Search in sources :

Example 16 with SharedStateContinuousSampler

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

the class CompositeSamplersTest method testSharedStateContinuousSamplerSamples.

/**
 * Test sampling is uniform across several SharedStateContinuousSampler samplers.
 */
@Test
void testSharedStateContinuousSamplerSamples() {
    final Builder<SharedStateContinuousSampler> builder = CompositeSamplers.newSharedStateContinuousSamplerBuilder();
    final UniformRandomProvider rng = RandomSource.PCG_RXS_M_XS_64_OS.create(0x567567345L);
    final int n = 11;
    final double min = -15.7;
    final double max = 123.4;
    addContinuousSamplers(builder, n, min, max, rng);
    assertContinuousSamplerSamples(builder.build(rng), min, max);
}
Also used : SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) ChiSquareTest(org.apache.commons.math3.stat.inference.ChiSquareTest) Test(org.junit.jupiter.api.Test)

Example 17 with SharedStateContinuousSampler

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

the class DynamicMultipleTargetTracingTest method testTraceMoleculesDisableLocalDiffusionModel.

/**
 * 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 testTraceMoleculesDisableLocalDiffusionModel() {
    // The test is not very robust and fails 20% of the time. A fixed seed corrects this.
    final UniformRandomProvider rng = RngUtils.create(0x12345L);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    // localisation precision (in pixels)
    final double s = 0.1;
    final SharedStateContinuousSampler intensity1 = SamplerUtils.createGaussianSampler(rng, 1000, 100);
    final SharedStateContinuousSampler intensity2 = SamplerUtils.createGaussianSampler(rng, 500, 50);
    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;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (i + gauss.sample() * s), (float) (i + gauss.sample() * s), (float) (intensity1.sample())));
    }
    // Second molecule is fixed in the centre with a lower intensity (allow
    // correct matching when tracks overlap)
    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) (intensity2.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, size, size, (float) (intensity1.sample())));
    // Move centre to centre each jump => sqrt(2 * 0.1^2) = 0.141 um or 0.02 um^2
    // MSD = 4D => D = 0.02 / 4 = 0.005
    final DmttConfiguration config = DmttConfiguration.newBuilder(0.005).setDisableLocalDiffusionModel(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));
    }
}
Also used : SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) DmttConfiguration(uk.ac.sussex.gdsc.smlm.results.DynamicMultipleTargetTracing.DmttConfiguration) CalibrationWriter(uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest) Test(org.junit.jupiter.api.Test)

Example 18 with SharedStateContinuousSampler

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

the class DynamicMultipleTargetTracingTest method testTraceMolecules.

/**
 * Test trace molecules using 2 molecules. One is fixed and the other moves across it. The tracing
 * should assign the fixed molecule correctly as it has a low local diffusion rate and different
 * intensity.
 */
@Test
void testTraceMolecules() {
    // The test is not very robust and fails 10% of the time. A fixed seed corrects this.
    final UniformRandomProvider rng = RngUtils.create(0x12345L);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    // localisation precision (in pixels)
    final double s = 0.1;
    final SharedStateContinuousSampler intensity1 = SamplerUtils.createGaussianSampler(rng, 1000, 100);
    final SharedStateContinuousSampler intensity2 = SamplerUtils.createGaussianSampler(rng, 500, 50);
    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;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (i + gauss.sample() * s), (float) (i + gauss.sample() * s), (float) (intensity1.sample())));
    }
    // Second molecule is fixed in the centre with a lower intensity (allow
    // correct matching when tracks overlap)
    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) (intensity2.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, size, size, (float) (intensity1.sample())));
    // 1 diffuses top-left to bottom-right.
    // 2 is fixed in the centre.
    // 3 is in the bottom-right for 1 frame.
    // 
    // 1
    // 1
    // 1
    // 12
    // 1
    // 1
    // 13
    // 
    // Molecule 3 can sometimes connect to the long lifetime molecules once they have been alive
    // long enough to create a local probability model. The default lifetime is 5 frames.
    // Setting this to 10 frames allows a better local model to be created.
    // Move centre to centre each jump => sqrt(2 * 0.1^2) = 0.141 um or 0.02 um^2
    // MSD = 4D => D = 0.02 / 4 = 0.005
    final DmttConfiguration config = DmttConfiguration.newBuilder(0.005).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));
    }
}
Also used : SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) DmttConfiguration(uk.ac.sussex.gdsc.smlm.results.DynamicMultipleTargetTracing.DmttConfiguration) CalibrationWriter(uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest) Test(org.junit.jupiter.api.Test)

Example 19 with SharedStateContinuousSampler

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

the class WeightedMeanFilterTest method filterPerformsWeightedMeanFiltering.

@SeededTest
void filterPerformsWeightedMeanFiltering(RandomSeed seed) {
    final DataFilter filter = createDataFilter();
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, 2, 0.2);
    final float[] offsets = getOffsets(filter);
    final int[] boxSizes = getBoxSizes(filter);
    final TDoubleArrayList l1 = new TDoubleArrayList();
    final FloatFloatBiPredicate equality = TestHelper.floatsAreClose(1e-6, 0);
    for (final int width : primes) {
        for (final int height : primes) {
            final float[] data = createData(width, height, rg);
            l1.reset();
            // Uniform weights
            final float[] w1 = new float[width * height];
            Arrays.fill(w1, 0.5f);
            // Weights simulating the variance of sCMOS pixels
            final float[] w2 = new float[width * height];
            for (int i = 0; i < w2.length; i++) {
                w2[i] = (float) (1.0 / Math.max(0.01, gs.sample()));
            }
            for (final int boxSize : boxSizes) {
                for (final float offset : offsets) {
                    for (final boolean internal : checkInternal) {
                        // For each pixel over the range around the pixel (vi).
                        // Mean filter: sum(vi) / n
                        // Weighted mean: sum(vi * wi) / sum(wi) == mean(vi * wi) / mean(wi)
                        filter.setWeights(null, width, height);
                        // Uniform weights
                        testfilterPerformsWeightedFiltering(filter, width, height, data, w1, boxSize, offset, internal, equality);
                        // Random weights.
                        testfilterPerformsWeightedFiltering(filter, width, height, data, w2, boxSize, offset, internal, equality);
                    }
                }
            }
        }
    }
}
Also used : TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) FloatFloatBiPredicate(uk.ac.sussex.gdsc.test.api.function.FloatFloatBiPredicate) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 20 with SharedStateContinuousSampler

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

the class PulseActivationAnalysis method showSimulationDialog.

private boolean showSimulationDialog() {
    final ExtendedGenericDialog gd = new ExtendedGenericDialog(title);
    final SimulationDistribution[] distributionValues = SimulationDistribution.values();
    final String[] distribution = SettingsManager.getNames((Object[]) distributionValues);
    // Random crosstalk if not set
    if (MathUtils.max(settings.ct) == 0) {
        // Have some crosstalk
        final SharedStateContinuousSampler sampler = ContinuousUniformSampler.of(getUniformRandomProvider(), 0.05, 0.15);
        for (int i = 0; i < settings.ct.length; i++) {
            settings.ct[i] = sampler.sample();
        }
    }
    // Three channel
    for (int c = 0; c < 3; c++) {
        final String ch = "_C" + (c + 1);
        gd.addNumericField("Molcules" + ch, settings.numberOfMolecules[c], 0);
        gd.addChoice("Distribution" + ch, distribution, distribution[settings.distribution[c].ordinal()]);
        gd.addNumericField("Precision_" + ch, settings.precision[c], 3);
        gd.addNumericField("Crosstalk_" + Settings.ctNames[2 * c], settings.ct[2 * c], 3);
        gd.addNumericField("Crosstalk_" + Settings.ctNames[2 * c + 1], settings.ct[2 * c + 1], 3);
    }
    gd.addHelp(HelpUrls.getUrl(helpKey));
    gd.showDialog();
    if (gd.wasCanceled()) {
        return false;
    }
    int count = 0;
    for (int c = 0; c < 3; c++) {
        settings.numberOfMolecules[c] = (int) Math.abs(gd.getNextNumber());
        if (settings.numberOfMolecules[c] > 0) {
            count++;
        }
        settings.distribution[c] = distributionValues[gd.getNextChoiceIndex()];
        settings.precision[c] = Math.abs(gd.getNextNumber());
        settings.ct[2 * c] = Math.abs(gd.getNextNumber());
        settings.ct[2 * c + 1] = Math.abs(gd.getNextNumber());
    }
    settings.save();
    if (gd.invalidNumber()) {
        return false;
    }
    if (count < 2) {
        IJ.error(title, "Simulation requires at least 2 channels");
        return false;
    }
    try {
        for (int i = 0; i < settings.ct.length; i += 2) {
            if (settings.numberOfMolecules[i / 2] > 0) {
                validateCrosstalk(i, i + 1);
            }
        }
    } catch (final IllegalArgumentException ex) {
        IJ.error(title, ex.getMessage());
        return false;
    }
    return true;
}
Also used : SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) NonBlockingExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)

Aggregations

SharedStateContinuousSampler (org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler)21 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)18 SeededTest (uk.ac.sussex.gdsc.test.junit5.SeededTest)9 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)5 Test (org.junit.jupiter.api.Test)4 ContinuousSampler (org.apache.commons.rng.sampling.distribution.ContinuousSampler)3 NormalizedGaussianSampler (org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler)3 CalibrationWriter (uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter)3 DmttConfiguration (uk.ac.sussex.gdsc.smlm.results.DynamicMultipleTargetTracing.DmttConfiguration)3 FloatFloatBiPredicate (uk.ac.sussex.gdsc.test.api.function.FloatFloatBiPredicate)3 FloatProcessor (ij.process.FloatProcessor)2 ImageProcessor (ij.process.ImageProcessor)2 Rectangle (java.awt.Rectangle)2 DiscreteSampler (org.apache.commons.rng.sampling.distribution.DiscreteSampler)2 ImageJImagePeakResults (uk.ac.sussex.gdsc.smlm.ij.results.ImageJImagePeakResults)2 ImagePlus (ij.ImagePlus)1 ImageStack (ij.ImageStack)1 File (java.io.File)1 BigDecimal (java.math.BigDecimal)1 MathContext (java.math.MathContext)1