Examples with CustomPowellOptimizer uk.ac.sussex.gdsc.smlm.math3.optim.nonlinear.scalar.noderiv.CustomPowellOptimizer used on opensource projects

Example 6 with CustomPowellOptimizer

use of uk.ac.sussex.gdsc.smlm.math3.optim.nonlinear.scalar.noderiv.CustomPowellOptimizer 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;
}

Example 7 with CustomPowellOptimizer

use of uk.ac.sussex.gdsc.smlm.math3.optim.nonlinear.scalar.noderiv.CustomPowellOptimizer 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 8 with CustomPowellOptimizer

use of uk.ac.sussex.gdsc.smlm.math3.optim.nonlinear.scalar.noderiv.CustomPowellOptimizer 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);
}