Examples with Mean org.apache.commons.math3.stat.descriptive.moment.Mean used on opensource projects

Example 51 with Mean

use of org.apache.commons.math3.stat.descriptive.moment.Mean in project GDSC-SMLM by aherbert.

the class FIRE method showQEstimationDialog.

private boolean showQEstimationDialog(final PrecisionHistogram histogram, final QPlot qplot, final FRCCurve frcCurve, final double nmPerPixel) {
    // This is used for the initial layout of windows
    final MyWindowOrganiser wo = new MyWindowOrganiser();
    // Use a simple workflow
    Workflow<WorkSettings, Object> workflow = new Workflow<WorkSettings, Object>();
    // Split the work to two children with a dummy initial worker
    int previous = workflow.add(new BaseWorker(wo));
    workflow.add(new HistogramWorker(wo, histogram), previous);
    workflow.add(new QPlotWorker(wo, qplot), previous);
    workflow.start();
    // The number of plots
    wo.expected = 4;
    String KEY_MEAN = "mean_estimate";
    String KEY_SIGMA = "sigma_estimate";
    String KEY_Q = "q_estimate";
    String macroOptions = Macro.getOptions();
    if (macroOptions != null) {
        // If inside a macro then just get the options and run the work
        double mean = Double.parseDouble(Macro.getValue(macroOptions, KEY_MEAN, Double.toString(histogram.mean)));
        double sigma = Double.parseDouble(Macro.getValue(macroOptions, KEY_SIGMA, Double.toString(histogram.sigma)));
        double qValue = Double.parseDouble(Macro.getValue(macroOptions, KEY_Q, Double.toString(qplot.qValue)));
        workflow.run(new WorkSettings(mean, sigma, qValue));
        workflow.shutdown(false);
    } else {
        // Draw the plots with the first set of work
        workflow.run(new WorkSettings(histogram.mean, histogram.sigma, qplot.qValue));
        // Build the dialog
        NonBlockingExtendedGenericDialog gd = new NonBlockingExtendedGenericDialog(TITLE);
        gd.addHelp(About.HELP_URL);
        double mu = histogram.mean / nmPerPixel;
        double sd = histogram.sigma / nmPerPixel;
        double plateauness = qplot.computePlateauness(qplot.qValue, mu, sd);
        gd.addMessage("Estimate the blinking correction parameter Q for Fourier Ring Correlation\n \n" + String.format("Initial estimate:\nPrecision = %.3f +/- %.3f\n", histogram.mean, histogram.sigma) + String.format("Q = %s\nCost = %.3f", Utils.rounded(qplot.qValue), plateauness));
        double mean10 = histogram.mean * 10;
        double sd10 = histogram.sigma * 10;
        double q10 = qplot.qValue * 10;
        gd.addSlider("Mean (x10)", Math.max(0, mean10 - sd10 * 2), mean10 + sd10 * 2, mean10);
        gd.addSlider("Sigma (x10)", Math.max(0, sd10 / 2), sd10 * 2, sd10);
        gd.addSlider("Q (x10)", 0, Math.max(50, q10 * 2), q10);
        gd.addCheckbox("Reset_all", false);
        gd.addMessage("Double-click a slider to reset");
        gd.addDialogListener(new FIREDialogListener(gd, histogram, qplot, workflow));
        // Show this when the workers have finished drawing the plots so it is on top
        try {
            long timeout = System.currentTimeMillis() + 5000;
            while (wo.size < wo.expected) {
                Thread.sleep(50);
                if (System.currentTimeMillis() > timeout)
                    break;
            }
        } catch (InterruptedException e) {
        // Ignore
        }
        gd.showDialog();
        // Finish the worker threads
        boolean cancelled = gd.wasCanceled();
        workflow.shutdown(cancelled);
        if (cancelled)
            return false;
    }
    // Store the Q value and the mean and sigma
    qValue = qplot.qValue;
    mean = qplot.mean;
    sigma = qplot.sigma;
    // Record the values for Macros since the NonBlockingDialog doesn't
    if (Recorder.record) {
        Recorder.recordOption(KEY_MEAN, Double.toString(mean));
        Recorder.recordOption(KEY_SIGMA, Double.toString(sigma));
        Recorder.recordOption(KEY_Q, Double.toString(qValue));
    }
    return true;
}

Example 52 with Mean

use of org.apache.commons.math3.stat.descriptive.moment.Mean in project GDSC-SMLM by aherbert.

the class FIRE method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.PlugIn#run(java.lang.String)
	 */
public void run(String arg) {
    extraOptions = Utils.isExtraOptions();
    SMLMUsageTracker.recordPlugin(this.getClass(), arg);
    // Require some fit results and selected regions
    int size = MemoryPeakResults.countMemorySize();
    if (size == 0) {
        IJ.error(TITLE, "There are no fitting results in memory");
        return;
    }
    if ("q".equals(arg)) {
        TITLE += " Q estimation";
        runQEstimation();
        return;
    }
    IJ.showStatus(TITLE + " ...");
    if (!showInputDialog())
        return;
    MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
    if (results == null || results.size() == 0) {
        IJ.error(TITLE, "No results could be loaded");
        return;
    }
    MemoryPeakResults results2 = ResultsManager.loadInputResults(inputOption2, false);
    results = cropToRoi(results);
    if (results.size() < 2) {
        IJ.error(TITLE, "No results within the crop region");
        return;
    }
    if (results2 != null) {
        results2 = cropToRoi(results2);
        if (results2.size() < 2) {
            IJ.error(TITLE, "No results2 within the crop region");
            return;
        }
    }
    initialise(results, results2);
    if (!showDialog())
        return;
    long start = System.currentTimeMillis();
    // Compute FIRE
    String name = results.getName();
    double fourierImageScale = SCALE_VALUES[imageScaleIndex];
    int imageSize = IMAGE_SIZE_VALUES[imageSizeIndex];
    if (this.results2 != null) {
        name += " vs " + results2.getName();
        FireResult result = calculateFireNumber(fourierMethod, samplingMethod, thresholdMethod, fourierImageScale, imageSize);
        if (result != null) {
            logResult(name, result);
            if (showFRCCurve)
                showFrcCurve(name, result, thresholdMethod);
        }
    } else {
        FireResult result = null;
        int repeats = (randomSplit) ? Math.max(1, FIRE.repeats) : 1;
        if (repeats == 1) {
            result = calculateFireNumber(fourierMethod, samplingMethod, thresholdMethod, fourierImageScale, imageSize);
            if (result != null) {
                logResult(name, result);
                if (showFRCCurve)
                    showFrcCurve(name, result, thresholdMethod);
            }
        } else {
            // Multi-thread this ... 			
            int nThreads = Maths.min(repeats, getThreads());
            ExecutorService executor = Executors.newFixedThreadPool(nThreads);
            TurboList<Future<?>> futures = new TurboList<Future<?>>(repeats);
            TurboList<FIREWorker> workers = new TurboList<FIREWorker>(repeats);
            setProgress(repeats);
            IJ.showProgress(0);
            IJ.showStatus(TITLE + " computing ...");
            for (int i = 1; i <= repeats; i++) {
                FIREWorker w = new FIREWorker(i, fourierImageScale, imageSize);
                workers.add(w);
                futures.add(executor.submit(w));
            }
            // Wait for all to finish
            for (int t = futures.size(); t-- > 0; ) {
                try {
                    // The future .get() method will block until completed
                    futures.get(t).get();
                } catch (Exception e) {
                    // This should not happen. 
                    // Ignore it and allow processing to continue (the number of neighbour samples will just be smaller).  
                    e.printStackTrace();
                }
            }
            IJ.showProgress(1);
            executor.shutdown();
            // Show a combined FRC curve plot of all the smoothed curves if we have multiples.
            LUT valuesLUT = LUTHelper.createLUT(LutColour.FIRE_GLOW);
            @SuppressWarnings("unused") LUT // Black at max value
            noSmoothLUT = LUTHelper.createLUT(LutColour.GRAYS).createInvertedLut();
            LUTHelper.DefaultLUTMapper mapper = new LUTHelper.DefaultLUTMapper(0, repeats);
            FrcCurve curve = new FrcCurve();
            Statistics stats = new Statistics();
            WindowOrganiser wo = new WindowOrganiser();
            boolean oom = false;
            for (int i = 0; i < repeats; i++) {
                FIREWorker w = workers.get(i);
                if (w.oom)
                    oom = true;
                if (w.result == null)
                    continue;
                result = w.result;
                if (!Double.isNaN(result.fireNumber))
                    stats.add(result.fireNumber);
                if (showFRCCurveRepeats) {
                    // Output each FRC curve using a suffix.
                    logResult(w.name, result);
                    wo.add(Utils.display(w.plot.getTitle(), w.plot));
                }
                if (showFRCCurve) {
                    int index = mapper.map(i + 1);
                    //@formatter:off
                    curve.add(name, result, thresholdMethod, LUTHelper.getColour(valuesLUT, index), Color.blue, //LUTHelper.getColour(noSmoothLUT, index)
                    null);
                //@formatter:on
                }
            }
            if (result != null) {
                wo.cascade();
                double mean = stats.getMean();
                logResult(name, result, mean, stats);
                if (showFRCCurve) {
                    curve.addResolution(mean);
                    Plot2 plot = curve.getPlot();
                    Utils.display(plot.getTitle(), plot);
                }
            }
            if (oom) {
                //@formatter:off
                IJ.error(TITLE, "ERROR - Parallel computation out-of-memory.\n \n" + TextUtils.wrap("The number of results will be reduced. " + "Please reduce the size of the Fourier image " + "or change the number of threads " + "using the extra options (hold down the 'Shift' " + "key when running the plugin).", 80));
            //@formatter:on
            }
        }
        // Only do this once
        if (showFRCTimeEvolution && result != null && !Double.isNaN(result.fireNumber))
            showFrcTimeEvolution(name, result.fireNumber, thresholdMethod, nmPerPixel / result.getNmPerPixel(), imageSize);
    }
    IJ.showStatus(TITLE + " complete : " + Utils.timeToString(System.currentTimeMillis() - start));
}

Example 53 with Mean

use of org.apache.commons.math3.stat.descriptive.moment.Mean in project GDSC-SMLM by aherbert.

the class FIRE method calculatePrecisionHistogram.

/**
	 * Calculate a histogram of the precision. The precision can be either stored in the results or calculated using the
	 * Mortensen formula. If the precision method for Q estimation is not fixed then the histogram is fitted with a
	 * Gaussian to create an initial estimate.
	 *
	 * @param precisionMethod
	 *            the precision method
	 * @return The precision histogram
	 */
private PrecisionHistogram calculatePrecisionHistogram() {
    boolean logFitParameters = false;
    String title = results.getName() + " Precision Histogram";
    // Check if the results has the precision already or if it can be computed.
    boolean canUseStored = canUseStoredPrecision(results);
    boolean canCalculatePrecision = canCalculatePrecision(results);
    // Set the method to compute a histogram. Default to the user selected option.
    PrecisionMethod m = null;
    if (canUseStored && precisionMethod == PrecisionMethod.STORED)
        m = precisionMethod;
    else if (canCalculatePrecision && precisionMethod == PrecisionMethod.CALCULATE)
        m = precisionMethod;
    if (m == null) {
        // We only have two choices so if one is available then select it.
        if (canUseStored)
            m = PrecisionMethod.STORED;
        else if (canCalculatePrecision)
            m = PrecisionMethod.CALCULATE;
        // If the user selected a method not available then log a warning
        if (m != null && precisionMethod != PrecisionMethod.FIXED) {
            IJ.log(String.format("%s : Selected precision method '%s' not available, switching to '%s'", TITLE, precisionMethod, m.getName()));
        }
        if (m == null) {
            // This does not matter if the user has provide a fixed input.
            if (precisionMethod == PrecisionMethod.FIXED) {
                PrecisionHistogram histogram = new PrecisionHistogram(title);
                histogram.mean = mean;
                histogram.sigma = sigma;
                return histogram;
            }
            // No precision
            return null;
        }
    }
    // We get here if we can compute precision.
    // Build the histogram 
    StoredDataStatistics precision = new StoredDataStatistics(results.size());
    if (m == PrecisionMethod.STORED) {
        for (PeakResult r : results.getResults()) {
            precision.add(r.getPrecision());
        }
    } else {
        final double nmPerPixel = results.getNmPerPixel();
        final double gain = results.getGain();
        final boolean emCCD = results.isEMCCD();
        for (PeakResult r : results.getResults()) {
            precision.add(r.getPrecision(nmPerPixel, gain, emCCD));
        }
    }
    //System.out.printf("Raw p = %f\n", precision.getMean());
    double yMin = Double.NEGATIVE_INFINITY, yMax = 0;
    // Set the min and max y-values using 1.5 x IQR 
    DescriptiveStatistics stats = precision.getStatistics();
    double lower = stats.getPercentile(25);
    double upper = stats.getPercentile(75);
    if (Double.isNaN(lower) || Double.isNaN(upper)) {
        if (logFitParameters)
            Utils.log("Error computing IQR: %f - %f", lower, upper);
    } else {
        double iqr = upper - lower;
        yMin = FastMath.max(lower - iqr, stats.getMin());
        yMax = FastMath.min(upper + iqr, stats.getMax());
        if (logFitParameters)
            Utils.log("  Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax);
    }
    if (yMin == Double.NEGATIVE_INFINITY) {
        int n = 5;
        yMin = Math.max(stats.getMin(), stats.getMean() - n * stats.getStandardDeviation());
        yMax = Math.min(stats.getMax(), stats.getMean() + n * stats.getStandardDeviation());
        if (logFitParameters)
            Utils.log("  Data range: %f - %f. Plotting mean +/- %dxSD: %f - %f", stats.getMin(), stats.getMax(), n, yMin, yMax);
    }
    // Get the data within the range
    double[] data = precision.getValues();
    precision = new StoredDataStatistics(data.length);
    for (double d : data) {
        if (d < yMin || d > yMax)
            continue;
        precision.add(d);
    }
    int histogramBins = Utils.getBins(precision, Utils.BinMethod.SCOTT);
    float[][] hist = Utils.calcHistogram(precision.getFloatValues(), yMin, yMax, histogramBins);
    PrecisionHistogram histogram = new PrecisionHistogram(hist, precision, title);
    if (precisionMethod == PrecisionMethod.FIXED) {
        histogram.mean = mean;
        histogram.sigma = sigma;
        return histogram;
    }
    // Fitting of the histogram to produce the initial estimate
    // Extract non-zero data
    float[] x = Arrays.copyOf(hist[0], hist[0].length);
    float[] y = hist[1];
    int count = 0;
    float dx = (x[1] - x[0]) * 0.5f;
    for (int i = 0; i < y.length; i++) if (y[i] > 0) {
        x[count] = x[i] + dx;
        y[count] = y[i];
        count++;
    }
    x = Arrays.copyOf(x, count);
    y = Arrays.copyOf(y, count);
    // Sense check to fitted data. Get mean and SD of histogram
    double[] stats2 = Utils.getHistogramStatistics(x, y);
    if (logFitParameters)
        Utils.log("  Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
    histogram.mean = stats2[0];
    histogram.sigma = stats2[1];
    // Standard Gaussian fit
    double[] parameters = fitGaussian(x, y);
    if (parameters == null) {
        Utils.log("  Failed to fit initial Gaussian");
        return histogram;
    }
    double newMean = parameters[1];
    double error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        Utils.log("  Failed to fit Gaussian: %f standard deviations from histogram mean", error);
        return histogram;
    }
    if (newMean < yMin || newMean > yMax) {
        Utils.log("  Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return histogram;
    }
    if (logFitParameters)
        Utils.log("  Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
    histogram.mean = parameters[1];
    histogram.sigma = parameters[2];
    return histogram;
}

Example 54 with Mean

use of org.apache.commons.math3.stat.descriptive.moment.Mean in project GDSC-SMLM by aherbert.

the class FIRE method fitGaussian.

/**
	 * Fit gaussian.
	 *
	 * @param x
	 *            the x
	 * @param y
	 *            the y
	 * @return new double[] { norm, mean, sigma }
	 */
private double[] fitGaussian(float[] x, float[] y) {
    WeightedObservedPoints obs = new WeightedObservedPoints();
    for (int i = 0; i < x.length; i++) obs.add(x[i], y[i]);
    Collection<WeightedObservedPoint> observations = obs.toList();
    GaussianCurveFitter fitter = GaussianCurveFitter.create().withMaxIterations(2000);
    GaussianCurveFitter.ParameterGuesser guess = new GaussianCurveFitter.ParameterGuesser(observations);
    double[] initialGuess = null;
    try {
        initialGuess = guess.guess();
        return fitter.withStartPoint(initialGuess).fit(observations);
    } catch (TooManyEvaluationsException e) {
        // Use the initial estimate
        return initialGuess;
    } catch (Exception e) {
        // Just in case there is another exception type, or the initial estimate failed
        return null;
    }
}

Example 55 with Mean

use of org.apache.commons.math3.stat.descriptive.moment.Mean in project GDSC-SMLM by aherbert.

the class FIRE method computeExpDecay.

private double[] computeExpDecay(double mean, double sigma, double[] q) {
    double[] hq = FRC.computeHq(q, mean, sigma);
    double[] exp_decay = new double[q.length];
    exp_decay[0] = 1;
    for (int i = 1; i < q.length; i++) {
        double sinc_q = sinc(Math.PI * q[i]);
        exp_decay[i] = sinc_q * sinc_q * hq[i];
    }
    return exp_decay;
}