Minor refactoring: add ability to MathUtils.normalizeFromLog10 to not go to linear domain but just substract max value from log values and return. Use this function in snp and indel GL computation.

This commit is contained in:
Guillermo del Angel 2011-09-25 21:18:56 -04:00
parent 8c4eb5827f
commit 9afccd11b1
3 changed files with 19 additions and 14 deletions

View File

@ -122,13 +122,10 @@ public class SNPGenotypeLikelihoodsCalculationModel extends GenotypeLikelihoodsC
aList.add(refAllele); aList.add(refAllele);
aList.add(altAllele); aList.add(altAllele);
double[] dlike = new double[]{likelihoods[refGenotype.ordinal()],likelihoods[hetGenotype.ordinal()],likelihoods[homGenotype.ordinal()]} ; double[] dlike = new double[]{likelihoods[refGenotype.ordinal()],likelihoods[hetGenotype.ordinal()],likelihoods[homGenotype.ordinal()]} ;
double maxElement = MathUtils.max(dlike[AlleleFrequencyCalculationModel.GenotypeType.AA.ordinal()],
dlike[AlleleFrequencyCalculationModel.GenotypeType.AB.ordinal()],dlike[AlleleFrequencyCalculationModel.GenotypeType.BB.ordinal()]);
for (int i=0; i < dlike.length; i++)
dlike[i] -= maxElement;
// normalize in log space so that max element is zero.
GLs.put(sample.getKey(), new MultiallelicGenotypeLikelihoods(sample.getKey(), GLs.put(sample.getKey(), new MultiallelicGenotypeLikelihoods(sample.getKey(),
aList, dlike, getFilteredDepth(pileup))); aList, MathUtils.normalizeFromLog10(dlike, false, true), getFilteredDepth(pileup)));
} }
return refAllele; return refAllele;

View File

@ -1041,20 +1041,14 @@ public class PairHMMIndelErrorModel {
double[] genotypeLikelihoods = new double[hSize*(hSize+1)/2]; double[] genotypeLikelihoods = new double[hSize*(hSize+1)/2];
int k=0; int k=0;
double maxElement = Double.NEGATIVE_INFINITY;
for (int j=0; j < hSize; j++) { for (int j=0; j < hSize; j++) {
for (int i=0; i <= j; i++){ for (int i=0; i <= j; i++){
genotypeLikelihoods[k++] = haplotypeLikehoodMatrix[i][j]; genotypeLikelihoods[k++] = haplotypeLikehoodMatrix[i][j];
if (haplotypeLikehoodMatrix[i][j] > maxElement)
maxElement = haplotypeLikehoodMatrix[i][j];
} }
} }
// renormalize // renormalize so that max element is zero.
for (int i=0; i < genotypeLikelihoods.length; i++) return MathUtils.normalizeFromLog10(genotypeLikelihoods, false, true);
genotypeLikelihoods[i] -= maxElement;
return genotypeLikelihoods;
} }
/** /**

View File

@ -444,11 +444,25 @@ public class MathUtils {
* @return a newly allocated array corresponding the normalized values in array, maybe log10 transformed * @return a newly allocated array corresponding the normalized values in array, maybe log10 transformed
*/ */
public static double[] normalizeFromLog10(double[] array, boolean takeLog10OfOutput) { public static double[] normalizeFromLog10(double[] array, boolean takeLog10OfOutput) {
double[] normalized = new double[array.length]; return normalizeFromLog10(array, takeLog10OfOutput, false);
}
public static double[] normalizeFromLog10(double[] array, boolean takeLog10OfOutput, boolean keepInLogSpace) {
// for precision purposes, we need to add (or really subtract, since they're // for precision purposes, we need to add (or really subtract, since they're
// all negative) the largest value; also, we need to convert to normal-space. // all negative) the largest value; also, we need to convert to normal-space.
double maxValue = Utils.findMaxEntry(array); double maxValue = Utils.findMaxEntry(array);
// we may decide to just normalize in log space with converting to linear space
if (keepInLogSpace) {
for (int i = 0; i < array.length; i++)
array[i] -= maxValue;
return array;
}
// default case: go to linear space
double[] normalized = new double[array.length];
for (int i = 0; i < array.length; i++) for (int i = 0; i < array.length; i++)
normalized[i] = Math.pow(10, array[i] - maxValue); normalized[i] = Math.pow(10, array[i] - maxValue);