Examples with Statistics uk.ac.sussex.gdsc.core.utils.Statistics used on opensource projects

Example 31 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionComputesSecondTargetGradientWith2Peaks.

private void functionComputesSecondTargetGradientWith2Peaks(int targetParameter) {
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) this.f2;
    final int gradientIndex = findGradientIndex(f2, targetParameter);
    final double[] dyda = new double[f2.getNumberOfGradients()];
    final double[] dyda2 = new double[dyda.length];
    double[] params;
    // Test fitting of second derivatives
    final ErfGaussian2DFunction f2a = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, maxx, maxy, flags, zModel);
    final ErfGaussian2DFunction f2b = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, maxx, maxy, flags, zModel);
    final Statistics s = new Statistics();
    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) {
                        for (final double[] w1 : testw1) {
                            for (final double angle1 : testangle1) {
                                // Peak 2
                                for (final double signal2 : testsignal2) {
                                    for (final double cx2 : testcx2) {
                                        for (final double cy2 : testcy2) {
                                            for (final double cz2 : testcz2) {
                                                for (final double[] w2 : testw2) {
                                                    for (final double angle2 : testangle2) {
                                                        params = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1, signal2, cx2, cy2, cz2, w2[0], w2[1], angle2);
                                                        f2.initialise2(params);
                                                        // Numerically solve gradient.
                                                        // Calculate the step size h to be an exact numerical representation
                                                        final double xx = params[targetParameter];
                                                        // Get h to minimise roundoff error
                                                        final double h = Precision.representableDelta(xx, stepH);
                                                        // Evaluate at (x+h) and (x-h)
                                                        params[targetParameter] = xx + h;
                                                        f2a.initialise1(params.clone());
                                                        params[targetParameter] = xx - h;
                                                        f2b.initialise1(params.clone());
                                                        for (final int x : testx) {
                                                            for (final int y : testy) {
                                                                final int i = y * maxx + x;
                                                                f2a.eval(i, dyda);
                                                                final double value2 = dyda[gradientIndex];
                                                                f2b.eval(i, dyda);
                                                                final double value3 = dyda[gradientIndex];
                                                                f2.eval2(i, dyda, dyda2);
                                                                final double gradient = (value2 - value3) / (2 * h);
                                                                final double error = DoubleEquality.relativeError(gradient, dyda2[gradientIndex]);
                                                                s.add(error);
                                                                Assertions.assertTrue((gradient * dyda2[gradientIndex]) >= 0, () -> gradient + " sign != " + dyda2[gradientIndex]);
                                                                // TestLog.fine(logger,"[%d,%d] %f == [%d] %f? (%g)", x, y,
                                                                // gradient,
                                                                // gradientIndex, dyda2[gradientIndex], error);
                                                                Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(gradient, dyda2[gradientIndex]), () -> gradient + " != " + dyda2[gradientIndex]);
                                                            }
                                                        }
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesSecondTargetGradient %s [%d] %s (error %s +/- %s)", f2.getClass().getSimpleName(), Gaussian2DFunction.getPeak(targetParameter), Gaussian2DFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 32 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class CubicSplineFunctionTest method functionComputesTargetGradient2.

private void functionComputesTargetGradient2(int targetParameter) {
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final Statistics s = new Statistics();
    final StandardGradient1Procedure p1a = new StandardGradient1Procedure();
    final StandardGradient1Procedure p1b = new StandardGradient1Procedure();
    final StandardGradient2Procedure p2 = new StandardGradient2Procedure();
    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) {
                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1);
                        // System.out.println(java.util.Arrays.toString(a));
                        f1.initialise2(a);
                        final boolean test = !f1.isNodeBoundary(gradientIndex);
                        // Comment out when printing errors
                        if (!test) {
                            continue;
                        }
                        // Evaluate all gradients
                        p2.getValues(f1, 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(f1, a);
                        a[targetParameter] = xx - h;
                        p1b.getValues(f1, 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.gradients[i][gradientIndex];
                                final double low = p1b.gradients[i][gradientIndex];
                                final double gradient = (high - low) / (2 * h);
                                final double d2yda2 = p2.gradients2[i][gradientIndex];
                                final double error = DoubleEquality.relativeError(gradient, d2yda2);
                                // gradient, gradientIndex, d2yda2, error);
                                if (test) {
                                    s.add(error);
                                    Assertions.assertTrue((gradient * d2yda2) >= 0, () -> String.format("%s sign != %s", gradient, d2yda2));
                                    // logger.fine(FunctionUtils.getSupplier("[%d,%d] %f == [%d] %f? (%g)", x, y,
                                    // gradient, gradientIndex, d2yda2, error);
                                    Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(gradient, d2yda2), () -> String.format("%s != %s", gradient, d2yda2));
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradient2 %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), CubicSplineFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 33 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class Gaussian2DFunctionTest method functionComputesTargetGradient.

private void functionComputesTargetGradient(int targetParameter) {
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final double[] dyda = new double[f1.gradientIndices().length];
    final double[] dyda2 = new double[dyda.length];
    double[] params;
    final Gaussian2DFunction f1a = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    final Gaussian2DFunction f1b = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    final Statistics s = new Statistics();
    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) {
                        for (final double[] w1 : testw1) {
                            for (final double angle1 : testangle1) {
                                params = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1);
                                f1.initialise(params);
                                // Numerically solve gradient.
                                // Calculate the step size h to be an exact numerical representation
                                final double xx = params[targetParameter];
                                // Get h to minimise roundoff error
                                final double h = Precision.representableDelta(xx, stepH);
                                // Evaluate at (x+h) and (x-h)
                                params[targetParameter] = xx + h;
                                f1a.initialise(params.clone());
                                params[targetParameter] = xx - h;
                                f1b.initialise(params.clone());
                                for (final int x : testx) {
                                    for (final int y : testy) {
                                        final int i = y * maxx + x;
                                        f1.eval(i, dyda);
                                        final double value2 = f1a.eval(i, dyda2);
                                        final double value3 = f1b.eval(i, dyda2);
                                        final double gradient = (value2 - value3) / (2 * h);
                                        final double error = DoubleEquality.relativeError(gradient, dyda2[gradientIndex]);
                                        s.add(error);
                                        if ((gradient * dyda2[gradientIndex]) < 0) {
                                            Assertions.fail(String.format("%s sign != %s", gradient, dyda2[gradientIndex]));
                                        }
                                        // gradient, gradientIndex, dyda[gradientIndex]);
                                        if (!eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex])) {
                                            Assertions.fail(String.format("%s != %s", gradient, dyda[gradientIndex]));
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradient %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), Gaussian2DFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 34 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class CmosAnalysis method computeError.

private void computeError(int slice, ImageStack simulationStack, WindowOrganiser wo) {
    final String label = simulationStack.getSliceLabel(slice);
    final float[] e = (float[]) simulationStack.getPixels(slice);
    final float[] o = (float[]) measuredStack.getPixels(slice);
    // Get the mean error
    final double[] error = new double[e.length];
    for (int i = e.length; i-- > 0; ) {
        error[i] = (double) o[i] - e[i];
    }
    final Statistics s = new Statistics();
    s.add(error);
    final StringBuilder result = new StringBuilder("Error ").append(label);
    result.append(" = ").append(MathUtils.rounded(s.getMean()));
    result.append(" +/- ").append(MathUtils.rounded(s.getStandardDeviation()));
    // Do statistical tests
    final double[] x = SimpleArrayUtils.toDouble(e);
    final double[] y = SimpleArrayUtils.toDouble(o);
    final PearsonsCorrelation c = new PearsonsCorrelation();
    result.append(" : R=").append(MathUtils.rounded(c.correlation(x, y)));
    // Plot these
    String title = TITLE + " " + label + " Simulation vs Measured";
    final Plot plot = new Plot(title, "simulated", "measured");
    plot.addPoints(e, o, Plot.DOT);
    plot.addLabel(0, 0, result.toString());
    ImageJUtils.display(title, plot, wo);
    // Histogram the error
    new HistogramPlotBuilder(TITLE + " " + label, DoubleData.wrap(error), "Error").setPlotLabel(result.toString()).show(wo);
    // Kolmogorov–Smirnov test that the distributions are the same
    double pvalue = TestUtils.kolmogorovSmirnovTest(x, y);
    result.append(" : Kolmogorov–Smirnov p=").append(MathUtils.rounded(pvalue)).append(' ').append(((pvalue < 0.001) ? REJECT : ACCEPT));
    if (slice == 3) {
        // Paired T-Test compares two related samples to assess whether their
        // population means differ.
        // T-Test is valid when the difference between the means is normally
        // distributed, e.g. gain
        pvalue = TestUtils.pairedTTest(x, y);
        result.append(" : Paired T-Test p=").append(MathUtils.rounded(pvalue)).append(' ').append(((pvalue < 0.001) ? REJECT : ACCEPT));
    } else {
        // Wilcoxon Signed Rank test compares two related samples to assess whether their
        // population mean ranks differ
        final WilcoxonSignedRankTest wsrTest = new WilcoxonSignedRankTest();
        pvalue = wsrTest.wilcoxonSignedRankTest(x, y, false);
        result.append(" : Wilcoxon Signed Rank p=").append(MathUtils.rounded(pvalue)).append(' ').append(((pvalue < 0.001) ? REJECT : ACCEPT));
    }
    ImageJUtils.log(result.toString());
}

Example 35 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.

the class DensityEstimator method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    // Require some fit results and selected regions
    if (MemoryPeakResults.countMemorySize() == 0) {
        IJ.error(TITLE, "There are no fitting results in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Currently this only supports pixel distance units
    final MemoryPeakResults results = ResultsManager.loadInputResults(settings.inputOption, false, DistanceUnit.PIXEL, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        IJ.showStatus("");
        return;
    }
    final long start = System.currentTimeMillis();
    IJ.showStatus("Calculating density ...");
    // Scale to um^2 from px^2
    final double scale = Math.pow(results.getDistanceConverter(DistanceUnit.UM).convertBack(1), 2);
    results.sort();
    final FrameCounter counter = results.newFrameCounter();
    final double localisationsPerFrame = (double) results.size() / (results.getLastFrame() - counter.currentFrame() + 1);
    final Rectangle bounds = results.getBounds(true);
    final double globalDensity = localisationsPerFrame / bounds.width / bounds.height;
    final int border = settings.border;
    final boolean includeSingles = settings.includeSingles;
    final int size = 2 * border + 1;
    final double minDensity = Math.pow(size, -2);
    ImageJUtils.log("%s : %s : Global density %s. Minimum density in %dx%d px = %s um^-2", TITLE, results.getName(), MathUtils.rounded(globalDensity * scale), size, size, MathUtils.rounded(minDensity * scale));
    final TIntArrayList x = new TIntArrayList();
    final TIntArrayList y = new TIntArrayList();
    final ExecutorService es = Executors.newFixedThreadPool(Prefs.getThreads());
    final LocalList<FrameDensity> densities = new LocalList<>();
    final LocalList<Future<?>> futures = new LocalList<>();
    results.forEach((PeakResultProcedure) (peak) -> {
        if (counter.advance(peak.getFrame())) {
            final FrameDensity fd = new FrameDensity(peak.getFrame(), x.toArray(), y.toArray(), border, includeSingles);
            densities.add(fd);
            futures.add(es.submit(fd));
            x.resetQuick();
            y.resetQuick();
        }
        x.add((int) peak.getXPosition());
        y.add((int) peak.getYPosition());
    });
    densities.add(new FrameDensity(counter.currentFrame(), x.toArray(), y.toArray(), border, includeSingles));
    futures.add(es.submit(densities.get(densities.size() - 1)));
    es.shutdown();
    // Wait
    ConcurrencyUtils.waitForCompletionUnchecked(futures);
    densities.sort((o1, o2) -> Integer.compare(o1.frame, o2.frame));
    final int total = densities.stream().mapToInt(fd -> fd.counts.length).sum();
    // Plot density
    final Statistics stats = new Statistics();
    final float[] frame = new float[total];
    final float[] density = new float[total];
    densities.stream().forEach(fd -> {
        for (int i = 0; i < fd.counts.length; i++) {
            final double d = (fd.counts[i] / fd.values[i]) * scale;
            frame[stats.getN()] = fd.frame;
            density[stats.getN()] = (float) d;
            stats.add(d);
        }
    });
    final double mean = stats.getMean();
    final double sd = stats.getStandardDeviation();
    final String label = String.format("Density = %s +/- %s um^-2", MathUtils.rounded(mean), MathUtils.rounded(sd));
    final Plot plot = new Plot("Frame vs Density", "Frame", "Density (um^-2)");
    plot.addPoints(frame, density, Plot.CIRCLE);
    plot.addLabel(0, 0, label);
    final WindowOrganiser wo = new WindowOrganiser();
    ImageJUtils.display(plot.getTitle(), plot, wo);
    // Histogram density
    new HistogramPlotBuilder("Local", StoredData.create(density), "Density (um^-2)").setPlotLabel(label).show(wo);
    wo.tile();
    // Log the number of singles
    final int singles = densities.stream().mapToInt(fd -> fd.singles).sum();
    ImageJUtils.log("Singles %d / %d (%s%%)", singles, results.size(), MathUtils.rounded(100.0 * singles / results.size()));
    IJ.showStatus(TITLE + " complete : " + TextUtils.millisToString(System.currentTimeMillis() - start));
}