fix #63, we update the extinction time at the end of the timestep to …

…ti-1 -> any species alive at the end of a timestep will be counted to be alive at the beginning of the next one

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: gen3sis
 Type: Package
 Title: General Engine for Eco-Evolutionary Simulations
-Version: 1.5.6
+Version: 1.5.7
 Authors@R: c(
   person("ETH Zürich", role=c("cph")),
   person("Oskar", "Hagen", email = "[email protected]", role = c("aut", "cre"), comment = "Landscape Ecology, WSL and ETH Zürich, Switzerland"),

R/gen3sis_main.R

@@ -247,7 +247,7 @@ run_simulation <- function(config = NA,
     val$vars$n_sp_alive <- sum(sapply(val$data$all_species, function(sp){ifelse(length(sp[["abundance"]]), 1, 0) }))
     val$vars$n_sp <- length(val$data$all_species)
 
-    val <- update_loop_steps_variable(val$config, val$data, val$vars)
+    val <- update_extinction_times(val$config, val$data, val$vars)
 
     if (verbose>=1){
       cat('step =', ti, ', species alive =', val$vars$n_sp_alive, ', species total =', val$vars$n_sp, '\n')  

R/internal_functions.R

@@ -41,15 +41,15 @@ get_geo_richness <- function(species_list, landscape){
 
 
 #' calculate the individual average traits for all given species
+#' currently deprecated
+# @param species_list a list
 #'
-#' @param species_list a list
-#'
-#' @return a matrix filled with average traits vs all species
-#' @noRd
-get_eco_by_sp <- function(species_list) {
-  averages <- t(sapply(species_list, function(sp) {colMeans(sp[["traits"]])} ))
-  return(invisible(averages))
-}
+# @return a matrix filled with average traits vs all species
+# @noRd
+#get_eco_by_sp <- function(species_list) {
+#  averages <- t(sapply(species_list, function(sp) {colMeans(sp[["traits"]])} ))
+#  return(invisible(averages))
+#}
 
 
 #' saves the current phylogeny in nex format(?)
@@ -72,7 +72,7 @@ write_nex <- function(phy, label="sp", output_file){
     phy$Ancestor[rootphy] <- 0
     phy$Speciation.Type <- as.character(phy$Speciation.Type)
     phy$Speciation.Type[rootphy] <- "COMB"
-    addroot <- c("0","0",phy$Speciation.Time[1],0, "ROOT")
+    addroot <- c("0","0",phy$Speciation.Time[1],val$config$gen3sis$general$end_time - 1, "ROOT")
     names(addroot) <- colnames(phy)
     phy <- rbind(addroot, phy)
     phy$Ancestor <- as.integer(phy$Ancestor)
@@ -90,7 +90,7 @@ write_nex <- function(phy, label="sp", output_file){
 
   if (nrow(phy_no_root)==0){
     #following TreeSimGM and TreeSim convertion
-    if (phy[1,"Extinction.Time"]==0){
+    if (phy[1,"Extinction.Time"]==val$config$gen3sis$general$end_time - 1){
       String_final <- "1" # tree with only root
     } else {
       String_final <- "0" # tree with only root that got extinct
@@ -148,7 +148,7 @@ write_nex <- function(phy, label="sp", output_file){
                            "begin trees;",
                            String_final,
                            "end;",
-                           "[!Tree generated with the General Allele Simulation Model (GaSM), following convention of TreeSim and TreeSimGM r-packages]",
+                           "[!Tree generated with the gen3sis R package, following convention of TreeSim and TreeSimGM r-packages]",
                            sep="\n"))
   } #end test if there is a tree
 

R/simulation.R

@@ -299,32 +299,30 @@ setup_distance_matrix <- function(config, data, vars) {
 }
 
 #-------------------------------------------#
-######## -> UPDATE LOOP STEPS VARIABLE  #######
-#--------------------------------------------#
-#' Updates the geo_richness, turnover and eco_by_sp objects
+######## -> UPDATE Extinction Times  #######
+#-------------------------------------------#
+#' Updates the extinction times 
 #'
 #' @param config the current config
 #' @param data the current data
 #' @param vars the current vars
 #'
 #' @return the general vals(config, data, vars) list
 #' @noRd
-update_loop_steps_variable <- function(config, data, vars) {
+update_extinction_times <- function(config, data, vars) {
   for (sp in data$all_species) {
     if(length(sp[["abundance"]])) {
       data$phy$"Extinction.Time"[as.integer(sp[["id"]])] <- vars$ti
     }
   }
-
+  return(list(config = config, data = data, vars = vars))
   # update turnover
-  data$turnover[toString(vars$ti),] <- c(vars$n_new_sp_ti, sum(data$phy$"Extinction.Time"==(vars$ti+1)), vars$n_sp_alive)
+  #data$turnover[toString(vars$ti),] <- c(vars$n_new_sp_ti, sum(data$phy$"Extinction.Time"==(vars$ti+1)), vars$n_sp_alive)
 
   # update geo_richness
-  data$geo_richness[rownames(data[["landscape"]][["coordinates"]]), as.character(vars$ti)] <- get_geo_richness(data$all_species, data[["landscape"]])
+  #data$geo_richness[rownames(data[["landscape"]][["coordinates"]]), as.character(vars$ti)] <- get_geo_richness(data$all_species, data[["landscape"]])
   # data$geo_richness <- update.geo.richness(geo_sp_ti=data$geo_sp_ti, ti=vars$ti, geo_richness = data$geo_richness )
-  data[["eco_by_sp"]] <- get_eco_by_sp(data$all_species)
-
-  return(list(config = config, data = data, vars = vars))
+  #data[["eco_by_sp"]] <- get_eco_by_sp(data$all_species)
 }
 
 #--------------------------------#

R/skeleton_config.R

@@ -6,6 +6,7 @@
 #' 
 #' @return compiled string
 #' @noRd
+
 skeleton_config <- function(){
 paste0(c('
 ######################################

R/summary_statistics.R

@@ -80,7 +80,11 @@ make_summary <- function(config, data, vars, total_runtime, save_file=TRUE){
   sgen3sis <- list("summary"= list(), "flag"=list(), "system"= list(), "parameters" = list())
 
   #summary
-  sgen3sis$summary <- c(data$summaries, list("richness-final"=data$geo_richness[,c('x','y',as.character(vars$ti))]))
+  final_richness <- cbind(data[["inputs"]][["coordinates"]], rep(0, nrow(data[["inputs"]][["coordinates"]])))
+  colnames(final_richness) <- append(colnames(data[["inputs"]][["coordinates"]]), toString(vars$ti))
+  richness <- get_geo_richness(data$all_species, data$landscape)
+  final_richness[as.integer(names(richness)), ncol(final_richness)] <- richness
+  sgen3sis$summary <- c(data$summaries, list("richness-final"= final_richness))
 
   #flag
   sgen3sis$flag <- vars$flag

inst/extdata/CaseStudy1/reference_saves/phy.nex

@@ -0,0 +1,5 @@
+#NEXUS
+begin trees;
+Tree tree = (((species1:5,species4:1):0,species3:5):0,species2:2):0;
+end;
+[!Tree generated with the gen3sis R package, following convention of TreeSim and TreeSimGM r-packages]

inst/extdata/CaseStudy1/reference_saves/phy.txt

@@ -0,0 +1,5 @@
+"Ancestor"	"Descendent"	"Speciation.Time"	"Extinction.Time"	"Speciation.Type"
+"1"	1	1	5	-1	"ROOT"
+"2"	1	2	5	2	"GENETIC"
+"3"	1	3	5	-1	"GENETIC"
+"4"	1	4	5	3	"GENETIC"

tests/testthat/test-gen3sis_main.R

@@ -11,7 +11,10 @@ test_that("run_simulation works", {
   config$gen3sis$general$start_time <- 5
   #re-set call of observer
   config$gen3sis$general$end_of_timestep_observer <- function(data, vars, config){}
+  tmp_output <- tempdir()
   s <- run_simulation(config = config, 
-                      landscape = file.path(datapath,"landscape"), output_directory = tempdir())
-  expect_true(!is.null(s))
+                      landscape = file.path(datapath,"landscape"), output_directory = tmp_output)
+  ref_summary <- readRDS(file.path(datapath, "reference_saves", "sgen3sis_summary.rds"))
+  expect_true(all.equal(ref_summary, s$summary))
+  expect_true(tools::md5sum(file.path(s$parameters$directories$output, "phy.nex")) == tools::md5sum(file.path(datapath, "reference_saves", "phy.nex")))
 })