text:
function finitegrad{T <: Number}(
    f::Function,
    x::T,
    direction::Symbol = :central,
)
    if direction == :forward
        ε = @forwardrule(x)
        return (f(x + ε) - f(x)) / ε
    elseif direction == :reverse
        ε = @forwardrule(x)
        return (f(x) - f(x - ε)) / ε
    elseif direction == :central
        ε = @centralrule(x)
        return (f(x + ε) - f(x - ε)) / (ε + ε)
    elseif direction == :complex
        ε = @complexrule(x)
        return imag(f(x + ε * im)) / ε
    else
        throw(
            ArgumentError(
                "direction must be :forward, :reverse, :central or :complex"
            )
        )
    end
end

function finitegrad!{T <: Number}(
    gr::Vector{T},
    f::Function,
    x::Vector{T},
    direction::Symbol = :central,
)
    n = length(x)
    if direction == :forward
        f_x = f(x)
        for i in 1:n
            ε = @forwardrule(x[i])
            x_i = x[i]
            x[i] = x_i + ε
            f_xpdx = f(x)
            gr[i] = (f_xpdx - f_x) / ε
            x[i] = x_i
        end
    elseif direction == :reverse
        f_x = f(x)
        for i in 1:n
            ε = @forwardrule(x[i])
            x_i = x[i]
            x[i] = x_i - ε
            f_xmdx = f(x)
            gr[i] = (f_x - f_xmdx) / ε
            x[i] = x_i
        end
    elseif direction == :central
        for i in 1:n
            ε = @centralrule(x[i])
            x_i = x[i]
            x[i] = x_i + ε
            f_xpdx = f(x)
            x[i] = x_i - ε
            f_xmdx = f(x)
            gr[i] = (f_xpdx - f_xmdx) / (ε + ε)
            x[i] = x_i
        end
    else
        error("direction must be :forward, :reverse or :central")
    end
    return gr
end

function finitegrad!{T <: Number}(
    gr::Vector{T},
    x_c::Vector{Complex{T}},
    f::Function,
    x::Vector{T},
)
    n = length(x)
    for i in 1:n
        x_c[i] = x[i]
    end
    for i in 1:n
        ε = @complexrule(x[i])
        x_c[i] = x[i] + ε * im
        gr[i] = imag(f(x_c)) / ε
        x_c[i] = x[i]
    end
    return gr
end

function finitegrad{T <: Real}(
    f::Function,
    x::Vector{T},
    direction::Symbol = :central,
)
    gr = similar(x)
    if direction == :complex
        x_c = Array(Complex{T}, length(x))
        finitegrad!(gr, x_c, f, x)
    else
        finitegrad!(gr, f, x, direction)
    end
    return gr
end

function taylor_finite_difference(
    f::Function,
    x::Real,
    direction::Symbol = :central,
    h::Real = 10e-4
)
    if direction == :forward
        f_x = f(x)
        d = 2^3 * (2^2 * (f(x + h) - f_x) - (f(x + 2 * h) - f_x))
        d += - (2^2 * (f(x + 2 * h) - f_x) - (f(x + 4 * h) - f_x))
        d /= 3 * 2^2 * h
    elseif direction == :central
        d = 4^5 * (2^3 * (f(x + h) - f(x - h)) - (f(x + 2 * h) - f(x - 2 * h)))
        d -= 2^3 * (f(x + 4 * h) - f(x - 4 * h)) - (f(x + 8 * h) - f(x - 8 * h))
        d /= (4^5 * (2^4 - 2^2) - (2^6 - 2^4)) * h
    else
        throw(ArgumentError("direction must be :forward or :central"))
    end
    return d
end


text:
# p34.jl - Allen-Cahn eq. u_t = eps*u_xx+u-u^3, u(-1)=-1, u(1)=1
#         (compare p6.jl and p32.jl)

# Differentiation matrix and initial data:
N = 20; (D,x) = cheb(N); D2 = D^2;     # use full-size matrix
D2[[1,N+1],:] = 0;                     # for convenience
eps = 0.01; dt = min(.01,50/(N^4*eps));
t = 0.0; v = @. .53*x + .47*sin(-1.5*pi*x);

# Solve PDE by Euler formula and plot results:
tmax = 100; tplot = 2; nplots = round(Int,tmax/tplot);
plotgap = round(Int,tplot/dt); dt = tplot/plotgap;
xx = -1:.025:1; vv = polyval(polyfit(x,v),xx);
plotdata = [vv zeros(length(xx),nplots)]; tdata = t;
for i = 1:nplots
    for n = 1:plotgap
        t = t+dt; v = v + dt*(eps*D2*(v-x) + v - v.^3);    # Euler
    end
    vv = polyval(polyfit(x,v),xx);
    plotdata[:,i+1] = vv; tdata = [tdata; t];
end
clf();
surf(xx,tdata,plotdata'); grid(true);
xlim(-1,1); ylim(0,tmax); zlim(-1,1);
gca()[:view_init](55,-150);
xlabel("x"); ylabel("t"); zlabel("u");

ext: jl


text:
module WireUp
    export wireup

    using Distributed
    using Sockets
    using PMI

    import Distributed: LPROC

    struct PMIManager <: Distributed.ClusterManager
        kvs::PMI.KVS
    end

    function start_worker(kvs, rank)
        Distributed.init_multi()
        close(stdin) # workers will not use it

        cookie = PMI.get(kvs, "cookie")

        init_worker(cookie)
        interface = IPv4(LPROC.bind_addr)
        if LPROC.bind_port == 0
            port_hint = 9000 + (getpid() % 1000)
            (port, sock) = listenany(interface, UInt16(port_hint))
            LPROC.bind_port = port
        else
            sock = listen(interface, LPROC.bind_port)
        end
        @async while isopen(sock)
            client = accept(sock)
            process_messages(client, client, true)
        end

        PMI.put!(kvs, string(rank), "$(LPROC.bind_port)#$(LPROC.bind_addr)")
        PMI.commit!(kvs)
        PMI.barrier()

        Sockets.nagle(sock, false)
        Sockets.quickack(sock, true)

        try
            # To prevent hanging processes on remote machines, newly launched workers exit if the
            # master process does not connect in time.
            Distributed.check_master_connect()
            while true; wait(); end
        catch err
            print(stderr, "unhandled exception on $(myid()): $(err)\nexiting.\n")
        end

        close(sock)
        PMI.finalize()
        exit(0)
    end

    function Distributed.launch(cm::PMIManager, params, launched, launch_ntfy)
        kvs = cm.kvs
        for wid in 1:(PMI.get_size()-1) # for all workers
            port, addr = split(PMI.get(kvs, string(wid)), "#")
            wc = WorkerConfig()
            wc.io = nothing
            wc.host = addr
            wc.bind_addr = addr
            wc.port = parse(Int, port)
            push!(launched, wc)
            notify(launch_ntfy)
        end
    end

    function Distributed.manage(::PMIManager, id, config, op)
    end

    function start_primary(kvs)
        PMI.barrier() # all workers have now written their port to the KVS, and are available for connection
        addprocs(PMIManager(kvs))
        atexit(PMI.finalize)
    end

    function wireup()
        PMI.init()
        rank = PMI.get_rank()
        kvs = PMI.KVS()
        if rank == 0
            PMI.put!(kvs, "cookie", cluster_cookie())
            PMI.commit!(kvs)
        end
        PMI.barrier()

        if rank == 0
            start_primary(kvs)
        else
            start_worker(kvs, rank)
        end
    end
end


text:
module ImageContrastAdjustment

using ColorVectorSpace
using ImageCore
using ImageTransformations: imresize
# Where possible we avoid a direct dependency to reduce the number of [compat] bounds
using ImageCore.MappedArrays
using Parameters: @with_kw # Same as Base.@kwdef but works on Julia 1.0

# TODO: port HistogramAdjustmentAPI to ImagesAPI
include("HistogramAdjustmentAPI/HistogramAdjustmentAPI.jl")
import .HistogramAdjustmentAPI: AbstractHistogramAdjustmentAlgorithm,
                                adjust_histogram, adjust_histogram!

# TODO Relax this to all image color types
const GenericGrayImage = AbstractArray{<:Union{Number, AbstractGray}}

include("core.jl")
include("build_histogram.jl")
include("algorithms/common.jl")
include("algorithms/adaptive_equalization.jl")
include("algorithms/equalization.jl")
include("algorithms/linear_stretching.jl")
include("algorithms/contrast_stretching.jl")
include("algorithms/gamma_correction.jl")
include("algorithms/matching.jl")
include("algorithms/midway_equalization.jl")
include("compat.jl")

export
    # main types and functions
    AdaptiveEqualization,
    Equalization,
    MidwayEqualization,
    Matching,
    GammaCorrection,
    LinearStretching,
    ContrastStretching,
    build_histogram,
    adjust_histogram,
    adjust_histogram!

end # module


text:
using DataFrames
using TBComponents

lorenz=DataFrame( pop=ones(10),inc=[10,11,12,13,14,15,16,17,30,50])
lorenz[!,:cumpop] = cumsum(lorenz.pop)
lorenz[!,:cuminc] = cumsum(lorenz.inc)

iq = makeinequality(lorenz,:pop,:inc)

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function init_timeteps(h5::HDF5.File, solution::Solution{DT,TT}) where {DT,TT}
    t = create_dataset(h5, "t", TT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    t[1] = timesteps(solution)[0]
end

init_timeteps(sio::SolutionHDF5, args...) = init_timeteps(hdf5(sio), args...)

function init_solution(h5::HDF5.File, solution::SolutionODE{DT,TT,1}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    q[1] = solution.q[0]
end

function init_solution(h5::HDF5.File, solution::SolutionODE{AT,TT,1}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    q[elaxes..., 1] = solution.q[0]
end

function init_solution(h5::HDF5.File, solution::SolutionODE{DT,TT,2}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    q[1, :] = solution.q[0, :]
end

function init_solution(h5::HDF5.File, solution::SolutionODE{AT,TT,2}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin, begin])
    elaxes = axes(solution.q[begin, begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.q, 2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end
end

function init_solution(h5::HDF5.File, solution::SolutionPODE{DT,TT,1}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    p = create_dataset(h5, "p", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    q[1] = solution.q[0]
    p[1] = solution.p[0]
end

function init_solution(h5::HDF5.File, solution::SolutionPODE{AT,TT,1}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    q[elaxes..., 1] = solution.q[0]

    elsize = size(solution.p[begin])
    elaxes = axes(solution.p[begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    p[elaxes..., 1] = solution.p[0]
end

function init_solution(h5::HDF5.File, solution::SolutionPODE{DT,TT,2}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    p = create_dataset(h5, "p", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    q[1, :] = solution.q[0, :]
    p[1, :] = solution.p[0, :]
end

function init_solution(h5::HDF5.File, solution::SolutionPODE{AT,TT,2}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin, begin])
    elaxes = axes(solution.q[begin, begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.q, 2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end

    elsize = size(solution.p[begin, begin])
    elaxes = axes(solution.p[begin, begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.p, 2)
        p[elaxes..., 1, k] = solution.p[0, k]
    end
end

function init_solution(h5::HDF5.File, solution::SolutionDAE{DT,TT,1}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    λ = create_dataset(h5, "λ", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    q[1] = solution.q[0]
    λ[1] = solution.λ[0]
end

function init_solution(h5::HDF5.File, solution::SolutionDAE{AT,TT,1}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    q[elaxes..., 1] = solution.q[0]

    elsize = size(solution.λ[begin])
    elaxes = axes(solution.λ[begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    λ[elaxes..., 1] = solution.λ[0]
end

function init_solution(h5::HDF5.File, solution::SolutionDAE{DT,TT,2}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    λ = create_dataset(h5, "λ", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    q[1, :] = solution.q[0, :]
    λ[1, :] = solution.λ[0, :]
end

function init_solution(h5::HDF5.File, solution::SolutionDAE{AT,TT,2}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin, begin])
    elaxes = axes(solution.q[begin, begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.q, 2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end

    elsize = size(solution.λ[begin, begin])
    elaxes = axes(solution.λ[begin, begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.λ, 2)
        λ[elaxes..., 1, k] = solution.λ[0, k]
    end
end

function init_solution(h5::HDF5.File, solution::SolutionPDAE{DT,TT,1}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    p = create_dataset(h5, "p", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    λ = create_dataset(h5, "λ", DT, ((solution.nt + 1,), (-1,)), chunk = (1,))
    q[1] = solution.q[0]
    p[1] = solution.p[0]
    λ[1] = solution.λ[0]
end

function init_solution(h5::HDF5.File, solution::SolutionPDAE{AT,TT,1}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin])
    elaxes = axes(solution.q[begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    q[elaxes..., 1] = solution.q[0]

    elsize = size(solution.p[begin])
    elaxes = axes(solution.p[begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    p[elaxes..., 1] = solution.p[0]

    elsize = size(solution.λ[begin])
    elaxes = axes(solution.λ[begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt + 1), (elsize..., -1)), chunk = (elsize..., 1))
    λ[elaxes..., 1] = solution.λ[0]
end

function init_solution(h5::HDF5.File, solution::SolutionPDAE{DT,TT,2}) where {DT<:Number,TT}
    q = create_dataset(h5, "q", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    p = create_dataset(h5, "p", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    λ = create_dataset(h5, "λ", DT, ((solution.nt + 1, solution.ni), (-1, solution.ni)), chunk = (1, 1))
    q[1, :] = solution.q[0, :]
    p[1, :] = solution.p[0, :]
    λ[1, :] = solution.λ[0, :]
end

function init_solution(h5::HDF5.File, solution::SolutionPDAE{AT,TT,2}) where {DT,AT<:Array{DT},TT}
    elsize = size(solution.q[begin, begin])
    elaxes = axes(solution.q[begin, begin])
    q = create_dataset(h5, "q", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.q, 2)
        q[elaxes..., 1, k] = solution.q[0, k]
    end

    elsize = size(solution.p[begin, begin])
    elaxes = axes(solution.p[begin, begin])
    p = create_dataset(h5, "p", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.p, 2)
        p[elaxes..., 1, k] = solution.p[0, k]
    end

    elsize = size(solution.λ[begin, begin])
    elaxes = axes(solution.λ[begin, begin])
    λ = create_dataset(h5, "λ", DT, ((elsize..., solution.nt + 1, solution.ni), (elsize..., -1, solution.ni)), chunk = (elsize..., 1, 1))
    for k in axes(solution.λ, 2)
        λ[elaxes..., 1, k] = solution.λ[0, k]
    end
end

init_solution(sio::SolutionHDF5, args...) = init_solution(hdf5(sio), args...)

function save_timeteps(h5::HDF5.File, sol::Solution, j1, j2, n1, n2)
    if size(h5["t"], 1) < j2
        HDF5.set_extent_dims(h5["t"], (j2,))
    end
    h5["t"][j1:j2] = timesteps(sol)[n1:n2]
end

save_timeteps(sio::SolutionHDF5, args...) = save_timeteps(hdf5(sio), args...)

function save_solution(h5::HDF5.File, solution::Union{SolutionODE{DT,TT,1},SolutionDAE{DT,TT,1}}, j1, j2, n1, n2) where {DT<:Number,TT}
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (j2,))
    end
    h5["q"][j1:j2] = solution.q[n1:n2]
end

function save_solution(h5::HDF5.File, solution::Union{SolutionODE{AT,TT,1},SolutionDAE{AT,TT,1}}, j1, j2, n1, n2) where {DT,AT<:Array{DT},TT}
    elaxes = axes(solution.q[begin])
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["q"][elaxes..., j] = solution.q[n]
    end
end

function save_solution(h5::HDF5.File, solution::Union{SolutionODE{DT,TT,2},SolutionDAE{DT,TT,2}}, j1, j2, n1, n2) where {DT<:Number,TT}
    if size(h5["q"], 1) < j2
        HDF5.set_extent_dims(h5["q"], (j2, size(h5["q"], 2)))
    end
    h5["q"][j1:j2, :] = solution.q[n1:n2, :]
end

function save_solution(h5::HDF5.File, solution::Union{SolutionODE{AT,TT,2},SolutionDAE{AT,TT,2}}, j1, j2, n1, n2) where {DT,AT<:Array{DT},TT}
    elaxes = axes(solution.q[begin, begin])
    if size(h5["q"])[end-1] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[begin:end-2]..., j2, size(h5["q"])[end]))
    end
    for k = 1:nsamples(solution.q)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["q"][elaxes..., j, k] = solution.q[n, k]
        end
    end
end

function save_solution(h5::HDF5.File, solution::Union{SolutionPODE{DT,TT,1},SolutionPDAE{DT,TT,1}}, j1, j2, n1, n2) where {DT<:Number,TT}
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (j2,))
    end
    if size(h5["p"])[end] < j2
        HDF5.set_extent_dims(h5["p"], (j2,))
    end
    h5["q"][j1:j2] = solution.q[n1:n2]
    h5["p"][j1:j2] = solution.p[n1:n2]
end

function save_solution(h5::HDF5.File, solution::Union{SolutionPODE{AT,TT,1},SolutionPDAE{AT,TT,1}}, j1, j2, n1, n2) where {DT,AT<:Array{DT},TT}
    elaxes = axes(solution.q[begin])
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["q"][elaxes..., j] = solution.q[n]
    end

    elaxes = axes(solution.p[begin])
    if size(h5["p"])[end] < j2
        HDF5.set_extent_dims(h5["p"], (size(h5["p"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["p"][elaxes..., j] = solution.p[n]
    end
end

function save_solution(h5::HDF5.File, solution::Union{SolutionPODE{DT,TT,2},SolutionPDAE{DT,TT,2}}, j1, j2, n1, n2) where {DT<:Number,TT}
    if size(h5["p"], 1) < j2
        HDF5.set_extent_dims(h5["p"], (j2, size(h5["p"], 2)))
    end
    if size(h5["p"], 1) < j2
        HDF5.set_extent_dims(h5["p"], (j2, size(h5["p"], 2)))
    end
    h5["q"][j1:j2, :] = solution.q[n1:n2, :]
    h5["p"][j1:j2, :] = solution.p[n1:n2, :]
end

function save_solution(h5::HDF5.File, solution::Union{SolutionPODE{AT,TT,2},SolutionPDAE{AT,TT,2}}, j1, j2, n1, n2) where {DT,AT<:Array{DT},TT}
    elaxes = axes(solution.q[begin, begin])
    if size(h5["q"])[end-1] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[begin:end-2]..., j2, size(h5["q"])[end]))
    end
    for k = 1:nsamples(solution.q)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["q"][elaxes..., j, k] = solution.q[n, k]
        end
    end

    elaxes = axes(solution.p[begin, begin])
    if size(h5["p"])[end-1] < j2
        HDF5.set_extent_dims(h5["p"], (size(h5["p"])[begin:end-2]..., j2, size(h5["p"])[end]))
    end
    for k = 1:nsamples(solution.p)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["p"][elaxes..., j, k] = solution.p[n, k]
        end
    end
end

save_solution(sio::SolutionHDF5, args...) = save_solution(hdf5(sio), args...)

function save_multiplier(h5::HDF5.File, solution::Union{SolutionDAE{DT,TT,1},SolutionPDAE{DT,TT,1}}, j1, j2, n1, n2) where {DT<:Number,TT}
    if size(h5["λ"])[end] < j2
        HDF5.set_extent_dims(h5["λ"], (j2,))
    end
    h5["λ"][j1:j2] = solution.λ[n1:n2]
end

function save_multiplier(h5::HDF5.File, solution::Union{SolutionDAE{AT,TT,1},SolutionPDAE{AT,TT,1}}, j1, j2, n1, n2) where {DT,AT<:Array{DT},TT}
    elaxes = axes(solution.λ[begin])
    if size(h5["λ"])[end] < j2
        HDF5.set_extent_dims(h5["λ"], (size(h5["λ"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["λ"][elaxes..., j] = solution.λ[n]
    end
end

function save_multiplier(h5::HDF5.File, solution::Union{SolutionDAE{DT,TT,2},SolutionPDAE{DT,TT,2}}, j1, j2, n1, n2) where {DT<:Number,TT}
    if size(h5["λ"], 1) < j2
        HDF5.set_extent_dims(h5["λ"], (j2, size(h5["λ"], 2)))
    end
    h5["λ"][j1:j2, :] = solution.λ[n1:n2, :]
end

function save_multiplier(h5::HDF5.File, solution::Union{SolutionDAE{AT,TT,2},SolutionPDAE{AT,TT,2}}, j1, j2, n1, n2) where {DT,AT<:Array{DT},TT}
    elaxes = axes(solution.λ[begin, begin])
    if size(h5["λ"])[end-1] < j2
        HDF5.set_extent_dims(h5["λ"], (size(h5["λ"])[begin:end-2]..., j2, size(h5["λ"])[end]))
    end
    for k = 1:nsamples(solution.λ)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["λ"][elaxes..., j, k] = solution.λ[n, k]
        end
    end
end

save_multiplier(::HDF5.File, ::Union{SolutionODE,SolutionPODE}, args...) = nothing
save_multiplier(sio::SolutionHDF5, args...) = save_multiplier(hdf5(sio), args...)


text:
using NeRCA
using Test

const ONLINEFILE = joinpath(@__DIR__, "data", "km3net_online.root")

hits = [Hit(1, 8, 100, 20, false),
        Hit(2, 9, 101, 21, true),
        Hit(3, 8, 112, 22, true),
        Hit(4, 8, 114, 23, false),
        Hit(5, 8, 134, 24, true),
        Hit(6, 8, 156, 25, false),
        Hit(7, 8, 133, 26, true),
        Hit(8, 8, 145, 26, false),
        Hit(9, 10, 900, 26, false),
        Hit(1, 11, 1000, 26, false),
        Hit(1, 11, 1001, 26, false),
        Hit(1, 11, 1002, 26, false)]

# triggered()
@testset "triggered()" begin
    thits = triggered(hits)
    @test 4 == length(thits)
    @test 9 == thits[1].dom_id
end

# nfoldhits()
@testset "nfoldhits()" begin
    twofoldhits = nfoldhits(hits, 10, 2)
    @test 4 == length(twofoldhits)
    threefoldhits = nfoldhits(hits, 15, 3)
    @test 3 == length(threefoldhits)
end

# domhits()
@testset "domhits()" begin
    dhits = domhits(hits)
    @test 7 == length(dhits[8])
    @test 20 == dhits[8][1].tot
    @test 1 == dhits[8][6].trigger_mask
end

# multiplicities
sorted_hits = sort(hits, by=h->h.t)
sort!(sorted_hits, by=h->h.dom_id)
mtps, mtp_ids = NeRCA.count_multiplicities(sorted_hits, 10)
@test (1, 1) == (mtps[1], mtp_ids[1])
@test (2, 2) == (mtps[2], mtp_ids[2])
@test (2, 2) == (mtps[3], mtp_ids[3])
@test (2, 3) == (mtps[4], mtp_ids[4])
@test (2, 3) == (mtps[5], mtp_ids[5])
@test (1, 4) == (mtps[6], mtp_ids[6])
@test (1, 5) == (mtps[7], mtp_ids[7])
@test (1, 6) == (mtps[8], mtp_ids[8])
@test (1, 7) == (mtps[9], mtp_ids[9])
@test (3, 8) == (mtps[10], mtp_ids[10])
@test (3, 8) == (mtps[11], mtp_ids[11])
@test (3, 8) == (mtps[12], mtp_ids[12])

# combine
event_shits = NeRCA.read_snapshot_hits(NeRCA.OnlineFile(ONLINEFILE))
event_thits = NeRCA.read_triggered_hits(NeRCA.OnlineFile(ONLINEFILE))
for (shits, thits) in zip(event_shits, event_thits)
    local hits = NeRCA.combine(shits, thits)
    @test length(hits) == length(shits)
    @test length(NeRCA.triggered(hits)) == length(thits)
end


text:
module RulesWishartOutTest

using Test
using ReactiveMP
using Random
using Distributions

import ReactiveMP: @test_rules

@testset "rules:Wishart:out" begin

    @testset "Belief Propagation: (m_ν::PointMass, m_S::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Wishart(:out, Marginalisation) [
            (input = (m_ν = PointMass(2.0), m_S = PointMass([ 1.0 0.0; 0.0 1.0 ])), output = Wishart(2.0, [ 1.0 0.0; 0.0 1.0 ])),
            (input = (m_ν = PointMass(3.0), m_S = PointMass([ 10.0 -1.0; -1.0 3.0 ])), output = Wishart(3.0, [ 10.0 -1.0; -1.0 3.0 ])),
            (input = (m_ν = PointMass(4.0), m_S = PointMass([ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])), output = Wishart(4.0, [ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])),
        ]

    end

    @testset "Variational: (m_ν::PointMass, q_S::Any)" begin

        @test_rules [ with_float_conversions = true ] Wishart(:out, Marginalisation) [
            (input = (m_ν = PointMass(2.0), q_S = PointMass([ 1.0 0.0; 0.0 1.0 ])), output = Wishart(2.0, [ 1.0 0.0; 0.0 1.0 ])),
            (input = (m_ν = PointMass(3.0), q_S = PointMass([ 10.0 -1.0; -1.0 3.0 ])), output = Wishart(3.0, [ 10.0 -1.0; -1.0 3.0 ])),
            (input = (m_ν = PointMass(4.0), q_S = PointMass([ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])), output = Wishart(4.0, [ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])),
        ]

    end

    @testset "Variational: (q_ν::Any, m_S::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Wishart(:out, Marginalisation) [
            (input = (q_ν = PointMass(2.0), m_S = PointMass([ 1.0 0.0; 0.0 1.0 ])), output = Wishart(2.0, [ 1.0 0.0; 0.0 1.0 ])),
            (input = (q_ν = PointMass(3.0), m_S = PointMass([ 10.0 -1.0; -1.0 3.0 ])), output = Wishart(3.0, [ 10.0 -1.0; -1.0 3.0 ])),
            (input = (q_ν = PointMass(4.0), m_S = PointMass([ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])), output = Wishart(4.0, [ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])),
        ]

    end

    @testset "Variational: (m_ν::Any, m_S::Any)" begin

        @test_rules [ with_float_conversions = true ] Wishart(:out, Marginalisation) [
            (input = (q_ν = PointMass(2.0), q_S = PointMass([ 1.0 0.0; 0.0 1.0 ])), output = Wishart(2.0, [ 1.0 0.0; 0.0 1.0 ])),
            (input = (q_ν = PointMass(3.0), q_S = PointMass([ 10.0 -1.0; -1.0 3.0 ])), output = Wishart(3.0, [ 10.0 -1.0; -1.0 3.0 ])),
            (input = (q_ν = PointMass(4.0), q_S = PointMass([ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])), output = Wishart(4.0, [ 1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0 ])),
        ]

    end

end

end


text:
using BenchmarkTools

module A

struct Triangle
    side1
    side2
    side3
    angle1
    angle2
    angle3
end

end # end of module A

module B

struct Triangle
    side1   :: Float64
    side2   :: Float64
    side3   :: Float64
    angle1  :: Float64
    angle2  :: Float64
    angle3  :: Float64
end

end # end of module B

module C

struct Triangle{T}
    side1   ::  T
    side2   ::  T
    side3   ::  T
    angle1  ::  T
    angle2  ::  T
    angle3  ::  T
end

end # end of module C

using .A
using .B
using .C

function main()

    a::Float64 = 5.0

    # t = A.Triangle(a, a, a, a, a, a)
    # t = B.Triangle(a, a, a, a, a, a)
    t = C.Triangle{Float64}(a, a, a, a, a, a)

end

@benchmark main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
using Pkg
Pkg.add("IJulia")
Pkg.add("Plots")

max_stars_repo_name: Montereytony/jupyterDesktop

max_issues_repo_name: Montereytony/jupyterDesktop

max_forks_repo_name: Montereytony/jupyterDesktop


text:
using MeshCat, MeshCatMechanisms, RigidBodyDynamics, GeometryTypes, CoordinateTransformations
using Zygote
using LinearAlgebra: I
using BatchedRoutines

"""
Visualize a trajectory in the browser

$(SIGNATURES)

**Arguments**:

- `qs :: Vector{Vector{Float64}}` Position vectors of the trajectory
- `ts` Timesteps of the trajectory
- `mechanism` Robot model as per `RigidBodyDynamics`
- `urdf` Path to URDF of the model
"""
function viz_traj(qs, ts, mechanism, urdf)
    vis = Visualizer()
    open(vis)
    mvis = MechanismVisualizer(mechanism, URDFVisuals(urdf), vis)
    settransform!(vis["/Cameras/default"], compose(Translation(0.0, 1.0, 1.0), LinearMap(RotZ(pi/2))));
    setanimation!(mvis, ts, qs)
end

"""
Matrix-matrix multiplication _batched_ and adjoint defined

$(SIGNATURES)
"""
function batchmul(A::AbstractArray{T,3},B::AbstractArray{V,3}) where {T, V}
    # XXX: To deal with Dual numbers from ForwardDiff
    if T ∉ (Float64, Float32) || V ∉ (Float64, Float32)
        return batchmul_slow(A, B)
    end
    batched_gemm('N','N', A, B)
end

function batchmul_slow(A, B)
    batchmm = [A[:,:,i] * B[:, :, i] for i in 1:size(A, 3)]
    C = reduce((xs,x)->cat(xs, x; dims=3),
               [reshape(batchmm[i], size(batchmm[i])..., 1) for i in 1:size(A)[end]])
    return C
end

Zygote.@adjoint function batchmul(A, B)
    batchmul(A, B), Δ -> (batched_gemm('N' ,'T', Δ , B), batched_gemm( 'T' ,'N', A, Δ))
end

eye(T, p) = Matrix{T}(I, p, p)
Zygote.@nograd eye


text:
module MDFlattenTests

if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using Base.Markdown
using Documenter.Utilities.MDFlatten

@testset "MDFlatten" begin
    @test mdflatten(Markdown.Paragraph("...")) == "..."
    @test mdflatten(Markdown.Header{1}("...")) == "..."

    # a simple test for blocks in top-level (each gets two newline appended to it)
    @test mdflatten(Markdown.parse("# Test\nTest")) == "Test\n\nTest\n\n"
    block_md = Markdown.parse("""
    # MDFlatten test

    ^^^ Ignoring extra whitespace.

    ```markdown
    code
    is forwarded as **is**
    ```
    """)
    block_text = """
    MDFlatten test

    ^^^ Ignoring extra whitespace.

    code
    is forwarded as **is**

    """
    @test mdflatten(block_md) == block_text

    # blocks
    @test mdflatten(Markdown.parse("> Test\n> Test\n\n> Test")) == "Test Test\n\nTest\n\n"
    @test mdflatten(Markdown.parse("HRs\n\n---\n\nto whitespace")) == "HRs\n\n\n\nto whitespace\n\n"
    @test mdflatten(Markdown.parse("HRs\n\n---\n\nto whitespace")) == "HRs\n\n\n\nto whitespace\n\n"
    @test mdflatten(Markdown.parse("HRs\n\n---\n\nto whitespace")) == "HRs\n\n\n\nto whitespace\n\n"

    # test some inline blocks
    @test mdflatten(Markdown.parse("`code` *em* normal **strong**")) == "code em normal strong\n\n"
    @test mdflatten(Markdown.parse("[link text *parsed*](link/itself/ignored)")) == "link text parsed\n\n"
    @test mdflatten(Markdown.parse("- a\n- b\n- c")) == "a\nb\nc\n\n"
    @test mdflatten(Markdown.parse("A | B\n---|---\naa|bb\ncc | dd")) == "A B\naa bb\ncc dd\n\n"

    # Math
    @test mdflatten(Markdown.parse("\$e=mc^2\$")) == "e=mc^2\n\n"
    if VERSION > v"0.5-"
        # backticks and blocks for math only in 0.5, i.e. these fail on 0.4
        @test mdflatten(Markdown.parse("``e=mc^2``")) == "e=mc^2\n\n"
        @test mdflatten(Markdown.parse("```math\n\\(m+n)(m-n)\nx=3\\sin(x)\n```")) == "(m+n)(m-n)\nx=3sin(x)\n\n"
    end

    # symbols in markdown
    @test mdflatten(Markdown.parse("A \$B C")) == "A B C\n\n"

    # linebreaks
    @test mdflatten(Markdown.parse("A\\\nB")) == "A\nB\n\n"

    # Only available on Julia 0.5.
    if isdefined(Base.Markdown, :Footnote)
        @test mdflatten(Markdown.parse("[^name]")) == "[name]\n\n"
        @test mdflatten(Markdown.parse("[^name]:**Strong** text.")) == "[name]: Strong text.\n\n"
    end

    if isdefined(Base.Markdown, :Admonition)
        @test mdflatten(Markdown.parse("!!! note \"Admonition Title\"\n    Test")) == "note: Admonition Title\nTest\n\n"
    end
end

end


text:
## Plot recipe
## define a hueristic to work around asymptotes
## just sort of succesful
@recipe function f(pq::AbstractRationalFunction{T}, a=nothing, b=nothing) where {T}

    xlims = get(plotattributes, :xlims, (nothing,nothing))
    ylims = get(plotattributes, :ylims, (nothing, nothing))
    rational_function_trim(pq, a, b, xlims, ylims)

end

isapproxreal(x::Real) = true
isapproxreal(x::Complex{T}) where {T} = imag(x) <= sqrt(eps(real(T)))
function toobig(pq)
    x -> begin
        y = pq(x)
        isinf(y) && return true
        isnan(y) && return true
        abs(y) > 1e8 && return true
        return false
    end
end

function rational_function_trim(pq, a, b, xlims, ylims)

    p,q = lowest_terms(//(pq...), method=:numerical)
    dpq = derivative(p//q)
    p′,q′ = lowest_terms(dpq)

    λs = Multroot.multroot(q).values
    λs = isempty(λs) ? λs : real.(filter(isapproxreal, λs))

    cps = Multroot.multroot(p′).values
    cps = isempty(cps) ? cps : real.(filter(isapproxreal, cps))
    cps = isempty(cps) ? cps : filter(!toobig(pq), cps)

    a = a == nothing ? xlims[1] : a
    b = b == nothing ? xlims[2] : b

    if a==nothing && b==nothing
        u= poly_interval(p)
        v= poly_interval(q)
        a,b = min(first(u), first(v)), max(last(u), last(v))

        if !isempty(λs)
            a,b = min(a, real(minimum(λs))), max(b, real(maximum(λs)))
        end
        if !isempty(cps)
            a,b = min(a, real(minimum(cps))), max(b, real(maximum(cps)))
        end
        a *= (a > 0 ? 1/1.5 : 1.25)
        b *= (b < 0 ? 1/1.5 : 1.25)
    end

    n = 601
    xs = range(a,stop=b, length=n)
    ys = pq.(xs)
    Mcps = isempty(cps) ? 5 : 3*maximum(abs, pq.(cps))
    M = max(5, Mcps, 1.25*maximum(abs, pq.((a,b))))

    lo = ylims[1] == nothing ? -M : ylims[1]
    hi = ylims[2] == nothing ?  M : ylims[2]
    ys′ = [lo <= yᵢ <= hi ? yᵢ : NaN for yᵢ ∈ ys]
    xs, ys′

end


text:
# Copyright 2013 Peter Chronz
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

# OFP_TYPE
# Immutable messages
const OFPT_HELLO = 0x00 # Symmetric message
const OFPT_ERROR = 0x01 # Symmetric message
const OFPT_ECHO_REQUEST = 0x02 # Symmetric message
const OFPT_ECHO_REPLY = 0x03 # Symmetric message
const OFPT_VENDOR = 0x04 # Symmetric message
# Switch configuration messages
const OFPT_FEATURES_REQUEST = 0x05 # Controller/switch message
const OFPT_FEATURES_REPLY = 0x06 # Controller/switch message
const OFPT_GET_CONFIG_REQUEST = 0x07 # Controller/switch message
const OFPT_GET_CONFIG_REPLY = 0x08 # Controller/switch message
const OFPT_SET_CONFIG = 0x09 # Controller/switch message
# Asynchronous messages
const OFPT_PACKET_IN = 0x0a # Async message
const OFPT_FLOW_REMOVED = 0x0b # Async message
const OFPT_PORT_STATUS = 0x0c # Async message
# Controller command messages
const OFPT_PACKET_OUT = 0x0d # Controller/switch message
const OFPT_FLOW_MOD = 0x0e # Controller/switch message
const OFPT_PORT_MOD = 0x0f # Controller/switch message
# Statistics messages
const OFPT_STATS_REQUEST = 0x10 # Controller/switch message
const OFPT_STATS_REPLY = 0x11 # Controller/switch message
# Barrier messages
const OFPT_BARRIER_REQUEST = 0x12 # Controller/switch message
const OFPT_BARRIER_REPLY = 0x13 # Controller/switch message
# Queue configuration messages
const OFPT_QUEUE_GET_CONFIG_REQUEST = 0x14 # Controller/switch message
const OFPT_QUEUE_GET_CONFIG_REPLY = 0x15 # Controller/switch message

# Port-related constants
# OFP_PORT_CONFIG Flags to indicate behavior of the physical port. These flags
    # are used in ofp_py_port to describe the current configuration. They are used
    # in the ofp_port_mod message to configure the port's behavior.
const OFPPC_PORT_DOWN = 1 << 0 # Port is administratively down
const OFPPC_NO_STP = 1 << 1 # Disable 802.1D spanning tree on port.
const OFPPC_NO_RECV = 1 << 2 # Drop all packets except 802.1D spanning tree packets.
const OFPPC_NO_RECV_STP = 1 << 3 # Drop received 802.1D STP packets.
const OFPPC_NO_FLOOD = 1 << 4 # Do not include this port when flooding.
const OFPPC_NO_FWD = 1 << 5 # Drop packets forwarded to port.
const OFPPC_NO_PACKET_IN = 1 << 6 # Do not send packet-in msgs fort port.
# OFP_PORT_STATE Current state of the physical port. These are not configurable
    # from the controller.
const OFPPS_LINK_DOWN = 1 << 0 # No physical link present.
# The OFPPS_STP_* bits have no effect on switch operation. The controller must
    # adjust OFPPC_NO_RECV, OFPPC_NO_FWD, and OFPPC_NO_PACKET_IN appropriately to
    # fully implement an 802.1D spanning tree.
const OFPPS_STP_LISTEN = 0 << 8 # Not learning or relaying frames.
const OFPPS_STP_LEARN = 1 << 8 # Learning but not relaying frames.
const OFPPS_STP_FORWARD = 2 << 8 # Learning and relaying frames.
const OFPPS_STP_BLOCK = 3 << 8 # Not part of spanning tree.
const OFPPS_STP_MASK = 3 << 8 # Bit mask for OFPPS_STP_* values.
# OFP_PORT Port numbering. Physical ports are numbered starting from 1.
# Maximum number of physical switch ports.
const OFPP_MAX = 0xff00
# Fake output "ports".
const OFPP_IN_PORT = 0xfff8 # Send the packet out the input port. This virtual
    # port must be explicitly used in order to send back out the input port.
const OFPP_TABLE = 0xfff9 # Perform actions in flow table. NB: This can only be
    # the destination port for packet-out messages.
const OFPP_NORMAL = 0xfffa # Process with normal L2/L3 switching.
const OFPP_FLOOD = 0xfffb # All physical ports except input port and those
# disabled by STP.
const OFPP_ALL = 0xfffc # All physical ports except input port.
const OFPP_CONTROLLER = 0xfffd # Send to controller.
const OFPP_LOCAL = 0xfffe # Local openflow "port".
const OFPP_NONE = 0xffff # Not associated with a physical port.
# OFP_PORT_FEATURES Features of physical ports available in a datapath.
const OFPPF_10MB_HD = 1 << 0 # 10 Mb half-duplex rate support.
const OFPPF_10MB_FD = 1 << 1 # 10 Mb full-duplex rate suport.
const OFPPF_100MB_HD = 1 << 2 # 100 Mb half-duplex rate support.
const OFPPF_100MB_FD = 1 << 3 # 100 Mb full-duplex rate support.
const OFPPF_1GB_HD = 1 << 4 # 1 Gb half-duplex rate support.
const OFPPF_1GB_FD = 1 << 5 # 1 Gb full-duplex rate support.
const OFPPF_10GB_FD = 1 << 6 # 10Gb full-duplex rate support.
const OFPPF_COPPER = 1 << 7 # Copper medium.
const OFPPF_FIBER = 1 << 8 # Fiber medium.
const OFPPF_AUTONEG = 1 << 9 # Auto-negotiation.
const OFPPF_PAUSE = 1 << 10 # Pause.
const OFPPF_PAUSE_ASYM = 1 << 11 # Asymmetric pause.

# OFP_CAPABILITIES Capabilities supported by datapath.
const OFPC_FLOW_STATS = 1 << 0 # Flow statistics.
const OFPC_TABLE_STATS = 1 << 1 # Table statistics.
const OFPC_PORT_STATS = 1 << 2 # Port statistics.
const OFPC_STP = 1 << 3 # 802.1d spanning tree.
const OFPC_RESERVED = 1 << 4 # Reserved, must be zero.
const OFPC_IP_REASM = 1 << 5 # Can reassemble IP fragments.
const OFPC_QUEUE_STATS = 1 << 6 # Queue statistics.
const OFPC_ARP_MATCH_IP = 1 << 7 # Match IP addressed in ARP pkts.

# Switch configuration flags
# OFP_CONFIG_FLAGS
# Handling of IP fragments.
const OFPC_FRAG_NORMAL = 0 # No special handling for fragments.
const OFPC_FRAG_DROP = 1 # Drop fragments.
const OFPC_FRAG_REASM = 2 # Reassemble (only if OFPC_IP_REASM set).
const OFPC_FRAG_MASK = 3

# Packet in reasons
# OFP_PACKET_IN_REASON
const OFPR_NO_MATCH =  0x00 # No matching flow.
const OFPR_ACTION = 0x00 # Action explicitly output to controller.

# Port status
# OFP_PORT_REASON
# What changed about the physical port.
const OFPPR_ADD = 0x00 # The port was added.
const OFPPR_DELETE = 0x01 # The port was removed.
const OFPPR_MODIFY = 0x02 # Some attribute of the port has changed.

# Flow wildcards
# OFP_FLOW_WILDCARDS
const OFPFW_IN_PORT = 1 << 0 # Switch input port.
const OFPFW_DL_VLAN = 1 << 1 # VLAN id.
const OFPFW_DL_SRC = 1 << 2 # Ethernet source address.
const OFPFW_DL_DST = 1 << 3 # Ethernet destination address.
const OFPFW_DL_TYPE = 1 << 4 # Ethernet frame type.
const OFPFW_NW_PROTO = 1 << 5 # IP protocol.
const OFPFW_TP_SRC = 1 << 6 # TCP/UDP source port.
const OFPFW_TP_DST = 1 << 7 # TCP/UDP destination port.
# IP source address wildcard bit count. 0 is exact match, 1 ignores the LSB, 2
# ignores the 2 least significant bits, ..., 32 and higher wildcard the entire
# field. This is the *opposite* of the usual convention where e.g. /24 indicates
# that 8 bits (not 24 bits) are wildcarded.
const OFPFW_NW_SRC_SHIFT = 8
const OFPFW_NW_SRC_BITS = 6
const OFPFW_NW_SRC_MASK = ((1 << OFPFW_NW_SRC_BITS) - 1) << OFPFW_NW_SRC_SHIFT
const OFPFW_NW_SRC_ALL = 32 << OFPFW_NW_SRC_SHIFT
# IP destination address wildcard bit count. Same format as source.
const OFPFW_NW_DST_SHIFT = 14
const OFPFW_NW_DST_BITS = 6
const OFPFW_NW_DST_MASK = ((1 << OFPFW_NW_DST_BITS) - 1) << OFPFW_NW_DST_SHIFT
const OFPFW_NW_DST_ALL = 32 << OFPFW_NW_DST_SHIFT
const OFPFW_DL_VLAN_PCP = 1 << 20 # VLAN priority.
const OFPFW_NW_TOS = 1 << 21 # IP ToS (DSCP field, 6 bits).
# Wildcards all fields.
const OFPFW_ALL = (1 << 22) - 1

# OFP_ACTION_TYPE
const OFPAT_OUTPUT = 0x00 # Output to switch port.
const OFPAT_SET_VLAN_VID = 0x01 # Set the 802.1q VLAN id.
const OFPAT_SET_VLAN_PCP = 0x02 # Set the 802.1q priority.
const OFPAT_STRIP_VLAN = 0x03 # Strip the 802.1q header.
const OFPAT_SET_DL_SRC = 0x04 # Ethernet source address.
const OFPAT_SET_DL_DST = 0x05 # Ethernet destination address.
const OFPAT_SET_NW_SRC = 0x06 # IP source address.
const OFPAT_SET_NW_DST = 0x07 # IP destination address.
const OFPAT_SET_NW_TOS = 0x08 # IP ToS (DSCP field, 6 bits).
const OFPAT_SET_TP_SRC = 0x09 # TCP/UDP source port.
const OFPAT_SET_TP_DST = 0x0a # TCP/UDP destination port.
const OFPAT_ENQUEUE = 0x0b # Output to queue.
const OFPAT_VENDOR = 0xffff

# OFP_FLOW_MOD_COMMAND
const OFPFC_ADD = 0x00 # New flow.
const OFPFC_MODIFY = 0x01 # Modify all matching flows.
const OFPFC_MODIFY_STRICT = 0x02 # Modify entry strictly matching wildcards.
const OFPFC_DELETE = 0x03 # Delete all matching flows.
const OFPFC_DELETE_STRICT = 0x04 # Strictly match wildcards and priority.

# OFP_FLOW_MOD_FLAGS
const OFPFF_SEND_FLOW_REM = 1 << 0 # Send flow removed message when flow expires
                                    # or is deleted
const OFPFF_CHECK_OVERLAP = 1 << 1 # Check for overlapping entries first.
const OFPFF_EMERG = 1 << 2 # Remark this is for emergency.

# Error messages.
# OFP_ERROR_TYPE
const OFPET_HELLO_FAILED = 0x00 # Hello protocol failed.
const OFPET_BAD_REQUEST = 0x01 # Request was not understood.
const OFPET_BAD_ACTION = 0x02 # Error in action description.
const OFPET_FLOW_MOD_FAILED = 0x03 # Problem modifying flow entry.
const OFPET_PORT_MOD_FAILED = 0x04 # Port mod request failed.
const OFPET_QUEUE_OP_FAILED = 0x05 # Queue operation failed.
# OFP_BAD_REQUEST_CODE
const OFPBRC_BAD_VERSION = 0x00 # ofp_header.version not supported.
const OFPBRC_BAD_TYPE = 0x01 # ofp_header.type not supported.
const OFPBRC_BAD_STAT = 0x02 # ofp_stats_request.type not supported.
const OFPBRC_BAD_VENDOR = 0x03 # Vendor not supported (in ofp_vendor_header or
    # ofp_stats_request or ofp_stats_reply).
const OFPBRC_BAD_SUBTYPE = 0x04 # Vendor subtype not supported.
const OFPBRC_EPERM = 0x05 # Permissions error.
const OFPBRC_BAD_LEN = 0x06 # Wrong request length for type.
const OFPBRC_BUFFER_EMPTY = 0x07 # Specified buffer has already been used.
const OFPBRC_BUFFER_UNKNOWN = 0x08 # Specified buffer does not exist.
# OFP_BAD_ACTION_CODE
# ofp_error_msg 'code' values for OFPET_BAD_ACTION. 'data' contains at least the
# first 64 bytes of the failed request.
const OFPBAC_BAD_TYPE = 0x00 # Unknown action type.
const OFPBAC_BAD_LEN = 0x01 # Length problem in actions.
const OFPBAC_BAD_VENDOR = 0x02 # Unkown vendor id specified.
const OFPBAC_BAD_VENDOR_TYPE = 0x03 # Unkown action type for vendor id.
const OFPBAC_BAD_OUT_PORT = 0x04 # Problem validating output action.
const OFPBAC_BAD_ARGUMENT = 0x05 # Bad action argument.
const OFPBAC_EPERM = 0x06 # Permissions error.
const OFPBAC_TOO_MANY = 0x07 # Can't handle this many actions.
const OFPBAC_BAD_QUEUE = 0x08 # Problem validating output queue.
# OFP_FLOW_MOD_FAILED_CODE
# ofp_error_msg 'code' values for OFPET_FLOW_MOD_FAILED. 'data' contains at
# least the first 64 bytes of the failed request.
const OFPFMFC_ALL_TABLES_FULL =  # Flow not added because of full tables.
const OFPFMFC_OVERLAP =  # Attempted to add overlapping flow with CHECK_OVERLAP
    # flag set.
const OFPFMFC_EPERM =  # Permissions error.
const OFPFMFC_BAD_EMERG_TIMEOUT =  # Flow not added because of non-zero
    # idle/hard timeout.
const OFPFMFC_BAD_COMMAND =  # Unknown command.
const OFPFMFC_UNSUPPORTED =  # Unsupported action list - cannot process in the
    # order specified.
# OFP_PORT_MOD_FAILED_CODE
# ofp_error_msg 'code' values for OFPET_PORT_MOD_FAILED> 'data' contains at
# least the first 64 bytes of the failed request.
const OFPMFC_BAD_PORT = 0x00 # Specified port does not exist.
const OFPMFC_BAD_HW_ADDR = 0x01 # Specified hardware address is wrong.
# OFP_QUEUE_OP_FAILED_CODE
# ofp_error msg 'code' values for OFPET_QUEUE_OP_FAILED. 'data' contains at
# least the first 64 bytes of the failed request.
const OFPQOFC_BAD_PORT = 0x00 # Invalid port (or port does not exist).
const OFPQOFC_BAD_QUEUE = 0x01 # Queue does not exist.
const OFPQOFC_EPERM = 0x02 # Permissions error.
# OfpQueueProperties
const OFPQT_NONE = 0 # No property defined for a queue (default).
const OFPQT_MIN_RATE = 1 # Minimum datarate guaranteed.
                           # Other rates should be added here (i.e. max rate,
                           # precedence, etc).
# OFP_STATS_TYPES
const OFPST_DESC = 0x00 # Description of this OpenFlow switch. The request body
    # is empty. The reply body is struct ofp_desc-stats.
const OFPST_FLOW = 0x01 # Individual flow statistics. The request body is struct
    # ofp_flow_stats_request. The reply body is an array of structu ofp_flow_stats.
const OFPST_AGGREGATE = 0x02 # Aggregate flow statistics. The request body is
    # struct ofp_aggregate_stats_request. The reply body is struct
    # ofp_aggregate_reply.
const OFPST_TABLE = 0x03 # Flow table statistics. The request body is empty. The
    # reply body is an array of struct ofp_table_stats.
const OFPST_PORT = 0x04 # Physical port statistics. The request body is struct
    # ofp_port_stats_request. The reply body is an array of struct ofp_port_stats.
const OFPST_QUEUE = 0x05 # Queue statistics for a port. The request body defines
    # the port. The reply body is an array of struct ofp_queue_stats.
const OFPST_VENDOR = 0xffff # Vendor extension. The request and reply bodies
    # begin witha 32-bit vendor ID, which takes the same form as in "struct
    # ofp_vendor_header". The request and reply bodies are otherwise
    # vendor-defined.
# OFP_FLOW_REMOVED_REASON
const OFPRR_IDLE_TIMEOUT = 0x00 # Flow idle time exceeded idle_timeout.
const OFPRR_HARD_TIMEOUT = 0x01 # Time exceeded hard_timeout.
const OFPRR_DELTE = 0x02 # Evicted by a DELETE flow mod.
# Other used constants
const OFP_MAX_ETH_ALEN = 6
const OFP_MAX_PORT_NAME_LEN = 16
const DESC_STR_LEN = 256
const SERIAL_NUM_LEN = 32
const OFP_MAX_TABLE_NAME_LEN = 32

max_stars_repo_name: pchronz/juliaflow

max_issues_repo_name: pchronz/juliaflow

max_forks_repo_name: pchronz/juliaflow


text:
using DiffEqBase, DiffEqJump, Test, StableRNGs

rng = StableRNG(12345)

# test for https://github.com/SciML/DiffEqJump.jl/issues/177
p   = [1., 2., 50.]
ns  = [[1 => 1], [1 => -1], [1 => 1]]
rs  = [[1 => 1], [1 => 1], Pair{Int64,Int64}[]]
maj = MassActionJump(p, rs, ns)
u0    = [5]
tspan = (0., 2e6)
dt    = tspan[2] / 1000
dprob = DiscreteProblem(u0, tspan, p)
jprob = JumpProblem(dprob, Direct(), maj, save_positions=(false,false), rng=rng)
sol   = solve(jprob, SSAStepper(), saveat=tspan[1]:dt:tspan[2])
@test length(unique(sol[(end-10):end][:])) > 1


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DelayEmbeddings
using DelimitedFiles
using StatsBase
using HypothesisTests

include("../../../src/pecuzal_method.jl")
include("../../../src/data_analysis_functions.jl")

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we look at real data of a chaotic electrochemical oscillator (Ni) from
# the Kiss group with two different resistances R₁ = 1,5kΩ, R₂ = 2,5kΩ leading to
# "more" chaoticity.

RQA_tde1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde1.csv")
RQA_tde2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde2.csv")
RQA_GA1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA1.csv")
RQA_GA2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA2.csv")
RQA_mdop1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop1.csv")
RQA_mdop2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop2.csv")
RQA_pec1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec1.csv")
RQA_pec2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec2.csv")

dim_tde1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde1.csv")
dim_tde2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde2.csv")
L_tde1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_tde1.csv")
L_tde2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_tde2.csv")
dim_GA1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA1.csv")
dim_GA2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA2.csv")
L_GA1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_GA1.csv")
L_GA2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_GA2.csv")
dim_mdop1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop1.csv")
dim_mdop2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop2.csv")
L_mdop1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_mdop1.csv")
L_mdop2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_mdop2.csv")
dim_pec1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec1.csv")
dim_pec2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec2.csv")
L_pec1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_pec1.csv")
L_pec2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/L_pec2.csv")

using Plots

p1 = histogram(RQA_tde1[1,:],title = "TDE 1 ENTR")
p2 = histogram(RQA_GA1[1,:],title = "GA 1 ENTR")
p3 = histogram(RQA_mdop1[1,:],title = "MDOP 1 ENTR")
p4 = histogram(RQA_pec1[1,:],title = "PEC 1 ENTR")
p5 = histogram(RQA_tde2[1,:],title = "TDE 2 ENTR")
p6 = histogram(RQA_GA2[1,:],title = "GA 2 ENTR")
p7 = histogram(RQA_mdop2[1,:],title = "MDOP 2 ENTR")
p8 = histogram(RQA_pec2[1,:],title = "PEC 2 ENTR")
plot(p1,p2,p3,p4,layout = (2, 2), legend = false)
plot(p5,p6,p7,p8,layout = (2, 2), legend = false)

p1 = histogram(RQA_tde1[2,:],title = "TDE 1 TRANS")
p2 = histogram(RQA_GA1[2,:],title = "GA 1 TRANS")
p3 = histogram(RQA_mdop1[2,:],title = "MDOP 1 TRANS")
p4 = histogram(RQA_pec1[2,:],title = "PEC 1 TRANS")
p5 = histogram(RQA_tde2[2,:],title = "TDE 2 TRANS")
p6 = histogram(RQA_GA2[2,:],title = "GA 2 TRANS")
p7 = histogram(RQA_mdop2[2,:],title = "MDOP 2 TRANS")
p8 = histogram(RQA_pec2[2,:],title = "PEC 2 TRANS")
plot(p1,p2,p3,p4,layout = (2, 2), legend = false)
plot(p5,p6,p7,p8,layout = (2, 2), legend = false)

med_TDE = zeros(2,4)
med_GA = zeros(2,4)
med_mdop = zeros(2,4)
med_pec = zeros(2,4)
std_TDE = zeros(2,4)
std_GA = zeros(2,4)
std_mdop = zeros(2,4)
std_pec = zeros(2,4)

Ltde1 = mean(L_tde1)
Lstd_tde1 = std(L_tde1)
Ltde2 = mean(L_tde2)
Lstd_tde2 = std(L_tde2)
dim_tde1 = mean(dim_tde1)
dim_tde2 = mean(dim_tde2)
dimstd_tde1 = std(dim_tde1)
dimstd_tde2 = std(dim_tde2)

LGA1 = mean(L_GA1)
Lstd_GA1 = std(L_GA1)
LGA2 = mean(L_GA2)
Lstd_GA2 = std(L_GA2)
dim_GA1 = mean(dim_GA1)
dim_GA2 = mean(dim_GA2)
dimstd_GA1 = std(dim_GA1)
dimstd_GA2 = std(dim_GA2)

Lmdop1 = mean(L_mdop1)
Lstd_mdop1 = std(L_mdop1)
Lmdop2 = mean(L_mdop2)
Lstd_mdop2 = std(L_mdop2)
dim_mdop1 = mean(dim_mdop1)
dim_mdop2 = mean(dim_mdop2)
dimstd_mdop1 = std(dim_mdop1)
dimstd_mdop2 = std(dim_mdop2)

Lpec1 = mean(L_pec1)
Lstd_pec1 = std(L_pec1)
Lpec2 = mean(L_pec2)
Lstd_pec2 = std(L_pec2)
dim_pec1 = mean(dim_pec1)
dim_pec2 = mean(dim_pec2)
dimstd_pec1 = std(dim_pec1)
dimstd_pec2 = std(dim_pec2)

for i = 1:4
    med_TDE[1,i] = mean(RQA_tde1[i,:])
    med_TDE[2,i] = mean(RQA_tde2[i,:])
    med_GA[1,i] = mean(RQA_GA1[i,:])
    med_GA[2,i] = mean(RQA_GA2[i,:])
    med_mdop[1,i] = mean(RQA_mdop1[i,:])
    med_mdop[2,i] = mean(RQA_mdop2[i,:])
    med_pec[1,i] = mean(RQA_pec1[i,:])
    med_pec[2,i] = mean(RQA_pec2[i,:])

    std_TDE[1,i] = std(RQA_tde1[i,:])
    std_TDE[2,i] = std(RQA_tde2[i,:])
    std_GA[1,i] = std(RQA_GA1[i,:])
    std_GA[2,i] = std(RQA_GA2[i,:])
    std_mdop[1,i] = std(RQA_mdop1[i,:])
    std_mdop[2,i] = std(RQA_mdop2[i,:])
    std_pec[1,i] = std(RQA_pec1[i,:])
    std_pec[2,i] = std(RQA_pec2[i,:])
end

RQA_measure = 4
# 1: ENTR
# 2: LAM
# 3: RTE
# 4: TRANS

wtest_tde = UnequalVarianceTTest(RQA_tde1[RQA_measure,:], RQA_tde2[RQA_measure,:])
wtest_GA = UnequalVarianceTTest(RQA_GA1[RQA_measure,:], RQA_GA2[RQA_measure,:])
wtest_mdop = UnequalVarianceTTest(RQA_mdop1[RQA_measure,:], RQA_mdop2[RQA_measure,:])
wtest_pec = UnequalVarianceTTest(RQA_pec1[RQA_measure,:], RQA_pec2[RQA_measure,:])

pvalue(wtest_tde)
pvalue(wtest_GA)
pvalue(wtest_mdop)
pvalue(wtest_pec)

wtest_tde = UnequalVarianceTTest(vec(L_tde1), vec(L_tde2))
wtest_GA = UnequalVarianceTTest(vec(L_GA1), vec(L_GA2))
wtest_mdop = UnequalVarianceTTest(vec(L_mdop1), vec(L_mdop2))
wtest_pec = UnequalVarianceTTest(vec(L_pec1), vec(L_pec2))

pvalue(wtest_tde)
pvalue(wtest_GA)
pvalue(wtest_mdop)
pvalue(wtest_pec)

max_stars_repo_path: scripts/Experimental data/Standard electrochemical chaos/analyze_standard_ecoscillator.jl


text:
import Base: size, length, reshape, hcat, vcat, fill

# Let the user get the `size` and `length` of `Node`s.
Base.size(x::Node, dims...) = size(unbox(x), dims...)
Base.length(x::Node) = length(unbox(x))

# Sensitivity for the first argument of `reshape`.
@explicit_intercepts reshape Tuple{∇Array, Vararg{Int}} [true, false]
@explicit_intercepts reshape Tuple{∇Array, Tuple{Vararg{Int}}} [true, false]
∇(::typeof(reshape), ::Type{Arg{1}}, _, y, ȳ, A::∇Array, args...) =
    reshape(ȳ, size(A)...)

@union_intercepts hcat Tuple{Vararg{∇Array}} Tuple{Vararg{AbstractArray}}
function Nabla.∇(
    ::typeof(hcat),
    ::Type{Arg{i}},
    _,
    y,
    ȳ,
    A::AbstractArray...
) where i
    l = sum([size(A[j], 2) for j in 1:(i - 1)])
    u = l + size(A[i], 2)
    # Using copy materializes the views returned by selectdim
    return copy(u > l + 1 ? selectdim(ȳ, 2, (l+1):u) : selectdim(ȳ, 2, u))
end

@union_intercepts vcat Tuple{Vararg{∇Array}} Tuple{Vararg{AbstractArray}}
function Nabla.∇(
    ::typeof(vcat),
    ::Type{Arg{i}},
    _,
    y,
    ȳ,
    A::AbstractArray...
) where i
    l = sum([size(A[j], 1) for j in 1:(i - 1)])
    u = l + size(A[i], 1)
    return copy(selectdim(ȳ, 1, (l+1):u))
end

@explicit_intercepts fill Tuple{Any, Tuple{Vararg{Integer}}} [true, false]
∇(::typeof(fill), ::Type{Arg{1}}, p, y, ȳ, value, dims...) = sum(ȳ)


text:
struct NewPolicy
    letter::Char
    first_position::Int
    second_position::Int
end

satisfies_policy(pw::AbstractString, p::NewPolicy) = (pw[p.first_position] == p.letter) != (pw[p.second_position] == p.letter)

passwords = split.(readlines("day2/input.txt"))
out = 0
for (range, char_text, password) in passwords
    first_pos, second_pos = parse.(Int, split(range, "-"))
    pw_policy = NewPolicy(char_text[1], first_pos, second_pos)
    global out += satisfies_policy(password, pw_policy)
end

print(out)


text:
import Base.Math: @horner

const FORWARD  =  true
const BACKWARD = false

const sqrtpi = 1.772453850905516027298
const edivsqrt2pi = 1.084437551419227546612

"""
Compute a typed 0.5.
"""
half(x::Number) = oftype(x,0.5)
half(x::Integer) = half(float(x))
half(::Type{T}) where {T<:Number} = convert(T,0.5)
half(::Type{T}) where {T<:Integer} = half(AbstractFloat)

"""
Compute a typed 2.
"""
two(x::Number) = oftype(x,2)
two(::Type{T}) where {T<:Number} = convert(T,2)

"""
The Kronecker ``\\delta`` function:

```math
\\delta_{k,j} = \\left\\{\\begin{array}{ccc} 1 & {\\rm for} & k = j,\\\\ 0 & {\\rm for} & k \\ne j.\\end{array}\\right.
```
"""
δ(k::Integer,j::Integer) = k == j ? 1 : 0

"""
Pochhammer symbol ``(x)_n = \\frac{\\Gamma(x+n)}{\\Gamma(x)}`` for the rising factorial.
"""
function pochhammer(x::Number,n::Integer)
    ret = one(x)
    if n≥0
        for i=0:n-1
            ret *= x+i
        end
    else
        ret /= pochhammer(x+n,-n)
    end
    ret
end

pochhammer(x::Number,n::Number) = isinteger(n) ? pochhammer(x,Int(n)) : ogamma(x)/ogamma(x+n)

function pochhammer(x::Number,n::UnitRange{T}) where T<:Real
    ret = Vector{promote_type(typeof(x),T)}(length(n))
    ret[1] = pochhammer(x,first(n))
    for i=2:length(n)
        ret[i] = (x+n[i]-1)*ret[i-1]
    end
    ret
end

ogamma(x::Number) = (isinteger(x) && x<0) ? zero(float(x)) : inv(gamma(x))

"""
Stirling's asymptotic series for ``\\Gamma(z)``.
"""
stirlingseries(z) = gamma(z)*sqrt((z/π)/2)*exp(z)/z^z

function stirlingseries(z::Float64)
    if z ≥ 3274.12075200175       # N =  4
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273)
    elseif z ≥ 590.1021805526798  # N =  5
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917)
    elseif z ≥ 195.81733962412835 # N =  6
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666)
    elseif z ≥ 91.4692823071966   # N =  7
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5)
    elseif z ≥ 52.70218954633605  # N =  8
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939)
    elseif z ≥ 34.84031591198865  # N =  9
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5)
    elseif z ≥ 25.3173982783047   # N = 10
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873)
    elseif z ≥ 19.685015283078513 # N = 11
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5)
    elseif z ≥ 16.088669099569266 # N = 12
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776)
    elseif z ≥ 13.655055978888104 # N = 13
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583)
    elseif z ≥ 11.93238782087875  # N = 14
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807)
    elseif z ≥ 10.668852439197263 # N = 15
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807,0.0005401647678926045)
    elseif z ≥ 9.715358216638403  # N = 16
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807,0.0005401647678926045,-0.02952788094569912)
    elseif z ≥ 8.979120323411497  # N = 17
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807,0.0005401647678926045,-0.02952788094569912,-0.002481743600264998)
    else
        gamma(z)*sqrt(z/2π)*exp(z)/z^z
    end
end

stirlingremainder(z::Number,N::Int) = (1+zeta(N))*gamma(N)/((2π)^(N+1)*z^N)/stirlingseries(z)

Aratio(n::Int,α::Float64,β::Float64) = exp((n/2+α+1/4)*log1p(-β/(n+α+β+1))+(n/2+β+1/4)*log1p(-α/(n+α+β+1))+(n/2+1/4)*log1p(α/(n+1))+(n/2+1/4)*log1p(β/(n+1)))
Aratio(n::Number,α::Number,β::Number) = (1+(α+1)/n)^(n+α+1/2)*(1+(β+1)/n)^(n+β+1/2)/(1+(α+β+1)/n)^(n+α+β+1/2)/(1+(zero(α)+zero(β))/n)^(n+1/2)

Cratio(n::Int,α::Float64,β::Float64) = exp((n+α+1/2)*log1p((α-β)/(2n+α+β+2))+(n+β+1/2)*log1p((β-α)/(2n+α+β+2))-log1p((α+β+2)/2n)/2)/sqrt(n)
Cratio(n::Number,α::Number,β::Number) = n^(-1/2)*(1+(α+1)/n)^(n+α+1/2)*(1+(β+1)/n)^(n+β+1/2)/(1+(α+β+2)/2n)^(2n+α+β+3/2)

Anαβ(n::Number,α::Number,β::Number) = 2^(α+β+1)/(2n+α+β+1)*exp(lgamma(n+α+1)-lgamma(n+α+β+1)+lgamma(n+β+1)-lgamma(n+1))
function Anαβ(n::Integer,α::Number,β::Number)
    if n==0
        2^(α+β+1)*beta(α+1,β+1)
    else
        val = Anαβ(0,α,β)
        for i=1:n
            val *= (i+α)*(i+β)/(i+α+β)/i*(2i+α+β-1)/(2i+α+β+1)
        end
        val
    end
end

function Anαβ(n::Integer,α::Float64,β::Float64)
    if n+min(α,β,α+β,0) ≥ 7.979120323411497
        2 ^ (α+β+1)/(2n+α+β+1)*stirlingseries(n+α+1)*Aratio(n,α,β)/stirlingseries(n+α+β+1)*stirlingseries(n+β+1)/stirlingseries(n+one(Float64))
    else
        (n+1)*(n+α+β+1)/(n+α+1)/(n+β+1)*Anαβ(n+1,α,β)*((2n+α+β+3)/(2n+α+β+1))
    end
end

"""
The Lambda function ``\\Lambda(z) = \\frac{\\Gamma(z+\\frac{1}{2})}{\\Gamma(z+1)}`` for the ratio of gamma functions.
"""
Λ(z::Number) = exp(lgamma(z+half(z))-lgamma(z+one(z)))
"""
For 64-bit floating-point arithmetic, the Lambda function uses the asymptotic series for ``\\tau`` in Appendix B of

I. Bogaert and B. Michiels and J. Fostier, 𝒪(1) computation of Legendre polynomials and Gauss–Legendre nodes and weights for parallel computing, *SIAM J. Sci. Comput.*, **34**:C83–C101, 2012.
"""
function Λ(x::Float64)
    if x > 9.84475
        xp = x+0.25
        @horner(inv(xp^2),1.0,-1.5625e-02,2.5634765625e-03,-1.2798309326171875e-03,1.343511044979095458984375e-03,-2.432896639220416545867919921875e-03,6.7542375336415716446936130523681640625e-03)/sqrt(xp)
    else
        (x+1.0)*Λ(x+1.0)/(x+0.5)
    end
end

"""
The Lambda function ``\\Lambda(z,λ₁,λ₂) = \\frac{\\Gamma(z+\\lambda_1)}{Γ(z+\\lambda_2)}`` for the ratio of gamma functions.
"""
Λ(z::Number,λ₁::Number,λ₂::Number) = exp(lgamma(z+λ₁)-lgamma(z+λ₂))
function Λ(x::Float64,λ₁::Float64,λ₂::Float64)
    if min(x+λ₁,x+λ₂) ≥ 8.979120323411497
        exp(λ₂-λ₁+(x-.5)*log1p((λ₁-λ₂)/(x+λ₂)))*(x+λ₁)^λ₁/(x+λ₂)^λ₂*stirlingseries(x+λ₁)/stirlingseries(x+λ₂)
    else
        (x+λ₂)/(x+λ₁)*Λ(x+1.,λ₁,λ₂)
    end
end

## TODO: deprecate when Lambert-W is supported in a mainstream repository such as SpecialFunctions.jl
"""
The principal branch of the Lambert-W function, defined by ``x = W_0(x) e^{W_0(x)}``, computed using Halley's method for ``x \\in [-e^{-1},\\infty)``.
"""
function lambertw(x::AbstractFloat)
    if x < -exp(-one(x))
        return throw(DomainError())
    elseif x == -exp(-one(x))
        return -one(x)
    elseif x < 0
        w0 = ℯ*x/(1+inv(inv(sqrt(2*ℯ*x+2))+inv(ℯ-1)-inv(sqrt(2))))
    else
        log1px = log1p(x)
        w0 = log1px*(1-log1p(log1px)/(2+log1px))
    end
    expw0 = exp(w0)
    w1 = w0 - (w0*expw0 - x)/((w0 + 1)*expw0 -
        (w0 + 2) * (w0*expw0 - x)/(2w0 + 2))
    while abs(w1/w0 - 1) > 2eps(typeof(x))
        w0 = w1
        expw0 = exp(w0)
        w1 = w0 - (w0*expw0 - x)/((w0 + 1)*expw0 -
            (w0 + 2) * (w0*expw0 - x)/(2w0 + 2))
    end
    return w1
end
lambertw(x::Real) = lambertw(float(x))

Cnλ(n::Integer,λ::Float64) = 2^λ/sqrtpi*Λ(n+λ)
Cnλ(n::Integer,λ::Number) = 2^λ/sqrt(oftype(λ,π))*Λ(n+λ)
function Cnλ(n::UnitRange{T},λ::Number) where T<:Integer
    ret = Vector{typeof(λ)}(undef, length(n))
    ret[1] = Cnλ(first(n),λ)
    for i=2:length(n)
        ret[i] = (n[i]+λ-half(λ))/(n[i]+λ)*ret[i-1]
    end
    ret
end

function Cnmλ(n::Integer,m::Integer,λ::Number)
    if m == 0
        Cnλ(n,λ)
    else
        (λ+m-1)/2/m*(m-λ)/(n+λ+m)*Cnmλ(n,m-1,λ)
    end
end

function Cnαβ(n::Integer,α::Number,β::Number)
    if n==0
        2^(α+β+1)*beta(α+1,β+1)/π
    else
        val = Cnαβ(0,α,β)
        for i=1:n
            val *= (i+α)*(i+β)/(i+(α+β+1)/2)/(i+(α+β)/2)
        end
        val
    end
end

function Cnαβ(n::Integer,α::Float64,β::Float64)
    if n+min(α,β) ≥ 7.979120323411497
        stirlingseries(n+α+1)/sqrtpi/stirlingseries(2n+α+β+2)*Cratio(n,α,β)*stirlingseries(n+β+1)
    else
        (n+(α+β+3)/2)/(n+β+1)*(n+(α+β+2)/2)/(n+α+1)*Cnαβ(n+1,α,β)
    end
end

function Cnmαβ(n::Integer,m::Integer,α::Number,β::Number)
    if m == 0
        Cnαβ(n,α,β)
    else
        Cnmαβ(n,m-1,α,β)/2(2n+α+β+m+1)
    end
end

function Cnmαβ(n::Integer,m::Integer,α::AbstractArray{T},β::AbstractArray{T}) where T<:Number
    shp = promote_shape(size(α),size(β))
    reshape([ Cnmαβ(n,m,α[i],β[i]) for i in eachindex(α,β) ], shp)
end

function absf(α::Number,β::Number,m::Int,θ::Number)
    ret = zero(θ)
    for l=0:m
        ret += pochhammer(half(α)+α,l)*pochhammer(half(α)-α,l)*pochhammer(half(β)+β,m-l)*pochhammer(half(β)-β,m-l)/factorial(l)/factorial(m-l)/sinpi(θ/2)^(l+α+half(α))/cospi(θ/2)^(m-l+β+half(β))
    end
    ret
end

function absf(α::AbstractArray{T},β::AbstractArray{T},m::Int,θ::Number) where T<:Number
    shp = promote_shape(size(α),size(β))
    reshape([ absf(α[i],β[i],m,θ) for i in eachindex(α,β) ], shp)
end

function absf(α::Number,β::Number,m::Int,θ::AbstractArray{T,1}) where T<:Number
    ret = zero(θ)
    cfs = zeros(T,m+1)
    for l=0:m
        @inbounds cfs[l+1] = pochhammer(half(α)+α,l)*pochhammer(half(α)-α,l)*pochhammer(half(β)+β,m-l)*pochhammer(half(β)-β,m-l)/factorial(l)/factorial(m-l)
    end
    @inbounds for i=1:length(θ),l=0:m
        ret[i] += cfs[l+1]/sinpi(θ[i]/2)^(l+α+half(α))/cospi(θ[i]/2)^(m-l+β+half(β))
    end
    ret
end
absf(α::Number,β::Number,m::Int,θ::AbstractArray{T,2}) where {T<:Number} = [ absf(α,β,m,θ[i,j]) for i=1:size(θ,1), j=1:size(θ,2) ]
absf(α::Number,β::Number,m::Int,θ::AbstractArray{T}) where {T<:Number} = reshape([ absf(α,β,m,θ[i]) for i in eachindex(θ) ], size(θ))

function compute_absf!(ret::Vector{T},cfs::Matrix{T},α::T,β::T,tempcos::Vector{T},tempsin::Vector{T},tempcosβsinα::Vector{T},m::Int) where T<:AbstractFloat
    @inbounds for i=1:length(ret)
        temp = inv(tempcos[i]^m*tempcosβsinα[i])
        ret[i] = cfs[m+1,1]*temp
        for l=1:m
            temp *= tempcos[i]/tempsin[i]
            ret[i] += cfs[m+1,l+1]*temp
        end
    end
    ret
end

function compute_absf!(ret::Vector{T},tempsin::Vector{T},tempsinλ::Vector{T},m::Int) where T<:AbstractFloat
    for i=1:length(ret)
        @inbounds ret[i] = inv(tempsin[i]^m*tempsinλ[i])
    end
    ret
end

function compute_umvm!(um::Vector{T},vm::Vector{T},cfs::Matrix{T},α::T,β::T,tempcos::Vector{T},tempsin::Vector{T},tempcosβsinα::Vector{T},m::Int,θ::Vector{T},ir::UnitRange{Int}) where T<:AbstractFloat
    @inbounds for i in ir
        temp = inv(tempcos[i]^m*tempcosβsinα[i])
        ϑ = (α+half(α))/2-(α+β+m+1)*θ[i]/2
        um[i] = cfs[m+1,1]*cospi(ϑ)*temp
        vm[i] = cfs[m+1,1]*sinpi(ϑ)*temp
        @inbounds for l=1:m
            temp *= tempcos[i]/tempsin[i]
            ϑ = (α+l+half(α))/2-(α+β+m+1)*θ[i]/2
            um[i] += cfs[m+1,l+1]*cospi(ϑ)*temp
            vm[i] += cfs[m+1,l+1]*sinpi(ϑ)*temp
        end
    end
end

function compute_umvm!(um::Vector{T},vm::Vector{T},λ::T,tempsinλm::Vector{T},m::Int,θ::Vector{T},ir::UnitRange{Int}) where T<:AbstractFloat
    @inbounds @simd for i in ir
        temp = inv(tempsinλm[i])
        ϑ = (m+λ)*(half(T)-θ[i])
        um[i] = cospi(ϑ)*temp
        vm[i] = sinpi(ϑ)*temp
    end
end

function findmindices!(Rαβjm::Vector{T},cfs::Matrix{T},α::T,β::T,j::Int,m::Int,tempcos::Vector{T},tempsin::Vector{T},tempcosβsinα::Vector{T}) where T<:AbstractFloat
    compute_absf!(Rαβjm,cfs,α,β,tempcos,tempsin,tempcosβsinα,m)
    rmul!(Rαβjm,Cnmαβ(j,m,α,β))
    rmin,imin = findmin(Rαβjm)
    if rmin < eps(T)
        i₁ = imin-1
        while i₁ ≥ 3
            if Rαβjm[i₁] < eps(T)
                i₁-=1
            else
                break
            end
        end
        i₂ = imin+1
        while i₂ ≤ length(Rαβjm)-2
            if Rαβjm[i₂] < eps(T)
                i₂+=1
            else
                break
            end
        end
        return i₁,i₂
    else
        return 1,0# error("There are no indices such that the interior asymptotic formula is valid.")# but adding an error is type-unstable :(.
    end
end

function findmindices!(Rαβjm::Vector{T},λ::T,j::Int,m::Int,tempsin::Vector{T},tempsinλ::Vector{T}) where T<:AbstractFloat
    compute_absf!(Rαβjm,tempsin,tempsinλ,m)
    rmul!(Rαβjm,Cnmλ(j,m,λ))
    rmin,imin = findmin(Rαβjm)
    if rmin < eps(T)
        i₁ = imin-1
        while i₁ ≥ 3
            if Rαβjm[i₁] < eps(T)
                i₁-=1
            else
                break
            end
        end
        i₂ = imin+1
        while i₂ ≤ length(Rαβjm)-2
            if Rαβjm[i₂] < eps(T)
                i₂+=1
            else
                break
            end
        end
        return i₁,i₂
    else
        return 1,0# error("There are no indices such that the interior asymptotic formula is valid.")# but adding an error is type-unstable :(.
    end
end

# initialization methods

function init_cfs(α::T,β::T,M::Int) where T<:AbstractFloat
    cfs = zeros(T,M+1,M+1)
    @inbounds for m=0:M,l=0:m
        cfs[m+1,l+1] = pochhammer(half(α)+α,l)*pochhammer(half(α)-α,l)*pochhammer(half(β)+β,m-l)*pochhammer(half(β)-β,m-l)/factorial(l)/factorial(m-l)
    end
    cfs
end

function init_c₁c₂!(c₁::Vector,c₂::Vector,a::Vector,b::Vector,j₁::Int,j₂::Int)
    @inbounds for j=1:j₁-1 c₁[j] = 0 end
    @inbounds for j=j₁:j₂ c₁[j] = a[j]*b[j] end
    @inbounds for j=j₂+1:length(c₁) c₁[j] = 0 end
    copyto!(c₂,c₁)
end

function init_c₁c₂!(c₁::Vector,c₂::Vector,u::Vector,v::Vector,c::Vector,i₁::Int,i₂::Int)
    @inbounds for i=1:i₁-1
        c₁[i] = 0
        c₂[i] = 0
    end
    @inbounds for i=i₁:i₂
        c₁[i] = u[i]*c[i]
        c₂[i] = v[i]*c[i]
    end
    @inbounds for i=i₂+1:length(c₁)
        c₁[i] = 0
        c₂[i] = 0
    end
end

"""
Modified Chebyshev moments of the first kind:

```math
    \\int_{-1}^{+1} T_n(x) {\\rm\\,d}x.
```
"""
function chebyshevmoments1(::Type{T}, N::Int) where T
    μ = zeros(T, N)
    for i = 0:2:N-1
        @inbounds μ[i+1] = two(T)/T(1-i^2)
    end
    μ
end

"""
Modified Chebyshev moments of the first kind with respect to the Jacobi weight:

```math
    \\int_{-1}^{+1} T_n(x) (1-x)^\\alpha(1+x)^\\beta{\\rm\\,d}x.
```
"""
function chebyshevjacobimoments1(::Type{T}, N::Int, α, β) where T
    μ = zeros(T, N)
    N > 0 && (μ[1] = 2 .^ (α+β+1)*beta(α+1,β+1))
    if N > 1
        μ[2] = μ[1]*(β-α)/(α+β+2)
        for i=1:N-2
            @inbounds μ[i+2] = (2(β-α)*μ[i+1]-(α+β-i+2)*μ[i])/(α+β+i+2)
        end
    end
    μ
end

"""
Modified Chebyshev moments of the first kind with respect to the logarithmic weight:

```math
    \\int_{-1}^{+1} T_n(x) \\log\\left(\\frac{1-x}{2}\\right){\\rm\\,d}x.
```
"""
function chebyshevlogmoments1(::Type{T}, N::Int) where T
    μ = zeros(T, N)
    N > 0 && (μ[1] = -two(T))
    if N > 1
        μ[2] = -one(T)
        for i=1:N-2
            cst = isodd(i) ? T(4)/T(i^2-4) : T(4)/T(i^2-1)
            @inbounds μ[i+2] = ((i-2)*μ[i]+cst)/(i+2)
        end
    end
    μ
end

"""
Modified Chebyshev moments of the second kind:

```math
    \\int_{-1}^{+1} U_n(x) {\\rm\\,d}x.
```
"""
function chebyshevmoments2(::Type{T}, N::Int) where T
    μ = zeros(T, N)
    for i = 0:2:N-1
        @inbounds μ[i+1] = two(T)/T(i+1)
    end
    μ
end

"""
Modified Chebyshev moments of the second kind with respect to the Jacobi weight:

```math
    \\int_{-1}^{+1} U_n(x) (1-x)^\\alpha(1+x)^\\beta{\\rm\\,d}x.
```
"""
function chebyshevjacobimoments2(::Type{T}, N::Int, α, β) where T
    μ = zeros(T, N)
    N > 0 && (μ[1] = 2 .^ (α+β+1)*beta(α+1,β+1))
    if N > 1
        μ[2] = 2μ[1]*(β-α)/(α+β+2)
        for i=1:N-2
            @inbounds μ[i+2] = (2(β-α)*μ[i+1]-(α+β-i)*μ[i])/(α+β+i+2)
        end
    end
    μ
end

"""
Modified Chebyshev moments of the second kind with respect to the logarithmic weight:

```math
    \\int_{-1}^{+1} U_n(x) \\log\\left(\\frac{1-x}{2}\\right){\\rm\\,d}x.
```
"""
function chebyshevlogmoments2(::Type{T}, N::Int) where T
    μ = chebyshevlogmoments1(T, N)
    if N > 1
        μ[2] *= two(T)
        for i=1:N-2
            @inbounds μ[i+2] = 2μ[i+2] + μ[i]
        end
    end
    μ
end

"""
Compute Jacobi expansion coefficients in ``P_n^{(\\alpha+1,\\beta)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\\alpha,\\beta)}(x)`` in-place.
"""
function incrementα!(c::AbstractVector,α,β)
    αβ,N = α+β,length(c)
    N > 1 && (c[1] -= (β+1)/(αβ+3)*c[2])
    @inbounds for i=2:N-1 c[i] = (αβ+i)/(αβ+2i-1)*c[i] - (β+i)/(αβ+2i+1)*c[i+1] end
    N > 1 && (c[N] *= (αβ+N)/(αβ+2N-1))
    c
end

"""
Compute Jacobi expansion coefficients in ``P_n^{(\\alpha,\\beta+1)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\\alpha,\\beta)}(x)`` in-place.
"""
function incrementβ!(c::AbstractVector,α,β)
    αβ,N = α+β,length(c)
    N > 1 && (c[1] += (α+1)/(αβ+3)*c[2])
    @inbounds for i=2:N-1 c[i] = (αβ+i)/(αβ+2i-1)*c[i] + (α+i)/(αβ+2i+1)*c[i+1] end
    N > 1 && (c[N] *= (αβ+N)/(αβ+2N-1))
    c
end

"""
Compute Jacobi expansion coefficients in ``P_n^{(\\alpha+1,\\alpha+1)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\\alpha,\\alpha)}(x)`` in-place.
"""
function incrementαβ!(c::AbstractVector,α,β)
    @assert α == β
    N = length(c)
    if N == 2
        c[2] *= (2α+2)/(2α+4)
    elseif N > 2
        c[1] -= (α+2)/(4α+10)*c[3]
        @inbounds for i=2:N-2 c[i] = (2α+i)*(2α+i+1)/(2α+2i-1)/(2α+2i)*c[i] - (α+i+1)/(4α+4i+6)*c[i+2] end
        c[N-1] *= (2α+N-1)*(2α+N)/(2α+2N-3)/(2α+2N-2)
        c[N] *= (2α+N)*(2α+N+1)/(2α+2N-1)/(2α+2N)
    end
    c
end

"""
Compute Jacobi expansion coefficients in ``P_n^{(\\alpha-1,\\beta)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\\alpha,\\beta)}(x)`` in-place.
"""
function decrementα!(c::AbstractVector,α,β)
    αβ,N = α+β,length(c)
    N > 1 && (c[N] *= (αβ+2N-2)/(αβ+N-1))
    @inbounds for i=N-1:-1:2 c[i] = (αβ+2i-2)/(αβ+i-1)*c[i] + (αβ+2i-2)/(αβ+2i)*(β+i)/(αβ+i-1)*c[i+1] end
    N > 1 && (c[1] += (β+1)/(αβ+2)*c[2])
    c
end

"""
Compute Jacobi expansion coefficients in ``P_n^{(\\alpha,\\beta-1)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\\alpha,\\beta)}(x)`` in-place.
"""
function decrementβ!(c::AbstractVector,α,β)
    αβ,N = α+β,length(c)
    N > 1 && (c[N] *= (αβ+2N-2)/(αβ+N-1))
    @inbounds for i=N-1:-1:2 c[i] = (αβ+2i-2)/(αβ+i-1)*c[i] - (αβ+2i-2)/(αβ+2i)*(α+i)/(αβ+i-1)*c[i+1] end
    N > 1 && (c[1] -= (α+1)/(αβ+2)*c[2])
    c
end

"""
Compute Jacobi expansion coefficients in ``P_n^{(\\alpha-1,\\alpha-1)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\\alpha,\\alpha)}(x)`` in-place.
"""
function decrementαβ!(c::AbstractVector,α,β)
    @assert α == β
    N = length(c)
    if N == 2
        c[2] *= (2α+2)/(2α)
    elseif N > 2
        c[N] *= (2α+2N-3)*(2α+2N-2)/(2α+N-2)/(2α+N-1)
        c[N-1] *= (2α+2N-5)*(2α+2N-4)/(2α+N-3)/(2α+N-2)
        @inbounds for i=N-2:-1:2 c[i] = (2α+2i-3)*(2α+2i-2)/(2α+i-2)/(2α+i-1)*(c[i] + (α+i)/(4α+4i+2)*c[i+2]) end
        c[1] += (α+1)/(4α+6)*c[3]
    end
    c
end

function modαβ(α)
    if -0.5 < α ≤ 0.5
        a = α
    else
        a = α%1
        a > 0.5 && (a-=1)
        a ≤ -0.5 && (a+=1)
    end
    a
end

function modλ(λ)
    if 0 ≤ λ < 1
        l = λ
    else
        l = λ%1
        l < 0 && (l+=1)
    end
    l
end

function tosquare!(ret::AbstractVector,α,β)
    a,b = modαβ(α),modαβ(β)
    A,B = α-a,β-b
    if α ≤ -0.5 && β ≤ -0.5
        incrementα!(ret,α,β)
        incrementβ!(ret,a,β)
    elseif α ≤ -0.5
        incrementα!(ret,α,β)
        for j=B:-1:1
            decrementβ!(ret,a,j+b)
        end
    elseif β ≤ -0.5
        incrementβ!(ret,α,β)
        for i=A:-1:1
            decrementα!(ret,i+a,b)
        end
    else
        for i=A:-1:1
            decrementα!(ret,i+a,β)
        end
        for j=B:-1:1
            decrementβ!(ret,a,j+b)
        end
    end
    ret
end

function fromsquare!(ret::AbstractVector,α,β)
    a,b = modαβ(α),modαβ(β)
    A,B = α-a,β-b
    if α ≤ -0.5 && β ≤ -0.5
        decrementα!(ret,a,b)
        decrementβ!(ret,α,b)
    elseif α ≤ -0.5
        decrementα!(ret,a,b)
        for j=0:B-1
            incrementβ!(ret,α,j+b)
        end
    elseif β ≤ -0.5
        decrementβ!(ret,a,b)
        for i=0:A-1
            incrementα!(ret,i+a,β)
        end
    else
        for i=0:A-1
            incrementα!(ret,i+a,b)
        end
        for j=0:B-1
            incrementβ!(ret,α,j+b)
        end
    end
    ret
end

function toline!(ret::AbstractVector,α,β)
    @assert α == β
    a,b = modαβ(α),modαβ(β)
    A,B = α-a,β-b
    if α ≤ -0.5
        incrementαβ!(ret,α,β)
    else
        for i=A:-1:1
            decrementαβ!(ret,i+a,i+a)
        end
    end
    ret
end

function fromline!(ret::AbstractVector,α,β)
    @assert α == β
    a,b = modαβ(α),modαβ(β)
    A,B = α-a,β-b
    if α ≤ -0.5
        decrementαβ!(ret,a,b)
    else
        for i=0:A-1
            incrementαβ!(ret,i+a,i+a)
        end
    end
    ret
end


text:
"""
    Send(buf, comm::Comm; dest::Integer, tag::Integer=0)

Perform a blocking send from the buffer `buf` to MPI rank `dest` of communicator
`comm` using the message tag `tag`.

    Send(obj, comm::Comm; dest::Integer, tag::Integer=0)

Complete a blocking send of an `isbits` object `obj` to MPI rank `dest` of
communicator `comm` using with the message tag `tag`.

# External links
$(_doc_external("MPI_Send"))
"""
Send(data, comm::Comm; dest::Integer, tag::Integer=Cint(0)) =
    Send(data, dest, tag, comm)

function Send(buf::Buffer, dest::Integer, tag::Integer, comm::Comm)
    # int MPI_Send(const void* buf, int count, MPI_Datatype datatype, int dest,
    #              int tag, MPI_Comm comm)
    @mpichk ccall((:MPI_Send, libmpi), Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm),
          buf.data, buf.count, buf.datatype, dest, tag, comm)
    return nothing
end
Send(arr::Union{Ref,AbstractArray}, dest::Integer, tag::Integer, comm::Comm) =
    Send(Buffer(arr), dest, tag, comm)
function Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T
    buf = Ref{T}(obj)
    Send(buf, dest, tag, comm)
end

"""
    send(obj, comm::Comm; dest::Integer, tag::Integer=0)

Complete a blocking send using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.
"""
send(obj, comm::Comm; dest::Integer, tag::Integer=0) =
    send(obj, dest, tag, comm)
function send(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Send(buf, dest, tag, comm)
end

"""
    Isend(data, comm::Comm; dest::Integer, tag::Integer=0)

Starts a nonblocking send of `data` to MPI rank `dest` of communicator `comm` using with
the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which [`Buffer_send`](@ref) is defined.

Returns the [`Request`](@ref) object for the nonblocking send.

# External links
$(_doc_external("MPI_Isend"))
"""
Isend(data, comm::Comm; dest::Integer, tag::Integer=0) =
    Isend(data, dest, tag, comm)

function Isend(buf::Buffer, dest::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Isend(const void* buf, int count, MPI_Datatype datatype, int dest,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Isend, libmpi), Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf.data, buf.count, buf.datatype, dest, tag, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return req
end
Isend(data, dest::Integer, tag::Integer, comm::Comm) =
    Isend(Buffer_send(data), dest, tag, comm)

"""
    isend(obj, comm::Comm; dest::Integer, tag::Integer=0)

Starts a nonblocking send of using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.

Returns the commication `Request` for the nonblocking send.
"""
isend(data, comm::Comm; dest::Integer, tag::Integer=0) =
    isend(data, dest, tag, comm)
function isend(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Isend(buf, dest, tag, comm)
end

"""
    data = Recv!(recvbuf, comm::Comm;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)
    data, status = Recv!(recvbuf, comm::Comm, MPI.Status;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)

Completes a blocking receive into the buffer `recvbuf` from MPI rank `source` of communicator
`comm` using with the message tag `tag`.

`recvbuf` can be a [`Buffer`](@ref), or any object for which `Buffer(recvbuf)` is defined.

Optionally returns the [`Status`](@ref) object of the receive.

# See also
- [`Recv`](@ref)
- [`recv`](@ref)

# External links
$(_doc_external("MPI_Recv"))
"""
Recv!(recvbuf, comm::Comm, status=nothing; source=MPI_ANY_SOURCE, tag=MPI_ANY_TAG) =
    Recv!(recvbuf, source, tag, comm, status)

function Recv!(recvbuf::Buffer, source::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status},Nothing})
    # int MPI_Recv(void* buf, int count, MPI_Datatype datatype, int source,
    #              int tag, MPI_Comm comm, MPI_Status *status)
    @mpichk ccall((:MPI_Recv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{Status}),
                  recvbuf.data, recvbuf.count, recvbuf.datatype, source, tag, comm, something(status, MPI_STATUS_IGNORE))
    return recvbuf.data
end
Recv!(recvbuf, source::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status},Nothing}) =
    Recv!(Buffer(recvbuf), source, tag, comm, status)
function Recv!(recvbuf, source::Integer, tag::Integer, comm::Comm, ::Type{Status})
    status = Ref{Status}()
    data = Recv!(recvbuf, source, tag, comm, status)
    return data, status[]
end

"""
    data = Recv(::Type{T}, comm::Comm;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)
    data, status = Recv(::Type{T}, comm::Comm, MPI.Status;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)

Completes a blocking receive of a single `isbits` object of type `T` from MPI rank `source` of
communicator `comm` using with the message tag `tag`.

Returns a tuple of the object of type `T` and optionally the [`Status`](@ref) of the receive.

# See also
- [`Recv!`](@ref)
- [`recv`](@ref)

# External links
$(_doc_external("MPI_Recv"))
"""
Recv(::Type{T}, comm::Comm, status=nothing; source=MPI_ANY_SOURCE, tag=MPI_ANY_TAG) where {T} =
    Recv(T, source, tag, comm, status)

function Recv(::Type{T}, source::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status}, Nothing}) where T
    data = Recv!(Ref{T}(), source, tag, comm, status)
    return data[]
end
function Recv(::Type{T}, source::Integer, tag::Integer, comm::Comm, ::Type{Status}) where T
    status = Ref{Status}()
    val = Recv(T, source, tag, comm, status)
    return val, status[]
end

"""
    obj = recv(comm::Comm;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)
    obj, status = recv(comm::Comm, MPI.Status;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)

Completes a blocking receive of a serialized object from MPI rank `source` of communicator
`comm` using with the message tag `tag`.

Returns the deserialized object and optionally the [`Status`](@ref) of the receive.
"""
recv(comm::Comm, status=nothing; source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG) =
    recv(source, tag, comm, status)
function recv(source::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status}, Nothing})
    stat = Probe(source, tag, comm)
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm, status)
    return MPI.deserialize(buf)
end
function recv(source::Integer, tag::Integer, comm::Comm, ::Type{Status})
    status = Ref{Status}()
    val = recv(source, tag, comm, status)
    return val, status[]
end

"""
    req = Irecv!(recvbuf, comm::Comm;
            source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG)

Starts a nonblocking receive into the buffer `data` from MPI rank `source` of communicator
`comm` using with the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which `Buffer(data)` is defined.

Returns the [`Request`](@ref) object for the nonblocking receive.

# External links
$(_doc_external("MPI_Irecv"))
"""
Irecv!(recvbuf, comm::Comm; source::Integer=MPI_ANY_SOURCE, tag::Integer=MPI_ANY_TAG) =
    Irecv!(recvbuf, source, tag, comm)
function Irecv!(buf::Buffer, source::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Irecv(void* buf, int count, MPI_Datatype datatype, int source,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Irecv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf.data, buf.count, buf.datatype, source, tag, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return req
end
Irecv!(data, source::Integer, tag::Integer, comm::Comm) =
    Irecv!(Buffer(data), source, tag, comm)

function irecv(src::Integer, tag::Integer, comm::Comm)
    (flag, stat) = Iprobe(src, tag, comm)
    if !flag
        return (false, nothing, nothing)
    end
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (true, MPI.deserialize(buf), stat)
end

"""
    data = Sendrecv!(sendbuf, recvbuf, comm;
            dest::Integer, sendtag::Integer=0, source::Integer=MPI_ANY_SOURCE, recvtag::Integer=MPI_ANY_TAG)
    data, status = Sendrecv!(sendbuf, recvbuf, comm, MPI.Status;
            dest::Integer, sendtag::Integer=0, source::Integer=MPI_ANY_SOURCE, recvtag::Integer=MPI_ANY_TAG)

Complete a blocking send-receive operation over the MPI communicator `comm`. Send
`sendcount` elements of type `sendtype` from `sendbuf` to the MPI rank `dest` using message
tag `tag`, and receive `recvcount` elements of type `recvtype` from MPI rank `source` into
the buffer `recvbuf` using message tag `tag`. Return a [`Status`](@ref) object.

If not provided, `sendtype`/`recvtype` and `sendcount`/`recvcount` are derived from the
element type and length of `sendbuf`/`recvbuf`, respectively.

# External links
$(_doc_external("MPI_Sendrecv"))
"""
Sendrecv!(sendbuf, recvbuf, comm::MPI.Comm, status=nothing; dest::Integer, sendtag::Integer=0, source::Integer=MPI_ANY_SOURCE, recvtag::Integer=MPI_ANY_TAG) =
    Sendrecv!(sendbuf, dest, sendtag, recvbuf, source, recvtag, comm, status)

function Sendrecv!(sendbuf::Buffer, dest::Integer, sendtag::Integer,
                   recvbuf::Buffer, source::Integer, recvtag::Integer,
                   comm::Comm, status::Union{Ref{Status}, Nothing})
    # int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, int dest,   int sendtag,
    #                        void *recvbuf, int recvcount, MPI_Datatype recvtype, int source, int recvtag,
    #                    MPI_Comm comm, MPI_Status *status)
    @mpichk ccall((:MPI_Sendrecv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPI_Comm, Ptr{Status}),
                  sendbuf.data, sendbuf.count, sendbuf.datatype, dest, sendtag,
                  recvbuf.data, recvbuf.count, recvbuf.datatype, source, recvtag,
                  comm, something(status, MPI_STATUS_IGNORE))
    return recvbuf.data
end
Sendrecv!(sendbuf, dest::Integer, sendtag::Integer, recvbuf, source::Integer, recvtag::Integer, comm::Comm, status::Union{Ref{Status}, Nothing}) =
    Sendrecv!(Buffer(sendbuf), dest, sendtag, Buffer(recvbuf), source, recvtag, comm, status)

function Sendrecv!(sendbuf, dest::Integer, sendtag::Integer, recvbuf, source::Integer, recvtag::Integer, comm::Comm, ::Type{Status})
    status = Ref{Status}()
    data = Sendrecv!(sendbuf, dest, sendtag, recvbuf, source, recvtag, comm, status)
    return data, status[]
end


text:
# Auxiliary tools to manipulate tuples
"""
    tsetdiff(a::Tuple, b::Tuple)

Assume `b` is completely contained in `a`, and return the tuple with elements of `a` that are
not in `b`.
"""
tsetdiff(a::Tuple, b::Tuple{}) = a
tsetdiff(a::Tuple{Any}, b::Tuple{Any}) = ()
tsetdiff(a::Tuple, b::Tuple{Any}) = a[1] == b[1] ? tail(a) : (a[1], tsetdiff(tail(a), b)...)
tsetdiff(a::Tuple, b::Tuple) = tsetdiff(tsetdiff(a, (b[1],)), tail(b))

@noinline tseterror() = throw(ArgumentError("tuples did not meet requirements"))

"""
    tunique(src::Tuple, dst::Tuple)

Assume that every element appears exactly twice in the combination of `src` and `dst`, and
each elements in `dst` only appears once in `dst`. Ruturn the list which contains every
element only once.
"""
tunique(src::Tuple) = tunique(src, ())
tunique(src::NTuple{N,Any}, dst::NTuple{N,Any}) where {N} = dst
tunique(src::Tuple, dst::Tuple) = src[1] in dst ? tunique((tail(src)..., src[1]), dst) : tunique(tail(src), (dst..., src[1]))

"""
    _findfirst(args...)

Reload the methods `findfirst(args...)` by returning `0` if nothing is found to make it
type stable.
"""
_findfirst(args...) = (i = findfirst(args...); i === nothing ? 0 : i)

"""
    _findnext(args...)

Reload the methods `findnext(args...)` by returning `0` if nothing is found to make it
type stable.
"""
_findnext(args...) = (i = findnext(args...); i === nothing ? 0 : i)

"""
    _findlast(args...)

Reload the methods `findlast(args...)` by returning `0` if nothing is found to make it
type stable.
"""
_findlast(args...) = (i = findlast(args...); i === nothing ? 0 : i)

# Auxiliary method to analyze trace indices
"""
    unique2(itr)

Returns an array containing only those elements that appear exactly once in itr,
and without any elements that appear more than once.
"""
function unique2(itr)
    out = reshape(collect(itr),length(itr))
    i = 1
    while i < length(out)
        inext = _findnext(isequal(out[i]), out, i+1)
        if inext == nothing || inext == 0
            i += 1
            continue
        end
        while !(inext == nothing || inext == 0)
            deleteat!(out,inext)
            inext = _findnext(isequal(out[i]), out, i+1)
        end
        deleteat!(out,i)
    end
    out
end

# Extract index information
"""
    add_indices(IA::NTuple{NA,Any}, IC::NTuple{NC,Any}) where {NA,NC}

Assume `IA` and `IC` has the same length, and `IA` can be obtained from `IC` by permutations.
Return the indices of `IA` where the elements of `IC` appears, i.e., the permutation.
"""
function add_indices(IA::NTuple{NA,Any}, IC::NTuple{NC,Any}) where {NA,NC}
    indCinA = map(l->_findfirst(isequal(l), IA), IC)
    (NA == NC && isperm(indCinA)) || throw(IndexError("invalid index specification: $IA to $IC"))
    return indCinA
end

"""
    trace_indices(IA::NTuple{NA,Any}, IC::NTuple{NC,Any}) where {NA,NC}

`IC` is obtained by tracing out the indices that appear twice in `IA` and then do a
permutation. Return the permutation and two lists which corresponds to the positions of the
first and second traced indices.
"""
function trace_indices(IA::NTuple{NA,Any}, IC::NTuple{NC,Any}) where {NA,NC}
    isodd(length(IA)-length(IC)) && throw(IndexError("invalid trace specification: $IA to $IC"))
    Itrace = tunique(tsetdiff(IA, IC)) # give the indices list that needed to be traced

    cindA1 = map(l->_findfirst(isequal(l), IA), Itrace)
    cindA2 = map(l->_findnext(isequal(l), IA, _findfirst(isequal(l), IA)+1), Itrace)
    indCinA = map(l->_findfirst(isequal(l), IA), IC)

    pA = (indCinA..., cindA1..., cindA2...)
    (isperm(pA) && length(pA) == NA) || throw(IndexError("invalid trace specification: $IA to $IC"))
    return indCinA, cindA1, cindA2
end

"""
    contract_indices(IA::NTuple{NA,Any}, IB::NTuple{NB,Any}, IC::NTuple{NC,Any}) where {NA,NB,NC}

`IC` is obtained by contracting the indices that are shared by `IA` and `IB`, and then do
a permutation. Return `oindA, cindA, oindB, cindB, indCinoAB`, where `oindA`, `oindB` are
positions of indices in `IA` and `IB` that are left open, `cindA`, `cindB` are positions of
indices in `IA` and `IB` that are contracted, `indCinoAB` is the permutation of the indices
`IC` to `(IopenA..., IopenB...)`.
"""
function contract_indices(IA::NTuple{NA,Any}, IB::NTuple{NB,Any}, IC::NTuple{NC,Any}) where {NA,NB,NC}
    # labels
    IAB = (IA..., IB...)
    isodd(length(IAB)-length(IC)) && throw(IndexError("invalid contraction pattern: $IA and $IB to $IC"))
    Icontract = tunique(tsetdiff(IAB, IC))
    IopenA = tsetdiff(IA, Icontract)
    IopenB = tsetdiff(IB, Icontract)

    # to indices
    cindA = map(l->_findfirst(isequal(l), IA), Icontract)
    cindB = map(l->_findfirst(isequal(l), IB), Icontract)
    oindA = map(l->_findfirst(isequal(l), IA), IopenA)
    oindB = map(l->_findfirst(isequal(l), IB), IopenB)
    indCinoAB = map(l->_findfirst(isequal(l), (IopenA..., IopenB...)), IC)

    if !isperm((oindA..., cindA...)) || !isperm((oindB..., cindB...)) || !isperm(indCinoAB)
        throw(IndexError("invalid contraction pattern: $IA and $IB to $IC"))
    end

    return oindA, cindA, oindB, cindB, indCinoAB
end