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

Example 26 with Statistics

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

the class ErfGaussian2DFunctionTest method functionComputesSecondTargetGradient.

private void functionComputesSecondTargetGradient(int targetParameter) {
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final double[] dyda = new double[f1.getNumberOfGradients()];
    final double[] dyda2 = new double[dyda.length];
    double[] params;
    // Test fitting of second derivatives
    final ErfGaussian2DFunction f1a = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    final ErfGaussian2DFunction f1b = (ErfGaussian2DFunction) 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.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;
                                f1a.initialise1(params.clone());
                                params[targetParameter] = xx - h;
                                f1b.initialise1(params.clone());
                                for (final int x : testx) {
                                    for (final int y : testy) {
                                        final int i = y * maxx + x;
                                        f1a.eval(i, dyda);
                                        final double value2 = dyda[gradientIndex];
                                        f1b.eval(i, dyda);
                                        final double value3 = dyda[gradientIndex];
                                        f1.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 %s (error %s +/- %s)", f1.getClass().getSimpleName(), Gaussian2DFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 27 with Statistics

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

the class TraceLengthAnalysis method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Load the results
    MemoryPeakResults results = ResultsManager.loadInputResults(settings.inputOption, false, null, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        return;
    }
    try {
        distanceConverter = results.getDistanceConverter(DistanceUnit.UM);
        timeConverter = results.getTimeConverter(TimeUnit.SECOND);
    } catch (final Exception ex) {
        IJ.error(TITLE, "Cannot convert units to um or seconds: " + ex.getMessage());
        return;
    }
    // Get the localisation error (4s^2) in raw units^2
    double precision = 0;
    try {
        final PrecisionResultProcedure p = new PrecisionResultProcedure(results);
        p.getPrecision();
        // Precision in nm using the median
        precision = new Percentile().evaluate(p.precisions, 50);
        // Convert from nm to um to raw units
        final double rawPrecision = distanceConverter.convertBack(precision / 1e3);
        // Get the localisation error (4s^2) in units^2
        error = 4 * rawPrecision * rawPrecision;
    } catch (final Exception ex) {
        ImageJUtils.log(TITLE + " - Unable to compute precision: " + ex.getMessage());
    }
    // Analyse the track lengths
    results = results.copy();
    results.sort(IdFramePeakResultComparator.INSTANCE);
    // Ensure the first result triggers an id change
    lastid = results.getFirst().getId() - 1;
    results.forEach(this::processTrackLength);
    // For the final track
    store();
    msds = msdList.toArray();
    lengths = lengthList.toArray();
    ids = idList.toArray();
    final int[] limits = MathUtils.limits(lengths);
    h1 = new int[limits[1] + 1];
    h2 = new int[h1.length];
    x1 = SimpleArrayUtils.newArray(h1.length, 0, 1f);
    y1 = new float[x1.length];
    y2 = new float[x1.length];
    // Sort by MSD
    final int[] indices = SimpleArrayUtils.natural(msds.length);
    SortUtils.sortIndices(indices, msds, false);
    final double[] msds2 = msds.clone();
    final int[] lengths2 = lengths.clone();
    final int[] ids2 = ids.clone();
    for (int i = 0; i < indices.length; i++) {
        msds[i] = msds2[indices[i]];
        lengths[i] = lengths2[indices[i]];
        ids[i] = ids2[indices[i]];
    }
    // Interactive analysis
    final NonBlockingExtendedGenericDialog gd = new NonBlockingExtendedGenericDialog(TITLE);
    ImageJUtils.addMessage(gd, "Split traces into fixed or moving using the track diffusion coefficient (D).\n" + "Localisation error has been subtracted from jumps (%s nm).", MathUtils.rounded(precision));
    final Statistics s = Statistics.create(msds);
    final double av = s.getMean();
    final String msg = String.format("Average D per track = %s um^2/s", MathUtils.rounded(av));
    gd.addMessage(msg);
    // Histogram the diffusion coefficients
    final WindowOrganiser wo = new WindowOrganiser();
    final HistogramPlot histogramPlot = new HistogramPlotBuilder("Trace diffusion coefficient", StoredData.create(msds), "D (um^2/s)").setRemoveOutliersOption(1).setPlotLabel(msg).build();
    histogramPlot.show(wo);
    final double[] xvalues = histogramPlot.getPlotXValues();
    final double min = xvalues[0];
    final double max = xvalues[xvalues.length - 1];
    // see if we can build a nice slider range from the histogram limits
    if (max - min < 5) {
        // Because sliders are used when the range is <5 and floating point
        gd.addSlider("D_threshold", min, max, settings.msdThreshold);
    } else {
        gd.addNumericField("D_threshold", settings.msdThreshold, 2, 6, "um^2/s");
    }
    gd.addCheckbox("Normalise", settings.normalise);
    gd.addDialogListener((gd1, event) -> {
        settings.msdThreshold = gd1.getNextNumber();
        settings.normalise = gd1.getNextBoolean();
        update();
        return true;
    });
    if (ImageJUtils.isShowGenericDialog()) {
        draw(wo);
        wo.tile();
    }
    gd.setOKLabel("Save datasets");
    gd.setCancelLabel("Close");
    gd.addHelp(HelpUrls.getUrl("trace-length-analysis"));
    gd.showDialog();
    if (gd.wasCanceled()) {
        return;
    }
    // Sort by ID
    final PeakResult[] list = results.toArray();
    Arrays.sort(list, IdFramePeakResultComparator.INSTANCE);
    createResults(results, "Fixed", 0, lastIndex, list);
    createResults(results, "Moving", lastIndex, msds.length, list);
}

Example 28 with Statistics

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

the class TraceMolecules method summarise.

private void summarise(Consumer<String> output, Trace[] traces, int filtered, double distanceThreshold, double timeThreshold) {
    IJ.showStatus("Calculating summary ...");
    final Statistics[] stats = new Statistics[NAMES.length];
    for (int i = 0; i < stats.length; i++) {
        stats[i] = (settings.getShowHistograms() || settings.getSaveTraceData()) ? new StoredDataStatistics() : new Statistics();
    }
    int singles = 0;
    for (final Trace trace : traces) {
        final int nBlinks = trace.getBlinks() - 1;
        stats[BLINKS].add(nBlinks);
        final int[] onTimes = trace.getOnTimes();
        final int[] offTimes = trace.getOffTimes();
        double timeOn = 0;
        for (final int t : onTimes) {
            stats[T_ON].add(t * exposureTime);
            timeOn += t * exposureTime;
        }
        stats[TOTAL_T_ON].add(timeOn);
        if (offTimes != null) {
            double timeOff = 0;
            for (final int t : offTimes) {
                stats[T_OFF].add(t * exposureTime);
                timeOff += t * exposureTime;
            }
            stats[TOTAL_T_OFF].add(timeOff);
        }
        final double signal = trace.getSignal() / results.getGain();
        stats[TOTAL_SIGNAL].add(signal);
        stats[SIGNAL_PER_FRAME].add(signal / trace.size());
        stats[DWELL_TIME].add((trace.getTail().getEndFrame() - trace.getHead().getFrame() + 1) * exposureTime);
        if (trace.size() == 1) {
            singles++;
        }
    }
    // Add to the summary table
    final StringBuilder sb = new StringBuilder();
    sb.append(results.getName()).append('\t');
    sb.append(outputName.equals("Cluster") ? getClusteringAlgorithm(settings.getClusteringAlgorithm()) : getTraceMode(settings.getTraceMode())).append('\t');
    sb.append(MathUtils.rounded(getExposureTimeInMilliSeconds(), 3)).append('\t');
    sb.append(MathUtils.rounded(distanceThreshold, 3)).append('\t');
    sb.append(MathUtils.rounded(timeThreshold, 3));
    if (settings.getSplitPulses()) {
        sb.append(" *");
    }
    sb.append('\t');
    sb.append(convertSecondsToFrames(timeThreshold)).append('\t');
    sb.append(traces.length).append('\t');
    sb.append(filtered).append('\t');
    sb.append(singles).append('\t');
    sb.append(traces.length - singles).append('\t');
    for (int i = 0; i < stats.length; i++) {
        sb.append(MathUtils.rounded(stats[i].getMean(), 3)).append('\t');
    }
    if (java.awt.GraphicsEnvironment.isHeadless()) {
        IJ.log(sb.toString());
        return;
    }
    output.accept(sb.toString());
    if (settings.getShowHistograms()) {
        IJ.showStatus("Calculating histograms ...");
        final WindowOrganiser windowOrganiser = new WindowOrganiser();
        final HistogramPlotBuilder builder = new HistogramPlotBuilder(pluginTitle).setNumberOfBins(settings.getHistogramBins());
        for (int i = 0; i < NAMES.length; i++) {
            if (pluginSettings.displayHistograms[i]) {
                builder.setData((StoredDataStatistics) stats[i]).setName(NAMES[i]).setIntegerBins(integerDisplay[i]).setRemoveOutliersOption((settings.getRemoveOutliers() || alwaysRemoveOutliers[i]) ? 2 : 0).show(windowOrganiser);
            }
        }
        windowOrganiser.tile();
    }
    if (settings.getSaveTraceData()) {
        saveTraceData(stats);
    }
    IJ.showStatus("");
}

Example 29 with Statistics

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

the class TraceDiffusion method summarise.

private void summarise(Trace[] traces, double[] fitMsdResult, int n, double[][] jdParams) {
    IJ.showStatus("Calculating summary ...");
    final Statistics[] stats = new Statistics[Settings.NAMES.length];
    for (int i = 0; i < stats.length; i++) {
        stats[i] = (clusteringSettings.getShowHistograms()) ? new StoredDataStatistics() : new Statistics();
    }
    for (final Trace trace : traces) {
        stats[Settings.T_ON].add(trace.getOnTime() * exposureTime);
        final double signal = trace.getSignal() / results.getGain();
        stats[Settings.TOTAL_SIGNAL].add(signal);
        stats[Settings.SIGNAL_PER_FRAME].add(signal / trace.size());
    }
    // Add to the summary table
    final StringBuilder sb = new StringBuilder(settings.title);
    sb.append('\t').append(createCombinedName());
    sb.append('\t');
    sb.append(MathUtils.rounded(exposureTime * 1000, 3)).append('\t');
    appendClusteringSettings(sb).append('\t');
    sb.append(clusteringSettings.getMinimumTraceLength()).append('\t');
    sb.append(clusteringSettings.getIgnoreEnds()).append('\t');
    sb.append(clusteringSettings.getTruncate()).append('\t');
    sb.append(clusteringSettings.getInternalDistances()).append('\t');
    sb.append(clusteringSettings.getFitLength()).append('\t');
    sb.append(clusteringSettings.getMsdCorrection()).append('\t');
    sb.append(clusteringSettings.getPrecisionCorrection()).append('\t');
    sb.append(clusteringSettings.getMle()).append('\t');
    sb.append(traces.length).append('\t');
    sb.append(MathUtils.rounded(precision, 4)).append('\t');
    // D
    double diffCoeff = 0;
    double precision = 0;
    if (fitMsdResult != null) {
        diffCoeff = fitMsdResult[0];
        precision = fitMsdResult[1];
    }
    sb.append(MathUtils.rounded(diffCoeff, 4)).append('\t');
    sb.append(MathUtils.rounded(precision * 1000, 4)).append('\t');
    sb.append(MathUtils.rounded(clusteringSettings.getJumpDistance() * exposureTime)).append('\t');
    sb.append(n).append('\t');
    sb.append(MathUtils.rounded(beta, 4)).append('\t');
    if (jdParams == null) {
        sb.append("\t\t\t");
    } else {
        sb.append(format(jdParams[0])).append('\t');
        sb.append(format(jdParams[1])).append('\t');
        sb.append(MathUtils.rounded(fitValue)).append('\t');
    }
    for (int i = 0; i < stats.length; i++) {
        sb.append(MathUtils.rounded(stats[i].getMean(), 3)).append('\t');
    }
    createSummaryTable().accept(sb.toString());
    if (java.awt.GraphicsEnvironment.isHeadless()) {
        return;
    }
    if (clusteringSettings.getShowHistograms()) {
        IJ.showStatus("Calculating histograms ...");
        for (int i = 0; i < Settings.NAMES.length; i++) {
            if (settings.displayHistograms[i]) {
                showHistogram((StoredDataStatistics) stats[i], Settings.NAMES[i], settings.alwaysRemoveOutliers[i], Settings.ROUNDED[i], false);
            }
        }
    }
    windowOrganiser.tile();
    IJ.showStatus("Finished " + TITLE);
}

Example 30 with Statistics

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

the class DynamicMultipleTargetTracing method createMatchedAction.

/**
 * Creates the action to process a connection between a trajectory and a localisation.
 *
 * <p>Currently the temporal window is occurrences and not an actual time span. Thus some local
 * statistics can be averaged over a long duration. However if the trajectory is for a blinking
 * particle with short on-times and long off-times the local statistics are no better or worse
 * than the global model.
 *
 * @param configuration the configuration
 * @return the matched action
 */
private BiConsumer<Trajectory, PeakResult> createMatchedAction(DmttConfiguration configuration) {
    // Optionally disable the local diffusion model
    final Consumer<Trajectory> updateLocalDiffusion;
    if (configuration.isDisableLocalDiffusionModel()) {
        updateLocalDiffusion = t -> {
        /* Do nothing */
        };
    } else {
        final double dMax = computeDMax(configuration);
        // D_max = r^2 / 4t ('t' is frame acquisition time)
        // dMax has been converted to frames
        final double r2max = dMax * 4;
        updateLocalDiffusion = t -> updateLocalDiffusion(configuration, r2max, t);
    }
    // Optionally disable the intensity model
    if (isDisableIntensityModel(configuration)) {
        return (t, r) -> {
            // Add the peak to the trajectory. Do not track on-frames for intensity.
            t.add(r);
            updateLocalDiffusion.accept(t);
        // Ignore local intensity
        };
    }
    final Statistics stats = new Statistics();
    return (t, r) -> {
        // Add the peak to the trajectory tracking on-frames
        final double intensity = r.getIntensity();
        final boolean on = t.isLocalIntensity ? isOn(intensity, t.meanI, t.sdI) : isOn(intensity, meanI, sdI);
        t.add(r, on);
        updateLocalDiffusion.accept(t);
        updateLocalIntensity(configuration, stats, t);
    };
}