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

Example 1 with NormalizedGaussianSampler

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

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

the class TraceExporter method export.

private void export(MemoryPeakResults results) {
    // Copy to allow manipulation
    results = results.copy();
    // Strip results with no trace Id
    results.removeIf(result -> result.getId() <= 0);
    // Sort by ID then time
    results.sort(IdFramePeakResultComparator.INSTANCE);
    // Split traces with big jumps and long tracks into smaller tracks
    results = splitTraces(results);
    results = splitLongTraces(results);
    // Count each ID and remove short traces
    int id = 0;
    int count = 0;
    int tracks = 0;
    int maxLength = 0;
    final TIntHashSet remove = new TIntHashSet();
    for (int i = 0, size = results.size(); i < size; i++) {
        final PeakResult result = results.get(i);
        if (result.getId() != id) {
            if (count < settings.minLength) {
                remove.add(id);
            } else {
                tracks++;
                maxLength = Math.max(maxLength, count);
            }
            count = 0;
            id = result.getId();
        }
        count += getLength(result);
    }
    // Final ID
    if (count < settings.minLength) {
        remove.add(id);
    } else {
        tracks++;
        maxLength = Math.max(maxLength, count);
    }
    if (!remove.isEmpty()) {
        results.removeIf(result -> remove.contains(result.getId()));
        results.sort(IdFramePeakResultComparator.INSTANCE);
    }
    if (settings.wobble > 0) {
        // Just leave any exceptions to trickle up and kill the plugin
        final TypeConverter<DistanceUnit> c = results.getDistanceConverter(DistanceUnit.NM);
        final double w = c.convertBack(settings.wobble);
        final UniformRandomProvider rng = UniformRandomProviders.create();
        final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
        final boolean is3D = results.is3D();
        results.forEach((PeakResultProcedure) peakResult -> {
            peakResult.setXPosition((float) (peakResult.getXPosition() + w * gauss.sample()));
            peakResult.setYPosition((float) (peakResult.getYPosition() + w * gauss.sample()));
            if (is3D) {
                peakResult.setZPosition((float) (peakResult.getZPosition() + w * gauss.sample()));
            }
        });
    }
    if (settings.format == 2) {
        exportVbSpt(results);
    } else if (settings.format == 1) {
        exportAnaDda(results);
    } else {
        exportSpotOn(results);
    }
    ImageJUtils.log("Exported %s: %s in %s", results.getName(), TextUtils.pleural(results.size(), "localisation"), TextUtils.pleural(tracks, "track"));
    if (settings.showTraceLengths) {
        // We store and index (count-1)
        final int[] h = new int[maxLength];
        id = 0;
        for (int i = 0, size = results.size(); i < size; i++) {
            final PeakResult result = results.get(i);
            if (result.getId() != id) {
                h[count - 1]++;
                count = 0;
                id = result.getId();
            }
            count += getLength(result);
        }
        h[count - 1]++;
        final String title = TITLE + ": " + results.getName();
        final Plot plot = new Plot(title, "Length", "Frequency");
        plot.addPoints(SimpleArrayUtils.newArray(h.length, 1, 1.0f), SimpleArrayUtils.toFloat(h), Plot.BAR);
        plot.setLimits(SimpleArrayUtils.findIndex(h, i -> i != 0), maxLength + 1, 0, Double.NaN);
        ImageJUtils.display(title, plot);
    }
}

Example 3 with NormalizedGaussianSampler

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

the class PulseActivationAnalysis method simulateActivations.

private int simulateActivations(UniformRandomProvider rng, DiscreteSampler bd, float[][] molecules, MemoryPeakResults results, int time, double precision) {
    if (bd == null) {
        return 0;
    }
    final int size = molecules.length;
    final int numberOfSamples = bd.sample();
    // Sample
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    final int[] samples = RandomUtils.sample(numberOfSamples, size, rng);
    for (final int index : samples) {
        final float[] xy = molecules[index];
        float x;
        float y;
        do {
            x = (float) (xy[0] + gauss.sample() * precision);
        } while (outOfBounds(x));
        do {
            y = (float) (xy[1] + gauss.sample() * precision);
        } while (outOfBounds(y));
        results.add(createResult(time, x, y));
    }
    return samples.length;
}

Example 4 with NormalizedGaussianSampler

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

the class DiffusionRateTest method simpleTest.

/**
 * Perform a simple diffusion test. This can be used to understand the distributions that are
 * generated during 3D diffusion.
 */
private void simpleTest() {
    if (!showSimpleDialog()) {
        return;
    }
    final StoredDataStatistics[] stats2 = new StoredDataStatistics[3];
    final StoredDataStatistics[] stats = new StoredDataStatistics[3];
    final NormalizedGaussianSampler[] gauss = new NormalizedGaussianSampler[3];
    for (int i = 0; i < 3; i++) {
        stats2[i] = new StoredDataStatistics(pluginSettings.simpleParticles);
        stats[i] = new StoredDataStatistics(pluginSettings.simpleParticles);
        gauss[i] = SamplerUtils.createNormalizedGaussianSampler(UniformRandomProviders.create());
    }
    final double scale = Math.sqrt(2 * pluginSettings.simpleD);
    final int report = Math.max(1, pluginSettings.simpleParticles / 200);
    for (int particle = 0; particle < pluginSettings.simpleParticles; particle++) {
        if (particle % report == 0) {
            IJ.showProgress(particle, pluginSettings.simpleParticles);
        }
        final double[] xyz = new double[3];
        if (pluginSettings.linearDiffusion) {
            final double[] dir = nextVector(gauss[0]);
            for (int step = 0; step < pluginSettings.simpleSteps; step++) {
                final double d = gauss[0].sample();
                for (int i = 0; i < 3; i++) {
                    xyz[i] += dir[i] * d;
                }
            }
        } else {
            for (int step = 0; step < pluginSettings.simpleSteps; step++) {
                for (int i = 0; i < 3; i++) {
                    xyz[i] += gauss[i].sample();
                }
            }
        }
        for (int i = 0; i < 3; i++) {
            xyz[i] *= scale;
        }
        double msd = 0;
        for (int i = 0; i < 3; i++) {
            msd += xyz[i] * xyz[i];
            stats2[i].add(msd);
            // Store the actual distances
            stats[i].add(xyz[i]);
        }
    }
    IJ.showProgress(1);
    for (int i = 0; i < 3; i++) {
        plotJumpDistances(TITLE, stats2[i], i + 1);
        // Save stats to file for fitting
        save(stats2[i], i + 1, "msd");
        save(stats[i], i + 1, "d");
    }
    windowOrganiser.tile();
}

Example 5 with NormalizedGaussianSampler

use of org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler 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);
}