# This R environment comes with all of CRAN preinstalled, as well as many other helpful packages
# The environment is defined by the kaggle/rstats docker image: https://github.com/kaggle/docker-rstats

library(h2o) # R API wrapper for h2o

# Input data files are available in the "../input/" directory.
# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory
#system("ls ../input")

# Any results you write to the current directory are saved as output.

# Note: the skeleton of workflow/submission preparation reused from 
# Mark Laundry's (https://www.kaggle.com/mlandry) script demostating H2O GBM application 
# to the current competition problem.
# George's contribution was to specifically build a framework to apply and validate H2O deep learning model 
# for this contest.
# The currently presented model does not perform well, however it could be a starting point for 
# future parameter tune-up and other advanced modelling techniques.
# Mark's original script can be reviewed at at https://www.kaggle.com/mlandry/grupo-bimbo-inventory-demand/h2o-gbm/


#################
## Set up Cluster (H2O is a Java ML Platform, with R/Python/Web/Java APIs)
#################
print(paste("Set up Cluster",Sys.time()))

## start a cluster; default is 2 cores for CRAN reasons; -1 = all
h2o.init(nthreads=-1,max_mem_size='7G')  

#############
## Load Data 
#############
print(paste("Load Data",Sys.time()))
## load the full training file, using all cores; load into H2O Frame known to R as train; H2O as train.hex
train<-h2o.uploadFile("..input/train.csv",destination_frame = "train.hex")

train[1:2,] ## take a look at a few rows
## train on the log of the target
train$target<-log(train$Demanda_uni_equil+1)

#################
## Partition Data
#################
print(paste("Partition Data",Sys.time()))
## This model will use three splits, partitioned by week: 
##   one to generate product averages, a second to fit a model, and a third to evaluate the model
dev<-train[train$Semana <= 5,]                    ## gets Semana 3,4,5
val<-train[train$Semana > 5 & train$Semana < 8,]  ## gets Semana 6, 7
val[1:5,]
final<-train[train$Semana == 8,]                  ## gets Semana 8
final[1:5,]

##############################
## Model: Product Groups & DL
##############################
print(paste("Model: Product Groups & DL",Sys.time()))
## use the fields available in test to predict; removing ID and Semana
##   for iteration, Semana would probably be a good field to use to control loops
##   in H2O you can directly use it as a "fold column" if you'd like
predictors<-c("Agencia_ID","Canal_ID","Ruta_SAK","Cliente_ID","Producto_ID")

## first part of model: use product averages, created on the dev set
##  this is the only time we will use the dev set
groups<-h2o.group_by(data=dev,by="Producto_ID",mean("target"))
groups[1:5,]

## apply groups back into dev and validation data sets as "mean_target"
## if there are NAs for this (new products), use a constant; used median of entire train target
newVal<-h2o.merge(x=val,y=groups,all.x = T)
newVal$mean_target[is.na(newVal$mean_target)]<-0.7
newVal[1:5,]
newFinal<-h2o.merge(x=final,y=groups,all.x = T)
newFinal$mean_target[is.na(newFinal$mean_target)]<-0.7
newFinal[1:5,]


## train a DL; use aggressive parameters to keep overall runtime within 20 minutes
## this model is fit on Semana 6 & 7, and evaluated on Semana 8.

hyper_params <- list(
  activation=c("Rectifier","Tanh","Maxout","RectifierWithDropout","TanhWithDropout","MaxoutWithDropout"),
  hidden=list(c(20,20),c(50,50),c(30,30,30),c(25,25,25,25)),
  input_dropout_ratio=c(0.04,0.05),
  l1=seq(0,1e-4,1e-6),
  l2=seq(0,1e-4,1e-6)
)
hyper_params

## Stop once the top 5 models are within 1% of each other (i.e., the windowed average varies less than 1%)
search_criteria <- list(strategy = "RandomDiscrete", 
           max_runtime_secs = 360, max_models = 100, seed=1234567, 
           stopping_rounds=5, stopping_tolerance=1e-2)

dl_random_grid <- h2o.grid(
  algorithm="deeplearning",
  grid_id="dl_grid_random",       ## internal H2O name for the grid
  training_frame=newVal,          ## H2O frame holding the training data
  validation_frame=newFinal,      ## extra holdout piece for three layer modeling
  x=predictors,                   ## this can be names or column numbers
  y="target",                     ## target: using the logged variable created earlier
  epochs=5,
  stopping_metric="MSE",
  stopping_tolerance=1e-2,        ## stop when mse does not improve by >=1% for 5 scoring events
  stopping_rounds=5,
  score_validation_samples=10000, ## downsample validation set for faster scoring
  score_duty_cycle=0.025,         ## don't score more than 2.5% of the wall time
  max_w2=10,                      ## can help improve stability for Rectifier
  hyper_params = hyper_params,
  search_criteria = search_criteria,
  nfolds=5,                       ## nfolds CV
  fold_assignment="Modulo",       # can be "AUTO", "Modulo", "Random" or "Stratified"
  offset_column = "mean_target"
)
dl_random_grid
                                
grid <- dl_random_grid

grid@summary_table[1,]
best_model <- h2o.getModel(grid@model_ids[[1]]) ## model with lowest MSE 

## look at model diagnostics
summary(best_model)
g <- best_model

# clean up frames no longer needed
h2o.rm(train)
h2o.rm(dev)
h2o.rm(val)
h2o.rm(newVal)

#####################
## Create Predictions
#####################
print(paste("Create Predictions",Sys.time()))
## load test file
test<-h2o.uploadFile("..input/test.csv",destination_frame = "test.hex")
test[1:2,] ## take a look at a few rows of the test data
## merge in the offset column, just as with val and final
newTest<-h2o.merge(x=test,y=groups,all.x = T)
newTest$mean_target[is.na(newTest$mean_target)]<-0.7
newTest[1:5,]
p<-h2o.predict(g,newTest)
p<-exp(p)-1
summary(p)

####################
## Create Submission
####################
print(paste("Create Submission",Sys.time()))
submissionFrame<-h2o.cbind(test$id,p)
colnames(submissionFrame)<-c("id","Demanda_uni_equil")