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

Example 6 with RandomDataGenerator

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

the class SCMOSLikelihoodWrapperTest method functionComputesTargetGradientPerDatum.

private void functionComputesTargetGradientPerDatum(Gaussian2DFunction f1, int targetParameter) {
    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.1;
        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_);
        for (int x : testx) for (int y : testy) {
            int i = y * maxx + x;
            a[targetParameter] = xx;
            ff1.likelihood(getVariables(indices, a), dyda, i);
            // Evaluate at (x+h) and (x-h)
            a[targetParameter] = xx + h;
            double value2 = ff1.likelihood(getVariables(indices, a), i);
            a[targetParameter] = xx - h;
            double value3 = ff1.likelihood(getVariables(indices, a), i);
            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 = eqPerDatum.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]);
            if (ok)
                count++;
            total++;
        }
    }
    double p = (100.0 * count) / total;
    logf("Per Datum %s : %s = %d / %d (%.2f)\n", f1.getClass().getSimpleName(), NAME[targetParameter], count, total, p);
    Assert.assertTrue(NAME[targetParameter] + " fraction too low per datum: " + p, p > 90);
}

Example 7 with RandomDataGenerator

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

the class SCMOSLikelihoodWrapperTest method canComputePValue.

private void canComputePValue(BaseNonLinearFunction nlf) {
    System.out.println(nlf.name);
    int n = maxx * maxx;
    double[] a = new double[] { 1 };
    // Simulate sCMOS camera
    nlf.initialise(a);
    RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    double[] k = Utils.newArray(n, 0, 1.0);
    for (int i = 0; i < n; i++) {
        double u = nlf.eval(i);
        if (u > 0)
            u = rdg.nextPoisson(u);
        k[i] = u * g[i] + rdg.nextGaussian(o[i], sd[i]);
    }
    SCMOSLikelihoodWrapper f = new SCMOSLikelihoodWrapper(nlf, a, k, n, var, g, o);
    double oll = f.computeObservedLikelihood();
    double oll2 = 0;
    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]);
    }
    System.out.printf("oll=%f, oll2=%f\n", oll, oll2);
    Assert.assertEquals("Observed Log-likelihood", oll2, oll, oll2 * 1e-10);
    double min = Double.POSITIVE_INFINITY;
    double mina = 0;
    for (int i = 5; i <= 15; i++) {
        a[0] = (double) i / 10;
        double ll = f.likelihood(a);
        double llr = f.computeLogLikelihoodRatio(ll);
        BigDecimal product = new BigDecimal(1);
        double ll2 = 0;
        for (int j = 0; j < n; j++) {
            double p1 = SCMOSLikelihoodWrapper.likelihood(nlf.eval(j), var[j], g[j], o[j], k[j]);
            ll2 -= Math.log(p1);
            double ratio = p1 / op[j];
            product = product.multiply(new BigDecimal(ratio));
        }
        double llr2 = -2 * Math.log(product.doubleValue());
        double q = f.computeQValue(ll);
        System.out.printf("a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f\n", a[0], ll, ll2, llr, llr2, product.round(new MathContext(4)).toString(), q);
        if (min > ll) {
            min = ll;
            mina = a[0];
        }
        // too small to store in a double.
        if (product.doubleValue() > 0)
            Assert.assertEquals("Log-likelihood", llr, llr2, llr * 1e-10);
    }
    Assert.assertEquals("min", 1, mina, 0);
}

Example 8 with RandomDataGenerator

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

the class CreateData method createBackground.

private float[] createBackground(RandomDataGenerator random) {
    float[] pixels2 = null;
    if (settings.background > 0) {
        if (random == null)
            random = new RandomDataGenerator(createRandomGenerator());
        createBackgroundPixels();
        pixels2 = Arrays.copyOf(backgroundPixels, backgroundPixels.length);
        // Add Poisson noise.
        if (uniformBackground) {
            // We can do N random samples thus ensuring the background average is constant.
            // Note: The number of samples must be Poisson distributed.
            // currently: pixels2[0] = uniform background level
            // => (pixels2[0] * pixels2.length) = amount of photons that fall on the image. 
            final int samples = (int) random.nextPoisson(pixels2[0] * pixels2.length);
            // Only do sampling if the number of samples is valid
            if (samples >= 1) {
                pixels2 = new float[pixels2.length];
                final int upper = pixels2.length - 1;
                for (int i = 0; i < samples; i++) {
                    pixels2[random.nextInt(0, upper)] += 1;
                }
            } else {
                // If using a uniform illumination then we can use a fixed Poisson distribution
                PoissonDistribution dist = new PoissonDistribution(random.getRandomGenerator(), pixels2[0], PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
                for (int i = 0; i < pixels2.length; i++) {
                    pixels2[i] = dist.sample();
                }
            }
        } else {
            for (int i = 0; i < pixels2.length; i++) {
                pixels2[i] = random.nextPoisson(pixels2[i]);
            }
        }
    } else {
        pixels2 = new float[settings.size * settings.size];
    }
    return pixels2;
}

Example 9 with RandomDataGenerator

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

the class CreateData method drawImage.

//StoredDataStatistics rawPhotons = new StoredDataStatistics();
//StoredDataStatistics drawPhotons = new StoredDataStatistics();
//	private synchronized void addRaw(double d)
//	{
//		//rawPhotons.add(d);
//	}
//
//	private synchronized void addDraw(double d)
//	{
//		//drawPhotons.add(d);
//	}
/**
	 * Create an image from the localisations using the configured PSF width. Draws a new stack
	 * image.
	 * <p>
	 * Note that the localisations are filtered using the signal. The input list of localisations will be updated.
	 * 
	 * @param localisationSets
	 * @return The localisations
	 */
private List<LocalisationModel> drawImage(final List<LocalisationModelSet> localisationSets) {
    if (localisationSets.isEmpty())
        return null;
    // Create a new list for all localisation that are drawn (i.e. pass the signal filters)
    List<LocalisationModelSet> newLocalisations = Collections.synchronizedList(new ArrayList<LocalisationModelSet>(localisationSets.size()));
    photonsRemoved = new AtomicInteger();
    t1Removed = new AtomicInteger();
    tNRemoved = new AtomicInteger();
    photonStats = new SummaryStatistics();
    // Add drawn spots to memory
    results = new MemoryPeakResults();
    Calibration c = new Calibration(settings.pixelPitch, settings.getTotalGain(), settings.exposureTime);
    c.setEmCCD((settings.getEmGain() > 1));
    c.setBias(settings.bias);
    c.setReadNoise(settings.readNoise * ((settings.getCameraGain() > 0) ? settings.getCameraGain() : 1));
    c.setAmplification(settings.getAmplification());
    results.setCalibration(c);
    results.setSortAfterEnd(true);
    results.begin();
    maxT = localisationSets.get(localisationSets.size() - 1).getTime();
    // Display image
    ImageStack stack = new ImageStack(settings.size, settings.size, maxT);
    final double psfSD = getPsfSD();
    if (psfSD <= 0)
        return null;
    ImagePSFModel imagePSFModel = null;
    if (imagePSF) {
        // Create one Image PSF model that can be copied
        imagePSFModel = createImagePSF(localisationSets);
        if (imagePSFModel == null)
            return null;
    }
    IJ.showStatus("Drawing image ...");
    // Multi-thread for speed
    // Note that the default Executors.newCachedThreadPool() will continue to make threads if
    // new tasks are added. We need to limit the tasks that can be added using a fixed size
    // blocking queue.
    // http://stackoverflow.com/questions/1800317/impossible-to-make-a-cached-thread-pool-with-a-size-limit
    // ExecutorService threadPool = Executors.newCachedThreadPool();
    ExecutorService threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
    List<Future<?>> futures = new LinkedList<Future<?>>();
    // Count all the frames to process
    frame = 0;
    totalFrames = maxT;
    // Collect statistics on the number of photons actually simulated
    // Process all frames
    int i = 0;
    int lastT = -1;
    for (LocalisationModelSet l : localisationSets) {
        if (Utils.isInterrupted())
            break;
        if (l.getTime() != lastT) {
            lastT = l.getTime();
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, i, lastT, createPSFModel(imagePSFModel), results, stack, poissonNoise, new RandomDataGenerator(createRandomGenerator()))));
        }
        i++;
    }
    // Finish processing data
    Utils.waitForCompletion(futures);
    futures.clear();
    if (Utils.isInterrupted()) {
        IJ.showProgress(1);
        return null;
    }
    // Do all the frames that had no localisations
    for (int t = 1; t <= maxT; t++) {
        if (Utils.isInterrupted())
            break;
        if (stack.getPixels(t) == null) {
            futures.add(threadPool.submit(new ImageGenerator(localisationSets, newLocalisations, maxT, t, null, results, stack, poissonNoise, new RandomDataGenerator(createRandomGenerator()))));
        }
    }
    // Finish
    Utils.waitForCompletion(futures);
    threadPool.shutdown();
    IJ.showProgress(1);
    if (Utils.isInterrupted()) {
        return null;
    }
    results.end();
    // Clear memory
    imagePSFModel = null;
    threadPool = null;
    futures.clear();
    futures = null;
    if (photonsRemoved.get() > 0)
        Utils.log("Removed %d localisations with less than %.1f rendered photons", photonsRemoved.get(), settings.minPhotons);
    if (t1Removed.get() > 0)
        Utils.log("Removed %d localisations with no neighbours @ SNR %.2f", t1Removed.get(), settings.minSNRt1);
    if (tNRemoved.get() > 0)
        Utils.log("Removed %d localisations with valid neighbours @ SNR %.2f", tNRemoved.get(), settings.minSNRtN);
    if (photonStats.getN() > 0)
        Utils.log("Average photons rendered = %s +/- %s", Utils.rounded(photonStats.getMean()), Utils.rounded(photonStats.getStandardDeviation()));
    //System.out.printf("rawPhotons = %f\n", rawPhotons.getMean());
    //System.out.printf("drawPhotons = %f\n", drawPhotons.getMean());
    //Utils.showHistogram("draw photons", drawPhotons, "photons", true, 0, 1000);
    // Update with all those localisation that have been drawn
    localisationSets.clear();
    localisationSets.addAll(newLocalisations);
    newLocalisations = null;
    IJ.showStatus("Displaying image ...");
    ImageStack newStack = stack;
    if (!settings.rawImage) {
        // Get the global limits and ensure all values can be represented
        Object[] imageArray = stack.getImageArray();
        float[] limits = Maths.limits((float[]) imageArray[0]);
        for (int j = 1; j < imageArray.length; j++) limits = Maths.limits(limits, (float[]) imageArray[j]);
        // Leave bias in place
        limits[0] = 0;
        // Check if the image will fit in a 16-bit range
        if ((limits[1] - limits[0]) < 65535) {
            // Convert to 16-bit
            newStack = new ImageStack(stack.getWidth(), stack.getHeight(), stack.getSize());
            // Account for rounding
            final float min = (float) (limits[0] - 0.5);
            for (int j = 0; j < imageArray.length; j++) {
                float[] image = (float[]) imageArray[j];
                short[] pixels = new short[image.length];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = (short) (image[k] - min);
                }
                newStack.setPixels(pixels, j + 1);
                // Free memory
                imageArray[j] = null;
                // Attempt to stay within memory (check vs 32MB)
                if (MemoryPeakResults.freeMemory() < 33554432L)
                    MemoryPeakResults.runGCOnce();
            }
        } else {
            // Keep as 32-bit but round to whole numbers
            for (int j = 0; j < imageArray.length; j++) {
                float[] pixels = (float[]) imageArray[j];
                for (int k = 0; k < pixels.length; k++) {
                    pixels[k] = Math.round(pixels[k]);
                }
            }
        }
    }
    // Show image
    ImagePlus imp = Utils.display(CREATE_DATA_IMAGE_TITLE, newStack);
    ij.measure.Calibration cal = new ij.measure.Calibration();
    String unit = "nm";
    double unitPerPixel = settings.pixelPitch;
    if (unitPerPixel > 100) {
        unit = "um";
        unitPerPixel /= 1000.0;
    }
    cal.setUnit(unit);
    cal.pixelHeight = cal.pixelWidth = unitPerPixel;
    imp.setCalibration(cal);
    imp.setDimensions(1, 1, newStack.getSize());
    imp.resetDisplayRange();
    imp.updateAndDraw();
    saveImage(imp);
    results.setSource(new IJImageSource(imp));
    results.setName(CREATE_DATA_IMAGE_TITLE + " (" + TITLE + ")");
    results.setConfiguration(createConfiguration((float) psfSD));
    results.setBounds(new Rectangle(0, 0, settings.size, settings.size));
    MemoryPeakResults.addResults(results);
    setBenchmarkResults(imp, results);
    if (benchmarkMode && benchmarkParameters != null)
        benchmarkParameters.setPhotons(results);
    List<LocalisationModel> localisations = toLocalisations(localisationSets);
    savePulses(localisations, results, CREATE_DATA_IMAGE_TITLE);
    // Saved the fixed and moving localisations into different datasets
    saveFixedAndMoving(results, CREATE_DATA_IMAGE_TITLE);
    return localisations;
}

Example 10 with RandomDataGenerator

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

the class EJMLLinearSolverTest method runInversionSpeedTest.

private void runInversionSpeedTest(int flags) {
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    int n = f0.size();
    final double[] y = new double[n];
    final TurboList<DenseMatrix64F> aList = new TurboList<DenseMatrix64F>();
    double[] testbackground = new double[] { 0.2, 0.7 };
    double[] testsignal1 = new double[] { 30, 100, 300 };
    double[] testcx1 = new double[] { 4.9, 5.3 };
    double[] testcy1 = new double[] { 4.8, 5.2 };
    double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    int np = f0.getNumberOfGradients();
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(np);
    final RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    //double lambda = 10;
    for (double background : testbackground) // Peak 1
    for (double signal1 : testsignal1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double w1 : testw1) {
        double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
        f0.initialise(p);
        f0.forEach(new ValueProcedure() {

            int i = 0;

            public void execute(double value) {
                // Poisson data 
                y[i++] = rdg.nextPoisson(value);
            }
        });
        double[][] alpha = new double[np][np];
        double[] beta = new double[np];
        //double ss = 
        calc.findLinearised(n, y, p, alpha, beta, f0);
        //System.out.printf("SS = %f\n", ss);
        // As per the LVM algorithm
        //for (int i = 0; i < np; i++)
        //	alpha[i][i] *= lambda;
        aList.add(EJMLLinearSolver.toA(alpha));
    }
    DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
    boolean[] ignore = new boolean[a.length];
    double[][] answer = new double[a.length][];
    int runs = 100000 / a.length;
    TimingService ts = new TimingService(runs);
    TurboList<InversionTimingTask> tasks = new TurboList<InversionTimingTask>();
    tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
    tasks.add(new LinearInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyLDLTInversionTimingTask(a, ignore, answer));
    tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
    tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
    for (InversionTimingTask task : tasks) if (!task.badSolver)
        ts.execute(task);
    ts.repeat();
    ts.report();
}