handle log

broadinstitute · Sep 5, 2015 · d3d2707 · d3d2707
1 parent 352e5fd
commit d3d2707
Show file tree

Hide file tree

Showing 6 changed files with 57 additions and 16 deletions.
diff --git a/src/main/java/org/broadinstitute/hellbender/utils/MathUtils.java b/src/main/java/org/broadinstitute/hellbender/utils/MathUtils.java
@@ -24,6 +24,13 @@ public final class MathUtils {
 
     public static final double LOG10_ONE_HALF = Math.log10(0.5);
 
+    /**
+     * Log10 of the e constant.
+     */
+    public static final double LOG10_OF_E = Math.log10(Math.E);
+
+    private static final double LN_10 = Math.log(10);
+
     /**
      * Private constructor.  No instantiating this class!
      */
@@ -125,6 +132,24 @@ public static int fastRound(final double d) {
         return (d > 0.0) ? (int) (d + 0.5d) : (int) (d - 0.5d);
     }
 
+    /**
+     * Converts LOG10 to LN
+     * @param log10 log10(x)
+     * @return ln(x)
+     */
+    public static double log10ToLog(final double log10){
+        return log10 * LN_10;
+    }
+
+    /**
+     * Converts LN to LOG10
+     * @param ln log(x)
+     * @return log10(x)
+     */
+    public static double logToLog10(final double ln) {
+        return ln * LOG10_OF_E;
+    }
+
     public static double approximateLogSumLog(final double[] vals) {
         return approximateLogSumLog(vals, vals.length);
     }

diff --git a/src/main/java/org/broadinstitute/hellbender/utils/QualityUtils.java b/src/main/java/org/broadinstitute/hellbender/utils/QualityUtils.java
@@ -185,7 +185,7 @@ public static double qualToLogErrorProb(final byte qual) {
      * The calculation is extremely efficient
      *
      * @param qual a phred-scaled quality score encoded as a double
-     * @return a probability (0.0-1.0)
+     * @return log of probability (0.0-1.0)
      */
     public static double qualToLogErrorProb(final double qual) {
         if ( qual < 0.0 ) throw new IllegalArgumentException("qual must be >= 0.0 but got " + qual);

diff --git a/.../hellbender/tools/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngineUnitTest.java b/.../hellbender/tools/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngineUnitTest.java
@@ -44,7 +44,7 @@ public void setup() {
 
     @Test
     public void testNormalizeDiploidLikelihoodMatrixFromLog10() {
-        double[][] likelihoodMatrix = {
+        double[][] log10_likelihoodMatrix = {
             {-90.2,     0,      0},
             {-190.1, -2.1,      0},
             {-7.0,  -17.5,  -35.9}
@@ -56,9 +56,9 @@ public void testNormalizeDiploidLikelihoodMatrixFromLog10() {
         };
 
 
-        Assert.assertTrue(compareDoubleArrays(normalizeDiploidLikelihoodMatrixFromLog10(likelihoodMatrix), normalizedMatrix));
+        Assert.assertTrue(compareDoubleArrays(normalizeDiploidLikelihoodMatrixFromLog10(log10_likelihoodMatrix), normalizedMatrix));
 
-        double[][] likelihoodMatrix2 = {
+        double[][] log10_likelihoodMatrix2 = {
                 {-90.2,     0,      0,        0},
                 {-190.1, -2.1,      0,        0},
                 {-7.0,  -17.5,  -35.9,        0},
@@ -70,14 +70,15 @@ public void testNormalizeDiploidLikelihoodMatrixFromLog10() {
                 {-4.9,  -15.4,  -33.8,        0},
                 {-4.9,  -15.4,  -33.8,   -997.9},
         };
-        Assert.assertTrue(compareDoubleArrays(normalizeDiploidLikelihoodMatrixFromLog10(likelihoodMatrix2), normalizedMatrix2));
+        Assert.assertTrue(compareDoubleArrays(normalizeDiploidLikelihoodMatrixFromLog10(log10_likelihoodMatrix2), normalizedMatrix2));
     }
 
-    static double[][] normalizeDiploidLikelihoodMatrixFromLog10( final double[][] likelihoodMatrix ) {
-        final double[] genotypeLikelihoods = copyLowerTriangleToArray(likelihoodMatrix);
-        final double[] normalizedGenotypeLikelihoods = MathUtils.normalizeFromLog10(genotypeLikelihoods, false, true);
-        copyArrayToLowerTriangle(normalizedGenotypeLikelihoods, likelihoodMatrix);
-        return likelihoodMatrix;
+    static double[][] normalizeDiploidLikelihoodMatrixFromLog10( final double[][] log10_likelihoodMatrix ) {
+        final double[] genotypeLog10Likelihoods = copyLowerTriangleToArray(log10_likelihoodMatrix);
+        //Note: normalizing from log10 is the same as normalizing from log
+        final double[] normalizedGenotypeLikelihoods = MathUtils.normalizeFromLog(genotypeLog10Likelihoods, false, true);
+        copyArrayToLowerTriangle(normalizedGenotypeLikelihoods, log10_likelihoodMatrix);
+        return log10_likelihoodMatrix;
     }
 
     private static double[] copyLowerTriangleToArray(final double[][] from) {
@@ -166,7 +167,7 @@ public void testComputeLikelihoods(){
         PairHMMLikelihoodCalculationEngine.writeLikelihoodsToFile = true;
 
         final ReadLikelihoodCalculationEngine lce = new PairHMMLikelihoodCalculationEngine((byte) SAMUtils.MAX_PHRED_SCORE,
-                PairHMM.Implementation.LOGLESS_CACHING, QualityUtils.qualToErrorProbLog10(LEAC.phredScaledGlobalReadMismappingRate),
+                PairHMM.Implementation.LOGLESS_CACHING, MathUtils.logToLog10(QualityUtils.qualToLogErrorProb(LEAC.phredScaledGlobalReadMismappingRate)),
                 PairHMMLikelihoodCalculationEngine.PCRErrorModel.CONSERVATIVE);
 
         final Map<String, List<GATKRead>> perSampleReadList= new HashMap<>();

diff --git a/src/test/java/org/broadinstitute/hellbender/utils/MathUtilsUnitTests.java b/src/test/java/org/broadinstitute/hellbender/utils/MathUtilsUnitTests.java
@@ -13,6 +13,8 @@
 import java.util.Random;
 
 import static java.lang.Math.log10;
+import static org.broadinstitute.hellbender.utils.MathUtils.log10ToLog;
+import static org.broadinstitute.hellbender.utils.MathUtils.logToLog10;
 
 /**
  * Basic unit test for MathUtils
@@ -300,4 +302,17 @@ public void testNormalizeFromReal(){
             Assert.assertEquals(expected[i], actual[i], error);
         }
     }
+
+    @Test
+    public void testLogLog10conversions(){
+        final double error = 1e-6;
+        for (final double x : new double[]{Math.E, 10.0, 0.5, 1, 1.5, 100.0}) {
+            final double logX = Math.log(x);
+            final double log10X = Math.log10(x);
+            Assert.assertEquals(log10ToLog(log10X), logX, error, "log10ToLog");
+            Assert.assertEquals(logToLog10(logX), log10X, error, "logToLog10");
+            Assert.assertEquals(logToLog10(log10ToLog(log10X)), log10X, error, "log10->log->log10");
+            Assert.assertEquals(log10ToLog(logToLog10(logX)), logX, error, "log->log10->log");
+        }
+    }
 }
diff --git a/...ava/org/broadinstitute/hellbender/utils/genotyper/PerReadAlleleLikelihoodMapUnitTest.java b/...ava/org/broadinstitute/hellbender/utils/genotyper/PerReadAlleleLikelihoodMapUnitTest.java
@@ -349,9 +349,9 @@ private static Pair<Allele, Allele> getMostLikelyDiploidAlleles(final PerReadAll
                 double haplotypeLikelihood = 0.0;
                 for( final Map<Allele,Double> alleleMap : map.getAlleleToLikelihoodMaps() ) {
                     // Compute log10(10^x1/2 + 10^x2/2) = log10(10^x1+10^x2)-log10(2)
-                    final double likelihood_i = alleleMap.get(allele_i);
-                    final double likelihood_j = alleleMap.get(allele_j);
-                    haplotypeLikelihood += MathUtils.approximateLog10SumLog10(likelihood_i, likelihood_j) + MathUtils.LOG10_ONE_HALF;
+                    final double log_likelihood_i = MathUtils.log10ToLog(alleleMap.get(allele_i));
+                    final double log_likelihood_j = MathUtils.log10ToLog(alleleMap.get(allele_j));
+                    haplotypeLikelihood += MathUtils.logToLog10(MathUtils.approximateLogSumLog(log_likelihood_i, log_likelihood_j) + MathUtils.LOG10_ONE_HALF);
 
                     // fast exit.  If this diploid pair is already worse than the max, just stop and look at the next pair
                     if ( haplotypeLikelihood < maxElement ) {

diff --git a/src/test/java/org/broadinstitute/hellbender/utils/pairhmm/PairHMMUnitTest.java b/src/test/java/org/broadinstitute/hellbender/utils/pairhmm/PairHMMUnitTest.java
@@ -376,7 +376,7 @@ void testReadSameAsHaplotype(final PairHMM hmm, final int readSize) {
         final byte gcp = 10;
         hmm.initialize(readBases.length, refBases.length);
         // running HMM with no haplotype caching. Should therefore pass null in place of nextRef bases
-        final double d = hmm.computeReadLikelihoodGivenHaplotypeLog10(refBases, readBases,
+        final double d = hmm.computeReadLogLikelihoodGivenHaplotype(refBases, readBases,
                 Utils.dupBytes(baseQual, readBases.length),
                 Utils.dupBytes(insQual, readBases.length),
                 Utils.dupBytes(delQual, readBases.length),
@@ -403,7 +403,7 @@ void testAllMatchingRead(final PairHMM hmm, final int readSize, final int refSiz
         Assert.assertEquals(d, expected, 1e-3, "Likelihoods should sum to just the error prob of the read " + String.format("readSize=%d refSize=%d", readSize, refSize));
     }
 
-    private static double getExpectedMathingLogLikelihood(final byte[] readBases, final byte[] refBases, final byte baseQual, final byte insQual) {
+    private static double getExpectedMatchingLogLikelihood(final byte[] readBases, final byte[] refBases, final byte baseQual, final byte insQual) {
         double expected =  0;
         final double initialCondition = ((double) Math.abs(refBases.length - readBases.length + 1))/refBases.length;
         if (readBases.length < refBases.length) {