Compute H (observed perplexity) for row i, and the row pij_i.

@param ignore Object to skip
@param it Distances list
@param p Output probabilities
@param mbeta {@code -1. / (2 * sigma * sigma)}
@return Observed perplexity
546968
protected static double computeH(DBIDRef ignore, DoubleDBIDListIter it, double[] p, double mbeta) {
    double sumP = 0.;
    // Skip point "i", break loop in two:
    it.seek(0);
    for(int j = 0; it.valid(); j++, it.advance()) {
      if(DBIDUtil.equal(ignore, it)) {
        p[j] = 0;
        continue;
      }
      sumP += (p[j] = FastMath.exp(it.doubleValue() * mbeta));
    }
    if(!(sumP > 0)) {
      // All pij are zero. Bad news.
      return Double.NEGATIVE_INFINITY;
    }
    final double s = 1. / sumP; // Scaling factor
    double sum = 0.;
    // While we could skip pi[i], it should be 0 anyway.
    it.seek(0);
    for(int j = 0; it.valid(); j++, it.advance()) {
      sum += it.doubleValue() * (p[j] *= s);
    }
    return FastMath.log(sumP) - mbeta * sum;
  }
400537
public static void nominateNeighbors(DBIDIter ignore, DBIDArrayIter di, double[] p, double norm, WritableDoubleDataStore scores) {
    for(di.seek(0); di.valid(); di.advance()) {
      if(DBIDUtil.equal(ignore, di)) {
        continue;
      }
      double v = p[di.getOffset()] * norm; // Normalize
      if(!(v > 0)) {
        break;
      }
      scores.putDouble(di, scores.doubleValue(di) * (1 - v));
    }
  }
957755
def _handle_outliers(self, p_o):
        
        if self.ignore_outliers:
            outliers = np.where(p_o.sum(axis=1)==0)[0]
            if outliers.size > 0:
                p_o[outliers, :] = 1.0
                self.found_outliers = True
        return p_o
230106
public static Measurement first(Iterable<Measurement> ms, Predicate<Measurement> p) {
    Iterator<Measurement> it = filter(ms, p).iterator();
    return it.hasNext() ? it.next() : null;
  }
1304808
private SimpleMatrix computeAffineH( List<AssociatedPair> observations ,
										 DMatrixRMaj H , DMatrixRMaj Hzero ) {

		SimpleMatrix A = new SimpleMatrix(observations.size(),3);
		SimpleMatrix b = new SimpleMatrix(A.numRows(),1);

		Point2D_F64 c = new Point2D_F64();
		Point2D_F64 k = new Point2D_F64();

		for( int i = 0; i < observations.size(); i++ ) {
			AssociatedPair a = observations.get(i);

			GeometryMath_F64.mult(Hzero, a.p1, k);
			GeometryMath_F64.mult(H,a.p2,c);

			A.setRow(i,0,k.x,k.y,1);
			b.set(i,0,c.x);
		}

		SimpleMatrix x = A.solve(b);

		SimpleMatrix Ha = SimpleMatrix.identity(3);
		Ha.setRow(0,0,x.getDDRM().data);

		return Ha;
	}
104758
public static void fe_neg(int[] h,int[] f)
{
  int f0 = f[0];
  int f1 = f[1];
  int f2 = f[2];
  int f3 = f[3];
  int f4 = f[4];
  int f5 = f[5];
  int f6 = f[6];
  int f7 = f[7];
  int f8 = f[8];
  int f9 = f[9];
  int h0 = -f0;
  int h1 = -f1;
  int h2 = -f2;
  int h3 = -f3;
  int h4 = -f4;
  int h5 = -f5;
  int h6 = -f6;
  int h7 = -f7;
  int h8 = -f8;
  int h9 = -f9;
  h[0] = (int)h0;
  h[1] = (int)h1;
  h[2] = (int)h2;
  h[3] = (int)h3;
  h[4] = (int)h4;
  h[5] = (int)h5;
  h[6] = (int)h6;
  h[7] = (int)h7;
  h[8] = (int)h8;
  h[9] = (int)h9;
}
110309
protected void add( int which , Homography2D_F64 H , DMatrixRMaj A ) {
		int idx = which*6*6;

		// Row 0
		A.data[idx++] = H.a11*H.a11-1;
		A.data[idx++] = 2*H.a11*H.a12;
		A.data[idx++] = 2*H.a11*H.a13;
		A.data[idx++] = H.a12*H.a12;
		A.data[idx++] = 2*H.a12*H.a13;
		A.data[idx++] = H.a13*H.a13;
		// Row 1
		A.data[idx++] = H.a11*H.a21;
		A.data[idx++] = H.a12*H.a21 + H.a11*H.a22-1;
		A.data[idx++] = H.a13*H.a21 + H.a11*H.a23;
		A.data[idx++] = H.a12*H.a22;
		A.data[idx++] = H.a13*H.a22 + H.a12*H.a23;
		A.data[idx++] = H.a13*H.a23;
		// Row 2
		A.data[idx++] = H.a11*H.a31;
		A.data[idx++] = H.a12*H.a31 + H.a11*H.a32;
		A.data[idx++] = H.a13*H.a31 + H.a11*H.a33-1;
		A.data[idx++] = H.a12*H.a32;
		A.data[idx++] = H.a13*H.a32 + H.a12*H.a33;
		A.data[idx++] = H.a13*H.a33;
		// Row 3
		A.data[idx++] = H.a21*H.a21;
		A.data[idx++] = 2*H.a21*H.a22;
		A.data[idx++] = 2*H.a21*H.a23;
		A.data[idx++] = H.a22*H.a22 - 1;
		A.data[idx++] = 2*H.a22*H.a23;
		A.data[idx++] = H.a23*H.a23;
		// Row 4
		A.data[idx++] = H.a21*H.a31;
		A.data[idx++] = H.a22*H.a31 + H.a21*H.a32;
		A.data[idx++] = H.a23*H.a31 + H.a21*H.a33;
		A.data[idx++] = H.a22*H.a32;
		A.data[idx++] = H.a23*H.a32 + H.a22*H.a33-1;
		A.data[idx++] = H.a23*H.a33;
		// Row 5
		A.data[idx++] = H.a31*H.a31;
		A.data[idx++] = 2*H.a31*H.a32;
		A.data[idx++] = 2*H.a31*H.a33;
		A.data[idx++] = H.a32*H.a32;
		A.data[idx++] = 2*H.a32*H.a33;
		A.data[idx  ] = H.a33*H.a33-1;
	}
538390
private INDArray x2p(final INDArray X, double tolerance, double perplexity) {
        int n = X.rows();
        final INDArray p = zeros(n, n);
        final INDArray beta = ones(n, 1);
        final double logU = Math.log(perplexity);

        INDArray sumX = pow(X, 2).sum(1);

        logger.debug("sumX shape: " + Arrays.toString(sumX.shape()));

        INDArray times = X.mmul(X.transpose()).muli(-2);

        logger.debug("times shape: " + Arrays.toString(times.shape()));

        INDArray prodSum = times.transpose().addiColumnVector(sumX);

        logger.debug("prodSum shape: " + Arrays.toString(prodSum.shape()));

        INDArray D = X.mmul(X.transpose()).mul(-2) // thats times
                        .transpose().addColumnVector(sumX) // thats prodSum
                        .addRowVector(sumX.transpose()); // thats D

        logger.info("Calculating probabilities of data similarities...");
        logger.debug("Tolerance: " + tolerance);
        for (int i = 0; i < n; i++) {
            if (i % 500 == 0 && i > 0)
                logger.info("Handled [" + i + "] records out of [" + n + "]");

            double betaMin = Double.NEGATIVE_INFINITY;
            double betaMax = Double.POSITIVE_INFINITY;
            int[] vals = Ints.concat(ArrayUtil.range(0, i), ArrayUtil.range(i + 1, n));
            INDArrayIndex[] range = new INDArrayIndex[] {new SpecifiedIndex(vals)};

            INDArray row = D.slice(i).get(range);
            Pair<Double, INDArray> pair = hBeta(row, beta.getDouble(i));
            //INDArray hDiff = pair.getFirst().sub(logU);
            double hDiff = pair.getFirst() - logU;
            int tries = 0;

            //while hdiff > tolerance
            while (Math.abs(hDiff) > tolerance && tries < 50) {
                //if hdiff > 0
                if (hDiff > 0) {
                    betaMin = beta.getDouble(i);
                    if (Double.isInfinite(betaMax))
                        beta.putScalar(i, beta.getDouble(i) * 2.0);
                    else
                        beta.putScalar(i, (beta.getDouble(i) + betaMax) / 2.0);
                } else {
                    betaMax = beta.getDouble(i);
                    if (Double.isInfinite(betaMin))
                        beta.putScalar(i, beta.getDouble(i) / 2.0);
                    else
                        beta.putScalar(i, (beta.getDouble(i) + betaMin) / 2.0);
                }

                pair = hBeta(row, beta.getDouble(i));
                hDiff = pair.getFirst() - logU;
                tries++;
            }
            p.slice(i).put(range, pair.getSecond());
        }


        //dont need data in memory after
        logger.info("Mean value of sigma " + sqrt(beta.rdiv(1)).mean(Integer.MAX_VALUE));
        BooleanIndexing.replaceWhere(p, 1e-12, Conditions.isNan());

        //set 0 along the diagonal
        INDArray permute = p.transpose();

        INDArray pOut = p.add(permute);

        pOut.divi(pOut.sumNumber().doubleValue() + 1e-6);

        pOut.muli(4);

        BooleanIndexing.replaceWhere(pOut, 1e-12, Conditions.lessThan(1e-12));
        //ensure no nans

        return pOut;
    }
544348
@Override
    public boolean isOrderDeterministic() {
        if ( ! super.isOrderDeterministic()) {
            return false;
        }
        IndexScanPlanNode index_scan = getInlineIndexScan();
        if ( ! index_scan.isOrderDeterministic()) {
            m_nondeterminismDetail = index_scan.m_nondeterminismDetail;
            return false;
        }
        return true;
    }
731834
public boolean computeRectifyingHomography( DMatrixRMaj H ) {
		H.reshape(4,4);

		// insert the results into H
		// H = [K 0;-p'*K 1 ]
		H.zero();
		for (int i = 0; i < 3; i++) {
			for (int j = i; j < 3; j++) {
				H.set(i,j,k.get(i,j));
			}
		}
		// p and k have different scales, fix that
		H.set(3,0, -(p.a1*k.a11 + p.a2*k.a21 + p.a3*k.a31));
		H.set(3,1, -(p.a1*k.a12 + p.a2*k.a22 + p.a3*k.a32));
		H.set(3,2, -(p.a1*k.a13 + p.a2*k.a23 + p.a3*k.a33));
		H.set(3,3,1);

		return true;
	}
584080
public static Iterable<Measurement> filter(
      final Iterable<Measurement> ms, final Predicate<Measurement> p) {
    return () -> new FilteredIterator<>(ms.iterator(), p);
  }
370157
private static void hqr(DenseMatrix A, double[] d, double[] e) {
        int n = A.nrows();
        int nn, m, l, k, j, its, i, mmin;
        double z, y, x, w, v, u, t, s, r = 0.0, q = 0.0, p = 0.0, anorm = 0.0;

        for (i = 0; i < n; i++) {
            for (j = Math.max(i - 1, 0); j < n; j++) {
                anorm += Math.abs(A.get(i, j));
            }
        }
        nn = n - 1;
        t = 0.0;
        while (nn >= 0) {
            its = 0;
            do {
                for (l = nn; l > 0; l--) {
                    s = Math.abs(A.get(l - 1, l - 1) + Math.abs(A.get(l, l)));
                    if (s == 0.0) {
                        s = anorm;
                    }
                    if (Math.abs(A.get(l, l - 1)) <= Math.EPSILON * s) {
                        A.set(l, l - 1, 0.0);
                        break;
                    }
                }
                x = A.get(nn, nn);
                if (l == nn) {
                    d[nn--] = x + t;
                } else {
                    y = A.get(nn - 1, nn - 1);
                    w = A.get(nn, nn - 1) * A.get(nn - 1, nn);
                    if (l == nn - 1) {
                        p = 0.5 * (y - x);
                        q = p * p + w;
                        z = Math.sqrt(Math.abs(q));
                        x += t;
                        if (q >= 0.0) {
                            z = p + Math.copySign(z, p);
                            d[nn - 1] = d[nn] = x + z;
                            if (z != 0.0) {
                                d[nn] = x - w / z;
                            }
                        } else {
                            d[nn] = x + p;
                            e[nn] = -z;
                            d[nn - 1] = d[nn];
                            e[nn - 1] = -e[nn];
                        }
                        nn -= 2;
                    } else {
                        if (its == 30) {
                            throw new IllegalStateException("Too many iterations in hqr");
                        }
                        if (its == 10 || its == 20) {
                            t += x;
                            for (i = 0; i < nn + 1; i++) {
                                A.sub(i, i, x);
                            }
                            s = Math.abs(A.get(nn, nn - 1)) + Math.abs(A.get(nn - 1, nn - 2));
                            y = x = 0.75 * s;
                            w = -0.4375 * s * s;
                        }
                        ++its;
                        for (m = nn - 2; m >= l; m--) {
                            z = A.get(m, m);
                            r = x - z;
                            s = y - z;
                            p = (r * s - w) / A.get(m + 1, m) + A.get(m, m + 1);
                            q = A.get(m + 1, m + 1) - z - r - s;
                            r = A.get(m + 2, m + 1);
                            s = Math.abs(p) + Math.abs(q) + Math.abs(r);
                            p /= s;
                            q /= s;
                            r /= s;
                            if (m == l) {
                                break;
                            }
                            u = Math.abs(A.get(m, m - 1)) * (Math.abs(q) + Math.abs(r));
                            v = Math.abs(p) * (Math.abs(A.get(m - 1, m - 1)) + Math.abs(z) + Math.abs(A.get(m + 1, m + 1)));
                            if (u <= Math.EPSILON * v) {
                                break;
                            }
                        }
                        for (i = m; i < nn - 1; i++) {
                            A.set(i + 2, i, 0.0);
                            if (i != m) {
                                A.set(i + 2, i - 1, 0.0);
                            }
                        }
                        for (k = m; k < nn; k++) {
                            if (k != m) {
                                p = A.get(k, k - 1);
                                q = A.get(k + 1, k - 1);
                                r = 0.0;
                                if (k + 1 != nn) {
                                    r = A.get(k + 2, k - 1);
                                }
                                if ((x = Math.abs(p) + Math.abs(q) + Math.abs(r)) != 0.0) {
                                    p /= x;
                                    q /= x;
                                    r /= x;
                                }
                            }
                            if ((s = Math.copySign(Math.sqrt(p * p + q * q + r * r), p)) != 0.0) {
                                if (k == m) {
                                    if (l != m) {
                                        A.set(k, k - 1, -A.get(k, k - 1));
                                    }
                                } else {
                                    A.set(k, k - 1, -s * x);
                                }
                                p += s;
                                x = p / s;
                                y = q / s;
                                z = r / s;
                                q /= p;
                                r /= p;
                                for (j = k; j < nn + 1; j++) {
                                    p = A.get(k, j) + q * A.get(k + 1, j);
                                    if (k + 1 != nn) {
                                        p += r * A.get(k + 2, j);
                                        A.sub(k + 2, j, p * z);
                                    }
                                    A.sub(k + 1, j, p * y);
                                    A.sub(k, j, p * x);
                                }
                                mmin = nn < k + 3 ? nn : k + 3;
                                for (i = l; i < mmin + 1; i++) {
                                    p = x * A.get(i, k) + y * A.get(i, k + 1);
                                    if (k + 1 != nn) {
                                        p += z * A.get(i, k + 2);
                                        A.sub(i, k + 2, p * r);
                                    }
                                    A.sub(i, k + 1, p * q);
                                    A.sub(i, k, p);
                                }
                            }
                        }
                    }
                }
            } while (l + 1 < nn);
        }
    }
1151260
public double get2DHV(List<? extends Solution<?>> solutionSet) {
    double hv = 0.0;
    if (solutionSet.size() > 0) {
      hv = Math.abs((solutionSet.get(0).getObjective(0) - referencePoint.getValue(0)) *
          (solutionSet.get(0).getObjective(1) - referencePoint.getValue(1)));

      for (int i = 1; i < solutionSet.size(); i++) {
        double tmp ;
        tmp = Math.abs((solutionSet.get(i).getObjective(0) - referencePoint.getValue(0)) * (solutionSet.get(i).getObjective(1) - solutionSet.get(i - 1).getObjective(1)));
        hv += tmp;
      }
    }
    return hv;
  }
275311
private void pldltptBK() throws Exception {
		//matrix S will be changed by the factorization, ad we do not want to change the matrix passed in by the client
		DoubleMatrix2D S = (rescaler==null)? this.Q.copy() : this.Q;
		int n = S.rows();
		DoubleMatrix2D A = S.copy();
		this.P = DoubleFactory2D.sparse.identity(n);
		this.D = DoubleFactory2D.sparse.make(n, n);
		this.L = DoubleFactory2D.sparse.make(n, n);
		DoubleMatrix2D LT = ALG.transpose(L);//@TODO: remove and work only with L
		double mu = (1 + Math.sqrt(17)) / 8;
		int s = 0;
		for (int j = 0; j < n; j++) {
			//log.debug("j: " + j);
			DoubleMatrix2D LPart = null;
			DoubleMatrix2D LTPart = null;
			DoubleMatrix2D APart = null;
			DoubleMatrix2D DPart = null;
			
			//Determine largest off-diagonal element in column j
			double lambda = -1;
			int r = -1;
			for(int rr = j+1; rr<n; rr++){
				double d = Math.abs(A.getQuick(rr, j)); 
				if(d > lambda){
					lambda = d;
					r = rr;
				}
			}
			double ajj = A.getQuick(j, j);
			if (Math.abs(ajj) > mu*lambda) {
				// Choose a 1 x 1 pivot,
				D.setQuick(j, j, ajj);
				s = 1;
			} else {
				//Determine largest off-diagonal element in row/column r
				double sigma = -1;
				for(int rr = 0; rr<n; rr++){
					if(rr==r){
						continue;
					}
					double d = Math.abs(A.getQuick(rr, r)); 
					if(d > sigma){
						sigma = d;
					}
				}
				if(sigma * Math.abs(ajj) > mu*lambda*lambda){
					//Choose a 1 x 1 pivot,
					D.setQuick(j, j, A.getQuick(j, j));
					s = 1;
				}else if(Math.abs(A.getQuick(r, r)) > mu*sigma){
					//Symmetrically permute row/column r to position j
				  A = ColtUtils.symmPermutation(A, r, j);
					P.setQuick(r, r, 0);
					P.setQuick(j, j, 0);
					P.setQuick(r, j, 1);
					P.setQuick(j, r, 1);
					//Choose a 1 x 1 pivot,
					D.setQuick(j, j, A.getQuick(j, j));
					s = 1;
				}else{
					//Symmetrically permute rows/columns r and p to positions j and j + 1 respectively
					int t = j+1;
					//int v = j+1;
					A = ColtUtils.symmPermutation(A, r, t);
					P.setQuick(r, r, 0);
					P.setQuick(t, t, 0);
					P.setQuick(r, t, 1);
					P.setQuick(t, r, 1);
//					A = Utils.symmPermutation(A, p, v);
//					P.setQuick(p, p, 0);
//					P.setQuick(v, v, 0);
//					P.setQuick(p, v, 1);
//					P.setQuick(v, p, 1);
					//Choose a 2 x pivot,
					//D(j:j+1)(j:j+1) =  A(j:j+1)(j:j+1)
					D.setQuick(j, j, A.getQuick(j, j));
					D.setQuick(j, j+1, A.getQuick(j, j+1));
					D.setQuick(j+1, j, A.getQuick(j+1, j));
					D.setQuick(j+1, j+1, A.getQuick(j+1, j+1));
					s = 2;
				}
			}
			//log.debug("s: " + s);

			// L(j:n)(j:j+s-1) = A(j:n)(j:j+s-1).DInv(j:j+s-1)(j:j+s-1)
			APart = A.viewPart(j, j, n - j, 1 + s - 1);
//			if(s==2){
//				try{
//					ColtUtils.invert2x2Matrix(D.viewPart(j, j, 1 + s - 1, 1 + s - 1));
//				}catch(Exception e){
//					throw e;
//				}
//			}
			DPart = ALG.inverse(D.viewPart(j, j, 1 + s - 1, 1 + s - 1));
			LPart = L.viewPart(j, j, n - j, 1 + s - 1);
			DoubleMatrix2D AD = ALG.mult(APart, DPart);
			for (int rr = 0; rr < LPart.rows(); rr++) {
				for (int cc = 0; cc < LPart.columns(); cc++) {
					LPart.setQuick(rr, cc, AD.getQuick(rr, cc));
				}
			}

			// A(j+s-1:n)(j+s-1:n) = A(j+s-1:n)(j+s-1:n) -
			// L(j+s-1:n)(j:j+s-1).D(j:j+s-1)(j:j+s-1).LT(j:j+s-1)(j+s-1:n)
			LPart = L.viewPart(j + s - 1, j, n - (j + s - 1), 1 + s - 1);
			DPart = D.viewPart(j, j, 1 + s - 1, 1 + s - 1);
			LTPart = LT.viewPart(j, j + s - 1, s, n - (j + s - 1));
			APart = A.viewPart(j + s - 1, j + s - 1, n - (j + s - 1), n - (j + s - 1));
			DoubleMatrix2D LDLT = ALG.mult(LPart, ALG.mult(DPart, LTPart));
			for (int rr = 0; rr < APart.rows(); rr++) {
				for (int cc = 0; cc < APart.columns(); cc++) {
					APart.setQuick(rr, cc, APart.getQuick(rr, cc) - LDLT.getQuick(rr, cc));
				}
			}
			//logger.debug("L: " + ArrayUtils.toString(L.toArray()));
			//logger.debug("A: " + ArrayUtils.toString(A.toArray()));

			j = j + s - 1;
		}
	}
231068
public static DMatrixSparseCSC permutationMatrix(int[] p, boolean inverse, int N, DMatrixSparseCSC P) {

        if( P == null )
            P = new DMatrixSparseCSC(N,N,N);
        else
            P.reshape(N,N,N);
        P.indicesSorted = true;
        P.nz_length = N;

        // each column should have one element inside of it
        if( !inverse ) {
            for (int i = 0; i < N; i++) {
                P.col_idx[i + 1] = i + 1;
                P.nz_rows[p[i]] = i;
                P.nz_values[i] = 1;
            }
        } else {
            for (int i = 0; i < N; i++) {
                P.col_idx[i + 1] = i + 1;
                P.nz_rows[i] = p[i];
                P.nz_values[i] = 1;
            }
        }

        return P;
    }
737142
private void primeProduct(List<InvPair> v, IAtomContainer atomContainer) {
        Iterator<InvPair> it = v.iterator();
        Iterator<IAtom> n;
        InvPair inv;
        IAtom a;
        long summ;
        while (it.hasNext()) {
            inv = (InvPair) it.next();
            List<IAtom> neighbour = atomContainer.getConnectedAtomsList(inv.getAtom());
            n = neighbour.iterator();
            summ = 1;
            while (n.hasNext()) {
                a = n.next();
                int next = ((InvPair) a.getProperty(InvPair.INVARIANCE_PAIR)).getPrime();
                summ = summ * next;
            }
            inv.setLast(inv.getCurr());
            inv.setCurr(summ);
        }
    }
1304155
public SortedClassProbability[] sortByDescendingClassProbability(BinomialModelPrediction p) {
    String[] domainValues = m.getDomainValues(m.getResponseIdx());
    double[] classProbabilities = p.classProbabilities;
    return sortByDescendingClassProbability(domainValues, classProbabilities);
  }