zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Moving the approximate summing of log10 vals to MathUtils; keeping the more efficient implementation of fast rounding.

This commit is contained in:
Eric Banks 2012-01-10 12:38:47 -05:00
parent 589397d611
commit 25d0d53d88
3 changed files with 53 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ExactAFCalculationModel.java
View File

@ -70,59 +70,6 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
return genotypeLikelihoods; 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. // 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 // finally, update the L(j,k) value
set.log10Likelihoods[j] = log10Max - logDenominator; 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 // 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++ ) { for ( int i = 0; i < set.ACcounts.getCounts().length; i++ ) {
int AC = set.ACcounts.getCounts()[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]; 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);
} }
} }
} }

2
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UGBoundAF.java
View File

@ -204,6 +204,6 @@ public class UGBoundAF extends RodWalker<VariantContext,Integer> {
return Math.log10(s_2 + (s_2 - s)/15.0); 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));
} }
} }

57
public/java/src/org/broadinstitute/sting/utils/MathUtils.java
View File

@ -56,17 +56,58 @@ public class MathUtils {
private MathUtils() { private MathUtils() {
} }
@Requires({"d > 0.0"}) // A fast implementation of the Math.round() method. This method does not perform
public static int fastPositiveRound(double d) { // under/overflow checking, so this shouldn't be used in the general case (but is fine
return (int) (d + 0.5); // 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) { public static double approximateLog10SumLog10(double[] vals) {
if (d > 0.0) {
return fastPositiveRound(d); final int maxElementIndex = MathUtils.maxElementIndex(vals);
} else { double approxSum = vals[maxElementIndex];
return -1 * fastPositiveRound(-1 * d); 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<Number> numbers) { public static double sum(Collection<Number> numbers) {