##VALIDATION-ACCURACY ASSESSMENT

#The interpolation is done first at the monthly add delta.

#AUTHOR: Benoit Parmentier                                                                        

#Change this to allow explicitly arguments...

#Arguments: 

#1)list of climatology files for all models...(365*nb of models)

#2)data_s:training

#3)data_v:testing

#6)dst: data at the monthly time scale

#Function used in the script

  calc_val_metrics<-function(x,y){

    #This functions calculates accurayc metrics on given two vectors.

    #Arguments: list of fitted models, raster stack of covariates

    #Output: spatial grid data frame of the subset of tiles

    #s_sgdf<-as(r_stack,"SpatialGridDataFrame") #Conversion to spatial grid data frame

    residuals<-x-y

    mae<-mean(abs(residuals),na.rm=T)

    rmse<-sqrt(mean((residuals)^2,na.rm=T))

    me<-mean(residuals,na.rm=T)

    r<-cor(x,y,use="complete")

    m50<-median(residuals,na.rm=T)

    metrics_dat<-as.data.frame(cbind(mae,rmse,me,r,m50))

    names(metrics_dat)<-c("mae","rmse","me","r","m50")

    metrics_obj<-list(metrics_dat,as.data.frame(residuals))

    names(metrics_obj)<-c("metrics_dat","residuals")

    return(metrics_obj)

  }

  calc_val_metrics_rast <-function(df,y_ref,pred_names){

    #1) df: data frame containing the observed and predicted variables (data_s or data_v)

    #2) y_ref: observed variable correspond to y_var_name??

    #3) pred_names: models run containig predicted values

    # library

    library(maptools)

    ## START SCRIPT

    list_metrics<-vector("list",length(pred_names))

    list_residuals<-vector("list",length(pred_names))

    names(list_metrics)<-pred_names

    names(list_residuals)<-pred_names

    for (j in 1:length(pred_names)){

      pred_var<-pred_names[j]

      metrics<-calc_val_metrics(df[[pred_var]],df[[y_ref]])

      list_metrics[[j]]<-metrics[[1]]

      list_residuals[[j]]<-metrics[[2]]

    }

    metrics_df<-do.call(rbind,list_metrics)

    metrics_df$pred_mod <- pred_names #adding name column

    residuals_df<-do.call(cbind,list_residuals) #creating data frame for residuals

    names(residuals_df)<-paste("res",pred_names,sep="_")

    accuracy_obj<-list(metrics_df,residuals_df) #output object

    names(accuracy_obj)<-c("metrics","residuals") 

    return(accuracy_obj)

  }  

  names_mod <- names(day_list) #names of the predicted variables

  data_s <- list_param$list_data_s[[i]]

  sampling_dat_day <- list_param$list_sampling_dat[[i]]

  ## Now create the stack

  rast_day_mod <- stack(day_list)

  names(rast_day_mod) <- names(day_list)

  nv <- nrow(data_v)

  #add sampling dat info...

  N=length(names_mod)

  if (nv > 0){

    run_info<-cbind(sampling_dat_day,n=nv)

    run_info[rep(seq_len(nrow(run_info)), each=N),] #repeating same row n times

    extract_data_v<-extract(rast_day_mod,data_v,df=TRUE)

    data_v <-spCbind(data_v,extract_data_v) #should match IDs before joining for good practice    

    data_v<-spCbind(data_v,metrics_v_obj$residuals)

  }

  extract_data_s<-extract(rast_day_mod,data_s,df=TRUE)  

  #prepare output object

  if (nv > 0){

    validation_obj<-list(metrics_s_df,metrics_v_df,data_s,data_v)

    names(validation_obj)<-c("metrics_s","metrics_v","data_s","data_v")

  }else{

    validation_obj<-list(metrics_s_df,data_s)

    names(validation_obj)<-c("metrics_s","data_s")

  }

  return(validation_obj)

}

extract_from_list_obj<-function(obj_list,list_name){

  list_tmp<-vector("list",length(obj_list))

  for (i in 1:length(obj_list)){

    tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame

    list_tmp[[i]]<-tmp

  }

  tb_list_tmp<-do.call(rbind,list_tmp) #long rownames

  return(tb_list_tmp) #this is  a data.frame

}

extract_list_from_list_obj<-function(obj_list,list_name){

  #Create a list of an object from a given list of object using a name prodived as input

  list_tmp<-vector("list",length(obj_list))

  for (i in 1:length(obj_list)){

    tmp<-obj_list[[i]][[list_name]] #double bracket to return data.frame

    list_tmp[[i]]<-tmp

  }

  return(list_tmp) #this is  a data.frame

}

  ### Start script

          #measure=mod_var, 

          id=c("date","pred_mod","prop"),

          na.rm=F)

  #avg_tb<-cast(t,pred_mod~variable,mean)

  #mod_names<-sort(unique(tb$pred_mod)) #kept for clarity

  tb_mod_list<-lapply(mod_names, function(k) subset(tb, pred_mod==k)) #this creates a list of 5 based on models names

    mod_metrics<-do.call(cbindX,tb_mod_list) #column bind the list??

  }else{

    mod_metrics<-tb_mod_list[[1]]

  }

  rows_total<-lapply(tb_mod_list,nrow)

    metric_ac<-metric_names[j]

    #browser()

    test<-mod_metrics[mod_var]

    #boxplot(test,outline=FALSE,horizontal=FALSE,cex=0.5,

    #        ylab=paste(metric_ac,"in degree C",sep=" "))

    axis(1, labels = FALSE)

    ## Create some text labels

    labels <- labels<- names(test)

    ## Plot x axis labels at default tick marks

    text(1:ncol(test), par("usr")[3] - 0.25, srt = 45, adj = 1,

         labels = labels, xpd = TRUE)

    axis(2)

    box()

    title(paste(metric_ac,"for",y_var_name,sep=" "),cex=0.8)

  }

  median_tb$n<-rows_total

}

  #Generate boxplot per month for models and accuracy metrics

  #Input parameters:

  #1) df: data frame containing accurayc metrics (RMSE etc.) per day)

  #2) metric_names: metrics used for validation

  #3) out_prefix

  #

  #################

  ## BEGIN

  tb_diagnostic$month<-strftime(date_f, "%m")          # current month of the date being processed

  mod_names<-sort(unique(tb_diagnostic$pred_mod)) #models that have accuracy metrics

  tb_mod_list<-lapply(mod_names, function(k) subset(tb_diagnostic, pred_mod==k)) #this creates a list of 5 based on models names

  names(tb_mod_list)<-mod_names

  t<-melt(tb_diagnostic,

          #measure=mod_var, 

          id=c("date","pred_mod","prop","month"),

          na.rm=F)

    mod_metrics <-tb_mod_list[[k]]

    current_mod_name<- mod_names[k]

    for (j in 1:length(metric_names)){    

      test<-mod_metrics[col_selected]

      #boxplot(test[[metric_ac]]~test[[c("month")]],outline=FALSE,horizontal=FALSE,cex=0.5,

      #        ylab=paste(metric_ac,"in degree C",sep=" "))

      axis(1, labels = FALSE)

      ## Create some text labels

      labels <- month.abb # abbreviated names for each month

      ## Plot x axis labels at default tick marks

      text(1:length(labels), par("usr")[3] - 0.25, srt = 45, adj = 1,

           labels = labels, xpd = TRUE)

      axis(2)

      box()

  names(summary_month_obj)<-c("tb_list","metric_month_avg","metric_month_sd")

  return(summary_month_obj)  

############ END OF SCRIPT #########