Examples with Well44497b org.apache.commons.math3.random.Well44497b used on opensource projects

Example 6 with Well44497b

use of org.apache.commons.math3.random.Well44497b 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
    RandomGenerator random = new Well44497b(System.currentTimeMillis() + System.identityHashCode(this));
    PoissonDistribution poisson = new PoissonDistribution(random, _photons, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    CustomGammaDistribution gamma = new CustomGammaDistribution(random, _photons, _gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    final int steps = simulationSize;
    int[] sample = new int[steps];
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0)
            IJ.showProgress(n, steps);
        // Poisson
        double d = poisson.sample();
        // Gamma
        if (d > 0) {
            gamma.setShapeUnsafe(d);
            d = gamma.sample();
        }
        // Gaussian
        d += _noise * random.nextGaussian();
        // Convert the sample to a count 
        sample[n] = (int) Math.round(d + _bias);
    }
    int max = Maths.max(sample);
    int[] h = new int[max + 1];
    for (int s : sample) h[s]++;
    return h;
}

Example 7 with Well44497b

use of org.apache.commons.math3.random.Well44497b in project GDSC-SMLM by aherbert.

the class EMGainAnalysis method simulateFromPDF.

/**
	 * Random sample from the cumulative probability distribution function that is used during fitting
	 * 
	 * @return The histogram
	 */
private int[] simulateFromPDF() {
    final double[] g = pdf(0, _photons, _gain, _noise, (int) _bias);
    // Debug this
    double[] x = Utils.newArray(g.length, 0, 1.0);
    Utils.display(TITLE + " PDF", new Plot(TITLE + " PDF", "ADU", "P", x, Arrays.copyOf(g, g.length)));
    // Get cumulative probability
    double sum = 0;
    for (int i = 0; i < g.length; i++) {
        final double p = g[i];
        g[i] += sum;
        sum += p;
    }
    for (int i = 0; i < g.length; i++) g[i] /= sum;
    // Ensure value of 1 at the end
    g[g.length - 1] = 1;
    // Randomly sample
    RandomGenerator random = new Well44497b(System.currentTimeMillis() + System.identityHashCode(this));
    int[] h = new int[g.length];
    final int steps = simulationSize;
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0)
            IJ.showProgress(n, steps);
        final double p = random.nextDouble();
        for (int i = 0; i < g.length; i++) if (p <= g[i]) {
            h[i]++;
            break;
        }
    }
    return h;
}

Example 8 with Well44497b

use of org.apache.commons.math3.random.Well44497b in project mixcr by milaboratory.

the class FullSeqAssemblerTest method testLargeCloneNoMismatches.

@Test
public void testLargeCloneNoMismatches() throws Exception {
    MasterSequence master = FullSeqAssemblerTest.masterSeq1WT;
    NSequenceWithQuality seq = new NSequenceWithQuality(master.getRange(-master.vPart + 10, 80), SequenceQuality.GOOD_QUALITY_VALUE);
    RunMiXCR.RunMiXCRAnalysis params0 = new RunMiXCR.RunMiXCRAnalysis(new SingleReadImpl(0, seq, ""));
    params0.cloneAssemblerParameters.setAssemblingFeatures(new GeneFeature[] { GeneFeature.VDJRegion });
    Clone largeClone = RunMiXCR.assemble(RunMiXCR.align(params0)).cloneSet.get(0);
    // ActionExportClonesPretty.outputCompact(System.out, largeClone);
    // System.exit(0);
    Well44497b rnd = new Well44497b(1234567889L);
    int nReads = 100_000;
    int readLength = 75, readGap = 150;
    // slice seq randomly
    PairedRead[] slicedReads = new PairedRead[nReads];
    for (int i = 0; i < nReads; ++i) {
        int r1from = rnd.nextInt(seq.size() - readLength - 1), r1to = r1from + readLength, r2from = r1from + 1 + rnd.nextInt(seq.size() - r1from - readLength - 1), r2to = r2from + readLength;
        assert r2from > r1from;
        slicedReads[i] = new PairedRead(new SingleReadImpl(i, seq.getRange(r1from, r1to), "" + i), new SingleReadImpl(i, seq.getRange(r2from, r2to).getReverseComplement(), "" + i));
    }
    RunMiXCR.RunMiXCRAnalysis params = new RunMiXCR.RunMiXCRAnalysis(slicedReads);
    // params.alignerParameters = VDJCParametersPresets.getByName("rna-seq");
    params.alignerParameters.setSaveOriginalReads(true);
    RunMiXCR.AlignResult align = RunMiXCR.align(params);
    RunMiXCR.AssembleResult assemble = RunMiXCR.assemble(align);
    // for (VDJCAlignments al : align.alignments) {
    // for (int i = 0; i < al.numberOfTargets(); i++) {
    // System.out.println(VDJCAlignmentsFormatter.getTargetAsMultiAlignment(al, i));
    // System.out.println();
    // }
    // System.out.println();
    // System.out.println(" ================================================ ");
    // System.out.println();
    // }
    Assert.assertEquals(1, assemble.cloneSet.size());
    Clone initialClone = assemble.cloneSet.get(0);
    NSequenceWithQuality cdr3 = initialClone.getFeature(GeneFeature.CDR3);
    List<VDJCAlignments> alignments = align.alignments.stream().filter(al -> cdr3.equals(al.getFeature(GeneFeature.CDR3))).collect(Collectors.toList());
    alignments.stream().filter(al -> Arrays.stream(al.getBestHit(GeneType.Variable).getAlignments()).filter(Objects::nonNull).anyMatch(a -> !a.getAbsoluteMutations().isEmpty())).filter(al -> al.getBestHit(GeneType.Variable).getGene().getName().contains("3-74")).forEach(al -> {
        for (int i = 0; i < al.numberOfTargets(); i++) {
            System.out.println(VDJCAlignmentsFormatter.getTargetAsMultiAlignment(al, i));
            System.out.println();
        }
        System.out.println();
        System.out.println(" ================================================ ");
        System.out.println();
    });
    // System.exit(0);
    System.out.println("=> Agg");
    CloneFactory cloneFactory = new CloneFactory(align.parameters.cloneAssemblerParameters.getCloneFactoryParameters(), align.parameters.cloneAssemblerParameters.getAssemblingFeatures(), align.usedGenes, align.parameters.alignerParameters.getFeaturesToAlignMap());
    FullSeqAssembler agg = new FullSeqAssembler(cloneFactory, DEFAULT_PARAMETERS, initialClone, align.parameters.alignerParameters);
    FullSeqAssembler.RawVariantsData prep = agg.calculateRawData(() -> CUtils.asOutputPort(alignments));
    List<Clone> clones = new ArrayList<>(new CloneSet(Arrays.asList(agg.callVariants(prep))).getClones());
    clones.sort(Comparator.comparingDouble(Clone::getCount).reversed());
    for (Clone clone : clones) {
        ActionExportClonesPretty.outputCompact(System.out, clone);
        System.out.println();
        System.out.println(" ================================================ ");
        System.out.println();
    }
}

Example 9 with Well44497b

use of org.apache.commons.math3.random.Well44497b in project mixcr by milaboratory.

the class PartialAlignmentsAssemblerAlignerTest method basicTest1.

@Test
public void basicTest1() throws Exception {
    Well44497b rg = new Well44497b(12312);
    VDJCAlignerParameters rnaSeqParams = VDJCParametersPresets.getByName("rna-seq");
    PartialAlignmentsAssemblerAligner aligner = new PartialAlignmentsAssemblerAligner(rnaSeqParams);
    VDJCLibrary lib = VDJCLibraryRegistry.getDefault().getLibrary("default", "hs");
    for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").getGenes()) if (gene.isFunctional())
        aligner.addGene(gene);
    TargetBuilder.VDJCGenes genes = new TargetBuilder.VDJCGenes(lib, "TRBV12-3*00", "TRBD1*00", "TRBJ1-3*00", "TRBC2*00");
    // | 305
    // 250V + 55CDR3 (20V 7N 10D 3N 15J) + 28J + 100C
    NucleotideSequence baseSeq = TargetBuilder.generateSequence(genes, "{CDR3Begin(-250)}V*270 NNNNNNN {DBegin(0)}D*10 NNN {CDR3End(-15):FR4End} {CBegin}C*100", rg);
    int len = 70;
    NucleotideSequence seq1 = baseSeq.getRange(0, len);
    NucleotideSequence seq2 = baseSeq.getRange(245, 245 + len);
    NucleotideSequence seq3 = baseSeq.getRange(320, 320 + len);
    VDJCAlignmentResult<VDJCMultiRead> alignment = aligner.process(MiXCRTestUtils.createMultiRead(seq1, seq2, seq3));
    VDJCAlignments al = alignment.alignment;
    Assert.assertNotNull(al);
    assertInHits(genes.v, al);
    assertInHits(genes.d, al);
    assertInHits(genes.j, al);
    assertInHits(genes.c, al);
    VDJCHit bestV = al.getBestHit(GeneType.Variable);
    VDJCHit bestD = al.getBestHit(GeneType.Diversity);
    VDJCHit bestJ = al.getBestHit(GeneType.Joining);
    VDJCHit bestC = al.getBestHit(GeneType.Constant);
    Assert.assertNotNull(bestV.getAlignment(0));
    Assert.assertNotNull(bestV.getAlignment(1));
    Assert.assertNull(bestV.getAlignment(2));
    Assert.assertNull(bestD.getAlignment(0));
    Assert.assertNotNull(bestD.getAlignment(1));
    Assert.assertNull(bestD.getAlignment(2));
    Assert.assertNull(bestJ.getAlignment(0));
    Assert.assertNotNull(bestJ.getAlignment(1));
    Assert.assertNotNull(bestJ.getAlignment(2));
    Assert.assertNull(bestC.getAlignment(0));
    Assert.assertNull(bestC.getAlignment(1));
    Assert.assertNotNull(bestC.getAlignment(2));
}

Example 10 with Well44497b

use of org.apache.commons.math3.random.Well44497b in project tetrad by cmu-phil.

the class RandomUtil method setSeed.

/**
 * Sets the seed to the given value.
 *
 * @param seed A long value. Once this seed is set, the behavior of the random number generator is deterministic, so
 *             setting the seed can be used to repeat previous behavior.
 */
public void setSeed(long seed) {
    // Do not change this generator; you will screw up innuerable unit tests!
    randomGenerator = new SynchronizedRandomGenerator(new Well44497b(seed));
    seedsToGenerators.put(seed, randomGenerator);
    normal = new NormalDistribution(randomGenerator, 0, 1);
    this.seed = seed;
}