import org.ejml.simple.SimpleEVD; //导入依赖的package包/类
@SuppressWarnings("rawtypes")
private static SimpleMatrix symmetricDecorrelation(SimpleMatrix x) {
	
	double d;
	SimpleMatrix QL;
	SimpleMatrix Q;
	
	SimpleEVD evd = x.transpose().mult(x).eig();
	int len = evd.getNumberOfEigenvalues();
	QL = new SimpleMatrix(len, len);
	Q = new SimpleMatrix(len, len);
	
	// Scale each column of the eigenvector matrix by the square root of 
	// the reciprocal of the associated eigenvalue
	for (int i = 0; i < len; i++) {
		d = evd.getEigenvalue(i).getReal();
		d = (d + Math.abs(d)) / 2;  // improve numerical stability by eliminating small negatives near singular matrix zeros
		QL.insertIntoThis(0, i, evd.getEigenVector(i).divide(Math.sqrt(d)));
		Q.insertIntoThis(0, i, evd.getEigenVector(i));
	}
	
	return x.mult(QL.mult(Q.transpose()));
}
 

import org.ejml.simple.SimpleEVD; //导入依赖的package包/类
/**
 * Calculates the eigen decomposition of a real matrix. The eigen
 * decomposition of matrix A is a set of two matrices: V and D such that A =
 * V × D × VT. A, V and D are all m × m matrices.
 *
 * @param a Given matrix.
 * @return Result W/V arrays.
 */
public static Array[] eigen_EJML(Array a) {
    int m = a.getShape()[0];
    Array Wa;
    Array Va = Array.factory(DataType.DOUBLE, new int[]{m, m});
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a);
    SimpleMatrix M = new SimpleMatrix(aa);
    SimpleEVD evd = M.eig();
    List<Complex_F64> evs = evd.getEigenvalues();
    boolean isComplex = evd.getEigenVector(0) == null;
    if (isComplex) {
        Wa = Array.factory(DataType.OBJECT, new int[]{m});
    } else {
        Wa = Array.factory(DataType.DOUBLE, new int[]{m});
    }
    for (int i = 0; i < m; i++) {
        if (isComplex) {
            Wa.setObject(i, new Complex(evs.get(i).real, evs.get(i).imaginary));
        } else {
            Wa.setDouble(i, evs.get(m - i - 1).real);
            SimpleBase v = evd.getEigenVector(m - i - 1);
            for (int j = 0; j < v.getNumElements(); j++) {
                Va.setDouble(j * m + i, v.get(j));
            }
        }
    }

    return new Array[]{Wa, Va};
}
 

import org.ejml.simple.SimpleEVD; //导入依赖的package包/类
/**
 * 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);
}
 

import org.ejml.simple.SimpleEVD; //导入依赖的package包/类
/**
 * 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);
}