text:
"The average reflection coefficients"
function reflection_coefficients(ωs::Union{T,AbstractVector{T}}, medium::Medium{T}, species::Vector{Specie{T}}; kws...) where T<:Number

    Rs = map(ωs) do ω
        k_effs = wavenumbers(ω, medium, species; kws...)
        w_effs = [EffectiveWave(ω, k_eff, medium, species; kws...) for k_eff in k_effs]
        reflection_coefficient(ω, w_effs, medium, species; kws...)
    end

    return Rs
end

"Calculates the reflection coefficient from each wave in waves and then sums the results."
reflection_coefficient(ω::T, waves::Vector{EffectiveWave{T}}, medium::Medium{T}, species::Vector{Specie{T}}; kws...) where T<:Number =
    sum(w -> reflection_coefficient(ω, w, medium, species; kws...), waves)

"Pairs each ω in ωs with each wave in waves to calculate the refleciton coefficients: reflection_coefficient(ω, wave)"
reflection_coefficients(ωs::AbstractVector{T}, waves::Vector{EffectiveWave{T}}, medium::Medium{T}, species::Vector{Specie{T}}; kws...) where T<:Number =
[ reflection_coefficient(ωs[i], waves[i], medium, species; kws...) for i in eachindex(ωs)]

"The average reflection coefficient"
function reflection_coefficient(ω::T, wave_eff::EffectiveWave{T}, medium::Medium{T}, species::Vector{Specie{T}};
        θin::T = zero(T), x::T = zero(T), kws...) where T<:Number

    k = ω/medium.c
    θ_ref = pi - wave_eff.θ_eff - θin
    S = length(species)
    ho = wave_eff.hankel_order

    kθ = (k*cos(θin) + wave_eff.k_eff*cos(wave_eff.θ_eff))
    R = 2.0im / (k*cos(θin) * kθ)
    R = R*sum(
        exp(im*n*θ_ref + im*x*kθ) * species[l].num_density *
        wave_eff.amplitudes[n+ho+1,l]
    for n=-ho:ho, l=1:S)

    return R
end

"The average reflection coefficient"
function wienerhopf_reflection_coefficient(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        tol::T = T(1e-7),
        θin::T = zero(T),
        radius_multiplier::T = 1.005,
        hankel_order::Int = 0,
        num_coefs::Int = 20000
    ) where T<:Number

    k = ω/medium.c
    ho = hankel_order

    t_vecs = t_vectors(ω, medium, species; hankel_order = ho)
    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]

    sToS(s,j::Int,l::Int) = (real(s) >= 0) ? sqrt(s^2 + (k*as[j,l]*sin(θin))^2) : -sqrt(s^2 + (k*as[j,l]*sin(θin))^2)

    function Ψ(s,j,l,m,n)
        (s^T(2) - (k*as[j,l]*cos(θin))^T(2)) * (n == m ? T(1) : T(0)) * (j == l ? T(1) : T(0)) +
        T(2) * as[j,l]^T(2) * pi * species[l].num_density * t_vecs[l][m+ho+1] *
        kernelN(n-m,k*as[j,l], sToS(s,j,l))
    end

    # kernelN(0,k*a12,Z) = k*a12*diffhankelh1(0,k*a12)*besselj(0,Z) - Z*hankelh1(0,k*a12)*diffbesselj(0,Z)
    q(s,j,l,m,n) = Ψ(s,j,l,m,n) / (s^T(2) - (k*as[j,l]*cos(θin))^T(2))

    Zs = LinRange(T(100),1/(10*tol),3000)
    maxZ = Zs[findfirst(Z -> abs(log(q(Z,1,1,0,0))) < 10*tol, Zs)]

    function Ψp(s, maxZ::T = maxZ, num_coefs::Int = num_coefs)
        Q(z) = log(q(z,1,1,0,0))/(z - s)
        xp = as[1,1]*k*cos(θin)*(-1.0+1.0im)
        q_pos = exp(
            (T(1.0)/(T(2)*pi*im)) * (
                sum(Fun(Q, Segment(-maxZ,xp), num_coefs)) +
                sum(Fun(Q, Segment(xp,-xp), num_coefs)) +
                sum(Fun(Q, Segment(-xp,maxZ), num_coefs))
            )
        )
        return (s + k*as[1,1]*cos(θin)) * q_pos
    end

    function Ψm(s, maxZ::T = maxZ, num_coefs::Int = num_coefs)
        Q(z) = log(q(z,1,1,0,0))/(z - s)
        xm = as[1,1]*k*cos(θin)*(-1.0+0.5im)
        a1 = T(0)
        q_neg = exp(
            -(T(1.0)/(T(2)*pi*im)) * (
                sum(Fun(Q, Segment(-maxZ,xm+a1), num_coefs)) +
                sum(Fun(Q, Segment(xm+a1,-xm+a1), num_coefs)) +
                sum(Fun(Q, Segment(-xm+a1, maxZ), num_coefs))
            )
        )
        return (s - k*as[1,1]*cos(θin)) * q_neg
    end

    x = -as[1,1]*k*cos(θin)*(-1.0+0.75im)

    # abs(Fp(x,maxZ,num_coefs) - Fp(x,maxZ, Int(round(num_coefs*1.1)))) / abs(Fp(x,maxZ,num_coefs))
    # abs(Fm(x,maxZ,num_coefs) - Fm(x,maxZ, Int(round(num_coefs*1.1)))) / abs(Fm(x,maxZ,num_coefs))

    x2 = as[1,1]*k*cos(θin)*(1.0+2.75im)

    err = abs(Ψp(x,maxZ,num_coefs) * Ψm(x,maxZ,num_coefs) - Ψ(x,1,1,0,0))/abs(Ψ(x,1,1,0,0))
    if err > tol
        @warn "Analytic split recovers original function with $err tolerance, instead of the specified tolernace: $tol"
    end

    R = Ψ(k*as[1,1]*cos(θin),1,1,0,0) / (Ψp(k*as[1,1]*cos(θin),maxZ,num_coefs))^2

    return R
end

# function F0(S,j,l,m,n)
#     (S^T(2) - (k*as[j,l])^T(2)) * (n == m ? T(1) : T(0)) * (j == l ? T(1) : T(0)) +
#     T(2) * as[j,l]^T(2) * pi*species[l].num_density*t_vecs[l][m+ho+1] * kernelN(n-m,k*as[j,l],S)
# end

# kernelN(0,k*a12,Z) = k*a12*diffhankelh1(0,k*a12)*besselj(0,Z) - Z*hankelh1(0,k*a12)*diffbesselj(0,Z)
# Q0(S,j,l,m,n) = F0(S,j,l,m,n) / (S^T(2) - (k*as[j,l])^T(2))

# function F0p(S, maxZ::T = maxZ, num_coefs::Int = num_coefs)
#     Q(Z) = log(Q0(Z,1,1,0,0))/(Z - S)
#     xp = as[1,1]*k*(-1.0+1.0im)
#     (S + k*as[1,1]) * exp(
#         (T(1.0)/(T(2)*pi*im)) * (
#             sum(Fun(Q, Segment(-maxZ, xp), num_coefs)) +
#             sum(Fun(Q, Segment(xp,-xp), num_coefs)) +
#             sum(Fun(Q, Segment(-xp,maxZ), num_coefs))
#         )
#     )
# end


text:
@testset "NumField/NfRel/Conjugates.jl" begin
  Qx, x = QQ["x"]
  K, a = NumberField(x^2 - 2, "a")
  Kt, t = K["t"]
  g = t^3 + a*t + 5
  L, b = NumberField(g, "b")
  r, s = @inferred signature(L)
  @test (r, s) == (2, 2)
  plc = @inferred infinite_places(L)
  @test count(isreal, plc) == 2
  @test count(!isreal, plc) == 2

  Labs, LabstoL, = absolute_field(L)
  c = conjugates(b, 64)
  cc = conjugates(LabstoL\b, 64)

  for _c in c
    @test count(_cc -> overlaps(_cc, _c), cc) == count(_cc -> overlaps(_cc, _c), c)
  end

  Ly, y = L["y"]
  M, m = number_field(y^2 + y + b - 5, "c")
  r, s = @inferred signature(M)
  @test (r, s) == (4, 4)

  plc = @inferred infinite_places(M)
  @test count(isreal, plc) == 4
  @test count(!isreal, plc) == 4

  Mabs, MabstoM, = absolute_field(M)
  Mabsabs, MabsabstoMabs = absolute_field(Mabs)

  o = rand(M, -1:1)

  c = conjugates(o, 64)
  cc = conjugates(MabsabstoMabs\(MabstoM\o), 64)

  for _c in c
    @test count(_cc -> overlaps(_cc, _c), cc) == count(_cc -> overlaps(_cc, _c), c)
  end
end

ext: jl


text:
module CommunalHelperChainedKevin

using ..Ahorn, Maple
using Ahorn.Cairo

@mapdef Entity "CommunalHelper/ChainedKevin" ChainedKevin(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    chillout::Bool=false,
    chainLength::Integer=64,
    direction::String="Right",
    chainOutline::Bool=true,
    centeredChain::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Kevin (Chained, $direction) (Communal Helper)" => Ahorn.EntityPlacement(
        ChainedKevin,
        "rectangle",
        Dict{String, Any}(
            "direction" => direction,
        ),
    ) for direction in Maple.move_block_directions
)

const smallFace = "objects/crushblock/idle_face"
const giantFace = "objects/crushblock/giant_block00"

const ChainedKevinColor = (98, 34, 43) ./ 255

Ahorn.editingOptions(entity::ChainedKevin) = Dict{String, Any}(
    "direction" => Maple.move_block_directions,
)

Ahorn.minimumSize(entity::ChainedKevin) = 24, 24
Ahorn.resizable(entity::ChainedKevin) = true, true

Ahorn.selection(entity::ChainedKevin) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ChainedKevin, room::Maple.Room)
    direction = String(get(entity.data, "direction", "Right"))
    chillout = get(entity.data, "chillout", false)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    face = height >= 48 && width >= 48 && chillout ? giantFace : smallFace
    frame = "objects/CommunalHelper/chainedKevin/block" * direction
    faceSprite = Ahorn.getSprite(face, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    chainLength = Int(get(entity.data, "chainLength", 64))
    fillColor = (0.0, 0.0, 0.0, 0.0)
    outlineColor = (1.0, 1.0, 1.0, 0.5)

    save(ctx)

    set_dash(ctx, [0.6, 0.2])
    set_antialias(ctx, 1)
    set_line_width(ctx, 1)

    if direction == "Up"
        Ahorn.drawRectangle(ctx, 1, -chainLength + 1, width - 1, height + chainLength - 1, fillColor, outlineColor)
    elseif direction == "Right"
        Ahorn.drawRectangle(ctx, 1, 1, width + chainLength - 1, height - 1, fillColor, outlineColor)
    elseif direction == "Left"
        Ahorn.drawRectangle(ctx, -chainLength + 1, 1, width + chainLength - 1, height - 1, fillColor, outlineColor)
    else
        Ahorn.drawRectangle(ctx, 1, 1, width - 1, height + chainLength - 1, fillColor, outlineColor)
    end

    restore(ctx)

    Ahorn.drawRectangle(ctx, 2, 2, width - 4, height - 4, ChainedKevinColor)
    Ahorn.drawImage(ctx, faceSprite, div(width - faceSprite.width, 2), div(height - faceSprite.height, 2))

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 24, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 24, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 24, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 24, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 24, 24, 8, 8)
end

end


text:
int_rules_8_8 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.8*Polylogarithm*function =#
    #= Antiderivative(PolyLog(2,(~a')*((~b')*(~x)^(~p'))^(~q')),~x) :=  x*PolyLog(2,a*(b*x^p)^q) + p*q*Antiderivative(log(1-a*(b*x^p)^q),x) <--  FreeQ([a,b,p,q],x) =#
    @apply_utils Antiderivative(PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')), ~x) => ~x * PolyLog(~n, ~a * (~b * (~x) ^ ~p) ^ ~q) - ~p * ~q * Antiderivative(PolyLog(~n - 1, ~a * (~b * (~x) ^ ~p) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~p, ~q], ~x) && GtQ(~n, 0)
    @apply_utils Antiderivative(PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')), ~x) => (~x * PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q)) / (~p * ~q) - (1 / (~p * ~q)) * Antiderivative(PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~p, ~q], ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative(PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')), ~x) => Unintegrable(PolyLog(~n, ~a * (~b * (~x) ^ ~p) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~n, ~p, ~q], ~x)
    @apply_utils Antiderivative(PolyLog(~n, ~(c') * (~(a') + ~(b') * ~x) ^ ~(p')) / (~(d') + ~(e') * ~x), ~x) => PolyLog(~n + 1, ~c * (~a + ~b * ~x) ^ ~p) / (~e * ~p) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x) && EqQ(~b * ~d, ~a * ~e)
    @apply_utils Antiderivative(PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')) / ~x, ~x) => PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q) / (~p * ~q) <-- FreeQ([~a, ~b, ~n, ~p, ~q], ~x)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')), ~x) => ((~d * ~x) ^ (~m + 1) * PolyLog(~n, ~a * (~b * (~x) ^ ~p) ^ ~q)) / (~d * (~m + 1)) - ((~p * ~q) / (~m + 1)) * Antiderivative((~d * ~x) ^ ~m * PolyLog(~n - 1, ~a * (~b * (~x) ^ ~p) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~d, ~m, ~p, ~q], ~x) && (NeQ(~m, -1) && GtQ(~n, 0))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')), ~x) => ((~d * ~x) ^ (~m + 1) * PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q)) / (~d * ~p * ~q) - ((~m + 1) / (~p * ~q)) * Antiderivative((~d * ~x) ^ ~m * PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~d, ~m, ~p, ~q], ~x) && (NeQ(~m, -1) && LtQ(~n, -1))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q')), ~x) => Unintegrable((~d * ~x) ^ ~m * PolyLog(~n, ~a * (~b * (~x) ^ ~p) ^ ~q), ~x) <-- FreeQ([~a, ~b, ~d, ~m, ~n, ~p, ~q], ~x)
    @apply_utils Antiderivative((log(~(c') * (~x) ^ ~(m')) ^ ~(r') * PolyLog(~n, ~(a') * (~(b') * (~x) ^ ~(p')) ^ ~(q'))) / ~x, ~x) => (log(~c * (~x) ^ ~m) ^ ~r * PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q)) / (~p * ~q) - ((~m * ~r) / (~p * ~q)) * Antiderivative((log(~c * (~x) ^ ~m) ^ (~r - 1) * PolyLog(~n + 1, ~a * (~b * (~x) ^ ~p) ^ ~q)) / ~x, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~q, ~r], ~x) && GtQ(~r, 0)
    @apply_utils Antiderivative(PolyLog(~n, ~(c') * (~(a') + ~(b') * ~x) ^ ~(p')), ~x) => (~x * PolyLog(~n, ~c * (~a + ~b * ~x) ^ ~p) - ~p * Antiderivative(PolyLog(~n - 1, ~c * (~a + ~b * ~x) ^ ~p), ~x)) + ~a * ~p * Antiderivative(PolyLog(~n - 1, ~c * (~a + ~b * ~x) ^ ~p) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && GtQ(~n, 0)
    @apply_utils Antiderivative(PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)) / (~(d') + ~(e') * ~x), ~x) => (log((1 - ~a * ~c) - ~b * ~c * ~x) * PolyLog(2, ~c * (~a + ~b * ~x))) / ~e + (~b / ~e) * Antiderivative(log((1 - ~a * ~c) - ~b * ~c * ~x) ^ 2 / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(~c * (~b * ~d - ~a * ~e) + ~e, 0)
    @apply_utils Antiderivative(PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)) / (~(d') + ~(e') * ~x), ~x) => (log(~d + ~e * ~x) * PolyLog(2, ~c * (~a + ~b * ~x))) / ~e + (~b / ~e) * Antiderivative((log(~d + ~e * ~x) * log((1 - ~a * ~c) - ~b * ~c * ~x)) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(~c * (~b * ~d - ~a * ~e) + ~e, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * PolyLog(2, ~c * (~a + ~b * ~x))) / (~e * (~m + 1)) + (~b / (~e * (~m + 1))) * Antiderivative(((~d + ~e * ~x) ^ (~m + 1) * log((1 - ~a * ~c) - ~b * ~c * ~x)) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && NeQ(~m, -1)
    #= Antiderivative(((~d')+(~e')*(~x))^(~m')*PolyLog((~n),(~c')*((~a')+(~b')*(~x))^(~p')),~x) :=   (d+e*x)^(m+1)*PolyLog(n,c*(a+b*x)^p)/(e*(m+1)) - b*p/(e*(m+1))*Antiderivative((d+e*x)^(m+1)*PolyLog(n-1,c*(a+b*x)^p)/(a+b*x),x)  <-- FreeQ([a,b,c,d,e,m,p],x) && GtQ(n,0) && IGtQ(m,0) =#
    @apply_utils Antiderivative((~x) ^ ~(m') * PolyLog(~n, ~(c') * (~(a') + ~(b') * ~x) ^ ~(p')), ~x) => (-(((~a) ^ (~m + 1) - (~b) ^ (~m + 1) * (~x) ^ (~m + 1))) * PolyLog(~n, ~c * (~a + ~b * ~x) ^ ~p)) / ((~m + 1) * (~b) ^ (~m + 1)) + (~p / ((~m + 1) * (~b) ^ ~m)) * Antiderivative(ExpandIntegrand(PolyLog(~n - 1, ~c * (~a + ~b * ~x) ^ ~p), ((~a) ^ (~m + 1) - (~b) ^ (~m + 1) * (~x) ^ (~m + 1)) / (~a + ~b * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (GtQ(~n, 0) && (IntegerQ(~m) && NeQ(~m, -1)))
    @apply_utils Antiderivative((~(g') + ~(h') * log(~(f') * (~(d') + ~(e') * ~x) ^ ~(n'))) * PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)), ~x) => (~x * (~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * PolyLog(2, ~c * (~a + ~b * ~x)) + ~b * Antiderivative((~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * log((1 - ~a * ~c) - ~b * ~c * ~x) * ExpandIntegrand(~x / (~a + ~b * ~x), ~x), ~x)) - ~e * ~h * ~n * Antiderivative(PolyLog(2, ~c * (~a + ~b * ~x)) * ExpandIntegrand(~x / (~d + ~e * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n], ~x)
    @apply_utils Antiderivative((log(1 + ~(e') * ~x) * PolyLog(2, ~(c') * ~x)) / ~x, ~x) => -(PolyLog(2, ~c * ~x) ^ 2) / 2 <-- FreeQ([~c, ~e], ~x) && EqQ(~c + ~e, 0)
    @apply_utils Antiderivative(((~g + ~(h') * log(1 + ~(e') * ~x)) * PolyLog(2, ~(c') * ~x)) / ~x, ~x) => ~g * Antiderivative(PolyLog(2, ~c * ~x) / ~x, ~x) + ~h * Antiderivative((log(1 + ~e * ~x) * PolyLog(2, ~c * ~x)) / ~x, ~x) <-- FreeQ([~c, ~e, ~g, ~h], ~x) && EqQ(~c + ~e, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(g') + ~(h') * log(~(f') * (~(d') + ~(e') * ~x) ^ ~(n'))) * PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)), ~x) => (((~x) ^ (~m + 1) * (~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * PolyLog(2, ~c * (~a + ~b * ~x))) / (~m + 1) + (~b / (~m + 1)) * Antiderivative(ExpandIntegrand((~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * log((1 - ~a * ~c) - ~b * ~c * ~x), (~x) ^ (~m + 1) / (~a + ~b * ~x), ~x), ~x)) - ((~e * ~h * ~n) / (~m + 1)) * Antiderivative(ExpandIntegrand(PolyLog(2, ~c * (~a + ~b * ~x)), (~x) ^ (~m + 1) / (~d + ~e * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n], ~x) && (IntegerQ(~m) && NeQ(~m, -1))
    @apply_utils Antiderivative(~Px * (~(g') + ~(h') * log(~(f') * (~(d') + ~(e') * ~x) ^ ~(n'))) * PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)), ~x) => With([u = IntHide(~Px, ~x)], (u * (~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * PolyLog(2, ~c * (~a + ~b * ~x)) + ~b * Antiderivative(ExpandIntegrand((~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * log((1 - ~a * ~c) - ~b * ~c * ~x), u / (~a + ~b * ~x), ~x), ~x)) - ~e * ~h * ~n * Antiderivative(ExpandIntegrand(PolyLog(2, ~c * (~a + ~b * ~x)), u / (~d + ~e * ~x), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n], ~x) && PolyQ(~Px, ~x)
    @apply_utils Antiderivative((~x) ^ ~m * ~Px * (~(g') + ~(h') * log(1 + ~(e') * ~x)) * PolyLog(2, ~(c') * ~x), ~x) => Coeff(~Px, ~x, -(~m) - 1) * Antiderivative(((~g + ~h * log(1 + ~e * ~x)) * PolyLog(2, ~c * ~x)) / ~x, ~x) + Antiderivative((~x) ^ ~m * (~Px - Coeff(~Px, ~x, -(~m) - 1) * (~x) ^ (-(~m) - 1)) * (~g + ~h * log(1 + ~e * ~x)) * PolyLog(2, ~c * ~x), ~x) <-- FreeQ([~c, ~e, ~g, ~h], ~x) && (PolyQ(~Px, ~x) && (ILtQ(~m, 0) && (EqQ(~c + ~e, 0) && NeQ(Coeff(~Px, ~x, -(~m) - 1), 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Px * (~(g') + ~(h') * log(~(f') * (~(d') + ~(e') * ~x) ^ ~(n'))) * PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)), ~x) => With([u = IntHide((~x) ^ ~m * ~Px, ~x)], (u * (~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * PolyLog(2, ~c * (~a + ~b * ~x)) + ~b * Antiderivative(ExpandIntegrand((~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * log((1 - ~a * ~c) - ~b * ~c * ~x), u / (~a + ~b * ~x), ~x), ~x)) - ~e * ~h * ~n * Antiderivative(ExpandIntegrand(PolyLog(2, ~c * (~a + ~b * ~x)), u / (~d + ~e * ~x), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~n], ~x) && (PolyQ(~Px, ~x) && IntegerQ(~m))
    @apply_utils Antiderivative((~x) ^ ~m * ~(Px') * (~(g') + ~(h') * log(~(f') * (~(d') + ~(e') * ~x) ^ ~(n'))) * PolyLog(2, ~(c') * (~(a') + ~(b') * ~x)), ~x) => Unintegrable((~x) ^ ~m * ~Px * (~g + ~h * log(~f * (~d + ~e * ~x) ^ ~n)) * PolyLog(2, ~c * (~a + ~b * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && PolyQ(~Px, ~x)
    @apply_utils Antiderivative(PolyLog(~n, ~(d') * ((~F) ^ (~(c') * (~(a') + ~(b') * ~x))) ^ ~(p')), ~x) => PolyLog(~n + 1, ~d * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~p) / (~b * ~c * ~p * log(~F)) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * PolyLog(~n, ~(d') * ((~F) ^ (~(c') * (~(a') + ~(b') * ~x))) ^ ~(p')), ~x) => ((~e + ~f * ~x) ^ ~m * PolyLog(~n + 1, ~d * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~p)) / (~b * ~c * ~p * log(~F)) - ((~f * ~m) / (~b * ~c * ~p * log(~F))) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * PolyLog(~n + 1, ~d * ((~F) ^ (~c * (~a + ~b * ~x))) ^ ~p), ~x) <-- FreeQ([~F, ~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && GtQ(~m, 0)
    @apply_utils Antiderivative(~u * PolyLog(~n, ~v), ~x) => With([w = DerivativeDivides(~v, ~u * ~v, ~x)], w * PolyLog(~n + 1, ~v) <-- Not(FalseQ(w))) <-- FreeQ(~n, ~x)
    @apply_utils Antiderivative(~u * log(~w) * PolyLog(~n, ~v), ~x) => With([z = DerivativeDivides(~v, ~u * ~v, ~x)], z * log(~w) * PolyLog(~n + 1, ~v) - Antiderivative(SimplifyIntegrand((z * D(~w, ~x) * PolyLog(~n + 1, ~v)) / ~w, ~x), ~x) <-- Not(FalseQ(z))) <-- FreeQ(~n, ~x) && InverseFunctionFreeQ(~w, ~x)
end

max_stars_repo_path: src/rules/8 Special functions/8.8 Polylogarithm function.jl


text:
include("../src/tools.jl")

# test insert_char
testString = insert_char("aaaaa", 2, "A")
@assert(testString == "aAaaaa")

# test build_strings
testSequence = [ "AABB", "ACA", "CC" ]
testChromRep = [ [4], [4, 5], [1, 3, 4] ]
@assert(build_strings(testSequence, testChromRep) == [ "AAB_B", "ACA__", "_C__C" ])

#test gap_counts
gap_counts, len = get_gap_count(testSequence, 1.2)
# need 1, 2, and 3 gaps respectively to reach a length of 5 for each string
# 5 because that is the ceiling of 4 * 1.2, the largest string in the sequence
@assert(gap_counts == [1, 2, 3])
@assert(len == 5)


text:
abstract type AbstractLeakyESN <: AbstractEchoStateNetwork end

struct ESN{T, S<:AbstractArray{T}, I, B, F, N} <: AbstractLeakyESN
    res_size::I
    in_size::I
    out_size::I
    train_data::S
    alpha::T
    nla_type::N
    activation::F
    W::S
    W_in::S
    states::S
    extended_states::B
end

function ESN(approx_res_size::Int,
        train_data::AbstractArray{T},
        degree::Int,
        radius::T;
        activation::Any = tanh,
        sigma::T = 0.1,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = Int(floor(approx_res_size/in_size)*in_size)
    W = init_reservoir_givendeg(res_size, radius, degree)
    W_in = init_input_layer(res_size, in_size, sigma)
    states = states_matrix(W, W_in, train_data, alpha, activation, extended_states)

    return ESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, activation, W, W_in, states, extended_states)
end

#reservoir matrix W given by the user
function ESN(W::AbstractArray{T},
        train_data::Array{T};
        activation::Any = tanh,
        sigma::T = 0.1,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = size(W, 1)
    W_in = init_input_layer(res_size, in_size, sigma)
    states = states_matrix(W, W_in, train_data, alpha, activation, extended_states)

    return ESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, activation, W, W_in, states, extended_states)
end

#input layer W_in given by the user
function ESN(approx_res_size::Int,
        train_data::AbstractArray{T},
        degree::Int,
        radius::T,
        W_in::AbstractArray{T};
        activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1) #needs to be different?
    res_size = Int(floor(approx_res_size/in_size)*in_size)
    W = init_reservoir_givendeg(res_size, radius, degree)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch("size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch("size(W_in, 2) must be equal to in_size"))
    end

    states = states_matrix(W, W_in, train_data, alpha, activation, extended_states)

    return ESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, activation, W, W_in, states, extended_states)
end

#reservoir matrix W and input layer W_in given by the user
"""
    ESN(W::AbstractArray{T}, train_data::AbstractArray{T}, W_in::AbstractArray{T}
    [, activation::Any, alpha::T, nla_type::NonLinearAlgorithm, extended_states::Bool])

Build an ESN struct given the input and reservoir matrices.
"""
function ESN(W::AbstractArray{T},
        train_data::AbstractArray{T},
        W_in::AbstractArray{T};
        activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = size(W, 1)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch("size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch("size(W_in, 2) must be equal to in_size"))
    end

    states = states_matrix(W, W_in, train_data, alpha, activation, extended_states)

    return ESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, activation, W, W_in, states, extended_states)
end

function states_matrix(W::AbstractArray{Float64},
        W_in::AbstractArray{Float64},
        train_data::AbstractArray{Float64},
        alpha::Float64,
        activation::Function,
        extended_states::Bool)

    train_len = size(train_data, 2)
    res_size = size(W, 1)
    in_size = size(train_data, 1)

    states = zeros(Float64, res_size, train_len)
    for i=1:train_len-1
        states[:, i+1] = leaky_fixed_rnn(activation, alpha, W, W_in, states[:, i], train_data[:, i])

    end

    if extended_states == true
        ext_states = vcat(states, hcat(zeros(Float64, in_size), train_data[:,1:end-1]))
        return ext_states
    else
        return states
    end
end

"""
    ESNpredict(esn::AbstractLeakyESN, predict_len::Int, W_out::AbstractArray{Float64})

Return the prediction for a given length of the constructed ESN struct.
"""
function ESNpredict(esn::AbstractLeakyESN,
    predict_len::Int,
    W_out::AbstractArray{Float64})

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
        end
    elseif esn.extended_states == true
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end
    return output
end

"""
    ESNpredict_h_steps(esn::AbstractLeakyESN, predict_len::Int, h_steps::Int,
    test_data::AbstractArray{Float64}, W_out::AbstractArray{Float64})

Return the prediction for h steps ahead of the constructed ESN struct.
"""
function ESNpredict_h_steps(esn::AbstractLeakyESN,
    predict_len::Int,
    h_steps::Int,
    test_data::AbstractArray{Float64},
    W_out::AbstractArray{Float64})

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            if mod(i, h_steps) == 0
                x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, test_data[:,i])
            else
                x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
            end
        end
    elseif esn.extended_states == true
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            if mod(i, h_steps) == 0
                x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], test_data[:,i]), test_data[:,i])
            else
                x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
            end
        end
    end
    return output
end

"""
    ESNfitted(esn::AbstractLeakyESN, W_out::Matrix; autonomous=false)

Return the prediction for the training data using the trained output layer. The autonomous trigger can be used to have have it return an autonomous prediction starting from the first point if true, or a point by point prediction if false.
"""

function ESNfitted(esn::AbstractLeakyESN, W_out::Matrix; autonomous=false)
    train_len = size(esn.train_data, 2)
    output = zeros(Float64, esn.in_size, train_len)
    x = zeros(size(esn.states, 1))

    if autonomous
        out = esn.train_data[:,1]
        return _fitted!(output, esn, x, train_len, W_out, out)
    else
        return _fitted!(output, esn, x, train_len, W_out, esn.train_data)
    end
end

function _fitted!(output, esn, state, train_len, W_out, vector::Vector)
    if esn.extended_states == false
        for i=1:train_len
            state = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, state, vector)
            x_new = nla(esn.nla_type, state)
            vector = (W_out*x_new)
            output[:, i] = vector
        end
    elseif esn.extended_states == true
        for i=1:train_len
            state = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, state[1:esn.res_size], vector), vector)
            x_new = nla(esn.nla_type, state)
            vector = (W_out*x_new)
            output[:, i] = vector
        end
    end
    return output
end

function _fitted!(output, esn, state, train_len, W_out, vector::Matrix)
    if esn.extended_states == false
        for i=1:train_len
            state = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, state, vector[:,i])
            x_new = nla(esn.nla_type, state)
            out = (W_out*x_new)
            output[:, i] = out
        end
    elseif esn.extended_states == true
        for i=1:train_len
            state = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, state[1:esn.res_size], vector[:,i]), vector[:,i])
            x_new = nla(esn.nla_type, state)
            out = (W_out*x_new)
            output[:, i] = out
        end
    end
    return output
end


text:
module TestForgettingWeights

using Compat
import Distributions
using Test

using TreeParzen

# When run with a single trial, ApFilterTrials.ap_filter_trials produces an obs_above and
# an obs_below, one of which will contain a single value and the other of which will contain
# nothing. However, due to the condition of `config.linear_forgetting < length(mus)` in
# Samplers.adaptive_parzen_normal, N can never be 0 even if config.linear_forgetting is 0.

# It also cannot be the case in HP.PChoice because although forgetting_weights is called in
# Samplers.categoricalindex, the user cannot submit 0 choices due to the PChoice's inner
# constructor.

@testset "Test a variety of N and lf" begin

    # N is a number of observations/mus, lf is based on config and indicates a positive linear
    # forgetting value beyond which older points are de-weighted, e.g. if lf == 25, only
    # the last 25 values will be weighted full. The rest will be linearly de-weighted.
    history_len = 50
    for N in 0:history_len, lf in 0:history_len
        case = TreeParzen.ForgettingWeights.forgetting_weights(N, lf)
        @test length(case) == N
        if !isempty(case)
            # There is no point putting in a history item that is to be completely forgotten
            @test minimum(case) > 0
            @test maximum(case) <= 1
        end

        # Test the fully-remembered part
        ones = [x for x in case if x == 1]
        # If don't you want to forget anything
        if N - lf <= 0
            @test length(ones) == N
        else
            @test length(ones) == lf
        end

        if !isempty(ones)
            @test unique(ones) == [1]
        end

        # Test the ramp
        ramp = [x for x in case if x != 1]
        # If don't you want to forget anything
        if N - lf <= 0
            @test isempty(ramp)
        else
            @test length(ramp) == N - lf
        end

        if !isempty(ramp)
            full_ramp = vcat(0, ramp, 1)
            # Check that the ramp is linear. It might be possible that one or
            # two additional values appear which are very similar due to
            # floating point errors.
            ramp_diff = unique(round.(diff(full_ramp), digits = 3))
            @test length(ramp_diff) == 1
        end
    end
end

end # module TestForgettingWeights
true


text:
export Configuration, Orbital, Kink, T2, T4, ImgTime, excite!, excite, Kinks, haskey, Kinks, times, excitations, keys, values, getindex, (==)

"""
Abstract type for single-particle basis states, implementation is required for each model.
"""
abstract type Orbital end

"""
Type alias for imaginary time.
`FixedPointNumbers` are used since these are stable for `==` and are thus stable as times in `Dict`.
"""
const ImgTime = Fixed{Int64,60}

@doc raw"""
    T2{T}

Parametric type representing a 1-particle excitation by specifying
a transition from one state (annihilator) to another (creator).
In the occupation number representation, this reads

$a^{\dagger}_i a_j$

with a creator orbital `i` and an annihilator orbital `j`.
The single-particle basis is represented by the type parameter `T`.
"""
struct T2{T}
  " creator "
  i :: T
  " annihilator "
  j :: T
end

@doc raw"""
    T4{T}

Parametric type representing a 2-particle excitation by specifying
a transition from two state (annihilators) to two other (creator) states.
In the occupation number representation, this reads

$a^{\dagger}_i a^{\dagger}_j a_k a_l$

with creator orbitals `i`, `j` and an annihilator orbitals `k`, `l`.
The single-particle basis is represented by the type parameter `T`.
"""
struct T4{T}
  " creator "
  i :: T
  j :: T
  " annihilator "
  k :: T
  l :: T
end

"""
    const Kink{T} = Union{T2{T}, T4{T}}

Parametric type representing either a one- or a two-particle scattering event
given by types 'T2{T}' and 'T4{T}', respectively.
This term is motivated by the geometrical interpretation of a path
in Configuration Path Integral Monte Carlo formulation of the partition function
in occupation number representation.
Discrete occupations are given by horizontal lines of (quasi-)particles
occupying the corresponding single-particle states and where transitions between
these states are thus indicated by vertical 'kinks' connecting the orbital lines
that are part of the transition.
"""
const Kink{T} = Union{T2{T}, T4{T}}

# outer constructor method to construct a T2 kink, inferring the type parameter from the arguments
Kink(i,j) = T2(i,j)
# outer constructor method to construct a T4 kink, inferring the type parameter from the arguments
Kink(i,j,k,l) = T4(i,j,k,l)

# outer constructor method to extract a kink from a pair where the second element is a kink.
Kink(p::Pair{S,T} where {T<:Kink} where {S}) = p[2]# first substitute S, then T

"""
    const Kinks{T} = Vector{Pair{ImgTime, Kink{T}}}

Structure for storing excitations and their imaginary times (kinks).

Can also be constructed by passing such pairs directly.

    Kinks(pairs::Pair{ImgTime,<:Kink{T}}...)
"""
const Kinks{T} = Vector{Pair{ImgTime, Kink{T}}}

Kinks(pairs::Pair{ImgTime,<:Kink{T}}...) where {T}  = reduce(push!, pairs, init=Kinks{T}())

"""
    excitations(ck::Kinks)

Return a list of the excitations of a `Kinks`-object, used to allow the use of dictionary syntax.
"""
excitations(ck::Kinks) = last.(ck)

"""
    times(ck::Kinks)

Return a list of the imaginary-times of a `Kinks`-object, used to allow the use of dictionary syntax.
"""
times(ck::Kinks) = first.(ck)

"""
    Base.haskey(ck::Kinks, key::ImgTime)

Check if a `Kinks`-object contains a kink at a specific time.
"""
Base.haskey(ck::Kinks, key::ImgTime) = in(key, times(ck))

"""
    Base.getindex(kinks::Kinks{T}, τ::ImgTime) where {T}

Allow the `kinks[τ]`-syntax to retrieve a `Kink` at a specific `ImgTime`.
"""
function Base.getindex(kinks::Kinks{T}, τ::ImgTime) where {T}
    τ_next, k = kinks[searchsortedfirst(kinks, by=first, τ)]
    if τ_next != τ
        throw(KeyError(τ))
    else
        return k
    end
end

function Base.setindex!(ck::Kinks, kink::Kink, τ::ImgTime)
    add!(ck::Kinks, τ::ImgTime, kink::Kink)
end

"""
Multi-particle trajectory in imaginary time.

`Configuration{T}` is a parametric type depending on the single-particle basis `T <: Orbital`.

**Fields**
- `occupations :: Set{T}`  -- orbitals which are occupied initially (`τ=0`)
- `kinks :: Kinks{T}`      -- collection of excitations and their imaginary times (kinks)

"""
mutable struct Configuration{T}
  occupations :: Set{T}
  kinks :: Kinks{T}
end

" outer constructor method for a configuration with occupations given by o and kinks given by k. "
Configuration(o::Set{T}, k) where {T} = Configuration(o, Kinks(k) )
" outer constructor method for a configuration with occupations given by o and no kinks. "
Configuration(o::Set{T}) where {T} = Configuration(o, Kinks{T}())
" outer constructor method for a configuriation with no occupations and kinks given by k. "
Configuration(k::Kinks{T}) where {T} = Configuration(Set{T}(), k)
" outer constructor method for a configuration with no occupations and kinks as given by varargs p... of Pair{ImgTime,<:Kink}, which are passed to the Kinks constructor "
Configuration(p::Pair{ImgTime,<:Kink{T}}...) where {T} = Configuration(Kinks(p...))
" outer constructor method for a configuration with occupations given by o and kinks as given by varargs p... of Pair{ImgTime,<:Kink}, which are passed to the Kinks constructor "
Configuration(o::Set{T}, p::Pair{ImgTime,<:Kink{T}}...) where {T} = Configuration(o, Kinks(p...))

" outer constructor method for empty Configurations{T} "
Configuration{T}() where T = Configuration(Set{T}())

Base.:(==)(c1::Configuration, c2::Configuration) = (c1.occupations == c2.occupations) && (c1.kinks == c2.kinks)

"""
    drop_orbs(occ, orbs)

Return a copy of the occupation `occ` without the orbitals in `orbs`.
"""
drop_orbs(occ, ::Nothing) = occ
drop_orbs(occ, orbs) = filter(o -> !in(o,orbs), occ)

"""
    add_orbs(occ, orbs)

Return a copy of the occupation `occ` with the orbitals in `orbs` added.
"""
add_orbs(occ, ::Nothing) = occ
add_orbs(occ, orbs) = union(occ, orbs)

"""
    drop_kinks(ck, s)

Return a copy of the kinks in `ck` without the kinks in `s`.
"""
drop_kinks(ck, ::Nothing) = ck
drop_kinks(ck, s) = filter(k -> !in(k,s), ck)

"""
    add_kinks(ck::Kinks, s)

Return a copy of the kinks in `ck` with the kinks in `s` added.
"""
add_kinks(ck::Kinks, ::Nothing) = ck

function add_kinks(ck::Kinks, s)
    ck_copy = deepcopy(ck)
    add_kinks!(ck_copy,s)
    return ck_copy
end

"""
    drop_orbs!(occ, s)

Remove the orbitals in `s` from the occupation `occ`.
"""
drop_orbs!(occ, ::Nothing) = nothing
drop_orbs!(occ, s) = filter!(o -> !in(o,s), occ)

"""
    add_orbs!(occ, s)

Add the orbitals in `s` to the occupation `occ`.
"""
add_orbs!(occ, ::Nothing) = nothing
add_orbs!(occ, s) = union!(occ, s)

"""
   drop_kinks!(ck, s)

Remove the kinks in `s` from `ck`.
"""
drop_kinks!(ck, ::Nothing) = nothing
drop_kinks!(ck, s) = filter!(k -> !in(k,s), ck)

"""
   add_kinks!(ck::Kinks, s)

Add the kinks in `s` to `ck`.
"""
add_kinks!(ck::Kinks, ::Nothing) = nothing
function add_kinks!(ck::Kinks, s)
    for k in s
        insert!(ck, searchsortedfirst(ck, by=first, first(k)), first(k) => last(k))
    end
end

"""
    basis(c::Configuration{T})

Return type parameter `T` of the configuration.
"""
basis(c::Configuration{T}) where T = T

"""
    orbs(::T2)

Return a Tuple of all orbitals that are affected by a `T2`-kink in the conventional ordering i, j.
"""
orbs(x::T2) = x.i, x.j

@doc raw"""
    orbs(::T4)

Return a Tuple of all orbitals that are affected by a `T4`-kink in the conventional ordering i, j, k, l.
This corresponds to the matrix element $w_{ijkl}$ in the same order.
"""
orbs(x::T4) = x.i, x.j, x.k, x.l

"""
    creators(::T4)

return a set of the two creators which are affected by a T4 kink
"""
creators(x::T4) = x.i, x.j

"""
    annihilators(x::T4)

return a set of the two annihilators which are affected by a T4 kink
"""
annihilators(x::T4) = x.k, x.l

"""
    ImgTime(t::Tuple{Nothing,Nothing})

Outer constructor method for a `Tuple{ImgTime,ImgTime}` from `Tuple{Nothing,Nothing}` which returns the bounds of the `ImgTime` interval as `Tuple (ImgTime(0), ImgTime(1))`
"""
ImgTime(t::Tuple{Nothing,Nothing}) = (ImgTime(0), ImgTime(1))

"""
    ImgTime(::Pair{ImgTime,<:Kink})

outer constructor method for an ImgTime from Pair{ImgTime,<:Kink}
this returns the first argument of the Pair
useful for automatic conversion
"""
ImgTime(p::Pair{ImgTime,<:Kink}) = first(p)

"""
    ImgTime(t::Tuple{Pair{ImgTime,<:Kink},Pair{ImgTime,<:Kink}})

Outer constructor method for a `Tuple{ImgTime,ImgTime}` from `Tuple{Pair{ImgTime,<:Kink},Pair{ImgTime,<:Kink}}` which returns the `ImgTime` of the Pairs as `Tuple{ImgTime,ImgTime}`
"""
ImgTime(t::Tuple{Pair{ImgTime,<:Kink},Pair{ImgTime,<:Kink}}) = (first(first(t)), first(last(t)))

@doc raw"""
    Δ(::ImgTime,::ImgTime)

return the length of the periodic interval between the two ::ImgTimes

the periodic distance of two imaginary times:

    `Δ(τ1,τ2) = τ2 - τ1` if `τ2 >= τ1` and
    `Δ(τ1,τ2) = 1 - (τ1 - τ2) = 1 + τ2 - τ1` else
"""
Δ(τ1::ImgTime,τ2::ImgTime) = τ1 > τ2 ? ImgTime(1) + τ2 - τ1 : τ2 - τ1

@doc raw"""
    excite(::Set{T}, ::T4{T})

Apply a kink to a set of basis states, i.e. return a set of basis states where
the states specified by the creators of the kink are added and the states specified
by the annihilators of the kink are dropped from the given set of basis states.
This has the physical meaning of a two-particle scattering event where two (quasi-)particles
in a many-body state change the single-particle states they occupy.
In occupation number representation this reads

$a^{\dagger}_i a^{\dagger}_j a_k a_l |\{n\}\rangle$

for creator orbitals `i` and `j` and annihilator orbitals `k` and `l`.
This function assumes fermionic particle statistics (Pauli principle),

$a^{\dagger}_i a^{\dagger}_i = a_i a_i = 0$

i.e. the target (creator) states must not be occupied and
the initial (annihilator) states must be occupied in the given set of states.
"""
function excite(o::Set{T}, κ::T4{T}) where T
  @assert ( in(κ.k, o) & in(κ.l, o) ) "Kink ($(κ.i),$(κ.j),$(κ.k),$(κ.l)) cannot be applied: one or two of the annihilators $(κ.k), $(κ.l) is not occupied."
  @assert ( !in(κ.i, o) & !in(κ.j, o) ) "Kink ($(κ.i),$(κ.j),$(κ.k),$(κ.l)) cannot be applied: one or two of the creators $(κ.i), $(κ.j) is already occupied. (Pauli principle)"
  union(setdiff(o, Set([κ.k,κ.l])), Set([κ.i, κ.j]))
end

"""
    excite(o::Set{T}, κ::Pair{ImgTime,T4{T}})

Apply a `T4`-kink to a set of basis states for a pair of a time and a `T4`-kink. This is useful for iterating over a `Kinks`-object when calculating the occupations at a specific time.
"""
excite(o::Set{T}, κ::Pair{Fixed{Int64,60},Union{T2{T}, T4{T}}}) where T = excite(o, last(κ))
excite(o::Set{T}, κ::Pair{ImgTime,T4{T}}) where T = excite(o, last(κ))

"""
    excite!(::Set{T}, i::T, j::T, k::T, l::T) where {T}
    excite!(o::Set{T}, κ::T4{T})

Apply an excitation in-place to a set of basis states that is given by creating orbitals i, j
and annihilating the orbitals k, l.
"""
function excite!(o::Set{T}, i::T, j::T, k::T, l::T) where {T}
    @assert (in(k, o) & in(l, o)) "Kink ($(i),$(j),$(k),$(l)) cannot be applied: one or two of the annihilators $(k), $(l) is not occupied. (Pauli-Principle)"
    @assert (!in(i, o) & !in(j, o)) "Kink ($(i),$(j),$(k),$(l)) cannot be applied: one or two of the creators $(i), $(j) is already occupied. (Pauli-Principle)"
    delete!(o, k)
    delete!(o, l)
    push!(o, i)
    push!(o, j)
end

excite!(o::Set{T}, κ::T4{T}) where T = excite!(o, κ.i, κ.j, κ.k, κ.l)

"""
    occupations_at(o::Set{T}, kinks::Kinks{T})

Return the occupied orbitals after applying all kinks to initial occupation.
"""
function occupations_at(o::Set{T}, kinks::Kinks{T}) :: Set{T} where {T}
  foldl(excite, kinks; init=o)
end

"""
    occupations_at(c::Configuration, τ::ImgTime)

Return the occupied orbitals to the right of τ.
"""
function occupations_at(c::Configuration, τ::ImgTime)
  occupations_at(c.occupations, filter(x -> first(x) <= τ, c.kinks))
end

"""
    next(a, τ)

Return index of the first kink after to the given `τ`. If there is no kink with `::ImgTime` larger than `τ`, return index of first kink.
Throw a `DomainError` if an empty list is passed as first argument.

This function assumes that the kinks in `a` are ordered with respect to their imaginary time!
"""
function next(a, τ)
    if isempty(a)
        throw(DomainError(" next(x::Kinks, τ) is not supported for empty x. It makes no sense to get the element next to some element for an empty collection. "))
    end
    index = searchsortedfirst(a, τ, by=first) # find first element x ∈ a with first(x) >= τ
    if index == length(a) + 1 # if there is no element x ∈ a with first(x) >= τ
        return 1#  return index of the first element in a
    end
    # catch equality:
    # if there there is a kink at time τ, `searchsortedfirst` returns the index of that kink
    if τ == first(a[index])
        index += 1 # increase index by one
    end
    # catch if kink at τ is the last kink in a:
    # then the index after the index of the kink at τ is length(a) + 1
    if index == length(a) + 1
        return 1 # return index of the first element
    else
        return index # return the index
    end
end

"""
     next_affecting(a::Any, orbs::Any, τ::Any)

Return the index of the closest kink to the right of τ
that affects one of the orbitals in orbs.
If no orbital in orbs is affected by any kink return 0.
"""
function next_affecting(a, orbs, τ)
    if isempty(a)
        return 0
    end
    i = next(a, τ)
    if any( orbs .∈ ( last(a[i]) ,) ) # if any orbs affected by the kink at i
        return i
    end
    # otherwise, iterate forward periodically until first affecting kink is reached
    for _ = 1:length(a) - 1 # iterate until before the same kink it reached. note: this should not occur if o is affected by any kink since there should always be at least two kinks for each affected orbital due to periodicity (i.e. all creators have to pair up with a corresponding annihilator)
        i += 1 # increment
        if i == length(a) + 1 # periodicity: if last index is reached, start from beginning of collection
            i = 1
        end
        if any( orbs .∈ ( last(a[i]) ,) )# if any orbs are affected by the kink at i
            return i # if o is affected by the kink at i
        end
    end
    return 0
end

"""
    prev(a, τ)

Return the index of the first kink earlier than `τ::ImgTime`. If there is no such kink the last kink is returned.
Throw a `DomainError` if an empty list is passed as first argument.

This function assumes that the kinks in `x` are ordered with respect to their imaginary time!
"""
function prev(a, τ)
    if isempty(a)
        throw(DomainError(" prev(x::Kinks, τ) is not supported for empty x. It makes no sense to get the element previous to some element for an empty collection. "))
    end
    index = searchsortedlast(a, τ, by=first) # find last element x ∈ a with first(x) <= τ
    if index == 0 # if there is no element x ∈ a with first(x) <= τ
        return length(a)#  return index of the last element in a
    end
    # catch equality:
    # if there there is a kink at time τ, `searchsortedlast` returns the index of that kink
    if τ == first(a[index])
        index -= 1 # reduce index by one
    end
    # catch if kink at τ is the first kink in a:
    # then the index before the index of the kink at τ is 0
    if index == 0
        return length(a) # return index of the last element
    else
        return index # return the index
    end
end

"""
     prev_affecting(a::Any, orbs::Any, τ::Any)

Return the index of the closest kink to the left of τ
that affects one of the orbitals in orbs.
If no orbital in orbs is affected by any kink return 0.
"""
function prev_affecting(a, orbs, τ)
    if isempty(a)
        return 0
    end
    i = prev(a, τ)
    if any( orbs .∈ ( last(a[i]), ) ) # if any orbs are affected by the kink at i
        return i
    end
    # otherwise, iterate backwards periodically until first affecting kink is reached
    for _ = 1:length(a) - 1 # iterate until before the same kink it reached. note: this should not occur if o is affected by any kink since there should always be at least two kinks for each affected orbital due to periodicity (i.e. all creators have to pair up with a corresponding annihilator)
        i -= 1 # decrement
        if i == 0 # periodicity: if first index is reached, start from end of collection
            i = length(a)
        end
        if any( orbs .∈ ( last(a[i]), ) ) # if any orbs are affected by the kink at i
            return i
        end
    end
    return 0
end

"""
     kinks_affecting_orbs(itr::Any, itr::Any)

Return all kinks in the first argument that affect one ore more of the orbitals the second argument.
The first itr is expected to contain tuples or at least types on which `last()` is defined
and `last(i)` is intended to contain a ::T4 for all elements `i ∈ itr` or at least a type
which has fields with the names `i`, `j`, `k`, `l`.
"""
kinks_affecting_orbs(ck, os) = filter( p -> any( orbs(last(p)) .∈ (os,) ), ck)
kinks_affecting_orbs(c::Configuration, os) = kinks_affecting_orbs(c.kinks, os)

"""
    adjacent_kinks(itr::Any, ::Any)

Return a tuple of the kink in `itr`
that is closest right of the time in the second argument.
This expects that methods for the functions `next()` and `prev()` are defined
for the given argument types.
Used for getting a tuple of the neighbouring kinks from some imaginary time.
"""
adjacent_kinks(ck::Any, τ::Any) = ck[prev(ck, τ)], ck[next(ck, τ)]

"""
    adjacent_kinks_affecting_orbs(::Any, ::Any, ::Any)

Return a tuple of the closest kink to the right and the closest kink to the left of some imaginary time
that affect one of the orbitals in some collection. If there are no such kinks, return (nothing,nothing).
The first argument is expected to be an iterable which contains pairs as elements, with
the first element of each pair containing an imaginary time and the second element of each pair containing a kink.
The second argument is expected to be a collection of orbitals, that are to be considered regarding affection with one of the
kinks in the first argument.
The third argument is expected to be an imaginary time from which the neighbouring kinks affected by any orbitals from the
second argument are to be determined."""
function adjacent_kinks_affecting_orbs(ck, os, τ)
    l = prev_affecting(ck, os, τ)

    if l == 0
        return (nothing,nothing)
    else
        return ck[l], ck[next_affecting(ck, os, τ)]
    end
end
adjacent_kinks_affecting_orbs(c::Configuration, os, τ) = adjacent_kinks_affecting_orbs(c.kinks, os, τ)

"""
     τ_prev_affecting(::Any, ::Any, ::Any)

Return the ImgTime of the closest kink to the left of τ
that affects one of the orbitals in os.
"""
function τ_prev_affecting(ck, os, τ)  # TODO: remove
    first(ck[prev_affecting(ck, os, τ)])
end

"""
     τ_next_affecting(::Any, ::Any, ::Any)

Return the ImgTime of the closest kink to the right of τ
that affects one of the orbitals in os.
"""
function τ_next_affecting(ck, os, τ)  # TODO: remove
    first(ck[next_affecting(ck, os, τ)])
end

"""
     τ_borders(::Kinks{T}, orbs, ::ImgTime) where {T <: Orbital}

Return a tuple of
the ImgTime of the closest kink to the right and
the ImgTime of the closest kink to the left of τ
that affect one of the orbitals in orbs.
If no orbital in os is affected by and kink from the collection in the first argument,
return a tuple of the interval bounds (ImgTime(0), ImgTime(1))."""
τ_borders(ck::Kinks{T}, orbs, τ::ImgTime) where {T <: Orbital} = ImgTime(adjacent_kinks_affecting_orbs(ck, orbs, τ))  # TODO: remove
τ_borders(c::Configuration{T}, orbs, τ::ImgTime) where {T <: Orbital} = τ_borders(c.kinks, orbs, τ)

"""
    isunaffected(Kinks, orb)

Return if an orbital is not affected by any kink.
"""
isunaffected(kinks, orb) = all(kink -> orb ∉ last(kink), kinks)

Base.in(i, k::T2) = i == k.i || i == k.j
Base.in(i, k::T4) = i == k.i || i == k.j || i == k.k || i == k.l

"""
    in_open_interval(τ, τ_first, τ_last)

Check if `τ` lies in the open `ImgTime`-interval `(τ_first,`τ_last)`, assuming that `τ_first` is the left border and `τ_last` the right one.
"""
in_open_interval(τ, τ_first, τ_last) = (τ_first != τ != τ_last != τ_first) && ((τ_first < τ_last) == ((τ < τ_first) != (τ < τ_last)))

"""
    isunaffected_in_interval(kinks, orb, τ_first::ImgTime, τ_last::ImgTime)

Return if an orbital is not affected by any of the kinks from `ck` in the open interval `(τ_first,τ_last)`.
"""
isunaffected_in_interval(kinks, orb, τ_first::ImgTime, τ_last::ImgTime) = all(kink -> ((orb ∉ last(kink)) || !in_open_interval(first(kink), τ_first, τ_last) ), kinks)

"""
    time_ordered_orbs(ck::Kinks{T}) where {T <: Orbital}

Get a list of orbitals that affect each kink in the time-ordering of the kinks and in the conventional ordering i, j, k, l.
"""
function time_ordered_orbs(ck::Kinks{T}) where {T <: Orbital}
    if isempty(ck)
        return Array{T,1}()
    else
        return collect( Iterators.flatten( orbs(k) for k in excitations(ck) ) )
    end
end

@doc raw"""
    ladder_operator_order_factor(sortedlist::Array{<:Orbital,1})

Returns $1$ or $-1$ depending on the order of all ladder operators as given by a list of orbitals.
"""
function ladder_operator_order_factor(sortedlist::Array{<:Orbital,1})
    phase_factor = 1
    while !isempty(sortedlist)
        index = 1
        op_orb = sortedlist[index]
        while sortedlist[index + 1] != op_orb
            sortedlist[index:index+1] = [sortedlist[index+1], sortedlist[index]]
            phase_factor *= -1
            index += 1
            @assert(index <= length(sortedlist))
        end
        deleteat!(sortedlist, [index, index+1])
    end
    return phase_factor
end

@doc raw"""
    ladder_operator_order_factor(ck::Kinks)

Returns $1$ or $-1$ depending on the order of all ladder operators as given by the orbitals that affect each kink in the time-ordering of the kinks and in the conventional ordering i, j, k, l,
used in the sign estimator.
"""
ladder_operator_order_factor(ck::Kinks) = ladder_operator_order_factor(time_ordered_orbs(ck))

"""
    is_type_1(left_kink::T4, right_kink::T4)

Returns True if left_kink and right_kink are entangled in a Type-1 way.
This does not check wether the two kinks are neighbouring."""
function is_type_1(left_kink::T4, right_kink::T4)
  if ((length(intersect(Set([left_kink.i, left_kink.j]),Set([right_kink.k, right_kink.l]))) == 1) &
        (length(intersect(Set([left_kink.k, left_kink.l]), Set([right_kink.i, right_kink.j]))) == 0)) ||
    ((length(intersect(Set([left_kink.i, left_kink.j]),Set([right_kink.k, right_kink.l]))) == 0) &
          (length(intersect(Set([left_kink.k, left_kink.l]), Set([right_kink.i, right_kink.j]))) == 1))
    return true
  else
    return false
  end
end

"""
    right_type_1_chain_length(ck, τ, count = 0)

Returns the length of the chain of type-1-entaglements starting with the Kink at τ counting to the right.
"""
function right_type_1_chain_length(ck::Kinks, τ, counted_τs = [])
    next_kink = ck[next_affecting(ck, orbs(ck[τ]), τ)]
    if is_type_1(ck[τ], last(next_kink)) & !in(τ, counted_τs)
        push!(counted_τs,τ)
        return right_type_1_chain_length(ck, first(next_kink), counted_τs)
    else
        return length(counted_τs)
    end
end

"""
    left_type_1_chain_length(ck::Kinks, τ, counted_τs = [])

Returns the length of the chain of type-1-entaglements starting with the Kink at τ counting to the left.
"""
function left_type_1_chain_length(ck::Kinks, τ, counted_τs = [])
    prev_kink = ck[prev(kinks_affecting_orbs(ck, orbs(ck[τ])), τ)]
    if is_type_1(last(prev_kink), ck[τ]) & !in(τ, counted_τs)
        push!(counted_τs,τ)
        return left_type_1_chain_length(ck, first(prev_kink), counted_τs)
    else
        return length(counted_τs)
    end
end

"""
    longest_type_1_chain_length(ck::Kinks)

Returns the longest chain of type-1-entaglements in ck.
"""
function longest_type_1_chain_length(ck::Kinks)
    longest_length = 0
    for (τ, kink) in ck
        longest_length = max(right_type_1_chain_length(ck, τ),
                                left_type_1_chain_length(ck, τ), longest_length)
    end
    return longest_length
end

"""
    right_type_1_count(ck::Kinks)

Returns the number of kinks that are type 1 entagled with their right neighbour.
"""
function right_type_1_count(ck::Kinks)
    count = 0
    for (τ, kink) in ck
        if is_type_1(kink, last(ck[next_affecting(ck, orbs(kink),τ)]))
            count += 1
        end
    end
    return count
end

"""
    random_shuffle(::Any, ::Any)

shuffle or not shuffle the given tuple with equal propability:
Return either the same tuple, or the tuple in reverse order,
with equal probability 0.5.
"""
function random_shuffle(i, j)
    if rand() < 0.5
        return i, j
    else
        return j, i
    end
end

"""
    shuffle_annihilators(::Any, ::Any, ::Any, ::Any)

Return a tuple where the last two arguments are randomly shuffled with equal probability 0.5
in comparison to the ordering of the input.
"""
shuffle_annihilators(i, j, k, l) = i, j, random_shuffle(k, l)...
shuffle_annihilators(κ::T4) = T4(shuffle_annihilators(κ.i, κ.j, κ.k, κ.l)...)

"""
    shuffle_creators(::Any, ::Any, ::Any, ::Any)

Return a tuple where the first two arguments are randomly shuffled with equal probability 0.5
in comparison to the ordering of the input.
"""
shuffle_creators(i, j, k, l) = random_shuffle(i, j)..., k, l
shuffle_creators(κ::T4) = T4(shuffle_creators(κ.i, κ.j, κ.k, κ.l)...)

"""
    shuffle_indices(::Any, ::Any, ::Any, ::Any)

Return a tuple where both the first two arguments and the last two arguments
are randomly shuffled with equal probability 0.5 respectively
in comparison to the ordering of the input.
"""
shuffle_indices(i, j, k, l) = random_shuffle(i, j)..., random_shuffle(k, l)...
shuffle_indices(κ::T4) = T4(shuffle_indices(κ.i, κ.j, κ.k, κ.l)...)

"""
    kinks_from_periodic_interval(::Kinks, τ1, τ2)

return kinks with τ ∈ (τ1,τ2) if τ1 < τ2 and τ ∈ (τ2,1) ∪ (0,τ1) if τ1 > τ2
"""
function kinks_from_periodic_interval(ck::Kinks, τ1, τ2)
    filter(x-> in_open_interval(first(x), τ1, τ2), ck)
end

"""
    times_from_periodic_interval(::Kinks, ::ImgTime, ::ImgTime)

return a list of all times of kinks with τ ∈ (τ1,τ2) if τ1 < τ2 or τ ∈ (τ2,1) ∪ (0,τ1) if τ1 > τ2
in the periodic ordering suggested by the relation of the first time-argument τ1 to the second time-argument τ2
"""
function times_from_periodic_interval(ck::Kinks, τ1::ImgTime, τ2::ImgTime)
    times(filter(x-> in_open_interval(first(x), τ1, τ2), ck))
end


text:
macro communative(ex)
    switched = copy(ex)
    @assert ex.head ∈ (:function, :(=))
    switched.args[1].args[2] = ex.args[1].args[3]
    switched.args[1].args[3] = ex.args[1].args[2]
    quote
        $(esc(ex))
        $(esc(switched))
    end
end


text:
include("PSCOPF.jl")

using Dates

using .PSCOPF

# load network
data_path = joinpath(@__DIR__, "..", "2buses")
network = PSCOPF.Data.pscopfdata2network(data_path)

created_instance_path = joinpath(data_path, "copied")

PSCOPF.PSCOPFio.write(created_instance_path, network)

# Generate uncertainties
uncertainties_distribution = PSCOPF.PSCOPFio.read_uncertainties_distributions(network, data_path)
nb_scenarios = 5

ts1 = Dates.DateTime("2015-01-01T11:00:00")
TS = PSCOPF.create_target_timepoints(ts1)
ECHs = []
for mode in [PSCOPF.PSCOPF_MODE_1, PSCOPF.PSCOPF_MODE_2, PSCOPF.PSCOPF_MODE_3]
    push!(ECHs, PSCOPF.generate_ech(network, TS, mode) )
end
horizon_timepoints = sort(unique(Iterators.flatten(ECHs)))

uncertainties = PSCOPF.generate_uncertainties(network, TS, horizon_timepoints,
                                            uncertainties_distribution, nb_scenarios)

PSCOPF.PSCOPFio.write(created_instance_path, uncertainties)

# Launch mode 1
mode = PSCOPF.PSCOPF_MODE_1

network = PSCOPF.Data.pscopfdata2network(created_instance_path)
un = PSCOPF.PSCOPFio.read_uncertainties(created_instance_path)
ts1 = Dates.DateTime("2015-01-01T11:00:00")
TS = PSCOPF.create_target_timepoints(ts1)
ECH = PSCOPF.generate_ech(network, TS, mode)

sequence = PSCOPF.generate_sequence(network, TS, ECH, mode)
exec_context = PSCOPF.PSCOPFContext(network, TS, mode)
PSCOPF.run!(exec_context, sequence)

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
####
#### Miscellaneous utils
####

ndigits(x::Int) = x == 0 ? 1 : floor(Int, log10(abs(x))) + 1
logrange(a::Real, b::Real, len::Int) = (r = exp.(range(log(a), log(b); length = len)); r[1] = a; r[end] = b; return r)
normcdf(x::T) where {T} = erfc(-x/sqrt(T(2)))/2 # Cumulative distribution for normal distribution
normccdf(x::T) where {T} = erfc(x/sqrt(T(2)))/2 # Compliment of normcdf, i.e. 1 - normcdf(x)

@inline mul_im(z::Complex) = Complex(-imag(z), real(z)) # optimized i*(a+b*i) = -b+a*i

@inline basisvector(::Type{SVector{D,T}}, i::Int) where {D,T} = SVector{D,T}(ntuple(d -> T(d == i), D))

function meshgrid(::Type{T}, iters...) where {T}
    A = [T(xs) for xs in Iterators.product(iters...)]
    reshape(A, length.(iters)...)
end
meshgrid(iters...) = meshgrid(Tuple, iters...)

@inline function SplitCartesianIndices(sz::NTuple{N,Int}, ::Val{M}) where {N,M}
    @assert 0 <= M <= N
    sz1 = sz[1:M]
    sz2 = sz[M+1:N]
    return CartesianIndices(sz1), CartesianIndices(sz2)
end
@inline SplitCartesianIndices(x::AbstractArray{<:Any,N}, ::Val{M}) where {N,M} = SplitCartesianIndices(size(x), Val(M))

@generated function fieldsof(::Type{T}, ::Type{C} = Tuple) where {T,C}
    fields = fieldnames(T) # fieldnames(T) allocates; hoist to generated function
    return C <: Tuple ?
        :($fields) : # default to returning tuple of field symbols
        :($(C(Symbol[fields...]))) # call container constructor on vector of symbols
end

@inline floattype(xs::Tuple) = float(promote_type(map(typeof, xs)...))
@inline floattype(xs::NamedTuple) = floattype(Tuple(xs))

function set_diag!(A::AbstractMatrix, val)
    @inbounds @simd ivdep for i in 1:min(size(A)...)
        A[i,i] = val
    end
    return A
end

function set_top!(A::AbstractArray, B::AbstractArray)
    @inbounds @simd ivdep for I in CartesianIndices(B)
        A[I] = B[I]
    end
    return A
end

function find_nearest(r::AbstractRange, x::Number)
    idx = x <= r[1] ? 1 :
          x >= r[end] ? length(r) :
          clamp(round(Int, 1 + (x - r[1]) / step(r)), 1, length(r))
    r[idx], idx # nearest value in r to x and corresponding index
end

function local_gridsearch(f, xs, i0)
    get(j) = @inbounds xs[clamp(j, firstindex(xs), lastindex(xs))]
    i = i0
    x⁻, x, x⁺ = get(i-1), get(i), get(i+1)
    y⁻, y, y⁺ = f(x⁻), f(x), f(x⁺)
    while !(y⁻ ≥ y ≤ y⁺) # search for local min
        if y⁻ < y
            i -= 1 # shift left
            x⁻, x, x⁺ = get(i-1), x⁻, x
            y⁻, y, y⁺ = f(x⁻), y⁻, y
        elseif y⁺ < y
            i += 1 # shift right
            x⁻, x, x⁺ = x, x⁺, get(i+1)
            y⁻, y, y⁺ = y, y⁺, f(x⁺)
        else
            break
        end
    end
    return (; x, y, i)
end

function mapfind(f, finder, xs)
    ys = map(f, xs)
    y, i = finder(ys)
    xs[i], y, i
end
mapfindmax(f, xs) = mapfind(f, findmax, xs)
mapfindmin(f, xs) = mapfind(f, findmin, xs)

# Threaded `foreach` construct, borrowing implementation from ThreadTools.jl:
#
#   https://github.com/baggepinnen/ThreadTools.jl/blob/55aaf2bbe735e52cefaad143e7614d4f00e312b0/src/ThreadTools.jl#L57
#
# Updated according to suggestions from the folks at DataFrames.jl:
#
#   https://github.com/jondeuce/DECAES.jl/issues/37
function tforeach(f, x::AbstractArray; blocksize::Integer = default_blocksize())
    nt = Threads.nthreads()
    len = length(x)
    if nt > 1 && len > blocksize
        @sync for p in split_indices(len, blocksize)
            Threads.@spawn begin
                @simd ivdep for i in p
                    f(@inbounds x[i])
                end
            end
        end
    else
        @simd ivdep for i in eachindex(x)
            f(@inbounds x[i])
        end
    end
    return nothing
end

function split_indices(len::Integer, basesize::Integer)
    len′ = Int64(len) # Avoid overflow on 32-bit machines
    np = max(1, div(len′, basesize))
    return (Int(1 + ((i - 1) * len′) ÷ np) : Int((i * len′) ÷ np) for i in 1:np)
end

default_blocksize() = 64

####
#### Timing utilities
####

tic() = time()
toc(t) = tic() - t

function hour_min_sec(t)
    hour = floor(Int, t / 3600)
    min = floor(Int, (t - 3600 * hour) / 60)
    sec = floor(Int, t - 3600 * hour - 60 * min)
    return (; hour, min, sec)
end

function pretty_time(t)
    if isnan(t) || isinf(t)
        "--h:--m:--s"
    else
        hh, mm, ss = hour_min_sec(t)
        lpad(hh, 2, "0") * "h:" * lpad(mm, 2, "0") * "m:" * lpad(ss, 2, "0") * "s"
    end
end

@with_kw mutable struct DECAESProgress
    progmeter::Progress
    io::IO
    iobuf::IOBuffer
    iolock::ReentrantLock = Threads.ReentrantLock()
    last_msg::AbstractString = ""
end

function DECAESProgress(io::IO, n::Int, desc::AbstractString; kwargs...)
    iobuf = IOBuffer()
    DECAESProgress(
        io = io,
        iobuf = iobuf,
        progmeter = Progress(n;
            dt = 0.0, desc = desc, color = :cyan, output = iobuf, barglyphs = BarGlyphs("[=> ]"),
            kwargs...
        )
    )
end
DECAESProgress(n::Int, desc::AbstractString; kwargs...) = DECAESProgress(stderr, n, desc; kwargs...)

function ProgressMeter.next!(p::DECAESProgress)
    next!(p.progmeter)
    msg = String(take!(p.iobuf)) # take!(::IOBuffer) is threadsafe
    if !isempty(msg)
        msg = replace(msg, "\r" => "")
        msg = replace(msg, "\u1b[K" => "")
        msg = replace(msg, "\u1b[A" => "")
        if msg != p.last_msg
            lock(p.iolock) do
                println(p.io, msg)
                p.last_msg = msg
            end
        end
    end
end

printheader(io, s) = (println(io, ""); printstyled(io, "* " * s * "\n"; color = :cyan))
printbody(io, s) = println(io, s)

# Macro for timing arbitrary code snippet and printing time
macro showtime(io, msg, ex)
    quote
        local io = $(esc(io))
        printheader(io, $(esc(msg)) * " ...")
        local val
        local t = @elapsed val = $(esc(ex))
        printheader(io, "Done ($(round(t; digits = 2)) seconds)")
        val
    end
end

####
#### Logging
####

# https://discourse.julialang.org/t/write-to-file-and-stdout/35042/3
struct Tee{T <: Tuple} <: IO
    streams::T
end
Tee(streams::IO...) = Tee(streams)
Base.flush(t::Tee) = tee(t)
Base.write(t::Tee, x::Array) = tee(Base.write, t, x) # needed to resolve method ambiguity
Base.write(t::Tee, args...) = tee(Base.write, t, args...)
Base.print(t::Tee, args...) = tee(Base.print, t, args...)
Base.println(t::Tee, args...) = tee(Base.println, t, args...)
Base.printstyled(t::Tee, args...; kwargs...) = tee(Base.printstyled, t, args...; kwargs...)

function tee(f, t::Tee, args...; kwargs...)
    for io in t.streams
        f(io, args...; kwargs...)
        flush(io)
    end
end
tee(t::Tee, args...; kwargs...) = tee(io -> nothing, t, args...; kwargs...)

function tee_capture(f; logfile = tempname(), suppress_terminal = false, suppress_logfile = false)
    open(suppress_logfile ? tempname() : logfile, "w+") do io
        io = Tee(suppress_terminal ? devnull : stderr, io)
        logger = ConsoleLogger(io)
        with_logger(logger) do
            f(io)
        end
    end
end

# https://discourse.julialang.org/t/redirect-stdout-and-stderr/13424/3
function redirect_to_files(f, outfile, errfile)
    open(outfile, "w") do out
        open(errfile, "w") do err
            redirect_stdout(out) do
                redirect_stderr(err) do
                    f()
                end
            end
        end
    end
end
redirect_to_tempfiles(f) = redirect_to_files(f, tempname() * ".log", tempname() * ".err")

function redirect_to_devnull(f)
    with_logger(ConsoleLogger(devnull)) do
        redirect_to_tempfiles() do
            f()
        end
    end
end

####
#### Optimizers
####

struct ADAM{N,T}
    η::T
    β::SVector{2,T}
    mt::SVector{N,T}
    vt::SVector{N,T}
    βp::SVector{2,T}
end
function ADAM{N,T}(η = 0.001, β = (0.9, 0.999)) where {N,T}
    @assert N >= 1
    S2 = SVector{2,T}
    SN = SVector{N,T}
    ADAM{N,T}(T(η), S2(β), zero(SN), zero(SN), ones(S2))
end

function update(∇::SVector{N,T}, o::ADAM{N,T}) where {N,T}
    @unpack η, β, mt, vt, βp = o

    ϵ  = T(1e-8)
    βp = @. βp * β
    ηt = η * √(1 - βp[2]) / (1 - βp[1])
    mt = @. β[1] * mt + (1 - β[1]) * ∇
    vt = @. β[2] * vt + (1 - β[2]) * ∇^2
    Δ  = @. ηt * mt / (√vt + ϵ)

    return Δ, ADAM{N,T}(η, β, mt, vt, βp)
end

@inline xform_periodic(t::S, lb::S, ub::S) where {N, T, S <: SVector{N,T}} = S(ntuple(i -> clamp(((lb[i] + ub[i])/2) + ((ub[i] - lb[i])/2) * sinpi(t[i]), lb[i], ub[i]), N))
@inline ∇xform_periodic(t::S, lb::S, ub::S) where {N, T, S <: SVector{N,T}} = S(ntuple(i -> ((ub[i] - lb[i])/2) * T(π) * cospi(t[i]), N))
@inline inv_xform_periodic(x::S, lb::S, ub::S) where {N, T, S <: SVector{N,T}} = S(ntuple(i -> asin(clamp((x[i] - ((lb[i] + ub[i])/2)) / ((ub[i] - lb[i])/2), -one(T), one(T))) / T(π), N))

function optimize(∇f, x0::SVector{N,T}, lb::SVector{N,T}, ub::SVector{N,T}, o::ADAM{N,T}; maxiter::Int = 1, xtol_rel = T(1e-3)) where {N,T}
    x = x0
    t = inv_xform_periodic(x, lb, ub)
    for i in 1:maxiter
        # Change of variables x->t
        x = xform_periodic(t, lb, ub)
        dxdt = ∇xform_periodic(t, lb, ub)
        dfdx = ∇f(x)
        dfdt = dfdx .* dxdt

        # Update in t-space
        Δt, o = update(dfdt, o)
        t -= Δt

        # Check for convegence in x-space
        xold, x = x, xform_periodic(t, lb, ub)
        maximum(abs.(x - xold)) < max(maximum(abs.(x)), maximum(abs.(xold))) * xtol_rel && break
    end
    return x, o
end

####
#### Generate (moderately) realistic mock images
####

# Mock CPMG image
function mock_t2map_opts(::Type{T} = Float64; kwargs...) where {T}
    T2mapOptions{T}(;
        MatrixSize = (2,2,2),
        TE = 10e-3,
        nTE = 32,
        T2Range = (10e-3, 2.0),
        nT2 = 40,
        Reg = "lcurve",
        SetFlipAngle = 165.0,
        SetRefConAngle = 150.0,
        kwargs...
    )
end

# Mock CPMG image
function mock_image(o::T2mapOptions{T} = mock_t2map_opts(Float64); kwargs...) where {T}
    oldseed = Random.seed!(0)

    @unpack MatrixSize, TE, nTE = T2mapOptions(o; kwargs...)
    SNR = T(50)
    eps = T(10^(-SNR/20))

    α = o.SetFlipAngle === nothing ? T(165.0) : o.SetFlipAngle
    β = o.SetRefConAngle === nothing ? T(150.0) : o.SetRefConAngle
    mag() = T(0.85) .* EPGdecaycurve(nTE, α, TE, T(65e-3), T(1), β) .+
            T(0.15) .* EPGdecaycurve(nTE, α, TE, T(15e-3), T(1), β) # bi-exponential signal with EPG correction
    noise(m) = abs(m[1]) .* eps .* randn(T, size(m)) # gaussian noise of size SNR relative to signal amplitude
    noiseysignal() = (m = mag(); sqrt.((m .+ noise(m)).^2 .+ noise(m).^2)) # bi-exponential signal with rician noise

    M = zeros(T, (MatrixSize..., nTE))
    @inbounds for I in CartesianIndices(MatrixSize)
        M[I,:] .= T(1e5 + 1e5*rand()) .* noiseysignal()
    end

    Random.seed!(oldseed)
    return M
end

# Mock T2 distribution, computed with default parameters
function mock_T2_dist(o::T2mapOptions = mock_t2map_opts(Float64); kwargs...)
    T2mapSEcorr(mock_image(o; kwargs...), T2mapOptions(o; kwargs..., Silent = true))[2]
end


text:
export hessian_dense

function hessian_dense(nlp, x)
    return Matrix(hess(nlp, x))
    #return Symmetric(Matrix(hess(nlp, x)))
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
Base.@kwdef mutable struct RKstruct
  nStage = nothing
  ARKE = nothing
  bRKE = nothing
  ARKI = nothing
  bRKI = nothing
  gRKI = nothing
end
function RungeKuttaMethod(Method)
  str = Method
  RK = RKstruct(;)
if str == "ARK32"
    RK.nStage=3;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.bRKE=zeros(1,RK.nStage);
    RK.ARKE[2,1]=2-sqrt(2);
    ARKE32=1/6*(3 + sqrt(2));
    RK.ARKE[3,1]=1-ARKE32;
    RK.ARKE[3,2]=ARKE32;
    RK.bRKE[1]=1/(2*sqrt(2));
    RK.bRKE[2]=1/(2*sqrt(2));
    RK.bRKE[3]=1-1/(sqrt(2));
    RK.ARKI=zeros(RK.nStage,RK.nStage);
    RK.bRKI=zeros(1,RK.nStage);
    RK.ARKI[2,1]=1-1/(sqrt(2));
    RK.ARKI[2,2]=1-1/(sqrt(2));
    RK.ARKI[3,1]=1/(2*sqrt(2));
    RK.ARKI[3,2]=1/(2*sqrt(2));
    RK.ARKI[3,3]=1-1/(sqrt(2));
    RK.bRKI[1]=1/(2*sqrt(2));
    RK.bRKI[2]=1/(2*sqrt(2));
    RK.bRKI[3]=1-1/(sqrt(2));
    RK.gRKI=1-1/(sqrt(2));
elseif str == "Trap2(2,3,2)"
    #Trap2(2,3,2)
    RK.nStage=4;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.bRKE=zeros(1,RK.nStage);
    RK.ARKE[2,1]=1;
    RK.ARKE[3,1]=1/2;
    RK.ARKE[3,2]=1/2;
    RK.ARKE[4,1]=1/2;
    RK.ARKE[4,3]=1/2;
    RK.bRKE[1]=1/2;
    RK.bRKE[2]=0;
    RK.bRKE[3]=1/2;
    RK.bRKE[4]=0;
elseif str == "RK1"
    RK.nStage=1;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.bRKE=zeros(1,RK.nStage);
    RK.bRKE[1]=1;
elseif str == "RK1I"
    RK.nStage=1;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.ARKE[1,1]=1;
    RK.bRKE=zeros(1,RK.nStage);
    RK.bRKE[1]=1;
elseif str == "RK4"
    RK.nStage=4;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.ARKE[2,1]=1/2;
    RK.ARKE[3,2]=1/2;
    RK.ARKE[4,3]=1;
    RK.bRKE=zeros(1,RK.nStage);
    RK.bRKE[1]=1/6;
    RK.bRKE[2]=1/3;
    RK.bRKE[3]=1/3;
    RK.bRKE[4]=1/6;
elseif str == "RK3"
    RK.nStage=3;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.ARKE[2,1]=1/3;
    RK.ARKE[3,2]=1/2;
    RK.bRKE=zeros(1,RK.nStage);
    RK.bRKE[1]=0;
    RK.bRKE[2]=0;
    RK.bRKE[3]=1;
elseif str == "DBM453"
    RK.nStage=5;
    RK.ARKE=zeros(RK.nStage,RK.nStage);
    RK.bRKE=zeros(1,RK.nStage);
    RK.ARKE[2,1]=0.10306208811591838;
    RK.ARKE[3,1]=-0.94124866143519894;
    RK.ARKE[3,2]=1.6626399742527356;
    RK.ARKE[4,1]=-1.3670975201437765;
    RK.ARKE[4,2]=1.3815852911016873;
    RK.ARKE[4,3]=1.2673234025619065;
    RK.ARKE[5,1]=-0.81287582068772448;
    RK.ARKE[5,2]=0.81223739060505738;
    RK.ARKE[5,3]=0.90644429603699305;
    RK.ARKE[5,4]=0.094194134045674111;

    RK.bRKE[1]=0.87795339639076672;
    RK.bRKE[2]=-0.72692641526151549;
    RK.bRKE[3]=0.7520413715737272;
    RK.bRKE[4]=-0.22898029400415090;
    RK.bRKE[5]=0.32591194130117247;

    RK.gRKI=0.32591194130117247;
    RK.ARKI=zeros(RK.nStage,RK.nStage);
    RK.bRKI=zeros(1,RK.nStage);
    RK.ARKI[2,1]=-0.22284985318525410;
    RK.ARKI[2,2]=RK.gRKI;
    RK.ARKI[3,1]=-0.46801347074080545;
    RK.ARKI[3,2]=0.86349284225716961;
    RK.ARKI[3,3]=RK.gRKI;
    RK.ARKI[4,1]=-0.46509906651927421;
    RK.ARKI[4,2]=0.81063103116959553;
    RK.ARKI[4,3]=0.61036726756832357;
    RK.ARKI[4,4]=RK.gRKI;
    RK.ARKI[5,1]=0.87795339639076675;
    RK.ARKI[5,2]=-0.72692641526151547;
    RK.ARKI[5,3]=0.75204137157372720;
    RK.ARKI[5,4]=-0.22898029400415088;
    RK.ARKI[5,5]=RK.gRKI;
    RK.bRKI=RK.bRKE;
end
return RK
end


text:
function output_init(ps_m::CanonicalModel,
                    devices::Array{PSD,1},
                    parameters::Bool) where {PSD <: PSY.ThermalGen}

    initial_conditions = Vector{InitialCondition}(undef, length(devices))

    idx = eachindex(devices)
    i, state = iterate(idx)
    for g in devices
        if !isnothing(g.tech.ramplimits)
            if parameters
                initial_conditions[i] = PSI.InitialCondition(g, PJ.add_parameter(ps_m.JuMPmodel, g.tech.activepower))
            else
                initial_conditions[i] = PSI.InitialCondition(g, g.tech.activepower)
            end
            y = iterate(idx, state)
            y === nothing && (i += 1; break)
            i, state = y
        end
    end

    deleteat!(initial_conditions, i:last(idx))

    ps_m.initial_conditions[:thermal_output] = initial_conditions

    return

end

function status_init(ps_m::CanonicalModel,
                    devices::Array{PSD,1},
                    parameters::Bool) where {PSD <: PSY.ThermalGen}

    initial_conditions = Vector{InitialCondition}(undef, length(devices))

    idx = eachindex(devices)
    i, state = iterate(idx)
    for g in devices
        if !isnothing(g.tech.ramplimits)
            if parameters
                initial_conditions[i] = PSI.InitialCondition(g, PJ.add_parameter(ps_m.JuMPmodel, 1.0*(g.tech.activepower > 0)))
            else
                initial_conditions[i] = PSI.InitialCondition(g, 1.0*(g.tech.activepower > 0))
            end
            y = iterate(idx, state)
            y === nothing && (i += 1; break)
            i, state = y
        end
    end

    deleteat!(initial_conditions, i:last(idx))

    ps_m.initial_conditions[:thermal_status] = initial_conditions

    return

end

function duration_init(ps_m::CanonicalModel,
                        devices::Array{PSD,1},
                        parameters::Bool) where {PSD <: PSY.ThermalGen}

    ini_cond_on = Vector{InitialCondition}(undef, length(devices))
    ini_cond_off = Vector{InitialCondition}(undef, length(devices))

    idx = eachindex(devices)
    i, state = iterate(idx)
    for g in devices
        if !isnothing(g.tech.ramplimits)
            if parameters
                ini_cond_on[i] = PSI.InitialCondition(g, PJ.add_parameter(ps_m.JuMPmodel, 999.0*(g.tech.activepower > 0)))
                ini_cond_off[i] = PSI.InitialCondition(g, PJ.add_parameter(ps_m.JuMPmodel, 999.0*(g.tech.activepower < 0)))
            else
                ini_cond_on[i] = PSI.InitialCondition(g, 999.0*(g.tech.activepower > 0))
                ini_cond_off[i] = PSI.InitialCondition(g, 999.0*(g.tech.activepower < 0))
            end
            y = iterate(idx, state)
            y === nothing && (i += 1; break)
            i, state = y
        end
    end

    deleteat!(ini_cond_on, i:last(idx))
    deleteat!(ini_cond_off, i:last(idx))

    ps_m.initial_conditions[:thermal_duration_on] = ini_cond_on
    ps_m.initial_conditions[:thermal_duration_off] = ini_cond_off

    return

end


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

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

module RationalTrees

using Nemo, NumberTheory, Counts

export ModuleRationalTrees
export EuclidTree, CalkinWilfTree, SchinzelSierpinskiEncoding

"""

Rational trees as understood here are binary trees enumerating the positive or
nonnegative rational numbers. Examples are the Euclid tree, the Kepler tree and the
Stern-Brocot tree (a.k.a. Farey tree). They are closely related to binary partitions
and to Stern's diatomic sequence or Dijkstra's fusc function.

Malter, Schleicher, Zagier, [New looks at old number theory](https://pdfs.semanticscholar.org/6d28/dcef911dd91f47e6ca4bd2c564c1f3099a05.pdf), Amer. Math. Monthly, 120(3), 2013, pp. 243-264.

* EuclidTree, CalkinWilfTree, SchinzelSierpinskiEncoding
"""
const ModuleRationalTrees = ""

# http://oeis.org/wiki/User:Peter_Luschny/SchinzelSierpinskiConjectureAndCalkinWilfTree
# A294442 A294446

"""

```
julia> for n ∈ 1:4 Println(EuclidTree(n)) end
[1//1]
[1//2, 2//1]
[1//3, 3//2, 2//3, 3//1]
[1//4, 4//3, 3//5, 5//2, 2//5, 5//3, 3//4, 4//1]
```
"""
function EuclidTree(n)

    function DijkstraFusc(m)
        a, b, k = 1, 0, m
        while k > 0
            k % 2 == 1 ? b += a : a += b
            k = k >> 1
        end
        b
    end

    DF = [DijkstraFusc(k) for k ∈ 2^(n-1):2^n]
    [fmpq(DF[j], DF[j+1]) for j ∈ 1:2^(n-1)]
end

"""

Alias for the (much better named) EuclidTree. See Malter, Schleicher, Zagier, New looks at old number theory, Amer. Math. Monthly, 120(3), 2013, pp. 243-264.
"""
CalkinWilfTree(n) = EuclidTree(n)

"""

Return the Schinzel-Sierpinski encoding of the positive rational number r.

```
julia> for n ∈ 1:4 println([SchinzelSierpinski(l) for l ∈ EuclidTree(n)]) end
[1//1]
[2//5, 5//2]
[3//11, 5//3, 3//5, 11//3]
[2//11, 3//2, 11//19, 19//7, 7//19, 19//11, 2//3, 11//2]
```
"""
function SchinzelSierpinskiEncoding(l, searchLimit = 500000)

    a, b = numerator(l), denominator(l)
    sgn = a < b ? -1 : 1
    p, q = 1, 2

    while q < searchLimit
        r = a * (q + 1)
        r % b == 0 && (p = div(r, b) - 1)
        isPrime(p) && return fmpq(p, q)
        q = NextPrime(q)
    end
    warn("Search limit reached for ", l)
    return 0
end

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

using Test, SeqUtils, PrimesIterator

function test()
    @testset "EuclidTree" begin
        S = [numerator(sum(r for r ∈ EuclidTree(n))) for n ∈ 1:9]
        L = fmpz[1, 5, 11, 23, 47, 95, 191, 383, 767] # A052940
        @test all(S[1:9] .== L[1:9])
    end
end

function demo()
    println("\nEuclidTree")
    for n ∈ 1:5
        Println(EuclidTree(n))
    end
    println()

    println("\nSchinzelSierpinskiEncoding of the Euclid tree")
    for n ∈ 1:5
        Println([SchinzelSierpinskiEncoding(l) for l ∈ EuclidTree(n)])
        # println(sum([SchinzelSierpinskiEncoding(l) for l ∈ EuclidTree(n)]))
    end
    println()

    println("\nSchinzelSierpinskiEncoding of primes")
    Println([SchinzelSierpinskiEncoding(fmpq(p, 1)) for p ∈ Primes(100)])
    println()
end

"""

EuclidTree(100) ::
    0.000005 seconds (4 allocations: 224 bytes)
[SchinzelSierpinskiEncoding(l) for l ∈ EuclidTree(100)] ::
    0.000007 seconds (6 allocations: 320 bytes)
"""
function perf()
    GC.gc()
    @time EuclidTree(100)
    GC.gc()
    @time [SchinzelSierpinskiEncoding(l) for l ∈ EuclidTree(100)]
end

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

main()

end # module

#=
# A062251 := n -> SchinzelSierpinskiPrimes(n,"numer"):
# A062251 = 2, 5, 11, 11, 19, 17, 41, 23, 53, 29, 43, 47, 103, 41,
# 59, 47, 67, 53, 113, 59, 83, 131, 137, 71, 149, 103, 107, 83,
# A060324 := n -> SchinzelSierpinskiPrimes(n,"denom"):
# A060324 = 2, 2, 3, 2, 3, 2, 5, 2, 5, 2, 3, 3, 7, 2, 3, 2,
# 3, 2, 5, 2, 3, 5, 5, 2, 5, 3, 3, 2, 5, 2, 13, 3
=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
include("Flux.jl")
include("CuArrays.jl")
include("Zygote.jl")
include("ReinforcementLearningBase.jl")


text:
const ALL_COLUMNS = (
    "auxHeat1",
    "auxHeat2",
    "auxHeat3",
    "compCool1",
    "compCool2",
    "compHeat1",
    "compHeat2",
    "dehumidifier",
    "dmOffset",
    "economizer",
    "fan",
    "humidifier",
    "hvacMode",
    "outdoorHumidity",
    "outdoorTemp",
    "sky",
    "ventilator",
    "wind",
    "zoneAveTemp",
    "zoneCalendarEvent",
    "zoneClimate",
    "zoneCoolTemp",
    "zoneHeatTemp",
    "zoneHumidity",
    "zoneHumidityHigh",
    "zoneHumidityLow",
    "zoneHvacMode",
    "zoneOccupancy"
  )

  const DEFAULT_COLUMNS = (
      "compCool1",
      "compHeat1",
      "fan",
      "hvacMode",
      "outdoorHumidity",
      "outdoorTemp",
      "zoneAveTemp",
      "zoneClimate",
      "zoneCoolTemp",
      "zoneHeatTemp",
      "zoneHumidity",
      "zoneHvacMode",
      "zoneOccupancy"
  )

const COLUMN_TITLE_LOOKUP = Dict(
    "auxHeat1" => "Aux Heat (sec)",
    "auxHeat2" => "Aux Heat Stage 2 (sec)",
    "auxHeat3" => "Aux Heat Stage 3 (sec)",
    "compCool1" => "Cool Stage 1 (sec)",
    "compCool2" => "Cool Stage 2 (sec)",
    "compHeat1" => "Heat Stage 1 (sec)",
    "compHeat2" => "Heat Stage 2 (sec)",
    "dehumidifier" => "Dehumidifier (sec)",
    "dmOffset" => "Demand Mgmt Offset (F)",
    "economizer" => "Economizer Runtime (sec)",
    "fan" => "Fan (sec)",
    "humidifier" => "Humidifier (sec)",
    "hvacMode" => "HVAC Mode",
    "outdoorHumidity" => "Outdoor Humidity (%)",
    "outdoorTemp" => "Outdoor Temp (F)",
    "sky" => "Sky Cover",
    "ventilator" => "Ventilator (sec)",
    "wind" => "Wind Speed (km/h)",
    "zoneAveTemp" => "Indoor Temp Avg (F)",
    "zoneCalendarEvent" => "Override Event",
    "zoneClimate" => "Climate Mode",
    "zoneCoolTemp" => "Zone Cool Temp",
    "zoneHeatTemp" => "Zone Heat Temp",
    "zoneHumidity" => "Humidity Avg (%)",
    "zoneHumidityHigh" => "Humidity High (%)",
    "zoneHumidityLow" => "Humidity Low (%)",
    "zoneHvacMode" => "HVAC System Mode",
    "zoneOccupancy" => "Zone Occupancy"
)


text:
module Absorption

using Parameters                # For constructing HitranTable with keywords
using DocStringExtensions       # For simplifying docstring
using Interpolations            # For interpolating in both lookup tables and qoft!
using JLD2                      # For saving and loading the interpolator
using ProgressMeter             # For showing progress, especially in creating interpolator
using KernelAbstractions        # For heterogeneous (GPU+CPU) programming
using CUDA                      # For GPU programming
using ForwardDiff, DiffResults  # For auto-differentiation
using ..Architectures           # For GPU/CPU convenience
using NetCDF                    # For loading NetCDF files with constants

include("Constants/constants.jl")                   # Scientific and mathematical constants
include("Types/types.jl")                           # All types used in this module
include("Hitran/hitran.jl")                         # HITRAN file-related functions
include("LineShapes/complex_error_functions.jl")    # CEFs used in line broadening
include("Constants/mol_weights.jl")                 # Molecular weights
include("Constants/TIPS_2017.jl")                   # Partition sums data
include("TableInterpolation/partition_sums.jl")     # Partition sums interpolator
include("TableInterpolation/cross_section_interpolator.jl") # CS interpolator functions
include("Hitran/absorption_cross_section.jl")       # Cross-section from HITRAN
include("Hitran/cross_section_autodiff.jl")         # Auto-differentiation

# Export the Cross Section models
export AbstractCrossSectionModel, HitranModel, InterpolationModel

# Export the absorption_cross_section functions
export compute_absorption_cross_section, absorption_cross_section

# Export the hitran functions from hitran.jl
export read_hitran, make_hitran_model

# Export the broadening function types
export Doppler, Lorentz, Voigt, AbstractBroadeningFunction

# Export the complex error functions
export AbstractComplexErrorFunction, HumlicekErrorFunction,
       HumlicekWeidemann32VoigtErrorFunction, HumlicekWeidemann32SDErrorFunction,
       CPF12ErrorFunction, ErfcHumliErrorFunctionVoigt, ErfcHumliErrorFunctionSD,
       ErfcErrorFunction

# Export the hitran table struct type
export HitranTable

# Export the interpolator functions
export make_interpolation_model, save_interpolation_model, load_interpolation_model

end


text:
module HarmonicPolya;

using LinearAlgebra
using MultivariatePolynomials
using DynamicPolynomials
using FastGaussQuadrature
using SpecialFunctions
using PolynomialBases
using FixedPolynomials

const MP = MultivariatePolynomials
const DP = DynamicPolynomials

import LinearAlgebra: length

include("HarmonicBases.jl")
include("SphericalQuadrature.jl")
include("SquaresOperator.jl")
include("FawziFangOperator.jl")
include("SphereMinimization.jl")

export harmonicdecomposition
export sphericalquadrature
export upperbound
export lowerboundfawzi
export lowerboundsquares

end


text:
struct analysis_lm_Gauss_NIP_Optim
    formula::FormulaTerm
    modelClass::String
    LikelihoodMod::String
    PriorMod::String
    Link::String
    ComputeMethod::String
    res
    fit
    LogLike::Float64
    AIC::Float64
    BIC::Float64
    R_sqr::Float64
    Adjusted_R_sqr::Float64
    sigma::Float64
    fittedResponse
    residuals
    Cooks_distance
end

function linear_reg(formula::FormulaTerm,data::DataFrame)

      modelClass = "LinearReg";
      LikelihoodMod="Gauss";
      PriorMod="NIP";
      Link = "Identity";
      ComputeMethod="Optimization";
      formla = formula;
      formula = apply_schema(formula, schema(formula, data));
      y, X = modelcols(formula, data);
      fm_frame=ModelFrame(formula,data);
      X=modelmatrix(fm_frame);
      p = size(X, 2);
      n = size(X, 1);

      res = lm(formula,data);
      fit = coeftable(res)
      R_sqr = GLM.r2(res);
      Adjusted_R_sqr = GLM.adjr2(res);

      logLike = GLM.loglikelihood(res)
      npar =(p+1);
      AIC = 2*npar - 2*logLike
      BIC = log(n)*npar - 2*logLike

      sse = GLM.deviance(res)
      degrees_of_freedon = GLM.dof_residual(res)
      mse = sse/degrees_of_freedon
      sigma = sqrt(mse)

      fittedResponse = predict(res);
      resid = residuals(res);

      cooks_distance = GLM.StatsBase.cooksdistance(res);
      ans = analysis_lm_Gauss_NIP_Optim(formla,modelClass,LikelihoodMod,PriorMod
      ,Link,ComputeMethod,res,fit, logLike, AIC, BIC, R_sqr
      ,Adjusted_R_sqr,sigma, fittedResponse,resid,cooks_distance)

      return ans

end

function linear_reg_predicts(obj,newdata::DataFrame)

    formula = obj.formula;
    fm_frame=ModelFrame(formula,newdata);
    X=modelmatrix(fm_frame);
    beta = coef(obj.res)
    y_pred = X*beta;
    y_pred
end

function linear_reg(formula::FormulaTerm,data::DataFrame,PriorMod::Prior_Ridge,h::Float64,sim_size::Int64)

    formula = apply_schema(formula, schema(formula, data));
    y, X = modelcols(formula, data);

    @model LinReg(X, y) = begin
      p=size(X, 2);

      #priors
      a0=0.1
      b0=0.1

      v ~ InverseGamma(h,h)
      σ ~ InverseGamma(a0,b0)
      α ~ Normal(0,v*σ)
      β ~ filldist(Normal(0,v*σ), p)

      #likelihood
      y ~ MvNormal(α .+ X * β, σ);
    end;
    LinReg_model=LinReg(X,y);
    chain = sample(CRRao_rng, LinReg_model, NUTS(), sim_size);
    summaries, quantiles = describe(chain);

    ans = MCMC_chain(chain,summaries,quantiles)
    ans
end

function linear_reg(formula::FormulaTerm,data::DataFrame,PriorMod::Prior_Laplace,h::Float64,sim_size::Int64)

    formula = apply_schema(formula, schema(formula, data));
    y, X = modelcols(formula, data);

    @model LinReg(X, y) = begin
      p=size(X, 2);

      #priors
      a0=0.1
      b0=0.1
      v ~ InverseGamma(h,h)
      σ ~ InverseGamma(a0,b0)
      α ~ Laplace(0,σ*v)
      β ~ filldist(Laplace(0,σ*v), p)

      #likelihood
      y ~ MvNormal(α .+ X * β, σ);
    end;
    LinReg_model=LinReg(X,y);
    chain = sample(CRRao_rng, LinReg_model, NUTS(), sim_size);
    summaries, quantiles = describe(chain);

    ans = MCMC_chain(chain,summaries,quantiles)
    ans

end

function linear_reg(formula::FormulaTerm,data::DataFrame,PriorMod::Prior_Cauchy,sim_size::Int64)

    formula = apply_schema(formula, schema(formula, data));
    y, X = modelcols(formula, data);

    @model LinReg(X, y) = begin
      p=size(X, 2);

      #priors
      σ ~ Truncated(TDist(1),0,Inf)
      α ~ TDist(1)*σ
      β ~ filldist(TDist(1)*σ, p)

      #likelihood
      y ~ MvNormal(α .+ X * β, σ);
    end;
    LinReg_model=LinReg(X,y);
    chain = sample(CRRao_rng, LinReg_model, NUTS(), sim_size);
    summaries, quantiles = describe(chain);

    ans = MCMC_chain(chain,summaries,quantiles)
    ans

end

function linear_reg(formula::FormulaTerm,data::DataFrame,PriorMod::Prior_TDist,h::Float64,sim_size::Int64)

    formula = apply_schema(formula, schema(formula, data));
    y, X = modelcols(formula, data);

    @model LinReg(X, y) = begin
      p=size(X, 2);

      #priors
      a0=0.1
      b0=0.1
      ν ~ InverseGamma(h,h)
      σ ~ InverseGamma(a0,b0)
      α ~ TDist(ν)*σ
      β ~ filldist(TDist(ν)*σ, p)

      #likelihood
      y ~ MvNormal(α .+ X * β, σ);
    end;
    LinReg_model=LinReg(X,y);
    chain = sample(CRRao_rng, LinReg_model, NUTS(), sim_size);
    summaries, quantiles = describe(chain);

    ans = MCMC_chain(chain,summaries,quantiles)
    ans

end

function linear_reg(formula::FormulaTerm,data::DataFrame,PriorMod::Prior_Uniform,h::Float64,sim_size::Int64)

    formula = apply_schema(formula, schema(formula, data));
    y, X = modelcols(formula, data);

    @model LinReg(X, y) = begin
      p=size(X, 2);

      #priors
      v=1/h;
      σ ~ Uniform(0,v)
      α ~ Uniform(-v*σ,v*σ)
      β ~ filldist(Uniform(-v,v), predictors)

      #likelihood
      y ~ MvNormal(α .+ X * β, σ);
    end;
    LinReg_model=LinReg(X,y);
    chain = sample(CRRao_rng, LinReg_model, NUTS(), sim_size);
    summaries, quantiles = describe(chain);

    ans = MCMC_chain(chain,summaries,quantiles)
    ans

end