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

Example 16 with StoredDataStatistics

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

the class ImageModel method createImage.

/**
 * Simulate an image of fluorophores. The total amount of time a fluorophore is on (i.e. sum of
 * tOn) is used to create the signal strength using the specified correlation.
 *
 * <p>A second set of fluorophores, independent of the first, are generated using
 * {@link #createFluorophores(List, int)}. The correlated on-times will be created by combining
 * the times using the provided correlation (r):
 *
 * <pre>
 * // X1 : Fluorophore total on-times
 * // X2 : Independently generated fluorophore total on-times
 * // r  : correlation
 * a = sqrt(1 - r * r)
 * newX = (r * X1 + a * X2) / (r + a)
 * </pre>
 *
 * <p>The signal is proportional to newly generated on-times (newX) with an average of the
 * specified photon budget.
 *
 * <p>The photon budget can either be distributed evenly over the fluorophore lifetime or per
 * frame (see {@link #isPhotonBudgetPerFrame()}). Each frame signal output will be subject to
 * Poisson sampling.
 *
 * <p>If the input correlation is zero then the number of photons will be sampled from the
 * configured photon distribution or, if this is null, will be uniform.
 *
 * <p>A random fraction of the fluorophores will move using a random walk with the diffusion
 * coefficient defined in the compound.
 *
 * @param compoundFluorophores The compounds containing the fluorophores
 * @param fixedFraction The fraction of molecules that are fixed
 * @param maxFrames The maximum frame for the simulation
 * @param photonBudget the average number of photons per frame/lifetime; see
 *        {@link #isPhotonBudgetPerFrame()}
 * @param correlation The correlation between the number of photons and the total on-time of the
 *        fluorophore
 * @param rotate Rotate the molecule per frame
 * @return the localisations
 */
public List<LocalisationModel> createImage(List<CompoundMoleculeModel> compoundFluorophores, double fixedFraction, int maxFrames, double photonBudget, double correlation, boolean rotate) {
    final List<LocalisationModel> localisations = new ArrayList<>();
    // Extract the fluorophores in all the compounds
    final ArrayList<FluorophoreSequenceModel> fluorophores = new ArrayList<>(compoundFluorophores.size());
    for (final CompoundMoleculeModel c : compoundFluorophores) {
        for (int i = c.getSize(); i-- > 0; ) {
            if (c.getMolecule(i) instanceof FluorophoreSequenceModel) {
                fluorophores.add((FluorophoreSequenceModel) c.getMolecule(i));
            }
        }
    }
    final int nMolecules = fluorophores.size();
    // Check the correlation bounds.
    // Correlation for photons per frame verses total on time should be negative.
    // Correlation for total photons verses total on time should be positive.
    double boundedCorrelation;
    if (photonBudgetPerFrame) {
        boundedCorrelation = MathUtils.clip(-1, 0, correlation);
    } else {
        boundedCorrelation = MathUtils.clip(0, 1, correlation);
    }
    // Create a photon budget for each fluorophore
    final double[] photons = new double[nMolecules];
    // Generate a second set of on times using the desired correlation
    if (boundedCorrelation != 0) {
        // Observations on real data show:
        // - there is a weak positive correlation between total photons and time
        // - There is a weak negative correlation between photons per frame and total on-time
        // How to generate a random correlation:
        // http://www.uvm.edu/~dhowell/StatPages/More_Stuff/CorrGen.html
        // http://stats.stackexchange.com/questions/13382/how-to-define-a-distribution-such-that-draws-from-it-correlate-with-a-draw-from
        // Create a second set of fluorophores. This is used to generate the correlated photon data
        final List<FluorophoreSequenceModel> fluorophores2 = new ArrayList<>();
        while (fluorophores2.size() < fluorophores.size()) {
            final FluorophoreSequenceModel f = createFluorophore(0, new double[] { 0, 0, 0 }, maxFrames);
            if (f != null) {
                fluorophores2.add(f);
            }
        }
        final double a = Math.sqrt(1 - boundedCorrelation * boundedCorrelation);
        // Q - How to generate a negative correlation?
        // Generate a positive correlation then invert the data and shift to the same range
        final boolean invert = (boundedCorrelation < 0);
        boundedCorrelation = Math.abs(boundedCorrelation);
        StoredDataStatistics correlatedOnTime = new StoredDataStatistics();
        for (int i = 0; i < nMolecules; i++) {
            final double X = getTotalOnTime(fluorophores.get(i));
            final double Z = getTotalOnTime(fluorophores2.get(i));
            correlatedOnTime.add((boundedCorrelation * X + a * Z) / (boundedCorrelation + a));
        }
        if (invert) {
            final double min = correlatedOnTime.getStatistics().getMin();
            final double range = correlatedOnTime.getStatistics().getMax() - min;
            final StoredDataStatistics newCorrelatedOnTime = new StoredDataStatistics();
            for (final double d : correlatedOnTime.getValues()) {
                newCorrelatedOnTime.add(range - d + min);
            }
            correlatedOnTime = newCorrelatedOnTime;
        }
        // Get the average on time from the correlated sample.
        // Using the population value (tOn * (1+blinks1)) over-estimates the on time.
        final double averageTotalTOn = correlatedOnTime.getMean();
        // Now create the localisations
        final double[] correlatedOnTimes = correlatedOnTime.getValues();
        for (int i = 0; i < nMolecules; i++) {
            // Generate photons using the correlated on time data
            final double p = photonBudget * correlatedOnTimes[i] / averageTotalTOn;
            photons[i] = p;
        }
    } else if (photonDistribution != null) {
        // Sample from the provided distribution. Do not over-write the class level distribution
        // to allow running again with a different shape parameter / photon budget.
        // Ensure the custom distribution is scaled to the correct photon budget
        final double photonScale = photonBudget / photonDistribution.getNumericalMean();
        // Generate photons sampling from the photon budget
        for (int i = 0; i < nMolecules; i++) {
            photons[i] = photonDistribution.sample() * photonScale;
        }
    } else {
        // No distribution so use the same number for all
        Arrays.fill(photons, photonBudget);
    }
    int photonIndex = 0;
    if (fixedFraction > 0) {
        for (final CompoundMoleculeModel c : compoundFluorophores) {
            final boolean fixed = random.nextDouble() < fixedFraction;
            photonIndex += createLocalisation(c, localisations, fixed, maxFrames, photons, photonIndex, !fixed && rotate);
        }
    } else {
        // No fixed molecules
        for (final CompoundMoleculeModel c : compoundFluorophores) {
            photonIndex += createLocalisation(c, localisations, false, maxFrames, photons, photonIndex, rotate);
        }
    }
    sortByTime(localisations);
    return localisations;
}

Example 17 with StoredDataStatistics

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

the class PoissonGammaGaussianFunctionTest method fasterThan.

private void fasterThan(RandomSeed seed, ConvolutionMode slow, ConvolutionMode fast) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    // Realistic EM-CCD parameters for speed test
    final double s = 7.16;
    final double g = 39.1;
    final PoissonGammaGaussianFunction f1 = new PoissonGammaGaussianFunction(1 / g, s);
    f1.setConvolutionMode(slow);
    final PoissonGammaGaussianFunction f2 = new PoissonGammaGaussianFunction(1 / g, s);
    f2.setConvolutionMode(fast);
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    // Generate realistic data from the probability mass function
    final double[][] samples = new double[photons.length][];
    for (int j = 0; j < photons.length; j++) {
        final int start = (int) (4 * -s);
        int mu = start;
        final StoredDataStatistics stats = new StoredDataStatistics();
        while (stats.getSum() < 0.995) {
            final double p = f1.likelihood(mu, photons[j]);
            stats.add(p);
            if (mu > 10 && p / stats.getSum() < 1e-6) {
                break;
            }
            mu++;
        }
        // Generate cumulative probability
        final double[] data = stats.getValues();
        for (int i = 1; i < data.length; i++) {
            data[i] += data[i - 1];
        }
        // Normalise
        for (int i = 0, end = data.length - 1; i < data.length; i++) {
            data[i] /= data[end];
        }
        // Sample
        final double[] sample = new double[1000];
        for (int i = 0; i < sample.length; i++) {
            final double p = rg.nextDouble();
            int x = 0;
            while (x < data.length && data[x] < p) {
                x++;
            }
            sample[i] = start + x;
        }
        samples[j] = sample;
    }
    // Warm-up
    run(f1, samples, photons);
    run(f2, samples, photons);
    long t1 = 0;
    for (int i = 0; i < 5; i++) {
        t1 += run(f1, samples, photons);
    }
    long t2 = 0;
    for (int i = 0; i < 5; i++) {
        t2 += run(f2, samples, photons);
    }
    logger.log(TestLogUtils.getTimingRecord(getName(f1), t1, getName(f2), t2));
}

Example 18 with StoredDataStatistics

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

the class PoissonGaussianConvolutionFunctionTest method pdfFasterThanPmf.

@SpeedTag
@SeededTest
void pdfFasterThanPmf(RandomSeed seed) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    // Realistic CCD parameters for speed test
    final double s = 7.16;
    final double g = 3.1;
    final PoissonGaussianConvolutionFunction f1 = PoissonGaussianConvolutionFunction.createWithStandardDeviation(1 / g, s);
    f1.setComputePmf(true);
    final PoissonGaussianConvolutionFunction f2 = PoissonGaussianConvolutionFunction.createWithStandardDeviation(1 / g, s);
    f2.setComputePmf(false);
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    // Generate realistic data from the probability mass function
    final double[][] samples = new double[photons.length][];
    for (int j = 0; j < photons.length; j++) {
        final int start = (int) (4 * -s);
        int mu = start;
        final StoredDataStatistics stats = new StoredDataStatistics();
        while (stats.getSum() < 0.995) {
            final double p = f1.likelihood(mu, photons[j]);
            stats.add(p);
            if (mu > 10 && p / stats.getSum() < 1e-6) {
                break;
            }
            mu++;
        }
        // Generate cumulative probability
        final double[] data = stats.getValues();
        for (int i = 1; i < data.length; i++) {
            data[i] += data[i - 1];
        }
        // Normalise
        for (int i = 0, end = data.length - 1; i < data.length; i++) {
            data[i] /= data[end];
        }
        // Sample
        final double[] sample = new double[1000];
        for (int i = 0; i < sample.length; i++) {
            final double p = rg.nextDouble();
            int x = 0;
            while (x < data.length && data[x] < p) {
                x++;
            }
            sample[i] = start + x;
        }
        samples[j] = sample;
    }
    // Warm-up
    run(f1, samples, photons);
    run(f2, samples, photons);
    long t1 = 0;
    for (int i = 0; i < 5; i++) {
        t1 += run(f1, samples, photons);
    }
    long t2 = 0;
    for (int i = 0; i < 5; i++) {
        t2 += run(f2, samples, photons);
    }
    logger.log(TestLogUtils.getTimingRecord("cdf", t1, "pdf", t2));
}

Example 19 with StoredDataStatistics

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

the class DiffusionRateTest method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    pluginSettings = Settings.load();
    pluginSettings.save();
    if (IJ.controlKeyDown()) {
        simpleTest();
        return;
    }
    extraOptions = ImageJUtils.isExtraOptions();
    if (!showDialog()) {
        return;
    }
    lastSimulation.set(null);
    final int totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
    conversionFactor = 1000000.0 / (settings.getPixelPitch() * settings.getPixelPitch());
    // Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
    final double diffusionRateInPixelsPerSecond = settings.getDiffusionRate() * conversionFactor;
    final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.getStepsPerSecond();
    final double precisionInPixels = myPrecision / settings.getPixelPitch();
    final boolean addError = myPrecision != 0;
    ImageJUtils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", MathUtils.rounded(settings.getDiffusionRate(), 4), MathUtils.rounded(myPrecision, 4));
    ImageJUtils.log("Mean-displacement per dimension = %s nm/sec", MathUtils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.getDiffusionRate()), 4));
    if (extraOptions) {
        ImageJUtils.log("Step size = %s, precision = %s", MathUtils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), MathUtils.rounded(precisionInPixels));
    }
    // Convert diffusion co-efficient into the standard deviation for the random walk
    final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
    final double diffusionSigma = ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
    ImageJUtils.log("Simulation step-size = %s nm", MathUtils.rounded(settings.getPixelPitch() * diffusionSigma, 4));
    // Move the molecules and get the diffusion rate
    IJ.showStatus("Simulating ...");
    final long start = System.nanoTime();
    final UniformRandomProvider random = UniformRandomProviders.create();
    final Statistics[] stats2D = new Statistics[totalSteps];
    final Statistics[] stats3D = new Statistics[totalSteps];
    final StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
    final StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
    for (int j = 0; j < totalSteps; j++) {
        stats2D[j] = new Statistics();
        stats3D[j] = new Statistics();
    }
    final SphericalDistribution dist = new SphericalDistribution(settings.getConfinementRadius() / settings.getPixelPitch());
    final Statistics asymptote = new Statistics();
    // Save results to memory
    final MemoryPeakResults results = new MemoryPeakResults(totalSteps);
    results.setCalibration(CalibrationHelper.create(settings.getPixelPitch(), 1, 1000.0 / settings.getStepsPerSecond()));
    results.setName(TITLE);
    results.setPsf(PsfHelper.create(PSFType.CUSTOM));
    int peak = 0;
    // Store raw coordinates
    final ArrayList<Point> points = new ArrayList<>(totalSteps);
    final StoredData totalJumpDistances1D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances2D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances3D = new StoredData(settings.getParticles());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
    for (int i = 0; i < settings.getParticles(); i++) {
        if (i % 16 == 0) {
            IJ.showProgress(i, settings.getParticles());
            if (ImageJUtils.isInterrupted()) {
                return;
            }
        }
        // Increment the frame so that tracing analysis can distinguish traces
        peak++;
        double[] origin = new double[3];
        final int id = i + 1;
        final MoleculeModel m = new MoleculeModel(id, origin.clone());
        if (addError) {
            origin = addError(origin, precisionInPixels, gauss);
        }
        if (pluginSettings.useConfinement) {
            // Note: When using confinement the average displacement should asymptote
            // at the average distance of a point from the centre of a ball. This is 3r/4.
            // See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM
            // The equivalent in 2D is 2r/3. However although we are plotting 2D distance
            // this is a projection of the 3D position onto the plane and so the particles
            // will not be evenly spread (there will be clustering at centre caused by the
            // poles)
            final double[] axis = (diffusionType == DiffusionType.LINEAR_WALK) ? nextVector(gauss) : null;
            for (int j = 0; j < totalSteps; j++) {
                double[] xyz = m.getCoordinates();
                final double[] originalXyz = xyz.clone();
                for (int n = pluginSettings.confinementAttempts; n-- > 0; ) {
                    if (diffusionType == DiffusionType.GRID_WALK) {
                        m.walk(diffusionSigma, random);
                    } else if (diffusionType == DiffusionType.LINEAR_WALK) {
                        m.slide(diffusionSigma, axis, random);
                    } else {
                        m.move(diffusionSigma, random);
                    }
                    if (!dist.isWithin(m.getCoordinates())) {
                        // Reset position
                        for (int k = 0; k < 3; k++) {
                            xyz[k] = originalXyz[k];
                        }
                    } else {
                        // The move was allowed
                        break;
                    }
                }
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
            asymptote.add(distance(m.getCoordinates()));
        } else if (diffusionType == DiffusionType.GRID_WALK) {
            for (int j = 0; j < totalSteps; j++) {
                m.walk(diffusionSigma, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else if (diffusionType == DiffusionType.LINEAR_WALK) {
            final double[] axis = nextVector(gauss);
            for (int j = 0; j < totalSteps; j++) {
                m.slide(diffusionSigma, axis, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else {
            for (int j = 0; j < totalSteps; j++) {
                m.move(diffusionSigma, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        }
        // Debug: record all the particles so they can be analysed
        // System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ());
        final double[] xyz = m.getCoordinates();
        double d2 = 0;
        totalJumpDistances1D.add(d2 = xyz[0] * xyz[0]);
        totalJumpDistances2D.add(d2 += xyz[1] * xyz[1]);
        totalJumpDistances3D.add(d2 += xyz[2] * xyz[2]);
    }
    final long nanoseconds = System.nanoTime() - start;
    IJ.showProgress(1);
    MemoryPeakResults.addResults(results);
    simulation = new SimulationData(results.getName(), myPrecision);
    // Convert pixels^2/step to um^2/sec
    final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    ImageJUtils.log("Raw data D=%s um^2/s, Precision = %s nm, N=%d, step=%s s, mean2D=%s um^2, " + "MSD 2D = %s um^2/s, mean3D=%s um^2, MSD 3D = %s um^2/s", MathUtils.rounded(settings.getDiffusionRate()), MathUtils.rounded(myPrecision), jumpDistances2D.getN(), MathUtils.rounded(results.getCalibrationReader().getExposureTime() / 1000), MathUtils.rounded(jumpDistances2D.getMean() / conversionFactor), MathUtils.rounded(msd2D), MathUtils.rounded(jumpDistances3D.getMean() / conversionFactor), MathUtils.rounded(msd3D));
    aggregateIntoFrames(points, addError, precisionInPixels, gauss);
    IJ.showStatus("Analysing results ...");
    if (pluginSettings.showDiffusionExample) {
        showExample(totalSteps, diffusionSigma, random);
    }
    // Plot a graph of mean squared distance
    final double[] xValues = new double[stats2D.length];
    final double[] yValues2D = new double[stats2D.length];
    final double[] yValues3D = new double[stats3D.length];
    final double[] upper2D = new double[stats2D.length];
    final double[] lower2D = new double[stats2D.length];
    final double[] upper3D = new double[stats3D.length];
    final double[] lower3D = new double[stats3D.length];
    final SimpleRegression r2D = new SimpleRegression(false);
    final SimpleRegression r3D = new SimpleRegression(false);
    final int firstN = (pluginSettings.useConfinement) ? pluginSettings.fitN : totalSteps;
    for (int j = 0; j < totalSteps; j++) {
        // Convert steps to seconds
        xValues[j] = (j + 1) / settings.getStepsPerSecond();
        // Convert values in pixels^2 to um^2
        final double mean2D = stats2D[j].getMean() / conversionFactor;
        final double mean3D = stats3D[j].getMean() / conversionFactor;
        final double sd2D = stats2D[j].getStandardDeviation() / conversionFactor;
        final double sd3D = stats3D[j].getStandardDeviation() / conversionFactor;
        yValues2D[j] = mean2D;
        yValues3D[j] = mean3D;
        upper2D[j] = mean2D + sd2D;
        lower2D[j] = mean2D - sd2D;
        upper3D[j] = mean3D + sd3D;
        lower3D[j] = mean3D - sd3D;
        if (j < firstN) {
            r2D.addData(xValues[j], yValues2D[j]);
            r3D.addData(xValues[j], yValues3D[j]);
        }
    }
    // TODO - Fit using the equation for 2D confined diffusion:
    // MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2)
    // s = localisation precision
    // R = confinement radius
    // D = 2D diffusion coefficient
    // t = time
    final PolynomialFunction fitted2D;
    final PolynomialFunction fitted3D;
    if (r2D.getN() > 0) {
        // Do linear regression to get diffusion rate
        final double[] best2D = new double[] { r2D.getIntercept(), r2D.getSlope() };
        fitted2D = new PolynomialFunction(best2D);
        final double[] best3D = new double[] { r3D.getIntercept(), r3D.getSlope() };
        fitted3D = new PolynomialFunction(best3D);
        // For 2D diffusion: d^2 = 4D
        // where: d^2 = mean-square displacement
        double diffCoeff = best2D[1] / 4.0;
        final String msg = "2D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")";
        IJ.showStatus(msg);
        ImageJUtils.log(msg);
        diffCoeff = best3D[1] / 6.0;
        ImageJUtils.log("3D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")");
    } else {
        fitted2D = fitted3D = null;
    }
    // Create plots
    plotMsd(totalSteps, xValues, yValues2D, lower2D, upper2D, fitted2D, 2);
    plotMsd(totalSteps, xValues, yValues3D, lower3D, upper3D, fitted3D, 3);
    plotJumpDistances(TITLE, jumpDistances2D, 2, 1);
    plotJumpDistances(TITLE, jumpDistances3D, 3, 1);
    // Show the total jump length for debugging
    // plotJumpDistances(TITLE + " total", totalJumpDistances1D, 1, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances2D, 2, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances3D, 3, totalSteps);
    windowOrganiser.tile();
    if (pluginSettings.useConfinement) {
        ImageJUtils.log("3D asymptote distance = %s nm (expected %.2f)", MathUtils.rounded(asymptote.getMean() * settings.getPixelPitch(), 4), 3 * settings.getConfinementRadius() / 4);
    }
}

Example 20 with StoredDataStatistics

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

the class BenchmarkSpotFit method summariseResults.

private void summariseResults(BenchmarkSpotFitResult spotFitResults, long runTime, final PreprocessedPeakResult[] preprocessedPeakResults, int uniqueIdCount, CandidateData candidateData, TIntObjectHashMap<List<Coordinate>> actualCoordinates) {
    // Summarise the fitting results. N fits, N failures.
    // Optimal match statistics if filtering is perfect (since fitting is not perfect).
    final StoredDataStatistics distanceStats = new StoredDataStatistics();
    final StoredDataStatistics depthStats = new StoredDataStatistics();
    // Get stats for all fitted results and those that match
    // Signal, SNR, Width, xShift, yShift, Precision
    createFilterCriteria();
    final StoredDataStatistics[][] stats = new StoredDataStatistics[3][filterCriteria.length];
    for (int i = 0; i < stats.length; i++) {
        for (int j = 0; j < stats[i].length; j++) {
            stats[i][j] = new StoredDataStatistics();
        }
    }
    final double nmPerPixel = simulationParameters.pixelPitch;
    double tp = 0;
    double fp = 0;
    int failCtp = 0;
    int failCfp = 0;
    int ctp = 0;
    int cfp = 0;
    final int[] singleStatus = new int[FitStatus.values().length];
    final int[] multiStatus = new int[singleStatus.length];
    final int[] doubletStatus = new int[singleStatus.length];
    final int[] multiDoubletStatus = new int[singleStatus.length];
    // Easier to materialise the values since we have a lot of non final variables to manipulate
    final TIntObjectHashMap<FilterCandidates> fitResults = spotFitResults.fitResults;
    final int[] frames = new int[fitResults.size()];
    final FilterCandidates[] candidates = new FilterCandidates[fitResults.size()];
    final int[] counter = new int[1];
    fitResults.forEachEntry((frame, candidate) -> {
        frames[counter[0]] = frame;
        candidates[counter[0]] = candidate;
        counter[0]++;
        return true;
    });
    for (final FilterCandidates result : candidates) {
        // Count the number of fit results that matched (tp) and did not match (fp)
        tp += result.tp;
        fp += result.fp;
        for (int i = 0; i < result.fitResult.length; i++) {
            if (result.spots[i].match) {
                ctp++;
            } else {
                cfp++;
            }
            final MultiPathFitResult fitResult = result.fitResult[i];
            if (singleStatus != null && result.spots[i].match) {
                // Debugging reasons for fit failure
                addStatus(singleStatus, fitResult.getSingleFitResult());
                addStatus(multiStatus, fitResult.getMultiFitResult());
                addStatus(doubletStatus, fitResult.getDoubletFitResult());
                addStatus(multiDoubletStatus, fitResult.getMultiDoubletFitResult());
            }
            if (noMatch(fitResult)) {
                if (result.spots[i].match) {
                    failCtp++;
                } else {
                    failCfp++;
                }
            }
            // We have multi-path results.
            // We want statistics for:
            // [0] all fitted spots
            // [1] fitted spots that match a result
            // [2] fitted spots that do not match a result
            addToStats(fitResult.getSingleFitResult(), stats);
            addToStats(fitResult.getMultiFitResult(), stats);
            addToStats(fitResult.getDoubletFitResult(), stats);
            addToStats(fitResult.getMultiDoubletFitResult(), stats);
        }
        // Statistics on spots that fit an actual result
        for (int i = 0; i < result.match.length; i++) {
            if (!result.match[i].isFitResult()) {
                // For now just ignore the candidates that matched
                continue;
            }
            final FitMatch fitMatch = (FitMatch) result.match[i];
            distanceStats.add(fitMatch.distance * nmPerPixel);
            depthStats.add(fitMatch.zdepth * nmPerPixel);
        }
    }
    if (tp == 0) {
        IJ.error(TITLE, "No fit results matched the simulation actual results");
        return;
    }
    // Store data for computing correlation
    final double[] i1 = new double[depthStats.getN()];
    final double[] i2 = new double[i1.length];
    final double[] is = new double[i1.length];
    int ci = 0;
    for (final FilterCandidates result : candidates) {
        for (int i = 0; i < result.match.length; i++) {
            if (!result.match[i].isFitResult()) {
                // For now just ignore the candidates that matched
                continue;
            }
            final FitMatch fitMatch = (FitMatch) result.match[i];
            final ScoredSpot spot = result.spots[fitMatch.index];
            i1[ci] = fitMatch.predictedSignal;
            i2[ci] = fitMatch.actualSignal;
            is[ci] = spot.spot.intensity;
            ci++;
        }
    }
    // We want to compute the Jaccard against the spot metric
    // Filter the results using the multi-path filter
    final ArrayList<MultiPathFitResults> multiPathResults = new ArrayList<>(fitResults.size());
    for (int i = 0; i < frames.length; i++) {
        final int frame = frames[i];
        final MultiPathFitResult[] multiPathFitResults = candidates[i].fitResult;
        final int totalCandidates = candidates[i].spots.length;
        final List<Coordinate> list = actualCoordinates.get(frame);
        final int nActual = (list == null) ? 0 : list.size();
        multiPathResults.add(new MultiPathFitResults(frame, multiPathFitResults, totalCandidates, nActual));
    }
    // Score the results and count the number returned
    final List<FractionalAssignment[]> assignments = new ArrayList<>();
    final TIntHashSet set = new TIntHashSet(uniqueIdCount);
    final FractionScoreStore scoreStore = set::add;
    final MultiPathFitResults[] multiResults = multiPathResults.toArray(new MultiPathFitResults[0]);
    // Filter with no filter
    final MultiPathFilter mpf = new MultiPathFilter(new SignalFilter(0), null, multiFilter.residualsThreshold);
    mpf.fractionScoreSubset(multiResults, NullFailCounter.INSTANCE, this.results.size(), assignments, scoreStore, CoordinateStoreFactory.create(0, 0, imp.getWidth(), imp.getHeight(), config.convertUsingHwhMax(config.getDuplicateDistanceParameter())));
    final double[][] matchScores = new double[set.size()][];
    int count = 0;
    for (int i = 0; i < assignments.size(); i++) {
        final FractionalAssignment[] a = assignments.get(i);
        if (a == null) {
            continue;
        }
        for (int j = 0; j < a.length; j++) {
            final PreprocessedPeakResult r = ((PeakFractionalAssignment) a[j]).peakResult;
            set.remove(r.getUniqueId());
            final double precision = Math.sqrt(r.getLocationVariance());
            final double signal = r.getSignal();
            final double snr = r.getSnr();
            final double width = r.getXSdFactor();
            final double xShift = r.getXRelativeShift2();
            final double yShift = r.getYRelativeShift2();
            // Since these two are combined for filtering and the max is what matters.
            final double shift = (xShift > yShift) ? Math.sqrt(xShift) : Math.sqrt(yShift);
            final double eshift = Math.sqrt(xShift + yShift);
            final double[] score = new double[8];
            score[FILTER_SIGNAL] = signal;
            score[FILTER_SNR] = snr;
            score[FILTER_MIN_WIDTH] = width;
            score[FILTER_MAX_WIDTH] = width;
            score[FILTER_SHIFT] = shift;
            score[FILTER_ESHIFT] = eshift;
            score[FILTER_PRECISION] = precision;
            score[FILTER_PRECISION + 1] = a[j].getScore();
            matchScores[count++] = score;
        }
    }
    // Add the rest
    set.forEach(new CustomTIntProcedure(count) {

        @Override
        public boolean execute(int uniqueId) {
            // This should not be null or something has gone wrong
            final PreprocessedPeakResult r = preprocessedPeakResults[uniqueId];
            if (r == null) {
                throw new IllegalArgumentException("Missing result: " + uniqueId);
            }
            final double precision = Math.sqrt(r.getLocationVariance());
            final double signal = r.getSignal();
            final double snr = r.getSnr();
            final double width = r.getXSdFactor();
            final double xShift = r.getXRelativeShift2();
            final double yShift = r.getYRelativeShift2();
            // Since these two are combined for filtering and the max is what matters.
            final double shift = (xShift > yShift) ? Math.sqrt(xShift) : Math.sqrt(yShift);
            final double eshift = Math.sqrt(xShift + yShift);
            final double[] score = new double[8];
            score[FILTER_SIGNAL] = signal;
            score[FILTER_SNR] = snr;
            score[FILTER_MIN_WIDTH] = width;
            score[FILTER_MAX_WIDTH] = width;
            score[FILTER_SHIFT] = shift;
            score[FILTER_ESHIFT] = eshift;
            score[FILTER_PRECISION] = precision;
            matchScores[count++] = score;
            return true;
        }
    });
    final FitConfiguration fitConfig = config.getFitConfiguration();
    // Debug the reasons the fit failed
    if (singleStatus != null) {
        String name = PeakFit.getSolverName(fitConfig);
        if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
            name += " Camera";
        }
        IJ.log("Failure counts: " + name);
        printFailures("Single", singleStatus);
        printFailures("Multi", multiStatus);
        printFailures("Doublet", doubletStatus);
        printFailures("Multi doublet", multiDoubletStatus);
    }
    final StringBuilder sb = new StringBuilder(300);
    // Add information about the simulation
    final double signal = simulationParameters.averageSignal;
    final int n = results.size();
    sb.append(imp.getStackSize()).append('\t');
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    sb.append(w).append('\t');
    sb.append(h).append('\t');
    sb.append(n).append('\t');
    final double density = ((double) n / imp.getStackSize()) / (w * h) / (simulationParameters.pixelPitch * simulationParameters.pixelPitch / 1e6);
    sb.append(MathUtils.rounded(density)).append('\t');
    sb.append(MathUtils.rounded(signal)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.pixelPitch)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.depth)).append('\t');
    sb.append(simulationParameters.fixedDepth).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.gain)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.readNoise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.background)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.noise)).append('\t');
    if (simulationParameters.fullSimulation) {
    // The total signal is spread over frames
    }
    sb.append(MathUtils.rounded(signal / simulationParameters.noise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd / simulationParameters.pixelPitch)).append('\t');
    sb.append(spotFilter.getDescription());
    // nP and nN is the fractional score of the spot candidates
    addCount(sb, (double) candidateData.countPositive + candidateData.countNegative);
    addCount(sb, candidateData.countPositive);
    addCount(sb, candidateData.countNegative);
    addCount(sb, candidateData.fractionPositive);
    addCount(sb, candidateData.fractionNegative);
    String name = PeakFit.getSolverName(fitConfig);
    if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
        name += " Camera";
    }
    add(sb, name);
    add(sb, config.getFitting());
    spotFitResults.resultPrefix = sb.toString();
    // Q. Should I add other fit configuration here?
    // The fraction of positive and negative candidates that were included
    add(sb, (100.0 * ctp) / candidateData.countPositive);
    add(sb, (100.0 * cfp) / candidateData.countNegative);
    // Score the fitting results compared to the original simulation.
    // Score the candidate selection:
    add(sb, ctp + cfp);
    add(sb, ctp);
    add(sb, cfp);
    // TP are all candidates that can be matched to a spot
    // FP are all candidates that cannot be matched to a spot
    // FN = The number of missed spots
    FractionClassificationResult match = new FractionClassificationResult(ctp, cfp, 0, simulationParameters.molecules - ctp);
    add(sb, match.getRecall());
    add(sb, match.getPrecision());
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // Score the fitting results:
    add(sb, failCtp);
    add(sb, failCfp);
    // TP are all fit results that can be matched to a spot
    // FP are all fit results that cannot be matched to a spot
    // FN = The number of missed spots
    add(sb, tp);
    add(sb, fp);
    match = new FractionClassificationResult(tp, fp, 0, simulationParameters.molecules - tp);
    add(sb, match.getRecall());
    add(sb, match.getPrecision());
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // Do it again but pretend we can perfectly filter all the false positives
    // add(sb, tp);
    match = new FractionClassificationResult(tp, 0, 0, simulationParameters.molecules - tp);
    // Recall is unchanged
    // Precision will be 100%
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // The mean may be subject to extreme outliers so use the median
    double median = distanceStats.getMedian();
    add(sb, median);
    final WindowOrganiser wo = new WindowOrganiser();
    String label = String.format("Recall = %s. n = %d. Median = %s nm. SD = %s nm", MathUtils.rounded(match.getRecall()), distanceStats.getN(), MathUtils.rounded(median), MathUtils.rounded(distanceStats.getStandardDeviation()));
    new HistogramPlotBuilder(TITLE, distanceStats, "Match Distance (nm)").setPlotLabel(label).show(wo);
    median = depthStats.getMedian();
    add(sb, median);
    // Sort by spot intensity and produce correlation
    double[] correlation = null;
    double[] rankCorrelation = null;
    double[] rank = null;
    final FastCorrelator fastCorrelator = new FastCorrelator();
    final ArrayList<Ranking> pc1 = new ArrayList<>();
    final ArrayList<Ranking> pc2 = new ArrayList<>();
    ci = 0;
    if (settings.showCorrelation) {
        final int[] indices = SimpleArrayUtils.natural(i1.length);
        SortUtils.sortData(indices, is, settings.rankByIntensity, true);
        correlation = new double[i1.length];
        rankCorrelation = new double[i1.length];
        rank = new double[i1.length];
        for (final int ci2 : indices) {
            fastCorrelator.add(Math.round(i1[ci2]), Math.round(i2[ci2]));
            pc1.add(new Ranking(i1[ci2], ci));
            pc2.add(new Ranking(i2[ci2], ci));
            correlation[ci] = fastCorrelator.getCorrelation();
            rankCorrelation[ci] = Correlator.correlation(rank(pc1), rank(pc2));
            if (settings.rankByIntensity) {
                rank[ci] = is[0] - is[ci];
            } else {
                rank[ci] = ci;
            }
            ci++;
        }
    } else {
        for (int i = 0; i < i1.length; i++) {
            fastCorrelator.add(Math.round(i1[i]), Math.round(i2[i]));
            pc1.add(new Ranking(i1[i], i));
            pc2.add(new Ranking(i2[i], i));
        }
    }
    final double pearsonCorr = fastCorrelator.getCorrelation();
    final double rankedCorr = Correlator.correlation(rank(pc1), rank(pc2));
    // Get the regression
    final SimpleRegression regression = new SimpleRegression(false);
    for (int i = 0; i < pc1.size(); i++) {
        regression.addData(pc1.get(i).value, pc2.get(i).value);
    }
    // final double intercept = regression.getIntercept();
    final double slope = regression.getSlope();
    if (settings.showCorrelation) {
        String title = TITLE + " Intensity";
        Plot plot = new Plot(title, "Candidate", "Spot");
        final double[] limits1 = MathUtils.limits(i1);
        final double[] limits2 = MathUtils.limits(i2);
        plot.setLimits(limits1[0], limits1[1], limits2[0], limits2[1]);
        label = String.format("Correlation=%s; Ranked=%s; Slope=%s", MathUtils.rounded(pearsonCorr), MathUtils.rounded(rankedCorr), MathUtils.rounded(slope));
        plot.addLabel(0, 0, label);
        plot.setColor(Color.red);
        plot.addPoints(i1, i2, Plot.DOT);
        if (slope > 1) {
            plot.drawLine(limits1[0], limits1[0] * slope, limits1[1], limits1[1] * slope);
        } else {
            plot.drawLine(limits2[0] / slope, limits2[0], limits2[1] / slope, limits2[1]);
        }
        ImageJUtils.display(title, plot, wo);
        title = TITLE + " Correlation";
        plot = new Plot(title, "Spot Rank", "Correlation");
        final double[] xlimits = MathUtils.limits(rank);
        double[] ylimits = MathUtils.limits(correlation);
        ylimits = MathUtils.limits(ylimits, rankCorrelation);
        plot.setLimits(xlimits[0], xlimits[1], ylimits[0], ylimits[1]);
        plot.setColor(Color.red);
        plot.addPoints(rank, correlation, Plot.LINE);
        plot.setColor(Color.blue);
        plot.addPoints(rank, rankCorrelation, Plot.LINE);
        plot.setColor(Color.black);
        plot.addLabel(0, 0, label);
        ImageJUtils.display(title, plot, wo);
    }
    add(sb, pearsonCorr);
    add(sb, rankedCorr);
    add(sb, slope);
    label = String.format("n = %d. Median = %s nm", depthStats.getN(), MathUtils.rounded(median));
    new HistogramPlotBuilder(TITLE, depthStats, "Match Depth (nm)").setRemoveOutliersOption(1).setPlotLabel(label).show(wo);
    // Plot histograms of the stats on the same window
    final double[] lower = new double[filterCriteria.length];
    final double[] upper = new double[lower.length];
    final double[] min = new double[lower.length];
    final double[] max = new double[lower.length];
    for (int i = 0; i < stats[0].length; i++) {
        final double[] limits = showDoubleHistogram(stats, i, wo, matchScores);
        lower[i] = limits[0];
        upper[i] = limits[1];
        min[i] = limits[2];
        max[i] = limits[3];
    }
    // Reconfigure some of the range limits
    // Make this a bit bigger
    upper[FILTER_SIGNAL] *= 2;
    // Make this a bit bigger
    upper[FILTER_SNR] *= 2;
    final double factor = 0.25;
    if (lower[FILTER_MIN_WIDTH] != 0) {
        // (assuming lower is less than 1)
        upper[FILTER_MIN_WIDTH] = 1 - Math.max(0, factor * (1 - lower[FILTER_MIN_WIDTH]));
    }
    if (upper[FILTER_MIN_WIDTH] != 0) {
        // (assuming upper is more than 1)
        lower[FILTER_MAX_WIDTH] = 1 + Math.max(0, factor * (upper[FILTER_MAX_WIDTH] - 1));
    }
    // Round the ranges
    final double[] interval = new double[stats[0].length];
    interval[FILTER_SIGNAL] = SignalFilter.DEFAULT_INCREMENT;
    interval[FILTER_SNR] = SnrFilter.DEFAULT_INCREMENT;
    interval[FILTER_MIN_WIDTH] = WidthFilter2.DEFAULT_MIN_INCREMENT;
    interval[FILTER_MAX_WIDTH] = WidthFilter.DEFAULT_INCREMENT;
    interval[FILTER_SHIFT] = ShiftFilter.DEFAULT_INCREMENT;
    interval[FILTER_ESHIFT] = EShiftFilter.DEFAULT_INCREMENT;
    interval[FILTER_PRECISION] = PrecisionFilter.DEFAULT_INCREMENT;
    interval[FILTER_ITERATIONS] = 0.1;
    interval[FILTER_EVALUATIONS] = 0.1;
    // Create a range increment
    final double[] increment = new double[lower.length];
    for (int i = 0; i < increment.length; i++) {
        lower[i] = MathUtils.floor(lower[i], interval[i]);
        upper[i] = MathUtils.ceil(upper[i], interval[i]);
        final double range = upper[i] - lower[i];
        // Allow clipping if the range is small compared to the min increment
        double multiples = range / interval[i];
        // Use 8 multiples for the equivalent of +/- 4 steps around the centre
        if (multiples < 8) {
            multiples = Math.ceil(multiples);
        } else {
            multiples = 8;
        }
        increment[i] = MathUtils.ceil(range / multiples, interval[i]);
        if (i == FILTER_MIN_WIDTH) {
            // Requires clipping based on the upper limit
            lower[i] = upper[i] - increment[i] * multiples;
        } else {
            upper[i] = lower[i] + increment[i] * multiples;
        }
    }
    for (int i = 0; i < stats[0].length; i++) {
        lower[i] = MathUtils.round(lower[i]);
        upper[i] = MathUtils.round(upper[i]);
        min[i] = MathUtils.round(min[i]);
        max[i] = MathUtils.round(max[i]);
        increment[i] = MathUtils.round(increment[i]);
        sb.append('\t').append(min[i]).append(':').append(lower[i]).append('-').append(upper[i]).append(':').append(max[i]);
    }
    // Disable some filters
    increment[FILTER_SIGNAL] = Double.POSITIVE_INFINITY;
    // increment[FILTER_SHIFT] = Double.POSITIVE_INFINITY;
    increment[FILTER_ESHIFT] = Double.POSITIVE_INFINITY;
    wo.tile();
    sb.append('\t').append(TextUtils.nanosToString(runTime));
    createTable().append(sb.toString());
    if (settings.saveFilterRange) {
        GUIFilterSettings filterSettings = SettingsManager.readGuiFilterSettings(0);
        String filename = (silent) ? filterSettings.getFilterSetFilename() : ImageJUtils.getFilename("Filter_range_file", filterSettings.getFilterSetFilename());
        if (filename == null) {
            return;
        }
        // Remove extension to store the filename
        filename = FileUtils.replaceExtension(filename, ".xml");
        filterSettings = filterSettings.toBuilder().setFilterSetFilename(filename).build();
        // Create a filter set using the ranges
        final ArrayList<Filter> filters = new ArrayList<>(4);
        // Create the multi-filter using the same precision type as that used during fitting.
        // Currently no support for z-filter as 3D astigmatism fitting is experimental.
        final PrecisionMethod precisionMethod = getPrecisionMethod((DirectFilter) multiFilter.getFilter());
        Function<double[], Filter> generator;
        if (precisionMethod == PrecisionMethod.POISSON_CRLB) {
            generator = parameters -> new MultiFilterCrlb(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        } else if (precisionMethod == PrecisionMethod.MORTENSEN) {
            generator = parameters -> new MultiFilter(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        } else {
            // Default
            generator = parameters -> new MultiFilter2(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        }
        filters.add(generator.apply(lower));
        filters.add(generator.apply(upper));
        filters.add(generator.apply(increment));
        if (saveFilters(filename, filters)) {
            SettingsManager.writeSettings(filterSettings);
        }
        // Create a filter set using the min/max and the initial bounds.
        // Set sensible limits
        min[FILTER_SIGNAL] = Math.max(min[FILTER_SIGNAL], 30);
        max[FILTER_SNR] = Math.min(max[FILTER_SNR], 10000);
        max[FILTER_PRECISION] = Math.min(max[FILTER_PRECISION], 100);
        // Make the 4-set filters the same as the 3-set filters.
        filters.clear();
        filters.add(generator.apply(min));
        filters.add(generator.apply(lower));
        filters.add(generator.apply(upper));
        filters.add(generator.apply(max));
        saveFilters(FileUtils.replaceExtension(filename, ".4.xml"), filters);
    }
    spotFitResults.min = min;
    spotFitResults.max = max;
}