# Data dictionary: https://support.google.com/analytics/answer/3437719?hl=en
# Google Merchandise Store: https://shop.googlemerchandisestore.com/?utm_source=Partners&utm_medium=affiliate&utm_campaign=Data%20Share%20Promo
# GA Demo Account: https://support.google.com/analytics/answer/6367342?hl=en
# GA Conversion Probability Report: https://support.google.com/analytics/answer/7551092?hl=en
#
# train_v2.csv - the updated training set - contains user transactions from August 1st 2016 to April 30th 2018.
# test_v2.csv - the updated test set - contains user transactions from May 1st 2018 to October 15th 2018 (5.5 months)
# predict for all users in test what their fu_logRevenue is December 1st 2018 to January 31st 2019

library(tidyverse)
library(glue)
library(data.table)
library(feather)
library(lubridate)
library(caret)
library(catboost)

CACHE_FILE = "fm.feather"
SUBMISSION = T
SKIP_PREP = F & file.exists(CACHE_FILE) & !SUBMISSION
REBALANCE = T & !SUBMISSION # Part of prep

# Time interval parameters and logic --------------------------------------
data_origin <- ymd("2016-08-01") #...to ymd("2018-10-15")
cut_offs <- ymd("2017-01-15") %m+% months(0:17)  # earliest we can start and have 5.5 months data (to make a consistent recent visitor pool)
train_cut_offs <- cut_offs[cut_offs <= max(cut_offs) %m-% months(2)] 
valid_cut_offs <- train_cut_offs %m+% months(2)
recent    <- function(cut_off) interval(cut_off %m-% months(5) - days(14), cut_off)
target    <- function(cut_off) interval(cut_off %m+% months(2) - days(14), cut_off %m+% months(4) - days(15))
expanding <- function(cut_off) interval(data_origin, cut_off)
stopifnot(int_end(target(train_cut_offs)) <= int_start(target(valid_cut_offs)) - days(1))
stopifnot(int_start(recent(min(cut_offs))) >= data_origin) # need enough history for all cut_offs, to create the pool of recent visitors of the last 5.5 months

if ( SKIP_PREP  ) { #$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  print("Skipping prep...")
  fm <- setDT(read_feather(CACHE_FILE), key = "fullVisitorId")
  print("Read cache.")
} else {

# Load data ---------------------------------------------------------------
print("Starting prep...")
dt <- setDT(read_feather("../input/gstore-1-flatten/df.feather"), key = "fullVisitorId")

# Rebalance classes -------------------------------------------------------
if ( REBALANCE ) {
  # Doing this before cleaning, saves some time
  purchasers <- dt[!is.na(transactionRevenue) & transactionRevenue > 0, .(dummy = 1), by = "fullVisitorId"][, dummy := NULL]
  set.seed(2018)
  non_purchasers <- dt[!purchasers, .(dummy = 1), by = "fullVisitorId"][, dummy := NULL] %>% sample_frac(0.013)
  visitors_subset <- bind_rows(purchasers, non_purchasers)
  dt <- setkey(dt[visitors_subset], "fullVisitorId")
  print("Rebalanced classes.")
} 

# Clean data --------------------------------------------------------------
dt[, date := as.Date(date)]
dt[is.na(pageviews), pageviews := 1] 
dt[is.na(transactionRevenue), transactionRevenue := 0]
dt[is.na(referralPath), referralPath := "NA"]
dt[is.na(campaignCode), campaignCode := "NA"]
dt[, hour := month(visitStartTime)]
dt[is.na(city), city := "NA"]

# Function to build a single feature window -------------------------------
calc_window_stats <- function(t) {
  
  # t = min(cut_offs) # For code walk through
  recent_visitors <- dt[date %within% recent(t), .(fullVisitorId)] %>% unique
  rv <- dt[recent_visitors] # recent visitors (for train) but ALL their data
  
  # Calculate target (then discard future data)
  future <- rv[date %within% target(t), .(fu_logRevenue = log1p(sum(transactionRevenue))), by = "fullVisitorId"][recent_visitors]
  rv <- subset(rv, date %within% expanding(t))

  # Summarise transactions per visitor
  tx <-
    rv[transactionRevenue > 0, 
       .(count_transactions = .N, total_revenue = sum(transactionRevenue), # Removed later
         first_date = min(date), last_date = max(date)),
         by = "fullVisitorId"
     ][, `:=`(last_transaction_scaled = as.numeric(t - last_date) / as.numeric(t - first_date), # Note, since first TRANSACTION, not visit
              transactions_per_day = (count_transactions - 1) / as.numeric(t - first_date), # Repeat transactions since first
              revenue_per_day = total_revenue / as.numeric(t - first_date),
              count_transactions = NULL, total_revenue = NULL, first_date = NULL, last_date = NULL)
     ][recent_visitors]
  # BTYD sufficient statistics are (first, last, count), plus optional monetary value (sum_tr)
  # We are using a subset of these (total_revenue/first, last/first, count_transactions/first), by incorporating "first"
  # I suspect this is estimating the rate parameter, lambda, of the poisson process: the 
  # mean number of repeat transactions since the first, under H0 where the customer is still "alive" at the cut off

  # Mean encode categoricals (at visit level)
  rv[, cr_operatingSystem := mean(transactionRevenue > 0), by = "operatingSystem"]
  rv[, cr_referralPath := mean(transactionRevenue > 0), by = "referralPath"]
  rv[, cr_medium := mean(transactionRevenue > 0), by = "medium"]
  rv[, cr_networkDomain := mean(transactionRevenue > 0), by = "networkDomain"]
  rv[, cr_campaignCode := mean(transactionRevenue > 0), by = "campaignCode"]
  rv[, cr_subContinent := mean(transactionRevenue > 0), by = "subContinent"]
  rv[, cr_country := mean(transactionRevenue > 0), by = "country"]
  rv[, cr_city := mean(transactionRevenue > 0), by = "city"]

  # Summarise past per visit
  past <- 
    rv[,
       .(first_date = min(date), last_date = max(date), # Removed later
         sum_pageviews = sum(pageviews),
         mean_pageviews = mean(pageviews),
         max_visitNumber = max(visitNumber), # very similar to count_visits
         cr_operatingSystem = median(cr_operatingSystem),
         cr_referralPath = median(cr_referralPath),
         cr_medium = median(cr_medium),
         cr_networkDomain = median(cr_networkDomain),
         cr_campaignCode = median(cr_campaignCode),
         cr_city = median(cr_city),
         cr_subContinent = median(cr_subContinent),
         cr_country = median(cr_country)),
         by = "fullVisitorId"
     ][, `:=`(first_visit = as.integer(t - first_date), # Done in two steps for speed
              last_visit = as.integer(t - last_date),
              first_date = NULL, last_date = NULL)
     ][recent_visitors] # Belt and braces, as already guaranteed one row per recent_visitors
  #BTYD sufficient statistics are (first_visit, last_visit, max_visitNumber), and "monetary" value (e.g. sum_pageviews & mean_pageviews)

  future[tx][past][is.na(fu_logRevenue), fu_logRevenue := 0][, cut_off := t]
}

# Build all windows
fm <- setDT(map_dfr(cut_offs, calc_window_stats), key = "fullVisitorId")

print("Finished prep.")
write_feather(fm, "fm.feather")
setDF(fm); write_rds(fm, "fm.rds", compress = "gz"); setDT(fm);
print("Saved cache.")
} #$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

# Set up data for CV -------------------------------------------------------
train <- lapply(train_cut_offs, function(t) which(fm$cut_off %within% interval(t %m-% months(12), t))) 
valid <- lapply(valid_cut_offs, function(t) which(fm$cut_off == t))
final <- which(fm$cut_off %within% interval(max(cut_offs) %m-% months(12), max(cut_offs)))

ignore_columns = c("cut_off", "fu_logRevenue", "fullVisitorId")
x <- setDF(copy(fm)) %>% select(-one_of(ignore_columns)) %>% mutate_if(is.character, as.factor)
y <- fm$fu_logRevenue

# Baseline model -----------------------------------------------------------
print("Baseline #1: always predict zero")
print(RMSE(0, fm$fu_logRevenue))

# CatBoost -----------------------------------------------------------------
grid <- expand.grid(depth = 6, #c(5,6,7),
                    learning_rate = 0.05, #c(0.2,0.1,0.05),
                    iterations = 100, #c(50,100,200)
                    l2_leaf_reg = 3,
                    rsm = 1.0, #c(0.8,1.0),
                    border_count = 64) #c(64,128,255))
report <-
  train(x, y, 
        method = catboost.caret,
        preProc = NULL,
        tuneGrid = grid, #tuneLength = 10, #limit the random search
        trControl = trainControl(index = train, indexOut = valid, indexFinal = final), #search = "random"
        # Params for CatBoost:
        logging_level = "Silent", #verbose = 10,
        allow_writing_files = FALSE,
        random_seed = 42,
        task_type = "CPU",
        thread_count = 4)

print("CV scores:")
print(glue("RMSE {format(report$results$RMSE, digits = 3)} +/- {format(report$results$RMSESD, digits = 3)}"))
print(glue("Rsquared {format(report$results$Rsquared, digits = 3)} +/- {format(report$results$RsquaredSD, digits = 3)}"))

scores <-
  report$resample %>% 
  mutate(Resample = NULL, MAE = NULL,
         `Training Cut Off` = valid_cut_offs,
         `Target Interval End` = int_end(target(`Training Cut Off`)))
#print(scores)
qplot(x = `Target Interval End`, y = Rsquared, data = scores, geom = "point", label = as.character(`Target Interval End`)) + ggrepel::geom_text_repel() 

fi <- catboost.get_feature_importance(report$finalModel, fstr_type = "FeatureImportance")
print(tibble(name = names(x), importance = round(fi[, 1], 2)) %>% arrange(desc(importance)) %>% as.data.frame)

#fin <- catboost.get_feature_importance(report$finalModel, fstr_type = "Interaction")
#print(tibble(f1 = names(x)[fin[, 1]+1], f2 = names(x)[fin[, 2]+1], importance = fin[, 3]) %>% arrange(-importance))

# Submission ---------------------------------------------------------------
if ( SUBMISSION ) {
  final_cut_off <- ymd("2018-10-15")
  final_fm <- setDT(calc_window_stats(final_cut_off), key = "fullVisitorId")
  y_test <- predict(report, newdata = final_fm)
  y_test <- pmax(y_test, 0) # some negative? not sure how this is possible with a tree...
  submit <- data.frame(fullVisitorId = final_fm$fullVisitorId, PredictedLogRevenue = y_test)
  stopifnot(length(unique(dt[(test), ]$fullVisitorId)) == nrow(submit))
  write_csv(submit, "submission.csv")
}