Examples with BoundedNonLinearConjugateGradientOptimizer org.apache.commons.math3.optim.nonlinear.scalar.gradient.BoundedNonLinearConjugateGradientOptimizer used on opensource projects

Example 1 with BoundedNonLinearConjugateGradientOptimizer

use of org.apache.commons.math3.optim.nonlinear.scalar.gradient.BoundedNonLinearConjugateGradientOptimizer in project GDSC-SMLM by aherbert.

the class MaximumLikelihoodFitter method computeFit.

/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.fitting.nonlinear.BaseFunctionSolver#computeFit(double[], double[], double[], double[])
	 */
public FitStatus computeFit(double[] y, double[] y_fit, double[] a, double[] a_dev) {
    final int n = y.length;
    LikelihoodWrapper maximumLikelihoodFunction = createLikelihoodWrapper((NonLinearFunction) f, n, y, a);
    @SuppressWarnings("rawtypes") BaseOptimizer baseOptimiser = null;
    try {
        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
            // This is as per the method in Numerical Recipes 10.5 (Direction Set (Powell's) method)
            // I could extend the optimiser and 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 = FastMath.sqrt(relativeThreshold);
            //final double lineAbs = FastMath.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;
            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
                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)));
                double[] solution = adapter.unboundedToBounded(optimum.getPointRef());
                optimum = new PointValuePair(solution, optimum.getValue());
            } else {
                if (powellFunction == null) {
                    // Python code by using the sqrt of the number of photons and background.
                    if (mapGaussian) {
                        Gaussian2DFunction gf = (Gaussian2DFunction) f;
                        // Re-map signal and background using the sqrt
                        int[] indices = gf.gradientIndices();
                        int[] map = new int[indices.length];
                        int count = 0;
                        // Background is always first
                        if (indices[0] == Gaussian2DFunction.BACKGROUND) {
                            map[count++] = 0;
                        }
                        // Look for the Signal in multiple peak 2D Gaussians
                        for (int i = 1; i < indices.length; i++) if (indices[i] % 6 == Gaussian2DFunction.SIGNAL) {
                            map[count++] = i;
                        }
                        if (count > 0) {
                            powellFunction = new MappedMultivariateLikelihood(maximumLikelihoodFunction, Arrays.copyOf(map, count));
                        }
                    }
                    if (powellFunction == null) {
                        powellFunction = new MultivariateLikelihood(maximumLikelihoodFunction);
                    }
                }
                // Update the maximum likelihood function in the Powell function wrapper
                powellFunction.fun = maximumLikelihoodFunction;
                OptimizationData positionChecker = null;
                // new org.apache.commons.math3.optim.PositionChecker(relativeThreshold, absoluteThreshold);
                SimpleBounds simpleBounds = null;
                if (powellFunction.isMapped()) {
                    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);
                    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.
            int numberOfInterpolationPoints = this.getNumberOfFittedParameters() + 2;
            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.
            // CMAESOptimiser based on Matlab code:
            // https://www.lri.fr/~hansen/cmaes.m
            // Take the defaults from the Matlab documentation
            //Double.NEGATIVE_INFINITY;
            double stopFitness = 0;
            boolean isActiveCMA = true;
            int diagonalOnly = 0;
            int checkFeasableCount = 1;
            RandomGenerator random = new Well19937c();
            boolean generateStatistics = false;
            // The sigma determines the search range for the variables. It should be 1/3 of the initial search region.
            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 = FastMath.min(sigma.length * sigma.length * 30, getMaxEvaluations() / 2);
            evaluations = 0;
            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
                    data[0] = new InitialGuess(optimum.getPointRef());
                    popSize *= 2;
                    data[1] = new CMAESOptimizer.PopulationSize(popSize);
                }
                CMAESOptimizer o = new CMAESOptimizer(getMaxIterations(), stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, new SimpleValueChecker(relativeThreshold, absoluteThreshold));
                baseOptimiser = o;
                PointValuePair result = o.optimize(data);
                iterations += o.getIterations();
                evaluations += o.getEvaluations();
                //		o.getEvaluations(), totalEvaluations);
                if (optimum == null || result.getValue() < optimum.getValue()) {
                    optimum = result;
                }
            }
            // Prevent incrementing the iterations again
            baseOptimiser = null;
        } else if (searchMethod == SearchMethod.BFGS) {
            // BFGS can use an approximate line search minimisation where as Powell and conjugate gradient
            // methods require a more accurate line minimisation. The BFGS search does not do a full 
            // minimisation but takes appropriate steps in the direction of the current gradient.
            // Do not use the convergence checker on the value of the function. Use the convergence on the 
            // point coordinate and gradient
            //BFGSOptimizer o = new BFGSOptimizer(new SimpleValueChecker(rel, abs));
            BFGSOptimizer o = new BFGSOptimizer();
            baseOptimiser = o;
            // Configure maximum step length for each dimension using the bounds
            double[] stepLength = new double[lower.length];
            for (int i = 0; i < stepLength.length; i++) {
                stepLength[i] = (upper[i] - lower[i]) * 0.3333333;
                if (stepLength[i] <= 0)
                    stepLength[i] = Double.POSITIVE_INFINITY;
            }
            // The GoalType is always minimise so no need to pass this in
            OptimizationData positionChecker = null;
            //new org.apache.commons.math3.optim.PositionChecker(relativeThreshold, absoluteThreshold);
            optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunctionGradient(new MultivariateVectorLikelihood(maximumLikelihoodFunction)), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), new InitialGuess(startPoint), new SimpleBounds(lowerConstraint, upperConstraint), new BFGSOptimizer.GradientTolerance(relativeThreshold), positionChecker, new BFGSOptimizer.StepLength(stepLength));
        } 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.
            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 l = 0;
            for (double d : gradient) l += d * d;
            final double bracketingStep = FastMath.min(0.001, ((l > 1) ? 1.0 / l : 1));
            //System.out.printf("Bracketing step = %f (length=%f)\n", bracketingStep, l);
            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));
        //maximumLikelihoodFunction.value(solution, gradient);
        //System.out.printf("Iter = %d, %g @ %s : %s\n", iterations, ll, Arrays.toString(solution),
        //		Arrays.toString(gradient));
        }
        final double[] solution = optimum.getPointRef();
        setSolution(a, solution);
        if (a_dev != null) {
            // Assume the Maximum Likelihood estimator returns the optimum fit (achieves the Cramer Roa
            // lower bounds) and so the covariance can be obtained from the Fisher Information Matrix.
            FisherInformationMatrix m = new FisherInformationMatrix(maximumLikelihoodFunction.fisherInformation(a));
            setDeviations(a_dev, m.crlb(true));
        }
        // Reverse negative log likelihood for maximum likelihood score
        value = -optimum.getValue();
    } catch (TooManyIterationsException e) {
        //e.printStackTrace();
        return FitStatus.TOO_MANY_ITERATIONS;
    } catch (TooManyEvaluationsException e) {
        //e.printStackTrace();
        return FitStatus.TOO_MANY_EVALUATIONS;
    } catch (ConvergenceException e) {
        //System.out.printf("Singular non linear model = %s\n", e.getMessage());
        return FitStatus.SINGULAR_NON_LINEAR_MODEL;
    } catch (BFGSOptimizer.LineSearchRoundoffException e) {
        //e.printStackTrace();
        return FitStatus.FAILED_TO_CONVERGE;
    } catch (Exception e) {
        //System.out.printf("Unknown error = %s\n", e.getMessage());
        e.printStackTrace();
        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 2 with BoundedNonLinearConjugateGradientOptimizer

use of org.apache.commons.math3.optim.nonlinear.scalar.gradient.BoundedNonLinearConjugateGradientOptimizer 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 3 with BoundedNonLinearConjugateGradientOptimizer

use of org.apache.commons.math3.optim.nonlinear.scalar.gradient.BoundedNonLinearConjugateGradientOptimizer 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 4 with BoundedNonLinearConjugateGradientOptimizer

use of org.apache.commons.math3.optim.nonlinear.scalar.gradient.BoundedNonLinearConjugateGradientOptimizer in project GDSC-SMLM by aherbert.

the class Image3DAligner method align.

/**
 * Align the image with the reference with sub-pixel accuracy. Compute the translation required to
 * move the target image onto the reference image for maximum correlation.
 *
 * @param target the target
 * @param refinements the maximum number of refinements for sub-pixel accuracy
 * @param error the error for sub-pixel accuracy (i.e. stop when improvements are less than this
 *        error)
 * @return [x,y,z,value]
 * @throws IllegalArgumentException If any dimension is less than 2, or if larger than the
 *         initialised reference
 */
private double[] align(DhtData target, int refinements, double error) {
    // Multiply by the reference. This allows the reference to be shared across threads.
    final DoubleDht3D correlation = target.dht.conjugateMultiply(reference.dht, buffer);
    // Store for reuse
    buffer = correlation.getData();
    correlation.inverseTransform();
    correlation.swapOctants();
    // Normalise:
    // ( Σ xiyi - nx̄ӯ ) / ( (Σ xi^2 - nx̄^2) (Σ yi^2 - nӯ^2) )^0.5
    // 
    // (sumXy - sumX*sumY/n) / sqrt( (sumXx - sumX^2 / n) * (sumYy - sumY^2 / n) )
    // Only do this over the range where at least half the original images overlap,
    // i.e. the insert point of one will be the middle of the other when shifted.
    int ix = Math.min(reference.ix, target.ix);
    int iy = Math.min(reference.iy, target.iy);
    int iz = Math.min(reference.iz, target.iz);
    int ixw = Math.max(reference.ix + reference.width, target.ix + target.width);
    int iyh = Math.max(reference.iy + reference.height, target.iy + target.height);
    int izd = Math.max(reference.iz + reference.depth, target.iz + target.depth);
    if (minimumDimensionOverlap > 0) {
        final double f = (1 - minimumDimensionOverlap) / 2;
        final int ux = (int) (Math.round(Math.min(reference.width, target.width) * f));
        final int uy = (int) (Math.round(Math.min(reference.height, target.height) * f));
        final int uz = (int) (Math.round(Math.min(reference.depth, target.depth) * f));
        ix += ux;
        ixw -= ux;
        iy += uy;
        iyh -= uy;
        iz += uz;
        izd -= uz;
    }
    cropDimensions = new int[] { ix, iy, iz, ixw - ix, iyh - iy, izd - iz };
    // The maximum correlation unnormalised. Since this is unnormalised
    // it will be biased towards the centre of the image. This is used
    // to restrict the bounds for finding the maximum of the normalised correlation
    // which should be close to this.
    int maxi = correlation.findMaxIndex(ix, iy, iz, cropDimensions[3], cropDimensions[4], cropDimensions[5]);
    // Check in the spatial domain
    checkCorrelation(target, correlation, maxi);
    // Compute sum from rolling sum using:
    // sum(x,y,z,w,h,d) =
    // + s(x+w-1,y+h-1,z+d-1)
    // - s(x-1,y+h-1,z+d-1)
    // - s(x+w-1,y-1,z+d-1)
    // + s(x-1,y-1,z+d-1)
    // /* Image above must be subtracted so reverse sign*/
    // - s(x+w-1,y+h-1,z-1)
    // + s(x-1,y+h-1,z-1)
    // + s(x+w-1,y-1,z-1)
    // - s(x-1,y-1,z-1)
    // Note:
    // s(i,j,k) = 0 when either i,j,k < 0
    // i = imax when i>imax
    // j = jmax when j>jmax
    // k = kmax when k>kmax
    // Note: The correlation is for the movement of the reference over the target
    final int nc_2 = nc / 2;
    final int nr_2 = nr / 2;
    final int ns_2 = ns / 2;
    final int[] centre = new int[] { nc_2, nr_2, ns_2 };
    // Compute the shift from the centre
    final int dx = nc_2 - ix;
    final int dy = nr_2 - iy;
    final int dz = ns_2 - iz;
    // For the reference (moved -dx,-dy,-dz over the target)
    int rx = -dx;
    int ry = -dy;
    int rz = -dz;
    // For the target (moved dx,dy,dz over the reference)
    int tx = dx;
    int ty = dy;
    int tz = dz;
    // Precompute the x-1,x+w-1,y-1,y+h-1
    final int nx = cropDimensions[3];
    final int[] rx1 = new int[nx];
    final int[] rxw1 = new int[nx];
    final int[] tx1 = new int[nx];
    final int[] txw1 = new int[nx];
    final int[] width = new int[nx];
    for (int c = ix, i = 0; c < ixw; c++, i++) {
        rx1[i] = Math.max(-1, rx - 1);
        rxw1[i] = Math.min(nc, rx + nc) - 1;
        rx++;
        tx1[i] = Math.max(-1, tx - 1);
        txw1[i] = Math.min(nc, tx + nc) - 1;
        tx--;
        width[i] = rxw1[i] - rx1[i];
    }
    final int ny = cropDimensions[4];
    final int[] ry1 = new int[ny];
    final int[] ryh1 = new int[ny];
    final int[] ty1 = new int[ny];
    final int[] tyh1 = new int[ny];
    final int[] h = new int[ny];
    for (int r = iy, j = 0; r < iyh; r++, j++) {
        ry1[j] = Math.max(-1, ry - 1);
        ryh1[j] = Math.min(nr, ry + nr) - 1;
        ry++;
        ty1[j] = Math.max(-1, ty - 1);
        tyh1[j] = Math.min(nr, ty + nr) - 1;
        ty--;
        h[j] = ryh1[j] - ry1[j];
    }
    final double[] rs = reference.sum;
    final double[] rss = reference.sumSq;
    final double[] ts = target.sum;
    final double[] tss = target.sumSq;
    final double[] rsum = new double[2];
    final double[] tsum = new double[2];
    final int size = Math.min(reference.size, target.size);
    final int minimumN = (int) (Math.round(size * minimumOverlap));
    int maxj = -1;
    double max = 0;
    for (int s = iz; s < izd; s++) {
        // Compute the z-1,z+d-1
        final int rz_1 = Math.max(-1, rz - 1);
        final int rz_d_1 = Math.min(ns, rz + ns) - 1;
        rz++;
        final int tz_1 = Math.max(-1, tz - 1);
        final int tz_d_1 = Math.min(ns, tz + ns) - 1;
        tz--;
        final int d = rz_d_1 - rz_1;
        for (int r = iy, j = 0; r < iyh; r++, j++) {
            final int base = s * nrByNc + r * nc;
            final int hd = h[j] * d;
            for (int c = ix, i = 0; c < ixw; c++, i++) {
                final double sumXy = buffer[base + c];
                compute(rx1[i], ry1[j], rz_1, rxw1[i], ryh1[j], rz_d_1, width[i], h[j], d, rs, rss, rsum);
                compute(tx1[i], ty1[j], tz_1, txw1[i], tyh1[j], tz_d_1, width[i], h[j], d, ts, tss, tsum);
                // Compute the correlation
                // (sumXy - sumX*sumY/n) / sqrt( (sumXx - sumX^2 / n) * (sumYy - sumY^2 / n) )
                final int n = width[i] * hd;
                final double numerator = sumXy - (rsum[X] * tsum[Y] / n);
                final double denominator1 = rsum[XX] - (rsum[X] * rsum[X] / n);
                final double denominator2 = tsum[YY] - (tsum[Y] * tsum[Y] / n);
                double corr;
                if (denominator1 == 0 || denominator2 == 0) {
                    // If there is data and all the variances are the same then correlation is perfect
                    if (rsum[XX] == tsum[YY] && rsum[XX] == sumXy && rsum[XX] > 0) {
                        corr = 1;
                    } else {
                        corr = 0;
                    }
                } else {
                    // Leave as raw for debugging, i.e. do not clip to range [-1:1]
                    corr = numerator / Math.sqrt(denominator1 * denominator2);
                }
                buffer[base + c] = corr;
                if (n < minimumN) {
                    continue;
                }
                // Check normalisation with some margin for error
                if (corr > 1.0001) {
                    // It is likely to occur at the bounds.
                    continue;
                }
                if (corr > max) {
                    max = corr;
                    maxj = base + c;
                } else if (corr == max) {
                    // Get shift from centre
                    final int[] xyz1 = correlation.getXyz(maxj);
                    final int[] xyz2 = correlation.getXyz(base + c);
                    int d1 = 0;
                    int d2 = 0;
                    for (int k = 0; k < 3; k++) {
                        d1 += MathUtils.pow2(xyz1[k] - centre[k]);
                        d2 += MathUtils.pow2(xyz2[k] - centre[k]);
                    }
                    if (d2 < d1) {
                        max = corr;
                        maxj = base + c;
                    }
                }
            }
        }
    }
    // The maximum correlation with normalisation
    // correlation.findMaxIndex(ix, iy, iz, iw - ix, ih - iy, id - iz);
    maxi = maxj;
    final int[] xyz = correlation.getXyz(maxi);
    // Report the shift required to move from the centre of the target image to the reference
    // @formatter:off
    final double[] result = new double[] { nc_2 - xyz[0], nr_2 - xyz[1], ns_2 - xyz[2], buffer[maxi] };
    if (refinements > 0) {
        // Create a cubic spline using a small region of pixels around the maximum
        if (calc == null) {
            calc = new CubicSplineCalculator();
        }
        // Avoid out-of-bounds errors. Only use the range that was normalised
        final int x = MathUtils.clip(ix, ixw - 4, xyz[0] - 1);
        final int y = MathUtils.clip(iy, iyh - 4, xyz[1] - 1);
        final int z = MathUtils.clip(iz, izd - 4, xyz[2] - 1);
        final DoubleImage3D crop = correlation.crop(x, y, z, 4, 4, 4, region);
        region = crop.getData();
        final CustomTricubicFunction f = CustomTricubicFunctionUtils.create(calc.compute(region));
        // Find the maximum starting at the current origin
        final int ox = xyz[0] - x;
        final int oy = xyz[1] - y;
        final int oz = xyz[2] - z;
        // Scale to the cubic spline dimensions of 0-1
        final double[] origin = new double[] { ox / 3.0, oy / 3.0, oz / 3.0 };
        // Simple condensing search
        if (searchMode == SearchMode.BINARY) {
            // Can this use the current origin as a start point?
            // Currently we evaluate 8-cube vertices. A better search
            // would evaluate 27 points around the optimum, pick the best then condense
            // the range.
            final double[] optimum = f.search(true, refinements, relativeThreshold, -1);
            final double value = optimum[3];
            if (value > result[3]) {
                result[3] = value;
                // Convert the maximum back with scaling
                for (int i = 0; i < 3; i++) {
                    result[i] -= (optimum[i] - origin[i]) * 3.0;
                }
                return result;
            }
        } else {
            // Gradient search
            try {
                final SplineFunction sf = new SplineFunction(f, origin);
                final BoundedNonLinearConjugateGradientOptimizer optimiser = new BoundedNonLinearConjugateGradientOptimizer(BoundedNonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES, // set the number of refinements
                new SimpleValueChecker(relativeThreshold, -1, refinements));
                final PointValuePair opt = optimiser.optimize(maxEvaluations, bounds, GoalType.MINIMIZE, new InitialGuess(origin), // Scale the error for the position check
                new PositionChecker(-1, error / 3.0), new ObjectiveFunction(sf::value), new ObjectiveFunctionGradient(point -> {
                    // This must be new each time
                    final double[] partialDerivative1 = new double[3];
                    sf.value(point, partialDerivative1);
                    return partialDerivative1;
                }));
                // Check it is higher. Invert since we did a minimisation.
                final double value = -opt.getValue();
                if (value > result[3]) {
                    result[3] = value;
                    // Convert the maximum back with scaling
                    final double[] optimum = opt.getPointRef();
                    for (int i = 0; i < 3; i++) {
                        result[i] -= (optimum[i] - origin[i]) * 3.0;
                    }
                    return result;
                }
            } catch (final Exception ex) {
            // Ignore this
            }
        }
    }
    return result;
}