text:
"""
    AeroState(alpha, beta, alpha_dot, tas, eas, cas, ias, qinf, mach)

Aerodynamic state struct storing relevant aerodynamic data.

# Constructors

Different constructors are provided:

    AeroState(alpha, beta, alpha_dot, tas, eas, cas, ias, qinf, mach)
    AeroState(state::State, env::Environment)
    AeroState(tas::Number, alpha::Number, beta::Number, height::Number)

"""
struct AeroState
    alpha :: Number  # AOA [rad]
    beta :: Number  # AOS [rad]

    alpha_dot::Number  # rad/s

    tas :: Number  # True Airspeed [m/s]
    eas :: Number  # Equivalent Airspeed [m/s]
    cas :: Number  # Calibrated Airspeed [m/s]
    ias :: Number  # Indicated Airspeed [m/s]

    qinf :: Number  # Dynamic Pressure [Pa]
    mach :: Number  # Mach number
end

function AeroState(state::State, env::Environment)

    p = get_pressure(env)
    ρ = get_density(env)
    a = get_sound_velocity(env)

    # coming from
    wind_hor = get_wind_NED(env)
    wind_body = hor2body(wind_hor..., get_euler_angles(state)...)

    aero_speed = get_body_velocity(state) + wind_body

    tas, alpha, beta = tas_alpha_beta_from_uvw(aero_speed...)

    cas = tas2cas(tas, ρ, p)
    ias = tas
    eas = tas2eas(tas, ρ)

    qinf = incompressible_qinf(tas, ρ)
    mach = tas / a

    AeroState(alpha, beta, 0, tas, eas, cas, ias, qinf, mach)
end

function AeroState(tas::Number, alpha::Number, beta::Number, height::Number)

    T, p, ρ, a = atmosphere_isa(height)

    cas = tas2cas(tas, ρ, p)
    ias = tas
    eas = tas2eas(tas, ρ)

    qinf = incompressible_qinf(tas, ρ)
    mach = tas / a

    AeroState(alpha, beta, 0, tas, eas, cas, ias, qinf, mach)
end

"""
    generate_state_aerostate(pos::Position, att::Attitude, tas::Number, alpha::Number,
                                beta::Number, env::Environment,
                                ang_vel=[0.0, 0.0, 0.0], accel=[0.0, 0.0, 0.0],
                                ang_accel=[0.0, 0.0, 0.0])

returns both state and aerostate provided the values of aircrft states, environment and
aerodynamic angles.
"""
function generate_state_aerostate(pos::Position, att::Attitude, tas::Number, alpha::Number,
                         beta::Number, env::Environment,
                         ang_vel=[0.0, 0.0, 0.0], accel=[0.0, 0.0, 0.0],
                         ang_accel=[0.0, 0.0, 0.0])

    aerostate = AeroState(tas, alpha, beta, get_height(pos))

    aero_vel_body = wind2body([tas, 0, 0]..., alpha, beta)
    wind_vel_body = hor2body(get_wind_NED(env)..., get_euler_angles(att)...)
    vel = aero_vel_body - wind_vel_body

    state = State(pos, att, vel, ang_vel, accel, ang_accel)
    return state, aerostate
end

"""
    check_state_aerostate_env_coherence(state::State, aerost::AeroState,
                                        env::Environment)

Check if aerodynamic velocity in aerostate is coherent with inertial velocity in
state and wind velocity in environment.
"""
function check_state_aerostate_env_coherence(state::State, aerost::AeroState,
                                             env::Environment)

    att = get_attitude(state)
    # wind vel from environment
    wind = get_wind_NED(env)
    # ground speed from state
    vel = get_horizon_velocity(state)
    # Composition (wind posititve blowing towards)
    aero_vel_estimation = vel + wind

    # Transform aerostate velocity to horizon
    tas, α, β = get_tas(aerost), get_alpha(aerost), get_beta(aerost)
    aero_vel = body2hor(wind2body(tas, 0, 0, α, β)..., get_euler_angles(state)...)

    isapprox(aero_vel_estimation, aero_vel)
end

get_alpha(aerost::AeroState) = aerost.alpha
get_beta(aerost::AeroState) = aerost.beta
get_aero_angles(aerost::AeroState) = [aerost.alpha, aerost.beta]
get_alpha_dot(aerost::AeroState) = aerost.alpha_dot

get_aero_speeds(aerost::AeroState) = [aerost.tas, aerost.eas, aerost.cas, aerost.ias]
get_tas(aerost::AeroState) = aerost.tas
get_eas(aerost::AeroState) = aerost.eas
get_cas(aerost::AeroState) = aerost.cas
get_ias(aerost::AeroState) = aerost.ias

get_qinf(aerost::AeroState) = aerost.qinf
get_mach(aerost::AeroState) = aerost.mach


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

using Markdown
using InteractiveUtils

# ╔═╡ 8ab9b70a-e98d-11ea-239c-73dc659722c2
using YaoLang, YaoPlots, ZXCalculus, Compose, Revise

# ╔═╡ 282bd90e-e9fd-11ea-0b79-174cade05905
using Yao, YaoExtensions

# ╔═╡ e0a83244-e98e-11ea-3213-fdf02f4df920
function Base.show(io::IO, mime::MIME"text/html", zx::Union{ZXDiagram, ZXGraph})
	g = plot(zx; size_x=13, size_y = 3, NODELABELSIZE=8pt)
	Compose.set_default_graphic_size(18cm, 5cm)
	Base.show(io, mime, compose(context(0.0, 0.0, 1.0, 1.0), g))
end

# ╔═╡ 227f7884-e99a-11ea-3a90-0beb697a2da6
md"# Construct a ZX diagram"

# ╔═╡ 49d6e6ac-e994-11ea-2ac5-27ab8242e297
z1 = ZXDiagram(4)

# ╔═╡ 51e72d20-e994-11ea-1a50-854039f728aa
push_gate!(z1, Val(:Z), 1, 3//2)

# ╔═╡ e1dbb828-e995-11ea-385d-fb20b58d1b49
let
	push_gate!(z1, Val(:H), 1)
	push_gate!(z1, Val(:H), 1)
	push_gate!(z1, Val(:Z), 1, 1//2)
	push_gate!(z1, Val(:H), 4)
	push_gate!(z1, Val(:CZ), 4, 1)
	push_gate!(z1, Val(:CNOT), 1, 4)
end

# ╔═╡ 60c59c0a-e994-11ea-02da-7360cbcf81f7
function load_graph()
	zxd = ZXDiagram(4)
	push_gate!(zxd, Val(:Z), 1, 3//2)
	push_gate!(zxd, Val(:H), 1)
    push_gate!(zxd, Val(:Z), 1, 1//2)
    push_gate!(zxd, Val(:H), 4)
    push_gate!(zxd, Val(:CZ), 4, 1)
    push_gate!(zxd, Val(:CNOT), 1, 4)
    push_gate!(zxd, Val(:H), 1)
    push_gate!(zxd, Val(:H), 4)
    push_gate!(zxd, Val(:Z), 1, 1//4)
    push_gate!(zxd, Val(:Z), 4, 3//2)
    push_gate!(zxd, Val(:X), 4, 1//1)
    push_gate!(zxd, Val(:H), 1)
    push_gate!(zxd, Val(:Z), 4, 1//2)
    push_gate!(zxd, Val(:X), 4, 1//1)
    push_gate!(zxd, Val(:Z), 2, 1//2)
    push_gate!(zxd, Val(:CNOT), 3, 2)
    push_gate!(zxd, Val(:H), 2)
    push_gate!(zxd, Val(:CNOT), 3, 2)
    push_gate!(zxd, Val(:Z), 2, 1//4)
    push_gate!(zxd, Val(:Z), 3, 1//2)
    push_gate!(zxd, Val(:H), 2)
    push_gate!(zxd, Val(:H), 3)
    push_gate!(zxd, Val(:Z), 3, 1//2)
    push_gate!(zxd, Val(:CNOT), 3, 2)
end

# ╔═╡ 90b83d5e-e99a-11ea-1fb2-95c907668262
md"# Simplify the ZX diagram"

# ╔═╡ 64bff9ec-e9b5-11ea-3b23-c51d2149697a
zxd = load_graph()

# ╔═╡ 5dbf9f96-e9a4-11ea-19d7-e15e7f2327c9
tcount(zxd)

# ╔═╡ db9a0d4e-e99e-11ea-22ab-1fead216dd07
zxg = ZXGraph(zxd)

# ╔═╡ 66eb6e1a-e99f-11ea-141c-a9017390524f
md"apply the `lc` rule recursively"

# ╔═╡ ce83c0be-e9a3-11ea-1a40-b1b5118a24bd
html"""

"""

# ╔═╡ a6b92942-e99a-11ea-227d-f9fe53f8a1cf
# simplify!(Rule{:lc}(), zxd)  #  this should not pass! use `DRule` and `GRule` to distinguish them?
simplify!(Rule{:lc}(), zxg)  # allow Rule(:lc) for simplicity.

# ╔═╡ 86475062-e99f-11ea-2f44-a3c270cc45e5
md"apply the p1 rule recursively"

# ╔═╡ f59bf644-e9a8-11ea-1944-b3843ef5d6c8
html"""

"""

# ╔═╡ b739540e-e99a-11ea-2a04-abd99889cf92
simplify!(Rule{:p1}(), zxg)  # does not have any effect?

# ╔═╡ 7af70558-e9b4-11ea-3aa9-3b73357f0a2a
srule!(sym::Symbol) = g -> simplify!(Rule{sym}(), g)

# ╔═╡ a5784394-e9b4-11ea-0e68-8d8211766409
srule_once!(sym::Symbol) = g -> replace!(Rule{sym}(), g)

# ╔═╡ c6f809e8-e9b4-11ea-2dcb-57c4a1d65bb7
zxg |> srule!(:lc) |> srule!(:p1) |> srule_once!(:pab)

# ╔═╡ 25d876b6-e9a9-11ea-2631-fd6f8934daa6
md"apply the `pab` rule once"

# ╔═╡ 3ce5329a-e9a9-11ea-2c7e-312416dd9483
html"""

"""

# ╔═╡ bd2b3364-e99a-11ea-06e7-4560cb873d2c
replace!(Rule{:pab}(), zxg)  # this naming is not explict, what about `simplify_recursive!` and `simplily!`.

# ╔═╡ c71cdf4c-e9b5-11ea-2aaf-5f4be0eb3e93
md"## To make life easier"

# ╔═╡ b6eaa762-e9b5-11ea-145e-7b5fa6b01872
zxd2 = load_graph()

# ╔═╡ 5464bc40-e9b5-11ea-2a2e-dfb0d31a33d3
ex_zxd = clifford_simplification(zxd2)

# ╔═╡ c0f046ea-e9b5-11ea-302c-c5fd4399d400
pt_zxd = phase_teleportation(zxd2)

# ╔═╡ c3e8b5b4-e99a-11ea-0e56-6b18757f94df
md"# Extract circuit"

# ╔═╡ c1b9828c-e99a-11ea-006a-013a2eab8cf3
zxd_simplified = circuit_extraction(zxg)

# ╔═╡ 6ddffea2-e9a4-11ea-1c32-0975a45aba7c
tcount(zxd_simplified)

# ╔═╡ 2082486e-e9fd-11ea-1a46-6395b4b34657
md"# Porting Yao"

# ╔═╡ 3a2c5d22-e9fd-11ea-1e48-a9be2e8b213b
zxd3 = ZXDiagram(4)

# ╔═╡ d341cde4-e9fd-11ea-1c3b-63f35fc648d5
begin
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::AbstractBlock)
	error("Block type `$c` is not supported.")
end
# rotation blocks
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,RotationGate{1,T,XGate}}) where {N,T}
	push_gate!(zxd, Val(:X), c.locs[1], c.content.theta)
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,RotationGate{1,T,ZGate}}) where {N,T}
	push_gate!(zxd, Val(:Z), c.locs[1], c.content.theta)
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,ShiftGate{T}}) where {N,T}
	push_gate!(zxd, Val(:Z), c.locs[1], c.content.theta)
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,HGate}) where {N}
	push_gate!(zxd, Val(:H), c.locs[1])
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::ChainBlock{N}) where {N}
	push_gate!.(Ref(zxd), subblocks(c))
	zxd
end

# constant block
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,RotationGate{1,T,XGate}}) where {N,T}
	push_gate!(zxd, Val(:X), c.locs[1], c.content.theta)
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,RotationGate{1,T,ZGate}}) where {N,T}
	push_gate!(zxd, Val(:Z), c.locs[1], c.content.theta)
end

# control blocks
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::ControlBlock{N,1,RotationGate{1,T,XGate}}) where {N,T}
	push_gate!(zxd, Val(:X), c.locs[1], c.content.theta)
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::ControBlock{N,1,RotationGate{1,T,ZGate}}) where {N,T}
	push_gate!(zxd, Val(:Z), c.locs[1], c.content.theta)
end
function ZXCalculus.push_gate!(zxd::ZXDiagram, c::ControlBlock{N,1,ShiftGate{T}}) where {N,T}
	push_gate!(zxd, Val(:Z), c.locs[1], c.content.theta)
end
end

# ╔═╡ 581e847c-e9fd-11ea-3fd0-6bbc0f6efd56
c = qft_circuit(4)

# ╔═╡ 2e84a5ce-e9fd-11ea-12d0-b3a3dd75a76f
push_gate!(zxd3, c)

# ╔═╡ 7b850816-ea02-11ea-183c-5db2d670be24
ZXDiagram(4) |> typeof

# ╔═╡ Cell order:
# ╠═8ab9b70a-e98d-11ea-239c-73dc659722c2
# ╟─e0a83244-e98e-11ea-3213-fdf02f4df920
# ╟─227f7884-e99a-11ea-3a90-0beb697a2da6
# ╠═49d6e6ac-e994-11ea-2ac5-27ab8242e297
# ╠═51e72d20-e994-11ea-1a50-854039f728aa
# ╠═e1dbb828-e995-11ea-385d-fb20b58d1b49
# ╠═60c59c0a-e994-11ea-02da-7360cbcf81f7
# ╟─90b83d5e-e99a-11ea-1fb2-95c907668262
# ╠═64bff9ec-e9b5-11ea-3b23-c51d2149697a
# ╠═5dbf9f96-e9a4-11ea-19d7-e15e7f2327c9
# ╠═db9a0d4e-e99e-11ea-22ab-1fead216dd07
# ╟─66eb6e1a-e99f-11ea-141c-a9017390524f
# ╟─ce83c0be-e9a3-11ea-1a40-b1b5118a24bd
# ╠═a6b92942-e99a-11ea-227d-f9fe53f8a1cf
# ╟─86475062-e99f-11ea-2f44-a3c270cc45e5
# ╟─f59bf644-e9a8-11ea-1944-b3843ef5d6c8
# ╠═b739540e-e99a-11ea-2a04-abd99889cf92
# ╠═7af70558-e9b4-11ea-3aa9-3b73357f0a2a
# ╠═a5784394-e9b4-11ea-0e68-8d8211766409
# ╠═c6f809e8-e9b4-11ea-2dcb-57c4a1d65bb7
# ╟─25d876b6-e9a9-11ea-2631-fd6f8934daa6
# ╟─3ce5329a-e9a9-11ea-2c7e-312416dd9483
# ╠═bd2b3364-e99a-11ea-06e7-4560cb873d2c
# ╟─c71cdf4c-e9b5-11ea-2aaf-5f4be0eb3e93
# ╠═b6eaa762-e9b5-11ea-145e-7b5fa6b01872
# ╠═5464bc40-e9b5-11ea-2a2e-dfb0d31a33d3
# ╠═c0f046ea-e9b5-11ea-302c-c5fd4399d400
# ╟─c3e8b5b4-e99a-11ea-0e56-6b18757f94df
# ╠═c1b9828c-e99a-11ea-006a-013a2eab8cf3
# ╠═6ddffea2-e9a4-11ea-1c32-0975a45aba7c
# ╠═2082486e-e9fd-11ea-1a46-6395b4b34657
# ╠═282bd90e-e9fd-11ea-0b79-174cade05905
# ╠═3a2c5d22-e9fd-11ea-1e48-a9be2e8b213b
# ╠═d341cde4-e9fd-11ea-1c3b-63f35fc648d5
# ╠═581e847c-e9fd-11ea-3fd0-6bbc0f6efd56
# ╠═2e84a5ce-e9fd-11ea-12d0-b3a3dd75a76f
# ╠═7b850816-ea02-11ea-183c-5db2d670be24

ext: jl


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

module TestAddmul

using Base: rtoldefault
using Test
using LinearAlgebra
using LinearAlgebra: AbstractTriangular
using SparseArrays
using Random

_rand(::Type{T}) where {T <: AbstractFloat} = T(randn())
_rand(::Type{T}) where {F, T <: Complex{F}} = T(_rand(F), _rand(F))
_rand(::Type{T}) where {T <: Integer} =
    T(rand(max(typemin(T), -10):min(typemax(T), 10)))
_rand(::Type{BigInt}) = BigInt(_rand(Int))

function _rand(A::Type{<:Array}, shape)
    T = eltype(A)
    data = T[_rand(T) for _ in 1:prod(shape)]
    return copy(reshape(data, shape))
end

constructor_of(::Type{T}) where T = getfield(parentmodule(T), nameof(T))

function _rand(A::Type{<: AbstractArray}, shape)
    data = _rand(Array{eltype(A)}, shape)
    T = constructor_of(A)
    if A <: Union{Bidiagonal, Hermitian, Symmetric}
        return T(data, rand([:U, :L]))
        # Maybe test with both :U and :L?
    end
    return T(data)
end

_rand(A::Type{<: SymTridiagonal{T}}, shape) where {T} =
    SymTridiagonal(_rand(Symmetric{T}, shape))

const FloatOrC = Union{AbstractFloat, Complex{<: AbstractFloat}}
const IntegerOrC = Union{Integer, Complex{<: Integer}}
const LTri = Union{LowerTriangular, UnitLowerTriangular, Diagonal}
const UTri = Union{UpperTriangular, UnitUpperTriangular, Diagonal}

needsquare(::Type{<:Matrix}) = false
needsquare(::Type) = true

testdata = []

sizecandidates = 1:4
floattypes = [
    Float64, Float32, ComplexF64, ComplexF32,  # BlasFloat
    BigFloat,
]
inttypes = [
    Int,
    BigInt,
]
# `Bool` can be added to `inttypes` but it's hard to handle
# `InexactError` bug that is mentioned in:
# https://github.com/JuliaLang/julia/issues/30094#issuecomment-440175887
alleltypes = [floattypes; inttypes]
celtypes = [Float64, ComplexF64, BigFloat, Int]

mattypes = [
    Matrix,
    Bidiagonal,
    Diagonal,
    Hermitian,
    LowerTriangular,
    SymTridiagonal,
    Symmetric,
    Tridiagonal,
    UnitLowerTriangular,
    UnitUpperTriangular,
    UpperTriangular,
]

isnanfillable(::AbstractArray) = false
isnanfillable(::Array{<:AbstractFloat}) = true
isnanfillable(A::AbstractArray{<:AbstractFloat}) = parent(A) isa Array

"""
Sample `n` elements from `S` on average but make sure at least one
element is sampled.
"""
function sample(S, n::Real)
    length(S) <= n && return S
    xs = randsubseq(S, n / length(S))
    return length(xs) > 0 ? xs : rand(S, 1)  # sample at least one
end

function inputeltypes(celt, alleltypes = alleltypes)
    # Skip if destination type is "too small"
    celt <: Bool && return []
    filter(alleltypes) do aelt
        celt <: Real && aelt <: Complex && return false
        !(celt <: BigFloat) && aelt <: BigFloat && return false
        !(celt <: BigInt) && aelt <: BigInt && return false
        celt <: IntegerOrC && aelt <: FloatOrC && return false
        if celt <: IntegerOrC && !(celt <: BigInt)
            typemin(celt) > typemin(aelt) && return false
            typemax(celt) < typemax(aelt) && return false
        end
        return true
    end
end
# Note: using `randsubseq` instead of `rand` to avoid repetition.

function inputmattypes(cmat, mattypes = mattypes)
    # Skip if destination type is "too small"
    cmat <: Union{Bidiagonal, Tridiagonal, SymTridiagonal,
                  UnitLowerTriangular, UnitUpperTriangular,
                  Hermitian, Symmetric} && return []
    filter(mattypes) do amat
        cmat <: Diagonal && (amat <: Diagonal || return false)
        cmat <: LowerTriangular && (amat <: LTri || return false)
        cmat <: UpperTriangular && (amat <: UTri || return false)
        return true
    end
end

n_samples = 1.5
# n_samples = Inf  # to try all combinations
for cmat in mattypes,
    amat in sample(inputmattypes(cmat), n_samples),
    bmat in sample(inputmattypes(cmat), n_samples),
    celt in celtypes,
    aelt in sample(inputeltypes(celt), n_samples),
    belt in sample(inputeltypes(celt), n_samples)

    push!(testdata, (cmat{celt}, amat{aelt}, bmat{belt}))
end

@testset "Alternative multiplication operators" begin
    for T in (Int, Float32, Float64, BigFloat)
        a = [T[1, 2], T[-3, 7]]
        b = [T[5, 11], T[-13, 17]]
        @test map(⋅, a, b) == map(dot, a, b) == [27, 158]
        @test map(⊙, a, b) == map(hadamard, a, b) == [a[1].*b[1], a[2].*b[2]]
        @test map(⊗, a, b) == map(tensor, a, b) == [a[1]*transpose(b[1]), a[2]*transpose(b[2])]
        @test hadamard!(fill(typemax(Int), 2), T[1, 2], T[-3, 7]) == [-3, 14]
        @test tensor!(fill(typemax(Int), 2, 2), T[1, 2], T[-3, 7]) == [-3 7; -6 14]
    end

    @test_throws DimensionMismatch [1,2] ⊙ [3]
    @test_throws DimensionMismatch hadamard!([0, 0, 0], [1,2], [-3,7])
    @test_throws DimensionMismatch hadamard!([0, 0], [1,2], [-3])
    @test_throws DimensionMismatch hadamard!([0, 0], [1], [-3,7])
    @test_throws DimensionMismatch tensor!(Matrix{Int}(undef, 2, 2), [1], [-3,7])
    @test_throws DimensionMismatch tensor!(Matrix{Int}(undef, 2, 2), [1,2], [-3])

    u, v = [2+2im, 3+5im], [1-3im, 7+3im]
    @test u ⋅ v == conj(u[1])*v[1] + conj(u[2])*v[2]
    @test u ⊙ v == [u[1]*v[1], u[2]*v[2]]
    @test u ⊗ v == [u[1]*v[1] u[1]*v[2]; u[2]*v[1] u[2]*v[2]]
    @test hadamard(u, v) == u ⊙ v
    @test tensor(u, v)   == u ⊗ v
    dest = similar(u)
    @test hadamard!(dest, u, v) == u ⊙ v
    dest = Matrix{Complex{Int}}(undef, 2, 2)
    @test tensor!(dest, u, v) == u ⊗ v

    for (A, B, b) in (([1 2; 3 4], [5 6; 7 8], [5,6]),
                      ([1+0.8im 2+0.7im; 3+0.6im 4+0.5im],
                       [5+0.4im 6+0.3im; 7+0.2im 8+0.1im],
                       [5+0.6im,6+0.3im]))
        @test A ⊗ b == cat(A*b[1], A*b[2]; dims=3)
        @test A ⊗ B == cat(cat(A*B[1,1], A*B[2,1]; dims=3),
                           cat(A*B[1,2], A*B[2,2]; dims=3); dims=4)
    end

    A, B = reshape(1:27, 3, 3, 3), reshape(1:4, 2, 2)
    @test A ⊗ B == [a*b for a in A, b in B]

    # Adjoint/transpose is a dual vector, not an AbstractMatrix
    v = [1,2]
    @test_throws ErrorException v ⊗ v'
    @test_throws ErrorException v ⊗ transpose(v)
    @test_throws ErrorException v' ⊗ v
    @test_throws ErrorException transpose(v) ⊗ v
    @test_throws ErrorException v' ⊗ v'
    @test_throws ErrorException transpose(v) ⊗ transpose(v)
    @test_throws ErrorException v' ⊗ transpose(v)
    @test_throws ErrorException transpose(v) ⊗ v'
    @test_throws ErrorException A ⊗ v'
    @test_throws ErrorException A ⊗ transpose(v)

    # Docs comparison to `kron`
    v, w = [1,2,3], [5,7]
    @test kron(v,w) == vec(w ⊗ v)
    @test w ⊗ v == reshape(kron(v,w), (length(w), length(v)))
end

@testset "mul!(::$TC, ::$TA, ::$TB, α, β)" for (TC, TA, TB) in testdata
    if needsquare(TA)
        na1 = na2 = rand(sizecandidates)
    else
        na1, na2 = rand(sizecandidates, 2)
    end
    if needsquare(TB)
        nb2 = na2
    elseif needsquare(TC)
        nb2 = na1
    else
        nb2 = rand(sizecandidates)
    end
    asize = (na1, na2)
    bsize = (na2, nb2)
    csize = (na1, nb2)

    @testset for α in Any[true, eltype(TC)(1), _rand(eltype(TC))],
                 β in Any[false, eltype(TC)(0), _rand(eltype(TC))]

        C = _rand(TC, csize)
        A = _rand(TA, asize)
        B = _rand(TB, bsize)

        # This is similar to how `isapprox` choose `rtol` (when
        # `atol=0`) but consider all number types involved:
        rtol = max(rtoldefault.(real.(eltype.((C, A, B))))...,
                   rtoldefault.(real.(typeof.((α, β))))...)

        Cc = copy(C)
        Ac = Matrix(A)
        Bc = Matrix(B)
        returned_mat = mul!(C, A, B, α, β)
        @test returned_mat === C
        @test collect(returned_mat) ≈ α * Ac * Bc + β * Cc  rtol=rtol

        y = C[:, 1]
        x = B[:, 1]
        yc = Vector(y)
        xc = Vector(x)
        returned_vec = mul!(y, A, x, α, β)
        @test returned_vec === y
        @test collect(returned_vec) ≈ α * Ac * xc + β * yc  rtol=rtol

        if TC <: Matrix
            @testset "adjoint and transpose" begin
                @testset for fa in [identity, adjoint, transpose],
                             fb in [identity, adjoint, transpose]
                    fa === fb === identity && continue

                    Af = fa === identity ? A : fa(_rand(TA, reverse(asize)))
                    Bf = fb === identity ? B : fb(_rand(TB, reverse(bsize)))

                    Ac = collect(Af)
                    Bc = collect(Bf)
                    Cc = collect(C)

                    returned_mat = mul!(C, Af, Bf, α, β)
                    @test returned_mat === C
                    @test collect(returned_mat) ≈ α * Ac * Bc + β * Cc  rtol=rtol
                end
            end
        end

        if isnanfillable(C)
            @testset "β = 0 ignores C .= NaN" begin
                parent(C) .= NaN
                Ac = Matrix(A)
                Bc = Matrix(B)
                returned_mat = mul!(C, A, B, α, zero(eltype(C)))
                @test returned_mat === C
                @test collect(returned_mat) ≈ α * Ac * Bc  rtol=rtol
            end
        end

        if isnanfillable(A)
            @testset "α = 0 ignores A .= NaN" begin
                parent(A) .= NaN
                Cc = copy(C)
                returned_mat = mul!(C, A, B, zero(eltype(A)), β)
                @test returned_mat === C
                @test collect(returned_mat) ≈ β * Cc  rtol=rtol
            end
        end
    end
end

end  # module


text:
# Object-based wrapper for BLIS utility routines.
#

macro blis_object_api_lvutil_mkspetial(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           a,      xObj)
    end
end

@blis_object_api_lvutil_mkspetial mkherm
@blis_object_api_lvutil_mkspetial mksymm
@blis_object_api_lvutil_mkspetial mktrim

macro blis_object_api_lvutil_norm(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           x,      xObj,
                           norm,   xObj)
    end
end

@blis_object_api_lvutil_norm norm1v
@blis_object_api_lvutil_norm normfv
@blis_object_api_lvutil_norm normiv
@blis_object_api_lvutil_norm norm1m
@blis_object_api_lvutil_norm normfm
@blis_object_api_lvutil_norm normim

macro blis_object_api_lvutil_rand(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           x,      xObj)
    end
end

@blis_object_api_lvutil_rand randv
@blis_object_api_lvutil_rand randnv
@blis_object_api_lvutil_rand randm
@blis_object_api_lvutil_rand randnm

macro blis_object_api_lvutil_print(funcname)
    return quote
        @blis_ccall_object($funcname,
                           Cvoid,
                           s1,     Ptr{Cchar},
                           x,      xObj,
                           format, Ptr{Cchar},
                           s2,     Ptr{Cchar})
    end
end

@blis_ccall_object(asumv,
                   Cvoid,
                   x,     xObj,
                   asum,  xObj)
@blis_ccall_object(sumsqv,
                   Cvoid,
                   x,     xObj,
                   scale, xObj,
                   sumsq, xObj)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # CBIOMES Global Model Output
#
# Here we retrieve model output from the [CBIOMES](https://cbiomes.org) project via two methods:
#
# - 1. in array format from the [MIT-CBIOMES opendap](http://engaging-opendap.mit.edu:8080/las/) (e.g. visit [this page](http://engaging-opendap.mit.edu:8080/las/UI.vm#panelHeaderHidden=false;differences=false;autoContour=false;xCATID=3C6AA795DF3E9F4E1208CEFE8341F298;xDSID=id-ab2a4e0c65;varid=MXLDEPTH-id-cdfa319965;imageSize=auto;over=xy;compute=Nonetoken;tlo=15-Jan-1992%2000:00;thi=15-Jan-1992%2000:00;catid=3C6AA795DF3E9F4E1208CEFE8341F298;dsid=id-ab2a4e0c65;varid=MXLDEPTH-id-cdfa319965;avarcount=0;xlo=-180;xhi=180;ylo=-90;yhi=90;operation_id=Plot_2D_XY_zoom;view=xy))
# - 2. in tabular format from the [Simons CMAP data base](https://cmap.readthedocs.io/en/latest/) (go to [this page](https://cmap.readthedocs.io/en/latest/catalog/datasets/Darwin_clim.html#darwin-clim))
#

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Using Opendap
#
# The `NCDatasets.jl` package readily supports lazy access to `opendap` data sets. The function below for example access Iron concentration for a chosen depth level and month. Only this two-dimensional slice is transferred over the network.
#
# _warning: this method has failed on mybinder in the past_
#
# ### ➥ `Array` (gridded data)
#
#

# + {"slideshow": {"slide_type": "subslide"}}
using NCDatasets, Plots

srv="http://engaging-opendap.mit.edu:8080/thredds/dodsC/las/"
fil="id-2e0ea5ca2c/data_usr_local_tomcat_content_cbiomes_20200206_17_Nutrients_FeT.0001.nc.jnl"

function test_opendapp(k,t)
    test=Dataset(srv*fil)
    lon=test["LON_C"]
    lat=test["LAT_C"]
    tmp=test["FeT"][:,:,k,t]
    tmp[findall(tmp.<=-0.99e34)].=NaN
    println( test["FeT"].attrib["long_name"]*" in "*test["FeT"].attrib["units"] )
    return vec(lon)[1:2:end],vec(lat)[1:2:end],transpose(tmp[1:2:end,1:2:end])
end

test=Dataset(srv*fil)

# + {"slideshow": {"slide_type": "subslide"}}
heatmap(test_opendapp(1,1))

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Using pycmap
#
# The `space_time` API call below returns the same data as before but as a python table object which is easily saved to a `csv` file.
#
# ### ➥ `DataFrame` (tabular data)

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# _Pre-requisites:_
#
# - 1. _install the [PyCmap](https://github.com/simonscmap/pycmap) python package and its dependencies using `pip`_
# - 2. _compile [PyCall.jl](https://github.com/simonscmap/pycmap) using external python distribution that installed `PyCmap`_
# - 3. _obtain your own API key from the [SimonsCMAP website](https://simonscmap.com) (free; takes `<30s`)_
# - 4. _import `pycmap` via `pycall`_

# + {"slideshow": {"slide_type": "skip"}}
if false
    run(`pip install pycmap`) #pycmap is used via PyCall later
    run(pipeline(`which python`,"whichpython.txt")) #external python path
    ENV["PYTHON"]=readline("whichpython.txt")
    import Pkg; Pkg.build("PyCall")
end

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# _Import `pycmap` Into julia:_
#
# [PyCmap](https://github.com/simonscmap/pycmap) is the `Python` API that we will use in `julia` (to query the `CMAP` data base) via the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package.
#
# _You may need to replace `your-own-API-key` (as outline below) with your own API key from the [SimonsCMAP website](https://simonscmap.com) and uncomment the command below._

# + {"slideshow": {"slide_type": "skip"}}
using PyCall
PyCmap = pyimport("pycmap")
#cmap = PyCmap.API(token="your-own-API-key")
cmap = PyCmap.API()

# + {"slideshow": {"slide_type": "skip"}, "cell_type": "markdown"}
# ### Get Data Catalog
#
# _See [SimonsCMAP website](https://simonscmap.com) for more information. The commented `df.to_csv` call below creates `catalog.csv` and writes the data from `df` to this file._

# + {"slideshow": {"slide_type": "skip"}}
df = cmap.get_catalog();
#df.to_csv("catalog.csv")

# + {"slideshow": {"slide_type": "subslide"}}
tables = ["tblDarwin_Nutrient"] # see catalog.csv  for the complete list of tables and variable names
variables = ["FeT"] # see catalog.csv  for the complete list of tables and variable names
startDate = "1994-01-03"
endDate = "1994-01-03"
lat1, lat2 = 10, 70
lon1, lon2 = -180, -80
depth1, depth2 = 0, 10

df = cmap.space_time(tables[1], variables[1], startDate, endDate,
    lat1, lat2, lon1, lon2, depth1, depth2)

pth="../samples/gradients/"
!isdir("$pth") ? mkdir("$pth") : nothing
df.to_csv("$pth"*"FeT.csv")

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# The data can then easily be read, manipulated, and plotted using Julia's `DataFrames.jl` package.
# -

using CSV, DataFrames
df = CSV.File("$pth"*"FeT.csv") |> DataFrame!

# + {"slideshow": {"slide_type": "subslide"}}
jan=findall( (df.time .== df.time[1]).&( (!ismissing).(df.FeT)) )
jan=rand(jan,1000)
scatter(df.lon[jan],df.lat[jan],zcolor=df.FeT[jan],title="random selection")
# -

max_stars_repo_name: gaelforget/MarineEcosystemNotebooks

max_forks_repo_name: gaelforget/MarineEcosystemNotebooks


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

const DEFAULT_DATASET_PATH = "../../../../src/datasets"

println("--> Dataset 1")
config = "test_dataset1_config.txt"
generations = 5000
st1 = run_NEAT_epochs(config, generations, false)
print_best_solution_perfomance(st1)

plot_neat_fitness(st1, imagename="fitness-dt1.png")
plot_neat_IP(st1, imagename="IP-dt1.png", nbins=10)
plot_output_distribution(st1, imagename="distribution-dt1.png")
plot_complexity(st1, imagename="solution-complexity-dt1.png")
plot_total_accuracy(st1, imagename="total-accuracy-dt1.png")
plot_output_distribution_test(st1, imagename="output-test-distribution-dt1.png")

#=
plot_neat_fitness(st1, imagename="fitness-dt1-pt.png", language=:pt)
plot_neat_IP(st1, imagename="IP-dt1-pt.png", nbins=10, language=:pt)
plot_output_distribution(st1, imagename="distribution-dt1-pt.png", language=:pt)
plot_complexity(st1, imagename="solution-complexity-dt1-pt.png", language=:pt)
plot_total_accuracy(st1, imagename="total-accuracy-dt1-pt.png", language=:pt)
plot_output_distribution_test(st1, imagename="output-test-distribution-dt1-pt.png", language=:pt)
=#

#=
st1 = get_simulation_state("checkpoint-dt1-limited-normalised-5000g-t2.jld2")
plot_neat_fitness(st1, imagename="fitness-dt1-crossover.png")
plot_neat_IP(st1, imagename="IP-dt1-crossover.png", nbins=10)
plot_neat_species(st1, imagename="species-dt1-crossover.png")
plot_output_distribution(st1, imagename="distribution-dt1-crossover.png")
plot_complexity(st1, imagename="solution-complexity-dt1-crossover.png")
plot_accuracy(st1, imagename="accuracy-dt1-crossover.png")
plot_total_accuracy(st1, imagename="total-accuracy-dt1-crossover.png")
plot_output_distribution_test(st1, imagename="output-test-distribution-best-NEAT-crossover.png")
plot_confusion_matrix(st1; imagename="confusion-matrix-testset-dt1-crossover.png", set=:test)
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
"""
    AbstractMoments{N}

A representation of statistical moments of a distribution

The following functions are supported
- `n_moments(m)`: get the number of recorded moments

The following getters return a single moment or
throw an error if the moment has not been recorded
- `mean(m)`: get the mean
- `var(m)`: get the variance
- `skewness(m)`: get the variance
- `kurtosis(m)`: get the variance
- `getindex(m,i)`: get the ith moment, i.e. indexing m[i]

The basic implementation `Moments` is immutable and
`convert(AbstractArray, m::Moments)` returns an `SArray{N,T}`.

# Examples
```jldoctest am; output = false
m = Moments(1,0.2);
n_moments(m) == 2
# output
true
```
```jldoctest am; output = false
var(m) == m[2]
# output
true
```
```julia
kurtosis(m) # throws error because its above 2nd moment
```
"""
abstract type AbstractMoments{N} end
n_moments(::Type{<:AbstractMoments{N}}) where N = N
n_moments(m::AbstractMoments{N}) where N = N
Distributions.mean(m::AbstractMoments) = n_moments(m) >= 1 ? m[1] :
    error("mean not recorded")
Distributions.var(m::AbstractMoments) = n_moments(m) >= 2 ? m[2] :
    error("variance not recorded")
Distributions.std(m::AbstractMoments) = n_moments(m) >= 2 ? sqrt(m[2]) :
    error("std not recorded")
Distributions.skewness(m::AbstractMoments) = n_moments(m) >= 3 ? m[3] :
    error("skewness not recorded")
Distributions.kurtosis(m::AbstractMoments) = n_moments(m) >= 4 ? m[4] :
    error("kurtosis not recorded")

struct Moments{N,T} <: AbstractMoments{N}
    all::SVector{N,T}
end
Moments(x...) = Moments(SVector{length(x)}(promote(x...)))
Moments() = Moments(SA[])
Base.getindex(m::Moments, i) = n_moments(m) >= i ? m.all[i] :
    error("$(i)th moment not recorded.")
Base.convert(::Type{AbstractArray}, m::Moments) = m.all

"""
    moments(D, ::Val{N} = Val(2))

Get the first N moments of a distribution.

See also type [`AbstractMoments`](@ref).

## Examples
```jldoctest; output = false, filter = r"Moments{4,"
moments(LogNormal(), Val(4))  # first four moments
moments(Normal())  # mean and variance

# output
Moments{2,Float64}([0.0, 1.0])
```
"""
function moments(d::Distribution, ::Val{N} = Val(2)) where N
    N > 4 && error("getting moments above 4 not yet implemented for distribution $(typeof(d)).")
    N == 4 && return(Moments(mean(d), var(d), skewness(d), kurtosis(d)))
    N == 3 && return(Moments(mean(d), var(d), skewness(d)))
    N == 2 && return(Moments(mean(d), var(d)))
    N == 1 && return(Moments(mean(d)))
    N == 0 && return(Moments())
    error("must provide a positive N")
end

"""
    fit(D, m)

Fit a statistical distribution of type `D` to given moments `m`.

# Arguments
- `D`: The type of distribution to fit
- `m`: The moments of the distribution

# Notes
This can be used to approximate one distribution by another.

See also [`AbstractMoments`](@ref), [`moments`](@ref).

# Examples
```jldoctest fm1; output = false
D = fit(LogNormal, Moments(3.2,4.6));
(mean(D), var(D)) .≈ (3.2,4.6)
# output
(true, true)
```
```jldoctest fm1; output = false
D = fit(LogNormal, moments(Normal(3,1.2)));
(mean(D), std(D)) .≈ (3,1.2)
# output
(true, true)
```
```julia
plot(D); lines(!Normal(3,1.2))
```
"""
function StatsBase.fit(::Type{<:Distribution}, m::AbstractMoments) end

struct QuantilePoint
    p
    q
    QuantilePoint(p,q) = 0 < p < 1 ? new(p,q) :
        error("p must be in (0,1)")
end
QuantilePoint(qp::QuantilePoint; p = qp.p, q = qp.q) = QuantilePoint(p,q)
Base.show(io::IO, qp::QuantilePoint) = print(io, "(p=$(qp.p),q=$(qp.q))")
function Base.isless(x::QuantilePoint,y::QuantilePoint)
    is_equal_q = (x.q == y.q)
    ((x.p == y.p) && !is_equal_q) && error("incompatible QuantilePoints: $x,$y")
    isless = (x.q < y.q)
    # for different p, q needs to be different
    (isless && (x.p > y.p)) && error("incompatible QuantilePoints: $(x),$(y)")
    (!isless && !is_equal_q && (x.p < y.p))  && error("incompatible QuantilePoints: $x,$y")
    return(isless)
end

macro qp(p,q) :(QuantilePoint($(esc(p)), $(esc(q)))) end
macro qp_ll(q0_025) :(QuantilePoint(0.025, $(esc(q0_025)))) end
macro qp_l(q0_05) :(QuantilePoint(0.05, $(esc(q0_05)))) end
macro qp_m(median) :(QuantilePoint(0.5, $(esc(median)))) end
macro qp_u(q0_95) :(QuantilePoint(0.95, $(esc(q0_95)))) end
macro qp_uu(q0_975) :(QuantilePoint(0.975, $(esc(q0_975)))) end

macro qs_cf90(q0_05,q0_95)
    :(Set([QuantilePoint(0.05,$(esc(q0_05))),QuantilePoint(0.95,$(esc(q0_95)))])) end
macro qs_cf95(q0_025,q0_975)
    :(Set([QuantilePoint(0.025,$(esc(q0_025))),QuantilePoint(0.975,$(esc(q0_975)))])) end

"""
    fit(D, lower, upper)

Fit a statistical distribution to a set of quantiles

# Arguments
- `D`: The type of the distribution to fit
- `lower`:  lower QuantilePoint (p,q)
- `upper`:  upper QuantilePoint (p,q)

# Notes
Several macros help to construct QuantilePoints
- `@qp(p,q)`    quantile at specified p: `QuantilePoint(p, q)`
- `@qp_ll(q0_025)`  quantile at very low p: `QuantilePoint(0.025, q0_025)`
- `@qp_l(q0_05)`    quantile at low p: `QuantilePoint(0.05, q0_05)`
- `@qp_m(median)`   quantile at median: `QuantilePoint(0.5, median)`
- `@qp_u(q0_95)`    quantile at high p: `QuantilePoint(0.95, q0_95)`
- `@qp_uu(q0_975)`  quantile at very high p: `QuantilePoint(0.975, q0_975)`

# Examples
```jldoctest; output = false
d = fit(LogNormal, @qp_m(3), @qp_uu(5));
quantile.(d, [0.5, 0.975]) ≈ [3,5]
# output
true
```
"""
function StatsBase.fit(::Type{<:Distribution}, lower::QuantilePoint, upper::QuantilePoint) end

"""
    fit(D, val, qp, ::Val{stats} = Val(:mean))

Fit a statistical distribution to a quantile and given statistics

# Arguments
- `D`: The type of distribution to fit
- `val`: The value of statistics
- `qp`: QuantilePoint(p,q)
- `stats` Which statistics to fit: defaults to `Val(:mean)`.
   Alternatives are: `Val(:mode)`, `Val(:median)`

# Examples
```jldoctest fm2; output = false
d = fit(LogNormal, 5, @qp_uu(14));
(mean(d),quantile(d, 0.975)) .≈ (5,14)
# output
(true, true)
```
```jldoctest fm2; output = false
d = fit(LogNormal, 5, @qp_uu(14), Val(:mode));
(mode(d),quantile(d, 0.975)) .≈ (5,14)
# output
(true, true)
```
"""
function StatsBase.fit(D::Type{<:Distribution}, val, qp::QuantilePoint, ::Val{stats} = Val(:mean)) where stats
    stats == :mean && return(fit_mean_quantile(D, val, qp))
    stats == :mode && return(fit_mode_quantile(D, val, qp))
    stats == :median && return(fit(D, @qp_m(val), qp))
    error("unknown stats: $stats")
end


text:
function lorentz_likelihood_hermitian_band_save(p,y::Vector{Float64},aex,al_small,indx,logdeta)
# Computes the likelihood of a sum of Exponential/Cosine kernels (which have a
# Lorentzian power spectrum) utilizing the approach of Ambikasaran (2015)
#   Numer. Linear Algebra Appl. 2015; 22:1102-1114 DOI: 10.1002/nla
# Uses (2.30) from Rasmussen & Williams (2006) to compute the log likelihood.
# The form of the kernel is:
#
#   K(t_i,t_j) = w_i \delta_{ij} + \sum_{k=1}^{p} \alpha_k \exp{-\beta_k |t_i-t_j|}
#
# This is represented by the matrix A_{ij} = K(t_i,t_j).
# The log likelihood is given by:
# ln(L) = -1/2 y^T A^{-1} y - 1/2 |det(A)|
# Note: if a \beta_k is complex, then its complex conjugate *must* be included as well.
#
# First, define band-diagonal matrix for kernel:
# [x_i,{r_{i,1},...,r_{i,p}},{l_{i,1},...,l_{i,p}}]

#tic()
n = length(y)
# There are p+1 sub-diagonals, p+1 super-diagonals + diagonal
# for a total of 2*p+3 non-zero diagonals:
nex = (2p+1)*n-2p
bex = zeros(Complex{Float64},nex)
for i=1:n
# Compute actual indices:
  irow =(i-1)*(1+2p)+1
  bex[irow] = y[i]
# Diagonal noise:
  for j=1:p
    if i < n
      bex[irow+  j] = 0.0
      bex[irow+p+j] = 0.0
    end
  end
end

# Specify the number of bands below & above the diagonal:
m1 = p+1
m2 = p+1
# Solve the equation A^{-1} y = b using band-diagonal LU back-substitution on
# the extended equations: A_{ex}^{-1} y_{ex} = b_{ex}:
banbks(aex,nex,m1,m2,al_small,indx,bex)
# Now select solution to compute the log likelihood:
# The equation A^{-1} y = b has been solved in bex (the extended vector).
# So, I need to pick out the b portion from bex, and take the dot product
# with y (which is the residuals of the data minus model, which is correlated
# noise that we are modeling with the multi-Lorentzian covariance function):
log_like = 0.0
for i=1:n
  log_like += real(bex[(i-1)*(2p+1)+1])*y[i]
end
# Convert this to log likelihood:
log_like = -0.5*log_like
return log_like
end


text:
abstract type Separator end

"""
ConstraintSeparator is a separator that represents a constraint defined directly
using `@constraint`.
"""
struct ConstraintSeparator{C, II, ex<:Union{Operation,Expr}} <: Separator
    variables::Vector{Symbol}
    constraint::II  # Interval or IntervalBox
    contractor::C
    expression::ex
end

ConstraintSeparator(constraint, contractor, expression) = ConstraintSeparator(contractor.variables, constraint, contractor, expression)

"""CombinationSeparator is a separator that is a combination (union, intersection,
or complement) of other separators.
"""
struct CombinationSeparator{F, ex<:Union{Operation,Expr}} <: Separator
    variables::Vector{Symbol}
    separator::F
    expression::ex
end

function (S::ConstraintSeparator)(X::IntervalBox)
    C = S.contractor
    a, b = S.constraint.lo, S.constraint.hi

    inner = C(IntervalBox(Interval(a, b)), X)

    local outer

    if a == -∞
        outer = C(IntervalBox(Interval(b, ∞)), X)

    elseif b == ∞
        outer = C(IntervalBox(Interval(-∞, a)), X)

    else

        outer1 = C(IntervalBox(Interval(-∞, a)), X)
        outer2 = C(IntervalBox(Interval(b, ∞)), X)

        outer = outer1 ∪ outer2
    end

    return (inner, outer)
end

"""`parse_comparison` parses comparisons like `x >= 10`
into the corresponding interval, expressed as `x ∈ [10,∞]`

Returns the expression and the constraint interval

TODO: Allow something like [3,4]' for the complement of [3,4]
"""

function parse_comparison(ex::Expr)
    expr, limits =
    @match ex begin
       ((a_ <= b_) | (a_ < b_) | (a_ ≤ b_))   => (a, (-∞, b))
       ((a_ >= b_) | (a_ > b_) | (a_ ≥ b_))   => (a, (b, ∞))
       ((a_ == b_) | (a_ = b_))   => (a, (b, b))

       ((a_ <= b_ <= c_)
        | (a_ < b_ < c_)
        | (a_ <= b_ < c)
        | (a_ < b_ <= c))         => (b, (a, c))

       ((a_ >= b_ >= c_)
       | (a_ > b_ > c_)
       | (a_ >= b_ > c_)
       | (a_ > b_ >= c))          => (b, (c, a))

       ((a_ ∈ [b_, c_])
       | (a_ in [b_, c_])
       | (a_ ∈ b_ .. c_)
       | (a_ in b_ .. c_))        => (a, (b, c))

       _                          => (ex, (-∞, ∞))

   end

   a, b = limits
   return (expr, a..b)   # expr ∈ [a,b]

end

function parse_comparison(ex::Operation)

    if isa(ex.args[1], ModelingToolkit.Constant)
        if ex.op == <
            a = ex.args[1].value
            b = Inf
        elseif ex.op == >
            a = -Inf
            b = ex.args[1].value
        end
        return (ex.args[2], a..b)
    elseif isa(ex.args[2], ModelingToolkit.Constant)
        if ex.op == >
            a = ex.args[2].value
            b = Inf
        elseif ex.op == <
            a = -Inf
            b = ex.args[2].value
        else
            a = ex.args[2].value
            b = ex.args[2].value
        end
        return (ex.args[1], a..b)
    end

end

function new_parse_comparison(ex)
    # @show ex
    if @capture ex begin
            (op_(a_, b_))
        end

        #return (op, a, b)
        # @show op, a, b

    elseif ex.head == :comparison
        println("Comparison")
        symbols = ex.args[1:2:5]
        operators = ex.args[2:2:4]

        # @show symbols
        # @show operators

    end
end

function make_constraint(expr, constraint, var =[])

    if isa(expr, Symbol)
        expr = :(1 * $expr)  # make into an expression!
    end

    contractor_name = make_symbol(:C)

    full_expr = Meta.quot(:($expr ∈ $constraint))

    contractor_code = make_contractor(expr, var)

    code = quote end
    push!(code.args, :($(esc(contractor_name)) = $(contractor_code)))

    push!(code.args, :(ConstraintSeparator($constraint, $(esc(contractor_name)), $full_expr)))

    code
end

make_constraint(expr::Variable, constraint) = make_constraint(Operation(expr), constraint)

function make_constraint(expr::Operation, constraint, var=[])
    C = Contractor(var, expr)
    ex = expr ∈ constraint
    ConstraintSeparator(constraint, C, ex)
end

"""Create a separator from a given constraint expression, written as
standard Julia code.

e.g. `C = @constraint x^2 + y^2 <= 1`

The variables (`x` and `y`, in this case) are automatically inferred.
External constants can be used as e.g. `\$a`:

```
a = 3
C = @constraint x^2 + y^2 <= \$a
```
"""
macro constraint(ex::Expr, variables = [])
    expr, constraint = parse_comparison(ex)
    isa(variables, Array) ? var = [] : var = variables.args
    make_constraint(expr, constraint, var)
end

"""
Create a separator without the use of macros using ModelingToolkit

e.g  vars = @variables x y z
S = Separator(vars, x^2+y^2<1)
X= IntervalBox(-0.5..1.5, -0.5..1.5, -0.5..1.5)
S(X)
"""
function Separator(variables, ex::Operation)
    expr, constraint = parse_comparison(ex)
    make_constraint(expr, constraint, variables)
end

Separator(ex::Operation) = Separator([], ex)

Separator(vars::Union{Vector{Operation}, Tuple{Vararg{Operation,N}}}, g::Function) where N = Separator(vars, g(vars...))

Separator(vars, f::Function) = Separator([Variable(Symbol(i))() for i in vars], f([Variable(Symbol(i))() for i in vars]...)) # if vars is not vector of variables

function show(io::IO, S::Separator)
    println(io, "Separator:")
    print(io, "  - variables: ")
    print(io, join(map(string, S.variables), ", "))
    println(io)
    print(io, "  - expression: ")
    print(io, S.expression)
end

(S::CombinationSeparator)(X) = S.separator(X)

"Unify the variables of two separators"
function unify_variables(vars1, vars2)

    variables = unique(sort(vcat(vars1, vars2)))
    numvars = length(variables)  # total number of variables

    indices1 = indexin(vars1, variables)
    indices2 = indexin(vars2, variables)

    where1 = zeros(Int, numvars)  # inverse of indices1
    where2 = zeros(Int, numvars)

    for (i, which) in enumerate(indices1)
        where1[which] = i
    end

    for (i, which) in enumerate(indices2)
        where2[which] = i
    end

    return variables, indices1, indices2, where1, where2
end

# TODO: when S1 and S2 have different variables -- amalgamate!

"""
    ∩(S1::Separator, S2::Separator)

Separator for the intersection of two sets given by the separators `S1` and `S2`.
Takes an iterator of intervals (`IntervalBox`, tuple, array, etc.), of length
equal to the total number of variables in `S1` and `S2`;
it returns inner and outer tuples of the same length
"""
function ∩(S1::Separator, S2::Separator)

    #=variables, indices1, indices2, where1, where2 = unify_variables(S1.variables, S2.variables)
    numvars = length(variables)=#
    variables = S1.variables   # as S1 and S2 have same variables
    f = X -> begin

       inner1, outer1 = S1(X)
       inner2, outer2 = S2(X)
       inner = inner1 ∩ inner2
       outer = outer1 ∪ outer2

        #=
        inner1, outer1 = S1(IntervalBox([X[i] for i in indices1]...))
        inner2, outer2 = S2(IntervalBox([X[i] for i in indices2]...))

        if any(isempty, inner1)
            inner1 = emptyinterval(IntervalBox(X))
        else
            inner1 = [i ∈ indices1 ? inner1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer1)
            outer1 = emptyinterval(IntervalBox(X))
        else
            outer1 = [i ∈ indices1 ? outer1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, inner2)
            inner2 = emptyinterval(IntervalBox(X))
        else
            inner2 = [i ∈ indices2 ? inner2[where2[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer2)
            outer2 = emptyinterval(IntervalBox(X))
        else
            outer2 = [i ∈ indices2 ? outer2[where2[i]] : X[i] for i in 1:numvars]
        end

        # Treat as if had X[i] in the other directions, except if empty

        inner = IntervalBox( [x ∩ y for (x,y) in zip(inner1, inner2) ]... )
        outer = IntervalBox( [x ∪ y for (x,y) in zip(outer1, outer2) ]... )
        =#

        return (inner, outer)

    end

    expression = :($(S1.expression) ∩ $(S2.expression))

    return CombinationSeparator(variables, f, expression)

end

function ∪(S1::Separator, S2::Separator)

    #=variables, indices1, indices2, where1, where2 = unify_variables(S1.variables, S2.variables)
    numvars = length(variables)=#
    variables = S1.variables    # S1 and S2 have same variables
    f = X -> begin

        inner1, outer1 = S1(X)
        inner2, outer2 = S2(X)
        inner = inner1 ∪ inner2
        outer = outer1 ∩ outer2

        #=
        inner1, outer1 = S1(IntervalBox([X[i] for i in indices1]...))
        inner2, outer2 = S2(IntervalBox([X[i] for i in indices2]...))

        if any(isempty, inner1)
            inner1 = emptyinterval(IntervalBox(X))
        else
            inner1 = [i ∈ indices1 ? inner1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer1)
            outer1 = emptyinterval(IntervalBox(X))
        else
            outer1 = [i ∈ indices1 ? outer1[where1[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, inner2)
            inner2 = emptyinterval(IntervalBox(X))
        else
            inner2 = [i ∈ indices2 ? inner2[where2[i]] : X[i] for i in 1:numvars]
        end

        if any(isempty, outer2)
            outer2 = emptyinterval(IntervalBox(X))
        else
            outer2 = [i ∈ indices2 ? outer2[where2[i]] : X[i] for i in 1:numvars]
        end

        inner = IntervalBox( [x ∪ y for (x,y) in zip(inner1, inner2) ]... )
        outer = IntervalBox( [x ∩ y for (x,y) in zip(outer1, outer2) ]... )
        =#

        return (inner, outer)
    end

    expression = :($(S1.expression) ∪ $(S2.expression))

    return CombinationSeparator(variables, f, expression)

end

function !(S::Separator)
    f = X -> begin
        inner, outer = S(X)
        return (outer, inner)
    end

    expression = :(!($(S.expression)))

    return CombinationSeparator(S.variables, f, expression)
end


text:
function _get_shifted_jacobian(M, o::Orbit; transform = x -> x, spline=true, tol=0.1,verbose=false)
    tm = o.coordinate.t
    tp = range(0.0,1.0,length=length(o.path))
    if spline
        # Create spline of orbit path
        A = hcat(o.path.energy,o.path.pitch,o.path.r,o.path.z)
        itp = extrapolate(scale(interpolate(A,(BSpline(Cubic(Periodic(OnGrid()))), NoInterp())), tp, 1:4), (Periodic(), Throw()))

        # Create new orbit path
        energy = itp.(tp .- tm, 1)
        pitch = itp.(tp .- tm, 2)
        r = itp.(tp .- tm, 3)
        z = itp.(tp .- tm, 4)

        dt = fill(step(tp)*o.tau_p,length(o.path))
        dt[end] = 0.0
        opath = OrbitPath(o.path.vacuum, o.path.drift, energy,pitch,r,z,r*0.0,dt)
    else
        os = get_orbit(M,o.coordinate,drift=o.path.drift,vacuum=o.path.vacuum)
        opath = os.path
    end

    # Calculate jacobian
    J = _get_jacobian(M, o.coordinate, opath, o.tau_p, transform, tol, verbose)

    # Shift jacobian
    tj = range(0.0,1.0,length=length(J))
    Jitp = extrapolate(scale(interpolate(J,BSpline(Linear())),tj), Periodic())
    Jshifted = max.(Jitp.(tp .+ tm),0.0)

    return Jshifted
end

function get_jacobian(M::AbstractEquilibrium, c::EPRCoordinate;
                      transform = x -> x, spline=true, tol=0.1,
                      verbose=false, kwargs...)
    o = get_orbit(M,c; kwargs...)
    if o.class == :invalid
        return _get_shifted_jacobian(M, o; transform=transform, spline=spline,tol=tol,verbose=verbose)
    end
    return _get_jacobian(M, o.coordinate, o.path, o.tau_p, transform, tol,verbose)
end

function get_jacobian(M::AbstractEquilibrium, o::Orbit;
                      transform = x -> x, spline=true, tol=0.1,
                      verbose=false)
    length(o.path) == 0 && return zeros(length(o.path))
    r0 = o.path.r[1]
    if r0 < o.coordinate.r && !isapprox(r0,o.coordinate.r,rtol=1e-4)
        return _get_shifted_jacobian(M, o; transform=transform, spline=spline,tol=tol,verbose=verbose)
    end
    return _get_jacobian(M, o.coordinate, o.path, o.tau_p, transform, tol,verbose)
end

function check_jacobian(J::Array; tol = 0.1)
    !(any(abs.(J .- circshift(J,-1))./median(J) .> tol))
end

function _get_jacobian(M::AbstractEquilibrium, c::Union{GCParticle,AbstractOrbitCoordinate},
                       o::OrbitPath, tau_p, transform, tol, verbose)
    np = length(o)
    detJ = zeros(np)
    if np == 0 || tau_p == 0.0
        return detJ
    end

    ## Do first half of orbit
    ed = Dual(o.energy[1],(1.0,0.0,0.0,0.0))
    pd = Dual(o.pitch[1], (0.0,1.0,0.0,0.0))
    rd = Dual(o.r[1],     (0.0,0.0,1.0,0.0))
    td = Dual(1e-30,      (0.0,0.0,0.0,1.0))
    zd = one(td)*o.z[1]

    gcp = GCParticle(ed,pd,rd,zd,c.m,c.q)
    path, stat = integrate(M, gcp; tmax=td*tau_p, interp_dt=0.0,
                           one_transit=false, store_path=true,
                           classify_orbit=false, drift=o.drift, vacuum=o.vacuum)

    x = transform([path.energy[end],path.pitch[end],path.r[end],path.z[end]])
    detJ[1] = max(abs(det(hcat((partials(xx) for xx in x)...))),0.0)

    ep = partials(path.energy[end])
    pp = partials(path.pitch[end])
    rp = partials(path.r[end])
    zp = partials(path.z[end])
    np1 = floor(Int,np/2)
    for i=2:np1
        ed = Dual(o.energy[i-1],ep)
        pd = Dual(o.pitch[i-1], pp)
        rd = Dual(o.r[i-1],     rp)
        zd = Dual(o.z[i-1],     zp)
        gcp = GCParticle(ed,pd,rd,zd,c.m,c.q)

        dt = o.dt[i-1]
        dt = dt + 1e-30*iszero(dt)
        path, stat = integrate(M, gcp; tmax=dt, interp_dt=0.0,
                               one_transit=false, store_path=true,
                               classify_orbit=false, drift=o.drift,vacuum=o.vacuum)

        x = transform([path.energy[end],path.pitch[end],path.r[end],path.z[end]])
        ep = partials(path.energy[end])
        pp = partials(path.pitch[end])
        rp = partials(path.r[end])
        zp = partials(path.z[end])
        detJ[i] = max(abs(det(hcat((partials(xx) for xx in x)...))),0.0)
    end

    ## Do second half of orbit
    ed = Dual(o.energy[1],(1.0,0.0,0.0,0.0))
    pd = Dual(o.pitch[1], (0.0,1.0,0.0,0.0))
    rd = Dual(o.r[1],     (0.0,0.0,1.0,0.0))
    td = Dual(-1e-30,     (0.0,0.0,0.0,1.0))
    zd = one(td)*o.z[1]
    gcp = GCParticle(ed,pd,rd,zd,c.m,c.q)
    path, stat = integrate(M, gcp; tmax=td*tau_p, interp_dt=0.0,
                           one_transit=false, store_path=true,
                           classify_orbit=false, drift=o.drift,vacuum=o.vacuum)

    x = transform([path.energy[end],path.pitch[end],path.r[end],path.z[end]])
    detJ[np] = max(abs(det(hcat((partials(xx) for xx in x)...))),0.0)

    ep = partials(path.energy[end])
    pp = partials(path.pitch[end])
    rp = partials(path.r[end])
    zp = partials(path.z[end])
    for i=(np-1):-1:(np1+1)
        ed = Dual(o.energy[i+1],ep)
        pd = Dual(o.pitch[i+1], pp)
        rd = Dual(o.r[i+1],     rp)
        zd = Dual(o.z[i+1],     zp)
        gcp = GCParticle(ed,pd,rd,zd,c.m,c.q)

        dt = o.dt[i]
        dt = -(dt + 1e-30*iszero(dt))
        path, stat = integrate(M, gcp; tmax=dt, interp_dt=0.0,
                               one_transit=false, store_path=true,
                               classify_orbit=false, drift=o.drift,vacuum=o.vacuum)

        x = transform([path.energy[end],path.pitch[end],path.r[end],path.z[end]])
        detJ[i] = max(abs(det(hcat((partials(xx) for xx in x)...))),0.0)

        ep = partials(path.energy[end])
        pp = partials(path.pitch[end])
        rp = partials(path.r[end])
        zp = partials(path.z[end])
    end

    if tol > 0.0 && !check_jacobian(detJ; tol=tol)
        verbose && @warn "Jacobian calculation failed. Setting to zero."
        detJ = zero(detJ)
    end

    return detJ
end

function transform(M::AbstractEquilibrium, o::Orbit, ::Type{EPRParticle}; kwargs...)
    J = _get_jacobian(M, o.coordinate, o.path, o.tau_p; kwargs...)[1]
    return EPRParticle(o.coordinate), J
end

function transform(M::AbstractEquilibrium, gcp::GCParticle, ::Type{EPRParticle}; kwargs...)
    o = get_orbit(M, gcp; kwargs...)
    return transform(M, o; kwargs...)
end

function transform(M::AbstractEquilibrium, c::EPRCoordinate, ::Type{HamiltonianCoordinate})
    KE_d = Dual(c.energy, (1.0, 0.0, 0.0))
    p_d =  Dual(c.pitch,  (0.0, 1.0, 0.0))
    r_d =  Dual(c.r,      (0.0, 0.0, 1.0))
    c_d = EPRCoordinate(KE_d, p_d, r_d, one(r_d)*c.z, one(r_d)*c.t, c.m, c.q)

    hc_d = HamiltonianCoordinate(M, c_d)

    hc = HamiltonianCoordinate(value(hc_d.energy), value(hc_d.mu), value(hc_d.p_phi), c.m, c.q)
    J = abs(det(hcat(partials(hc_d.energy), partials(hc_d.mu), partials(hc_d.p_phi))))

    return hc, J
end


text:
function Base.show(io::IO, ::MIME"text/plain", G::GSDist)
    α, g, k, γ = params(G)
    println(io, "Generalized S-Distribution ($(G.dist))")
    println(io, " α: $α")
    println(io, " g: $g")
    println(io, " k: $k")
    print(io, " γ: $γ")
end


text:
#=
 = This module returns a dictionary of symbols and alphanumeric values, as given by a single
 = ASCII line from an International Maritime Meteorological Archive file (IMMA Release 2.5;
 = following http://rda.ucar.edu/datasets/ds540.0/docs/R2.5-imma_short.pdf).  This format
 = consists of a core entry and up to six attachments, some of which are still in the process
 = of being standardized.  Note that the dictionary returns only defined values (that is, the
 = original string if a scale parameter is zero, a scaled float if the scale is positive, and
 = a float converted from a different base if the scale is negative) and for the time being,
 = attachments 5 and 6 are not parsed (at least until 5 is formalized) - RD March 2016.
 =#

module ICOADSDict                                       # first define the values contained in the core
export imma                                             # and attachments by symbol, location, and scale

const symb0 = [:YR :MO :DY  :HR :LAT :LON :IM :ATTC :TI :LI :DS :VS :NID :II :ID :C1 :DI :D :WI  :W :VI :VV :WW :W1 :SLP :A :PPP :IT :AT :WBTI :WBT :DPTI :DPT :SI :SST :N :NH :CL :HI :H :CM :CH :WD :WP :WH :SD :SP :SH]'
const leng0 = [  4   2   2    4    5    6   2     1   1   1   1   1    2   2   9   2   1  3   1   3   1   2   2   1    5  1    3   1   4     1    4     1    4   2    4  1   1   1   1  1   1   1   2   2   2   2   2   2]'
const scal0 = [  1   1   1 0.01 0.01 0.01   1     1   1   1   1   1    1   1   0   0   1  1   1 0.1   1   1   1   1  0.1  1  0.1   1 0.1     1  0.1     1  0.1   1  0.1  1   1   0   1  0   0   0   1   1   1   1   1   1]'

const symb1 = [:ATTI1 :ATTL1 :BSI :B10 :B1 :DCK :SID :PT :DUPS :DUPC :TC :PB :WX :SX :C2 :SQZ :SQA :AQZ :AQA :UQZ :UQA :VQZ :VQA :PQZ :PQA :DQZ :DQA :ND :SF :AF :UF :VF :PF :RF :ZNC :WNC :BNC :XNC :YNC :PNC :ANC :GNC :DNC :SNC :CNC :ENC :FNC :TNC :QCE :LZ :QCZ]'
const leng1 = [     2      2    1    3   2    3    3   2     2     1   1   1   1   1   2    1    1    1    1    1    1    1    1    1    1    1    1   1   1   1   1   1   1   1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    2   1    2]'
const scal1 = [     1      1    0    1   1    1    1   1     1     1   1   1   1   1   1  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36   1 -36 -36 -36 -36 -36 -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36  -36    1   1    1]'

const symb2 = [:ATTI2 :ATTL2 :OS :OP :FM :IX :W2 :SGN :SGT :SGH :WMI :SD2 :SP2 :SH2 :IS :ES :RS :IC1 :IC2 :IC3 :IC4 :IC5 :IR :RRR :TR :QCI :QI1 :QI2 :QI3 :QI4 :QI5 :QI6 :QI7 :QI8 :QI9 :QI10 :QI11 :QI12 :QI13 :QI14 :QI15 :QI16 :QI17 :QI18 :QI19 :QI20 :QI21 :HDG :COG :SOG :SLL :SLHH :RWD :RWS]'
const leng2 = [     2      2   1   1   2   1   1    1    1    2    1    2    2    2   1   2   1    1    1    1    1    1   1    3   1    1    1    1    1    1    1    1    1    1    1     1     1     1     1     1     1     1     1     1     1     1     1    3    3    2    2     3    3    3]'
const scal2 = [     1      1   1   1   1   1   1    1    0    1    1    1    1    1   1   1   1    0    0    0    0    0   1    1   1    1    1    1    1    1    1    1    1    1    1     1     1     1     1     1     1     1     1     1     1     1     1    1    1    1    1     1    1    1]'

const symb3 = [:ATTI3 :ATTL3 :CCCC :BUID :BMP :BSWU :SWU :BSWV :SWV :BSAT :BSRH :SRH :SIX :BSST :MST :MSH :BY :BM :BD :BH :BFL]'
const leng3 = [     2      2     4     6    5     4    4     4    4     4     3    3    1     4    1    3   4   2   2   2    2]'
const scal3 = [     1      1     0     0  0.1   0.1  0.1   0.1  0.1   0.1     1    1    1   0.1    1    1   1   1   1   1    1]'

const symb4 = [:ATTI4 :ATTL4 :C1M :OPM :KOV :COR :TOB :TOT :EOT :LOT :TOH :EOH :SIM :LOV :DOS :HOP :HOT :HOB :HOA :SMF :SME :SMV]'
const leng4 = [     2      2    2    2    2    2    3    3    2    2    1    2    3    3    2    3    3    3    3    5    5    2]'
const scal4 = [     1      1    0    1    0    0    0    0    0    0    0    0    0    1    1    1    1    1    1    1    1    1]'

const symb5 = [:ATTI5 :ATTL5 :WFI :WF :XWI :XW :XDI :XD :SLPI :TAI :TA :XNI :XN :OTHERSTUFF5]'
const leng5 = [     2      2    1   2    1   3    1   2     1    1   4    1   2           71]'
const scal5 = [     1      1    0   1    0   1    0   0     0    0   1    0   0            0]'

const symb6 = [:ATTI6 :ATTL6 :ATTE :OTHERSTUFF6]'
const leng6 = [     2      2     1          999]'
const scal6 = [     1      1     0            0]'

function imma(line::AbstractString)
  vals = Dict{Symbol, Any}()                                                  # initialize a "vals" dictionary
  linlen = length(line)

  subsec = 0                                                                  # then loop through the core and
  substa = 1                                                                  # all available attachments until
  while substa < linlen                                                       # we're past the end of the line
    if substa > 1                                                             # (sta/end are section start/end)
      subsec = parse(Int, line[substa:substa+1])
    end

    symb = eval(symbol("symb", subsec))                                       # use the relevant section arrays
    leng = eval(symbol("leng", subsec))
    scal = eval(symbol("scal", subsec))
    sumlen = sum(leng)
    subend = substa + sumlen - 1

    if subsec < 5                                                             # ignore historical/supplemental
      if subend <= linlen                                                     # isolate this section and pad if
        sublin =      line[substa:subend]                                     # at the end of a truncated line
      else
        sublin = rpad(line[substa:   end], sumlen, ' ')
      end
#     println(subsec, " :", sublin, ": ")

      index = [0; cumsum(leng)] + 1                                           # add any valid entries from each
      for (a, sym) in enumerate(symb)                                         # section and return the dictionary
        tmpval = sublin[index[a]:index[a]+leng[a]-1]
        if tmpval != " "^leng[a]
          if scal[a] > 0
            vals[sym] = float(tmpval) * scal[a]
          elseif scal[a] == 0
            vals[sym] = strip(tmpval)
          else
            vals[sym] = float(parse(Int, tmpval, 36))
          end
#         println(sym, " ", vals[sym])
        end
      end
    end
    substa = subend + 1
  end
  return vals
end

end                                                                           # end ICOADSDict module


text:
function documentedmodules(md = Main)
    modules = Module[]
    for obj in names(md)
        if isa(getfield(md, obj), Module) && isdefined(md.(obj), METADATA)
            push!(modules, md.(obj))
        end
    end
    modules
end

@doc """
Full text search of all docstrings curently available.

**Example:**

```julia
query("query")
```
""" {
    :tags       => ["search", "fulltext"],
    :section    => "Queries and Searches",
    :parameters => [
        (:str, "text to search for in all available documentation")
        ]
    } ->
function query(str::String)
    results = (Any,Entry)[]
    for md in documentedmodules()
        for (obj, ent) in getfield(md, METADATA).entries
            searchmarkdown(ent.docs, str) && push!(results, (obj, ent))
        end
    end
    results
end

function searchmarkdown(md::Markdown.Content, str::String)
    found = false
    for block in md.content
        found = found ||
            if :text in (ns = names(block);)
                contains(block.text, str)
            elseif :code in ns
                contains(block.code, str)
            elseif :content in ns
                searchmarkdown(block, str)
            end
    end
    found
end

@doc """
Search the module `m` for documentation on object `q`.

If `q` is not provided search for documentation related to the module
`m` itself.

**Examples:**

```julia
query(Docile)                     # docs for the module `Docile`
query(query)                      # function and method docs for `query`
query(query, Docile; all = false) # just the function docs in module `Docile`.
```
""" {
    :section => "Queries and Searches",
    :parameters => [
        (:m, "the target module to search through"),
        (:q, "optional object to search the module for"),
        (:all, "optional switch to show all methods associated with a function")
        ]
    } ->
function query(q, ms = documentedmodules(); all = true)
    results = (Any,Entry)[]
    for m in [ms]
        if haskey(getfield(m, METADATA).entries, q)
            push!(results, (q, getfield(m, METADATA).entries[q]))
        end
        if isa(q, Function) && all
            for mt in q.env
                if haskey(getfield(m, METADATA).entries, mt)
                    push!(results, (mt, getfield(m, METADATA).entries[mt]))
                end
            end
        end
    end
    results
end

@doc "Search for documentation in Docile-generated docstrings." -> query

@doc """
Search for documentation related to the method that would be called if
the expression `q` was evaluated. Behavour is similar to that of
`Base.@which`, but returns the documentation as well.

**Examples:**

```julia
@query doctest(Docile)
@query @query
```
""" {
    :section    => "Queries and Searches",
    :parameters => [
        (:q, "the expression to query")
        ]
    } ->
macro query(q)
    if isa(q, Symbol)
        Expr(:call, query, Expr(:quote, symbol(q)))
    elseif isa(q, Expr)
        if q.head == :macrocall
            Expr(:call, query, Expr(:quote, symbol("$(q.args[1])")))
        else
            mt = Expr(:macrocall, symbol("@which"), q)
            esc(:(query($(mt), $(mt).func.code.module, )))
        end
    end
end


text:
#######################################
#### Canopy-Atmosphere Decoupling  ####
#######################################

#' Canopy-Atmosphere Decoupling Coefficient
#'
#' The canopy-atmosphere decoupling coefficient 'Omega'.
#'
#' - data        DataFrame or matrix containing all required input variables
#' - Tair        Air temperature (deg C)
#' - pressure    Atmospheric pressure (kPa)
#' - Ga          Aerodynamic conductance to heat/water vapor (m s-1)
#' - Gs          Surface conductance (m s-1)
#' - approach    Approach used to calculate omega. Either `"Jarvis&McNaughton_1986"` (default)
#'                    or `"Martin_1989"`.
#' - LAI         Leaf area index (m2 m-2), only used if `approach = "Martin_1989"`.
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants   Kelvin - conversion degree Celsius to Kelvin
#'                    cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                    eps - ratio of the molecular weight of water vapor to dry air (-)
#'                    sigma - Stefan-Boltzmann constant (W m-2 K-4)
#'                    Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Details
 The decoupling coefficient Omega ranges from 0 to 1 and quantifies the
#'          linkage of the conditions (foremost humidity and temperature) at the canopy surface
#'          to the ambient air. Values close to 0 indicate well coupled conditions
#'          characterized by high physiological (i.e. stomatal) control on transpiration
#'          and similar conditions at the canopy surface compared to the atmosphere above
#'          the canopy. Values close to 1 indicate the opposite, i.e. decoupled conditions and
#'          a low stomatal control on transpiration (Jarvis & McNaughton 1986).
#'          The `"Jarvis&McNaughton_1986"` approach (default option) is the original
#'          formulation for the decoupling coefficient, given by (for an amphistomatous
#'          canopy):
#'
#'          ``\\Omega = \frac{\epsilon + 1``{\epsilon + 1 + \frac{Ga}{Gc}}}{%
#'          \\Omega = (\epsilon + 1) / ( \epsilon + 1 + Ga/Gc)}
#'
#'          where ``\epsilon = \frac{s``{\\gamma}}{\epsilon = s/\\gamma} is a dimensionless coefficient
#'          with s being the slope of the saturation vapor pressure curve (Pa K-1), and ``\\gamma`` the
#'          psychrometric constant (Pa K-1).
#'
#'          The approach `"Martin_1989"` by Martin 1989 additionally takes radiative coupling
#'          into account:
#'
#'          ``\\Omega = \frac{\epsilon + 1 + \frac{Gr``{Ga}}{\epsilon + (1 + \frac{Ga}{Gs}) (1 + \frac{Gr}{Ga})}}{%
#'          \\Omega = (\epsilon + 1 + Gr/Ga) / (\epsilon + (1 + Ga/Gs) (1 + Gr/Ga))}
#'
#' # Value
 - ``\\Omega`` -: the decoupling coefficient Omega (-)
#'
#' #References
#' Jarvis P_G., McNaughton K_G., 1986: Stomatal control of transpiration:
#'             scaling up from leaf to region. Advances in Ecological Research 15, 1-49.
#'
#'             Martin P., 1989: The significance of radiative coupling between
#'             vegetation and the atmosphere. Agricultural and Forest Meteorology 49, 45-53.
#'
#' #See also
#' [`aerodynamic_conductance!`](@ref), [`surface_conductance`](@ref),
#'          [`equilibrium_imposed_ET`](@ref)
#'
#' ```@example; output = false
#' ```
#' # Omega calculated following Jarvis & McNaughton 1986
#' set_seed(3)
#' df = DataFrame(Tair=rnorm(20,25,1),pressure=100,Ga_h=rnorm(20,0.06,0.01),
#'                  Gs_ms=rnorm(20,0.005,0.001))
#' decoupling(df,approach="Jarvis&McNaughton_1986")
#'
#' # Omega calculated following Martin 1989 (requires LAI)
#' decoupling(df,approach="Martin_1989",LAI=4)
#'
"""
"""
function decoupling(data,Tair="Tair",pressure="pressure",Ga="Ga_h",Gs="Gs_ms",
                       approach=c("Jarvis&McNaughton_1986","Martin_1989"),
                       LAI,Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                       constants=BigleafConstants())

  approach    = match_arg(approach)

  check_input(data,list(Tair,pressure,Ga,Gs))

  Delta   = Esat_slope(Tair,Esat_formula,constants)[,"Delta"]
  gamma   = psychrometric_constant(Tair,pressure,constants)
  epsilon = Delta/gamma

  if (approach == "Jarvis&McNaughton_1986")

    Omega = (epsilon + 1) / (epsilon + 1 + Ga/Gs)

elseif (approach == "Martin_1989")

    if (isnothing(LAI))

      stop("LAI is not provided!")

else

      Gr    = longwave_conductance(Tair,LAI,constants)
      Omega = (epsilon + 1 + Gr/Ga) / (epsilon + 1 + Ga/Gs + Gr/Gs + Gr/Ga)

end
end

  return(Omega)

end

#' Longwave Radiative Transfer Conductance of the Canopy
#'
#' - Tair      Air temperature (deg C)
#' - LAI       Leaf area index (m2 m-2)
#' - constants Kelvin - conversion degree Celsius to Kelvin
#'                  sigma - Stefan-Boltzmann constant (W m-2 K-4)
#'                  cp - specific heat of air for constant pressure (J K-1 kg-1)
#'
#' # Details
 the following Esat_formula is used (Martin, 1989):
#'
#'          ``Gr = 4 \sigma Tair^3 LAI / cp``
#'
#' # Value
 - Gr -: longwave radiative transfer conductance of the canopy (m s-1)
#'
#' #References
#' Martin P., 1989: The significance of radiative coupling between
#'             vegetation and the atmosphere. Agricultural and Forest Meteorology 49, 45-53.
#'
#' ```@example; output = false
#' ```
#' longwave_conductance(25,seq(1,8,1))
#'
#' @export
function longwave_conductance(Tair,LAI,constants=BigleafConstants())

  Tair = Tair + constants.Kelvin

  Gr = 4 * constants.sigma] * Tair^3 * LAI / constants[:cp

  return(Gr)
end


text:
module Sequencing

export
    # Molecules
    Molecules,
    amplify,
    fragment,
    select,
    subsample,
    tag,
    flip,
    # Reads
    paired_reads,
    unpaired_reads,
    nreads,
    ClearSubstitutions,
    FixedProbSubstitutions,
    edit_substitutions,
    make_substitutions,
    generate,
    # Utils
    needed_sample_size,
    expected_coverage,
    sequence,
    X, x, bp,
    makereads,
    coverage_report,
    uncovered_positions,
    uncovered_regions

using BioSequences, FASTX, Random

struct ExpectedCoverage
    val::Int
end

const X = ExpectedCoverage
const x = ExpectedCoverage

Base.:*(val::Int, ::Type{ExpectedCoverage}) = ExpectedCoverage(val)
Base.show(io::IO, v::ExpectedCoverage) = println(io, v.val, "X sequencing coverage")

struct SequenceLength
    val::Int
end

const bp = SequenceLength

Base.:*(val::Int, ::Type{SequenceLength}) = SequenceLength(val)
Base.show(io::IO, x::SequenceLength) = println(io, x.val, " base pairs")

Base.:+(x::SequenceLength, y::SequenceLength) = SequenceLength(x.val + y.val)

Base.:*(cov::ExpectedCoverage, len::SequenceLength) = SequenceLength(cov.val * len.val)

function Base.:*(val::Int, rl::SequenceLength)
    @assert val == 2
    return (rl, rl)
end

Base.div(x::SequenceLength, y::SequenceLength) = div(x.val, y.val)
Base.iterate(x::SequenceLength) = (x, nothing)
Base.iterate(x::SequenceLength, ::Any) = nothing

function needed_sample_size(coverage::ExpectedCoverage, genome_len::SequenceLength, lens)
    return div(coverage * genome_len, sum(lens))
end

expected_coverage(G::Int, L::Int, N::Int) = div(L * N, G)

include("sequencing_views.jl")
include("Molecules.jl")
include("Reads.jl")
include("DSL.jl")

struct CoverageReport
    genome::Vector{LongSequence{DNAAlphabet{2}}}
    covs::Vector{Vector{UInt}}
end

function coverage_report(x)
    covs = [zeros(UInt, length(g)) for g in genome(x)]
    vs = views(x)
    for v in vs
        c = covs[seqid(v)]
        @inbounds for i in min(first(v), last(v)):max(first(v), last(v))
            c[i] = c[i] + 1
        end
    end
    return CoverageReport(genome(x), covs)
end

function uncovered_positions(cr::CoverageReport, thresh = 0)
    return [[i for i in eachindex(c) if c[i] < thresh] for c in cr.covs]
end

function uncovered_regions(cr::CoverageReport, thresh = 0)
    upos = uncovered_positions(cr, thresh)
    out = [Vector{UnitRange}() for i in 1:length(upos)]
    for i in eachindex(upos)
        v = upos[i]
        o = out[i]
        if !isempty(v)
            start = stop = first(v)
            for j in 2:lastindex(v)
                vj = v[j]
                if vj == stop + 1
                    stop = vj
                else
                    push!(o, start:stop)
                    start = stop = vj
                end
            end
            push!(o, start:stop)
        end
    end
    return out
end

function _summarize_cov(V::Vector{UInt64})
    min = typemax(UInt64)
    max = typemin(UInt64)
    sum = typemin(UInt64)
    for v in V
        min = ifelse(v < min, v, min)
        max = ifelse(v > max, v, max)
        sum = sum + v
    end
    return (min, max, sum)
end

function summarize(cr::Sequencing.CoverageReport, bychrom = false)
    stats = [_summarize_cov(c) for c in cr.covs]
    if bychrom
        println("Coverage Summary:")
        for i in eachindex(stats)
            min, max, total = stats[i]
            mean = total / length(cr.covs[i])
            println("Chromosome $i:")
            println("\tmin: $min")
            println("\tmean: $mean")
            println("\tmax: $max")
        end
    else
        min = minimum(x[1] for x in stats)
        total = sum(x[3] for x in stats)
        max = maximum(x[2] for x in stats)
        mean = total / sum(length(x) for x in cr.covs)
        println("Coverage Summary:")
        println("\tmin: $min")
        println("\tmean: $mean")
        println("\tmax: $max")
    end
end

end # module


text:
# ---
# title: 1668. Maximum Repeating Substring
# id: problem1668
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/maximum-repeating-substring/description/>
# hidden: true
# ---
#
# For a string `sequence`, a string `word` is **`k`-repeating** if `word`
# concatenated `k` times is a substring of `sequence`. The `word`'s
# **maximum`k`-repeating value** is the highest value `k` where `word` is
# `k`-repeating in `sequence`. If `word` is not a substring of `sequence`,
# `word`'s maximum `k`-repeating value is `0`.
#
# Given strings `sequence` and `word`, return _the **maximum`k`-repeating
# value** of `word` in `sequence`_.
#
#
#
# **Example 1:**
#
#
#
#     Input: sequence = "ababc", word = "ab"
#     Output: 2
#     Explanation: "abab" is a substring in " _abab_ c".
#
#
# **Example 2:**
#
#
#
#     Input: sequence = "ababc", word = "ba"
#     Output: 1
#     Explanation: "ba" is a substring in "a _ba_ bc". "baba" is not a substring in "ababc".
#
#
# **Example 3:**
#
#
#
#     Input: sequence = "ababc", word = "ac"
#     Output: 0
#     Explanation: "ac" is not a substring in "ababc".
#
#
#
#
# **Constraints:**
#
#   * `1 <= sequence.length <= 100`
#   * `1 <= word.length <= 100`
#   * `sequence` and `word` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Compat

@static if Compat.Sys.iswindows()
    using WinRPM
    WinRPM.install("libxml2-2", yes=true)
end


text:
@reset
@set printstyle none

Econv = [
-76.335767822597347
-56.422272522003723
-231.188695053088594
-154.616795317070142
-114.180708994251702
-279.415437830677774
-40.448675124014166
-378.532023978714108
-228.312964998766091
-109.342963378094467
]

Edf = [
-76.33577146094731
-56.42227793333907
-231.18868648219075
-154.61678476272598
-114.18069925474244
-279.41549876125123
-40.44867306182801
-378.53202120936407
-228.31293477911845
-109.34269895623703
]

@testset "CCSD" begin
    @testset "Conventional" begin
        Fermi.Options.set("df", false)

        for i = eachindex(molecules)
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*molecules[i]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", basis[i])

            wf = @energy ccsd
            @test isapprox(wf.energy, Econv[i], rtol=tol) # Energy from Psi4
        end
    end

    @testset "Density Fitted" begin
        Fermi.Options.set("df", true)
        Fermi.Options.set("jkfit", "cc-pvqz-jkfit")
        Fermi.Options.set("rifit", "cc-pvqz-rifit")

        for i = eachindex(molecules)
            # Read molecule
            path = joinpath(@__DIR__, "xyz/"*molecules[i]*".xyz")
            mol = open(f->read(f,String), path)

            # Define options
            Fermi.Options.set("molstring", mol)
            Fermi.Options.set("basis", basis[i])

            wf = @energy ccsd
            @test isapprox(wf.energy, Edf[i], rtol=tol) # Energy from Psi4
        end
    end
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using Test
include("simpson_integrate.jl")

testfunc(x) = 4.0x^2 + 2.0x - 1.0
stepsize = 1e-2
xrange = 0:stepsize:1
funcvals = testfunc.(xrange)
integral = simpson_integrate(funcvals,stepsize)
exactintegral = 4.0/3
@test integral ≈ exactintegral

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
"""
    StochasticProcesses

A lightweight, efficient package for simulating stochastic processes on various domains.

##### Abstract types
- AbstractMarkovChain
- AbstractStochasticProcess

##### Concrete types
- MarkovChain
- MultiRandomWalk
- PoissonProcess
- RandomWalk

##### Methods
- gridwalk
- nexpectedarrivals
- nwalkers
- plot
- rand
- randomwalkmc
"""
module StochasticProcesses

export

# abstract types
AbstractMarkovChain, AbstractStochasticProcess,

# concrete types
MarkovChain, MultiRandomWalk, PoissonProcess, RandomWalk,

# methods
gridwalk, nexpectedarrivals, nwalkers, plot, rand, randomwalkmc

# include submodules
include("Samplers/Samplers.jl")
include("MarkovChains/MarkovChains.jl")
include("RandomWalks/RandomWalks.jl")
include("GridWalks/GridWalks.jl")
include("PoissonProcesses/PoissonProcesses.jl")
include("PlotStochasticProcesses/PlotStochasticProcesses.jl")

# use submodules
using .Samplers
using .MarkovChains
using .RandomWalks
using .GridWalks
using .PoissonProcesses
using .PlotStochasticProcesses
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
using DataDeps

register(DataDep(
  "balloons-b",
  "http://archive.ics.uci.edu/ml/datasets/Balloons",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/balloons/adult+stretch.data",
  "56f8efed1e7f2e5e21a7f78b1cb2df75860642667b145330b3a0ea196b2e3274",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=5,
        feature_indices=1:4,
        categoric_indices=1:4,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using PtFEM

ProjDir = dirname(@__FILE__)

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 4, 8, 4, 4, :y, Quadrilateral(8, 2), false),
  :properties => [100.0 1.0e5 0.3;],
  :x_coords => [0.0, 1.0, 2.0, 3.0, 4.0, 5.5, 7.0, 9.0, 12.0],
  :y_coords => [0.0, -1.25, -2.5, -3.75, -5.0],
  :support => [
    (  1, [0 1]), (  2, [0 1]), (  3, [0 1]), (  4, [0 1]), (  5, [0 1]), (  6, [0 1]),
    (  7, [0 1]), (  8, [0 1]), (  9, [0 0]), ( 14, [0 0]), ( 23, [0 0]), ( 28, [0 0]),
    ( 37, [0 0]), ( 42, [0 0]), ( 51, [0 0]), ( 56, [0 0]), ( 65, [0 0]), ( 70, [0 0]),
    ( 79, [0 0]), ( 84, [0 0]), ( 93, [0 0]), ( 98, [0 0]), (107, [0 0]), (112, [0 0]),
    (113, [0 0]), (114, [0 0]), (115, [0 0]), (116, [0 0]), (117, [0 0]), (118, [0 0]),
    (119, [0 0]), (120, [0 0]), (121, [0 0])
  ],
  :loaded_nodes => [
    ( 1, [0.0 -0.166667]), (10, [0.0 -0.666667]), (15, [0.0 -0.333333]),
    (24, [0.0 -0.666667]), (29, [0.0 -0.166667])
  ],
  :tol => 0.001,
  :limit => 250,
  :qincs => [200.0, 100.0, 50.0, 50.0, 50.0, 30.0, 20.0, 10.0, 5.0, 1.0, 0.5],
  :cg_tol => 0.0001,
  :cg_limit => 100
)

@time m = p62a(data)
println()
m |> display


text:
## ## parametric copula without rotation
## ##-----------------------------------

## function rand(cop::ParamPC)
##     fam, params = getFamAndParams(cop)

##     @mput fam params
##     simVal = @matlab PairCopulaRand(fam, 1, params)
##     return simVal
## end

## ## multiple values
## function rand(cop::ParamPC, n::Int)
##     fam, params = getFamAndParams(cop)
##     nFloat = float(n)

##     @mput fam params nFloat
##     simVal = @matlab PairCopulaRand(fam, nFloat, params)
##     return simVal
## end

## ## parametric copula with rotation
## ##--------------------------------

## function rand(cop::ModPC)
##     return rand(cop.cop, cop.mod)
## end

## function rand(cop::ParamPC, rot::Rotation)
##     fam, params = getFamAndParams(cop)
##     degr = rot.degrees

##     @mput fam params degrees
##     simVal = @matlab PairCopulaRand(fam, 1, params, degr)
##     return simVal
## end

## ## multiple values
## function rand(cop::ModPC, n::Int)
##     return rand(cop.cop, cop.mod, n)
## end

## function rand(cop::ParamPC, rot::Rotation, n::Int)
##     fam, params = getFamAndParams(cop)
##     nFloat = float(n)
##     degr = rot.degrees

##     @mput fam params nFloat degr
##     simVal = @matlab PairCopulaRand(fam, nFloat, params, degr)
##     return simVal
## end


text:
import .AST: ASTMDP, ASTMetrics

"""
Collection of failure related metrics.
- iteration of first failure
- number of failures
- failure rate
- highest log-likelihood of failure
"""
@with_kw struct FailureMetrics
    num_terminals = NaN
    first_failure = NaN
    num_failures = NaN
    failure_rate = NaN
    highest_loglikelihood = NaN
end

"""
Calculate highest log-likelihood of failure.
"""
highest_loglikelihood_of_failure(planner) = highest_loglikelihood_of_failure(planner.mdp.metrics)
function highest_loglikelihood_of_failure(metrics::ASTMetrics)
    # TODO: get this index from the returned `action_trace` itself
    failure_loglikelihood = NaN
    if any(metrics.event)
        # Failures were found.
        failure_loglikelihood = maximum(metrics.logprob[metrics.event])
    end
    return failure_loglikelihood
end

"""
    failure_metrics(planner)
    failure_metrics(mdp::ASTMDP)
    failure_metrics(metrics::ASTMetrics)

Collect failure metrics including:
- iteration of first failure
- number of failures
- failure rate
- highest log-likelihood of failure
"""
failure_metrics(planner; verbose=true) = failure_metrics(planner.mdp.metrics; verbose=verbose)
failure_metrics(mdp::ASTMDP; verbose=true) = failure_metrics(mdp.metrics; verbose=verbose)
function failure_metrics(metrics::ASTMetrics; verbose=true)
    E = metrics.event
    terminal_idx = findall(metrics.terminal)

    if isempty(terminal_idx)
        verbose && @warn "No terminal states. Something may be incorrect in simulation."
    else
        num_terminals = length(terminal_idx)
        if findfirst(E) === nothing
            verbose && @info "No failures recorded."
            return FailureMetrics(num_terminals, NaN, 0, 0, 0)
        else
            possible_failures = metrics.event[terminal_idx] # 0 or 1 if the terminal state was a failure.
            first_failure = findfirst(possible_failures)
            num_failures = sum(possible_failures)
            failure_rate = num_failures/num_terminals * 100
            highest_loglikelihood = highest_loglikelihood_of_failure(metrics)
            failure_metrics = FailureMetrics(num_terminals, first_failure, num_failures, failure_rate, highest_loglikelihood)
            return failure_metrics
        end
    end
end

"""
    print_metrics(planner)
    print_metrics(mdp::ASTMDP)
    print_metrics(metrics::ASTMetrics)

Print failure metrics including:
- iteration of first failure
- number of failures
- failure rate
- highest log-likelihood of failure
"""
print_metrics(planner) = print_metrics(planner.mdp.metrics)
print_metrics(mdp::ASTMDP) = print_metrics(mdp.metrics)
function print_metrics(metrics::ASTMetrics)
    fail_metrics = failure_metrics(metrics)
    if fail_metrics isa FailureMetrics
        println("First failure: ", fail_metrics.first_failure, " of ", fail_metrics.num_terminals)
        println("Number of failures: ", fail_metrics.num_failures)
        println("Failure rate: ", round(fail_metrics.failure_rate, digits=5), "%")
        println("Highest log-likelihood of failure: ", round(fail_metrics.highest_loglikelihood, digits=5))
    end
    return fail_metrics
end

"""
Display failure metrics in a LaTeX enviroment.
Useful in Pluto.jl notebooks.
"""
function latex_metrics(metrics::FailureMetrics; digits=5)
    # Note indenting is important here to render correctly.
    return Markdown.parse(string("
\$\$\\begin{align}",
"p(\\text{fail}) &=", round(metrics.failure_rate/100, digits=digits), "\\\\",
"\\text{first failure index} &= ", metrics.first_failure, "\\text{ of }", metrics.num_terminals, "\\\\",
"\\text{highest log-likelihood of failure} &= ", round(metrics.highest_loglikelihood, digits=digits), "\\\\",
"\\end{align}\$\$"))
end

function latex_metrics(metrics_mean::FailureMetrics, metrics_std::FailureMetrics; digits=5)
    # Note indenting is important here to render correctly.
    return Markdown.parse(string("
\$\$\\begin{align}",
"p(\\text{fail}) &=", round(metrics_mean.failure_rate/100, digits=digits), " \\pm ", round(metrics_std.failure_rate/100, digits=digits), "\\\\",
"\\text{first failure index} &= ", round(metrics_mean.first_failure, digits=digits), " \\pm ", round(metrics_std.first_failure, digits=digits), "\\\\",
"\\text{num. episodes} &= ", round(metrics_mean.num_terminals, digits=digits), " \\pm ", round(metrics_std.num_terminals, digits=digits), "\\\\",
"\\text{highest log-likelihood of failure} &= ", round(metrics_mean.highest_loglikelihood, digits=digits), " \\pm ", round(metrics_std.highest_loglikelihood, digits=digits), "\\\\",
"\\end{align}\$\$"))
end

function Statistics.mean(fms::Vector{FailureMetrics})
    num_terminals = mean([m.num_terminals for m in fms])
    first_failure = mean(filter(!isnan, [m.first_failure for m in fms]))
    num_failures = mean([m.num_failures for m in fms])
    failure_rate = mean([m.failure_rate for m in fms])
    highest_loglikelihood = mean([m.highest_loglikelihood for m in fms])

    return FailureMetrics(
        num_terminals,
        first_failure,
        num_failures,
        failure_rate,
        highest_loglikelihood)
end

function Statistics.std(fms::Vector{FailureMetrics})
    num_terminals = std([m.num_terminals for m in fms])
    first_failure = std(filter(!isnan, [m.first_failure for m in fms]))
    num_failures = std([m.num_failures for m in fms])
    failure_rate = std([m.failure_rate for m in fms])
    highest_loglikelihood = std([m.highest_loglikelihood for m in fms])

    return FailureMetrics(
        num_terminals,
        first_failure,
        num_failures,
        failure_rate,
        highest_loglikelihood)
end


text:
"""
plots a timeseries
"""
function plot(jl_dateTimes::Array{DateTime},y;options=Dict())
    #Make everything a string
    dateTimeStrings = string.(jl_dateTimes)
    #Not sure if this line is really necessary, but change format to matlab format
    dateTimeStrings = replace.(dateTimeStrings,"T"=>" ")
    lengthDTStrings = length.(dateTimeStrings)
    #Ensure Consistent formatting
    locationsOfShortStrings = findall(lengthDTStrings.==19)
    dateTimeStrings[locationsOfShortStrings].=dateTimeStrings[locationsOfShortStrings].*".000"
    MATLAB.@mput dateTimeStrings
    #Convert to MATLAB DateTime
    MATLAB.eval_string("x = datetime(dateTimeStrings,\"InputFormat\",\"yyyy-MM-dd HH:mm:ss.SSS\")")
    MATLAB.@mput y
    #Plot!
    optionsString = namevaluepairs(options)
    MATLAB.eval_string("plot(x,y"*optionsString*")")
end

"""
plots a timeseries
"""
function plot(jl_dateTimes::Array{DateTime},y,lineSpec;options=Dict())
    #Make everything a string
    dateTimeStrings = string.(jl_dateTimes)
    #Not sure if this line is really necessary, but change format to matlab format
    dateTimeStrings = replace.(dateTimeStrings,"T"=>" ")
    lengthDTStrings = length.(dateTimeStrings)
    #Ensure Consistent formatting
    locationsOfShortStrings = findall(lengthDTStrings.==19)
    dateTimeStrings[locationsOfShortStrings].=dateTimeStrings[locationsOfShortStrings].*".000"
    MATLAB.@mput dateTimeStrings
    #Convert to MATLAB DateTime
    MATLAB.eval_string("x = datetime(dateTimeStrings,\"InputFormat\",\"yyyy-MM-dd HH:mm:ss.SSS\")")
    MATLAB.@mput y
    #Plot!
    optionsString = namevaluepairs(options)
    MATLAB.eval_string("plot(x,y"*",\"$lineSpec\""*optionsString*")")
end

"""
plots a timeseries
"""
function plot_ts(jl_dateTimes,y,lineSpec;options=Dict())
    #Make everything a string
    dateTimeStrings = string.(jl_dateTimes)
    #Not sure if this line is really necessary, but change format to matlab format
    dateTimeStrings = replace.(dateTimeStrings,"T"=>" ")
    lengthDTStrings = length.(dateTimeStrings)
    #Ensure Consistent formatting
    locationsOfShortStrings = findall(lengthDTStrings.==19)
    dateTimeStrings[locationsOfShortStrings].=dateTimeStrings[locationsOfShortStrings].*".000"
    MATLAB.@mput dateTimeStrings
    #Convert to MATLAB DateTime
    MATLAB.eval_string("x = datetime(dateTimeStrings,\"InputFormat\",\"yyyy-MM-dd HH:mm:ss.SSS\")")
    MATLAB.@mput y
    #Plot!
    optionsString = namevaluepairs(options)
    MATLAB.eval_string("plot(x,y"*",\"$lineSpec\""*optionsString*")")
end


text:
module KnuthMorrisPratt
  export KMP_search

  #=
    @fn     word search using KMP
    @param  [1] (text::String) run text to search
            [2] (word::String) search word
    @return (Int64) Found     -> position of word in text
                    Not Found -> length(text) + 1
  =#
  function KMP_search(text::String, word::String) :: Int64
    table = generate_KMP_table(word)
    p = 0
    i = 1
    # search
    while p+i <= length(text)
      if word[i] == text[p+i]
        i += 1
        if i == length(word)+1
          return p+1
        end
      else
        p += i - table[i]
        if i > 1
          i = table[i]
        end
      end
    end
    # not found
    return length(text)+1
  end

  #=
    @fn     Generate partial match table for KMP
    @param  (word::String) search word
    @return (Array{Int64, 1}) partial match table
  =#
  function generate_KMP_table(word::String) :: Array{Int64, 1}
    table = ones(Int64, length(word))
    table[1] = 0
    i = 3
    j = 1
    while i <= length(word)
      if word[i-1] == word[j]
        table[i] = j+1
        i += 1
        j += 1
      elseif j > 1
        j = table[j]
      else
        table[i] = 1
        i += 1
      end
    end
    return table
  end

end


text:
## -----
## refinement

# For [a,b] an isolating interval for a square free polynomial
# p(a) and p(b) have different signs.
# this finds a refined isolating interval [a′, b′] with |b′ - a′| < 2^{-L}
# if found within 64 refinement steps.
function refine_interval(p, a, b, L=min(53,maximum(precision, (a,b))))
    n = length(p) - 1
    n′ = sum(divrem(n,2))

    p′ = [i*p[i+1] for i ∈ 1:n]

    ϵ = epsᴸ(L)
    I = 𝐈(a,b)
    N = 1
    bisect_depth = 0

    while true

        if bisect_depth > 64
            return I
        end

        val, J = newtontest(p, p′, I, N)
        if val
            diff(J) < ϵ  && return J
            λ = diff(I)/diff(J)
            I = J
            N += 1
            continue
        end

        val, J = boundarytest(p, I, N)
        if val
            diff(J) < ϵ  && return J
            I = J
            continue
        end

        Iₗ, Iᵣ = bisect_interval(p, I, n′, L)

        val = signtest(p, Iₗ)
        if val == nothing
            val = signtest(p, Iᵣ)
            val != nothing && (val = !val)
        end

        if val == nothing
            dₗ,dᵣ = diff(Iₗ), diff(Iᵣ)
            if dₗ <= dᵣ
                dₗ ≤ ϵ && return Iₗ
            else
                dᵣ ≤ ϵ && return Iᵣ
            end
            error("I did not cleanly divide")
        end

        I = val ? Iₗ : Iᵣ
        diff(I) ≤ ϵ && return I
        N = max(0, N-1)
        bisect_depth += 1
    end
end

# refine intervals to size `L`, return midpoint
refine_roots(st::State) = midpoint.(refine_interval(st.p, I...) for I ∈ st)


text:
"""
    function c_sys_case4(; kwargs...)

Constructs the system for the example Case 4. The example based from [1], note that the
values of the elements are modified.

"""
function c_sys_case4(; kwargs...)
    sys_kwargs = filter_kwargs(; kwargs...)
    # Buses of the System
    nodes = nodes4_tcrd()
    # Build the System
    c_sys4_tcrd = PSY.System(
        100.0,
        nodes,
        thermal_generators4_tcrd(nodes),
        loads4_tcrd(nodes),
        branches4_tcrd(nodes);
        time_series_in_memory = get(sys_kwargs, :time_series_in_memory, true),
        sys_kwargs...,
    )

    # Dates for the System
    dates = collect(
        DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
        "1/1/2024  23:00:00",
        "d/m/y  H:M:S",
        ),
    )

    # Load pattern for 24 hrs
    #time series per zone
    loadz1_ts = ones(Float64, 24)

    # Loads of the system
    loads4_tcrd(nodes4_tcrd())
    timeseries_DA4_tcrd = [
        TimeArray(dates, loadz1_ts),
        TimeArray(dates, loadz1_ts),
        ]

    # Forecasts of the System
    if get(kwargs, :add_forecasts, true)
        forecast_data = SortedDict{Dates.DateTime, TimeSeries.TimeArray}()
        for (ix, l) in enumerate(PSY.get_components(PowerLoad, c_sys4_tcrd))
            ini_time = TimeSeries.timestamp(timeseries_DA4_tcrd[ix])[1]
            forecast_data[ini_time] = timeseries_DA4_tcrd[ix]
            PSY.add_time_series!(
                c_sys4_tcrd,
                l,
                PSY.Deterministic("max_active_power", forecast_data),
                )
        end
    end
    return c_sys4_tcrd
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LSODA_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LSODA")
JLLWrappers.@generate_main_file("LSODA", UUID("aae0fff6-70f8-51e1-979a-8150c7cc74bd"))
end  # module LSODA_jll


text:
export AbstractBag
"""
    AbstractBag{BT, N}<:TagBlock{BT, N}

Abstract `Bag` is a wrapper of a block that conserves all properties.
Including `mat`, `apply!`, `ishermitian`, `isreflexive`, `isunitary`,
`occupied_locs`, `apply_back!` and `mat_back!`.
"""
abstract type AbstractBag{BT, N}<:TagBlock{BT, N} end

Yao.mat(::Type{T}, bag::AbstractBag{N}) where {T,N} = mat(T, content(bag))
Yao.apply!(reg::AbstractRegister, bag::AbstractBag) = apply!(reg, content(bag))
Yao.ishermitian(bag::AbstractBag) = ishermitian(content(bag))
Yao.isreflexive(bag::AbstractBag) = isreflexive(content(bag))
Yao.isunitary(bag::AbstractBag) = isunitary(content(bag))
Yao.occupied_locs(bag::AbstractBag) = occupied_locs(content(bag))

function Yao.AD.apply_back!(state, b::AbstractBag, collector)
    Yao.AD.apply_back!(state, content(b), collector)
end
function Yao.AD.mat_back!(::Type{T}, b::AbstractBag, adjy, collector) where T
    Yao.AD.mat_back!(T, content(b), adjy, collector)
end

export Bag, enable_block!, disable_block!, setcontent!, isenabled
"""
    Bag{N}<:TagBlock{AbstractBlock, N}

A bag is a trivil container, but can
    * `setcontent!(bag, content)`
    * `disable_block!(bag)`
    * `enable_block!(bag)`
"""
mutable struct Bag{N}<:AbstractBag{AbstractBlock, N}
    content::AbstractBlock{N}
    mask::Bool
end
Bag(b::AbstractBlock) = Bag(b,true)

Yao.content(bag::Bag{N}) where N = bag.mask ? bag.content : put(N, 1=>I2)
Yao.chcontent(bag::Bag, content) = Bag(content)
setcontent!(bag::Bag, content) = (bag.content = content; bag)
disable_block!(b::Bag) = (b.mask = false; b)
enable_block!(b::Bag) = (b.mask = true; b)
isenabled(b::Bag) = b.mask

function Yao.print_annotation(io::IO, bag::Bag)
    printstyled(io, isenabled(bag) ? "[⊙] " : "[⊗] "; bold=true, color=isenabled(bag) ? :green : :red)
end

function Base.show(io::IO, ::MIME"plain/text", blk::Bag)
    return print_tree(io, blk; title=false, compact=false)
end


text:
module Message

struct Msgs
  start_stmt::String
  check_drinks_stmt::String
  mode_statement::String
  mode1_statement::String
  mode2_statement::String

end

const start_st = "Please enter multiple drinks in chronological order. The order can be in any order. The format should be 'YYYY/MM/DD H:M', such as '2020/12/01 21:34'."
const check_drinks_st = "How many drinks do you want to calculate?"
const mode_stmt = "Mode 1 is calculated by the amount of caffeine per 100ml, Mode 2 is calculated by the total amount of the caffeine.\n Choose 1 or 2."
const mode1_stmt = "Enter first the amount of the drink and then the amount of caffeine per 100 ml, separated by a space"
const mode2_stmt = "Enter the total amount of caffeine in mg."

const messages = Msgs(start_st, check_drinks_st, mode_stmt, mode1_stmt, mode2_stmt)

export Msgs

end

max_stars_repo_name: Kaikei-e/caffeCalJuli

max_issues_repo_name: Kaikei-e/caffeCalJuli

max_forks_repo_name: Kaikei-e/caffeCalJuli