Java Code Examples for org.ejml.simple.SimpleMatrix#set()

/**
	 * 用textRank算法计算权重矩阵。
	 * 
	 * @param matrix
	 * @return
	 */
	protected SimpleMatrix buildWeightVector(SimpleMatrix matrix) {
		SimpleMatrix vector = new SimpleMatrix(matrix.numCols(), 1);
		vector.set(1);
		
		SimpleMatrix vecDamp = new SimpleMatrix(matrix.numCols(), 1);
		vecDamp.set(1- this.damp);
		
		double diff = 1;
		while(diff > this.threshold){
			SimpleMatrix next = matrix.mult(vector);
			//next = (1-damp)+damp * next;
			next = vecDamp.plus(this.damp, next);			
			diff = next.minus(vector).normF();
			vector = next;
//			System.out.println("weight==========");
//			System.out.println(vector);
		}
		return vector;
	}
 

/**
 * 建立相似度矩阵
 * 
 * @param sentences
 * @return
 */
protected SimpleMatrix buildSimilarityMatrix(List<Sentence> sentences) {
	for (Sentence sentence : sentences) {
		List<Term> terms = this.tokenizer.tokenize(sentence.toString());
		((SentenceWrapper) sentence).setTerms(terms);
	}
	SimpleMatrix matrix = new SimpleMatrix(sentences.size(),
			sentences.size());
	matrix.set(0);
	for (int i = 0; i < sentences.size(); i++)
		for (int j = i + 1; j < sentences.size(); j++) {
			// 相似度+1，消除0值；
			double similarity = this.similarity(sentences.get(i),
					sentences.get(j)) + 1;
			matrix.set(i, j, similarity);
			matrix.set(j, i, similarity);
		}
	return matrix;
}
 

private static SimpleMatrix transformMatrix(SimpleMatrix trnsF, SimpleMatrix matrix, Double[] validElements, int countValidElements) {
	// forward transform the pdf and remove non-valid eigendirections
	//if(matrix.numRows()!= countValidElements)
	//	System.out.println("projection was necessary!!!");
	SimpleMatrix tmp = trnsF.mult(matrix).mult(trnsF.transpose());
	SimpleMatrix trnsMatrix = new SimpleMatrix(countValidElements, countValidElements);
	int row=0, column=0;
	for(int i=0; i< validElements.length; i++){
		for(int j=0; j< validElements.length; j++){
			if(validElements[i] == 1 && validElements[j] == 1)
				trnsMatrix.set(row, column++, tmp.get(i, j));
		}
		column = 0;
		row++;
	}
	return trnsMatrix;
}
 

/**
 * Splits a single component distribution into two components as described in the oKDE-paper.
 * @return a TwoComponentDistribution
 */
public TwoComponentDistribution split(double parentWeight){
	SimpleSVD<?> svd = mGlobalCovariance.svd(true);
	SimpleMatrix S = svd.getW();
	SimpleMatrix V = svd.getV();
	SimpleMatrix d = S.extractDiag();
	double max = MatrixOps.maxVectorElement(d);
	int maxIndex = MatrixOps.maxVectorElementIndex(d);
	int len = mGlobalCovariance.numRows();
	SimpleMatrix M = new SimpleMatrix(len,1);
	M.set(maxIndex, 0, 1.0d);
	SimpleMatrix dMean = V.mult(M).scale(0.5*Math.sqrt(max));
	SimpleMatrix meanSplit1 = mGlobalMean.plus(dMean);
	SimpleMatrix meanSplit2 = mGlobalMean.minus(dMean);
	
	SimpleMatrix dyadMean = mGlobalMean.mult(mGlobalMean.transpose());
	SimpleMatrix dyadMeanSplit1 = meanSplit1.mult(meanSplit1.transpose());
	SimpleMatrix dyadMeanSplit2 = meanSplit2.mult(meanSplit2.transpose());
	SimpleMatrix covSplit = mGlobalCovariance.plus(dyadMean).minus(dyadMeanSplit1.plus(dyadMeanSplit2).scale(0.5));
	
	SimpleMatrix[] means = {meanSplit1, meanSplit2};
	SimpleMatrix[] covariances = {covSplit, covSplit};
	double[] weights = {0.5, 0.5};
	TwoComponentDistribution splitDist = null;
	try {
		splitDist = new TwoComponentDistribution(weights, means, covariances, mBandwidthMatrix);
		splitDist.setGlobalWeight(parentWeight*mGlobalWeight);
		splitDist.setGlobalCovariance(mGlobalCovariance);
		splitDist.setGlobalMean(mGlobalMean);
	} catch (TooManyComponentsException e) {
		// cant be thrown
	}
	return splitDist;
}
 

/**
 * Factory method for converting RxPos to a SimpleMatrix
 * @param rxPos rxPos Coordinates object
 * @return RxPos as a 4z1 vector
 */
public static SimpleMatrix getRxPosAsVector(Coordinates rxPos){
    SimpleMatrix rxPosSimpleVector = new SimpleMatrix(4, 1);
    rxPosSimpleVector.set(0, rxPos.getX());
    rxPosSimpleVector.set(1, rxPos.getY());
    rxPosSimpleVector.set(2, rxPos.getZ());
    rxPosSimpleVector.set(3, 0);

    return rxPosSimpleVector;
}
 

private static SimpleMatrix projectBandwidthToOriginalSpace(SampleModel distribution, double bandwidthFactor) {
	SimpleMatrix bandwidth = SimpleMatrix.identity(distribution.getGlobalCovariance().numCols());
	SimpleMatrix subSpaceBandwidth = distribution.getSubspaceGlobalCovariance().scale(Math.pow(bandwidthFactor, 2));
	ArrayList<Integer> subspace = distribution.getmSubspace();
	for (int i = 0; i < subSpaceBandwidth.numRows(); i++) {
		for (int j = 0; j < subSpaceBandwidth.numCols(); j++) {
			if (subspace.contains(new Integer(i)) && subspace.contains(new Integer(j)))
				bandwidth.set(i, j, subSpaceBandwidth.get(i, j));
		}
	}
	SimpleMatrix invSubspaceCov = distribution.getSubspaceInverseCovariance();
	bandwidth = invSubspaceCov.transpose().mult(bandwidth).mult(invSubspaceCov);
	return bandwidth;
}
 

public static SimpleMatrix ones(int rows, int cols) {
	SimpleMatrix matrix = new SimpleMatrix(rows, cols);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < rows; j++) {
			matrix.set(i, j, 1);
		}
	}
	return matrix;
}
 

public static SimpleMatrix deleteElementsFromVector(SimpleMatrix vector, List<Double> elements, int vectorSize) {
	SimpleMatrix newVector = new SimpleMatrix(vectorSize, 1);
	int j = 0;
	for (int i = 0; i < vector.numRows(); i++)
		if (elements.get(i) == 1)
			newVector.set(j++, 0, vector.get(i));
	return newVector;
}
 

public static SimpleMatrix elemPow(SimpleMatrix matrix, double p) {
	for (int i = 0; i < matrix.numRows(); i++) {
		for (int j = 0; j < matrix.numCols(); j++) {
			matrix.set(i, j, Math.pow(matrix.get(i, j), p));
		}
	}
	return matrix;
}
 

public static SimpleMatrix elemSqrt(SimpleMatrix matrix) {
	for (int i = 0; i < matrix.numRows(); i++) {
		for (int j = 0; j < matrix.numCols(); j++) {
			matrix.set(i, j, Math.sqrt(matrix.get(i, j)));
		}
	}
	return matrix;
}
 

public static SimpleMatrix abs(SimpleMatrix matrix) {
	for (int i = 0; i < matrix.numRows(); i++) {
		for (int j = 0; j < matrix.numCols(); j++) {
			matrix.set(i, j, Math.abs(matrix.get(i, j)));
		}
	}
	return matrix;
}
 

private static SimpleMatrix setVectorElements(SimpleMatrix v1, SimpleMatrix v2, Double[] elementsInV1) {
	int j = 0;
	for (int i = 0; i < v1.numRows(); i++)
		if (elementsInV1[i] == 1)
			v1.set(i, 0, v2.get(j++));
	return v1;
}
 

private static SimpleMatrix backTransformMatrix(SimpleMatrix matrix, SimpleMatrix matrixToTransform, Double[] validElements) {
	// add removed eigendirections and backwards transform
	int row=0, column=0;
	for(int i=0; i< validElements.length; i++){
		for(int j=0; j< validElements.length; j++){
			if(validElements[i] == 1 && validElements[j] == 1)
				matrix.set(i, j, matrixToTransform.get(row, column++));
		}
		column = 0;
		row++;
	}
	
	/*
	
	int a=0,b=0;
	for (int row = 0; row < matrix.numRows(); row++) {
		for (int column = 0; column < matrix.numCols(); column++) {
			if (validElements[row] == 1)
				matrix.set(row, column, matrixToTransform.get(row, column));
			else if (validElements[column] == 1)
				matrix.set(row, column, matrixToTransform.get(row, column));
			else
				matrix.set(row, column, 0);
		}
	}*/
	return matrix;
}
 

/**
 * The team assignment matrix is often referred to as the "A" matrix.
 * It's a matrix whose rows represent the players and the columns
 * represent teams. At Matrix[row, column] represents that player[row]
 * is on team[col] Positive values represent an assignment and a
 * negative value means that we subtract the value of the next team
 * since we're dealing with pairs. This means that this matrix always
 * has teams - 1 columns. The only other tricky thing is that values
 * represent the play percentage.
 * <p>
 * For example, consider a 3 team game where team1 is just player1, team
 * 2 is player 2 and player 3, and team3 is just player 4. Furthermore,
 * player 2 and player 3 on team 2 played 25% and 75% of the time (e.g.
 * partial play), the A matrix would be:
 * <p><pre>
 * A = this 4x2 matrix:
 * |  1.00  0.00 |
 * | -0.25  0.25 |
 * | -0.75  0.75 |
 * |  0.00 -1.00 |
 * </pre>
 */
private static SimpleMatrix CreatePlayerTeamAssignmentMatrix(List<ITeam> teamAssignmentsList, int totalPlayers) {

    List<List<Double>> playerAssignments = new ArrayList<List<Double>>();
    int totalPreviousPlayers = 0;

    for (int i = 0; i < teamAssignmentsList.size() - 1; i++) {
        ITeam currentTeam = teamAssignmentsList.get(i);

        // Need to add in 0's for all the previous players, since they're not
        // on this team
        List<Double> currentRowValues = new ArrayList<Double>();
        for(int j = 0; j < totalPreviousPlayers; j++) currentRowValues.add(0.);
        playerAssignments.add(currentRowValues);

        for(IPlayer player: currentTeam.keySet()) {
            currentRowValues.add(PartialPlay.getPartialPlayPercentage(player));
            // indicates the player is on the team
            totalPreviousPlayers++;
        }

        ITeam nextTeam = teamAssignmentsList.get(i + 1);
        for(IPlayer nextTeamPlayer : nextTeam.keySet()) {
            // Add a -1 * playing time to represent the difference
            currentRowValues.add(-1 * PartialPlay.getPartialPlayPercentage(nextTeamPlayer));
        }
    }

    SimpleMatrix playerTeamAssignmentsMatrix = new SimpleMatrix(totalPlayers, teamAssignmentsList.size() - 1);
    for(int i=0; i < playerAssignments.size(); i++)
        for(int j=0; j < playerAssignments.get(i).size(); j++)
            playerTeamAssignmentsMatrix.set(j, i, playerAssignments.get(i).get(j));

    return playerTeamAssignmentsMatrix;
}
 

/**
 * 将matrix归一化处理。
 * 
 * @param matrix
 */
protected void normalize(SimpleMatrix matrix) {
	SimpleMatrix one = new SimpleMatrix(matrix.numCols(), matrix.numRows());
	one.set(1);
	SimpleMatrix sum = matrix.mult(one);
	//CommonOps.elementDiv(matrix.getMatrix(), sum.getMatrix());
	matrix.set(matrix.elementDiv(sum));
	//CommonOps.transpose(matrix.getMatrix());
	matrix.set(matrix.transpose());
}
 

/**
 * This method searches a local maximum by gradient-quadratic search. First a direct leap to the maximum by 
 * quadratic optimization is tried. Then gradient search is used to refine the result in case of an overshoot.
 * Uses means, covariances and component weights given as parameters.
 * 
 * This algorithm was motivated by this paper: 
 * Miguel A. Carreira-Perpinan (2000): "Mode-finding for mixtures of
 * Gaussian distributions", IEEE Trans. on Pattern Analysis and
 * Machine Intelligence 22(11): 1318-1323.
 * 
 * @param start Defines the starting point for the search.
 * @return The serach result containing the point and the probability value at that point.
 */
public static SearchResult gradQuadrSearch(SimpleMatrix start, ArrayList<SimpleMatrix> means, ArrayList<SimpleMatrix> covs, ArrayList<Double> weights, SampleModel model){

	SimpleMatrix gradient = new SimpleMatrix(2,1);
	SimpleMatrix hessian = new SimpleMatrix(2,2);
	double n = means.get(0).numRows();
	double a = Math.pow(Math.sqrt(2 * Math.PI), n);
	
	SimpleMatrix x = new SimpleMatrix(2,1);
	x.set(0,0,start.get(start.numRows()-2,0));
	x.set(1,0,start.get(start.numRows()-1,0));
	ArrayList<Double> mahalanobisDistances;
	double step = START_STEP_SIZE;
	double probability = 0;
	SimpleMatrix gradStep = null;
	int count =0;
	do {
		mahalanobisDistances = mahalanobis(x, means, covs);
		double prob = 0;
		// this loop calculates gradient and hessian as well as probability at x
		for (int i = 0; i < means.size(); i++) {
			// check whether the component actually contributes to to the density at given point by mahalanobis distance
			if(mahalanobisDistances.get(i) < MAX_MAHALANOBIS_DIST) {
				SimpleMatrix m = means.get(i);
				SimpleMatrix dm = m.minus(x);
				SimpleMatrix c = covs.get(i);
				SimpleMatrix invC = c.invert();
				double w = weights.get(i);
				//probability p(x,m) under component m
				double p = ((1 / (a * Math.sqrt(c.determinant()))) * Math.exp((-0.5d) * mahalanobisDistances.get(i))) * w;
				prob += p; 
				// gradient at x
				gradient = gradient.plus( invC.mult(dm).scale(p) );
				// hessian at x
				hessian = hessian.plus( invC.mult( dm.mult(dm.transpose()).minus(c) ).mult(invC).scale(p) );
			}


		}
		// save x
		SimpleMatrix xOld = new SimpleMatrix(x);
		double tst = evaluate(xOld, means, covs, weights);
		// check if hessian is negative definite
		SimpleEVD hessianEVD = hessian.eig();
		int maxEVIndex = hessianEVD.getIndexMax();
		// try a direct leap by quadratic optimization
		if(hessianEVD.getEigenvalue(maxEVIndex).getReal() < 0){
			gradStep = hessian.invert().mult(gradient);
			x = xOld.minus(gradStep);
		}
		double prob1 = 	evaluate(x, means, covs, weights);
		// if quadratic optimization did not work try gradient ascent
		if( prob1 <= prob | hessianEVD.getEigenvalue(maxEVIndex).getReal() >= 0) {
			gradStep = gradient.scale(step);
			x = xOld.plus(gradStep);
			// if still not ok decrease step size iteratively
			while(evaluate(x, means, covs, weights) < prob){
				step = step/2;
				gradStep = gradient.scale(step);
				x = xOld.plus(gradStep);
			}
		}
		probability =	model.evaluate(x, means, covs, weights);
		count++;
		// continue until the last step is sufficiently small or
		// a predefined amount of steps was performed
	}while(gradStep.elementMaxAbs() > STOP_STEP_SIZE && count<10);

	// return results
	return new SearchResult(x, probability);
}
 

public static SimpleMatrix doubleListToMatrix(List<Double> valueList) {
	SimpleMatrix m = new SimpleMatrix(1, valueList.size());
	for (int i = 0; i < valueList.size(); i++)
		m.set(0, i, valueList.get(i));
	return m;
}
 

/**
 * This method derives the conditional distribution of the sample model kde with distribution p(x).
 * It marginalizes on the first n dimensions. The number n is given as a parameter.
 * @param firstDimensions
 * @return
 */
public ConditionalDistribution getMarginalDistribution(int n){
	ArrayList<SimpleMatrix> means = this.getSubMeans();
	ArrayList<SimpleMatrix> marginalMeans = new ArrayList<SimpleMatrix>();

	ArrayList<SimpleMatrix> covs = this.getSubSmoothedCovariances();
	ArrayList<SimpleMatrix> marginalCovs = new ArrayList<SimpleMatrix>();
	
	ArrayList<Double> weights = this.getSubWeights();
	ArrayList<Double> marginalWeights = new ArrayList<Double>();

	ConditionalDistribution result = null;
	
	double a = Math.pow(Math.sqrt(2 * Math.PI), n);

	
	for(int i=0; i<means.size(); i++) {
		SimpleMatrix c = covs.get(i);
		SimpleMatrix m = means.get(i);
		SimpleMatrix m1 = new SimpleMatrix(n,1);
		
		// extract all elements from covariance that correspond only to m1
		// that means extract the block in the left top corner with height=width=n
		SimpleMatrix newC1 = new SimpleMatrix(n,n);
		for(int j=0; j<n; j++) {
			for(int k=0; k<n; k++) {
				newC1.set(j, k, c.get(j,k) );
			}
		}
		//extract last rows from mean to m1
		for(int j=0; j<n; j++) {
			m1.set(j,0,m.get(j,0));
		}
		
		marginalMeans.add(m1);
		marginalCovs.add(newC1);

	}
	result = new ConditionalDistribution(marginalMeans, marginalCovs, weights);
	return result;
}
 

/**
 * This method derives the conditional distribution of the actual sample model kde with distribution p(x).
 * It takes a condition parameter that is a vector c of dimension m. Using this vector
 * it finds the conditional distribution p(x*|c) where c=(x_0,...,x_m), x*=(x_m+1,...,x_n).
 * For detailed description see:
 * @param condition A vector that defines c in p(x*|c)
 * @return The conditional distribution of this sample model under the given condition
 */
public ConditionalDistribution getConditionalDistribution(SimpleMatrix condition){
	int lenCond = condition.numRows();
	
	ArrayList<SimpleMatrix> means = this.getSubMeans();
	ArrayList<SimpleMatrix> conditionalMeans = new ArrayList<SimpleMatrix>();

	ArrayList<SimpleMatrix> covs = this.getSubSmoothedCovariances();
	ArrayList<SimpleMatrix> conditionalCovs = new ArrayList<SimpleMatrix>();
	
	ArrayList<Double> weights = this.getSubWeights();
	ArrayList<Double> conditionalWeights = new ArrayList<Double>();

	ConditionalDistribution result = null;
	
	double n = condition.numRows();
	double a = Math.pow(Math.sqrt(2 * Math.PI), n);

	
	for(int i=0; i<means.size(); i++) {
		SimpleMatrix c = covs.get(i);
		SimpleMatrix invC = c.invert();
		SimpleMatrix m = means.get(i);
		int lenM1 = m.numRows()-lenCond;
		SimpleMatrix m1 = new SimpleMatrix(lenM1,1);
		SimpleMatrix m2 = new SimpleMatrix(lenCond,1);
		
		// extract all elements from inverse covariance that correspond only to m1
		// that means extract the block in the right bottom corner with height=width=lenM1
		SimpleMatrix newC1 = new SimpleMatrix(lenM1,lenM1);
		for(int j=0; j<lenM1; j++) {
			for(int k=0; k<lenM1; k++) {
				newC1.set(j, k, invC.get(j+lenCond,k+lenCond) );
			}
		}
		// extract all elements from inverse covariance that correspond to m1 and m2
		// from the the block in the left bottom corner with height=width=lenM1
		SimpleMatrix newC2 = new SimpleMatrix(lenM1,lenCond);
		for(int j=0; j<lenM1; j++) {
			for(int k=0; k<lenCond; k++) {
				newC2.set(j, k, invC.get(j+lenCond,k) );
			}
		}		
		
		//extract first rows from mean to m2
		for(int j=0; j<lenCond; j++) {
			m2.set(j,0,m.get(j,0));
		}
		//extract last rows from mean to m1
		for(int j=0; j<lenM1; j++) {
			m1.set(j,0,m.get(j+lenCond,0));
		}
		SimpleMatrix invNewC1 = newC1.invert();
		// calculate new mean and new covariance of conditional distribution
		SimpleMatrix condMean = m1.minus( invNewC1.mult(newC2).mult( condition.minus(m2) ) );
		SimpleMatrix condCovariance = invNewC1;
		conditionalMeans.add(condMean);
		conditionalCovs.add(condCovariance);
		
		// calculate new weights
		
		// extract all elements from inverse covariance that correspond only to m2
		// that means extract the block in the left top corner with height=width=lenCond
		SimpleMatrix newC22 = new SimpleMatrix(lenCond,lenCond);
		for(int j=0; j<lenCond; j++) {
			for(int k=0; k<lenCond; k++) {
				newC22.set(j, k, c.get(j,k) );
			}
		}
		double mahalanobisDistance = condition.minus(m2).transpose().mult(newC22.invert()).mult(condition.minus(m2)).trace();
		double newWeight = ((1 / (a * Math.sqrt(newC22.determinant()))) * Math.exp((-0.5d) * mahalanobisDistance))* weights.get(i);
		conditionalWeights.add(newWeight);
	}
	// normalize weights
	double weightSum = 0;
	for(int i=0; i<conditionalWeights.size(); i++) {
		weightSum += conditionalWeights.get(i);
	}
	for(int i=0; i<conditionalWeights.size(); i++) {
		double weight = conditionalWeights.get(i);
		weight = weight /weightSum;
		conditionalWeights.set(i,weight);
	}
	result = new ConditionalDistribution(conditionalMeans, conditionalCovs, conditionalWeights);
	return result;
}
 

/**
 * Find Maximum by gradient-quadratic search.
 * First a conditional distribution is derived from the kde.
 * @param start
 * @return
 */
public SearchResult gradQuadrSearch(SimpleMatrix start){
	
	
	SimpleMatrix condVector = new SimpleMatrix(4,1);
	for(int i=0; i<condVector.numRows(); i++){
		condVector.set(i,0,start.get(i,0));
	}
	ConditionalDistribution conditionalDist = getConditionalDistribution(condVector);
	
	ArrayList<SimpleMatrix> means = conditionalDist.conditionalMeans;
	ArrayList<SimpleMatrix> covs = conditionalDist.conditionalCovs;
	ArrayList<Double> weights = conditionalDist.conditionalWeights;

	SimpleMatrix gradient = new SimpleMatrix(2,1);
	SimpleMatrix hessian = new SimpleMatrix(2,2);
	double n = means.get(0).numRows();
	double a = Math.pow(Math.sqrt(2 * Math.PI), n);
	
	SimpleMatrix x = new SimpleMatrix(2,1);
	x.set(0,0,start.get(start.numRows()-2,0));
	x.set(1,0,start.get(start.numRows()-1,0));
	ArrayList<Double> mahalanobisDistances;
	double step = 1;
	double probability = 0;
	SimpleMatrix gradStep = null;
	do {
		mahalanobisDistances = mahalanobis(x, means, covs);
		//calculate gradient and hessian:
		double prob = 0;
		for (int i = 0; i < means.size(); i++) {
			// check wether the component actually contributes to to the density at given point 
			if(mahalanobisDistances.get(i) < MAX_MAHALANOBIS_DIST) {
				SimpleMatrix m = means.get(i);

				SimpleMatrix dm = m.minus(x);
				SimpleMatrix c = covs.get(i);

				
				SimpleMatrix invC = c.invert();
				double w = weights.get(i);
				//probability p(x,m)
				double p = ((1 / (a * Math.sqrt(c.determinant()))) * Math.exp((-0.5d) * mahalanobisDistances.get(i))) * w;
				prob += p; 
				gradient = gradient.plus( invC.mult(dm).scale(p) );
				hessian = hessian.plus( invC.mult( dm.mult(dm.transpose()).minus(c) ).mult(invC).scale(p) );
			}


		}
		// save x
		SimpleMatrix xOld = new SimpleMatrix(x);
		SimpleEVD<?> hessianEVD = hessian.eig();
		int maxEVIndex = hessianEVD.getIndexMax();
		if(hessianEVD.getEigenvalue(maxEVIndex).getReal() < 0){
			gradStep = hessian.invert().mult(gradient);
			x = xOld.minus(gradStep);
		}
		double prob1 = 	evaluate(x, means, covs, weights);
		if( prob1 <= prob | hessianEVD.getEigenvalue(maxEVIndex).getReal() >= 0) {
			gradStep = gradient.scale(step);
			x = xOld.plus(gradStep);
			while(evaluate(x, means, covs, weights) < prob){
				step = step/2;
				gradStep = gradient.scale(step);
				x = xOld.plus(gradStep);
			}
		}
		probability =	evaluate(x, means, covs, weights); 
	}while(gradStep.elementMaxAbs() > 1E-10);
	
	return new SearchResult(x, probability);
}