Solves the kriging system as a vectorized operation. This method
        can take a lot of memory for large grids and/or large datasets.
970884
def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids):
        """"""

        npt = bd.shape[0]
        n = self.X_ADJUSTED.shape[0]
        zero_index = None
        zero_value = False

        a_inv = scipy.linalg.inv(a)

        if np.any(np.absolute(bd) <= self.eps):
            zero_value = True
            zero_index = np.where(np.absolute(bd) <= self.eps)

        if self.UNBIAS:
            b = np.zeros((npt, n_withdrifts+1, 1))
        else:
            b = np.zeros((npt, n_withdrifts, 1))
        b[:, :n, 0] = - self.variogram_function(self.variogram_model_parameters, bd)
        if zero_value:
            b[zero_index[0], zero_index[1], 0] = 0.0

        i = n
        if self.regional_linear_drift:
            b[:, i, 0] = xyz[:, 2]
            i += 1
            b[:, i, 0] = xyz[:, 1]
            i += 1
            b[:, i, 0] = xyz[:, 0]
            i += 1
        if self.specified_drift:
            for spec_vals in spec_drift_grids:
                b[:, i, 0] = spec_vals.flatten()
                i += 1
        if self.functional_drift:
            for func in self.functional_drift_terms:
                b[:, i, 0] = func(xyz[:, 2], xyz[:, 1], xyz[:, 0])
                i += 1
        if i != n_withdrifts:
            warnings.warn("Error in setting up kriging system. "
                          "Kriging may fail.", RuntimeWarning)
        if self.UNBIAS:
            b[:, n_withdrifts, 0] = 1.0

        if (~mask).any():
            mask_b = np.repeat(mask[:, np.newaxis, np.newaxis],
                               n_withdrifts+1, axis=1)
            b = np.ma.array(b, mask=mask_b)

        if self.UNBIAS:
            x = np.dot(a_inv, b.reshape((npt, n_withdrifts+1)).T).reshape((1, n_withdrifts+1, npt)).T
        else:
            x = np.dot(a_inv, b.reshape((npt, n_withdrifts)).T).reshape((1, n_withdrifts, npt)).T
        kvalues = np.sum(x[:, :n, 0] * self.VALUES, axis=1)
        sigmasq = np.sum(x[:, :, 0] * -b[:, :, 0], axis=1)

        return kvalues, sigmasq
49472
def solve(self,problem):
        
    
        # Local vars
        norm2 = self.norm2
        norminf = self.norminf
        parameters = self.parameters
        
        # Parameters
        feastol = parameters['feastol']
        optol = parameters['optol']
        maxiter = parameters['maxiter']
        quiet = parameters['quiet']
        sigma = parameters['sigma']
        eps = parameters['eps']
        ls_maxiter = parameters['line_search_maxiter']

        # Problem
        problem = cast_problem(problem)
        self.problem = problem
        
        # Linsolver
        self.linsolver = new_linsolver(parameters['linsolver'],'symmetric')

        # Reset
        self.reset()

        # Checks
        if not np.all(problem.l <= problem.u):
            raise OptSolverError_NoInterior(self)

        # Constants
        self.A = problem.A
        self.AT = problem.A.T
        self.b = problem.b
        self.u = problem.u+feastol/10.
        self.l = problem.l-feastol/10.
        self.n = problem.get_num_primal_variables()
        self.m1 = problem.get_num_linear_equality_constraints()
        self.m2 = problem.get_num_nonlinear_equality_constraints()
        self.e = np.ones(self.n)
        self.I = eye(self.n,format='coo')
        self.Omm1 = coo_matrix((self.m1,self.m1))
        self.Omm2 = coo_matrix((self.m2,self.m2))

        # Initial primal
        if problem.x is None:
            self.x = (self.u + self.l)/2.
        else:
            self.x = np.maximum(np.minimum(problem.x,problem.u),problem.l)

        # Initial duals
        if problem.lam is None:
            self.lam = np.zeros(problem.get_num_linear_equality_constraints())
        else:
            self.lam = problem.lam.copy()
        if problem.nu is None:
            self.nu = np.zeros(problem.get_num_nonlinear_equality_constraints())
        else:
            self.nu = problem.nu.copy()
        self.mu = np.minimum(1./(self.u-self.x), 1.)
        self.pi = np.minimum(1./(self.x-self.l), 1.)

        # Init vector
        self.y = np.hstack((self.x,self.lam,self.nu,self.mu,self.pi))

        # Average violation of complementarity slackness
        self.eta_mu = (np.dot(self.mu,self.u-self.x)/self.x.size) if self.x.size else 0.
        self.eta_pi = (np.dot(self.pi,self.x-self.l)/self.x.size) if self.x.size else 0.

        # Objective scaling
        fdata = self.func(self.y)
        self.obj_sca = np.maximum(norminf(problem.gphi)/10.,1.)
        fdata = self.func(self.y)

        # Header
        if not quiet:
            print('\nSolver: inlp')
            print('------------')
                                   
        # Outer
        s = 0.
        self.k = 0
        while True:

            # Average violation of complementarity slackness
            self.eta_mu = (np.dot(self.mu,self.u-self.x)/self.x.size) if self.x.size else 0.
            self.eta_pi = (np.dot(self.pi,self.x-self.l)/self.x.size) if self.x.size else 0.
            
            # Init eval
            fdata = self.func(self.y)

            # Target
            tau = sigma*norminf(fdata.GradF)
           
            # Header
            if not quiet:
                if self.k > 0:
                    print('')
                print('{0:^3s}'.format('iter'),end=' ')
                print('{0:^9s}'.format('phi'),end=' ')
                print('{0:^9s}'.format('pres'),end=' ')
                print('{0:^9s}'.format('dres'),end=' ')
                print('{0:^9s}'.format('gmax'),end=' ')
                print('{0:^8s}'.format('cu'),end=' ')
                print('{0:^8s}'.format('cl'),end=' ')
                print('{0:^8s}'.format('alpha'))
 
            # Inner
            while True:
                
                # Eval
                fdata = self.func(self.y)
                pres = norminf(np.hstack((fdata.rp1,fdata.rp2)))
                dres = norminf(np.hstack((fdata.rd,fdata.ru,fdata.rl)))
                gmax = norminf(fdata.GradF)                
                compu = norminf(self.mu*(self.u-self.x))
                compl = norminf(self.pi*(self.x-self.l))
                phi = problem.phi
                
                # Show progress
                if not quiet:
                    print('{0:^3d}'.format(self.k),end=' ')
                    print('{0:^9.2e}'.format(phi),end=' ')
                    print('{0:^9.2e}'.format(pres),end=' ')
                    print('{0:^9.2e}'.format(dres),end=' ')
                    print('{0:^9.2e}'.format(gmax),end=' ')
                    print('{0:^8.1e}'.format(compu),end=' ')
                    print('{0:^8.1e}'.format(compl),end=' ')
                    print('{0:^8.1e}'.format(s))

                # Done
                if self.k > 0 and pres < feastol and dres < optol and sigma*np.maximum(self.eta_mu,self.eta_pi) < optol:
                    self.set_status(self.STATUS_SOLVED)
                    self.set_error_msg('')
                    return
                    
                # Done
                if gmax < tau: 
                    break

                # Done 
                if pres < feastol and dres < optol and np.maximum(compu,compl) < optol:
                    break

                # Maxiters
                if self.k >= maxiter:
                    raise OptSolverError_MaxIters(self)
                
                # Search direction
                ux = self.u-self.x
                xl = self.x-self.l
                D1 = spdiags(self.mu/ux,0,self.n,self.n,format='coo')
                D2 = spdiags(self.pi/xl,0,self.n,self.n,format='coo')
                fbar = np.hstack((-fdata.rd+fdata.ru/ux-fdata.rl/xl,fdata.rp1,fdata.rp2))
                Hbar = coo_matrix((np.concatenate((problem.Hphi.data/self.obj_sca,
                                                   problem.H_combined.data,
                                                   D1.data,
                                                   D2.data)),
                                   (np.concatenate((problem.Hphi.row,
                                                    problem.H_combined.row,
                                                    D1.row,
                                                    D2.row)),
                                    np.concatenate((problem.Hphi.col,
                                                    problem.H_combined.col,
                                                    D1.col,
                                                    D2.col)))))
                Jbar = bmat([[Hbar,None,None],
                             [-self.A,self.Omm1,None],
                             [-problem.J,None,self.Omm2]],
                            format='coo')
                try:
                    if not self.linsolver.is_analyzed():
                        self.linsolver.analyze(Jbar)
                    pbar = self.linsolver.factorize_and_solve(Jbar,fbar)
                except RuntimeError:
                    raise OptSolverError_BadLinSystem(self)
                px = pbar[:self.n]
                pmu = (-fdata.ru + self.mu*px)/ux
                ppi = (-fdata.rl - self.pi*px)/xl
                p = np.hstack((pbar,pmu,ppi))

                # Steplength bounds
                indices = px > 0
                s1 = np.min(np.hstack(((self.u-self.x)[indices]/px[indices],np.inf)))
                indices = px < 0
                s2 = np.min(np.hstack(((self.l-self.x)[indices]/px[indices],np.inf)))
                indices = pmu < 0
                s3 = np.min(np.hstack((-self.mu[indices]/pmu[indices],np.inf)))
                indices = ppi < 0
                s4 = np.min(np.hstack((-self.pi[indices]/ppi[indices],np.inf)))
                smax = (1.-eps)*np.min([s1,s2,s3,s4])
                spmax = (1.-eps)*np.min([s1,s2])
                sdmax = (1.-eps)*np.min([s3,s4])
                
                # Line search
                try:
                    s, fdata = self.line_search(self.y, p, fdata.F, fdata.GradF, self.func, smax=smax, maxiter=ls_maxiter)

                    # Update point
                    self.y += s*p
                    self.x, self.lam, self.nu, self.mu, self.pi = self.extract_components(self.y)
                    
                except OptSolverError_LineSearch:
                    sp = np.minimum(1., spmax)
                    sd = np.minimum(1., sdmax)
                    s = np.minimum(sp,sd)

                    # Update point
                    self.x += sp*px
                    self.lam += sd*pbar[self.x.size:self.x.size+self.lam.size]
                    self.nu += sd*pbar[self.x.size+self.lam.size:]
                    self.mu += sd*pmu
                    self.pi += sd*ppi
                    self.y = np.hstack((self.x,self.lam,self.nu,self.mu,self.pi))

                # Update iters
                self.k += 1                

                # Check
                try:
                    assert(np.all(self.x < self.u))
                    assert(np.all(self.x > self.l))
                    assert(np.all(self.mu > 0))
                    assert(np.all(self.pi > 0))
                except AssertionError:
                    raise OptSolverError_Infeasibility(self)

            # Update iters
            self.k += 1
41626
private void computeMatrixTopoGraphHalfEdges_(int geometry_a, int geometry_b) {
		boolean bRelationKnown = false;

		int id_a = m_topo_graph.getGeometryID(geometry_a);
		int id_b = m_topo_graph.getGeometryID(geometry_b);

		m_visited_index = m_topo_graph.createUserIndexForHalfEdges();

		for (int cluster = m_topo_graph.getFirstCluster(); cluster != -1; cluster = m_topo_graph
				.getNextCluster(cluster)) {
			int first_half_edge = m_topo_graph.getClusterHalfEdge(cluster);
            if (first_half_edge == -1)
            {
                if (m_predicates_cluster != -1)
                {
                    // Treat cluster as an interior point
                    switch (m_predicates_cluster)
                    {
                        case Predicates.AreaPointPredicates:
                            bRelationKnown = areaPointPredicates_(cluster, id_a, id_b);
                            break;
                        case Predicates.LinePointPredicates:
                            bRelationKnown = linePointPredicates_(cluster, id_a, id_b);
                            break;
                        default:
                            throw GeometryException.GeometryInternalError();
                    }
                }

                continue;
            }

			int next_half_edge = first_half_edge;

			do {
				int half_edge = next_half_edge;
				int visited = m_topo_graph.getHalfEdgeUserIndex(half_edge,
						m_visited_index);

				if (visited != 1) {
					do {
						// Invoke relational predicates
						switch (m_predicates_half_edge) {
						case Predicates.AreaAreaPredicates:
							bRelationKnown = areaAreaPredicates_(half_edge,
									id_a, id_b);
							break;
						case Predicates.AreaLinePredicates:
							bRelationKnown = areaLinePredicates_(half_edge,
									id_a, id_b);
							break;
						case Predicates.LineLinePredicates:
							bRelationKnown = lineLinePredicates_(half_edge,
									id_a, id_b);
							break;
						default:
							throw GeometryException.GeometryInternalError();
						}

						if (bRelationKnown)
							break;

						m_topo_graph.setHalfEdgeUserIndex(half_edge,
								m_visited_index, 1);
						half_edge = m_topo_graph.getHalfEdgeNext(half_edge);
					} while (half_edge != next_half_edge && !bRelationKnown);
				}

				if (bRelationKnown)
					break;

				next_half_edge = m_topo_graph.getHalfEdgeNext(m_topo_graph
						.getHalfEdgeTwin(half_edge));
			} while (next_half_edge != first_half_edge);

			if (bRelationKnown)
				break;
		}

		if (!bRelationKnown)
			setRemainingPredicatesToFalse_();

		m_topo_graph.deleteUserIndexForHalfEdges(m_visited_index);
	}
1352829
private HashSet<Integer> vector2set(int y[]) {
        HashSet<Integer> Y = new HashSet<Integer>();
        for (int j = 0; j < y.length; j++) {
            if (y[j] > 0) {
                Y.add(j);
            }
        }
        return Y;
    }
775333
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        int BnumCols = B.numCols;
        
        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                a.data[i*2]   = B.data[indexB];
                a.data[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                double []u = QR[n];

                double realVV = u[n*2];
                double imagVV = u[n*2+1];

                u[n*2]   = 1;
                u[n*2+1] = 0;

                QrHelperFunctions_ZDRM.rank1UpdateMultR(a, u,0, gammas[n], 0, n, numRows, temp.data);

                u[n*2]   = realVV;
                u[n*2+1] = imagVV;
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(R.data, a.data, numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                X.data[indexB]   = a.data[i*2];
                X.data[indexB+1] = a.data[i*2+1];
            }
        }
    }
93854
public DoubleMatrix2D solve(DoubleMatrix2D B) {
		if (B.rows() != dim) {
			log.error("wrong dimension of vector b: expected " + dim + ", actual "	+ B.rows());
			throw new RuntimeException("wrong dimension of vector b: expected " + dim	+ ", actual " + B.rows());
		}

	// with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z
		if (this.rescaler != null) {
			// B = ALG.mult(this.U, B);
			B = ColtUtils.diagonalMatrixMult(this.U, B);
		}

		int nOfColumns = B.columns();

		// Solve LY = B (the same as L.Yc = Bc for each column Yc e Bc)
		final double[][] Y = new double[dim][nOfColumns];
		for (int i = 0; i < diagonalLength; i++) {
			double LII = LData[0][i];
			double[] YI = Y[i];
			DoubleMatrix1D BI = B.viewRow(i);
			for (int col = 0; col < nOfColumns; col++) {
				YI[col] = BI.getQuick(col) / LII;
			}
		}
		for (int i = diagonalLength; i < dim; i++) {
			double[] LI = LData[i - diagonalLength + 1];
			double LII = LI[i];
			double[] sum = new double[nOfColumns];
			for (int j = 0; j < i; j++) {
				double LIJ = LI[j];
				double[] YJ = Y[j];
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LIJ * YJ[col];
				}
			}
			double[] YI = Y[i];
			DoubleMatrix1D BI = B.viewRow(i);
			for (int col = 0; col < nOfColumns; col++) {
				YI[col] = (BI.getQuick(col) - sum[col]) / LII;
			}
		}

		// Solve L[T].X = Y (the same as L[T].Xc = Yc for each column)
		final DoubleMatrix2D X = F2.make(dim, nOfColumns);
		for (int i = dim - 1; i > diagonalLength - 1; i--) {
			double LII = LData[i - diagonalLength + 1][i];
			double[] sum = new double[nOfColumns];
			for (int j = dim - 1; j > i; j--) {
				double[] LJ = LData[j - diagonalLength + 1];
				double LJI = LJ[i];
				DoubleMatrix1D XJ = X.viewRow(j);
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LJI * XJ.getQuick(col);
				}
			}
			DoubleMatrix1D XI = X.viewRow(i);
			double[] YI = Y[i];
			for (int col = 0; col < nOfColumns; col++) {
				XI.setQuick(col, (YI[col] - sum[col]) / LII);
			}
		}
		for (int i = diagonalLength - 1; i > -1; i--) {
			double LII = LData[0][i];
			double[] sum = new double[nOfColumns];
			for (int j = dim - 1; j > diagonalLength - 1; j--) {
				double[] LJ = LData[j - diagonalLength + 1];
				double LJI = LJ[i];
				DoubleMatrix1D XJ = X.viewRow(j);
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LJI * XJ.getQuick(col);
				}
			}
			DoubleMatrix1D XI = X.viewRow(i);
			double[] YI = Y[i];
			for (int col = 0; col < nOfColumns; col++) {
				XI.setQuick(col, (YI[col] - sum[col]) / LII);
			}
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, X);
			return ColtUtils.diagonalMatrixMult(this.U, X);
		} else {
			return X;
		}
	}
856584
def main(dimension, iterations):
    
    objective_function = minimize(functions.sphere)
    stopping_condition = max_iterations(iterations)
    (solution, metrics) = optimize(objective_function=objective_function,
                                   domain=Domain(-5.12, 5.12, dimension),
                                   stopping_condition=stopping_condition,
                                   parameters={'seed': 3758117674,
                                               'rho': 1.0, 'e_s': 15, 'e_f': 5},
                                   velocity_update=gc_velocity_update,
                                   parameter_update=update_rho,
                                   measurements=[fitness_measurement])
    return solution
139995
public static void recomputeCovarianceMatrices(Relation<? extends NumberVector> relation, WritableDataStore<double[]> probClusterIGivenX, List<? extends EMClusterModel<?>> models, double prior) {
    final int k = models.size();
    boolean needsTwoPass = false;
    for(EMClusterModel<?> m : models) {
      m.beginEStep();
      needsTwoPass |= m.needsTwoPass();
    }
    // First pass, only for two-pass models.
    if(needsTwoPass) {
      for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
        double[] clusterProbabilities = probClusterIGivenX.get(iditer);
        NumberVector instance = relation.get(iditer);
        for(int i = 0; i < clusterProbabilities.length; i++) {
          final double prob = clusterProbabilities[i];
          if(prob > 1e-10) {
            models.get(i).firstPassE(instance, prob);
          }
        }
      }
      for(EMClusterModel<?> m : models) {
        m.finalizeFirstPassE();
      }
    }
    double[] wsum = new double[k];
    for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
      double[] clusterProbabilities = probClusterIGivenX.get(iditer);
      NumberVector instance = relation.get(iditer);
      for(int i = 0; i < clusterProbabilities.length; i++) {
        final double prob = clusterProbabilities[i];
        if(prob > 1e-10) {
          models.get(i).updateE(instance, prob);
        }
        wsum[i] += prob;
      }
    }
    for(int i = 0; i < models.size(); i++) {
      // MLE / MAP
      final double weight = prior <= 0. ? wsum[i] / relation.size() : (wsum[i] + prior - 1) / (relation.size() + prior * k - k);
      models.get(i).finalizeEStep(weight, prior);
    }
  }
938951
public void process(Clustering<?> result1, Clustering<?> result2) {
    // Get the clusters
    final List<? extends Cluster<?>> cs1 = result1.getAllClusters();
    final List<? extends Cluster<?>> cs2 = result2.getAllClusters();

    // Initialize
    size1 = cs1.size();
    size2 = cs2.size();
    contingency = new int[size1 + 2][size2 + 2];
    noise1 = BitsUtil.zero(size1);
    noise2 = BitsUtil.zero(size2);

    // Fill main part of matrix
    {
      final Iterator<? extends Cluster<?>> it2 = cs2.iterator();
      for(int i2 = 0; it2.hasNext(); i2++) {
        final Cluster<?> c2 = it2.next();
        if(c2.isNoise()) {
          BitsUtil.setI(noise2, i2);
        }
        contingency[size1 + 1][i2] = c2.size();
        contingency[size1 + 1][size2] += c2.size();
      }
    }
    final Iterator<? extends Cluster<?>> it1 = cs1.iterator();
    for(int i1 = 0; it1.hasNext(); i1++) {
      final Cluster<?> c1 = it1.next();
      if(c1.isNoise()) {
        BitsUtil.setI(noise1, i1);
      }
      final DBIDs ids = DBIDUtil.ensureSet(c1.getIDs());
      contingency[i1][size2 + 1] = c1.size();
      contingency[size1][size2 + 1] += c1.size();

      final Iterator<? extends Cluster<?>> it2 = cs2.iterator();
      for(int i2 = 0; it2.hasNext(); i2++) {
        final Cluster<?> c2 = it2.next();
        int count = DBIDUtil.intersectionSize(ids, c2.getIDs());
        contingency[i1][i2] = count;
        contingency[i1][size2] += count;
        contingency[size1][i2] += count;
        contingency[size1][size2] += count;
      }
    }
  }
21709
private void flood() throws Exception {

        /* define directions */
        // Initialise the vector to a supposed dimension, if the number of
        // unresolved pixel overload the vector there are a method which resized
        // the vectors.
        pitsStackSize = (int) (nCols * nRows * 0.1);
        pstack = pitsStackSize;
        dn = new int[pitsStackSize];
        currentPitRows = new int[pitsStackSize];
        currentPitCols = new int[pitsStackSize];
        ipool = new int[pstack];
        jpool = new int[pstack];
        firstCol = 0;
        firstRow = 0;
        lastCol = nCols;
        lastRow = nRows;

        setdf();

    }
805918
public static INDArray col2im(INDArray col, int sy, int sx, int ph, int pw, int h, int w) {
        //number of images
        long n = col.size(0);
        //number of columns
        long c = col.size(1);
        //kernel height
        long kh = col.size(2);
        //kernel width
        long kw = col.size(3);
        //out height
        long outH = col.size(4);
        //out width
        long outW = col.size(5);

        INDArray img = Nd4j.create(n, c, h + 2 * ph + sy - 1, w + 2 * pw + sx - 1);
        for (int i = 0; i < kh; i++) {
            //iterate over the kernel rows
            long iLim = i + sy * outH;
            for (int j = 0; j < kw; j++) {
                //iterate over the kernel columns
                long jLim = j + sx * outW;
                INDArrayIndex[] indices = new INDArrayIndex[] {NDArrayIndex.all(), NDArrayIndex.all(),
                                NDArrayIndex.interval(i, sy, iLim), NDArrayIndex.interval(j, sx, jLim)};

                INDArray get = img.get(indices);

                INDArray colAdd = col.get(NDArrayIndex.all(), NDArrayIndex.all(), NDArrayIndex.point(i),
                                NDArrayIndex.point(j), NDArrayIndex.all(), NDArrayIndex.all());
                get.addi(colAdd);
                img.put(indices, get);

            }
        }

        //return the subset of the padded image relative to the height/width of the image and the padding width/height
        return img.get(NDArrayIndex.all(), NDArrayIndex.all(), NDArrayIndex.interval(ph, ph + h),
                        NDArrayIndex.interval(pw, pw + w));
    }
644168
public DoubleVector centerOfMass() {
        // Handle lazy initialization
        if (centroid == null) {
            if (indices.size() == 1)
                centroid = sumVector;
            else {
                // Update the centroid by normalizing by the number of elements.
                // We expect that the centroid vector might be compared with
                // other vectors multiple times. If we used a ScaledVector here,
                // we would be re-doing this multiplication each time, which is
                // wasted. The centerOfMass is already lazily instantiated, so
                // we know that if we do the computation here we'll be using the
                // results at least once.  Therefore do the normalization here
                // once to save cost.
                int length = sumVector.length();
                double d = 1d / indices.size();
                if (sumVector instanceof SparseVector) {
                    centroid = new SparseHashDoubleVector(length);
                    SparseVector sv = (SparseVector)sumVector;
                    for (int nz : sv.getNonZeroIndices())
                        centroid.set(nz, sumVector.get(nz) * d);
                }
                else {
                    centroid = new DenseVector(length);
                    for (int i = 0; i < length; ++i) 
                        centroid.set(i, sumVector.get(i) * d);                    
                }
            }
        }
        return centroid;
    }
673899
private static void clusterIteration(Matrix matrix,
                                         int numClusters,
                                         KMeansSeed seedType,
                                         CriterionFunction criterion) {
        DoubleVector[] centers = seedType.chooseSeeds(numClusters, matrix);

        // Compute the initial set of assignments for each data point based on
        // the initial assignments.
        int[] initialAssignments = new int[matrix.rows()];

        // If there is to be only one cluster, then everything will be auto
        // assigned to the first cluster.  This is just a special case that only
        // comes up when comparing other solutions to the non-solution.
        if (numClusters != 1) {
            int nc = 0;
            for (int i = 0; i < matrix.rows(); ++i) {

                DoubleVector vector = matrix.getRowVector(i);
                double bestSimilarity = 0;
                for (int c = 0; c < numClusters; ++c) {
                    double similarity = Similarity.cosineSimilarity(
                            centers[c], vector);
                    nc++;
                    if (similarity >= bestSimilarity) {
                        bestSimilarity = similarity;
                        initialAssignments[i] = c;
                    }
                }
            }
        }

        // Setup the criterion function with it's meta data.
        criterion.setup(matrix, initialAssignments, numClusters);

        // Iteratively swap each data point to a better cluster if such an
        // assignment exists.
        List<Integer> indices = new ArrayList<Integer>(matrix.rows());
        for (int i = 0; i < matrix.rows(); ++i)
            indices.add(i);

        // Iterate through each data point in random order.  For each data
        // point, try to assign it to a new cluster.  If no data point is moved
        // in an iteration, end the iterations.
        boolean changed = true;
        while (changed) {
            changed = false;
            Collections.shuffle(indices);
            for (int index : indices)
                changed |= criterion.update(index);
        }
    }
208523
private double optimize(int k0, int kmax, double sumx, double sumx2, double xp, double yp, double sumxy, double sumy) {
    int k = kmax - k0 + 1;
    return (sumxy - xp * sumy - yp * sumx + k * xp * yp) / (sumx2 - 2 * sumx * xp + k * xp * xp);
    // return (-xp * yp * + yp * sumx - sumxy + xp * sumy) / (-xp * xp *
    // kmax - sumx2 + 2 * xp * sumx);
  }
847270
function vectorize(fun) {
  return function(x) {
    if (!isArray(x)) {
      return fun(x);
    }
    return x.map(function(x) {
      return fun(x);
    });
  };
}
500024
public DoubleMatrix1D solve(DoubleMatrix1D b) {
		if (b.size() != dim) {
			log.error("wrong dimension of vector b: expected " + dim + ", actual "	+ b.size());
			throw new RuntimeException("wrong dimension of vector b: expected " + dim	+ ", actual " + b.size());
		}

		// with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z
		if (this.rescaler != null) {
			// b = ALG.mult(this.U, b);
			b = ColtUtils.diagonalMatrixMult(this.U, b);
		}

		// Solve L.y = b
		final double[] y = new double[dim];
		for (int i = 0; i < dim; i++) {
			double[] LI = LData[i];
			double sum = 0;
			for (int j = 0; j < i; j++) {
				sum += LI[j] * y[j];
			}
			y[i] = (b.getQuick(i) - sum) / LI[i];
		}

		// logger.debug("b: " + ArrayUtils.toString(b));
		// logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y)));

		// Solve D.z = y
		final double[] z = new double[dim];
		for (int i = 0; i < dim; i++) {
			z[i] = y[i] / DData[i];
		}

		// Solve L[T].x = z
		final DoubleMatrix1D x = F1.make(dim);
		for (int i = dim - 1; i > -1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > i; j--) {
				sum += LData[j][i] * x.getQuick(j);
			}
			x.setQuick(i, (z[i] - sum) / LData[i][i]);
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, x);
			return ColtUtils.diagonalMatrixMult(this.U, x);
		} else {
			return x;
		}
	}
1166769
private void precomputeSlope() {
        float invDirX = 1.0f / dirX;
        float invDirY = 1.0f / dirY;
        float invDirZ = 1.0f / dirZ;
        s_yx = dirX * invDirY;
        s_xy = dirY * invDirX;
        s_zy = dirY * invDirZ;
        s_yz = dirZ * invDirY;
        s_xz = dirZ * invDirX;
        s_zx = dirX * invDirZ;
        c_xy = originY - s_xy * originX;
        c_yx = originX - s_yx * originY;
        c_zy = originY - s_zy * originZ;
        c_yz = originZ - s_yz * originY;
        c_xz = originZ - s_xz * originX; // <- original paper had a bug here. It switched originZ/originX
        c_zx = originX - s_zx * originZ; // <- original paper had a bug here. It switched originZ/originX
        int sgnX = signum(dirX);
        int sgnY = signum(dirY);
        int sgnZ = signum(dirZ);
        classification = (byte) ((sgnZ+1) << 4 | (sgnY+1) << 2 | (sgnX+1));
    }