From 3b67aa8aeeeb2631b98b04cd7ebf26e35ae5a46a Mon Sep 17 00:00:00 2001
From: Mark DePristo <depristo@broadinstitute.org>
Date: Mon, 11 Feb 2013 12:42:50 -0800
Subject: [PATCH] Final edge case bug fixes to QualityUtil routines

-- log10 functions in QualityUtils allow -Infinity to allow log10(0.0) values
-- Fix edge condition of log10OneMinusX failing with Double.MIN_VALUE
-- Fix another edge condition of log10OneMinusX failing with a small but not min_value double
---
 .../broadinstitute/sting/utils/MathUtils.java | 22 +++++++++++++++----
 .../sting/utils/QualityUtils.java             |  6 +++--
 .../sting/utils/QualityUtilsUnitTest.java     | 12 ++++++++++
 3 files changed, 34 insertions(+), 6 deletions(-)

diff --git a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java
index 4db55b275..2459c1d36 100644
--- a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java
+++ b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java
@@ -1244,11 +1244,23 @@ public class MathUtils {
     /**
      * Checks that the result is a well-formed log10 probability
      *
-     * @param result a supposedly well-formed log10 probability value
+     * @param result a supposedly well-formed log10 probability value.  By default allows
+     *               -Infinity values, as log10(0.0) == -Infinity.
      * @return true if result is really well formed
      */
     public static boolean goodLog10Probability(final double result) {
-        return result <= 0.0 && ! Double.isInfinite(result) && ! Double.isNaN(result);
+        return goodLog10Probability(result, true);
+    }
+
+    /**
+     * Checks that the result is a well-formed log10 probability
+     *
+     * @param result a supposedly well-formed log10 probability value
+     * @param allowNegativeInfinity should we consider a -Infinity value ok?
+     * @return true if result is really well formed
+     */
+    public static boolean goodLog10Probability(final double result, final boolean allowNegativeInfinity) {
+        return result <= 0.0 && result != Double.POSITIVE_INFINITY && (allowNegativeInfinity || result != Double.NEGATIVE_INFINITY) && ! Double.isNaN(result);
     }
 
     /**
@@ -1779,7 +1791,9 @@ public class MathUtils {
             return Double.NEGATIVE_INFINITY;
         else if ( x == 0.0 )
             return 0.0;
-        else
-            return Math.log10(1 / x - 1) + Math.log10(x);
+        else {
+            final double d = Math.log10(1 / x - 1) + Math.log10(x);
+            return Double.isInfinite(d) || d > 0.0 ? 0.0 : d;
+        }
     }
 }
diff --git a/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java b/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java
index dd958cbb0..1dcd5a9ae 100644
--- a/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java
+++ b/public/java/src/org/broadinstitute/sting/utils/QualityUtils.java
@@ -312,7 +312,8 @@ public class QualityUtils {
      * This is a very generic method, that simply computes a phred-scaled double quality
      * score given an error rate.  It has the same precision as a normal double operation
      *
-     * @param trueRateLog10 the probability of being right (0.0-1.0)
+     * @param trueRateLog10 the log10 probability of being right (0.0-1.0).  Can be -Infinity to indicate
+     *                      that the result is impossible in which MIN_PHRED_SCALED_QUAL is returned
      * @return a phred-scaled version of the error rate implied by trueRate
      */
     @Ensures("result >= 0.0")
@@ -340,7 +341,8 @@ public class QualityUtils {
      * This is a very generic method, that simply computes a phred-scaled double quality
      * score given an error rate.  It has the same precision as a normal double operation
      *
-     * @param errorRateLog10 the log10 probability of being wrong (0.0-1.0)
+     * @param errorRateLog10 the log10 probability of being wrong (0.0-1.0).  Can be -Infinity, in which case
+     *                       the result is MIN_PHRED_SCALED_QUAL
      * @return a phred-scaled version of the error rate
      */
     @Ensures("result >= 0.0")
diff --git a/public/java/test/org/broadinstitute/sting/utils/QualityUtilsUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/QualityUtilsUnitTest.java
index 1efce3cb0..f5c7a14df 100644
--- a/public/java/test/org/broadinstitute/sting/utils/QualityUtilsUnitTest.java
+++ b/public/java/test/org/broadinstitute/sting/utils/QualityUtilsUnitTest.java
@@ -101,6 +101,18 @@ public class QualityUtilsUnitTest extends BaseTest {
         }
     }
 
+    @Test
+    public void testTrueProbWithMinDouble() {
+        final byte actual = QualityUtils.trueProbToQual(Double.MIN_VALUE);
+        Assert.assertEquals(actual, 1, "Failed to convert true prob of min double to 1 qual");
+    }
+
+    @Test
+    public void testTrueProbWithVerySmallValue() {
+        final byte actual = QualityUtils.trueProbToQual(1.7857786272673852E-19);
+        Assert.assertEquals(actual, 1, "Failed to convert true prob of very small value 1.7857786272673852E-19 to 1 qual");
+    }
+
     @Test
     public void testQualCaches() {
         Assert.assertEquals(QualityUtils.qualToErrorProb((byte) 20), 0.01, 1e-6);