First API code version.

using ArgParse

s = ArgParseSettings()

@add_arg_table s begin

    "model"

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

        arg_type = ASCIIString

        required = true

    "features"

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

        arg_type = ASCIIString

        required = true

    "train-features"

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

        arg_type = ASCIIString

        required = true

    "train-labels"

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

        arg_type = ASCIIString

        required = true

    "source_dirs"

        help = "Location of the raw data csv files"

        arg_type = ASCIIString

        nargs = '+'

        required = true

    "--num-samples"

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

        arg_type = Int64

        default = 100

    "--ntree-limit"

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

        arg_type = Int64

        default = 10000

    "--quiet", "-q"

        help = "Don't print anything"

        action = :store_true

end

args = parse_args(s)

using Prescience

using JSON

using HttpServer

using XGBoost

using TimeZones

using StreamingTimeSeries

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

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

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

# pre-define our date parsers

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

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

baseRate = sum(trainLabels)/length(trainLabels)

# compute groups of features that consider as units

featureGroups = Dict()

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

    else

        featureGroups[feature] = [i]

    end

end

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

# this causes the sum of all features to be correct

d = ones(args["num-samples"])

d[1] = 1e8

W = spdiagm(d);

function explained_pred(procId, currTime)

    # load the events

    events = nothing

    for d in args["source_dirs"]

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

        if isfile(fname)

            events = open(dataevents, fname)

            break

        end

    end

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

    # find which events fall before the given time

    currPos = 1

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

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

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

    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)

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

    Dict(

        "baseRate" => baseRate,

        "pred" => origPred,

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

    )

end

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

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

    procId = parts[2]

    time = dparse_utc(parts[3])

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

end

server = Server(http)

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