Cleanup and simplify

realtimeBaseRate = sum(realtimeLabels)/length(realtimeLabels)

preopBaseRate = sum(preopLabels)/length(preopLabels)

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

        if matches != nothing

            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

    end

    featureGroups

end

realtimeFeatureGroups = build_feature_groups(realtimeFeatureNames)

preopFeatureGroups = build_feature_groups(preopFeatureNames)

# get feature groups

realtimeFeatureGroups = feature_groups(realtimeFeatureNames)

preopFeatureGroups = feature_groups(preopFeatureNames)

# pre-allocate space for synthetic samples

preopSampleSpace = samplespace(sparse(preopKmedians), args["num-samples"], 0.5) # assume no more than 50% density

realtimeSampleSpace = samplespace(sparse(realtimeKmedians), args["num-samples"], 0.5) # assume no more than 50% density

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

# 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::SparseMatrixCSC{Float32, Int64})

    handle = Ref{Ptr{Void}}()

    @xgboost_ccall(

        :XGDMatrixCreateFromCSR,

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

        convert(Array{UInt64, 1}, data.colptr - 1),

        convert(Array{UInt32, 1}, data.rowval - 1), data.nzval,

        convert(UInt64, size(data.colptr)[1]),

        convert(UInt64, nnz(data)),

        handle

    )

    return handle[]

end

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

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

            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)

        println("currTimeNormalized = $currTimeNormalized")

    end

    # find which events fall before the given time

    currPos = 1

    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], currTimeNormalized) : 0 for k in featureNames]

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

    predictFunction = data->begin

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

function explained_pred(data, model, featureGroups, X, ntree, nsamples, baseRate, weights, synthSampleSpace, predType)

    predictFunction = x->begin

        XGBoost.predict(model, DMatrix(x, true), ntree_limit=ntree)

    end

    # find the Shapley values of the features

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

        synthSampleSpace=synthSampleSpace,

        featureGroups=featureGroupInds,

        nsamples=args["num-samples"],

        nsamples=nsamples,

        weights=weights

    )

    Dict(

        "baseRate" => baseRate,

        "pred" => logit(predictFunction(sparse(reshape(data, length(data), 1)))[1]),

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

        "predType" => predType,

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

    )

end

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

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

    procId = parts[2]

    println(parts)

    if parts[3] == "preop"

        fdata = build_features(procId, "preop", preopFeatureNames, args["raw-data-dir"])

        res.data = json(explained_pred(

            procId, "preop", preopModel, preopFeatureNames,

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

            preopBaseRate, preopKmediansWeights, preopSampleSpace

            fdata, preopModel, preopFeatureGroups,

            preopKmedians, args["preop-ntree-limit"], args["num-samples"],

            preopBaseRate, preopKmediansWeights, preopSampleSpace, "preop"

        ))

    else

        time = dparse_utc(parts[3])

        time = dparse_utc(URIParser.unescape(parts[3]))

        fdata = build_features(procId, time, realtimeFeatureNames, args["raw-data-dir"])

        res.data = json(explained_pred(

            procId, time, realtimeModel, realtimeFeatureNames,

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

            realtimeBaseRate, realtimeKmediansWeights, realtimeSampleSpace

            fdata, realtimeModel, realtimeFeatureGroups,

            realtimeKmedians, args["realtime-ntree-limit"], args["num-samples"],

            realtimeBaseRate, realtimeKmediansWeights, realtimeSampleSpace, "realtime"

        ))

    end

end

#! /bin/sh

julia index.jl \

    --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-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.0179307 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' \

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

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

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

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

    --preop-ntree-limit 3631 \

    --raw-data-dir '/scratch/slund1/prescience/data/merged/train' \

    --raw-data-dir '/scratch/slund1/prescience/data/merged/validation' \

    --raw-data-dir '/scratch/slund1/prescience/data/merged/test1' \

    --raw-data-dir '/scratch/slund1/prescience/data/merged/test2' \

    --num-samples 1000 \

    --port 5023