Merge branch 'feature/ibdcalccpp' into develop

nickbrazeau · May 5, 2024 · 1690c9f · 1690c9f
2 parents 76a3b1a + ca5fff9
commit 1690c9f
Show file tree

Hide file tree

Showing 2 changed files with 25 additions and 20 deletions.
diff --git a/src/main.cpp b/src/main.cpp
@@ -281,7 +281,7 @@ Rcpp::List get_bvibd_cpp(Rcpp::List args) {
   vector<int> haplo_index = rcpp_to_vector_int(args["haplo_index"]);
   int L = recomb[0][0][0].size();
   int tlim = recomb.size();
-  int coi1 = coi_final[host_index[0]]; // COI from first host 
+  int coi1 = coi_final[host_index[0]]; // COI from first host
   int n = host_index.size(); // number of haplotypes within hosts (cumulative)
 
   // objects for storing results
@@ -321,23 +321,19 @@ Rcpp::List get_bvibd_cpp(Rcpp::List args) {
           sample_haplo[i] = parent_haplo2[t][this_host][this_haplo];
         }
 
-      }
+      } // end sample loop
 
       // check for between host coalescence
       for (int i = 0; i < coi1; ++i) {
-        if (coalesced[i]) {
-          continue;
-        }
         for (int j = coi1; j < n; ++j) {
-          if (coalesced[j]) {
-            continue;
-          }
-          if (sample_host[j] == sample_host[i]) {
+          if (sample_host[i] == sample_host[j]) {
+            if (sample_haplo[i] == sample_haplo[j]) {
               // ibd event
               ibd_target_store[l] = true;
+            }
           }
-        } // end j loop
-      }  // end i loop
+        } // end j coal loop
+      }  // end i coal loop
 
     }  // end t loop
   }  // end l loop

diff --git a/tests/testthat/test-cppIBDcalc.R b/tests/testthat/test-cppIBDcalc.R
@@ -1,7 +1,9 @@
 test_that("test Cpp IBD calculation by hand", {
   # NB, cpp IBD calculation already checked in migration matrix form multiple times
   # here handwriting ARG to confirm IBD
-  swf <- polySimIBD::sim_swf(pos = c(0,1e3,2e3), 
+  # NB, cpp IBD calculation already checked in migration matrix form multiple times
+  # here handwriting ARG to confirm IBD
+  swf <- polySimIBD::sim_swf(pos = c(0,1001,2001,3001), 
                              N = 2, 
                              m = 0.5, 
                              mean_coi = 1, 
@@ -10,18 +12,25 @@ test_that("test Cpp IBD calculation by hand", {
                              tlim = 2)
   # call bvibd
   bvIBD <- polySimIBD::get_bvibd(swf = swf, host_index = c(1,2))
+
+
   # calculate by hand from ARG
   ARG <- polySimIBD::get_arg(swf = swf, host_index = c(1,2))
-  coi <- swf$coi
-  bvtree1 <- ARG[[1]]@c
-  bvtree2 <- ARG[[1]]@c 
-  bvtree3 <- ARG[[1]]@c 
+  (coi <- swf$coi)
+  # always burn first loci
+  bvtree1 <- ARG[[2]]@c
+  bvtree2 <- ARG[[3]]@c 
+  bvtree3 <- ARG[[4]]@c 
 
   # calculate bvtree1 - always look left 
-  h1 <- sum( bvtree1[coi[2]] %in% (coi[1]:coi[2] - 1) ) 
-  h2 <- sum( bvtree2[coi[2]] %in% (coi[1]:coi[2] - 1) ) 
-  h3 <- sum( bvtree2[coi[2]] %in% (coi[1]:coi[2] - 1) ) 
-  handIBD <- sum(h1,h2,h3)/sum(coi)
+  h1 <- any( bvtree1[(coi[1]+1):sum(coi)] %in% (1:coi[1] - 1) ) 
+  h2 <- any( bvtree2[(coi[1]+1):sum(coi)] %in% (1:coi[1] - 1) ) 
+  h3 <- any( bvtree3[(coi[1]+1):sum(coi)] %in% (1:coi[1] - 1) ) 
+  # under SNP vs PSMC (Li/Durbin model) don't know begin and end, so treat as missing info - ie burn first loci
+  wi <- diff(swf$pos)/sum(diff(swf$pos))
+  # weighted average (each loci, denom is 1)
+  handIBD <- sum(c(h1,h2,h3) * wi) # 3 loci, equidistant and weighted thus
+
 
   # confirm
   testthat::expect_equal(bvIBD, handIBD)