cern.colt.matrix.linalg.Algebra Java Examples

@Override
public void aggregate(double[][][] values) {
	result = null;
	int offset = super.calculateOffset(values,featureNameIndecis);
	int[][] featureIndecis = super.collapseFeatures(values,featureNameIndecis);
	result[0] = 0.0;
	windowLength = featureNameIndecis.length-offset;
	featureLength = featureIndecis[0].length;
	for (int i=offset;i<values.length;++i){
		for(int j=0;j<featureIndecis.length;++j){
			result[0] += values[i][featureIndecis[j][0]][featureIndecis[j][1]];
		}
	}		
	terms = new DenseDoubleMatrix2D(xDim*yDim,windowLength*featureLength);
	z = new DenseDoubleMatrix2D(1,featureLength);
	calcTerms(terms);
	result = ((new Algebra()).solve(terms,z)).viewRow(0).toArray();
}
 

/**
 */
public static void doubleTest34() {
	double[][] data = 
	{ 
		{ 3, 0, 0, 0 },
		{ 0, 4, 2, 0 },
		{ 0, 0, 0, 0 },
		{ 0, 0, 0, 0 },
	};
	
	DoubleMatrix2D A = new DenseDoubleMatrix2D(data);
	Property.DEFAULT.generateNonSingular(A);
	DoubleMatrix2D inv = Algebra.DEFAULT.inverse(A);


	System.out.println("\n\n\n"+A);
	System.out.println("\n"+inv);
	DoubleMatrix2D B = A.zMult(inv,null);
	System.out.println(B);
	if (!(B.equals(DoubleFactory2D.dense.identity(A.rows)))) {
		throw new InternalError();
	}
}
 

/**
 */
public static void doubleTest34() {
	double[][] data = 
	{ 
		{ 3, 0, 0, 0 },
		{ 0, 4, 2, 0 },
		{ 0, 0, 0, 0 },
		{ 0, 0, 0, 0 },
	};
	
	DoubleMatrix2D A = new DenseDoubleMatrix2D(data);
	Property.DEFAULT.generateNonSingular(A);
	DoubleMatrix2D inv = Algebra.DEFAULT.inverse(A);


	System.out.println("\n\n\n"+A);
	System.out.println("\n"+inv);
	DoubleMatrix2D B = A.zMult(inv,null);
	System.out.println(B);
	if (!(B.equals(DoubleFactory2D.dense.identity(A.rows)))) {
		throw new InternalError();
	}
}
 

public PepPlaneLinModel(Residue res1, Residue res2){
    //construct from the current conformations
    //of the residue involved in this peptide plane
    
    if(res1.template.name.equalsIgnoreCase("PRO") || res2.template.name.equalsIgnoreCase("PRO"))
        throw new RuntimeException("ERROR: CATS doesn't handle proline");
    
    double CA1[] = res1.getCoordsByAtomName("CA");
    double CA2[] = res2.getCoordsByAtomName("CA");
    double N[] = res2.getCoordsByAtomName("N");
    double C[] = res1.getCoordsByAtomName("C");
    double O[] = res1.getCoordsByAtomName("O");
    double H[] = res2.getCoordsByAtomName("H");
    
    DoubleMatrix2D M = makeAxisMatrix(CA1,N,CA2);
    DoubleMatrix2D vec = DoubleFactory2D.dense.make(3,1);
    
    vec.viewColumn(0).assign(VectorAlgebra.subtract(C, CA1));
    CCoeffs = Algebra.DEFAULT.solve(M, vec).viewColumn(0);
    
    vec.viewColumn(0).assign(VectorAlgebra.subtract(O, CA1));
    OCoeffs = Algebra.DEFAULT.solve(M, vec).viewColumn(0);
    
    vec.viewColumn(0).assign(VectorAlgebra.subtract(H, CA1));
    HCoeffs = Algebra.DEFAULT.solve(M, vec).viewColumn(0);
}
 

public double[] transformMatrix(double[][] inputMatrix, int dim) {
    Algebra algebra = new Algebra();

    DoubleMatrix2D H = new DenseDoubleMatrix2D(inputMatrix);
    RobustEigenDecomposition decomposition = new RobustEigenDecomposition(H);
    DoubleMatrix1D eigenvalues = decomposition.getRealEigenvalues();

    normalizeEigenvalues(eigenvalues);

    DoubleMatrix2D V = decomposition.getV();

    transformEigenvalues(eigenvalues);

    double[][] negativeHessianInverse = algebra.mult(
            algebra.mult(V, DoubleFactory2D.dense.diagonal(eigenvalues)),
            algebra.inverse(V)
    ).toArray();

    double[] massArray = new double[dim * dim];
    for (int i = 0; i < dim; i++) {
        System.arraycopy(negativeHessianInverse[i], 0, massArray, i * dim, dim);
    }
    return massArray;
}
 

@Override
public double train(List<Minimizer.Result> trainingSet, List<Minimizer.Result> testSet) {

	// make sure all the samples have the right dimension
	for (Minimizer.Result sample : trainingSet) {
		if (sample.dofValues.size() != numDofs) {
			throw new IllegalArgumentException("samples have wrong number of dimensions");
		}
	}

	LinearSystem trainingSystem = new LinearSystem(trainingSet);
	LinearSystem testSystem = new LinearSystem(testSet);

	// solve Ax = b in least squares sense
	coefficients.assign(new QRDecomposition(trainingSystem.A).solve(trainingSystem.b).viewColumn(0));

	// calculate the residual (Ax - b) for the test set
	DoubleMatrix1D residual = new Algebra()
		.mult(testSystem.A, coefficients)
		.assign(testSystem.b.viewColumn(0), (ri, bi) -> Math.abs(ri - bi));

	// calculate the max error
	maxe = 0.0;
	for (int i=0; i<residual.size(); i++) {
		assert (residual.get(i) >= 0.0);
		maxe = Math.max(maxe, residual.get(i));
	}

	return maxe;
}
 

@Test
void givenTwoMatrices_whenMultiply_thenMultiplicatedMatrix() {
    DoubleFactory2D doubleFactory2D = DoubleFactory2D.dense;

    DoubleMatrix2D firstMatrix = doubleFactory2D.make(
      new double[][] {
        new double[] {1d, 5d},
        new double[] {2d, 3d},
        new double[] {1d ,7d}
      }
    );

    DoubleMatrix2D secondMatrix = doubleFactory2D.make(
      new double[][] {
        new double[] {1d, 2d, 3d, 7d},
        new double[] {5d, 2d, 8d, 1d}
      }
    );

    DoubleMatrix2D expected = doubleFactory2D.make(
      new double[][] {
        new double[] {26d, 12d, 43d, 12d},
        new double[] {17d, 10d, 30d, 17d},
        new double[] {36d, 16d, 59d, 14d}
      }
    );

    Algebra algebra = new Algebra();
    DoubleMatrix2D actual = algebra.mult(firstMatrix, secondMatrix);

    assertThat(actual).isEqualTo(expected);
}
 

@Benchmark
public Object coltMatrixMultiplication(MatrixProvider matrixProvider) {
    DoubleFactory2D doubleFactory2D = DoubleFactory2D.dense;

    DoubleMatrix2D firstMatrix = doubleFactory2D.make(matrixProvider.getFirstMatrix());
    DoubleMatrix2D secondMatrix = doubleFactory2D.make(matrixProvider.getSecondMatrix());

    Algebra algebra = new Algebra();
    return algebra.mult(firstMatrix, secondMatrix);
}
 

@Benchmark
public Object coltMatrixMultiplication(BigMatrixProvider matrixProvider) {
    DoubleFactory2D doubleFactory2D = DoubleFactory2D.dense;

    DoubleMatrix2D firstMatrix = doubleFactory2D.make(matrixProvider.getFirstMatrix());
    DoubleMatrix2D secondMatrix = doubleFactory2D.make(matrixProvider.getSecondMatrix());

    Algebra algebra = new Algebra();
    return algebra.mult(firstMatrix, secondMatrix);
}
 

/**
 */
public static void doubleTest27() {

	
System.out.println("\n\n");
System.out.println("initializing...");

int rows=51;
int columns=10;
double[][] trainingSet = new double[columns][rows];
for (int i=columns; --i >= 0; ) trainingSet[i][i]=2.0;

int patternIndex        = 0;
int unitIndex           = 0;

DoubleMatrix2D patternMatrix       = null;
DoubleMatrix2D transposeMatrix     = null;
DoubleMatrix2D QMatrix             = null;
DoubleMatrix2D inverseQMatrix      = null;
DoubleMatrix2D pseudoInverseMatrix = null;
DoubleMatrix2D weightMatrix        = null;

// form a matrix with the columns as training vectors
patternMatrix = DoubleFactory2D.dense.make (rows, columns);

// copy the patterns into the matrix
for (patternIndex=0;patternIndex<columns;patternIndex++) {
	for (unitIndex=0;unitIndex<rows;unitIndex++) {
		patternMatrix.setQuick (unitIndex, patternIndex, trainingSet[patternIndex][unitIndex]);
	}
}

transposeMatrix     = Algebra.DEFAULT.transpose (patternMatrix);
QMatrix             = Algebra.DEFAULT.mult (transposeMatrix,patternMatrix);
inverseQMatrix      = Algebra.DEFAULT.inverse (QMatrix);
pseudoInverseMatrix = Algebra.DEFAULT.mult (inverseQMatrix,transposeMatrix);
weightMatrix        = Algebra.DEFAULT.mult (patternMatrix,pseudoInverseMatrix);
System.out.println("done.");
	
}
 

/**
 */
public static void doubleTest26(int size) {

	
System.out.println("\n\n");
System.out.println("initializing...");
boolean dense = true;
DoubleMatrix2D A;
DoubleFactory2D factory;
if (dense) 
	factory = Factory2D.dense;
else 
	factory = Factory2D.sparse;
	
double value = 0.5;
A = factory.make(size,size,value);
Property.generateNonSingular(A);
cern.colt.Timer timer = new cern.colt.Timer().start();

DoubleMatrix2DComparator fun = new DoubleMatrix2DComparator() {
	public int compare(DoubleMatrix2D a, DoubleMatrix2D b) {
		return a.zSum() == b.zSum() ? 1 : 0;
	}
};

System.out.println(A);
System.out.println(Algebra.ZERO.inverse(A));

timer.stop().display();

System.out.println("done.");
	
}
 

/**
 */
public static void doubleTest25(int size) {

	
System.out.println("\n\n");
System.out.println("initializing...");
boolean dense = true;
DoubleMatrix2D A;
DoubleFactory2D factory;
if (dense) 
	factory = Factory2D.dense;
else 
	factory = Factory2D.sparse;
	
double value = 0.5;
A = factory.make(size,size,value);
Property.generateNonSingular(A);
cern.colt.Timer timer = new cern.colt.Timer().start();

System.out.println(A);
System.out.println(Algebra.ZERO.inverse(A));

timer.stop().display();

System.out.println("done.");
	
}
 

public static void main(String[] args) {
	DoubleMatrix1D x1 = DoubleFactory1D.dense
			.make(new double[] { 1.0, 2.0});
	DoubleMatrix1D x2 = DoubleFactory1D.dense
			.make(new double[] { 1.0, -2.0});
	DoubleMatrix1D x3 = DoubleFactory1D.dense.make(new double[] { -1.0, -2.0});

	System.out.println(Algebra.DEFAULT.normInfinity(x1));
	System.out.println(Algebra.DEFAULT.normInfinity(x2));
	System.out.println(Algebra.DEFAULT.normInfinity(x3));
}
 

public EPolyPC( EPoly template, int fullOrder, int PCOrder, double PCFac) {
    //set up the principal component basis and the DOF information
    //coeffs will be added later
    //since we will be using this EPolyPC object when performing the fitting
    
    super(template.numDOFs,template.DOFmax,template.DOFmin,template.center,
            template.minE,null,fullOrder,template.DOFNames);
    
    //get the coordinate transformation into the eigenbasis of the template's Hessian
    //so we can use the high-eigenvalue directions as our principcal components
    DoubleMatrix2D hess = SeriesFitter.getHessian(template.coeffs,numDOFs,false);
    EigenvalueDecomposition edec = new EigenvalueDecomposition(hess);

    DoubleMatrix1D eigVals = edec.getRealEigenvalues();
    DoubleMatrix2D eigVecs = edec.getV();//The columns of eigVec are the eigenvectors

    invAxisCoeffs = eigVecs;//axisCoeffs' inverse is its transpose, since it's orthogonal
    axisCoeffs = Algebra.DEFAULT.transpose(eigVecs);
    

    this.fullOrder = fullOrder;
    this.PCOrder = PCOrder;


    //Now figure out what PC's to use
    //Let's try the criterion to use all PCs
    //whose eigenvalues' absolute values are within a factor of PCFac
    //of the biggest

    double maxEigVal = 0;
    for(int n=0; n<numDOFs; n++)
        maxEigVal = Math.max( Math.abs(eigVals.get(n)), maxEigVal );

    isPC = new boolean[numDOFs];

    for(int n=0; n<numDOFs; n++){
        if( Math.abs(eigVals.get(n)) >= maxEigVal*PCFac )
            isPC[n] = true;
    }
}
 

double getScaledDistFromCenter(DoubleMatrix1D x){
    //How far are we from the center of the ellipse, relative to the boundary?
    
    if(x==null)
        System.out.println("woah...");
    
    DoubleMatrix1D xrel = x.copy().assign(nc,Functions.plus);
    return xrel.zDotProduct( Algebra.DEFAULT.mult(A, xrel) );
}
 

public double[] getCartesianCoords(double[] polars) {
	if (polars.length==0) { return new double[0]; }
	
	EigenvalueDecomposition evd = new EigenvalueDecomposition(A);
	DoubleMatrix2D Q = evd.getV();
	DoubleMatrix2D L = evd.getD();
	DoubleMatrix2D qT = Q.viewDice().copy();
	Algebra alg = new Algebra();		
	
	int n = polars.length;		
	double radius = polars[0];
	double[] phis = new double[n-1];
	for (int i=1; i<n; i++) { phis[i-1] = polars[i]; }
	
	double[] cartCoords = new double[n];
	for (int i=0; i<n; i++) {
		double prod = 1;
		for (int j=0; j<i; j++) { prod = prod * Math.sin(phis[j]); }
		if (i<n-1) { prod = prod * Math.cos(phis[i]); } 			
		cartCoords[i] = radius*prod;
	}
	
	// transform cartesian coordinates back!
	
	double[] u = new double[cartCoords.length];
	for (int i=0; i<u.length; i++) {
		u[i] = cartCoords[i]*Math.sqrt(L.get(i, i));
	}
	double[] s = alg.mult(Q, DoubleFactory1D.dense.make(u)).toArray();
	double[] x = new double[s.length];
	for (int i=0; i<x.length; i++) { x[i] = s[i] + c.get(i); }
	
	return x;
}
 

public static void main(String[] args) {
	DoubleMatrix1D x1 = DoubleFactory1D.dense
			.make(new double[] { 1.0, 2.0});
	DoubleMatrix1D x2 = DoubleFactory1D.dense
			.make(new double[] { 1.0, -2.0});
	DoubleMatrix1D x3 = DoubleFactory1D.dense.make(new double[] { -1.0, -2.0});

	System.out.println(Algebra.DEFAULT.normInfinity(x1));
	System.out.println(Algebra.DEFAULT.normInfinity(x2));
	System.out.println(Algebra.DEFAULT.normInfinity(x3));
}
 

/**
 */
public static void doubleTest25(int size) {

	
System.out.println("\n\n");
System.out.println("initializing...");
boolean dense = true;
DoubleMatrix2D A;
DoubleFactory2D factory;
if (dense) 
	factory = Factory2D.dense;
else 
	factory = Factory2D.sparse;
	
double value = 0.5;
A = factory.make(size,size,value);
Property.generateNonSingular(A);
cern.colt.Timer timer = new cern.colt.Timer().start();

System.out.println(A);
System.out.println(Algebra.ZERO.inverse(A));

timer.stop().display();

System.out.println("done.");
	
}
 

/**
 */
public static void doubleTest26(int size) {

	
System.out.println("\n\n");
System.out.println("initializing...");
boolean dense = true;
DoubleMatrix2D A;
DoubleFactory2D factory;
if (dense) 
	factory = Factory2D.dense;
else 
	factory = Factory2D.sparse;
	
double value = 0.5;
A = factory.make(size,size,value);
Property.generateNonSingular(A);
cern.colt.Timer timer = new cern.colt.Timer().start();

DoubleMatrix2DComparator fun = new DoubleMatrix2DComparator() {
	public int compare(DoubleMatrix2D a, DoubleMatrix2D b) {
		return a.zSum() == b.zSum() ? 1 : 0;
	}
};

System.out.println(A);
System.out.println(Algebra.ZERO.inverse(A));

timer.stop().display();

System.out.println("done.");
	
}
 

/**
 */
public static void doubleTest27() {

	
System.out.println("\n\n");
System.out.println("initializing...");

int rows=51;
int columns=10;
double[][] trainingSet = new double[columns][rows];
for (int i=columns; --i >= 0; ) trainingSet[i][i]=2.0;

int patternIndex        = 0;
int unitIndex           = 0;

DoubleMatrix2D patternMatrix       = null;
DoubleMatrix2D transposeMatrix     = null;
DoubleMatrix2D QMatrix             = null;
DoubleMatrix2D inverseQMatrix      = null;
DoubleMatrix2D pseudoInverseMatrix = null;
DoubleMatrix2D weightMatrix        = null;

// form a matrix with the columns as training vectors
patternMatrix = DoubleFactory2D.dense.make (rows, columns);

// copy the patterns into the matrix
for (patternIndex=0;patternIndex<columns;patternIndex++) {
	for (unitIndex=0;unitIndex<rows;unitIndex++) {
		patternMatrix.setQuick (unitIndex, patternIndex, trainingSet[patternIndex][unitIndex]);
	}
}

transposeMatrix     = Algebra.DEFAULT.transpose (patternMatrix);
QMatrix             = Algebra.DEFAULT.mult (transposeMatrix,patternMatrix);
inverseQMatrix      = Algebra.DEFAULT.inverse (QMatrix);
pseudoInverseMatrix = Algebra.DEFAULT.mult (inverseQMatrix,transposeMatrix);
weightMatrix        = Algebra.DEFAULT.mult (patternMatrix,pseudoInverseMatrix);
System.out.println("done.");
	
}
 

private DoubleMatrix2D makeAxisMatrix(double[] CA1, double[] N, double CA2[]){
    //matrix of the axes used to expand our atom coordinates
    double axis1[] = VectorAlgebra.subtract(N, CA1);
    double axis2[] = VectorAlgebra.subtract(CA2, CA1);
    double axis3[] = VectorAlgebra.cross(axis1, axis2);
            
    DoubleMatrix2D M = DoubleFactory2D.dense.make(new double[][] {axis1,axis2,axis3});
    return Algebra.DEFAULT.transpose(M);
}
 

public LinearMultivariateRealFunction[] getJoptPolytope(RCTuple conf, ArrayList<DegreeOfFreedom> DOFs,
        double ineqTol){
    //get the polytope for conf as non-redundant LinearMultivariateRealFunctions
    //relax them all by ineqTol
    final double redundancyTol = 1e-6;//DEBUG!! for checking redundancy of constraints.  From SQPMinimizer
    ArrayList<LinearConstraint> constrList = getFullStericPolytope(conf, DOFs);
    
    ArrayList<LinearMultivariateRealFunction> ajlnv = new ArrayList<>();
    for(LinearConstraint c : constrList){
        LinearMultivariateRealFunction f = LPChecks.toLinearMultivariateRealFunction(c,ineqTol);
        boolean redundant = false;
        for(LinearMultivariateRealFunction g : ajlnv){
            if(Math.abs(f.getR()-g.getR())<redundancyTol){
                if(Algebra.DEFAULT.norm2(f.getQ().copy().assign(g.getQ(),Functions.minus))<redundancyTol*redundancyTol){
                    redundant = true;
                    break;
                }
                //DEBUG!!  This is somewhat redundant with constr cache in getFullStericPolytope
                //but it will get non-voxel redundant constr (intra constr will appears as (at1,at2) and (at2,at1))
            }
        }
        if(!redundant)
            ajlnv.add(f);
    }
    
    LinearMultivariateRealFunction constr[] = ajlnv.toArray(new LinearMultivariateRealFunction[ajlnv.size()]);
    return constr;
}
 

public HashMap<DegreeOfFreedom,Double> calcSpecialCenter(){
    //center for some degrees of freedom may, due to non-box constr, 
    //differ from center of box constr
    //NOTE: the map only contains those degrees of freedom where center is special
    //(not just middle of lb,ub)
    //DEBUG!!! This assumes the particular form of lin constr (pairs of parallel constr,
    //DOFs not in constrained space assumed to center at 0) used in basis big CATS stuff
    
    HashMap<DegreeOfFreedom,Double> ans = new HashMap<>();
    int numSpecialDOFConstr = voxLinConstr.size()/2;
    if(numSpecialDOFConstr==0)
        return ans;
    
    HashSet<Integer> specialDOFSet = new HashSet<>();
    for(int spesh=0; spesh<numSpecialDOFConstr; spesh++){
        RealVector coeff = voxLinConstr.get(2*spesh).getCoefficients();
        if(coeff.getDistance(voxLinConstr.get(2*spesh+1).getCoefficients()) > 0)
            throw new RuntimeException("ERROR: voxLinConstr not paired like expected");
        
        for(int d=0; d<dofIntervals.size(); d++){
            if(coeff.getEntry(d)!=0)
                specialDOFSet.add(d);
        }
    }
    
    int numSpecialDOFs = specialDOFSet.size();
    ArrayList<Integer> specialDOFList = new ArrayList<>();
    specialDOFList.addAll(specialDOFSet);
    
    DoubleMatrix2D specialConstr = DoubleFactory2D.dense.make(numSpecialDOFConstr, numSpecialDOFs);
    //we will constrain values first of the lin combs of DOFs given in voxLinConstr,
    //then of lin combs orth to those
    DoubleMatrix2D target = DoubleFactory2D.dense.make(numSpecialDOFs,1);
    for(int spesh=0; spesh<numSpecialDOFConstr; spesh++){
        RealVector coeffs = voxLinConstr.get(2*spesh).getCoefficients();
        for(int d=0; d<numSpecialDOFs; d++)
            specialConstr.set(spesh, d, coeffs.getEntry(specialDOFList.get(d)));
        target.set( spesh, 0, 0.5*(voxLinConstr.get(2*spesh).getValue()+voxLinConstr.get(2*spesh+1).getValue()) );
    }
    
    //remaining targets (orth components to voxLinConstr) are 0
    DoubleMatrix2D orthComponents = getOrthogVectors(specialConstr);
    DoubleMatrix2D M = DoubleFactory2D.dense.compose(
            new DoubleMatrix2D[][] {
                new DoubleMatrix2D[] {specialConstr},
                new DoubleMatrix2D[] {Algebra.DEFAULT.transpose(orthComponents)}
            }
        );
    
    DoubleMatrix1D specialDOFVals = Algebra.DEFAULT.solve(M, target).viewColumn(0);
    
    for(int d=0; d<numSpecialDOFs; d++){
        ans.put(dofIntervals.get(specialDOFList.get(d)).dof, specialDOFVals.get(d));
    }
    
    return ans;
}
 

/**
 * Calculates based on windows of magnitude spectrum. Encompasses portion of
 * Moments class, but has a delay of lengthOfWindow windows before any
 * results are calculated.
 * 
 * @param samples
 *            The samples to extract the feature from.
 * @param sampling_rate
 *            The sampling rate that the samples are encoded with.
 * @param other_feature_values
 *            The values of other features that are needed to calculate this
 *            value. The order and offsets of these features must be the
 *            same as those returned by this class's getDependencies and
 *            getDependencyOffsets methods respectively. The first indice
 *            indicates the feature/window and the second indicates the
 *            value.
 * @return The extracted feature value(s).
 * @throws Exception
 *             Throws an informative exception if the feature cannot be
 *             calculated.
 */
public double[] extractFeature(double[] samples, double sampling_rate,
		double[][] other_feature_values) throws Exception {
	if((featureLength != other_feature_values[0].length)||(windowLength != other_feature_values.length)){
		terms = new DenseDoubleMatrix2D(k*l,windowLength*featureLength);
		z = new DenseDoubleMatrix2D(1,featureLength*windowLength);
		calcTerms(terms);
	}
	for(int i=0;i<windowLength;++i){
		for(int j=0;j<featureLength;++j){
			z.set(0,featureLength*i+j,other_feature_values[i][j]);
		}
	}
	DoubleMatrix2D retMatrix = (new Algebra()).solve(terms,z);
	return retMatrix.viewRow(0).toArray();
}
 

/**
 * Calculates based on windows of magnitude spectrum. Encompasses portion of
 * Moments class, but has a delay of lengthOfWindow windows before any
 * results are calculated.
 * 
 * @param samples
 *            The samples to extract the feature from.
 * @param sampling_rate
 *            The sampling rate that the samples are encoded with.
 * @param other_feature_values
 *            The values of other features that are needed to calculate this
 *            value. The order and offsets of these features must be the
 *            same as those returned by this class's getDependencies and
 *            getDependencyOffsets methods respectively. The first indice
 *            indicates the feature/window and the second indicates the
 *            value.
 * @return The extracted feature value(s).
 * @throws Exception
 *             Throws an informative exception if the feature cannot be
 *             calculated.
 */
public double[] extractFeature(double[] samples, double sampling_rate,
		double[][] other_feature_values) throws Exception {
	if((featureLength != other_feature_values[0].length)||(windowLength != other_feature_values.length)){
		terms = new DenseDoubleMatrix2D(k*l,windowLength*featureLength);
		z = new DenseDoubleMatrix2D(1,featureLength*windowLength);
		calcTerms(terms);
	}
	for(int i=0;i<windowLength;++i){
		for(int j=0;j<featureLength;++j){
			z.set(0,featureLength*i+j,other_feature_values[i][j]);
		}
	}
	DoubleMatrix2D retMatrix = (new Algebra()).solve(terms,z);
	return retMatrix.viewRow(0).toArray();
}
 

/**
 * Calculates based on windows of magnitude spectrum. Encompasses portion of
 * Moments class, but has a delay of lengthOfWindow windows before any
 * results are calculated.
 * 
 * @param samples
 *            The samples to extract the feature from.
 * @param sampling_rate
 *            The sampling rate that the samples are encoded with.
 * @param other_feature_values
 *            The values of other features that are needed to calculate this
 *            value. The order and offsets of these features must be the
 *            same as those returned by this class's getDependencies and
 *            getDependencyOffsets methods respectively. The first indice
 *            indicates the feature/window and the second indicates the
 *            value.
 * @return The extracted feature value(s).
 * @throws Exception
 *             Throws an informative exception if the feature cannot be
 *             calculated.
 */
public double[] extractFeature(double[] samples, double sampling_rate,
		double[][] other_feature_values) throws Exception {
	if((featureLength != other_feature_values[0].length)||(windowLength != other_feature_values.length)){
		terms = new DenseDoubleMatrix2D(windowLength*featureLength,k*l);
		z = new DenseDoubleMatrix2D(1,featureLength*windowLength);
		calcTerms(terms);
	}
	for(int i=0;i<windowLength;++i){
		for(int j=0;j<featureLength;++j){
			z.set(featureLength*i+j,0,other_feature_values[i][j]);
		}
	}
	DoubleMatrix2D retMatrix = (new Algebra()).solve(terms,z);
	return retMatrix.viewRow(0).toArray();
}
 

public static void main(String args[]) {

       // For COLT
       DoubleMatrix2D xmatrix,ymatrix,zmatrix;

       DoubleFactory2D myfactory;
       myfactory = DoubleFactory2D.dense;
       xmatrix = myfactory.make(8,2);
       ymatrix = myfactory.make(8,1);

       xmatrix.set(0,0,1);
       xmatrix.set(1,0,1);
       xmatrix.set(2,0,1);
       xmatrix.set(3,0,1);
       xmatrix.set(4,0,1);
       xmatrix.set(5,0,1); 
       xmatrix.set(6,0,1);
       xmatrix.set(7,0,1);

       xmatrix.set(0,1,80);
       xmatrix.set(1,1,220);
       xmatrix.set(2,1,140);
       xmatrix.set(3,1,120);
       xmatrix.set(4,1,180);
       xmatrix.set(5,1,100);
       xmatrix.set(6,1,200);
       xmatrix.set(7,1,160);

       ymatrix.set(0,0,0.6);
       ymatrix.set(1,0,6.70);
       ymatrix.set(2,0,5.30);
       ymatrix.set(3,0,4.00);
       ymatrix.set(4,0,6.55);
       ymatrix.set(5,0,2.15);
       ymatrix.set(6,0,6.60);
       ymatrix.set(7,0,5.75);

       Algebra myAlgebra = new Algebra();
       zmatrix = myAlgebra.solve(xmatrix,ymatrix);
       System.err.println(xmatrix);
       System.err.println(ymatrix);
       System.err.println(zmatrix);

       /*
       // For JAMA
       Matrix amatrix,bmatrix,cmatrix;
       amatrix = new Matrix(8,2);
       bmatrix = new Matrix(8,1);

       amatrix.set(0,0,1);
       amatrix.set(1,0,1);
       amatrix.set(2,0,1);
       amatrix.set(3,0,1);
       amatrix.set(4,0,1);
       amatrix.set(5,0,1);
       amatrix.set(6,0,1);
       amatrix.set(7,0,1);

       amatrix.set(0,1,80);
       amatrix.set(1,1,220);
       amatrix.set(2,1,140);
       amatrix.set(3,1,120);
       amatrix.set(4,1,180);
       amatrix.set(5,1,100);
       amatrix.set(6,1,200);
       amatrix.set(7,1,160);

       bmatrix.set(0,0,0.6);
       bmatrix.set(1,0,6.70);
       bmatrix.set(2,0,5.30);
       bmatrix.set(3,0,4.00);
       bmatrix.set(4,0,6.55);
       bmatrix.set(5,0,2.15);
       bmatrix.set(6,0,6.60);
       bmatrix.set(7,0,5.75);

       cmatrix = amatrix.solve(bmatrix);
       amatrix.print(8,5);
       bmatrix.print(8,5);
       cmatrix.print(8,5);
		*/
    }
 

/**
 */
public static void testLU() {
double[][] vals = {
	{-0.074683,  0.321248,-0.014656, 0.286586,0},
	{-0.344852, -0.16278 , 0.173711, 0.00064 ,0},
	{-0.181924, -0.092926, 0.184153, 0.177966,1},
	{-0.166829, -0.10321 , 0.582301, 0.142583,0},
	{ 0       , -0.112952,-0.04932 ,-0.700157,0},
	{ 0       , 0        ,0        ,0        ,0}
};
 
DoubleMatrix2D H = new DenseDoubleMatrix2D( vals ); // see values below...
System.out.println("\nHplus="+H.viewDice().zMult(H,null));

DoubleMatrix2D Hplus = Algebra.DEFAULT.inverse(H.viewDice().zMult( H,null )).zMult(H.viewDice(),null);
Hplus.assign(cern.jet.math.Functions.round(1.0E-10));
System.out.println("\nHplus="+Hplus);

		/*
DoubleMatrix2D HtH = new DenseDoubleMatrix2D( 5, 5 );
DoubleMatrix2D Hplus = new DenseDoubleMatrix2D( 5, 6 );
LUDecompositionQuick LUD = new LUDecompositionQuick();
		//H.zMult( H, HtH, 1, 0, true, false );
		//DoubleMatrix2D res = Algebra.DEFAULT.inverse(HtH).zMult(H,null,1,0,false,true);
		LUD.decompose( HtH );
		// first fill Hplus with the transpose of H...
		for (int i = 0; i < 6; i++ ) {
			for ( int j = 0; j < 5; j++ ) {
				Hplus.set( j, i, H.get( i, j ) );
			}
		}
		LUD.solve( Hplus );

		DoubleMatrix2D perm = Algebra.DEFAULT.permute(Hplus, null,LUD.getPivot());
		DoubleMatrix2D inv = Algebra.DEFAULT.inverse(HtH);//.zMult(H,null,1,0,false,true);
		*/

		// in matlab...
		// Hplus = inv(H' * H) * H'

//System.out.println("\nLU="+LUD);
//System.out.println("\nHplus="+Hplus);
//System.out.println("\nperm="+perm);
//System.out.println("\ninv="+inv);
//System.out.println("\nres="+res);
}
 

/**
 */
public static void testLU() {
double[][] vals = {
	{-0.074683,  0.321248,-0.014656, 0.286586,0},
	{-0.344852, -0.16278 , 0.173711, 0.00064 ,0},
	{-0.181924, -0.092926, 0.184153, 0.177966,1},
	{-0.166829, -0.10321 , 0.582301, 0.142583,0},
	{ 0       , -0.112952,-0.04932 ,-0.700157,0},
	{ 0       , 0        ,0        ,0        ,0}
};
 
DoubleMatrix2D H = new DenseDoubleMatrix2D( vals ); // see values below...
System.out.println("\nHplus="+H.viewDice().zMult(H,null));

DoubleMatrix2D Hplus = Algebra.DEFAULT.inverse(H.viewDice().zMult( H,null )).zMult(H.viewDice(),null);
Hplus.assign(cern.jet.math.Functions.round(1.0E-10));
System.out.println("\nHplus="+Hplus);

		/*
DoubleMatrix2D HtH = new DenseDoubleMatrix2D( 5, 5 );
DoubleMatrix2D Hplus = new DenseDoubleMatrix2D( 5, 6 );
LUDecompositionQuick LUD = new LUDecompositionQuick();
		//H.zMult( H, HtH, 1, 0, true, false );
		//DoubleMatrix2D res = Algebra.DEFAULT.inverse(HtH).zMult(H,null,1,0,false,true);
		LUD.decompose( HtH );
		// first fill Hplus with the transpose of H...
		for (int i = 0; i < 6; i++ ) {
			for ( int j = 0; j < 5; j++ ) {
				Hplus.set( j, i, H.get( i, j ) );
			}
		}
		LUD.solve( Hplus );

		DoubleMatrix2D perm = Algebra.DEFAULT.permute(Hplus, null,LUD.getPivot());
		DoubleMatrix2D inv = Algebra.DEFAULT.inverse(HtH);//.zMult(H,null,1,0,false,true);
		*/

		// in matlab...
		// Hplus = inv(H' * H) * H'

//System.out.println("\nLU="+LUD);
//System.out.println("\nHplus="+Hplus);
//System.out.println("\nperm="+perm);
//System.out.println("\ninv="+inv);
//System.out.println("\nres="+res);
}
 

public static void main(String args[]) {

       // For COLT
       DoubleMatrix2D xmatrix,ymatrix,zmatrix;

       DoubleFactory2D myfactory;
       myfactory = DoubleFactory2D.dense;
       xmatrix = myfactory.make(8,2);
       ymatrix = myfactory.make(8,1);

       xmatrix.set(0,0,1);
       xmatrix.set(1,0,1);
       xmatrix.set(2,0,1);
       xmatrix.set(3,0,1);
       xmatrix.set(4,0,1);
       xmatrix.set(5,0,1); 
       xmatrix.set(6,0,1);
       xmatrix.set(7,0,1);

       xmatrix.set(0,1,80);
       xmatrix.set(1,1,220);
       xmatrix.set(2,1,140);
       xmatrix.set(3,1,120);
       xmatrix.set(4,1,180);
       xmatrix.set(5,1,100);
       xmatrix.set(6,1,200);
       xmatrix.set(7,1,160);

       ymatrix.set(0,0,0.6);
       ymatrix.set(1,0,6.70);
       ymatrix.set(2,0,5.30);
       ymatrix.set(3,0,4.00);
       ymatrix.set(4,0,6.55);
       ymatrix.set(5,0,2.15);
       ymatrix.set(6,0,6.60);
       ymatrix.set(7,0,5.75);

       Algebra myAlgebra = new Algebra();
       zmatrix = myAlgebra.solve(xmatrix,ymatrix);
       System.err.println(xmatrix);
       System.err.println(ymatrix);
       System.err.println(zmatrix);

       /*
       // For JAMA
       Matrix amatrix,bmatrix,cmatrix;
       amatrix = new Matrix(8,2);
       bmatrix = new Matrix(8,1);

       amatrix.set(0,0,1);
       amatrix.set(1,0,1);
       amatrix.set(2,0,1);
       amatrix.set(3,0,1);
       amatrix.set(4,0,1);
       amatrix.set(5,0,1);
       amatrix.set(6,0,1);
       amatrix.set(7,0,1);

       amatrix.set(0,1,80);
       amatrix.set(1,1,220);
       amatrix.set(2,1,140);
       amatrix.set(3,1,120);
       amatrix.set(4,1,180);
       amatrix.set(5,1,100);
       amatrix.set(6,1,200);
       amatrix.set(7,1,160);

       bmatrix.set(0,0,0.6);
       bmatrix.set(1,0,6.70);
       bmatrix.set(2,0,5.30);
       bmatrix.set(3,0,4.00);
       bmatrix.set(4,0,6.55);
       bmatrix.set(5,0,2.15);
       bmatrix.set(6,0,6.60);
       bmatrix.set(7,0,5.75);

       cmatrix = amatrix.solve(bmatrix);
       amatrix.print(8,5);
       bmatrix.print(8,5);
       cmatrix.print(8,5);
		*/
    }