Search in sources :

Example 31 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class PoissonGaussianFunction2Test method probabilityMatchesLogProbability.

private static void probabilityMatchesLogProbability(final double gain, double mu, final double sd, final boolean usePicard) {
    // Note: The input s parameter is pre-gain.
    final PoissonGaussianFunction2 f = PoissonGaussianFunction2.createWithStandardDeviation(1.0 / gain, sd * gain);
    f.setUsePicardApproximation(usePicard);
    // Evaluate an initial range.
    // Gaussian should have >99% within +/- s
    // Poisson will have mean mu with a variance mu.
    // At large mu it is approximately normal so use 3 sqrt(mu) for the range added to the mean
    final int[] range = PoissonGaussianFunctionTest.getRange(gain, mu, sd);
    final int min = range[0];
    final int max = range[1];
    // Note: The input mu parameter is pre-gain.
    final double e = mu;
    final Supplier<String> msg = () -> String.format("g=%f, mu=%f, s=%f", gain, mu, sd);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-3, 0);
    for (int x = min; x <= max; x++) {
        final double p = f.likelihood(x, e);
        if (p == 0) {
            continue;
        }
        final double logP = f.logLikelihood(x, e);
        TestAssertions.assertTest(Math.log(p), logP, predicate, msg);
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)

Example 32 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class PoissonLikelihoodWrapperTest method cumulativeProbabilityIsOneFromLikelihood.

private void cumulativeProbabilityIsOneFromLikelihood(final double mu) {
    // Determine upper limit for a Poisson
    final int limit = new PoissonDistribution(mu).inverseCumulativeProbability(0.999);
    // Expand to allow for the gain
    final int n = (int) Math.ceil(limit / alpha);
    // Evaluate all values from zero to large n
    final double[] k = SimpleArrayUtils.newArray(n, 0, 1.0);
    final double[] a = new double[1];
    final double[] g = new double[1];
    final NonLinearFunction nlf = new DummyNonLinearFunction(mu);
    PoissonLikelihoodWrapper func = new PoissonLikelihoodWrapper(nlf, a, Arrays.copyOf(k, n), n, alpha);
    // Keep evaluating up until no difference
    final double changeTolerance = 1e-6;
    double total = 0;
    double pvalue = 0;
    int index = 0;
    while (index < n) {
        final double nll = func.computeLikelihood(index);
        final double nll2 = func.computeLikelihood(g, index);
        index++;
        Assertions.assertEquals(nll, nll2, 1e-10, "computeLikelihood(i)");
        total += nll;
        final double pp = StdMath.exp(-nll);
        // logger.fine(FunctionUtils.getSupplier("mu=%f, o=%f, i=%d, pp=%f", mu, mu / alpha, i, pp);
        pvalue += pp;
        if (pvalue > 0.5 && pp / pvalue < changeTolerance) {
            break;
        }
    }
    logger.log(TestLogUtils.getRecord(LOG_LEVEL, "mu=%f, limit=%d, p=%f", mu, limit, pvalue));
    Assertions.assertEquals(1, pvalue, 0.02, () -> String.format("mu=%f", mu));
    // Check the function can compute the same total
    func = new PoissonLikelihoodWrapper(nlf, a, k, index, alpha);
    final double sum = func.computeLikelihood();
    final double sum2 = func.computeLikelihood(g);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    TestAssertions.assertTest(total, sum, predicate, "computeLikelihood");
    TestAssertions.assertTest(total, sum2, predicate, "computeLikelihood with gradient");
}
Also used : PoissonDistribution(org.apache.commons.math3.distribution.PoissonDistribution) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)

Example 33 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionCanComputeIntegralWith2Peaks.

@Test
void functionCanComputeIntegralWith2Peaks() {
    Assumptions.assumeTrue(null != f2);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-8, 0);
    double[] params;
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        for (final double[] w1 : testw1) {
                            for (final double angle1 : testangle1) {
                                // Peak 2
                                for (final double signal2 : testsignal2) {
                                    for (final double cx2 : testcx2) {
                                        for (final double cy2 : testcy2) {
                                            for (final double cz2 : testcz2) {
                                                for (final double[] w2 : testw2) {
                                                    for (final double angle2 : testangle2) {
                                                        params = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1, signal2, cx2, cy2, cz2, w2[0], w2[1], angle2);
                                                        final double e = new IntegralValueProcedure().getIntegral(f2, params);
                                                        final double o = f2.integral(params);
                                                        TestAssertions.assertTest(e, o, predicate);
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) IntegralValueProcedure(uk.ac.sussex.gdsc.smlm.function.IntegralValueProcedure) Gaussian2DFunctionTest(uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunctionTest) Test(org.junit.jupiter.api.Test)

Example 34 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionCanComputeIntegral.

@Test
void functionCanComputeIntegral() {
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-8, 0);
    double[] params;
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        for (final double[] w1 : testw1) {
                            for (final double angle1 : testangle1) {
                                params = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1);
                                final double e = new IntegralValueProcedure().getIntegral(f1, params);
                                final double o = f1.integral(params);
                                TestAssertions.assertTest(e, o, predicate);
                            }
                        }
                    }
                }
            }
        }
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) IntegralValueProcedure(uk.ac.sussex.gdsc.smlm.function.IntegralValueProcedure) Gaussian2DFunctionTest(uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunctionTest) Test(org.junit.jupiter.api.Test)

Example 35 with DoubleDoubleBiPredicate

use of uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate in project GDSC-SMLM by aherbert.

the class PrecisionTest method functionsComputeSameValue.

private static void functionsComputeSameValue(int maxx, SinglePrecision f1, DoublePrecision f2, final double precision) {
    f1.setMaxX(maxx);
    f2.setMaxX(maxx);
    final float[] p1 = params1.clone();
    final double[] p2 = params2.clone();
    p1[Gaussian.X_POSITION] = (float) (p2[Gaussian.X_POSITION] = (float) (0.123 + maxx / 2));
    p1[Gaussian.Y_POSITION] = (float) (p2[Gaussian.Y_POSITION] = (float) (0.789 + maxx / 2));
    f1.initialise(p1);
    f2.initialise(p2);
    final int n = p1.length;
    final float[] g1 = new float[n];
    final double[] g2 = new double[n];
    double t1 = 0;
    double t2 = 0;
    final double[] tg1 = new double[n];
    final double[] tg2 = new double[n];
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(precision, 0);
    for (int i = 0, limit = maxx * maxx; i < limit; i++) {
        final float v1 = f1.eval(i);
        t1 += v1;
        final double v2 = f2.eval(i);
        t2 += v2;
        TestAssertions.assertTest(v2, v1, predicate, "Different values");
        final float vv1 = f1.eval(i, g1);
        final double vv2 = f2.eval(i, g2);
        Assertions.assertEquals(v1, vv1, "Different f1 values");
        Assertions.assertEquals(v2, vv2, "Different f2 values");
        for (int j = 0; j < n; j++) {
            tg1[j] += g1[j];
            tg2[j] += g2[j];
        }
        TestAssertions.assertArrayTest(g2, toDouble(g1), predicate, "Different gradients");
    }
    TestAssertions.assertArrayTest(tg2, tg1, predicate, "Different total gradients");
    TestAssertions.assertTest(t2, t1, predicate, "Different totals");
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)

Aggregations

DoubleDoubleBiPredicate (uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)63 SeededTest (uk.ac.sussex.gdsc.test.junit5.SeededTest)27 Test (org.junit.jupiter.api.Test)22 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)12 PoissonDistribution (org.apache.commons.math3.distribution.PoissonDistribution)6 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)4 ArrayList (java.util.ArrayList)4 MixtureMultivariateGaussianDistribution (uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution)4 MultivariateGaussianDistribution (uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution.MultivariateGaussianDistribution)4 UnivariateFunction (org.apache.commons.math3.analysis.UnivariateFunction)3 SimpsonIntegrator (org.apache.commons.math3.analysis.integration.SimpsonIntegrator)3 ArrayRealVector (org.apache.commons.math3.linear.ArrayRealVector)3 RealMatrix (org.apache.commons.math3.linear.RealMatrix)3 RealVector (org.apache.commons.math3.linear.RealVector)3 ParameterizedTest (org.junit.jupiter.params.ParameterizedTest)3 BigDecimal (java.math.BigDecimal)2 MathContext (java.math.MathContext)2 MixtureMultivariateNormalDistribution (org.apache.commons.math3.distribution.MixtureMultivariateNormalDistribution)2 MultivariateNormalDistribution (org.apache.commons.math3.distribution.MultivariateNormalDistribution)2 MultivariateJacobianFunction (org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction)2