Another kind of laplacian normalization, used in the matlab PVF code.
    Uses the formula: L = I - D^{-1} * W
348364
def ave_laplacian(self):
    ''''''
    W = self.matrix('dense')
    # calculate -inv(D)
    Dinv = W.sum(axis=0)
    mask = Dinv!=0
    Dinv[mask] = -1./Dinv[mask]
    # calculate -inv(D) * W
    lap = (Dinv * W.T).T
    # add I
    lap.flat[::W.shape[0]+1] += 1
    # symmetrize
    return (lap + lap.T) / 2.0
420853
def _laplace_fit(self,obj_type):
        

        # Get Mode and Inverse Hessian information
        y = self.fit(method='PML',printer=False)

        if y.ihessian is None:
            raise Exception("No Hessian information - Laplace approximation cannot be performed")
        else:

            self.latent_variables.estimation_method = 'Laplace'

            theta, Y, scores, states, states_var, X_names = self._categorize_model_output(self.latent_variables.get_z_values())

            # Change this in future
            try:
                latent_variables_store = self.latent_variables.copy()
            except:
                latent_variables_store = self.latent_variables

            return LaplaceResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
                model_type=self.model_type, latent_variables=latent_variables_store,data=Y,index=self.index,
                multivariate_model=self.multivariate_model,objective_object=obj_type, 
                method='Laplace',ihessian=y.ihessian,signal=theta,scores=scores,
                z_hide=self._z_hide,max_lag=self.max_lag,states=states,states_var=states_var)
753585
def f(self, v, t):
        
        ret = np.zeros(8)
        ret[0] = v[4]
        ret[4] = (self.He1*self.ke1*(self.g1*self.S(v[1]-v[2]) + self.u_mean) -
                  2*self.ke1*v[4] - self.ke1*self.ke1*v[0])

        ret[1] = v[5]
        ret[5] = (self.He2*self.ke2*(self.g2*self.S(v[0]) + self.p_mean) -
                  2*self.ke2*v[5] - self.ke2*self.ke2*v[1])

        ret[2] = v[6]
        ret[6] = (self.Hi*self.ki*self.g4*self.S(v[3]) - 2*self.ki*v[6] -
                  self.ki*self.ki*v[2])

        ret[3] = v[7]
        ret[7] = (self.He3*self.ke3*self.g3*self.S(v[1]-v[2]) -
                  2*self.ke3*v[7] - self.ke3*self.ke3*v[3])
        return ret
561948
function spfa(graph, s) {
  const distance = {};
  const previous = {};
  const queue = {};
  const isInQue = {};
  const cnt = {};
  let head = 0;
  let tail = 1;
  // initialize
  distance[s] = 0;
  queue[0] = s;
  isInQue[s] = true;
  cnt[s] = 1;
  graph.vertices.forEach(v => {
    if (v !== s) {
      distance[v] = Infinity;
      isInQue[v] = false;
      cnt[v] = 0;
    }
  });

  let currNode;
  while (head !== tail) {
    currNode = queue[head++];
    isInQue[currNode] = false;
    const neighbors = graph.neighbors(currNode);
    for (let i = 0; i < neighbors.length; i++) {
      const v = neighbors[i];
      // relaxation
      const newDistance = distance[currNode] + graph.edge(currNode, v);
      if (newDistance < distance[v]) {
        distance[v] = newDistance;
        previous[v] = currNode;
        if (!isInQue[v]) {
          queue[tail++] = v;
          isInQue[v] = true;
          cnt[v]++;
          if (cnt[v] > graph.vertices.size) {
            // indicates negative-weighted cycle
            return {
              distance: {}
            };
          }
        }
      }
    }
  }

  return {
    distance,
    previous
  };
}
818633
def _normalize(self):
        
        # Get length
        L = self.norm()
        if not L:
            raise ValueError('Quaternion cannot have 0-length.')
        # Correct
        self.w /= L
        self.x /= L
        self.y /= L
        self.z /= L
1049357
def compute_laplacian_matrix(affinity_matrix, method='auto', **kwargs):
    """"""
    if method == 'auto':
        method = 'geometric'
    return Laplacian.init(method, **kwargs).laplacian_matrix(affinity_matrix)
478444
def _lmder1_helical_valley():
    """"""
    tpi = 2 * np.pi

    def func(params, vec):
        if params[0] == 0:
            tmp1 = np.copysign(0.25, params[1])
        elif params[0] > 0:
            tmp1 = np.arctan(params[1] / params[0]) / tpi
        else:
            tmp1 = np.arctan(params[1] / params[0]) / tpi + 0.5

        tmp2 = np.sqrt(params[0]**2 + params[1]**2)

        vec[0] = 10 * (params[2] - 10 * tmp1)
        vec[1] = 10 * (tmp2 - 1)
        vec[2] = params[2]

    def jac(params, jac):
        temp = params[0]**2 + params[1]**2
        tmp1 = tpi * temp
        tmp2 = np.sqrt(temp)
        jac[0,0] = 100 * params[1] / tmp1
        jac[0,1] = 10 * params[0] / tmp2
        jac[0,2] = 0
        jac[1,0] = -100 * params[0] / tmp1
        jac[1,1] = 10 * params[1] / tmp2
        jac[2,0] = 10
        jac[2,1] = 0
        jac[1,2] = 0
        jac[2,2] = 1

    guess = np.asfarray([-1, 0, 0])

    _lmder1_driver(3, func, jac, guess,
                   50., 0.993652310343e-16,
                   [0.100000000000e+01, -0.624330159679e-17, 0.000000000000e+00])
    _lmder1_driver(3, func, jac, guess * 10,
                   0.102956301410e+03, 0.104467885065e-18,
                   [0.100000000000e+01, 0.656391080516e-20, 0.000000000000e+00])
    _lmder1_driver(3, func, jac, guess * 100,
                   0.991261822124e+03, 0.313877781195e-28,
                   [0.100000000000e+01, -0.197215226305e-29, 0.000000000000e+00])
84636
def vertex_normals(vertices, faces):
    """"""

    def normalize_v3(arr):
        ''' Normalize a numpy array of 3 component vectors shape=(n,3) '''
        lens = np.sqrt(arr[:, 0]**2 + arr[:, 1]**2 + arr[:, 2]**2)
        arr /= lens[:, np.newaxis]

    tris = vertices[faces]
    facenorms = np.cross(tris[::, 1] - tris[::, 0], tris[::, 2] - tris[::, 0])
    normalize_v3(facenorms)

    norm = np.zeros(vertices.shape, dtype=vertices.dtype)
    norm[faces[:, 0]] += facenorms
    norm[faces[:, 1]] += facenorms
    norm[faces[:, 2]] += facenorms
    normalize_v3(norm)
    return norm
179573
def ppf(q, df, loc=0.0, scale=1.0, gamma = 1.0):
        
        result = np.zeros(q.shape[0])
        probzero = Skewt.cdf(x=np.zeros(1),loc=np.zeros(1),df=df,gamma=gamma)
        result[q<probzero] = 1.0/gamma*ss.t.ppf(((np.power(gamma,2) + 1.0) * q[q<probzero])/2.0,df)
        result[q>=probzero] = gamma*ss.t.ppf((1.0 + 1.0/np.power(gamma,2))/2.0*(q[q >= probzero] - probzero) + 0.5, df)
        return result
1286802
def tf_smooth_eig_vec(self):
    """"""
    _, matrix_m = self.dual_object.get_full_psd_matrix()
    # Easier to think in terms of max so negating the matrix
    [eig_vals, eig_vectors] = tf.self_adjoint_eig(-matrix_m)
    exp_eig_vals = tf.exp(tf.divide(eig_vals, self.smooth_placeholder))
    scaling_factor = tf.reduce_sum(exp_eig_vals)
    # Multiplying each eig vector by exponential of corresponding eig value
    # Scaling factor normalizes the vector to be unit norm
    eig_vec_smooth = tf.divide(
        tf.matmul(eig_vectors, tf.diag(tf.sqrt(exp_eig_vals))),
        tf.sqrt(scaling_factor))
    return tf.reshape(
        tf.reduce_sum(eig_vec_smooth, axis=1),
        shape=[eig_vec_smooth.shape[0].value, 1])
804188
def normalize(self):
        
        cm = self.centerOfMass()
        coords = self.coordinates()
        if not len(coords):
            return
        pts = coords - cm
        xyz2 = np.sum(pts * pts, axis=0)
        scale = 1 / np.sqrt(np.sum(xyz2) / len(pts))
        t = vtk.vtkTransform()
        t.Scale(scale, scale, scale)
        t.Translate(-cm)
        tf = vtk.vtkTransformPolyDataFilter()
        tf.SetInputData(self.poly)
        tf.SetTransform(t)
        tf.Update()
        return self.updateMesh(tf.GetOutput())
543285
def Svd(a, uv, full):
    
    u, s, v = np.linalg.svd(a, full_matrices=full, compute_uv=uv)
    return s, u, v
1363622
def normalized(self):
    

    norm = self.magnitude()
    return Vector(self.x / norm, self.y / norm, self.z / norm)
1067945
def cublasDsymv(handle, uplo, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasDsymv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], n, 
                                       ctypes.byref(ctypes.c_double(alpha)), 
                                       int(A), lda, int(x), incx, 
                                       ctypes.byref(ctypes.c_double(beta)), 
                                       int(y), incy)
    cublasCheckStatus(status)
1204488
function normalize(vec3) {
    if (vec3.x === 0 && vec3.y === 0 && vec3.z === 0) {
        return vec3;
    }
    const length = Math.sqrt(vec3.x * vec3.x + vec3.y * vec3.y + vec3.z * vec3.z);
    const invScalar = 1 / length;
    vec3.x *= invScalar;
    vec3.y *= invScalar;
    vec3.z *= invScalar;
    return vec3;
}
695968
def PCA(Y, components):
	

	N,D = Y.shape
	sv = linalg.svd(Y, full_matrices=0);
	[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
	v = s0.std(axis=0)
	s0 /= v;
	w0 *= v;
	return [s0, w0]

	if N>D:
		sv = linalg.svd(Y, full_matrices=0);
		[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
		v = s0.std(axis=0)
		s0 /= v;
		w0 *= v;
		return [s0, w0]
	else:
		K=np.cov(Y)
		sv = linalg.eigh(K)
		std_var = np.sqrt(sv[0])
		pc = sv[1]*std_var[np.newaxis(),0]
		#import ipdb
		#ipdb.set_trace()
		return [pc,std_var]
1288174
def interpolate2dUnstructuredIDW(x, y, v, grid, power=2):
    '''
    
    '''
    n = len(v)
    gx = grid.shape[0]
    gy = grid.shape[1]
    for i in range(gx):
        for j in range(gy):
            overPx = False  # if pixel position == point position
            sumWi = 0.0
            value = 0.0

            for k in range(n):
                xx = x[k]
                yy = y[k]
                vv = v[k]
                if xx == i and yy == j:
                    grid[i, j] = vv
                    overPx = True
                    break
                # weight from inverse distance:
                wi = 1 / ((xx - i)**2 + (yy - j)**2)**(0.5 * power)
                sumWi += wi
                value += wi * vv
            if not overPx:
                grid[i, j] = value / sumWi
    return grid