Spatial Critter Swarming Simulation

[This article was first published on Econometrics by Simulation, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# I am interested in how small bits of individualized instructions can create collective action.

# In this simulation I will give a single instruction to each individual in the swarm.

# Choose another individual who is not too close, then accelerate towards that individual.

# I also control momentum causing the previous movement and direction to only decay at a small rate.

# TO SEE Original Script

# Critters are initially distributed randomly on a 1 x 1 grid.

ncritters = 40

xypos = matrix(runif(ncritters*2),ncol=2)
plot(xypos, main=”Critters are Initially Distributed Randomly”
          , xlab=”X”, ylab=”Y”)

# Now let’s imagine that each critter has an ideal safe distance from each other critter.

safe.dist = .3

critter.speed = .001

# If another critter is not at that safe distance than the critter will move towards the closest nearby critter.

# Let’s see how this works.

# First let’s check how close each critter is to each other critter.
# We will accomplish this by going through each critter and checking how far away each other critter is.
distances = NULL
for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

# In order to prevent critters from always chasing whatever is closest to them (and themselves) we drop anything which is closer than the safe.distance.
distances[abs(distances)closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]
  # The apply command will apply the order command to each row whiel the [1,] selects only the critter that is closes.

# Plot the
plot(xypos, xlab = “X”, ylab = “Y”)
for (i in 1:ncritters) arrows(x0=xypos[i,1], y0=xypos[i,2],

# This calculates the difference between the current position of each critter and that of the closest critter.
ab = xypos-xypos[closest,]

# To see how this is calculated, see my previous post simulating a werewolf attack.

# Now calculate the difference in the horizontal and vertical axes that the critters will move as a projection into the direction of the closest critter outside of the safe zone. = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5 = (critter.speed^^2)^.5

# This corrects the movement to ensure that the critters are flying at each other rather than away from each other.
movement = cbind( * sign(ab[,2]), * sign(ab[,1]))
between = function(xy1,xy2,point) (point>xy1&pointxy2&pointmovement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)

# Set the new xypos
xypos1 = xypos+movement

points(xypos1, col=”red”)

# ——————————————————
# Let’s turn this into an animation.


# loopnum = 100; ncritters=40; inertia = .5; show.grid=T; ani.pause=F; plot.fixed=F; plot.centered=F; brownian = F; arrow = T
flocking <- ani.pause="F," arrow="T)" brownian="F," function="" inertia=".5," loopnum="100," ncritters="40," p="" plot.centered="F," plot.fixed="F," show.grid="T,">
  # Generate xy initial positions.
  # xypos will hold the current critter position while
  # xypos0 will hold the position of the critters the previous time.
  xypos = xypos0 = matrix(runif(ncritters*2),ncol=2)-.5

  movement0 = 0

  # Loop though all of the loops.
  for (i in 1:loopnum) {

  # This specifies the range to be graphed.
  if (plot.fixed) rangex=rangey = -.5:.5
  if (!plot.fixed) {
    rangex = c(min(xypos[,1]), max(xypos[,1]))
    rangey = c(min(xypos[,2]), max(xypos[,2]))

  #  This handles the grid size when
  if (plot.centered&!plot.fixed) {
    rangex=c(-max(abs(xypos[,1])), max(abs(xypos[,1])))
    rangey=c(-max(abs(xypos[,2])), max(abs(xypos[,2])))

  # This centers the plot at the middle (0,0) if the plot width is also set to be fixed.
  if (plot.centered&plot.fixed) rangex=rangex-mean(rangex)
  if (plot.centered&plot.fixed) rangey=rangey-mean(rangey)

  # Draw critters
  plot(xypos, main=”Swarming Animation”, xlab=”X”, ylab=”Y”, axes=F, ylim=rangey, xlim=rangex, type=”p”)
  # Draw arrows.  The start of the arrows is the previous periods location.
  if (arrow&i>1) arrows(x0=xypos0[,1],y0=xypos0[,2],x1=xypos[,1],y1=xypos[,2], length = .1)

  # Show the grid in the background.
  if (show.grid) {

  # Show the origin
  text(0,0, “(0,0)”)

  # Calculate each critters distance from each other
  distances = NULL
  for (i in 1:ncritters) distances = rbind(distances, apply((xypos[i,]-t(xypos))^2,2,sum))

  # Drop those within the safe zone.
  # This selects the critter closest to the selected critter.
  closest =  matrix(1:ncritters, ncol=ncritters, nrow=ncritters)[apply(abs(distances), 1, order)[1,]]

#   distances[as.apply(!apply(distances, 1,,1,sum)==0,]=0

  # As done above
  ab = xypos-xypos[closest,] = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5 = (critter.speed^^2)^.5

  movement = cbind( * sign(ab[,2]), * sign(ab[,1]))

  between = function(xy1,xy2,point) (point>xy1&pointxy2&point
  movement = movement*(-1)^between(xypos,xypos[closest,], xypos-movement)


  movement0 = movement0*inertia + movement

  # This fancy dodad allows half of the change in movement to be due to random variation.
  if (brownian) movement0=movement0+matrix(rnorm(ncritters*2),ncol=2)*critter.speed/2

  # Set the previous round’s xy position to be equal to the current round’s.
  xypos0 = xypos

  # Update the current round’s.
  xypos = xypos+movement0

  # This is only used in the event that the animate package is in use.
  if (ani.pause) ani.pause()


# This generates a GIF animation demonstrating smoothly how these GIFs can be incorper
  ani.options(ani.width=400, ani.height=400, interval=.1)

# You must have imagemagick installed for this to work.
  saveGIF(flocking(300,100,.999, ani.pause=T), = “Swarming.gif”, replace=T)

# Here are two different graphs generated by the previous command(though the one on the bottom uses 200 frames while the one on the top uses 300)

# Let’s see how this works.
flocking(400,100,.99, plot.fixed=T)

To leave a comment for the author, please follow the link and comment on their blog: Econometrics by Simulation. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)