Examples with Max org.apache.commons.math3.stat.descriptive.rank.Max used on opensource projects

Example 51 with Max

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

the class SCMOSLikelihoodWrapperTest method instanceLikelihoodMatches.

private void instanceLikelihoodMatches(final double mu, boolean test) {
    // Determine upper limit for a Poisson
    int limit = new PoissonDistribution(mu).inverseCumulativeProbability(P_LIMIT);
    // Map to observed values using the gain and offset
    double max = limit * G;
    double step = 0.1;
    int n = (int) Math.ceil(max / step);
    // Evaluate all values from (zero+offset) to large n
    double[] k = Utils.newArray(n, O, step);
    double[] a = new double[0];
    double[] gradient = new double[0];
    float[] var = newArray(n, VAR);
    float[] g = newArray(n, G);
    float[] o = newArray(n, O);
    NonLinearFunction nlf = new NonLinearFunction() {

        public void initialise(double[] a) {
        }

        public int[] gradientIndices() {
            return new int[0];
        }

        public double eval(int x, double[] dyda, double[] w) {
            return 0;
        }

        public double eval(int x) {
            return mu;
        }

        public double eval(int x, double[] dyda) {
            return mu;
        }

        public boolean canComputeWeights() {
            return false;
        }

        public double evalw(int x, double[] w) {
            return 0;
        }

        public int getNumberOfGradients() {
            return 0;
        }
    };
    SCMOSLikelihoodWrapper f = new SCMOSLikelihoodWrapper(nlf, a, k, n, var, g, o);
    double total = 0, p = 0;
    double maxp = 0;
    int maxi = 0;
    for (int i = 0; i < n; i++) {
        double nll = f.computeLikelihood(i);
        double nll2 = f.computeLikelihood(gradient, i);
        double nll3 = SCMOSLikelihoodWrapper.negativeLogLikelihood(mu, var[i], g[i], o[i], k[i]);
        total += nll;
        Assert.assertEquals("computeLikelihood @" + i, nll3, nll, nll * 1e-10);
        Assert.assertEquals("computeLikelihood+gradient @" + i, nll3, nll2, nll * 1e-10);
        double pp = FastMath.exp(-nll);
        if (maxp < pp) {
            maxp = pp;
            maxi = i;
        //System.out.printf("mu=%f, e=%f, k=%f, pp=%f\n", mu, mu * G + O, k[i], pp);
        }
        p += pp * step;
    }
    // Expected max of the distribution is the mode of the Poisson distribution.
    // This has two modes for integer input counts. We take the mean of those.
    // https://en.wikipedia.org/wiki/Poisson_distribution
    // Note that the shift of VAR/(G*G) is a constant applied to both the expected and
    // observed values and consequently cancels when predicting the max, i.e. we add
    // a constant count to the observed values and shift the distribution by the same
    // constant. We can thus compute the mode for the unshifted distribution.
    double lambda = mu;
    double mode1 = Math.floor(lambda);
    double mode2 = Math.ceil(lambda) - 1;
    // Scale to observed values
    double kmax = ((mode1 + mode2) * 0.5) * G + O;
    //System.out.printf("mu=%f, p=%f, maxp=%f @ %f  (expected=%f  %f)\n", mu, p, maxp, k[maxi], kmax, kmax - k[maxi]);
    Assert.assertEquals("k-max", kmax, k[maxi], kmax * 1e-3);
    if (test) {
        Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
    }
    // Check the function can compute the same total
    double delta = total * 1e-10;
    double sum, sum2;
    sum = f.computeLikelihood();
    sum2 = f.computeLikelihood(gradient);
    Assert.assertEquals("computeLikelihood", total, sum, delta);
    Assert.assertEquals("computeLikelihood with gradient", total, sum2, delta);
    // Check the function can compute the same total after duplication
    f = f.build(nlf, a);
    sum = f.computeLikelihood();
    sum2 = f.computeLikelihood(gradient);
    Assert.assertEquals("computeLikelihood", total, sum, delta);
    Assert.assertEquals("computeLikelihood with gradient", total, sum2, delta);
}

Example 52 with Max

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

the class PoissonLikelihoodWrapperTest method cumulativeProbabilityIsOneWithRealData.

private void cumulativeProbabilityIsOneWithRealData(final double mu, int max) {
    double p = 0;
    SimpsonIntegrator in = new SimpsonIntegrator();
    p = in.integrate(20000, new UnivariateFunction() {

        public double value(double x) {
            return PoissonLikelihoodWrapper.likelihood(mu, x);
        }
    }, 0, max);
    System.out.printf("mu=%f, p=%f\n", mu, p);
    Assert.assertEquals(String.format("mu=%f", mu), 1, p, 0.02);
}

Example 53 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max 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 54 with Max

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

the class EMGainAnalysis method simulateFromPoissonGammaGaussian.

/**
	 * Randomly generate a histogram from poisson-gamma-gaussian samples
	 * 
	 * @return The histogram
	 */
private int[] simulateFromPoissonGammaGaussian() {
    // Randomly sample
    RandomGenerator random = new Well44497b(System.currentTimeMillis() + System.identityHashCode(this));
    PoissonDistribution poisson = new PoissonDistribution(random, _photons, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    CustomGammaDistribution gamma = new CustomGammaDistribution(random, _photons, _gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    final int steps = simulationSize;
    int[] sample = new int[steps];
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0)
            IJ.showProgress(n, steps);
        // Poisson
        double d = poisson.sample();
        // Gamma
        if (d > 0) {
            gamma.setShapeUnsafe(d);
            d = gamma.sample();
        }
        // Gaussian
        d += _noise * random.nextGaussian();
        // Convert the sample to a count 
        sample[n] = (int) Math.round(d + _bias);
    }
    int max = Maths.max(sample);
    int[] h = new int[max + 1];
    for (int s : sample) h[s]++;
    return h;
}

Example 55 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max 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;
}