# Spatial Critter Swarming Simulation

May 10, 2013
By

(This article was first published on Econometrics by Simulation, and kindly contributed to R-bloggers)

# 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.

# 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],
x1=xypos[closest,][i,1],
y1=xypos[closest,][i,2],
length=.1)

# 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.
a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
b.prime = (critter.speed^2-a.prime^2)^.5

# This corrects the movement to ensure that the critters are flying at each other rather than away from each other.
movement = cbind(a.prime * sign(ab[,2]), b.prime * 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.

library(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) {
abline(v=seq(-10,10,.1))
abline(h=seq(-10,10,.1))
}

# 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.
distances[abs(distances)
# 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, is.na),1,sum)==0,]=0

# As done above
ab = xypos-xypos[closest,]

a.prime = critter.speed/(1 + (ab[,2]^2)/(ab[,1]^2))^.5
b.prime = (critter.speed^2-a.prime^2)^.5

movement = cbind(a.prime * sign(ab[,2]), b.prime * sign(ab[,1]))

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

movement[is.na(movement)]=0

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), movie.name = “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()
flocking(400,100,.99)
flocking(400,100,.99, plot.fixed=T)

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...