Now handles both pre-op and real-time predictions

using ArgParse

s = ArgParseSettings()

@add_arg_table s begin

    "model"

    "--realtime-model"

        help = "Location of XGBoost model to use for prediction"

        arg_type = ASCIIString

        required = true

    "features"

    "--realtime-features"

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

        arg_type = ASCIIString

        required = true

    "train-features"

    "--preop-model"

        help = "Location of XGBoost model to use for prediction"

        arg_type = ASCIIString

        required = true

    "--preop-features"

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

    "train-labels"

    "--preop-train-labels"

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

        arg_type = ASCIIString

        required = true

    "source_dirs"

    "--raw-data-dir"

        help = "Location of the raw data csv files"

        arg_type = ASCIIString

        nargs = '+'

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

        arg_type = Int64

        default = 100

    "--ntree-limit"

    "--mask-rate"

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

        arg_type = Float64

        default = 0.05

    "--realtime-ntree-limit"

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

        arg_type = Int64

        default = 10000

    "--preop-ntree-limit"

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

        arg_type = Int64

        default = 10000

using TimeZones

using StreamingTimeSeries

@assert all(map(isdir, args["source_dirs"])) "The given source_dirs were not all found!"

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

model = Booster(model_file = args["model"])

featureNames = ASCIIString[strip(l) for l in open(readlines, args["features"])]

# load the models

realtimeModel = Booster(model_file = args["realtime-model"])

realtimeFeatureNames = ASCIIString[strip(l) for l in open(readlines, args["realtime-features"])]

preopModel = Booster(model_file = args["preop-model"])

preopFeatureNames = ASCIIString[strip(l) for l in open(readlines, args["preop-features"])]

# pre-define our date parsers

# pre-define our date parser

dformat_raw = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss")

utcTime = TimeZone("UTC")

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

# pre-load the training data we will sample from

trainX = open(deserialize, args["train-features"])

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

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

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

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

baseRate = sum(trainLabels)/length(trainLabels)

realtimeBaseRate = sum(realtimeLabels)/length(realtimeLabels)

preopBaseRate = sum(preopLabels)/length(preopLabels)

# compute groups of features that consider as units

# compute groups of features that we consider as units

function build_feature_groups(featureNames)

    featureGroups = Dict()

    for (i,feature) in enumerate(featureNames)

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

            featureGroups[feature] = [i]

        end

    end

    featureGroups

end

realtimeFeatureGroups = build_feature_groups(realtimeFeatureNames)

preopFeatureGroups = build_feature_groups(preopFeatureNames)

# build a weighting matrix that forces the model to match the first sample

# this causes the sum of all features to be correct

d[1] = 1e8

W = spdiagm(d);

function explained_pred(procId, currTime)

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

    # load the events

    summaryInfo = nothing

    events = nothing

    for d in args["source_dirs"]

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

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

        if isfile(fname)

            events = open(dataevents, fname)

            f = open(fname)

            summaryInfo = Prescience.SummaryInfoDatum(readline(f))

            seekstart(f)

            events = dataevents(f)

            close(f)

            break

        end

    end

    @assert events != nothing "$procId not found!"

    # find a time before this procedure began if this is a preop prediction

    currTimeNormalized = currTime

    if currTime == "preop"

        nowZoned = TimeZones.ZonedDateTime(now(), utcTime)

        startTime = min(

            get(summaryInfo.procStartTime, nowZoned),

            get(summaryInfo.anesthesiaStartTime, nowZoned),

            get(summaryInfo.inRoomTime, nowZoned)

        )

        if startTime == nowZoned || startTime < summaryInfo.dos

            warn("Refusing to predict for a procedure with an invalid start time.")

            return

        end

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

    end

    # find which events fall before the given time

    currPos = 1

    while currPos <= length(events) && events[currPos].time <= currTime

    while currPos <= length(events) && events[currPos].time <= currTimeNormalized

        currPos += 1

    end

    # compute the features

    features = Dict()

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

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

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

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

    # find the model's output for this example

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

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

    # 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)) .> 0.01;

    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[filteredGroupNames[j]] for j in collect(1:length(filteredGroupNames))[vec(!randMasks[i,:])]]...)

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

        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=args["ntree-limit"], output_margin=true) - logit(baseRate)

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

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

    Dict(

        "baseRate" => baseRate,

        "pred" => origPred,

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

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

    )

end

http = HttpHandler() do req::Request, res::Response

    parts = split(req.resource, '/', limit=3)

    procId = parts[2]

    if parts[3] == "preop"

        res.data = json(explained_pred(

            procId, "preop", preopModel, preopFeatureNames,

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

            preopBaseRate

        ))

    else

        time = dparse_utc(parts[3])

    res.data = json(explained_pred(procId, time))

        res.data = json(explained_pred(

            procId, time, realtimeModel, realtimeFeatureNames,

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

            realtimeBaseRate

        ))

    end

end

server = Server(http)

#! /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 \

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

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

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

    --num-samples 10000