about summary refs log tree commit diff
path: root/src/lmm.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/lmm.cpp')
-rw-r--r--src/lmm.cpp50
1 files changed, 18 insertions, 32 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 129cbd6..2b730f3 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -2012,23 +2012,17 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
     // continue;
 
     auto tup = fetch_snp(t); // use the callback
-    auto success = get<0>(tup);
-    if (!success)
+    auto state = get<0>(tup);
+    if (state == SKIP)
       continue;
-    // typedef tuple< bool,MarkerInfo,vector<double> > SnpNameValues2;
-    // auto marker = get<1>(tup);
-    // auto chr = get<2>(tup);
-    // auto mpos = get<3>(tup);
+    if (state == LAST)
+      break; // marker loop because of LOCO
+
     auto markerinfo = get<1>(tup);
     auto gs = get<2>(tup);
 
     markers.push_back(markerinfo);
 
-    // check whether SNP is included in gwasnps (used by LOCO)
-    /*
-    if (process_gwasnps && gwasnps.count(snp) == 0)
-      continue;
-    */
     // drop missing idv and plug mean values for missing geno
     double x_total = 0.0; // sum genotype values to compute x_mean
     uint vpos = 0;        // position in target vector
@@ -2067,7 +2061,6 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
     gsl_vector_safe_memcpy(&Xlarge_col.vector, x);
     c++; // count markers going in
 
-
     if (c % msize == 0) {
       batch_compute(msize,markers);
       markers.clear();
@@ -2212,10 +2205,6 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
   mdb_stat(rtxn, geno_mdb, &stat);
   auto num_markers = stat.ms_entries;
 
-  // fetch_snp is a callback function for every SNP row
-  // returns typedef std::tuple<bool,string,std::vector<double> > SnpNameValues;
-
-  // size_t prev_line = 0;
   auto mdb_fetch = MDB_FIRST;
 
   auto cursor = lmdb::cursor::open(rtxn, geno_mdb);
@@ -2227,21 +2216,16 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
 
     string_view key,value;
 
-    /*
-    while (prev_line <= num) {
-      cursor.get(key, value, MDB_NEXT);
-      prev_line++;
-    }
-    */
-    auto success = cursor.get(key, value, mdb_fetch);
+    auto mdb_success = cursor.get(key, value, mdb_fetch);
     mdb_fetch = MDB_NEXT;
 
     // uint8_t chr;
     vector<double> gs;
     MarkerInfo markerinfo;
 
-    if (success) {
+    if (mdb_success) {
       size_t size = 0;
+      // ---- Depending on the format we get different buffers - currently float and byte are supported:
       if (format == "Gb") {
         size_t num_bytes = value.size() / sizeof(uint8_t);
         assert(num_bytes == ni_total);
@@ -2270,13 +2254,6 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
       // "S>L>L>"
       const uint8_t* data = reinterpret_cast<const uint8_t*>(key.data());
       auto chr = static_cast<uint8_t>(data[1]);
-
-      // printf("%#02x %#02x\n", chr, loco_chr);
-
-      if (is_loco && loco_chr != chr) {
-        return make_tuple(false, markerinfo, gs);
-      }
-
       // Extract big-endian uint32
       // uint32_t rest = static_cast<uint32_t>(data[2]);
       uint32_t pos =  (data[2] << 24) | (data[3] << 16) |
@@ -2285,6 +2262,15 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
       uint32_t num = (data[6] << 24) | (data[7] << 16) |
         (data[8] << 8) | data[9];
 
+      // printf("%#02x %#02x\n", chr, loco_chr);
+
+      if (is_loco && loco_chr != chr) {
+        if (chr > loco_chr)
+            return make_tuple(LAST, MarkerInfo { .name="", .chr=chr, .pos=pos } , gs);
+          else
+            return make_tuple(SKIP, MarkerInfo { .name="", .chr=chr, .pos=pos } , gs);
+      }
+
       string_view value2;
       marker_mdb.get(rtxn,key,value2);
       auto marker = string(value2);
@@ -2299,7 +2285,7 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
 
       // cout << "!!!!" << size << marker << ": af" << maf << " " << gs[0] << "," << gs[1] << "," << gs[2] << "," << gs[3] << endl;
     }
-    return make_tuple(success, markerinfo, gs);
+    return make_tuple(COMPUTE, markerinfo, gs);
   };
   LMM::mdb_analyze(fetch_snp,U,eval,UtW,Uty,W,y,num_markers);