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

Create the temp storage for calculating cell values just once as per Mark's TODO

This commit is contained in:
Eric Banks 2012-01-12 10:27:10 -05:00
parent f5f5ed5dcd
commit e7fe9910f7
2 changed files with 37 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -177,12 +177,18 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
ACqueue.add(zeroSet);
indexesToACset.put(zeroSet.ACcounts, zeroSet);
// optimization: create the temporary storage for computing L(j,k) just once
final int maxPossibleDependencies = numAlternateAlleles + (numAlternateAlleles * (numAlternateAlleles + 1) / 2) + 1;
final double[][] tempLog10ConformationLikelihoods = new double[numSamples+1][maxPossibleDependencies];
for ( int i = 0; i < maxPossibleDependencies; i++ )
tempLog10ConformationLikelihoods[0][i] = Double.NEGATIVE_INFINITY;
// keep processing while we have AC conformations that need to be calculated
double maxLog10L = Double.NEGATIVE_INFINITY;
while ( !ACqueue.isEmpty() ) {
// compute log10Likelihoods
final ExactACset set = ACqueue.remove();
final double log10LofKs = calculateAlleleCountConformation(set, genotypeLikelihoods, maxLog10L, numChr, preserveData, ACqueue, indexesToACset, log10AlleleFrequencyPriors, result);
final double log10LofKs = calculateAlleleCountConformation(set, genotypeLikelihoods, maxLog10L, numChr, preserveData, ACqueue, indexesToACset, log10AlleleFrequencyPriors, result, tempLog10ConformationLikelihoods);
// adjust max likelihood seen if needed
maxLog10L = Math.max(maxLog10L, log10LofKs);
@ -197,13 +203,14 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
final Queue<ExactACset> ACqueue,
final HashMap<ExactACcounts, ExactACset> indexesToACset,
final double[][] log10AlleleFrequencyPriors,
final AlleleFrequencyCalculationResult result) {
final AlleleFrequencyCalculationResult result,
final double[][] tempLog10ConformationLikelihoods) {
//if ( DEBUG )
// System.out.printf(" *** computing LofK for set=%s%n", set.ACcounts);
// compute the log10Likelihoods
computeLofK(set, genotypeLikelihoods, indexesToACset, log10AlleleFrequencyPriors, result);
computeLofK(set, genotypeLikelihoods, indexesToACset, log10AlleleFrequencyPriors, result, tempLog10ConformationLikelihoods);
// clean up memory
if ( !preserveData ) {
@ -309,7 +316,8 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
final ArrayList<double[]> genotypeLikelihoods,
final HashMap<ExactACcounts, ExactACset> indexesToACset,
final double[][] log10AlleleFrequencyPriors,
final AlleleFrequencyCalculationResult result) {
final AlleleFrequencyCalculationResult result,
final double[][] tempLog10ConformationLikelihoods) {
set.log10Likelihoods[0] = 0.0; // the zero case
final int totalK = set.getACsum();
@ -321,10 +329,6 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
}
// k > 0 for at least one k
else {
// all possible likelihoods for a given cell from which to choose the max
final int numPaths = set.ACsetIndexToPLIndex.size() + 1;
final double[][] log10ConformationLikelihoods = new double[set.log10Likelihoods.length][numPaths]; // TODO can be created just once, since you initialize it
// deal with the non-AA possible conformations
int conformationIndex = 1;
for ( Map.Entry<ExactACcounts, Integer> mapping : set.ACsetIndexToPLIndex.entrySet() ) {
@ -337,10 +341,10 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
if ( totalK <= 2*j ) { // skip impossible conformations
final double[] gl = genotypeLikelihoods.get(j);
log10ConformationLikelihoods[j][conformationIndex] =
tempLog10ConformationLikelihoods[j][conformationIndex] =
determineCoefficient(mapping.getValue(), j, set.ACcounts.getCounts(), totalK) + dependent.log10Likelihoods[j-1] + gl[mapping.getValue()];
} else {
log10ConformationLikelihoods[j][conformationIndex] = Double.NEGATIVE_INFINITY;
tempLog10ConformationLikelihoods[j][conformationIndex] = Double.NEGATIVE_INFINITY;
}
}
@ -348,17 +352,18 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
}
// finally, deal with the AA case (which depends on previous cells in this column) and then update the L(j,k) value
final int numPaths = set.ACsetIndexToPLIndex.size() + 1;
for ( int j = 1; j < set.log10Likelihoods.length; j++ ) {
if ( totalK < 2*j-1 ) {
final double[] gl = genotypeLikelihoods.get(j);
log10ConformationLikelihoods[j][0] = MathUtils.log10Cache[2*j-totalK] + MathUtils.log10Cache[2*j-totalK-1] + set.log10Likelihoods[j-1] + gl[HOM_REF_INDEX];
tempLog10ConformationLikelihoods[j][0] = MathUtils.log10Cache[2*j-totalK] + MathUtils.log10Cache[2*j-totalK-1] + set.log10Likelihoods[j-1] + gl[HOM_REF_INDEX];
} else {
log10ConformationLikelihoods[j][0] = Double.NEGATIVE_INFINITY;
tempLog10ConformationLikelihoods[j][0] = Double.NEGATIVE_INFINITY;
}
final double logDenominator = MathUtils.log10Cache[2*j] + MathUtils.log10Cache[2*j-1];
final double log10Max = MathUtils.approximateLog10SumLog10(log10ConformationLikelihoods[j]);
final double log10Max = MathUtils.approximateLog10SumLog10(tempLog10ConformationLikelihoods[j], numPaths);
set.log10Likelihoods[j] = log10Max - logDenominator;
}
}

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

@ -63,14 +63,18 @@ public class MathUtils {
return (d > 0) ? (int)(d + 0.5d) : (int)(d - 0.5d);
}
public static double approximateLog10SumLog10(double[] vals) {
public static double approximateLog10SumLog10(final double[] vals) {
return approximateLog10SumLog10(vals, vals.length);
}
final int maxElementIndex = MathUtils.maxElementIndex(vals);
public static double approximateLog10SumLog10(final double[] vals, final int endIndex) {
final int maxElementIndex = MathUtils.maxElementIndex(vals, endIndex);
double approxSum = vals[maxElementIndex];
if ( approxSum == Double.NEGATIVE_INFINITY )
return approxSum;
for ( int i = 0; i < vals.length; i++ ) {
for ( int i = 0; i < endIndex; i++ ) {
if ( i == maxElementIndex || vals[i] == Double.NEGATIVE_INFINITY )
continue;
@ -582,11 +586,15 @@ public class MathUtils {
return normalizeFromLog10(array, false);
}
public static int maxElementIndex(double[] array) {
public static int maxElementIndex(final double[] array) {
return maxElementIndex(array, array.length);
}
public static int maxElementIndex(final double[] array, final int endIndex) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
int maxI = -1;
for (int i = 0; i < array.length; i++) {
for (int i = 0; i < endIndex; i++) {
if (maxI == -1 || array[i] > array[maxI])
maxI = i;
}
@ -594,11 +602,15 @@ public class MathUtils {
return maxI;
}
public static int maxElementIndex(int[] array) {
public static int maxElementIndex(final int[] array) {
return maxElementIndex(array, array.length);
}
public static int maxElementIndex(final int[] array, int endIndex) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
int maxI = -1;
for (int i = 0; i < array.length; i++) {
for (int i = 0; i < endIndex; i++) {
if (maxI == -1 || array[i] > array[maxI])
maxI = i;
}