Search in sources :

Example 41 with RandomSeed

use of uk.ac.sussex.gdsc.test.junit5.RandomSeed in project GDSC-SMLM by aherbert.

the class PoissonGradientProcedureTest method crlbIsHigherWithPrecomputed.

@SeededTest
void crlbIsHigherWithPrecomputed(RandomSeed seed) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final int n = func.getNumberOfGradients();
    // Get a background
    final double[] b = new double[func.size()];
    for (int i = 0; i < b.length; i++) {
        b[i] = nextUniform(rng, 1, 2);
    }
    for (int i = 0; i < iter; i++) {
        a[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
        a[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
        a[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
        a[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a[Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        final PoissonGradientProcedure p1 = PoissonGradientProcedureUtils.create(func);
        p1.computeFisherInformation(a);
        final PoissonGradientProcedure p2 = PoissonGradientProcedureUtils.create(OffsetGradient1Function.wrapGradient1Function(func, b));
        p2.computeFisherInformation(a);
        final FisherInformationMatrix m1 = new FisherInformationMatrix(p1.getLinear(), n);
        final FisherInformationMatrix m2 = new FisherInformationMatrix(p2.getLinear(), n);
        final double[] crlb1 = m1.crlb();
        final double[] crlb2 = m2.crlb();
        Assertions.assertNotNull(crlb1);
        Assertions.assertNotNull(crlb2);
        // Arrays.toString(crlb2));
        for (int j = 0; j < n; j++) {
            Assertions.assertTrue(crlb1[j] < crlb2[j]);
        }
    }
}
Also used : ErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction) FisherInformationMatrix(uk.ac.sussex.gdsc.smlm.fitting.FisherInformationMatrix) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 42 with RandomSeed

use of uk.ac.sussex.gdsc.test.junit5.RandomSeed in project GDSC-SMLM by aherbert.

the class LogFactorialTest method testLogFactorialDouble.

/**
 * Test the factorial of a fractional number against Commons Math logGamma(1+n).
 */
@SeededTest
void testLogFactorialDouble(RandomSeed seed) {
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final DoubleDoubleBiPredicate tol = TestHelper.doublesAreClose(1e-14);
    for (int i = 0; i < 200; i++) {
        final double n = rng.nextDouble() * 200;
        final double expected = n <= 1.5 ? Gamma.logGamma1p(n) : Gamma.logGamma(1 + n);
        TestAssertions.assertTest(expected, LogFactorial.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 43 with RandomSeed

use of uk.ac.sussex.gdsc.test.junit5.RandomSeed in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method instanceAndFastMethodIsApproximatelyEqualToStaticMethod.

@SeededTest
void instanceAndFastMethodIsApproximatelyEqualToStaticMethod(RandomSeed seed) {
    final DoubleEquality eq = new DoubleEquality(3e-4, 0);
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    // Test for different x. The calculator approximation begins
    final int n = 100;
    final double[] u = new double[n];
    final double[] x = new double[n];
    double expected;
    double observed;
    for (final double testx : new double[] { Math.nextDown(PoissonCalculator.APPROXIMATION_X), PoissonCalculator.APPROXIMATION_X, Math.nextUp(PoissonCalculator.APPROXIMATION_X), PoissonCalculator.APPROXIMATION_X * 1.1, PoissonCalculator.APPROXIMATION_X * 2, PoissonCalculator.APPROXIMATION_X * 10 }) {
        final String X = Double.toString(testx);
        Arrays.fill(x, testx);
        final PoissonCalculator pc = new PoissonCalculator(x);
        expected = PoissonCalculator.maximumLogLikelihood(x);
        observed = pc.getMaximumLogLikelihood();
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "[%s] Instance MaxLL = %g vs %g (error = %g)", X, expected, observed, DoubleEquality.relativeError(expected, observed)));
        Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(expected, observed), () -> "Instance Max LL not equal: x=" + X);
        observed = PoissonCalculator.fastMaximumLogLikelihood(x);
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "[%s] Fast MaxLL = %g vs %g (error = %g)", X, expected, observed, DoubleEquality.relativeError(expected, observed)));
        Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(expected, observed), () -> "Fast Max LL not equal: x=" + X);
        // Generate data around the value
        for (int i = 0; i < n; i++) {
            u[i] = x[i] + rg.nextDouble() - 0.5;
        }
        expected = PoissonCalculator.logLikelihood(u, x);
        observed = pc.logLikelihood(u);
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "[%s] Instance LL = %g vs %g (error = %g)", X, expected, observed, DoubleEquality.relativeError(expected, observed)));
        Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(expected, observed), () -> "Instance LL not equal: x=" + X);
        observed = PoissonCalculator.fastLogLikelihood(u, x);
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "[%s] Fast LL = %g vs %g (error = %g)", X, expected, observed, DoubleEquality.relativeError(expected, observed)));
        Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(expected, observed), () -> "Fast LL not equal: x=" + X);
        expected = PoissonCalculator.logLikelihoodRatio(u, x);
        observed = pc.getLogLikelihoodRatio(observed);
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "[%s] Instance LLR = %g vs %g (error = %g)", X, expected, observed, DoubleEquality.relativeError(expected, observed)));
        Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(expected, observed), () -> "Instance LLR not equal: x=" + X);
    }
}
Also used : DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 44 with RandomSeed

use of uk.ac.sussex.gdsc.test.junit5.RandomSeed in project GDSC-SMLM by aherbert.

the class FastMleGradient2ProcedureTest method gradientCalculatorComputesGradient.

@SeededTest
void gradientCalculatorComputesGradient(RandomSeed seed) {
    gradientCalculatorComputesGradient(seed, new SingleFreeCircularErfGaussian2DFunction(blockWidth, blockWidth));
    // Use a reasonable z-depth function from the Smith, et al (2010) paper (page 377)
    final double sx = 1.08;
    final double sy = 1.01;
    final double gamma = 0.389;
    final double d = 0.531;
    final double Ax = -0.0708;
    final double Bx = -0.073;
    final double Ay = 0.164;
    final double By = 0.0417;
    final HoltzerAstigmatismZModel zModel = HoltzerAstigmatismZModel.create(sx, sy, gamma, d, Ax, Bx, Ay, By);
    gradientCalculatorComputesGradient(seed, new SingleAstigmatismErfGaussian2DFunction(blockWidth, blockWidth, zModel));
}
Also used : SingleAstigmatismErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.SingleAstigmatismErfGaussian2DFunction) HoltzerAstigmatismZModel(uk.ac.sussex.gdsc.smlm.function.gaussian.HoltzerAstigmatismZModel) SingleFreeCircularErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 45 with RandomSeed

use of uk.ac.sussex.gdsc.test.junit5.RandomSeed in project GDSC-SMLM by aherbert.

the class FastMleGradient2ProcedureTest method gradientProcedureComputesSameWithPrecomputed.

@SeededTest
void gradientProcedureComputesSameWithPrecomputed(RandomSeed seed) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a1 = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] a2 = new double[1 + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] x = new double[f1.size()];
    final double[] b = new double[f1.size()];
    for (int i = 0; i < iter; i++) {
        final int ii = i;
        a2[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
        a2[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a2[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 3, 5);
        a2[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 3, 5);
        a2[Gaussian2DFunction.Z_POSITION] = nextUniform(rng, -2, 2);
        a2[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_POSITION] = nextUniform(rng, 5, 7);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 5, 7);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Z_POSITION] = nextUniform(rng, -3, 1);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        // Simulate Poisson data
        f2.initialise0(a2);
        f1.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                if (value > 0) {
                    x[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
                } else {
                    x[index++] = 0;
                }
            }
        });
        // Precompute peak 2 (no background)
        a1[Gaussian2DFunction.BACKGROUND] = 0;
        for (int j = 1; j < 7; j++) {
            a1[j] = a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + j];
        }
        f1.initialise0(a1);
        f1.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                b[index++] = value;
            }
        });
        // Reset to peak 1
        for (int j = 0; j < 7; j++) {
            a1[j] = a2[j];
        }
        // Compute peak 1+2
        final FastMleGradient2Procedure p12 = FastMleGradient2ProcedureUtils.create(x, f2);
        p12.computeSecondDerivative(a2);
        final double[] d11 = Arrays.copyOf(p12.d1, f1.getNumberOfGradients());
        final double[] d21 = Arrays.copyOf(p12.d2, f1.getNumberOfGradients());
        // Compute peak 1+(precomputed 2)
        final FastMleGradient2Procedure p1b2 = FastMleGradient2ProcedureUtils.create(x, OffsetGradient2Function.wrapGradient2Function(f1, b));
        p1b2.computeSecondDerivative(a1);
        final double[] d12 = p1b2.d1;
        final double[] d22 = p1b2.d2;
        Assertions.assertArrayEquals(p12.u, p1b2.u, 1e-10, () -> " Result: Not same @ " + ii);
        Assertions.assertArrayEquals(d11, d12, 1e-10, () -> " D1: Not same @ " + ii);
        Assertions.assertArrayEquals(d21, d22, 1e-10, () -> " D2: Not same @ " + ii);
        final double[] v1 = p12.computeValue(a2);
        final double[] v2 = p1b2.computeValue(a1);
        Assertions.assertArrayEquals(v1, v2, 1e-10, () -> " Value: Not same @ " + ii);
        final double[] d1 = Arrays.copyOf(p12.computeFirstDerivative(a2), f1.getNumberOfGradients());
        final double[] d2 = p1b2.computeFirstDerivative(a1);
        Assertions.assertArrayEquals(d1, d2, 1e-10, () -> " 1st derivative: Not same @ " + ii);
    }
}
Also used : ValueProcedure(uk.ac.sussex.gdsc.smlm.function.ValueProcedure) SingleFreeCircularErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction) SingleAstigmatismErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.SingleAstigmatismErfGaussian2DFunction) ErfGaussian2DFunction(uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Aggregations

SeededTest (uk.ac.sussex.gdsc.test.junit5.SeededTest)161 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)142 DoubleDoubleBiPredicate (uk.ac.sussex.gdsc.test.api.function.DoubleDoubleBiPredicate)17 Rectangle (java.awt.Rectangle)12 SpeedTag (uk.ac.sussex.gdsc.test.junit5.SpeedTag)12 TimingService (uk.ac.sussex.gdsc.test.utils.TimingService)11 SharedStateContinuousSampler (org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler)7 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)6 FloatProcessor (ij.process.FloatProcessor)5 SingleFreeCircularErfGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction)5 TimingResult (uk.ac.sussex.gdsc.test.utils.TimingResult)5 ArrayList (java.util.ArrayList)4 ErfGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction)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 DoubleConvolutionValueProcedure (uk.ac.sussex.gdsc.smlm.utils.Convolution.DoubleConvolutionValueProcedure)4 FloatFloatBiPredicate (uk.ac.sussex.gdsc.test.api.function.FloatFloatBiPredicate)4 MixtureMultivariateNormalDistribution (org.apache.commons.math3.distribution.MixtureMultivariateNormalDistribution)3 MultivariateNormalMixtureExpectationMaximization (org.apache.commons.math3.distribution.fitting.MultivariateNormalMixtureExpectationMaximization)3 DoubleEquality (uk.ac.sussex.gdsc.core.utils.DoubleEquality)3