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);
}
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);
}
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]));
}
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]));
}
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 });
}
Aggregations