03_UncertaintyAnalysis.Rmd

---
title: "03_UncertaintyAnalysis"
author: "Nia Bartolucci; Cameron Reimer; Kangjoon Cho; Zhenpeng Zuo"
date: "4/27/2021"
output: html_document
---

```{r}
#environment load
source('00C_Library+Directory_Setting.R')
```


```{r}
#NOAA data load
# Download NOAA climate forecasts (hourly) and downsample to daily scale
source("00B_NOAAconversion.R")

# define site names
site_names <- c("BART","KONZ","OSBS","SRER")

####If you don't have NOAA data, run this code
#for (S in site_names){
#  download_noaa_files_s3(siteID = S, 
#                         date = "2021-03-01", 
#                         cycle = "00", 
#                         local_directory <- paste0(basePath,"/drivers/"))
#}

NOAA_Driver = noaa_gefs_read(paste0(basePath,"/drivers/noaa/NOAAGEFS_1hr"), "2021-03-01", "00", "KONZ")

predict_time = subset(NOAA_Driver, ensemble==1)
predict_time = predict_time$time

## Driver data preparation

#shortwave radiance
sw_driver = subset(NOAA_Driver, ensemble!=0)
sw_driver = sw_driver$surface_downwelling_shortwave_flux_in_air
sw_driver = matrix(sw_driver, nrow=30 ,byrow = TRUE)

#longwave radiance
lw_driver = subset(NOAA_Driver, ensemble!=0)
lw_driver = lw_driver$surface_downwelling_longwave_flux_in_air
lw_driver = matrix(lw_driver, nrow=30 ,byrow = TRUE)

#air temperature
temp_driver = subset(NOAA_Driver, ensemble!=0)
temp_driver = temp_driver$air_temperature
temp_driver = matrix(temp_driver, nrow=30 ,byrow = TRUE)
tmp = matrix(273.15,30,841)
temp_driver = temp_driver - tmp  # conversion kelvin degree to celcius degree (-273.15)

#precipitation flux
precip_driver = subset(NOAA_Driver, ensemble!=0)
precip_driver = precip_driver$precipitation_flux
precip_driver = matrix(precip_driver, nrow=30 ,byrow = TRUE)
precip_driver = 1800 * precip_driver # unit conversion (30 min -> 1800 sec)

#storage to make 30 min interval driver data
sw_driver_gf = matrix(0, nrow=30, ncol=1679)
lw_driver_gf = matrix(0, nrow=30, ncol=1679)
temp_driver_gf = matrix(0, nrow=30, ncol=1679)
precip_driver_gf = matrix(0, nrow=30, ncol=1679)

## filling gap (interpolation using average)
for(i in 1:839){
  sw_driver_gf[,2*i-1]=sw_driver[,i]
  sw_driver_gf[,2*i]=(sw_driver[,i]+sw_driver[,i+1])/2
  lw_driver_gf[,2*i-1]=lw_driver[,i]
  lw_driver_gf[,2*i]=(lw_driver[,i]+lw_driver[,i+1])/2
  temp_driver_gf[,2*i-1]=temp_driver[,i]
  temp_driver_gf[,2*i]=(temp_driver[,i]+temp_driver[,i+1])/2
  precip_driver_gf[,2*i-1]=precip_driver[,i]
  precip_driver_gf[,2*i]=(precip_driver[,i]+precip_driver[,i+1])/2
}
sw_driver_gf[,1679]=sw_driver[,840]
lw_driver_gf[,1679]=lw_driver[,840]
temp_driver_gf[,1679]=temp_driver[,840]
precip_driver_gf[,1679]=precip_driver[,840]

```


```{r}

loadFilename <- sprintf("%s.Rdata","joint_burn")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

loadFilename <- sprintf("%s.Rdata","Target_time")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

loadFilename <- sprintf("%s.Rdata","Target_LE")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

loadFilename <- sprintf("%s.Rdata","Target_NEE")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

loadFilename <- sprintf("%s.Rdata","Target_VSWC")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

## Forward Simulation
### settings
Nmc = 1000         ## set number of Monte Carlo draws
N.cols <- c("red","green","blue") ## set colors
trans <- 0.8       ## set transparancy
time = 1:length(time_KONZ)+1679    ## total time (1yr + 35 days)
time1 = 1:length(time_KONZ)       ## calibration period
time2 = (length(time_KONZ)+1):(length(time_KONZ)+1679)   ## forecast period
timeN_predict = length(time2)
tmp = matrix(0,1,length(time))
#ylim = c(-500,700)
```


```{r}
##` @param IC    Initial Conditions
##` @param r     Intrinsic growth rate
##` @param Kg    Across-site ('global') mean carrying capacity
##` @param alpha Site random effect
##` @param beta  Slope of precipitation effect on K
##` @param ppt   Precipitation forecast
##` @param Q     Process error (default = 0 for deterministic runs)
##` @param n     Size of Monte Carlo ensemble
forecastN <- function(IC_NEE,IC_LE,IC_VSWC,beta_NEE,beta_LE,beta_VSWC,beta_NL,beta_NV,beta_LV,beta_LN,beta_VN,beta_VL,
                        beta_NEEI,beta_LEI,beta_VSWCI,beta_sw1,beta_sw2,beta_lw,beta_temp,beta_precip,sw,lw,temp,precip
                        ,Q=0,n=Nmc){
  NEE <- matrix(NA,n,timeN_predict)  ## storage
  LE <- matrix(NA,n,timeN_predict)
  VSWC <- matrix(NA,n,timeN_predict)
  Nprev_NEE <- IC_NEE           ## initialize
  Nprev_LE <- IC_LE
  Nprev_VSWC <- IC_VSWC
  
  
  for(t in 1:timeN_predict){
    mu_nee = (1+beta_NEE)*Nprev_NEE+beta_NEEI+beta_NL*Nprev_LE+beta_NV*Nprev_VSWC+beta_sw1*sw[,t] +beta_temp*temp[,t]
    mu_le = (1+beta_LE)*Nprev_LE+beta_LEI+beta_LN*Nprev_NEE+beta_LV*Nprev_VSWC+beta_sw2*sw[,t] +beta_lw*lw[,t]
    mu_vswc = (1+beta_VSWC)*Nprev_VSWC+beta_VSWCI+beta_VN*Nprev_NEE+beta_VL*Nprev_LE+beta_precip*precip[,t]
    
    NEE[,t] <- rnorm(n,mu_nee,Q)                         ## predict next step
    LE[,t] <- rnorm(n,mu_le,Q)
    VSWC[,t] <- rnorm(n,mu_vswc,Q)
    Nprev_NEE <- NEE[,t]                                 ## update IC
    Nprev_LE <- LE[,t]
    Nprev_VSWC <- VSWC[,t]
  }
  N <- list(NEE_predict=NEE, LE_predict=LE, VSWC_predict=VSWC)
  return(N)
}

```

```{r}

## calculate mean of all inputs
sw.mean <- matrix(apply(sw_driver_gf,2,mean),1,timeN_predict) ## driver
lw.mean <- matrix(apply(lw_driver_gf,2,mean),1,timeN_predict)
temp.mean <- matrix(apply(temp_driver_gf,2,mean),1,timeN_predict)
precip.mean <- matrix(apply(precip_driver_gf,2,mean),1,timeN_predict)
## parameters
params_mat <- as.matrix(joint_out$params)
param.mean <- apply(params_mat,2,mean)
## initial conditions
IC <- as.matrix(joint_out$predict)
x.cols <- grep("^NEE",colnames(IC))
IC_NEE <- IC[,x.cols]
x.cols <- grep("^LE",colnames(IC))
IC_LE <- IC[,x.cols]
x.cols <- grep("^VSWC",colnames(IC))
IC_VSWC <- IC[,x.cols]
```

```{r}

##Deterministic prediction
N.det <- forecastN(IC_NEE=mean(IC_NEE[,"NEE[17520]"]),
                   IC_LE=mean(IC_LE[,"LE[17520]"]),
                   IC_VSWC=mean(IC_VSWC[,"VSWC[17520]"]),
                   beta_NEE=param.mean["beta_NEE"],
                   beta_LE=param.mean["beta_LE"],
                   beta_VSWC=param.mean["beta_VSWC"],
                   beta_NL=param.mean["beta_NL"],
                   beta_NV=param.mean["beta_NV"],
                   beta_LV=param.mean["beta_LV"],
                   beta_LN=param.mean["beta_LN"],
                   beta_VN=param.mean["beta_VN"],
                   beta_VL=param.mean["beta_VL"],
                   beta_NEEI=param.mean["beta_NEEI"],
                   beta_LEI=param.mean["beta_LEI"],
                   beta_VSWCI=param.mean["beta_VSWCI"],
                   beta_sw1=param.mean["beta_sw1"],
                   beta_sw2=param.mean["beta_sw2"],
                   beta_lw=param.mean["beta_lw"],
                   beta_temp=param.mean["beta_temp"],
                   beta_precip=param.mean["beta_precip"],
                   sw=sw.mean,
                   lw=lw.mean,
                   temp=temp.mean,
                   precip=precip.mean,
                   Q=0,
                   n=1)

```

## Monte Carlo Error Propagation
## I couldn't figure out why weird sentence appear, anyway we could get result N.I
```{r}
Nmc = 1000
#@ sample parameter rows from previous analysis
prow = sample.int(nrow(params_mat),Nmc,replace=TRUE)

N.I <- forecastN(IC_NEE=IC_NEE[prow,"NEE[17520]"],
                   IC_LE=IC_LE[prow,"LE[17520]"],
                   IC_VSWC=IC_VSWC[prow,"VSWC[17520]"],
                   beta_NEE=param.mean["beta_NEE"],
                   beta_LE=param.mean["beta_LE"],
                   beta_VSWC=param.mean["beta_VSWC"],
                   beta_NL=param.mean["beta_NL"],
                   beta_NV=param.mean["beta_NV"],
                   beta_LV=param.mean["beta_LV"],
                   beta_LN=param.mean["beta_LN"],
                   beta_VN=param.mean["beta_VN"],
                   beta_VL=param.mean["beta_VL"],
                   beta_NEEI=param.mean["beta_NEEI"],
                   beta_LEI=param.mean["beta_LEI"],
                   beta_VSWCI=param.mean["beta_VSWCI"],
                   beta_sw1=param.mean["beta_sw1"],
                   beta_sw2=param.mean["beta_sw2"],
                   beta_lw=param.mean["beta_lw"],
                   beta_temp=param.mean["beta_temp"],
                   beta_precip=param.mean["beta_precip"],
                   sw=sw.mean,
                   lw=lw.mean,
                   temp=temp.mean,
                   precip=precip.mean,
                   Q=0,
                   n=Nmc)
```

## Parameter Uncertainty
## I couldn't figure out why weird sentence appear, anyway we could get result N.IP
```{r}

N.IP <- forecastN(IC_NEE=IC_NEE[prow,"NEE[17520]"],
                   IC_LE=IC_LE[prow,"LE[17520]"],
                   IC_VSWC=IC_VSWC[prow,"VSWC[17520]"],
                   beta_NEE=params_mat[prow,"beta_NEE"],
                   beta_LE=params_mat[prow,"beta_LE"],
                   beta_VSWC=params_mat[prow,"beta_VSWC"],
                   beta_NL=params_mat[prow,"beta_NL"],
                   beta_NV=params_mat[prow,"beta_NV"],
                   beta_LV=params_mat[prow,"beta_LV"],
                   beta_LN=params_mat[prow,"beta_LN"],
                   beta_VN=params_mat[prow,"beta_VN"],
                   beta_VL=params_mat[prow,"beta_VL"],
                   beta_NEEI=params_mat[prow,"beta_NEEI"],
                   beta_LEI=params_mat[prow,"beta_LEI"],
                   beta_VSWCI=params_mat[prow,"beta_VSWCI"],
                   beta_sw1=params_mat[prow,"beta_sw1"],
                   beta_sw2=params_mat[prow,"beta_sw2"],
                   beta_lw=params_mat[prow,"beta_lw"],
                   beta_temp=params_mat[prow,"beta_temp"],
                   beta_precip=params_mat[prow,"beta_precip"],
                   sw=sw.mean,
                   lw=lw.mean,
                   temp=temp.mean,
                   precip=precip.mean,
                   Q=0,
                   n=Nmc)

```

##Driver Uncertainty
```{r}
#sample driver rows
drow_sw = sample.int(nrow(sw_driver_gf),Nmc,replace = TRUE)
drow_lw = sample.int(nrow(lw_driver_gf),Nmc,replace = TRUE)
drow_temp = sample.int(nrow(temp_driver_gf),Nmc,replace = TRUE)
drow_precip = sample.int(nrow(precip_driver_gf),Nmc,replace = TRUE)

N.IPD <- forecastN(IC_NEE=IC_NEE[prow,"NEE[17520]"],
                   IC_LE=IC_LE[prow,"LE[17520]"],
                   IC_VSWC=IC_VSWC[prow,"VSWC[17520]"],
                   beta_NEE=params_mat[prow,"beta_NEE"],
                   beta_LE=params_mat[prow,"beta_LE"],
                   beta_VSWC=params_mat[prow,"beta_VSWC"],
                   beta_NL=params_mat[prow,"beta_NL"],
                   beta_NV=params_mat[prow,"beta_NV"],
                   beta_LV=params_mat[prow,"beta_LV"],
                   beta_LN=params_mat[prow,"beta_LN"],
                   beta_VN=params_mat[prow,"beta_VN"],
                   beta_VL=params_mat[prow,"beta_VL"],
                   beta_NEEI=params_mat[prow,"beta_NEEI"],
                   beta_LEI=params_mat[prow,"beta_LEI"],
                   beta_VSWCI=params_mat[prow,"beta_VSWCI"],
                   beta_sw1=params_mat[prow,"beta_sw1"],
                   beta_sw2=params_mat[prow,"beta_sw2"],
                   beta_lw=params_mat[prow,"beta_lw"],
                   beta_temp=params_mat[prow,"beta_temp"],
                   beta_precip=params_mat[prow,"beta_precip"],
                   sw=sw_driver_gf[drow_sw,],
                   lw=lw_driver_gf[drow_lw,],
                   temp=temp_driver_gf[drow_temp,],
                   precip=precip_driver_gf[drow_precip,],
                   Q=0,
                   n=Nmc)
```

```{r}
## process error samples
#Qmc <- 1/sqrt(params_mat[prow,"Q"])  ## convert from precision to standard deviation

#N.IPDE <- forecastN(IC_NEE=IC_NEE[prow,"NEE[17520]"],
#                   IC_LE=IC_LE[prow,"LE[17520]"],
#                   IC_VSWC=IC_VSWC[prow,"VSWC[17520]"],
#                   beta_NEE=params_mat[prow,"beta_NEE"],
#                   beta_LE=params_mat[prow,"beta_LE"],
#                   beta_VSWC=params_mat[prow,"beta_VSWC"],
#                   beta_NL=params_mat[prow,"beta_NL"],
#                   beta_NV=params_mat[prow,"beta_NV"],
#                   beta_LV=params_mat[prow,"beta_LV"],
#                   beta_LN=params_mat[prow,"beta_LN"],
#                   beta_VN=params_mat[prow,"beta_VN"],
#                   beta_VL=params_mat[prow,"beta_VL"],
#                   beta_NEEI=params_mat[prow,"beta_NEEI"],
#                   beta_LEI=params_mat[prow,"beta_LEI"],
#                   beta_VSWCI=params_mat[prow,"beta_VSWCI"],
#                   beta_sw1=params_mat[prow,"beta_sw1"],
#                   beta_sw2=params_mat[prow,"beta_sw2"],
#                   beta_lw=params_mat[prow,"beta_lw"],
#                   beta_temp=params_mat[prow,"beta_temp"],
#                   beta_precip=params_mat[prow,"beta_precip"],
#                   sw=sw_driver_gf[drow_sw,],
#                   lw=lw_driver_gf[drow_lw,],
#                   temp=temp_driver_gf[drow_temp,],
#                   precip=precip_driver_gf[drow_precip,],
#                   Q=Qmc,
#                   n=Nmc)

## Plot run
#plot.run()
#N.IPDE.ci = apply(N.IPDE,2,quantile,c(0.025,0.5,0.975))
#ecoforecastR::ciEnvelope(time2,N.IPDE.ci[1,],N.IPDE.ci[3,],col=col.alpha(N.cols[4],trans))
#ecoforecastR::ciEnvelope(time2,N.IPD.ci[1,],N.IPD.ci[3,],col=col.alpha(N.cols[3],trans))
#ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
#ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
#lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
## Random effect samples
#tau.mc <- 1/sqrt(params_mat[prow,"tau_site"]) ## convert from precision to std deviation
#aNew.mc <- rnorm(Nmc,0,tau.mc)            ## draw out-of-sample predictions of alpha at a new site

#N.IPDE <- forecastN(IC_NEE=IC_NEE[prow,"NEE[17520]"],
#                   IC_LE=IC_LE[prow,"LE[17520]"],
#                   IC_VSWC=IC_VSWC[prow,"VSWC[17520]"],
#                   beta_NEE=params_mat[prow,"beta_NEE"],
#                   beta_LE=params_mat[prow,"beta_LE"],
#                   beta_VSWC=params_mat[prow,"beta_VSWC"],
#                   beta_NL=params_mat[prow,"beta_NL"],
#                   beta_NV=params_mat[prow,"beta_NV"],
#                   beta_LV=params_mat[prow,"beta_LV"],
#                   beta_LN=params_mat[prow,"beta_LN"],
#                   beta_VN=params_mat[prow,"beta_VN"],
#                   beta_VL=params_mat[prow,"beta_VL"],
#                   beta_NEEI=params_mat[prow,"beta_NEEI"],
#                   beta_LEI=params_mat[prow,"beta_LEI"],
#                   beta_VSWCI=params_mat[prow,"beta_VSWCI"],
#                   beta_sw1=params_mat[prow,"beta_sw1"],
#                   beta_sw2=params_mat[prow,"beta_sw2"],
#                   beta_lw=params_mat[prow,"beta_lw"],
#                   beta_temp=params_mat[prow,"beta_temp"],
#                   beta_precip=params_mat[prow,"beta_precip"],
#                   sw=sw_driver_gf[drow_sw,],
#                   lw=lw_driver_gf[drow_lw,],
#                   temp=temp_driver_gf[drow_temp,],
#                   precip=precip_driver_gf[drow_precip,],
#                   alpha = aNew.mc,
#                   Q=Qmc,
#                   n=Nmc)

## Plot run
#plot.run()
#N.IPDEA.ci = apply(N.IPDEA,2,quantile,c(0.025,0.5,0.975))
#ecoforecastR::ciEnvelope(time2,N.IPDEA.ci[1,],N.IPDEA.ci[3,],col=col.alpha(N.cols[5],trans))
#ecoforecastR::ciEnvelope(time2,N.IPDE.ci[1,],N.IPDE.ci[3,],col=col.alpha(N.cols[4],trans))
#ecoforecastR::ciEnvelope(time2,N.IPD.ci[1,],N.IPD.ci[3,],col=col.alpha(N.cols[3],trans))
#ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
#ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
#lines(time2,N.I.ci[2,],lwd=0.5)
```



```{r}
plot.run_NEE <- function(){
  plot(time,tmp,type='n',ylim=range(NEE_KONZ,na.rm=TRUE),ylab="NEE")
  ecoforecastR::ciEnvelope(time1,ci[1,],ci[3,],col=col.alpha("lightBlue",0.6))
  lines(time1,ci[2,],col="blue")
  points(time1,NEE_BART)
}
```

```{r}
plot.run_LE <- function(){
  plot(time,tmp,type='n',ylim=range(LE_KONZ,na.rm=TRUE),ylab="LE")
  ecoforecastR::ciEnvelope(time1,ci[1,],ci[3,],col=col.alpha("lightBlue",0.6))
  lines(time1,ci[2,],col="blue")
  points(time1,LE_BART)
}
```

```{r}
plot.run_VSWC <- function(){
  plot(time,tmp,type='n',ylim=range(VSWC_KONZ,na.rm=TRUE),ylab="VSWC")
  ecoforecastR::ciEnvelope(time1,ci[1,],ci[3,],col=col.alpha("lightBlue",0.6))
  lines(time1,ci[2,],col="blue")
  points(time1,VSWC_BART)
}
```

### Plot for NEE

```{r,echo=FALSE}
out <- as.matrix(joint_out$predict)
x.cols <- grep("^NEE",colnames(out)) ## grab all columns that start with the letter x
ci <- apply(out[,x.cols],2,quantile,c(0.025,0.5,0.975)) 
plot.run_NEE()
```


```{r}
## Deterministic prediction
plot.run_NEE()
lines(time2,N.det$NEE_predict,col="purple",lwd=3)

```

```{r}
#Monte Carlo plot
plot.run_NEE()
N.I.ci = apply(N.I$NEE_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
## Parameter Uncertainty
plot.run_NEE()
N.IP.ci = apply(N.IP$NEE_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
## Driver Uncertainty
plot.run_NEE()
N.IPD.ci = apply(N.IPD$NEE_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.IPD.ci[1,],N.IPD.ci[3,],col=col.alpha(N.cols[3],trans))
ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
### calculation of variances
varI     <- apply(N.I$NEE_predict,2,var)
varIP    <- apply(N.IP$NEE_predict,2,var)
varIPD   <- apply(N.IPD$NEE_predict,2,var)
#varIPDE  <- apply(N.IPDE,2,var)
#varIPDEA <- apply(N.IPDEA,2,var)
varMat   <- rbind(varI,varIP,varIPD)

## out-of-sample stacked area plot
V.pred.rel <- apply(varMat,2,function(x) {x/max(x)})
plot(time2,V.pred.rel[1,],ylim=c(0,1),type='n',main="Relative Variance: Out-of-Sample",ylab="Proportion of Variance",xlab="time")
ciEnvelope(time2,rep(0,ncol(V.pred.rel)),V.pred.rel[1,],col=N.cols[1])
ciEnvelope(time2,V.pred.rel[1,],V.pred.rel[2,],col=N.cols[2])
ciEnvelope(time2,V.pred.rel[2,],V.pred.rel[3,],col=N.cols[3])
ciEnvelope(time2,V.pred.rel[3,],V.pred.rel[3,],col=N.cols[3])
#ciEnvelope(time2,V.pred.rel[4,],V.pred.rel[5,],col=N.cols[5])
legend("topright",legend=c("Driver","Parameter","InitCond"),col=rev(N.cols),lty=1,lwd=3)

## in-sample stacked area plot
#V.pred.rel.in <- apply(varMat[-5,],2,function(x) {x/max(x)})
#plot(time2,V.pred.rel.in[1,],ylim=c(0,1),type='n',main="Relative Variance: In-Sample",ylab="Proportion of Variance",xlab="time")
#ciEnvelope(time2,rep(0,ncol(V.pred.rel.in)),V.pred.rel.in[1,],col=N.cols[1])
#ciEnvelope(time2,V.pred.rel.in[1,],V.pred.rel.in[2,],col=N.cols[2])
#ciEnvelope(time2,V.pred.rel.in[2,],V.pred.rel.in[3,],col=N.cols[3])
#ciEnvelope(time2,V.pred.rel.in[3,],V.pred.rel.in[4,],col=N.cols[4])
#legend("topright",legend=c("Driver","Parameter","InitCond"),col=rev(N.cols[-5]),lty=1,lwd=5)
```

### Plot for LE

```{r,echo=FALSE}
out <- as.matrix(joint_out$predict)
x.cols <- grep("^LE",colnames(out)) ## grab all columns that start with the letter x
ci <- apply(out[,x.cols],2,quantile,c(0.025,0.5,0.975)) 
plot.run_LE()
```


```{r}
## Deterministic prediction
plot.run_LE()
lines(time2,N.det$LE_predict,col="purple",lwd=3)

```

```{r}
#Monte Carlo plot
plot.run_LE()
N.I.ci = apply(N.I$LE_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
## Parameter Uncertainty
plot.run_LE()
N.IP.ci = apply(N.IP$LE_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
## Driver Uncertainty
plot.run_LE()
N.IPD.ci = apply(N.IPD$LE_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.IPD.ci[1,],N.IPD.ci[3,],col=col.alpha(N.cols[3],trans))
ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
### calculation of variances
varI     <- apply(N.I$LE_predict,2,var)
varIP    <- apply(N.IP$LE_predict,2,var)
varIPD   <- apply(N.IPD$LE_predict,2,var)
#varIPDE  <- apply(N.IPDE,2,var)
#varIPDEA <- apply(N.IPDEA,2,var)
varMat   <- rbind(varI,varIP,varIPD)

## out-of-sample stacked area plot
V.pred.rel <- apply(varMat,2,function(x) {x/max(x)})
plot(time2,V.pred.rel[1,],ylim=c(0,1),type='n',main="Relative Variance: Out-of-Sample",ylab="Proportion of Variance",xlab="time")
ciEnvelope(time2,rep(0,ncol(V.pred.rel)),V.pred.rel[1,],col=N.cols[1])
ciEnvelope(time2,V.pred.rel[1,],V.pred.rel[2,],col=N.cols[2])
ciEnvelope(time2,V.pred.rel[2,],V.pred.rel[3,],col=N.cols[3])
#ciEnvelope(time2,V.pred.rel[3,],V.pred.rel[4,],col=N.cols[4])
#ciEnvelope(time2,V.pred.rel[4,],V.pred.rel[5,],col=N.cols[5])
legend("topright",legend=c("Driver","Parameter","InitCond"),col=rev(N.cols),lty=1,lwd=3)

## in-sample stacked area plot
#V.pred.rel.in <- apply(varMat[-5,],2,function(x) {x/max(x)})
#plot(time2,V.pred.rel.in[1,],ylim=c(0,1),type='n',main="Relative Variance: In-Sample",ylab="Proportion of Variance",xlab="time")
#ciEnvelope(time2,rep(0,ncol(V.pred.rel.in)),V.pred.rel.in[1,],col=N.cols[1])
#ciEnvelope(time2,V.pred.rel.in[1,],V.pred.rel.in[2,],col=N.cols[2])
#ciEnvelope(time2,V.pred.rel.in[2,],V.pred.rel.in[3,],col=N.cols[3])
#ciEnvelope(time2,V.pred.rel.in[3,],V.pred.rel.in[4,],col=N.cols[4])
#legend("topright",legend=c("Driver","Parameter","InitCond"),col=rev(N.cols[-5]),lty=1,lwd=5)
```

### Plot for Soil Moisture
```{r,echo=FALSE}
out <- as.matrix(joint_out$predict)
x.cols <- grep("^VSWC",colnames(out)) ## grab all columns that start with the letter x
ci <- apply(out[,x.cols],2,quantile,c(0.025,0.5,0.975)) 
plot.run_VSWC()
```


```{r}
## Deterministic prediction
plot.run_VSWC()
lines(time2,N.det$VSWC_predict,col="purple",lwd=3)

```

```{r}
#Monte Carlo plot
plot.run_VSWC()
N.I.ci = apply(N.I$VSWC_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```


```{r}
## Parameter Uncertainty
plot.run_VSWC()
N.IP.ci = apply(N.IP$VSWC_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
## Driver Uncertainty
plot.run_VSWC()
N.IPD.ci = apply(N.IPD$VSWC_predict,2,quantile,c(0.025,0.5,0.975),na.rm=TRUE)
ecoforecastR::ciEnvelope(time2,N.IPD.ci[1,],N.IPD.ci[3,],col=col.alpha(N.cols[3],trans))
ecoforecastR::ciEnvelope(time2,N.IP.ci[1,],N.IP.ci[3,],col=col.alpha(N.cols[2],trans))
ecoforecastR::ciEnvelope(time2,N.I.ci[1,],N.I.ci[3,],col=col.alpha(N.cols[1],trans))
lines(time2,N.I.ci[2,],lwd=0.5)
```

```{r}
### calculation of variances
varI     <- apply(N.I$VSWC_predict,2,var)
varIP    <- apply(N.IP$VSWC_predict,2,var)
varIPD   <- apply(N.IPD$VSWC_predict,2,var)
#varIPDE  <- apply(N.IPDE,2,var)
#varIPDEA <- apply(N.IPDEA,2,var)
varMat   <- rbind(varI,varIP,varIPD)

## out-of-sample stacked area plot
V.pred.rel <- apply(varMat,2,function(x) {x/max(x)})
plot(time2,V.pred.rel[1,],ylim=c(0,1),type='n',main="Relative Variance: Out-of-Sample",ylab="Proportion of Variance",xlab="time")
ciEnvelope(time2,rep(0,ncol(V.pred.rel)),V.pred.rel[1,],col=N.cols[1])
ciEnvelope(time2,V.pred.rel[1,],V.pred.rel[2,],col=N.cols[2])
ciEnvelope(time2,V.pred.rel[2,],V.pred.rel[3,],col=N.cols[3])
ciEnvelope(time2,V.pred.rel[3,],V.pred.rel[3,],col=N.cols[3])
#ciEnvelope(time2,V.pred.rel[4,],V.pred.rel[5,],col=N.cols[5])
legend("topright",legend=c("Driver","Parameter","InitCond"),col=rev(N.cols),lty=1,lwd=3)

## in-sample stacked area plot
#V.pred.rel.in <- apply(varMat[-5,],2,function(x) {x/max(x)})
#plot(time2,V.pred.rel.in[1,],ylim=c(0,1),type='n',main="Relative Variance: In-Sample",ylab="Proportion of Variance",xlab="time")
#ciEnvelope(time2,rep(0,ncol(V.pred.rel.in)),V.pred.rel.in[1,],col=N.cols[1])
#ciEnvelope(time2,V.pred.rel.in[1,],V.pred.rel.in[2,],col=N.cols[2])
#ciEnvelope(time2,V.pred.rel.in[2,],V.pred.rel.in[3,],col=N.cols[3])
#ciEnvelope(time2,V.pred.rel.in[3,],V.pred.rel.in[4,],col=N.cols[4])
#legend("topright",legend=c("Driver","Parameter","InitCond"),col=rev(N.cols[-5]),lty=1,lwd=5)
```