Computes a basis (the principle components) from the most dominant eigenvectors.
3698
public void computeBasis() {
		if (sampleIndex != numSamples)
			throw new IllegalArgumentException("Not all the data has been added");
		if (numComponents > numSamples)
			throw new IllegalArgumentException(
					"More data needed to compute the desired number of components");

		means = new DenseMatrix64F(sampleSize, 1);
		// compute the mean of all the samples
		for (int i = 0; i < numSamples; i++) {
			for (int j = 0; j < sampleSize; j++) {
				double val = means.get(j);
				means.set(j, val + A.get(i, j));
			}
		}
		for (int j = 0; j < sampleSize; j++) {
			double avg = means.get(j) / numSamples;
			means.set(j, avg);
		}

		// subtract the mean from the original data
		for (int i = 0; i < numSamples; i++) {
			for (int j = 0; j < sampleSize; j++) {
				A.set(i, j, A.get(i, j) - means.get(j));
			}
		}

		// compute SVD and save time by not computing U
		SingularValueDecomposition<DenseMatrix64F> svd = DecompositionFactory.svd(numSamples, sampleSize,
				false, true, compact);
		if (!svd.decompose(A))
			throw new RuntimeException("SVD failed");

		V_t = svd.getV(null, true);
		W = svd.getW(null);

		// singular values are in an arbitrary order initially and need to be sorted in descending order
		SingularOps.descendingOrder(null, false, W, V_t, true);

		// strip off unneeded components and find the basis
		V_t.reshape(numComponents, sampleSize, true);

	}
209907
def A(self):
        """"""
        wishart = self._wishart_cdf

        base = np.zeros([wishart.size, wishart.size])
        if wishart.n_min % 2:
            # If matrix has odd number of elements, we need to append a
            # row and a col, in order for the pfaffian algorithm to work
            base = self._make_size_even(base)
        return base
1158312
private static void compute(Kernel kernel, int[] in, int[] out, int width, int height, boolean alpha, int edge)
    {
        final float[] matrix = kernel.getMatrix();
        final int cols = kernel.getWidth();
        final int cols2 = cols / 2;

        for (int y = 0; y < height; y++)
        {
            int index = y;
            final int ioffset = y * width;
            for (int x = 0; x < width; x++)
            {
                compute(matrix, in, out, x, index, ioffset, cols2, width, alpha, edge);
                index += height;
            }
        }
    }
1267290
def _to_diagonally_dominant(mat):
    """"""
    mat += np.diag(np.sum(mat != 0, axis=1) + 0.01)
    return mat
776057
private void buildHierarchy(Relation<V> database, Clustering<SubspaceModel> clustering, List<Cluster<SubspaceModel>> clusters, int dimensionality) {
    StringBuilder msg = LOG.isDebugging() ? new StringBuilder() : null;
    final int db_dim = RelationUtil.dimensionality(database);
    Hierarchy<Cluster<SubspaceModel>> hier = clustering.getClusterHierarchy();

    for(int i = 0; i < clusters.size() - 1; i++) {
      Cluster<SubspaceModel> c_i = clusters.get(i);
      final Subspace s_i = c_i.getModel().getSubspace();
      int subspaceDim_i = dimensionality - s_i.dimensionality();
      NumberVector ci_centroid = ProjectedCentroid.make(s_i.getDimensions(), database, c_i.getIDs());
      long[] pv1 = s_i.getDimensions();

      for(int j = i + 1; j < clusters.size(); j++) {
        Cluster<SubspaceModel> c_j = clusters.get(j);
        final Subspace s_j = c_j.getModel().getSubspace();
        int subspaceDim_j = dimensionality - s_j.dimensionality();

        if(subspaceDim_i < subspaceDim_j) {
          if(msg != null) {
            msg.append("\n l_i=").append(subspaceDim_i).append(" pv_i=[").append(BitsUtil.toStringLow(s_i.getDimensions(), db_dim)).append(']') //
                .append("\n l_j=").append(subspaceDim_j).append(" pv_j=[").append(BitsUtil.toStringLow(s_j.getDimensions(), db_dim)).append(']');
          }

          // noise level reached
          if(s_j.dimensionality() == 0) {
            // no parents exists -> parent is noise
            if(hier.numParents(c_i) == 0) {
              clustering.addChildCluster(c_j, c_i);
              if(msg != null) {
                msg.append("\n [").append(BitsUtil.toStringLow(s_j.getDimensions(), db_dim)) //
                    .append("] is parent of [").append(BitsUtil.toStringLow(s_i.getDimensions(), db_dim)).append(']');
              }
            }
          }
          else {
            NumberVector cj_centroid = ProjectedCentroid.make(c_j.getModel().getDimensions(), database, c_j.getIDs());
            long[] pv2 = s_j.getDimensions();
            long[] commonPreferenceVector = BitsUtil.andCMin(pv1, pv2);
            int subspaceDim = subspaceDimensionality(ci_centroid, cj_centroid, pv1, pv2, commonPreferenceVector);

            double d = weightedDistance(ci_centroid, cj_centroid, commonPreferenceVector);
            if(msg != null) {
              msg.append("\n dist = ").append(subspaceDim);
            }

            if(subspaceDim == subspaceDim_j) {
              if(msg != null) {
                msg.append("\n d = ").append(d);
              }
              if(d <= 2 * epsilon) {
                // no parent exists or c_j is not a parent of the already
                // existing parents
                if(hier.numParents(c_i) == 0 || !isParent(database, c_j, hier.iterParents(c_i), db_dim)) {
                  clustering.addChildCluster(c_j, c_i);
                  if(msg != null) {
                    msg.append("\n [").append(BitsUtil.toStringLow(s_j.getDimensions(), db_dim)) //
                        .append("] is parent of [") //
                        .append(BitsUtil.toStringLow(s_i.getDimensions(), db_dim)).append(']');
                  }
                }
              }
              else {
                throw new RuntimeException("Should never happen: d = " + d);
              }
            }
          }
        }
      }
    }
    if(msg != null) {
      LOG.debug(msg.toString());
    }
  }
1056857
def compute_svd(X, n_components, n_iter, random_state, engine):
    """"""

    # Determine what SVD engine to use
    if engine == 'auto':
        engine = 'sklearn'

    # Compute the SVD
    if engine == 'fbpca':
        if FBPCA_INSTALLED:
            U, s, V = fbpca.pca(X, k=n_components, n_iter=n_iter)
        else:
            raise ValueError('fbpca is not installed; please install it if you want to use it')
    elif engine == 'sklearn':
        U, s, V = extmath.randomized_svd(
            X,
            n_components=n_components,
            n_iter=n_iter,
            random_state=random_state
        )
    else:
        raise ValueError("engine has to be one of ('auto', 'fbpca', 'sklearn')")

    U, V = extmath.svd_flip(U, V)

    return U, s, V
1327131
private void _symmetricTred2 ()
  {
    // This is derived from the Algol procedures tred2 by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
    }

    // Householder reduction to tridiagonal form.

    for (int i = m_nDim - 1; i > 0; i--)
    {
      // Scale to avoid under/overflow.
      double scale = 0.0;
      double h = 0.0;
      for (int k = 0; k < i; k++)
        scale += MathHelper.abs (m_aEVd[k]);
      if (scale == 0.0)
      {
        m_aEVe[i] = m_aEVd[i - 1];
        for (int j = 0; j < i; j++)
        {
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
          m_aEigenVector[j][i] = 0.0;
        }
      }
      else
      {
        // Generate Householder vector.

        for (int k = 0; k < i; k++)
        {
          m_aEVd[k] /= scale;
          h += m_aEVd[k] * m_aEVd[k];
        }
        double f = m_aEVd[i - 1];
        double g = Math.sqrt (h);
        if (f > 0)
        {
          g = -g;
        }
        m_aEVe[i] = scale * g;
        h -= f * g;
        m_aEVd[i - 1] = f - g;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] = 0.0;
        }

        // Apply similarity transformation to remaining columns.
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          m_aEigenVector[j][i] = f;
          g = m_aEVe[j] + m_aEigenVector[j][j] * f;
          for (int k = j + 1; k <= i - 1; k++)
          {
            g += m_aEigenVector[k][j] * m_aEVd[k];
            m_aEVe[k] += m_aEigenVector[k][j] * f;
          }
          m_aEVe[j] = g;
        }
        f = 0.0;
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] /= h;
          f += m_aEVe[j] * m_aEVd[j];
        }
        final double hh = f / (h + h);
        for (int j = 0; j < i; j++)
        {
          m_aEVe[j] -= hh * m_aEVd[j];
        }
        for (int j = 0; j < i; j++)
        {
          f = m_aEVd[j];
          g = m_aEVe[j];
          for (int k = j; k <= i - 1; k++)
          {
            m_aEigenVector[k][j] -= (f * m_aEVe[k] + g * m_aEVd[k]);
          }
          m_aEVd[j] = m_aEigenVector[i - 1][j];
          m_aEigenVector[i][j] = 0.0;
        }
      }
      m_aEVd[i] = h;
    }

    // Accumulate transformations.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      m_aEigenVector[m_nDim - 1][i] = m_aEigenVector[i][i];
      m_aEigenVector[i][i] = 1.0;
      final double h = m_aEVd[i + 1];
      if (h != 0.0)
      {
        for (int k = 0; k <= i; k++)
        {
          m_aEVd[k] = m_aEigenVector[k][i + 1] / h;
        }
        for (int j = 0; j <= i; j++)
        {
          double g = 0.0;
          for (int k = 0; k <= i; k++)
          {
            g += m_aEigenVector[k][i + 1] * m_aEigenVector[k][j];
          }
          for (int k = 0; k <= i; k++)
          {
            m_aEigenVector[k][j] -= g * m_aEVd[k];
          }
        }
      }
      for (int k = 0; k <= i; k++)
      {
        m_aEigenVector[k][i + 1] = 0.0;
      }
    }
    for (int j = 0; j < m_nDim; j++)
    {
      m_aEVd[j] = m_aEigenVector[m_nDim - 1][j];
      m_aEigenVector[m_nDim - 1][j] = 0.0;
    }
    m_aEigenVector[m_nDim - 1][m_nDim - 1] = 1.0;
    m_aEVe[0] = 0.0;
  }
86038
def eigenvalues_hadamard(matrix1, matrix2):
    

    matrix1 = array(matrix1)  # as arrays
    matrix2 = array(matrix2)

    eig_a = eigvalsh(matrix1)  # eigenvalues (optimized for Hermitian matrices)
    eig_b = eigvalsh(matrix2)

    return eig_a, eig_b
735329
private void hqr2() {
    // FIXME: does this fail on NaN/inf values?

    // 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
    final int nn = this.n, low = 0, high = nn - 1;

    // Store roots isolated by balanc and compute matrix norm
    double norm = 0.;
    for(int i = 0; i < nn; i++) {
      // If eventually allowing low,high to be set, use this:
      // if(i < low || i > high) { d[i] = H[i][i]; e[i] = 0.; }
      for(int j = (i > 0 ? i - 1 : 0); j < nn; j++) {
        norm += abs(H[i][j]);
      }
    }

    // Outer loop over eigenvalue index
    {
      double exshift = 0.0;
      // double p = 0, q = 0, r = 0; // , s = 0, z = 0;
      int iter = 0;
      for(int n = nn - 1; n >= low;) {
        // Look for single small sub-diagonal element
        int l = n;
        for(; l > low; --l) {
          double s = abs(H[l - 1][l - 1]) + abs(H[l][l]);
          if(abs(H[l][l - 1]) < 0x1P-52 * (s == 0. ? norm : s)) {
            break;
          }
        }

        // Check for convergence
        if(l == n) {
          // One root found
          d[n] = (H[n][n] += exshift);
          e[n] = 0.;
          n--;
          iter = 0;
          continue;
        }
        final double[] Hn = H[n], Hnm1 = H[n - 1];
        if(l == n - 1) {
          // Two roots found
          double w = Hn[n - 1] * Hnm1[n];
          double p = (Hnm1[n - 1] - Hn[n]) * 0.5, q = p * p + w;
          double x = (Hn[n] += exshift), z = FastMath.sqrt(abs(q));
          Hnm1[n - 1] += exshift;

          if(q >= 0) {
            // Real pair
            z = (p >= 0) ? p + z : p - z;
            d[n - 1] = x + z;
            d[n] = (z != 0.) ? x - w / z : d[n - 1];
            e[n] = e[n - 1] = 0.;
            x = Hn[n - 1];
            double s = abs(x) + abs(z);
            p = x / s;
            q = z / s;
            double r = FastMath.hypot(p, q);
            p /= r;
            q /= r;

            // Row modification
            for(int j = n - 1; j < nn; j++) {
              double tmp = Hnm1[j];
              Hnm1[j] = q * tmp + p * Hn[j];
              Hn[j] = q * Hn[j] - p * tmp;
            }

            // Column modification
            for(int i = 0; i <= n; i++) {
              modifyQP(H[i], n, p, q);
            }

            // Accumulate transformations
            for(int i = low; i <= high; i++) {
              modifyQP(V[i], n, p, q);
            }
          }
          else {
            // Complex pair
            d[n] = d[n - 1] = x + p;
            e[n] = -(e[n - 1] = z);
          }
          n -= 2;
          iter = 0;
          continue;
        }
        // No convergence yet

        // Form shift
        double x = Hn[n], y = 0., w = 0.;
        y = Hnm1[n - 1];
        w = Hn[n - 1] * Hnm1[n];

        // Wilkinson's original ad hoc shift
        if(iter == 10) {
          exshift += x;
          for(int i = low; i <= n; i++) {
            H[i][i] -= x;
          }
          double s = abs(Hn[n - 1]) + abs(Hnm1[n - 2]);
          x = y = 0.75 * s;
          w = -0.4375 * s * s;
        }

        // MATLAB's new ad hoc shift
        if(iter == 30) {
          double s = (y - x) * 0.5;
          s = s * s + w;
          if(s > 0) {
            s = FastMath.sqrt(s);
            s = (y < x) ? -s : s;
            s = x - w / ((y - x) * 0.5 + s);
            for(int i = low; i <= n; i++) {
              H[i][i] -= s;
            }
            exshift += s;
            x = y = w = 0.964;
          }
        }

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

        double p = 0., q = 0., r = 0.;
        // Look for two consecutive small sub-diagonal elements
        int m = n - 2;
        for(; m >= l; m--) {
          final double[] Hm = H[m], Hmp1 = H[m + 1];
          final double z = Hm[m], xz = x - z;
          p = (xz * (y - z) - w) / Hmp1[m] + Hm[m + 1];
          q = Hmp1[m + 1] - xz - y;
          r = H[m + 2][m + 1];
          double tmp = abs(p) + abs(q) + abs(r);
          p /= tmp;
          q /= tmp;
          r /= tmp;
          if(m == l || abs(Hm[m - 1]) * (abs(q) + abs(r)) < 0x1P-52 * (abs(p) * (abs(H[m - 1][m - 1]) + abs(z) + abs(Hmp1[m + 1])))) {
            break;
          }
        }

        for(int i = m + 2; i <= n; i++) {
          final double[] Hi = H[i];
          Hi[i - 2] = 0.;
          if(i > m + 2) {
            Hi[i - 3] = 0.;
          }
        }

        // Double QR step involving rows l:n and columns m:n
        for(int k = m, last = n - 1; k <= last; k++) {
          boolean notlast = (k != last);
          final double[] Hk = H[k], Hkp1 = H[k + 1],
              Hkp2 = notlast ? H[k + 2] : null;
          if(k != m) {
            p = Hk[k - 1];
            q = Hkp1[k - 1];
            r = notlast ? Hkp2[k - 1] : 0.;
            x = abs(p) + abs(q) + abs(r);
            if(x == 0.0) {
              continue; // Jama 1.0.3 fix
            }
            p /= x;
            q /= x;
            r /= x;
          }
          double s = FastMath.hypot(p, q, r);
          s = (p < 0) ? -s : s;
          if(s != 0) {
            if(k != m) {
              Hk[k - 1] = -s * x;
            }
            else if(l != m) {
              Hk[k - 1] = -Hk[k - 1];
            }
            p += s;
            x = p / s;
            y = q / s;
            double z = r / s;
            q /= p;
            r /= p;

            // Row modification
            for(int j = k; j < nn; j++) {
              double tmp = Hk[j] + q * Hkp1[j];
              if(notlast) {
                tmp += r * Hkp2[j];
                Hkp2[j] -= tmp * z;
              }
              Hk[j] -= tmp * x;
              Hkp1[j] -= tmp * y;
            }

            // Column modification
            for(int i = 0; i <= Math.min(n, k + 3); i++) {
              modifyQR(H[i], k, notlast, q, r, x, y, z);
            }

            // Accumulate transformations
            for(int i = low; i <= high; i++) {
              modifyQR(V[i], k, notlast, q, r, x, y, z);
            }
          } // (s != 0)
        } // k loop
        // check convergence
      } // while (n >= low)
    }

    if(norm == 0.0) {
      return;
    }

    // Backsubstitute to find vectors of upper triangular form
    for(int n = nn - 1; n >= 0; n--) {
      double p = d[n], q = e[n];
      if(q == 0) {
        hqr2BacksubstituteReal(n, p, norm);
      }
      else if(q < 0) {
        hqr2BacksubstituteComplex(n, p, q, norm);
      }
    }

    // Vectors of isolated roots
    // Only matters if the user can modify low, high:
    // for(int i = 0; i < nn; i++) {
    // if(i < low || i > high) { System.arraycopy(H[i], i, V[i], i, nn - i); }
    // }

    // Back transformation to get eigenvectors of original matrix
    hqr2BackTransformation(nn, low, high);
  }
946169
def eigendecompose(self, normalise=False):
        
        self.eigenvalues, self.eigenvectors = utils.eigendecompose(self.tensor,
                                                                   normalise)
        return self.eigenvalues, self.eigenvectors
56306
private boolean computeEigenValues() {
       // make a copy of the internal tridiagonal matrix data for later use
       diagSaved = helper.copyDiag(diagSaved);
       offSaved = helper.copyOff(offSaved);

       vector.setQ(null);
       vector.setFastEigenvalues(true);

       // extract the eigenvalues
       if( !vector.process(-1,null,null) )
           return false;

       // save a copy of them since this data structure will be recycled next
       values = helper.copyEigenvalues(values);
       return true;
   }
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);
            }
        }
    }
596760
void constructLinearSystem(Planar<I> input , Se3_F32 g ) {
		int numBands = imageType.getNumBands();

		// first precompute everything that does not depend on pixel values
		inboundsPixels = 0;
		for (int i = 0; i < keypixels.size(); i++) {
			Pixel p = keypixels.data[i];

			// Apply the known warp
			SePointOps_F32.transform(g, p.p3, S);

			if( S.z <= 0 ) {
				p.valid = false;
				continue;
			}

			// Compute projected warped pixel coordinate on image I_1
			p.proj.x = (S.x / S.z) * fx + cx;
			p.proj.y = (S.y / S.z) * fy + cy;

			// make sure it's in the bounds
			if (p.proj.x < 0 || p.proj.x > input.width - 1 || p.proj.y < 0 || p.proj.y > input.height - 1) {
				p.valid = false;
				continue;
			} else {
				p.valid = true;
			}
			inboundsPixels++;

			// pi matrix derivative relative to t at S
			float ZZ = S.z * S.z;

			p.dP11 = fx / S.z;
			p.dP13 = -S.x * fx / ZZ;
			p.dP22 = fy / S.z;
			p.dP23 = -S.y * fy / ZZ;
		}

		// how compute the components which require
		errorOptical = 0;
		int row = 0;
		for (int band = 0; band < numBands; band++) {
			interpDX.setImage(derivX.getBand(band));
			interpDY.setImage(derivY.getBand(band));
			interpI.setImage(input.getBand(band));

			for (int i = 0; i < keypixels.size(); i++) {
				Pixel p = keypixels.data[i];

				if( !p.valid )
					continue;

				// Apply the known warp
				// TODO precompute?
				SePointOps_F32.transform(g, p.p3, S);

				// sample pixel values at warped location in I_1
				// NOTE: This could be highly optimized.  Compute and save interpolation weights once per input
				//       instead of for each band and image (current,dx,dy)
				// TODO create a special bilinear class for this?
				float current = interpI.get( p.proj.x, p.proj.y);
				float dx      = interpDX.get(p.proj.x, p.proj.y);
				float dy      = interpDY.get(p.proj.x, p.proj.y);

				// B = grad^T * dPI/dt = shape(1,3)
				float b1 = dx*p.dP11;
				float b2 = dy*p.dP22;
				float b3 = dx*p.dP13 + dy*p.dP23;

				// C * A(S'(x)) = shape(1,6)
				int indexA = row * 6;
				A.data[indexA++] = -b2*S.z + b3*S.y;
				A.data[indexA++] =  b1*S.z - b3*S.x;
				A.data[indexA++] = -b1*S.y + b2*S.x;
				A.data[indexA++] = b1;
				A.data[indexA++] = b2;
				A.data[indexA]   = b3;

				float error = -(current - p.bands[band]);
				y.data[row] = error;

				errorOptical += Math.abs(error);
				row += 1;
			}
		}
		errorOptical /= row;

		A.numRows = row;
		y.numRows = row;
	}
767896
def eigenvalue_band_properties(self):
        
        vbm = -float("inf")
        vbm_kpoint = None
        cbm = float("inf")
        cbm_kpoint = None
        for spin, d in self.eigenvalues.items():
            for k, val in enumerate(d):
                for (eigenval, occu) in val:
                    if occu > self.occu_tol and eigenval > vbm:
                        vbm = eigenval
                        vbm_kpoint = k
                    elif occu <= self.occu_tol and eigenval < cbm:
                        cbm = eigenval
                        cbm_kpoint = k
        return max(cbm - vbm, 0), cbm, vbm, vbm_kpoint == cbm_kpoint
1131690
def get_principal_axes(self):
        
        # Perform eigendecomposition - returns in order P, B, T
        _ = self.eigendecompose(normalise=True)
        principal_axes = GCMTPrincipalAxes()
        # Eigenvalues
        principal_axes.p_axis = {'eigenvalue': self.eigenvalues[0]}
        principal_axes.b_axis = {'eigenvalue': self.eigenvalues[1]}
        principal_axes.t_axis = {'eigenvalue': self.eigenvalues[2]}
        # Eigen vectors
        # 1) P axis
        azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 0], True)
        principal_axes.p_axis['azimuth'] = azim
        principal_axes.p_axis['plunge'] = plun
        # 2) B axis
        azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 1], True)
        principal_axes.b_axis['azimuth'] = azim
        principal_axes.b_axis['plunge'] = plun
        # 3) T axis
        azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 2], True)
        principal_axes.t_axis['azimuth'] = azim
        principal_axes.t_axis['plunge'] = plun
        return principal_axes
619428
private void hqr2 (double[][] V, double[][] H) {
   
      //  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 = this.n;
      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,q=0,r=0,s=0,z=0,t,w,x,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[i] = H[i][i];
            e[i] = 0.0;
         }
         for (int j = Math.max(i-1,0); j < nn; j++) {
            norm = norm + Math.abs(H[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[l-1][l-1]) + Math.abs(H[l][l]);
            if (s == 0.0) {
               s = norm;
            }
            if (Math.abs(H[l][l-1]) < eps * s) {
               break;
            }
            l--;
         }
       
         // Check for convergence
         // One root found
   
         if (l == n) {
            H[n][n] = H[n][n] + exshift;
            d[n] = H[n][n];
            e[n] = 0.0;
            n--;
            iter = 0;
   
         // Two roots found
   
         } else if (l == n-1) {
            w = H[n][n-1] * H[n-1][n];
            p = (H[n-1][n-1] - H[n][n]) / 2.0;
            q = p * p + w;
            z = Math.sqrt(Math.abs(q));
            H[n][n] = H[n][n] + exshift;
            H[n-1][n-1] = H[n-1][n-1] + exshift;
            x = H[n][n];
   
            // Real pair
   
            if (q >= 0) {
               if (p >= 0) {
                  z = p + z;
               } else {
                  z = p - z;
               }
               d[n-1] = x + z;
               d[n] = d[n-1];
               if (z != 0.0) {
                  d[n] = x - w / z;
               }
               e[n-1] = 0.0;
               e[n] = 0.0;
               x = H[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[n-1][j];
                  H[n-1][j] = q * z + p * H[n][j];
                  H[n][j] = q * H[n][j] - p * z;
               }
   
               // Column modification
   
               for (int i = 0; i <= n; i++) {
                  z = H[i][n-1];
                  H[i][n-1] = q * z + p * H[i][n];
                  H[i][n] = q * H[i][n] - p * z;
               }
   
               // Accumulate transformations
   
               for (int i = low; i <= high; i++) {
                  z = V[i][n-1];
                  V[i][n-1] = q * z + p * V[i][n];
                  V[i][n] = q * V[i][n] - p * z;
               }
   
            // Complex pair
   
            } else {
               d[n-1] = x + p;
               d[n] = x + p;
               e[n-1] = z;
               e[n] = -z;
            }
            n = n - 2;
            iter = 0;
   
         // No convergence yet
   
         } else {
   
            // Form shift
   
            x = H[n][n];
            y = 0.0;
            w = 0.0;
            if (l < n) {
               y = H[n-1][n-1];
               w = H[n][n-1] * H[n-1][n];
            }
   
            // Wilkinson's original ad hoc shift
   
            if (iter == 10) {
               exshift += x;
               for (int i = low; i <= n; i++) {
                  H[i][i] -= x;
               }
               s = Math.abs(H[n][n-1]) + Math.abs(H[n-1][n-2]);
               x = y = 0.75 * s;
               w = -0.4375 * s * s;
            }

            // 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[i][i] -= 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[m][m];
               r = x - z;
               s = y - z;
               p = (r * s - w) / H[m+1][m] + H[m][m+1];
               q = H[m+1][m+1] - z - r - s;
               r = H[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[m][m-1]) * (Math.abs(q) + Math.abs(r)) <
                  eps * (Math.abs(p) * (Math.abs(H[m-1][m-1]) + Math.abs(z) +
                  Math.abs(H[m+1][m+1])))) {
                     break;
               }
               m--;
            }
   
            for (int i = m+2; i <= n; i++) {
               H[i][i-2] = 0.0;
               if (i > m+2) {
                  H[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[k][k-1];
                  q = H[k+1][k-1];
                  r = (notlast ? H[k+2][k-1] : 0.0);
                  x = Math.abs(p) + Math.abs(q) + Math.abs(r);
                  if (x != 0.0) {
                     p = p / x;
                     q = q / x;
                     r = r / x;
                  }
               }
               if (x == 0.0) {
                  break;
               }
               s = Math.sqrt(p * p + q * q + r * r);
               if (p < 0) {
                  s = -s;
               }
               if (s != 0) {
                  if (k != m) {
                     H[k][k-1] = -s * x;
                  } else if (l != m) {
                     H[k][k-1] = -H[k][k-1];
                  }
                  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[k][j] + q * H[k+1][j];
                     if (notlast) {
                        p = p + r * H[k+2][j];
                        H[k+2][j] = H[k+2][j] - p * z;
                     }
                     H[k][j] = H[k][j] - p * x;
                     H[k+1][j] = H[k+1][j] - p * y;
                  }
   
                  // Column modification
   
                  for (int i = 0; i <= Math.min(n,k+3); i++) {
                     p = x * H[i][k] + y * H[i][k+1];
                     if (notlast) {
                        p = p + z * H[i][k+2];
                        H[i][k+2] = H[i][k+2] - p * r;
                     }
                     H[i][k] = H[i][k] - p;
                     H[i][k+1] = H[i][k+1] - p * q;
                  }
   
                  // Accumulate transformations
   
                  for (int i = low; i <= high; i++) {
                     p = x * V[i][k] + y * V[i][k+1];
                     if (notlast) {
                        p = p + z * V[i][k+2];
                        V[i][k+2] = V[i][k+2] - p * r;
                     }
                     V[i][k] = V[i][k] - p;
                     V[i][k+1] = V[i][k+1] - 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[n];
         q = e[n];
   
         // Real vector
   
         if (q == 0) {
            int l = n;
            H[n][n] = 1.0;
            for (int i = n-1; i >= 0; i--) {
               w = H[i][i] - p;
               r = 0.0;
               for (int j = l; j <= n; j++) {
                  r = r + H[i][j] * H[j][n];
               }
               if (e[i] < 0.0) {
                  z = w;
                  s = r;
               } else {
                  l = i;
                  if (e[i] == 0.0) {
                     if (w != 0.0) {
                        H[i][n] = -r / w;
                     } else {
                        H[i][n] = -r / (eps * norm);
                     }
   
                  // Solve real equations
   
                  } else {
                     x = H[i][i+1];
                     y = H[i+1][i];
                     q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
                     t = (x * s - z * r) / q;
                     H[i][n] = t;
                     if (Math.abs(x) > Math.abs(z)) {
                        H[i+1][n] = (-r - w * t) / x;
                     } else {
                        H[i+1][n] = (-s - y * t) / z;
                     }
                  }
   
                  // Overflow control
   
                  t = Math.abs(H[i][n]);
                  if ((eps * t) * t > 1) {
                     for (int j = i; j <= n; j++) {
                        H[j][n] = H[j][n] / t;
                     }
                  }
               }
            }
   
         // Complex vector
   
         } else if (q < 0) {
            int l = n-1;

            // Last vector component imaginary so matrix is triangular
   
            if (Math.abs(H[n][n-1]) > Math.abs(H[n-1][n])) {
               H[n-1][n-1] = q / H[n][n-1];
               H[n-1][n] = -(H[n][n] - p) / H[n][n-1];
            } else {
               cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q);
               H[n-1][n-1] = cdivr;
               H[n-1][n] = cdivi;
            }
            H[n][n-1] = 0.0;
            H[n][n] = 1.0;
            for (int i = n-2; i >= 0; i--) {
               double ra,sa,vr,vi;
               ra = 0.0;
               sa = 0.0;
               for (int j = l; j <= n; j++) {
                  ra = ra + H[i][j] * H[j][n-1];
                  sa = sa + H[i][j] * H[j][n];
               }
               w = H[i][i] - p;
   
               if (e[i] < 0.0) {
                  z = w;
                  r = ra;
                  s = sa;
               } else {
                  l = i;
                  if (e[i] == 0) {
                     cdiv(-ra,-sa,w,q);
                     H[i][n-1] = cdivr;
                     H[i][n] = cdivi;
                  } else {
   
                     // Solve complex equations
   
                     x = H[i][i+1];
                     y = H[i+1][i];
                     vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
                     vi = (d[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));
                     }
                     cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
                     H[i][n-1] = cdivr;
                     H[i][n] = cdivi;
                     if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
                        H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x;
                        H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x;
                     } else {
                        cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q);
                        H[i+1][n-1] = cdivr;
                        H[i+1][n] = cdivi;
                     }
                  }
   
                  // Overflow control

                  t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n]));
                  if ((eps * t) * t > 1) {
                     for (int j = i; j <= n; j++) {
                        H[j][n-1] = H[j][n-1] / t;
                        H[j][n] = H[j][n] / 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[i][j] = H[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[i][k] * H[k][j];
            }
            V[i][j] = z;
         }
      }
   }
1111518
private Set<HiCSSubspace> calculateSubspaces(Relation<? extends NumberVector> relation, ArrayList<ArrayDBIDs> subspaceIndex, Random random) {
    final int dbdim = RelationUtil.dimensionality(relation);

    FiniteProgress dprog = LOG.isVerbose() ? new FiniteProgress("Subspace dimensionality", dbdim, LOG) : null;
    if(dprog != null) {
      dprog.setProcessed(2, LOG);
    }

    TreeSet<HiCSSubspace> subspaceList = new TreeSet<>(HiCSSubspace.SORT_BY_SUBSPACE);
    TopBoundedHeap<HiCSSubspace> dDimensionalList = new TopBoundedHeap<>(cutoff, HiCSSubspace.SORT_BY_CONTRAST_ASC);
    FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Generating two-element subsets", (dbdim * (dbdim - 1)) >> 1, LOG) : null;
    // compute two-element sets of subspaces
    for(int i = 0; i < dbdim; i++) {
      for(int j = i + 1; j < dbdim; j++) {
        HiCSSubspace ts = new HiCSSubspace();
        ts.set(i);
        ts.set(j);
        calculateContrast(relation, ts, subspaceIndex, random);
        dDimensionalList.add(ts);
        LOG.incrementProcessed(prog);
      }
    }
    LOG.ensureCompleted(prog);

    IndefiniteProgress qprog = LOG.isVerbose() ? new IndefiniteProgress("Testing subspace candidates", LOG) : null;
    for(int d = 3; !dDimensionalList.isEmpty(); d++) {
      if(dprog != null) {
        dprog.setProcessed(d, LOG);
      }
      // result now contains all d-dimensional sets of subspaces

      ArrayList<HiCSSubspace> candidateList = new ArrayList<>(dDimensionalList.size());
      for(Heap<HiCSSubspace>.UnorderedIter it = dDimensionalList.unorderedIter(); it.valid(); it.advance()) {
        subspaceList.add(it.get());
        candidateList.add(it.get());
      }
      dDimensionalList.clear();
      // candidateList now contains the *m* best d-dimensional sets
      Collections.sort(candidateList, HiCSSubspace.SORT_BY_SUBSPACE);

      // TODO: optimize APRIORI style, by not even computing the bit set or?
      for(int i = 0; i < candidateList.size() - 1; i++) {
        for(int j = i + 1; j < candidateList.size(); j++) {
          HiCSSubspace set1 = candidateList.get(i), set2 = candidateList.get(j);

          HiCSSubspace joinedSet = new HiCSSubspace();
          joinedSet.or(set1);
          joinedSet.or(set2);
          if(joinedSet.cardinality() != d) {
            continue;
          }

          calculateContrast(relation, joinedSet, subspaceIndex, random);
          dDimensionalList.add(joinedSet);
          LOG.incrementProcessed(qprog);
        }
      }
      // Prune
      for(HiCSSubspace cand : candidateList) {
        for(Heap<HiCSSubspace>.UnorderedIter it = dDimensionalList.unorderedIter(); it.valid(); it.advance()) {
          if(it.get().contrast > cand.contrast) {
            subspaceList.remove(cand);
            break;
          }
        }
      }
    }
    LOG.setCompleted(qprog);
    if(dprog != null) {
      dprog.setProcessed(dbdim, LOG);
      dprog.ensureCompleted(LOG);
    }
    return subspaceList;
  }