16 "git.arvados.org/arvados.git/sdk/go/arvados"
17 "github.com/james-bowman/nlp"
18 "github.com/kshedden/gonpy"
19 log "github.com/sirupsen/logrus"
20 "gonum.org/v1/gonum/mat"
23 type pythonPCA struct{}
25 func (cmd *pythonPCA) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
29 fmt.Fprintf(stderr, "%s\n", err)
32 flags := flag.NewFlagSet("", flag.ContinueOnError)
33 flags.SetOutput(stderr)
34 projectUUID := flags.String("project", "", "project `UUID` for output data")
35 inputFilename := flags.String("i", "-", "input `file`")
36 priority := flags.Int("priority", 500, "container request priority")
37 err = flags.Parse(args)
38 if err == flag.ErrHelp {
41 } else if err != nil {
45 runner := arvadosContainerRunner{
46 Name: "lightning pca",
47 Client: arvados.NewClientFromEnv(),
48 ProjectUUID: *projectUUID,
53 err = runner.TranslatePaths(inputFilename)
57 runner.Prog = "python3"
58 runner.Args = []string{"-c", `import sys
60 from sklearn.decomposition import PCA
61 scipy.save(sys.argv[2], PCA(n_components=4).fit_transform(scipy.load(sys.argv[1])))`, *inputFilename, "/mnt/output/pca.npy"}
63 output, err = runner.Run()
67 fmt.Fprintln(stdout, output+"/pca.npy")
73 func (cmd *goPCA) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
77 fmt.Fprintf(stderr, "%s\n", err)
80 flags := flag.NewFlagSet("", flag.ContinueOnError)
81 flags.SetOutput(stderr)
82 pprof := flags.String("pprof", "", "serve Go profile data at http://`[addr]:port`")
83 runlocal := flags.Bool("local", false, "run on local host (default: run in an arvados container)")
84 projectUUID := flags.String("project", "", "project `UUID` for output data")
85 priority := flags.Int("priority", 500, "container request priority")
86 inputFilename := flags.String("i", "-", "input `file`")
87 outputFilename := flags.String("o", "-", "output `file`")
88 components := flags.Int("components", 4, "number of components")
89 onehot := flags.Bool("one-hot", false, "recode tile variants as one-hot")
90 err = flags.Parse(args)
91 if err == flag.ErrHelp {
94 } else if err != nil {
100 log.Println(http.ListenAndServe(*pprof, nil))
105 if *outputFilename != "-" {
106 err = errors.New("cannot specify output file in container mode: not implemented")
109 runner := arvadosContainerRunner{
110 Name: "lightning pca-go",
111 Client: arvados.NewClientFromEnv(),
112 ProjectUUID: *projectUUID,
113 RAM: 100000000000, // maybe 10x input size?
117 err = runner.TranslatePaths(inputFilename)
121 runner.Args = []string{"pca-go", "-local=true", fmt.Sprintf("-one-hot=%v", *onehot), "-i", *inputFilename, "-o", "/mnt/output/pca.npy"}
123 output, err = runner.Run()
127 fmt.Fprintln(stdout, output+"/pca.npy")
131 var input io.ReadCloser
132 if *inputFilename == "-" {
133 input = ioutil.NopCloser(stdin)
135 input, err = os.Open(*inputFilename)
142 tilelib := tileLibrary{
144 compactGenomes: map[string][]tileVariantID{},
146 err = tilelib.LoadGob(context.Background(), input, strings.HasSuffix(*inputFilename, ".gz"), nil)
155 log.Print("converting cgs to array")
156 data, rows, cols := cgs2array(tilelib.compactGenomes)
158 log.Printf("recode one-hot: %d rows, %d cols", rows, cols)
159 data, _, cols = recodeOnehot(data, cols)
162 log.Printf("creating matrix backed by array: %d rows, %d cols", rows, cols)
163 mtx := array2matrix(rows, cols, data).T()
166 transformer := nlp.NewPCA(*components)
168 log.Printf("transforming")
169 mtx, err = transformer.Transform(mtx)
175 rows, cols = mtx.Dims()
176 log.Printf("copying result to numpy output array: %d rows, %d cols", rows, cols)
177 out := make([]float64, rows*cols)
178 for i := 0; i < rows; i++ {
179 for j := 0; j < cols; j++ {
180 out[i*cols+j] = mtx.At(i, j)
184 var output io.WriteCloser
185 if *outputFilename == "-" {
186 output = nopCloser{stdout}
188 output, err = os.OpenFile(*outputFilename, os.O_CREATE|os.O_WRONLY, 0777)
194 bufw := bufio.NewWriter(output)
195 npw, err := gonpy.NewWriter(nopCloser{bufw})
199 npw.Shape = []int{rows, cols}
200 log.Printf("writing numpy: %d rows, %d cols", rows, cols)
201 npw.WriteFloat64(out)
214 func array2matrix(rows, cols int, data []uint16) mat.Matrix {
215 floatdata := make([]float64, rows*cols)
216 for i, v := range data {
217 floatdata[i] = float64(v)
219 return mat.NewDense(rows, cols, floatdata)