R/plothist.R

In stranda/holoSimCell: Linked demographic/genetic simulation

#' plot the population history
#'
#' Plots the history of colonization events from a spatially-explicit forward simulation of population expansion.
#'
#' @param ph population history object (direct output from \code{gepophist2.cells})
#' @param maxtime the maximum number of generations to plot, sets the end of the period where colonization history is plotted (i.e., \code{maxtime=100} plots the first 100 generations of the forward simulation).
#'
#' @details 
#' \itemize{
#' This function takes a fine-grained spatio-temporal history of a species generated by holoSimCell's forward-time simulator and produces a map that shows population spread through time. The input object includes the following named components:
#' \item{\code{pophist}} {a data frame with population ID, cell row, cell column, timing of colonization, and source population associated with colonists.}
#' \item{\code{Nvecs}} {a data frame with number of generations columns and number of populations rows that provides population size through time from the forward demographic simulation.}
#' \item{\code{tmat}} {a matrix of pairwise migration rates between populations, used in coalescent simulations.}
#' \item{\code{struct}} {a vector summarizing the structure of the simulated landscape from the forward simulation, x and y dimensions of the grid, number of generations, disperal parameters, etc.}
#' \item{\code{hab_suit}} {a matrix with number of generations rows and number of populations columns that provides habitat suitability through time.}
#' \item{\code{coalhist}} {a data frame specifying historical events for the coalescent simulation, includes time of event (time), source (src) and sink (snk) IDs, and fractional colonization (prop).}
#' \item{\code{popslst}} {a list object with one element per population from the forward-time demographic simulation.  Each element includes information on location on the landscape, colonization time, source populations, and population size through time.}
#' }
#'
#' @return
#' Does not return an object, but creates a multi-panel plot illustrating the history of colonization of cells on the landscape (colored arrows, red arrows are associated with older events, yellow and white with more recent colonization), histogram illustrating the distribution of colonization times, and a simplified plot of centroid (mass-weighted range centroid - 'centroid' in \code{enmSdm::bioticVelocity metrics} - and centroid of the northern portion of the range - 'nCentroid' in \code{enmSdm::bioticVelocity metrics}) biotic velocity through time (generation-by-generation).
#'
#' @examples
#' library(holoSimCell)
#' parms <- drawParms(control = system.file("extdata/ashpaper","Ash_priors.csv",package="holoSimCell"))
#' load(file=paste0(system.file(package="holoSimCell"),"/extdata/landscapes/",pollenPulls[[1]]$file))
#' refpops <- pollenPulls[[1]]$refs
#' avgCellsz <- mean(c(res(landscape$sumrast)))
#'
#' ph = getpophist2.cells(h = landscape$details$ncells, xdim = landscape$details$x.dim, ydim = landscape$details$y.dim,
#'                        landscape=landscape,
#'                        refs=refpops,   
#'                        refsz=parms$ref_Ne,
#'                        lambda=parms$lambda,
#'                        mix=parms$mix,  
#'                        shortscale=parms$shortscale*avgCellsz,  
#'                        shortshape=parms$shortshape, 
#'                        longmean=parms$longmean*avgCellsz,  
#'                        ysz=res(landscape$sumrast)[2], 
#'                        xsz=res(landscape$sumrast)[1], 
#'                        K = parms$Ne) 
#'
#' plothist(ph)
#'
#' @seealso \code{\link{getpophist2.cells}}, \code{\link{plotHistSlice}}, \code{\link[enmSdm]{bioticVelocity}}
#'
#' @export
#'
plothist <- function(ph, maxtime=NULL)
{
    ch=ph$coalhist
    pops=ph$pophist
    pops=data.frame(do.call(rbind,lapply(unique(pops$pop),function(p) 
    {
        if (length(which(pops==p))>1)
        {
            tmp <- pops[pops$pop==p,]
            tmp <- tmp[order(-tmp$arrive),][1,]
        } else pops[pops$pop==p,]
    })))
    
    ch=merge(ch,pops,by.x="src",by.y="pop")[,c("src","time","snk","col","row")]
    ch=ch[with(ch,order(-time,src,snk)),]
    ch=rbind(data.frame(src=pops$pop[(!is.na(pops$arrive))&(pops$arrive==0)],
                        time=pops$arrive[(!is.na(pops$arrive))&(pops$arrive==0)],
                        snk=pops$source[(!is.na(pops$arrive))&(pops$arrive==0)],
                        col=pops$col[(!is.na(pops$arrive))&(pops$arrive==0)],
                        row=pops$row[(!is.na(pops$arrive))&(pops$arrive==0)]),
             ch)
                        
   
    layout(matrix(c(2,1,1,1,
                    2,1,1,1,
                    2,1,1,1,
                    0,3,4,5),nrow=4,ncol=4,byrow=T),
           widths=c(0.22,0.26,0.26,0.26)) 

    
    rows <- floor(min(pops$row)):ceiling(max(pops$row))
    rw <- ceiling(max(rows))
    
    cols <- floor(min(pops$col)):ceiling(max(pops$col))
    cl <- ceiling(max(cols))
    
    ch$coldist=NA

    
    mai = par("mai")
    par(mai=c(0.1,0.1,0.1,0.1))
    
    plot(1,type="n",ylab="",xlab="",ylim=range(rows),xlim=range(cols),axes=F,asp=1)
    rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],col = "gray") 

    ((pops$pop-1) %/% cl) +1->row
    ((pops$pop-1) %% cl) +1 ->col
    
    if ((!is.null(ph$Nvecs))&(min(dim(ph$Nvecs))>0))
    {
        harmmean <- function(x){1/mean(1/x)}
        sz <- apply((ph$Nvecs+1)[pops$pop,],1,harmmean)-1
        pops$sz <- sz
    } else {sz <- log(dim(pops)[1])/6}

    for (i in 1:dim(ch)[1])
    {
        if (!is.na(ch$snk[i]))
            {
                x1=ch[i,"col"]
                y1=ch[i,"row"]
                x0= ch[ch$src==ch[i,"snk"],"col"][1]
                y0= ch[ch$src==ch[i,"snk"],"row"][1]
                if (is.null(maxtime))
                    clrs = heat.colors(max(ch$time,na.rm=T))[ch$time[i]]  else
                    clrs = heat.colors(maxtime)[ch$time[i]]
                
                arrows(x0,y0,x1,y1,
                       col=clrs,
                       lwd=2,
                       length=0.1
                       )
                ch$coldist[i] <- sqrt((y0-y1)^2 + (x0-x1)^2)
            }
    }
    
    points(row~col,pops,pch=16,cex=log(sz+1))

    mxt <- ifelse(is.null(maxtime),max(ch$time,na.rm=T),maxtime)
steps=5
    lgd <- round(seq(1,mxt,length.out=steps))

    par(mai=c(0.5,0.6,0.5,0.6))
    
    barplot(matrix(rep(diff(lgd)[1],steps),ncol=1),col=heat.colors(mxt)[lgd],axes=F,
            main="Event Age \n(from start)")
    for (lvl in 1:steps)
    {
        text(x=0.5,y=lgd[lvl]+(diff(lgd)[1]/2),labels=lgd[lvl],cex=2,adj=0.5)
        }

par(mai=mai)

    hist(ch$coldist,xlab="Colonization distance",main="Realized colonization distances")
    hist(ch$time,xlab="Colonization time",main="Realized colonization times")

    recentph <- ph
    recentph$pophist <- pops 
    recentph$Nvecs <- recentph$Nvecs[,-1]
    pharray <- pophistToArray(recentph, times = seq(-21000,0,by=30))
    bV <- bioticVelocity(pharray$pophistAsArray, metrics = c("centroid","nCentroid"), times = pharray$times, longitude = pharray$longitude, latitude = pharray$latitude, atTimes = seq(-21000,0,by=990))
    cent <- bV$centroidVelocity
    ncent <- bV$nCentroidVelocity
    both <- c(cent,ncent)
    both <- both[!is.na(both)|is.finite(both)]

    plot(cent,type="l",
         xlab="Time (bigger more recent)",ylab="Centroid velocity",
         main="Range movement",ylim=range(c(both)))
    points(ncent,type="l",col="red")
    legend(x=0.4*length(cent),y=0.6*max(both),legend=c("Center","North margin"),col=c(1,2),lty=c(1,1))

    layout(matrix(1))
    
}