text:
module LegendrePolynomials

using OffsetArrays
using SpecialFunctions

export Pl
export collectPl
export collectPl!
export Plm
export collectPlm
export collectPlm!
export dnPl
export collectdnPl
export collectdnPl!

# cache values of logfactorial to avoid expensive calls
const lfdictInt = Vector{Dict{Int, Float64}}()
const lfdictBigInt = Vector{Dict{BigInt, BigFloat}}()

function __init__()
    empty!(lfdictInt)
    empty!(lfdictBigInt)
    resize!(lfdictInt, Threads.nthreads())
    resize!(lfdictBigInt, Threads.nthreads())
    for i in 1:Threads.nthreads()
        lfdictInt[i] = Dict{Int, Float64}()
        lfdictBigInt[i] = Dict{BigInt, BigFloat}()
    end
end

checkdomain(x) = true
function checkdomain(x::Number)
    f = abs(x) > 1
    if f isa Bool
        f && throw(DomainError(x,
        "Legendre Polynomials are defined for arguments lying in -1 ⩽ x ⩽ 1"))
    end
    return nothing
end

assertnonnegative(l::Integer) = l >= 0 || throw(ArgumentError("l must be >= 0, received " * string(l)))

checklm(l::Integer, m::Integer) = -l <= m <= l || throw(ArgumentError("m must satisfy -l <= m <= l"))

function _checkvalues(x, l, m, n)
    checkdomain(x)
    assertnonnegative(l)
    checklm(l, m)
    n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
end

function checklength(arr, minlength)
    length(arr) >= minlength || throw(ArgumentError(
        "array is not large enough to store all values, require a minimum length of " * string(minlength)))
end

@inline function Pl_recursion(::Type{T}, ℓ, (Plm1, Plm2, α), x) where {T}
  # relation is valid from ℓ = 1
  Pl = ((2ℓ-1) * x * Plm1 - (ℓ-1) * Plm2)/ℓ
  convert(T, Pl), α
end

@inline function Plm_recursion(::Type{T}, ℓ, m, (Plm1, Plm2, α_ℓm1_m), x) where {T}
    # relation is valid from ℓ = m+2
    α_ℓ_m = coeff(ℓ, m)
    Pl = α_ℓ_m * (x * Plm1 - inv(α_ℓm1_m) * Plm2)
    return convert(T, Pl), α_ℓ_m
end

@inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, P_n_lm2, x) where {T}
  P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1) - (ℓ - 1) * P_n_lm2)/ℓ
  convert(T, P_n_l)
end

# special cases
# case 1 : l == n, in which case there are no (l-1,n) and (l-2,n) terms
@inline function dPl_recursion(::Type{T}, ℓ, n, ::Nothing, P_nm1_lm1, ::Nothing, x) where {T}
    P_n_l = (2ℓ-1) * P_nm1_lm1
    convert(T, P_n_l)
end
# case 1 : l == n + 1, in which case there's no (l-2,n) term
@inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, ::Nothing, x) where {T}
    P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1))/ℓ
    convert(T, P_n_l)
end

# Allow memoization
_logfactorial(n) = logfactorial(n)
_logfactorial(n::Int) = get!(lfdictInt[Threads.threadid()], n) do
                            logfactorial(n)
                        end
_logfactorial(n::BigInt) = get!(lfdictBigInt[Threads.threadid()], n) do
                            logfactorial(n)
                        end

function logplm_norm(l, m)
    T = promote_type(typeof(l), typeof(m))
    (log(T(2)) - log(2T(l)+1) + _logfactorial(l + m) - _logfactorial(l - m))/2
end
function _maybebigexp(t)
    if t < log(prevfloat(typemax(t)))
        return exp(t)
    else
        return exp(big(t))
    end
end
function plm_norm(l, m)
    t = logplm_norm(l, m)
    _maybebigexp(t)
end
function logabspll_prefactor(l, T = typeof(l))
    lT = T(l)
    a = (logfactorial(2lT-1) + log(2lT+1) - log(lT))/2
    b = l*log(T(2)) + _logfactorial(lT-1)
    a - b
end
neg1pow(l) = iseven(l) ? 1 : -1
function pll_prefactor(l, T = typeof(l); csphase::Bool = true)
    l == 0 && return sqrt(oftype(_logfactorial(T(l)), 1/2))
    t = logabspll_prefactor(l, T)
    (csphase ? neg1pow(l) : 1) * exp(t)
end
function coeff(l, m, T = float(promote_type(typeof(l), typeof(m))))
    √((2T(l) - 1)/((T(l)-m)*(T(l)+m))*(2T(l)+1))
end

polytype(x) = typeof(float.(x))
polytype(x::Array{<:Number}) = Array{promote_type(eltype(x), Float64), ndims(x)}
polytype(x::Number) = float(typeof(x))
polytype(x::Number, m::Number) = promote_type(polytype(x), polytype(m))
polytype(x, m::Number) = typeof(float.(x) .* m)
polytype(x::Array{<:Number}, m::Number) = Array{promote_type(eltype(x), Float64, typeof(m)), ndims(x)}
polytype(x, ::Nothing) = polytype(x)

struct LegendrePolynomialIterator{T, M<:Union{Integer, Nothing}, V}
    x :: V
    m :: M
    csphase :: Bool
    function LegendrePolynomialIterator{T,M,V}(x::V, m::M, csphase::Bool) where {T, M<:Union{Integer, Nothing}, V}
        checkdomain(x)
        new{T,M,V}(x, m, csphase)
    end
end
_checkm(::Nothing) = true
_checkm(m::Integer) = m >= 0
function LegendrePolynomialIterator(x, m::Union{Integer, Nothing} = nothing; csphase::Bool = true)
    @assert _checkm(m) "m must be non-negative"
    LegendrePolynomialIterator{polytype(x, m), typeof(m), typeof(x)}(x, m, csphase)
end

Base.eltype(::Type{<:LegendrePolynomialIterator{T}}) where {T} = T
Base.IteratorSize(::Type{<:LegendrePolynomialIterator}) = Base.IsInfinite()

# Iteration for Legendre polynomials
# starting value, l == 0 for m == 0
function Base.iterate(iter::LegendrePolynomialIterator{T,Nothing}) where {T}
    Pl = convert(T, one(iter.x))
    Plm1 = convert(T, zero(iter.x))
    nextstate = (Pl, Plm1, nothing)
    return Pl, (1, nextstate)
end
function Base.iterate(iter::LegendrePolynomialIterator{T,Nothing}, state) where {T}
    l, prevstate = state
    x = iter.x
    # standard-norm Bonnet iteration
    Pl, α = Pl_recursion(T, l, prevstate, x)
    Plm1 = first(prevstate)
    nextstate = (Pl, Plm1, α)
    return Pl, (l+1, nextstate)
end

# Iteration for associated Legendre polynomials
# starting value, l == m
function Base.iterate(iter::LegendrePolynomialIterator{T,<:Integer}) where {T}
    m = iter.m
    x = iter.x
    if m == 0
        Pl = convert(T, one(x))
    else
        f = pll_prefactor(m, csphase = iter.csphase)
        t = f * (√(one(x) - x^2))^m
        Pl = convert(T, t)
    end
    Plm1 = convert(T, zero(x))
    α_ℓ_m = Inf
    return Pl, (m+1, (Pl, Plm1, α_ℓ_m))
end

function Base.iterate(iter::LegendrePolynomialIterator{T,<:Integer}, state) where {T}
    l, prevstate = state
    x = iter.x
    m = iter.m
    if m == 0
        Pl, α_ℓ_m = Pl_recursion(T, l, prevstate, x)
    else
        Pl, α_ℓ_m = Plm_recursion(T, l, m, prevstate, x)
    end
    Plm1 = first(prevstate)
    nextstate = (Pl, Plm1, α_ℓ_m)
    return Pl, (l+1, nextstate)
end

throw_normerr(norm) =
    throw(ArgumentError("norm = $norm undefined,"*
        " valid norms are Val(:standard), Val(:normalized), Val(:schmidtquasi) or Val(:schmidt)"))

maybenormalize(P, l, ::Val{:normalized}) = oftype(P, P * √(l+1/2))
maybenormalize(P, l, ::Union{Val{:standard}, Val{:schmidtquasi}}) = P
maybenormalize(P, l, ::Val{:schmidt}) = oftype(P, P * √(2l+1))
maybenormalize(P, l, norm) = throw_normerr(norm)

plm_norm(::Val{:normalized}, lm...) = 1.0
plm_norm(::Val{:standard}, lm...) = plm_norm(lm...)
plm_norm(::Val{:schmidtquasi}, l, m) = √(2*(2 - iszero(m))/(2l + 1))
function plm_norm(::Val{:schmidt}, l, m)
    two = first(promote(2, l, m))
    √(two*(two - iszero(m)))
end
plm_norm(norm, args...) = throw_normerr(norm)

function maybenormalize(P, l, m, norm, csphase = true)
    if m == 0
        return maybenormalize(P, l, norm)
    else
        return (m < 0 ? neg1pow(m) : 1) * plm_norm(norm, l, m) * P
    end
end

"""
    Pl(x, l::Integer; [norm = Val(:standard)])

Compute the Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and the degree `l`.

The default norm is chosen to be `Val(:standard)`, in which case the polynomials
satisfy `Pl(1, l) == 1` for all `l`.
Optionally, the `norm` may be set to `Val(:normalized)` to obtain normalized Legendre polynomials.
These have an L2 norm of `1`.

# Examples
```jldoctest
julia> Pl(1, 2)
1.0

julia> Pl(0.5, 4) ≈ -37/128 # analytically obtained value
true

julia> Pl(0.5, 20, norm = Val(:normalized))
-0.21895188261094017
```
"""
function Pl(x, l::Integer; norm = Val(:standard))
    _checkvalues(x, l, 0, 0)
    iter_full = LegendrePolynomialIterator(x)
    iter = Iterators.drop(iter_full, l)
    Pl = first(iter)
    return maybenormalize(Pl, l, norm)
end

function doublefactorial(T, n)
    p = convert(T, one(n))
    for i in n:-2:1
        p *= convert(T, i)
    end
    convert(T, p)
end

"""
    Plm(x, l::Integer, m::Integer; [kwargs...])

Compute the associated Legendre polynomial ``P_\\ell^m(x)`` for degree `l`
and order ``m`` lying in `-l:l`.

For `m == 0` this function returns Legendre polynomials.

# Optional keyword arguments
* `norm`: The normalization used in the function. Possible
    options are `Val(:standard)` (default), `Val(:normalized)`, `Val(:schmidtquasi)` or `Val(:schmidt)`.
    * The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
    * The functions obtained with `norm = Val(:schmidtquasi)` have an L2 norm of ``\\sqrt{\\frac{2(2-\\delta_{m0})}{2l+1}}``.
    * The functions obtained with `norm = Val(:schmidt)` have an L2 norm of ``\\sqrt{2(2-\\delta_{m0})}``.
* `csphase::Bool`: The Condon-shortley phase ``(-1)^m``, which is included by default.

# Examples

```jldoctest
julia> Plm(0.5, 3, 2) ≈ 45/8 # analytically obtained value
true

julia> Plm(0.5, 4, 0) == Pl(0.5, 4)
true
```
"""
function Plm(x, l::Integer, m::Integer; norm = Val(:standard), csphase::Bool = true)
    _checkvalues(x, l, m, 0)
    lT, mT = promote(l, m)
    iter_full = LegendrePolynomialIterator(x, abs(mT); csphase)
    iter = Iterators.drop(iter_full, l - abs(m))
    Plm = first(iter)
    return maybenormalize(Plm, lT, mT, norm, csphase)
end

Base.@propagate_inbounds function _unsafednPl!(cache, x, l, n)
    # unsafe, assumes 1-based indexing
    checklength(cache, l - n + 1)
    if n == l # may short-circuit this
        cache[1] = doublefactorial(eltype(cache), 2l-1)
    else
        collectPl!(cache, x, lmax = l - n)

        for ni in 1:n
            # We denote the terms as P_ni_li

            # li == ni
            P_nim1_nim1 = cache[1]
            P_ni_ni = dPl_recursion(eltype(cache), ni, ni, nothing, P_nim1_nim1, nothing, x)
            cache[1] = P_ni_ni

            # li == ni + 1
            P_nim1_ni = cache[2]
            P_ni_nip1 = dPl_recursion(eltype(cache), ni + 1, ni, P_ni_ni, P_nim1_ni, nothing, x)
            cache[2] = P_ni_nip1

            for li in ni+2:min(l, l - n + ni)
                P_ni_lim2 = cache[li - ni - 1]
                P_ni_lim1 = cache[li - ni]
                P_nim1_lim1 = cache[li - ni + 1]
                P_ni_li = dPl_recursion(eltype(cache), li, ni, P_ni_lim1, P_nim1_lim1, P_ni_lim2, x)
                cache[li - ni + 1] = P_ni_li
            end
        end
    end
    nothing
end

"""
    dnPl(x, l::Integer, n::Integer, [cache::AbstractVector])

Compute the ``n``-th derivative ``d^n P_\\ell(x)/dx^n`` of the Legendre polynomial ``P_\\ell(x)``.
Optionally a pre-allocated vector `cache` may be provided, which must have a minimum length of `l - n + 1`
and may be overwritten during the computation.

The order of the derivative `n` must be non-negative. For `n == 0` this function just returns
Legendre polynomials.

# Examples

```jldoctest
julia> dnPl(0.5, 3, 2) # second derivative of P3(x) at x = 0.5
7.5

julia> dnPl(0.5, 4, 0) == Pl(0.5, 4) # zero-th order derivative == Pl(x)
true
```
"""
Base.@propagate_inbounds function dnPl(x, l::Integer, n::Integer,
    A = begin
        _checkvalues(x, l, 0, n)
        # do not allocate A if the value is trivially zero
        if l < n
            return zero(polytype(x))
        end
        zeros(polytype(x), l - n + 1)
    end
    )

    _checkvalues(x, l, 0, n)
    # check if the value is trivially zero in case A is provided in the function call
    if l < n
        return zero(eltype(A))
    end

    cache = OffsetArrays.no_offset_view(A)
    # function barrier, as no_offset_view may be type-unstable
    _unsafednPl!(cache, x, l, n)

    return cache[l - n + 1]
end

"""
    collectPl!(v::AbstractVector, x; [lmin::Integer = 0], [lmax::Integer = length(v) - 1 + lmin], [kwargs...])

Compute the Legendre Polynomials ``P_\\ell(x)`` for the argument `x` and degrees `l = lmin:lmax`,
and store them in `v`.

At output `v[firstindex(v) + l] == Pl(x, l; kwargs...)`.

# Optional keyword arguments
* `norm`: The normalization used in the function. Possible
    options are `Val(:standard)` (default) and `Val(:normalized)`.
    The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.

# Examples
```jldoctest
julia> v = zeros(4);

julia> collectPl!(v, 0.5);

julia> all(zip(0:3, v)) do (l, vl); vl ≈ Pl(0.5, l); end
true

julia> collectPl!(v, 0.5, norm = Val(:normalized));

julia> all(zip(0:3, v)) do (l,vl); vl ≈ Pl(0.5, l, norm = Val(:normalized)); end
true

julia> v = zeros(0:4);

julia> collectPl!(v, 0.5, lmax = 3) # only l from 0 to 3 are populated
5-element OffsetArray(::Vector{Float64}, 0:4) with eltype Float64 with indices 0:4:
  1.0
  0.5
 -0.125
 -0.4375
  0.0
```
"""
function collectPl!(v::AbstractVector, x; lmin::Integer = 0,
        lmax::Integer  = length(v) - 1 + lmin, norm = Val(:standard))
    _checkvalues(x, lmin, 0, 0)
    lmin <= lmax || throw(ArgumentError("lmin must be less than or equal to lmax"))
    n = lmax - lmin + 1
    checklength(v, n)

    lmaxT, lminT = promote(lmax, lmin)
    iter = Iterators.take(Iterators.drop(LegendrePolynomialIterator(x), lmin), n)

    inds = (firstindex(v)-1) .+ (1:n)
    v_section = @view v[inds]

    for (ind, Pl) in enumerate(iter)
        l = ind - 1 + lminT
        v_section[ind] = maybenormalize(Pl, l, norm)
    end

    return v
end

"""
    collectPl(x; lmax::Integer, [lmin::Integer = 0], [kwargs...])

Compute the Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and degrees `l = lmin:lmax`.
Return `P` with indices `lmin:lmax`, where `P[l] == Pl(x, l; kwargs...)`

# Optional keyword arguments
* `norm`: The normalization used in the function. Possible
    options are `Val(:standard)` (default) and `Val(:normalized)`.
    The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.

# Examples
```jldoctest
julia> v = collectPl(0.5, lmax = 4);

julia> all(l -> v[l] ≈ Pl(0.5, l), 0:4)
true

julia> v = collectPl(0.5, lmax = 4, norm = Val(:normalized));

julia> all(l -> v[l] ≈ Pl(0.5, l, norm = Val(:normalized)), 0:4)
true
```
"""
function collectPl(x; lmax::Integer, lmin::Integer = 0, norm = Val(:standard))
    _checkvalues(x, lmin, 0, 0)
    lmin <= lmax || throw(ArgumentError("lmin must be less than or equal to lmax"))
    Tnorm = promote_type(typeof(lmax), typeof(lmin))
    T = promote_type(polytype(x), Tnorm)
    v_offset = OffsetVector{T}(undef, lmin:lmax)
    v = parent(v_offset)
    collectPl!(v, x; lmax, lmin, norm)
    return v_offset
end

"""
    collectPlm(x; m::Integer, lmax::Integer, [lmin::Integer = abs(m)], [kwargs...])

Compute the associated Legendre polynomials ``P_\\ell^m(x)`` for the argument `x`,
degrees `l = lmin:lmax` and order `m` lying in `-l:l`.

Returns `v` with indices `lmin:lmax`, with `v[l] == Plm(x, l, m; kwargs...)`.

# Optional keyword arguments
* `norm`: The normalization used in the function. Possible
    options are `Val(:standard)` (default), `Val(:normalized)`, `Val(:schmidtquasi)` or `Val(:schmidt)`.
    * The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
    * The functions obtained with `norm = Val(:schmidtquasi)` have an L2 norm of ``\\sqrt{\\frac{2(2-\\delta_{m0})}{2l+1}}``.
    * The functions obtained with `norm = Val(:schmidt)` have an L2 norm of ``\\sqrt{2(2-\\delta_{m0})}``.
* `csphase::Bool`: The Condon-shortley phase ``(-1)^m``, which is included by default.
* `Tnorm`: Type of the normalization factor, which is dynamicall inferred by default,
    and is used in allocating an appropriate container.
    In case it is known that the norm may be expressed as a certain type without overflow
    (eg. `Float64`), this may be provided. Care must be taken to avoid overflows if `Tnorm`
    is passed as an argument.

# Examples

```jldoctest
julia> v = collectPlm(0.5, lmax = 3, m = 2);

julia> all(l -> v[l] ≈ Plm(0.5, l, 2), 2:3)
true

julia> v = collectPlm(0.5, lmax = 3, m = 2, norm = Val(:normalized));

julia> all(l -> v[l] ≈ Plm(0.5, l, 2, norm = Val(:normalized)), 2:3)
true
```
"""
function collectPlm(x; lmax::Integer, m::Integer, norm = Val(:standard),
        lmin::Integer = abs(m),
        Tnorm = begin
            norm === Val(:standard) ? begin
            _checkvalues(x, lmax, m, 0)
            typeof(plm_norm(lmax, m))
            end : float(mapreduce(typeof, promote_type, (m, lmax, lmin)))
        end,
        csphase::Bool = true,
        )
    _checkvalues(x, lmin, m, 0)
    lmin <= lmax || throw(ArgumentError("lmin must be less than or equal to lmax"))
    T = promote_type(polytype(x), Tnorm)
    v = OffsetVector{T}(undef, lmin:lmax)
    collectPlm!(parent(v), x; lmin, lmax, m, norm, csphase)
    v
end

"""
    collectPlm!(v::AbstractVector, x; m::Integer, [lmin::Integer = abs(m)], [lmax::Integer = length(v) - 1 + lmin], [kwargs...])

Compute the associated Legendre polynomials ``P_\\ell^m(x)`` for the argument `x`,
degrees `l = lmin:lmax` and order `m` lying in `-l:l`, and store the result in `v`.

At output, `v[l + firstindex(v)] == Plm(x, l, m; kwargs...)` for `l = lmin:lmax`.

# Optional keyword arguments
* `norm`: The normalization used in the function. Possible
    options are `Val(:standard)` (default), `Val(:normalized)`, `Val(:schmidtquasi)` or `Val(:schmidt)`.
    * The functions obtained with `norm = Val(:normalized)` have an L2 norm of ``1``.
    * The functions obtained with `norm = Val(:schmidtquasi)` have an L2 norm of ``\\sqrt{\\frac{2(2-\\delta_{m0})}{2l+1}}``.
    * The functions obtained with `norm = Val(:schmidt)` have an L2 norm of ``\\sqrt{2(2-\\delta_{m0})}``.
* `csphase::Bool`: The Condon-shortley phase ``(-1)^m``, which is included by default.

# Examples

```jldoctest
julia> v = zeros(2);

julia> collectPlm!(v, 0.5, lmax = 3, m = 2);

julia> all(zip(2:3, v)) do (l, vl); vl ≈ Plm(0.5, l, 2); end
true
```
"""
function collectPlm!(v, x; m::Integer, lmin::Integer = abs(m),
        lmax::Integer = length(v) - 1 + lmin,
        norm = Val(:standard),
        csphase::Bool = true,
        )

    lmin <= lmax || throw(ArgumentError("lmin must be less than or equal to lmax"))
    checklm(lmin, m)
    n = lmax - lmin + 1
    checklength(v, n)

    lminT, lmaxT, mT = promote(lmin, lmax, m)
    iter = Iterators.drop(LegendrePolynomialIterator(x, abs(mT); csphase), lmin - abs(m))

    inds = (firstindex(v)-1) .+ (1:n)
    v_section = @view v[inds]

    for (ind, Plm) in enumerate(iter)
        l = ind - 1 + lminT
        v_section[ind] = maybenormalize(Plm, l, mT, norm, csphase)
        ind == n && break
    end

    return v
end

"""
    collectdnPl(x; lmax::Integer, n::Integer)

Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`.

The order of the derivative `n` must be greater than or equal to zero.

Returns `v` with indices `0:lmax`, where `v[l] == dnPl(x, l, n)`.

# Examples

```jldoctest
julia> collectdnPl(0.5, lmax = 3, n = 2)
4-element OffsetArray(::Vector{Float64}, 0:3) with eltype Float64 with indices 0:3:
 0.0
 0.0
 3.0
 7.5
```
"""
function collectdnPl(x; lmax::Integer, n::Integer)
    assertnonnegative(lmax)
    n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
    v = zeros(polytype(x), 0:lmax)
    if lmax >= n
        collectdnPl!(parent(v), x; lmax = lmax, n = n)
    end
    v
end

"""
    collectdnPl!(v::AbstractVector, x; lmax::Integer, n::Integer)

Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`,
and store the result in `v`.

The order of the derivative `n` must be greater than or equal to zero.

At output, `v[l + firstindex(v)] == dnPl(x, l, n)` for `l = 0:lmax`.

# Examples

```jldoctest
julia> v = zeros(4);

julia> collectdnPl!(v, 0.5, lmax = 3, n = 2)
4-element Vector{Float64}:
 0.0
 0.0
 3.0
 7.5
```
"""
function collectdnPl!(v, x; lmax::Integer, n::Integer)
    assertnonnegative(lmax)
    n >= 0 || throw(ArgumentError("order of derivative n must be >= 0"))
    checklength(v, lmax + 1)

    # trivially zero for l < n
    fill!((@view v[(0:n-1) .+ firstindex(v)]), zero(eltype(v)))
    # populate the other elements
    @inbounds dnPl(x, lmax, n, @view v[(n:lmax) .+ firstindex(v)])

    v
end

end


text:
export normpdf, normcdfCody1990, erfCody, erfcCody, erfcxCody

const a = [3.1611237438705656, 113.864154151050156, 377.485237685302021, 3209.37758913846947, 0.185777706184603153]
const b = [23.6012909523441209, 244.024637934444173, 1282.61652607737228, 2844.23683343917062]
const c__ = [
    0.564188496988670089,
    8.88314979438837594,
    66.1191906371416295,
    298.635138197400131,
    881.95222124176909,
    1712.04761263407058,
    2051.07837782607147,
    1230.33935479799725,
    2.15311535474403846e-8,
]
const d__ = [
    15.7449261107098347,
    117.693950891312499,
    537.181101862009858,
    1621.38957456669019,
    3290.79923573345963,
    4362.61909014324716,
    3439.36767414372164,
    1230.33935480374942,
]
const p = [
    0.305326634961232344,
    0.360344899949804439,
    0.125781726111229246,
    0.0160837851487422766,
    6.58749161529837803e-4,
    0.0163153871373020978,
]
const q =
    [2.56852019228982242, 1.87295284992346047, 0.527905102951428412, 0.0605183413124413191, 0.00233520497626869185]

const sqrpi = 0.56418958354775628695
const thresh = 0.46875
const sixten = 16.0
# Explanation of machine-dependent constants */
#   XMIN   = the smallest positive floating-point number. */
#   XINF   = the largest positive finite floating-point number. */
#   XNEG   = the largest negative argument acceptable to ERFCX */
#            the negative of the solution to the equation */
#            2*exp(x*x) = XINF. */
#   XSMALL = argument below which erf(x) may be represented by */
#            2*x/sqrt(pi)  and above which  x*x  will not underflow. */
#            A conservative value is the largest machine number X */
#           such that   1.0 + X = 1.0   to machine precision. */
#   XBIG   = largest argument acceptable to ERFC  solution to */
#            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where */
#            W(x) = Math.exp(-x*x)/[x*sqrt(pi)]. */
#   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to */
#            machine precision.  A conservative value is */
#            1/[2*sqrt(XSMALL)] */
#   XMAX   = largest acceptable argument to ERFCX the minimum */
#           of XINF and 1/[sqrt(pi)*XMIN]. */
# The numbers below were preselected for IEEE .
const xinf = 1.79e308
const xneg = -26.628
const xsmall = 1.11e-16
const xbig = 26.543
const xhuge = 6.71e7
const xmax = 2.53e307

const norm_cdf_asymptotic_expansion_first_threshold = -10.0
const norm_cdf_asymptotic_expansion_second_threshold = -1 / sqrt(eps())

function normpdf(x::T)::T where {T}
    sqrtTwoPi = sqrt(2*T(pi))
    return exp(-x^2 / 2) /sqrtTwoPi
end
function normpdf(x::Float64)::Float64
    return exp(-x^2 / 2) * OneOverSqrt2Pi
end
function normcdfCody1990(z::Float64)::Float64
    if z <= norm_cdf_asymptotic_expansion_first_threshold
        # Asymptotic expansion for very negative z following (26.2.12) on page 408
        # in M. Abramowitz and A. Stegun, Pocketbook of Mathematical Functions, ISBN 3-87144818-4.
        sum = 1.0
        if z >= norm_cdf_asymptotic_expansion_second_threshold
            zsqr = z * z
            i = 1
            g = 1.0
            x = floatmax()
            y = floatmax()
            a = floatmax()
            lasta = 0.0
            while (true)
                lasta = a
                x = Float64(4 * i - 3) / zsqr
                y = x * (Float64(4 * i - 1) / zsqr)
                a = g * (x - y)
                sum -= a
                g *= y
                i += 1
                a = abs(a)
                if lasta > a && a >= abs(sum * eps())
                    break
                end
            end
        end
        return -normpdf(z) * sum / z
    end
    return 0.5 * erfc(-z * OneOverSqrtTwo)
end

function d_int(x::Float64)::Float64
    if x > 0
        return floor(x)
    else
        return floor(-x)
    end
end

function calerf(x::Float64, jint::Int)::Float64

    y = 0.0
    del = 0.0
    ysq = 0.0
    xden = 0.0
    xnum = 0.0
    result = 0.0

    y = abs(x)
    if y <= thresh

        ysq = 0.0
        if y > xsmall
            ysq = y * y
        end
        xnum = a[5] * ysq
        xden = ysq
        for i__ = 1:3
            xnum = (xnum + a[i__]) * ysq
            xden = (xden + b[i__]) * ysq
        end
        result = x * (xnum + a[4]) / (xden + b[4])
        if jint != 0
            result = 1.0 - result
        end
        if jint == 2
            result = exp(ysq) * result
        end

        return result

    elseif y <= 4.0
        xnum = c__[9] * y
        xden = y
        for i__ = 1:7
            xnum = (xnum + c__[i__]) * y
            xden = (xden + d__[i__]) * y

        end
        result = (xnum + c__[8]) / (xden + d__[8])
        if jint != 2
            d__1 = y * sixten
            ysq = d_int(d__1) / sixten
            del = (y - ysq) * (y + ysq)
            d__1 = exp(-ysq * ysq) * exp(-del)
            result = d__1 * result
        end
    else
        result = 0.0
        if y >= xbig
            if jint != 2 || y >= xmax
                return l300(jint, x, result)
            end
            if y >= xhuge
                result = sqrpi / y
                return l300(jint, x, result)
            end
        end
        ysq = 1.0 / (y * y)
        xnum = p[6] * ysq
        xden = ysq
        for i__ = 1:4
            xnum = (xnum + p[i__]) * ysq
            xden = (xden + q[i__]) * ysq
        end
        result = ysq * (xnum + p[5]) / (xden + q[5])
        result = (sqrpi - result) / y
        if jint != 2
            d__1 = y * sixten
            ysq = d_int(d__1) / sixten
            del = (y - ysq) * (y + ysq)
            d__1 = exp(-ysq * ysq) * exp(-del)
            result = d__1 * result
        end
    end

    return l300(jint, x, result)
end

function l300(jint::Int, x::Float64, result::Float64)::Float64
    if jint == 0
        result = (0.5 - result) + 0.5
        if x < 0.0
            result = -(result)
        end
    elseif jint == 1
        if x < 0.0
            result = 2.0 - result
        end
    else
        if x < 0.0
            if x < xneg
                result = xinf
            else
                d__1 = x * sixten
                ysq = d_int(d__1) / sixten
                del = (x - ysq) * (x + ysq)
                y = exp(ysq * ysq) * exp(del)
                result = y + y - result
            end
        end
    end
    return result
end

function erfCody(x::Float64)::Float64
    # -------------------------------------------------------------------- */
    #This subprogram computes approximate values for erf(x). */
    #   (see comments heading CALERF). */
    #   Author/date: W. J. Cody, January 8, 1985 */

    return calerf(x, 0)
end

function erfcCody(x::Float64)::Float64

    return calerf(x, 1)
end

function erfcxCody(x::Float64)::Float64
    # ------------------------------------------------------------------ */
    # This subprogram computes approximate values for Math.exp(x*x) * erfc(x). */
    #   (see comments heading CALERF). */
    #   Author/date: W. J. Cody, March 30, 1987 */
    # ------------------------------------------------------------------ */

    return calerf(x, 2)
end

ext: jl


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

using Test
using Behavior.Selection
using Behavior.Gherkin
using Behavior.Gherkin: AbstractScenario

@testset "Selection            " begin
    # These tests check that the tag selector is parsed and that the expressions are used
    # properly. It uses only tags on scenarios to check this, for simplicity. Tests in a
    # section below check that tags are inherited properly between features and scenarios.
    @testset "Tag selector" begin
        @testset "Feature has tag @foo and no scenarios; Selecting @foo returns no scenarios" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @foo and one scenario; selecting @foo returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags and no scenarios; select @foo returns no scenarios" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], String[])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @bar; Select for @foo returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@bar"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @foo; Select for @bar returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@bar")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tags @bar and @foo; Select for @foo returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@bar", "@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags; Selector is (not @foo); select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tag @foo; Selector is (not @foo); select returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @foo; Empty selector; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags; Empty selector; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags; Selector is only whitespace; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("   ")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tag @baz; Selector is only whitespace; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("   ")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end
    end

    # These tests check that tags are inherited properly between features and scenarios.
    # It assumes that the parsing of tag selectors works as intended, tested above.
    @testset "Feature and Scenario selection" begin
        @testset "Feature has tag @foo and one scenario, no tags; select returns the feature with one scenario" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[])
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Only one scenario has tag @bar; select returns the feature with only that scenario" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String["@bar"], ScenarioStep[]),
                Scenario("Some other scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@bar")
            newfeature = select(selector, feature)

            # Assert
            @test length(newfeature.scenarios) == 1
            @test newfeature.scenarios[1] == feature.scenarios[1]
        end

        @testset "One scenario has tag @ignore; selecting on (not @ignore) returns only the feature without that tag" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String["@ignore"], ScenarioStep[]),
                Scenario("Some other scenario", String["@other"], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @ignore")
            newfeature = select(selector, feature)

            # Assert
            @test length(newfeature.scenarios) == 1
            @test newfeature.scenarios[1] == feature.scenarios[2]
        end
    end

    @testset "Or expressions" begin
        @testset "Feature has tags @foo and one scenario; selecting @foo,@bar returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo,@bar")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tags @foo and one scenario; selecting @bar,@foo returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@bar,@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tag @foo; Selector is (not @foo,@bar); select returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo,@bar")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @bar; Selector is (not @foo,@bar); select returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@bar"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo,@bar")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end
    end
end

@testset "New Selection        " begin
    #
    # A feature template is created that shows which scenarios are expected
    # to be selected by the tag expression. Each `Yes` scenario is expected
    # to be selected, and each `No` scenario is expected to _not_ be selected.
    # This is all done for readability, so it's clear what each test actually tests.
    #
    # The `FeatureTemplate` struct contains all information needed to create an
    # actual Feature and its Scenarios, to send into the selection method.
    #

    abstract type SelectedOrNot end

    struct Yes <: SelectedOrNot
        tags::Vector{String}
    end

    struct No <: SelectedOrNot
        tags::Vector{String}
    end

    struct FeatureTemplate
        tags::Vector{String}
        selections::Vector{<:SelectedOrNot}
    end

    struct TV
        description::String
        expression::String
        featuretemplate::FeatureTemplate
    end

    function makescenarios(tv::TV)
        AbstractScenario[
            Scenario("Some scenario: $(i)", selectedornot.tags, ScenarioStep[])
            for (i, selectedornot) in enumerate(tv.featuretemplate.selections)
        ]
    end

    function makeyesscenarios(tv::TV)
        AbstractScenario[
            Scenario("Some scenario: $(i)", selectedornot.tags, ScenarioStep[])
            for (i, selectedornot) in enumerate(tv.featuretemplate.selections)
            if isa(selectedornot, Yes)
        ]
    end

    function makefeature(tv::TV)
        scenarios = makescenarios(tv)
        Feature(
            FeatureHeader("Some feature", String[], tv.featuretemplate.tags),
            scenarios)
    end

    testvectors = [
        TV(
            "Expression @foo: Feature has @foo, matches all scenarios",
            "@foo",
            FeatureTemplate(["@foo"],
                [
                    Yes(String[])
                ]
            )
        ),

        TV(
            "Expression @foo: Feature has @bar, matches no scenarios",
            "@foo",
            FeatureTemplate(["@bar"],
                [
                    No(String[])
                ]
            )
        )
    ]

    for tv in testvectors
        @testset "$(tv.description)" begin
            # Arrange
            feature = makefeature(tv)

            # Act
            selector = newparsetagselector(tv.expression)
            newfeature = select(selector, feature)

            # Assert
            expectedscenarios = makeyesscenarios(tv)
            @test newfeature.scenarios == expectedscenarios
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayBackendHealthHttpSettings <: SwaggerModel
    backendHttpSettings::Any # spec type: Union{ Nothing, ApplicationGatewayBackendHttpSettings } # spec name: backendHttpSettings
    servers::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayBackendHealthServer} } # spec name: servers

    function ApplicationGatewayBackendHealthHttpSettings(;backendHttpSettings=nothing, servers=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendHealthHttpSettings, Symbol("backendHttpSettings"), backendHttpSettings)
        setfield!(o, Symbol("backendHttpSettings"), backendHttpSettings)
        validate_property(ApplicationGatewayBackendHealthHttpSettings, Symbol("servers"), servers)
        setfield!(o, Symbol("servers"), servers)
        o
    end
end # type ApplicationGatewayBackendHealthHttpSettings

const _property_map_ApplicationGatewayBackendHealthHttpSettings = Dict{Symbol,Symbol}(Symbol("backendHttpSettings")=>Symbol("backendHttpSettings"), Symbol("servers")=>Symbol("servers"))
const _property_types_ApplicationGatewayBackendHealthHttpSettings = Dict{Symbol,String}(Symbol("backendHttpSettings")=>"ApplicationGatewayBackendHttpSettings", Symbol("servers")=>"Vector{ApplicationGatewayBackendHealthServer}")
Base.propertynames(::Type{ ApplicationGatewayBackendHealthHttpSettings }) = collect(keys(_property_map_ApplicationGatewayBackendHealthHttpSettings))
Swagger.property_type(::Type{ ApplicationGatewayBackendHealthHttpSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendHealthHttpSettings[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendHealthHttpSettings }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendHealthHttpSettings[property_name]

function check_required(o::ApplicationGatewayBackendHealthHttpSettings)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendHealthHttpSettings }, name::Symbol, val)
end


text:
@testset "NACA 4" begin
function read_airfoil_data(fn::AbstractString)
    xs = Float64[]
    ys = Float64[]
    open(fn, "r") do io
        title = readline(io)
        while true
            line = strip(readline(io))
            if isempty(line)
                break
            end
            x, y = parse.(Float64, split(line))
            push!(xs, x)
            push!(ys, y)
        end
    end
    if (iseven(length(xs)))
        error("File Broken")
    end
    n = (length(xs) + 1) >> 1
    xu = xs[n:-1:1]
    yu = ys[n:-1:1]
    xl = xs[n:end]
    yl = ys[n:end]
    xu, yu, xl, yl
end
xu, yu, xl, yl = read_airfoil_data("NACA0015.dat")
naca0015 = NACA"0015"
yuaf = y_upper(naca0015, xu)
ylaf = y_lower(naca0015, xl)
@test all(isapprox.(yu, yuaf, atol=1e-4))
@test all(isapprox.(yl, ylaf, atol=1e-4))
naca4312s = NACA"4312"s
xu, yu, xl, yl = read_airfoil_data("NACA4312S.dat")
xtest = collect(0:0.1:1)
xuaf = x_upper(naca4312s, xtest)
yuaf = y_upper(naca4312s, xtest)
xlaf = x_lower(naca4312s, xtest)
ylaf = y_lower(naca4312s, xtest)
@test all(isapprox.(xu, xuaf, atol=1e-4))
@test all(isapprox.(yu, yuaf, atol=1e-4))
@test all(isapprox.(xl, xlaf, atol=1e-4))
@test all(isapprox.(yl, ylaf, atol=1e-4))
x1, y1 = gen_airfoil(naca4312s, 65, orient=:CCW)
@test y1[end-1] < y1[2]
x2, y2 = gen_airfoil(naca4312s, 65, orient=:CW)
@test x1 ≈ reverse(x2) && y1 ≈ reverse(y2)
@test_throws ArgumentError gen_airfoil(naca4312s, 65, orient=:OTHER)
end


text:
using StaticPolynomials
using Base.Test
using BenchmarkTools

@polyvar x y z

@testset "monomials" begin
    m = x^2
    @test isa(m, Monomial)
    @test variables(m) == (x,)
    @test exponents(m) == (2,)

    m = x * y
    @test isa(m, Monomial)
    @test variables(m) == (x, y)
    @test exponents(m) == (1, 1)

    m = y * x
    @test isa(m, Monomial)
    @test variables(m) == (x, y)
    @test exponents(m) == (1, 1)

    m = y * x^2 * y^3
    @test isa(m, Monomial)
    @test variables(m) == (x, y)
    @test exponents(m) == (2, 4)

    @test x * y < x^2
    @test (@benchmark $(x * y) < $(x^2)).allocs == 0
end

@testset "terms" begin
    t = @inferred 3x
    @test isa(t, Term)
    @test variables(t) == (x,)
    @test coefficient(t) == 3
    @test monomial(t) == convert(Monomial, x)

    t = @inferred x * 3
    @test isa(t, Term)
    @test variables(t) == (x,)
    @test coefficient(t) == 3
    @test monomial(t) == convert(Monomial, x)

    t = @inferred (2x) * y * x
    @test isa(t, Term)
    @test variables(t) == (x, y)
    @test coefficient(t) == 2
    @test exponents(t) == (2, 1)
end

@testset "polynomials" begin
    p = @inferred x + 1
    @test isa(p, Polynomial)
    @test p.terms[2] == Term(x)
    @test variables(p) == (x,)

    p = @inferred x^2 + x * x
    @test isa(p, Polynomial)
    @test length(terms(p)) == 1
    @test coefficient(terms(p)[1]) == 2

    p = @inferred (1 + x) + (y + 2)
    @test coefficient(terms(p)[1]) == 3
    @test exponents(terms(p)[1]) == (0, 0)
    @test variables(p) == (x, y)
    @test coefficient(terms(p)[2]) == 1
    @test exponents(terms(p)[2]) == (0, 1)
    @test coefficient(terms(p)[3]) == 1
    @test exponents(terms(p)[3]) == (1, 0)

    p = @inferred x^2 + y + x * x + 3 * x * y + x * y^2
    @test length(terms(p)) == 4
    @test coefficient.(terms(p)) == [1, 3, 2, 1]
    @test exponents.(terms(p)) == [(0, 1), (1, 1), (2, 0), (1, 2)]
end

function testmonomials(var, degree)
    [var^i for i in 0:degree]
end

@testset "monomial vector" begin
    ms = @inferred testmonomials(x, 3)
    @test eltype(ms) == Monomial{1, (x,)}
    @test (@benchmark testmonomials($x, 3)).allocs == 1
end

include("substitution.jl")


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 9da308cd-880e-45ca-9c9f-79315cb575c6
using InteractiveUtils, PlutoUI, TickTock, ConfParser, Dates, Random, StaticArrays, Distributed, OrdinaryDiffEq, JLD2, Plots, LoopVectorization, BenchmarkTools

# ╔═╡ fe4cd7b2-b5d2-11eb-29f5-81ec8d2a9a9e
md"
## Particle Tracing Simulation of Non-Linear Resonant Electron Scattering by Chorus Waves
#### General notes

 - default ingestion of data will be via `.jld2` files in a particular subdirectory

#### There's plenty of plots

Lost electrons are cool.

"

# ╔═╡ 9834a489-ea0f-45b7-a299-74454b4b3ace
md"
conf file = $(@bind conffile TextField())

directory = $(@bind dirname TextField())

basename  = $(@bind basename TextField())

batches   = $(@bind num_batches NumberField(1:100; default=1))
"

# ╔═╡ 8fa0edf6-5276-436f-bb76-e4e90429d0f2
begin
	const Re   = 6370e3;        # Earth radius, f64
	const c    = 3e8;           # speedo lite, f64
	const Beq  = 3.e-5;         # B field at equator (T), f64

	# Parsing the conf file
	conf = ConfParse(dirname*"/"*conffile)
	parse_conf!(conf)
	numParticles = parse(Int64, retrieve(conf, "numberOfParticles"));
	startTime = parse(Float64, retrieve(conf, "startTime"));
	endTime = parse(Float64, retrieve(conf, "endTime"));
	tspan = (startTime, endTime); # integration time
	lossConeAngle = parse(Float64, retrieve(conf, "lossConeAngle"));
	saveDecimation = parse(Float64, retrieve(conf, "saveDecimation"));
	L = parse(Float64, retrieve(conf, "L"));
	omegam = parse(Float64, retrieve(conf, "omegam"));
	Omegape = parse(Float64, retrieve(conf, "Omegape"));
	z0 = parse(Float64, retrieve(conf, "z0"));
	lambda0 = parse(Float64, retrieve(conf, "lambda0"));
	waveAmplitudeModifier = parse(Float64, retrieve(conf, "waveAmplitudeModifier"));
	ELo = parse(Float64, retrieve(conf, "ELo"));
	EHi = parse(Float64, retrieve(conf, "EHi"));
	Esteps = parse(Float64, retrieve(conf, "Esteps"));
	PALo = parse(Float64, retrieve(conf, "PALo"));
	PAHi = parse(Float64, retrieve(conf, "PAHi"));
	PAsteps = parse(Float64, retrieve(conf, "PAsteps"));
	ICrange = [ELo, EHi, Esteps, PALo, PAHi, PAsteps];
	batches = parse(Int64, retrieve(conf, "batches"));
	numThreads = parse(Int64, retrieve(conf, "numberOfThreads"))
	@info "Parsed Config file: $conffile"
end

# ╔═╡ 59c01c3b-effc-4260-be8f-981feda851c8
begin

	calcb(lambda::Float64)  = sqrt(1+3*sin(lambda)^2)/(cos(lambda)^6)

	calcdb(lambda::Float64) = (3*(27*sin(lambda)-5*sin(3*lambda))) / (cos(lambda)^8*(4+12*sin(lambda)^2))

	calcGamma(pz::Float64,mu::Float64,b::Float64) = sqrt(1 + pz^2 + 2*mu*b)

	calcK(b::Float64,lambda::Float64) = (Omegape * (cos(lambda)^-(5/2)))/sqrt(b/omegam - 1)

	calcAlpha(mu::Float64, gamma::Float64) = rad2deg(asin(sqrt((2*mu)/(gamma^2 - 1))))

	function loadData(directory::String, basename::String, num_batches::Int64)
		#=
		Loads in the JLD2 file and creates timeseries of all the data.
		Takes in filename in string, returns 5 TxN vectors for time,
		z, pz, pitch angle, and energy timeseries for all N particles.
		=#

		allZ = Vector{Vector{Float64}}();
		allPZ = Vector{Vector{Float64}}();
		allE = Vector{Vector{Float64}}();
		allPA = Vector{Vector{Float64}}();
		allT = Vector{Vector{Float64}}();

		Threads.@threads for i in 1:num_batches
			JLD2.@load jldname = directory*"/"*basename*"_$i.jld2" sol

			for traj in sol # TODO make this multithreaded
				vars = hcat(traj.u...) # pulls out the canonical position/momentum
				# z        = vars[1,:];
				# pz       = vars[2,:];
				# zeta     = vars[3,:];
				# mu       = vars[4,:];
				# lambda   = vars[5,:];
				b = @views @avx @. sqrt(1+3*sin(vars[5,:])^2)/(cos(vars[5,:])^6)
				gamma = @views @avx @. sqrt(1 + vars[2,:]^2 + 2*vars[4,:]*b)

				@views push!(allT, traj.t);
				@views push!(allZ, vars[1,:]);
				@views push!(allPZ, vars[2,:]);
				@views push!(allPA, @. rad2deg( @avx asin(sqrt((2*vars[4,:])/(gamma^2 - 1)))));
				@views push!(allE, @. (511*(gamma - 1)));
			end
			@info "$(length(sol)) particles loaded in from $(basename*"_$i.jld2")"
		end
		@info "Loaded total $(length(allT)) particles!"
		return allZ, allPZ, allT, allPA, allE
	end

	function countLostParticles(allT)
		```
		Based on time vectors, counts which ones were lost
		and at what time. Returns a Nx2 array where the first
		column is the time at which the particle was lost, and
		the 2nd column denotes the number lost at that point.
		```
		lossCounter = []; # initialize vector

		for ntime in allT # loop thru each particle
			if maximum(ntime) != endTime # particle lost if ended early
				push!(lossCounter, maximum(ntime)); # the final entry in the time vector is when the particle got lost
			end
		end
		if isempty(lossCounter) # if particle wasn't lost, then throw in a NaN
			push!(lossCounter, NaN)
			@warn "All particles trapped!" # since all particles trapped
		end
		lostParticles = [sort(lossCounter) collect(1:length(lossCounter))] # sort it by time, so you can see when particles got lost

		if maximum(lostParticles[:,1]) != endTime # this adds in a final hline from last particle lost to end of simulation
			lostParticles = vcat(lostParticles, [endTime (maximum(lostParticles[:,2]))]);
		end
		#lostParticles[1:end .!= 6,: ] # clever 1 liner to rid a row
		@info "Total of $(lostParticles[end,end]) particles lost during sim"
		return lostParticles
	end

end

# ╔═╡ de04d879-e805-4d27-8383-6398fb44872d
allZ, allPZ, allT, allPA, allE = loadData(dirname, basename, num_batches);

# ╔═╡ 7af09807-4012-46e2-bb2e-d307fa092b28
# @time countLostParticles(allT)

# ╔═╡ dbd59420-c97c-48e0-81e4-57bd72a5e52b

# ╔═╡ Cell order:
# ╟─fe4cd7b2-b5d2-11eb-29f5-81ec8d2a9a9e
# ╠═9da308cd-880e-45ca-9c9f-79315cb575c6
# ╟─9834a489-ea0f-45b7-a299-74454b4b3ace
# ╠═8fa0edf6-5276-436f-bb76-e4e90429d0f2
# ╠═59c01c3b-effc-4260-be8f-981feda851c8
# ╠═de04d879-e805-4d27-8383-6398fb44872d
# ╠═7af09807-4012-46e2-bb2e-d307fa092b28
# ╠═dbd59420-c97c-48e0-81e4-57bd72a5e52b
# ╠═27be5594-8ad4-48e4-8a3d-c06c9e712311

max_stars_repo_name: ethantsai/nlwhistlers

max_issues_repo_name: ethantsai/nlwhistlers

max_forks_repo_name: ethantsai/nlwhistlers


text:
const ageTypes = Dict( "age" => :age, "tenure" => :tenure,
    "time in node" => :timeInNode )
const ageTypesS = Dict( :age => "age", :tenure => "tenure",
    :timeInNode => "time in current node" )

function writeInitialPopData( report::DataFrame, mpSim::MPsim,
    timeFactor::Float64, reportGenerationTime::Float64, filename::String,
    overwrite::Bool )

    tStart = now()
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[1]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header.
        ws["A1"] = "Simulation length"
        ws["B1"] = mpSim.simLength / timeFactor
        ws["C1"] = "years"
        ws["A2"] = "Data generation time"
        ws["B2"] = reportGenerationTime
        ws["C2"] = "seconds"
        ws["A3"] = "Excel generation time"
        ws["C3"] = "seconds"

        if isempty( report )
            ws["A5"] = "No initial population in simulation"
            excelGenerationTime = ( now() - tStart ).value / 1000.0
            ws["B3"] = excelGenerationTime
            return
        end  # if isempty( report )

        # Initial population report.
        ws["A5"] = "Node"
        ws["B5"] = "# in initial population"
        nNodes = size( report, 1 )

        for ii in 1:nNodes
            ws[ii + 5, 1] = report[ii, 1]
            ws[ii + 5, 2] = report[ii, 2]
        end  # for ii in 1:nNodes

        ws[nNodes + 6, 1] = "Total"
        ws[nNodes + 6, 2] = sum( report[:, 2] )

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws["B3"] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # writeInitialPopData( report, mpSim, timeFactor, reportGenerationTime,
     #   filename, overwrite )

function writeInitialPopAgeData( report::DataFrame, mpSim::MPsim,
    ageType::Symbol, timeFactor::Float64, reportGenerationTime::Float64,
    filename::String, overwrite::Bool )

    tStart = now()
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[1]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header.
        ws["A1"] = "Simulation length"
        ws["B1"] = mpSim.simLength / timeFactor
        ws["C1"] = "years"
        ws["A2"] = "Report on"
        ws["B2"] = ageTypesS[ageType]
        ws["A3"] = "Data generation time"
        ws["B3"] = reportGenerationTime
        ws["C3"] = "seconds"
        ws["A4"] = "Excel generation time"
        ws["C4"] = "seconds"

        if isempty( report )
            ws["A6"] = "No initial population in simulation"
            excelGenerationTime = ( now() - tStart ).value / 1000.0
            ws["B4"] = excelGenerationTime
            return
        end  # if isempty( report )

        ws["A5"] = "Initial population size"
        ws["B5"] = sum( report[1:(end-5), 2] )

        ws["A7"] = "Summary statistics"
        ws["A8"] = "Mean"
        ws["A9"] = "Median"
        ws["A10"] = "Std. Dev."
        ws["A11"] = "Minimum"
        ws["A12"] = "Maximum"

        for ii in 1:5
            val = report[end - (5-ii), 2]
            ws[ii + 7, 2] = val isa Missing ? "N/A" : val / timeFactor
        end  # for ii in 1:5

        ws["A14"] = "Age"
        ws["A15"] = "Amount"
        nBrackets = size( report, 1 ) - 5

        for ii in 1:nBrackets
            ws[14, ii + 1] = report[ii, 1] / timeFactor
            ws[15, ii + 1] = report[ii, 2]
        end  # for ii in 1:nBrackets

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws["B4"] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # writeInitialPopAgeData( report, mpSim, ageType, timeFactor,
     #   reportGenerationTime, filename, overwrite )


text:
@testset "get_interface_addrs" begin
    results = AWSClusterManagers.get_interface_addrs()

    @test results isa Vector{AWSClusterManagers.InterfaceAddress}
    @test length(results) > 0
end


text:
using MLJ
import RDatasets: dataset
using PrettyPrinting
import DataFrames: DataFrame, select, Not
DTC = @load DecisionTreeClassifier pkg=DecisionTree

carseats = dataset("ISLR", "Carseats")

first(carseats, 3) |> pretty

High = ifelse.(carseats.Sales .<= 8, "No", "Yes") |> categorical;
carseats[!, :High] = High;

X = select(carseats, Not([:Sales, :High]))
y = carseats.High;

HotTreeClf = OneHotEncoder() |> DTC()

mdl = HotTreeClf
mach = machine(mdl, X, y)
fit!(mach);

ypred = predict_mode(mach, X)
misclassification_rate(ypred, y)

train, test = partition(eachindex(y), 0.5, shuffle=true, rng=333)
fit!(mach, rows=train)
ypred = predict_mode(mach, rows=test)
misclassification_rate(ypred, y[test])

r_mpi = range(mdl, :(decision_tree_classifier.max_depth), lower=1, upper=10)
r_msl = range(mdl, :(decision_tree_classifier.min_samples_leaf), lower=1, upper=50)

tm = TunedModel(
    model=mdl,
    ranges=[r_mpi, r_msl],
    tuning=Grid(resolution=8),
    resampling=CV(nfolds=5, rng=112),
    operation=predict_mode,
    measure=misclassification_rate
)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

ypred = predict_mode(mtm, rows=test)
misclassification_rate(ypred, y[test])

fitted_params(mtm).best_model.decision_tree_classifier

DTR = @load DecisionTreeRegressor pkg=DecisionTree

boston = dataset("MASS", "Boston")

y, X = unpack(boston, ==(:MedV))

train, test = partition(eachindex(y), 0.5, shuffle=true, rng=551);

scitype(X)

X = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))

dtr_model = DTR()
dtr = machine(dtr_model, X, y)

fit!(dtr, rows=train)

ypred = MLJ.predict(dtr, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

r_depth = range(dtr_model, :max_depth, lower=2, upper=20)

tm = TunedModel(
    model=dtr_model,
    ranges=[r_depth],
    tuning=Grid(resolution=10),
    resampling=CV(nfolds=5, rng=1254),
    measure=rms
)

mtm = machine(tm, X, y)

fit!(mtm, rows=train)

ypred = MLJ.predict(mtm, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

RFR = @load RandomForestRegressor pkg=ScikitLearn

rf_mdl = RFR()
rf = machine(rf_mdl, X, y)
fit!(rf, rows=train)

ypred = MLJ.predict(rf, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

XGBR = @load XGBoostRegressor

xgb_mdl = XGBR(num_round=10, max_depth=10)
xgb = machine(xgb_mdl, X, y)
fit!(xgb, rows=train)

ypred = MLJ.predict(xgb, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using netgengen
using Base.Test

teststr =
"""
solid unioni2 = not orthobrick (0.0, 0.0, 0.0; 1.0, 1.0, 1.0) or
	cylinder (0.0, 0.0, 0.0; 1.0, 1.0, 1.0; 1.0);
tlo unioni2 -col=[1,0,0];
solid laatikko = orthobrick (0.0, 0.0, 0.0; 1.0, 1.0, 1.0);
solid loota = orthobrick (0.0, 0.0, 0.0; 1.0, 0.5, 2.0);
solid unioni = laatikko or
	loota;
tlo unioni -col=[0,1,0];
"""

laatikko = brick("laatikko", [0.0,0.0,0.0],[1.0,1.0,1.0])
loota = brick("loota", [0.0,0.0,0.0],[1.0,0.5,2.0])
unioni2 = csgunion("unioni2", [not(brick([0.0,0.0,0.0],[1.0,1.0,1.0])), cylinder([0.0,0.0,0.0],[1.0,1.0,1.0], 1)])
unioni = csgunion("unioni", [laatikko, loota])

io=IOBuffer()

tlo(unioni2, io, col=[1, 0, 0])
tlo(unioni, io, col=[0, 1, 0])
outstr = takebuf_string(io)

@test outstr == teststr


text:
"""
`process_pokes`
"""

function process_pokes(filepath::String)
    curr_data = FileIO.load(filepath) |> DataFrame
    if !in(:Poke,names(curr_data))
        DataFrames.rename!(curr_data, names(curr_data)[1] => :Poke) #change poke counter name
    end
    if in(:delta,names(curr_data))
        DataFrames.rename!(curr_data, :delta => :Delta)
    end
    curr_data[!,:Poke] = curr_data[!,:Poke].+1
    start_time = curr_data[1,:PokeIn]
    curr_data[!,:PokeIn] = curr_data[!,:PokeIn] .- start_time
    curr_data[!,:PokeOut] = curr_data[!,:PokeOut] .- start_time
    curr_data[!,:PokeDur] = curr_data[!,:PokeOut] - curr_data[!,:PokeIn]
    if !iscolumn(curr_data,:Wall)
        curr_data[:Wall] = zeros(size(curr_data,1))
    end
    booleans=[:Reward,:Side,:SideHigh,:Stim,:Wall]#columns to convert to Bool
    for x in booleans
        if curr_data[1,x] isa AbstractString
                curr_data[!,x] = parse.(Bool,curr_data[!,x])
        elseif curr_data[1,x] isa Real
            curr_data[!,x] = Bool.(curr_data[:,x])
        end
    end
    curr_data[!,:Side] = [a ? "L" : "R" for a in curr_data[!,:Side]]
    if iscolumn(curr_data,:ProbVec0)
        integers=[:Protocollo,:ProbVec0,:ProbVec1,:GamVec0,:GamVec1,:Delta]; #columns to convert to Int64
        for x in integers
            curr_data[!,x] = Int64.(curr_data[!,x])
        end
        curr_data[!,:Protocol] = Flipping.get_protocollo(curr_data)
        for x in[:ProbVec0,:ProbVec1,:GamVec0,:GamVec1,:Protocollo]
            DataFrames.select!(curr_data,DataFrames.Not(x))
        end
        if !iscolumn(curr_data,:StimFreq)
            curr_data[!,:StimFreq] .= 10000
        end
        curr_data[!,:StimFreq] = [a == 10000 ? 25 : a  for a in curr_data[!,:StimFreq]]
        curr_data[!,:Box] .= "Box0"
    elseif iscolumn(curr_data,:Prwd)
        curr_data[!,:Protocol] = string.(curr_data[!,:Prwd],'/',curr_data[!,:Ptrs])
        curr_data[!,:Box] = "Box".*string.(curr_data[:,:Box])
    end
    mouse, day, daily_session, session = session_info(filepath)
    curr_data[!,:MouseID] .= mouse
    curr_data[!,:Day] .= parse(Int64,day)
    curr_data[!,:Daily_Session] .= daily_session
    curr_data[!,:Session] .= session
    curr_data[!,:Gen] = Flipping.gen.(curr_data[:,:MouseID])
    curr_data[!,:Drug] = Flipping.pharm.(curr_data[:,:Day])
    curr_data[!,:Stim_Day] .= length(findall(curr_data[:,:Stim])) == 0 ? false : true
    curr_data[!,:Streak] = count_sequence(curr_data[:,:Side])
    curr_data[!,:ReverseStreak] = reverse(curr_data[:,:Streak])
    curr_data[!,:Poke_within_Streak] .= 0
    curr_data[!,:Poke_Hierarchy] .= 0.0
    curr_data[!,:Poke_within_Streak] = Vector{Union{Float64,Missing}}(undef,size(curr_data,1))
    curr_data[!,:Pre_Interpoke] = Vector{Union{Float64,Missing}}(undef,size(curr_data,1))
    curr_data[!,:Post_Interpoke] = Vector{Union{Float64,Missing}}(undef,size(curr_data,1))
    curr_data[!,:LastPoke] .= false
    by(curr_data,:Streak) do dd
        dd[:,:Poke_within_Streak] = count_sequence(dd[!,:Poke])
        dd[:,:Pre_Interpoke] =  dd[!,:PokeIn] .-lag(dd[!,:PokeOut],default = missing)
        dd[:,:Post_Interpoke] = lead(dd[!,:PokeIn],default = missing).- dd[!,:PokeOut]
        dd[:,:Poke_Hierarchy] = Flipping.get_hierarchy(dd[!,:Reward])
        dd.LastPoke[end] = true
    end
    curr_data[!,:Block] = count_sequence(curr_data[!,:Wall])
    curr_data[!,:Streak_within_Block] .= 0
    by(curr_data,:Block) do dd
        dd[:,:Streak_within_Block] = count_sequence(dd[!,:Side])
    end
    curr_data[!,:SideHigh] = [x ? "L" : "R" for x in curr_data[!,:SideHigh]]
    curr_data[!,:Correct] = curr_data[!,:Side] .== curr_data[!,:SideHigh]
    return curr_data
end


text:
using LoopVectorization, Test

const GAMMA = 1.4
@testset "Reject Unroll" begin
  function r_manual_inline!(r, u)
    @inbounds @fastmath begin
      m = size(u, 2)
      n = size(u, 3)
      fill!(r, 0.0)
      for j = 2:n-1
        for i = 1:m-1
          rl = u[1, i, j]
          rul = u[2, i, j]
          rvl = u[3, i, j]
          rEl = u[4, i, j]

          rr = u[1, i+1, j]
          rur = u[2, i+1, j]
          rvr = u[3, i+1, j]
          rEr = u[4, i+1, j]

          rr1 = 1.0 / rr
          ur = rur * rr1
          vr = rvr * rr1
          Er = rEr * rr1
          u2r = ur * ur + vr * vr
          pr = (GAMMA - 1.0) * (rEr - 0.5 * rr * u2r)
          hr = Er + pr * rr1
          unr = ur

          rhr = rEr + pr

          rl1 = 1.0 / rl
          ul = rul * rl1
          vl = rvl * rl1
          El = rEl * rl1
          u2l = ul * ul + vl * vl
          pl = (GAMMA - 1.0) * (rEl - 0.5 * rl * u2l)
          hl = El + pl * rl1
          unl = ul

          rhl = rEl + pl

          di = sqrt(rr * rl1)
          d1 = 1.0 / (di + 1.0)
          ui = (di * ur + ul) * d1
          vi = (di * vr + vl) * d1
          hi = (di * hr + hl) * d1
          ci2 = (GAMMA - 1.0) * (hi - 0.5 * (ui * ui + vi * vi))
          ci = sqrt(ci2)
          af = 0.5 * (ui * ui + vi * vi)
          uni = ui

          dr = rr - rl
          dru = rur - rul
          drv = rvr - rvl
          drE = rEr - rEl

          rlam1 = abs(uni + ci)
          rlam2 = abs(uni - ci)
          rlam3 = abs(uni)

          s1 = 0.5 * (rlam1 + rlam2)
          s2 = 0.5 * (rlam1 - rlam2)
          al1x = (GAMMA - 1.0) * (af * dr - ui * dru - vi * drv + drE)
          al2x = -uni * dr + dru
          cc1 = ((s1 - rlam3) * al1x / ci2) + (s2 * al2x / ci)
          cc2 = (s2 * al1x / ci) + (s1 - rlam3) * al2x

          f1 = 0.5 * (rr * unr + rl * unl) - 0.5 * (rlam3 * dr + cc1)
          f2 =
            0.5 * (rur * unr + rul * unl) + 0.5 * (pr + pl) -
            0.5 * (rlam3 * dru + cc1 * ui + cc2)
          f3 = 0.5 * (rvr * unr + rvl * unl) - 0.5 * (rlam3 * drv + cc1 * vi)
          f4 = 0.5 * (rhr * unr + rhl * unl) - 0.5 * (rlam3 * drE + cc1 * hi + cc2 * uni)

          r[1, i, j] -= f1
          r[2, i, j] -= f2
          r[3, i, j] -= f3
          r[4, i, j] -= f4

          r[1, i+1, j] += f1
          r[2, i+1, j] += f2
          r[3, i+1, j] += f3
          r[4, i+1, j] += f4
        end
      end
    end
    r
  end

  function r_turbo!(r, u)
    @inbounds begin
      m = size(u, 2)
      n = size(u, 3)
      fill!(r, 0.0)
      @turbo for j = 2:n-1
        for i = 1:m-1
          rl = u[1, i, j]
          rul = u[2, i, j]
          rvl = u[3, i, j]
          rEl = u[4, i, j]

          rr = u[1, i+1, j]
          rur = u[2, i+1, j]
          rvr = u[3, i+1, j]
          rEr = u[4, i+1, j]

          rr1 = 1.0 / rr
          ur = rur * rr1
          vr = rvr * rr1
          Er = rEr * rr1
          u2r = ur * ur + vr * vr
          pr = (GAMMA - 1.0) * (rEr - 0.5 * rr * u2r)
          hr = Er + pr * rr1
          unr = ur

          rhr = rEr + pr

          rl1 = 1.0 / rl
          ul = rul * rl1
          vl = rvl * rl1
          El = rEl * rl1
          u2l = ul * ul + vl * vl
          pl = (GAMMA - 1.0) * (rEl - 0.5 * rl * u2l)
          hl = El + pl * rl1
          unl = ul

          rhl = rEl + pl

          di = sqrt(rr * rl1)
          d1 = 1.0 / (di + 1.0)
          ui = (di * ur + ul) * d1
          vi = (di * vr + vl) * d1
          hi = (di * hr + hl) * d1
          ci2 = (GAMMA - 1.0) * (hi - 0.5 * (ui * ui + vi * vi))
          ci = sqrt(ci2)
          af = 0.5 * (ui * ui + vi * vi)
          uni = ui

          dr = rr - rl
          dru = rur - rul
          drv = rvr - rvl
          drE = rEr - rEl

          rlam1 = abs(uni + ci)
          rlam2 = abs(uni - ci)
          rlam3 = abs(uni)

          s1 = 0.5 * (rlam1 + rlam2)
          s2 = 0.5 * (rlam1 - rlam2)
          al1x = (GAMMA - 1.0) * (af * dr - ui * dru - vi * drv + drE)
          al2x = -uni * dr + dru
          cc1 = ((s1 - rlam3) * al1x / ci2) + (s2 * al2x / ci)
          cc2 = (s2 * al1x / ci) + (s1 - rlam3) * al2x

          f1 = 0.5 * (rr * unr + rl * unl) - 0.5 * (rlam3 * dr + cc1)
          f2 =
            0.5 * (rur * unr + rul * unl) + 0.5 * (pr + pl) -
            0.5 * (rlam3 * dru + cc1 * ui + cc2)
          f3 = 0.5 * (rvr * unr + rvl * unl) - 0.5 * (rlam3 * drv + cc1 * vi)
          f4 = 0.5 * (rhr * unr + rhl * unl) - 0.5 * (rlam3 * drE + cc1 * hi + cc2 * uni)

          r[1, i, j] -= f1
          r[2, i, j] -= f2
          r[3, i, j] -= f3
          r[4, i, j] -= f4

          r[1, i+1, j] += f1
          r[2, i+1, j] += f2
          r[3, i+1, j] += f3
          r[4, i+1, j] += f4
        end
      end
    end
    r
  end

  u = rand(4, 128, 128)
  u[[1, 4], :, :] .= 1.0
  u[[2, 3], :, :] .= 0.1
  @test r_manual_inline!(similar(u), u) ≈ r_turbo!(similar(u), u)

end


text:
export PolyZeros
import Polynomials: roots, degree

"""
    PolyZeros(z, mult)

represents a set of distinct (potential) roots z[j] of a polynomial with their
multiplicities. If the same root value appears repeatedly, the surplus occurrences
are discarded and the multiplicity addapted accordingly.
"""
struct PolyZeros{S<:Number}
    z::AbstractVector{S}		# the distinct root values
    mult::AbstractArray{Int}	# the multiplicities of the roots

    function PolyZeros(zz::AbstractVector{S}, ll::AbstractVector{<:Integer} = Int[]) where S<:Number
        zz, ll = filter_duplicates(zz, ll)
        new{S}(zz, ll)
    end
end

function PolyZeros(A::AbstractMatrix{T}) where T
    m, n = size(A)
    n == 2 || throw(ArgumentError("need exactly 2 colums"))
    zz = A[:,1]
    ll = Int.(A[:,2])
    PolyZeros(zz, ll)
end

Base.isreal(z::PolyZeros) = z.real_coefficients
multiplicities(z::PolyZeros) = z.mult
roots(z::PolyZeros) = z.z
degree(z::PolyZeros) = sum(abs.(z.mult))

"""
  `lessabs(a,b) -> Bool`

Order by absolute value of a and b.
For equal absolute values first negative real part, then positive imaginary part first.

For example 0 < -1 < 1 < -4+3im < -4-3im < -3-4im < 3+4im < 5

"""
function lessabs(a::S, b::T) where {S<:Number,T<:Number}
    diff = abs2(a) - abs2(b)
    Z = zero(diff)
    if diff != Z || ( isreal(a) && isreal(b) )
      return diff < Z || ( diff == Z && real(a) < real(b) )
    end
    real(a) < real(b) || ( real(a) == real(b) && imag(a) > imag(b) )
end

function less_real_first(a::S, b::T) where {S<:Number, T<:Number}
    ra, rb = real(a), real(b)
    ia, ib = imag(a), imag(b)
    ra < rb ? true : ra > rb ? false : abs(ia) < abs(ib) ? true : abs(ia) > abs(ib) ? false : ia < ib
end

"""
    zerostructure(zz::PolyZero) -> Type, Vector{Int}
    zerostructure(z::Vector, mult::Vector)

Return the needed coefficients type (real or complex) of polynomial with given zeros,
and for each zero an integer indication with the following meaning:

```
== 0:  real root
 < 0:  complex root without corresponding conjugate root of same multiplicity
 > 0:  complex root with corresponding conjugate root at position given by indication
```
 """
function zerostructure(zz::PolyZeros{<:Number})
    zerostructure(zz.z, zz.mult)
end
function zerostructure(z::AbstractVector{S}, ll::AbstractVector{<:Integer}) where S<:Number
    n = length(z)
    res = zeros(Int, n)
    isreal(z) && return real(S), res
    for k = 1:n
        res[k] != 0 && continue
        zk = z[k]
        isreal(zk) && continue
        res[k] = -1
        lk = ll[k]
        for j = k+1:n
            zj = z[j]
            lj = ll[j]
            if lj == lk && zk == conj(zj)
                res[k] = j
                res[j] = k
                break
            end
        end
    end
    T = all(res .>= 0) ? real(S) : S
    return T, res
end

"""
    `filter_duplicates(roots, multiplicities)`

Remove roots with multiplicities <= 0 and join duplicates
(add multiplicities for equal root values).
"""
function filter_duplicates(z::AbstractVector{<:Number}, ll::AbstractVector{<:Integer})
    n = length(z)
    if length(ll) == 0
        ll = ones(Int, n)
    else
        length(ll) == n ||
            throw(ArgumentError("roots and multiplicities have differing sizes"))
        all( ll .> 0 ) || throw(ArgumentError("only positive multiplicities allowed"))
    end
    lorig = ll
    for k = 1:n
        zk = z[k]
        for j = k+1:n
            zj = z[j]
            lj = ll[j]
            if zk == z[j] && lj > 0
                if ll === lorig
                    ll = copy(ll)
                end
                ll[k] += lj
                ll[j] = 0
            end
        end
    end
    if ll !== lorig
        jj = findall(iszero, ll)
        z = deleteat!(copy(z), jj)
        deleteat!(ll, jj)
    end
    z, ll
end

#
# calculate product over j of (x - z[j])^ll[j]
#
function evaluate(z::PolyZeros{S}, x::T, ll::AbstractVector{Int}) where {S,T}

  l = ll
	s = one(promote_type(T,S))
	n = length(ll)
	i = 1
	zz = z.z
	while i <= n
	    zi = x - zz[i]
		li = l[i]
	    if i < n && zz[i+1] == conj(zz[i])
		    m = min(li, l[i+1])
			if m > 0
		        s = s * (real(zi)^2 + imag(zi)^2) ^ m
			end
			if li > m
			    s = s * zi ^ (li - m)
			else
			    s = s * conj(zi) ^ (l[i+1] - m)
			end
			i += 2
		else
	        s = s * zi ^ li
			i += 1
		end
    end
    isreal(s) ? real(s) : s
end
evaluate(z::PolyZeros{S}, x::T) where {S,T} = evaluate(z, x, z.mult)

"""
"""
function Base.sort(z::PolyZeros)
    m = length(z.z)
    flist = [
             a::Int->z.mult[a],
             a::Int->abs(z.z[a]),
             a::Int->real(z.z[a]),
             a::Int->imag(z.z[a])
            ]
    wlist = [
             1000,
             1.0,
             0.1,
             0.01
            ]
    ev(f, a) = f(a)
    function lt(a::Int, b::Int)
        dot(ev.(flist, a), wlist) < dot(ev.(flist, b), wlist)
    end
    jj = sortperm(1:m, lt=lt)
    PolyZeros(z.z[jj], z.mult[jj])
end

function diffmult(za::PolyZeros, zb::PolyZeros)
    maximum(abs, za.mult - zb.mult)
end
function diffzeros(za::PolyZeros, zb::PolyZeros)
    norm(za.z - zb.z, Inf)
end

"""
    subsetapprox(p::Vector{S}, q::Vector{S}, dist::Function)

Let `length(p) <= length(q)` and `S` have a distance function `dist(::S,::S)`.
Find a vector `x` of indices into `q` with `length(jj) == length(p)`
which minimizes the term `maximum(dist(p[i], q[x[i]) for i in 1:length(p))`.
"""
function subsetapprox(a::AbstractVector{T}, b::AbstractVector{T}, dist::Function) where T
    A = dist.(b, transpose(a))
    subsetapprox(A)
end

function subsetapprox(A::AbstractMatrix{S}) where S<:Real
    fi = zeros(Int, size(A, 2))
    branch_and_bound(A, fi, typemax(S))
end

function branch_and_bound2(A, fi, inbound)
    println("called b&b($fi <= $inbound)")
    res = branch_and_bound1(A, fi, inbound)
    v, x = res
    println("returning $v")
    if minimum(x) < 1 || maximum(x) > size(A,1)
        println("wrong x after b&b:")
        display(x')
    end
    if false
        vv, xx = result(A, x)
        display([x xx]')
    end
    res
end

function result(A, y)
    z = getindex.(Ref(A), y, axes(y)[1])
    target(z), z
end

function target(z)
    norm(z, 2)
end

function fix!(fj, fi, j, i)
    if !iszero(fi[i])
        fj[fi[i]] = 0
    end
    fi[i] = j
    fj[j] = i
end

function checkfi(f)
    g = f[f .> 0]
    length(g) == length(unique(g))
end

function copyfi(fi, n)
    fi = copy(fi)
    fj = zeros(eltype(fi), n)
    for k = 1:length(fi)
        if fi[k] != 0
            fj[fi[k]] = k
        end
    end
    fi, fj
end

function branch_and_bound(A::Matrix{S}, fi, inbound::S) where S<:Real
    isopt, bound, y, ix = heuristics(A, fi, inbound)
    isopt == 1 && return bound, y
    isopt == 2 && return Inf, y
    n, m = size(A)
    fi, fj = copyfi(fi, n)
    fix!(fj, fi, y[ix], ix)
    ci, xi = branch_and_bound(A, fi, bound)
    if ci >= bound
        ci, xi = bound, y
    end
    for j = sortperm(view(A, :, ix))
        iszero(fj[j]) || continue
        cij = A[j,ix]
        cij > ci && continue
        fix!(fj, fi, j, ix)
        bij, xx = branch_and_bound(A, fi, ci)
        isinf(bij) && continue
        if bij < ci
            ci = bij
            xi = xx
        end
    end
    ci, xi
end

function heuristics(A::AbstractMatrix{S}, fi::V, inbounds::S) where {S<:Real,V<:AbstractVector{Int}}

    n, m = size(A)
    m <= n || throw(ArgumentError("first vector longer than second vector"))
    function findmin2(i, fj)
        cm, jm = typemax(S), 0
        for j = axes(fj)[1]
            if iszero(fj[j]) && A[j,i] < cm
                cm = A[j,i]
                jm = j
            end
        end
        cm, jm
    end
    findmin2(i, fj, fi) = fi[i] == 0 ? findmin2(i, fj) : (A[fi[i],i], fi[i])
    function findmax2(c, fi)
        cm, ii = zero(S), 0
        for i = axes(fi)[1]
            if iszero(fi[i]) && c[i] > cm
                cm = c[i]
                ii = i
            end
        end
        cm, ii
    end

    fi, fj = copyfi(fi, n)
    c = similar(fi, S)
    y = similar(fi)

    for i = 1:m
        c[i], y[i] = findmin2(i, fj, fi)
    end
    ci = target(c)
    if ci > inbounds # the most optimistic estimation can't be better than inbounds
        return 2, ci, y, 0
    end
    ix = 0
    bb = zero(S)
    isopt = 1
    free = count(iszero, fi)
    while free > 0
        ci, ic = findmax2(c, fi)
        jc = y[ic]
        if ci > bb
            bb = ci
            ix = ic
        end
        fix!(fj, fi, jc, ic)
        free -= 1
        for i in 1:m
            if fi[i] == 0 && y[i] == jc
                c[i], y[i] = findmin2(i, fj)
                if c[i] > bb
                    bb = ci
                    ix = i
                end
                isopt = 0
            end
        end
    end
    cc = target(c)
    isopt, cc, y, ix
end

# brute-force calculation of minimum
function bf(A)
    c = typemax(eltype(A))
    x = Int[]
    for jj in allcomb(size(A,2), size(A,1))
        cjj = target(A[jj[i], i] for i in axes(A,2))
        if cjj < c
            c = cjj
            x = jj
        end
    end
    c, x
end
function allcomb(a::Integer, b::Integer)
    allcomb.(b-a+1, b, 0:binomial(b, a) * factorial(a)-1)
end
function allcomb(m::Integer, n::Integer, z::Integer)
    function op((a, b), k)
        x, y = divrem(b, k)
        [a; y], x
    end
    function op2((a,b), k)
        j = 0
        for l = 1:length(b)
            if b[l] == 0
                if j == k
                    b[l] = 1
                    return [a; l], b
                else
                    j += 1
                end
            end
        end
        (a, b)
    end
    foldl(op2, foldl(op, n:-1:m, init = (Int[], z))[1], init = (Int[], zeros(n)))[1]
end


text:
# Hello world example, similar to the Boost.Python hello world

using Base.Test
using Compat
using CxxWrap

# Wrap the functions defined in C++
wrap_modules(CxxWrap._l_types)

using CppTypes
using CppTypes.World

# Stress test
for i in 1:1000000
  d = CppTypes.DoubleData()
end

# Default constructor
@test World <: CxxWrap.CppAny
@test supertype(World) == CxxWrap.CppAny
w = World()
println("Dumping type w...")
dump(w)
@test CppTypes.greet(w) == "default hello"

@show fw = CppTypes.world_factory()
@test CppTypes.greet(fw) == "factory hello"

swf = CppTypes.shared_world_factory()
@test CppTypes.greet_shared(swf) == "shared factory hello"
@test CppTypes.greet_shared_const(swf) == "shared factory hello"
@test CppTypes.greet(swf[]) == "shared factory hello" # Explicit dereference
@test CppTypes.greet(swf) == "shared factory hello" # Automatic conversion
swf2 = CppTypes.smart_world_factory()
@test CppTypes.greet_smart(swf2) == "smart factory hello"
@test CppTypes.greet_smart(swf) == "shared factory hello" # auto-convert between pointers
@test CppTypes.greet(swf2[]) == "smart factory hello" # Explicit dereference
@test CppTypes.greet(swf2) == "smart factory hello" # Automatic conversion
@test CppTypes.greet_weak(swf) == "shared factory hello"
@test_throws ErrorException CppTypes.greet_weak(swf2) == "shared factory hello"

swfr = CppTypes.shared_world_ref()
@test CppTypes.greet(swfr) == "shared factory hello ref"
CppTypes.reset_shared_world!(swfr, "reset shared pointer")
@test CppTypes.greet(swfr) == "reset shared pointer"

@test CppTypes.greet(CppTypes.boxed_world_factory()) == "boxed world"
@test CppTypes.greet(CppTypes.boxed_world_pointer_factory()) == "boxed world pointer"
@test CppTypes.greet(CppTypes.world_ref_factory()) == "reffed world"

@show uwf = CppTypes.unique_world_factory()
@test CppTypes.greet(uwf) == "unique factory hello"

byval = CppTypes.world_by_value()
@test CppTypes.greet(byval) == "world by value hello"

CppTypes.set(w, "hello")
@show CppTypes.greet(w)
@test CppTypes.greet(w) == "hello"
@test CppTypes.greet_lambda(w) == "hello"

w = World("constructed")
@test CppTypes.greet(w) == "constructed"

w_assigned = w
w_deep = deepcopy(w)

@test w_assigned == w
@test w_deep != w

# Destroy w: w and w_assigned should be dead, w_deep alive
finalize(w)
println("finalized w")
if !(is_windows() && Sys.WORD_SIZE == 32)
  @test_throws ErrorException CppTypes.greet(w)
  println("throw test 1 passed")
  @test_throws ErrorException CppTypes.greet(w_assigned)
  println("throw test 2 passed")
end
@test CppTypes.greet(w_deep) == "constructed"
println("completed deepcopy test")

wnum = World(1)
@test CppTypes.greet(wnum) == "NumberedWorld"
finalize(wnum)
@test CppTypes.greet(wnum) == "NumberedWorld"

noncopyable = CppTypes.NonCopyable()
other_noncopyable = deepcopy(noncopyable)
@test other_noncopyable.cpp_object == noncopyable.cpp_object

@test value(value(ReturnConstRef())) == 42

wptr = World()
@test CppTypes.greet(ConstPtrConstruct(wptr)) == "default hello"

call_op = CallOperator()
@test call_op() == 43
@test call_op(42) == 42

type JuliaTestType
    a::Float64
    b::Float64
end

function julia_test_func(data)
  println("a: ", data.a, ", b: ", data.b)
  @test data.a == 2.
  @test data.b == 3.
end

CppTypes.call_testype_function()

@test enum_to_int(CppTypes.EnumValA) == 0
@test enum_to_int(CppTypes.EnumValB) == 1
@test get_enum_b() == CppTypes.EnumValB
@test CppTypes.EnumValA + CppTypes.EnumValB == CppTypes.EnumValB
@test CppTypes.EnumValA | CppTypes.EnumValB == CppTypes.EnumValB

using CppTypes: Foo, name, data, print_foo_array

foovec = Any[Foo("a", [1.0, 2.0, 3.0]), Foo("b", [11.0, 12.0, 13.0])] # Must be Any because of the boxing
@show name(foovec[1])
@show data(foovec[1])
print_foo_array(foovec)


text:
using ResumableFunctions, SimJulia, Distributions, BenchmarkTools

@resumable function exp_source(sim::Simulation, lambd::Float64, server::Resource, mu::Float64)
  while true
    dt = rand(Exponential(1 / lambd))
    @yield timeout(sim, dt)
    @process customer(sim, server, mu)
  end
end

@resumable function customer(sim::Simulation, server::Resource, mu::Float64)
  @yield request(server)
  dt = rand(Exponential(1 / mu))
  @yield timeout(sim, dt)
  @yield release(server)
end

function test_mm1(n::Float64)
  sim = Simulation()
  server = Resource(sim)
  @process exp_source(sim, 1.0, server, 1.1)
  run(sim, n)
end

test_mm1(100.0)
@btime test_mm1(100.0)


text:
#=
  Euler #19 in Julia.

  Problem 19:
  """
  You are given the following information, but you may prefer
  to do some research for yourself.

  * 1 Jan 1900 was a Monday.
  * Thirty days has September,
    April, June and November.
    All the rest have thirty-one,
    Saving February alone,
    Which has twenty-eight, rain or shine.
    And on leap years, twenty-nine.
  * A leap year occurs on any year evenly divisible by 4, but not
    on a century unless it is divisible by 400.

  How many Sundays fell on the first of the month during the
  twentieth century (1 Jan 1901 to 31 Dec 2000)?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

using Dates

# 0.00001321s
function euler19a()
    count = 0;
    for year in 1901:2000
        for month in 1:12
            if Dates.dayofweek(Date(year, month, 1)) == 7
                count += 1
            end
        end
    end

    return count
end

run_euler(euler19a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
#
# Kelley's method with lazy hyperplanes.
#
"""
    KelleyMethod

The parameters for solving a Lagrangian dual with Kelley's method.

# Arguments
* initialbound   Initial upper or lower bound for the Lagrangian dual
* tol            Tolerance: we stop when the gap between our approximation of the function and the actual function value is less than `tol`
* maxit          To terminate the method
"""
mutable struct KelleyMethod{T<:Tolerance} <: AbstractLagrangianMethod
    initialbound::Float64   # starting bound for the Lagrangian dual problem
    tol::T                  # tolerance for terminating
    maxit::Int              # a cap on iterations
end
function KelleyMethod(; initialbound=0.0, tol=Unit(1e-6), maxit=10_000)
    KelleyMethod(initialbound, tol, maxit)
end

"""
    lagrangian_method!(lp::LinearProgramData{KelleyMethod}, m::JuMP.Model, π::Vector{Float64})

Kelley's method with a lazy callback approach.

# Arguments
* lp        Information about the primal problem
* m         The primal problem
* π         Initial iterate

# Returns
* status, objective, and modifies π
"""
function lagrangian_method!(lp::LinearProgramData{KelleyMethod{T}}, m::JuMP.Model, π::Vector{Float64}) where T

    kelleys = lp.method
    N      = length(π)
    tol    = kelleys.tol

    # Gap between approximate model and true function each iteration
    gap = Inf
    # Let's make new storage for the best multiplier found so far
    bestmult = copy(π)
    # Dual problem has the opposite sense to the primal
    dualsense = getdualsense(m)
    fdash = zeros(N)

    # The approximate model will be a made from linear hyperplanes
    approx_model = Model(solver=m.solver)

    @variables approx_model begin
        θ                   # The objective of the approximate model
        x[i=1:N]            # The x variables in the approximate model are the Lagrangian duals.
    end
    # There are sign restrictions on some duals
    for (i, sense) in enumerate(lp.senses)
        if sense == :ge
            setupperbound(x[i], 0.0)
        elseif sense == :le
            setlowerbound(x[i], 0.0)
        end
    end
    # Let's not be unbounded from the beginning
    if dualsense == :Min
        setlowerbound(θ, kelleys.initialbound)
        best_actual, f_actual, f_approx = Inf, Inf, -Inf
    else
        setupperbound(θ, kelleys.initialbound)
        best_actual, f_actual, f_approx = -Inf, -Inf, Inf
    end

    iteration = 0

    while iteration < kelleys.maxit
        iteration += 1
        # Evaluate the real function and a subgradient
        f_actual, fdash = solve_primal(m, lp, π)

        # Improve the model and update best function value so far
        if dualsense == :Min
            @constraint(approx_model, θ >= f_actual + dot(fdash, x - π))
            if f_actual < best_actual
                best_actual = f_actual
                bestmult .= π
            end
        else
            @constraint(approx_model, θ <= f_actual + dot(fdash, x - π))
            if f_actual > best_actual
                best_actual = f_actual
                bestmult .= π
            end
        end
        # Get a bound from the approximate model
        @objective(approx_model, dualsense, θ)
        @assert solve(approx_model) == :Optimal
        f_approx = getobjectivevalue(approx_model)::Float64
        # Check the gap
        gap = abs(best_actual - f_approx)
        #=
            Stop if best_actual ≈ f_approx
            Note: if fdash ≈ ̃0, then we expect that best_actual ≈ f_approx
            anyway. We still check for this condition in case the solver allows
            a component of ̃x to get very large, so that dot(fdash, x) is nonzero
            and f_approx is incorrect.
        =#
        if closetozero(gap, best_actual, f_approx, tol) || isclose(norm(fdash), 0.0, tol)
            π .= bestmult
            if dualsense == :Min
                π .*= -1 # bestmult not the same as getvalue(x), approx_model may have just gotten lucky
            end
            return :Optimal, best_actual::Float64
        end

        # Get the next iterate
        π .= getvalue(x)
    end
    warn("Lagrangian relaxation did not solve properly.")
    return :IterationLimit, f_approx::Float64

end