text:
struct FINITO_lbfgs_adaptive_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg, TH, Ts} <: CIAO_iterable
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                    # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    η::R                    # ls division parameter for γ
    β::R                    # ls division parameter for τ
    sweeping::Int8          # to only use one stepsize γ
    batch::Int              # batch size
    α::R                    # in (0, 1), e.g.: 0.99
    H::TH
    adaptive::Bool          # adaptive stepsizes
    tol_b::R
    S::Ts                   # functions s_i for calculating the Hassian of f_i s
end

mutable struct FINITO_lbfgs_adaptive_state{R<:Real,Tx, TH}
    γ::Tx                   # stepsize parameter, plays the rolde of hat gamma
    hat_γ::R                # average γ
    av::Tx                  # the running average
    ind::Array{Array{Int}}  # running index set
    d::Int                  # number of batches
    H::TH                   # Hessian approx

    # some extra placeholders
    sum_z::Tx
    sum_fx::R
    sum_fv::R
    sum_nrmx::R             # sum |x_i|^2/γ_i
    sum_nabla::Tx           # sum \nabla f(x_i)
    sum_innprod::R          # sum  < \nabla f(x_i), x_i >
    val_fg::R               # value of g
    z::Tx
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx
    zbar::Tx                # bar z
    zbar_prev::Maybe{Tx}    # bar z previous
    vbar::Tx                # bar v
    vbar_prev::Maybe{Tx}    # bar v previous
    dir::Tx                 # direction
    xbar::Tx                # linesearch candidate
    inds::Array{Int}        # needed for shuffled only!
    τ::Float64              # number of epochs
end

function FINITO_lbfgs_adaptive_state(γ::Tx, hat_γ::R, av::Tx, ind, d, H::TH, sum_z, sum_nrmx, sum_nabla, sum_innprod, sum_fx) where {R,Tx,TH}
    return FINITO_lbfgs_adaptive_state{R,Tx,TH}(
        γ,
        hat_γ,
        av,
        ind,
        d,
        H,
        sum_z,
        sum_fx,
        R(0),
        sum_nrmx,
        sum_nabla,
        sum_innprod,
        R(0),
        copy(av),
        copy(av),
        copy(av),
        copy(av),
        nothing,
        copy(av),
        nothing,
        copy(av),
        copy(av),
        collect(1:d),
        1.0
        )
end

function Base.iterate(iter::FINITO_lbfgs_adaptive_iterable{R}) where {R}
    N = iter.N
    r = iter.batch # batch size
    # create index sets
    ind = Vector{Vector{Int}}(undef, 0)
    d = Int(floor(N / r))
    for i = 1:d
        push!(ind, collect(r*(i-1)+1:i*r))
    end
    r * d < N && push!(ind, collect(r*d+1:N))

    #initializing the vectors
    sum_nabla = zero(iter.x0)
    sum_innprod = R(0)
    sum_fx = R(0)
    for i = 1:N  # nabla f(x0)
        ∇f, fi_x = gradient(iter.F[i], iter.x0)
        sum_innprod += real(dot(∇f, iter.x0))
        sum_nabla .+= ∇f
        sum_fx += fi_x
    end

    # updating the stepsize
    γ = iter.γ
    if iter.γ === nothing ##  Lip. of f(x) = 1/N sum nabla f_i(x)
        if iter.L === nothing ||  iter.adaptive
           xeps = iter.x0 .+ one(R)
           av_eps = zero(iter.x0)
           for i in 1:N
                ∇f, ~ = gradient(iter.F[i], xeps)
                av_eps .+= ∇f
            end
            nmg = norm(sum_nabla - av_eps)
            t = 1
            while nmg < eps(R)  # in case xeps has the same gradient
                println("initial upper bound for L is too small")
                xeps = iter.x0 .+ rand(t * [-1, 1], size(iter.x0))
                av_eps = zero(iter.x0)
                for i in 1:N
                    ∇f, ~ = gradient(iter.F[i], xeps)
                    av_eps .+= ∇f
                end
                # grad_f_xeps, f_xeps = gradient(iter.F[i], xeps)
                nmg = norm(sum_nabla - av_eps)
                t *= 2
            end
            ###### decide how to initialize!
            L_int = nmg / (t * sqrt(length(iter.x0)))
            L_int /= iter.N # to account for 1/N
            γ = iter.α / (L_int)
            # γ = iter.S(iter.x0)
            println("gamma specified by L_int: ", γ)
        else
            println("gamma specified by L")
            isa(iter.L, R) ? (γ = iter.α * R(iter.N) / iter.L) :
                (γ = iter.α * R(N) / (maximum(iter.L)) )
        end
    end
    γ = isa(γ, R) ? N * fill(γ, (N,)) : γ
    hat_γ = 1 / sum(1 ./ γ)
    # update the average
    av = copy(sum_nabla) .* (-hat_γ / N)
    av .+= iter.x0

    sum_nrmx = norm(iter.x0)^2 / hat_γ
    sum_z = iter.x0 / hat_γ
    state = FINITO_lbfgs_adaptive_state(γ, hat_γ, av, ind, cld(N, r), iter.H, sum_z, sum_nrmx, sum_nabla, sum_innprod, sum_fx)

    return state, state
end

function Base.iterate(
    iter::FINITO_lbfgs_adaptive_iterable{R},
    state::FINITO_lbfgs_adaptive_state{R},
) where {R}

    if state.zbar_prev === nothing
        state.zbar_prev = zero(state.z)
        state.vbar_prev = zero(state.z)
    end

    γ_init = state.hat_γ

    # full z update
     while true
        if state.hat_γ < iter.tol_b / iter.N
            @warn "parameter `γ` became too small at ls 1, ($(state.hat_γ))"
            return nothing
        end
        prox!(state.zbar, iter.g, state.av, state.hat_γ)

        sum_fz = 0
        for i in 1:iter.N
            sum_fz += iter.F[i](state.zbar)
        end

        ls_temp_l = sum_fz
        ls_temp_l -= real(dot(state.sum_nabla, state.zbar))
        ls_temp_l -= state.sum_fx
        ls_temp_l += state.sum_innprod
        temp = (norm(state.zbar)^2)/(2*state.hat_γ) + (state.sum_nrmx)/2
        ls_temp_r =  (1 + 10^-6) * temp - real(dot(state.sum_z, state.zbar))
        ls_temp_r *= iter.N * iter.α

        tol = 10^(-6)  * (1 + abs(ls_temp_r))
        # break
        R(ls_temp_l) <= ls_temp_r + tol && break
        println("ls 1")
        println("x: ", ls_temp_l)
        println("y: ", ls_temp_r)
        println("hat_γ: ", state.hat_γ)
        println("hat_calc_γ: ", 1 / sum(1 ./ state.γ))
        println("norm: ", state.sum_nrmx)
        println("p1: ", (norm(state.zbar)^2)/(2*state.hat_γ) + (state.sum_nrmx)/2)
        println("p2: ", real(dot(state.sum_z, state.zbar)))

        γ_prev = state.hat_γ
        state.hat_γ *= iter.η
        # update av
        state.av .-=  ((state.hat_γ - γ_prev) / iter.N) .* state.sum_nabla
        reset!(state.H)
    end

    state.sum_fx = R(0)
    state.sum_innprod = R(0)
    state.sum_nabla = zero(state.av)
    for i = 1:iter.N
        state.∇f_temp, fi_z = gradient(iter.F[i], state.zbar) # update the gradient
        state.sum_fx += fi_z
        state.sum_nabla .+= state.∇f_temp
    end
    state.sum_innprod += real(dot(state.sum_nabla,state.zbar))
    state.av .= state.zbar .- (state.hat_γ / iter.N) .* state.sum_nabla
    state.sum_nrmx = norm(state.zbar)^2 * iter.N

    # full v update
    while true
        if state.hat_γ < iter.tol_b / iter.N
            @warn "parameter `γ` became too small at ls 2, ($(state.hat_γ))"
            return nothing
        end
        state.val_fg = prox!(state.vbar, iter.g, state.av, state.hat_γ) # bar v update
        state.sum_fv = 0
        for i in 1:iter.N
            state.sum_fv += iter.F[i](state.vbar)
        end

        f_model_cnst = state.sum_nrmx / (2 * state.hat_γ) + state.sum_fx - state.sum_innprod

        f_model_z = real(dot(state.vbar, state.av))
        f_model_z -=  norm(state.vbar)^2 / 2
        f_model_z *= iter.N / state.hat_γ

        f_model_z += state.sum_fv
        tol = 10^(-6)  * (1 + abs(f_model_cnst))
        # tol = eps(R)
        # break
        R(f_model_z) <= f_model_cnst + tol && break
        println("ls 2")

        γ_prev = state.hat_γ
        state.hat_γ *= iter.η
        # update av
        state.av .-=  ((state.hat_γ - γ_prev) / iter.N) .* state.sum_nabla
        reset!(state.H)
    end
    state.vbar .-= state.zbar # \bar v- \bar z

    # prepare for linesearch: compute varphi(bar z)
    envVal = state.sum_fx / iter.N
    envVal += state.val_fg
    envVal += real(dot(state.sum_nabla, state.vbar)) / iter.N
    envVal += norm(state.vbar)^2 / (2 *  state.hat_γ)

    # update lbfgs
    update!(state.H, state.zbar - state.zbar_prev, -state.vbar +  state.vbar_prev)
    mul!(state.dir, state.H, state.vbar)
    # store vectors for next update
    copyto!(state.zbar_prev, state.zbar)
    copyto!(state.vbar_prev, state.vbar)

    state.τ = 1
    γ_prev = state.hat_γ
    # while true  ### z_line should be removed using some other placeholder
    for i=1:5
        state.hat_γ = γ_prev
        state.xbar .=  state.zbar .+ (1- state.τ) .* state.vbar + state.τ * state.dir

        # compute varphi(xbar)
        state.sum_nabla = zero(state.temp)
        state.sum_fx = 0
        for i = 1:iter.N
            state.∇f_temp, fi_z = gradient(iter.F[i], state.xbar) # update the gradient
            state.sum_fx += fi_z / iter.N
            state.sum_nabla .+= state.∇f_temp
        end

        while true
            if state.hat_γ < iter.tol_b / iter.N
                @warn "parameter `γ` became too small at ls 3, ($(state.hat_γ))"
                return nothing
            end
            state.av .= state.xbar
            state.av .-= (state.hat_γ / iter.N) .* state.sum_nabla
            state.val_fg = prox!(state.z, iter.g, state.av, state.hat_γ)

            temp = 0
            for i = 1:iter.N
                ~, fi_z = gradient(iter.F[i], state.z) # update the gradient
                temp += fi_z / iter.N
            end

            envVal_trial = 0
            envVal_trial += state.sum_fx
            state.z .-= state.xbar # r_γ(z^trial)
            envVal_trial += real(dot(state.sum_nabla, state.z)) / iter.N
            envVal_trial += norm(state.z)^2 / (2 *  state.hat_γ)

            tol = 10^(-6)  * (1 + abs(envVal_trial))
            # tol = eps(R)
            # break
            if temp <= envVal_trial + tol
                # println("ls condition met")
                break
            end
            state.hat_γ *= iter.η
            println("ls 3")
        end
        envVal_trial += state.val_fg

        tol = 10^(-6) * abs(envVal)
        envVal_trial <= envVal + tol && break
        state.τ *= iter.β      ##### bug prone: change in reporting if you change this! ######
        println("ls on τ")
        println("envVal_trial: ", envVal_trial)
        println("envVal      : ", envVal)
    end
    state.zbar .= state.xbar # of line search
    ###### at the momoent I am computing the next prox twice (it is already computed in the ls)
    iter.sweeping == 3 && (state.inds = randperm(state.d)) # shuffled
    state.temp .= state.av

    state.sum_fx = R(0)
    state.sum_z = zero(state.temp)
    state.sum_innprod =  R(0)
    state.sum_nrmx =  R(0)
    while_flag = false
    flag = true

    state.γ *= (state.hat_γ / γ_init)
    state.zbar .= state.xbar

    for j in state.inds
        for i in state.ind[j]
            # backtrack γ (warn if γ gets too small)

            # println("γ_i: ", iter.S(state.zbar, i)/iter.N)
            #--------------------------- restarting section ---------------------------
            # hat_γ_prev = state.hat_γ
            # γ_prev = state.γ[i]
            # # gradient!(state.∇f_temp, iter.F[i], state.zbar)
            # # state.γ[i] = iter.S(state.zbar, state.∇f_temp, i)
            # state.γ[i] = iter.S(state.zbar, i)
            # state.hat_γ = 1/(1/state.hat_γ - 1/γ_prev + 1/state.γ[i])
            # # state.hat_γ = 1 / sum(1 ./ state.γ)
            # state.av .+= (hat_γ_prev - state.hat_γ) / iter.N * state.sum_nabla + (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev) * state.xbar
            #--------------------------- restarting section ---------------------------

            # print("current gamma: $(state.hat_γ) - current index: $(state.γ[i])")
            while true
                if state.hat_γ < iter.tol_b / iter.N
                    @warn "parameter `γ` became too small in ls 4 ($(state.hat_γ)) - for index $(state.γ[i])"
                    return nothing
                end
                prox!(state.zbar, iter.g, state.av, state.hat_γ)

                fi_x = gradient!(state.∇f_temp, iter.F[i], state.xbar) # update the gradient
                state.val_fg = iter.F[i](state.zbar)  ##### change to function value! #####

                fi_model = fi_x + real(dot(state.∇f_temp, state.zbar .- state.xbar)) +
                    (0.5 * iter.α * iter.N/ state.γ[i]) * (norm(state.zbar .- state.xbar)^2)

                tol = 10^(-6)  * (1 + abs(fi_model))
                # break
                R(state.val_fg) <= fi_model + tol && break

                println("ls 4 γ: $(state.hat_γ) for sample $(i)")
                println("x: ", state.val_fg)
                println("y: ", fi_model)

                hat_γ_prev = state.hat_γ
                γ_prev = state.γ[i]
                state.γ[i] *= iter.η # detach point!!
                state.hat_γ = 1 / sum(1 ./ state.γ)
                state.av .+= (hat_γ_prev - state.hat_γ) / iter.N * state.sum_nabla + (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev) * state.xbar

                reset!(state.H)
            end

            # iterates
            state.av .+= (state.hat_γ / iter.N) .* state.∇f_temp
            state.sum_nabla .-= state.∇f_temp # update sum nabla f_i for next iter
            gradient!(state.∇f_temp, iter.F[i], state.zbar) # update the gradient
            state.sum_nabla .+= state.∇f_temp # update sum nabla f_i for next iter
            state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
            state.av .+= state.hat_γ * (state.zbar .- state.xbar) ./ state.γ[i]

            # updates for the next linesearch
            state.sum_z += state.zbar / state.γ[i]
            state.sum_fx += state.val_fg
            state.sum_nrmx   += norm(state.zbar)^2 / state.γ[i]
            state.sum_innprod += real(dot(state.∇f_temp, state.zbar))
        end
    end
    # println("current stepsize is $(state.hat_γ)")
    return state, state
end
solution(state::FINITO_lbfgs_adaptive_state) = state.zbar
epoch_count(state::FINITO_lbfgs_adaptive_state) = state.τ   # number of epochs is 2+ 1/tau , where 1/tau is from ls

#### TODO:
    # careful with the minibatch and adaptive stepsizes!
    # stepsizes if gamma is given...
    # add LS parameter as iter field


text:
function optimal_simple(d::Divergence, m::DRO)
    p = zero(m.q)
    p[m.Imax] .= m.q[m.Imax]
    return p./sum(m.q[m.Imax])
end

# ----------------------------------------------------------------------------------------------------------
# Kullback-Leibler divergence
# ----------------------------------------------------------------------------------------------------------
"""
    KullbackLeibler

An empty structure representing Kullback-Leibler divergence.
"""
struct KullbackLeibler <: Divergence end

"""
    generate(d::KullbackLeibler)

Returns the generating function `ϕ(t) = t⋅log(t)` of the Kullback-Leibler divergence.
The generating function for `t == 0` is defined as `0`.
"""
generate(d::KullbackLeibler) = ϕ(t) = iszero(t) ? zero(t) : t*log(t)

"""
    name(d::KullbackLeibler)

Returns the full name of the Kullback-Leibler divergence.
"""
name(d::KullbackLeibler) = "Kullback-Leibler divergence"

"""
    bounds(d::KullbackLeibler, m::DRO)

Returns lower and upper bound for finding the root of the function `h` using the bisection method.
"""
bounds(d::KullbackLeibler, m::DRO) = (0, (m.cmax - m.cmin)/m.ε)

function h(d::KullbackLeibler, m::DRO, λ::Real)
    p_hat = m.q.*exp.(m.c./λ)
    val   = p_hat' * (m.c./λ .- log(sum(p_hat)) .- m.ε)
    return isnan(val) || val == - typemax(val) ? typemax(val) : val
end

"""
    function optimal(s::Our, d::KullbackLeibler, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with Kullback-Leibler divergence.
"""
function optimal(s::Our, d::KullbackLeibler, m::DRO; kwargs...)
    p = optimal_simple(d, m)
    ϕ = generate(d)

    if sum(m.q .* ϕ.(p ./ m.q)) <= m.ε
        return p
    else
        λ = bisection(d, m; kwargs...)
        p = m.q.*exp.(m.c./λ)
        return p./sum(p)
    end
end

# ----------------------------------------------------------------------------------------------------------
# Burg entropy
# ----------------------------------------------------------------------------------------------------------
"""
    Burg

An empty structure representing Burg entropy.
"""
struct Burg <: Divergence end

"""
    generate(d::Burg, t)

Returns the generating function `ϕ(t) = - log(t)` of the Burg entropy.
"""
generate(d::Burg) = ϕ(t) = - log(t)

"""
    name(d::Burg)

Returns the full name of the Burg entropy.
"""
name(d::Burg) = "Burg entropy"

"""
    bounds(d::Burg, m::DRO)

Returns lower and upper bound for finding the root of the function `h` using the bisection method.
"""
bounds(d::Burg, m::DRO) = (m.cmax, m.cmax + (m.cmax - m.cmin)/m.ε)

function h(d::Burg, m::DRO, λ::Real)
    val = sum(m.q.*log.(λ .- m.c)) + log(sum(m.q./(λ .- m.c))) - m.ε
    return isnan(val) ? typemax(val) : val
end

"""
    function optimal(s::Our, d::Burg, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with Burg entropy.
"""
function optimal(s::Our, d::Burg, m::DRO; kwargs...)
    p = optimal_simple(d, m)
    ϕ = generate(d)

    if sum(m.q .* ϕ.(p ./ m.q)) <= m.ε
        return p
    else
        λ = bisection(d, m; kwargs...)
        p = m.q./(λ .- m.c)
        return p./sum(p)
    end
end

# ----------------------------------------------------------------------------------------------------------
# Hellinger distance
# ----------------------------------------------------------------------------------------------------------
"""
    Hellinger

An empty structure representing Hellinger distance.
"""
struct Hellinger <: Divergence end

"""
    generate(d::Hellinger, t)

Returns the generating function `ϕ(t) = (√t - 1)²`  of the Hellinger distance.
"""
generate(d::Hellinger) = ϕ(t) = (sqrt(t) - 1)^2

"""
    name(d::Hellinger)

Returns the full name of the Hellinger distance.
"""
name(d::Hellinger) = "Hellinger distance"

"""
    bounds(d::Hellinger, m::DRO)

Returns lower and upper bound for finding the root of the function `h` using the bisection method.
"""
bounds(d::Hellinger, m::DRO) = (m.cmax, m.cmax + (2 - m.ε)*(m.cmax - m.cmin)/m.ε)

function h(d::Hellinger, m::DRO, λ::Real)
    val =  2*sum(m.q./(λ .- m.c)) - (2 - m.ε)*sqrt(sum(m.q./((λ .- m.c).^2)))
    return isnan(val) ? - typemax(val) : val
end

"""
    function optimal(s::Our, d::Hellinger, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with Hellinger distance.
"""
function optimal(s::Our, d::Hellinger, m::DRO; kwargs...)
    p = optimal_simple(d, m)
    ϕ = generate(d)

    if sum(m.q .* ϕ.(p ./ m.q)) <= m.ε
        return p
    else
        λ = bisection(d, m; kwargs...)
        p = m.q./(λ .- m.c).^2
        return p./sum(p)
    end
end

# ----------------------------------------------------------------------------------------------------------
# χ²-distance
# ----------------------------------------------------------------------------------------------------------
"""
    ChiSquare

An empty structure representing χ²-distance.
"""
struct ChiSquare <: Divergence end

"""
    generate(d::ChiSquare, t)

Returns the generating function `ϕ(t) = (t - 1)²/t`  of the χ²-distance.
"""
generate(d::ChiSquare) = ϕ(t) = ((t - 1)^2)/t

"""
    name(d::ChiSquare)

Returns the full name of the χ²-distance.
"""
name(d::ChiSquare) = "χ²-distance"

"""
    initial(d::ChiSquare, m::DRO)

Returns the initial point for finding the root of the function `h` using the newton method.
"""
function initial(d::ChiSquare, m::DRO)
    f(μ)  = h(d, m, μ)
    λ_min = m.cmax
    q     = 0.01
    f_val = f(λ_min + q)
    add_eval()

    while f_val <= 0
        f == 0 && return λ_min + q

        q    /= 10
        f_val = f(λ_min + q)
        add_eval()
    end
    return λ_min + q
end

function h(d::ChiSquare, m::DRO, λ::Real)
    val = sum(m.q .* sqrt.(λ .- m.c))*sum(m.q./sqrt.(λ .- m.c)) - 1 - m.ε
    return isnan(val) ? - typemax(val) : val
end

function ∇h(d::ChiSquare, m::DRO, λ::Real)
    λc = λ .- m.c
    return (sum(m.q ./ sqrt.(λc))^2 - sum(m.q .* sqrt.(λc))*sum(m.q ./ (λc.^(3/2))))/2
end

"""
    function optimal(s::Our, d::ChiSquare, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with χ²-distance.
"""
function optimal(s::Our, d::ChiSquare, m::DRO; kwargs...)
    p = optimal_simple(d, m)
    ϕ = generate(d)

    if sum(m.q .* ϕ.(p ./ m.q)) <= m.ε
        return p
    else
        λ = newton(d, m; kwargs...)
        p = m.q./sqrt.(λ .- m.c)
        return p./sum(p)
    end
end

# ----------------------------------------------------------------------------------------------------------
# Modified χ²-distance
# ----------------------------------------------------------------------------------------------------------
"""
    ModifiedChiSquare

An empty structure representing modified χ²-distance.
"""
struct ModifiedChiSquare <: Divergence end

"""
    generate(d::ModifiedChiSquare, t)

Returns the generating function `ϕ(t) = (t - 1)²`  of the modified χ²-distance.
"""
generate(d::ModifiedChiSquare) = ϕ(t) = (t - 1)^2

"""
    name(d::ModifiedChiSquare)

Returns the full name of the modified χ²-distance.
"""
name(d::ModifiedChiSquare) = "Modified χ²-distance"

"""
    initial(d::ModifiedChiSquare, m::DRO)

Returns the initial point for finding the root of the function `h` using the newton method.
"""
function initial(d::ModifiedChiSquare, m::DRO)
    f(λ)  = h(d, m, λ)
    λ     = - m.cmax + 10
    f_val = f(λ)
    add_eval()

    while f_val >= 0
        f == 0 && return λ

        λ    *= 10
        f_val = f(λ)
        add_eval()
    end
    return λ
end

function h(d::ModifiedChiSquare, m::DRO, λ::Real)
    val = sum(m.q .* max.(λ .+ m.c, 0).^2) - (1 + m.ε) * (sum(m.q .* max.(λ .+ m.c, 0)))^2
    return isnan(val) ? - typemax(val) : val
end

function ∇h(d::ModifiedChiSquare, m::DRO, λ::Real)
    I = findall(λ .+ m.c .>= 0)
    return 2 * sum(m.q[I] .* (λ .+ m.c[I])) * (1 - (1 + m.ε) * sum(m.q[I]))
end

"""
    function optimal(s::Our, d::ModifiedChiSquare, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with modified χ²-distance.
"""
function optimal(s::Our, d::ModifiedChiSquare, m::DRO; kwargs...)
    p = optimal_simple(d, m)
    ϕ = generate(d)

    if sum(m.q .* ϕ.(p ./ m.q)) <= m.ε
        return p
    else
        λ = newton(d, m; kwargs...)
        p = m.q.*max.(λ .+ m.c, 0)
        return p./sum(p)
    end
end

ext: jl

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
# ---
# title: 235. Lowest Common Ancestor of a Binary Search Tree
# id: problem235
# author: zhwang
# date: 2022-02-23
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/lowest-common-ancestor-of-a-binary-search-tree/description/>
# hidden: true
# ---
#
# Given a binary search tree (BST), find the lowest common ancestor (LCA) of two
# given nodes in the BST.
#
# According to the [definition of LCA on
# Wikipedia](https://en.wikipedia.org/wiki/Lowest_common_ancestor): "The lowest
# common ancestor is defined between two nodes p and q as the lowest node in T
# that has both p and q as descendants (where we allow **a node to be a
# descendant of itself** )."
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/14/binarysearchtree_improved.png)
#
#
#
#     Input: root = [6,2,8,0,4,7,9,null,null,3,5], p = 2, q = 8
#     Output: 6
#     Explanation: The LCA of nodes 2 and 8 is 6.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/14/binarysearchtree_improved.png)
#
#
#
#     Input: root = [6,2,8,0,4,7,9,null,null,3,5], p = 2, q = 4
#     Output: 2
#     Explanation: The LCA of nodes 2 and 4 is 2, since a node can be a descendant of itself according to the LCA definition.
#
#
# **Example 3:**
#
#
#
#     Input: root = [2,1], p = 2, q = 1
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[2, 105]`.
#   * `-109 <= Node.val <= 109`
#   * All `Node.val` are **unique**.
#   * `p != q`
#   * `p` and `q` will exist in the BST.
#
#
## @lc code=start
using LeetCode

function lowest_common_ancestor_235(root::TreeNode, p::TreeNode, q::TreeNode)::TreeNode
    lv, gv = p.val < q.val ? (p.val, q.val) : (q.val, p.val)
    while true
        lv <= root.val <= gv && return root
        root = root.val < lv ? root.right : root.left
    end
end
## @lc code=end

## @lc test=start
@testset "235.lowest-common-ancestor-of-a-binary-search-tree.jl" begin
    root = TreeNode{Int}([6,2,8,0,4,7,9,nothing,nothing,3,5])
    @test lowest_common_ancestor_235(root, root.left, root.right) == root
    @test lowest_common_ancestor_235(root, root.left, root.left.right) == root.left
    root = TreeNode{Int}([2, 1])
    @test lowest_common_ancestor_235(root, root, root.left) == root
end
## @lc test=end

## @lc info=start
# link: [solution to 235](https://leetcode-cn.com/problems/lowest-common-ancestor-of-a-binary-search-tree/solution/235-er-cha-sou-suo-shu-de-zui-jin-gong-g-yne3/)
## @lc info=end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

using ProgressMeter, Plots
using GEMPIC

# # Strong Landau Damping
#
# Electrostatic example of strong Landau damping
# $$
# f(x,v) =\frac{1}{2\pi\sigma^2}  \exp
# \Big( - \frac{v_1^2 + v_2^2}{2\sigma^2} \Big)
# ( 1+\alpha \cos(k x)􏰁),
# $$
#
# The physical parameters

# +
σ, μ = 1.0, 0.0
kx, α = 0.5, 0.5
xmin, xmax = 0, 2π/kx
∆t = 0.05
nx = 32
n_particles = 100000
mesh = Mesh( xmin, xmax, nx)
spline_degree = 3

df = CosSumGaussian{1,2}([[kx]], [α], [[σ,σ]], [[μ,μ]] )

mass, charge = 1.0, 1.0
particle_group = ParticleGroup{1,2}(n_particles)
sampler = ParticleSampler{1,2}( :sobol, true, n_particles)

sample!(particle_group, sampler, df, mesh)

xp = view(particle_group.array, 1, :)
vp = view(particle_group.array, 2:3, :)
wp = view(particle_group.array, 4, :);
# -

p = plot(layout=(3,1))
histogram!(p[1,1], xp, weights=wp, normalize=true, bins = 100, lab = "")
plot!(p[1,1], x-> (1+α*cos(kx*x))/(2π/kx), 0., 2π/kx, lab="")
histogram!(p[2,1], vp[1,:], weights=wp, normalize=true, bins = 100, lab = "")
plot!(p[2,1], v-> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="")
histogram!(p[3,1], vp[2,:], weights=wp, normalize=true, bins = 100, lab = "")
plot!(p[3,1], v-> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="")
savefig("histograms.png")

# Initialize the arrays for the spline coefficients of the fields

kernel_smoother1 = ParticleMeshCoupling( mesh, n_particles, spline_degree-1, :galerkin)
kernel_smoother0 = ParticleMeshCoupling( mesh, n_particles, spline_degree, :galerkin)

# Initialize the field solver
rho = zeros(Float64, nx)
# efield by Poisson
efield_poisson = zeros(Float64, nx)
maxwell_solver = Maxwell1DFEM(mesh, spline_degree)
solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )
xg = LinRange(xmin, xmax, nx)
sval = eval_uniform_periodic_spline_curve(spline_degree-1, rho)
p = plot(layout=(1,2))
plot!(p[1,1], xg, sval )
sval = eval_uniform_periodic_spline_curve(spline_degree-1, efield_poisson)
plot!(p[1,2], xg, sval )

# +
function run( steps)

    σ, μ = 1.0, 0.0
    kx, α = 0.5, 0.5
    xmin, xmax = 0, 2π/kx
    dt = 0.02
    nx = 32
    n_particles = 200000
    mesh = Mesh( xmin, xmax, nx)
    spline_degree = 3

    df = CosSumGaussian{1,2}([[kx]], [α], [[σ,σ]], [[μ,μ]] )

    mass, charge = 1.0, 1.0
    particle_group = ParticleGroup{1,2}(n_particles)
    sampler = ParticleSampler{1,2}( :sobol, true, n_particles)

    sample!(particle_group, sampler, df, mesh)

    kernel_smoother1 = ParticleMeshCoupling( mesh, n_particles, spline_degree-1, :galerkin)
    kernel_smoother0 = ParticleMeshCoupling( mesh, n_particles, spline_degree, :galerkin)

    rho = zeros(Float64, nx)
    efield_poisson = zeros(Float64, nx)
    # Initialize the field solver
    maxwell_solver = Maxwell1DFEM(mesh, spline_degree)
    # efield by Poisson
    solve_poisson!( efield_poisson, particle_group, kernel_smoother0, maxwell_solver, rho )

    efield_dofs = [efield_poisson, zeros(Float64, nx)]
    bfield_dofs = zeros(Float64, nx)

    propagator = HamiltonianSplitting( maxwell_solver,
                                       kernel_smoother0,
                                       kernel_smoother1,
                                       particle_group,
                                       efield_dofs,
                                       bfield_dofs)

    efield_dofs_n = propagator.e_dofs

    thdiag = TimeHistoryDiagnostics( particle_group, maxwell_solver,
                            kernel_smoother0, kernel_smoother1 )

    @showprogress 1 for j = 1:steps # loop over time

        # Strang splitting
        strang_splitting!(propagator, dt, 1)

        # Diagnostics
        solve_poisson!( efield_poisson, particle_group,
                        kernel_smoother0, maxwell_solver, rho)

        write_step!(thdiag, j * dt, spline_degree,
                        efield_dofs,  bfield_dofs,
                        efield_dofs_n, efield_poisson)

    end

    return thdiag.data

end
# -

@time results = run(2500) # change number of steps

plot(results[!,:Time], log.(results[!,:PotentialEnergyE1]))

savefig("potential_energy.png")


text:
"""
    function addUpdateLocal_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions related with Local Updates.
"""
function addUpdateLocal_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if "UpdateLocal" in keys(p.agent.declaredUpdates)

        code = p.agent.declaredUpdates["UpdateLocal"]

        #Construct functions
        f = simpleFirstLoopWrapInFunction_(platform,:locStep_!,code)
        f = vectorize_(p.agent,f,p)

        push!(p.declareF.args,
            f)

        push!(p.execInloop.args,
                :(@platformAdapt locStep_!(ARGS_))
            )

    end

    return nothing
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
const MONOKAY = Dict(
            "COLOR_ASPHALT"              => colorant"0x708090",
            "COLOR_LANE_MARKINGS_WHITE"  => colorant"0xDFDFDF",
            "COLOR_LANE_MARKINGS_YELLOW" => colorant"0xFFEF00",
            "COLOR_CAR_EGO"              => colorant"0x00FFFF",
            "COLOR_CAR_OTHER"            => colorant"0xFF007F",
            "COLOR_ROAD_BOUNDARIES"      => colorant"white",
            "foreground" => colorant"0xCFBFAD",
            "background" => colorant"black",
            "color1"     => colorant"0x52E3F6", # light blue
            "color2"     => colorant"0xA7EC21", # light green
            "color3"     => colorant"0xFF007F", # red
            "color4"     => colorant"0xF9971F", # orange
            "color5"     => colorant"0x79ABFF", # cobalt
    )

const OFFICETHEME = Dict(
            "COLOR_ASPHALT"              => colorant"0x708090",
            "COLOR_LANE_MARKINGS_WHITE"  => colorant"0xDFDFDF",
            "COLOR_LANE_MARKINGS_YELLOW" => colorant"0xFFEF00",
            "COLOR_CAR_EGO"              => colorant"#0090bc",
            "COLOR_CAR_OTHER"            => colorant"#ce0300" ,
            "COLOR_ROAD_BOUNDARIES"      => colorant"white",
            "foreground" => colorant"0xCFBFAD",
            "background" => colorant"transparent",
            "color1"     => colorant"#0090bc", # light blue
            "color2"     => colorant"#00a321", # light green
            "color3"     => colorant"#ce0300", # red
    )

const LIGHTTHEME = Dict(
            "COLOR_ASPHALT"              => colorant"transparent",
            "COLOR_LANE_MARKINGS_WHITE"  => colorant"black",
            "COLOR_LANE_MARKINGS_YELLOW" => colorant"0xFFEF00",
            "COLOR_CAR_EGO"              => colorant"#0090bc",
            "COLOR_CAR_OTHER"            => colorant"#ce0300" ,
            "COLOR_ROAD_BOUNDARIES"      => colorant"white",
            "foreground" => colorant"0xCFBFAD",
            "background" => colorant"transparent",
            "color1"     => colorant"#0090bc", # light blue
            "color2"     => colorant"#00a321", # light green
            "color3"     => colorant"#ce0300", # red
    )

global colortheme = MONOKAY

"""
    set_color_theme(colortheme)
Change the color theme of the package
"""
function set_color_theme(ct)
    global colortheme
    colortheme = ct
end

COLOR_ASPHALT       = colortheme["COLOR_ASPHALT"]
COLOR_LANE_MARKINGS_WHITE = colortheme["COLOR_LANE_MARKINGS_WHITE"]
COLOR_LANE_MARKINGS_YELLOW = colortheme["COLOR_LANE_MARKINGS_YELLOW"]
COLOR_CAR_EGO       = colortheme["COLOR_CAR_EGO"]
COLOR_CAR_OTHER     = colortheme["COLOR_CAR_OTHER"]

function Vec.lerp(a::Colorant, b::Colorant, t::Real)

    ra = red(a)
    rb = red(b)
    ga = green(a)
    gb = green(b)
    ba = blue(a)
    bb = blue(b)
    aa = alpha(a)
    ab = alpha(b)

    r = ra + (rb - ra)*t
    g = ga + (gb - ga)*t
    b = ba + (bb - ba)*t
    a = aa + (ab - aa)*t

    RGBA(r,g,b,a)
end

function get_pastel_car_colors(scene::EntityFrame{S,D,I}; saturation::Float64=0.85, value::Float64=0.85) where {S,D,I}
    retval = Dict{I,Colorant}()
    n = length(scene)
    for (i,veh) in enumerate(scene)
        retval[veh.id] = convert(RGB, HSV(180*(i-1)/max(n-1,1), saturation, value))
    end
    return retval
end

"""
    id_to_color(id)

Random color based on hash code of the ID
see https://stackoverflow.com/questions/11120840/hash-string-into-rgb-color
"""
function id_to_color(id)
    idhash = hash(id)
    color = RGB(
        .3 + .7*((idhash & 0xFF0000) >> 16)/255,
        .3 + .7*((idhash & 0x00FF00) >> 8)/255,
        .3 + .7*((idhash & 0x0000FF))/255,
    )
    return color
end


text:
using DynamicalSystems

# Short time series of two coupled logistic maps (x drives y)
sys = logistic2_unidir(c_xy = 0.5)
npts = 40
x, y = columns(trajectory(sys, npts, Ttr = 200))

# average over a few different binning schemes, use a few different
# prediction lags and use the transfer operator grid estimator on the point
# cloud generated from the invariant measure over the triangulation.
binnings = [RectangularBinning(i) for i = 2:3]
estimator = TransferOperatorGrid()
ηs = 1:3

# Perform causality test, both from x to y and y to x
test = ApproximateSimplexIntersectionTest(binning = binnings, estimator = estimator, ηs = ηs)

@test causality(x, y, test) isa Vector{T} where {T}


text:
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n
using .MyPrintf

const DIGIT = 100
const MAXN = 500

function viete(n)
    numerator = sqrt(BigFloat(2))
    result = BigFloat(1)

    for i = 1:n
        result *= (numerator / BigFloat(2))
        numerator = sqrt(BigFloat(2) + numerator)
    end

    return BigFloat(2) / result
end

function main()
    n = GetDigit_n.getn(viete, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(viete, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(viete, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, viete(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
function check_err(err)
    err[1] != C_NULL && throw(RocksDBException(error(unsafe_string(err[1]))))
    nothing
end

function open_db(file_path, create_if_missing)
    options = @threadcall( (:rocksdb_options_create, librocksdb), Ptr{Nothing}, ())
    if create_if_missing
        @threadcall( (:rocksdb_options_set_create_if_missing, librocksdb), Nothing,
              (Ptr{Nothing}, UInt8), options, 1)
    end
    err = Ptr{UInt8}[0]
    db = @threadcall( (:rocksdb_open, librocksdb), Ptr{Nothing},
               (Ptr{Nothing}, Ptr{UInt8}, Ptr{Ptr{UInt8}}) , options, file_path, err)

    @threadcall( (:rocksdb_options_destroy, librocksdb), Nothing, (Ptr{Nothing},), options)

    db == C_NULL && check_err(err)

    return db
end

function close_db(db)
    @threadcall( (:rocksdb_close, librocksdb), Nothing, (Ptr{Nothing},), db)
end

function db_put(db, key, value; raw=false)
    options = @threadcall( (:rocksdb_writeoptions_create, librocksdb), Ptr{Nothing}, ())
    k = byte_array(key, endian_conv=true)
    v = raw ? value : byte_array(value)
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_put, librocksdb), Nothing,
          (Ptr{Nothing}, Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt, Ptr{Ptr{UInt8}} ),
          db, options,k, length(k), v, length(v), err)
    check_err(err)
end

function db_put_sync(db, key, value; raw=false)
    options = @threadcall( (:rocksdb_writeoptions_create, librocksdb), Ptr{Nothing}, ())
    @threadcall( (:rocksdb_writeoptions_set_sync, librocksdb), Nothing,
                 (Ptr{Nothing}, UInt8), options, 1)
    k = byte_array(key, endian_conv=true)
    v = raw ? value : byte_array(value)
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_put, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt, Ptr{Ptr{UInt8}} ),
                 db, options,k, length(k), v, length(v), err)
    check_err(err)
end

function db_get(db, key; raw=false)
    # rocksdb_get will allocate the buffer for return value
    options = @threadcall( (:rocksdb_readoptions_create, librocksdb), Ptr{Nothing}, ())
    err = Ptr{UInt8}[0]
    val_len = Csize_t[0]
    k = byte_array(key, endian_conv=true)
    value = @threadcall( (:rocksdb_get, librocksdb), Ptr{UInt8},
          (Ptr{Nothing}, Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{Csize_t},  Ptr{Ptr{UInt8}} ),
          db, options, k, length(k), val_len, err)
    check_err(err)

    s = unsafe_wrap(Array, value, (val_len[1],), own=true)
    return val_len[1] == 0 ? nothing : (raw ? s : array_to_type(s))
end

function db_delete(db, key)
    options = @threadcall( (:rocksdb_writeoptions_create, librocksdb), Ptr{Nothing}, ())
    err = Ptr{UInt8}[0]
    k = byte_array(key, endian_conv=true)
    @threadcall( (:rocksdb_delete, librocksdb), Nothing,
          (Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, UInt, Ptr{Ptr{UInt8}} ),
          db, options, k, length(k), err)
    check_err(err)
end

function create_write_batch()
    batch = @threadcall( (:rocksdb_writebatch_create, librocksdb), Ptr{Nothing},())
    return batch
end

function batch_put(batch, key, value; raw=false)
    k = byte_array(key, endian_conv=true)
    v = raw ? value : byte_array(value)
    @threadcall( (:rocksdb_writebatch_put, librocksdb), Nothing,
          (Ptr{UInt8}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt),
          batch, k, length(k), v, length(v))
end

function write_batch(db, batch)
    options = @threadcall( (:rocksdb_writeoptions_create, librocksdb), Ptr{Nothing}, ())
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_write, librocksdb), Nothing,
          (Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing},  Ptr{Ptr{UInt8}} ),
          db, options, batch, err)
    check_err(err)
end

function create_iter(db::Ptr{Nothing}, options::Ptr{Nothing})
    @threadcall( (:rocksdb_create_iterator, librocksdb), Ptr{Nothing},
                 (Ptr{Nothing}, Ptr{Nothing}),
                 db, options)
end

function iter_valid(it::Ptr{Nothing})
    @threadcall( (:rocksdb_iter_valid, librocksdb), UInt8,
                 (Ptr{Nothing},),
                 it) != 0
end

function iter_destroy(it::Ptr{Nothing})
    @threadcall( (:rocksdb_iter_destroy, librocksdb), UInt8,
                 (Ptr{Nothing},),
                 it) != 0
end

function iter_key(it::Ptr{Nothing})
  k_len = Csize_t[0]
  key = @threadcall( (:rocksdb_iter_key, librocksdb), Ptr{UInt8},
    (Ptr{Nothing}, Ptr{Csize_t}),
    it, k_len)
  r = unsafe_wrap(Array, key, (k_len[1],), own=false)
end

function iter_value(it::Ptr{Nothing})
  v_len = Csize_t[0]
  value = @threadcall( (:rocksdb_iter_value, librocksdb), Ptr{UInt8},
    (Ptr{Nothing}, Ptr{Csize_t}),
    it, v_len)
  r = unsafe_wrap(Array, value, (v_len[1],), own=false)
end

function iter_seek(it::Ptr{Nothing}, key)
  @threadcall( (:rocksdb_iter_seek, librocksdb), Nothing,
    (Ptr{Nothing}, Ptr{UInt8}, UInt),
    it, key, length(key))
end

function iter_next(it::Ptr{Nothing})
  @threadcall( (:rocksdb_iter_next, librocksdb), Nothing,
    (Ptr{Nothing},),
    it)
end

abstract type AbstractRange end

function range_close(range::AbstractRange)
    if !range.destroyed
        range.destroyed = true
        @threadcall( (:rocksdb_iter_destroy, librocksdb), Nothing,
                     (Ptr{Nothing},),
                     range.iter)
        @threadcall( (:rocksdb_readoptions_destroy, librocksdb), Nothing,
                     (Ptr{Nothing},),
                     range.options)
    end
end

mutable struct Range <: AbstractRange
    iter::Ptr{Nothing}    # RocksDB iterator
    options::Ptr{Nothing} # RocksDB options
    key_start::Vector{UInt8}  # Starting key
    key_end::Vector{UInt8}    # Ending key
    raw::Bool          # Read the value without de-serializing to julia type
    destroyed::Bool    # RocksDB iterator and options are destroyed
end

mutable struct KeyRange <: AbstractRange
    iter::Ptr{Nothing}
    options::Ptr{Nothing}
    key_start::Vector{UInt8}
    key_end::Vector{UInt8}
    destroyed::Bool
end

function db_range(db, key_start, key_end; raw=false)
    options = @threadcall( (:rocksdb_readoptions_create, librocksdb), Ptr{Nothing}, ())
    ks = byte_array(key_start, endian_conv=true)
    ke = byte_array(key_end, endian_conv=true)
    iter = create_iter(db, options)
    range = Range(iter, options, ks, ke, raw, false)
    finalizer(range_close, range)
    range
end

function Base.iterate(range::AbstractRange)
    iter_seek(range.iter, range.key_start)
    iterate(range, nothing)
end

function Base.iterate(range::Range, state)
    if !iter_valid(range.iter) # No more valid keys
        return nothing
    end
    r = iter_key(range.iter)
    k = length(r) == 0 ? nothing : array_to_type(r, endian_conv=true)
    r = iter_value(range.iter)
    v = length(r) == 0 ? nothing : (range.raw ? r : array_to_type(r))
    iter_next(range.iter)
    (k == nothing || k == array_to_type(range.key_end, endian_conv=true)) && return nothing
    ((k, v), nothing)
end

#=
function Base.start(range::AbstractRange)
    iter_seek(range.iter, range.key_start)
end

function Base.done(range::AbstractRange, state=nothing)
    if range.destroyed
        return true
    end
    isdone = iter_valid(range.iter) ? iter_key(range.iter) > range.key_end : true
    if isdone
        range_close(range)
    end
    isdone
end

function next(range::Range, state=nothing)
    r = iter_key(range.iter)
    k = length(r) == 0 ? nothing : array_to_type(r, endian_conv=true)
    r = iter_value(range.iter)
    v = length(r) == 0 ? nothing : (range.raw ? r : array_to_type(r))
    iter_next(range.iter)
    ((k, v), nothing)
end
=#

"""
    db_key_range(db, key_start, key_end)
A range that returns only keys, no values
"""
function db_key_range(db, key_start, key_end)
    options = @threadcall( (:rocksdb_readoptions_create, librocksdb), Ptr{Nothing}, ())
    ks = byte_array(key_start, endian_conv=true)
    ke = byte_array(key_end, endian_conv=true)
    iter = create_iter(db, options)
    KeyRange(iter, options, ks, ke, false)
end

#=
function next(range::KeyRange, state=nothing)
    r = iter_key(range.iter)
    k = length(r) == 0 ? nothing : array_to_type(r, endian_conv=true)
    iter_next(range.iter)
    (k, nothing)
end
=#

"""
    db_snap_key_range(db, snap, key_start, key_end)
A range that returns only keys, no values
"""
function db_snap_key_range(db, snap, key_start, key_end)
    options = @threadcall( (:rocksdb_readoptions_create, librocksdb), Ptr{Nothing}, ())
    @threadcall( (:rocksdb_readoptions_set_snapshot, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{Nothing}), options, snap)
    ks = byte_array(key_start, endian_conv=true)
    ke = byte_array(key_end, endian_conv=true)
    iter = create_iter(db, options)
    Range(iter, options, ks, ke, false)
end

"""
    db_compact_range(db, start_key, end_key)
Compact (free storage) in the given range start_key:end_key
Caution: This call will delete ranges in snapshots.
"""
function db_compact_range(db, skey, ekey)
    err = Ptr{UInt8}[0]
    sk = byte_array(skey, endian_conv=true)
    ek = byte_array(ekey, endian_conv=true)
    # Delete files in range to free up storage
    @threadcall( (:rocksdb_delete_file_in_range, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt, Ptr{Ptr{UInt8}} ),
                 db, sk, length(sk), ek, length(ek), err)
    check_err(err)
    # Compact range to free up storage
    @threadcall( (:rocksdb_compact_range, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{UInt8}, UInt, Ptr{UInt8}, UInt),
                 db, sk, length(sk), ek, length(ek))
end

"""
    db_delete_range(db, start_key, limit_key)
Delete keys in the range start_key:limit_key including start key and limit_key.
"""
function db_delete_range(db, start_key, limit_key)
    sk = byte_array(start_key, endian_conv=true)
    lk = byte_array(limit_key, endian_conv=true)
    for k in RocksDB.db_key_range(db, sk, lk)
        db_delete(db, k)
    end
    #compact_range(db, start_key, limit_key)
end

"""
    db_backup_open(backup_dir)
Returns an instance of backup engine which backs up in *backup_dir*.
"""
function db_backup_open(backup_dir)
    options = @threadcall( (:rocksdb_options_create, librocksdb), Ptr{Nothing}, ())
    err = Ptr{UInt8}[0]
    be = @threadcall( (:rocksdb_backup_engine_open, librocksdb), Ptr{Nothing},
                 (Ptr{Nothing}, Ptr{UInt8}, Ptr{Ptr{UInt8}}),
                 options, backup_dir, err)
    check_err(err)
    return be
end

"""
    db_backup_create(be, db)
Given the backup engine *be*, backs up the database *db* in the directory
opened by *be*.
"""
function db_backup_create(be, db)
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_backup_engine_create_new_backup, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{Nothing}, Ptr{Ptr{UInt8}}),
                 be, db, err)
    check_err(err)
end

function db_backup_close(be)
    @threadcall( (:rocksdb_backup_engine_close, librocksdb), Nothing, (Ptr{Nothing},), be)
end

"""
    db_backup_purge(be, num_to_keep)
Keep *num_to_keep* back ups and purge the older ones.
"""
function db_backup_purge(be, num_to_keep)
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_backup_engine_purge_old_backups, librocksdb), Nothing,
                 (Ptr{Nothing}, UInt, Ptr{Ptr{UInt8}}),
                 be, num_to_keep, err)
    check_err(err)
end

"""
    db_backup_restore(be, db_dir)
Restore the latest backup in *be* to *db_dir*. The restored db can be opened
using *db_open*.
"""
function db_backup_restore(be, db_dir)
    options = @threadcall( (:rocksdb_restore_options_create, librocksdb), Ptr{Nothing}, ())
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_backup_engine_restore_db_from_latest_backup, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}, Ptr{Nothing}, Ptr{Ptr{UInt8}}),
                 be, db_dir, db_dir, options, err)
    check_err(err)
    @threadcall( (:rocksdb_restore_options_destroy, librocksdb), Nothing,
                 (Ptr{Nothing},), options)
end

"""
    db_create_snapshot(db)
For the given *db* return a snapshot object *snap* which can be iterated on with
*db_snap_key_range*. Snapshots are lightweight in RocksDB due to its LSM origins.
Snapshots are transient because they cannot be recovered after restarting the
process.
"""
function db_create_snapshot(db)
    @threadcall( (:rocksdb_create_snapshot, librocksdb), Ptr{Nothing},
                 (Ptr{Nothing},), db)
end

"""
    db_release_snapshot(db, snap)
Release storage for the given snapshot *snap* associated with *db*.
"""
function db_release_snapshot(db, snap)
    @threadcall( (:rocksdb_create_snapshot, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{Nothing}), db, snap)
end

"""
    db_create_checkpoint(db, path)
For the given *db* create a checkpoint in *path*. *path* must be an absolute
path to a directory where the checkpoint will be created. The last component of
*path* should not exist when making this call. Checkpoints are persistent.
They can be opened for reads and writes, as a regular db would, by passing
the *path*.

Checkpoints are more lightweight than backups. While backups copy the contents
of the database checkpoints create links to the original files. Checkpoints are
deleted by deleting the directory.
"""
function db_create_checkpoint(db, path)
    err = Ptr{UInt8}[0]
    @threadcall( (:rocksdb_create_checkpoint, librocksdb), Nothing,
                 (Ptr{Nothing}, Ptr{UInt8}, Ptr{Ptr{UInt8}}),
                 db, path, err)
    check_err(err)
end


text:
function repr_get_obj_fields_display(obj)
    field_names = fieldnames(typeof(obj))
    if length(field_names) == 0
        theme = term_theme[]
        return RenderableText(
            "$obj{$(theme.repr_type_style)}::$(typeof(obj)){/$(theme.repr_type_style)}",
        )
        return nothing
    end
    field_types = map(f -> "::" * string(f), typeof(obj).types)
    _values = map(f -> getfield(obj, f), field_names)

    fields = map(
        ft -> RenderableText(
            apply_style(string(ft[1]), term_theme[].repr_accent_style) *
            apply_style(string(ft[2]), term_theme[].repr_type_style),
        ),
        zip(field_names, field_types),
    )

    values = []
    for val in _values
        val = truncate(string(val), 45)
        push!(values, RenderableText.(val; style = term_theme[].repr_values_style))
    end

    line = vLine(length(fields); style = term_theme[].repr_line_style)
    space = Spacer(1, length(fields))

    content = rvstack(fields...) * line * space * lvstack(values...)
    return content
end

"""
    typename(typedef::Expr)

Get the name of a type as an expression
"""
function typename(typedef::Expr)
    return if typedef.args[2] isa Symbol
        typedef.args[2]
    elseif typedef.args[2].args[1] isa Symbol
        typedef.args[2].args[1]
    elseif typedef.args[2].args[1].args[1] isa Symbol
        typedef.args[2].args[1].args[1]
    else
        error("Could not parse type-head from: $typedef")
    end
end

function repr_panel(
    obj,
    content,
    subtitle;
    width = console_width() - 10,
    justify = :center,
    kwargs...,
)
    p = Panel(
        content;
        fit = true,
        title = isnothing(obj) ? obj : escape_brackets(string(typeof(obj))),
        title_justify = :left,
        width = width,
        justify = justify,
        style = term_theme[].repr_panel_style,
        title_style = term_theme[].repr_name_style,
        padding = (2, 1, 1, 1),
        subtitle = subtitle,
        subtitle_justify = :right,
        kwargs...,
    )

    return p
end

function vec_elems2renderables(v::Union{Tuple,AbstractVector}, N, max_w)
    shortsting(x) = x isa AbstractRenderable ? info(x) : truncate(string(x), max_w)
    out = highlight.(shortsting.(v[1:N]))

    length(v) > N && push!(out, "⋮";)
    return out
end

function matrix2content(mtx::AbstractMatrix; max_w = 12, max_items = 100, max_D = 10)
    max_D = console_width() < 150 ? 5 : max_D
    N = min(max_items, size(mtx, 1))
    D = min(max_D, size(mtx, 2))

    columns = [vec_elems2renderables(mtx[:, i], N, max_w) for i in 1:D]
    counts = RenderableText.("(" .* string.(1:N) .* ")"; style = "dim")
    top_counts = RenderableText.("(" .* string.(1:D) .* ")"; style = "dim white bold")

    space1, space2 = Spacer(3, length(counts)), Spacer(2, length(counts))

    content = ("" / "" / rvstack(counts...)) * space1
    for i in 1:(D - 1)
        content *= cvstack(top_counts[i], "", lvstack(columns[i])) * space2
    end
    content *= cvstack(top_counts[end], "", lvstack(columns[end]))

    if D < size(mtx, 2)
        content *= "" / "" / vstack((" {bold}⋯{/bold}" for i in 1:N)...)
    end
    return content
end

function vec2content(vec::Union{Tuple,AbstractVector})
    max_w = 88
    max_items = 100
    N = min(max_items, length(vec))

    if N == 0
        return "{bright_blue}empty vector{/bright_blue}"
    end

    vec_items = vec_elems2renderables(vec, N, max_w)
    counts = "(" .* string.(1:length(vec_items)) .* ")"

    content = Table(
        [counts vec_items];
        show_header = false,
        columns_justify = [:right, :left],
        columns_style = ["dim", "default"],
        columns_widths = [12, 60],
        vpad = 0,
        hpad = 2,
        compact = true,
        box = :NONE,
    )

    # content = rvstack(counts...) * Spacer(3, length(counts)) * cvstack(vec_items)
    return content
end


text:
@plugin struct TemplateFiles <: Plugin
    files::Array{String} = []
    dests::Array{String} = []
    vars::Dict{String, String} = Dict()
end

view(p::TemplateFiles, t::Template, pkg::AbstractString) = merge(p.vars, Dict(
    "PKG" => pkg,
))

# validate() = ...

function hook(p::TemplateFiles, t::Template, pkg_dir::AbstractString)
    pkg = basename(pkg_dir)

    for (file, dest) in zip(p.files, p.dests)
        dest = joinpath(pkg_dir, dest)
        f = render_file(file, combined_view(p, t, pkg), tags(p))
        gen_file(dest, f)
    end
end


text:
function convert_arguments(LS::Type{<:LineSegments}, aab::AABox{2})
    return convert_arguments(LS, edges(aab))
end

function convert_arguments(LS::Type{<:LineSegments}, R::Vector{<:AABox})
    return convert_arguments(LS, vcat([convert_arguments(LS, aab)[1] for aab ∈ R]...))
end

function convert_arguments(LS::Type{<:LineSegments}, R::SVector{N, <:AABox}) where {N}
    return convert_arguments(LS, collect(R))
end

function convert_arguments(M::Type{<:GLMakieMesh}, aab::AABox{2})
    vertices = [coordinates(v) for v in vertices(aab)]
    faces = [1 2 3;
             3 4 1]
    return convert_arguments(M, vertices, faces)
end

function convert_arguments(M::Type{<:GLMakieMesh}, aab::AABox{3})
    return convert_arguments(M, faces(aab))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#================================================

Shared testing functions for AbstractSkipList subtypes

================================================#

using Random, SkipLists, Test
using Base.Iterators: partition
using SkipLists: AbstractSkipList

#=========================
AbstractSkipList{T,:List} subtype tests
=========================#

@generated function test_construct_list(::Type{L}) where {L <: AbstractSkipList}
    quote
        @testset "Construct new $($L)" begin
            list = $L{Int64}()
            @test height(list) == 1
            @test length(list) == 0
            @test eltype(list) == Int64

            buff = IOBuffer()
            show(buff, list)
            seekstart(buff)

            :(
                if L == :SkipList
                    @test String(read(buff)) == "SkipList{Int64}(length = 0, height = 1)"
                end
            )

            # An error should be raised if we attempt to construct a skip list in an
            # invalid mode
            @test_throws ErrorException $L{Int64,:Foo}()

            # Try constructing a set instead
            set = $L{String,:Set}()
            @test height(set) == 1
            @test length(set) == 0
            @test eltype(set) == String

            buff = IOBuffer()
            show(buff, set)
            seekstart(buff)

            :(
                if L == :SkipList
                    @test String(read(buff)) == "SkipListSet{String}(length = 0, height = 1)"
                end
            )
        end
    end
end

@generated function test_insert_into_list(::Type{L}) where {L <: AbstractSkipList}
    quote
        @testset "Insert into $($L)" begin
            # Insert sorted values
            list = $L{Int64}()
            result_success = true
            for ii = 1:20
                result = insert!(list, ii)
                result_success = result_success && (result == Some(ii))
            end

            @test collect(list) == collect(1:20)
            @test isa(collect(list), Vector{Int64})
            @test length(list) == 20
            @test result_success

            insert!(list, 1)

            @test length(list) == 21
            @test collect(list) == cat([1], 1:20; dims=1)

            # Insert shuffled values
            list = $L{Int64}()
            result_success = true
            for ii in shuffle(1:20)
                result = insert!(list, ii)
                result_success = result_success && (result == Some(ii))
            end

            @test collect(list) == collect(1:20)
            @test isa(collect(list), Vector{Int64})
            @test length(list) == 20
            @test result_success
        end
    end
end

@generated function test_iterate_over_list(::Type{L}) where {L <: AbstractSkipList}
    quote
        @testset "Iterate over $($L)" begin
            vals = shuffle(1:100)
            list = $L{Int64}()
            for ii in vals
                insert!(list, ii)
            end

            success = true
            for (x1, x2) in zip(sort(vals), list)
                success = success && x1 == x2
            end

            @test success
        end
    end
end

@generated function test_list_membership(::Type{L}) where {L <: AbstractSkipList}
    quote
        @testset "Test membership in $($L)" begin
            list = $L{Int64}()
            @test 1 ∉ list

            insert!(list, 1)
            @test 1 ∈ list
        end
    end
end

@generated function test_add_duplicate_elements_to_list(::Type{L}) where {L <: AbstractSkipList}
    quote
        @testset "Add duplicate elements to $($L)" begin
            list = $L{Int64}()
            for ii = 1:2
                insert!(list, 1)
                insert!(list, 2)
            end

            @test length(list) == 4
            @test collect(list) == [1, 1, 2, 2]
            @test 1 ∈ list && 2 ∈ list

            delete!(list, 1)
            delete!(list, 2)
            @test length(list) == 2
            @test collect(list) == [1, 2]
            @test 1 ∈ list && 2 ∈ list

            delete!(list, 1)
            delete!(list, 2)
            @test length(list) == 0
            @test collect(list) == []
        end
    end
end

@generated function test_delete_from_list(::Type{L}) where {L <: AbstractSkipList}
    quote
        @testset "Delete from $($L)" begin
            list = $L{Int64}()
            insert!(list, 1)
            insert!(list, 2)
            insert!(list, 3)

            @test delete!(list, 1) == Some(1)
            @test length(list) == 2
            @test 1 ∉ list
            @test collect(list) == collect(2:3)

            @test delete!(list, 2) == Some(2)
            @test length(list) == 1
            @test 2 ∉ list
            @test collect(list) == collect(3:3)

            @test delete!(list, 3) == Some(3)
            @test length(list) == 0
            @test 3 ∉ list

            @test delete!(list, 0) == nothing
            @test collect(list) == []
            @test length(list) == 0
        end
    end
end

#=========================
SkipListSet tests
=========================#

@generated function test_insert_into_skip_list_set(::Type{S}) where {S <: AbstractSkipList{T,:Set} where T}
    quote
        @testset "Insert into $($S)" begin
            set = $S{Int64}()
            result_success = true
            for ii = 1:10
                result = insert!(set, ii)
                result_success = result_success && (result == Some(ii))
            end

            @test length(set) == 10
            @test collect(set) == 1:10
            @test result_success

            # If we now try to insert a duplicate element into the set, it shouldn't
            # have any effect
            result_success = true
            for ii = shuffle(1:10)
                result = insert!(set, ii)
                result_success = result_success && (result == nothing)
            end

            @test length(set) == 10
            @test collect(set) == 1:10
            @test result_success
        end
    end
end

@generated function test_delete_from_skip_list_set(::Type{S}) where {S <: AbstractSkipList{T,:Set} where T}
    quote
        @testset "Remove from $($S)" begin
            set = $S{Int64}()
            orig = 1:100

            result_success = true
            for ii in shuffle(orig)
                # Insert every element twice. The first insertion should return
                # Some(val); the second insertion fails and thus should return Nothing.
                result_success = result_success && (insert!(set, ii) == Some(ii))
                result_success = result_success && (insert!(set, ii) == nothing)
            end

            @test length(set) == length(orig)
            @test collect(set) == sort(orig)
            @test result_success

            # Remove all of the even elements
            to_remove = filter(iseven, orig)
            remaining = filter(isodd, orig) |> sort
            result_success = true
            for ii in shuffle(to_remove)
                # Delete the element twice. The first delete should be successful and
                # return Some(ii). The second delete should fail and return Nothing.
                result_success = result_success && (delete!(set, ii) == Some(ii))
                result_success = result_success && (delete!(set, ii) == nothing)
            end

            @test length(set) == length(remaining)
            @test collect(set) == remaining
            @test result_success

            # Test membership of remaining elements
            success = true
            for ii in remaining
                success = success && ii ∈ set
            end
            @test success
        end
    end
end

@generated function test_mixed_insert_delete_from_skip_list_set(::Type{S}) where {S <: AbstractSkipList{T,:Set} where T}
    quote
        @testset "Mixed insertion / deletion from $($S)" begin
            set = $S{Int64}()
            N = 10_000
            vals = rand(Int64, 2N)

            for val in vals
                insert!(set, val)
            end

            @test length(set) == 2N
            @test collect(set) == sort(vals)

            # Delete the first half of the elements from the vals array, and
            # simultaneously insert new elements into the array
            new_vals = rand(Int64, N)
            vals_to_delete = vals[1:N]

            insert_ops = zip(new_vals, repeated(:insert))
            delete_ops = zip(vals_to_delete, repeated(:delete))
            ops = cat(collect(insert_ops), collect(delete_ops); dims=1)

            success = true
            for (val, op) in ops
                if !success
                    break
                end

                if op == :insert
                    insert!(set, val)
                    success = val ∈ set
                else
                    delete!(set, val)
                    success = val ∉ set
                end
            end

            expected_vals = cat(vals[N+1:end], new_vals; dims=1)

            @test length(set) == 2N
            @test collect(set) == sort(expected_vals)
            @test success
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 77311686-7a2b-11eb-2bc7-95b2a2211969
begin
  using Pkg
  Pkg.activate("MLJ_env", shared=true)

  # using MLJ
  using CSV
  using DataFrames
  using PlutoUI
  using Random
  using Test
  using Printf
  using Plots
end

# ╔═╡ e9f8e942-7b29-11eb-1d64-d7e7ff02159c
begin
  include("./utils.jl");
  include("./dt_utils.jl");
end

# ╔═╡ 5435ff66-7a2b-11eb-1474-b96dcad21315
md"""
## C03w05: Boosting a decision stump

The goal of this notebook is to implement our own boosting module.

  - Use SFrames to do some feature engineering.
  - Modify the decision trees to incorporate weights.
  - Implement Adaboost ensembling.
  - Use your implementation of Adaboost to train a boosted decision stump ensemble.
  - Evaluate the effect of boosting (adding more decision stumps) on performance of the model.
  - Explore the robustness of Adaboost to overfitting.
"""

# ╔═╡ 81b7a9e0-7cc8-11eb-1d64-d7e7ff02159c
const DF= Union{DataFrame, SubDataFrame}

# ╔═╡ 771130dc-7a2b-11eb-170f-934c11fceede
md"""
#### Load LendingClub dataset
"""

# ╔═╡ 76f62e40-7a2b-11eb-095e-f36e207f06a2
begin
  loans =  CSV.File(("../../ML_UW_Spec/C03/data/lending-club-data.csv");
  header=true) |> DataFrame;
  size(loans)
end

# ╔═╡ 76de16d4-7a2b-11eb-09c6-e36b6c893c1f
md"""
Like the previous assignment, we reassign the labels to have +1 for a safe loan, and -1 for a risky (bad) loan.
"""

# ╔═╡ 76c17c4a-7a2b-11eb-2c47-97106bd58f99
insertcols!(loans, :safe_loans => ifelse.(loans.bad_loans .== 0, 1, -1),
  makeunique=true);

# ╔═╡ 76a83684-7a2b-11eb-1961-cb133db8c164
md"""
#### Extracting the target and the feature columns

We will now repeat some of the feature processing steps that we saw in the previous assignment:

Next, we select four categorical features:
  1. grade of the loan
  2. the length of the loan term
  3. the home ownership status: own, mortgage, rent
  4. number of years of employment.
"""

# ╔═╡ 768d6b08-7a2b-11eb-388a-8763c8da0a51
begin
    const Features = [
    	:grade,           # grade of the loan
        :emp_length,      # number of years of employment
        :home_ownership,  # home_ownership status: own, mortgage or rent
        :term,            # the term of the loan
    ]

  Target = :safe_loans # prediction target (y) (+1 means safe, -1 is risky)

  # Extract the feature columns and target column
  select!(loans, [Features..., Target]);
  length(names(loans)), names(loans)
end

# ╔═╡ 7673d8e4-7a2b-11eb-02b4-2f837181b4e9
md"""
#### Subsample dataset to make sure classes are balanced

Just as we did in the previous assignment, we will undersample the larger class (safe loans) in order to balance out our dataset. This means we are throwing away many data points.
"""

# ╔═╡ 3773355e-7cc8-11eb-2890-7577dda5d4de
begin
  function loan_ratios(safe_, risky_; n=size(loans)[1])
      (size(safe_)[1] / n, size(risky_)[1] / n)
  end

  function partition_sr(df::DF; target=Target)
    """Partition between safe and risky loan"""
    (df[df[!, target] .== 1, :], df[df[!, target] .!= 1, :])  # was -1
  end

  function resample!(safe_loans::DF, risky_loans::DF)
    """
    Since there are fewer risky loans than safe loans, find the ratio of
    the sizes and use that percentage to downsample the safe loans.
    """
    ## => deps on sampling(...)
    n = size(safe_loans)[1]
    perc = size(risky_loans)[1] / n

    ## down sample safe loans
    safe_loans = safe_loans[sampling(n, perc; seed=42), :];
    [risky_loans; safe_loans]
  end

end

# ╔═╡ 542db682-7a2e-11eb-35ad-25c82d6919e0
begin
  safe_loans₀, risky_loans₀ = partition_sr(loans)
  p_safe_loans₀, p_risky_loans₀ = loan_ratios(safe_loans₀, risky_loans₀);

  with_terminal() do
    @printf("Number of safe loans  : %d\n", size(safe_loans₀)[1])
    @printf("Number of risky loans : %d\n", size(risky_loans₀)[1])
    @printf("Percentage of safe loans:  %1.2f%%\n", p_safe_loans₀ * 100.)
    @printf("Percentage of risky loans: %3.2f%%\n", p_risky_loans₀ * 100.)
  end
end

# ╔═╡ 6327985a-7cc8-11eb-071c-0fb7dd06e885
loans_data = resample!(safe_loans₀, risky_loans₀);

# ╔═╡ 7641d404-7a2b-11eb-042d-bbc2c500ac59
begin
  safe_loans₁, risky_loans₁ = partition_sr(loans_data)
  p_safe_loans₁, p_risky_loans₁ = loan_ratios(safe_loans₁, risky_loans₁;
    n=size(loans_data)[1]);

  with_terminal() do
    @printf("Number of safe loans  : %d\n", size(safe_loans₁)[1])
    @printf("Number of risky loans : %d\n", size(risky_loans₁)[1])
    @printf("Percentage of safe loans:  %1.2f%%\n", p_safe_loans₁ * 100.)
    @printf("Percentage of risky loans: %3.2f%%\n", p_risky_loans₁ * 100.)
  end
end

# ╔═╡ 760c171a-7a2b-11eb-1341-19ba5a1824a1
md"""
#### Transform categorical data into binary features

In this assignment, we will implement **binary decision trees** (decision trees for binary features, a specific case of categorical variables taking on two values, e.g., true/false). Since all of our features are currently categorical features, we want to turn them into binary features.

For instance, the `home_ownership` feature represents the home ownership status of the loanee, which is either `own`, `mortgage` or `rent`. For example, if a data point has the feature
```
   {'home_ownership': 'RENT'}
```
we want to turn this into three features:
```
 {
   'home_ownership = OWN'      : 0,
   'home_ownership = MORTGAGE' : 0,
   'home_ownership = RENT'     : 1
 }
```
"""

# ╔═╡ 75c1a64e-7a2b-11eb-21f6-55dd3a8ce164
begin
  hot_encode!(loans_data)
  loans_data[1:5, :]
end

# ╔═╡ 75a52c44-7a2b-11eb-2f02-2312f3a3dac9
begin
  features = names(loans_data) |> a_ -> Symbol.(a_);
  deleteat!(features,
    findall(x -> x == Target, features))
end

# ╔═╡ 75736a9c-7a2b-11eb-2278-ab37447cf22d
## Expect 25 features after hot-encoding
@test length(features) == 25

# ╔═╡ 116eb5f8-7a33-11eb-2c4c-5f937540f987
length(loans_data.grade_A)
# vs 46508 in original notebook

# ╔═╡ a5c56d68-7a32-11eb-2ef6-e58ed65dc561
md"""
**Checkpoint:** Make sure the following answers match up.
"""

# ╔═╡ a5a1cdf6-7a32-11eb-18dd-0dfee8eb885c
@test sum(loans_data.grade_A) == 6537
# vs 6422 in original notebook

# ╔═╡ a56e44e0-7a32-11eb-2587-6f5f8b8c21f3
md"""
#### Train-test split

We split the data into a train test split with 80% of the data in the training set and 20% of the data in the test set.
"""

# ╔═╡ 6fad8afe-7a33-11eb-3208-6b289c3ccd23
begin
  train_data, test_data = train_test_split(loans_data; split=0.8, seed=70);
  size(train_data), size(test_data)
    ## vs ((37224, 26), (9284, 26)) ...
end

# ╔═╡ 6f88e4f6-7a33-11eb-10cf-e3e8608e616b
md"""
### Weighted decision trees

Let's modify our decision tree code from Module 5 to support weighting of individual data points.

#### Weighted error definition

Consider a model with $N$ data points with:
  -  Predictions $\hat{y}_1 ... \hat{y}_n$
  -  Target $y_1 ... y_n$
  -  Data point weights $\alpha_1 ... \alpha_n$.

Then the **weighted error** is defined by:

$$\mathrm{E}(\mathbf{\alpha}, \mathbf{\hat{y}}) = \frac{\sum_{i=1}^{n} \alpha_i \times 1[y_i \neq \hat{y_i}]}{\sum_{i=1}^{n} \alpha_i}$$

where $1[y_i \neq \hat{y_i}]$ is an indicator function that is set to $1$ if $y_i \neq \hat{y_i}$.

#### Write a function to compute weight of mistakes

Write a function that calculates the weight of mistakes for making the "weighted-majority" predictions for a dataset. The function accepts two inputs:
  -  `labels_in_node`: Targets $y_1 ... y_n$
  -  `data_weights`: Data point weights $\alpha_1 ... \alpha_n$

We are interested in computing the (total) weight of mistakes, i.e.

$$\mathrm{WM}(\mathbf{\alpha}, \mathbf{\hat{y}}) = \sum_{i=1}^{n} \alpha_i \times 1[y_i \neq \hat{y_i}].$$

This quantity is analogous to the number of mistakes, except that each mistake now carries different weight. It is related to the weighted error in the following way:

$$\mathrm{E}(\mathbf{\alpha}, \mathbf{\hat{y}}) = \frac{\mathrm{WM}(\mathbf{\alpha}, \mathbf{\hat{y}})}{\sum_{i=1}^{n} \alpha_i}$$

The function `inter_node_weighted_mistakes` should first compute two weights:
   - weight of mistakes $\mathrm{WM}_{-1}$, when all predictions are
$\hat{y}_i = -1,\ i.e\ \mathrm{WM}(\mathbf{\alpha}, \mathbf{-1})$

   -  weight of mistakes $\mathrm{WM}_{+1}$ when all predictions are
$\hat{y}_i = +1,\ i.e\ \mbox{WM}(\mathbf{\alpha}, \mathbf{+1})$

 where $\mathbf{-1}$ and $\mathbf{+1}$ are vectors where all values are -1 and +1 respectively.

After computing $\mathrm{WM}_{-1}$ and $\mathrm{WM}_{+1}$, the function `inter_node_weighted_mistakes` should return the lower of the two weights of mistakes, along with the class associated with that weight.
"""

# ╔═╡ 49f07e5c-7b2b-11eb-34ff-9b3613ca036c
function inter_node_weighted_mistakes(labels_in_node::AbstractVector{T},
                                      data_weights) where {T <: Real}

  ## Sum the weights of all +1 entries ≡ Weight of mistakes for predicting all -1's
  w_mistakes_all_neg = sum(data_weights[labels_in_node .== 1])

  ## Sum the weights of all -1 entries ≡ Weight of mistakes for predicting all +1's
  w_mistakes_all_pos = sum(data_weights[labels_in_node .== -1])

  ## Return the tuple (weight, class_label) representing the lower of the two weights.
  ## class_label should be an integer of value +1 or -1. tie => (tot_weight_po, +1)
  return w_mistakes_all_pos ≤ w_mistakes_all_neg ? (w_mistakes_all_pos, 1) :
  (w_mistakes_all_neg, -1)
end

# ╔═╡ 0179604e-7ccc-11eb-37bc-7df8411fb645
md"""
**Checkpoint**: Test your `inter_node_weighted_mistakes` function, run the following cell:
"""

# ╔═╡ 0163185c-7ccc-11eb-30dd-677a45cd9920
begin
  ex1_labels = Int[-1, -1, 1, 1, 1]
  ex1_data_weights = Float64[1., 2., .5, 1., 1.]

  @test inter_node_weighted_mistakes(ex1_labels, ex1_data_weights) == (2.5, -1)
  @test inter_node_weighted_mistakes(ex1_labels, ex1_data_weights)[2] == -1
end

# ╔═╡ 0149a662-7ccc-11eb-071c-0fb7dd06e885
md"""
Recall that the **classification error** is defined as follows:

$$\mbox{classification error} = \frac{\mbox{\# mistakes}}{\mbox{\# all data points}}$$

**Quiz Question:** If we set the weights $\mathbf{\alpha} = 1$ for all data points, how is the weight of mistakes $\mbox{WM}(\mathbf{\alpha}, \mathbf{\hat{y}})$ related to the `classification error`?
"""

# ╔═╡ 0131e5be-7ccc-11eb-2890-7577dda5d4de
## Answer: N × classifcation Error, N num. of data points
begin
  nex_labels = Int[-1, -1, 1, 1, 1]
  nex_data_weights = Float64[1., 1., 1., 1., 1.]  ## All 1s
  inter_node_weighted_mistakes(nex_labels, nex_data_weights)
end

# ╔═╡ 074dee62-7ccd-11eb-3b98-bd0e15b4de11
md"""
#### Function to pick best feature to split on

We continue modifying our decision tree code from the earlier assignment to incorporate weighting of individual data points. The next step is to pick the best feature to split on.

The `best_splitting_feature` function is similar to the one from the earlier assignment with two minor modifications:

  1. The function `best_splitting_feature` should now accept an extra parameter `data_weights` to take account of weights of data points.

  2. Instead of computing the number of mistakes in the left and right side of the split, we compute the weight of mistakes for both sides, add up the two weights, and divide it by the total weight of the data.
"""

# ╔═╡ 0731dd4e-7ccd-11eb-12fa-4d42c67876b2
function best_splitting_feature(df::DF, features::Vector{Symbol},
                                target::Symbol, data_weights)

  best_feature, best_error = nothing, Base.Inf

  for f ∈ features
    l_split = df[df[!, f] .== 0, :]
    r_split = df[df[!, f] .== 1, :]

    ## Apply the same filtering to data_weights to create l_data_weights,
    ## r_data_weights
    l_data_weights = data_weights[df[!, f] .== 0, :]
    r_data_weights = data_weights[df[!, f] .== 1, :]

    ## Calculate the weight of mistakes for left and right sides
    l_w_mistakes, _ = inter_node_weighted_mistakes(l_split[!, target], l_data_weights)
    r_w_mistakes, _ = inter_node_weighted_mistakes(r_split[!, target], r_data_weights)

    # Compute weighted error
    s_1, s_2 = sum(l_data_weights), sum(r_data_weights)
    error = (l_w_mistakes + r_w_mistakes) / (s_1 + s_2)

    if error < best_error
      best_error = error
      best_feature = f
    end
  end

  best_feature
end

# ╔═╡ 0715ebb6-7ccd-11eb-20d9-87257af20d51
md"""
**Checkpoint**: Now, we have another checkpoint to make sure you are on the right track
"""

# ╔═╡ 06e27130-7ccd-11eb-0a51-851d7325ad47
begin
  tr_data_weights = fill(1.5, size(train_data)[1])

  @test string(best_splitting_feature(train_data, features, Target,
    tr_data_weights)) == "term_ 60 months"
end

# ╔═╡ 06b1c8aa-7ccd-11eb-08e3-635d3737c018
md"""
**Very Optional**. Relationship between weighted error and weight of mistakes

By definition, the weighted error is the weight of mistakes divided by the weight of all data points, so

$$\mathrm{E}(\mathbf{\alpha}, \mathbf{\hat{y}}) = \frac{\sum_{i=1}^{n} \alpha_i \times 1[y_i \neq \hat{y_i}]}{\sum_{i=1}^{n} \alpha_i} = \frac{\mathrm{WM}(\mathbf{\alpha}, \mathbf{\hat{y}})}{\sum_{i=1}^{n} \alpha_i}.$$

In the code above, we obtain $\mathrm{E}(\mathbf{\alpha}, \mathbf{\hat{y}})$ from the two weights of mistakes from both sides, $\mathrm{WM}(\mathbf{\alpha}_{\mathrm{left}}, \mathbf{\hat{y}}_{\mathrm{left}})$ and $\mathrm{WM}(\mathbf{\alpha}_{\mathrm{right}}, \mathbf{\hat{y}}_{\mathrm{right}})$. First, notice that the overall weight of mistakes $\mathrm{WM}(\mathbf{\alpha}, \mathbf{\hat{y}})$ can be broken into two weights of mistakes over either side of the split:

```math
\mathrm{WM}(\mathbf{\alpha}, \mathbf{\hat{y}})
= \sum_{i=1}^{n} \alpha_i \times 1[y_i \neq \hat{y_i}]
= \sum_{\mathrm{left}} \alpha_i \times 1[y_i \neq \hat{y_i}]
 + \sum_{\mathrm{right}} \alpha_i \times 1[y_i \neq \hat{y_i}]\\
= \mathrm{WM}(\mathbf{\alpha}_{\mathrm{left}}, \mathbf{\hat{y}}_{\mathrm{left}}) + \mathrm{WM}(\mathbf{\alpha}_{\mathrm{right}}, \mathbf{\hat{y}}_{\mathrm{right}})
```

We then divide through by the total weight of all data points to obtain $\mathrm{E}(\mathbf{\alpha}, \mathbf{\hat{y}})$:

```math
\mathrm{E}(\mathbf{\alpha}, \mathbf{\hat{y}})
= \frac{\mathrm{WM}(\mathbf{\alpha}_{\mathrm{left}}, \mathbf{\hat{y}}_{\mathrm{left}}) + \mathrm{WM}(\mathbf{\alpha}_{\mathrm{right}}, \mathbf{\hat{y}}_{\mathrm{right}})}{\sum_{i=1}^{n} \alpha_i}
```

"""

# ╔═╡ 0694786a-7ccd-11eb-0a69-8fbbd0859786
md"""
#### Building the tree

With the above functions implemented correctly, we are now ready to build our decision tree. Recall from the previous assignments that each node in the decision tree is represented as a dictionary which contains the following keys:

    {
       'is_leaf'            : True/False.
       'prediction'         : Prediction at the leaf node.
       'left'               : (dictionary corresponding to the left tree).
       'right'              : (dictionary corresponding to the right tree).
       'features_remaining' : List of features that are posible splits.
    }

Let us start with a function that creates a leaf node given a set of target values:
"""

# ╔═╡ 068154f6-7ccd-11eb-05df-8b16a1439d68
function create_leaf(target_values, data_weights;
                     splitting_feature=nothing, left=nothing,
                     right=nothing, is_leaf=true)

  ## Create a leaf node
  leaf = Dict{Symbol, Any}(
                           :is_leaf => is_leaf,
                           :splitting_feature => splitting_feature,
                           :left => left,
                           :right => right
                           )

  if is_leaf
    ## Compute weight of mistakes.
    _, best_class = inter_node_weighted_mistakes(target_values, data_weights)
    ## Store the predicted class (1 or -1) in leaf['prediction']
    leaf[:prediction] = best_class
  end

  leaf
end

# ╔═╡ 0100bf38-7ccc-11eb-34ff-9b3613ca036c
md"""
And now the function that learns a weighted decision tree recursively and implements 3 stopping conditions:

  1. All data points in a node are from the same class.
  2. No more features to split on.
  3. Stop growing the tree when the tree depth reaches max_depth.

"""

# ╔═╡ b6fc15ae-7a37-11eb-071c-0fb7dd06e885
function weighted_decision_tree_create(df::DF, features, target, data_weights; 				curr_depth=1, max_depth=10, verbose=false, ϵ=1e-15)

  rem_features = copy(features)
  target_vals = df[!, target]

  if verbose
    println("-------------------------------------------------------------")
    @printf("Subtree, depth = %s (%s data points).\n", curr_depth,
            size(target_vals)[1])
  end

  if inter_node_weighted_mistakes(target_vals, data_weights)[1] ≤ ϵ
    ## Stopping cond1: no mistakes at current node?
    verbose && @printf("Stopping condition 1 reached.")
    return create_leaf(target_vals, data_weights)

  elseif length(rem_features) == 0
    ## Stopping cond2: no remaining features to consider splitting on?
    verbose && @printf("Stopping condition 2 reached.")
    return create_leaf(target_vals, data_weights)

  elseif curr_depth > max_depth
    ## Additional stopping condition (limit tree depth)
    verbose && @printf("Early stopping condition 1 reached. Maximum depth.")
    return create_leaf(target_vals, data_weights)
  end

  ## Now Find the best splitting feature
  splitting_feature = best_splitting_feature(df, features, target, data_weights)

  # Split on the best feature that we found.
  l_split = df[df[!, splitting_feature] .== 0, :]
  r_split = df[df[!, splitting_feature] .== 1, :]
  l_dweights = data_weights[df[!, splitting_feature] .== 0, :]
  r_dweights = data_weights[df[!, splitting_feature] .== 1, :]

  deleteat!(rem_features,
            findall(x -> x == splitting_feature, rem_features))

  verbose && @printf("Split on feature %s. (%s, %s)\n", splitting_feature,
                     size(l_split)[1], size(r_split)[1])

  # Create a leaf node if the split is "perfect"
  if size(l_split)[1] == size(df)[1]
    verbose && println("Creating leaf node.")
    return create_leaf(l_split[!, target], data_weights)
  end

  if size(r_split)[1] == size(df)[1]
    verbose && println("Creating leaf node.")
    return create_leaf(r_split[!, target], data_weights)
  end

  # Recurse on left and right subtrees
  ltree = weighted_decision_tree_create(l_split, rem_features, target,
		l_dweights; curr_depth=curr_depth + 1, max_depth, verbose, ϵ)

  rtree = weighted_decision_tree_create(r_split, rem_features, target,
		r_dweights; curr_depth=curr_depth + 1, max_depth, verbose, ϵ)

  return create_leaf(target_vals, data_weights;
                     splitting_feature=string(splitting_feature),
                     left=ltree,
                     right=rtree,
                     is_leaf=false)
end

# ╔═╡ cdb69d9c-7b2b-11eb-12fa-4d42c67876b2
md"""

Run the following test code to check the implementation.
"""

# ╔═╡ cd9e2b4a-7b2b-11eb-20d9-87257af20d51
begin
  small_data_weights = ones(Float64, size(train_data)[1])

  small_dt = weighted_decision_tree_create(train_data, features, Target,      small_data_weights; max_depth=2, verbose=true)

  @test count_nodes(small_dt) == 7
end

# ╔═╡ dfb94004-7d62-11eb-071c-0fb7dd06e885
small_dt

# ╔═╡ c872733c-7b34-11eb-05df-8b16a1439d68
md"""
#### Making Predictions with a weighted decision tree

we will use our `classify` function (cf. include).

### Evaluating the tree

Now, we will write a function to evaluate a decision tree by computing the classification error of the tree on the given dataset.

Again, recall that the **classification error** is defined as follows:
$$\mbox{classification error} = \frac{\mbox{\# mistakes}}{\mbox{\# all data points}}$$

we will use our `eval_classification_error` function (cf. include).
"""

# ╔═╡ 75fc81fc-7d63-11eb-0a69-8fbbd0859786
round(eval_classification_error(small_dt, test_data, Target),
  digits=4)

# ╔═╡ 75da5262-7d63-11eb-05df-8b16a1439d68
md"""
#### Example: Training a weighted decision tree

To build intuition on how weighted data points affect the tree being built, consider the following:

Suppose we only care about making good predictions for the **first 10 and last 10 items** in `train_data`, we assign weights:
  - 1 to the last 10 items
  - 1 to the first 10 items
  - and 0 to the rest.

Let us fit a weighted decision tree with `max_depth = 2`.
"""

# ╔═╡ 75c205c2-7d63-11eb-37bc-7df8411fb645
begin
  	## Assign weights
  	ex_data_weights = Float64[ones(Float64, 10)...,
		zeros(Float64, size(train_data)[1] - 20)..., ones(Float64, 10)...]

	## Train a weighted decision tree model.
	small_data_dt_subset20 = weighted_decision_tree_create(train_data, features,
		Target, ex_data_weights, max_depth=2)
end

# ╔═╡ 75a67960-7d63-11eb-30dd-677a45cd9920
md"""
Now, we will compute the classification error on the subset20, i.e. the subset of data points whose weight is 1 (namely the first and last 10 data points).
"""

# ╔═╡ a56f646c-7d64-11eb-071c-0fb7dd06e885
begin
  subset20 = vcat(first(train_data, 10), last(train_data, 10))
  round(eval_classification_error(small_data_dt_subset20, subset20, Target),
    digits=4)
end

# ╔═╡ a5394d28-7d64-11eb-34ff-9b3613ca036c
md"""
Now, let us compare the classification error of the model `small_data_dt_subset20` on the entire train\_data set:
"""

# ╔═╡ a5213ba2-7d64-11eb-1d64-d7e7ff02159c
round(eval_classification_error(small_data_dt_subset20, train_data, Target),
  digits=4)

# ╔═╡ 1892ec06-7d66-11eb-37bc-7df8411fb645
md"""
The model `small_data_dt_subset20` performs **a lot** better on `subset20` than on `train_data`.

So, what does this mean?
* The points with higher weights are the ones that are more important during the training process of the weighted decision tree.
* The points with zero weights are basically ignored during training.

**Quiz Question**: Will you get the same model as `small_data_dt_subset20` if you trained a decision tree with only the 20 data points with non-zero weights from the set of points in `subset_20`?
  - Yes
"""

# ╔═╡ a57ffa0a-7d66-11eb-05df-8b16a1439d68
md"""
### Implementing our own Adaboost (on decision stumps)

Now that we have a weighted decision tree working, it takes only a bit of work to implement Adaboost. For the sake of simplicity, let us stick with **decision tree stumps** by training trees with **`max_depth=1`**.

Recall from the lecture the procedure for Adaboost:

1\. Start with unweighted data with $\alpha_j = 1$

2\. For t = 1,...T:
  * Learn $f_t(x)$ with data weights $\alpha_j$
  * Compute coefficient $\hat{w}_t$:
     $$\hat{w}_t = \frac{1}{2}\ln{\left(\frac{1- \mbox{E}(\mathbf{\alpha}, \mathbf{\hat{y}})}{\mbox{E}(\mathbf{\alpha}, \mathbf{\hat{y}})}\right)}$$

  * Re-compute weights $\alpha_j$:
     $$\alpha_j \gets \begin{cases}
     \alpha_j \exp{(-\hat{w}_t)} & \text{ if }f_t(x_j) = y_j\\
     \alpha_j \exp{(\hat{w}_t)} & \text{ if }f_t(x_j) \neq y_j
     \end{cases}$$

  * Normalize weights $\alpha_j$:
      $$\alpha_j \gets \frac{\alpha_j}{\sum_{i=1}^{N}{\alpha_i}}$$

Write the function `adaboost_with_tree_stumps` (given Python skeleton).

"""

# ╔═╡ f962b6f0-7d66-11eb-0a69-8fbbd0859786
function adaboost_with_tree_stumps(df, features, target;
                                   num_tree_stumps=2, verbose=false)
  α = ones(Float64, size(df)[1])  ## start with unweighted data
  weights, tree_stumps = Float64[], []
  target_vals = select(df, target)[:, target]

  for _t ∈ 1:num_tree_stumps
    if verbose
      println("=====================================================")
      @printf("Adaboost Iteration %d\n", _t)
      println("=====================================================")
    end

    ## Learn a weighted decision tree stump. Use max_depth=1
    tree_stump = weighted_decision_tree_create(df, features, target, α;
                                               max_depth=1, verbose=verbose)
    push!(tree_stumps, tree_stump)
    preds = map(x -> classify(tree_stump, x), eachrow(df))

    ## Boolean array indicating whether each data point was correctly classified
    is_correct = preds .== target_vals
    is_wrong   = preds .≠ target_vals

    weighted_err = sum(α[is_wrong]) / sum(α)  ## Compute weighted error

    ## Compute model coefficient using weighted error
    weight = .5 * log((1 - weighted_err) / weighted_err)
    push!(weights, weight)

    ## Adjust weights on data point
    adjustment = ifelse.(is_correct, exp(-weight), exp(weight))

    ## Scale α by: × by adjustment, then normalize data points weights
    α = α .* adjustment
    α /= sum(α)
  end
  weights, tree_stumps
end

# ╔═╡ a55ef7c4-7d66-11eb-37bc-7df8411fb645
md"""
### Checking our Adaboost code

Train an ensemble of **two** tree stumps and see which features those stumps split on. We will run the algorithm with the following parameters:
  - `train_data`
  - `features`
  - `target`
  - `num_tree_stumps = 2`
"""

# ╔═╡ a516c486-7d66-11eb-30dd-677a45cd9920
stump_weights, tree_stumps = adaboost_with_tree_stumps(train_data, features, Target;	num_tree_stumps=2, verbose=true)

# ╔═╡ a4f7d690-7d66-11eb-071c-0fb7dd06e885
function  print_stump(tree, name=:root)
    split_name = tree[:splitting_feature] # ex. 'term. 36 months'

    if isnothing(split_name)
        print("(leaf, label: $(tree[:prediction]))")
        return nothing
  end
    println("                       $(name)")
    println("         |---------------|----------------|")
    println("         |                                |")
    println("         |                                |")
    println("         |                                |")
    println("     [$(split_name) == 0]              [$(split_name) == 1]")
    println("         |                                |")
    println("         |                                |")
    println("         |                                |")
    @printf("     (%s)                     (%s)\n",
        tree[:left][:is_leaf] ?
      string("leaf, label: ", string(tree[:left][:prediction])) : "subtree",
    tree[:right][:is_leaf] ?
          string("leaf, label: ", string(tree[:right][:prediction])) : "subtree")
end

# ╔═╡ a4dd0c0a-7d66-11eb-2890-7577dda5d4de
with_terminal() do
  print_stump(tree_stumps[1])
end

# ╔═╡ a4c39572-7d66-11eb-34ff-9b3613ca036c
with_terminal() do
  print_stump(tree_stumps[2])
end

# ╔═╡ a4aa3a96-7d66-11eb-1d64-d7e7ff02159c
md"""
If our Adaboost is correctly implemented, the following things should be true:

  - `tree_stumps[1]` should split on **term. 60 months** with the prediction -1 on the left and +1 on the right.
  - `tree_stumps[2]` should split on **grade.A** with the prediction -1 on the left and +1 on the right.
  - Weights should be approximately `[0.151, 0.187]`

**Reminders**
  - Stump weights ($\mathbf{\hat{w}}$) and data point weights ($\mathbf{\alpha}$) are two different concepts.
  - Stump weights ($\mathbf{\hat{w}}$) tell you how important each stump is while making predictions with the entire boosted ensemble.
  - Data point weights ($\mathbf{\alpha}$) tell you how important each data point is while training a decision stump.
"""

# ╔═╡ 28bccc2c-7d80-11eb-2c80-6b7e9646f76e
md"""
### Training a boosted ensemble of 10 stumps

Let us train an ensemble of 10 decision tree stumps with Adaboost. We run the `adaboost_with_tree_stumps` function with the following parameters:
  - `train_data`
  - `features`
  - `target`
  - `num_tree_stumps = 10`
"""

# ╔═╡ 4653cf88-7d80-11eb-354e-f54c38b7781c
stump_weights10, tree_stumps10 = adaboost_with_tree_stumps(train_data, features, Target,
                                                           num_tree_stumps=10)

# ╔═╡ 28a27e62-7d80-11eb-2d86-2ff7694e596e
md"""
#### Making predictions

Recall from the lecture that in order to make predictions, we use the following formula:

$$\hat{y} = sign\left(\sum_{t=1}^T \hat{w}_t f_t(x)\right)$$

We need to do the following things:
- Compute the predictions $f_t(x)$ using the $t$-th decision tree
- Compute $\hat{w}_t f_t(x)$ by multiplying the `stump_weights` with the predictions $f_t(x)$ from the decision trees
- Sum the weighted predictions over each stump in the ensemble.

Complete the following skeleton for making predictions:
"""

# ╔═╡ 2885d134-7d80-11eb-3039-85b1d81f2b21
function predict_adaboost(stump_weights, tree_stumps, df)

  scores = zeros(Float64, size(df)[1])

  for (ix, tree_stump) ∈ enumerate(tree_stumps)
	preds = map(x -> classify(tree_stump, x), eachrow(df))
    scores += stump_weights[ix] .* preds
  end

  ifelse.(scores .> zero(eltype(scores)), 1, -1)
end

# ╔═╡ 69dc69d2-7d81-11eb-2afd-216aca94000f
begin
  preds10 = predict_adaboost(stump_weights, tree_stumps, test_data)
  acc10 = sum(test_data[!, Target] .== preds10) / size(test_data)[1]

  with_terminal() do
    @printf("\nAccuracy of 10-component ensemble = %1.4f\n", acc10)
  end
end

# ╔═╡ 4205b9e0-7d81-11eb-283e-174f221087c0
stump_weights10

# ╔═╡ 525dffa0-7d81-11eb-346f-a78d0b34253e
md"""
**Quiz Question:** Are the weights monotonically decreasing, monotonically increasing, or neither?
  - Neither (from iter1 to iter 2: increase, after that decrease)

**Reminder**: Stump weights ($\mathbf{\hat{w}}$) tell you how important each stump is while making predictions with the entire boosted ensemble.
"""

# ╔═╡ 4214e72a-7d82-11eb-1d64-d7e7ff02159c
md"""
#### Performance plots

In this section, we will try to reproduce some of the performance plots dicussed in the lecture.

### How does accuracy change with adding stumps to the ensemble?

We will now train an ensemble with:
  - `train_data`
  - `features`
  - `target`
  - `num_tree_stumps = 30`

Once we are done with this, we will then do the following:
  - Compute the classification error at the end of each iteration.
  - Plot a curve of classification error vs iteration.

First, let's train the model.
"""

# ╔═╡ 62a40a70-7d82-11eb-34ff-9b3613ca036c
stump_weights30, tree_stumps30 = adaboost_with_tree_stumps(train_data, features, Target, num_tree_stumps=30)

# ╔═╡ 72b17ef2-7d82-11eb-30dd-677a45cd9920
md"""
##### Computing training error at the end of each iteration

Now, we will compute the classification error on the `train_data` and see how it is reduced as trees are added.
"""

# ╔═╡ 728ff304-7d82-11eb-071c-0fb7dd06e885
function calc_error(df, label, stump_weights, tree_stumps;
		target=Target, n=30)
    err_all = Float64[]

    for ix ∈ 1:n
      preds = predict_adaboost(stump_weights[1:ix], tree_stumps[1:ix], df)
      err = 1.0 - sum(df[!, target] .== preds) / size(df)[1]
      push!(err_all, err)
      @printf("Iteration %2d, %s error = %2.5f\n", ix, label, err_all[ix])
  end

  err_all
end

# ╔═╡ 72750690-7d82-11eb-2890-7577dda5d4de
error_all = calc_error(train_data, :training, stump_weights30, tree_stumps30)

# ╔═╡ 56fc7b38-7d84-11eb-05df-8b16a1439d68
test_error_all = calc_error(test_data, :test, stump_weights30, tree_stumps30)

# ╔═╡ 05a54422-7d84-11eb-37bc-7df8411fb645
begin
  plot(collect(1:30), error_all, linewidth=3.0, label="Training error")
  plot!(collect(1:30), test_error_all, linewidth=3.0, label="Testing error")
end

# ╔═╡ Cell order:
# ╟─5435ff66-7a2b-11eb-1474-b96dcad21315
# ╠═77311686-7a2b-11eb-2bc7-95b2a2211969
# ╠═81b7a9e0-7cc8-11eb-1d64-d7e7ff02159c
# ╠═e9f8e942-7b29-11eb-1d64-d7e7ff02159c
# ╟─771130dc-7a2b-11eb-170f-934c11fceede
# ╠═76f62e40-7a2b-11eb-095e-f36e207f06a2
# ╟─76de16d4-7a2b-11eb-09c6-e36b6c893c1f
# ╠═76c17c4a-7a2b-11eb-2c47-97106bd58f99
# ╟─76a83684-7a2b-11eb-1961-cb133db8c164
# ╠═768d6b08-7a2b-11eb-388a-8763c8da0a51
# ╟─7673d8e4-7a2b-11eb-02b4-2f837181b4e9
# ╠═3773355e-7cc8-11eb-2890-7577dda5d4de
# ╠═542db682-7a2e-11eb-35ad-25c82d6919e0
# ╠═6327985a-7cc8-11eb-071c-0fb7dd06e885
# ╠═7641d404-7a2b-11eb-042d-bbc2c500ac59
# ╟─760c171a-7a2b-11eb-1341-19ba5a1824a1
# ╠═75c1a64e-7a2b-11eb-21f6-55dd3a8ce164
# ╠═75a52c44-7a2b-11eb-2f02-2312f3a3dac9
# ╠═75736a9c-7a2b-11eb-2278-ab37447cf22d
# ╠═116eb5f8-7a33-11eb-2c4c-5f937540f987
# ╟─a5c56d68-7a32-11eb-2ef6-e58ed65dc561
# ╠═a5a1cdf6-7a32-11eb-18dd-0dfee8eb885c
# ╟─a56e44e0-7a32-11eb-2587-6f5f8b8c21f3
# ╠═6fad8afe-7a33-11eb-3208-6b289c3ccd23
# ╟─6f88e4f6-7a33-11eb-10cf-e3e8608e616b
# ╠═49f07e5c-7b2b-11eb-34ff-9b3613ca036c
# ╟─0179604e-7ccc-11eb-37bc-7df8411fb645
# ╠═0163185c-7ccc-11eb-30dd-677a45cd9920
# ╟─0149a662-7ccc-11eb-071c-0fb7dd06e885
# ╠═0131e5be-7ccc-11eb-2890-7577dda5d4de
# ╟─074dee62-7ccd-11eb-3b98-bd0e15b4de11
# ╠═0731dd4e-7ccd-11eb-12fa-4d42c67876b2
# ╟─0715ebb6-7ccd-11eb-20d9-87257af20d51
# ╠═06e27130-7ccd-11eb-0a51-851d7325ad47
# ╟─06b1c8aa-7ccd-11eb-08e3-635d3737c018
# ╟─0694786a-7ccd-11eb-0a69-8fbbd0859786
# ╠═068154f6-7ccd-11eb-05df-8b16a1439d68
# ╟─0100bf38-7ccc-11eb-34ff-9b3613ca036c
# ╠═b6fc15ae-7a37-11eb-071c-0fb7dd06e885
# ╟─cdb69d9c-7b2b-11eb-12fa-4d42c67876b2
# ╠═cd9e2b4a-7b2b-11eb-20d9-87257af20d51
# ╠═dfb94004-7d62-11eb-071c-0fb7dd06e885
# ╟─c872733c-7b34-11eb-05df-8b16a1439d68
# ╠═75fc81fc-7d63-11eb-0a69-8fbbd0859786
# ╟─75da5262-7d63-11eb-05df-8b16a1439d68
# ╠═75c205c2-7d63-11eb-37bc-7df8411fb645
# ╟─75a67960-7d63-11eb-30dd-677a45cd9920
# ╠═a56f646c-7d64-11eb-071c-0fb7dd06e885
# ╟─a5394d28-7d64-11eb-34ff-9b3613ca036c
# ╠═a5213ba2-7d64-11eb-1d64-d7e7ff02159c
# ╟─1892ec06-7d66-11eb-37bc-7df8411fb645
# ╟─a57ffa0a-7d66-11eb-05df-8b16a1439d68
# ╠═f962b6f0-7d66-11eb-0a69-8fbbd0859786
# ╟─a55ef7c4-7d66-11eb-37bc-7df8411fb645
# ╠═a516c486-7d66-11eb-30dd-677a45cd9920
# ╠═a4f7d690-7d66-11eb-071c-0fb7dd06e885
# ╠═a4dd0c0a-7d66-11eb-2890-7577dda5d4de
# ╠═a4c39572-7d66-11eb-34ff-9b3613ca036c
# ╟─a4aa3a96-7d66-11eb-1d64-d7e7ff02159c
# ╟─28bccc2c-7d80-11eb-2c80-6b7e9646f76e
# ╠═4653cf88-7d80-11eb-354e-f54c38b7781c
# ╟─28a27e62-7d80-11eb-2d86-2ff7694e596e
# ╠═2885d134-7d80-11eb-3039-85b1d81f2b21
# ╠═69dc69d2-7d81-11eb-2afd-216aca94000f
# ╠═4205b9e0-7d81-11eb-283e-174f221087c0
# ╟─525dffa0-7d81-11eb-346f-a78d0b34253e
# ╟─4214e72a-7d82-11eb-1d64-d7e7ff02159c
# ╠═62a40a70-7d82-11eb-34ff-9b3613ca036c
# ╟─72b17ef2-7d82-11eb-30dd-677a45cd9920
# ╠═728ff304-7d82-11eb-071c-0fb7dd06e885
# ╠═72750690-7d82-11eb-2890-7577dda5d4de
# ╠═56fc7b38-7d84-11eb-05df-8b16a1439d68
# ╠═05a54422-7d84-11eb-37bc-7df8411fb645