## Boot up lam with lamhosts10 and go into xmpi, or ## boot up bvm with xpvm and 10 hosts ## ## Parallel Kriging ## library(sgeostat) example("krige") ## to get stuff set up # kriging on a grid (slow!) grid <- list(x=seq(min(maas$x),max(maas$x),by=100), y=seq(min(maas$y),max(maas$y),by=100)) grid$xr <- range(grid$x) grid$xs <- grid$xr[2] - grid$xr[1] grid$yr <- range(grid$y) grid$ys <- grid$yr[2] - grid$yr[1] grid$max <- max(grid$xs, grid$ys) grid$xy <- data.frame(cbind(c(matrix(grid$x, length(grid$x), length(grid$y))), c(matrix(grid$y, length(grid$x), length(grid$y), byrow=TRUE)))) colnames(grid$xy) <- c("x", "y") grid$point <- point(grid$xy) data(maas.bank) grid$krige <- krige(grid$point,maas.point,'zinc',maas.vmod, maxdist=1000,extrap=FALSE,border=maas.bank) op <- par(no.readonly = TRUE) par(pty="s") plot(grid$xy, type="n", xlim=c(grid$xr[1], grid$xr[1]+grid$max), ylim=c(grid$yr[1], grid$yr[1]+grid$max)) image(grid$x,grid$y, matrix(grid$krige$zhat,length(grid$x),length(grid$y)), add=TRUE) contour(grid$x,grid$y, matrix(grid$krige$zhat,length(grid$x),length(grid$y)), add=TRUE) data(maas.bank) lines(maas.bank$x,maas.bank$y,col="blue") par(op) parKrige <- function(cl, s, ...) { # split the prediction points s idx <- clusterSplit(cl, 1: dim(s)[1]) ssplt <- lapply(idx, function(i) s[i,]) # compute the predictions in parallel v <- clusterApply(cl, ssplt, krige, ...) # assemble and return the results merge <- function(x, f) do.call("c", lapply(x, f)) s.o <- point(s) s.o$do <- merge(v, function(y) y$do) s.o$zhat <- merge(v, function(y) y$zhat) s.o$sigma2hat <- merge(v, function(y) y$sigma2hat) return(s.o) } system.time(k <- krige(grid$point,maas.point,'zinc',maas.vmod, maxdist=1000,extrap=FALSE,border=maas.bank)) clusterEvalQ(cl, library(sgeostat)) system.time(pk <- parKrige(cl,grid$point,maas.point,'zinc',maas.vmod, maxdist=1000,extrap=FALSE,border=maas.bank)) identical(k, pk) ## Load balanced version parKrigeLB <- function(cl, s, ..., LBF = 4) { # split the prediction points s idx <- splitIndices(dim(s)[1], length(cl) * LBF) ssplt <- lapply(idx, function(i) s[i,]) # compute the predictions in parallel v <- clusterApplyLB(cl, ssplt, krige, ...) # assemble and return the results merge <- function(x, f) do.call("c", lapply(x, f)) s.o <- point(s) s.o$do <- merge(v, function(y) y$do) s.o$zhat <- merge(v, function(y) y$zhat) s.o$sigma2hat <- merge(v, function(y) y$sigma2hat) return(s.o) } system.time(pklb <- parKrigeLB(cl,grid$point,maas.point,'zinc',maas.vmod, maxdist=1000,extrap=FALSE,border=maas.bank)) identical(pklb, pk) ## Making the plot grid$krige <- pk op <- par(no.readonly = TRUE) par(pty="s") plot(grid$xy, type="n", xlim=c(grid$xr[1], grid$xr[1]+grid$max), ylim=c(grid$yr[1], grid$yr[1]+grid$max)) image(grid$x,grid$y, matrix(grid$krige$zhat,length(grid$x),length(grid$y)), add=TRUE) contour(grid$x,grid$y, matrix(grid$krige$zhat,length(grid$x),length(grid$y)), add=TRUE) data(maas.bank) lines(maas.bank$x,maas.bank$y,col="blue") par(op)