Examples with RealMatrix org.apache.commons.math3.linear.RealMatrix used on opensource projects

Example 76 with RealMatrix

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

the class IntegerCopyNumberTransitionProbabilityCacheCollectionUnitTest method performCompleteTest.

private void performCompleteTest(final IntegerCopyNumberTransitionProbabilityCacheCollection cache, final boolean padded) {
    for (final String sexGenotype : HOMO_SAPIENS_SEX_GENOTYPES) {
        for (final int dist : DISTANCES) {
            for (final String contig : HOMO_SAPIENS_ALL_CONTIGS) {
                /* set the per-base transition matrix according to contig and sex genotype */
                final RealMatrix perBaseTransitionMatrix;
                if (HOMO_SAPIENS_AUTOSOMAL_CONTIGS.contains(contig)) {
                    perBaseTransitionMatrix = IntegerCopyNumberTransitionMatrixUnitTest.HOMO_SAPIENS_COPY_NUMBER_TRANSITION_AUTOSOMAL_TRUTH;
                } else if (contig.equals("X")) {
                    if (sexGenotype.equals("SEX_XX")) {
                        perBaseTransitionMatrix = IntegerCopyNumberTransitionMatrixUnitTest.HOMO_SAPIENS_COPY_NUMBER_TRANSITION_XX_X_TRUTH;
                    } else {
                        /* SEX_XY */
                        perBaseTransitionMatrix = IntegerCopyNumberTransitionMatrixUnitTest.HOMO_SAPIENS_COPY_NUMBER_TRANSITION_XY_X_TRUTH;
                    }
                } else {
                    /* contig = Y */
                    if (sexGenotype.equals("SEX_XX")) {
                        perBaseTransitionMatrix = IntegerCopyNumberTransitionMatrixUnitTest.HOMO_SAPIENS_COPY_NUMBER_TRANSITION_XX_Y_TRUTH;
                    } else {
                        /* SEX_XY */
                        perBaseTransitionMatrix = IntegerCopyNumberTransitionMatrixUnitTest.HOMO_SAPIENS_COPY_NUMBER_TRANSITION_XY_Y_TRUTH;
                    }
                }
                /* the list of copy number states */
                final List<IntegerCopyNumberState> integerCopyNumberStates = IntStream.range(0, cache.getMaxCopyNumber(sexGenotype, contig) + 1).mapToObj(IntegerCopyNumberState::new).collect(Collectors.toList());
                /* calculate the log transition matrix directly */
                final RealMatrix transitionMatrixDirect = getDirectMatrixPower(perBaseTransitionMatrix, dist);
                /* calculate the log transition matrix using the cache class */
                final RealMatrix transitionMatrixFromCache = getTransitionMatrix(dist, sexGenotype, contig, integerCopyNumberStates, integerCopyNumberStates, cache);
                /* if the collection is padded, we must pad the truth as well */
                final RealMatrix transitionMatrixDirectPadded;
                if (padded) {
                    transitionMatrixDirectPadded = MatrixUtils.createRealIdentityMatrix(MAX_COPY_NUMBER + 1);
                    for (int i = 0; i < transitionMatrixDirect.getRowDimension(); i++) {
                        for (int j = 0; j < transitionMatrixDirect.getColumnDimension(); j++) {
                            transitionMatrixDirectPadded.setEntry(i, j, transitionMatrixDirect.getEntry(i, j));
                        }
                    }
                } else {
                    transitionMatrixDirectPadded = transitionMatrixDirect;
                }
                assertEqualMatrices(transitionMatrixDirectPadded, transitionMatrixFromCache);
            }
        }
    }
}

Example 77 with RealMatrix

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

the class IntegerCopyNumberTransitionProbabilityCacheUnitTest method testBasicSoundness.

@Test
public void testBasicSoundness() {
    for (final RealMatrix transitionMatrix : TRANSITION_MATRICES) {
        final IntegerCopyNumberTransitionProbabilityCache cache = new IntegerCopyNumberTransitionProbabilityCache(new IntegerCopyNumberTransitionMatrix(transitionMatrix, 0));
        for (final int dist : DISTANCES) {
            final RealMatrix transitionMatrixExponentiated = cache.getTransitionProbabilityMatrix(dist);
            /* assert positivity */
            Assert.assertTrue(Arrays.stream(transitionMatrixExponentiated.getData()).flatMapToDouble(Arrays::stream).allMatch(d -> d >= 0));
            /* assert conservation of probability */
            for (int c = 0; c < transitionMatrix.getColumnDimension(); c++) {
                Assert.assertEquals(Arrays.stream(transitionMatrixExponentiated.getColumn(c)).sum(), 1.0, EPSILON);
            }
            /* assert correctness, T(2*d) = T(d)*T(d) */
            assertEqualMatrices(cache.getTransitionProbabilityMatrix(2 * dist), transitionMatrixExponentiated.multiply(transitionMatrixExponentiated));
        }
        /* assert loss of initial state over long distances, i.e. all columns must be equal */
        final RealMatrix longRangeTransitionMatrix = cache.getTransitionProbabilityMatrix(Integer.MAX_VALUE);
        final double[] firstColumn = longRangeTransitionMatrix.getColumn(0);
        final RealMatrix syntheticLongRangeTransitionMatrix = new Array2DRowRealMatrix(firstColumn.length, firstColumn.length);
        for (int i = 0; i < firstColumn.length; i++) {
            syntheticLongRangeTransitionMatrix.setColumn(i, firstColumn);
        }
        assertEqualMatrices(longRangeTransitionMatrix, syntheticLongRangeTransitionMatrix);
        final double[] stationary = cache.getStationaryProbabilityVector().toArray();
        ArrayAsserts.assertArrayEquals(stationary, firstColumn, EPSILON);
    }
}

Example 78 with RealMatrix

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

the class IntegerCopyNumberTransitionMatrixUnitTest method testPadding.

@Test
public void testPadding() {
    final IntegerCopyNumberTransitionMatrix data = new IntegerCopyNumberTransitionMatrix(new Array2DRowRealMatrix(new double[][] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } }), 2);
    final RealMatrix expected = new Array2DRowRealMatrix(new double[][] { { 1.0 / 12, 2.0 / 15, 3.0 / 18, 0, 0 }, { 4.0 / 12, 5.0 / 15, 6.0 / 18, 0, 0 }, { 7.0 / 12, 8.0 / 15, 9.0 / 18, 0, 0 }, { 0, 0, 0, 1, 0 }, { 0, 0, 0, 0, 1 } });
    Assert.assertEquals(data.getTransitionMatrix().subtract(expected).getNorm(), 0, 1e-12);
}

Example 79 with RealMatrix

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

the class SomaticGenotypingEngine method getAsRealMatrix.

//convert a likelihood matrix of alleles x reads into a RealMatrix
public static RealMatrix getAsRealMatrix(final LikelihoodMatrix<Allele> matrix) {
    final RealMatrix result = new Array2DRowRealMatrix(matrix.numberOfAlleles(), matrix.numberOfReads());
    result.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {

        @Override
        public double visit(int row, int column, double value) {
            return matrix.get(row, column);
        }
    });
    return result;
}

Example 80 with RealMatrix

use of org.apache.commons.math3.linear.RealMatrix 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);
    }
}