use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class ResultsMatchCalculator method getRoiCoordinates.
@Nullable
private static double[] getRoiCoordinates(String line) {
// Extract the startT and x,y coordinates from the first available point in the pair
final int[] index = { 1, 4 };
final String[] fields = line.split("\t");
final int startT = Integer.parseInt(fields[0]);
for (final int i : index) {
if (i < fields.length) {
if (fields[i].equals("-")) {
continue;
}
final double x = Double.parseDouble(fields[i]);
final double y = Double.parseDouble(fields[i + 1]);
return new double[] { startT, x, y };
}
}
return null;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable 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.
*/
@Nullable
private double[][] doFitJumpDistancesMle(double[] jumpDistances, double estimatedD, int n) {
final MaxEval maxEval = new MaxEval(20000);
final 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
final SimpleBounds bounds = new SimpleBounds(function.getLowerBounds(), function.getUpperBounds(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY));
final PointValuePair solution = powellOptimizer.optimize(maxEval, new ObjectiveFunction(function), new InitialGuess(function.guess()), bounds, new CustomPowellOptimizer.BasisStep(function.step()), GoalType.MAXIMIZE);
final double[] fitParams = solution.getPointRef();
ll = solution.getValue();
lastFitValue = fitValue = MathUtils.getAkaikeInformationCriterion(ll, 1);
final double[] coefficients = fitParams;
final double[] fractions = new double[] { 1 };
LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=1) : %s, MLE = %s, Akaike IC = %s (%d evaluations)", formatD(fitParams[0]), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), powellOptimizer.getEvaluations());
return new double[][] { coefficients, fractions };
} catch (final TooManyEvaluationsException ex) {
LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many evaluation (%d)", powellOptimizer.getEvaluations());
} catch (final TooManyIterationsException ex) {
LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : Too many iterations (%d)", powellOptimizer.getIterations());
} catch (final ConvergenceException ex) {
LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=1) : %s", ex.getMessage());
}
return null;
}
final MixedJumpDistanceFunction function = new MixedJumpDistanceFunction(jumpDistances, estimatedD, n);
final 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();
LoggerUtils.log(logger, Level.FINE, "Powell optimiser fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
} catch (final TooManyEvaluationsException ex) {
LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many evaluation (%d)", n, powellOptimizer.getEvaluations());
} catch (final TooManyIterationsException ex) {
LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : Too many iterations (%d)", n, powellOptimizer.getIterations());
} catch (final ConvergenceException ex) {
LoggerUtils.log(logger, Level.INFO, "Powell optimiser failed to fit (N=%d) : %s", n, ex.getMessage());
}
if (constrainedSolution == null) {
LoggerUtils.log(logger, Level.INFO, "Trying CMAES optimiser with restarts ...");
final double[] uB = function.getUpperBounds();
final 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.
final double[] s = new double[lB.length];
for (int i = 0; i < s.length; i++) {
s[i] = (uB[i] - lB[i]) / 3;
}
final OptimizationData sigma = new CMAESOptimizer.Sigma(s);
final OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(function.x.length))));
// Iterate this for stability in the initial guess
final CMAESOptimizer cmaesOptimizer = createCmaesOptimizer();
for (int i = 0; i <= fitRestarts; i++) {
// Try from the initial guess
try {
final 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;
LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%da] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
}
} catch (final TooManyEvaluationsException ex) {
// No solution
}
if (constrainedSolution == null) {
continue;
}
// Try from the current optimum
try {
final 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;
LoggerUtils.log(logger, Level.FINE, "CMAES optimiser [%db] fit (N=%d) : MLE = %f (%d evaluations)", i, n, solution.getValue(), evaluations);
}
} catch (final TooManyEvaluationsException ex) {
// No solution
}
}
if (constrainedSolution != null) {
try {
// Re-optimise with Powell?
final 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;
LoggerUtils.log(logger, Level.INFO, "Powell optimiser re-fit (N=%d) : MLE = %f (%d evaluations)", n, constrainedSolution.getValue(), powellOptimizer.getEvaluations());
}
} catch (final TooManyEvaluationsException ex) {
// No solution
} catch (final TooManyIterationsException ex) {
// No solution
} catch (final ConvergenceException ex) {
// No solution
}
}
}
if (constrainedSolution == null) {
LoggerUtils.log(logger, Level.INFO, "Failed to fit N=%d", n);
return null;
}
final 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
fitValue = MathUtils.getAkaikeInformationCriterion(ll, fitParams.length - 1);
final double[] d = new double[n];
final 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);
final double[] coefficients = d;
final double[] fractions = f;
LoggerUtils.log(logger, Level.INFO, "Fit Jump distance (N=%d) : %s (%s), MLE = %s, Akaike IC = %s (%d evaluations)", n, formatD(d), format(f), MathUtils.rounded(ll, 4), MathUtils.rounded(fitValue, 4), evaluations);
if (isValid(d, f)) {
lastFitValue = fitValue;
return new double[][] { coefficients, fractions };
}
return null;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class GaussianFit method fitSingle.
/**
* Fits a 2D Gaussian to the given data. Fits all the specified peaks.
*
* <p>Data must be arranged in yx block order, i.e. height rows of width.
*
* @param gf the gf
* @param data the data
* @param width the width
* @param height the height
* @param index Index of the data to fit
* @param estimatedHeight Estimated height for the peak (input from smoothed data)
* @return Array containing the fitted curve data: The first value is the Background. The
* remaining values are Amplitude, PosX, PosY, StdDevX, StdDevY for each fitted peak. Null
* if no fit is possible.
*/
@Nullable
private double[] fitSingle(Gaussian2DFitter gf, float[] data, int width, int height, int index, double estimatedHeight) {
this.fitResult = gf.fit(SimpleArrayUtils.toDouble(data), width, height, new int[] { index }, new double[] { estimatedHeight });
if (fitResult.getStatus() == FitStatus.OK) {
chiSquared = fitResult.getError();
final double[] params = fitResult.getParameters();
convertParameters(params);
// Check the fit is within the data
if (params[Gaussian2DFunction.X_POSITION] < 0 || params[Gaussian2DFunction.X_POSITION] > width || params[Gaussian2DFunction.Y_POSITION] < 0 || params[Gaussian2DFunction.Y_POSITION] > height) {
fitResult = new FitResult(FitStatus.OUTSIDE_FIT_REGION, fitResult.getDegreesOfFreedom(), fitResult.getError(), fitResult.getInitialParameters(), fitResult.getParameters(), fitResult.getParameterDeviations(), fitResult.getNumberOfPeaks(), fitResult.getNumberOfFittedParameters(), fitResult.getStatusData(), fitResult.getIterations(), fitResult.getEvaluations());
return null;
}
return params;
}
return null;
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class ImageJ3DResultsViewer method createSphereSizeFromPrecision.
@Nullable
private static Point3f[] createSphereSizeFromPrecision(MemoryPeakResults results) {
final PrecisionResultProcedure p = new PrecisionResultProcedure(results);
try {
final PrecisionMethod m = p.getPrecision();
IJ.log("Using precision method " + FitProtosHelper.getName(m));
final Point3f[] size = new Point3f[results.size()];
for (int i = 0, j = 0; i < p.precisions.length; i++) {
// Precision is in NM which matches the rendering
final float v = (float) p.precisions[i];
size[j++] = new Point3f(v, v, v);
}
return size;
} catch (final DataException ex) {
IJ.error(TITLE, "The results have no precision: " + ex.getMessage());
return null;
}
}
use of uk.ac.sussex.gdsc.core.annotation.Nullable in project GDSC-SMLM by aherbert.
the class ImageJ3DResultsViewer method createAlphaFromSize.
@Nullable
private static float[] createAlphaFromSize(double minA, double maxA, Point3f[] sphereSize) {
if (sphereSize == null) {
return null;
}
final double[] d = new double[sphereSize.length];
for (int i = 0; i < d.length; i++) {
final Point3f p = sphereSize[i];
if (p.x == p.y && p.y == p.z) {
d[i] = 3.0 * p.x * p.x;
} else {
// Use the squared distance. This is the equivalent of the area of the shape projected to
// 2D.
d[i] = (double) p.x * p.x + p.y * p.y + p.z * p.z;
}
// Use the average radius. This is the equivalent of the mean radius of an enclosing
// ellipsoid.
// d[i] = Math.sqrt(d[i] / 3);
}
final double[] limits = MathUtils.limits(d);
final double min = limits[0];
final double max = limits[1];
if (min == max) {
ImageJUtils.log("No per-item transparency as size is fixed");
return null;
}
final double range = (maxA - minA) / (max - min);
final float[] alpha = new float[d.length];
for (int i = 0; i < alpha.length; i++) {
// Largest distance has lowest alpha (more transparent)
alpha[i] = (float) (minA + range * (max - d[i]));
}
return alpha;
}
Aggregations