text:
"""
    reconstruct!(KS::SolverSet, ctr::AbstractArray)

Reconstruct solutions in cells
"""
function reconstruct!(KS::SolverSet, ctr::T) where {T<:AbstractArray{ControlVolume1D,1}}

    if KS.set.interpOrder == 1
        return
    end

    dx = KS.ps.dx

    if first(eachindex(KS.pSpace.x)) < 1
        idx0 = 1
        idx1 = KS.pSpace.nx
    else
        idx0 = 2
        idx1 = KS.pSpace.nx - 1
    end

    if ctr[1].w isa Number
        @inbounds Threads.@threads for i = idx0:idx1
            ctr[i].sw = reconstruct3(
                ctr[i-1].w,
                ctr[i].w,
                ctr[i+1].w,
                0.5 * (dx[i-1] + dx[i]),
                0.5 * (dx[i] + dx[i+1]),
                Symbol(KS.set.limiter),
            )
        end
    else
        @inbounds Threads.@threads for i = idx0:idx1
            reconstruct3!(
                ctr[i].sw,
                ctr[i-1].w,
                ctr[i].w,
                ctr[i+1].w,
                0.5 * (dx[i-1] + dx[i]),
                0.5 * (dx[i] + dx[i+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

    return nothing

end

function reconstruct!(KS::SolverSet, ctr::T) where {T<:AbstractArray{ControlVolume1D1F,1}}

    if KS.set.interpOrder == 1
        return
    end

    dx = KS.ps.dx

    if first(eachindex(KS.pSpace.x)) < 1
        idx0 = 1
        idx1 = KS.pSpace.nx
    else
        idx0 = 2
        idx1 = KS.pSpace.nx - 1
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sw,
            ctr[i-1].w,
            ctr[i].w,
            ctr[i+1].w,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sf,
            ctr[i-1].f,
            ctr[i].f,
            ctr[i+1].f,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

end

function reconstruct!(KS::SolverSet, ctr::T) where {T<:AbstractArray{ControlVolume1D2F,1}}

    if KS.set.interpOrder == 1
        return
    end

    dx = KS.ps.dx

    if first(eachindex(KS.pSpace.x)) < 1
        idx0 = 1
        idx1 = KS.pSpace.nx
    else
        idx0 = 2
        idx1 = KS.pSpace.nx - 1
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sw,
            ctr[i-1].w,
            ctr[i].w,
            ctr[i+1].w,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sh,
            ctr[i-1].h,
            ctr[i].h,
            ctr[i+1].h,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
        reconstruct3!(
            ctr[i].sb,
            ctr[i-1].b,
            ctr[i].b,
            ctr[i+1].b,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

end

function reconstruct!(KS::SolverSet, ctr::T) where {T<:AbstractArray{ControlVolume1D3F,1}}

    if KS.set.interpOrder == 1
        return
    end

    dx = KS.ps.dx

    if first(eachindex(KS.pSpace.x)) < 1
        idx0 = 1
        idx1 = KS.pSpace.nx
    else
        idx0 = 2
        idx1 = KS.pSpace.nx - 1
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sw,
            ctr[i-1].w,
            ctr[i].w,
            ctr[i+1].w,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sh0,
            ctr[i-1].h0,
            ctr[i].h0,
            ctr[i+1].h0,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
        reconstruct3!(
            ctr[i].sh1,
            ctr[i-1].h1,
            ctr[i].h1,
            ctr[i+1].h1,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
        reconstruct3!(
            ctr[i].sh2,
            ctr[i-1].h2,
            ctr[i].h2,
            ctr[i+1].h2,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

end

function reconstruct!(KS::SolverSet, ctr::T) where {T<:AbstractArray{ControlVolume1D4F,1}}

    if KS.set.interpOrder == 1
        return
    end

    dx = KS.ps.dx

    if first(eachindex(KS.pSpace.x)) < 1
        idx0 = 1
        idx1 = KS.pSpace.nx
    else
        idx0 = 2
        idx1 = KS.pSpace.nx - 1
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sw,
            ctr[i-1].w,
            ctr[i].w,
            ctr[i+1].w,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

    @inbounds Threads.@threads for i = idx0:idx1
        reconstruct3!(
            ctr[i].sh0,
            ctr[i-1].h0,
            ctr[i].h0,
            ctr[i+1].h0,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
        reconstruct3!(
            ctr[i].sh1,
            ctr[i-1].h1,
            ctr[i].h1,
            ctr[i+1].h1,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
        reconstruct3!(
            ctr[i].sh2,
            ctr[i-1].h2,
            ctr[i].h2,
            ctr[i+1].h2,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
        reconstruct3!(
            ctr[i].sh3,
            ctr[i-1].h3,
            ctr[i].h3,
            ctr[i+1].h3,
            0.5 * (dx[i-1] + dx[i]),
            0.5 * (dx[i] + dx[i+1]),
            Symbol(KS.set.limiter),
        )
    end

end

# ------------------------------------------------------------
# 2D1F
# ------------------------------------------------------------
function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{ControlVolume2D,2}}

    if KS.set.interpOrder == 1
        return
    end

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    #--- conservative variables ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        swL = extract_last(ctr[1, j].sw, 1, mode = :view)
        reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j]))

        swR = extract_last(ctr[nx, j].sw, 1, mode = :view)
        reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dy[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        swD = extract_last(ctr[i, 1].sw, 2, mode = :view)
        reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2]))

        swU = extract_last(ctr[i, ny].sw, 2, mode = :view)
        reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny]))
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sw = extract_last(ctr[i, j].sw, 1, mode = :view)
            reconstruct3!(
                sw,
                ctr[i-1, j].w,
                ctr[i, j].w,
                ctr[i+1, j].w,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sw = extract_last(ctr[i, j].sw, 2, mode = :view)
            reconstruct3!(
                sw,
                ctr[i, j-1].w,
                ctr[i, j].w,
                ctr[i, j+1].w,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

end

function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{ControlVolume2D1F,2}}

    if KS.set.interpOrder == 1
        return
    end

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    #--- conservative variables ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        swL = extract_last(ctr[1, j].sw, 1, mode = :view)
        reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j]))

        swR = extract_last(ctr[nx, j].sw, 1, mode = :view)
        reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dx[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        swD = extract_last(ctr[i, 1].sw, 2, mode = :view)
        reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2]))

        swU = extract_last(ctr[i, ny].sw, 2, mode = :view)
        reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny]))
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sw = extract_last(ctr[i, j].sw, 1, mode = :view)
            reconstruct3!(
                sw,
                ctr[i-1, j].w,
                ctr[i, j].w,
                ctr[i+1, j].w,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sw = extract_last(ctr[i, j].sw, 2, mode = :view)
            reconstruct3!(
                sw,
                ctr[i, j-1].w,
                ctr[i, j].w,
                ctr[i, j+1].w,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

    #--- particle distribution function ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        sfL = extract_last(ctr[1, j].sf, 1, mode = :view)
        reconstruct2!(sfL, ctr[1, j].f, ctr[2, j].f, 0.5 * (dx[1, j] + dx[2, j]))

        sfR = extract_last(ctr[nx, j].sf, 1, mode = :view)
        reconstruct2!(sfR, ctr[nx-1, j].f, ctr[nx, j].f, 0.5 * (dx[nx-1, j] + dx[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        sfD = extract_last(ctr[i, 1].sf, 2, mode = :view)
        reconstruct2!(sfD, ctr[i, 1].f, ctr[i, 2].f, 0.5 * (dy[i, 1] + dy[i, 2]))

        sfU = extract_last(ctr[i, ny].sf, 2, mode = :view)
        reconstruct2!(sfU, ctr[i, ny-1].f, ctr[i, ny].f, 0.5 * (dy[i, ny-1] + dy[i, ny]))
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sf = extract_last(ctr[i, j].sf, 1, mode = :view)
            reconstruct3!(
                sf,
                ctr[i-1, j].f,
                ctr[i, j].f,
                ctr[i+1, j].f,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sf = extract_last(ctr[i, j].sf, 2, mode = :view)
            reconstruct3!(
                sf,
                ctr[i, j-1].f,
                ctr[i, j].f,
                ctr[i, j+1].f,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

end

function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{ControlVolume2D2F,2}}

    if KS.set.interpOrder == 1
        return
    end

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    #--- conservative variables ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        swL = extract_last(ctr[1, j].sw, 1, mode = :view)
        reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j]))

        swR = extract_last(ctr[nx, j].sw, 1, mode = :view)
        reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dx[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        swD = extract_last(ctr[i, 1].sw, 2, mode = :view)
        reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2]))

        swU = extract_last(ctr[i, ny].sw, 2, mode = :view)
        reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny]))
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sw = extract_last(ctr[i, j].sw, 1, mode = :view)
            reconstruct3!(
                sw,
                ctr[i-1, j].w,
                ctr[i, j].w,
                ctr[i+1, j].w,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sw = extract_last(ctr[i, j].sw, 2, mode = :view)
            reconstruct3!(
                sw,
                ctr[i, j-1].w,
                ctr[i, j].w,
                ctr[i, j+1].w,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

    #--- particle distribution function ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        shL = extract_last(ctr[1, j].sh, 1, mode = :view)
        reconstruct2!(shL, ctr[1, j].h, ctr[2, j].h, 0.5 * (dx[1, j] + dx[2, j]))
        sbL = extract_last(ctr[1, j].sb, 1, mode = :view)
        reconstruct2!(sbL, ctr[1, j].b, ctr[2, j].b, 0.5 * (dx[1, j] + dx[2, j]))

        shR = extract_last(ctr[nx, j].sh, 1, mode = :view)
        reconstruct2!(shR, ctr[nx-1, j].h, ctr[nx, j].h, 0.5 * (dx[nx-1, j] + dx[nx, j]))
        sbR = extract_last(ctr[nx, j].sb, 1, mode = :view)
        reconstruct2!(sbR, ctr[nx-1, j].b, ctr[nx, j].b, 0.5 * (dx[nx-1, j] + dx[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        shD = extract_last(ctr[i, 1].sh, 2, mode = :view)
        reconstruct2!(shD, ctr[i, 1].h, ctr[i, 2].h, 0.5 * (dy[i, 1] + dy[i, 2]))
        sbD = extract_last(ctr[i, 1].sb, 2, mode = :view)
        reconstruct2!(sbD, ctr[i, 1].b, ctr[i, 2].b, 0.5 * (dy[i, 1] + dy[i, 2]))

        shU = extract_last(ctr[i, ny].sh, 2, mode = :view)
        reconstruct2!(shU, ctr[i, ny-1].h, ctr[i, ny].h, 0.5 * (dy[i, ny-1] + dy[i, ny]))
        sbU = extract_last(ctr[i, ny].sb, 2, mode = :view)
        reconstruct2!(sbU, ctr[i, ny-1].b, ctr[i, ny].b, 0.5 * (dy[i, ny-1] + dy[i, ny]))
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sh = extract_last(ctr[i, j].sh, 1, mode = :view)
            reconstruct3!(
                sh,
                ctr[i-1, j].h,
                ctr[i, j].h,
                ctr[i+1, j].h,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )

            sb = extract_last(ctr[i, j].sb, 1, mode = :view)
            reconstruct3!(
                sb,
                ctr[i-1, j].b,
                ctr[i, j].b,
                ctr[i+1, j].b,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sh = extract_last(ctr[i, j].sh, 2, mode = :view)
            reconstruct3!(
                sh,
                ctr[i, j-1].h,
                ctr[i, j].h,
                ctr[i, j+1].h,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )

            sb = extract_last(ctr[i, j].sb, 2, mode = :view)
            reconstruct3!(
                sb,
                ctr[i, j-1].b,
                ctr[i, j].b,
                ctr[i, j+1].b,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

end

function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractArray{ControlVolume2D3F,2}}

    if KS.set.interpOrder == 1
        return
    end

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    #--- macroscopic variables ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        swL = extract_last(ctr[1, j].sw, 1, mode = :view)
        reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j]))

        swR = extract_last(ctr[nx, j].sw, 1, mode = :view)
        reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dx[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        swD = extract_last(ctr[i, 1].sw, 2, mode = :view)
        reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2]))

        swU = extract_last(ctr[i, ny].sw, 2, mode = :view)
        reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny]))
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sw = extract_last(ctr[i, j].sw, 1, mode = :view)
            reconstruct3!(
                sw,
                ctr[i-1, j].w,
                ctr[i, j].w,
                ctr[i+1, j].w,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sw = extract_last(ctr[i, j].sw, 2, mode = :view)
            reconstruct3!(
                sw,
                ctr[i, j-1].w,
                ctr[i, j].w,
                ctr[i, j+1].w,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

    #--- particle distribution function ---#
    # boundary
    @inbounds Threads.@threads for j = 1:ny
        s0L = extract_last(ctr[1, j].sh0, 1, mode = :view)
        reconstruct2!(s0L, ctr[1, j].h0, ctr[2, j].h0, 0.5 * (dx[1, j] + dx[2, j]))
        s1L = extract_last(ctr[1, j].sh1, 1, mode = :view)
        reconstruct2!(s1L, ctr[1, j].h1, ctr[2, j].h1, 0.5 * (dx[1, j] + dx[2, j]))
        s2L = extract_last(ctr[1, j].sh2, 1, mode = :view)
        reconstruct2!(s2L, ctr[1, j].h2, ctr[2, j].h2, 0.5 * (dx[1, j] + dx[2, j]))

        s0R = extract_last(ctr[nx, j].sh0, 1, mode = :view)
        reconstruct2!(s0R, ctr[nx-1, j].h0, ctr[nx, j].h0, 0.5 * (dx[nx-1, j] + dx[nx, j]))
        s1R = extract_last(ctr[nx, j].sh1, 1, mode = :view)
        reconstruct2!(s1R, ctr[nx-1, j].h1, ctr[nx, j].h1, 0.5 * (dx[nx-1, j] + dx[nx, j]))
        s2R = extract_last(ctr[nx, j].sh2, 1, mode = :view)
        reconstruct2!(s2R, ctr[nx-1, j].h2, ctr[nx, j].h2, 0.5 * (dx[nx-1, j] + dx[nx, j]))
    end

    @inbounds Threads.@threads for i = 1:nx
        s0D = extract_last(ctr[i, 1].sh0, 2, mode = :view)
        reconstruct2!(s0D, ctr[i, 1].h0, ctr[i, 2].h0, 0.5 * (dy[i, 1] + dy[i, 2]))
        s1D = extract_last(ctr[i, 1].sh1, 2, mode = :view)
        reconstruct2!(s1D, ctr[i, 1].h1, ctr[i, 2].h1, 0.5 * (dy[i, 1] + dy[i, 2]))
        s2D = extract_last(ctr[i, 1].sh2, 2, mode = :view)
        reconstruct2!(s2D, ctr[i, 1].h2, ctr[i, 2].h2, 0.5 * (dy[i, 1] + dy[i, 2]))

        s0U = extract_last(ctr[i, ny].sh0, 2, mode = :view)
        reconstruct2!(s0U, ctr[i, ny-1].h0, ctr[i, ny].h0, 0.5 * (dy[i, ny-1] + dy[i, ny]))
        s1U = extract_last(ctr[i, KS.pSpace.ny].sh1, 2, mode = :view)
        reconstruct2!(
            s1U,
            ctr[i, KS.pSpace.ny-1].h1,
            ctr[i, KS.pSpace.ny].h1,
            0.5 * (dy[i, ny-1] + dy[i, ny]),
        )
        s2U = extract_last(ctr[i, KS.pSpace.ny].sh2, 2, mode = :view)
        reconstruct2!(
            s2U,
            ctr[i, KS.pSpace.ny-1].h2,
            ctr[i, KS.pSpace.ny].h2,
            0.5 * (dy[i, ny-1] + dy[i, ny]),
        )
    end

    # inner
    @inbounds Threads.@threads for j = 1:ny
        for i = 2:nx-1
            sh0 = extract_last(ctr[i, j].sh0, 1, mode = :view)
            reconstruct3!(
                sh0,
                ctr[i-1, j].h0,
                ctr[i, j].h0,
                ctr[i+1, j].h0,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )

            sh1 = extract_last(ctr[i, j].sh1, 1, mode = :view)
            reconstruct3!(
                sh1,
                ctr[i-1, j].h1,
                ctr[i, j].h1,
                ctr[i+1, j].h1,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )

            sh2 = extract_last(ctr[i, j].sh2, 1, mode = :view)
            reconstruct3!(
                sh2,
                ctr[i-1, j].h2,
                ctr[i, j].h2,
                ctr[i+1, j].h2,
                0.5 * (dx[i-1, j] + dx[i, j]),
                0.5 * (dx[i, j] + dx[i+1, j]),
                Symbol(KS.set.limiter),
            )
        end
    end

    @inbounds Threads.@threads for j = 2:ny-1
        for i = 1:nx
            sh0 = extract_last(ctr[i, j].sh0, 2, mode = :view)
            reconstruct3!(
                sh0,
                ctr[i, j-1].h0,
                ctr[i, j].h0,
                ctr[i, j+1].h0,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )

            sh1 = extract_last(ctr[i, j].sh1, 2, mode = :view)
            reconstruct3!(
                sh1,
                ctr[i, j-1].h1,
                ctr[i, j].h1,
                ctr[i, j+1].h1,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )

            sh2 = extract_last(ctr[i, j].sh2, 2, mode = :view)
            reconstruct3!(
                sh2,
                ctr[i, j-1].h2,
                ctr[i, j].h2,
                ctr[i, j+1].h2,
                0.5 * (dy[i, j-1] + dy[i, j]),
                0.5 * (dy[i, j] + dy[i, j+1]),
                Symbol(KS.set.limiter),
            )
        end
    end

end

function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractVector{ControlVolumeUS}}

    if KS.set.interpOrder == 1
        return
    end

    @inbounds Threads.@threads for i in eachindex(ctr)
        ids = ps.cellNeighbors[i, :]
        deleteat!(ids, findall(x -> x == -1, ids))
        id1, id2 = ids[1:2]

        swx = KitBase.extract_last(ctr[i].sw, 1, mode = :view)
        swy = KitBase.extract_last(ctr[i].sw, 2, mode = :view)

        wL = ctr[id1].w
        wR = ctr[id2].w

        dxL = ctr[i].x[1] - ctr[id1].x[1]
        dxR = ctr[id2].x[1] - ctr[i].x[1]
        dyL = ctr[i].x[2] - ctr[id1].x[2]
        dyR = ctr[id2].x[2] - ctr[i].x[2]

        reconstruct3!(swx, wL, ctr[i].w, wR, dxL, dxR, Symbol(KS.set.limiter))
        reconstruct3!(swy, wL, ctr[i].w, wR, dyL, dyR, Symbol(KS.set.limiter))
    end

end

function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractVector{ControlVolumeUS1F}}

    if KS.set.interpOrder == 1
        return
    end

    @inbounds Threads.@threads for i in eachindex(ctr)
        ids = ps.cellNeighbors[i, :]
        deleteat!(ids, findall(x -> x == -1, ids))
        id1, id2 = ids[1:2]

        swx = KitBase.extract_last(ctr[i].sw, 1, mode = :view)
        swy = KitBase.extract_last(ctr[i].sw, 2, mode = :view)
        sfx = KitBase.extract_last(ctr[i].sf, 1, mode = :view)
        sfy = KitBase.extract_last(ctr[i].sf, 2, mode = :view)

        wL = ctr[id1].w
        wR = ctr[id2].w
        fL = ctr[id1].f
        fR = ctr[id2].f

        dxL = ctr[i].x[1] - ctr[id1].x[1]
        dxR = ctr[id2].x[1] - ctr[i].x[1]
        dyL = ctr[i].x[2] - ctr[id1].x[2]
        dyR = ctr[id2].x[2] - ctr[i].x[2]

        reconstruct3!(swx, wL, ctr[i].w, wR, dxL, dxR, Symbol(KS.set.limiter))
        reconstruct3!(swy, wL, ctr[i].w, wR, dyL, dyR, Symbol(KS.set.limiter))
        reconstruct3!(sfx, fL, ctr[i].f, fR, dxL, dxR, Symbol(KS.set.limiter))
        reconstruct3!(sfy, fL, ctr[i].f, fR, dyL, dyR, Symbol(KS.set.limiter))
    end

end

function reconstruct!(
    KS::X,
    ctr::Y,
) where {X<:AbstractSolverSet,Y<:AbstractVector{ControlVolumeUS2F}}

    if KS.set.interpOrder == 1
        return
    end

    @inbounds Threads.@threads for i in eachindex(ctr)
        ids = KS.pSpace.cellNeighbors[i, :]
        if -1 in ids
            deleteat!(ids, findall(x -> x == -1, ids))
            id1, id2 = ids[1:2]

            swx = KitBase.extract_last(ctr[i].sw, 1, mode = :view)
            swy = KitBase.extract_last(ctr[i].sw, 2, mode = :view)
            shx = KitBase.extract_last(ctr[i].sh, 1, mode = :view)
            shy = KitBase.extract_last(ctr[i].sh, 2, mode = :view)
            sbx = KitBase.extract_last(ctr[i].sb, 1, mode = :view)
            sby = KitBase.extract_last(ctr[i].sb, 2, mode = :view)

            wL = ctr[id1].w
            wR = ctr[id2].w
            hL = ctr[id1].h
            hR = ctr[id2].h
            bL = ctr[id1].b
            bR = ctr[id2].b

            dxL = ctr[i].x[1] - ctr[id1].x[1]
            dxR = ctr[id2].x[1] - ctr[i].x[1]
            dyL = ctr[i].x[2] - ctr[id1].x[2]
            dyR = ctr[id2].x[2] - ctr[i].x[2]

            reconstruct3!(swx, wL, ctr[i].w, wR, dxL, dxR, Symbol(KS.set.limiter))
            reconstruct3!(swy, wL, ctr[i].w, wR, dyL, dyR, Symbol(KS.set.limiter))
            reconstruct3!(shx, hL, ctr[i].h, hR, dxL, dxR, Symbol(KS.set.limiter))
            reconstruct3!(shy, hL, ctr[i].h, hR, dyL, dyR, Symbol(KS.set.limiter))
            reconstruct3!(sbx, bL, ctr[i].b, bR, dxL, dxR, Symbol(KS.set.limiter))
            reconstruct3!(sby, bL, ctr[i].b, bR, dyL, dyR, Symbol(KS.set.limiter))
        else
            id1, id2, id3 = ids

            swx = KitBase.extract_last(ctr[i].sw, 1, mode = :view)
            swy = KitBase.extract_last(ctr[i].sw, 2, mode = :view)
            shx = KitBase.extract_last(ctr[i].sh, 1, mode = :view)
            shy = KitBase.extract_last(ctr[i].sh, 2, mode = :view)
            sbx = KitBase.extract_last(ctr[i].sb, 1, mode = :view)
            sby = KitBase.extract_last(ctr[i].sb, 2, mode = :view)

            w1 = ctr[id1].w
            w2 = ctr[id2].w
            w3 = ctr[id3].w
            h1 = ctr[id1].h
            h2 = ctr[id2].h
            h3 = ctr[id3].h
            b1 = ctr[id1].b
            b2 = ctr[id2].b
            b3 = ctr[id3].b

            reconstruct4!(
                swx,
                ctr[i].w,
                w1,
                w2,
                w3,
                ctr[i].x[1] - ctr[id1].x[1],
                ctr[i].x[1] - ctr[id2].x[1],
                ctr[i].x[1] - ctr[id3].x[1],
                Symbol(KS.set.limiter),
            )
            reconstruct4!(
                swy,
                ctr[i].w,
                w1,
                w2,
                w3,
                ctr[i].x[2] - ctr[id1].x[2],
                ctr[i].x[2] - ctr[id2].x[2],
                ctr[i].x[2] - ctr[id3].x[2],
                Symbol(KS.set.limiter),
            )
            reconstruct4!(
                shx,
                ctr[i].h,
                h1,
                h2,
                h3,
                ctr[i].x[1] - ctr[id1].x[1],
                ctr[i].x[1] - ctr[id2].x[1],
                ctr[i].x[1] - ctr[id3].x[1],
                Symbol(KS.set.limiter),
            )
            reconstruct4!(
                shy,
                ctr[i].h,
                h1,
                h2,
                h3,
                ctr[i].x[2] - ctr[id1].x[2],
                ctr[i].x[2] - ctr[id2].x[2],
                ctr[i].x[2] - ctr[id3].x[2],
                Symbol(KS.set.limiter),
            )
            reconstruct4!(
                sbx,
                ctr[i].b,
                b1,
                b2,
                b3,
                ctr[i].x[1] - ctr[id1].x[1],
                ctr[i].x[1] - ctr[id2].x[1],
                ctr[i].x[1] - ctr[id3].x[1],
                Symbol(KS.set.limiter),
            )
            reconstruct4!(
                sby,
                ctr[i].b,
                b1,
                b2,
                b3,
                ctr[i].x[2] - ctr[id1].x[2],
                ctr[i].x[2] - ctr[id2].x[2],
                ctr[i].x[2] - ctr[id3].x[2],
                Symbol(KS.set.limiter),
            )

            #=w = (ctr[id1].w, ctr[id2].w, ctr[id3].w)
            h = (ctr[id1].h, ctr[id2].h, ctr[id3].h)
            b = (ctr[id1].b, ctr[id2].b, ctr[id3].b)

            s = [
                (ctr[i].w[end] .- w[1][end]),
                (ctr[i].w[end] .- w[2][end]),
                (ctr[i].w[end] .- w[3][end]),
            ]

            # x
            idxs = deepcopy(ids)
            Δx = [
                (ctr[i].x[1] - ctr[id1].x[1]),
                (ctr[i].x[1] - ctr[id2].x[1]),
                (ctr[i].x[1] - ctr[id3].x[1]),
            ]
            maxid = findmax(abs.(s ./ Δx))[2]
            rmid = findall(x -> x == maxid, idxs)
            deleteat!(idxs, rmid)
            deleteat!(Δx, rmid)

            reconstruct3!(
                swx,
                ctr[idxs[1]].w,
                ctr[i].w,
                ctr[idxs[2]].w,
                Δx[1],
                -Δx[2],
                Symbol(KS.set.limiter),
            )
            reconstruct3!(
                shx,
                ctr[idxs[1]].h,
                ctr[i].h,
                ctr[idxs[2]].h,
                Δx[1],
                -Δx[2],
                Symbol(KS.set.limiter),
            )
            reconstruct3!(
                sbx,
                ctr[idxs[1]].b,
                ctr[i].b,
                ctr[idxs[2]].b,
                Δx[1],
                -Δx[2],
                Symbol(KS.set.limiter),
            )

            # y
            idys = deepcopy(ids)
            Δy = [
                (ctr[i].x[2] - ctr[id1].x[2]),
                (ctr[i].x[2] - ctr[id2].x[2]),
                (ctr[i].x[2] - ctr[id3].x[2]),
            ]
            maxid = findmax(abs.(s ./ Δy))[2]
            rmid = findall(x -> x == maxid, idys)
            deleteat!(idys, rmid)
            deleteat!(Δy, rmid)

            reconstruct3!(
                swy,
                ctr[idys[1]].w,
                ctr[i].w,
                ctr[idys[2]].w,
                Δy[1],
                -Δy[2],
                Symbol(KS.set.limiter),
            )
            reconstruct3!(
                shy,
                ctr[idys[1]].h,
                ctr[i].h,
                ctr[idys[2]].h,
                Δy[1],
                -Δy[2],
                Symbol(KS.set.limiter),
            )
            reconstruct3!(
                sby,
                ctr[idys[1]].b,
                ctr[i].b,
                ctr[idys[2]].b,
                Δy[1],
                -Δy[2],
                Symbol(KS.set.limiter),
            )=#
        end
    end

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

module Logging

using Base.Threads

@enum LogLevel LOG_NONE=0 LOG_ERROR=1 LOG_WARNING=2 LOG_INFO=3 LOG_DEBUG=4 LOG_TRACE=5 LOG_ALL=100
export LogLevel
export LOG_NONE, LOG_ERROR, LOG_WARNING, LOG_INFO, LOG_DEBUG, LOG_TRACE, LOG_ALL

import Base: +,-

function +(level::LogLevel, i::Integer)
    l = Int(level)
    l_min = Int(LOG_ERROR)
    l_max = Int(LOG_TRACE)
    ifelse(
        l_min <= l <= l_max,
        LogLevel(clamp(l + i, l_min, l_max)),
        level
    )
end

-(level::LogLevel, i::Integer) = level + (-i)

const log_colors = Dict(
    LOG_ERROR => Base.error_color(),
    LOG_WARNING => Base.warn_color(),
    LOG_INFO => Base.info_color(),
    LOG_DEBUG => :light_green,
    LOG_TRACE => :green,
)

const log_prefix = Dict(
    LOG_ERROR => "ERROR",
    LOG_WARNING => "WARNING",
    LOG_INFO => "INFO",
    LOG_DEBUG => "DEBUG",
    LOG_TRACE => "TRACE",
)

const _global_lock = RecursiveSpinLock()

const _output_io = Ref{IO}()

function output_logging_msg(level::LogLevel, msg...)
    io = _output_io[]
    color = log_colors[level]
    prefix = log_prefix[level]
    Base.print_with_color(color, io, prefix, " ($(myid()), $(threadid())): "; bold = true)
    Base.println_with_color(color, io, chomp(string(msg...)))
    Base.flush(io)
    nothing
end

export @enable_logging

function enable_logging_macro()
    quote
        _log_level = Base.Threads.Atomic{Int}(BAT.Logging.LOG_INFO)
    end
end

macro enable_logging()
    esc(enable_logging_macro())
end

# get_log_level_macro() = :(LogLevel(_log_level[]))
#
# macro get_log_level()
#     esc(get_log_level_macro())
# end
#
# set_log_level_macro(level) = :(_log_level[] = $level; $(get_log_level_macro()))
#
# macro set_log_level!(level)
#     esc(set_log_level_macro(level))
# end

export get_log_level, set_log_level!

get_log_level(m::Module) = m._log_level[]

function set_log_level!(m::Module, level::LogLevel)
    atomic_level = m._log_level
    atomic_level[] = level
    LogLevel(atomic_level[])
end

export @log_msg
export @log_error, @log_warning, @log_info, @log_debug, @log_trace

function logging_macro(level, msg)
    l_var = gensym("l_")
    m_var = gensym("m_")
    quote
        let $l_var = $level
            if _log_level[] >= Int($l_var) && $l_var != LOG_NONE
                # lock early in case message-generation code is not thread-safe
                lock(BAT.Logging._global_lock)
                try
                    let $m_var = $msg
                        BAT.Logging.output_logging_msg($l_var, $m_var)
                    end
                finally
                    unlock(BAT.Logging._global_lock)
                end
            end
        end
    end
end

macro log_msg(log_level, msg)
    esc(logging_macro(log_level, msg))
end

for (macro_name, log_level) in (
    (:log_error, LOG_ERROR),
    (:log_warning, LOG_WARNING),
    (:log_info, LOG_INFO),
    (:log_debug, LOG_DEBUG),
    (:log_trace, LOG_TRACE)
)

    @eval macro $macro_name(msg)
        esc(logging_macro($log_level, msg))
    end
end

function __init__()
    # STDERR object at runtime differs from precompilation:
    _output_io[] = STDERR
end

end # module

ext: jl


text:
using JuMP, Cbc, CSV

lightModel = Model(Cbc.Optimizer)

@variable(lightModel, fixtures[1:5], binary=true)
@variable(lightModel, 0 <= luminosity[1:5] <= 300)

@objective(lightModel, Min, 450*sum(fixtures)+2*sum(luminosity))

corridors = [1 2; 1 3; 2 3; 2 4; 3 5; 4 5]
minLuminosity = [250,300,150,200,350,180]

for y = 1:size(fixtures,1)
    @constraint(lightModel, luminosity[y] <= 300*fixtures[y])
end

for c = 1:size(corridors,1)
    @constraint(lightModel, luminosity[corridors[c,1]]+luminosity[corridors[c,2]] >= minLuminosity[c])
end

result = JuMP.optimize!(lightModel)

JuMP.value.(fixtures)
JuMP.value.(luminosity)

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
module test_Plot1

import ModiaMath

# Desired:
#   using Test
#
# In order that Test needs not to be defined in the user environment, it is included via ModiaMath:
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
    t = linspace(0.0, 10.0, 100)
else
    using ModiaMath.Test
    t = range(0.0, stop=10.0, length=100)
end

result = Dict{AbstractString,Any}()

result["time"] = t
result["phi"]  = sin.(t)

ModiaMath.plot(result, :phi, heading="Sine(time)")

# Test also that warn works for non-existing variables
println("... Next plot should give a warning:")
ModiaMath.plot(result, :SignalNotDefined, heading="Sine(time)", figure=2)

end


text:
"""
    get_signals(x)

Signals can be recursively retrieved from any type containing widgets.
"""
function get_signals end

get_signals(xs::Tuple) = begin
    sig_array = []
    for x in xs
        if x != nothing && x != tuple()
            sig = get_signals(x)
            if is_usable(sig)
                push!(sig_array, sig)
            end
        end
    end
    mapsig(sig_array)
end
get_signals(x::OrderedDict) = begin sigs = []
    used_keys = []
    k = keys(x)
    for key in k
        if x != nothing && x != tuple()
            sig = get_signals(x[key])
            if is_usable(sig)
                push!(sigs, sig)
                push!(used_keys, key)
            end
        end
    end
    if length(sigs) > 0
        return map(s::Vararg -> OrderedDict(zip(used_keys, s)), sigs...)
    else
        return nothing
    end
end
get_signals(xs::AbstractArray) =
    map(s::Vararg -> reshape([s...], size(xs)), get_signals.(xs)...)
@require Interact begin
    get_signals(x::Signal) = x
    get_signals(x::Interact.InputWidget) = signal(x)
    is_usable(sig) = isa(sig, Signal)
    mapsig(sig_array) = map(sigs::Vararg -> sigs, map(s -> signal(s), sig_array)...)
end
@require InteractNext begin
    get_signals(x::Observable) = x
    get_signals(x::WebIO.Scope) = observe(x)
    mapsig(sig_array) = begin
        map(sigs::Vararg -> sigs, sig_array[1], sig_array...)
    end
    is_usable(sig) = isa(sig, Observable)
end

"""
    apply_signals!(data, signal)

Signals can be recursively applied back to data.
"""
function apply_signals! end

apply_signals!(data::Associative, signal::Associative) = begin
    for key in keys(signal)
        data[key] = apply_signals!(data[key], signal[key])
    end
    return data
end
apply_signals!(data::Tuple, signal::Tuple) = apply_signals!.(data, signal)
@require StaticArrays begin
apply_signals!(data::FieldVector, signal) = data .= signal
end
apply_signals!(data, signal::Associative) = begin
    for key in keys(signal)
        setfield!(data, key, apply_signals!(getfield(data, key), signal[key]))
    end
    return data
end
apply_signals!(data, signal) = signal


text:
using BlockMatching
using Images, TestImages
using Documenter, DemoCards

format = Documenter.HTML(;
    prettyurls=get(ENV, "CI", "false") == "true",
    canonical="https://johnnychen94.github.io/BlockMatching.jl",
    assets=String[],
)

makedocs(;
    modules=[BlockMatching],
    sitename="BlockMatching.jl",
    format=format,
    pages=[
        "Home" => "index.md",
        "References" => "reference.md"
    ],
)

deploydocs(;
    repo="github.com/johnnychen94/BlockMatching.jl",
)


text:
module ConstrainedParameter

using Bijectors
using Bijectors: Bijector, Identity, Logit
import Functors: functor

export AbstractParameter, Parameter, Constant, value, num_of_param

include("parameters.jl")
include("idset.jl")
include("parameter_count.jl")

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Elfutils"
version = v"0.182"

# Collection of sources required to build Elfutils
sources = [
    ArchiveSource("https://sourceware.org/elfutils/ftp/$(version.major).$(version.minor)/elfutils-$(version.major).$(version.minor).tar.bz2",
                  "ecc406914edf335f0b7fc084ebe6c460c4d6d5175bfdd6688c1c78d9146b8858"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/elfutils-*/
if [[ ${target} = *-musl* ]] ; then
    for patchfile in $WORKSPACE/srcdir/patches/*; do
        atomic_patch -p1 $patchfile
    done
    cp $WORKSPACE/srcdir/error.h src/
    cp $WORKSPACE/srcdir/error.h lib/

    apk add bsd-compat-headers
    # /usr/include isn't in search path of cross-cc, so copy cdefs.h
    mkdir -p $prefix/include/sys
    # Skip warning macro at top of file
    tail -n +2 /usr/include/sys/cdefs.h >$prefix/include/sys/cdefs.h
    autoreconf -vif
fi
export CC=gcc
export CXX=g++
CFLAGS="-Wno-error=unused-result" CPPFLAGS="-I${prefix}/include" ./configure \
    --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-debuginfod \
    --disable-libdebuginfod
make -j${nproc}
make install
rm -f "${includedir}/sys/cdefs.h"
install_license COPYING*
"""

# Only build for Linux
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libasm", :libasm),
    LibraryProduct("libdw", :libdw),
    LibraryProduct("libelf", :libelf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    # Future versions of bzip2 should allow a more relaxed compat because the
    # soname of the macOS library shouldn't change at every patch release.
    Dependency("Bzip2_jll"; compat="1.0.8"),
    Dependency("XZ_jll"),
    Dependency("argp_standalone_jll"),
    Dependency("fts_jll"),
    Dependency("obstack_jll"; compat="~1.2.2"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
InputFolder = './Images/TissueSimpleMovie/Movie1/';
OutputFolder = './Results/Images/TissueSimpleMovie/Movie1/';

@iA = '*.tif';

@fxg_mWaterTiles [iA] > [M];
params.GRad = 4;
params.ExtendedMinThr = 1;
/endf

/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
include("distribution.jl")
include("const.jl")

# data[w', w, a] -> |MI(W';W|A)|
function MC_W_CMI(data::Matrix{Int64}, w_cardinality::Int64; base=2)
  @assert (size(data)[2] == 2) "MC_W required two columns (W,A)."
  w2      = data[2:end,1]
  w1      = data[1:end-1,1]
  a1      = data[1:end-1,2]
  pw2w1a1 = fe3p(hcat(w2,w1,a1)) # add possibility to use other estimators
  return CMI(pw2w1a1, base=base) / log(base, w_cardinality) # |MI(W';W|A)|
end

# data[w', w, a] -> 1 - |MI(W';A|W)|
function MC_A_CMI(data::Matrix{Int64}, w_cardinality::Int64; base=2)
  @assert (size(data)[2] == 2) "MC_W required two columns (W,A)."
  w2      = data[2:end,1]
  w1      = data[1:end-1,1]
  a1      = data[1:end-1,2]
  pw2a1w1 = fe3p(hcat(w2,a1,w1)) # add possibility to use other estimators
  return 1.0 - CMI(pw2a1w1, base=base)/log(base, w_cardinality) # 1.0 - |MI(W';A|W)|
end

# data[w', w, a] -> |MI(W';W|A)|
function MC_W(data::Matrix{Int64}, w_cardinality::Int64; base=2)
  @assert (size(data)[2] == 2) "MC_W required two columns (W,A)."

  w2         = data[2:end,1]
  w1         = data[1:end-1,1]
  a1         = data[1:end-1,2]

  w2_max     = maximum(w2);
  w1_max     = maximum(w1);
  a1_max     = maximum(a1);

  pw2w1a1    = fe3p(hcat(w2,w1,a1)) # add possibility to use other estimators
  pw1a1      = sum(pw2w1a1, 1)
  pw2a1      = sum(pw2w1a1, 2)
  pa1        = sum(pw1a1,   2)

  pw2_c_a1   = zeros(w2_max, a1_max);
  pw2_c_w1a1 = zeros(w2_max, w1_max, a1_max);

  for w2_index = 1:w2_max
    for a1_index = 1:a1_max
      if pa1[1,1,a1_index] > 0.0
        pw2_c_a1[w2_index, a1_index] = pw2a1[w2_index, 1, a1_index] / pa1[1, 1, a1_index]
      end
    end
  end

  for w2_index = 1:w2_max
    for w1_index = 1:w1_max
      for a1_index = 1:a1_max
        pw2_c_w1a1[w2_index, w1_index, a1_index] = 0.0
        if pw2a1[w2_index, 1, a1_index] > 0.0
          pw2_c_w1a1[w2_index, w1_index, a1_index] = pw2w1a1[w2_index, w1_index, a1_index] / pw2a1[w2_index, 1, a1_index]
        end
      end
    end
  end

  r = 0.0
  for w2_index = 1:w2_max
    for w1_index = 1:w1_max
      for a1_index = 1:a1_max
        if pw2w1a1[w2_index, w1_index, a1_index]    > 0.0 &&
           pw2_c_w1a1[w2_index, w1_index, a1_index] > 0.0 &&
           pw2_c_a1[w2_index, a1_index]          > 0.0
          r = r + pw2w1a1[w2_index, w1_index, a1_index] * (log(base, pw2_c_w1a1[w2_index, w1_index, a1_index]) - log(base, pw2_c_a1[w2_index, a1_index]))
        end
      end
    end
  end

  return r / log(base, w_cardinality) # |MI(W';W|A)|
end

# data[w', w, a] -> 1 - |MI(W';A|W)|
# D(p(w'|w,a) || p(w'|w))
function MC_A(data::Matrix{Int64}, w_cardinality::Int64; base=2)
  @assert (size(data)[2] == 2) "MC_W required two columns (W,A)."

  w2         = data[2:end,1]
  w1         = data[1:end-1,1]
  a1         = data[1:end-1,2]

  w2_max     = maximum(w2);
  w1_max     = maximum(w1);
  a1_max     = maximum(a1);

  pw2w1a1    = fe3p(hcat(w2,w1,a1)) # add possibility to use other estimators
  pw1a1      = sum(pw2w1a1, 1)
  pw2w1      = sum(pw2w1a1, 3)
  pw1        = sum(pw1a1,   3)

  pw2_c_w1   = zeros(w2_max, w1_max);
  pw2_c_w1a1 = zeros(w2_max, w1_max, a1_max);

  for w2_index = 1:w2_max
    for w1_index = 1:a1_max
      if pw1[1,w1_index,1] > 0.0
        pw2_c_w1[w2_index, w1_index] = pw2w1[w2_index, w1_index, 1] / pw1[1, w1_index, 1]
      end
    end
  end

  for w2_index = 1:w2_max
    for w1_index = 1:w1_max
      for a1_index = 1:a1_max
        if pw1a1[1, w1_index, a1_index] > 0.0
          pw2_c_w1a1[w2_index, w1_index, a1_index] = pw2w1a1[w2_index, w1_index, a1_index] / pw1a1[1, w1_index, a1_index]
        end
      end
    end
  end

  r = 0.0
  for w2_index = 1:w2_max
    for w1_index = 1:w1_max
      for a1_index = 1:a1_max
        if pw2w1a1[w2_index, w1_index, a1_index]   > 0.0 &&
          pw2_c_w1a1[w2_index, w1_index, a1_index] > 0.0 &&
          pw2_c_w1[w2_index, w1_index, 1]          > 0.0
          r = r + pw2w1a1[w2_index, w1_index, a1_index] * (log(base, pw2_c_w1a1[w2_index, w1_index, a1_index]) - log(base, pw2_c_w1[w2_index, w1_index]))
        end
      end
    end
  end

  return 1.0 - r/log(base, w_cardinality) # 1.0 - |MI(W';A|W)|
end

# CA = ∑ p(s,a) p(s'|do(a)) log (p(s'|do(a)) / p(s'|do(s)))
function C_A(data::Matrix{Int64}, s_cardinality::Int64; base=2)
  @assert (size(data)[2] == 2) "C_A required two columns (S,A)."

  s_max   = maximum(data[:,1])
  a_max   = maximum(data[:,2])

  s2      = data[2:end,1]
  s1      = data[1:end-1,1]
  a1      = data[1:end-1,2]

  ps2s1a1 = fe3p(hcat(s2,s1,a1)) # add possibility to use other estimators

  ps1a1   = sum(ps2s1a1, 1)
  pa1     = sum(ps2s1a1, (1,2))
  ps1     = sum(ps2s1a1, (1,3))

  # p(s'|s,a) = p(s',s,a) / p(s,a)
  ps2_c_s1a1 = zeros(size(ps2s1a1))
  for s2_index = 1:s_max
    for s1_index = 1:s_max
      for a1_index = 1:a_max
        if abs(ps1a1[1, s1_index, a1_index]) > ϵ
          ps2_c_s1a1[s2_index, s1_index, a1_index] = ps2s1a1[s2_index, s1_index, a1_index] / ps1a1[1, s1_index, a1_index]
        end
      end
    end
  end

  # p(s'| do(a)) = ∑ p(s'|s,a) * p(s)
  ps2_do_a1 = zeros(s_max, a_max)
  for s2_index = 1:s_max
    for a1_index = 1:a_max
      for s1_index = 1:s_max
        ps2_do_a1[s2_index, a1_index] = ps2_do_a1[s2_index, a1_index] + ps2_c_s1a1[s2_index, s1_index, a1_index] * ps1[1,s1_index, 1]
      end
    end
  end

  # p(a|s) = p(s,a) / p(s)
  pa1_c_s1 = zeros(a_max, s_max)
  for a1_index = 1:a_max
    for s1_index = 1:s_max
      if abs(ps1[s1_index]) > ϵ
        pa1_c_s1[a1_index, s1_index] = ps1a1[1, s1_index, a1_index] / ps1[1, s1_index, 1]
      end
    end
  end

  # p(s'|do(s)) = ∑ p(a|s) p(s'|do(a))
  ps2_do_s1 = zeros(s_max, s_max)
  for s2_index = 1:s_max
    for s1_index = 1:s_max
      for a1_index = 1:a_max
        ps2_do_s1[s2_index, s1_index] = ps2_do_s1[s2_index, s1_index] + pa1_c_s1[a1_index, s1_index] * ps2_do_a1[s2_index, a1_index]
      end
    end
  end

  # CA = ∑ p(s,a) p(s'|do(a)) log (p(s'|do(a)) / p(s'|do(s)))
  r = 0.0
  for s2_index = 1:s_max
    for s1_index = 1:s_max
      for a1_index = 1:a_max
        if abs(ps1a1[1, s1_index, a1_index]) > ϵ && abs(ps2_do_a1[s2_index, a1_index]) > ϵ && abs(ps2_do_s1[s2_index, s1_index]) > ϵ
          r = r + ps1a1[1, s1_index, a1_index] * ps2_do_a1[s2_index, a1_index] * (log(base, ps2_do_a1[s2_index, a1_index]) - log(base, ps2_do_s1[s2_index, s1_index]))
        end
      end
    end
  end
  return r / log(base, s_cardinality)
end

# data[s', s, a] -> D_KL(p(s'|s,a)||p(s'|a))
function C_W(data::Matrix{Int64}, s_cardinality::Int64; base=2)
  @assert (size(data)[2] == 2) "C_A required two columns (S,A)."
  s_max   = maximum(data[:,1])
  a_max   = maximum(data[:,2])
  s2      = data[2:end,1]
  s1      = data[1:end-1,1]
  a1      = data[1:end-1,2]
  ps2s1a1 = fe3p(hcat(s2,s1,a1)) # add possibility to use other estimators

  ps1     = sum(ps2s1a1, (1,3))
  pa1     = sum(ps2s1a1, (1,2))
  ps2s1   = sum(ps2s1a1, 3)
  ps1a1   = sum(ps2s1a1, 1)
  ps2a1   = sum(ps2s1a1, 2)

  # p(s'|s) = p(s',s) / p(s)
  ps2_c_s1 = zeros(s_max, s_max)
  for s2_index = 1:s_max
    for s1_index = 1:s_max
      if abs(ps1[s1_index]) > ϵ
        ps2_c_s1[s2_index, s1_index] = ps2s1[s2_index, s1_index] / ps1[s1_index]
      end
    end
  end

  # p(s'|a) = p(s',a) / p(a)
  ps2_c_a1 = zeros(s_max, a_max)
  for s2_index = 1:s_max
    for a1_index = 1:a_max
      if abs(pa1[a1_index]) > ϵ
        ps2_c_a1[s2_index, a1_index] = ps2a1[s2_index, 1, a1_index] / pa1[a1_index]
      end
    end
  end

  # p(a|s) = p(a,s) / p(s)
  pa1_c_s1 = zeros(a_max, s_max)
  for a1_index = 1:a_max
    for s1_index = 1:s_max
      if abs(ps1[1, s1_index, 1]) > ϵ
        pa1_c_s1[a1_index, s1_index] = ps1a1[1, s1_index, a1_index] / ps1[1, s1_index, 1]
      end
    end
  end

  #̂̂p^(s'|s) = ∑ p(s'|a) p(a|s)
  phat_s2_c_s1 = zeros(s_max, s_max)
  for s2_index = 1:s_max
    for s1_index = 1:s_max
      for a1_index = 1:a_max
        phat_s2_c_s1[s2_index, s1_index] = ps2_c_a1[s2_index, a1_index] * pa1_c_s1[a1_index, s1_index]
      end
    end
  end

  # C_W = ∑ p(s',s) log p(s'|s) / phat(s'|s)
  r = 0
  for s1_index = 1:s_max
    for s2_index = 1:s_max
      if abs(ps2s1[s2_index, s1_index, 1]) > ϵ && abs(ps2_c_s1[s2_index, s1_index]) > ϵ && abs(phat_s2_c_s1[s2_index, s1_index]) > ϵ
        r = r + ps2s1[s2_index, s1_index, 1] * (log(base, ps2_c_s1[s2_index, s1_index]) - log(base, phat_s2_c_s1[s2_index, s1_index]))
      end
    end
  end

  return r / log(base, s_cardinality)
end

# Morphological Computation
function quantify(data::Matrix{Int64}; base=2, mode="MC A", pseudocount=0)
  modes      = ["MC A",
  "MC W",
  "ASOC A",
  "ASOC W",
  "C A",
  "C W"
  ]
  umodes     = map(x->uppercase(x), modes)
  known_mode = uppercase(mode) in umodes
  pmodes     = map(x->", $x",modes)
  pmodes     = foldl(*, pmodes[1][3:end], pmodes[2:end])
  @assert known_mode "Mode may be any of the following: [$pmodes]."

  if uppercase(mode) in umodes[1:2]
    return MC_A(data, base)
  elseif uppercase(mode) in umodes[3:4]
    return MC_W(data, base)
  elseif uppercase(mode) in umodes[5]
    return MC_A(data, base)
  elseif uppercase(mode) in umodes[6]
    return MC_W(data, base)
  elseif uppercase(mode) in umodes[7]
    return C_A(data, base)
  elseif uppercase(mode) in umodes[8]
    return C_W(data, base)
  end
end


text:
function SeisSmooth1(A;L=7)
	a = isodd(L) ? 0 : 1
	kernel = ones(eltype(A),L)/L;
	A = conv(A,kernel)
	A = A[floor(Int,L/2) + 1:end-floor(Int,L/2) + a]
	return A
end

function SeisSmooth1(A,h::Array{Header,1};L=7)
	for itrace = 1 : size(A[:,:],2)
		A[:,itrace] = SeisSmooth1(A[:,itrace],L=L)
	end
	return A,h
end

function SeisSmooth1(in::AbstractString,out::AbstractString;L=7)
	@compat parameters = Dict(:L=>L)
	SeisProcess(in,out,[SeisSmooth1],[parameters];group="some",ntrace=100000)
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
mutable struct EconometricsMLJModel{M, C, F, L} <: Probabilistic
    contrasts::C
    formula::F
    label::L
end

function EconometricsMLJModel(;model::Type{M},
                               contrasts::C = nothing,
                               formula::F,
                               label::L) where M where C where F where L
    return EconometricsMLJModel{M, C, F, L}(contrasts, formula, label)
end

function _merge_X_and_y(X, # X must be a table
                        y::AbstractVector,
                        label)
    return merge(columntable(X), NamedTuple{(label,)}((y,)))
end

function MLJModelInterface.fit(model::EconometricsMLJModel{M},
                               verbosity::Integer,
                               X, # X must be a table
                               y::AbstractVector,
                               w = nothing,
                               ) where M
    contrasts = model.contrasts
    formula = model.formula
    label = model.label
    data = _merge_X_and_y(X, y, model.label)
    if contrasts === nothing
        fitresult = fit(M, formula, data)
    else
        fitresult = fit(M, formula, data; contrasts = contrasts)
    end
    cache = nothing
    report = NamedTuple{}()
    return (fitresult, cache, report)
end

function clean!(model::EconometricsMLJModel)
    warning = ""
    return warning
end

# TODO: implement the `MLJModelInterface.predict` method
function MLJModelInterface.predict(model::EconometricsMLJModel,
                                   fitresult,
                                   Xnew, # Xnew must be a table
                                   )
    throw(ErrorException("The predict method has not yet been implemented."))
end


text:
using Bolt
using ForwardDiff
using PyPlot
using BenchmarkTools

# Cₗ as a function of baryon density
function clb(Ω_b::DT, ells) where DT
    par = CosmoParams{DT}(Ω_b=Ω_b)
    bg = Background(par)
    ih = IonizationHistory(Peebles(), par, bg)
    k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 500)
    sf = source_grid(par, bg, ih, k_grid, BasicNewtonian())
    return cltt(ells, par, bg, ih, sf)
end

ells = 10:2:1200
f(Ω_b) = clb(Ω_b, ells)
@time cl = f(0.046)
@time ∂cl = ForwardDiff.derivative(f, 0.046)  # you can just ForwardDiff the whole thing

##
clf()
plt.figure(figsize=(10,5))
plot(ells, cl .* ells.^2, "-", label=raw"$C_{\ell}$")
plot(ells, ∂cl .* ells.^2 / 10, "-",
    label=raw"$\partial C_{\ell}/\partial\Omega_b / 10$")

# uncomment to see finite differences
# Δ = 1e-3
# @time finitediff_∂cl = (f(0.046 + Δ) .- f(0.046 - Δ)) ./ 2Δ
# plt.plot(ells, finitediff_∂cl .* ells.^2 / 10, "--",
#     label=raw"$\Delta C_{\ell}/\Delta\Omega_b / 10$")

ylabel(raw"$\ell^2 C_{\ell}^{TT}$")
xlabel(raw"$\ell$")
legend()
ylim(-0.3, 0.5)
tight_layout()
# savefig("docs/assets/example_spectrum.png")
gcf()

##
sf

##


text:
# Variational approximations
# (Excepted for loglikelihood(Normal, Normal))

"""
Expectation of the log-likelihood of a gaussian wrt. another gaussian
(Penny and Roberts, 2000)
E_{a} L(b)
"""
function loglikelihood(a::Normal, b::Normal)
    if iszero(b.σ)
        (a.μ == b.μ) ? Inf : 0
    else
        -(1/2) * (log2π + log(b.σ^2) + (a.σ^2 / b.σ^2) + ((b.μ - a.μ)^2 / b.σ^2))
    end
end

function loglikelihood_va(a::MixtureModel, b::MixtureModel)
    I, J = ncomponents(a), ncomponents(b)

    # Optimal lower-bound to the expected log-likelihood
    lb = 0.0

    # Expected log-likelihoods
    logl = Matrix{Float64}(undef, (I, J))

    # η[i,j]: probability that an observation from
    # component i of a corresponds to component j of b
    logη = Matrix{Float64}(undef, (I, J))

    # Type annotation is important for performane
    # otherwiser the compiler infers `Any`.
    logp::Vector{Float64} = log.(b.prior.p::Vector{Float64})

    for (i, ai) in enumerate(a.components)
        for (j, bj) in enumerate(b.components)
            logl[i,j] = loglikelihood(ai, bj)
            logη[i,j] = logp[j] + logl[i,j]
        end
        logη[i,:] .-= logsumexp(logη[i,:])
    end

    for i in OneTo(I)
        p::Float64 = (a.prior.p::Vector{Float64})[i]
        for j in OneTo(J)
            lb += p * exp(logη[i,j]) * (logp[j] - logη[i,j] + logl[i,j])
        end
    end

    lb, logη
end

"""
Variational lower bound to the expected log-likelihood.

Coviello 14, Hershey 07.
"""
function loglikelihood_va(hmm1::HMM, hmm2::HMM, lgmm::Matrix{Float64}, τ::Integer)
    @argcheck size(hmm1, 1) == size(lgmm, 1)
    @argcheck size(hmm2, 1) == size(lgmm, 2)
    @argcheck τ >= 0

    K, L = size(hmm1, 1), size(hmm2, 1)

    # Optimal lower-bound to the expected log-likelihood
    lhmm = 0.0

    # logl[t,βp,ρp]
    logl = zeros(τ+1, K, L)

    # logϕ[t,β,ρp,ρ]
    logϕ = zeros(τ, K, L, L)

    # logϕ1[β,ρ]
    # TODO: Merge with logϕ (ignore previous state for t = 1) ?
    logϕ1 = zeros(K, L)

    # Initial probabilities (π in the paper)
    # Transition matrices (a or A in the paper)
    loga2 = log.(hmm2.a)
    logA2 = log.(hmm2.A)

    # Appendix B. (p. 737)
    # β  = β_t
    # βp = β_{t-1}

    for t in τ:-1:2
        for β in 1:K
            for ρp in 1:L
                for ρ in 1:L
                    logϕ[t,β,ρp,ρ] = logA2[ρp,ρ] + lgmm[β,ρ] + logl[t+1,β,ρ]
                end

                norm = logsumexp(logϕ[t,β,ρp,:])
                logϕ[t,β,ρp,:] .-= norm

                for βp in 1:K
                    logl[t,βp,ρp] += hmm1.A[βp,β] * norm
                end
            end
        end
    end

    for β in 1:K
        for ρ in 1:L
            logϕ1[β,ρ] = loga2[ρ] + lgmm[β,ρ] + logl[2,β,ρ]
        end

        norm = logsumexp(logϕ1[β,:])
        logϕ1[β,:] .-= norm

        lhmm += hmm1.a[β] * norm
    end

    lhmm, logϕ, logϕ1
end

function loglikelihood_va(hmm1::HMM, hmm2::HMM, τ::Integer)
    Ki, Kj = size(hmm1, 1), size(hmm2, 1)

    logη = Matrix{Matrix{Float64}}(undef, Ki, Kj)
    lgmm = zeros(Ki, Kj)

    for (β, Miβ) in enumerate(hmm1.B)
        for (ρ, Mjρ) in enumerate(hmm2.B)
            lgmm[β,ρ], logη[β,ρ] = loglikelihood_va(Miβ, Mjρ)
        end
    end

    logη, lgmm, loglikelihood_va(hmm1, hmm2, lgmm, τ)
end


text:
typealias  Graph Dict{Int, Set{Int}}

function addnode!(g::Graph, u::Int)
    if !haskey(g, u)
        g[u] = valtype(g)()
    end
    nothing
end

function delnode!(g::Graph, u::Int)
    if haskey(g, u)
        for v in g[u]
            delete!(g[v], u)
        end
        delete!(g, u)
    end
    nothing
end

function addedge!(g::Graph, u::Int, v::Int)
    addnode!(g, u)
    addnode!(g, v)
    push!(g[u], v)
    push!(g[v], u)
    nothing
end

function deledge!(g::Graph, u::Int, v::Int)
    if haskey(g, u) && haskey(g, v)
        delete!(g[u], v)
        delete!(g[v], u)
    end
    nothing
end

function ne(g::Graph)
    m = 0
    for u in keys(g)
        m += length(g[u])
    end
    return div(m, 2)
end

neighbors(g::Graph, u::Int) = g[u]

nodes(g::Graph) = keys(g)

hasnode(g::Graph, u::Int) = haskey(g, u)

hasedge(g::Graph, u::Int, v::Int) = hasnode(g, u) && in(v, g[u])


text:
using Base.Test
using JSON
using NMR

addprocs(Sys.CPU_CORES)
@everywhere using NMR

include("test_numberofflights.jl")
include("test_listofflights.jl")
include("test_numpassengers.jl")
include("test_lineofsight.jl")


text:
abstract type AbstractStack{T} end

### Interface

# Push an item onto the top of the stack.
push!(s::AbstractStack, item) = nothing

# Remove and return top element of the stack.
pop!(s::AbstractStack) = nothing

# Return the top element.
topvalue(s::AbstractStack) = nothing

# Return the number of elements in the stack.
length(s::AbstractStack) = nothing

# Check if stack is empty.
isempty(s::AbstractStack) = nothing

### Generic

eltype(s::AbstractStack{T}) where T = T


text:
using PyCall
using FromPython: AST
using FromPython

ast = pyimport("ast")

@testset "BinOp" begin
    for (op, x, y, jl_expr) in [
        (ast.Add(), 1, 2, :(1 + 2)),
        (ast.Sub(), 1, 2, :(1 - 2)),
        (ast.Mult(), 4, 5, :(4 * 5)),
        (ast.Div(), 4, 2, :(4 / 2)),

        (ast.Mod(), 15, 4, :(15 % 4)),
        (ast.Pow(), 2, 3, :(2 ^ 3)),
        (ast.LShift(), 2, 1, :(2 << 1)),
        (ast.RShift(), 2, 1, :(2 >> 1)),
        (ast.BitOr(),  4, 1, :(4 | 1)),
        # (ast.BitXor(), 4, 1, :(xor(4, 1)),
        (ast.BitAnd(), 4, 1, :(4 & 1)),

        # (ast.FloorDiv(),=> floor ∘ (/)
        # (ast.MatMult(), => @not_implemented_yet
    ]
        binop = ast.BinOp(op=op, left=ast.Num(x), right=ast.Num(y))
        @test AST.dispatch(Expr, binop) == jl_expr
    end
end

# @testset "Comparisons" begin
#     (op, x, y, jl_expr) in [
#     # (ast.Eq(), => (==),
#     # (ast.NotEq(), => (!=),
#     # (ast.Lt(), => (<),
#     # (ast.LtE(), => (<=),
#     # (ast.Gt(),  => (>),
#     # (ast.GtE(), => (>=),
#     # (ast.Is(),  => (===),
#     # (ast.IsNot(), => (!==),
#     # (ast.In(),    => (in),
#     #  (ast.NotIn(), => (!in),
#       ]
# end

# @testset "Assign" begin
#     assgn = ast.parse("x = 1").body[1]
#     tup_assgn = ast.parse("y, z = 2, 3").body[1]
#     @test AST.dispatch(Expr, assgn) == :(x = 1)
#     @test AST.dispatch(Expr, tup_assgn) == :((y, z) = (2, 3))
# end

@testset "Expr" begin
    binop = ast.BinOp(op=ast.Add, left=ast.Num(1), right=ast.Num(2))
    pyexpr = ast.Expr(value=binop)
    @test AST.dispatch(Expr, pyexpr) == :(1 + 2)
end

@testset "Module" begin
    binop = ast.BinOp(op=ast.Add, left=ast.Num(1), right=ast.Num(2))
    pyexpr = ast.Expr(value=binop)
    modl = ast.Module(body=[pyexpr])
    @test AST.dispatch(Expr, modl) == [:(1 + 2)]
end

# @testset "Call" begin
#     call_exp = ast.parse("numpy.dot(x, y)").body[1].value
#     @test AST.dispatch(Expr, call_exp) == [:(numpy.dot(x, y))]

#     FromPython.register!("numpy.dot",
#                             (ctx, func, args, keywords) -> :(dot($(AST.dispatch(ctx, args[1])),
#                                                                  $(AST.dispatch(ctx, args[2])))))
#     call_exp = ast.parse("numpy.dot(x, y)").body[1].value
#     @test AST.dispatch(Expr, call_exp) == [:(dot(x, y))]
# end


text:
#!/usr/bin/env julia

"""
Tests for the [`Sequences`](@ref) common module.
"""
module TestSequences

include(joinpath("..", "src", "Sequences.jl"))

using Test

@testset "Sequences" begin
    @testset "Fibonacci" begin
        @test first(Sequences.Fibonacci()) == 1
        @test first(Sequences.Fibonacci(UInt8)) == 1
        @test typeof(first(Sequences.Fibonacci(UInt8))) == UInt8
        @test first(Sequences.Fibonacci(BigInt)) == 1
        @test typeof(first(Sequences.Fibonacci(BigInt))) == BigInt

        fib_seq = [1, 2, 3, 5, 8, 13]
        for (i, n) in enumerate(Sequences.Fibonacci())
            if i > length(fib_seq)
                break
            end
            @test n == fib_seq[i]
        end

        for (i, n) in enumerate(Sequences.Fibonacci(BigInt))
            if i <= length(fib_seq)
                @test n == fib_seq[i]
            elseif i == 13
                @test n == 377
            elseif i == 37
                @test n == 39088169
            elseif i == 170
                @test (
                    n == parse(BigInt, "244006547798191185585064349218729154")
                )
            elseif i > 170
                break
            end
        end
    end

    @testset "arithmetic_series" begin
        @test Sequences.arithmetic_series(1, 4) == 10
        @test Sequences.arithmetic_series(2, 6, -1) == -3
        @test Sequences.arithmetic_series(3, 5, 4) == 55
        @test Sequences.arithmetic_series(4, 1) == 4
        @test Sequences.arithmetic_series(4, 1, 1) == 4
        @test Sequences.arithmetic_series(4, 1, 2) == 4
        @test Sequences.arithmetic_series(1, 2, -1) == 1
        @test Sequences.arithmetic_series(19, 5, 4) == 135
        @test Sequences.arithmetic_series(7, 8, -3) == -28
        @test Sequences.arithmetic_series(-14, 8, 3) == -28
    end

    @testset "fibonacci" begin
        @test Sequences.fibonacci(0) == 1
        @test Sequences.fibonacci(1) == 1
        @test Sequences.fibonacci(2) == 2
        @test Sequences.fibonacci(3) == 3
        @test Sequences.fibonacci(4) == 5
        @test Sequences.fibonacci(5) == 8
        @test Sequences.fibonacci(6) == 13
        @test Sequences.fibonacci(13) == 377
        @test Sequences.fibonacci(37) == 39088169
        @test (
            Sequences.fibonacci(BigInt, 200)
                == parse(BigInt, "453973694165307953197296969697410619233826")
        )
    end

    @testset "sum_of_squares" begin
        @test Sequences.sum_of_squares(0) == 0
        @test Sequences.sum_of_squares(1) == 1
        @test Sequences.sum_of_squares(2) == 5
        @test Sequences.sum_of_squares(3) == 14
        @test Sequences.sum_of_squares(4) == 30
        @test Sequences.sum_of_squares(5) == 55
        @test Sequences.sum_of_squares(6) == 91
        @test Sequences.sum_of_squares(7) == 140
        @test Sequences.sum_of_squares(401) == 21574201
        @test Sequences.sum_of_squares(34594) == 13800661997045
    end

    @testset "triangular" begin
        tri_nums = [1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120]
        for (i, n) in enumerate(tri_nums)
            @test Sequences.triangular(i) == n
        end
        @test Sequences.triangular(292) == 42778
        @test Sequences.triangular(38483) == 740489886
        @test Sequences.triangular(4946666) == 12234754731111
    end
end

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
# Abstract type for Global Spatial Autocorrelation
abstract type AbstractGlobalSpatialAutocorrelation end

function Base.show(io::IO, x::AbstractGlobalSpatialAutocorrelation)

    println(io, testname(x) * " test of Global Spatial Autocorrelation")
    println(io, "--------------------------------------------")
    println(io, "")
    println(io, testname(x) * ": ", round(score(x), digits = 7))
    println(io, "Expectation: ", round(expected(x), digits = 7))
    println(io, "")
    println(io, "Randomization test with ", length(scoreperms(x)), " permutations.")
    println(io, " Mean: ", round(mean(x), digits = 7))
    println(io, " Std Error: ", round(std(x), digits = 7))
    println(io, " zscore: ", round(zscore(x), digits = 7))
    println(io, " p-value: ", round(pvalue(x), digits = 7))

end


text:
function permanova(data)
    X = data[:X]
    Y = data[:Y]
    Z = data[:Z]
    R"""
    suppressMessages(suppressWarnings(library(vegan)))
    z <- data.frame($Z)
    z[] <- lapply(z, as.factor)
    x <- data.frame($X)
    x[] <- lapply(x, as.factor)
    sim <- adonis($Y ~ $X + $Z, method="bray", nperm=9999)
    print(sim)
    pval <- sim$aov.tab$`Pr(>F)`[1]
    """
    @rget pval
    return Dict{Symbol, Any}(:pval => pval)
end

max_stars_repo_name: reich-group/MIMIX

max_issues_repo_name: reich-group/MIMIX

max_forks_repo_name: reich-group/MIMIX


text:
const data      = "../data"
const simulated = joinpath(data, "simulated")
const real      = joinpath(data, "real")
const realbam   = joinpath(real,"bam")
const simbam    = joinpath(simulated,"bam")
const ref       = joinpath(data, "ref")

const varfile = joinpath(data,"variants.bed")
const inbam = joinpath(realbam, "069_gdna_sort.bam")
const hg19 = joinpath(ref, "ucsc.hg19.fasta")
const outbam = joinpath(simbam, "069_gdna_sort_sim.bam")

const snvfrac = 1
const mutfrac = 0.5
const numsnv  = 10
const coveragefile = joinpath(data,"coverage.bed") #todo
const haplosize = 2
const procs = 50
const picardjar = "/home/guo/haplxo/Github/picard-tools-1.131/picard.jar"
const mindepth = 5
const maxdepth = 2000
const minmutreads = 4 #important threshold
#const avoidreads = 100
const maxopen = 1024
const aligner = "bwa"
const alignopts = ""

const varfile1 = joinpath(data,"variants_1.bed")
const outbam1 = joinpath(simbam, "069_gdna_sort_sim1.bam")

cmd1 = `addsnv.py -v $varfile1 -f $inbam -r $hg19 -o $outbam1
                 --tmpdir "/haplox/tmp" `
#=
cmd = `addsnv.py -v $varfile -f $inbam -r $hg19 -o $outbam -s $snvfrac
                 -m $mutfrac -n $numsnv -d $coveragefile -z $haplosize -p $procs
                 --picardjar $picardjar --mindepth $mindepth --maxdepth $maxdepth
                 --minmutreads $minmutreads --avoidreads $avoidreads --nomut
                 --ignoresnps  --maxopen $maxopen --requirepaired
                 --skipmerge --ignorepileup --aligner $aligner --alignopts $alignopts
                 --tmpdir "/haplox/tmp" `
=#
println("comand for addsnv.py is:")
println(cmd1)

run(cmd1)


text:
import Markdown

"""
    look(s[, k])

Look up information about system or variable. Both system type `S` and instance `s` are accepted. For looking up a variable, the name of variable `k` needs to be specified in a symbol.

See also: [`@look`](@ref), [`dive`](@ref)

# Arguments
- `s::Union{System,Type{<:System}}`: target system.
- `k::Symbol`: name of variable.

# Examples
```julia-repl
julia> "my system"
       @system S(Controller) begin
           "a param"
           a => 1 ~ preserve(parameter)
       end;

julia> s = instance(S);

julia> look(s)
[doc]
  my system

[system]
S
  context = <Context>
  config = <Config>
  a = 1.0
julia> look(s, :a)
[doc]
  a param

[code]
  a => 1 ~ preserve(parameter)

[value]
1.0
```

"""
look(s::System; kw...) = look(stdout, s; kw...)
look(S::Type{<:System}; kw...) = look(stdout, S; kw...)
look(s::System, k::Symbol; kw...) = look(stdout, s, k; kw...)
look(S::Type{<:System}, k::Symbol; kw...) = look(stdout, S, k; kw...)
look(m::Module, s::Symbol; kw...) = look(getfield(m, s); kw...)

look(io::IO, s::Union{S,Type{S}}; header=true, doc=true, system=true, kw...) where {S<:System} = begin
    print(io, join(filter(!isempty, strip.([
        doc ? buf2str(io -> lookdoc(io, s; header, kw...)) : "",
        system ? buf2str(io -> looksystem(io, s; header, kw...)) : "",
    ])), "\n\n"))
end
look(io::IO, S::Type{<:System}, k::Symbol; header=true, doc=true, code=true, kw...) = begin
    print(io, join(filter(!isempty, strip.([
        doc ? buf2str(io -> lookdoc(io, S, k; header, kw...)) : "",
        code ? buf2str(io -> lookcode(io, S, k; header, kw...)) : "",
    ])), "\n\n"))
end
look(io::IO, s::S, k::Symbol; header=true, value=true, kw...) where {S<:System} = begin
    print(io, join(filter(!isempty, strip.([
        buf2str(io -> look(io, S, k; header, kw...)),
        value ? buf2str(io -> lookvalue(io, s, k; header, kw...)) : "",
    ])), "\n\n"))
end

lookdoc(io::IO, ::Union{S,Type{S}}; header=false, kw...) where {S<:System} = begin
    header && printstyled(io, "[doc]\n", color=:light_black)
    try
        #HACK: mimic REPL.doc(b) with no dynamic concatenation
        md = Docs.formatdoc(fetchdocstr(S))
        show(io, MIME("text/plain"), md)
    catch
    end
end
looksystem(io::IO, s::Union{S,Type{S}}; header=false, kw...) where {S<:System} = begin
    header && printstyled(io, "[system]\n", color=:light_black)
    printstyled(io, namefor(S), color=:light_magenta)
    for (n, a) in fieldnamesalias(S)
        print(io, "\n  ")
        printstyled(io, uncanonicalname(n), color=:light_blue)
        !isnothing(a) && printstyled(io, " (", uncanonicalname(a), ")", color=:light_black)
        s isa Type && continue
        printstyled(io, " = ", color=:light_black)
        print(io, labelstring(s[n]))
    end
end

lookdoc(io::IO, ::Union{S,Type{S}}, k::Symbol; header=false, excerpt=false, kw...) where {S<:System} = begin
    header && printstyled(io, "[doc]\n", color=:light_black)
    #HACK: mimic REPL.fielddoc(b, k) with no default description
    docstr = fetchdocstr(S)
    isnothing(docstr) && return
    n = canonicalname(k, S)
    ds = get(docstr.data[:fields], n, nothing)
    isnothing(ds) && return
    md = ds isa Markdown.MD ? ds : Markdown.parse(ds)
    s = if excerpt
        ts = Markdown.terminline_string(io, md)
        split(strip(ts), '\n')[1] |> Text
    else
        md
    end
    show(io, MIME("text/plain"), s)
end
lookcode(io::IO, ::Union{S,Type{S}}, k::Symbol; header=false, kw...) where {S<:System} = begin
    header && printstyled(io, "[code]\n", color=:light_black)
    d = dependency(S)
    n = canonicalname(k, S)
    v = d.M[n]
    Highlights.highlight(io, MIME("text/ansi"), "  " * string(v.line), Highlights.Lexers.JuliaLexer)
end
lookvalue(io::IO, s::System, k::Symbol; header=false, kw...) = begin
    header && printstyled(io, "[value]\n", color=:light_black)
    n = canonicalname(k, s)
    show(io, MIME("text/plain"), s[n])
end

buf2str(f; color=true, kw...) = begin
    d = Dict(:color => color, kw...)
    b = IOBuffer()
    io = IOContext(b, d...)
    f(io)
    String(take!(b))
end

using MacroTools: @capture

"""
    @look ex
    @look s[, k]

Macro version of `look` supports a convenient way of accessing variable without relying on symbol. Both `@look s.a` and `@look s a` work the same as `look(s, :a)`.

See also: [`look`](@ref)

# Examples
```julia-repl
julia> "my system"
       @system S(Controller) begin
           "a param"
           a => 1 ~ preserve(parameter)
       end;

julia> @look S.a
[doc]
  a param

[code]
  a => 1 ~ preserve(parameter)
```

"""
macro look(ex)
    if @capture(ex, s_.k_)
        :(Cropbox.look($(esc(s)), $(Meta.quot(k))))
    else
        :(Cropbox.look($(esc(ex))))
    end
end

macro look(s, k)
    :(Cropbox.look($(esc(s)), $(Meta.quot(k))))
end

fetchdocstr(S::Type{<:System}) = begin
    b = Docs.Binding(scopeof(S), nameof(typefor(S)))
    for m in Docs.modules
        d = Docs.meta(m)
        haskey(d, b) && return d[b].docs[Union{}]
    end
    nothing
end

getdoc(S::Type{<:System}) = begin
    docstr = fetchdocstr(S)
    isnothing(docstr) && return Markdown.parse("""
    No documentation found.

    Type `@look $(scopeof(S)).$(namefor(S))` for more information.
    """)
    b = IOBuffer()
    io = IOContext(b, :color => true)
    header = false
    doc = true
    look(io, S; header, doc, system=false)
    fields = docstr.data[:fields]
    if !isempty(fields)
        for (n, a) in fieldnamesalias(S)
            !haskey(fields, n) && continue
            println(io, '\n')
            entry = isnothing(a) ? "- `$n`" : "- `$n` (`$a`)"
            show(io, MIME("text/plain"), Markdown.parse(entry))
            print(io, ": ")
            look(io, S, n; header, doc, excerpt=true, code=false)
        end
    end
    String(take!(b)) |> Text
end

export look, @look


text:
type DistSparseMatrix{T} <: ElementalMatrix{T}
    obj::Ptr{Void}
end

for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:Complex64, :c),
                    (:Complex128, :z))
    @eval begin
        function DistSparseMatrix(::Type{$elty}, comm::ElComm = CommWorld)
            obj = Ref{Ptr{Void}}(C_NULL)
            err = ccall(($(string("ElDistSparseMatrixCreate_", ext)), libEl), Cuint,
                (Ref{Ptr{Void}}, ElComm),
                obj, comm)
            err == 0 || throw(ElError(err))
            return DistSparseMatrix{$elty}(obj[])
        end

        function comm(A::DistSparseMatrix{$elty})
            cm = Ref{ElComm}()
            err = ccall(($(string("ElDistSparseMatrixComm_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElComm}),
                A.obj, cm)
            err == 0 || throw(ElError(err))
            return cm[]
        end

        function destroy(A::DistSparseMatrix{$elty})
            err = ccall(($(string("ElDistSparseMatrixDestroy_", ext)), libEl), Cuint,
                (Ptr{Void},),
                A.obj)
            err == 0 || throw(ElError(err))
            return 0
        end

        function globalRow(A::DistSparseMatrix{$elty}, iLoc::Integer)
            i = Ref{ElInt}(0)
            err = ccall(($(string("ElDistSparseMatrixGlobalRow_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, Ref{ElInt}),
                A.obj, iLoc-1, i)
            err == 0 || throw(ElError(err))
            return i[]+1
        end

        function height(A::DistSparseMatrix{$elty})
            i = Ref{ElInt}(0)
            err = ccall(($(string("ElDistSparseMatrixHeight_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                A.obj, i)
            err == 0 || throw(ElError(err))
            return i[]
        end

        function localHeight(A::DistSparseMatrix{$elty})
            i = Ref{ElInt}(0)
            err = ccall(($(string("ElDistSparseMatrixLocalHeight_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                A.obj, i)
            err == 0 || throw(ElError(err))
            return i[]
        end

        function numLocalEntries(A::DistSparseMatrix{$elty})
            n = Ref{ElInt}(0)
            err = ccall(($(string("ElDistSparseMatrixNumLocalEntries_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                A.obj, n)
            err == 0 || throw(ElError(err))
            return n[]
        end

        function processQueues(A::DistSparseMatrix{$elty})
            err = ccall(($(string("ElDistSparseMatrixProcessQueues_", ext)), libEl), Cuint,
                (Ptr{Void},),
                A.obj)
            err == 0 || throw(ElError(err))
            return nothing
        end

        function queueLocalUpdate(A::DistSparseMatrix{$elty}, localRow::Integer, col::Integer, value::$elty)
            err = ccall(($(string("ElDistSparseMatrixQueueLocalUpdate_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt, $elty),
                A.obj, localRow-1, col-1, value)
            err == 0 || throw(ElError(err))
            return nothing
        end

        function queueUpdate(A::DistSparseMatrix{$elty}, row::Integer, col::Integer, value::$elty, passive::Bool = true)
            err = ccall(($(string("ElDistSparseMatrixQueueUpdate_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt, $elty, Bool),
                A.obj, row-1, col-1, value, passive)
            err == 0 || error("something is wrong here!")
            return nothing
        end

        function reserve(A::DistSparseMatrix{$elty}, numLocalEntries::Integer, numRemoteEntries::Integer = 0)
            err = ccall(($(string("ElDistSparseMatrixReserve_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt),
                A.obj, numLocalEntries, numRemoteEntries)
            err == 0 || throw(ElError(err))
            return nothing
        end

        function resize!(A::DistSparseMatrix{$elty}, height::Integer, width::Integer = 1) # to mimic vector behavior
            err = ccall(($(string("ElDistSparseMatrixResize_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt),
                A.obj, height, width)
            err == 0 || throw(ElError(err))
            return A
        end

        function width(A::DistSparseMatrix{$elty})
            i = Ref{ElInt}(0)
            err = ccall(($(string("ElDistSparseMatrixWidth_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                A.obj, i)
            err == 0 || throw(ElError(err))
            return i[]
        end
    end
end

# The other constructors don't have a version with dimensions. Should they, or should this one go?
function DistSparseMatrix{T}(::Type{T}, m::Integer, n::Integer, comm::ElComm = CommWorld)
    A = DistSparseMatrix(T, comm)
    resize!(A, m, n)
    return A
end

# Julia convenience

countnz(A::DistSparseMatrix) = MPI.allreduce(Cint(numLocalEntries(A)), +, comm(A))

function showarray(io::IO, A::DistSparseMatrix; kwargs...)
    print(io, "$(size(A, 1))x$(size(A, 2)) $(typeof(A))")
end


text:
#
# xsduri may URL as well as local file
#
function parse_xsd(xsduri)

	xsddoc = parse_file(xsduri)

	xsdroot = root(xsddoc)
	if name(xsdroot) != "schema"
		error("URI does not reference an XSD: the root is a not a schema element")
	end

	astroot = ASTNode(:xsd)

	if has_attribute(xsdroot, "targetNamespace")
		astroot.args[:targetNamespace] = attribute(xsdroot, "targetNamespace")
	end

	# TODO: LightXML does not currently have support for namespaces, so cannot retrieve additional namespace information

	parse_xsd_elements(xsdroot, astroot, xsduri)

	astroot

end

function parse_xsd_attribute(element::XMLElement, node::ASTNode, attributename)
		if has_attribute(element, attributename)
 			node.args[Symbol(attributename)] = attribute(element, attributename)
		end
end

function parse_xsd_elements(parentelement::XMLElement, parentnode::ASTNode, xsduri)

	for element in child_elements(parentelement)

		elementname = name(element)

		if elementname in ["element", "attribute", "group", "attributeGroup",
														"simpleType", "complexType",
														"simpleContent", "complexContent",
														"union", "restriction", "extension", "list",
														"all", "choice", "sequence",
														"enumeration", "pattern", "fractionDigits", "length", "maxExclusive", "maxInclusive",
															"maxLength", "minExclusive", "minInclusive", "minLength", "totalDigits","whiteSpace",
														"any", "anyAttribute"]

			node = ASTNode(Symbol(elementname))

			if elementname in ["element", "attribute", "group", "attributeGroup", "simpleType", "complexType"]
				parse_xsd_attribute(element, node, "name")
			end

			if elementname in ["element"]
				parse_xsd_attribute(element, node, "abstract")
				parse_xsd_attribute(element, node, "substitutionGroup")
			end

			if elementname in ["element", "attribute", "group", "attributeGroup"]
				parse_xsd_attribute(element, node, "ref")
			end

			if elementname in ["element", "attribute"]
				parse_xsd_attribute(element, node, "type")
			end

			if elementname in ["element", "group", "all", "sequence", "choice"]
				parse_xsd_attribute(element, node, "minOccurs")
				parse_xsd_attribute(element, node, "maxOccurs")
			end

			if elementname in ["element", "attribute"]
				parse_xsd_attribute(element, node, "fixed")
			end

			if elementname in ["attribute"]
				parse_xsd_attribute(element, node, "use")
			end

			if elementname in ["complexType", "complexContent"]
				parse_xsd_attribute(element, node, "mixed")
			end

			if elementname in ["union"]
				parse_xsd_attribute(element, node, "memberTypes")
			end

			if elementname in ["restriction"]
				parse_xsd_attribute(element, node, "base")
			end

			if elementname in ["extension"]
				parse_xsd_attribute(element, node, "base")
			end

			if elementname in ["enumeration", "pattern", "fractionDigits", "length", "maxExclusive", "maxInclusive",
				"maxLength", "minExclusive", "minInclusive", "minLength", "totalDigits","whiteSpace"]
				parse_xsd_attribute(element, node, "value")
			end

			if elementname in ["list"]
				parse_xsd_attribute(element, node, "itemType")
			end

			if elementname in ["any", "anyAttribute"]
				parse_xsd_attribute(element, node, "namespace")
			end

			srcelement = new_element(elementname)
			set_attributes(srcelement, attributes_dict(element))
			node.source = split(string(srcelement),'\n')

			push!(parentnode.children, node)

			parse_xsd_elements(element, node, xsduri)

		elseif elementname == "annotation"

			parentnode.source = [parentnode.source; split(string(element),'\n')]

		elseif elementname == "include"

			schemalocation = attribute(element, "schemaLocation")

			# Some simple logic regarding absolute and relative URIs
			# But needs to be improved (perhaps URIParser package etc.)

			newslashpos = rsearch(schemalocation,'/')
			if newslashpos > 0
				# assume absolute
				newxsduri = schemalocation
			else
				# assume relative
				oldslashpos = rsearch(xsduri,'/')
				if oldslashpos == 0
					includexsduri = schemalocation
				else
					includexsduri = xsduri[1:oldslashpos] * schemalocation
				end
			end

			info("including schema from $includexsduri")

			# parse new URI
			includedastroot = parse_xsd(includexsduri)

			# include all elements from under root in the AST
			parentnode.children = [parentnode.children; includedastroot.children]

		else

			error("Do not know how to process XSD element $(elementname)")

		end

	end

end


text:
module Test3D

using Test
using Trixi2Img

include("test_trixi2img.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "3D" begin
  run_trixi(joinpath("3d", "elixir_advection_extended.jl"), maxiters=1)

  @testset "uniform mesh as PNG" begin
    @testset "default slice" begin
      test_trixi2img("solution_000000.h5", outdir)
        # To be able to compare hashes, we would need to make sure that exactly the same versions
        # of all required libraries are installed everywhere, which does not seem to be a good option
        # right now.
    end

    @testset "x-axis slice" begin
      test_trixi2img("solution_000000.h5", outdir, slice_axis=:x, slice_axis_intercept=0.6)
    end

    @testset "y-axis slice" begin
      test_trixi2img("solution_000000.h5", outdir, slice_axis=:y, slice_axis_intercept=-0.35)
    end

    @testset "z-axis slice" begin
      test_trixi2img("solution_000000.h5", outdir, slice_axis=:z, slice_axis_intercept=-0.89)
    end

    @testset "negative border slice" begin
      test_trixi2img("solution_000000.h5", outdir, slice_axis=:y, slice_axis_intercept=-1)
    end

    @testset "positive border slice" begin
      test_trixi2img("solution_000000.h5", outdir, slice_axis=:z, slice_axis_intercept=1)
    end

    @testset "outside negative border slice" begin
      slice_axis_intercept=-1.01
      slice_axis=:x
      @test_throws ErrorException(string(
          "Slice plane $slice_axis = $slice_axis_intercept outside of domain. ",
          "$slice_axis must be between -1.0 and 1.0")) trixi2img(
          joinpath(outdir, "solution_000000.h5"); output_directory=outdir,
          slice_axis=slice_axis, slice_axis_intercept=slice_axis_intercept)
    end

    @testset "outside positive border slice" begin
      slice_axis_intercept=1.005
      slice_axis=:y
      @test_throws ErrorException(string(
          "Slice plane $slice_axis = $slice_axis_intercept outside of domain. ",
          "$slice_axis must be between -1.0 and 1.0")) trixi2img(
          joinpath(outdir, "solution_000000.h5"); output_directory=outdir,
          slice_axis=slice_axis, slice_axis_intercept=slice_axis_intercept)
    end
  end

  # PDF tests do not work on Windows
  if !Sys.iswindows()
    @testset "uniform mesh as PDF with grid lines" begin
      test_trixi2img("solution_000000.h5", outdir,
          format=:pdf, grid_lines=true)
          # To be able to compare hashes, we would need to make sure that exactly the same versions
          # of all required libraries are installed everywhere, which does not seem to be a good option
          # right now.
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end


text:
# -----------------------------------------------------------------------
# Initiate population randomly accordingly to a distribution
# -----------------------------------------------------------------------
# Required packages
using Interpolations, StatsBase
include("stepper_move.jl")

"""
    randweights(coord::AbstractArray, p::AbstractVector, α = 1.0)

## Description
Find the weights for a random sampling considering one dimension only (used for the upstream cylinder position)
## Arguments
- coord: coordinate array
- p: pressure array measured at coord
- α: magnification factor (p^α)
## Output
A weight function (linearly interpolated) is returned
## Example
```jldoctest
julia> x = 0:15:90;
julia> p = [-95,-90,-80,-72,-80,-90,-95];

julia> itpx = randweights(x, p, α=2.0);
julia> itpx(50)
0.4156138645934563
```
"""
function randweights(coord::AbstractArray, p::AbstractVector; α = 1.0)

    val = p .- maximum(p)
    val = val .- minimum(val)
    prob = val.^α ./ sum(val.^α)
    itp = LinearInterpolation(coord, prob)

    return itp
end

"""
    randpos(n::Int, lims::AbstractArray)

## Description
Function to generate random position of OneCylType
## Arguments
- n: number of positions
- lims: array with limits
- weights: array with probabilities for the upstream cylinder position (x and y positions)
## Example
```jldoctest
julia> p = randompos(5, [0:90,-25:25])
5-element Array{OneCylType,1}:
 OneCylType(50.0, -23.0)
 OneCylType(76.0, -22.0)
 OneCylType(85.0, 18.0)
 OneCylType(48.0, 15.0)
 OneCylType(6.0, 25.0)
```
"""
function randpos(n::Int, lims::AbstractArray)

    if length(lims) == 2
        pos = OneCylType.(sample(lims[1],n), sample(lims[2],n))
    elseif length(lims) == 4
        pos = TwoCylType.(sample(lims[1],n), sample(lims[2],n), sample(lims[3],n), sample(lims[4],n))
    elseif length(lims) == 6
        pos = ThreeCylType.(sample(lims[1],n), sample(lims[2],n), sample(lims[3],n), sample(lims[4],n),
                            sample(lims[5],n), sample(lims[6],n))
    else
        throw(BoundsError)
    end
end

function randpos(n::Int, lims::AbstractArray, weigths::AbstractArray)

    if length(lims) == 2
        pos = OneCylType.(sample(lims[1],Weights(weigths[1]),n), sample(lims[2],Weights(weigths[2]),n))
    elseif length(lims) == 4
        pos = TwoCylType.(sample(lims[1],Weights(weigths[1]),n), sample(lims[2],Weights(weigths[2]),n),
                        sample(lims[3],n), sample(lims[4],n))
    elseif length(lims) == 6
        pos = ThreeCylType.(sample(lims[1],Weights(weigths[1]),n), sample(lims[2],Weights(weigths[2]),n),
                        sample(lims[3],n), sample(lims[4],n), sample(lims[5],n), sample(lims[6],n))
    else
        throw(BoundsError)
    end
end

"""
    randindividual(nc=3; size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4], d=5.0)

## Description
Random initialization
## Arguments
- nc: number of cylinders
- size: size of the search 2^size (mm)
- coord_zero: initial coordinates (offset)
- d: cylinder diameter + offset (mm)
## Example
```jldoctest
julia> randindividual(1)
```
"""
function randindividual(nc=3; size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4], d=5.0)
    xs = []
    ys = []

    if length(size) != length(coord_zero)
        error("size and coord_zero must have the same length!")
    end

    for i in 1:2:length(size)
        push!(xs, rand(0:2^(size[i])-1) + coord_zero[i])
        k = i+1
        push!(ys, rand(0:2^(size[k])-1) + coord_zero[k])
    end

    if nc == 1
        p = OneCylType(xs[1], ys[1])
    elseif nc == 2
        p = TwoCylType(xs[1], ys[1], xs[2]+d, ys[2])
    elseif nc == 3
        p = ThreeCylType(xs[1], ys[1], xs[2]+d, ys[2], xs[3]+d, ys[3])
    else
        error("Number of cylinders not defined")
    end
    return p
end

"""
    population_init(n, nc=3; size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4], d=5.0)

## Description
Create initial population with random individuals.
## Arguments
- n: number os individuals
- nc: number of cylinders
- size: size of the search 2^size (mm)
- coord_zero: initial coordinates (offset)
- d: cylinder diameter + offset (mm)
## Example
```jldoctest
# Population with 2 control cylinders and 20 individuals
julia> p = population_init(20,2);
```
"""
function population_init(n, nc=3; size=[5,4,3,3,3,3], coord_zero=[0,-5,0,-4,0,-4], d=5.0)
    P = [randindividual(nc, size=size, coord_zero=coord_zero, d=d) for _ in 1:n]
end


text:
#using LinearAlgebraicRepresentation
#Lar = LinearAlgebraicRepresentation

include("./utilities.jl")
include("./planar_arrangement.jl")
include("./dimension_travel.jl")
include("./largrid.jl")

include("./interface.jl")

include("./struct.jl")
include("./mapper.jl")
include("./integr.jl")
include("./simplexn.jl")
include("./refactoring.jl")
include("./fenvs.jl")


text:
# functions to generate and manipulate masks

""" Read mask containing air wavelengths in csv format.
   input: filename for mask with columns: `lambda` in air and `depth`
   output: DataFrame containing the line wavlengths in vacuum and line depths
"""
function read_mask_air(fn::String)
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda","depth"])
    @assert issorted(df[!,:lambda])
    df[!,:lambda] .= λ_air_to_vac.(df[!,:lambda]) #convert air to vacuum wavelength
    return df
end

""" Read mask containing vacuum wavelengths in csv format.
   input: filename for mask with columns: `lambda` in vacuum and `depth`
   output: DataFrame containing the line wavlengths in vacuum and line depths
"""
function read_mask_vacuum(fn::String)
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda","depth"])
    @assert issorted(df[!,:lambda])
    return df
end

""" Read mask in csv format. Filename must include the substring _AIR" or "_VACUUM" denoting wavelengths medium
   input: filename for mask with columns: `lambda` in vacuum and `depth`
   output: DataFrame containing the line wavlengths in vacuum and line depths
"""
function read_mask(fn::String)
    @assert occursin("_VACUUM",fn) | occursin("_AIR",fn)
    local df = CSV.read(fn,DataFrame,threaded=false,header=["lambda","depth"])
    @assert issorted(df[!,:lambda])
    if occursin("_AIR",fn)
      df[!,:lambda] .= λ_air_to_vac.(df[!,:lambda]) #convert air to vacuum wavelength
    end
    return df
end

function get_airVacString(inst::Module)
   local airVacString
   if inst in (NEID, EXPRES)
      airVacString="_VACUUM"
   elseif inst in (HARPS, HARPSN)
      airVacString="_AIR"
   end
   return airVacString
end

function binMask(mask::DataFrame, nbin::Int; Inst::Module = EXPRES, orders_to_use::UnitRange{Int64} = 43:72, binParam::Symbol = :weight, depthPercentile::Bool = true)
    bins = Dict()
    if depthPercentile #this should be updated to include blaze/normalization/some sort of SNR
        pSorted = sortperm(mask[:,binParam])
        dCumSum = cumsum(mask[pSorted,:weight])
        for i in 1:nbin
            dLower = dCumSum[end]/nbin * (i-1)
            dUpper = dCumSum[end]/nbin * i
            push!(bins, i=>pSorted[(dCumSum .>= dLower) .& (dCumSum .<= dUpper)])
        end
    else
        for i in 1:nbin
            pLower = percentile(mask[:,binParam], 100/nbin * (i-1))
            pUpper = percentile(mask[:,binParam], 100/nbin * i)
            push!(bins, i=>findall((mask[:,binParam] .>= pLower) .& (mask[:,binParam] .< pUpper)))
            if i==nbin
                append!(bins[i],argmax(mask[:,binParam]))
                sort!(bins[8])
            end
        end
    end
    return [sort(mask[bins[i],:], :lambda) for i in 1:nbin]
end

function get_param_range(mask::DataFrame, param::Symbol, param_min::Number, param_max::Number)
   @assert String(param) in names(mask)
   @assert param_min < param_max
   indices = sortperm(mask[!,param])
   index_min = findfirst(mask[indices,param] .> param_min)
   index_max = findlast(mask[indices,param] .< param_max)
   view(mask,sort(indices[index_min:index_max]), :)
end

function get_depth_range(mask::DataFrame, depth_min::Number, depth_max::Number)
   get_param_range(mask, :depth, depth_min, depth_max)
end

function get_lambda_range(mask::DataFrame, lambda_min::Number, lambda_max::Number)
   get_param_range(mask, :lambda, lambda_min, lambda_max)
end

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

using FEMBase
using FEMTruss
using FEMBase: get_formulation_type

using Base.Test
# thsi is the same for all tests testing K and orientation
function make_test_problem(coords::Vector{Tuple{Int64, Vector{Float64}}}, dim::Int64, t::Float64=0.0)
    elem_id = map(x->x[1], coords)
    coord_dict=Dict(coords)
    elem1 = Element(Seg2, elem_id) # connects to nodes 1, 2
    elem1.id = 1
    update!(elem1, "geometry", coord_dict)
    update!(elem1, "youngs modulus", 288.0)
    update!(elem1, "cross section area", 0.1)
    p1 = Problem(Truss, "my truss problem", dim) # 1 dofs/node for now;
    empty!(p1.assembly)
    assemble!(p1.assembly, p1, elem1, t)
    return p1
end

@testset "FEMTrussElements.jl" begin

zeros_2d = zeros(1,2)
zeros_3d = zeros(1,3)
K_oracle_base =  [2.88 -2.88;-2.88 2.88]
#truss test
X = [(1 , [0.0]),(2, [10.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)

@test get_unknown_field_name(p1)=="displacement"

#This one fails
#Test threw an exception of type UndefVarError
#  Expression: get_formulation_type(p1)
#  UndefVarError: get_formulation_type not defined
@test get_formulation_type(p1) == :total

K_truss = full(p1.assembly.K)
K_oracle =  K_oracle_base
@test isapprox(K_truss, K_oracle)

# now for the 2d x case
X = [(1 , [0.0,0.0]),(2, [10.0,0.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [1 0]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# lets do y as well
X = [(1 , [0.0,0.0]),(2, [0.0,10.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [0 1]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# now for the 2d x case reverted
X = [(2 , [0.0,0.0]),(1, [10.0,0.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [-1 0]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# lets do y as well reverted
X = [(2 , [0.0,0.0]),(1, [0.0,10.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [0 -1]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# now for the 3d case
X = [(1 , [0.0,0.0,0.0]),(2, [10.0,0.0,0.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_3d = [1 0 0]
trans = [trans_3d zeros_3d;zeros_3d trans_3d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

end

@testset "test nodal elements" begin
    element = Element(Poi1, [1])
    update!(element, "fixed displacement 1", 1.0)
    update!(element, "nodal force 2", 2.0)
    problem = Problem(Truss, "test problem", 2)
    add_elements!(problem, [element])
    assemble!(problem, 0.0)
    C1 = full(problem.assembly.C1, 2, 2)
    C2 = full(problem.assembly.C2, 2, 2)
    f = full(problem.assembly.f, 2, 1)
    g = full(problem.assembly.g, 2, 1)
    @test isapprox(C1, [1.0 0.0; 0.0 0.0])
    @test isapprox(C1, C2)
    @test isapprox(f, [0.0, 2.0])
    @test isapprox(g, [1.0, 0.0])
end


text:
using TensorBoardLogger #import the TensorBoardLogger package
using Logging #import Logging package

logger = TBLogger("scalarlogs", tb_append) #create tensorboard logger

################log scalars example: y = x²################
#using logger interface
with_logger(logger) do
    for x in 1:20
        @info "scalar/loggerinterface" y = x*x
    end
end
#using explicit function interface
for x in 1:20
    log_value(logger, "scalar/explicitinterface", x*x, step = x)
end

################log scalar example: y = x-xi################
with_logger(logger) do
    for x in 1:10
        z = x-x*im
        @info "scalar/complex" y = z
    end
end