bd.R

# Device ask should be true if device is interactive
# options(device.ask.default=grDevices::dev.interactive(orNone = TRUE))

respPlot <- function(pop, b, d, lpos, ylab, plab="Population size", title, logscale=FALSE, legendSize=1, fontSize=1){
	ymin = ifelse(logscale, min(c(b, d)), 0)
	ymax = max(c(b,min(d)))
	logPar <- ifelse(logscale, "y", "")

	plot(pop, b,
		cex.axis = fontSize, cex.lab=fontSize,
		ylim = c(ymin, ymax),
		xlab = plab,
		ylab = ylab,
		type = "l", lwd=2, col="blue", main=title, log=logPar
	)
	lines(pop, d, lty=2, lwd=2)
	if (legendSize>0) legend(lpos, cex=legendSize,
		legend = c("Birth rate", "Death rate"),
		col = c("blue", "black"),
		lty = c(1, 2)
	)
}

## Plot just a rate (difference between birth and death)
ratePlot <- function(
	pop, y, lpos, ylab
	, plab="Population size", title, legendSize=1, fontSize=1
	, xrange=NULL
){
	if (is.null(xrange)) xrange<-range(c(0, pop))
	use <- (pop>xrange[[1]]) & (pop<xrange[[2]])
	yrange <- range(c(0, y[use]))
	plot(pop, y
		, cex.axis = fontSize, cex.lab=fontSize
		, xlab = plab
		, xlim = xrange
		, ylim = yrange
		, ylab = ylab
		, type = "l", lwd=2, col="blue", main=title
	)
	abline(h=0, lty=2, lwd=1.5)
}

# Take populations generated by rfun and plot either per-capita or total plots of birth rates and death rates, possibly doing some complicated stuff to add arrows and elines
bdplots = function(pop, b, d, reportTotal=FALSE, title="Birth-death plot", fontSize, legendSize, plab, elines=FALSE, arrows=TRUE, apos=NULL, awidth=0.3, timeUnit="yr"){
	ylab <- paste0("Per capita rate (1/", timeUnit, ")")
	lpos <- "topright"

	if(reportTotal) {
		ylab <- paste0("Total rate (pop/", timeUnit, ")")
		lpos <- "bottomright"
		b <- b*pop
		d <- d*pop
	}

	if(elines) legendSize <- 0

	par(cex=1.6)
	respPlot(pop, b, d, lpos, ylab, title, fontSize=fontSize,
	legendSize=legendSize, plab=plab)
	# abline(v=40)
	if(elines){
		r <- b-d
		eq <- numeric(0)
		print(eq)
		if(reportTotal)
			eq <- 0
		for(j in 2:(length(pop)-1)){
			if (r[[j]]==0)
				eq <- c(eq, pop[[j]])
			else
				if (r[[j]]*r[[j+1]]<0)
					eq <- c(eq,
					(r[j+1]*pop[j]-r[j]*pop[j+1])/(r[j+1]-r[j]))
		}
		abline(v=eq, lty=3)
		if (arrows){
			af <- approxfun(pop, r)
			if (is.null(apos)) apos <- mean(b)/10
			div <- c(eq, pop[[length(pop)]])
			if (!reportTotal) div <- c(pop[1], div)
			for (j in 1:(length(div)-1)){
				pm <- (div[j]+div[j+1])/2
				pw <- div[j+1]-div[j]
				dir <- sign(af(pm))
				arrows(pm-dir*awidth*pw, apos, pm+dir*awidth*pw, apos)
			}
		}
	}
}
rate <- function(popMax=100, b0=1, bDD=NULL, bAllee=NULL, d0=0.5, dDD=NULL, dAllee=NULL, reportTotal=TRUE, popSteps=100, fontSize=1, legendSize=1, title="", tlab = "Time (years)", plab="Population size", elines=FALSE, arrows=TRUE, xrange=NULL){
	pop <- 1:popSteps*(popMax/popSteps)
	b <- rfun(b0, bDD, bAllee, pop, TRUE)
	d <- rfun(d0, dDD, dAllee, pop, FALSE)

	rplots(pop, b, d, reportTotal, title, fontSize, legendSize, plab, elines, arrows, xrange=xrange)
}

## Plot r or dN/dt; this is a mess and should be integrated with the bd stuff
## Done in Lukang, prior to HYH1
rplots = function(pop, b, d
	, reportTotal=FALSE, title="Rate of change"
	, fontSize, legendSize, plab, elines=FALSE
	, arrows=TRUE, apos=NULL, awidth=0.3, xrange=NULL
){
	ylab <- "Per capita rate (1/t)"
	lpos <- "topright"

	if(reportTotal) {
		ylab <- "Total rate (pop/t)"
		lpos <- "bottomright"
		b <- b*pop
		d <- d*pop
	}
	y = b-d

	if(elines) legendSize <- 0

	par(cex=1.6)
	ratePlot(pop, y, lpos, ylab, title, fontSize=fontSize,
		legendSize=legendSize, plab=plab, xrange=xrange)

	if(elines){
		if(reportTotal)
			eq <- 0
		for(j in 2:(length(pop)-1)){
			if (y[[j]]==0)
				eq <- c(eq, pop[[j]])
			else
				if (y[[j]]*y[[j+1]]<0)
					eq <- c(eq,
					(y[j+1]*pop[j]-y[j]*pop[j+1])/(y[j+1]-y[j]))
		}
		abline(v=eq, lty=3)
		if (arrows){
			af <- approxfun(pop, y)
			if (is.null(apos)) apos <- mean(b)/10
			div <- c(eq, pop[[length(pop)]])
			if (!reportTotal) div <- c(pop[1], div)
			for (j in 1:(length(div)-1)){
				pm <- (div[j]+div[j+1])/2
				pw <- div[j+1]-div[j]
				dir <- sign(af(pm))
				arrows(pm-dir*awidth*pw, apos, pm+dir*awidth*pw, apos)
			}
		}
	}
}

rfun <- function(r0, DD, Allee, pop, birth=TRUE, divOffset=1/2, mmax=1000){
	mult <- 1 + 0*pop
	if (!is.null(DD)) mult <- mult*exp(pop/DD)
	if (!is.null(Allee))
		mult <- mult*exp((Allee+divOffset)/(pop+divOffset))
	mult <- mmax*mult/(mmax+mult)
	if (birth) {mult <- 1/mult}
	return(r0*mult)
}

ndot = function(time, vars, parms){
	ndot <- with(as.list(c(vars, parms)),
		rfun(b0, bDD, bAllee, exp(n), TRUE) 
			- rfun(d0, dDD, dAllee, exp(n), FALSE) 
	)
	list(c(ndot))
}

popSim = function (N0, MaxTime, steps, parms){
	sim <- as.data.frame(lsoda(
		y = c(n=log(N0)),
		times = (0:steps)*MaxTime/steps,
		func = ndot,
		parms
	))
	sim$N <- exp(sim$n)
	return(sim)
}

bd <- function(N0=NULL, MaxTime=20, steps=100, popMax=100, b0=1, bDD=NULL, bAllee=NULL, d0=0.5, dDD=NULL, dAllee=NULL, reportPcTotal="b", popSteps=100, fontSize=1, legendSize=1, title="", tlab = "Time (years)", plab="Population size", elines=FALSE, arrows=TRUE){
	pop <- 1:popSteps*(popMax/popSteps)

	b <- rfun(b0, bDD, bAllee, pop, TRUE)
	d <- rfun(d0, dDD, dAllee, pop, FALSE)

	if (reportPcTotal == "p" | reportPcTotal == "b")
		bdplots(pop, b, d, reportTotal=FALSE, title, fontSize=fontSize, 
			legendSize=legendSize, plab=plab, elines=elines, arrows=arrows)

	if (reportPcTotal == "t" | reportPcTotal == "b")
		bdplots(pop, b, d, reportTotal=TRUE, title, fontSize=fontSize, 
			legendSize=legendSize, plab=plab, elines=elines, arrows=arrows)

	if(!is.null(N0)){
		sim <- bdsim(N0, MaxTime, steps, b0, bDD, bAllee, d0, dDD, dAllee)
		# print(data.frame(sim$time, sim$N))
		plot(sim$time, sim$N,
			cex.lab=fontSize, cex.axis=fontSize,
			main=title, xlab = tlab, ylab = "Population",
			type = "l", ylim = c(1, max(sim$N)), log="y"
		)
		plot(sim$time, sim$N,
			cex.lab=fontSize, cex.axis=fontSize,
			main=title, xlab = tlab, ylab = "Population",
			type = "l", ylim = c(0, max(sim$N))
		)
	}
}

bdsim <- function(N0=1, MaxTime=20, steps=100, b0=1, bDD=NULL, bAllee=NULL, d0=0.5, dDD=NULL, dAllee=NULL){

	parms = list(b0=b0, bDD=bDD, bAllee=bAllee, d0=d0, dDD=dDD, dAllee=dAllee)
	return(popSim(N0, MaxTime, steps, parms))
}

args(bd)