text:
using Polynomials, Unitful, UnitfulAstro

# Convenience function for wavelength conversion
aa_to_invum(wave::Real) = 10000 / wave
aa_to_invum(wave::Quantity) = aa_to_invum(ustrip(u"angstrom", wave))

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

# Optical coefficients
const ccm89_ca = Poly([1.0, 0.17699, -0.50447, -0.02427, 0.72085, 0.01979, -0.7753, 0.32999])
const ccm89_cb = Poly([0.0, 1.41338, 2.28305, 1.07233, -5.38434, -0.62251, 5.3026, -2.09002])

const od94_ca = Poly([1.0, 0.104, -0.609, 0.701, 1.137, -1.718, -0.827, 1.647, -0.505])
const od94_cb = Poly([0.0, 1.952, 2.908, -3.989, -7.985, 11.102, 5.491, -10.805, 3.347])

"""
    ccm89(λ::Real, Rv=3.1)
    ccm89(λ::Quantity, Rv=3.1)

Clayton, Cardelli and Mathis (1989) dust law.

Returns E(B-V) in magnitudes at the given wavelength relative to the extinction
at 5494.5 Å. `λ` is the wavelength in Å and has support over [1000, 33333].
Outside of that range this will return 0. `Rv` is the selective extinction
and is valid over [2, 6]. A typical value for the Milky Way is 3.1

If `λ` is a `Unitful.Quantity` it will be automatically converted to Å and the
returned value will be `UnitfulAstro.mag`.

# References
[Clayton,Cardelli and Mathis (1989)](https://ui.adsabs.harvard.edu/abs/1989ApJ...345..245C)
"""
function ccm89(λ::Real, Rv = 3.1)
    x = aa_to_invum(λ)
    return ccm89_invum(x, Rv, ccm89_ca, ccm89_cb)
end

ccm89(λ::Quantity, Rv::Real = 3.1) = ccm89(ustrip(u"angstrom", λ), Rv) * u"mag"

"""
    od94(λ::Real, Rv=3.1)
    od94(λ::Quantity, Rv=3.1)

O'Donnell (1994) dust law.

This is identical to the Clayton, Cardelli and Mathis (1989) dust law, except
for different coefficients used in the optical (3030.3 Å to 9090.9 Å).

If `λ` is a `Unitful.Quantity` it will be automatically converted to Å and
the returned value will be `UnitfulAstro.mag`.

# References
[O'Donnell (1994)](https://ui.adsabs.harvard.edu/abs/1994ApJ...422..158O)

# See Also
[`ccm89`](@ref)
"""
function od94(λ::Real, Rv = 3.1)
    x = aa_to_invum(λ)
    return ccm89_invum(x, Rv, od94_ca, od94_cb)
end

od94(λ::Quantity, Rv = 3.1) = ccm89(ustrip(u"angstrom", λ), Rv) * u"mag"

function ccm89_invum(x::Real, Rv::Real, c_a::Poly{<:Real}, c_b::Poly{<:Real})
    if x < 0.3
        return 0.0x
    elseif x < 1.1  # Near IR
        y = x^1.61
        a = 0.574y
        b = -0.527y
    elseif x < 3.3  # Optical
        y = x - 1.82
        yn = 1.0
        a = c_a(y)
        b = c_b(y)
    elseif x < 8.0  # NUV
        a =  1.752 - 0.316x - (0.104 / ((x - 4.67)^2 + 0.341))
        b = -3.090 + 1.825x + (1.206 / ((x - 4.62)^2 + 0.263))
        if x > 5.9 # Far NUV
            y = x - 5.9
            a += Poly([0.0, 0.0, -0.04473, -0.009779])(y)
            b += Poly([0.0, 0.0, 0.213, 0.1207])(y)
        end
    elseif x ≤ 10.0 # FUV
        y = x - 8.0
        a = Poly([-1.073, -0.628, 0.137, -0.07])(y)
        b = Poly([13.67, 4.257, -0.42, 0.374])(y)
    else
        return 0.0x
    end

    return a + b / Rv
end

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

"""
    cal00(λ::Real, Rv=4.05)
    cal00(λ::Quantity, Rv=4.05)

Calzetti et al. (2000) Dust Law.

Returns E(B-V) in magnitudes at the given wavelength. `λ` is the wavelength in Å
 and has support over [1200, 22000]. Outside of that range this will return 0.

Calzetti et al. (2000) developed a recipe for dereddening the spectra of
galaxies where massive stars dominate the radiation output. They found the best
 fit value for such galaxies was 4.05±0.80.

If `λ` is a `Unitful.Quantity` it will be automatically converted to Å and the
returned value will be `UnitfulAstro.mag`.

# References
[Calzetti et al. (2000)](https://ui.adsabs.harvard.edu/abs/2000ApJ...533..682C)
"""
function cal00(λ::Real, Rv = 3.1)
    # Convert to inverse-um
    x = aa_to_invum.(λ)
    return cal00_invum(x, Rv)
end

cal00(λ::Quantity, Rv::Real = 3.1) = cal00(ustrip(u"angstrom", λ), Rv) * u"mag"

function cal00_invum(x::Real, Rv::Real)

    if x > 1 / 0.12
        return 0.0x
    elseif x > 1 / 0.63
        k = 2.659 * (((0.011x - 0.198)x + 1.509)x - 2.156)
    elseif x >  1 / 2.2
        k = 2.659 * (1.040x - 1.857)
    else
        return 0.0x
    end

    return 1.0 + k / Rv

end


text:
function read_IMA(dt::Date)

    @assert BusinessDays.isbday(:BRSettlement, dt) "$dt should be a business day in Brazil"

    url = "https://www.anbima.com.br/informacoes/ima/ima-sh-down.asp"
    str_dt = Dates.format(dt, "dd%2Fmm%2Fyyyy")

    header = [
        # POST /informacoes/ima/ima-sh-down.asp HTTP/1.1
        "Host" =>  "www.anbima.com.br"
        "User-Agent" =>  "Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:80.0) Gecko/20100101 Firefox/80.0"
        "Accept" =>  "text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8"
        "Accept-Language" =>  "pt-BR,pt;q=0.8,en-US;q=0.5,en;q=0.3"
        "Accept-Encoding" =>  "gzip, deflate, br"
        "Content-Type" =>  "application/x-www-form-urlencoded"
        "Content-Length" =>  "94"
        "Origin" =>  "https://www.anbima.com.br"
        "Connection" =>  "keep-alive"
        "Referer" =>  "https://www.anbima.com.br/informacoes/ima/ima-sh.asp"
        "Cookie" =>  "lumClientId=8A2AB290749450B201749709905D1504; JSESSIONID=4B569D6B9E949CB734EBD02313131CD1.LumisProdA; lumUserLocale=pt_BR; lumMonUid=GB9v_UzD1CZbwjVzKKYvn0QTc3CJxsuf; AWSELB=359F6F8906C702E7D2BBD3C04A400F8539C6D91FD238C2952FBF4A0F2BB3C549A6A51A2676854C5B9E1FDC2EB533CB8BAA69DCEB7EAD83CF8BAD191790A634B39952D64C276ECAB45E8107BE5B4A03390F9EA29851; ASPSESSIONIDQESTDBAS=BECPHPJAMIHJIAGBHAMGCGCB; BIGipServerPool_ANBSPCLD-WEB02_SSL=204005386.47873.0000; __utma=234609614.893389923.1600261823.1600261823.1600275110.2; __utmc=234609614; __utmz=234609614.1600261823.1.1.utmcsr=google|utmccn=(organic)|utmcmd=organic|utmctr=(not%20provided); _ga=GA1.3.893389923.1600261823; _gid=GA1.3.220011512.1600261823; _fbp=fb.2.1600261824060.1044754934; _hjTLDTest=1; _hjid=67d810ba-8ef2-40f5-80d2-63d7dbac00b7; __trf.src=encoded_eyJmaXJzdF9zZXNzaW9uIjp7InZhbHVlIjoiMjM0NjA5NjE0LjE2MDAyNjE4MjMuMS4xLnV0bWNzcj1nb29nbGV8dXRtY2NuPShvcmdhbmljKXx1dG1jbWQ9b3JnYW5pY3x1dG1jdHI9KG5vdCBwcm92aWRlZCkiLCJleHRyYV9wYXJhbXMiOnt9fSwiY3VycmVudF9zZXNzaW9uIjp7InZhbHVlIjoiaHR0cHM6Ly93d3cuZ29vZ2xlLmNvbS8iLCJleHRyYV9wYXJhbXMiOnt9fSwiY3JlYXRlZF9hdCI6MTYwMDI2NDEzMjA4OX0=; rdtrk=%7B%22id%22%3A%221d5766e6-4883-442c-930a-d1d22efb296f%22%7D; ADRUM=s=1600275093544&r=https%3A%2F%2Fwww.anbima.com.br%2Fpt_br%2Finformar%2Fprecos-e-indices%2Findices%2Fima.htm%3F0; lumUserSessionId=7yqmdlOCqWvJfHnbGi6641Sq1QT_xFq2; lumUserName=Guest; lumIsLoggedUser=false"
        "Upgrade-Insecure-Requests" =>  "1"
    ]

    # why u not work??
    # body = HTTP.Form(Dict([
    #   "Tipo" => ""
    #   "DataRef" => ""
    #   "Pai" => "ima"
    #   "escolha" => "2"
    #   "Idioma" => "EN"
    #   "saida" => "csv"
    #   "Dt_Ref_Ver" => "20200908"
    #   "Dt_Ref" => "03%2F08%2F2020"
    # ]))

    body = "Tipo=&DataRef=&Pai=ima&escolha=2&Idioma=PT&saida=csv&Dt_Ref_Ver=20200908&Dt_Ref=$str_dt"

    r = HTTP.post(url, header, body)

    response_body = String(r.body)
    out = _parse_ima_result(response_body)
    for (k, ima) in out.elements
        @assert ima.date == dt "Invalid date for $(k): expected $(dt) got $(ima.date)"
    end
    return out

end

function _parse_ima_result(s::String)
    lines = split(s, "\r\n")
    if length(lines) <= 2
        error("nao ha dados de quadro resumo")
    end
    v_elements = Vector{IMA}()
    for line in lines[3:(end-1)]
        vals = split(line, ";")
        @assert length(vals) == 18 "esperava 18 campos na linha mas recebeu $(length(vals)): $line"
        element = IMA(
            vals[1],
            Date(vals[2], "dd/mm/yyyy"),
            parse_value(Float64, vals[3]),
            parse_value(Float64, vals[4]),
            parse_value(Float64, vals[5]),
            parse_value(Float64, vals[6]),
            parse_value(Float64, vals[7]),
            parse_value(Float64, vals[8]),
            parse_value(Float64, vals[9]),
            parse_value(Float64, vals[10]),
            parse_value(Int64, vals[11]),
            parse_value(Float64, vals[12]),
            parse_value(Float64, vals[13]),
            parse_value(Float64, vals[14]),
            parse_value(Float64, vals[15]),
            parse_value(Float64, vals[16]),
            parse_value(Float64, vals[17]),
            parse_value(Float64, vals[18])
        )
        push!(v_elements, element)
    end
    return IMAResult(v_elements)
end

ext: jl


text:
function update!(s::Simulation{T}, event::AbstractEvent) where T <: EpidemicModel
  if _time(event) < s.time
    @error "Time of a new event must be >= that of the previous event"
  elseif _time(event) < Inf
    update!(s.events, event)
    update!(s.disease_states, event)
    update!(s.transmission_network,
            generate(Transmission,
                     s.transmission_rates,
                     event))
    update!(s.transmission_rates,
            s.event_rates,
            event,
            s.disease_states,
            s.population,
            s.risk_functions,
            s.risk_parameters)
    # Count the iteration as having occurred.
    s.iterations += 1
  end
  # Update simulation time
  s.time = _time(event)
  # Return updated Simulation object
  return s
end


text:
import Base.isempty

export AbstractCentrallySymmetric,
       center,
       an_element

"""
    AbstractCentrallySymmetric{N<:Real} <: LazySet{N}

Abstract type for centrally symmetric sets.

### Notes

Every concrete `AbstractCentrallySymmetric` must define the following functions:

- `center(::AbstractCentrallySymmetric{N})::Vector{N}` -- return the center
    point

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractCentrallySymmetric)
3-element Array{Any,1}:
 Ball2
 Ballp
 Ellipsoid
```
"""
abstract type AbstractCentrallySymmetric{N<:Real} <: LazySet{N} end

"""
    dim(S::AbstractCentrallySymmetric)::Int

Return the ambient dimension of a centrally symmetric set.

### Input

- `S` -- set

### Output

The ambient dimension of the set.
"""
@inline function dim(S::AbstractCentrallySymmetric)::Int
    return length(center(S))
end

"""
    isbounded(S::AbstractCentrallySymmetric)::Bool

Determine whether a centrally symmetric set is bounded.

### Input

- `S` -- centrally symmetric set

### Output

`true` (since a set with a unique center must be bounded).
"""
function isbounded(::AbstractCentrallySymmetric)::Bool
    return true
end

"""
    an_element(S::AbstractCentrallySymmetric{N})::Vector{N} where {N<:Real}

Return some element of a centrally symmetric set.

### Input

- `S` -- centrally symmetric set

### Output

The center of the centrally symmetric set.
"""
function an_element(S::AbstractCentrallySymmetric{N})::Vector{N} where {N<:Real}
    return center(S)
end

"""
    isempty(S::AbstractCentrallySymmetric)::Bool

Return if a centrally symmetric set is empty or not.

### Input

- `S` -- centrally symmetric set

### Output

`false`.
"""
function isempty(::AbstractCentrallySymmetric)::Bool
    return false
end


text:
@testset "NEQ Kirschmer" begin
  R,x = FlintZZ["x"]
  K,a = number_field(x^2+3x+1,"a")
  kt,t = K["t"]
  E, b = number_field( t^2 + (a+5)* (-2*a+2) )
  @test Hecke.is_norm( E, -a*(-2*a+2) )[1]

  K, a = NumberField(x^2 - 2, "a")
  Kt, t = K["t"]
  L, b = number_field(t^2 - 7)
  @test is_norm(L, K(324))[1]

  K,a = number_field(x^2+3*x+1)
  kt,t=K["t"]
  E,b=number_field(t^2-(a-1),"b")
  @test is_norm(E, a+3)[1]

  K, a = quadratic_field(-5)
  kt,t=K["t"]
  theta = 11*(4*a-22)
  E,b=number_field(t^2+theta,"b")

  @test is_norm(E, K(-11))[1]
end


text:
struct RuntimeDevice{D}
    device::D
end
default_device() = RuntimeDevice(default_device(RUNTIME[]))
default_device(::typeof(HSA_rt)) = HSARuntime.get_default_agent()

struct RuntimeQueue{Q}
    queue::Q
end
default_queue(device) = RuntimeQueue(default_queue(RUNTIME[], device))
default_queue(::typeof(HSA_rt), device) =
    HSARuntime.get_default_queue(device.device)
get_device(queue::RuntimeQueue{HSAQueue}) = RuntimeDevice(queue.queue.agent)

default_isa(device::RuntimeDevice{HSAAgent}) =
    HSARuntime.get_first_isa(device.device)

struct RuntimeEvent{E}
    event::E
end
create_event() = RuntimeEvent(create_event(RUNTIME[]))
create_event(::typeof(HSA_rt)) = HSASignal()
Base.wait(event::RuntimeEvent; kwargs...) = wait(event.event; kwargs...)

struct RuntimeExecutable{E}
    exe::E
end
create_executable(device, entry, obj; globals=()) =
    RuntimeExecutable(create_executable(RUNTIME[], device, entry, obj; globals=globals))
function create_executable(::typeof(HSA_rt), device, entry, obj; globals=())
    # link with ld.lld
    ld_path = HSARuntime.ld_lld_path
    @assert ld_path != "" "ld.lld was not found; cannot link kernel"
    path_exe = mktemp() do path_o, io_o
        write(io_o, obj)
        flush(io_o)
        path_exe = path_o*".exe"
        run(`$ld_path -shared -o $path_exe $path_o`)
        path_exe
    end
    data = read(path_exe)
    rm(path_exe)

    return HSAExecutable(device.device, data, entry; globals=globals)
end
HSARuntime.get_global(exe::RuntimeExecutable, sym::Symbol) =
    HSARuntime.get_global(exe.exe, sym)

struct RuntimeKernel{K}
    kernel::K
end
create_kernel(device, exe, entry, args) =
    RuntimeKernel(create_kernel(RUNTIME[], device, exe, entry, args))
create_kernel(::typeof(HSA_rt), device, exe, entry, args) =
    HSAKernelInstance(device.device, exe.exe, entry, args)
launch_kernel(queue, kern, event; kwargs...) =
    launch_kernel(RUNTIME[], queue, kern, event; kwargs...)
launch_kernel(::typeof(HSA_rt), queue, kern, event;
              groupsize=nothing, gridsize=nothing) =
    HSARuntime.launch!(queue.queue, kern.kernel, event.event;
                       workgroup_size=groupsize, grid_size=gridsize)


text:
@views function fit(estimator::LinearModelEstimators,
					X::AbstractMatrix,
					y::AbstractVector,
					z::AbstractMatrix,
					Z::AbstractMatrix,
					wts::AbstractVector)
	if !isa(wts, Weights)
        wts = FrequencyWeights(wts)
    else
        isa(wts, FrequencyWeights) || throw(ArgumentError("Only frequency weights are supported."))
    end
    X = transform(estimator, X, wts)
    y = transform(estimator, y, wts)
    z = transform(estimator, z, wts)
    Z = transform(estimator, Z, wts)
    wts = transform(estimator, FrequencyWeights(wts))
    if !isempty(z)
        Z̃ = hcat(X, Z)
        F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true, check = false)
        bkr = count(x -> abs(x) ≥ √eps(), F.D)
		if bkr < size(F, 2)
            lin_ind = sort!(invperm(F.p)[1:bkr])
            count(x -> lin_ind > size(X, 2) + 1) ≥ size(z, 2) ||
                throw(ArgumentError("Insufficient number of instruments."))
            Z̃ = Z̃[:,lin_ind]
            F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true)
        else
            lin_ind = collect(1:size(F, 2))
        end
        γ = F \ (Z̃' * Diagonal(wts) * z)
        X̃ = hcat(X, Z̃ * γ)
    else
        X̃ = X
    end
    F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true, check = false)
    bkr = count(x -> abs(x) ≥ √eps(), F.D)
    if bkr < size(F, 2)
        lin_ind = sort!(invperm(F.p)[1:bkr])
        X̃ = convert(Matrix{Float64}, X̃[:,lin_ind])
        F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true)
    else
        lin_ind = collect(1:size(F, 2))
    end
    β = F \ (X̃' * Diagonal(wts) * y)
    Ψ = Hermitian(inv(F))
    ŷ = isempty(z) ? X̃ * β : hcat(X, z)[:,lin_ind] * β
	X̃, y, β, Ψ, ŷ, wts, lin_ind
end
@views function fit!(obj::EconometricModel{<:LinearModelEstimators})
	@unpack estimator, X, y, z, Z, wts = obj
    if !isa(wts, Weights)
        wts = FrequencyWeights(wts)
    else
        isa(wts, FrequencyWeights) || throw(ArgumentError("Only frequency weights are supported."))
    end
    X = transform(estimator, X, wts)
    y = transform(estimator, y, wts)
    z = transform(estimator, z, wts)
    Z = transform(estimator, Z, wts)
    wts = transform(estimator, FrequencyWeights(wts))
    if !isempty(z)
        Z̃ = hcat(X, Z)
        F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true, check = false)
        bkr = count(x -> abs(x) ≥ √eps(), F.D)
		if bkr < size(F, 2)
            lin_ind = sort!(invperm(F.p)[1:bkr])
            count(x -> lin_ind > size(X, 2) + 1) ≥ size(z, 2) ||
                throw(ArgumentError("Insufficient number of instruments."))
            Z̃ = Z̃[:,lin_ind]
            F = bunchkaufman(Hermitian(Z̃' * Diagonal(w) * Z̃), true)
        else
            lin_ind = collect(1:size(F, 2))
        end
        γ = F \ (Z̃' * Diagonal(wts) * z)
        X̃ = hcat(X, Z̃ * γ)
    else
        X̃ = X
    end
    F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true, check = false)
    bkr = count(x -> abs(x) ≥ √eps(), F.D)
    if bkr < size(F, 2)
        lin_ind = sort!(invperm(F.p)[1:bkr])
        X̃ = convert(Matrix{Float64}, X̃[:,lin_ind])
        F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true)
    else
        lin_ind = collect(1:size(F, 2))
    end
    β = F \ (X̃' * Diagonal(wts) * y)
    Ψ = Hermitian(inv(F))
    ŷ = isempty(z) ? X̃ * β : hcat(X, z)[:,lin_ind] * β
	X = X̃
	wts = FrequencyWeights(collect(wts))
	@pack! obj = X, y, β, Ψ, ŷ, wts
	obj.vars = (obj.vars[1], obj.vars[2][lin_ind])
	obj
end
@views function fit!(obj::EconometricModel{<:ContinuousResponse})
	@unpack estimator, X, y, z, Z, wts = obj
    if !isa(wts, Weights)
        wts = FrequencyWeights(wts)
    else
        isa(wts, FrequencyWeights) || throw(ArgumentError("Only frequency weights are supported."))
    end
    y₀ = y
    X = transform(estimator, X, wts)
    y = transform(estimator, y, wts)
    z = transform(estimator, z, wts)
    Z = transform(estimator, Z, wts)
    w = transform(estimator, FrequencyWeights(wts))
    if !isempty(z)
        Z̃ = hcat(X, Z)
        F = bunchkaufman(Hermitian(Z̃' * Diagonal(w) * Z̃), true, check = false)
        bkr = count(x -> abs(x) ≥ √eps(), F.D)
		if bkr < size(F, 2)
            lin_ind = sort!(invperm(F.p)[1:bkr])
            count(x -> lin_ind > size(X, 2) + 1) ≥ size(z, 2) ||
                throw(ArgumentError("Insufficient number of instruments."))
            Z̃ = Z̃[:,lin_ind]
            F = bunchkaufman(Hermitian(Z̃' * Diagonal(w) * Z̃), true)
        else
            lin_ind = collect(1:size(F, 2))
        end
        γ = F \ (Z̃' * Diagonal(w) * z)
        X̃ = hcat(X, Z̃ * γ)
    else
        X̃ = X
    end
    F = bunchkaufman(Hermitian(X̃' * Diagonal(w) * X̃), true, check = false)
    bkr = count(x -> abs(x) ≥ √eps(), F.D)
    if bkr < size(F, 2)
        lin_ind = sort!(invperm(F.p)[1:bkr])
        X̃ = convert(Matrix{Float64}, X̃[:,lin_ind])
        F = bunchkaufman(Hermitian(X̃' * Diagonal(w) * X̃), true)
    else
        lin_ind = collect(1:size(F, 2))
    end
    β = F \ (X̃' * Diagonal(w) * y)
    Ψ = Hermitian(inv(F))
    ŷ = isempty(z) ? X̃ * β : hcat(X, z)[:,lin_ind] * β
    û = y - ŷ
    ŷ .= y₀ .- û
	X = X̃
	y = y₀
	@pack! obj = X, y, β, Ψ, ŷ
	obj.vars = (obj.vars[1], obj.vars[2][lin_ind])
	obj
end
"""
    obtain_Ω(A::AbstractMatrix{<:Real},
             μ::AbstractMatrix{<:Real},
             ω::FrequencyWeights)
Obtain Ω for a multinomial regression by building the matrix by blocks.
"""
@views function obtain_Ω(A, μ, wts)
  Aᵀ = transpose(A)
  k = size(μ, 2)
  p = size(A, 2)
  Σ = zeros(p * (k - 1), p * (k - 1))
  for (idx₀, col) ∈ enumerate(1:p:(p * (k - 1)))
    for (idx₁, row) ∈ enumerate(col:p:(p * (k - 1)))
      rows = col:col + p - 1
      cols = row:row + p - 1
      if col == row
        Σ[rows, cols] = Aᵀ * Diagonal(μ[:,idx₀ + 1] .* (1 .- μ[:,idx₀ + 1]) .* wts) * A
      else
        Σ[rows, cols] = -(Aᵀ * Diagonal(μ[:,idx₀ + 1] .* μ[:,idx₁ + 1] .* wts) * A)
      end
    end
  end
  Hermitian(Σ)
end
@views function fit!(obj::EconometricModel{<:NominalResponse})
	@unpack estimator, X, y, z, Z, wts = obj
	@unpack categories = estimator
    @assert isempty(z) && isempty(Z) "Nominal response models can only contain exogenous features"
    b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
    F = qr(X, Val(true))
    qrr = count(x -> abs(x) ≥ √eps(), diag(F.R))
    if qrr < size(F, 2)
        lin_ind = sort!(invperm(F.p)[1:qrr])
        X = convert(Matrix{Float64}, X[:,lin_ind])
        F = qr(X, Val(true))
    else
        lin_ind = collect(1:size(F, 2))
    end
    Q = Matrix(F.Q)
    Q⊤ = transpose(Q)
    m, p = size(F)
    k = size(b, 2)
    η = zeros(m, k)
    μ, μ′, w = zero(η), zero(η), zero(η)
    β = zeros(p, k)
    ℓℓ = [Inf, Inf, Inf]
    max_iteration = 250
    tol = 1e-9
    converged = false
    for iteration ∈ 1:max_iteration
      converged && break
      for row ∈ 1:m
        μ[row,:] .= softmax(η[row,:])
      end
      μ′ .= max.(μ .* (one(Float64) .- μ), √eps())
      w .= wts .* μ′
      ℓℓ[1] = ℓℓ[2]
      ℓℓ[2] = sum(wts[idx[1]] * logpdf(Categorical(collect(μ[idx[1],:])), idx[2]) for idx in findall(b))
      if isone(iteration)
          ℓℓ₀ = ℓℓ[2]
      end
      η .+= (b .- μ) ./ μ′
      for idx ∈ 2:size(b, 2)
        C = cholesky!(Hermitian(Q⊤ * Diagonal(w[:,idx]) * Q)).factors
        β[:,idx] = LowerTriangular(transpose(C)) \ (Q⊤ * Diagonal(w[:,idx]) * η[:,idx])
        β[:,idx] = UpperTriangular(C) \ β[:,idx]
      end
      η .= Q * β
      ℓℓ[3] = abs(ℓℓ[2] - ℓℓ[1])
      converged = ℓℓ[3] < 1e-9 || iszero(ℓℓ[1])
    end
    converged || throw(ConvergenceException(max_iteration, ℓℓ[3], tol))
    β .= F \ η
    Ψ = Hermitian(inv(bunchkaufman!(obtain_Ω(X, μ, wts))))
    ŷ = X * β
    β = collect(vec(β)[size(X, 2) + 1:end])
	@pack! obj = X, y, β, Ψ, ŷ
	obj.vars = (obj.vars[1], obj.vars[2][lin_ind])
	obj
end
@views function fit!(obj::EconometricModel{<:OrdinalResponse})
	@unpack estimator, X, y, z, Z, wts = obj
	@unpack categories = estimator
    @assert isempty(z) && isempty(Z) "Ordinal response models can only contain exogenous features"
    F = qr(X, Val(true))
    qrr = count(x -> abs(x) ≥ √eps(), diag(F.R))
    if qrr < size(F, 2)
        lin_ind = sort!(invperm(F.p)[1:qrr])
        X = convert(Matrix{Float64}, X[:,lin_ind])
    else
        lin_ind = collect(1:size(F, 2))
    end
    m, n = size(X)
    xs = 1:n
    ks = length(categories) - 1
    ζs = n + 1:(n + ks)
    δ₀ = mapreduce(x -> (1:ks .== x - 1)', vcat, y)
    δ₁ = mapreduce(x -> (1:ks .== x)', vcat, y)
    η = Vector{Union{Float64,Dual}}(undef, m)
    ζ = vcat(-Inf, similar(η, ks), Inf)
    y₁ = y .+ 1
    Y₀ = similar(η)
    Y₁ = similar(η)
    F₀ = similar(η)
    F₁ = similar(η)
    f₀ = similar(η)
    f₁ = similar(η)
    pr = similar(η)
    θ₁ = ones(ks)
    function f(β)
      η .= X * β[xs]
      ζ[2:end - 1] .= cumsum(vcat(β[n + 1], exp.(β[n + 2:end])))
      Y₀ .= max.(-100, ζ[y] - η)
      Y₁ .= min.(100, ζ[y₁] - η)
      F₀ .= cdf.(Logistic(), Y₀)
      F₁ .= cdf.(Logistic(), Y₁)
      pr .= F₁ .- F₀
      -sum(w * log(p) for (w, p) ∈ zip(wts, pr))
    end
    function g!(G, β)
      ζ[2:end - 1] .= cumsum(vcat(β[n + 1], exp.(β[n + 2:end])))
      η .= X * β[xs]
      Y₀ .= max.(-100, ζ[y] - η)
      Y₁ .= min.(100, ζ[y₁] - η)
      F₀ .= cdf.(Logistic(), Y₀)
      F₁ .= cdf.(Logistic(), Y₁)
      pr .= F₁ .- F₀
      f₀ .= pdf.(Logistic(), Y₀)
      f₁ .= pdf.(Logistic(), Y₁)
      G[xs] .= X' * (wts .* (f₁ .- f₀) ./ pr)
      g2 = - (δ₁ .* f₁ .- δ₀ .* f₀)' * (wts ./ pr)
      θ₁[2:end] .= exp.(β[ζs[2:end]])
      G[ζs] .= diagm((d => θ₁[1:ks - d] for d ∈ 0:ks - 1)...) * g2
    end
    β₀ = vcat(zero(xs), log.(1:ks))
    td = TwiceDifferentiable(f, g!, β₀, autodiff = :forwarddiff)
    β = optimize(td, β₀) |> minimizer
    H = hessian!(td, β)
    Ψ = inv(bunchkaufman!(Hermitian(H)))
    A = Matrix(1.0I, size(Ψ)...)
    A[ζs, ζs] .= diagm((d => θ₁[1:ks + d] for d ∈ 0:-1:1 - ks)...)
    Ψ = Hermitian(A * Ψ * A')
    ŷ = X * β[1:size(X, 2)]
    β[n + 1:end] .= cumsum(vcat(β[n + 1], exp.(β[n + 2:end])))
	@pack! obj = X, y, β, Ψ, ŷ
	obj.vars = (obj.vars[1], obj.vars[2][lin_ind .+ 1])
	obj
end


text:
module InferenceBase

include("util.jl")
include("kernels.jl")
include("auto.jl")

end # module

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
#
###### Data transformations
#

# Lag data takes in a matrix. If we had columns a b c, then we will have a matrix with
# a(-lags) b(-lags) c(-lags) a(-lags + 1) b(-lags + 1) c(-lags + 1) ... ... ... a(-1) b(-1) c(-1).
# Maybe the columns should be reordered to have the closest lags at the front? changed!

function lagData(data, lags::Int64, obs::Int64)
	VARdata = data[lags:(obs-1),:]
	for i=(lags-1):-1:1
		VARdata = [VARdata data[i:(obs-lags+i-1),:]]
	end
	return VARdata
end

# dataMatrix is a function that creates a data matrix that goes into the estimation
# args:
#	- data: matrix of numbers
# 	- lags: number of lags in the estimation
# 	- typ: ASCIIString describing the type of VAR to be fitted, can be either of the following
# 			- None, do not include neither constant nor trend
# 			- Const, do include a constant
# 			- Trend, do include a trend
# 			- Const and trend, do include both constant and trend
# Returns a data matrix as from lagData with first columns corresponding to type.

function dataMatrix(data, lags, typ)
	obs = size(data)[1]
	data = lagData(data, lags, obs)

	typ=="None" && return(data)
	typ=="Const" && return([ones(obs-lags) data])
	typ=="Const and trend" && return([ones(obs-lags) 1:(obs-lags) data])
	typ=="Trend" && return([1:(obs-lags) data])
end

#
######
#


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

using Markdown
using InteractiveUtils

# ╔═╡ aaad96c0-801e-11ec-0b38-399812dced78
begin
	using PlutoUI
	using Test
	bigbreak = HTML("<br/>"^4)
end

# ╔═╡ aaad96c0-801e-11ec-27fe-bb16ccca80d2
md"""
# Strings

Topics:
1. How to get a string
2. String interpolation
3. String concatenation

"""

# ╔═╡ 10e3ed4e-53ee-4042-8538-a9ec15a416a4
md"""## How to get a string
Enclose your characters in `" "` or `\"\"\" \"\"\"`.

_(Single quotes  `' '` are exclusive for characters, not strings)_ """

# ╔═╡ 71fff504-34f1-4b9e-9d99-ca5e7a245f12
s1 = "I am a string."

# ╔═╡ 97ea92ed-d68f-40b7-8e3f-bd9e27d3a27c
s2 = """I am also a string."""

# ╔═╡ 0cd07f4a-8ab4-4e81-a2ee-0f7b64456232
md"
There are a couple functional differences between strings enclosed in single and triple quotes.
One difference is that, in the latter case, you can use quotation marks within your string.
"

# ╔═╡ d0001d61-e7f4-41ff-9095-4b6d5e8c15f9
"Here, we get an "error" because it's ambiguous where this string ends "

# ╔═╡ 6676d221-4a5a-4748-8f28-d865e910bfcd
"Look, Mom, no \"errors\"!!! "

# ╔═╡ 40791229-8e16-4eae-b1ce-59fc76d8b5b1
""" No "errors" here as well 😜 """

# ╔═╡ 2f819995-9494-40b4-b5be-42978dd7a207
md"**Note that ' ' define a character, but NOT a string!**"

# ╔═╡ 5f10ef9c-36da-4222-8c11-cf00f3d7c051
typeof('a')

# ╔═╡ b52aece1-3652-45e4-b292-0c7db7a37e05
bigbreak

# ╔═╡ 7c478f7a-e3cc-4c89-9339-31abe7bf358f
md"""
## String interpolation

We can use the `$` sign to insert existing variables into a string and to evaluate expressions within a string.
Below is an example that contains some highly sensitive personal information.
"""

# ╔═╡ 47c641a4-cb3f-4000-bebb-af762426dc7d
begin
	name = "Jane"
	num_fingers = 10
	num_toes = 10
end

# ╔═╡ 499b267a-7c8c-4def-a5ea-9f7e4ccae197
"Hello, my name is $name."

# ╔═╡ 702f1673-94b2-406a-b4d9-0aa7f8c2ad7d
"I have $num_fingers fingers and $num_toes toes."

# ╔═╡ 5c78e280-a5f8-495c-90ce-79d02d6bc635
"That is $(num_fingers + num_toes) digits in all!!"

# ╔═╡ 0e0fe89a-748b-4f0f-a189-52b5697dbc49
bigbreak

# ╔═╡ 1e4877f9-2e54-4fe0-be09-ccfed3daf301
md"""
## String concatenation

Below are three ways we can concatenate strings!

The first way is to use the `string()` function.
`string()` converts non-string inputs to strings.
"""

# ╔═╡ b9324e14-da0a-465e-b58c-ebd54eaf49e7
begin
	s3 = "How many cats ";
	s4 = "is too many cats?";
	😺 = 10
end

# ╔═╡ b50fb9f6-427e-4a39-9ca3-ba9e37c59b91
string(s3, s4)

# ╔═╡ a1b26473-0b27-467a-80f1-50a8121caa79
string("I don't know, but ", 😺, " is too few.")

# ╔═╡ 3f5915e4-7564-4395-beae-1074d86d3845
md""

# ╔═╡ 75b87fec-3acb-4442-9eeb-242198b2fa28
md"We can also use * for concatenation!"

# ╔═╡ 0fec7dc0-8ccf-4e25-abfd-30c995cfb3ae
s3*s4

# ╔═╡ 49e637b3-15a6-4220-8ffa-477532e61269
bigbreak

# ╔═╡ 183b4690-dbd1-4aa9-8afa-28f3c9e30e7b
md"""
## Exercises
2.1

Create a string that says "hi" 1000 times, first with `repeat` and then with the exponentiation operator, which can call * under the hood. Assign it the variable hi below.
"""

# ╔═╡ a5602b8b-c16a-493e-acc2-607938edd8f5
hi_repeat = "hi"

# ╔═╡ 5e6c11fe-e294-4480-b33c-1e19eb1afdb0
hi_power = "hi"

# ╔═╡ 2b387e43-7c33-4d65-82b5-3d63ca748e4c
@testset "2.1" begin
	result = "hihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihihi"
	@test hi_repeat == result
	@test hi_power == result
end

# ╔═╡ 3fc75a0f-cf3e-4893-b037-cbc9463d5222
md"""
2.2

Declare two variables

```julia
a = 3
b = 4
```
and use them to create two strings:
- `"3 + 4"`
- `"7"`
and store the results in `c` and `d` respectively
"""

# ╔═╡ 172769c7-9cfb-4fd4-ae91-d882e6c6ea5e
begin
	a = missing
	b = missing
	c = missing
	d = missing
end

# ╔═╡ 3dd734e2-6c5d-4f70-85da-7760eb4b64c1
@testset "2.2" begin
	@test c == "3 + 4"
	@test d == "7"
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[compat]
PlutoUI = "~0.7.32"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─aaad96c0-801e-11ec-0b38-399812dced78
# ╟─aaad96c0-801e-11ec-27fe-bb16ccca80d2
# ╟─10e3ed4e-53ee-4042-8538-a9ec15a416a4
# ╠═71fff504-34f1-4b9e-9d99-ca5e7a245f12
# ╠═97ea92ed-d68f-40b7-8e3f-bd9e27d3a27c
# ╟─0cd07f4a-8ab4-4e81-a2ee-0f7b64456232
# ╠═d0001d61-e7f4-41ff-9095-4b6d5e8c15f9
# ╠═6676d221-4a5a-4748-8f28-d865e910bfcd
# ╠═40791229-8e16-4eae-b1ce-59fc76d8b5b1
# ╟─2f819995-9494-40b4-b5be-42978dd7a207
# ╠═5f10ef9c-36da-4222-8c11-cf00f3d7c051
# ╟─b52aece1-3652-45e4-b292-0c7db7a37e05
# ╟─7c478f7a-e3cc-4c89-9339-31abe7bf358f
# ╠═47c641a4-cb3f-4000-bebb-af762426dc7d
# ╠═499b267a-7c8c-4def-a5ea-9f7e4ccae197
# ╠═702f1673-94b2-406a-b4d9-0aa7f8c2ad7d
# ╠═5c78e280-a5f8-495c-90ce-79d02d6bc635
# ╟─0e0fe89a-748b-4f0f-a189-52b5697dbc49
# ╟─1e4877f9-2e54-4fe0-be09-ccfed3daf301
# ╠═b9324e14-da0a-465e-b58c-ebd54eaf49e7
# ╠═b50fb9f6-427e-4a39-9ca3-ba9e37c59b91
# ╠═a1b26473-0b27-467a-80f1-50a8121caa79
# ╟─3f5915e4-7564-4395-beae-1074d86d3845
# ╟─75b87fec-3acb-4442-9eeb-242198b2fa28
# ╠═0fec7dc0-8ccf-4e25-abfd-30c995cfb3ae
# ╟─49e637b3-15a6-4220-8ffa-477532e61269
# ╟─183b4690-dbd1-4aa9-8afa-28f3c9e30e7b
# ╠═a5602b8b-c16a-493e-acc2-607938edd8f5
# ╠═5e6c11fe-e294-4480-b33c-1e19eb1afdb0
# ╠═2b387e43-7c33-4d65-82b5-3d63ca748e4c
# ╟─3fc75a0f-cf3e-4893-b037-cbc9463d5222
# ╠═172769c7-9cfb-4fd4-ae91-d882e6c6ea5e
# ╠═3dd734e2-6c5d-4f70-85da-7760eb4b64c1
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: slides/1 - Strings.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/1 - Strings.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/1 - Strings.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
"""
    load_OPFproblems(input_type, instance_path::String)

Read provided files in `instance_path` function of `input_type` and create a structure ds of type DataSource for each scenario. Also create empty structures gs and mp for each scenario
For each scenario, treat data in ds to complete gs : from bus and links in ds, generate node_linksin and node_linksout
For each scenario, complete mp structure
Return a dictonary Scenario => Scenario structure

# Arguments
- `input_type::T` : type of entry format for example MatpowerInput, GOCInput
- `instance_path::String` : path to the instance to treat, provided the instance path must be associated to an instance of format input_type

# Output
- `OPFpbs`

# Examples
```
julia > OPFpbs = load_OPFproblems(MatpowerInput, "instances\\matpower\\WB2.m")
julia > OPFpbs
Dict{String,Scenario} with 1 entry:
  "BaseCase" => Scenario
```
"""

function load_OPFproblems(input_type, instance_path::String)
  OPFpbs = read_input(input_type, instance_path)
  # baseMVA, bus, link, node_elems, link_elems, elem_to_typekey = read_input(input_type, filenames)
  # ds = DataSource(baseMVA, bus, link)

  for scenario in keys(OPFpbs)
    node_linksin, node_linksout = Dict{String, Set{Link}}(), Dict{String, Set{Link}}()

    for (linkname, _) in OPFpbs[scenario].ds.link
      if !haskey(node_linksin, linkname.dest)
        node_linksin[linkname.dest] = Set{Link}()
      end
      union!(node_linksin[linkname.dest], [linkname])

      if !haskey(node_linksout, linkname.orig)
        node_linksout[linkname.orig] = Set{Link}()
      end
      union!(node_linksout[linkname.orig], [linkname])
    end

    OPFpbs[scenario].gs.node_linksin = node_linksin
    OPFpbs[scenario].gs.node_linksout = node_linksout

    node_vars=Dict{String, Dict{String, Variable}}(busname => Dict{String, Variable}() for busname in keys(OPFpbs[scenario].ds.bus))
    link_vars=Dict{Link, Dict{String, Variable}}(link => Dict{String, Variable}() for link in keys(OPFpbs[scenario].ds.link))

    node_formulations = Dict{String, Dict{String, Symbol}}()
    link_formulations = Dict{Link, Dict{String, Symbol}}()
    for bus in keys(OPFpbs[scenario].ds.bus)
      node_formulations[bus] = Dict{String, Symbol}(elem => :NbMinVar for elem in keys(OPFpbs[scenario].ds.bus[bus]))
    end
    for link in keys(OPFpbs[scenario].ds.link)
      link_formulations[link] = Dict{String, Symbol}(elem => :NbMinVar for elem in keys(OPFpbs[scenario].ds.link[link]))
    end

    OPFpbs[scenario].mp = MathematicalProgramming(node_formulations, link_formulations, node_vars, link_vars)
  end

  return OPFpbs
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
include("../utils.jl")

@testset "isdef, isanydef" begin
    @test G.isdef(1)
    @test !G.isdef(nothing)

    struct Foo_isdef
        a::G.Option{Int}
        b::G.Option{Int}
    end
    Foo = Foo_isdef

    f1 = Foo(1,2)
    f2 = Foo(1,nothing)
    f3 = Foo(nothing,2)
    f4 = Foo(nothing,nothing)
    @test G.isanydef(f1) && G.isanydef(f2) && G.isanydef(f3)
    @test !G.isanydef(f4)
end

@testset "reset!" begin
    mutable struct Foo_reset
        a::Int
        b::Int
    end
    Foo_reset(; bar=1) = Foo_reset(1, bar)
    Foo = Foo_reset

    f  = Foo(3, 2)
    f2 = Foo(3, 2)
    f3 = Foo(3, 2)
    G.reset!(f)
    G.reset!(f2, exclude=[:b])
    G.reset!(f3, bar=7)
    @test f.a == 1 == f.b
    @test f2.a == 1 && f2.b == 2
    @test f3.a == 1 && f3.b == 7
end

@testset "nvec" begin
    mutable struct Foo_nvec
        a::Int
    end
    Foo_nvec() = Foo_nvec(1)
    Foo = Foo_nvec

    v = G.nvec(2, Foo)
    @test length(v) == 2
    @test eltype(v) == Foo_nvec
    @test v[1].a == v[2].a == 1
end


text:
using Random, StatsBase, Parameters, LinearAlgebra, Distributions
import Base: isempty
include("simulate.jl")

Random.seed!(0)

struct Job
    unique_id::Int64
end

# function Job(unique_id::Int64)
#     return Job(unique_id)
# end

struct Server
    unique_id::Int
    buffer::Queue{Job}           # First element in buffer is currently served job
    max_buffer_size::Int         # Buffer size before overflow
    overflow_vector::Vector{Float64} # Probabilities of where an overflown job will go
    overflow_leave::Float64
    move_vector::Vector{Float64}     # Probabilities of where a finished job will go
    move_leave::Float64
    service_rate::Float64        # Mean service rate of this server
end

function Server(unique_id::Int, max_buffer_size::Int,
    overflow_vector::Vector{Float64},
    move_vector::Vector{Float64},
    service_rate::Float64)
    return Server(unique_id::Int, Queue{Job}(), max_buffer_size, overflow_vector, 1-sum(overflow_vector), move_vector, 1-sum(move_vector), service_rate)
end

struct ArrivalEvent <: Event
    job::Job
end

"""
    ArrivalEvent(unique_id::Int64)

Creates a new `ArrivalEvent` and populates it with a new `Job`; the `Job`'s `unique_id` is as given
"""
function ArrivalEvent(unique_id::Int64)
    return ArrivalEvent(Job(unique_id))
end

abstract type TransferEvent <: Event end

struct OverflowEvent <: TransferEvent
    transfer_to::Server
    job::Job
end

struct MoveEvent <: TransferEvent
    transfer_to::Server
    job::Job
end

struct LeaveEvent <: Event
    job::Job
end

# This event is called when the server finishes serving
struct ServedEvent <: Event
    server::Server
end

#The @with_kw macro comes from the Parameters package
@with_kw struct NetworkParameters
    L::Int
    gamma_scv::Float64 #This is constant for all scenarios at 3.0
    λ::Float64 #This is undefined for the scenarios since it is varied
    η::Float64 #This is assumed constant for all scenarios at 4.0
    μ_vector::Vector{Float64} #service rates
    P::Matrix{Float64} #routing matrix
    Q::Matrix{Float64} #overflow matrix
    p_e::Vector{Float64} #external arrival distribution
    K::Vector{Int} #-1 means infinity
end

mutable struct DiscreteState <: State
    servers::Vector{Server}     # Enumerated list of servers
    left_jobs::Vector{Job}      # Enumerated list of jobs that have left the system
    moving_jobs::Vector{Job}    # Enumerated list of jobs that are moving across servers
    moving_jobs_num::Int        # Number of moving jobs (cheap memory-wise to keep them here)
    last_job::Int64             # Index of last-added job
    jobs_in_system::Int64       # Number of total jobs in sys (also cheap memory-wise)
    is_debug::Bool
    params::NetworkParameters   # Hold our params in the state (is this a good idea? we'll find out)
end

"""
    DiscreteState(params::NetworkParameters)

Creates a DiscreteState, populates it with `servers` and `params` based on the given. First `job` index is 1, and all `job` storages start empty
"""
function DiscreteState(params::NetworkParameters)
    servers = [Server(i, params.K[i], Vector{Float64}(vec(params.Q[i:params.L:end])), Vector{Float64}(vec(params.P[i:params.L:end])), params.μ_vector[i]) for i in 1:params.L]
    return DiscreteState(servers, Vector{Job}(), Vector{Job}(), 0, 1, 0, false, params)
end

"""
Pull out the important info from our state - we don't need `params`
"""
function copy_state_items(state::DiscreteState)
    return deepcopy(state.servers), deepcopy(state.left_jobs)
end

"""
Returns whether a `Server`'s buffer is full
"""
function isfull(server::Server)
    return server.max_buffer_size != -1 && length(server.buffer) > server.max_buffer_size
end

"""
Returns whether a `Server` has no current `Job` and empty buffer
"""
function isempty(server::Server)
    return isempty(server.buffer)
end

"""
Returns index of `Server` a `Job` should transfer to, or 0 if the `Job` should leave
"""
function sample_probabilities(transfer_probs::Vector{Float64}, chance_leave::Float64)::Int
    return rand() < chance_leave ? 0 : sample(1:length(transfer_probs), Weights(transfer_probs))
end

"""
Returns server index a job should overflow to, or 0 if job should leave
"""
function get_overflow_index_from(server::Server)::Int
    return sample_probabilities(server.overflow_vector, server.overflow_leave)
end

"""
Returns server index a job should move to, or 0 if job should leave
"""
function get_move_index_from(server::Server)::Int
    return sample_probabilities(server.move_vector, server.move_leave)
end

"""
Returns the Server to transfer (overflow, move) Job to, or nothing if Job is to leave
"""
function get_transfer_from_index(state::DiscreteState, index::Int)
    return state.servers[index]
end

"""
Process a job leaving the system. Simply adds it to the state's array of left jobs
"""
function leave_job(state::DiscreteState, job::Job)
    state.jobs_in_system -= 1
    # push!(state.left_jobs, job)
    # # state.is_debug && println("Job #", job.unique_id, " left")
end

"""
Handles moving a Job between Servers.
If Job is to leave:
 - Processes its departure, returns nothing
If it's to move:
 - Returns the appropriate MoveEvent or OverflowEvent (determined by is_overflow flag)
"""
function process_transfer_from(state::DiscreteState, server::Server, job::Job, is_overflow::Bool)
    ind = is_overflow ? get_overflow_index_from(server) : get_move_index_from(server)
    if ind == 0
        leave_job(state, job)
        return nothing
    else
        state.moving_jobs_num += 1
        return is_overflow ? OverflowEvent(get_transfer_from_index(state, ind), job) : MoveEvent(get_transfer_from_index(state, ind), job)
    end
end

"""
Handles adding a Job to a Server.
If buffer is full:
 - Return OverflowEvent if moving to different buffer
 - Return nothing if leaving system
If buffer is empty:
 - Add Job to Server, return ServedEvent
Else:
 - Add Job to Server buffer, return nothing
"""
function add_job(state::DiscreteState, server::Server, job::Job)
    if isfull(server)
        return process_transfer_from(state, server, job, true)
    else
        ret = isempty(server) ? ServedEvent(server) : nothing
        enqueue!(server.buffer, job)
        return ret
    end
end

"""
For when a server finishes processing its job
"""
function process_job(server::Server)
    return dequeue!(server.buffer)
end

"""
    process_event(time::Float64, state::DiscreteState, event::ArrivalEvent)

Processes a Job's arrival:
 - Schedules the next Job's arrival
 - Attempts to add this job to a server
"""
function process_event(time::Float64, state::DiscreteState, event::ArrivalEvent)

    # # state.is_debug && println("Job #", state.last_job, " Arrived")
    state.last_job += 1
    state.jobs_in_system += 1
    new_timed_events = TimedEvent[]

    # Prepare next arrival
    push!(new_timed_events,generate_timed_event(time, state, ArrivalEvent(state.last_job)))

    server = state.servers[sample_probabilities(state.params.p_e, .0)]  # We will need to assert sum(p_e) = 1

    new_event = add_job(state, server, event.job)
    !isnothing(new_event) && push!(new_timed_events, generate_timed_event(time,state, new_event))
    # # If this is the only job on the server
    # state.number_in_system == 1 && push!(new_timed_events,TimedEvent(EndOfServiceEvent(), time + 1/μ))
    # Increase number in system

    return new_timed_events
end

"""
    process_event(time::Float64, state::DiscreteState, event::ServedEvent)

Handles a `Job`'s completion:
- Removes `Job` from Server
 - If buffer's not empty, create new `ServedEvent` for `Server`
- Handles `Job`'s subsequent path
 - If it moves to another `Server`, create `TransferEvent`
 - Else, do nothing
"""
function process_event(time::Float64, state::DiscreteState, event::ServedEvent)
    job = process_job(event.server)
    # # state.is_debug && println("Served Job #", job.unique_id, " at Server #", event.server.unique_id)
    new_timed_events = TimedEvent[]
    !isempty(event.server) && push!(new_timed_events, generate_timed_event(time, state, ServedEvent(event.server)))

    new_event = process_transfer_from(state, event.server, job, false)

    # Decide which server to send the job to, based on the event Vector
        # If sent out of system, return that event instead
    !isnothing(new_event) && append!(new_timed_events, [generate_timed_event(time, state, new_event)])
    # Decide whether to move normally or overflow
        # Return the appropriate event
    return new_timed_events
end

"""
    process_event(time::Float64, state::DiscreteState, event::TransferEvent)

Attempts to add a `Job` to a `Server`'s buffer, handles results
"""
function process_event(time::Float64, state::DiscreteState, event::TransferEvent)
    job = event.job
    state.moving_jobs_num -= 1
    server_to = event.transfer_to
    # # state.is_debug && println("Transferring Job #", job.unique_id, " to Server #", server_to.unique_id)
    new_event = add_job(state, server_to, job)
    return !isnothing(new_event) ? TimedEvent[generate_timed_event(time, state, new_event)] : TimedEvent[]
end

"""
    generate_timed_event(time::Float64, state::DiscreteState, event::TransferEvent)

Generates a `TimedEvent` for a `TransferEvent`. Its `time` is the current + global gamma distributed `η`
"""
function generate_timed_event(time::Float64, state::DiscreteState, event::TransferEvent)::TimedEvent
    return TimedEvent(event, time + generate_time(state.params.η, state.params.gamma_scv))
end

"""
    generate_timed_event(time::Float64, state::DiscreteState, event::ServedEvent)

Generates a `TimedEvent` for a `ServedEvent`. Its `time` is the current + gamma distributed `Server`'s rate
"""
function generate_timed_event(time::Float64, state::DiscreteState, event::ServedEvent)::TimedEvent
    return TimedEvent(event, time + generate_time(event.server.service_rate, state.params.gamma_scv))
end

"""
    generate_timed_event(time::Float64, state::DiscreteState, event::ArrivalEvent)

Generates a `TimedEvent` for an `ArrivalEvent`. Its `time` is the current + global arrival rate (will not be constant)
"""
function generate_timed_event(time::Float64, state::DiscreteState, event::ArrivalEvent)::TimedEvent
    return TimedEvent(event, time + generate_time(state.params.λ, state.params.gamma_scv))
end

"""
A convenience function to make a Gamma distribution with desired rate (inverse of shape) and SCV.
"""
function rate_scv_gamma(desired_rate::Float64, desired_scv::Float64)
    return Gamma(1/desired_scv, desired_scv/desired_rate)
end

"""
Generates a time based on mean rate μ and SCV gamma_scv
"""
function generate_time(μ::Float64, gamma_scv::Float64)
    rand(rate_scv_gamma(μ, gamma_scv))
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SetPartitionsTypeM

using Nemo, IntPartitions, OrderedSetPartitionsTypeM, GeneralizedBinomial
using DataStructures, Products, Polynomials

export ModuleSetPartitionsOfTypeM
export ShapePartitionsOfTypeM, SetPartitionsOfTypeM
# by shape
export L036040, L257490, L327003, L327004
# by length # X008284
export L156289, L291451, L291452
# diagonal
export L001147, L025035, L025036
# row sum
export L005046, L291973, L291975
# alt row sum
export V153881, L000587, L260884, V291974, V291976
# central
export L007820, L327416, L327417, L327418
# see also A260876.
export L014606, V291973, V291975, V005046, V260884, V000587, L291976, L291974

"""

For example consider the case n = 4. There are five integer partitions of 4:

* P = [[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]]. The shapes are m times the parts of the integer partitions: S(m) = [[4m], [3m, m], [2m, 2m], [2m, m, m], [m, m, m, m]].

* In the case m = 1 we look at set partitions of {1, 2, 3, 4} with sizes ∈  [[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]] which gives rise to [1, 4, 3, 6, 1] with sum 15.

* In the case m = 2 we look at set partitions of {1, 2, .., 8} with sizes ∈ [[8], [6, 2], [4, 4], [4, 2, 2], [2, 2, 2, 2]] which gives rise to [1, 28, 35, 210, 105] with sum 379.

* In the case m = 0 we look at set partitions of {} with sizes ∈ [[0], [0, 0], [0, 0], [0, 0, 0], [0, 0, 0, 0]] which gives rise to [1, 1, 1, 1, 1] with sum 5 (because the only partition of the empty set is the set that contains the empty set, thus from the definition T(0,4) = Sum_{S(0)} card({0}) = A000041(4) = 5).

* If n runs through 0, 1, 2,... then the result is an irregular triangle ∈ which the n-th row lists multinomials for partitions of [m*n] which have only parts which are multiples of m. These are the triangles A080575 (m = 1), A257490 (m = 2), A327003 (m = 3), A327004 (m = 4). In the case m = 0 the triangle is A000012 subdivided into rows of length A000041. See the  references below how this case integrates into the full picture.

| type  | m = 0 | m = 1 | m = 2 | m = 3 | m = 4 |
|-------|-------|-------|-------|-------|-------|
| by shape | [A000012](https://oeis.org/A000012) | [A036040](https://oeis.org/A036040) | [A257490](https://oeis.org/A257490) | [A327003](https://oeis.org/A327003) | [A327004](https://oeis.org/A327004) |
| by length | [A008284](https://oeis.org/A008284) | [A048993](https://oeis.org/A048993) | [A156289](https://oeis.org/A156289) | [A291451](https://oeis.org/A291451) | [A291452](https://oeis.org/A291452) |
| diagonal | [A000012](https://oeis.org/A000012) | [A000012](https://oeis.org/A000012) | [A001147](https://oeis.org/A001147) | [A025035](https://oeis.org/A025035) | [A025036](https://oeis.org/A025036) |
| row sum | [A000041](https://oeis.org/A000041) | [A000110](https://oeis.org/A000110) | [A005046](https://oeis.org/A005046) | [A291973](https://oeis.org/A291973) | [A291975](https://oeis.org/A291975) |
| alt row sum | [A081362](https://oeis.org/A081362) | [A000587](https://oeis.org/A000587) | [A260884](https://oeis.org/A260884) | [A291974](https://oeis.org/A291974) | [A291976](https://oeis.org/A291976) |
| central | [A000041](https://oeis.org/A000041) | [A007820](https://oeis.org/A007820) | [A327416](https://oeis.org/A327416) | [A327417](https://oeis.org/A327417) | [A327418](https://oeis.org/A327418) |

See also [A260876](https://oeis.org/A260876).
"""
const ModuleSetPartitionsOfTypeM = ""

"""
--
"""
function ShapePartitionsOfTypeM(m, n, k)
    shapes = (map(x -> x*m, p) for p ∈ IntegerPartitions(n, k))
    isempty(shapes) && return [fmpz(0)]
    return [CardinalityOfShapePartitions(s) for s ∈ shapes]
end

"""

Return the number of shape-partitions of type m of an mn-set.
# Examples
```julia-repl
for n in 0:5 ShapePartitionsOfTypeM(2, n) |> Println end
[1]
[0, 1]
[0, 1, 3]
[0, 1, 15, 15]
[0, 1, 28, 35, 210, 105]
[0, 1, 45, 210, 630, 1575, 3150, 945]
```
"""
function ShapePartitionsOfTypeM(m, n)
    reduce(vcat, [ShapePartitionsOfTypeM(m, n, k) for k ∈ 0:n])
end

"""

Return the list of the number of set partitions of a mn set into k parts which are m-sized.
"""
function SetPartitionsOfTypeM(m, n)
    m == 0 && return [k == 1 ? 1 : 0 for k in 0:n]
    EgfCoeffs(OrderedSetPolynomials(m, n))
end

"""

Return the number of partitions of ``{1, 2, ..., mn}`` into sets of size ``m``.
"""
SetPartitionsOfSizeTypeM(m, len) = [SetPartitionsOfTypeM(m, n)[n+1] for n in 0:len-1]

# A060540 T(n,k) = (n*k)!/(k!^n*n!), (n>=1, k>=1).
# the number of ways of dividing nk labeled items into n unlabeled boxes with k items in each box.

"""

Return the number of partitions of ``{1, 2, ..., 2n}`` into sets of size ``2``.
# Examples
```julia-repl
julia> L001147(7)
7-element Array{fmpz,1}:
[1, 1, 3, 15, 105, 945, 10395]
```
"""
L001147(len) = SetPartitionsOfSizeTypeM(2, len)

"""

Return the number of partitions of ``{1, 2, ..., 3n}`` into sets of size ``3``.
# Examples
```julia-repl
julia> L025035(7)
7-element Array{fmpz,1}:
[1, 1, 10, 280, 15400, 1401400, 190590400]
```
"""
L025035(len) = SetPartitionsOfSizeTypeM(3, len)

"""

Return the number of partitions of ``{1, 2, ..., 4n}`` into sets of size ``4``.
# Examples
```julia-repl
julia> L025036(7)
7-element Array{fmpz,1}:
[1, 1, 35, 5775, 2627625, 2546168625, 4509264634875]
```
"""
L025036(len) = SetPartitionsOfSizeTypeM(4, len)

"""

Return the number of set partitions of a mn set into k parts which are m-sized.
"""
SetPartitionsOfTypeM(m, n, k) = SetPartitionsOfTypeM(m, n)[k+1]

"""

Return the number of set partitions of n of type 0 ordered by length .
# Examples
```julia-repl
julia> L008284(10)
7-element Array{fmpz,1}:
[0, 1, 43690, 7128576, 99379280, 285885600, 190590400]
```
"""
L008284(n) = SetPartitionsOfTypeM(0, n)

"""

Return the number of end rhyme patterns of a poem of an even number of lines.
# Examples
```julia-repl
julia> L156289(6)
7-element Array{fmpz,1}:
[ 0, 1, 1023, 21120, 65835, 51975, 10395]
```
"""
L156289(n) = SetPartitionsOfTypeM(2, n)

"""

Return the number of set partitions of type 3.
# Examples
```julia-repl
julia> L291451(10)
7-element Array{fmpz,1}:
[0, 1, 43690, 7128576, 99379280, 285885600, 190590400]
```
"""
L291451(n) = SetPartitionsOfTypeM(3, n)

"""

``exp(x (cos(z) + cosh(z) - 2)/2)``
# Examples
```julia-repl
L291452(6)
7-element Array{fmpz,1}:
[0, 1, 2098175, 2941884000, 181262956875, 1932541986375, 4509264634875]
```
"""
L291452(n) = SetPartitionsOfTypeM(4, n)

"""

Return the number of partitions of a 2n-set into even blocks.
# Examples
```julia-repl
julia> V005046(30)
5215270659830794008813481789355446576694940667184287549
```
"""
V005046(n) = sum(SetPartitionsOfTypeM(2, n))

"""

Return a list of the number of partitions of even sized sets into even blocks.
# Examples
```julia-repl
julia> L005046(10)
10-element Array{fmpz,1}:
[ 1, 1, 4, 31, 379, 6556, 150349, 4373461, 156297964, 6698486371]
```
"""
L005046(len) = [V005046(n) for n ∈ 0:len-1]

"""

``(3n)! [z^(3n)] exp(exp(z)/3 + 2exp(-z/2)cos(z sqrt(3)/2)/3 - 1).``
# Examples
```julia-repl
julia> V291973(9)
31728742163212641
```
"""
V291973(n) = sum(SetPartitionsOfTypeM(3, n))

"""

Return the n-th complementary Bell number.
# Examples
```julia-repl
julia> V000587(50)
-29113173035759403920216141265491160286912
```
"""
V000587(n) = sum(AltEgfCoeffs(OrderedSetPolynomials(1, n)))

SignedShapePartitions(m, n) = reduce(vcat, [(-1)^k*ShapePartitionsOfTypeM(m, n, k) for k ∈ 0:n])

"""

Return ...
"""
V153881(n) = sum(SignedShapePartitions(0, n))

"""

Return a list of complementary Bell numbers of length ``len``.
# Examples
```julia-repl
julia> L000587(10)
10-element Array{fmpz,1}:
[1, -1, 0, 1, 1, -2, -9, -9, 50, 267]
```
"""
L000587(len) = [V000587(n) for n ∈ 0:len-1]

"""

Return the number of set partitions of a 2n-set into even blocks which have even length minus the number of partitions into even blocks which have odd length.
# Examples
```julia-repl
julia> V260884(19)
5097857816569586800024019
```
"""
V260884(n) = sum(AltEgfCoeffs(OrderedSetPolynomials(2, n)))

"""

Return a list of the first ``len`` terms V260884(n).
# Examples
```julia-repl
julia> L260884(10)
10-element Array{fmpz,1}:
[1, -1, 2, -1, -43, 254, 4157, -70981, -1310398, 40933619]
```
"""
L260884(len) = [V260884(n) for n ∈ 0:len-1]

"""

(3n)! [z^(3n)] exp(-(exp(z)/3 + 2exp(-z/2) cos(z sqrt(3)/2)/3 - 1)).
# Examples
```julia-repl
julia> V291974(9)
-3166484321001
```
"""
V291974(n) = sum(AltEgfCoeffs(OrderedSetPolynomials(3, n)))

"""

# Examples
```julia-repl
julia> L291974(6)
6-element Array{fmpz,1}:
[1, -1, 9, -197, 6841, -254801]
```
"""
L291974(len) = [V291974(n) for n ∈ 0:len-1]

"""

Return ordered the number of set partitions of type 3.
# Examples
```julia-repl
julia> V291975(9)
926848347928901638652131
```
"""
V291975(n) = sum(SetPartitionsOfTypeM(4, n))

"""

# Examples
```julia-repl
L291975(5)
5-element Array{fmpz,1}:
[1, 1, 36, 6271, 3086331]
```
"""
L291975(len) = [V291975(n) for n ∈ 0:len-1]

"""

Return the number of ordered set partitions of type 4.
# Examples
```julia-repl
julia> V291976(7)
-6440372006517541
```
"""
V291976(n) = sum(AltEgfCoeffs(OrderedSetPolynomials(4, n)))

"""

Return the number of ordered set partitions of type 4.
# Examples
```julia-repl
L291976(7)
7-element Array{fmpz,1}:
[ 1, -1, 34, -5281, 2185429, -1854147586, 2755045819549]
```
"""
L291976(len) = [V291976(n) for n ∈ 0:len-1]

"""

# Examples
```julia-repl
L014606(7)
7-element Array{fmpz,1}:
[1, 1, 20, 1680, 369600, 168168000, 137225088000]
```
"""
L014606(len) = Diagonal(n -> OrderedSetPolynomials(3, n), len)

"""

Return
# Examples
```julia-repl
julia> L291973(6)
6-element Array{fmpz,1}:
[1, 1, 11, 365, 25323, 3068521]
```
"""
L291973(len) = [V291973(n) for n ∈ 0:len-1]

"""

Return
# Examples
```julia-repl
for n in 0:5 SetPartitionsOfTypeM(0, n) |> Println end
[1]
[0, 1]
[0, 1, 0]
[0, 1, 0, 0]
[0, 1, 0, 0, 0]
[0, 1, 0, 0, 0, 0]
```
"""
Lnix(n) = SetPartitionsOfTypeM(0, n)

"""

Return the number of set partitions of n with block sizes given by the k-th integer partition of n when the integer partitions are sorted in Hindenburg order. This is row n of the triangle of multinomial coefficients.
# Examples
```julia-repl
for n in 0:5 SetPartitionsOfTypeM(1, n) |> Println end
[1]
[0, 1]
[0, 1, 1]
[0, 1, 3, 1]
[0, 1, 7, 6, 1]
[0, 1, 15, 25, 10, 1]
```
"""
L036040(n) = SetPartitionsOfTypeM(1, n)

"""

Return the number of set partitions of 2n with block sizes given by the even parts of integer partition of 2n when the integer partitions are sorted in Hindenburg order.
# Examples
```julia-repl
for n in 0:5 SetPartitionsOfTypeM(2, n) |> Println end
[1]
[0, 1, 3]
[0, 1, 15, 15]
[0, 1, 63, 210, 105]
[0, 1, 255, 2205, 3150, 945]
```
"""
L257490(n) = SetPartitionsOfTypeM(2, n)

"""

Return the number of set partitions of 3n with block sizes given by the parts of integer partition of 3n which are multiples of 3 and when the integer partitions are sorted in Hindenburg order.
# Examples
```julia-repl
for n in 0:5 SetPartitionsOfTypeM(3, n) |> Println end
[1]
[0, 1]
[0, 1, 10]
[0, 1, 84, 280]
[0, 1, 682, 9240, 15400]
[0, 1, 5460, 260260, 1401400, 1401400]
```
"""
L327003(n) = SetPartitionsOfTypeM(3, n)

"""

Return the number of set partitions of 4n with block sizes given by the parts of integer partition of 4n which are multiples of 4 and when the integer partitions are sorted in Hindenburg order.
# Examples
```julia-repl
for n in 0:5 SetPartitionsOfTypeM(4, n) |> Println end
[1]
[0, 1]
[0, 1, 35]
[0, 1, 495, 5775]
[0, 1, 8255, 450450, 2627625]
[0, 1, 130815, 35586525, 727476750, 2546168625]
```
"""
L327004(n) = SetPartitionsOfTypeM(4, n)

"""

Return the central column of the triangle of set partitions of 1-type. (Stirling set numbers).
"""
L007820(len) = [SetPartitionsOfTypeM(1, 2n)[n+1]  for n in 0:len-1]

"""

Return the central column of the triangle of set partitions of 2-type.
"""
L327416(len) = [SetPartitionsOfTypeM(2, 2n)[n+1]  for n in 0:len-1]

"""

Return the central column of the triangle of set partitions of 3-type.
"""
L327417(len) = [SetPartitionsOfTypeM(3, 2n)[n+1]  for n in 0:len-1]

"""

Return the central column of the triangle of set partitions of 4-type.
"""
L327418(len) = [SetPartitionsOfTypeM(4, 2n)[n+1]  for n in 0:len-1]

#START-TEST-########################################################

using Test, SeqTests, SeqUtils

function test()

    @testset "SetPartMType" begin
        println("SetPartMType")
        L036040(6) |> println
        L257490(6) |> println
        L327003(6) |> println
        L327004(6) |> println
        # X008284, X048993,
        L156289(6) |> println
        @test L156289(5) == [0, 1, 255, 2205, 3150, 945]
        @test L156289(6) == [0, 1, 1023, 21120, 65835, 51975, 10395]
        L291451(6) |> println
        @test L291451(5) == [0, 1, 5460, 260260, 1401400, 1401400]
        L291452(6) |> println
        @test L291452(5) == [0, 1, 130815, 35586525, 727476750, 2546168625]
        V005046(6) |> println
        L005046(6) |> println
        @test V005046(7) == 4373461
        @test L005046(7) == [1, 1, 4, 31, 379, 6556, 150349]
        V291973(6) |> println
        L291973(6) |> println
        @test L291973(7) == [1, 1, 11, 365, 25323, 3068521, 583027547]
        V291975(6) |> println
        L291975(6) |> println

        V000587(6) |> println
        L000587(6) |> println
        V260884(6) |> println
        L260884(6) |> println
        V291974(6) |> println
        V291976(6) |> println

        @test L007820(6) == [1, 1, 7, 90, 1701, 42525]
        @test L327416(6) == [1, 1, 63, 21120, 20585565, 44025570225]
        @test L327417(6) == [1, 1, 682, 7128576, 429120851544, 94066556834970720]
        @test L327418(6) == [1, 1, 8255, 2941884000, 11957867341948125, 294040106448733743008625]

        L014606(6) |> println
        @test L291976(6) == [1, -1, 34, -5281, 2185429, -1854147586]
        @test L260884(7) == [1, -1, 2, -1, -43, 254, 4157]
        @test V000587(7) == -9
        @test L000587(10) == [1, -1, 0, 1, 1, -2, -9, -9, 50, 267]
        @test L291974(7) == [1, -1, 9, -197, 6841, -254801, -3000807]
        @test L291975(6) == [1, 1, 36, 6271, 3086331, 3309362716]

        @test L001147(6) == [1, 1, 3, 15, 105, 945]
        @test L025035(6) == [1, 1, 10, 280, 15400, 1401400]
        @test L025036(6) == [1, 1, 35, 5775, 2627625, 2546168625]

    end

    @testset "SetPartByShape" begin
        @test ShapePartitionsOfTypeM(3, 4) == [0, 1, 220, 462, 9240, 15400]
    end
end

function demo()

    for n ∈ 0:4
        [ShapePartitionsOfTypeM(3, n, k) for k ∈ 0:n] |> println
    end

    for m ∈ 0:4
        println("\nUnordered ShapePartitions m = $m")
        for n ∈ 0:6
            ShapePartitionsOfTypeM(m, n) |> Println
        end
    end

    for m ∈ 0:4
        println("\nUnordered SetPartitions m = $m")
        for n ∈ 0:6
            SetPartitionsOfTypeM(m, n) |> Println
        end
    end

end

function perf() end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
"""
 Return the Stirling numbers of 2nd kind.
 Exported from Module StirlingLahNumbers
"""
L048993(n) = CoeffExp(OrderedSetPolynomials(1, n))
=#

#=
"""

Return the number of set partitions of type 1.
# Examples
```julia-repl
julia> LA048993(10)
7-element Array{fmpz,1}:
[0, 1, 43690, 7128576, 99379280, 285885600, 190590400]
```
"""
L048993(n) = SetPartitionsOfTypeM(1, n)
=#

#=

('--', 0, '---------------------')

by Shape    A000012
[0] [1]
[1] [1]
[2] [1, 1]
[3] [1, 1, 1]
[4] [1, 1, 1, 1, 1]
[5] [1, 1, 1, 1, 1, 1, 1]
[6] [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

by Length   A008284
[0] [1]
[1] [0, 1]
[2] [0, 1, 1]
[3] [0, 1, 1, 1]
[4] [0, 1, 2, 1, 1]
[5] [0, 1, 2, 2, 1, 1]
[6] [0, 1, 3, 3, 2, 1, 1]

[1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 2, 1, 1, 0, 1, 2, 2, 1, 1]

RowSum      A000041 [1, 1, 2, 3, 5, 7, 11]
AlterRowSum A081362 [1, -1, 0, -1, 1, -1, 1]
Diagonal    A000012 [1, 1, 1, 1, 1, 1, 1]
Central     A000041 [1, 1, 2, 3, 5, 7]

('--', 1, '---------------------')

by Shape    A036040
[0] [1]
[1] [1]
[2] [1, 1]
[3] [1, 3, 1]
[4] [1, 4, 3, 6, 1]
[5] [1, 5, 10, 10, 15, 10, 1]
[6] [1, 6, 15, 10, 15, 60, 15, 20, 45, 15, 1]

[1, 1, 1, 1, 1, 3, 1, 1, 4, 3, 6, 1]

by Length   A048993
[0] [1]
[1] [0, 1]
[2] [0, 1, 1]
[3] [0, 1, 3, 1]
[4] [0, 1, 7, 6, 1]
[5] [0, 1, 15, 25, 10, 1]
[6] [0, 1, 31, 90, 65, 15, 1]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 1]
[3] [0, 1, 3, 1]
[4] [0, 1, 7, 6, 1]
[5] [0, 1, 15, 25, 10, 1]
[6] [0, 1, 31, 90, 65, 15, 1]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 2]
[3] [0, 1, 6, 6]
[4] [0, 1, 14, 36, 24]
[5] [0, 1, 30, 150, 240, 120]
[6] [0, 1, 62, 540, 1560, 1800, 720]

[1, 0, 1, 0, 1, 1, 0, 1, 3, 1, 0, 1, 7, 6, 1, 0, 1, 15, 25, 10, 1]

RowSum      A000110 [1, 1, 2, 5, 15, 52, 203]
AlterRowSum A000587 [1, -1, 0, 1, 1, -2, -9]
Diagonal    A000012 [1, 1, 1, 1, 1, 1, 1]
Central     A007820 [1, 1, 7, 90, 1701, 42525]

('--', 2, '---------------------')

by Shape    A257490
[0] [1]
[1] [1]
[2] [1, 3]
[3] [1, 15, 15]
[4] [1, 28, 35, 210, 105]
[5] [1, 45, 210, 630, 1575, 3150, 945]
[6] [1, 66, 495, 462, 1485, 13860, 5775, 13860, 51975, 51975, 10395]

[1, 1, 1, 3, 1, 15, 15, 1, 28, 35, 210, 105]

by Length   A156289
[0] [1]
[1] [0, 1]
[2] [0, 1, 3]
[3] [0, 1, 15, 15]
[4] [0, 1, 63, 210, 105]
[5] [0, 1, 255, 2205, 3150, 945]
[6] [0, 1, 1023, 21120, 65835, 51975, 10395]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 3]
[3] [0, 1, 15, 15]
[4] [0, 1, 63, 210, 105]
[5] [0, 1, 255, 2205, 3150, 945]
[6] [0, 1, 1023, 21120, 65835, 51975, 10395]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 6]
[3] [0, 1, 30, 90]
[4] [0, 1, 126, 1260, 2520]
[5] [0, 1, 510, 13230, 75600, 113400]
[6] [0, 1, 2046, 126720, 1580040, 6237000, 7484400]

[1, 0, 1, 0, 1, 3, 0, 1, 15, 15, 0, 1, 63, 210, 105, 0, 1, 255, 2205, 3150, 945]

RowSum      A005046 [1, 1, 4, 31, 379, 6556, 150349]
AlterRowSum A260884 [1, -1, 2, -1, -43, 254, 4157]
Diagonal    A001147 [1, 1, 3, 15, 105, 945, 10395]
Central     A327416 [1, 1, 63, 21120, 20585565, 44025570225]

('--', 3, '---------------------')

by Shape    A327003
[0] [1]
[1] [1]
[2] [1, 10]
[3] [1, 84, 280]
[4] [1, 220, 462, 9240, 15400]
[5] [1, 455, 5005, 50050, 210210, 1401400, 1401400]
[6] [1, 816, 18564, 24310, 185640, 4084080, 2858856, 13613600, 85765680, 285885600, 190590400]

[1, 1, 1, 10, 1, 84, 280, 1, 220, 462, 9240, 15400]

by Length   A291451
[0] [1]
[1] [0, 1]
[2] [0, 1, 10]
[3] [0, 1, 84, 280]
[4] [0, 1, 682, 9240, 15400]
[5] [0, 1, 5460, 260260, 1401400, 1401400]
[6] [0, 1, 43690, 7128576, 99379280, 285885600, 190590400]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 10]
[3] [0, 1, 84, 280]
[4] [0, 1, 682, 9240, 15400]
[5] [0, 1, 5460, 260260, 1401400, 1401400]
[6] [0, 1, 43690, 7128576, 99379280, 285885600, 190590400]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 20]
[3] [0, 1, 168, 1680]
[4] [0, 1, 1364, 55440, 369600]
[5] [0, 1, 10920, 1561560, 33633600, 168168000]
[6] [0, 1, 87380, 42771456, 2385102720, 34306272000, 137225088000]

[1, 0, 1, 0, 1, 10, 0, 1, 84, 280, 0, 1, 682, 9240, 15400, 0, 1, 5460, 260260, 1401400, 1401400]

RowSum      A291973 [1, 1, 11, 365, 25323, 3068521, 583027547]
AlterRowSum A291974 [1, -1, 9, -197, 6841, -254801, -3000807]
Diagonal    A025035 [1, 1, 10, 280, 15400, 1401400, 190590400]
Central     A327417 [1, 1, 682, 7128576, 429120851544, 94066556834970720]

('--', 4, '---------------------')

by Shape    A327004
[0] [1]
[1] [1]
[2] [1, 35]
[3] [1, 495, 5775]
[4] [1, 1820, 6435, 450450, 2627625]
[5] [1, 4845, 125970, 4408950, 31177575, 727476750, 2546168625]
[6] [1, 10626, 735471, 1352078, 25741485, 1338557220, 1577585295, 15616500900, 165646455975, 1932541986375, 4509264634875]

[1, 1, 1, 35, 1, 495, 5775, 1, 1820, 6435, 450450, 2627625]

by Length   A291452
[0] [1]
[1] [0, 1]
[2] [0, 1, 35]
[3] [0, 1, 495, 5775]
[4] [0, 1, 8255, 450450, 2627625]
[5] [0, 1, 130815, 35586525, 727476750, 2546168625]
[6] [0, 1, 2098175, 2941884000, 181262956875, 1932541986375, 4509264634875]

by Recur/k!
[0] [1]
[1] [0, 1]
[2] [0, 1, 35]
[3] [0, 1, 495, 5775]
[4] [0, 1, 8255, 450450, 2627625]
[5] [0, 1, 130815, 35586525, 727476750, 2546168625]
[6] [0, 1, 2098175, 2941884000, 181262956875, 1932541986375, 4509264634875]

by Recurrence
[0] [1]
[1] [0, 1]
[2] [0, 1, 70]
[3] [0, 1, 990, 34650]
[4] [0, 1, 16510, 2702700, 63063000]
[5] [0, 1, 261630, 213519150, 17459442000, 305540235000]
[6] [0, 1, 4196350, 17651304000, 4350310965000, 231905038365000, 3246670537110000]

[1, 0, 1, 0, 1, 35, 0, 1, 495, 5775, 0, 1, 8255, 450450, 2627625, 0, 1, 130815, 35586525, 727476750, 2546168625]

RowSum      A291975 [1, 1, 36, 6271, 3086331, 3309362716, 6626013560301]
AlterRowSum A291976 [1, -1, 34, -5281, 2185429, -1854147586, 2755045819549]
Diagonal    A025036 [1, 1, 35, 5775, 2627625, 2546168625, 4509264634875]
Central     A327418 [1, 1, 8255, 2941884000, 11957867341948125, 294040106448733743008625]

=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
#Todo:
#plotting via timeseries
#PkgBenchmark
#HAC s.e.s from CovariancesMatrices.jl?
#how to export arch?
#Forecasting
#Float16/32 don't seem to work anymore. Problem in Optim?
#support missing data? timeseries?
#a simulated AM should probably contain a (zero) intercept, so that fit! is consistent with fit.
#the constructor for UnivariateARCHModel should make a copy of its args
#implement lrtest
#allow uninititalized constructors for UnivariateVolatilitySpec, MeanSpec and StandardizedDistribution? If so, then be consistent with how they are defined
#  (change for meanspec and dist ), document, and test. Also, NaN is prob. safer than undef.
#logconst needs to return the correct type
"""
The ARCHModels package for Julia. For documentation, see https://s-broda.github.io/ARCHModels.jl/dev.
"""
module ARCHModels
using Reexport
using Requires
@reexport using StatsBase
using StatsFuns: normcdf, normccdf, normlogpdf, norminvcdf, log2π, logtwo, RFunctions.tdistrand, RFunctions.tdistinvcdf, RFunctions.gammarand, RFunctions.gammainvcdf
using SpecialFunctions: beta, lgamma, gamma
using Optim
using ForwardDiff
using Distributions
using HypothesisTests
using Roots
using LinearAlgebra
using DataStructures: CircularBuffer
using DelimitedFiles
using Statistics: cov

import Distributions: quantile
import Base: show, showerror, eltype
import Statistics: mean
import Random: rand, AbstractRNG
import HypothesisTests: HypothesisTest, testname, population_param_of_interest, default_tail, show_params, pvalue
import StatsBase: StatisticalModel, stderror, loglikelihood, nobs, fit, fit!, confint, aic,
                  bic, aicc, dof, coef, coefnames, coeftable, CoefTable,
				  informationmatrix, islinear, score, vcov, residuals, predict
export ARCHModel, UnivariateARCHModel, UnivariateVolatilitySpec, StandardizedDistribution, Standardized, MeanSpec,
       simulate, simulate!, selectmodel, StdNormal, StdT, StdGED, Intercept, Regression,
       NoIntercept, ARMA, AR, MA, BG96, volatilities, mean, quantile, VaRs, pvalue, means, VolatilitySpec,
	   MultivariateVolatilitySpec, MultivariateStandardizedDistribution, MultivariateARCHModel, MultivariateStdNormal,
	   EGARCH, ARCH, GARCH, TGARCH, ARCHLMTest, DQTest,
	   DOW29, DCC, CCC, covariances, correlations

include("utils.jl")
include("general.jl")
include("univariatearchmodel.jl")
include("meanspecs.jl")
include("univariatestandardizeddistributions.jl")
include("EGARCH.jl")
include("TGARCH.jl")
include("tests.jl")
include("multivariatearchmodel.jl")
include("multivariatestandardizeddistributions.jl")
include("DCC.jl")
function __init__()
	@require GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" begin
		using .GLM
		import .StatsModels: TableRegressionModel
		function fit(vs::Type{VS}, lm::TableRegressionModel{<:LinearModel}; kwargs...) where VS<:UnivariateVolatilitySpec
			fit(vs, response(lm.model); meanspec=Regression(modelmatrix(lm.model); coefnames=coefnames(lm)), kwargs...)
		end
		function fit(vs::Type{VS}, lm::LinearModel; kwargs...) where VS<:UnivariateVolatilitySpec
			fit(vs, response(lm); meanspec=Regression(modelmatrix(lm)), kwargs...)
		end
	end
end
end#module


text:
import AlphaZero.GI
using AlphaZero: AbstractSchedule, ConstSchedule
using GraphNeuralNetworks
using Random
using Crayons
using JSON3
using StructTypes

o2ji(o, M, N) = CartesianIndices((M, N))[o]
ji2o(j, i, M, N) = LinearIndices((M, N))[j, i]

struct GameSpec <: GI.AbstractGameSpec
  M::Pair{AbstractSchedule,AbstractSchedule}
  N::Pair{AbstractSchedule,AbstractSchedule}
  A::Pair{AbstractSchedule,AbstractSchedule}
  K::Pair{AbstractSchedule,AbstractSchedule}
  P::Pair{UInt8,UInt8}
  T::Pair{UInt8,UInt8}
end

mutable struct AdaptiveNodes
  src::Matrix{UInt8}
  tar::Matrix{UInt8}
  info::Matrix{UInt8}
  done_time::Matrix{UInt16}
end

Base.copy(s::AdaptiveNodes) = AdaptiveNodes(copy(s.src), copy(s.tar), copy(s.info), copy(s.done_time))
Base.isequal(a::AdaptiveNodes, b::AdaptiveNodes) = isequal((a.src, a.tar, a.info, a.done_time), (b.src, b.tar, b.info, b.done_time))
function Base.hash(s::AdaptiveNodes, h::UInt)
  fs = (getfield(s, k) for k in fieldnames(AdaptiveNodes))
  return foldl((h, f) -> hash(f, h), fs, init=hash(AdaptiveNodes, h))
end

function generate_conjuctive_edges(num_operations, N_OPP, T, S)
  N = length(num_operations)
  target = [collect(UInt8, 2:N_OPP*2+1); repeat([T], N+1); zeros(UInt8, N)]
  index = 1
  for (i, num) in enumerate(num_operations)
    target[S+i-1] = index
    index += num * 2
    target[index-1] = T - i
  end
  return target
end

function gen_done_time(p_time, t_time, conj_tar, start_edges, S, T)
  done_time = zeros(UInt16, S)
  for o in start_edges
    last_done_time = 0
    last_m = size(t_time, 1)#M+2
    while (true)#propagate expected done time
      t = conj_tar[o]
      o = conj_tar[t]
      if (o == T)
        done_time[t] = last_done_time + t_time[last_m, end]
        done_time[T] = max(done_time[T], done_time[t])
        break
      end
      m_time, m = findmin(p_time[o÷2, :])
      done_time[t] = last_done_time + t_time[last_m, m]
      last_done_time = done_time[o] = done_time[t] + m_time
      last_m = m
    end
  end
  return done_time
end

function gen_action_values(p_time, t_time, conj_tar, start_edges, K, S)
  nodes = AdaptiveNodes(Matrix{UInt8}(undef, 0, 5), Matrix{UInt8}(undef, 0, 5), Matrix{UInt8}(undef, 0, 4), Matrix{UInt8}(undef, 0, 2))
  for (i, e) in enumerate(start_edges)
    t = conj_tar[e]
    o = conj_tar[t]
    for (m, m_time) in enumerate(p_time[o÷2, :])
      m_time == 0xff && continue
      for k in 1:K
        t_node_id = S + length(nodes.done_time) + 1
        m_node_id = t_node_id + 1
        nodes.src = [nodes.src; [S S t_node_id m_node_id S]]
        nodes.tar = [nodes.tar; [t_node_id t_node_id m_node_id conj_tar[o] m_node_id]]
        nodes.info = [nodes.info; [m k o i]]
        nodes.done_time = [nodes.done_time; [t_time[end, m] t_time[end, m] + m_time]]
      end
    end
  end
  return nodes
end

mutable struct GameEnv <: GI.AbstractGameEnv
  #Problem instance
  process_time::Matrix{UInt8} #Index in (Mxn)*M
  transport_time::Matrix{UInt8} #M+2xM+2
  M::UInt8
  N::UInt8
  K::UInt8
  N_OPP::UInt16
  T::UInt16
  S::UInt16
  UB::UInt16
  LB::UInt16
  conj_src::Vector{UInt8}
  conj_tar::Vector{UInt8}
  #State
  disj_src::Vector{UInt8}
  disj_tar::Vector{UInt8}
  is_done::Vector{Bool}
  done_time::Vector{UInt16}
  adaptive_nodes::AdaptiveNodes
  #Info
  prev_operation::Matrix{UInt8}
  prev_machine::Vector{UInt8}
  prev_vehicle::Matrix{UInt8}
  start_machine::Vector{UInt8}
  start_vehicle::Vector{UInt8}
end

function GI.init(spec::GameSpec, itc::Int, rng::AbstractRNG)
  N = rand(rng, spec.N.first[itc]:spec.N.second[itc])
  M = rand(rng, spec.M.first[itc]:spec.M.second[itc])
  K = rand(rng, spec.K.first[itc]:spec.K.second[itc])
  num_operations = [rand(rng, 1:M) for _ in 1:N]
  N_OPP = sum(num_operations)
  T = N_OPP * 2 + N + 1
  S = N_OPP * 2 + N + 2
  p_time = rand(rng, spec.P.first:spec.P.second, N_OPP, M)
  for o in 1:N_OPP
    ind = randperm(rng, M)[1:M-rand(rng, spec.A.first[itc]:spec.A.second[itc])]
    p_time[o, ind] .= 0xff
  end
  t_time = rand(rng, spec.T.first:spec.T.second, M + 1, M + 1)
  for m in 1:M+1
    t_time[m,m] = 0x00
  end
  conj_tar = generate_conjuctive_edges(num_operations, N_OPP, T, S) #fix
  start_edges_n = collect(0:N-1) .+ S
  node_done = gen_done_time(p_time, t_time, conj_tar, start_edges_n, S, T) #fix
  adaptive_nodes = gen_action_values(p_time, t_time, conj_tar, start_edges_n, K, S)
  return GameEnv(
    #Problem instance
    p_time,
    t_time,
    M,
    N,
    K,
    N_OPP,
    T,
    S,
    sum(maximum(replace(p_time, 0xff => 0x00), dims=2)) + maximum(t_time) * 2 * N_OPP, #All worst operations in a row,
    node_done[T],
    [collect(1:N_OPP*2 + N); T; S * ones(N)],
    conj_tar,
    #State
    [],
    [],
    [falses(N_OPP * 2 + N); false; true], #[Transport, Operations, returns, T, S]
    node_done,
    adaptive_nodes,
    #Info
    repeat([S (M + 1)], N),
    repeat([S], M),
    repeat([S (M + 1)], K),
    zeros(UInt8, M),
    zeros(UInt8, K)
  )
end

struct FJSPTInstance
  num_vehicles::Int
  num_machines::Int
  num_operations::Vector{Int}
  process_time::Vector{Int}
  transport_time::Vector{Int}
end

StructTypes.StructType(::Type{FJSPTInstance}) = StructTypes.Struct()

function GI.init(spec::GameSpec, instance_string::String)
  instance = JSON3.read(instance_string, FJSPTInstance)
  p_time = reshape(instance.process_time, :, instance.num_machines)
  t_time = reshape(instance.transport_time, :, instance.num_machines + 1)
  t_time = [t_time[2:end, 2:end] t_time[2:end, 1]; t_time[1, 2:end]... 0]#reshape so LU is last
  N = length(instance.num_operations)
  M = instance.num_machines
  K = instance.num_vehicles
  N_OPP = sum(instance.num_operations)
  T = N_OPP * 2 + N + 1
  S = N_OPP * 2 + N + 2
  conj_tar = generate_conjuctive_edges(instance.num_operations, N_OPP, T, S)
  start_edges_n = collect(0:N-1) .+ S
  node_done = gen_done_time(p_time, t_time, conj_tar, start_edges_n, S, T)
  adaptive_nodes = gen_action_values(p_time, t_time, conj_tar, start_edges_n, K, S)
  return GameEnv(
    p_time,
    t_time,
    instance.num_machines,
    N,
    instance.num_vehicles,
    N_OPP,
    T,
    S,
    sum(maximum(replace(p_time, 0xff => 0x00), dims=2)) + maximum(t_time) * 2 * N_OPP, #All worst operations in a row,
    node_done[T],
    [collect(1:N_OPP*2 + N); T; S * ones(N)],
    conj_tar,
    #State
    [],
    [],
    [falses(N_OPP * 2 + N); false; true],
    node_done,
    adaptive_nodes,
    #Info
    repeat([S (M + 1)], N),
    repeat([S], M),
    repeat([S (M + 1)], K),
    zeros(UInt8, M),
    zeros(UInt8, K)
  )
end

GI.init(spec::GameSpec, s) = GameEnv(
  #Static values
  s.process_time,
  s.transport_time,
  s.M,
  s.N,
  s.K,
  s.N_OPP,
  s.T,
  s.S,
  s.UB,
  s.LB,
  s.conj_src,
  s.conj_tar,
  #Mutable values
  copy(s.disj_src),
  copy(s.disj_tar),
  copy(s.is_done),
  copy(s.done_time),
  copy(s.adaptive_nodes),
  copy(s.prev_operation),
  copy(s.prev_machine),
  copy(s.prev_vehicle),
  copy(s.start_machine),
  copy(s.start_vehicle)
)

GI.spec(g::GameEnv) = GameSpec(ConstSchedule(g.M) => ConstSchedule(g.M), ConstSchedule(g.N) => ConstSchedule(g.N), ConstSchedule(2) => ConstSchedule(2), ConstSchedule(2) => ConstSchedule(2), 1 => 5, 1 => 5)

GI.two_players(::GameSpec) = false

GI.state_dim(spec::GameSpec) = (2, (spec.M.second[1] * spec.N.second[1] + 2))#Opperations + source and sink

function GI.set_state!(g::GameEnv, s)
  g.process_time = s.process_time
  g.transport_time = s.transport_time
  g.M = s.M
  g.N = s.N
  g.K = s.K
  g.N_OPP = s.N_OPP
  g.T = s.T
  g.S = s.S
  g.UB = s.UB
  g.LB = s.LB
  g.conj_src = s.conj_src
  g.conj_tar = s.conj_tar
  g.disj_src = copy(s.disj_src)
  g.disj_tar = copy(s.disj_tar)
  g.is_done = copy(s.is_done)
  g.done_time = copy(s.done_time)
  g.adaptive_nodes = copy(s.adaptive_nodes)
  g.prev_operation = copy(s.prev_operation)
  g.prev_machine = copy(s.prev_machine)
  g.prev_vehicle = copy(s.prev_vehicle)
  g.start_machine = copy(s.start_machine)
  g.start_vehicle = copy(s.start_vehicle)
end

GI.current_state(g::GameEnv) = (
  process_time=g.process_time,
  transport_time=g.transport_time,
  M=g.M,
  N=g.N,
  K=g.K,
  N_OPP=g.N_OPP,
  T=g.T,
  S=g.S,
  UB=g.UB,
  LB=g.LB,
  conj_src=g.conj_src,
  conj_tar=g.conj_tar,
  disj_src=copy(g.disj_src),
  disj_tar=copy(g.disj_tar),
  is_done=copy(g.is_done),
  done_time=copy(g.done_time),
  adaptive_nodes=copy(g.adaptive_nodes),
  prev_operation=copy(g.prev_operation),
  prev_machine=copy(g.prev_machine),
  prev_vehicle=copy(g.prev_vehicle),
  start_machine=copy(g.start_machine),
  start_vehicle=copy(g.start_vehicle)
)

GI.white_playing(g::GameEnv) = true

GI.game_terminated(g::GameEnv) = g.is_done[g.T]

GI.actions(spec::GameSpec) = collect(1:(spec.M.second[1]*spec.N.second[1]+2))

GI.available_actions(g::GameEnv) = collect(g.S+2:2:g.S+length(g.adaptive_nodes.done_time))

GI.actions_mask(g::GameEnv) = [falses(g.S); repeat([false, true], size(g.adaptive_nodes.done_time, 1))]#fix

function GI.play!(g::GameEnv, action)
  #mark operation scheduled
  m, k, o, i = g.adaptive_nodes.info[(action-g.S)÷2, :]
  t = o - 1
  @assert g.is_done[o] == false
  g.is_done[o] = true
  g.is_done[t] = true
  #update previous operation and machine
  p = g.prev_machine[m] #previous operation done on machine of todo operation
  l, m_l = g.prev_operation[i, :] #previous operation done in job of todo operation
  u, m_u = g.prev_vehicle[k, :]
  #update info vectors
  g.prev_machine[m] = o
  g.prev_operation[i, :] = [o, m]
  isequal(m_l, m) || (g.prev_vehicle[k, :] = [t, m])#don't update if no transport needed
  iszero(g.start_machine[m]) && (g.start_machine[m] = o)
  iszero(g.start_vehicle[k]) && (g.start_vehicle[k] = t)
  #add disjunctive link
  if(!isequal(m_l, m))
    append!(g.disj_src, u)
    append!(g.disj_tar, t)
  end
  append!(g.disj_src, p)
  append!(g.disj_tar, o)
  #set done time
  g.done_time[t] = max(g.done_time[p], max(g.done_time[l], g.done_time[u] + g.transport_time[m_u, m_l]) + g.transport_time[m_l, m])
  g.done_time[o] = max(g.done_time[t], g.done_time[p]) + g.process_time[o÷2, m]
  #remove old actions
  mask = g.adaptive_nodes.info[:, 3] .!== o
  g.adaptive_nodes.src = g.adaptive_nodes.src[mask, :]
  g.adaptive_nodes.tar = g.adaptive_nodes.tar[mask, :]
  g.adaptive_nodes.info = g.adaptive_nodes.info[mask, :]
  g.adaptive_nodes.done_time = g.adaptive_nodes.done_time[mask, :]
  #fix disjunctive edges
  g.adaptive_nodes.src[g.adaptive_nodes.info[:, 2].==k, 2] .= t
  g.adaptive_nodes.src[g.adaptive_nodes.info[:, 1].==m, 5] .= o
  #fix node ids
  num_adaptive_nodes = size(g.adaptive_nodes.done_time, 1) * 2
  t_ids = collect(g.S+1:2:g.S+num_adaptive_nodes)
  m_ids = t_ids .+ 1
  g.adaptive_nodes.tar[:, 1] = t_ids
  g.adaptive_nodes.tar[:, 2] = t_ids
  g.adaptive_nodes.src[:, 3] = t_ids
  g.adaptive_nodes.src[:, 4] = m_ids
  g.adaptive_nodes.tar[:, 3] = m_ids
  g.adaptive_nodes.tar[:, 5] = m_ids
  #add new actions
  next_t = g.conj_tar[o]
  next_o = g.conj_tar[next_t]
  next_next_t = g.conj_tar[next_o]
  if (!isequal(next_o, g.T))
    for (p_m, p_time) in enumerate(g.process_time[next_o÷2, :])
      p_time == 0xff && continue
      for k in 1:g.K
        t_node_id = g.S + length(g.adaptive_nodes.done_time) + 1
        m_node_id = t_node_id + 1
        g.adaptive_nodes.src = [g.adaptive_nodes.src; [o g.prev_vehicle[k] t_node_id m_node_id g.prev_machine[p_m]]]
        g.adaptive_nodes.tar = [g.adaptive_nodes.tar; [t_node_id t_node_id m_node_id next_next_t m_node_id]]
        g.adaptive_nodes.info = [g.adaptive_nodes.info; [p_m k next_o i]]
        transport_done = g.done_time[o] + g.transport_time[m, p_m]
        g.adaptive_nodes.done_time = [g.adaptive_nodes.done_time; [transport_done transport_done + p_time]]
      end
    end
  else
    g.is_done[next_t] = true
    append!(g.disj_src, t)
    append!(g.disj_tar, next_t)
    g.prev_vehicle[k, :] = [next_t, g.M+1]
  end
  #mark last node as done
  isempty(g.adaptive_nodes.done_time) && (g.is_done[g.T] = true)
  #propagate expected done time
  last_done_time = g.done_time[o]
  last_m = m
  while (true)
    t = g.conj_tar[o]
    o = g.conj_tar[t]
    if (o == g.T)
      g.done_time[t] = last_done_time + g.transport_time[last_m, end]
      g.done_time[o] = max(g.done_time[o], g.done_time[t])
      return
    end
    m_time, m = findmin(g.process_time[o÷2, :])
    g.done_time[t] = last_done_time + g.transport_time[last_m, m]
    last_done_time = g.done_time[o] = g.done_time[t] + m_time
    last_m = m
  end
end

function GI.white_reward(g::GameEnv)
  if (all(g.is_done))
    #return -convert(Float32, g.done_time[g.T])
    return ((g.UB - g.done_time[g.T]) / (g.UB - g.LB))
  end
  return 0
end

function GI.vectorize_state(::GameSpec, state)
  num_actions = length(state.adaptive_nodes.done_time)
  return GNNGraph([state.conj_src; state.disj_src; vec(state.adaptive_nodes.src)],
    [state.conj_tar; state.disj_tar; vec(state.adaptive_nodes.tar)],
    num_nodes=state.S + num_actions,
    ndata=Float32.([[state.done_time; vec(state.adaptive_nodes.done_time)] [state.is_done; zeros(num_actions)]]'))
end

#####
##### Interaction APIc
#####

function GI.action_string(spec::GameSpec, o)
  mn = o2ji(o, spec.M.second, spec.N.second)
  return string("job: ", mn[2], " machine: ", mn[1])
end

function GI.parse_action(spec::GameSpec, str)
  try
    s = split(str)
    @assert length(s) == 2
    n = parse(Int, s[1])
    m = parse(Int, s[2])
    return ji2o(m, n, spec.M.second, spec.N.second)
  catch e
    return nothing
  end
end

function GI.read_state(::GameSpec)#input problem
  return nothing
end

function GI.render(g::GameEnv)
  m_opperation = zeros(UInt8, g.N_OPP)
  for (m, o) in enumerate(g.start_machine)
    print(crayon"white", "m", m, ":")
    last_done = 0
    while (!iszero(o))
      i = count(>(g.N_OPP*2), g.conj_tar[1:o-1])+1
      idle = (crayon"dark_gray", repeat("-", g.done_time[o] - g.process_time[o÷2, m] - last_done))
      active = (Crayon(foreground=i), repeat("=", g.process_time[o÷2, m]))
      print(idle..., active...)
      last_done = g.done_time[o]
      m_opperation[o÷2] = m
      #get next operation for machine
      index = findfirst(isequal(o), g.disj_src)
      if(isnothing(index))
        break
      end
      o = g.disj_tar[index]
    end
    fill = (crayon"dark_gray", repeat("-", g.done_time[g.T] - last_done))
    print(fill...)
    println(crayon"reset")
  end
  for (k, t) in enumerate(g.start_vehicle)
    print(crayon"white", "k", k, ":")
    last_done = 0
    location = g.M+1
    while (!iszero(t))
      if(t<=g.N_OPP*2)
        i = count(>(g.N_OPP*2), g.conj_tar[1:t])+1
        dropoff = m_opperation[(t+1)÷2]
      else
        i = Int64(g.T - t)
        dropoff = g.M + 1
      end
      pickup_index = findfirst(isequal(t), g.conj_tar)
      pickup = pickup_index>g.N_OPP*2 ? g.M+1 : m_opperation[pickup_index÷2]
      empty_time = g.transport_time[location, pickup]
      loaded_time = g.transport_time[pickup, dropoff]
      idle = (crayon"dark_gray", repeat("-", g.done_time[t] - last_done - empty_time  - loaded_time))
      empty = (crayon"dark_gray", repeat("=", empty_time))
      loaded = (Crayon(foreground=i), repeat("=", loaded_time))
      print(idle..., empty..., loaded...)
      last_done = g.done_time[t]
      location = dropoff
      #get next transport for vehicle
      index = findfirst(isequal(t), g.disj_src)
      isnothing(index) && break
      t = g.disj_tar[index]
    end
    fill = (crayon"dark_gray", repeat("-", g.done_time[g.T] - last_done))
    print(fill...)
    println(crayon"reset")
  end
  for n in 1:g.N
    print(Crayon(foreground=n), "n", n, " ")
  end
  print(crayon"white", "makespan: ", g.done_time[g.T])
  println(crayon"reset")
end


text:
"""
    WhiteNoise(amplitude)
Defining WhiteNoise struct, values generated by it are all between -amplitude and +amplitude

# Example
```julia.repl
julia> w = WhiteNoise(1)
WhiteNoise(1)
julia> w(1)
```
"""
struct WhiteNoise

    amplitude::Number

end

"""
    GaussianNoise(amplitude)
Defining GaussianNoise struct, values are generated by normal distribution function randn()

# Example
```julia.repl
julia> w = WhiteNoise(1)
WhiteNoise(1)
julia> w(1)
```
"""
struct GaussianNoise

    amplitude::Number

end

struct BrownianNoise

    amplitude::Number

end

struct PinkNoise

    amplitude::Number
    beta::Number

end

function (wn::WhiteNoise)(t, fr)

    time = 0:inv(fr):t
    amps = zeros(length(time)) + wn

    return Wave(
        amps,
        0:inv(fr):t,
        fr,
        Signal(fr, wn.amplitude, x->0, 0, Dict("name" => "Noise")),
    )

end

function (gn::GaussianNoise)(t, fr)

    time = 0:inv(fr):t
    amps = zeros(length(time)) + gn

    return Wave(
        amps,
        0:inv(fr):t,
        fr,
        Signal(fr, gn.amplitude, x -> 0, 0, Dict("name" => "Noise")),
    )

end

function (pn::PinkNoise)(t, fr)

    time = 0:inv(fr):t
    amps = zeros(length(time)) + pn

    return Wave(
        amps,
        0:inv(fr):t,
        fr,
        Signal(fr, pn.amplitude, x->0,0, Dict("name" => "Noise")),
    )

end

# TO Do check type interference!
function +(v, wn::WhiteNoise)

    a =  wn.amplitude

    r = rand(-a:eps(Float64):a,size(v))

    return v + r

end

function +(v, gn::GaussianNoise)

    a =  gn.amplitude

    r = randn(size(v)) .* a

    return v + r

end

function +(v, pn::PinkNoise)

    s = length(v)

    wn = WhiteNoise(pn.amplitude)(1 - inv(s),s)
    swn = wn |> Spectrum

    denominator = swn.fs .^(pn.beta/2)

    denominator[1] = 1

    pink_amps = swn.amps ./ denominator

    spec = Spectrum(pink_amps,inv(s),abs2.(pink_amps), Signal(100,1),wn.framerate ,length(wn.amps))

    w = Wave(spec)

    amps = w.amps

    amps = amps .- mean(amps)

    amps = amps .* (pn.amplitude / maximum(abs,amps))

    return v + amps

end

function +(v, bn::BrownianNoise)

    a = bn.amplitude

    r = rand(-a:eps(Float64):a,size(v))

    cumr = cumsum(r)

    cumr = cumr .- mean(cumr) # make signal zero mean

    cumr = cumr * bn.amplitude/maximum(abs,cumr) # normalize

    return v + cumr

end

function (bn::BrownianNoise)(t, fr)

    time = 0:inv(fr):t
    n = length(time)
    v = zeros(n)

    brownian_amps = v + bn

    return Wave(
        brownian_amps,
        time,
        fr,
        Signal(fr, bn.amplitude,x -> 0, 0, Dict("name" => "Noise")),
    )

end

function +(s::Signal, wn::WhiteNoise)

    s_new = deepcopy(s)

    s_new.info["WhiteNoise"] = wn.amplitude

    if haskey(s_new.info, "name")

        s_new.info["name"] *= " with White Noise"

    end

    return s_new

end

function +(s::Signal, gn::GaussianNoise)

    s_new = deepcopy(s)

    s_new.info["GaussianNoise"] = gn.amplitude

    if haskey(s_new.info, "name")

        s_new.info["name"] *= " with Gaussian Noise"

    end

    return s_new

end

function +(s::Signal, bn::BrownianNoise)

    s_new = deepcopy(s)

    s_new.info["BrownianNoise"] = bn.amplitude

    if haskey(s_new.info, "name")

        s_new.info["name"] *= " with Brownian Noise"

    end

    return s_new

end

function +(s::Signal, pn::PinkNoise)

    s_new = deepcopy(s)

    s_new.info["PinkNoise"] = (pn.amplitude, pn.beta)

    if haskey(s_new.info, "name")

        s_new.info["name"] *= " with Pink Noise"

    end

    return s_new

end