#Analyses, figures, tables and data are also produced in the script.

#AUTHOR: Benoit Parmentier 

#CREATED ON: 03/23/2014  

#Version: 5

#PROJECT: Environmental Layers project     

#COMMENTS: Initial commit, script based on part 2 of assessment, will be modified further for overall assessment 

#in_dir <- "" #PARAM 0

#y_var_name <- "dailyTmax" #PARAM1

#interpolation_method <- c("gam_CAI") #PARAM2

#out_suffix <- "output_run10_1000x3000_global_analyses_02102015"

#out_suffix <- "run10_1500x4500_global_analyses_pred_1992_10052015" #PARAM3

#out_dir <- "/data/project/layers/commons/NEX_data/output_run10_1500x4500_global_analyses_pred_1992_10052015" #PARAM4

  setwd(out_dir)

  list_outfiles_names <- vector("list", length=23) #collect names of output files

  counter_fig <- 0 #index of figure to collect outputs

  list_tb_fname <- list.files(path=file.path(in_dir,in_dir_list),"tb_diagnostic_v_NA_.*._run_global_analyses_pred_.*._reg4.txt",full.names=T)

  list_df_fname <- list.files(path=file.path(in_dir,in_dir_list),"df_tile_processed_.*._run_global_analyses_pred_.*._reg4.txt",full.names=T)

  list_summary_metrics_v_fname <- list.files(path=file.path(in_dir,in_dir_list),"summary_metrics_v2_NA_.*._run_global_analyses_pred_.*._reg4.txt",full.names=T)

  list_tb <- lapply(list_tb_fname,function(x){read.table(x,stringsAsFactors=F,sep=",")})

  tb <- do.call(rbind,list_tb)

  list_df <- lapply(list_df_fname,function(x){read.table(x,stringsAsFactors=F,sep=",")})

  df_tile_processed <- do.call(rbind,list_df)  

  list_summary_metrics_v <- lapply(list_summary_metrics_v_fname,function(x){read.table(x,stringsAsFactors=F,sep=",")})

  summary_metrics_v <- do.call(rbind,list_summary_metrics_v)  

  tb <-  read.table(list_tb[1],stringsAsFactors=F,sep=",")

  try(tb$reg <- tb$reg.x)

  try(summary_metrics_v$year_predicted <- summary_metrics_v$year_predicted.x)

  try(summary_metrics_v$reg <- summary_metrics_v$reg.x)  

  #tb_all <- tb

  #summary_metrics_v_all <- summary_metrics_v 

  #df_tiled_processed <- na.omit(df_tile_processed) #remove other list of folders irrelevant

  #list_shp_reg_files <- df_tiled_processed$shp_files

  #list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,

  #          as.character(df_tile_processed$tile_coord),"shapefiles",basename(list_shp_reg_files))

  #list_shp_reg_files <- file.path("/data/project/layers/commons/NEX_data/",out_dir,

  table(tb$pred_mod)

  table(tb$index_d)

  table(subset(tb,pred_mod=="mod1")$index_d)

  #aggregate()

  tb$predicted <- 1

  test <- aggregate(predicted~pred_mod+tile_id,data=tb,sum)

  as.character(unique(test$tile_id)) #141 tiles

  ##### Figure 8: Breaking down accuracy by regions!! #####

  #summary_metrics_v <- merge(summary_metrics_v,df_tile_processed,by="tile_id")

  ## Figure 8a

  res_pix <- 480

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  p<- bwplot(rmse~pred_mod | reg, data=tb,

  print(p)

  dev.off()

  png(filename=paste("Figure8b_boxplot_overall_separated_by_region_scaling_","_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  p<- bwplot(rmse~pred_mod | reg, data=tb,ylim=c(0,5),

  print(p)

  dev.off()

  list_outfiles[[counter_fig+1]] <- paste("Figure8b_boxplot_overall_separated_by_region_scaling_",model_name[i],"_",out_suffix,".png",sep="")

  counter_fig <- counter_fig + 1

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

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

  #### Figure 9: plotting boxplot by year and regions ###########

  res_pix <- 480

  col_mfrow <- 1

  row_mfrow <- 1

  png(filename=paste("Figure9a_boxplot_overall_separated_by_region_year_with_oultiers_",model_name[i],"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  p<- bwplot(rmse~pred_mod | reg + year_predicted, data=tb,

             main="RMSE per model and region over all tiles")

  print(p)

  png(filename=paste("Figure8b_boxplot_overall_separated_by_region_year_scaling_",model_name[i],"_",out_suffix,".png",sep=""),

      width=col_mfrow*res_pix,height=row_mfrow*res_pix)

  p<- bwplot(rmse~pred_mod | reg + year_predicted, data=tb,ylim=c(0,5),

             main="RMSE per model and region over all tiles")

  print(p)