text:
using SyntaxTree
using SymEngine
using CUDA
include("types.jl")

function GenStencilKernel(st_sym, max_radius,prev_time_coeff, bdim=32)

    function search_edges(stencil::Basic, max_radius::Integer, dim::Integer)
        radius = max_radius
        Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2*max_radius+1),
                       j = 1:(2*max_radius+1), k=1:(2*max_radius+1)]
        Doff = [[i,j,k] for i = -radius:radius, j=-radius:radius, k=-radius:radius]
        mn = 100
        mx = -100
        for i = 1:(2*max_radius+1), j = 1:(2*max_radius+1),k = 1:(2*max_radius+1)
            cf = SymEngine.coeff(stencil,Dsym[i,j,k])
            of = Doff[i,j,k]
            if cf != 0
                mn > of[dim] ? (mn= of[dim]) : nothing
                mx < of[dim] ? (mx= of[dim]) : nothing
            end
        end
        return (mn,mx)
    end
    #bdim = 32
    mnx, mxx = search_edges(st_sym, max_radius, 1)
    mny, mxy = search_edges(st_sym, max_radius, 2)
    mnz, mxz = search_edges(st_sym, max_radius, 3)
    shmem_radius = max_radius#max( abs.((mnx, mxx, mny, mxy))... )
    v_info, nv_info = GenerateStencil(st_sym, shmem_radius)
    q = CreateInitPart(shmem_radius, bdim, v_info, nv_info)
    l,b1 = CreateLoopPart(shmem_radius, bdim)
    b2,cinit = CreateCalcsPart(v_info, nv_info, shmem_radius,prev_time_coeff=prev_time_coeff)
    b3 = CreateRotateRegsPart(v_info, nv_info)
    f = CreateFinalStorePart(v_info, nv_info,prev_time_coeff=prev_time_coeff)
    for i in b2.args
        push!(b1.args, i)
    end
    for i in b3.args
        push!(b1.args, i)
    end
    push!(l.args, b1)
    push!(q.args, cinit)
    push!(q.args, l)
    for i in f.args
        push!(q.args,i)
    end
    return linefilter!(q), v_info.exists, abs(mxz), abs(mnz)
end

"""
    search_edges(stencil, max_radius, dim)

    For a `stencil`'s symbolic math expression, find how much it extends
    to the given dimesnion `dim`. `max_radius` is used to create the array of
    math symbols.

    Returns the minimum and maximum index
"""
function search_edges(stencil::Basic, max_radius::Integer, dim::Integer)
    radius = max_radius
    Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2*max_radius+1),
                   j = 1:(2*max_radius+1), k=1:(2*max_radius+1)]
    Doff = [[i,j,k] for i = -radius:radius, j=-radius:radius, k=-radius:radius]
    mn = 100
    mx = -100
    for i = 1:(2*max_radius+1), j = 1:(2*max_radius+1),k = 1:(2*max_radius+1)
        cf = SymEngine.coeff(stencil,Dsym[i,j,k])
        of = Doff[i,j,k]
        # if dim == 3 && cf != 0
        #     @show cf, of
        # end
        if cf != 0
            mn > of[dim] ? (mn= of[dim]) : nothing
            mx < of[dim] ? (mx= of[dim]) : nothing
        end
    end
    return (mn,mx)
end

function GenerateStencil(st_sym::Basic, radius::Integer)
    ## Check if it uses vsq array
    check_vsq = coeff(st_sym, v) != 0
    vsq_stencil = nothing
    non_vsq_stencil = st_sym
    ## Split vsq stencil and non vsq stencil
    if check_vsq
        vsq_stencil = coeff(st_sym,v)
        non_vsq_stencil = expand(st_sym - v*vsq_stencil)
    end
    ## ----------------------------------------
    #           FIND INFO FOR EACH STENCIL
    ## ----------------------------------------
    vsq_info = stencil_info()
    if !check_vsq
        vsq_info.exists = false
    else
        vsq_info.exists = true
        vsq_info.sym_math = vsq_stencil
        vsq_info.max_radius = radius
        mnz, mxz = search_edges(vsq_stencil, radius, 3)
        # front cannot be negative
        vsq_info.front_max = mnz > 0 ? 0 : abs(mnz)
        # behind can't be positive
        vsq_info.behind_max = mxz < 0 ? 0 : mxz
    end
    non_vsq_info = stencil_info()
    if non_vsq_stencil == 0
        non_vsq_info.exists = false
    else
        non_vsq_info.exists = true
        non_vsq_info.sym_math = non_vsq_stencil
        non_vsq_info.max_radius = radius
        mnz, mxz = search_edges(non_vsq_stencil, radius, 3)
        # front cannot be negative
        non_vsq_info.front_max = mnz > 0 ? 0 : abs(mnz)
        # behind can't be positive
        non_vsq_info.behind_max = mxz < 0 ? 0 : mxz
    end

    ## FIX front/behind_max to match both stencils
    if vsq_info.exists == non_vsq_info.exists == true
        # data have to be passed from front to behind
        # until they are stored. So both stencil must
        # have the same amount of behind registers.
        if non_vsq_info.behind_max > vsq_info.behind_max
            vsq_info.behind_max = non_vsq_info.behind_max
        else
            non_vsq_info.behind_max = vsq_info.behind_max
        end
    end
    return (vsq_info, non_vsq_info)
end

# function convert_to_sym(st::StencilDefinition, radius::Integer)
#     if st.expr isa Nothing
#         ## For coef stencil
#         coefs = copy(st.coefs)
#         if st.max_radius < radius
#             pad = radius - st.max_radius
#             coefs = zeros(Float64, 2*radius+1, 2*radius+1, 2*radius+1)
#             # for i in  length(coefs)
#             #     coefs[i] = (0)
#             #     println(coefs)
#             # end
#             # println(coefs)
#             #println(size(coefs))
#             for i = 1:(2st.max_radius+1), j = 1:(2st.max_radius+1), k = 1:(2st.max_radius+1)
#                 coefs[pad+i, pad+j,pad+k] = st.coefs[i,j,k]
#                 #println(size(coefs))
#             end
#         end
#         #println(size(coefs))
#         global D
#         D = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2radius+1),
#                                             j = 1:(2radius+1),
#                                             k = 1:(2radius+1)]
#         st_sym = sum(coefs.*D)
#         if st.usus_vsq
#             global v
#             st_sym *= v
#         end
#         return st_sym
#     else
#         ## For expression stencil
#         coefs = copy(st.coefs[:])
#         expr = copy(st.expr)
#         for i in [:x :y :z]
#             expr = replace_symbol(i, 0, expr)
#         end
#         if length(coefs) < radius
#             pad = radius - length(coefs)
#             for i = 1:pad
#                 push!(coefs, 0)
#             end
#         end
#         global c = [reverse(coefs[2:end]); coefs[:]]
#         expr = fix_indices(expr, radius+1)
#         st_sym = eval(expr)
#         return st_sym
#     end
# end

function ConvertToFunction(ex::Expr, uses_vsq::Bool)
    if uses_vsq
        return SyntaxTree.genfun(ex, [:g_input, :g_output, :offz_f, :offz_b, :offx, :offy, :g_vsq])
    else
        return SyntaxTree.genfun(ex, [:g_input, :g_output, :offz_f, :offz_b, :offx, :offy])
    end
end

function CreateInitPart(shmem_radius::Integer, bdim::Integer,
                vsq_st::stencil_info, non_vsq_st::stencil_info)
    radius = shmem_radius
    temp = bdim + 2*radius
    q = quote
        dimz = size(g_input,3)
        dimzend = dimz-$(radius)
        bdimx = $(CUDA).blockDim().x
        bdimy = $(CUDA).blockDim().y
        bx = ($(CUDA).blockIdx().x-1)*$(bdim) + 1 + offx
        by = ($(CUDA).blockIdx().y-1)*$(bdim) + 1 + offy
        tx = $(CUDA).threadIdx().x
        ty = $(CUDA).threadIdx().y
        txr = $(CUDA).threadIdx().x+$(radius)
        tyr = $(CUDA).threadIdx().y+$(radius)
        tile = $(CUDA).@cuDynamicSharedMem(Float32,
                                   ($(temp), $(temp)))

        @inbounds l_input = @view g_input[bx:(bx+$(temp-1)),
                               by:(by+$(temp-1)), :]
        @inbounds l_output = @view g_output[bx:(bx+$(temp-1)),
                               by:(by+$(temp-1)), :]
    end
    if non_vsq_st.exists
        for i = non_vsq_st.behind_max:(-1):1
            sym = Symbol(string("behind", i))
            push!(q.args, :($(sym) = Float32(0.0)))
        end
        push!(q.args, :(current = Float32(0.0)))
        for i = 1:non_vsq_st.front_max
            sym = Symbol(string("infront", i))
            push!(q.args, :($(sym) = Float32(0.0)))
        end
    end
    """Add the init for non vsq registers
        vsq_infront max is at least current
        vsq_behind max is as normal begind max
    """
    if vsq_st.exists
        push!(q.args, :(@inbounds l_vsq = @view g_vsq[bx:(bx+$(bdim-1)),
                                   by:(by+$(bdim-1)), :]))
        for i = vsq_st.behind_max:(-1):1
            sym = Symbol(string("vsq_behind", i))
            push!(q.args, :($(sym) = Float32(0.0)))
        end
        push!(q.args, :(vsq_current = Float32(0.0)))
        for i = 1:vsq_st.front_max
            sym = Symbol(string("vsq_infront", i))
            push!(q.args, :($(sym) = Float32(0.0)))
        end
    end
    return linefilter!(q)
end

function CreateCalcsPart(vsq_st, non_vsq_st, radius; prev_time_coeff=0)
    coeff_arr = Vector{Float32}(undef,0)
    function check_coeff(c_arr, c)
        ind = findfirst(x-> x==Float32(c), c_arr)
        if ind isa Nothing
            push!(c_arr, Float32(c))
            ind = length(c_arr)
        end
        return :(cfarr[$(ind)])

    end
    function init_coeff_arr(cf)
        ex = :(cfarr = $(StaticArrays).SA_F32[])
        for i in cf
            push!(ex.args[2].args , i)
        end
        return ex
    end

    max_radius = radius
    Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2*max_radius+1),
                   j = 1:(2*max_radius+1), k=1:(2*max_radius+1)]

    Doff = [[i,j,k] for i = -radius:radius, j=-radius:radius, k=-radius:radius]
    q = Expr(:block)
    # ## Determine Calculation for the 2d xy Tile
    # for  k = radius+1,i = 1:size(Dsym,1), j = 1:size(Dsym,2)
    #     pos = Doff[i,j,k]
    #     sym = Dsym[i,j,k]
    #     if vsq_st.exists
    #         cf = Float64(coeff(vsq_st.sym_math, sym))
    #         if cf != 0
    #             ex = :()
    #             if cf == 1
    #                 ex = :(@inbounds vsq_current += tile[txr+$(pos[1]), tyr + $(pos[2])])
    #             else
    #                 ex = :(@inbounds vsq_current += Float32($(cf))*tile[txr+$(pos[1]), tyr + $(pos[2])])
    #             end
    #             push!(q.args, ex)
    #         end
    #     end
    #     if non_vsq_st.exists
    #         cf = Float64(coeff(non_vsq_st.sym_math, sym))
    #         if cf != 0
    #             ex = :()
    #             if cf == 1
    #                 ex = :(@inbounds current += tile[txr+$(pos[1]), tyr + $(pos[2])])
    #             else
    #                 ex = :(@inbounds current += Float32($(cf))*tile[txr+$(pos[1]), tyr + $(pos[2])])
    #             end
    #             push!(q.args, ex)
    #         end
    #     end
    # end

    ## Determine Calculation for contribution at infront and behind
    for x = 1:(2*max_radius+1), y = 1:(2*max_radius+1)
        of = Doff[x,y,1]
        read_to_temp = :(@inbounds temp = tile[txr+$(of[1]), tyr+$(of[2])])
        flag_read_temp = true

        # Current
        pos = Doff[x,y,max_radius+1]
        sym = Dsym[x,y,max_radius+1]
        if vsq_st.exists
            cf = Float64(coeff(vsq_st.sym_math, sym))
            if cf != 0
                if flag_read_temp
                    push!(q.args, read_to_temp)
                    flag_read_temp = false
                end
                ex = :()
                if cf == 1
                    ex = :(vsq_current += temp)
                else
                    csym = check_coeff(coeff_arr, Float32(cf))
                    ex = :(vsq_current += $(csym)*temp)
                end
                push!(q.args, ex)
            end
        end
        if non_vsq_st.exists
            cf = Float64(coeff(non_vsq_st.sym_math, sym))
            if cf != 0
                if flag_read_temp
                    push!(q.args, read_to_temp)
                    flag_read_temp = false
                end
                ex = :()
                if cf == 1
                    ex = :(current += temp)
                else
                    csym = check_coeff(coeff_arr, Float32(cf))
                    ex = :(current += $(csym)*temp)
                end
                push!(q.args, ex)
            end
        end

        for z = 1:max_radius
        ## FRONT
            var = Dsym[x,y,z]
            off = Doff[x,y,z]
            if vsq_st.exists
                sym = Symbol(string("vsq_infront",abs(off[3])))
                cf = Float64(coeff(vsq_st.sym_math, var))
                if cf != 0
                    if flag_read_temp
                        push!(q.args, read_to_temp)
                        flag_read_temp = false
                    end
                    if cf == 1
                        push!(q.args, :($(sym) += temp))
                    else
                        csym = check_coeff(coeff_arr, Float32(cf))
                        push!(q.args, :($(sym) += $(csym)*temp))
                    end
                end
            end

            if non_vsq_st.exists
                sym = Symbol(string("infront",abs(off[3])))
                cf = Float64(coeff(non_vsq_st.sym_math, var))
                if cf != 0
                    if flag_read_temp
                        push!(q.args, read_to_temp)
                        flag_read_temp = false
                    end
                    if cf == 1
                        push!(q.args, :($(sym) += temp))
                    else
                        csym = check_coeff(coeff_arr, Float32(cf))
                        push!(q.args, :($(sym) += $(csym)*temp))
                    end
                end
            end
        end
        for z = 1:max_radius
        ## Behind
            var = Dsym[x,y,2*max_radius + 2 - z]
            off = Doff[x,y,2*max_radius + 2 - z]
            if vsq_st.exists
                sym = Symbol(string("vsq_behind",abs(off[3])))
                cf = Float64(coeff(vsq_st.sym_math, var))
                if cf != 0
                    if flag_read_temp
                        push!(q.args, read_to_temp)
                        flag_read_temp = false
                    end
                    if cf == 1
                        push!(q.args, :($(sym) += temp))
                    else
                        csym = check_coeff(coeff_arr, Float32(cf))
                        push!(q.args, :($(sym) += $(csym)*temp))
                    end
                end
            end

            if non_vsq_st.exists
                sym = Symbol(string("behind",abs(off[3])))
                cf = Float64(coeff(non_vsq_st.sym_math, var))
                if cf != 0
                    if flag_read_temp
                        push!(q.args, read_to_temp)
                        flag_read_temp = false
                    end
                    if cf == 1
                        push!(q.args, :($(sym) += temp))
                    else
                        csym = check_coeff(coeff_arr, Float32(cf))
                        push!(q.args, :($(sym) += $(csym)*temp))
                    end
                end
            end
        end
    end

    ## SAVE RESULTS
    if non_vsq_st.exists == vsq_st.exists == true
        if non_vsq_st.behind_max == 0
            if prev_time_coeff == 0
                push!(q.args, :(@inbounds l_output[txr,tyr,z] =
                                         vsq_current*l_vsq[tx,ty,z] + current))
            else
                push!(q.args, :(@inbounds l_output[txr,tyr,z] =
                                    vsq_current*l_vsq[tx,ty,z] +
                                    current + Float32($(prev_time_coeff))*l_output[txr,tyr,z]))
            end
        else
            bmax = non_vsq_st.behind_max
            snv = Symbol(string("behind", non_vsq_st.behind_max))
            sv = Symbol(string("vsq_behind", non_vsq_st.behind_max))
            if prev_time_coeff == 0
                push!(q.args, :(if z > $(bmax)
                                @inbounds l_output[txr,tyr,z-$(bmax)] = $(sv)*l_vsq[tx,ty,z-$(bmax)] + $(snv)
                            end))
            else
                push!(q.args, :(if z > $(bmax)
                                @inbounds l_output[txr,tyr,z-$(bmax)] = $(sv)*l_vsq[tx,ty,z-$(bmax)] + $(snv) + Float32($(prev_time_coeff)*l_output[txr,tyr,z-$(bmax)])
                            end))
            end

        end

    else
        if non_vsq_st.exists
            if non_vsq_st.behind_max == 0
                if prev_time_coeff == 0
                    push!(q.args, :(
                                    @inbounds l_output[txr,tyr,z] = current))
                else
                    push!(q.args, :(
                                    @inbounds l_output[txr,tyr,z] = current + Float32($(prev_time_coeff))*l_output[txr,tyr,z]
                                ))
                end
            else
                bmax = non_vsq_st.behind_max
                sym = Symbol(string("behind", non_vsq_st.behind_max))
                if prev_time_coeff == 0
                    push!(q.args, :(if z > $(bmax)
                                    @inbounds l_output[txr,tyr,z-$(bmax)] = $(sym)
                                end))
                else
                    push!(q.args, :(if z > $(bmax)
                                    @inbounds l_output[txr,tyr,z-$(bmax)] = $(sym) + Float32($(prev_time_coeff))*l_output[txr,tyr,z-$(bmax)]
                                end))

                end
            end
        end

        if vsq_st.exists
            if vsq_st.behind_max == 0
                if prev_time_coeff == 0
                    push!(q.args, :(
                                    @inbounds l_output[txr,tyr,z] = vsq_current
                                ))
                else
                    push!(q.args, :(
                                    @inbounds l_output[txr,tyr,z] = vsq_current + Float32($(prev_time_coeff))*l_output[txr,tyr,z]
                                ))
                end
            else
                bmax = vsq_st.behind_max
                sym = Symbol(string("vsq_behind", vsq_st.behind_max))
                if prev_time_coeff == 0
                    push!(q.args, :(if z > $(bmax)
                                    @inbounds l_output[txr,tyr,z-$(bmax)] = $(sym)*l_vsq[tx,ty,z-$(bmax)]
                                end))
                else
                    push!(q.args, :(if z > $(bmax)
                                    @inbounds l_output[txr,tyr,z-$(bmax)] = $(sym)*l_vsq[tx,ty,z-$(bmax)] + Float32($(prev_time_coeff))*l_output[txr,tyr,z-$(bmax)]
                                end))
                end
            end
        end

    end
    return linefilter!(q), init_coeff_arr(coeff_arr)
end

function CreateLoopPart(radius, bdim)
    #radius = max(non_vsq_st.behind_max, vsq_st.behind_max)
    l = Expr(:for)
    push!(l.args, :(z = ((1+offz_f)):((dimz-offz_b))))

    q = Expr(:block)
    push!(q.args, :($(CUDA).sync_threads()))

    len = bdim + 2*radius
    for i = 0:bdim:len
        for j = 0:bdim:len
            tq = :(
                 if tx+$(i) <= $(len) && ty+$(j) <= $(len)
                     @inbounds tile[tx+$(i), ty+$(j)] = l_input[tx+$(i), ty+$(j),z]
                 end
            )

            push!(q.args, tq)
        end
    end
    push!(q.args, :(($CUDA).sync_threads()))
    #push!(l.args, q)
    return l,q
end

function CreateRotateRegsPart(vsq_st, non_vsq_st)
    q = Expr(:block)
    if non_vsq_st.exists
        regs = []

        for i = non_vsq_st.behind_max:-1:1
            push!(regs, Symbol(string("behind", i)))
        end
        push!(regs, :current)
        for i = 1:non_vsq_st.front_max
            push!(regs, Symbol(string("infront",i )))
        end
        for i = 1:(length(regs) - 1)
            push!(q.args, :($(regs[i]) = $(regs[i+1])))
        end
        sym = non_vsq_st.front_max == 0 ? :current : Symbol(string("infront", non_vsq_st.front_max))
        push!(q.args, :($(sym) = Float32(0)))
    end
    if vsq_st.exists
        regs = []

        for i = vsq_st.behind_max:-1:1
            push!(regs, Symbol(string("vsq_behind", i)))
        end
        push!(regs, :vsq_current)
        for i = 1:vsq_st.front_max
            push!(regs, Symbol(string("vsq_infront",i )))
        end
        for i = 1:(length(regs) - 1)
            push!(q.args, :($(regs[i]) = $(regs[i+1])))
        end
        sym = vsq_st.front_max == 0 ? :vsq_current : Symbol(string("vsq_infront", vsq_st.front_max))
        push!(q.args, :($(sym) = Float32(0)))
    end
    return q
end

function CreateFinalStorePart(vsq_st, non_vsq_st; prev_time_coeff = 0)
    q = Expr(:block)
    behindz = max(vsq_st.behind_max, non_vsq_st.behind_max)
    behindz == 0 && return :(return nothing)

    for i = 1:behindz
        ## SAVE RESULTS
        if non_vsq_st.exists == vsq_st.exists == true
            bmax = i
            snv = Symbol(string("behind", i))
            sv = Symbol(string("vsq_behind", i))
            if prev_time_coeff == 0
                push!(q.args, :(@inbounds l_output[txr,tyr,dimz-$(bmax-1)-offz_b] = $(sv)*l_vsq[tx,ty,dimz-$(bmax-1)-offz_b] + $(snv)))
            else
                push!(q.args, :(@inbounds l_output[txr,tyr,dimz-$(bmax-1)-offz_b] = $(sv)*l_vsq[tx,ty,dimz-$(bmax-1)-offz_b] + $(snv) + Float32($(prev_time_coeff))*l_output[txr,tyr,dimz-$(bmax-1)-offz_b]))
            end
        else
            if non_vsq_st.exists
                bmax = i#non_vsq_st.behind_max
                sym = Symbol(string("behind", i))
                if prev_time_coeff == 0
                    push!(q.args, :(@inbounds l_output[txr,tyr,dimz-$(bmax-1)-offz_b] = $(sym)))
                else
                    push!(q.args, :(@inbounds l_output[txr,tyr,dimz-$(bmax-1)-offz_b] = $(sym) + Float32($(prev_time_coeff)*l_output[txr,tyr,dimz-$(bmax-1)-offz_b])))

                end
            end

            if vsq_st.exists
                bmax = i#non_vsq_st.behind_max
                sym = Symbol(string("vsq_behind", i))
                if prev_time_coeff == 0
                    push!(q.args, :(@inbounds l_output[txr,tyr,dimz-$(bmax-1)-offz_b] = $(sym)*l_vsq[tx,ty,dimz-$(bmax-1)-offz_b]))
                else
                    push!(q.args, :(@inbounds l_output[txr,tyr,dimz-$(bmax-1)-offz_b] = $(sym)*l_vsq[tx,ty,dimz-$(bmax-1)-offz_b] + Float32($(prev_time_coeff))*l_output[txr,tyr,dimz-$(bmax-1)-offz_b]))
                end
            end

        end
    end

    push!(q.args, :(return nothing))
    return q
end


text:
using GalacticOptim, Optim, Test, ReverseDiff, Tracker

x0 = zeros(2)
rosenbrock(x, p=nothing) =  (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
l1 = rosenbrock(x0)

function g!(G, x)
    G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    G[2] = 200.0 * (x[2] - x[1]^2)
end

function h!(H, x)
    H[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
    H[1, 2] = -400.0 * x[1]
    H[2, 1] = -400.0 * x[1]
    H[2, 2] = 200.0
end

G1 = Array{Float64}(undef,2)
G2 = Array{Float64}(undef,2)
H1 = Array{Float64}(undef, 2, 2)
H2 = Array{Float64}(undef, 2, 2)

g!(G1, x0)
h!(H1, x0)

optf = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
optprob = GalacticOptim.instantiate_function(optf,x0,GalacticOptim.AutoForwardDiff(),nothing)
optprob.grad(G2, x0)
@test G1 == G2
optprob.hess(H2, x0)
@test H1 == H2

prob = OptimizationProblem(optprob, x0)

sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

sol = solve(prob, Newton())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.KrylovTrustRegion())
@test 10*sol.minimum < l1

optf = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())
optprob = GalacticOptim.instantiate_function(optf,x0,GalacticOptim.AutoZygote(),nothing)
optprob.grad(G2, x0)
@test G1 == G2
optprob.hess(H2, x0)
@test H1 == H2

prob = OptimizationProblem(optprob, x0)

sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

sol = solve(prob, Newton())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.KrylovTrustRegion())
@test 10*sol.minimum < l1

optf = OptimizationFunction(rosenbrock, GalacticOptim.AutoReverseDiff())
optprob = GalacticOptim.instantiate_function(optf,x0,GalacticOptim.AutoReverseDiff(),nothing)
optprob.grad(G2, x0)
@test G1 == G2
optprob.hess(H2, x0)
@test H1 == H2

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

sol = solve(prob, Newton())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.KrylovTrustRegion())
@test 10*sol.minimum < l1

optf = OptimizationFunction(rosenbrock, GalacticOptim.AutoTracker())
optprob = GalacticOptim.instantiate_function(optf,x0,GalacticOptim.AutoTracker(),nothing)
optprob.grad(G2, x0)
@test G1 == G2
@test_throws ErrorException optprob.hess(H2, x0)

prob = OptimizationProblem(optprob, x0)

sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

@test_throws ErrorException solve(prob, Newton())

optf = OptimizationFunction(rosenbrock, GalacticOptim.AutoFiniteDiff())
optprob = GalacticOptim.instantiate_function(optf,x0,GalacticOptim.AutoFiniteDiff(),nothing)
optprob.grad(G2, x0)
@test G1 ≈ G2 rtol=1e-6
optprob.hess(H2, x0)
@test H1 ≈ H2 rtol=1e-6

prob = OptimizationProblem(optprob, x0)
sol = solve(prob, BFGS())
@test 10*sol.minimum < l1

sol = solve(prob, Newton())
@test 10*sol.minimum < l1

sol = solve(prob, Optim.KrylovTrustRegion())
@test sol.minimum < l1 #the loss doesn't go below 5e-1 here

ext: jl


text:
import OpenNI
using OpenNI

function distance(pnt1, pnt2, norm=2)
    sum = 0
    for i=1:length(pnt1)
        sum += (pnt1[i] - pnt2[i]) ^ norm
    end
    return sum ^ (1 / norm)
end

function angleBetweenVectors(pnt1, pnt2, pnt3)
    c = distance(pnt3, pnt1)
    a = distance(pnt1, pnt2)
    b = distance(pnt2, pnt3)

    try
        return acos((a ^ 2 + b ^ 2 - c ^ 2) /
         (a * b * 2))
    catch
        return pi
    end
end

function convertKinectPntToSimPnt(kinectPnt::Vector{NIFloat}, origin::Vector{NIFloat})
    for i=1:length(origin)
        kinectPnt[i] -= origin[i]
        #now scale down, first convert to inches, then to game points
        kinectPnt[i] *= 0.0393701 / .27
    end
    kinectPnt[2] += 70 #add to y? TODO way to automate
    kinectPnt[3] *= -1 #flip z
    return kinectPnt
end

function convertSkeletonToSim(jointPos::Matrix{NIFloat}, origin::Vector{NIFloat})
    for i=1:length(jointPos[:, 1])
        jointPos[i, :] = convertKinectPntToSimPnt(jointPos[i, :], origin)
    end
end

function delimVector(pnt)
    str = ""
    for i=1:length(pnt)
        str = "$str$(pnt[i]),"
    end
    return str
end

###Record only once skeleton is detected and then the user puts their right hand by torso

g_elapseSinceLastEvent = 0
g_userFound = false
g_startPlaying = false
g_stopPlaying = false
g_origin = Vector{NIFloat}()

#g_outputStream = STDOUT
g_outputStream = open(joinpath(@__DIR__, "../../data/csv/$(ARGS[1]).csv"), true, true, true, true, false)

#Init without any preloaded nodes
print("Creating context...")

context = NIContext()
InitContext(context)

println("Done!")
print("Creating Player...")

player = NIPlayer()
OpenContextFileRecording(player, context,  joinpath(@__DIR__, "../../data/$(ARGS[1]).oni"))

println("Done!")
print("Attaching generators...")

#Now create a UserGenerator node
usergen = NIUserGenerator()
CreateUserGenerator(usergen, context)

#Function for callbacks (new and lost user)
function userFound()
    println("User recognized registered!")
end

function userLost()
    g_stopPlaying = true
    println("User lost. Stopping recording")
end

RegisterUserCallbacks(usergen, userFound, userLost)

println("Done!")

#Begin generators
StartGeneratingAll(context)

print("Starting playback...")

#init array of user IDs
userIDs = Vector{UInt32}(2)

tic()
#start program loop
while !g_stopPlaying
    #each frame update, but wait for usergen to have new data
    WaitOneUpdateAll(context, usergen)

    userIDs = GetUsers(usergen, userIDs)

    for i=1:length(userIDs)
        #If tracking, print head coordinates
        if (IsTrackingUser(usergen, userIDs[i]))

            skeleton = NISkeleton(usergen, userIDs[i])
            jntHandID = Int(OpenNI.NI_SKEL_LEFT_HAND) #mirrored
            jntTorsoID = Int(OpenNI.NI_SKEL_TORSO)
            dist = distance(skeleton.jointPos[jntHandID, :], skeleton.jointPos[jntTorsoID, :])

            g_elapseSinceLastEvent += toq()
            tic()

            if (!g_userFound && dist < 250)
                g_userFound = true
                print("User struck pose, beginning motion in 2 seconds...")

                #calculate origin
                g_origin = deepcopy(skeleton.jointPos[jntHandID, :])

                #adjust the y? TODO
                #adjust the z
                g_origin[3] -= 55 / 0.0393701 * .27

                println("Origin:", g_origin)

                g_elapseSinceLastEvent = 0
                tic()
            end

            if (g_userFound && !g_startPlaying && g_elapseSinceLastEvent > 2)
                println("Starting motion...")
                g_startPlaying = true
                g_elapseSinceLastEvent = 0
                tic()
            end

            if (g_startPlaying)
                println(g_outputStream, delimVector(convertKinectPntToSimPnt(deepcopy(skeleton.jointPos[jntHandID, :]), g_origin)))
            end

            if (g_startPlaying && g_elapseSinceLastEvent > 5)
                println("Stopping playback")
                g_stopPlaying = true
            end
        end #IsTrackingUser
    end #for
end

println("Finished playback")


text:
abstract type MultiplicativeStructure{Twisted,I} <: AbstractMatrix{I} end

_istwisted(::MultiplicativeStructure{T}) where {T} = T
_product(ms::MultiplicativeStructure, g, h) = _product(Val(_istwisted(ms)), g, h)

_product(::Val{false}, g, h) = g * h
_product(::Val{true}, g, h) = star(g) * h

struct ProductNotDefined <: Exception
    i::Any
    j::Any
    msg::Any
end

function Base.showerror(io::IO, ex::ProductNotDefined)
    print(io, "Product of elements $(ex.i) and $(ex.j) is not defined on the basis")
    print(io, " or the multiplicative structure could not be completed")
    if isdefined(ex, :msg)
        print(io, ": $(ex.msg)")
    end
    print(io, ".")
end

struct TrivialMStructure{Tw,I,B<:AbstractBasis} <: MultiplicativeStructure{Tw,I}
    basis::B
end

TrivialMStructure{Tw}(basis::AbstractBasis{T,I}) where {Tw,T,I} =
    TrivialMStructure{Tw,I,typeof(basis)}(basis)

basis(mstr::TrivialMStructure) = mstr.basis
Base.size(mstr::TrivialMStructure) = (l = length(basis(mstr)); (l, l))

Base.@propagate_inbounds function Base.getindex(
    mstr::TrivialMStructure,
    i::Integer,
    j::Integer,
)
    @boundscheck checkbounds(mstr, i, j)
    b = basis(mstr)
    g, h = b[i], b[j]
    gh = _product(mstr, g, h)
    gh in b || throw(ProductNotDefined(i, j, "$g · $h = $gh"))
    return b[gh]
end


text:
module MicrostructureNoise

using SpecialFunctions, Distributions

using Statistics

include("microstructure.jl")

end # module


text:
"""
    eliminate!(g::SimpleGraph{Int64}, v::Int64, labels)::Nothing
Connect neighbors of v (in graph g) together if they are not already connected, delete v
and update the labels array. Labels is pop and element v updated.

# Arguments
- `g::SimpleGraph{Int64}` the graph to consider.
- `v::Int64` the vertex to eliminate.
- `labels::Vector{Int64}` labels to keep track of nodes id in original graph before elimination.

# Return
- Nothing.
"""
function eliminate!(g::SimpleGraph{Int64}, v::Int64, labels::Vector{Int64})::Nothing
    ns = neighbors(g, v)
    len = length(ns)
    for i = 1:(len-1)
        for j = (i+1):len
            add_edge!(g, ns[i], ns[j])
        end
    end
    rem_vertex!(g, v)
    labels[v] = labels[end]
    pop!(labels)
    nothing
end

"""
    treewidth_by_elimination!(g::SimpleGraph{Int64}, order::Vector{Int64})
Computes the treewidth by eliminating in order vertices of the graph and keep track of the max degree.

# Arguments
- `g::SimpleGraph{Int64}` the graph to consider.
- `order::Vector{Int64}` the vertex elimination order to perform.

# Return
- `treewidth::Int64` the max degree of vertices during elimination.
"""
function treewidth_by_elimination!(g::SimpleGraph{Int64}, order::Vector{Int64})::Int64
    treewidth = 0
    n = nv(g)
    labels = collect(1:n)
    for rm_origin_id in order
        toremove = findfirst(id -> id == rm_origin_id, labels)
        treewidth = max(treewidth, degree(g, toremove))
        eliminate!(g, toremove, labels)
    end
    return treewidth
end

"""
    tree_order!(graph, nodes)
Computes the optimal elimination order for nodes in a tree. Graph is modified, all its edges are deleted.

# Arguments
- `graph::SimpleGraph{Int64}` the subtree to consider.
- `nodes::Vector{Int64}` indices of nodes in the root greaph.

# Return
- `order::Vector{Int64}` vertex elimination order (indices in the root graph) to get treewidth of 1 in the subtree.
"""
function tree_order!(graph::SimpleGraph{Int64}, nodes::Vector{Int64})::Vector{Int64}
    n = length(nodes)
    eliminated = 1
    order = zeros(Int64, n)
    leafs = filter(node -> length(neighbors(graph, node)) == 1, 1:n)

    # trick to get the last node
    lastnode = sum(nodes)
    while eliminated < n
        parents = unique(map(it -> neighbors(graph, it)[1], leafs))
        for rm in leafs
            order[eliminated] = nodes[rm]
            lastnode -= nodes[rm]
            if ne(graph) != 0
                rem_edge!(graph, rm, neighbors(graph, rm)[1])
            end
            eliminated += 1
        end
        # new leafs are parents with 1 neighbors
        leafs = filter(
            node -> length(neighbors(graph, node)) == 1,
            parents
        )
    end
    # fix last node have no neighbor because all edges are deleted
    if order[n] == 0
        order[n] = lastnode
    end
    return order
end

"""
    graph_from_gr(filename::String)

Read a graph from the provided gr file.
"""
function graph_from_gr(filename::String)::SimpleGraph{Int64}
    lines = readlines(filename)

    # Create a Graph with the correct number of vertices.
    num_vertices, num_edges = parse.(Int, split(lines[1], ' ')[3:end])
    G = SimpleGraph(num_vertices)

    # Add an edge to the graph for every other line in the file.
    for line in lines[2:end]
        src, dst = parse.(Int, split(line, ' '))
        add_edge!(G, src, dst)
    end

    G
end

"""
    square_lattice_graph(n::Int64)
Return a square lattice graph of dimension n. The dimension means
a big square made of n x n squares so (n+1)² nodes in total.

# Example

n = 2 give a graph like this:
    1 ---- 2 ---- 3
    |      |      |
    4 ---- 5 ---- 6
    |      |      |
    7 ---- 8 ---- 9

"""
function square_lattice_graph(n::Int64)::SimpleGraph{Int64}
    nrow = (n + 1)
    nvertices =  nrow * nrow
    g = SimpleGraph{Int64}(nvertices)
    for i in 1:nvertices
        if i % nrow != 0
            add_edge!(g, i, i + 1)
        end
        if trunc(i / nrow) < nrow - 1
            add_edge!(g, i, i + nrow)
        end
    end
    return g
end


text:
# Exceptions
struct IndexOperationException <: Exception
    op::String
    itype::String
end

Base.showerror(io::IO, e::IndexOperationException) =
    print(io, "Failed call `$(e.op)` on `$(e.itype)` index.")

# Abstract types
abstract type AbstractIndex end

"""
    knn_search(index, point, k, keep)

Searches for the `k` nearest neighbors of `point` in data contained in
the `index`. The index may vary from a simple wrapper inside a matrix
to more complex structures such as k-d trees, etc. Only neighbors
present in `keep` are returned.
"""
function knn_search(index::AbstractIndex,
                    point::AbstractVector,
                    k::Integer,
                    keep::AbstractVector;
                    kwargs...)
    throw(IndexOperationException("knn_search", string(typeof(index))))
end

"""
    length(index)

Returns the number of points indexed in `index`.
"""
function length(index::AbstractIndex)
    throw(IndexOperationException("length", string(typeof(index))))
end

### Utility function for indexes
"""
    densify(array)

Transforms sparse arrays into dense ones.
"""
densify(m::AbstractMatrix) = m
densify(m::AbstractSparseMatrix{T,I}) where {T,I} = Matrix{T}(m)
densify(v::AbstractVector) = v
densify(v::AbstractSparseVector{T,I}) where {T,I} = Vector{T}(v)


text:
using PolyaGammaSamplers
using Distributions
using Random
using Statistics

@testset begin
    s1 = PolyaGammaPSWSampler(2, 3.0)
    s2 = PolyaGammaPSWSampler(1, 1.0)
    @test mean(s1) ≈ 0.3017160845482888
    @test mean(s2) ≈ 0.23105857863000487
    @test var(s2) ≈ 0.03444664538852302
end

@testset begin
    t = 0.64
    tol = 1e-3
    rng = MersenneTwister(1)
    @test abs(PolyaGammaSamplers.mass_texpon(0.0, t) - 0.5776972) < tol
    @test abs(PolyaGammaSamplers.mass_texpon(1.0, t) - 0.4605903) < tol
    @test abs(PolyaGammaSamplers.mass_texpon(2.0, t) - 0.2305365) < tol
    x = [PolyaGammaSamplers.rtigauss(rng, 1.0, t) for _ in 1:10000]
    @test abs(mean(x) - .372498) < .005
end

@testset begin
    tol = 5
    n = 10^6

    s = PolyaGammaPSWSampler(1, 1.0)
    @test abs(mean(s) - mean(rand(s, n))) < tol * √(var(s) / n)
    @test abs(var(s) - var(rand(s, n))) < tol * √(2 * var(s)^2 / n)

    s = PolyaGammaPSWSampler(1, 5.0)
    @test abs(mean(s) - mean(rand(s, n))) < tol * √(var(s) / n)
    @test abs(var(s) - var(rand(s, n))) < tol * √(2 * var(s)^2 / n)
end


text:
# This file was modified from part of Julia's Base:
# `base/summarysize.jl`

module MemorySummarySize

using Core: SimpleVector
using Base: unsafe_convert, isbitsunion, unwrap_unionall, isdeprecated, gc_alignment

Base.@kwdef mutable struct FieldResult
    # The unique id for the FieldResultalue itself
    objectid::UInt = 0x0

    size::Int = 0
    type::Type
    is_collection::Bool = false
    skipped_self_reference::Bool = false
    children::Dict{String, FieldResult} = Dict()
    parent::Union{FieldResult,Nothing} = nothing
end
function Base.show(io::IO, fr::FieldResult)
    # Print everything but parent
    print(io,"FieldResult($(fr.objectid), $(fr.size), $(fr.type), $(fr.is_collection), $(fr.children))")
end

mutable struct FrontierNode
    x
    i::Int
    parent_result::FieldResult
end

struct SummarySize
    seen::IdDict{Any,Any}
    fieldresult::FieldResult
    frontier::Vector{FrontierNode}
    exclude::Any
    chargeall::Any
    SummarySize(obj, exclude, chargeall) = new(IdDict(), FieldResult(type=typeof(obj)), [], exclude, chargeall)
end

"""
    Base.summarysize(obj; exclude=Union{...}, chargeall=Union{...}) -> Int

Compute the amount of memory, in bytes, used by all unique objects reachable from the argument.

# Keyword Arguments
- `exclude`: specifies the types of objects to exclude from the traversal.
- `chargeall`: specifies the types of objects to always charge the size of all of their
  fields, even if those fields would normally be excluded.
"""
function summarysize(obj;
                     parentname="obj",
                     exclude = Union{DataType, Core.TypeName, Core.MethodInstance},
                     chargeall = Union{Core.TypeMapEntry, Method})
    @nospecialize obj exclude chargeall
    ss = SummarySize(obj, exclude, chargeall)
    parent = ss.fieldresult
    size::Int = ss(parent, obj)
    while !isempty(ss.frontier)
        # DFS heap traversal of everything without a specialization
        # BFS heap traversal of anything with a specialization
        node = ss.frontier[end]
        x = node.x
        i = node.i
        parent_result = node.parent_result
        val = nothing
        name = parentname
        if isa(x, SimpleVector)
            nf = length(x)
            if isassigned(x, i)
                val = x[i]
                name = "$i"
            end
        elseif isa(x, Array)
            nf = length(x)
            if ccall(:jl_array_isassigned, Cint, (Any, UInt), x, i - 1) != 0
                val = x[i]
                name = "$i"
            end
        else
            nf = nfields(x)
            ft = typeof(x).types
            if !isbitstype(ft[i]) && isdefined(x, i)
                val = getfield(x, i)
                name = "$(fieldname(typeof(x), i))"
            end
        end
        if nf > i
            ss.frontier[end].i = i + 1
        else
            pop!(ss.frontier)
        end
        if val !== nothing && !isa(val, Module) && (!isa(val, ss.exclude) || isa(x, ss.chargeall))
            fieldresult = get!(parent_result.children, name, FieldResult(type=typeof(val),parent=parent_result))
            valsize = ss(fieldresult, val)::Int
            _finish_fieldresult(fieldresult, val, valsize)
            p = fieldresult.parent
            while p !== nothing
                p.size += valsize
                p = p.parent
            end
            size += valsize
        end
    end
    _finish_fieldresult(ss.fieldresult, obj, size)
    return ss.fieldresult
end
function _finish_fieldresult(fieldresult, val, size)
    fieldresult.objectid = Base.objectid(val)
    fieldresult.size = size
    # Mark Collection types not handled below
    if val isa Tuple
        fieldresult.is_collection = true
    end
end

(ss::SummarySize)(fieldresult, @nospecialize obj) = _summarysize(ss, fieldresult, obj)
# define the general case separately to make sure it is not specialized for every type
@noinline function _summarysize(ss::SummarySize, fieldresult, @nospecialize obj)
    isdefined(typeof(obj), :instance) && return 0
    # NOTE: this attempts to discover multiple copies of the same immutable value,
    # and so is somewhat approximate.
    key = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), obj)
    if handle_seen(fieldresult, ss, key) return 0 end
    if nfields(obj) > 0
        push!(ss.frontier, FrontierNode(obj, 1, fieldresult))
    end
    if isa(obj, UnionAll) || isa(obj, Union)
        # black-list of items that don't have a Core.sizeof
        sz = 2 * sizeof(Int)
    else
        sz = Core.sizeof(obj)
    end
    if sz == 0
        # 0-field mutable structs are not unique
        return gc_alignment(0)
    end
    return sz
end
function handle_seen(fieldresult, ss, obj)
    if haskey(ss.seen, obj)
        fieldresult.skipped_self_reference = true
        return true
    else
        ss.seen[obj] = true
        return false
    end
end

(::SummarySize)(fieldresult, obj::Symbol) = 0
(::SummarySize)(fieldresult, obj::SummarySize) = 0

function (ss::SummarySize)(f, obj::String)
    key = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), obj)
    if handle_seen(f, ss, obj) return 0 end
    return Core.sizeof(Int) + Core.sizeof(obj)
end

function (ss::SummarySize)(f, obj::DataType)
    key = pointer_from_objref(obj)
    if handle_seen(f, ss, obj) return 0 end
    size::Int = 7 * Core.sizeof(Int) + 6 * Core.sizeof(Int32)
    size += 4 * nfields(obj) + ifelse(Sys.WORD_SIZE == 64, 4, 0)
    size += ss(f, obj.parameters)::Int
    if isdefined(obj, :types)
        size += ss(f, obj.types)::Int
    end
    return size
end

function (ss::SummarySize)(fieldresult, obj::Core.TypeName)
    key = pointer_from_objref(obj)
    if handle_seen(fieldresult, ss, obj) return 0 end
    return Core.sizeof(obj) + (isdefined(obj, :mt) ? ss(fieldresult, obj.mt) : 0)
end

function (ss::SummarySize)(fieldresult, obj::Array)
    fieldresult.is_collection = true
    if handle_seen(fieldresult, ss, obj) return 0 end
    headersize = 4*sizeof(Int) + 8 + max(0, ndims(obj)-2)*sizeof(Int)
    size::Int = headersize
    datakey = unsafe_convert(Ptr{Cvoid}, obj)
    if !haskey(ss.seen, datakey)
        ss.seen[datakey] = true
        dsize = Core.sizeof(obj)
        if isbitsunion(eltype(obj))
            # add 1 union selector byte for each element
            dsize += length(obj)
        end
        size += dsize
        if !isempty(obj) && !Base.allocatedinline(eltype(obj))
            push!(ss.frontier, FrontierNode(obj, 1, fieldresult))
        end
    end
    return size
end

function (ss::SummarySize)(fieldresult, obj::SimpleVector)
    fieldresult.is_collection = true
    key = pointer_from_objref(obj)
    if handle_seen(fieldresult, ss, obj) return 0 end
    size::Int = Core.sizeof(obj)
    if !isempty(obj)
        push!(ss.frontier, FrontierNode(obj, 1, fieldresult))
    end
    return size
end

function (ss::SummarySize)(fieldresult, obj::Module)
    if handle_seen(fieldresult, ss, obj) return 0 end
    size::Int = Core.sizeof(obj)
    for binding in names(obj, all = true)
        if isdefined(obj, binding) && !isdeprecated(obj, binding)
            value = getfield(obj, binding)
            if !isa(value, Module) || parentmodule(value) === obj
                size += ss(fieldresult, value)::Int
                if isa(value, UnionAll)
                    value = unwrap_unionall(value)
                end
                if isa(value, DataType) && value.name.module === obj && value.name.name === binding
                    # charge a TypeName to its module (but not to the type)
                    size += ss(fieldresult, value.name)::Int
                end
            end
        end
    end
    return size
end

function (ss::SummarySize)(fieldresult, obj::Task)
    if handle_seen(fieldresult, ss, obj) return 0 end
    size::Int = Core.sizeof(obj)
    if isdefined(obj, :code)
        size += ss(fieldresult, obj.code)::Int
    end
    size += ss(fieldresult, obj.storage)::Int
    size += ss(fieldresult, obj.backtrace)::Int
    size += ss(fieldresult, obj.donenotify)::Int
    size += ss(fieldresult, obj.exception)::Int
    size += ss(fieldresult, obj.result)::Int
    # TODO: add stack size, and possibly traverse stack roots
    return size
end

end # module


text:
#=
A type to solve the career / job choice model due to Derek Neal.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-05

References
----------

http://quant-econ.net/jl/career.html

[Neal1999] Neal, D. (1999). The Complexity of Job Mobility among Young Men,
Journal of Labor Economics, 17(2), 237-261.
=#

"""
Career/job choice model fo Derek Neal (1999)

##### Fields

- `beta::Real` : Discount factor in (0, 1)
- `N::Int` : Number of possible realizations of both epsilon and theta
- `B::Real` : upper bound for both epsilon and theta
- `theta::AbstractVector` : A grid of values on [0, B]
- `epsilon::AbstractVector` : A grid of values on [0, B]
- `F_probs::AbstractVector` : The pdf of each value associated with of F
- `G_probs::AbstractVector` : The pdf of each value associated with of G
- `F_mean::Real` : The mean of the distribution F
- `G_mean::Real` : The mean of the distribution G

"""
type CareerWorkerProblem
    beta::Real
    N::Int
    B::Real
    theta::AbstractVector
    epsilon::AbstractVector
    F_probs::AbstractVector
    G_probs::AbstractVector
    F_mean::Real
    G_mean::Real
end

"""
Constructor with default values for `CareerWorkerProblem`

##### Arguments

- `beta::Real(0.95)` : Discount factor in (0, 1)
- `B::Real(5.0)` : upper bound for both epsilon and theta
- `N::Real(50)` : Number of possible realizations of both epsilon and theta
- `F_a::Real(1), F_b::Real(1)` : Parameters of the distribution F
- `G_a::Real(1), G_b::Real(1)` : Parameters of the distribution F

##### Notes

There is also a version of this function that accepts keyword arguments for
each parameter
"""
# use key word argument
function CareerWorkerProblem(;beta::Real=0.95, B::Real=5.0, N::Real=50,
                             F_a::Real=1, F_b::Real=1, G_a::Real=1,
                             G_b::Real=1)
    theta = linspace(0, B, N)
    epsilon = copy(theta)
    F_probs::Vector{Float64} = pdf(BetaBinomial(N-1, F_a, F_b))
    G_probs::Vector{Float64} = pdf(BetaBinomial(N-1, G_a, G_b))
    F_mean = sum(theta .* F_probs)
    G_mean = sum(epsilon .* G_probs)
    CareerWorkerProblem(beta, N, B, theta, epsilon, F_probs, G_probs,
                        F_mean, G_mean)
end

"""
Apply the Bellman operator for a given model and initial value.

##### Arguments

- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output
- `;ret_policy::Bool(false)`: Toggles return of value or policy functions

##### Returns

None, `out` is updated in place. If `ret_policy == true` out is filled with the
policy function, otherwise the value function is stored in `out`.

"""
function update_bellman!(cp::CareerWorkerProblem, v::Array, out::Array;
                           ret_policy=false)
    # new life. This is a function of the distribution parameters and is
    # always constant. No need to recompute it in the loop
    v3 = (cp.G_mean + cp.F_mean + cp.beta .*
          cp.F_probs' * v * cp.G_probs)[1]  # don't need 1 element array

    for j=1:cp.N
        for i=1:cp.N
            # stay put
            v1 = cp.theta[i] + cp.epsilon[j] + cp.beta * v[i, j]

            # new job
            v2 = (cp.theta[i] .+ cp.G_mean .+ cp.beta .*
                  v[i, :]'*cp.G_probs)[1]  # don't need a single element array

            if ret_policy
                if v1 > max(v2, v3)
                    action = 1
                elseif v2 > max(v1, v3)
                    action = 2
                else
                    action = 3
                end
                out[i, j] = action
            else
                out[i, j] = max(v1, v2, v3)
            end
        end
    end
end

function update_bellman(cp::CareerWorkerProblem, v::Array; ret_policy=false)
    out = similar(v)
    update_bellman!(cp, v, out, ret_policy=ret_policy)
    return out
end

"""
Extract the greedy policy (policy function) of the model.

##### Arguments

- `cp::CareerWorkerProblem` : Instance of `CareerWorkerProblem`
- `v::Matrix`: Current guess for the value function
- `out::Matrix` : Storage for output

##### Returns

None, `out` is updated in place to hold the policy function

"""
function get_greedy!(cp::CareerWorkerProblem, v::Array, out::Array)
    update_bellman!(cp, v, out, ret_policy=true)
end

function get_greedy(cp::CareerWorkerProblem, v::Array)
    update_bellman(cp, v, ret_policy=true)
end


text:
using KungIChi
using Test

function multiply_by_two(y)
    @test true
    return 2y
end

@schema struct Foo
    x::Int = 1
    y::Int

    @validate function (x)
        @test true
        @assert x < 10
        return x
    end

    @validate multiply_by_two(y)

    @check multiple_of(y, 1.0)
    @check multiple_of(x, 2.0)
    @check function (x)
        x > 0
    end
end

@testset "@validate/@check" begin
    @test Foo(2, 2).y == 4
    @test_throws ValidationError Foo(1, 2)
    @test_throws ValidationError Foo(-2, 2)
end


text:
using Distributions
using Plots

unicodeplots()

Σ₁₁ = 1.5
Σ₂₂ = 3.0
Σ₁₂ = Σ₂₁ = 0.7

μ₁ = 3.0
μ₂ = -1.0

μ = [
    μ₁
    μ₂
]

Σ = [
    Σ₁₁ Σ₁₂
    Σ₂₁ Σ₂₂
]

function g(x, μ, σ)
    σ² = σ^2
    return exp(-0.5 * (x - μ)^2 / σ²) / √(2π * σ²)
end

N = 10000
d = MvNormal(μ, Σ)

data = rand(d, N)
xs, ys = data[1, :], data[2, :]

p1 = plot(size = (300, 300))
histogram!(p1, xs, norm = true)
plot!(p1, x -> g(x, μ₁, √Σ₁₁))

p2 = plot(size = (300, 300))
histogram!(p2, ys, norm = true)
plot!(p2, x -> g(x, μ₂, √Σ₂₂))

plot(p1, p2)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Documenter, ElectronDisplay

makedocs(
	modules=[ElectronDisplay],
	sitename="ElectronDisplay.jl",
	analytics="UA-132838790-1",
	pages=[
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo="github.com/queryverse/ElectronDisplay.jl.git"
)


text:
using Combinatorics

const HOLES = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H']
const PEGS = [1, 2, 3, 4, 5, 6, 7, 8]
const EDGES = [('A', 'C'), ('A', 'D'), ('A', 'E'),
               ('B', 'D'), ('B', 'E'), ('B', 'F'),
               ('C', 'G'), ('C', 'D'), ('D', 'G'),
               ('D', 'E'), ('D', 'H'), ('E', 'F'),
               ('E', 'G'), ('E', 'H'), ('F', 'H')]

goodperm(p) = all(e->abs(p[e[1]-'A'+1] - p[e[2]-'A'+1]) > 1, EDGES)

goodplacements() = [p for p in permutations(PEGS) if goodperm(p)]

const BOARD = raw"""
        A   B
       /|\ /|\
      / | X | \
     /  |/ \|  \
    C - D - E - F
     \  |\ /|  /
      \ | X | /
       \|/ \|/
        G   H
"""

function printsolutions()
    solutions = goodplacements()
    println("Found $(length(solutions)) solutions.")
    for soln in solutions
        board = BOARD
        for (i, n) in enumerate(soln)
            board = replace(board, string('A' + i - 1) => string(n))
        end
        println(board); exit(1) # remove this exit for all solutions
    end
end

printsolutions()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# - Parse Json_dict to Struct

function json2ps_struct(data::Dict{String,Any})
    """
    Takes a PowerSystems dictionary and return an array of PowerSystems struct for Bus, Generator, Branch and load
    """
    if haskey(data, "bus")
        Buses = bus_json_parse(data["bus"])
    else
        @error "Key error : key 'bus' not found in PowerSystems dictionary, this will result in an empty Bus array"
        Buses =[]
    end
    if haskey(data, "gen")
        Generators, Storage = gen_json_parser(data["gen"])
    else
        @error "Key error : key 'gen' not found in PowerSystems dictionary, this will result in an empty Generators and Storage array"
        Generators =[]
        Storage = []
    end
    if haskey(data, "branch")
        Branches = branch_json_parser(data["branch"])
    else
        @warn "Key error : key 'branch' not found in PowerSystems dictionary, this will result in an empty Branches array"
        Branches =[]
    end
    if haskey(data, "load")
        Loads = load_json_parser(data["load"])
    else
        @error "Key error : key 'load'  not found in PowerSystems dictionary, this will result in an empty Loads array"
        Loads =[]
    end
    return Buses, Generators, Storage, Branches, Loads
end

function bus_json_parse(dict::Dict{String,Any})
    Buses = Array{Bus}(undef, 0)
    for (bus_key,bus_dict) in dict
        push!(Buses,Bus(bus_dict["number"],
                                    bus_dict["name"],
                                    bus_dict["bustype"],
                                    bus_dict["angle"],
                                    bus_dict["voltage"],
                                    (min =bus_dict["voltagelimits"]["min"],max=bus_dict["voltagelimits"]["max"]),
                                    bus_dict["basevoltage"]
                                    ))
    end
    return Buses
end

function gen_json_parser(dict::Dict{String,Any})
    Generators = Array{Generator}(undef, 0)
    Storage_gen = Array{Storage}(undef, 0)
    for (gen_type_key,gen_type_dict) in dict
        if gen_type_key =="Thermal"
            for (thermal_key,thermal_dict) in gen_type_dict
                push!(Generators,ThermalStandard(thermal_dict["name"],
                                                            thermal_dict["available"],
                                                            Bus(thermal_dict["bus"]["number"],
                                                                            thermal_dict["bus"]["name"],
                                                                            thermal_dict["bus"]["bustype"],
                                                                            thermal_dict["bus"]["angle"],
                                                                            thermal_dict["bus"]["voltage"],
                                                                            (min =thermal_dict["bus"]["voltagelimits"]["min"],max=thermal_dict["bus"]["voltagelimits"]["max"]),
                                                                            thermal_dict["bus"]["basevoltage"] ),
                                                            TechThermal(thermal_dict["tech"]["activepower"],
                                                                        (min =thermal_dict["tech"]["activepowerlimits"]["min"],max =thermal_dict["tech"]["activepowerlimits"]["max"]),
                                                                        thermal_dict["tech"]["reactivepower"],
                                                                        (min =thermal_dict["tech"]["reactivepowerlimits"]["min"],max =thermal_dict["tech"]["reactivepowerlimits"]["min"]),
                                                                        (up=thermal_dict["tech"]["ramplimits"]["up"],down=thermal_dict["tech"]["ramplimits"]["down"]),
                                                                        (up=thermal_dict["tech"]["timelimits"]["up"],down=thermal_dict["tech"]["timelimits"]["down"])),
                                                            ThreePartCost(json_var_cost(thermal_dict["econ"]["variable"]),
                                                                          thermal_dict["econ"]["fixedcost"],
                                                                          thermal_dict["econ"]["startup"],
                                                                          thermal_dict["econ"]["shutdn"])
                            ))
            end
        elseif gen_type_key =="Hydro"
            for (hydro_key,hydro_dict) in gen_type_dict
                push!(Generators,HydroDispatch(hydro_dict["name"],
                                                            hydro_dict["available"],
                                                            Bus(hydro_dict["bus"]["number"],
                                                                            hydro_dict["bus"]["name"],
                                                                            hydro_dict["bus"]["bustype"],
                                                                            hydro_dict["bus"]["angle"],
                                                                            hydro_dict["bus"]["voltage"],
                                                                            (min =hydro_dict["bus"]["voltagelimits"]["min"],max=hydro_dict["bus"]["voltagelimits"]["max"]),
                                                                            hydro_dict["bus"]["basevoltage"] ),
                                                            TechHydro(  hydro_dict["tech"]["rating"],
                                                                        hydro_dict["tech"]["activepower"],
                                                                        (min =hydro_dict["tech"]["activepowerlimits"]["min"],max =hydro_dict["tech"]["activepowerlimits"]["max"]),
                                                                        hydro_dict["tech"]["reactivepower"],
                                                                        (min =hydro_dict["tech"]["reactivepowerlimits"]["min"],max =hydro_dict["tech"]["reactivepowerlimits"]["max"]),
                                                                        (up=hydro_dict["tech"]["ramplimits"]["up"],down=hydro_dict["tech"]["ramplimits"]["down"]),
                                                                        (up=hydro_dict["tech"]["timelimits"]["up"],down=hydro_dict["tech"]["timelimits"]["down"])),
                                                            hydro_dict["econ"]["curtailcost"],
                                                            dict_to_timearray(hydro_dict["scalingfactor"])
                            ))
            end
        elseif gen_type_key =="Renewable"
            for (ren_key,ren_dict) in  gen_type_dict
                if ren_key == "PV"
                    for (pv_key,pv_dict) in ren_dict
                        push!(Generators,RenewableDispatch(pv_dict["name"],
                                                                    pv_dict["available"],
                                                                    Bus(pv_dict["bus"]["number"],
                                                                            pv_dict["bus"]["name"],
                                                                            pv_dict["bus"]["bustype"],
                                                                            pv_dict["bus"]["angle"],
                                                                            pv_dict["bus"]["voltage"],
                                                                            (min =pv_dict["bus"]["voltagelimits"]["min"],max=pv_dict["bus"]["voltagelimits"]["max"]),
                                                                            pv_dict["bus"]["basevoltage"] ),
                                                                    pv_dict["tech"]["rating"],
                                                                    TwoPartCost(pv_dict["econ"]["curtailcost"],
                                                                                pv_dict["econ"]["interruptioncost"]),
                                                                    dict_to_timearray(pv_dict["scalingfactor"])
                                    ))
                    end
                elseif ren_key == "RTPV"
                    for (rtpv_key,rtpv_dict) in ren_dict
                        push!(Generators,RenewableFix(rtpv_dict["name"],
                                                                    rtpv_dict["available"],
                                                                    Bus(rtpv_dict["bus"]["number"],
                                                                            rtpv_dict["bus"]["name"],
                                                                            rtpv_dict["bus"]["bustype"],
                                                                            rtpv_dict["bus"]["angle"],
                                                                            rtpv_dict["bus"]["voltage"],
                                                                            (min =rtpv_dict["bus"]["voltagelimits"]["min"],max=rtpv_dict["bus"]["voltagelimits"]["max"]),
                                                                            rtpv_dict["bus"]["basevoltage"] ),
                                                                    rtpv_dict["tech"]["rating"],
                                                                    dict_to_timearray(rtpv_dict["scalingfactor"])
                                    ))
                    end
                elseif ren_key == "WIND"
                    for (wind_key,wind_dict) in ren_dict
                        push!(Generators,RenewableDispatch(wind_dict["name"],
                                                                    wind_dict["available"],
                                                                    Bus(wind_dict["bus"]["number"],
                                                                            wind_dict["bus"]["name"],
                                                                            wind_dict["bus"]["bustype"],
                                                                            wind_dict["bus"]["angle"],
                                                                            wind_dict["bus"]["voltage"],
                                                                            (min =wind_dict["bus"]["voltagelimits"]["min"],max=wind_dict["bus"]["voltagelimits"]["max"]),
                                                                            wind_dict["bus"]["basevoltage"] ),
                                                                    wind_dict["tech"]["rating"],
                                                                    TwoPartCost(wind_dict["econ"]["curtailcost"],
                                                                                wind_dict["econ"]["interruptioncost"]),
                                                                    dict_to_timearray(wind_dict["scalingfactor"])
                                    ))
                    end
                end
            end
        elseif gen_type_key =="Storage"
            for (storage_key,storage_dict) in  gen_type_dict
                push!(Storage_gen,GenericBattery(storage_dict["name"],
                                                            storage_dict["available"],
                                                            Bus(storage_dict["bus"]["number"],
                                                                            storage_dict["bus"]["name"],
                                                                            storage_dict["bus"]["bustype"],
                                                                            storage_dict["bus"]["angle"],
                                                                            storage_dict["bus"]["voltage"],
                                                                            (min =storage_dict["bus"]["voltagelimits"]["min"],max=storage_dict["bus"]["voltagelimits"]["max"]),
                                                                            storage_dict["bus"]["basevoltage"] ),
                                                            storage_dict["energy"],
                                                            (min=storage_dict["capacity"]["min"],max=storage_dict["capacity"]["max"]),
                                                            storage_dict["activepower"],
                                                            storage_dict["inputactivepowerlimit"],
                                                            storage_dict["outputactivepowerlimit"],
                                                            (in =storage_dict["efficiency"]["in"],out=storage_dict["efficiency"]["out"]),
                                                            storage_dict["reactivepower"],
                                                            (min=storage_dict["reactivepowerlimits"]["min"],max=storage_dict["reactivepowerlimits"]["max"])
                            ))
            end
        end
    end
    return Generators, Storage_gen
end

function branch_json_parser(dict)
    Branches = Array{Branch}(undef, 0)
    for (branch_key,branch_dict) in dict
        if branch_key == "Transformers"
            for (trans_key,trans_dict) in branch_dict
                bus_f =Bus(trans_dict["arc"]["from"]["number"],
                                        trans_dict["arc"]["from"]["name"],
                                        trans_dict["arc"]["from"]["bustype"],
                                        trans_dict["arc"]["from"]["angle"],
                                        trans_dict["arc"]["from"]["voltage"],
                                        (min =trans_dict["arc"]["from"]["voltagelimits"]["min"],max=trans_dict["arc"]["from"]["voltagelimits"]["max"]),
                                        trans_dict["arc"]["from"]["basevoltage"] )
                bus_t =Bus(trans_dict["arc"]["to"]["number"],
                                        trans_dict["arc"]["to"]["name"],
                                        trans_dict["arc"]["to"]["bustype"],
                                        trans_dict["arc"]["to"]["angle"],
                                        trans_dict["arc"]["to"]["voltage"],
                                        (min =trans_dict["arc"]["to"]["voltagelimits"]["min"],max=trans_dict["arc"]["to"]["voltagelimits"]["max"]),
                                        trans_dict["arc"]["to"]["basevoltage"] )
                if trans_dict["tap"] ==1.0
                    push!(Branches,Transformer2W(trans_dict["name"],
                                                trans_dict["available"],
                                                (from = bus_f,to = bus_t),
                                                trans_dict["r"],
                                                trans_dict["x"],
                                                trans_dict["primaryshunt"],
                                                trans_dict["rate"]
                                                ))
                elseif trans_dict["tap"] !=1.0
                    push!(Branches,TapTransformer(trans_dict["name"],
                                                trans_dict["available"],
                                                (from = bus_f,to = bus_t),
                                                trans_dict["r"],
                                                trans_dict["x"],
                                                trans_dict["primaryshunt"],
                                                trans_dict["tap"],
                                                trans_dict["rate"]
                                                ))
                end
            end
        elseif branch_key == "Lines"
            for (line_key,line_dict) in branch_dict
                bus_t =Bus(line_dict["arc"]["to"]["number"],
                                        line_dict["arc"]["to"]["name"],
                                        line_dict["arc"]["to"]["bustype"],
                                        line_dict["arc"]["to"]["angle"],
                                        line_dict["arc"]["to"]["voltage"],
                                        (min =line_dict["arc"]["to"]["voltagelimits"]["min"],max=line_dict["arc"]["to"]["voltagelimits"]["max"]),
                                        line_dict["arc"]["to"]["basevoltage"] )
                bus_f =Bus(line_dict["arc"]["from"]["number"],
                                        line_dict["arc"]["from"]["name"],
                                        line_dict["arc"]["from"]["bustype"],
                                        line_dict["arc"]["from"]["angle"],
                                        line_dict["arc"]["from"]["voltage"],
                                        (min =line_dict["arc"]["from"]["voltagelimits"]["min"],max=line_dict["arc"]["from"]["voltagelimits"]["max"]),
                                        line_dict["arc"]["from"]["basevoltage"] )
                push!(Branches,Line(line_dict["name"],
                                    line_dict["available"],
                                    (from = bus_f, to = bus_t),
                                    line_dict["r"],
                                    line_dict["x"],
                                    (from =line_dict["b"]["from"],to=line_dict["b"]["to"]),
                                    line_dict["rate"],
                                    line_dict["anglelimits"]
                                    ))
            end
        end
    end
    return Branches
end

function load_json_parser(dict)
    Loads = Array{ElectricLoad}(undef, 0)
    for (load_key,load_dict) in dict
        push!(Loads, PowerLoad(load_dict["name"],
                load_dict["available"],
                Bus(load_dict["bus"]["number"],
                                load_dict["bus"]["name"],
                                load_dict["bus"]["bustype"],
                                load_dict["bus"]["angle"],
                                load_dict["bus"]["voltage"],
                                (min =load_dict["bus"]["voltagelimits"]["min"],max=load_dict["bus"]["voltagelimits"]["max"]),
                                load_dict["bus"]["basevoltage"] ),
                load_dict["model"],
                load_dict["maxactivepower"],
                load_dict["maxreactivepower"],
                dict_to_timearray(load_dict["scalingfactor"])
                ))
    end
    return Loads
end

# Write dict to json file
function dict_to_json(dict,filename)
    stringdata =JSON2.write(dict)
    open("$filename.json", "w") do f
        write(f, stringdata)
    end
end

# Read json save file as a dict
function json_parser(filename)
    data =Dict{String,Any}()
    if isfile("../data/CDM/RTS/JSON/RTS-GMLC_Test_Case.json")
        temp =Dict{}()
        open("../data/CDM/RTS/JSON/RTS-GMLC_Test_Case.json", "r") do f
        global temp
        dicttxt = readstring(f)  # file information to string
        temp = JSON2.read(dicttxt, Dict{Any,Array{Dict}})  # parse and transform data
        data = temp
        end
    else
        @warn "JSON file doesn't exist"
    end
    Buses = bus_json_parse(data["Buses"])
    Generators,Storages = gen_json_parser(data["Generators"])
    Branches = branch_json_parser(data["Branches"])
    Loads = load_json_parser(data["Loads"]);
    return Buses,Generators,Branches,Loads
end

# Read back array of Variable cost from json dict
function json_var_cost(array)
    var_cost = [float(array[i][1]) for i in 1:length(array)]
    return var_cost
end

# Read back timeseries time array from json dict
function dict_to_timearray(dict)
    ts_string =[dict["timestamp"][i] for i in 1:length(dict["timestamp"])]
    ts_values =[dict["values"][i] for i in 1:length(dict["values"])]
    ts_raw =TimeSeries.TimeArray(Dates.DateTime(ts_string),ts_values)
    return ts_raw
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Options <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Options(; kwargs...)
        obj = new(meta(Options), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Options
const __meta_Options = Ref{ProtoMeta}()
function meta(::Type{Options})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Options)
            __meta_Options[] = target = ProtoMeta(Options)
            allflds = Pair{Symbol,Union{Type,String}}[:infer => Bool, :trace_inference => Bool, :explain => Bool, :parallel => Bool, :prefetch_size => Int32, :prefetch => Bool, :session_idle_timeout_millis => Int32, :schema_lock_acquire_timeout_millis => Int32, :read_any_replica => Bool]
            oneofs = Int[1,2,3,4,5,6,7,8,9]
            oneof_names = Symbol[Symbol("infer_opt"),Symbol("trace_inference_opt"),Symbol("explain_opt"),Symbol("parallel_opt"),Symbol("prefetch_size_opt"),Symbol("prefetch_opt"),Symbol("session_idle_timeout_opt"),Symbol("schema_lock_acquire_timeout_opt"),Symbol("read_any_replica_opt")]
            meta(target, Options, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_Options[]
    end
end
function Base.getproperty(obj::Options, name::Symbol)
    if name === :infer
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :trace_inference
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :explain
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :parallel
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :prefetch_size
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :prefetch
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :session_idle_timeout_millis
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :schema_lock_acquire_timeout_millis
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :read_any_replica
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

export Options


text: abstract type Drift{R,W} <: CellRule{R,W} end


text:
using OrdinaryDiffEq, DiffEqBase, Base.Test
srand(100)
setprecision(400)

f = (t,u) -> (2u)
analytic = (t,u0) -> u0*exp(t)
prob = ODETestProblem(f,1/2,analytic)

sol3 =solve(prob,RK4(),dt=1/2^(6),abstol=1,reltol=0)

prob = ODETestProblem(f,BigInt(1)//BigInt(2),analytic,(BigInt(0)//BigInt(1),BigInt(1)//BigInt(1)))

integrator = init(prob,RK4(),dt=BigInt(1)//BigInt(2)^(6),abstol=1,reltol=0)

sol =solve(prob,RK4(),dt=BigInt(1)//BigInt(2)^(6),abstol=1,reltol=0)
sol2 =solve(prob,RK4(),dt=BigInt(1)/BigInt(2)^(6),abstol=1,reltol=0)

sol.u
sol2.u
sol3.u
@test eltype(sol.u) == Rational{BigInt}
@test eltype(sol2.u) == Rational{BigInt}
@test eltype(sol3.u) == Float64

sol4 =solve(prob,DP5(),dt=BigInt(1)//BigInt(2)^(3),adaptive=false)

@test eltype(sol4.u) == Rational{BigInt}

tabalg = ExplicitRK(tableau = constructDormandPrince(Rational{BigInt}))

integrator = init(prob,tabalg,dt=BigInt(1)//BigInt(2)^(3),abstol=1,reltol=0,adaptive=false)
sol5 =solve(prob,tabalg,dt=BigInt(1)//BigInt(2)^(3),abstol=1,reltol=0,adaptive=false)

@test typeof(sol5.u[end]) == Rational{BigInt}

prob = ODETestProblem(f,1/2+1/2im,analytic)


text:
using Test, ImageFeatures, Images
import ImageFeatures.trilinear_interpolate!

@testset "HOG Feature" begin
    img = rand(128, 64)
    @test length(create_descriptor(img, HOG())) == 3780

    img = rand(RGB, 128, 64)
    @test length(create_descriptor(img, HOG())) == 3780

    img = zeros(128, 64)
    img[64, :] .= 1
    @test length(create_descriptor(img, HOG())) == 3780

    #tests for function trilinear_interpolate!
    rows = 12
    cols = 12
    cell_size = 4
    cell_rows::Int = rows/cell_size
    cell_cols::Int = cols/cell_size
    orientations = 9

    hist = zeros(9, 3, 3)
    trilinear_interpolate!(hist, 1, 170, orientations, CartesianIndex(1, 1), cell_size, cell_rows, cell_cols, rows, cols)
    @test hist[1, 1, 1] == hist[9, 1, 1] == 0.5

    hist = zeros(9, 3, 3)
    trilinear_interpolate!(hist, 1, 150, orientations, CartesianIndex(1, 4), cell_size, cell_rows, cell_cols, rows, cols)
    @test hist[8, 1, 1] == hist[9, 1, 1] == hist[8, 1, 2] == hist[9, 1, 2] == 0.25

    hist = zeros(9, 3, 3)
    trilinear_interpolate!(hist, 1, 145, orientations, CartesianIndex(4, 4), cell_size, cell_rows, cell_cols, rows, cols)
    @test hist[8, 1, 1] == hist[8, 1, 2] == hist[8, 2, 1] == hist[8, 2, 2] == 0.1875
    @test hist[9, 1, 1] == hist[9, 1, 2] == hist[9, 2, 1] == hist[9, 2, 2] == 0.0625
end


text:
using Random

Base.rand(X::Interval{T}) where {T <: AbstractFloat} =
    X.lo + rand(T) * (X.hi - X.lo)

Random.gentype(::Type{Interval{T}}) where T = T
Random.rand(rng::AbstractRNG, X::Random.SamplerTrivial{Interval{T}}) where T =
    rand(X[])

# Base.eltype(::Type{IntervalBox{2,Float64}}) = SArray{Tuple{2},Float64,1,2}
# Random.rand(rng::AbstractRNG,X::Random.SamplerTrivial{IntervalBox{2,Float64}}) = rand.(X[])


text:
"""
和费米面相关的功能
"""
module Surface

using ..Fermi
using ..Basics

export @onsurface
export const_energy_line, const_energy_line_in_patches
export const_energy_triangle

"""
获取等能面
"""
function const_energy_line(
    ltris::Vector{T},
    ladjs::Vector{Tuple{Union{Missing, T}, Union{Missing, T}, Union{Missing, T}}},
    eng::Float64,
    disp::Function) where T <: Basics.AbstractTriangle
    #计算每个位置上面的能量
    eng_dict = Dict()
    for tri in ltris
        #能量这里减去需要寻找的能量大小
        eng_dict[tri] = disp(tri.center.x, tri.center.y) - eng
    end
    #
    edges::Vector{Basics.Segment} = []
    for (idx, tri) in enumerate(ltris)
        #寻找从负数到正数的边界
        #跳过所有的正数
        if eng_dict[tri] > 0.
            continue
        end
        #查看相邻的几个有没有大于零的
        for (eidx, adj) in enumerate(ladjs[idx])
            #已经到头的时候是不继续添加的
            if ismissing(adj)
                continue
            end
            if eng_dict[adj] > 0
                push!(edges, tri.edges[eidx])
            end
        end
    end
    return edges
end

"""
获取费米面，但是要分好属于哪一组
"""
function const_energy_line_in_patches(
    ltris::Vector{T},
    ladjs::Vector{Tuple{Union{Missing, T}, Union{Missing, T}, Union{Missing, T}}},
    lpats::Vector{Int64},
    eng::Float64,
    disp::Function) where T <: Basics.AbstractTriangle
    eng_dict = Dict()
    for tri in ltris
        #能量这里减去需要寻找的能量大小
        eng_dict[tri] = disp(tri.center.x, tri.center.y) - eng
    end
    #
    edges::Vector{Basics.Segment} = []
    epidx::Vector{Int64} = []
    for (idx, tri) in enumerate(ltris)
        #寻找从负数到正数的边界
        if eng_dict[tri] > 0.
            continue
        end
        #查看相邻的几个有没有大于零的
        for (eidx, adj) in enumerate(ladjs[idx])
            #已经到头的时候是不继续添加的
            if ismissing(adj)
                continue
            end
            if eng_dict[adj] > 0
                push!(edges, tri.edges[eidx])
                #如果大于0和小于0的不在一个patch，会产生一些混淆。但是没明显的影响。
                push!(epidx, lpats[idx])
            end
        end
    end
    return edges, epidx
end

"""
判断小三角是否穿过费米面
"""
macro onsurface(tri, disp, eng)
    return esc(quote
        ver1 = ($tri).vertex[1]
        sgn1 = sign($disp(ver1.x, ver1.y)-$eng)
        ver2 = ($tri).vertex[2]
        sgn2 = sign($disp(ver2.x, ver2.y)-$eng)
        ver3 = ($tri).vertex[3]
        sgn3 = sign($disp(ver3.x, ver3.y)-$eng)
        sgn1 != sgn2 || sgn1 != sgn3
    end)
end

"""
穿过费米面的小三角
"""
function const_energy_triangle(
    ltris::Vector{T},
    eng::Float64,
    disp::Function) where T <: Basics.AbstractTriangle
    #
    #如果一个小三角型的三个顶点的能量符号不同，那就是穿过了费米面
    #
    edges::Vector{T} = []
    for tri in ltris
        if (@onsurface tri disp eng)
            push!(edges, tri)
        end
    end
    return edges
end

end # end module


text:
# Sane behavior when run from the REPL
source_dir = typeof(Base.source_dir()) == Void ? joinpath(Pkg.dir("Crispulator"), "src") : Base.source_dir()
"""
List the custom simulation scripts that are available
"""
function ls()
    files = readdir(joinpath(source_dir, "exps"))
    filter(file -> file != "common.jl", files)
end

"""
bootstrap_exp(analysis_file::ASCIIString, output_dir::ASCIIString, proc_count::Int, debug::Bool)

Run custom simulation using one of the scripts in the `exps/` folder
"""
function bootstrap_exp(analysis_file::Compat.String,
                       output_dir::Compat.String,
                       proc_count::Int,
                       debug::Bool)

    # add additional threads
    addprocs(proc_count)
    println("Using $(nprocs()) threads")

    # load simulation code on all cores
    include(joinpath(source_dir, "simulation", "load.jl"))
    print("Loading analysis code...")
    analysis_full_path = joinpath(source_dir, "exps", analysis_file)

    (!isfile(analysis_full_path)) && error("Please provide a valid analysis file to run")
    (!isdir(dirname(abspath(output_dir)))) && error("Please provide a valid directory for output")

    # load common analysis code
    include(joinpath(source_dir, "exps", "common.jl"))

    # load specific analysis
    include(analysis_full_path)
    println("Done.")

    # compute output filename
    output_file = compute_name(output_dir)

    # Run
    println("Running $analysis_file and saving output in $output_file")
    func_name = Symbol(splitext(analysis_file)[1])
    getfield(Main, func_name)(output_file, debug=debug)
end

"""
bootstrap_config(config_file::String, output_dir::String, suppress_graph::Bool)

Run simulation using parameters supplied in a YAML configuration file
"""
function bootstrap_config(config_file::Compat.String,
                          output_dir::Compat.String,
                          suppress_graph::Bool)
    # load config
    config = YAML.load_file(config_file)

    if !isdir(output_dir)
        info("Directory $output_dir does not exist, attempting to create")
        try
            mkdir(output_dir)
        catch ex
            if isa(ex, SystemError)
                error("Unable to create $output_dir")
            end
        end
    end

    info("Using $(nprocs()) thread(s)")

    # load analysis files and simulation
    include(joinpath(source_dir, "simulation", "load.jl"))
    include(joinpath(source_dir, "utils", "parse_yaml.jl"))
    include(joinpath(source_dir, "utils", "runconfig.jl"))

    # run
    Base.invokelatest(runconfig, Base.invokelatest(parse, config)..., output_dir, suppress_graph)
end


text:
# ---
# title: 224. Basic Calculator
# id: problem224
# author: Qling
# date: 2021-03-10
# difficulty: Hard
# categories: Math, Stack
# link: <https://leetcode.com/problems/basic-calculator/description/>
# hidden: true
# ---
#
# Implement a basic calculator to evaluate a simple expression string.
#
# The expression string may contain open `(` and closing parentheses `)`, the
# plus `+` or minus sign `-`, **non-negative** integers and empty spaces ` `.
#
# **Example 1:**
#
#
#
#     Input: "1 + 1"
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: " 2-1 + 2 "
#     Output: 3
#
# **Example 3:**
#
#
#
#     Input: "(1+(4+5+2)-3)+(6+8)"
#     Output: 23
#
# **Note:**
#
#   * You may assume that the given expression is always valid.
#   * **Do not** use the `eval` built-in library function.
#
#
## @lc code=start
using LeetCode

function calculate(s::String)::Int
    res, num, sign, stk = 0, 0, 1, []
    for c in s
        if isdigit(c)
            num = 10 * num + parse(Int, c)
        elseif c in ['-', '+']
            res += num * sign
            num = 0
            sign = (c == '-') ? -1 : 1
        elseif c == '('
            push!(stk, res)
            push!(stk, sign)
            res, sign = 0, 1
        elseif c == ')'
            res += num * sign
            res *= pop!(stk)
            res += pop!(stk)
            num = 0
        end
    end

    return res + num * sign
end
## @lc code=end


text:
__precompile__(true)
module GenomicAnnotations

using DataFrames
using BioSequences
using GZip

export Chromosome, Gene, AbstractGene, GeneDataView, Locus
export readgbk, sequence, iscomplement, addgene!, pushproperty!, printgbk
export @genes

include("types.jl")
include("genedataview.jl")
include("readgbk.jl")
include("macro.jl")

end #module


text:
################################################################################

function calcNestPosForEquator(z, z_abs, scaled_phi)
    jp = floor(Integer, NSIDE_MAX * (0.5 + scaled_phi - z * 0.75))
    jm = floor(Integer, NSIDE_MAX * (0.5 + scaled_phi + z * 0.75))

    idfp = div(jp, NSIDE_MAX) # in {0,4}
    idfm = div(jm, NSIDE_MAX)

    local face_num
    if idfp == idfm
        face_num = (idfp % 4) + 4
    elseif idfp < idfm
        face_num = (idfp % 4)
    else
        face_num = (idfm % 4) + 8
    end

    ix = mod(jm, NSIDE_MAX)
    iy = NSIDE_MAX - mod(jp, NSIDE_MAX) - 1

    (ix, iy, face_num)
end

################################################################################

function calcNestPosForPole(z, z_abs, scaled_phi)
    ntt = floor(Integer, scaled_phi)
    if ntt >= 4
        ntt = 3
    end

    tp = scaled_phi - ntt
    tmp = sqrt(3 * (1 - z_abs)) # in ]0,1]

    jp = floor(Integer, NSIDE_MAX * tp * tmp)
    jm = floor(Integer, NSIDE_MAX * (1 - tp) * tmp)

    # Clip jp and jm
    jp = jp < NSIDE_MAX - 1 ? jp : NSIDE_MAX - 1
    jm = jm < NSIDE_MAX - 1 ? jm : NSIDE_MAX - 1

    local ix, iy, face_num
    if z >= 0
        face_num = ntt # in {0,3}
        ix = NSIDE_MAX - jm - 1
        iy = NSIDE_MAX - jp - 1
    else
        face_num = ntt + 8 # in {8,11} */
        ix =  jp
        iy =  jm
    end

    (ix, iy, face_num)
end

################################################################################

"""
    ang2pixNest(resol::Resolution, theta, phi) -> Integer

Return the index of the pixel which contains the point with
coordinates (`theta`, the colatitude, and `phi`, the longitude), in
radians, for a Healpix map with pixels in nested order. Note that
pixel indexes are 1-based (this is Julia)!
"""
function ang2pixNest(resol::Resolution, theta, phi)

    nside = resol.nside
    local ix, iy, face_num

    z          = cos(theta)
    z_abs      = abs(z)
    scaled_phi = mod2pi(phi) / (π / 2) # in [0,4[

    if z_abs ≤ 2//3
        (ix, iy, face_num) = calcNestPosForEquator(z, z_abs, scaled_phi)
    else
        (ix, iy, face_num) = calcNestPosForPole(z, z_abs, scaled_phi)
    end

    (ix_hi, ix_low) = divrem(ix, 128)
    (iy_hi, iy_low) = divrem(iy, 128)

    ipf = ((x2pix[ix_hi + 1] + y2pix[iy_hi + 1]) * (128^2)
           + (x2pix[ix_low + 1] + y2pix[iy_low + 1]))
    ipf = floor(Integer, ipf / ((NSIDE_MAX / nside) ^ 2))

    # Add 1 to have a 1-based index
    ipf + face_num * nside * nside + 1

end

################################################################################

function calcRingPosForEquator(resol::Resolution, z, z_abs, tt)

    jp = floor(Integer, resol.nside * (0.5 + tt - z * 0.75))
    jm = floor(Integer, resol.nside * (0.5 + tt + z * 0.75))

    ir = resol.nside + 1 + jp - jm
    kshift = (mod(ir, 2) == 0) ? 1 : 0

    nl4 = resol.nsideTimesFour;

    local ip = div(jp + jm - resol.nside + kshift + 1, 2) + 1
    if ip > nl4
	    ip = ip - nl4
    end

    resol.ncap + nl4 * (ir - 1) + ip
end

################################################################################

function calcRingPosForPole(resol::Resolution, z, z_abs, tt)

    tp = tt - floor(tt)
    tmp = sqrt(3. * (1. - z_abs))

    jp = floor(Integer, resol.nside * tp * tmp )
    jm = floor(Integer, resol.nside * (1 - tp) * tmp)

    ir = jp + jm + 1
    ip = floor(Integer, tt * ir) + 1
    if ip > 4ir
	    ip -= 4ir
    end

    if z ≤ 0
	    resol.numOfPixels - 2ir * (ir + 1) + ip
    else
        2ir * (ir - 1) + ip
    end
end

################################################################################

"""
    ang2vec(theta, phi) -> Array{Float64}

Given a direction in the sky with colatitude `theta` and longitude
`phi` (in radians), return an array of 3 elements containing the
`x`, `y`, and `z` components of the one-length vector pointing to
that direction.
"""
function ang2vec(theta, phi)
    (0 ≤ theta ≤ π) || throw(DomainError(theta, "Invalid value of theta"))

    sintheta = sin(theta)
    return [sintheta * cos(phi), sintheta * sin(phi), cos(theta)]
end

################################################################################

"""
    vec2ang(x, y, z) -> (Number, Number)

Given a vector (not necessarily normalized) whose Cartesian components are `x`,
`y`, and `z`, return a pair (`theta`, `phi`) containing the colatitude `theta`
and the longitude `phi` (in radians) of the direction in the sky the vector is
pointing at.
"""
function vec2ang(x, y, z)
    norm = sqrt(x^2 + y^2 + z^2)
    theta = acos(z / norm)
    phi = atan(y, x)
    if phi < 0
        phi += 2π
    end

    return (theta, phi)
end

################################################################################

"""
    ang2pixRing(resol::Resolution, theta, phi) -> Integer

Return the index of the pixel which contains the point with
coordinates (`theta`, the colatitude, and `phi`, the longitude), in
radians, for a Healpix map with pixels in ring order. Note that pixel
indexes are 1-based (this is Julia)!
"""
function ang2pixRing(resol::Resolution, theta, phi)

    z = cos(theta)
    z_abs = abs(z)
    scaled_phi = mod2pi(phi) / (π / 2) # in [0,4[

    if z_abs ≤ 2//3
        calcRingPosForEquator(resol, z, z_abs, scaled_phi)
    else
        calcRingPosForPole(resol, z, z_abs, scaled_phi)
    end
end

################################################################################

"""
    pix2angNest(resol::Resolution, pixel) -> (Float64, Float64)

Given the (1-based) index of a pixel in a Healpix map in nested
order, return a pair containing the (`colatitude`, `longitude`) angles
corresponding to its center, both expressed in radians.
"""
function pix2angNest(resol::Resolution, pixel)

    jrll = [ 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 ]
    jpll = [ 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 ]

    floatNside = Float64(resol.nside)
    fact1 = 1. / (3.0 * floatNside^2)
    fact2 = 2. / (3.0 * floatNside)

    # face number in {0,11} and pixel number within the face
    faceNum, ipf = divrem(pixel - 1, resol.pixelsPerFace)

    ip_trunc, ip_low = divrem(ipf, 1024)
    ip_hi, ip_med = divrem(ip_trunc, 1024)

    ix = 1024 * pix2x[ip_hi+1] + 32 * pix2x[ip_med+1] + pix2x[ip_low+1]
    iy = 1024 * pix2y[ip_hi+1] + 32 * pix2y[ip_med+1] + pix2y[ip_low+1]

    # Transforms this in (horizontal, vertical) coordinates
    jrt = ix + iy # 'vertical' in {0,2*(nside-1)}
    jpt = ix - iy # 'horizontal' in {-nside+1,nside-1}

    jr = jrll[faceNum+1] * resol.nside - jrt - 1
    local nr = resol.nside # Equatorial region (the most frequent)
    z = (2resol.nside - jr) * fact2
    local kshift = Int(mod(jr - resol.nside, 2))
    if jr < resol.nside
        # North polar cap
	nr = jr
	z = 1. - nr^2 * fact1
	kshift = 0
    elseif jr > 3 * resol.nside
        # South polar cap
	nr = resol.nsideTimesFour - jr
	z = - 1. + nr^2 * fact1
	kshift = 0
    end

    local jp = div(jpll[faceNum+1] * nr + jpt + 1 + kshift, 2)
    if jp > resol.nsideTimesFour
        jp = jp - resol.nsideTimesFour
    end
    if jp < 1
        jp = jp + resol.nsideTimesFour
    end

    (acos(z), (jp - (kshift + 1) * 0.5) * (π / (2nr)))
end

################################################################################

"""
    pix2ringpos(resol::Resolution, pixel)

Given the (1-based) index of a pixel in a Healpix map in ring
order, return a pair of numbers (n, i, j) whose meaning is the following:

- `n` can be one of the symbols `:northcap`, `:equator`, or
  `:southcap`, representing the region of the sky
- `i` is the ring index, from 1 to 4NSIDE - 1
- `j` is the pixel-in-ring index
"""
function pix2ringpos(resol::Resolution, pixel)
    # Any reference to equations in this routine refers to Gorski et al. (2005)
    if pixel ≤ resol.ncap
        # North polar cap

        p_h = pixel / 2 # Defined in Gorsky et al. (2005) before Eq. (2)
        floor_p_h = floor(p_h)
        # Eq. (2)
        i = floor(Int, sqrt(p_h - sqrt(floor_p_h))) + 1
        # Eq. (3)
        j = pixel - 2i * (i - 1)

        (:northcap, i, j)
    elseif pixel ≤ resol.nsideTimesTwo * (5resol.nside + 1)
        ip    = pixel - resol.ncap - 1
        # Eq. (6) - ring counts from the North pole; resol.nside is
        # the number of pixels in the North Polar ring
        i = floor(Int, ip / resol.nsideTimesFour) + resol.nside
        # Eq. (7) - zero-based index of the pixel within this ring
        j  = Int(mod(ip, resol.nsideTimesFour)) + 1
        (:equator, i, j)
    else
        # South polar cap

        ip = resol.numOfPixels - pixel + 1
        p_h = ip / 2
        floor_p_h = floor(p_h)
        i = floor(Int, sqrt(p_h - sqrt(floor_p_h))) + 1 # counted from S. pole
        j = Int(4 * i + 1 - (ip - 2i * (i - 1)))

        (:southcap, i, j)
    end
end

################################################################################

"""
    pix2angRing(resol::Resolution, pixel) -> (Float64, Float64)

Given the (1-based) index of a pixel in a Healpix map in ring
order, return a pair containing the (`colatitude`, `longitude`) angles
corresponding to its center, both expressed in radians.
"""
function pix2angRing(resol::Resolution, pixel)
    fact1 = 1.5 * resol.nside
    fact2 = 3.0 * resol.pixelsPerFace

    # Any reference to equations in this routine refers to Gorski et al. (2005)
    cap, i, j = pix2ringpos(resol, pixel)
    if cap == :northcap
        # Colatitude: Eq. (4); longitude: Eq. (5)
        return (acos(1 - i^2 / fact2),
                    (Float64(j) - 0.5) * π / (2i))
    elseif pixel ≤ resol.nsideTimesTwo * (5resol.nside + 1)
        # Equatorial belt

        # Eq. (9) - this equals 1 if i + resol.nside is odd, 1/2
        # otherwise. It is used to convert j into a longitude (since
        # pixel centers in odd rings are shifted with respect to
        # centers in even rings)
        s_half = 0.5 * (1 + mod(Float64(i + resol.nside), 2))

        # Colatitude: Eq. (8) in disguise, latitude: Eq. (9)
        return (acos((resol.nsideTimesTwo - i) / fact1),
                    (Float64(j) - s_half) * π / (2resol.nside))
    else
        # South Polar cap

        # The pixels in this cap are handled like the ones in the
        # North Polar cap, except that we must flip the value of "ip".
        return (acos(-1 + i^2 / fact2),
                    (float(j) - 0.5) * π / (2i))
    end
end

################################################################################

vec2pixNest(res::Resolution, x, y, z) = ang2pixNest(res, vec2ang(x, y, z)...)
vec2pixRing(res::Resolution, x, y, z) = ang2pixRing(res, vec2ang(x, y, z)...)
pix2vecNest(res::Resolution, pixel) = ang2vec(pix2angNest(res, pixel)...)
pix2vecRing(res::Resolution, pixel) = ang2vec(pix2angRing(res, pixel)...)

################################################################################

function ringWeightPath(datadir, nside)
    joinpath(datadir, @sprintf("weight_ring_n%05d.fits", nside))
end

function readWeightRing(fileName, nside)
    f = FITSIO.fits_open_table(fileName)
    try
        weights = Array(Float64, 2 * nside)
        FITSIO.fits_read_col(f, 1, 1, 1, weights)

        return weights
    finally
        FITSIO.fits_close_file(f)
    end
end

################################################################################

function pixelWindowPath(datadir, nside)
    joinpath(datadir, @sprintf("pixel_window_n%04d.fits", nside))
end

function readPixelWindowT(fileName, nside)
    f = FITSIO.fits_open_table(fileName)
    try
        pixwin = Array(Float64, FITSIO.fits_get_num_rows(f))
        FITSIO.fits_read_col(f, 1, 1, 1, pixwin)

        return pixwin
    finally
        FITSIO.fits_close_file(f)
    end
end

function readPixelWindowP(fileName, nside)
    local pixwinT
    local pixwinP

    f = FITSIO.fits_open_table(fileName)
    try
        pixwinT = Array(Float64, FITSIO.fits_get_num_rows(f))
        pixwinP = Array(Float64, FITSIO.fits_get_num_rows(f))
        FITSIO.fits_read_col(f, 1, 1, 1, pixwinT)
        FITSIO.fits_read_col(f, 2, 1, 1, pixwinP)

        return (pixwinT, pixwinP)
    finally
        FITSIO.fits_close_file(f)
    end
end

################################################################################

"""
    ang2pix{T, O <: Order}(map::Map{T, O}, theta::Real, phi::Real)

Convert the direction specified by the colatitude `theta` (∈ [0, π])
and the longitude `phi` (∈ [0, 2π]) into the index of the pixel in the
Healpix map `map`.
"""
function ang2pix(map::OrderedMap{T, Ring}, theta, phi) where {T}
    ang2pixRing(map.resolution, Float64(theta), Float64(phi))
end

function ang2pix(map::OrderedMap{T, Nest}, theta, phi) where {T}
    ang2pixNest(map.resolution, Float64(theta), Float64(phi))
end

################################################################################

"""
    pix2ang{T, O <: Order}(map::Map{T, O}, ipix) -> (Float64, Float64)

Return the pair (`theta`, `phi`), where `theta` is the colatitude and
`phi` the longitude of the direction of the pixel center with index
`ipix`.
"""
function pix2ang(map::OrderedMap{T, Ring}, ipix) where {T}
    pix2angRing(map.resolution, ipix)
end

function pix2ang(map::OrderedMap{T, Nest}, ipix) where {T}
    pix2angNest(map.resolution, ipix)
end

################################################################################