Examples with RealMatrix org.apache.commons.math3.linear.RealMatrix used on opensource projects

Example 1 with RealMatrix

use of org.apache.commons.math3.linear.RealMatrix in project jstructure by JonStargaryen.

the class MultiDimensionalScaling method computeEmbedding.

public List<double[]> computeEmbedding(RealMatrix distanceMap, int targetDimension) {
    numberOfDataPoints = distanceMap.getRowDimension();
    this.targetDimension = targetDimension;
    if (this.targetDimension > numberOfDataPoints) {
        throw new IllegalArgumentException("target dimension must not exceed number of data points");
    }
    this.distanceMap = distanceMap;
    proximityMap = computeSquaredProximityMap(this.distanceMap);
    centeringMap = computeConfiguration(proximityMap);
    normalizedEigenvectors = new ArrayList<>();
    embedding = new ArrayList<>();
    EigenDecomposition evd = new EigenDecomposition(centeringMap);
    // we are looking for the m biggest eigenvalues - they are at the last elements of the matrix
    RealMatrix eigenvectors = evd.getV();
    // Matrix eigenvalues = evd.getD();
    double[] eigenvalues = evd.getRealEigenvalues();
    Map<Integer, Double> eigenvalueMap = new HashMap<>();
    for (int eigenvalueIndex = 0; eigenvalueIndex < eigenvalues.length; eigenvalueIndex++) {
        eigenvalueMap.put(eigenvalueIndex, eigenvalues[eigenvalueIndex]);
    }
    List<Entry<Integer, Double>> sortedEigenvalues = entriesSortedByValues(eigenvalueMap).subList(0, targetDimension);
    // normalize eigenvectors
    for (Entry<Integer, Double> sortedEigenvalue : sortedEigenvalues) {
        if (sortedEigenvalue.getValue() <= 0) {
            throw new IllegalArgumentException("eigenvalue is negative: " + sortedEigenvalue.getValue());
        }
        double[] eigenvector = eigenvectors.getColumn(sortedEigenvalue.getKey());
        normalizedEigenvectors.add(normalize(eigenvector, Math.sqrt(sortedEigenvalue.getValue())));
    }
    // compose embedded data points from normalized eigenvectors
    for (int dataPointIndex = 0; dataPointIndex < numberOfDataPoints; dataPointIndex++) {
        double[] dataPoint = new double[this.targetDimension];
        for (int dataPointDimension = 0; dataPointDimension < this.targetDimension; dataPointDimension++) {
            dataPoint[dataPointDimension] = normalizedEigenvectors.get(dataPointDimension)[dataPointIndex];
        }
        embedding.add(dataPoint);
    }
    return embedding;
}

Example 2 with RealMatrix

use of org.apache.commons.math3.linear.RealMatrix in project jstructure by JonStargaryen.

the class MultiDimensionalScaling method computeConfiguration.

private RealMatrix computeConfiguration(RealMatrix proximityMap) {
    RealMatrix centeringMap = new Array2DRowRealMatrix(numberOfDataPoints, numberOfDataPoints);
    double[] rowAverage = new double[numberOfDataPoints];
    double[] columnAverage = new double[numberOfDataPoints];
    double overallAverage = 0;
    // assess rows and overall average
    for (int row = 0; row < numberOfDataPoints; row++) {
        double tempRowAverage = 0;
        for (int column = 0; column < numberOfDataPoints; column++) {
            double entry = proximityMap.getEntry(row, column);
            tempRowAverage += entry;
            overallAverage += entry;
        }
        rowAverage[row] = tempRowAverage / numberOfDataPoints;
    }
    overallAverage /= numberOfDataPoints * numberOfDataPoints;
    // assess columns
    for (int column = 0; column < numberOfDataPoints; column++) {
        double tempColumnAverage = 0;
        for (int row = 0; row < numberOfDataPoints; row++) {
            tempColumnAverage += proximityMap.getEntry(row, column);
        }
        columnAverage[column] = tempColumnAverage / numberOfDataPoints;
    }
    for (int row = 0; row < numberOfDataPoints; row++) {
        for (int column = 0; column < numberOfDataPoints; column++) {
            // b_ij = a_ij - a_i* - a_j* + a_**
            centeringMap.setEntry(row, column, proximityMap.getEntry(row, column) - rowAverage[row] - columnAverage[column] + overallAverage);
        }
    }
    return centeringMap;
}

Example 3 with RealMatrix

use of org.apache.commons.math3.linear.RealMatrix in project jstructure by JonStargaryen.

the class SVDSuperimposer method align.

@Override
public StructureAlignmentResult align(AtomContainer reference, AtomContainer query) {
    AtomContainer originalReference = reference;
    AtomContainer originalCandidate = query;
    Pair<GroupContainer, GroupContainer> atomContainerPair = AbstractAlignmentAlgorithm.comparableGroupContainerPair(reference, query, minimalSetOfAtomNames, maximalSetOfAtomNames);
    reference = atomContainerPair.getLeft();
    query = atomContainerPair.getRight();
    // calculate centroids
    double[] centroid1 = reference.calculate().centroid().getValue();
    double[] centroid2 = query.calculate().centroid().getValue();
    // center atoms
    reference.calculate().center();
    query.calculate().center();
    // compose covariance matrix and calculate SVD
    RealMatrix matrix1 = convertToMatrix(reference);
    RealMatrix matrix2 = convertToMatrix(query);
    RealMatrix covariance = matrix2.transpose().multiply(matrix1);
    SingularValueDecomposition svd = new SingularValueDecomposition(covariance);
    // R = (V * U')'
    RealMatrix ut = svd.getU().transpose();
    RealMatrix rotationMatrix = svd.getV().multiply(ut).transpose();
    // check if reflection
    if (new LUDecomposition(rotationMatrix).getDeterminant() < 0) {
        RealMatrix v = svd.getV().transpose();
        v.setEntry(2, 0, (0 - v.getEntry(2, 0)));
        v.setEntry(2, 1, (0 - v.getEntry(2, 1)));
        v.setEntry(2, 2, (0 - v.getEntry(2, 2)));
        rotationMatrix = v.transpose().multiply(ut).transpose();
    }
    double[][] rotation = rotationMatrix.getData();
    // calculate translation
    double[] translation = LinearAlgebra.on(centroid1).subtract(LinearAlgebra.on(centroid2).multiply(rotation)).getValue();
    logger.trace("rotation matrix\n{}\ntranslation vector\n{}", Arrays.deepToString(rotationMatrix.getData()), Arrays.toString(translation));
    /* transform 2nd atom select - employ neutral translation (3D vector of zeros), because the atoms are already
        * centered and calculate RMSD */
    query.calculate().transform(new Transformation(rotation));
    double rmsd = calculateRmsd(reference, query);
    // return alignment
    return new StructureAlignmentResult(originalReference, originalCandidate, query, rmsd, translation, rotation);
}

Example 4 with RealMatrix

use of org.apache.commons.math3.linear.RealMatrix in project GDSC-SMLM by aherbert.

the class ApacheLVMFitter method computeFit.

public FitStatus computeFit(double[] y, final double[] y_fit, double[] a, double[] a_dev) {
    int n = y.length;
    try {
        // Different convergence thresholds seem to have no effect on the resulting fit, only the number of
        // iterations for convergence
        final double initialStepBoundFactor = 100;
        final double costRelativeTolerance = 1e-10;
        final double parRelativeTolerance = 1e-10;
        final double orthoTolerance = 1e-10;
        final double threshold = Precision.SAFE_MIN;
        // Extract the parameters to be fitted
        final double[] initialSolution = getInitialSolution(a);
        // TODO - Pass in more advanced stopping criteria.
        // Create the target and weight arrays
        final double[] yd = new double[n];
        final double[] w = new double[n];
        for (int i = 0; i < n; i++) {
            yd[i] = y[i];
            w[i] = 1;
        }
        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
        //@formatter:off
        LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd).weight(new DiagonalMatrix(w));
        if (f instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) f).canComputeValuesAndJacobian()) {
            // Compute together, or each individually
            builder.model(new ValueAndJacobianFunction() {

                final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) f;

                public Pair<RealVector, RealMatrix> value(RealVector point) {
                    final double[] p = point.toArray();
                    final Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
                    return new Pair<RealVector, RealMatrix>(new ArrayRealVector(result.getFirst(), false), new Array2DRowRealMatrix(result.getSecond(), false));
                }

                public RealVector computeValue(double[] params) {
                    return new ArrayRealVector(fun.computeValues(params), false);
                }

                public RealMatrix computeJacobian(double[] params) {
                    return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
                }
            });
        } else {
            // Compute separately
            builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) f, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) f, a, n));
        }
        LeastSquaresProblem problem = builder.build();
        Optimum optimum = optimizer.optimize(problem);
        final double[] parameters = optimum.getPoint().toArray();
        setSolution(a, parameters);
        iterations = optimum.getIterations();
        evaluations = optimum.getEvaluations();
        if (a_dev != null) {
            try {
                double[][] covar = optimum.getCovariances(threshold).getData();
                setDeviationsFromMatrix(a_dev, covar);
            } catch (SingularMatrixException e) {
                // Matrix inversion failed. In order to return a solution 
                // return the reciprocal of the diagonal of the Fisher information 
                // for a loose bound on the limit 
                final int[] gradientIndices = f.gradientIndices();
                final int nparams = gradientIndices.length;
                GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams);
                double[][] alpha = new double[nparams][nparams];
                double[] beta = new double[nparams];
                calculator.findLinearised(nparams, y, a, alpha, beta, (NonLinearFunction) f);
                FisherInformationMatrix m = new FisherInformationMatrix(alpha);
                setDeviations(a_dev, m.crlb(true));
            }
        }
        // Compute function value
        if (y_fit != null) {
            Gaussian2DFunction f = (Gaussian2DFunction) this.f;
            f.initialise0(a);
            f.forEach(new ValueProcedure() {

                int i = 0;

                public void execute(double value) {
                    y_fit[i] = value;
                }
            });
        }
        // As this is unweighted then we can do this to get the sum of squared residuals
        // This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
        // just computes the dot product anyway.
        value = optimum.getResiduals().dotProduct(optimum.getResiduals());
    } catch (TooManyEvaluationsException e) {
        return FitStatus.TOO_MANY_EVALUATIONS;
    } catch (TooManyIterationsException e) {
        return FitStatus.TOO_MANY_ITERATIONS;
    } catch (ConvergenceException e) {
        // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
        return FitStatus.SINGULAR_NON_LINEAR_MODEL;
    } catch (Exception e) {
        // TODO - Find out the other exceptions from the fitter and add return values to match. 
        return FitStatus.UNKNOWN;
    }
    return FitStatus.OK;
}

Example 5 with RealMatrix

use of org.apache.commons.math3.linear.RealMatrix in project incubator-systemml by apache.

the class LibCommonsMath method computeEigen.

/**
	 * Function to perform Eigen decomposition on a given matrix.
	 * Input must be a symmetric matrix.
	 * 
	 * @param in matrix object
	 * @return array of matrix blocks
	 * @throws DMLRuntimeException if DMLRuntimeException occurs
	 */
private static MatrixBlock[] computeEigen(MatrixObject in) throws DMLRuntimeException {
    if (in.getNumRows() != in.getNumColumns()) {
        throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
    }
    Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
    EigenDecomposition eigendecompose = new EigenDecomposition(matrixInput);
    RealMatrix eVectorsMatrix = eigendecompose.getV();
    double[][] eVectors = eVectorsMatrix.getData();
    double[] eValues = eigendecompose.getRealEigenvalues();
    //Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
    int n = eValues.length;
    for (int i = 0; i < n; i++) {
        int k = i;
        double p = eValues[i];
        for (int j = i + 1; j < n; j++) {
            if (eValues[j] < p) {
                k = j;
                p = eValues[j];
            }
        }
        if (k != i) {
            eValues[k] = eValues[i];
            eValues[i] = p;
            for (int j = 0; j < n; j++) {
                p = eVectors[j][i];
                eVectors[j][i] = eVectors[j][k];
                eVectors[j][k] = p;
            }
        }
    }
    MatrixBlock mbValues = DataConverter.convertToMatrixBlock(eValues, true);
    MatrixBlock mbVectors = DataConverter.convertToMatrixBlock(eVectors);
    return new MatrixBlock[] { mbValues, mbVectors };
}