From e7fe9910f7d8e3df40f2848231dfadf301f09424 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Thu, 12 Jan 2012 10:27:10 -0500 Subject: [PATCH] Create the temp storage for calculating cell values just once as per Mark's TODO --- .../genotyper/ExactAFCalculationModel.java | 31 +++++++++++-------- .../broadinstitute/sting/utils/MathUtils.java | 26 +++++++++++----- 2 files changed, 37 insertions(+), 20 deletions(-) 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 986b2a800..1594c92cb 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 @@ -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 ACqueue, final HashMap 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 genotypeLikelihoods, final HashMap 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 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; } } diff --git a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java index 5ffd634cc..408faf61a 100644 --- a/public/java/src/org/broadinstitute/sting/utils/MathUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/MathUtils.java @@ -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; }