19566: Add constant, check GLM results against Python.
[lightning.git] / slice.go
1 // Copyright (C) The Lightning Authors. All rights reserved.
2 //
3 // SPDX-License-Identifier: AGPL-3.0
4
5 package lightning
6
7 import (
8         "bufio"
9         "encoding/gob"
10         "errors"
11         "flag"
12         "fmt"
13         "io"
14         "net/http"
15         _ "net/http/pprof"
16         "os"
17         "path/filepath"
18         "runtime"
19         "strings"
20         "sync"
21         "sync/atomic"
22
23         "git.arvados.org/arvados.git/sdk/go/arvados"
24         "github.com/klauspost/pgzip"
25         log "github.com/sirupsen/logrus"
26 )
27
28 type slicecmd struct{}
29
30 func (cmd *slicecmd) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
31         var err error
32         defer func() {
33                 if err != nil {
34                         fmt.Fprintf(stderr, "%s\n", err)
35                 }
36         }()
37         flags := flag.NewFlagSet("", flag.ContinueOnError)
38         flags.SetOutput(stderr)
39         pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
40         runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
41         projectUUID := flags.String("project", "", "project `UUID` for output data")
42         priority := flags.Int("priority", 500, "container request priority")
43         outputDir := flags.String("output-dir", "./out", "output `directory`")
44         tagsPerFile := flags.Int("tags-per-file", 50000, "tags per file (nfiles will be ~10M÷x)")
45         err = flags.Parse(args)
46         if err == flag.ErrHelp {
47                 err = nil
48                 return 0
49         } else if err != nil {
50                 return 2
51         }
52         inputDirs := flags.Args()
53         if len(inputDirs) == 0 {
54                 err = errors.New("no input dirs specified")
55                 return 2
56         }
57
58         if *pprof != "" {
59                 go func() {
60                         log.Println(http.ListenAndServe(*pprof, nil))
61                 }()
62         }
63
64         if !*runlocal {
65                 runner := arvadosContainerRunner{
66                         Name:        "lightning slice",
67                         Client:      arvados.NewClientFromEnv(),
68                         ProjectUUID: *projectUUID,
69                         RAM:         500000000000,
70                         VCPUs:       64,
71                         Priority:    *priority,
72                         KeepCache:   2,
73                         APIAccess:   true,
74                 }
75                 for i := range inputDirs {
76                         err = runner.TranslatePaths(&inputDirs[i])
77                         if err != nil {
78                                 return 1
79                         }
80                 }
81                 runner.Args = append([]string{"slice", "-local=true",
82                         "-pprof", ":6060",
83                         "-output-dir", "/mnt/output",
84                 }, inputDirs...)
85                 var output string
86                 output, err = runner.Run()
87                 if err != nil {
88                         return 1
89                 }
90                 fmt.Fprintln(stdout, output)
91                 return 0
92         }
93
94         err = Slice(*tagsPerFile, *outputDir, inputDirs)
95         if err != nil {
96                 return 1
97         }
98         return 0
99 }
100
101 // Read tags+tiles+genomes from srcdir, write to dstdir with (up to)
102 // the specified number of tags per file.
103 func Slice(tagsPerFile int, dstdir string, srcdirs []string) error {
104         var infiles []string
105         for _, srcdir := range srcdirs {
106                 files, err := allFiles(srcdir, matchGobFile)
107                 if err != nil {
108                         return err
109                 }
110                 infiles = append(infiles, files...)
111         }
112         // dirNamespace[dir] is an int in [0,len(dirNamespace)), used below to
113         // namespace variant numbers from different dirs.
114         dirNamespace := map[string]tileVariantID{}
115         for _, path := range infiles {
116                 dir, _ := filepath.Split(path)
117                 if _, ok := dirNamespace[dir]; !ok {
118                         dirNamespace[dir] = tileVariantID(len(dirNamespace))
119                 }
120         }
121         namespaces := tileVariantID(len(dirNamespace))
122
123         var (
124                 tagset     [][]byte
125                 tagsetOnce sync.Once
126                 fs         []*os.File
127                 bufws      []*bufio.Writer
128                 gzws       []*pgzip.Writer
129                 encs       []*gob.Encoder
130
131                 countTileVariants int64
132                 countGenomes      int64
133                 countReferences   int64
134         )
135
136         throttle := throttle{Max: runtime.GOMAXPROCS(0)}
137         for _, infile := range infiles {
138                 infile := infile
139                 throttle.Go(func() error {
140                         f, err := open(infile)
141                         if err != nil {
142                                 return err
143                         }
144                         defer f.Close()
145                         dir, _ := filepath.Split(infile)
146                         namespace := dirNamespace[dir]
147                         log.Printf("reading %s (namespace %d)", infile, namespace)
148                         return DecodeLibrary(f, strings.HasSuffix(infile, ".gz"), func(ent *LibraryEntry) error {
149                                 if err := throttle.Err(); err != nil {
150                                         return err
151                                 }
152                                 if len(ent.TagSet) > 0 {
153                                         tagsetOnce.Do(func() {
154                                                 tagset = ent.TagSet
155                                                 var err error
156                                                 fs, bufws, gzws, encs, err = openOutFiles(dstdir, len(ent.TagSet), tagsPerFile)
157                                                 if err != nil {
158                                                         throttle.Report(err)
159                                                         return
160                                                 }
161                                                 for _, enc := range encs {
162                                                         err = enc.Encode(LibraryEntry{TagSet: tagset})
163                                                         if err != nil {
164                                                                 throttle.Report(err)
165                                                                 return
166                                                         }
167                                                 }
168                                         })
169                                 }
170                                 if err := throttle.Err(); err != nil {
171                                         return err
172                                 }
173                                 atomic.AddInt64(&countTileVariants, int64(len(ent.TileVariants)))
174                                 for _, tv := range ent.TileVariants {
175                                         tv.Variant = tv.Variant*namespaces + namespace
176                                         fileno := 0
177                                         if !tv.Ref {
178                                                 fileno = int(tv.Tag) / tagsPerFile
179                                         }
180                                         err := encs[fileno].Encode(LibraryEntry{
181                                                 TileVariants: []TileVariant{tv},
182                                         })
183                                         if err != nil {
184                                                 return err
185                                         }
186                                 }
187                                 // Here, each output file gets a
188                                 // CompactGenome entry for each
189                                 // genome, even if there are no
190                                 // variants in the relevant range.
191                                 // Easier for downstream code.
192                                 atomic.AddInt64(&countGenomes, int64(len(ent.CompactGenomes)))
193                                 for _, cg := range ent.CompactGenomes {
194                                         for i, v := range cg.Variants {
195                                                 if v > 0 {
196                                                         cg.Variants[i] = v*namespaces + namespace
197                                                 }
198                                         }
199                                         for i, enc := range encs {
200                                                 start := i * tagsPerFile
201                                                 end := start + tagsPerFile
202                                                 if max := len(cg.Variants)/2 + int(cg.StartTag); end > max {
203                                                         end = max
204                                                 }
205                                                 if start < int(cg.StartTag) {
206                                                         start = int(cg.StartTag)
207                                                 }
208                                                 var variants []tileVariantID
209                                                 if start < end {
210                                                         variants = cg.Variants[(start-int(cg.StartTag))*2 : (end-int(cg.StartTag))*2]
211                                                 }
212                                                 err := enc.Encode(LibraryEntry{CompactGenomes: []CompactGenome{{
213                                                         Name:     cg.Name,
214                                                         Variants: variants,
215                                                         StartTag: tagID(start),
216                                                         EndTag:   tagID(start + tagsPerFile),
217                                                 }}})
218                                                 if err != nil {
219                                                         return err
220                                                 }
221                                         }
222                                 }
223                                 // Write all ref seqs to the first
224                                 // slice. Easier for downstream code.
225                                 atomic.AddInt64(&countReferences, int64(len(ent.CompactSequences)))
226                                 if len(ent.CompactSequences) > 0 {
227                                         for _, cs := range ent.CompactSequences {
228                                                 for _, tseq := range cs.TileSequences {
229                                                         for i, libref := range tseq {
230                                                                 tseq[i].Variant = libref.Variant*namespaces + namespace
231                                                         }
232                                                 }
233                                         }
234                                         err := encs[0].Encode(LibraryEntry{CompactSequences: ent.CompactSequences})
235                                         if err != nil {
236                                                 return err
237                                         }
238                                 }
239                                 return nil
240                         })
241                 })
242         }
243         throttle.Wait()
244         if throttle.Err() != nil {
245                 closeOutFiles(fs, bufws, gzws, encs)
246                 return throttle.Err()
247         }
248         defer log.Printf("Total %d tile variants, %d genomes, %d reference sequences", countTileVariants, countGenomes, countReferences)
249         return closeOutFiles(fs, bufws, gzws, encs)
250 }
251
252 func openOutFiles(dstdir string, tags, tagsPerFile int) (fs []*os.File, bufws []*bufio.Writer, gzws []*pgzip.Writer, encs []*gob.Encoder, err error) {
253         nfiles := (tags + tagsPerFile - 1) / tagsPerFile
254         fs = make([]*os.File, nfiles)
255         bufws = make([]*bufio.Writer, nfiles)
256         gzws = make([]*pgzip.Writer, nfiles)
257         encs = make([]*gob.Encoder, nfiles)
258         for i := 0; i*tagsPerFile < tags; i++ {
259                 fs[i], err = os.Create(dstdir + fmt.Sprintf("/library%04d.gob.gz", i))
260                 if err != nil {
261                         return
262                 }
263                 bufws[i] = bufio.NewWriterSize(fs[i], 1<<26)
264                 gzws[i] = pgzip.NewWriter(bufws[i])
265                 encs[i] = gob.NewEncoder(gzws[i])
266         }
267         return
268 }
269
270 func closeOutFiles(fs []*os.File, bufws []*bufio.Writer, gzws []*pgzip.Writer, encs []*gob.Encoder) error {
271         var firstErr error
272         for _, gzw := range gzws {
273                 if gzw != nil {
274                         err := gzw.Close()
275                         if err != nil && firstErr == nil {
276                                 firstErr = err
277                         }
278                 }
279         }
280         for _, bufw := range bufws {
281                 if bufw != nil {
282                         err := bufw.Flush()
283                         if err != nil && firstErr == nil {
284                                 firstErr = err
285                         }
286                 }
287         }
288         for _, f := range fs {
289                 if f != nil {
290                         err := f.Close()
291                         if err != nil && firstErr == nil {
292                                 firstErr = err
293                         }
294                 }
295         }
296         return firstErr
297 }