03_Module_1.Rmd

---
title: 'D-place FARM documentation: Module 1'
author: "Ty Tuff, Bruno Vilela, and Carlos Botero"
date: 'project began: 15 May 2016, document updated: `r strftime(Sys.time(), format
  = "%d %B %Y")`'
output:
  html_notebook: default
  html_document: default
  pdf_document: default
  word_document: default
bibliography: FARM package.bib
---

# Module 1: Simulation to produce a world and a tree 

  After the spatial network was built to host the simulations, the machinery for tracking simulation progress was in place, and the input parameters were assigned, then each simulation could launch according to a set of initial launch rules and then play out according to a separate set of continuing behavioral rules. The initial launch rules begin by stipulating that the first human culture will originate within modern-day Ethiopia and radiate outward from that point. All of these original cultures use foraging as their only mode of subsistence until at least 80% of the available nodes are occupied by a culture and then agriculture can start to arise by foragers switching to agriculturalists. The ability to switch from forager to agriculturalist without any outside influence is restricted to within known origins of agriculture as defined by Larson et al. (2014).  Switching from forager to agriculturalist can only occur on nodes within those known origins, but agriculturalists can switch to foragers with equal probability on any node. Nodes can only interact via edges so, edges define neighboring nodes. These initial launch rules are the same for all simulations, regardless of the hypothesized mechanism they are simulating. 

  When simulations were terminated at 30,000 time steps, they will have produced a spatial pattern of subsistence modes distributed across the geographic network and a phylogeny describing the relatedness of those societies through time. These two objects contain a great deal of information about the trajectory of particular replicates but that complexity makes it difficult to compare replicates directly. No test currently exists to compare a combination of spatial and phylogenetic patterns between replicated simulations. To simplify these outputs and allow them to be compared directly to each other, we calculated a suite of 12 summary statistics to describe different features of both the spatial distribution and the phylogeny. During preliminary tests of the simulations, we used 19 summary statistics but later eliminated 7 of those original statistics because they didn't stabilize over time to where they could be described by a linear model. The 12 remaining statistics stabilized by either increasing regularly through time or came to an asymptote as the simulation reached equilibrium (SI figure…).  After all of the simulations were complete, these summary statistics were concatenated into a single dataset, labelled with the rule set (mechanism) used to create each one, and then passed on to the random forest algorithm for analysis.


```{r eval=FALSE}
# Install the most recent version of FARM from a .zip file
install.packages(file.choose(), repos=NULL) 
```

```{r}
library(FARM)
ls("package:FARM")
```

## Inputs
 The specific behavior of any particular replicate simulation is controlled by 17 input parameters. All of these values are unique to each replicate simulation, but don't change from the beginning to the end of a simulation.  These parameters are all constrained between 0 and 1 and fall within five general categories: speciation, extinction, cultural diffusion, diffusion by takeover, and arisal. The first two categories, speciation and extinction, each require four parameters to describe the different ways that farmers and foragers can interact with environments that are suitable or unsuitable for agriculture. Diffusion by takeover also requires four parameters to describe the four ways that farmers or foragers can displace their neighbors to expand their subsistence mode. Cultural diffusion is simpler than diffusion by takeover and only requires two input parameters: diffusion from farmers to foragers and diffusion from foragers to farmers. Arisal follows the same convention as the cultural diffusion inputs but the probability of switching from foragers to farmers is constraining to within the known origins of agriculture (see Figure 2). 

Most of the input parameters are drawn randomly from a uniform distribution constrained between 0 and 1, but some of these uniform random draws are constrained more narrowly. Parameter values in the speciation or extinction categories are ordered so that farmers in environments suitable for farming will have the highest speciation rate and lowest extinction rate, farmers in environments unsuitable for agriculture will have the lowest speciation rate and highest extinction rate, and all foragers will have an intermediate probability between them. The first value selected during a random draw is assigned to the high probability input, the second is constrained to be smaller than the first draw and assigned to the low probability input, the third is constrained between the first two values. Arisal is constrained so that the probability of switching from forager to farmer is always higher than switching from farmer to forager.  

The input parameters are also where the model type is assigned to each simulation. The hypothesized mechanism prescribed to each simulation are relayed to the algorithm by setting specific sets of parameters to 0 so certain events have no chance of happening. The full ensemble model (+culture +takeover) assigns a value to all 17 parameters, the cultural diffusion only model (+culture) sets all takeover values to 0, the diffusion by takeover model (+takeover) sets all cultural diffusion values to 0, and the basic model sets all culture and takeover values to 0.  


## Module 1 functions
#### The first set of RunSim functions are the default pipeline where only one output is saved at the end of the simulation. 


This first function controls error messages coming from the primary function below. 
```{r eval=FALSE}
# Run the simulation function skiping the erros and atributing NA if it occurs
RunSimUltimate <- function(myWorld, P.extinction, P.speciation,
                           P.diffusion, P.Arisal, P.TakeOver, nbs, independent,
                           N.steps, multiplier,
                           silent = TRUE, start = NULL) {

  result <- try(RunSim(myWorld, P.extinction, P.speciation,
                       P.diffusion, P.Arisal, P.TakeOver, nbs,
                       independent, N.steps,
                       multiplier, start = start), silent = silent)
  if (class(result) == "try-error") {
    result <- NA
  }
  return(result)
}

```

This is the primary function running the simulation.
```{r eval=FALSE}
#==================================================================
# SimulationFunctions.R
#
# Contains a function for simulation of cultural evolution in space and time
# Allows for (1) Vertical Transmission (phylogenetic inheritance); (2) Horizontal
# Transmission (cultural diffusion); (3) Ecological selection (Both speciation and
# extinction are determined by the match between the state of a binary trait and the
# environment a societuy occupies).
#
# 7 Jun 2016
# Carlos A. Botero, Bruno Vilela & Ty Tuff
# Washington University in Saint Louis
#==================================================================
RunSim <- function(myWorld, P.extinction, P.speciation,
                   P.diffusion, P.Arisal, P.TakeOver, nbs, independent,
                   N.steps, multiplier, start) {
  # myWorld = The hexagonal world created with the function BuildWorld
  # P.extinction = Probability matrix of extinction
  # P.speciation = Probability matrix of speciation
  # P.diffusion = Probability matrix of diffusion
  # P.Arisal = Probability matrix of arisal
  # P.TakeOver = Probability matrix of takeover
  # N.steps = Number of steps in the model
  # multiplier = The number that will multiply the probabilities according
  # to environmetal fitness.
  # start = the point ID in 'myWorld' that will give risen to humans.
  # (humans origin will be in one of the existing positions)

  world.size <- nrow(myWorld)
  # Initialize parameters we will use later to build the phylogeny
  rootnode <-  world.size + 1 # standard convention for root node number

  # set the seed for simulation
  if (is.null(start)) {
  start <- sample(1:world.size, 1)
  }

  myWorld[start, 4:6] <- c(0, 0, 1) # Setting root(0), time(0), ancestral(1, forager)

  mytree <- TheOriginOfSpecies(world.size, start) # Empty tree
  myT <- 0 # Time starts at zero

  # Common input and output for all the internal modules
  input <- list(P.speciation, P.Arisal, P.diffusion, P.extinction, P.TakeOver,
                myWorld, mytree, myT, multiplier, nbs, independent)

  # Functions order to be randomized
  rand_order_func_run <- list("Extinction", "Diffusion", "SpeciationTakeOver", "Arisal")

  cat("0% [") # Time count

  for (steps in 1:N.steps) { # Starts the loop with 'n' steps

    if (steps %% round((N.steps / 10)) == 0) { # Time count
      cat('-') # Time count
    }# Time count
    if (steps == N.steps) { # Time count
      cat("] 100 %\n")# Time count
    }# Time count

    # Randomize functions order
    rand_order <- sample(rand_order_func_run)
    # Run the functions
    input <- do.call(rand_order[[1]], list(input = input))
    input <- do.call(rand_order[[2]], list(input = input))
    input <- do.call(rand_order[[3]], list(input = input))
    input <- do.call(rand_order[[4]], list(input = input))

  }
  # Trunsform the input/output into the final result and return it
  myWorld <- as.data.frame(input[[6]])
  myWorld[, 8] <- paste0("t", myWorld[, 8])
  mytree <- makePhy(input[[7]])
  mytree$edge.length <- mytree$edge.length / N.steps
  return(list('mytree' = mytree, 'myWorld' = myWorld))
}

```


We track each simulation from start to finish using two data storage objects, a spatial object storing the current subsistence mode of each node and a phylogenetic object tracking the relationship between all nodes through time. For each action taken during the simulation, the algorithm first checks the current state of these two objects, then prescribes a change to a single node within the network according to a set of rules, and then modifies that node in both storage objects before moving on to the next node. This process is repeated for each node within each time step and randomized across time steps to create an entire simulation. The phylogenetic object is built under standard evolutionary assumptions necessary for many of the summary statistics used later, so the tree is always bifurcating, non-reticulated, and ultrametric. 



## Push versions 

```{r eval=FALSE}
# Run the simulation function skiping the erros and atributing NA if it occurs
RunSimUltimate.push <- function(myWorld, P.extinction, P.speciation,
                                P.diffusion, P.Arisal, P.TakeOver, nbs, independent,
                                N.steps, multiplier,
                                silent = TRUE, start = NULL) {

  result <- try(RunSim.push(myWorld, P.extinction, P.speciation,
                            P.diffusion, P.Arisal, P.TakeOver, nbs,
                            independent, N.steps,
                            multiplier, start = start), silent = silent)
  if (class(result) == "try-error") {
    result <- NA
  }
  return(result)
}


```





```{r eval=FALSE}
#==================================================================
# SimulationFunctions.R
#
# Contains a function for simulation of cultural evolution in space and time
# Allows for (1) Vertical Transmission (phylogenetic inheritance); (2) Horizontal
# Transmission (cultural diffusion); (3) Ecological selection (Both speciation and
# extinction are determined by the match between the state of a binary trait and the
# environment a societuy occupies).
#
# 7 Jun 2016
# Carlos A. Botero, Bruno Vilela & Ty Tuff
# Washington University in Saint Louis
#==================================================================
RunSim.push <- function(myWorld, P.extinction, P.speciation,
                   P.diffusion, P.Arisal, P.TakeOver, nbs, independent,
                   N.steps, multiplier, start) {
  # myWorld = The hexagonal world created with the function BuildWorld
  # P.extinction = Probability matrix of extinction
  # P.speciation = Probability matrix of speciation
  # P.diffusion = Probability matrix of diffusion
  # P.Arisal = Probability matrix of arisal
  # P.TakeOver = Probability matrix of takeover
  # N.steps = Number of steps in the model
  # multiplier = The number that will multiply the probabilities according
  # to environmetal fitness.
  # start = the point ID in 'myWorld' that will give risen to humans.
  # (humans origin will be in one of the existing positions)

  world.size <- nrow(myWorld)
  # Initialize parameters we will use later to build the phylogeny
  rootnode <-  world.size + 1 # standard convention for root node number

  # set the seed for simulation
  if (is.null(start)) {
    start <- sample(1:world.size, 1)
  }

  myWorld[start, 4:6] <- c(0, 0, 1) # Setting root(0), time(0), ancestral(1, forager)

  mytree <- TheOriginOfSpecies(world.size, start) # Empty tree
  myT <- 0 # Time starts at zero

  # Common input and output for all the internal modules
  input <- list(P.speciation, P.Arisal, P.diffusion, P.extinction, P.TakeOver,
                myWorld, mytree, myT, multiplier, nbs, independent)

  # Functions order to be randomized
  rand_order_func_run <- list("Extinction", "Diffusion",
                              "SpeciationTakeOver.push", "Arisal")

  cat("0% [") # Time count

  for (steps in 1:N.steps) { # Starts the loop with 'n' steps

    if (steps %% round((N.steps / 10)) == 0) { # Time count
      cat('-') # Time count
    }# Time count
    if (steps == N.steps) { # Time count
      cat("] 100 %\n")# Time count
    }# Time count

    # Randomize functions order
    rand_order <- sample(rand_order_func_run)
    # Run the functions
    input <- do.call(rand_order[[1]], list(input = input))
    input <- do.call(rand_order[[2]], list(input = input))
    input <- do.call(rand_order[[3]], list(input = input))
    input <- do.call(rand_order[[4]], list(input = input))

  }
  # Trunsform the input/output into the final result and return it
  myWorld <- as.data.frame(input[[6]])
  myWorld[, 8] <- paste0("t", myWorld[, 8])
  mytree <- makePhy(input[[7]])
  mytree$edge.length <- mytree$edge.length / N.steps
  return(list('mytree' = mytree, 'myWorld' = myWorld))
}

```







































#### The second set of RunSim functions save an output each timestep if we want to look at trends through time. We use this to make videos of the simulation running. 
```{r}
RunSimUltimate2 <- function (myWorld, P.extinction, P.speciation, P.diffusion, P.Arisal, 
    P.TakeOver, nbs, independent, N.steps, multiplier, silent = TRUE, 
    count, resolution = seq(1, N.steps, 100), P.Arisal0, start = NULL) 
{
    result <- try(RunSim2(myWorld, P.extinction, P.speciation, 
        P.diffusion, P.Arisal, P.TakeOver, nbs, independent, 
        N.steps, multiplier, count = count, resolution = resolution, 
        P.Arisal0 = P.Arisal0, start), silent = silent)
    if (class(result) == "try-error") {
        result <- NA
    }
    return(result)
}
```


```{r}
RunSim2 <- function (myWorld, P.extinction, P.speciation, P.diffusion, P.Arisal, 
    P.TakeOver, nbs, independent, N.steps, multiplier, count, 
    resolution, P.Arisal0, start = NULL) 
{
    folder <- paste0("./Module_1_outputs/myOut_rep_", formatC(count, 
        width = 2, flag = 0), "_combo_", formatC(count, width = 2, 
        flag = 0), "_", "params", "_P.speciation_", paste(formatC(P.speciation, 
        width = 2, flag = 0), collapse = "_"), "_P.extinction_", 
        paste(formatC(P.extinction, width = 2, flag = 0), collapse = "_"), 
        "_P.diffusion_", paste(formatC(P.diffusion, width = 2, 
            flag = 0), collapse = "_"), "_P.TO_", paste(formatC(P.TakeOver, 
            width = 2, flag = 0), collapse = "_"), "_P.Arisal_", 
        paste(formatC(P.Arisal0, width = 2, flag = 0), collapse = "_"), 
        "_timesteps_", N.steps)
    world.size <- nrow(myWorld)
    rootnode <- world.size + 1
    if (is.null(start)) {
        start <- sample(1:world.size, 1)
    }
    myWorld[start, 4:6] <- c(0, 0, 1)
    mytree <- TheOriginOfSpecies(world.size, start)
    myT <- 0
    input <- list(P.speciation, P.Arisal, P.diffusion, P.extinction, 
        P.TakeOver, myWorld, mytree, myT, multiplier, nbs, independent)
    rand_order_func_run <- list("Extinction", "Diffusion", "SpeciationTakeOver", 
        "Arisal")
    cat("0% [")
    for (steps in 1:N.steps) {
        if (steps%%round((N.steps/10)) == 0) {
            cat("-")
        }
        if (steps == N.steps) {
            cat("] 100 %\n")
        }
        rand_order <- sample(rand_order_func_run)
        input <- do.call(rand_order[[1]], list(input = input))
        input <- do.call(rand_order[[2]], list(input = input))
        input <- do.call(rand_order[[3]], list(input = input))
        input <- do.call(rand_order[[4]], list(input = input))
        if (steps %in% resolution) {
            myWorld <- as.data.frame(input[[6]])
            myWorld[, 8] <- paste0("t", myWorld[, 8])
            if (nrow(na.omit(input[[7]])) > 1) {
                mytree <- makePhy(input[[7]])
            }
            else {
                mytree <- NA
            }
            myOut <- list(mytree = mytree, myWorld = myWorld)
            save(myOut, file = paste0(folder, "_", formatC(steps, 
                10, flag = 0), ".Rdata"))
            stats <- Module_2(myOut)
            save(stats, file = paste0(folder, "_", formatC(steps, 
                10, flag = 0), "_stats", ".Rdata"))
        }
    }
    myWorld <- as.data.frame(input[[6]])
    myWorld[, 8] <- paste0("t", myWorld[, 8])
    mytree <- makePhy(input[[7]])
    mytree$edge.length <- mytree$edge.length/N.steps
    return(list(mytree = mytree, myWorld = myWorld))
}
```


#..And the push version of saving each time step


```{r}
# Run the simulation function skiping the erros and atributing NA if it occurs
RunSimUltimate2.push <- function(myWorld, P.extinction, P.speciation,
                            P.diffusion, P.Arisal, P.TakeOver, nbs, independent,
                            N.steps, multiplier,
                            silent = TRUE, count, resolution = seq(1, N.steps, 100),
                            P.Arisal0, start = NULL) {

  result <- try(RunSim2.push(myWorld, P.extinction, P.speciation,
                        P.diffusion, P.Arisal, P.TakeOver, nbs,
                        independent, N.steps,
                        multiplier, count = count, resolution = resolution,
                        P.Arisal0 = P.Arisal0, start),
                silent = silent)
  if (class(result) == "try-error") {
    result <- NA
  }
  return(result)
}

```


```{r}
#==================================================================
RunSim2.push <- function(myWorld, P.extinction, P.speciation,
                    P.diffusion, P.Arisal, P.TakeOver, nbs, independent,
                    N.steps, multiplier, count, resolution, P.Arisal0,
                    start = NULL) {
  # myWorld = The hexagonal world created with the function BuildWorld
  # P.extinction = Probability matrix of extinction
  # P.speciation = Probability matrix of speciation
  # P.diffusion = Probability matrix of diffusion
  # P.Arisal = Probability matrix of arisal
  # P.TakeOver = Probability matrix of takeover
  # N.steps = Number of steps in the model
  # multiplier = The number that will multiply the probabilities according
  # to environmetal fitness.
  # start = the point ID in 'myWorld' that will give risen to humans.
  # (humans origin will be in one of the existing positions)
  folder <- paste0("./Module_1_outputs/myOut_rep_",
                   formatC(count, width = 2,flag = 0),
                   "_combo_",
                   formatC(count, width = 2,flag = 0),
                   "_","params", "_P.speciation_",
                   paste(formatC(P.speciation, width = 2,flag = 0),
                         collapse="_"),"_P.extinction_",
                   paste(formatC(P.extinction, width = 2,flag = 0),
                         collapse="_"), "_P.diffusion_",
                   paste(formatC(P.diffusion, width = 2,flag = 0),
                         collapse="_"), "_P.TO_",
                   paste(formatC(P.TakeOver, width = 2,flag = 0),
                         collapse="_"),"_P.Arisal_",
                   paste(formatC(P.Arisal0, width = 2,flag = 0),
                         collapse="_"), "_timesteps_",
                   N.steps)
  world.size <- nrow(myWorld)
  # Initialize parameters we will use later to build the phylogeny
  rootnode <-  world.size + 1 # standard convention for root node number

  # set the seed for simulation
  if (is.null(start)) {
    start <- sample(1:world.size, 1)
  }

  myWorld[start, 4:6] <- c(0, 0, 1) # Setting root(0), time(0), ancestral(1, forager)

  mytree <- TheOriginOfSpecies(world.size, start) # Empty tree
  myT <- 0 # Time starts at zero

  # Common input and output for all the internal modules
  input <- list(P.speciation, P.Arisal, P.diffusion, P.extinction, P.TakeOver,
                myWorld, mytree, myT, multiplier, nbs, independent)

  # Functions order to be randomized
  rand_order_func_run <- list("Extinction", "Diffusion",
                              "SpeciationTakeOver.push", "Arisal")

  cat("0% [") # Time count

  for (steps in 1:N.steps) { # Starts the loop with 'n' steps

    if (steps %% round((N.steps / 10)) == 0) { # Time count
      cat('-') # Time count
    }# Time count
    if (steps == N.steps) { # Time count
      cat("] 100 %\n")# Time count
    }# Time count

    # Randomize functions order
    rand_order <- sample(rand_order_func_run)
    # Run the functions
    input <- do.call(rand_order[[1]], list(input = input))
    input <- do.call(rand_order[[2]], list(input = input))
    input <- do.call(rand_order[[3]], list(input = input))
    input <- do.call(rand_order[[4]], list(input = input))
    # Save
    if(steps %in% resolution) {
      myWorld <- as.data.frame(input[[6]])
      myWorld[, 8] <- paste0("t", myWorld[, 8])
      if(nrow(na.omit(input[[7]])) > 1) {
        mytree <- makePhy(input[[7]])
      } else {
        mytree <- NA
      }
      myOut <- list('mytree' = mytree, 'myWorld' = myWorld)
      save(myOut, file= paste0(folder,"_", formatC(steps, 10, flag = 0), ".Rdata"))
      stats <- Module_2(myOut)
      save(stats, file= paste0(folder,"_", formatC(steps, 10, flag = 0),
                               "_stats", ".Rdata"))
    }
  }
  # Trunsform the input/output into the final result and return it
  myWorld <- as.data.frame(input[[6]])
  myWorld[, 8] <- paste0("t", myWorld[, 8])
  mytree <- makePhy(input[[7]])
  mytree$edge.length <- mytree$edge.length / N.steps
  return(list('mytree' = mytree, 'myWorld' = myWorld))
}


```