Last commit before Shapley switch

    "--num-samples"

        help = "How many samples should be used to estimate the expectation for each held-out feature."

        arg_type = Int64

        default = 100

    "--mask-rate"

        help = "The fraction of held-out features in each sample."

        arg_type = Float64

        default = 0.05

        default = 1000

    "--realtime-ntree-limit"

        help = "Restrict the prediction to use at most this many trees."

        arg_type = Int64

        help = "Restrict the prediction to use at most this many trees."

        arg_type = Int64

        default = 10000

    "--port"

        help = "Which port number to listen on."

        arg_type = Int64

        default = 5023

    "--quiet", "-q"

        help = "Don't print anything"

        action = :store_true

using XGBoost

using TimeZones

using StreamingTimeSeries

using GLM

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

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

# 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

    for (i,feature) in enumerate(featureNames)

        matches = match(r"^(.*)_(ema|emv|decay|ts)[0-9\.]*$", feature)

        if matches != nothing

            featureGroups[matches[1]] = vcat(get(featureGroups, matches[1], Any[]), [i])

            name = matches[1]

            if name in ("TIDALVOLUME", "TV")

                name = "tidalVolume"

            elseif name in ("PEAK", "PEAKPRESSURE", "PIP")

                name = "peakPressure"

            elseif name in ("ETSEV", "ETSEVO")

                name = "ETSEV"

            end

            featureGroups[name] = vcat(get(featureGroups, name, Any[]), [i])

        elseif startswith(feature, "gender_")

            featureGroups["gender"] = vcat(get(featureGroups, "gender", Any[]), [i])

        elseif feature in ("heightInches", "weightPounds", "bmi")

            featureGroups["heightWeight"] = vcat(get(featureGroups, "heightWeight", Any[]), [i])

        elseif feature in ("TIDALVOLUME", "TV")

            println("TV")

            featureGroups["tidalVolume"] = vcat(get(featureGroups, "tidalVolume", Any[]), [i])

        elseif feature in ("PEAK", "PEAKPRESSURE", "PIP")

            featureGroups["peakPressure"] = vcat(get(featureGroups, "peakPressure", Any[]), [i])

        elseif feature in ("ETSEV", "ETSEVO")

            featureGroups["ETSEV"] = vcat(get(featureGroups, "ETSEV", Any[]), [i])

        else

            featureGroups[feature] = [i]

        end

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)

        shuffle!(inds)

        for j in 1:numMasked

            maskMemory[inds[j],i] = 1.0

        end

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

    end

end

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]

            end

        end

    end

end

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

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

                        synthSamplesMemory=nothing, weightsMemory=nothing,

                        maskyMemory=nothing)

    # 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))

    # 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)"

    # mask features randomly for our samples

    rand_mask!(maskMemory)

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

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

    # 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

    # 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]

end

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

    # load the events

    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)))

    # find the model's output for this example

    origPred = XGBoost.predict(model, data, ntree_limit=ntree, output_margin=true)[1]

        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)

    # mask each feature group randomly for our samples

    randMasks = rand(args["num-samples"], length(filteredGroupNames)) .> args["mask-rate"];

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

    # randomly replace entries in the data matrix that have been not masked

    dataMatrix = vcat([data for i in 1:args["num-samples"]]...)

    for i in 1:args["num-samples"]

        inds = vcat([featureGroups[g] for g in filteredGroupNames[vec(!randMasks[i,:])]]...)

        inds2 = vcat([featureGroups[filteredGroupNames[j]] for j in collect(1:length(filteredGroupNames))[vec(!randMasks[i,:])]]...)

        @assert all(inds .== inds2)

        dataMatrix[i,inds] = X[rand(1:size(X)[1]),inds]

    end

    # compute the difference in model predictions from the base rate for each sample and solve for the additive effects

    y = XGBoost.predict(model, dataMatrix, ntree_limit=ntree, output_margin=true) - logit(baseRate)

    vals = inv(randMasks'*W*randMasks + 1e-8*I)*randMasks'*W*y

    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"]

    )

    println(vals)

    Dict(

        "baseRate" => baseRate,

        "pred" => origPred,

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

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

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

    )

        ))

    else

        time = dparse_utc(parts[3])

        println("realtimeBaseRate: $realtimeBaseRate")

        res.data = json(explained_pred(

            procId, time, realtimeModel, realtimeFeatureNames,

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

end

server = Server(http)

run(server, host=IPv4(127,0,0,1), port=5023)

run(server, host=IPv4(127,0,0,1), port=args["port"])

#! /bin/sh

julia index.jl \

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

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

    --realtime-train-data '../notebooks/data/X-train-2-standard-desat_bool92_10_5_5-51min.jls' \

    --realtime-train-labels '../notebooks/data/y-train-2-standard-desat_bool92_10_5_5-51min.jls' \

    --realtime-ntree-limit 1170 \

    --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-ntree-limit 1200 \

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

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

    --num-samples 10000

    --num-samples 10000 \

    --port 5024