Examples with MaxEval org.apache.commons.math3.optim.MaxEval used on opensource projects

Example 26 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project GDSC-SMLM by aherbert.

the class PcPalmFitting method runBoundedOptimiser.

private PointValuePair runBoundedOptimiser(double[] initialSolution, double[] lowerB, double[] upperB, SumOfSquaresModelFunction function) {
    // Create the functions to optimise
    final ObjectiveFunction objective = new ObjectiveFunction(new SumOfSquaresMultivariateFunction(function));
    final ObjectiveFunctionGradient gradient = new ObjectiveFunctionGradient(new SumOfSquaresMultivariateVectorFunction(function));
    final boolean debug = false;
    // Try a gradient optimiser since this will produce a deterministic solution
    PointValuePair optimum = null;
    boundedEvaluations = 0;
    final MaxEval maxEvaluations = new MaxEval(2000);
    MultivariateOptimizer opt = null;
    for (int iteration = 0; iteration <= settings.fitRestarts; iteration++) {
        try {
            final double relativeThreshold = 1e-6;
            opt = new BoundedNonLinearConjugateGradientOptimizer(BoundedNonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES, new SimpleValueChecker(relativeThreshold, -1));
            optimum = opt.optimize(maxEvaluations, gradient, objective, GoalType.MINIMIZE, new InitialGuess((optimum == null) ? initialSolution : optimum.getPointRef()), new SimpleBounds(lowerB, upperB));
            if (debug) {
                System.out.printf("Bounded Iter %d = %g (%d)\n", iteration, optimum.getValue(), opt.getEvaluations());
            }
        } catch (final RuntimeException ex) {
            // No need to restart
            break;
        } finally {
            if (opt != null) {
                boundedEvaluations += opt.getEvaluations();
            }
        }
    }
    // Try a CMAES optimiser which is non-deterministic. To overcome this we perform restarts.
    // CMAESOptimiser based on Matlab code:
    // https://www.lri.fr/~hansen/cmaes.m
    // Take the defaults from the Matlab documentation
    final double stopFitness = 0;
    final boolean isActiveCma = true;
    final int diagonalOnly = 0;
    final int checkFeasableCount = 1;
    final RandomGenerator random = new RandomGeneratorAdapter(UniformRandomProviders.create());
    final boolean generateStatistics = false;
    final ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10);
    // The sigma determines the search range for the variables. It should be 1/3 of the initial
    // search region.
    final double[] range = new double[lowerB.length];
    for (int i = 0; i < lowerB.length; i++) {
        range[i] = (upperB[i] - lowerB[i]) / 3;
    }
    final OptimizationData sigma = new CMAESOptimizer.Sigma(range);
    final OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(initialSolution.length))));
    final SimpleBounds bounds = new SimpleBounds(lowerB, upperB);
    opt = new CMAESOptimizer(maxEvaluations.getMaxEval(), stopFitness, isActiveCma, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
    // Restart the optimiser several times and take the best answer.
    for (int iteration = 0; iteration <= settings.fitRestarts; iteration++) {
        try {
            // Start from the initial solution
            final PointValuePair constrainedSolution = opt.optimize(new InitialGuess(initialSolution), objective, GoalType.MINIMIZE, bounds, sigma, popSize, maxEvaluations);
            if (debug) {
                System.out.printf("CMAES Iter %d initial = %g (%d)\n", iteration, constrainedSolution.getValue(), opt.getEvaluations());
            }
            boundedEvaluations += opt.getEvaluations();
            if (optimum == null || constrainedSolution.getValue() < optimum.getValue()) {
                optimum = constrainedSolution;
            }
        } catch (final TooManyEvaluationsException | TooManyIterationsException ex) {
        // Ignore
        } finally {
            boundedEvaluations += maxEvaluations.getMaxEval();
        }
        if (optimum == null) {
            continue;
        }
        try {
            // Also restart from the current optimum
            final PointValuePair constrainedSolution = opt.optimize(new InitialGuess(optimum.getPointRef()), objective, GoalType.MINIMIZE, bounds, sigma, popSize, maxEvaluations);
            if (debug) {
                System.out.printf("CMAES Iter %d restart = %g (%d)\n", iteration, constrainedSolution.getValue(), opt.getEvaluations());
            }
            if (constrainedSolution.getValue() < optimum.getValue()) {
                optimum = constrainedSolution;
            }
        } catch (final TooManyEvaluationsException | TooManyIterationsException ex) {
        // Ignore
        } finally {
            boundedEvaluations += maxEvaluations.getMaxEval();
        }
    }
    return optimum;
}

Example 27 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project GDSC-SMLM by aherbert.

the class PcPalmMolecules method optimiseSimplex.

private double[] optimiseSimplex(float[] x, float[] y, double[] initialSolution) {
    // Simplex optimisation
    final SkewNormalMultivariateFunction sn2 = new SkewNormalMultivariateFunction(initialSolution);
    sn2.addData(x, y);
    final NelderMeadSimplex simplex = new NelderMeadSimplex(4);
    final SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
    final PointValuePair solution = opt.optimize(new MaxEval(1000), new InitialGuess(initialSolution), simplex, new ObjectiveFunction(sn2), GoalType.MINIMIZE);
    return solution.getPointRef();
}

Example 28 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project GDSC-SMLM by aherbert.

the class BinomialFitter method fitBinomial.

/**
 * Fit the binomial distribution (n,p) to the cumulative histogram. Performs fitting assuming a
 * fixed n value and attempts to optimise p.
 *
 * @param histogram The input histogram
 * @param mean The histogram mean (used to estimate p). Calculated if NaN.
 * @param n The n to evaluate
 * @param zeroTruncated True if the model should ignore n=0 (zero-truncated binomial)
 * @return The best fit (n, p)
 * @throws IllegalArgumentException If any of the input data values are negative
 * @throws IllegalArgumentException If any fitting a zero truncated binomial and there are no
 *         values above zero
 */
public PointValuePair fitBinomial(double[] histogram, double mean, int n, boolean zeroTruncated) {
    if (Double.isNaN(mean)) {
        mean = getMean(histogram);
    }
    if (zeroTruncated && histogram[0] > 0) {
        log("Fitting zero-truncated histogram but there are zero values - " + "Renormalising to ignore zero");
        double cumul = 0;
        for (int i = 1; i < histogram.length; i++) {
            cumul += histogram[i];
        }
        if (cumul == 0) {
            throw new IllegalArgumentException("Fitting zero-truncated histogram but there are no non-zero values");
        }
        histogram[0] = 0;
        for (int i = 1; i < histogram.length; i++) {
            histogram[i] /= cumul;
        }
    }
    final int nFittedPoints = Math.min(histogram.length, n + 1) - ((zeroTruncated) ? 1 : 0);
    if (nFittedPoints < 1) {
        log("No points to fit (%d): Histogram.length = %d, n = %d, zero-truncated = %b", nFittedPoints, histogram.length, n, zeroTruncated);
        return null;
    }
    // The model is only fitting the probability p
    // For a binomial n*p = mean => p = mean/n
    final double[] initialSolution = new double[] { Math.min(mean / n, 1) };
    // Create the function
    final BinomialModelFunction function = new BinomialModelFunction(histogram, n, zeroTruncated);
    final double[] lB = new double[1];
    final double[] uB = new double[] { 1 };
    final SimpleBounds bounds = new SimpleBounds(lB, uB);
    // Fit
    // CMAESOptimizer or BOBYQAOptimizer support bounds
    // CMAESOptimiser based on Matlab code:
    // https://www.lri.fr/~hansen/cmaes.m
    // Take the defaults from the Matlab documentation
    final int maxIterations = 2000;
    final double stopFitness = 0;
    final boolean isActiveCma = true;
    final int diagonalOnly = 0;
    final int checkFeasableCount = 1;
    final RandomGenerator random = new RandomGeneratorAdapter(UniformRandomProviders.create());
    final boolean generateStatistics = false;
    final ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10);
    // The sigma determines the search range for the variables. It should be 1/3 of the initial
    // search region.
    final OptimizationData sigma = new CMAESOptimizer.Sigma(new double[] { (uB[0] - lB[0]) / 3 });
    final OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(2))));
    try {
        PointValuePair solution = null;
        boolean noRefit = maximumLikelihood;
        if (n == 1 && zeroTruncated) {
            // No need to fit
            solution = new PointValuePair(new double[] { 1 }, 0);
            noRefit = true;
        } else {
            final GoalType goalType = (maximumLikelihood) ? GoalType.MAXIMIZE : GoalType.MINIMIZE;
            // Iteratively fit
            final CMAESOptimizer opt = new CMAESOptimizer(maxIterations, stopFitness, isActiveCma, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
            for (int iteration = 0; iteration <= fitRestarts; iteration++) {
                try {
                    // Start from the initial solution
                    final PointValuePair result = opt.optimize(new InitialGuess(initialSolution), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2));
                    // opt.getEvaluations());
                    if (solution == null || result.getValue() < solution.getValue()) {
                        solution = result;
                    }
                } catch (final TooManyEvaluationsException | TooManyIterationsException ex) {
                // No solution
                }
                if (solution == null) {
                    continue;
                }
                try {
                    // Also restart from the current optimum
                    final PointValuePair result = opt.optimize(new InitialGuess(solution.getPointRef()), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2));
                    // opt.getEvaluations());
                    if (result.getValue() < solution.getValue()) {
                        solution = result;
                    }
                } catch (final TooManyEvaluationsException | TooManyIterationsException ex) {
                // No solution
                }
            }
            if (solution == null) {
                return null;
            }
        }
        if (noRefit) {
            // Although we fit the log-likelihood, return the sum-of-squares to allow
            // comparison across different n
            final double p = solution.getPointRef()[0];
            double ss = 0;
            final double[] obs = function.pvalues;
            final double[] exp = function.getP(p);
            for (int i = 0; i < obs.length; i++) {
                ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
            }
            return new PointValuePair(solution.getPointRef(), ss);
        // We can do a LVM refit if the number of fitted points is more than 1.
        } else if (nFittedPoints > 1) {
            // Improve SS fit with a gradient based LVM optimizer
            final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
            try {
                final BinomialModelFunctionGradient gradientFunction = new BinomialModelFunctionGradient(histogram, n, zeroTruncated);
                // @formatter:off
                final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(solution.getPointRef()).target(gradientFunction.pvalues).weight(new DiagonalMatrix(gradientFunction.getWeights())).model(gradientFunction, gradientFunction::jacobian).build();
                // @formatter:on
                final Optimum lvmSolution = optimizer.optimize(problem);
                // Check the pValue is valid since the LVM is not bounded.
                final double p = lvmSolution.getPoint().getEntry(0);
                if (p <= 1 && p >= 0) {
                    // True if the weights are 1
                    final double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
                    // ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
                    if (ss < solution.getValue()) {
                        // MathUtils.rounded(100 * (solution.getValue() - ss) / solution.getValue(), 4));
                        return new PointValuePair(lvmSolution.getPoint().toArray(), ss);
                    }
                }
            } catch (final TooManyIterationsException ex) {
                log("Failed to re-fit: Too many iterations: %s", ex.getMessage());
            } catch (final ConvergenceException ex) {
                log("Failed to re-fit: %s", ex.getMessage());
            } catch (final Exception ex) {
            // Ignore this ...
            }
        }
        return solution;
    } catch (final RuntimeException ex) {
        log("Failed to fit Binomial distribution with N=%d : %s", n, ex.getMessage());
    }
    return null;
}

Example 29 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project GDSC-SMLM by aherbert.

the class MaximumLikelihoodFitter method computeFit.

@Override
public FitStatus computeFit(double[] y, double[] fx, double[] a, double[] parametersVariance) {
    final int n = y.length;
    final LikelihoodWrapper maximumLikelihoodFunction = createLikelihoodWrapper((NonLinearFunction) function, n, y, a);
    @SuppressWarnings("rawtypes") BaseOptimizer baseOptimiser = null;
    try {
        final double[] startPoint = getInitialSolution(a);
        PointValuePair optimum = null;
        if (searchMethod == SearchMethod.POWELL || searchMethod == SearchMethod.POWELL_BOUNDED || searchMethod == SearchMethod.POWELL_ADAPTER) {
            // Non-differentiable version using Powell Optimiser
            // Background: see Numerical Recipes 10.5 (Direction Set (Powell's) method).
            // The optimiser could be extended to implement bounds on the directions moved.
            // However the mapping adapter seems to work OK.
            final boolean basisConvergence = false;
            // Perhaps these thresholds should be tighter?
            // The default is to use the sqrt() of the overall tolerance
            // final double lineRel = Math.sqrt(relativeThreshold);
            // final double lineAbs = Math.sqrt(absoluteThreshold);
            // final double lineRel = relativeThreshold * 1e2;
            // final double lineAbs = absoluteThreshold * 1e2;
            // Since we are fitting only a small number of parameters then just use the same tolerance
            // for each search direction
            final double lineRel = relativeThreshold;
            final double lineAbs = absoluteThreshold;
            final CustomPowellOptimizer o = new CustomPowellOptimizer(relativeThreshold, absoluteThreshold, lineRel, lineAbs, null, basisConvergence);
            baseOptimiser = o;
            OptimizationData maxIterationData = null;
            if (getMaxIterations() > 0) {
                maxIterationData = new MaxIter(getMaxIterations());
            }
            if (searchMethod == SearchMethod.POWELL_ADAPTER) {
                // Try using the mapping adapter for a bounded Powell search
                final MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(new MultivariateLikelihood(maximumLikelihoodFunction), lower, upper);
                optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded(startPoint)));
                final double[] solution = adapter.unboundedToBounded(optimum.getPointRef());
                optimum = new PointValuePair(solution, optimum.getValue());
            } else {
                if (powellFunction == null) {
                    powellFunction = new MultivariateLikelihood(maximumLikelihoodFunction);
                }
                // Update the maximum likelihood function in the Powell function wrapper
                powellFunction.fun = maximumLikelihoodFunction;
                final OptimizationData positionChecker = null;
                // new org.apache.commons.math3.optim.PositionChecker(relativeThreshold,
                // absoluteThreshold);
                SimpleBounds simpleBounds = null;
                if (powellFunction.isMapped()) {
                    final MappedMultivariateLikelihood adapter = (MappedMultivariateLikelihood) powellFunction;
                    if (searchMethod == SearchMethod.POWELL_BOUNDED) {
                        simpleBounds = new SimpleBounds(adapter.map(lower), adapter.map(upper));
                    }
                    optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(powellFunction), GoalType.MINIMIZE, new InitialGuess(adapter.map(startPoint)), positionChecker, simpleBounds);
                    final double[] solution = adapter.unmap(optimum.getPointRef());
                    optimum = new PointValuePair(solution, optimum.getValue());
                } else {
                    if (searchMethod == SearchMethod.POWELL_BOUNDED) {
                        simpleBounds = new SimpleBounds(lower, upper);
                    }
                    optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(powellFunction), GoalType.MINIMIZE, new InitialGuess(startPoint), positionChecker, simpleBounds);
                }
            }
        } else if (searchMethod == SearchMethod.BOBYQA) {
            // Differentiable approximation using Powell's BOBYQA algorithm.
            // This is slower than the Powell optimiser and requires a high number of evaluations.
            final int numberOfInterpolationpoints = this.getNumberOfFittedParameters() + 2;
            final BOBYQAOptimizer o = new BOBYQAOptimizer(numberOfInterpolationpoints);
            baseOptimiser = o;
            optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new InitialGuess(startPoint), new SimpleBounds(lower, upper));
        } else if (searchMethod == SearchMethod.CMAES) {
            // TODO - Understand why the CMAES optimiser does not fit very well on test data. It appears
            // to converge too early and the likelihood scores are not as low as the other optimisers.
            // The CMAESOptimiser is based on Matlab code:
            // https://www.lri.fr/~hansen/cmaes.m
            // Take the defaults from the Matlab documentation
            final double stopFitness = 0;
            final boolean isActiveCma = true;
            final int diagonalOnly = 0;
            final int checkFeasableCount = 1;
            final RandomGenerator random = new RandomGeneratorAdapter(UniformRandomProviders.create());
            final boolean generateStatistics = false;
            // The sigma determines the search range for the variables. It should be 1/3 of the initial
            // search region.
            final double[] sigma = new double[lower.length];
            for (int i = 0; i < sigma.length; i++) {
                sigma[i] = (upper[i] - lower[i]) / 3;
            }
            int popSize = (int) (4 + Math.floor(3 * Math.log(sigma.length)));
            // The CMAES optimiser is random and restarting can overcome problems with quick
            // convergence.
            // The Apache commons documentations states that convergence should occur between 30N and
            // 300N^2
            // function evaluations
            final int n30 = Math.min(sigma.length * sigma.length * 30, getMaxEvaluations() / 2);
            evaluations = 0;
            final OptimizationData[] data = new OptimizationData[] { new InitialGuess(startPoint), new CMAESOptimizer.PopulationSize(popSize), new MaxEval(getMaxEvaluations()), new CMAESOptimizer.Sigma(sigma), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new SimpleBounds(lower, upper) };
            // Iterate to prevent early convergence
            int repeat = 0;
            while (evaluations < n30) {
                if (repeat++ > 1) {
                    // Update the start point and population size
                    if (optimum != null) {
                        data[0] = new InitialGuess(optimum.getPointRef());
                    }
                    popSize *= 2;
                    data[1] = new CMAESOptimizer.PopulationSize(popSize);
                }
                final CMAESOptimizer o = new CMAESOptimizer(getMaxIterations(), stopFitness, isActiveCma, diagonalOnly, checkFeasableCount, random, generateStatistics, new SimpleValueChecker(relativeThreshold, absoluteThreshold));
                baseOptimiser = o;
                final PointValuePair result = o.optimize(data);
                iterations += o.getIterations();
                evaluations += o.getEvaluations();
                if (optimum == null || result.getValue() < optimum.getValue()) {
                    optimum = result;
                }
            }
            // Prevent incrementing the iterations again
            baseOptimiser = null;
        } else {
            // The line search algorithm often fails. This is due to searching into a region where the
            // function evaluates to a negative so has been clipped. This means the upper bound of the
            // line cannot be found.
            // Note that running it on an easy problem (200 photons with fixed fitting (no background))
            // the algorithm does sometimes produces results better than the Powell algorithm but it is
            // slower.
            final BoundedNonLinearConjugateGradientOptimizer o = new BoundedNonLinearConjugateGradientOptimizer((searchMethod == SearchMethod.CONJUGATE_GRADIENT_FR) ? Formula.FLETCHER_REEVES : Formula.POLAK_RIBIERE, new SimpleValueChecker(relativeThreshold, absoluteThreshold));
            baseOptimiser = o;
            // Note: The gradients may become unstable at the edge of the bounds. Or they will not
            // change direction if the true solution is on the bounds since the gradient will always
            // continue towards the bounds. This is key to the conjugate gradient method. It searches
            // along a vector until the direction of the gradient is in the opposite direction (using
            // dot products, i.e. cosine of angle between them)
            // NR 10.7 states there is no advantage of the variable metric DFP or BFGS methods over
            // conjugate gradient methods. So I will try these first.
            // Try this:
            // Adapt the conjugate gradient optimiser to use the gradient to pick the search direction
            // and then for the line minimisation. However if the function is out of bounds then clip
            // the variables at the bounds and continue.
            // If the current point is at the bounds and the gradient is to continue out of bounds then
            // clip the gradient too.
            // Or: just use the gradient for the search direction then use the line minimisation/rest
            // as per the Powell optimiser. The bounds should limit the search.
            // I tried a Bounded conjugate gradient optimiser with clipped variables:
            // This sometimes works. However when the variables go a long way out of the expected range
            // the gradients can have vastly different magnitudes. This results in the algorithm
            // stalling since the gradients can be close to zero and the some of the parameters are no
            // longer adjusted. Perhaps this can be looked for and the algorithm then gives up and
            // resorts to a Powell optimiser from the current point.
            // Changed the bracketing step to very small (default is 1, changed to 0.001). This improves
            // the performance. The gradient direction is very sensitive to small changes in the
            // coordinates so a tighter bracketing of the line search helps.
            // Tried using a non-gradient method for the line search copied from the Powell optimiser:
            // This also works when the bracketing step is small but the number of iterations is higher.
            // 24.10.2014: I have tried to get conjugate gradient to work but the gradient function
            // must not behave suitably for the optimiser. In the current state both methods of using a
            // Bounded Conjugate Gradient Optimiser perform poorly relative to other optimisers:
            // Simulated : n=1000, signal=200, x=0.53, y=0.47
            // LVM : n=1000, signal=171, x=0.537, y=0.471 (1.003s)
            // Powell : n=1000, signal=187, x=0.537, y=0.48 (1.238s)
            // Gradient based PR (constrained): n=858, signal=161, x=0.533, y=0.474 (2.54s)
            // Gradient based PR (bounded): n=948, signal=161, x=0.533, y=0.473 (2.67s)
            // Non-gradient based : n=1000, signal=151.47, x=0.535, y=0.474 (1.626s)
            // The conjugate optimisers are slower, under predict the signal by the most and in the case
            // of the gradient based optimiser, fail to converge on some problems. This is worse when
            // constrained fitting is used and not tightly bounded fitting.
            // I will leave the code in as an option but would not recommend using it. I may remove it
            // in the future.
            // Note: It is strange that the non-gradient based line minimisation is more successful.
            // It may be that the gradient function is not accurate (due to round off error) or that it
            // is simply wrong when far from the optimum. My JUnit tests only evaluate the function
            // within the expected range of the answer.
            // Note the default step size on the Powell optimiser is 1 but the initial directions are
            // unit vectors.
            // So our bracketing step should be a minimum of 1 / average length of the first gradient
            // vector to prevent the first step being too large when bracketing.
            final double[] gradient = new double[startPoint.length];
            maximumLikelihoodFunction.likelihood(startPoint, gradient);
            double length = 0;
            for (final double d : gradient) {
                length += d * d;
            }
            final double bracketingStep = Math.min(0.001, ((length > 1) ? 1.0 / length : 1));
            o.setUseGradientLineSearch(gradientLineMinimisation);
            optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunctionGradient(new MultivariateVectorLikelihood(maximumLikelihoodFunction)), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new InitialGuess(startPoint), new SimpleBounds(lowerConstraint, upperConstraint), new BoundedNonLinearConjugateGradientOptimizer.BracketingStep(bracketingStep));
        }
        if (optimum == null) {
            return FitStatus.FAILED_TO_CONVERGE;
        }
        final double[] solution = optimum.getPointRef();
        setSolution(a, solution);
        if (parametersVariance != null) {
            // Compute assuming a Poisson process.
            // Note:
            // If using a Poisson-Gamma-Gaussian model then these will be incorrect.
            // However the variance for the position estimates of a 2D PSF can be
            // scaled by a factor of 2 as in Mortensen, et al (2010) Nature Methods 7, 377-383, SI 4.3.
            // Since the type of function is unknown this cannot be done here.
            final FisherInformationMatrix m = new FisherInformationMatrix(maximumLikelihoodFunction.fisherInformation(solution));
            setDeviations(parametersVariance, m);
        }
        // Reverse negative log likelihood for maximum likelihood score
        value = -optimum.getValue();
    } catch (final TooManyIterationsException ex) {
        return FitStatus.TOO_MANY_ITERATIONS;
    } catch (final TooManyEvaluationsException ex) {
        return FitStatus.TOO_MANY_EVALUATIONS;
    } catch (final ConvergenceException ex) {
        // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
        return FitStatus.SINGULAR_NON_LINEAR_MODEL;
    } catch (final Exception ex) {
        Logger.getLogger(getClass().getName()).log(Level.SEVERE, "Failed to fit", ex);
        return FitStatus.UNKNOWN;
    } finally {
        if (baseOptimiser != null) {
            iterations += baseOptimiser.getIterations();
            evaluations += baseOptimiser.getEvaluations();
        }
    }
    // Check this as likelihood functions can go wrong
    if (Double.isInfinite(value) || Double.isNaN(value)) {
        return FitStatus.INVALID_LIKELIHOOD;
    }
    return FitStatus.OK;
}

Example 30 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project GDSC-SMLM by aherbert.

the class JumpDistanceAnalysis method doFitJumpDistanceHistogram.

/**
 * Fit the jump distance histogram using a cumulative sum with the given number of species.
 *
 * <p>Results are sorted by the diffusion coefficient ascending.
 *
 * @param jdHistogram The cumulative jump distance histogram. X-axis is um^2, Y-axis is cumulative
 *        probability. Must be monototic ascending.
 * @param estimatedD The estimated diffusion coefficient
 * @param n The number of species in the mixed population
 * @return Array containing: { D (um^2), Fractions }. Can be null if no fit was made.
 */
private double[][] doFitJumpDistanceHistogram(double[][] jdHistogram, double estimatedD, int n) {
    calibrated = isCalibrated();
    if (n == 1) {
        // Fit using a single population model
        final LevenbergMarquardtOptimizer lvmOptimizer = new LevenbergMarquardtOptimizer();
        try {
            final JumpDistanceCumulFunction function = new JumpDistanceCumulFunction(jdHistogram[0], jdHistogram[1], estimatedD);
            // @formatter:off
            final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, function::jacobian).build();
            // @formatter:on
            final Optimum lvmSolution = lvmOptimizer.optimize(problem);
            final double[] fitParams = lvmSolution.getPoint().toArray();
            // True for an unweighted fit
            ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
            // ss = calculateSumOfSquares(function.getY(), function.value(fitParams));
            lastFitValue = fitValue = MathUtils.getAdjustedCoefficientOfDetermination(ss, MathUtils.getTotalSumOfSquares(function.getY()), function.x.length, 1);
            final double[] coefficients = fitParams;
            final double[] fractions = new double[] { 1 };
            LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=1) : %s, SS = %s, Adjusted R^2 = %s (%d evaluations)", formatD(fitParams[0]), MathUtils.rounded(ss, 4), MathUtils.rounded(fitValue, 4), lvmSolution.getEvaluations());
            return new double[][] { coefficients, fractions };
        } catch (final TooManyIterationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "LVM optimiser failed to fit (N=1) : Too many iterations : %s", ex.getMessage());
        } catch (final ConvergenceException ex) {
            LoggerUtils.log(logger, Level.INFO, "LVM optimiser failed to fit (N=1) : %s", ex.getMessage());
        }
    }
    // Uses a weighted sum of n exponential functions, each function models a fraction of the
    // particles.
    // An LVM fit cannot restrict the parameters so the fractions do not go below zero.
    // Use the CustomPowell/CMEASOptimizer which supports bounded fitting.
    final MixedJumpDistanceCumulFunctionMultivariate function = new MixedJumpDistanceCumulFunctionMultivariate(jdHistogram[0], jdHistogram[1], estimatedD, n);
    final double[] lB = function.getLowerBounds();
    int evaluations = 0;
    PointValuePair constrainedSolution = null;
    final MaxEval maxEval = new MaxEval(20000);
    final CustomPowellOptimizer powellOptimizer = createCustomPowellOptimizer();
    try {
        // The Powell algorithm can use more general bounds: 0 - Infinity
        constrainedSolution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MINIMIZE);
        evaluations = powellOptimizer.getEvaluations();
        LoggerUtils.log(logger, Level.FINE, "Powell optimiser fit (N=%d) : SS = %f (%d evaluations)", n, constrainedSolution.getValue(), evaluations);
    } catch (final TooManyEvaluationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many evaluations (%d)", n, powellOptimizer.getEvaluations());
    } catch (final TooManyIterationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
    } catch (final ConvergenceException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : %s", n, ex.getMessage());
    }
    if (constrainedSolution == null) {
        LoggerUtils.log(logger, Level.INFO, "Trying CMAES optimiser with restarts ...");
        final double[] uB = function.getUpperBounds();
        final SimpleBounds bounds = new SimpleBounds(lB, uB);
        // The sigma determines the search range for the variables. It should be 1/3 of the initial
        // search region.
        final double[] s = new double[lB.length];
        for (int i = 0; i < s.length; i++) {
            s[i] = (uB[i] - lB[i]) / 3;
        }
        final OptimizationData sigma = new CMAESOptimizer.Sigma(s);
        final OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
        // Iterate this for stability in the initial guess
        final CMAESOptimizer cmaesOptimizer = createCmaesOptimizer();
        for (int i = 0; i <= fitRestarts; i++) {
            // Try from the initial guess
            try {
                final PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(function.guess()), new ObjectiveFunction(function), GoalType.MINIMIZE, bounds, sigma, popSize, maxEval);
                if (constrainedSolution == null || solution.getValue() < constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%da] fit (N=%d) : SS = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
            if (constrainedSolution == null) {
                continue;
            }
            // Try from the current optimum
            try {
                final PointValuePair solution = cmaesOptimizer.optimize(new InitialGuess(constrainedSolution.getPointRef()), new ObjectiveFunction(function), GoalType.MINIMIZE, bounds, sigma, popSize, maxEval);
                if (solution.getValue() < constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%db] fit (N=%d) : SS = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
        }
        if (constrainedSolution != null) {
            // Re-optimise with Powell?
            try {
                final PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(constrainedSolution.getPointRef()), new SimpleBounds(lB, function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY)), new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MINIMIZE);
                if (solution.getValue() < constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.INFO, "Powell optimiser re-fit (N=%d) : SS = %f (%d evaluations)", n, constrainedSolution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            } catch (final TooManyIterationsException ex) {
            // No solution
            } catch (final ConvergenceException ex) {
            // No solution
            }
        }
    }
    if (constrainedSolution == null) {
        LoggerUtils.log(logger, Level.INFO, "Failed to fit N=%d", n);
        return null;
    }
    double[] fitParams = constrainedSolution.getPointRef();
    ss = constrainedSolution.getValue();
    // TODO - Try a bounded BFGS optimiser
    // Try and improve using a LVM fit
    final MixedJumpDistanceCumulFunctionGradient functionGradient = new MixedJumpDistanceCumulFunctionGradient(jdHistogram[0], jdHistogram[1], estimatedD, n);
    Optimum lvmSolution;
    final LevenbergMarquardtOptimizer lvmOptimizer = new LevenbergMarquardtOptimizer();
    try {
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(fitParams).target(functionGradient.getY()).weight(new DiagonalMatrix(functionGradient.getWeights())).model(functionGradient, functionGradient::jacobian).build();
        // @formatter:on
        lvmSolution = lvmOptimizer.optimize(problem);
        // True for an unweighted fit
        final double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
        // All fitted parameters must be above zero
        if (ss < this.ss && MathUtils.min(lvmSolution.getPoint().toArray()) > 0) {
            LoggerUtils.log(logger, Level.INFO, "  Re-fitting improved the SS from %s to %s (-%s%%)", MathUtils.rounded(this.ss, 4), MathUtils.rounded(ss, 4), MathUtils.rounded(100 * (this.ss - ss) / this.ss, 4));
            fitParams = lvmSolution.getPoint().toArray();
            this.ss = ss;
            evaluations += lvmSolution.getEvaluations();
        }
    } catch (final TooManyIterationsException ex) {
        LoggerUtils.log(logger, Level.WARNING, "Failed to re-fit : Too many iterations : %s", ex.getMessage());
    } catch (final ConvergenceException ex) {
        LoggerUtils.log(logger, Level.WARNING, "Failed to re-fit : %s", ex.getMessage());
    }
    // Since the fractions must sum to one we subtract 1 degree of freedom from the number of
    // parameters
    fitValue = MathUtils.getAdjustedCoefficientOfDetermination(ss, MathUtils.getTotalSumOfSquares(function.getY()), function.x.length, fitParams.length - 1);
    final double[] d = new double[n];
    final double[] f = new double[n];
    double sum = 0;
    for (int i = 0; i < d.length; i++) {
        f[i] = fitParams[i * 2];
        sum += f[i];
        d[i] = fitParams[i * 2 + 1];
    }
    for (int i = 0; i < f.length; i++) {
        f[i] /= sum;
    }
    // Sort by coefficient size
    sort(d, f);
    final double[] coefficients = d;
    final double[] fractions = f;
    LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=%d) : %s (%s), SS = %s, Adjusted R^2 = %s (%d evaluations)", n, formatD(d), format(f), MathUtils.rounded(ss, 4), MathUtils.rounded(fitValue, 4), evaluations);
    if (isValid(d, f)) {
        lastFitValue = fitValue;
        return new double[][] { coefficients, fractions };
    }
    return null;
}