Search in sources :

Example 11 with UnivariateFunction

use of org.apache.commons.math3.analysis.UnivariateFunction 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;
    UnivariateIntegrator integrator = new SimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
    UnivariateFunction f = new UnivariateFunction() {

        public double value(double x) {
            //return AiryPattern.intensity(x) * 2 * Math.PI * x / (4 * Math.PI);
            return AiryPattern.intensity(x) * 0.5 * x;
        }
    };
    // Integrate up to a set number of dark rings
    int samples = 1000;
    final double step = RINGS[SAMPLE_RINGS] / samples;
    double to = 0, from = 0;
    r = new double[samples + 1];
    sum = new double[samples + 1];
    for (int i = 1; i < sum.length; i++) {
        from = to;
        r[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 RuntimeException("Failed to create the Airy distribution");
    SplineInterpolator si = new SplineInterpolator();
    spline = si.interpolate(sum, r);
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) UnivariateIntegrator(org.apache.commons.math3.analysis.integration.UnivariateIntegrator) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction) SplineInterpolator(org.apache.commons.math3.analysis.interpolation.SplineInterpolator)

Example 12 with UnivariateFunction

use of org.apache.commons.math3.analysis.UnivariateFunction in project GDSC-SMLM by aherbert.

the class PoissonFunctionTest method cumulativeProbabilityIsOneWithRealAbove4.

private void cumulativeProbabilityIsOneWithRealAbove4(final double gain, final double mu, int min, int max) {
    final double o = mu * gain;
    final PoissonFunction f = new PoissonFunction(1.0 / gain, true);
    double p = 0;
    SimpsonIntegrator in = new SimpsonIntegrator(3, 30);
    try {
        p = in.integrate(Integer.MAX_VALUE, new UnivariateFunction() {

            public double value(double x) {
                return f.likelihood(x, o);
            }
        }, min, max);
        System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
        Assert.assertEquals(String.format("g=%f, mu=%f", gain, mu), 1, p, 0.02);
    } catch (TooManyEvaluationsException e) {
        //double inc = max / 20000.0;
        //for (double x = 0; x <= max; x += inc)
        //{
        //	final double pp = f.likelihood(x, o);
        //	//System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, pp);
        //	p += pp;
        //}
        //System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
        Assert.assertFalse(e.getMessage(), true);
    }
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) TooManyEvaluationsException(org.apache.commons.math3.exception.TooManyEvaluationsException) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction)

Example 13 with UnivariateFunction

use of org.apache.commons.math3.analysis.UnivariateFunction in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method cumulativeProbabilityIsOneWithRealData.

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

        public double value(double x) {
            double v;
            v = PoissonCalculator.likelihood(mu, x);
            //System.out.printf("x=%f, v=%f\n", x, v);
            return v;
        }
    }, min, max);
    System.out.printf("mu=%f, p=%f\n", mu, p);
    if (test) {
        Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
    }
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) UnivariateIntegrator(org.apache.commons.math3.analysis.integration.UnivariateIntegrator) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction)

Example 14 with UnivariateFunction

use of org.apache.commons.math3.analysis.UnivariateFunction in project GDSC-SMLM by aherbert.

the class SCMOSLikelihoodWrapperTest method cumulativeProbabilityIsOneWithRealData.

private void cumulativeProbabilityIsOneWithRealData(final double mu, double min, double max, boolean test) {
    double p = 0;
    // Test using a standard Poisson-Gaussian convolution
    //min = -max;
    //final PoissonGaussianFunction pgf = PoissonGaussianFunction.createWithVariance(1, 1, VAR);
    UnivariateIntegrator in = new SimpsonIntegrator();
    p = in.integrate(20000, new UnivariateFunction() {

        public double value(double x) {
            double v;
            v = SCMOSLikelihoodWrapper.likelihood(mu, VAR, G, O, x);
            //System.out.printf("x=%f, v=%f\n", x, v);
            return v;
        }
    }, min, max);
    //System.out.printf("mu=%f, p=%f\n", mu, p);
    if (test) {
        Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
    }
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) UnivariateIntegrator(org.apache.commons.math3.analysis.integration.UnivariateIntegrator) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction)

Example 15 with UnivariateFunction

use of org.apache.commons.math3.analysis.UnivariateFunction 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);
}
Also used : SimpsonIntegrator(org.apache.commons.math3.analysis.integration.SimpsonIntegrator) UnivariateFunction(org.apache.commons.math3.analysis.UnivariateFunction)

Aggregations

UnivariateFunction (org.apache.commons.math3.analysis.UnivariateFunction)17 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)7 SimpsonIntegrator (org.apache.commons.math3.analysis.integration.SimpsonIntegrator)6 NoBracketingException (org.apache.commons.math3.exception.NoBracketingException)6 FastMath (org.apache.commons.math3.util.FastMath)6 VisibleForTesting (com.google.common.annotations.VisibleForTesting)4 List (java.util.List)4 Collectors (java.util.stream.Collectors)4 IntStream (java.util.stream.IntStream)4 Nonnull (javax.annotation.Nonnull)4 Nullable (javax.annotation.Nullable)4 ImmutablePair (org.apache.commons.lang3.tuple.ImmutablePair)4 ImmutableTriple (org.apache.commons.lang3.tuple.ImmutableTriple)4 UnivariateIntegrator (org.apache.commons.math3.analysis.integration.UnivariateIntegrator)4 RobustBrentSolver (org.broadinstitute.hellbender.tools.coveragemodel.math.RobustBrentSolver)4 Utils (org.broadinstitute.hellbender.utils.Utils)4 ParamUtils (org.broadinstitute.hellbender.utils.param.ParamUtils)4 INDArray (org.nd4j.linalg.api.ndarray.INDArray)4 Nd4j (org.nd4j.linalg.factory.Nd4j)4 NDArrayIndex (org.nd4j.linalg.indexing.NDArrayIndex)4