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

/**
 * Returns 2n+k sigma points starting with mean as the first point
 * 
 * @param mean
 * @param cov
 * @param no
 * @param k
 * @return
 */
private static List<SimpleMatrix> getSigmaPoints(SimpleMatrix mean, SimpleMatrix cov, int no, int k) {
	List<SimpleMatrix> resultVectors = new ArrayList<SimpleMatrix>();

	int n = cov.numRows();
	SimpleSVD<?> svd = cov.svd(true);
	SimpleMatrix U = svd.getU();
	SimpleMatrix S = svd.getW();

	S = U.mult(MatrixOps.elemSqrt(S)).scale(Math.sqrt(n + k));

	for (int i = 0; i < S.numCols(); i++) {
		SimpleMatrix columnVector = S.extractVector(false, i);
		SimpleMatrix negColumnVector = S.extractVector(false, i).scale(-1);
		resultVectors.add(columnVector.plus(mean));
		resultVectors.add(negColumnVector.plus(mean));
	}
	if (k != 0)
		resultVectors.add(mean);

	return resultVectors;
}
 

/**
 * Evaluate a gaussian mixture defined by the given means, covariances and component weights at a given point x.
 * 
 * @param x The point where the mixture shall be evaluated.
 * @param means The component means.
 * @param covs The component covariances.
 * @param weights The component weights.
 * @return The probability at point x.
 */
private static double evaluate(SimpleMatrix x, List<SimpleMatrix> means, List<SimpleMatrix> covs, List<Double> weights){
	ArrayList<Double> mahalanobisDistances = mahalanobis(x, means, covs);
	double n = x.numRows();
	double a = Math.pow(Math.sqrt(2 * Math.PI), n);

	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 c = covs.get(i);
			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; 
		}
	}
	return prob;
}
 

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;
}
 

@Override
public double[] value(SimpleMatrix xx) {
    int n = xx.numRows();
    double [] retval = new double[n];
    for (int i = 0; i < n; i++){
            retval[i] = value(xx.extractVector(true, i).getMatrix().getData());
    }
		
    return retval;    
}
 

@Override
public double[] value(SimpleMatrix xx) {
    int n = xx.numRows();
    double [] retval = new double[n];
    for (int i = 0; i < n; i++){
            retval[i] = value(xx.extractVector(true, i).getMatrix().getData());
    }
		
    return retval;    
}
 

@Override
public double[] value(SimpleMatrix xx) {
    int n = xx.numRows();
    double [] retval = new double[n];
    for (int i = 0; i < n; i++){
            retval[i] = value(xx.extractVector(true, i).getMatrix().getData());
    }
		
    return retval;    
}
 

@Override
public double[] value(SimpleMatrix xx) {
    int n = xx.numRows();
    double [] retval = new double[n];
    for (int i = 0; i < n; i++){
            retval[i] = value(xx.extractVector(true, i).getMatrix().getData());
    }
		
    return retval;    
}
 

@Override
public double[] value(SimpleMatrix xx) {
    int n = xx.numRows();
    double [] retval = new double[n];
    for (int i = 0; i < n; i++){
            retval[i] = value(xx.extractVector(true, i).getMatrix().getData());
    }
		
    return retval;    
}
 

@Override
public double[] value(SimpleMatrix xx){
    
    int n = xx.numRows();
    double [] retval = new double[n];
    for (int i = 0; i < n; i++){
            retval[i] = value(xx.extractVector(true, i).getMatrix().getData());
    }
		
    return retval;

}
 

/**
 * 将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());
}
 

public static int maxVectorElementIndex(SimpleMatrix matrix){
	double d = Double.MIN_VALUE;
	int row = 0;
	for (int i = 0; i < matrix.numRows(); i++) {
			if(matrix.get(i,0)>d){
				d = matrix.get(i,0);
				row = i;
			}
	}
	return row;
}
 

public static double maxVectorElement(SimpleMatrix matrix){
	double d = Double.MIN_VALUE;
	for (int i = 0; i < matrix.numRows(); i++) {
			if(matrix.get(i,0)>d)
				d = matrix.get(i,0);
	}
	return d;
}
 

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 double getIntSquaredHessian(SimpleMatrix[] means, Double[] weights, SimpleMatrix[] covariance, SimpleMatrix F, SimpleMatrix g) {
	long time = System.currentTimeMillis();
	long d = means[0].numRows();
	long N = means.length;
	// normalizer
	double constNorm = Math.pow((1d / (2d * Math.PI)), (d / 2d));

	// test if F is identity for speedup
	SimpleMatrix Id = SimpleMatrix.identity(F.numCols());
	double deltaF = F.minus(Id).elementSum();

	double w1, w2, m, I = 0, eta, f_t, c;
	SimpleMatrix s1, s2, mu1, mu2, dm, ds, B, b, C;
	for (int i1 = 0; i1 < N; i1++) {
		s1 = covariance[i1].plus(g);
		mu1 = means[i1];
		w1 = weights[i1];
		for (int i2 = i1; i2 < N; i2++) {
			s2 = covariance[i2];
			mu2 = means[i2];
			w2 = weights[i2];
			SimpleMatrix A = s1.plus(s2).invert();
			dm = mu1.minus(mu2);

			// if F is not identity
			if (deltaF > 1e-3) {
				ds = dm.transpose().mult(A);
				b = ds.transpose().mult(ds);
				B = A.minus(b.scale(2));
				C = A.minus(b);
				f_t = constNorm * Math.sqrt(A.determinant()) * Math.exp(-0.5 * ds.mult(dm).trace());
				c = 2 * F.mult(A).mult(F).mult(B).trace() + Math.pow(F.mult(C).trace(), 2);
			} else {
				m = dm.transpose().mult(A).mult(dm).get(0);
				f_t = constNorm * Math.sqrt(A.determinant()) * Math.exp(-0.5 * m);

				DenseMatrix64F A_sqr = new DenseMatrix64F(A.numRows(), A.numCols());
				CommonOps.elementMult(A.getMatrix(), A.transpose().getMatrix(), A_sqr);
				double sum = CommonOps.elementSum(A_sqr);
				c = 2d * sum * (1d - 2d * m) + Math.pow((1d - m), 2d) * Math.pow(A.trace(), 2);
			}

			// determine the weight of the current term
			if (i1 == i2)
				eta = 1;
			else
				eta = 2;
			I = I + f_t * c * w2 * w1 * eta;
		}
	}
	/*time = System.currentTimeMillis()-time;
	if((time) > 100)
		System.out.println("Time for IntSqrdHessian: "+ ((double)time/1000)+"s"+"  loopcount: "+N);*/
	return I;
}
 

/**
 * 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);
}