use of uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.
the class MaximumLikelihoodFitter method computeFisherInformationMatrix.
@Override
protected FisherInformationMatrix computeFisherInformationMatrix(double[] y, double[] a) {
final int n = y.length;
final LikelihoodWrapper maximumLikelihoodFunction = createLikelihoodWrapper((NonLinearFunction) function, n, y, a);
final double[] solution = getInitialSolution(a);
// Since the type of function is unknown this cannot be done here.
return new FisherInformationMatrix(maximumLikelihoodFunction.fisherInformation(solution));
}
use of uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix 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;
}
use of uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.
the class WLseLvmSteppingFunctionSolver method computeFunctionFisherInformationMatrix.
@Override
protected FisherInformationMatrix computeFunctionFisherInformationMatrix(double[] y, double[] a) {
// Compute using the scaled Hessian as per the above method.
// Use a dedicated procedure that omits computing beta.
final WPoissonGradientProcedure p = WPoissonGradientProcedureUtils.create(y, getWeights(y.length), (Gradient1Function) function);
p.computeFisherInformation(a);
if (p.isNaNGradients()) {
throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
}
return new FisherInformationMatrix(p.getLinear(), p.numberOfGradients);
}
use of uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.
the class ApacheLvmFitter method computeFit.
@Override
public FitStatus computeFit(double[] y, final double[] fx, double[] a, double[] parametersVariance) {
final int n = y.length;
try {
// Different convergence thresholds seem to have no effect on the resulting fit, only the
// number of
// iterations for convergence
final double initialStepBoundFactor = 100;
final double costRelativeTolerance = 1e-10;
final double parRelativeTolerance = 1e-10;
final double orthoTolerance = 1e-10;
final double threshold = Precision.SAFE_MIN;
// Extract the parameters to be fitted
final double[] initialSolution = getInitialSolution(a);
// TODO - Pass in more advanced stopping criteria.
// Create the target and weight arrays
final double[] yd = new double[n];
// final double[] w = new double[n];
for (int i = 0; i < n; i++) {
yd[i] = y[i];
// w[i] = 1;
}
final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
// @formatter:off
final LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd);
if (function instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) function).canComputeValuesAndJacobian()) {
// Compute together, or each individually
builder.model(new ValueAndJacobianFunction() {
final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) function;
@Override
public Pair<RealVector, RealMatrix> value(RealVector point) {
final double[] p = point.toArray();
final org.apache.commons.lang3.tuple.Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
return new Pair<>(new ArrayRealVector(result.getKey(), false), new Array2DRowRealMatrix(result.getValue(), false));
}
@Override
public RealVector computeValue(double[] params) {
return new ArrayRealVector(fun.computeValues(params), false);
}
@Override
public RealMatrix computeJacobian(double[] params) {
return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
}
});
} else {
// Compute separately
builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) function, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) function, a, n));
}
final LeastSquaresProblem problem = builder.build();
final Optimum optimum = optimizer.optimize(problem);
final double[] parameters = optimum.getPoint().toArray();
setSolution(a, parameters);
iterations = optimum.getIterations();
evaluations = optimum.getEvaluations();
if (parametersVariance != null) {
// Set up the Jacobian.
final RealMatrix j = optimum.getJacobian();
// Compute transpose(J)J.
final RealMatrix jTj = j.transpose().multiply(j);
final double[][] data = (jTj instanceof Array2DRowRealMatrix) ? ((Array2DRowRealMatrix) jTj).getDataRef() : jTj.getData();
final FisherInformationMatrix m = new FisherInformationMatrix(data);
setDeviations(parametersVariance, m);
}
// Compute function value
if (fx != null) {
final ValueFunction function = (ValueFunction) this.function;
function.initialise0(a);
function.forEach(new ValueProcedure() {
int index;
@Override
public void execute(double value) {
fx[index++] = value;
}
});
}
// As this is unweighted then we can do this to get the sum of squared residuals
// This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
// just computes the dot product anyway.
value = optimum.getResiduals().dotProduct(optimum.getResiduals());
} catch (final TooManyEvaluationsException ex) {
return FitStatus.TOO_MANY_EVALUATIONS;
} catch (final TooManyIterationsException ex) {
return FitStatus.TOO_MANY_ITERATIONS;
} 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) {
// TODO - Find out the other exceptions from the fitter and add return values to match.
return FitStatus.UNKNOWN;
}
return FitStatus.OK;
}
use of uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix in project GDSC-SMLM by aherbert.
the class FastMleSteppingFunctionSolver method computeLastFisherInformationMatrix.
@Override
protected FisherInformationMatrix computeLastFisherInformationMatrix(double[] fx) {
Gradient2Function f2 = (Gradient2Function) function;
// Capture the y-values if necessary
if (fx != null && fx.length == f2.size()) {
f2 = new Gradient2FunctionValueStore(f2, fx);
}
// Add the weights if necessary
if (obsVariances != null) {
f2 = OffsetGradient2Function.wrapGradient2Function(f2, obsVariances);
}
// The fisher information is that for a Poisson process
final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(f2);
initialiseAndRun(p);
if (p.isNaNGradients()) {
throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
}
return new FisherInformationMatrix(p.getLinear(), p.numberOfGradients);
}
Aggregations