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

Updated to work exclusively in log10 space 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4069 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
fromer 2010-08-19 21:31:07 +00:00
parent 3af4e618cc
commit 1c4784999a
2 changed files with 25 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
java/src/org/broadinstitute/sting/playground/gatk/walkers/ReadBackedPhasingWalker.java
View File

@ -282,15 +282,16 @@ public class ReadBackedPhasingWalker extends LocusWalker<PhasingStatsAndOutput,
logger.debug("\nPhasing table [AFTER NORMALIZATION]:\n" + sampleHaps + "\n"); logger.debug("\nPhasing table [AFTER NORMALIZATION]:\n" + sampleHaps + "\n");
PhasingTable.PhasingTableEntry maxEntry = sampleHaps.maxEntry(); PhasingTable.PhasingTableEntry maxEntry = sampleHaps.maxEntry();
double posteriorProb = maxEntry.getScore().getValue();
// convert posteriorProb to PHRED scale, but do NOT cap the quality as in QualityUtils.probToQual(posteriorProb):
PreciseNonNegativeDouble sumErrorProbs = new PreciseNonNegativeDouble(ZERO); PreciseNonNegativeDouble sumErrorProbs = new PreciseNonNegativeDouble(ZERO);
for (PhasingTable.PhasingTableEntry pte : sampleHaps) { for (PhasingTable.PhasingTableEntry pte : sampleHaps) {
if (pte != maxEntry) if (pte != maxEntry)
sumErrorProbs.plusEqual(pte.getScore()); sumErrorProbs.plusEqual(pte.getScore());
} }
// log10(x) = log(x) / log(10): double phaseQuality = -10.0 * (sumErrorProbs.getLog10Value());
double phaseQuality = -10.0 * (sumErrorProbs.getLogValue() / Math.log(10.0)); // convert to PHRED scale, do NOT cap the quality as in QualityUtils.probToQual(posteriorProb)
double posteriorProb = maxEntry.getScore().getValue();
logger.debug("MAX hap:\t" + maxEntry.getHaplotypeClass() + "\tposteriorProb:\t" + posteriorProb + "\tphaseQuality:\t" + phaseQuality); logger.debug("MAX hap:\t" + maxEntry.getHaplotypeClass() + "\tposteriorProb:\t" + posteriorProb + "\tphaseQuality:\t" + phaseQuality);
if (statsWriter != null) if (statsWriter != null)

40
java/src/org/broadinstitute/sting/utils/PreciseNonNegativeDouble.java
View File

@ -9,43 +9,43 @@ package org.broadinstitute.sting.utils;
*/ */
/* PreciseNonNegativeDouble permits arithmetic operations on NON-NEGATIVE double values /* PreciseNonNegativeDouble permits arithmetic operations on NON-NEGATIVE double values
with precision (prevents underflow by representing in log space). with precision (prevents underflow by representing in log10 space).
*/ */
public class PreciseNonNegativeDouble implements Comparable<PreciseNonNegativeDouble> { public class PreciseNonNegativeDouble implements Comparable<PreciseNonNegativeDouble> {
private static double EQUALS_THRESH = 1e-6; private static double EQUALS_THRESH = 1e-6;
private static double INFINITY = Double.POSITIVE_INFINITY; private static double INFINITY = Double.POSITIVE_INFINITY;
private double logValue; private double log10Value;
public PreciseNonNegativeDouble(double d) { public PreciseNonNegativeDouble(double d) {
this(d, false); this(d, false);
} }
public PreciseNonNegativeDouble(double d, boolean isLog) { public PreciseNonNegativeDouble(double d, boolean isLog10) {
if (isLog) { if (isLog10) {
this.logValue = d; this.log10Value = d;
} }
else { else {
if (d < 0) if (d < 0)
throw new IllegalArgumentException("non-log PreciseNonNegativeDouble argument must be non-negative"); throw new IllegalArgumentException("non-log PreciseNonNegativeDouble argument must be non-negative");
this.logValue = Math.log(d); this.log10Value = Math.log10(d);
} }
} }
public PreciseNonNegativeDouble(PreciseNonNegativeDouble pd) { public PreciseNonNegativeDouble(PreciseNonNegativeDouble pd) {
this.logValue = pd.logValue; this.log10Value = pd.log10Value;
} }
public double getValue() { public double getValue() {
return Math.exp(logValue); return Math.pow(10, log10Value);
} }
public double getLogValue() { public double getLog10Value() {
return logValue; return log10Value;
} }
public PreciseNonNegativeDouble setEqual(PreciseNonNegativeDouble other) { public PreciseNonNegativeDouble setEqual(PreciseNonNegativeDouble other) {
logValue = other.logValue; log10Value = other.log10Value;
return this; return this;
} }
@ -62,12 +62,12 @@ public class PreciseNonNegativeDouble implements Comparable<PreciseNonNegativeDo
} }
public PreciseNonNegativeDouble absDiff(PreciseNonNegativeDouble other) { public PreciseNonNegativeDouble absDiff(PreciseNonNegativeDouble other) {
return new PreciseNonNegativeDouble(absSubLog(this.logValue, other.logValue), true); return new PreciseNonNegativeDouble(absSubLog(this.log10Value, other.log10Value), true);
} }
public int compareTo(PreciseNonNegativeDouble other) { public int compareTo(PreciseNonNegativeDouble other) {
// Since log is monotonic: e^a R e^b <=> a R b, where R is one of: >, <, == // Since log is monotonic: e^a R e^b <=> a R b, where R is one of: >, <, ==
double logValDiff = this.logValue - other.logValue; double logValDiff = this.log10Value - other.log10Value;
if (Math.abs(logValDiff) <= EQUALS_THRESH) if (Math.abs(logValDiff) <= EQUALS_THRESH)
return 0; // this.equals(other) return 0; // this.equals(other)
@ -87,23 +87,23 @@ public class PreciseNonNegativeDouble implements Comparable<PreciseNonNegativeDo
} }
public PreciseNonNegativeDouble plusEqual(PreciseNonNegativeDouble other) { public PreciseNonNegativeDouble plusEqual(PreciseNonNegativeDouble other) {
logValue = addInLogSpace(logValue, other.logValue); log10Value = addInLogSpace(log10Value, other.log10Value);
return this; return this;
} }
public PreciseNonNegativeDouble timesEqual(PreciseNonNegativeDouble other) { public PreciseNonNegativeDouble timesEqual(PreciseNonNegativeDouble other) {
logValue += other.logValue; log10Value += other.log10Value;
return this; return this;
} }
public PreciseNonNegativeDouble divEqual(PreciseNonNegativeDouble other) { public PreciseNonNegativeDouble divEqual(PreciseNonNegativeDouble other) {
logValue -= other.logValue; log10Value -= other.log10Value;
return this; return this;
} }
// If x = log(a), y = log(b), returns log(a+b) // If x = log(a), y = log(b), returns log(a+b)
public static double addInLogSpace(double x, double y) { public static double addInLogSpace(double x, double y) {
if (x == INFINITY || y == INFINITY) return INFINITY; //log( e^INFINITY + e^y ) = INFINITY if (x == INFINITY || y == INFINITY) return INFINITY; // log(e^INFINITY + e^y) = INFINITY
if (x == -INFINITY) return y; if (x == -INFINITY) return y;
if (y == -INFINITY) return x; if (y == -INFINITY) return x;
@ -119,7 +119,7 @@ public class PreciseNonNegativeDouble implements Comparable<PreciseNonNegativeDo
} }
// x + log(1+e^(y-x)) = log(a) + log(1+e^(log(b)-log(a))) = log(a) + log(1+b/a) = log(a+b) // x + log(1+e^(y-x)) = log(a) + log(1+e^(log(b)-log(a))) = log(a) + log(1+b/a) = log(a+b)
return maxVal + Math.log(1.0 + Math.exp(negDiff)); return maxVal + Math.log10(1.0 + Math.pow(10, negDiff));
} }
// If x = log(a), y = log(b), returns log |a-b| // If x = log(a), y = log(b), returns log |a-b|
@ -129,9 +129,9 @@ public class PreciseNonNegativeDouble implements Comparable<PreciseNonNegativeDo
return -INFINITY; return -INFINITY;
} }
else if (x >= y) // x + log(1-e^(y-x)) = log(a) + log(1-e^(log(b)-log(a))) = log(a) + log(1-b/a) = a - b = |a-b|, since x >= y else if (x >= y) // x + log(1-e^(y-x)) = log(a) + log(1-e^(log(b)-log(a))) = log(a) + log(1-b/a) = a - b = |a-b|, since x >= y
return x + Math.log(1 - Math.exp(y-x)); return x + Math.log10(1 - Math.pow(10, y-x));
else else
return y + Math.log(1 - Math.exp(x-y)); return y + Math.log10(1 - Math.pow(10, x-y));
} }
public String toString() { public String toString() {