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

```{r include = FALSE}
  library(knitr)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  tidy=T,
  fig.align='center',
  tidy.opts = list(width.cutoff=80),
  results='hold'
)
```

```{r tidy=T, out.width='70%', echo=F, fig.cap='*gen3sis*: **gen**eral **e**ngine for **e**co-**e**volutionary **si**mulation**s**', fig.margin=F}
knitr::include_graphics("../inst/logo/gen3sis_logo.png")
```

In this vignette you will be introduced to *gen3sis*, an engine for simulating various spatial eco-evolutionary models. Here will create a virtual planet, biodiversify it and examine virtual species.


In particular, we will go through the following steps: 

[1. **Set-up**](#setup) *gen3sis* package and check if we are set with all the necessary input data.

[2. **Run**](#run) one simple example simulation

[3. **Visualize**](#visualize) the outputs

[4. **Analyze**](#analyze) the outputs

<br>



The first step for running a simulation consists of setting up the folder structure that should contain input data (i.e. config and landscape). 



## 1. Set-up {#setup}
First of all, load the *gen3sis* package. 
```{r setup, message=F}
library(gen3sis)
```


From now on, if you run into problems, refer to *gen3sis* the for the full help documentation.

Let's check the version of the *gen3sis* package we are using

```{r eval=TRUE}
print(paste("gen3sis version:", packageVersion("gen3sis")))
```



All *gen3sis* simulations need a **landscape** object which define the environment over which the simulation will take place and a **configuration** object which defines the mechanisms of speciation, dispersal, trait evolution and ecological interactions. Here we will use a published example[^1]. We will consider the continent of South America with a time span from 65Ma to 0Ma (1 time-step = 1 myr) and a spatial resolution of 1&deg;.

All the data we will use here (i.e. the **landscape** and the **configuration** object) is available in the [simulation repository](https://github.com/project-gen3sis/Simulations), so you need to download the data (i.e. landscape, configs and outputs) and set the correct *datapath*. For more information about these input data, please check the associated metadata.This is stored in the *METADATA.txt* file inside the landscape folder and,  for the configuration, as a comment in the beginning of the config R script.

To access this dataset, define the path to the data contained inside the package.

```{r eval=TRUE, echo=FALSE}
datapath <- system.file(file.path("extdata", "SouthAmerica"), package="gen3sis")
```

The full [landscape](https://github.com/project-Gen3sis/Simulations) can be downloaded and stored locally for running the experiment. An experiment folder should look like this:

|>EXPERIMENT_FOLDER     |
|:-------------------------|
|&nbsp;&nbsp;&nbsp;**>landscape**|
|&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;>distances_local|
|&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;>landscapes.rds|
|&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;>METADATA.txt|
|&nbsp;&nbsp;&nbsp;**>config**|
|&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;>config.R|




### Landscape
The **landscape** is the virtual environment in which the processes of speciation, dispersal, evolution and ecology take place (e.g. temperature and aridity data that might be spatially and temporally dynamic). Landscape objects used in *gen3sis* are generated based on temporal sequences of landscapes in the form of raster files, which are summarized in the form of two classes. 

The first landscape class contains: (i) the geographic coordinates of the landscape site, (ii) the corresponding information on which sites are generally suitable for the organisms modeled (e.g. land or ocean), and (iii) the environmental conditions relevant for the species’ ecology (e.g. temperature and aridity). The landscape may be simplified into a single geographic axis for theoretical experiments, or it may consider realistic configurations aimed at reproducing real local or global landscapes. 

The second landscape class defines the permeability of the landscape for species movement between sites via dispersal. Connection cost between sites is computed for each time step from the gridded landscape data based on harvesine geographic distances modified by a user-defined cost function and stored as sparse distance matrices. Distance matrices, containing the connection costs, are provided at every time step as either: (i) a pre-computed full distance matrix, containing all habitable sites in the landscape (faster simulations but more storage required); or (ii) a local distance matrix, computed from neighboring site distances up to a user-defined range limit (slower simulations but less storage required).


**NOTE** the *gen3sis* landscape contains several files. The landscape.rds file contains the landscape changing over time with its environmental variables. The distances folder stores the information on how the landscape is connected (connections costs) which will interact with the disperal ability of the modeled organisms. To see how to create a landscape object, refer to [**design_landscape vignette**](design_landscape.html) and [**create_input_landscape vignette**](create_input_landscape.html).




```{r eval=T, fig.width=7, fig.height=4, message=F, echo=F, warning=T, results='hide', fig.cap='This figure shows the temperature and aridity of the lanscape used in this vignette at 65, 30 and 0Ma.'}
library(raster)
# read landscapes.rds
landscapes <- readRDS(file.path(datapath, "landscape/landscapes.rds"))
# create rasters for each timestep
temp_65 <- rasterFromXYZ(landscapes$temp[, c(1,2,68)])
temp_30 <- rasterFromXYZ(landscapes$temp[, c(1,2,33)])
temp_0 <- rasterFromXYZ(landscapes$temp[, c(1,2,3)])   

# TEMPERATURE
oldpar <- par(no.readonly = TRUE)
par(mar=c(1,2,3,1), oma=c(0,0,3,0))
layout(matrix(c(1,1,1,4,
                2,2,2,4,
                3,3,3,4), ncol=3))
maxtemp <- 33
mintemp <- -27
temp_breaks <- seq(mintemp, maxtemp, by=1)
temp_colors <- rev(heat.colors(length(temp_breaks)-1))
image(temp_65, col=temp_colors, breaks=temp_breaks, main='65 Ma')
image(temp_30, col=temp_colors, breaks=temp_breaks, main='30 Ma')
title(main=expression('Temperature [\u00B0C]'), line=1, outer=T, cex.main=2)
image(temp_0, col=temp_colors, breaks=temp_breaks, main='0 Ma')

plot.new()
legend_df <- as.data.frame(cbind(seq(0, length(temp_breaks)-1, length.out=(length(temp_breaks))), rep(0.25, (length(temp_breaks))), temp_breaks))
legend_image <- rasterFromXYZ(legend_df, res=0.01)
plot(legend_image, legend.only=T, col=temp_colors, horizontal=T, smallplot=c(0.2, 0.8, 0.55, 0.7), 
   axis.args=list(at=seq(mintemp, maxtemp, 5),labels=seq(mintemp, maxtemp, 5)))

# ARIDITY
arid_300 <- rasterFromXYZ(landscapes$arid[, c(1,2,68)])
arid_65 <- rasterFromXYZ(landscapes$arid[, c(1,2,33)])
arid_0 <- rasterFromXYZ(landscapes$arid[, c(1,2,3)])   

par(mar=c(1,2,3,1), oma=c(0,0,3,0))
layout(matrix(c(1,1,1,4,
                2,2,2,4,
                3,3,3,4), ncol=3))
maxarid <- 1
minarid <- 0
arid_breaks <- seq(minarid, maxarid, by=0.01)
arid_colors <- colorRampPalette(c('grey95', 'peru'))(length(arid_breaks)-1)
image(arid_300, col=arid_colors, breaks=arid_breaks, main='65 Ma')
image(arid_65, col=arid_colors, breaks=arid_breaks, main='30 Ma')
title(main=expression('Aridity [index]'), line=1, outer=T, cex.main=2)
image(arid_0, col=arid_colors, breaks=arid_breaks, main='0 Ma')

plot.new()
legend_df <- as.data.frame(cbind(seq(0, length(arid_breaks)-1, length.out=(length(arid_breaks))), rep(0.25, (length(arid_breaks))), arid_breaks))
legend_image <- rasterFromXYZ(legend_df, res=0.01)
plot(legend_image, legend.only=T, col=arid_colors, horizontal=T, smallplot=c(0.2, 0.8, 0.55, 0.7), 
axis.args=list(at=seq(minarid, maxarid, 0.2),labels=seq(minarid, maxarid, 0.2)))
par(oldpar)
# # AREA
# area_300 <- rasterFromXYZ(landscapes$area[, c(1,2,68)])
# area_65 <- rasterFromXYZ(landscapes$area[, c(1,2,38)])
# area_0 <- rasterFromXYZ(landscapes$area[, c(1,2,3)])   
# 
# par(mar=c(1,2,3,1), oma=c(0,0,3,0))
# layout(matrix(c(1,1,1,4,
#                 2,2,2,4,
#                 3,3,3,4), ncol=3))
# maxarea <- 197000
# minarea <- 151000
# area_breaks <- seq(minarea, maxarea, by=1000)
# area_colors <- rev(gray(seq(0.05, 0.9, length.out=length(area_breaks)-1)))
# image(area_300, col=area_colors, breaks=area_breaks, main='65 Ma')
# image(area_65, col=area_colors, breaks=area_breaks, main='30 Ma')
# title(main=expression('Area [' ~km^2~ ']'), line=1, outer=T, cex.main=2)
# image(area_0, col=area_colors, breaks=area_breaks, main='0 Ma')
#  
# plot.new()
# legend_df <- as.data.frame(cbind(seq(0, length(area_breaks)-1, length.out=(length(area_breaks))), rep(0.25, (length(area_breaks))), area_breaks))
# legend_image <- rasterFromXYZ(legend_df, res=0.01)
# plot(legend_image, legend.only=T, col=area_colors, horizontal=T, smallplot=c(0.2, 0.8, 0.55, 0.7), 
# axis.args=list(at=seq(minarea, maxarea, 7500),labels=seq(minarea, maxarea, 7500)))
``` 
   

### Configuration
The **configuration** are the rules of our virtual planet, it is where we define all the ecological and evolutionary equations and parameters.

The configuration object includes information on the model initialization, the observer function and the biological functions, namely speciation, dispersal, evolution, and ecology. The possibility of the user to customize these functions confers the flexibility and generality of *gen3sis*. The *initialization* function creates the ancestor species at the start of the simulation. Users can define the number of species, their distribution within the landscape and their trait values. The *observer* function allows recording any abiotic or biotic information of the virtual world as it changes by defining the outputs that are saved at specified time-steps. The configuration further allows to define the *speciation* function, the *dispersal* function, the *evolution* function and the *ecology* function, which are defined in the section on core processes below. Altogether, those six functions are specified in the configuration object and interact with a set of core objects defined in the simulation engine. The configuration object further lists the ecological traits considered in the simulation, sets a random seed allowing simulation reproducibility, start and end times of the simulation, as well as simulation aborting rules including the limits for a total or local number of species.

In our example, we start our simulation with the 40$^{th}$ time-step and end it a the latest available time-step of the landscape, i.e. 40 to 0Ma. We set ten different ancestor species, each colonizing one random cell in the landscape, with initial optimum temperature equal to the local temperature of the site. The dispersal of all species is the same and stochastic, following a Weibull distribution with shape 1.5 and scale 133. Speciation takes place after every second time-step (2 Myr) of isolation. Clustered populations (exchanging genes) had their trait homogenized and weighted by abundance, meaning that a temperature optimum of populations that are doing well in a site will contribute more to the average trait of a cluster. Temperature optimum evolves based on a normal distribution with standard deviation 0.001, possibly increasing or decreasing species optimum temperature. Abundances in a site are based on how close the species population optimum temperature 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. If abundances are below 1, species are considered extinct, if total abundance in a site (sum of all species 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. For more details regarding the configuration used here and how to modify it, please refer to [**create_config vignette**](create_config.html).


[^1]: O. Hagen, B. Flück, F. Fopp, J.S. Cabral, F. Hartig, M. Pontarp, T.F. Rangel, L. Pellissier. (2020). GENƎSIS: the GENeral Engine for Eco-Evolutionary SImulationS on the origins of biodiversity. (in review)


## 2. Run {#run}


Now that we have set up the general folder structure, the input data and configuration file, we can run a simulation using the *run_simulation* function. This function will:

1) Read in the config and prepare the output directories, initial landscape and call *create_ancestor_species* from the user config to create the initial species for the simulation.

2) Start the main loop over the time-steps. For every time-step it loads the appropriate landscape, removes all sites that became uninhabitable in the new time-step, and calls the main steps of any iteration.

3) At the end of every time-step, if desired, the simulation saves the species, landscapes, species richness patterns, etc... by calling the *end_of_timestep_observer* from the user config. In this function the user can implement customized observer functions, for example calculating summary statistics or creating species patterns plots.

4) When the simulation reaches the end, a phylogeny, runtime_information, the world starting conditions and the used config are saved at the output folder. If the output folder is not specified, it is deduced from the location of the landscape object. The function will return a summary object. In this vignette we will save it as *sim*.

Additionally to the main functions, the package provides several convenience functions to generate input data, configuration files, and plots. Moreover, all functions are accessible to the observer functions, which requests the variables for calculation during the model runs. This allows to transform the output of the simulation model into a format that can be readily analyzed. 


To launch a simulation you need to call the *run_simulation* function. We only store one intermediate step between the starting and end time-steps by setting "call_observer=1", this means one time-step equally distributed between the starting and ending time-steps. By doing so, we will have a look at our virtual world at 40, 20 and 0 Ma. 

```{r eval=FALSE, echo=TRUE}
#we set verbose to 0 to avoid a large console outputs of how the simulation is developing
sim <- run_simulation(config = file.path(datapath, "config/config_southamerica.R"), 
               landscape = file.path(datapath, "landscape"),
               output_directory = tempdir(),
               call_observer = 1,
               verbose=0)
```


```{r eval=T, echo=F}
sim <- readRDS(file.path(datapath, "output/config_southamerica/sgen3sis.rds"))
```

After the simulation is finished, a plot is created inside the SouthAmerica folder to show the state of the initial world. They are automatically stored in a sub folder called **output** inside our SouthAmerica folder and inside another folder named after our configuration file (i.e. config). It should look like this:

|>EXPERIMENT_FOLDER     |
|:-------------------------|
|&nbsp;&nbsp;&nbsp;>landscape|
|&nbsp;&nbsp;&nbsp;>config|
|&nbsp;&nbsp;&nbsp;**>output**|

## 3. Visualize {#visualize}

After having run the simulation, we can load and plot the species richness at different time-steps. Here, we will load them at the time-steps 40, 20 and 0 corresponding to 40, 20 and 0 Ma (1 time-step = 1 myr).

```{r eval=F, echo=T, fig.width=7, fig.height=2}
timesteps <- c(40, 20, 0)
oldpar <- par(no.readonly = TRUE)
par(mfrow=c(1,3))
for(i in timesteps){
  landscape_i <- readRDS(file.path(datapath, paste0('output/config_southamerica/landscapes/landscape_t_', i ,'.rds')))
  species_i <- readRDS(file.path(datapath, paste0('output/config_southamerica/species/species_t_', i ,'.rds')))
  plot_richness(species_i, landscape_i)
}
par(oldpar)
```


```{r eval=F, echo=F, fig.width=7, fig.height=2}
#GENERATE OFFLINE PLOT
timesteps <- c(40, 20, 0)
png("../inst/logo/richness_plot.png", height=500, width=1750, pointsize=28)
oldpar <- par(no.readonly = TRUE)
par(mfrow=c(1,3))
for(i in timesteps){
  landscape_i <- readRDS(file.path(datapath, paste0('output/config_southamerica/landscapes/landscape_t_', i ,'.rds')))
  species_i <- readRDS(file.path(datapath, paste0('output/config_southamerica/species/species_t_', i ,'.rds')))
  plot_richness(species_i, landscape_i)
}
par(oldpar)
dev.off()
```

```{r eval=T, echo=F, out.width='700px', height='200px'}
#PLOT OFFLINE PLOT
knitr::include_graphics("../inst/logo/richness_plot.png")
```


*gen3sis* also returns a summary object at the end of the simulation storing default summary statistics and important data over time. One of the default observer summaries is the species richness over time. 
**NOTE** You can also create your own observer function. 

We can quickly visualize some main results using the *plot_summary function*.
Here we stored the summary object (sgen3sis class) as *sim*.

```{r eval=TRUE, fig.width=7, fig.height=7}
plot_summary(sim)
```


There are a few visualization tools already included in the package, but you are free to explore and check the outputs with your favorite plotting functions and colors.

## 4. Analyze {#analyze}

Now that we have run the simulation, we are ready to perform some analysis to investigate the model behavior. Because we set a limit to the total abundances lower in colder and in more arid environments, we expect that the number of species that can co-exist also decreases in those environments.
We will perform a few statistical analyses to investigate those patterns based on the species richness pattern that emerged at the end of the simulation. To do so, we first load the landscape and combine it with the simulated species richness.

```{r eval=TRUE, fig.width=8, fig.height=5, fig.retina=2}
landscapes <- readRDS(file.path(datapath, "landscape", "landscapes.rds")) #get the input landscape
landscape_t0 <- as.data.frame(cbind(landscapes$temp[, 1:2], temp=landscapes$temp[, 3], arid=landscapes$arid[,3], area=landscapes$area[,3])) #get landscape at last time-step
landscape_t0 <- cbind(landscape_t0, rich=sim$summary$`richness-final`[,3]) #add richness to the dataframe
landscape_t0 <- na.omit(landscape_t0)
```



Much is possible here, from looking at phylogenies to spatial correlations. For this introduction we will keep things simple and investigate the relationship between richness and environmental variables at the final time-step of the simulation. For this, we will fit two univariate models, where each explains the relationship between species richness and one of the environmental variables.

For this, we first fit a generalized linear model between richness and temperature.
```{r eval=TRUE, echo=TRUE}
glm.uni <- glm(rich ~ poly(temp, 2), data=landscape_t0, family=poisson)
cor(landscape_t0$temp, landscape_t0$rich)
```

Second, we plot the response curve.
```{r eval=T, fig.width=7, fig.height=6, fig.retina=2, message=F, results='hide'}
# prepare data with temperature and predicted richness from our model
data_plot <- data.frame(cbind(landscape_t0$temp, predict(glm.uni,type = "response")))
# sort data for plotting and ommit NA's
data_plot <- na.omit(data_plot[order(data_plot[,1], decreasing = FALSE),])
# get the number of observations
n <- paste0('observations (n = ', length(landscape_t0$rich), ')')
# plot model curve
plot(data_plot[,1],data_plot[,2], xlab="Temperature [\u00B0C]", ylab=expression(paste(alpha," richness")), frame.plot=F, type="l", col='red', lwd=2, xlim=c(min(landscape_t0$temp), max(landscape_t0$temp)), ylim=c(min(landscape_t0$rich), max(landscape_t0$rich)))
# add observed points
points(landscape_t0$temp, landscape_t0$rich, col=rgb(0.5,0.5,0.5, alpha=0.4), pch=16) 
# add legend
legend(-20,30, col=c(rgb(0.5,0.5,0.5,0.4), 'red'), legend=c(n, 'model fit'), pch=c(16, NA), lty=c(NA, 1), lwd=c(NA, 2), bty='n')
```


By looking at the model fit of the generalized linear model we can see that the relationship between the two variables is rather hump-shaped, indicating that species richness is highest at intermediate temperature levels. First the higher richness found above 10° is consistent with the ecological rule that we imposed. At higher temperature, we have a higher energy and capacity facilitating coexist. This is consistent with our model expectations. Moreover, we can observe a large spread of the points at high temperature. This is likely because high temperature regions frequently have high aridity. Thus we will next investigate the association of aridity and species richness by fitting a generalized linear model between richness and aridity.


```{r eval=TRUE, fig.width=7, fig.height=6, fig.retina=2, results='hide', echo=F}
glm.uni <- glm(rich ~ poly(arid, 2), data=landscape_t0, family=poisson)
cor(landscape_t0$arid, landscape_t0$rich)

#Plot the response curve
data_plot <- data.frame(cbind(landscape_t0$arid, predict(glm.uni,type = "response")))
sorted <- na.omit(data_plot[order(data_plot[,1], decreasing = FALSE),])
data_plot <- data.frame(cbind(sorted[,1],sorted[,2]))
plot(data_plot[,1],data_plot[,2], xlab="Aridity", ylab=expression(paste(alpha," richness")), frame.plot=F, type="l", col='blue', 
     lwd=2, xlim=c(min(landscape_t0$arid), max(landscape_t0$arid)), ylim=c(min(landscape_t0$rich), max(landscape_t0$rich)))
points(landscape_t0$arid, landscape_t0$rich, col=rgb(0.5,0.5,0.5, alpha=0.4), pch=16) 
legend(0.4,30, col=c(rgb(0.5,0.5,0.5,0.4), 'blue'), legend=c(n, 'model fit'), pch=c(16, NA), lty=c(NA, 1), lwd=c(NA, 2), bty='n')
```


The number of species is negatively related to aridity. The explained deviance of this relationship is D2=`r round((glm.uni$null.deviance - glm.uni$deviance) / glm.uni$null.deviance,3)`. This relationship emerge because we set a lower carrying capacity in sites with higher level of aridity. Hence, the models behaves as expected from the configuration function.


Finally, we fit a multivariate model for explaining species richness by taking into account the two predictors together.

```{r eval=TRUE}
richness_model <- glm(rich ~ poly(temp, 2) + poly(arid, 2), family='quasipoisson', data=landscape_t0)
summary(richness_model)$coefficients
```

Above, we report the explained deviance of the glm model on richness considering temperature and aridity. The explained deviance of this relationship is D2=`r round((richness_model$null.deviance - richness_model$deviance) / richness_model$null.deviance,3)`. This explained deviance will not be total since there are other effects including dispersal limitations, that are not considered in this statistical summary. It is important to consider that we have now only been analyzing the species richness at one time-step, namely the last one (0 Ma). The emerged patterns are more likely not only to be explained by the environmental conditions of that time but also dependent on historical conditions and multiple processes interactions tackled by *gen3sis*. 

For more details on how to create a landscape and prepare it to use in *gen3sis* see [**design_landscape vignette**](design_landscape.html) and [**create_input_landscape vignette**](create_input_landscape.html). For more details on how to create or modify a configuration object see [**create_config vignette**](create_config.html). It's now up to you to explore the "virtual" world.
