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

Example 1 with ContinuousSampler

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

the class BaseFunctionSolverTest method createData.

private static double[][] createData(RandomSeed source) {
    // Per observation read noise.
    // This is generated once so create the randon generator here.
    final UniformRandomProvider rg = RngUtils.create(source.getSeed());
    final ContinuousSampler ed = SamplerUtils.createExponentialSampler(rg, variance);
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, gain, gainSD);
    final double[] w = new double[size * size];
    final double[] n = new double[size * size];
    for (int i = 0; i < w.length; i++) {
        final double pixelVariance = ed.sample();
        final double pixelGain = Math.max(0.1, gs.sample());
        // weights = var / g^2
        w[i] = pixelVariance / (pixelGain * pixelGain);
        // Convert to standard deviation for simulation
        n[i] = Math.sqrt(w[i]);
    }
    return new double[][] { w, n };
}

Example 2 with ContinuousSampler

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

the class GaussianFilterTest method filter1IsSameAsFilter2.

private void filter1IsSameAsFilter2(RandomSeed seed, GFilter f1, GFilter f2, boolean weighted, double tolerance) {
    final UniformRandomProvider rand = RngUtils.create(seed.getSeed());
    final float[] data = createData(rand, size, size);
    float[] weights = null;
    if (weighted) {
        final ContinuousSampler ed = SamplerUtils.createExponentialSampler(rand, 57);
        weights = new float[data.length];
        for (int i = 0; i < weights.length; i++) {
            weights[i] = (float) (1.0 / Math.max(0.01, ed.sample()));
        }
        // w[i] = (float) (1.0 / Math.max(0.01, rand.nextGaussian() * 0.2 + 2));
        // w[i] = 0.5f;
        f1.setWeights(weights);
        f2.setWeights(weights);
    }
    for (final double sigma : sigmas) {
        final float[] e = data.clone();
        f2.run(e, sigma);
        final float[] o = data.clone();
        f1.run(o, sigma);
        double max = 0;
        for (int i = 0; i < e.length; i++) {
            final double d = DoubleEquality.relativeError(e[i], o[i]);
            if (max < d) {
                max = d;
            }
        }
        logger.fine(FunctionUtils.getSupplier("%s vs %s w=%b @ %.1f = %g", f1.getName(), f2.getName(), weighted, sigma, max));
        Assertions.assertTrue(max < tolerance);
    }
}

Example 3 with ContinuousSampler

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

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

the class CmosAnalysis method simulate.

private void simulate() throws IOException {
    // Create the offset, variance and gain for each pixel
    final int n = settings.size * settings.size;
    final float[] pixelOffset = new float[n];
    final float[] pixelVariance = new float[n];
    final float[] pixelGain = new float[n];
    IJ.showStatus("Creating random per-pixel readout");
    final long start = System.currentTimeMillis();
    final UniformRandomProvider rg = UniformRandomProviders.create();
    final DiscreteSampler pd = PoissonSamplerUtils.createPoissonSampler(rg, settings.offset);
    final ContinuousSampler ed = SamplerUtils.createExponentialSampler(rg, settings.variance);
    final SharedStateContinuousSampler gauss = SamplerUtils.createGaussianSampler(rg, settings.gain, settings.gainStdDev);
    Ticker ticker = ImageJUtils.createTicker(n, 0);
    for (int i = 0; i < n; i++) {
        // Q. Should these be clipped to a sensible range?
        pixelOffset[i] = pd.sample();
        pixelVariance[i] = (float) ed.sample();
        pixelGain[i] = (float) gauss.sample();
        ticker.tick();
    }
    IJ.showProgress(1);
    // Save to the directory as a stack
    final ImageStack simulationStack = new ImageStack(settings.size, settings.size);
    simulationStack.addSlice("Offset", pixelOffset);
    simulationStack.addSlice("Variance", pixelVariance);
    simulationStack.addSlice("Gain", pixelGain);
    simulationImp = new ImagePlus("PerPixel", simulationStack);
    // Only the info property is saved to the TIFF file
    simulationImp.setProperty("Info", String.format("Offset=%s; Variance=%s; Gain=%s +/- %s", MathUtils.rounded(settings.offset), MathUtils.rounded(settings.variance), MathUtils.rounded(settings.gain), MathUtils.rounded(settings.gainStdDev)));
    IJ.save(simulationImp, new File(settings.directory, "perPixelSimulation.tif").getPath());
    // Create thread pool and workers
    final int threadCount = getThreads();
    final ExecutorService executor = Executors.newFixedThreadPool(threadCount);
    final LocalList<Future<?>> futures = new LocalList<>(numberOfThreads);
    // Simulate the exposure input.
    final int[] photons = settings.getPhotons();
    // For saving stacks
    final int blockSize = 10;
    int numberPerThread = (int) Math.ceil((double) settings.frames / numberOfThreads);
    // Convert to fit the block size
    numberPerThread = (int) Math.ceil((double) numberPerThread / blockSize) * blockSize;
    final Pcg32 rng = Pcg32.xshrs(start);
    // Note the bias is increased by 3-fold so add 2 to the length
    ticker = Ticker.createStarted(new ImageJTrackProgress(true), (long) (photons.length + 2) * settings.frames, threadCount > 1);
    for (final int p : photons) {
        ImageJUtils.showStatus(() -> "Simulating " + TextUtils.pleural(p, "photon"));
        // Create the directory
        final Path out = Paths.get(settings.directory, String.format("photon%03d", p));
        Files.createDirectories(out);
        // Increase frames for bias image
        final int frames = settings.frames * (p == 0 ? 3 : 1);
        for (int from = 0; from < frames; ) {
            final int to = Math.min(from + numberPerThread, frames);
            futures.add(executor.submit(new SimulationWorker(ticker, rng.split(), out.toString(), simulationStack, from, to, blockSize, p)));
            from = to;
        }
        ConcurrencyUtils.waitForCompletionUnchecked(futures);
        futures.clear();
    }
    final String msg = "Simulation time = " + TextUtils.millisToString(System.currentTimeMillis() - start);
    IJ.showStatus(msg);
    ImageJUtils.clearSlowProgress();
    executor.shutdown();
    ImageJUtils.log(msg);
}

Example 5 with ContinuousSampler

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

the class StandardFluorophoreSequenceModel method init.

private void init(double startT, double onTime, double offTime, double offTime2, double blinks1, double blinks2, boolean useGeometricBlinkingDistribution, UniformRandomProvider rand) {
    // Model two dark states: short and long. The second offTime and blinks1 is for the long dark
    // state:
    // 
    // ++-+-+++-+.................+-+--++-+................................+--+++-+
    // 
    // + = on
    // - = Short dark state
    // . = Long dark state
    // Note: 1+blinks1 is the number of on-states
    final TDoubleArrayList sequence = new TDoubleArrayList();
    // The exponential distribution is just scaled by the mean
    final ContinuousSampler sampler = SamplerUtils.createExponentialSampler(rand);
    // Perform a set number of long blinks
    final int nLongBlinks = getBlinks(useGeometricBlinkingDistribution, rand, blinks2);
    double time = startT;
    for (int n = 0; n <= nLongBlinks; n++) {
        // For each burst between long blinks perform a number of short blinks
        final int nShortBlinks = getBlinks(useGeometricBlinkingDistribution, rand, blinks1);
        // Starts on the current time
        sequence.add(time);
        // Stops after the on-time
        time += sampler.sample() * onTime;
        sequence.add(time);
        // Remaining bursts
        for (int i = 0; i < nShortBlinks; i++) {
            // Next burst starts after the short off-time
            time += sampler.sample() * offTime;
            sequence.add(time);
            // Stops after the on-time
            time += sampler.sample() * onTime;
            sequence.add(time);
        }
        // Add the long dark state if there are more bursts.
        time += sampler.sample() * offTime2;
    }
    // Convert the sequence to the burst sequence array
    setBurstSequence(sequence.toArray());
}