R snippets

Synchronous vs. asynchronous agent activation example

2012-02-24T03:09:00.001-08:00

This time I have implemented NetLogo Voting model to verify how agent activation scheme influences the results.
The code executing the simulation is given below. It simulates two types of voter preferences encoded as 1 and -1. In this way average preference equal to 0 indicates 50/50 split. Voters are arranged on square grid with vertical and horizontal wrapping.

nei8 <- function(x, y, size) {

base.x <- c(-1, -1, -1, 0, 0, 1, 1, 1)

base.y <- c(-1, 0, 1, -1, 1, -1, 0, 1)

1 + ((cbind(x + base.x, y + base.y) - 1) %% size)

}

step.syn <- function(space) {

size <- nrow(space)

new.space <- space

for (x in 1:size) {

for (y in 1:size) {

nei.pref <- sum(space[nei8(x, y, size)])

if (nei.pref > 0) { new.space[x, y] <- 1 }

if (nei.pref < 0) { new.space[x, y] <- -1 }

}

return(new.space)

}

step.asyn <- function(space) {

size <- nrow(space)

old.space <- space

all.x <- rep(1:size, size)

all.y <- rep(1:size, each = size)

dec.seq <- sample.int(length(all.x))

for (i in dec.seq) {

x <- all.x[i]

y <- all.y[i]

nei.pref <- sum(space[nei8(x, y, size)])

if (nei.pref > 0) { space[x, y] <- 1 }

if (nei.pref < 0) { space[x, y] <- -1 }

}

return(space)

}

simulate <- function (size, is.syn, do.plot) {

x <- rep(1:size, size)

y <- rep(1:size, each = size)

space <- 2 * matrix(rbinom(size ^ 2, 1, 0.5), nrow = size) - 1

rep <- 0

while (TRUE) {

rep <- rep + 1

old.space <- space

if (is.syn) {

space <- step.syn(space)

} else {

space <- step.asyn(space)

}

if (do.plot) {

par(pin=c(3,3))

plot(x , y, axes = FALSE, xlab="", ylab="",

col = space + 2, pch = 15, cex = 40 / size)

}

if (all(old.space == space)) {

return(c(rep, mean(space)))

}

if (rep > 100) {

return(c(NA, mean(space)))

}

It uses either step.syn function which assumes that all voters make decision at the same time in each time step and step.asyn where all voters are activated in random order.

To reproduce graph similar to NetLogo original the do.plot option should be set to TRUE. For example simulate(32, FALSE, TRUE) generates the following picture:

In order to compare synchronous and asynchronous voter activation regimes I have run the following code:

s.r <- replicate (1024, simulate(32, TRUE, FALSE))

as.r <- replicate (1024, simulate(32, FALSE, FALSE))

par(mfrow = c(1, 2))

boxplot(cbind("synchronous"=s.r[1,],"asynchronous"=as.r[1,]),

main = "time")

boxplot(cbind("synchronous"=s.r[2,],"asynchronous"=as.r[2,]),

main = "mean preference")

As shown on the picture below the distribution of mean voter preferences is very similar. However convergence speed of both methods varies greatly. Moreover under synchronous activation around 0.5% of simulations do not reach stable state.

Generating directed Watts-Strogatz network

2012-02-11T03:38:00.000-08:00

There are two limitations of Watts-Strogatz network generator in igraph package: (1) it works only for undirected graphs and (2) rewiring algorithm can produce loops or multiple edges.

You can use simplify function of such a graph, but then number of edges in the graph is reduced.
Below I give ws.graph function that generates directed graph without these problems:

library(igraph)
library(colorspace)

resample <- function(x, ...) {
    x[sample.int(length(x), ...)]
}

ws.graph <- function(n, nei, p) {
    stopifnot(nei < n)
    edge.list <- vector("list", n)
    for (v in 0:(n-1)) {
        edge.end <- union((v + 1:nei) %% n,
                          (v + (-1:-nei)) %% n)
        rewire <- (runif(length(edge.end)) < p)
        edge.end <- edge.end[!rewire]
        rewired <- resample(setdiff(0 : (n-1),
             c(edge.end, v)), sum(rewire))
        edges <- rep(v, 4 * nei)
        edges[c(F, T)] <- c(edge.end, rewired)
        edge.list[[v + 1]] <- edges
    }
    graph(unlist(edge.list))
}

n <- 8
nei <- 2
p.levels <- c(0, 0.25, 0.5, 1)
reps <- 2^16
m <- matrix(0, nrow = n, ncol = n)
m <- list(m, m, m, m)
for (i in 1:reps) {
    for (j in seq_along(p.levels)) {
        g <- ws.graph(n, nei, p.levels[j])
        m[[j]] <- m[[j]] + get.adjacency(g)
    }
}

x <- rep(1:n, n)
y <- rep(1:n, each = n)
par(mfrow = c(2, 2), mar= c(5, 5, 2, 2))
for (i in 1:4) {
    mc <- as.vector(m[[i]]) / reps
    mc <- cbind(mc, mc, mc)
    mc <- 1 - mc
    plot(x, y, col = hex(RGB(mc)),  pch = 19, ylab = "",
         xlab = paste("p =", round(p.levels[i], 4)), cex = 1.5)
}

This is the resulting plot:

As expected increasing rewiring probability to 1 makes edge probability distribution more uniform.

Implementing Circles example

2012-02-04T08:01:00.000-08:00

This week I reimplemented part of Conic Sections 1 model from NetLogo. In the model turtles seek to to be in target distance from center.

My code takes only one center point, so only circles can be obtained. Apart from turtle location plot given in NetLogo implementation I added:

plot showing maximal difference between turtle distance and target distance;
decreasing turtle step size.

Here is the plot showing final simulation state, but it is also nice to watch the simulation run:

Below is the code generating the simulation:

# n: number of turtles
# p.x, p.y: location of center
# range: turtles have random position from [0,range]
#        and will move in random angle a
# step: how fast turtles move
# target: target distance from center
# time: simulation time
init <- function(n, p.x, p.y, range, step, target, time) {
    sim <- list(
        turtles = data.frame(x = runif(n, max = range),
                             y = runif(n, max = range),
                             a = runif(n, max = 2 * pi)),
        p.x = p.x, p.y = p.y, step = step, target = target,
        time = time, max.dist = rep(NA, time))

    # Calculate turtle distance from center
    sim$turtles$dist <- sqrt((sim$turtles$x - p.x) ^ 2 +
                             (sim$turtles$y - p.y) ^ 2)
    return(sim)
}

step <- function(sim) {
    x <- sim$turtles$x
    y <- sim$turtles$y

    # Remember last distance and save current distance
    o.dist <- sim$turtles$dist
    n.dist <- sqrt((x - sim$p.x) ^ 2 + (y - sim$p.y) ^ 2)
    sim$turtles$dist <- n.dist

    # For turtles that are too far and are moving out
    # or too close and are moving in randomly change direction
    w.dist <- ((n.dist < o.dist) & (n.dist < sim$target)) |
               ((n.dist > o.dist) & (n.dist > sim$target))
    sim$turtles$a[w.dist] <- runif(sum(w.dist), max = 2 * pi)

    sim$turtles$x <- x + sin(sim$turtles$a) * sim$step
    sim$turtles$y <- y + cos(sim$turtles$a) * sim$step
    return(sim)
}

do.plot <- function(sim) {
    rng <- quantile(c(sim$turtles$x, sim$turtles$y),
                    c(0.05, 0.95))
    rng <- round(rng, -1) + c(-10, 10)
    par(mai = rep(0.5, 4), mfrow = c(1, 2))
    plot(sim$turtles$x, sim$turtles$y, pch = ".",
         xlim = rng, ylim = rng, xlab = "", ylab = "",
         main = "Turtle location")
    points(sim$p.x, sim$p.y, col = "red", pch = 20, cex = 2)
    plot(sim$max.dist, type = "l",
         ylim = c(0, max(sim$max.dist, na.rm = TRUE) + 5),
         xlab = "", ylab = "", main = "Max difference from target")
}

run <- function(sim) {
    for (i in 1:sim$time) {
        sim <- step(sim)
        sim$step <- sim$step * 127 / 128
        sim$max.dist[i] <- max(sim$turtles$dist) - sim$target
        do.plot(sim)
    }
}

sim <- init(4096, 128, 128, 256, 2, 128, 512)
set.seed(0)
run(sim)

Replicating NetLogo Fire model

2012-01-27T10:27:00.000-08:00

While preparing for the new semester I have started reimplementing standard NetLogo examples in R. The first is Fire model.
The simulation in R is presented here:

# Forest matrix trees encoding:

# 3 – green, 2 – burning, 1 – burnt, 0 – no tree

simulate <- function (forest.density = 0.6, size = 251) {

forest <- matrix(3 * rbinom(size ^ 2, 1, forest.density),

nrow = size)

forest[1,] <- 2

while (sum(forest == 2) > 0) {

list.red <- which(forest == 2)

# find spots touching burning trees

ignite.x <- 1 + floor((list.red - 1) / size)

ignite.x <- rep(ignite.x, each = 4) +

rep(c(1, -1, 0,0), length(ignite.x))

ignite.y <- 1 + ((list.red - 1) %% size)

ignite.y <- rep(ignite.y, each = 4) +

rep(c(0, 0, 1,-1), length(ignite.y))

include <- (ignite.x > 0) & (ignite.x <= size) &

(ignite.y > 0) & (ignite.y <= size)

ignite.n <- (ignite.x[include] - 1) * size +

ignite.y[include]

# ignite green trees

is.green <- (forest[ignite.n] == 3)

forest[ignite.n[is.green]] <- 2

# trees burn for one period

forest[list.red] <- 1

}

return(1 - sum(forest[-1,] == 3) / sum(forest[-1,] > 0))

}

It is interesting to notice that it is quite compact. The only exception is handling of finding neighbors which is much easier in NetLogo.

The only thing left was to check that the results produced by NetLogo and R implementations are identical. I have simulated both models 1000 times for forest density equal to 60%. Here is the plot comparing percent burned density estimates and Kolmogorov-Smirnov test p-value.

It can be concluded that for the parametrization I have selected models produce similar results.