---
title: "Creating input configurations in gen3sis"
output: rmarkdown::html_vignette
date: "`r format(Sys.time(),  '%d.%m.%Y')`"
vignette: >
  %\VignetteIndexEntry{Creating Configuration}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

The configuration object defines the core processes of the system modeled and define some run settings. It is essential to define the model to be simulated.

- Core processes:  include a customizable initialization, observer, speciation, dispersal, evolution and ecology function. Altogether, these six functions are applied as defined in the simulation engine. The possibility to customize these functions confers the high flexibility of *gen3sis* in terms of including a wide range of theoretical knowledge.
- Settings: include the ecological traits considered in the simulation: whether a random seed is used; allowing simulation reproducibility; start and end times of the simulation; and rules about aborting the simulation, including the maximum global or local species number permitted.
Additionally, the configuration object lists the 

This vignette will guide you trough the config input object. It can be either an R object or a file, that can be easily shared and modified.

Here we will:

[1. **Create an empty config**](#create) using the *write_config_skeleton* function. 

[2. **Explain an example config**](#explore) using the config from [**introduction vignette**](introduction.html).

[3. **Modify a config**](#modify) from within R.




```{r setup}
library(gen3sis)
```

## 1. Create an empty config {#create}


In order to create an empty config, we use the function *write_config_skeleton* and define the output location. For that we have to define where the empty config will be stored.

```{r eval=TRUE}
# get data path
datapath <- system.file(file.path("extdata", "EmptyConfig"), package="gen3sis")
# set config_empty.R file path
config_file_path <- file.path(tempdir(), "config_empty.R")
```

And then create an empty config.

```{r eval=FALSE, message=FALSE}
#writes out a config skeleton
write_config_skeleton(config_file_path)
```
This is how the config_empty.R looks like. We present them by main areas and each part is described briefly. For more information and a hand's on example see next session.


::: {.infobox .caution data-latex="{caution}"}
> ## **Content of config_empty.R**
>
>***Metadata***
>```{r eval=FALSE}
# Version: 1.0
# Author:
# Date:
# Landscape:
# Publications:
# Description:
```
>***Settings***
>```{r eval=FALSE}
# set the random seed for the simulation.
random_seed = NA
# set the starting time step or leave NA to use the earliest/highest time-step.
start_time = NA
# set the end time step or leave as NA to use the latest/lowest time-step (0).
end_time = NA
# maximum total number of species in the simulation before it is aborted.
max_number_of_species = 25000
# maximum number of species within one cell before the simulation is aborted.
max_number_of_coexisting_species = 2500
# a list of traits to include with each species.
# a "dispersal" trait is implicitly added in any case.
trait_names = c("dispersal")
# ranges to scale the input environments with:
# not listed variable:         no scaling takes place
# listed, set to NA:           env. variable will be scaled from [min, max] to [0, 1]
# listed with a given range r: env. variable will be scaled from [r1, r2] to [0, 1]
environmental_ranges = list( )
```
>***Observer***
>```{r eval=FALSE}
end_of_timestep_observer = function(data, vars, config){
  # the list of all species can be found in data$all_species.
  # the current landscape can be found in data$landscape.
}
```
>***Initialization***
>```{r eval=FALSE}
# the initial abundance of a newly colonized cell, both during setup and later when 
# colonizing a cell during the dispersal.
initial_abundance = 1
# place species in the landscape:
create_ancestor_species <- function(landscape, config) {
 stop("create the initial species here")
}
```
>***Core Processes***
>```{r eval=FALSE}
### Dispersal
# the maximum range to consider when calculating the distances from local distance inputs.
max_dispersal <- Inf
# returns n dispersal values.
get_dispersal_values <- function(n, species, landscape, config) {
  stop("calculate dispersal values here")
}
### Speciation
# threshold for genetic distance after which a speciation event takes place.
divergence_threshold = NULL
# factor by which the divergence is increased between geographically isolated population.
# can also be a matrix between the different population clusters.
get_divergence_factor <- function(species, cluster_indices, landscape, config) {
  stop("calculate divergence factor here")
}
### Evolution
# mutate the traits of a species and return the new traits matrix.
apply_evolution <- function(species, cluster_indices, landscape, config) {
  stop("mutate species traits here")
}
### Ecology
# called for every cell with all occurring species, this function calculates abundances and/or who survives for each sites
# returns a vector of abundances.
# set the abundance to 0 for every species supposed to die.
apply_ecology <- function(abundance, traits, environment, config) {
  stop("calculate species abundances and deaths here")
}
```
:::

Note that the main core functions have stop warnings by default, this forces the configuration to be edited according to the rules and assumptions (model) that is to be created. You can do this from scratch or by modifying an existing config.


## 2. Explain example config {#explore}

Here we will go bit by bit of the config used at the [**introduction vignette**](introduction.html) vignette. We will follow the structure presented above.


### **Content of config_.empty_worldcenter.R**

### ***Metadata***
The first section of the config file contains metadata. It informs the user when and by whom the file has been created, which landscape it was made for and informs about associated publications.

>```{r eval=FALSE}
# Version: 1.0
# Author: Oskar Hagen
# Date: 1.7.2020
# Landscape: WorldCenter
# Publications: R-package gen3sis
# Description: Example config used at the introduction vignette and similar to case study global configs in Hagen et al. 2020.
# O. Hagen, B. Flück, F. Fopp, J.S. Cabral, F. Hartig, M. Pontarp, T.F. Rangel, L. Pellissier. gen3sis: The GENeral Engine for Eco-Evolutionary SImulationS on the origins of biodiversity.
```


### ***Settings***
We start our simulation with the first avaiable timestep and end it a the latest available time-step at the landscape, thus *start_time* and *end_time* are set to NA. We limit the maximum total number of species alive in a simulation to 50000 with *max_number_of_species* and the maximum number of species within one site to 10000 with *max_number_of_coexisting_species*. Like so, we stop simulations that could possible go off hand and generate too many species. We define which traits we will consider in our simulation with *traits_names*, in our example, species only have a dispersal ability and a optimum temperature traits. For our model to work, we normalize the environmental values, according to the range informed at *environmental_ranges*.

>```{r eval=FALSE}
random_seed = 666
start_time = NA
end_time = NA
max_number_of_species = 50000
max_number_of_coexisting_species = 10000
trait_names = c("temp",  "dispersal")
environmental_ranges = list("temp" = c(-45, 55), "area"=c(101067, 196949), "prec"=c(1,0.5))
```

### ***Observer***
With the observer function, changes over time in any abiotic or biotic information of the virtual world can be recorded by defining the outputs that are saved at specified time steps, and results can be saved and plotted in real-time as the model runs. 
In our observer function we save the species and plot the richness. Like so, plots are printed as the simulation progresses and also saved inside the plot folder inside output.

>```{r eval=FALSE}
end_of_timestep_observer = function(data, vars, config){
  save_species()
  plot_richness(data$all_species, data$landscape)
}
```

Combined plots can be generated inside the observer function as for example:

>```{r eval=FALSE}
end_of_timestep_observer = function(data, vars, config){
    par(mfrow=c(2,3))
    plot_raster_single(data$landscape$environment[,"temp"], data$landscape, "temp", NA)
    plot_raster_single(data$landscape$environment[,"prec"], data$landscape, "prec", NA)
    plot_raster_single(data$landscape$environment[,"area"], data$landscape, "area", NA)
    plot_richness(data$all_species, data$landscape)
    plot_species_presence(data$all_species[[1]], data$landscape)
    plot(0,type='n',axes=FALSE,ann=FALSE)
    mtext("STATUS",1)
}
```

### ***Initialization***
The initialization function creates the ancestor species at the start of the simulation. Users can define the number of ancestor species, their distribution within the ancient landscape and their trait values.
In our example, the initial abundance of the ancestor species as well as the abundance of species on newly colonized sites is 1. We set one single ancestor species spread across the entire globe, i.e. between longitudes [-180, 180] and latitudes [-90, 90]. We set the **temp** trait of each population of the initial species to the environmental temperature. This means, that each population is adapted to the local conditions, as we will see later on the ecology function. Dispersal is set to one, since it is not ranging. The function *create_ancestor_species* is expected to return a list of new species, i.e. the ancestor(s) one(s).

>```{r eval=FALSE}
initial_abundance = 1
create_ancestor_species <- function(landscape, config) {
  range <- c(-180, 180, -90, 90)
  co <- landscape$coordinates
  selection <- co[, "x"] >= 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"]
  new_species$traits[ , "dispersal"] <- 1
  return(list(new_species))
}
```

### ***Core Processes***

#### ***Dispersal***
The dispersal function iterates over all species populations and determines the connectivity between sites and the colonization of new sites in the grid cell. In our example, the dispersal of all species is the same and is stochastic, following a weibull distribution with shape 6 and scale 999. For different dispersal between species, the dispersal trait should be used. Note that dispersal can be made deterministic by making the function return a fix value. 

>```{r eval=FALSE}
get_dispersal_values <- function(n, species, landscape, config) {
  values <- rweibull(n, shape = 6, scale = 999)
  return(values) }
```

The histogram of 100 draws from the dispersal function used above.

>```{r eval=TRUE, fig.width=6, fig.height=3.2}
n <- 100
hist(rweibull(n, shape = 6, scale = 999), col="black")
```

### ***Speciation***
The speciation iterates over every species separately, registers populations’ geographic occupancy (species range), and determines when geographic isolation between population clusters is higher than a user-defined threshold, triggering a lineage splitting event of cladogenesis. The clustering of occupied sites is based on the species’ dispersal capacity and the landscape connection costs. Over time, disconnected clusters gradually accumulate incompatibility, analogous to genetic differentiation. When the divergence between clusters is above the speciation threshold, those clusters become two or more distinct species, and a divergence matrix reset follows. On the other hand, if geographic clusters come into secondary contact before the speciation occurs, they coalesce and incompatibilities are gradually reduced to zero. In our example, speciation takes place after 12 time-steps of isolation and the divergence increase is the same for all species as indicated by *get_divergence_threshold*.
>```{r eval=FALSE}
divergence_threshold = 12 #this is 2Myrs
get_divergence_factor <- function(species, cluster_indices, landscape, config) {
  return(1) }
```

### ***Evolution***

In the function bellow, clustered populations (exchanging genes) had their trait homogenized and weighted by abundance, meaning that a trait of a population that is doing well in a site, as dictated by the ecology function, will contribute more to the average trait of a cluster. Later, populations mutate based on a normal distribution with standard deviation 0.001, possibly increasing or decreasing species optimum temperature.

>```{r eval=FALSE}
# mutate the traits of a species and return the new traits matrix
apply_evolution <- function(species, cluster_indices, landscape, config) {
  trait_evolutionary_power <- 0.001
  traits <- species[["traits"]]
  cells <- rownames(traits)
  #homogenize trait based on abundance
  for(cluster_index in unique(cluster_indices)){
    cells_cluster <- cells[which(cluster_indices == cluster_index)]
    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)
  }
  #mutations
  mutation_deltas <-rnorm(length(traits[, "temp"]), mean=0, sd=trait_evolutionary_power)
  traits[, "temp"] <- traits[, "temp"] + mutation_deltas
  return(traits)
}
```

###  ***Ecology***
The ecology function determines the abundance or presence of the populations in occupied sites of each species. The function iterates over all occupied sites and updates the species population abundances or presences on the basis of local environmental values, updated co-occurrence patterns and species traits. The function takes as input the species abundance, species trait, species divergence and clusters, and the landscape values. In our example, we calculate abundances in a site based on how close the species is to the site temperature. We scale the values to avoid small numbers and apply a carrying capacity based on aridity and temperature corrected by the area of a site. If abundances are below 1, species are considered extinct, if total abundance in a site is above the carrying capacity, small abundances are removed progressively and randomly distributed across the present species until total abundance is smaller or equal to the carrying capacity.

>```{r eval=FALSE}
apply_ecology <- function(abundance, traits, landscape, config) {
  abundance_scale = 10
  abundance_threshold = 1
  #abundance treashold
  survive <- abundance>=abundance_threshold
  abundance[!survive] <- 0
  abundance <- (( 1-abs( traits[, "temp"] - landscape[, "temp"]))*abundance_scale)*as.numeric(survive)
  #abundance thhreashold
  abundance[abundance<abundance_threshold] <- 0
  k <- ((landscape[,"area"]*(landscape[,"prec"]+0.1)*(landscape[,"temp"]+0.1))*abundance_scale^2)
  total_ab <- sum(abundance)
  subtract <- total_ab-k
  if (subtract > 0) {
    # print(paste("should:", k, "is:", total_ab, "DIFF:", round(subtract,0) ))
    while (total_ab>k){
      alive <- abundance>0
      loose <- sample(1:length(abundance[alive]),1)
      abundance[alive][loose] <- abundance[alive][loose]-1
      total_ab <- sum(abundance)
    }
    #set negative abundances to zero
    abundance[!alive] <- 0
  }
  return(abundance)
}
```

## 3. Modify a config {#modify}

A configuration object can be modified by editing a file or the config object in R. Here we will load the config_worldcenter.R as an R object and modify it's random seed.

```{r eval=TRUE}
# get data path
datapath <- system.file(file.path("extdata", "WorldCenter"), package="gen3sis")
# creates config object from config file
config_object <- create_input_config(file.path(datapath, "config/config_worldcenter.R"))
# modify random seed
config_object$gen3sis$general$random_seed <- 2020
```

Like so, functions and parameters can be conveniently changed and simulated in gen3sis. Remember that *run_simulation* works either with a path to a config file or config object, so we can run the old and the new modified config.

Before using a new config, you can test if it is a valid one. If the function *verify_config* returns TRUE, you are good to go.

```{r}
verify_config(config_object)
```

```{r eval=F}
# run simulation by indicating config file path
sim_old <- run_simulation(config = file.path(datapath, "config/config_worldcenter.R"), 
               landscape = file.path(datapath, "landscape"),
               output_directory=tempdir())

# run simulation by indicating config object
sim_new <- run_simulation(config = config_object, 
               landscape = file.path(datapath, "landscape"),
               output_directory=tempdir())

```