use of org.apache.commons.math3.optim.PointValuePair 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;
}
use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.
the class FIRE method runQEstimation.
private void runQEstimation() {
IJ.showStatus(TITLE + " ...");
if (!showQEstimationInputDialog())
return;
MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
if (results == null || results.size() == 0) {
IJ.error(TITLE, "No results could be loaded");
return;
}
if (results.getCalibration() == null) {
IJ.error(TITLE, "The results are not calibrated");
return;
}
results = cropToRoi(results);
if (results.size() < 2) {
IJ.error(TITLE, "No results within the crop region");
return;
}
initialise(results, null);
// We need localisation precision.
// Build a histogram of the localisation precision.
// Get the initial mean and SD and plot as a Gaussian.
PrecisionHistogram histogram = calculatePrecisionHistogram();
if (histogram == null) {
IJ.error(TITLE, "No localisation precision available.\n \nPlease choose " + PrecisionMethod.FIXED + " and enter a precision mean and SD.");
return;
}
StoredDataStatistics precision = histogram.precision;
//String name = results.getName();
double fourierImageScale = SCALE_VALUES[imageScaleIndex];
int imageSize = IMAGE_SIZE_VALUES[imageSizeIndex];
// Create the image and compute the numerator of FRC.
// Do not use the signal so results.size() is the number of localisations.
IJ.showStatus("Computing FRC curve ...");
FireImages images = createImages(fourierImageScale, imageSize, false);
// DEBUGGING - Save the two images to disk. Load the images into the Matlab
// code that calculates the Q-estimation and make this plugin match the functionality.
//IJ.save(new ImagePlus("i1", images.ip1), "/scratch/i1.tif");
//IJ.save(new ImagePlus("i2", images.ip2), "/scratch/i2.tif");
FRC frc = new FRC();
frc.progress = progress;
frc.setFourierMethod(fourierMethod);
frc.setSamplingMethod(samplingMethod);
frc.setPerimeterSamplingFactor(perimeterSamplingFactor);
FRCCurve frcCurve = frc.calculateFrcCurve(images.ip1, images.ip2, images.nmPerPixel);
if (frcCurve == null) {
IJ.error(TITLE, "Failed to compute FRC curve");
return;
}
IJ.showStatus("Running Q-estimation ...");
// Note:
// The method implemented here is based on Matlab code provided by Bernd Rieger.
// The idea is to compute the spurious correlation component of the FRC Numerator
// using an initial estimate of distribution of the localisation precision (assumed
// to be Gaussian). This component is the contribution of repeat localisations of
// the same molecule to the numerator and is modelled as an exponential decay
// (exp_decay). The component is scaled by the Q-value which
// is the average number of times a molecule is seen in addition to the first time.
// At large spatial frequencies the scaled component should match the numerator,
// i.e. at high resolution (low FIRE number) the numerator is made up of repeat
// localisations of the same molecule and not actual structure in the image.
// The best fit is where the numerator equals the scaled component, i.e. num / (q*exp_decay) == 1.
// The FRC Numerator is plotted and Q can be determined by
// adjusting Q and the precision mean and SD to maximise the cost function.
// This can be done interactively by the user with the effect on the FRC curve
// dynamically updated and displayed.
// Compute the scaled FRC numerator
double qNorm = (1 / frcCurve.mean1 + 1 / frcCurve.mean2);
double[] frcnum = new double[frcCurve.getSize()];
for (int i = 0; i < frcnum.length; i++) {
FRCCurveResult r = frcCurve.get(i);
frcnum[i] = qNorm * r.getNumerator() / r.getNumberOfSamples();
}
// Compute the spatial frequency and the region for curve fitting
double[] q = FRC.computeQ(frcCurve, false);
int low = 0, high = q.length;
while (high > 0 && q[high - 1] > maxQ) high--;
while (low < q.length && q[low] < minQ) low++;
// Require we fit at least 10% of the curve
if (high - low < q.length * 0.1) {
IJ.error(TITLE, "Not enough points for Q estimation");
return;
}
// Obtain initial estimate of Q plateau height and decay.
// This can be done by fitting the precision histogram and then fixing the mean and sigma.
// Or it can be done by allowing the precision to be sampled and the mean and sigma
// become parameters for fitting.
// Check if we can sample precision values
boolean sampleDecay = precision != null && FIRE.sampleDecay;
double[] exp_decay;
if (sampleDecay) {
// Random sample of precision values from the distribution is used to
// construct the decay curve
int[] sample = Random.sample(10000, precision.getN(), new Well19937c());
final double four_pi2 = 4 * Math.PI * Math.PI;
double[] pre = new double[q.length];
for (int i = 1; i < q.length; i++) pre[i] = -four_pi2 * q[i] * q[i];
// Sample
final int n = sample.length;
double[] hq = new double[n];
for (int j = 0; j < n; j++) {
// Scale to SR pixels
double s2 = precision.getValue(sample[j]) / images.nmPerPixel;
s2 *= s2;
for (int i = 1; i < q.length; i++) hq[i] += FastMath.exp(pre[i] * s2);
}
for (int i = 1; i < q.length; i++) hq[i] /= n;
exp_decay = new double[q.length];
exp_decay[0] = 1;
for (int i = 1; i < q.length; i++) {
double sinc_q = sinc(Math.PI * q[i]);
exp_decay[i] = sinc_q * sinc_q * hq[i];
}
} else {
// Note: The sigma mean and std should be in the units of super-resolution
// pixels so scale to SR pixels
exp_decay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
}
// Smoothing
double[] smooth;
if (loessSmoothing) {
// Note: This computes the log then smooths it
double bandwidth = 0.1;
int robustness = 0;
double[] l = new double[exp_decay.length];
for (int i = 0; i < l.length; i++) {
// Original Matlab code computes the log for each array.
// This is equivalent to a single log on the fraction of the two.
// Perhaps the two log method is more numerically stable.
//l[i] = Math.log(Math.abs(frcnum[i])) - Math.log(exp_decay[i]);
l[i] = Math.log(Math.abs(frcnum[i] / exp_decay[i]));
}
try {
LoessInterpolator loess = new LoessInterpolator(bandwidth, robustness);
smooth = loess.smooth(q, l);
} catch (Exception e) {
IJ.error(TITLE, "LOESS smoothing failed");
return;
}
} else {
// Note: This smooths the curve before computing the log
double[] norm = new double[exp_decay.length];
for (int i = 0; i < norm.length; i++) {
norm[i] = frcnum[i] / exp_decay[i];
}
// Median window of 5 == radius of 2
MedianWindow mw = new MedianWindow(norm, 2);
smooth = new double[exp_decay.length];
for (int i = 0; i < norm.length; i++) {
smooth[i] = Math.log(Math.abs(mw.getMedian()));
mw.increment();
}
}
// Fit with quadratic to find the initial guess.
// Note: example Matlab code frc_Qcorrection7.m identifies regions of the
// smoothed log curve with low derivative and only fits those. The fit is
// used for the final estimate. Fitting a subset with low derivative is not
// implemented here since the initial estimate is subsequently optimised
// to maximise a cost function.
Quadratic curve = new Quadratic();
SimpleCurveFitter fit = SimpleCurveFitter.create(curve, new double[2]);
WeightedObservedPoints points = new WeightedObservedPoints();
for (int i = low; i < high; i++) points.add(q[i], smooth[i]);
double[] estimate = fit.fit(points.toList());
double qValue = FastMath.exp(estimate[0]);
//System.out.printf("Initial q-estimate = %s => %.3f\n", Arrays.toString(estimate), qValue);
// This could be made an option. Just use for debugging
boolean debug = false;
if (debug) {
// Plot the initial fit and the fit curve
double[] qScaled = FRC.computeQ(frcCurve, true);
double[] line = new double[q.length];
for (int i = 0; i < q.length; i++) line[i] = curve.value(q[i], estimate);
String title = TITLE + " Initial fit";
Plot2 plot = new Plot2(title, "Spatial Frequency (nm^-1)", "FRC Numerator");
String label = String.format("Q = %.3f", qValue);
plot.addPoints(qScaled, smooth, Plot.LINE);
plot.setColor(Color.red);
plot.addPoints(qScaled, line, Plot.LINE);
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
Utils.display(title, plot, Utils.NO_TO_FRONT);
}
if (fitPrecision) {
// Q - Should this be optional?
if (sampleDecay) {
// If a sample of the precision was used to construct the data for the initial fit
// then update the estimate using the fit result since it will be a better start point.
histogram.sigma = precision.getStandardDeviation();
// Normalise sum-of-squares to the SR pixel size
double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
}
// Do a multivariate fit ...
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
PointValuePair p = null;
MultiPlateauness f = new MultiPlateauness(frcnum, q, low, high);
double[] initial = new double[] { histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, qValue };
p = findMin(p, opt, f, scale(initial, 0.1));
p = findMin(p, opt, f, scale(initial, 0.5));
p = findMin(p, opt, f, initial);
p = findMin(p, opt, f, scale(initial, 2));
p = findMin(p, opt, f, scale(initial, 10));
if (p != null) {
double[] point = p.getPointRef();
histogram.mean = point[0] * images.nmPerPixel;
histogram.sigma = point[1] * images.nmPerPixel;
qValue = point[2];
}
} else {
// If so then this should be optional.
if (sampleDecay) {
if (precisionMethod != PrecisionMethod.FIXED) {
histogram.sigma = precision.getStandardDeviation();
// Normalise sum-of-squares to the SR pixel size
double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
}
exp_decay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
}
// Estimate spurious component by promoting plateauness.
// The Matlab code used random initial points for a Simplex optimiser.
// A Brent line search should be pretty deterministic so do simple repeats.
// However it will proceed downhill so if the initial point is wrong then
// it will find a sub-optimal result.
UnivariateOptimizer o = new BrentOptimizer(1e-3, 1e-6);
Plateauness f = new Plateauness(frcnum, exp_decay, low, high);
UnivariatePointValuePair p = null;
p = findMin(p, o, f, qValue, 0.1);
p = findMin(p, o, f, qValue, 0.2);
p = findMin(p, o, f, qValue, 0.333);
p = findMin(p, o, f, qValue, 0.5);
// Do some Simplex repeats as well
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
p = findMin(p, opt, f, qValue * 0.1);
p = findMin(p, opt, f, qValue * 0.5);
p = findMin(p, opt, f, qValue);
p = findMin(p, opt, f, qValue * 2);
p = findMin(p, opt, f, qValue * 10);
if (p != null)
qValue = p.getPoint();
}
QPlot qplot = new QPlot(frcCurve, qValue, low, high);
// Interactive dialog to estimate Q (blinking events per flourophore) using
// sliders for the mean and standard deviation of the localisation precision.
showQEstimationDialog(histogram, qplot, frcCurve, images.nmPerPixel);
IJ.showStatus(TITLE + " complete");
}
use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.
the class JumpDistanceAnalysis method createCMAESOptimizer.
private CMAESOptimizer createCMAESOptimizer() {
double rel = 1e-8;
double abs = 1e-10;
int maxIterations = 2000;
//Double.NEGATIVE_INFINITY;
double stopFitness = 0;
boolean isActiveCMA = true;
int diagonalOnly = 20;
int checkFeasableCount = 1;
RandomGenerator random = new Well19937c();
boolean generateStatistics = false;
ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(rel, abs);
// Iterate this for stability in the initial guess
return new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
}
use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.
the class PCPALMMolecules method optimiseSimplex.
private double[] optimiseSimplex(float[] x, float[] y, double[] initialSolution) {
// Simplex optimisation
SkewNormalMultivariateFunction sn2 = new SkewNormalMultivariateFunction(initialSolution);
sn2.addData(x, y);
NelderMeadSimplex simplex = new NelderMeadSimplex(4);
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
PointValuePair solution = opt.optimize(new MaxEval(1000), new InitialGuess(initialSolution), simplex, new ObjectiveFunction(sn2), GoalType.MINIMIZE);
double[] skewParameters2 = solution.getPointRef();
return skewParameters2;
}
use of org.apache.commons.math3.optim.PointValuePair in project GDSC-SMLM by aherbert.
the class BFGSOptimizer method doOptimize.
/** {@inheritDoc} */
@Override
protected PointValuePair doOptimize() {
final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
double[] p = getStartPoint();
// Assume minimisation
sign = -1;
LineStepSearch lineSearch = new LineStepSearch();
// In case there are no restarts
if (restarts <= 0)
return bfgsWithRoundoffCheck(checker, p, lineSearch);
PointValuePair lastResult = null;
PointValuePair result = null;
//int lastConverge = 0;
int iteration = 0;
//int[] count = new int[3];
while (iteration <= restarts) {
iteration++;
result = bfgsWithRoundoffCheck(checker, p, lineSearch);
if (converged == GRADIENT) {
// If no gradient remains then we cannot move anywhere so return
break;
}
if (lastResult != null) {
// Check if the optimum was improved using the convergence criteria
if (checker != null && checker.converged(getIterations(), lastResult, result)) {
break;
}
if (positionChecker.converged(lastResult.getPointRef(), result.getPointRef())) {
break;
}
}
// Store the new optimum and repeat
lastResult = result;
//lastConverge = converged;
p = lastResult.getPointRef();
}
return result;
}
Aggregations