04_IterativeParticleFilter.Rmd

---
title: "04_IterativeParticleFilter"
author: "Nia Bartolucci"
date: "4/26/2021"
output: html_document
---

Potential issues with this iterative forecast:
1. Initial conditions/model parameters sampled independently from the chain rather than jointly. I think it's better to sample them jointly (with replacement) because they could be correlated under the chain.
2. Standard deviation for the particle filter (Chunk 9) is currently just 1. This should be obtained by the data as suggested by Kangjoon.
3. Due to issues with the numbers, some questionable math is used. We just normalized the weights so they all lie between 10^-20, 10^20. Maybe one should resample instead.
4. This is not really a problem with the forecast, but rather with the model. The model thinks the NEE is periodic with a similar amplitude year round, which makes it very inaccurate in the winter.


```{r}
## Package check and load

source('00C_Library+Directory_Setting.R')

```

```{r}

# load the data file [30 min Target data]
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)
loadFilename <- sprintf("%s.Rdata","Target_time")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

```

```{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}
# load MCMC output
loadFilename <- sprintf("%s.Rdata","joint_burn")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)


#SET ENSEMBLE RUNS
ne = 500         #needs to stay 30 unless we also sample (with replacement) the noaa driver ensembles
#yes, we will want way more than 30, so will have to sample with replacement

re_sw_driver_gf=sw_driver_gf[sample(nrow(sw_driver_gf),size=ne,replace=TRUE),]
re_lw_driver_gf=lw_driver_gf[sample(nrow(sw_driver_gf),size=ne,replace=TRUE),]
re_temp_driver_gf=temp_driver_gf[sample(nrow(sw_driver_gf),size=ne,replace=TRUE),]
re_precip_driver_gf=precip_driver_gf[sample(nrow(sw_driver_gf),size=ne,replace=TRUE),]

#AVERAGE CHAINS
params = joint_out$params
params_re = matrix(NA,nrow=1503,ncol=23)
params_re[1:501,] = params[[1]]
params_re[502:1002,] = params[[2]]
params_re[1003:1503,] = params[[3]]

#parsing MCMC output 
beta_LE = sample(params_re[,1], ne, replace=TRUE)
beta_LEI = sample(params_re[,2], ne, replace = TRUE)
beta_LN = sample(params_re[,3], ne, replace = TRUE)
beta_LV = sample(params_re[,4], ne, replace = TRUE)
beta_NEE = sample(params_re[,5], ne, replace = TRUE)
beta_NEEI = sample(params_re[,6], ne, replace = TRUE)
beta_NL = sample(params_re[,7], ne, replace = TRUE)
beta_NV = sample(params_re[,8], ne, replace = TRUE)
beta_VL = sample(params_re[,9], ne, replace = TRUE)
beta_VN = sample(params_re[,10], ne, replace = TRUE)
beta_VSWC = sample(params_re[,11], ne, replace = TRUE)
beta_VSWCI = sample(params_re[,12], ne, replace = TRUE)
beta_lw = sample(params_re[,13], ne, replace = TRUE)
beta_precip = sample(params_re[,14], ne, replace = TRUE)
beta_sw1 = sample(params_re[,15], ne, replace = TRUE)
beta_sw2 = sample(params_re[,16], ne, replace = TRUE)
beta_temp = sample(params_re[,17], ne, replace = TRUE)
tau_le_add = sample(params_re[,18], ne, replace = TRUE)
tau_le_obs = sample(params_re[,19], ne, replace = TRUE)
tau_nee_add = sample(params_re[,20], ne, replace = TRUE)
tau_nee_obs = sample(params_re[,21], ne, replace = TRUE)
tau_vswc_add = sample(params_re[,22], ne, replace = TRUE)
tau_vswc_obs = sample(params_re[,23], ne, replace = TRUE)
#why sample from chain 2? 

#I feel that sampling from each of these independently is not correct, rather we should sample 30 random
#iterates of the 500*3 different iterates from the chains. 
#For now we can leave it this way.


#Initial conditions: starting from last observed value <-- is this a bad idea? -- yes :P
#qa_nee = joint_out$data$NEE_obs[!is.na(joint_out$data$NEE_obs)]
#IC_NEE = rnorm(ne, mean = qa_nee[length(qa_nee)], sd = 0.1)
#rm(qa_nee)


#qa_le = joint_out$data$LE_obs[!is.na(joint_out$data$LE_obs)]
#IC_LE = rnorm(ne, mean = qa_le[length(qa_le)], sd = 0.1)
#rm(qa_le)

#qa_vswc = joint_out$data$VSWC_obs[!is.na(joint_out$data$VSWC_obs)]   #remember outlier 
#IC_VSWC = rnorm(ne, mean = qa_vswc[length(qa_vswc)], sd = 0.1)
#rm(qa_vswc)

#the above seems to start at some arbitrary value, and add normal noise. 
#Instead we should sample from the chain output. 
#also arbitrarily sample form the second chain; could improve by sampling from all three.
#17520, 35040, 52560 are the last predictions for LE, NEE, VSWC repectively.

predict = joint_out$predict
predict_re = matrix(NA,nrow=1503,ncol=52560)
predict_re[1:501,] = predict[[1]]
predict_re[502:1002,] = predict[[2]]
predict_re[1003:1503,] = predict[[3]]

IC_LE = sample(predict_re[,17520], ne, replace = TRUE)
IC_NEE = sample(predict_re[,35040], ne, replace = TRUE)
IC_VSWC = sample(predict_re[,52560], ne, replace = TRUE)


```

```{r}
ensembleforecast <- 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){
  
  
  Nprev_NEE <- IC_NEE           
  Nprev_LE <- IC_LE
  Nprev_VSWC <- IC_VSWC
  NEE = (1+beta_NEE)*Nprev_NEE+beta_NEEI+beta_NL*Nprev_LE+beta_NV*Nprev_VSWC+beta_sw1*sw +beta_temp*temp
  LE = (1+beta_LE)*Nprev_LE+beta_LEI+beta_LN*Nprev_NEE+beta_LV*Nprev_VSWC+beta_sw2*sw +beta_lw*lw
  VSWC = (1+beta_VSWC)*Nprev_VSWC+beta_VSWCI+beta_VN*Nprev_NEE+beta_VL*Nprev_LE+beta_precip*precip
  return(cbind(NEE=NEE, LE=LE, VSWC=VSWC))
                     }
```




```{r}
#Initial Forecast

nt = 1679
#nt = 35 * 48                           ## 35 days of 30min; production run should be nrow(inputs) *********
output = array(NA, c(ne, nt, 3))     ## output storage [time step,ensembles,variables]

## forward ensemble simulation
for(t in 1:nt){
  output[,t , ] <- ensembleforecast(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,
                     re_sw_driver_gf[,t],re_lw_driver_gf[,t],re_temp_driver_gf[,t],re_precip_driver_gf[,t])  
  #reset initial conditions
  IC_NEE = output[,t ,1]
  IC_LE = output[,t ,2]
  IC_VSWC = output[,t ,3]
  #X <- output[t, , 1:3]                          ## set most recent prediction to be the next IC
  #if((t %% 336) == 0) print(t / 336)             ## counter: weeks elapsed (7*48 = 1 week)
}
```


```{r}
## Forward Simulation
### settings
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))
out <- as.matrix(joint_out$predict)

time_plot = time_KONZ
for (i in 1:1679){
  time_plot[17520+i]<-time_plot[17519+i]+1800
} 
time_predict = time_plot[17521:19199]
#ylim = c(-500,700)
```


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

```{r}
plot.run_NEE2 <- function(){
  plot(time_plot,tmp,type='n',ylim=range(NEE_KONZ,na.rm=TRUE),ylab="NEE",main = "NEE Ensemble Forecast", 
       xlab = "time", xlim=c(as.POSIXct("2021-03-10",tz="UTC"),as.POSIXct("2021-03-20",tz="UTC")))
  ecoforecastR::ciEnvelope(time_KONZ,ci[1,],ci[3,],col=col.alpha("lightBlue",0.6))
  lines(time_KONZ,ci[2,],col="blue")
  points(time_KONZ,NEE_BART)
}
```

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

```{r}
plot.run_LE2 <- function(){
  plot(time_plot,tmp,type='n',ylim=range(LE_KONZ,na.rm=TRUE),ylab="LE",main = "LE Ensemble Forecast", 
       xlab = "time", xlim=c(as.POSIXct("2021-03-10",tz="UTC"),as.POSIXct("2021-03-20",tz="UTC")))
  ecoforecastR::ciEnvelope(time_KONZ,ci[1,],ci[3,],col=col.alpha("lightBlue",0.6))
  lines(time_KONZ,ci[2,],col="blue")
  points(time_KONZ,NEE_BART)
}
```

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

```{r}
plot.run_VSWC2 <- function(){
  plot(time_plot,tmp,type='n',ylim=range(ci,na.rm=TRUE),ylab="VSWC",main = "VSWC Ensemble Forecast", 
       xlab = "time", xlim=c(as.POSIXct("2021-03-10",tz="UTC"),as.POSIXct("2021-03-20",tz="UTC")))
  ecoforecastR::ciEnvelope(time_KONZ,ci[1,],ci[3,],col=col.alpha("lightBlue",0.6))
  lines(time_KONZ,ci[2,],col="blue")
  points(time_KONZ,VSWC_BART)
}
```

### Plot for NEE

```{r,echo=FALSE}
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()
ci = apply(output[, , 1], 2, quantile, c(0.025, 0.5, 0.975),na.rm=TRUE)   ## calculate CI over ensemble members
ciEnvelope(time_predict, ci[1, ], ci[3, ], col = col.alpha(N.cols[1], 0.5)) ## plot interval
lines(time_predict,ci[2,],lwd=0.5)
```

### Plot for NEE (subset)

```{r,echo=FALSE}
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_NEE2()
ci = apply(output[, , 1], 2, quantile, c(0.025, 0.5, 0.975),na.rm=TRUE)   ## calculate CI over ensemble members
ciEnvelope(time_predict, ci[1, ], ci[3, ], col = col.alpha(N.cols[1], 0.5)) ## plot interval
lines(time_predict,ci[2,],lwd=0.5)
```

### Plot for NEE

```{r,echo=FALSE}
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()
ci = apply(output[, , 2], 2, quantile, c(0.025, 0.5, 0.975),na.rm=TRUE)   ## calculate CI over ensemble members
ciEnvelope(time_predict, ci[1, ], ci[3, ], col = col.alpha(N.cols[2], 0.5)) ## plot interval
lines(time_predict,ci[2,],lwd=0.5)
```

### Plot for LE (subset)

```{r,echo=FALSE}
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_LE2()
ci = apply(output[, , 2], 2, quantile, c(0.025, 0.5, 0.975),na.rm=TRUE)   ## calculate CI over ensemble members
ciEnvelope(time_predict, ci[1, ], ci[3, ], col = col.alpha(N.cols[2], 0.5)) ## plot interval
lines(time_predict,ci[2,],lwd=0.5)
```

### Plot for VSWC

```{r,echo=FALSE}
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()
ci = apply(output[, , 3], 2, quantile, c(0.025, 0.5, 0.975),na.rm=TRUE)   ## calculate CI over ensemble members
ciEnvelope(time_predict, ci[1, ], ci[3, ], col = col.alpha(N.cols[3], 0.5)) ## plot interval
lines(time_predict,ci[2,],lwd=0.5)
```

### Plot for VSWC (subset)

```{r,echo=FALSE}
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_VSWC2()
ci = apply(output[, , 3], 2, quantile, c(0.025, 0.5, 0.975),na.rm=TRUE)   ## calculate CI over ensemble members
ciEnvelope(time_predict, ci[1, ], ci[3, ], col = col.alpha(N.cols[3], 0.5)) ## plot interval
lines(time_predict,ci[2,],lwd=0.5)
```



```{r}
#get new data into the right format
loadFilename <- sprintf("%s.Rdata","Dataframe_updated")
loadFilename <- paste(dataPath, loadFilename, sep="", collapse = NULL)
load(file = loadFilename)

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

NEEnew = data_up$NEE_obs
LEnew = data_up$LE_obs
VSWCnew = data_up$VSWC_obs

```



```{r}
## calculate the cumulative likelihoods
## to be used as PF weights

NEEcast = output[,,1]
LEcast = output[,,2]
VSWCcast = output[,,3]
NEElike = array(NA, dim(NEEcast))  ## storage, same size as model [ensemble, NEE new time points]
LElike = array(NA, dim(LEcast))
VSWClike = array(NA, dim(VSWCcast))
#NEEcalc = matrix(NA, nrow = nrow(NEEcast), ncol =length(NEEnew))
#print(dim(NEEcalc))
sel = 1:length(NEEnew)
for(i in 1:ne){
  #print(length(NEEcast[i, sel]))
  #print(length(NEEnew))
  NEElike[i, sel] = dnorm(NEEcast[i, sel], NEEnew, 1, log = TRUE)  ## calculate log likelihoods 
  LElike[i,sel] = dnorm(LEcast[i,sel],LEnew,1,log = TRUE)
  VSWClike[i,sel] = dnorm(VSWCcast[i,sel],VSWCnew,1,log=TRUE)
  
  #for now just using standard dev 1 from lack of standard dev data!! Need to fix this.
  
  NEElike[i, is.na(NEElike[i, ])] = 0       ## missing data as weight 1; log(1)=0
  LElike[i, is.na(LElike[i, ])] = 0 
  VSWClike[i, is.na(VSWClike[i, ])] = 0 
  #print(NEElike)
  #NEElike[i, ] = exp(cumsum(NEElike[i, ]))  ## convert to cumulative log likelihood and take out of log-space
  NEElike[i, ] = cumsum(NEElike[i, ])
  LElike[i,] = cumsum(LElike[i,])
  VSWClike[i,]=cumsum(VSWClike[i,])
  #set each column to have sum zero.
  #print(NEElike[i,1:500])
  #problem: there are way too many timepoints, so the exponentials of cumulative sums of likelihoods will be incredibly tiny. Need a more reasonable re-weighting method to make the quantiles for the particle filter. 
}
#set each column to have sum zero
NEEsum = colSums(NEElike)
LEsum = colSums(LElike)
VSWCsum = colSums(VSWClike)
NEEshift = array(1,dim(NEElike))
NEEshift = NEEshift%*%diag(NEEsum)
NEEshift = NEEshift/ne
NEElike = NEElike - NEEshift
LEshift = array(1,dim(LElike))
LEshift = LEshift%*%diag(LEsum)
LEshift = LEshift/ne
LElike = LElike - LEshift
VSWCshift = array(1,dim(VSWClike))
VSWCshift = VSWCshift%*%diag(VSWCsum)
VSWCshift = VSWCshift/ne
VSWClike = VSWClike - VSWCshift


#make sure the max and min are not too large

#NEEmax = apply(NEElike,2,max)
NEEmax = max(NEElike)
LEmax = max(LElike)
VSWCmax = max(VSWClike)
NEElike = 30*NEElike/NEEmax
LElike = 30*LElike/LEmax
VSWClike = 30*VSWClike/VSWCmax


#take out of log space
NEElike = exp(NEElike)
LElike = exp(LElike)
VSWClike = exp(VSWClike)

hist(NEElike[,ncol(NEElike)],main="Final Ensemble Weights(NEE)", xlab = "NEElike for the ensemble")
hist(LElike[,ncol(LElike)],main="Final Ensemble Weights(LE)", xlab = "NEElike for the ensemble")
hist(VSWClike[,ncol(VSWClike)],main="Final Ensemble Weights(Soil Moisture)", xlab = "NEElike for the ensemble")
#not showing anything because it is all zero.
``` 
Plotting the non-resampling particle filter

```{r}
## Non-resampling Particle Filter
## calculation of CI
#nobs = ncol(LAIlike)                     ## number of observations
NEEpf = matrix(NA, 3, nt)              ## storage [intervals,time]
LEpf = matrix(NA,3,nt)
VSWCpf = matrix(NA,3,nt)
wbar = apply(NEElike, 2, mean)           ## mean weight at each time point
wbar2 = apply(LElike,2,mean)
wbar3 = apply(VSWClike,2,mean)
for(i in 1:nt){
  ## calculate weighted median and CI
  vec1 = NEEcast[,i]
  vec2 = LEcast[,i]
  vec3 = VSWCcast[,i]
  wt = NEElike[, i] / wbar[i]
  wt2 = LElike[,i] / wbar2[i]
  wt3 = VSWClike[,i] / wbar3[i]
  NEEpf[, i] = wtd.quantile(NEEcast[, i], NEElike[, i] / wbar[i], c(0.025, 0.5, 0.975))  
  LEpf[, i] = wtd.quantile(LEcast[, i], LElike[, i] / wbar2[i], c(0.025, 0.5, 0.975)) 
  VSWCpf[, i] = wtd.quantile(VSWCcast[, i], VSWClike[, i] / wbar3[i], c(0.025, 0.5, 0.975)) 
}
```

```{r}

## plot original ensemble and PF with data (NEE)
col.pf   = c(col.alpha("lightGrey", 0.5), col.alpha("lightBlue", 0.5), 
             col.alpha("lightGreen", 0.5))                      ## color sequence
names.pf = c("ensemble", "non-resamp PF", "resamp PF")          ## legend names

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()
ci = apply(output[, , 1], 2, quantile, c(0.025, 0.5, 0.975))   ## calculate CI over ensemble members
ciEnvelope(time_predict,ci[1, ], ci[3, ], col = col.alpha("lightGrey", 0.5)) ## plot interval
lines(ci[2, ])


ciEnvelope(time_predict,NEEpf[1, ], NEEpf[3, ], col = col.pf[2])      ## non-resampling Particle Filter
lines(time_predict,NEEpf[2, ],col="green")

points(data_up$time,NEEnew,col="red")
#points(Mtime, LAIr)                                                   ## observations
#for(i in 1:length(LAIr)){                                             ## observation uncertainty
#  if(!is.na(QC[i])){
#    lines(rep(Mtime[i], 2), LAIr[i]+c(-1, 1) * LAIr.sd[i])            ## data is +/- 1 SD; NOT 95%
#  }
#}
legend("topleft", legend = names.pf[1:2], col = col.pf[1:2], lwd = 5)
```



### Plot for NEE (subset)
```{r}
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_NEE2()
ci = apply(output[, , 1], 2, quantile, c(0.025, 0.5, 0.975))   ## calculate CI over ensemble members
ciEnvelope(time_predict,ci[1, ], ci[3, ], col = col.alpha("lightGrey", 0.5)) ## plot interval
lines(ci[2, ])
ciEnvelope(time_predict,NEEpf[1, ], NEEpf[3, ], col = col.pf[2])      ## non-resampling Particle Filter
lines(time_predict,NEEpf[2, ],col="green")
points(data_up$time,NEEnew,col="red")
legend("topleft", legend = names.pf[1:2], col = col.pf[1:2], lwd = 5)
```

```{r}

## plot original ensemble and PF with data (LE)
col.pf   = c(col.alpha("lightGrey", 0.5), col.alpha("lightBlue", 0.5), 
             col.alpha("lightGreen", 0.5))                      ## color sequence
names.pf = c("ensemble", "non-resamp PF", "resamp PF")          ## legend names

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()
ci = apply(output[, , 2], 2, quantile, c(0.025, 0.5, 0.975))   ## calculate CI over ensemble members
ciEnvelope(time_predict,ci[1, ], ci[3, ], col = col.alpha("lightGrey", 0.5)) ## plot interval
lines(ci[2, ])


ciEnvelope(time_predict,LEpf[1, ], LEpf[3, ], col = col.pf[2])      ## non-resampling Particle Filter
lines(time_predict,LEpf[2, ],col="green")

points(data_up$time,LEnew,col="red")
#points(Mtime, LAIr)                                                   ## observations
#for(i in 1:length(LAIr)){                                             ## observation uncertainty
#  if(!is.na(QC[i])){
#    lines(rep(Mtime[i], 2), LAIr[i]+c(-1, 1) * LAIr.sd[i])            ## data is +/- 1 SD; NOT 95%
#  }
#}
legend("topleft", legend = names.pf[1:2], col = col.pf[1:2], lwd = 5)
```

### Plot for LE (subset)
```{r}
col.pf   = c(col.alpha("lightGrey", 0.5), col.alpha("lightBlue", 0.5), 
             col.alpha("lightGreen", 0.5))                      ## color sequence
names.pf = c("ensemble", "non-resamp PF", "resamp PF")          ## legend names
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_LE2()
ci = apply(output[, , 2], 2, quantile, c(0.025, 0.5, 0.975))   ## calculate CI over ensemble members
ciEnvelope(time_predict,ci[1, ], ci[3, ], col = col.alpha("lightGrey", 0.5)) ## plot interval
lines(ci[2, ])
ciEnvelope(time_predict,LEpf[1, ], LEpf[3, ], col = col.pf[2])      ## non-resampling Particle Filter
lines(time_predict,LEpf[2, ],col="green")
points(data_up$time,LEnew,col="red")
legend("topleft", legend = names.pf[1:2], col = col.pf[1:2], lwd = 5)
```

```{r}

## plot original ensemble and PF with data (VSWC)
col.pf   = c(col.alpha("lightGrey", 0.5), col.alpha("lightBlue", 0.5), 
             col.alpha("lightGreen", 0.5))                      ## color sequence
names.pf = c("ensemble", "non-resamp PF", "resamp PF")          ## legend names

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()
ci = apply(output[, , 3], 2, quantile, c(0.025, 0.5, 0.975))   ## calculate CI over ensemble members
ciEnvelope(time_predict,ci[1, ], ci[3, ], col = col.alpha("lightGrey", 0.5)) ## plot interval
lines(ci[2, ])


ciEnvelope(time_predict,VSWCpf[1, ], VSWCpf[3, ], col = col.pf[2])      ## non-resampling Particle Filter
lines(time_predict,VSWCpf[2, ],col="green")

points(data_up$time,VSWCnew,col="red")
#points(Mtime, LAIr)                                                   ## observations
#for(i in 1:length(LAIr)){                                             ## observation uncertainty
#  if(!is.na(QC[i])){
#    lines(rep(Mtime[i], 2), LAIr[i]+c(-1, 1) * LAIr.sd[i])            ## data is +/- 1 SD; NOT 95%
#  }
#}
legend("topleft", legend = names.pf[1:2], col = col.pf[1:2], lwd = 5)
```

### Plot for VSWC (subset)
```{r}
## plot original ensemble and PF with data (VSWC)
col.pf   = c(col.alpha("lightGrey", 0.5), col.alpha("lightBlue", 0.5), 
             col.alpha("lightGreen", 0.5))                      ## color sequence
names.pf = c("ensemble", "non-resamp PF", "resamp PF")          ## legend names
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_VSWC2()
ci = apply(output[, , 3], 2, quantile, c(0.025, 0.5, 0.975))   ## calculate CI over ensemble members
ciEnvelope(time_predict,ci[1, ], ci[3, ], col = col.alpha("lightGrey", 0.5)) ## plot interval
lines(ci[2, ])
ciEnvelope(time_predict,VSWCpf[1, ], VSWCpf[3, ], col = col.pf[2])      ## non-resampling Particle Filter
lines(time_predict,VSWCpf[2, ],col="green")
points(data_up$time,VSWCnew,col="red")
legend("topleft", legend = names.pf[1:2], col = col.pf[1:2], lwd = 5)
```

```{r}
## save data

Fantastic4casters <- data.frame(time = time_predict, statistic = matrix("mean",1679,1), siteID = matrix("KONZ",1679,1), 
                                forecast = matrix(1,1679,1), data_assimilation=matrix(0,1679,1), nee=NEEpf[2,], le=LEpf[2,],
                                vswc=VSWCpf[2,])

newFilename <- sprintf("%s.Rdata","Prediction")
newFilename <- paste(dataPath, newFilename, sep="", collapse = NULL)
save(Fantastic4casters,file=newFilename)

```