Example usage for org.apache.commons.math3.optim.nonlinear.vector.jacobian LevenbergMarquardtOptimizer LevenbergMarquardtOptimizer

Introduction

In this page you can find the example usage for org.apache.commons.math3.optim.nonlinear.vector.jacobian LevenbergMarquardtOptimizer LevenbergMarquardtOptimizer.

Prototype

public LevenbergMarquardtOptimizer(double initialStepBoundFactor, double costRelativeTolerance,
        double parRelativeTolerance, double orthoTolerance, double threshold) 

Source Link

Document

The arguments control the behaviour of the default convergence checking procedure.

Usage

From source file:gdsc.smlm.fitting.nonlinear.ApacheLVMFitter.java

public FitStatus fit(int n, double[] y, double[] y_fit, double[] a, double[] a_dev, double[] error,
        double noise) {
    numberOfFittedPoints = n;//from  w ww  . j  a va 2  s  . c o m

    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);
        PointVectorValuePair optimum = optimizer.optimize(new MaxIter(getMaxEvaluations()),
                new MaxEval(Integer.MAX_VALUE),
                new ModelFunctionJacobian(new MultivariateMatrixFunctionWrapper(f, a, n)),
                new ModelFunction(new MultivariateVectorFunctionWrapper(f, a, n)), new Target(yd),
                new Weight(w), new InitialGuess(initialSolution));

        final double[] parameters = optimum.getPoint();
        setSolution(a, parameters);
        iterations = optimizer.getIterations();
        evaluations = optimizer.getEvaluations();
        if (a_dev != null) {
            double[][] covar = optimizer.computeCovariances(parameters, threshold);
            setDeviations(a_dev, covar);
        }
        // Compute fitted function if desired 
        if (y_fit != null) {
            f.initialise(a);
            for (int i = 0; i < n; i++)
                y_fit[i] = f.eval(i);
        }

        residualSumOfSquares = error[0] = optimizer.getChiSquare();
        totalSumOfSquares = getSumOfSquares(n, y);
    } catch (TooManyEvaluationsException e) {
        return FitStatus.FAILED_TO_CONVERGE;
    } 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;
}

---

From source file:gdsc.smlm.ij.plugins.BlinkEstimator.java

/**
 * Fit the dark time to counts of molecules curve. Only use the first n fitted points.
 * <p>/*ww w .  ja  va 2s .c  o  m*/
 * Calculates:<br/>
 * N = The number of photoblinking molecules in the sample<br/>
 * nBlink = The average number of blinks per flourophore<br/>
 * tOff = The off-time
 * 
 * @param td
 *            The dark time
 * @param ntd
 *            The counts of molecules
 * @param nFittedPoints
 * @param log
 *            Write the fitting results to the ImageJ log window
 * @return The fitted parameters [N, nBlink, tOff], or null if no fit was possible
 */
public double[] fit(double[] td, double[] ntd, int nFittedPoints, boolean log) {
    blinkingModel = new BlinkingFunction();
    blinkingModel.setLogging(true);
    for (int i = 0; i < nFittedPoints; i++)
        blinkingModel.addPoint(td[i], ntd[i]);

    // Different convergence thresholds seem to have no effect on the resulting fit, only the number of
    // iterations for convergence
    double initialStepBoundFactor = 100;
    double costRelativeTolerance = 1e-6;
    double parRelativeTolerance = 1e-6;
    double orthoTolerance = 1e-6;
    double threshold = Precision.SAFE_MIN;
    LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor,
            costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
    try {
        double[] obs = blinkingModel.getY();

        PointVectorValuePair optimum = optimizer.optimize(new MaxIter(1000), new MaxEval(Integer.MAX_VALUE),
                new ModelFunctionJacobian(new MultivariateMatrixFunction() {
                    public double[][] value(double[] point) throws IllegalArgumentException {
                        return blinkingModel.jacobian(point);
                    }
                }), new ModelFunction(blinkingModel), new Target(obs), new Weight(blinkingModel.getWeights()),
                new InitialGuess(new double[] { ntd[0], 0.1, td[1] }));

        blinkingModel.setLogging(false);

        double[] parameters = optimum.getPoint();

        double[] exp = optimum.getValue();
        double mean = 0;
        for (double d : obs)
            mean += d;
        mean /= obs.length;
        double ssResiduals = 0, ssTotal = 0;
        for (int i = 0; i < obs.length; i++) {
            ssResiduals += (obs[i] - exp[i]) * (obs[i] - exp[i]);
            ssTotal += (obs[i] - mean) * (obs[i] - mean);
        }

        r2 = 1 - ssResiduals / ssTotal;
        adjustedR2 = getAdjustedCoefficientOfDetermination(ssResiduals, ssTotal, obs.length, parameters.length);

        if (log) {
            Utils.log("  Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, Utils.rounded(r2, 4),
                    Utils.rounded(adjustedR2, 4));
            Utils.log("  N=%s, nBlink=%s, tOff=%s (%s frames)", Utils.rounded(parameters[0], 4),
                    Utils.rounded(parameters[1], 4), Utils.rounded(parameters[2], 4),
                    Utils.rounded(parameters[2] / msPerFrame, 4));
        }

        return parameters;
    } catch (TooManyIterationsException e) {
        if (log)
            Utils.log("  Failed to fit %d points: Too many iterations (%d)", blinkingModel.size(),
                    optimizer.getIterations());
        return null;
    } catch (ConvergenceException e) {
        if (log)
            Utils.log("  Failed to fit %d points", blinkingModel.size());
        return null;
    }
}

---