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

##########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:300

######Training Data############
set.seed(1234)
tr_raw <- fread("../input/train.csv")


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)

#Collapse to one record per Id
tr <- tr_raw[, .(
        target = log1p(mean(Expected, na.rm = T)),
        ref = mean(Ref, na.rm = T),
        ref1 = mean(RefComposite, na.rm = T),
        mp = sum(mp, na.rm = T),
        rd = mean(radardist_km, na.rm = T),
        records = .N,
        naCounts = sum(is.na(Ref)),
        mean_RhoHV = mean(RhoHV,na.rm=T),
        meanRef_5x5_10th = mean(Ref_5x5_10th,na.rm=T),
        meanRef_5x5_50th = mean(Ref_5x5_50th,na.rm=T),
        meanRef_5x5_90th = mean(Ref_5x5_90th,na.rm=T),
        meanRhoHV_5x5_10th=mean(RhoHV_5x5_10th,na.rm=T),
        meanRhoHV_5x5_50th=mean(RhoHV_5x5_50th,na.rm=T),
        meanRhoHV_5x5_90th=mean(RhoHV_5x5_90th,na.rm=T),
        mean_RefComposite_5x5_10th=mean(RefComposite_5x5_10th,na.rm=T),
        mean_RefComposite_5x5_50th=mean(RefComposite_5x5_50th,na.rm=T),
        mean_RefComposite_5x5_90th=mean(RefComposite_5x5_90th,na.rm=T),
        mean_Zdr = mean(Zdr,na.rm=T),
        mean_Zdr_5x5_10th=mean(Zdr_5x5_10th,na.rm=T),
        mean_Zdr_5x5_50th=mean(Zdr_5x5_50th,na.rm=T),
        mean_Zdr_5x5_90th=mean(Zdr_5x5_90th,na.rm=T),
        mean_Kdp = mean(Kdp,na.rm=T),
        median_Ref=median(dt*Ref,na.rm=T),
        sum_Ref=sum(dt*Ref,na.rm=T)
), Id]
tr<-subset(tr,records-naCounts>0)
tr<-subset(tr,target<expm1(70))

print("training model...")
y<-tr$target
tr<-as.data.frame(tr)
tr<-tr[,3:25]
MAE<- function(preds, dtrain) {
        labels <- getinfo(dtrain, "label")
        elab<-expm1(as.numeric(labels))
        epreds<-expm1(as.numeric(preds))
        err <- mean(abs(epreds-elab))
        return(list(metric = "MAE", value = err))
}
nrow(tr)
h<-sample(nrow(tr),10000)

dval<-xgb.DMatrix(data=data.matrix(tr[h,]),label=y[h],missing=NA) #before missing=0
dtrain<-xgb.DMatrix(data=data.matrix(tr[-h,]),label=y[-h],missing=NA) #before missing=0
watchlist<-list(val=dval,train=dtrain)


param0 <- list("objective"  = "reg:linear" 
               ##, "eval_metric" = "rmse"
               , "booster" = "gbtree"
               , "eta" = 0.010281327 # 0.010281327 0.1
               , "subsample" = 0.8 # 0.8 0.85
               , "min_child_weight" =10    
               , "max_depth" = 8
               #, "colsample_bytree" = 0.6
               , "nthreads" = 4
)
#xgtrain = xgb.DMatrix(as.matrix(tr), label = y, missing = NA)
#gc()

x.mod.t  <- xgb.train(params = param0, data = dtrain , nrounds =2246,
                      verbose             = 1,
                      watchlist = watchlist,
                      maximize            = FALSE,
                      early.stop.round    = 30,
                      feval=MAE)