19073: Fix dup tag detection.
[lightning.git] / tilelib.go
index e396e1c9d05ed08c630b0221522c223ae63ec8fa..df6cb73b3be2b9e16418f6d5eb7c5135b7686efe 100644 (file)
@@ -61,6 +61,7 @@ type tileLibrary struct {
        retainNoCalls       bool
        skipOOO             bool
        retainTileSequences bool
+       useDups             bool
 
        taglib         *tagLibrary
        variant        [][][blake2b.Size256]byte
@@ -71,6 +72,8 @@ type tileLibrary struct {
        variants       int64
        // if non-nil, write out any tile variants added while tiling
        encoder *gob.Encoder
+       // set Ref flag when writing new variants to encoder
+       encodeRef bool
 
        onAddTileVariant func(libref tileLibRef, hash [blake2b.Size256]byte, seq []byte) error
        onAddGenome      func(CompactGenome) error
@@ -120,7 +123,7 @@ func (tilelib *tileLibrary) loadTileVariants(tvs []TileVariant, variantmap map[t
        for _, tv := range tvs {
                // Assign a new variant ID (unique across all inputs)
                // for each input variant.
-               variantmap[tileLibRef{Tag: tv.Tag, Variant: tv.Variant}] = tilelib.getRef(tv.Tag, tv.Sequence).Variant
+               variantmap[tileLibRef{Tag: tv.Tag, Variant: tv.Variant}] = tilelib.getRef(tv.Tag, tv.Sequence, tv.Ref).Variant
        }
        return nil
 }
@@ -234,7 +237,7 @@ func (tilelib *tileLibrary) loadCompactSequences(cseqs []CompactSequence, varian
        return nil
 }
 
-func allGobFiles(path string) ([]string, error) {
+func allFiles(path string, re *regexp.Regexp) ([]string, error) {
        var files []string
        f, err := open(path)
        if err != nil {
@@ -249,21 +252,24 @@ func allGobFiles(path string) ([]string, error) {
                if fi.Name() == "." || fi.Name() == ".." {
                        continue
                } else if child := path + "/" + fi.Name(); fi.IsDir() {
-                       add, err := allGobFiles(child)
+                       add, err := allFiles(child, re)
                        if err != nil {
                                return nil, err
                        }
                        files = append(files, add...)
-               } else if strings.HasSuffix(child, ".gob") || strings.HasSuffix(child, ".gob.gz") {
+               } else if re == nil || re.MatchString(child) {
                        files = append(files, child)
                }
        }
+       sort.Strings(files)
        return files, nil
 }
 
+var matchGobFile = regexp.MustCompile(`\.gob(\.gz)?$`)
+
 func (tilelib *tileLibrary) LoadDir(ctx context.Context, path string) error {
        log.Infof("LoadDir: walk dir %s", path)
-       files, err := allGobFiles(path)
+       files, err := allFiles(path, matchGobFile)
        if err != nil {
                return err
        }
@@ -307,7 +313,7 @@ func (tilelib *tileLibrary) LoadDir(ctx context.Context, path string) error {
                                        mtx.Unlock()
                                }
                                for _, tv := range ent.TileVariants {
-                                       variantmap[tileLibRef{Tag: tv.Tag, Variant: tv.Variant}] = tilelib.getRef(tv.Tag, tv.Sequence).Variant
+                                       variantmap[tileLibRef{Tag: tv.Tag, Variant: tv.Variant}] = tilelib.getRef(tv.Tag, tv.Sequence, tv.Ref).Variant
                                }
                                cgs = append(cgs, ent.CompactGenomes...)
                                cseqs = append(cseqs, ent.CompactSequences...)
@@ -538,15 +544,16 @@ func (tilelib *tileLibrary) dump(out io.Writer) {
 }
 
 type importStats struct {
-       InputFile              string
-       InputLabel             string
-       InputLength            int
-       InputCoverage          int
-       PathLength             int
-       DroppedOutOfOrderTiles int
+       InputFile             string
+       InputLabel            string
+       InputLength           int
+       InputCoverage         int
+       PathLength            int
+       DroppedRepeatedTags   int
+       DroppedOutOfOrderTags int
 }
 
-func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader, matchChromosome *regexp.Regexp) (tileSeq, []importStats, error) {
+func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader, matchChromosome *regexp.Regexp, isRef bool) (tileSeq, []importStats, error) {
        ret := tileSeq{}
        type jobT struct {
                label string
@@ -554,6 +561,7 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader, matchChro
        }
        todo := make(chan jobT, 1)
        scanner := bufio.NewScanner(rdr)
+       scanner.Buffer(make([]byte, 256), 1<<29) // 512 MiB, in case fasta does not have line breaks
        go func() {
                defer close(todo)
                var fasta []byte
@@ -599,14 +607,33 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader, matchChro
                        log.Warnf("%s %s no tags found", filelabel, job.label)
                }
 
-               skipped := 0
+               droppedDup := 0
+               if !tilelib.useDups {
+                       // Remove any tags that appeared more than once
+                       dup := map[tagID]bool{}
+                       for _, ft := range found {
+                               _, dup[ft.tagid] = dup[ft.tagid]
+                       }
+                       dst := 0
+                       for _, ft := range found {
+                               if !dup[ft.tagid] {
+                                       found[dst] = ft
+                                       dst++
+                               }
+                       }
+                       droppedDup = len(found) - dst
+                       log.Infof("%s %s dropping %d non-unique tags", filelabel, job.label, droppedDup)
+                       found = found[:dst]
+               }
+
+               droppedOOO := 0
                if tilelib.skipOOO {
-                       log.Infof("%s %s keeping longest increasing subsequence", filelabel, job.label)
                        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)
+                       droppedOOO = len(found) - len(keep)
+                       log.Infof("%s %s dropping %d out-of-order tags", filelabel, job.label, droppedOOO)
                        found = found[:len(keep)]
                }
 
@@ -630,7 +657,7 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader, matchChro
                                } else {
                                        endpos = found[i+1].pos + taglen
                                }
-                               path[i] = tilelib.getRef(f.tagid, job.fasta[startpos:endpos])
+                               path[i] = tilelib.getRef(f.tagid, job.fasta[startpos:endpos], isRef)
                                if countBases(job.fasta[startpos:endpos]) != endpos-startpos {
                                        atomic.AddInt64(&lowquality, 1)
                                }
@@ -645,14 +672,15 @@ func (tilelib *tileLibrary) TileFasta(filelabel string, rdr io.Reader, matchChro
                ret[job.label] = pathcopy
 
                basesIn := countBases(job.fasta)
-               log.Infof("%s %s fasta in %d coverage in %d path len %d low-quality %d skipped-out-of-order %d", filelabel, job.label, len(job.fasta), basesIn, len(path), lowquality, skipped)
+               log.Infof("%s %s fasta in %d coverage in %d path len %d low-quality %d", filelabel, job.label, len(job.fasta), basesIn, len(path), lowquality)
                stats = append(stats, importStats{
-                       InputFile:              filelabel,
-                       InputLabel:             job.label,
-                       InputLength:            len(job.fasta),
-                       InputCoverage:          basesIn,
-                       PathLength:             len(path),
-                       DroppedOutOfOrderTiles: skipped,
+                       InputFile:             filelabel,
+                       InputLabel:            job.label,
+                       InputLength:           len(job.fasta),
+                       InputCoverage:         basesIn,
+                       PathLength:            len(path),
+                       DroppedOutOfOrderTags: droppedOOO,
+                       DroppedRepeatedTags:   droppedDup,
                })
 
                totalPathLen += len(path)
@@ -667,7 +695,7 @@ func (tilelib *tileLibrary) Len() int64 {
 
 // Return a tileLibRef for a tile with the given tag and sequence,
 // adding the sequence to the library if needed.
-func (tilelib *tileLibrary) getRef(tag tagID, seq []byte) tileLibRef {
+func (tilelib *tileLibrary) getRef(tag tagID, seq []byte, usedByRef bool) tileLibRef {
        dropSeq := false
        if !tilelib.retainNoCalls {
                for _, b := range seq {
@@ -788,6 +816,7 @@ func (tilelib *tileLibrary) getRef(tag tagID, seq []byte) tileLibRef {
                tilelib.encoder.Encode(LibraryEntry{
                        TileVariants: []TileVariant{{
                                Tag:      tag,
+                               Ref:      usedByRef,
                                Variant:  variant,
                                Blake2b:  seqhash,
                                Sequence: saveSeq,