text:
# load code
include("channel.jl")
include("vertex.jl")

abstract type Action end

# get interpolated / extrapolated self energy for general action
function get_Σ(
    w :: Float64,
    m :: Mesh,
    a :: Action
    ) :: Float64

    # init value
    val = 0.0

    # check if in bounds, otherwise extrapolate as 1 / w
    if abs(w) <= m.σ[end]
        p   = get_param(abs(w), m.σ)
        val = sign(w) * (p.lower_weight * a.Σ[p.lower_index] + p.upper_weight * a.Σ[p.upper_index])
    else
        val = m.σ[end] * a.Σ[end] / w
    end

    return val
end

# get interpolated vertex component for general action
function get_Γ_comp(
    comp        :: Int64,
    site        :: Int64,
    bs          :: Buffer,
    bt          :: Buffer,
    bu          :: Buffer,
    r           :: Reduced_lattice,
    a           :: Action,
    apply_flags :: Function
    ;
    ch_s        :: Bool = true,
    ch_t        :: Bool = true,
    ch_u        :: Bool = true
    )           :: Float64

    # init with bare value
    val = a.Γ[comp].bare[site]

    # add s channel
    if ch_s
        # check for site exchange
        site_s = site

        if bs.exchange_flag
            site_s = r.exchange[site_s]
        end

        # apply other flags
        comp_s        = comp
        sgn_s, comp_s = apply_flags(bs, comp_s)

        # check for mapping to u channel and interpolate
        if bs.map_flag
            val += sgn_s * get_vertex(site_s, bs, a.Γ[comp_s], 3)
        else
            val += sgn_s * get_vertex(site_s, bs, a.Γ[comp_s], 1)
        end
    end

    # add t channel
    if ch_t
        # check for site exchange
        site_t = site

        if bt.exchange_flag
            site_t = r.exchange[site_t]
        end

        # apply other flags
        comp_t        = comp
        sgn_t, comp_t = apply_flags(bt, comp_t)

        # interpolate
        val += sgn_t * get_vertex(site_t, bt, a.Γ[comp_t], 2)
    end

    # add u channel
    if ch_u
        # check for site exchange
        site_u = site

        if bu.exchange_flag
            site_u = r.exchange[site_u]
        end

        # apply other flags
        comp_u        = comp
        sgn_u, comp_u = apply_flags(bu, comp_u)

        # check for mapping to s channel and interpolate
        if bu.map_flag
            val += sgn_u * get_vertex(site_u, bu, a.Γ[comp_u], 1)
        else
            val += sgn_u * get_vertex(site_u, bu, a.Γ[comp_u], 3)
        end
    end

    return val
end

# get interpolated vertex component on all lattice sites
function get_Γ_comp_avx!(
    comp        :: Int64,
    r           :: Reduced_lattice,
    bs          :: Buffer,
    bt          :: Buffer,
    bu          :: Buffer,
    a           :: Action,
    apply_flags :: Function,
    temp        :: SubArray{Float64, 1, Array{Float64, 3}}
    ;
    ch_s        :: Bool = true,
    ch_t        :: Bool = true,
    ch_u        :: Bool = true
    )           :: Nothing

    # init with bare value
    @turbo temp .= a.Γ[comp].bare

    # add s channel
    if ch_s
        # apply flags
        comp_s        = comp
        sgn_s, comp_s = apply_flags(bs, comp_s)

        # check for mapping to u channel and interpolate
        if bs.map_flag
            get_vertex_avx!(r, bs, a.Γ[comp_s], 3, temp, bs.exchange_flag, sgn_s)
        else
            get_vertex_avx!(r, bs, a.Γ[comp_s], 1, temp, bs.exchange_flag, sgn_s)
        end
    end

    # add t channel
    if ch_t
        # apply flags
        comp_t        = comp
        sgn_t, comp_t = apply_flags(bt, comp_t)

        # interpolate
        get_vertex_avx!(r, bt, a.Γ[comp_t], 2, temp, bt.exchange_flag, sgn_t)
    end

    # add u channel
    if ch_u
        # apply flags
        comp_u        = comp
        sgn_u, comp_u = apply_flags(bu, comp_u)

        # check for mapping to s channel and interpolate
        if bu.map_flag
            get_vertex_avx!(r, bu, a.Γ[comp_u], 1, temp, bu.exchange_flag, sgn_u)
        else
            get_vertex_avx!(r, bu, a.Γ[comp_u], 3, temp, bu.exchange_flag, sgn_u)
        end
    end

    return nothing
end

# load saving and reading for channels and vertices
include("disk.jl")

# load specialized code for different symmetries
include("action_lib/action_su2.jl")  ; include("checkpoint_lib/checkpoint_su2.jl")
include("action_lib/action_u1_dm.jl"); include("checkpoint_lib/checkpoint_u1_dm.jl")

# interface function to replace action with another action (except for bare)
function replace_with!(
    a1 :: Action,
    a2 :: Action
    )  :: Nothing

    # replace self energy
    @turbo a1.Σ .= a2.Σ

    # replace vertices
    for i in eachindex(a1.Γ)
        replace_with!(a1.Γ[i], a2.Γ[i])
    end

    return nothing
end

# interface function to replace action with another action only on the vertex level (except for bare)
function replace_with_Γ!(
    a1 :: Action,
    a2 :: Action
    )  :: Nothing

    # replace vertices
    for i in eachindex(a1.Γ)
        replace_with!(a1.Γ[i], a2.Γ[i])
    end

    return nothing
end

# interface function to multiply action with factor (except for bare)
function mult_with!(
    a   :: Action,
    fac :: Float64
    )   :: Nothing

    # multiply self energy
    @turbo a.Σ .*= fac

    # multiply vertices
    for i in eachindex(a.Γ)
        mult_with!(a.Γ[i], fac)
    end

    return nothing
end

# interface function to multiply action with factor only on the vertex level (except for bare)
function mult_with_Γ!(
    a   :: Action,
    fac :: Float64
    )   :: Nothing

    # multiply vertices
    for i in eachindex(a.Γ)
        mult_with!(a.Γ[i], fac)
    end

    return nothing
end

# interface function to reset an action to zero (except for bare)
function reset!(
    a :: Action
    ) :: Nothing

    mult_with!(a, 0.0)

    return nothing
end

# interface function to reset an action to zero only on the vertex level (except for bare)
function reset_Γ!(
    a :: Action
    ) :: Nothing

    mult_with_Γ!(a, 0.0)

    return nothing
end

# interface function to multiply action with some factor and add to other action (except for bare)
function mult_with_add_to!(
    a2  :: Action,
    fac :: Float64,
    a1  :: Action
    )   :: Nothing

    # multiply add for the self energy
    @turbo a1.Σ .+= fac .* a2.Σ

    # multiply add for the vertices
    for i in eachindex(a1.Γ)
        mult_with_add_to!(a2.Γ[i], fac, a1.Γ[i])
    end

    return nothing
end

# interface function to multiply action with some factor and add to other action only on the vertex level (except for bare)
function mult_with_add_to_Γ!(
    a2  :: Action,
    fac :: Float64,
    a1  :: Action
    )   :: Nothing

    # multiply add for the vertices
    for i in eachindex(a1.Γ)
        mult_with_add_to!(a2.Γ[i], fac, a1.Γ[i])
    end

    return nothing
end

# interface function to add two actions (except for bare)
function add_to!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to!(a2, 1.0, a1)

    return nothing
end

# interface function to add two actions only on the vertex level (except for bare)
function add_to_Γ!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to_Γ!(a2, 1.0, a1)

    return nothing
end

# interface function to subtract two actions (except for bare)
function subtract_from!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to!(a2, -1.0, a1)

    return nothing
end

# interface function to subtract two actions only on the vertex level (except for bare)
function subtract_from_Γ!(
    a2 :: Action,
    a1 :: Action
    )  :: Nothing

    mult_with_add_to_Γ!(a2, -1.0, a1)

    return nothing
end

"""
    get_abs_max(
        a :: Action
        ) :: Float64

Returns maximum absolute value of an action across all vertex components.
"""
function get_abs_max(
    a :: Action
    ) :: Float64

    abs_max_Γ = zeros(Float64, length(a.Γ))

    for i in eachindex(a.Γ)
        abs_max_Γ[i] = get_abs_max(a.Γ[i])
    end

    abs_max = maximum(abs_max_Γ)

    return abs_max
end

# scan cut through channel, where x is assumed to be generated by the get_mesh function with linear fraction p0
# in order to optimize the frequency mesh we use the following heuristics:
#   a) the linear extent δ should be large enough to capture sharp peaks at finite frequencies, that is, δ >= 'point where channel has decayed to 10 percent of peak value'
#   b) the linear extent δ should be large enough to capture sign changes at finite frequencies, that is, δ >= zero crossing
#   c) the linear extent δ should be small enough to capture sharp peaks at the origin, that is, Δ should fulfill p1 < Δ < p2,
#      where Δ is the relative deviation between the values at the origin and the first finite frequency
#   d) the linear extent δ should be bounded to avoid overambitious shrinking or broadening of the mesh, that is, num_lin * p3 <= δ <= p4,
#      where num_lin is the number of linear frequencies
function scan(
    Λ  :: Float64,
    x  :: Vector{Float64},
    y  :: Vector{Float64},
    p0 :: Float64,
    p1 :: Float64,
    p2 :: Float64,
    p3 :: Float64,
    p4 :: Float64
    )  :: Float64

    # determine current mesh layout
    num_lin = ceil(Int64, p0 * (length(x) - 1))

    # a) check for peaks at finite frequencies and set δa as 'point where channel has decayed to 10 percent of peak value'
    δa  = Inf
    idx = argmax(abs.(y))

    if idx > 1
        for i in idx + 1 : length(x)
            if abs(y[i] / y[idx]) <= 0.1
                δa = x[i]
                break
            end
        end
    end

    # b) check for sign change at finite frequencies and set δb accordingly
    δb = Inf

    for i in 2 : length(y)
        if sign(y[i]) != sign(y[1])
            δb = x[i]
            break
        end
    end

    # c) check if the value at the origin is finite. If so, set δc such that p1 < Δ < p2
    δc = num_lin * x[2]

    if abs(y[1]) / abs(y[idx]) > 1e-3
        # compute Δ
        Δ = abs(y[2] - y[1]) / max(abs(y[2]), abs(y[1]))

        # adjust δc if Δ is out of required bounds
        while (p1 < Δ < p2) == false
            # if Δ is too large decrease δc by two percent
            if Δ >= p2
                δc *= 0.98
            # if Δ is too small increase δc by two percent
            elseif Δ <= p1
                δc *= 1.02
            end

            # ensure δc is within required bounds
            if δc < num_lin * p3
                break
            elseif δc > p4
                break
            end

            # generate new reference data
            xp = get_mesh(δc, x[end], length(x) - 1, p0)
            yp = similar(y)

            for i in eachindex(yp)
                p     = get_param(xp[i], x)
                yp[i] = p.lower_weight * y[p.lower_index] + p.upper_weight * y[p.upper_index]
            end

            # recompute Δ
            Δ = abs(yp[2] - yp[1]) / max(abs(yp[2]), abs(yp[1]))
        end
    end

    # d) merge results from heuristics a), b) and c) and perform sanity checks
    δ = min(max(num_lin * p3, min(δa, δb, δc)), p4)

    return δ
end

# auxiliary function to scan a single channel
function scan_channel(
    Λ   :: Float64,
    p_Ω :: NTuple{5, Float64},
    p_ν :: NTuple{5, Float64},
    Ω   :: Vector{Float64},
    ν   :: Vector{Float64},
    ch  :: Channel
    )   :: NTuple{2, Float64}

    # deref data
    q3 = ch.q3

    # determine position of the maximum
    idxs = argmax(abs.(q3))

    # get cuts through the maximum
    q3_Ω   = q3[idxs[1], :, idxs[3], idxs[4]]
    q3_ν_1 = Float64[q3[idxs[1], idxs[2],      x,       x] - q3[idxs[1], idxs[2],     end,     end] for x in eachindex(ν)]
    q3_ν_2 = Float64[q3[idxs[1], idxs[2],      x, idxs[4]] - q3[idxs[1], idxs[2],     end, idxs[4]] for x in eachindex(ν)]
    q3_ν_3 = Float64[q3[idxs[1], idxs[2], idxs[3],      x] - q3[idxs[1], idxs[2], idxs[3],     end] for x in eachindex(ν)]

    # scan bosonic cut
    Ω_lin = scan(Λ, Ω, q3_Ω, p_Ω[1], p_Ω[2], p_Ω[3], p_Ω[4] * Λ, p_Ω[5] * Λ)

    # scan fermionic cuts
    ν_lin_1 = scan(Λ, ν, q3_ν_1, p_ν[1], p_ν[2], p_ν[3], p_ν[4] * Λ, p_ν[5] * Λ)
    ν_lin_2 = scan(Λ, ν, q3_ν_2, p_ν[1], p_ν[2], p_ν[3], p_ν[4] * Λ, p_ν[5] * Λ)
    ν_lin_3 = scan(Λ, ν, q3_ν_3, p_ν[1], p_ν[2], p_ν[3], p_ν[4] * Λ, p_ν[5] * Λ)
    ν_lin   = min(ν_lin_1, ν_lin_2, ν_lin_3)

    return Ω_lin, ν_lin
end

# resample an action to new meshes via scanning and trilinear interpolation
function resample_from_to(
    Λ      :: Float64,
    p_σ    :: NTuple{2, Float64},
    p_Ωs   :: NTuple{5, Float64},
    p_νs   :: NTuple{5, Float64},
    p_Ωt   :: NTuple{5, Float64},
    p_νt   :: NTuple{5, Float64},
    p_χ    :: NTuple{5, Float64},
    lins   :: NTuple{5, Float64},
    bounds :: NTuple{5, Float64},
    m_old  :: Mesh,
    a_old  :: Action,
    a_new  :: Action,
    χ      :: Vector{Matrix{Float64}}
    )      :: Mesh

    # scale linear bounds
    σ_lin  = lins[2] * Λ
    Ωs_lin = lins[3] * Λ
    νs_lin = lins[4] * Λ
    Ωt_lin = lins[3] * Λ
    νt_lin = lins[4] * Λ
    χ_lin  = lins[5] * Λ

    # adjust meshes via scanning once required scale is reached
    if Λ < lins[1]
        # scan self energy
        σ_lin = p_σ[2] * m_old.σ[argmax(abs.(a_old.Σ))]

        # scan the channels
        Ωs_lins, νs_lins = zeros(Float64, length(a_old.Γ)), zeros(length(a_old.Γ))
        Ωt_lins, νt_lins = zeros(Float64, length(a_old.Γ)), zeros(length(a_old.Γ))

        for i in eachindex(a_old.Γ)
            Ωs_lins[i], νs_lins[i] = scan_channel(Λ, p_Ωs, p_νs, m_old.Ωs, m_old.νs, a_old.Γ[i].ch_s)
            Ωt_lins[i], νt_lins[i] = scan_channel(Λ, p_Ωt, p_νt, m_old.Ωt, m_old.νt, a_old.Γ[i].ch_t)
        end

        Ωs_lin, νs_lin = minimum(Ωs_lins), minimum(νs_lins)
        Ωt_lin, νt_lin = minimum(Ωt_lins), minimum(νt_lins)

        # scan the correlations
        χ_lins = zeros(Float64, length(χ))

        for i in eachindex(χ_lins)
            χ_lins[i] = scan(Λ, m_old.χ, χ[i][argmax(abs.(χ[i]))[1], :], p_χ[1], p_χ[2], p_χ[3], p_χ[4] * Λ, p_χ[5] * Λ)
        end

        χ_lin = minimum(χ_lins)
    end

    # build new frequency meshes according to scanning results
    σ     = get_mesh( σ_lin, bounds[2] * max(Λ, bounds[1]), m_old.num_σ - 1, p_σ[1])
    Ωs    = get_mesh(Ωs_lin, bounds[3] * max(Λ, bounds[1]), m_old.num_Ω - 1, p_Ωs[1])
    νs    = get_mesh(νs_lin, bounds[4] * max(Λ, bounds[1]), m_old.num_ν - 1, p_νs[1])
    Ωt    = get_mesh(Ωt_lin, bounds[3] * max(Λ, bounds[1]), m_old.num_Ω - 1, p_Ωt[1])
    νt    = get_mesh(νt_lin, bounds[4] * max(Λ, bounds[1]), m_old.num_ν - 1, p_νt[1])
    χ     = get_mesh( χ_lin, bounds[5] * max(Λ, bounds[1]), m_old.num_χ - 1, p_χ[1])
    m_new = Mesh(m_old.num_σ, m_old.num_Ω, m_old.num_ν, m_old.num_χ, σ, Ωs, νs, Ωt, νt, χ)

    # resample self energy
    for w in eachindex(m_new.σ)
        a_new.Σ[w] = get_Σ(m_new.σ[w], m_old, a_old)
    end

    # resample vertices
    for i in eachindex(a_new.Γ)
        resample_from_to!(m_old, a_old.Γ[i], m_new, a_new.Γ[i])
    end

    return m_new
end

# generate action dummy
function get_action_empty(
    symmetry :: String,
    r        :: Reduced_lattice,
    m        :: Mesh
    ;
    S        :: Float64 = 0.5
    )        :: Action

    if symmetry == "su2"
        return get_action_su2_empty(S, r, m)
    elseif symmetry == "u1-dm"
        return get_action_u1_dm_empty(r, m)
    end
end

"""
    read_checkpoint(
        file     :: HDF5.File,
        Λ        :: Float64
        )        :: Tuple{Float64, Float64, Mesh, Action}

Read checkpoint of FRG calculation from HDF5 file.
Returns cutoff Λ, ODE stepwidth dΛ, frequency meshes (wrapped in Mesh struct) and vertices (wrapped in Action struct).
"""
function read_checkpoint(
    file     :: HDF5.File,
    Λ        :: Float64
    )        :: Tuple{Float64, Float64, Mesh, Action}

    # read symmetry group from file
    symmetry = read(file, "symmetry")

    if symmetry == "su2"
        return read_checkpoint_su2(file, Λ)
    elseif symmetry == "u1-dm"
        return read_checkpoint_u1_dm(file, Λ)
    end
end

# load tests and timers
include("test.jl")
include("timers.jl")


text:
"""
    SIS

A type representing a Susceptible-Infected-Susceptible model with infection rate `β` and recovery rate `γ`.
"""
struct SIS <: AbstractCompartimentalModel
    β::Real
    γ::Real
end

SI(sis::SIS) = SI(sis.β)

num_states(::SIS) = 2

# Metapopulation

function rate(i, state, mp::Metapopulation{SIS})
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    h = mp.h
    β*state[i,1]*state[i,2] + γ*state[i,2] + sum(D .* state[i,:]) * (outdegree(h,i) > 0)
end

function update_state_and_rates!(state, a, k, mp::Metapopulation{SIS})
    N = nv(mp.h)
    D = mp.D
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    h = mp.h
    p = vcat(β*state[k,1]*state[k,2], γ*state[k,2], D.*state[k,:]*(outdegree(h,k)>0))
    j = sample(ProbabilityWeights(p, a[k]))
    if j == 1
        state[k,1] -= 1
        state[k,2] += 1
        a[k] = rate(k, state, mp)
    elseif j == 2
        state[k,1] += 1
        state[k,2] -= 1
        a[k] = rate(k, state, mp)
    else
        μ = k
        ν = rand(outneighbors(h, μ))
        state[μ,j-2] -= 1 # NB. j-2 is the state of the individual that migrates
        state[ν,j-2] += 1 # NB. because the first 2 elements of p are the reactions
        a[μ] = rate(μ, state, mp)
        a[ν] = rate(ν, state, mp)
    end
end

function meanfield_fun(mp::Metapopulation{SIS})
    L = normalized_laplacian(mp.h)
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    f! = function(dx, x, p, t)
        tmp = (β.*x[:,1] .- γ).*x[:,2]
        dx[:,1] .= .-tmp .- D[1].*(L'*x[:,1])
        dx[:,2] .= tmp .- D[2].*(L'x[:,2])
    end
    return f!
end

# Contact Process

function rate(k, state, cp::ContactProcess{SIS})
    if state[k] == 1
        l = length(filter(i->state[i]==2, inneighbors(cp.g, k)))
        return cp.dynamics.β * l
    else
        return cp.dynamics.γ
    end
end

function update_state_and_rates!(state, a, k, cp::ContactProcess{SIS})
    g = cp.g
    β = cp.dynamics.β
    if state[k] == 1
        state[k] = 2
        for i in Iterators.filter(j->state[j]==1, outneighbors(g, k))
            a[i] += β
        end
    else
        state[k] = 1
        for i in Iterators.filter(j->state[j]==1, outneighbors(g, k))
            a[i] -= β
        end
    end
    a[k] = rate(k, state, cp)
end

function update_output!(output, state, n, k, cp::ContactProcess{SIS})
    if state[k] == 2 # newly infected node
        output[n] = output[n-1] + [-1, 1]
    else
        output[n] = output[n-1] + [1, -1]
    end
end

function init_state_mf(cp::ContactProcess{SIS}, x0)
    return float(x0 .- 1)
end

function meanfield_fun(cp::ContactProcess{SIS})
    A = adjacency_matrix(cp.g)
    β = cp.dynamics.β
    γ = cp.dynamics.γ
    f! = function(dx, x, p, t)
        dx .= β.*(1.0 .- x).*(A*x) .- γ.*x
    end
    return f!
end

function finalize_meanfield(cp::ContactProcess{SIS}, sol)
    N = nv(cp.g)
    n = length(sol.u)
    output = Array{Float64,2}(undef, n, N)
    for i in 1:n
        output[i,:] .= sol.u[i]
    end
    return sol.t, output
end

# Metaplex

function rate(k, state, mpx::Metaplex{SIS})
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    od = outdegree(mpx.h, x_μ[k]) > 0
    if x_i[k] == 1
        l = length(filter(i->x_i[i]==2 && x_μ[i]==x_μ[k], inneighbors(mpx.g, k)))
        return β*l + D[1]*od
    else
        return γ + D[2]*od
    end
end

function update_state_and_rates!(state, a, k, mpx::Metaplex{SIS})
    N = nv(mpx.g)
    M = nv(mpx.h)
    g = mpx.g
    h = mpx.h
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    popcounts = state[:popcounts]
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    μ = x_μ[k]
    if x_i[k] == 1
        if rand()*a[k] < D[1] # susceptible node migrates
            ν = rand(outneighbors(h, μ))
            x_μ[k] = ν
            popcounts[μ,1] -= 1
            popcounts[ν,1] += 1
            a[k] = rate(k, state, mpx)
        else # susceptible node becomes infected
            x_i[k] = 2
            popcounts[μ,1] -= 1
            popcounts[μ,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g,k))
                a[i] += β
            end
        end
    else
        if rand()*a[k] < D[2] # infected node migrates
            ν = rand(outneighbors(h,μ))
            x_μ[k] = ν
            popcounts[μ,2] -= 1
            popcounts[ν,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1, outneighbors(g,k))
                if x_μ[i] == μ
                    a[i] -= β
                elseif x_μ[i] == ν
                    a[i] += β
                end
            end
        else # infected node recovers
            x_i[k] = 1
            popcounts[μ,2] -= 1
            popcounts[μ,1] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g,k))
                a[i] -= β
            end
        end
    end
end

function meanfield_fun(mpx::Metaplex{SIS})
    N = nv(mpx.g)
    M = nv(mpx.h)
    A = adjacency_matrix(mpx.g)
    L = normalized_laplacian(mpx.h)
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    f! = function(dx, x, p, t)
        for i in 1:N
            dx[1,i,:] .= .-D[1]*L'*x[1,i,:]
            dx[2,i,:] .= .-D[2]*L'*x[2,i,:]
        end
        for μ in 1:M
            infection = β.*x[1,:,μ].*(A*x[2,:,μ]) .- γ.*x[2,:,μ]
            dx[1,:,μ] .-= infection
            dx[2,:,μ] .+= infection
        end
    end
    return f!
end

# Heterogeneous Metaplex

function rate(k, state, mpx::HeterogeneousMetaplex{SIS})
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    μ = x_μ[k]
    od = outdegree(mpx.h, x_μ[k]) > 0
    if x_i[k] == 1
        l = length(filter(i->x_i[i]==2 && x_μ[i]==x_μ[k], inneighbors(mpx.g[μ], k)))
        return β*l + D[1]*od
    else
        return γ + D[2]*od
    end
end

function update_state_and_rates!(state, a, k, mpx::HeterogeneousMetaplex{SIS})
    g = mpx.g
    h = mpx.h
    x_i = state[:state_i]
    x_μ = state[:state_μ]
    popcounts = state[:popcounts]
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    μ = x_μ[k]
    if x_i[k] == 1
        if rand()*a[k] < D[1] # susceptible node migrates
            ν = rand(outneighbors(h, μ))
            x_μ[k] = ν
            popcounts[μ,1] -= 1
            popcounts[ν,1] += 1
            a[k] = rate(k, state, mpx)
        else # susceptible node becomes infected
            x_i[k] = 2
            popcounts[μ,1] -= 1
            popcounts[μ,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g[μ],k))
                a[i] += β
            end
        end
    else
        if rand()*a[k] < D[2] # infected node migrates
            ν = rand(outneighbors(h,μ))
            x_μ[k] = ν
            popcounts[μ,2] -= 1
            popcounts[ν,2] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1, outneighbors(g[μ],k))
                if x_μ[i] == μ
                    a[i] -= β
                elseif x_μ[i] == ν
                    a[i] += β
                end
            end
        else # infected node recovers
            x_i[k] = 1
            popcounts[μ,2] -= 1
            popcounts[μ,1] += 1
            a[k] = rate(k, state, mpx)
            for i in Iterators.filter(j->x_i[j]==1 && x_μ[j]==μ, outneighbors(g[μ],k))
                a[i] -= β
            end
        end
    end
end

function meanfield_fun(mpx::HeterogeneousMetaplex{SIS})
    N = nv(mpx.g[1])
    M = nv(mpx.h)
    A = adjacency_matrix.(mpx.g)
    L = normalized_laplacian(mpx.h)
    β = mpx.dynamics.β
    γ = mpx.dynamics.γ
    D = mpx.D
    f! = function(dx, x, p, t)
        for i in 1:N
            dx[1,i,:] .= .-D[1]*L'*x[1,i,:]
            dx[2,i,:] .= .-D[2]*L'*x[2,i,:]
        end
        for μ in 1:M
            infection = β.*x[1,:,μ].*(A[μ]*x[2,:,μ]) .- γ.*x[2,:,μ]
            dx[1,:,μ] .-= infection
            dx[2,:,μ] .+= infection
        end
    end
    return f!
end

# HeterogeneousMetapopulation

function rate(i, state, mp::HeterogeneousMetapopulation{SIS})
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    h = mp.h
    ks = mp.ks
    N = mp.N
    βᵢ = β * ks[i] / N
    return βᵢ*state[i,1]*state[i,2] + γ*state[i,2] + sum(D .* state[i,:]) * (outdegree(h,i) > 0)
end

function update_state_and_rates!(state, a, k, mp::HeterogeneousMetapopulation{SIS})
    M = nv(mp.h)
    N = mp.N
    ks = mp.ks
    D = mp.D
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    h = mp.h
    βₖ = β * ks[k] / N
    p = vcat(βₖ*state[k,1]*state[k,2], γ*state[k,2], D.*state[k,:]*(outdegree(h,k)>0))
    j = sample(ProbabilityWeights(p, a[k]))
    if j == 1
        state[k,1] -= 1
        state[k,2] += 1
        a[k] = rate(k, state, mp)
    elseif j == 2
        state[k,1] += 1
        state[k,2] -= 1
        a[k] = rate(k, state, mp)
    else
        μ = k
        ν = rand(outneighbors(h, μ))
        state[μ,j-2] -= 1 # NB. j-2 is the state of the individual that migrates
        state[ν,j-2] += 1 # NB. because the first 2 elements of p are the reactions
        a[μ] = rate(μ, state, mp)
        a[ν] = rate(ν, state, mp)
    end
end

function meanfield_fun(mp::HeterogeneousMetapopulation{SIS})
    #L = normalized_laplacian(mp.h)
    L = laplacian_matrix(mp.h)
    ks = mp.ks
    N = mp.N
    β = mp.dynamics.β
    γ = mp.dynamics.γ
    D = mp.D
    βs = β * ks / N
    f! = function(dx, x, p, t)
        tmp = (βs .* x[:,1] .- γ).*x[:,2]
        dx[:,1] .= .-tmp .- D[1].*(L'*x[:,1])
        dx[:,2] .= tmp .- D[2].*(L'x[:,2])
    end
    return f!
end

ext: jl


text:
import ..UncertainDatasets:
	AbstractUncertainValueDataset,
	UncertainValueDataset,
	UVAL_COLLECTION_TYPES

import ..SamplingConstraints:
	NoConstraint,
	TruncateLowerQuantile,
	TruncateUpperQuantile,
	TruncateQuantiles

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint) -> Vector{T} where T

Resample `x` (a collection of uncertain values) once, applying the provided sampling `constraint`.

Returns a `length(x)`-element vector. The `i`-th element of this vector is generated by
truncating the `i`-th uncertain value by the sampling `constraint`, then drawing a single random
number from the truncated value.

See also [`UVAL_COLLECTION_TYPES`](@ref).

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate each distribution at +- 0.5 standard deviations, then resample.
resample(uvals, TruncateStd(0.5))
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint)
	[resample(uv.values[i], constraint) for i in 1:length(uv)]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint}) -> Vector{T} where T

Resample `x` (a collection of uncertain values) once, applying the provided sampling `constraint`s.
The number of constraints must match the number of elements in `x`.

Returns a `length(x)`-element vector. The `i`-th element of this vector is generated by
truncating the `i`-th uncertain value by the `i`-th sampling `constraint`, then drawing a single random
number from the truncated value.

See also [`UVAL_COLLECTION_TYPES`](@ref).

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate each distribution at +- 0.5 standard deviations, then resample.
resample(uvals, TruncateStd(0.5))
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint})
    [resample(uv.values[i], constraint[i]) for i in 1:length(uv)]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint, n::Int) -> Vector{Vector{T}} where T

Resample `x` (a collection of uncertain values) `n` times, applying the provided sampling `constraint`.

Returns an `n`-element vector of `length(x)`-element vectors. Each of these vectors is an independent
draw from `x`. The `i`-th element of each draw is generated by truncating the `i`-th uncertain value by
the sampling `constraint`, then drawing a single random number from the truncated value.

See also [`UVAL_COLLECTION_TYPES`](@ref).

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate the first 50 elements at the 90th percentile range, and the
# last 50 elements at the 40th percentile range.
constraints = [i <= 50 ? TruncateQuantiles(0.05, 0.95) : TruncateQuantiles(0.3, 0.7) for i = 1:100]

# Truncate the distributions, then draw ten independent realisations of the collection subject
# to the provided constraints.
resample(uvals, constraints, 10)
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint, n::Int)
	[[resample(uv.values[i], constraint) for i in 1:length(uv)] for k = 1:n]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint}, n::Int) -> Vector{Vector{T}} where T

Resample `x` (a collection of uncertain values) `n` times, applying the provided sampling `constraint`s.

Returns an `n`-element vector of `length(x)`-element vectors. Each of these vectors is an independent
draw from `x`. The `i`-th element of each draw is generated by truncating the `i`-th uncertain value by
the `i`-th sampling `constraint`, then drawing a single random number from the truncated value.

See also [`UVAL_COLLECTION_TYPES`](@ref).

## Example

```julia
# Generate some uncertain values where the `i`-th value is given by a normal
# distribution with mean `i` and a standard deviation drawn from a uniform
# distribution on `[0, 1]`.
uvals = [UncertainValue(Normal(i, rand())) for i = 1:100]

# Truncate the first 50 elements at `± 0.5` standard deviations, and the
# last 50 elements at `± 1.2` standar deviations.
constraints = [i <= 50 ? TruncateStd(0.5) : TruncateStd(1.2) for i = 1:100]

# Apply the constraints element-wise, then draw ten independent realisations
# of the collection subject to those constraints.
resample(uvals, constraints, 10)
```
"""
function resample(uv::UVAL_COLLECTION_TYPES, constraint::Vector{<:SamplingConstraint}, n::Int)
	[[resample(uv.values[i], constraint[i]) for i in 1:length(uv)] for k = 1:n]
end

"""
	resample(x::UVAL_COLLECTION_TYPES, constraint::SamplingConstraint, n::Int) -> Vector{Vector{T}} where T

Resample `x` (a collection of uncertain values) once by drawing a single random number from
each of the uncertain values in `x`.

See also [`UVAL_COLLECTION_TYPES`](@ref).

## Example

```julia
# Generate some uncertain values represented by gamma distributions
uvals = [UncertainValue(Gamma(i, rand())) for i = 1:100]

# Resample the collection once
resample(uvals)
```
"""
function resample(uvd::UVAL_COLLECTION_TYPES)
	L = length(uvd)
	[resample(uvd.values[i]) for i in 1:L]
end

"""
	resample(uvd::UVAL_COLLECTION_TYPES, n::Int) -> Vector{Vector{T}}

Draw `n` realisations of an uncertain value dataset according to the distributions
of the uncertain values comprising it.

See also [`UVAL_COLLECTION_TYPES`](@ref).

## Example

```julia
# Generate some uncertain values represented by gamma distributions
uvals = [UncertainValue(Gamma(i, rand())) for i = 1:100]

# Resample the collection once
resample(uvals)
```
"""
function resample(uvd::UVAL_COLLECTION_TYPES, n::Int)
	L = length(uvd)
	[[resample(uvd.values[i]) for i in 1:L] for k in 1:n]
end


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
include("faberschauder.jl")
using ZigZagBoomerang, SparseArrays, LinearAlgebra, Revise
# Drift
α = 1.5
β = 0.1
b(x) = -β*x + α*sin(2pi*x)
# First derivative
b′(x) = -β + 2*α*2pi*cos(2pi*x)
# Second derivative
b″(x) = -2*α*(2pi)^2*sin(2pi*x)
# length diffusion bridge
T = 2.0

L = 7
n = (2 << L) + 1
ξ0 = 0randn(n)
u, v = 1.5, -0.5  # initial and fianl point
ξ0[1] = u/sqrt(T)
ξ0[end] = v/sqrt(T)
c = ones(n)
c[end] = c[1] = 0.0
θ0 = rand((-1.0, 1.0), n)
θ0[end] = θ0[1] =  0.0 # fix final point
T′ = 5000.0 # final clock of the pdmp

Γ = sparse(1.0I, n, n)
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ!, 0.0, ξ0, θ0, T, 10.0, Boomerang(Γ, ξ0*0, 0.1; ρ=0.9), 1, L, adapt=false);
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ, 0.0, ξ0, rand((-1.0, 1.0), n), T, 40.0*ones(n), ZigZag(Γ, ξ0*0), 5, L, adapt=false);
trace, (t, ξ, θ), (acc, num), c = @time spdmp(∇ϕmoving, 0.0, ξ0,
    θ0, T′, c, ZigZag(Γ, ξ0*0), SelfMoving(), L, T, adapt=true);
#trace, (t, ξ, θ), (acc, num) = @time pdmp(∇ϕ, 0.0, ξ0, rand((-1.0,1.0), n), T, 100.0*ones(n), FactBoomerang(Γ, ξ0*0, 0.1), 5, L, adapt=false);
ts, ξs = splitpairs(discretize(trace, T′/n))
S = T*(0:n)/(n+1)

### change to false is you want to Plot the results
if true
    error("Stop before plotting the results...")
end

using CairoMakie
p1 = lines(S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
for ξ in ξs[1:5:end]
    lines!(S, [dotψ(ξ, s, L, T) for s in S], linewidth=0.3)
end
display(p1)
save("./diffbridges.png", p1)


text:
using Documenter, MicrobiomePlots

makedocs(
    format = :html,
    sitename = "MicrobiomePlots.jl",
    pages = [
        "Home" => "index.md",
        "Recipes" => "recipes.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD"
)

deploydocs(
    repo = "github.com/BioJulia/MicrobiomePlots.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
Base.:(==)(t1::CXType, t2::CXType) = clang_equalTypes(t1, t2) != 0

"""
    getCanonicalType(t::CXType) -> CXType
    getCanonicalType(t::CLType) -> CLType
Return the canonical type for a CXType.

Clang's type system explicitly models typedefs and all the ways a specific type can be
represented. The canonical type is the underlying type with all the "sugar" removed.
For example, if 'T' is a typedef for 'int', the canonical type for 'T' would be 'int'.
Wrapper for libclang's [`clang_getCanonicalType`](@ref).
"""
getCanonicalType(t::CXType) = clang_getCanonicalType(t)
getCanonicalType(t::CLType)::CLType = clang_getCanonicalType(t)

"""
    isConstQualifiedType(t::Union{CXType,CLType}) -> Bool
Determine whether a CXType has the "const" qualifier set, without looking through typedefs
that may have added "const" at a different level.
Wrapper for libclang's [`clang_isConstQualifiedType`](@ref).
"""
isConstQualifiedType(t::Union{CXType,CLType})::Bool = clang_isConstQualifiedType(t)

"""
    isVolatileQualifiedType(t::Union{CXType,CLType}) -> Bool
Determine whether a CXType has the "volatile" qualifier set, without looking through typedefs
that may have added "volatile" at a different level.
Wrapper for libclang's [`clang_isVolatileQualifiedType`](@ref).
"""
isVolatileQualifiedType(t::CXType)::Bool = clang_isVolatileQualifiedType(t)

"""
    isRestrictQualifiedType(t::Union{CXType,CLType}) -> Bool
Determine whether a CXType has the "restrict" qualifier set, without looking through typedefs
that may have added "restrict" at a different level.
Wrapper for libclang's [`clang_isRestrictQualifiedType`](@ref).
"""
isRestrictQualifiedType(t::Union{CXType,CLType})::Bool = clang_isRestrictQualifiedType(t)

"""
    getAddressSpace(t::Union{CXType,CLType})
Returns the address space of the given type.
Wrapper for libclang's [`clang_getAddressSpace`](@ref).
"""
getAddressSpace(t::Union{CXType,CLType}) = clang_getAddressSpace(t)

"""
    getTypedefName(t::Union{CXType,CLType}) -> String
Return the typedef name of the given type.
Wrapper for libclang's [`clang_getTypedefName`](@ref).
"""
function getTypedefName(t::Union{CXType,CLType})
    cxstr = clang_getTypedefName(t)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    getPointeeType(t::CXType) -> CXType
    getPointeeType(t::CLType) -> CLType
Return the type of the pointee for pointer types.
Wrapper for libclang's [`clang_getPointeeType`](@ref).
"""
getPointeeType(t::CXType) = clang_getPointeeType(t)
getPointeeType(t::CLType)::CLType = clang_getPointeeType(t)

"""
    getTypeDeclaration(t::CXType) -> CXCursor
    getTypeDeclaration(t::CLType) -> CLCursor
Return the cursor for the declaration of the given type. To get the type of the cursor,
see [`type`](@ref). Wrapper for libclang's [`clang_getTypeDeclaration`](@ref).
"""
getTypeDeclaration(t::CXType) = clang_getTypeDeclaration(t)
getTypeDeclaration(t::CLType)::CLCursor = clang_getTypeDeclaration(t)

"""
    getFunctionTypeCallingConv(t::Union{CXType,CLFunctionNoProto,CLFunctionProto}) -> CXCallingConv
Return the calling convention associated with a function type.
Wrapper for libclang's [`clang_getFunctionTypeCallingConv`](@ref).
"""
getFunctionTypeCallingConv(t::Union{CXType,CLFunctionNoProto,CLFunctionProto})::CXCallingConv = clang_getFunctionTypeCallingConv(t)

"""
    getCursorResultType(t::CXType) -> CXType
    getCursorResultType(t::Union{CLFunctionNoProto,CLFunctionProto}) -> CLType
Return the return type associated with a function type.
Wrapper for libclang's [`clang_getResultType`](@ref).
"""
getCursorResultType(t::CXType) = clang_getResultType(t)
getCursorResultType(t::Union{CLFunctionNoProto,CLFunctionProto})::CLType = clang_getResultType(t)

## TODO:
# clang_getExceptionSpecificationType

"""
    getNumArguments(t::Union{CXType,CLFunctionNoProto,CLFunctionProto}) -> Int
Return the number of non-variadic parameters associated with a function type.
Wrapper for libclang's [`clang_getNumArgTypes`](@ref).
"""
getNumArguments(t::Union{CXType,CLFunctionNoProto,CLFunctionProto,CLUnexposed})::Int = clang_getNumArgTypes(t)

"""
    getArgType(t::CXType, i::Unsigned) -> CXType
    getArgType(t::Union{CLFunctionNoProto,CLFunctionProto}, i::Integer) -> CLType
Return the type of a parameter of a function type.
Wrapper for libclang's [`clang_getArgType`](@ref).
"""
getArgType(t::CXType, i::Unsigned) = clang_getArgType(t, Unsigned(i))
getArgType(t::Union{CLFunctionNoProto,CLFunctionProto,CLUnexposed}, i::Integer)::CLType = getArgType(t.type, Unsigned(i))

"""
    isVariadic(t::Union{CXType,CLType}) -> Bool
Return true if the CXType is a variadic function type.
Wrapper for libclang's [`clang_isFunctionTypeVariadic`](@ref).
"""
isVariadic(t::Union{CXType,CLType})::Bool = clang_isFunctionTypeVariadic(t)

"""
    isPODType(t::Union{CXType,CLType}) -> Bool
Return true if the CXType is a plain old data type.
Wrapper for libclang's [`clang_isPODType`](@ref).
"""
isPODType(t::CXType)::Bool = clang_isPODType(t)

"""
    getElementType(t::CXType) -> CXType
    getElementType(t::Union{CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray,CLComplex}) -> CLType
Return the element type of an array, complex, or vector type.
Wrapper for libclang's [`clang_getElementType`](@ref).
"""
getElementType(t::CXType) = clang_getElementType(t)
getElementType(t::Union{CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray,CLComplex})::CLType = clang_getElementType(t)
# old API: clang_getArrayElementType

"""
    getNumElements(t::Union{CXType,CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray}) -> Int
Return the number of elements of an array or vector type.
Wrapper for libclang's [`clang_getNumElements`](@ref).
"""
getNumElements(t::Union{CXType,CLVector,CLConstantArray,CLIncompleteArray,CLVariableArray,CLDependentSizedArray})::Int = clang_getNumElements(t)
# old API: clang_getArraySize

"""
    getNamedType(t::CXType) -> CXType
    getNamedType(t::CLElaborated) -> CLType
Return the type named by the qualified-id.
Wrapper for libclang's [`clang_Type_getNamedType`](@ref).
"""
getNamedType(t::CXType) = clang_Type_getNamedType(t)
getNamedType(t::CLElaborated)::CLType = clang_Type_getNamedType(t)

"""
    getSizeOf(t::CXType) -> Int
    getSizeOf(t::CLType) -> Int
Return the size of a type in bytes as per C++[expr.sizeof] standard,

It returns a minus number for layout errors, please convert the result to
a [`CXTypeLayoutError`](@ref) to see what the error is.
"""
getSizeOf(t::CXType)::Int = clang_Type_getSizeOf(t)
getSizeOf(t::CLType) = getSizeOf(t.type)

"""
    getAlignOf(t::CXType) -> Int
    getAlignOf(t::CLType) -> Int
Return the alignment of a type in bytes as per C++[expr.alignof] standard.

It returns a minus number for layout errors, please convert the result to
a [`CXTypeLayoutError`](@ref) to see what the error is.
"""
getAlignOf(t::CXType)::Int = clang_Type_getAlignOf(t)
getAlignOf(t::CLType) = getAlignOf(t.type)

"""
    getOffsetOf(t::CXType, s) -> Int
    getOffsetOf(t::CLType, s::AbstractString) -> Int
Return the offset of a field named S in a record of type T in bits as it
would be returned by __offsetof__ as per C++11[18.2p4].

It returns a minus number for layout errors, please convert the result to
a [`CXTypeLayoutError`](@ref) to see what the error is.
"""
getOffsetOf(t::CXType, s)::Int = clang_Type_getOffsetOf(t, s)
getOffsetOf(t::CLType, s::AbstractString) = getOffsetOf(t.type, s)

"""
    isInvalid(t::CXType) -> Bool
    isInvalid(t::CLType) -> Bool
Return true if the type is a valid type.
"""
isInvalid(t::CXType) = kind(t) != CXType_Invalid
isInvalid(t::CLType) = isInvalid(t.type)

## TODO:
# clang_Type_isTransparentTagTypedef (ObjectiveC)
# clang_Type_getClassType (C++)
# clang_Type_getNumTemplateArguments (C++)
# clang_Type_getTemplateArgumentAsType (C++)
# clang_Type_getCXXRefQualifier (C++)
#

# helper
"""
    kind(t::CXType) -> CXTypeKind
    kind(t::CLType) -> CXTypeKind
Return the kind of the given type.
"""
kind(t::CXType) = t.kind
kind(t::CLType) = kind(t.type)

"""
    spelling(t::Union{CXType,CLType}) -> String
Pretty-print the underlying type using the rules of the language of the translation unit
from which it came. If the type is invalid, an empty string is returned.
Wrapper for libclang's [`clang_getTypeSpelling`](@ref).
"""
function spelling(t::Union{CXType,CLType})
    cxstr = clang_getTypeSpelling(t)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    spelling(kind::CXTypeKind) -> String
Return the spelling of a given CXTypeKind.
"""
function spelling(kind::CXTypeKind)
    cxstr = clang_getTypeKindSpelling(kind)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

# visitor
function ct_field_visitor(cursor::CXCursor, list)::CXVisitorResult
    push!(list, cursor)
    return CXVisit_Continue
end

"""
    fields(ty::CXType) -> Vector{CXCursor}
    fields(ty::CLType) -> Vector{CLCursor}
Return a child cursor vector of the given cursor.
"""
function fields(ty::CXType)
    list = CXCursor[]
    ct_visitor_cb = @cfunction(
        ct_field_visitor, CXVisitorResult, (CXCursor, Ref{Vector{CXCursor}})
    )
    GC.@preserve list ccall(
        (:clang_Type_visitFields, LibClang.libclang),
        UInt32,
        (CXType, CXFieldVisitor, Any),
        ty,
        ct_visitor_cb,
        list,
    )
    return list
end
fields(ty::CLType)::Vector{CLCursor} = fields(ty.type)

is_elaborated(ty::CXType) = kind(ty) == CXType_Elaborated
is_elaborated(ty::CLType) = is_elaborated(ty.type)

"""
    has_elaborated_reference(ty::CLType) -> Bool
    has_elaborated_reference(ty::CXType) -> Bool
Return true if the type is an elaborated type or the type indirectly refers to an
elaborated type.
"""
function has_elaborated_reference(ty::CXType)
    if kind(ty) == CXType_Pointer
        ptreety = getPointeeType(ty)
        return has_elaborated_reference(ptreety)
    elseif kind(ty) == CXType_ConstantArray || kind(ty) == CXType_IncompleteArray
        elty = getElementType(ty)
        return has_elaborated_reference(elty)
    else
        return is_elaborated(ty)
    end
end
has_elaborated_reference(ty::CLType) = has_elaborated_reference(ty.type)

"""
    get_elaborated_cursor(ty::CLType) -> CLCursor
    get_elaborated_cursor(ty::CXType) -> CXCursor
Return the cursor of the elaborated type that is referenced by the input type.
The input type can be a pointer/array.
This function returns [`clang_getNullCursor`](@ref) if the input type is not refer to an
elaborated type.
"""
function get_elaborated_cursor(ty::CXType)
    if kind(ty) == CXType_Pointer
        ptreety = getPointeeType(ty)
        return get_elaborated_cursor(ptreety)
    elseif kind(ty) == CXType_ConstantArray || kind(ty) == CXType_IncompleteArray
        elty = getElementType(ty)
        return get_elaborated_cursor(elty)
    elseif is_elaborated(ty)
        return getTypeDeclaration(ty)
    else
        return getNullCursor()
    end
end
get_elaborated_cursor(ty::CLType)::CLCursor = get_elaborated_cursor(ty.type)

function is_function(ty::CXType)
    canonical_ty = kind(ty) == CXType_Typedef ? getCanonicalType(ty) : ty
    k = kind(canonical_ty)
    return k == CXType_FunctionProto || k == CXType_FunctionNoProto
end
is_function(ty::CLType) = is_function(ty.type)

"""
    has_function_reference(ty::CLType) -> Bool
    has_function_reference(ty::CXType) -> Bool
Return true if the type is a function or the type indirectly refers to a function.
"""
function has_function_reference(ty::CXType)
    if kind(ty) == CXType_Pointer
        ptreety = getPointeeType(ty)
        return has_function_reference(ptreety)
    elseif kind(ty) == CXType_ConstantArray || kind(ty) == CXType_IncompleteArray
        elty = getElementType(ty)
        return has_function_reference(elty)
    else
        return is_function(ty)
    end
end
has_function_reference(ty::CLType) = has_function_reference(ty.type)


text:
using Gnuplot
x = 0:0.01:2π
files = Vector{String}()
for i in 1:50
    @gp xlab = "x" ylab = "sin"  x sin.(x .+ i / 10.0) "w l lc 'black' notit"
    save(term="pngcairo size 600,400", output="advancedAnim$(i).png")
    push!(files, "advancedAnim$(i).png")
end
run(`convert -limit memory 1 -limit map 1 $files advancedAnim.gif`)
run(`rm -f $files`)

# https://github.com/gcalderone/Gnuplot.jl/issues/9

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
include("./params.jl")
include("./utils.jl")

function compute_optimal_bw{T<:Real}(C::Function, ws::AbstractArray{T,1})::T

    ################# WRITE SOMETHING #################

    return wopt

end

function get_costfunc(data)
    # Cost function

    return C
end

function sskernel(data; ws=nothing)
    # Example usage:
    # x = 0.5-0.5*log(rand(1,1e3));
    # ws = linspace(0.0001, 1.0, 100);
    # rs = sskernel(x, ws=ws); # or simply sskernel(x)
    #
    # Input argument
    # data: Spike-rate.
    # w:    Bandwidth for radial basis function
    #
    # Output argument
    # rs:   Estimated intensity function

    ################# WRITE SOMETHING #################

    return rs
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
using Millboard
using Base.Test

board = ()
@test """
+==+""" == table(board) |> string

board = ([],)
@test """
+===+===+
|   | 1 |
+===+===+
| 1 |   |
+---+---+""" == table(board) |> string

board = ((1,2,3),)
@test """
+===+=========+
|   |       1 |
+===+=========+
| 1 | (1,2,3) |
+---+---------+""" == table(board) |> string

board = ([1,2,3],)
@test """
+===+===+
|   | 1 |
+===+===+
|   | 1 |
|   | 2 |
| 1 | 3 |
+---+---+""" == table(board) |> string

board = [()]
@test """
+===+====+
|   |  1 |
+===+====+
| 1 | () |
+---+----+""" == table(board) |> string

board = [(1,2,3)]
@test """
+===+=========+
|   |       1 |
+===+=========+
| 1 | (1,2,3) |
+---+---------+""" == table(board) |> string

board = [(1,2,3) (4,5,6)]
@test """
+===+=========+=========+
|   |       1 |       2 |
+===+=========+=========+
| 1 | (1,2,3) | (4,5,6) |
+---+---------+---------+""" == table(board) |> string


text:
function expa(stoichiometric_matrix::Array{Float64,2})

    # initialize -
	(ℳ,ℛ) = size(stoichiometric_matrix)
    T = [Matrix{Float64}(I,ℛ,ℛ) transpose(stoichiometric_matrix)];

    # outer loop -
    for col_index = 1:ℳ

        # adjust the col index (from the org matrix -)
	    local_col_index = col_index + ℛ;

        # get the current RHS -
        RHS = T[:,(ℛ + 1):end] # this gets the cols belonging to the stoichiometric matrix -

        # Pick a col -
	    TMP_COL = RHS[:,col_index]; # this is grabbing a metabolite -

        # how many zeros, + and - in this col -
        idx_zero_array = findall(x->x==0.0,TMP_COL)
        idx_positive_array = findall(x->x>0.0,TMP_COL)
        idx_negative_array = findall(x->x<0.0,TMP_COL)

        # make new T -
        TNEW = T[idx_zero_array,:]

        # zero out non-zero elements -
        TMP_ARR = Array{Array{Float64,1},1}()
        for idx_positive ∈ idx_positive_array
            for idx_negative ∈ idx_negative_array

                # setup α and β -
			    α = abs(T[idx_positive,local_col_index]);
			    β = abs(T[idx_negative,local_col_index]);

                # compute a tmp row -
                TMP_ROW = α*T[idx_negative,:] .+ β*T[idx_positive,:];

                # grab the tmp row -
                push!(TMP_ARR,TMP_ROW)
            end
        end

		# ok, so we have a new set of possible pathways, check for independent pathways -
        P = transpose(hcat(TMP_ARR...))
		if (isempty(P) == false)
			# Update T -
			T = vcat(TNEW,P);
		end
    end # main -

	# ok, check for indepedent pathways -
	rows_to_keep_set = compute_convex_independent_rows(T)

	# grab: which rows should we keep?
    T_IND = T[(rows_to_keep_set |> collect),:]

    # return -
    return T_IND

end # function -

function compute_convex_independent_rows(P)::Set{Int64}

	# ok, so we have a new set of possible pathways, check for independent pathways -
	(NRP,_) = size(P)
	IPA = Matrix{Float64}(I,NRP,NRP)
	for i ∈ 1:NRP

		# get rᵢ -
		rᵢ = P[i,:]

		# compare rᵢ against rⱼ i neq j
		for j ∈ 1:NRP
			if (i != j)

				# get rⱼ -
				rⱼ = P[j,:]

				# check: is convex independent?
				if (is_convex_independent(rᵢ,rⱼ) == true)
					IPA[i,j] = 1.0
				else
					IPA[i,j] = 0.0
				end
			end
		end # inner
	end # outer

	# rows we should keep -
	rows_to_keep_set = Set{Int64}()
	for i ∈ 1:NRP
		sᵢ = sum(IPA[i,:])
		if (sᵢ == NRP)
			push!(rows_to_keep_set,i)
		end
	end

	# return -
	return rows_to_keep_set
end

function is_convex_independent(rᵢ,rⱼ)::Bool
	idx_zero_rᵢ = findall(x->x==0.0,rᵢ)
	idx_zero_rⱼ = findall(x->x==0.0,rⱼ)
	return issubset(idx_zero_rᵢ,idx_zero_rⱼ) ? false : true
end


text:
# gravity
get_gravity(g::T) where T <: Real = SVector{3,T}([0.0; 0.0; g])
get_gravity(g::Vector{T}) where T = SVector{3,T}(g)
get_gravity(g::SVector) = g


text:
doc"""
    FDist(ν1,ν2)

The *F distribution* has probability density function

$f(x; \nu_1, \nu_2) = \frac{1}{x B(\nu_1/2, \nu_2/2)}
\sqrt{\frac{(\nu_1 x)^{\nu_1} \cdot \nu_2^{\nu_2}}{(\nu_1 x + \nu_2)^{\nu_1 + \nu_2}}},
\quad x>0$

It is related to the [`Chisq`](:func:`Chisq`) distribution via the property that if $X_1
\sim \operatorname{Chisq}(\nu_1)$ and $X_2 \sim \operatorname{Chisq}(\nu_2)$, then
$(X_1/\nu_1) / (X_2 / \nu_2) \sim FDist(\nu_1, \nu_2)`.

```julia
FDist(d1, d2)     # F-Distribution with parameters d1 and d2

params(d)         # Get the parameters, i.e. (d1, d2)
```

External links

* [F distribution on Wikipedia](http://en.wikipedia.org/wiki/F-distribution)
    """
immutable FDist{T<:Real} <: ContinuousUnivariateDistribution
    ν1::T
    ν2::T

    function FDist(ν1::T, ν2::T)
        @check_args(FDist, ν1 > zero(ν1) && ν2 > zero(ν2))
        new(ν1, ν2)
    end
end

FDist{T<:Real}(ν1::T, ν2::T) = FDist{T}(ν1, ν2)
FDist(ν1::Integer, ν2::Integer) = FDist(Float64(ν1), Float64(ν2))
FDist(ν1::Real, ν2::Real) = FDist(promote(ν1, ν2)...)

@distr_support FDist 0.0 Inf

#### Conversions
function convert{T <: Real, S <: Real}(::Type{FDist{T}}, ν1::S, ν2::S)
    FDist(T(ν1), T(ν2))
end
function convert{T <: Real, S <: Real}(::Type{FDist{T}}, d::FDist{S})
    FDist(T(d.ν1), T(d.ν2))
end

#### Parameters

params(d::FDist) = (d.ν1, d.ν2)
@inline partype{T<:Real}(d::FDist{T}) = T

#### Statistics

mean{T<:Real}(d::FDist{T}) = (ν2 = d.ν2; ν2 > 2 ? ν2 / (ν2 - 2) : T(NaN))

function mode{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    ν1 > 2 ? ((ν1 - 2)/ν1) * (ν2 / (ν2 + 2)) : zero(T)
end

function var{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    ν2 > 4 ? 2ν2^2 * (ν1 + ν2 - 2) / (ν1 * (ν2 - 2)^2 * (ν2 - 4)) : T(NaN)
end

function skewness{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    if ν2 > 6
        return (2ν1 + ν2 - 2) * sqrt(8(ν2 - 4)) / ((ν2 - 6) * sqrt(ν1 * (ν1 + ν2 - 2)))
    else
        return T(NaN)
    end
end

function kurtosis{T<:Real}(d::FDist{T})
    (ν1, ν2) = params(d)
    if ν2 > 8
        a = ν1 * (5ν2 - 22) * (ν1 + ν2 - 2) + (ν2 - 4) * (ν2 - 2)^2
        b = ν1 * (ν2 - 6) * (ν2 - 8) * (ν2 - 2)
        return 12a / b
    else
        return T(NaN)
    end
end

function entropy(d::FDist)
    (ν1, ν2) = params(d)
    hν1 = ν1/2
    hν2 = ν2/2
    hs = (ν1 + ν2)/2
    return log(ν2 / ν1) + lgamma(hν1) + lgamma(hν2) - lgamma(hs) +
        (1 - hν1) * digamma(hν1) + (-1 - hν2) * digamma(hν2) +
        hs * digamma(hs)
end

#### Evaluation & Sampling

@_delegate_statsfuns FDist fdist ν1 ν2

rand(d::FDist) = StatsFuns.RFunctions.fdistrand(d.ν1, d.ν2)


text:
"""
    solve_hadc(problem::HADC, gradient::ScalarGradient)

Compute the ADC using a homogenized ADC model (HADC). This is currently only implemented for
scalar gradients.
"""
function solve(problem::HADC, gradient::ScalarGradient)
    (; model, matrices, odesolver, reltol, abstol) = problem
    (; mesh, D) = model
    (; M_cmpts, S_cmpts, G, volumes) = matrices

    f = gradient.profile
    dir = gradient.dir
    TE = echotime(f)

    # Deduce sizes
    ncompartment = length(D)

    # Initial conditions
    ω₀ = zeros.(size.(mesh.points, 2))

    # Time dependent ODE function
    function Mdω!(dω, ω, p, t)
        # @show t
        (; mS, f, surfint) = p
        mul!(dω, mS, ω)
        dω .+= integral(f, t) .* surfint
    end

    # Create ODE function and Jacobian from matrices
    Jac!(J, _, p, t) = (J .= p.mS)

    # Allocate output array
    adc = zeros(ncompartment)

    # Iterate over compartments and gradient sequences and directions
    for icmpt = 1:ncompartment
        # Free diffusivity in gradient direction
        D₀ = dir' * D[icmpt] * dir

        # Surface integrals in gradient direction
        surfint = G[icmpt] * (D[icmpt] * dir)

        p = (; mS = -S_cmpts[icmpt], f, surfint)

        odefunction = ODEFunction(
            Mdω!,
            jac = Jac!,
            jac_prototype = p.mS,
            mass_matrix = M_cmpts[icmpt],
        )
        odeproblem = ODEProblem(odefunction, ω₀[icmpt], (0, TE), p, progress = false)

        # Solve ODE, keep all time steps (for integral)
        sol = OrdinaryDiffEq.solve(odeproblem, odesolver, reltol = reltol, abstol = abstol)

        # Integral over compartment boundary
        a, = quadgk(t -> integral(f, t) * (surfint' * sol(t)), 0, TE)

        # HADC (free diffusivity minus correction)
        adc[icmpt] = D₀ - a / volumes[icmpt] / int_F²(f)
    end

    adc
end


text:
@doc doc"""
    GroupOperationAction(group::AbstractGroupManifold, AD::ActionDirection = LeftAction())

Action of a group upon itself via left or right translation.
"""
struct GroupOperationAction{G,AD} <: AbstractGroupAction{AD}
    group::G
end

function GroupOperationAction(
    G::AbstractGroupManifold,
    ::TAD = LeftAction(),
) where {TAD<:ActionDirection}
    return GroupOperationAction{typeof(G),TAD}(G)
end

function show(io::IO, A::GroupOperationAction)
    print(io, "GroupOperationAction($(A.group), $(direction(A)))")
end

base_group(A::GroupOperationAction) = A.group

g_manifold(A::GroupOperationAction) = A.group

function switch_direction(A::GroupOperationAction)
    return GroupOperationAction(A.group, switch_direction(direction(A)))
end

apply(A::GroupOperationAction, a, x) = translate(A.group, a, x, direction(A))

apply!(A::GroupOperationAction, y, a, x) = translate!(A.group, y, a, x, direction(A))

function inverse_apply(A::GroupOperationAction, a, x)
    return inverse_translate(A.group, a, x, direction(A))
end

function inverse_apply!(A::GroupOperationAction, y, a, x)
    return inverse_translate!(A.group, y, a, x, direction(A))
end

function apply_diff(A::GroupOperationAction, a, x, v)
    return translate_diff(A.group, a, x, v, direction(A))
end

function apply_diff!(A::GroupOperationAction, vout, a, x, v)
    return translate_diff!(A.group, vout, a, x, v, direction(A))
end

function inverse_apply_diff(A::GroupOperationAction, a, x, v)
    return inverse_translate_diff(A.group, a, x, v, direction(A))
end

function inverse_apply_diff!(A::GroupOperationAction, vout, a, x, v)
    return inverse_translate_diff!(A.group, vout, a, x, v, direction(A))
end

function optimal_alignment(A::GroupOperationAction, x1, x2)
    return inverse_apply(switch_direction(A), x1, x2)
end

function optimal_alignment!(A::GroupOperationAction, y, x1, x2)
    return inverse_apply!(switch_direction(A), y, x1, x2)
end

function center_of_orbit(
    A::GroupOperationAction,
    pts::AbstractVector,
    q,
    mean_method::AbstractEstimationMethod,
)
    m = mean(A.group, pts, mean_method)
    return inverse_apply(switch_direction(A), q, m)
end
function center_of_orbit(A::GroupOperationAction, pts::AbstractVector, q)
    m = mean(A.group, pts)
    return inverse_apply(switch_direction(A), q, m)
end


text:
#
# Advent of Code
# Day 4: Security Through Obscurity
#

println("Advent of Code")
println("Day 4: Security Through Obscurity")
println("")

function pairSorter(pair::Pair{Char,Int})
  fullString = string(pair[2]) * string(pair[1])
  return fullString
end

function parseRooms(inputDataLines)
  total = 0

  for line in inputDataLines
    letterDict = Dict{Char, Int}()

    inSector = false
    sectorID = ""

    inChecksum = false
    checksum = ""

    # evaluate all the letters one by one
    for letter in line
      # check if we have [A-za-z]
      if isalpha(letter)
        # if we've started the checksum, add the letter to the checksum
        if inChecksum
          checksum = checksum * string(letter)
        # if not, and we haven't seen this letter before, start the count
        elseif !haskey(letterDict, letter)
          letterDict[letter] = 1
        # otherwise, increment the count
        else
          letterDict[letter] += 1
        end
      # but if we have a number, we must be in the sectorID
      elseif isnumber(letter)
        # check if we've been in the sectorID yet
        if !inSector
          # we're just starting building the sector string
          sectorID = string(letter)
          inSector = true
        else
          # or, add to the sector string
          sectorID = sectorID * string(letter)
        end
      # If it's not a letter or number, maybe we're starting the checksum
      elseif letter == '['
        if !inChecksum
          inChecksum = true
        end
      end
    end

    # generate the calculated checksum
    generatedChecksum = ""
    letterList = Vector{Pair{Char,Int}}()

    # we need to reverse the numbers for the sort
    lineLength = length(line)

    # get the key, value pairs into a list
    for kv in letterDict
      kv = Pair(kv[1], lineLength - kv[2])
      push!(letterList, kv)
    end

    # sort the list by the values, then the key
    sort!(letterList,by=pairSorter)

    # grab the first five (if they exist)
    if length(letterList) >= 5
      for item in letterList[1:5]
        generatedChecksum = generatedChecksum * string(item[1])
      end

      # compare the generatedChecksum to the checksum
      # and if it matches, then grab the value of the sectorID
      # and then add it to the total

      if generatedChecksum == checksum
        shift = parse(sectorID)

        total += shift

        # part 2 is to decrypt the room name with a shift cypher
        roomName = ""

        for letter in line
          if isnumber(letter)
            break
          end

          if letter == '-'
            newLetter = ' '
          else
            # all letters are lowercase, thank goodness
            newLetter = Char(rem((Int(letter) - 97 + shift), 26) + 97)
          end

          roomName *= string(newLetter)
        end

        println("Room: $roomName ($sectorID)")

      end
    end
  end

  return total

end

# Data file
# Input data has been saved to this file
inputDataFileName = "day04-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

#inputData = "aaaaa-bbb-z-y-x-123[abxyz]\n"
#inputData *= "a-b-c-d-e-f-g-h-987[abcde]\n"
#inputData *= "not-a-real-room-404[oarel]\n"
#inputData *= "totally-real-room-200[decoy]\n"

inputDataLines = split(strip(inputData), '\n')

part1Total = parseRooms(inputDataLines)

println("Part 1 total: $part1Total")


text:
module JuliaOrgs

using GitHub
using TOML
#using ProgressMeter

export get_orglist

ghauth = GitHub.authenticate(ENV["GITHUB_TOKEN"])

function get_orglist()
    orgnames = Dict{String, Any}()

    registrypath = ENV["HOME"] * "/.julia/registries/General"
    packages = TOML.parsefile(registrypath * "/Registry.toml")["packages"]

    for p in packages
        package_path = p[2]["path"]
        repo_url = TOML.parsefile(registrypath * "/" * package_path * "/Package.toml")["repo"]
        name = split(repo_url, "/")[4]
        orgnames[name] = nothing
    end

    function ghowner(x)
        try
            GitHub.owner(x, auth=myauth)
        catch
            x
        end
    end

    return asyncmap(x->ghowner(x), collect(keys(orgnames)))

end

end


text:
""" Generalized algebraic theories (GATs) in Julia.

GATs generalize (multisorted) algebraic theories by incorporating a fragment of
dependent type theory. They allow type and term constructors to be partially
defined. GATs provide a convenient formal syntax for categorical structures.

At present, this module only supports defining the *signature* of a GAT. In the
future we may allow *axioms* to be expressed as well. Regardless, we will
persist in calling this module "GAT". Signatures are defined using the
`@signature` macro.

There are several useful things to do with a GAT signature:

1. Define *instances* of the signature (models of the theory) in Julia code
   using the `@instance` macro
2. Define *syntax* systems using the `@syntax` macro (see `Syntax` module)
3. Define models of the theory by finite *presentation* using a syntax system
   and the `@present` macro (see `Present` module)

References:

- (Cartmell, 1986, "Generalized algebraic theories and contextual categoies")
- (Pitts, 1995, "Categorical logic", Sec 6: "Dependent types")
"""
module GAT
export @signature, @instance, invoke_term

using Base.Meta: ParseError
using Compat
using AutoHashEquals
using DataStructures: OrderedDict
using Match

using ..Meta

# Data types
############

""" Base type for method stubs in GAT signature.
"""
abstract type Stub end

const Context = OrderedDict{Symbol,Expr0}

""" Type constructor in a GAT.
"""
@auto_hash_equals struct TypeConstructor
  name::Symbol
  params::Vector{Symbol}
  context::Context
  doc::Union{String,Nothing}

  function TypeConstructor(name::Symbol, params::Vector,
                           context::Context, doc=nothing)
    new(name, params, context, doc)
  end
end

""" Term constructor in a GAT.
"""
@auto_hash_equals struct TermConstructor
  name::Symbol
  params::Vector{Symbol}
  typ::Expr0
  context::Context
  doc::Union{String,Nothing}

  function TermConstructor(name::Symbol, params::Vector, typ::Expr0,
                           context::Context, doc=nothing)
    new(name, params, typ, context, doc)
  end
end

""" Signature for a generalized algebraic theory (GAT).
"""
@auto_hash_equals struct Signature
  types::Vector{TypeConstructor}
  terms::Vector{TermConstructor}
end

""" Typeclass = GAT signature + Julia-specific content.
"""
struct Typeclass
  name::Symbol
  type_params::Vector{Symbol}
  signature::Signature
  functions::Vector{JuliaFunction}
end

struct SignatureBinding
  name::Symbol
  params::Vector{Symbol}
end
struct SignatureHead
  main::SignatureBinding
  base::Vector{SignatureBinding}
  SignatureHead(main, base=[]) = new(main, base)
end

# Signatures
############

""" Define a signature of a generalized algebraic theory (GAT).

Three kinds of things can go in the signature body:

1. Type constructors, indicated by the special type `TYPE`, e.g.
   `Hom(X::Ob,Y::Ob)::TYPE`
2. Term constructors, e.g.,
   `id(X::Ob)::Hom(X,X)`
3. Julia functions operating on the term constructors to provide additional
   functionality

A signature can extend existing signatures (at present only one).
"""
macro signature(head, body)
  # Parse signature header.
  head = parse_signature_head(head)
  @assert length(head.base) <= 1 "Multiple signature extension not supported"
  if length(head.base) == 1
    base_name, base_params = head.base[1].name, head.base[1].params
    @assert all(p in head.main.params for p in base_params)
  else
    base_name, base_params = nothing, []
  end

  # Parse signature body: GAT types/terms and extra Julia functions.
  types, terms, functions = parse_signature_body(body)
  signature = Signature(types, terms)
  class = Typeclass(head.main.name, head.main.params, signature, functions)

  # We must generate and evaluate the code at *run time* because the base
  # signature, if specified, is not available at *parse time*.
  expr = :(signature_code($class, $(esc(base_name)), $base_params))
  Expr(:block,
    Expr(:call, esc(:eval), expr),
    :(Core.@__doc__ $(esc(head.main.name))))
end
function signature_code(main_class, base_mod, base_params)
  # Add types/terms/functions from base class, if provided.
  if isnothing(base_mod)
    class = main_class
  else
    base_class = base_mod.class()
    bindings = Dict(zip(base_class.type_params, base_params))
    main_sig = main_class.signature
    base_sig = replace_types(bindings, base_class.signature)
    sig = Signature([base_sig.types; main_sig.types],
                    [base_sig.terms; main_sig.terms])
    functions = [ [ replace_symbols(bindings, f) for f in base_class.functions ];
                  main_class.functions ]
    class = Typeclass(main_class.name, main_class.type_params, sig, functions)
  end
  signature = class.signature

  # Generate module with stub types.
  mod = Expr(:module, true, class.name,
    Expr(:block, [
      # Prevents error about export not being at toplevel.
      # https://github.com/JuliaLang/julia/issues/28991
      LineNumberNode(0);
      Expr(:export, [cons.name for cons in signature.types]...);
      map(gen_abstract_type, signature.types);
      :(class() = $(class));
    ]...))

  # Generate method stubs.
  # (We put them outside the module, so the stub type names must be qualified.)
  bindings = Dict(cons.name => Expr(:(.), class.name, QuoteNode(cons.name))
                  for cons in signature.types)
  fns = interface(class)
  toplevel = [ generate_function(replace_symbols(bindings, f)) for f in fns ]

  # Modules must be at top level:
  # https://github.com/JuliaLang/julia/issues/21009
  Expr(:toplevel, mod, toplevel...)
end

function parse_signature_head(expr::Expr)::SignatureHead
  parse = parse_signature_binding
  @match expr begin
    (Expr(:call, [:(=>), Expr(:tuple, bases), main])
      => SignatureHead(parse(main), map(parse, bases)))
    Expr(:call, [:(=>), base, main]) => SignatureHead(parse(main), [parse(base)])
    _ => SignatureHead(parse(expr))
  end
end

function parse_signature_binding(expr::Expr)::SignatureBinding
  @match expr begin
    Expr(:call, [name::Symbol, params...]) => SignatureBinding(name, params)
    _ => throw(ParseError("Ill-formed signature binding $expr"))
  end
end

""" Parse the body of a GAT signature declaration.
"""
function parse_signature_body(expr::Expr)
  @assert expr.head == :block
  types = OrderedDict{Symbol,TypeConstructor}()
  terms = TermConstructor[]
  funs = JuliaFunction[]
  for elem in strip_lines(expr).args
    head = last(parse_docstring(elem)).head
    if head in (:(::), :call)
      cons = parse_constructor(elem)
      if isa(cons, TypeConstructor)
        if haskey(types, cons.name)
          throw(ParseError("Duplicate type constructor $elem"))
        else
          types[cons.name] = cons
        end
      else
        push!(terms, cons)
      end
    elseif head in (:(=), :function)
      push!(funs, parse_function(elem))
    else
      throw(ParseError("Ill-formed signature element $elem"))
    end
  end
  return (collect(values(types)), terms, funs)
end

""" Julia functions for type parameter accessors.
"""
function accessors(sig::Signature)::Vector{JuliaFunction}
  vcat(map(accessors, sig.types)...)
end
function accessors(cons::TypeConstructor)::Vector{JuliaFunction}
  [ JuliaFunction(Expr(:call, param, Expr(:(::), cons.name)),
                  strip_type(cons.context[param]))
    for param in cons.params ]
end

""" Julia functions for term constructors of GAT.
"""
function constructors(sig::Signature)::Vector{JuliaFunction}
  [ constructor(cons, sig) for cons in sig.terms ]
end
function constructor(cons::TermConstructor, sig::Signature)::JuliaFunction
  arg_names = cons.params
  arg_types = [ strip_type(cons.context[name]) for name in arg_names ]
  args = [ Expr(:(::), name, typ) for (name,typ) in zip(arg_names, arg_types) ]
  return_type = strip_type(cons.typ)

  call_expr = Expr(:call, cons.name, args...)
  if !any(has_type(sig, typ) for typ in arg_types)
    call_expr = add_type_dispatch(call_expr, return_type)
  end
  JuliaFunction(call_expr, return_type)
end

""" Complete set of Julia functions for a type class.
"""
function interface(class::Typeclass)::Vector{JuliaFunction}
  [ accessors(class.signature);
    constructors(class.signature);
    class.functions ]
end

""" Get type constructor by name.

Unlike term constructors, type constructors cannot be overloaded, so there is at
most one type constructor with a given name.
"""
function get_type(sig::Signature, name::Symbol)::TypeConstructor
  indices = findall(cons -> cons.name == name, sig.types)
  @assert length(indices) == 1
  sig.types[indices[1]]
end
function has_type(sig::Signature, name::Symbol)::Bool
  findfirst(cons -> cons.name == name, sig.types) != nothing
end

""" Add a type-valued first argument to a Julia function signature.

We need this to avoid ambiguity in method dispatch when a term constructor has
no arguments (e.g., `munit()`) or more generally has no arguments that are types
in the signature (e.g., object generators in a category).

The fundamental reason why these cases must be treated specially is that Julia
does not (yet?) support
[dispatching on return type](https://github.com/JuliaLang/julia/issues/19206).
"""
function add_type_dispatch(call_expr::Expr, type_expr::Expr0)::Expr
  @match call_expr begin
    (Expr(:call, [name, args...]) =>
      Expr(:call, name, :(::Type{$type_expr}), args...))
    _ => throw(ParseError("Ill-formed call expression $call_expr"))
  end
end

# GAT expressions
#################

""" Parse a raw expression in a GAT.

A "raw expression" is a just composition of function and constant symbols.
"""
function parse_raw_expr(expr)
  @match expr begin
    Expr(:call, args) => map(parse_raw_expr, args)
    head::Symbol => nothing
    _ => throw(ParseError("Ill-formed raw expression $expr"))
  end
  expr # Return the expression unmodified. This function just checks syntax.
end

""" Parse context for term or type in a GAT.
"""
function parse_context(expr::Expr)::Context
  context = Context()
  args = expr.head == :tuple ? expr.args : [ expr ]
  for arg in args
    name, typ = @match arg begin
      Expr(:(::), [name::Symbol, typ]) => (name, parse_raw_expr(typ))
      _ => throw(ParseError("Ill-formed context expression $expr"))
    end
    push_context!(context, name, typ)
  end
  return context
end
function push_context!(context, name, expr)
  if haskey(context, name)
    throw(ParseError("Name $name already defined"))
  end
  context[name] = expr
end

""" Parse type or term constructor in a GAT.
"""
function parse_constructor(expr::Expr)::Union{TypeConstructor,TermConstructor}
  # Context is optional.
  doc, expr = parse_docstring(expr)
  cons_expr, context = @match expr begin
    Expr(:call, [:<=, inner, context]) => (inner, parse_context(context))
    _ => (expr, Context())
  end

  # Allow abbreviated syntax where tail of context is included in parameters.
  function parse_param(param::Expr0)::Symbol
    @match param begin
      Expr(:(::), [name::Symbol, typ]) => begin
        push_context!(context, name, parse_raw_expr(typ))
        name
      end
      name::Symbol => name
      _ => throw(ParseError("Ill-formed type/term parameter $param"))
    end
  end

  @match cons_expr begin
    (Expr(:(::), [name::Symbol, :TYPE])
      => TypeConstructor(name, [], context, doc))
    (Expr(:(::), [Expr(:call, [name::Symbol, params...]), :TYPE])
      => TypeConstructor(name, map(parse_param, params), context, doc))
    (Expr(:(::), [Expr(:call, [name::Symbol, params...]), typ])
      => TermConstructor(name, map(parse_param, params), parse_raw_expr(typ),
                         context, doc))
    _ => throw(ParseError("Ill-formed type/term constructor $cons_expr"))
  end
end

""" Generate abstract type definition from a GAT type constructor.
"""
function gen_abstract_type(cons::TypeConstructor)::Expr
  stub_name = GlobalRef(GAT, :Stub)
  expr = :(abstract type $(cons.name) <: $stub_name end)
  generate_docstring(expr, cons.doc)
end

""" Replace names of type constructors in a GAT.
"""
function replace_types(bindings::Dict, sig::Signature)::Signature
  Signature([ replace_types(bindings, t) for t in sig.types ],
            [ replace_types(bindings, t) for t in sig.terms ])
end
function replace_types(bindings::Dict, cons::TypeConstructor)::TypeConstructor
  TypeConstructor(replace_symbols(bindings, cons.name), cons.params,
                  replace_types(bindings, cons.context), cons.doc)
end
function replace_types(bindings::Dict, cons::TermConstructor)::TermConstructor
  TermConstructor(cons.name, cons.params,
                  replace_symbols(bindings, cons.typ),
                  replace_types(bindings, cons.context), cons.doc)
end
function replace_types(bindings::Dict, context::Context)::Context
  GAT.Context(((name => @match expr begin
    (Expr(:call, [sym::Symbol, args...]) =>
      Expr(:call, replace_symbols(bindings, sym), args...))
    sym::Symbol => replace_symbols(bindings, sym)
  end) for (name, expr) in context))
end

""" Remove type parameters from dependent type.
"""
function strip_type(expr)::Symbol
  @match expr begin
    Expr(:call, [head::Symbol, args...]) => head
    sym::Symbol => sym
  end
end

# GAT expressions in a signature
################################

""" Expand context variables that occur implicitly in an expression.

Reference: (Cartmell, 1986, Sec 10: 'Informal syntax').
"""
function expand_in_context(expr, params::Vector{Symbol},
                           context::Context, sig::Signature)
  @match expr begin
    Expr(:call, [name::Symbol, args...]) => begin
      expanded = [expand_in_context(e, params, context, sig) for e in args]
      Expr(:call, name, expanded...)
    end
    name::Symbol => begin
      if name in params
        name
      elseif haskey(context, name)
        expand_symbol_in_context(name, params, context, sig)
      else
        error("Name $name missing from context $context")
      end
    end
    _ => throw(ParseError("Ill-formed raw expression $expr"))
  end
end
function expand_symbol_in_context(sym::Symbol, params::Vector{Symbol},
                                  context::Context, sig::Signature)
  # This code expands symbols that occur as direct arguments to type
  # constructors. If there are term constructors in between, it does not work:
  # indeed, it cannot work in general because the term constructors are not
  # necessarily injective. For example, we can expand :X in
  #   (:X => :Ob, :f => :(Hom(X)))
  # but not in
  #   (:X => :Ob, :Y => :Ob, :f => :(Hom(otimes(X,Y))))
  names = collect(keys(context))
  start = findfirst(names .== sym)
  for name in names[start+1:end]
    expr = context[name]
    if isa(expr, Expr) && expr.head == :call && sym in expr.args[2:end]
      cons = get_type(sig, expr.args[1])
      accessor = cons.params[findfirst(expr.args[2:end] .== sym)]
      expanded = Expr(:call, accessor, name)
      return expand_in_context(expanded, params, context, sig)
    end
  end
  error("Name $sym does not occur explicitly among $params in context $context")
end

""" Expand context variables that occur implicitly in the type expression
of a term constructor.
"""
function expand_term_type(cons::TermConstructor, sig::Signature)
  isa(cons.typ, Symbol) ? cons.typ :
    expand_in_context(cons.typ, cons.params, cons.context, sig)
end

""" Implicit equations defined by a context.

This function allows a generalized algebraic theory (GAT) to be expressed as
an essentially algebraic theory, i.e., as partial functions whose domains are
defined by equations.

References:
 - (Cartmell, 1986, Sec 6: "Essentially algebraic theories and categories with
    finite limits")
 - (Freyd, 1972, "Aspects of topoi")
"""
function equations(context::Context, sig::Signature)::Vector{Pair}
  # The same restrictions as `expand_symbol_in_context` apply here.
  eqs = Pair[]
  names = collect(keys(context))
  for (start, var) in enumerate(names)
    for name in names[start+1:end]
      expr = context[name]
      if isa(expr, Symbol) && !has_type(sig, expr)
        # If the constructor is a symbol and there isn't a matching type in
        # the signature, assume it's a Julia type. For now, these are
        # completely ignored by the syntax system.
        continue
      end
      expr = isa(expr, Symbol) ? Expr(:call, expr) : expr
      cons = get_type(sig, expr.args[1])
      accessors = cons.params[findall(expr.args[2:end] .== var)]
      append!(eqs, (Expr(:call, a, name) => var for a in accessors))
    end
  end
  eqs
end

""" Implicit equations defined by context, allowing for implicit variables.
"""
function equations(params::Vector{Symbol}, context::Context,
                   sig::Signature)::Vector{Pair}
  eqs = [ (expand_in_context(lhs, params, context, sig) =>
           expand_in_context(rhs, params, context, sig))
          for (lhs, rhs) in equations(context, sig) ]
  # Remove tautologies (expr == expr) resulting from expansions.
  # FIXME: Should we worry about redundancies from the symmetry of equality,
  # i.e., (expr1 == expr2) && (expr2 == expr1)?
  filter(eq -> eq.first != eq.second, eqs)
end

""" Implicit equations for term constructor.
"""
function equations(cons::TermConstructor, sig::Signature)::Vector{Pair}
  equations(cons.params, cons.context, sig)
end

# Instances
###########

""" Define an *instance* of a generalized algebraic theory (GAT).

These are perfectly analagous to instances of a type class in Haskell. See also
the Typeclass.jl library for Julia.
"""
macro instance(head, body)
  # Parse the instance definition.
  head = parse_signature_binding(head)
  functions = parse_instance_body(body)

  # We must generate and evaluate the code at *run time* because the signature
  # module is not defined at *parse time*.
  # Also, we "throw away" any docstring.
  # FIXME: Is there a better place to put the docstring?
  expr = :(instance_code($(esc(head.name)), $(esc(head.params)), $functions))
  Expr(:block,
    Expr(:call, esc(:eval), expr),
    :(Core.@__doc__ abstract type $(esc(gensym(:instance_doc))) end)) # /dev/null
end
function instance_code(mod, instance_types, instance_fns)
  code = Expr(:block)
  class = mod.class()
  bindings = Dict(zip(class.type_params, instance_types))
  bound_fns = [ replace_symbols(bindings, f) for f in interface(class) ]
  bound_fns = OrderedDict(parse_function_sig(f) => f for f in bound_fns)
  instance_fns = Dict(parse_function_sig(f) => f for f in instance_fns)
  for (sig, f) in bound_fns
    if haskey(instance_fns, sig)
      f_impl = instance_fns[sig]
    elseif !isnothing(f.impl)
      f_impl = f
    else
      error("Method $(f.call_expr) not implemented in $(class.name) instance")
    end
    push!(code.args, generate_function(f_impl))
  end
  return code
end

""" Parse the body of a GAT instance definition.
"""
function parse_instance_body(expr::Expr)::Vector{JuliaFunction}
  @match strip_lines(expr) begin
    Expr(:block, args) => map(parse_function, args)
    _ => throw(ParseEror("Ill-formed instance definition"))
  end
end

""" Invoke a term constructor by name on an instance.

This method provides reflection for GAT signatures. In everyday use the generic
method for the constructor should be called directly, not through this function.

Cf. Julia's builtin `invoke()` function.
"""
function invoke_term(signature_module::Module, instance_types::Tuple,
                     constructor_name::Symbol, args...)
  # Get the corresponding Julia method from the parent module.
  method = getfield(parentmodule(signature_module), constructor_name)
  args = collect(Any, args)

  # Add dispatch on return type, if necessary.
  if !any(typeof(arg) <: typ for typ in instance_types for arg in args)
    # Case 1: Name refers to type constructor, e.g., generator constructor
    # in syntax system.
    signature = signature_module.class().signature
    index = findfirst(cons -> cons.name == constructor_name, signature.types)
    if isnothing(index)
      # Case 2: Name refers to term constructor.
      # FIXME: Terms constructors can be overloaded, so there may be multiple
      # term constructors with the same name. Distinguishing them requires type
      # inference. I am punting on that right now.
      constructor = signature.terms[
        findfirst(cons -> cons.name == constructor_name, signature.terms)]
      return_name = strip_type(constructor.typ)
      index = findfirst(cons -> cons.name == return_name, signature.types)
    end
    insert!(args, 1, instance_types[index])
  end

  # Invoke the method!
  method(args...)
end

end


text:
# Autogenerated wrapper script for mimalloc_jll for x86_64-linux-gnu
export libmimalloc

JLLWrappers.@generate_wrapper_header("mimalloc")
JLLWrappers.@declare_library_product(libmimalloc, "libmimalloc.so.2.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmimalloc,
        "lib/libmimalloc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# call with julia --track-allocation=user test/runprofiler.jl

using VortexCollisions

include("test_functions.jl")

function run_profiler()
    M = 64
    N = 64

    ℳcut = M
    𝒩cut = N

    grid = Grid2d(M,N)
    ft   = FourierTransform(grid; ℳcut=ℳcut, 𝒩cut=𝒩cut)
    op   = FokkerPlanckOperator(grid, ft)

    u    = get_field(grid)
    divJ = zeros(u)

    evaluate_function_on_grid(grid, u_test, u)

    collision_operator!(op, u, divJ)

    Profile.clear_malloc_data()

    collision_operator!(op, u, divJ)

end

run_profiler()


text:
function init_face_data(elem::Tri, N; quad_rule_face = gauss_quad(0,0,N))
    # Nodes on faces, and face node coordinate
    r1D, w1D = quad_rule_face
    e = ones(size(r1D)) # vector of all ones
    z = zeros(size(r1D)) # vector of all zeros
    rf,sf = map_face_nodes(elem,r1D)
    wf = vec(repeat(w1D,3,1));
    nrJ = [z; e; -e]
    nsJ = [-e; e; z]
    return rf,sf,wf,nrJ,nsJ
end

function init_face_data(elem::Quad, N; quad_rule_face = gauss_quad(0,0,N))
    r1D,w1D = quad_rule_face
    e = ones(size(r1D))
    z = zeros(size(r1D))
    rf,sf = map_face_nodes(elem,r1D)
    wf = vec(repeat(w1D,4,1)); # 4 faces
    nrJ = [z; e; z; -e]
    nsJ = [-e; z; e; z]
    return rf,sf,wf,nrJ,nsJ
end

function init_face_data(elem::Hex, N; quad_rule_face = quad_nodes(Quad(), N))
    rquad, squad, wquad = quad_rule_face
    e = ones(size(rquad))
    zz = zeros(size(rquad))
    rf,sf,tf = map_face_nodes(elem,rquad,squad)
    Nfaces = 6
    wf = vec(repeat(wquad,Nfaces,1));
    nrJ = [-e; e; zz;zz; zz;zz]
    nsJ = [zz;zz; -e; e; zz;zz]
    ntJ = [zz;zz; zz;zz; -e; e]
    return rf,sf,tf,wf,nrJ,nsJ,ntJ
end

function init_face_data(elem::Tet, N; quad_rule_face = quad_nodes(Tri(), N))
    rquad,squad,wquad = quad_rule_face
    e = ones(size(rquad))
    zz = zeros(size(rquad))
    rf,sf,tf = map_face_nodes(elem,rquad,squad)
    Nfaces = 4
    wf = vec(repeat(wquad,Nfaces,1));
    nrJ = [zz; e; -e; zz]
    nsJ = [-e; e; zz; zz]
    ntJ = [zz; e; zz; -e]
    return rf,sf,tf,wf,nrJ,nsJ,ntJ
end

"""
    RefElemData(elem::Line, N;
                quad_rule_vol = quad_nodes(elem, N+1))
    RefElemData(elem::Union{Tri, Quad}, N;
                 quad_rule_vol = quad_nodes(elem, N),
                 quad_rule_face = gauss_quad(0, 0, N))
    RefElemData(elem::Union{Hex, Tet}, N;
                 quad_rule_vol = quad_nodes(elem, N),
                 quad_rule_face = quad_nodes(Quad(), N))
    RefElemData(elem; N, kwargs...) # version with keyword args

Constructor for `RefElemData` for different element types.
"""
function RefElemData(elem::Line, approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N+1),
                     Nplot = 10)

    fv = face_vertices(elem)

    # Construct matrices on reference elements
    r = nodes(elem,N)
    Fmask = [1 N+1]
    VDM = vandermonde(elem, N, r)
    Dr = grad_vandermonde(elem, N, r)/VDM

    V1 = vandermonde(elem,1,r)/vandermonde(elem,1,[-1;1])

    rq,wq = quad_rule_vol
    Vq = vandermonde(elem,N,rq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))

    rf = [-1.0;1.0]
    nrJ = [-1.0;1.0]
    wf = [1.0;1.0]
    Vf = vandermonde(elem,N,rf)/VDM
    LIFT = M\(Vf') # lift matrix

    # plotting nodes
    rp = equi_nodes(elem,Nplot)
    Vp = vandermonde(elem,N,rp)/VDM

    return RefElemData(elem,approxType,N,fv,V1,
                       tuple(r),VDM,vec(Fmask),
                       Nplot,tuple(rp),Vp,
                       tuple(rq),wq,Vq,
                       tuple(rf),wf,Vf,tuple(nrJ),
                       M,Pq,tuple(Dr),LIFT)
end

function RefElemData(elem::Union{Tri, Quad},  approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N),
                     quad_rule_face = quad_nodes(face_type(elem), N),
                     Nplot = 10)

    fv = face_vertices(elem) # set faces for triangle

    # Construct matrices on reference elements
    r,s = nodes(elem,N)
    Fmask = hcat(find_face_nodes(elem,r,s)...)

    VDM,Vr,Vs = basis(elem,N,r,s)
    Dr = Vr/VDM
    Ds = Vs/VDM

    # low order interpolation nodes
    r1,s1 = nodes(elem,1)
    V1 = vandermonde(elem,1,r,s)/vandermonde(elem,1,r1,s1)

    rf,sf,wf,nrJ,nsJ = init_face_data(elem, N, quad_rule_face = quad_rule_face)

    rq,sq,wq = quad_rule_vol
    Vq = vandermonde(elem,N,rq,sq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))

    Vf = vandermonde(elem,N,rf,sf)/VDM # interpolates from nodes to face nodes
    LIFT = M\(Vf'*diagm(wf)) # lift matrix used in rhs evaluation

    # plotting nodes
    rp, sp = equi_nodes(elem,Nplot)
    Vp = vandermonde(elem,N,rp,sp)/VDM

    Drs = (Dr,Ds)

    return RefElemData(elem, approxType, N, fv, V1,
                       tuple(r, s), VDM, vec(Fmask),
                       Nplot, tuple(rp, sp), Vp,
                       tuple(rq, sq), wq, Vq,
                       tuple(rf, sf), wf, Vf, tuple(nrJ, nsJ),
                       M, Pq, Drs, LIFT)
end

function RefElemData(elem::Tet, approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N),
                     quad_rule_face = quad_nodes(face_type(elem), N),
                     Nplot = 10)

    fv = face_vertices(elem)

    # Construct matrices on reference elements
    r,s,t = nodes(elem,N)
    Fmask = hcat(find_face_nodes(elem,r,s,t)...)
    VDM,Vr,Vs,Vt = basis(elem,N,r,s,t)
    Dr,Ds,Dt = (A->A/VDM).((Vr,Vs,Vt))

    # low order interpolation nodes
    r1,s1,t1 = nodes(elem,1)
    V1 = vandermonde(elem,1,r,s,t)/vandermonde(elem,1,r1,s1,t1)

    # Nodes on faces, and face node coordinate
    rf,sf,tf,wf,nrJ,nsJ,ntJ = init_face_data(elem, N, quad_rule_face = quad_rule_face)

    # quadrature nodes - build from 1D nodes.
    rq,sq,tq,wq = quad_rule_vol
    Vq = vandermonde(elem,N,rq,sq,tq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))

    Vf = vandermonde(elem,N,rf,sf,tf)/VDM
    LIFT = M\(Vf'*diagm(wf))

    # plotting nodes
    rp,sp,tp = equi_nodes(elem,Nplot)
    Vp = vandermonde(elem,N,rp,sp,tp)/VDM

    Drst = (Dr,Ds,Dt)

    return RefElemData(elem, approxType, N, fv, V1,
                       tuple(r, s, t), VDM, vec(Fmask),
                       Nplot, tuple(rp, sp, tp), Vp,
                       tuple(rq, sq, tq), wq, Vq,
                       tuple(rf, sf, tf), wf, Vf, tuple(nrJ, nsJ, ntJ),
                       M, Pq, Drst, LIFT)
end

# specialize constructor for `Hex` to allow for higher polynomial degrees `N`
function RefElemData(elem::Hex, approxType::Polynomial, N;
                     quad_rule_vol = quad_nodes(elem, N),
                     quad_rule_face = quad_nodes(face_type(elem), N),
                     Nplot = 10)

    fv = face_vertices(elem)

    # Construct matrices on reference elements
    r, s, t = nodes(elem, N)
    Fmask = hcat(find_face_nodes(elem, r, s, t)...)

    # construct 1D operator for faster Kronecker solves
    r1D = nodes(Line(), N)
    rq1D, wq1D = quad_nodes(Line(), N)
    VDM_1D = vandermonde(Line(), N, r1D)
    Vq1D = vandermonde(Line(), N, rq1D) / VDM_1D
    invVDM_1D = inv(VDM_1D)
    invM_1D = VDM_1D * VDM_1D'
    M1D = Vq1D' * diagm(wq1D) * Vq1D

    # form kronecker products of multidimensional matrices to invert/multiply
    M = kronecker(M1D, M1D, M1D)
    VDM = kronecker(VDM_1D, VDM_1D, VDM_1D)
    invVDM = kronecker(invVDM_1D, invVDM_1D, invVDM_1D)
    invM = kronecker(invM_1D, invM_1D, invM_1D)

    _, Vr, Vs, Vt = basis(elem, N, r, s, t)
    Dr, Ds, Dt = (A -> A * invVDM).((Vr, Vs, Vt))

    # low order interpolation nodes
    r1, s1, t1 = nodes(elem, 1)
    V1 = vandermonde(elem, 1, r, s, t) / vandermonde(elem, 1, r1, s1, t1)

    # Nodes on faces, and face node coordinate
    rf, sf, tf, wf, nrJ, nsJ, ntJ = init_face_data(elem, N, quad_rule_face=quad_rule_face)

    # quadrature nodes - build from 1D nodes.
    rq, sq, tq, wq = quad_rule_vol
    Vq = vandermonde(elem, N, rq, sq, tq) * invVDM
    # M = Vq' * diagm(wq) * Vq
    Pq = invM * (Vq' * diagm(wq))

    Vf = vandermonde(elem, N, rf, sf, tf) * invVDM
    LIFT = invM * (Vf' * diagm(wf))

    # plotting nodes
    rp, sp, tp = equi_nodes(elem, Nplot)
    # Vp = vandermonde(elem, N, rp, sp, tp) * invVDM
    rp1D = LinRange(-1, 1, Nplot + 1)
    Vp1D = vandermonde(Line(), N, rp1D) / VDM_1D
    Vp = kronecker(Vp1D, Vp1D, Vp1D)

    Drst = (Dr, Ds, Dt)

    return RefElemData(elem, approxType, N, fv, V1,
                       tuple(r, s, t), VDM, vec(Fmask),
                       Nplot, tuple(rp, sp, tp), Vp,
                       tuple(rq, sq, tq), wq, Vq,
                       tuple(rf, sf, tf), wf, Vf, tuple(nrJ, nsJ, ntJ),
                       M, Pq, Drst, LIFT)
end


text:
module GitHub

import Base.show

using Requests
import JSON
using HttpCommon

abstract GitHubType
abstract Owner <: GitHubType

# types
export User,
       Organization,
       Repo,
       Issue,
       HttpError,
       AuthError,
       StatsError

# methods
export authenticate,
       set_api_endpoint,
       set_web_endpoint,
       user,
       star,
       unstar,
       stargazers,
       starred,
       forks,
       fork,
       contributors,
       commit_activity,
       code_frequency,
       participation,
       punch_card,
       collaborators,
       iscollaborator,
       add_collaborator,
       remove_collaborator,
       watchers,
       watched,
       watching,
       watch,
       unwatch,
       followers,
       following,
       org,
       orgs,
       repo,
       issue,
       create_issue,
       edit_issue,
       issues

include("utils.jl")
include("endpoint.jl")
include("error.jl")
include("auth.jl")
include("users.jl")
include("organizations.jl")
include("repos.jl")
include("issues.jl")
include("starring.jl")
include("forks.jl")
include("statistics.jl")
include("collaborators.jl")
include("watching.jl")

end


text:
using Documenter
using SymmetricSparseMatrices

makedocs(
    sitename = "SymmetricSparseMatrices",
    format = Documenter.HTML(),
    modules = [SymmetricSparseMatrices]
)

deploydocs(repo = "github.com/victorsndvg/SymmetricSparseMatrices.jl.git")


text:
for abstract type AbstractVector{T}
where T
  length(using) -> Int
  index(using, Int) -> T
end


text:
using Flux, Gym, Printf, Zygote
using Zygote: @adjoint
using Flux.Optimise: update!
using Statistics: mean
#using CuArrays

import Base.sign

#Load game environment

env = make("CartPole-v0")
reset!(env)

#ctx = Ctx(env)

#display(ctx.s)
#using Blink# when not on Juno
#body!(Blink.Window(), ctx.s)

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(env._env.state)
ACTION_SIZE = length(env._env.action_space)
MAX_TRAIN_REWARD = env._env.x_threshold * env._env.θ_threshold_radians
SEQ_LEN = 8

# Optimiser params
η = 3f-2
# ------------------------------ Model Architecture ----------------------------
sign(x) = Base.sign.(x)
@adjoint sign(x) = sign(x), x̄ -> (x̄,)

model = Chain(Dense(STATE_SIZE, 24, relu),
              Dense(24, 48, relu),
              Dense(48, 1, tanh), x->sign(x)) |> gpu

opt = ADAM(η)

action(state) = state |> model |> (model_output) -> (3 .+ model_output) / 2

loss(rewards) = Flux.mse(rewards, MAX_TRAIN_REWARD)

# ----------------------------- Helper Functions -------------------------------

function train_reward(env::EnvWrapper)
    s = env._env.state
    x, ẋ, θ, θ̇  = s
    # Custom reward for training
    # Product of Triangular function over x-axis and θ-axis
    # Min reward = 0, Max reward = env.x_threshold * env.θ_threshold_radians
    x_upper = env._env.x_threshold - x
    x_lower = env._env.x_threshold + x

    r_x     = max(0f0, min(x_upper, x_lower))

    θ_upper = env._env.θ_threshold_radians - θ
    θ_lower = env._env.θ_threshold_radians + θ

    r_θ     = max(0f0, min(θ_upper, θ_lower))

    return r_x * r_θ
end

function μEpisode(env::EnvWrapper)
    l = 0
    for frames ∈ 1:SEQ_LEN
        #render(env, ctx)
        #sleep(0.01)
        a = action(env._env.state)
        s′, r, done, _ = step!(env, a)

        if trainable(env)
            l += loss(train_reward(env))
        end
        game_over(env) && break
    end
    return l
end

function episode!(env::EnvWrapper)
    reset!(env)
    while !game_over(env)
        if trainable(env)
            grads = gradient(()->μEpisode(env), params(model))
            update!(opt, params(model), grads)
        else
            μEpisode(env)
        end
    end
    env.total_reward
end

# -------------------------------- Testing -------------------------------------

function test(env::EnvWrapper)
    score_mean = 0f0
    testmode!(env)
    for _=1:100
        total_reward = episode!(env)
        score_mean += total_reward / 100
    end
    testmode!(env, false)
    return score_mean
end

# ------------------------------ Training --------------------------------------

e = 1
while true
    global e
    total_reward = @sprintf "%6.2f" episode!(env)
    print("Episode: $e | Score: $total_reward | ")

    score_mean = test(env)
    score_mean_str = @sprintf "%6.2f" score_mean
    print("Mean score over 100 test episodes: " * score_mean_str)

    println()

    if score_mean > env.reward_threshold
        println("CartPole-v0 solved!")
    break
    end
    e += 1
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
import Base.show

function show(io::IO, x::Select)
    write(io, "SELECT ",
        join(UTF8String[sprint(io->show(io, x)) for x in x.relation.attrs], ", "))
    show(io, x.from)
end

show(io::IO, x::DataStreamRelation) = fieldset(x).name

function show(io::IO, x::Attr)
    write(io, string(name(x)))
end

function show(io::IO, x::QualifiedAttr)
    write(io, string(x.relation.name), ".", string(name(x.attr)))
end

function show(io::IO, x::Where)
    show(io, x.source)
    write(io, " WHERE ")
    show(io, x.predicate)
end

function show(io::IO, x::GroupBy)
    show(io, x.source)
    write(io, " GRPOUPBY ")
    show(io, x.grouping)
end

function show(io::IO, x::Aggregate)
    show(io, x.source)
    write(io, ' ')
    show(io, x.aggregator)
end

function show(io::IO, agg::Count)
    write(io, "count()")
end

function show(io::IO, agg::Avg)
    write(io, "avg(")
    show(io, agg.attr)
    write(io, ")")
end

function show(io::IO, agg::Sum)
    write(io, "sum(")
    show(io, agg.attr)
    write(io, ")")
end

function show(io::IO, a::QualifiedAttr)
    write(io, string(a.relation.name))
    write(io, '.')
    write(io, string(name(a)))
end

function show{op}(io::IO, a::BinOp{op})
    write(io, '(')
    show(io, a.l)
    write(io, ' ')
    write(io, op)
    write(io, ' ')
    show(io, a.r)
    write(io, ')')
end


text:
function bisection(f::Function, a::Number, b::Number;
    tol::AbstractFloat=1e-5, maxiter::Integer=100)
    fa = f(a)
    fa*f(b) <= 0 || error("No real root in [a,b]")
    i = 0
    local c
    while b-a > tol
        i += 1
        i != maxiter || error("Max iteration exceeded")
        c = (a+b)/2
        fc = f(c)
        if fc == 0
            break
        elseif fa*fc > 0
            a = c  # Root is in the right half of [a,b].
            fa = fc
        else
            b = c  # Root is in the left half of [a,b].
        end
    end
    return c
end


text:
using LinearTransforms, BenchmarkTools

m, n = 100, 100
A = randn(m, n)
prod = v -> A*v
tprod = w -> A'*w
op = LinearTransform(m, n, prod, tprod)
x = randn(n)
y = randn(m)

println("Forward native multiply")
@benchmark A*x
println("Forward operator multiply")
@benchmark op*x

println("Ajoint native multiply")
@benchmark A'*y
println("Adjoint operator multiply")
@benchmark op'*y


text:
# Autogenerated wrapper script for Libgpg_error_jll for aarch64-apple-darwin
export libgpg_error

JLLWrappers.@generate_wrapper_header("Libgpg_error")
JLLWrappers.@declare_library_product(libgpg_error, "@rpath/libgpg-error.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgpg_error,
        "lib/libgpg-error.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()