simplify usage

Author: brentp

README.md

@@ -19,12 +19,11 @@ is an attempt to provide an API so that users can write their own tools with lit
 Design
 ------
 
-+ Separate (sorted) data-sources are merged into a single stream (but each interval knows it's original source)
 + data-sources must support the *Relatable* Interface.
 + a user-defined function returns true if 2 *Relatable*'s are related. (only a small number of interval-pairs
-  are sent to be tested--this is handled automatically by *irelate*.).
-+ i.Related() gives access to all of the related intervals (after they are added internally by *irelated*)
-+ the API is a for loop
+  are sent to be tested--this is handled automatically by `IRelate`.).
++ i.Related() gives access to all of the related intervals (after they are added internally by `IRelate`)
++ the "API" is a for loop
 
 Example
 -------
@@ -44,21 +43,21 @@ func CheckRelatedByOverlap(a Relatable, b Relatable) bool {
 // a and b are channels that send Relatables.
 a := ScanToRelatable('intervals.bed', IntervalFromBedLine)
 b := BamToRelatable('some.bam')
-merged := Merge(a, b)
-for interval := range IRelate(merged, CheckRelatedByOverlap, false, 0) {
+for interval := range IRelate(CheckRelatedByOverlap, false, 0, a, b) {
     fmt.Fprintf("%s\t%d\t%d\t%d\n", interval.Chrom(), interval.Start(), interval.End(), len(interval.Related()))
 }
 ```
 
-The 3rd argument to *IRelate* determines if intervals from the same source (file) should be
-related (almost always false). The 4th argument determines the *query* set of intervals. So,
+The 2nd argument to *IRelate* determines if intervals from the same source (file) should be
+related (almost always false). The 3rd argument determines the *query* set of intervals. So,
 only intervals from `a` (the 0th) source will be sent from IRelate. If this is set to -1, then
-all intervals from all sources will be sent.
+all intervals from all sources will be sent. After this, any number of interval streams
+can be passed to `IRelate`
 
 If we only want to count alignments with a given mapping quality, the loop becomes:
 
 ```go
-for interval := range IRelate(merged, CheckRelatedByOverlap, false, 0) {
+for interval := range IRelate(CheckRelatedByOverlap, false, 0, a, b) {
     n := 0
     for _, b := range interval.Related() {
          // cast to a bam to ge the mapping quality.
@@ -82,7 +81,7 @@ CheckRelatedByOverlap) and a for loop, it is easy to create custom applications.
 For example, here is the function to relate all intervals within 2KB:
 ```go
 // CheckRelatedBy2KB returns true if intervals are within 2KB.
-func CheckRelatedByOverlap(a Relatable, b Relatable) bool {
+func CheckRelatedBy2KB(a Relatable, b Relatable) bool {
         distance := uint32(2000)
         // note with distance == 0 this just overlap.
         return (b.Start()-distance < a.End()) && (b.Chrom() == a.Chrom())

bam_test.go

@@ -7,7 +7,7 @@ import (
 func TestBam(t *testing.T) {
 	var g RelatableChannel
 	g = BamToRelatable("data/ex.bam")
-	for i := range IRelate(g, CheckRelatedByOverlap, false, 0) {
+	for i := range IRelate(CheckRelatedByOverlap, false, 0, g) {
 		if len(i.Related()) != 0 {
 			t.Errorf("should have another relation: %d", len(i.Related()))
 

bench.sh

@@ -0,0 +1 @@
+go test -cpu=1,2,4 -bench . -run NOTHING -benchtime 3s -cpuprofile cpu.prof

bench_test.go

@@ -0,0 +1,25 @@
+package irelate
+
+import (
+	"testing"
+)
+
+func benchmarkStreams(nStreams int, b *testing.B) {
+
+	for n := 0; n < b.N; n++ {
+		streams := make([]RelatableChannel, 0)
+		f := "data/test.bed.gz"
+
+		for i := 0; i < nStreams; i++ {
+			streams = append(streams, Streamer(f))
+		}
+
+		for a := range IRelate(CheckRelatedByOverlap, false, 0, streams...) {
+			a.Start()
+		}
+
+	}
+}
+
+func Benchmark2Streams(b *testing.B) { benchmarkStreams(2, b) }
+func Benchmark3Streams(b *testing.B) { benchmarkStreams(3, b) }

data/test.bed.gz

@@ -7,7 +7,7 @@ import (
 func TestGff(t *testing.T) {
 	var g RelatableChannel
 	g = GFFToRelatable("data/ex.gff")
-	for i := range IRelate(g, CheckRelatedByOverlap, true, 0) {
+	for i := range IRelate(CheckRelatedByOverlap, true, 0, g) {
 		if len(i.Related()) != 1 {
 			t.Errorf("should have another relation: %d", len(i.Related()))
 

gff_test.go

@@ -42,9 +42,10 @@ func relate(a Relatable, b Relatable, includeSameSourceRelations bool, relativeT
 
 // CheckRelatedByOverlap returns true if Relatables overlap.
 func CheckRelatedByOverlap(a Relatable, b Relatable) bool {
-	distance := uint32(0)
+	return (b.Start() < a.End()) && (b.Chrom() == a.Chrom())
 	// note with distance == 0 this just overlap.
-	return (b.Start()-distance < a.End()) && (b.Chrom() == a.Chrom())
+	//distance := uint32(0)
+	//return (b.Start()-distance < a.End()) && (b.Chrom() == a.Chrom())
 }
 
 // CheckKNN relates an interval to its k-nearest neighbors.
@@ -97,21 +98,23 @@ func sendSortedRelatables(sendQ *relatableQueue, cache []Relatable, out chan Rel
 // it is assumed that no other `b` Relatables could possibly be related to `a`
 // and so `a` is sent to the returnQ. It is likely that includeSameSourceRelations
 // will only be set to true if one is doing something like a merge.
-func IRelate(stream RelatableChannel,
-	checkRelated func(a Relatable, b Relatable) bool,
+// streams are a variable number of channels that send intervals.
+func IRelate(checkRelated func(a Relatable, b Relatable) bool,
 	includeSameSourceRelations bool,
-	relativeTo int) chan Relatable {
+	relativeTo int,
+	streams ...RelatableChannel) chan Relatable {
 
+	stream := Merge(streams...)
 	out := make(chan Relatable, 64)
 	go func() {
 
 		// use the cache to keep relatables to test against.
-		cache := make([]Relatable, 1, 256)
+		cache := make([]Relatable, 1, 1024)
 		cache[0] = <-stream
 
 		// Use sendQ to make sure we output in sorted order.
 		// We know we can print something when sendQ.minStart < cache.minStart
-		sendQ := make(relatableQueue, 0, 256)
+		sendQ := make(relatableQueue, 0, 1024)
 		nils := 0
 
 		// TODO:if we know the ends are sorted (in addition to start) then we have some additional
@@ -122,6 +125,9 @@ func IRelate(stream RelatableChannel,
 
 			for i, c := range cache {
 				// tried using futures for checkRelated to parallelize... got slower
+				if c == nil {
+					continue
+				}
 				if checkRelated(c, interval) {
 					relate(c, interval, includeSameSourceRelations, relativeTo)
 				} else {
@@ -134,20 +140,20 @@ func IRelate(stream RelatableChannel,
 			}
 
 			// only do this when we have a lot of nils as it's expensive to create a new slice.
-			if nils > 0 {
+			if nils > 1 {
 				// remove nils from the cache (must do this before sending)
 				cache, nils = filter(cache, nils), 0
 				// send the elements from cache in order.
 				// use heuristic to minimize the sending.
-				if len(sendQ) > 128 {
+				if len(sendQ) > 12 {
 					sendSortedRelatables(&sendQ, cache, out)
 				}
 			}
 			cache = append(cache, interval)
 
 		}
 		for _, c := range filter(cache, nils) {
-			if relativeTo == -1 || c.Source() == uint32(relativeTo) {
+			if c.Source() == uint32(relativeTo) || relativeTo == -1 {
 				heap.Push(&sendQ, c)
 			}
 		}
@@ -179,7 +185,6 @@ func Merge(streams ...RelatableChannel) RelatableChannel {
 			interval = heap.Pop(&q).(Relatable)
 			source := interval.Source()
 			ch <- interval
-			// need the case/select stmt here to handle end of each stream
 			// pull the next interval from the same source.
 			next_interval, ok := <-streams[source]
 			if ok {

irelate.go

@@ -47,9 +47,8 @@ func main() {
 
 	buf := bufio.NewWriter(os.Stdout)
 
-	merged := I.Merge(streams...)
 	//for interval := range I.IRelate(merged, I.CheckRelatedByOverlap) {
-	for interval := range I.IRelate(merged, I.CheckRelatedByOverlap, false, 0) {
+	for interval := range I.IRelate(I.CheckRelatedByOverlap, false, 0, streams...) {
 		// for bam output:
 		// bam := *(interval).(*I.Bam)
 		fmt.Fprintf(buf, "%s\t%d\t%d\t%d\n", interval.Chrom(), interval.Start(), interval.End(), len(interval.Related()))

main/main.go

@@ -60,14 +60,11 @@ func ScanToRelatable(file string, fn func(line string) Relatable) RelatableChann
 	scanner, fh := OpenScanFile(file)
 	ch := make(chan Relatable, 32)
 	go func() {
-		var i Relatable
-		defer fh.Close()
 		for scanner.Scan() {
-			line := scanner.Text()
-			i = fn(line)
-			ch <- i
+			ch <- fn(scanner.Text())
 		}
 
+		fh.Close()
 		close(ch)
 	}()
 	return ch