diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ExactAFCalculationModel.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ExactAFCalculationModel.java index 295d3f9f0..8da72ef7a 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ExactAFCalculationModel.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ExactAFCalculationModel.java @@ -70,59 +70,6 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { return genotypeLikelihoods; } - final static double approximateLog10SumLog10(double[] vals) { - - final int maxElementIndex = MathUtils.maxElementIndex(vals); - double approxSum = vals[maxElementIndex]; - if ( approxSum == Double.NEGATIVE_INFINITY ) - return approxSum; - - for ( int i = 0; i < vals.length; i++ ) { - if ( i == maxElementIndex || vals[i] == Double.NEGATIVE_INFINITY ) - continue; - - final double diff = approxSum - vals[i]; - if ( diff < MathUtils.MAX_JACOBIAN_TOLERANCE ) { - final int ind = fastRound(diff / MathUtils.JACOBIAN_LOG_TABLE_STEP); // hard rounding - approxSum += MathUtils.jacobianLogTable[ind]; - } - } - - return approxSum; - } - - final static double approximateLog10SumLog10(double small, double big) { - // make sure small is really the smaller value - if ( small > big ) { - final double t = big; - big = small; - small = t; - } - - if ( small == Double.NEGATIVE_INFINITY || big == Double.NEGATIVE_INFINITY ) - return big; - - final double diff = big - small; - if ( diff >= MathUtils.MAX_JACOBIAN_TOLERANCE ) - return big; - - // OK, so |y-x| < tol: we use the following identity then: - // we need to compute log10(10^x + 10^y) - // By Jacobian logarithm identity, this is equal to - // max(x,y) + log10(1+10^-abs(x-y)) - // we compute the second term as a table lookup with integer quantization - // we have pre-stored correction for 0,0.1,0.2,... 10.0 - final int ind = fastRound(diff / MathUtils.JACOBIAN_LOG_TABLE_STEP); // hard rounding - return big + MathUtils.jacobianLogTable[ind]; - } - - // A fast implementation of the Math.round() method. This method does not perform - // under/overflow checking, so this shouldn't be used in the general case (but is fine - // here because we already make those checks before calling in to the rounding). - final static int fastRound(double d) { - return (d > 0) ? (int)(d + 0.5d) : (int)(d - 0.5d); - } - // ------------------------------------------------------------------------------------- // // Multi-allelic implementation. @@ -403,7 +350,7 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { } } - final double log10Max = approximateLog10SumLog10(log10ConformationLikelihoods); + final double log10Max = MathUtils.approximateLog10SumLog10(log10ConformationLikelihoods); // finally, update the L(j,k) value set.log10Likelihoods[j] = log10Max - logDenominator; @@ -427,10 +374,10 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { // update the likelihoods/posteriors vectors which are collapsed views of each of the various ACs for ( int i = 0; i < set.ACcounts.getCounts().length; i++ ) { int AC = set.ACcounts.getCounts()[i]; - result.log10AlleleFrequencyLikelihoods[i][AC] = approximateLog10SumLog10(result.log10AlleleFrequencyLikelihoods[i][AC], log10LofK); + result.log10AlleleFrequencyLikelihoods[i][AC] = MathUtils.approximateLog10SumLog10(result.log10AlleleFrequencyLikelihoods[i][AC], log10LofK); final double prior = log10AlleleFrequencyPriors[nonRefAlleles-1][AC]; - result.log10AlleleFrequencyPosteriors[i][AC] = approximateLog10SumLog10(result.log10AlleleFrequencyPosteriors[i][AC], log10LofK + prior); + result.log10AlleleFrequencyPosteriors[i][AC] = MathUtils.approximateLog10SumLog10(result.log10AlleleFrequencyPosteriors[i][AC], log10LofK + prior); } } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UGBoundAF.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UGBoundAF.java index e40054c9f..99d55bc69 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UGBoundAF.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UGBoundAF.java @@ -204,6 +204,6 @@ public class UGBoundAF extends RodWalker { return Math.log10(s_2 + (s_2 - s)/15.0); } - return ExactAFCalculationModel.approximateLog10SumLog10(simpAux(likelihoods,a,c,eps/2,s_l,fa,fc,fd,cap-1),simpAux(likelihoods, c, b, eps / 2, s_r, fc, fb, fe, cap - 1)); + return MathUtils.approximateLog10SumLog10(simpAux(likelihoods,a,c,eps/2,s_l,fa,fc,fd,cap-1),simpAux(likelihoods, c, b, eps / 2, s_r, fc, fb, fe, cap - 1)); } } diff --git a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java index 737f4bb5f..5ffd634cc 100644 --- a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java @@ -56,17 +56,58 @@ public class MathUtils { private MathUtils() { } - @Requires({"d > 0.0"}) - public static int fastPositiveRound(double d) { - return (int) (d + 0.5); + // A fast implementation of the Math.round() method. This method does not perform + // under/overflow checking, so this shouldn't be used in the general case (but is fine + // if one is already make those checks before calling in to the rounding). + public static int fastRound(double d) { + return (d > 0) ? (int)(d + 0.5d) : (int)(d - 0.5d); } - public static int fastRound(double d) { - if (d > 0.0) { - return fastPositiveRound(d); - } else { - return -1 * fastPositiveRound(-1 * d); + public static double approximateLog10SumLog10(double[] vals) { + + final int maxElementIndex = MathUtils.maxElementIndex(vals); + double approxSum = vals[maxElementIndex]; + if ( approxSum == Double.NEGATIVE_INFINITY ) + return approxSum; + + for ( int i = 0; i < vals.length; i++ ) { + if ( i == maxElementIndex || vals[i] == Double.NEGATIVE_INFINITY ) + continue; + + final double diff = approxSum - vals[i]; + if ( diff < MathUtils.MAX_JACOBIAN_TOLERANCE ) { + // See notes from the 2-inout implementation below + final int ind = fastRound(diff / MathUtils.JACOBIAN_LOG_TABLE_STEP); // hard rounding + approxSum += MathUtils.jacobianLogTable[ind]; + } + } + + return approxSum; + } + + public static double approximateLog10SumLog10(double small, double big) { + // make sure small is really the smaller value + if ( small > big ) { + final double t = big; + big = small; + small = t; } + + if ( small == Double.NEGATIVE_INFINITY || big == Double.NEGATIVE_INFINITY ) + return big; + + final double diff = big - small; + if ( diff >= MathUtils.MAX_JACOBIAN_TOLERANCE ) + return big; + + // OK, so |y-x| < tol: we use the following identity then: + // we need to compute log10(10^x + 10^y) + // By Jacobian logarithm identity, this is equal to + // max(x,y) + log10(1+10^-abs(x-y)) + // we compute the second term as a table lookup with integer quantization + // we have pre-stored correction for 0,0.1,0.2,... 10.0 + final int ind = fastRound(diff / MathUtils.JACOBIAN_LOG_TABLE_STEP); // hard rounding + return big + MathUtils.jacobianLogTable[ind]; } public static double sum(Collection numbers) {