#+TITLE: dispersal.org #+AUTHOR: Rainer M Krug #+EMAIL: rkrug@ecolmod #+DATE: 2010-11-05 Fri #+DESCRIPTION: #+KEYWORDS: #+LANGUAGE: en #+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t #+OPTIONS: TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc #+INFOJS_OPT: view:nil toc:nil ltoc:t mouse:underline buttons:0 path:http://orgmode.org/org-info.js #+EXPORT_SELECT_TAGS: export #+EXPORT_EXCLUDE_TAGS: noexport #+LINK_UP: #+LINK_HOME: #+XSLT: #+BABEL: :session *R* :results output :tangle dispersal.R * Calculate 2D dispersal kernel Calculates the 2D dispersal kernel based on the 1D dispersal kernel (here: pweibull(q, shape=0.44, scale=9.72, log=FALSE) #+BEGIN_SRC R dispProb2D <- function() { ## This can be changed cumDispProb <- function(q) {pweibull(q, shape=0.44, scale=9.72, log=FALSE)} ## relative strengt of wind in direction wind <- data.frame( dir = (1:360), strength = 1.2 + sin(pi*2*(1:360)/360) ) #### maxDist <- 2000 yres <- 100 xres <- 100 f <- max(abs(wind[[2]])) maxX <- ceiling( f * maxDist / xres ) maxY <- ceiling( f * maxDist / yres ) x <- c(-maxX:maxX) y <- c(-maxY:maxY) dm <- matrix(y, nrow=length(x), ncol=length(y), byrow=TRUE) dir <- round( atan2(dm, -y) * 180 / pi ) dir[dir<0] <- (360 + dir[dir<0]) f <- apply( dir, c(1,2), function(d) {ifelse( d==0, 1, wind$strength[[d]])} ) colnames(f) <- y rownames(f) <- x dm <- sqrt((dm*yres)^2 + (x*xres)^2) colnames(dm) <- y rownames(dm) <- x dm <- dm / f vm.min <- cumDispProb(dm-(min(xres, yres)/2)) vm.max <- cumDispProb(dm+(min(xres, yres)/2)) vm <- vm.max - vm.min vt <- vm / sum(vm) return(vt) } #+END_SRC #+results: * Dispersal in 2D Dispersal of seeds in seed, using sd2D ** To be optimized *** dispRorg This is the original function #+BEGIN_SRC R :tangle no dispRorg <- function( maxx = maxx, maxy = maxy, seeds = seeds, sd2D = sd2D, dx = dx, dy = dy, dispSeeds = dispSeeds # the return value ){ for (x in 1:maxx) { cat(x, " ") for (y in 1:maxy) { s <- seeds[x:(x+2*dx), y:(y+2*dy)][!is.na(seeds[x:(x+2*dx), y:(y+2*dy)])] dispSeeds[x, y] <- ifelse( length(s) > 0, sum( rbinom( seeds[x:(x+2*dx), y:(y+2*dy)], seeds[x:(x+2*dx), y:(y+2*dy)], sd2D ), na.rm=TRUE ), NA ) } } cat("\n") return(dispSeeds) } #+END_SRC #+results: *** dispRnew #+BEGIN_SRC R :tangle no dispRnew <- function( maxx = maxx, maxy = maxy, seeds = seeds, sd2D = sd2D, dx = dx, dy = dy, dispSeeds = dispSeeds # the return value ){ ## if ("package:jit" %in% search()) { jit(2,1) } for (x in 1:maxx) { cat(x, " ") for (y in 1:maxy) { s <- seeds[x:(x+2*dx), y:(y+2*dy)] dispSeeds[x, y] <- ifelse( all(is.na(s)), NA, sum( rbinom( s, s, sd2D ), na.rm=TRUE ) ) } } ## if ("package:jit" %in% search()) { jit(0) } cat("\n") return(dispSeeds) } #+END_SRC *** dispRpar #+BEGIN_SRC R :tangle no dispRpar <- function( maxx = maxx, maxy = maxy, seeds = seeds, sd2D = sd2D, dx = dx, dy = dy, dispSeeds = dispSeeds # the return value ) { # jit(1,1) dx2 <- 2*dx dy2 <- 2*dy return( sapply( 1:maxy, function(y){ cat(y, " ") sapply( 1:maxx, function(x){ s <- seeds[x:(x+dx2), y:(y+dy2)] if (all(is.na(s))) { return(NA) } else { return( sum( rbinom( s, s, sd2D ), na.rm = TRUE ) ) } } ) } ) ) } #+END_SRC *** dispRparFor #+BEGIN_SRC R dispRparFor <- function( maxx = maxx, maxy = maxy, seeds = seeds, sd2D = sd2D, dx = dx, dy = dy, dispSeeds = dispSeeds # the return value ) { dx2 <- 2*dx dy2 <- 2*dy for ( y in 1:maxy ) { cat(y, " ") for (x in 1:maxx) { s <- seeds[x:(x+dx2), y:(y+dy2)] if (all(is.na(s))) { dispSeeds[x,y] <- NA } else { dispSeeds[x,y] <- sum( rbinom( s, s, sd2D ), na.rm = TRUE ) } } } return(dispSeeds) } #+END_SRC *** dispRparForSugarOrg #+BEGIN_SRC R fx <- cxxfunction( sig = signature(DX2 = "integer", DY2 = "integer", SD2D = "numeric", SEEDS = "integer", DISPSEEDS = "integer"), body = ' // The input parameter Rcpp::IntegerVector dx2 (DX2); Rcpp::IntegerVector dy2 (DY2); Rcpp::NumericVector sd2D (SD2D); Rcpp::IntegerMatrix seeds (SEEDS); Rcpp::IntegerMatrix dispSeeds (DISPSEEDS); // internal variables Rcpp::NumericVector s (SD2D); // int res; int indS; //counter for matrix seeds // for( int y=0; y0) { res += rbinom( 1, s[i], sd2D[i] )[0]; } } dispSeeds(x, y) = res; } } return wrap( dispSeeds ); ', plugin = "Rcpp", verbose = TRUE ) #+END_SRC #+BEGIN_SRC R dispRparForSugarOrg <- function( maxx = maxx, maxy = maxy, seeds = seeds, sd2D = sd2D, dx = dx, dy = dy, dispSeeds = dispSeeds # the return value ) { library(inline) library(Rcpp) dx2 <- 2*dx dy2 <- 2*dy return(fx(dx2, dy2, sd2D, seeds, dispSeeds)) } #+END_SRC #+results: * The main dispersal function #+BEGIN_SRC R disperse <- function(seeds, sd2D, dispFun=dispRorg) { dispSeeds <- seeds * 0 # matrix to hold the dispersed seeds dx <- (ncol(sd2D) - 1) / 2 dy <- (nrow(sd2D) - 1) / 2 maxx <- nrow(dispSeeds) maxy <- ncol(dispSeeds) # return( dispFun(maxx, maxy, seeds, sd2D, dx, dy, dispSeeds) ) ## buffer for dispersal into cells at the edge buffer <- matrix(NA, nrow=nrow(seeds), ncol=dx) seeds <- cbind(buffer, seeds, buffer) buffer <- matrix(NA, ncol=ncol(seeds), nrow=dy) seeds <- rbind(buffer, seeds, buffer) ## Disperse seeds from seeds.m into dispSeeds ## use rbinom to calculate the number of seeds dispersed INTO the cell dispSeeds[x,y] maxx <- nrow(dispSeeds) maxy <- ncol(dispSeeds) dispSeeds <- dispFun(maxx, maxy, seeds, sd2D, dx, dy, dispSeeds) return(dispSeeds) } #+END_SRC #+results: * Preparation Preparation of tests #+BEGIN_SRC R # seeds.m <- as.matrix(read.table("dispSeeds.m.txt")) # seeds to be dispersed set.seed <- 123 seeds.m <- matrix(ceiling(runif(100*100)*5000), ncol=200) x[x<2] <- NA sd2D <- dispProb2D() # 2D dispersal kernel #+END_SRC * For Illustration #+BEGIN_SRC R ## to execute: print(range(sd2D)) print(system.time(dispSeeds <- disperse(seeds.m, sd2D, dispRparForSugarOrg))) print(range(sd2D)) ## par(mfcol=c(2,2)) ## image(sd2D) ## image(seeds.m) ## dispSeeds[dispSeeds==0] <- NA ## image(dispSeeds) #+END_SRC