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

Example 61 with UniformRandomProvider

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

the class OffsetFunctionTest method offsetGradient2FunctionWrapsPrecomputedValues.

@SeededTest
void offsetGradient2FunctionWrapsPrecomputedValues(RandomSeed seed) {
    final int n = 3;
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    final Gradient2Function 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 Gradient2Function f1 = OffsetGradient2Function.wrapGradient2Function(f0, b1);
    final Gradient2Function f2 = OffsetGradient2Function.wrapGradient2Function(f1, b2);
    final double[] p = new double[n];
    for (int i = 0; i < n; i++) {
        p[i] = r.nextDouble();
    }
    final double[] d0 = new double[n];
    final double[] d1 = new double[n];
    final double[] d2 = new double[n];
    final double[] d20 = new double[n];
    final double[] d21 = new double[n];
    final double[] d22 = new double[n];
    final double[] v0 = evaluateGradient2Function(f0, p, d0, d20);
    final double[] v1 = evaluateGradient2Function(f1, p, d1, d21);
    final double[] v2 = evaluateGradient2Function(f2, p, d2, d22);
    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");
    }
    Assertions.assertArrayEquals(d0, d1, "d1");
    Assertions.assertArrayEquals(d0, d2, "d2");
    Assertions.assertArrayEquals(d20, d21, "d21");
    Assertions.assertArrayEquals(d20, d22, "d22");
}

Example 62 with UniformRandomProvider

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

the class ScmosLikelihoodWrapperTest method functionComputesTargetGradient.

private void functionComputesTargetGradient(RandomSeed seed, Gaussian2DFunction f1, int targetParameter, double threshold) {
    final int[] indices = f1.gradientIndices();
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final double[] dyda = new double[indices.length];
    double[] params;
    ScmosLikelihoodWrapper ff1;
    final int n = maxx * maxx;
    int count = 0;
    int total = 0;
    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);
    for (final double background : testbackground) {
        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);
                                // Create y as a function we would want to move towards
                                final double[] a2 = params.clone();
                                a2[targetParameter] *= 1.3;
                                f1.initialise(a2);
                                final double[] data = new double[n];
                                for (int i = 0; i < n; i++) {
                                    // Simulate sCMOS camera
                                    final double u = f1.eval(i);
                                    data[i] = GdscSmlmTestUtils.createPoissonSampler(r, u).sample() * g[i] + o[i] + gs.sample() * sd[i];
                                }
                                ff1 = new ScmosLikelihoodWrapper(f1, params, data, n, var, g, o);
                                // Numerically solve gradient.
                                // Calculate the step size h to be an exact numerical representation
                                final double xx = params[targetParameter];
                                // Get h to minimise roundoff error
                                final double h = Precision.representableDelta(xx, stepH);
                                ff1.likelihood(getVariables(indices, params), dyda);
                                // Evaluate at (x+h) and (x-h)
                                params[targetParameter] = xx + h;
                                final double value2 = ff1.likelihood(getVariables(indices, params));
                                params[targetParameter] = xx - h;
                                final double value3 = ff1.likelihood(getVariables(indices, params));
                                final double gradient = (value2 - value3) / (2 * h);
                                boolean ok = Math.signum(gradient) == Math.signum(dyda[gradientIndex]) || Math.abs(gradient - dyda[gradientIndex]) < 0.1;
                                // dyda[gradientIndex]));
                                if (!ok) {
                                    Assertions.fail(NAME[targetParameter] + ": " + gradient + " != " + dyda[gradientIndex]);
                                }
                                ok = eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]);
                                if (ok) {
                                    count++;
                                }
                                total++;
                            }
                        }
                    }
                }
            }
        }
    }
    final double p = (100.0 * count) / total;
    logger.log(TestLogUtils.getRecord(Level.INFO, "%s : %s = %d / %d (%.2f)", f1.getClass().getSimpleName(), NAME[targetParameter], count, total, p));
    Assertions.assertTrue(p > threshold, FunctionUtils.getSupplier("%s fraction too low: %s", NAME[targetParameter], p));
}

Example 63 with UniformRandomProvider

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

Example 64 with UniformRandomProvider

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

the class ScmosLikelihoodWrapperTest method createData.

private static Object createData(RandomSeed source) {
    final int n = maxx * maxx;
    final SCcmosLikelihoodWrapperTestData data = new SCcmosLikelihoodWrapperTestData();
    data.var = new float[n];
    data.gain = new float[n];
    data.offset = new float[n];
    data.sd = new float[n];
    final UniformRandomProvider rg = RngUtils.create(source.getSeed());
    final DiscreteSampler pd = GdscSmlmTestUtils.createPoissonSampler(rg, O);
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, G, G_SD);
    final ContinuousSampler ed = SamplerUtils.createExponentialSampler(rg, VAR);
    for (int i = 0; i < n; i++) {
        data.offset[i] = pd.sample();
        data.var[i] = (float) ed.sample();
        data.sd[i] = (float) Math.sqrt(data.var[i]);
        data.gain[i] = (float) gs.sample();
    }
    return data;
}

Example 65 with UniformRandomProvider

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

the class ScmosLikelihoodWrapperTest method functionComputesTargetGradientPerDatum.

private void functionComputesTargetGradientPerDatum(RandomSeed seed, Gaussian2DFunction f1, int targetParameter) {
    final int[] indices = f1.gradientIndices();
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final double[] dyda = new double[indices.length];
    double[] params;
    ScmosLikelihoodWrapper ff1;
    final int n = maxx * maxx;
    int count = 0;
    int total = 0;
    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);
    for (final double background : testbackground) {
        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);
                                // Create y as a function we would want to move towards
                                final double[] a2 = params.clone();
                                a2[targetParameter] *= 1.1;
                                f1.initialise(a2);
                                final double[] data = new double[n];
                                for (int i = 0; i < n; i++) {
                                    // Simulate sCMOS camera
                                    final double u = f1.eval(i);
                                    data[i] = GdscSmlmTestUtils.createPoissonSampler(r, u).sample() * g[i] + o[i] + gs.sample() * sd[i];
                                }
                                ff1 = new ScmosLikelihoodWrapper(f1, params, data, n, var, g, o);
                                // Numerically solve gradient.
                                // Calculate the step size h to be an exact numerical representation
                                final double xx = params[targetParameter];
                                // Get h to minimise roundoff error
                                final double h = Precision.representableDelta(xx, stepH);
                                for (final int x : testx) {
                                    for (final int y : testy) {
                                        final int i = y * maxx + x;
                                        params[targetParameter] = xx;
                                        ff1.likelihood(getVariables(indices, params), dyda, i);
                                        // Evaluate at (x+h) and (x-h)
                                        params[targetParameter] = xx + h;
                                        final double value2 = ff1.likelihood(getVariables(indices, params), i);
                                        params[targetParameter] = xx - h;
                                        final double value3 = ff1.likelihood(getVariables(indices, params), i);
                                        final double gradient = (value2 - value3) / (2 * h);
                                        boolean ok = Math.signum(gradient) == Math.signum(dyda[gradientIndex]) || Math.abs(gradient - dyda[gradientIndex]) < 0.1;
                                        // dyda[gradientIndex]);
                                        if (!ok) {
                                            Assertions.fail(NAME[targetParameter] + ": " + gradient + " != " + dyda[gradientIndex]);
                                        }
                                        ok = eqPerDatum.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]);
                                        if (ok) {
                                            count++;
                                        }
                                        total++;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    final double p = (100.0 * count) / total;
    logger.log(TestLogUtils.getRecord(Level.INFO, "Per Datum %s : %s = %d / %d (%.2f)", f1.getClass().getSimpleName(), NAME[targetParameter], count, total, p));
    Assertions.assertTrue(p > 90, () -> NAME[targetParameter] + " fraction too low per datum: " + p);
}