List of usage examples for org.apache.commons.math3.linear EigenDecomposition getRealEigenvalues
public double[] getRealEigenvalues()
From source file:com.datumbox.framework.machinelearning.featureselection.continuous.PCA.java
@Override protected void estimateModelParameters(Dataset originaldata) { ModelParameters modelParameters = knowledgeBase.getModelParameters(); int n = originaldata.size(); int d = originaldata.getColumnSize(); //convert data into matrix Map<Object, Integer> feature2ColumnId = modelParameters.getFeature2ColumnId(); MatrixDataset matrixDataset = MatrixDataset.newInstance(originaldata, false, feature2ColumnId); RealMatrix X = matrixDataset.getX(); //calculate means and subtract them from data double[] meanValues = new double[d]; for (Map.Entry<Object, Integer> entry : feature2ColumnId.entrySet()) { Object feature = entry.getKey(); Integer columnId = entry.getValue(); meanValues[columnId] = Descriptives .mean(originaldata.extractColumnValues(feature).toFlatDataCollection()); for (int row = 0; row < n; ++row) { X.addToEntry(row, columnId, -meanValues[columnId]); //inplace subtraction!!! }/*from ww w .j a v a 2 s . com*/ } modelParameters.setMean(meanValues); RealMatrix components; double[] eigenValues; /* if(d>n) { //turned off because of the algorithm could not be validated //Karhunen Lowe Transform to speed up calculations //nxn matrix RealMatrix covarianceNN = (X.multiply(X.transpose())).scalarMultiply(1.0/(n-1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceNN); eigenValues = decomposition.getRealEigenvalues(); RealMatrix eigenVectors = decomposition.getV(); double[] sqrtInverseEigenValues = new double[eigenValues.length]; for(int i=0;i<eigenValues.length;++i) { if(eigenValues[i]==0.0) { sqrtInverseEigenValues[i] = 0.0; } else { sqrtInverseEigenValues[i] = 1.0/Math.sqrt(eigenValues[i]); } } components = X.transpose().multiply(eigenVectors); //Components = X'*V*L^-0.5; To whiten them we multiply with L^0.5 which //cancels out the previous multiplication. So below we multiply by //L^-0.5 ONLY if we don't whiten. if(!knowledgeBase.getTrainingParameters().isWhitened()) { components = components.multiply(new DiagonalMatrix(sqrtInverseEigenValues)); } } else { //Normal PCA goes here } */ //dxd matrix RealMatrix covarianceDD = (X.transpose().multiply(X)).scalarMultiply(1.0 / (n - 1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceDD); eigenValues = decomposition.getRealEigenvalues(); components = decomposition.getV(); //Whiten Components W = U*L^0.5; To whiten them we multiply with L^0.5. if (knowledgeBase.getTrainingParameters().isWhitened()) { double[] sqrtEigenValues = new double[eigenValues.length]; for (int i = 0; i < eigenValues.length; ++i) { sqrtEigenValues[i] = Math.sqrt(eigenValues[i]); } components = components.multiply(new DiagonalMatrix(sqrtEigenValues)); } //the eigenvalues and their components are sorted by descending order no need to resort them Integer maxDimensions = knowledgeBase.getTrainingParameters().getMaxDimensions(); Double varianceThreshold = knowledgeBase.getTrainingParameters().getVarianceThreshold(); if (varianceThreshold != null && varianceThreshold <= 1) { double sum = 0.0; double totalVariance = StatUtils.sum(eigenValues); int varCounter = 0; for (double l : eigenValues) { sum += l / totalVariance; ++varCounter; if (sum >= varianceThreshold) { break; } } if (maxDimensions == null || maxDimensions > varCounter) { maxDimensions = varCounter; } } if (maxDimensions != null && maxDimensions < d) { //keep only the maximum selected eigenvalues double[] newEigenValues = new double[maxDimensions]; System.arraycopy(eigenValues, 0, newEigenValues, 0, maxDimensions); eigenValues = newEigenValues; //keep only the maximum selected eigenvectors components = components.getSubMatrix(0, components.getRowDimension() - 1, 0, maxDimensions - 1); } modelParameters.setRows(components.getRowDimension()); modelParameters.setCols(components.getColumnDimension()); modelParameters.setEigenValues(eigenValues); modelParameters.setComponents(components.getData()); }
From source file:lanchester.MultiArena2.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2);/*from w w w.ja v a 2 s.c o m*/ } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); // for(int i1=0;i1<eVals.length;i1++){ // System.out.println("eVals["+i1+"] = "+eVals[i1]); // } if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { AthenaForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getForceSize() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; int cntr = 0; int numGone; do { MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; if (forces.get(i1).getForceSize() > lb) { for (int i2 = 0; i2 < numFoes; i2++) { // updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * timeStep); updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); if (updatedForce < 1.) { updatedForce = 0.; numGone++; } // updatedForce+=coeffs[i2]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i2]*timeStep); // updatedForce+=coeffs[i1]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i1]*timeStep); } } else { updatedForce = lb / 2.; numGone++; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); } history.add(currentStep); aboveFloor = checkAboveFloors(); cntr++; } while (aboveFloor && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2); } } } // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:com.datumbox.framework.core.machinelearning.featureselection.PCA.java
/** {@inheritDoc} */ @Override/*from w ww.j a va2s .c o m*/ protected void _fit(Dataframe trainingData) { ModelParameters modelParameters = knowledgeBase.getModelParameters(); int n = trainingData.size(); int d = trainingData.xColumnSize(); //convert data into matrix Map<Object, Integer> featureIds = modelParameters.getFeatureIds(); DataframeMatrix matrixDataset = DataframeMatrix.newInstance(trainingData, false, null, featureIds); RealMatrix X = matrixDataset.getX(); //calculate means and subtract them from data RealVector meanValues = new OpenMapRealVector(d); for (Integer columnId : featureIds.values()) { double mean = 0.0; for (int row = 0; row < n; row++) { mean += X.getEntry(row, columnId); } mean /= n; for (int row = 0; row < n; row++) { X.addToEntry(row, columnId, -mean); } meanValues.setEntry(columnId, mean); } modelParameters.setMean(meanValues); //dxd matrix RealMatrix covarianceDD = (X.transpose().multiply(X)).scalarMultiply(1.0 / (n - 1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceDD); RealVector eigenValues = new ArrayRealVector(decomposition.getRealEigenvalues(), false); RealMatrix components = decomposition.getV(); //Whiten Components W = U*L^0.5; To whiten them we multiply with L^0.5. if (knowledgeBase.getTrainingParameters().isWhitened()) { RealMatrix sqrtEigenValues = new DiagonalMatrix(d); for (int i = 0; i < d; i++) { sqrtEigenValues.setEntry(i, i, FastMath.sqrt(eigenValues.getEntry(i))); } components = components.multiply(sqrtEigenValues); } //the eigenvalues and their components are sorted by descending order no need to resort them Integer maxDimensions = knowledgeBase.getTrainingParameters().getMaxDimensions(); Double variancePercentageThreshold = knowledgeBase.getTrainingParameters().getVariancePercentageThreshold(); if (variancePercentageThreshold != null && variancePercentageThreshold <= 1) { double totalVariance = 0.0; for (int i = 0; i < d; i++) { totalVariance += eigenValues.getEntry(i); } double sum = 0.0; int varCounter = 0; for (int i = 0; i < d; i++) { sum += eigenValues.getEntry(i) / totalVariance; varCounter++; if (sum >= variancePercentageThreshold) { break; } } if (maxDimensions == null || maxDimensions > varCounter) { maxDimensions = varCounter; } } if (maxDimensions != null && maxDimensions < d) { //keep only the maximum selected eigenvalues eigenValues = eigenValues.getSubVector(0, maxDimensions); //keep only the maximum selected eigenvectors components = components.getSubMatrix(0, components.getRowDimension() - 1, 0, maxDimensions - 1); } modelParameters.setEigenValues(eigenValues); modelParameters.setComponents(components); }
From source file:com.datumbox.framework.core.machinelearning.featureselection.continuous.PCA.java
/** {@inheritDoc} */ @Override/* w w w. ja va2 s . co m*/ protected void _fit(Dataframe originalData) { ModelParameters modelParameters = kb().getModelParameters(); int n = modelParameters.getN(); int d = modelParameters.getD(); //convert data into matrix Map<Object, Integer> featureIds = modelParameters.getFeatureIds(); MatrixDataframe matrixDataset = MatrixDataframe.newInstance(originalData, false, null, featureIds); RealMatrix X = matrixDataset.getX(); //calculate means and subtract them from data double[] meanValues = new double[d]; for (Integer columnId : featureIds.values()) { meanValues[columnId] = 0.0; for (double v : X.getColumn(columnId)) { meanValues[columnId] += v; } meanValues[columnId] /= n; for (int row = 0; row < n; row++) { X.addToEntry(row, columnId, -meanValues[columnId]); } } modelParameters.setMean(meanValues); RealMatrix components; double[] eigenValues; //dxd matrix RealMatrix covarianceDD = (X.transpose().multiply(X)).scalarMultiply(1.0 / (n - 1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceDD); eigenValues = decomposition.getRealEigenvalues(); components = decomposition.getV(); //Whiten Components W = U*L^0.5; To whiten them we multiply with L^0.5. if (kb().getTrainingParameters().isWhitened()) { double[] sqrtEigenValues = new double[eigenValues.length]; for (int i = 0; i < eigenValues.length; i++) { sqrtEigenValues[i] = Math.sqrt(eigenValues[i]); } components = components.multiply(new DiagonalMatrix(sqrtEigenValues)); } //the eigenvalues and their components are sorted by descending order no need to resort them Integer maxDimensions = kb().getTrainingParameters().getMaxDimensions(); Double variancePercentageThreshold = kb().getTrainingParameters().getVariancePercentageThreshold(); if (variancePercentageThreshold != null && variancePercentageThreshold <= 1) { double sum = 0.0; double totalVariance = StatUtils.sum(eigenValues); int varCounter = 0; for (double l : eigenValues) { sum += l / totalVariance; varCounter++; if (sum >= variancePercentageThreshold) { break; } } if (maxDimensions == null || maxDimensions > varCounter) { maxDimensions = varCounter; } } if (maxDimensions != null && maxDimensions < d) { //keep only the maximum selected eigenvalues double[] newEigenValues = new double[maxDimensions]; System.arraycopy(eigenValues, 0, newEigenValues, 0, maxDimensions); eigenValues = newEigenValues; //keep only the maximum selected eigenvectors components = components.getSubMatrix(0, components.getRowDimension() - 1, 0, maxDimensions - 1); } modelParameters.setRows(components.getRowDimension()); modelParameters.setCols(components.getColumnDimension()); modelParameters.setEigenValues(eigenValues); modelParameters.setComponents(components.getData()); }
From source file:MultivariateNormalDistribution.java
/** * Creates a multivariate normal distribution with the given mean vector and * covariance matrix.//from w w w . j av a 2s. co m * <br/> * The number of dimensions is equal to the length of the mean vector * and to the number of rows and columns of the covariance matrix. * It is frequently written as "p" in formulae. * * @param rng Random Number Generator. * @param means Vector of means. * @param covariances Covariance matrix. * @throws DimensionMismatchException if the arrays length are * inconsistent. * @throws SingularMatrixException if the eigenvalue decomposition cannot * be performed on the provided covariance matrix. * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is * negative. */ public MultivariateNormalDistribution(RandomGenerator rng, final double[] means, final double[][] covariances) throws SingularMatrixException, DimensionMismatchException, NonPositiveDefiniteMatrixException { super(rng, means.length); final int dim = means.length; if (covariances.length != dim) { throw new DimensionMismatchException(covariances.length, dim); } for (int i = 0; i < dim; i++) { if (dim != covariances[i].length) { throw new DimensionMismatchException(covariances[i].length, dim); } } this.means = MathArrays.copyOf(means); covarianceMatrix = new Array2DRowRealMatrix(covariances); // Covariance matrix eigen decomposition. final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix); // Compute and store the inverse. covarianceMatrixInverse = covMatDec.getSolver().getInverse(); // Compute and store the determinant. covarianceMatrixDeterminant = covMatDec.getDeterminant(); // Eigenvalues of the covariance matrix. final double[] covMatEigenvalues = covMatDec.getRealEigenvalues(); for (int i = 0; i < covMatEigenvalues.length; i++) { if (covMatEigenvalues[i] < 0) { throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0); } } // Matrix where each column is an eigenvector of the covariance matrix. final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim); for (int v = 0; v < dim; v++) { final double[] evec = covMatDec.getEigenvector(v).toArray(); covMatEigenvectors.setColumn(v, evec); } final RealMatrix tmpMatrix = covMatEigenvectors.transpose(); // Scale each eigenvector by the square root of its eigenvalue. for (int row = 0; row < dim; row++) { final double factor = FastMath.sqrt(covMatEigenvalues[row]); for (int col = 0; col < dim; col++) { tmpMatrix.multiplyEntry(row, col, factor); } } samplingMatrix = covMatEigenvectors.multiply(tmpMatrix); }
From source file:lanchester.MultiArena3.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; currentCeiling = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; currentCeiling[i1][i2] = 100.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2);// w ww.j a v a2s . c om setCeiling(i1, i2); } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { MultiForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getNumber() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; boolean belowCeiling = true; int cntr = 0; int numGone; do { timeStep = determineTimeStep(); MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; if (stillAlive[i1]) { for (int i2 = 0; i2 < numFoes; i2++) { updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); } if (updatedForce < 1.) { updatedForce = lb; stillAlive[i1] = false; numGone++; } } else { numGone++; updatedForce = lb; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); } history.add(currentStep); aboveFloor = checkAboveFloors(); belowCeiling = checkBelowCeilings(); cntr++; } while (aboveFloor && belowCeiling && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setFloor(i1, i2); } } } // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:lanchester.MultiArena4.java
public void step() { boolean aboveFloor = true; double currentCycle = 0.; int numFoes = forces.size(); if (isMyTurn == null) { isMyTurn = new boolean[numFoes][numFoes]; stanceArray = new int[numFoes][numFoes]; currentFloor = new double[numFoes][numFoes]; currentCeiling = new double[numFoes][numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { // int ind1 = forceMap.get(forces.get(i1)); for (int i2 = 0; i2 < numFoes; i2++) { // int ind2 = forceMap.get(forces.get(i2)); isMyTurn[i1][i2] = true; if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; currentCeiling[i1][i2] = 1000000.; } else { stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); // setStances(i1,i2); setFloor(i1, i2);// ww w .j ava 2 s .c om setCeiling(i1, i2); } } } } Array2DRowRealMatrix mat = getMat(); EigenDecomposition eigen = new EigenDecomposition(mat); double det = eigen.getDeterminant(); double[] eVals = eigen.getRealEigenvalues(); if (eigen.hasComplexEigenvalues()) { System.out.println("Complex eigenvalues"); for (int i1 = 0; i1 < forces.size(); i1++) { MultiForce f = forces.get(i1); System.out.println(f.getName() + " has " + f.getNumber() + " forces remaining"); } } double[] initialNums = getInitialNumbers(forces); Array2DRowRealMatrix eVectors = (Array2DRowRealMatrix) eigen.getV(); LUDecomposition lu = new LUDecomposition(eVectors); double det2 = lu.getDeterminant(); double[] coeffs = new double[numFoes]; for (int i1 = 0; i1 < numFoes; i1++) { Array2DRowRealMatrix tmpMat = (Array2DRowRealMatrix) eVectors.copy(); tmpMat.setColumn(i1, initialNums); LUDecomposition tmpLU = new LUDecomposition(tmpMat); double tmpDet = tmpLU.getDeterminant(); coeffs[i1] = tmpDet / det2; } aboveFloor = true; boolean belowCeiling = true; int cntr = 0; int numGone; do { timeStep = determineTimeStep(); MultiTimeStep currentStep = new MultiTimeStep(numFoes); currentTime += timeStep; currentCycle += timeStep; currentStep.setTime(currentTime); numGone = 0; for (int i1 = 0; i1 < numFoes; i1++) { double updatedForce = 0.; // if (forces.get(i1).getForceSize() > lb) { if (stillAlive[i1]) { for (int i2 = 0; i2 < numFoes; i2++) { // updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * timeStep); updatedForce += coeffs[i2] * eVectors.getEntry(i1, i2) * Math.exp(eVals[i2] * currentCycle); // updatedForce+=coeffs[i2]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i2]*timeStep); // updatedForce+=coeffs[i1]*eVectors.getEntry(i2, i1)*Math.exp(eVals[i1]*timeStep); } if (updatedForce < 1.) { updatedForce = lb; stillAlive[i1] = false; numGone++; } } else { updatedForce = lb; numGone++; } forces.get(i1).updateForce(updatedForce); currentStep.setForceNumber(updatedForce, i1); // for (int i2 = 0; i2 < numFoes; i1++) { // if (i1 != i2) { // // } // } } history.add(currentStep); aboveFloor = checkAboveFloors(); belowCeiling = checkBelowCeilings(); cntr++; } while (aboveFloor && belowCeiling && cntr < 2000 && (numFoes - numGone) > 1); for (int i1 = 0; i1 < numFoes; i1++) { for (int i2 = 0; i2 < numFoes; i2++) { if (i1 == i2) { stanceArray[i1][i2] = AthenaConstants.ALLIED_POSTURE; currentFloor[i1][i2] = 0.; } else { // stanceArray[i1][i2] = initializeStance(forces.get(i1), forces.get(i2)); setStances(i1, i2); setFloor(i1, i2); setCeiling(i1, i2); } } } // eVectors. // this.currentTime++; // Truncator truncator = new Truncator(); if (numFoes - numGone == 1) { loneSurvivor = true; // System.out.println("time = " + time); } }
From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * Pre-calculate intermediate matrices, perform setup, etc. */// w w w. j av a2 s. c om private void initialize(Test test, RealMatrix FEssence, RealMatrix FtFinverse, int perGroupN, FixedRandomMatrix CFixedRand, RealMatrix U, RealMatrix thetaNull, RealMatrix beta, RealMatrix sigmaError, RealMatrix sigmaG, boolean exact) throws PowerException { debug("entering initialize"); // reset member variables this.T1 = null; this.FT1 = null; this.S = null; this.mzSq = null; this.H0 = 0; this.sStar = 0; // cache inputs this.test = test; this.FEssence = FEssence; this.FtFinverse = FtFinverse; this.perGroupN = perGroupN; this.CFixedRand = CFixedRand; this.U = U; this.thetaNull = thetaNull; this.beta = beta; this.sigmaError = sigmaError; this.sigmaG = sigmaG; // calculate intermediate matrices // RealMatrix FEssence = params.getDesignEssence().getFullDesignMatrixFixed(); // TODO: do we ever get here with values that can cause integer overflow, // and if so, does it matter? this.N = (double) FEssence.getRowDimension() * perGroupN; this.exact = exact; try { // TODO: need to calculate H0, need to adjust H1 for Unirep // get design matrix for fixed parameters only qF = FEssence.getColumnDimension(); // a = CFixedRand.getCombinedMatrix().getRowDimension(); // get fixed contrasts RealMatrix Cfixed = CFixedRand.getFixedMatrix(); RealMatrix CGaussian = CFixedRand.getRandomMatrix(); // build intermediate terms h1, S if (FtFinverse == null) { FtFinverse = new LUDecomposition(FEssence.transpose().multiply(FEssence)).getSolver().getInverse(); debug("FEssence", FEssence); debug("FtFinverse = (FEssence transpose * FEssence) inverse", FtFinverse); } else { debug("FtFinverse", FtFinverse); } RealMatrix PPt = Cfixed.multiply(FtFinverse.scalarMultiply(1 / (double) perGroupN)) .multiply(Cfixed.transpose()); debug("Cfixed", Cfixed); debug("n = " + perGroupN); debug("PPt = Cfixed * FtF inverse * (1/n) * Cfixed transpose", PPt); T1 = forceSymmetric(new LUDecomposition(PPt).getSolver().getInverse()); debug("T1 = PPt inverse", T1); FT1 = new CholeskyDecomposition(T1).getL(); debug("FT1 = Cholesky decomposition (L) of T1", FT1); // calculate theta difference // RealMatrix thetaNull = params.getTheta(); RealMatrix C = CFixedRand.getCombinedMatrix(); // RealMatrix beta = params.getScaledBeta(); // RealMatrix U = params.getWithinSubjectContrast(); // thetaHat = C * beta * U RealMatrix thetaHat = C.multiply(beta.multiply(U)); debug("C", C); debug("beta", beta); debug("U", U); debug("thetaHat = C * beta * U", thetaHat); // thetaDiff = thetaHat - thetaNull RealMatrix thetaDiff = thetaHat.subtract(thetaNull); debug("thetaNull", thetaNull); debug("thetaDiff = thetaHat - thetaNull", thetaDiff); // TODO: specific to HLT or UNIREP RealMatrix sigmaStarInverse = getSigmaStarInverse(U, sigmaError, test); debug("sigmaStarInverse", sigmaStarInverse); RealMatrix H1matrix = thetaDiff.transpose().multiply(T1).multiply(thetaDiff).multiply(sigmaStarInverse); debug("H1matrix = thetaDiff transpose * T1 * thetaDiff * sigmaStarInverse", H1matrix); H1 = H1matrix.getTrace(); debug("H1 = " + H1); // Matrix which represents the non-centrality parameter as a linear combination of chi-squared r.v.'s. S = FT1.transpose().multiply(thetaDiff).multiply(sigmaStarInverse).multiply(thetaDiff.transpose()) .multiply(FT1).scalarMultiply(1 / H1); debug("S = FT1 transpose * thetaDiff * sigmaStar inverse * thetaDiff transpose * FT1 * (1/H1)", S); // We use the S matrix to generate the F-critical, numerical df's, and denominator df's // for a central F distribution. The resulting F distribution is used as an approximation // for the distribution of the non-centrality parameter. // See formulas 18-21 and A8,A10 from Glueck & Muller (2003) for details. EigenDecomposition sEigenDecomp = new EigenDecomposition(S); sEigenValues = sEigenDecomp.getRealEigenvalues(); // calculate H0 if (sEigenValues.length > 0) H0 = H1 * (1 - sEigenValues[0]); if (H0 <= 0) H0 = 0; // count the # of positive eigen values for (double value : sEigenValues) { if (value > 0) sStar++; } // TODO: throw error if sStar is <= 0 // TODO: NO: throw error if sStar != sEigenValues.length instead??? double stddevG = Math.sqrt(sigmaG.getEntry(0, 0)); RealMatrix svec = sEigenDecomp.getVT(); mzSq = svec.multiply(FT1.transpose()).multiply(CGaussian).scalarMultiply(1 / stddevG); for (int i = 0; i < mzSq.getRowDimension(); i++) { for (int j = 0; j < mzSq.getColumnDimension(); j++) { double entry = mzSq.getEntry(i, j); mzSq.setEntry(i, j, entry * entry); // TODO: is there an apache function to do this? } } debug("exiting initialize normally"); } catch (RuntimeException e) { LOGGER.warn("exiting initialize abnormally", e); throw new PowerException(e.getMessage(), PowerErrorEnum.INVALID_DISTRIBUTION_NONCENTRALITY_PARAMETER); } }
From source file:com.trickl.math.lanczos.LanczosSolver.java
public Pair<DoubleMatrix2D, Info> getEigenvectorsWithInfo(double[] eigenvalues, RandomGenerator randomGenerator, int maxIterations) { List<DoubleMatrix1D> eigvectors = new LinkedList<>(); // contains ritz vectors. List<List<Double>> Tvectors = new LinkedList<>(); // contains // calculation of eigen vectors of T matrix(consists of alphas & betas): int start = 0; int end = eigenvalues.length; int n1 = alpha.length; double mamax, error, lambda; Pair<Double, Double> a_and_b; int ma = 0, deltam; int nth, maMax = 0, count; List<Double> eigenval_a = new LinkedList<>(); List<Double> residuum = new LinkedList<>(); List<ErrorInfo> status = new LinkedList<>(); findM1M2(eigenvalues);/*from www .j a va 2 s .c o m*/ int M1_itr = 0; int M2_itr = 0; for (int pos = start; pos < end; ++pos) { int maxcount = 10; lambda = 0; count = 0; ErrorInfo errInf = ErrorInfo.OK; // calculation of ma starts: if (M1[M1_itr] != 0 && M2[M2_itr] != 0) { ma = (3 * (M1[M1_itr]) + M2[M2_itr]) / 4 + 1; deltam = ((3 * (M1[M1_itr]) + 5 * (M2[M2_itr])) / 8 + 1 - ma) / 10 + 1; } else if (M1[M1_itr] != 0 && M2[M2_itr] == 0) { ma = (5 * (M1[M1_itr])) / 4 + 1; mamax = Math.min((11 * n1) / 8 + 12, (13 * (M1[M1_itr])) / 8 + 1); deltam = (int) ((mamax - ma) / 10) + 1; if (maxIterations > 0) { maxcount = maxIterations / deltam; } } else { errInf = ErrorInfo.NO_EIGENVALUE; deltam = 0; ma = 0; } // calculation of ma ends. eigvectors.add(DoubleFactory1D.dense.make(startvector.size())); // new ritz vector is being added in eigvectors. List<Double> Tvector = new LinkedList<>(); Tvectors.add(Tvector); // new T matrix vector is being added in Tvectors. if (ma == 0) { eigvectors.get(eigvectors.size() - 1).assign(0); } if (ma != 0) { double[] eval; double[] nthEigenvector = null; do { if (ma > alpha.length) { // size of T matrix is to be increased. LanczosIteration iter = new LanczosIterationFixed(ma); generateTMatrix(iter); } count++; // on return, z contains all orthonormal eigen vectors of T matrix. EigenDecomposition eigenDecomposition = new EigenDecomposition(Arrays.copyOfRange(alpha, 0, ma), Arrays.copyOfRange(beta, 0, ma)); eval = eigenDecomposition.getRealEigenvalues(); Arrays.sort(eval); // Consistent with IETL // search for the value of nth starts, where nth is the nth eigen vector in z. for (nth = ma - 1; nth >= 0; nth--) { if (Math.abs(eval[nth] - eigenvalues[pos]) <= threshold) { break; } } // search for the value of ith ends, where ith is the ith eigen vector in z. if (nth == -1) { error = 0; ma = 0; eigvectors.get(eigvectors.size() - 1).assign(0); errInf = ErrorInfo.NO_EIGENVALUE; } else { nthEigenvector = eigenDecomposition.getEigenvector(ma - 1 - nth).toArray(); error = Math.abs(beta[ma - 1] * nthEigenvector[ma - 1]); // beta[ma - 1] = betaMplusOne. if (error > errorTolerance) { ma += deltam; } } // end of else } while (error > errorTolerance && count < maxcount); if (error > errorTolerance) { eigvectors.get(eigvectors.size() - 1).assign(0); errInf = ErrorInfo.NOT_CALCULATED; } else { // if error is small enough. if (ma != 0) { for (int i = 0; i < ma; i++) { (Tvectors.get(Tvectors.size() - 1)).add(nthEigenvector[i]); } if (ma > maMax) { maMax = ma; } lambda = eval[nth]; } } } eigenval_a.add(lambda); // for Info object. Ma.add(ma); // for Info object. status.add(errInf); M1_itr++; M2_itr++; } // end of while(in_eigvals_start != in_eigvals_end) // basis transformation of eigen vectors of T. These vectors are good // approximation of eigen vectors of actual matrix. int eigenvectors_itr; int Tvectors_itr; a_and_b = makeFirstStep(randomGenerator); if (a_and_b.getFirst() != alpha[0] || a_and_b.getSecond() != beta[0]) { throw new RuntimeException("T-matrix problem at first step"); } eigenvectors_itr = 0; Tvectors_itr = 0; while (eigenvectors_itr != end) { if (!Tvectors.get(Tvectors_itr).isEmpty()) { DoubleMatrix1D eigvector = eigvectors.get(eigenvectors_itr); eigvector.assign(0); eigvector.assign(startvector, PlusMult.plusMult(Tvectors.get(Tvectors_itr).get(0))); eigvector.assign(vec2, PlusMult.plusMult(Tvectors.get(Tvectors_itr).get(1))); } eigenvectors_itr++; Tvectors_itr++; } vec2index = 2; for (int j = 2; j < maMax; j++) { a_and_b = makeStep(j - 1); if (a_and_b.getFirst() != alpha[j - 1] || a_and_b.getSecond() != beta[j - 1]) { throw new RuntimeException("T-matrix problem"); } ++vec2index; eigenvectors_itr = 0; Tvectors_itr = 0; while (eigenvectors_itr != end) { if (Tvectors.get(Tvectors_itr).size() > j) { DoubleMatrix1D eigvector = eigvectors.get(eigenvectors_itr); eigvector.assign(vec2, PlusMult.plusMult(Tvectors.get(Tvectors_itr).get(j))); } // vec2 is being added in one vector of eigvectors. eigenvectors_itr++; Tvectors_itr++; } } // end of basis transformation. // end of basis transformation. // copying to the output iterator & residuum calculation starts: int i = 0; DoubleMatrix2D eigenvectors = DoubleFactory2D.dense.make(startvector.size(), end - start); for (eigenvectors_itr = start; eigenvectors_itr != end; eigenvectors_itr++) { DoubleMatrix1D eigvector = eigvectors.get(eigenvectors_itr); eigenvectors.viewColumn(i).assign(eigvector); matrix.zMult(eigvector, vec3); vec3.assign(eigvector, PlusMult.minusMult(eigenval_a.get(i))); // now vec3 is (A*v - eigenval_a*v); *eigenvectors_itr) is being added in vec3. residuum.add(norm2(vec3)); i++; } // copying to the output iterator ends. Info info = new Info(M1, M2, Ma, eigenval_a, residuum, status); return new Pair<>(eigenvectors, info); }
From source file:com.analog.lyric.dimple.solvers.core.parameterizedMessages.MultivariateNormalParameters.java
private void validateMatrix(double[][] m) { final int n = m.length; boolean allZero = true; for (double[] row : m) { for (double value : row) { if (value != value) { throw new DimpleException("Matrix contains a NaN value"); }/*from w w w. ja v a 2s . co m*/ if (value != 0.0) { allZero = false; } } } if (allZero) { return; } boolean infiniteDiagonal = true; for (int i = 0; i < n; ++i) { final double[] row = m[i]; if (row.length != n) { throw new DimpleException("Matrix is not square"); } for (int j = 0; j < n; ++j) { final double vij = row[j]; if (j == i) { infiniteDiagonal &= (vij == Double.POSITIVE_INFINITY); } else { final double vji = m[j][i]; if (Math.abs(vji - vij) > 1e-10) { throw new DimpleException("Matrix is not symmetric at entry (%d,%d)", i, j); } } } } if (infiniteDiagonal) { return; } if (n > 0) { EigenDecomposition eig = new EigenDecomposition(wrapRealMatrix(m)); for (double value : eig.getRealEigenvalues()) { if (value <= 0) { throw notPositiveDefinite(); } } } }