Search in sources :

Example 16 with Statistics

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());
    }
}
Also used : StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics)

Example 17 with Statistics

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;
}
Also used : IJImageSource(uk.ac.sussex.gdsc.smlm.ij.IJImageSource) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) Point(java.awt.Point)

Example 18 with Statistics

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;
}
Also used : Trace(uk.ac.sussex.gdsc.smlm.results.Trace) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) PeakResult(uk.ac.sussex.gdsc.smlm.results.PeakResult)

Example 19 with Statistics

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;
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) SingleEllipticalGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.SingleEllipticalGaussian2DFunction) ArrayList(java.util.ArrayList) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) SingleFreeCircularGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.SingleFreeCircularGaussian2DFunction) SingleEllipticalGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.SingleEllipticalGaussian2DFunction) SingleNbFixedGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.SingleNbFixedGaussian2DFunction) EllipticalGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction) SingleFixedGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.SingleFixedGaussian2DFunction) Gaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction) SingleCircularGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.SingleCircularGaussian2DFunction) EjmlLinearSolver(uk.ac.sussex.gdsc.smlm.fitting.linear.EjmlLinearSolver) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 20 with Statistics

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()));
    });
}
Also used : StandardGradient1Procedure(uk.ac.sussex.gdsc.smlm.function.StandardGradient1Procedure) StandardValueProcedure(uk.ac.sussex.gdsc.smlm.function.StandardValueProcedure) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics)

Aggregations

Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)46 StoredDataStatistics (uk.ac.sussex.gdsc.core.utils.StoredDataStatistics)16 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)11 Rectangle (java.awt.Rectangle)10 ArrayList (java.util.ArrayList)10 WindowOrganiser (uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser)10 LocalList (uk.ac.sussex.gdsc.core.utils.LocalList)10 Plot (ij.gui.Plot)8 PeakResult (uk.ac.sussex.gdsc.smlm.results.PeakResult)8 PeakResultProcedure (uk.ac.sussex.gdsc.smlm.results.procedures.PeakResultProcedure)7 ImagePlus (ij.ImagePlus)6 ExecutorService (java.util.concurrent.ExecutorService)6 Future (java.util.concurrent.Future)6 HistogramPlotBuilder (uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder)6 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)6 Ticker (uk.ac.sussex.gdsc.core.logging.Ticker)6 DistanceUnit (uk.ac.sussex.gdsc.smlm.data.config.UnitProtos.DistanceUnit)6 Trace (uk.ac.sussex.gdsc.smlm.results.Trace)6 TIntArrayList (gnu.trove.list.array.TIntArrayList)5 ImageStack (ij.ImageStack)5