text:
# compute the full right side of the BSEs for all channels
function compute_Γ!(
    Λ      :: Float64,
    r      :: Reduced_lattice,
    m      :: Mesh,
    a1     :: Action,
    a2     :: Action,
    tbuffs :: Vector{NTuple{3, Matrix{Float64}}},
    temps  :: Vector{Array{Float64, 3}},
    corrs  :: Array{Float64, 3},
    eval   :: Int64,
    Γ_tol  :: NTuple{2, Float64}
    )      :: Nothing

    @sync begin
        for w1 in 1 : m.num_Ω
            # compute q1
            Threads.@spawn begin
                compute_channel_s_BSE!(Λ, 1, w1, -1, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                compute_channel_t_BSE!(Λ, 1, w1, -1, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                compute_channel_u_BSE!(Λ, 1, w1, -1, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
            end

            for w3 in 1 : m.num_ν
                # compute q2_1
                Threads.@spawn begin
                    compute_channel_s_BSE!(Λ, 2, w1, w3, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_BSE!(Λ, 2, w1, w3, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_BSE!(Λ, 2, w1, w3, -1, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                # compute q2_2
                Threads.@spawn begin
                    compute_channel_s_BSE!(Λ, 3, w1, -1, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_BSE!(Λ, 3, w1, -1, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_BSE!(Λ, 3, w1, -1, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                for w2 in w3 : m.num_ν
                    # compute q3
                    Threads.@spawn begin
                        compute_channel_s_BSE!(Λ, 4, w1, w2, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_BSE!(Λ, 4, w1, w2, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_BSE!(Λ, 4, w1, w2, w3, r, m, a1, a2, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end
                end
            end
        end
    end

    symmetrize!(r, a2)

    return nothing
end

# compute the full right side of the Katanin truncated flow equations for all channels
function compute_dΓ_1l!(
    Λ      :: Float64,
    r      :: Reduced_lattice,
    m      :: Mesh,
    a      :: Action,
    da     :: Action,
    tbuffs :: Vector{NTuple{3, Matrix{Float64}}},
    temps  :: Vector{Array{Float64, 3}},
    corrs  :: Array{Float64, 3},
    eval   :: Int64,
    Γ_tol  :: NTuple{2, Float64}
    )      :: Nothing

    # compute boundary corrections
    compute_corrs_kat!(Λ, m, a, da, corrs, Γ_tol)

    @sync begin
        for w1 in 1 : m.num_Ω
            # compute q1
            Threads.@spawn begin
                compute_channel_s_kat!(Λ, 1, w1, -1, -1, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                compute_channel_t_kat!(Λ, 1, w1, -1, -1, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                compute_channel_u_kat!(Λ, 1, w1, -1, -1, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
            end

            for w3 in 1 : m.num_ν
                # compute q2_1
                Threads.@spawn begin
                    compute_channel_s_kat!(Λ, 2, w1, w3, -1, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_kat!(Λ, 2, w1, w3, -1, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_kat!(Λ, 2, w1, w3, -1, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                # compute q2_2
                Threads.@spawn begin
                    compute_channel_s_kat!(Λ, 3, w1, -1, w3, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_kat!(Λ, 3, w1, -1, w3, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_kat!(Λ, 3, w1, -1, w3, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                for w2 in w3 : m.num_ν
                    # compute q3
                    Threads.@spawn begin
                        compute_channel_s_kat!(Λ, 4, w1, w2, w3, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_kat!(Λ, 4, w1, w2, w3, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_kat!(Λ, 4, w1, w2, w3, r, m, a, da, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end
                end
            end
        end
    end

    symmetrize!(r, da)

    return nothing
end

# compute the full right side of the two loop truncated flow equations for all channels
function compute_dΓ_2l!(
    Λ      :: Float64,
    r      :: Reduced_lattice,
    m      :: Mesh,
    a      :: Action,
    da     :: Action,
    da_l   :: Action,
    tbuffs :: Vector{NTuple{3, Matrix{Float64}}},
    temps  :: Vector{Array{Float64, 3}},
    corrs  :: Array{Float64, 3},
    eval   :: Int64,
    Γ_tol  :: NTuple{2, Float64}
    )      :: Nothing

    # compute one loop
    compute_dΓ_1l!(Λ, r, m, a, da, tbuffs, temps, corrs, eval, Γ_tol)

    # compute boundary corrections
    compute_corrs!(Λ, m, a, corrs, Γ_tol)

    @sync begin
        for w1 in 1 : m.num_Ω
            # compute q1
            Threads.@spawn begin
                compute_channel_s_left!(Λ, 1, w1, -1, -1, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                compute_channel_t_left!(Λ, 1, w1, -1, -1, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                compute_channel_u_left!(Λ, 1, w1, -1, -1, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
            end

            for w3 in 1 : m.num_ν
                # compute q2_1
                Threads.@spawn begin
                    compute_channel_s_left!(Λ, 2, w1, w3, -1, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_left!(Λ, 2, w1, w3, -1, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_left!(Λ, 2, w1, w3, -1, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                # compute q2_2
                Threads.@spawn begin
                    compute_channel_s_left!(Λ, 3, w1, -1, w3, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_left!(Λ, 3, w1, -1, w3, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_left!(Λ, 3, w1, -1, w3, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                for w2 in 1 : m.num_ν
                    # compute q3
                    Threads.@spawn begin
                        compute_channel_s_left!(Λ, 4, w1, w2, w3, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_left!(Λ, 4, w1, w2, w3, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_left!(Λ, 4, w1, w2, w3, r, m, a, da, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end
                end
            end
        end
    end

    symmetrize_add_to!(r, da_l, da)

    return nothing
end

# compute the full right side of the multiloop truncated flow equations for all channels
function compute_dΓ_ml!(
    Λ       :: Float64,
    r       :: Reduced_lattice,
    m       :: Mesh,
    loops   :: Int64,
    a       :: Action,
    da      :: Action,
    da_l    :: Action,
    da_c    :: Action,
    da_temp :: Action,
    da_Σ    :: Action,
    tbuffs  :: Vector{NTuple{3, Matrix{Float64}}},
    temps   :: Vector{Array{Float64, 3}},
    corrs   :: Array{Float64, 3},
    eval    :: Int64,
    Γ_tol   :: NTuple{2, Float64}
    )       :: Nothing

    # compute two loop
    compute_dΓ_2l!(Λ, r, m, a, da, da_l, tbuffs, temps, corrs, eval, Γ_tol)

    # update temporary buffer and reset terms for self energy corrections
    reset_Γ!(da_temp)
    reset_Γ!(da_Σ)
    symmetrize_add_to!(r, da_l, da_temp)

    for loop in 3 : loops
        @sync begin
            for w1 in 1 : m.num_Ω
                # compute q1
                Threads.@spawn begin
                    compute_channel_s_central!(Λ, 1, w1, -1, -1, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_central!(Λ, 1, w1, -1, -1, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_central!(Λ, 1, w1, -1, -1, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                for w3 in 1 : m.num_ν
                     # compute q2_1
                    Threads.@spawn begin
                        compute_channel_s_central!(Λ, 2, w1, w3, -1, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_central!(Λ, 2, w1, w3, -1, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_central!(Λ, 2, w1, w3, -1, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end

                     # compute q2_2
                    Threads.@spawn begin
                        compute_channel_s_central!(Λ, 3, w1, -1, w3, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_central!(Λ, 3, w1, -1, w3, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_central!(Λ, 3, w1, -1, w3, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end

                    for w2 in w3 : m.num_ν
                         # compute q3
                        Threads.@spawn begin
                            compute_channel_s_central!(Λ, 4, w1, w2, w3, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                            compute_channel_t_central!(Λ, 4, w1, w2, w3, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                            compute_channel_u_central!(Λ, 4, w1, w2, w3, r, m, a, da_l, da_c, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        end
                    end
                end
            end
        end

        @sync begin
            for w1 in 1 : m.num_Ω
                # compute q1
                Threads.@spawn begin
                    compute_channel_s_left!(Λ, 1, w1, -1, -1, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_t_left!(Λ, 1, w1, -1, -1, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    compute_channel_u_left!(Λ, 1, w1, -1, -1, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                end

                for w3 in 1 : m.num_ν
                    # compute q2_1
                    Threads.@spawn begin
                        compute_channel_s_left!(Λ, 2, w1, w3, -1, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_left!(Λ, 2, w1, w3, -1, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_left!(Λ, 2, w1, w3, -1, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end

                    # compute q2_2
                    Threads.@spawn begin
                        compute_channel_s_left!(Λ, 3, w1, -1, w3, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_t_left!(Λ, 3, w1, -1, w3, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        compute_channel_u_left!(Λ, 3, w1, -1, w3, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                    end

                    # compute q3
                    for w2 in 1 : m.num_ν
                        Threads.@spawn begin
                            compute_channel_s_left!(Λ, 4, w1, w2, w3, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                            compute_channel_t_left!(Λ, 4, w1, w2, w3, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                            compute_channel_u_left!(Λ, 4, w1, w2, w3, r, m, a, da_temp, da_l, tbuffs[Threads.threadid()], temps[Threads.threadid()], corrs, eval, Γ_tol)
                        end
                    end
                end
            end
        end

        # update temporary buffer
        symmetrize!(r, da_c)
        replace_with_Γ!(da_temp, da_c)
        symmetrize_add_to!(r, da_l, da_temp)

        # update self energy corrections and flow
        add_to_Γ!(da_c, da_Σ)
        add_to_Γ!(da_temp, da)
    end

    return nothing
end


text:
# Autogenerated wrapper script for xsum_jll for x86_64-linux-musl
export libxsum

JLLWrappers.@generate_wrapper_header("xsum")
JLLWrappers.@declare_library_product(libxsum, "libxsum.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libxsum,
        "lib/libxsum.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
include("admm_common.jl")

type RegressorSelectionState <: NonconvexADMMState
    rho::Float64
    K::Int
    A::Matrix{Float64}
    b::Vector{Float64}
    m::Int
    n::Int
    At::Matrix{Float64}
    Atb::Vector{Float64}
    x::Vector{Float64}
    z::Vector{Float64}
    u::Vector{Float64}
    L::Matrix{Float64}
    U::Matrix{Float64}
    function RegressorSelectionState(A::Matrix{Float64}, b::Vector{Float64},
        K::Int, rho::Float64)
        m, n = size(A)
        At = A'
        Atb = At * b            # save a matrix-vector multiply
        L, U = factor(A, rho)   # cache the factorization
        new(rho, iround(K), A, b, m, n, At, Atb,
            zeros(n), zeros(n), zeros(n), L, U)
    end
end

function regressor_selection_admm(A, b, K, rho::Float64 = 1.0; options...)
    as = RegressorSelectionState(A, vec(b), iround(K), rho)
    admm_loop(as; options...)
end

function update_xzu!(s::RegressorSelectionState, iter::Int)
    # x-update
    q = s.Atb + s.rho*(s.z - s.u)    # temporary value
    if ( s.m >= s.n )    # if skinny
        x = s.U \ (s.L \ q)
    else            # if fat
        x = q/s.rho - (s.At*(s.U \ ( s.L \ (s.A*q) )))/s.rho^2
    end
    s.x = vec(x)

    # z-update with relaxation
    s.z = keep_largest!(s.x + s.u, s.K)

    # u-update
    s.u = s.u + (s.x - s.z)
end

solution(s::RegressorSelectionState) = s.z

function objective(s::RegressorSelectionState)
    sum( (s.A*s.x - s.b).^2 )
end

function keep_largest!(z::Vector{Float64}, K::Int)
    pos = sortperm(abs(z), rev=true)
    z[pos[(K+1):end]] = 0
    z
end

import Base.factor
function factor(A, rho)
    m, n = size(A)
    if ( m >= n )  # if skinny
        L = chol( A'*A + rho*speye(n), :L )
    else            # if fat
        L = chol( speye(m) + 1/rho*(A*A'), :L )
    end

    return L, L'
end


text:
@testset "Cycler" begin
    cycl = Cycler([:a, :b, :c])
    @test length(cycl) == 3
    @test !isempty(cycl)
    @test cycl[1] == :a
    @test cycl[5] == :b
    @test cycl[3:7] == [:c, :a, :b, :c, :a]
    @test cycl[0] == :c
    @test cycl[-1] == :b
    @test cycl[-2:1] == [:a, :b, :c, :a]
end

ps1 = M.Style(layout=M.Layout(font_size=10))
ps2 = M.Style(layout=M.Layout(font_family="Helvetica"))

@testset "Style(::Style, ::Style)" begin
    ps3 = M.Style(ps1, ps2)
    @test ps3.layout[:font] == Dict(:family=>"Helvetica", :size=>10)
    @test isempty(ps3.global_trace)
    @test isempty(ps3.trace)

    gg = M.ggplot_style()
    new_font = M.Style(layout=M.Layout(font_family="Source Code Pro"))
    ps4 = M.Style(gg, new_font)
    ps5 = M.Style(new_font, gg)

    # make sure things were set properly
    for (k, v) in gg.layout
        if k == :font
            want = gg.layout[:font]
            want[:family] = "Source Code Pro"
            @test ps4.layout[k] == ps5.layout[k] == want
        else
            @test ps4.layout[k] == ps5.layout[k] == gg.layout[k]
        end
    end

end

@testset "setting plot attributes" begin

    goofy = Style(
        global_trace=attr(
            marker_color="red", line_dash=Cycler(["dash", "dot"]),
        ),
        trace=Dict(:scatter => attr(mode="markers")),
    )

    p1 = Plot(scatter(y=1:3, mode="lines", marker_symbol="square"), style=goofy)
    p2 = Plot(scatter(y=1:3, marker_color="green", line_dash="dot"), style=goofy)

    # now call JSON.lower to enforce style
    PlotlyBase.JSON.lower(p1)
    PlotlyBase.JSON.lower(p2)

    @test p1.data[1]["marker_color"] == "red"
    @test p2.data[1]["marker_color"] == "green"

    @test p1.data[1]["line_dash"] == "dash"
    @test p2.data[1]["line_dash"] == "dot"

    p3 = Plot(rand(5, 6), style=goofy)
    restyle!(p3, 4, line_dash="something_else")
    PlotlyBase.JSON.lower(p3)
    @test p3.data[1]["line_dash"] == "dash"
    @test p3.data[2]["line_dash"] == "dot"
    @test p3.data[3]["line_dash"] == "dash"
    @test p3.data[4]["line_dash"] == "something_else"
    @test p3.data[5]["line_dash"] == "dot"
    @test p3.data[6]["line_dash"] == "dash"
end


text:
#=
Implementation of 32 bit xxHash by Yann Collet

https://github.com/Cyan4973/xxHash
=#
const prime1 = UInt32(2654435761)
const prime2 = UInt32(2246822519)
const prime3 = UInt32(3266489917)
const prime4 = UInt32( 668265263)
const prime5 = UInt32( 374761393)

@inline rotl(x::UInt32, r::Int)::UInt32 = (x << r) | (x >> (32 - r))

function _xxhash32(data, seed::UInt32)::UInt32
    len = length(data)
    acc::UInt32 = UInt32(0)
	i = n = 1
    bytes = Ptr{UInt8}(pointer(data))
    words = Ptr{UInt32}(pointer(data))

	if len >= 16
        limit = len - 15
        acc1::UInt32 = seed + prime1 + prime2
        acc2::UInt32 = seed + prime2
        acc3::UInt32 = seed
        acc4::UInt32 = seed - prime1

        @inbounds for _ = 1:16:limit
            acc1 = rotl(acc1 + unsafe_load(words, n  ) * prime2, 13) * prime1
            acc2 = rotl(acc2 + unsafe_load(words, n+1) * prime2, 13) * prime1
            acc3 = rotl(acc3 + unsafe_load(words, n+2) * prime2, 13) * prime1
            acc4 = rotl(acc4 + unsafe_load(words, n+3) * prime2, 13) * prime1
            i += 16
            n += 4
        end

        acc = rotl(acc1, 1) + rotl(acc2, 7) + rotl(acc3, 12) + rotl(acc4, 18)
    else
        acc = seed + prime5
    end

    acc += len

    @inbounds for _ = i:4:len-3
        acc = acc + unsafe_load(words, n) * prime3
        acc = rotl(acc, 17) * prime4
        i += 4
        n += 1
    end

    @inbounds for _ = i:len
        acc = acc + unsafe_load(bytes, i) * prime5
        acc = rotl(acc, 11) * prime1
        i += 1
    end

    acc = acc ⊻ acc >> 15
    acc = acc * prime2
    acc = acc ⊻ acc >> 13
    acc = acc * prime3
    acc = acc ⊻ acc >> 16
end

const Hashable = Union{AbstractString,Array{UInt8}}
xxhash32(key::Hashable)::UInt32 = _xxhash32(key, UInt32(0))
xxhash32(key::Hashable, seed::UInt32)::UInt32 = _xxhash32(key, seed)
xxhash32(key::Hashable, seed::Int)::UInt32 = _xxhash32(key, UInt32(seed))

"""
    xxhash32(key [, seed])

32 bit hash of `key` using xxHash, optionally seeded with `seed`.
"""
xxhash32


text:
struct SummaryAlign{T<:Real,S<:Integer}
    Pij::Matrix{T}
    Pi::Vector{T}
    Pijpc::Matrix{T}
    Pipc::Vector{T}
    Header::Vector{String}
    Seqs::Vector{String}
    Z::Matrix{S}
    W::Vector{T}
    M::Int
    N::Int
    q::Int
end

function SummaryAlign(
    Pij::Matrix{T},
    Pi::Vector{T},
    Pijpc::Matrix{T},
    Pipc::Vector{T},
    Header::Vector{String},
    Seqs::Vector{String},
    Z::Matrix{S},
    W::Vector{T},
) where {T<:Real} where {S<:Integer}
    M1, N1 = size(Pij)
    N2 = length(Pi)
    N, M = size(Z)
    if M1 != N1 || N1 != N2
        throw(DimensionMismatch("size(Pij) = $M1,$N1 length(Pi) = $N2. They should be all equal"))
    end

    M2, N2 = size(Pijpc)
    N3 = length(Pipc)

    (M1 == M2 && N1 == N2) ||
        throw(DimensionMismatch("Pij and Pijpc have different sizes"))
    N3 == N2 || throw(DimensionMismatch("Pi and Pij have diffeent sizes"))
    Nz, Mz = size(Z)
    Mw = length(W)
    Mz == Mw || throw(DimensionMismatch("W and Z have different sizes"))
    MH = length(Header)
    MH == Mz || throw(DimensionMismatch("Header and Z have different sizes"))
    MS = length(Seqs)
    MS == MH || throw(DimensionMismatch("Seqs and Z have different sizes"))
    qmin, qmax = extrema(Z)
    qmin < 0 && error("qmin = $qmin < 0")
    qmax > 21 && error("qmax = $qmax > 21")
    q = div(M1, N)
    qmin ≤ q ≤ qmax || error("q = $qmin  ∉ [$qmin,$qmax]")
    return SummaryAlign(Pij, Pi, Pijpc, Pipc, Header, Seqs, Z, W, M1, N, q + 1)
end

"""
    function summary_align(Z::Matrix{Int8};
        pseudocount::Real = 0.8,
        mode::Symbol = :nodiag,
        theta = :auto,
        max_gap_fraction::Real = 0.9,
        score::Symbol = :frob,
        min_separation::Integer = 5,
        remove_dups::Bool = false) -> Corr
Compute one and two points statistics. Returns a `Corr` structure containing
the `N(q-1) × N(q-1)` two point correlation matrix `Pij`, `Pi`, and `Z` the
`q × M` alignment data translated in `Int8`.
"""
function summary_align(
    Z::Array{Int8,2},
    Headers::Array{String,1},
    Seqs::Array{String,1};
    pseudocount::Real = 0.0,
    mode::Symbol = :nodiag,
    theta = 0.0,
    max_gap_fraction::Real = 1.0,
    score::Symbol = :frob,
    min_separation::Integer = 1,
    remove_dups::Bool = true,
)

    0 <= pseudocount <= 1 ||
        throw(DomainError("pseudocount = $pseudocount should be in the interval [0,1]"))

    println("pseudocount = $pseudocount")
    GaussDCA.use_threading(true)
    if remove_dups
        Z, _ = GaussDCA.remove_duplicate_seqs(Z)
    end
    N, M = size(Z)
    q = Int(maximum(Z))
    q > 32 && error("parameter q=$q is too big (max 32 is allowed)")

    Pi_true, Pij_true, Meff, W = GaussDCA.compute_new_frequencies(Z, q, theta)
    N == div(length(Pi_true), q - 1) ||
        throw(DimensionMismatch("something wrong here"))
    W .= W / sum(W)
    Pi, Pij = GaussDCA.add_pseudocount(Pi_true, Pij_true, pseudocount, N, q)

    return SummaryAlign(Pij_true, Pi_true, Pij, Pi, Headers, Seqs, Z, W)
end

"""
    function summary_align(nomefile::String;
        pseudocount::Real = 0.8,
        mode::Symbol = :nodiag,
        theta = :auto,
        max_gap_fraction::Real = 0.9,
        score::Symbol = :frob,
        min_separation::Integer = 5,
        remove_dups::Bool = false) -> SummaryAlign

Compute one and two points statistics of a `FASTA` msa contained in `nomefile`.
Returns a `Corr` structure containingthe `N(q-1) × N(q-1)` two point correlation
matrix `Pij`, `Pi`, and `Z` the`q × M` alignment data translated in `Int8`.
"""
summary_align(nomefile::String; args...) =
    summary_align(read_fasta_alignment(nomefile, 1.0)...; args...)

function read_fasta_alignment(filename::AbstractString, max_gap_fraction::Real)

    f = FastaReader(filename)

    max_gap_fraction = Float64(max_gap_fraction)

    # pass 1

    seqs = Int[]
    inds = Int[]
    fseqlen = 0

    for (name, seq) in f
        ngaps = 0
        if f.num_parsed == 1
            ls = length(seq)
            resize!(inds, ls)
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    fseqlen += 1
                    inds[fseqlen] = i
                    c == '-' && (ngaps += 1)
                end
            end
        else
            ls = length(seq)
            ls == length(inds) || error("inputs are not aligned")
            tstfseqlen = 0
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    tstfseqlen += 1
                    inds[tstfseqlen] == i || error("inconsistent inputs")
                    c == '-' && (ngaps += 1)
                end
            end
            tstfseqlen == fseqlen || error("inconsistent inputs")
        end
        ngaps / fseqlen <= max_gap_fraction && push!(seqs, f.num_parsed)
    end

    length(seqs) > 0 ||
        error("Out of $(f.num_parsed) sequences, none passed the filter (max_gap_fraction=$max_gap_fraction)")

    # pass 2
    Z = Array{Int8}(undef, fseqlen, length(seqs))
    Header = Array{String}(undef, length(seqs))
    Seqs = Array{String}(undef,length(seqs))
    seqid = 1
    for (name, seq) in f
        seqs[end] < f.num_parsed && break
        seqs[seqid] == f.num_parsed || continue
        for i = 1:fseqlen
            c = seq[inds[i]]
            Z[i, seqid] = letter2num(c)
        end
        Header[seqid] = name
        Seqs[seqid] = seq
        seqid += 1
    end
    @assert seqid == length(seqs) + 1

    close(f)

    return Z, Header, Seqs
end

let alphabet = [
        1,
        21,
        2,
        3,
        4,
        5,
        6,
        7,
        8,
        21,
        9,
        10,
        11,
        12,
        21,
        13,
        14,
        15,
        16,
        17,
        21,
        18,
        19,
        21,
        20,
    ]
    # A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y
    global letter2num
    function letter2num(c::Union{Char,UInt8})
        i = UInt8(c) - 0x40
        1 <= i <= 25 && return alphabet[i]
        return 21
    end
end

"""
    function read_stockholm(iFile::String) -> Dict{String,String}

Read stokholm formatted `iFile` and return a `Dict{String,String}` of the alignment:
`keys -> header`, `values -> sequence`.

"""
function read_stockholm(iFile::String)

    fasta_dict = Dict{String,String}()
    header_order = String[]
    open(iFile) do iFid
        for line in eachline(iFid)
            line = chomp(line)
            line == "" && continue
            line[1] == '#' && continue
            line == "//" && continue
            header, seq = split(line)
            if haskey(fasta_dict, header)
                fasta_dict[header] = string(fasta_dict[header], seq)
            else
                fasta_dict[header] = seq
                push!(header_order, header)
            end
        end
    end
    fasta_dict
end

function extract_ins(iFile::String)

    f = FastaReader(iFile)
    fasta_dict = Dict{String, String}()
    for (name,line) in f
        newline = replace(line, r"[.]" => s"")
        f_match = 0
        l_match = length(newline)
        i = 1
        idx_s = 0
        idx_e = 0
        count = 0
        while i <= length(newline)
            if newline[i] != '-'
                count = count + 1
            end
            if (newline[i] == '-' || isuppercase(newline[i])) && idx_s == 0
                idx_s = i
            end
            if isuppercase(newline[i]) && f_match == 0
                f_match = count
            end
            i = i + 1
        end
        i = length(newline)
        count = count + 1
        while i >= 1
            if newline[i] != '-'
                count = count - 1
            end
            if (newline[i] == '-' || isuppercase(newline[i])) && idx_e == 0
                idx_e = i
            end
            if isuppercase(newline[i]) && l_match == length(newline)
                l_match = count
            end
            i = i - 1
        end
        if occursin("/", name)
            newname = name #pfam
        else
            newname = name * "/" * string(f_match) * "-" * string(l_match) #mafft
        end
        fasta_dict[newname] = newline[idx_s:idx_e]
    end

    return fasta_dict

end

function extract_align(f)

    fasta_dict = Dict{String, String}()
    L = 0
    for (name,line) in f
        line = replace(line, r"[.]" => s"")
        newline = ""
        for i = 1:length(line)
            if !islowercase(line[i])
                newline = newline * line[i]
            end
        end
        L = length(newline)
        fasta_dict[name] = newline
    end

    return fasta_dict, L

end
"""
    linsi(iFile::String, oFile::String; nthread::Int=1 eFile::Union{String,Nothing}=nothing)

Runs linsi (mafft accurate alignment strategy) on input fasta file `iFile` and writes
the output in `oFile` (aligned fasta file). Optional arguments are, `nthread::Int`
(default = 1) and `eFile::String` the file to collect errors (default = `nothing`,
i.e it writes on `stderr`)
"""
function linsi(
    iFile::String,
    oFile::String;
    nthreads::Int = 1,
    eFile::Union{String,Nothing} = nothing,
)
    if eFile === nothing
        run(pipeline(
            `linsi --thread $nthreads --anysymbol $iFile`;
            stdout = oFile,
            stderr = string(oFile, ".stderr"),
        ))
    else
        run(pipeline(
            `linsi --thread $nthreads --anysymbol $iFile`;
            stdout = oFile,
            stderr = eFile,
        ))
    end
end

function hmmbuild(iFile::String, oFile::String, hmmFile::String;
                  nthreads::Int=1)
    run(pipeline(`hmmbuild -O $oFile --cpu $nthreads $hmmFile $iFile`, devnull))
end

function hmmalign(hmmFile::String,iFile::String, oFile::String)
    run(pipeline(`hmmalign -o $oFile $hmmFile $iFile`, devnull))
end


text:
# Linear Least Square

function _llsq{T<:BlasReal}(trans::Bool, X::StridedMatrix{T}, y::StridedVecOrMat{T}, bias::T, method::Symbol)

	X_ = if bias == zero(T)
		copy(X)
	else
		trans ? augment_cols(X, bias) :
				augment_rows(X, bias)
	end::Matrix{T}
	ch = trans ? 'T' : 'N'
	y_ = copy(y)
	rcond = eps(convert(T, length(X)))

	r = if method == :qrlq
		gels!(ch, X_, y_)[2]
	elseif method == :orth
		trans ? gelsy!(X_', y_, rcond)[1] : gelsy!(X_, y_, rcond)[1]
	elseif method == :svd
		trans ? gelsd!(X_', y_, rcond)[1] : gelsd!(X_, y_, rcond)[1]
	else
		error("Invalid method value: $method.")
	end

	return r
end

function llsq{T<:BlasReal}(X::StridedMatrix{T}, y::StridedVecOrMat{T};
						   bias::Real=zero(T),
						   method::Symbol=:qrlq,
						   trans::Bool=false)

	_llsq(trans, X, y, convert(T, bias)::T, method)
end

# Ridge regression

function _ridgereg{T<:BlasReal}(trans::Bool, X::StridedMatrix{T}, y::StridedVecOrMat{T}, bias::T, r::T)
	X_ = if bias == zero(T)
		X
	else
		trans ? augment_cols(X, bias) :
				augment_rows(X, bias)
	end

	H = trans ? A_mul_Bt(X_, X_) : At_mul_B(X_, X_)
	n = trans ? size(X, 1) : size(X, 2)
	n_ = bias == zero(T) ? n : n+1
	@assert size(H, 1) == size(H, 2) == n_
	for i = 1:n
		H[i,i] += r
	end
	rhs = trans ? X_ * y : X_'y
	cholfact!(H) \ rhs
end

function ridgereg{T<:BlasReal}(X::StridedMatrix{T}, y::StridedVecOrMat{T}, r::Real;
							   bias::Real=zero(T),
							   trans::Bool=false)

	_ridgereg(trans, X, y, convert(T, bias)::T, convert(T, r)::T)
end


text:
using EzXML: Node, readxml
using Vulkan_Headers_jll: vk_xml

const xdoc = readxml(vk_xml)
const xroot = xdoc.root

const extension_types = [
    :Display,
    :VisualID,
    :Window,
    :RROutput,
    :wl_display,
    :wl_surface,
    :HINSTANCE,
    :HWND,
    :HMONITOR,
    :HANDLE,
    :SECURITY_ATTRIBUTES,
    :DWORD,
    :LPCWSTR,
    :xcb_connection_t,
    :xcb_visualid_t,
    :xcb_window_t,
    :IDirectFB,
    :IDirectFBSurface,
    :zx_handle_t,
    :GgpStreamDescriptor,
    :GgpFrameToken,
    :MirConnection,
    :MirSurface,
    :ANativeWindow,
    :AHardwareBuffer,
    :CAMetalLayer,
    :_screen_context,
    :_screen_window,
]

include("spec/utils.jl")
include("spec/aliases.jl")
include("spec/types.jl")
include("spec/specs.jl")
include("spec/structure_types.jl")

for name ∈ names(@__MODULE__, all = true)
    if any(startswith.(string(name), ["spec_", "Spec"]))
        @eval export $name
    end
end

export
    # Alias manipulation
    follow_alias,
    isalias,
    hasalias,

    # Specifications
    ### Utility
    default,
    func_by_name,
    struct_by_name,
    union_by_name,
    handle_by_name,
    bitmask_by_name,
    flag_by_name,
    enum_by_name,
    constant_by_name,
    create_func,
    create_funcs,
    create_func_no_batch,
    destroy_func,
    destroy_funcs,
    follow_constant,
    children,
    parent_spec,
    has_length,
    has_computable_length,
    is_length,
    is_arr,
    is_size,
    length_chain,
    is_fn_ptr,
    is_version,
    wrappable_constructors,

    ### Render passes
    RenderPassRequirement,
    RenderPassInside,
    RenderPassOutside,
    render_pass_compatibiltiy_map,

    ### Queues
    QueueType,
    QueueCompute,
    QueueGraphics,
    QueueTransfer,
    QueueSparseBinding,

    ### Structures & functions
    PARAM_REQUIREMENT,
    OPTIONAL,
    REQUIRED,
    POINTER_OPTIONAL,
    POINTER_REQUIRED,
    STRUCT_TYPE,
    CREATE_INFO,
    ALLOCATE_INFO,
    GENERIC_INFO,
    DATA,
    PROPERTY,
    FUNC_TYPE,
    CREATE,
    DESTROY,
    ALLOCATE,
    FREE,
    COMMAND,
    QUERY,
    OTHER


text:
#Attempts to read every .mtx file in the test directory
using MatrixMarket
using Compat.Test
using Compat

include("dl-matrixmarket.jl")

num_errors = 0
num_pass = 0

@testset "Idempotent IO" begin
    @testset "read and write $filename" for filename in filter(t -> endswith(t, ".mtx"), readdir())
        new_filename = "$(filename)_"
        A = MatrixMarket.mmread(filename)
        Compat.@info("$(typeof(A))  $(size(A))")

        # verify mmread(mmwrite(A)) == A
        MatrixMarket.mmwrite(new_filename, A)
        new_A = MatrixMarket.mmread(new_filename)
        @test new_A == A

        rm(filename)
        rm(new_filename)
    end
end


text:
# Packages required for testing
using Test
using Logging

# Package Under Test
using WorldMapData

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Define package tests
@time @testset "WorldMapData Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "WorldMapData.LoadData" begin
        include(joinpath(testdir, "test_load_data.jl"))
    end
end


text:
# ========== Single Step Autocorrelation Wavelet Transform ==========
## 1D ##
"""
    acdwt_step(v, j, h, g)

Performs one level of the autocorrelation discrete wavelet transform (acdwt) on the vector
`v`, which is the j-th level scaling coefficients (Note the 0th level scaling coefficients
is the raw signal). The vectors `h` and `g` are the detail and scaling filters.

# Arguments
- `v::AbstractVector{T} where T<:Number`: Vector of coefficients from a node at level `d`.
- `j::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::Vector{T}`: Output from the low pass filter.
- `w₂::Vector{T}`: Output from the high pass filter.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
ACWT.acdwt_step(v, 0, h, g)
```

**See also:** [`acdwt`](@ref), [`iacdwt`](@ref)
"""
function acdwt_step(v::AbstractVector{T}, d::Integer, h::Array{S,1}, g::Array{S,1}) where
                   {T<:Number, S<:Number}
    n = length(v)
    w₁ = zeros(T,n)
    w₂ = zeros(T,n)

    acdwt_step!(w₁, w₂, v, d, h, g)
    return w₁, w₂
end

"""
    acdwt_step!(w₁, w₂, v, j, h, g)

Same with `acdwt_step` but without array allocation.

# Arguments
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.
- `v::AbstractVector{T} where T<:Number`: Vector of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::Vector{T}`: Output from the low pass filter.
- `w₂::Vector{T}`: Output from the high pass filter.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
w₁ = similar(v)
w₂ = similar(v)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
ACWT.acdwt_step!(w₁, w₂, v, 0, h, g)
```

**See also:** [`acdwt_step!`](@ref), [`acdwt`](@ref), [`iacdwt`](@ref)
"""
function acdwt_step!(w₁::AbstractVector{T},
                     w₂::AbstractVector{T},
                     v::AbstractVector{T},
                     j::Integer,
                     h::Array{T,1},
                     g::Array{T,1}) where {T<:Number, S<:Number}
    # Sanity check
    @assert length(w₁) == length(w₂) == length(v)
    @assert length(h) == length(g)

    # Setup
    N = length(v)
    L = length(h)

    # One step of autocorrelation transform
    for i in 1:N
        t = i+(1<<j) |> t -> t>N ? mod1(t,N) : t
        i = mod1(i + (L÷2+1) * 2^j,N) # Need to shift by half the filter size because of periodicity assumption
        @inbounds w₁[i] = g[1] * v[t]
        @inbounds w₂[i] = h[1] * v[t]
        for n in 2:L
            t = t+(1<<j) |> t -> t>N ? mod1(t,N) : t
            @inbounds w₁[i] += g[n] * v[t]
            @inbounds w₂[i] += h[n] * v[t]
        end
    end
    return w₁, w₂
end

"""
    iacdwt_step(w₁, w₂)

Perform one level of the inverse autocorrelation discrete wavelet transform (IACDWT) on the
vectors `w₁` and `w₂`, which are the `j+1`-th level scaling coefficients (Note that the 0th
level scaling coefficients is the raw signal).

# Arguments
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.

# Returns
- `v::Vector{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
w₁, w₂ = ACWT.acdwt_step(v, 0, h, g)

# One step of IACDWT
v̂ = ACWT.iacdwt_step(w₁, w₂)
```

**See also:** [`iacdwt_step!`](@ref), [`acdwt_step`](@ref), [`iacdwt`](@ref)
"""
function iacdwt_step(w₁::AbstractVector{T}, w₂::AbstractVector{T}) where T<:Number
    v = similar(w₁)
    iacdwt_step!(v, w₁, w₂)
    return v
end

"""
    iacdwt_step!(v, w₁, w₂)

Same as `iacdwt_step` but without array allocation.

# Arguments
- `v::AbstractVector{T} where T<:Number`: Vector allocation for reconstructed coefficients.
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.

# Returns
- `v::Vector{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.make_acreverseqmfpair(wt)

# One step of ACDWT
w₁, w₂ = ACWT.acdwt_step(v, 0, h, g)

# One step of IACDWT
v̂ = similar(v)
ACWT.iacdwt_step!(v̂, w₁, w₂)
```

**See also:** [`iacdwt_step`](@ref), [`acdwt_step`](@ref), [`iacdwt`](@ref)
"""
function iacdwt_step!(v::AbstractVector{T},
                      w₁::AbstractVector{T},
                      w₂::AbstractVector{T}) where T<:Number
    @assert length(v) == length(w₁) == length(w₂)
    for i in eachindex(v)
        @inbounds v[i] = (w₁[i]+w₂[i]) / √2
    end
end


text:
function rand(internal::TwoInternalStates{T})::T where {T<:IUI32}
    s1, s2 = state(internal)

	nextvalue = s1 * 0x9E3779BB
	nextvalue = rotl(nextvalue, 5)
	nextvalue += nextvalue << 2

	s2 = s2 ⊻ s1
	s1 = rotl(s1, 26)
	s1 = s1 ⊻ (s2 ⊻ (s2 << 9))
	s2 = rotl(s2, 13)

	state(internal, s1, s2)

	return nextvalue
end

function rand(internal::TwoInternalStates{T}) where {T<:IUI64}
    s1, s2 = state(internal)

	nextvalue =  s1
	nextvalue += s1 << 2
	nextvalue =  rotl(nextvalue, 7)
	nextvalue += nextvalue << 3

	s2 ⊻= s1
	s1 = rotl(s1, 24)
	s1 ⊻= s2 ⊻ (s2 << 16)
	s2 = rotl(s2, 37)

	state(internal, s1, s2)

	return nextvalue
end

function rand(internal::FourInternalStates{T}) where {T<:IUI64}
    s1, s2, s3, s4 = state(internal)

	nextvalue =  s2
	nextvalue += s2 << 2
	nextvalue =  rotl(nextvalue, 7)
	nextvalue += nextvalue << 3

	t = s2 << 17

	s3 ⊻= s1
	s4 ⊻= s2
	s2 ⊻= s3
	s1 ⊻= s4

	s3 ⊻= t
	s4 = rotl(s4, 45)

	state(internal, s1, s2, s3, s4)

	return nextvalue
end

function rand(internal::EightInternalStates{T}) where {T<:IUI64}
	s1, s2, s3, s4, s5, s6, s7, s8 = state(internal)

	nextvalue =  s2
	nextvalue += s2 << 2
	nextvalue =  rotl(nextvalue, 7)
	nextvalue += nextvalue << 3

	t = s2 << 11

	s3 ⊻= s1
	s6 ⊻= s2
	s2 ⊻= s3
	s8 ⊻= s4
	s4 ⊻= s5
	s5 ⊻= s6
	s1 ⊻= s7
	s7 ⊻= s8

	s7 ⊻= t
	s8 = rotl(s8, 45)

	state(internal, s1, s2, s3, s4, s5, s6, s7, s8)

	return nextvalue
end


text:
using LightGraphs

function betweenness_centrality(c::Int, n::Int, k::Int)
    bc, elapsed = 0.0, 0.0
    for i in 1:c
        graph = barabasi_albert(n, k)
        elapsed += @elapsed begin
            bc += mean(LightGraphs.betweenness_centrality(graph))
            #  bc += mean(parallelbetweenness_centrality(graph))
        end
    end
    bc /= c
    elapsed /= c
    println("Betweenness Centrality:  $bc ($(elapsed)s/call)")
end

function clustering_coefficient(c::Int, n::Int, k::Int)
    cc, elapsed = 0.0, 0.0
    for i in 1:c
        graph = barabasi_albert(n, k)
        elapsed += @elapsed begin
            cc += mean(local_clustering_coefficient(graph))
        end
    end
    cc /= c
    elapsed /= c
    println("Average Clustering Coefficient:  $cc ($(elapsed)s/call)")
end

function main()
    n = parse(Int, ARGS[1])
    k = parse(Int, ARGS[2])
    if length(ARGS) >= 3
        srand(parse(Int, ARGS[3]))
    end

    betweenness_centrality(5, n, k)
    clustering_coefficient(10, n, k)
end

main()

max_stars_repo_name: ELIFE-ASU/meet-julia

max_issues_repo_name: ELIFE-ASU/meet-julia

max_forks_repo_name: ELIFE-ASU/meet-julia


text:
@with_kw mutable struct RewardModPOMDP{S, A, O} <: POMDP{S,A,O}
    pomdp::POMDP{S,A,O}
    reward_model
end

POMDPs.initialstate(mdp::RewardModPOMDP) = initialstate(mdp.pomdp)
POMDPs.initialobs(mdp::RewardModPOMDP, s) = initialobs(mdp.pomdp, s)

POMDPs.actions(mdp::RewardModPOMDP) = actions(mdp.pomdp)
POMDPs.actionindex(mdp::RewardModPOMDP, a) = actionindex(mdp.pomdp)

POMDPs.isterminal(mdp::RewardModPOMDP, s) = isterminal(mdp.pomdp, s)
POMDPs.discount(mdp::RewardModPOMDP) = discount(mdp.pomdp)

POMDPs.observation(mdp::RewardModPOMDP, s) = observation(mdp.pomdp, s)

render(mdp::RewardModPOMDP, s, a = nothing; kwargs...) = render(mdp.pomdp, s, a; kwargs...)
render(mdp::RewardModPOMDP; kwargs...) = render(mdp.pomdp; kwargs...)

function POMDPs.gen(mdp::RewardModPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; info=Dict())
    sp, o, r = gen(mdp.pomdp, s, a, rng, info=info)
    rnew = mdp.reward_model((s=s, a=a, sp=sp, o=o, r=r))
    return (sp=sp, o=o, r=rnew)
end

max_stars_repo_name: ancorso/POMDPGym

max_issues_repo_name: ancorso/POMDPGym

max_forks_repo_name: ancorso/POMDPGym


text:
"""
    implemented(function, Tuple{Arg1Type, Arg2Type})

Check whether there is an implementation available that will return a suitable value.
"""
implemented(f::Function, TT::TupleType) = hasmethod(f, TT)

"""
    @implemented function(::Arg1Type, ::Arg2Type)

Check whether there is an implementation available that will return a suitable value.
"""
macro implemented(ex)
    tplex = convert_req(ex)
    return quote
        implemented($(esc(tplex))...)
    end
end

"""
    get_requirements(f::Function, args::Tuple)

Return a RequirementSet for the function f and arguments args.
"""
get_requirements(f::Function, args::Tuple) = Unspecified((f, typeof(args)))

"""
    @get_requirements f(arg1, arg2)

Call get_requirements(f, (arg1,arg2)).
"""
macro get_requirements(call)
    return quote get_requirements($(esc(convert_call(call)))...) end
end

"""
    @POMDP_require solve(s::CoolSolver, p::POMDP) begin
        PType = typeof(p)
        @req states(::PType)
        @req actions(::PType)
        @req transition(::PType, ::S, ::A)
        s = first(states(p))
        a = first(actions(p))
        t_dist = transition(p, s, a)
        @req rand(::AbstractRNG, ::typeof(t_dist))
    end

Create a get_requirements implementation for the function signature and the requirements block.
"""
macro POMDP_require(typedcall, block)
    fname, args, types = unpack_typedcall(typedcall)
    tconstr = Expr[:($(Symbol(:T,i))<:$(esc(C))) for (i,C) in enumerate(types)] # oh snap
    ts = Symbol[Symbol(:T,i) for i in 1:length(types)]
    req_spec = :(($fname, Tuple{$(types...)}))
    fimpl = quote
        function POMDPs.get_requirements(f::typeof($(esc(fname))), args::Tuple{$(ts...)}) where {$(tconstr...)} # dang
            ($([esc(a) for a in args]...),) = args # whoah
            return $(pomdp_requirements(req_spec, block))
        end
    end
    return fimpl
end

"""
    reqs = @POMDP_requirements CoolSolver begin
        PType = typeof(p)
        @req states(::PType)
        @req actions(::PType)
        @req transition(::PType, ::S, ::A)
        s = first(states(p))
        a = first(actions(p))
        t_dist = transition(p, s, a)
        @req rand(::AbstractRNG, ::typeof(t_dist))
    end

Create a RequirementSet object.
"""
macro POMDP_requirements(name, block)
    return pomdp_requirements(name, block)
end

"""
    @warn_requirements solve(solver, problem)

Print a warning if there are unmet requirements.
"""
macro warn_requirements(call::Expr)
    quote
        reqs = get_requirements($(esc(convert_call(call)))...)
        if !check_requirements(reqs)
            show_requirements(reqs)
        end
    end
end

"""
    @show_requirements solve(solver, problem)

Print a a list of requirements for a function call.
"""
macro show_requirements(call::Expr)
    quote
        reqs = get_requirements($(esc(convert_call(call)))...)
        show_requirements(reqs)
    end
end

"""
    @requirements_info ASolver() [YourPOMDP()]

Print information about the requirements for a solver.
"""
macro requirements_info(exprs...)
    quote
        requirements_info($([esc(ex) for ex in exprs]...))
    end
end

"""
    requirements_info(s::Solver, p::Union{POMDP,MDP}, ...)

Print information about the requirement for solver s.
"""
function requirements_info(s::Union{Solver,Simulator})
    stype = typeof(s)
    try
        stype = stype.name.name
    catch ex
        # do nothing
    end
    println("""Please supply a POMDP as a second argument to requirements_info.
            e.g. `@requirements_info $(stype)() YourPOMDP()`
            """)
end
function requirements_info(s::Union{Solver,Simulator}, p::Union{POMDP,MDP}, args...)
    reqs = get_requirements(solve, (s, p))
    show_requirements(reqs)
end

"""
    @req f( ::T1, ::T2)

Convert a `f( ::T1, ::T2)` expression to a `(f, Tuple{T1,T2})::Req` for pushing to a `RequirementSet`.

If in a `@POMDP_requirements` or `@POMDP_require` block, marks the requirement for including in the set of requirements.
"""
macro req(ex)
    return esc(convert_req(ex))
end

"""
    @subreq f(arg1, arg2)

In a `@POMDP_requirements` or `@POMDP_require` block, include the requirements for `f(arg1, arg2)` as a child argument set.
"""
macro subreq(ex)
    return quote
        get_requirements($(esc(convert_call(ex)))...)
    end
end

"""
    check_requirements(r::AbstractRequirementSet)

Check whether the methods in `r` have implementations with `implemented()`. Return true if all methods have implementations.
"""
function check_requirements(r::AbstractRequirementSet)
    analyzed = Set()
    return recursively_check(r, analyzed)
end

"""
    show_requirements(r::AbstractRequirementSet)

Check whether the methods in `r` have implementations with `implemented()` and print out a formatted list showing which are missing. Return true if all methods have implementations.
"""
function show_requirements(r::AbstractRequirementSet)
    buf = stdout
    reported = Set{Req}()
    analyzed = Set()

    show_heading(buf, r.requirer)
    println(buf)

    allthere, first_exception = recursively_show(buf, r, analyzed, reported)

    if !allthere
        println("Note: Missing methods are often due to incorrect importing. You must explicitly import POMDPs functions to add new methods.")
        println()
    end

    if first_exception != nothing
        print("Throwing the first exception (from processing ")
        printstyled(handle_method(first_exception.requirer), color=:blue)
        println(" requirements):\n")
        rethrow(first_exception.exception)
    end

    return allthere
end


text:
"""
    LanguageServerInstance(pipe_in, pipe_out, env="", depot="", err_handler=nothing, symserver_store_path=nothing)

Construct an instance of the language server.

Once the instance is `run`, it will read JSON-RPC from `pipe_out` and
write JSON-RPC from `pipe_in` according to the [language server
specification](https://microsoft.github.io/language-server-protocol/specifications/specification-3-14/).
For normal usage, the language server can be instantiated with
`LanguageServerInstance(stdin, stdout, false, "/path/to/environment")`.

# Arguments
- `pipe_in::IO`: Pipe to read JSON-RPC from.
- `pipe_out::IO`: Pipe to write JSON-RPC to.
- `env::String`: Path to the
  [environment](https://docs.julialang.org/en/v1.2/manual/code-loading/#Environments-1)
  for which the language server is running. An empty string uses julia's
  default environment.
- `depot::String`: Sets the
  [`JULIA_DEPOT_PATH`](https://docs.julialang.org/en/v1.2/manual/environment-variables/#JULIA_DEPOT_PATH-1)
  where the language server looks for packages required in `env`.
- `err_handler::Union{Nothing,Function}`: If not `nothing`, catch all errors and pass them to an error handler
  function with signature `err_handler(err, bt)`. Mostly used for the VS Code crash reporting implementation.
- `symserver_store_path::Union{Nothing,String}`: if `nothing` is passed, the symbol server cash is stored in
  a folder in the package. If an absolute path is passed, the symbol server will store the cache files in that
  path. The path must exist on disc before this is called.
"""
mutable struct LanguageServerInstance
    jr_endpoint::Union{JSONRPC.JSONRPCEndpoint,Nothing}
    workspaceFolders::Set{String}
    _documents::Dict{URI2,Document}

    env_path::String
    depot_path::String
    symbol_server::SymbolServer.SymbolServerInstance
    symbol_results_channel::Channel{Any}
    global_env::StaticLint.ExternalEnv
    roots_env_map::Dict{Document,StaticLint.ExternalEnv}
    symbol_store_ready::Bool

    format_options::DocumentFormat.FormatOptions
    runlinter::Bool
    lint_options::StaticLint.LintOptions
    lint_missingrefs::Symbol
    lint_disableddirs::Vector{String}
    completion_mode::Symbol

    combined_msg_queue::Channel{Any}

    err_handler::Union{Nothing,Function}

    status::Symbol

    number_of_outstanding_symserver_requests::Int
    symserver_use_download::Bool

    current_symserver_progress_token::Union{Nothing,String}

    clientcapability_window_workdoneprogress::Bool
    clientcapability_workspace_didChangeConfiguration::Bool
    # Can probably drop the above 2 and use the below.
    clientCapabilities::Union{ClientCapabilities,Missing}
    clientInfo::Union{InfoParams,Missing}

    shutdown_requested::Bool

    function LanguageServerInstance(pipe_in, pipe_out, env_path="", depot_path="", err_handler=nothing, symserver_store_path=nothing, download=true)
        new(
            JSONRPC.JSONRPCEndpoint(pipe_in, pipe_out, err_handler),
            Set{String}(),
            Dict{URI2,Document}(),
            env_path,
            depot_path,
            SymbolServer.SymbolServerInstance(depot_path, symserver_store_path),
            Channel(Inf),
            StaticLint.ExternalEnv(deepcopy(SymbolServer.stdlibs), SymbolServer.collect_extended_methods(SymbolServer.stdlibs), collect(keys(SymbolServer.stdlibs))),
            Dict(),
            false,
            DocumentFormat.FormatOptions(),
            true,
            StaticLint.LintOptions(),
            :all,
            LINT_DIABLED_DIRS,
            :qualify, # options: :import or :qualify, anything else turns this off
            Channel{Any}(Inf),
            err_handler,
            :created,
            0,
            download,
            nothing,
            false,
            false,
            missing,
            missing,
            false
        )
    end
end
function Base.display(server::LanguageServerInstance)
    println("Root: ", server.workspaceFolders)
    for d in getdocuments_value(server)
        display(d)
    end
end

function hasdocument(server::LanguageServerInstance, uri::URI2)
    return haskey(server._documents, uri)
end

function getdocument(server::LanguageServerInstance, uri::URI2)
    return server._documents[uri]
end

function getdocuments_key(server::LanguageServerInstance)
    return keys(server._documents)
end

function getdocuments_pair(server::LanguageServerInstance)
    return pairs(server._documents)
end

function getdocuments_value(server::LanguageServerInstance)
    return values(server._documents)
end

function setdocument!(server::LanguageServerInstance, uri::URI2, doc::Document)
    server._documents[uri] = doc
end

function deletedocument!(server::LanguageServerInstance, uri::URI2)
    doc = getdocument(server, uri)
    StaticLint.clear_meta(getcst(doc))
    delete!(server._documents, uri)

    for d in getdocuments_value(server)
        if getroot(d) === doc
            setroot(d, d)
            semantic_pass(getroot(d))
        end
    end
end

function create_symserver_progress_ui(server)
    if server.clientcapability_window_workdoneprogress
        token = string(uuid4())
        server.current_symserver_progress_token = token
        JSONRPC.send(server.jr_endpoint, window_workDoneProgress_create_request_type, WorkDoneProgressCreateParams(token))

        JSONRPC.send(
            server.jr_endpoint,
            progress_notification_type,
            ProgressParams(token, WorkDoneProgressBegin("Julia Language Server", missing, "Indexing packages...", missing))
        )
    end
end

function destroy_symserver_progress_ui(server)
    if server.clientcapability_window_workdoneprogress && server.current_symserver_progress_token !== nothing
        progress_token = server.current_symserver_progress_token
        server.current_symserver_progress_token = nothing
        JSONRPC.send(
            server.jr_endpoint,
            progress_notification_type,
            ProgressParams(progress_token, WorkDoneProgressEnd(missing))
        )
    end
end

function trigger_symbolstore_reload(server::LanguageServerInstance)
    server.symbol_store_ready = false
    if server.number_of_outstanding_symserver_requests == 0 && server.status == :running
        create_symserver_progress_ui(server)
    end
    server.number_of_outstanding_symserver_requests += 1

    if server.symserver_use_download
        @debug "Will download symbol server caches for this instance."
    end

    @async try
        # TODO Add try catch handler that links into crash reporting
        ssi_ret, payload = SymbolServer.getstore(
            server.symbol_server,
            server.env_path,
            function (i)
                if server.clientcapability_window_workdoneprogress && server.current_symserver_progress_token !== nothing
                    JSONRPC.send(
                        server.jr_endpoint,
                        progress_notification_type,
                        ProgressParams(server.current_symserver_progress_token, WorkDoneProgressReport(missing, "Indexing $i...", missing))
                    )
                else
                    @info "Indexing $i..."
                end
            end,
            server.err_handler,
            download = server.symserver_use_download
        )

        server.number_of_outstanding_symserver_requests -= 1

        if server.number_of_outstanding_symserver_requests == 0
            destroy_symserver_progress_ui(server)
        end

        if ssi_ret == :success
            push!(server.symbol_results_channel, payload)
        elseif ssi_ret == :failure
            error_payload = Dict(
                "command" => "symserv_crash",
                "name" => "LSSymbolServerFailure",
                "message" => payload === nothing ? "" : String(take!(payload)),
                "stacktrace" => "")
            JSONRPC.send(
                server.jr_endpoint,
                telemetry_event_notification_type,
                error_payload
            )
        elseif ssi_ret == :package_load_crash
            error_payload = Dict(
                "command" => "symserv_pkgload_crash",
                "name" => payload.package_name,
                "message" => payload.stderr === nothing ? "" : String(take!(payload.stderr)))
            JSONRPC.send(
                server.jr_endpoint,
                telemetry_event_notification_type,
                error_payload
            )
        end
        server.symbol_store_ready = true
    catch err
        bt = catch_backtrace()
        if server.err_handler !== nothing
            server.err_handler(err, bt)
        else
            Base.display_error(stderr, err, bt)
        end
    end
end

function request_wrapper(func, server::LanguageServerInstance)
    return function (conn, params)
        if server.shutdown_requested
            # it's fine to always return a value here, even for notifications, because
            # JSONRPC discards it anyways in that case
            return JSONRPC.JSONRPCError(
                -32600,
                "LS shutdown was requested.",
                nothing
            )
        end
        func(params, server, conn)
    end
end

"""
    run(server::LanguageServerInstance)

Run the language `server`.
"""
function Base.run(server::LanguageServerInstance)
    server.status = :started

    run(server.jr_endpoint)

    trigger_symbolstore_reload(server)

    @async try
        @debug "LS: Starting client listener task."
        while true
            msg = JSONRPC.get_next_message(server.jr_endpoint)
            put!(server.combined_msg_queue, (type = :clientmsg, msg = msg))
        end
    catch err
        bt = catch_backtrace()
        if server.err_handler !== nothing
            server.err_handler(err, bt)
        else
            io = IOBuffer()
            Base.display_error(io, err, bt)
            print(stderr, String(take!(io)))
        end
    finally
        if isopen(server.combined_msg_queue)
            put!(server.combined_msg_queue, (type = :close,))
            close(server.combined_msg_queue)
        end
        @debug "LS: Client listener task done."
    end

    @async try
        @debug "LS: Starting symbol server listener task."
        while true
            msg = take!(server.symbol_results_channel)
            put!(server.combined_msg_queue, (type = :symservmsg, msg = msg))
        end
    catch err
        bt = catch_backtrace()
        if server.err_handler !== nothing
            server.err_handler(err, bt)
        else
            io = IOBuffer()
            Base.display_error(io, err, bt)
            print(stderr, String(take!(io)))
        end
    finally
        if isopen(server.combined_msg_queue)
            put!(server.combined_msg_queue, (type = :close,))
            close(server.combined_msg_queue)
        end
        @debug "LS: Symbol server listener task done."
    end

    msg_dispatcher = JSONRPC.MsgDispatcher()

    msg_dispatcher[textDocument_codeAction_request_type] = request_wrapper(textDocument_codeAction_request, server)
    msg_dispatcher[workspace_executeCommand_request_type] = request_wrapper(workspace_executeCommand_request, server)
    msg_dispatcher[textDocument_completion_request_type] = request_wrapper(textDocument_completion_request, server)
    msg_dispatcher[textDocument_signatureHelp_request_type] = request_wrapper(textDocument_signatureHelp_request, server)
    msg_dispatcher[textDocument_definition_request_type] = request_wrapper(textDocument_definition_request, server)
    msg_dispatcher[textDocument_formatting_request_type] = request_wrapper(textDocument_formatting_request, server)
    msg_dispatcher[textDocument_references_request_type] = request_wrapper(textDocument_references_request, server)
    msg_dispatcher[textDocument_rename_request_type] = request_wrapper(textDocument_rename_request, server)
    msg_dispatcher[textDocument_prepareRename_request_type] = request_wrapper(textDocument_prepareRename_request, server)
    msg_dispatcher[textDocument_documentSymbol_request_type] = request_wrapper(textDocument_documentSymbol_request, server)
    msg_dispatcher[textDocument_documentHighlight_request_type] = request_wrapper(textDocument_documentHighlight_request, server)
    msg_dispatcher[textDocument_semanticTokens_request_type] = request_wrapper(textDocument_semanticTokens_request, server)
    msg_dispatcher[textDocument_semanticTokens_full_request_type] = request_wrapper(textDocument_semanticTokens_full_request, server)
    msg_dispatcher[julia_getModuleAt_request_type] = request_wrapper(julia_getModuleAt_request, server)
    msg_dispatcher[julia_getDocAt_request_type] = request_wrapper(julia_getDocAt_request, server)
    msg_dispatcher[textDocument_hover_request_type] = request_wrapper(textDocument_hover_request, server)
    msg_dispatcher[initialize_request_type] = request_wrapper(initialize_request, server)
    msg_dispatcher[initialized_notification_type] = request_wrapper(initialized_notification, server)
    msg_dispatcher[shutdown_request_type] = request_wrapper(shutdown_request, server)
    msg_dispatcher[exit_notification_type] = request_wrapper(exit_notification, server)
    msg_dispatcher[cancel_notification_type] = request_wrapper(cancel_notification, server)
    msg_dispatcher[setTrace_notification_type] = request_wrapper(setTrace_notification, server)
    msg_dispatcher[setTraceNotification_notification_type] = request_wrapper(setTraceNotification_notification, server)
    msg_dispatcher[julia_getCurrentBlockRange_request_type] = request_wrapper(julia_getCurrentBlockRange_request, server)
    msg_dispatcher[julia_activateenvironment_notification_type] = request_wrapper(julia_activateenvironment_notification, server)
    msg_dispatcher[textDocument_didOpen_notification_type] = request_wrapper(textDocument_didOpen_notification, server)
    msg_dispatcher[textDocument_didClose_notification_type] = request_wrapper(textDocument_didClose_notification, server)
    msg_dispatcher[textDocument_didSave_notification_type] = request_wrapper(textDocument_didSave_notification, server)
    msg_dispatcher[textDocument_willSave_notification_type] = request_wrapper(textDocument_willSave_notification, server)
    msg_dispatcher[textDocument_willSaveWaitUntil_request_type] = request_wrapper(textDocument_willSaveWaitUntil_request, server)
    msg_dispatcher[textDocument_didChange_notification_type] = request_wrapper(textDocument_didChange_notification, server)
    msg_dispatcher[workspace_didChangeWatchedFiles_notification_type] = request_wrapper(workspace_didChangeWatchedFiles_notification, server)
    msg_dispatcher[workspace_didChangeConfiguration_notification_type] = request_wrapper(workspace_didChangeConfiguration_notification, server)
    msg_dispatcher[workspace_didChangeWorkspaceFolders_notification_type] = request_wrapper(workspace_didChangeWorkspaceFolders_notification, server)
    msg_dispatcher[workspace_symbol_request_type] = request_wrapper(workspace_symbol_request, server)
    msg_dispatcher[julia_refreshLanguageServer_notification_type] = request_wrapper(julia_refreshLanguageServer_notification, server)
    msg_dispatcher[julia_getDocFromWord_request_type] = request_wrapper(julia_getDocFromWord_request, server)
    msg_dispatcher[textDocument_selectionRange_request_type] = request_wrapper(textDocument_selectionRange_request, server)

    while true
        message = take!(server.combined_msg_queue)

        if message.type == :close
            @info "Shutting down server instance."
            return
        elseif message.type == :clientmsg
            msg = message.msg
            JSONRPC.dispatch_msg(server.jr_endpoint, msg_dispatcher, msg)
        elseif message.type == :symservmsg
            @info "Received new data from Julia Symbol Server."

            server.global_env.symbols = message.msg
            server.global_env.extended_methods = SymbolServer.collect_extended_methods(server.global_env.symbols)
            server.global_env.project_deps = collect(keys(server.global_env.symbols))

            # redo roots_env_map
            for (root, _) in server.roots_env_map
                @debug "resetting get_env_for_root"
                newenv = get_env_for_root(root, server)
                if newenv === nothing
                    delete!(server.roots_env_map, root)
                else
                    server.roots_env_map[root] = newenv
                end
            end

            @debug "starting re-lint of everything"
            relintserver(server)
            @debug "re-lint done"
        end
    end
end

function relintserver(server)
    roots = Set{Document}()
    documents = getdocuments_value(server)
    for doc in documents
        StaticLint.clear_meta(getcst(doc))
        set_doc(getcst(doc), doc)
    end
    for doc in documents
        # only do a pass on documents once
        root = getroot(doc)
        if !(root in roots)
            push!(roots, root)
            semantic_pass(root)
        end
    end
    for doc in documents
        lint!(doc, server)
    end
end


text:
using Requires

struct D5DW_Domainwall_operator{Dim,T,fermion,wilsonfermion} <: Dirac_operator{Dim}   where  T <: AbstractGaugefields
    U::Array{T,1}
    wilsonoperator::Union{Wilson_Dirac_operator{Dim,T,wilsonfermion},Wilson_Dirac_operator_faster{Dim,T,wilsonfermion} }
    mass::Float64
    _temporary_fermi::Array{fermion,1}
    L5::Int64
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    verbose_print::Verbose_print
    _temporary_fermion_forCG::Vector{fermion}
    boundarycondition::Vector{Int8}

end

function D5DW_Domainwall_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters,mass) where  {NC,Dim}
    @assert haskey(parameters,"L5") "parameters should have the keyword L5"
    L5 = parameters["L5"]
    if L5 != x.L5
        @warn "L5 in Dirac operator and fermion fields is not same. Now L5 = $L5 and x.L5 = $(x.L5)"
        #@assert L5 == x.L5 "L5 in Dirac operator and fermion fields should be same. Now L5 = $L5 and x.L5 = $(x.L5)"
    end

    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end

    #boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    T = eltype(U)

    #Parameters for Wilson operator
        r = 1
        M = check_parameters(parameters,"M",-1)
        #@assert haskey(parameters,"M") "parameters should have the keyword M"
        #M = parameters["M"]
        κ_wilson = 1/(2*Dim*r+2M)
        parameters_wilson = Dict()
        parameters_wilson["κ"] = κ_wilson
        parameters_wilson["numtempvec"] = 4
        parameters_wilson["numtempvec_CG"] = 1
        parameters_wilson["boundarycondition"] = boundarycondition
        x_wilson  = x.w[1]
        fasterversion = check_parameters(parameters,"faster version",false)
        if fasterversion
            wilsonoperator = Wilson_Dirac_operator_faster(U,x_wilson,parameters_wilson)
        else
            wilsonoperator = Wilson_Dirac_operator(U,x_wilson,parameters_wilson)
        end
    #--------------------------------

    xtype = typeof(x)
    num = 2
    _temporary_fermi = Array{xtype,1}(undef,num)
    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    #println("eps_CG = ",eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)
    verbose_print = Verbose_print(verbose_level)

    method_CG = check_parameters(parameters,"method_CG","bicg")
    #println("xtype ",xtype)
    #println("D ", typeof(wilsonoperator))

    return D5DW_Domainwall_operator{Dim,T,xtype,typeof(x_wilson)}(
        U,
        wilsonoperator,
        mass,
        _temporary_fermi,
        L5,
        eps_CG,
        MaxCGstep,
        verbose_level,
        method_CG,
        verbose_print,
        _temporary_fermion_forCG,
        boundarycondition
        )
end

function (D::D5DW_Domainwall_operator{Dim,T,fermion,wilsonfermion})(U) where {Dim,T,fermion,wilsonfermion}
    return D5DW_Domainwall_operator{Dim,T,fermion,wilsonfermion}(
        U,
        D.wilsonoperator(U),
        D.mass,
        D._temporary_fermi,
        D.L5,
        D.eps_CG,
        D.MaxCGstep,
        D.verbose_level,
        D.method_CG,
        D.verbose_print,
        D._temporary_fermion_forCG,
        D.boundarycondition
    )
end

struct Adjoint_D5DW_Domainwall_operator{T} <: Adjoint_Dirac_operator
    parent::T
end

struct Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion} <: Dirac_operator{Dim}  where T <: AbstractGaugefields
    U::Array{T,1}
    D5DW::D5DW_Domainwall_operator{Dim,T,fermion,wilsonfermion}
    D5DW_PV::D5DW_Domainwall_operator{Dim,T,fermion,wilsonfermion}
    mass::Float64
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    verbose_print::Verbose_print
    boundarycondition::Vector{Int8}

end

function Domainwall_Dirac_operator(U::Array{<:AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
    @assert haskey(parameters,"mass") "parameters should have the keyword mass"
    mass = parameters["mass"]
    D5DW = D5DW_Domainwall_operator(U,x,parameters,mass)
    D5DW_PV = D5DW_Domainwall_operator(U,x,parameters,1)

    #boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])

    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    #println("eps_CG = ",eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)
    verbose_print = Verbose_print(verbose_level)

    method_CG = check_parameters(parameters,"method_CG","bicg")

    return Domainwall_Dirac_operator{Dim,eltype(U),typeof(x),typeof(x.w[1])}(
            U,
            D5DW,D5DW_PV,mass,
            eps_CG,
            MaxCGstep,
            verbose_level,
            method_CG,
            verbose_print,
            boundarycondition
    )
end

function (D::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion})(U) where {Dim,T,fermion,wilsonfermion}
    return Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}(
        U,
        D.D5DW(U),D.D5DW_PV(U),
        D.mass,
        D.eps_CG,
        D.MaxCGstep,
        D.verbose_level,
        D.method_CG,
        D.verbose_print,
        D.boundarycondition
    )
end

function get_temporaryvectors(A::T) where T <: Domainwall_Dirac_operator
    return A.D5DW._temporary_fermi
end

function get_temporaryvectors_forCG(A::T) where T <: Domainwall_Dirac_operator
    return A.D5DW._temporary_fermion_forCG
end

struct Adjoint_Domainwall_operator{Dim,T,fermion,wilsonfermion} <: Adjoint_Dirac_operator
    parent::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}
end

function Base.adjoint(A::D5DW_Domainwall_operator)
    Adjoint_D5DW_Domainwall_operator(A)
end

function Base.adjoint(A::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
    Adjoint_Domainwall_operator{Dim,T,fermion,wilsonfermion}(A)
end

struct DdagD_Domainwall_operator{Dim,T,fermion,wilsonfermion} <: DdagD_operator
    dirac::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}
    function DdagD_Domainwall_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
        return new{Dim,eltype(U),typeof(x),typeof(x.w[1])}(Domainwall_Dirac_operator(U,x,parameters))
    end

    function DdagD_Domainwall_operator(D::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
        return new{Dim,T,fermion,wilsonfermion}(D)
    end
end

struct D5DWdagD5DW_Wilson_operator{T} <: DdagD_operator
    dirac::T
    function D5DWdagD5DW_Wilson_operator(U::Array{T,1},x,parameters,mass) where  T <: AbstractGaugefields
        return new{T}(D5DW_Domainwall_operator(U,x,parameters,mass))
    end

    function D5DWdagD5DW_Wilson_operator(D::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
        dtype = typeof(D.D5DW)
        return new{dtype}(D.D5DW)
    end

    function D5DWdagD5DW_Wilson_operator(D::D5DW_Domainwall_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
        dtype = typeof(D)
        return new{dtype}(D)
    end
end

struct Reflection_matrixJ
end

struct Permutation_matrixP
end

struct Adjoint_Permutation_matrixP
end

function Base.adjoint(A::Permutation_matrixP)
    return Adjoint_Permutation_matrixP()
end

function Base.adjoint(A::Adjoint_Permutation_matrixP)
    return Permutation_matrixP()
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Reflection_matrixJ, T3 <:AbstractFermionfields}
    apply_J!(y,x)
    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Permutation_matrixP, T3 <:AbstractFermionfields}
    apply_P!(y,x)
    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_Permutation_matrixP, T3 <:AbstractFermionfields}
    apply_Pdag!(y,x)
    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: D5DWdagD5DW_Wilson_operator, T3 <:AbstractFermionfields}
    temp = A.dirac._temporary_fermi[1]
    mul!(temp,A.dirac,x)
    mul!(y,A.dirac',temp)

    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: D5DW_Domainwall_operator, T3 <:AbstractFermionfields}
    D5DWx!(y,A.U,x,A.mass,A.wilsonoperator,A.L5)
    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_D5DW_Domainwall_operator, T3 <:AbstractFermionfields}
    D5DWdagx!(y,A.parent.U,x,A.parent.mass,A.parent.wilsonoperator,A.parent.L5)
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_Domainwall_operator, T3 <:AbstractFermionfields}

    #A = D5DW(m)*D5DW(m=1)^(-1)
    #=
    A^+ = [D5DW(m)*D5DW(m=1)^(-1)]^+
        = D5DW(m=1)^(-1)^+ D5DW(m)^+
        = [D5DW(m=1)^+]^(-1) D5DW(m)^+
    y = A^+*x = [D5DW(m=1)^+]^(-1) D5DW(m)^+*x
    =#
    mul!(A.parent.D5DW_PV._temporary_fermi[1],A.parent.D5DW',x)
    bicg(y,A.parent.D5DW_PV',A.parent.D5DW_PV._temporary_fermi[1],
        eps=A.parent.eps_CG,maxsteps = A.parent.MaxCGstep,verbose = A.parent.verbose_print)
        #maxsteps = 10000)
        #verbose = Verbose_3())

    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Domainwall_Dirac_operator, T3 <:AbstractFermionfields}
    #A = D5DW(m)*D5DW(m=1))^(-1)
    #y = A*x = D5DW(m)*D5DW(m=1))^(-1)*x
    #println("x ",x.w[1][1,1,1,1])
    bicg(A.D5DW_PV._temporary_fermi[1],A.D5DW_PV,x,
        eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)
    mul!(y,A.D5DW,A.D5DW_PV._temporary_fermi[1])
    #println("y ",y.w[1][1,1,1,1])

    #bicg(A.D5DW._temporary_fermi[1],A.D5DW,y,
    #eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)
    #mul!(A.D5DW._temporary_fermi[2],A.D5DW_PV,A.D5DW._temporary_fermi[1])
    #println("tmp2 ",A.D5DW._temporary_fermi[2].w[1][1,1,1,1])

    #error("Do not use Domainwall_operator directory. Use D5DW_Domainwall_operator M = D5DW(m)*D5DW(-1)^{-1}")
    #D5DWx!(y,A.U,x,A.m,A.wilsonoperator._temporary_fermi)
    return
end

include("./DomainwallFermion_5d.jl")
include("./DomainwallFermion_5d_wing.jl")
include("./DomainwallFermion_3d_wing.jl")

function Initialize_DomainwallFermion(u::AbstractGaugefields{NC,Dim},L5;nowing=false) where {NC,Dim}
    _,_,NN... = size(u)
    return Initialize_DomainwallFermion(L5,NC,NN...,nowing = nowing)
end

function Initialize_DomainwallFermion(L5,NC,NN...;nowing = false)
    Dim = length(NN)
    if Dim == 4
        if nowing
            fermion = DomainwallFermion_5D(L5,NC,NN...,nowing = nowing)
        else
            fermion = DomainwallFermion_5D_wing(L5,NC,NN...)
        end
        #fermion = DomainwallFermion_5D_wing(L5,NC,NN...)
        #fermion = WilsonFermion_4D_wing{NC}(NN...)
        #fermion = WilsonFermion_4D_wing(NC,NN...)
    elseif Dim == 2
        fermion = DomainwallFermion_3D_wing(L5,NC,NN...)
        #fermion = WilsonFermion_2D_wing{NC}(NN...)
    else
        error("Dimension $Dim is not supported")
    end
    return fermion
end

function bicg(x,A::Domainwall_Dirac_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1))^(-1)
    #A^-1 = D5DW(m=1)*DsDW(m)^-1
    #x = A^-1*b = D5DW(m=1)*DsDW(m)^-1*b
    #println("b ",b.w[1][1,1,1,1])
    bicg(A.D5DW._temporary_fermi[1],A.D5DW,b;eps=eps,maxsteps = maxsteps,verbose = verbose)

    #mul!(x,A.D5DW,A.D5DW._temporary_fermi[1])
    #println("x ",x.w[1][1,1,1,1])

    mul!(x,A.D5DW_PV,A.D5DW._temporary_fermi[1])
end

function bicg(x,A::Adjoint_Domainwall_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) where T #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1)^(-1)
    #A' = (D5DW(m=1)^+)^(-1) D5DW(m)^+
    #A'^-1 = D5DW(m)^+^-1 D5DW(m=1)^+
    #x = A'^-1*b =D5DW(m)^+^-1 D5DW(m=1)^+*b
    mul!(A.parent.D5DW._temporary_fermi[1],A.parent.D5DW_PV',b)
    bicg(x,A.parent.D5DW',A.parent.D5DW._temporary_fermi[1];eps=eps,maxsteps = maxsteps,verbose = verbose)

end

function bicgstab(x,A::Domainwall_Dirac_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1))^(-1)
    #A^-1 = D5DW(m=1)*DsDW(m)^-1
    #x = A^-1*b = D5DW(m=1)*DsDW(m)^-1*b

    bicgstab(A.D5DW._temporary_fermi[1],A.D5DW,b;eps=eps,maxsteps = maxsteps,verbose = verbose)
    mul!(x,A.D5DW_PV,A.D5DW._temporary_fermi[1])
end

function bicgstab(x,A::Adjoint_Domainwall_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) where T #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1)^(-1)
    #A' = (D5DW(m=1)^+)^(-1) D5DW(m)^+
    #A'^-1 = D5DW(m)^+^-1 D5DW(m=1)^+
    #x = A'^-1*b =D5DW(m)^+^-1 D5DW(m=1)^+*b
    mul!(A.parent.D5DW._temporary_fermi[1],A.parent.D5DW_PV',b)
    bicgstab(x,A.parent.D5DW',A.parent.D5DW._temporary_fermi[1];eps=eps,maxsteps = maxsteps,verbose = verbose)

end

function cg(x,A::DdagD_Domainwall_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2))
    #=
    A^-1 = ( (D5DW(m)*D5DW(m=1))^(-1))^+ D5DW(m)*D5DW(m=1))^(-1) )^-1
      = ( D5DW(m=1)^+)^(-1) D5DW(m)^+ D5DW(m)*D5DW(m=1))^(-1) )^-1
      = D5DW(m=1) (  D5DW(m)^+ D5DW(m) )^(-1) )^-1 D5DW(m=1)^+
    x = A^-1*b = D5DW(m=1) (  D5DW(m)^+ D5DW(m) )^(-1)  D5DW(m=1)^+*b
    =#
    mul!(A.dirac.D5DW_PV._temporary_fermi[1],A.dirac.D5DW_PV',b) #D5DW(m=1)^+*b

    temp = A.dirac.D5DW_PV._temporary_fermi[1]
    cg(A.dirac.D5DW._temporary_fermi[1],A.DdagD,temp;eps=eps,maxsteps = maxsteps,verbose = verbose) #(  D5DW(m)^+ D5DW(m) )^(-1)  D5DW(m=1)^+*b
    mul!(x,A.dirac.D5DW_PV,A.dirac.D5DW._temporary_fermi[1])
end


text:
immutable WeightDim
    rows::Int
    cols::Int
    size::Int

    WeightDim(rows, cols) = new(rows, cols, rows * cols)
end

immutable WeightPos{S}
    seq::S
    dims::Tuple{DimT, DimT}
end

WeightPos{S}(seq::S, dim1, dim2) = WeightPos{S}(seq, (dim1, dim2))

Base.show(io::IO, pos::WeightPos) = print(io, "< Size: $(pos.dims), Seq: $((seqBegin(pos.seq),seqEnd(pos.seq))) >")

immutable MLPData{C, A}
    ctx::C
    numLayers::Int
    weightCount::Int
    empty::A
    weightsArray::A
    signal::Vector{A}
    weights::Vector{WeightPos}
end

function MLPData{C}(ctx::C, weightDims::Vector{WeightDim})
    numLayers = length(weightDims) + 1
    empty = array(ctx)
    A = typeof(empty)
    weights = Vector{WeightPos}()
    signal = Vector{A}()
    weightCount = map(d -> d.size, weightDims) |> sum
    weightsArray = array(ctx, Float32, weightCount)
    wIdx = 0
    for i in 1:numLayers
        push!(signal, array(ctx))
        if i != numLayers
            dims = weightDims[i]
            push!(weights, WeightPos(seq(ctx, wIdx, wIdx + dims.size - 1), dims.rows, dims.cols))
            wIdx += dims.size
        end
    end
    MLPData{C, A}(ctx, numLayers, weightCount, empty, weightsArray, signal, weights)
end

Base.show(io::IO, mlp::MLPData) = print(io, "Weight Count: $(mlp.weightCount)\nWeight Layout: $(mlp.weights)")

function release!(mlp::MLPData)
    release!(mlp.weightsArray)
    for s in mlp.signal
        release!(s)
    end
end

function addBias!(mlp::MLPData, input)
    joinArrays(1, constant(mlp.ctx, 1.0f0, dims(input, 0)), input)
end

function predict!(impl, mlp::MLPData, input)
    mlp.signal[1][] = input
    const eof = mlp.numLayers - 1
    for i in 1:eof
        scope!(mlp.ctx) do this
            inVec = makeInputVector(impl, i)
            weights = mlp.weights[i]
            currWeights = moddims(mlp.weightsArray[weights.seq], weights.dims...)
            outVec = matmul(inVec, currWeights)
            if i == eof
                mlp.signal[i + 1][] = sigmoid(outVec)
            else
                mlp.signal[i + 1][] = sin(outVec ./ 4f0)
            end
        end
    end
    mlp.signal[mlp.numLayers]
end

function setWeights!(mlp::MLPData, weights)
    @assert dims(weights) == dims(mlp.weightsArray)
    mlp.weightsArray[] = weights
end

reset!(mlp::MLPData) = for s in mlp.signal
    s[] = mlp.empty
end

arrayType{C, A}(mlp::MLPData{C, A}) = A

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using ForwardDiff
using ForwardDiff: Dual

seed(x::Real, ::Val) = Dual(x, true)

function seed(x, ::Val{N}, offset = 0) where N
  map(x, reshape(1:length(x), size(x))) do x, i
    Dual(x, ntuple(j -> j+offset == i, Val(N)))
  end
end

extract(x::ForwardDiff.Dual) = x.value, [x.partials...]

function extract(xs::AbstractArray{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
  J = similar(xs, V, N, length(xs))
  for i = 1:length(xs), j = 1:N
    J[j, i] = xs[i].partials.values[j]
  end
  return map(x -> x.value, xs), J
end

function forward_jacobian(f, x, ::Val{N}) where N
  y, _J = extract(f(seed(x, Val(N))))
  J = similar(_J, length(x), length(y))
  J[1:N,:] = _J
  offset = 0
  while offset + N < length(x)
    offset += N
    _, _J = extract(f(seed(x, Val(N), offset)))
    range = (1+offset):min(N+offset,length(x))
    J[range,:] = @view _J[range.-offset,:]
  end
  return y, J
end

function forward_jacobian(f, x)
  len = min(length(x), ForwardDiff.DEFAULT_CHUNK_THRESHOLD)
  forward_jacobian(f, x, Val(len))
end

vec_scalar(x) = vec(x)
vec_scalar(x::Real) = [x]
reshape_scalar(x, y) = reshape(y, size(x))
reshape_scalar(x::Real, y) = y[]

"""
    forwarddiff(f, x) -> f(x)

Runs `f(x)` as usual, but instructs Zygote to differentiate `f` using forward
mode, rather than the usual reverse mode.

Forward mode takes time linear in `length(x)` but only has constant memory
overhead, and is very efficient for scalars, so in some cases this can be a
useful optimisation.

```julia
julia> function pow(x, n)
         r = one(x)
         for i = 1:n
           r *= x
         end
         return r
       end
pow (generic function with 1 method)

julia> gradient(5) do x
         forwarddiff(x) do x
           pow(x, 2)
         end
       end
(10,)
```

Note that the function `f` will *drop gradients* for any closed-over values.

```julia
julia> gradient(2, 3) do a, b
         forwarddiff(a) do a
           a*b
         end
       end
(3, nothing)
```

This can be rewritten by explicitly passing through `b`, i.e.

```julia
gradient(2, 3) do a, b
  forwarddiff([a, b]) do (a, b)
    a*b
  end
end
```
"""
forwarddiff(f, x) = f(x)

@adjoint function forwarddiff(f, x)
  y, J = forward_jacobian(f, x)
  return y, ȳ -> (nothing, reshape_scalar(x, J*vec_scalar(ȳ)))
end


text:
for i in 1:10

    a = now()
    b = now()
    println(i)
    sleep(3)

end

DateTime(time())

DateTime(2018, 9, 9, 14, 0, 0)

a = now()

b = now()

println()

floor(DateTime(now()), Dates.second)

second(now())

second(now()) + 10

sleep(10)

DateTime(Year(now()), Month(now()), , 1, 1, 1)

Time(now())

# if endTime kleiner als now, ergebnis ist Null

x = true

startTime = Time(Hour(now()), minute(now()), second(now()))
endTime = Time(Hour(now()), minute(now()), second(now()))

d = endTime - startTime

second(DateTime(d))

a = Time(now())
b = Time(now())

Dates.value(b-a)

Int64(round(Float64(b-a)/1000000000))

Int64(round((b-a)/1000000000))

while x == true

    if endTime < Time(now())
        println("Zeit vorbei du looser!")

        x = false
    end

end

using Dates

while endTime < Time(now())

    if typeof(a) < String a = readline() end

end

function Eingabe() x = readline() end

x = true
endTime = Time(Hour(now()), minute(now()), second(now()) + 10)
while endTime > Time(now())

    readline()

end

Dates.value(Time(now()))


text:
function ClassicalRungeKuttaIntegrator(arg0::jdouble)
    return ClassicalRungeKuttaIntegrator((jdouble,), arg0)
end

function get_a(obj::ClassicalRungeKuttaIntegrator)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::ClassicalRungeKuttaIntegrator)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::ClassicalRungeKuttaIntegrator)
    return jcall(obj, "getC", Vector{jdouble}, ())
end

function integrate(obj::RungeKuttaIntegrator, arg0::ExpandableODE, arg1::ODEState, arg2::jdouble)
    return jcall(obj, "integrate", ODEStateAndDerivative, (ExpandableODE, ODEState, jdouble), arg0, arg1, arg2)
end

function single_step(obj::RungeKuttaIntegrator, arg0::OrdinaryDifferentialEquation, arg1::jdouble, arg2::Vector{jdouble}, arg3::jdouble)
    return jcall(obj, "singleStep", Vector{jdouble}, (OrdinaryDifferentialEquation, jdouble, Vector{jdouble}, jdouble), arg0, arg1, arg2, arg3)
end


text:
module TreeTools

using Random

## Iteration
import Base: eltype, iterate, IteratorEltype, IteratorSize, length
## Indexing
import Base: eachindex, firstindex, get!, getindex, lastindex, setindex!

## Others
import Base: ==, cat, convert, copy, count, hash, in, intersect
import Base: isempty, isequal, keys, map!, setdiff, show, size
import Base: union, union!, unique, unique!

##
include("objects.jl")
export Tree, TreeNode, TreeNodeData, MiscData
export isleaf, isroot, ancestor, children, branch_length, label

include("methods.jl")
export lca, node2tree, node2tree!, node_depth, distance, divtime, share_labels

include("iterators.jl")
export POT, POTleaves, nodes, leaves, internals

include("prunegraft.jl")
export delete_node!, graftnode!, prunenode!, prunenode, prunesubtree!

include("reading.jl")
export parse_newick_string, read_tree

include("writing.jl")
export write_newick, write_fasta

include("misc.jl")
export print_tree, check_tree, print_tree_ascii

include("splits.jl")
export Split, SplitList
export arecompatible, iscompatible

end


text:
using DSGE, Dates
path = dirname(@__FILE__)

mb_empty =  MeansBands()
@show mb_empty
@test isempty(mb_empty)
mb_full = load("$path/../reference/MeansBands.jld2", "mb")

@testset "Test that you can construct MeansBands objects and do stuff with them" begin
     # If one is empty, just return the non-empty
    @test cat(mb_empty, mb_full).means == mb_full.means
    @test cat(mb_full, mb_empty).means == mb_full.means
    @test get_class(mb_full) == :obs
    @test get_class(mb_empty) == :none
    @test get_product(mb_full) == :forecast
    @test get_product(mb_empty) == :none
    @test DSGE.get_cond_type(mb_full) == :none
    @test DSGE.get_cond_type(mb_empty) == :none
    @test DSGE.get_para(mb_full) == :full
    @test DSGE.get_para(mb_empty) == :none
    @test DSGE.n_vars_means(mb_full) == 3
    @test DSGE.n_vars_means(mb_empty) == 1
    @test DSGE.get_vars_means(mb_empty) == [:none]
    @test DSGE.get_vars_means(mb_full) == [:obs_gdp, :obs_cpi, :obs_nominalrate]
    @test DSGE.n_periods_means(mb_full) == 60
    @test DSGE.n_periods_means(mb_empty) == 1
    @test DSGE.startdate_means(mb_full) == quartertodate("2015-Q4")
    @test DSGE.startdate_means(mb_empty) == Dates.Date(0000, 1, 1)
    @test DSGE.enddate_means(mb_full) == quartertodate("2030-Q3")
    @test DSGE.enddate_means(mb_empty) == Dates.Date(0000, 1, 1)

    @test DSGE.n_vars_bands(mb_full) == 3
    @test DSGE.n_vars_bands(mb_empty) == 1
    @test DSGE.n_periods_bands(mb_full) == 60
    @test DSGE.n_periods_bands(mb_empty) == 1
    @test DSGE.startdate_bands(mb_full) == quartertodate("2015-Q4")
    @test DSGE.startdate_bands(mb_empty) == Dates.Date(0000, 1, 1)
    @test DSGE.enddate_bands(mb_full) == quartertodate("2030-Q3")
    @test DSGE.enddate_bands(mb_empty) == Dates.Date(0000, 1, 1)

    # Can only call on shockdec or irf MeansBands objects
    @test_throws AssertionError DSGE.get_shocks(mb_empty)
    @test_throws AssertionError DSGE.get_variables(mb_empty)
    @test_throws AssertionError DSGE.get_scenario_key(mb_empty)

    @test which_density_bands(mb_full) == ["90.0% LB", "80.0% LB", "70.0% LB", "60.0% LB", "50.0% LB", "50.0% UB", "60.0% UB", "70.0% UB", "80.0% UB", "90.0% UB"]
    @test which_density_bands(mb_full, uniquify = true) == ["50.0%", "60.0%", "70.0%", "80.0%", "90.0%"]
    @test which_density_bands(mb_empty) == String[]

    DSGE.prepare_meansbands_table_timeseries(mb_full, :obs_gdp)

    DSGE.parse_transform(Symbol("DSGE.loggrowthtopct_annualized"))([1.0]) == loggrowthtopct_annualized([1.0])

    # Ideally would also test prepare_means_table_shockdec...
end

@testset "Test meansbands_to_matrix works" begin
    meansbands_to_matrix(mb_full)
end


text:
"""
    couple_models(aero::Wagner, stru::TypicalSection, flap::SimpleFlap)

Create an aerostructural model using an unsteady aerodynamic model based on
Wagner's function, a two-degree of freedom typical section model, and a linear
steady-state control surface model.  This model introduces the freestream
velocity ``U_\\infty``, air density ``\\rho_\\infty``, and flap deflection
``\\delta`` as additional parameters.
"""
couple_models(aero::Wagner, stru::TypicalSection, flap::SimpleFlap) = (aero, stru, flap)

# --- Traits --- #

number_of_additional_parameters(::Type{<:Wagner}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = 3

coupling_inplaceness(::Type{<:Wagner}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = OutOfPlace()

coupling_rate_jacobian_type(::Type{<:Wagner}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Linear()
coupling_state_jacobian_type(::Type{<:Wagner}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_parameter_jacobian_type(::Type{<:Wagner}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Nonlinear()
coupling_time_gradient_type(::Type{<:Wagner}, ::Type{TypicalSection}, ::Type{SimpleFlap}) = Zeros()

# --- Methods --- #

function get_coupling_inputs(aero::Wagner, stru::TypicalSection, flap::SimpleFlap, dx, x, p, t)
    # extract rate variables
    dλ1, dλ2, dh, dθ, dhdot, dθdot = dx
    # extract state variables
    λ1, λ2, h, θ, hdot, θdot = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # extract model constants
    C1 = aero.C1
    C2 = aero.C2
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)
    # calculate aerodynamic loads (except contribution from state rates)
    Na, Aa, Ma = wagner_loads(a, b, ρ, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # loads per unit span
    L = Na + Nf
    M = Ma + Mf
    # return inputs
    return SVector(u, v, ω, L, M)
end

# --- Performance Overloads --- #

function get_coupling_rate_jacobian(aero::Wagner, stru::TypicalSection,
    flap::SimpleFlap, dx, x, p, t)
    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # calculate loads
    N_hddot, A_hddot, M_hddot = wagner_loads_vdot(a, b, ρ)
    N_θddot, A_θddot, M_θddot = wagner_loads_ωdot(a, b, ρ)
    # construct submatrices
    Mda = @SMatrix [0 0; 0 0; 0 0]
    Mds = @SMatrix [0 0 0 0; 0 0 0 0; 0 0 0 0]
    Mra = @SMatrix [0 0; 0 0]
    Mrs = @SMatrix [0 0 N_hddot N_θddot; 0 0 M_hddot M_θddot]
    # assemble mass matrix
    return [Mda Mds; Mra Mrs]
end

function get_coupling_state_jacobian(aero::Wagner, stru::TypicalSection,
    flap::SimpleFlap, dx, x, p, t) where {N,TF,SV,SA}
    # extract state variables
    λ1, λ2, h, θ, hdot, θdot = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p
    # extract model constants
    C1 = aero.C1
    C2 = aero.C2
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    u_θ, v_θ, ω_θ = section_velocities_θ(U)
    u_hdot, v_hdot, ω_hdot = section_velocities_hdot()
    u_θdot, v_θdot, ω_θdot = section_velocities_θdot()
    # calculate loads
    r_λ = wagner_loads_λ(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)
    N_λ, A_λ, M_λ = r_λ[1,:], r_λ[2,:], r_λ[3,:]
    N_θ, M_θ = wagner_loads_θ(a, b, ρ, a0, C1, C2, U)
    N_hdot, A_hdot, M_hdot = wagner_loads_v(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)
    N_θdot, A_θdot, M_θdot = wagner_loads_ω(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)
    # compute jacobian sub-matrices
    Jda = @SMatrix [0 0; 0 0; 0 0]
    Jds = @SMatrix [0 0 0 0; 0 v_θ v_hdot 0; 0 0 0 ω_θdot]
    Jra = vcat(N_λ', M_λ')
    Jrs = @SMatrix [0 N_θ N_hdot N_θdot; 0 M_θ M_hdot M_θdot]
    # return jacobian
    return [Jda Jds; Jra Jrs]
end

# TODO: Parameter jacobian

# --- Convenience Methods --- #

function set_additional_parameters!(padd, aero::Wagner, stru::TypicalSection,
    flap::SimpleFlap; U, rho, delta)

    padd[1] = U
    padd[2] = rho
    padd[3] = delta

    return padd
end

function separate_additional_parameters(aero::Wagner, stru::TypicalSection,
    flap::SimpleFlap, padd)

    return (U = padd[1], rho = padd[2], delta = padd[3])
end

# --- Plotting --- #

@recipe function f(aero::Wagner, stru::TypicalSection, flap::SimpleFlap, dx, x, y, p, t)

    framestyle --> :origin
    grid --> false
    xlims --> (-1.0, 1.0)
    ylims --> (-1.5, 1.5)
    label --> @sprintf("t = %6.3f", t)

    # extract state variables
    λ1, λ2, h, θ, hdot, θdot = x

    # extract parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, cnδ, caδ, cmδ, U, ρ, δ = p

    xplot = [-(0.5 + a*b)*cos(θ),    (0.5 - a*b)*cos(θ)]
    yplot = [ (0.5 + a*b)*sin(θ)-h, -(0.5 - a*b)*sin(θ)-h]

    return xplot, yplot
end


text:
using Documenter, StingerGraphs

makedocs()

deploydocs(
    julia = "nightly",
    repo = "github.com/stingergraph/StingerGraphs.jl.git"
)