#MODIFIED ON: 04/06/2014            

function_analyses_paper2 <-"multi_timescales_paper_interpolation_functions_03182014.R"

#Table 2: Covariates models and functional form (not produced in R)

#Table 3: Combinations of models tested (not produced in R)

#Table 4. Average RMSE for all models acroos year 2010 wiht 30% daily hold out

##### Table 4: Contribution of covariates using validation accuracy metrics

table4_paper <- as.data.frame(do.call(rbind,list(table_gam,table_kriging,table_gwr)))    

table4_paper <- round(table4_paper,digit=3) #roundto three digits teh differences

table4_paper$Methods <- c(rep("gam",7),

table4_paper$Formulas<-rep(list_formulas,3)                             

table4_paper<-table4_paper[(c(5,4,1,2,3))]  #reordering columns                           

###Table 4, writeout

names_table_col <-c("Models","Method","Single time scale","Multiple time scale CAI","Multiple time scale FSS") # 

names(table4_paper)<- names_table_col

print(table4_paper) #show resulting table

#Now write out table 4

file_name<-paste("table4_multi_timescale_paper","_",out_prefix,".txt",sep="")

write.table(table4_paper,file=file_name,sep=",")

################ Table 6 ################

table6 <- do.call(rbind,list_table_method_holdout)

table6_paper <- table6[,c(9,1,2,3,4,5,6,7,8)]

file_name<-paste("table6_multi_timescale_paper","_",out_prefix,".txt",sep="")

write.table(table6_paper,file=file_name,sep=",")

            main="RME by monthly hold-out proportions and  interpolation methods",

##Do analyses using density of stations per holdout

##DO lm modeling with factor  variable being hold out proportion

tb_mv <- do.call(rbind.fill,tb_mv_list) #[c("gam_CAI","kriging_CAI","gwr_CAI")])

tb_mv <- subset(tb_mv,tb_mv$method_interp %in%c("gam_CAI","kriging_CAI","gwr_CAI") )

#names(tb_mv)

#reg_area_r <- rasterize(interp_area,elev,field=1,fun=mean)

#tot_area <-freq(reg_area_r)[1,2]*prod(res(reg_area_r)/1000)

#table(tb_mv$pred_mod)

#unique(tb_mv$n) #note that

#station_density_prop <- unique(tb_mv$n)/tot_area

#station_density_prop[1:7]*1000

extract_table_term_factor <- function(i,list_param){

  #This function generate a linear model for proportion of hold out effect on accuracy

  #Add option to choose MAE later

  #First parse arguments

  interpolation_method <- list_param$list_interp_method[[i]]

  tb_mv <- list_param$tb_mv

  #Begin script:

  list_pred_mod_name <- unique(tb_mv$pred_mod)

  list_prop_cat <- unique(tb_mv$prop_month)  

  list_pred_mod_name <- grep(c("mod_kr"),list_pred_mod_name,value=TRUE,invert=TRUE)

  list_mod_table <- vector("list",length=length(list_pred_mod_name))

  tb_dat<- subset(tb_mv,tb_mv$method_interp==interpolation_method) 

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

      mod <-lm(rmse~ as.factor(prop_month),

               data=subset(tb_dat,tb_dat$pred_mod==list_pred_mod_name[[j]]))      

      term_table <- summary(mod)$coefficients

      list_mod_table[[j]] <- term_table  

  }

  names(list_mod_table) <- list_pred_mod_name

  tx <- lapply(list_mod_table,function(x){x[,c(1,4)]})

  tx <-lapply(tx, round,digit=3)

  #tx <- format(tx,digits=3)

  #column_tab <- paste(format(tx[[1]][,1],digits=3)," ",

  #             "(",format(tx[[1]][,2],digits=3),")",sep="")

  #column_tab <-lapply(tx,function(x){paste(format(x[,1],digits=3)," ",

  #             "(",format(x[,2],digits=3),")",sep="")})

  column_tab <-lapply(tx,function(x){as.data.frame(paste(format(x[,1],digits=3)," ",

               "(",format(x[,2],digits=3),")",sep=""))})

  table_method <- as.data.frame(do.call(cbindX,column_tab))

  names(table_method) <- list_pred_mod_name

  table_method$method_interp <- rep(interpolation_method,nrow(table_method))

  table_method$prop <- list_prop_cat

  mod_table_obj<- list(list_mod_table,table_method)

  names(mod_table_obj) <- c("list_mod_table","table_method")

  return(mod_table_obj)

}  

  table4_paper <- as.data.frame(do.call(rbind,list(table_gam,table_kriging,table_gwr)))    

table4_paper <- round(table4_paper,digit=3) #roundto three digits teh differences

table4_paper$Methods <- c(rep("gam",7),

                          rep("kriging",7),

                          rep("gwr",7))    

}

#(res(reg_area_r))

#### DATA MONTH: correlation between Tmax and LST comparison...

plot(1:12,x_cor_tmp1[[1]][1:12,13],ylim=c(-0.8,0.5),

### NOW EXMAINE LST-TMAx and elev cor

plot(TMax~mm_01,data=subset(data_month_all,data_month_all$month==1))

plot(TMax~mm_02,data=subset(data_month_all,data_month_all$month==1))

plot(TMax~mm_06,data=subset(data_month_all,data_month_all$month==6))

plot(TMax~mm_07,data=subset(data_month_all,data_month_all$month==7))

cor(data_month)

names_tmp<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08",

             "mm_09","mm_10","mm_11","mm_12")#,"LST","elev_s")

list_cor_x <- vector("list",length=12)

for(i in 1:12){

  data_month<-list_data_month_s[[i]]

  data_x <- data_month[,c("TMax","elev_s",names_tmp[i],"LC1","LC6","LC9")]

  data_x <- as.data.frame(data_x)[,1:3]

  cor_x<- cor(data_x,use="complete.obs")

  list_cor_x[[i]] <- cor_x

}

cor_tab<-lapply(list_cor_x,function(x){x[1,2:5]})

cor_tab<-lapply(list_cor_x,function(x){x[1,2:4]})

cor_tab<-do.call(rbind,cor_tab)

#cor()

data

#plot(TMax~mm_01,data_month_all)

#         for(i in 1:12){

#  data_month<-list_data_month_s[[i]]

#  data_x <- data_month[,c("TMax","elev_s",names_tmp[i],"LC1","LC6","LC9")]

#  data_x <- as.data.frame(data_x)[,1:3]

#  cor_x<- cor(data_x,use="complete.obs")

#  list_cor_x[[i]] <- cor_x

#}