# Calculate growth rate using regression of dbh on time in years
growthRegression = function(fulldata, x, graphit= TRUE) 
{
  onetree=subset(fulldata, tag==x & dbh>0)   ## Filter all records of one tag and having a dbh
  if(nrow(onetree)<=1) return(NA)   ## In case there is only one dbh, or none, bail out and return NA
      
  onetree$time = onetree$julian/365.25
  linear_model = lm(formula = onetree$dbh ~ onetree$time)


  if(graphit)   ## Only show the graph it requested by setting graphit=TRUE
   { plot(dbh~time,data=onetree, xlab = 'Julian Time',
          ylab = 'Diameter at breast height/dbh (cm)',
          main = "Relative Growth Rates for Tropical Tree Species",
          pch = 16, col="black")
     abline(linear_model)}
  
  growth_rate = linear_model$coefficients[2]  ## This is where lm returns the slope coefficient
  mean_dbh = mean(onetree$dbh)## Store the mean dbh as well
  initial_dbh = (onetree$dbh[1])
  s = summary(linear_model)
  
  return(c(mean_dbh, initial_dbh, growth_rate, trt, s$r.squared))
}

## Loop through all tree tags to calculate growth rate of each, and save species name, plot, growth, and mean dbh in dataframe
growthAllTree=function(fulldata)
{
  Ntree = nrow(fulldata)
  alltree = unique(fulldata[,c('sp_code','tag','plot')])
  alltree$rate = alltree$meandbh = alltree$initialdbh = alltree$r2 = NA
  
  for(i in 1:Ntree)
      alltree[i,c('meandbh', 'initialdbh','rate', 'r2')]=growthRegression(fulldata,alltree$tag[i],FALSE)

  alltree$logdbh = log(alltree$initialdbh)
  return(data.frame(alltree))
}

#growthAllTree(fulldata)
#growthrate = growthAllTree(fulldata)
#head(growthrate)

Nplot = c(1, 10, 12, 13, 17, 22)
NPplot = c(3, 5, 11, 18, 20, 23)
Pplot = c(4, 6, 9, 15, 16, 21)
Cplot = c(2, 7, 8, 14, 19, 24)
treatments <- cbind(Nplot, NPplot, Pplot, Cplot)

growthSubset = function(gdata=growthrate, onesp = NULL, gmax= 4, gmin=-1 ,centerdata = TRUE) {
  goodgrowth = subset(gdata, rate<=gmax & rate>=gmin)
  
  if(is.null(onesp)) usegrowth = goodgrowth
  else usegrowth=subset(goodgrowth, sp_code==onesp)

  usegrowth$treatment = 'C'
  usegrowth$treatment[usegrowth$plot%in%Nplot] = 'N'
  usegrowth$treatment[usegrowth$plot%in%NPplot] = 'NP'
  usegrowth$treatment[usegrowth$plot%in%Pplot] = 'P'

  usegrowth$ctrdbh = usegrowth$logdbh - mean(usegrowth$logdbh)
  if(centerdata) usegrowth$usedbh=usegrowth$ctrdbh
  else usegrowth$useddbh=usegrowth$logdbh

  gmodel = lm(rate ~ usedbh + treatment, data = usegrowth)
  return(usegrowth) }


growthModel = function(usegrowth)
{gmodel = lm(rate ~ ctrdbh + treatment, data = usegrowth)
return(summary(gmodel))}

graphModel = function(gdata = growthrate, choosesp = NULL, gmax=4, gmin=-2) {
  alldata = growthSubset(gdata, onesp= choosesp, gmax, gmin)
  modall = growthModel(alldata)
  plot(rate ~ usedbh, data=alldata, pch=16, cex=0.5)
  coef = modall$coefficients[,1]
  abline(coef[1], coef[2])
  treatcoef = coef[1] + coef[3:5]
  abline(treatcoef[1], coef[2], col='blue')
  abline(treatcoef[2], coef[2], col='green')
  abline(treatcoef[3], coef[2], col='red')
  return(modall)
}

#alldata= growthSubset(gdata, onesp=choosesp, gmax, gmin)



barGraphmodel = function(gdata = growthrate, choosesp = NULL, gmax=4, gmin=-2, plotnum= FALSE, jig = .1, wid = .1, yr = NULL, title = "title") {
    alldata = growthSubset(gdata, onesp= choosesp, gmax, gmin)
    modall = growthModel(alldata)
    
    gplot = tapply(alldata$rate, list(alldata$plot, alldata$treatment), mean)
    N = nrow(gplot)
    
    coef = modall$coefficients[,1]
    treatcoef = c(coef[1], coef[1] + coef[3:5])
    treatcoef.sd = modall$coefficients[-2,2]
    treatname = colnames(gplot)
  
    
    if(is.null(yr)) yr = range(alldata$rate)
    
    bp = barplot(treatcoef, names.arg = treatname, ylim = yr, border = NA, col = "white", ylab= 'Growth Rate', xlab = 'Treatment', main= title, font.main=4, font.lab=1 )
    arrows(bp, treatcoef - treatcoef.sd, bp, treatcoef + treatcoef.sd, code = 3, angle = 90, length = wid, lwd = 2)
    segments(x0 =0, x1=5, y0=yr[1], y1=yr[1])
    
    clrs = c('black', 'blue', 'green', 'red')
    
    par(xpd=FALSE)
    for(i in 1:4)
    {
      onetreat = subset(alldata, treatment == treatname[i])
      x=rep(bp[i,1]-jig,nrow(onetreat))
    points(x, onetreat$rate, pch =16 , cex= 0.5, col='gray')
    segments(x0=bp[i,1]-2*wid,x1=bp[i,1]+2*wid,y0=treatcoef[i],y1=treatcoef[i],lwd=2)
      }
     for(i in 1:4)
     {
       if(plotnum)
       text(x=rep(bp[i,1] + jig,N), y=gplot[,i], lab=rownames(gplot))
       else(points(x=rep(bp[i,1] + jig,N), y=gplot[,i]))
     }
    #return(gplot)
    
}


#cbind(treatcoef.M, treatcoef.sd, plot.M)

#head(fulldata)