diff --git a/R/privateMutations.R b/R/privateMutations.R
index a23e28247..2fe1f8ab3 100644
--- a/R/privateMutations.R
+++ b/R/privateMutations.R
@@ -16,13 +16,46 @@ f <- function(d, acs) {
legend("topright", legend=lapply(acs, function(x) paste("AC =", x)), fill=cols, title="Sub-population")
}
-expected <- function(ps, N, eps) {
- co = 2 * N / ( 1 - eps )
- v = co * ((1 - ps)/(1-eps))^(2 * N - 1)
- v / sum(v)
+expected <- function(maxAN, N, eps, ac1scale = F) {
+ scale = 10
+
+ f <- function(ps, N) {
+ co = 2 * N / ( 1 - eps )
+ co * ((1 - ps)/(1-eps))^(2 * N - 1)
+ }
+
+ # these are the points that we'll actually show, but we need to do the calculation
+ # special for the AC = 1 given the equation actually fits an infinite population
+ # not a discrete population with max chromosomes
+ ps = 1:maxAN / maxAN
+ v = f(ps, N)
+ v = v / sum(v)
+
+ if ( ac1scale ) {
+ subps = seq(1, maxAN*scale) / (maxAN * scale)
+ #print(subps)
+ subv = f(subps, N)
+ #print(subv)
+ #print(v[1:10])
+ pBelowAC1 = sum(subv[1:scale] / sum(subv))
+ #print(list(pBelowAC1=pBelowAC1, v1=v[1]))
+ v[1] = v[1] + pBelowAC1
+ }
+
+ list(ps = ps, pr = v)
}
f(afs, c(1,2,3,5,10,50))
-x = 1:MAX_AC / 200000
-points(x, expected(x,1000,1e-8),type="l",lty=3,lwd=3)
+
+if ( F ) {
+scale = 100
+ex1 = expected(200000, 1000, 1e-8)
+ex2 = expected(200000*scale, 1000, 1e-8)
+i = 1:(200000*scale) %% scale == 1
+plot(ex2$ps[i], cumsum(ex1$pr), type="l",lty=3,lwd=3, log="x", col="red")
+points(ex2$ps[i], cumsum(ex2$pr)[i], type="l",lty=3,lwd=3, log="x")
+}
+
+ex = expected(200000, 1000, 1e-8, T)
+points(ex$ps, ex$pr, type="l",lty=3,lwd=3)
diff --git a/build.xml b/build.xml
index ee7cb2851..ee6c4debc 100644
--- a/build.xml
+++ b/build.xml
@@ -519,10 +519,8 @@
-
+
+
diff --git a/java/src/org/broadinstitute/sting/utils/threading/GenomeLocProcessingTracker.java b/java/src/org/broadinstitute/sting/utils/threading/GenomeLocProcessingTracker.java
new file mode 100644
index 000000000..4d3507e4c
--- /dev/null
+++ b/java/src/org/broadinstitute/sting/utils/threading/GenomeLocProcessingTracker.java
@@ -0,0 +1,104 @@
+package org.broadinstitute.sting.utils.threading;
+
+import org.broadinstitute.sting.utils.GenomeLoc;
+
+import java.util.List;
+
+/**
+ * Created by IntelliJ IDEA.
+ * User: depristo
+ * Date: 1/13/11
+ * Time: 9:38 AM
+ * To change this template use File | Settings | File Templates.
+ */
+public abstract class GenomeLocProcessingTracker {
+ /**
+ * Information about processing locations and their owners
+ */
+ public static final class ProcessingLoc {
+ GenomeLoc loc;
+ String owner;
+
+ /**
+ * Create an unowned loc
+ * @param loc
+ */
+ public ProcessingLoc(GenomeLoc loc) {
+ this(loc, null);
+
+ }
+
+ /**
+ * Create a loc that's already owned
+ * @param loc
+ * @param owner
+ */
+ public ProcessingLoc(GenomeLoc loc, String owner) {
+ this.loc = loc;
+ this.owner = owner;
+ }
+
+ public GenomeLoc getLoc() {
+ return loc;
+ }
+
+ public String getOwner() {
+ return owner;
+ }
+
+ public void setOwner(String owner) {
+ this.owner = owner;
+ }
+
+ public boolean isOwned() { return owner != null; }
+ }
+
+ // --------------------------------------------------------------------------------
+ //
+ // Code to claim intervals for processing and query for their ownership
+ //
+ // --------------------------------------------------------------------------------
+ /**
+ * Queries the current database if a location is owned. Does not guarantee that the
+ * loc can be owned in a future call, though.
+ *
+ * @param loc
+ * @return
+ */
+ public boolean locIsOwned(GenomeLoc loc) {
+ return findOwner(loc) != null;
+ }
+
+ public ProcessingLoc findOwner(GenomeLoc loc) {
+ for ( ProcessingLoc l : getProcessingLocs() ) {
+ if ( l.getLoc().equals(loc) )
+ return l;
+ }
+
+ return null;
+ }
+
+ /**
+ * The workhorse routine. Attempt to claim processing ownership of loc, with my name.
+ * This is an atomic operation -- other threads / processes will wait until this function
+ * returns. The return result is the ProcessingLoc object describing who owns this
+ * location. If the location isn't already claimed and we now own the location, the pl owner
+ * will be myName. Otherwise, the name of the owner can found in the pl.
+ *
+ * @param loc
+ * @param myName
+ * @return
+ */
+ public abstract ProcessingLoc claimOwnership(GenomeLoc loc, String myName);
+
+ /**
+ * Returns the list of currently owned locations, updating the database as necessary.
+ * DO NOT MODIFY THIS LIST! As with all parallelizing data structures, the list may be
+ * out of date immediately after the call returns, or may be updating on the fly.
+ *
+ * This is really useful for printing, counting, etc. operations that aren't mission critical
+ *
+ * @return
+ */
+ protected abstract List getProcessingLocs();
+}
diff --git a/java/src/org/broadinstitute/sting/utils/threading/SharedMemoryGenomeLocProcessingTracker.java b/java/src/org/broadinstitute/sting/utils/threading/SharedMemoryGenomeLocProcessingTracker.java
new file mode 100644
index 000000000..a57921c4a
--- /dev/null
+++ b/java/src/org/broadinstitute/sting/utils/threading/SharedMemoryGenomeLocProcessingTracker.java
@@ -0,0 +1,31 @@
+package org.broadinstitute.sting.utils.threading;
+
+import org.broadinstitute.sting.utils.GenomeLoc;
+
+import java.util.ArrayList;
+import java.util.List;
+
+/**
+ * For algorithmic testing purposes only. Uses synchronization to keep a consistent
+ * processing list in shared memory.
+ */
+public class SharedMemoryGenomeLocProcessingTracker extends GenomeLocProcessingTracker {
+ List processingLocs = new ArrayList();
+
+ public synchronized ProcessingLoc claimOwnership(GenomeLoc loc, String myName) {
+ // processingLocs is a shared memory synchronized object, and this
+ // method is synchonized, so we can just do our processing
+ ProcessingLoc owner = super.findOwner(loc);
+
+ if ( owner == null ) { // we are unowned
+ owner = new ProcessingLoc(loc, myName);
+ processingLocs.add(owner);
+ }
+
+ return owner;
+ }
+
+ protected synchronized List getProcessingLocs() {
+ return processingLocs;
+ }
+}