# gam_daily, kriging_daily,gwr_daily,gam_CAI,gam_fusion,kriging_fusion,gwr_fusion and other options added.

#For multiple time scale methods, the interpolation is done first at the monthly time scale then delta surfaces are added.

#AUTHOR: Benoit Parmentier                                                                        

#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#568--     

#

# TO DO:

  #5)infile_covariates: raster covariate brick, tif file

  #6)covar_names: covar_names #remove at a later stage...

  #7)var: variable being interpolated-TMIN or TMAX

  #9)CRS_locs_WGS84

  #10)screen_data_training

  #

  #19)prop_minmax_month

  #20)dates_selected

  #

  #21)list_models: model formulas in character vector

  #22)lst_avg: LST climatology name in the brick of covariate--change later

  #25)script_path: path to script

  #26)interpolation_method: c("gam_fusion","gam_CAI") #other otpions to be added later

  #27) use_clim_image

  #28) join_daily

  #29)list_models2: models' formulas as string vector for daily devation

  #30)interp_method2: intepolation method for daily devation step

  ###Loading R library and packages     

  interp_method2 <- list_param_raster_prediction$interp_method2

  lst_avg<-list_param_raster_prediction$lst_avg

  script_path<-list_param_raster_prediction$script_path

  interpolation_method<-list_param_raster_prediction$interpolation_method

  screen_data_training <-list_param_raster_prediction$screen_data_training

  use_clim_image <- list_param_raster_prediction$use_clim_image # use predicted image as a base...rather than average Tmin at the station for delta

  join_daily <- list_param_raster_prediction$join_daily # join monthly and daily station before calucating delta

  ###################### START OF THE SCRIPT ########################

  cat("Starting script process time:",file=log_fname,sep="\n",append=TRUE)

  cat(as.character(time1),file=log_fname,sep="\n",append=TRUE)    

  ############### READING INPUTS: DAILY STATION DATA AND OTHER DATASETS  #################

  #Should clean this up, reduce the number of if

  if (dates_selected==""){

  names(s_raster)<-covar_names               #Assigning names to the raster layers: making sure it is included in the extraction

  #Reading monthly data

  #construct date based on input end_year !!!

  day_tmp <- rep("15",length=nrow(dst))

  year_tmp <- rep(as.character(end_year),length=nrow(dst))

  sampling_month_obj <- sampling_training_testing(list_param_sampling)

  ########### PREDICT FOR MONTHLY SCALE  ##################

  #First predict at the monthly time scale: climatology

  cat("Predictions at monthly scale:",file=log_fname,sep="\n", append=TRUE)

  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),

      file=log_fname,sep="\n")

  t1<-proc.time()

  j=12

  #browser() #Missing out_path for gam_fusion list param!!!

  #if (interpolation_method=="gam_fusion"){

  if (interpolation_method %in% c("gam_fusion","kriging_fusion","gwr_fusion")){

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

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

      tmp<-clim_method_mod_obj[[i]]$clim

      list_tmp[[i]]<-tmp

    }

    clim_yearlist<-list_tmp

  }

  if (interpolation_method %in% c("gam_CAI","kriging_CAI", "gwr_CAI")){

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

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

      tmp<-clim_method_mod_obj[[i]]$clim

  }

  t2<-proc.time()-t1

  cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)

  ################## PREDICT AT DAILY TIME SCALE #################

  #Predict at daily time scale from single time scale or multiple time scale methods: 2 methods availabe now

  #put together list of clim models per month...

  #rast_clim_yearlist<-list_tmp

  #Second predict at the daily time scale: delta

  #method_mod_obj<-mclapply(1:1, runGAMFusion,mc.preschedule=FALSE,mc.cores = 1) #This is the end bracket from mclapply(...) statement

  cat("Predictions at the daily scale:",file=log_fname,sep="\n", append=TRUE)

  t1<-proc.time()

  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),

    #test <- run_prediction_daily_deviation(1,list_param=list_param_run_prediction_daily_deviation) #This is the end bracket from mclapply(...) statement

    #test <- mclapply(1:9,list_param=list_param_run_prediction_daily_deviation,run_prediction_daily_deviation,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

  }

  #TODO : Same call for all functions!!! Replace by one "if" for all daily single time scale methods...

  if (interpolation_method=="gam_daily"){

    #input a list:note that ghcn.subsets is not sampling_obj$data_day_ghcn

    names(list_param_run_prediction_gam_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","screen_data_training","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")

    #test <- runGAM_day_fun(1,list_param_run_prediction_gam_daily)

    #method_mod_obj<-mclapply(1:22,list_param=list_param_run_prediction_gam_daily,runGAM_day_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

    list_param_run_prediction_kriging_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_kriging_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")

    #test <- runKriging_day_fun(1,list_param_run_prediction_kriging_daily)

    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

    save(method_mod_obj,file= file.path(out_path,paste("method_mod_obj_",interpolation_method,"_",y_var_name,out_prefix,".RData",sep="")))

    list_param_run_prediction_gwr_daily <-list(i,s_raster,covar_names,lst_avg,list_models,dst,var,y_var_name, sampling_obj,interpolation_method,out_prefix,out_path)

    names(list_param_run_prediction_gwr_daily)<-c("list_index","covar_rast","covar_names","lst_avg","list_models","dst","var","y_var_name","sampling_obj","interpolation_method","out_prefix","out_path")

    #test <- run_interp_day_fun(1,list_param_run_prediction_gwr_daily)

    #method_mod_obj<-mclapply(1:22,list_param=list_param_run_prediction_gwr_daily,run_interp_day_fun,mc.preschedule=FALSE,mc.cores = 11) #This is the end bracket from mclapply(...) statement

    #method_mod_obj<-mclapply(1:18,list_param=list_param_run_prediction_kriging_daily,runKriging_day_fun,mc.preschedule=FALSE,mc.cores = 9) #This is the end bracket from mclapply(...) statement

  ############### NOW RUN VALIDATION #########################

  #SIMPLIFY THIS PART: one call

  cat("Validation step:",file=log_fname,sep="\n", append=TRUE)

  t1<-proc.time()

  cat(paste("Local Date and Time: ",as.character(Sys.time()),sep=""),

    #debug(calculate_accuracy_metrics)

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation)

    save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep="")))

    t2<-proc.time()-t1

    cat(as.character(t2),file=log_fname,sep="\n", append=TRUE)

  }

  ### Run monthly validation if multi-time scale methods and add information to daily...

  if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){

    multi_time_scale <- TRUE

    i<-1

    #debug(calculate_accuracy_metrics)

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation)

    save(validation_mod_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_obj_",y_var_name,out_prefix,".RData",sep="")))

    ### monthly time scale

    #debug(calculate_accuracy_metrics)    

    #test_val2 <-calculate_accuracy_metrics(1,list_param_validation_month)

    #test_val<-calculate_accuracy_metrics(1,list_param_validation)

    save(validation_mod_month_obj,file= file.path(out_path,paste(interpolation_method,"_validation_mod_month_obj_",y_var_name,out_prefix,".RData",sep="")))

    tb_month_diagnostic_s$method_interp <- interpolation_method #add type of interpolation...out_prefix too??

  }

  #################### ASSESSMENT OF PREDICTIONS: PLOTS OF ACCURACY METRICS ###########

  ##Create data.frame with validation and fit metrics for a full year/full numbe of runs

  tb_diagnostic_v<-extract_from_list_obj(validation_mod_obj,"metrics_v") 

  #tb_diagnostic_v contains accuracy metrics for models sample and proportion for every run...if full year then 365 rows maximum

  ################### PREPARE RETURN OBJECT ###############

  #Will add more information to be returned

  if (interpolation_method %in% c("gam_CAI","kriging_CAI","gwr_CAI","gam_fusion","kriging_fusion","gwr_fusion")){

    raster_prediction_obj<-list(clim_method_mod_obj,method_mod_obj,validation_mod_obj,validation_mod_month_obj, tb_diagnostic_v,

                                tb_diagnostic_s,tb_month_diagnostic_v,tb_month_diagnostic_s,summary_metrics_v,summary_month_metrics_v)

    save(raster_prediction_obj,file= file.path(out_path,paste("raster_prediction_obj_",interpolation_method,"_", y_var_name,out_prefix,".RData",sep="")))

  }

  return(raster_prediction_obj)

}

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