Search in sources :

Example 21 with DoubleDoubleBiPredicate

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

the class ScmosLikelihoodWrapperTest method canComputePValue.

private static void canComputePValue(RandomSeed seed, BaseNonLinearFunction nlf) {
    logger.log(TestLogUtils.getRecord(Level.INFO, nlf.name));
    final int n = maxx * maxx;
    final double[] a = new double[] { 1 };
    // Simulate sCMOS camera
    nlf.initialise(a);
    final SCcmosLikelihoodWrapperTestData testData = (SCcmosLikelihoodWrapperTestData) dataCache.computeIfAbsent(seed, ScmosLikelihoodWrapperTest::createData);
    final float[] var = testData.var;
    final float[] g = testData.gain;
    final float[] o = testData.offset;
    final float[] sd = testData.sd;
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 1);
    final double[] k = SimpleArrayUtils.newArray(n, 0, 1.0);
    for (int i = 0; i < n; i++) {
        double mean = nlf.eval(i);
        if (mean > 0) {
            mean = GdscSmlmTestUtils.createPoissonSampler(r, mean).sample();
        }
        k[i] = mean * g[i] + o[i] + gs.sample() * sd[i];
    }
    final ScmosLikelihoodWrapper f = new ScmosLikelihoodWrapper(nlf, a, k, n, var, g, o);
    final double oll = f.computeObservedLikelihood();
    double oll2 = 0;
    final double[] op = new double[n];
    for (int j = 0; j < n; j++) {
        op[j] = ScmosLikelihoodWrapper.likelihood((k[j] - o[j]) / g[j], var[j], g[j], o[j], k[j]);
        oll2 -= Math.log(op[j]);
    }
    logger.log(TestLogUtils.getRecord(Level.INFO, "oll=%f, oll2=%f", oll, oll2));
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    TestAssertions.assertTest(oll2, oll, predicate, "Observed Log-likelihood");
    final TDoubleArrayList list = new TDoubleArrayList();
    final int imin = 5;
    final int imax = 15;
    for (int i = imin; i <= imax; i++) {
        a[0] = (double) i / 10;
        final double ll = f.likelihood(a);
        list.add(ll);
        final double llr = f.computeLogLikelihoodRatio(ll);
        BigDecimal product = new BigDecimal(1);
        double ll2 = 0;
        for (int j = 0; j < n; j++) {
            final double p1 = ScmosLikelihoodWrapper.likelihood(nlf.eval(j), var[j], g[j], o[j], k[j]);
            ll2 -= Math.log(p1);
            final double ratio = p1 / op[j];
            product = product.multiply(new BigDecimal(ratio));
        }
        final double llr2 = -2 * Math.log(product.doubleValue());
        final double q = f.computeQValue(ll);
        logger.log(TestLogUtils.getRecord(Level.INFO, "a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f", a[0], ll, ll2, llr, llr2, product.round(new MathContext(4)).toString(), q));
        // too small to store in a double.
        if (product.doubleValue() > 0) {
            TestAssertions.assertTest(llr, llr2, predicate, "Log-likelihood");
        }
    }
    // Find min using quadratic fit
    final double[] data = list.toArray();
    int index = SimpleArrayUtils.findMinIndex(data);
    final double mina = (double) (imin + index) / 10;
    double fita = mina;
    try {
        if (index == 0) {
            index++;
        }
        if (index == data.length - 1) {
            index--;
        }
        final int i1 = index - 1;
        final int i2 = index;
        final int i3 = index + 1;
        fita = QuadraticUtils.findMinMax((double) (imin + i1) / 10, data[i1], (double) (imin + i2) / 10, data[i2], (double) (imin + i3) / 10, data[i3]);
    } catch (final DataException ex) {
    // Ignore
    }
    // Allow a tolerance as the random data may alter the p-value computation.
    // Should allow it to be less than 2 increment either side of the answer.
    logger.log(TestLogUtils.getRecord(Level.INFO, "min fit = %g => %g", mina, fita));
    Assertions.assertEquals(1, fita, 0.199, "min");
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) BigDecimal(java.math.BigDecimal) MathContext(java.math.MathContext) DataException(uk.ac.sussex.gdsc.core.data.DataException) TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider)

Example 22 with DoubleDoubleBiPredicate

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

the class PoissonGaussianConvolutionFunctionTest method probabilityMatchesLogProbability.

private static void probabilityMatchesLogProbability(final double gain, double mu, final double sd, boolean computePmf) {
    // Note: The input s parameter is pre-gain.
    final PoissonGaussianConvolutionFunction f = PoissonGaussianConvolutionFunction.createWithStandardDeviation(1.0 / gain, sd * gain);
    f.setComputePmf(computePmf);
    // 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, erf=%b", gain, mu, sd, computePmf);
    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 23 with DoubleDoubleBiPredicate

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

the class FactorialTest method testFactorialDouble.

/**
 * Test the factorial of a fractional number against Commons Math gamma(1+n).
 */
@SeededTest
void testFactorialDouble(RandomSeed seed) {
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final DoubleDoubleBiPredicate tol = TestHelper.doublesAreClose(5e-15).or(TestHelper.doublesEqual());
    for (int i = 0; i < 100; i++) {
        final double n = rng.nextDouble() * 180;
        final double expected = n < 1.5 ? 1 / (1 + Gamma.invGamma1pm1(n)) : Gamma.gamma(1 + n);
        TestAssertions.assertTest(expected, Factorial.value(n), tol, () -> Double.toString(n));
    }
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 24 with DoubleDoubleBiPredicate

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

the class FactorialTest method testFactorial.

/**
 * Test the factorial of a fractional number against reference values.
 */
@ParameterizedTest
@CsvSource({ // @formatter:off
"0.25, 0.906402477055477077939", "0.75, 0.919062526848883233914", "1.125, 1.05946053733091417382", "1.75, 1.60835942198554565977", "2.75, 4.42298841046025056293", "3.25, 8.28508514183522016498", "4.75, 78.7844810613232131788", "6.25, 1155.3810139199896867", "12.75, 3255990905.16494153091", "17.125, 508919418354999.991547", "27.75, 1.32099626069721824877e+29", "57.75, 8.50381139447823243123e+77", "97.75, 2.99327649630298942593e+153", "137.75, 2.01698432357024766181e+236", "154.25, 1.08954758738660925068e+272", "164.25, 1.17747769003558155385e+294", "170.25, 2.62296687867878940345e+307", "170.624, 1.79421175992481035917e+308", // Limit
"170.6243769563027, 1.7976931348622301e+308", // Infinite
"170.62437695630274, 1.79769313486249261288e+308", "170.625, 1.80346206550076047216e+308", "171.5, 1.62639753771045308624e+310" // @formatter:on
})
void testFactorial(double n, double expected) {
    final DoubleDoubleBiPredicate tol = TestHelper.doublesAreClose(1e-15).or(TestHelper.doublesEqual());
    TestAssertions.assertTest(expected, Factorial.value(n), tol, () -> Double.toString(n));
}
Also used : DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) CsvSource(org.junit.jupiter.params.provider.CsvSource) ParameterizedTest(org.junit.jupiter.params.ParameterizedTest)

Example 25 with DoubleDoubleBiPredicate

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

the class MultivariateGaussianMixtureExpectationMaximizationTest method canCreateMixedMultivariateGaussianDistribution.

@SeededTest
void canCreateMixedMultivariateGaussianDistribution(RandomSeed seed) {
    // Will be normalised
    final double[] weights = { 1, 1 };
    final double[][] means = new double[2][];
    final double[][][] covariances = new double[2][][];
    final double[][] data1 = { { 1, 2 }, { 2.5, 1.5 }, { 3.5, 1.0 } };
    final double[][] data2 = { { 4, 2 }, { 3.5, -1.5 }, { -3.5, 1.0 } };
    means[0] = getColumnMeans(data1);
    covariances[0] = getCovariance(data1);
    means[1] = getColumnMeans(data2);
    covariances[1] = getCovariance(data2);
    // Create components. This does not have to be zero based.
    final LocalList<double[]> list = new LocalList<>();
    list.addAll(Arrays.asList(data1));
    list.addAll(Arrays.asList(data2));
    final double[][] data = list.toArray(new double[0][]);
    final int[] components = { -1, -1, -1, 3, 3, 3 };
    // Randomise the data
    for (int n = 0; n < 3; n++) {
        final long start = n + seed.getSeedAsLong();
        // This relies on the shuffle being the same
        RandomUtils.shuffle(data, RngUtils.create(start));
        RandomUtils.shuffle(components, RngUtils.create(start));
        final MixtureMultivariateGaussianDistribution dist = MultivariateGaussianMixtureExpectationMaximization.createMixed(data, components);
        Assertions.assertArrayEquals(new double[] { 0.5, 0.5 }, dist.getWeights());
        final MultivariateGaussianDistribution[] distributions = dist.getDistributions();
        Assertions.assertEquals(weights.length, distributions.length);
        final DoubleDoubleBiPredicate test = TestHelper.doublesAreClose(1e-8);
        for (int i = 0; i < means.length; i++) {
            TestAssertions.assertArrayTest(means[i], distributions[i].getMeans(), test);
            TestAssertions.assertArrayTest(covariances[i], distributions[i].getCovariances(), test);
        }
    }
}
Also used : LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate) MixtureMultivariateGaussianDistribution(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution) MultivariateGaussianDistribution(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution.MultivariateGaussianDistribution) MixtureMultivariateGaussianDistribution(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

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