Calculate standard vertical modes. Simply
            the eigenvectors of the stretching matrix S
961416
def vertical_modes(self):
        

        evals,evecs = np.linalg.eig(-self.S)

        asort = evals.argsort()

        # deformation wavenumbers and radii
        self.kdi2 = evals[asort]
        self.radii = np.zeros_like(self.kdi2)
        self.radii[0] = self.g*self.H/np.abs(self.f) # barotropic def. radius
        self.radii[1:] = 1./np.sqrt(self.kdi2[1:])

        # eigenstructure
        self.pmodes = evecs[:,asort]

        # normalize to have unit L2-norm
        Ai = (self.H / (self.Hi[:,np.newaxis]*(self.pmodes**2)).sum(axis=0))**0.5
        self.pmodes = Ai[np.newaxis,:]*self.pmodes
495244
def index_search(right_eigenvectors):
    

    num_micro, num_eigen = right_eigenvectors.shape

    index = np.zeros(num_eigen, 'int')

    # first vertex: row with largest norm
    index[0] = np.argmax(
            [norm(right_eigenvectors[i]) for i in range(num_micro)])

    ortho_sys = right_eigenvectors - np.outer(np.ones(num_micro),
                                              right_eigenvectors[index[0]])

    for j in range(1, num_eigen):
        temp = ortho_sys[index[j - 1]].copy()
        for l in range(num_micro):
            ortho_sys[l] -= temp * dot(ortho_sys[l], temp)

        dist_list = np.array([norm(ortho_sys[l]) for l in range(num_micro)])

        index[j] = np.argmax(dist_list)

        ortho_sys /= dist_list.max()

    return index
1300178
def diagonalize(self):
        ''''''
        self.eigvals, self.eigvecs = np.linalg.eig(
            (self.tensor.transpose() + self.tensor) / 2.0)
        self.eigvals = np.diag(np.dot(
            np.dot(self.eigvecs.transpose(), self.tensor), self.eigvecs))
623514
def row_standard_coordinates(self, X):
        
        utils.validation.check_is_fitted(self, 's_')
        return self.row_coordinates(X).div(self.eigenvalues_, axis='columns')
417135
@Override
	protected void vertical( GrayF32 intensity ) {
		float[] hXX = horizXX.data;
		float[] hXY = horizXY.data;
		float[] hYY = horizYY.data;
		final float[] inten = intensity.data;

		final int imgHeight = horizXX.getHeight();
		final int imgWidth = horizXX.getWidth();

		final int kernelWidth = radius * 2 + 1;

		final int startX = radius;
		final int endX = imgWidth - radius;

		final int backStep = kernelWidth * imgWidth;

		BoofConcurrency.loopBlocks(radius,imgHeight-radius,(y0,y1)->{
		final float[] tempXX = work.pop();
		final float[] tempXY = work.pop();
		final float[] tempYY = work.pop();
		for (int x = startX; x < endX; x++) {
			// defines the A matrix, from which the eigenvalues are computed
			int srcIndex = x + (y0-radius)*imgWidth;
			int destIndex = imgWidth * y0 + x;
			float totalXX = 0, totalXY = 0, totalYY = 0;

			int indexEnd = srcIndex + imgWidth * kernelWidth;
			for (; srcIndex < indexEnd; srcIndex += imgWidth) {
				totalXX += hXX[srcIndex];
				totalXY += hXY[srcIndex];
				totalYY += hYY[srcIndex];
			}

			tempXX[x] = totalXX;
			tempXY[x] = totalXY;
			tempYY[x] = totalYY;

			// compute the eigen values
			inten[destIndex] = this.intensity.compute(totalXX,totalXY,totalYY);
			destIndex += imgWidth;
		}

		// change the order it is processed in to reduce cache misses
		for (int y = y0 + 1; y < y1; y++) {
			int srcIndex = (y + radius) * imgWidth + startX;
			int destIndex = y * imgWidth + startX;

			for (int x = startX; x < endX; x++, srcIndex++, destIndex++) {
				float totalXX = tempXX[x] - hXX[srcIndex - backStep];
				tempXX[x] = totalXX += hXX[srcIndex];
				float totalXY = tempXY[x] - hXY[srcIndex - backStep];
				tempXY[x] = totalXY += hXY[srcIndex];
				float totalYY = tempYY[x] - hYY[srcIndex - backStep];
				tempYY[x] = totalYY += hYY[srcIndex];

				inten[destIndex] = this.intensity.compute(totalXX,totalXY,totalYY);
			}
		}
		work.recycle(tempXX);
		work.recycle(tempXY);
		work.recycle(tempYY);
		});
	}
1130572
def eigenvalue(nodes, node=None):
    
    if node is None:
        return sorted([eigenvalue(nodes, n) for n in nodes])[0]  # return lowest eigenvalue
    else:
        return len([1 for n in nodes if hops(node, n)])
1345409
function matrixToGravityArray(mat) {
  const rows = mat.length;
  const cols = mat[0].length;
  const rowTip = rows - 1;
  const colTip = cols - 1;
  const gravityArr = [];

  let row, col, sum, topSum = 0;
  for (row = 0; row < rows; row++) {
    for (col = 0; col < cols; col++) {
      sum = mat[row][col];
      if (!sum) continue;

      topSum = Math.max(topSum, sum);

      gravityArr.push({ row, col, sum });
    }
  }

  return {
    gravityArr,
    topSum
  };
}
497121
def _spectral_norm(self):
        
        w = self.params.get('weight').data(self.ctx)
        w_mat = nd.reshape(w, [w.shape[0], -1])

        _u = self.u.data(self.ctx)
        _v = None

        for _ in range(POWER_ITERATION):
            _v = nd.L2Normalization(nd.dot(_u, w_mat))
            _u = nd.L2Normalization(nd.dot(_v, w_mat.T))

        sigma = nd.sum(nd.dot(_u, w_mat) * _v)
        if sigma == 0.:
            sigma = EPSILON

        with autograd.pause():
            self.u.set_data(_u)

        return w / sigma
928914
def _svd(cls, matrix, num_concepts=5):
        
        u, s, v = svds(matrix, k=num_concepts)
        return u, s, v
848895
def check_psd(matrix, tolerance=1e-6):
    

    hermitian = (matrix + matrix.T.conjugate()) / 2
    eigenvalues = np.linalg.eigh(hermitian)[0]
    return (eigenvalues > -tolerance).all()
1181849
def compute_eig(self):
        
        self.eig1 = _SHGrid.from_array(_np.zeros_like(self.vxx.data),
                                       grid='DH')
        self.eig2 = _SHGrid.from_array(_np.zeros_like(self.vxx.data),
                                       grid='DH')
        self.eig3 = _SHGrid.from_array(_np.zeros_like(self.vxx.data),
                                       grid='DH')

        for i in range(self.nlat):
            for j in range(self.nlon):
                a = _np.array([[self.vxx.data[i, j],
                                self.vxy.data[i, j],
                                self.vxz.data[i, j]],
                               [self.vyx.data[i, j],
                                self.vyy.data[i, j],
                                self.vyz.data[i, j]],
                               [self.vzx.data[i, j],
                                self.vzy.data[i, j],
                                self.vzz.data[i, j]]])

                eigs = _eigvalsh(a)

                self.eig1.data[i, j] = eigs[2]
                self.eig2.data[i, j] = eigs[1]
                self.eig3.data[i, j] = eigs[0]
1219686
def _check_scalar_vertical_extents(self, ds, z_variable):
        '''
        
        '''
        vert_min = ds.geospatial_vertical_min
        vert_max = ds.geospatial_vertical_max
        msgs = []
        total = 2

        zvalue = ds.variables[z_variable][:].item()
        if not np.isclose(vert_min, vert_max):
            msgs.append("geospatial_vertical_min != geospatial_vertical_max for scalar depth values, %s != %s" % (
                vert_min,
                vert_max
            ))

        if not np.isclose(vert_max, zvalue):
            msgs.append("geospatial_vertical_max != %s values, %s != %s" % (
                z_variable,
                vert_max,
                zvalue
            ))

        return Result(BaseCheck.MEDIUM,
                      (total - len(msgs), total),
                      'geospatial_vertical_extents_match',
                      msgs)
1028725
private static EigenPair[] processDecomposition(EigenvalueDecomposition evd) {
    double[] eigenvalues = evd.getRealEigenvalues();
    double[][] eigenvectors = evd.getV();

    EigenPair[] eigenPairs = new EigenPair[eigenvalues.length];
    for(int i = 0; i < eigenvalues.length; i++) {
      double e = Math.abs(eigenvalues[i]);
      double[] v = VMath.getCol(eigenvectors, i);
      eigenPairs[i] = new EigenPair(v, e);
    }
    Arrays.sort(eigenPairs, Comparator.reverseOrder());
    return eigenPairs;
  }
324750
def eigendecompose(tensor, normalise=False):
    '''
    
    '''
    if normalise:
        tensor, tensor_norm = normalise_tensor(tensor)
    else:
        tensor_norm = 1.

    eigvals, eigvects = np.linalg.eigh(tensor, UPLO='U')

    isrt = np.argsort(eigvals)
    eigenvalues = eigvals[isrt] * tensor_norm
    eigenvectors = eigvects[:, isrt]
    return eigenvalues, eigenvectors
586688
def LaplaceCentreWeight(self):
        """"""

        sz = [1,] * self.S.ndim
        for ax in self.axes:
            sz[ax] = self.S.shape[ax]
        lcw = 2*len(self.axes)*np.ones(sz, dtype=self.dtype)
        for ax in self.axes:
            lcw[(slice(None),)*ax + ([0, -1],)] -= 1.0
        return lcw
140830
def _project_eigenvectors(self):
        """"""

        self._p_eigenvectors = []
        for vecs_q in self._eigenvectors:
            p_vecs_q = []
            for vecs in vecs_q.T:
                p_vecs_q.append(np.dot(vecs.reshape(-1, 3),
                                       self._projection_direction))
            self._p_eigenvectors.append(np.transpose(p_vecs_q))
        self._p_eigenvectors = np.array(self._p_eigenvectors)
1046616
private void hqr2(Matrix H, Matrix V, Vector d, Vector e) {

        // This is derived from the Algol procedure hqr2,
        // by Martin and Wilkinson, Handbook for Auto. Comp.,
        // Vol.ii-Linear Algebra, and the corresponding
        // Fortran subroutine in EISPACK.

        // Initialize

        int nn = e.length();
        int n = nn - 1;
        int low = 0;
        int high = nn - 1;
        double eps = Math.pow(2.0, -52.0);
        double exshift = 0.0;
        double p = 0;
        double q = 0;
        double r = 0;
        double s = 0;
        double z = 0;
        double t;
        double w;
        double x;
        double y;

        // Store roots isolated by balanc and compute matrix norm

        double norm = 0.0;
        for (int i = 0; i < nn; i++) {
            if (i < low | i > high) {
                d.set(i, H.get(i, i));
                e.set(i, 0.0);
            }
            for (int j = Math.max(i - 1, 0); j < nn; j++) {
                norm = norm + Math.abs(H.get(i, j));
            }
        }

        // Outer loop over eigenvalue index

        int iter = 0;
        while (n >= low) {

            // Look for single small sub-diagonal element

            int l = n;
            while (l > low) {
                s = Math.abs(H.get(l - 1, l - 1)) 
                    + Math.abs(H.get(l, l));

                if (s == 0.0) {
                    s = norm;
                }
                if (Math.abs(H.get(l, l - 1)) < eps * s) {
                    break;
                }
                l--;
            }

            // Check for convergence
            // One root found

            if (l == n) {
                H.updateAt(n, n, Matrices.asPlusFunction(exshift));
                d.set(n, H.get(n, n));
                e.set(n, 0.0);
                n--;
                iter = 0;

                // Two roots found

            } else if (l == n - 1) {
                w = H.get(n, n - 1) * H.get(n - 1, n);
                p = (H.get(n - 1, n - 1) - H.get(n, n)) / 2.0;
                q = p * p + w;
                z = Math.sqrt(Math.abs(q));
                H.updateAt(n, n, Matrices.asPlusFunction(exshift));
                H.updateAt(n - 1, n - 1, Matrices.asPlusFunction(exshift));
                x = H.get(n, n);

                // Real pair

                if (q >= 0) {
                    if (p >= 0) {
                        z = p + z;
                    } else {
                        z = p - z;
                    }
                    d.set(n - 1, x + z);
                    d.set(n, d.get(n - 1));
                    if (z != 0.0) {
                        d.set(n, x - w / z);
                    }
                    e.set(n - 1, 0.0);
                    e.set(n, 0.0);
                    x = H.get(n, n - 1);
                    s = Math.abs(x) + Math.abs(z);
                    p = x / s;
                    q = z / s;
                    r = Math.sqrt(p * p + q * q);
                    p = p / r;
                    q = q / r;

                    // Row modification

                    for (int j = n - 1; j < nn; j++) {
                        z = H.get(n - 1, j);
                        H.set(n - 1, j, q * z + p * H.get(n, j));
                        H.set(n, j, q * H.get(n, j) - p * z);
                    }

                    // Column modification

                    for (int i = 0; i <= n; i++) {
                        z = H.get(i, n - 1);
                        H.set(i, n - 1, q * z + p * H.get(i, n));
                        H.set(i, n, q * H.get(i, n) - p * z);
                    }

                    // Accumulate transformations

                    for (int i = low; i <= high; i++) {
                        z = V.get(i, n - 1);
                        V.set(i, n - 1, q * z + p * V.get(i, n));
                        V.set(i, n, q * V.get(i, n) - p * z);
                    }

                    // Complex pair

                } else {
                    d.set(n - 1, x + p);
                    d.set(n, x + p);
                    e.set(n - 1, z);
                    e.set(n, -z);
                }
                n = n - 2;
                iter = 0;

                // No convergence yet

            } else {

                // Form shift

                x = H.get(n, n);
                y = 0.0;
                w = 0.0;
                if (l < n) {
                    y = H.get(n - 1, n - 1);
                    w = H.get(n, n - 1) * H.get(n - 1, n);
                }

                // Wilkinson's original ad hoc shift

                if (iter == 10) {
                    exshift += x;
                    for (int i = low; i <= n; i++) {
                        H.updateAt(i, i, Matrices.asMinusFunction(x));
                    }
                    s = Math.abs(H.get(n, n - 1)) 
                        + Math.abs(H.get(n - 1, n - 2));
                    x = y = 0.75 * s; // WTF ???
                    w = -0.4375 * s * s; // Are you kidding me???
                }

                // MATLAB's new ad hoc shift

                if (iter == 30) {
                    s = (y - x) / 2.0;
                    s = s * s + w;
                    if (s > 0) {
                        s = Math.sqrt(s);
                        if (y < x) {
                            s = -s;
                        }
                        s = x - w / ((y - x) / 2.0 + s);
                        for (int i = low; i <= n; i++) {
                            H.updateAt(i, i, Matrices.asMinusFunction(s));
                        }
                        exshift += s;
                        x = y = w = 0.964;
                    }
                }

                iter = iter + 1; // (Could check iteration count here.)

                // Look for two consecutive small sub-diagonal elements

                int m = n - 2;
                while (m >= l) {
                    z = H.get(m, m);
                    r = x - z;
                    s = y - z;
                    p = (r * s - w) / H.get(m + 1, m) 
                                      + H.get(m, m + 1);
                    q = H.get(m + 1, m + 1) - z - r - s;
                    r = H.get(m + 2, m + 1);
                    s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                    p = p / s;
                    q = q / s;
                    r = r / s;
                    if (m == l) {
                        break;
                    }
                    if (Math.abs(H.get(m, m - 1)) * (Math.abs(q) + Math.abs(r)) < eps
                            * (Math.abs(p) * (Math.abs(H.get(m - 1, m - 1))
                                    + Math.abs(z) + Math.abs(H.get(m + 1, m + 1))))) {
                        break;
                    }
                    m--;
                }

                for (int i = m + 2; i <= n; i++) {
                    H.set(i, i - 2, 0.0);
                    if (i > m + 2) {
                        H.set(i, i - 3, 0.0);
                    }
                }

                // Double QR step involving rows l:n and columns m:n

                for (int k = m; k <= n - 1; k++) {
                    boolean notlast = (k != n - 1);
                    if (k != m) {
                        p = H.get(k, k - 1);
                        q = H.get(k + 1, k - 1);
                        r = (notlast ? H.get(k + 2, k - 1) : 0.0);
                        x = Math.abs(p) + Math.abs(q) + Math.abs(r);
                        if (x == 0.0) {
                            continue;
                        }
                        p = p / x;
                        q = q / x;
                        r = r / x;
                    }

                    s = Math.sqrt(p * p + q * q + r * r);
                    if (p < 0) {
                        s = -s;
                    }
                    if (s != 0) {
                        if (k != m) {
                            H.set(k, k - 1, -s * x);
                        } else if (l != m) {
                            H.updateAt(k, k - 1, Matrices.INV_FUNCTION);
                        }
                        p = p + s;
                        x = p / s;
                        y = q / s;
                        z = r / s;
                        q = q / p;
                        r = r / p;

                        // Row modification

                        for (int j = k; j < nn; j++) {
                            p = H.get(k, j) + q * H.get(k + 1, j);
                            if (notlast) {
                                p = p + r * H.get(k + 2, j);
                                H.updateAt(k + 2, j,
                                    Matrices.asMinusFunction(p * z));
                            }
                            H.updateAt(k, j, Matrices.asMinusFunction(p * x));
                            H.updateAt(k + 1, j, Matrices.asMinusFunction(p * y));
                        }

                        // Column modification

                        for (int i = 0; i <= Math.min(n, k + 3); i++) {
                            p = x * H.get(i, k) + y 
                                * H.get(i, k + 1);
                            if (notlast) {
                                p = p + z * H.get(i, k + 2);
                                H.updateAt(i, k + 2,
                                    Matrices.asMinusFunction(p * r));
                            }
                            H.updateAt(i, k, Matrices.asMinusFunction(p));
                            H.updateAt(i, k + 1, Matrices.asMinusFunction(p * q));
                        }

                        // Accumulate transformations

                        for (int i = low; i <= high; i++) {
                            p = x * V.get(i, k) + y 
                                * V.get(i, k + 1);
                            if (notlast) {
                                p = p + z * V.get(i, k + 2);
                                V.updateAt(i, k + 2,
                                    Matrices.asMinusFunction(p * r));
                            }
                            V.updateAt(i, k, Matrices.asMinusFunction(p));
                            V.updateAt(i, k + 1, Matrices.asMinusFunction(p * q));
                        }
                    } // (s != 0)
                } // k loop
            } // check convergence
        } // while (n >= low)

        // Backsubstitute to find vectors of upper triangular form

        if (norm == 0.0) {
            return;
        }

        for (n = nn - 1; n >= 0; n--) {
            p = d.get(n);
            q = e.get(n);

            // Real vector

            if (q == 0) {
                int l = n;
                H.set(n, n, 1.0);
                for (int i = n - 1; i >= 0; i--) {
                    w = H.get(i, i) - p;
                    r = 0.0;
                    for (int j = l; j <= n; j++) {
                        r = r + H.get(i, j) * H.get(j, n);
                    }
                    if (e.get(i) < 0.0) {
                        z = w;
                        s = r;
                    } else {
                        l = i;
                        if (e.get(i) == 0.0) {
                            if (w != 0.0) {
                                H.set(i, n, -r / w);
                            } else {
                                H.set(i, n, -r / (eps * norm));
                            }

                            // Solve real equations

                        } else {
                            x = H.get(i, i + 1);
                            y = H.get(i + 1, i);
                            q = (d.get(i) - p) * (d.get(i) - p) 
                                + e.get(i) * e.get(i);
                            t = (x * s - z * r) / q;
                            H.set(i, n, t);
                            if (Math.abs(x) > Math.abs(z)) {
                                H.set(i + 1, n, (-r - w * t) / x);
                            } else {
                                H.set(i + 1, n, (-s - y * t) / z);
                            }
                        }

                        // Overflow control

                        t = Math.abs(H.get(i, n));
                        if ((eps * t) * t > 1) {
                            for (int j = i; j <= n; j++) {
                                H.updateAt(j, n, Matrices.asDivFunction(t));
                            }
                        }
                    }
                }

                // Complex vector

            } else if (q < 0) {
                int l = n - 1;

                // Last vector component imaginary so matrix is triangular

                if (Math.abs(H.get(n, n - 1)) 
                        > Math.abs(H.get(n - 1, n))) {
                    H.set(n - 1, n - 1, q / H.get(n, n - 1));
                    H.set(n - 1, n, -(H.get(n, n) - p) 
                                            / H.get(n, n - 1));
                } else {
                    double[] cdiv = cdiv(0.0, -H.get(n - 1, n),
                         H.get(n - 1, n - 1) - p, q);
                    H.set(n - 1, n - 1, cdiv[0]);
                    H.set(n - 1, n, cdiv[1]);
                }
                H.set(n, n - 1, 0.0);
                H.set(n, n, 1.0);
                for (int i = n - 2; i >= 0; i--) {
                    double ra;
                    double sa;
                    double vr;
                    double vi;
                    ra = 0.0;
                    sa = 0.0;
                    for (int j = l; j <= n; j++) {
                        ra = ra + H.get(i, j) * H.get(j, n - 1);
                        sa = sa + H.get(i, j) * H.get(j, n);
                    }
                    w = H.get(i, i) - p;

                    if (e.get(i) < 0.0) {
                        z = w;
                        r = ra;
                        s = sa;
                    } else {
                        l = i;
                        if (e.get(i) == 0) {
                            double[] cdiv = cdiv(-ra, -sa, w, q);
                            H.set(i, n - 1, cdiv[0]);
                            H.set(i, n, cdiv[1]);
                        } else {

                            // Solve complex equations

                            x = H.get(i, i + 1);
                            y = H.get(i + 1, i);
                            vr = (d.get(i) - p) * (d.get(i) - p) 
                                 + e.get(i) * e.get(i) - q * q;
                            vi = (d.get(i) - p) * 2.0 * q;
                            if (vr == 0.0 & vi == 0.0) {
                                vr = eps
                                        * norm
                                        * (Math.abs(w) + Math.abs(q)
                                                + Math.abs(x) + Math.abs(y) + Math
                                                    .abs(z));
                            }
                            double[] cdiv = cdiv(x * r - z * ra + q * sa,
                                                 x * s - z * sa - q * ra, vr, vi);
                            H.set(i, n - 1, cdiv[0]);
                            H.set(i, n, cdiv[1]);
                            if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
                                H.set(i + 1, n - 1, (-ra - w 
                                             * H.get(i, n - 1) + q
                                             * H.get(i, n)) / x);
                                H.set(i + 1, n, (-sa - w 
                                             * H.get(i, n) - q 
                                             * H.get(i, n - 1)) / x);
                            } else {
                                cdiv = cdiv(-r - y 
                                     * H.get(i, n - 1), -s - y 
                                     * H.get(i, n), z, q);
                                H.set(i + 1, n - 1, cdiv[0]);
                                H.set(i + 1, n, cdiv[1]);
                            }
                        }

                        // Overflow control

                        t = Math.max(Math.abs(H.get(i, n - 1)), 
                                     Math.abs(H.get(i, n)));
                        if ((eps * t) * t > 1) {
                            for (int j = i; j <= n; j++) {
                                H.updateAt(j, n - 1, Matrices.asDivFunction(t));
                                H.updateAt(j, n, Matrices.asDivFunction(t));
                            }
                        }
                    }
                }
            }
        }

        // Vectors of isolated roots

        for (int i = 0; i < nn; i++) {
            if (i < low | i > high) {
                for (int j = i; j < nn; j++) {
                    V.set(i, j, H.get(i, j));
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix

        for (int j = nn - 1; j >= low; j--) {
            for (int i = low; i <= high; i++) {
                z = 0.0;
                for (int k = low; k <= Math.min(j, high); k++) {
                    z = z + V.get(i, k) * H.get(k, j);
                }
                V.set(i, j, z);
            }
        }
    }