Examples with SharedStateContinuousSampler org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler used on opensource projects

Example 6 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project GDSC-SMLM by aherbert.

the class WeightedFilterTest method filterDoesNotAlterFilteredImageMean.

@SeededTest
void filterDoesNotAlterFilteredImageMean(RandomSeed seed) {
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    // ExponentialDistribution ed = new ExponentialDistribution(rand, 57,
    // ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    final DataFilter filter = createDataFilter();
    final float[] offsets = getOffsets(filter);
    final int[] boxSizes = getBoxSizes(filter);
    final TDoubleArrayList l1 = new TDoubleArrayList();
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, 2, 0.2);
    for (final int width : primes) {
        for (final int height : primes) {
            final float[] data = createData(width, height, rg);
            l1.reset();
            filter.setWeights(null, width, height);
            for (final int boxSize : boxSizes) {
                for (final float offset : offsets) {
                    for (final boolean internal : checkInternal) {
                        l1.add(getMean(data, width, height, boxSize - offset, internal, filter));
                    }
                }
            }
            final double[] e = l1.toArray();
            int ei = 0;
            // Uniform weights
            final float[] w = new float[width * height];
            Arrays.fill(w, 0.5f);
            filter.setWeights(w, width, height);
            for (final int boxSize : boxSizes) {
                for (final float offset : offsets) {
                    for (final boolean internal : checkInternal) {
                        testMean(data, width, height, boxSize - offset, internal, filter, "w=0.5", e[ei++], 1e-5);
                    }
                }
            }
            // Weights simulating the variance of sCMOS pixels
            for (int i = 0; i < w.length; i++) {
                // w[i] = (float) (1.0 / Math.max(0.01, ed.sample()));
                w[i] = (float) (1.0 / Math.max(0.01, gs.sample()));
            // w[i] = 0.5f;
            }
            ei = 0;
            filter.setWeights(w, width, height);
            for (final int boxSize : boxSizes) {
                for (final float offset : offsets) {
                    for (final boolean internal : checkInternal) {
                        testMean(data, width, height, boxSize - offset, internal, filter, "w=?", e[ei++], 5e-2);
                    }
                }
            }
        }
    }
}

Example 7 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler 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 8 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler 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 9 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler 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 10 with SharedStateContinuousSampler

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