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

Example 51 with MaxEval

use of org.apache.commons.math3.optim.MaxEval in project tetrad by cmu-phil.

the class GeneralizedSemEstimator method optimize.

private double[] optimize(MultivariateFunction function, double[] values, int optimizer) {
    PointValuePair pair;
    if (optimizer == 1) {
        // 0.01, 0.000001
        // 2.0D * FastMath.ulp(1.0D), 1e-8
        MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7);
        pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000));
    } else if (optimizer == 2) {
        MultivariateOptimizer search = new SimplexOptimizer(1e-7, 1e-7);
        pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000), new NelderMeadSimplex(values.length));
    } else if (optimizer == 3) {
        int dim = values.length;
        int additionalInterpolationPoints = 0;
        final int numIterpolationPoints = 2 * dim + 1 + additionalInterpolationPoints;
        BOBYQAOptimizer search = new BOBYQAOptimizer(numIterpolationPoints);
        pair = search.optimize(new MaxEval(100000), new ObjectiveFunction(function), GoalType.MINIMIZE, new InitialGuess(values), SimpleBounds.unbounded(dim));
    } else if (optimizer == 4) {
        MultivariateOptimizer search = new CMAESOptimizer(3000000, .05, false, 0, 0, new MersenneTwister(), false, new SimplePointChecker<PointValuePair>(0.5, 0.5));
        pair = search.optimize(new MaxEval(30000), new ObjectiveFunction(function), GoalType.MINIMIZE, new InitialGuess(values), new CMAESOptimizer.Sigma(new double[values.length]), new CMAESOptimizer.PopulationSize(1000));
    } else if (optimizer == 5) {
        // 0.01, 0.000001
        // 2.0D * FastMath.ulp(1.0D), 1e-8
        MultivariateOptimizer search = new PowellOptimizer(.05, .05);
        pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000));
    } else if (optimizer == 6) {
        MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7);
        pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MAXIMIZE, new MaxEval(10000));
    } else {
        throw new IllegalStateException();
    }
    return pair.getPoint();
}

Example 52 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, SimplexOptimizer optimiser, MultivariateFunction func, double qvalue) {
    try {
        final NelderMeadSimplex simplex = new NelderMeadSimplex(1);
        final double[] initialSolution = { qvalue };
        final PointValuePair solution = optimiser.optimize(new MaxEval(1000), new InitialGuess(initialSolution), simplex, new ObjectiveFunction(func), GoalType.MINIMIZE);
        final UnivariatePointValuePair next = (solution == null) ? null : new UnivariatePointValuePair(solution.getPointRef()[0], solution.getValue());
        if (next == null) {
            return current;
        }
        if (current != null) {
            return (next.getValue() < current.getValue()) ? next : current;
        }
        return next;
    } catch (final Exception ex) {
        return current;
    }
}

Example 53 with MaxEval

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

the class EmGainAnalysis method fit.

/**
 * Fit the EM-gain distribution (Gaussian * Gamma).
 *
 * @param histogram The distribution
 */
private void fit(int[] histogram) {
    final int[] limits = limits(histogram);
    final double[] x = getX(limits);
    final double[] y = getY(histogram, limits);
    Plot plot = new Plot(TITLE, "ADU", "Frequency");
    double yMax = MathUtils.max(y);
    plot.setLimits(limits[0], limits[1], 0, yMax);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    ImageJUtils.display(TITLE, plot);
    // Estimate remaining parameters.
    // Assuming a gamma_distribution(shape,scale) then mean = shape * scale
    // scale = gain
    // shape = Photons = mean / gain
    double mean = getMean(histogram) - settings.bias;
    // Note: if the bias is too high then the mean will be negative. Just move the bias.
    while (mean < 0) {
        settings.bias -= 1;
        mean += 1;
    }
    double photons = mean / settings.gain;
    if (settings.settingSimulate) {
        ImageJUtils.log("Simulated bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.settingBias, MathUtils.rounded(settings.settingGain), MathUtils.rounded(settings.settingNoise), MathUtils.rounded(settings.settingPhotons));
    }
    ImageJUtils.log("Estimate bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    final int max = (int) x[x.length - 1];
    double[] pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot.setColor(Color.blue);
    plot.addPoints(x, pg, Plot.LINE);
    ImageJUtils.display(TITLE, plot);
    // Perform a fit
    final CustomPowellOptimizer o = new CustomPowellOptimizer(1e-6, 1e-16, 1e-6, 1e-16);
    final double[] startPoint = new double[] { photons, settings.gain, settings.noise, settings.bias };
    int maxEval = 3000;
    final String[] paramNames = { "Photons", "Gain", "Noise", "Bias" };
    // Set bounds
    final double[] lower = new double[] { 0, 0.5 * settings.gain, 0, settings.bias - settings.noise };
    final double[] upper = new double[] { 2 * photons, 2 * settings.gain, settings.gain, settings.bias + settings.noise };
    // Restart until converged.
    // TODO - Maybe fix this with a better optimiser. This needs to be tested on real data.
    PointValuePair solution = null;
    for (int iter = 0; iter < 3; iter++) {
        IJ.showStatus("Fitting histogram ... Iteration " + iter);
        try {
            // Basic Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(fun), GoalType.MINIMIZE, new InitialGuess((solution == null) ? startPoint : solution.getPointRef()));
            if (solution == null || optimum.getValue() < solution.getValue()) {
                final double[] point = optimum.getPointRef();
                // Check the bounds
                for (int i = 0; i < point.length; i++) {
                    if (point[i] < lower[i] || point[i] > upper[i]) {
                        throw new ComputationException(String.format("Fit out of of estimated range: %s %f", paramNames[i], point[i]));
                    }
                }
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
        try {
            // Bounded Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(fun, lower, upper);
            PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded((solution == null) ? startPoint : solution.getPointRef())));
            final double[] point = adapter.unboundedToBounded(optimum.getPointRef());
            optimum = new PointValuePair(point, optimum.getValue());
            if (solution == null || optimum.getValue() < solution.getValue()) {
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Bounded Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
    }
    ImageJUtils.finished();
    if (solution == null) {
        ImageJUtils.log("Failed to fit the distribution");
        return;
    }
    final double[] point = solution.getPointRef();
    photons = point[0];
    settings.gain = point[1];
    settings.noise = point[2];
    settings.bias = (int) Math.round(point[3]);
    final String label = String.format("Fitted bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    ImageJUtils.log(label);
    if (settings.settingSimulate) {
        ImageJUtils.log("Relative Error bias=%s, gain=%s, noise=%s, photons=%s", MathUtils.rounded(relativeError(settings.bias, settings.settingBias)), MathUtils.rounded(relativeError(settings.gain, settings.settingGain)), MathUtils.rounded(relativeError(settings.noise, settings.settingNoise)), MathUtils.rounded(relativeError(photons, settings.settingPhotons)));
    }
    // Show the PoissonGammaGaussian approximation
    double[] approxValues = null;
    if (settings.showApproximation) {
        approxValues = new double[x.length];
        final PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / settings.gain, settings.noise);
        final double expected = photons * settings.gain;
        for (int i = 0; i < approxValues.length; i++) {
            approxValues[i] = fun.likelihood(x[i] - settings.bias, expected);
        }
        yMax = MathUtils.maxDefault(yMax, approxValues);
    }
    // Replot
    pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot = new Plot(TITLE, "ADU", "Frequency");
    plot.setLimits(limits[0], limits[1], 0, yMax * 1.05);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    plot.setColor(Color.red);
    plot.addPoints(x, pg, Plot.LINE);
    plot.addLabel(0, 0, label);
    if (settings.showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, approxValues, Plot.LINE);
    }
    ImageJUtils.display(TITLE, plot);
}

Example 54 with MaxEval

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

the class EmGainAnalysis method getFunction.

private static MultivariateFunction getFunction(final int[] limits, final double[] y, final int max, final int maxEval) {
    return new MultivariateFunction() {

        int eval;

        @Override
        public double value(double[] point) {
            IJ.showProgress(++eval, maxEval);
            if (ImageJUtils.isInterrupted()) {
                throw new TooManyEvaluationsException(maxEval);
            }
            // Compute the sum of squares between the two functions
            final double photons = point[0];
            final double gain = point[1];
            final double noise = point[2];
            final int bias = (int) Math.round(point[3]);
            final double[] g = pdf(max, photons, gain, noise, bias);
            double ss = 0;
            for (int c = limits[0]; c <= limits[1]; c++) {
                final double d = g[c] - y[c];
                ss += d * d;
            }
            return ss;
        }
    };
}