Use the new Shapley values

using XGBoost

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

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

realtimeKmedians = sparse(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"))

preopKmedians = sparse(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

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 = ones(args["num-samples"])

d[1] = 1e8

W = spdiagm(d);

# 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

# 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

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

        end

    end

end

function XGDMatrixCreateFromCSR(data::XGMatrixCSR)

function XGDMatrixCreateFromCSR(data::SparseMatrixCSC{Float32, Int64})

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

        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

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

            # 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

                end

                addsample!(synthSampleCache, s1)

                addsample!(synthSampleCache, s2)

                pos += 2

            end

        end

    end

    # run the provided function

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

    # 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

"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

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

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

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

    # 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

        # update our estimates for a block of samples

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

        # 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

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

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

    # load the events

    summaryInfo = nothing

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

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

        synthSampleSpace=synthSampleSpace,

        featureGroups=featureGroupInds,

        nsamples=args["num-samples"],

        sampleWeights=sampleWeights

        weights=weights

    )

    Dict(

        "baseRate" => baseRate,

        "pred" => logit(fdata),

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

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

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

    )

        res.data = json(explained_pred(

            procId, "preop", preopModel, preopFeatureNames,

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

            preopBaseRate, preopKmediansWeights, synthSamples

            preopBaseRate, preopKmediansWeights, preopSampleSpace

        ))

    else

        time = dparse_utc(parts[3])

        res.data = json(explained_pred(

            procId, time, realtimeModel, realtimeFeatureNames,

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

            realtimeBaseRate, realtimeKmediansWeights, synthSamples

            realtimeBaseRate, realtimeKmediansWeights, realtimeSampleSpace

        ))

    end

end

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

    --num-samples 2000 \

    --num-samples 1000 \

    --port 5023