Examples with Well19937c org.apache.commons.math3.random.Well19937c used on opensource projects

Example 16 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method canComputeLogLikelihoodRatio.

private void canComputeLogLikelihoodRatio(BaseNonLinearFunction nlf) {
    System.out.println(nlf.name);
    int n = maxx * maxx;
    double[] a = new double[] { 1 };
    // Simulate Poisson process
    nlf.initialise(a);
    RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    double[] x = new double[n];
    double[] u = new double[n];
    for (int i = 0; i < n; i++) {
        u[i] = nlf.eval(i);
        if (u[i] > 0)
            x[i] = rdg.nextPoisson(u[i]);
    }
    double ll = PoissonCalculator.logLikelihood(u, x);
    double mll = PoissonCalculator.maximumLogLikelihood(x);
    double llr = -2 * (ll - mll);
    double llr2 = PoissonCalculator.logLikelihoodRatio(u, x);
    System.out.printf("llr=%f, llr2=%f\n", llr, llr2);
    Assert.assertEquals("Log-likelihood ratio", llr, llr2, llr * 1e-10);
    double[] op = new double[x.length];
    for (int i = 0; i < n; i++) op[i] = PoissonCalculator.maximumLikelihood(x[i]);
    double max = Double.NEGATIVE_INFINITY;
    double maxa = 0;
    //TestSettings.setLogLevel(gdsc.smlm.TestSettings.LogLevel.DEBUG);
    int df = n - 1;
    ChiSquaredDistributionTable table = ChiSquaredDistributionTable.createUpperTailed(0.05, df);
    ChiSquaredDistributionTable table2 = ChiSquaredDistributionTable.createUpperTailed(0.001, df);
    System.out.printf("Chi2 = %f (q=%.3f), %f (q=%.3f)  %f %b  %f\n", 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));
    for (int i = 5; i <= 15; 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);
        llr = PoissonCalculator.logLikelihoodRatio(u, x);
        BigDecimal product = new BigDecimal(1);
        double ll2 = 0;
        for (int j = 0; j < n; j++) {
            double p1 = PoissonCalculator.likelihood(u[j], x[j]);
            ll2 += Math.log(p1);
            double ratio = p1 / op[j];
            product = product.multiply(new BigDecimal(ratio));
        }
        llr2 = -2 * Math.log(product.doubleValue());
        double p = ChiSquaredDistributionTable.computePValue(llr, df);
        double q = ChiSquaredDistributionTable.computeQValue(llr, df);
        TestSettings.info("a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f, q=%f (reject=%b @ %.3f, reject=%b @ %.3f)\n", 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());
        if (max < ll) {
            max = ll;
            maxa = a[0];
        }
        // too small to store in a double.
        if (product.doubleValue() > 0) {
            Assert.assertEquals("Log-likelihood", ll, ll2, Math.abs(ll2) * 1e-10);
            Assert.assertEquals("Log-likelihood ratio", llr, llr2, Math.abs(llr) * 1e-10);
        }
    }
    Assert.assertEquals("max", 1, maxa, 0);
}

Example 17 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method instanceAndFastMethodIsApproximatelyEqualToStaticMethod.

@Test
public void instanceAndFastMethodIsApproximatelyEqualToStaticMethod() {
    DoubleEquality eq = new DoubleEquality(3e-4, 0);
    RandomGenerator rg = new Well19937c(30051977);
    // Test for different x. The calculator approximation begins
    int n = 100;
    double[] u = new double[n];
    double[] x = new double[n];
    double e, o;
    for (double testx : new double[] { Math.nextAfter(PoissonCalculator.APPROXIMATION_X, -1), PoissonCalculator.APPROXIMATION_X, Math.nextUp(PoissonCalculator.APPROXIMATION_X), PoissonCalculator.APPROXIMATION_X * 1.1, PoissonCalculator.APPROXIMATION_X * 2, PoissonCalculator.APPROXIMATION_X * 10 }) {
        String X = Double.toString(testx);
        Arrays.fill(x, testx);
        PoissonCalculator pc = new PoissonCalculator(x);
        e = PoissonCalculator.maximumLogLikelihood(x);
        o = pc.getMaximumLogLikelihood();
        System.out.printf("[%s] Instance MaxLL = %g vs %g (error = %g)\n", X, e, o, DoubleEquality.relativeError(e, o));
        Assert.assertTrue("Instance Max LL not equal", eq.almostEqualRelativeOrAbsolute(e, o));
        o = PoissonCalculator.fastMaximumLogLikelihood(x);
        System.out.printf("[%s] Fast MaxLL = %g vs %g (error = %g)\n", X, e, o, DoubleEquality.relativeError(e, o));
        Assert.assertTrue("Fast Max LL not equal", eq.almostEqualRelativeOrAbsolute(e, o));
        // Generate data around the value
        for (int i = 0; i < n; i++) u[i] = x[i] + rg.nextDouble() - 0.5;
        e = PoissonCalculator.logLikelihood(u, x);
        o = pc.logLikelihood(u);
        System.out.printf("[%s] Instance LL = %g vs %g (error = %g)\n", X, e, o, DoubleEquality.relativeError(e, o));
        Assert.assertTrue("Instance LL not equal", eq.almostEqualRelativeOrAbsolute(e, o));
        o = PoissonCalculator.fastLogLikelihood(u, x);
        System.out.printf("[%s] Fast LL = %g vs %g (error = %g)\n", X, e, o, DoubleEquality.relativeError(e, o));
        Assert.assertTrue("Fast LL not equal", eq.almostEqualRelativeOrAbsolute(e, o));
        e = PoissonCalculator.logLikelihoodRatio(u, x);
        o = pc.getLogLikelihoodRatio(o);
        System.out.printf("[%s] Instance LLR = %g vs %g (error = %g)\n", X, e, o, DoubleEquality.relativeError(e, o));
        Assert.assertTrue("Instance LLR not equal", eq.almostEqualRelativeOrAbsolute(e, o));
    }
}

Example 18 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.

the class PoissonGradientProcedureTest method gradientProcedureIsFasterUnrolledThanGradientProcedure.

private void gradientProcedureIsFasterUnrolledThanGradientProcedure(final int nparams, final boolean precomputed) {
    org.junit.Assume.assumeTrue(speedTests || TestSettings.RUN_SPEED_TESTS);
    final int iter = 100;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    final ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    createFakeParams(nparams, iter, paramsList);
    // Remove the timing of the function call by creating a dummy function
    FakeGradientFunction f = new FakeGradientFunction(blockWidth, nparams);
    final Gradient1Function func = (precomputed) ? PrecomputedGradient1Function.wrapGradient1Function(f, Utils.newArray(f.size(), 0.1, 1.3)) : f;
    for (int i = 0; i < paramsList.size(); i++) {
        PoissonGradientProcedure p1 = new PoissonGradientProcedure(func);
        p1.computeFisherInformation(paramsList.get(i));
        p1.computeFisherInformation(paramsList.get(i));
        PoissonGradientProcedure p2 = PoissonGradientProcedureFactory.create(func);
        p2.computeFisherInformation(paramsList.get(i));
        p2.computeFisherInformation(paramsList.get(i));
        // Check they are the same
        Assert.assertArrayEquals("M " + i, p1.getLinear(), p2.getLinear(), 0);
    }
    // Realistic loops for an optimisation
    final int loops = 15;
    // Run till stable timing
    Timer t1 = new Timer() {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                PoissonGradientProcedure p1 = new PoissonGradientProcedure(func);
                for (int j = loops; j-- > 0; ) p1.computeFisherInformation(paramsList.get(k++ % iter));
            }
        }
    };
    long time1 = t1.getTime();
    Timer t2 = new Timer(t1.loops) {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                PoissonGradientProcedure p2 = PoissonGradientProcedureFactory.create(func);
                for (int j = loops; j-- > 0; ) p2.computeFisherInformation(paramsList.get(k++ % iter));
            }
        }
    };
    long time2 = t2.getTime();
    log("Precomputed=%b : Standard %d : Unrolled %d = %d : %fx\n", precomputed, time1, nparams, time2, (1.0 * time1) / time2);
    Assert.assertTrue(time2 < time1);
}

Example 19 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.

the class ChiSquaredDistributionTableTest method canPerformChiSquaredTest.

@Test
public void canPerformChiSquaredTest() {
    ChiSquareTest test = new ChiSquareTest();
    for (int n : new int[] { 10, 50, 100 }) {
        double[] x = Utils.newArray(n, 0.5, 1.0);
        long[] l = new long[x.length];
        RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
        for (int i = 0; i < x.length; i++) l[i] = rdg.nextPoisson(x[i]);
        double chi2 = test.chiSquare(x, l);
        double ep = test.chiSquareTest(x, l);
        int df = x.length - 1;
        double o = ChiSquaredDistributionTable.computeQValue(chi2, df);
        Assert.assertEquals(ep, o, 1e-10);
        ChiSquaredDistributionTable upperTable = ChiSquaredDistributionTable.createUpperTailed(o, df);
        double upper = chi2 * 1.01;
        double lower = chi2 * 0.99;
        Assert.assertTrue("Upper did not reject higher", upperTable.reject(upper, df));
        Assert.assertFalse("Upper did not reject actual value", upperTable.reject(o, df));
        Assert.assertFalse("Upper did not accept lower", upperTable.reject(lower, df));
    }
}

Example 20 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.

the class ErfTest method erfxxHasLowError.

private void erfxxHasLowError(BaseErf erf, double expected) {
    RandomGenerator rg = new Well19937c(30051977);
    int range = 3;
    double max = 0;
    for (int xi = -range; xi <= range; xi++) {
        for (int xi2 = -range; xi2 <= range; xi2++) {
            for (int i = 0; i < 5; i++) {
                double x = xi + rg.nextDouble();
                for (int j = 0; j < 5; j++) {
                    double x2 = xi2 + rg.nextDouble();
                    double o = erf.erf(x, x2);
                    double e = org.apache.commons.math3.special.Erf.erf(x, x2);
                    double error = Math.abs(o - e);
                    if (max < error)
                        max = error;
                    //System.out.printf("x=%f, x2=%f, e=%f, o=%f, error=%f\n", x, x2, e, o, error);
                    Assert.assertTrue(error < expected);
                }
            }
        }
    }
    System.out.printf("erfxx %s max error = %g\n", erf.name, max);
}