## 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] 
lai <- ans$lai.cum[idx] 

agdd <- ans$AGDD[idx] 

plot1 <- xyplot(ag ~ agdd, main = "Above ground biomass",
                ylab = "Above gorund biomass",
                xlab = "Growing degree days",
                ylim = c(0,70),
                panel = function(x,y,...){
                  panel.xyplot(x,y,col="blue",type="p", cex=1.5)
                  panel.xyplot(ans$AGDD, ans$AG.cum, type="l", col="blue")
                })

plot2 <- xyplot(lai ~ agdd, main = "Leaf Area Index",
                ylab = "Leaf Area Index",
                xlab = "Growing degree days",
                panel = function(x,y,...){
                  panel.xyplot(x,y,col="purple",type="p", cex=1.5)
                  panel.xyplot(ans$AGDD, ans$lai.cum, type="l", col="purple")
                })

pdf('../figs/obs-sim-model.pdf', width = 10)
print(plot1, c(0,0,0.5,1), more = TRUE)
print(plot2, c(0.5,0,1,1))
dev.off()

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

rss <- function(obs, sim){

  if(length(obs) != length(sim))
    stop("length of obs should be equal to length sim")
  
  resid.sq <- (obs - sim)^2

  rss <- sum(resid.sq)

  rss

}

objfun <- function(cfs, obs.ag, obs.lai, sim.idx, wetdat){

  rue <- cfs[1]
  lai.c <- cfs[2]

  sim <- RUEmod(wetdat[,3],wetdat[,4], lai.c = lai.c, rue = rue)

  sim.ag <- sim$AG.cum[sim.idx]
  sim.lai <- sim$lai.cum[sim.idx]
  
  val1 <- rss(obs.ag, sim.ag)

  val2 <- rss(obs.lai, sim.lai)

  val <- val1 + val2

}

op <- optim(c(2,0.01), objfun,
            obs.ag = ag,
            obs.lai = lai,
            sim.idx = idx,
            wetdat = cmi,
            control = list(trace = TRUE))



## Can we recover data simulated with error ?

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

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

agdd <- ans$AGDD[idx] 

op <- optim(c(2,0.01), objfun,
            obs.ag = ag,
            obs.lai = lai,
            sim.idx = idx,
            wetdat = cmi)