use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class BenchmarkFit method runAnalysis.
private void runAnalysis() {
benchmarkParameters = benchmarkResults.getFirst().benchmarkParameters;
final double sa = getSa();
// The fitting could have used centre-of-mass or not making the number of points different.
// Find the shortest array (this will be the one where the centre-of-mass was not used)
int length = Integer.MAX_VALUE;
for (final BenchmarkResult benchmarkResult : benchmarkResults) {
if (length > benchmarkResult.results.length) {
length = benchmarkResult.results.length;
}
}
// Build a list of all the frames which have results
final int[] valid = new int[length];
int index = 0;
final int[] counts = new int[benchmarkResults.size()];
for (final BenchmarkResult benchmarkResult : benchmarkResults) {
int count = 0;
for (int i = 0; i < valid.length; i++) {
if (benchmarkResult.results[i] != null) {
count++;
valid[i]++;
}
}
counts[index++] = count;
}
final int target = benchmarkResults.size();
// Check that we have data
if (!validData(valid, target)) {
IJ.error(TITLE, "No frames have fitting results from all methods");
return;
}
// Get the number of start points normalised for all the results
final int totalFrames = benchmarkParameters.frames;
final double numberOfStartPointsNormalisation = (double) length / totalFrames;
createAnalysisTable();
// Create the results using only frames where all the fitting methods were successful
index = 0;
for (final BenchmarkResult benchmarkResult : benchmarkResults) {
final double[] answer = benchmarkResult.answer;
final Statistics[] stats = new Statistics[NAMES.length];
for (int i = 0; i < stats.length; i++) {
stats[i] = new Statistics();
}
for (int i = 0; i < valid.length; i++) {
if (valid[i] < target) {
continue;
}
addResult(stats, answer, benchmarkParameters.framePhotons[i % totalFrames], sa, benchmarkResult.results[i], benchmarkResult.resultsTime[i]);
}
final StringBuilder sb = new StringBuilder(benchmarkResult.parameters);
// Now output the actual results ...
sb.append('\t');
final double recall = (stats[0].getN() / numberOfStartPointsNormalisation) / benchmarkParameters.getMolecules();
sb.append(MathUtils.rounded(recall));
// Add the original recall
sb.append('\t');
final double recall2 = (counts[index++] / numberOfStartPointsNormalisation) / benchmarkParameters.getMolecules();
sb.append(MathUtils.rounded(recall2));
// Convert to units of the image (ADUs and pixels)
final double[] convert = benchmarkResult.convert;
for (int i = 0; i < stats.length; i++) {
if (convert[i] != 0) {
sb.append('\t').append(MathUtils.rounded(stats[i].getMean() * convert[i], 6)).append('\t').append(MathUtils.rounded(stats[i].getStandardDeviation() * convert[i]));
} else {
sb.append("\t0\t0");
}
}
analysisTable.append(sb.toString());
}
}
use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class SpotAnalysis method extractSpotProfile.
private double[][] extractSpotProfile(ImagePlus imp, Rectangle bounds, ImageStack rawSpot) {
final int nSlices = imp.getStackSize();
final IJImageSource rawSource = new IJImageSource(imp);
final double[][] profile = new double[2][nSlices];
for (int n = 0; n < nSlices; n++) {
IJ.showProgress(n, nSlices);
final float[] data = rawSource.next(bounds);
rawSpot.setPixels(data, n + 1);
final Statistics stats = Statistics.create(data);
profile[0][n] = stats.getMean() / gain;
profile[1][n] = stats.getStandardDeviation() / gain;
}
return profile;
}
use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class SpotAnalysis method saveSpot.
private void saveSpot() {
if (onFrames.isEmpty() || !updated) {
return;
}
createResultsWindow();
id++;
double signal = 0;
Trace trace = null;
final float psfWidth = Float.parseFloat(widthTextField.getText());
final float cx = areaBounds.x + areaBounds.width / 2.0f;
final float cy = areaBounds.y + areaBounds.height / 2.0f;
for (final Spot s : onFrames) {
// Get the signal again since the background may have changed.
final double spotSignal = getSignal(s.frame);
signal += spotSignal;
final float[] params = Gaussian2DPeakResultHelper.createOneAxisParams(0, (float) (spotSignal), cx, cy, 0, psfWidth);
final PeakResult result = new PeakResult(s.frame, (int) cx, (int) cy, 0, 0, 0, 0, params, null);
if (trace == null) {
trace = new Trace(result);
} else {
trace.add(result);
}
}
if (trace == null) {
return;
}
final Statistics tOn = Statistics.create(trace.getOnTimes());
final Statistics tOff = Statistics.create(trace.getOffTimes());
resultsWindow.append(String.format("%d\t%.1f\t%.1f\t%s\t%s\t%s\t%d\t%s\t%s\t%s", id, cx, cy, MathUtils.rounded(signal, 4), MathUtils.rounded(tOn.getSum() * msPerFrame, 3), MathUtils.rounded(tOff.getSum() * msPerFrame, 3), trace.getBlinks() - 1, MathUtils.rounded(tOn.getMean() * msPerFrame, 3), MathUtils.rounded(tOff.getMean() * msPerFrame, 3), imp.getTitle()));
// Save the individual on/off times for use in creating a histogram
traces.put(id, trace);
updated = false;
}
use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class GradientCalculatorSpeedTest method gradientCalculatorComputesSameOutputWithBias.
@SeededTest
void gradientCalculatorComputesSameOutputWithBias(RandomSeed seed) {
final Gaussian2DFunction func = new SingleEllipticalGaussian2DFunction(blockWidth, blockWidth);
final int nparams = func.getNumberOfGradients();
final GradientCalculator calc = new GradientCalculator(nparams);
final int n = func.size();
final int iter = 50;
final ArrayList<double[]> paramsList = new ArrayList<>(iter);
final ArrayList<double[]> yList = new ArrayList<>(iter);
final ArrayList<double[][]> alphaList = new ArrayList<>(iter);
final ArrayList<double[]> betaList = new ArrayList<>(iter);
final ArrayList<double[]> xList = new ArrayList<>(iter);
// Manipulate the background
final double defaultBackground = background;
final boolean report = logger.isLoggable(Level.INFO);
try {
background = 1e-2;
createData(RngUtils.create(seed.getSeed()), 1, iter, paramsList, yList, true);
final EjmlLinearSolver solver = new EjmlLinearSolver(1e-5, 1e-6);
for (int i = 0; i < paramsList.size(); i++) {
final double[] y = yList.get(i);
final double[] a = paramsList.get(i);
final double[][] alpha = new double[nparams][nparams];
final double[] beta = new double[nparams];
calc.findLinearised(n, y, a, alpha, beta, func);
alphaList.add(alpha);
betaList.add(beta.clone());
for (int j = 0; j < nparams; j++) {
if (Math.abs(beta[j]) < 1e-6) {
logger.info(FunctionUtils.getSupplier("[%d] Tiny beta %s %g", i, func.getGradientParameterName(j), beta[j]));
}
}
// Solve
if (!solver.solve(alpha, beta)) {
throw new AssertionError();
}
xList.add(beta);
// System.out.println(Arrays.toString(beta));
}
final double[][] alpha = new double[nparams][nparams];
final double[] beta = new double[nparams];
final Statistics[] rel = new Statistics[nparams];
final Statistics[] abs = new Statistics[nparams];
for (int i = 0; i < nparams; i++) {
rel[i] = new Statistics();
abs[i] = new Statistics();
}
final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
// for (double b : new double[] { -500, -100, -10, -1, -0.1, 0.1, 1, 10, 100, 500 })
for (final double b : new double[] { -10, -1, -0.1, 0.1, 1, 10 }) {
if (report) {
for (int i = 0; i < nparams; i++) {
rel[i].reset();
abs[i].reset();
}
}
for (int i = 0; i < paramsList.size(); i++) {
final double[] y = add(yList.get(i), b);
final double[] a = paramsList.get(i).clone();
a[0] += b;
calc.findLinearised(n, y, a, alpha, beta, func);
final double[][] alpha2 = alphaList.get(i);
final double[] beta2 = betaList.get(i);
final double[] x2 = xList.get(i);
TestAssertions.assertArrayTest(beta2, beta, predicate, "Beta");
TestAssertions.assertArrayTest(alpha2, alpha, predicate, "Alpha");
// Solve
solver.solve(alpha, beta);
Assertions.assertArrayEquals(x2, beta, 1e-10, "X");
if (report) {
for (int j = 0; j < nparams; j++) {
rel[j].add(DoubleEquality.relativeError(x2[j], beta[j]));
abs[j].add(Math.abs(x2[j] - beta[j]));
}
}
}
if (report) {
for (int i = 0; i < nparams; i++) {
logger.info(FunctionUtils.getSupplier("Bias = %.2f : %s : Rel %g +/- %g: Abs %g +/- %g", b, func.getGradientParameterName(i), rel[i].getMean(), rel[i].getStandardDeviation(), abs[i].getMean(), abs[i].getStandardDeviation()));
}
}
}
} finally {
background = defaultBackground;
}
}
use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class CubicSplineFunctionTest method functionComputesTargetGradient1With2Peaks.
private void functionComputesTargetGradient1With2Peaks(int targetParameter) {
final int gradientIndex = findGradientIndex(f2, targetParameter);
final Statistics s = new Statistics();
final StandardValueProcedure p1a = new StandardValueProcedure();
final StandardValueProcedure p1b = new StandardValueProcedure();
final StandardGradient1Procedure p2 = new StandardGradient1Procedure();
for (final double background : testbackground) {
// Peak 1
for (final double signal1 : testsignal1) {
for (final double cx1 : testcx1) {
for (final double cy1 : testcy1) {
for (final double cz1 : testcz1) {
// Peak 2
for (final double signal2 : testsignal2) {
for (final double cx2 : testcx2) {
for (final double cy2 : testcy2) {
for (final double cz2 : testcz2) {
final double[] a = createParameters(background, signal1, cx1, cy1, cz1, signal2, cx2, cy2, cz2);
// System.out.println(java.util.Arrays.toString(a));
// Evaluate all gradients
p2.getValues(f2, a);
// 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
final double h = Precision.representableDelta(xx, stepH);
// Evaluate at (x+h) and (x-h)
a[targetParameter] = xx + h;
p1a.getValues(f2, a);
a[targetParameter] = xx - h;
p1b.getValues(f2, a);
// Only test close to the XY centre
for (final int x : testx) {
for (final int y : testy) {
final int i = y * maxx + x;
final double high = p1a.values[i];
final double low = p1b.values[i];
final double gradient = (high - low) / (2 * h);
final double dyda = p2.gradients[i][gradientIndex];
final double error = DoubleEquality.relativeError(gradient, dyda);
s.add(error);
Assertions.assertTrue((gradient * dyda) >= 0, () -> String.format("%s sign != %s", gradient, dyda));
// logger.fine(FunctionUtils.getSupplier("[%d,%d] %f == [%d] %f? (%g)", x,
// y, gradient, gradientIndex, dyda, error);
Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(gradient, dyda), () -> String.format("%s != %s", gradient, dyda));
}
}
}
}
}
}
}
}
}
}
}
logger.info(() -> {
return String.format("functionComputesTargetGradient1With2Peaks %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), CubicSplineFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
});
}
Aggregations