fig.dim <- 5 knitr::opts_chunk$set(fig.width=2*fig.dim,fig.height=fig.dim,fig.align='center') library(Matrix) library(raster) library(landsim) # helper plotting function pl <- function (z,zlim=range(0,z,finite=TRUE),...) { hab <- do.call(stack,list(habitat)[rep(1,NCOL(z))]); values(hab)[habitable] <- z; plot(hab,nr=1,zlim=zlim,legend.width=2,legend.mar=12,...) }
We'd like to run long-distance migration on a disconnected habitat where much of the range is inaccessible, like this one:
set.seed(42) habitat <- raster::raster(xmn=-1e4, xmx=1e4, ymn=-1e4, ymx=1e4, resolution=100, crs="+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") raster::values(habitat) <- pmax(-10,pmin(20,rcauchy(raster::ncell(habitat)))) habitat <- 20*migrate_raster( habitat, kern="gaussian", sigma=300, radius=1500 ) values(habitat)[values(habitat) < (-5)] <- NA plot(habitat)
We'll use negative values as accessible but inhospitable locations, and NA values as inaccessible locations.
The problem is that the number of nonzero entries in the migration matrix grows quadratically with the radius of the migration kernel. To deal with this, we can use a smaller-radius migration matrix, applied many times. To do this we will
In general, the migration operator we want is [ x \mapsto \sum_{n=0}^m w_n M^n , ] with where $\sum w_n = 1$, and $M$ is some kind of fairly short-distance kernel. Here's the setup:
accessible <- ( !is.na(values(habitat)) ) habitable <- ( accessible & values(habitat) > 0 ) # positions of habitable cells in vector of accessible cells habitable.inds <- cumsum(accessible)[habitable] M <- migration_matrix( habitat, kern="gaussian", sigma=100, radius=400, from=which(accessible), to=which(accessible) ) x <- rpois( sum(habitable), values(habitat)[habitable] ) weighted_migrate <- function (n.weights) { Mnx <- numeric(sum(accessible)) Mnx[habitable.inds] <- x out <- n.weights[1]*Mnx[habitable.inds] for (n in seq_along(n.weights)[-1]) { Mnx <- M %*% Mnx out <- out + n.weights[n]*Mnx[habitable.inds] } return(out) }
With $w_n = (0,0,\ldots,1)$, this should look just like a Gaussian smooth, with standard deviation equal to $\sigma \sqrt{m}$:
n.weights <- c(rep(0,10),1) pl( cbind( x=x, Mx=weighted_migrate(n.weights) ), main=c("original","smoothed") )
If we take the smoothing operator as a geometric sum, [ x \mapsto \frac{p}{1-(1-p)^{n=1}} \sum_{n=0}^m (1-p)^n M^n , ] with $p$ fairly small, then we get a long-tailed distribution:
p <- 1/10 n.weights <- (1-p)^(0:50) n.weights <- n.weights/sum(n.weights) pl( cbind( x=x, Mx=weighted_migrate(n.weights) ), main=c("original","smoothed") )
This is implemented via the n.weights
parameter in migration
objects
(note that the first, "zero migration" index is implied):
p <- 1/10 n.weights <- (1-p)^(0:50) n.weights <- n.weights/sum(n.weights) pop <- population( habitat, accessible = ( !is.na(values(habitat)) ), habitable = ( ( !is.na(values(habitat)) ) & ( values(habitat) > 0 ) ), genotypes = 'aa' ) pop$N[] <- rpois( sum(pop$habitable), values(pop$habitat)[pop$habitable] ) migr <- migration( kern="gaussian", sigma=100, radius=400, n.weights = n.weights[-1] ) migr <- setup_migration( migr, pop ) pl( cbind( x=pop$N[,"aa"], Mx=migrate( pop$N, migr ) ), main=c("original","smoothed") )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.