19566: Precompute cov glm.

Arvados-DCO-1.1-Signed-off-by: Tom Clegg <tom@curii.com>

diff --git a/glm.go b/glm.go
index 6ae98df928ac02a64d5381ff0a4f5bd8bfb0bfe5..b68bdd1898419ffab3c15d5516c8259ad4646073 100644 (file)
--- a/glm.go
+++ b/glm.go
@@ -35,7 +35,7 @@ func normalize(a []float64) {
 // onehot is the observed outcome, in same order as sampleInfo, but
 // shorter because it only has entries for samples with
 // isTraining==true.
-func pvalueGLM(sampleInfo []sampleInfo, onehot []bool, nPCA int) (p float64) {
+func glmPvalueFunc(sampleInfo []sampleInfo, nPCA int) func(onehot []bool) float64 {
        pcaNames := make([]string, 0, nPCA)
        data := make([][]statmodel.Dtype, 0, nPCA)
        for pca := 0; pca < nPCA; pca++ {
@@ -50,17 +50,11 @@ func pvalueGLM(sampleInfo []sampleInfo, onehot []bool, nPCA int) (p float64) {
                pcaNames = append(pcaNames, fmt.Sprintf("pca%d", pca))
        }
 
-       variant := make([]statmodel.Dtype, 0, len(sampleInfo))
        outcome := make([]statmodel.Dtype, 0, len(sampleInfo))
        constants := make([]statmodel.Dtype, 0, len(sampleInfo))
        row := 0
        for _, si := range sampleInfo {
                if si.isTraining {
-                       if onehot[row] {
-                               variant = append(variant, 1)
-                       } else {
-                               variant = append(variant, 0)
-                       }
                        if si.isCase {
                                outcome = append(outcome, 1)
                        } else {
@@ -70,27 +64,50 @@ func pvalueGLM(sampleInfo []sampleInfo, onehot []bool, nPCA int) (p float64) {
                        row++
                }
        }
-       data = append(data, variant, outcome, constants)
-       dataset := statmodel.NewDataset(data, append(pcaNames, "variant", "outcome", "constants"))
+       data = append([][]statmodel.Dtype{outcome, constants}, data...)
+       names := append([]string{"outcome", "constants"}, pcaNames...)
+       dataset := statmodel.NewDataset(data, names)
 
-       defer func() {
-               if recover() != nil {
-                       // typically "matrix singular or near-singular with condition number +Inf"
-                       p = math.NaN()
-               }
-       }()
-       model, err := glm.NewGLM(dataset, "outcome", append([]string{"constants"}, pcaNames...), glmConfig)
+       model, err := glm.NewGLM(dataset, "outcome", names[1:], glmConfig)
        if err != nil {
-               return math.NaN()
+               log.Printf("%s", err)
+               return func([]bool) float64 { return math.NaN() }
        }
        resultCov := model.Fit()
        logCov := resultCov.LogLike()
-       model, err = glm.NewGLM(dataset, "outcome", append([]string{"constants", "variant"}, pcaNames...), glmConfig)
-       if err != nil {
-               return math.NaN()
+
+       return func(onehot []bool) (p float64) {
+               defer func() {
+                       if recover() != nil {
+                               // typically "matrix singular or near-singular with condition number +Inf"
+                               p = math.NaN()
+                       }
+               }()
+
+               variant := make([]statmodel.Dtype, 0, len(sampleInfo))
+               row := 0
+               for _, si := range sampleInfo {
+                       if si.isTraining {
+                               if onehot[row] {
+                                       variant = append(variant, 1)
+                               } else {
+                                       variant = append(variant, 0)
+                               }
+                               row++
+                       }
+               }
+
+               data := append([][]statmodel.Dtype{data[0], variant}, data[1:]...)
+               names := append([]string{"outcome", "variant"}, names[1:]...)
+               dataset := statmodel.NewDataset(data, names)
+
+               model, err := glm.NewGLM(dataset, "outcome", names[1:], glmConfig)
+               if err != nil {
+                       return math.NaN()
+               }
+               resultComp := model.Fit()
+               logComp := resultComp.LogLike()
+               dist := distuv.ChiSquared{K: 1}
+               return dist.Survival(-2 * (logCov - logComp))
        }
-       resultComp := model.Fit()
-       logComp := resultComp.LogLike()
-       dist := distuv.ChiSquared{K: 1}
-       return dist.Survival(-2 * (logCov - logComp))
 }

diff --git a/glm_test.go b/glm_test.go
index 9761b7b4f0f10388adcfe7cabff32276a474f6d7..62c907f55e656c127ce4573578688f5402a0c660 100644 (file)
--- a/glm_test.go
+++ b/glm_test.go
@@ -174,7 +174,7 @@ func (s *glmSuite) TestPvalueRealDataVsPython(c *check.C) {
                }
        }
 
-       pGo := pvalueGLM(samples, onehot, nPCA)
+       pGo := glmPvalueFunc(samples, nPCA)(onehot)
        c.Logf("pGo = %g", pGo)
 
        var pydata bytes.Buffer
@@ -263,7 +263,7 @@ func (s *glmSuite) TestPvalue(c *check.C) {
                return
        }
 
-       c.Check(pvalueGLM(csv2test(`
+       samples, onehot, npca := csv2test(`
 # case=1, onehot=1, pca1, pca2, pca3
 0, 0, 1, 1.21, 2.37
 0, 0, 2, 1.22, 2.38
@@ -274,7 +274,8 @@ func (s *glmSuite) TestPvalue(c *check.C) {
 1, 1, 1, 1.23, 2.36
 1, 1, 2, 1.22, 2.32
 1, 1, 3, 1.21, 2.31
-`)), check.Equals, 0.002789665435066107)
+`)
+       c.Check(glmPvalueFunc(samples, npca)(onehot), check.Equals, 0.002789665435066107)
 }
 
 var benchSamples, benchOnehot = func() ([]sampleInfo, []bool) {
@@ -300,7 +301,7 @@ var benchSamples, benchOnehot = func() ([]sampleInfo, []bool) {
 
 func (s *glmSuite) BenchmarkPvalue(c *check.C) {
        for i := 0; i < c.N; i++ {
-               p := pvalueGLM(benchSamples, benchOnehot, len(benchSamples[0].pcaComponents))
+               p := glmPvalueFunc(benchSamples, len(benchSamples[0].pcaComponents))(benchOnehot)
                c.Check(p, check.Equals, 0.0)
        }
 }

diff --git a/slicenumpy.go b/slicenumpy.go
index fb3ae95187c5289e4a525caf9fc58e410838f666..1dbe328686f7d84deca7c9f6fe1e0a21678f0251 100644 (file)
--- a/slicenumpy.go
+++ b/slicenumpy.go
@@ -52,6 +52,7 @@ type sliceNumpy struct {
        samples         []sampleInfo
        trainingSet     []int // samples index => training set index, or -1 if not in training set
        trainingSetSize int
+       pvalue          func(onehot []bool) float64
 }
 
 func (cmd *sliceNumpy) RunCommand(prog string, args []string, stdin io.Reader, stdout, stderr io.Writer) int {
@@ -187,6 +188,9 @@ func (cmd *sliceNumpy) run(prog string, args []string, stdin io.Reader, stdout,
                if err != nil {
                        return err
                }
+               if len(cmd.samples[0].pcaComponents) > 0 {
+                       cmd.pvalue = glmPvalueFunc(cmd.samples, cmd.pcaComponents)
+               }
        } else if *caseControlOnly {
                return fmt.Errorf("-case-control-only does not make sense without -samples")
        }
@@ -280,6 +284,9 @@ func (cmd *sliceNumpy) run(prog string, args []string, stdin io.Reader, stdout,
                                cmd.trainingSet[i] = -1
                        }
                }
+               cmd.pvalue = func(onehot []bool) float64 {
+                       return pvalue(onehot, cmd.chi2Cases)
+               }
        }
        if cmd.filter.MinCoverage == 1 {
                // In the generic formula below, floating point
@@ -1602,12 +1609,7 @@ func (cmd *sliceNumpy) tv2homhet(cgs map[string]CompactGenome, maxv tileVariantI
                if col < 4 && !cmd.includeVariant1 {
                        continue
                }
-               var p float64
-               if len(cmd.samples[0].pcaComponents) > 0 {
-                       p = pvalueGLM(cmd.samples, obs[col], cmd.pcaComponents)
-               } else {
-                       p = pvalue(obs[col], cmd.chi2Cases)
-               }
+               p := cmd.pvalue(obs[col])
                if cmd.chi2PValue < 1 && !(p < cmd.chi2PValue) {
                        continue
                }