Commit before overhaul

        help = "Location of the feature names in the model"

        arg_type = ASCIIString

        required = true

    "--realtime-train-data"

        help = "Location of training feature data for the model"

        arg_type = ASCIIString

        required = true

    # "--realtime-train-data"

    #     help = "Location of training feature data for the model"

    #     arg_type = ASCIIString

    #     required = true

    "--realtime-train-labels"

        help = "Location of training label data for the model"

        arg_type = ASCIIString

        required = true

    "--preop-train-data"

        help = "Location of training feature data for the model"

        arg_type = ASCIIString

        required = true

    # "--preop-train-data"

    #     help = "Location of training feature data for the model"

    #     arg_type = ASCIIString

    #     required = true

    "--preop-train-labels"

        help = "Location of training label data for the model"

        arg_type = ASCIIString

using TimeZones

using StreamingTimeSeries

using GLM

using ShapleyValues

@assert all(map(isdir, args["raw-data-dir"])) "The given raw-data-dir directories were not all found!"

utcTime = TimeZone("UTC")

dparse_utc(str::AbstractString) = ZonedDateTime(DateTime(str[1:end-1], dformat_raw), utcTime)

# load the K-medians we use to approximate the data

realtimeKmedians = convert(Array{Float32, 2}, readdlm(args["realtime-model"][1:end-5]*"20medians.txt"))

realtimeKmediansWeights = convert(Array{Float32, 1}, vec(readdlm(args["realtime-model"][1:end-5]*"20mediansWeights.txt")))

preopKmedians = convert(Array{Float32, 2}, readdlm(args["preop-model"][1:end-5]*"19medians.txt"))

preopKmediansWeights = convert(Array{Float32, 1}, vec(readdlm(args["preop-model"][1:end-5]*"19mediansWeights.txt")))

# pre-load the training data we will sample from

realtimeX = open(deserialize, args["realtime-train-data"])'

#realtimeX = open(deserialize, args["realtime-train-data"])'

realtimeLabels = open(deserialize, args["realtime-train-labels"])

preopX = open(deserialize, args["preop-train-data"])'

#preopX = open(deserialize, args["preop-train-data"])'

preopLabels = open(deserialize, args["preop-train-labels"])

# compute the base occurence rate so we can compare the predicted risks with it

d[1] = 1e8

W = spdiagm(d);

function rand_mask!(maskMemory)

    maskMemory[:,1] = 1.0 # the first sample enforces correct output when all features are included

    maskMemory[:,2] = 0.0 # the second sample enforces correct output when all features are excluded

    maskMemory[end,2] = 1.0

    inds = collect(1:size(maskMemory)[1]-1)

    for i in 3:size(maskMemory)[2]

        numMasked = rand(inds)

# define a data type to hold our synthetic samples

type XGMatrixCSR

    colptr::Array{UInt64,1}

    colptrLen::Int64

    rowval::Array{UInt32,1}

    rowvalLen::Int64

    nzval::Array{Float32,1}

    nzvalLen::Int64

end

XGMatrixCSR(maxLength, maxDensity) = XGMatrixCSR(

    Array(UInt64, maxLength+1), 0,

    Array(UInt32, maxLength*maxDensity), 0,

    Array(Float32, maxLength*maxDensity), 0

)

function reset!(x::XGMatrixCSR)

    x.colptrLen = 1

    x.colptr[1] = 0 # zero-based for easy

    x.rowvalLen = 0

    x.nzvalLen = 0

end

# define a function to map the sythetic samples to an XGBoost DMatrix

macro xgboost_ccall(f, argTypes, args...)

    argTypes = eval(argTypes)

    return quote

        err = ccall(($f, XGBoost._xgboost), Int64, ($(argTypes...),), $(args...))

        if err != 0

            errMsg = bytestring(ccall((:XGBGetLastError, _xgboost), Ptr{UInt8}, ()))

            error("Call to XGBoost C function "*string($f)*" failed: $errMsg")

        end

    end

end

function XGDMatrixCreateFromCSR(data::XGMatrixCSR)

    handle = Ref{Ptr{Void}}()

    @xgboost_ccall(

        :XGDMatrixCreateFromCSR,

        (Ptr{UInt64}, Ptr{UInt32}, Ptr{Float32}, UInt64, UInt64, Ref{Ptr{Void}}),

        data.colptr, data.rowval, data.nzval,

        convert(UInt64, data.colptrLen),

        convert(UInt64, data.nzvalLen),

        handle

    )

    return handle[]

end

# pre-allocate space for synthetic samples

synthLength = 2*args["num-samples"]*size(realtimeKmedians)[2]

synthSamples = XGMatrixCSR(synthLength, round(Int64, size(realtimeKmedians)[1]/2)) # assume no more than 50% density of realtime x

"Distributes the given number of samples proportionally."

function allocate_samples(proportions, nsamples)

    counts = round(Int, nsamples*proportions/sum(proportions))

    total = sum(counts)

    for ind in randperm(length(counts))

        total != nsamples || break

        if total < nsamples

            counts[ind] += 1

            total += 1

        elseif counts[ind] > 0

            counts[ind] -= 1

            total -= 1

        end

    end

    counts

end

# http://www.nowozin.net/sebastian/blog/streaming-mean-and-variance-computation.html

type MeanVarianceAccumulator

    sumw::Float64

    wmean::Float64

    t::Float64

    n::Int

    function MeanVarianceAccumulator()

        new(0.0, 0.0, 0.0, 0)

    end

end

function observe!(mvar::MeanVarianceAccumulator, value, weight)

    @assert weight >= 0.0

    q = value - mvar.wmean

    temp_sumw = mvar.sumw + weight

    r = q*weight / temp_sumw

    mvar.wmean += r

    mvar.t += q*r*mvar.sumw

    mvar.sumw = temp_sumw

    mvar.n += 1

    nothing

end

count(mvar::MeanVarianceAccumulator) = mvar.n

Base.mean(mvar::MeanVarianceAccumulator) = mvar.wmean

var(mvar::MeanVarianceAccumulator) = (mvar.t*mvar.n)/(mvar.sumw*(mvar.n-1))

std(mvar::MeanVarianceAccumulator) = sqrt(var(mvar))

"Push a new sample into a preallocated XGMatrixCSR."

function addsample!(s::XGMatrixCSR, x::Array{Float32,1})

    s.colptrLen += 1

    s.colptr[s.colptrLen] = s.colptr[s.colptrLen-1]

    for i in 1:length(x)

        if x[i] != zero(Float32)

            s.colptr[s.colptrLen] += 1

            s.rowvalLen += 1

            s.rowval[s.rowvalLen] = i-1 # zero based to help XGBoost later

            s.nzvalLen += 1

            s.nzval[s.nzvalLen] = x[i]

        end

    end

end

"The core method that updates the Shapley value estimates designed for XGMatrixCSR caches."

function update_estimates!(deltas, x, f, X::Array{Float32,2}, sampleWeights::Array{Float32,1}, g,

                           featureGroups::Array{Array{Int64,1}}, sampleCounts::Array{Int64,1},

                           synthSampleCache::XGMatrixCSR)

    P,N = size(X)

    M = length(featureGroups)

    @assert length(sampleCounts) == M "sampleCounts should be an array of counts for each feature group!"

    # build the synthentic samples

    inds = collect(1:M)

    s1 = Array(Float32,P)

    s2 = Array(Float32,P)

    reset!(synthSampleCache)

    pos = 1

    for i in 1:M

        for j in 1:sampleCounts[i]

            shuffle!(inds)

        for j in 1:numMasked

            maskMemory[inds[j],i] = 1.0

            # find where in the permutation we are

            ind = findfirst(inds, i)

            for k in 1:N

                # save two synthetic samples with and without the current group replaced

                copy!(s1, x)

                copy!(s2, x)

                for l in ind:M

                    for m in featureGroups[inds[l]]

                        l != ind && (s1[m] = X[m,k])

                        s2[m] = X[m,k]

                    end

        maskMemory[end,i] = 1.0 # the constant offset parameter

                end

                addsample!(synthSampleCache, s1)

                addsample!(synthSampleCache, s2)

                pos += 2

            end

        end

    end

    # run the provided function

    y::Array{Float32,1} = vec(f(synthSampleCache))

function synth_samples!(synthSamplesMemory, x, X, mask, featureGroups)

    for i in 1:size(synthSamplesMemory)[2]

        randPos = rand(1:size(X)[2])

        synthSamplesMemory[:,i] = x

        for j in 1:length(featureGroups)

            mask[j,i] == 0.0 || continue

            for ind in featureGroups[j]

                synthSamplesMemory[ind,i] = X[ind,randPos]

    # sum the totals and keep an estimate of the variance differences

    pos = 1

    withiExp = zero(Float32)

    for i in 1:M

        for j in 1:sampleCounts[i]

            withiExp = zero(Float32)

            withoutiExp = zero(Float32)

            sumw = zero(Float32)

            for k in 1:N

                w = sampleWeights[k]

                withiExp += w*y[pos]

                withoutiExp += w*y[pos+1]

                sumw += w

                pos += 2

            end

            withiExp /= sumw

            withoutiExp /= sumw

            observe!(deltas[i], g(withiExp) - g(withoutiExp), 1)

        end

    end

end

function explain_model3(x, X, basePred, model, linkFunction;

                        nsamples=1000, featureGroups=nothing, maskMemory=nothing,

                        synthSamplesMemory=nothing, weightsMemory=nothing,

                        maskyMemory=nothing)

"Identify which feature groups vary enough in the data to merit estimation."

function varying_groups(x, X, sampleWeights, featureGroups, weightThreshold)

    varyingWeights = zeros(length(featureGroups))

    for (i,inds) in enumerate(featureGroups)

        varyingWeights[i] = sum(sampleWeights[vec(sum(x[inds] .== X[inds,:],1) .!= length(inds))])

    end

    find(varyingWeights .> weightThreshold)

end

    # fill in default values

    featureGroups != nothing || (featureGroups = Array{Int64,1}[Int64[i] for i in 1:length(x)])

    maskMemory != nothing || (maskMemory = zeros(length(featureGroups)+1, nsamples))

    synthSamplesMemory != nothing || (synthSamplesMemory = zeros(length(x), nsamples))

    weightsMemory != nothing || (weightsMemory = zeros(nsamples))

    maskyMemory != nothing || (maskyMemory = zeros(nsamples))

"Designed to determine the Shapley values (importance) of each feature for f(x)."

function shapleyvalues{T}(x, f::Function, X, g::Function=identity, synthSampleCache::T=nothing; featureGroups=nothing,

                          sampleWeights=nothing, nsamples=1000) # maxStdDevFraction=0.02

    P,N = size(X)

    # make sure dimensions all match

    @assert length(featureGroups)+1 == size(maskMemory)[1] "Provided mask memory must match length(featureGroups)+1"

    @assert nsamples == size(maskMemory)[2] "Provided mask memory must match nsamples ($nsamples)"

    @assert length(x) == size(synthSamplesMemory)[1] "Provided synthSamples memory must match length(x)"

    @assert nsamples == size(synthSamplesMemory)[2] "Provided synthSamples memory must match nsamples ($nsamples)"

    @assert nsamples == length(weightsMemory) "Provided weightsMemory memory must match nsamples ($nsamples)"

    @assert nsamples == length(maskyMemory) "Provided maskyMemory memory must match nsamples ($nsamples)"

    # give default values to omitted arguments

    sampleWeights != nothing || (sampleWeights = ones(N))

    featureGroups != nothing || (featureGroups = Array{Int64,1}[Int64[i] for i in 1:length(x)])

    featureGroups = convert(Array{Array{Int64,1},1}, featureGroups)

    # mask features randomly for our samples

    rand_mask!(maskMemory)

    # find the feature groups we will test. If a feature rarely changes from its

    # current value then we know it doesn't have a large impact on the model

    varyingInds = varying_groups(x, X, sampleWeights, featureGroups, 0.01)

    varyingFeatureGroups = featureGroups[varyingInds]

    M = length(varyingFeatureGroups)

    # create samples by randomly replacing entries in the data matrix that have been masked

    synth_samples!(synthSamplesMemory, x, X, maskMemory, featureGroups)

    # loop through the estimation process focusing samples on groups with high variance

    nextSamples = allocate_samples(ones(M), min(10M, nsamples))

    deltas = [MeanVarianceAccumulator() for i in 1:M]

    totalSamples = 0

    while true

    # compute the model predictions for each sample

    println(size(synthSamplesMemory'))

    println(size(vec(model(synthSamplesMemory'))))

    maskyMemory[:] = convert(Array{Float64,1}, vec(model(synthSamplesMemory')))

    maskyMemory[2] = basePred # first entry is model(x), second is base rate

        # update our estimates for a block of samples

        update_estimates!(deltas, x, f, X, sampleWeights, g, varyingFeatureGroups, nextSamples, synthSampleCache)

    # solve for the additive effects

    weightsMemory[:] = 1.0

    weightsMemory[1] = nsamples*10

    weightsMemory[2] = nsamples*10

    m = fit(GeneralizedLinearModel, maskMemory', maskyMemory, Normal(), linkFunction, wts=weightsMemory, minStepFac=1e-20, maxIter=1000)

    coef(m)[1:end-1]

        # keep track of our samples and optimize their allocation to minimize variance (Neyman allocation)

        totalSamples += sum(nextSamples)

        if totalSamples < nsamples

            vs = [var(a) for a in deltas]

            nextSamples = allocate_samples(vs, min(round(Int, nsamples/3), nsamples-totalSamples))

        else break end

    end

function explained_pred(procId, currTime, model, featureNames, featureGroups, X, ntree, baseRate)

    # compute the Shapley values along with estimated variances of the estimates

    φ = zeros(length(featureGroups))

    φ[varyingInds] = [mean(a) for a in deltas]

    φVar = zeros(length(featureGroups))

    φVar[varyingInds] = [var(a)/count(a) for a in deltas]

    # φCounts = zeros(length(featureGroups))

    # φCounts[varyingInds] = [count(a) for a in deltas]

    # find f(x) and E_x[f(x)]

    reset!(synthSampleCache)

    addsample!(synthSampleCache, x)

    fx = f(synthSampleCache)[1]

    reset!(synthSampleCache)

    for i in 1:N addsample!(synthSampleCache, X[:,i]) end

    fnull = sum(f(synthSampleCache).*sampleWeights)

    # We ensure that the total of all features equals f(x)

    trueSum = g(fx) - g(fnull)

    tmp = inv(φ*φ' + I*(trueSum*1e-8))

    β = inv(φ*φ' + I*(trueSum*1e-8))*φ*(sum(φ) - trueSum)

    println(maximum(β))

    φ .-= β.*φ

    @assert all(β .<= 1) "Rescaling failed! (indicates poor Shapley value estimates)"

    # return the Shapley values along with estimated variances of the estimates

    φ,φVar,fx#,φCounts

end

println("DIR ", args["raw-data-dir"], " ", args["num-samples"])

function explained_pred(procId, currTime, model, featureNames, featureGroups, X, ntree, baseRate, sampleWeights, synthSampleCache)

    # load the events

    summaryInfo = nothing

    events = nothing

    for d in args["raw-data-dir"]

        println(d)

        fname = d*"/"*procId*".csv"

        if isfile(fname)

            println("found! $fname")

            f = open(fname)

            summaryInfo = Prescience.SummaryInfoDatum(readline(f))

            seekstart(f)

            return

        end

        currTimeNormalized = max(startTime - Dates.Hour(4), summaryInfo.dos)

        println("currTimeNormalized = $currTimeNormalized")

    end

    # find which events fall before the given time

    features = Dict()

    updatefeatures!(features, events[1:currPos-1])

    data = Float32[haskey(features, k) ? valueat(features[k], currTimeNormalized) : 0 for k in featureNames]

    data = reshape(data, 1, length(data))

    display(find(data))

    println(data)

    predictFunction = data->begin

        println(sum(abs(data),2))

        println("ASDF ", size(sum(data,2)))

        println(minimum(sum(abs(data),2)))

        out = XGBoost.predict(model, sparse(data), ntree_limit=ntree)

    end

    # only check the importance of features that are present in this example (this is way faster)

    groupNames = collect(keys(featureGroups))

    filteredGroupNames = filter(x->!all(data[featureGroups[x]] .== 0), groupNames)

    filteredGroups = [featureGroups[g] for g in filteredGroupNames]

    println(size(filteredGroupNames))

    println(filteredGroups)

    println(args["num-samples"])

    println(baseRate)

    println(X[1,:])

    vals = explain_model3(data, X, baseRate, predictFunction, LogitLink(),

        featureGroups=filteredGroups,

        nsamples=args["num-samples"]

        XGBoost.predict(model, XGDMatrixCreateFromCSR(data), ntree_limit=ntree)

    end

    # find the Shapley values of the features

    featureGroupNames = collect(keys(featureGroups))

    featureGroupInds = Array{Int64,1}[featureGroups[g] for g in keys(featureGroups)]

    vals,vars,fdata = shapleyvalues(data, predictFunction, X, logit, synthSampleCache,

        featureGroups=featureGroupInds,

        nsamples=args["num-samples"],

        sampleWeights=sampleWeights

    )

    println(vals)

    Dict(

        "baseRate" => baseRate,

        "pred" => linkinv(LogitLink(), predictFunction(data)[1]),

        "pred" => logit(fdata),

        "predType" => currTime == "preop" ? "preop" : "realtime",

        "featureEffects" => Dict(zip(filteredGroupNames, vals))

        "featureEffects" => Dict(zip(featureGroupNames, vals))

    )

end

    if parts[3] == "preop"

        res.data = json(explained_pred(

            procId, "preop", preopModel, preopFeatureNames,

            preopFeatureGroups, preopX, args["preop-ntree-limit"],

            preopBaseRate

            preopFeatureGroups, preopKmedians, args["preop-ntree-limit"],

            preopBaseRate, preopKmediansWeights, synthSamples

        ))

    else

        time = dparse_utc(parts[3])

        println("realtimeBaseRate: $realtimeBaseRate")

        res.data = json(explained_pred(

            procId, time, realtimeModel, realtimeFeatureNames,

            realtimeFeatureGroups, realtimeX, args["realtime-ntree-limit"],

            realtimeBaseRate

            realtimeFeatureGroups, realtimeKmedians, args["realtime-ntree-limit"],

            realtimeBaseRate, realtimeKmediansWeights, synthSamples

        ))

    end

end

#! /bin/sh

julia index.jl \

    --realtime-model '../notebooks/data/xgb-2-standard-desat_bool92_5_unfiltered-51min nrounds=3100 subsample=0.5 min_child_weight=10 max_depth=5 base_score=0.026031412 eta=0.01 gamma=1.0.model' \

    --realtime-features '../notebooks/data/features-train-2-standard-desat_bool92_5_unfiltered-51min.txt' \

    --realtime-train-data '../notebooks/data/Xsubset-train-2-standard-desat_bool92_5_unfiltered-51min.jls' \

    --realtime-train-labels '../notebooks/data/ysubset-train-2-standard-desat_bool92_5_unfiltered-51min.jls' \

    --realtime-ntree-limit 3100 \

    --preop-model '../notebooks/data/xgb-2-standard-anydesat_bool92_10_5_5-51min nrounds=2000 subsample=0.5 min_child_weight=1 max_depth=4 base_score=0.078426786 eta=0.01 gamma=1.0.model' \

    --preop-features '../notebooks/data/features-train-2-standard-anydesat_bool92_10_5_5-51min.txt' \

    --preop-train-data '../notebooks/data/X-train-2-standard-anydesat_bool92_10_5_5-51min.jls' \

    --preop-train-labels '../notebooks/data/y-train-2-standard-anydesat_bool92_10_5_5-51min.jls' \

    --preop-ntree-limit 1200 \

    --realtime-model '../notebooks/data/xgb_train_validation-2-standard-desat_bool92_5_nodesat-100min nrounds=2000 subsample=0.5 min_child_weight=10 max_depth=6 base_score=0.017933857 eta=0.02 gamma=1.0.model' \

    --realtime-features '../notebooks/data/features-train_validation-2-standard-desat_bool92_5_nodesat-100min.txt' \

    --realtime-train-labels '../notebooks/data/y-train_validation-2-standard-desat_bool92_5_nodesat-100min.jls' \

    --realtime-ntree-limit 2000 \

    --preop-model '../notebooks/data/xgb_train_validation-2-standard-anydesat_bool92_10_5_5-50min nrounds=2000 subsample=0.5 min_child_weight=1 max_depth=4 base_score=0.07841552 eta=0.01 gamma=1.0.model' \

    --preop-features '../notebooks/data/features-train_validation-2-standard-anydesat_bool92_10_5_5-50min.txt' \

    --preop-train-labels '../notebooks/data/y-train_validation-2-standard-anydesat_bool92_10_5_5-50min.jls' \

    --preop-ntree-limit 1500 \

    --raw-data-dir '/scratch/slund1/prescience/data/merged/merged.rawtrain' \

    --raw-data-dir '/scratch/slund1/prescience/data/merged/merged.rawvalidation' \

    --num-samples 10000 \

    --port 5024

    --raw-data-dir '/scratch/slund1/prescience/data/merged/merged.rawtest1' \

    --num-samples 2000 \

    --port 5023