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

# read in species 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 parameter values:
growth <- read.csv(file="GrowthCoef.csv",header=T)
mort <- read.csv(file="MortCoef.csv",header=T)
density <- read.csv(file="DensCoef.csv",header=T)

############################################################################
# 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)
}

#########################################################################################################
############## calculate self-thinning relationship between growth and morality (eq. 3 in manuscript):

# function to calculate mortality rate from growth rate assuming all other terms are constant
gm1 <- function(G,a1,b1,a2,b2,c2,a,b,c,d,diam) {
y1 <- double.exponential.function(diam,a,b,c,d) #observed mortality fitted with double exponential
dfdt <- 0 - double.exponential.function(diam,a,b,c,d) *  weibull.function(diam,a2,b2,c2)  + delta.term2.weibull(diam,a1,b1,a2,b2,c2) # pop change rate
delta2.term2.weibull <- function(diam,a1,b1,a2,b2,c2,G) 
{   (G) * (weibull.function.prime(diam,a2,b2,c2)) + (weibull.function(diam,a2,b2,c2)) * (power.function.prime(diam,a1,b1))  } 
M <- -1/(weibull.function(diam,a2,b2,c2)) * (dfdt + delta2.term2.weibull(diam,a1,b1,a2,b2,c2,G))
return(M)
}

# function to calculate mortality rate from growth rate assuming all other terms are constant except that df/dt=0
gm2 <- function(G,a1,b1,a2,b2,c2,a,b,c,d,diam) {
y1 <- double.exponential.function(diam,a,b,c,d) #observed mortality fitted with double exponential
dfdt <- 0 # pop change rate
delta2.term2.weibull <- function(diam,a1,b1,a2,b2,c2,G) 
{   (G) * (weibull.function.prime(diam,a2,b2,c2)) + (weibull.function(diam,a2,b2,c2)) * (power.function.prime(diam,a1,b1))  } 
M <- -1/(weibull.function(diam,a2,b2,c2)) * (dfdt + delta2.term2.weibull(diam,a1,b1,a2,b2,c2,G))
return(M)
}

###############
ming<-0.2
maxg<-1.0
sp1 <- "All"
a1 <- subset(growth$a,growth$Species==sp1)
b1 <- subset(growth$b,growth$Species==sp1)
a2 <- subset(density$a,density$Species==sp1)
b2 <- subset(density$b,density$Species==sp1)
c2 <- subset(density$c,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)
xs <- seq(ming,maxg,0.01)
ys <- gm1(xs,a1,b1,a2,b2,c2,a,b,c,d,diam)
selfthin <- data.frame(x=xs,y=ys)
save(selfthin,file="selfthinning.Rda")
xs <- seq(ming,maxg,0.01)
ys <- gm2(xs,a1,b1,a2,b2,c2,a,b,c,d,diam)
equilthin <- data.frame(x=xs,y=ys)
save(equilthin,file="equilibrium.Rda")

###################################################################################################
############# plot relationship between mortality and growth associated with self-thinning: 

diameter <- 30 #do calcs for modeled 30 cm tree
setwd("/Users/cmarks/Documents/Data/demography/predictions/SI")
#take data for common species:
succession.species <- c("SANI","PODE3","ACNE2","CAOV2","ACSA2","JUCI","PLOC","ULAM","QUPA2","FRPE","ACRU","FRNI","TIAM","CACO15","ABBA", "ACSA3", "BELE", "FAGR", "FRAM2", "PIST", "PRSE2", "QURU", "TSCA", "ULRU")
sppnum<-length(succession.species)
sp.results <- data.frame(matrix( rep(NA, sppnum), nrow=sppnum, ncol=20))
sp.results[,1] <- succession.species
ming<-0.2
maxg<-1.0
#look up mortality and growth rates:
for (i in 1:sppnum) {
	sp1 <- as.character(succession.species[i])
	load(paste("Results for ", as.character(sp1), ".rda",sep = ""))
	closestdiam <- min(abs(results$diam - diameter)) 
	sp.results[i,2:20] <- subset(results,abs(results$diam - diameter)==closestdiam)
}
colnames(sp.results) <- c("Species", colnames(results))

#set positions for labels in plot:
textpos <- data.frame(Species=sp.results$Species,growth=sp.results$growth,mort=sp.results$mort)
textpos[textpos[,1]=="FRAM2",2] <- textpos[textpos[,1]=="FRAM2",2]-0.075
textpos[textpos[,1]=="ACSA3",3] <- textpos[textpos[,1]=="ACSA3",3]-0.005
textpos[textpos[,1]=="ACSA3",2] <- textpos[textpos[,1]=="ACSA3",2]-0.035
textpos[textpos[,1]=="QURU",2] <- textpos[textpos[,1]=="QURU",2]-0.025
textpos[textpos[,1]=="QURU",3] <- textpos[textpos[,1]=="QURU",3]-0.005
textpos[textpos[,1]=="CACO15",2] <- textpos[textpos[,1]=="CACO15",2]-0.025
textpos[textpos[,1]=="CACO15",3] <- textpos[textpos[,1]=="CACO15",3]+0.005
textpos[textpos[,1]=="FRPE",2] <- textpos[textpos[,1]=="FRPE",2]-0.025
textpos[textpos[,1]=="FRPE",3] <- textpos[textpos[,1]=="FRPE",3]+0.005
textpos[textpos[,1]=="BELE",2] <- textpos[textpos[,1]=="BELE",2]-0.025
textpos[textpos[,1]=="BELE",3] <- textpos[textpos[,1]=="BELE",3]+0.005

par(oma = c(0., 0., 0, 0.))
plot(sp.results[,'growth'], sp.results[,'mort'], xlim=c(ming,maxg), pch=16, ylab="Probability of mortality (/year)", xlab="DBH growth rate (cm/year)", cex.axis=1.0,cex.lab=1.25)
arrows(sp.results[,'growth'], sp.results[,'mort.l'], sp.results[,'growth'], sp.results[,'mort.u'], col="grey", length=0.05, angle=90, code=3)
arrows(sp.results[,'growth.l'], sp.results[,'mort'], sp.results[,'growth.u'], sp.results[,'mort'],  col="grey", length=0.05, angle=90, code=3)
points(sp.results[,'growth'], sp.results[,'mort'], pch=16)
text(textpos[,2], textpos[,3], sp.results[,1], cex=0.8, pos=4, col=1) 
#text(0.2, 0.17, paste("r =", as.character(round(cor(sp.results[,2],sp.results[,3]), digits = 2))), cex=1.25, pos=4, col=1)
lines(equilthin$x,equilthin$y,lwd=2, lty=3)
lines(selfthin$x,selfthin$y,lwd=2, lty=3)
savePlot(file="Mortality-growth relationship across species.png",type="png")

########################################################################################################