# Required packages

rm(list=ls())
gc(TRUE)
Sys.getpid()

library(mgcv)
library(Matrix)
library(MASS)
library(BayesX)


# Initialize some characteristics of the model
min_year     <- 1980
fixed_s      <- FALSE
expectile1   <- 0.90
isotropy     <- "GMRF"
delta        <- 1e-4
reltol       <- 1e-5
initial_w    <- TRUE
step_max     <- 10
fitting_type <- "gcv"
dev_schall   <- "gcv"
bounce_off   <- TRUE

# Load expectile regression based on the mgcv package
source("function_expectreg_mgcv_bam11.R")

# Define name for saving
fixed_s_text <- "fixed"
if(!fixed_s) {
    fixed_s_text <- "optim"
}


delta_text <- as.character(delta)
delta_text <- gsub(x=delta_text,pattern=".",replacement="_",fixed=T)
delta_text <- paste(delta_text,step_max,sep="_")
reltol_text <- as.character(reltol)
reltol_text <- gsub(x=reltol_text,pattern=".",replacement="_",fixed=T)
reltol_text <- paste(reltol_text,step_max,sep="_")

expect_text <- as.character(expectile1*100)
initial_w_text <- "T"
if(!initial_w) initial_w_text <- "F"


NameSave <- paste("SpatTemp3d",min_year,isotropy,fixed_s_text,
                  expect_text,fitting_type,dev_schall,delta_text,reltol_text,initial_w_text,
                  ".RData",sep="_")

# Load data set
homstart <- read.table("Table_Temperatur_new3d.txt",sep=";",dec=".",header=T)

# Load Regions map
BND <- read.bnd("DEU_ROR.bnd")
neighbors <- read.gra("DEU_ROR.gra")

# Data cleaning
homstart <- subset(homstart, year >= min_year)
homstart <- subset(homstart, elevation <= 1850)
homstart$raw[homstart$raw < -100] <- NA
homstart$y <- homstart$raw

homstart <- homstart[!is.na(homstart$y),]
homstart$ROR <- as.factor(homstart$ROR)

homstart <- homstart[order(homstart$STATIONS_ID,homstart$MESS_DATUM),]

homstart <- homstart[,c("y","elevation","day","lon","lat","year","ROR")]

xt2 <- list()
xt2[[2]] <- list(polys=BND,penalty=neighbors)
xt1 <- list(polys=BND)

# Define formula
if(isotropy == "aniso") {
    f_exp <- y ~ s(elevation, bs="ps", k=7) + s(year, bs="ps", k=15, pc=1980) + 
        ti(day, bs="cp",k=15, np=F) + 
        ti(lon, bs="ps",k=6, np=F) + ti(lat, bs="ps",k=9, np=F) + 
        ti(lon,lat, bs=c("ps","ps"),d=c(1,1),k=c(6,9), np=F) + 
        ti(lon,day, bs=c("ps","cp"),d=c(1,1),k=c(6,15), np=F) + 
        ti(day,lat, bs=c("cp","ps"),d=c(1,1),k=c(15,9), np=F) +
        ti(day,lon,lat, bs=c("cp","ps","ps"),d=c(1,1,1),k=c(15,6,9), np=F)
}
if(isotropy == "GMRF") {
    f_exp <- y ~ s(elevation, bs="ps", k=7) + s(year, bs="ps", k=15, pc=1980) + 
        ti(day, bs="cp",k=15, np=F) + 
        ti(ROR, bs="mrf",xt=xt1,k=95, np=F)  +
        ti(day,ROR, bs=c("cp","mrf"),d=c(1,1),k=c(15,95), np=F, xt=xt2)
}

# Run intial model
# Note: The computational time is approximatly 1h
Sys.time()
set.seed(111)
rainmodel1 <- bam(f_exp, data = homstart, family = gaussian(), 
                  weight=rep(0.5,times=nrow(homstart)),method="GCV.Cp")
print(object.size(rainmodel1),units="Mb")
Sys.time()
gc(TRUE)
save.image(NameSave)
Sys.time()

# Extract initial smoothing parameter
sm_par_vec2 <- rainmodel1$sp

print(sm_par_vec2)

Sys.time()

quietly <- FALSE

# Run desired model
# Warning: The computational time is several days!!!
Sys.time()
gc(TRUE)
save.image(NameSave)
set.seed(111)
rainmodel2 <- expectreg_bam_smooth(formula=f_exp, data=homstart, expectiles=expectile1,
                                   sm_par_vec = sm_par_vec2, fixed=fixed_s,
                                   delta=delta, step_max=step_max, reltol = reltol,
                                   initial_w = initial_w, opt_type = fitting_type,
                                   dev_schall = dev_schall, quietly=quietly,
                                   bounce_off=bounce_off)

# Save workspace and quit process
Sys.time()
gc(TRUE)
save.image(NameSave)

if(.Platform$OS.type == "unix"){
    q("no")
}

