Examples with StoredDataStatistics gdsc.core.utils.StoredDataStatistics used on opensource projects

Example 6 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class DoubletAnalysis method summariseResults.

/**
	 * Summarise results.
	 *
	 * @param results
	 *            the results
	 * @param runTime
	 * @param density2
	 */
private void summariseResults(ArrayList<DoubletResult> results, double density, long runTime) {
    // Store results in memory for later analysis
    doubletResults = results;
    // If we are only assessing results with no neighbour candidates
    // TODO - Count the number of actual results that have no neighbours
    numberOfMolecules = this.results.size() - ignored.get();
    // Store details we want in the analysis table
    StringBuilder sb = new StringBuilder();
    sb.append(Utils.rounded(density)).append("\t");
    sb.append(Utils.rounded(getSa())).append("\t");
    sb.append(config.getRelativeFitting()).append("\t");
    sb.append(fitConfig.getFitFunction().toString());
    sb.append(":").append(PeakFit.getSolverName(fitConfig));
    if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
        sb.append(":Camera\t");
        // Add details of the noise model for the MLE
        sb.append("EM=").append(fitConfig.isEmCCD());
        sb.append(":A=").append(Utils.rounded(fitConfig.getAmplification()));
        sb.append(":N=").append(Utils.rounded(fitConfig.getReadNoise()));
        sb.append("\t");
    } else
        sb.append("\t\t");
    analysisPrefix = sb.toString();
    // -=-=-=-=-
    showResults(results, showResults);
    createSummaryTable();
    sb.setLength(0);
    final int n = countN(results);
    // Create the benchmark settings and the fitting settings
    sb.append(numberOfMolecules).append("\t");
    sb.append(n).append("\t");
    sb.append(Utils.rounded(density)).append("\t");
    sb.append(Utils.rounded(simulationParameters.minSignal)).append("\t");
    sb.append(Utils.rounded(simulationParameters.maxSignal)).append("\t");
    sb.append(Utils.rounded(simulationParameters.signalPerFrame)).append("\t");
    sb.append(Utils.rounded(simulationParameters.s)).append("\t");
    sb.append(Utils.rounded(simulationParameters.a)).append("\t");
    sb.append(Utils.rounded(getSa() * simulationParameters.a)).append("\t");
    sb.append(Utils.rounded(simulationParameters.gain)).append("\t");
    sb.append(Utils.rounded(simulationParameters.readNoise)).append("\t");
    sb.append(Utils.rounded(simulationParameters.b)).append("\t");
    // Compute the noise
    double noise = Math.sqrt((simulationParameters.b * ((simulationParameters.emCCD) ? 2 : 1)) / simulationParameters.gain + simulationParameters.readNoise * simulationParameters.readNoise);
    sb.append(Utils.rounded(noise)).append("\t");
    sb.append(Utils.rounded(simulationParameters.signalPerFrame / noise)).append("\t");
    sb.append(config.getRelativeFitting()).append("\t");
    sb.append(fitConfig.getFitFunction().toString());
    sb.append(":").append(PeakFit.getSolverName(fitConfig));
    if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
        sb.append(":Camera\t");
        // Add details of the noise model for the MLE
        sb.append("EM=").append(fitConfig.isEmCCD());
        sb.append(":A=").append(Utils.rounded(fitConfig.getAmplification()));
        sb.append(":N=").append(Utils.rounded(fitConfig.getReadNoise()));
        sb.append("\t");
    } else
        sb.append("\t\t");
    // Now output the actual results ...
    // Show histograms as cumulative to avoid problems with bin width
    // Residuals scores 
    // Iterations and evaluations where fit was OK
    StoredDataStatistics[] stats = new StoredDataStatistics[NAMES2.length];
    for (int i = 0; i < stats.length; i++) stats[i] = new StoredDataStatistics();
    // For Jaccard scoring we need to count the score with no residuals threshold,
    // i.e. Accumulate the score accepting all doublets that were fit 
    double tp = 0;
    double fp = 0;
    double bestTp = 0, bestFp = 0;
    ArrayList<DoubletBonus> data = new ArrayList<DoubletBonus>(results.size());
    for (DoubletResult result : results) {
        final double score = result.getMaxScore();
        // Filter the singles that would be accepted
        if (result.good1) {
            // Filter the doublets that would be accepted
            if (result.good2) {
                final double tp2 = result.tp2a + result.tp2b;
                final double fp2 = result.fp2a + result.fp2b;
                tp += tp2;
                fp += fp2;
                if (result.tp2a > 0.5) {
                    bestTp += result.tp2a;
                    bestFp += result.fp2a;
                }
                if (result.tp2b > 0.5) {
                    bestTp += result.tp2b;
                    bestFp += result.fp2b;
                }
                // Store this as a doublet bonus
                data.add(new DoubletBonus(score, result.getAvScore(), tp2 - result.tp1, fp2 - result.fp1));
            } else {
                // No doublet fit so this will always be the single fit result
                tp += result.tp1;
                fp += result.fp1;
                if (result.tp1 > 0.5) {
                    bestTp += result.tp1;
                    bestFp += result.fp1;
                }
            }
        }
        // Build statistics
        final int c = result.c;
        // True results, i.e. where there was a choice between single or doublet
        if (result.valid) {
            stats[c].add(score);
        }
        // Of those where the fit was good, summarise the iterations and evaluations
        if (result.good1) {
            stats[3].add(result.iter1);
            stats[4].add(result.eval1);
            // about the iteration increase for singles that are not doublets.
            if (c != 0 && result.good2) {
                stats[5].add(result.iter2);
                stats[6].add(result.eval2);
            }
        }
    }
    // Debug the counts
    //		double tpSingle = 0;
    //		double fpSingle = 0;
    //		double tpDoublet = 0;
    //		double fpDoublet = 0;
    //		int nSingle = 0, nDoublet = 0;
    //		for (DoubletResult result : results)
    //		{
    //			if (result.good1)
    //			{
    //				if (result.good2)
    //				{
    //					tpDoublet += result.tp2;
    //					fpDoublet += result.fp2;
    //					nDoublet++;
    //				}
    //				tpSingle += result.tp1;
    //				fpSingle += result.fp1;
    //				nSingle++;
    //			}
    //		}
    //		System.out.printf("Single %.1f,%.1f (%d) : Doublet %.1f,%.1f (%d)\n", tpSingle, fpSingle, nSingle, tpDoublet, fpDoublet, nDoublet*2);
    // Summarise score for true results
    Percentile p = new Percentile(99);
    for (int c = 0; c < stats.length; c++) {
        double[] values = stats[c].getValues();
        // Sorting is need for the percentile and the cumulative histogram so do it once 
        Arrays.sort(values);
        sb.append(Utils.rounded(stats[c].getMean())).append("+/-").append(Utils.rounded(stats[c].getStandardDeviation())).append(" (").append(stats[c].getN()).append(") ").append(Utils.rounded(p.evaluate(values))).append('\t');
        if (showHistograms && displayHistograms[c + NAMES.length])
            showHistogram(values, NAMES2[c]);
    }
    sb.append(MATCHING[matching]).append('\t');
    // Plot a graph of the additional results we would fit at all score thresholds.
    // This assumes we just pick the the doublet if we fit it (NO FILTERING at all!)
    // Initialise the score for residuals 0
    // Store this as it serves as a baseline for the filtering analysis
    computeScores(data, tp, fp, numberOfMolecules, true);
    _residualsScoreMax = this.residualsScore;
    computeScores(data, tp, fp, numberOfMolecules, false);
    _residualsScoreAv = this.residualsScore;
    residualsScore = (useMaxResiduals) ? _residualsScoreMax : _residualsScoreAv;
    if (showJaccardPlot)
        plotJaccard(residualsScore, null);
    String bestJaccard = Utils.rounded(bestTp / (bestFp + numberOfMolecules)) + '\t';
    analysisPrefix += bestJaccard;
    sb.append(bestJaccard);
    addJaccardScores(sb);
    sb.append("\t").append(Utils.timeToString(runTime / 1000000.0));
    summaryTable.append(sb.toString());
}

Example 7 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class TraceDiffusion method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    jumpDistanceParameters = null;
    extraOptions = Utils.isExtraOptions();
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    ArrayList<MemoryPeakResults> allResults = new ArrayList<MemoryPeakResults>();
    // Option to pick multiple input datasets together using a list box.
    if ("multi".equals(arg)) {
        if (!showMultiDialog(allResults))
            return;
    }
    // This shows the dialog for selecting trace options
    if (!showTraceDialog(allResults))
        return;
    if (// Sense check
    allResults.isEmpty())
        return;
    Utils.log(TITLE + "...");
    // This optionally collects additional datasets then gets the traces:
    // - Trace each single dataset (and store in memory)
    // - Combine trace results held in memory
    Trace[] traces = getTraces(allResults);
    // This still allows a zero entry in the results table.
    if (traces.length > 0)
        if (!showDialog())
            return;
    int count = traces.length;
    double[] fitMSDResult = null;
    int n = 0;
    double[][] jdParams = null;
    if (count > 0) {
        calculatePrecision(traces, allResults.size() > 1);
        //--- MSD Analysis ---
        // Conversion constants
        final double px2ToUm2 = results.getCalibration().getNmPerPixel() * results.getCalibration().getNmPerPixel() / 1e6;
        final double px2ToUm2PerSecond = px2ToUm2 / exposureTime;
        // Get the maximum trace length
        int length = settings.minimumTraceLength;
        if (!settings.truncate) {
            for (Trace trace : traces) {
                if (length < trace.size())
                    length = trace.size();
            }
        }
        // Get the localisation error (4s^2) in um^2
        final double error = (settings.precisionCorrection) ? 4 * precision * precision / 1e6 : 0;
        // Pre-calculate MSD correction factors. This accounts for the fact that the distance moved 
        // in the start/end frames is reduced due to the averaging of the particle location over the 
        // entire frame into a single point. The true MSD may be restored by applying a factor.
        // Note: These are used for the calculation of the diffusion coefficients per molecule and 
        // the MSD passed to the Jump Distance analysis. However the error is not included in the 
        // jump distance analysis so will be subtracted from the fitted D coefficients later.
        final double[] factors;
        if (settings.msdCorrection) {
            factors = new double[length];
            for (int t = 1; t < length; t++) factors[t] = JumpDistanceAnalysis.getConversionfactor(t);
        } else {
            factors = Utils.newArray(length, 0.0, 1.0);
        }
        // Extract the mean-squared distance statistics
        Statistics[] stats = new Statistics[length];
        for (int i = 0; i < stats.length; i++) stats[i] = new Statistics();
        ArrayList<double[]> distances = (saveTraceDistances || displayTraceLength) ? new ArrayList<double[]>(traces.length) : null;
        // Store all the jump distances at the specified interval
        StoredDataStatistics jumpDistances = new StoredDataStatistics();
        final int jumpDistanceInterval = settings.jumpDistance;
        // Compute squared distances
        StoredDataStatistics msdPerMoleculeAllVsAll = new StoredDataStatistics();
        StoredDataStatistics msdPerMoleculeAdjacent = new StoredDataStatistics();
        for (Trace trace : traces) {
            ArrayList<PeakResult> results = trace.getPoints();
            // Sum the MSD and the time
            final int traceLength = (settings.truncate) ? settings.minimumTraceLength : trace.size();
            // Get the mean for each time separation
            double[] sumDistance = new double[traceLength + 1];
            double[] sumTime = new double[sumDistance.length];
            // Do the distances to the origin (saving if necessary)
            {
                final float x = results.get(0).getXPosition();
                final float y = results.get(0).getYPosition();
                if (distances != null) {
                    double[] msd = new double[traceLength - 1];
                    for (int j = 1; j < traceLength; j++) {
                        final int t = j;
                        final double d = distance2(x, y, results.get(j));
                        msd[j - 1] = px2ToUm2 * d;
                        if (t == jumpDistanceInterval)
                            jumpDistances.add(msd[j - 1]);
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                    distances.add(msd);
                } else {
                    for (int j = 1; j < traceLength; j++) {
                        final int t = j;
                        final double d = distance2(x, y, results.get(j));
                        if (t == jumpDistanceInterval)
                            jumpDistances.add(px2ToUm2 * d);
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                }
            }
            if (settings.internalDistances) {
                // Do the internal distances
                for (int i = 1; i < traceLength; i++) {
                    final float x = results.get(i).getXPosition();
                    final float y = results.get(i).getYPosition();
                    for (int j = i + 1; j < traceLength; j++) {
                        final int t = j - i;
                        final double d = distance2(x, y, results.get(j));
                        if (t == jumpDistanceInterval)
                            jumpDistances.add(px2ToUm2 * d);
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                }
                // Add the average distance per time separation to the population
                for (int t = 1; t < traceLength; t++) {
                    // Note: (traceLength - t) == count
                    stats[t].add(sumDistance[t] / (traceLength - t));
                }
            } else {
                // Add the distance per time separation to the population
                for (int t = 1; t < traceLength; t++) {
                    stats[t].add(sumDistance[t]);
                }
            }
            // Fix this for the precision and MSD adjustment.
            // It may be necessary to:
            // - sum the raw distances for each time interval (this is sumDistance[t])
            // - subtract the precision error
            // - apply correction factor for the n-frames to get actual MSD
            // - sum the actual MSD
            double sumD = 0, sumD_adjacent = Math.max(0, sumDistance[1] - error) * factors[1];
            double sumT = 0, sumT_adjacent = sumTime[1];
            for (int t = 1; t < traceLength; t++) {
                sumD += Math.max(0, sumDistance[t] - error) * factors[t];
                sumT += sumTime[t];
            }
            // Calculate the average displacement for the trace (do not simply use the largest 
            // time separation since this will miss moving molecules that end up at the origin)
            msdPerMoleculeAllVsAll.add(px2ToUm2PerSecond * sumD / sumT);
            msdPerMoleculeAdjacent.add(px2ToUm2PerSecond * sumD_adjacent / sumT_adjacent);
        }
        StoredDataStatistics dPerMoleculeAllVsAll = null;
        StoredDataStatistics dPerMoleculeAdjacent = null;
        if (saveTraceDistances || (settings.showHistograms && displayDHistogram)) {
            dPerMoleculeAllVsAll = calculateDiffusionCoefficient(msdPerMoleculeAllVsAll);
            dPerMoleculeAdjacent = calculateDiffusionCoefficient(msdPerMoleculeAdjacent);
        }
        if (saveTraceDistances) {
            saveTraceDistances(traces.length, distances, msdPerMoleculeAllVsAll, msdPerMoleculeAdjacent, dPerMoleculeAllVsAll, dPerMoleculeAdjacent);
        }
        if (displayTraceLength) {
            StoredDataStatistics lengths = calculateTraceLengths(distances);
            showHistogram(lengths, "Trace length (um)");
        }
        if (displayTraceSize) {
            StoredDataStatistics sizes = calculateTraceSizes(traces);
            showHistogram(sizes, "Trace size", true);
        }
        // Plot the per-trace histogram of MSD and D
        if (settings.showHistograms) {
            if (displayMSDHistogram) {
                showHistogram(msdPerMoleculeAllVsAll, "MSD/Molecule (all-vs-all)");
                showHistogram(msdPerMoleculeAdjacent, "MSD/Molecule (adjacent)");
            }
            if (displayDHistogram) {
                showHistogram(dPerMoleculeAllVsAll, "D/Molecule (all-vs-all)");
                showHistogram(dPerMoleculeAdjacent, "D/Molecule (adjacent)");
            }
        }
        // Calculate the mean squared distance (MSD)
        double[] x = new double[stats.length];
        double[] y = new double[x.length];
        double[] sd = new double[x.length];
        // Intercept is the 4s^2 (in um^2)
        y[0] = 4 * precision * precision / 1e6;
        for (int i = 1; i < stats.length; i++) {
            x[i] = i * exposureTime;
            y[i] = stats[i].getMean() * px2ToUm2;
            //sd[i] = stats[i].getStandardDeviation() * px2ToUm2;
            sd[i] = stats[i].getStandardError() * px2ToUm2;
        }
        String title = TITLE + " MSD";
        Plot2 plot = plotMSD(x, y, sd, title);
        // Fit the MSD using a linear fit
        fitMSDResult = fitMSD(x, y, title, plot);
        // Jump Distance analysis
        if (saveRawData)
            saveStatistics(jumpDistances, "Jump Distance", "Distance (um^2)", false);
        // Calculate the cumulative jump-distance histogram
        double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(jumpDistances.getValues());
        // Always show the jump distance histogram
        jdTitle = TITLE + " Jump Distance";
        jdPlot = new Plot2(jdTitle, "Distance (um^2)", "Cumulative Probability", jdHistogram[0], jdHistogram[1]);
        display(jdTitle, jdPlot);
        // Fit Jump Distance cumulative probability
        n = jumpDistances.getN();
        jumpDistanceParameters = jdParams = fitJumpDistance(jumpDistances, jdHistogram);
    }
    summarise(traces, fitMSDResult, n, jdParams);
}

Example 8 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class TraceDiffusion method calculateTraceLengths.

public StoredDataStatistics calculateTraceLengths(ArrayList<double[]> distances) {
    StoredDataStatistics lengths = new StoredDataStatistics();
    for (double[] trace : distances) {
        double sum = 0;
        for (double d : trace) {
            sum += Math.sqrt(d);
        }
        lengths.add(sum);
    }
    return lengths;
}

Example 9 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class SpotInspector method plotHistogram.

private void plotHistogram(float[] data, double yMin, double yMax) {
    if (plotHistogram) {
        String title = TITLE + " Histogram";
        Utils.showHistogram(title, new StoredDataStatistics(data), SORT_ORDER[sortOrderIndex], 0, (removeOutliers) ? 1 : 0, histogramBins);
    }
}

Example 10 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class CreateData method createPhotonDistribution.

/**
	 * @return A photon distribution loaded from a file of floating-point values with the specified population mean.
	 */
private RealDistribution createPhotonDistribution() {
    if (PHOTON_DISTRIBUTION[PHOTON_CUSTOM].equals(settings.photonDistribution)) {
        // Get the distribution file
        String filename = Utils.getFilename("Photon_distribution", settings.photonDistributionFile);
        if (filename != null) {
            settings.photonDistributionFile = filename;
            try {
                InputStream is = new FileInputStream(new File(settings.photonDistributionFile));
                BufferedReader in = new BufferedReader(new UnicodeReader(is, null));
                StoredDataStatistics stats = new StoredDataStatistics();
                try {
                    String str = null;
                    double val = 0.0d;
                    while ((str = in.readLine()) != null) {
                        val = Double.parseDouble(str);
                        stats.add(val);
                    }
                } finally {
                    in.close();
                }
                if (stats.getSum() > 0) {
                    // Update the statistics to the desired mean.
                    double scale = (double) settings.photonsPerSecond / stats.getMean();
                    double[] values = stats.getValues();
                    for (int i = 0; i < values.length; i++) values[i] *= scale;
                    // TODO - Investigate the limits of this distribution. 
                    // How far above and below the input data will values be generated.
                    // Create the distribution using the recommended number of bins
                    final int binCount = stats.getN() / 10;
                    EmpiricalDistribution dist = new EmpiricalDistribution(binCount, createRandomGenerator());
                    dist.load(values);
                    return dist;
                }
            } catch (IOException e) {
            // Ignore
            } catch (NullArgumentException e) {
            // Ignore 
            } catch (NumberFormatException e) {
            // Ignore
            }
        }
        Utils.log("Failed to load custom photon distribution from file: %s. Default to fixed.", settings.photonDistributionFile);
    } else if (PHOTON_DISTRIBUTION[PHOTON_UNIFORM].equals(settings.photonDistribution)) {
        if (settings.photonsPerSecond < settings.photonsPerSecondMaximum) {
            UniformRealDistribution dist = new UniformRealDistribution(createRandomGenerator(), settings.photonsPerSecond, settings.photonsPerSecondMaximum);
            return dist;
        }
    } else if (PHOTON_DISTRIBUTION[PHOTON_GAMMA].equals(settings.photonDistribution)) {
        final double scaleParameter = settings.photonsPerSecond / settings.photonShape;
        GammaDistribution dist = new GammaDistribution(createRandomGenerator(), settings.photonShape, scaleParameter, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
        return dist;
    } else if (PHOTON_DISTRIBUTION[PHOTON_CORRELATED].equals(settings.photonDistribution)) {
        // No distribution required
        return null;
    }
    settings.photonDistribution = PHOTON_DISTRIBUTION[PHOTON_FIXED];
    return null;
}