From f8f4578e1c6e1e5a9bb5eeea42f504ee3ee4ed48 Mon Sep 17 00:00:00 2001 From: Oskar Hagen Date: Mon, 12 Dec 2022 15:55:45 +0100 Subject: [PATCH] update --- code/day1_OH.R | 28 ++-- code/source.R | 13 +- data/configs/SA_1d/config_browser.R | 16 +- .../preservation/config_preservation.R | 68 ++++++--- .../preservation/config_preservation_2.R | 141 ++++++++++++++++++ 5 files changed, 218 insertions(+), 48 deletions(-) create mode 100644 data/configs/preservation/config_preservation_2.R diff --git a/code/day1_OH.R b/code/day1_OH.R index 37adb5e..4ca76f6 100644 --- a/code/day1_OH.R +++ b/code/day1_OH.R @@ -37,26 +37,20 @@ install.packages("gen3sis") -library(raster) -library(gen3sis) - ### [] download data -------- ### [] store it into working directory ------------ ### [] source.R ------------ -getwd() -source("C:/Users/am92guke/Documents/iDiv/Teaching/gen3sis workshop iDiv/202211_14-18_WS_MS3BM/code/source.R") -getwd() +source("./source.R") ### [] Reflection -------- ### [] Questions: Sourcing, WD, relative and absolute paths clear? # What does it means, if the bellow does not work? ------------ source(file.path(getwd(), "source.R")) ### [] Questions: Did you noticed the variables in the global environment? - - +#get version of gen3sis. Who is running the dev. version? print(paste("gen3sis version:", packageVersion("gen3sis"))) ### attention! Download the data If using South America!!! @@ -87,7 +81,7 @@ lc$temp[1:10, 1:10] colnames(lc$temp) -plot(rasterFromXYZ(lc$temp[,c("x", "y", "0")])) +plot(rasterFromXYZ(lc$temp[ ,c("x", "y", "0")])) plot(rasterFromXYZ(lc$temp[,c("x", "y", "65")])) # overlay @@ -128,7 +122,7 @@ plot(rl0, col=rgb(0,0,1,0.5,1), add=T) } } - plots(lc$temp) + plots(lc$arid) } @@ -320,7 +314,7 @@ if(run_slow==T){ # time series is constant! t_steps <- get_time_steps(dataset = "Krapp2021") - plot(t_steps) # [] whatfor this plot? ----- + plot(t_steps) # [] what for this plot? ----- # make timesteps at 10k selected_t_steps <- seq(-790000, 0, 10000) @@ -501,7 +495,7 @@ sim <- run_simulation(config = conf_n, # OH NO! Somethings is not working.... -#resolution +#resolution... conf_n$gen3sis$initialization$create_ancestor_species <- function(landscape, config) { # browser() range <- c(-50, -45, -30, 20) @@ -522,12 +516,12 @@ conf_n$gen3sis$initialization$create_ancestor_species <- function(landscape, con # you can see what objects are in the internal environment with ls(). # timestep at the landscape object can help you stop at your convenience: -# stop_time <- "10" -# get_dispersal_values <- function(n, species, landscape, config) { - # if(landscape$timestep== stop_time){browser()} +stop_time <- "64" +get_dispersal_values <- function(n, species, landscape, config) { +if(landscape$timestep== stop_time){browser()} # also found on vars: -# vars <- dynGet("vars", inherits = TRUE) -# if(vars$ti== stop_time){browser()} +vars <- dynGet("vars", inherits = TRUE) +if(vars$ti== stop_time){browser()} # get your error stats # options()$error diff --git a/code/source.R b/code/source.R index eb4eaf8..3c51156 100644 --- a/code/source.R +++ b/code/source.R @@ -8,12 +8,16 @@ ## Author: Hagen (oskar@hagen.bio), Skeels and Rosenbaum ##=======================================================================## -library(gen3sis) -# example course +# Here we have a list of packages we want to install +lop <- c("gen3sis", "raster") #list.of.packages +# And finally we install the missing packages, including their dependency. +new.packages <- lop[!(lop %in% installed.packages()[,"Package"])] +if(length(new.packages)) install.packages(new.packages) +# After the installation process completes, we load all packages. +sapply(lop,require,character=TRUE) # set working directory ----------- -wd <- "C:/Users/am92guke/Documents/iDiv/Teaching/gen3sis workshop iDiv/202211_14-18_WS_MS3BM/code" -setwd(wd) +setwd(getwd()) ### [] create variables for directory relative location ------- # dd= data directory @@ -27,7 +31,6 @@ cat( "# During this prac we won't run the code in the run_slow sections (will take too long!) \n # but so you get an idea of the code, please take a look and see whats going on in these sections \n # keep this as FALSE \n") - run_slow <- FALSE diff --git a/data/configs/SA_1d/config_browser.R b/data/configs/SA_1d/config_browser.R index 75ae06e..93b1c54 100644 --- a/data/configs/SA_1d/config_browser.R +++ b/data/configs/SA_1d/config_browser.R @@ -35,14 +35,14 @@ trait_names = c("temp", "niche_width", "dispersal") # "prec", # a place to inspect the internal state of the simulation and collect additional information if desired end_of_timestep_observer = function(data, vars, config){ - # oldpar <- par(no.readonly = TRUE) - # on.exit(par(oldpar)) - # browser() #first,browse# - # par(mfrow=c(2,2)) + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar)) + #browser() #first,browse# + par(mfrow=c(2,2)) plot_richness(data$all_species, data$landscape) - # plot_species_presence(data$all_species[[1]], data$landscape) - # plot_species_presence(data$all_species[[2]], data$landscape) - # plot_species_presence(data$all_species[[3]], data$landscape) + plot_species_presence(data$all_species[[1]], data$landscape) + plot_species_presence(data$all_species[[2]], data$landscape) + plot_species_presence(data$all_species[[3]], data$landscape) #save # save_species() @@ -90,7 +90,7 @@ create_ancestor_species <- function(landscape, config) { ################# # returns n dispersal values get_dispersal_values <- function(n, species, landscape, config) { - # browser() #first,browse# + browser() #first,browse# # if(landscape$timestep==stop_time){browser()} #second,browse# # hint, look at #78 dispersal.R num_draws <- length(free_cells) * length(presence_spi_ti). Start looking for: get_dispersal_values With Ctrl+Shift+F # at the species level... e.g. plot_species_presence(species, landscape) diff --git a/data/configs/preservation/config_preservation.R b/data/configs/preservation/config_preservation.R index 89c9b3d..ee1ccf9 100644 --- a/data/configs/preservation/config_preservation.R +++ b/data/configs/preservation/config_preservation.R @@ -19,7 +19,7 @@ max_number_of_coexisting_species = 1000 # a list of traits to include with each species # a "dispersion" trait is implictly added in any case #trait_names = c("t_min", "a_min", "competition", "dispersion") -trait_names = c("temp", "dispersal", "foss") # "prec", +trait_names = c("opt_temp", "dispersal", "foss") # "prec", # ranges to scale the input environemts with: # not listed variable: no scaling takes place @@ -30,10 +30,22 @@ environmental_ranges = list() #"temp" = c(-45, 55), "area"=c(6895.094, 196948.4) # a place to inspect the internal state of the simulation and collect additional information if desired end_of_timestep_observer = function(data, vars, config){ #plot - plot_richness(data$all_species, data$landscape) + #plot_richness(data$all_species, data$landscape) #save + + # oldpar <- par(no.readonly = TRUE) + # on.exit(par(oldpar)) + #browser() #first,browse# + #par(mfrow=c(1,2)) + plot_richness(data$all_species, data$landscape) + # plot_species_presence(data$all_species[[1]], data$landscape) + # plot_species_presence(data$all_species[[2]], data$landscape) + # plot_species_presence(data$all_species[[3]], data$landscape) + + save_species() save_landscape() + save_phylogeny() } @@ -45,6 +57,9 @@ end_of_timestep_observer = function(data, vars, config){ initial_abundance = 1 # place species within rectangle, our case entire globe create_ancestor_species <- function(landscape, config) { + #### BROWSER ! ---------- + # browser() + #### BROWSER ! ---------- range <- c(-180, 180, -90, 90) co <- landscape$coordinates selection <- co[, "x"] >= range[1] & @@ -52,13 +67,25 @@ create_ancestor_species <- function(landscape, config) { co[, "y"] >= range[3] & co[, "y"] <= range[4] initial_cells <- rownames(co)[selection] - new_species <- create_species(initial_cells, config) - #set local adaptation to max optimal temp equals local temp - new_species$traits[ , "temp"] <- landscape$environment[,"temp"] - #set fossilization to be randomly selected in a uniform prior - new_species$traits[ , "foss"] <- runif(1,0,1) - new_species$traits[ , "dispersal"] <- 1 - return(list(new_species)) + n_s <- 2 + new_speciez <- list() + opt_temps <- c(17,22) + for (sp_i in 1:n_s){ + new_species <- create_species(initial_cells, config) + #set local adaptation to max optimal temp equals local temp + new_species$traits[ , "opt_temp"] <- opt_temps[sp_i] #landscape$environment[,"temp"] + #set fossilization to be randomly selected in a uniform prior + new_species$traits[ , "foss"] <- runif(1,0,1) + new_species$traits[ , "dispersal"] <- 1 + new_speciez[[sp_i]] <- new_species + } + # new_species <- create_species(initial_cells, config) + # #set local adaptation to max optimal temp equals local temp + # new_species$traits[ , "opt_temp"] <- 22 #landscape$environment[,"temp"] + # #set fossilization to be randomly selected in a uniform prior + # new_species$traits[ , "foss"] <- runif(1,0,1) + # new_species$traits[ , "dispersal"] <- 1 + return(new_speciez) } @@ -67,6 +94,9 @@ create_ancestor_species <- function(landscape, config) { ################# # returns n dispersal values get_dispersal_values <- function(n, species, landscape, config) { + #### BROWSER ! ---------- + browser() + #### BROWSER ! ---------- values <- rweibull(n, shape =3, scale =781) ### VARY return(values) } @@ -100,12 +130,12 @@ apply_evolution <- function(species, cluster_indices, landscape, config) { # hist(traits[cells_cluster, "temp"], main="before") mean_abd <- mean(species$abundance[cells_cluster]) weight_abd <- species$abundance[cells_cluster]/mean_abd - traits[cells_cluster, "temp"] <- mean(traits[cells_cluster, "temp"]*weight_abd) + traits[cells_cluster, "opt_temp"] <- mean(traits[cells_cluster, "opt_temp"]*weight_abd) # hist(traits[cells_cluster, "temp"], main="after") } #mutations - mutation_deltas <-rnorm(length(traits[, "temp"]), mean=0, sd=trait_evolutionary_power) - traits[, "temp"] <- traits[, "temp"] + mutation_deltas + mutation_deltas <-rnorm(length(traits[, "opt_temp"]), mean=0, sd=trait_evolutionary_power) + traits[, "opt_temp"] <- traits[, "opt_temp"] + mutation_deltas mutation_deltas <-rnorm(length(traits[, "foss"]), mean=0, sd=trait_evolutionary_power) traits[, "foss"] <- traits[, "foss"] + mutation_deltas # rang between 0 and 1 @@ -124,17 +154,19 @@ apply_evolution <- function(species, cluster_indices, landscape, config) { # returns a vector of abundances # set the abundance to 0 for every species supposed to die apply_ecology <- function(abundance, traits, landscape, config, abundance_scale = 10, abundance_threshold = 1) { - - fg <- function(x,a,b,c){ - v <- a*exp(-((x-b)^2/(2*c^2))) - return(v) - } + #### BROWSER ! ---------- + # browser() + #### BROWSER ! ---------- + # fg <- function(x,a,b,c){ + # v <- a*exp(-((x-b)^2/(2*c^2))) + # return(v) + # } # # plot(fg(x=seq(0,1,0.01), a=10, b=0.5, c=0.3), type='l') # c ranges from 0.001 to 0.3 (very wide niche) # abline(h=1) # gaussian - abundance <- (abundance_scale*exp(-((traits[, "temp"] - landscape[, "temp"])**2/(21.232**2))))*(landscape[,"prec"]) + abundance <- (abundance_scale*exp(-((traits[, "opt_temp"] - landscape[, "temp"])**2/(21.232**2))))*(landscape[,"arid"]) #abundance thhreashold abundance[abundance= range[1] & + co[, "x"] <= range[2] & + co[, "y"] >= range[3] & + co[, "y"] <= range[4] + initial_cells <- rownames(co)[selection] + new_species <- create_species(initial_cells, config) + #set local adaptation to max optimal temp equals local temp + new_species$traits[ , "temp"] <- landscape$environment[,"temp"] + #set fossilization to be randomly selected in a uniform prior + new_species$traits[ , "foss"] <- runif(1,0,1) + new_species$traits[ , "dispersal"] <- 1 + return(list(new_species)) +} + + +################# +### Dispersal ### +################# +# returns n dispersal values +get_dispersal_values <- function(n, species, landscape, config) { + values <- rweibull(n, shape =3, scale =781) ### VARY + return(values) +} + + +################## +### Speciation ### +################## +# threshold for genetic distance after which a speciation event takes place +divergence_threshold =50 + +# factor by which the divergence is increased between geographicaly isolated population +# can also be a matrix between the different population clusters +get_divergence_factor <- function(species, cluster_indices, landscape, config) { + return(1) +} + + +################ +### Mutation ### +################ +# mutate the traits of a species and return the new traits matrix +apply_evolution <- function(species, cluster_indices, landscape, config) { + trait_evolutionary_power <-0.085 ### VARY + traits <- species[["traits"]] + cells <- rownames(traits) + #homogenize trait based on abundance + for(cluster_index in unique(cluster_indices)){ + # cluster_index <- 1 + cells_cluster <- cells[which(cluster_indices == cluster_index)] + # hist(traits[cells_cluster, "temp"], main="before") + mean_abd <- mean(species$abundance[cells_cluster]) + weight_abd <- species$abundance[cells_cluster]/mean_abd + traits[cells_cluster, "temp"] <- mean(traits[cells_cluster, "temp"]*weight_abd) + # hist(traits[cells_cluster, "temp"], main="after") + } + #mutations + mutation_deltas <-rnorm(length(traits[, "temp"]), mean=0, sd=trait_evolutionary_power) + traits[, "temp"] <- traits[, "temp"] + mutation_deltas + mutation_deltas <-rnorm(length(traits[, "foss"]), mean=0, sd=trait_evolutionary_power) + traits[, "foss"] <- traits[, "foss"] + mutation_deltas + # rang between 0 and 1 + tf <- traits[,"foss"] + tf[tf>1] <- 1 + tf[tf<0] <- 0 + traits[,"foss"] <- tf + return(traits) +} + + +############### +### Ecology ### +############### +# called for every cell with all occuring species, this functin calculates the who survives in the current cells +# returns a vector of abundances +# set the abundance to 0 for every species supposed to die +apply_ecology <- function(abundance, traits, landscape, config, abundance_scale = 10, abundance_threshold = 1) { + + fg <- function(x,a,b,c){ + v <- a*exp(-((x-b)^2/(2*c^2))) + return(v) + } + # + # plot(fg(x=seq(0,1,0.01), a=10, b=0.5, c=0.3), type='l') # c ranges from 0.001 to 0.3 (very wide niche) + # abline(h=1) + + # gaussian + abundance <- (abundance_scale*exp(-((traits[, "temp"] - landscape[, "temp"])**2/(21.232**2))))*(landscape[,"arid"]) + #abundance thhreashold + abundance[abundance