text:
"""
$(TYPEDEF)

Interpolant-basis weighted sum-of-squares polynomial (of dimension `U`) epigraph
of `ℓ₂` norm (of dimension `R`) cone, parametrized by vector of matrices `Ps`
derived from interpolant basis and polynomial domain constraints.

    $(FUNCTIONNAME){T}(R::Int, U::Int, Ps::Vector{Matrix{T}}, use_dual::Bool = false)
"""
mutable struct WSOSInterpEpiNormEucl{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    R::Int
    U::Int
    Ps::Vector{Matrix{T}}
    nu::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    mat::Vector{Matrix{T}}
    matfact::Vector
    ΛLi_Λ::Vector{Vector{Matrix{T}}}
    Λ11::Vector{Matrix{T}}
    tempLL::Vector{Matrix{T}}
    tempLU::Vector{Matrix{T}}
    tempLU2::Vector{Matrix{T}}
    tempLU_vec::Vector{Vector{Matrix{T}}}
    tempLU_vec2::Vector{Vector{Matrix{T}}}
    tempLU_vec3::Vector{Matrix{T}}
    tempLL_vec::Vector{Vector{Matrix{T}}}
    tempLL_vec2::Vector{Matrix{T}}
    tempLRUR::Vector{Matrix{T}}
    tempUU::Matrix{T}
    ΛLiPs_edge::Vector{Vector{Matrix{T}}}
    PΛiPs::Vector{Matrix{T}}
    PΛiP_blocks_U::Vector
    Λ11LiP::Vector{Matrix{T}} # also equal to the block on the diagonal of ΛLiP
    matLiP::Vector{Matrix{T}} # also equal to block (1, 1) of ΛLiP
    PΛ11iP::Vector{Matrix{T}}
    Λfact::Vector
    point_views::Vector
    Ps_times::Vector{Float64}
    Ps_order::Vector{Int}

    function WSOSInterpEpiNormEucl{T}(
        R::Int,
        U::Int,
        Ps::Vector{Matrix{T}};
        use_dual::Bool = false,
        ) where {T <: Real}
        for Pj in Ps
            @assert size(Pj, 1) == U
        end
        cone = new{T}()
        cone.use_dual_barrier = !use_dual # using dual barrier
        cone.dim = U * R
        cone.R = R
        cone.U = U
        cone.Ps = Ps
        cone.nu = 2 * sum(size(Pk, 2) for Pk in Ps)
        return cone
    end
end

function setup_extra_data!(cone::WSOSInterpEpiNormEucl{T}) where {T <: Real}
    U = cone.U
    R = cone.R
    Ps = cone.Ps
    K = length(Ps)
    Ls = [size(Pk, 2) for Pk in cone.Ps]
    cone.mat = [zeros(T, L, L) for L in Ls]
    cone.matfact = Vector{Any}(undef, K)
    cone.ΛLi_Λ = [[zeros(T, L, L) for _ in 1:(R - 1)] for L in Ls]
    cone.Λ11 = [zeros(T, L, L) for L in Ls]
    cone.tempLL = [zeros(T, L, L) for L in Ls]
    cone.tempLU = [zeros(T, L, U) for L in Ls]
    cone.tempLU2 = [zeros(T, L, U) for L in Ls]
    cone.tempLU_vec = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.tempLU_vec2 = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.tempLU_vec3 = [zeros(T, L, U * R) for L in Ls]
    cone.tempLL_vec = [[zeros(T, L, L) for _ in 1:(R - 1)] for L in Ls]
    cone.tempLL_vec2 = [zeros(R * L, L) for L in Ls]
    cone.tempLRUR = [zeros(T, L * R, U * R) for L in Ls]
    cone.tempUU = zeros(T, U, U)
    cone.ΛLiPs_edge = [[zeros(T, L, U) for _ in 1:(R - 1)] for L in Ls]
    cone.matLiP = [zeros(T, L, U) for L in Ls]
    cone.PΛiPs = [zeros(T, R * U, R * U) for _ in eachindex(Ls)]
    cone.Λ11LiP = [zeros(T, L, U) for L in Ls]
    cone.PΛ11iP = [zeros(T, U, U) for _ in eachindex(Ps)]
    cone.PΛiP_blocks_U = [[view(PΛiPk, block_idxs(U, r), block_idxs(U, s))
        for r in 1:R, s in 1:R] for PΛiPk in cone.PΛiPs]
    cone.Λfact = Vector{Any}(undef, K)
    cone.point_views = [view(cone.point, block_idxs(U, i)) for i in 1:R]
    cone.Ps_times = zeros(K)
    cone.Ps_order = collect(1:K)
    return cone
end

function set_initial_point!(arr::AbstractVector, cone::WSOSInterpEpiNormEucl)
    @views arr[1:cone.U] .= 1
    @views arr[(cone.U + 1):end] .= 0
    return arr
end

function update_feas(cone::WSOSInterpEpiNormEucl)
    @assert !cone.feas_updated
    U = cone.U
    Λfact = cone.Λfact
    matfact = cone.matfact
    point_views = cone.point_views

    # order the Ps by how long it takes to check feasibility, to improve efficiency
    sortperm!(cone.Ps_order, cone.Ps_times, initialized = true) # stochastic

    cone.is_feas = true
    for k in cone.Ps_order
        cone.Ps_times[k] = @elapsed @inbounds begin
            Pk = cone.Ps[k]
            Λ11k = cone.Λ11[k]
            LUk = cone.tempLU[k]
            ΛLi_Λ = cone.ΛLi_Λ[k]
            mat = cone.mat[k]

            # first lambda
            @. LUk = Pk' * point_views[1]'
            mul!(Λ11k, LUk, Pk)
            copyto!(mat, Λ11k)
            Λfact[k] = cholesky!(Symmetric(Λ11k, :U), check = false)
            if !isposdef(Λfact[k])
                cone.is_feas = false
                break
            end

            # subtract others
            uo = U + 1
            for r in 1:(cone.R - 1)
                @. LUk = Pk' * point_views[r + 1]'
                mul!(ΛLi_Λ[r], LUk, Pk)
                ldiv!(Λfact[k].L, ΛLi_Λ[r])
                mul!(mat, ΛLi_Λ[r]', ΛLi_Λ[r], -1, true)
                uo += U
            end

            matfact[k] = cholesky!(Symmetric(mat, :U), check = false)
            if !isposdef(matfact[k])
                cone.is_feas = false
                break
            end
        end
    end

    cone.feas_updated = true
    return cone.is_feas
end

function update_grad(cone::WSOSInterpEpiNormEucl)
    @assert cone.is_feas
    U = cone.U
    R = cone.R
    Λfact = cone.Λfact
    matfact = cone.matfact
    grad = cone.grad

    grad .= 0
    @inbounds for k in eachindex(cone.Ps)
        Pk = cone.Ps[k]
        Λ11LiP = cone.Λ11LiP[k]
        ΛLiP_edge = cone.ΛLiPs_edge[k]
        matLiP = cone.matLiP[k]
        ΛLi_Λ = cone.ΛLi_Λ[k]

        ldiv!(Λ11LiP, cone.Λfact[k].L, Pk') # TODO may be more efficient to do
        # ldiv(fact.U', B) than ldiv(fact.L, B) here and elsewhere since the factorizations are of symmetric :U matrices

        # prep PΛiP halves
        ldiv!(matLiP, matfact[k].L, Pk')
        # top edge of ΛLiP
        for r in 1:(R - 1)
            mul!(ΛLiP_edge[r], ΛLi_Λ[r]', Λ11LiP, -1, false)
            ldiv!(matfact[k].L, ΛLiP_edge[r])
        end

        @views for u in 1:U
            grad[u] -= sum(abs2, Λ11LiP[:, u]) + sum(abs2, matLiP[:, u])
            idx = U + u
            for r in 2:R
                grad[idx] -= 2 * dot(ΛLiP_edge[r - 1][:, u], matLiP[:, u])
                grad[u] -= sum(abs2, ΛLiP_edge[r - 1][:, u])
                idx += U
            end
        end
    end

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::WSOSInterpEpiNormEucl)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    U = cone.U
    R = cone.R
    R2 = R - 2
    UU = cone.tempUU
    matfact = cone.matfact

    H .= 0
    @inbounds for k in eachindex(cone.Ps)
        PΛiPs = cone.PΛiP_blocks_U[k]
        PΛ11iP = cone.PΛ11iP[k]
        ΛLiP_edge = cone.ΛLiPs_edge[k]
        matLiP = cone.matLiP[k]
        Λ11LiP = cone.Λ11LiP[k]

        # prep PΛiPs
        # P * inv(Λ_11) * P' for (1, 1) hessian block and adding to PΛiPs[r][r]
        mul!(PΛ11iP, Λ11LiP', Λ11LiP)
        # block-(1,1) is P * inv(mat) * P'
        mul!(PΛiPs[1, 1], matLiP', matLiP)
        # get all the PΛiPs that are in row one or on the diagonal
        for r in 2:R
            mul!(PΛiPs[r, 1], ΛLiP_edge[r - 1]', matLiP)
            mul!(PΛiPs[r, r], ΛLiP_edge[r - 1]', ΛLiP_edge[r - 1])
            @. PΛiPs[r, r] += PΛ11iP
            for r2 in 2:(r - 1)
                mul!(PΛiPs[r, r2], ΛLiP_edge[r - 1]', ΛLiP_edge[r2 - 1])
            end
        end
        copytri!(cone.PΛiPs[k], 'L')

        for i in 1:U, k in 1:i
            H[k, i] -= abs2(PΛ11iP[k, i]) * R2
        end

        @. @views H[1:U, 1:U] += abs2(PΛiPs[1, 1])
        for r in 2:R
            idxs = block_idxs(U, r)
            for s in 1:(r - 1)
                # block (1,1)
                @. UU = abs2(PΛiPs[r, s])
                @. @views H[1:U, 1:U] += UU
                @. @views H[1:U, 1:U] += UU'
                # blocks (1,r)
                @. @views H[1:U, idxs] += PΛiPs[s, 1] * PΛiPs[s, r]
            end
            # block (1,1)
            @. @views H[1:U, 1:U] += abs2(PΛiPs[r, r])
            # blocks (1,r)
            @. @views H[1:U, idxs] += PΛiPs[r, 1] * PΛiPs[r, r]
            # blocks (1,r)
            for s in (r + 1):R
                @. @views H[1:U, idxs] += PΛiPs[s, 1] * PΛiPs[s, r]
            end

            # blocks (r, r2)
            # for H[idxs, idxs], UU are symmetric
            @. UU = PΛiPs[r, 1] * PΛiPs[r, 1]'
            @. @views H[idxs, idxs] += UU
            @. UU = PΛiPs[1, 1] * PΛiPs[r, r]
            @. @views H[idxs, idxs] += UU
            for r2 in (r + 1):R
                idxs2 = block_idxs(U, r2)
                @. UU = PΛiPs[r, 1] * PΛiPs[r2, 1]'
                @. @views H[idxs, idxs2] += UU
                @. UU = PΛiPs[1, 1] * PΛiPs[r2, r]'
                @. @views H[idxs, idxs2] += UU
            end
        end
    end
    @. @views H[:, (U + 1):cone.dim] *= 2

    cone.hess_updated = true
    return cone.hess
end

function dder3(cone::WSOSInterpEpiNormEucl, dir::AbstractVector)
    @assert cone.hess_updated
    dder3 = cone.dder3
    dder3 .= 0
    R = cone.R
    U = cone.U

    @inbounds for k in eachindex(cone.Ps)
        Pk = cone.Ps[k]
        LUk = cone.tempLU[k]
        LLk = cone.tempLL[k]
        ΛFLPk = cone.Λ11LiP[k]

        @views mul!(LUk, ΛFLPk, Diagonal(dir[1:U]))
        mul!(LLk, LUk, ΛFLPk')
        mul!(LUk, Hermitian(LLk), ΛFLPk)
        @views for u in 1:U
            dder3[u] += sum(abs2, LUk[:, u])
        end
    end
    @. @views dder3[1:U] *= 2 - R

    @inbounds for k in eachindex(cone.Ps)
        L = size(cone.Ps[k], 2)
        ΛLiP_edge = cone.ΛLiPs_edge[k]
        matLiP = cone.matLiP[k]
        Λ11LiP = cone.Λ11LiP[k]
        scaled_row = cone.tempLU_vec[k]
        scaled_col = cone.tempLU_vec2[k]
        ΛLiP_edge_ctgs = cone.tempLU_vec3[k]
        ΛLiP_D_ΛLiPt_row = cone.tempLL_vec[k]
        ΛLiP_D_ΛLiPt_col = cone.tempLL_vec2[k]
        ΛLiP_D_ΛLiPt_diag = cone.tempLL[k]
        dder3_half = cone.tempLRUR[k]
        dder3_half .= 0

        # get ΛLiP * D * PΛiP where D is diagonalized dir scattered in an arrow
        # and ΛLiP is half an arrow
        # ΛLiP * D is an arrow matrix but row edge doesn't equal column edge
        @views scaled_diag = mul!(cone.tempLU[k], Λ11LiP, Diagonal(dir[1:U]))
        @views scaled_pt = mul!(cone.tempLU2[k], matLiP, Diagonal(dir[1:U]))
        @views for r in 2:R
            mul!(scaled_pt, ΛLiP_edge[r - 1], Diagonal(dir[block_idxs(U, r)]),
                true, true)
            mul!(scaled_row[r - 1], matLiP, Diagonal(dir[block_idxs(U, r)]))
            mul!(scaled_row[r - 1], ΛLiP_edge[r - 1], Diagonal(dir[1:U]),
                true, true)
            mul!(scaled_col[r - 1], Λ11LiP, Diagonal(dir[block_idxs(U, r)]))
        end

        # ΛLiP is half an arrow, multiplying arrow by its transpose gives arrow
        @views mul!(ΛLiP_D_ΛLiPt_col[1:L, :], scaled_pt, matLiP')
        @views for r in 2:R
            mul!(ΛLiP_D_ΛLiPt_col[1:L, :], scaled_row[r - 1], ΛLiP_edge[r - 1]',
                true, true)
            mul!(ΛLiP_D_ΛLiPt_row[r - 1], scaled_row[r - 1], Λ11LiP')
            mul!(ΛLiP_D_ΛLiPt_col[block_idxs(L, r), :], scaled_col[r - 1], matLiP')
            mul!(ΛLiP_D_ΛLiPt_col[block_idxs(L, r), :], scaled_diag,
                ΛLiP_edge[r - 1]', true, true)
            mul!(ΛLiP_D_ΛLiPt_diag, scaled_diag, Λ11LiP')
        end

        # TODO check that this is worthwhile
        @. @views ΛLiP_edge_ctgs[:, 1:U] = matLiP
        for r in 2:R
            @. @views ΛLiP_edge_ctgs[:, block_idxs(U, r)] = ΛLiP_edge[r - 1]
        end

        # multiply by ΛLiP
        mul!(dder3_half, ΛLiP_D_ΛLiPt_col, ΛLiP_edge_ctgs)
        @views for r in 2:R
            mul!(dder3_half[1:L, block_idxs(U, r)], ΛLiP_D_ΛLiPt_row[r - 1],
                Λ11LiP, true, true)
            mul!(dder3_half[block_idxs(L, r), block_idxs(U, r)], ΛLiP_D_ΛLiPt_diag,
                Λ11LiP, true, true)
        end

        @views for u in 1:U
            dder3[u] += sum(abs2, dder3_half[:, u])
            idx = U + u
            for r in 2:R
                dder3[idx] += 2 * dot(dder3_half[:, idx], dder3_half[:, u])
                dder3[u] += sum(abs2, dder3_half[:, idx])
                idx += U
            end
        end

    end

    return dder3
end


text:
using Random
using Gaugefields
using LinearAlgebra

function calc_action(gauge_action,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action,U)/NC #evaluate_Gauge_action(gauge_action,U) = tr(evaluate_Gaugeaction_untraced(gauge_action,U))
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(gauge_action,U,p,MDsteps,Dim,Uold)
    Δτ = 1/MDsteps
    gauss_distribution!(p)
    Sold = calc_action(gauge_action,U,p)
    substitute_U!(Uold,U)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,gauge_action)

        P_update!(U,p,1.0,Δτ,Dim,gauge_action)

        U_update!(U,p,0.5,Δτ,Dim,gauge_action)
    end
    Snew = calc_action(gauge_action,U,p)
    println("Sold = $Sold, Snew = $Snew")
    println("Snew - Sold = $(Snew-Sold)")
    ratio = min(1,exp(Snew-Sold))
    if rand() > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)

    end
end

function P_update!(U,p,ϵ,Δτ,Dim,gauge_action) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    dSdUμ = temps[end]
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,gauge_action,μ,U)
        mul!(temps[1],U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function HMC_test_4D(NX,NY,NZ,NT,NC,β)
    Dim = 4
    Nwing = 1

    Random.seed!(123)

    #=
    u1 = RandomGauges(NC,Nwing,NX,NY,NZ,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = RandomGauges(NC,Nwing,NX,NY,NZ,NT)
    end
    =#

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot",randomnumber="Reproducible")
    #"Reproducible"

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end

    factor = 1/(comb*U[1].NV*U[1].NC)

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("0 plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("0 polyakov loop = $(real(poly)) $(imag(poly))")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = β/2
    push!(gauge_action,β,plaqloop)

    #show(gauge_action)

    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients.
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 100
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0

    numtrj = 10
    for itrj = 1:numtrj
        accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold)
        numaccepted += ifelse(accepted,1,0)

        #plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        #println("$itrj plaq_t = $plaq_t")

        if itrj % 10 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
            println("acceptance ratio ",numaccepted/itrj)
        end
    end

    return plaq_t,numaccepted/numtrj

end

function HMC_test_2D(NX,NT,NC)
    Dim = 2
    Nwing = 1

    #=
    u1 = RandomGauges(NC,Nwing,NX,NT)
    U = Array{typeof(u1),1}(undef,Dim)
    U[1] = u1
    for μ=2:Dim
        U[μ] = RandomGauges(NC,Nwing,NX,NT)
    end
    =#

    U = Initialize_Gaugefields(NC,Nwing,NX,NT,condition = "hot",randomnumber="Reproducible")

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimension $Dim is not supported")
    end

    factor = 1/(comb*U[1].NV*U[1].NC)

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("0 plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("0 polyakov loop = $(real(poly)) $(imag(poly))")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette",Dim=Dim)
    append!(plaqloop,plaqloop')
    β = 5.7*(NC/3)/2
    push!(gauge_action,β,plaqloop)

    #show(gauge_action)

    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients.
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 100
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0

    numtrj = 10
    for itrj = 1:numtrj
        accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold)
        numaccepted += ifelse(accepted,1,0)

        #plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        #println("$itrj plaq_t = $plaq_t")

        if itrj % 10 == 0
            @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            println("$itrj plaq_t = $plaq_t")
            poly = calculate_Polyakov_loop(U,temp1,temp2)
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
            println("acceptance ratio ",numaccepted/itrj)
        end
    end
    return plaq_t,numaccepted/numtrj

end

println("2D system")
@testset "2D" begin
    NX = 4
    #NY = 4
    #NZ = 4
    NT = 4
    Nwing = 1

    @testset "NC=1" begin
        β = 2.3
        NC = 1
        println("NC = $NC")
        #val =0.6414596466929057
        val =  0.9768786716327604
        @time plaq_t,ratio = HMC_test_2D(NX,NT,NC)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =0.6414596466929057
        val =  0.9768786716327604
        @time plaq_t,ratio = HMC_test_2D(NX,NT,NC)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        #val = 0.5779454661484242
        val = 0.9656356864814539
        @time plaq_t,ratio = HMC_test_2D(NX,NT,NC)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        #val  =0.19127260002797497
        val = 0.8138836242603148
        @time plaq_t,ratio = HMC_test_2D(NX,NT,NC)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

end

println("4D system")
@testset "4D" begin
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 1

    @testset "NC=2" begin
        β = 2.3
        NC = 2
        println("NC = $NC")
        #val =0.6414596466929057
        #val = 0.5920897445000382
        val  =0.9440125563836135
        @time plaq_t,ratio = HMC_test_4D(NX,NY,NZ,NT,NC,β)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=3" begin
        β = 5.7
        NC = 3
        println("NC = $NC")
        #val = 0.5779454661484242
        #val  =0.9440125563836135
        #val = 0.5385142466966718
        val = 0.8786515255315753
        @time plaq_t,ratio = HMC_test_4D(NX,NY,NZ,NT,NC,β)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

    @testset "NC=4" begin
        β = 5.7
        NC = 4
        println("NC = $NC")
        #val  =0.19127260002797497
        #val = 0.1904815857904191
        val = 0.7301232810349298
        @time plaq_t,ratio  =HMC_test_4D(NX,NY,NZ,NT,NC,β)
        @test ratio > 0.5
        #@test abs(plaq_t-val)/abs(val) < eps
    end

end

ext: jl


text:
module Weave
using Compat

#Contains report global properties
type Report <: Display
  source::String
  documentationmode::Bool
  cwd::String
  basename::String
  formatdict
  pending_code::String
  cur_result::String
  fignum::Int
  figures::Array
  term_state::Symbol
  cur_chunk::Dict{Symbol, Any}

  function Report()
        new("", false, "", "",  Any[], "", "", 1, Any[], :text, @compat Dict{Symbol, Any}() )
  end

end

const report = Report()

const supported_mime_types =
    [MIME"image/png",
     MIME"text/plain"]

function Base.display(doc::Report, data)
    for m in supported_mime_types
        if mimewritable(m(), data)
            display(doc, m(), data)
            brea
        end
    end
end

function list_out_formats()
  for format = keys(formats)
      println(string(format,": ",  formats[format].description))
  end
end

#module ReportSandBox
#end

function tangle(source ; out_path=:doc, informat="noweb")
    cwd, fname = splitdir(abspath(source))
    basename = splitext(fname)[1]

    #Set the output directory
    if out_path == :pwd
        cwd = pwd()
    elseif out_path != :doc
        cwd = out_path
    end

    outname = "$(cwd)/$(basename).jl"
    open(outname, "w") do io
        for chunk in read_document(source, informat)
            if chunk[:type] == "code"
                write(io, chunk[:content]*"\n")
            end
        end
    end

    info("Writing to file $(basename).jl")
end

function weave(source ; doctype = "pandoc", plotlib="Gadfly", informat="noweb", out_path=:doc, fig_path = "figures", fig_ext = nothing)

    cwd, fname = splitdir(abspath(source))
    basename = splitext(fname)[1]
    formatdict = formats[doctype].formatdict
    if fig_ext == nothing
        rcParams[:chunk_defaults][:fig_ext] = formatdict[:fig_ext]
    else
        rcParams[:chunk_defaults][:fig_ext] = fig_ext
    end

    #Set the output directory
    if out_path == :doc
        report.cwd = cwd
    elseif out_path == :pwd
        report.cwd = pwd()
    else
        report.cwd = out_path
    end

    report.source = source
    report.basename = basename
    rcParams[:chunk_defaults][:fig_path] = fig_path
    report.formatdict = formatdict

    if plotlib == nothing
        rcParams[:chunk_defaults][:fig] = false
    else
        l_plotlib = lowercase(plotlib)
        rcParams[:chunk_defaults][:fig] = true
        if l_plotlib == "winston"
            eval(parse("""include(Pkg.dir("Weave","src","winston.jl"))"""))
            rcParams[:plotlib] = "Winston"
        elseif l_plotlib == "pyplot"
            eval(Expr(:using, :PyPlot))
            rcParams[:plotlib] = "PyPlot"
        elseif l_plotlib == "gadfly"
            eval(parse("""include(Pkg.dir("Weave","src","gadfly.jl"))"""))
            rcParams[:plotlib] = "Gadfly"
        end
    end

    pushdisplay(report)
    parsed = read_document(source, informat)
    executed = run(parsed)
    popdisplay(report)
    formatted = format(executed, doctype)
    outname = "$(report.cwd)/$(report.basename).$(formatdict[:extension])"
    open(outname, "w") do io
        write(io, join(formatted, "\n"))
    end

    info("Report weaved to $(report.basename).$(formatdict[:extension])")

end

function run_block(code_str)
    oldSTDOUT = STDOUT
    result = ""

    rw, wr = redirect_stdout()
    #If there is nothing to read code will hang
    println()

    try
      n = length(code_str)
      pos = 1 #The first character is extra line end
      while pos < n
          oldpos = pos
          code, pos = parse(code_str, pos)
          s = eval(ReportSandBox, code)
          if rcParams[:plotlib] == "Gadfly"
              s != nothing && display(s)
          end
      end
    finally

      redirect_stdout(oldSTDOUT)
      close(wr)
      result = readall(rw)
      close(rw)
    end

    return string("\n", result)
end

function run_term(code_str)
    prompt = "\njulia> "
    codestart = "\n\n"*report.formatdict[:codestart]

    if haskey(report.formatdict, :indent)
        prompt = indent(prompt, report.formatdict[:indent])
    end

    #Emulate terminal
    n = length(code_str)
    pos = 2 #The first character is extra line end
    while pos < n
        oldpos = pos
        code, pos = parse(code_str, pos)

        report.term_state == :fig && (report.cur_result*= codestart)
	    prompts = string(prompt, rstrip(code_str[oldpos:(pos-1)]), "\n")
	    report.cur_result *= prompts
        report.term_state = :text
        s = eval(ReportSandBox, code)
        s != nothing && display(s)
    end

    return string(report.cur_result)
end

function run(parsed)
    #Clear sandbox for each document
    #Raises a warning, couldn't find a "cleaner"
    #way to do it.
    eval(parse("module ReportSandBox\nend"))

    i = 1
    for chunk in copy(parsed)
        if chunk[:type] == "code"
            #print(chunk["content"])
            info("Weaving chunk $(chunk[:number]) from line $(chunk[:start_line])")
            defaults = copy(rcParams[:chunk_defaults])
            options = copy(chunk[:options])
            try
                options = merge(rcParams[:chunk_defaults], options)
            catch
                options = rcParams[:chunk_defaults]
                warn("Invalid format for chunk options line: $(chunk[:start_line])")
            end

            merge!(chunk, options)
            delete!(chunk, :options)

            if !chunk[:eval]
                chunk[:result] = ""
                chunk[:fig] = false
                parsed[i] = copy(chunk)
                i += 1
                continue #Do nothing if eval is false
            end

            report.fignum = 1
            report.cur_result = ""
            report.figures = String[]
            report.cur_chunk = chunk
            report.term_state = :text
            if haskey(report.formatdict, :out_width) && chunk[:out_width] == nothing
                chunk[:out_width] = report.formatdict[:out_width]
            end

            if chunk[:term]
                chunk[:result] = run_term(chunk[:content])
                chunk[:term_state] = report.term_state
            else
                chunk[:result] = run_block(chunk[:content])
            end

            if rcParams[:plotlib] == "PyPlot"
                chunk[:fig] && (chunk[:figure] = savefigs(chunk))
            else
                chunk[:fig] && (chunk[:figure] = copy(report.figures))
            end
        end
        parsed[i] = copy(chunk)
        i += 1
    end
    parsed
end

function savefigs(chunk)
    l_plotlib = lowercase(rcParams[:plotlib])
    if l_plotlib == "pyplot"
      return savefigs_pyplot(chunk)
    end
end

function savefigs_pyplot(chunk)
    fignames = String[]
    ext = report.formatdict[:fig_ext]
    figpath = joinpath(report.cwd, chunk[:fig_path])
    isdir(figpath) || mkdir(figpath)
    chunkid = (chunk[:name] == nothing) ? chunk[:number] : chunk[:name]
    #Iterate over all open figures, save them and store names
    for fig = plt.get_fignums()
        full_name, rel_name = get_figname(report, chunk, fignum=fig)
        savefig(full_name, dpi=chunk[:dpi])
        push!(fignames, rel_name)
        plt.draw()
        plt.close()
    end
    return fignames
end

function Base.display(report::Report, m::MIME"text/plain", data)
  if report.term_state == :fig #Catch Winston plot command output
    report.cur_result *= "\n" * report.formatdict[:codestart] * "\n"
  end

  s = reprmime(m, data)
  haskey(report.formatdict, :indent) && (s = indent(s, report.formatdict[:indent]))

  report.cur_result *= s * "\n"

  if report.term_state == :fig #Catch Winston plot command output
    report.cur_result *= "\n" * report.formatdict[:codeend] * "\n"
  end
end

function get_figname(report::Report, chunk; fignum = nothing)
    figpath = joinpath(report.cwd, chunk[:fig_path])
    isdir(figpath) || mkdir(figpath)
    ext = chunk[:fig_ext]
    fignum == nothing && (fignum = report.fignum)

    chunkid = (chunk[:name] == nothing) ? chunk[:number] : chunk[:name]
    full_name = joinpath(report.cwd, chunk[:fig_path], "$(report.basename)_$(chunkid)_$(fignum)$ext")
    rel_name = "$(chunk[:fig_path])/$(report.basename)_$(chunkid)_$(fignum)$ext" #Relative path is used in output
    return full_name, rel_name
end

export weave, list_out_formats, tangle

include("config.jl")
include("readers.jl")
include("formatters.jl")
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for writing the curtailment values of the different variable renewable resources.
"""
function write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	dfGen = inputs["dfGen"]
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	dfCurtailment = DataFrame(Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone], AnnualSum = Array{Union{Missing,Float32}}(undef, G))
	for i in 1:G
		if i in inputs["VRE"]
			dfCurtailment[!,:AnnualSum][i] = sum(inputs["omega"].*(inputs["pP_Max"][i,:]).*value.(EP[:eTotalCap])[i,:].- inputs["omega"].*value.(EP[:vP])[i,:])
		else
			dfCurtailment[!,:AnnualSum][i] = 0
		end
	end
	if setup["ParameterScale"] ==1
		dfCurtailment.AnnualSum = dfCurtailment.AnnualSum * ModelScalingFactor
		dfCurtailment = hcat(dfCurtailment, DataFrame(( ModelScalingFactor * (inputs["pP_Max"]).*value.(EP[:eTotalCap]).- value.(EP[:vP])), :auto))
	else
		dfCurtailment = hcat(dfCurtailment, DataFrame(((inputs["pP_Max"]).*value.(EP[:eTotalCap]).- value.(EP[:vP])), :auto))
	end

	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
	rename!(dfCurtailment,auxNew_Names)
	total = DataFrame(["Total" 0 sum(dfCurtailment[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
	for t in 1:T
		total[:,t+3] .= sum(dfCurtailment[:,Symbol("t$t")][1:G])
	end
	rename!(total,auxNew_Names)
	dfCurtailment = vcat(dfCurtailment, total)
	CSV.write(joinpath(path, "curtail.csv"), dftranspose(dfCurtailment, false), writeheader=false)
	return dfCurtailment
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
#=
Please see
https://computationalmindset.com/en/neural-networks/ordinary-differential-equation-solvers.html#ode2
for details
=#

#=Original problem
x′′ + x′ + 2x = 0
x(0)=1
x′(0)=0
=#

#=Esplicit form of equation
x′′ = -x′ - 2x
=#

#=Analytical solution
x(t) = e^(-t/2) (cos(√7 t / 2) + sin(√7 t / 2) / √7)
=#

using DifferentialEquations
using Plots

function ode_fn(dx,x,p,t)
    -dx -2.0 * x
end

an_sol(t) = exp(-t/2.0) *
            (cos(sqrt(7.0) * t / 2.0) + sin(sqrt(7.0) * t / 2.0)/sqrt(7.0))

t_begin=0.0
t_end=12.0
tspan = (t_begin,t_end)
x_init=1.0
dxdt_init=0.0

prob = SecondOrderODEProblem(ode_fn, dxdt_init, x_init, tspan)
num_sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
x_num_sol = [u[2] for u in num_sol.u]

plot(num_sol.t, an_sol.(num_sol.t),
    linewidth=2, ls=:dash,
    title="ODE 2nd order IVP solved by D.E. package",
    xaxis="t", yaxis="x",
    label="analytical",
    legend=true)
plot!(num_sol.t, x_num_sol,
    linewidth=1,
    label="numerical")

max_stars_repo_name: ettoremessina/differential-equations

max_issues_repo_name: ettoremessina/differential-equations

max_forks_repo_name: ettoremessina/differential-equations


text:
type CNS # central nervous system
    # sympathetic, parasympathetic activations
    np::Vector{Float64}
    ns::Vector{Float64}
    # heart rate and maximum ventriclular elastance
    H::Vector{Float64}
    Emaxlv::Vector{Float64}
    Emaxrv::Vector{Float64}
    # lower/upper body compartmental resistances, compliances, unstressed vols.
    R2L::Vector{Float64}
    R3L::Vector{Float64}
    R2U::Vector{Float64}
    R3U::Vector{Float64}
    C4L::Vector{Float64}
    C5L::Vector{Float64}
    C4U::Vector{Float64}
    C5U::Vector{Float64}
    V4L::Vector{Float64}
    V5L::Vector{Float64}
    V4U::Vector{Float64}
    V5U::Vector{Float64}
    # feedback model weights (alpha:ns, beta:np, gamma:bias)
    alphaH::Float64
    betaH::Float64
    gammaH::Float64
    alphalv::Float64
    gammalv::Float64
    alpharv::Float64
    gammarv::Float64
    alphaR2L::Float64
    gammaR2L::Float64
    alphaR3L::Float64
    gammaR3L::Float64
    alphaR2U::Float64
    gammaR2U::Float64
    alphaR3U::Float64
    gammaR3U::Float64
    alphaC4L::Float64
    gammaC4L::Float64
    alphaC5L::Float64
    gammaC5L::Float64
    alphaC4U::Float64
    gammaC4U::Float64
    alphaC5U::Float64
    gammaC5U::Float64
    alphaV4L::Float64
    gammaV4L::Float64
    alphaV5L::Float64
    gammaV5L::Float64
    alphaV4U::Float64
    gammaV4U::Float64
    alphaV5U::Float64
    gammaV5U::Float64
    # time constants
    tauH::Float64
    tauEmax::Float64
    tauC::Float64
    tauR::Float64
    tauV::Float64
    # activation function parameters
    Ptarget::Float64
    Paverage::Vector{Float64}
    baselinefirerate::Float64
    a::Float64
    b::Float64
    c::Float64
    # indices for baroreflex pressure averaging
    avgindexstart::Int32
    avgindexend::Int32

    function CNS()
        this = new()
        this.ns = Vector{Float64}[];
        this.np = Vector{Float64}[];
        this.H = Vector{Float64}[];
        this.Emaxlv = Vector{Float64}[];
        this.Emaxrv = Vector{Float64}[];
        this.R2L = Vector{Float64}[];
        this.R3L = Vector{Float64}[];
        this.R2U = Vector{Float64}[];
        this.R3U = Vector{Float64}[];
        this.C4L = Vector{Float64}[];
        this.C5L = Vector{Float64}[];
        this.C4U = Vector{Float64}[];
        this.C5U = Vector{Float64}[];
        this.V4L = Vector{Float64}[];
        this.V5L = Vector{Float64}[];
        this.V4U = Vector{Float64}[];
        this.V5U = Vector{Float64}[];
        this.tauH = 2;
        this.tauEmax = 2;
        this.tauC = 20;
        this.tauR = 6;
        this.tauV = 20;
        this.Ptarget = 94.6*mmHgToPa;
        this.Paverage = Vector{Float64}[];
        this.avgindexstart = 1;
        this.avgindexend = 0;
        this.baselinefirerate = 0.25;
        this.a = 1;
        this.b = -log(this.baselinefirerate);
        this.c = 6;
        return this
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module ActivationFunctions

export sigmoid, relu, softmax

"""
    sigmoid(x)

Computes the sigmoid activation of scalar `x`, returning a scalar.
"""
sigmoid(x::Number)::Number = (1.0 + exp(-x))^(-1)

"""
    sigmoid(x)

Computes the sigmoid activation of vector `x`, returning a vector.
"""
sigmoid(x::AbstractVector{<:Number})::AbstractVector{<:Number} = sigmoid.(x)

"""
    relu(x)

Computes the rectified linear unit activation of scalar `x`, returning a scalar.
"""
relu(x::Number)::Number = max(0.0, x)

"""
    relu(x)

Computes the rectified linear unit activation of vector `x`, returning a vector.
"""
relu(x::AbstractVector{<:Number})::AbstractVector{<:Number} = max.(0.0, x)

"""
    softmax(x)

Computes the softmax activation of vector `x`, returning a vector.

Unlike some other activation functions, softmax does not treat vector elements
independently.
"""
softmax(x::AbstractVector{<:Number})::AbstractVector{<:Number} =
    exp.(x .- maximum(x)) .* sum(exp.(x .- maximum(x)))^(-1)

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
module  TestManipulate

using PopGenCore
using DataFrames
using StatsBase
using Test

cats = @nancycats;
sharks = @gulfsharks;

@testset "Manipulate.jl" begin
    @testset "locations" begin
        x = rand(length(samplenames(cats))) ; y = rand(length(samplenames(cats)))
        locationdata!(cats, longitude = x, latitude = y)
        @test cats.sampleinfo.longitude == Float32.(x)
        @test cats.sampleinfo.latitude == Float32.(y)
        @test locationdata(cats).longitude == Float32.(x)
        @test locationdata(cats).latitude == Float32.(y)
    end

    @testset "decimal-minutes locations" begin
        x = fill("11 22.33W", length(samplenames(cats))) ; y = fill("-41 31.52", length(samplenames(cats)))
        locationdata!(cats, longitude = x, latitude = y)
        @test all(cats.sampleinfo.longitude .== Float32(-11.3722))
        @test all(cats.sampleinfo.latitude .== Float32(-41.5253))
    end

    @testset "loci and genotypes" begin
        @test length(loci(cats)) == 9
        @test genotypes(cats, "N115" => "fca8") == (135,135)
        N115 = genotypes(cats, "N115")
        @test length(N115) == 9
        @test N115 isa GenoArray
        @test size(genotypes(cats, ["N115", "N7"] => "fca8")) == (2,4)
        @test size(genotypes(cats, "N115" => ["fca8", "fca37"])) == (2,4)
        @test size(genotypes(cats, ["N115", "N7"] => ["fca8", "fca37"])) == (4,4)
        @test length(genotypes(sharks, "contig_475")) == 212
    end

    @testset "populations" begin
        @test length(populations(cats)) == 17
        @test populations(cats, counts = true).population isa Vector{String}
        @test populations(cats, counts = true).count isa Vector{Int}

        rn_dict = Dict("1" => "one", "2" => "two")
        populations!(cats, rn_dict)
        @test "one" ∈ cats.sampleinfo.population && "1" ∉ cats.sampleinfo.population
        @test "two" ∈ cats.sampleinfo.population && "2" ∉ cats.sampleinfo.population

        rn_vect = string.(1:17)
        populations!(cats, rn_vect)
        @test "one" ∉ cats.sampleinfo.population && "1" ∈ cats.sampleinfo.population
        @test "two" ∉ cats.sampleinfo.population && "2" ∈ cats.sampleinfo.population
        rn_vect = string.(1:237)
        populations!(cats, rn_vect)
        @test cats.metadata.populations == 237
        @test cats.sampleinfo.population == rn_vect
        @test_throws DimensionMismatch populations!(cats, string.(1:22))

        rn_vect_ii = ["N215", "N297"]
        rn_vect_iin = ["one", "seventeen"]
        populations!(cats, rn_vect_ii, rn_vect_iin)
        @test "one" ∈ cats.sampleinfo.population && "seventeen" ∈ cats.sampleinfo.population
        @test cats.sampleinfo[cats.sampleinfo.name .== "N215", :population] == ["one"]
        @test cats.sampleinfo[cats.sampleinfo.name .== "N297", :population] == ["seventeen"]
    end

    @testset "exclusion" begin
        cats = @nancycats ;
        tmp = exclude(cats, name = "N100", population = ["1", "15"])
        @test length(samplenames(tmp)) == 215
        @test length(populations(tmp)) == 15

        tmp = exclude(cats, name = "N7", locus= "fca8", population = "3")
        @test length(loci(tmp)) == 8
        @test length(populations(tmp)) == 16
        @test length(samplenames(tmp)) == 224

        tmp = keep(cats, name = "N104")
        @test length(samplenames(tmp)) == 1
        tmp = keep(cats, name = ["N104", "N105"])
        @test length(samplenames(tmp)) == 2
        tmp = keep(cats, locus = "fca8")
        @test length(loci(tmp)) == 1
        tmp = keep(cats, locus = ["fca8","fca37"])
        @test length(loci(tmp)) == 2
        tmp = keep(cats, population = "1")
        @test length(populations(tmp)) == 1
        tmp = keep(cats, population = 1:2)
        @test length(populations(tmp)) == 2
    end
end

end # module


text:
const COLOR_SCHEMES = [
    (0xffffff, 0x000000, 0xff0000, 0x00ff00, 0x0000ff, 0x00ffff, 0xffff00, 0xff00ff),
    (0x282c34, 0xd7dae0, 0xcb4e42, 0x99c27c, 0x85a9fc, 0x5ab6c1, 0xd09a6a, 0xc57bdb),
    (0xfdf6e3, 0x657b83, 0xdc322f, 0x859900, 0x268bd2, 0x2aa198, 0xb58900, 0xd33682),
    (0x002b36, 0x839496, 0xdc322f, 0x859900, 0x268bd2, 0x2aa198, 0xb58900, 0xd33682)
]

const BASIC_COLORS = [ 0, 1, 984, 987, 989, 983, 994, 988 ]

const SERIES_COLORS = [989, 982, 980, 981, 996, 983, 995, 988, 986, 990, 991, 984, 992, 993, 994, 987, 985, 997, 998, 999]

const COLOR_INDICES = Dict{Symbol, Int}(
    :scheme => 0,
    :colormap => GR.COLORMAP_VIRIDIS
)

function resetcolors()
    for k in keys(COLOR_INDICES)
        if k ∉ (:scheme, :colormap)
            COLOR_INDICES[k] = 0
        end
    end
    return nothing
end

# Color codes

"""
    rgb(color)

Return the normalized RGB values (between 0 and 1)
corresponding to a given hexadecimal color value.

# Examples

```jldoctest
julia> GRUtils.rgb(0xff3300)
(1.0, 0.2, 0.0)
```
"""
function rgb(color::Integer)
    r = float((color >> 16) & 0xff) / 255.0
    g = float((color >> 8)  & 0xff) / 255.0
    b = float( color        & 0xff) / 255.0
    return (r, g, b)
end

"""
    color(r, g, b)

Return the hexadecimal integer corresponding to the given RGB values.

# Examples

```jldoctest
julia> GRUtils.color(1, 0.2, 0)
0x00ff3300
```
"""
function color(r, g, b)
    rint = round(UInt32, r * 255)
    gint = round(UInt32, g * 255)
    bint = round(UInt32, b * 255)
    return rint << 16 + gint << 8 + bint
end

"""
    switchbytes(hexcolor::UInt32)

Switch the R and B channels of hexadecimal colors codes.

Use this function to turn a byte-ordered hexadecimal color code (ARGB, little-endian)
into a word-ordered color code (RGBA, big endian) and vice versa.

# Examples

```jldoctest
julia> GRUtils.switchbytes(0xfe12ac34)
0xfe34ac12
```
"""
switchbytes(c::UInt32) = (c & 0xff00ff00) + (c & 0x00ff0000) >> 16 + (c & 0x000000ff) << 16

"""
    colorindex(hexcolor[, byteorder=false])

Define a color given by an hexadecimal code and return its index in GR's color map.

This is a wrapper to `GR.inqcolorfromrgb` for hexadecimal color codes.
By default it considers that RGB values are ordered from the most to the least
significant bytes. If it is a byte-ordered code, set the second argument to `true`.
"""
function colorindex(hexcolor, byteorder=false)::Int
    if byteorder
        hexcolor = switchbytes(hexcolor)
    end
    GR.inqcolorfromrgb(rgb(hexcolor)...)
end

# Colormaps

const COLORMAPS = Dict( k => i-1 for (i, k) in enumerate((
    "uniform", "temperature", "grayscale", "glowing", "rainbowlike",
    "geologic", "greenscale", "cyanscale", "bluescale", "magentascale",
    "redscale", "flame", "brownscale", "pilatus", "autumn", "bone",
    "cool", "copper", "gray", "hot", "hsv", "jet", "pink", "spectral",
    "spring", "summer", "winter", "gistearth", "gistheat", "gistncar",
    "gistrainbow", "giststern", "afmhot", "brg", "bwr", "coolwarm",
    "cmrmap", "cubehelix", "gnuplot", "gnuplot2", "ocean", "rainbow",
    "seismic", "terrain", "viridis", "inferno", "plasma", "magma"
)))

"""
    colormap(cmap)
    colormap([cmap,] T::Type{<:Integer}[, byteorder::Bool=false])

Set or retrieve the values of the current colormap.

Give the identifier `cmap` to set the colormap. This identifier can be
the number or a string with the name of any of the
[GR built-in colormaps](https://gr-framework.org/colormaps.html).

Give an integer data type `T` to return the set of colors defined
in the current colormap, as a vector of hexadecimal color codes.
By default those codes represent the word-ordered RGB values
in a little-endian system. The optional argument `byteorder`
can be set to `true` in order to retrieve the codes as byte-ordered codes.

At least one of the arguments `cmap::Union{Integer, AbstractString}`
or `T::Type{<:Integer}` are required. If `cmap` not given, the previous
colormap is left as current. If `T` is not given, the function returns
`nothing`.

# Examples

```julia
$(_example("colormap"))
```
"""
function colormap(cmap::Integer)
    GR.setcolormap(cmap)
    COLOR_INDICES[:colormap] = cmap
    return nothing
end

function colormap(cmap::AbstractString)
    cmap = lowercase(replace(cmap, (' ', '_') => ""))
    colormap(lookup(cmap, COLORMAPS))
end

function colormap(::Type{T}, byteorder::Bool=false) where {T <: Integer}
    if byteorder
        return [T(GR.inqcolor(i)) for i ∈ 1000:1255]
    else
        return [GR.inqcolor(i) |> UInt32 |> switchbytes |> T for i ∈ 1000:1255]
    end
end

function colormap(cmap::Union{Integer, AbstractString}, T::DataType, args...)
    colormap(cmap)
    colormap(T, args...)
end

"""
    to_rgba(value[, alpha, cmap])

Calculate the hexadecimal color code of a normalized value between 0 and 1
in the current colormap (byte-ordered, big-endian RGBA).

Optionally, the level of the `alpha` channel can be passed as a value
between 0 (transparent) and 1 (opaque), and a custom colormap can be
defined by `cmap` with the corresponding RGB values
(as byte-ordered hexadecimal colors or triplets of normalized float values).
"""
function to_rgba(value, alpha, cmap=colormap(UInt32, true)::Vector{<:Integer})
    isnan(value) && return zero(UInt32)
    color = UInt32(cmap[round(Int, value * 255 + 1)])
    return color + round(UInt32, alpha * 255) << 24
end

function to_rgba(value)
    isnan(value) && return zero(UInt32)
    colormap(UInt32, true)[round(Int, value * 255 + 1)] + 0xff000000
end

# Color schemes

"""
    applycolorscheme(scheme)

Apply the given color scheme, coded as an integer number between 0 and 4.
See [`colorscheme`](@ref) for the values of the color schemes.
"""
function applycolorscheme(scheme)
    scheme == 0 && return nothing
    # Replace the basic color indices
    for colorind in 1:8
        color = COLOR_SCHEMES[scheme][colorind]
        r, g, b = rgb(color)
        # Set 1-8
        GR.setcolorrep(colorind - 1, r, g, b)
        # In set 980-999 (except for the first scheme)
        if scheme ≠ 1
            GR.setcolorrep(BASIC_COLORS[colorind], r, g, b)
        end
    end
    # Background RGB values
    r, g, b = rgb(COLOR_SCHEMES[scheme][1])
    # Difference between foreground and background
    rdiff, gdiff, bdiff = rgb(COLOR_SCHEMES[scheme][2]) .- (r, g, b)
    # replace the 12 "grey" shades
    for (colorind, f) in enumerate(LinRange(1, 0, 12))
        GR.setcolorrep(79 + colorind, r + f*rdiff, g + f*gdiff, b + f*bdiff)
    end
    return nothing
end

"""
    colorscheme(scheme)

Set the color scheme for subsequent plots.

The value of the scheme can be one of the following numbers or strings:

* 0: `"none"`
* 1: `"light"`
* 2: `"dark"`
* 3: `"solarized light"`
* 4: `"solarized dark"`

The scheme applies to all plots that are drawn after calling
`colorscheme`, even if they have been created beforehand.
To apply a particular scheme to some plot or figure, use
[`colorscheme!`](@ref).

If the scheme is `"none"` (`0`), the standard scheme is set, which
is the same as `"light"` (1), except that the default background of the figure
is transparent.

# Examples

```julia
$(_example("colorscheme"))
```
"""
function colorscheme(scheme::Int)
    COLOR_INDICES[:scheme] = scheme
    return nothing
end

function colorscheme(scheme::AbstractString)
    scheme = replace(scheme, " " => "")
    scheme = lowercase(scheme)
    scheme_dict = Dict("none" => 0, "light" => 1, "dark" => 2,
        "solarizedlight" => 3, "solarizeddark" => 4)
    colorscheme(scheme_dict[scheme])
end


text:
#=

   Dice problem.

   From Think Bayes, page 21.
   """
   Suppose I have a box of dice that contains a 4-sided die, a 6-sided die, an
   8-sided die, a 12-sided die, and a 20-sided die. If you have ever played
   Dungeons & Dragons, you know what I am talking about.

   Suppose I select a die from the box at random, roll it, and get a 6. What is
   the probability that I rolled each die?

   ...

   What if we roll a few more times and get 6, 8, 7, 7, 5, and 4?
   """

   As we get more information (i.e. that all values are <= 8) we get more and
   more sure that it's the 8-sided die.

   throws: [6]
   Distributions of variable dice (num:0)
   6.00000 =>    3881  (0.388100)
   8.00000 =>    2965  (0.296500)
   12.00000 =>    1992  (0.199200)
   20.00000 =>    1162  (0.116200)

   throws: [6, 6]
   Distributions of variable dice (num:0)
   6.00000 =>    5217  (0.521700)
   8.00000 =>    3008  (0.300800)
   12.00000 =>    1321  (0.132100)
   20.00000 =>     454  (0.045400)

   throws: [6, 6, 8]
   Distributions of variable dice (num:0)
   8.00000 =>    7438  (0.743800)
   12.00000 =>    2102  (0.210200)
   20.00000 =>     460  (0.046000)

   throws: [6, 6, 8, 7]
   Distributions of variable dice (num:0)
   8.00000 =>    8183  (0.818300)
   12.00000 =>    1615  (0.161500)
   20.00000 =>     202  (0.020200)

   throws: [6, 6, 8, 7, 7]
   Distributions of variable dice (num:0)
   8.00000 =>    8766  (0.876600)
   12.00000 =>    1130  (0.113000)
   20.00000 =>     104  (0.010400)

   throws: [6, 6, 8, 7, 7, 5]
   Distributions of variable dice (num:0)
   8.00000 =>    9163  (0.916300)
   12.00000 =>     792  (0.079200)
   20.00000 =>      45  (0.004500)

   throws: [6, 6, 8, 7, 7, 5, 4]
   Distributions of variable dice (num:0)
   8.00000 =>    9413  (0.941300)
   12.00000 =>     575  (0.057500)
   20.00000 =>      12  (0.001200)

   However, if the next throw is 11 then all previous bets are off:

   throws: [6, 6, 8, 7, 7, 5, 4, 11]
   Distributions of variable dice (num:0)
   12.00000 =>    9851  (0.985100)
   20.00000 =>     149  (0.014900)

   Cf ~/blog/dice_problem.blog
      ~/psi/dice_problem.psi
      ~/webppl/dice_problem.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function dice_problem()
    throws_len = length(throws)

    values=[4,6,8,12,20]
    values_len = length(values)

    # Index of the values array
    ix ~ Categorical([1/values_len for _ in 1:values_len])
    # The value of the die
    dice ~ Dirac(values[ix])

    throws ~ filldist(DiscreteUniform(1,dice),throws_len)

end

throws=[6,6,8,7,7,5,4,11]

function run_dice_problem(throws)
    println("throws: ", throws)
    model = dice_problem(throws)

    num_chains = 4

    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    chns = sample(model, SMC(10), 10_000)
    # chns = sample(model, IS(), 100_000)

    # display(chns)
    # display(plot(chns))

    show_var_dist_pct(chns,:dice)

end

for i in 1:length(throws)
    run_dice_problem(throws[1:i])
    println()
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Stack models together

export Model2Fit_Stacked

@doc raw"""

Define a composite model, where each submodel is responsible of a given ROI

| (1)        | (2)       | ...
|------------|-----------| ...
| models[1]  | model[2]  | ...

Domain (1) range is defined as 1:ROI_sizes[1]
Domain (2) range is defined as (ROI_sizes[1]+1):ROI_sizes[2]
etc...
"""
struct Model2Fit_Stacked{MODEL2FIT<:Abstract_Model2Fit} <: Abstract_Model2Fit
    _models::Vector{MODEL2FIT}
    _X_ranges::Vector{UnitRange{Int}}
    _θ_ranges::Vector{UnitRange{Int}}

    _ROI_lengths::Vector{Int}

    function Model2Fit_Stacked(models::ABSVECT_MODEL2FIT,
                               ROI_lengths::AbstractVector{Int}) where {MODEL2FIT<:Abstract_Model2Fit,
                                                                        ABSVECT_MODEL2FIT<:AbstractVector{MODEL2FIT}}

        @assert length(models) == length(ROI_lengths)

        # Compute θ and X ranges associated to each submodel
        #
        n_models = length(models)

        X_ranges = Vector{UnitRange{Int}}(undef,n_models)
        θ_ranges = Vector{UnitRange{Int}}(undef,n_models)

        if n_models>0

            X_range_first = 1
            θ_range_first = 1

            resize!(X_ranges,0)
            resize!(θ_ranges,0)

            for (model,roi_length) in zip(models,ROI_lengths)

                X_range = X_range_first:(X_range_first + roi_length -1)
                push!(X_ranges,X_range)

                θ_length = parameter_size(model)
                θ_range = θ_range_first:(θ_range_first + θ_length - 1)
                push!(θ_ranges,θ_range)

                X_range_first += roi_length
                θ_range_first += θ_length
            end

        end

        new{MODEL2FIT}(models,X_ranges,θ_ranges)

    end
end

# Helper  ================
#
_get_XY(model::Model2Fit_Stacked,submodel_idx::Int,XY::AbstractVector) = @view XY[model._X_ranges[submodel_idx]]
_get_θ(model::Model2Fit_Stacked,submodel_idx::Int,θ::AbstractVector) = @view θ[model._θ_ranges[submodel_idx]]

# Visit  ================
#
visit_submodel_size(model::Model2Fit_Stacked) = length(model._models)

visit_get_submodel(model::Model2Fit_Stacked,submodel_idx::Int) = model._models[submodel_idx]
visit_get_Y(model::Model2Fit_Stacked,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = _get_XY(model,submodel_idx,Y)
visit_get_X(model::Model2Fit_Stacked,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = _get_XY(model,submodel_idx,X)
visit_get_θ(model::Model2Fit_Stacked,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = _get_θ(model,submodel_idx,θ)

# Interface ================
#
parameter_size(m::Model2Fit_Stacked) = last(last(m._θ_ranges))

function accumulate_y!(model::Model2Fit_Stacked,Y::AbstractVector,X::AbstractVector,θ::AbstractVector)
    n_submodel = visit_submodel_size(model)

    for i in 1:n_submodel
        submodel = visit_get_submodel(model,i)
        subY = _get_XY(model,i,Y)
        subX = _get_XY(model,i,X)
        subθ = _get_θ(model,i,θ)

        accumulate_y!(submodel,subY,subX,subθ)
    end

    Y
end


text:
# # ImageInspector

# ImageInspector is a small package for educational purposes. Its main goal is not presenting functionality, but presenting package structure. This is its short documentation created in the package [Literate](https://fredrikekre.github.io/Literate.jl/v2) which uses the [Markdown](https://www.markdownguide.org/cheat-sheet) syntax.

# To use the package, we need to load first the required packages.

using ImageInspector
using Plots

# ## Grayscale images

# As a test example, we create the real matrix `img1` representing a circle. We first discretize the domain $[-1,1]$ in `xs`. We assign black colour whenever $x^2 + y^2 \le 1$. Since the white colour is represented by `[1; 1; 1]` and the black colour by `[0; 0; 0]`, we can do it by the following code:

xs = -1:0.001:1
img1 = [x^2+y^2>1 for x in xs, y in xs];

# This is a two-dimensional matrix, which represents a grayscale image. We convert it to an image by calling `image` and then we plot it.

plot(image(img1); axis = nothing, border = :none)

# We can convert it to a three-dimensinal representation by the `reshape` function. The resulting plot is the same even though a different method for the `image` function is called.

img2 = reshape(img1, (size(img1)..., 1))

plot(image(img2); axis = nothing, border = :none)

# ## Coloured images

# We can repeat a similar procedure for coloured images. We create again a three-dimensional matrix, convert it to an image and then plot it.

img2 = ones(size(img1)..., 3)
img2[:,:,2:3] .= [x^2+y^2>1 for x in xs, y in xs]

plot(image(img2); axis = nothing, border = :none)


text:
# # Linear regression

# We estimate simple linear regression model with a half-T prior.
# First, we load the packages we use.

using TransformVariables, LogDensityProblems, DynamicHMC, MCMCDiagnostics,
    Parameters, Statistics, Distributions

# Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

"""
Linear regression model ``y ∼ Xβ + ϵ``, where ``ϵ ∼ N(0, σ²)`` IID.

Flat prior for `β`, half-T for `σ`.
"""
struct LinearRegressionProblem{TY <: AbstractVector, TX <: AbstractMatrix,
                               Tν <: Real}
    "Observations."
    y::TY
    "Covariates"
    X::TX
    "Degrees of freedom for prior."
    ν::Tν
end

# Then make the type callable with the parameters *as a single argument*.

function (problem::LinearRegressionProblem)(θ)
    @unpack y, X, ν = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    loglikelihood(Normal(0, σ), y .- X*β) + logpdf(TDist(ν), σ)
end

# We should test this, also, this would be a good place to benchmark and
# optimize more complicated problems.

N = 100
X = hcat(ones(N), randn(N, 2));
β = [1.0, 2.0, -1.0]
σ = 0.5
y = X*β .+ randn(N) .* σ;
p = LinearRegressionProblem(y, X, 1.0);
p((β = β, σ = σ))

# For this problem, we write a function to return the transformation (as it varies with the number of covariates).

problem_transformation(p::LinearRegressionProblem) =
    as((β = as(Array, size(p.X, 2)), σ = asℝ₊))

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(problem_transformation(p), p)
∇P = FluxGradientLogDensity(P);

# Finally, we sample from the posterior. `chain` holds the chain (positions and
# diagnostic information), while the second returned value is the tuned sampler
# which would allow continuation of sampling.

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

# We use the transformation to obtain the posterior from the chain.

posterior = transform.(Ref(∇P.transformation), get_position.(chain));

# Extract the parameter posterior means: `β`,

posterior_β = mean(first, posterior)

# then `σ`:

posterior_σ = mean(last, posterior)

# Effective sample sizes (of untransformed draws)

ess = mapslices(effective_sample_size,
                get_position_matrix(chain); dims = 1)

# NUTS-specific statistics

NUTS_statistics(chain)


text:
#
# Autogenerated by Thrift Compiler (1.0.0-dev)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

# types encapsulating arguments and return values of method OpenSession

type OpenSession_args
  req::TOpenSessionReq
  OpenSession_args() = (o=new(); fillunset(o); o)
end # type OpenSession_args

type OpenSession_result
  success::TOpenSessionResp
  OpenSession_result() = (o=new(); fillunset(o); o)
  OpenSession_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type OpenSession_result
meta(t::Type{OpenSession_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method CloseSession

type CloseSession_args
  req::TCloseSessionReq
  CloseSession_args() = (o=new(); fillunset(o); o)
end # type CloseSession_args

type CloseSession_result
  success::TCloseSessionResp
  CloseSession_result() = (o=new(); fillunset(o); o)
  CloseSession_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type CloseSession_result
meta(t::Type{CloseSession_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetInfo

type GetInfo_args
  req::TGetInfoReq
  GetInfo_args() = (o=new(); fillunset(o); o)
end # type GetInfo_args

type GetInfo_result
  success::TGetInfoResp
  GetInfo_result() = (o=new(); fillunset(o); o)
  GetInfo_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetInfo_result
meta(t::Type{GetInfo_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method ExecuteStatement

type ExecuteStatement_args
  req::TExecuteStatementReq
  ExecuteStatement_args() = (o=new(); fillunset(o); o)
end # type ExecuteStatement_args

type ExecuteStatement_result
  success::TExecuteStatementResp
  ExecuteStatement_result() = (o=new(); fillunset(o); o)
  ExecuteStatement_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type ExecuteStatement_result
meta(t::Type{ExecuteStatement_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetTypeInfo

type GetTypeInfo_args
  req::TGetTypeInfoReq
  GetTypeInfo_args() = (o=new(); fillunset(o); o)
end # type GetTypeInfo_args

type GetTypeInfo_result
  success::TGetTypeInfoResp
  GetTypeInfo_result() = (o=new(); fillunset(o); o)
  GetTypeInfo_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetTypeInfo_result
meta(t::Type{GetTypeInfo_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetCatalogs

type GetCatalogs_args
  req::TGetCatalogsReq
  GetCatalogs_args() = (o=new(); fillunset(o); o)
end # type GetCatalogs_args

type GetCatalogs_result
  success::TGetCatalogsResp
  GetCatalogs_result() = (o=new(); fillunset(o); o)
  GetCatalogs_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetCatalogs_result
meta(t::Type{GetCatalogs_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetSchemas

type GetSchemas_args
  req::TGetSchemasReq
  GetSchemas_args() = (o=new(); fillunset(o); o)
end # type GetSchemas_args

type GetSchemas_result
  success::TGetSchemasResp
  GetSchemas_result() = (o=new(); fillunset(o); o)
  GetSchemas_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetSchemas_result
meta(t::Type{GetSchemas_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetTables

type GetTables_args
  req::TGetTablesReq
  GetTables_args() = (o=new(); fillunset(o); o)
end # type GetTables_args

type GetTables_result
  success::TGetTablesResp
  GetTables_result() = (o=new(); fillunset(o); o)
  GetTables_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetTables_result
meta(t::Type{GetTables_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetTableTypes

type GetTableTypes_args
  req::TGetTableTypesReq
  GetTableTypes_args() = (o=new(); fillunset(o); o)
end # type GetTableTypes_args

type GetTableTypes_result
  success::TGetTableTypesResp
  GetTableTypes_result() = (o=new(); fillunset(o); o)
  GetTableTypes_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetTableTypes_result
meta(t::Type{GetTableTypes_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetColumns

type GetColumns_args
  req::TGetColumnsReq
  GetColumns_args() = (o=new(); fillunset(o); o)
end # type GetColumns_args

type GetColumns_result
  success::TGetColumnsResp
  GetColumns_result() = (o=new(); fillunset(o); o)
  GetColumns_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetColumns_result
meta(t::Type{GetColumns_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetFunctions

type GetFunctions_args
  req::TGetFunctionsReq
  GetFunctions_args() = (o=new(); fillunset(o); o)
end # type GetFunctions_args

type GetFunctions_result
  success::TGetFunctionsResp
  GetFunctions_result() = (o=new(); fillunset(o); o)
  GetFunctions_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetFunctions_result
meta(t::Type{GetFunctions_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetOperationStatus

type GetOperationStatus_args
  req::TGetOperationStatusReq
  GetOperationStatus_args() = (o=new(); fillunset(o); o)
end # type GetOperationStatus_args

type GetOperationStatus_result
  success::TGetOperationStatusResp
  GetOperationStatus_result() = (o=new(); fillunset(o); o)
  GetOperationStatus_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetOperationStatus_result
meta(t::Type{GetOperationStatus_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method CancelOperation

type CancelOperation_args
  req::TCancelOperationReq
  CancelOperation_args() = (o=new(); fillunset(o); o)
end # type CancelOperation_args

type CancelOperation_result
  success::TCancelOperationResp
  CancelOperation_result() = (o=new(); fillunset(o); o)
  CancelOperation_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type CancelOperation_result
meta(t::Type{CancelOperation_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method CloseOperation

type CloseOperation_args
  req::TCloseOperationReq
  CloseOperation_args() = (o=new(); fillunset(o); o)
end # type CloseOperation_args

type CloseOperation_result
  success::TCloseOperationResp
  CloseOperation_result() = (o=new(); fillunset(o); o)
  CloseOperation_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type CloseOperation_result
meta(t::Type{CloseOperation_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetResultSetMetadata

type GetResultSetMetadata_args
  req::TGetResultSetMetadataReq
  GetResultSetMetadata_args() = (o=new(); fillunset(o); o)
end # type GetResultSetMetadata_args

type GetResultSetMetadata_result
  success::TGetResultSetMetadataResp
  GetResultSetMetadata_result() = (o=new(); fillunset(o); o)
  GetResultSetMetadata_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetResultSetMetadata_result
meta(t::Type{GetResultSetMetadata_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method FetchResults

type FetchResults_args
  req::TFetchResultsReq
  FetchResults_args() = (o=new(); fillunset(o); o)
end # type FetchResults_args

type FetchResults_result
  success::TFetchResultsResp
  FetchResults_result() = (o=new(); fillunset(o); o)
  FetchResults_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type FetchResults_result
meta(t::Type{FetchResults_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method GetDelegationToken

type GetDelegationToken_args
  req::TGetDelegationTokenReq
  GetDelegationToken_args() = (o=new(); fillunset(o); o)
end # type GetDelegationToken_args

type GetDelegationToken_result
  success::TGetDelegationTokenResp
  GetDelegationToken_result() = (o=new(); fillunset(o); o)
  GetDelegationToken_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type GetDelegationToken_result
meta(t::Type{GetDelegationToken_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method CancelDelegationToken

type CancelDelegationToken_args
  req::TCancelDelegationTokenReq
  CancelDelegationToken_args() = (o=new(); fillunset(o); o)
end # type CancelDelegationToken_args

type CancelDelegationToken_result
  success::TCancelDelegationTokenResp
  CancelDelegationToken_result() = (o=new(); fillunset(o); o)
  CancelDelegationToken_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type CancelDelegationToken_result
meta(t::Type{CancelDelegationToken_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# types encapsulating arguments and return values of method RenewDelegationToken

type RenewDelegationToken_args
  req::TRenewDelegationTokenReq
  RenewDelegationToken_args() = (o=new(); fillunset(o); o)
end # type RenewDelegationToken_args

type RenewDelegationToken_result
  success::TRenewDelegationTokenResp
  RenewDelegationToken_result() = (o=new(); fillunset(o); o)
  RenewDelegationToken_result(success) = (o=new(); fillset(o, :success); o.success=success; o)
end # type RenewDelegationToken_result
meta(t::Type{RenewDelegationToken_result}) = meta(t, Symbol[:success], Int[0], Dict{Symbol,Any}())

# Processor for TCLIService service (to be used in server implementation)
type TCLIServiceProcessor <: TProcessor
  tp::ThriftProcessor
  function TCLIServiceProcessor()
    p = new(ThriftProcessor())
    handle(p.tp, ThriftHandler("OpenSession", _OpenSession, OpenSession_args, OpenSession_result))
    handle(p.tp, ThriftHandler("CloseSession", _CloseSession, CloseSession_args, CloseSession_result))
    handle(p.tp, ThriftHandler("GetInfo", _GetInfo, GetInfo_args, GetInfo_result))
    handle(p.tp, ThriftHandler("ExecuteStatement", _ExecuteStatement, ExecuteStatement_args, ExecuteStatement_result))
    handle(p.tp, ThriftHandler("GetTypeInfo", _GetTypeInfo, GetTypeInfo_args, GetTypeInfo_result))
    handle(p.tp, ThriftHandler("GetCatalogs", _GetCatalogs, GetCatalogs_args, GetCatalogs_result))
    handle(p.tp, ThriftHandler("GetSchemas", _GetSchemas, GetSchemas_args, GetSchemas_result))
    handle(p.tp, ThriftHandler("GetTables", _GetTables, GetTables_args, GetTables_result))
    handle(p.tp, ThriftHandler("GetTableTypes", _GetTableTypes, GetTableTypes_args, GetTableTypes_result))
    handle(p.tp, ThriftHandler("GetColumns", _GetColumns, GetColumns_args, GetColumns_result))
    handle(p.tp, ThriftHandler("GetFunctions", _GetFunctions, GetFunctions_args, GetFunctions_result))
    handle(p.tp, ThriftHandler("GetOperationStatus", _GetOperationStatus, GetOperationStatus_args, GetOperationStatus_result))
    handle(p.tp, ThriftHandler("CancelOperation", _CancelOperation, CancelOperation_args, CancelOperation_result))
    handle(p.tp, ThriftHandler("CloseOperation", _CloseOperation, CloseOperation_args, CloseOperation_result))
    handle(p.tp, ThriftHandler("GetResultSetMetadata", _GetResultSetMetadata, GetResultSetMetadata_args, GetResultSetMetadata_result))
    handle(p.tp, ThriftHandler("FetchResults", _FetchResults, FetchResults_args, FetchResults_result))
    handle(p.tp, ThriftHandler("GetDelegationToken", _GetDelegationToken, GetDelegationToken_args, GetDelegationToken_result))
    handle(p.tp, ThriftHandler("CancelDelegationToken", _CancelDelegationToken, CancelDelegationToken_args, CancelDelegationToken_result))
    handle(p.tp, ThriftHandler("RenewDelegationToken", _RenewDelegationToken, RenewDelegationToken_args, RenewDelegationToken_result))
    p
  end
  _OpenSession(inp::OpenSession_args) = OpenSession_result(OpenSession(inp.req))
  _CloseSession(inp::CloseSession_args) = CloseSession_result(CloseSession(inp.req))
  _GetInfo(inp::GetInfo_args) = GetInfo_result(GetInfo(inp.req))
  _ExecuteStatement(inp::ExecuteStatement_args) = ExecuteStatement_result(ExecuteStatement(inp.req))
  _GetTypeInfo(inp::GetTypeInfo_args) = GetTypeInfo_result(GetTypeInfo(inp.req))
  _GetCatalogs(inp::GetCatalogs_args) = GetCatalogs_result(GetCatalogs(inp.req))
  _GetSchemas(inp::GetSchemas_args) = GetSchemas_result(GetSchemas(inp.req))
  _GetTables(inp::GetTables_args) = GetTables_result(GetTables(inp.req))
  _GetTableTypes(inp::GetTableTypes_args) = GetTableTypes_result(GetTableTypes(inp.req))
  _GetColumns(inp::GetColumns_args) = GetColumns_result(GetColumns(inp.req))
  _GetFunctions(inp::GetFunctions_args) = GetFunctions_result(GetFunctions(inp.req))
  _GetOperationStatus(inp::GetOperationStatus_args) = GetOperationStatus_result(GetOperationStatus(inp.req))
  _CancelOperation(inp::CancelOperation_args) = CancelOperation_result(CancelOperation(inp.req))
  _CloseOperation(inp::CloseOperation_args) = CloseOperation_result(CloseOperation(inp.req))
  _GetResultSetMetadata(inp::GetResultSetMetadata_args) = GetResultSetMetadata_result(GetResultSetMetadata(inp.req))
  _FetchResults(inp::FetchResults_args) = FetchResults_result(FetchResults(inp.req))
  _GetDelegationToken(inp::GetDelegationToken_args) = GetDelegationToken_result(GetDelegationToken(inp.req))
  _CancelDelegationToken(inp::CancelDelegationToken_args) = CancelDelegationToken_result(CancelDelegationToken(inp.req))
  _RenewDelegationToken(inp::RenewDelegationToken_args) = RenewDelegationToken_result(RenewDelegationToken(inp.req))
end # type TCLIServiceProcessor
process(p::TCLIServiceProcessor, inp::TProtocol, outp::TProtocol) = process(p.tp, inp, outp)
distribute(p::TCLIServiceProcessor) = distribute(p.tp)

# Server side methods to be defined by user:
# function OpenSession(req::TOpenSessionReq)
#     # returns TOpenSessionResp
# function CloseSession(req::TCloseSessionReq)
#     # returns TCloseSessionResp
# function GetInfo(req::TGetInfoReq)
#     # returns TGetInfoResp
# function ExecuteStatement(req::TExecuteStatementReq)
#     # returns TExecuteStatementResp
# function GetTypeInfo(req::TGetTypeInfoReq)
#     # returns TGetTypeInfoResp
# function GetCatalogs(req::TGetCatalogsReq)
#     # returns TGetCatalogsResp
# function GetSchemas(req::TGetSchemasReq)
#     # returns TGetSchemasResp
# function GetTables(req::TGetTablesReq)
#     # returns TGetTablesResp
# function GetTableTypes(req::TGetTableTypesReq)
#     # returns TGetTableTypesResp
# function GetColumns(req::TGetColumnsReq)
#     # returns TGetColumnsResp
# function GetFunctions(req::TGetFunctionsReq)
#     # returns TGetFunctionsResp
# function GetOperationStatus(req::TGetOperationStatusReq)
#     # returns TGetOperationStatusResp
# function CancelOperation(req::TCancelOperationReq)
#     # returns TCancelOperationResp
# function CloseOperation(req::TCloseOperationReq)
#     # returns TCloseOperationResp
# function GetResultSetMetadata(req::TGetResultSetMetadataReq)
#     # returns TGetResultSetMetadataResp
# function FetchResults(req::TFetchResultsReq)
#     # returns TFetchResultsResp
# function GetDelegationToken(req::TGetDelegationTokenReq)
#     # returns TGetDelegationTokenResp
# function CancelDelegationToken(req::TCancelDelegationTokenReq)
#     # returns TCancelDelegationTokenResp
# function RenewDelegationToken(req::TRenewDelegationTokenReq)
#     # returns TRenewDelegationTokenResp

# Client implementation for TCLIService service
type TCLIServiceClient <: TCLIServiceClientBase
  p::TProtocol
  seqid::Int32
  TCLIServiceClient(p::TProtocol) = new(p, 0)
end # type TCLIServiceClient

# Client callable method for OpenSession
function OpenSession(c::TCLIServiceClientBase, req::TOpenSessionReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "OpenSession", Thrift.MessageType.CALL, c.seqid)
  inp = OpenSession_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, OpenSession_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function OpenSession

# Client callable method for CloseSession
function CloseSession(c::TCLIServiceClientBase, req::TCloseSessionReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "CloseSession", Thrift.MessageType.CALL, c.seqid)
  inp = CloseSession_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, CloseSession_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function CloseSession

# Client callable method for GetInfo
function GetInfo(c::TCLIServiceClientBase, req::TGetInfoReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetInfo", Thrift.MessageType.CALL, c.seqid)
  inp = GetInfo_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetInfo_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetInfo

# Client callable method for ExecuteStatement
function ExecuteStatement(c::TCLIServiceClientBase, req::TExecuteStatementReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "ExecuteStatement", Thrift.MessageType.CALL, c.seqid)
  inp = ExecuteStatement_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, ExecuteStatement_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function ExecuteStatement

# Client callable method for GetTypeInfo
function GetTypeInfo(c::TCLIServiceClientBase, req::TGetTypeInfoReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetTypeInfo", Thrift.MessageType.CALL, c.seqid)
  inp = GetTypeInfo_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetTypeInfo_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetTypeInfo

# Client callable method for GetCatalogs
function GetCatalogs(c::TCLIServiceClientBase, req::TGetCatalogsReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetCatalogs", Thrift.MessageType.CALL, c.seqid)
  inp = GetCatalogs_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetCatalogs_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetCatalogs

# Client callable method for GetSchemas
function GetSchemas(c::TCLIServiceClientBase, req::TGetSchemasReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetSchemas", Thrift.MessageType.CALL, c.seqid)
  inp = GetSchemas_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetSchemas_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetSchemas

# Client callable method for GetTables
function GetTables(c::TCLIServiceClientBase, req::TGetTablesReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetTables", Thrift.MessageType.CALL, c.seqid)
  inp = GetTables_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetTables_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetTables

# Client callable method for GetTableTypes
function GetTableTypes(c::TCLIServiceClientBase, req::TGetTableTypesReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetTableTypes", Thrift.MessageType.CALL, c.seqid)
  inp = GetTableTypes_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetTableTypes_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetTableTypes

# Client callable method for GetColumns
function GetColumns(c::TCLIServiceClientBase, req::TGetColumnsReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetColumns", Thrift.MessageType.CALL, c.seqid)
  inp = GetColumns_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetColumns_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetColumns

# Client callable method for GetFunctions
function GetFunctions(c::TCLIServiceClientBase, req::TGetFunctionsReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetFunctions", Thrift.MessageType.CALL, c.seqid)
  inp = GetFunctions_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetFunctions_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetFunctions

# Client callable method for GetOperationStatus
function GetOperationStatus(c::TCLIServiceClientBase, req::TGetOperationStatusReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetOperationStatus", Thrift.MessageType.CALL, c.seqid)
  inp = GetOperationStatus_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetOperationStatus_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetOperationStatus

# Client callable method for CancelOperation
function CancelOperation(c::TCLIServiceClientBase, req::TCancelOperationReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "CancelOperation", Thrift.MessageType.CALL, c.seqid)
  inp = CancelOperation_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, CancelOperation_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function CancelOperation

# Client callable method for CloseOperation
function CloseOperation(c::TCLIServiceClientBase, req::TCloseOperationReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "CloseOperation", Thrift.MessageType.CALL, c.seqid)
  inp = CloseOperation_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, CloseOperation_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function CloseOperation

# Client callable method for GetResultSetMetadata
function GetResultSetMetadata(c::TCLIServiceClientBase, req::TGetResultSetMetadataReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetResultSetMetadata", Thrift.MessageType.CALL, c.seqid)
  inp = GetResultSetMetadata_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetResultSetMetadata_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetResultSetMetadata

# Client callable method for FetchResults
function FetchResults(c::TCLIServiceClientBase, req::TFetchResultsReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "FetchResults", Thrift.MessageType.CALL, c.seqid)
  inp = FetchResults_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, FetchResults_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function FetchResults

# Client callable method for GetDelegationToken
function GetDelegationToken(c::TCLIServiceClientBase, req::TGetDelegationTokenReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "GetDelegationToken", Thrift.MessageType.CALL, c.seqid)
  inp = GetDelegationToken_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, GetDelegationToken_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function GetDelegationToken

# Client callable method for CancelDelegationToken
function CancelDelegationToken(c::TCLIServiceClientBase, req::TCancelDelegationTokenReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "CancelDelegationToken", Thrift.MessageType.CALL, c.seqid)
  inp = CancelDelegationToken_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, CancelDelegationToken_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function CancelDelegationToken

# Client callable method for RenewDelegationToken
function RenewDelegationToken(c::TCLIServiceClientBase, req::TRenewDelegationTokenReq)
  p = c.p
  c.seqid = (c.seqid < (2^31-1)) ? (c.seqid+1) : 0
  Thrift.writeMessageBegin(p, "RenewDelegationToken", Thrift.MessageType.CALL, c.seqid)
  inp = RenewDelegationToken_args()
  Thrift.set_field!(inp, :req, req)
  Thrift.write(p, inp)
  Thrift.writeMessageEnd(p)
  Thrift.flush(p.t)

  (fname, mtype, rseqid) = Thrift.readMessageBegin(p)
  (mtype == Thrift.MessageType.EXCEPTION) && throw(Thrift.read(p, Thrift.TApplicationException()))
  outp = Thrift.read(p, RenewDelegationToken_result())
  Thrift.readMessageEnd(p)
  (rseqid != c.seqid) && throw(Thrift.TApplicationException(ApplicationExceptionType.BAD_SEQUENCE_ID, "response sequence id $rseqid did not match request ($(c.seqid))"))
  Thrift.has_field(outp, :success) && (return Thrift.get_field(outp, :success))
  throw(Thrift.TApplicationException(Thrift.ApplicationExceptionType.MISSING_RESULT, "retrieve failed: unknown result"))
end # function RenewDelegationToken


text:
export Jackknife, JackknifeSet
export jackknife

mutable struct Jackknife <: ScalarObservable
    xs :: Vector{Float64}
end

Jackknife(jk::Jackknife, f::Function) = Jackknife(map(f,jk.xs))

function jk_helper(xs)
    s = sum(xs)
    n = length(xs)-1
    ret = similar(xs)
    for i in 1:n+1
        ret[i] = (s-xs[i])/n
    end
    return ret
end

function Jackknife(o::ScalarObservable)
    return count(o) > 0 ? Jackknife(jk_helper(o.bins)) : Jackknife(zeros(0))
end
function Jackknife(b::BinningObservable)
    if count(b) > 0
        return Jackknife(jk_helper(ifelse(b.lastbin<b.binsize,
                                          b.bins[1:end-1],
                                          b.bins)))
    else
        return Jackknife(zeros(0))
    end
end

count(jk::Jackknife) = length(jk.xs)

mean(jk::Jackknife) = mean(jk.xs)
function stderror(jk::Jackknife)
    n = count(jk)
    if n < 2
        return NaN
    else
        m2 = sum(abs2, jk.xs)
        m2 /= n
        m = mean(jk)
        sigma2 = m2 - m*m
        sigma2 *= n-1
        sigma2 = maxzero(sigma2)
        return sqrt(sigma2)
    end
end
function var(jk::Jackknife)
    n = count(jk)
    if n < 2
        return NaN
    else
        m = mean(jk)
        return sum(abs2, jk.xs.-m)*(n-1)
    end
end
stddev(jk::Jackknife) = sqrt(var(jk))

function confidence_interval(jk::Jackknife, confidence_rate::Real)
    q = 0.5+0.5*confidence_rate
    correction = quantile( TDist(count(jk)-1), q)
    serr = stderror(jk)
    return correction * serr
end
function confidence_interval(jk::Jackknife, confidence_rate_symbol::Symbol = :sigma1)
    n = parsesigma(confidence_rate_symbol)
    return confidence_interval(jk, erf(0.5n*sqrt(2.0)))
end

unary_functions = (
                   :-,
                   :sin, :cos, :tan,
                   :sind, :cosd, :tand,
                   :sinpi, :cospi,
                   :sinh, :cosh, :tanh,
                   :asin, :acos, :atan,
                   :asind, :acosd, :atand,
                   :sec, :csc, :cot,
                   :secd, :cscd, :cotd,
                   :asec, :acsc, :acot,
                   :asecd, :acscd, :acotd,
                   :sech, :csch, :coth,
                   :asinh, :acosh, :atanh,
                   :asech, :acsch, :acoth,
                   :sinc, :cosc,
                   :log, :log2, :log10, :log1p,
                   :exp, :exp2, :exp10, :expm1,
                   :abs, :abs2,
                   :sqrt, :cbrt,
                  )

for op in unary_functions
    @eval Base.$op(jk::Jackknife) = Jackknife(jk, $op)
end

binary_functions = (
                    :+, :-, :*, :/, :\,
                   )

for op in binary_functions
    @eval Base.$op(jk::Jackknife, rhs::Real) = Jackknife(jk, lhs->($op)(lhs,rhs))
    @eval Base.$op(lhs::Real, jk::Jackknife) = Jackknife(jk, rhs->($op)(lhs,rhs))
    @eval Base.$op(lhs::Jackknife, rhs::Jackknife) = Jackknife(broadcast($op, lhs.xs, rhs.xs))
end

import Base.^
^(lhs::Jackknife, rhs::Real) = Jackknife((lhs.xs).^rhs)
^(lhs::Jackknife, rhs::Integer) = Jackknife((lhs.xs).^rhs)
^(lhs::Real, rhs::Jackknife) = Jackknife(lhs.^(rhs.xs))
^(lhs::Integer, rhs::Jackknife) = Jackknife(lhs.^(rhs.xs))
^(lhs::Jackknife, rhs::Jackknife) = Jackknife((lhs.xs).^(rhs.xs))

const JackknifeSet = MCObservableSet{Jackknife}

jackknife(obs::ScalarObservable) = Jackknife(obs)
function jackknife(obsset :: MCObservableSet)
    JK = JackknifeSet()
    for (k,v) in obsset
        JK[k] = Jackknife(v)
    end
    return JK
end


text:
using ClangCompiler
using Test

if haskey(ENV, "CI")
    include("clang/status.jl")
    include("platform/env.jl")
end

include("clang/diagnostic.jl")
include("clang/source.jl")
include("clang/instance.jl")

include("lookup.jl")
if Sys.isapple()
    include("call.jl")
    include("execution.jl")
end
include("types.jl")

if haskey(ENV, "CLANGCOMPILER_ENABLE_BOOT")
    include("boot.jl")
end


text:
using FuzzySets, Test

x = [1, 4.5]
y = [3.5, 4]
@test FuzzySets.distance(x, y) == 2.5495097567963922


text:
# Load Julia packages (libraries) needed

using StatisticalRethinking

# ### snippet 2.1

ways  = [0, 3, 8, 9, 0];

#-

ways/sum(ways)

# ### snippet 2.2

# Create a distribution with n = 9 (e.g. tosses) and p = 0.5.

d = Binomial(9, 0.5)

# Probability density for 6 `waters` holding n = 9 and p = 0.5.

pdf(d, 6)

# End of `02/clip-01-02.jl`


text:
struct DistanceFromSphereSurface{M,R1,S,m} <: TaskMap{M,R1,S}
    center::SVector{m,S}
    radius::S
end
DistanceFromSphereSurface{M,R1}(c::SVector{m,S}, r::S) where {M,R1,S,m} =
    DistanceFromSphereSurface{M,R1,S,embdim(M)}(c, r)

struct DistanceFromSphereSurfaceT{M,R1,S,m} <: TaskMapT{M,R1,S}
    base_map::DistanceFromSphereSurface{M,R1,S,m}
end
DistanceFromSphereSurfaceT{M,R1}(c::SVector{m,S}, r::S) where {M,R1,S,m} =
    DistanceFromSphereSurfaceT{M,R1,S,m}(DistanceFromSphereSurface{M,R1}(c, r))

task_map_emb(::EmbRep, ::EmbRep, xme, task_map::DistanceFromSphereSurface{ℝ{m},R1}) where m =
    SA[abs(norm(xme - task_map.center) - task_map.radius)]


text:
using HMCExamples
# The main purpose is for tracing the compiler, not for exection speed or results.  Keep the number of samples small
HMCExamples.main_rbc_1_kalman(["--num_samples", "100"])
HMCExamples.main_rbc_2_joint(["--num_samples", "100"])