Friday 31 October 2014

Conflict data for October

Here's a quick animation of the Armed Conflict Location Event Database's realtime data for October.




I created this with the code:

acl <- read.csv("/Users/david/ACLEDoct27.csv")
acl$date <- as.Date(acl$EVENT_DATE, "%d-%b-%y")
acl <- acl[acl$date > as.Date("2014-10-01"),]


library(anim.plots)
library(maps)

tmp <- anim.points(LATITUDE ~ LONGITUDE + as.numeric(date),  

      data=acl , speed=1, col=rgb(1,0,0,.7), pch=19,
      cex=log(TOTAL_FATALITIES+1), show=F)

saveGIF(
  replay(tmp,
  before=map('world', fill=T, col=terrain.colors(5)[2:5],

     xlim=c(-20, 60), ylim=c(-40, 40)),
  after=legend("bottomleft", col="red", pt.cex=log(1+c(1,10,100)) ,
      legend=c(1,10,100), pch=19, y.intersp=1.5, x.intersp=1.5, 

      bty="n", title="Fatalities")),
  "conflict.gif")


Easy animations in R using anim.plots


After the last two R posts, here's something a bit more useful.

Creating animations in R was hard. Then Yihui Xie wrote the animation package. animation lets you plot individual frames, then record them.

I've been working on a little R package called anim.plots which adds some syntactic sugar to this. It provides a few commands, similar to the standard R graphics, to create animated plots. Here's a quick demo. (The speed is quite slow: this is because I converted the files to animated GIFs. The original animations go reasonably fast.)
library(anim.plots)
anim.plot(1:5, 1:5, col="green") 
To control what gets plotted when, use the times parameter.

x <- rep(1:100/10, 20)
times <- rep(1:20, each=100) # twenty frames with 100 points each
y <- sin(x*times/4)
waves <- anim.plot(x,y,times, type="l", col="orange", lwd=2, speed=2)
You can make incremental animations using the window parameter. Here's a printout of the first 20 plot symbols:
anim.plot(rep(1:10,2), rep(2:1, each=10), window=1:t, pch=1:20, ylim=c(0,3),cex=2, col=1:5, xlab=paste("Plot symbols", 1:20))

The above code also shows how parameters get recycled where appropriate, either to the number of points, or to the number of frames. For more complex parameters, you might have to use a matrix. Here we zoom in on a distribution of points by changing the xlim and ylim parameters.
 x <- rnorm(4000)
 y <- rnorm(4000)
 x <- rep(x, 40)
 y <- rep(y, 40)
 xlims <- 4*2^(-(1:40/10))
 ylims <- xlims <- rbind(xlims, -xlims)
 anim.plot(x, y, times=40, speed=5, xlim=xlims, ylim=ylims,
       col=rgb(0,0,0,.3), pch=19, bg="white")
The window argument is quite powerful. You can use it to create moving plots over time.


## discoveries 1860-1959
xlim <- rbind(1860:1959,1870:1969)
anim.plot(1860:1959, discoveries, times=1:100, xlim=xlim, col="red", xaxp=rbind(xlim, 10), window=t:(t+10), type="h", lwd=8, speed=5)

There's also a formula interface. Here's a plot of some chicks being fed one of four different diets:

data(ChickWeight)
ChickWeight$chn <- as.numeric(as.factor(ChickWeight$Chick))
 
tmp <- anim.plot(weight ~ chn + Time, data=ChickWeight, col=as.numeric(Diet), pch=as.numeric(Diet), speed=3)



Sometimes you need to run extra plotting commands after each frame. You can do this using the replay command:

replay(tmp, after=legend("topleft", legend=paste("Diet", 1:4), pch=1:4, col=1:4))
 You aren't limited to the standard plot function. Here's a histogram with increasingly fine bins:
anim.hist(rep(rnorm(5000), 7), times=rep(1:7, each=5000), breaks=c(5,10,20,50,100,200, 500, 1000))
 And here's how to animate a plot of a mathematical expression:
anim.curve(x^t, times=10:50/10, n=20)
 You can animate contour plots. Here's the map of the Maunga Whau volcano in Auckland, emerging from a sphere. (Well, I had to come up with something!)
data(volcano)
# create a circle:
tmp <- volcano
tmp[] <- 200 - ((row(tmp) - 43)^2 + (col(tmp) - 30)^2)/20
# an animated contour needs a 3D array:
cplot <- array(NA, dim=c(87,61,20))
cplot[,,1] <- tmp
cplot[,,20] <- volcano
# morph the volcano into a circle:
cplot <- apply(cplot, 1:2, function(x) seq(x[1], x[20], length.out=20))
cplot <- aperm(cplot, c(2,3,1))
anim.contour(z=cplot, times=1:20, speed=3, levels=80 + 1:12*10, lty=c(1,2,2))
Finally, here are last week's earthquakes on the world map.
eq = read.table(
     file="http://earthquake.usgs.gov/earthquakes/catalogs/eqs7day-M1.txt",
     fill=TRUE, sep=",", header=T)
     eq$time <- as.numeric(strptime(eq$Datetime, "%A, %B %d, %Y %X UTC"))
     eq <- eq[-1,]
library(maps)
map('world')
tmp <- anim.points(Lat ~ Lon + time, data=eq, cex=Magnitude, col=rgb(
       1-Depth/200, 0, Depth/200,.7), pch=19, speed=3600*12, show=FALSE)  
replay(tmp, before=map('world', fill=TRUE, col="wheat"))
Code is available at http://github.com/hughjonesd/anim.plots.  If you have devtools installed you can get it with:
install_github("hughjonesd/anim.plots")


Saturday 25 October 2014

Screw UKIP

... and Cameron too!

I'm sick of this anti-European bullshit, and I'm going to stand up against it the only way I know how.

par(bg=rgb(0,51/255,153/255))
symbols(sin(1:12/12*2*pi), cos(1:12/12*2*pi),
  stars=matrix(rep(c(0.52,1,0.52, 0.4)*1/7, 10), nrow=12, ncol=20, byrow=T),  
  fg=rgb(255/255,204/255,0), bg=rgb(255/255,204/255,0), inches=F, axes=F, ann=F, crt=90)

The proportions are guesses - couldn't be bothered to do the math for what is described here. The colours are right though. I had to use a 10 pointed star to make sure the stars were the right way up.

Tuesday 21 October 2014

Simulating infections


Here's a little R simulation of fashion/disease infection/group identity in a population. There are two groups, red and white. Individuals move around at random. Sometimes, they change group. Whether they change depends on the number of their neighbours who are from a different group. Also, individuals who have just changed are "infectious" (or, trendsetters). Individuals who have not changed for a long time are not very infectious, and are also easy to infect.



It's fun to play with the parameters. It would also be fun to implement this with real people, as e.g. a mobile app, using real time movement.

##############################################
N <- 500
maxdist <- .1 # distance within which "neighbours" are counted
maxt <- 240   # number of repetitions
delay <- 0    # time delay between reps
step <- FALSE # ask after each step?
speed <- 0.01 # speed indivs move
decay <- 1    # speed decay
cols <- c("red", "white", "pink", "yellow") # last two for changers
##############################################
 prop2changeprob <- function(props) {
  props^2
}

infectivity <- function(lifespan) 1/lifespan
dist <- function(ind, oths) sqrt((ind$x-oths$x)^2+(ind$y-oths$y)^2)
par(bg="gray")

inds <- data.frame(x=runif(N), y=runif(N), col=sample(1:2, N, replace=TRUE))
inds$oldx <- inds$x
inds$oldy <- inds$y
inds$lifespan <- 0 # since last change

for (t in 1:(maxt)) {
  inds$lifespan <- inds$lifespan + 1
  propother <- sapply(1:nrow(inds), function(i) {
    i <- inds[i,]
    close <- dist(i, inds) < maxdist
    if (any(close)) {
      wt <- infectivity(inds$lifespan)
      wt[!close] <- 0
      sum(wt * (inds$col!=i$col)) / sum(wt) # weighted by
    } else 0
  })
  change <- runif(N) < prop2changeprob(propother)
  displaycol <- inds$col
  inds$col <- ifelse(change, 3-inds$col, inds$col)
  inds$lifespan[change] <- 0
  displaycol <- ifelse(displaycol==inds$col, displaycol, inds$col+2)
  inds$oldx <- inds$x
  inds$oldy <- inds$y
  inds$x <- inds$x + runif(N, -speed,speed) + decay*(inds$x-inds$oldx)
  inds$y <- inds$y + runif(N, -speed,speed)+ decay*(inds$y-inds$oldy)
  #inds$x <- pmax(0, pmin(1, inds$x)) # square
  #inds$y <- pmax(0, pmin(1, inds$y))
  inds$x <- inds$x %% 1 # torus
  inds$y <- inds$y %% 1
 
  plot(inds$x, inds$y, col=cols[displaycol], pch=19, cex=.7, bg="grey",
    xlim=0:1, ylim=0:1)
  Sys.sleep(delay)
  if (step) {
    rl <- readline()
    if (rl=="q") stop("Quitting")
  }
}