FASTQ directory is gone

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4986 348d0f76-0448-11de-a6fe-93d51630548a

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) {
+expected <- function(maxAN, N, eps, ac1scale = F) {
-	co = 2 * N / ( 1 - eps )
+	scale = 10
-	v = co * ((1 - ps)/(1-eps))^(2 * N - 1)
+	
-	v / sum(v)
+	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)
 

build.xml

@@ -519,10 +519,8 @@
                 <jvmarg value="-Xmx${test.maxmemory}" />
                 <jvmarg value="-Djava.awt.headless=true" />
                 <jvmarg value="-Dpipeline.run=${pipeline.run}" />
-<!--
+<!--                 <jvmarg value="-Xdebug"/> -->
-                <jvmarg value="-Xdebug"/>
+<!--                 <jvmarg value="-Xrunjdwp:transport=dt_socket,server=y,suspend=y,address=5005"/> -->
-                <jvmarg value="-Xrunjdwp:transport=dt_socket,server=y,suspend=y,address=5005"/>
--->
                 <classpath>
                     <path refid="external.dependencies" />
                     <pathelement location="${java.classes}" />

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<ProcessingLoc> getProcessingLocs();
+}

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<ProcessingLoc> processingLocs = new ArrayList<ProcessingLoc>();
+
+    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<ProcessingLoc> getProcessingLocs() {
+        return processingLocs;
+    }
+}