Jaywalking Simulation

library(caTools)         # external package providing write.gif() function

# sets up color palette funciton
jet.colors <- colorRampPalette(c("#000000", "#00FF00", "#AAAA00",
                                 "#FF0000", "#FFFFFF"))
height <- 400 # height of gif in pixels
width <- 200  # width of gif in pixels
steps <- 1000 # number of times steps

accel <- -2 # acceleration; negative is up
fric <- 0.5 # friction; resistance to movement

# initializes position array
# columns are y, x, -vy, -vx, state
# x is measured in pixels from the left
# y is measured in pixels from the top
# state is whether a person is jaywalking
the.pos <- array(c(height, 100, -4, 0, 0), dim=c(1,5))

# initializes array to hold frame informaiton
Z <- array(0, dim=c(height, width))

# draws road markers
Z[c(95:105,295:305),] <- 0.5
Z[c(195:205),(1:width)[(1:width + 10) %% 40 > 20]] <- 0.5
Z.base <- Z # sets this as the starting state for each frame

# initializes array to hold gif information
X <- array(0, dim=c(height, width, steps/2))

# draws people on the frame
draw.people <- function(pos.array, Z) {
  pos.array <- matrix(round(pos.array), ncol = 2)
  for (i in 1:dim(pos.array)[1]) {
    # gets position of person i
    y.pos <- as.integer(pos.array[i, 1])
    x.pos <- as.integer(pos.array[i, 2])

    # gets boundaries of person i's box
    top <- max(c(y.pos-1, 1))
    bottom <- min(c(y.pos+1), height)
    left <- max(c(x.pos-1), 1)
    right <- min(c(x.pos+1), height)

    # draws person i
    Z[top:bottom,left:right] <- 1
  }

  return(Z)
}

# the potential each person sees
main.der.potential <- function(x, state, t) {
  if ((t - 15) %% 500 >= 400) { # if green light, the potential is zero
    result <- 0
  } else if (x >= height - 40) { # no potential to some point before the road
    result <- 0
  } else if (x >= (height - 100) && state == 0) { # increasing potential to curb
    result <- (-x + (height - 40))/ 16
  } else if (x >= (height - 300)) { # decreasing potential to opposite side
    result <- -1
  } else { # no potential on opposite side
    result <- 0
  }
  return(result)
}

# causes people to jaywalk
change.state <- function(pos.array) {
  for (i in 1:dim(pos.array)[1]) {
    if (pos.array[i, 1] < (height - 70) && pos.array[i, 5] == 0) {
      cur.prob <- (sum(abs(the.pos[, 1] - 250) <= 50) * 100 + 1) * 0.0008
      pos.array[i, 5] <- as.integer(runif(1) < cur.prob)
    }
  }
  return(pos.array)
}

# move people depending on their speed
move.people <- function(pos.array, t) {
  pos.array <- matrix(pos.array, ncol = 5)
  for (i in 1:dim(pos.array)[1]) {
    y <- pos.array[i, 1]
    x <- pos.array[i, 2]
    dy <- pos.array[i, 3]
    dx <- pos.array[i, 4]
    dy <- dy + accel - fric * sign(dy) * sqrt(dx^2 + dy^2) +
      main.der.potential(pos.array[i, 1], pos.array[i, 5], t)
    dx <- dx - fric * sign(dx) * sqrt(dx^2 + dy^2)
    y <- y + dy
    x <- x + dx
    pos.array[i, 1:4] <- c(y, x, dy, dx)
  }
  return(pos.array)
}

# adds a new person
add.person <- function(pos.array) {
  taken <- round(pos.array[, 2])
  taken <- c(taken, taken + 1, taken + 2, taken + 3, taken - 1 , taken - 2, taken - 3)
  taken <- c(1, 2, 3, width - 2, width - 1, width, taken)
  numpos <- length((1:width)[-taken])

  if (numpos > 50) {
    new.y <- ((1:width)[-taken])[as.integer(round(runif(1) * numpos))]
    pos.array <- rbind(pos.array, c(height, new.y, -4, 0, 0))
  }

  return(pos.array)
}

# remove people who have moved off the gif
remove.past <- function(pos.array) {
  pos.array <- matrix(pos.array[pos.array[, 1] > 0, ], ncol=5)
  return(pos.array)
}

# adds the green/red sign
add.sign <- function(Z, t) {
  if ((t - 15) %% 500 >= 400) {
    Z[10:20, 180:190] <- 0.25
  } else {
    Z[20:30, 180:190] <- 0.75
  }
  return(Z)
}

for (k in 1:steps) {        # iterate through the loop
  Z <- Z.base # initializes the frame
  Z <- draw.people(the.pos[, 1:2], Z) # draws people
  Z <- add.sign(Z, k) # adds the sign
  X[,,as.integer(k/2)] <- Z # store the frame
  the.pos <- move.people(the.pos, k) # moves people for next frame
  the.pos <- remove.past(the.pos) # remove people out of the frame
  if ((dim(the.pos)[1] == 0) || runif(1) < 0.03) { # randomly or if no people
    the.pos <- add.person(the.pos) # add a new person
  }
  the.pos <- change.state(the.pos) # try to change states
}
write.gif(X, "Jaywalking.gif", col=jet.colors, delay=8)

a <- 400:1
b0 <- vapply(a, main.der.potential, double(1), t=200, state=0)
c0 <- double(400)
for (i in 2:400) {
  c0[i] <- c0[i-1] + b0[i-1]
}

b1 <- vapply(a, main.der.potential, double(1), t=200, state=1)
c1 <- double(400)
for (i in 2:400) {
  c1[i] <- c1[i-1] + b1[i-1]
}

png(filename="jaywalkpotential.png", width=400, height=400, units="px")
plot(1:400,c1, type="l", ylim=c(-275,125), col="blue", ylab="Potential", xlab="Position", main="Potential vs. Position")
lines(c0, col="red")
abline(v=100, lty=2)
abline(v=300, lty=2)
legend("topright", legend=c("normal", "jaywalking"), fill=c("red", "blue"))
dev.off()