BML_Model_WORKING_FUNCTIONS.R

###Working Functions for BML Traffic Model

##initialization matrix##
#input is r=num rows, c=num cols,p=density of matrix (% of cells with a car in it)
bml.init <- function(r, c, p){
  m<-matrix(0,nrow=r,ncol=c)
  m<-apply(m,c(1,2),function(x) ifelse(sample(c(TRUE,FALSE),1,prob=c(p,1-p)),x<-sample(c(1,2),1),x<-0))
  return(m)
}

#in case need to assign for testing purposes   assign("m",m,pos=1)


##Eastward movement (1s)##
#move.e moves can move cars INTO the border column, but cannot move border cars to "wrap around".
move.e<-function(m){
  for(i in 1:nrow(m)){
    for(j in (ncol(m)-1):1){
      if(m[i,j] == 1 & m[i,j+1] == 0){
        m[i,j]=0
        m[i,j+1]=1
      }
    }
  }
  return(m)
}

#border.e will move cars on the border to the wraparound edge, but only those cars that were originally there
#in the refrence matrix (otherwise we could theoretically have a car be moved twice if it was in the 2nd to
#last column, once by move.e and then again by border.e)
border.e<-function(m,m.ref){
  for(i in 1:nrow(m)){
    if(m.ref[i,ncol(m)]==1 & m[i,1]==0){
      m[i,ncol(m)]=0
      m[i,1]=1
    }
  }
  return(m)
}

move.e.2<-function(m,m.ref){
  for(i in 1:nrow(m)){
    for(j in (ncol(m)-1):1){
      if(m.ref[i,j] == 1 & m[i,j] == 1 & m[i,j+1] == 0){
        m[i,j]=0
        m[i,j+1]=1
      }
    }
  }
  return(m)
}

#order of operations. move.e -> border.e (m is output of move.e) -> move.e.deux
test<-matrix(cbind(c(1,0,0,1,1),c(1,1,0,0,0),c(0,0,0,1,1),c(0,0,1,1,1),c(1,0,0,0,1)),nrow=5,byrow=TRUE)

bml.step.e<-function(m){
  m.ref<-m
  m<-move.e(m)
  m<-border.e(m,m.ref)
  m<-move.e.2(m,m.ref)
  return(m)
}

##Northward movement (2s)##
#move.n can move cars north into the top border (first row), but cannot wrap them around. SECOND VERSION WORKS
move.n<-function(m){
  for(j in 1:ncol(m)){
    for(i in 2:nrow(m)){
      if(m[i,j] == 2 & m[i-1,j] == 0){
        m[i,j]=0
        m[i-1,j]=2
      }
    }
  }
  return(m)
}

#border.n wraps around those cars on the border only from the refrence matrix (should be called after move.n)
border.n<-function(m,m.ref){
  for(j in 1:ncol(m)){
    if(m.ref[1,j]==2 & m[nrow(m),j]==0){
      m[1,j]=0
      m[nrow(m),j]=2
    }
  }
  return(m)
}

move.n.2<-function(m,m.ref){
  for(j in 1:ncol(m)){
    for(i in 2:nrow(m)){
      if(m.ref[i,j] == 2 & m[i,j] == 2 & m[i-1,j] == 0){
        m[i,j]=0
        m[i-1,j]=2
      }
    }
  }
  return(m)
}

#Sundays version of the NORTHWARD movement function.
#order of operations. move.n -> border.n (m is output of move.n) -> move.n.deux
#m.ref is defined at very beginiing of the bml.step.n fucntion, see bml.step.e for refrence.

test.n<-matrix(cbind(c(2,2,0,0,2),c(2,0,0,2,2),c(2,2,0,0,0),c(2,2,2,0,0),c(2,0,0,0,2)),nrow=5,byrow=FALSE)

bml.step.n<-function(m){
  m.ref<-m
  m<-move.n(m)
  m<-border.n(m,m.ref)
  m<-move.n.2(m,m.ref)
  return(m)
}

##One complete Movement Step (1s and 2s)## 
bml.step<-function(m){
  m.grid<-m
  m<-bml.step.e(m)
  m<-bml.step.n(m)
  grid.new<-identical(m.grid,m)
  return(list(m,grid.new))
}

#Write the simulation function. The simulation function takes as its inputs the parameters of the initial
#system state (eg; size and proportion of cells with cars in them). It then outputs the initial grid, the
#final grid, and the number of steps needed to reach gridlock. It iterates through 5000 steps; if gridlock
#is reached before that, it will indicate how many steps were taken. If not, it will return 5000.

bml.sim <- function(r, c, p){
  m.initial<-bml.init(r,c,p)
  m.1<-bml.step(m.initial)
  m<-m.1[[1]]
  gridlock<-m.1[[2]]
  iterations<-5000
  for(i in 2:iterations){
    m.1<-bml.step(m)
    m<-m.1[[1]]
    gridlock<-m.1[[2]]
    if(gridlock == TRUE) break
  }
  if(i == 5000){
    i = "System gridlock free after 5000 iterations"
  }
  else if(i != 5000){
    i = paste(paste("System reached gridlock after", i, sep=" "),"iterations",sep=" ")
  }
  return(list(m.initial,m,i))
}