+ if *onehotSingle {
+ fnm := fmt.Sprintf("%s/onehot.npy", *outputDir)
+ err = writeNumpyUint32(fnm, onehot, 2, nzCount)
+ if err != nil {
+ return err
+ }
+ fnm = fmt.Sprintf("%s/onehot-columns.npy", *outputDir)
+ err = writeNumpyInt32(fnm, onehotXref2int32(xrefs), 5, len(xrefs))
+ if err != nil {
+ return err
+ }
+ fnm = fmt.Sprintf("%s/stats.json", *outputDir)
+ j, err := json.Marshal(map[string]interface{}{
+ "pvalueCallCount": cmd.pvalueCallCount,
+ })
+ if err != nil {
+ return err
+ }
+ err = os.WriteFile(fnm, j, 0777)
+ if err != nil {
+ return err
+ }
+ }
+ if *onlyPCA {
+ cols := 0
+ for _, c := range onehot[nzCount:] {
+ if int(c) >= cols {
+ cols = int(c) + 1
+ }
+ }
+ if cols == 0 {
+ return fmt.Errorf("cannot do PCA: one-hot matrix is empty")
+ }
+ log.Printf("have %d one-hot cols", cols)
+ stride := 1
+ for *maxPCATiles > 0 && cols > *maxPCATiles*2 {
+ cols = (cols + 1) / 2
+ stride = stride * 2
+ }
+ if cols%2 == 1 {
+ // we work with pairs of columns
+ cols++
+ }
+ log.Printf("creating full matrix (%d rows) and training matrix (%d rows) with %d cols, stride %d", len(cmd.cgnames), cmd.trainingSetSize, cols, stride)
+ mtxFull := mat.NewDense(len(cmd.cgnames), cols, nil)
+ mtxTrain := mat.NewDense(cmd.trainingSetSize, cols, nil)
+ for i, c := range onehot[nzCount:] {
+ if int(c/2)%stride == 0 {
+ outcol := int(c/2)/stride*2 + int(c)%2
+ mtxFull.Set(int(onehot[i]), outcol, 1)
+ if trainRow := cmd.trainingSet[int(onehot[i])]; trainRow >= 0 {
+ mtxTrain.Set(trainRow, outcol, 1)
+ }
+ }
+ }
+ log.Print("fitting")
+ transformer := nlp.NewPCA(cmd.pcaComponents)
+ transformer.Fit(mtxTrain.T())
+ log.Printf("transforming")
+ pca, err := transformer.Transform(mtxFull.T())
+ if err != nil {
+ return err
+ }
+ pca = pca.T()
+ outrows, outcols := pca.Dims()
+ log.Printf("copying result to numpy output array: %d rows, %d cols", outrows, outcols)
+ out := make([]float64, outrows*outcols)
+ for i := 0; i < outrows; i++ {
+ for j := 0; j < outcols; j++ {
+ out[i*outcols+j] = pca.At(i, j)
+ }
+ }
+ fnm := fmt.Sprintf("%s/pca.npy", *outputDir)
+ log.Printf("writing numpy: %s", fnm)
+ output, err := os.OpenFile(fnm, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
+ if err != nil {
+ return err
+ }
+ npw, err := gonpy.NewWriter(nopCloser{output})
+ if err != nil {
+ return fmt.Errorf("gonpy.NewWriter: %w", err)
+ }
+ npw.Shape = []int{outrows, outcols}
+ err = npw.WriteFloat64(out)
+ if err != nil {
+ return fmt.Errorf("WriteFloat64: %w", err)
+ }
+ err = output.Close()
+ if err != nil {
+ return err
+ }
+ log.Print("done")
+
+ log.Print("copying pca components to sampleInfo")
+ for i := range cmd.samples {
+ cmd.samples[i].pcaComponents = make([]float64, outcols)
+ for c := 0; c < outcols; c++ {
+ cmd.samples[i].pcaComponents[i] = pca.At(i, c)
+ }
+ }
+ log.Print("done")
+
+ err = writeSampleInfo(cmd.samples, *outputDir)
+ if err != nil {
+ return err
+ }
+ }
+ }
+ if !*mergeOutput && !*onehotChunked && !*onehotSingle && !*onlyPCA {
+ tagoffsetFilename := *outputDir + "/chunk-tag-offset.csv"
+ log.Infof("writing tag offsets to %s", tagoffsetFilename)
+ var f *os.File
+ f, err = os.Create(tagoffsetFilename)