# Modelling volume and height development using the Wenk model, Ap - 01.06.2011. Last updated on 20.03.2013.

rm(list = ls())

cppFile <- "/Users/afs80b/Dropbox/Rcpp/"

library(WenkRegression)
help(WenkRegression)

# Constants
c3 <- 0.4
c4 <- 0.0
trans <- 10

options(digits = 14, width = 100)

# Select data 
dataFile <- "Gwy1_486"
tdata <- read.table(paste(cppFile, "StemAnalyses/", dataFile, ".txt", sep = ""), header = TRUE)
names(tdata)

# tdata <- tdata[tdata$age > 29,] # You can select data

# Calculate form factor
tdata$ff <- tdata$v / (pi * (tdata$dbh/200)^2 * tdata$h)

# Calculate observed volume, height, diameter and form factor growth multipliers and the corresponding relative increments
tdata$mv <- calcTrueMultipliers(tdata$v)
tdata$mh <- calcTrueMultipliers(tdata$h)
tdata$md <- calcTrueMultipliers(tdata$d)
tdata$mf <- calcTrueMultipliers(tdata$ff)

# Calculate observed relative volume, height, diameter and form factor increments
tdata$pv <- calcTrueRelIncrements(tdata$v)
tdata$ph <- calcTrueRelIncrements(tdata$h)
tdata$pd <- calcTrueRelIncrements(tdata$d)
tdata$pf <- calcTrueRelIncrements(tdata$ff)

# Calculate observed allometric coefficients (volume -> height and height -> diameter)
tdata$m1 <- log(1 - tdata$pv) / log(1 - tdata$ph)
tdata$m1[1] <- 0
tdata$m2 <- log(1 - tdata$ph) / log(1 - tdata$pd)
tdata$m2[1] <- 0


# Regression c1, c2, m1 and m2 (allometric coefficients) for volume, height and diameter
weight.v <- 1.0 / 3.0
weight.h <- 1.0 / 3.0
weight.d <- 1.0 / 3.0
sum(weight.v, weight.h, weight.d)
abdn0 <- c(0.17, 1.46, 2.56, 2.56)
abdn.L2 <- optim(par = abdn0, fn = calcEnergy, c3 = c3, c4 = c4, trans = 10, tree = tdata, weightV = weight.v, weightH = weight.h, weightD = weight.d, control = list(maxit = 300, trace = TRUE, REPORT = 3000, temp = 0.01, tmax = 100))#, method = "BFGS") 

# Show regression parameters
abdn.L2$par

# If regression is unsuccessful and the results cannot be improved by modifying the starting values use the following routine
abdn.L2$par <- simulatedAnnealing(tdata, c(0.05, 0.1, 0.1, 0.1, 1.0, 5.0, 5.0, 5.0), c(0.05, 1000000, 0, 0.7), c3, c4, trans, weight.v, weight.h, weight.d)

# Calculate statistics
dn <- calcAllRelativeIncrements(tdata$age, tdata$pf, abdn.L2$par[1], abdn.L2$par[2], c3, c4, abdn.L2$par[3], abdn.L2$par[4], trans)

(vvarres <- var(tdata$pv[-1] - dn$pvx[-1], na.rm = TRUE))

(vbias <- mean(tdata$pv[-1] - dn$pvx[-1], na.rm = TRUE))

(vrmse <- sqrt(vvarres + vbias^2))

(rss.v <- sum((tdata$pv[-1] - dn$pvx[-1])^2, na.rm = TRUE))

# (trss.v <- (length(tdata$pv[-1]) - 1) * var(tdata$pv[-1], na.rm = TRUE))
(trss.v <- length(tdata$pv[-1]) * var(tdata$pv[-1], na.rm = TRUE))

(rv2 <- 1 - rss.v / trss.v)

(hvarres <- var(tdata$ph[-1] - dn$phx[-1], na.rm = TRUE))

(hbias <- mean(tdata$ph[-1] - dn$phx[-1], na.rm = TRUE))

(hrmse <- sqrt(hvarres + hbias^2))

(rss.h <- sum((tdata$ph[-1] - dn$phx[-1])^2, na.rm = TRUE))

# (trss.h <- (length(tdata$ph[-1]) - 1) * var(tdata$ph[-1], na.rm = TRUE))
(trss.h <- length(tdata$ph[-1]) * var(tdata$ph[-1], na.rm = TRUE))

(rh2 <- 1 - rss.h / trss.h)

(dvarres <- var(tdata$pd[-1] - dn$pdx1[-1], na.rm = TRUE))

(dbias1 <- mean(tdata$pd[-1] - dn$pdx1[-1], na.rm = TRUE))

(drmse1 <- sqrt(dvarres + dbias1^2))

(dvarres <- var(tdata$pd[-1] - dn$pdx2[-1], na.rm = TRUE))

(rss.d <- sum((tdata$pd[-1] - dn$pdx1[-1])^2, na.rm = TRUE))

# (trss.d <- (length(tdata$pd[-1]) - 1) * var(tdata$pd[-1], na.rm = TRUE))
(trss.d <- length(tdata$pd[-1]) * var(tdata$pd[-1], na.rm = TRUE))

(rd2.1 <- 1 - rss.d / trss.d)

(dbias2 <- mean(tdata$pd[-1] - dn$pdx2[-1], na.rm = TRUE))

(drmse2 <- sqrt(dvarres + dbias2^2))

(rss.d <- sum((tdata$pd[-1] - dn$pdx2[-1])^2, na.rm = TRUE))

# (trss.d <- (length(tdata$pd[-1]) - 1) * var(tdata$pd[-1], na.rm = TRUE))
(trss.d <- length(tdata$pd[-1]) * var(tdata$pd[-1], na.rm = TRUE))

(rd2.2 <- 1 - rss.d / trss.d)
  
# Collect regression results of multiple sessions in output file 
outputF <- paste(cppFile, "StemAnalyses/", "Parameters.txt", sep = "")
if(file.exists(outputF)) {
  parameters <- read.table(outputF, header = TRUE)
  parameters <- rbind(parameters, data.frame(site = tdata$site[1], tree = tdata$tree_no[1], species = tdata$species[1], c1 = abdn.L2$par[1], c2 = abdn.L2$par[2], c3 = c3, c4 = c4, m1 = abdn.L2$par[3], m2 = abdn.L2$par[4], vbias = vbias, vrmse = vrmse, rv2 = rv2, hbias = hbias, hrmse = hrmse, rh2 = rh2, dbias1 = dbias1, drmse1 = drmse1, rd2.1 = rd2.1, dbias2 = dbias2, drmse2 = drmse2, rd2.2 = rd2.2))
  write.table(parameters, row.names = F, file = outputF)
} else {
  parameters <- NULL
  parameters$site <- tdata$site[1]
  parameters$tree <- tdata$tree_no[1]
  parameters$species <- tdata$species[1]
  parameters$c1 <- abdn.L2$par[1]
  parameters$c2 <- abdn.L2$par[2]
  parameters$c3 <- c3
  parameters$c4 <- c4
  parameters$m1 <- abdn.L2$par[3]
  parameters$m2 <- abdn.L2$par[4]
  parameters$vbias <- vbias
  parameters$vrmse <- vrmse
  parameters$rv2 <- rv2
  parameters$hbias <- hbias
  parameters$hrmse <- hrmse
  parameters$rh2 <- rh2
  parameters$dbias1 <- dbias1
  parameters$drmse1 <- drmse1
  parameters$rd2.1 <- rd2.1
  parameters$dbias2 <- dbias2
  parameters$drmse2 <- drmse2
  parameters$rd2.2 <- rd2.2
  write.table(parameters, file = outputF, row.names = F)
}

# Optimise growth parameter c1 for every time step
tdata$c1 <- 0
for (i in 2:length(tdata$age)) {
    abdn.L2.1 <- optimize(calcEnergy1, c(0.05, 2.5), abdn.L2$par[2], c3, c4, abdn.L2$par[3], abdn.L2$par[4], tdata$age[i - 1], tdata$age[i] - tdata$age[i - 1], 10, tdata$pv[i], tdata$ph[i], tdata$pd[i], weight.v, weight.h, weight.d) 
    tdata$c1[i] <- abdn.L2.1$minimum
}
mv1 <- 0
for (i in 2:length(tdata$age)) 
    mv1[i] <- interpolateMultiplier(tdata$age[i - 1], trans, tdata$age[i] - tdata$age[i - 1], as.double(tdata$c1[i]), abdn.L2$par[2], c3, c4) 
v1 <- calcEstimatedGrowthVariable(tdata$v[1], length(tdata$age), mv1) 

par(mar = c(4, 4, 2, 2))
plot(tdata$age, tdata$v, las = 1, ylab = "Volume", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age, v1, col = "red", lwd = 2)

par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], tdata$c1[-1], las = 1, ylim = c(0, 0.6), ylab = "c1", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
curve(abdn.L2$par[1] + x * 0, from = tdata$age[2], to = tdata$age[length(tdata$age)], lty = 1, lwd = 2, col = "red", add = TRUE)

# Calculate estimated volume, height and diameter
dn <- calcEstimatedMultipliers(abdn.L2$par[1], abdn.L2$par[2], c3, c4, abdn.L2$par[3], abdn.L2$par[4], tdata$age, tdata$pf, trans)
tdata <- cbind(tdata, dn)
names(tdata)

tdata$v.est <- calcEstimatedGrowthVariable(tdata$v[1], length(tdata$age), tdata$mv.est) 
tdata$h.est <- calcEstimatedGrowthVariable(tdata$h[1], length(tdata$age), tdata$mh.est) 
tdata$d1.est <- calcEstimatedGrowthVariable(tdata$dbh[1], length(tdata$age), tdata$md1.est) 
tdata$d2.est <- calcEstimatedGrowthVariable(tdata$dbh[1], length(tdata$age), tdata$md1.est) 

# Visualise results
# pdf(file = paste(cppFile, "Gwydyr.pdf", sep = ""))
par(mfrow = c(5, 2))
par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], 1 - tdata$pv[-1], las = 1, ylim = c(0, 1), ylab = "1 - pv", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age[-1], 1 / tdata$mv.est[-1], col = "red", lwd = 2)

par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], 1 - tdata$ph[-1], las = 1, ylim = c(0, 1), ylab = "1 - ph", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age[-1], 1 / tdata$mh.est[-1], col = "red", lwd = 2)

par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], 1 - tdata$pd[-1], las = 1, ylim = c(0, 1), ylab = "1 - pd", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age[-1], 1 / tdata$md1.est[-1], col = "red", lwd = 2)

par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], 1 - tdata$pd[-1], las = 1, ylim = c(0, 1), ylab = "1 - pd", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age[-1], 1 / tdata$md2.est[-1], col = "red", lwd = 2)

# pdf(file = paste(cppFile, "Volume.pdf", sep = ""))
# par(mar = c(3, 3, 2, 2))
par(mar = c(4, 4, 2, 2))
plot(tdata$age, tdata$v, las = 1, ylab = "Volume", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)#, ylim = c(0, 10), xlim = c(10, 70))  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age, tdata$v.est, col = "red", lwd = 2)
# dev.off()

par(mar = c(4, 4, 2, 2))
plot(tdata$age, tdata$h, las = 1, ylab = "Height", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age, tdata$h.est, col = "red", lwd = 2)

par(mar = c(4, 4, 2, 2))
plot(tdata$age, tdata$dbh, las = 1, ylab = "Diameter", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age, tdata$d1.est, col = "red", lwd = 2)

par(mar = c(4, 4, 2, 2))
plot(tdata$age, tdata$dbh, las = 1, ylab = "Diameter", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
lines(tdata$age, tdata$d2.est, col = "red", lwd = 2)
# dev.off()

par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], tdata$m1[-1], las = 1, ylim = c(0, 5), ylab = "m1", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
curve(abdn.L2$par[3] + x * 0, from = tdata$age[2], to = tdata$age[length(tdata$age)], lty = 1, lwd = 2, col = "red", add = TRUE)

par(mar = c(4, 4, 2, 2))
plot(tdata$age[-1], tdata$m2[-1], las = 1, ylim = c(0, 2), ylab = "m2", xlab = "Age", cex = .9, col = "black", pch = 19, axes = FALSE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
curve(abdn.L2$par[4] + x * 0, from = tdata$age[2], to = tdata$age[length(tdata$age)], lty = 1, lwd = 2, col = "red", add = TRUE)


# Visualise relation between c1 and c2 of all trees analysed (collective output file)
xdata <- read.table(paste(cppFile, "StemAnalyses/", "Parameters", ".txt", sep = ""), header = TRUE)
names(xdata)
par(mar = c(4, 4, 2, 2))
library(calibrate)
plot(xdata$c1, xdata$c2, las = 1, ylab = "c2", xlab = "c1", cex = .9, col = "black", pch = 19, axes = FALSE)  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
box(lwd = 2)
textxy(xdata$c1, xdata$c2, xdata$tree, cx = 0.9)
# [-length(xdata$c1)]