text:
"""
$(SIGNATURES)

Calculate moments of Gaussian distribution `G = (λ / π)^(D / 2) * exp[-λ(c^2 + ξ^2)]`
"""
function gauss_moments(prim::T) where {T<:AA{<:Real,1}}
    if eltype(prim) <: Int
        MuL = OffsetArray(similar(prim, Float64, 7), 0:6)
    else
        MuL = OffsetArray(similar(prim, 7), 0:6)
    end
    MuR = similar(MuL)
    Mu = similar(MuL)

    MuL[0] = 0.5 * SpecialFunctions.erfc(-sqrt(prim[end]) * prim[2])
    MuL[1] = prim[2] * MuL[0] + 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    MuR[0] = 0.5 * SpecialFunctions.erfc(sqrt(prim[end]) * prim[2])
    MuR[1] = prim[2] * MuR[0] - 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    for i = 2:6
        MuL[i] = prim[2] * MuL[i-1] + 0.5 * (i - 1) * MuL[i-2] / prim[end]
        MuR[i] = prim[2] * MuR[i-1] + 0.5 * (i - 1) * MuR[i-2] / prim[end]
    end
    @. Mu = MuL + MuR

    if length(prim) == 3

        return Mu, MuL, MuR

    elseif length(prim) == 4

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        return Mu, Mv, MuL, MuR

    elseif length(prim) == 5

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        Mw = similar(MuL)
        Mw[0] = 1.0
        Mw[1] = prim[4]
        for i = 2:6
            Mw[i] = prim[4] * Mw[i-1] + 0.5 * (i - 1) * Mw[i-2] / prim[end]
        end

        return Mu, Mv, Mw, MuL, MuR

    end
end

"""
$(SIGNATURES)

A more general function dealing with internal energy
"""
function gauss_moments(prim::T, inK) where {T<:AA{<:Real,1}}
    if eltype(prim) <: Int
        MuL = OffsetArray(similar(prim, Float64, 7), 0:6)
    else
        MuL = OffsetArray(similar(prim, 7), 0:6)
    end
    MuR = similar(MuL)
    Mu = similar(MuL)

    MuL[0] = 0.5 * SpecialFunctions.erfc(-sqrt(prim[end]) * prim[2])
    MuL[1] = prim[2] * MuL[0] + 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    MuR[0] = 0.5 * SpecialFunctions.erfc(sqrt(prim[end]) * prim[2])
    MuR[1] = prim[2] * MuR[0] - 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    for i = 2:6
        MuL[i] = prim[2] * MuL[i-1] + 0.5 * (i - 1) * MuL[i-2] / prim[end]
        MuR[i] = prim[2] * MuR[i-1] + 0.5 * (i - 1) * MuR[i-2] / prim[end]
    end
    @. Mu = MuL + MuR

    if length(prim) == 3

        Mxi = similar(MuL, 0:2)
        Mxi[0] = 1.0
        Mxi[1] = 0.5 * inK / prim[end]
        Mxi[2] = (inK^2 + 2.0 * inK) / (4.0 * prim[end]^2)

        return Mu, Mxi, MuL, MuR

    elseif length(prim) == 4

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        Mxi = similar(MuL, 0:2)
        Mxi[0] = 1.0
        Mxi[1] = 0.5 * inK / prim[end]
        Mxi[2] = (inK^2 + 2.0 * inK) / (4.0 * prim[end]^2)

        return Mu, Mv, Mxi, MuL, MuR

    elseif length(prim) == 5

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        Mw = similar(MuL)
        Mw[0] = 1.0
        Mw[1] = prim[4]
        for i = 2:6
            Mw[i] = prim[4] * Mw[i-1] + 0.5 * (i - 1) * Mw[i-2] / prim[end]
        end

        return Mu, Mv, Mw, MuL, MuR

    end
end

"""
$(SIGNATURES)

Calculate conservative moments of particle distribution
"""
moments_conserve(Mu::T, alpha::I) where {T<:OffsetArray{<:FN,1},I<:Integer} = Mu[alpha]

"""
$(SIGNATURES)
"""
function moments_conserve(
    Mu::T,
    Mxi::T,
    alpha::I,
    delta::I,
) where {T<:OffsetArray{<:FN,1},I<:Integer}

    uv = similar(Mu, 3)
    uv[1] = Mu[alpha] * Mxi[delta÷2]
    uv[2] = Mu[alpha+1] * Mxi[delta÷2]
    uv[3] = 0.5 * (Mu[alpha+2] * Mxi[delta÷2] + Mu[alpha] * Mxi[(delta+2)÷2])

    return uv

end

"""
$(SIGNATURES)
"""
function moments_conserve(
    Mu::T,
    Mv::T,
    Mw::T,
    alpha::I,
    beta::I,
    delta::I,
) where {T<:OffsetArray{<:FN,1},I<:Integer}

    if length(Mw) == 3 # internal motion
        uv = similar(Mu, 4)
        uv[1] = Mu[alpha] * Mv[beta] * Mw[delta÷2]
        uv[2] = Mu[alpha+1] * Mv[beta] * Mw[delta÷2]
        uv[3] = Mu[alpha] * Mv[beta+1] * Mw[delta÷2]
        uv[4] =
            0.5 * (
                Mu[alpha+2] * Mv[beta] * Mw[delta÷2] +
                Mu[alpha] * Mv[beta+2] * Mw[delta÷2] +
                Mu[alpha] * Mv[beta] * Mw[(delta+2)÷2]
            )
    else
        uv = similar(Mu, 5)
        uv[1] = Mu[alpha] * Mv[beta] * Mw[delta]
        uv[2] = Mu[alpha+1] * Mv[beta] * Mw[delta]
        uv[3] = Mu[alpha] * Mv[beta+1] * Mw[delta]
        uv[4] = Mu[alpha] * Mv[beta] * Mw[delta+1]
        uv[5] =
            0.5 * (
                Mu[alpha+2] * Mv[beta] * Mw[delta] +
                Mu[alpha] * Mv[beta+2] * Mw[delta] +
                Mu[alpha] * Mv[beta] * Mw[delta+2]
            )
    end

    return uv

end

"""
$(SIGNATURES)

Discrete moments of conservative variables

1F1V
"""
function moments_conserve(f::X, u::T, ω::T) where {X<:AA{<:FN,1},T<:AA{<:FN,1}}
    w = similar(f, 3)
    w[1] = discrete_moments(f, u, ω, 0)
    w[2] = discrete_moments(f, u, ω, 1)
    w[3] = 0.5 * discrete_moments(f, u, ω, 2)

    return w
end

"""
$(SIGNATURES)

2F1V
"""
function moments_conserve(h::X, b::X, u::T, ω::T) where {X<:AA{<:FN,1},T<:AA{<:FN,1}}
    w = similar(h, 3)
    w[1] = discrete_moments(h, u, ω, 0)
    w[2] = discrete_moments(h, u, ω, 1)
    w[3] = 0.5 * (discrete_moments(h, u, ω, 2) + discrete_moments(b, u, ω, 0))

    return w
end

"""
$(SIGNATURES)

1F2V
"""
function moments_conserve(f::X, u::T, v::T, ω::T) where {X<:AA{<:FN,2},T<:AA{<:FN,2}}
    w = similar(f, 4)
    w[1] = discrete_moments(f, u, ω, 0)
    w[2] = discrete_moments(f, u, ω, 1)
    w[3] = discrete_moments(f, v, ω, 1)
    w[4] = 0.5 * (discrete_moments(f, u, ω, 2) + discrete_moments(f, v, ω, 2))

    return w
end

"""
$(SIGNATURES)

2F2V
"""
function moments_conserve(h::X, b::X, u::T, v::T, ω::T) where {X<:AA{<:FN,2},T<:AA{<:FN,2}}
    w = similar(h, 4)
    w[1] = discrete_moments(h, u, ω, 0)
    w[2] = discrete_moments(h, u, ω, 1)
    w[3] = discrete_moments(h, v, ω, 1)
    w[4] =
        0.5 * (
            discrete_moments(h, u, ω, 2) +
            discrete_moments(h, v, ω, 2) +
            discrete_moments(b, u, ω, 0)
        )

    return w
end

"""
$(SIGNATURES)

3F2V
"""
function moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AA{<:FN,2},T<:AA{<:FN,2}}
    w = similar(h0, 5)
    w[1] = discrete_moments(h0, u, ω, 0)
    w[2] = discrete_moments(h0, u, ω, 1)
    w[3] = discrete_moments(h0, v, ω, 1)
    w[4] = discrete_moments(h1, u, ω, 0)
    w[5] =
        0.5 * (
            discrete_moments(h0, u, ω, 2) +
            discrete_moments(h0, v, ω, 2) +
            discrete_moments(h2, u, ω, 0)
        )

    return w
end

"""
$(SIGNATURES)

1F3V
"""
function moments_conserve(f::X, u::T, v::T, w::T, ω::T) where {X<:AA{<:FN,3},T<:AA{<:FN,3}}
    moments = similar(f, 5)

    moments[1] = discrete_moments(f, u, ω, 0)
    moments[2] = discrete_moments(f, u, ω, 1)
    moments[3] = discrete_moments(f, v, ω, 1)
    moments[4] = discrete_moments(f, w, ω, 1)
    moments[5] =
        0.5 * (
            discrete_moments(f, u, ω, 2) +
            discrete_moments(f, v, ω, 2) +
            discrete_moments(f, w, ω, 2)
        )

    return moments
end

"""
$(SIGNATURES)

4F1V
"""
function moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    h3::X,
    u::T,
    ω::T,
) where {X<:AA{<:FN,1},T<:AA{<:FN,1}}
    moments = similar(h0, 5)

    moments[1] = discrete_moments(h0, u, ω, 0)
    moments[2] = discrete_moments(h0, u, ω, 1)
    moments[3] = discrete_moments(h1, u, ω, 0)
    moments[4] = discrete_moments(h2, u, ω, 0)
    moments[5] = 0.5 * discrete_moments(h0, u, ω, 2) + 0.5 * discrete_moments(h3, u, ω, 0)

    return moments
end

"""
$(SIGNATURES)

Calculate conservative moments of diatomic particle distribution
"""
function diatomic_moments_conserve(
    h::X,
    b::X,
    r::X,
    u::T,
    ω::T,
) where {X<:AA{<:FN,1},T<:AA{<:FN,1}}
    w = similar(h, 4)
    w[1] = discrete_moments(h, u, ω, 0)
    w[2] = discrete_moments(h, u, ω, 1)
    w[3] =
        0.5 * (
            discrete_moments(h, u, ω, 2) +
            discrete_moments(b, u, ω, 0) +
            discrete_moments(r, u, ω, 0)
        )
    w[4] = 0.5 * discrete_moments(r, u, ω, 0)

    return w
end

"""
$(SIGNATURES)
"""
function diatomic_moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AA{<:FN,2},T<:AA{<:FN,2}}
    w = similar(h0, 5)
    w[1] = discrete_moments(h0, u, ω, 0)
    w[2] = discrete_moments(h0, u, ω, 1)
    w[3] = discrete_moments(h0, v, ω, 1)
    w[4] =
        0.5 * (
            discrete_moments(h0, u, ω, 2) +
            discrete_moments(h0, v, ω, 2) +
            discrete_moments(h1, u, ω, 0) +
            discrete_moments(h2, u, ω, 0)
        )
    w[5] = 0.5 * discrete_moments(h2, u, ω, 0)

    return w
end

"""
$(SIGNATURES)

Calculate slope-related conservative moments
`a = a1 + u * a2 + 0.5 * u^2 * a3`

"""
moments_conserve_slope(a, Mu::T, alpha::I) where {T<:OffsetArray{<:Real,1},I<:Int} =
    a * moments_conserve(Mu, alpha)

"""
$(SIGNATURES)
"""
moments_conserve_slope(
    a::X,
    Mu::Y,
    Mxi::Y,
    alpha::I,
) where {X<:AA{<:Real,1},Y<:OffsetArray{<:Real,1},I<:Int} =
    a[1] .* moments_conserve(Mu, Mxi, alpha + 0, 0) .+
    a[2] .* moments_conserve(Mu, Mxi, alpha + 1, 0) .+
    0.5 * a[3] .* moments_conserve(Mu, Mxi, alpha + 2, 0) .+
    0.5 * a[3] .* moments_conserve(Mu, Mxi, alpha + 0, 2)

"""
$(SIGNATURES)
"""
function moments_conserve_slope(
    a::X,
    Mu::Y,
    Mv::Y,
    Mw::Y,
    alpha::I,
    beta::I,
    delta = 0::I,
) where {X<:AA{<:Real,1},Y<:OffsetArray{<:Real,1},I<:Integer}

    if length(a) == 4
        return a[1] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, 0) .+
               a[2] .* moments_conserve(Mu, Mv, Mw, alpha + 1, beta + 0, 0) .+
               a[3] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 1, 0) .+
               0.5 * a[4] .* moments_conserve(Mu, Mv, Mw, alpha + 2, beta + 0, 0) .+
               0.5 * a[4] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 2, 0) .+
               0.5 * a[4] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, 2)
    elseif length(a) == 5
        return a[1] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, delta + 0) .+
               a[2] .* moments_conserve(Mu, Mv, Mw, alpha + 1, beta + 0, delta + 0) .+
               a[3] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 1, delta + 0) .+
               a[4] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, delta + 1) .+
               0.5 * a[5] .* moments_conserve(Mu, Mv, Mw, alpha + 2, beta + 0, delta + 0) .+
               0.5 * a[5] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 2, delta + 0) .+
               0.5 * a[5] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, delta + 2)
    end

end

"""
$(SIGNATURES)

Calculate conservative moments from microscopic moments
"""
function flux_conserve!(fw, args...)
    if length(fw) == 3
        return flux_conserve_1d!(fw, args...)
    elseif length(fw) == 4
        return flux_conserve_2d!(fw, args...)
    elseif length(fw) == 5
        return flux_conserve_3d!(fw, args...)
    end

    return nothing
end

#--- 1f1v ---#
function flux_conserve_1d!(fw, ff, u, ω)
    fw[1] = sum(ω .* ff)
    fw[2] = sum(u .* ω .* ff)
    fw[end] = 0.5 * sum(u .^ 2 .* ω .* ff)
end

#--- 2f1v ---#
function flux_conserve_1d!(fw, fh, fb, u, ω)
    fw[1] = sum(ω .* fh)
    fw[2] = sum(u .* ω .* fh)
    fw[end] = 0.5 * (sum(u .^ 2 .* ω .* fh) + sum(ω .* fb))
end

#--- 1f2v ---#
function flux_conserve_2d!(fw, ff, u, v, ω)
    fw[1] = sum(ω .* ff)
    fw[2] = sum(u .* ω .* ff)
    fw[3] = sum(v .* ω .* ff)
    fw[end] = 0.5 * sum((u .^ 2 .+ v .^ 2) .* ω .* ff)
end

#--- 2f2v ---#
function flux_conserve_2d!(fw, fh, fb, u, v, ω)
    fw[1] = sum(ω .* fh)
    fw[2] = sum(u .* ω .* fh)
    fw[3] = sum(v .* ω .* fh)
    fw[end] = 0.5 * (sum((u .^ 2 .+ v .^ 2) .* ω .* fh) + sum(ω .* fb))
end

#--- 1f3v ---#
function flux_conserve_3d!(fw, ff, u, v, w, ω)
    fw[1] = sum(ω .* ff)
    fw[2] = sum(u .* ω .* ff)
    fw[3] = sum(v .* ω .* ff)
    fw[4] = sum(w .* ω .* ff)
    fw[end] = 0.5 * sum((u .^ 2 .+ v .^ 2 + w .^ 2) .* ω .* ff)
end

"""
$(SIGNATURES)

Discrete moments of particle distribution
"""
discrete_moments(f, ω) = sum(@. ω * f)

"""
$(SIGNATURES)
"""
discrete_moments(f, u, ω, n) = sum(@. ω * u^n * f)

"""
$(SIGNATURES)

Calculate pressure from particle distribution function
"""
pressure(f, prim, u, ω) = sum(@. ω * (u - prim[2])^2 * f)

"""
$(SIGNATURES)
"""
pressure(h, b, prim, u, ω, K) =
    (sum(@. ω * (u - prim[2])^2 * h) + sum(@. ω * b)) / (K + 1.0)

"""
$(SIGNATURES)
"""
pressure(h, b, prim, u, v, ω, K) =
    (sum(@. ω * ((u - prim[2])^2 + (v - prim[3])^2) * h) + sum(@. ω * b)) / (K + 2.0)

"""
$(SIGNATURES)

Calculate stress tensor from particle distribution function
"""
stress(f, prim, u, ω) = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f)

"""
$(SIGNATURES)
"""
function stress(f, prim, u, v, ω)
    P = similar(prim, 2, 2)

    P[1, 1] = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f)
    P[1, 2] = sum(@. ω * (u - prim[2]) * (v - prim[3]) * f)
    P[2, 1] = P[1, 2]
    P[1, 2] = sum(@. ω * (v - prim[3]) * (v - prim[3]) * f)

    return P
end

"""
$(SIGNATURES)

Calculate heat flux from particle distribution function

Multiple dispatch doesn't consider unstructured multi-dimensional velocity space.
In that case a new method needs to be defined.

1F1V
"""
heat_flux(f, prim, u, ω) = 0.5 * sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * f) # 1F1V

"""
$(SIGNATURES)

2F1V
"""
heat_flux(
    h::X,
    b::X,
    prim::Y,
    u::Z,
    ω::Z,
) where {X<:AA{<:FN,1},Y<:AA{<:Real,1},Z<:AA{<:FN,1}} =
    0.5 * (sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * h) + sum(@. ω * (u - prim[2]) * b))

"""
$(SIGNATURES)

3F1V Rykov model
"""
function heat_flux(
    h::X,
    b::X,
    r::X,
    prim::Y,
    u::Z,
    ω::Z,
) where {X<:AA{<:FN,1},Y<:AA{<:Real,1},Z<:AA{<:FN,1}}

    q = similar(h, 2)

    q[1] =
        0.5 *
        (sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * h) + sum(@. ω * (u - prim[2]) * b))
    q[2] = 0.5 * (sum(@. ω * (u - prim[2]) * r))

    return q

end

"""
$(SIGNATURES)

1F2V
"""
function heat_flux(
    h::X,
    prim::Y,
    u::Z,
    v::Z,
    ω::Z,
) where {X<:AA{<:FN,2},Y<:AA{<:Real,1},Z<:AA{<:FN,2}}

    q = similar(h, 2)
    q[1] = 0.5 * sum(@. ω * (u - prim[2]) * ((u - prim[2])^2 + (v - prim[3])^2) * h)
    q[2] = 0.5 * sum(@. ω * (v - prim[3]) * ((u - prim[2])^2 + (v - prim[3])^2) * h)

    return q

end

"""
$(SIGNATURES)

2F2V
"""
function heat_flux(
    h::X,
    b::X,
    prim::Y,
    u::Z,
    v::Z,
    ω::Z,
) where {X<:AA{<:FN,2},Y<:AA{<:Real,1},Z<:AA{<:FN,2}}

    q = similar(h, 2)

    q[1] =
        0.5 * (
            sum(@. ω * (u - prim[2]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (u - prim[2]) * b)
        )
    q[2] =
        0.5 * (
            sum(@. ω * (v - prim[3]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (v - prim[3]) * b)
        )

    return q

end

"""
$(SIGNATURES)

3F2V Rykov model
"""
function heat_flux(
    h::X,
    b::X,
    r::X,
    prim::Y,
    u::Z,
    v::Z,
    ω::Z,
) where {X<:AA{<:FN,2},Y<:AA{<:Real,1},Z<:AA{<:FN,2}}

    q = similar(h, 4)

    q[1] =
        0.5 * (
            sum(@. ω * (u - prim[2]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (u - prim[2]) * b)
        )
    q[2] =
        0.5 * (
            sum(@. ω * (v - prim[3]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (v - prim[3]) * b)
        )
    q[3] = 0.5 * sum(@. ω * (u - prim[2]) * r)
    q[4] = 0.5 * sum(@. ω * (v - prim[3]) * r)

    return q

end

"""
$(SIGNATURES)

1F3V
"""
function heat_flux(
    f::X,
    prim::Y,
    u::Z,
    v::Z,
    w::Z,
    ω::Z,
) where {X<:AA{<:FN,3},Y<:AA{<:Real,1},Z<:AA{<:FN,3}}

    q = similar(f, 3)

    q[1] =
        0.5 * sum(
            @. ω *
               (u - prim[2]) *
               ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) *
               f
        )
    q[2] =
        0.5 * sum(
            @. ω *
               (v - prim[3]) *
               ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) *
               f
        )
    q[3] =
        0.5 * sum(
            @. ω *
               (w - prim[4]) *
               ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) *
               f
        )

    return q

end


text:
using BAT
using Test

using Distributions
using ValueShapes
using IntervalSets
using LinearAlgebra: Diagonal, ones

@testset "ahmi_integration" begin
    function test_integration(algorithm::IntegrationAlgorithm, title::String,
                              dist::Distribution; val_expected::Real=1.0,
                              val_rtol::Real=3.5, err_max::Real=0.2)
        @testset "$title" begin
            samplingalg = MCMCSampling(
                mcalg = MetropolisHastings(),
                trafo = NoDensityTransform(),
                nsteps = 2*10^5,
                burnin = MCMCMultiCycleBurnin(nsteps_per_cycle = 10^5, max_ncycles = 60)
            )
            sample = bat_sample(dist, samplingalg).result
            sample_integral = bat_integrate(sample, algorithm).result

            @test isapprox(sample_integral.val, val_expected, atol=val_rtol*sample_integral.err)
            @test sample_integral.err < err_max
        end
    end
    test_integration(AHMIntegration(), "funnel distribution", BAT.FunnelDistribution(), val_rtol = 15)
    test_integration(AHMIntegration(), "multimodal student-t distribution", BAT.MultimodalStudentT(), val_rtol = 50)
    test_integration(AHMIntegration(), "Gaussian shell", BAT.GaussianShell(), val_rtol = 15)
    test_integration(AHMIntegration(), "MvNormal", MvNormal(Diagonal(ones(5))), val_rtol = 15)

    @testset "ahmi_integration_defaults" begin
        bsample = @inferred(bat_sample(product_distribution([Normal() for i in 1:3]), IIDSampling(nsamples = 10^4))).result
        @test isapprox(bat_integrate(bsample).result.val, 1.0, rtol=15)
        eff_sample_size_dsv = @inferred(bat_eff_sample_size(bsample)).result
        eff_sample_size_aosa = @inferred(bat_eff_sample_size(bsample.v)).result
        for i in eachindex(@inferred(bat_eff_sample_size(bsample)).result)
            @test 8000 <= eff_sample_size_dsv[i] <= 10^4
            @test eff_sample_size_dsv[i] == eff_sample_size_aosa[i]
        end
    end
end

ext: jl


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

__precompile__(true)

module H3DDetectorSystems

using Sockets

using ArraysOfArrays
using BufferedStreams
using Unitful

include("polaris_data.jl")
include("datastream.jl")
include("polaris_detector.jl")

end # module


text:
"""
    dualgraph(dist::Matrix{Float64}; method::Symbol = :inverse, σ::Float64 = 1.0)

build the dual graph's weight matrix based on the given non-trivial eigenvector
metric.

# Input Arguments
- `dist::Matrix{Float64}`: eigenvector distance matrix
- `method::Symbol`: default is by taking inverse of the distance between
    eigenvectors. Ways to build the dual graph edge weights. Option: `:inverse`,
    `:gaussian`.
- `σ::Float64`: default is `1.0`. Gaussian variance parameter.

# Output Argument
- `G_star::GraphSig`: A `GraphSig` object containing the weight matrix of the
    dual graph.

"""
function dualgraph(dist::Matrix{Float64}; method::Symbol = :inverse, σ::Float64 = 1.0)
    N = Base.size(dist, 1)
    W_star = zeros(N, N)
    if method == :inverse
        for i = 1:(N - 1), j = (i + 1):N
            W_star[i, j] = 1 / dist[i, j]
        end
    elseif method == :gaussian
        for i = 1:(N - 1), j = (i + 1):N
            W_star[i, j] = exp(-dist[i, j] / σ^2)
        end
    else
        error("method must be :inverse or :gaussian.")
    end
    W_star = W_star + W_star'
    G_star = GraphSig(sparse(W_star); name = "dual graph")
    return G_star
end


text:
try
    using ArgParse
catch e
    println("Package ArgParse is required, but not installed. Install now? [y/n]")
    if readline() == "y"
        import Pkg
        Pkg.add("ArgParse")
        using ArgParse
    else
        exit()
    end
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--nprocs", "-p"
            help = "number of worker processors"
            arg_type = Int
            default = 1
        "--verbose", "-v"
            help = "print the solver outputs"
            action = :store_true
        "--gcoff", "-g"
            help = "turn of the julia garbage collector"
            action = :store_true
        "--quick", "-q"
            help = "run tests with reduced number of instances"
            action = :store_true
        "--decode", "-d"
            help = "decode the cutest instances"
            action = :store_true
        "testsets"
            help = "testsets for benchmark (separated by comma). possible values: cutest, power"
            required = true
        "solvers"
            help = "solvers to benchmark (separated by comma). possible values: current, master, ipopt"
            required = true
    end

    return parse_args(s)
end

function main()
    pargs = parse_commandline()

    CLASSES = split(pargs["testsets"],",")
    SOLVERS = split(pargs["solvers"],",")

    # sanity check
    issubset(CLASSES,["cutest","power"]) || error("argument testsets is incorrect")
    issubset(SOLVERS,["current","master","ipopt","knitro"]) || error("argument solvers is incorrect")

    PROJECT_PATH = dirname(@__FILE__)

    cp(
        joinpath(PROJECT_PATH, ".Project.toml"),
        joinpath(PROJECT_PATH, "Project.toml"),
        force=true
    )

    for class in CLASSES
        for solver in SOLVERS
            launch_script = joinpath(PROJECT_PATH, "benchmark-$class.jl")
            run(`julia --project=$PROJECT_PATH $launch_script $(pargs["nprocs"]) $solver $(pargs["verbose"]) $(pargs["quick"]) $(pargs["gcoff"]) $(pargs["decode"])`)
        end
    end

    run(Cmd([["julia","--project=$PROJECT_PATH","plot.jl"];String.(CLASSES);String.(SOLVERS)]))

end

main()


text:
# Autogenerated wrapper script for TauDEM_jll for i686-w64-mingw32-cxx03
export aread8, areadinf, d8flowdir, d8flowpathextremeup, d8hdisttostrm, dinfavalanche, dinfconclimaccum, dinfdecayaccum, dinfdistdown, dinfdistup, dinfflowdir, dinfrevaccum, dinftranslimaccum, dinfupdependence, dropanalysis, gagewatershed, gridnet, lengtharea, moveoutletstostrm, peukerdouglas, pitremove, slopearea, slopearearatio, slopeavedown, streamnet, threshold

using GDAL_jll
using MPICH_jll
using MicrosoftMPI_jll
JLLWrappers.@generate_wrapper_header("TauDEM")
JLLWrappers.@declare_executable_product(aread8)
JLLWrappers.@declare_executable_product(areadinf)
JLLWrappers.@declare_executable_product(d8flowdir)
JLLWrappers.@declare_executable_product(d8flowpathextremeup)
JLLWrappers.@declare_executable_product(d8hdisttostrm)
JLLWrappers.@declare_executable_product(dinfavalanche)
JLLWrappers.@declare_executable_product(dinfconclimaccum)
JLLWrappers.@declare_executable_product(dinfdecayaccum)
JLLWrappers.@declare_executable_product(dinfdistdown)
JLLWrappers.@declare_executable_product(dinfdistup)
JLLWrappers.@declare_executable_product(dinfflowdir)
JLLWrappers.@declare_executable_product(dinfrevaccum)
JLLWrappers.@declare_executable_product(dinftranslimaccum)
JLLWrappers.@declare_executable_product(dinfupdependence)
JLLWrappers.@declare_executable_product(dropanalysis)
JLLWrappers.@declare_executable_product(gagewatershed)
JLLWrappers.@declare_executable_product(gridnet)
JLLWrappers.@declare_executable_product(lengtharea)
JLLWrappers.@declare_executable_product(moveoutletstostrm)
JLLWrappers.@declare_executable_product(peukerdouglas)
JLLWrappers.@declare_executable_product(pitremove)
JLLWrappers.@declare_executable_product(slopearea)
JLLWrappers.@declare_executable_product(slopearearatio)
JLLWrappers.@declare_executable_product(slopeavedown)
JLLWrappers.@declare_executable_product(streamnet)
JLLWrappers.@declare_executable_product(threshold)
function __init__()
    JLLWrappers.@generate_init_header(GDAL_jll, MPICH_jll, MicrosoftMPI_jll)
    JLLWrappers.@init_executable_product(
        aread8,
        "bin\\aread8.exe",
    )

    JLLWrappers.@init_executable_product(
        areadinf,
        "bin\\areadinf.exe",
    )

    JLLWrappers.@init_executable_product(
        d8flowdir,
        "bin\\d8flowdir.exe",
    )

    JLLWrappers.@init_executable_product(
        d8flowpathextremeup,
        "bin\\d8flowpathextremeup.exe",
    )

    JLLWrappers.@init_executable_product(
        d8hdisttostrm,
        "bin\\d8hdisttostrm.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfavalanche,
        "bin\\dinfavalanche.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfconclimaccum,
        "bin\\dinfconclimaccum.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfdecayaccum,
        "bin\\dinfdecayaccum.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfdistdown,
        "bin\\dinfdistdown.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfdistup,
        "bin\\dinfdistup.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfflowdir,
        "bin\\dinfflowdir.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfrevaccum,
        "bin\\dinfrevaccum.exe",
    )

    JLLWrappers.@init_executable_product(
        dinftranslimaccum,
        "bin\\dinftranslimaccum.exe",
    )

    JLLWrappers.@init_executable_product(
        dinfupdependence,
        "bin\\dinfupdependence.exe",
    )

    JLLWrappers.@init_executable_product(
        dropanalysis,
        "bin\\dropanalysis.exe",
    )

    JLLWrappers.@init_executable_product(
        gagewatershed,
        "bin\\gagewatershed.exe",
    )

    JLLWrappers.@init_executable_product(
        gridnet,
        "bin\\gridnet.exe",
    )

    JLLWrappers.@init_executable_product(
        lengtharea,
        "bin\\lengtharea.exe",
    )

    JLLWrappers.@init_executable_product(
        moveoutletstostrm,
        "bin\\moveoutletstostrm.exe",
    )

    JLLWrappers.@init_executable_product(
        peukerdouglas,
        "bin\\peukerdouglas.exe",
    )

    JLLWrappers.@init_executable_product(
        pitremove,
        "bin\\pitremove.exe",
    )

    JLLWrappers.@init_executable_product(
        slopearea,
        "bin\\slopearea.exe",
    )

    JLLWrappers.@init_executable_product(
        slopearearatio,
        "bin\\slopearearatio.exe",
    )

    JLLWrappers.@init_executable_product(
        slopeavedown,
        "bin\\slopeavedown.exe",
    )

    JLLWrappers.@init_executable_product(
        streamnet,
        "bin\\streamnet.exe",
    )

    JLLWrappers.@init_executable_product(
        threshold,
        "bin\\threshold.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module api
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("api_pb.jl")
end


text:
weights = SimpleWeightedGraphs.weights

"""
Find the shortest path between two nodes of an AbstractGraph by using the
\'a star\' method. Returns a vector of AbstractEdge. If no path is found,
returns an empty vector.
"""
function shortest_path(graph::AbstractGraph, src::Int64, dst::Int64)
    path = a_star(SimpleDiGraph(weights(graph)), src, dst, weights(graph))
    return path
end
# Change SimpleDiGraph to SimpleWeightedDiGraph

"""
Finds the total weight/length of a path for a given AbstractGraph.
Not to be confused with get_pathcv(), which finds the vector of costs for a
path in a QNetwork.

!!! warning
    For whatever reason, SimpleWeightedDiGraph indexes weights with [dest, src].
    Keep in mind this path_length will likely not work if you pass in anything
    but a SimpleWeightedDiGraph.
"""
function path_length(graph::AbstractGraph, path::Vector{Tuple{Int64, Int64}})::Float64
    dist = 0.0
    for edge in path
        # SimpleWeightedDiGraph indexes by [dst, src]. Be careful here!
        dist += graph.weights[edge[2], edge[1]]
    end
    return dist
end

function path_length(graph::AbstractGraph, path::Vector{LightGraphs.SimpleGraphs.SimpleEdge{Int64}})::Float64
    dist = 0.0
    for edge in path
        # SimpleWeightedDiGraph indexes by [dst, src]. Be careful here!
        dist += graph.weights[edge.dst, edge.src]
    end
    return dist
end

"""
Removes the shortest path between two nodes of an abstract graph and returns the
path along with its length. If no path exists, returns nothing for both.
"""
function remove_shortest_path!(graph::AbstractGraph, src::Int64, dst::Int64)

    # Get shortest path. If no path exists, return nothing for path and length
    path = shortest_path(graph, src, dst)
    if length(path) == 0
        return nothing
    end

    path_len = path_length(graph, path)

    # Remove path
    for edge in path
        status = hard_rem_edge!(graph, edge.src, edge.dst)
        if status == false
            error("Failed to remove edge in path")
        end
    end

    return path, path_len
end

"""
Removes the shortest path between two nodes of a QNetwork *with respect to a
given cost* and returns the path with its corresponding cost vector.
"""
function remove_shortest_path!(network::QNetwork, cost_id::String, src::Int64, dst::Int64)
    @assert cost_id in keys(zero_costvector()) "Invalid cost"

    # Find shortest path in terms of the given cost
    g = network.graph[cost_id]
    path = shortest_path(g, src, dst)
    # If no path exists, return nothing for path and path_cv
    if length(path) == 0
        return nothing, nothing
    end

    # Find the costs associated with the path
    path_cv = get_pathcv(network, path)

    # Remove the path from the network graphs
    for cost_id in keys(zero_costvector())
        graph = network.graph[cost_id]
        for edge in path
            status = hard_rem_edge!(graph, edge.src, edge.dst)
            if status == false
                error("Failed to remove edge in path")
            end
        end
    end
    return path, path_cv
end

"""
Remove the shortest_path of a TemporalGraph, making sure not to remove the edge
between the asynchronous node and its temporal counterpart. When specifying
src and dst, specify them as you would if the graph were non-temporal. E.G. if
you have a 1x2 network with 2 timesteps and you wanted to find the best path
between the two nodes, you would specify src:1, dst:2
"""
function remove_shortest_path!(tempnet::QuNet.TemporalGraph, cost_id::String,
    src::Int64, dst::Int64)

    @assert cost_id in keys(zero_costvector()) "Invalid cost"

    # Check if src or dst nodes are asynchronus
    async_src = false
    async_dst = false

    if src > tempnet.nv * tempnet.steps
        async_src = true
    end
    if dst > tempnet.nv * tempnet.steps
        async_dst = true
    end

    # TODO: Remove this
    #     # Check that src and dst are nodes within range(1, size-of-static-graph)
    #     if (src > tempnet.nv && dst > tempnet.nv)
    #     @assert 1==0 "(src::$src) and (dst::$dst) async_pair = true: nodes must be asynchronous"
    #     # Reindex src and dst so that they point to their asynchronus counterparts
    #     src += tempnet.nv * tempnet.steps
    #     dst += tempnet.nv * tempnet.steps
    # end

    # Find shortest path in terms of the given cost
    t = tempnet.graph[cost_id]
    path = shortest_path(t, src, dst)

    # If length(path) == 0, no path exists. Return nothing
    if length(path) == 0
        return nothing, nothing
    end

    # If either src or dst was asynchronus:
    # Pop the asynchronus edge so it doesn't get removed.
    if async_src == true
        popfirst!(path)
    end
    if async_dst == true
        pop!(path)
    end

    # Find the costs associated with the path:
    path_cv = get_pathcv(tempnet, path)

    # Remove the shortest path
    for cost_id in keys(zero_costvector())
        graph = tempnet.graph[cost_id]
        for edge in path
            status = hard_rem_edge!(graph, edge.src, edge.dst)
            if status == false
                error("Failed to remove edge in path")
            end
        end
    end
    return path, path_cv
end

"""
greedy_multi_path! is an entanglement routing strategy for a quantum network
with n end user pairs.

1. Each user pair has an associated "bundle" of paths that starts out empty.
For each user pair, find the shortest path between them, and add its cost to the
bundle. If no path exists, add "nothing" to the bundle.

2. When all the bundles are filled up to the maximum number "maxpaths", use an
entanglement purification method between the path costs. Add the resulting cost
vector to the array pur_paths. If no paths exist for a given bundle, increment
the collision count by one, and add "nothing" to pur_paths

3. Return the array of purified cost vectors and the collision count.
"""
function greedy_multi_path!(network::QNetwork, purification_method,
    users, maxpaths::Int64=3)

    # List of paths for each userpair
    pathset = [Vector() for i in 1:length(users)]
    # List of path costs for each userpair
    path_costs = [Vector{Dict{Any,Any}}() for i in 1:length(users)]

    for i in 1:maxpaths
        for (userid, user) in enumerate(users)
            src = user[1]; dst = user[2]
            # Remove the shortest path in terms of Z-dephasing and get its cost vector.
            path, path_cv = remove_shortest_path!(network, "Z", src, dst)
            # If pathcv is nothing, no path was found.
            if path_cv == nothing
                break
            else
                push!(pathset[userid], path)
                push!(path_costs[userid], path_cv)
            end
        end
    end

    # A tally of the number of paths each end-user purifies together.
    pathuse_count = [0 for i in 0:maxpaths]
    # An array of purified cost vectors for each end-user.
    pur_paths = []

    # Purify end-user paths
    for userpaths in path_costs
        # Increment the tally for the number of paths beting purified
        len = length(userpaths)
        pathuse_count[len + 1] += 1
        # If len == 0, no paths were found between the end-user.
        if len == 0
            push!(pur_paths, nothing)
        # Otherwise, purify the paths
        else
            purcost::Dict{Any, Any} = purification_method(userpaths)
            # Convert purcost from decibels to metric form
            purcost = convert_costs(purcost)
            push!(pur_paths, purcost)
        end
    end
    return pathset, pur_paths, pathuse_count
end

function greedy_multi_path!(network::QuNet.TemporalGraph, purification_method,
    users, maxpaths::Int64=3)

    # For TemporalGraph, paths must arrive at identical time_depth
    # route_layer_known is false until first path is routed.
    route_layer_known = false

    # List of paths for each userpair
    pathset = [Vector() for i in 1:length(users)]
    # List of path costs for each userpair
    path_costs = [Vector{Dict{Any,Any}}() for i in 1:length(users)]

    for i in 1:maxpaths
        for (userid, user) in enumerate(users)
            src = user[1]; dst = user[2]
            # Remove the shortest path in terms of Z-dephasing and get its cost vector.
            path, path_cv = remove_shortest_path!(network, "Z", src, dst)
            # If pathcv is nothing, no path was found.
            if path_cv == nothing
                break
            else
                push!(pathset[userid], path)
                push!(path_costs[userid], path_cv)
            end

            # If we found a path, and we haven't fixed a routing time:
            if path_cv != nothing && route_layer_known == false
                route_layer_known = true
                last_edge = last(path)
                # If last node of the path is asynchronus:
                    # Remove all async edges except for the one at T = depth
                    # This means all future paths will have to route to the same time
                if last_edge.dst > network.nv * network.steps
                    node = last_edge.src
                    depth = QuNet.node_timedepth(T.nv, T.steps, node)
                    # Remove all asynchronus edges except for that depth
                    QuNet.fix_async_nodes_in_time(T, [node])
                end
            end

        end
    end

    # A tally of the number of paths each end-user purifies together.
    pathuse_count = [0 for i in 0:maxpaths]
    # An array of purified cost vectors for each end-user.
    pur_paths = []

    # Purify end-user paths
    for userpaths in path_costs
        # Increment the tally for the number of paths beting purified
        len = length(userpaths)
        pathuse_count[len + 1] += 1
        # If len == 0, no paths were found between the end-user.
        if len == 0
            push!(pur_paths, nothing)
        # Otherwise, purify the paths
        else
            purcost::Dict{Any, Any} = purification_method(userpaths)
            # Convert purcost from decibels to metric form
            purcost = convert_costs(purcost)
            push!(pur_paths, purcost)
        end
    end
    return pathset, pur_paths, pathuse_count
end

"""
Find the maximum timedepth reached by a given pathset
"""
function max_timedepth(pathset, T)
    max_depth = 1
    for bundle in pathset
        for path in bundle
            edge = last(path)
            node = edge.dst
            # Check if node is temporal. If it is, use 2nd last node in path instead
            if node > T.nv * T.steps
                node = edge.src
            end
            # use node - 1 here because if node % T.nv == 0, depth is off by one
            depth = QuNet.node_timedepth(T.nv, T.steps, node)
            if depth > max_depth
                max_depth = depth
            end
        end
    end
    return max_depth
end

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
### A Pluto.jl notebook ###
# v0.17.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 43e5cbb6-5557-492e-9aeb-f7cbde220e15
begin
    import Pkg
    Pkg.activate() # disable Pluto's package management
end

# ╔═╡ e2d1cb38-2698-11eb-2e59-632faa201d6f
begin
    using DelimitedFiles
    using PlutoUI

    ENV["GRDISPLAY"] = "pluto"
    using GR
end

# ╔═╡ 49d120ea-2699-11eb-203d-fddbb33bb08b
begin
    url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
    download(url, "covid19.csv")
end;

# ╔═╡ 4b5f800a-2699-11eb-2528-b5917c0461cf
begin
    data = readdlm("covid19.csv", ',')
    ncountries, ncols = size(data)
    ndays = ncols - 4
end;

# ╔═╡ 1c4aec86-269a-11eb-2204-098352e8267c
@bind countries MultiSelect(["Germany", "Austria", "Belgium", "Netherlands", "France", "Italy", "Spain", "US"], default=["Germany"])

# ╔═╡ 75ad9770-2699-11eb-0db7-5948da39d6e1
begin
    cummulated = Dict()
    for i in 1:ncountries
        country = data[i,2]
        if country in countries
            if !haskey(cummulated, country) cummulated[country] = zeros(ndays) end
            cummulated[country] .+= collect(data[i,5:end])
        end
    end
end

# ╔═╡ 86c70334-2699-11eb-1b52-cfa50b314e2b
begin
    day = collect(Float64, 1:ndays);
    confirmed = hcat([cummulated[country] for country in countries]...)
end;

# ╔═╡ 918f640a-2699-11eb-02d8-998ed76c614a
plot(day, confirmed, xlim=(0, ndays+1), ylim=(10, 100_000_000), ylog=true,
     title="Confirmed SARS–CoV–2 infections", xlabel="Day", ylabel="Confirmed",
     labels=countries, location=4)

# ╔═╡ Cell order:
# ╠═43e5cbb6-5557-492e-9aeb-f7cbde220e15
# ╠═e2d1cb38-2698-11eb-2e59-632faa201d6f
# ╠═49d120ea-2699-11eb-203d-fddbb33bb08b
# ╠═4b5f800a-2699-11eb-2528-b5917c0461cf
# ╠═75ad9770-2699-11eb-0db7-5948da39d6e1
# ╠═86c70334-2699-11eb-1b52-cfa50b314e2b
# ╠═1c4aec86-269a-11eb-2204-098352e8267c
# ╠═918f640a-2699-11eb-02d8-998ed76c614a


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    BSplineCurve{P,S,N,M} <: Spline{P,S,N,M}

`N` is the spatial dimension of the curve.
`M` is the curve order, i.e., the highest power of the parameterizing variable, `u`.
All curve segments are assumed to be of the same order.

```julia
BSplineCurve{P,S,N,M}(knots::KnotVector{S}, controlpoints::AbstractArray{MVector{N,S},1})
```
"""
struct BSplineCurve{P, S, N, M} <: Spline{P, S, N, M}
    knotvector::KnotVector{S}
    controlpolygon::Array{MVector{N, S}, 1} # may have to allow different orders on different curve segments

    function BSplineCurve{P, S, N, M}(
        knots::KnotVector{S},
        controlpoints::AbstractArray{MVector{N, S}, 1},
    ) where {P, S, N, M}
        bsplineinvariant(numknots(knots), length(controlpoints), M)
        newpoints = [MVector{N, S}(point) for point in controlpoints]
        return new{P, S, N, M}(knots, newpoints)
    end
end

numknots(curve::BSplineCurve) = numknots(curve.knotvector)

function findspan(curve::BSplineCurve{T, S, N, M}, u) where {T, S, N, M}
    return findspan(curve.knotvector, M, u)
end

function numspans(curve::BSplineCurve{T, S, N, M}) where {T, S, N, M}
    return numknots(curve) - M - 1 # number of spans = numknots - curveorder - 1
end

spatialdimension(::Spline{T, S, N, M}) where {T, S, N, M} = N
curveorder(::Spline{T, S, N, M}) where {T, S, N, M} = M

function point(
    curve::BSplineCurve{P, S, N, M},
    u::T,
)::SVector{N, S} where {T <: Real, P, S, N, M}
    # returns the raw point type of the curve - for Homogeneous curve types this will be a homogenous point
    span = findspan(curve, u)
    bases = basisfunctions(curve.knotvector, u, M)
    c = zeros(MVector{N, S})
    for i in 1:(M + 1)
        pt = curve.controlpolygon[span - (M + 1) + i]
        c .= c .+ bases[i] .* pt
    end
    return SVector{N, S}(c)
end

function euclideanpoint(
    curve::BSplineCurve{Euclidean, S, N, M},
    ::T,
)::SVector where {T <: Real, S, N, M}
    # might seem weird to have this function for Euclidean curves. But in user code can use euclidean point without worrying about whether curve is homogeneous and still get correct results.
    return point(curve)
end

function euclideanpoint(
    curve::BSplineCurve{Rational, S, N, M},
    u::T,
)::SVector where {T <: Real, S, N, M}
    return toeuclidean(point(curve, u))
end

"""
    BSplineSurface{P,S,N,M} <: SplineSurface{P,S,N,M}

Curve order is the same in the u and v direction and fixed over all spans.
u and v knot vectors are allowed to be different - _may change this to make them both the same_.

Control points in the u direction correspond to columns, with the lowest value of u corresponding to row 1.
Control points in the v direction correspond to rows, with the lowest value of v corresponding to col 1.

!!! danger
    This surface does not create a valid half-space, requires updates to function correctly.

```julia
BSplineSurface{P,S,N,M}(knots::KnotVector{S}, controlpoints::AbstractArray{<:AbstractArray{S,1},2})
BSplineSurface{P,S,N,M}(uknots::KnotVector{S}, vknots::KnotVector{S}, controlpoints::AbstractArray{<:AbstractArray{S,1},2})
```
"""
struct BSplineSurface{P, S, N, M} <: SplineSurface{P, S, N, M}
    uknotvector::KnotVector{S}
    vknotvector::KnotVector{S}
    controlpolygon::Array{MVector{N, S}, 2} # may have to allow different orders on different curve segments

    function BSplineSurface{P, S, N, M}(
        knots::KnotVector{S},
        controlpoints::AbstractArray{<:AbstractArray{S, 1}, 2},
    ) where {P, S, N, M}
        return BSplineSurface(knots, knots, controlpoints)
    end

    function BSplineSurface{P, S, N, M}(
        uknots::KnotVector{S},
        vknots::KnotVector{S},
        controlpoints::AbstractArray{<:AbstractArray{S, 1}, 2},
    ) where {P, S, N, M}
        # m == n + M + 1
        # relation between number of knots = m + 1, number of control points = n + 1, and curve order M
        n_u, n_v = size(controlpoints)
        bsplineinvariant(numknots(uknots), n_u, M)
        bsplineinvariant(numknots(vknots), n_v, M)
        newpoints = [MVector{N, S}(point) for point in controlpoints]
        return new{P, S, N, M}(uknots, vknots, newpoints)
    end
end

function numknots(surf::BSplineSurface)
    return (numknots(surf.uknotvector), numknots(surf.vknotvector))
end

function uvrange(surface::BSplineSurface)
    return (urange(surface.uknotvector), urange(surface.vknotvector))
end

function findspan(surface::BSplineSurface{T, S, N, M}, u, v) where {T, S, N, M}
    return (
        findspan(surface.uknotvector, M, u),
        findspan(surface.vknotvector, M, v),
    )
end

function insertknot(
    knots::AbstractArray{S, 1},
    knotindex::Int,
    controlpoints::Array{MVector{N, S}, 1},
    curveorder::Int,
) where {N, S}
    knotvalue = knots[knotindex]
    M = curveorder
    K = [zeros(MVector{N, S}) for _ in 1:(M + 1)]
    newknots = vcat(knots[1:knotindex], [knotvalue], knots[(knotindex + 1):end])
    newcontrolpoints = vcat(
        controlpoints[1:(knotindex - M)],
        K,
        controlpoints[(knotindex + 1):end],
    )

    for i in (knotindex - M + 1):knotindex
        αᵢ = (knotvalue - knots[i]) / (knots[i + M] - knots[i])
        @. newcontrolpoints[i] =
            αᵢ * controlpoints[i] + (1 - αᵢ) * controlpoints[i - 1]
    end
    for i in (knotindex + 1):length(newcontrolpoints)
        newcontrolpoints[i] .= controlpoints[i - 1]
    end

    return (newknots, newcontrolpoints)
end

function insertknot(
    curve::BSplineCurve{T, S, N, M},
    knotindex::Int,
) where {T, S, N, M}
    newknots, newcontrolpoints =
        insertknot(curve.knotvector.knots, knotindex, curve.controlpolygon, M)
    return BSplineCurve{T, S, N, M}(KnotVector{S}(newknots), newcontrolpoints)
end

function insertknots(curve::BSplineCurve{T, S, N, M}) where {T, S, N, M}
    newknots, newcontrolpoints =
        insertknots(curve.knotvector.knots, curve.controlpolygon, M)
    return BSplineCurve{T, S, N, M}(newknots, newcontrolpoints)
end

function insertknots(surf::BSplineSurface{T, S, N, M}) where {T, S, N, M}
    # TODO avoid using '
    newvknots, temppoints =
        expandcontrolpoints(surf.vknotvector.knots, surf.controlpolygon, M)
    newuknots, temppoints =
        expandcontrolpoints(surf.uknotvector.knots, copy(temppoints'), M)
    # TODO this adjoint is a pain for messing up all the types, ideally we wouldn't need to deepcopy
    # should have an arg for expandcontrolpoints that controls direction or something
    return BSplineSurface{T, S, N, M}(
        KnotVector{S}(newuknots),
        KnotVector{S}(newvknots),
        copy(temppoints'),
    )
end

function expandcontrolpoints(
    vknots::AbstractArray{S, 1},
    controlpolygon::Array{MVector{N, S}, 2},
    curveorder::Int,
) where {N, S}
    upts, vpts = size(controlpolygon)
    temp = Array{Array{MVector{N, S}, 1}, 1}(undef, 0) # TODO static array
    let newknots = vknots
        for i in 1:upts
            newknots, newcontrolpoints =
                insertknots(vknots, view(controlpolygon, i, :), curveorder)
            # not the most efficient since new knot arrays are being computed for reach row when all rows of knots are the same. Optimize later if necessary.
            # println("temp $temp newpts $newcontrolpoints")
            push!(temp, newcontrolpoints)
        end
        newcontrolpolygon = hcat(temp[:]...) # TODO more efficient way
        # make new 2D control point array
        return newknots, newcontrolpolygon
    end
end

function insertknots(
    knots::AbstractArray{S, 1},
    controlpoints::AbstractArray{MVector{N, S}, 1},
    curveorder::Int,
) where {N, S}
    insertionknots = knotstoinsert(knots, curveorder)
    newknots = copy(knots)
    newcontrolpoints = Array{MVector{N, S}, 1}(controlpoints)
    offset = 0
    for insertion in insertionknots
        for i in 1:insertion[2]
            newknots, newcontrolpoints = insertknot(
                newknots,
                insertion[1] + offset,
                newcontrolpoints,
                curveorder,
            )
        end
        offset += insertion[2]  # inserting knots changes the index where the next set of knots needs to be inserted.
    end
    return (newknots, newcontrolpoints)
end

function knotstoinsert(curve::Spline{T, S, N, M}) where {T, S, N, M}
    return knotstoinsert(curve.knotvector)
end

function knotstoinsert(knots::AbstractArray{S, 1}, curveorder::Int) where {S}
    numknots = length(knots)
    index = curveorder + 2
    knotcounts = Array{Tuple{Int64, Int64}, 1}(undef, 0)

    while index < numknots - curveorder
        stop = index + 1
        while knots[index] == knots[stop]
            stop += 1       # stop will have the index after the last repeated knot in this sequence
        end
        numtoinsert = (curveorder) - (stop - index)
        if numtoinsert != 0
            push!(knotcounts, (index, numtoinsert))
        end
        # println("index $index stop $stop")
        index = stop
    end
    return knotcounts
end

function tobeziersegments(curve::BSplineCurve{T, S, N, M}) where {T, S, N, M}
    # returns an array of arrays of groups of M+1 Bezier control points, where M is the curve order.
    return tobeziersegments(curve.knotvector.knots, curve.controlpolygon, M)
end

function tobeziersegments(
    knots::AbstractArray{S, 1},
    controlpoints::Array{MVector{N, S}, 1},
    curveorder::Int,
) where {N, S}
    # returns an array of arrays of groups of M+1 Bezier control points, where M is the curve order. Inefficient because it creates a new curve for every knot insertion. Optimize when necessary
    M = curveorder
    _, newcontrolpoints = insertknots(knots, controlpoints, curveorder)
    # each knot is now inserted M times, except for the first and last, which are repeated M+1 times.
    numsegments = length(knots) - 2 * (M + 1) + 1
    # println(newcurve.controlpolygon)
    # println("numsegments $numsegments")
    segments = Array{Array{MVector{N, S}, 1}}(undef, numsegments)

    for i in 1:numsegments
        start = (i - 1) * M + 1
        segments[i] = newcontrolpoints[start:(start + M)]
    end

    return segments
end

function tobeziersegments(surf::BSplineSurface{P, S, N, M}) where {P, S, N, M}
    newsurf = insertknots(surf)
    newcontrolpoints = newsurf.controlpolygon
    # each knot is now inserted M times, except for the first and last, which are repeated M+1 times.
    uknots, vknots = numknots(newsurf)
    numusegments = (uknots - 2 * (M + 1)) ÷ M + 1
    numvsegments = (vknots - 2 * (M + 1)) ÷ M + 1
    # println(newcurve.controlpolygon)
    # println("numsegments $numsegments")
    segments =
        Array{BezierSurface{P, S, N, M}}(undef, numusegments, numvsegments)
    # println("usges $numusegments vseg $numvsegments")
    for i in 1:numusegments
        for j in 1:numvsegments
            startu = (i - 1) * M + 1
            startv = (j - 1) * M + 1
            segments[i, j] = BezierSurface{P, S, N, M}(
                newcontrolpoints[startu:(startu + M), startv:(startv + M)],
            )
        end
    end

    return segments
end

function point(
    surface::BSplineSurface{P, S, N, M},
    u::T,
    v::T,
)::SVector{N, S} where {T <: Real, P, S, N, M}
    spanu, spanv = findspan(surface, u, v)

    ubases = basisfunctions(surface.uknotvector, u, M)
    vbases = basisfunctions(surface.vknotvector, v, M)
    sum = zeros(MVector{N, S})
    c = zeros(MVector{N, S})

    for j in 1:(M + 1)
        for i in 1:(M + 1)
            pt =
                surface.controlpolygon[spanu - (M + 1) + i, spanv - (M + 1) + j]
            c .= c .+ ubases[i] .* pt
        end
        sum .= sum .+ vbases[j] .* c
        c .= 0
    end

    return SVector{N, S}(sum)
end

function bsplineinvariant(numknots::Int, numcontrolpoints::Int, curveorder::Int)
    m = numknots - 1
    n = numcontrolpoints - 1
    @assert m == n + curveorder + 1 "This invariant should hold: m = n + curveorder + 1. Actual values: m=$m n=$n curveorder=$curveorder"
end

# Note TO SELF: need code to extract curves in the u and v directions from BSpline surfaces to work with the conversion to Bezier form.
# "converts a NURBS curve into Bezier curve segments, represented as Bezier control points. For an m degree NURBS with n knots this will return n - (m+1) + 1 Bezier segments, each with m+1 control points."
# function toBezier(curve::BSplineCurve{T,S,N,M}) where {T,S,N,M}
#     #this code follows that in the NURBS book closely. Use OffsetArrays to index by zero to avoid the nightmare of converting to indexing by 1.
#     U = curve.knots
#     Pw = curve.controlpolygon

#     Qw = OffsetArray{Array{OffestArray{Array{SVector{N,S},1}},1}}(undef,0)
#     temppoints = Array{SVector{N,S},1}(undef,0)

#     m = M + numknots(curve) + 1
#     a = M
#     b = M+1
#     nb = 0

#     for i in 0:M
#         push!(temppoints,Pw[i])
#     end
#     Qw[nb] = temppoints

#     while b < m
#         i = b
#         while b < m && U[b+1] == U[b]
#             b += 1;
#         end
#         mult = b-i+1
#         if mult < M
#             numerator = U[b] - U[a] #numerator of alpha
#             #compute and store alphas

#             for j in p:-1:mult+1
#                 alphas[j-mult-1] = numerator/(U[a+j] - U[a])
#             end

#             r = p - mult #insert knot r times
#             for j in 1:r
#                 save = r-j
#                 s = mult+j # This many new points
#                 for k in p:-1:s
#                     alpha = alphas[k-s]
#                     Qw[nb][k] = alpha*Qw[nb][k] + (1-alpha)*Qw[nb][k-1]
#                 end

#                 if b < m #control point of next segment
#                     Qw[nb+1][save] = Qw[nb][p]
#                 end
#             end
#         end

#         nb += 1 #Bezier segment completed
#         if b < m  #initialize for next segment
#             for i in p-mult:pairs
#                 Qw[nb][i] = Pw[b-p+1]
#             end
#             a = b
#             b = b+1
#         end
#     end
# end


text:
s = open("input.txt") do file
    readlines(file);
end

function testInput(inp, dims)
    return reshape(replace(inp, "\n"=>"")|>collect, dims)
end

function GetOccupiedNeighbors(i, j, board)
    count = 0;
    directions = [
        [0, 1], [1, 1], [1, 0], [1, -1], [0, -1], [-1, -1], [-1, 0], [-1, 1]
    ];

    for dir in directions
        # println(dir);
        loc = [i, j] .+ dir;
        while 0 < loc[1] <= size(board, 1) && 0 < loc[2] <= size(board, 2)
            # println(board[loc[1], loc[2]]);
            if board[loc[1], loc[2]] == '#'
                count += 1;
                break;
            elseif board[loc[1], loc[2]] == 'L'
                break;
            end

            loc .+= dir;
        end
    end

    return count;
end

function SimulateSeating(startState)
    simulating = true;
    curState = startState;
    nextState = Array{Char}(undef, 96, 91);
    nextState .= curState;
    lastChanges = -1;
    while simulating
        changes = 0;
        for i = 1:size(startState,1)
            for j = 1:size(startState,2)
                if curState[i, j] == '.'
                    continue
                end

                nCount = GetOccupiedNeighbors(i, j, curState);

                if curState[i, j] == 'L' && nCount == 0
                    nextState[i, j] = '#';
                    changes += 1;
                elseif curState[i, j] == '#' && nCount >= 5
                    nextState[i, j] = 'L';
                    changes += 1;
                else
                    nextState[i,j] = curState[i, j];
                end
            end
        end
        println(changes);

        if all(nextState .== curState)
            simulating = false;
            break;
        end

        curState .= nextState;
    end

    println("=================");
    println(sum(map(x->x=='#', curState)));
end

rowLen = length(s[1]);
seating = vcat(map(collect, s)...);
seatingV = reshape(seating, (rowLen,91));
SimulateSeating(seatingV)


text:
if lowercase(get(ENV, "CI", "false")) == "true"
    include("install_dependencies.jl")
end

using Simplices
using Simplices: SimplexSplitting

using Test
using LinearAlgebra
using Printf

println("Running tests for Simplices.jl ....")

include("even_sampling.jl")
include("test_delaunay.jl")
include("test_intersectionvol.jl")
include("tensordecomposition.jl")

println("Tests finished ....")


text:
module ColorMath

import Base: +, -, *, /, ==, hex

export RGB, HSB, CMYK, rgb255, red, green, blue, clear, alpha
export white, black, gray
export cyan, magenta, yellow
#export hue, saturation, bright, cyan, magenta, yellow, black
#export gray, purple, brown

# package code goes here
include("types.jl")
include("colors.jl")

function +(lhs::RGB, rhs::RGB)
    RGB(lhs.r + rhs.r, lhs.g + rhs.g, lhs.b + rhs.b, lhs.a)
end

function -(lhs::RGB, rhs::RGB)
    RGB(lhs.r - rhs.r, lhs.g - rhs.g, lhs.b - rhs.b, lhs.a)
end

*(lhs::RGB, m::Int) = RGB(m * lhs.r, m * lhs.g, m * lhs.b, lhs.a)
*(m::Int, lhs::RGB) = *(lhs, m)

/(lhs::RGB, m::Int) = RGB(lhs.r/m, lhs.g/m, lhs.b/m, lhs.a)

function ==(lhs::RGB, rhs::RGB)
    isapprox(lhs.r, rhs.r) && isapprox(lhs.g, rhs.g) && isapprox(lhs.b, rhs.b) && isapprox(lhs.a, rhs.a)
end

function ==(lhs::CMYK, rhs::CMYK)
    isapprox(lhs.c, rhs.c) && isapprox(lhs.m, rhs.m) && isapprox(lhs.y, rhs.y) && isapprox(lhs.k, rhs.k)
end

function hex(lhs::RGB)
    function rounddown(n::Float16)
        n = n<0 ? 0.0 : n
        uppercase(hex(round(Int, 0xffn, RoundDown), 2))
    end
    string("#", rounddown(lhs.r), rounddown(lhs.g), rounddown(lhs.b))
end

function rgb255(r::Int, g::Int, b::Int, a::Float64 = 1.0)
    RGB(r/0xff, g/0xff, b/0xff, a)
end

alpha(lhs::RGB, a::FloatInt) = RGB(lhs.r, lhs.g, lhs.b, a)

end # module

max_stars_repo_name: wookay/ColorMath

max_issues_repo_name: wookay/ColorMath

max_forks_repo_name: wookay/ColorMath


text:
function _pwf2pm_corrections_PV!(pm_data::Dict, pwf_data::Dict, software::Organon)
    for (i, bus) in pm_data["bus"]
        if bus_type_num_to_str[bus["bus_type"]] == "PV"
            filters = [
                gen -> element_status[gen["gen_status"]] == "ON",
                gen -> gen["qmin"] == gen["qmax"]
            ]
            if !isempty(generators_from_bus(pm_data, parse(Int, i); filters = filters))
                bus["bus_type"] = bus_type_str_to_num["PQ"]
                # if there is no load in this PV bus, create one to
                # later allocate this generation
                if isempty(load_from_bus(pm_data, parse(Int, i)))
                    _pwf2pm_load!(pm_data, pwf_data, parse(Int,i))
                end
                @warn "Active generator with QMIN = QMAX found in a PV bus number $i. Changing bus type from PV to PQ."
            end
        end
    end
end

function _pwf2pm_corrections_PQ!(pm_data::Dict, software::Organon)
    for (i, bus) in pm_data["bus"]
        if bus_type_num_to_str[bus["bus_type"]] == "PQ"
            filters = [
                gen -> element_status[gen["gen_status"]] == "ON", #1
                gen -> gen["qmin"] < gen["qmax"],                 #2
                gen -> gen["qmin"] == gen["qmax"]                 #3
            ]
            gen_keys_case1 = generators_from_bus(pm_data, parse(Int, i); filters = filters[[1,2]])
            gen_keys_case2 = generators_from_bus(pm_data, parse(Int, i); filters = filters[[1,3]])

            if !isempty(gen_keys_case1)
                # change bus type to PV
                bus["bus_type"] = bus_type_str_to_num["PV"]
                @warn "Active generator with QMIN < QMAX found in a PQ bus. Changing bus $i type to PV."
            elseif !isempty(gen_keys_case2)
                # change generator status to off and sum load power with gen power
                Pg, Qg = sum_generators_power_and_turn_off(pm_data, gen_keys_case2)
                load_key = load_from_bus(pm_data, parse(Int, i))
                @assert length(load_key) == 1
                # sum load power with the negative of generator power
                pm_data["load"][load_key[1]]["pd"] += - Pg
                pm_data["load"][load_key[1]]["qd"] += - Qg
                @warn "Active generator with QMIN = QMAX found in PQ bus $i. Adding generator power " *
                    "to load power and changing generator status to off."
            end
        end
    end
    return
end

function _pwf2pm_corrections_shunt!(pm_data::Dict, software::Organon; add_control_data::Bool=false)
    for (s, shunt) in pm_data["shunt"]
        bus = pm_data["bus"]["$(shunt["shunt_bus"])"]
        bus_type = bus["bus_type"]

        if add_control_data
            if !(bus_type_num_to_str[bus_type] == "PQ" && shunt["control_data"]["shunt_control_type"] == 2) # Discrete
            _fix_shunt_voltage_bounds(shunt, pm_data)
            end
        end
    end
end

function handle_min_max_value(element::Dict, spec::String, min::String, max::String)
    spec_value = element[spec]
    min_value  = element[min]
    max_value  = element[max]

    if (min_value == max_value) && min_value == 0
        return spec_value, spec_value
    end

    return min_value, max_value
end

function _pwf2pm_corrections_gen!(pm_data::Dict, pwf_data::Dict, software::Organon)
    for (i, gen) in pm_data["gen"]
        gen["qmin"], gen["qmax"] = handle_min_max_value(gen, "qg", "qmin", "qmax")
    end
end


text:
# Defines BLIS Objects and related methods.
#

export BliObj, BliObjBase

" BliObjBase resembles obj_t of BLIS' object-based interface. "
mutable struct BliObjBase
    # Basic fields.
    root::Ptr{BliObjBase}
    off::NTuple{2, BliDim}
    dim::NTuple{2, BliDim}
    diag_off::BliDoff
    info1::BliObjBits
    info2::BliObjBits
    enum_size::BliSiz
    buffer::Ptr{Nothing}
    rs::BliInc
    cs::BliInc
    is::BliInc

    # Bufferless scalar.
    scalar::BliAtomic

    # Pack-related fields.
    m_padded::BliDim
    n_padded::BliDim
    ps::BliInc
    pd::BliInc
    m_panel::BliDim
    n_panel::BliDim

    # Empty inner-constructor only,
    #  resembling a plain C struct.
    BliObjBase() = begin
        obj = new(Ptr{Nothing}(),
                  (0, 0),
                  (0, 0),
                  0, 0, 0, 0,
                  Ptr{Nothing}(), 0, 0, 0,
                  0,
                  0, 0, 0, 0,
                  0, 0)
    end
end

"""
Julia object BliObj packs obj_t with a reference to buffer array.
This is to preserve the base array from being recycled by GC module.
"""
struct BliObj
    obj::BliObjBase
    # Referee of obj.buffer, which itself is also a reference.
    # Purpose is to keep data alive.
    data::RefValue
end

" Object from strided matrix (view). "
BliObj(A::StridedMatrix{T}) where {T<:BliCompatibleType} = begin
    # Query data format.
    mA, nA = size(A)
    rsA, csA = strides(A)
    dtA = ctype_to_bli_num[T]

    # Initialize object with reference.
    objA = BliObj(BliObjBase(), Ref(A))
    # Pass address of A to object initializer.
    bli_obj_create_with_attached_buffer!(dtA, mA, nA,
                                         unsafe_convert(Ptr{Nothing}, A), rsA, csA,
                                         objA.obj)

    objA
end

" Object from strided vector. "
BliObj(x::StridedVector{T}) where {T<:BliCompatibleType} = begin
    mx, = size(x)
    incx, = strides(x)
    dt = ctype_to_bli_num[T]

    objx = BliObj(BliObjBase(), Ref(x))
    # Creation of vector buffer is roughly the same as matrix.
    bli_obj_create_with_attached_buffer!(dt, mx, 1,
                                         unsafe_convert(Ptr{Nothing}, x), incx, mx*incx,
                                         objx.obj)

    objx
end

" Object from scalar number. "
BliObj(γ::T) where {T<:BliCompatibleType} = begin
    dtγ = ctype_to_bli_num[T]

    objγ = BliObj(BliObjBase(), Ref(nothing))
    # bli_obj_create_1x1_with_attached_buffer!(dtγ,
    #                                          unsafe_convert(Ptr{Nothing}, objγ.data),
    #                                          objγ.obj)
    # Instead of something like above, use internal scalar storage.
    bli_obj_scalar_init_detached!(dtγ, objγ.obj)
    bli_setsc!(real(γ), imag(γ), objγ.obj)

    objγ
end

#----------------------------------------------------------
# Object-creation methods are exposed for the in-principle
#  consistency with BLIS C API and for advanced usages.
# Usually object constructors should be enough for Julia.
#
# NOTE: bli_obj_create, bli_obj_free, bli_obj_alloc_buffer,
#       bli_obj_create_1x1 and bli_obj_create_conf_to.
#       are not interfaced as Julia is a GC language.
#

bli_obj_create_without_buffer!(dt::BliNum,
                               m::BliDim,
                               n::BliDim,
                               obj::BliObjBase) = ccall(dlsym(libblis, :bli_obj_create_without_buffer),
                                                        Nothing,
                                                        (BliNum,
                                                         BliDim,
                                                         BliDim,
                                                         Ptr{BliObjBase}),
                                                        dt, m, n, Ref(obj))

bli_obj_attach_buffer!(p::Ptr{Nothing},
                       rs::BliInc,
                       cs::BliInc,
                       is::BliInc,
                       obj::BliObjBase) = ccall(dlsym(libblis, :bli_obj_attach_buffer),
                                                Nothing,
                                                (Ptr{Nothing},
                                                 BliInc,
                                                 BliInc,
                                                 BliInc,
                                                 Ptr{BliObjBase}),
                                                p, rs, cs, is, Ref(obj))

bli_obj_create_with_attached_buffer!(dt::BliNum,
                                     m::BliDim,
                                     n::BliDim,
                                     p::Ptr{Nothing},
                                     rs::BliInc,
                                     cs::BliInc,
                                     obj::BliObjBase) = ccall(dlsym(libblis, :bli_obj_create_with_attached_buffer),
                                                              Nothing,
                                                              (BliNum,
                                                               BliDim,
                                                               BliDim,
                                                               Ptr{Nothing},
                                                               BliInc,
                                                               BliInc,
                                                               Ptr{BliObjBase}),
                                                              dt, m, n, p, rs, cs, Ref(obj))

bli_obj_create_1x1_with_attached_buffer!(dt::BliNum,
                                         p::Ptr{Nothing},
                                         obj::BliObjBase) = ccall(dlsym(libblis, :bli_obj_create_1x1_with_attached_buffer),
                                                                  Nothing,
                                                                  (BliNum,
                                                                   Ptr{Nothing},
                                                                   Ptr{BliObjBase}),
                                                                  dt, p, Ref(obj))

bli_obj_scalar_init_detached!(dt::BliNum,
                              obj::BliObjBase) = ccall(dlsym(libblis, :bli_obj_scalar_init_detached),
                                                       Nothing,
                                                       (BliNum,
                                                        Ptr{BliObjBase}),
                                                       dt, Ref(obj))

bli_setsc!(zeta_r::Real,
           zeta_i::Real,
           obj::BliObjBase) = ccall(dlsym(libblis, :bli_setsc),
                                    Nothing,
                                    (Float64,
                                     Float64,
                                     Ptr{BliObjBase}),
                                    Float64(zeta_r), Float64(zeta_i), Ref(obj))

#----------------------------------------------------------------
# These methods replicates BLIS' C marcos for setting obj_t.info.
# They are exposed in a minimum fashion.

bli_obj_set_info_bits!(dest::BliObjBits,
                       src::BliObjBits,
                       obj::BliObjBase) = begin
    obj.info1 = UInt32( ( obj.info1 & ~dest ) | src )

    nothing
end

bli_obj_set_conjtrans!(t::BliTrans,
                       obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_CONJTRANS_BITS,
                                                                 BliObjBits(t.enum),
                                                                 obj)

bli_obj_set_onlytrans!(t::BliTrans,
                       obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_TRANS_BIT,
                                                                 BliObjBits(t.enum),
                                                                 obj)

bli_obj_set_conj!(c::BliConj,
                  obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_CONJ_BIT,
                                                            BliObjBits(c.enum),
                                                            obj)

bli_obj_set_uplo!(ul::BliUpLo,
                  obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_UPLO_BITS,
                                                            BliObjBits(ul.enum),
                                                            obj)

bli_obj_set_diag!(d::BliDiag,
                  obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_UNIT_DIAG_BIT,
                                                            BliObjBits(d.enum),
                                                            obj)

bli_obj_set_invert_diag!(id::BliInvDiag,
                         obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_INVERT_DIAG_BIT,
                                                                   BliObjBits(id.enum),
                                                                   obj)

bli_obj_set_struc!(st::BliStruc,
                   obj::BliObjBase) = bli_obj_set_info_bits!(BLIS_STRUC_BITS,
                                                             BliObjBits(st.enum),
                                                             obj)


text:
############################################################################################
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    x::Tuple
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
## Obtain flatten/unflatten constructor for each element
    x_vecs_and_backs = map(x) do xᵢ
        construct_flatten(T, flattentype, unflattentype, xᵢ)
    end
    _flatten, _unflatten = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
## Run flatten once to obtain dimensions of individual elements
    x_vecs = map(_flatten, x) do flat, xᵢ
         return flat(xᵢ)
    end
    lngth = map(length, x_vecs)
    sz = cumsum(lngth)
## Return flatten/unflatten for Tuple
    function construct_flatten_Tuple(x::Tuple)
        reduce(vcat,
            map(_flatten, x) do flat, xᵢ
                flat(xᵢ)
            end
        )
    end
    function construct_unflatten_Tuple(v::Union{R,AbstractVector{R}}) where {R<:Real}
        map(_unflatten, lngth, sz) do unflat, l, s
             return unflat(view(v, (s - l + 1):s))
        end
    end
    return construct_flatten_Tuple, construct_unflatten_Tuple
end

############################################################################################
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    constraint::Tuple,
    x::Tuple
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
## Obtain flatten/unflatten constructor for each element
    x_vecs_and_backs = map(x, constraint) do xᵢ, constraintᵢ
        construct_flatten(T, flattentype, unflattentype, constraintᵢ, xᵢ)
    end
    _flatten, _unflatten = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
## Run flatten once to obtain dimensions of individual elements
    x_vecs = map(_flatten, x) do flat, xᵢ
         return flat(xᵢ)
    end
    lngth = map(length, x_vecs)
    sz = cumsum(lngth)
## Return flatten/unflatten for Tuple
    function construct_flatten_Tuple(x::Tuple)
        reduce(vcat,
            map(_flatten, x) do flat, xᵢ
                flat(xᵢ)
            end
        )
    end
    function construct_unflatten_Tuple(v::Union{R,AbstractVector{R}}) where {R<:Real}
        map(_unflatten, lngth, sz) do unflat, l, s
             return unflat(view(v, (s - l + 1):s))
        end
    end
    return construct_flatten_Tuple, construct_unflatten_Tuple
end

############################################################################################
#!NOTE: If we map over collections, need to extend this to params.jl functions:
# _check, construct_transform, constrain, unconstrain, log_abs_det_jac

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::Tuple,
    x::Tuple
)
    return all(
            map(constraint, x) do cᵢ, xᵢ
                _check(_rng, cᵢ, xᵢ)
            end
    )
end

############################################################################################
function construct_transform(
    constraint::Tuple,
    x::Tuple
)
## Obtain transform/inversetransform constructor for each element
    x_transforms = map(x, constraint) do xᵢ, constraintᵢ
        construct_transform(constraintᵢ, xᵢ)
    end
    _transform, _inversetransform = first.(x_transforms), last.(x_transforms)
## Return flatten/unflatten for AbstractArray
    return _transform, _inversetransform
end

############################################################################################
function constrain(
    b⁻¹ᵥ::Tuple,
    θᵤ::Tuple
)
    return map(constrain, b⁻¹ᵥ, θᵤ)
end

############################################################################################
function unconstrain(
    bᵥ::Tuple,
    θ::Tuple
)
    return map(unconstrain, bᵥ, θ)
end

############################################################################################
function log_abs_det_jac(
    bᵥ::Tuple,
    θ::Tuple
)
    return sum(map(log_abs_det_jac, bᵥ, θ))
end

############################################################################################
#Export
export
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
"""
    FieldTimeSeries(filepath, name; architecture=CPU(), backend=InMemory())

Returns a `Dict` containing a `FieldTimeSeries` for each field in the JLD2 file located at `filepath`.
Note that model output must have been saved with halos. The `InMemory` backend will store the data
fully in memory as a 4D multi-dimensional array while the `OnDisk` backend will lazily load field time
snapshots when the `FieldTimeSeries` is indexed linearly.
"""
function FieldDataset(filepath; architecture=CPU(), backend=InMemory())
    file = jldopen(filepath)

    field_names = keys(file["timeseries"])
    filter!(k -> k != "t", field_names)  # Time is not a field.

    ds = Dict{String, FieldTimeSeries}(
        name => FieldTimeSeries(filepath, name; architecture, backend)
        for name in field_names
    )

    close(file)

    return ds
end


text:
import MPI
using PyCall
np=pyimport("numpy")

include("../types/types.jl")
include("../setting/constants.jl")
include("../utils/readfiles.jl")
include("../utils/kdtree.jl")
include("../utils/parse_commandline_ppm.jl")

function generate_profile_points(latnpts,lonnpts, vnpts, lon1, lat1, dep1, lon2, lat2, dep2,comm,latproc,lonproc)
    # MPI parameters
    rank = MPI.Comm_rank(comm)
    size_mpi = MPI.Comm_size(comm)

    # get ranges for the three directions
    rank_lat=rank%latproc+1
    rank_lon=div(rank,latproc)+1
    coor_lat=np[:array_split](1:latnpts,latproc)[rank_lat]
    coor_lon=np[:array_split](1:lonnpts,lonproc)[rank_lon]

    # init rθϕ_new
    ngll_new_this_rank=length(coor_lat)*length(coor_lon)*vnpts

    deplatlon_new_this_rank=zeros(Float64, 3, ngll_new_this_rank)
    # fill in deplatlon_new of this rank
    latnpts_this_rank=length(coor_lat)
    lonnpts_this_rank=length(coor_lon)
    for (vindex,dep) in enumerate(range(dep1,stop=dep2,length=vnpts))
        for (latindex,lat) in enumerate(range(lat1,stop=lat2,length=latnpts)[coor_lat])
            for (lonindex,lon) in enumerate(range(lon1,stop=lon2,length=lonnpts)[coor_lon])
                id=(vindex-1)*latnpts_this_rank*lonnpts_this_rank+(latindex-1)*lonnpts_this_rank+lonindex
                deplatlon_new_this_rank[:,id]=[dep,lat,lon]
            end
        end
    end

    xyz_new = zeros(Float64, 3, ngll_new_this_rank)
    # convert dep,lat,lon to x,y,z
    for id in 1:ngll_new_this_rank
        (dep,lat,lon)=deplatlon_new_this_rank[:,id]
        R_EARTH_KM = 6371.e0

        # normalized r
        r=(R_EARTH_KM-dep)/R_EARTH_KM
        θ=90-lat
        ϕ=lon

        # get x,y,z
        z=r*cosd(θ)
        h=r*sind(θ)
        x=h*cosd(ϕ)
        y=h*sind(ϕ)

        xyz_new[:,id]=[x,y,z]
    end

    return xyz_new
end

function run_interp( command_args::Dict{String,Any},comm::MPI.Comm)
    # MPI
    rank=MPI.Comm_rank(comm)
    isroot=(rank==0)

    # * init some variables
    nproc_old = command_args["nproc_old"]
    old_mesh_dir = command_args["old_mesh_dir"]
    old_model_dir = command_args["old_model_dir"]
    model_tags = command_args["model_tags"]
    output_file = command_args["output_file"]
    region=command_args["region"]
    npts=command_args["npts"]
    nproc=command_args["nproc"]
    (lon1,lat1,lon2,lat2,dep1,dep2)=parse.(Float64,split(region,"/"))
    (latnpts,lonnpts,vnpts)=parse.(Int64,split(npts,"/"))
    (latproc,lonproc)=parse.(Int64,split(nproc,"/"))

    # * init some variables
    mesh_old = sem_mesh_data()

    # * parse model tags
    model_names = String.(split(model_tags, ","))
    nmodel = length(model_names)

    if isroot
        @info "# nmodel=$(nmodel)"
        @info "# model names=$(model_names)"
    end

    # * interpolate old mesh/model onto new mesh

    # typical element size at surface for old mesh
    # read from constants.jl for old mesh
    typical_size = deg2rad(max(ANGULAR_WIDTH_XI_IN_DEGREES_VAL / NEX_XI_VAL, ANGULAR_WIDTH_ETA_IN_DEGREES_VAL / NEX_ETA_VAL)) * R_UNIT_SPHERE
    max_search_dist = 10.0 * typical_size
    max_misloc = typical_size / 4.0

    # get xyz_new
    xyz_new=generate_profile_points(latnpts,lonnpts, vnpts, lon1, lat1, dep1, lon2, lat2, dep2,comm,latproc,lonproc)

    ngll_new_this_rank = size(xyz_new)[2]
    idoubling_new = zeros(Int64, ngll_new_this_rank)
    idoubling_new.=IFLAG_DUMMY

    # * initialize variables for interpolation
    location_1slice = Array{sem_mesh_location}(undef, ngll_new_this_rank)
    stat_final = zeros(Int64, ngll_new_this_rank)
    misloc_final = zeros(Float64, ngll_new_this_rank)

    stat_final .= -1
    misloc_final .= HUGEVAL

    model_interp_this_rank = zeros(nmodel, ngll_new_this_rank)

    # * loop each slices of the old mesh
    flag = true

    # * loop for all points in xyz_new
    for iproc_old = 0:nproc_old - 1
        @info "[$rank]# iproc_old=$(iproc_old)"
        # read old mesh slice
        mesh_old = sem_mesh_read(old_mesh_dir, iproc_old)
        nspec_old = mesh_old.nspec
        # test if the new and old mesh slices are separated apart
        min_dist = HUGEVAL
        for ispec = 1:nspec_old
            iglob = mesh_old.ibool[MIDX,MIDY,MIDZ,ispec]
            old_x=mesh_old.xyz_glob[1,iglob]
            old_y=mesh_old.xyz_glob[2,iglob]
            old_z=mesh_old.xyz_glob[3,iglob]
            dist_this_spec=sqrt(minimum(@. (xyz_new[1,:]-old_x)^2+(xyz_new[2,:]-old_y)^2+(xyz_new[3,:]-old_z)^2))
            min_dist=min(min_dist,dist_this_spec)
        end
        if min_dist>max_search_dist
            @info "[$rank]#$(iproc_old) slices too far away, skip"
            continue
        end

        # read old model
        model_gll_old = zeros(Float64, nmodel, NGLLX, NGLLY, NGLLZ, nspec_old)
        sem_io_read_gll_file_n!(old_model_dir, iproc_old, model_names, nmodel, model_gll_old)

        # locate points in this mesh slice
        nnearest = 10
        location_1slice = sem_mesh_locate_kdtree2!(mesh_old, ngll_new_this_rank, xyz_new, idoubling_new, nnearest, max_search_dist, max_misloc, iproc_old)

        for igll = 1:ngll_new_this_rank
            if (stat_final[igll] == 1 && location_1slice[igll].stat == 1)
                @info "[$rank]# multi-located, $(xyz_new[:,igll])"
                continue
            end
            # for point located inside one element in the first time or closer to one element than located before
            if location_1slice[igll].stat == 1 || (location_1slice[igll].stat == 0 && location_1slice[igll].misloc < misloc_final[igll])
                for imodel = 1:nmodel
                    model_interp_this_rank[imodel,igll] = sum(location_1slice[igll].lagrange .* model_gll_old[imodel,:,:,:,location_1slice[igll].eid])
                end
                stat_final[igll] = location_1slice[igll].stat
                misloc_final[igll] = location_1slice[igll].misloc
            end
        end
    end

    # gather all ngll_new_this_rank to ngll_new

    # * write out gll files for this new mesh slice
    # get latnpts_this_rank and lonnpts_this_rank
    rank = MPI.Comm_rank(comm)
    size_mpi = MPI.Comm_size(comm)

    rank_lat=rank%latproc+1
    rank_lon=div(rank,latproc)+1
    coor_lat=np[:array_split](1:latnpts,latproc)[rank_lat]
    coor_lon=np[:array_split](1:lonnpts,lonproc)[rank_lon]

    latnpts_this_rank=length(coor_lat)
    lonnpts_this_rank=length(coor_lon)

    # output files
    if isroot
        run(`rm -rf $output_file`)
        run(`mkdir -p $output_file`)
    end
    MPI.Barrier(comm)

    open(output_file*"/$rank","w") do io
        for (vindex,dep) in enumerate(range(dep1,stop=dep2,length=vnpts))
            for (latindex,lat) in enumerate(range(lat1,stop=lat2,length=latnpts)[coor_lat])
                for (lonindex,lon) in enumerate(range(lon1,stop=lon2,length=lonnpts)[coor_lon])
                    id=(vindex-1)*latnpts_this_rank*lonnpts_this_rank+(latindex-1)*lonnpts_this_rank+lonindex
                    thesize=size(model_interp_this_rank)[1]
                    if thesize==1
                        write(io,"$lon $lat $dep $(model_interp_this_rank[1,id]) \n")
                    else
                        out_string="$lon $lat $dep "
                        for item in model_interp_this_rank[:,id]
                            out_string*="$item "
                        end
                        out_string*="\n"
                        write(io,out_string)
                    end
                end
            end
        end
    end

end

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD

    command_args = parse_commandline()
    run_interp(command_args,comm)

    MPI.Finalize()
end

main()

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
isinstalled(a) = isa(Pkg.installed(a), VersionNumber)
if isinstalled("PyPlot")
	using PyPlot

    function axisij()
        a = axis()
        axis((a[1],a[2],a[4],a[3]))
    end

	function plotshape(a, args...; kargs...)
		if size(a,1)==2
			plot(vec(a[2,:]), vec(a[1,:]), args...; kargs...)
            axisij()
			xlabel("N")
			ylabel("M")
		elseif size(a,1)==3
			plot3D(vec(a[1,:]), vec(a[2,:]), vec(a[3,:]), args...; kargs...)
			xlabel("M")
			ylabel("N")
			zlabel("O")
		else
			error("Can't plot data of size $(size(a))")
		end
		axis("equal")
	end

	function plotshapes(a)
		if size(a,1)==3
			for i = 1:last(size(a))
				plotshape(a[:,:,i])
				hold(true)
			end
			hold(false)
		else
			gridplot(a)
		end
	end

	at(a,i) = slicedim(a,ndims(a),i)

	function gridplot(a)
		N = last(size(a))
		sm = floor(sqrt(N))
		sn = ceil(N/sm)
		for m = 1:sm, n = 1:sn
			datai = (n-1)*sm+m
			ploti = (m-1)*sn+n
			if datai <= N
				subplot(sm,sn,ploti)
				plotshape(at(a,round(Int,datai)))
			end
		end
	end

end


text:
mutable struct MarkedQueued
    queued::Bool
end
MarkedQueued() = MarkedQueued(false)
CachedArrays.markqueued!(x::MarkedQueued) = (x.queued = true)
CachedArrays.isqueued(x::MarkedQueued) = (x.queued)

@testset "Testing FreeBuffer" begin
    buffer = CachedArrays.FreeBuffer{MarkedQueued}()
    @test CachedArrays.candrain(buffer) == false

    a = MarkedQueued()
    @test CachedArrays.isqueued(a) == false
    push!(buffer, a)
    @test CachedArrays.isqueued(a) == true
    @test CachedArrays.candrain(buffer) == true
    CachedArrays.unsafe_swap!(buffer)

    buffer = CachedArrays.unsafe_get(buffer)
    @test buffer == [a]
end


text:
using QuantumStateBase
using Test
using LinearAlgebra
using ClassicalOrthogonalPolynomials

const QSB = QuantumStateBase
const DIM = QSB.DIM

Base.:(==)(s1::StateVector, s2::StateVector) = (s1.v == s2.v) && (s1.dim == s2.dim)
Base.:(==)(s1::StateMatrix, s2::StateMatrix) = (s1.𝛒 == s2.𝛒) && (s1.dim == s2.dim)
Base.:(≈)(s1::StateVector, s2::StateVector) = (s1.v ≈ s2.v) && (s1.dim == s2.dim)
Base.:(≈)(s1::StateMatrix, s2::StateMatrix) = (s1.𝛒 ≈ s2.𝛒) && (s1.dim == s2.dim)

@testset "QuantumStateBase.jl" begin
    # polynomial
    include("polynomial.jl")

    # state
    include("representation.jl")
    include("basis.jl")
    include("operator.jl")
    include("state.jl")

    # wigner
    include("wigner_util.jl")
    include("wigner.jl")

    # pdf of quadrature
    # include("quadrature_pdf.jl")
    # include("sampler.jl")
end


text:
abstract type ELMTypes end
struct Regression <: ELMTypes end
struct Classification <: ELMTypes
    nclass::Int
end

mutable struct ELM{T<:BaseType}
    W::Matrix{T}
    b::Matrix{T}
    beta::Matrix{T}
    L::Int
    mf::MappingFunctions
    elmtype::ELMTypes

    function ELM{T}(L::Int;
                 mf::MappingFunctions=Sigmoid(),
                 elmtype::ELMTypes=Regression()) where T<:BaseType
        net = new()
        net.L = L
        net.mf = mf
        net.elmtype = elmtype
        return net
    end
end

function train!(net::ELM{T}, X::AbstractMatrix{T}, Y::AbstractMatrix{S};
                tm::TrainingMethods=Batch(),
                init::Initializers=Uniform()) where {T<:BaseType, S<:Union{BaseType,Integer}}
    @assert(size(Y,2) == size(X,2))

    net.W = rand_arr(T, (net.L, size(X,1)), init)
    net.b = rand_arr(T, (net.L, 1), init)

    X, Y = preprocess(X, Y, net.elmtype)
    net.beta = train_core(net.W, net.b, X, Y, net.mf, tm)

    return net
end

function predict(net::ELM{T}, X::AbstractMatrix{T}) where {T<:BaseType}
    H = mapping(net.W, net.b, X, net.mf)
    Y = net.beta*H
    return postprocess(Y, net.elmtype)
end

function preprocess(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, elmtype::Regression) where {T<:BaseType}
    return X, Y
end

function postprocess(Y::AbstractMatrix{T}, elmtype::Regression) where {T<:BaseType}
    return Y
end

function preprocess(X::AbstractMatrix{T}, Y::AbstractMatrix{S}, elmtype::Classification) where {T<:BaseType, S<:Integer}
    @assert(maximum(Y) < elmtype.nclass, "Class labels must be in {0,1,...,nclass-1}.")
    if elmtype.nclass == 2
        Yc = zeros(T, 1, size(Y,2))
        for i=1:size(Y,2)
            Yc[1,i] = Y[1,i] > zero(S) ? one(T) : -one(T)
        end
    else
        Yc = zeros(T, elmtype.nclass, size(Y,2))
        for i=1:size(Y,2)
            Yc[Y[1,i]+one(S),i] = one(T)
        end
    end
    return X, Yc
end

function postprocess(Yc::AbstractMatrix{T}, elmtype::Classification) where {T<:BaseType}
    Y = zeros(Int, 1, size(Yc,2))
    if elmtype.nclass == 2
        for i=1:size(Yc,2)
            Y[1,i] = real(Yc[i]) > zero(Int) ? one(Int) : zero(Int)
        end
    else
        for i=1:size(Yc,2)
            Y[1,i] = indmax(real(Yc[:,i])) - one(Int)
        end
    end
    return Y
end