Write out K

5 files changed, 17 insertions, 32 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index 6b42763..1cced77 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -1915,7 +1915,16 @@ void GEMMA::BatchRun(PARAM &cPar) {
     cout << "Calculating Relatedness Matrix ... " << endl;
 
     if (cPar.is_mdb) {
-      cPar.MdbCalcKin();
+      time_start = clock();
+      auto K = cPar.MdbCalcKin();
+      cPar.time_G = (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
+      if (cPar.a_mode == M_KIN) {
+        cPar.WriteMatrix(K, "cXX");
+      } else { // M_KIN2
+        cPar.WriteMatrix(K, "sXX");
+      }
+
+      gsl_matrix_safe_free(K);
       return;
     }
 
diff --git a/src/gemma_io.cpp b/src/gemma_io.cpp
index 8fa0eab..baa44b9 100644
--- a/src/gemma_io.cpp
+++ b/src/gemma_io.cpp
@@ -1477,7 +1477,7 @@ void ReadFile_eigenD(const string &file_kd, bool &error, gsl_vector *eval) {
 
 
 // Read bimbam mean genotype file and calculate kinship matrix.
-int mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode, const vector<int> &indicator_snp) {
+gsl_matrix *mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode, const vector<int> &indicator_snp) {
   checkpoint("mdb-kin",file_geno);
 
   /* Create and open the LMDB environment: */
@@ -1537,20 +1537,10 @@ int mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode, const v
     if (t % display_pace == 0) {
       ProgressBar("Reading SNPs", t, num_markers - 1);
     }
-
-    // auto tokens = tokenize_whitespace(line,ni_total+3,"");
-    // auto token_i = tokens.begin();
-    // const char *snp = tokens[0]; // first field
-    // token_i++; // skip snp name
-    // token_i++; // skip nucleotide field
-    // token_i++; // skip nucleotide field
-
     string_view key, value;
     cursor.get(key, value, mdb_fetch);
     mdb_fetch = MDB_NEXT;
-    // cout << " key: " << key << endl ;
     size_t num_floats = value.size() / sizeof(float);
-    // cout << " num: " << num_floats << endl ;
     assert(num_floats == ni_total);
     const float* gs = reinterpret_cast<const float*>(value.data());
     size_t n_miss = 0;
@@ -1559,23 +1549,11 @@ int mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode, const v
     gsl_vector_set_all(geno_miss, 0);
 
     for (size_t i = 0; i < ni_total; ++i) {
-      // enforce_str(token_i != tokens.end(), line + " number of fields");
-      // auto field = *token_i;
-      // auto sfield = std::string(field);
-      // cout << i << ":" << sfield << "," << endl;
-      // if (strncmp(field,"NA",2)==0) { // missing value
-      //  gsl_vector_set(geno_miss, i, 0);
-      //  n_miss++;
-      // } else {
       double d = gs[i];
-      // if (is_strict_mode() && d == 0.0)
-      //   enforce_is_float(std::string(field));  // rule out non NA and non-float fields
       gsl_vector_set(geno, i, d);
       gsl_vector_set(geno_miss, i, 1);
       geno_mean += d;
       geno_var += d * d;
-      // }
-      // token_i++;
     }
 
     geno_mean /= (double)(ni_total - n_miss);
@@ -1617,7 +1595,7 @@ int mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode, const v
     fast_eigen_dgemm("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
     write(matrix_kin,"K updated");
   }
-  ProgressBar("Reading SNPs", 100, 100);
+  ProgressBar("Reading SNPs", num_markers, num_markers);
   cout << endl;
 
   enforce_gsl(gsl_matrix_scale(matrix_kin, 1.0 / (double)ns_test));
@@ -1625,9 +1603,8 @@ int mdb_calc_kin(const string file_geno, bool is_loco, const int k_mode, const v
   gsl_vector_free(geno);
   gsl_vector_free(geno_miss);
   gsl_matrix_free(Xlarge);
-  gsl_matrix_free(matrix_kin);
 
-  return true;
+  return matrix_kin;
 }
 
 // Read bimbam mean genotype file and calculate kinship matrix.
diff --git a/src/gemma_io.h b/src/gemma_io.h
index ad13379..de74ada 100644
--- a/src/gemma_io.h
+++ b/src/gemma_io.h
@@ -87,7 +87,7 @@ void ReadFile_mk(const string &file_mk, vector<int> &indicator_idv,
 void ReadFile_eigenU(const string &file_u, bool &error, gsl_matrix *U);
 void ReadFile_eigenD(const string &file_d, bool &error, gsl_vector *eval);
 
-int mdb_calc_kin(const string file_geno, const bool is_loco, const int mode, const vector<int> &indicator_snp);
+gsl_matrix *mdb_calc_kin(const string file_geno, const bool is_loco, const int mode, const vector<int> &indicator_snp);
 bool BimbamKin(const string file_geno, const set<string> ksnps,
                vector<int> &indicator_snp, const int k_mode,
                const int display_pace, gsl_matrix *matrix_kin,
diff --git a/src/param.cpp b/src/param.cpp
index 3b74e53..1adee9e 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -1293,9 +1293,8 @@ void PARAM::ReadBIMBAMGenotypes(gsl_matrix *UtX, gsl_matrix *K, const bool calc_
   return;
 }
 
-void PARAM::MdbCalcKin() {
-  error = !mdb_calc_kin(file_geno, is_loco(), a_mode - 20, indicator_snp);
-  return;
+gsl_matrix *PARAM::MdbCalcKin() {
+  return mdb_calc_kin(file_geno, is_loco(), a_mode - 20, indicator_snp);
 }
 
 void PARAM::CalcKin(gsl_matrix *matrix_kin) {
diff --git a/src/param.h b/src/param.h
index aa08b5d..7cf4a6c 100644
--- a/src/param.h
+++ b/src/param.h
@@ -347,7 +347,7 @@ public:
   void ProcessCvtPhen();
   void CopyCvtPhen(gsl_matrix *W, gsl_vector *y, size_t flag);
   void CopyCvtPhen(gsl_matrix *W, gsl_matrix *Y, size_t flag);
-  void MdbCalcKin();
+  gsl_matrix *MdbCalcKin();
   void CalcKin(gsl_matrix *matrix_kin);
   void CalcS(const map<string, double> &mapRS2wA,
              const map<string, double> &mapRS2wK, const gsl_matrix *W,