Examples with MarsagliaTsangGammaSampler uk.ac.sussex.gdsc.core.utils.rng.MarsagliaTsangGammaSampler used on opensource projects

Example 1 with MarsagliaTsangGammaSampler

use of uk.ac.sussex.gdsc.core.utils.rng.MarsagliaTsangGammaSampler in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method drawGaussian.

/**
 * Draw a Gaussian with Poisson shot noise and Gaussian read noise.
 *
 * @param params The Gaussian parameters
 * @param noise The read noise
 * @param noiseModel the noise model
 * @param rg the random generator
 * @return The data
 */
static double[] drawGaussian(double[] params, double[] noise, NoiseModel noiseModel, UniformRandomProvider rg) {
    final int n = params.length / Gaussian2DFunction.PARAMETERS_PER_PEAK;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(n, size, size, flags, null);
    final double[] data = f.computeValues(params);
    // Poisson noise
    for (int i = 0; i < data.length; i++) {
        if (data[i] > 0) {
            data[i] = GdscSmlmTestUtils.createPoissonSampler(rg, data[i]).sample();
        }
    }
    // Simulate EM-gain
    if (noiseModel == NoiseModel.EMCCD) {
        // Use a gamma distribution
        // Since the call random.nextGamma(...) creates a Gamma distribution
        // which pre-calculates factors only using the scale parameter we
        // create a custom gamma distribution where the shape can be set as a property.
        final MarsagliaTsangGammaSampler gd = new MarsagliaTsangGammaSampler(rg, 1, emGain);
        for (int i = 0; i < data.length; i++) {
            if (data[i] > 0) {
                gd.setAlpha(data[i]);
                // The sample will amplify the signal so we remap to the original scale
                data[i] = gd.sample() / emGain;
            }
        }
    }
    // Read-noise
    final NormalizedGaussianSampler gs = SamplerUtils.createNormalizedGaussianSampler(rg);
    if (noise != null) {
        for (int i = 0; i < data.length; i++) {
            data[i] += gs.sample() * noise[i];
        }
    }
    // uk.ac.sussex.gdsc.core.ij.Utils.display("Spot", data, size, size);
    return data;
}

Example 2 with MarsagliaTsangGammaSampler

use of uk.ac.sussex.gdsc.core.utils.rng.MarsagliaTsangGammaSampler in project GDSC-SMLM by aherbert.

the class CameraModelAnalysis method simulatePoissonGammaGaussian.

/**
 * Randomly generate a histogram from poisson-gamma-gaussian samples.
 *
 * @param settings the settings
 * @param random the random
 * @return The histogram
 */
private static IntHistogram simulatePoissonGammaGaussian(CameraModelAnalysisSettings settings, UniformRandomProvider random) {
    final int[] poissonSample = simulatePoisson(settings, random);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
    final double gain = getGain(settings);
    // Note that applying a separate EM-gain and then the camera gain later
    // is the same as applying the total gain in the gamma distribution and no camera gain
    // later, i.e. the Gamma distribution is just squeezed.
    // Thus we sample from a Gamma distribution with a fixed gain (the scale) and the
    // number of photons is the shape.
    final double noise = getReadNoise(settings);
    final int samples = settings.getSamples();
    final int emSamples = settings.getEmSamples();
    final int noiseSamples = (noise > 0) ? settings.getNoiseSamples() : 1;
    final int[] sample = new int[samples * emSamples * noiseSamples];
    int count = 0;
    final Round round = getRound(settings);
    if (gain < 1) {
        throw new IllegalStateException("EM-gain must be >= 1: " + gain);
    }
    final MarsagliaTsangGammaSampler gamma = new MarsagliaTsangGammaSampler(random, 1, gain);
    for (int n = poissonSample.length; n-- > 0; ) {
        if (poissonSample[n] != 0) {
            // Gamma
            gamma.setAlpha(poissonSample[n]);
            // Over-sample the Gamma
            for (int k = emSamples; k-- > 0; ) {
                // Do rounding to simulate a discrete PMF.
                final double d2 = round.round(gamma.sample());
                // Over-sample the Gaussian
                for (int j = noiseSamples; j-- > 0; ) {
                    // Convert the sample to a count
                    sample[count++] = round.round(d2 + noise * gauss.sample());
                }
            }
        } else {
            // Still over-sample the Gamma even though it was zero
            for (int k = emSamples; k-- > 0; ) {
                // Over-sample the Gaussian
                for (int j = noiseSamples; j-- > 0; ) {
                    // Convert the sample to a count
                    sample[count++] = round.round(noise * gauss.sample());
                }
            }
        }
    }
    return createHistogram(sample, count);
}

Example 3 with MarsagliaTsangGammaSampler

use of uk.ac.sussex.gdsc.core.utils.rng.MarsagliaTsangGammaSampler in project GDSC-SMLM by aherbert.

the class EmGainAnalysis method simulateFromPoissonGammaGaussian.

/**
 * Randomly generate a histogram from poisson-gamma-gaussian samples.
 *
 * @return The histogram
 */
private int[] simulateFromPoissonGammaGaussian() {
    // Randomly sample
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final PoissonSampler poisson = new PoissonSampler(rng, settings.settingPhotons);
    final MarsagliaTsangGammaSampler gamma = new MarsagliaTsangGammaSampler(rng, settings.settingPhotons, settings.settingGain);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    final int steps = settings.simulationSize;
    final int[] samples = new int[steps];
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0) {
            IJ.showProgress(n, steps);
        }
        // Poisson
        double sample = poisson.sample();
        // Gamma
        if (sample > 0) {
            gamma.setAlpha(sample);
            sample = gamma.sample();
        }
        // Gaussian
        sample += settings.settingNoise * gauss.sample();
        // Convert the sample to a count
        samples[n] = (int) Math.round(sample + settings.settingBias);
    }
    final int max = MathUtils.max(samples);
    final int[] histogram = new int[max + 1];
    for (final int s : samples) {
        histogram[s]++;
    }
    return histogram;
}