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

Example 31 with MaxEval

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

the class JumpDistanceAnalysis method doFitJumpDistancesMle.

/**
 * Fit the jump distances using a maximum likelihood estimation with the given number of species.
 *
 * <p>Results are sorted by the diffusion coefficient ascending.
 *
 * @param jumpDistances The jump distances (in um^2)
 * @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.
 */
@Nullable
private double[][] doFitJumpDistancesMle(double[] jumpDistances, double estimatedD, int n) {
    final MaxEval maxEval = new MaxEval(20000);
    final CustomPowellOptimizer powellOptimizer = createCustomPowellOptimizer();
    calibrated = isCalibrated();
    if (n == 1) {
        try {
            final JumpDistanceFunction function = new JumpDistanceFunction(jumpDistances, estimatedD);
            // The Powell algorithm can use more general bounds: 0 - Infinity
            final SimpleBounds bounds = new SimpleBounds(function.getLowerBounds(), function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY));
            final PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), bounds, new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
            final double[] fitParams = solution.getPointRef();
            ll = solution.getValue();
            lastFitValue = fitValue = MathUtils.getAkaikeInformationCriterion(ll, 1);
            final double[] coefficients = fitParams;
            final double[] fractions = new double[] { 1 };
            LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=1) : %s, MLE = %s, Akaike IC = %s (%d evaluations)", formatD(fitParams[0]), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), powellOptimizer.getEvaluations());
            return new double[][] { coefficients, fractions };
        } catch (final TooManyEvaluationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many evaluation (%d)", powellOptimizer.getEvaluations());
        } catch (final TooManyIterationsException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many iterations (%d)", powellOptimizer.getIterations());
        } catch (final ConvergenceException ex) {
            LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : %s", ex.getMessage());
        }
        return null;
    }
    final MixedJumpDistanceFunction function = new MixedJumpDistanceFunction(jumpDistances, estimatedD, n);
    final double[] lB = function.getLowerBounds();
    int evaluations = 0;
    PointValuePair constrainedSolution = null;
    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.MAXIMIZE);
        evaluations = powellOptimizer.getEvaluations();
        LoggerUtils.log(logger, Level.FINE, "Powell optimiser fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
    } catch (final TooManyEvaluationsException ex) {
        LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many evaluation (%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);
        // Try a bounded CMAES optimiser since the Powell optimiser appears to be
        // sensitive to the order of the parameters. It is not good when the fast particle
        // is the minority fraction. Could this be due to too low an upper bound?
        // 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.MAXIMIZE, 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) : MLE = %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.MAXIMIZE, 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) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
                }
            } catch (final TooManyEvaluationsException ex) {
            // No solution
            }
        }
        if (constrainedSolution != null) {
            try {
                // Re-optimise with Powell?
                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.MAXIMIZE);
                if (solution.getValue() > constrainedSolution.getValue()) {
                    evaluations = cmaesOptimizer.getEvaluations();
                    constrainedSolution = solution;
                    LoggerUtils.log(logger, Level.INFO, "Powell optimiser re-fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
                }
            } 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;
    }
    final double[] fitParams = constrainedSolution.getPointRef();
    ll = constrainedSolution.getValue();
    // Since the fractions must sum to one we subtract 1 degree of freedom from the number of
    // parameters
    fitValue = MathUtils.getAkaikeInformationCriterion(ll, 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), MLE = %s, Akaike IC = %s (%d evaluations)", n, formatD(d), format(f), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), evaluations);
    if (isValid(d, f)) {
        lastFitValue = fitValue;
        return new double[][] { coefficients, fractions };
    }
    return null;
}

Example 32 with MaxEval

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

the class PoissonGammaGaussianFisherInformation method findMaximum.

/**
 * Find the maximum of the integrated function A^2/P. P is the Poisson-Gamma convolution, A is the
 * partial gradient.
 *
 * <p>When both A and P are convolved with a Gaussian kernel, the integral of this function - 1 is
 * the Fisher information.
 *
 * <p>This method is used to determine the maximum of the function. The integral of A should equal
 * 1. The range can be determined using a fraction of the maximum, or integrating until the sum is
 * 1.
 *
 * @param theta the Poisson mean
 * @param rel Relative threshold
 * @return [max,max value,evaluations]
 */
public double[] findMaximum(final double theta, double rel) {
    if (theta < MIN_MEAN) {
        throw new IllegalArgumentException("Poisson mean must be positive");
    }
    final double dirac = PoissonGammaFunction.dirac(theta);
    final UnivariateFunction f = new UnivariateFunction() {

        double[] gradient = new double[1];

        @Override
        public double value(double x) {
            // Return F without the dirac. This is so that the maximum of the function
            // is known for relative height analysis.
            double value;
            if (x == 0) {
                gradient[0] = dirac / gain;
                value = gradient[0] * theta;
            } else {
                value = PoissonGammaFunction.unscaledPoissonGammaPartial(x, theta, gain, gradient);
            }
            return getF(value, gradient[0]);
        }
    };
    // Determine the extent of the function: A^2/P
    // Without convolution this will have a maximum that is above, and
    // converges to [Poisson mean * amplification].
    final double mean = theta * gain;
    final int iterations = 10;
    final QuickConvergenceChecker checker = new QuickConvergenceChecker(iterations);
    final BrentOptimizer opt = new BrentOptimizer(rel, Double.MIN_VALUE, checker);
    UnivariatePointValuePair pair;
    try {
        final double start = (theta < 1) ? 0 : mean;
        pair = opt.optimize(new SearchInterval(0, mean * 1.5, start), GoalType.MAXIMIZE, new MaxEval(50000), new UnivariateObjectiveFunction(f));
    } catch (final TooManyEvaluationsException ex) {
        // This should not occur with the quick checker fixed iterations
        // - just get the current best
        pair = checker.getBest();
    }
    final double[] max = new double[3];
    max[0] = pair.getPoint();
    max[1] = pair.getValue();
    max[2] = opt.getEvaluations();
    return max;
}

Example 33 with MaxEval

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

the class Fire method findMin.

private static PointValuePair findMin(PointValuePair current, SimplexOptimizer optimiser, MultivariateFunction func, double[] initialSolution) {
    try {
        final NelderMeadSimplex simplex = new NelderMeadSimplex(initialSolution.length);
        final PointValuePair next = optimiser.optimize(new MaxEval(1000), new InitialGuess(initialSolution), simplex, new ObjectiveFunction(func), GoalType.MINIMIZE);
        if (next == null) {
            return current;
        }
        if (current != null) {
            return (next.getValue() < current.getValue()) ? next : current;
        }
        return next;
    } catch (final Exception ex) {
        return current;
    }
}

Example 34 with MaxEval

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

the class Fire method findMin.

private static UnivariatePointValuePair findMin(UnivariatePointValuePair current, UnivariateOptimizer optimiser, UnivariateFunction func, double qvalue, double factor) {
    try {
        final BracketFinder bracket = new BracketFinder();
        bracket.search(func, GoalType.MINIMIZE, qvalue * factor, qvalue / factor);
        final UnivariatePointValuePair next = optimiser.optimize(GoalType.MINIMIZE, new MaxEval(3000), new SearchInterval(bracket.getLo(), bracket.getHi(), bracket.getMid()), new UnivariateObjectiveFunction(func));
        if (next == null) {
            return current;
        }
        if (current != null) {
            return (next.getValue() < current.getValue()) ? next : current;
        }
        return next;
    } catch (final Exception ex) {
        return current;
    }
}

Example 35 with MaxEval

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

the class BoundedNonLinearConjugateGradientOptimizer method doOptimize.

@Override
protected PointValuePair doOptimize() {
    final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
    final double[] point = getStartPoint();
    final GoalType goal = getGoalType();
    final int n = point.length;
    sign = (goal == GoalType.MINIMIZE) ? -1 : 1;
    double[] unbounded = point.clone();
    applyBounds(point);
    double[] gradient = computeObjectiveGradient(point);
    checkGradients(gradient, unbounded);
    if (goal == GoalType.MINIMIZE) {
        for (int i = 0; i < n; i++) {
            gradient[i] = -gradient[i];
        }
    }
    // Initial search direction.
    double[] steepestDescent = preconditioner.precondition(point, gradient);
    double[] searchDirection = steepestDescent.clone();
    double delta = 0;
    for (int i = 0; i < n; ++i) {
        delta += gradient[i] * searchDirection[i];
    }
    // Used for non-gradient based line search
    LineSearch line = null;
    double rel = 1e-6;
    double abs = 1e-10;
    if (getConvergenceChecker() instanceof SimpleValueChecker) {
        rel = ((SimpleValueChecker) getConvergenceChecker()).getRelativeThreshold();
        abs = ((SimpleValueChecker) getConvergenceChecker()).getRelativeThreshold();
    }
    line = new LineSearch(Math.sqrt(rel), Math.sqrt(abs));
    PointValuePair current = null;
    int maxEval = getMaxEvaluations();
    for (; ; ) {
        incrementIterationCount();
        final double objective = computeObjectiveValue(point);
        final PointValuePair previous = current;
        current = new PointValuePair(point, objective);
        if (previous != null && checker.converged(getIterations(), previous, current)) {
            // We have found an optimum.
            return current;
        }
        double step;
        if (useGradientLineSearch) {
            // Classic code using the gradient function for the line search:
            // Find the optimal step in the search direction.
            final UnivariateFunction lsf = new LineSearchFunction(point, searchDirection);
            final double uB;
            try {
                uB = findUpperBound(lsf, 0, initialStep);
                // Check if the bracket found a minimum. Otherwise just move to the new point.
                if (noBracket) {
                    step = uB;
                } else {
                    // XXX Last parameters is set to a value close to zero in order to
                    // work around the divergence problem in the "testCircleFitting"
                    // unit test (see MATH-439).
                    // System.out.printf("Bracket %f - %f - %f\n", 0., 1e-15, uB);
                    step = solver.solve(maxEval, lsf, 0, uB, 1e-15);
                    // Subtract used up evaluations.
                    maxEval -= solver.getEvaluations();
                }
            } catch (final MathIllegalStateException ex) {
                // System.out.printf("Failed to bracket %s @ %s\n", Arrays.toString(point),
                // Arrays.toString(searchDirection));
                // Line search without gradient (as per Powell optimiser)
                final UnivariatePointValuePair optimum = line.search(point, searchDirection);
                step = optimum.getPoint();
            // throw ex;
            }
        } else {
            // Line search without gradient (as per Powell optimiser)
            final UnivariatePointValuePair optimum = line.search(point, searchDirection);
            step = optimum.getPoint();
        }
        // System.out.printf("Step = %f x %s\n", step, Arrays.toString(searchDirection));
        for (int i = 0; i < point.length; ++i) {
            point[i] += step * searchDirection[i];
        }
        unbounded = point.clone();
        applyBounds(point);
        gradient = computeObjectiveGradient(point);
        checkGradients(gradient, unbounded);
        if (goal == GoalType.MINIMIZE) {
            for (int i = 0; i < n; ++i) {
                gradient[i] = -gradient[i];
            }
        }
        // Compute beta.
        final double deltaOld = delta;
        final double[] newSteepestDescent = preconditioner.precondition(point, gradient);
        delta = 0;
        for (int i = 0; i < n; ++i) {
            delta += gradient[i] * newSteepestDescent[i];
        }
        if (delta == 0) {
            return new PointValuePair(point, computeObjectiveValue(point));
        }
        final double beta;
        switch(updateFormula) {
            case FLETCHER_REEVES:
                beta = delta / deltaOld;
                break;
            case POLAK_RIBIERE:
                double deltaMid = 0;
                for (int i = 0; i < gradient.length; ++i) {
                    deltaMid += gradient[i] * steepestDescent[i];
                }
                beta = (delta - deltaMid) / deltaOld;
                break;
            default:
                // Should never happen.
                throw new MathInternalError();
        }
        steepestDescent = newSteepestDescent;
        // Compute conjugate search direction.
        if (getIterations() % n == 0 || beta < 0) {
            // Break conjugation: reset search direction.
            searchDirection = steepestDescent.clone();
        } else {
            // Compute new conjugate search direction.
            for (int i = 0; i < n; ++i) {
                searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
            }
        }
        // The gradient has already been adjusted for the search direction
        checkGradients(searchDirection, unbounded, -sign);
    }
}