library(data.table)
library(fastmatch)
library(zoo)
library(xgboost)

##########FUNCTIONS###########
#Fast %in%
`%fin%` <- function(x, lkup) {
  fmatch(x, lkup, nomatch = 0L) > 0L
}

#Get the time differences between each measure
time_difference <- function(times, num_per_segment = 60) {
  n <- length(times)
  valid_time <- vector(mode="numeric", length = n)
  valid_time[1] <- times[1]
  valid_time[-1] <- diff(times, 1)
  valid_time[n] <- valid_time[n] + num_per_segment - sum(valid_time)
  valid_time <- valid_time / num_per_segment
  valid_time
}

#Convert reflectivity (dbz) to mm/hr
marshall_palmer <- function(dbz) {
  ((10**(dbz/10))/200) ** 0.625
}



##Valid values based on 0.01in measurements
valid_vals <- 0.254 * 1:276

######Training Data############
tr_raw <- fread("../input/train.csv",  select = c(
    "Id", 
    "minutes_past", 
    "radardist_km",
    "Ref", 
    "Ref_5x5_50th",
    "Ref_5x5_90th",
    "RefComposite", 
    "RefComposite_5x5_50th",
    "RefComposite_5x5_90th",
    "Expected")
    )

#Cut off Ref values < 0
tr_raw$Ref_5x5_50th[which(tr_raw$Ref_5x5_50th < 0)] <- NA
tr_raw$Ref_5x5_90th[which(tr_raw$Ref_5x5_90th < 0)] <- NA
tr_raw$RefComposite[which(tr_raw$RefComposite < 0)] <- NA
tr_raw$RefComposite_5x5_50th[which(tr_raw$RefComposite_5x5_50th < 0)] <- NA
tr_raw$Ref[which(tr_raw$Ref < 0)] <- NA
tr_raw$RefComposite_5x5_90th[which(tr_raw$RefComposite_5x5_90th < 0)] <- NA

tr_raw <- tr_raw[round(Expected, 4) %fin% valid_vals]
tr_raw$dt <- time_difference(tr_raw$minutes_past)
tr_raw$mp <- marshall_palmer(tr_raw$Ref)

#cor(tr_raw,use = "pairwise.complete.obs")

#Collapse to one record per Id
tr <- tr_raw[, .(
    target = log1p(mean(Expected, na.rm = T)),
    ref = mean(dt * Ref, na.rm = T),
    ref5 = mean(dt * Ref_5x5_50th, na.rm = T),
    ref9 = mean(dt * Ref_5x5_90th, na.rm = T),
    ref1 = mean(dt * RefComposite, na.rm = T),
    refc5 = mean(dt * RefComposite_5x5_50th, na.rm = T),
    refc9 = mean(dt * RefComposite_5x5_90th, na.rm = T),
    mp = sum(dt * mp, na.rm = T),
    rd = mean(radardist_km, na.rm = T),
    records = .N,
    naCounts = sum(is.na(Ref))
    ), Id]

   
print("training model...")
cs <- c("ref", "ref1",   "mp", "rd", "records")
y<-tr$target
tr<-as.data.frame(tr)
tr<-tr[,cs]
param0 <- list("objective"  = "reg:linear" 
  , "eval_metric" = "rmse"
  , "eta" = 0.007
  , "subsample" = 0.7
  , "min_child_weight" =10    
  , "max_depth" = 8
  , "nthreads" = 4
)
xgtrain = xgb.DMatrix(as.matrix(tr), label = y, missing = NA)
rm(tr,tr_raw)
gc()

x.mod.t  <- xgb.train(params = param0, data = xgtrain , nrounds =1200)


#model <- xgb.dump(x.mod.t, with.stats = T)
#model[1:10]
 features <- dimnames(xgtrain)[[2]]
 
importanceMatrix <- xgb.importance(features, model = x.mod.t)
xgb.plot.importance(importanceMatrix)




print("Processing test data...")
te_raw<-fread("../input/test.csv", select=c(    
    "Id", 
    "minutes_past", 
    "radardist_km",
    "Ref",
    "Ref_5x5_50th",
    "Ref_5x5_90th",
    "RefComposite", 
    "RefComposite_5x5_50th",
    "RefComposite_5x5_90th",
    "Expected")
    )
    
    #Cut off Ref values < 0
te_raw$Ref_5x5_50th[which(te_raw$Ref_5x5_50th < 0)] <- NA
te_raw$Ref_5x5_90th[which(te_raw$Ref_5x5_90th < 0)] <- NA
te_raw$RefComposite[which(te_raw$RefComposite < 0)] <- NA
te_raw$RefComposite_5x5_50th[which(te_raw$RefComposite_5x5_50th < 0)] <- NA
te_raw$Ref[which(te_raw$Ref < 0)] <- NA
te_raw$RefComposite_5x5_90th[which(te_raw$RefComposite_5x5_90th < 0)] <- NA

te_raw$dt <- time_difference(te_raw$minutes_past)
te_raw$mp <- marshall_palmer(te_raw$Ref)
    
te <- te_raw[, .(
    ref = mean(dt * Ref, na.rm = T),
    ref5 = mean(dt * Ref_5x5_50th, na.rm = T),
    ref9 = mean(dt * Ref_5x5_90th, na.rm = T),
    ref1 = mean(dt * RefComposite, na.rm = T),
    refc5 = mean(dt * RefComposite_5x5_50th, na.rm = T),
    refc9 = mean(dt * RefComposite_5x5_90th, na.rm = T),
    mp = sum(dt * mp, na.rm = T),
    rd = mean(radardist_km),
    records = .N
    ),Id]
te<-as.data.frame(te)
  
xgtest = xgb.DMatrix(as.matrix(te[,cs ]), missing = NA)

pr  <- predict(x.mod.t,xgtest)
sample_sol <-fread("../input/sample_solution.csv")

xgb_prediction <- expm1(pr)

res <- data.frame(
  Id = te$Id,
  Expected =  xgb_prediction 
)

#convert expected values to 0.01in values
res$Expected <- round(res$Expected / 0.254) * 0.254

summary(res)

write.csv(res, "xgbpvnr218.csv", row.names = FALSE, col.names = TRUE)