Examples with CustomSimpsonIntegrator uk.ac.sussex.gdsc.smlm.math3.analysis.integration.CustomSimpsonIntegrator used on opensource projects

Example 1 with CustomSimpsonIntegrator

use of uk.ac.sussex.gdsc.smlm.math3.analysis.integration.CustomSimpsonIntegrator in project GDSC-SMLM by aherbert.

the class CameraModelAnalysis method cumulativeDistribution.

private static double[][] cumulativeDistribution(CameraModelAnalysisSettings settings, double[][] cdf, final LikelihoodFunction fun) {
    // Q. How to match this is the discrete cumulative histogram using the continuous
    // likelihood function:
    // 1. Compute integral up to the value
    // 2. Compute integral up to but not including the next value using trapezoid integration
    // 3. Compute integral up to but not including the next value using flat-top integration
    // Since the function will be used on continuous float data when fitting PSFs the best
    // match for how it will perform in practice is a continuous (trapezoid) integration.
    // The simplest is a flat-top integration.
    // Compute the probability at each value
    final double e = settings.getPhotons();
    double[] x = cdf[0];
    double[] y = new double[x.length];
    for (int i = 0; i < x.length; i++) {
        y[i] = fun.likelihood(x[i], e);
    }
    // Add more until the probability change is marginal
    double sum = MathUtils.sum(y);
    final TDoubleArrayList list = new TDoubleArrayList(y);
    for (int o = (int) x[x.length - 1] + 1; ; o++) {
        final double p = fun.likelihood(o, e);
        list.add(p);
        if (p == 0) {
            break;
        }
        sum += p;
        if (p / sum < PROBABILITY_DELTA) {
            break;
        }
    }
    final TDoubleArrayList list2 = new TDoubleArrayList(10);
    for (int o = (int) x[0] - 1; ; o--) {
        final double p = fun.likelihood(o, e);
        list2.add(p);
        if (p == 0) {
            break;
        }
        sum += p;
        if (p / sum < PROBABILITY_DELTA) {
            break;
        }
    }
    // Insert at start
    double start = x[0];
    if (!list2.isEmpty()) {
        start -= list2.size();
        list2.reverse();
        list.insert(0, list2.toArray());
    }
    y = list.toArray();
    x = SimpleArrayUtils.newArray(y.length, start, 1.0);
    if (settings.getSimpsonIntegration()) {
        int c0 = -1;
        double dirac = 0;
        int minc = 0;
        int maxc = 0;
        CachingUnivariateFunction uf = null;
        if (settings.getMode() == MODE_EM_CCD && isPoissonGammaLikelihoodFunction(settings)) {
            // A spike is expected at c=0 due to the Dirac delta contribution.
            // This breaks integration, especially when noise < 0.1.
            // Fix by integrating around c=0 fully then integrating the rest.
            c0 = Arrays.binarySearch(x, 0);
            final double noise = getReadNoise(settings);
            final double p = settings.getPhotons();
            if (noise == 0) {
                // Pure Poisson-Gamma. Just subtract the delta, do the simple integration
                // below and add the delta back. Only functions that support noise==0
                // will be allowed so this solution works.
                dirac = PoissonGammaFunction.dirac(p);
                if (c0 != -1) {
                    y[c0] -= dirac;
                }
            } else {
                // Fix integration around c=0 using the range of the Gaussian
                minc = (int) Math.max(x[0], Math.floor(-5 * noise));
                maxc = (int) Math.min(x[x.length - 1], Math.ceil(5 * noise));
                uf = new CachingUnivariateFunction(fun, p);
            }
        }
        // Use Simpson's integration with n=4 to get the integral of the probability
        // over the range of each count.
        // Note the Poisson-Gamma function cannot be integrated with the
        // Dirac delta function at c==0
        // Compute the extra function points
        final double[] f = new double[y.length * SIMPSON_N + 1];
        int index = f.length;
        final int mod;
        if (settings.getRoundDown()) {
            // Do this final value outside the loop as y[index/n] does not exist
            mod = 0;
            index--;
            f[index] = fun.likelihood(start + index * SIMPSON_H, e);
        } else {
            // Used when computing for rounding up/down
            mod = SIMPSON_N_2;
            start -= SIMPSON_N_2 * SIMPSON_H;
        }
        while (index-- > 0) {
            if (index % SIMPSON_N == mod) {
                f[index] = y[index / SIMPSON_N];
            } else {
                f[index] = fun.likelihood(start + index * SIMPSON_H, e);
            }
        }
        // Compute indices for the integral
        final TIntArrayList tmp = new TIntArrayList(SIMPSON_N);
        for (int j = 1; j <= SIMPSON_N_2 - 1; j++) {
            tmp.add(2 * j);
        }
        final int[] i2 = tmp.toArray();
        tmp.resetQuick();
        for (int j = 1; j <= SIMPSON_N_2; j++) {
            tmp.add(2 * j - 1);
        }
        final int[] i4 = tmp.toArray();
        // Compute integral
        for (int i = 0; i < y.length; i++) {
            final int a = i * SIMPSON_N;
            final int b = a + SIMPSON_N;
            sum = f[a] + f[b];
            for (int j = i2.length; j-- > 0; ) {
                sum += 2 * f[a + i2[j]];
            }
            for (int j = i4.length; j-- > 0; ) {
                sum += 4 * f[a + i4[j]];
            }
            sum *= SIMPSON_H / 3;
            // System.out.printf("y[%d] = %f => %f\n", i, y[i], s);
            y[i] = sum;
        }
        // Fix Poisson-Gamma ...
        if (c0 != -1) {
            if (uf != null) {
                // Convolved Poisson-Gamma. Fix in the range of the Gaussian around c=0
                final CustomSimpsonIntegrator in = new CustomSimpsonIntegrator(SIMPSON_RELATIVE_ACCURACY, SIMPSON_ABSOLUTE_ACCURACY, SIMPSON_MIN_ITERATION_COUNT, CustomSimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
                final double lower = (settings.getRoundDown()) ? 0 : -0.5;
                final double upper = lower + 1;
                // Switch from c<=maxc to c<maxc. Avoid double computation at minc==maxc
                if (maxc != minc) {
                    maxc++;
                }
                maxc++;
                for (int c = minc, i = Arrays.binarySearch(x, minc); c < maxc; c++, i++) {
                    uf.reset();
                    try {
                        y[i] = in.integrate(2000, uf, c + lower, c + upper);
                    } catch (final TooManyEvaluationsException ex) {
                        // Q. Is the last sum valid?
                        if (in.getLastSum() > 0) {
                            y[i] = in.getLastSum();
                        } else {
                            // Otherwise use all the cached values to compute a sum
                            // using the trapezoid rule. This will underestimate the sum.
                            // Note: The Simpson integrator will have computed the edge values
                            // as the first two values in the cache.
                            final double[] g = uf.list.toArray();
                            final double dx = (g[3] - g[1]) / in.getN();
                            final int total = 1 + 2 * ((int) in.getN());
                            sum = 0;
                            for (int j = 4; j < total; j += 2) {
                                sum += g[j];
                            }
                            y[i] = (g[0] + g[2] + 2 * sum) / dx;
                        }
                    }
                }
            } else {
                // Pure Poisson-Gamma. Just add back the delta.
                y[c0] += dirac;
            }
        }
    }
    // Simple flat-top integration
    sum = 0;
    for (int i = 0; i < y.length; i++) {
        sum += y[i];
        y[i] = sum;
    }
    return new double[][] { x, y };
}

Example 2 with CustomSimpsonIntegrator

use of uk.ac.sussex.gdsc.smlm.math3.analysis.integration.CustomSimpsonIntegrator in project GDSC-SMLM by aherbert.

the class AiryPsfModel method createAiryDistribution.

private static synchronized void createAiryDistribution() {
    if (spline != null) {
        return;
    }
    final double relativeAccuracy = 1e-4;
    final double absoluteAccuracy = 1e-8;
    final int minimalIterationCount = 3;
    final int maximalIterationCount = 32;
    final UnivariateIntegrator integrator = new CustomSimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
    // The pattern profile is in one dimension.
    // Multiply by the perimeter of a circle to convert to 2D volume then normalise by 4 pi
    // return AiryPattern.intensity(x) * 2 * Math.PI * x / (4 * Math.PI)
    final UnivariateFunction f = x -> AiryPattern.intensity(x) * 0.5 * x;
    // Integrate up to a set number of dark rings
    final int samples = 1000;
    final double step = RINGS[SAMPLE_RINGS] / samples;
    double to = 0;
    double from = 0;
    final double[] radius = new double[samples + 1];
    final double[] sum = new double[samples + 1];
    for (int i = 1; i < sum.length; i++) {
        from = to;
        radius[i] = to = step * i;
        sum[i] = integrator.integrate(2000, f, from, to) + sum[i - 1];
    }
    if (DoubleEquality.relativeError(sum[samples], POWER[SAMPLE_RINGS]) > 1e-3) {
        throw new ComputationException("Failed to create the Airy distribution");
    }
    final SplineInterpolator si = new SplineInterpolator();
    spline = si.interpolate(sum, radius);
}

Example 3 with CustomSimpsonIntegrator

use of uk.ac.sussex.gdsc.smlm.math3.analysis.integration.CustomSimpsonIntegrator in project GDSC-SMLM by aherbert.

the class PoissonGammaGaussianFunction method createIntegrator.

private UnivariateIntegrator createIntegrator() {
    UnivariateIntegrator in = integrator;
    if (in == null) {
        // This is the integrator for the Poisson-Gamma when observed count x>=1
        // i.e. not at the boundary x=0.
        final double relativeAccuracy = 1e-4;
        final double absoluteAccuracy = 1e-16;
        int minimalIterationCount;
        switch(convolutionMode) {
            case SIMPSON_PDF:
                // This is a CustomSimpsonIntegrator that computes 1 refinement
                // on the first iteration.
                // Number of function evaluations = 2^(iteration+1) + 1
                // => 5 for 1 iterations
                // => 9 for 2 iterations
                minimalIterationCount = 1;
                in = new CustomSimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minimalIterationCount, CustomSimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT);
                break;
            case LEGENDRE_GAUSS_PDF:
                // Not sure how to configure this.
                // The integration points are used for each sub-interval.
                // Function evaluations = integrationpoints * intervals.
                // The intervals start at 1,2 and increase by at least 4 at each stage after that.
                // At least 1 stage is done thus 3 * integrationpoints functions evaluations
                // will be done for minimalIterationCount=1.
                minimalIterationCount = 1;
                final int maximalIterationCount = 32;
                final int integrationpoints = 8;
                in = new IterativeLegendreGaussIntegrator(integrationpoints, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
                break;
            default:
                throw new IllegalStateException();
        }
        integrator = in;
    }
    return in;
}