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

Example 11 with SharedStateContinuousSampler

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

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

the class MultivariateGaussianMixtureExpectationMaximizationTest method create.

/**
 * Creates a n-dimension array with values in the given range.
 *
 * @param n the array dimension
 * @param rng the random generator
 * @param lower the lower
 * @param upper the upper
 * @return the data
 */
private static double[] create(int n, UniformRandomProvider rng, double lower, double upper) {
    final SharedStateContinuousSampler sampler = ContinuousUniformSampler.of(rng, lower, upper);
    final double[] data = new double[n];
    for (int i = 0; i < n; i++) {
        data[i] = sampler.sample();
    }
    return data;
}

Example 13 with SharedStateContinuousSampler

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

the class DynamicMultipleTargetTracingTest method testTraceMoleculesDisableIntensityModel.

/**
 * Test trace molecules using 2 molecules. One is fixed and the other moves past it. The tracing
 * should assign the fixed molecule correctly as it has a low local diffusion rate.
 */
@Test
void testTraceMoleculesDisableIntensityModel() {
    final UniformRandomProvider rng = RngUtils.create(125631236L);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    // localisation precision (in pixels)
    final double s = 0.1;
    final SharedStateContinuousSampler intensity = SamplerUtils.createGaussianSampler(rng, 1000, 100);
    final MemoryPeakResults results = new MemoryPeakResults(100);
    final CalibrationWriter writer = results.getCalibrationWriterSafe();
    // 0.1 um pixels, 1 second exposure time
    writer.setDistanceUnit(DistanceUnit.PIXEL);
    writer.setNmPerPixel(100);
    writer.setExposureTime(1000);
    results.setCalibration(writer.getCalibration());
    // First molecule diffuses roughly across the field from top-left to bottom-right.
    // 5 frames is the default for local stats, 15 frames for trajectory removal.
    // Use 20 so we build local stats and can expire a trajectory.
    final int size = 20;
    final float x1 = size / 2 + 0.5f;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (x1 + gauss.sample() * s), (float) (i + gauss.sample() * s), (float) (intensity.sample())));
    }
    // Second molecule is fixed in the centre with a same intensity
    final int x = size / 2;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (x + gauss.sample() * s), (float) (x + gauss.sample() * s), (float) (intensity.sample())));
    }
    // Add a single molecule that will not connect to anything in the second frame.
    // This should create a trajectory that will expire.
    results.add(new PeakResult(1, x1, size, (float) (intensity.sample())));
    // Move centre to centre each jump => 0.1 um or 0.01 um^2
    // MSD = 4D => D = 0.01 / 4 = 0.0025
    final DmttConfiguration config = DmttConfiguration.newBuilder(0.0025).setDisableIntensityModel(true).setTemporalWindow(10).build();
    final List<Trace> traces = new DynamicMultipleTargetTracing(results).traceMolecules(config);
    // Should have 3 traces
    Assertions.assertEquals(3, traces.size());
    // Assert ids start from 1
    for (int i = 0; i < traces.size(); i++) {
        Assertions.assertEquals(i + 1, traces.get(i).getId());
    }
    // Traces should be 2 full length and 1 single peak
    Assertions.assertEquals(size, traces.get(0).size());
    Assertions.assertEquals(size, traces.get(1).size());
    Assertions.assertEquals(1, traces.get(2).size());
    // Do an analysis on the actual tracks.
    // One should be based in the centre and the other should have parts close to position (i,i)
    // for each frame i.
    final PeakResult[] peaks = results.toArray();
    // Assume traces are initially created using the input order of the results.
    final Trace t1 = traces.get(0);
    final Trace t2 = traces.get(1);
    for (int i = 0; i < size; i++) {
        Assertions.assertSame(peaks[i], t1.get(i));
        Assertions.assertSame(peaks[i + size], t2.get(i));
    }
}

Example 14 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project gdsc by aherbert.

the class FindFociTest method combine.

private static short[] combine(UniformRandomProvider rg, float[] data1, float[] data2, float[] data3) {
    // Combine images and add a bias and read noise
    final SharedStateContinuousSampler g = SamplerUtils.createGaussianSampler(rg, BIAS, 5);
    final short[] data = new short[data1.length];
    for (int i = 0; i < data.length; i++) {
        final double mu = data1[i] + data2[i] + data3[i];
        double value = g.sample();
        if (mu != 0) {
            value += new PoissonSampler(rg, mu).sample();
        }
        if (value < 0) {
            value = 0;
        } else if (value > 65535) {
            value = 65535;
        }
        data[i] = (short) value;
    }
    return data;
}

Example 15 with SharedStateContinuousSampler

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

the class CompositeSamplersTest method testSharedStateContinuousSampler.

/**
 * Test the SharedStateSampler implementation for the composite
 * SharedStateContinuousSampler.
 *
 * @param customFactory Set to true to use a custom discrete sampler factory that does not
 * support a shared stated sampler.
 */
private static void testSharedStateContinuousSampler(boolean customFactory) {
    final UniformRandomProvider rng1 = RandomSource.SPLIT_MIX_64.create(0L);
    final UniformRandomProvider rng2 = RandomSource.SPLIT_MIX_64.create(0L);
    final Builder<SharedStateContinuousSampler> builder = CompositeSamplers.newSharedStateContinuousSamplerBuilder();
    if (customFactory) {
        addFactoryWithNoSharedStateSupport(builder);
    }
    // Sample within the ranges between the ticks
    final double[] ticks = { 7.89, 13.99, 21.7, 35.6, 45.5 };
    for (int i = 1; i < ticks.length; i++) {
        final DoubleRangeSampler sampler = new DoubleRangeSampler(rng1, ticks[i - 1], ticks[i]);
        // Weight using the range
        builder.add(sampler, sampler.range());
    }
    final SharedStateContinuousSampler sampler1 = builder.build(rng1);
    final SharedStateContinuousSampler sampler2 = sampler1.withUniformRandomProvider(rng2);
    RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}