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

Example 56 with Well19937c

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

the class FRCTest method canComputeMirrored.

@Test
public void canComputeMirrored() {
    // Sample lines through an image to create a structure.
    int size = 1024;
    double[][] data = new double[size * 2][];
    RandomGenerator r = new Well19937c(30051977);
    for (int x = 0, y = 0, y2 = size, i = 0; x < size; x++, y++, y2--) {
        data[i++] = new double[] { x + r.nextGaussian() * 5, y + r.nextGaussian() * 5 };
        data[i++] = new double[] { x + r.nextGaussian() * 5, y2 + r.nextGaussian() * 5 };
    }
    // Create 2 images
    Rectangle bounds = new Rectangle(0, 0, size, size);
    IJImagePeakResults i1 = createImage(bounds);
    IJImagePeakResults i2 = createImage(bounds);
    int[] indices = Utils.newArray(data.length, 0, 1);
    MathArrays.shuffle(indices, r);
    for (int i : indices) {
        IJImagePeakResults image = i1;
        i1 = i2;
        i2 = image;
        image.add((float) data[i][0], (float) data[i][1], 1);
    }
    i1.end();
    i2.end();
    ImageProcessor ip1 = i1.getImagePlus().getProcessor();
    ImageProcessor ip2 = i2.getImagePlus().getProcessor();
    // Test
    FRC frc = new FRC();
    FloatProcessor[] fft1, fft2;
    fft1 = frc.getComplexFFT(ip1);
    fft2 = frc.getComplexFFT(ip2);
    float[] dataA1 = (float[]) fft1[0].getPixels();
    float[] dataB1 = (float[]) fft1[1].getPixels();
    float[] dataA2 = (float[]) fft2[0].getPixels();
    float[] dataB2 = (float[]) fft2[1].getPixels();
    float[] numeratorE = new float[dataA1.length];
    float[] absFFT1E = new float[dataA1.length];
    float[] absFFT2E = new float[dataA1.length];
    FRC.compute(numeratorE, absFFT1E, absFFT2E, dataA1, dataB1, dataA2, dataB2);
    Assert.assertTrue("numeratorE", FRC.checkSymmetry(numeratorE, size));
    Assert.assertTrue("absFFT1E", FRC.checkSymmetry(absFFT1E, size));
    Assert.assertTrue("absFFT2E", FRC.checkSymmetry(absFFT2E, size));
    float[] numeratorA = new float[dataA1.length];
    float[] absFFT1A = new float[dataA1.length];
    float[] absFFT2A = new float[dataA1.length];
    FRC.computeMirrored(size, numeratorA, absFFT1A, absFFT2A, dataA1, dataB1, dataA2, dataB2);
    //for (int y=0, i=0; y<size; y++)
    //	for (int x=0; x<size; x++, i++)
    //	{
    //		System.out.printf("[%d,%d = %d] %f ?= %f\n", x, y, i, numeratorE[i], numeratorA[i]);
    //	}
    Assert.assertArrayEquals("numerator", numeratorE, numeratorA, 0);
    Assert.assertArrayEquals("absFFT1", absFFT1E, absFFT1A, 0);
    Assert.assertArrayEquals("absFFT2", absFFT2E, absFFT2A, 0);
    FRC.computeMirroredFast(size, numeratorA, absFFT1A, absFFT2A, dataA1, dataB1, dataA2, dataB2);
    // Check this.
    for (int y = 1; y < size; y++) for (int x = 1, i = y * size + 1; x < size; x++, i++) {
        Assert.assertEquals("numerator", numeratorE[i], numeratorA[i], 0);
        Assert.assertEquals("absFFT1", absFFT1E[i], absFFT1A[i], 0);
        Assert.assertEquals("absFFT2", absFFT2E[i], absFFT2A[i], 0);
    }
}

Example 57 with Well19937c

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

the class FRCTest method computeMirroredIsFaster.

@Test
public void computeMirroredIsFaster() {
    // Sample lines through an image to create a structure.
    final int size = 2048;
    double[][] data = new double[size * 2][];
    RandomGenerator r = new Well19937c(30051977);
    for (int x = 0, y = 0, y2 = size, i = 0; x < size; x++, y++, y2--) {
        data[i++] = new double[] { x + r.nextGaussian() * 5, y + r.nextGaussian() * 5 };
        data[i++] = new double[] { x + r.nextGaussian() * 5, y2 + r.nextGaussian() * 5 };
    }
    // Create 2 images
    Rectangle bounds = new Rectangle(0, 0, size, size);
    IJImagePeakResults i1 = createImage(bounds);
    IJImagePeakResults i2 = createImage(bounds);
    int[] indices = Utils.newArray(data.length, 0, 1);
    MathArrays.shuffle(indices, r);
    for (int i : indices) {
        IJImagePeakResults image = i1;
        i1 = i2;
        i2 = image;
        image.add((float) data[i][0], (float) data[i][1], 1);
    }
    i1.end();
    i2.end();
    ImageProcessor ip1 = i1.getImagePlus().getProcessor();
    ImageProcessor ip2 = i2.getImagePlus().getProcessor();
    // Test
    FRC frc = new FRC();
    FloatProcessor[] fft1, fft2;
    fft1 = frc.getComplexFFT(ip1);
    fft2 = frc.getComplexFFT(ip2);
    final float[] dataA1 = (float[]) fft1[0].getPixels();
    final float[] dataB1 = (float[]) fft1[1].getPixels();
    final float[] dataA2 = (float[]) fft2[0].getPixels();
    final float[] dataB2 = (float[]) fft2[1].getPixels();
    final float[] numerator = new float[dataA1.length];
    final float[] absFFT1 = new float[dataA1.length];
    final float[] absFFT2 = new float[dataA1.length];
    TimingService ts = new TimingService(10);
    ts.execute(new MyTimingTask("compute") {

        public Object run(Object data) {
            FRC.compute(numerator, absFFT1, absFFT2, dataA1, dataB1, dataA2, dataB2);
            return null;
        }
    });
    ts.execute(new MyTimingTask("computeMirrored") {

        public Object run(Object data) {
            FRC.computeMirrored(size, numerator, absFFT1, absFFT2, dataA1, dataB1, dataA2, dataB2);
            return null;
        }
    });
    ts.execute(new MyTimingTask("computeMirroredFast") {

        public Object run(Object data) {
            FRC.computeMirroredFast(size, numerator, absFFT1, absFFT2, dataA1, dataB1, dataA2, dataB2);
            return null;
        }
    });
    ts.repeat(ts.getSize());
    ts.report();
}

Example 58 with Well19937c

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

the class SCMOSLikelihoodWrapperTest method functionComputesTargetGradient.

private void functionComputesTargetGradient(Gaussian2DFunction f1, int targetParameter, double threshold) {
    int[] indices = f1.gradientIndices();
    int gradientIndex = findGradientIndex(f1, targetParameter);
    double[] dyda = new double[indices.length];
    double[] a;
    SCMOSLikelihoodWrapper ff1;
    int n = maxx * maxx;
    int count = 0, total = 0;
    RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    for (double background : testbackground) for (double signal1 : testsignal1) for (double angle1 : testangle1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
        a = createParameters(background, signal1, angle1, cx1, cy1, w1[0], w1[1]);
        // Create y as a function we would want to move towards
        double[] a2 = a.clone();
        a2[targetParameter] *= 1.3;
        f1.initialise(a2);
        double[] data = new double[n];
        for (int i = 0; i < n; i++) {
            // Simulate sCMOS camera
            double u = f1.eval(i);
            data[i] = rdg.nextPoisson(u) * g[i] + rdg.nextGaussian(o[i], sd[i]);
        }
        ff1 = new SCMOSLikelihoodWrapper(f1, a, data, n, var, g, o);
        // Numerically solve gradient. 
        // Calculate the step size h to be an exact numerical representation
        final double xx = a[targetParameter];
        // Get h to minimise roundoff error
        double h = Precision.representableDelta(xx, h_);
        ff1.likelihood(getVariables(indices, a), dyda);
        // Evaluate at (x+h) and (x-h)
        a[targetParameter] = xx + h;
        double value2 = ff1.likelihood(getVariables(indices, a));
        a[targetParameter] = xx - h;
        double value3 = ff1.likelihood(getVariables(indices, a));
        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)
            Assert.assertTrue(NAME[targetParameter] + ": " + gradient + " != " + dyda[gradientIndex], ok);
        ok = eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]);
        if (ok)
            count++;
        total++;
    }
    double p = (100.0 * count) / total;
    logf("%s : %s = %d / %d (%.2f)\n", f1.getClass().getSimpleName(), NAME[targetParameter], count, total, p);
    Assert.assertTrue(NAME[targetParameter] + " fraction too low: " + p, p > threshold);
}

Example 59 with Well19937c

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

the class PrecomputedFunctionTest method precomputedGradient1FunctionWrapsPrecomputedValues.

@Test
public void precomputedGradient1FunctionWrapsPrecomputedValues() {
    int n = 3;
    RandomGenerator r = new Well19937c(30051977);
    Gradient1Function f0 = new FakeGradientFunction(3, n);
    int size = f0.size();
    double[] b1 = new PseudoRandomGenerator(size, r).getSequence();
    double[] b2 = new PseudoRandomGenerator(size, r).getSequence();
    Gradient1Function f1 = PrecomputedGradient1Function.wrapGradient1Function(f0, b1);
    Gradient1Function f2 = PrecomputedGradient1Function.wrapGradient1Function(f1, b2);
    double[] p = new double[n];
    for (int i = 0; i < n; i++) p[i] = r.nextDouble();
    double[] d0 = new double[n];
    double[] d1 = new double[n];
    double[] d2 = new double[n];
    double[] v0 = evaluateGradient1Function(f0, p, d0);
    double[] v1 = evaluateGradient1Function(f1, p, d1);
    double[] v2 = evaluateGradient1Function(f2, p, d2);
    for (int i = 0; i < v0.length; i++) {
        double e = v0[i] + b1[i] + b2[i];
        double o1 = v1[i] + b2[i];
        double o2 = v2[i];
        Assert.assertEquals("o1", e, o1, 0);
        Assert.assertEquals("o2", e, o2, 1e-6);
    }
    Assert.assertArrayEquals("d1", d0, d1, 0);
    Assert.assertArrayEquals("d2", d0, d2, 0);
}

Example 60 with Well19937c

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

the class PrecomputedFunctionTest method precomputedValueFunctionWrapsPrecomputedValues.

@Test
public void precomputedValueFunctionWrapsPrecomputedValues() {
    int n = 3;
    RandomGenerator r = new Well19937c(30051977);
    ValueFunction f0 = new FakeGradientFunction(3, n);
    int size = f0.size();
    double[] b1 = new PseudoRandomGenerator(size, r).getSequence();
    double[] b2 = new PseudoRandomGenerator(size, r).getSequence();
    ValueFunction f1 = PrecomputedValueFunction.wrapValueFunction(f0, b1);
    ValueFunction f2 = PrecomputedValueFunction.wrapValueFunction(f1, b2);
    double[] p = new double[n];
    for (int i = 0; i < n; i++) p[i] = r.nextDouble();
    double[] v0 = evaluateValueFunction(f0, p);
    double[] v1 = evaluateValueFunction(f1, p);
    double[] v2 = evaluateValueFunction(f2, p);
    for (int i = 0; i < v0.length; i++) {
        double e = v0[i] + b1[i] + b2[i];
        double o1 = v1[i] + b2[i];
        double o2 = v2[i];
        Assert.assertEquals("o1", e, o1, 0);
        Assert.assertEquals("o2", e, o2, 1e-6);
    }
}