Search in sources :

Example 1 with DefaultRealMatrixPreservingVisitor

use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk-protected by broadinstitute.

the class NormalizeSomaticReadCountsIntegrationTest method assertBetaHatsRobustToOutliers.

/**
     * Asserts that the calculation of beta hats is not significantly affected by zero-coverage outlier counts
     * We perform this check by randomly setting some coverages to zero in copy ratio space (-infinity in log space).
     * betaHats imputes 0 in log space (1 in copy ratio space) whenever coverage is below a certain low threshold
     * and should thus be robust to this type of noise.
     */
private void assertBetaHatsRobustToOutliers(final ReadCountCollection preTangentNormalized, final File ponFile) {
    try (final HDF5File ponReader = new HDF5File(ponFile)) {
        final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
        final List<String> ponTargets = pon.getPanelTargetNames();
        final RealMatrix input = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);
        // randomly set some entries to zero in copy-ratio space (-infinity in log space)
        final Random random = new Random(13);
        final double noiseProportion = 0.01;
        final RealMatrix noisyInput = input.copy();
        noisyInput.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {

            @Override
            public double visit(final int row, final int column, final double value) {
                return random.nextDouble() < noiseProportion ? Double.NEGATIVE_INFINITY : value;
            }
        });
        final RealMatrix betaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), input, PCATangentNormalizationUtils.EPSILON);
        final RealMatrix noisyBetaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), noisyInput, PCATangentNormalizationUtils.EPSILON);
        final RealMatrix difference = betaHats.subtract(noisyBetaHats);
        difference.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {

            @Override
            public void visit(final int row, int column, double value) {
                Assert.assertEquals(value, 0, 0.01);
            }
        });
    }
}
Also used : HDF5PCACoveragePoN(org.broadinstitute.hellbender.tools.pon.coverage.pca.HDF5PCACoveragePoN) DefaultRealMatrixPreservingVisitor(org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) RealMatrix(org.apache.commons.math3.linear.RealMatrix) DefaultRealMatrixChangingVisitor(org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor) PCACoveragePoN(org.broadinstitute.hellbender.tools.pon.coverage.pca.PCACoveragePoN) HDF5PCACoveragePoN(org.broadinstitute.hellbender.tools.pon.coverage.pca.HDF5PCACoveragePoN) HDF5File(org.broadinstitute.hdf5.HDF5File)

Example 2 with DefaultRealMatrixPreservingVisitor

use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk by broadinstitute.

the class CopyNumberTriStateTransitionProbabilityCacheUnitTest method markovianPropertiesTest.

//test various properties of a transition matrix
@Test(dataProvider = "meanEventSizeAndEventStartProbability")
public void markovianPropertiesTest(final double meanEventSize, final double eventStartProbability) {
    final CopyNumberTriStateTransitionProbabilityCache cache = new CopyNumberTriStateTransitionProbabilityCache(meanEventSize, eventStartProbability);
    for (final int d : DISTANCES) {
        //check symmetries -- these are part of the model and need not be true in the future
        final RealMatrix transitionMatrix = cache.getAsMatrixInProbabilitySpace(d);
        assertSymmetries(transitionMatrix);
        //check that columns sums equal 1
        for (int column = 0; column < transitionMatrix.getColumnDimension(); column++) {
            Assert.assertEquals(MathUtils.sum(transitionMatrix.getColumn(column)), 1, EPSILON);
        }
        //check that all elements are positive
        transitionMatrix.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {

            @Override
            public void visit(int row, int column, double value) {
                Assert.assertTrue(value >= 0);
            }
        });
        //check that T(2d) = T(d)*T(d)
        assertEqualMatrices(cache.getAsMatrixInProbabilitySpace(2 * d), transitionMatrix.multiply(transitionMatrix));
        //check that the largest eigenvalue of the transition matrix is 1 (this corresponds to the asymptotic stationary state)
        Assert.assertEquals(MathUtils.arrayMax(new EigenDecomposition(transitionMatrix).getRealEigenvalues()), 1, EPSILON);
    }
    // check that at long distances memory of the initial state is lost and all initial distributions tend toward
    // the same asymptotic stationary distribution.  That is, all columns of the large-distance transition matrix are equal
    final RealMatrix asymptoticMatrix = cache.getAsMatrixInProbabilitySpace(HUGE_DISTANCE);
    for (int column = 1; column < asymptoticMatrix.getColumnDimension(); column++) {
        Assert.assertEquals(asymptoticMatrix.getColumnVector(0).subtract(asymptoticMatrix.getColumnVector(column)).getL1Norm(), 0, EPSILON);
    }
}
Also used : DefaultRealMatrixPreservingVisitor(org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor) EigenDecomposition(org.apache.commons.math3.linear.EigenDecomposition) RealMatrix(org.apache.commons.math3.linear.RealMatrix) Test(org.testng.annotations.Test)

Example 3 with DefaultRealMatrixPreservingVisitor

use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk by broadinstitute.

the class NormalizeSomaticReadCountsIntegrationTest method assertBetaHatsRobustToOutliers.

/**
     * Asserts that the calculation of beta hats is not significantly affected by zero-coverage outlier counts
     * We perform this check by randomly setting some coverages to zero in copy ratio space (-infinity in log space).
     * betaHats imputes 0 in log space (1 in copy ratio space) whenever coverage is below a certain low threshold
     * and should thus be robust to this type of noise.
     */
private void assertBetaHatsRobustToOutliers(final ReadCountCollection preTangentNormalized, final File ponFile) {
    try (final HDF5File ponReader = new HDF5File(ponFile)) {
        final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
        final List<String> ponTargets = pon.getPanelTargetNames();
        final RealMatrix input = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);
        // randomly set some entries to zero in copy-ratio space (-infinity in log space)
        final Random random = new Random(13);
        final double noiseProportion = 0.01;
        final RealMatrix noisyInput = input.copy();
        noisyInput.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {

            @Override
            public double visit(final int row, final int column, final double value) {
                return random.nextDouble() < noiseProportion ? Double.NEGATIVE_INFINITY : value;
            }
        });
        final RealMatrix betaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), input, PCATangentNormalizationUtils.EPSILON);
        final RealMatrix noisyBetaHats = PCATangentNormalizationUtils.calculateBetaHats(pon.getReducedPanelPInverseCounts(), noisyInput, PCATangentNormalizationUtils.EPSILON);
        final RealMatrix difference = betaHats.subtract(noisyBetaHats);
        difference.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {

            @Override
            public void visit(final int row, int column, double value) {
                Assert.assertEquals(value, 0, 0.01);
            }
        });
    }
}
Also used : HDF5PCACoveragePoN(org.broadinstitute.hellbender.tools.pon.coverage.pca.HDF5PCACoveragePoN) DefaultRealMatrixPreservingVisitor(org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) RealMatrix(org.apache.commons.math3.linear.RealMatrix) DefaultRealMatrixChangingVisitor(org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor) PCACoveragePoN(org.broadinstitute.hellbender.tools.pon.coverage.pca.PCACoveragePoN) HDF5PCACoveragePoN(org.broadinstitute.hellbender.tools.pon.coverage.pca.HDF5PCACoveragePoN) HDF5File(org.broadinstitute.hdf5.HDF5File)

Example 4 with DefaultRealMatrixPreservingVisitor

use of org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor in project gatk-protected by broadinstitute.

the class CopyNumberTriStateTransitionProbabilityCacheUnitTest method markovianPropertiesTest.

//test various properties of a transition matrix
@Test(dataProvider = "meanEventSizeAndEventStartProbability")
public void markovianPropertiesTest(final double meanEventSize, final double eventStartProbability) {
    final CopyNumberTriStateTransitionProbabilityCache cache = new CopyNumberTriStateTransitionProbabilityCache(meanEventSize, eventStartProbability);
    for (final int d : DISTANCES) {
        //check symmetries -- these are part of the model and need not be true in the future
        final RealMatrix transitionMatrix = cache.getAsMatrixInProbabilitySpace(d);
        assertSymmetries(transitionMatrix);
        //check that columns sums equal 1
        for (int column = 0; column < transitionMatrix.getColumnDimension(); column++) {
            Assert.assertEquals(MathUtils.sum(transitionMatrix.getColumn(column)), 1, EPSILON);
        }
        //check that all elements are positive
        transitionMatrix.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {

            @Override
            public void visit(int row, int column, double value) {
                Assert.assertTrue(value >= 0);
            }
        });
        //check that T(2d) = T(d)*T(d)
        assertEqualMatrices(cache.getAsMatrixInProbabilitySpace(2 * d), transitionMatrix.multiply(transitionMatrix));
        //check that the largest eigenvalue of the transition matrix is 1 (this corresponds to the asymptotic stationary state)
        Assert.assertEquals(MathUtils.arrayMax(new EigenDecomposition(transitionMatrix).getRealEigenvalues()), 1, EPSILON);
    }
    // check that at long distances memory of the initial state is lost and all initial distributions tend toward
    // the same asymptotic stationary distribution.  That is, all columns of the large-distance transition matrix are equal
    final RealMatrix asymptoticMatrix = cache.getAsMatrixInProbabilitySpace(HUGE_DISTANCE);
    for (int column = 1; column < asymptoticMatrix.getColumnDimension(); column++) {
        Assert.assertEquals(asymptoticMatrix.getColumnVector(0).subtract(asymptoticMatrix.getColumnVector(column)).getL1Norm(), 0, EPSILON);
    }
}
Also used : DefaultRealMatrixPreservingVisitor(org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor) EigenDecomposition(org.apache.commons.math3.linear.EigenDecomposition) RealMatrix(org.apache.commons.math3.linear.RealMatrix) Test(org.testng.annotations.Test)

Aggregations

DefaultRealMatrixPreservingVisitor (org.apache.commons.math3.linear.DefaultRealMatrixPreservingVisitor)4 RealMatrix (org.apache.commons.math3.linear.RealMatrix)4 Array2DRowRealMatrix (org.apache.commons.math3.linear.Array2DRowRealMatrix)2 DefaultRealMatrixChangingVisitor (org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor)2 EigenDecomposition (org.apache.commons.math3.linear.EigenDecomposition)2 HDF5File (org.broadinstitute.hdf5.HDF5File)2 HDF5PCACoveragePoN (org.broadinstitute.hellbender.tools.pon.coverage.pca.HDF5PCACoveragePoN)2 PCACoveragePoN (org.broadinstitute.hellbender.tools.pon.coverage.pca.PCACoveragePoN)2 Test (org.testng.annotations.Test)2