ts read and write (double cross over breakpoints included)

.github/workflows/R-CMD-check.yaml

@@ -43,7 +43,7 @@ jobs:
 
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::rcmdcheck
+          extra-packages: github::HighlanderLab/RcppTskit/RcppTskit, any::rcmdcheck
           needs: check
 
       - uses: r-lib/actions/check-r-package@v2

.github/workflows/document.yaml

@@ -29,7 +29,7 @@ jobs:
       - name: Install dependencies
         uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::roxygen2
+          extra-packages: github::HighlanderLab/RcppTskit/RcppTskit, any::roxygen2
           needs: roxygen2
 
       - name: Document

.github/workflows/pkgdown.yaml

@@ -30,7 +30,7 @@ jobs:
 
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::pkgdown, local::.
+          extra-packages: github::HighlanderLab/RcppTskit/RcppTskit, any::pkgdown, local::.
           needs: website
 
       - name: Build site

.github/workflows/test-coverage.yaml

@@ -36,7 +36,7 @@ jobs:
 
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::covr, any::xml2
+          extra-packages: github::HighlanderLab/RcppTskit/RcppTskit, any::covr, any::xml2
           needs: coverage
 
       - name: Test coverage

.gitignore

@@ -30,3 +30,5 @@ src/Makevars
 src/Makevars.win
 vignettes/*.html
 vignettes/*.pdf
+.idea/
+src/.idea/

R/Class-Pop.R

@@ -686,7 +686,7 @@ newPop = function(rawPop,ploidy=NULL,simParam=NULL,nThreads=NULL,...){
 .newPop = function(rawPop, id=NULL, mother=NULL, father=NULL,
                    iMother=NULL, iFather=NULL, isDH=NULL,
                    femaleParentPop=NULL, maleParentPop=NULL,
-                   hist=NULL, simParam=NULL, nThreads=NULL,...){
+                   hist=NULL, histGen=NULL, simParam=NULL, nThreads=NULL,...){
   if(is.null(simParam)){
     simParam = get("SP",envir=.GlobalEnv)
   }
@@ -813,7 +813,7 @@ newPop = function(rawPop,ploidy=NULL,simParam=NULL,nThreads=NULL,...){
 
   if(simParam$isTrackPed){
     if(simParam$isTrackRec){
-      simParam$addToRec(lastId,id,iMother,iFather,isDH,hist,output@ploidy)
+      simParam$addToRec(lastId,id,iMother,iFather,isDH,hist,histGen,output@ploidy) #Jinyang modified
     }else{
       simParam$addToPed(lastId,id,iMother,iFather,isDH)
     }

R/Class-SimParam.R

@@ -111,6 +111,8 @@ SimParam = R6Class(
       private$.pedigree = matrix(NA_integer_,nrow=0,ncol=3)
       private$.isTrackRec = FALSE
       private$.recHist = list()
+      private$.isTrackRecGen = FALSE
+      private$.recHistGen = list() # Jinyang added
       private$.varA = numeric()
       private$.varG = numeric()
       private$.varE = numeric()
@@ -187,6 +189,25 @@ SimParam = R6Class(
       invisible(self)
     },
 
+    #' @description Sets genetic-coordinate recombination tracking for the simulation. Jinyang added.
+    #' By default this is turned off. When turned on, it will also turn on pedigree tracking.
+    #'
+    #' @param isTrackRecGen should genetic-coordinate recombination tracking be on.
+    #' @param force should the check for a running simulation be ignored.
+    setTrackRecGen = function(isTrackRecGen, force=FALSE){
+      stopifnot(is.logical(isTrackRecGen))
+      if(!force){
+        private$.isRunning()
+      }
+      private$.isTrackRecGen = isTrackRecGen
+      if(isTrackRecGen){
+        private$.isTrackPed = TRUE
+        private$.isTrackRec = TRUE
+      }
+      invisible(self)
+    },
+
+
     #' @description Resets the internal lastId, the pedigree
     #' and recombination tracking (if in use) to the
     #' supplied lastId. Be careful using this function because
@@ -217,6 +238,10 @@ SimParam = R6Class(
       if(private$.isTrackRec){
         private$.recHist = private$.recHist[0:lastId]
       }
+      # Jinyang added
+      if(private$.isTrackRecGen){
+        private$.recHistGen = private$.recHistGen[0:lastId]
+      }
       invisible(self)
     },
 
@@ -2139,9 +2164,12 @@ SimParam = R6Class(
     #' @param father vector of father iids
     #' @param isDH indicator for DH lines
     #' @param hist new recombination history
+    #' @param histGen new recombination history (genetic coordinate)
     #' @param ploidy ploidy level
     addToRec = function(lastId,id,mother,father,isDH,
-                        hist,ploidy){
+                        hist,
+                        histGen=NULL, # Jinyang added
+                        ploidy){
       nNewInd = lastId-private$.lastId
       stopifnot(nNewInd>0)
       if(length(isDH)==1) isDH = rep(isDH,nNewInd)
@@ -2172,12 +2200,24 @@ SimParam = R6Class(
         names(newRecHist) = id
         private$.recHist = c(private$.recHist, newRecHist)
         private$.lastHaplo = tmpLastHaplo
+
+        # Jinyang added
+        if(private$.isTrackRecGen){
+          #newRecHistGen = vector("list", nNewInd)
+          #names(newRecHistGen) = id
+          private$.recHistGen = c(private$.recHistGen, newRecHist)
+        }
       }else{
         # Add hist to recombination history
         private$.hasHap = c(private$.hasHap, rep(FALSE, nNewInd))
         private$.isFounder = c(private$.isFounder, rep(FALSE, nNewInd))
         names(hist) = id
         private$.recHist = c(private$.recHist, hist)
+        # Jinyang added
+        if(private$.isTrackRecGen){
+          names(histGen) = id
+          private$.recHistGen = c(private$.recHistGen, histGen)
+        }
       }
       private$.pedigree = rbind(private$.pedigree, tmp)
       private$.lastId = lastId
@@ -2292,6 +2332,8 @@ SimParam = R6Class(
     .pedigree="matrix",
     .isTrackRec="logical",
     .recHist="list",
+    .isTrackRecGen = "logical",
+    .recHistGen = "list", #Jinyang added
     .varA="numeric",
     .varG="numeric",
     .varE="numeric",
@@ -2735,6 +2777,25 @@ SimParam = R6Class(
       }
     },
 
+    #' @field isTrackRecGen is recombination being tracked. Jinyang added.
+    isTrackRecGen = function(value){
+      if(missing(value)){
+        private$.isTrackRecGen
+      }else{
+        stop("`$isTrackRecGen` is read only",call.=FALSE)
+      }
+    },
+
+    #' @field recHistGen list of historic recombination events. Jinyang added.
+    recHistGen = function(value){
+      if(missing(value)){
+        private$.recHistGen
+      }else{
+        stop("`$recHistGen` is read only",call.=FALSE)
+      }
+    },
+
+
     #' @field haplotypes list of computed IBD haplotypes
     haplotypes=function(value){
       if(missing(value)){

R/RcppExports.R

@@ -292,8 +292,8 @@ createIbdMat <- function(ibd, chr, nLoci, ploidy, nThreads) {
     .Call(`_AlphaSimR_createIbdMat`, ibd, chr, nLoci, ploidy, nThreads)
 }
 
-cross <- function(motherGeno, mother, fatherGeno, father, femaleMap, maleMap, trackRec, motherPloidy, fatherPloidy, v, p, motherCentromere, fatherCentromere, quadProb, nThreads) {
-    .Call(`_AlphaSimR_cross`, motherGeno, mother, fatherGeno, father, femaleMap, maleMap, trackRec, motherPloidy, fatherPloidy, v, p, motherCentromere, fatherCentromere, quadProb, nThreads)
+cross <- function(motherGeno, mother, fatherGeno, father, femaleMap, maleMap, trackRec, motherPloidy, fatherPloidy, v, p, motherCentromere, fatherCentromere, quadProb, nThreads, trackRecGen) {
+    .Call(`_AlphaSimR_cross`, motherGeno, mother, fatherGeno, father, femaleMap, maleMap, trackRec, motherPloidy, fatherPloidy, v, p, motherCentromere, fatherCentromere, quadProb, nThreads, trackRecGen)
 }
 
 createDH2 <- function(geno, nDH, genMap, v, p, trackRec, nThreads) {