analysis_chl.Rmd

---
title: "Lake Geneva DO Analysis"
author: "Ethan Deyle, Damien Bouffard, et al."
date: '`r format(Sys.Date(), "%B %d, %Y")`'
output:
  html_document: 
    toc: true
    number_sections: true
---

# Intro and Setup

```{r Packages, warning=FALSE}
library(readr)
library(tidyverse)
library(ggplot2)
library(gridExtra)
library(lubridate)

library("rEDM",lib.loc = "./LIB")

rm(list = ls())
```

```{r Help Functions}
source('./FUNCTIONS/F_EDM_analysis.R')
source('./FUNCTIONS/F_data_and_plots.R')
```

```{r Import Data}
data_file <- './DATA/lake_geneva_proc_extended_2019_07_30.Rdata'
```

## Basic Plots

### Univariate EDM Plots for Chl

```{r Plot rho vs tp}
load(data_file)

target_var <- "chl"

data_univar <- lake_geneva_interp %>% filter(date>="1976-05-08") %>% pull(chl)

out_simplex <- simplex( data_univar,E=1:15,tau=1,tp=1)

ggplot(out_simplex,aes(x=E,y=rho)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(1,10),breaks = seq(2,10,2))

out_tp <- do.call(bind_rows,
                  lapply(1:48,function(tp) {
                    simplex(lake_geneva_interp[,target_var],
                            pred = c(1,NROW(lake_geneva_interp)-48+tp),
                            tau=1,tp=tp,E=12,
                            exclusion_radius = 12)
                  }))

out_smap <- s_map(data_univar,tau=1,tp=1,E=12)

ggplot(out_smap,aes(x=theta,y=rho)) + geom_line() + theme_bw() + scale_x_continuous(limits=c(0,8),breaks = seq(2,10,2))
ggplot(out_tp,aes(x=tp,y=rho)) + geom_line()
```

### Scatterplots vs. chl

```{r Plot scatter with chl}
load(data_file)

data_long <- tidyr::gather(lake_geneva_interp, variable, value, -c(date,chl))
data_plot <- ggplot(data_long, aes(y=chl, x=value)) + 
  geom_path(color = 'grey70', lwd = 0.1,alpha = .8) +  
  geom_point(color = "red", size = 0.2) +
  facet_wrap(~variable, scales = "free_x") + 
  theme_bw()
# coord_flip()

print(data_plot)
```
### Time Scales

Do basic statistical analysis of the candidate time series to identify the de-correlation time (first time in the auto-correlation function drops below 1/e ~ .37). This helps indicate if the time series are over-sampled and should be treat with tau/tp > 1.

```{r Plot chl ACF}
acf.chl <- acf(lake_geneva_interp$chl,na.action=na.pass,lag.max = 48,plot = FALSE)

plot(acf.chl,axes=FALSE,main = "",xlab= "")
axis(1,at=seq(0,48,by=12))
axis(2)
title(main='Chl ACF',xlab='months')
box()
```

## Multivariate EDM

Get three baselines of forecast skill for the target variable: univariate EDM, seasonal, and multiview.

```{r mEDM Setup}
load(data_file)

block_chl_mEDM_raw <- lake_geneva_interp %>%
    mutate(year_sine=sin(2*pi*yday(date)/365)) %>%
    mutate(year_cosine=cos(2*pi*yday(date)/365)) %>%
    filter(date>="1976-05-08") %>%
  select(-pp)

block_chl_mEDM <- lake_geneva_interp %>%
    mutate(year_sine=sin(2*pi*yday(date)/365)) %>%
    mutate(year_cosine=cos(2*pi*yday(date)/365)) %>%
    filter(date>="1976-05-08") %>%
    mutate_at(vars(-date),funs(./sd(.,na.rm=TRUE))) %>%
  select(-pp)
```

```{r mEDM Baselines}
file_baselines <- "./OUTPUTS/mEDM_chl_baselines.Rdata"

if(!file.exists(file_baselines)){

results_mEDM_multiview_fullfit <-  block_chl_mEDM %>% 
        select(-T_surf,-Ptot_lake,-Ptot_epi,-T_air,-year_cosine) %>%
        multiview(lib=c(1,nrow(block_chl_mEDM)),pred=c(1,nrow(block_chl_mEDM)),
                  target_column="chl",
                  max_lag=3,E=4,
                  first_column_time = T)

results_mEDM_seasonal <- do_mEDM_models(block=block_chl_mEDM,  list( c('year_sine','year_cosine') ),target_var='chl',tp=1)$multi_stats[[1]]

E.univar <- 12
results_univar <- s_map(data_univar,tau=1,tp=1,E=E.univar)

save(results_univar,
     results_mEDM_seasonal,
     results_mEDM_multiview_fullfit,
     file=file_baselines)
}
```

Forecast chl using greedy approach from the two focal variables, nutrient loading (PO4_lake) and temperature (T_surf_model).

```{r mEDM Greedy}
file_greedy <- "./OUTPUTS/mEDM_chl_greedy.Rdata"

if(!file.exists(file_greedy)){
  L_candidate_variables <- c('year_sine','chl','rhone','Z_thermo_model','PO4_epi')
  
  results_mEDM_greedy <- do_mEDM_greedy(block_chl_mEDM,embed0=c('T_surf_model','PO4_lake'),L_candidate_variables,target_col="chl",max_E=5)
  save(results_mEDM_greedy,file=file_greedy)
}
```

Plot the mEDM model selected with the greedy algorithm along with the baseline models.

```{r}
load('./OUTPUTS/mEDM_chl_baselines.Rdata') # results_univar,results_mEDM_multiview_fullfit,results_mEDM_multiview_nfold
load('./OUTPUTS/mEDM_chl_greedy.Rdata') # results_mEDM_greedy

ggplot(data=results_mEDM_greedy,aes(x=theta,y=rho,color=embedding_label)) + 
    geom_line(lwd=1) +
    geom_line(data=results_univar,aes(color="univariate")) +
    geom_line(data=results_mEDM_seasonal,aes(color="seasonal")) +
    geom_hline(aes(yintercept=results_mEDM_multiview_fullfit$rho,color="multiview"))

```

Same with pp

```{r}
load(data_file)

lake_geneva_interp <- lake_geneva_interp %>%
    filter(date >= ymd("1983-01-01")) %>%
    filter(date <= ymd("2012-01-01")) %>%
    mutate(year_sine=sin(2*pi*yday(date)/365)) %>%
    mutate(year_cosine=cos(2*pi*yday(date)/365)) %>%
  select(-date)

L_models <- list( c('PO4_epi','PO4_lake','T_air','pp'),
                  c('PO4_lake','T_air','pp'),
                  c('PO4_lake','T_surf','pp'),
                  c('T_air','pp'),
                  c('year_sine','year_cosine'),
                  c('year_sine','year_cosine','PO4_lake'))

target_var <- 'pp'

lake_geneva_interp_norm <- lake_geneva_interp %>%
  mutate(chl_lag_1 = lag(chl,1)) %>%
  mutate_all(funs(./sd(.,na.rm=TRUE)))

results_mEDM_exp1 <- do_mEDM_models(block=lake_geneva_interp_norm, L_models,target_var=target_var,tp=1)

g_multi <- map(results_mEDM_exp1$multi_preds,function(df){
  ggplot(df,aes(x=time)) + geom_line(aes(y=obs,col='obs')) + geom_line(aes(y=pred,col='pred'))
})

g_errors_vs_oxy <- map(results_mEDM_exp1$multi_preds,function(df){
    df %>% bind_cols(lake_geneva_interp) %>%
  ggplot(aes(x=oxydeep,y=(pred-obs))) + geom_point(size=.75) + stat_smooth(method="lm") + ylab('Error') 
})

g_errors_vs_deltaoxy <- map(results_mEDM_exp1$multi_preds,function(df){
    df %>% bind_cols(lake_geneva_interp) %>%
  ggplot(aes(x=obs,y=(pred-obs))) + geom_point(size=.75) + stat_smooth() + ylab('Error') 
})

do.call(bind_rows,results_mEDM_exp1$multi_stats) %>%
    mutate(embedding=label_embeddings(embedding,
                                      names(lake_geneva_interp_norm),
                                      dictionary = dictionary_extended) %>% as.factor()) %>%
ggplot(aes(x=theta,y=rho,color=embedding)) + geom_line(lwd=1) + theme_bw()

do.call(bind_rows,results_mEDM_exp1$multi_stats) %>%
    mutate(embedding=label_embeddings(embedding,
                                      names(lake_geneva_interp_norm),
                                      dictionary = dictionary_extended) %>% as.factor()) %>%
ggplot(aes(x=theta,y=rmse,color=embedding)) + geom_line(lwd=1) + theme_bw()
```

### Look at forecast decay

We look at S-map forecast skill as a function of the prediction interval (number of time-steps).

```{r}
load(data_file)

lake_geneva_interp <- lake_geneva_interp %>%
    mutate(year_sine=sin(2*pi*yday(date)/365)) %>%
    mutate(year_cosine=cos(2*pi*yday(date)/365)) %>%
    filter(date>="1976-05-08") %>%
    mutate(T_diff_model=T_air-T_surf_model)


L_models <- list( c('T_surf_model','PO4_lake','year_sine'),
                  c('T_surf_model','PO4_lake','year_sine','chl'),
                  c('year_cosine','year_sine'))

sim_col <- "chl"
phys_vars <- c("T_surf","Z_thermo","wind_speed","T_air","rhone","T_surf_model","Z_thermo_model","year_sine","year_cosine")
tp_max <- 3*12

lake_geneva_sim <- lake_geneva_interp %>%
    rename(time=date) %>%
    # mutate(oxydeep_delta = c(NA,diff(oxydeep))) %>%
    mutate_at(vars(phys_vars),funs(lead(.,1))) %>%
    filter(complete.cases(.))

I_start_points <- 1:( NROW(lake_geneva_sim) - tp_max)



results_sim_tp_smap <- map_df(L_models,function(L_model_i){
    
    block_model_i <- lake_geneva_sim[,c('time',union(L_model_i,sim_col))]
    
    results_sim <- map_df(I_start_points,function(i_start){
        
        block_sim_i <- block_model_i[i_start + 0:tp_max,]
        
        out_sim_i <- sim_EDM_smap(block_model_i,
                                       block_sim_i,
                                       tp_max,
                                       tp = 1,
                                       lib_train = c(1,NROW(block_model_i)),
                                       pred_sim = c(1,NROW(block_model_i)),
                                       sim_col = sim_col, # could set to which columns are NAs
                                       normalize = TRUE,
                                       predictor_col=L_model_i,
                                       theta = c(0,.1,.5,1,2,3,4,5,6,7,8),
                                       keep_trajectories = TRUE,
                                       first_column_time = TRUE,
                                       bootstrap = FALSE,
                                       exclusion_radius = 90
                                       ) %>%
            rename(date=time) %>%
            group_by(type) %>%
            mutate(tp=row_number()-1) %>%
            ungroup()
    }
    
    )
    
    results_sim_chl <- full_join( results_sim %>%
                                      filter(type=='sim') %>%
                                      select(date,tp,chl) %>%
                                      rename(pred=chl) ,
                                  results_sim %>%
                                      filter(type=='true') %>%
                                      select(date,tp,chl) %>%
                                      rename(obs=chl),
                                  by = c("date","tp")) %>%
                        mutate(embedding=label_embeddings(paste(1:length(L_model_i),collapse=","),
                                                          L_model_i,
                                                          dictionary=dictionary_extended))
    
}) # map(L_models_)


results_sim_tp_smap %>%
    group_by(tp,embedding) %>%
    summarise(compute_stats(obs,pred) %>% list()) %>%
    unnest() %>%
    ggplot(aes(x=tp,y=rho)) + geom_line(aes(color=embedding))



results_sim_tp_smap %>%
    group_by(tp,embedding) %>%
    summarise(compute_stats(obs,pred) %>% list()) %>%
    unnest() %>%
    ggplot(aes(x=tp,y=mae)) + geom_line(aes(color=embedding))
```

We can visualize time series of 6m, 1yr, 18m, 2yr forecasts to check that the S-map models are producing reasonable behavior even as forecast skill decreases.

```{r}
results_sim_tp_exp <- "./OUTPUTS/mEDM_chl_sim_tp.Rdata"
load(data_file)

lake_geneva_interp <- lake_geneva_interp %>%
    mutate(year_sine=sin(2*pi*yday(date)/365)) %>%
    mutate(year_cosine=cos(2*pi*yday(date)/365)) %>%
    filter(date>="1976-05-08") %>%
    mutate(T_diff_model=T_air-T_surf_model)

L_models <- list( c('T_surf_model','PO4_lake','year_sine'),
                  c('T_surf_model','PO4_lake','year_sine','chl'),
                  c('year_cosine','year_sine'))

sim_col <- "chl"

v_plot_times <- c(6,24,36)
tp_max <- max(v_plot_times)+1

lake_geneva_sim <- lake_geneva_interp %>%
    rename(time=date) %>%
    mutate(oxydeep = lag(oxydeep,1)) %>%
    filter(complete.cases(.))

I_start_points <- 1:( NROW(lake_geneva_sim) - tp_max)

if(!file.exists(results_sim_tp_exp)){
  results_sim_tp_exp <- map_df(L_models,function(L_model_i){
    
    block_model_i <- lake_geneva_sim[,c('time',union(sim_col,L_model_i))]
    
    results_sim <- map_df(I_start_points,function(i_start){
      
      block_sim_i <- block_model_i[i_start + 0:tp_max,]
      
      out_sim_i <- sim_EDM_smap(block_model_i,
                                block_sim_i,
                                tp_max,
                                tp = 1,
                                lib_train = c(1,NROW(block_model_i)),
                                pred_sim = c(1,NROW(block_model_i)),
                                sim_col = sim_col, # could set to which columns are NAs
                                predictor_col=L_model_i,
                                normalize = TRUE,
                                theta = c(0,.1,.5,1,2,3,4,5,6,7,8),
                                keep_trajectories = TRUE,
                                first_column_time = TRUE,
                                bootstrap = FALSE,
                                exclusion_radius = 90) %>%
        
        
        rename(date=time) %>%
        group_by(type) %>%
        mutate(tp=row_number()-1)
    }
    
    )
    
    results_sims <- results_sim %>% ungroup() %>%
      filter(type=='sim') %>%
      select(date,tp,chl) %>%
      mutate(embedding=label_embeddings(paste(1:length(L_model_i),collapse=","),
                                        L_model_i,
                                        dictionary=dictionary_extended)) %>%
      filter(tp %in% v_plot_times)
    
  }) # map(L_models_)
  
  save(results_sim_tp_exp,file=file_sim_tp)
} # if(file.exists())

load(file_sim_tp)

results_sim_tp_exp %>%
    ggplot(aes(x=date,y=chl)) + geom_line(data=lake_geneva_interp,aes(color="OBSERVED")) + geom_line(aes(color=embedding)) + facet_grid(tp~embedding)

results_sim_tp_exp %>% left_join(lake_geneva_interp %>% select(date,chl),by="date") %>% group_by(tp,embedding) %>% summarise(rho=cor(chl.x,chl.y))    
```