import org.jblas.Singular; //导入依赖的package包/类
public Matrix[] calc(Matrix source) {
	final DoubleMatrix matrix;
	if (source instanceof JBlasDenseDoubleMatrix2D) {
		matrix = ((JBlasDenseDoubleMatrix2D) source).getWrappedObject();
	} else if (source instanceof HasColumnMajorDoubleArray1D) {
		matrix = new JBlasDenseDoubleMatrix2D(MathUtil.longToInt(source.getRowCount()), MathUtil.longToInt(source
				.getColumnCount()), ((HasColumnMajorDoubleArray1D) source).getColumnMajorDoubleArray1D())
				.getWrappedObject();
	} else {
		matrix = new JBlasDenseDoubleMatrix2D(source).getWrappedObject();
	}
	final DoubleMatrix[] svd = Singular.fullSVD(matrix);
	final Matrix u = new JBlasDenseDoubleMatrix2D(svd[0]);
	DoubleMatrix sVector = svd[1];
	final Matrix s = new JBlasDenseDoubleMatrix2D(MathUtil.longToInt(source.getRowCount()),
			MathUtil.longToInt(source.getColumnCount()));
	for (int i = 0; i < sVector.rows; i++) {
		s.setAsDouble(sVector.get(i), i, i);
	}
	final Matrix v = new JBlasDenseDoubleMatrix2D(svd[2]);
	return new Matrix[] { u, s, v };
}
 

import org.jblas.Singular; //导入依赖的package包/类
protected DoubleMatrix[] computeTruncatedSVD(DoubleMatrix mat, Double singValTol, int maxSize)
{
	DoubleMatrix[] svdResult = Singular.sparseSVD(mat);
	
	notifyPSRObservers(svdResult[1].toArray());
	
	int singIter;
	for(singIter = 0; singIter < svdResult[1].getRows(); singIter++)
	{
		if(svdResult[1].get(singIter,0) < singValTol)
			break;
	}
	
	int numToKeep = Math.min(singIter, maxSize);
	maxDim = numToKeep;
			
	DoubleMatrix u = new DoubleMatrix(mat.getRows(),numToKeep);
	u = svdResult[0].getRange(0,svdResult[0].getRows(), 0, numToKeep);
	DoubleMatrix s = DoubleMatrix.diag(svdResult[1].getRange(0, numToKeep, 0, 1));
	DoubleMatrix v = new DoubleMatrix(mat.getColumns(), numToKeep);
	v = svdResult[2].getRange(0, svdResult[2].getRows(), 0, numToKeep);
	
	DoubleMatrix[] truncatedSVDResult = {u.transpose(),s,v};
	return truncatedSVDResult;
}
 

import org.jblas.Singular; //导入依赖的package包/类
private DoubleMatrix computeMapping(DoubleMatrix x) {
    DoubleMatrix alphas = MatrixFunctions.pow(
            MatrixFunctions.pow(
                    this.highDimControlPoints.subRowVector(x),
                    2
            ).rowSums(),
            -1
    );
    double alpha = alphas.sum();

    DoubleMatrix xTilde = this.highDimControlPoints.mulColumnVector(alphas).columnSums().div(alpha);
    DoubleMatrix yTilde = this.lowDimControlPoints.mulColumnVector(alphas).columnSums().div(alpha);
    DoubleMatrix xHats = this.highDimControlPoints.subRowVector(xTilde);
    @SuppressWarnings("SuspiciousNameCombination")
    DoubleMatrix yHats = this.lowDimControlPoints.subRowVector(yTilde);

    DoubleMatrix sqrtAlphas = MatrixFunctions.sqrt(alphas);
    DoubleMatrix A = xHats.mulColumnVector(sqrtAlphas);
    DoubleMatrix B = yHats.mulColumnVector(sqrtAlphas);

    DoubleMatrix[] svdComposition = Singular.sparseSVD(A.transpose().mmul(B));
    DoubleMatrix U = svdComposition[0];
    DoubleMatrix V = svdComposition[2];

    DoubleMatrix M = U.mmul(V);

    return x.sub(xTilde).mmul(M).add(yTilde);
}
 

import org.jblas.Singular; //导入依赖的package包/类
/**
 * 
 * 
 * @param A
 * @param dim
 * @return
 */
public DoubleMatrix pca(DoubleMatrix A, int dim) {
	
	logger.info("Type 'pca' started (target dimensions = "+dim+").");
	
	A = centerData(A);
	save("A", A);
	
	logger.info("Computing SVD...");
	DoubleMatrix[] usv = Singular.fullSVD(A);
	DoubleMatrix U = usv[0];
	DoubleMatrix S = usv[1];
	
	save("U", U);
	save("S", S);
	
	// 
	logger.info("Reducing U to Uk...");
	DoubleMatrix Uk = new DoubleMatrix(U.rows, dim);
	for(int i=0; i<dim; i++)
		Uk.putColumn(i, U.getColumn(i));
	save("Uk", Uk);
	
	// build S matrix
	logger.info("Reducing S to Sk...");
	DoubleMatrix Sk = new DoubleMatrix(dim, dim);
	for (int i = 0; i < dim; i++) {
		Sk.put(i, i, S.get(i));
	}
	save("Sk", Sk);
	
	// calculate principal component matrix...
	logger.info("Computing principal components...");
	DoubleMatrix B = Uk.mmul(Sk);
	save("B", B);
	
	return B;
	
}
 

import org.jblas.Singular; //导入依赖的package包/类
/**
 * @param A
 * @param k
 * @return
 */
public DoubleMatrix compress(DoubleMatrix A, int k) {
	
	logger.info("Type 'compress' started...");

	A = centerData(A);
	
	DoubleMatrix[] usv = Singular.fullSVD(A);
	// n x n
	DoubleMatrix U = usv[0];
	// n x p
	DoubleMatrix S = usv[1];
	// p x p (straight)
	DoubleMatrix V = usv[2];
			
	// k x k
	DoubleMatrix Sk = new DoubleMatrix(U.columns, V.columns);
	for (int i = 0; i < k; i++) {
		Sk.put(i, i, S.get(i));
	}
	
	// 
	DoubleMatrix Aapprox = U.mmul(Sk).mmul(V.transpose());
	save("Aapprox", Aapprox);
	
	return Aapprox;
	
}
 

import org.jblas.Singular; //导入依赖的package包/类
public DoubleMatrix reconstruct2(DoubleMatrix A, int k) {
		
		logger.info("Type 'reconstruct2' started...");

		A = centerData(A);
		save("A", A);
		
		DoubleMatrix[] usv = Singular.fullSVD(A);
		// n x n
		DoubleMatrix U = usv[0];
		// n x p
		DoubleMatrix S = usv[1];
		// p x p (straight)
		DoubleMatrix V = usv[2];
		
		save("U", U);
		save("S", S);
		save("V", V);
		
		for(int i = k; i < U.columns; i++)
			U.putColumn(i, DoubleMatrix.zeros(U.rows, 1));
		save("Uk", U);
		
		DoubleMatrix Sm = new DoubleMatrix(A.rows, A.columns);
		for (int i = 0; i < k; i++) {
			Sm.put(i, i, S.get(i));
		}
		save("Sk", Sm);
		
//		for (int i = k; i < V.rows; i++) {
//			V.putColumn(i, DoubleMatrix.zeros(1, V.columns));
//		}
//		visual("Vk", V);
 		
		// 
		DoubleMatrix Aapprox = U.mmul(Sm).mmul(V.transpose());
		save("Aapprox", Aapprox);
				
		return Aapprox;
	}
 

import org.jblas.Singular; //导入依赖的package包/类
/**
 * 初始化PCA
 * 
 * @param x 特征矩阵
 */
public PCA(FloatMatrix x) {
	super();
	this.x = x;

	this.m = x.rows;
	FloatMatrix Sigma = x.transpose().mmul(x).div(m);
	FloatMatrix[] result = Singular.fullSVD(Sigma);

	this.u = result[0].neg();// 与matlab保持一致
	this.s = result[1];
}
 

import org.jblas.Singular; //导入依赖的package包/类
public static DoubleMatrix pinv(final DoubleMatrix A) {
  final DoubleMatrix[] usv = Singular.fullSVD(A);
  final DoubleMatrix s_inv = new DoubleMatrix(usv[1].length);
  for (int i = 0; i < usv[1].length; i++) {
    s_inv.put(i, Math.abs(usv[1].get(i)) < TOLERANCE ? 0.0 : 1.0 / usv[1].get(i));
  }
  return usv[2].mmul(DoubleMatrix.diag(s_inv)).mmul(usv[0].transpose());
}
 

import org.jblas.Singular; //导入依赖的package包/类
public static DoubleMatrix[] svd(DoubleMatrix src)
{
    DoubleMatrix[] usv = Singular.fullSVD(src);
       DoubleMatrix S = usv[1];
       usv[1] = convertDiagMatrix(S, src.rows, src.columns);
       return usv;
}
 

import org.jblas.Singular; //导入依赖的package包/类
@Test
public void testSVD()
{
	double[][] testData = { 
			{ 36, 49, 47, 11 }, 
			{ 2, 68, 27, 42 }, 
			{ 42, 25, 38, 3 } 
		};
	RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
	SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
       System.out.println(svd.getU());
       System.out.println(svd.getS());
       System.out.println(svd.getV());
	
	DoubleMatrix[] usv = Singular.fullSVD(new DoubleMatrix(testData));
       System.out.println(usv[0]);
       System.out.println(usv[1]);
       System.out.println(usv[2]);

       DoubleMatrix U = usv[0];
       DoubleMatrix S = usv[1];
       DoubleMatrix V = usv[2];
       DoubleMatrix mt = new DoubleMatrix(3, 4);
       for (int i = 0; i < S.length; i++)
       {
           mt.put(i, i, S.get(i));
       }
       System.out.println(mt.toString().replace(";", "\n"));
       DoubleMatrix src = U.mmul(mt).mmul(V.transpose());
       System.out.println(src.toString().replace(";", "\n"));
       mt = Solve.pinv(mt);
       System.out.println(mt.toString().replace(";", "\n"));
}
 

import org.jblas.Singular; //导入依赖的package包/类
/**
 * @param A
 * @param k
 * @return
 */
public DoubleMatrix reconstruct(DoubleMatrix A, int k) {
	
	logger.info("Type 'reconstruct' started...");

	A = centerData(A);
	save("A", A);
	
	DoubleMatrix[] usv = Singular.fullSVD(A);
	// n x n
	DoubleMatrix U = usv[0];
	// n x p
	DoubleMatrix S = usv[1];
	// p x p (straight)
	DoubleMatrix V = usv[2];
	
	save("U", U);
	save("S", S);
	save("V", V);
	
	// n x k
	DoubleMatrix Uk = new DoubleMatrix(U.rows, k);
	for(int i=0; i<k; i++)
		Uk.putColumn(i, U.getColumn(i));
	save("Uk", Uk);
	
	// k x k
	DoubleMatrix Sk = new DoubleMatrix(k, k);
	for (int i = 0; i < k; i++) {
		Sk.put(i, i, S.get(i));
	}
	save("Sk", Sk);
	
	// p x k (straight)
	DoubleMatrix Vk = new DoubleMatrix(A.columns, k);
	for (int i = 0; i < k; i++) {
		Vk.putColumn(i, V.getColumn(i));
	}
	save("Vk", Vk);
		
	// 
	DoubleMatrix Aapprox = Uk.mmul(Sk).mmul(Vk.transpose());
	save("Aapprox", Aapprox);
	
	DoubleMatrix Areduced = new DoubleMatrix(A.rows, k);
	for(int i=0; i<k; i++)
		Areduced.putColumn(i, Aapprox.getColumn(i));
	save("C"+k, Areduced);
	
	return Areduced;
	
}
 

import org.jblas.Singular; //导入依赖的package包/类
public void train(float[][] xraw, float[][] yraw) {
	this.xOrig = new float[xraw.length][];
	for (int i=0; i<xraw.length; ++i) {
		this.xOrig[i] = a.append(1.0f, xraw[i]);
	}
	this.yOrig = yraw;
	
	long start = System.nanoTime();
	this.xMean = a.scale(a.sum(a.transpose(xOrig)), 1.0f / xOrig.length);
	FloatMatrix meanMat = new FloatMatrix(xMean);
	FloatMatrix XMat = new FloatMatrix(xOrig);
	FloatMatrix XMatSubMean = XMat.subRowVector(meanMat);
	FloatMatrix covXMat = XMatSubMean.transpose().mmul(XMatSubMean);
	FloatMatrix WTransMat = Singular.fullSVD(covXMat)[0].transpose();
	System.out.println("PCA time: "+(System.nanoTime() - start) / 1e6);

	start = System.nanoTime();
	this.projDirection = Arrays.copyOfRange(WTransMat.toArray2(), 0, K);
	this.proj = new float[K][];
	this.x = new float[K][][];
	this.y = new float[K][][];
	for (int k=0; k<K; ++k) {
		final int kfinal = k;
		this.proj[k] = new float[xOrig.length];
		for (int i=0; i<xOrig.length; ++i) {
			proj[k][i] = a.innerProd(projDirection[k], a.comb(xOrig[i], 1.0f, xMean, -1.0f));
		}
		Integer[] indices = new Integer[xOrig.length];
		for (int i=0; i<xOrig.length; ++i) indices[i] = i;
		Arrays.sort(indices, new Comparator<Integer>() {
			public int compare(Integer i1, Integer i2) {
				if (proj[kfinal][i1] < proj[kfinal][i2]) {
					return -1;
				} else if (proj[kfinal][i1] > proj[kfinal][i2]) {
					return 1;
				} else {
					return 0;
				}
			}
		});
		float[] projOld = proj[k];
		x[k] = new float[xOrig.length][];
		y[k] = new float[yOrig.length][];
		proj[k] = new float[projOld.length];
		for (int i=0; i<xOrig.length; ++i) {
			x[k][i] = xOrig[indices[i]];
			y[k][i] = yOrig[indices[i]];
			proj[k][i] = projOld[indices[i]];
		}
	}
	System.out.println("Project / sort time: "+(System.nanoTime() - start) / 1e6);
	
	this.xMean_d = Matrix.build(1, xMean.length, xMean);
	this.xMean_d.setDontFree(true);
	this.x_d = new Matrix[K];
	this.y_d = new Matrix[K];
	this.xSqrNorms_d = new Matrix[K];
	this.projDirection_d = new Matrix[K];
	this.proj_d = new Matrix[K];
	for (int k=0; k<K; ++k) {
		this.x_d[k] = Matrix.build(x[k]);
		this.x_d[k].setDontFree(true);
		this.y_d[k] = Matrix.build(y[k]);
		this.y_d[k].setDontFree(true);
		this.xSqrNorms_d[k] = x_d[k].sqr().colSum();
		this.xSqrNorms_d[k].setDontFree(true);
		this.projDirection_d[k] = Matrix.build(projDirection[k].length, 1, projDirection[k]);
		this.projDirection_d[k].setDontFree(true);
		this.proj_d[k] = Matrix.build(proj[k].length, 1, proj[k]);
		this.proj_d[k].setDontFree(true);
	}
}
 

import org.jblas.Singular; //导入依赖的package包/类
public void train(float[][] xraw, float[][] yraw) {
	this.x = new float[xraw.length][];
	for (int i=0; i<xraw.length; ++i) {
		this.x[i] = a.append(1.0f, xraw[i]);
	}
	this.y = yraw;
	
	long start = System.nanoTime();
	this.xMean = a.scale(a.sum(a.transpose(x)), 1.0f / x.length);
	FloatMatrix meanMat = new FloatMatrix(xMean);
	FloatMatrix XMat = new FloatMatrix(x);
	FloatMatrix XMatSubMean = XMat.subRowVector(meanMat);
	FloatMatrix covXMat = XMatSubMean.transpose().mmul(XMatSubMean);
	FloatMatrix WTransMat = Singular.fullSVD(covXMat)[0].transpose();
	System.out.println("PCA time: "+(System.nanoTime() - start) / 1e6);

	start = System.nanoTime();
	this.projDirection = WTransMat.toArray2()[0];
	this.proj = new float[x.length];
	for (int i=0; i<x.length; ++i) {
		proj[i] = a.innerProd(projDirection, a.comb(x[i], 1.0f, xMean, -1.0f));
	}
	Integer[] indices = new Integer[x.length];
	for (int i=0; i<x.length; ++i) indices[i] = i;
	Arrays.sort(indices, new Comparator<Integer>() {
		public int compare(Integer i1, Integer i2) {
			if (proj[i1] < proj[i2]) {
				return -1;
			} else if (proj[i1] > proj[i2]) {
				return 1;
			} else {
				return 0;
			}
		}
	});
	float[][] xold = x;
	float[][] yold = y;
	float[] projOld = proj;
	x = new float[xold.length][];
	y = new float[yold.length][];
	proj = new float[projOld.length];
	for (int i=0; i<x.length; ++i) {
		x[i] = xold[indices[i]];
		y[i] = yold[indices[i]];
		proj[i] = projOld[indices[i]];
	}
	System.out.println("Project / sort time: "+(System.nanoTime() - start) / 1e6);

	this.x_d = Matrix.build(x);
	this.y_d = Matrix.build(y);
	this.xSqrNorms_d = x_d.sqr().colSum();
	this.xMean_d = Matrix.build(1, xMean.length, xMean);
	this.projDirection_d = Matrix.build(projDirection.length, 1, projDirection);
	this.proj_d = Matrix.build(proj.length, 1, proj);
	
	this.x_d.setDontFree(true);
	this.y_d.setDontFree(true);
	this.xSqrNorms_d.setDontFree(true);
	this.xMean_d.setDontFree(true);
	this.projDirection_d.setDontFree(true);
	this.proj_d.setDontFree(true);
}
 

import org.jblas.Singular; //导入依赖的package包/类
public CPSR(TrainingDataSet trainData, double minSingularVal, int maxSVDDim, int projDim, int maxHistLen, ProjType projType, boolean histCompress, boolean randStart)
{
	super(trainData, minSingularVal, maxSVDDim, maxHistLen);
	this.projType = projType;
	this.histCompress = histCompress;
	this.projDim = projDim;

	tPhi = constructRandomMatrix(projType, projDim, tests.size(), projDim);
	if(!histCompress)
	{
		hPhi = constructRandomMatrix(ProjType.Eye,histories.size()+1,histories.size()+1, 1.0);
	}
	else
	{
		hPhi = constructRandomMatrix(projType,histories.size()+1,projDim, projDim);
	}

	DoubleMatrix one = DoubleMatrix.zeros(histories.size()+1, 1);
	this.randStart = randStart;

	one.put(0,0, 1.0);
	if(randStart)
	{
		for(int i = 0; i < one.getRows(); i++)
			one.put(i,0,1.0);
	}

	if(histCompress)
	{
		DoubleMatrix[] hPhiSVD = Singular.sparseSVD(hPhi);
		DoubleMatrix sigma = DoubleMatrix.diag(hPhiSVD[1]);


		e = hPhiSVD[2].mmul(Solve.pinv(sigma)).mmul(hPhiSVD[0].transpose());
		
		e = e.mmul(one);
	}
	else
	{
		e = one;
	}
	
	
}
 

import org.jblas.Singular; //导入依赖的package包/类
protected DoubleMatrix[] updateSVD(DoubleMatrix newTHData)
{
	DoubleMatrix biggerV = DoubleMatrix.concatVertically(svdResults[2],
			DoubleMatrix.zeros(newTHData.getColumns()-svdResults[2].getRows(), svdResults[2].getColumns()));
	
	DoubleMatrix m = svdResults[0].mmul(newTHData);
	DoubleMatrix p = newTHData.sub((svdResults[0].transpose()).mmul(m));
	
	DoubleMatrix pBase = Singular.sparseSVD(p)[0];
	pBase = DoubleMatrix.concatHorizontally(pBase, DoubleMatrix.zeros(pBase.getRows(), p.getColumns()-pBase.getColumns()));
	DoubleMatrix rA = (pBase.transpose()).mmul(p);
	
	DoubleMatrix n = (biggerV.transpose()).mmul(DoubleMatrix.eye(newTHData.getColumns()));
	DoubleMatrix q = DoubleMatrix.eye(newTHData.getColumns()).sub(biggerV.mmul(n));
	
	DoubleMatrix qBase = Singular.sparseSVD(q)[0];
	DoubleMatrix rB = (qBase.transpose()).mmul(q);
	
	DoubleMatrix z = DoubleMatrix.zeros(m.getRows(), m.getColumns());
	DoubleMatrix z2 = DoubleMatrix.zeros(m.getColumns(),m.getColumns());
	
	DoubleMatrix top = DoubleMatrix.concatHorizontally(svdResults[1],z);
	DoubleMatrix bottom = DoubleMatrix.concatHorizontally(z.transpose(),z2);
	DoubleMatrix comb = DoubleMatrix.concatVertically(top, bottom);
	
	DoubleMatrix mRA = DoubleMatrix.concatVertically(m, rA);
	DoubleMatrix nRB = DoubleMatrix.concatVertically(n, rB);

	DoubleMatrix mRANRB = mRA.mmul(nRB.transpose());
	DoubleMatrix k = comb.add(mRANRB);
	
	DoubleMatrix[] newBases = computeTruncatedSVD(k, Double.MIN_VALUE, maxDim);
			
	DoubleMatrix uP = (DoubleMatrix.concatHorizontally(svdResults[0].transpose(), p)).mmul(newBases[0].transpose());
	DoubleMatrix vP = (DoubleMatrix.concatHorizontally(biggerV,q)).mmul(newBases[2]);
	
	DoubleMatrix[] updatedSVDResults = new DoubleMatrix[3];
	
	updatedSVDResults[0] = uP.transpose();
	updatedSVDResults[1] = newBases[1];
	updatedSVDResults[2] = vP;
	
	double[] singVals = new double[updatedSVDResults[1].getRows()];
	for(int i = 0; i < updatedSVDResults[1].getRows(); i++)
		singVals[i] = updatedSVDResults[1].get(i,i);
			
	return updatedSVDResults;
}