Examples with PoissonDistribution org.apache.commons.math3.distribution.PoissonDistribution used on opensource projects

Example 1 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class CreateData method createBackground.

private float[] createBackground(RandomDataGenerator random) {
    float[] pixels2 = null;
    if (settings.background > 0) {
        if (random == null)
            random = new RandomDataGenerator(createRandomGenerator());
        createBackgroundPixels();
        pixels2 = Arrays.copyOf(backgroundPixels, backgroundPixels.length);
        // Add Poisson noise.
        if (uniformBackground) {
            // We can do N random samples thus ensuring the background average is constant.
            // Note: The number of samples must be Poisson distributed.
            // currently: pixels2[0] = uniform background level
            // => (pixels2[0] * pixels2.length) = amount of photons that fall on the image. 
            final int samples = (int) random.nextPoisson(pixels2[0] * pixels2.length);
            // Only do sampling if the number of samples is valid
            if (samples >= 1) {
                pixels2 = new float[pixels2.length];
                final int upper = pixels2.length - 1;
                for (int i = 0; i < samples; i++) {
                    pixels2[random.nextInt(0, upper)] += 1;
                }
            } else {
                // If using a uniform illumination then we can use a fixed Poisson distribution
                PoissonDistribution dist = new PoissonDistribution(random.getRandomGenerator(), pixels2[0], PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
                for (int i = 0; i < pixels2.length; i++) {
                    pixels2[i] = dist.sample();
                }
            }
        } else {
            for (int i = 0; i < pixels2.length; i++) {
                pixels2[i] = random.nextPoisson(pixels2[i]);
            }
        }
    } else {
        pixels2 = new float[settings.size * settings.size];
    }
    return pixels2;
}

Example 2 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class CMOSAnalysis method simulate.

private void simulate() {
    // Create the offset, variance and gain for each pixel
    int n = size * size;
    float[] pixelOffset = new float[n];
    float[] pixelVariance = new float[n];
    float[] pixelGain = new float[n];
    IJ.showStatus("Creating random per-pixel readout");
    long start = System.currentTimeMillis();
    RandomGenerator rg = new Well19937c();
    PoissonDistribution pd = new PoissonDistribution(rg, offset, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    ExponentialDistribution ed = new ExponentialDistribution(rg, variance, ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    totalProgress = n;
    stepProgress = Utils.getProgressInterval(totalProgress);
    for (int i = 0; i < n; i++) {
        if (i % n == 0)
            IJ.showProgress(i, n);
        // Q. Should these be clipped to a sensible range?
        pixelOffset[i] = (float) pd.sample();
        pixelVariance[i] = (float) ed.sample();
        pixelGain[i] = (float) (gain + rg.nextGaussian() * gainSD);
    }
    IJ.showProgress(1);
    // Avoid all the file saves from updating the progress bar and status line
    Utils.setShowStatus(false);
    Utils.setShowProgress(false);
    JLabel statusLine = Utils.getStatusLine();
    progressBar = Utils.getProgressBar();
    // Save to the directory as a stack
    ImageStack simulationStack = new ImageStack(size, 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", Utils.rounded(offset), Utils.rounded(variance), Utils.rounded(gain), Utils.rounded(gainSD)));
    IJ.save(simulationImp, new File(directory, "perPixelSimulation.tif").getPath());
    // Create thread pool and workers
    ExecutorService executor = Executors.newFixedThreadPool(getThreads());
    TurboList<Future<?>> futures = new TurboList<Future<?>>(nThreads);
    // Simulate the zero exposure input.
    // Simulate 20 - 200 photon images.
    int[] photons = new int[] { 0, 20, 50, 100, 200 };
    totalProgress = photons.length * frames;
    stepProgress = Utils.getProgressInterval(totalProgress);
    progress = 0;
    progressBar.show(0);
    // For saving stacks
    int blockSize = 10;
    int nPerThread = (int) Math.ceil((double) frames / nThreads);
    // Convert to fit the block size
    nPerThread = (int) Math.ceil((double) nPerThread / blockSize) * blockSize;
    long seed = start;
    for (int p : photons) {
        statusLine.setText("Simulating " + Utils.pleural(p, "photon"));
        // Create the directory
        File out = new File(directory, String.format("photon%03d", p));
        if (!out.exists())
            out.mkdir();
        for (int from = 0; from < frames; ) {
            int to = Math.min(from + nPerThread, frames);
            futures.add(executor.submit(new SimulationWorker(seed++, out.getPath(), simulationStack, from, to, blockSize, p)));
            from = to;
        }
        // Wait for all to finish
        for (int t = futures.size(); t-- > 0; ) {
            try {
                // The future .get() method will block until completed
                futures.get(t).get();
            } catch (Exception e) {
                // This should not happen. 
                e.printStackTrace();
            }
        }
        futures.clear();
    }
    Utils.setShowStatus(true);
    Utils.setShowProgress(true);
    IJ.showProgress(1);
    executor.shutdown();
    Utils.log("Simulation time = " + Utils.timeToString(System.currentTimeMillis() - start));
}

Example 3 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class PoissonCalculatorTest method canComputeLikelihoodForIntegerData.

@Test
public void canComputeLikelihoodForIntegerData() {
    for (double u : photons) {
        PoissonDistribution pd = new PoissonDistribution(u);
        for (int x = 0; x < 100; x++) {
            double e = pd.probability(x);
            double o = PoissonCalculator.likelihood(u, x);
            if (e > 1e-100)
                Assert.assertEquals(e, o, e * 1e-10);
            e = pd.logProbability(x);
            o = PoissonCalculator.logLikelihood(u, x);
            Assert.assertEquals(e, o, Math.abs(e) * 1e-10);
        }
    }
}

Example 4 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project gatk-protected by broadinstitute.

the class TargetCoverageSexGenotypeCalculator method calculateSexGenotypeData.

/**
     * Calculates the likelihood of different sex genotypes for a given sample index
     *
     * @param sampleIndex sample index
     * @return an instance of {@link SexGenotypeData}
     */
private SexGenotypeData calculateSexGenotypeData(final int sampleIndex) {
    final int[] allosomalReadCounts = Arrays.stream(processedReadCounts.getColumnOnSpecifiedTargets(sampleIndex, allosomalTargetList, true)).mapToInt(n -> (int) n).toArray();
    final double readDepth = getSampleReadDepthFromAutosomalTargets(sampleIndex);
    logger.debug("Read depth for " + processedReadCounts.columnNames().get(sampleIndex) + ": " + readDepth);
    final List<Double> logLikelihoods = new ArrayList<>();
    final List<String> sexGenotypesList = new ArrayList<>();
    final int numAllosomalTargets = allosomalTargetList.size();
    for (final String genotypeName : sexGenotypeIdentifiers) {
        /* calculate log likelihood */
        final int[] currentAllosomalTargetPloidies = allosomalTargetPloidies.get(genotypeName);
        final double[] poissonParameters = IntStream.range(0, numAllosomalTargets).mapToDouble(i -> readDepth * (currentAllosomalTargetPloidies[i] > 0 ? currentAllosomalTargetPloidies[i] : baselineMappingErrorProbability)).toArray();
        final double currentLogLikelihood = IntStream.range(0, numAllosomalTargets).mapToDouble(i -> {
            final PoissonDistribution pois = new PoissonDistribution(poissonParameters[i]);
            return pois.logProbability(allosomalReadCounts[i]);
        }).sum();
        sexGenotypesList.add(genotypeName);
        logLikelihoods.add(currentLogLikelihood);
    }
    /* infer the most likely sex genotype */
    final Integer[] indices = new Integer[sexGenotypesList.size()];
    IntStream.range(0, sexGenotypesList.size()).forEach(i -> indices[i] = i);
    Arrays.sort(indices, (li, ri) -> Double.compare(logLikelihoods.get(ri), logLikelihoods.get(li)));
    return new SexGenotypeData(processedReadCounts.columnNames().get(sampleIndex), sexGenotypesList.get(indices[0]), sexGenotypesList, logLikelihoods.stream().mapToDouble(d -> d).toArray());
}

Example 5 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project gatk-protected by broadinstitute.

the class CoverageModelCopyRatioEmissionProbabilityCalculator method logLikelihoodPoisson.

/**
     * Calculate emission log probability directly using Poisson distribution
     *
     * TODO github/gatk-protected issue #854
     *
     * @implNote This is a naive implementation where the variance of log bias ($\Psi$) is
     * ignored altogether. This routine must be improved by performing a few-point numerical
     * integration of:
     *
     *      \int_{-\infty}^{+\infty} db Poisson(n | \lambda = d*c*exp(b) + eps*d)
     *                                                           * Normal(b | \mu, \Psi)
     *
     * @param emissionData an instance of {@link CoverageModelCopyRatioEmissionData}
     * @param copyRatio copy ratio on which the emission probability is conditioned on
     * @return a double value
     */
private double logLikelihoodPoisson(@Nonnull CoverageModelCopyRatioEmissionData emissionData, double copyRatio) {
    final double multBias = FastMath.exp(emissionData.getMu());
    final double readDepth = emissionData.getCopyRatioCallingMetadata().getSampleCoverageDepth();
    final double err = emissionData.getMappingErrorProbability();
    final double poissonMean = readDepth * ((1 - err) * copyRatio * multBias + err);
    return new PoissonDistribution(null, poissonMean, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS).logProbability(emissionData.getReadCount());
}