## First load the RUE model

source("RUEmod.r")

## We load the weather data file

cmi <- read.table("cmiWet-1990.txt" , header=TRUE)

## Simulate some data

ans <- RUEmod(cmi$Rad, cmi$Temp, rue = 2.2, lai.c = 0.0152)
#ans <- RUEmod(cmi$Rad, cmi$Temp)

idx <- c(23, 46, 72, 98, 124, 135)
set.seed(12345)
ag <- ans$AG.cum[idx] + rnorm(6,0, sd = 3)
lai <- ans$lai.cum[idx] + rnorm(6, 0 , sd = 2)

agdd <- ans$AGDD[idx] 

## How do we optimize two parameters given the data?

## Need to construct an objective function

## Let's go by parts

## Let's create a rss function

## Parameters rue 2.1998 , lai.c 0.0152

sim <- RUEmod(cmi$Rad, cmi$Temp, rue = 2.2, lai.c = 0.0152) 

sim.ag <- sim$AG.cum[idx]
sim.lai <- sim$lai.cum[idx]

## Comparing graphically

pdf('../figs/obs-sim.pdf')
xyplot(sim.ag ~ ag,
       ylab = "Simulated Above Ground Biomass",
       xlab = "Observed Above Ground Biomass",
       cex = 1.7,
       panel = function(x, y, ...){
         panel.xyplot(x,y, ...)
         panel.abline(0,1)
       })
dev.off()

## What are possible measures of agreement?

## The R-squared is often reported, but I discourage it

fit <- lm(sim.ag ~ ag)
summary(fit)

## R-squared 0.983, adjusted 0.979

## Plotting residuals

resid <- ag - sim.ag

pdf('../figs/residuals.pdf')
xyplot(resid ~ ag,
       ylab = "Residuals",
       xlab = "Observed",
       ylim = c(-6,6),
       cex = 1.7,
       panel = function(x,y,...){
         panel.xyplot(x,y,...)
         panel.abline(h=0)
       })
dev.off()

## The difference between observed and simulated is one of the most basic measures of dis agreement

resid <- ag - sim.ag

## From here the mean Bias can be easily computed

mbias <- mean(resid)

## Why do we have a negative bias? Why over prediction dominates?

## The Mean Squared Error is also just one more step from the bias

mse <- mean(resid^2)

## The root Mean Squared Error

Rmse <- sqrt(mse) ## Same units as the Y variable

## Mean absolute error

mae <- mean(abs(resid))

## Another approach is to make these measures of disagreement in relative terms

## For example the Relative Root Mean Squared Error

mY <- mean(ag)

rRmse <- Rmse/mY

## Relative mean absolute error

rmae <- mean(abs(resid)/abs(ag))

## Small function for model efficiency

mef <- function(obs, sim){

    num <- sum((obs - sim)^2)
    den <- sum((obs - mean(obs))^2)
    ans <- 1 - num/den
    ans
    
}