torus.width <- 50 torus.height <- 50 agent.count <- 10 # Distance agents move per tick (second or whatever) agent.speed.per.tick <- 1 payoff.collide <- 0 payoff.move <- 1 collision.radius <- 1; # Discount factor lambda <- .90 k <- 8 n <- (2*pi)/k A <- seq(0, 2*pi - n, n) # We give all agents some random uniform starting position along the torus agent.coordinates <- matrix ( data = c(runif(agent.count, 0, torus.width), runif(agent.count, 0, torus.height)), ncol = agent.count, nrow = 2, dimnames = list( c("x","y"), c()), byrow = T) # Mean payoffs for all directions and agents theta <- matrix ( data = 1, ncol = agent.count, nrow = length(A) ) # Look in the file if you're interested source("Functions.R") # For now repeat the process 100 times max.ticks <- 100 for (payoff.move in seq(10, 100, by = 10) ){ for(epsilon in seq(0.1, 0.5, by = 0.1)) { print(paste(payoff.move, epsilon)) plotData <- matrix(data = NA, ncol = length(A), nrow = max.ticks) for(tick in 1:max.ticks) { # We use the weighted random method directions <- apply(theta, 2, function(propensities) { sample(A, size = 1, prob = propensities) }) # Give us the new coordinates for all agents if they would have moved coords.new <- t(move(torus.width, torus.height, agent.coordinates["x",], agent.coordinates["y",], agent.speed.per.tick, directions)) collisions <- sapply(1:agent.count, function(agent.id) { # Note that we use the entire column for x2 and y2 so we get a matrix of the correct size # We fix this later on collision (torus.width, torus.height, collision.radius, x1 = coords.new["x", agent.id], x2 = agent.coordinates["x",], y1 = coords.new["y", agent.id], y2 = agent.coordinates["y",]) }) # Fix what we did above and ignore agent collisions with itself diag(collisions) <- F # Reduce the matrix to a logical vector, handy for later on collisions <- apply(collisions, 2, function(collision) { any(collision) }) # Create the unit matrix with 1's for each direction chosen per agent e <- matrix(data = 0, nrow = length(A), ncol = agent.count) e[cbind(match(directions, A), 1:length(directions))] <- 1 # Create the payoff vector by looking at collisions and giving the appropriate payoff u <- diag(ifelse(collisions, payoff.collide, payoff.move)) # Update theta using a linear updating model theta <- round(linearUpdate(lambda, theta, u, e, 0), 4) plotData[tick, ] <- rowMeans(theta) agent.coordinates <- sapply(1:agent.count, function(agent.id) { if(collisions[agent.id]) { agent.coordinates[,agent.id] } else { coords.new[,agent.id] } }) } r.sums <- rowSums(theta) ratios <- r.sums / sum(r.sums) dir.index <- which(ratios >= 0.9) if(length(dir.index) > 0) { write.csv(theta, file= paste(epsilon, payoff.move, "theta.csv", sep = "-"), quote = FALSE, row.names = FALSE) write.csv(plotData, file= paste(epsilon, payoff.move, "plotData.csv", sep = "-"), quote = FALSE, row.names = FALSE) matplot(plotData, type = "l", xlab = "Time", ylab = "Mean propensity", main = "Linear Updating Model") } } }