about summary refs log tree commit diff
path: root/src/gemma.cpp
diff options
context:
space:
mode:
authorPjotr Prins2025-11-23 14:09:24 +0100
committerPjotr Prins2025-11-23 14:09:24 +0100
commit8229e4e18dddae7bf570ca37767ff39726eede45 (patch)
tree17986d394dbc3a03da9bebe6eba96119dff571dd /src/gemma.cpp
parentff196955e3ca5bdd671a44f852f42323d3828be2 (diff)
downloadpangemma-8229e4e18dddae7bf570ca37767ff39726eede45.tar.gz
Comments HEAD master
Diffstat (limited to 'src/gemma.cpp')
-rw-r--r--src/gemma.cpp30
1 files changed, 15 insertions, 15 deletions
diff --git a/src/gemma.cpp b/src/gemma.cpp
index 4ef2059..b30c3c3 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2572,19 +2572,19 @@ void GEMMA::BatchRun(PARAM &cPar) {
       cPar.a_mode == M_LMM4 || cPar.a_mode == M_LMM5 || cPar.a_mode == M_LMM9 ||
       cPar.a_mode == M_EIGEN) { // Fit LMM or mvLMM or eigen
     write(cPar.a_mode, "Mode");
-    gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph);
+    gsl_matrix *Y = gsl_matrix_safe_alloc(cPar.ni_test, cPar.n_ph); // Y is ind x phenotypes
     enforce_msg(Y, "allocate Y"); // just to be sure
-    gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);
-    gsl_matrix *B = gsl_matrix_safe_alloc(Y->size2, W->size2); // B is a d by c
+    gsl_matrix *W = gsl_matrix_safe_alloc(Y->size1, cPar.n_cvt);    // W is ind x covariates
+    gsl_matrix *B = gsl_matrix_safe_alloc(Y->size2, W->size2);      // B is a d by c (effect size)
                                                           // matrix
-    gsl_matrix *se_B = gsl_matrix_safe_alloc(Y->size2, W->size2);
-    gsl_matrix *K = gsl_matrix_safe_alloc(Y->size1, Y->size1);
+    gsl_matrix *se_B = gsl_matrix_safe_alloc(Y->size2, W->size2);   // se
+    gsl_matrix *K = gsl_matrix_safe_alloc(Y->size1, Y->size1);      // Kinship aka as G/GRM in the equation ind x ind
     gsl_matrix *U = gsl_matrix_safe_alloc(Y->size1, Y->size1);
-    gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);
-    gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2);
-    gsl_vector *eval = gsl_vector_calloc(Y->size1);
-    gsl_vector *env = gsl_vector_safe_alloc(Y->size1);
-    gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);
+    gsl_matrix *UtW = gsl_matrix_calloc(Y->size1, W->size2);        // Transformed ind x covariates
+    gsl_matrix *UtY = gsl_matrix_calloc(Y->size1, Y->size2);        // Transforemd ind x phenotypes
+    gsl_vector *eval = gsl_vector_calloc(Y->size1);                 // eigen values
+    gsl_vector *env = gsl_vector_safe_alloc(Y->size1);              // E environment
+    gsl_vector *weight = gsl_vector_safe_alloc(Y->size1);           // weights
     debug_msg("Started on LMM");
     assert_issue(is_issue(26), UtY->data[0] == 0.0);
 
@@ -2696,8 +2696,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
       LMM cLmm;
       cLmm.CopyFromParam(cPar);
 
-      gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
-      gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
+      gsl_vector_view Y_col = gsl_matrix_column(Y, 0); // get a single phenotype column
+      gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); // and transformed column
 
       cLmm.AnalyzeGene(U, eval, UtW, &UtY_col.vector, W,
                        &Y_col.vector); // y is the predictor, not the phenotype
@@ -2814,8 +2814,8 @@ void GEMMA::BatchRun(PARAM &cPar) {
 
           // if (is_check_mode()) disable_segfpe(); // disable fast NaN checking for now
 
-          gsl_vector_view Y_col = gsl_matrix_column(Y, 0);
-          gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0);
+          gsl_vector_view Y_col = gsl_matrix_column(Y, 0); // get pheno column
+          gsl_vector_view UtY_col = gsl_matrix_column(UtY, 0); // and its transformed
           write(&Y_col.vector, "Y_col");
 
           if (!cPar.file_bfile.empty()) {
@@ -2840,7 +2840,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
                                     &Y_col.vector, env);
             }
           }
-          cLmm.WriteFiles();
+          cLmm.WriteFiles(); // write output files
           cLmm.CopyToParam(cPar);
         } else {
           debug_msg("fit mvLMM (multiple phenotypes)");