19566: glm one column at a time.
authorTom Clegg <tom@curii.com>
Tue, 29 Nov 2022 16:10:32 +0000 (11:10 -0500)
committerTom Clegg <tom@curii.com>
Tue, 29 Nov 2022 16:10:32 +0000 (11:10 -0500)
Arvados-DCO-1.1-Signed-off-by: Tom Clegg <tom@curii.com>

glm.go
glm_test.go
slicenumpy.go

diff --git a/glm.go b/glm.go
index cc06f39e2ce4104c5fabe24dd1ab121ad7e6cbc2..35ed37891898f989dc2619c13f17e8a74c1d6ae0 100644 (file)
--- a/glm.go
+++ b/glm.go
@@ -18,7 +18,7 @@ var glmConfig = &glm.Config{
        ConcurrentIRLS: 1000,
 }
 
-func pvalueGLM(sampleInfo []sampleInfo, onehotPair [][]bool) float64 {
+func pvalueGLM(sampleInfo []sampleInfo, onehot []bool) float64 {
        nPCA := len(sampleInfo[0].pcaComponents)
        pcaNames := make([]string, 0, nPCA)
        data := make([][]statmodel.Dtype, 0, nPCA)
@@ -37,7 +37,7 @@ func pvalueGLM(sampleInfo []sampleInfo, onehotPair [][]bool) float64 {
        outcome := make([]statmodel.Dtype, 0, len(sampleInfo))
        for row, si := range sampleInfo {
                if si.isTraining {
-                       if onehotPair[0][row] {
+                       if onehot[row] {
                                variant = append(variant, 1)
                        } else {
                                variant = append(variant, 0)
index 875a96c1d5bb4d53b7eae13cafeea4371638855a..8cd06233f2dd352550040d6ad85410b70fba3946 100644 (file)
@@ -21,9 +21,8 @@ func (s *glmSuite) TestPvalue(c *check.C) {
                {id: "sample2", isCase: false, isTraining: true, pcaComponents: []float64{7, -1.2, 2}},
                {id: "sample3", isCase: true, isTraining: true, pcaComponents: []float64{7, -1.2, 2}},
                {id: "sample4", isCase: true, isTraining: true, pcaComponents: []float64{-4, 1.1, -2}},
-       }, [][]bool{
-               {false, false, true, true},
-               {false, false, true, true},
+       }, []bool{
+               false, false, true, true,
        }), check.Equals, 0.09589096738494937)
 
        c.Check(pvalueGLM([]sampleInfo{
@@ -36,9 +35,8 @@ func (s *glmSuite) TestPvalue(c *check.C) {
                {id: "sample7", isCase: true, isTraining: true, pcaComponents: []float64{1, 1.23, 2.36}},
                {id: "sample8", isCase: true, isTraining: true, pcaComponents: []float64{2, 1.22, 2.32}},
                {id: "sample9", isCase: true, isTraining: true, pcaComponents: []float64{3, 1.21, 2.31}},
-       }, [][]bool{
-               {false, false, false, false, false, true, true, true, true},
-               {false, false, false, false, false, true, true, true, true},
+       }, []bool{
+               false, false, false, false, false, true, true, true, true,
        }), check.Equals, 0.001028375654911555)
 
        c.Check(pvalueGLM([]sampleInfo{
@@ -61,16 +59,15 @@ func (s *glmSuite) TestPvalue(c *check.C) {
                {id: "sample17", isCase: true, isTraining: true, pcaComponents: []float64{2.008, 8.09, -2.09}},
                {id: "sample18", isCase: true, isTraining: true, pcaComponents: []float64{2.009, 8.10, -2.10}},
                {id: "sample19", isCase: true, isTraining: true, pcaComponents: []float64{2.010, 8.11, -2.11}},
-       }, [][]bool{
-               {false, false, false, false, false, false, false, false, false, true, true, true, true, true, true, true, true, true, true},
-               {false, false, false, false, false, false, false, false, false, true, true, true, true, true, true, true, true, true, true},
+       }, []bool{
+               false, false, false, false, false, false, false, false, false, true, true, true, true, true, true, true, true, true, true,
        }), check.Equals, 0.9999944849940106)
 }
 
-var benchSamples, benchOnehot = func() ([]sampleInfo, [][]bool) {
+var benchSamples, benchOnehot = func() ([]sampleInfo, []bool) {
        pcaComponents := 10
        samples := []sampleInfo{}
-       onehot := make([][]bool, 2)
+       onehot := []bool{}
        r := make([]float64, pcaComponents)
        for j := 0; j < 10000; j++ {
                for i := 0; i < len(r); i++ {
@@ -83,8 +80,7 @@ var benchSamples, benchOnehot = func() ([]sampleInfo, [][]bool) {
                        isTraining:    true,
                        pcaComponents: append([]float64(nil), r...),
                })
-               onehot[0] = append(onehot[0], j%2 == 0)
-               onehot[1] = append(onehot[1], j%2 == 0)
+               onehot = append(onehot, j%2 == 0)
        }
        return samples, onehot
 }()
index c555635585f2fc3a2f6aaead4227ff3a8aee104f..25fc466c116c4b9cc0a9508378ff33db91247c0b 100644 (file)
@@ -1603,7 +1603,7 @@ func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantI
                }
                var p float64
                if len(cmd.samples[0].pcaComponents) > 0 {
-                       p = pvalueGLM(cmd.samples, obs[col:col+2])
+                       p = pvalueGLM(cmd.samples, obs[col])
                } else {
                        p = pvalue(obs[col], cmd.chi2Cases)
                }