Fix panic on blank line in fasta.

[lightning.git] / tilelib.go

diff --git a/tilelib.go b/tilelib.go
index 339bc215e3e5a82704b03ba81f93d0f8377e70d7..49f5169665bbd19f433cf1d2da63a9313d7fde94 100644 (file)
--- a/tilelib.go
+++ b/tilelib.go
@@ -20,7 +20,7 @@ type tileLibRef struct {
 
 type tileSeq map[string][]tileLibRef
 
-func (tseq tileSeq) Variants() []tileVariantID {
+func (tseq tileSeq) Variants() ([]tileVariantID, int, int) {
        maxtag := 0
        for _, refs := range tseq {
                for _, ref := range refs {
@@ -30,12 +30,18 @@ func (tseq tileSeq) Variants() []tileVariantID {
                }
        }
        vars := make([]tileVariantID, maxtag+1)
+       var kept, dropped int
        for _, refs := range tseq {
                for _, ref := range refs {
+                       if vars[int(ref.tag)] != 0 {
+                               dropped++
+                       } else {
+                               kept++
+                       }
                        vars[int(ref.tag)] = ref.variant
                }
        }
-       return vars
+       return vars, kept, dropped
 }
 
 type tileLibrary struct {
@@ -63,7 +69,7 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq,
                var seqlabel string
                for scanner.Scan() {
                        buf := scanner.Bytes()
-                       if len(buf) == 0 || buf[0] == '>' {
+                       if len(buf) > 0 && buf[0] == '>' {
                                todo <- jobT{seqlabel, fasta}
                                seqlabel, fasta = string(buf[1:]), nil
                                log.Debugf("%s %s reading fasta", filelabel, seqlabel)
@@ -80,6 +86,7 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq,
        }
        found := make([]foundtag, 2000000)
        path := make([]tileLibRef, 2000000)
+       totalFoundTags := 0
        totalPathLen := 0
        skippedSequences := 0
        for job := range todo {
@@ -95,19 +102,21 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq,
                tilelib.taglib.FindAll(job.fasta, func(tagid tagID, pos, taglen int) {
                        found = append(found, foundtag{pos: pos, tagid: tagid, taglen: taglen})
                })
+               totalFoundTags += len(found)
+
+               skipped := 0
                path = path[:0]
                last := foundtag{tagid: -1}
-               for i, f := range found {
-                       if tilelib.skipOOO {
-                               if f.tagid < last.tagid+1 {
-                                       // e.g., last=B, this=A
-                                       continue
-                               }
-                               if f.tagid > last.tagid+1 && i+1 < len(found) && found[i+1].tagid <= f.tagid {
-                                       // e.g., last=A, this=C, next=B
-                                       continue
-                               }
+               if tilelib.skipOOO {
+                       keep := longestIncreasingSubsequence(len(found), func(i int) int { return int(found[i].tagid) })
+                       for i, x := range keep {
+                               found[i] = found[x]
                        }
+                       skipped = len(found) - len(keep)
+                       found = found[:len(keep)]
+               }
+               for i, f := range found {
+                       log.Tracef("%s %s found[%d] == %#v", filelabel, job.label, i, f)
                        if last.taglen > 0 {
                                path = append(path, tilelib.getRef(last.tagid, job.fasta[last.pos:f.pos+f.taglen]))
                        }
@@ -120,10 +129,10 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader) (tileSeq,
                pathcopy := make([]tileLibRef, len(path))
                copy(pathcopy, path)
                ret[job.label] = pathcopy
-               log.Debugf("%s %s tiled with path len %d", filelabel, job.label, len(path))
+               log.Debugf("%s %s tiled with path len %d, skipped %d", filelabel, job.label, len(path), skipped)
                totalPathLen += len(path)
        }
-       log.Printf("%s tiled with total path len %d in %d sequences (skipped %d sequences with '_' in name)", filelabel, totalPathLen, len(ret), skippedSequences)
+       log.Printf("%s tiled with total path len %d in %d sequences (skipped %d sequences with '_' in name, skipped %d out-of-order tags)", filelabel, totalPathLen, len(ret), skippedSequences, totalFoundTags-totalPathLen)
        return ret, scanner.Err()
 }