List of usage examples for org.apache.commons.math3.linear RealMatrix scalarMultiply
RealMatrix scalarMultiply(double d);
From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * Pre-calculate intermediate matrices, perform setup, etc. *//* w ww . ja v a2 s. c o m*/ 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:org.moeaframework.core.operator.real.AdaptiveMetropolisTest.java
/** * This is sample output from running Jasper Vrugt's AMALGAM codes in * MATLAB that is positive definite./*from w w w . ja v a2 s . c om*/ */ @Test public void testVrugtSample() { RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] { { 0.0000, 0.0002, -0.0006, 0.0008, 0.0001, 0.0018, 0.0001, 0.0003, 0.0006, -0.0035 }, { 0.0934, -0.0010, -0.0062, 0.0001, 0.0024, 0.0017, 0.0019, 0.0004, 0.0008, -0.0067 }, { 0.1435, -0.0047, 0.0007, -0.0051, -0.0049, 0.0054, 0.0052, 0.0018, 0.0016, 0.0014 }, { 0.0026, 0.0001, -0.0001, 0.0002, 0.0001, 0.0004, -0.0002, -0.0002, 0.0000, -0.0010 }, { 0.0262, 0.0028, -0.0003, 0.0019, -0.0013, 0.0002, -0.0009, -0.0004, -0.0012, -0.0021 }, { 0.0432, 0.0005, -0.0006, 0.0008, 0.0001, 0.0018, -0.0015, 0.0003, 0.0006, -0.0035 }, { 0.1220, -0.0028, -0.0082, 0.0006, 0.0033, -0.0044, -0.0019, -0.0033, -0.0029, 0.0010 }, { 0.0657, -0.0049, -0.0005, 0.0004, 0.0016, 0.0020, -0.0010, 0.0047, -0.0047, -0.0033 }, { 0.1594, 0.0001, -0.0001, 0.0026, 0.0058, 0.0004, -0.0001, -0.0002, 0.0000, -0.0018 }, { 0.1329, 0.0057, -0.0010, -0.0053, 0.0026, 0.0048, 0.0061, 0.0019, 0.0023, -0.0020 }, { 0.0998, 0.0029, 0.0036, -0.0014, 0.0015, 0.0018, -0.0035, 0.0006, -0.0013, 0.0083 }, { 0.0160, -0.0003, -0.0009, 0.0008, 0.0006, -0.0003, -0.0010, 0.0008, -0.0005, -0.0033 }, { 0.0061, -0.0001, -0.0001, -0.0001, 0.0001, 0.0001, 0.0001, 0.0000, 0.0000, -0.0003 } }); double jumpRate = Math.pow(2.4 / Math.sqrt(10), 2.0); RealMatrix actual = new CholeskyDecomposition( new Covariance(matrix.scalarMultiply(jumpRate)).getCovarianceMatrix()).getLT(); RealMatrix expected = MatrixUtils.createRealMatrix(new double[][] { { 0.0335, -0.0001, -0.0003, -0.0006, 0.0005, 0.0003, 0.0006, 0.0000, 0.0001, 0.0005 }, { 0.0000, 0.0017, 0.0005, -0.0002, 0.0004, 0.0002, 0.0001, -0.0002, 0.0005, 0.0003 }, { 0.0000, 0.0000, 0.0016, -0.0003, -0.0006, 0.0008, -0.0000, 0.0006, 0.0001, 0.0010 }, { 0.0000, 0.0000, 0.0000, 0.0012, 0.0008, -0.0007, -0.0010, -0.0003, -0.0004, -0.0001 }, { 0.0000, 0.0000, 0.0000, 0.0000, 0.0009, -0.0001, -0.0001, 0.0003, -0.0002, -0.0001 }, { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0008, 0.0006, 0.0005, 0.0004, -0.0015 }, { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0007, -0.0001, 0.0004, -0.0006 }, { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0005, -0.0005, -0.0003 }, { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0003, -0.0003 }, { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0003 } }); TestUtils.assertEquals(expected, actual, new AbsoluteError(0.0001)); }
From source file:org.moeaframework.core.operator.real.AdaptiveMetropolisTest.java
/** * Returns the covariance matrix for the specified solutions. * // w w w.ja v a2 s. com * @param cluster the cluster * @param jumpRateCoefficient the jump rate coefficient * @return the covariance matrix for the specified solutions */ private RealMatrix getCovariance(Solution[] parents, double jumpRateCoefficient) { RealMatrix rm = MatrixUtils.createRealMatrix(parents.length, 2); for (int i = 0; i < parents.length; i++) { rm.setRow(i, EncodingUtils.getReal(parents[i])); } rm = rm.scalarMultiply(Math.pow(jumpRateCoefficient / Math.sqrt(2), 2.0)); return new Covariance(rm).getCovarianceMatrix(); }
From source file:org.rhwlab.BHCnotused.GaussianGIWPrior.java
@Override public RealMatrix getPrecision() { RealMatrix ret = new Array2DRowRealMatrix(m.getDimension(), m.getDimension()); int n = 0;//from www.j a v a2 s . c om for (LabeledFieldVector v : data) { int label = v.getLabel(); RealVector[] vs = source.getClusterVectors(label); for (int i = 0; i < vs.length; ++i) { RealVector del = vs[i].subtract(getMean()); ret = ret.add(del.outerProduct(del)); ++n; } } ret = ret.scalarMultiply(1.0 / n); LUDecomposition lud = new LUDecomposition(ret); RealMatrix prec = lud.getSolver().getInverse(); return prec; }
From source file:org.rhwlab.dispim.datasource.ClusteredDataSource.java
public static void main(String[] args) throws Exception { ClusteredDataSource source = new ClusteredDataSource( "/nfs/waterston/pete/Segmentation/Cherryimg_SimpleSegmentation0075.xml"); GaussianMixture gm = new GaussianMixture(); gm.setSource(source);/* ww w .j a v a2 s. c o m*/ gm.setIterationMax(600); gm.setAlpha0(0.001); gm.setBeta0(0.000001); gm.setNu0(400.0); gm.setMu0(source.getDataMean()); RealMatrix W0 = MatrixUtils.createRealIdentityMatrix(source.getD()); W0 = W0.scalarMultiply(0.00001); gm.setW0(W0); gm.init(source.getClusterCount()); gm.run(); gm.saveAsXML("/nfs/waterston/pete/Segmentation/Cherryimg75_SimpleSegmentation_GMM0075.xml"); }
From source file:org.rhwlab.dispim.datasource.GaussianComponent.java
public RealMatrix precision(RealVector mu) { RealMatrix ret = new Array2DRowRealMatrix(mu.getDimension(), mu.getDimension()); ret.scalarMultiply(0.0); // make sure it is zero int n = 0;//from w ww. j av a 2s . c om for (Integer index : indexes) { RealVector x = source.get(index).coords; RealVector del = x.subtract(mu); ret = ret.add(del.outerProduct(del)); ++n; } ret = ret.scalarMultiply(1.0 / n); LUDecomposition lud = new LUDecomposition(ret); RealMatrix prec = lud.getSolver().getInverse(); return prec; }
From source file:org.rhwlab.dispim.datasource.MicroCluster.java
License:asdf
public static RealMatrix precision(List<MicroCluster> data, RealVector mu) { RealMatrix ret = new Array2DRowRealMatrix(mu.getDimension(), mu.getDimension()); RealVector v = new ArrayRealVector(mu.getDimension()); long n = 0;// w ww. ja v a2s . c o m for (MicroCluster micro : data) { for (int p = 0; p < micro.points.length; ++p) { for (int d = 0; d < mu.getDimension(); ++d) { v.setEntry(d, micro.points[p][d]); } RealVector del = v.subtract(mu); ret = ret.add(del.outerProduct(del)); ++n; } } ret = ret.scalarMultiply(1.0 / n); LUDecomposition lud = new LUDecomposition(ret); RealMatrix prec = null; if (lud.getSolver().isNonSingular()) { prec = lud.getSolver().getInverse(); } return prec; }
From source file:org.rhwlab.dispim.datasource.MicroClusterDataSource.java
public RealMatrix getDataVariance() { RealMatrix ret = new Array2DRowRealMatrix(getD(), getD()); for (int k = 0; k < micros.length; ++k) { RealVector v = micros[k].asRealVector(); ret = ret.add(v.outerProduct(v)); }/* w w w.ja v a 2s .c om*/ return ret.scalarMultiply(1.0 / micros.length); }
From source file:org.rhwlab.dispim.nucleus.NamedNucleusFile.java
License:asdf
static public RealMatrix rotationMatrix(Vector3D A, Vector3D B) { Vector3D a = A.normalize();/*from ww w . jav a 2 s. c o m*/ Vector3D b = B.normalize(); Vector3D v = a.crossProduct(b); double s = v.getNormSq(); double c = a.dotProduct(b); RealMatrix vx = MatrixUtils.createRealMatrix(3, 3); vx.setEntry(1, 0, v.getZ()); vx.setEntry(0, 1, -v.getZ()); vx.setEntry(2, 0, -v.getY()); vx.setEntry(0, 2, v.getY()); vx.setEntry(2, 1, v.getX()); vx.setEntry(1, 2, -v.getX()); RealMatrix vx2 = vx.multiply(vx); RealMatrix scaled = vx2.scalarMultiply((1.0 - c) / s); RealMatrix ident = MatrixUtils.createRealIdentityMatrix(3); RealMatrix sum = vx.add(scaled); RealMatrix ret = ident.add(sum); return ret; }
From source file:org.rhwlab.variationalbayesian.GaussianMixture.java
public void statistics() { // calculate N calculateN();/* w w w . j ava2 s. com*/ // calculate xBar for (int k = 0; k < K; ++k) { if (inModel[k]) { xBar[k].set(0.0); for (int n = 0; n < X.getN(); ++n) { Voxel vox = X.get(n); double intensity = vox.getAdjusted(); RealVector x = vox.coords; // double rnk = (intensity-background)*r.getEntry(n, k); double rnk = intensity * r.getEntry(n, k); xBar[k] = xBar[k].add(x.mapMultiply(rnk)); } xBar[k].mapDivideToSelf(N[k]); } } // calculate the S for (int k = 0; k < K; ++k) { if (inModel[k]) { if (k == 31) { int iusahdfis = 90; } for (int row = 0; row < X.getD(); ++row) { for (int col = 0; col < X.getD(); ++col) { S[k].setEntry(row, col, 0.0); } } for (int n = 0; n < X.getN(); ++n) { Voxel vox = X.get(n); RealVector x = vox.coords; RealVector del = x.subtract(xBar[k]); RealMatrix mat = del.outerProduct(del); double rnk = r.getEntry(n, k); if (rnk > 0) { int iuhsafuid = 0; } // S[k] = S[k].add(mat.scalarMultiply((intensity-background)*rnk)); S[k] = S[k].add(mat.scalarMultiply(vox.getAdjusted() * rnk)); } S[k] = S[k].scalarMultiply(1.0 / N[k]); } } int uisahdfius = 0; }