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

Example 1 with RealVector

use of org.apache.commons.math3.linear.RealVector 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 2 with RealVector

use of org.apache.commons.math3.linear.RealVector in project gatk by broadinstitute.

the class CoveragePoNQCUtils method hasSuspiciousContigs.

/**
     *  Given a single sample tangent normalization (or other coverage profile), determine whether any contig looks like
     *   it has an arm level event (defined as 25% (or more) of the contig amplified/deleted)
     *
     * @param singleSampleTangentNormalized Tangent normalized data for a single sample.
     * @return never {@code null}
     */
private static Boolean hasSuspiciousContigs(final ReadCountCollection singleSampleTangentNormalized, final Map<String, Double> contigToMedian) {
    final List<String> allContigsPresent = retrieveAllContigsPresent(singleSampleTangentNormalized);
    for (String contig : allContigsPresent) {
        final ReadCountCollection oneContigReadCountCollection = singleSampleTangentNormalized.subsetTargets(singleSampleTangentNormalized.targets().stream().filter(t -> t.getContig().equals(contig)).collect(Collectors.toSet()));
        final RealVector counts = oneContigReadCountCollection.counts().getColumnVector(0);
        for (int i = 0; i < 4; i++) {
            final RealVector partitionCounts = counts.getSubVector(i * counts.getDimension() / 4, counts.getDimension() / 4);
            final double[] partitionArray = DoubleStream.of(partitionCounts.toArray()).map(d -> Math.pow(2, d)).sorted().toArray();
            double median = new Median().evaluate(partitionArray);
            final double medianShiftInCRSpace = contigToMedian.getOrDefault(contig, 1.0) - 1.0;
            median -= medianShiftInCRSpace;
            if ((median > AMP_THRESHOLD) || (median < DEL_THRESHOLD)) {
                logger.info("Suspicious contig: " + singleSampleTangentNormalized.columnNames().get(0) + " " + contig + " (" + median + " -- " + i + ")");
                return true;
            }
        }
    }
    return false;
}

Example 3 with RealVector

use of org.apache.commons.math3.linear.RealVector in project gatk-protected by broadinstitute.

the class Nd4jApacheAdapterUtilsUnitTest method assertINDArrayToApacheVectorCorrectness.

public void assertINDArrayToApacheVectorCorrectness(final INDArray arr) {
    final RealVector result = Nd4jApacheAdapterUtils.convertINDArrayToApacheVector(arr);
    final RealVector expected = new ArrayRealVector(arr.dup().data().asDouble());
    Assert.assertEquals(result.getDimension(), expected.getDimension());
    ArrayAsserts.assertArrayEquals(result.toArray(), expected.toArray(), EPS);
}

Example 4 with RealVector

use of org.apache.commons.math3.linear.RealVector in project knime-core by knime.

the class BinaryNominalSplitsPCA method calculateWeightedCovarianceMatrix.

/**
 * Calculates the weighted covariance matrix of the class probability vectors of the CombinedAttributeValues in
 * attVals
 *
 * @param attVals
 * @param meanClassProbabilityVec
 * @param totalWeight
 * @param numTargetVals
 * @return The weighted covariance matrix of the class probability vectors of the CombinedAttributeValues
 */
static RealMatrix calculateWeightedCovarianceMatrix(final CombinedAttributeValues[] attVals, final RealVector meanClassProbabilityVec, final double totalWeight, final int numTargetVals) {
    RealMatrix weightedCovarianceMatrix = MatrixUtils.createRealMatrix(numTargetVals, numTargetVals);
    for (CombinedAttributeValues attVal : attVals) {
        RealVector diff = attVal.m_classProbabilityVector.subtract(meanClassProbabilityVec);
        weightedCovarianceMatrix = weightedCovarianceMatrix.add(diff.outerProduct(diff).scalarMultiply(attVal.m_totalWeight));
    }
    weightedCovarianceMatrix = weightedCovarianceMatrix.scalarMultiply(1 / (totalWeight - 1));
    return weightedCovarianceMatrix;
}

Example 5 with RealVector

use of org.apache.commons.math3.linear.RealVector in project knime-core by knime.

the class BinaryNominalSplitsPCA method calculatePCAOrdering.

/**
 * Applies the PCA algorithm ("Partitioning Nominal Attributes in Decision Trees", Coppersmith et al. (1999)) to
 * order the attribute values according to their principal component score.
 *
 * @param attVals
 * @param totalWeight
 * @param numTargetVals
 * @return a CombinedAttributeValues array ordered by the principal component score of its elements.
 */
static CombinedAttributeValues[] calculatePCAOrdering(final CombinedAttributeValues[] attVals, final double totalWeight, final int numTargetVals) {
    RealVector meanClassProbabilityVec = calculateMeanClassProbabilityVector(attVals, totalWeight, numTargetVals);
    RealMatrix weightedCovarianceMatrix = calculateWeightedCovarianceMatrix(attVals, meanClassProbabilityVec, totalWeight, numTargetVals);
    // eigenvalue decomposition
    EigenDecomposition eig;
    try {
        eig = new EigenDecomposition(weightedCovarianceMatrix);
    } catch (Exception e) {
        return null;
    }
    // get principal component
    RealVector principalComponent = eig.getEigenvector(argMax(eig.getRealEigenvalues()));
    // calculate principal component scores
    for (CombinedAttributeValues attVal : attVals) {
        attVal.m_principalComponentScore = principalComponent.dotProduct(attVal.m_classProbabilityVector);
    }
    // sort attribute values list ascending by principal component score
    Comparator<CombinedAttributeValues> pcScoreComparator = new Comparator<CombinedAttributeValues>() {

        @Override
        public int compare(final CombinedAttributeValues arg0, final CombinedAttributeValues arg1) {
            if (arg0.m_principalComponentScore < arg1.m_principalComponentScore) {
                return -1;
            } else if (arg0.m_principalComponentScore > arg1.m_principalComponentScore) {
                return 1;
            } else {
                return 0;
            }
        }
    };
    CombinedAttributeValues[] result = Arrays.copyOf(attVals, attVals.length);
    Arrays.sort(result, pcScoreComparator);
    return result;
}