Examples with UniformRandomProvider org.apache.commons.rng.UniformRandomProvider used on opensource projects

Example 56 with UniformRandomProvider

use of org.apache.commons.rng.UniformRandomProvider in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method canComputeLogLikelihoodRatio.

private static void canComputeLogLikelihoodRatio(RandomSeed seed, BaseNonLinearFunction nlf) {
    logger.log(TestLogUtils.getRecord(LOG_LEVEL, nlf.name));
    final int n = maxx * maxx;
    final double[] a = new double[] { 1 };
    // Simulate Poisson process
    nlf.initialise(a);
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final double[] x = new double[n];
    final double[] u = new double[n];
    for (int i = 0; i < n; i++) {
        u[i] = nlf.eval(i);
        if (u[i] > 0) {
            x[i] = GdscSmlmTestUtils.createPoissonSampler(rng, u[i]).sample();
        }
    }
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    double ll = PoissonCalculator.logLikelihood(u, x);
    final double mll = PoissonCalculator.maximumLogLikelihood(x);
    double llr = -2 * (ll - mll);
    double llr2 = PoissonCalculator.logLikelihoodRatio(u, x);
    logger.log(TestLogUtils.getRecord(LOG_LEVEL, "llr=%f, llr2=%f", llr, llr2));
    TestAssertions.assertTest(llr, llr2, predicate, "Log-likelihood ratio");
    final double[] op = new double[x.length];
    for (int i = 0; i < n; i++) {
        op[i] = PoissonCalculator.maximumLikelihood(x[i]);
    }
    // TestSettings.setLogLevel(uk.ac.sussex.gdsc.smlm.TestLog.Level.FINE);
    final int df = n - 1;
    final ChiSquaredDistributionTable table = ChiSquaredDistributionTable.createUpperTailed(0.05, df);
    final ChiSquaredDistributionTable table2 = ChiSquaredDistributionTable.createUpperTailed(0.001, df);
    if (logger.isLoggable(LOG_LEVEL)) {
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "Chi2 = %f (q=%.3f), %f (q=%.3f)  %f %b  %f", table.getCrititalValue(df), table.getSignificanceValue(), table2.getCrititalValue(df), table2.getSignificanceValue(), ChiSquaredDistributionTable.computeQValue(24, 2), ChiSquaredDistributionTable.createUpperTailed(0.05, 2).reject(24, 2), ChiSquaredDistributionTable.getChiSquared(1e-6, 2)));
    }
    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;
        nlf.initialise(a);
        for (int j = 0; j < n; j++) {
            u[j] = nlf.eval(j);
        }
        ll = PoissonCalculator.logLikelihood(u, x);
        list.add(ll);
        llr = PoissonCalculator.logLikelihoodRatio(u, x);
        BigDecimal product = new BigDecimal(1);
        double ll2 = 0;
        for (int j = 0; j < n; j++) {
            final double p1 = PoissonCalculator.likelihood(u[j], x[j]);
            ll2 += Math.log(p1);
            final double ratio = p1 / op[j];
            product = product.multiply(new BigDecimal(ratio));
        }
        llr2 = -2 * Math.log(product.doubleValue());
        final double p = ChiSquaredDistributionTable.computePValue(llr, df);
        final double q = ChiSquaredDistributionTable.computeQValue(llr, df);
        logger.log(TestLogUtils.getRecord(LOG_LEVEL, "a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f, q=%f " + "(reject=%b @ %.3f, reject=%b @ %.3f)", a[0], ll, ll2, llr, llr2, product.round(new MathContext(4)).toString(), p, q, table.reject(llr, df), table.getSignificanceValue(), table2.reject(llr, df), table2.getSignificanceValue()));
        // too small to store in a double.
        if (product.doubleValue() > 0) {
            TestAssertions.assertTest(ll, ll2, predicate, "Log-likelihood");
            TestAssertions.assertTest(llr, llr2, predicate, "Log-likelihood ratio");
        }
    }
    // Find max using quadratic fit
    final double[] data = list.toArray();
    int index = SimpleArrayUtils.findMaxIndex(data);
    final double maxa = (double) (imin + index) / 10;
    double fita = maxa;
    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(LOG_LEVEL, "max fit = %g => %g", maxa, fita));
    Assertions.assertEquals(1, fita, 0.199, "max");
}

Example 57 with UniformRandomProvider

use of org.apache.commons.rng.UniformRandomProvider in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method cannotSubtractConstantBackgroundAndComputeLogLikelihoodRatio.

private static void cannotSubtractConstantBackgroundAndComputeLogLikelihoodRatio(RandomSeed seed, BaseNonLinearFunction nlf1, BaseNonLinearFunction nlf2, BaseNonLinearFunction nlf3) {
    // System.out.println(nlf1.name);
    final int n = maxx * maxx;
    final double[] a = new double[] { 1 };
    // Simulate Poisson process of 3 combined functions
    nlf1.initialise(a);
    nlf2.initialise(a);
    nlf3.initialise(a);
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    double[] x = SimpleArrayUtils.newArray(n, 0, 1.0);
    final double[] u = new double[x.length];
    final double[] b1 = new double[x.length];
    final double[] b2 = new double[x.length];
    final double[] b3 = new double[x.length];
    for (int i = 0; i < n; i++) {
        b1[i] = nlf1.eval(i);
        b2[i] = nlf2.eval(i);
        b3[i] = nlf3.eval(i);
        u[i] = b1[i] + b2[i] + b3[i];
        if (u[i] > 0) {
            x[i] = GdscSmlmTestUtils.createPoissonSampler(rng, u[i]).sample();
        }
    }
    // x is the target data
    // b1 is the already computed background
    // b2 is the first function to add to the model
    // b3 is the second function to add to the model
    // Adjust x to ensure that x - background is positive.
    double[] xb = subtract(x, b1);
    SimpleArrayUtils.apply(xb, z -> z < 0 ? 0 : z);
    x = add(xb, b1);
    // Compute the LLR of adding b3 to b2 given we already have b1 modelling data x
    final double[] b12 = add(b1, b2);
    final double ll1a = PoissonCalculator.logLikelihood(b12, x);
    final double ll2a = PoissonCalculator.logLikelihood(add(b12, b3), x);
    final double llra = -2 * (ll1a - ll2a);
    // logger.fine(FunctionUtils.getSupplier("x|(a+b+c) ll1=%f, ll2=%f, llra=%f", ll1a, ll2a, llra);
    // Compute the LLR of adding b3 to b2 given we already have x minus b1
    final double ll1b = PoissonCalculator.logLikelihood(b2, xb);
    final double ll2b = PoissonCalculator.logLikelihood(add(b2, b3), xb);
    final double llrb = -2 * (ll1b - ll2b);
    // logger.fine(FunctionUtils.getSupplier("x-a|(b+c) : ll1=%f, ll2=%f, llrb=%f", ll1b, ll2b,
    // llrb);
    // logger.fine(FunctionUtils.getSupplier("llr=%f (%g), llr2=%f (%g)", llra,
    // PoissonCalculator.computePValue(llra, 1), llrb,
    // PoissonCalculator.computePValue(llrb, 1));
    Assertions.assertThrows(AssertionError.class, () -> {
        TestAssertions.assertTest(llra, llrb, TestHelper.doublesAreClose(1e-10, 0), "Log-likelihood ratio");
    });
}

Example 58 with UniformRandomProvider

use of org.apache.commons.rng.UniformRandomProvider in project GDSC-SMLM by aherbert.

the class BoundedLvmSteppingFunctionSolverTest method fitSingleGaussianWithoutBias.

private void fitSingleGaussianWithoutBias(RandomSeed seed, boolean applyBounds, int clamping) {
    Assumptions.assumeTrue(runTests);
    final double bias = 100;
    final SteppingFunctionSolver solver = getSolver(clamping, false);
    final SteppingFunctionSolver solver2 = getSolver(clamping, false);
    final String name = getLvmName(applyBounds, clamping, false);
    final int loops = 5;
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    final StoredDataStatistics[] stats = new StoredDataStatistics[6];
    for (final double s : signal) {
        final double[] expected = createParams(1, s, 0, 0, 1);
        double[] lower = null;
        double[] upper = null;
        if (applyBounds) {
            lower = createParams(0, s * 0.5, -0.2, -0.2, 0.8);
            upper = createParams(3, s * 2, 0.2, 0.2, 1.2);
            solver.setBounds(lower, upper);
        }
        final double[] expected2 = addBiasToParams(expected, bias);
        if (applyBounds) {
            final double[] lower2 = addBiasToParams(lower, bias);
            final double[] upper2 = addBiasToParams(upper, bias);
            solver2.setBounds(lower2, upper2);
        }
        for (int loop = loops; loop-- > 0; ) {
            final double[] data = drawGaussian(expected, rg);
            final double[] data2 = data.clone();
            for (int i = 0; i < data.length; i++) {
                data2[i] += bias;
            }
            for (int i = 0; i < stats.length; i++) {
                stats[i] = new StoredDataStatistics();
            }
            for (final double db : base) {
                for (final double dx : shift) {
                    for (final double dy : shift) {
                        for (final double dsx : factor) {
                            final double[] p = createParams(db, s, dx, dy, dsx);
                            final double[] p2 = addBiasToParams(p, bias);
                            final double[] fp = fitGaussian(solver, data, p, expected);
                            final double[] fp2 = fitGaussian(solver2, data2, p2, expected2);
                            // The result should be the same without a bias
                            Assertions.assertEquals(solver.getEvaluations(), solver2.getEvaluations(), () -> name + " Iterations");
                            fp2[0] -= bias;
                            Assertions.assertArrayEquals(fp, fp2, 1e-6, () -> name + " Solution");
                        }
                    }
                }
            }
        }
    }
}

Example 59 with UniformRandomProvider

use of org.apache.commons.rng.UniformRandomProvider 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);
    }
}

Example 60 with UniformRandomProvider

use of org.apache.commons.rng.UniformRandomProvider in project GDSC-SMLM by aherbert.

the class OffsetFunctionTest method offsetValueFunctionWrapsPrecomputedValues.

@SeededTest
void offsetValueFunctionWrapsPrecomputedValues(RandomSeed seed) {
    final int n = 3;
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    final ValueFunction f0 = new FakeGradientFunction(3, n);
    final int size = f0.size();
    final double[] b1 = GdscSmlmTestUtils.generateDoubles(size, r);
    final double[] b2 = GdscSmlmTestUtils.generateDoubles(size, r);
    final ValueFunction f1 = OffsetValueFunction.wrapValueFunction(f0, b1);
    final ValueFunction f2 = OffsetValueFunction.wrapValueFunction(f1, b2);
    final double[] p = new double[n];
    for (int i = 0; i < n; i++) {
        p[i] = r.nextDouble();
    }
    final double[] v0 = evaluateValueFunction(f0, p);
    final double[] v1 = evaluateValueFunction(f1, p);
    final double[] v2 = evaluateValueFunction(f2, p);
    for (int i = 0; i < v0.length; i++) {
        final double e = v0[i] + b1[i] + b2[i];
        final double o1 = v1[i] + b2[i];
        final double o2 = v2[i];
        Assertions.assertEquals(e, o1, "o1");
        Assertions.assertEquals(e, o2, 1e-6, "o2");
    }
}