text:
# Orbit Space Covariance

function get_global_covariance(c1::EPRCoordinate, ot1::Symbol,
                               c2::EPRCoordinate, ot2::Symbol,
                               p::AbstractVector, norms::AbstractVector)

    x = S3(c1.energy, c1.pitch, c1.r)
    y = S3(c2.energy, c2.pitch, c2.r)
    mu = (x .- y)./norms

    GO = Set((:trapped,:potato,:stagnation,:co_passing))
    same_ot = ((ot1 in GO) && (ot2 in GO)) || (ot1 == ot2)

    same_orb = c1 == c2

    np = length(p)
    if np == 2
        k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,mu)/(p[2]^2)) + same_orb*1e-6
    elseif np == 3
        k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,mu)/(p[2]^2)) + same_orb*p[3]^2
    elseif np == 4
        Σ_rho_inv = Diagonal(p[2:4].^(-2))
        k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,Σ_rho_inv*mu)) + same_orb*1e-6
    elseif np == 5
        Σ_rho_inv = Diagonal(p[2:4].^(-2))
        k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,Σ_rho_inv*mu)) + same_orb*p[5]^2
    else
        error("Unsupported number of parameters")
    end

    return k
end

function get_global_covariance(o1::Orbit, o2::Orbit, p::Vector; norms=S3(1.0,1.0,1.0))
    c1 = o1.coordinate
    c2 = o2.coordinate
    ot1 = o1.class
    ot2 = o2.class
    return get_global_covariance(c1,ot1,c2,ot2,p,norms)
end

function get_global_covariance_matrix(orbs, p::Vector; norms=S3(1.0,1.0,1.0))
    n = length(orbs)
    Σ = zeros(n,n)
    for i=1:n
        oi = orbs[i]
        for j = i:n
            oj = orbs[j]
            k = get_global_covariance(oi.coordinate, oi.class,
                                      oj.coordinate, oj.class,
                                      p, norms)
            Σ[i,j] = k
            Σ[j,i] = k
        end
    end
    return Σ
end

function get_global_covariance_matrix(orbs1::Vector, orbs2::Vector, p::Vector; norms=S3(1.0,1.0,1.0))
    n1 = length(orbs1)
    n2 = length(orbs2)
    Σ = zeros(n1,n2)
    for i=1:n1
        oi = orbs1[i]
        for j = 1:n2
            oj = orbs2[j]
            k = get_global_covariance(oi.coordinate, oi.class,
                                      oj.coordinate, oj.class,
                                      p, norms)
            Σ[i,j] = k
        end
    end
    return Σ
end

function epr_cov(energy, pitch, r, orbit_type, p; independent=false, norm = ones(3))
    norb = length(energy)
    np = length(p)

    GO = Set((:trapped,:potato,:stagnation,:co_passing))
    Sigma_F = zeros(norb,norb)
    for i=1:norb
        x = [energy[i], pitch[i], r[i]]
        ox = orbit_type[i]
        for j=i:norb
            y = [energy[j], pitch[j], r[j]]
            oy = orbit_type[j]
            if independent
                same_ot = ox == oy
            else
                same_ot = ((ox in GO) && (oy in GO)) || (ox == oy)
            end
            mu = (x .- y)./norm
            if np == 2
                k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,mu)/(p[2]^2)) + (i==j)*1e-6
            elseif np == 3
                k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,mu)/(p[2]^2)) + (i==j)*p[3]^2
            elseif np == 4
                Σ_rho_inv = Diagonal(p[2:4].^(-2))
                k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,Σ_rho_inv*mu)) + (i==j)*1e-6
            elseif np == 5
                Σ_rho_inv = Diagonal(p[2:4].^(-2))
                k = same_ot*(p[1]^2)*exp(-0.5*dot(mu,Σ_rho_inv*mu)) + (i==j)*p[5]^2
            else
                error("Unsupported number of parameters")
            end

	    Sigma_F[i,j] = k
            Sigma_F[j,i] = k
        end
    end

    return Sigma_F
end

# EPRZ Covariance
struct RepeatedBlockDiagonal{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    data::A
    n::Int
end

Base.size(A::RepeatedBlockDiagonal) = size(A.data).*A.n

function Base.getindex(A::RepeatedBlockDiagonal,I::Int...)
    i, j = I
    ni, nj = size(A.data)
    ii = mod(i-1, ni) + 1
    jj = mod(j-1, nj) + 1

    return div(i-1,ni) == div(j-1,nj) ? A.data[ii,jj] : zero(eltype(A.data))
end

Base.getindex(A::RepeatedBlockDiagonal,I::Int) = A[CartesianIndices(A)[I]]

function RepeatedBlockDiagonal(A::Matrix,n::Int)
    RepeatedBlockDiagonal{eltype(A),2,typeof(A)}(A,n)
end

Base.:(*)(A::RepeatedBlockDiagonal,a::T) where {T<:Number} = RepeatedBlockDiagonal(a*A.data,A.n)
Base.:(*)(a::T, A::RepeatedBlockDiagonal) where {T<:Number} = RepeatedBlockDiagonal(a*A.data,A.n)
Base.:(/)(A::RepeatedBlockDiagonal,a::T) where {T<:Number} = RepeatedBlockDiagonal(A.data/a,A.n)

function Base.:(*)(A::RepeatedBlockDiagonal , a::Union{SparseVector{S,T},SparseMatrixCSC{S,T}}) where {S,T}
    nr,nc = size(A.data)
    colptr = vcat(1:nr:(nr*nc),fill(nr*nc+1,nr*A.n + 1 - nc))
    colptr2 = Array{Int}(undef,length(colptr))
    rowval = repeat(1:nr,nc)
    rowval2 = Array{Int}(undef,length(rowval))
    K = SparseMatrixCSC(size(A)..., colptr, rowval, vec(A.data))

    p = Progress(A.n)
    channel = RemoteChannel(()->Channel{Bool}(A.n), 1)
    B = fetch(@sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end
        @async begin
            B = @distributed (+) for i=1:A.n
                circshift!(colptr2,colptr,(i-1)*nc)
                rowval2 .= rowval .+ (i-1)*nr
                K = SparseMatrixCSC(K.m,K.n,colptr2,rowval2, K.nzval)
                Ka = K*a
                put!(channel,true)
                Ka
            end
            put!(channel, false)
            B
        end
    end)
    return B
end

function Base.:(*)(a::Union{SparseVector{S,T},SparseMatrixCSC{S,T}}, A::RepeatedBlockDiagonal) where {S,T}
    nr,nc = size(A.data)
    colptr = vcat(1:nr:(nr*nc),fill(nr*nc+1,nr*A.n + 1 - nc))
    colptr2 = Array{Int}(undef,length(colptr))
    rowval = repeat(1:nr,nc)
    rowval2 = Array{Int}(undef,length(rowval))
    K = SparseMatrixCSC(size(A)..., colptr, rowval, vec(A.data))

    p = Progress(A.n)
    channel = RemoteChannel(()->Channel{Bool}(A.n), 1)
    B = fetch(@sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end
        @async begin
            B = @distributed (+) for i=1:A.n
                circshift!(colptr2,colptr,(i-1)*nc)
                rowval2 .= rowval .+ (i-1)*nr
                K = SparseMatrixCSC(K.m,K.n,colptr2,rowval2, K.nzval)
                aK = a*K
                put!(channel,true)
                aK
            end
            put!(channel, false)
            B
        end
    end)
    return B
end

function transform_eprz_cov(C,R)
    nr,nc = size(C.data)
    colptr = vcat(1:nr:(nr*nc),fill(nr*nc+1,nr*C.n + 1 - nc))
    colptr2 = Array{Int}(undef,length(colptr))
    rowval = repeat(1:nr,nc)
    rowval2 = Array{Int}(undef,length(rowval))
    K = SparseMatrixCSC(size(C)..., colptr, rowval, vec(C.data))
    S = @distributed (+) for i=1:C.n
        circshift!(colptr2,colptr,(i-1)*nc)
        rowval2 .= rowval .+ (i-1)*nr
        K = SparseMatrixCSC(K.m,K.n,colptr2,rowval2, K.nzval)
        Array((R*(R*K)'))
    end

    return S
end

function transform_eprz_cov2(C,R)
    nr,nc = size(C.data)
    Rnr = size(R)[1]
    Cs = sparse(C.data)

    S = @distributed (+) for i=1:nr:(C.n*nr)
        Ri = R[:,i:(i+nr-1)]
        if nnz(Ri) != 0
            Si = Ri*(Ri*Cs)'
        else
            Si = spzeros(Rnr,Rnr)
        end
        Si
    end

    return Array(S)
end

function transform_eprz_cov3(C,R)
    nr,nc = size(C.data)
    Rnr = size(R)[1]
    Cs = sparse(C.data)

    Ris = SparseMatrixCSC{Float64,Int64}[]
    for i=1:nr:(C.n*nr)
        Ri = R[:,i:(i+nr-1)]
        nnz(Ri) == 0 && continue
        push!(Ris,Ri)
    end

    Si = [zeros(Rnr,Rnr) for i=1:Threads.nthreads()]
    Threads.@threads for Ri in Ris
        id = Threads.threadid()
        Si[id] .= Si[id] .+ Ri*(Ri*Cs)'
    end

    return sum(Si)
end

function LinearAlgebra.inv(A::RepeatedBlockDiagonal)
    RepeatedBlockDiagonal(inv(A.data),A.n)
end

function ep_cov(energy, pitch, sigE, sigp)

    nenergy = length(energy)
    npitch = length(pitch)
    nep = nenergy*npitch

    Σ_ep = zeros(typeof(sigE),nep,nep)
    Σ_p_inv = SMatrix{2,2}(inv(Diagonal([sigE,sigp].^2)))
    subs = CartesianIndices((nenergy,npitch))
    @inbounds for i=1:nep
        ie,ip = Tuple(subs[i])
        x = SVector{2}(energy[ie],pitch[ip])
        for j=i:nep
            je,jp = Tuple(subs[j])
            y = SVector{2}(energy[je],pitch[jp])
            d = x .- y
            l = d'*Σ_p_inv*d
            Σ_ep[i,j] = exp(-0.5*l)
            Σ_ep[j,i] = Σ_ep[i,j]
        end
    end

    return Σ_ep
end

function ep_cov(energy, pitch, p::Vector)
    if length(p) == 2
        return ep_cov(energy, pitch, p[1], p[2])
    elseif length(p) == 3
        return (p[1]^2)*ep_cov(energy, pitch, p[2], p[3])
    else
        error("Unsupported number of parameters")
    end
end

function eprz_cov(energy, pitch, r, z, p::Vector)
    nr = length(r)
    nz = length(z)
    if length(p) == 2
        Σ_ep = ep_cov(energy, pitch, p[1], p[2])
    elseif length(p) == 3
        Σ_ep = (p[1]^2)*ep_cov(energy, pitch, p[2], p[3])
    else
        error("Unsupported number of parameters")
    end
    return RepeatedBlockDiagonal(Σ_ep,nr*nz)
end

function eprz_kernel(x,y,Σ_inv)
    theta_x = acos(clamp(x[2],-1.0,1.0))
    theta_y = acos(clamp(y[2],-1.0,1.0))
    d = S4(x[1] - y[1], # dE
           (theta_x - theta_y)/(2pi), # dtheta/2pi
           x[3] - y[3], # dr
           x[4] - y[4]) # dz
    #d = x .- y
    l = d'*Σ_inv*d
    return exp(-0.5*l)
end

function make_jacobian_spline(J::Vector{T}, t) where T<:Number
    if length(t) == 0
        Js = scale(interpolate(zeros(2),BSpline(Linear())),range(0.0,1.0,length=2))
    else
        Js = scale(interpolate(J, BSpline(Linear())), t)
    end
    return Js
end

make_jacobian_spline(J::T, t) where T<:AbstractInterpolation = J

function get_covariance(o1::OrbitSpline, J1::T, o2::OrbitSpline, J2::T, Σ_inv::S44, atol) where {T<:AbstractInterpolation}
    (length(o1) == 0 || length(o2) == 0) && return 0.0
    K, kerr = hcubature(x -> J1(x[1])*J2(x[2])*eprz_kernel(o1(x[1]), o2(x[2]), Σ_inv),
                        (0.0,0.0), (1.0,1.0), atol=atol)
    return abs(K)
end

function get_covariance(o1::OrbitSpline, J1::T, o2::OrbitSpline, J2::T, sigma::AbstractVector, atol) where {T<:AbstractInterpolation}
    Σ_inv = S44(inv(Diagonal(sigma.^2)))
    return get_covariance(o1, J1, o2, J2, Σ_inv, atol)
end

function get_covariance(M::AxisymmetricEquilibrium, o, sigma;
                        J = Float64[], atol = 1e-3, kwargs...)

    n = length(o)
    t = range(0.0, 1.0, length = n)
    oi = OrbitSpline(o, t)
    if isempty(J)
        J1i = make_jacobian_spline(get_jacobian(M, o; kwargs...),t)
    else
        J1i = make_jacobian_spline(J, t)
    end

    return get_covariance(oi, J1i, oi, J1i, sigma, atol)
end

function get_covariance(M::AxisymmetricEquilibrium, o1, o2, sigma;
                        J1 = Float64[], J2 = Float64[], atol=1e-3, kwargs...)

    n1 = length(o1)
    n2 = length(o2)
    if n1 == 0 || n2 == 0
        return 0.0
    end

    t1 = range(0.0, 1.0, length=n1)
    o1i = OrbitSpline(o1, t1)

    t2 = range(0.0, 1.0, length=n2)
    o2i = OrbitSpline(o2, t2)

    if isempty(J1)
        J1i = make_jacobian_spline(get_jacobian(M, o1; kwargs...), t1)
    else
        J1i = make_jacobian_spline(J1, t1)
    end

    if isempty(J2)
        J2i = make_jacobian_spline(get_jacobian(M, o2; kwargs...), t2)
    else
        J2i = make_jacobian_spline(J2, t2)
    end

    return get_covariance(o1i, J1i, o2i, J2i, sigma, atol)
end

function get_covariance(M::AxisymmetricEquilibrium, c1::T, c2::T, sigma;
                         k1::Float64 = 0.0, k2::Float64 = 0.0,
                         J1 = Float64[], J2 = Float64[],
                        atol = 1e-3, kwargs...) where T <: Union{AbstractOrbitCoordinate,AbstractParticle}
    o1 = get_orbit(M, c1; kwargs...)
    o2 = get_orbit(M, c2; kwargs...)
    return get_covariance(M, o1, o2, sigma; J1 = J1, k1 = k1, J2 = J2, k2 = k2, atol = atol, kwargs...)
end

function get_correlation(M::AxisymmetricEquilibrium, o1, o2, sigma;
                         k1::Float64 = 0.0, k2::Float64 = 0.0,
                         J1 = Float64[], J2 = Float64[],
                         atol = 1e-3, kwargs...)

    n1 = length(o1)
    n2 = length(o2)
    if n1 == 0 || n2 == 0
        return 0.0
    end

    t1 = range(0.0, 1.0, length=n1)
    o1i = OrbitSpline(o1, t1)

    t2 = range(0.0, 1.0, length=n2)
    o2i = OrbitSpline(o2, t2)

    if isempty(J1)
        J1i = make_jacobian_spline(get_jacobian(M, o1; kwargs...), t1)
    else
        J1i = make_jacobian_spline(J1, t1)
    end

    if isempty(J2)
        J2i = make_jacobian_spline(get_jacobian(M, o2; kwargs...), t2)
    else
        J2i = make_jacobian_spline(J2, t2)
    end

    if k1 == 0.0
        k11 = get_covariance(o1i, J1i, o1i, J1i, sigma, atol)
    else
        k11 = k1
    end

    if k2 == 0.0
        k22 = get_covariance(o2i, J2i, o2i, J2i, sigma, atol)
    else
        k22 = k2
    end

    if k11 == 0.0 || k22 == 0.0
        return k11*k22
    end

    k12 = get_covariance(o1i, J1i, o2i, J2i, sigma, atol)

    return k12/sqrt(k11*k22)
end

#1-arg Threaded Covariance
function compute_covariance_matrix(orbs, Js, Σ_inv, atol)
    n = length(orbs)
    N = Threads.nthreads()
    Σ = zeros(n,n)
    Threads.@threads for id=1:N
        @inbounds for i = id:N:n
            oi = orbs[i]
            Ji = Js[i]
            for j=i:n
                oj = orbs[j]
                Jj = Js[j]
                s = get_covariance(oi, Ji, oj, Jj, Σ_inv, atol)
                Σ[i,j] = s
                Σ[j,i] = s
            end
        end
    end

    return Σ
end

#1-arg Threaded Sparse Covariance
function compute_covariance_matrix(orbs, Js, Σ_inv, atol, minval)
    n = length(orbs)
    N = Threads.nthreads()
    Σs = [spzeros(n,n) for i=1:N]
    Threads.@threads for id=1:N
        @inbounds for i=id:N:n
            oi = orbs[i]
            Ji = Js[i]
            Σs[id][i,i] = 1.0
            for j=(i+1):n
                oj = orbs[j]
                Jj = Js[j]
                s = get_covariance(oi, Ji, oj, Jj, Σ_inv, atol)
                if s > minval
                    Σs[id][i,j] = s
                    Σs[id][j,i] = s
                end
            end
        end
    end
    Σ = sum(Σs)

    return Σ
end

#1-arg Distributed Covariance
function compute_covariance_matrix(orbs, Js, Σ_inv, atol, pool::AbstractWorkerPool; batch_size = 0)
    n = length(orbs)
    indices = [(i,j) for i=1:n, j=1:n if i >= j]

    if batch_size == 0
        batch_size = round(Int, length(indices)/(5*nprocs()))
    end

    c = pmap(x -> get_covariance(orbs[x[1]],Js[x[1]],orbs[x[2]],Js[x[2]], Σ_inv, atol),
             pool, indices, on_error=x->0.0, retry_delays=zeros(3), batch_size=batch_size)

    Σ = zeros(n,n)
    for (ii,I) in enumerate(indices)
        i,j = I
        Σ[i,j] = c[ii]
        Σ[j,i] = c[ii]
    end

    return Σ
end

#1-arg Distributed Sparse Covariance
function compute_covariance_matrix(orbs, Js, Σ_inv, atol, minval, pool::AbstractWorkerPool; batch_size=0)
    Σ = compute_covariance_matrix(orbs, Js, Σ_inv, atol, pool, batch_size=batch_size)
    Σ[Σ .< minval] .= zero(eltype(Σ))

    return sparse(Σ)
end

#1-arg Covariance
function get_covariance_matrix(M::AxisymmetricEquilibrium, orbits::Vector, sigma::AbstractVector;
                               Js::Vector{Vector{Float64}} = Vector{Float64}[],
                               sparse::Bool = false, atol::Float64 = 1e-3,
                               minval::Float64 = atol, distributed = false,
                               batch_size = 0, kwargs...)
    n = length(orbits)
    ns = length.(orbits)
    ts = [range(0.0, 1.0, length=nn) for nn in ns]
    orbs = [OrbitSpline(o) for o in orbits]

    if isempty(Js)
        J = Array{Vector{Float64}}(undef, n)
        @inbounds Threads.@threads for i=1:n
            oi = orbits[i]::Orbit
            Ji = get_jacobian(M, oi; kwargs...)
            J[i] = Ji
        end
    else
        J = Js
    end
    Jis = [make_jacobian_spline(jj,tt) for (jj, tt) in zip(J,ts)]

    Σ_inv = S44(inv(Diagonal(sigma.^2)))

    if distributed
        pool = CachingPool(workers())
        if sparse
            Σ = compute_covariance_matrix(orbs, Jis, Σ_inv, atol, minval, pool; batch_size=batch_size)
        else
            Σ = compute_covariance_matrix(orbs, Jis, Σ_inv, atol, pool; batch_size=batch_size)
        end
        clear!(pool)
        return Σ
    else
        if sparse
            return compute_covariance_matrix(orbs, Jis, Σ_inv, atol, minval)
        end
        return compute_covariance_matrix(orbs, Jis, Σ_inv, atol)
    end
end

#1-arg Correlation
function get_correlation_matrix(cov)
    nr,nc = size(cov)
    nr != nc && error("Covariance matrix must be square")
    K = sqrt.(diag(cov))
    return cov./(K*K')
end

function get_correlation_matrix(M::AxisymmetricEquilibrium, orbits::Vector, sigma::AbstractVector;
                                Js::Vector{Vector{Float64}} = Vector{Float64}[],
                                sparse::Bool = false, atol::Float64 = 1e-3,
                                minval::Float64 = atol, distributed = false,
                                batch_size = 0, kwargs...)

    Σ = get_covariance_matrix(M, orbits, sigma; Js=Js, sparse=sparse,
                              atol=atol, minval=minval, distributed=distributed,
                              batch_size=batch_size, kwargs...)

    return get_correlation_matrix(Σ)
end

#2-arg Threaded Covariance
function compute_covariance_matrix(orbs1, Js1, orbs2, Js2, Σ_inv, atol)
    n1 = length(orbs1)
    n2 = length(orbs2)
    Σ = zeros(n1,n2)
    N = Threads.nthreads()
    Threads.@threads for id=1:N
        @inbounds for j=id:N:n2
            oj = orbs2[j]
            Jj = Js2[j]
            for i=1:n1
                oi = orbs1[i]
                Ji = Js1[i]
                Σ[i,j] = get_covariance(oi, Ji, oj, Jj, Σ_inv, atol)
            end
        end
    end
    return Σ
end

#2-arg Threaded Sparse Covariance
function compute_covariance_matrix(orbs1, Js1, orbs2, Js2, Σ_inv, atol, minval)
    n1 = length(orbs1)
    n2 = length(orbs2)
    N = Threads.nthreads()
    Σs = [spzeros(n1,n2) for i=1:N]
    Threads.@threads for id=1:N
        @inbounds for j=id:N:n2
            oj = orbs2[j]
            Jj = Js2[j]
            for i=1:n1
                oi = orbs1[i]
                Ji = Js1[i]
                s = get_covariance(oi, Ji, oj, Jj, Σ_inv, atol)
                if s > minval
                    Σs[id][i,j] = s
                end
            end
        end
    end
    Σ = sum(Σs)
    return Σ
end

#2-arg Distributed Covariance
function compute_covariance_matrix(orbs1, Js1, orbs2, Js2, Σ_inv, atol, pool::AbstractWorkerPool; batch_size=0)
    n1 = length(orbs1)
    n2 = length(orbs2)
    indices = vec([(i,j) for i=1:n1, j=1:n2])

    if batch_size == 0
        batch_size = round(Int, length(indices)/(5*nprocs()))
    end

    Σ = pmap(x -> get_covariance(orbs1[x[1]],Js1[x[1]],orbs2[x[2]],Js2[x[2]], Σ_inv, atol),
             pool, indices, on_error=x->0.0, retry_delays=zeros(3), batch_size=batch_size)

    return reshape(Σ,(n1,n2))
end

#2-arg Distributed Sparse Covariance
function compute_covariance_matrix(orbs1, Js1, orbs2, Js2, Σ_inv, atol, minval, pool::AbstractWorkerPool; batch_size=0)
    Σ = compute_covariance_matrix(orbs1, Js1, orbs2, Js2, Σ_inv, atol, pool; batch_size=batch_size)
    Σ[Σ .< minval] .= zero(eltype(Σ))
    return sparse(Σ)
end

#2-arg Covariance
function get_covariance_matrix(M::AxisymmetricEquilibrium, orbits_1::Vector, orbits_2::Vector, sigma::AbstractVector;
                                Js_1::Vector{Vector{Float64}} = Vector{Float64}[],
                                Js_2::Vector{Vector{Float64}} = Vector{Float64}[],
                                sparse::Bool = false, atol::Float64 = 1e-3,
                                minval::Float64 = atol, distributed=false,
                                batch_size=0, kwargs...)

    Σ_inv = S44(inv(Diagonal(sigma.^2)))

    n1 = length(orbits_1)
    ns1 = length.(orbits_1)
    ts1 = [range(0.0, 1.0, length=nn) for nn in ns1]
    orbs1 = [OrbitSpline(o) for o in orbits_1]

    if isempty(Js_1)
        J1 = Array{Vector{Float64}}(undef, n1)
        @inbounds Threads.@threads for i=1:n1
            oi = orbits_1[i]::Orbit
            Ji = get_jacobian(M, oi; kwargs...)
            J1[i] = Ji
        end
    else
        J1 = Js_1
    end
    Jis1 = [make_jacobian_spline(jj,tt) for (jj, tt) in zip(J1,ts1)]

    n2 = length(orbits_2)
    ns2 = length.(orbits_2)
    ts2 = [range(0.0, 1.0, length=nn) for nn in ns2]
    orbs2 = [OrbitSpline(o) for o in orbits_2]

    if isempty(Js_2)
        J2 = Array{Vector{Float64}}(undef, n2)
        @inbounds Threads.@threads for i=1:n2
            oi = orbits_2[i]::Orbit
            Ji = get_jacobian(M, oi; kwargs...)
            J2[i] = Ji
        end
    else
        J2 = Js_2
    end
    Jis2 = [make_jacobian_spline(jj,tt) for (jj, tt) in zip(J2,ts2)]

    if distributed
        pool = CachingPool(workers())
        if sparse
            Σ = compute_covariance_matrix(orbs1, Jis1, orbs2, Jis2, Σ_inv, atol, minval, pool; batch_size=batch_size)
        else
            Σ = compute_covariance_matrix(orbs1, Jis1, orbs2, Jis2, Σ_inv, atol, pool; batch_size=batch_size)
        end
        clear!(pool)
        return Σ
    else
        if sparse
            return compute_covariance_matrix(orbs1, Jis1, orbs2, Jis2, Σ_inv, atol, minval)
        end
        return compute_covariance_matrix(orbs1, Jis1, orbs2, Jis2, Σ_inv, atol)
    end
end

#2-arg Threaded Correlation
function compute_correlation_matrix(orbs1, Js1, Ks1, orbs2, Js2, Ks2, Σ_inv, atol)
    n1 = length(orbs1)
    n2 = length(orbs2)
    Σ = zeros(n1,n2)
    N = Threads.nthreads()
    Threads.@threads for id=1:N
        @inbounds for j=id:N:n2
            oj = orbs2[j]
            Jj = Js2[j]
            kj = Ks2[j]
            for i=1:n1
                oi = orbs1[i]
                Ji = Js1[i]
                ki = Ks1[i]
                Σ[i,j] = get_covariance(oi, Ji, oj, Jj, Σ_inv, atol)/sqrt(ki*kj)
            end
        end
    end
    return Σ
end

#2-arg Threaded Sparse Correlation
function compute_correlation_matrix(orbs1, Js1, Ks1, orbs2, Js2, Ks2, Σ_inv, atol, minval)
    n1 = length(orbs1)
    n2 = length(orbs2)
    N = Threads.nthreads()
    Σs = [spzeros(n1,n2) for i=1:N]
    Threads.@threads for id=1:N
        @inbounds for j=id:N:n2
            oj = orbs2[j]
            Jj = Js2[j]
            kj = Ks2[j]
            for i=1:n1
                oi = orbs1[i]
                Ji = Js1[i]
                ki = Ks1[i]
                s = get_covariance(oi, Ji, oj, Jj, Σ_inv, atol)/sqrt(ki*kj)
                if s > minval
                    Σs[id][i,j] = s
                end
            end
        end
    end
    Σ = sum(Σs)
    return Σ
end

#2-arg Distributed Correlation
function compute_correlation_matrix(orbs1, Js1, Ks1, orbs2, Js2, Ks2, Σ_inv, atol, pool::AbstractWorkerPool; batch_size=0)
    n1 = length(orbs1)
    n2 = length(orbs2)
    indices = vec([(i,j) for i=1:n1, j=1:n2])

    if batch_size == 0
        batch_size = round(Int, length(indices)/(5*nprocs()))
    end

    Σ = pmap(x -> get_covariance(orbs1[x[1]],Js1[x[1]],orbs2[x[2]],Js2[x[2]], Σ_inv, atol)/sqrt(Ks1[x[1]]*Ks2[x[2]]),
             pool, indices, on_error=x->0.0, retry_delays=zeros(3), batch_size=batch_size)

    return reshape(Σ,(n1,n2))
end

#2-arg Distributed Sparse Correlation
function compute_correlation_matrix(orbs1, Js1, Ks1, orbs2, Js2, Ks2, Σ_inv, atol, minval, pool::AbstractWorkerPool; batch_size=0)
    Σ = compute_correlation_matrix(orbs1, Js1, Ks1, orbs2, Js2, Ks2, Σ_inv, atol, pool; batch_size=batch_size)
    Σ[Σ .< minval] .= zero(eltype(Σ))
    return sparse(Σ)
end

function get_correlation_matrix(M::AxisymmetricEquilibrium, orbits_1::Vector, orbits_2::Vector, sigma::AbstractVector;
                                ks_1::Vector{Float64} = Float64[],
                                Js_1::Vector{Vector{Float64}} = Vector{Float64}[],
                                ks_2::Vector{Float64} = Float64[],
                                Js_2::Vector{Vector{Float64}} = Vector{Float64}[],
                                sparse::Bool = false, atol::Float64 = 1e-3,
                                minval::Float64 = atol, distributed=false,
                                batch_size = 0, kwargs...)

    Σ_inv = S44(inv(Diagonal(sigma.^2)))

    n1 = length(orbits_1)
    ns1 = length.(orbits_1)
    ts1 = [range(0.0, 1.0, length=nn) for nn in ns1]
    orbs1 = [OrbitSpline(o) for o in orbits_1]

    if isempty(Js_1)
        J1 = Array{Vector{Float64}}(undef, n1)
        @inbounds Threads.@threads for i=1:n1
            oi = orbits_1[i]::Orbit
            Ji = get_jacobian(M, oi; kwargs...)
            J1[i] = Ji
        end
    else
        J1 = Js_1
    end
    Jis1 = [make_jacobian_spline(jj,tt) for (jj, tt) in zip(J1,ts1)]

    if isempty(ks_1)
        K1 = zeros(n1)
        @inbounds Threads.@threads for i=1:n1
            oi = orbs1[i]
            Ji = Jis1[i]
            K1[i] = get_covariance(oi, Ji, oi, Ji, Σ_inv, atol)
        end
    else
        K1 = ks_1
    end

    n2 = length(orbits_2)
    ns2 = length.(orbits_2)
    ts2 = [range(0.0, 1.0, length=nn) for nn in ns2]
    orbs2 = [OrbitSpline(o) for o in orbits_2]

    if isempty(Js_2)
        J2 = Array{Vector{Float64}}(undef, n2)
        @inbounds Threads.@threads for i=1:n2
            oi = orbits_2[i]::Orbit
            Ji = get_jacobian(M, oi; kwargs...)
            J2[i] = Ji
        end
    else
        J2 = Js_2
    end
    Jis2 = [make_jacobian_spline(jj,tt) for (jj, tt) in zip(J2,ts2)]

    if isempty(ks_2)
        K2 = zeros(n2)
        @inbounds Threads.@threads for i=1:n2
            oi = orbs2[i]
            Ji = Jis2[i]
            K2[i] = get_covariance(oi, Ji, oi, Ji, Σ_inv, atol)
        end
    else
        K2 = ks_2
    end

    if distributed
        pool = CachingPool(workers())
        if sparse
            Σ = compute_correlation_matrix(orbs1, Jis1, K1, orbs2, Jis2, K2, Σ_inv, atol, minval, pool, batch_size=batch_size)
        else
            Σ = compute_correlation_matrix(orbs1, Jis1, K1, orbs2, Jis2, K2, Σ_inv, atol, pool, batch_size=batch_size)
        end
        clear!(pool)
        return Σ
    else
        if sparse
            return compute_correlation_matrix(orbs1, Jis1, K1, orbs2, Jis2, K2, Σ_inv, atol, minval)
        end
        return compute_correlation_matrix(orbs1, Jis1, K1, orbs2, Jis2, K2, Σ_inv, atol)
    end

end


text:
######### CmdStan program example  ###########

using CmdStan, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) #do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  datatheta = rel_path_cmdstan("..", "examples", "BernoulliInitTheta", "bernoulli.data.R")
  inittheta = rel_path_cmdstan("..", "examples", "BernoulliInitTheta", "bernoulli.init.R")

  global stanmodel
  stanmodel = Stanmodel(name="bernoulli3", model=bernoullimodel,
    output_format=:array, random=CmdStan.Random(seed=-1));

  global rc, sim
  rc, sim, cnames = stan(stanmodel, datatheta, ProjDir,
    init=inittheta, CmdStanDir=CMDSTAN_HOME, summary=false)

  if rc == 0
    println()
    println("Test 0.2 <= mean(theta[1]) <= 0.5)")
    @test 0.2 <= round.(mean(sim[:,8,:]), digits=1) <= 0.5
  end

  #end # cd

ext: jl


text:
InputFolder = './Images/Yeasts/';
OutputFolder = './Results/Images/Yeasts/';
Rescale = 4;

@iA = '*.tif';

@fxg_mRadStdThr [iA] > [M];
params.GRad = 1;
params.RSet = [16 24 32];
params.Fracstd = 0.5;
params.SegThr = 0.6;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
using Plotly;
include("../interps.jl")
f(x) = 2016x^7 + 12x^5 +x^4 + x^3 + x^2 + x + 1

xs = linspace(-1,1, 6)
Newton = interpNewton(xs,map(f,xs))

xss = linspace(-1,1, 10000)
trace1 = scatter(;x=xss, y= map(Newton,xss), name="Newton")
trace2 = scatter(;x=xss, y=  map(f        ,xss), name="f")
plot([trace1, trace2])

#zera czebyszewa
k = 6
xs = []
for i = 0:(k-1)
push!(xs, cos(pi*(2*i+1)/(2*k) ))
end
xs = sort(xs)
@show xs

Newton = interpNewton(xs,map(f,xs))

xss = linspace(-1,1, 10000)
trace1 = scatter(;x=xss, y= map(Newton,xss), name="Newton")
trace2 = scatter(;x=xss, y=  map(f        ,xss), name="f")
plot([trace1, trace2])

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
"""
    calc_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{FreeSlip}, ::ModelSetup)
apply free slip boundary condition for velocity
sets reflective ghost point
"""
@inline function calc_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{FreeSlip},
    ::Union{ModelSetup},
    ::Union{AbstractDiffusion, Nothing},
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρ = state⁻.ρ

    ρu⁻ = state⁻.ρu
    state⁺.ρu = ρu⁻ - 2 * n⁻ ⋅ ρu⁻ .* SVector(n⁻)

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{NoSlip}, ::ModelSetup)
apply no slip boundary condition for velocity
sets reflective ghost point
"""
@inline function calc_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{NoSlip},
    ::Union{ModelSetup},
    ::Union{AbstractDiffusion, Nothing},
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρ = state⁻.ρ
    state⁺.ρu = -state⁻.ρu

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxFirstOrder, ::Insulating, ::HBModel)

apply insulating boundary condition for temperature
sets transmissive ghost point
"""
@inline function calc_boundary_state!(
    ::Union{NumericalFluxFirstOrder},
    ::Union{Insulating, CoupledPrimaryBoundary},
    ::Union{ModelSetup},
    ::Union{AbstractDiffusion, Nothing},
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρθ = state⁻.ρθ

    nothing
end

calc_boundary_state!(
    ::Union{NumericalFluxFirstOrder},
    ::Union{Insulating, CoupledSecondaryBoundary},
    ::Union{SlabLandModelSetup},
    _...,
) = nothing

# """
#     calc_boundary_state!(::NumericalFluxFirstOrder, ::Penetrable{FreeSlip}, ::ModelSetup)
# no mass boundary condition for penetrable
# """
# @inline calc_boundary_state!(
#     ::NumericalFluxFirstOrder,
#     ::Penetrable{FreeSlip},
#     ::ModelSetup,
#     ::ConstantViscosity,
#     _...,
# ) = nothing

# """
#     calc_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{MomentumFlux}, ::HBModel)
# apply kinematic stress boundary condition for velocity
# applies free slip conditions for first-order and gradient fluxes
# """
# @inline function calc_boundary_state!(
#     nf::NumericalFluxFirstOrder,
#     ::Impenetrable{<:MomentumFlux},
#     model::ModelSetup,
#     diff::AbstractDiffusion,
#     args...,
# )
#     return calc_boundary_state!(
#         nf,
#         Impenetrable(FreeSlip()),
#         model,
#         diff,
#         args...,
#     )
# end

# """
#     calc_boundary_state!(::NumericalFluxFirstOrder, ::Penetrable{MomentumFlux}, ::HBModel)
# apply kinematic stress boundary condition for velocity
# applies free slip conditions for first-order and gradient fluxes
# """
# @inline function calc_boundary_state!(
#     nf::NumericalFluxFirstOrder,
#     ::Penetrable{<:MomentumFlux},
#     model::ModelSetup,
#     diff::AbstractDiffusion,
#     args...,
# )
#     return calc_boundary_state!(
#         nf,
#         Penetrable(FreeSlip()),
#         model,
#         diff,
#         args...,
#     )
# end

#     return nothing
# end

# """
#     calc_boundary_state!(::NumericalFluxFirstOrder, ::TemperatureFlux, ::HBModel)

# apply temperature flux boundary condition for velocity
# applies insulating conditions for first-order and gradient fluxes
# """
# @inline function calc_boundary_state!(
#     nf::NumericalFluxFirstOrder,
#     ::TemperatureFlux,
#     model::ModelSetup,
#     args...,
# )
#     return calc_boundary_state!(nf, Insulating(), model, args...)
# end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
function courses_test()
    @testset "Courses" begin
        dCourses = 7.5;
        mCourses = data_to_model_courses(dCourses);
        @test isa(mCourses, Integer)
        dCourses2 = model_to_data_courses(mCourses);
        @test dCourses2 == dCourses;
    end
end

function dollars_test()
	@testset "Dollars" begin
		dDollars = 1000.0;
		mDollars = dollars_data_to_model(dDollars, :perYear);
        dDollars2 = dollars_model_to_data(mDollars, :perYear);
        @test isapprox(dDollars, dDollars2)

		mDollars3 = dollars_data_to_model(dDollars, :perDay);
		@test mDollars3 > 300 * mDollars
		dDollars3 = dollars_model_to_data(mDollars3, :perDay);
		@test dDollars3 ≈ dDollars
	end
end

function time_units_test()
    @testset "Time units" begin
        mtu = [0.7, 0.2];
        validate_mtu(mtu);
        hours = mtu_to_hours_per_week(mtu);
        mtu2 = hours_per_week_to_mtu(hours);
        @test mtu ≈ mtu2

        modelHours = [0.2, 0.7];
        for timeFactor ∈ [:weeksPerYear, :daysPerYear, :hoursPerYear, :hoursPerWeek]
            hd = hours_mtu_to_data(modelHours, timeFactor);
            @test all(hd .> 0.0)
            @test all(hd .< CollegeStratBase.hoursPerYear)
            hm = hours_data_to_mtu(hd, timeFactor);
            @test all(hm .≈ modelHours)
        end

        hd = hours_mtu_to_data(modelHours, :hoursPerWeek);
        hd2 = mtu_to_hours_per_week(modelHours);
        @test all(hd .≈ hd2)
    end
end

@testset "Unit conversion" begin
    courses_test();
    dollars_test();
    time_units_test();
end

# -------------


text:
function flip(direction::Char, boxes::Array{Int})::Array{Int}
    sort!(boxes, rev=direction=='L')
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
@testset "find_traversal & list_traversal" for t in TEST_TREES
    for l in list_traversal(t)
        res = find_traversal(t, walk(t, l))
        @test all(c -> walk(t, c) === walk(t, l), res)
    end
    @test isempty(find_traversal(t, Leaf(-1)))
    @test isempty(find_traversal(t, [t, t]))
    @test isempty(find_traversal(t, BinaryVertex(-1, t, t)))
end


text:
function my_seqsample(n::Int, k::Int)
    k <= n || error("length(x) should not exceed length(a)")

    i = 0
    j = 0
    out = Int[]
    while k > 0
        u = rand()
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        push!(out, i+=1)
        n -= 1
        k -= 1
    end

    return out
end

@time my_seqsample(19,1)

x = collect(1:10)

@time rand_delete!(x, 2)

@time rand_insert!(x, 2, 1:10)

x


text:
module MyPlayground

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) MyPlayground

end # module


text:
# some basic profiling

using Profile
using ProfileView
using NearestNeighborDescent
using MLDatasets

mnist_x = MNIST.convert2features(MNIST.traintensor(Float64))
DescentGraph(rand(5, 100), 10) # precompile
Profile.clear_malloc_data()
DescentGraph(mnist_x, 10)


text: # TODO: write tests for physic


text:
"""
    MPI.Errhandler

An MPI error handler object. Currently only two are supported:

- `ERRORS_ARE_FATAL` (default): program will immediately abort
- `ERRORS_RETURN`: program will throw an `MPIError`.
"""
mutable struct Errhandler
    val::MPI_Errhandler
end
Base.:(==)(a::Errhandler, b::Errhandler) = a.val == b.val
Base.cconvert(::Type{MPI_Errhandler}, errhandler::Errhandler) = errhandler
Base.unsafe_convert(::Type{MPI_Errhandler}, errhandler::Errhandler) = errhandler.val
Base.unsafe_convert(::Type{Ptr{MPI_Errhandler}}, errhandler::Errhandler) = convert(Ptr{MPI_Errhandler}, pointer_from_objref(errhandler))

const ERRHANDLER_NULL  = Errhandler(Consts.MPI_ERRHANDLER_NULL[])
const ERRORS_ARE_FATAL = Errhandler(Consts.MPI_ERRORS_ARE_FATAL[])
const ERRORS_RETURN    = Errhandler(Consts.MPI_ERRORS_RETURN[]   )
add_load_time_hook!(() -> ERRHANDLER_NULL.val  = Consts.MPI_ERRHANDLER_NULL[] )
add_load_time_hook!(() -> ERRORS_ARE_FATAL.val = Consts.MPI_ERRORS_ARE_FATAL[])
add_load_time_hook!(() -> ERRORS_RETURN.val    = Consts.MPI_ERRORS_RETURN[]   )

Errhandler() = Errhandler(ERRHANDLER_NULL.val)

function free(errh::Errhandler)
    if errh != ERRHANDLER_NULL && errh != ERRORS_ARE_FATAL && errh != ERRORS_RETURN && !Finalized()
        # int MPI_Errhandler_free(MPI_Errhandler *errhandler)
        @mpichk ccall((:MPI_Errhandler_free, libmpi), Cint, (Ptr{MPI_Errhandler},), errh)
    end
    return nothing
end

function set_default_error_handler_return()
    set_errorhandler!(COMM_SELF, ERRORS_RETURN)
    set_errorhandler!(COMM_WORLD, ERRORS_RETURN)
end

"""
    MPI.get_errorhandler(comm::MPI.Comm)
    MPI.get_errorhandler(win::MPI.Win)
    MPI.get_errorhandler(file::MPI.File.FileHandle)

Get the current [`Errhandler`](@ref) for the relevant MPI object.

# See also
- [`set_errorhandler!`](@ref)
"""
function get_errorhandler(comm::Comm)
    errh = Errhandler()
    @mpichk ccall((:MPI_Comm_get_errhandler, libmpi), Cint, (MPI_Comm, Ptr{MPI_Errhandler}), comm, errh)
    return errh
end
function get_errorhandler(win::Win)
    errh = Errhandler()
    @mpichk ccall((:MPI_Win_get_errhandler, libmpi), Cint, (MPI_Win, Ptr{MPI_Errhandler}), win, errh)
    return errh
end
function get_errorhandler(file::File.FileHandle)
    errh = Errhandler()
    @mpichk ccall((:MPI_File_get_errhandler, libmpi), Cint, (MPI_File, Ptr{MPI_Errhandler}), file, errh)
    return errh
end

"""
    MPI.set_errorhandler!(comm::MPI.Comm, errh::Errhandler)
    MPI.set_errorhandler!(win::MPI.Win, errh::Errhandler)
    MPI.set_errorhandler!(file::MPI.File.FileHandle, errh::Errhandler)

Set the [`Errhandler`](@ref) for the relevant MPI object.

# See also
- [`get_errorhandler`](@ref)
"""
function set_errorhandler!(comm::Comm, errh::Errhandler)
    @mpichk ccall((:MPI_Comm_set_errhandler, libmpi), Cint, (MPI_Comm, MPI_Errhandler), comm, errh)
    return nothing
end
function set_errorhandler!(win::Win, errh::Errhandler)
    @mpichk ccall((:MPI_Win_set_errhandler, libmpi), Cint, (MPI_Win, MPI_Errhandler), win, errh)
    return nothing
end
function set_errorhandler!(file::File.FileHandle, errh::Errhandler)
    @mpichk ccall((:MPI_File_set_errhandler, libmpi), Cint, (MPI_File, MPI_Errhandler), file, errh)
    return nothing
end


text:
using LinearAlgebra: I, diagm
import KalmanFilters.UnscentedKalmanFilter: unscented_transform

@testset """
Test unscented transform gives original mean and covariance
if σ-points not changed.
""" begin
    for n = 1:5
        x = Float64[1.0 for i = 1:n]
        P = Matrix{Float64}(I, n, n) * 5.0
        σ_parameters = MerweScaled(n=n, α=1e-3, β=2.0, κ=3.0 - n)
        ukf = UKFState(dim_x=n, dim_z=1, σ_parameters=σ_parameters)

        Σ = calculate_σ_points(σ_parameters=σ_parameters, x=x, P=P)
        xn, Pn = unscented_transform(
            ukf=ukf, Σ=Σ, residual_f=ukf.residual_x, mean_f=ukf.mean_x,
        )
        @test all(abs.(x - xn) .< 1e-6)
        @test all(abs.(P - Pn) .< 1e-6)
    end
end

@testset "Test linear" begin
    function fx(
        ;x::AbstractArray{Float64, 1}, δt::Float64,
    )::AbstractArray{Float64, 1}
        F = Float64[[1 δt 0 0];[0 1 0 0];[0 0 1 δt];[0 0 0 1]]
        F * x
    end

    hx(;x::AbstractArray{Float64, 1})::AbstractArray{Float64, 1} = x[[1, 3]]

    δt = 1.0
    σ_parameters = MerweScaled(n=4, α=0.1, β=2.0, κ=-1.0)
    ukf = UKFState(dim_x=4, dim_z=2, σ_parameters=σ_parameters)
    ukf.x .= [0, 0, 0, 0]
    ukf.P .*= 0.01
    Q = nothing
    R = diagm([0.01 ^ 2, 0.01 ^ 2])

    n = 10
    target_state = Float64[n, 1., n, 1.]
    measurements = vcat([Float64[i i] for i = 1:n]...)
    for i = 1:n
        predict!(ukf=ukf, δt=δt, Q=Q, fx=fx)
        update!(ukf=ukf, z=measurements[i, :], R=R, hx=hx)
    end
    @test all(abs.(ukf.x - target_state) .< 1e-2)
end

max_stars_repo_name: pxl-th/KalmanFilters

max_issues_repo_name: pxl-th/KalmanFilters

max_forks_repo_name: pxl-th/KalmanFilters


text:
@inline function Double64(str::S) where {S<:AbstractString}
  oldprec = precision(BigFloat)
  setprecision(BigFloat, 254)
  result = Double64(BigFloat(str))
  setprecision(BigFloat, oldprec)
  return result
end

macro df64_str(val::AbstractString)
  :(Double64($val))
end

@inline function Double32(str::S) where {S<:AbstractString}
  oldprec = precision(BigFloat)
  setprecision(BigFloat, 126)
  result = Double32(BigFloat(str))
  setprecision(BigFloat, oldprec)
  return result
end

macro df32_str(val::AbstractString)
  :(Double32($val))
end

@inline function Double16(str::S) where {S<:AbstractString}
  oldprec = precision(BigFloat)
  setprecision(BigFloat, 62)
  result = Double16(BigFloat(str))
  setprecision(BigFloat, oldprec)
  return result
end

macro df16_str(val::AbstractString)
  :(Double16($val))
end

function tryparse(::Type{DoubleFloat{Float64}}, str::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(BigFloat, str, base=base)
    return ifelse(isnothing(x), nothing, DoubleFloat{Float64}(x))
end

function tryparse(::Type{DoubleFloat{Float32}}, str::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(BigFloat, str, base=base)
    return ifelse(isnothing(x), nothing, DoubleFloat{Float32}(x))
end

function tryparse(::Type{DoubleFloat{Float16}}, str::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(BigFloat, str, base=base)
    return ifelse(isnothing(x), nothing, DoubleFloat{Float16}(x))
end

function tryparse(::Type{Complex{DoubleFloat{Float64}}}, str::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(Complex{BigFloat}, str)
    return ifelse(isnothing(x), nothing, Complex{DoubleFloat{Float64}}(x))
end

function tryparse(::Type{Complex{DoubleFloat{Float32}}}, str::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(Complex{BigFloat}, str)
    return ifelse(isnothing(x), nothing, Complex{DoubleFloat{Float32}}(x))
end

function tryparse(::Type{Complex{DoubleFloat{Float16}}}, str::S; base::Int=10) where {S<:AbstractString}
    x = tryparse(Complex{BigFloat}, str)
    return ifelse(isnothing(x), nothing, Complex{DoubleFloat{Float16}}(x))
end

@inline function Complex{Double64}(str::S) where {S<:AbstractString}
  oldprec = precision(BigFloat)
  setprecision(BigFloat, 508)
  cplxbf = tryparse(Complex{BigFloat}, str)
  result = Complex{Double64}(Double64(real(cplxbf)), Double64(imag(cplxbf)))
  setprecision(BigFloat, oldprec)
  return result
end

@inline function Complex{Double32}(str::S) where {S<:AbstractString}
  oldprec = precision(BigFloat)
  setprecision(BigFloat, 254)
  cplxbf = tryparse(Complex{BigFloat}, str)
  result = Complex{Double32}(Double32(real(cplxbf)), Double32(imag(cplxbf)))
  setprecision(BigFloat, oldprec)
  return result
end

@inline function Complex{Double16}(str::S) where {S<:AbstractString}
  oldprec = precision(BigFloat)
  setprecision(BigFloat, 120)
  cplxbf = tryparse(Complex{BigFloat}, str)
  result = Complex{Double16}(Double16(real(cplxbf)), Double16(imag(cplxbf)))
  setprecision(BigFloat, oldprec)
  return result
end


text:
# This file was generated, do not modify it. # hide
Xs = source()
ys = source()

atom = @load DecisionTreeRegressor
atom.n_subfeatures = 4 # to ensure diversity among trained atomic models

machines = (machine(atom, Xs, ys) for i in 1:100)


text:
function getattribute(tiffdata::EzXML.Node, ::Type{Int}, attr::String)
    try
        parse(Int, tiffdata[attr])+1
    catch e
        (!isa(e, KeyError)) && rethrow(e)
        -1
    end
end

function getattribute(tiffdata::EzXML.Node, ::Type{Float64}, attr::String)
    try
        parse(Float64, tiffdata[attr])
    catch e
        (!isa(e, KeyError)) && rethrow(e)
        NaN
    end
end

"""
    ifdindex!(ifds, image, dims, imageidx)

Update the master ifd index list, `ifds`, using the TiffData's in `image`.
`dims` is a NamedTuple of the size of each dimension of `image` in the order
specified by the `DimensionOrder` parameter. `imageidx` is positive integer
corresponding to the index of the current image in the OME-TIFF file.

When we read the `TiffFile` we'll know what indices in the 6D matrix each IFD
belongs to.
"""
function ifdindex!(ifd_index::Array{Union{NTuple{4, Int}, Nothing}},
                   ifd_files::Array{Union{Tuple{String, String}, Nothing}},
                   obs_filepaths::Set{String},
                   image::EzXML.Node,
                   dims::NamedTuple,
                   filepath::String,
                   imageidx::Int)

    tiffdatas = findall(".//ns:TiffData", image, ["ns"=>namespace(image)])

    ifd = 1
    # this is an offset value since multiple ifds can share the same index if
    # they are split across files, IFD1 (File1), IFD1 (File2), etc
    file_ifd_offset = 1
    for tiffdata in tiffdatas
        try # if this tiffdata specifies the corresponding IFD
            ifd = parse(Int, tiffdata["IFD"]) + 1
        catch
            ifd = 1
        end

        uuid_node = findfirst("./ns:UUID", tiffdata, ["ns"=>namespace(tiffdata)])
        if uuid_node != nothing
            uuid = nodecontent(uuid_node)
            filepath = joinpath(dirname(filepath), uuid_node["FileName"])
            # if this file isn't one we've observed before, increment the offset
            if !in(filepath, obs_filepaths)
                ifd = file_ifd_offset
                file_ifd_offset += 1
                push!(obs_filepaths, filepath)
            end
            ifd_files[ifd] = (uuid, filepath)
        end

        # get Z, C, T indices (in order specified by `dims`)
        indices = Tuple(getattribute(tiffdata, Int, "First$x") for x in keys(dims)[3:5])
        # how many ifds does this tiffdata correspond to
        p = getattribute(tiffdata, Int, "PlaneCount") - 1

        # if none of the Z, C, T indices are specified then we'll assume the
        # indices starting with the inner dimension, etc
        if all(indices .< 0)
            # index in the master ifd list
            idx = ifd
            # reverse iterate since we cycle the inner dimension the most
            for k=1:dims[5], j=1:dims[4], i=1:dims[3]
                ifd_index[idx] = (i, j, k, imageidx)
                ifd_files[idx] = nothing
                # if this tiffdata applies to multiple ifds then check that we
                # don't exceed the specified number of ifds
                (p > 1 && idx >= p+ifd-1) && break
                idx += 1
            end
        # if any of the indices are specified in the tiffdata then use these
        else
            indices = (indices..., imageidx) # add the position index
            # all the indices that are not specified, we assume the first index
            ifd_index[ifd] = Tuple(pos > 0 ? pos : 1 for pos in indices)
        end
    end
end

"""
    get_unitful_axis(image, dimsize, stepsize, units)

Attempts to return a unitful axis with a length of `dimsize`. `stepsize` and
`units` should be the XML tags in `image`.
"""
function get_unitful_axis(image::EzXML.Node, dimsize::Int, stepsize::String, units::String)
    try
        # OME-XML stores the step size
        increment = parse(Float64, image[stepsize])

        # This is an ugly hack to convert the unit string into Unitful.Unit till
        # https://github.com/PainterQubits/Unitful.jl/issues/214 gets fixed
        unittype = @eval @u_str $(image[units])

        # Create a unitful range
        return 0*unittype:increment*unittype:increment*(dimsize-1)*unittype
    catch
        # alternatively, just return the bare unitless range
        return 1:dimsize
    end
end

const axis_name_mapping = (X = :x, Y = :y, Z=:z, T=:time, C=:channel, P=:position)
"""
    build_axes(omexml::EzXML.Node)

Returns an array of ints with dimension sizes and an array of `AxisArrays.Axis`
objects both in XYZCT order given the Pixels node of the OME-XML document
"""
function build_axes(image::EzXML.Node)
    order = "YX"*join(replace(split(image["DimensionOrder"], ""), "X"=>"", "Y"=>""))
    order = Tuple(Symbol(dim) for dim in order)
    dims = NamedTuple{order, NTuple{5, Int}}(Tuple(parse(Int, image["Size$(x)"]) for x in order))
    dims = merge(dims, [:P=>1])

    # extract channel names
    channel_names = nodecontent.(findall("ns:Channel/@Name", image, ["ns"=>namespace(image)]))
    if isempty(channel_names)
        channel_names = ["C$x" for x in 1:dims[:C]]
    end

    vals = (
        X=get_unitful_axis(image, dims[:X], "PhysicalSizeX", "PhysicalSizeXUnit"),
        Y=get_unitful_axis(image, dims[:Y], "PhysicalSizeY", "PhysicalSizeYUnit"),
        Z=get_unitful_axis(image, dims[:Z], "PhysicalSizeZ", "PhysicalSizeZUnit"),
        C=to_symbol.(channel_names),
        T=get_unitful_axis(image, dims[:T], "TimeIncrement", "TimeIncrementUnit"),
        P=1
    )

    axes = [Axis{axis_name_mapping[key]}(vals[key]) for key in keys(dims)]

    dims, axes
end

"""
    get_elapsed_times(containers, master_dims, masteraxis; default_unit)

Extracts the actual acquisition times from the OME-XML data. Takes `containers`, a vector of the XML nodes
corresponding to the root of each image.
"""
function get_elapsed_times(containers::Vector{EzXML.Node}, master_dims::NamedTuple, masteraxis::Vector{Axis}; default_unit=Unitful.s)

    # get the used dims that aren't x or y (since each xy plane has one elapsed time)
    nonxydims = NamedTuple{Tuple(k for (k, v) in pairs(master_dims) if k != :X && k != :Y && v > 1)}(master_dims)
    elapsed_times = fill(NaN*default_unit, values(nonxydims))

    # OMETIFF stores the Z, T, C info for each plane in TheZ, TheT, and TheC attributes
    attrnames = ["The$k" for k in keys(nonxydims) if k != :P]
    didx = fill(1, length(nonxydims))

    # we have to track position independently since it isn't stored in the plane tags
    # get which index corresponds to position (usually it's last)
    posdim = findfirst(x->x==:P, keys(nonxydims))
    unitstr = string(default_unit)

    for (pos, container) in enumerate(containers)
        planes = findall("ns:Plane[@DeltaT]", container, ["ns"=>namespace(container)])

        for plane in planes
            for d in 1:length(didx)
                # if this is the position dimension, get it from the current Pixel object since it's
                # not stored in the plane itself
                if d == posdim
                    didx[d] = pos
                else
                    didx[d] = getattribute(plane, Int, attrnames[d])
                end
            end
            try
                unitstr = plane["DeltaTUnit"]
            catch
            end
            unittype = @eval @u_str $unitstr
            elapsed_times[didx...] = getattribute(plane, Float64, "DeltaT")*unittype
        end
    end
    used_axes = (masteraxis[i] for i in findall(x->x in keys(nonxydims), keys(master_dims)))
    AxisArray(elapsed_times, used_axes...)
end


text:
push!(LOAD_PATH, "./models")
push!(LOAD_PATH, ".")
module Vessel_network
# Module that implements methods to facilitate the generation and
# manipulation of blood vessel structures.

using Interpolations
using Roots
using QuadGK

export Tree, addchild, Absolute_Vessel, displacement, evaluation, Uniform_sampler

type TreeNode
    parent::Int
    children::Vector{Int}
    curveData::Array{Any,1}
end
type Tree
    nodes::Vector{TreeNode}
end
Tree(data) = Tree([TreeNode(0, Vector{Int}(), data)])
function addchild(tree::Tree, id::Int, curveData::Array{Any,1})
    1 <= id <= length(tree.nodes) || throw(BoundsError(tree, id))
    push!(tree.nodes, TreeNode(id, Vector{}(), curveData))
    child = length(tree.nodes)
    push!(tree.nodes[id].children, child)
    child
end
children(tree, id) = tree.nodes[id].children
parent(tree,id) = tree.nodes[id].parent
data(tree,id) = tree.nodes[id].curveData

function angle_direction(direction)
    return pi/2 - atan(direction[2],direction[1])
end

function segment_generator(Vessel, id)
    parameters = data(Vessel,id)
    x0 = parameters[1]
    tf = parameters[2]
    angle = parameters[3]
    speed = parameters[4]
    curve_type = parameters[5]
    params = parameters[6]
    if curve_type == "straight"
        f  = t -> [0 ; t]
        df = t -> [0 ; 1]
    elseif curve_type == "ellipse"
        if length(params)!= 3
            error("Incorrect number of input parameters")
        end
        side = params[1]
        x_semiax = params[2]
        y_semiax = params[3]
        g(t) = [x_semiax*(1-cos(t))*side; y_semiax*sin(t)]
        dg(t) = [x_semiax*side*sin(t); y_semiax*cos(t) ]
        # now we need to obtain the natural parametrization
            norm_derivative(t) = sqrt(dg(t)[1]^2 + dg(t)[2]^2)
            arclength(t) = quadgk(norm_derivative, 0, t)
            # invert the arc length
            max_sample = 1
            while arclength(max_sample)[1]-tf < 0
                max_sample = max_sample*2
            end
            t_sample = [t for t = linspace(0,max_sample,1000)]
            arclength_sample = [arclength(t)[1] for t in t_sample]
            inv_interp = LinearInterpolation(arclength_sample, t_sample)
        f = t -> g(inv_interp(t))
        df = t-> dg(inv_interp(t))/sqrt( dg(inv_interp(t))[1]^2 + dg(inv_interp(t))[2]^2 )
    else
        error("Not a valid curve type, find the correct names")
    end
    # Rotation matrix
    rotation_mat = [cos(angle) -sin(angle); sin(angle) cos(angle)]
    # Translation, speed scalation, rotation
    segment(t) = rotation_mat*f(t*speed) + x0
    dsegment(t) = rotation_mat*df(t*speed)
    return [segment , dsegment]
end

function Absolute_Vessel(Vessel,id, x0, angle)
    datta = data(Vessel,id)
    datta[2] = datta[2]-angle
    prepend!(datta,[x0])
    Vessel.nodes[id].curveData = datta
    segt, dsegt = segment_generator(Vessel, id)
    new_angle = angle_direction(dsegt(datta[2]))
    new_position = segt(datta[2])
    for j in children(Vessel,id)
        Absolute_Vessel(Vessel,j, new_position, new_angle)
    end
end

# Function to move along the vessels, from a starting position
# (id, t) and a displacement dt
function displacement(Vessel, id, t, dt)
    # This method is not robust to big displacements !
    tf = data(Vessel,id)[2]
    if 0 <= t+dt && t+dt <= tf
        return (id, t+dt)
    elseif t+dt > tf
        next_vessel = children(Vessel,id)
        if isempty(next_vessel)
            return NaN
        else
        # Choose some children at random, assuming at most 2
            new_id = next_vessel[rand(1:2)]
            return (new_id, t+dt-tf)
        end
    else
        # backward case, t+dt < 0
        new_id = parent(Vessel,id)
        if new_id == 0
            return NaN
        else
            new_tf = data(Vessel,new_id)[2]
            return (new_id, new_tf+t+dt)
        end
    end
end

#function to convert the position (id,t) to a point in R2.
function evaluation(Vessel, id, t)
    segt, _ = segment_generator(Vessel, id)
    return segt(t)
end

# Uniform samples on the network
function Uniform_sampler(Vessel)
    times = []
    for id =1:length(Vessel.nodes)
        append!(times, data(Vessel,id)[2])
    end
    tot_time = sum(times)
    selector = rand()*tot_time
    for i=1:length(times)
        if selector < times[i]
            return (i, selector)
        else
            selector = selector - times[i]
        end
    end
end
end


text:
using Test
using carnivalparty

@testset "My test" begin
  @test 1 + 1 ≈ 2
end

max_stars_repo_name: bischtob/carnivalparty

max_issues_repo_name: bischtob/carnivalparty

max_forks_repo_name: bischtob/carnivalparty


text:
@testset "Rational elliptic curves" begin

  curves_to_test_tor_short = [
  ([1, 2], 4),
  ([2, 3], 2),
  ([0, 1], 6),
  ([1, 0], 2),
  ([1, 1], 1)
  ]

  curves_to_test_tor_struc =
  [
  ([0, -1, 1, -7820, -263580], [1]),
  ([1, 0, 1, -2731, -55146], [2]),
  ([0, 1, 1, -9, -15], [3]),
  ([1, 1, 1, -80, 242], [4]),
  ([0, -1, 1, -10, -20], [5]),
  ([1, 0, 1, -36, -70], [6]),
  ([1, -1, 1, -3, 3], [7]),
  ([1, 1, 1, 35, -28], [8]),
  ([1, -1, 1, -14, 29], [9]),
  ([1, 0, 0, -45, 81], [10]),
  ([1, -1, 1, -122, 1721], [12]),
  ([1, 1, 1, -135, -660], [2, 2]),
  ([1, 1, 1, -10, -10], [2, 4]),
  ([1, 0, 1, -19, 26], [2, 6]),
  ([1, 0, 0, -1070, 7812], [2, 8]),
  ]

  @testset "Point order computation" begin
    E = EllipticCurve([0, -1, 1, -7820, -263580])
    @test 1 == @inferred order(infinity(E))

    E = EllipticCurve([1, 0, 1, -2731, -55146])
    @test 2 == @inferred order(E([fmpq(-121, 4), fmpq(117, 8)]))

    E = EllipticCurve([0, 1, 1, -9, -15])
    @test 3 == @inferred order(E([5, 9]))

    E = EllipticCurve([1, 1, 1, -80, 242])
    @test 4 == @inferred order(E([5, -2]))

    E = EllipticCurve([0, -1, 1, -10, -20])
    @test 5 == @inferred order(E([5, 5]))

    E = EllipticCurve([1, 0, 1, -36, -70])
    @test 6 == @inferred order(E([-4, 5]))

    E = EllipticCurve([1, -1, 1, -3, 3])
    @test 7 == @inferred order(E([1, 0]))

    E = EllipticCurve([1, 1, 1, 35, -28])
    @test 8 == @inferred order(E([2, 6]))

    E = EllipticCurve([1, -1, 1, -14, 29])
    @test 9 == @inferred order(E([-3, 7]))

    E = EllipticCurve([1, 0, 0, -45, 81])
    @test 10 == @inferred order(E([0, 9]))

    E = EllipticCurve([1, -1, 1, -122, 1721])
    @test 12 == @inferred order(E([-9, 49]))
  end

  @testset "Torsion test" begin
    E = EllipticCurve([1, -1, 1, -19353, 958713])
    @test @inferred istorsion_point(E([103,172]))
    @test @inferred !istorsion_point(E([-121,1292]))
  end

  @testset "Torsion points (Lutz-Nagell)" begin
    for c in curves_to_test_tor_short
      E = EllipticCurve(c[1])
      T = @inferred torsion_points_lutz_nagell(E)
      @test c[2] == length(T)
    end
  end

  @testset "Torsion points (division polynomials" begin
    for c in curves_to_test_tor_short
      E = EllipticCurve(c[1])
      T = @inferred torsion_points_division_poly(E)
      @test c[2] == length(T)
    end

    for c in curves_to_test_tor_struc
      E = EllipticCurve(c[1])
      T = @inferred torsion_points_division_poly(E)
      @test prod(c[2]) == length(T)
    end
  end

  @testset "Torsion points" begin
    for c in curves_to_test_tor_short
      E = EllipticCurve(c[1])
      T = @inferred torsion_points(E)
      @test c[2] == length(T)
      # test the caching
      T = @inferred torsion_points(E)
      @test c[2] == length(T)
    end

    for c in curves_to_test_tor_struc
      E = EllipticCurve(c[1])
      T = @inferred torsion_points(E)
      @test prod(c[2]) == length(T)
      # test the caching
      T = @inferred torsion_points(E)
      @test prod(c[2]) == length(T)
    end
  end

  @testset "Torsion point structure" begin
    for c in curves_to_test_tor_struc
      E = EllipticCurve(c[1])
      T = @inferred torsion_structure(E)
      @test T[1] == c[2]
      @test order(T[2][1]) == T[1][1]
      if length(T[1]) == 2
        @test order(T[2][2]) == T[1][2]
      end
      # test the caching
      T = @inferred torsion_structure(E)
      @test T[1] == c[2]
      @test order(T[2][1]) == T[1][1]
      if length(T[1]) == 2
        @test order(T[2][2]) == T[1][2]
      end
    end
  end

  @testset "Minimal model (Laska-Kraus-Connell" begin
    E = EllipticCurve([1,2,3,4,5])
    EE = @inferred laska_kraus_connell(E)
    @test EE.coeff == [1, -1, 0, 4, 3]

    E = EllipticCurve([625, -15625, 19531250, -2929687500, -34332275390625])
    EE = @inferred laska_kraus_connell(E)
    @test EE.coeff == [1, -1, 0, 4, 3]
  end

  @testset "Tates algorithm" begin
    E = EllipticCurve([625, -15625, 19531250, -2929687500, -34332275390625])
    EE = @inferred tates_algorithm_global(E)
    @test EE.coeff == [1, -1, 0, 4, 3]

    E = EllipticCurve([1,2,3,4,5])
    EE = @inferred tates_algorithm_global(E)
    @test EE.coeff == [1, -1, 0, 4, 3]

    #  25350.a1
    E = EllipticCurve([1, 1, 0, 40050, 7557750])
    Ep, K, f, c = tates_algorithm_local(E, 2)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "I1"
    @test f == 1
    @test c == 1

    Ep, K, f, c = tates_algorithm_local(E, 3)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "I2"
    @test f == 1
    @test c == 2

    Ep, K, f, c = tates_algorithm_local(E, 5)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "III*"
    @test f == 2
    @test c == 2

    Ep, K, f, c = tates_algorithm_local(E, 13)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "IV*"
    @test f == 2
    @test c == 1

    # 150.a1
    E = EllipticCurve([1, 1, 0, -20700, 1134000])
    Ep, K, f, c = tates_algorithm_local(E, 2)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "I5"
    @test f == 1
    @test c == 1

    Ep, K, f, c = tates_algorithm_local(E, 3)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "I10"
    @test f == 1
    @test c == 2

    Ep, K, f, c = tates_algorithm_local(E, 5)
    @test tidy_model(Ep).coeff == E.coeff
    @test K == "III*"
    @test f == 2
    @test c == 2
  end
end


text:
using ClimateMachine.Mesh.Filters
using KernelAbstractions
using KernelAbstractions.Extras: @unroll

import ClimateMachine.Mesh.Filters: apply_async!
import ClimateMachine.Mesh.Filters: AbstractFilterTarget
import ClimateMachine.Mesh.Filters:
    number_state_filtered, vars_state_filtered, compute_filter_argument!, compute_filter_result!

function modified_filter_matrix(r, Nc, σ)
    N = length(r) - 1
    T = eltype(r)

    @assert N >= 0
    @assert 0 <= Nc

    a, b = GaussQuadrature.legendre_coefs(T, N)
    V = (N == 0 ? ones(T, 1, 1) : GaussQuadrature.orthonormal_poly(r, a, b))

    Σ = ones(T, N + 1)
    if Nc ≤ N
        Σ[(Nc:N) .+ 1] .= σ.(((Nc:N) .- Nc) ./ (N - Nc))
    end

    V * Diagonal(Σ) / V
end

abstract type AbstractMassPreservingSpectralFilter <: AbstractSpectralFilter end
"""
    CutoffFilter(grid, Nc=polynomialorders(grid))
Returns the spectral filter that zeros out polynomial modes greater than or
equal to `Nc`.
"""
struct MassPreservingCutoffFilter{FM} <: AbstractMassPreservingSpectralFilter
    "filter matrices in all directions (tuple of filter matrices)"
    filter_matrices::FM

    function MassPreservingCutoffFilter(grid, Nc = polynomialorders(grid))
        dim = dimensionality(grid)

        # Support different filtering thresholds in different
        # directions (default behavior is to apply the same threshold
        # uniformly in all directions)
        if Nc isa Integer
            Nc = ntuple(i -> Nc, dim)
        elseif Nc isa NTuple{2} && dim == 3
            Nc = (Nc[1], Nc[1], Nc[2])
        end
        @assert length(Nc) == dim

        # Tuple of polynomial degrees (N₁, N₂, N₃)
        N = polynomialorders(grid)
        # In 2D, we assume same polynomial order in the horizontal
        @assert dim == 2 || N[1] == N[2]
        @assert all(0 .<= Nc)

        σ(η) = 0

        AT = arraytype(grid)
        ξ = referencepoints(grid)
        filter_matrices = ntuple(i -> AT(modified_filter_matrix(ξ[i], Nc[i], σ)), dim)
        new{typeof(filter_matrices)}(filter_matrices)
    end
end

function apply_async!(
    Q,
    target::AbstractFilterTarget,
    grid::DiscontinuousSpectralElementGrid,
    filter::MassPreservingCutoffFilter;
    dependencies,
    state_auxiliary = nothing,
    direction = EveryDirection(),
)
    topology = grid.topology

    device = typeof(Q.data) <: Array ? CPU() : CUDADevice()

    dim = dimensionality(grid)
    N = polynomialorders(grid)
    # Currently only support same polynomial in both horizontal directions
    @assert dim == 2 || N[1] == N[2]
    Nq = N .+ 1
    Nq1 = Nq[1]
    Nq2 = Nq[2]
    Nq3 = dim == 2 ? 1 : Nq[dim]

    nrealelem = length(topology.realelems)
    # parallel sum info
    nreduce = 2^ceil(Int, log2(Nq1 * Nq2 * Nq3))
    event = dependencies

    if direction isa EveryDirection || direction isa HorizontalDirection
        @assert dim == 2 || Nq1 == Nq2
        filtermatrix = filter.filter_matrices[1]
        event = kernel_apply_mp_filter!(device, (Nq1, Nq2, Nq3))(
            Val(nreduce),
            Val(dim),
            Val(N),
            Val(vars(Q)),
            Val(isnothing(state_auxiliary) ? nothing : vars(state_auxiliary)),
            HorizontalDirection(),
            Q.data,
            isnothing(state_auxiliary) ? nothing : state_auxiliary.data,
            target,
            filtermatrix,
            grid.vgeo,
            ndrange = (nrealelem * Nq1, Nq2, Nq3),
            dependencies = event,
        )
    end
    if direction isa EveryDirection || direction isa VerticalDirection
        filtermatrix = filter.filter_matrices[end]
        event = kernel_apply_mp_filter!(device, (Nq1, Nq2, Nq3))(
            Val(nreduce),
            Val(dim),
            Val(N),
            Val(vars(Q)),
            Val(isnothing(state_auxiliary) ? nothing : vars(state_auxiliary)),
            VerticalDirection(),
            Q.data,
            isnothing(state_auxiliary) ? nothing : state_auxiliary.data,
            target,
            filtermatrix,
            grid.vgeo,
            ndrange = (nrealelem * Nq1, Nq2, Nq3),
            dependencies = event,
        )
    end
    return event
end

const _M = Grids._M
# mp for masspreserving
"""
    kernel_apply_mp_filter!(::Val{dim}, ::Val{N}, direction,
                         Q, state_auxiliary, target, filtermatrix
                        ) where {dim, N}
Computational kernel: Applies the `filtermatrix` to `Q` given a
custom target `target`.
The `direction` argument is used to control if the filter is applied in the
"""
@kernel function kernel_apply_mp_filter!(
    ::Val{nreduce},
    ::Val{dim},
    ::Val{N},
    ::Val{vars_Q},
    ::Val{vars_state_auxiliary},
    direction,
    Q,
    state_auxiliary,
    target::AbstractFilterTarget,
    filtermatrix,
    vgeo,
) where {nreduce, dim, N, vars_Q, vars_state_auxiliary}
    @uniform begin
        FT = eltype(Q)

        Nqs = N .+ 1
        Nq1 = Nqs[1]
        Nq2 = Nqs[2]
        Nq3 = dim == 2 ? 1 : Nqs[dim]

        if direction isa EveryDirection
            filterinξ1 = filterinξ2 = true
            filterinξ3 = dim == 2 ? false : true
        elseif direction isa HorizontalDirection
            filterinξ1 = true
            filterinξ2 = dim == 2 ? false : true
            filterinξ3 = false
        elseif direction isa VerticalDirection
            filterinξ1 = false
            filterinξ2 = dim == 2 ? true : false
            filterinξ3 = dim == 2 ? false : true
        end

        nstates = varsize(vars_Q)
        nfilterstates = number_state_filtered(target, FT)
        nfilteraux = isnothing(state_auxiliary) ? 0 : varsize(vars_state_auxiliary)

        # ugly workaround around problems with @private
        # hopefully will be soon fixed in KA
        l_Q2 = MVector{nstates, FT}(undef)
        l_Qfiltered2 = MVector{nfilterstates, FT}(undef)
    end

    s_Q = @localmem FT (Nq1, Nq2, Nq3, nfilterstates) # element local
    s_MQᴮ = @localmem FT (Nq1 * Nq2 * Nq3, nstates) # before applying filter
    s_MQᴬ = @localmem FT (Nq1 * Nq2 * Nq3, nstates) # after applying filter
    s_M = @localmem FT (Nq1 * Nq2 * Nq3) # local mass matrix

    l_Q = @private FT (nstates,)
    l_Qfiltered = @private FT (nfilterstates,) # scratch space for storing mat mul
    l_aux = @private FT (nfilteraux,)

    e = @index(Group, Linear)
    i, j, k = @index(Local, NTuple)
    ijk = @index(Local, Linear)

    @inbounds begin

        @unroll for s in 1:nstates
            l_Q[s] = Q[ijk, s, e]
        end

        @unroll for s in 1:nfilteraux
            l_aux[s] = state_auxiliary[ijk, s, e]
        end

        # Load mass weighted quantities to shared memory
        s_M[ijk] = vgeo[ijk, _M, e]
        @unroll for s in 1:nstates
            s_MQᴮ[ijk, s] = s_M[ijk] * l_Q[s]
        end

        fill!(l_Qfiltered2, -zero(FT))

        compute_filter_argument!(
            target,
            Vars{vars_state_filtered(target, FT)}(l_Qfiltered2),
            Vars{vars_Q}(l_Q[:]),
            Vars{vars_state_auxiliary}(l_aux[:]),
        )

        @unroll for fs in 1:nfilterstates
            l_Qfiltered[fs] = zero(FT)
        end

        @unroll for fs in 1:nfilterstates
            s_Q[i, j, k, fs] = l_Qfiltered2[fs]
        end

        if filterinξ1
            @synchronize
            @unroll for n in 1:Nq1
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += filtermatrix[i, n] * s_Q[n, j, k, fs]
                end
            end

            if filterinξ2 || filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    s_Q[i, j, k, fs] = l_Qfiltered[fs]
                    l_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ2
            @synchronize
            @unroll for n in 1:Nq2
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += filtermatrix[j, n] * s_Q[i, n, k, fs]
                end
            end

            if filterinξ3
                @synchronize
                @unroll for fs in 1:nfilterstates
                    s_Q[i, j, k, fs] = l_Qfiltered[fs]
                    l_Qfiltered[fs] = zero(FT)
                end
            end
        end

        if filterinξ3
            @synchronize
            @unroll for n in 1:Nq3
                @unroll for fs in 1:nfilterstates
                    l_Qfiltered[fs] += filtermatrix[k, n] * s_Q[i, j, n, fs]
                end
            end
        end

        @unroll for s in 1:nstates
            l_Q2[s] = l_Q[s]
        end

        compute_filter_result!(
            target,
            Vars{vars_Q}(l_Q2),
            Vars{vars_state_filtered(target, FT)}(l_Qfiltered[:]),
            Vars{vars_state_auxiliary}(l_aux[:]),
        )
        # Store result
        @unroll for s in 1:nstates
            l_Q[s] = l_Q2[s]
            s_MQᴬ[ijk, s] = s_M[ijk] * l_Q[s]
        end

        @synchronize
        @unroll for n in 11:-1:1
            if nreduce ≥ (1 << n)
                ijkshift = ijk + (1 << (n - 1))
                if ijk ≤ (1 << (n - 1)) && ijkshift ≤ Nq1 * Nq2 * Nq3
                    s_M[ijk] += s_M[ijkshift]
                    @unroll for s in 1:nstates
                        s_MQᴮ[ijk, s] += s_MQᴮ[ijkshift, s]
                        s_MQᴬ[ijk, s] += s_MQᴬ[ijkshift, s]
                    end
                end
                @synchronize
            end
        end

        @synchronize
        M⁻¹ = FT(1) / s_M[1]
        @unroll for s in 1:nstates
            Q[ijk, s, e] = l_Q[s] + M⁻¹ * (s_MQᴮ[1, s] - s_MQᴬ[1, s])
        end

        @synchronize
    end

end


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test

@testset "Wigner Chains W^3j_k" begin
    @testset "k = 1" begin
        @testset "ell = $ell" for ell=0:10
            w3j_1_slow = SFB.WignerChains.calc_w3j_k([ell], [0], [0])
            w3j_1_fast = SFB.WignerChains.calc_w3j_1([ell], [0], [0])
            @test w3j_1_slow ≈ (-1)^ell * √(2*ell+1) rtol=eps(10.0)
            @test w3j_1_fast ≈ (-1)^ell * √(2*ell+1) rtol=eps(10.0)

            for L=1:10, M=-L:L
                w3j_1_slow = SFB.WignerChains.calc_w3j_k([ell], [L], [M])
                w3j_1_fast = SFB.WignerChains.calc_w3j_1([ell], [L], [M])
                @test w3j_1_slow ≈ 0 atol=eps(1.0)
                @test w3j_1_fast ≈ 0 atol=eps(1.0)
            end
        end
        k = 1
        @show k
    end

    @testset "k = 2" begin
        lmax = 5
        SFB.WignerChains.calc_w3j_k([0,0], [0,0], [0,0]) # compile
        SFB.WignerChains.calc_w3j_2([0,0], [0,0], [0,0]) # compile
        SFB.WignerChains.calc_w3j_2_simple([0,0], [0,0], [0,0]) # compile
        tslow = 0.0
        tfast = 0.0
        tsimple = 0.0
        @time for l1=0:lmax, l2=0:lmax, L1=0:lmax, L2=0:lmax, M1=-L1:L1, M2=-L2:L2
            ell = [l1, l2]
            L = [L1, L2]
            M = [M1, M2]
            #@show ell,L,M
            tslow += @elapsed w3j_2_slow = SFB.WignerChains.calc_w3j_k(ell, L, M)
            tfast += @elapsed w3j_2_fast = SFB.WignerChains.calc_w3j_2(ell, L, M)
            tsimple += @elapsed w3j_2_simple = SFB.WignerChains.calc_w3j_2_simple(ell, L, M)
            @test w3j_2_fast ≈ w3j_2_simple atol=eps(10.0)
            @test w3j_2_slow ≈ w3j_2_simple atol=eps(10.0)
        end
        k = 2
        @show k tsimple tslow tfast
    end

    @testset "k = 3" begin
        lmax = 3
        SFB.WignerChains.calc_w3j_k([0,1,0], [0,0,0], [0,0,0]) # compile
        SFB.WignerChains.calc_w3j_3([0,1,0], [0,0,0], [0,0,0]) # compile
        tslow = 0.0
        tfast = 0.0
        @time for l1=0:lmax, l2=0:lmax, l3=0:lmax, L1=0:lmax, L2=0:lmax, L3=0:lmax
            for M1=-L1:L1, M2=-L2:L2, M3=-L3:L3
                ell = [l1, l2, l3]
                L = [L1, L2, L3]
                M = [M1, M2, M3]
                #@show ell,L,M
                tslow += @elapsed w3j_3_slow = SFB.WignerChains.calc_w3j_k(ell, L, M)
                tfast += @elapsed w3j_3_fast = SFB.WignerChains.calc_w3j_3(ell, L, M)
                @test w3j_3_slow ≈ w3j_3_fast atol=eps(10.0)
            end
        end
        k = 3
        @show k tslow tfast
    end

    @testset "k = 4" begin
        lmax = 2
        SFB.WignerChains.calc_w3j_k([0,1,0,0], [0,0,0,0], [0,0,0,0]) # compile
        SFB.WignerChains.calc_w3j_4([0,1,0,0], [0,0,0,0], [0,0,0,0]) # compile
        tslow = 0.0
        tfast = 0.0
        @time for l1=0:lmax, l2=0:lmax, l3=0:lmax, l4=0:lmax
            for L1=0:lmax, L2=0:lmax, L3=0:lmax, L4=0:lmax
                for M1=-L1:L1, M2=-L2:L2, M3=-L3:L3, M4=-L4:L4
                    ell = [l1, l2, l3, l4]
                    L = [L1, L2, L3, L4]
                    M = [M1, M2, M3, M4]
                    #@show ell,L,M
                    tslow += @elapsed w3j_4_slow = SFB.WignerChains.calc_w3j_k(ell, L, M)
                    tfast += @elapsed w3j_4_fast = SFB.WignerChains.calc_w3j_4(ell, L, M)
                    @test w3j_4_slow ≈ w3j_4_fast atol=eps(10.0)
                end
            end
        end
        k = 4
        @show k tslow tfast
    end

    @testset "k = 5" begin
        lmax = 1
        SFB.WignerChains.calc_w3j_k([0,1,0,0,0], [0,0,0,0,0], [0,0,0,0,0]) # compile
        SFB.WignerChains.calc_w3j_5([0,1,0,0,0], [0,0,0,0,0], [0,0,0,0,0]) # compile
        tslow = 0.0
        tfast = 0.0
        @time for l1=0:lmax, l2=0:lmax, l3=0:lmax, l4=0:lmax, l5=0:lmax
            for L1=0:lmax, L2=0:lmax, L3=0:lmax, L4=0:lmax, L5=0:lmax
                for M1=-L1:L1, M2=-L2:L2, M3=-L3:L3, M4=-L4:L4, M5=-L5:L5
                    ell = [l1, l2, l3, l4, l5]
                    L = [L1, L2, L3, L4, L5]
                    M = [M1, M2, M3, M4, M5]
                    #@show ell,L,M
                    tslow += @elapsed w3j_5_slow = SFB.WignerChains.calc_w3j_k(ell, L, M)
                    tfast += @elapsed w3j_5_fast = SFB.WignerChains.calc_w3j_5(ell, L, M)
                    @test w3j_5_slow ≈ w3j_5_fast atol=eps(10.0)
                end
            end
        end
        k = 5
        @show k tslow tfast
    end
end

# vim: set sw=4 et sts=4 :


text:
module InferenceUtilities
using Base.Test

export @isinferred_noneval, @isinferred, @return_types

"""
    @isinferred_noneval f(x)

Checks whether the call expression `f(x)` is type stable without evaluating `f`
(but it does evaluate the argument(s) x) and returns a boolean value.
`@isinferred` is similar to `@inferred`, but does not evaluate `f` and returns
a boolean rather than throwing an exception. It also does not capture some
corner cases that `@inferred` (and hence `@isinferred`, which is built on
`@inferred`, captures).
This macro is useful for unit tests (`@test @isinferred_noneval f(x)`).

`f(x)` can be any call expression.

```jldoctest
julia> using Base.Test

julia> f(a,b,c) = b > 1 ? 1 : 1.0
f (generic function with 1 method)

julia> @code_warntype f(1,2,3)
Variables:
  #self#::#f
  a::Int64
  b::Int64
  c::Int64

Body:
  begin
      unless (Base.slt_int)(1, b::Int64)::Bool goto 3
      return 1
      3:
      return 1.0
  end::UNION{FLOAT64, INT64}

julia> @isinferred_noneval f(1,2,3)
false

julia> @isinferred_noneval max(1,2)
false
```
"""
macro isinferred_noneval(ex)
    if Meta.isexpr(ex, :ref)
        ex = Expr(:call, :getindex, ex.args...)
    end
    Meta.isexpr(ex, :call)|| error("@inferred requires a call expression")

    quote
        let
            args = ($([esc(ex.args[i]) for i = 2:length(ex.args)]...),)
            inftypes = Base.return_types($(esc(ex.args[1])), Base.typesof(args...))
            length(inftypes) == 1 && isleaftype(first(inftypes))
        end
    end
end

"""
    @isinferred f(x)
Checks whether the call expression `f(x)` is type stable and returns a
boolean value. `@isinferred` is similar to `@inferred`,
but returns a boolean rather than throwing an exception.
This macro is useful for unit tests (`@test @isinferred f(x)`).
`f(x)` can be any call expression.
```jldoctest
julia> using Base.Test
julia> f(a,b,c) = b > 1 ? 1 : 1.0
f (generic function with 1 method)
julia> @code_warntype f(1,2,3)
Variables:
  #self#::#f
  a::Int64
  b::Int64
  c::Int64
Body:
  begin
      unless (Base.slt_int)(1, b::Int64)::Bool goto 3
      return 1
      3:
      return 1.0
  end::UNION{FLOAT64, INT64}
julia> @isinferred f(1,2,3)
false
julia> @isinferred max(1,2)
false
```
"""
macro isinferred(ex)
    quote
        try
            @inferred $(esc(ex))
            true
        catch err
            isa(err, ErrorException) ? false : rethrow(err)
        end
    end
end

"""
  @return_types f(x)

  Evaluates arguments `x` and returns the inferred return type of `f(x)` without
  evaluating `f(x)`

```jldoctests
julia> f(a,b,c) = b > 1 ? 1 : 1.0
f (generic function with 1 method)

julia> @return_types f(1,2,3)
1-element Array{Any,1}:
 Union{Float64, Int64}
```
"""
macro return_types(ex)
    f = ex.args[1]
    f_args = [:(typeof($i)) for i in ex.args[2:end]]
    :(Base.return_types($(esc(f)), Base.typesof($(f_args...))))
end
end # module


text:
if VERSION >= v"1.3" && !(haskey(ENV, "JULIA_IPOPT_LIBRARY_PATH") && haskey(ENV, "JULIA_IPOPT_EXECUTABLE_PATH"))
    exit()  # Use Ipopt_jll instead.
end

using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libipopt"], :libipopt),
    ExecutableProduct(prefix, "ipopt", :amplexe),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/Ipopt_jll.jl/releases/download/Ipopt-v3.13.2+0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-gnu-libgfortran3-cxx03.tar.gz", "bed359bbf63128215d5adb1dc53044ed97566cad757f72b440275932e74fcd2c"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-gnu-libgfortran3-cxx11.tar.gz", "a236e3fdfae8c397fc2400c60d57f31f57541b6f22650a8df9b8a0f7b198a7bb"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-gnu-libgfortran4-cxx03.tar.gz", "3c3507a9b7f7c54881d720f5dd20f39c871407d0196055cdcffbc7f4059f0cd7"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-gnu-libgfortran4-cxx11.tar.gz", "9b0f28c8987b6f2c5c39122a37e9b5accda4019edeacb3be3b8ebe0de320e693"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-gnu-libgfortran5-cxx03.tar.gz", "aa06b89cc0aa4b9157a99260f0adb0a9e85a18c20c5cac4b781498d6ed2d1234"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-gnu-libgfortran5-cxx11.tar.gz", "f6e18d9fb562aa44ed9ddbab1cb2122fcd2cb0dc034686fb5af8a58884c7554a"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-musl-libgfortran3-cxx03.tar.gz", "5c1bfc631ddd0e6de19223558d482aec733f53f99d7af26436753de2a94090cd"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-musl-libgfortran3-cxx11.tar.gz", "e9dcad5a9e8d72cd1250168ffbbb5c68759869ba84e057b8877abf9146940f07"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-musl-libgfortran4-cxx03.tar.gz", "0bb93a5db8fd0b9aa0bdc41db68025bfdebe8ef468a4c6b88b83c6cb4ad3e5a1"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-musl-libgfortran4-cxx11.tar.gz", "65373153d5f8d14504b1e7d0eaf967baf8408db3e749e6f342dda7d755468847"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-musl-libgfortran5-cxx03.tar.gz", "329dd313ba6da542c0ae06a236402f0279e0c45f4d8a9aa1373a8bf82a869c4e"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.aarch64-linux-musl-libgfortran5-cxx11.tar.gz", "6f36bc0fccbd6292f764108c218b26e4449e891a8f063d4b7958700005424f55"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-gnueabihf-libgfortran3-cxx03.tar.gz", "269d9d86c8d23285f864fec19b65d8f6c618ab2f217ea23971f527b3b776a450"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-gnueabihf-libgfortran3-cxx11.tar.gz", "ed2aa1c432f53adf0e5c64524e1b401fbb65a99e4a3d16d9561484f0e5d385bf"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-gnueabihf-libgfortran4-cxx03.tar.gz", "2217cd2ed481ba19de2ce49b13b7f85b0fc87151f207b4846806dc7d076157e5"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-gnueabihf-libgfortran4-cxx11.tar.gz", "e83987f3a60cf400411bdcd295a8fe18398a77ce160e6a09ec92eeb79089e4da"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-gnueabihf-libgfortran5-cxx03.tar.gz", "a9cc1439fc05657ba2dbf6b5f2b548e13364e00940d4f7fd0aad0835a0a1bf91"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-gnueabihf-libgfortran5-cxx11.tar.gz", "6add6a40fb8c6aeee961b7abadf2a3749a821ec68aa5c52cbc0fcede97865da8"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-musleabihf-libgfortran3-cxx03.tar.gz", "e7f93eb56c0408e21381a2efebf7d197e1cd8b423e5768159a541fd386d4c3fa"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-musleabihf-libgfortran3-cxx11.tar.gz", "1d9055d6db6bb0899201506ad885fd6398de2c541489b84e6f34a42abd8475d9"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-musleabihf-libgfortran4-cxx03.tar.gz", "68cce661e19725c087edce05c20ddcf9e78ad697df63af3a2088660a57ebcac9"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-musleabihf-libgfortran4-cxx11.tar.gz", "86fefe99701772d0b45a4759f1382a49c80963cbcd47266e7fe85fc8f7b0632d"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-musleabihf-libgfortran5-cxx03.tar.gz", "398122631c32ccd584bb1761f3f8be9630b3d9f4aa9649d2af503f2e7055b756"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.armv7l-linux-musleabihf-libgfortran5-cxx11.tar.gz", "aa5166ffcd1ad94783beccdc3bb171321974268ef813d65192851a6a2c82f656"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-gnu-libgfortran3-cxx03.tar.gz", "f597e24bc521fee19bd9547347a859b4201deb0c0ca1bc457e7c41cb113fbff8"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-gnu-libgfortran3-cxx11.tar.gz", "ad5278604795f5990f11e27fd1b3927094599c2b458e97aec92e93cd240139b7"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-gnu-libgfortran4-cxx03.tar.gz", "43ac1f94e559016fc117404aa0f4d33cada2cc5a15a308564eb20c06d8c6c3c7"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-gnu-libgfortran4-cxx11.tar.gz", "5c1527018aa2039f6f4251b0eb3c9287aead603f521f9a9563d6baa2807c5eca"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-gnu-libgfortran5-cxx03.tar.gz", "3c352a9b6772ed9d9745efd5bf61a7ddb8e742609db8aa5509cc8ffce42c1c06"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-gnu-libgfortran5-cxx11.tar.gz", "694d7a0191aca91c7d04fc1459cd4d170b051b70bc6c0bb966347e0325d08147"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-musl-libgfortran3-cxx03.tar.gz", "d5c08f899e90e357f69d886731334f8e57e162d66d8bddb4d54f25b4a86f7a63"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-musl-libgfortran3-cxx11.tar.gz", "a88d8ea6d26065b9b478ffe7d05013a043bb15ba8ecd0f08972b78caed6a4e1d"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-musl-libgfortran4-cxx03.tar.gz", "ed0b8242669dec58f797b5647f21fb045000b85edf0ddd48ee08a0a88adf858a"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-musl-libgfortran4-cxx11.tar.gz", "d357ebe823f76e29b6380b72818a0a9597c0e1fe9a3a5a9be752ba9daf5255ae"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-musl-libgfortran5-cxx03.tar.gz", "2baea86a2b9343fe8257accd6fb7cf563a79c971e2263d76c98a3a0186de737a"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-linux-musl-libgfortran5-cxx11.tar.gz", "07197d3aee44396bd289030d7b9ef31399d8959ebc6ff16cb58cc9ef15e12d8d"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-w64-mingw32-libgfortran3-cxx03.tar.gz", "d67d2be408fa63ee67fb8b443f3a52c905d7295284d7c7d6c7fb30bbb650035e"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-w64-mingw32-libgfortran3-cxx11.tar.gz", "7417c88af4cd026d967e44447035ba3c248ae730ab590a36a55a6004b7fc4fdf"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-w64-mingw32-libgfortran4-cxx03.tar.gz", "1bbe178ea16c014e6294bb60c3a88027c4a4511bb30eb57ce54a3bd5ae0384f6"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-w64-mingw32-libgfortran4-cxx11.tar.gz", "4ce9232e7f902a8b9999f3cfa0c934600110176a40a55d6c3921814bf64a35da"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.i686-w64-mingw32-libgfortran5-cxx03.tar.gz", "69fd33000a7085e3713cf7a50f2feca1ec7a0351567f89309cf05b06c9cea7b0"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.i686-w64-mingw32-libgfortran5-cxx11.tar.gz", "92672ea65b1599e47ef77eab638b497e383d52d4447b2037572e598fd59229f7"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.powerpc64le-linux-gnu-libgfortran3-cxx03.tar.gz", "f33dd8f30eaf99ad6849b04988dabd4796ebd92c73866d2c1d1a99647aa54330"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.powerpc64le-linux-gnu-libgfortran3-cxx11.tar.gz", "bd82361e489b530c388fe647d541d8debdc1a050d92c4796f82ba2db221e24a4"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.powerpc64le-linux-gnu-libgfortran4-cxx03.tar.gz", "53c9f9de27a94d3240294ea16c1ca39b72d34bb246cfb17acbff4a93ba9def40"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.powerpc64le-linux-gnu-libgfortran4-cxx11.tar.gz", "e2a057666807d274643c3f9eaa4d3ce867a4e6a74537f40afd458af2d7c44d06"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.powerpc64le-linux-gnu-libgfortran5-cxx03.tar.gz", "67a5fcbc1cf469f35237457537ab2d6878544b6e10e4467a64165c6709a98536"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.powerpc64le-linux-gnu-libgfortran5-cxx11.tar.gz", "8c8c9a6ee63e849c783d7f0d034f79700ddb44ce30da3b7c7b71096fa3cc4c4a"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-apple-darwin14-libgfortran3-cxx03.tar.gz", "a905edd95dbbc38a766408215df90607c5de1149713e54829134c8f205fe189b"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-apple-darwin14-libgfortran3-cxx11.tar.gz", "6b1954a890099cf07d76ad403a89535f80583d58a6be891097dfd57d3a2199b1"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-apple-darwin14-libgfortran4-cxx03.tar.gz", "57ec817ca75ff6befc5b016bfe4796af3cc6d74acbda491d6f4a822ae32e2616"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-apple-darwin14-libgfortran4-cxx11.tar.gz", "2223d9702ead939291ccebd37b65974561c49ba1f08cb826294c2bfd820eca63"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-apple-darwin14-libgfortran5-cxx03.tar.gz", "6c844ed7212af30073964b57400241e24de99ca2262e85090f5d34a21a9f972f"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-apple-darwin14-libgfortran5-cxx11.tar.gz", "72c3d4ae1acdd4815129f329e660008c1767a291847605c09a6d87fe1277cb58"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-gnu-libgfortran3-cxx03.tar.gz", "d496a355c79221228b69dac865567bda121d75772e06f9fbe5bd0b0ab13ceab1"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-gnu-libgfortran3-cxx11.tar.gz", "c9199f5d4bf79121208d3a48f216bbe643dccdac6ee31c7ceff9b727b58c71b9"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-gnu-libgfortran4-cxx03.tar.gz", "4837bc4de9dbb328c171938eceba24f44fd6989a434f9981918348b3c42b8aeb"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-gnu-libgfortran4-cxx11.tar.gz", "4583675e767cd3dee64ecc3cd54a65b943a1769972d805fd2eb1dc4b3cda6bf3"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-gnu-libgfortran5-cxx03.tar.gz", "f46bf6a3a42a4cb175feb73102b2ec62ede908a05d5c2e871d1f8e1dbf4aceec"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-gnu-libgfortran5-cxx11.tar.gz", "be3577870b855476c7b22f134950e98ee27b9bcd73ab456a2d626c8f6f95cc76"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-musl-libgfortran3-cxx03.tar.gz", "a039c6ea305705daada0fffea56ba12f089cd773ab8973fb690ea35d2cba34fb"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-musl-libgfortran3-cxx11.tar.gz", "bd4cae26147fd6d870d81ef052894c19db1f85094e22cf5bdf7842dc62d3585f"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-musl-libgfortran4-cxx03.tar.gz", "1a2de72395c925ea03f9f6d1f84f3a526d6c66df3d8b361a6cee41744190278f"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-musl-libgfortran4-cxx11.tar.gz", "c9ee7a579cca559554f7503fa797d9c08008af4a417945c29ec1fe73be889bf1"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-musl-libgfortran5-cxx03.tar.gz", "89ed817d417d5552a77618efd28e8efc6a3a27dfa75c2dfd14cdd325b331109a"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-linux-musl-libgfortran5-cxx11.tar.gz", "b64ccf90b27f73cb4c0831817176f80520be69c5dcd6da9f06093097cebfb4d7"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-unknown-freebsd11.1-libgfortran3-cxx03.tar.gz", "bcdaee0b1a7dd796f7bcf130740e7d62aba6da3bbaae93028789279ec3215d5c"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-unknown-freebsd11.1-libgfortran3-cxx11.tar.gz", "165f086abd60da2b4ac1f08b02958c69ad1c64880b45441c58081bd0f9e5abbf"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-unknown-freebsd11.1-libgfortran4-cxx03.tar.gz", "24899c094690436654e54f1648cdb637c7699d64f268ef843c3cf74e91cddf2b"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-unknown-freebsd11.1-libgfortran4-cxx11.tar.gz", "c3ef1366203443ed9079d51b8932f4c5b38b0fac2c1f6cc1d1c2c3276670b2e9"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-unknown-freebsd11.1-libgfortran5-cxx03.tar.gz", "2a4d281f72f6b9ea67456daed3920273e06ba1e80f51040ebf76a1e4497b76ea"),
    FreeBSD(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-unknown-freebsd11.1-libgfortran5-cxx11.tar.gz", "9cc2d63dfd3ae768aab4b20c9e2641c6804dc9a2b5d27f41908b3c9370b7d4a3"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-w64-mingw32-libgfortran3-cxx03.tar.gz", "70cae55c61a6bd8eaabec96c5c2085951862cc5052b23cff8f87ae05080a032b"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc4, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-w64-mingw32-libgfortran3-cxx11.tar.gz", "f01e4120877ac722514a7eded94a351477e2f2f29f5f9b4a94d4c24bd7aa139e"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-w64-mingw32-libgfortran4-cxx03.tar.gz", "a837bda68c6059dff057a1968743c9fbc3957ce8568ad226d7517c364646682f"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-w64-mingw32-libgfortran4-cxx11.tar.gz", "5fe62d97efd77b35ed7ac428a102c7201a2d859c6852547944380a5eef9b26ab"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx03)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-w64-mingw32-libgfortran5-cxx03.tar.gz", "fe339c6d2c1371267ba8f7519de6afa665ff77a1337bdd7566cfdc373464436d"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8, :cxx11)) => ("$bin_prefix/Ipopt.v3.13.2.x86_64-w64-mingw32-libgfortran5-cxx11.tar.gz", "c173c91afe765288658ae9eb1e55e2ebd72604c63274a2aa9c3a552f617ddf2e"),
)

# To fix gcc4 bug in Windows
# https://sourceforge.net/p/mingw-w64/bugs/727/
this_platform = platform_key_abi()
if typeof(this_platform)==Windows && this_platform.compiler_abi.gcc_version == :gcc4
   this_platform = Windows(arch(this_platform), libc=libc(this_platform), compiler_abi=CompilerABI(:gcc6))
end

function update_product(product::LibraryProduct, library_path, binary_path)
    product = LibraryProduct(library_path, product.libnames, product.variable_name)
    # compatible with dll generate by coinbrew & msvc
    if ~satisfied(product)
        product = LibraryProduct(library_path, "ipopt", product.variable_name)
    end
    return product
end

function update_product(product::ExecutableProduct, library_path, binary_path)
    ExecutableProduct(joinpath(binary_path, basename(product.path)), product.variable_name)
end

custom_library = false
if haskey(ENV,"JULIA_IPOPT_LIBRARY_PATH") && haskey(ENV,"JULIA_IPOPT_EXECUTABLE_PATH")
    custom_products = [update_product(product, ENV["JULIA_IPOPT_LIBRARY_PATH"], ENV["JULIA_IPOPT_EXECUTABLE_PATH"]) for product in products]
    if all(satisfied(p; verbose=verbose) for p in custom_products)
        products = custom_products
        custom_library = true
    else
        error("Could not install custom libraries from $(ENV["JULIA_IPOPT_LIBRARY_PATH"]) and $(ENV["JULIA_IPOPT_EXECUTABLE_PATH"]).\nTo fall back to BinaryProvider call delete!(ENV,\"JULIA_IPOPT_LIBRARY_PATH\");delete!(ENV,\"JULIA_IPOPT_EXECUTABLE_PATH\") and run build again.")
    end
end

if !custom_library
    # Install unsatisfied or updated dependencies:
    # We added `, isolate=true` as otherwise, it would segfault when closing `OpenBLAS32`,
    # probably because it is conflicting with Julia openblas.
    unsatisfied = any(!satisfied(p; verbose=verbose, isolate=true) for p in products)

    dl_info = choose_download(download_info, platform_key_abi())
    if dl_info === nothing && unsatisfied
        # If we don't have a compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something even more ambitious here.
        error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
    end

    # If we have a download, and we are unsatisfied (or the version we're
    # trying to install is not itself installed) then load it up!
    if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
        # Download and install binaries
        evalfile("build_CompilerSupportLibraries.v0.3.3.jl")
        evalfile("build_ASL.v0.1.1.jl")
        evalfile("build_METIS.v5.1.0.jl")
        evalfile("build_OpenBLAS32.v0.3.9.jl")
        evalfile("build_MUMPS_seq.v5.2.1.jl")
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)