Cellular automaton - Game of Life using {stars} - generated using ChatGPT

game_of_life_stars.R

# Load the necessary library
library(stars)

# Set seed for reproducibility
set.seed(123)

# Define grid dimensions
nx <- 100  # Number of columns (x-axis)
ny <- 100  # Number of rows (y-axis)

# Define the number of simulation steps
num_steps <- 100

# Initialize the grid with random states (0: dead, 1: alive)
initial_matrix <- matrix(
  sample(c(0, 1), nx * ny, replace = TRUE, prob = c(0.7, 0.3)),
  nrow = ny, ncol = nx
)  # rows correspond to y-axis

# Initialize a 3D array to store all steps
# Dimensions: y, x, step
# Initialize all steps to 0
all_steps_array <- array(
  0,
  dim = c(ny, nx, num_steps + 1),
  dimnames = list(
    y = 1:ny,
    x = 1:nx,
    step = 0:num_steps
  )
)

# Assign the initial state to the first slice (step 0)
all_steps_array[ , , 1] <- initial_matrix

# Function to shift the matrix in a specified direction
shift_matrix <- function(mat, direction) {
  nr <- nrow(mat)
  nc <- ncol(mat)
  
  if (direction == "up") {
    shifted <- rbind(mat[-1, ], rep(0, nc))
  } else if (direction == "down") {
    shifted <- rbind(rep(0, nc), mat[-nr, ])
  } else if (direction == "left") {
    shifted <- cbind(mat[, -1], rep(0, nr))
  } else if (direction == "right") {
    shifted <- cbind(rep(0, nr), mat[, -nc])
  } else if (direction == "up_left") {
    shifted <- rbind(mat[-1, -1], rep(0, nc - 1))
    shifted <- cbind(shifted, rep(0, nr))
  } else if (direction == "up_right") {
    shifted <- rbind(mat[-1, -nc], rep(0, nc - 1))
    shifted <- cbind(rep(0, nr), shifted)
  } else if (direction == "down_left") {
    shifted <- rbind(rep(0, nc - 1), mat[-nr, -1])
    shifted <- cbind(shifted, rep(0, nr))
  } else if (direction == "down_right") {
    shifted <- rbind(rep(0, nc - 1), mat[-nr, -nc])
    shifted <- cbind(rep(0, nr), shifted)
  } else {
    stop("Invalid direction")
  }
  
  shifted
}

# Function to count alive neighbors
count_alive_neighbors <- function(mat) {
  # Shift the matrix in all eight directions
  up        <- shift_matrix(mat, "up")
  down      <- shift_matrix(mat, "down")
  left      <- shift_matrix(mat, "left")
  right     <- shift_matrix(mat, "right")
  up_left   <- shift_matrix(mat, "up_left")
  up_right  <- shift_matrix(mat, "up_right")
  down_left <- shift_matrix(mat, "down_left")
  down_right<- shift_matrix(mat, "down_right")
  
  # Sum all the shifted matrices to get the count of alive neighbors
  neighbor_count <- up + down + left + right + up_left + up_right + down_left + down_right
  neighbor_count
}

# Simulation loop
for (step in 1:num_steps) {
  
  # Extract the current matrix from the previous step
  current_matrix <- all_steps_array[ , , step]
  
  # Count alive neighbors
  neighbor_count <- count_alive_neighbors(current_matrix)
  
  # Apply Game of Life rules to determine the next state
  # Rule 1: Any live cell with two or three live neighbors survives.
  survive <- (current_matrix == 1) & (neighbor_count == 2 | neighbor_count == 3)
  
  # Rule 2: Any dead cell with exactly three live neighbors becomes a live cell.
  birth <- (current_matrix == 0) & (neighbor_count == 3)
  
  # Rule 3: All other live cells die in the next generation. Similarly, all other dead cells stay dead.
  next_matrix <- matrix(0, nrow = ny, ncol = nx)
  next_matrix[survive | birth] <- 1
  
  # Assign the next state to the corresponding slice in the array
  all_steps_array[ , , step + 1] <- next_matrix
  
  # Optional: Print progress
  if (step %% 10 == 0) {
    cat("Completed step", step, "\n")
  }
}

# Convert the 3D array into a stars object with dimensions x, y, step
ca_stars <- st_as_stars(all_steps_array)

# Assign proper dimension names and units if necessary
# (Optional: Modify as per your requirements)
st_dimensions(ca_stars) <- st_dimensions(
  x = 1:nx,
  y = 1:ny,
  step = 0:num_steps
)

# Assign dimension names
names(st_dimensions(ca_stars)) <- c("x", "y", "step")

# Plot all steps
plot(ca_stars, col = c("white", "black"), axes = FALSE)

# Animate the simulation using a loop to plot each step
for (step in 0:num_steps) {
  plot(ca_stars[,,,step + 1], main = paste("Step", step), col = c("white", "black"), axes = FALSE)
  Sys.sleep(0.1)  # Adjust the pause duration as needed
}