List of usage examples for org.apache.commons.math3.linear RealMatrix multiply
RealMatrix multiply(RealMatrix m) throws DimensionMismatchException;
From source file:org.cirdles.squid.algorithms.weightedMeans.WeightedMeanCalculators.java
/** * * @param values/* w w w .j ava 2s. co m*/ * @param varCov * @return */ public static WtdAvCorrResults wtdAvCorr(double[] values, double[][] varCov) { // assume varCov is variance-covariance matrix (i.e. SigRho = false) WtdAvCorrResults results = new WtdAvCorrResults(); int n = varCov.length; RealMatrix omegaInv = new BlockRealMatrix(varCov); RealMatrix omega = MatrixUtils.inverse(omegaInv); double numer = 0.0; double denom = 0.0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { numer += (values[i] + values[j]) * omega.getEntry(i, j); denom += omega.getEntry(i, j); } } // test denom if (denom > 0.0) { double meanVal = numer / denom / 2.0; double meanValSigma = Math.sqrt(1.0 / denom); double[][] unwtdResidsArray = new double[n][1]; for (int i = 0; i < n; i++) { unwtdResidsArray[i][0] = values[i] - meanVal; } RealMatrix unwtdResids = new BlockRealMatrix(unwtdResidsArray); RealMatrix transUnwtdResids = unwtdResids.transpose(); RealMatrix product = transUnwtdResids.multiply(omega); RealMatrix sumWtdResids = product.multiply(unwtdResids); double mswd = 0.0; double prob = 0.0; if (n > 1) { mswd = sumWtdResids.getEntry(0, 0) / (n - 1); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 1), 1E9); prob = 1.0 - fdist.cumulativeProbability(mswd); } results.setBad(false); results.setMeanVal(meanVal); results.setSigmaMeanVal(meanValSigma); results.setMswd(mswd); results.setProb(prob); } return results; }
From source file:org.easotope.shared.math.ExponentialDecayFitter.java
private void calculateIfNecessary() { if (!needsCalculation) { return;//from www . j a va2 s . c o m } a = b = c = Double.NaN; if (pointList.size() < 2) { return; } Point[] point = pointList.toArray(new Point[pointList.size()]); Arrays.sort(point, new PointComparator()); double[] s = new double[point.length]; s[0] = 0.0; for (int k = 1; k < point.length; k++) { s[k] = s[k - 1] + 0.5 * (point[k].y + point[k - 1].y) * (point[k].x - point[k - 1].x); } // calculate m RealMatrix m = MatrixUtils.createRealMatrix(2, 2); for (int k = 0; k < point.length; k++) { m.setEntry(0, 0, m.getEntry(0, 0) + (point[k].x - point[0].x) * (point[k].x - point[0].x)); m.setEntry(0, 1, m.getEntry(0, 1) + (point[k].x - point[0].x) * s[k]); m.setEntry(1, 1, m.getEntry(1, 1) + s[k] * s[k]); } m.setEntry(1, 0, m.getEntry(0, 1)); // invert m try { m = MatrixUtils.inverse(m); } catch (Exception e) { return; } // calculate n RealMatrix n = MatrixUtils.createRealMatrix(2, 1); for (int k = 0; k < point.length; k++) { n.setEntry(0, 0, n.getEntry(0, 0) + (point[k].y - point[0].y) * (point[k].x - point[0].x)); n.setEntry(1, 0, n.getEntry(1, 0) + (point[k].y - point[0].y) * s[k]); } // calculate c c = m.multiply(n).getEntry(1, 0); // calculate m m = MatrixUtils.createRealMatrix(2, 2); m.setEntry(0, 0, point.length); for (int k = 0; k < point.length; k++) { m.setEntry(0, 1, m.getEntry(0, 1) + Math.exp(c * point[k].x)); m.setEntry(1, 1, m.getEntry(1, 1) + Math.exp(2.0 * c * point[k].x)); } m.setEntry(1, 0, m.getEntry(0, 1)); // invert m try { m = MatrixUtils.inverse(m); } catch (Exception e) { return; } // calculate n n = MatrixUtils.createRealMatrix(2, 1); for (int k = 0; k < point.length; k++) { n.setEntry(0, 0, n.getEntry(0, 0) + point[k].y); n.setEntry(1, 0, n.getEntry(1, 0) + point[k].y * Math.exp(c * point[k].x)); } // calculate a & b RealMatrix r = m.multiply(n); a = r.getEntry(0, 0); b = r.getEntry(1, 0); }
From source file:org.knime.al.util.noveltydetection.knfst.KNFST.java
public static RealMatrix projection(final RealMatrix kernelMatrix, final String[] labels) throws KNFSTException { final ClassWrapper[] classes = ClassWrapper.classes(labels); // check labels if (classes.length == 1) { throw new IllegalArgumentException( "not able to calculate a nullspace from data of a single class using KNFST (input variable \"labels\" only contains a single value)"); }/* ww w . j av a2 s . c o m*/ // check kernel matrix if (!kernelMatrix.isSquare()) { throw new IllegalArgumentException("The KernelMatrix must be quadratic!"); } // calculate weights of orthonormal basis in kernel space final RealMatrix centeredK = centerKernelMatrix(kernelMatrix); final EigenDecomposition eig = new EigenDecomposition(centeredK); final double[] eigVals = eig.getRealEigenvalues(); final ArrayList<Integer> nonZeroEigValIndices = new ArrayList<Integer>(); for (int i = 0; i < eigVals.length; i++) { if (eigVals[i] > 1e-12) { nonZeroEigValIndices.add(i); } } int eigIterator = 0; final RealMatrix eigVecs = eig.getV(); RealMatrix basisvecs = null; try { basisvecs = MatrixUtils.createRealMatrix(eigVecs.getRowDimension(), nonZeroEigValIndices.size()); } catch (final Exception e) { throw new KNFSTException("Something went wrong. Try different parameters or a different kernel."); } for (final Integer index : nonZeroEigValIndices) { final double normalizer = 1 / Math.sqrt(eigVals[index]); final RealVector basisVec = eigVecs.getColumnVector(eigIterator).mapMultiply(normalizer); basisvecs.setColumnVector(eigIterator++, basisVec); } // calculate transformation T of within class scatter Sw: // T= B'*K*(I-L) and L a block matrix final RealMatrix L = kernelMatrix.createMatrix(kernelMatrix.getRowDimension(), kernelMatrix.getColumnDimension()); int start = 0; for (final ClassWrapper cl : classes) { final int count = cl.getCount(); L.setSubMatrix(MatrixFunctions.ones(count, count).scalarMultiply(1.0 / count).getData(), start, start); start += count; } // need Matrix M with all entries 1/m to modify basisvecs which allows // usage of // uncentered kernel values (eye(size(M)).M)*basisvecs final RealMatrix M = MatrixFunctions .ones(kernelMatrix.getColumnDimension(), kernelMatrix.getColumnDimension()) .scalarMultiply(1.0 / kernelMatrix.getColumnDimension()); final RealMatrix I = MatrixUtils.createRealIdentityMatrix(M.getColumnDimension()); // compute helper matrix H final RealMatrix H = ((I.subtract(M)).multiply(basisvecs)).transpose().multiply(kernelMatrix) .multiply(I.subtract(L)); // T = H*H' = B'*Sw*B with B=basisvecs final RealMatrix T = H.multiply(H.transpose()); // calculate weights for null space RealMatrix eigenvecs = MatrixFunctions.nullspace(T); if (eigenvecs == null) { final EigenDecomposition eigenComp = new EigenDecomposition(T); final double[] eigenvals = eigenComp.getRealEigenvalues(); eigenvecs = eigenComp.getV(); final int minId = MatrixFunctions.argmin(MatrixFunctions.abs(eigenvals)); final double[] eigenvecsData = eigenvecs.getColumn(minId); eigenvecs = MatrixUtils.createColumnRealMatrix(eigenvecsData); } // System.out.println("eigenvecs:"); // test.printMatrix(eigenvecs); // calculate null space projection final RealMatrix proj = ((I.subtract(M)).multiply(basisvecs)).multiply(eigenvecs); return proj; }
From source file:org.micromanager.plugins.magellan.coordinates.AffineCalibrator.java
private AffineTransform computeAffine(Point2D.Double pix1, Point2D.Double pix2, Point2D.Double pix3, Point2D.Double stage1, Point2D.Double stage2, Point2D.Double stage3) { //solve system x' = xA for A, x is stage points and x' is pixel points //6x6 matrix// w w w . j ava 2 s . c o m RealMatrix pixPoints = new Array2DRowRealMatrix(new double[][] { { pix1.x, pix1.y, 1, 0, 0, 0 }, { 0, 0, 0, pix1.x, pix1.y, 1 }, { pix2.x, pix2.y, 1, 0, 0, 0 }, { 0, 0, 0, pix2.x, pix2.y, 1 }, { pix3.x, pix3.y, 1, 0, 0, 0 }, { 0, 0, 0, pix3.x, pix3.y, 1 } }); //6x1 matrix RealMatrix stagePoints = new Array2DRowRealMatrix( new double[] { stage1.x, stage1.y, stage2.x, stage2.y, stage3.x, stage3.y }); //invert stagePoints matrix RealMatrix stagePointsInv = new LUDecomposition(pixPoints).getSolver().getInverse(); RealMatrix A = stagePointsInv.multiply(stagePoints); AffineTransform transform = new AffineTransform( new double[] { A.getEntry(0, 0), A.getEntry(3, 0), A.getEntry(1, 0), A.getEntry(4, 0) }); return transform; }
From source file:org.nd4j.linalg.Nd4jTestsComparisonFortran.java
@Test public void testGemvApacheCommons() { int[] rowsArr = new int[] { 4, 4, 4, 8, 8, 8 }; int[] colsArr = new int[] { 2, 1, 10, 2, 1, 10 }; for (int x = 0; x < rowsArr.length; x++) { int rows = rowsArr[x]; int cols = colsArr[x]; List<Pair<INDArray, String>> matrices = NDArrayCreationUtil.getAllTestMatricesWithShape(rows, cols, 12345);// ww w. j a va 2s . co m List<Pair<INDArray, String>> vectors = NDArrayCreationUtil.getAllTestMatricesWithShape(cols, 1, 12345); for (int i = 0; i < matrices.size(); i++) { for (int j = 0; j < vectors.size(); j++) { Pair<INDArray, String> p1 = matrices.get(i); Pair<INDArray, String> p2 = vectors.get(j); String errorMsg = getTestWithOpsErrorMsg(i, j, "mmul", p1, p2); INDArray m = p1.getFirst(); INDArray v = p2.getFirst(); RealMatrix rm = new BlockRealMatrix(m.rows(), m.columns()); for (int r = 0; r < m.rows(); r++) { for (int c = 0; c < m.columns(); c++) { double d = m.getDouble(r, c); rm.setEntry(r, c, d); } } RealMatrix rv = new BlockRealMatrix(cols, 1); for (int r = 0; r < v.rows(); r++) { double d = v.getDouble(r, 0); rv.setEntry(r, 0, d); } INDArray gemv = m.mmul(v); RealMatrix gemv2 = rm.multiply(rv); assertArrayEquals(new int[] { rows, 1 }, gemv.shape()); assertArrayEquals(new int[] { rows, 1 }, new int[] { gemv2.getRowDimension(), gemv2.getColumnDimension() }); //Check entries: for (int r = 0; r < rows; r++) { double exp = gemv2.getEntry(r, 0); double act = gemv.getDouble(r, 0); assertEquals(errorMsg, exp, act, 1e-5); } } } } }
From source file:org.ojalgo.benchmark.contestant.ACM.java
@Override public BenchmarkContestant<RealMatrix>.MatrixMultiplier getMatrixMultiplier() { return new MatrixMultiplier() { @Override// w ww . j a v a 2s.c om public RealMatrix apply(final RealMatrix left, final RealMatrix right) { return left.multiply(right); } }; }
From source file:org.ojalgo.benchmark.contestant.ACM.java
@Override public BenchmarkContestant<RealMatrix>.TransposedMultiplier getTransposedMultiplier() { return new TransposedMultiplier() { @Override//from w ww . j a va 2s . c o m public RealMatrix apply(final RealMatrix left, final RealMatrix right) { return left.multiply(right.transpose()); } }; }
From source file:org.ojalgo.benchmark.lab.library.ACM.java
@Override protected RealMatrix multiply(final RealMatrix... factors) { RealMatrix retVal = factors[0]; for (int f = 1; f < factors.length; f++) { retVal = retVal.multiply(factors[f]); }/*from w w w . ja v a 2 s . c om*/ return retVal; }
From source file:org.pmad.gmm.MyEDC.java
/** * Computes the square-root of the matrix. * This implementation assumes that the matrix is symmetric and positive * definite.//from w ww .ja v a2 s .c o m * * @return the square-root of the matrix. * @throws MathUnsupportedOperationException if the matrix is not * symmetric or not positive definite. * @since 3.1 */ public RealMatrix getSquareRoot() { if (!isSymmetric) { throw new MathUnsupportedOperationException(); } final double[] sqrtEigenValues = new double[realEigenvalues.length]; for (int i = 0; i < realEigenvalues.length; i++) { final double eigen = realEigenvalues[i]; if (eigen <= 0) { throw new MathUnsupportedOperationException(); } sqrtEigenValues[i] = FastMath.sqrt(eigen); } final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues); final RealMatrix v = getV(); final RealMatrix vT = getVT(); return v.multiply(sqrtEigen).multiply(vT); }
From source file:org.pmad.gmm.MyMixMNDEM.java
/** * Fit a mixture model to the data supplied to the constructor. * * The quality of the fit depends on the concavity of the data provided to * the constructor and the initial mixture provided to this function. If the * data has many local optima, multiple runs of the fitting function with * different initial mixtures may be required to find the optimal solution. * If a SingularMatrixException is encountered, it is possible that another * initialization would work./*from w ww . j a v a2 s . c om*/ * * @param initialMixture Model containing initial values of weights and * multivariate normals * @param maxIterations Maximum iterations allowed for fit * @param threshold Convergence threshold computed as difference in * logLikelihoods between successive iterations * @throws SingularMatrixException if any component's covariance matrix is * singular during fitting * @throws NotStrictlyPositiveException if numComponents is less than one * or threshold is less than Double.MIN_VALUE * @throws DimensionMismatchException if initialMixture mean vector and data * number of columns are not equal */ public void fit(final MyMixMND initialMixture, final int maxIterations, final double threshold) throws SingularMatrixException, NotStrictlyPositiveException, DimensionMismatchException { if (maxIterations < 1) { throw new NotStrictlyPositiveException(maxIterations); } if (threshold < Double.MIN_VALUE) { throw new NotStrictlyPositiveException(threshold); } final int n = data.length; // Number of data columns. Jagged data already rejected in constructor, // so we can assume the lengths of each row are equal. final int numCols = data[0].length; final int k = initialMixture.getComponents().size(); final int numMeanColumns = initialMixture.getComponents().get(0).getSecond().getMeans().length; if (numMeanColumns != numCols) { throw new DimensionMismatchException(numMeanColumns, numCols); } int numIterations = 0; double previousLogLikelihood = 0d; logLikelihood = Double.NEGATIVE_INFINITY; // Initialize model to fit to initial mixture. fittedModel = new MyMixMND(initialMixture.getComponents()); while (numIterations++ <= maxIterations && FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { System.out.println(numIterations); previousLogLikelihood = logLikelihood; double sumLogLikelihood = 0d; // Mixture components final List<Pair<Double, MyMND>> components = fittedModel.getComponents(); // Weight and distribution of each component final double[] weights = new double[k]; final MyMND[] mvns = new MyMND[k]; for (int j = 0; j < k; j++) { weights[j] = components.get(j).getFirst(); mvns[j] = components.get(j).getSecond(); } // E-step: compute the data dependent parameters of the expectation // function. // The percentage of row's total density between a row and a // component final double[][] gamma = new double[n][k]; // Sum of gamma for each component final double[] gammaSums = new double[k]; // for (double d : weights) { // System.out.print(d+" "); // } // System.out.println(); // Sum of gamma times its row for each each component final double[][] gammaDataProdSums = new double[k][numCols]; for (int i = 0; i < n; i++) { final double rowDensity = fittedModel.density(data[i]); sumLogLikelihood += FastMath.log(rowDensity); for (int j = 0; j < k; j++) { gamma[i][j] = weights[j] * mvns[j].density(data[i]); if (rowDensity == 0) { if (gamma[i][j] == 0) { gamma[i][j] = weights[j]; } else { System.err.println("ele"); } } else { gamma[i][j] /= rowDensity; } gammaSums[j] += gamma[i][j]; for (int col = 0; col < numCols; col++) { gammaDataProdSums[j][col] += gamma[i][j] * data[i][col]; } } } if (Utils.hasNaN(gamma)) { System.out.println("gamma has NaN"); } logLikelihood = sumLogLikelihood / n; // M-step: compute the new parameters based on the expectation // function. final double[] newWeights = new double[k]; final double[][] newMeans = new double[k][numCols]; for (int j = 0; j < k; j++) { newWeights[j] = gammaSums[j] / n; for (int col = 0; col < numCols; col++) { newMeans[j][col] = gammaDataProdSums[j][col] / gammaSums[j]; } } // Compute new covariance matrices final RealMatrix[] newCovMats = new RealMatrix[k]; for (int j = 0; j < k; j++) { newCovMats[j] = new Array2DRowRealMatrix(numCols, numCols); } for (int i = 0; i < n; i++) { for (int j = 0; j < k; j++) { final RealMatrix vec = new Array2DRowRealMatrix(MathArrays.ebeSubtract(data[i], newMeans[j])); final RealMatrix dataCov = vec.multiply(vec.transpose()).scalarMultiply(gamma[i][j]); newCovMats[j] = newCovMats[j].add(dataCov); } } // Converting to arrays for use by fitted model final double[][][] newCovMatArrays = new double[k][numCols][numCols]; for (int j = 0; j < k; j++) { newCovMats[j] = newCovMats[j].scalarMultiply(1d / gammaSums[j]); newCovMatArrays[j] = newCovMats[j].getData(); // if (Utils.hasNaN(newCovMatArrays[j])) { // System.out.println("covmet nan"); // } } // if (Utils.hasNaN(newMeans)) { // System.out.println("neams nan"); // } // Update current model fittedModel = new MyMixMND(newWeights, newMeans, newCovMatArrays); } if (FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { // Did not converge before the maximum number of iterations // throw new ConvergenceException(); System.out.println("Did not converge"); } }