use of org.apache.commons.rng.sampling.distribution.DiscreteSampler in project GDSC-SMLM by aherbert.
the class BinomialFitterTest method createData.
private static int[] createData(UniformRandomProvider rg, int n, double p, boolean zeroTruncated) {
final BinomialDistribution bd = new BinomialDistribution(null, n, p);
final DiscreteSampler sampler = new InverseTransformDiscreteSampler(rg, pvalue -> bd.inverseCumulativeProbability(pvalue));
final int[] data = new int[2000];
if (zeroTruncated) {
if (p <= 0) {
throw new RuntimeException("p must be positive");
}
for (int i = 0; i < data.length; i++) {
int count;
do {
count = sampler.sample();
} while (count == 0);
data[i] = count;
}
} else {
for (int i = 0; i < data.length; i++) {
data[i] = sampler.sample();
}
}
return data;
}
use of org.apache.commons.rng.sampling.distribution.DiscreteSampler 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;
}
use of org.apache.commons.rng.sampling.distribution.DiscreteSampler 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);
}
Aggregations