Before we begin modeling we will add some functions that can be re-used, in R studio open a new R script, copy and past the following code and save as my_lib.R.
library(plot.matrix)
library(ggplot2)
library(dplyr)
library(reshape2)
library(sf)
simple_fill <- function(mat, pct) {
mat[] <- 1
for (m in 1:nrow(mat)){
for (n in 1:ncol(mat)){
if (sample(1:100, 1) <= pct){
mat[m, n] <- 2
}
}
}
mat
}
complex_fill <- function(mat, factors, pct) {
if(sum(pct) > 100){
stop("Asked for > 100%% feature coverage")
}
pct2 <- c(100-sum(pct), pct)
pct2 <- cumsum(pct2)
mat[] <- 1
for (n in 1:ncol(mat)){
for (m in 1:nrow(mat)){
rnd = sample(1:100, 1)
for(i in length(pct2):1)
if (rnd <= as.integer(pct2[i])){
mat[m, n] <- i
}
}
}
mat
}
plot_grid <- function(grid, labels, colours) {
labels <- c('blank', labels)
colours <- c('white', colours)
m <- melt(grid) %%>%% mutate(feat = labels[value])
pal <- data.frame(t(colours))
colnames(pal) <- labels
ggplot(m, aes(x = Var1, y = Var2)) + geom_raster(aes(fill = feat)) + scale_fill_manual(values = pal)
}
get_neighbors <- function(tmp_grid, x, y, r){
min_x <- max(1, x-r)
max_x <- min(nrow(tmp_grid), x+r)
min_y <- max(1, y-r)
max_y <- min(ncol(tmp_grid), y+r)
neighbors <- list()
for (n in min_y:max_y){
for (m in min_x:max_x){
if (!(m == x && n == y)){
neighbors <- cbind(neighbors, tmp_grid[m, n])
}
}
}
unlist(neighbors)
}
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
apply_rule <- function(mat, fun) {
grid <- mat
for (n in 1:ncol(mat)){
for (m in 1:nrow(mat)){
grid[m, n] <- fun(mat, m, n)
}
}
grid
}
create_agents <- function(features, colours, num_agents){
a <- data.frame()
for (i in 1:length(num_agents)){
a <- rbind(a, data.frame(replicate(2, sample(1:100, num_agents[i], rep=T)), features[i]))
}
a <- cbind(a, seq(1:sum(num_agents)))
colnames(a) <- c('x', 'y', 'type', 'id')
agents <- st_as_sf(a, coords = c('x', 'y'))
}
drop_agent <- function(agents, agent_id){
agents <- agents %%>%% filter(id != agent_id)
}
add_agent <- function(agents, target_type){
new_agent <- data.frame(x = sample(1:100, 1), y = sample(1:100, 1), type = target_type, id = max(agents$id) + 1)
agents <- rbind(agents, st_as_sf(new_agent, coords = c('x', 'y')))
}
plot_agents <- function(agents, labels, colours) {
ggplot(agents) + geom_sf(aes(colour = type)) + xlim(c(0, 100)) + ylim(c(0, 100)) + scale_colour_manual(values = colours)
}
get_nearest <- function(agents, target_id){
target <- agents %%>%% filter(id == target_id)
others <- agents %%>%% filter(id != target_id)
near <- st_nearest_feature(target, others)
others[near, ]$id
}
get_nearest_type <- function(agents, target_id, target_type){
target <- agents %%>%% filter(id == target_id)
others <- agents %%>%% filter(id != target_id, type == target_type)
near <- st_nearest_feature(target, others)
others[near, ]$id
}
get_distance <- function(agents, a, b){
A <- agents %%>%% filter(id == a)
B <- agents %%>%% filter(id == b)
as.numeric(st_distance(A, B))
}
move_agent <- function(agents, target_id, distance, direction){
target <- agents[agents$id == target_id, ]
x <- st_coordinates(target)[[1]]
y <- st_coordinates(target)[[2]]
x <- x + (distance * sin(direction * (pi/180)))
y <- y + (distance * cos(direction * (pi/180)))
if (x > 100){
x = 100
}
if (x < 0){
x = 0
}
if (y > 100){
y = 100
}
if (y < 0){
y = 0
}
agents[agents$id == target_id, ]$geometry <- st_geometry(st_point(c(x, y)))
agents
}
apply_agent_rule <- function(agen, fun) {
for (id in agen$id){
agen <- fun(agen, id)
}
agen
}
With the custom library added start Your modeling script; in this file instead of loading libraries we will load the R script just created with the source command.
source('~/L/GEOG413/lab10/my_lib.R')
Cell Based Automata
We will start with a cell based modeling approach, with the framework defined by the code we just added build a grid with 2 types of features, you can name and colour the feature types, then build a grid by size.
The grid will start empty, simple_fill will set a a percentage of a matrix to the value 2 (1 is blank so 2 is our first class).
# Label Classes, and pick colour for display, must be same length
features <- c('sheep', 'wolf')
colours <- c('blue', 'red')
# Define a grid size, first arg is fill balue, always defaults to blank
grid <- matrix(1, ncol=50, nrow=50)
# Simple fill resets grid to feature 1, fills with desired percent feature 2
grid <- simple_fill(mat=grid, pct=10) plot_grid(grid, features, colours)
# Example of how to show current model
plot_grid(grid, features, colours)
When using more than one feature there is also the option of complex fill. In this case pct is a list with the percentage of each feature desired, the remainder will be blank, in this case 50% sheep, 30% wolf, leaving 20% blank tiles.
# Complex fill takes a grid to be filled, and the percentage of each class grid <- complex_fill(mat=grid, pct=c(50, 30)) plot_grid(grid, features, colours)
With a grid defined and populated with data the next step is to define the rules of the model.
# Reduces each feature value by one, modulo for wrap around, add 1 to solve 0 index
rule1 <- function(grid, x, y){
(grid[x, y] + 1) %%%% (length(features)+1) + 1
}
# Example Stepping through rule 1
grid <- apply_rule(grid, rule1)
plot_grid(grid, features, colours)
Rules can also be dependent on neighbors; the get_neighbors() function takes the overall grid, as well as x, and y of the current cell, the final argument is the search radius. A radius of 1 will return up to 8 neighbors, while 2 would return up to 24 neighbors.
# Makes each cell become like average neighbor
rule2 <- function(grid, x, y){
# Note Mode has capital not part of base R.
Mode(get_neighbors(grid, x, y, 1))
}
# Example stepping through rule 2
grid <- apply_rule(grid, rule2)
plot_grid(grid, features, colours)
We can run the simulation for more than a single round by using a loop; sys.sleep is used to give the plot a chance to render before the next iteration.
grid <- complex_fill(mat=grid, pct=c(50, 30))
for (i in 1:10){
grid <- apply_rule(grid, rule2)
print(plot_grid(grid, features, colours))
Sys.sleep(1)
}
Agent Based Modeling
Using cell based automata can be limiting in terms of resolution and rules. Agents are able to have more degrees of movement, as well as the ability for different agents to play by different sets of rules. These differences may be based on agent type or could be random traits within a population of agents.
We will start by making 200 sheep and 5 wolves
# Create new set of agents
features <- c('sheep', 'wolf')
colours <- c('blue', 'red')
num_agents <- c(200, 5)
agents <- create_agents(features, colours, num_agents)
# Plot Agents
plot_agents(agents, features, colours)
There are a number of functions made available by my_lib.R
# Add new Agent agents <- add_agent(agents, 'sheep') # Delet Agent by ID agents <- drop_agent(agents, 160) # Plot Agents plot_agents(agents, features, colours) # Move Agent by ID # List of agents, ID, distance, direction in degrees agents <- move_agent(agents, 3, 15, 270) # Get information about neighbors get_distance(agents, 1, 2) # Returns Distance between agents 1 & 2 get_nearest(agents, 1) # Returns ID of agent closest to agent 1 get_nearest_type(agents, 1, 'wolf') #Returns closes wolf to agent 1
A simple simulation is to move agents a random distance and direction
# Example Rule
move_random <- function(agents, t_id){
move_agent(agents, t_id, sample(1:10, 1), sample(1:360, 1))
}
# Apply_agent applies rule to all agents, initative is based on id
agents <- apply_agent_rule(agents, move_random)
plot_agents(agents, features, colours)
A more complex example
# More complex Rule
eat_sheep <- function(agent, t_id){
if (agent[agent$id == t_id, ]$type == 'wolf'){
nearest <- get_nearest_type(agent, t_id, 'sheep')
distance <- get_distance(agent, t_id, nearest)
if (distance < 15){
agent <- drop_agent(agent, nearest)
}
}
agent
}
# Dinner Time
for (frame in 1:10){
for (sub_frame in 1:5){ # Run five cycles between plotting for faster performance
agents <- apply_agent_rule(agents, eat_sheep)
agents <- apply_agent_rule(agents, move_random)
}
print(plot_agents(agents, features, colours))
Sys.sleep(1)
}
Assignment
Question 1 (1.5 points): Implement Conway’s Game of life. The rules for the Game of life are as follows
1. A cell with less than 2 neighbors will die of loneliness
2. A cell with 2 or 3 neighbors will come alive
2. A cell with more than 4 neighbors will die of crowding
Hint: Assuming live cells are the first feature (2 after accounting for blank) the following will get the number of live neighbors.
n <- sum(get_neighbors(grid, x, y, 1)==2)
Question 2 (3.5 points): Building upon the previous model extend with these additional rules
- Add a 3rd agent type of grass
- Add 10 grass per cycle
- If a wolf does not eat it dies
- If two wolves are within 10 distance add another wolf (this will result in 2 wolves being added due to agent order rules)
- If two sheep are within 5 distance and within 5 of grass, eat the grass and add a sheep