Search in sources :

Example 41 with Array2DRowRealMatrix

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

the class SomaticLikelihoodsEngineUnitTest method testAlleleFractionsPosterior.

@Test
public void testAlleleFractionsPosterior() {
    //likelihoods completely favor allele 0 over allele 1 for every read, so
    // we should get no counts for allele 1
    final double[] prior1 = new double[] { 1, 1 };
    final RealMatrix mat1 = new Array2DRowRealMatrix(2, 4);
    mat1.setRow(0, new double[] { 0, 0, 0, 0 });
    mat1.setRow(1, new double[] { -10, -10, -10, -10 });
    final double[] posterior1 = SomaticLikelihoodsEngine.alleleFractionsPosterior(mat1, prior1);
    final double[] expectedCounts1 = new double[] { 4, 0 };
    final double[] expectedPosterior1 = MathArrays.ebeAdd(prior1, expectedCounts1);
    Assert.assertArrayEquals(posterior1, expectedPosterior1, 1.0e-6);
    //prior is extremely strong and outweighs ambiguous likelihoods
    final double[] prior2 = new double[] { 1e8, 1 };
    final RealMatrix mat2 = new Array2DRowRealMatrix(2, 4);
    mat2.setRow(0, new double[] { 0, 0, 0, 0 });
    mat2.setRow(1, new double[] { 0, 0, 0, 0 });
    final double[] posterior2 = SomaticLikelihoodsEngine.alleleFractionsPosterior(mat2, prior2);
    final double[] expectedCounts2 = new double[] { 4, 0 };
    final double[] expectedPosterior2 = MathArrays.ebeAdd(prior2, expectedCounts2);
    Assert.assertArrayEquals(posterior2, expectedPosterior2, 1.0e-6);
    //prior is extremely weak and likelihoods speak for themselves
    final double[] prior3 = new double[] { 1e-6, 1e-6 };
    final RealMatrix mat3 = new Array2DRowRealMatrix(2, 4);
    mat3.setRow(0, new double[] { 0, 0, 0, -10 });
    mat3.setRow(1, new double[] { -10, -10, -10, 0 });
    final double[] posterior3 = SomaticLikelihoodsEngine.alleleFractionsPosterior(mat3, prior3);
    final double[] expectedCounts3 = new double[] { 3, 1 };
    final double[] expectedPosterior3 = MathArrays.ebeAdd(prior3, expectedCounts3);
    Assert.assertArrayEquals(posterior3, expectedPosterior3, 1.0e-6);
    // test convergence
    final double[] prior4 = new double[] { 0.2, 1.7 };
    final RealMatrix mat4 = new Array2DRowRealMatrix(2, 4);
    mat4.setRow(0, new double[] { 0.1, 5.2, 0.5, 0.2 });
    mat4.setRow(1, new double[] { 2.6, 0.6, 0.5, 0.4 });
    final double[] posterior4 = SomaticLikelihoodsEngine.alleleFractionsPosterior(mat4, prior4);
    final double[] counts4 = MathArrays.ebeSubtract(posterior4, prior4);
    Assert.assertArrayEquals(counts4, SomaticLikelihoodsEngine.getEffectiveCounts(mat4, posterior4), 1.0e-3);
}
Also used : Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) RealMatrix(org.apache.commons.math3.linear.RealMatrix) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 42 with Array2DRowRealMatrix

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

the class SomaticLikelihoodsEngineUnitTest method testEvidence.

@Test
public void testEvidence() {
    // one exact limit for the evidence is when the likelihoods of each read are so peaked (i.e. the most likely allele
    // of each read is much likelier than all other alleles) that the sum over latent read-to-allele assignments
    // (that is, over the indicator z in the notes) is dominated by the max-likelihood allele configuration
    // and thus the evidence reduces to exactly integrating out the Dirichlet allele fractions
    final double[] prior = new double[] { 1, 2 };
    final RealMatrix log10Likelihoods = new Array2DRowRealMatrix(2, 4);
    log10Likelihoods.setRow(0, new double[] { 0.1, 4.0, 3.0, -10 });
    log10Likelihoods.setRow(1, new double[] { -12, -9, -5.0, 0.5 });
    final double calculatedLog10Evidence = SomaticLikelihoodsEngine.log10Evidence(log10Likelihoods, prior);
    final double[] maxLikelihoodCounts = new double[] { 3, 1 };
    final double expectedLog10Evidence = SomaticLikelihoodsEngine.log10DirichletNormalization(prior) - SomaticLikelihoodsEngine.log10DirichletNormalization(MathArrays.ebeAdd(prior, maxLikelihoodCounts)) + new IndexRange(0, log10Likelihoods.getColumnDimension()).sum(read -> log10Likelihoods.getColumnVector(read).getMaxValue());
    Assert.assertEquals(calculatedLog10Evidence, expectedLog10Evidence, 1e-5);
    // when there's just one read we can calculate the likelihood exactly
    final double[] prior2 = new double[] { 1, 2 };
    final RealMatrix log10Likelihoods2 = new Array2DRowRealMatrix(2, 1);
    log10Likelihoods2.setRow(0, new double[] { 0.1 });
    log10Likelihoods2.setRow(1, new double[] { 0.5 });
    final double calculatedLog10Evidence2 = SomaticLikelihoodsEngine.log10Evidence(log10Likelihoods2, prior2);
    final double[] delta0 = new double[] { 1, 0 };
    final double[] delta1 = new double[] { 0, 1 };
    final double expectedLog10Evidence2 = MathUtils.log10SumLog10(log10Likelihoods2.getEntry(0, 0) + SomaticLikelihoodsEngine.log10DirichletNormalization(prior2) - SomaticLikelihoodsEngine.log10DirichletNormalization(MathArrays.ebeAdd(prior2, delta0)), +log10Likelihoods2.getEntry(1, 0) + SomaticLikelihoodsEngine.log10DirichletNormalization(prior2) - SomaticLikelihoodsEngine.log10DirichletNormalization(MathArrays.ebeAdd(prior2, delta1)));
    Assert.assertEquals(calculatedLog10Evidence2, expectedLog10Evidence2, 0.05);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) MathUtils(org.broadinstitute.hellbender.utils.MathUtils) Assert(org.testng.Assert) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) MathArrays(org.apache.commons.math3.util.MathArrays) RealMatrix(org.apache.commons.math3.linear.RealMatrix) Test(org.testng.annotations.Test) Assert(org.junit.Assert) IndexRange(org.broadinstitute.hellbender.utils.IndexRange) IndexRange(org.broadinstitute.hellbender.utils.IndexRange) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) RealMatrix(org.apache.commons.math3.linear.RealMatrix) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 43 with Array2DRowRealMatrix

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

the class GATKProtectedMathUtilsTest method testRowNaN.

@Test
public void testRowNaN() {
    final double[][] array = { { 1, 2, Double.NaN }, { 5, 5, 5 }, { 7, 8, 9 }, { -15, 2, 12 } };
    final double[] guess = GATKProtectedMathUtils.rowSums(new Array2DRowRealMatrix(array));
    final double[] gt = { Double.NaN, 15, 24, -1 };
    Assert.assertEquals(guess.length, 4);
    Assert.assertEquals(guess[1], gt[1]);
    Assert.assertEquals(guess[2], gt[2]);
    Assert.assertEquals(guess[3], gt[3]);
    Assert.assertTrue(Double.isNaN(guess[0]));
}
Also used : Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) Test(org.testng.annotations.Test)

Example 44 with Array2DRowRealMatrix

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

the class GATKProtectedMathUtilsTest method testRowNegInf.

@Test
public void testRowNegInf() {
    final double[][] array = { { 1, 2, Double.NEGATIVE_INFINITY }, { 5, 5, 5 }, { 7, 8, 9 }, { -15, 2, 12 } };
    final double[] guess = GATKProtectedMathUtils.rowSums(new Array2DRowRealMatrix(array));
    final double[] gt = { Double.NEGATIVE_INFINITY, 15, 24, -1 };
    Assert.assertEquals(guess.length, 4);
    Assert.assertEquals(guess[1], gt[1]);
    Assert.assertEquals(guess[2], gt[2]);
    Assert.assertEquals(guess[3], gt[3]);
    Assert.assertTrue(Double.isInfinite(guess[0]));
}
Also used : Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) Test(org.testng.annotations.Test)

Example 45 with Array2DRowRealMatrix

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

the class GATKProtectedMathUtilsTest method testRowMeans.

@Test
public void testRowMeans() {
    final double[][] array = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
    final double[] rowMeans = GATKProtectedMathUtils.rowMeans(new Array2DRowRealMatrix(array));
    Assert.assertEquals(rowMeans, new double[] { 2, 5, 8 });
}
Also used : Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) Test(org.testng.annotations.Test)

Aggregations

Array2DRowRealMatrix (org.apache.commons.math3.linear.Array2DRowRealMatrix)141 RealMatrix (org.apache.commons.math3.linear.RealMatrix)101 Test (org.testng.annotations.Test)60 IntStream (java.util.stream.IntStream)31 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)28 File (java.io.File)27 Collectors (java.util.stream.Collectors)25 ArrayList (java.util.ArrayList)24 Assert (org.testng.Assert)24 List (java.util.List)22 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)22 Target (org.broadinstitute.hellbender.tools.exome.Target)18 java.util (java.util)15 Random (java.util.Random)14 ReadCountCollection (org.broadinstitute.hellbender.tools.exome.ReadCountCollection)14 ParamUtils (org.broadinstitute.hellbender.utils.param.ParamUtils)14 DataProvider (org.testng.annotations.DataProvider)14 Stream (java.util.stream.Stream)13 Arrays (java.util.Arrays)12 DoubleStream (java.util.stream.DoubleStream)12