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

Example 81 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project gatk by broadinstitute.

the class CoverageDropoutDetectorTest method getUnivariateGaussianTargets.

private Object[][] getUnivariateGaussianTargets(final double sigma) {
    Random rng = new Random(337);
    final RandomGenerator randomGenerator = RandomGeneratorFactory.createRandomGenerator(rng);
    NormalDistribution n = new NormalDistribution(randomGenerator, 1, sigma);
    final int numDataPoints = 10000;
    final int numEventPoints = 2000;
    final List<ReadCountRecord.SingleSampleRecord> targetList = new ArrayList<>();
    for (int i = 0; i < (numDataPoints - numEventPoints); i++) {
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2 * i, 101 + 2 * i)), n.sample()));
    }
    for (int i = (numDataPoints - numEventPoints); i < numDataPoints; i++) {
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2 * i, 101 + 2 * i)), 0.5 + n.sample()));
    }
    HashedListTargetCollection<ReadCountRecord.SingleSampleRecord> targets = new HashedListTargetCollection<>(targetList);
    List<ModeledSegment> segments = new ArrayList<>();
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 100, 16050), 8000, 1));
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 16100, 20200), 2000, 1.5));
    return new Object[][] { { targets, segments } };
}

Example 82 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project gatk-protected by broadinstitute.

the class AdaptiveMetropolisSampler method sample.

public double sample(final RandomGenerator rng, final Function<Double, Double> logPDF) {
    Utils.nonNull(rng);
    Utils.nonNull(logPDF);
    final AbstractRealDistribution normal = new NormalDistribution(rng, 0, 1);
    final double proposal = xCurrent + stepSize * normal.sample();
    final double acceptanceProbability = (proposal < lowerBound || upperBound < proposal) ? 0 : Math.min(1, Math.exp(logPDF.apply(proposal) - logPDF.apply(xCurrent)));
    //adjust stepSize larger/smaller to decrease/increase the acceptance rate
    final double correctionFactor = (acceptanceProbability - optimalAcceptanceRate) * adjustmentRate * (timeScale / (timeScale + iteration));
    stepSize *= Math.exp(correctionFactor);
    iteration++;
    return rng.nextDouble() < acceptanceProbability ? proposal : xCurrent;
}

Example 83 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project gatk-protected by broadinstitute.

the class HMM method generateHiddenStateChain.

default default List<S> generateHiddenStateChain(final List<T> positions) {
    final RandomGenerator rg = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED_FOR_CHAIN_GENERATION));
    final List<S> hiddenStates = hiddenStates();
    final List<S> result = new ArrayList<>(positions.size());
    final S initialState = GATKProtectedMathUtils.randomSelect(hiddenStates, s -> Math.exp(logPriorProbability(s, positions.get(0))), rg);
    result.add(initialState);
    IntStream.range(1, positions.size()).forEach(n -> result.add(GATKProtectedMathUtils.randomSelect(hiddenStates, s -> Math.exp(logTransitionProbability(result.get(n - 1), positions.get(n - 1), s, positions.get(n))), rg)));
    return result;
}

Example 84 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project gatk-protected by broadinstitute.

the class CoverageDropoutDetectorTest method getUnivariateGaussianTargets.

private Object[][] getUnivariateGaussianTargets(final double sigma) {
    Random rng = new Random(337);
    final RandomGenerator randomGenerator = RandomGeneratorFactory.createRandomGenerator(rng);
    NormalDistribution n = new NormalDistribution(randomGenerator, 1, sigma);
    final int numDataPoints = 10000;
    final int numEventPoints = 2000;
    final List<ReadCountRecord.SingleSampleRecord> targetList = new ArrayList<>();
    for (int i = 0; i < (numDataPoints - numEventPoints); i++) {
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2 * i, 101 + 2 * i)), n.sample()));
    }
    for (int i = (numDataPoints - numEventPoints); i < numDataPoints; i++) {
        targetList.add(new ReadCountRecord.SingleSampleRecord(new Target("arbitrary_name", new SimpleInterval("chr1", 100 + 2 * i, 101 + 2 * i)), 0.5 + n.sample()));
    }
    HashedListTargetCollection<ReadCountRecord.SingleSampleRecord> targets = new HashedListTargetCollection<>(targetList);
    List<ModeledSegment> segments = new ArrayList<>();
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 100, 16050), 8000, 1));
    segments.add(new ModeledSegment(new SimpleInterval("chr1", 16100, 20200), 2000, 1.5));
    return new Object[][] { { targets, segments } };
}

Example 85 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project gatk-protected by broadinstitute.

the class JointAFCRSegmenterUnitTest method testSegmentation.

@Test
public void testSegmentation() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
    // probability that a datum is a het i.e. #hets / (#hets + #targets)
    final double hetProportion = 0.25;
    final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
    final double[] trueMinorAlleleFractions = new double[] { 0.12, 0.32, 0.5 };
    final double[] trueLog2CopyRatios = new double[] { -2.0, 0.0, 1.7 };
    final List<AFCRHiddenState> trueJointStates = IntStream.range(0, trueLog2CopyRatios.length).mapToObj(n -> new AFCRHiddenState(trueMinorAlleleFractions[n], trueLog2CopyRatios[n])).collect(Collectors.toList());
    final double trueMemoryLength = 1e5;
    final double trueCauchyWidth = 0.2;
    final int initialNumCRStates = 20;
    final int initialNumAFStates = 20;
    final AlleleFractionGlobalParameters trueAFParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);
    final JointAFCRHMM trueJointModel = new JointAFCRHMM(trueJointStates, trueWeights, trueMemoryLength, trueAFParams, AllelicPanelOfNormals.EMPTY_PON, trueCauchyWidth);
    // generate joint truth
    final int chainLength = 10000;
    final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
    final List<Integer> trueHiddenStates = trueJointModel.generateHiddenStateChain(positions);
    final List<AFCRHiddenState> trueAFCRSequence = trueHiddenStates.stream().map(trueJointModel::getHiddenStateValue).collect(Collectors.toList());
    final double[] trueCopyRatioSequence = trueAFCRSequence.stream().mapToDouble(AFCRHiddenState::getLog2CopyRatio).toArray();
    final double[] trueAlleleFractionSequence = trueAFCRSequence.stream().mapToDouble(AFCRHiddenState::getMinorAlleleFraction).toArray();
    // generate separate af and cr data
    final GammaDistribution biasGenerator = AlleleFractionSegmenterUnitTest.getGammaDistribution(trueAFParams, rng);
    final double outlierProbability = trueAFParams.getOutlierProbability();
    final AllelicCountCollection afData = new AllelicCountCollection();
    final List<Double> crData = new ArrayList<>();
    final List<Target> crTargets = new ArrayList<>();
    for (int n = 0; n < positions.size(); n++) {
        final SimpleInterval position = positions.get(n);
        final AFCRHiddenState jointState = trueAFCRSequence.get(n);
        final double minorFraction = jointState.getMinorAlleleFraction();
        final double log2CopyRatio = jointState.getLog2CopyRatio();
        if (rng.nextDouble() < hetProportion) {
            // het datum
            afData.add(AlleleFractionSegmenterUnitTest.generateAllelicCount(minorFraction, position, rng, biasGenerator, outlierProbability));
        } else {
            //target datum
            crTargets.add(new Target(position));
            crData.add(CopyRatioSegmenterUnitTest.generateData(trueCauchyWidth, log2CopyRatio, rng));
        }
    }
    final ReadCountCollection rcc = new ReadCountCollection(crTargets, Arrays.asList("SAMPLE"), new Array2DRowRealMatrix(crData.stream().mapToDouble(x -> x).toArray()));
    final JointAFCRSegmenter segmenter = JointAFCRSegmenter.createJointSegmenter(initialNumCRStates, rcc, initialNumAFStates, afData, AllelicPanelOfNormals.EMPTY_PON);
    final TargetCollection<SimpleInterval> tc = new HashedListTargetCollection<>(positions);
    final List<Pair<SimpleInterval, AFCRHiddenState>> segmentation = segmenter.findSegments();
    final List<ACNVModeledSegment> jointSegments = segmentation.stream().map(pair -> {
        final SimpleInterval position = pair.getLeft();
        final AFCRHiddenState jointState = pair.getRight();
        final PosteriorSummary crSummary = PerformJointSegmentation.errorlessPosterior(jointState.getLog2CopyRatio());
        final PosteriorSummary afSummary = PerformJointSegmentation.errorlessPosterior(jointState.getMinorAlleleFraction());
        return new ACNVModeledSegment(position, crSummary, afSummary);
    }).collect(Collectors.toList());
    final double[] segmentCopyRatios = jointSegments.stream().flatMap(s -> Collections.nCopies(tc.targetCount(s.getInterval()), s.getSegmentMeanPosteriorSummary().getCenter()).stream()).mapToDouble(x -> x).toArray();
    final double[] segmentMinorFractions = jointSegments.stream().flatMap(s -> Collections.nCopies(tc.targetCount(s.getInterval()), s.getMinorAlleleFractionPosteriorSummary().getCenter()).stream()).mapToDouble(x -> x).toArray();
    final double averageMinorFractionError = Arrays.stream(MathArrays.ebeSubtract(trueAlleleFractionSequence, segmentMinorFractions)).map(Math::abs).average().getAsDouble();
    final double averageCopyRatioError = Arrays.stream(MathArrays.ebeSubtract(trueCopyRatioSequence, segmentCopyRatios)).map(Math::abs).average().getAsDouble();
    Assert.assertEquals(averageMinorFractionError, 0, 0.04);
    Assert.assertEquals(averageCopyRatioError, 0, 0.04);
}