use of org.apache.commons.math3.fraction.Fraction in project GDSC-SMLM by aherbert.
the class PSFCreator method run.
/*
* (non-Javadoc)
*
* @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor)
*/
public void run(ImageProcessor ip) {
loadConfiguration();
BasePoint[] spots = getSpots();
if (spots.length == 0) {
IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
return;
}
ImageStack stack = getImageStack();
final int width = imp.getWidth();
final int height = imp.getHeight();
final int currentSlice = imp.getSlice();
// Adjust settings for a single maxima
config.setIncludeNeighbours(false);
fitConfig.setDuplicateDistance(0);
ArrayList<double[]> centres = new ArrayList<double[]>(spots.length);
int iterations = 1;
LoessInterpolator loess = new LoessInterpolator(smoothing, iterations);
// TODO - The fitting routine may not produce many points. In this instance the LOESS interpolator
// fails to smooth the data very well. A higher bandwidth helps this but perhaps
// try a different smoothing method.
// For each spot
Utils.log(TITLE + ": " + imp.getTitle());
Utils.log("Finding spot locations...");
Utils.log(" %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
StoredDataStatistics averageSd = new StoredDataStatistics();
StoredDataStatistics averageA = new StoredDataStatistics();
Statistics averageRange = new Statistics();
MemoryPeakResults allResults = new MemoryPeakResults();
allResults.setName(TITLE);
allResults.setBounds(new Rectangle(0, 0, width, height));
MemoryPeakResults.addResults(allResults);
for (int n = 1; n <= spots.length; n++) {
BasePoint spot = spots[n - 1];
final int x = (int) spot.getX();
final int y = (int) spot.getY();
MemoryPeakResults results = fitSpot(stack, width, height, x, y);
allResults.addAllf(results.getResults());
if (results.size() < 5) {
Utils.log(" Spot %d: Not enough fit results %d", n, results.size());
continue;
}
// Get the results for the spot centre and width
double[] z = new double[results.size()];
double[] xCoord = new double[z.length];
double[] yCoord = new double[z.length];
double[] sd = new double[z.length];
double[] a = new double[z.length];
int i = 0;
for (PeakResult peak : results.getResults()) {
z[i] = peak.getFrame();
xCoord[i] = peak.getXPosition() - x;
yCoord[i] = peak.getYPosition() - y;
sd[i] = FastMath.max(peak.getXSD(), peak.getYSD());
a[i] = peak.getAmplitude();
i++;
}
// Smooth the amplitude plot
double[] smoothA = loess.smooth(z, a);
// Find the maximum amplitude
int maximumIndex = findMaximumIndex(smoothA);
// Find the range at a fraction of the max. This is smoothed to find the X/Y centre
int start = 0, stop = smoothA.length - 1;
double limit = smoothA[maximumIndex] * amplitudeFraction;
for (int j = 0; j < smoothA.length; j++) {
if (smoothA[j] > limit) {
start = j;
break;
}
}
for (int j = smoothA.length; j-- > 0; ) {
if (smoothA[j] > limit) {
stop = j;
break;
}
}
averageRange.add(stop - start + 1);
// Extract xy centre coords and smooth
double[] smoothX = new double[stop - start + 1];
double[] smoothY = new double[smoothX.length];
double[] smoothSd = new double[smoothX.length];
double[] newZ = new double[smoothX.length];
for (int j = start, k = 0; j <= stop; j++, k++) {
smoothX[k] = xCoord[j];
smoothY[k] = yCoord[j];
smoothSd[k] = sd[j];
newZ[k] = z[j];
}
smoothX = loess.smooth(newZ, smoothX);
smoothY = loess.smooth(newZ, smoothY);
smoothSd = loess.smooth(newZ, smoothSd);
// Since the amplitude is not very consistent move from this peak to the
// lowest width which is the in-focus spot.
maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
// Find the centre at the amplitude peak
double cx = smoothX[maximumIndex] + x;
double cy = smoothY[maximumIndex] + y;
int cz = (int) newZ[maximumIndex];
double csd = smoothSd[maximumIndex];
double ca = smoothA[maximumIndex + start];
// The average should weight the SD using the signal for each spot
averageSd.add(smoothSd[maximumIndex]);
averageA.add(ca);
if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) {
Utils.log(" Spot %d was ignored", n);
continue;
}
// Store result - it may have been moved interactively
maximumIndex += this.slice - cz;
cz = (int) newZ[maximumIndex];
csd = smoothSd[maximumIndex];
ca = smoothA[maximumIndex + start];
Utils.log(" Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f\n", n, cx, cy, cz, csd, ca);
centres.add(new double[] { cx, cy, cz, csd, n });
}
if (interactiveMode) {
imp.setSlice(currentSlice);
imp.setOverlay(null);
// Hide the amplitude and spot plots
Utils.hide(TITLE_AMPLITUDE);
Utils.hide(TITLE_PSF_PARAMETERS);
}
if (centres.isEmpty()) {
String msg = "No suitable spots could be identified centres";
Utils.log(msg);
IJ.error(TITLE, msg);
return;
}
// Find the limits of the z-centre
int minz = (int) centres.get(0)[2];
int maxz = minz;
for (double[] centre : centres) {
if (minz > centre[2])
minz = (int) centre[2];
else if (maxz < centre[2])
maxz = (int) centre[2];
}
IJ.showStatus("Creating PSF image");
// Create a stack that can hold all the data.
ImageStack psf = createStack(stack, minz, maxz, magnification);
// For each spot
Statistics stats = new Statistics();
boolean ok = true;
for (int i = 0; ok && i < centres.size(); i++) {
double progress = (double) i / centres.size();
final double increment = 1.0 / (stack.getSize() * centres.size());
IJ.showProgress(progress);
double[] centre = centres.get(i);
// Extract the spot
float[][] spot = new float[stack.getSize()][];
Rectangle regionBounds = null;
for (int slice = 1; slice <= stack.getSize(); slice++) {
ImageExtractor ie = new ImageExtractor((float[]) stack.getPixels(slice), width, height);
if (regionBounds == null)
regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius);
spot[slice - 1] = ie.crop(regionBounds);
}
int n = (int) centre[4];
final float b = getBackground(n, spot);
if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
Utils.log(" Spot %d was ignored", n);
continue;
}
stats.add(b);
// Adjust the centre using the crop
centre[0] -= regionBounds.x;
centre[1] -= regionBounds.y;
// This takes a long time so this should track progress
ok = addToPSF(maxz, magnification, psf, centre, spot, regionBounds, progress, increment, centreEachSlice);
}
if (interactiveMode) {
Utils.hide(TITLE_INTENSITY);
}
IJ.showProgress(1);
if (threadPool != null) {
threadPool.shutdownNow();
threadPool = null;
}
if (!ok || stats.getN() == 0)
return;
final double avSd = getAverage(averageSd, averageA, 2);
Utils.log(" Average background = %.2f, Av. SD = %s px", stats.getMean(), Utils.rounded(avSd, 4));
normalise(psf, maxz, avSd * magnification, false);
IJ.showProgress(1);
psfImp = Utils.display("PSF", psf);
psfImp.setSlice(maxz);
psfImp.resetDisplayRange();
psfImp.updateAndDraw();
double[][] fitCom = new double[2][psf.getSize()];
Arrays.fill(fitCom[0], Double.NaN);
Arrays.fill(fitCom[1], Double.NaN);
double fittedSd = fitPSF(psf, loess, maxz, averageRange.getMean(), fitCom);
// Compute the drift in the PSF:
// - Use fitted centre if available; otherwise find CoM for each frame
// - express relative to the average centre
double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / magnification);
double[] slice = Utils.newArray(psf.getSize(), 1, 1.0);
String title = TITLE + " CoM Drift";
Plot2 plot = new Plot2(title, "Slice", "Drift (nm)");
plot.addLabel(0, 0, "Red = X; Blue = Y");
//double[] limitsX = Maths.limits(com[0]);
//double[] limitsY = Maths.limits(com[1]);
double[] limitsX = getLimits(com[0]);
double[] limitsY = getLimits(com[1]);
plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1]));
plot.setColor(Color.red);
plot.addPoints(slice, com[0], Plot.DOT);
plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(slice, com[1], Plot.DOT);
plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE);
Utils.display(title, plot);
// TODO - Redraw the PSF with drift correction applied.
// This means that the final image should have no drift.
// This is relevant when combining PSF images. It doesn't matter too much for simulations
// unless the drift is large.
// Add Image properties containing the PSF details
final double fwhm = getFWHM(psf, maxz);
psfImp.setProperty("Info", XmlUtils.toXML(new PSFSettings(maxz, nmPerPixel / magnification, nmPerSlice, stats.getN(), fwhm, createNote())));
Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, Utils.rounded(nmPerPixel / magnification, 3), Utils.rounded(nmPerSlice, 3), stats.getN(), Utils.rounded(fittedSd * nmPerPixel, 4), Utils.rounded(fwhm));
createInteractivePlots(psf, maxz, nmPerPixel / magnification, fittedSd * nmPerPixel);
IJ.showStatus("");
}
use of org.apache.commons.math3.fraction.Fraction 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.fraction.Fraction 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.fraction.Fraction in project GDSC-SMLM by aherbert.
the class BenchmarkFilterAnalysis method filterAnalysis.
private int filterAnalysis(FilterSet filterSet, int setNumber, DirectFilter currentOptimum, double rangeReduction) {
// Check if the filters are the same so allowing optimisation
final boolean allSameType = filterSet.allSameType();
this.ga_resultsList = resultsList;
Chromosome<FilterScore> best = null;
String algorithm = "";
// All the search algorithms use search dimensions.
// Create search dimensions if needed (these are used for testing if the optimum is at the limit).
ss_filter = null;
ss_lower = null;
ss_upper = null;
FixedDimension[] originalDimensions = null;
boolean rangeInput = false;
boolean[] disabled = null;
double[][] seed = null;
boolean nonInteractive = false;
if (allSameType) {
// There should always be 1 filter
ss_filter = (DirectFilter) filterSet.getFilters().get(0);
int n = ss_filter.getNumberOfParameters();
// Option to configure a range
rangeInput = filterSet.getName().contains("Range");
double[] range = new double[n];
if (rangeInput && filterSet.size() == 4) {
originalDimensions = new FixedDimension[n];
// This is used as min/lower/upper/max
final Filter minF = ss_filter;
final Filter lowerF = filterSet.getFilters().get(1);
final Filter upperF = filterSet.getFilters().get(2);
final Filter maxF = filterSet.getFilters().get(3);
for (int i = 0; i < n; i++) {
double min = minF.getParameterValue(i);
double lower = lowerF.getParameterValue(i);
double upper = upperF.getParameterValue(i);
range[i] = upper - lower;
double max = maxF.getParameterValue(i);
double minIncrement = ss_filter.getParameterIncrement(i);
try {
originalDimensions[i] = new FixedDimension(min, max, minIncrement, lower, upper);
} catch (IllegalArgumentException e) {
Utils.log(TITLE + " : Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i));
originalDimensions = null;
rangeInput = false;
break;
}
}
}
if (rangeInput && (filterSet.size() == 3 || filterSet.size() == 2)) {
originalDimensions = new FixedDimension[n];
// This is used as lower/upper[/increment]
final Filter lowerF = ss_filter;
final Filter upperF = filterSet.getFilters().get(1);
for (int i = 0; i < n; i++) {
// Do not disable if the increment is not set. This is left to the user to decide.
// if (incF.getParameterValue(i) == incF.getDisabledParameterValue(i) ||
// Double.isInfinite(incF.getParameterValue(i)))
// {
// // Not enabled
// dimensions[i] = new SearchDimension(incF.getDisabledParameterValue(i));
// continue;
// }
double lower = lowerF.getParameterValue(i);
double upper = upperF.getParameterValue(i);
range[i] = upper - lower;
ParameterType type = ss_filter.getParameterType(i);
double min = BenchmarkSpotFit.getMin(type);
double max = BenchmarkSpotFit.getMax(type);
double minIncrement = ss_filter.getParameterIncrement(i);
try {
originalDimensions[i] = new FixedDimension(min, max, minIncrement, lower, upper);
} catch (IllegalArgumentException e) {
Utils.log(TITLE + " : Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i));
originalDimensions = null;
rangeInput = false;
break;
}
}
}
// Get the dimensions from the filters
if (originalDimensions == null) {
originalDimensions = new FixedDimension[n];
// Allow inputing a filter set (e.g. saved from previous optimisation)
// Find the limits in the current scores
final double[] lower = ss_filter.getParameters().clone();
final double[] upper = lower.clone();
// Allow the SearchSpace algorithms to be seeded with an initial population
// for the first evaluation of the optimum. This is done when the input filter
// set is not a range.
seed = new double[filterSet.size()][];
int c = 0;
for (Filter f : filterSet.getFilters()) {
final double[] point = f.getParameters();
seed[c++] = point;
for (int j = 0; j < lower.length; j++) {
if (lower[j] > point[j])
lower[j] = point[j];
if (upper[j] < point[j])
upper[j] = point[j];
}
}
// Min/max must be set using values from BenchmarkSpotFit.
for (int i = 0; i < n; i++) {
if (lower[i] == upper[i]) {
// Not enabled
originalDimensions[i] = new FixedDimension(lower[i]);
continue;
}
ParameterType type = ss_filter.getParameterType(i);
double min = BenchmarkSpotFit.getMin(type);
double max = BenchmarkSpotFit.getMax(type);
double minIncrement = ss_filter.getParameterIncrement(i);
if (min > lower[i])
min = lower[i];
if (max < upper[i])
max = upper[i];
try {
originalDimensions[i] = new FixedDimension(min, max, minIncrement, lower[i], upper[i]);
} catch (IllegalArgumentException e) {
Utils.log(TITLE + " : Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i));
originalDimensions = null;
break;
}
}
if (originalDimensions == null) {
// Failed to work out the dimensions. No optimisation will be possible.
// Sort so that the filters are in a nice order for reporting
filterSet.sort();
// This will not be used when the dimensions are null
seed = null;
}
}
if (originalDimensions != null) {
// Use the current optimum if we are doing a range optimisation
if (currentOptimum != null && rangeInput && currentOptimum.getType().equals(ss_filter.getType()) && evolve != 0) {
// Suppress dialogs and use the current settings
nonInteractive = true;
double[] p = currentOptimum.getParameters();
// Range search uses SearchDimension and we must centre on the optimum after creation.
for (int i = 0; i < originalDimensions.length; i++) {
double centre = p[i];
double r = 0;
if (originalDimensions[i].isActive()) {
// Set the range around the centre.
// This uses the range for each param when we read the filters.
r = range[i];
// Optionally reduce the width of the dimensions.
if (rangeReduction > 0 && rangeReduction < 1)
r *= rangeReduction;
}
double lower = centre - r * 0.5;
double upper = centre + r * 0.5;
originalDimensions[i] = originalDimensions[i].create(lower, upper);
}
}
// Store the dimensions so we can do an 'at limit' check
disabled = new boolean[originalDimensions.length];
ss_lower = new double[originalDimensions.length];
ss_upper = new double[originalDimensions.length];
for (int i = 0; i < disabled.length; i++) {
disabled[i] = !originalDimensions[i].isActive();
ss_lower[i] = originalDimensions[i].lower;
ss_upper[i] = originalDimensions[i].upper;
}
}
} else {
// Sort so that the filters are in a nice order for reporting
filterSet.sort();
}
analysisStopWatch = StopWatch.createStarted();
if (evolve == 1 && originalDimensions != null) {
// Collect parameters for the genetic algorithm
pauseFilterTimer();
// Remember the step size settings
double[] stepSize = stepSizeMap.get(setNumber);
if (stepSize == null || stepSize.length != ss_filter.length()) {
stepSize = ss_filter.mutationStepRange().clone();
for (int j = 0; j < stepSize.length; j++) stepSize[j] *= delta;
// See if the same number of parameters have been optimised in other algorithms
boolean[] enabled = searchRangeMap.get(setNumber);
if (enabled != null && enabled.length == stepSize.length) {
for (int j = 0; j < stepSize.length; j++) if (!enabled[j])
stepSize[j] *= -1;
}
}
GenericDialog gd = null;
int[] indices = ss_filter.getChromosomeParameters();
boolean runAlgorithm = nonInteractive;
if (!nonInteractive) {
// Ask the user for the mutation step parameters.
gd = new GenericDialog(TITLE);
String prefix = setNumber + "_";
gd.addMessage("Configure the genetic algorithm for [" + setNumber + "] " + filterSet.getName());
gd.addNumericField(prefix + "Population_size", populationSize, 0);
gd.addNumericField(prefix + "Failure_limit", failureLimit, 0);
gd.addNumericField(prefix + "Tolerance", tolerance, -1);
gd.addNumericField(prefix + "Converged_count", convergedCount, 0);
gd.addSlider(prefix + "Mutation_rate", 0.05, 1, mutationRate);
gd.addSlider(prefix + "Crossover_rate", 0.05, 1, crossoverRate);
gd.addSlider(prefix + "Mean_children", 0.05, 3, meanChildren);
gd.addSlider(prefix + "Selection_fraction", 0.05, 0.5, selectionFraction);
gd.addCheckbox(prefix + "Ramped_selection", rampedSelection);
gd.addCheckbox(prefix + "Save_option", saveOption);
gd.addMessage("Configure the step size for each parameter");
for (int j = 0; j < indices.length; j++) {
// Do not mutate parameters that were not expanded, i.e. the input did not vary them.
final double step = (originalDimensions[indices[j]].isActive()) ? stepSize[j] : 0;
gd.addNumericField(getDialogName(prefix, ss_filter, indices[j]), step, 2);
}
gd.showDialog();
runAlgorithm = !gd.wasCanceled();
}
if (runAlgorithm) {
// Used to create random sample
FixedDimension[] dimensions = Arrays.copyOf(originalDimensions, originalDimensions.length);
if (!nonInteractive) {
populationSize = (int) Math.abs(gd.getNextNumber());
if (populationSize < 10)
populationSize = 10;
failureLimit = (int) Math.abs(gd.getNextNumber());
tolerance = gd.getNextNumber();
// Allow negatives
convergedCount = (int) gd.getNextNumber();
mutationRate = Math.abs(gd.getNextNumber());
crossoverRate = Math.abs(gd.getNextNumber());
meanChildren = Math.abs(gd.getNextNumber());
selectionFraction = Math.abs(gd.getNextNumber());
rampedSelection = gd.getNextBoolean();
saveOption = gd.getNextBoolean();
for (int j = 0; j < indices.length; j++) {
stepSize[j] = gd.getNextNumber();
}
// Store for repeat analysis
stepSizeMap.put(setNumber, stepSize);
}
for (int j = 0; j < indices.length; j++) {
// A zero step size will keep the parameter but prevent range mutation.
if (stepSize[j] < 0) {
dimensions[indices[j]] = new FixedDimension(ss_filter.getDisabledParameterValue(indices[j]));
disabled[indices[j]] = true;
}
}
// // Reset negatives to zero
// stepSize = stepSize.clone();
// for (int j = 0; j < stepSize.length; j++)
// if (stepSize[j] < 0)
// stepSize[j] = 0;
// Create the genetic algorithm
RandomDataGenerator random = new RandomDataGenerator(new Well44497b());
SimpleMutator<FilterScore> mutator = new SimpleMutator<FilterScore>(random, mutationRate);
// Override the settings with the step length, a min of zero and the configured upper
double[] upper = ss_filter.upperLimit();
mutator.overrideChromosomeSettings(stepSize, new double[stepSize.length], upper);
Recombiner<FilterScore> recombiner = new SimpleRecombiner<FilterScore>(random, crossoverRate, meanChildren);
SelectionStrategy<FilterScore> selectionStrategy;
// If the initial population is huge ensure that the first selection culls to the correct size
final int selectionMax = (int) (selectionFraction * populationSize);
if (rampedSelection)
selectionStrategy = new RampedSelectionStrategy<FilterScore>(random, selectionFraction, selectionMax);
else
selectionStrategy = new SimpleSelectionStrategy<FilterScore>(random, selectionFraction, selectionMax);
ToleranceChecker<FilterScore> ga_checker = new InterruptChecker(tolerance, tolerance * 1e-3, convergedCount);
// Create new random filters if the population is initially below the population size
List<Filter> filters = filterSet.getFilters();
if (filterSet.size() < populationSize) {
filters = new ArrayList<Filter>(populationSize);
// Add the existing filters if they are not a range input file
if (!rangeInput)
filters.addAll(filterSet.getFilters());
// Add current optimum to seed
if (nonInteractive)
filters.add(currentOptimum);
// The GA does not use the min interval grid so sample without rounding
double[][] sample = SearchSpace.sampleWithoutRounding(dimensions, populationSize - filters.size(), null);
filters.addAll(searchSpaceToFilters(sample));
}
ga_population = new Population<FilterScore>(filters);
ga_population.setPopulationSize(populationSize);
ga_population.setFailureLimit(failureLimit);
selectionStrategy.setTracker(this);
// Evolve
algorithm = EVOLVE[evolve];
ga_statusPrefix = algorithm + " [" + setNumber + "] " + filterSet.getName() + " ... ";
ga_iteration = 0;
ga_population.setTracker(this);
createGAWindow();
resumeFilterTimer();
best = ga_population.evolve(mutator, recombiner, this, selectionStrategy, ga_checker);
if (best != null) {
// In case optimisation was stopped
IJ.resetEscape();
// The GA may produce coordinates off the min interval grid
best = enumerateMinInterval(best, stepSize, indices);
// Now update the filter set for final assessment
filterSet = new FilterSet(filterSet.getName(), populationToFilters(ga_population.getIndividuals()));
// Option to save the filters
if (saveOption)
saveFilterSet(filterSet, setNumber, !nonInteractive);
}
} else
resumeFilterTimer();
}
if ((evolve == 2 || evolve == 4) && originalDimensions != null) {
// Collect parameters for the range search algorithm
pauseFilterTimer();
boolean isStepSearch = evolve == 4;
// The step search should use a multi-dimension refinement and no range reduction
SearchSpace.RefinementMode myRefinementMode = SearchSpace.RefinementMode.MULTI_DIMENSION;
// Remember the enabled settings
boolean[] enabled = searchRangeMap.get(setNumber);
int n = ss_filter.getNumberOfParameters();
if (enabled == null || enabled.length != n) {
enabled = new boolean[n];
Arrays.fill(enabled, true);
// See if the same number of parameters have been optimised in other algorithms
double[] stepSize = stepSizeMap.get(setNumber);
if (stepSize != null && enabled.length == stepSize.length) {
for (int j = 0; j < stepSize.length; j++) if (stepSize[j] < 0)
enabled[j] = false;
}
}
GenericDialog gd = null;
boolean runAlgorithm = nonInteractive;
if (!nonInteractive) {
// Ask the user for the search parameters.
gd = new GenericDialog(TITLE);
String prefix = setNumber + "_";
gd.addMessage("Configure the " + EVOLVE[evolve] + " algorithm for [" + setNumber + "] " + filterSet.getName());
gd.addSlider(prefix + "Width", 1, 5, rangeSearchWidth);
if (!isStepSearch) {
gd.addCheckbox(prefix + "Save_option", saveOption);
gd.addNumericField(prefix + "Max_iterations", maxIterations, 0);
String[] modes = SettingsManager.getNames((Object[]) SearchSpace.RefinementMode.values());
gd.addSlider(prefix + "Reduce", 0.01, 0.99, rangeSearchReduce);
gd.addChoice("Refinement", modes, modes[refinementMode]);
}
gd.addNumericField(prefix + "Seed_size", seedSize, 0);
// Add choice of fields to optimise
for (int i = 0; i < n; i++) gd.addCheckbox(getDialogName(prefix, ss_filter, i), enabled[i]);
gd.showDialog();
runAlgorithm = !gd.wasCanceled();
}
if (runAlgorithm) {
SearchDimension[] dimensions = new SearchDimension[n];
if (!nonInteractive) {
rangeSearchWidth = (int) gd.getNextNumber();
if (!isStepSearch) {
saveOption = gd.getNextBoolean();
maxIterations = (int) gd.getNextNumber();
refinementMode = gd.getNextChoiceIndex();
rangeSearchReduce = gd.getNextNumber();
}
seedSize = (int) gd.getNextNumber();
for (int i = 0; i < n; i++) enabled[i] = gd.getNextBoolean();
searchRangeMap.put(setNumber, enabled);
}
if (!isStepSearch)
myRefinementMode = SearchSpace.RefinementMode.values()[refinementMode];
for (int i = 0; i < n; i++) {
if (enabled[i]) {
try {
dimensions[i] = originalDimensions[i].create(rangeSearchWidth);
dimensions[i].setPad(true);
// Prevent range reduction so that the step search just does a single refinement step
dimensions[i].setReduceFactor((isStepSearch) ? 1 : rangeSearchReduce);
// Centre on current optimum
if (nonInteractive)
dimensions[i].setCentre(currentOptimum.getParameterValue(i));
} catch (IllegalArgumentException e) {
IJ.error(TITLE, String.format("Unable to configure dimension [%d] %s: " + e.getMessage(), i, ss_filter.getParameterName(i)));
return -1;
}
} else {
dimensions[i] = new SearchDimension(ss_filter.getDisabledParameterValue(i));
}
}
for (int i = 0; i < disabled.length; i++) disabled[i] = !dimensions[i].active;
// Check the number of combinations is OK
long combinations = SearchSpace.countCombinations(dimensions);
if (!nonInteractive && combinations > 10000) {
gd = new GenericDialog(TITLE);
gd.addMessage(String.format("%d combinations for the configured dimensions.\n \nClick 'Yes' to optimise.", combinations));
gd.enableYesNoCancel();
gd.hideCancelButton();
gd.showDialog();
if (!gd.wasOKed()) {
combinations = 0;
}
}
if (combinations == 0) {
resumeFilterTimer();
} else {
algorithm = EVOLVE[evolve] + " " + rangeSearchWidth;
ga_statusPrefix = algorithm + " [" + setNumber + "] " + filterSet.getName() + " ... ";
ga_iteration = 0;
es_optimum = null;
SearchSpace ss = new SearchSpace();
ss.setTracker(this);
if (seedSize > 0) {
double[][] sample;
// Add current optimum to seed
if (nonInteractive) {
sample = new double[1][];
sample[0] = currentOptimum.getParameters();
seed = merge(seed, sample);
}
int size = (seed == null) ? 0 : seed.length;
// Sample without rounding as the seed will be rounded
sample = SearchSpace.sampleWithoutRounding(dimensions, seedSize - size, null);
seed = merge(seed, sample);
}
// Note: If we have an optimum and we are not seeding this should not matter as the dimensions
// have been centred on the current optimum
ss.seed(seed);
ConvergenceChecker<FilterScore> checker = new InterruptConvergenceChecker(0, 0, maxIterations);
createGAWindow();
resumeFilterTimer();
SearchResult<FilterScore> optimum = ss.search(dimensions, this, checker, myRefinementMode);
if (optimum != null) {
// In case optimisation was stopped
IJ.resetEscape();
best = ((SimpleFilterScore) optimum.score).r.filter;
if (seedSize > 0) {
// Not required as the search now respects the min interval
// The optimum may be off grid if it was from the seed
//best = enumerateMinInterval(best, enabled);
}
// Now update the filter set for final assessment
filterSet = new FilterSet(filterSet.getName(), searchSpaceToFilters((DirectFilter) best, ss.getSearchSpace()));
// Option to save the filters
if (saveOption)
saveFilterSet(filterSet, setNumber, !nonInteractive);
}
}
} else
resumeFilterTimer();
}
if (evolve == 3 && originalDimensions != null) {
// Collect parameters for the enrichment search algorithm
pauseFilterTimer();
boolean[] enabled = searchRangeMap.get(setNumber);
int n = ss_filter.getNumberOfParameters();
if (enabled == null || enabled.length != n) {
enabled = new boolean[n];
Arrays.fill(enabled, true);
// See if the same number of parameters have been optimised in other algorithms
double[] stepSize = stepSizeMap.get(setNumber);
if (stepSize != null && enabled.length == stepSize.length) {
for (int j = 0; j < stepSize.length; j++) if (stepSize[j] < 0)
enabled[j] = false;
}
}
GenericDialog gd = null;
boolean runAlgorithm = nonInteractive;
if (!nonInteractive) {
// Ask the user for the search parameters.
gd = new GenericDialog(TITLE);
String prefix = setNumber + "_";
gd.addMessage("Configure the enrichment search algorithm for [" + setNumber + "] " + filterSet.getName());
gd.addCheckbox(prefix + "Save_option", saveOption);
gd.addNumericField(prefix + "Max_iterations", maxIterations, 0);
gd.addNumericField(prefix + "Converged_count", convergedCount, 0);
gd.addNumericField(prefix + "Samples", enrichmentSamples, 0);
gd.addSlider(prefix + "Fraction", 0.01, 0.99, enrichmentFraction);
gd.addSlider(prefix + "Padding", 0, 0.99, enrichmentPadding);
// Add choice of fields to optimise
for (int i = 0; i < n; i++) gd.addCheckbox(getDialogName(prefix, ss_filter, i), enabled[i]);
gd.showDialog();
runAlgorithm = !gd.wasCanceled();
}
if (runAlgorithm) {
FixedDimension[] dimensions = Arrays.copyOf(originalDimensions, originalDimensions.length);
if (!nonInteractive) {
saveOption = gd.getNextBoolean();
maxIterations = (int) gd.getNextNumber();
convergedCount = (int) gd.getNextNumber();
enrichmentSamples = (int) gd.getNextNumber();
enrichmentFraction = gd.getNextNumber();
enrichmentPadding = gd.getNextNumber();
for (int i = 0; i < n; i++) enabled[i] = gd.getNextBoolean();
searchRangeMap.put(setNumber, enabled);
}
for (int i = 0; i < n; i++) {
if (!enabled[i])
dimensions[i] = new FixedDimension(ss_filter.getDisabledParameterValue(i));
}
for (int i = 0; i < disabled.length; i++) disabled[i] = !dimensions[i].active;
algorithm = EVOLVE[evolve];
ga_statusPrefix = algorithm + " [" + setNumber + "] " + filterSet.getName() + " ... ";
ga_iteration = 0;
es_optimum = null;
SearchSpace ss = new SearchSpace();
ss.setTracker(this);
// Add current optimum to seed
if (nonInteractive) {
double[][] sample = new double[1][];
sample[0] = currentOptimum.getParameters();
seed = merge(seed, sample);
}
ss.seed(seed);
ConvergenceChecker<FilterScore> checker = new InterruptConvergenceChecker(0, 0, maxIterations, convergedCount);
createGAWindow();
resumeFilterTimer();
SearchResult<FilterScore> optimum = ss.enrichmentSearch(dimensions, this, checker, enrichmentSamples, enrichmentFraction, enrichmentPadding);
if (optimum != null) {
// In case optimisation was stopped
IJ.resetEscape();
best = ((SimpleFilterScore) optimum.score).r.filter;
// Not required as the search now respects the min interval
// Enumerate on the min interval to produce the final filter
//best = enumerateMinInterval(best, enabled);
// Now update the filter set for final assessment
filterSet = new FilterSet(filterSet.getName(), searchSpaceToFilters((DirectFilter) best, ss.getSearchSpace()));
// Option to save the filters
if (saveOption)
saveFilterSet(filterSet, setNumber, !nonInteractive);
}
} else
resumeFilterTimer();
}
IJ.showStatus("Analysing [" + setNumber + "] " + filterSet.getName() + " ...");
// Do not support plotting if we used optimisation
double[] xValues = (best != null || isHeadless || (plotTopN == 0)) ? null : new double[filterSet.size()];
double[] yValues = (xValues == null) ? null : new double[xValues.length];
SimpleFilterScore max = null;
// It can just assess the top 1 required for the summary.
if (best != null) {
// Only assess the top 1 filter for the summary
List<Filter> list = new ArrayList<Filter>();
list.add((DirectFilter) best);
filterSet = new FilterSet(filterSet.getName(), list);
}
// Score the filters and report the results if configured.
FilterScoreResult[] scoreResults = scoreFilters(setUncomputedStrength(filterSet), showResultsTable);
if (scoreResults == null)
return -1;
analysisStopWatch.stop();
for (int index = 0; index < scoreResults.length; index++) {
final FilterScoreResult scoreResult = scoreResults[index];
if (xValues != null) {
xValues[index] = scoreResult.filter.getNumericalValue();
yValues[index] = scoreResult.score;
}
final SimpleFilterScore result = new SimpleFilterScore(scoreResult, allSameType, scoreResult.criteria >= minCriteria);
if (result.compareTo(max) < 0) {
max = result;
}
}
if (showResultsTable) {
BufferedTextWindow tw = null;
if (resultsWindow != null) {
tw = new BufferedTextWindow(resultsWindow);
tw.setIncrement(Integer.MAX_VALUE);
}
for (int index = 0; index < scoreResults.length; index++) addToResultsWindow(tw, scoreResults[index].text);
if (resultsWindow != null)
resultsWindow.getTextPanel().updateDisplay();
}
// Check the top filter against the limits of the original dimensions
char[] atLimit = null;
if (allSameType && originalDimensions != null) {
DirectFilter filter = max.r.filter;
int[] indices = filter.getChromosomeParameters();
atLimit = new char[indices.length];
StringBuilder sb = new StringBuilder(200);
for (int j = 0; j < indices.length; j++) {
atLimit[j] = ComplexFilterScore.WITHIN;
final int p = indices[j];
if (disabled[p])
continue;
final double value = filter.getParameterValue(p);
double lowerLimit = originalDimensions[p].getLower();
double upperLimit = originalDimensions[p].getUpper();
int c1 = Double.compare(value, lowerLimit);
if (c1 <= 0) {
atLimit[j] = ComplexFilterScore.FLOOR;
sb.append(" : ").append(filter.getParameterName(p)).append(' ').append(atLimit[j]).append('[').append(Utils.rounded(value));
if (c1 == -1) {
atLimit[j] = ComplexFilterScore.BELOW;
sb.append("<").append(Utils.rounded(lowerLimit));
}
sb.append("]");
} else {
int c2 = Double.compare(value, upperLimit);
if (c2 >= 0) {
atLimit[j] = ComplexFilterScore.CEIL;
sb.append(" : ").append(filter.getParameterName(p)).append(' ').append(atLimit[j]).append('[').append(Utils.rounded(value));
if (c2 == 1) {
atLimit[j] = ComplexFilterScore.ABOVE;
sb.append(">").append(Utils.rounded(upperLimit));
}
sb.append("]");
}
}
}
if (sb.length() > 0) {
if (max.criteriaPassed) {
Utils.log("Warning: Top filter (%s @ %s|%s) [%s] at the limit of the expanded range%s", filter.getName(), Utils.rounded((invertScore) ? -max.score : max.score), Utils.rounded((invertCriteria) ? -minCriteria : minCriteria), limitFailCount + limitRange, sb.toString());
} else {
Utils.log("Warning: Top filter (%s @ -|%s) [%s] at the limit of the expanded range%s", filter.getName(), Utils.rounded((invertCriteria) ? -max.criteria : max.criteria), limitFailCount + limitRange, sb.toString());
}
}
}
// Note that max should never be null since this method is not run with an empty filter set
// We may have no filters that pass the criteria
String type = max.r.filter.getType();
if (!max.criteriaPassed) {
Utils.log("Warning: Filter does not pass the criteria: %s : Best = %s using %s", type, Utils.rounded((invertCriteria) ? -max.criteria : max.criteria), max.r.filter.getName());
return 0;
}
// This could be an option?
boolean allowDuplicates = true;
// XXX - Commented out the requirement to be the same type to store for later analysis.
// This may break the code, however I think that all filter sets should be able to have a best filter
// irrespective of whether they were the same type or not.
//if (allSameType)
//{
ComplexFilterScore newFilterScore = new ComplexFilterScore(max.r, atLimit, algorithm, analysisStopWatch.getTime(), "", 0);
addBestFilter(type, allowDuplicates, newFilterScore);
// Add spacer at end of each result set
if (isHeadless) {
if (showResultsTable && filterSet.size() > 1)
IJ.log("");
} else {
if (showResultsTable && filterSet.size() > 1)
resultsWindow.append("");
if (plotTopN > 0 && xValues != null) {
// Check the xValues are unique. Since the filters have been sorted by their
// numeric value we only need to compare adjacent entries.
boolean unique = true;
for (int ii = 0; ii < xValues.length - 1; ii++) {
if (xValues[ii] == xValues[ii + 1]) {
unique = false;
break;
}
}
String xAxisName = filterSet.getValueName();
if (unique) {
// Check the values all refer to the same property
for (Filter filter : filterSet.getFilters()) {
if (!xAxisName.equals(filter.getNumericalValueName())) {
unique = false;
break;
}
}
}
if (!unique) {
// If not unique then renumber them and use an arbitrary label
xAxisName = "Filter";
for (int ii = 0; ii < xValues.length; ii++) xValues[ii] = ii + 1;
}
String title = filterSet.getName();
// Check if a previous filter set had the same name, update if necessary
NamedPlot p = getNamedPlot(title);
if (p == null)
plots.add(new NamedPlot(title, xAxisName, xValues, yValues));
else
p.updateValues(xAxisName, xValues, yValues);
if (plots.size() > plotTopN) {
Collections.sort(plots);
p = plots.remove(plots.size() - 1);
}
}
}
return 0;
}
use of org.apache.commons.math3.fraction.Fraction in project GDSC-SMLM by aherbert.
the class SCMOSLikelihoodWrapperTest method functionComputesTargetGradientPerDatum.
private void functionComputesTargetGradientPerDatum(Gaussian2DFunction f1, int targetParameter) {
int[] indices = f1.gradientIndices();
int gradientIndex = findGradientIndex(f1, targetParameter);
double[] dyda = new double[indices.length];
double[] a;
SCMOSLikelihoodWrapper ff1;
int n = maxx * maxx;
int count = 0, total = 0;
RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
for (double background : testbackground) for (double signal1 : testsignal1) for (double angle1 : testangle1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
a = createParameters(background, signal1, angle1, cx1, cy1, w1[0], w1[1]);
// Create y as a function we would want to move towards
double[] a2 = a.clone();
a2[targetParameter] *= 1.1;
f1.initialise(a2);
double[] data = new double[n];
for (int i = 0; i < n; i++) {
// Simulate sCMOS camera
double u = f1.eval(i);
data[i] = rdg.nextPoisson(u) * g[i] + rdg.nextGaussian(o[i], sd[i]);
}
ff1 = new SCMOSLikelihoodWrapper(f1, a, data, n, var, g, o);
// Numerically solve gradient.
// Calculate the step size h to be an exact numerical representation
final double xx = a[targetParameter];
// Get h to minimise roundoff error
double h = Precision.representableDelta(xx, h_);
for (int x : testx) for (int y : testy) {
int i = y * maxx + x;
a[targetParameter] = xx;
ff1.likelihood(getVariables(indices, a), dyda, i);
// Evaluate at (x+h) and (x-h)
a[targetParameter] = xx + h;
double value2 = ff1.likelihood(getVariables(indices, a), i);
a[targetParameter] = xx - h;
double value3 = ff1.likelihood(getVariables(indices, a), i);
double gradient = (value2 - value3) / (2 * h);
boolean ok = Math.signum(gradient) == Math.signum(dyda[gradientIndex]) || Math.abs(gradient - dyda[gradientIndex]) < 0.1;
// dyda[gradientIndex]);
if (!ok)
Assert.assertTrue(NAME[targetParameter] + ": " + gradient + " != " + dyda[gradientIndex], ok);
ok = eqPerDatum.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]);
if (ok)
count++;
total++;
}
}
double p = (100.0 * count) / total;
logf("Per Datum %s : %s = %d / %d (%.2f)\n", f1.getClass().getSimpleName(), NAME[targetParameter], count, total, p);
Assert.assertTrue(NAME[targetParameter] + " fraction too low per datum: " + p, p > 90);
}
Aggregations