use of org.apache.commons.math3.optim.OptimizationData 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
LevenbergMarquardtOptimizer lvmOptimizer = new LevenbergMarquardtOptimizer();
try {
final JumpDistanceCumulFunction function = new JumpDistanceCumulFunction(jdHistogram[0], jdHistogram[1], estimatedD);
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(function.guess()).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return function.jacobian(point);
}
}).build();
//@formatter:on
Optimum lvmSolution = lvmOptimizer.optimize(problem);
double[] fitParams = lvmSolution.getPoint().toArray();
// True for an unweighted fit
ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
//ss = calculateSumOfSquares(function.getY(), function.value(fitParams));
lastIC = ic = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.x.length, 1);
double[] coefficients = fitParams;
double[] fractions = new double[] { 1 };
logger.info("Fit Jump distance (N=1) : %s, SS = %s, IC = %s (%d evaluations)", formatD(fitParams[0]), Maths.rounded(ss, 4), Maths.rounded(ic, 4), lvmSolution.getEvaluations());
return new double[][] { coefficients, fractions };
} catch (TooManyIterationsException e) {
logger.info("LVM optimiser failed to fit (N=1) : Too many iterations : %s", e.getMessage());
} catch (ConvergenceException e) {
logger.info("LVM optimiser failed to fit (N=1) : %s", e.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.
MixedJumpDistanceCumulFunctionMultivariate function = new MixedJumpDistanceCumulFunctionMultivariate(jdHistogram[0], jdHistogram[1], estimatedD, n);
double[] lB = function.getLowerBounds();
int evaluations = 0;
PointValuePair constrainedSolution = null;
MaxEval maxEval = new MaxEval(20000);
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();
logger.debug("Powell optimiser fit (N=%d) : SS = %f (%d evaluations)", n, constrainedSolution.getValue(), evaluations);
} catch (TooManyEvaluationsException e) {
logger.info("Powell optimiser failed to fit (N=%d) : Too many evaluations (%d)", n, powellOptimizer.getEvaluations());
} catch (TooManyIterationsException e) {
logger.info("Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
} catch (ConvergenceException e) {
logger.info("Powell optimiser failed to fit (N=%d) : %s", n, e.getMessage());
}
if (constrainedSolution == null) {
logger.info("Trying CMAES optimiser with restarts ...");
double[] uB = function.getUpperBounds();
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.
double[] s = new double[lB.length];
for (int i = 0; i < s.length; i++) s[i] = (uB[i] - lB[i]) / 3;
OptimizationData sigma = new CMAESOptimizer.Sigma(s);
OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
// Iterate this for stability in the initial guess
CMAESOptimizer cmaesOptimizer = createCMAESOptimizer();
for (int i = 0; i <= fitRestarts; i++) {
// Try from the initial guess
try {
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;
logger.debug("CMAES optimiser [%da] fit (N=%d) : SS = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
}
} catch (TooManyEvaluationsException e) {
}
if (constrainedSolution == null)
continue;
// Try from the current optimum
try {
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;
logger.debug("CMAES optimiser [%db] fit (N=%d) : SS = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
}
} catch (TooManyEvaluationsException e) {
}
}
if (constrainedSolution != null) {
// Re-optimise with Powell?
try {
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;
logger.info("Powell optimiser re-fit (N=%d) : SS = %f (%d evaluations)", n, constrainedSolution.getValue(), evaluations);
}
} catch (TooManyEvaluationsException e) {
} catch (TooManyIterationsException e) {
} catch (ConvergenceException e) {
}
}
}
if (constrainedSolution == null) {
logger.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;
LevenbergMarquardtOptimizer lvmOptimizer = new LevenbergMarquardtOptimizer();
try {
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(fitParams).target(functionGradient.getY()).weight(new DiagonalMatrix(functionGradient.getWeights())).model(functionGradient, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return functionGradient.jacobian(point);
}
}).build();
//@formatter:on
lvmSolution = lvmOptimizer.optimize(problem);
// True for an unweighted fit
double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
// All fitted parameters must be above zero
if (ss < this.ss && Maths.min(lvmSolution.getPoint().toArray()) > 0) {
logger.info(" Re-fitting improved the SS from %s to %s (-%s%%)", Maths.rounded(this.ss, 4), Maths.rounded(ss, 4), Maths.rounded(100 * (this.ss - ss) / this.ss, 4));
fitParams = lvmSolution.getPoint().toArray();
this.ss = ss;
evaluations += lvmSolution.getEvaluations();
}
} catch (TooManyIterationsException e) {
logger.error("Failed to re-fit : Too many iterations : %s", e.getMessage());
} catch (ConvergenceException e) {
logger.error("Failed to re-fit : %s", e.getMessage());
}
// Since the fractions must sum to one we subtract 1 degree of freedom from the number of parameters
ic = Maths.getAkaikeInformationCriterionFromResiduals(ss, function.x.length, fitParams.length - 1);
double[] d = new double[n];
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);
double[] coefficients = d;
double[] fractions = f;
logger.info("Fit Jump distance (N=%d) : %s (%s), SS = %s, IC = %s (%d evaluations)", n, formatD(d), format(f), Maths.rounded(ss, 4), Maths.rounded(ic, 4), evaluations);
if (isValid(d, f)) {
lastIC = ic;
return new double[][] { coefficients, fractions };
}
return null;
}
use of org.apache.commons.math3.optim.OptimizationData 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.
*/
private double[][] doFitJumpDistancesMLE(double[] jumpDistances, double estimatedD, int n) {
MaxEval maxEval = new MaxEval(20000);
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
SimpleBounds bounds = new SimpleBounds(function.getLowerBounds(), function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY));
PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), bounds, new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
double[] fitParams = solution.getPointRef();
ll = solution.getValue();
lastIC = ic = Maths.getAkaikeInformationCriterion(ll, jumpDistances.length, 1);
double[] coefficients = fitParams;
double[] fractions = new double[] { 1 };
logger.info("Fit Jump distance (N=1) : %s, MLE = %s, IC = %s (%d evaluations)", formatD(fitParams[0]), Maths.rounded(ll, 4), Maths.rounded(ic, 4), powellOptimizer.getEvaluations());
return new double[][] { coefficients, fractions };
} catch (TooManyEvaluationsException e) {
logger.info("Powell optimiser failed to fit (N=1) : Too many evaluation (%d)", powellOptimizer.getEvaluations());
} catch (TooManyIterationsException e) {
logger.info("Powell optimiser failed to fit (N=1) : Too many iterations (%d)", powellOptimizer.getIterations());
} catch (ConvergenceException e) {
logger.info("Powell optimiser failed to fit (N=1) : %s", e.getMessage());
}
return null;
}
MixedJumpDistanceFunction function = new MixedJumpDistanceFunction(jumpDistances, estimatedD, n);
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();
logger.debug("Powell optimiser fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
} catch (TooManyEvaluationsException e) {
logger.info("Powell optimiser failed to fit (N=%d) : Too many evaluation (%d)", n, powellOptimizer.getEvaluations());
} catch (TooManyIterationsException e) {
logger.info("Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
} catch (ConvergenceException e) {
logger.info("Powell optimiser failed to fit (N=%d) : %s", n, e.getMessage());
}
if (constrainedSolution == null) {
logger.info("Trying CMAES optimiser with restarts ...");
double[] uB = function.getUpperBounds();
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.
double[] s = new double[lB.length];
for (int i = 0; i < s.length; i++) s[i] = (uB[i] - lB[i]) / 3;
OptimizationData sigma = new CMAESOptimizer.Sigma(s);
OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
// Iterate this for stability in the initial guess
CMAESOptimizer cmaesOptimizer = createCMAESOptimizer();
for (int i = 0; i <= fitRestarts; i++) {
// Try from the initial guess
try {
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;
logger.debug("CMAES optimiser [%da] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
}
} catch (TooManyEvaluationsException e) {
}
if (constrainedSolution == null)
continue;
// Try from the current optimum
try {
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;
logger.debug("CMAES optimiser [%db] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
}
} catch (TooManyEvaluationsException e) {
}
}
if (constrainedSolution != null) {
try {
// Re-optimise with Powell?
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;
logger.info("Powell optimiser re-fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
}
} catch (TooManyEvaluationsException e) {
} catch (TooManyIterationsException e) {
} catch (ConvergenceException e) {
}
}
}
if (constrainedSolution == null) {
logger.info("Failed to fit N=%d", n);
return null;
}
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
ic = Maths.getAkaikeInformationCriterion(ll, jumpDistances.length, fitParams.length - 1);
double[] d = new double[n];
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);
double[] coefficients = d;
double[] fractions = f;
logger.info("Fit Jump distance (N=%d) : %s (%s), MLE = %s, IC = %s (%d evaluations)", n, formatD(d), format(f), Maths.rounded(ll, 4), Maths.rounded(ic, 4), evaluations);
if (isValid(d, f)) {
lastIC = ic;
return new double[][] { coefficients, fractions };
}
return null;
}
use of org.apache.commons.math3.optim.OptimizationData in project GDSC-SMLM by aherbert.
the class BinomialFitter method fitBinomial.
/**
* Fit the binomial distribution (n,p) to the cumulative histogram. Performs fitting assuming a fixed n value and
* attempts to optimise p.
*
* @param histogram
* The input histogram
* @param mean
* The histogram mean (used to estimate p). Calculated if NaN.
* @param n
* The n to evaluate
* @param zeroTruncated
* True if the model should ignore n=0 (zero-truncated binomial)
* @return The best fit (n, p)
* @throws IllegalArgumentException
* If any of the input data values are negative
* @throws IllegalArgumentException
* If any fitting a zero truncated binomial and there are no values above zero
*/
public PointValuePair fitBinomial(double[] histogram, double mean, int n, boolean zeroTruncated) {
if (Double.isNaN(mean))
mean = getMean(histogram);
if (zeroTruncated && histogram[0] > 0) {
log("Fitting zero-truncated histogram but there are zero values - Renormalising to ignore zero");
double cumul = 0;
for (int i = 1; i < histogram.length; i++) cumul += histogram[i];
if (cumul == 0)
throw new IllegalArgumentException("Fitting zero-truncated histogram but there are no non-zero values");
histogram[0] = 0;
for (int i = 1; i < histogram.length; i++) histogram[i] /= cumul;
}
int nFittedPoints = Math.min(histogram.length, n + 1) - ((zeroTruncated) ? 1 : 0);
if (nFittedPoints < 1) {
log("No points to fit (%d): Histogram.length = %d, n = %d, zero-truncated = %b", nFittedPoints, histogram.length, n, zeroTruncated);
return null;
}
// The model is only fitting the probability p
// For a binomial n*p = mean => p = mean/n
double[] initialSolution = new double[] { FastMath.min(mean / n, 1) };
// Create the function
BinomialModelFunction function = new BinomialModelFunction(histogram, n, zeroTruncated);
double[] lB = new double[1];
double[] uB = new double[] { 1 };
SimpleBounds bounds = new SimpleBounds(lB, uB);
// Fit
// CMAESOptimizer or BOBYQAOptimizer support bounds
// CMAESOptimiser based on Matlab code:
// https://www.lri.fr/~hansen/cmaes.m
// Take the defaults from the Matlab documentation
int maxIterations = 2000;
//Double.NEGATIVE_INFINITY;
double stopFitness = 0;
boolean isActiveCMA = true;
int diagonalOnly = 0;
int checkFeasableCount = 1;
RandomGenerator random = new Well19937c();
boolean generateStatistics = false;
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.
OptimizationData sigma = new CMAESOptimizer.Sigma(new double[] { (uB[0] - lB[0]) / 3 });
OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(2))));
try {
PointValuePair solution = null;
boolean noRefit = maximumLikelihood;
if (n == 1 && zeroTruncated) {
// No need to fit
solution = new PointValuePair(new double[] { 1 }, 0);
noRefit = true;
} else {
GoalType goalType = (maximumLikelihood) ? GoalType.MAXIMIZE : GoalType.MINIMIZE;
// Iteratively fit
CMAESOptimizer opt = new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
for (int iteration = 0; iteration <= fitRestarts; iteration++) {
try {
// Start from the initial solution
PointValuePair result = opt.optimize(new InitialGuess(initialSolution), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2));
// opt.getEvaluations());
if (solution == null || result.getValue() < solution.getValue()) {
solution = result;
}
} catch (TooManyEvaluationsException e) {
} catch (TooManyIterationsException e) {
}
if (solution == null)
continue;
try {
// Also restart from the current optimum
PointValuePair result = opt.optimize(new InitialGuess(solution.getPointRef()), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2));
// opt.getEvaluations());
if (result.getValue() < solution.getValue()) {
solution = result;
}
} catch (TooManyEvaluationsException e) {
} catch (TooManyIterationsException e) {
}
}
if (solution == null)
return null;
}
if (noRefit) {
// Although we fit the log-likelihood, return the sum-of-squares to allow
// comparison across different n
double p = solution.getPointRef()[0];
double ss = 0;
double[] obs = function.p;
double[] exp = function.getP(p);
for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
return new PointValuePair(solution.getPointRef(), ss);
} else // We can do a LVM refit if the number of fitted points is more than 1
if (nFittedPoints > 1) {
// Improve SS fit with a gradient based LVM optimizer
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
try {
final BinomialModelFunctionGradient gradientFunction = new BinomialModelFunctionGradient(histogram, n, zeroTruncated);
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(solution.getPointRef()).target(gradientFunction.p).weight(new DiagonalMatrix(gradientFunction.getWeights())).model(gradientFunction, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return gradientFunction.jacobian(point);
}
}).build();
//@formatter:on
Optimum lvmSolution = optimizer.optimize(problem);
// Check the pValue is valid since the LVM is not bounded.
double p = lvmSolution.getPoint().getEntry(0);
if (p <= 1 && p >= 0) {
// True if the weights are 1
double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
// ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
if (ss < solution.getValue()) {
// Utils.rounded(100 * (solution.getValue() - ss) / solution.getValue(), 4));
return new PointValuePair(lvmSolution.getPoint().toArray(), ss);
}
}
} catch (TooManyIterationsException e) {
log("Failed to re-fit: Too many iterations: %s", e.getMessage());
} catch (ConvergenceException e) {
log("Failed to re-fit: %s", e.getMessage());
} catch (Exception e) {
// Ignore this ...
}
}
return solution;
} catch (Exception e) {
log("Failed to fit Binomial distribution with N=%d : %s", n, e.getMessage());
}
return null;
}
use of org.apache.commons.math3.optim.OptimizationData 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;
}
use of org.apache.commons.math3.optim.OptimizationData in project GDSC-SMLM by aherbert.
the class PCPALMFitting method runBoundedOptimiser.
private PointValuePair runBoundedOptimiser(double[][] gr, double[] initialSolution, double[] lB, double[] uB, SumOfSquaresModelFunction function) {
// Create the functions to optimise
ObjectiveFunction objective = new ObjectiveFunction(new SumOfSquaresMultivariateFunction(function));
ObjectiveFunctionGradient gradient = new ObjectiveFunctionGradient(new SumOfSquaresMultivariateVectorFunction(function));
final boolean debug = false;
// Try a BFGS optimiser since this will produce a deterministic solution and can respect bounds.
PointValuePair optimum = null;
boundedEvaluations = 0;
final MaxEval maxEvaluations = new MaxEval(2000);
MultivariateOptimizer opt = null;
for (int iteration = 0; iteration <= fitRestarts; iteration++) {
try {
opt = new BFGSOptimizer();
final double relativeThreshold = 1e-6;
// Configure maximum step length for each dimension using the bounds
double[] stepLength = new double[lB.length];
for (int i = 0; i < stepLength.length; i++) stepLength[i] = (uB[i] - lB[i]) * 0.3333333;
// The GoalType is always minimise so no need to pass this in
optimum = opt.optimize(maxEvaluations, gradient, objective, new InitialGuess((optimum == null) ? initialSolution : optimum.getPointRef()), new SimpleBounds(lB, uB), new BFGSOptimizer.GradientTolerance(relativeThreshold), new BFGSOptimizer.StepLength(stepLength));
if (debug)
System.out.printf("BFGS Iter %d = %g (%d)\n", iteration, optimum.getValue(), opt.getEvaluations());
} catch (TooManyEvaluationsException e) {
// No need to restart
break;
} catch (RuntimeException e) {
// No need to restart
break;
} finally {
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
//Double.NEGATIVE_INFINITY;
double stopFitness = 0;
boolean isActiveCMA = true;
int diagonalOnly = 0;
int checkFeasableCount = 1;
//Well19937c();
RandomGenerator random = new Well44497b();
boolean generateStatistics = false;
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.
double[] range = new double[lB.length];
for (int i = 0; i < lB.length; i++) range[i] = (uB[i] - lB[i]) / 3;
OptimizationData sigma = new CMAESOptimizer.Sigma(range);
OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(initialSolution.length))));
SimpleBounds bounds = new SimpleBounds(lB, uB);
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 <= fitRestarts; iteration++) {
try {
// Start from the initial solution
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 (TooManyEvaluationsException e) {
} catch (TooManyIterationsException e) {
} finally {
boundedEvaluations += maxEvaluations.getMaxEval();
}
if (optimum == null)
continue;
try {
// Also restart from the current optimum
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 (TooManyEvaluationsException e) {
} catch (TooManyIterationsException e) {
} finally {
boundedEvaluations += maxEvaluations.getMaxEval();
}
}
return optimum;
}
Aggregations