#This extract a data.frame object from raster prediction obj and combine them in one data.frame 

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

}

calc_stat_prop_tb <-function(names_mod,raster_prediction_obj,testing=TRUE){

  #add for testing??

  if (testing==TRUE){

    tb <-raster_prediction_obj$tb_diagnostic_v #use testing accuracy information

  }else{

    tb <-raster_prediction_obj$tb_diagnostic_s #use training accuracy information

  }

  t<-melt(subset(tb,pred_mod==names_mod),

          measure=c("mae","rmse","r","m50"), 

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

          na.rm=T)

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

  sd_tb<-cast(t,pred_mod+prop~variable,sd)

  n_tb<-cast(t,pred_mod+prop~variable,length)

  #n_NA<-cast(t,dst_cat1~variable,is.na)

  #### prepare returning object

  prop_obj<-list(tb,avg_tb,sd_tb,n_tb)

  names(prop_obj) <-c("tb","avg_tb","sd_tb","n_tb")

  return(prop_obj)

}

#ploting 

plot_prop_metrics <-function(list_param){

  #

  #list_dist_obj: list of dist object 

  #col_t: list of color for each 

  #pch_t: symbol for line

  #legend_text: text for line and symbol

  #mod_name: selected models

  #

  ## BEGIN ##

  list_obj<-list_param$list_prop_obj

  col_t <-list_param$col_t 

  pch_t <- list_param$pch_t 

  legend_text <- list_param$legend_text

  list_mod_name<-list_param$mod_name

  metric_name<-list_param$metric_name

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

    l<-list_obj[[i]]

    mod_name<-list_mod_name[i]

    avg_tb<-subset(l$avg_tb,pred_mod==mod_name,select=metric_name) #selecte relevant accuarcy metric

    n_tb<-subset(l$n_tb,pred_mod==mod_name,select=metric_name) 

    sd_tb<-subset(l$sd_tb,pred_mod==mod_name,select=metric_name) #l$sd_abs_tb[,metric_name]

    xlab_text<-"holdout proportion"

    no <- unlist(as.data.frame(n_tb))

    y <- unlist(as.data.frame(avg_tb))

    x<- l$avg_tb$prop

    y_sd <- unlist(as.data.frame(sd_tb)) #sd_tb

    ciw <-y_sd

    #ciw2   <- qt(0.975, n) * y_sd2 / sqrt(n)

    #plotCI(y=y, x=x, uiw=ciw, col="red", main=paste(" MAE for ",mod_name,sep=""), barcol="blue", lwd=1,

    #       ylab="RMSE (C)", xlab=xlab_text)

    ciw   <- qt(0.975, no) * y_sd / sqrt(no)

    if(i==1){

      plotCI(y=y, x=x, uiw=ciw, col=col_t[i], main=paste(" Comparison of ",metric_name," in ",mod_name,sep=""), barcol="blue", lwd=1,

             ylab="RMSE (C)", xlab=xlab_text)

      lines(y~x, col=col_t[i])

    }else{

      lines(y~x, col=col_t[i])

    }

  }

  legend("topleft",legend=legend_text, 

         cex=1.2, pch=pch_t,col=col_t,lty=1,bty="n")

  #axis(1,at=1:length(stat_tb[,1]),labels=stat_tb[,1])

}

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

in_dir5<- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_03142013/output_data_365d_gam_daily_mults10_lst_comb3_08082013"

raster_obj_file_5 <- "raster_prediction_obj_gam_daily_dailyTmax_365d_gam_daily_mults10_lst_comb3_08082013.RData"

raster_obj_file_6 <- "raster_prediction_obj_kriging_daily_dailyTmax_365d_kriging_daily_mults10_lst_comb3_08062013.RData"

raster_obj_file_7 <- "raster_prediction_obj_gwr_daily_dailyTmax_365d_gwr_daily_mults10_lst_comb3_08072013.RData"

#raster_prediction_obj_gam_daily_dailyTmax_365d_gam_daily_mults10_lst_comb3_08082013.RData

raster_prediction_obj_7 <-load_obj(file.path(in_dir7,raster_obj_file_7)) #kriging daily multisampling 

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

####### Now create figures #############

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

################### Figure 3

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

#########Figure 4. RMSE and MAE, mulitisampling and hold out for FSS and GAM.

#names_mod<-c("mod1","mod2","mod3","mod4","mod5","mod6","mod7","mod8","mod9")

names_mod<-c("mod1")

debug(calc_stat_prop_tb)

prop_obj_gam<-calc_stat_prop_tb(names_mod,raster_prediction_obj_5)

prop_obj_kriging<-calc_stat_prop_tb(names_mod,raster_prediction_obj_6)

prop_obj_gwr<-calc_stat_prop_tb(names_mod,raster_prediction_obj_7)

list_prop_obj<-list(prop_obj_gam,prop_obj_kriging,prop_obj_gwr)

col_t<-c("red","blue","black")

pch_t<- 1:length(col_t)

legend_text <- c("GAM","Kriging","GWR")

mod_name<-c("mod1","mod1","mod1")#selected models

metric_name<-"rmse"

#png_names<- 

#png(file.path(out_path,paste("clim_surface_",y_var_name,"_",sampling_dat$date[i],"_",sampling_dat$prop[i],

#                             "_",sampling_dat$run_samp[i],out_prefix,".png", sep="")))

list_param_plot<-list(list_prop_obj,col_t,pch_t,legend_text,mod_name,metric_name)

names(list_param_plot)<-c("list_prop_obj","col_t","pch_t","legend_text","mod_name","metric_name")

#debug(plot_prop_metrics)

plot_prop_metrics(list_param_plot)

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

#########Figure 5. Overtraining tendency

tb7 <-raster_prediction_obj_7$tb_diagnostic_v  #gwr daily methods

prop_obj_gam_s<-calc_stat_prop_tb(names_mod,raster_prediction_o  bj_5,testing=FALSE)

prop_obj_kriging_s<-calc_stat_prop_tb(names_mod,raster_prediction_obj_6,testing=FALSE)

prop_obj_gwr_s<-calc_stat_prop_tb(names_mod,raster_prediction_obj_7,testing=FALSE)

prop_obj_gam$avg_tb - prop_obj_gam_s$avg_tb

plot(prop_obj_gam_s$avg_tb$rmse ~ prop_obj_gam_s$avg_tb$prop, type="b",)

y_range<-range(prop_obj_gam$avg_tb$rmse,prop_obj_gam_s$avg_tb$rmse)

plot(prop_obj_gam$avg_tb$rmse ~ prop_obj_gam$avg_tb$prop, type="b",ylim=y_range)

lines(prop_obj_gam_s$avg_tb$rmse ~ prop_obj_gam_s$avg_tb$prop, type="b",ylim=y_range,col=c("red"))

lines(prop_obj_gam$avg_tb$rmse ~ prop_obj_gam$avg_tb$prop, type="b",ylim=y_range,col=c("red"),lty=2)

lines(prop_obj_gwr_s$avg_tb$rmse ~ prop_obj_gwr_s$avg_tb$prop, type="b",ylim=y_range,col=c("black"))

lines(prop_obj_gwr$avg_tb$rmse ~ prop_obj_gam$avg_tb$prop, type="b",ylim=y_range,col=c("black"),lty=2)

y_range<-range(prop_obj_kriging$avg_tb$rmse,prop_obj_kriging_s$avg_tb$rmse)

plot(prop_obj_kriging$avg_tb$rmse ~ prop_obj_kriging$avg_tb$prop, type="b",ylim=y_range,col=c("blue"),lty=2)

lines(prop_obj_kriging_s$avg_tb$rmse ~ prop_obj_kriging_s$avg_tb$prop, type="b",ylim=y_range,col=c("blue"))

## Calculate average difference for RMSE for all three methods

#read in relevant data:

tb1_s<-extract_from_list_obj(raster_prediction_obj_1$validation_mod_obj,"metrics_s")

rownames(tb1_s)<-NULL #remove row names

tb1_s$method_interp <- "gam_daily" #add type of interpolation...out_prefix too??

tb3_s<-extract_from_list_obj(raster_prediction_obj_4$validation_mod_obj,"metrics_s")

rownames(tb1_s)<-NULL #remove row names

tb3_s$method_interp <- "kriging_daily" #add type of interpolation...out_prefix too??

tb4_s<-extract_from_list_obj(raster_prediction_obj_3$validation_mod_obj,"metrics_s")

rownames(tb4_s)<-NULL #remove row names

tb4_s$method_interp <- "gwr_daily" #add type of interpolation...out_prefix too??

#tb1_s <-raster_prediction_obj_1$tb_diagnostic_s  #gam dailycontains the accuracy metrics for each run...

#tb3_s <-raster_prediction_obj_3$tb_diagnostic_s  #Kriging daily methods

#tb4_s <-raster_prediction_obj_4$tb_diagnostic_s  #gwr daily methods

tb1 <-raster_prediction_obj_1$tb_diagnostic_v  #gam dailycontains the accuracy metrics for each run...

tb3 <-raster_prediction_obj_3$tb_diagnostic_v  #Kriging daily methods

tb4 <-raster_prediction_obj_4$tb_diagnostic_v  #gwr daily methods

diff_df<-function(tb_s,tb_v,list_metric_names){

  tb_diff<-vector("list", length(list_metric_names))

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

    metric_name<-list_metric_names[i]

    tb_diff[[i]] <-tb_s[,c(metric_name)] - tb_v[,c(metric_name)]

  }

  names(tb_diff)<-list_metric_names

  tb_diff<-as.data.frame(do.call(cbind,tb_diff))

  return(tb_diff)

#debug(diff_df)

diff_tb1 <-diff_df(tb1_s[tb1_s$pred_mod=="mod1",],tb1[tb1$pred_mod=="mod1",],c("mae","rmse")) #select differences for mod1

diff_tb3 <-diff_df(tb3_s[tb3_s$pred_mod=="mod1",],tb3[tb3$pred_mod=="mod1",],c("mae","rmse"))

diff_tb4 <-diff_df(tb4_s[tb4_s$pred_mod=="mod1",],tb4[tb4$pred_mod=="mod1",],c("mae","rmse"))

x<-data.frame(gam=diff_tb1$mae,gwr=diff_tb3$mae,kriging=diff_tb4$mae)

boxplot(x) #plot differences in training and testing accuracies for three methods

mae_tmp<- data.frame(gam=tb1[tb1$pred_mod=="mod1",c("mae")],

                     kriging=tb3[tb3$pred_mod=="mod1",c("mae")],

                     gwr=tb4[tb4$pred_mod=="mod1",c("mae")])

plot(mae_tmp$gam,col=c("red"),type="b")

lines(mae_tmp$kriging,col=c("blue"),type="b")

lines(mae_tmp$gwr,col=c("black"),type="b")

legend("topleft",legend=legend_text, 

       cex=1.2, pch=pch_t,col=col_t,lty=1,bty="n")

max(mae_tmp$gam)

x2<-tb1[tb1$pred_mod=="mod1",c("mae","date")]

arrange(x2,desc(mae))

pmax(x2$mae,x2$date)

kriging=tb3[tb3$pred_mod=="mod1",c("mae")],

                     gwr=tb4[tb4$pred_mod=="mod1",c("mae")])

##### MONTHLY AVERAGES

tb1_month<-raster_prediction_obj_1$summary_month_metrics_v[[1]] #note that this is for model1

tb3_month<-raster_prediction_obj_3$summary_month_metrics_v[[1]]

tb4_month<- raster_prediction_obj_4$summary_month_metrics_v[[1]]

y_range<-range(tb1_month$mae,tb3_month$mae,tb4_month$mae)

plot(1:12,tb1_month$mae,col=c("red"),type="b",ylim=y_range)

lines(1:12,tb3_month$mae,col=c("blue"),type="b")

lines(1:12,tb4_month$mae,col=c("black"),type="b")

date<-strptime(tb1$date, "%Y%m%d")   # interpolation date being processed

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

tb1$month<-month

x3<-tb1[tb1$pred_mod=="mod1",]

(plot(x3[month=="01",c("mae")]))

median(x3[x3$month=="03",c("mae")],na.rm=T)

mean(x3[x3$month=="03",c("mae")],na.rm=T)

plot_dst_spat_fun<-function(stat_tb,names_var,cat_val){

  range_y<-range(as.vector(unlist(stat_tb[,names_var])),na.rm=T) #flatten data.frame

  col_t<-rainbow(length(names_var))

  pch_t<- 1:length(names_var)

  plot(stat_tb[,names_var[1]], ylim=range_y,pch=pch_t[1],col=col_t[1],type="b",

       yla="MAE (in degree C)",xlab="",xaxt="n")

  #points((stat_tb[,names_var[k]], ylim=range_y,pch=pch_t[1],col=col_t[1]),type="p")

  for (k in 2:length(names_var)){

    lines(stat_tb[,names_var[k]], ylim=range_y,pch=pch_t[k],col=col_t[k],type="b",

          xlab="",axes=F)

    #points((stat_tb[,names_var[k]], ylim=range_y,pch=pch_t[k],col=col_t[k]),type="p")

  }

  legend("topleft",legend=names_var, 

         cex=1.2, pch=pch_t,col=col_t,lty=1,bty="n")

  axis(1,at=1:length(stat_tb[,1]),labels=stat_tb[,1])

}

plot_prop_metrics <-function(list_param){

  list_obj<-list_param$list_prop_obj

  col_t <-list_param$col_t 

  pch_t <- list_param$pch_t 

  metric_name<-list_param$metric_name

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

    l<-list_obj[[i]]

    avg_tb<-subset(l$avg_tb,pred_mod==mod_name,select=metric_name) #selecte relevant accuarcy metric

    n_tb<-subset(l$n_tb,pred_mod==mod_name,select=metric_name) 

    sd_tb<-subset(l$sd_tb,pred_mod==mod_name,select=metric_name) #l$sd_abs_tb[,metric_name]

    xlab_text<-"holdout proportion"

    no <- unlist(as.data.frame(n_tb))

    y <- unlist(as.data.frame(avg_tb))

    x<- l$avg_tb$prop

    y_sd <- unlist(as.data.frame(sd_tb)) #sd_tb

    ciw   <- qt(0.975, no) * y_sd / sqrt(no)

      plotCI(y=y, x=x, uiw=ciw, col=col_t[i], main=paste(" Comparison of ",metric_name," in ",mod_name,sep=""), barcol="blue", lwd=1,