########################################################################################################
# Title: r-script to compute baseline mortality rate and population trend as a function of size 
# Notes: includes ploting support intervals
# Author: Christian O. Marks (The Nature Conservancy)
# Date: January 21, 2015
########################################################################################################
library(gtools)

# read in the data:
setwd("/Users/cmarks/Documents/Data/demography")
sppnames <- read.csv(file="tblSpecies.csv",header=T)
spplist <- read.csv(file="spplist.csv",header=T)
spptable <- merge(spplist,sppnames,by.x="Species", by.y="Species_code")#includes common species names, maximum measured size, and replication
#read in demographic function parameters:
growth <- read.csv(file="GrowthCoef.csv",header=T)
mort <- read.csv(file="MortCoef.csv",header=T)
density <- read.csv(file="DensCoef.csv",header=T)
beaver <- read.csv(file="BeavMortCoef.csv",header=T)

#set maximum diameters manually for species where support intervals get too big for largest diameters in cm (otherwise keep largest measured diameter):
spptable$maxdiam[spptable$Species == "PODE3"] <- 105
spptable$maxdiam[spptable$Species == "PIST"] <- 100
spptable$maxdiam[spptable$Species == "ACSA3"] <- 100
spptable$maxdiam[spptable$Species == "CAOV2"] <- 80
spptable$maxdiam[spptable$Species == "FRAM2"] <- 100
spptable$maxdiam[spptable$Species == "FRNI"] <- 35
spptable$maxdiam[spptable$Species == "FRPE"] <- 105
spptable$maxdiam[spptable$Species == "JUCI"] <- 60
spptable$maxdiam[spptable$Species == "PRSE2"] <- 60
spptable$maxdiam[spptable$Species == "QURU"] <- 100

###################################################################################################################
# define demographic functions:

neg.exponential.function <- function(diam,a,b) { a*exp(-b*diam) }

neg.exponential.function.prime <- function(diam,a,b) { -b*a*exp(-b*diam) }

neg.power.function <- function(diam,a,b) { a*(diam^(-b)) }

neg.power.function.prime <- function(diam,a,b) { -b*a*(diam^(-b-1)) }

power.function <- function(diam,a,b) { a/10*(diam^b) } # growth function

power.function.prime <- function(diam,a,b) { b*a/10*(diam^(b-1)) }

double.exponential.function <- function(diam,a,b,c,d)
{   1 - ((1-a*exp(b*diam/100)) * (exp(c*((diam/100)^d))))}

weibull.function <- function(diam,a,b,c)
{    (a*(diam)^(c-1))*(exp(-b*(diam)^c))  }

weibull.function.prime <- function(diam,a,b,c)
{    (a*(c-1)*(diam)^(c-2))*(exp(-b*(diam)^c)) + (a*(diam)^(c-1))*(exp(-b*(diam)^c))*(-b*c*(diam)^(c-1))  }

basemort.negexp <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) / (neg.exponential.function(diam,a2,b2)) * (neg.exponential.function.prime(diam,a2,b2)) - (power.function.prime(diam,a1,b1))  } 

basemort.negpow <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) / (neg.power.function(diam,a2,b2)) * (neg.power.function.prime(diam,a2,b2)) - (power.function.prime(diam,a1,b1))  } 

basemort.weibull <- function(diam,a1,b1,a2,b2,c2) 
{    - (power.function(diam,a1,b1)) / (weibull.function(diam,a2,b2,c2)) * (weibull.function.prime(diam,a2,b2,c2)) - (power.function.prime(diam,a1,b1))  } 

delta.term2.negexp <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) * (neg.exponential.function.prime(diam,a2,b2)) - (neg.exponential.function(diam,a2,b2)) * (power.function.prime(diam,a1,b1))  } 

delta.term2.negpow <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) * (neg.power.function.prime(diam,a2,b2)) - (neg.power.function(diam,a2,b2)) * (power.function.prime(diam,a1,b1))  } 

delta.term2.weibull <- function(diam,a1,b1,a2,b2,c2) 
{    - (power.function(diam,a1,b1)) * (weibull.function.prime(diam,a2,b2,c2)) - (weibull.function(diam,a2,b2,c2)) * (power.function.prime(diam,a1,b1))  } 

###################################################################################################################
# prepare functions for computing support intervals:

#function to replace high (1) and low (0) from permutation with actual upper and lower support interval fitted values
valrep <- function(v,l,h){  
v <- ifelse(v==0, l, h)  
return(v) 
}

#function to get all permutations of upper and lower support interval values for fitted coefficients:
permhighlow <- function(r){ # r is the number of coefficients in equation to permute
bv <- permutations(n=2,r=r,v=c(0, 1),repeats.allowed=T)
return(bv)
}

#######################################################################################################################
# make plots with data from all species combined:
#Set plot margins:
par(mar=c(4.5, 4.5, 2.5, 0.5))

# find fitted equation parameters for all species combined:
sp1 <- "All"
maxdiam <- max(spptable$maxdiam)
x <- seq(10/pi,maxdiam,0.25)
results <- data.frame(diam=x,growth=NA,growth.l=NA,growth.u=NA,mort=NA,mort.l=NA,mort.u=NA,basemort=NA,basemort.l=NA,basemort.u=NA,density=NA,density.l=NA,density.u=NA,change=NA,change.l=NA,change.u=NA)
a1 <- subset(growth$a,growth$Species==sp1)
b1 <- subset(growth$b,growth$Species==sp1)
al1 <- subset(growth$al,growth$Species==sp1)
bl1 <- subset(growth$bl,growth$Species==sp1)
au1 <- subset(growth$au,growth$Species==sp1)
bu1 <- subset(growth$bu,growth$Species==sp1)
a2 <- subset(density$a,density$Species==sp1)
b2 <- subset(density$b,density$Species==sp1)
c2 <- subset(density$c,density$Species==sp1)
al2 <- subset(density$al,density$Species==sp1)
bl2 <- subset(density$bl,density$Species==sp1)
cl2 <- subset(density$cl,density$Species==sp1)
au2 <- subset(density$au,density$Species==sp1)
bu2 <- subset(density$bu,density$Species==sp1)
cu2 <- subset(density$cu,density$Species==sp1)
a <- subset(mort$a, mort$Species==sp1)
b <- subset(mort$b, mort$Species==sp1)
c <- subset(mort$c, mort$Species==sp1)
d <- subset(mort$d, mort$Species==sp1)
al <- subset(mort$al, mort$Species==sp1)
bl <- subset(mort$bl, mort$Species==sp1)
cl <- subset(mort$cl, mort$Species==sp1)
dl <- subset(mort$dl, mort$Species==sp1)
au <- subset(mort$au, mort$Species==sp1)
bu <- subset(mort$bu, mort$Species==sp1)
cu <- subset(mort$cu, mort$Species==sp1)
du <- subset(mort$du, mort$Species==sp1)

#observed mortality fitted with double exponential: 
y1 <- double.exponential.function(x,a,b,c,d)
bv <- permhighlow(4)
coefs <- data.frame(a=valrep(bv[,1],al,au),b=valrep(bv[,2],bl,bu),c=valrep(bv[,3],cl,cu),d=valrep(bv[,4],dl,du)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- double.exponential.function(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4])
}
yl1 <- apply(vals,1,min)
yu1 <- apply(vals,1,max)
results$mort <- y1
results$mort.l <- yl1
results$mort.u <- yu1

#baseline mortality calculated from fitted growth and density functions:
y0 <- basemort.weibull(x,a1,b1,a2,b2,c2) 
bv <- permhighlow(5)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2),c2=valrep(bv[,5],cl2,cu2)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- basemort.weibull(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4],coefs[i,5])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results$basemort <- y0
results$basemort.l <- yl0
results$basemort.u <- yu0

# size-dependent growth rate:
y0 <- power.function(x,a1,b1)
bv <- permhighlow(2)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- power.function(x,coefs[i,1],coefs[i,2])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results$growth <- y0
results$growth.l <- yl0
results$growth.u <- yu0

# population size structure:
y0 <- weibull.function(x,a2,b2,c2)
bv <- permhighlow(3)
coefs <- data.frame(a2=valrep(bv[,1],al2,au2),b2=valrep(bv[,2],bl2,bu2),c2=valrep(bv[,3],cl2,cu2)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- weibull.function(x,coefs[i,1],coefs[i,2],coefs[i,3])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results$density <- y0
results$density.l <- yl0
results$density.u <- yu0

# modeled trend in population structure as is :
y0 <- -double.exponential.function(x,a,b,c,d) * weibull.function(x,a2,b2,c2) + delta.term2.weibull(x,a1,b1,a2,b2,c2) 
bv <- permhighlow(9)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2),c2=valrep(bv[,5],cl2,cu2),a=valrep(bv[,6],al,au),b=valrep(bv[,7],bl,bu),c=valrep(bv[,8],cl,cu),d=valrep(bv[,9],dl,du)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- -double.exponential.function(x,coefs[i,6],coefs[i,7],coefs[i,8],coefs[i,9]) * weibull.function(x,coefs[i,3],coefs[i,4],coefs[i,5]) + delta.term2.weibull(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4],coefs[i,5])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results$change <- y0
results$change.l <- yl0
results$change.u <- yu0

save(results, file="AllSpecies.rda")

##############################################################################################################
# make multi-panel plot of demographic rates for data from all species combined (figure 3 in manuscript):

load("AllSpecies.rda")
par(mfcol=c(2,2), oma = c(0.5, 0.5, 0, 0.5))

# plot size-dependent growth rate:
par(mar = c(1,5,4,0))
plot(results$diam,results$growth, type="l",xlab="", xaxt="n",ylim=c(0,0.85), lwd=2,col="blue",lty=1,ylab="DBH growth rate (cm/year)",cex.axis=1.0,cex.lab=1.25)
lines(results$diam,results$growth.l,lwd=2,col="blue",lty=2)
lines(results$diam,results$growth.u,lwd=2,col="blue",lty=2)

# plot size-dependent mortality rates:
par(mar = c(4,5,1,0))
plot(results$diam,results$basemort, type="l",xlab="DBH (cm)",ylim=c(0,0.1), lwd=2,col="red",lty=1,ylab="Probability of mortality (/year)",cex.axis=1.0,cex.lab=1.25)
lines(results$diam,results$mort,lwd=2,col="blue")
lines(results$diam,results$basemort.l,lwd=2,col="red",lty=2)
lines(results$diam,results$basemort.u,lwd=2,col="red",lty=2)
lines(results$diam,results$mort.l,lwd=2,col="blue",lty=2)
lines(results$diam,results$mort.u,lwd=2,col="blue",lty=2)
legend("topright",legend=c("observed","baseline"),lty=c(1,1),lwd=3,col=c("blue","red"),bty="n",cex=1.25)

# plot population size structure:
par(mar = c(1,5,4,0))
plot(results$diam,results$density, type="l",lwd=2,xaxt="n", col="blue",lty=1,xlab="",ylab="Population density (/ha)",cex.axis=1.0,cex.lab=1.25)
lines(results$diam,results$density.l,lwd=2,col="blue",lty=2)
lines(results$diam,results$density.u,lwd=2,col="blue",lty=2)

# plot modeled trend in population structure as is :
par(mar = c(4,5,1,0))
plot(results$diam,results$change, type="l",lwd=2,col="blue",lty=1,xlab="DBH (cm)",ylab="Density change (/ha/year)",cex.axis=1.0,cex.lab=1.25)
lines(results$diam,c(rep(0, length(results$diam))),lty=3,col="black")  
lines(results$diam,results$change.l,lwd=2,col="blue",lty=2)
lines(results$diam,results$change.u,lwd=2,col="blue",lty=2)

savePlot(file="All species multi-panel.png",type="png")

##############################################################################################################
# plot modeled trends for each of the common species individually:

#take data for common species individually:
common.species <- subset(spptable,spptable$N>=50)
sppnum <- dim(common.species)[1]

#Select data and set up functions to fit to the data:

for (j in 1:sppnum) {

# find fitted equation parameters for selected species:
sp1 <- as.character(common.species[j,"Species"])
Dtype <- subset(density$Function,density$Species==sp1)
maxdiam <- subset(spptable$maxdiam, spptable$Species==sp1)
x <- seq(10/pi,maxdiam,0.25)
results <- data.frame(diam=x,growth=NA,growth.l=NA,growth.u=NA,mort=NA,mort.l=NA,mort.u=NA,basemort=NA,basemort.l=NA,basemort.u=NA,density=NA,density.l=NA,density.u=NA,change=NA,change.l=NA,change.u=NA,beaver=NA,beaver.l=NA,beaver.u=NA)
a1 <- subset(growth$a,growth$Species==sp1)
b1 <- subset(growth$b,growth$Species==sp1)
al1 <- subset(growth$al,growth$Species==sp1)
bl1 <- subset(growth$bl,growth$Species==sp1)
au1 <- subset(growth$au,growth$Species==sp1)
bu1 <- subset(growth$bu,growth$Species==sp1)
a2 <- subset(density$a,density$Species==sp1)
b2 <- subset(density$b,density$Species==sp1)
al2 <- subset(density$al,density$Species==sp1)
bl2 <- subset(density$bl,density$Species==sp1)
au2 <- subset(density$au,density$Species==sp1)
bu2 <- subset(density$bu,density$Species==sp1)
a <- subset(mort$a, mort$Species==sp1)
b <- subset(mort$b, mort$Species==sp1)
c <- subset(mort$c, mort$Species==sp1)
d <- subset(mort$d, mort$Species==sp1)
al <- subset(mort$al, mort$Species==sp1)
bl <- subset(mort$bl, mort$Species==sp1)
cl <- subset(mort$cl, mort$Species==sp1)
dl <- subset(mort$dl, mort$Species==sp1)
au <- subset(mort$au, mort$Species==sp1)
bu <- subset(mort$bu, mort$Species==sp1)
cu <- subset(mort$cu, mort$Species==sp1)
du <- subset(mort$du, mort$Species==sp1)
a3 <- subset(beaver$Beaver, beaver$Species==sp1)
au3 <- subset(beaver$Beaverh, beaver$Species==sp1)
al3 <- subset(beaver$Beaverl, beaver$Species==sp1)

#observed mortality fitted with double exponential: 
y1 <- double.exponential.function(x,a,b,c,d)
bv <- permhighlow(4)
coefs <- data.frame(a=valrep(bv[,1],al,au),b=valrep(bv[,2],bl,bu),c=valrep(bv[,3],cl,cu),d=valrep(bv[,4],dl,du)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- double.exponential.function(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4])
}
yl1 <- apply(vals,1,min)
yu1 <- apply(vals,1,max)
results$mort <- y1
results$mort.l <- yl1
results$mort.u <- yu1

#baseline mortality calculated from fitted growth and density functions:
if(Dtype=="exponential") { y0 <-  basemort.negexp(x,a1,b1,a2,b2) } else  {  y0 <-  basemort.negpow(x,a1,b1,a2,b2)  } 
bv <- permhighlow(4)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
if(Dtype=="exponential") { vals[,i] <-  basemort.negexp(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4]) } else  {  vals[,i] <-  basemort.negpow(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4])  }
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results$basemort <- y0
results$basemort.l <- yl0
results$basemort.u <- yu0

#windows()
maxy<-max(max(y1),max(yu1),max(yu0),max(y0))
plot(x,y0, type="l",lwd=2,col="red",ylim=c(0,maxy),xlab="DBH (cm)",ylab="Probability of mortality (/yr)",cex.axis=1.5,cex.lab=1.75)
lines(x,yl1,lwd=2,col="blue",lty=2)
lines(x,yu1,lwd=2,col="blue",lty=2)
lines(x,y1,lwd=2,col="blue")
lines(x,yl0,lwd=2,col="red",lty=2)
lines(x,yu0,lwd=2,col="red",lty=2)
title(main = bquote(italic(.(as.character(common.species[j,"Scientific_name"])))*plain(" N=")*plain(.(as.character(common.species[j,"N"])))))
legend("topright",legend=c("baseline","observed"),lty=c(1,1),lwd=3,col=c("red","blue"),bty="n",cex=1.5)
savePlot(file=paste("Baseline and observed mortality rate for ", as.character(sp1), ".png",sep = ""),type="png")

# plot size-dependent growth rate:
y0 <- power.function(x,a1,b1)
bv <- permhighlow(2)
coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1)) 
vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
for (i in 1:nrow(bv)){
vals[,i] <- power.function(x,coefs[i,1],coefs[i,2])
}
yl0 <- apply(vals,1,min)
yu0 <- apply(vals,1,max)
results$growth <- y0
results$growth.l <- yl0
results$growth.u <- yu0

maxy<-max(max(y0),max(yu0))
plot(x,y0, type="l",lwd=2,col="red",ylim=c(0,maxy),lty=1,xlab="DBH (cm)",ylab="DBH growth rate (cm/yr)",cex.axis=1.5,cex.lab=1.75)
lines(x,yl0,lwd=2,col="red",lty=2)
lines(x,yu0,lwd=2,col="red",lty=2)
title(main = bquote(italic(.(as.character(common.species[j,"Scientific_name"])))*plain(" N=")*plain(.(as.character(common.species[j,"N"])))))
savePlot(file=paste("Growth rate for ", as.character(sp1), ".png",sep = ""),type="png")

# plot modeled trend in population structure as is and with the effect of beaver removed:
if(Dtype=="exponential"){ 
	# calc pop trend with total mortality:
	y0 <- -double.exponential.function(x,a,b,c,d) *  neg.exponential.function(x,a2,b2) + delta.term2.negexp(x,a1,b1,a2,b2)
 	bv <- permhighlow(8)
	coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2),a=valrep(bv[,5],al,au),b=valrep(bv[,6],bl,bu),c=valrep(bv[,7],cl,cu),d=valrep(bv[,8],dl,du)) 
	vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv)))
	for (i in 1:nrow(bv)){
	vals[,i] <- -double.exponential.function(x,coefs[i,5],coefs[i,6],coefs[i,7],coefs[i,8]) *  neg.exponential.function(x,coefs[i,3],coefs[i,4]) + delta.term2.negexp(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4])
	}
	yl0 <- apply(vals,1,min)
	yu0 <- apply(vals,1,max)
	# calc pop trend with effect of beaver subtracted out: 
	y1 <- y0 + (a3) *  neg.exponential.function(x,a2,b2)
 	bv <- permhighlow(9)
	coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2),a=valrep(bv[,5],al,au),b=valrep(bv[,6],bl,bu),c=valrep(bv[,7],cl,cu),d=valrep(bv[,8],dl,du),a3=valrep(bv[,9],al3,au3)) 
	vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv)))
	for (i in 1:nrow(bv)){
	vals[,i] <- -double.exponential.function(x,coefs[i,5],coefs[i,6],coefs[i,7],coefs[i,8]) *  neg.exponential.function(x,coefs[i,3],coefs[i,4]) + delta.term2.negexp(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4]) + (coefs[i,9]) *  neg.exponential.function(x,coefs[i,3],coefs[i,4])
	}
	yl1 <- apply(vals,1,min)
	yu1 <- apply(vals,1,max)
} else { 
	# calc pop trend with total mortality:
	y0 <- -double.exponential.function(x,a,b,c,d) *  neg.power.function(x,a2,b2) + delta.term2.negpow(x,a1,b1,a2,b2)
 	bv <- permhighlow(8)
	coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2),a=valrep(bv[,5],al,au),b=valrep(bv[,6],bl,bu),c=valrep(bv[,7],cl,cu),d=valrep(bv[,8],dl,du)) 
	vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv)))
	for (i in 1:nrow(bv)){
	vals[,i] <- -double.exponential.function(x,coefs[i,5],coefs[i,6],coefs[i,7],coefs[i,8]) *  neg.power.function(x,coefs[i,3],coefs[i,4]) + delta.term2.negpow(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4])
	}
	yl0 <- apply(vals,1,min)
	yu0 <- apply(vals,1,max)
	# calc pop trend with effect of beaver subtracted out:  
	y1 <- y0 + (a3) *  neg.power.function(x,a2,b2)
 	bv <- permhighlow(9)
	coefs <- data.frame(a1=valrep(bv[,1],al1,au1),b1=valrep(bv[,2],bl1,bu1),a2=valrep(bv[,3],al2,au2),b2=valrep(bv[,4],bl2,bu2),a=valrep(bv[,5],al,au),b=valrep(bv[,6],bl,bu),c=valrep(bv[,7],cl,cu),d=valrep(bv[,8],dl,du),a3=valrep(bv[,9],al3,au3)) 
	vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv)))
	for (i in 1:nrow(bv)){
	vals[,i] <- -double.exponential.function(x,coefs[i,5],coefs[i,6],coefs[i,7],coefs[i,8]) *  neg.power.function(x,coefs[i,3],coefs[i,4]) + delta.term2.negpow(x,coefs[i,1],coefs[i,2],coefs[i,3],coefs[i,4]) + (coefs[i,9]) *  neg.power.function(x,coefs[i,3],coefs[i,4])
	}
	yl1 <- apply(vals,1,min)
	yu1 <- apply(vals,1,max)
}
results$change <- y0
results$change.l <- yl0
results$change.u <- yu0
results$beaver <- y1
results$beaver.l <- yl1
results$beaver.u <- yu1

#plot pop trend with total mort
#windows()
maxy <- max(max(y0),0,max(yu0))
miny <- min(min(y0),0,min(yl0))
plot(x,y0, type="l",lwd=2,col="red",ylim=c(miny,maxy),lty=1,xlab="DBH (cm)",ylab="Population density change (/ha/yr)",cex.axis=1.5,cex.lab=1.75)
lines(x,c(rep(0, length(x))),lty=2,col="black") 
lines(x,yl0,lwd=2,col="red",lty=2)
lines(x,yu0,lwd=2,col="red",lty=2)
title(main = bquote(italic(.(as.character(common.species[j,"Scientific_name"])))*plain(" N=")*plain(.(as.character(common.species[j,"N"])))))
savePlot(file=paste("Population density change rate for ", as.character(sp1), ".png",sep = ""),type="png")

#plot pop trend with beaver mort subtracted out
#windows()
maxy <- max(max(y1),max(y0),max(yu1),max(yu0),0)
miny <- min(min(y1),min(y0),min(yl1),min(yl0),0)
plot(x,y0, type="l",lwd=2,col="red",ylim=c(miny,maxy),lty=1,xlab="DBH (cm)",ylab="Population density change (/ha/yr)",cex.axis=1.5,cex.lab=1.75)
lines(x,c(rep(0, length(x))),lty=2,col="black")  
lines(x,y1,lwd=2,col="blue",lty=1)
lines(x,yl1,lwd=2,col="blue",lty=2)
lines(x,yu1,lwd=2,col="blue",lty=2)
lines(x,yl0,lwd=2,col="red",lty=2)
lines(x,yu0,lwd=2,col="red",lty=2)
title(main = bquote(italic(.(as.character(common.species[j,"Scientific_name"])))*plain(" N=")*plain(.(as.character(common.species[j,"N"])))))
legend("bottomright",legend=c("observed","no beaver"),lty=c(1,1),lwd=3,col=c("red","blue"),bty="n",cex=1.5)
savePlot(file=paste("Population density change rate with beaver effect removed for ", as.character(sp1), ".png",sep = ""),type="png")

# plot population size structure:
if(Dtype=="exponential"){ 
	y0 <- neg.exponential.function(x,a2,b2) 
	bv <- permhighlow(2)
	coefs <- data.frame(a2=valrep(bv[,1],al2,au2),b2=valrep(bv[,2],bl2,bu2)) 
	vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
	for (i in 1:nrow(bv)){
	vals[,i] <- neg.exponential.function(x,coefs[i,1],coefs[i,2])
	}
	yl0 <- apply(vals,1,min)
	yu0 <- apply(vals,1,max) 
} else { 
	y0 <- neg.power.function(x,a2,b2) 
	bv <- permhighlow(2)
	coefs <- data.frame(a2=valrep(bv[,1],al2,au2),b2=valrep(bv[,2],bl2,bu2)) 
	vals <- data.frame(matrix( rep(NA, length(x)), nrow=length(x), ncol=nrow(bv))) 
	for (i in 1:nrow(bv)){
	vals[,i] <- neg.power.function(x,coefs[i,1],coefs[i,2])
	}
	yl0 <- apply(vals,1,min)
	yu0 <- apply(vals,1,max)
}
results$density <- y0
results$density.l <- yl0
results$density.u <- yu0

maxy <- max(max(y0),max(yu0))
plot(x,y0, type="l",lwd=2,col="blue",ylim=c(0,maxy),lty=1,xlab="DBH (cm)",ylab="Population density (/ha)",cex.axis=1.5,cex.lab=1.75)
lines(x,yl0,lwd=2,col="blue",lty=2)
lines(x,yu0,lwd=2,col="blue",lty=2)
title(main = bquote(italic(.(as.character(common.species[j,"Scientific_name"])))*plain(" N=")*plain(.(as.character(common.species[j,"N"])))))
savePlot(file=paste("Density for ", as.character(sp1), ".png",sep = ""),type="png")

save(results, file=paste("Results for ", as.character(sp1), ".rda",sep = ""))
}
#########################################################################################################