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 aa743f86f..295d3f9f0 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,15 +70,25 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { return genotypeLikelihoods; } - final static double approximateLog10SumLog10(double[] vals) { - if ( vals.length < 2 ) - throw new ReviewedStingException("Passing array with fewer than 2 values when computing approximateLog10SumLog10"); - double approx = approximateLog10SumLog10(vals[0], vals[1]); - for ( int i = 2; i < vals.length; i++ ) - approx = approximateLog10SumLog10(approx, vals[i]); - return approx; + 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) { @@ -89,27 +99,29 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel { small = t; } - if (small == Double.NEGATIVE_INFINITY || big == Double.NEGATIVE_INFINITY ) + if ( small == Double.NEGATIVE_INFINITY || big == Double.NEGATIVE_INFINITY ) return big; - if (big >= small + MathUtils.MAX_JACOBIAN_TOLERANCE) + 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 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 = (int)(((big-small)/JACOBIAN_LOG_TABLE_STEP)); // hard rounding - int ind = (int)(Math.round((big-small)/MathUtils.JACOBIAN_LOG_TABLE_STEP)); // hard rounding - - //double z =Math.log10(1+Math.pow(10.0,-diff)); - //System.out.format("x: %f, y:%f, app: %f, true: %f ind:%d\n",x,y,t2,z,ind); + 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); + } // ------------------------------------------------------------------------------------- // diff --git a/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java b/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java index 7c0dba558..646ede836 100755 --- a/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java @@ -28,7 +28,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { public void testMultiSamplePilot1() { WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec( baseCommand + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10,022,000-10,025,000", 1, - Arrays.asList("66ed60c6c1190754abd8a0a9d1d8d61e")); + Arrays.asList("d61c7055bd09024abb8902bde6bd3960")); executeTest("test MultiSample Pilot1", spec); } @@ -294,7 +294,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest { WalkerTest.WalkerTestSpec spec4 = new WalkerTest.WalkerTestSpec( baseCommandIndelsb37 + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + validationDataLocation + "ALL.wgs.union_v2_chr20_100_110K.20101123.indels.sites.vcf -I " + validationDataLocation + "phase1_GBR_realigned.chr20.100K-110K.bam -o %s -L 20:100,000-110,000", 1, - Arrays.asList("1e2a4aab26e9ab0dae709d33a669e036")); + Arrays.asList("69bfc9bec43a4fdd85dda3b947e6a98e")); executeTest("test MultiSample Phase1 indels with complicated records", spec4); }