From bcba1ff424a440b6d641eb2f8c5f3eba1da6e374 Mon Sep 17 00:00:00 2001
From: jmaguire <jmaguire@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Sun, 19 Apr 2009 17:52:24 +0000
Subject: [PATCH] Fix a minor rounding bug and putz around with fractional
 counts in the pooled caller.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@472 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../gatk/walkers/AlleleFrequencyWalker.java   | 42 +++++++------------
 .../gatk/walkers/PoolCallingExperiment.java   | 11 +++--
 .../utils/AlleleFrequencyEstimate.java        |  9 +++-
 3 files changed, 31 insertions(+), 31 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
index a78e042cc..998e71a17 100755
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
@@ -16,6 +16,7 @@ import java.util.List;
 import java.util.Arrays;
 import java.util.Random;
 import java.io.PrintStream;
+import java.io.DataInputStream;
 
 public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate, String>// implements AllelicVariant
 {
@@ -144,25 +145,6 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
         return alleleFreq;
     }
 
-    /*
-    static public String getBases (LocusContext context)
-    {
-        // Convert bases to CharArray
-        int numReads = context.getReads().size(); //numReads();
-        char[] bases = new char[numReads];
-        //int refnum = nuc2num[ref];
-
-        List<SAMRecord> reads = context.getReads();
-        List<Integer> offsets = context.getOffsets();
-        for (int i =0; i < numReads; i++ ) {
-            SAMRecord read = reads.get(i);
-            int offset = offsets.get(i);
-            bases[i] = read.getReadString().charAt(offset);
-        }
-        return new String(bases);
-    }
-    */
-
     /*
     static public String[] getIndels(LocusContext context)
     {
@@ -325,6 +307,8 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
             bestMixtures[0] = new MixtureLikelihood(posterior_null_hyp, 0.0, 0.0);
             posteriors[0] = posterior_null_hyp;
 
+            //System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", location, 0.0, 0.0, P_D_q(bases, quals, 0.0, refnum, altnum), P_q_G(bases, N, 0.0, 0, 0), P_G(N, 0), prior_alt_frequency, posterior_null_hyp, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
+
             assert(! Double.isNaN(P_D_q(bases, quals, 0.0, refnum, altnum)));
             assert(! Double.isNaN(P_q_G(bases, N, 0.0, 0, 0)));
             assert(! Double.isNaN(P_G(N, 0)));
@@ -332,8 +316,6 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
             assert(! Double.isNaN(posteriors[0]));
             assert(! Double.isInfinite(posteriors[0]));
 
-            //System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", location, 0.0, 0.0, P_D_q(bases, quals, 0.0, refnum, altnum), P_q_G(bases, N, 0.0, 0, 0), P_G(N, 0), posterior_null_hyp, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
-
             // qstar - true allele balances
             //for (qstar = epsilon + ((1.0 - 2*epsilon)/N), qstar_N = 1; qstar <= 1.0; qstar += (1.0 - 2*epsilon)/N, qstar_N++) 
             for (qstar_N = 1; qstar_N <= N; qstar_N += 1)
@@ -347,8 +329,8 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
                 bestMixtures[qstar_N] = new MixtureLikelihood(posterior, qstar, q);
                 posteriors[qstar_N] = posterior;
 
-                //System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", 
-                //        location, q, qstar, pDq, pqG, pG, posterior, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
+                //System.out.format("DBG %s %.2f %.2f %5.2f %5.2f %5.2f %5.2f %5.2f %d %d %s\n", 
+                        location, q, qstar, pDq, pqG, pG, prior_alt_frequency, posterior, (int)(q*bases.length), (int)((1.0-q)*bases.length), new String(bases));
             }
         }
 
@@ -512,7 +494,10 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
         {
             // For small n and the edges, compute it directly.
             ans = Math.log10((double)nchoosek(n, k)) + Math.log10(Math.pow(p, k)) + Math.log10(Math.pow(1-p, n-k));
-            //System.out.printf("DBG1: %d %d %f %f\n", k, n, p, ans);
+            //System.out.printf("DBG1: %d %d %f %f %f %f %f\n", 
+            //                        k, n, p, 
+            //                        nchoosek(n,k), Math.pow(p,k), Math.pow(1-p, n-k),
+            //                        ans);
         }
         else
         {
@@ -537,7 +522,7 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
             ans = ans_1;
         }
 
-        //System.out.printf("DBG3: %f\n", ans);
+        //System.out.printf("DBG3: %d %d %f %f\n", n, k, p, ans);
 
         return ans;
     }
@@ -555,7 +540,8 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
         for (long i = 1; i <= k; i++)
             accum = accum * (n-k+i) / i;
 
-        return accum + 0.5; // avoid rounding error
+        //return accum + 0.5; // avoid rounding error
+        return accum; // avoid rounding error
 
                             /*
         long m = n - k;
@@ -664,6 +650,9 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
     private double prior_alt_frequency = -1.0;
     public void setAlleleFrequencyPrior(double frequency)
     {
+        assert(! Double.isNaN(frequency) );
+        assert(! Double.isInfinite(frequency));
+        if (frequency == 0) { frequency = 1e-5; } // how many epsilons do we need!? This is worrisome.
         this.prior_alt_frequency = frequency;
     }
 
@@ -755,6 +744,7 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
         AlleleFrequencyWalker.binomialProb(1500,2965,0.508065);
         AlleleFrequencyWalker.binomialProb(0,197,0.5);
         AlleleFrequencyWalker.binomialProb(13,197,0.5);
+        AlleleFrequencyWalker.binomialProb(0, 7, 1e-3);
         System.exit(0);
 
         {
diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCallingExperiment.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCallingExperiment.java
index c667ea1d1..b768deb5c 100644
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCallingExperiment.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCallingExperiment.java
@@ -183,11 +183,12 @@ public class PoolCallingExperiment extends LocusWalker<AlleleFrequencyEstimate,
                     {
 	                    for (int j = 0; j <= shallow_calls[i].N; j++)
 	                    {
-	                        //System.out.printf("DBG3: %d %f %d\n", j, shallow_calls[i].posteriors[j], shallow_calls[i].N);
-	                        EM_sum += shallow_calls[i].posteriors[j] * (double)j;
-	                        EM_sum += shallow_calls[i].N;
+                            if (Double.isInfinite(shallow_calls[i].posteriors[j])) { shallow_calls[i].posteriors[j] = -10000; }
+	                        System.out.printf("DBG3: %d %f %d\n", j, shallow_calls[i].posteriors[j], shallow_calls[i].N);
+	                        EM_sum += Math.pow(10,shallow_calls[i].posteriors[j]) * (double)j;
+	                        EM_N   += shallow_calls[i].N;
 	                    }
-                    }
+                    }   
 
 	            }
 	        }
@@ -195,6 +196,8 @@ public class PoolCallingExperiment extends LocusWalker<AlleleFrequencyEstimate,
 	        shallow_calls_fraction_correct = correct_shallow_calls / total_shallow_calls;
             trajectory[iterations+1] = EM_alt_freq;
             likelihood_trajectory[iterations+1] = likelihood/(double)total_shallow_calls;
+
+            System.out.printf("DBGTRAJ %f %f %f %f\n", EM_sum, EM_N, trajectory[iterations], trajectory[iterations+1]);
         }
 
         // 7. Compare to estimation from the pool.
diff --git a/java/src/org/broadinstitute/sting/playground/utils/AlleleFrequencyEstimate.java b/java/src/org/broadinstitute/sting/playground/utils/AlleleFrequencyEstimate.java
index 95020cb1a..1f7de6a86 100755
--- a/java/src/org/broadinstitute/sting/playground/utils/AlleleFrequencyEstimate.java
+++ b/java/src/org/broadinstitute/sting/playground/utils/AlleleFrequencyEstimate.java
@@ -46,7 +46,14 @@ public class AlleleFrequencyEstimate {
         if( Double.isNaN(pBest)) { System.out.printf("pBest is NaN\n"); }
         if( Double.isNaN(pRef)) { System.out.printf("pRef is NaN\n"); }
 
-        if( Double.isInfinite(lodVsRef)) { System.out.printf("lodVsRef is Infinite\n"); }
+        if( Double.isInfinite(lodVsRef)) 
+        { 
+            System.out.printf("lodVsRef is Infinite: %c %s\n", ref, bases); 
+            for (int i = 0; i < posteriors.length; i++)
+            {
+                System.out.printf("POSTERIOR %d %f\n", i, posteriors[i]);
+            }
+        }
         if( Double.isInfinite(lodVsNextBest)) { System.out.printf("lodVsNextBest is Infinite\n"); }
         if( Double.isInfinite(qhat)) { System.out.printf("qhat is Infinite\n"); }
         if( Double.isInfinite(qstar)) { System.out.printf("qstar is Infinite\n"); }