Examples with Evaluation org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation used on opensource projects

Example 6 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation in project GDSC-SMLM by aherbert.

the class BFGSOptimizer method bfgs.

protected PointValuePair bfgs(ConvergenceChecker<PointValuePair> checker, double[] p, LineStepSearch lineSearch) {
    final int n = p.length;
    final double EPS = epsilon;
    double[] hdg = new double[n];
    double[] xi = new double[n];
    double[][] hessian = new double[n][n];
    // Get the gradient for the the bounded point
    applyBounds(p);
    double[] g = computeObjectiveGradient(p);
    checkGradients(g, p);
    // Initialise the hessian and search direction
    for (int i = 0; i < n; i++) {
        hessian[i][i] = 1.0;
        xi[i] = -g[i];
    }
    PointValuePair current = null;
    while (true) {
        incrementIterationCount();
        // Get the value of the point
        double fp = computeObjectiveValue(p);
        if (checker != null) {
            PointValuePair previous = current;
            current = new PointValuePair(p, fp);
            if (previous != null && checker.converged(getIterations(), previous, current)) {
                // We have found an optimum.
                converged = CHECKER;
                return current;
            }
        }
        // Move along the search direction.
        final double[] pnew;
        try {
            pnew = lineSearch.lineSearch(p, fp, g, xi);
        } catch (LineSearchRoundoffException e) {
            // This can happen if the Hessian is nearly singular or non-positive-definite.
            // In this case the algorithm should be restarted.
            converged = ROUNDOFF_ERROR;
            //System.out.printf("Roundoff error, iter=%d\n", getIterations());
            return new PointValuePair(p, fp);
        }
        // We assume the new point is on/within the bounds since the line search is constrained
        double fret = lineSearch.f;
        // Test for convergence on change in position
        if (positionChecker.converged(p, pnew)) {
            converged = POSITION;
            return new PointValuePair(pnew, fret);
        }
        // Update the line direction
        for (int i = 0; i < n; i++) {
            xi[i] = pnew[i] - p[i];
        }
        p = pnew;
        // Save the old gradient
        double[] dg = g;
        // Get the gradient for the new point
        g = computeObjectiveGradient(p);
        checkGradients(g, p);
        // If necessary recompute the function value. 
        // Doing this after the gradient evaluation allows the value to be cached when 
        // computing the objective gradient
        fp = fret;
        // Test for convergence on zero gradient.
        double test = 0;
        for (int i = 0; i < n; i++) {
            final double temp = Math.abs(g[i]) * FastMath.max(Math.abs(p[i]), 1);
            //final double temp = Math.abs(g[i]);
            if (test < temp)
                test = temp;
        }
        // Compute the biggest gradient relative to the objective function
        test /= FastMath.max(Math.abs(fp), 1);
        if (test < gradientTolerance) {
            converged = GRADIENT;
            return new PointValuePair(p, fp);
        }
        for (int i = 0; i < n; i++) dg[i] = g[i] - dg[i];
        for (int i = 0; i < n; i++) {
            hdg[i] = 0.0;
            for (int j = 0; j < n; j++) hdg[i] += hessian[i][j] * dg[j];
        }
        double fac = 0, fae = 0, sumdg = 0, sumxi = 0;
        for (int i = 0; i < n; i++) {
            fac += dg[i] * xi[i];
            fae += dg[i] * hdg[i];
            sumdg += dg[i] * dg[i];
            sumxi += xi[i] * xi[i];
        }
        if (fac > Math.sqrt(EPS * sumdg * sumxi)) {
            fac = 1.0 / fac;
            final double fad = 1.0 / fae;
            for (int i = 0; i < n; i++) dg[i] = fac * xi[i] - fad * hdg[i];
            for (int i = 0; i < n; i++) {
                for (int j = i; j < n; j++) {
                    hessian[i][j] += fac * xi[i] * xi[j] - fad * hdg[i] * hdg[j] + fae * dg[i] * dg[j];
                    hessian[j][i] = hessian[i][j];
                }
            }
        }
        for (int i = 0; i < n; i++) {
            xi[i] = 0.0;
            for (int j = 0; j < n; j++) xi[i] -= hessian[i][j] * g[j];
        }
    }
}

Example 7 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation in project gatk-protected by broadinstitute.

the class CoverageDropoutDetector method retrieveGaussianMixtureModelForFilteredTargets.

/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p>
     * <p>Filters targets to populations where more than the minProportion lie in a single segment.</p>
     * <p>Returns null if no pass filtering.  Please note that in these cases,
     * in the rest of this class, we use this to assume that a GMM is not a good model.</p>
     *
     * @param segments  -- segments with segment mean in log2 copy ratio space
     * @param targets -- targets with a log2 copy ratio estimate
     * @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used
     *                      in the evaluation
     * @param numComponents -- number of components to use in the GMM.  Usually, this is 2.
     * @return  never {@code null}.  Fitting result with indications whether it converged or was even attempted.
     */
private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments, final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents) {
    // For each target in a segment that contains enough targets, normalize the difference against the segment mean
    //  and collapse the filtered targets into the copy ratio estimates.
    final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion);
    if (filteredTargetsSegDiff.size() < numComponents) {
        return new MixtureMultivariateNormalFitResult(null, false, false);
    }
    // Assume that Apache Commons wants data points in the first dimension.
    // Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API.
    final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2];
    // Convert the filtered targets into 2d array (even if second dimension is length 1).  The second dimension is
    //  uncorrelated Gaussian.  This is only to get around funny API in Apache Commons, which will throw an
    //  exception if the length of the second dimension is < 2
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
    final NormalDistribution nd = new NormalDistribution(rng, 0, .1);
    for (int i = 0; i < filteredTargetsSegDiff.size(); i++) {
        filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i);
        filteredTargetsSegDiff2d[i][1] = nd.sample();
    }
    final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents);
    final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d);
    try {
        multivariateNormalMixtureExpectationMaximization.fit(estimateEM0);
    } catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) {
        //  did not converge.  Include the model as it was when the exception was thrown.
        return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true);
    }
    return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true);
}

Example 8 with Evaluation

use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation in project gatk by broadinstitute.

the class CoverageDropoutDetector method retrieveGaussianMixtureModelForFilteredTargets.

/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p>
     * <p>Filters targets to populations where more than the minProportion lie in a single segment.</p>
     * <p>Returns null if no pass filtering.  Please note that in these cases,
     * in the rest of this class, we use this to assume that a GMM is not a good model.</p>
     *
     * @param segments  -- segments with segment mean in log2 copy ratio space
     * @param targets -- targets with a log2 copy ratio estimate
     * @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used
     *                      in the evaluation
     * @param numComponents -- number of components to use in the GMM.  Usually, this is 2.
     * @return  never {@code null}.  Fitting result with indications whether it converged or was even attempted.
     */
private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments, final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents) {
    // For each target in a segment that contains enough targets, normalize the difference against the segment mean
    //  and collapse the filtered targets into the copy ratio estimates.
    final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion);
    if (filteredTargetsSegDiff.size() < numComponents) {
        return new MixtureMultivariateNormalFitResult(null, false, false);
    }
    // Assume that Apache Commons wants data points in the first dimension.
    // Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API.
    final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2];
    // Convert the filtered targets into 2d array (even if second dimension is length 1).  The second dimension is
    //  uncorrelated Gaussian.  This is only to get around funny API in Apache Commons, which will throw an
    //  exception if the length of the second dimension is < 2
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
    final NormalDistribution nd = new NormalDistribution(rng, 0, .1);
    for (int i = 0; i < filteredTargetsSegDiff.size(); i++) {
        filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i);
        filteredTargetsSegDiff2d[i][1] = nd.sample();
    }
    final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents);
    final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d);
    try {
        multivariateNormalMixtureExpectationMaximization.fit(estimateEM0);
    } catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) {
        //  did not converge.  Include the model as it was when the exception was thrown.
        return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true);
    }
    return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true);
}