Midpoint_Report_Code.Rmd

This file represents the coded portion of our project, which will be combined with our written portion for the final submission


Just loading any necessary libraries
```{r}
library(tidyverse)

```

For proof-of-concept we are currently manually creating one specific world and set of parameters to demonstrate how the functions will run on this example
```{r}
world <- list("1" = list("n1" = 10, "n2" = 10), "2" = list("n1" = 0, "n2" = 0), "3" =  list("n1" = 0, "n2" = 0))
params <- list("n1" = list("r" = 0.1, "k" = 1000), "n2" = list("r" = 0.9, "k" = 100))
```
Now that the world is created, we have made a set of functions and algorithms to run our simulation.


The following will determine how many individuals will migrate from a patch based on a binomial draw
```{r}
migrants_number <- function(n){
  return(rbinom(n = 1, size = n, prob = 0.5))
}
```

The following will randomly determine which direction a population will migrate to.
```{r}
direction_generator <-function(){
  return(sample(c(-1, 1), 1))
}
direction_generator()

```

The following is the population growth, that calculates by how much a given population will grow based on the logistic growth equation 
```{r}
growth <- function(r, k, n){
  return((r*n*(1-(n/k))))
}
growth(r = 0.5, k = 100, n = 10)
```

The following chunk simulates the the growth and migration for 10 timesteps, note that there are a lot of things that we have hard coded a lot of the steps here instead of creating functions for everything so that we have a proof of concept that works. It first applies the logistic growth rate on all the populations in all the patches, then applies the movement function on each population in each patch after they have grown for 1 timestep, as this meets our assumptions more accurately.
the plan is to have each of these steps converted to their unique functions of growth and movement across any world of any size. Then they will be combined to one "simulate" function so that eventually we can include the checker that will stop the simulation when one of the populations wins. 

please ignore the comments which were used for troubleshooting

also note that currently the _ceiling()_ function is being called on the n parameter when determining the number of migrants, as not doing so causes rbinom to take a float, that was generated as a result of us using the logistic growth function, and rbinom cannot work on floats, it requires whole numbers/integers
```{r}
i <- 0
while(i<10){
  for(patch in seq(1, length(world))){
    #print(seq(1, length(world)))
    #print(patch)
    for(pop in seq(1, length(world[[patch]]))){
      #print(seq(1, length(patch)))
      #print(pop)
      curr_pop_name <- names(world[[patch]][pop])
      #print(curr_pop_name)
      curr_pop <- world[[patch]][[curr_pop_name]]
      #print(curr_pop)
      #print(growth(r = params[[curr_pop_name]][["r"]], k = params[[curr_pop_name]][["k"]], n = curr_pop))
      world[[patch]][[curr_pop_name]] <- curr_pop + growth(r = params[[curr_pop_name]][["r"]], k = params[[curr_pop_name]][["k"]], n = curr_pop) -> curr_pop
      #print(curr_pop)
    }}
  for(patch in seq(1, length(world))){
    #print(seq(1, length(world)))
    #print(patch)
    for(pop in seq(1, length(world[[patch]]))){
      #print(seq(1, length(patch)))
      #print(pop)
      curr_pop_name <- names(world[[patch]][pop])
     #print(curr_pop_name)
      curr_pop <- world[[patch]][[curr_pop_name]]
      #print(curr_pop)
      migrants <- migrants_number(ceiling(curr_pop))
      direction <- direction_generator()
      destination <- patch + direction
      #print(destination)
      if(all(destination > 0 & destination < length(world))){
        world[[patch]][[curr_pop_name]] <- world[[patch]][[curr_pop_name]] - migrants
        world[[destination]][[curr_pop_name]] <- world[[destination]][[curr_pop_name]] + migrants
    }
   }
  }
  i <- i + 1
}    

```

From here, we plan on tracking the winners, and seeing frequently each population wins as their parameters(r and k) change, and plot our results

Potential plottign code taken from 'Improving the Presentation of Results of Logistic Regression with R' by Marcelino de la Cruz Rot code will be modified to better suit our purpose.

```{r]
datafake <- data.frame('n1_wins' = c(35, 45, 50), 'n2_wins'= c(65, 55, 50), 
                       'r1' = c(0.45, 0.5, 0.55), 'r2' = c(0.9, 0.85, 0.8), 
                       'k1' = c(1000, 900, 800),  'k2' = c(500, 550, 600))


plot.logi.hist <- function(independ, depend, logi.mod = 1, type = "dit", boxp = TRUE, rug = FALSE, las.h = 1, ...) {
  
  # Get the label for the x-axis
  xlabel <- paste(deparse(substitute(independ)))
  
  # Define the functions:
  
  # Scatter plot function
  logi.scater <- function(independ, depend, scater = "n", x.lab = xlabel, las = las.h) {
    plot(independ, depend, cex = 1, type = scater, ylab = "Predicted probability", xlab = x.lab, cex.lab = 1.5, las = las)
  }
  
  # Add rug plot if desired
  logi.rug <- function(independ, depend, pch.rug = 16, cex.rug = 1) {
    points(independ, depend, pch = pch.rug, cex = cex.rug)
  }
  
  # Box plot function
  logi.box <- function(independ, depend, col.box = "gray", x.lab = xlabel, las = las.h) {
    plot(independ, depend, cex = 1, type = "n", ylim = c(-0.1, 1.1), ylab = "Predicted probability", xlab = x.lab, cex.lab = 1.5, las = las)
    indep.1 <- independ[depend == 1]
    indep.0 <- independ[depend == 0]
    boxplot(indep.1, horizontal = TRUE, add = TRUE, at = 1.05, boxwex = 0.1, col = col.box, notch = TRUE)
    boxplot(indep.0, horizontal = TRUE, add = TRUE, at = -0.05, boxwex = 0.1, col = col.box, notch = TRUE)
  }
  
  # Fit binomial GLM and add predicted curve
  logi.curve <- function(independ, depend, mod = logi.mod, col.cur = "red", lwd.cur = 4) {
    if (mod == 1) mod3 <- glm(depend ~ independ, family = binomial)
    if (mod == 2) mod3 <- glm(depend ~ independ + I(independ^2), family = binomial)
    x.new <- seq(min(independ), max(independ), len = 100)
    y.new <- predict(mod3, data.frame(independ = x.new), type = "response")
    lines(x.new, y.new, lwd = lwd.cur, col = col.cur)
  }
  
  # Add dit plot
  logi.dit <- function(independ, depend, cex.p = 1, pch.dit = 1, incre = 0.02) {
    indep.0 <- independ[depend == 0]
    indep.1 <- independ[depend == 1]
    uni.plot.0 <- function(x) length(which(indep.0 == x))
    uni.plot.1 <- function(x) length(which(indep.1 == x))
    # Get the number of repeated values of “independ”:
    cosa.0 <- apply(as.matrix(unique(indep.0)), 1, uni.plot.0)
    cosa.1 <- apply(as.matrix(unique(indep.1)), 1, uni.plot.1)
    # Start plotting:
    points(independ, depend, pch = pch.dit, cex = cex.p)
    for (i in 1:max(cosa.0)) {
      for (j in 1:i) {
        points(unique(indep.0)[which(cosa.0 == i + 1)], rep(0 + incre * j, length(which(cosa.0 == i + 1))), pch = pch.dit, cex = cex.p)
      }
    }
    for (i in 1:max(cosa.1)) {
      for (j in 1:i) {
        points(unique(indep.1)[which(cosa.1 == i + 1)], rep(1 - incre * j, length(which(cosa.1 == i + 1))), pch = pch.dit, cex = cex.p)
      }
    }
  }
  
  # Add histograms and frequency axes
  logi.hist <- function(independ, depend, scale.hist = 5, col.hist = gray(0.7), count.hist = FALSE, intervalo = 0, las.h1 = las.h) {
    h.br <- hist(independ, plot = F)$br
    if (intervalo > 0) h.br <- seq(from = range(h.br)[1], to = range(h.br)[2], by = intervalo)
    h.x <- hist(independ[depend == 0], breaks = h.br, plot = F)$mid
    h.y0 <- hist(independ[depend == 0], breaks = h.br, plot = F)$counts
    h.y1 <- hist(independ[depend == 1], breaks = h.br, plot = F)$counts
    h.y0n <- h.y0 / (max(c(h.y0, h.y1)) * scale.hist)
    h.y1n <- 1 - h.y1 / (max(c(h.y0, h.y1)) * scale.hist)
    # Draw bottom histogram:
    for (i in 1:length(h.y0n)) {
      if (h.y0n[i] > 0) polygon(c(rep(h.br[i], 2), rep(h.br[i + 1], 2)), c(0, rep(h.y0n[i], 2), 0), col = col.hist)
    }
    # Draw top histogram:
    for (i in 1:length(h.y1n)) {
      if (h.y1n[i] < 1) polygon(c(rep(h.br[i], 2), rep(h.br[i + 1], 2)), c(h.y1n[i], 1, 1, h.y1n[i]), col = col.hist)
    }
    # Add counts to bins if required:
    if (count.hist == TRUE) {
      for (i in 1:length(h.x)) {
        text(h.x[i], h.y1n[i], h.y1[i], cex = 1, pos = 1)
        text(h.x[i], h.y0n[i], h.y0[i], cex = 1, pos = 3)
      }
    }
    # Plot the axes of histograms
    axis.hist <- function(h.y0, h.y1, scale.hist, las = las.h1) {
      tope <- max(c(h.y0, h.y1))
      label.down <- c(0, (ceiling(tope / 10)) * 5, (ceiling(tope / 10)) * 10)
      label.up <- c((ceiling(tope / 10)) * 10, (ceiling(tope / 10)) * 5, 0)
      at.down <- label.down / (tope * scale.hist)
      at.up <- 1 - (label.up / (tope * scale.hist))
      at.hist <- c(at.down, at.up)
      label.hist <- c(label.down, label.up)
      axis(side = 4, at = at.hist, labels = label.hist, las = las)
      mtext("Frequency", side = 4, line = 2, cex = 1.5)
    }
    axis.hist(h.y0, h.y1, scale.hist)
    axis(side = 2, las = las.h1)
  }
  
  # Set the margins of plot area
  old.mar <- par()$mar
  par(mar = c(5.1, 4.1, 4.1, 4.1))
  
  # Plot the combined graph
  if (boxp == TRUE) logi.box(independ, depend)
  if (boxp == FALSE) logi.scater(independ, depend)
  if (type != "dit") logi.hist(independ, depend, ...)
  if (rug == TRUE) logi.rug(independ, depend)
  logi.curve(independ, depend)
  if (type == "dit") logi.dit(independ, depend)
  
  # Reset the margins to old margins
  par(mar = old.mar)
}


depend <- ifelse(datafake$n1_wins > median(datafake$n1_wins), 1, 0)
datafake$depend <- depend

plot.logi.hist(independ = datafake$r1, depend = datafake$depend,   
  logi.mod = 1,  type = "dit",  boxp = TRUE, rug = TRUE, las.h = 1                   
)

```