CSCI-540: Chapter 8: C code

.gitignore

*.so
*.o
Rplots.pdf

README.md

C code necessary to run BML model

git clone https://gist.github.com/6174858c3ec9076e130b.git BML
cd BML
make
make run

The BML.R file should contain all the code necessary to run some simulations.

Code Results: Video on Screencast.com

BML.c

#include "BML.h"

unsigned int
bml_move(int *grid, int *dims, int *locations, int *num, Direction dir, int speed)
{
	int i;
	int r, c, rprime, cprime;
	unsigned int numMoved = 0;

	for (i = 0; i < *num; i++) {
		r = locations[i] - 1;
		c = locations[i + *num] - 1;

		if (dir == RED) {
			rprime = r;
			cprime = c + speed > dims[1] - 1 ? c + speed - dims[1] : c + speed; // (c + speed) % dims[1];
		} else {
			rprime = r + speed > dims[0] - 1 ? r + speed - dims[0] : r + speed; // (r + speed) % dims[0];
			cprime = c;
		}

		if (grid[rprime + dims[0] * cprime] == 0) {
			grid[r + dims[0] * c] = 0;
			grid[rprime + dims[0] * cprime] = dir;
			if (dir == RED)
				locations[i + *num] = cprime + 1;
			else
				locations[i] = rprime + 1;
			numMoved++;
		}
	}

	return numMoved;
}

BML.h

// These are symbolic constants for specifying the direction.
typedef enum {RED = 1, BLUE = 2} Direction;

// This is the function that does the actual moving of the cars for a given time step.
unsigned int bml_move(int *grid, int *dims, int *locations, int *num, Direction dir, int speed);

#include <stdio.h>
#include <string.h>

BML.R

getCarLocations = function(g) {
	w = (g != "")
	i = row(g)[w]
	j = col(g)[w]
	pos = cbind(i, j)
	structure(pos, dimnames = list(g[pos], c("i", "j")))
}

createGrid = function(dims = c(100, 100), numCars = .3) {
	if (length(dims) == 1) dims = rep(dims, 2)
	if (length(numCars) == 1 && numCars < 1)
		numCars = rep(prod(dims) * numCars/2, 2)
	if (length(numCars) == 1) numCars = rep(numCars, 2)

	grid = matrix("", dims[1], dims[2])
	pos = sample(1:prod(dims), sum(numCars))
	grid[pos] = sample(rep(c("red", "blue"), ceiling(numCars)))[seq(along = pos)]
	class(grid) = c("BMLGrid", class(grid))
	grid
}

plot.BMLGrid = function(x, ...) {
	if (typeof(x) == "character")
	  z = matrix(
	    match(x, c("", "red", "blue")),
	    nrow(x), ncol(x)
	  )
	else
		z = x
  image(
    t(z), col=c("white", "red", "blue"),
    axes=FALSE, xlab="", ylab="", ...
  )
  box()
}

dyn.load("BML.so")

crunBML = function(grid, numIter = 100L) {
	stopifnot(is.loaded("R_BML"))
	z = class(grid)
	gi = matrix(match(grid, c("red", "blue"), 0L), nrow(grid), ncol(grid))
	pos = getCarLocations(grid)
	red = pos[ rownames(pos) == "red", ]
	blue = pos[ rownames(pos) == "blue", ]

	velocity = matrix(0L, as.integer(numIter), 2L, dimnames=list(NULL, c("red", "blue")))

	val = .C("R_BML", grid = gi, dim(gi),
					red = red, nrow(red),
					blue = blue, nrow(blue),
					as.integer(numIter), TRUE,
					velocity = velocity, c(1L, 1L))
	val = val[c("grid", "velocity")]
	class(val$grid) = z
	val
}

# --------------------------------
#    NOT THE HOMEWORK PROBLEM
# --------------------------------

runGrid = function(dims, num, iter=1000, plot=TRUE) {
	grid = createGrid(dims, num)
	gout = crunBML(grid, iter)
	if (plot) {
		plot(grid)
		plot(gout$grid)
	}
	invisible(list(initial = grid,
		final = gout$grid,
		velocity = gout$velocity))
}


set.seed(13123)

runs = lapply(
  c(.25, .33, .38, .38, .55, .65),
  function(density)
    runGrid(1024, density, 5000, FALSE)
)

plot(runs[[1]]$final, main="rho .25")
plot(runs[[2]]$final, main="rho .33")
plot(runs[[3]]$final, main="rho .38")
plot(runs[[4]]$final, main="rho .38")
plot(runs[[5]]$final, main="rho .55")
plot(runs[[6]]$final, main="rho .65")

# --------------------------------
#     BEGIN HOMEWORK PROBLEM
# --------------------------------

runGrid = function(dims, num, iter=1000, plot=TRUE) {
	grid = createGrid(dims, num)
	g.out = crunBML(grid, iter)
	if (plot) {
		plot(grid)
		plot(g.out$grid)
	}
	tmp = sum(colMeans(g.out$velocity))
	norm = tmp / length(g.out$grid[g.out$grid != "0"])
	c(num, norm)
}


library("parallel")

size = 1024

# Normal execution
run.slow.time = system.time(run.slow <- lapply(
  rep(seq(0.3, 0.4, length.out=20), 8),
  function(density) runGrid(size, density, 64000, FALSE)
))

# Parallel execution
run.fast.time = system.time(run.fast <- mclapply(
  rep(seq(0.3, 0.4, length.out=20), 8),
  function(density) runGrid(size, density, 64000, FALSE),
  mc.cores = 8L
))

plot(matrix(unlist(run.slow), ncol=2, byrow=TRUE), xlab="Density", ylab="Mean % of cars moving", main="Series: L = 64")
plot(matrix(unlist(run.fast), ncol=2, byrow=TRUE), xlab="Density", ylab="Mean % of cars moving", main="Parallel: L = 64")

Makefile

build: clean
	@R CMD SHLIB -o BML.so BML.c RBML.c

clean:
	@rm -f *.o
	@rm -f *.so

run:
	@R --no-save --no-restore < BML.R

RBML.c

#include "BML.h"

/* This is the routine that is called from R and takes the same parameters
   as the move() routine for the most part, but takes also does the looping for
   numSteps of time and specifies the values in slightly different form.
grid - an R _integer_ matrix which is of dimensions r x c (specified in the dims argument)
       and whose values are either 0 for no car, 1 for blue and 2 for red. (Or the other way round for 1 and 2?)
dims - an int array containing 2 elements and giving the number of rows and columns of grid.
redLocations - an integer matrix with as many rows as there are red cars and with 2 columns.
               Each row in the matrix gives the row and column of that red car
numReds - a single integer giving the number of red cars.
blueLocations -  same as redLocations, but for blue cars
numBlues  -  same as numReds, but for blue.
numSteps - an integer indicating how many time steps/iterations to  perform
verbose - a logical value (TRUE or FALSE) indicating whether to print out the iteration
          number or not.
ans - an integer vector or matrix with 2 * numSteps elements into which the number
      of red and then blue cars are inserted.
speeds - an integer vector of length 2 giving the number of cells a red car and a blue
        car respectively moves in a single jump/time interval. This is c(1L, 1L) typically in R.
*/
void
R_BML(int *grid, int *dims, int *redLocations, int *numReds,
      int *blueLocations, int *numBlues, int *numSteps,
      int *verbose, int *ans, int *speeds)
{
	unsigned int t = 0;
	for (t = 0 ; t < *numSteps ; t++) {
		if (*verbose && (t % 100) == 99)
			fprintf(stderr, "%d\n", (int) (t + 1));

		ans[t] = bml_move(grid, dims, redLocations, numReds, RED, speeds[0]);
		ans[t + *numSteps] = bml_move(grid, dims, blueLocations, numBlues, BLUE, speeds[1]);
		if (0 == ans[t] && 0 == ans[t + *numSteps]) {
			fprintf(stderr, "Deatlock at %d\n", t + 1);
			break;
		}
	}
	for (t = t + 1; t < *numSteps ; t++) {
		ans[t] = 0;
		ans[t + *numSteps] = 0;
	}
}