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::LynxJinyangii/RcppTskit/RcppTskit@add-multiple-functions-on-pr-131, 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::LynxJinyangii/RcppTskit/RcppTskit@add-multiple-functions-on-pr-131, 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::LynxJinyangii/RcppTskit/RcppTskit@add-multiple-functions-on-pr-131, 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::LynxJinyangii/RcppTskit/RcppTskit@add-multiple-functions-on-pr-131, any::covr, any::xml2
           needs: coverage
 
       - name: Test coverage

.gitignore

@@ -10,6 +10,7 @@
 .Rhistory
 .Rproj.user/
 .Ruserdata
+.Rproj.user.bak/
 .httr-oauth
 /*.Rcheck/
 /.*tar.gz
@@ -30,3 +31,20 @@ src/Makevars
 src/Makevars.win
 vignettes/*.html
 vignettes/*.pdf
+.idea/
+src/.idea/
+.Rlib/
+*tar.gz
+__pycache__/
+testData/out_phase1_useMacsMut_FALSE/*.trees
+testData/out_msprime_from_macs/*.trees
+dev/testData/out_phase1_useMacsMut_FALSE/*.trees
+dev/testData/out_msprime_from_macs/*.trees
+dev/testData/*.trees
+dev/testData/*_pairs.json
+dev/testData/hotspots*.txt
+dev/testData/*_hotspots_for_window_counts.txt
+dev/testData/out_phase1_useMacsMut_FALSE/*.debug.txt
+dev/testData/out_phase1_useMacsMut_FALSE/*.demography.png
+dev/*.html
+dev/tmp_tree.*

DESCRIPTION

@@ -43,8 +43,7 @@ Depends: R (>= 4.0.0)
 Imports: Rcpp (>= 0.12.7), Rdpack, methods, R6, RcppTskit (>= 0.3.0)
 RdMacros: Rdpack
 LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH, dqrng (>= 0.4.1), RcppTskit
-RoxygenNote: 7.3.3
 Suggests: knitr, rmarkdown, testthat
 VignetteBuilder: knitr
 NeedsCompilation: true
-Config/roxygen2/version: 8.0.0
+Config/roxygen2/version: 8.1.0

NAMESPACE

@@ -146,31 +146,37 @@ import(Rcpp)
 import(RcppTskit)
 importFrom(R6,R6Class)
 importFrom(Rdpack,reprompt)
-importFrom(methods,.hasSlot)
-importFrom(methods,classLabel)
-importFrom(methods,is)
-importFrom(methods,new)
-importFrom(methods,show)
-importFrom(methods,slot)
-importFrom(methods,validObject)
-importFrom(stats,aggregate)
-importFrom(stats,coef)
-importFrom(stats,cov2cor)
-importFrom(stats,dnorm)
-importFrom(stats,lm)
-importFrom(stats,model.matrix)
-importFrom(stats,na.omit)
-importFrom(stats,pnorm)
-importFrom(stats,qgamma)
-importFrom(stats,qnorm)
-importFrom(stats,rbinom)
-importFrom(stats,rnorm)
-importFrom(stats,rpois)
-importFrom(stats,runif)
-importFrom(stats,setNames)
-importFrom(stats,var)
-importFrom(utils,combn)
-importFrom(utils,packageDescription)
-importFrom(utils,read.table)
-importFrom(utils,write.table)
+importFrom(methods,
+  .hasSlot,
+  classLabel,
+  is,
+  new,
+  show,
+  slot,
+  validObject
+)
+importFrom(stats,
+  aggregate,
+  coef,
+  cov2cor,
+  dnorm,
+  lm,
+  model.matrix,
+  na.omit,
+  pnorm,
+  qgamma,
+  qnorm,
+  rbinom,
+  rnorm,
+  rpois,
+  runif,
+  setNames,
+  var
+)
+importFrom(utils,
+  combn,
+  packageDescription,
+  read.table,
+  write.table
+)
 useDynLib(AlphaSimR, .registration = TRUE)

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)
     }
@@ -823,6 +823,14 @@ newPop = function(rawPop,ploidy=NULL,simParam=NULL,nThreads=NULL,...){
 
   output = simParam$finalizePop(output, simParam=simParam, ...)
 
+  ts_forward_attrs <- c("tsForwardSource", "tsForwardPosMeta")
+  for (nm in ts_forward_attrs) {
+    val <- attr(rawPop, nm, exact = TRUE)
+    if (!is.null(val)) {
+      attr(output, nm) <- val
+    }
+  }
+
   return(output)
 }
 
@@ -1534,4 +1542,4 @@ unnameMultiPop = function(x, level = Inf) {
   }
   validObject(mp)
   return(mp)
-}
+}

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()
       private$.varA = numeric()
       private$.varG = numeric()
       private$.varE = numeric()
@@ -187,6 +189,51 @@ SimParam = R6Class(
       invisible(self)
     },
 
+    #' @description Sets genetic-coordinate recombination tracking for the simulation.
+    #' 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 Sets forward tree-sequence tracking on/off.
+    #' When enabled, this also enables genetic-coordinate recombination
+    #' tracking and can initialize the forward TS recorder immediately.
+    #'
+    #' @param isTrackTs should forward tree-sequence tracking be on.
+    #' @param founderPop founder/sample population used to seed TS node maps.
+    #' Required on enable if no recorder is already attached.
+    #' @param force should the check for a running simulation be ignored.
+    setTrackTs = function(isTrackTs, founderPop=NULL, force=FALSE){
+      stopifnot(is.logical(isTrackTs), length(isTrackTs)==1L, !is.na(isTrackTs))
+      if(isTRUE(isTrackTs)){
+        self$setTrackRecGen(TRUE, force=TRUE)
+        hasRecorder = !is.null(.tsForwardGetRecorder(self))
+        if(is.null(founderPop)){
+          if(!hasRecorder){
+            stop("setTrackTs(TRUE) requires founderPop when no recorder is attached", call. = FALSE)
+          }
+        } else {
+          tsForwardInitOnSimParam(self, founderPop=founderPop)
+        }
+      }else{
+        .tsForwardSetRecorder(self, NULL)
+      }
+      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 +264,9 @@ SimParam = R6Class(
       if(private$.isTrackRec){
         private$.recHist = private$.recHist[0:lastId]
       }
+      if(private$.isTrackRecGen){
+        private$.recHistGen = private$.recHistGen[0:lastId]
+      }
       invisible(self)
     },
 
@@ -2139,9 +2189,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,
+                        ploidy){
       nNewInd = lastId-private$.lastId
       stopifnot(nNewInd>0)
       if(length(isDH)==1) isDH = rep(isDH,nNewInd)
@@ -2153,6 +2206,8 @@ SimParam = R6Class(
                 length(isDH)==nNewInd)
       tmp = cbind(mother,father,isDH)
       rownames(tmp) = id
+      keepRecHistGen = isTRUE(getOption("AlphaSimR.tsForwardKeepRecHistGen", FALSE))
+      storeRecHistGen = keepRecHistGen
       if(is.null(hist)){
         newRecHist = vector("list",nNewInd)
         tmpLastHaplo = private$.lastHaplo
@@ -2172,12 +2227,35 @@ SimParam = R6Class(
         names(newRecHist) = id
         private$.recHist = c(private$.recHist, newRecHist)
         private$.lastHaplo = tmpLastHaplo
+
+        if(private$.isTrackRecGen){
+          if(storeRecHistGen){
+            private$.recHistGen = c(private$.recHistGen, newRecHist)
+          }else{
+            emptyRecHistGen = vector("list", nNewInd)
+            names(emptyRecHistGen) = id
+            private$.recHistGen = c(private$.recHistGen, emptyRecHistGen)
+          }
+        }
       }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)
+        if(private$.isTrackRecGen){
+          if(storeRecHistGen){
+            if(is.null(histGen)){
+              histGen = vector("list", nNewInd)
+            }
+            names(histGen) = id
+            private$.recHistGen = c(private$.recHistGen, histGen)
+          }else{
+            emptyRecHistGen = vector("list", nNewInd)
+            names(emptyRecHistGen) = id
+            private$.recHistGen = c(private$.recHistGen, emptyRecHistGen)
+          }
+        }
       }
       private$.pedigree = rbind(private$.pedigree, tmp)
       private$.lastId = lastId
@@ -2292,6 +2370,8 @@ SimParam = R6Class(
     .pedigree="matrix",
     .isTrackRec="logical",
     .recHist="list",
+    .isTrackRecGen = "logical",
+    .recHistGen = "list",
     .varA="numeric",
     .varG="numeric",
     .varE="numeric",
@@ -2735,6 +2815,25 @@ SimParam = R6Class(
       }
     },
 
+    #' @field isTrackRecGen is genetic-coordinate recombination being tracked.
+    isTrackRecGen = function(value){
+      if(missing(value)){
+        private$.isTrackRecGen
+      }else{
+        stop("`$isTrackRecGen` is read only",call.=FALSE)
+      }
+    },
+
+    #' @field recHistGen list of historic genetic-coordinate recombination events.
+    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, returnRecHistGen, recordTsDirect, returnTsSegGen, directAppendFnSEXP) {
+    .Call(`_AlphaSimR_cross`, motherGeno, mother, fatherGeno, father, femaleMap, maleMap, trackRec, motherPloidy, fatherPloidy, v, p, motherCentromere, fatherCentromere, quadProb, nThreads, trackRecGen, returnRecHistGen, recordTsDirect, returnTsSegGen, directAppendFnSEXP)
 }
 
 createDH2 <- function(geno, nDH, genMap, v, p, trackRec, nThreads) {
@@ -399,6 +399,10 @@ MaCS <- function(args, maxSites, inbred, ploidy, nThreads, seed) {
     .Call(`_AlphaSimR_MaCS`, args, maxSites, inbred, ploidy, nThreads, seed)
 }
 
+MaCSTS <- function(args, nChr, inbred, ploidy, nThreads, seed, usePhysicalPositions = FALSE, useMacsMut = FALSE, Nref = NA_real_, expandInbredSamples = TRUE) {
+    .Call(`_AlphaSimR_MaCSTS`, args, nChr, inbred, ploidy, nThreads, seed, usePhysicalPositions, useMacsMut, Nref, expandInbredSamples)
+}
+
 #' @title Summarise `tskit` table collection
 #' @param tc an external pointer to a \code{tsk_table_collection_t} object.
 #' @return A list.
@@ -426,3 +430,27 @@ rtsk_treeseq_get_num_individuals2 <- function(ts) {
     .Call(`_AlphaSimR_rtsk_treeseq_get_num_individuals2`, ts)
 }
 
+tsMutateTableCollection <- function(tc, theta, seed) {
+    invisible(.Call(`_AlphaSimR_tsMutateTableCollection`, tc, theta, seed))
+}
+
+tsFinalizeInbredTableCollection <- function(tc, ploidy) {
+    invisible(.Call(`_AlphaSimR_tsFinalizeInbredTableCollection`, tc, ploidy))
+}
+
+tsForwardNodeTableAddRows <- function(tc, flags, time, population, individual) {
+    .Call(`_AlphaSimR_tsForwardNodeTableAddRows`, tc, flags, time, population, individual)
+}
+
+tsForwardNodeTableAddRowsWithMetadata <- function(tc, flags, time, population, individual, nodeKey) {
+    .Call(`_AlphaSimR_tsForwardNodeTableAddRowsWithMetadata`, tc, flags, time, population, individual, nodeKey)
+}
+
+tsForwardEdgeTableAddRows <- function(tc, left, right, parent, child) {
+    invisible(.Call(`_AlphaSimR_tsForwardEdgeTableAddRows`, tc, left, right, parent, child))
+}
+
+tsForwardSetSampleFlags <- function(tc, samples, clearExisting = TRUE) {
+    invisible(.Call(`_AlphaSimR_tsForwardSetSampleFlags`, tc, samples, clearExisting))
+}
+