This repository has been archived by the owner on Dec 22, 2024. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathProject_Script.Rmd
156 lines (129 loc) · 6.63 KB
/
Project_Script.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
Just loading any necessary libraries
```{r libraries}
library(tidyverse)
library(ggplot2)
```
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 worlds and params generator}
my_world <- list("1" = list("n1" = 10, "n2" = 10),
"2" = list("n1" = 0, "n2" = 0),
"3" = list("n1" = 0, "n2" = 0),
"4" = list("n1" = 0, "n2" = 0),
"5" = list("n1" = 0, "n2" = 0),
"6" = list("n1" = 0, "n2" = 0),
"7" = list("n1" = 0, "n2" = 0),
"8" = list("n1" = 0, "n2" = 0),
"9" = list("n1" = 0, "n2" = 0),
"10" = list("n1" = 0, "n2" = 0),
"11" = list("n1" = 0, "n2" = 0),
"12" = list("n1" = 0, "n2" = 0),
"13" = list("n1" = 0, "n2" = 0),
"14" = list("n1" = 0, "n2" = 0),
"15" = list("n1" = 0, "n2" = 0),
"16" = list("n1" = 0, "n2" = 0),
"17" = list("n1" = 0, "n2" = 0),
"18" = list("n1" = 0, "n2" = 0),
"19" = list("n1" = 0, "n2" = 0),
"20" = list("n1" = 0, "n2" = 0))
my_params <- list("n1" = list("r" = 0.9, "k" = 900),
"n2" = list("r" = 0.1, "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 generator}
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}
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}
growth <- function(r, k, n){
return((r*(1-min(n,k)/k))*min(n,k))
}
growth(r = 10, 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 THE SIMULATOR}
simulate <- function(world, params){
i <- 0
while(i<1000){
#growth on all patches
for(patch in seq(1, length(world))){
for(pop in seq(1, length(world[[patch]]))){
curr_pop_name <- names(world[[patch]][pop])
curr_pop <- world[[patch]][[curr_pop_name]]
world[[patch]][[curr_pop_name]] <- world[[patch]][[curr_pop_name]] + round(growth(r = params[[curr_pop_name]][["r"]], k = params[[curr_pop_name]][["k"]], n = curr_pop)) -> curr_pop
}
}
#movement on all patches
for(patch in seq(1, length(world))){
for(pop in seq(1, length(world[[patch]]))){
curr_pop_name <- names(world[[patch]][pop])
curr_pop <- world[[patch]][[curr_pop_name]]
migrants <- migrants_number(floor(curr_pop))
direction <- direction_generator()
destination <- patch + direction
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
}
}
}
#checking if one of them won
if(sum(unlist(world$"20")) > 0){
result <- list(world)
edge_patch <- world$"20"
#checking which one said winner is
the_winner <- character()
ifelse(edge_patch$n1 > 0 && edge_patch$n2 == 0, the_winner <- "n1", ifelse(edge_patch$n1 == 0 && edge_patch$n2 > 0, the_winner <- "n2", the_winner <- NA))
#put the results in a new row in simulation_data
simulation_data %>% add_row(worlds = result, timestep = i, parameters = list(params), winner = the_winner, r1 = params$n1$r, r2 = params$n2$r, k1 = params$n1$k, k2 = params$n2$k) -> simulation_data
return(simulation_data)
}
i <- i + 1
}
}
```
```{r clearing dataframe}
#making the dataframe empty
simulation_data <- tibble(worlds = list(), timestep = numeric(), parameters = list(), winner = character(), r1 = numeric(), r2 = numeric(), k1 = numeric(), k2 = numeric())
```
```{r running the simulation}
for(i in seq(1, 900, 1)){
simulation_data <- simulate(my_world, my_params)
}
```
```{r 1 sim isolated}
simulation_data %>% filter(!is.na(winner)) %>% group_by(winner, r1, k1, r2, k2) %>% filter(r1 == 0.5, k1 == 900, r2 == 0.5, k2 == 100) %>%
ggplot(aes(x = timestep))+
scale_fill_ordinal()+
facet_wrap(~r1+r2+k1+k2, labeller = labeller(.multi_line = FALSE, .default = label_both))+
geom_histogram(aes(x = timestep, fill = winner), binwidth = 1, colour = "black")+
labs(x = "Win Time", y = "# of Wins", title = "Comparison of Win Speed at Different r and K Values")+
theme(text = element_text(family = "Tahoma"),
axis.title = element_text(face = "bold"),
panel.background = element_rect(fill = "white"),
panel.grid = element_line(colour = "lightgrey"))
```
```{r combined}
simulation_data %>% filter(!is.na(winner)) %>% group_by(winner, r1, k1, r2, k2) %>%
ggplot(aes(x = timestep))+
scale_fill_ordinal()+
facet_wrap(~r1+r2+k1+k2, labeller = labeller(.multi_line = FALSE, .default = label_both))+
geom_histogram(aes(x = timestep, fill = winner), binwidth = 1, colour = "black")+
labs(x = "Win Time", y = "# of Wins", title = "Comparison of Win Speed at Different r and K Values")+
theme(text = element_text(family = "Tahoma"),
axis.title = element_text(face = "bold"),
panel.background = element_rect(fill = "white"),
panel.grid = element_line(colour = "lightgrey"))
```