############################  Script for the experimentation of GAM model fitting ##############################

#The purpose is to test different parameters for the GAM function to allow fitting in area with low stations count.

##We experiment with:

# - gamma:  multiplying factor that acts on the smoothing parameter lambda

# - k: dimension for each term which affects the number of knots and complextiy of hte relationship

#

#The code generate a table of diagnostic that includes: k-index, rmse,mae, aic, gcv.

#

#MODIFIED ON: 08/07/2014            

#PROJECT: Environmental Layers project  NCEAS-NASA

# - check that step are increased by 1 for k

# - k starting  value can  be controled by increase or decrease

# - k initial value can be set

test_k_gam <- function(formula,l_k,data_training,l_k_min=NULL){

  if(is.null(l_k_min)){

    l_k_tmp <- as.numeric(l_k)/2  #if  no starting l_k then use half the l_k value (which correspond to the max)

  }

  if(!is.null(l_k_min)){

    l_k_tmp <- l_k_min  #if l_k then use l_k_min

  }

    #l_k_tmp <- l_k_tmp - rep(1,length(l_k_tmp)) #can modify this option to go -1 instead of plus 1?

fit_gam_model_with_diagnostics <- function(l_k,data_training,formula,names_mod,y_var_name,out_suffix,out_path,l_k_min=NULL){

  #formula: model formul without the k values!! eg y_var ~ s(lat,lon) + s(elev_s)

  #y_var_name: variable being model eg  dailyTmax

  #out_suffix: ouput suffix added toe the file  saved

  #out_path: output directory

  #TODO: 

  # -allow default val for out_suffix, y_var_name and out_path

  #### library and param

  ########## START  SCRIPT #############

  l_k_obj <- test_k_gam(formula,l_k,data_training,l_k_min)

  out_prefix  <- out_suffix

  filename_obj <- file.path(out_path,paste("diagnostics_obj_gam_fitting","_", y_var_name,out_prefix,".RData",sep=""))

  save(diagnostics_obj,file=filename_obj )

#First, copy data from NEX to Atlas:

##Second, load data from Northen Africa

in_dir <- "/data/project/layers/commons/NEX_data/output_regions/15.0_20.0" #Atlas location

raster_obj_infile <- "raster_prediction_obj_gam_CAI_dailyTmax15.0_20.0.RData" #Raster prediction object

out_path <- in_dir #output dir can  be set

out_suffix  <-"_08062014" #is the same as "out_prefix" in the workflow 

y_var_name <- "dailyTmax" #see workflow

setwd(out_path)

raster_obj <- load_obj(raster_obj_infile) #load raster predction object containing all the information

nb_models <- length((raster_obj$clim_method_mod_obj[[1]]$formulas)) #number of models predicted

list_formulas <- (raster_obj$clim_method_mod_obj[[1]]$formulas) #list of formula with models predicted

#Models used/predicted:

pred_mod <- paste("mod",c(1:nb_models,"_kr"),sep="") #names of models

nb_models <- length(pred_mod) #4 instead of 3

#Select one month to test and explored model fitting parameters:

clim_method_mod_obj <- raster_obj$clim_method_mod_obj[[j]] #object containing model fitted at the monthly time scale

#Object is a list made of "clim",data_month, data_month_v , sampling_month_dat, mod and formulas

names(clim_method_mod_obj)

clim_method_mod_obj$data_month #fitting  data

clim_method_mod_obj$data_month_v #no hold out for monthly data

l_mod <- clim_method_mod_obj$mod # R model object for GAM

data_training <- clim_method_mod_obj$data_month #no hold out for monthly data

plot(reg_rast,y=15) #plot SRTM (elev)

dim(reg_rast) #7,200 ; 3,600 ; 15

### Which model to test? Use model 2 because it is not fitted.

formula <-list_formulas[[k]] #model 2: y_var ~ s(lat, lon) + s(elev_s) + s(LST)

#### TRY with different gamma parameter value

mod_g1 <- try(gam(y_var ~ s(lat, lon) + s(elev_s),gamma=1.4, data=data_training)) #change to any model!!

mod_g3<- try(gam(y_var ~ s(lat, lon) + s(elev_s),gamma=100, data=data_training)) #change to any model!!

gam.check(mod_g3) #increase in k-index!!

#Now try with model2: y_var ~ s(lat, lon) + s(elev_s) + s(LST)

mod<- try(gam(y_var ~ s(lat, lon) + s(elev_s) + s(LST),gamma=10, data=data_training)) #does not fit!!

mod<- try(gam(y_var ~ s(lat, lon) + s(elev_s) + s(LST),gamma=100, data=data_training)) #does not fit!!

## check for residual pattern, removeable by increasing `k'

## typically `k', below, chould be substantially larger than 

## the original, `k' but certainly less than n/2.

### Visualization of predicted values

vis.gam(mod_t1)

vis.gam(mod_t1,view=c("lat","lon"),theta= 35) # plot against by variable

x_data <- mod_t1$model #model matrix also known as design matrix, note that it is stored as a data.frame

plot(mod_t1)

y<- fitted.values(mod_t1)

x<- x_data$lat

plot(y ~ x) #this is predicted versus lat

points(x_data$y_var ~ x,col="red")

y<- fitted.values(mod_t1)

x<- x_data$LST

plot(y~x) #this is predicted vs LST

names_mod <- paste("mod",k,sep="")

l_k <- c(30,10,10) #default values : max val

l_k_min <- c(10,4,4) #minimum values, default is null

l_k_obj <- test_k_gam(formula,l_k,data_training,l_k_min)

list_mod <- l_k_obj$list_mod 

list_mod<- list_mod[unlist(lapply(list_mod,FUN=function(x){!inherits(x,"try-error")}))]

k <- 2 #select model 3 with LST

data_training <- clim_method_mod_obj$data_month #this is for July

l_k_min <- c(10,4,4) #minimum values, default is null

#y_var_name, out_suffix, out_path are all set at the beginning of the script

#test_obj <- fit_gam_model_with_diagnostics(l_k,data_training,formula,model_name)

out_suffix_s  <- paste(7,"_",out_suffix,sep="") #7 for month of july

#could transform the input arguments to be a list??

#use l_k_min=NULL to start with half of the maximum k values...

test1_gam_fit_obj <- fit_gam_model_with_diagnostics(l_k,data_training,formula,names_mod,y_var_name,out_suffix_s,out_path,l_k_min=NULL)

#use l_k_min=c(10,4,4) 

out_suffix_s  <- paste(7,"_","lk_min_",out_suffix,sep="") #7 for month of july

test2_gam_fit_obj <- fit_gam_model_with_diagnostics(l_k,data_training,formula,names_mod,y_var_name,out_suffix_s,out_path,l_k_min)

names(test2_gam_fit_obj)

df_diagnostics <- test2_gam_fit_obj$df_diagnostics #extract the diagnostic  table

names(df_diagnostics)

df_diagnostics$fit_no

df_diagnostics$rmse

plot(df_diagnostics$rmse)

#plot(df_diagnostics$rmse~df_diagnostics$fit_no,data=df_diagnostics)

#Now do this over 12 months for each model??

#loop here or lapply