text:
#----------------------------------------------------------------------------------------------------------------------------------------
# Implements solution to the diffusion equation in a N-layered finite cylinder as given in [1].
# Solutions are given in the spatial, frequency, and time domains for a point source located in the middle of the cylinder top.
#
# [1] André Liemert and Alwin Kienle, "Light diffusion in a turbid cylinder. II. Layered case," Opt. Express 18, 9266-9279 (2010)
#----------------------------------------------------------------------------------------------------------------------------------------

@with_kw struct Nlayer_cylinder{T <: Real} <: DiffusionParameters
    μsp::Vector{T} = [10.0, 10.0, 10.0, 10.0]               # reduced scattering coefficient (1/cm)
    μa::Vector{T} = [0.1, 0.1, 0.1, 0.1]                    # absorption coefficient (1/cm)
    n_ext::T = 1.0                                          # surrounding index of refraction
    n_med::Vector{T} = [1.0, 1.0, 1.0, 1.0]                 # layers index of refraction

    l::Vector{T} = [0.5, 0.8, 1.0, 5.0]                     # length of cylinder layers (cm)
    ρ::Union{T, AbstractVector{T}} = 1.0                    # source-detector separation (cm)
    a::T = 5.0                                              # radius of cylinder (cm)
    z::T = 0.0                                              # detector depth (cm)

    ω::T = 0.0                                              # modulation frequency
    bessels::AbstractVector{T} = besselroots[1:1000]        # roots of J0
end

#-------------------------------------------
# Steady-State Fluence
#-------------------------------------------
"""
    fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots)

Compute the steady-state fluence in an N-layered cylinder. Source is assumed to be located on the top middle of the cylinder.
No checks on arguments are made except enforcing that the source is located in the top layer.
ρ, μa, and z should be >= 0.0.
μsp, n_ext, n_med, a should be > 0.0.
Each layer in l should be > 0.0.
In general, ρ should always be < a.
All arguments should have the same types and the length of μsp, μa, n_med, l should be equal.
It is advised to use the Nlayer_cylinder() structure format for inputs.

# Arguments
- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: layer thicknesses (cm)
- `a`: cylinder radius (cm)
- `z`: source depth within cylinder
- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0

# Examples
julia> `fluence_DA_Nlay_cylinder_CW(1.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)`
"""
function fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots)
    D = D_coeff.(μsp)
    N = length(D)
    A = A_coeff.(n_med / n_ext)
    z0 = z0_coeff(μsp[1])
    zb = zb_coeff.(A, D)
    @assert z0 < l[1]

    if z < l[1]
        return _kernel_fluence_DA_Nlay_cylinder(ρ, D, μa, a, zb, z, z0, l, n_med, besselroots, _green_Nlaycylin_top, N) / (π * (a + zb[1])^2)
    elseif z > sum(l[1:end - 1])
        return _kernel_fluence_DA_Nlay_cylinder(ρ, D, μa, a, zb, z, z0, l, n_med, besselroots, _green_Nlaycylin_bottom, N) / (π * (a + zb[1])^2)
    end
end
"""
    fluence_DA_Nlay_cylinder_CW(data)

Wrapper to fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots) with inputs given as a structure (data).

# Examples
julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0)
julia> `fluence_DA_Nlay_cylinder_CW(data)`
"""
function fluence_DA_Nlay_cylinder_CW(data)
    return fluence_DA_Nlay_cylinder_CW(data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, data.bessels)
end

#-------------------------------------------
# Steady-State Flux
#-------------------------------------------
"""
    flux_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots)

Compute the steady-state flux using Fick's law D[1]*∂ϕ(ρ)/∂z for z = 0 (reflectance) and -D[end]*∂ϕ(ρ)/∂z for z = sum(l) (transmittance) in an N-layered cylinder.
Source is assumed to be located on the top middle of the cylinder. z must be equal to 0 or the total length sum(l) of cylinder.

# Arguments
- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: layer thicknesses (cm)
- `a`: cylinder radius (cm)
- `z`: source depth within cylinder: must be equal to 0 or sum(l)
- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0

# Examples
julia> `flux_DA_Nlay_cylinder_CW(1.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)`
"""
function flux_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots)
    @assert z == zero(eltype(z)) || z == sum(l)
    D = D_coeff.(μsp)
    if z == zero(eltype(z))
        return D[1] * ForwardDiff.derivative(dz -> fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, dz, besselroots), z)
    elseif z == sum(l)
        return -D[end] * ForwardDiff.derivative(dz -> fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, dz, besselroots), z)
    end
end
"""
    flux_DA_Nlay_cylinder_CW(data, besselroots)

Wrapper to flux_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots) with inputs given as a structure (data).

# Examples
julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0)
julia> `flux_DA_Nlay_cylinder_CW(data, besselroots)`
"""
function flux_DA_Nlay_cylinder_CW(data)
    return flux_DA_Nlay_cylinder_CW(data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, data.bessels)
end

#-------------------------------------------
# Time-Domain Fluence
#-------------------------------------------
"""
    fluence_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24, ILT = hyper_fixed)

Compute the time-domain fluence in an N-layered cylinder. Source is assumed to be located on the top middle of the cylinder.
See notes on fluence_DA_Nlay_cylinder_CW for input checks with the initial constraint that all values of t should be > 0.0 and ordered least to greatest.
Do not have any value of t be zero or negative.

It is best to use `hyper_fixed` as the inverse laplace transform if t consists of many time points. Utilize `hyperbola` for a single time point.
The value of N should be proportional to the dynamic range of the time-domain signal needed. For later times you will need a larger N.
The lowest fluence value you can compute will be no less than the machine precision you utilize. Float64 values are limited to fluence > ~2.2-16.

# Arguments
- `t`: time point or vector (ns)
- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: layer thicknesses (cm)
- `a`: cylinder radius (cm)
- `z`: source depth within cylinder
- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0
- `N`: number of Hankel-Laplace calculations
- `ILT`: inverse laplace transform function

# Examples
julia> `fluence_DA_Nlay_cylinder_TD(0.1:0.1:2.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)`
"""
function fluence_DA_Nlay_cylinder_TD(t::AbstractFloat, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 18, ILT = hyperbola)
    return ILT(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N)
end
function fluence_DA_Nlay_cylinder_TD(t::AbstractArray, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24, ILT = hyper_fixed)
    T = promote_type(eltype(ρ), eltype(μa), eltype(μsp), eltype(z), eltype(l))
    return ILT(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N, T = T)
end
"""
    fluence_DA_Nlay_cylinder_TD(t, data; N = 24)

Wrapper to fluence_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24, ILT = hyper_fixed) with inputs given as a structure (data).

# Examples
julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0)
julia> `fluence_DA_Nlay_cylinder_TD(0.5:1.0:2.5, data)`
"""
function fluence_DA_Nlay_cylinder_TD(t, data; N = 24)
    return fluence_DA_Nlay_cylinder_TD(t, data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, data.bessels, N = N)
end

#-------------------------------------------
# Time-Domain Flux
#-------------------------------------------
"""
    flux_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24)

Compute the time-domain flux using Fick's law D[1]*∂ϕ(ρ, t)/∂z for z = 0 (reflectance) and -D[end]*∂ϕ(ρ, t)/∂z for z = sum(l) (transmittance) in an N-layered cylinder.
Source is assumed to be located on the top middle of the cylinder. z must be equal to 0 or the total length sum(l) of cylinder.
Uses the hyperbola contour if t is an AbstractFloat and the fixed hyperbola contour if t is an AbstractVector.
∂ϕ(ρ, t)/∂z is calculated using forward mode auto-differentiation with ForwardDiff.jl

# Arguments
- `t`: time point or vector (ns)
- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: layer thicknesses (cm)
- `a`: cylinder radius (cm)
- `z`: source depth within cylinder
- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0
- `N`: number of Hankel-Laplace calculations

# Examples
julia> `fluence_DA_Nlay_cylinder_TD(0.1:0.1:2.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, besselroots)`
"""
function flux_DA_Nlay_cylinder_TD(t::AbstractVector, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24)
    @assert z == zero(eltype(z)) || z == sum(l)
    D = D_coeff.(μsp)

    ## need to know the type of z to preallocate vectors coorectly in hyper_fixed for autodiff
    function _ILT(z::T) where {T}
        return hyper_fixed(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots), t, N = N, T = T)
    end
    if z == zero(eltype(z))
        return D[1]*ForwardDiff.derivative(_ILT, z)
    elseif z == sum(l)
        return -D[end]*ForwardDiff.derivative(_ILT, z)
    end
end
function flux_DA_Nlay_cylinder_TD(t::AbstractFloat, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24)
    @assert z == zero(eltype(z)) || z == sum(l)
    D = D_coeff.(μsp)
    if z == zero(eltype(z))
        return D[1] * ForwardDiff.derivative(dz -> hyperbola(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, dz, s, besselroots), t, N = N), 0.0)
    elseif z == sum(l)
        return -D[end] * ForwardDiff.derivative(dz -> hyperbola(s -> _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, dz, s, besselroots), t, N = N), sum(l))
    end
end

"""
    flux_DA_Nlay_cylinder_TD(t, data; N = 24)

Wrapper to flux_DA_Nlay_cylinder_TD(t, ρ, μa, μsp, n_ext, n_med, l, a, z, besselroots; N = 24) with inputs given as a structure (data).

# Examples
julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0)
julia> `flux_DA_Nlay_cylinder_TD(0.5:1.0:2.5, data)`
"""
function flux_DA_Nlay_cylinder_TD(t, data;  N = 24)
    return flux_DA_Nlay_cylinder_TD(t, data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, data.bessels, N = N)
end

#-------------------------------------------
# Frequency-Domain Fluence
#-------------------------------------------
"""
    fluence_DA_Nlay_cylinder_FD(ρ, μa, μsp, n_ext, n_med, l, a, z, ω, besselroots)

Compute the frequency modulated fluence in an N-layered cylinder. Source is assumed to be located on the top middle of the cylinder.

# Arguments
- `ρ`: source-detector separation in cylindrical coordinates (distance from middle z-axis of cylinder) (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `l`: layer thicknesses (cm)
- `a`: cylinder radius (cm)
- `z`: source depth within cylinder
- `ω`: source modulation frequency (1/ns)
- `besselroots`: roots of bessel function of first kind zero order J0(x) = 0

# Examples
julia> `fluence_DA_Nlay_cylinder_CW(1.0, [0.1, 0.1], [10.0, 10.0], 1.0, [1.0, 1.0], [4.5, 4.5], 10.0, 0.0, 1.0, besselroots)`
"""
function fluence_DA_Nlay_cylinder_FD(ρ, μa, μsp, n_ext, n_med, l, a, z, ω, besselroots)
    ν = ν_coeff.(n_med)
    μa_complex = μa .+ ω * im ./ ν

    return fluence_DA_Nlay_cylinder_CW(ρ, μa_complex, μsp, n_ext, n_med, l, a, z, besselroots)
end
"""
    fluence_DA_Nlay_cylinder_FD(data)

Wrapper to fluence_DA_Nlay_cylinder_FD(ρ, μa, μsp, n_ext, n_med, l, a, z, ω, besselroots) with inputs given as a structure (data).

# Examples
julia> data = Nlayer_cylinder(a = 10.0, l = [1.0, 1.0, 1.0, 2.0], z = 5.0, ω = 1.0)
julia> `fluence_DA_Nlay_cylinder_FD(data)`
"""
function fluence_DA_Nlay_cylinder_FD(data)
    return fluence_DA_Nlay_cylinder_FD(data.ρ, data.μa, data.μsp, data.n_ext, data.n_med, data.l, data.a, data.z, data.ω, data.bessels)
end

#-------------------------------------------------------------------------------
# this function is the base for the laplace transform in the time-domain
#-------------------------------------------------------------------------------
function _fluence_DA_Nlay_cylinder_Laplace(ρ, μa, μsp, n_ext, n_med, l, a, z, s, besselroots)
    ν = ν_coeff.(n_med)
    μa_complex = μa .+ s ./ ν

    return fluence_DA_Nlay_cylinder_CW(ρ, μa_complex, μsp, n_ext, n_med, l, a, z, besselroots)
end

#-------------------------------------------------------------------------------
# _kernel_fluence_DA_Nlay_cylinder provides the main kernel of the program.
# The inputs are similar to fluence_DA_Nlay_cylinder_CW.
# D is the diffusion coefficient and N is the number of layers.
# green is the Green's function for either the first or bottom layer (below).
#-------------------------------------------------------------------------------
function _kernel_fluence_DA_Nlay_cylinder(ρ::AbstractFloat, D, μa, a, zb, z, z0, l, n_med, besselroots, green, N)
    ϕ = zero(eltype(μa))
    ϕ_tmp = zero(eltype(μa))
    α = zeros(eltype(μa), N)

    for ind in eachindex(besselroots)
        ϕ_tmp = green(α, besselroots[ind] / (a + zb[1]), μa, D, z, z0, zb, l, n_med, N)
        ϕ_tmp *= besselj0(besselroots[ind] / (a + zb[1]) * ρ)
        ϕ_tmp /= (besselj1(besselroots[ind]))^2
        ϕ += ϕ_tmp
    end

    return ϕ
end

function _kernel_fluence_DA_Nlay_cylinder(ρ::AbstractVector, D, μa, a, zb, z, z0, l, n_med, besselroots, green, N)
    ϕ = zeros(eltype(μa), length(ρ))
    ϕ_tmp = zero(eltype(μa))
    α = zeros(eltype(μa), N)

    for ind in eachindex(besselroots)
        tmp = besselroots[ind] / (a + zb[1])
        ϕ_tmp = green(α, tmp, μa, D, z, z0, zb, l, n_med, N)
        ϕ_tmp /= (besselj1(besselroots[ind]))^2
        for ρ_ind in eachindex(ρ)
            ϕ[ρ_ind] += ϕ_tmp * besselj0(tmp * ρ[ρ_ind])
        end
    end

    return ϕ
end
#-------------------------------------------------------------------------------
# Calculates the Green's function in the first (top) and last (bottom) layer
# sinh and cosh have been expanded as exponentials.
# A common term has been factored out. This cancels for G1 but not GN (see _βγN_correction)
# For N = 2, 3, 4 coefficients are explicitly calculated.
# For N > 4, β and γ are calculated recursively using eqn. 17 & 18.
#-------------------------------------------------------------------------------
@inline function α_coeff!(α, μa, D, sn)
    @inbounds for ind in 1:length(μa)
        α[ind] = sqrt(μa[ind] / D[ind] + sn^2)
    end
    return α
end
@inline function _green_Nlaycylin_top(α, sn, μa, D, z, z0, zb, l, n, N)
    α = α_coeff!(α, μa, D, sn)

    if N == 4
        β, γ = _get_βγ4(α, D, n, zb, l)
    elseif N == 3
        β, γ = _get_βγ3(α, D, n, zb, l)
    elseif N == 2
        β, γ = _get_βγ2(α, D, n, zb, l)
    elseif N > 4
        β, γ = _get_βγk(α, D, n, zb, l)
    end

    tmp1 = D[1] * α[1] * n[1]^2 * β
    tmp2 = D[2] * α[2] * n[2]^2 * γ
    tmp3 = exp(-2 * α[1] * (l[1] + zb[1]))

    g  = (exp(-α[1] * abs(z - z0)) - exp(-α[1] * (z + z0 + 2 * zb[1])))
    g1 = exp(α[1] * (z + z0 - 2 * l[1])) * (1 - exp(-2 * α[1] * (z0 + zb[1]))) * (1 - exp(-2 * α[1] * (z + zb[1])))
    g1 *= (tmp1 - tmp2)
    g1 /= (tmp1 * (1 + tmp3) + tmp2 * (1 - tmp3))

    return (g + g1) / (2 * D[1] * α[1])
end
@inline function _green_Nlaycylin_bottom(α, sn, μa, D, z, z0, zb, l, n, N)
    α = α_coeff!(α, μa, D, sn)

    if N == 4
        β, γ = _get_βγ4(α, D, n, zb, l)
        βγ_correction = _βγ4_correction(α, zb, l)
    elseif N == 3
        β, γ = _get_βγ3(α, D, n, zb, l)
        βγ_correction = _βγ3_correction(α, zb, l)
    elseif N == 2
        β, γ = _get_βγ2(α, D, n, zb, l)
        βγ_correction = _βγ2_correction(α, zb, l)
    elseif N > 4
        β, γ = _get_βγk(α, D, n, zb, l)
        βγ_correction = _βγN_correction(α, zb, l)
    end

    tmp = 1.0
    for ind in (N - 1):-1:2
        tmp *= D[ind] * α[ind] * n[ind]^2
    end

    gN = n[end]^2 * tmp * 2^(N - 1) / 2
    gN *= exp(α[1] * (z0 - l[1]) + α[end] * (sum(l) + zb[end] - z) - βγ_correction)
    gN /= D[1] * α[1] * n[1]^2 * β * (1 + exp(-2 * α[1] * (l[1] + zb[1]))) + D[2] * α[2] * n[2]^2 * γ * (1 - exp(-2 * α[1] * (l[1] + zb[1])))
    gN *= (1 - exp(-2 * α[1] * (z0 + zb[1]))) * (1 - exp(-2 * α[end] * (sum(l) + zb[end] - z)))

    return gN
 end

#-------------------------------------------------------------------------------
# Calculate β and γ coefficients with eqn. 17 in [1].
# sinh and cosh have been expanded as exponentials.
# A common term has been factored out. This cancels for G1 but not GN (see _βγN_correction)
# For N = 2, 3, 4 coefficients are explicitly calculated.
# For N > 4, β and γ are calculated recursively using eqn. 17 & 18.
#-------------------------------------------------------------------------------
@inline function _get_βγ2(α, D, n, zb, l)
    tmp1 = exp(-2 * α[2] * (l[2] + zb[2]))

    β = (1 - tmp1)
    γ = (1 + tmp1)

    return β, γ
end
@inline function _get_βγ3(α, D, n, zb, l)
    tmp1 = D[2] * α[2] * n[2]^2
    tmp2 = D[3] * α[3] * n[3]^2
    tmp3 = exp(-2 * α[2] * l[2])
    tmp4 = exp(-2 * α[3] * (l[3] + zb[2]))

    β  = tmp1 * (1 + tmp3) * (1 - tmp4)
    β += tmp2 * (1 - tmp3) * (1 + tmp4)

    γ  = tmp1 * (1 - tmp3) * (1 - tmp4)
    γ += tmp2 * (1 + tmp3) * (1 + tmp4)

    return β, γ
end
@inline function _get_βγ4(α, D, n, zb, l)
    tmp5 = D[2] * α[2] * n[2]^2
    tmp1 = D[3] * α[3] * n[3]^2
    tmp4 = D[4] * α[4] * n[4]^2

    tmp6 = exp(-2 * α[2] * l[2])
    tmp2 = exp(-2 * α[3] * l[3])
    tmp3 = exp(-2 * α[4] * (l[4] + zb[2]))

    β_4  = tmp1 * (1 + tmp2) * (1 - tmp3)
    β_4 += tmp4 * (1 - tmp2) * (1 + tmp3)

    γ_4  = tmp1 * (1 - tmp2) * (1 - tmp3)
    γ_4 += tmp4 * (1 + tmp2) * (1 + tmp3)

    β = (tmp5 * β_4 * (1 + tmp6) + tmp1 * γ_4 * (1 - tmp6))
    γ = (tmp5 * β_4 * (1 - tmp6) + tmp1 * γ_4 * (1 + tmp6))

    return β, γ
end
@inline function _get_βγk(α, D, n, zb, l)
    βN, γN = _get_βγN(α, D, n, zb, l)

    for k in length(α):-1:4
        tmp1 = D[k - 2] * α[k - 2] * n[k - 2]^2 * βN
        tmp2 = D[k - 1] * α[k - 1] * n[k - 1]^2 * γN
        tmp3 = exp(-2 * α[k - 2] * l[k - 2])

        βN  =  tmp1 * (1 + tmp3)
        βN +=  tmp2 * (1 - tmp3)
        γN  =  tmp1 * (1 - tmp3)
        γN +=  tmp2 * (1 + tmp3)

    end

    return βN, γN
end
@inline function _get_βγN(α, D, n, zb, l)
    tmp1 = D[end - 1] * α[end - 1] * n[end - 1]^2
    tmp2 = D[end] * α[end] * n[end]^2
    tmp3 = exp(-2 * α[end - 1] * l[end - 1])
    tmp4 = exp(-2 * α[end] * (l[end] + zb[2]))

    βN  =   tmp1 * (1 + tmp3) *  (1 - tmp4)
    βN +=   tmp2 * (1 - tmp3) *  (1 + tmp4)

    γN =  tmp1 * (1 - tmp3) *  (1 - tmp4)
    γN +=  tmp2 * (1 + tmp3) *  (1 + tmp4)

    return βN, γN
end

#-------------------------------------------------------------------------------
# Calculate βγ correction for N layer green's function.
# This is done because a common term was factored out from general expressions
# They don't cancel for the N layer so have to divide them out.
# For N = 2, 3, 4 coefficients are explicitly calculated.
# For N > 4, β and γ are calculated recursively
#-------------------------------------------------------------------------------
@inline function _βγ2_correction(α, zb, l)
    return α[2] * (l[2] + zb[2])
end
function _βγ3_correction(α, zb, l)
    return α[2] * l[2] + α[3] * (l[3] + zb[2])
end
@inline function _βγ4_correction(α, zb, l)
    return α[2] * l[2] + α[3] * l[3] + α[4] * (l[4] + zb[2])
end
@inline function _βγ5_correction(α, zb, l)
    return α[2] * l[2] + α[3] * l[3] + α[4] * l[4] + α[5] * (l[5] + zb[2])
end
@inline function _βγN_correction(α, zb, l)
    out = α[end] * (l[end] + zb[2])
    for ind in (length(α) - 1):-1:2
        out += α[ind] * l[ind]
    end
    return out
end
#-------------------------------------------------------------------------------

max_stars_repo_name: heltonmc/LightPropagation

max_issues_repo_name: heltonmc/LightPropagation


text:
# ASuMa May, 2018
# This script takes a directory with corpus files and one or more
# AdaGram model(s) to annotate the corpus.
# It outputs the corpus files annotated with word senses, if there is more
# than one sense above the threshold for a given word

# usage:
#   annotate_corpus model(file or dir) corpus-dir
#                   output-dir [--joiner=@] [--min-prob=0.3] [--window=4]
# see arg table for meaning of parameters

#push!(LOAD_PATH, "./src/")

function annotate_word(outfile, separator, threshold, vm, dict, word, context)
"""
    Writes word in the outfile.
    If the word has more than one sense above threshold in vm model,
    its context is used to get the current sense, and
    it's annotated with separator and sense number
"""
    letters = "abcdefghjklmnopqrstuvwxyz" # "i" skipped on purpose bc LG dict
    output_word = word
    id = get(dict.word2id, word, -1)
    if id != -1
        priors = expected_pi(vm, id)
        if length(find(priors .> threshold)) > 1
            probs = disambiguate(vm, dict, word, context, true, threshold)
            best_sense = findmax(probs)[2]
            tag = string(letters[best_sense % length(letters)])
            output_word = word * separator * tag
        end
    end
    write(outfile, output_word * " ")
end

function annotate_file(corpus, outfile, vm, dict, separator, min_prob, win)
"""
   Goes through lines in corpus and calls annotate_word to re-write corpus
"""
    open(corpus, "r") do fi
        open(outfile, "w") do fo
            for line in eachline(fi)
                if line == "\n"
                    write(fo, "\n")
                end
                split_line = split(line)
                #println(split_line)
                for i in enumerate(split_line)
                    llim = max(1, i[1] - win)
                    ulim = min(length(split_line), i[1] + win)
                    context = split_line[llim:ulim]
                    deleteat!(context, i[1] + 1 - llim)
                    rm_index = []
                    for word in enumerate(context)
                        if get(dict.word2id, word[2], -1) == -1
                            push!(rm_index, word[1])
                        end
                    end
                    deleteat!(context, rm_index)
                    #println(context)
                    annotate_word(fo, separator, min_prob, vm, dict, i[2], context)
                end
                seek(fo, position(fo) - 1)
                write(fo, "\n")
            end
        end
    end
end

using ArgParse
using AdaGram

s = ArgParseSettings()

@add_arg_table s begin
  "model"
    help = "File or directory with AdaGram model(s)"
    arg_type = AbstractString
    required = true
  "corpus-dir"
    help = "Directory with files to annotate"
    arg_type = AbstractString
    required = true
  "output-dir"
    help = "Directory to output annotated corpus"
    arg_type = AbstractString
    required = true
  "--joiner"
    help = "Symbol to annotate word senses"
    arg_type = AbstractString
    default = "@"
  "--min-prob"
    help = "Min probability of second sense to consider a word ambiguous"
    arg_type = Float64
    default = 0.3
  "--window"
    help = "Size of window to look for context"
    arg_type = Int64
    default = 4
end

args = parse_args(ARGS, s)

separator = args["joiner"]
min_prob = args["min-prob"]
win = args["window"]
corpus_dir = args["corpus-dir"]
output_dir = args["output-dir"]
model = args["model"]
if !isdir(output_dir)
    mkdir(output_dir)
end

temp = "temp_folder"
# if model is single file, copy it in a temp folder
if !isdir(model)
    mkdir(temp)
    cp(model, temp * "/" * model)
    model = temp
end

for curr_model in readdir(model)

    vm, dict = load_model(model * "/" * curr_model);

    for file in readdir(corpus_dir)
        println("Annotating " * file * " with model " * curr_model)
        corpus_file = corpus_dir * "/" * file
        output_file = output_dir * "/" * file * curr_model * "_disamb"
        annotate_file(corpus_file, output_file, vm, dict, separator, min_prob, win)
    end

end

# if temp folder was created, delete it
if isdir(temp)
    rm(temp, recursive=true)
end

ext: jl

max_stars_repo_name: vsbogd/language-learning

max_issues_repo_name: vsbogd/language-learning

max_forks_repo_name: vsbogd/language-learning


text:
using Gadgetron
using Test

@testset "Gadgetron.jl" begin
    # Write your tests here.
end


text:
using MLDatasets
using HDF5

type Dataset
    traindata
    testdata
    worddict
    chardict
    tagdict
end

type Token
    word::Int
    chars::Vector{Int}
end

function setup_data()
    h5file = joinpath(dirname(@__FILE__), "wordembeds_nyt100.h5")
    isfile(h5file) || download("https://cl.naist.jp/~shindo/wordembeds_nyt100.h5", h5file)
    words = h5read(h5file, "s")
    #wordembeds = Lookup(h5read(h5file,"v"))
    #charembeds = Lookup(Float32,100,10)

    worddict = IntDict(words)
    chardict = IntDict{String}()
    tagdict = IntDict{String}()

    traindata = UD_English.traindata()
    testdata = UD_English.testdata()
    #traindata = CoNLL.read(".data/wsj_00-18.conll")
    #testdata = CoNLL.read(".data/wsj_22-24.conll")
    info("# sentences of train data: $(length(traindata))")
    info("# sentences of test data: $(length(testdata))")

    traindata = encode(traindata, worddict, chardict, tagdict, true)
    testdata = encode(testdata, worddict, chardict, tagdict, false)
    info("# words: $(length(worddict))")
    info("# chars: $(length(chardict))")
    info("# tags: $(length(tagdict))")

    Dataset(traindata, testdata, worddict, chardict, tagdict)
end

function encode(data::Vector, worddict, chardict, tagdict, append::Bool)
    data_x, data_y = Vector{Token}[], Vector{Int}[]
    unkword = worddict["UNKNOWN"]
    for sent in data
        push!(data_x, Token[])
        push!(data_y, Int[])
        for items in sent
            word, tag = items[2], items[5]
            word0 = replace(word, r"[0-9]", '0')
            wordid = get(worddict, lowercase(word0), unkword)

            chars = Vector{Char}(word0)
            if append
                charids = map(c -> push!(chardict,string(c)), chars)
            else
                charids = map(c -> get(chardict,string(c),0), chars)
            end
            tagid = push!(tagdict, tag)
            token = Token(wordid, charids)
            push!(data_x[end], token)
            push!(data_y[end], tagid)
        end
    end
    data_x, data_y
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# This includes a few helper variables and functions that provide information about the
# test environment (command line flags, current module, etc).
# This file can be included multiple times in the same module if necessary,
# which can happen with unisolated test runs.

if !@isdefined(testenv_defined)
    const testenv_defined = true
    if haskey(ENV, "JULIA_TEST_EXEFLAGS")
        const test_exeflags = `$(Base.shell_split(ENV["JULIA_TEST_EXEFLAGS"]))`
    else
        const test_exeflags = Base.julia_cmd()
        filter!(test_exeflags.exec) do c
            return !(startswith(c, "--depwarn") || startswith(c, "--check-bounds"))
        end
        push!(test_exeflags.exec, "--check-bounds=yes")
        push!(test_exeflags.exec, "--startup-file=no")
        push!(test_exeflags.exec, "--depwarn=error")
    end

    if haskey(ENV, "JULIA_TEST_EXENAME")
        popfirst!(test_exeflags.exec)
        const test_exename = `$(Base.shell_split(ENV["JULIA_TEST_EXENAME"]))`
    else
        const test_exename = popfirst!(test_exeflags.exec)
    end

    addprocs_with_testenv(X; kwargs...) = addprocs(X; exename=test_exename, exeflags=test_exeflags, kwargs...)

    const curmod = @__MODULE__
    const curmod_name = fullname(curmod)
    const curmod_str = curmod === Main ? "Main" : join(curmod_name, ".")
    const curmod_prefix = "$(["$m." for m in curmod_name]...)"
end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
@testset "Utils" begin
    @testset "edge encoding/decoding" begin
        # not is_bidirected
        n = 5
        s = [1,1,2,3,3,4,5]
        t = [1,3,1,1,2,5,5]

        # directed=true
        idx, maxid = GNNGraphs.edge_encoding(s, t, n)
        @test maxid == n^2
        @test idx == [1, 3, 6, 11, 12, 20, 25]

        sdec, tdec = GNNGraphs.edge_decoding(idx, n)
        @test sdec == s
        @test tdec == t

        n1, m1 = 10, 30
        g = rand_graph(n1, m1)
        s1, t1 = edge_index(g)
        idx, maxid = GNNGraphs.edge_encoding(s1, t1, n1)
        sdec, tdec = GNNGraphs.edge_decoding(idx, n1)
        @test sdec == s1
        @test tdec == t1

        # directed=false
        idx, maxid = GNNGraphs.edge_encoding(s, t, n, directed=false)
        @test maxid == n*(n+1)÷2
        @test idx == [1, 3, 2, 3, 7, 14, 15]

        mask = s .> t
        snew = copy(s)
        tnew = copy(t)
        snew[mask] .= t[mask]
        tnew[mask] .= s[mask]
        sdec, tdec = GNNGraphs.edge_decoding(idx, n, directed=false)
        @test sdec == snew
        @test tdec == tnew

        n1, m1 = 6, 8
        g = rand_graph(n1, m1)
        s1, t1 = edge_index(g)
        idx, maxid = GNNGraphs.edge_encoding(s1, t1, n1, directed=false)
        sdec, tdec = GNNGraphs.edge_decoding(idx, n1, directed=false)
        mask = s1 .> t1
        snew = copy(s1)
        tnew = copy(t1)
        snew[mask] .= t1[mask]
        tnew[mask] .= s1[mask]
        @test sdec == snew
        @test tdec == tnew
    end
end


text:
module Varpro
using LinearAlgebra
using SparseArrays

include("VarproTypes.jl")

using NL2sol

export varpro, FitContext, NL2SOL, NL2SNO

"""
  # Description

    varpro solves a separable nonlinear least squares problem.

        varpro(ctx::FitContext)
        Returns: [alpha, c, wresid, wresid_norm, y_est, Regression]

        ctx is constructed by the FitContext() constructor.  See types.jl
        for more info

    Given a set of m observations y[1],...,y[m]
    this program computes a weighted least squares fit using the model

        eta(alpha,c,t) =
            c_1 * phi_1 (alpha,t) + ...  + c_n * phi_n (alpha,t)
        (possibly with an extra term  + phi_{n+1} (alpha,t) ).

    This program determines optimal values of the q nonlinear parameters
    alpha and the n linear parameters c, given observations y at m
    different values of the "time" t and given evaluation of phi and
    (optionally) derivatives of phi.

    Varpro calls either NL2sol, NL2sno (fd Jacobian) which solves a
    non-linear least squares problem.

    What distinguishes varpro from levenberg_marquardt is that, for efficiency and
    reliability, varpro causes the nonlinear solver NL2sol to only iterate on the
    nonlinear parameters.  Given the information in phi and dphi, this
    requires an intricate but inexpensive computation of partial
    derivatives, and this is handled by the varpro function formJacobian.

    nl2sol is in the NL2sol.jl package.

    The original Fortran implementation of the variable projection
    algorithm (ref. 2) was modified in 1977 by John Bolstad
    Computer Science Department, Serra House, Stanford University,
    using ideas of Linda Kaufman (ref. 5) to speed up the
    computation of derivatives.  He also allowed weights on
    the observations, and computed the covariance matrix.

    The orginal MATLAB version of the Varpro program is documented in
    Dianne P. O'Leary and Bert W. Rust,
    "Variable Projection for Nonlinear Least Squares Problems",
    Computational Optimization and Applications (2012)
    doi 10.1007/s10589-012-9492-9.

    US National Inst. of Standards and Technology, 2010.

    Main reference:

        0.  Gene H. Golub and V. Pereyra, 'Separable nonlinear least
            squares: the variable projection method and its applications,'
            Inverse Problems 19, R1-R26 (2003).

    See also these papers, cited in John Bolstad's Fortran code:

    1.  Gene H. Golub and V. Pereyra, 'The differentiation of
        pseudo-inverses and nonlinear least squares problems whose
        variables separate,' SIAM J. Numer. Anal. 10, 413-432
        (1973).
    2.  ------, same title, Stanford C.S. Report 72-261, Feb. 1972.
    3.  Michael R. Osborne, 'Some aspects of non-linear least
        squares calculations,' in Lootsma, Ed., 'Numerical Methods
        for Non-Linear Optimization,' Academic Press, London, 1972.
    4.  Fred Krogh, 'Efficient implementation of a variable projection
        algorithm for nonlinear least squares problems,'
        Comm. ACM 17:3, 167-169 (March, 1974).
    5.  Linda Kaufman, 'A variable projection method for solving
        separable nonlinear least squares problems', B.I.T. 15,
        49-57 (1975).
    6.  C. Lawson and R. Hanson, Solving Least Squares Problems,
        Prentice-Hall, Englewood Cliffs, N. J., 1974.

    These books discuss the statistical background:

    7.  David A. Belsley, Edwin Kuh, and Roy E. Welsch, Regression
        Diagnostics, John Wiley & Sons, New York, 1980, Chap. 2.
    8.  G.A.F. Seber and C.J. Wild, Nonlinear Regression,
        John Wiley & Sons, New York, 1989, Sec. 2.1, 5.1, and 5.2.

    Dianne P. O'Leary, NIST and University of Maryland, February 2011.
    Bert W. Rust,      NIST                             February 2011.
    Comments updated 07-2012.

    MODIFIED: May 2016
    Ported to Julia, extended to handle complex residuals of complex
    parameters (both linear and non-linear), restructured to use context
    struct, style changes and added ability to call nl2sol as the nonlinear
    solver

    Don MacMillen

# Input

    On Input and stored in the FitContext struct ctx are:

    y    m x 1   vector containing the m observations (may be complex)
    w    m x 1   vector of weights used in the least squares
        fit.  We minimize the norm of the weighted residual
        vector r, where, for i=1:m,

        r[i] = w[i] * (y[i] - eta(alpha, c, t[i,:])).

        Therefore, w[i] should be set to 1 divided by
        the standard deviation in the measurement y[i].
        If this number is unknown, set w[i] = 1.
        Only real weights are supported

    alpha q x 1  initial estimates of the parameters alpha.
        If alpha = [], Varpro assumes that the problem
        is linear and returns estimates of the c parameters.
        May be complex
    n            number of linear parameters c
    ada          a function handle, described below.
    gada         a function handle, described below.
    lb    q x 1  lower bounds on the parameters alpha.
        (Optional)   (Omit this argument or use [] if there are
        no lower bounds.)  May be complex but currently unimplmented
    ub    q x 1  upper bounds on the parameters alpha.
        (Optional)   (Omit this argument or use [] if there are
        no upper bounds.)  May be complex but currently unimplemented

# Output

    alpha  q      estimates of the nonlinear parameters.
    c      n      estimates of the linear parameters.
    wresid m      weighted residual vector, with i-th component
        w(i) * (y(i) - eta(alpha, c, t(i,:))).
    wresid_norm   norm of wresid.
    y_est  m      the model estimates = eta(alpha, c, t(i,:)))
    regression    a structure containing diagnostics about the model

                **************************************************
                *                C a u t i o n:                  *
                *   The theory that makes statistical            *
                *   diagnostics useful is derived for            *
                *   linear regression, with no upper- or         *
                *   lower-bounds on variables.                   *
                *   The relevance of these quantities to our     *
                *   nonlinear model is determined by how well    *
                *   the linearized model (Taylor series model)   *
                *         eta(alpha_true, c_true)                *
                *            +  phi * (c  - c_true)              *
                *            + dphi * (alpha - alpha_true)       *
                *   fits the data in the neighborhood of the     *
                *   true values for alpha and c, where phi       *
                *   and dphi contain the partial derivatives     *
                *   of the model with respect to the c and       *
                *   alpha parameters, respectively, and are      *
                *   defined in ada.                              *
                **************************************************

    regression.results:
        This structure includes information on the solution
        process, including the number of iterations,
        termination criterion, and exitflag from lsqnonlin.

    regression.rank is the computed rank of the
        matrix for the linear subproblem.  If this equals
        n, then the linear coefficients are well-determined.
        If it is less than n, then although the model might
        fit the data well, other linear coefficients might
        give just as good a fit.

    regression.sigma:
        The estimate of the standard deviation is the
        weighted residual norm divided by the square root
        of the number of degrees of freedom.
        This is also called the "regression standard error"
        or the square-root of the weighted SSR (sum squared
        residual) divided by the square root of the
        number of degrees of freedom.

    regression.rms:
       The "residual mean square" is equal to sigma^2:
       RMS = wresid_norm^2 / (m-n+q)

    regression.coef_determ:
        The "coefficient of determination" for the fit,
        also called the square of the multiple correlation
        coefficient, is sometimes called R^2.
        It is computed as 1 - wresid_norm^2/CTSS,
        where the "corrected total sum of squares"
        CTSS is the norm-squared of W*(y-y_bar),
        and the entries of y_bar are all equal to
        (the sum of W_i^2 y_i) divided by (the sum of W_i^2).
        A value of .95, for example, indicates that 95 per
        cent of the CTTS is accounted for by the fit.

    regression.covmx: (n+q) x (n+q)
        This is the estimated variance/covariance matrix for
        the parameters.  The linear parameters c are ordered
        first, followed by the nonlinear parameters alpha.
        This is empty if dphi is not computed by ada.

    regression.cormx: (n+q) x (n+q)
        This is the estimated correlation matrix for the
        parameters.  The linear parameters c are ordered
        first, followed by the nonlinear parameters alpha.
        This is empty if dphi is not computed by ada.

    regression.std_param: (n+q) x 1
        This vector contains the estimate of the standard
        deviation for each parameter.
        The k-th element is the square root of the k-th main
        diagonal element of Regression.CovMx
        This is empty if dphi is not computed by ada.

    regression.t_ratio:   (n+q)
        The t-ratio for each parameter is equal to the
        parameter estimate divided by its standard deviation.
        (linear parameters c first, followed by alpha)
        This is empty if dphi is not computed by ada.

    regression.standardized_wresid:
        The k-th component of the "standardized weighted
        residual" is the k-th component of the weighted
        residual divided by its standard deviation.
        This is empty if dphi is not computed by ada.

    Specification of the function ada, which computes information
    related to phi:

    function  phi, dphi, ind = ada(alpha)

    This function computes phi and returns it.  Additionally, it
    computes more information that is needed in formJacobian and
    stores it in the context

    On Input:

        alpha q     contains the current value of the alpha parameters.

        Note:  If more input arguments are needed, use the standard
               Matlab syntax to accomplish this.  For example, if
               the input arguments to ada are t, z, and alpha, then
               before calling varpro, initialize t and z, and in calling
               varpro, replace "@ada" by "@(alpha)ada(t,z,alpha)".

     On Output:

        phi   m x n1   where phi(i,j) = phi_j(alpha,t_i).
                       (n1 = n if there is no extra term;
                        n1 = n+1 if an extra term is used)
        dphi  m x p    where the columns contain partial derivative
                       information for phi and p is the number of
                       columns in Ind
                       (or dphi = [] if derivatives are not available).
        Ind   2 x p    Column k of dphi contains the partial
                       derivative of phi_j with respect to alpha_i,
                       evaluated at the current value of alpha,
                       where j = Ind(1,k) and i = Ind(2,k).
                       Columns of dphi that are always zero, independent
                       of alpha, need not be stored.
        Example:  if  phi_1 is a function of alpha_2 and alpha_3,
                  and phi_2 is a function of alpha_1 and alpha_2, then
                  we can set
                          Ind = [ 1 1 2 2
                                  2 3 1 2 ]
                  In this case, the p=4 columns of dphi contain
                          d phi_1 / d alpha_2,
                          d phi_1 / d alpha_3,
                          d phi_2 / d alpha_1,
                          d phi_2 / d alpha_2,
                  evaluated at each t_i.
                  There are no restrictions on how the columns of
                  dphi are ordered, as long as Ind correctly specifies
                  the ordering.

        If derivatives dphi are not available, then set dphi = Ind = [].

"""
function varpro(ctx)
    y, W, q, m, n1, n =  ctx.y, ctx.w, ctx.q, ctx.m, ctx.n1, ctx.n

    T = eltype(ctx.alpha)
    nl2_msg = ""

    #
    # Solve the least squares problem using NL2sol or, if there
    # are no nonlinear parameters, using the SVD procedure in formJacobian.
    #
    regression = Regression(T, q, n, m)

    if q > 0  # The problem is nonlinear.

        if ctx.iscomplex
            alpha_real = [real(ctx.alpha); imag(ctx.alpha)]
            mreal = 2 * ctx.m
        else
            alpha_real = ctx.alpha
            mreal = ctx.m
        end

        f(a, r) = f_lsq(a, r, ctx)
        g(a, j) = g_lsq(a, j, ctx)
        iv, v = nl2_set_defaults(mreal, length(alpha_real))
        iv[MXFCAL] = ctx.mxfcal
        iv[MXITER] = ctx.mxiter
        iv[PRUNIT] = 0
        if ctx.opto == NL2SOL
            results = nl2sol(f, g, alpha_real, mreal, iv, v)
            ctx.verbose && println("\nNL2sol return code: ", return_code[iv[1]])
        else
            results = nl2sno(f, alpha_real, mreal, iv, v)
            ctx.verbose && println("\nNL2sno return code: ", return_code[iv[1]])
        end

        ctx.verbose && println(nl2_msg)
        ctx.verbose && @show(results)
        alpha_real .= results.minimum
        if ctx.iscomplex
            ctx.alpha .= complex.(alpha_real[1:ctx.q], alpha_real[ctx.q+1:end])  ## cwt barks here
        end
        wresid_norm2 = results.f_minimum
        f_lsq(alpha_real, ctx.wresid_real, ctx)
        g_lsq(alpha_real, ctx.jac_real,  ctx)
        r, Jacobian, phi, dphi, y_est = ctx.wresid, ctx.jac, ctx.phi, ctx.dphi, ctx.y_est

        wresid = r
        wresid_norm = sqrt(wresid_norm2)
        regression.results = results
        regression.rank = ctx.rank

    else       # The problem is linear.

        ctx.verbose && println("VARPRO problem is linear, since length(alpha) = 0")

        Jacobian = formJacobian(ctx)
        c = ctx.c
        wresid = ctx.wresid
        y_est = ctx.y_est
        regression.report.rank = ctx.rank
        wresid_norm = norm(wresid)
        wresid_norm2 = wresid_norm * wresid_norm
    end

    # Calculate sample variance,  the norm-squared of the residual
    #    divided by the number of degrees of freedom.
    sigma2 = wresid_norm2 / float(m - n - q)

    # Compute  Regression.sigma:
    #     square-root of weighted residual norm squared divided
    #     by number of degrees of freedom.
    regression.sigma = sqrt(sigma2)

    # Compute Regression.coef_determ:
    #     The coeficient of determination for the fit, also called the square
    #     of the multiple correlation coefficient, or R^2. It is computed as
    #     1 - wresid_norm^2/CTSS, where the corrected total sum of squares
    #     CTSS is the norm-squared of W*(y-y_bar), and the entries of y_bar
    #     are all equal to (the sum of W_i y_i) divided by (the sum of W_i).

    w = diag(W)
    y_bar = sum(w.*y) / sum(w)
    CTTS = norm(W * (y .- y_bar)) ^ 2
    regression.coef_determ = 1 - wresid_norm^2 / CTTS

    # Compute  regression.RMS = sigma^2:
    #     the weighted residual norm divided by the number of degrees of
    #     freedom. RMS = wresid_norm / sqrt(m-n+q)

    regression.rms = sigma2

    if !isempty(dphi) && !ctx.skip_stats
        # Calculate the covariance matrix CovMx, which is sigma^2 times the
        # inverse of H'*H, where
        #              H = W*[phi,J]
        # contains the partial derivatives of  wresid  with
        # respect to the parameters in alpha and c.

        if ctx.iscomplex
            dphi = dphi[:,1:q]
        end
        xx, pp = size(dphi)
        J = zeros(eltype(ctx.alpha), m, q)
        for kk = 1:pp
            j = ctx.ind[1, kk]
            i = ctx.ind[2, kk]
            if j > n
                J[:, i] = J[:, i] + dphi[:, kk]
            else
                J[:, i] = J[:, i] + ctx.c[j] * dphi[:, kk]
            end
        end
        # Uses compact pivoted QR.
        F = qr(W*[phi[:, 1:n] J], Val(true))
        Qj = F.Q
        Rj = F.R
        Pj = F.p
        T2 = Rj \ diagm(ones((size(Rj, 1))))
        covmx = sigma2 * T2 * T2'
        regression.covmx[Pj, Pj] = covmx  # Undo the pivoting permutation.

        # Compute  regression.CorMx:
        #     estimated correlation matrix (n+q) x (n+q) for the
        #     parameters.  The linear parameters are ordered
        #     first, followed by the nonlinear parameters.
        d = 1 ./ sqrt.(diag(regression.covmx))
        D = spdiagm(0 => d)
        regression.cormx = D * regression.covmx * D

        # Compute  regression.std_param:
        #     The k-th element is the square root of the k-th main
        #     diagonal element of CovMx.
        regression.std_param = sqrt.(diag(regression.covmx))

        # Compute  regression.t_ratio:
        #     parameter estimates divided by their standard deviations.
        regression.t_ratio = [ctx.c; ctx.alpha] .* d

        # Compute  regression.standardized_wresid:
        #     The k-th component is the k-th component of the
        #     weighted residual, divided by its standard deviation.
        #     Let X = W*[phi, J],
        #        h[k] = k-th main diagonal element of covariance
        #               matrix for wresid
        #             = k-th main diagonal element of X*inv(X'*X)*X'
        #             = k-th main diagonal element of Qj*Qj'
        #     Then the standard deviation is estimated by sigma*sqrt(1-h[k])

        temp = zeros(eltype(ctx.alpha), m)
        for k = 1:m
            temp[k] = (Qj[k, :] * Qj[k, :]')[1]
        end
        regression.std_wresid = wresid ./(regression.sigma*sqrt.(1 .- temp))
     end

    if ctx.verbose
        println(" ")
        println("VARPRO Results:")
        println(" Linear Parameters: ", ctx.c)

        println(" Nonlinear Parameters: ", ctx.alpha)

        println(" Norm-squared of weighted residual  = ", wresid_norm2)
        println(" Norm-squared of data vector        = ", norm(w.*y)^2)
        println(" Norm         of weighted residual  = ", wresid_norm)
        println(" Norm         of data vector        = ", norm(w.*y))
        println(" Expected error of observations     = ", regression.sigma)
        println(" Coefficient of determination       = ", regression.coef_determ)
    end

    return ctx.alpha, ctx.c, wresid, wresid_norm, y_est, regression
end

function update_alpha!(alpha_trial, ctx)
    if ctx.iscomplex
        # We do it this way since alpha_trial can be a NL2Array and we haven't
        # implemented slicing for that type.
        for i = 1:ctx.q
            ctx.alpha_real[i] = alpha_trial[i]
            ctx.alpha_real[i+ctx.q] = alpha_trial[i+ctx.q]
        end
        ctx.alpha = complex.(ctx.alpha_real[1:ctx.q], ctx.alpha_real[ctx.q+1:end])
    else
        # again, cannot use a colon assignment here
        for i = 1:ctx.q
            ctx.alpha[i] = alpha_trial[i]
        end
    end
    return nothing
end

"""
# Description
    This function is used by a nonlinear least squares solver NL2sol
    to compute wresid, the current estimate of the weighted residual
    (but it is wrapped in a closure first so it is only a function
     of alpha_trial and wresid).

# Input
    alpha_trial  trial vector of non-linear parameters
    ctx: a varpro context
    function f_lsq(alpha_trial, ctx)

# Output
    wresid.  Also calculates current estimates of c, phi, y_est,
        and myrank and stores in the context struct

# NB
    It uses the user-supplied function ada

    An added complication is that both alpha_trial and r must be
    real for the optimizer.  However both of these values may be
    complex.  In that case, we check ctx to see and then recombine
    into the complex values for ada and then blast back to real
    before returning
"""
function f_lsq(alpha_trial, r, ctx)

    # Convert back to complex if necessary
    update_alpha!(alpha_trial, ctx)
    W = ctx.w

    # We determine the optimal linear parameters c for
    # the given values of alpha, the resulting residual and the y estimate
    #
    # We use the singular value decomposition to solve the linear least
    # squares problem
    #
    #    min_{c} || W resid ||.
    #       resid =  y - phi * c.
    #
    # If W*phi has any singular value less than m * eps() its largest singular value,
    # these singular values are set to zero.
    # The following values on stored on the FitContext:
    #      c        n x 1 the optimal linear parameters for this choice of alpha.
    #      wresid   m x 1 the weighted residual = W(y - phi * c)
    #      y_est    m x 1 the model estimates = phi * c)
    #      rank     1 x 1 the rank of the matrix W*phi.

    ctx.phi = ctx.ada(ctx.alpha, ctx)
    U, S, V = svd(W * ctx.phi[:, 1:ctx.n], full=true)
    ctx.U = U
    ctx.s = S
    ctx.V = V

    if ctx.n > 0
        s = S
    else    #  n = 0

        ctx.c = ctx.iscomplex ? Base.Complex{Float64}[] : Float64[]
        ctx.y_est .= ctx.phi
        ctx.wresid .= W * (y - ctx.y_est)
        r .= ctx.wresid
        ctx.rank = 1
        return ctx.wresid
    end

    # Time to "regularize"
    tol = ctx.m * eps()
    ctx.rank = sum(s .> tol*s[1] ) # number of singular values > tol*norm(W*phi)
    s = s[1:ctx.rank]

    if ctx.rank < ctx.n && ctx.verbose
        println("Warning from VARPRO:")
        println("   The linear parameters are currently not well-determined.")
        println("   The rank of the matrix in the subproblem is $(ctx.rank)")
        println("   which is less than the $(ctx.n) linear parameters.")
    end

    yuse = copy(ctx.y)
    if ctx.n < ctx.n1
        yuse  = ctx.y - ctx.phi[:, n1]  # extra function phi(:,n+1)
    end
    temp  = U[:, 1:ctx.rank]' * (W * yuse)
    ctx.c = V[:, 1:ctx.rank] * (temp ./ s)
    ctx.y_est = ctx.phi[:, 1:ctx.n] * ctx.c
    ctx.wresid = W * (yuse - ctx.y_est)

    # This is most certainly wrong now, but we don't use the case n < n1
    if ctx.n < ctx.n1
        ctx.y_est += ctx.phi[:, ctx.n1]
    end
    if ctx.iscomplex
        ctx.wresid_real[1:ctx.m] = real(ctx.wresid)
        ctx.wresid_real[ctx.m+1:end] = imag(ctx.wresid)
    else
        ctx.wresid_real .= ctx.wresid
    end

    # Set r[:] for NL2sol. This is unrolled because r might be an NL2Array
    for i in 1:length(ctx.wresid_real)
        r[i] = ctx.wresid_real[i]
    end

    return ctx.wresid_real
end

function g_lsq(alpha_trial, jnl2, ctx)
    # function g_lsq(alpha_trial, ctx)
    #
    # Returns: jacobian, but also calculates dphi and stores it in ctx
    #
    # This is an 'extra' function that does not appear in the MATLAB version
    # of Varpro.  It is needed by levenberg_marquardt solver to compute the
    # the Jacobian matrix for the nonlinear parameters
    #
    # It also computes dphi the partial derivatives of phi (if available).
    # It uses the user-supplied function gada and the Varpro function formJacobian.
    #
    # We depend on f_lsq to update ctx.alpha and ctx.alpha_real.  In effect,
    # we assume that we never ask for a gradient unless we have first asked for
    # for the residual.
    f_lsq(alpha_trial, ctx.wresid_real, ctx)

    ctx.gada(ctx.alpha, ctx)
    formJacobian(ctx)

    if ctx.iscomplex
        # This is true because of Cauchy-Riemann Equations
        ctx.jac_real = [real(ctx.jac)  -imag(ctx.jac);
                        imag(ctx.jac)   real(ctx.jac)]
    else
        ctx.jac_real = ctx.jac
    end

    # Now must set the memory that nl2sol uses.
    m, n = size(ctx.jac_real)
    for j = 1:n
        for i in 1:m
            jnl2[i, j] = ctx.jac_real[i, j]
        end
    end

    return ctx.jac_real
end

"""
    function formJacobian(ctx)

    Returns: Jacobian

    This function computes the Jacobian and dphi
    It is used by the functions Varpro and f_lsq.

    Notation: there are m data observations
                        n1 basis functions (columns of phi)
                        n linear parameters c
                            (n = n1, or n = n1 - 1)
                        q nonlinear parameters alpha
                        p nonzero columns of partial derivatives of phi

    Input: (on ctx)
        alpha  q x 1   the nonlinear parameters,
        phi    m x n1  the basis functions phi(alpha),
        dphi   m x p   the partial derivatives of phi

    The variables W, y, q, m, n1, and n are also used.

    Output:

        Jacobian  m x p the Jacobian matrix, with J[i, k] = partial
                    derivative of W resid[i] with respect to alpha[k].
"""
function formJacobian(ctx)
    if ctx.q == 0 || isempty(ctx.dphi)
        ctx.jac = []
        return
    end

    # Compute the Jacobian.
    # There are two pieces, which we call Jac1 and Jac2,
    # with Jacobian = - (Jac1 + Jac2).
    #
    # The formula for Jac1 is (P D(W*phi) pinv(W*phi) y,
    #             and Jac2 is ((W*phi)^+})^T (P D(W*phi))^T y.
    # where  P           is the projection onto the orthogonal complement
    #                       of the range of W*phi,
    #        D(W*phi)    is the m x n x q tensor of derivatives of W*phi,
    #        pinv(W*phi) is the pseudo-inverse of W*phi.
    #  (See Golub&Pereyra (1973) equation (5.4).  We use their notational
    #   conventions for multiplications by tensors.)
    #
    # Golub & Pereyra (2003), p. R5 break these formulas down by columns:
    #     The j-th column of Jac1 is P D_j pinv(W*phi) y
    #                             =  P D_j c
    # and the j-th column of Jac2 is (P D_j pinv(W*phi))^T y,
    #                             =  (pinv(W*phi))^T D_j^T P^T y
    #                             =  (pinv(W*phi))^T D_j^T wresid.
    # where D_j is the m x n matrix of partial derivatives of W*phi
    #     with respect to alpha(j).
    #
    # We begin the computation by precomputing
    #       Wdphi, which contains the derivatives of W*phi, and
    #       Wdphi_r, which contains Wdphi' * wresid.

    Wdphi = ctx.w * ctx.dphi
    Wdphi_r = Wdphi' * ctx.wresid
    T2 = zeros(eltype(ctx.alpha), ctx.n1, ctx.q)
    ctemp = ctx.c
    if ctx.n1 > ctx.n
        ctemp = [ctemp; 1]
    end

    #   Now we work column-by-column, for j=1:q.
    #
    #   We form Jac1 = D(W*phi) ctemp.
    #   After the loop, this matrix is multiplied by
    #        P = U[:, ctx.rank+1:m] * (U[:, ctx.rank+1:m]' * Jac1)
    #   to complete the computation.
    #
    #   We also form  T2 = (D_j(W*phi))^T wresid  by unpacking
    #   the information in Wdphi_r, using ind.
    #   After the loop, T2 is multiplied by the pseudoinverse
    #       (pinv(W*phi))^T = U[:,1:ctx.rank] * diagm(1./s) * (V[:, 1:ctx.rank]'
    #   to complete the computation of Jac2.
    #   Note: if n1 > n, last row of T2 is not needed.
    for j = 1:ctx.q                            # for each nonlinear parameter alpha[j]
        range = findall(ctx.ind[2, :] .== j)   # columns of Wdphi relevant to alpha[j]
        indrows = vec(ctx.ind[1, range])       # relevant rows of ctemp, need as vector
        ctx.jac1[:, j] = Wdphi[:, range] * ctemp[indrows]
        T2[indrows, j] = Wdphi_r[range]
    end

    ctx.jac1 .= ctx.U[:, ctx.rank+1:ctx.m] * (ctx.U[:, ctx.rank+1:ctx.m]' * ctx.jac1)
    T2 = diagm(0 => 1 ./ ctx.s[1:ctx.rank]) * (ctx.V[:, 1:ctx.rank]' * T2[1:ctx.n, :])
    ctx.jac2 .= ctx.U[:, 1:ctx.rank] * T2

    ctx.jac .= -(ctx.jac1 + ctx.jac2)

    if ctx.debug
        if ctx.neglect
            ctx.verbose && println("VARPRO norm(neglected Jacobian)/norm(Jacobian) = ",
                                   norm(ctx.jac2,"fro")/norm(ctx.jac,"fro"))
            ctx.verbose && println("neglecting Jac2")
            ctx.jac .= -ctx.jac1
        end
    end
    return ctx.jac
end

end # module Varpro


text:
"""
`get_lambda(fmhz)`

Returns wavelength in meters at frequency `fmhz`.
Argument `fmhz` is in mega-hertz.
"""
function get_lambda(fmhz::Real)
    c = 299792458.0       # m/s
    return c / (fmhz * 1e6)
end

# if delta is the length of a segment, and lambda is the wavelength
# note that lambda / delta is segments per wavelength
# we want delta < 0.05 * lambda (according to NEC-2 guide)
# Or, we want lambda / delta > 20 => at least 20 segments per wavelength
# we'll call it 25 to be safe
#
# length is L meters
# wavelength is lambda meters/wave
# then L / lambda is how many wavelengths long this is
# So...
#  25 * L / lambda is minimum number of segments.
# Obviously, round this up; then make sure it's odd (easier to split)
#
# other things I want to do:
#  add a minimum number of segments but ensure they aren't too small
function divide(L::Real, fmhz::Real)
    lambda = get_lambda(fmhz)
    n_segments = round(Int, 25.0 * L / lambda, RoundUp)

    # make sure this is an odd number
    if iseven(n_segments)
        n_segments += 1
    end

    # TODO:
    # make sure we have some minimum number of segments but not too many

    return n_segments
end

"""
`awg(n::Int)`

Returns radius in meters of No. `n` AWG wire.
"""
awg(n::Int) = 0.5 * 0.001 * 0.127 * 92.0 ^ ((36.0-n) / 39.0)

""" Radius in meters of size 1 AWG (American Wire Gauge) """
const AWG1 = awg(1)
""" Radius in meters of size 2 AWG (American Wire Gauge) """
const AWG2 = awg(2)
""" Radius in meters of size 3 AWG (American Wire Gauge) """
const AWG3 = awg(3)
""" Radius in meters of size 4 AWG (American Wire Gauge) """
const AWG4 = awg(4)
""" Radius in meters of size 5 AWG (American Wire Gauge) """
const AWG5 = awg(5)
""" Radius in meters of size 6 AWG (American Wire Gauge) """
const AWG6 = awg(6)
""" Radius in meters of size 7 AWG (American Wire Gauge) """
const AWG7 = awg(7)
""" Radius in meters of size 8 AWG (American Wire Gauge) """
const AWG8 = awg(8)
""" Radius in meters of size 9 AWG (American Wire Gauge) """
const AWG9 = awg(9)
""" Radius in meters of size 10 AWG (American Wire Gauge) """
const AWG10 = awg(10)
""" Radius in meters of size 11 AWG (American Wire Gauge) """
const AWG11 = awg(11)
""" Radius in meters of size 12 AWG (American Wire Gauge) """
const AWG12 = awg(12)
""" Radius in meters of size 13 AWG (American Wire Gauge) """
const AWG13 = awg(13)
""" Radius in meters of size 14 AWG (American Wire Gauge) """
const AWG14 = awg(14)
""" Radius in meters of size 15 AWG (American Wire Gauge) """
const AWG15 = awg(15)
""" Radius in meters of size 16 AWG (American Wire Gauge) """
const AWG16 = awg(16)
""" Radius in meters of size 17 AWG (American Wire Gauge) """
const AWG17 = awg(17)
""" Radius in meters of size 18 AWG (American Wire Gauge) """
const AWG18 = awg(18)
""" Radius in meters of size 19 AWG (American Wire Gauge) """
const AWG19 = awg(19)
""" Radius in meters of size 20 AWG (American Wire Gauge) """
const AWG20 = awg(20)


text:
function exampleHBF(n=1)
    if n == 1
        f1, q1, f2, q2 = HBFex1()
    elseif n == 2
        f1, q1, f2, q2 = HBFex2()
    elseif n == 3
        f1, q1, f2, q2 = Saramaki88()
    elseif n == 4
        f1, q1, f2, q2 = Saramaki90()
    else
        error("No code for example $n available.")
    end
    F1 = convertForm(f1, q1)
    F2 = convertForm(f2, q2)
    return F1, F2
end

function convertForm(f, q)
    n = length(f)
    val = f
    csd = Array{Matrix{Int}}(undef, n)
    for i = 1:n
        csdrep = q[(2*i-1):(2*i),:]
        csd[i] = csdrep[:, csdrep[1,:] .!= 0]
    end
    return val, csd
end

function HBFex1()
    q1 = [0 -4 -5; 1 1 1; 0 -4 -5; -1 -1 -1; -1 -3 -5; 1 1 1; -3 -5 0; -1 -1 0]
    f1 = bunquantize(q1)
    q2 = [-1 -3 -7; 1 1 1; -2 -5 -7; -1 1 1; -3 -8 0; 1 -1 0; -4 -6 -8; -1 -1 -1;
          -4 -8 0; 1 -1 0; -5 -7 -8; -1 -1 -1; -5 0 0; 1 0 0; -6 -7 0; -1 -1 0;
          -7 -8 0; -1 -1 0; -6 0 0; 1 0 0]
    f2 = bunquantize(q2)
    return f1, q1, f2, q2
end

function HBFex2()
    q1 = [0 -4 -7; 1 -1 1; -1 -3 -6; -1 -1 -1; -2 -4 -7; 1 -1 1]
    f1 = bunquantize(q1)
    q2 = [-1 -3 -8; 1 1 -1; -2 -4 -9; -1 1 -1; -3 -5 -9; 1 -1 1; -4 -7 -8;
          -1 1 1; -5 -8 -11; 1 -1 -1; -6 -9 -11; -1 1 -1]
    f2 = bunquantize(q2)
    return f1, q1, f2, q2
end

function Saramaki88()
    # Coefficients from "Efficient VLSI-Realizable Decimators for Sigma-Delta
    # Analog-to-Digital Converters," T. Saramaki and H. Tenhunen, ISCAS 1988,
    # pp 1525-2528.
    q1 = [0 -2; 1 -1; -2 -9; -1 1]
    f1 = bunquantize(q1)
    q2 = [-1 -3 0; 1 1 0; -3 -5 -6; -1 -1 -1; -4 -6 0; 1 1 0]
    f2 = bunquantize(q2)
    return f1, q1, f2, q2
end

function Saramaki90()
    # Coefficients from "Multiplier-Free Decimation Algorithms for Superresolution
    # Oversampled Converters," T. Saramaki, T. Karema, T. Ritoniemi, and H. Tenhunen
    # ISCAS 1990, vol. 4. pp 1525-2528.
    q1 = [0 -4 -5; 1 1 1; 0 -4 -5; -1 -1 -1; -1 -3 -5; 1 1 1; -3 -5 -20; -1 -1  1]
    f1 = bunquantize(q1)
    q2 = [-1 -3 -7; 1 1 1; -2 -5 -7; -1 1 1; -3 -8 0; 1 -1 0; -4 -6 -8; -1 -1 -1;
          -4 -8 0; 1 -1 0; -5 -7 -8; -1 -1 -1; -5 0 0; 1 0 0; -6 -7 0; -1 -1 0;
          -6 0 0; 1 0 0; -7 -8 0; -1 -1 0; -6 0 0; 1 0 0]
    f2 = bunquantize(q2)
    return f1, q1, f2, q2
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4"]);
gb = SignatureGB.f5([x1, x2, x3, x4]);F = [841*x1^4-5510*x1^3*x2+9025*x1^2*x2^2-638*x1^2*x2*x3+2843*x1^2*x3^2+2726*x1^2*x3*x4+2092*x1*x2^2*x3-9310*x1*x2*x3^2-8820*x1*x2*x3*x4-56*x1*x3*x4^2+8465*x2^4-1104*x2^3*x3-13506*x2^3*x4-13091*x2^2*x3^2-5920*x2^2*x3*x4+8445*x2^2*x4^2-214*x2*x3^3+10066*x2*x3^2*x4+5476*x2*x3*x4^2-3080*x2*x4^3+7585*x3^4+9934*x3^3*x4+3578*x3^2*x4^2+784*x4^4-13800*x1*x2^2+1090*x1*x2*x3+11100*x1*x2*x4+10800*x1*x3^2+5550*x1*x3*x4+190*x2^3+10450*x2^2*x4-5320*x2*x4^2+6089*x1^2-1520*x1*x2+154*x1*x3+9179*x2^2-176*x2*x3+8470*x2*x4+784*x3^2+752*x3*x4-4312*x4^2+14630*x2+5993, -5510*x1^3+18050*x1^2*x2-638*x1^2*x3+4184*x1*x2*x3-9310*x1*x3^2-8820*x1*x3*x4+33860*x2^3-3312*x2^2*x3-40518*x2^2*x4-26182*x2*x3^2-11840*x2*x3*x4+16890*x2*x4^2-214*x3^3+10066*x3^2*x4+5476*x3*x4^2-3080*x4^3-27600*x1*x2+1090*x1*x3+11100*x1*x4+570*x2^2+20900*x2*x4-5320*x4^2-1520*x1+18358*x2-176*x3+8470*x4+14630, -638*x1^2*x2+5686*x1^2*x3+2726*x1^2*x4+2092*x1*x2^2-18620*x1*x2*x3-8820*x1*x2*x4-56*x1*x4^2-1104*x2^3-26182*x2^2*x3-5920*x2^2*x4-642*x2*x3^2+20132*x2*x3*x4+5476*x2*x4^2+30340*x3^3+29802*x3^2*x4+7156*x3*x4^2+1090*x1*x2+21600*x1*x3+5550*x1*x4+154*x1-176*x2+1568*x3+752*x4, 2726*x1^2*x3-8820*x1*x2*x3-112*x1*x3*x4-13506*x2^3-5920*x2^2*x3+16890*x2^2*x4+10066*x2*x3^2+10952*x2*x3*x4-9240*x2*x4^2+9934*x3^3+7156*x3^2*x4+3136*x4^3+11100*x1*x2+5550*x1*x3+10450*x2^2-10640*x2*x4+8470*x2+752*x3-8624*x4];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
using DiffEqSensitivity, ForwardDiff, Distributions, OrdinaryDiffEq,
      LinearAlgebra, Test

function fiip(du,u,p,t)
    du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
    du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
function g(sol)
    J = extract_local_sensitivities(sol,true)[2]
    det(J'*J)
end

u0 = [1.0,1.0]
p = [1.5,1.0,3.0,1.0]
prob = ODEForwardSensitivityProblem(fiip,u0,(0.0,10.0),p,saveat=0:10)
sol = solve(prob, Tsit5())
u0_dist = [Uniform(0.9,1.1), 1.0]
p_dist = [1.5, truncated(Normal(1.5,.1),1.1, 1.9),3.0,1.0]
u0_dist_extended = vcat(u0_dist,zeros(length(p)*length(u0)))

function fiip_expe_SciML_forw_sen_SciML()
    prob = ODEForwardSensitivityProblem(fiip,u0,(0.0,10.0),p,saveat=0:10)
    prob_func = function (prob, i, repeat)
        _prob = remake(prob, u0=[isa(ui,Distribution) ? rand(ui) : ui for ui in u0_dist], p=[isa(pj,Distribution) ? rand(pj) : pj for pj in p_dist])
        _prob
    end
    output_func = function (sol, i)
        (g(sol), false)
    end
    monte_prob = EnsembleProblem(prob;output_func=output_func,prob_func=prob_func)
    sol = solve(monte_prob,Tsit5(),EnsembleSerial(),trajectories=100_000)
    mean(sol.u)
end

@test fiip_expe_SciML_forw_sen_SciML() ≈ 3.56e6 rtol=4e-2


text:
# =========================
# Conversion
# =========================

@inline nodim_msg() = throw(ArgumentError("to use the `static` option you should pass " *
                                          "the system's dimension argument `dim=...`"))

# no-op
_reconvert(Ω0::Zonotope{N, Vector{N}, Matrix{N}}, static::Val{false}, dim, ngens) where {N} = Ω0
_reconvert(Ω0::Zonotope{N, <:SVector, <:SMatrix}, static::Val{true}, dim) where {N} = Ω0

# convert any zonotope to be represented wih regular arrays
_reconvert(Ω0::Zonotope, static::Val{false}, dim, ngens) = Zonotope(Vector(Ω0.center), Matrix(Ω0.generators))

# convert any zonotope to be represented with static arrays
function _reconvert(Ω0::Zonotope{N, VN, MN}, static::Val{true}, dim::Missing, ngens::Missing) where {N, VN, MN}
    n, m = size(Ω0.generators)
    _reconvert(Ω0, static, Val(n), Val(m))
end

function _reconvert(Ω0::Zonotope{N, VN, MN}, static::Val{true}, dim::Val{n}, ngens::Val{p}) where {N, VN, MN, n, p}
    G = Ω0.generators
    m = size(G, 2)
    c = SVector{n, N}(Ω0.center)

    if m == p
        return Zonotope(c, SMatrix{n, p}(G))

    elseif m < p
        # extend with zeros
        Gext = hcat(SMatrix{n, m}(G), zeros(MMatrix{n, p-m, N}))
        return Zonotope(c, Gext)

    else
        throw(ArgumentError("can't reconvert a zontope with $m generators to a " *
                            "zonotope with $p generators; you should reduce it first"))
    end
end

# no-op
_reconvert(Ω0::Hyperrectangle{N, Vector{N}, Vector{N}}, static::Val{false}, dim::Missing) where {N} = Ω0
_reconvert(Ω0::Hyperrectangle{N, Vector{N}, Vector{N}}, static::Val{true}, dim::Missing) where {N} = nodim_msg()
_reconvert(Ω0::Hyperrectangle{N, <:SVector, <:SVector}, static::Val{true}, dim) where {N} = Ω0

# convert any Hyperrectangle to be represented wih regular arrays
function _reconvert(Ω0::Hyperrectangle, static::Val{false}, dim::Missing)
    Ω0 = Hyperrectangle(Vector(Ω0.center), Matrix(Ω0.radius), check_bounds=false)
end

# convert any Hyperrectangle to be represented with static arrays
function _reconvert(Ω0::Hyperrectangle{N, VNC, VNR}, static::Val{true}, dim::Val{n}) where {N, VNC, VNR, n}
    #n = length(Ω0.center) # dimension
    Ω0 = Hyperrectangle(SVector{n, N}(Ω0.center), SVector{n, N}(Ω0.radius), check_bounds=false)
end

# no-op
_reconvert(Φ::Matrix{N}, static::Val{false}, dim) where {N} = Φ
_reconvert(Φ::IntervalMatrix{N}, static::Val{false}, dim) where {N} = Φ
_reconvert(Φ::AbstractMatrix, static::Val{false}, dim) = Matrix(Φ)
_reconvert(Φ::SMatrix, static::Val{true}, dim) = Φ
_reconvert(Φ::AbstractMatrix, static::Val{true}, dim::Missing) = _reconvert(Φ, static, Val(size(Φ, 1)))

function _reconvert(Φ::AbstractMatrix{N}, static::Val{true}, dim::Val{n}) where {N, n}
    #n = size(Φ, 1)
    Φ = SMatrix{n, n, N, n*n}(Φ)
end

function _reconvert(Φ::IntervalMatrix{N, IN, Matrix{IN}}, static::Val{true}, dim::Val{n}) where {N, IN, n}
    #n = size(Φ, 1)
    Φ = IntervalMatrix(SMatrix{n, n, IN, n*n}(Φ))
end

# AbstractVPolygon is not yet available in LazySets
const VPOLY{N, VN} = Union{VPolygon{N, VN}, VPolytope{N, VN}}

# no-op
_reconvert(V::VPOLY{N, VN}, static::Val{false}, dim) where {N, VN<:AbstractVector{N}} = V
_reconvert(V::VPOLY{N, VN}, static::Val{true}, dim) where {N, VN<:SVector{N}} = V

function _reconvert(V::VPOLY{N, VN}, static::Val{true}, dim::Val{n}) where {N, VN<:AbstractVector{N}, n}
    VP = n == 2 ? VPolygon : VPolytope
    return VP([SVector{n, N}(vi) for vi in vertices_list(V)])
end

# dimension is missing
function _reconvert(V::VPOLY{N, VN}, static::Val{true}, dim::Missing) where {N, VN<:AbstractVector{N}}
    _reconvert(V, static, Val(LazySets.dim(V)))
end

# fallback implementation for conversion (if applicable) or overapproximation
function _convert_or_overapproximate(T::Type{<:AbstractPolytope}, X::LazySet)
    if applicable(convert, T, X)
        return convert(T, X)
    elseif applicable(overapproximate, X, T)
        return overapproximate(X, T)
    else
        return convert(T, overapproximate(X, Hyperrectangle))
    end
end

function _convert_or_overapproximate(X::LazySet, T::Type{<:AbstractPolytope})
    return _convert_or_overapproximate(T, X)
end

function _overapproximate(X::Hyperrectangle, T::Type{HPolytope{N, VT}}) where {N, VT}
    # TODO create overapproximation using VT directly
    Y = overapproximate(X, BoxDirections(dim(X)))
    return convert(T, Y)
end

Base.convert(::Type{HPolytope{N, VT}}, P::HPolytope{N, VT}) where {N, VT} = P
Base.convert(::Type{HPolytope{N, VT}}, P) where {N, VT} = HPolytope([HalfSpace(VT(c.a), c.b) for c in constraints_list(P)])

function Base.convert(::Type{Hyperrectangle{N, Vector{N}, Vector{N}}},
                      H::Hyperrectangle{N, SVector{L, N}, SVector{L, N}}) where {N, L}
    return Hyperrectangle(Vector(H.center), Vector(H.radius))
end

function Base.convert(::Type{Singleton},
                      cp::CartesianProduct{N, S1, S2}) where{N, S1<:Singleton{N}, S2<:Singleton{N}}
    x = element(cp.X)
    y = element(cp.Y)
    return Singleton(vcat(x, y))
end

# overload sampling functions
LazySets._default_sampler(X::IntervalArithmetic.Interval) = LazySets._default_sampler(convert(Interval, X))
LazySets._default_sampler(X::IntervalArithmetic.IntervalBox) = LazySets._default_sampler(convert(Hyperrectangle, X))
LazySets.sample(X::IntervalArithmetic.Interval, d::Integer; kwargs...) = sample(convert(Interval, X), d; kwargs...)
LazySets.sample(X::IntervalArithmetic.IntervalBox, d::Integer; kwargs...) = sample(convert(Hyperrectangle, X), d; kwargs...)

# ------------------------------------------------
# Functions to handle splitting of IntervalBoxes
# TODO refactor to LazySets
# See also: LazySets#2651, IntervalArithmetic#444
# ------------------------------------------------

function Base.convert(HT::Type{Hyperrectangle{N, Vector{N}, Vector{N}}}, H::AbstractHyperrectangle) where {N}
    c = convert(Vector{N}, LazySets.center(H))
    r = convert(Vector{N}, radius_hyperrectangle(H))
    return Hyperrectangle(c, r)
end

function Base.convert(HT::Type{Hyperrectangle{N, Vector{N}, Vector{N}}}, B::IntervalBox{D, N}) where {D, N}
    H = convert(Hyperrectangle, B)
    return convert(Hyperrectangle{N, Vector{N}, Vector{N}}, H)
end

function LazySets.split(B::IntervalBox{D, N}, partition::AbstractVector{Int}) where {D, N}
    H = convert(Hyperrectangle{N, Vector{N}, Vector{N}}, B)
    return split(H, partition)
end

# =========================
# Projection
# =========================

# extend LazySets concrete projection for other arg fomats
LazySets.project(X::LazySet, vars::NTuple{D, <:Integer}) where {D} = project(X, collect(vars))
LazySets.project(X::LazySet; vars) = project(X, vars)

# extend LazySets lazy projection for other arg fomats
LazySets.Projection(X::LazySet, vars::NTuple{D, <:Integer}) where {D} = Projection(X, collect(vars))
LazySets.Projection(X::LazySet; vars) = Projection(X, vars)

# ===============================
# Decompositions and partitions
# ===============================

#const Partition{}

# concrete decomposition using a uniform block partition
#using LazySets.Arrays: projection_matrix

#const Partition{N, VT} = AbstractVector{VT} where {VT<:AbstractVector{Int}}

function _decompose(X::LazySet{N},
                    blocks, # ::AbstractVector{<:AbstractVector{Int}}
                    set_type::Type{ST}) where {N, ST<:LazySet}
    n = dim(X)
    result = Vector{ST}(undef, length(blocks))

    @inbounds for (i, bi) in enumerate(blocks)
        πX = Projection(X, bi)
        result[i] = overapproximate(πX, ST)
    end
    return CartesianProductArray(result)
end

# split the symmetric box [-1, 1]^n in nparts in each dimension
function _split_symmetric_box(n::Int, nparts::Int)
    return IA.mince(symBox(n), nparts)
end

# split the symmetric box [-1, 1]^n in partition[i] parts for each dimension
function _split_symmetric_box(D::Int, partition::Vector{Int})
    S = BallInf(zeros(D), 1.0)
    Sp = split(S, partition)
    return convert.(IA.IntervalBox, Sp)
end

# =========================
# Overapproximation
# =========================

# compared to LazySets.Approximations._overapproximate_hparallelotope,
# this function does inv(Matrix(Γ))
function _overapproximate_hparallelotope(Z::AbstractZonotope, indices=1:dim(Z))
    length(indices) == dim(Z) || throw(ArgumentError("the number of generator indices is $(length(indices)), " *
                                                     "but it was expected to be $(dim(Z))"))

    p, n = ngens(Z), dim(Z)
    if p == n
        return Z
    elseif p < n
        error("the zonotope order is $(order(Z)) but it should be at least 1")
    end

    G = genmat(Z)
    Γ = G[:, indices]
    □Γ⁻¹Z = box_approximation(linear_map(inv(Matrix(Γ)), Z))
    return linear_map(Γ, □Γ⁻¹Z)
end

function _overapproximate(lm::LinearMap{N, <:AbstractZonotope{N}, NM, <:AbstractIntervalMatrix{NM}},
                            ::Type{<:Zonotope}) where {N<:Real, NM}

    Mc, Ms = _split(matrix(lm))
    Z = LazySets.set(lm)
    c = center(Z)
    G = genmat(Z)
    _overapproximate_interval_linear_map(Mc, Ms, c, G)
end

function _overapproximate_interval_linear_map(Mc::AbstractMatrix{N},
                                              Ms::AbstractMatrix{N},
                                              c::AbstractVector{N},
                                              G::AbstractMatrix{N}) where {N}
    n = length(c)
    m = size(G, 2) # number of generators
    c_oa = Mc * c
    Ggens = Mc * G

    dvec = zeros(N, n)
    @inbounds for i in 1:n
        dvec[i] = abs(c[i])
        for j in 1:m
            dvec[i] += abs(G[i, j])
        end
    end
    q = size(Mc, 1)
    α = Ms * dvec # vector of length q
    αnz = findall(!iszero, α)
    DV = zeros(N, q, length(αnz))
    @inbounds for (j, idx) in enumerate(αnz)
        DV[j, idx] = α[idx]
    end
    G_oa = hcat(Ggens, DV)
    return Zonotope(c_oa, G_oa)
end

function _overapproximate_interval_linear_map(Mc::SMatrix{n, n, N, LM},
                                              Ms::SMatrix{n, n, N, LM},
                                              c::SVector{n, N},
                                              G::SMatrix{n, m, N, LG}) where {n, N, LM, m, LG}
    c_oa = Mc * c
    Ggens = Mc * G

    dvec = zeros(N, n)
    @inbounds for i in 1:n
        dvec[i] = abs(c[i])
        for j in 1:m
            dvec[i] += abs(G[i, j])
        end
    end
    q = size(Mc, 1)
    α = Ms * dvec # vector of length q
    αnz = findall(!iszero, α)
    DV = zeros(MMatrix{q, q, N}) # NOTE: sole difference with regular arrays, may refactor
    @inbounds for (j, idx) in enumerate(αnz)
        DV[j, idx] = α[idx]
    end
    G_oa = hcat(Ggens, DV)
    return Zonotope(c_oa, G_oa)
end

function _split_fallback!(A::IntervalMatrix{T}, C, S) where {T}
    m, n = size(A)
    @inbounds for j in 1:n
        for i in 1:m
            itv = A[i, j]
            radius = (sup(itv) - inf(itv)) / T(2)
            C[i, j] = inf(itv) + radius
            S[i, j] = radius
        end
    end
    return C, S
end

function _split(A::IntervalMatrix{T, IT, MT}) where {T, IT, MT<:AbstractMatrix{IT}}
    m, n = size(A)
    C = Matrix{T}(undef, m, n)
    S = Matrix{T}(undef, m, n)
    _split_fallback!(A, C, S)
    return C, S
end

function _split(A::IntervalMatrix{T, IT, MT}) where {T, IT, ST, MT<:StaticArray{ST, IT}}
    m, n = size(A)
    # TODO: use MMatrix and convert to SMatrix afterwards?
    C = Matrix{T}(undef, m, n)
    S = Matrix{T}(undef, m, n)
    _split_fallback!(A, C, S)
    return SMatrix{m, n, T}(C), SMatrix{m, n, T}(S)
end

_symmetric_interval_hull(x::Interval) = LazySets.symmetric_interval_hull(x)
_symmetric_interval_hull(x::Hyperrectangle) = LazySets.symmetric_interval_hull(x)

# type-stable version
function _symmetric_interval_hull(S::LazySet{N}) where {N}
    # fallback returns a hyperrectangular set
    (c, r) = LazySets.Approximations.box_approximation_helper(S)
    #if r[1] < 0
    #    return EmptySet{N}(dim(S))
    #end
    return Hyperrectangle(zeros(N, length(c)), abs.(c) .+ r)
end

# type-stable version
function _overapproximate(S::LazySet{N}, ::Type{<:Hyperrectangle}) where {N}
    c, r = LazySets.Approximations.box_approximation_helper(S)
    #if r[1] < 0
    #    return EmptySet{N}(dim(S))
    #end
    return Hyperrectangle(c, r)
end

# TEMP
function LazySets.Approximations.box_approximation(x::IntervalArithmetic.Interval)
    return convert(Hyperrectangle, Interval(x))
end

# TEMP
function LazySets.Approximations.box_approximation(x::IntervalArithmetic.IntervalBox)
    return convert(Hyperrectangle, x)
end

LazySets.box_approximation(S::UnionSetArray) = overapproximate(S, Hyperrectangle)

function LazySets.overapproximate(S::UnionSetArray{N}, ::Type{<:Hyperrectangle}) where {N}
    c, r = box_approximation_helper(S)
    if r[1] < 0
        return EmptySet{N}(dim(S))
    end
    return Hyperrectangle(c, r)
end

@inline function box_approximation_helper(S::UnionSetArray{N}) where {N}
    zero_N = zero(N)
    one_N = one(N)
    n = dim(S)
    c = Vector{N}(undef, n)
    r = Vector{N}(undef, n)
    d = zeros(N, n)
    @inbounds for i in 1:n
        d[i] = one_N
        htop = ρ(d, S)
        d[i] = -one_N
        hbottom = -ρ(d, S)
        d[i] = zero_N
        c[i] = (htop + hbottom) / 2
        r[i] = (htop - hbottom) / 2
        if r[i] < 0
            # contradicting bounds => set is empty
            # terminate with first radius entry being negative
            r[1] = r[i]
            break
        end
    end
    return c, r
end

# overapproximate a hyperrectangular set with a polytope
# TODO clean-up
_overapproximate(H::AbstractHyperrectangle, T::Type{<:HPolytope}) = _overapproximate_hyperrectangle(H, T)
_overapproximate(H::Hyperrectangle, T::Type{<:HPolytope}) = _overapproximate_hyperrectangle(H, T)
function _overapproximate_hyperrectangle(H, ::Type{<:HPolytope})
    P = overapproximate(H, Hyperrectangle)
    HPolytope([HalfSpace(Vector(c.a), c.b) for c in constraints_list(H)])
end

function Base.:(*)(M::AbstractMatrix, X::UnionSetArray{N, <:AbstractSingleton{N}}) where {N}
    return UnionSetArray([linear_map(M, p) for p in array(X)])
end

"""
    relative_error(x, x_ref)

Compute the relative error between interval `x` and a reference interval `xref`.

### Input

- `x`    -- interval
- `xref` -- reference interval

### Output

An interval representing the relative error (in percentage) of `x` with respect to
the reference interval `xref`.

### Algorithm

If ``x = [x_L, x_H]``` and ``xref = [xref_L, xref_H]``, the output is the interval
``y = 100 * [y_L, y_H]`` computed as ``y_L = -(x_L - xref_L) / den`` and
``y_H = (x_H - xref_H) / den``, where ``den = xref_H - xref_L``.

This function measures the relative error between an interval `x` and a reference
interval `x_ref` accounting for it the lower and the upper range bounds separately
(see  Eq. (20) in [1]).

### References

- [1] Althoff, Matthias, Dmitry Grebenyuk, and Niklas Kochdumper.
      "Implementation of Taylor models in CORA 2018."
      Proc. of the 5th International Workshop on Applied Verification for Continuous
      and Hybrid Systems. 2018. [pdf](https://easychair.org/publications/paper/9Tz3)
"""
function relative_error(x, x_ref)
    x_low, x_high = inf(x), sup(x)
    x_ref_low, x_ref_high = inf(x_ref), sup(x_ref)
    denom = x_ref_high - x_ref_low
    rel_low = -(x_low - x_ref_low) / denom
    rel_high = (x_high - x_ref_high) / denom
    return 100 * IntervalArithmetic.Interval(rel_low, rel_high)
end

# ==================================
# Zonotope order reduction methods
# ==================================

abstract type AbstractReductionMethod end

# These methods split the given zonotope Z into two zonotopes, K and L, where
# K contains the most "representative" generators and L contains the generators
# that are reduced, Lred using a box overapproximation
struct GIR05 <: AbstractReductionMethod end
struct COMB03 <: AbstractReductionMethod end

const _COMB03 = COMB03()
const _GIR05 = GIR05()

# algorithm selection
_reduce_order(Z::Zonotope, r::Number; force_reduction::Bool=false) = _reduce_order_GIR05(Z, r, force_reduction=force_reduction) # default
_reduce_order(Z::Zonotope, r::Number, ::GIR05; force_reduction::Bool=false) = _reduce_order_GIR05(Z, r, force_reduction=force_reduction)
_reduce_order(Z::Zonotope, r::Number, ::COMB03; force_reduction::Bool=false) = _reduce_order_COMB03(Z, r, force_reduction=force_reduction)

# zonotope with mixed static array types
function _reduce_order(Z::Zonotope{N, SVector{n, N}, MMatrix{n, p, N, L}},
                       r::Number, alg::GIR05; force_reduction::Bool=false) where {n, N, p, L}
    c = Z.center
    G = SMatrix(Z.generators)
    return _reduce_order(Zonotope(c, G), r, alg, force_reduction=force_reduction)
end

# return the indices of the generators in G (= columns) sorted according to the COMB03 method
# the generator index with highest score goes first
function _weighted_gens!(indices, G::AbstractMatrix{N}, ::COMB03) where {N}
    p = size(G, 2)
    weights = Vector{N}(undef, p)
    @inbounds for j in 1:p
        v = view(G, :, j)
        weights[j] = norm(v, 2)
    end
    sortperm!(indices, weights, rev=true, initialized=false)
    return indices
end

# return the indices of the generators in G (= columns) sorted according to the GIR05 method
# the generator index with highest score goes first
function _weighted_gens!(indices, G::AbstractMatrix{N}, ::GIR05) where {N}
    n, p = size(G)
    weights = Vector{N}(undef, p)
    @inbounds for j in 1:p
        aux_norm_1 = zero(N)
        aux_norm_inf = zero(N)
        for i in 1:n
            abs_Gij = abs(G[i, j])
            aux_norm_1 += abs_Gij
            if aux_norm_inf < abs_Gij
                aux_norm_inf = abs_Gij
            end
        end
        weights[j] = aux_norm_1 - aux_norm_inf
    end
    sortperm!(indices, weights, rev=true, initialized=false)
    return indices
end

# compute interval hull of the generators of G (= columns) corresponding to `indices`
function _interval_hull(G::AbstractMatrix{N}, indices) where {N}
    n, p = size(G)
    Lred = zeros(N, n, n)
    @inbounds for i in 1:n
        for j in indices
            Lred[i, i] += abs(G[i, j])
        end
    end
    return Lred
end

# implementation for static arrays
function _interval_hull(G::SMatrix{n, p, N, L}, indices) where {n, p, N, L}
    Lred = zeros(MMatrix{n, n, N})
    @inbounds for i in 1:n
        for j in indices
            Lred[i, i] += abs(G[i, j])
        end
    end
    return SMatrix{n, n}(Lred)
end

# given an n x p matrix G and a vector of m integer indices with m <= p,
# concatenate the columns of G given by `indices` with the matrix Lred
function _hcat_KLred(G::AbstractMatrix, indices, Lred::AbstractMatrix)
    K = view(G, :, indices)
    return hcat(K, Lred)
end

# implementation for static arrays
function _hcat_KLred(G::SMatrix{n, p, N, L1}, indices, Lred::SMatrix{n, n, N, L2}) where {n, p, N, L1, L2}
    m = length(indices)
    K = SMatrix{n, m}(view(G, :, indices))
    return hcat(K, Lred)
end

# Implements zonotope order reduction method from [COMB03]
# We follow the notation from [YS18]
function _reduce_order_COMB03(Z::Zonotope{N}, r::Number; force_reduction::Bool=false) where {N}
    r >= 1 || throw(ArgumentError("the target order should be at least 1, but it is $r"))
    c = Z.center
    G = Z.generators
    n, p = size(G)

    # r is bigger than the order of Z => don't reduce
    if !force_reduction
        (r * n >= p) && return Z
    end

    # this algorithm sort generators by decreasing 2-norm
    indices = Vector{Int}(undef, p)
    _weighted_gens!(indices, G, _COMB03)

    # the first m generators have greatest weight
    m = floor(Int, n * (r - 1))

    # compute interval hull of L
    Lred = _interval_hull(G, view(indices, (m+1):p))

    isone(r) && return Zonotope(c, Lred)

    Gred = _hcat_KLred(G, view(indices, 1:m), Lred)
    return Zonotope(c, Gred)
end

# Implements zonotope order reduction method from [GIR05]
# We follow the notation from [YS18]
function _reduce_order_GIR05(Z::Zonotope{N}, r::Number; force_reduction::Bool=false) where {N}
    r >= 1 || throw(ArgumentError("the target order should be at least 1, but it is $r"))
    c = Z.center
    G = Z.generators
    n, p = size(G)

    # r is bigger than the order of Z => don't reduce
    if !force_reduction
        (r * n >= p) && return Z
    end

    # this algorithm sorts generators by ||⋅||₁ - ||⋅||∞ difference
    indices = Vector{Int}(undef, p)
    _weighted_gens!(indices, G, _GIR05)

    # the first m generators have greatest weight
    m = floor(Int, n * (r - 1))

    # compute interval hull of L
    Lred = _interval_hull(G, view(indices, (m+1):p))

    isone(r) && return Zonotope(c, Lred)

    Gred = _hcat_KLred(G, view(indices, 1:m), Lred)
    return Zonotope(c, Gred)
end

# ====================================
# Disjointness checks
# ====================================

abstract type AbstractDisjointnessMethod end

# we pass the sets to the disjointness check algorithm without pre-processing
struct FallbackDisjointness <: AbstractDisjointnessMethod end

const NoEnclosure = FallbackDisjointness

# we overapproximate the reach-set with a zonotope, then make the disjointness check
struct ZonotopeEnclosure <: AbstractDisjointnessMethod end

# we overapproximate the reach-set with a hyperrectangle, then make the disjointness check
struct BoxEnclosure <: AbstractDisjointnessMethod end

# this is a dummy disjointness check which returns "false" irrespective of the value of its arguments
struct Dummy <: AbstractDisjointnessMethod end

# --------------------------------------------------------------------
# Methods to evaluate disjointness
# --------------------------------------------------------------------

# fallbacks
_is_intersection_empty(X::LazySet, Y::LazySet, ::FallbackDisjointness) = isdisjoint(X, Y)
_is_intersection_empty(X, Y) = _is_intersection_empty(X, Y, FallbackDisjointness())

# -----------------------------------------------
# Disjointness checks between specific set types
# -----------------------------------------------

using LazySets: _geq, _leq

# H : {x : ax <= b}, one-dimensional with a != 0
@commutative function _is_intersection_empty(X::Interval, H::HalfSpace)
    a = H.a[1]
    b = H.b
    N = promote_type(eltype(X), eltype(H))
    if a > zero(N)
        return !_leq(min(X), b/a)
    else
        return !_geq(max(X), b/a)
    end
end

# H : {x : ax = b}, one-dimensional with a != 0
@commutative function _is_intersection_empty(X::Interval, H::Hyperplane)
    q = H.b / H.a[1]
    return !_geq(q, min(X)) || !_leq(q, max(X))
end

# if X is polyhedral and Y is the set union of half-spaces, X ∩ Y is empty iff
# X ∩ Hi is empty for each half-space Hi in Y
# NOTE the algorithm below solves an LP for each X ∩ Hi; however, we can proceed
# more efficiently using support functions
# see LazySets.is_intersection_empty_helper_halfspace
@commutative function isdisjoint(X::AbstractPolytope{N},
                                 Y::UnionSetArray{N, <:HalfSpace{N}}) where {N}
    if dim(X) == 2 # use vrep in 2D
        Xp = convert(VPolygon, X)
        return all(Yi -> isdisjoint(Xp, Yi), array(Y))
    end

    clist_X = constraints_list(X)
    for ci in Y.array
        # using support functions
        #!(-ρ(-hs.a, X) > hs.b) && return false # TODO use robust check

        # using LP
        Y_ci = vcat(clist_X, ci)
        remove_redundant_constraints!(Y_ci) && return false
    end
    return true
end

# ====================================================================
# Concrete intersection
#
#  Methods to compute the intersection between two or more sets
# ====================================================================

# -----------------------
# Auxiliary functions
# -----------------------

# converts the normal vector of a list of half-spaces to be a Vector
const VECH{N, VT} = Vector{HalfSpace{N, VT}}
_to_vec(c::HalfSpace{N, Vector{N}}) where {N} = c
_to_vec(c::HalfSpace{N, VT}) where {N, VT<:AbstractVector{N}} = HalfSpace(Vector(c.a), c.b)
_to_vec(x::VECH{N, Vector{N}}) where {N} = x
_to_vec(x::VECH{N, VT}) where {N, VT<:AbstractVector{N}} = broadcast(_to_vec, x)

# concatenates lists of half-spaces such that the normal vectors of the final list
# are all Vector
_vcat(args::VECH...) = vcat(map(_to_vec, args)...)

# ------------------------
# FallbackIntersection
# ------------------------

abstract type AbstractIntersectionMethod end

# "fallback" implementation that uses LazySets intersection(X, Y)
struct FallbackIntersection{T} <: AbstractIntersectionMethod
     backend::T
end

_intersection(X::LazySet, Y::LazySet, ::FallbackIntersection) = intersection(X, Y)
_intersection(X, Y) = _intersection(X, Y, FallbackIntersection())

FallbackIntersection() = FallbackIntersection(nothing)

has_backend(alg::FallbackIntersection) = !isnothing(alg.backend)

function _intersection(X::AbstractPolyhedron{N}, Y::AbstractPolyhedron{N}, alg::FallbackIntersection) where {N}
    if has_backend(alg)
        return intersection(X, Y, backend=alg.backend)
    else
        return intersection(X, Y)
    end
end

# ------------------------
# HRepIntersection
# ------------------------

# evaluate X ∩ Y exactly using the constraint representation of X and Y
# evaluate X₁ ∩ ⋯ ∩ Xₖ using the constraint representation of each Xᵢ
#
# TODO Annotate normal vector types?
# struct HRepIntersection{SX, SY} <: AbstractIntersectionMethod end
# setrep(int_method::HRepIntersection{SX<:AbstractPolytope}, SY<:AbstractPolyhedron}) = SX
#
struct HRepIntersection <: AbstractIntersectionMethod
    #
end

setrep(::HRepIntersection) = HPolytope{Float64, Vector{Float64}}

function _intersection(X::AbstractPolyhedron, Y::AbstractPolyhedron, ::HRepIntersection)
    clist_X = constraints_list(X)
    clist_Y = constraints_list(Y)
    out = _vcat(clist_X, clist_Y)
    success = remove_redundant_constraints!(out)
    return (success, out)
end

function _intersection(X::AbstractPolyhedron, Y::AbstractPolyhedron, Z::AbstractPolyhedron, ::HRepIntersection)
    clist_X = constraints_list(X)
    clist_Y = constraints_list(Y)
    clist_Z = constraints_list(Z)
    out = _vcat(clist_X, clist_Y, clist_Z)
    success = remove_redundant_constraints!(out)
    return (success, out)
end

function _intersection(X::LazySet, Y::AbstractPolyhedron, Z::AbstractPolyhedron, W::AbstractPolyhedron, ::HRepIntersection)
    clist_X = constraints_list(X)
    clist_Y = constraints_list(Y)
    clist_Z = constraints_list(Z)
    clist_W = constraints_list(W)
    out = _vcat(clist_X, clist_Y, clist_Z, clist_W)
    success = remove_redundant_constraints!(out)
    return (success, out)
end

# ------------------------
# BoxIntersection
# ------------------------

struct BoxIntersection <: AbstractIntersectionMethod
    #
end

setrep(::BoxIntersection) = Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}

# ------------------------
# TemplateHullIntersection
# ------------------------

# evaluate X ∩ Y approximately using support function evaluations
#
# if lazy = false (default) use the property that the support function
# of X ∩ Y along direction d is not greater
# ρ(d, X ∩ Y) <= min(ρ(d, X), ρ(d, Y))
# by the same token, compute X₁ ∩ ⋯ ∩ Xₖ approximately using the same property
# if the template is provided, we have TN<:AbstractDirections{N, VN}
# otherwise, the constraints of X and Y are used and TN is Missing
#
# if lazy = true, use specialized approximation of lazy intersections
# assuming that Y is polyhedral
struct TemplateHullIntersection{N, VN, TN, L} <: AbstractIntersectionMethod
    dirs::TN
    lazy::L
end

# constructor with template directions provided
function TemplateHullIntersection(dirs::TN; lazy=false) where {N, VN, TN<:AbstractDirections{N, VN}}
    lazy_val = Val(lazy)
    TemplateHullIntersection{N, VN, TN, typeof(lazy_val)}(dirs, lazy_val)
end

# constructor without template directions => directions are missing until evaluated
function TemplateHullIntersection{N, VN}(; lazy=false) where {N, VN<:AbstractVector{N}}
    lazy_val = Val(lazy)
    TemplateHullIntersection{N, VN, Missing, typeof(lazy_val)}(missing, lazy_val)
end

function TemplateHullIntersection(; lazy=false)
    TemplateHullIntersection{Float64, Vector{Float64}}(; lazy=lazy)
end

setrep(::TemplateHullIntersection{N, VN}) where {N, VN} = HPolytope{N, VN}
setrep(::TemplateHullIntersection{N, SEV}) where {N, SEV<:SingleEntryVector{N}} = Union{HPolytope{N, SEV}, HPolytope{N, Vector{N}}}
setrep(::TemplateHullIntersection{N, SP}) where {N, SP<:SparseVector{N}} = Union{HPolytope{N, SP}, HPolytope{N, Vector{N}}}

# if the template directions is missing => use the constraints of both X and Y
# doesn't remove redundant constraints
function _intersection(X::LazySet, Y::LazySet, method::TemplateHullIntersection{N, VN, Missing, Val{false}}) where {N, VN}
    clist_X = constraints_list(X)
    clist_Y = constraints_list(Y)

    out = Vector{HalfSpace{N, VN}}()
    sizehint!(out, length(clist_X) + length(clist_Y))

    @inbounds for (i, c) in enumerate(clist_X)
        d = convert(VN, c.a) # normal vector
        b = min(c.b, ρ(d, Y)) # we know that ρ(d, X) = b
        push!(out, HalfSpace(d, b))
    end

    @inbounds for (i, c) in enumerate(clist_Y)
        d = convert(VN, c.a) # normal vector
        b = min(ρ(d, X), c.b) # we know that ρ(d, Y) = b
        push!(out, HalfSpace(d, b))
    end

    return HPolytope(out)
end

# use ρ(d, X∩Y) ≤ min(ρ(d, X), ρ(d, Y)) for each d in the template
function _intersection(X::LazySet, Y::LazySet, method::TemplateHullIntersection{N, VN, TN, Val{false}}) where {N, VN, TN<:AbstractDirections{N, VN}}
    dirs = method.dirs
    out = Vector{HalfSpace{N, VN}}(undef, length(dirs))
    @inbounds for (i, d) in enumerate(dirs)
        d = convert(VN, d)
        b = min(ρ(d, X), ρ(d, Y))
        out[i] = HalfSpace(d, b)
    end
    return HPolytope(out)
end

# compute the min of X ∩ Hi for each Hi in Y (assuming the second set is polyhedral,
# requires the list of constraints of Y) for each template direction d
function _intersection(X::LazySet, Y::LazySet, method::TemplateHullIntersection{N, VN, TN, Val{true}}) where {N, VN, TN<:AbstractDirections{N, VN}}
    dirs = method.dirs
    out = Vector{HalfSpace{N, VN}}(undef, length(dirs))
    H = constraints_list(Y)

    @inbounds for (i, d) in enumerate(dirs)
        d = convert(VN, d)
        b = minimum(ρ(d, X ∩ H[j]) for j in 1:length(H))
        out[i] = HalfSpace(d, b)
    end
    return HPolytope(out)
end

# =====================================
# Methods for checking inclusion
# =====================================

# given two reach-sets A, B check whether f(A) ⊆ g(B) where
# A ⊆ f(A) and B ⊆ g(B)
abstract type AbstractInclusionMethod end

# no pre-processing of the sets
struct FallbackInclusion <: AbstractInclusionMethod
#
end

_iscontained(X::LazySet, Y::LazySet, ::FallbackInclusion) = LazySets.⊆(X, Y)
_iscontained(X, Y) = _iscontained(X, Y, FallbackInclusion())


text:
#
# Potential energy function
#

function upair(r,data)
  r6 = r^6
  r12 = r6^2
  return data.ε4*(data.σ12/r12 - data.σ6/r6)
end


text:
JULIA_ENV = realpath(joinpath(@__DIR__, "..", ".."))
data_root = realpath(joinpath(@__DIR__, "..", "..", "data"))
data_input_root = joinpath(data_root, "input", "synthetic")
data_parameter_file = joinpath(data_input_root, "parameters.jser")
data_output_root = joinpath(data_root, "output")

worker_list = [("localhost", 1)]
exeflags = `--project="$JULIA_ENV"`
sshflags= `-i path/to/ssh/key/file`
SAMPLE_SIZE_LIMIT = 10000

fmt_string = "[{name} | {date} | {level}]: {msg}"
loglevel = "debug"

### data directories ###
data_dirs = ["clusters", "dims", "sizes"]

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation