text:
# FibonacciAnyons
"""
    struct FibonacciAnyon <: Sector
    FibonacciAnyon(s::Symbol)

Represents the anyons (isomorphism classes of simple objects) of the Fibonacci fusion
category. It can take two values, corresponding to the unit object `FibonacciAnyon(:I)`
and the non-trivial object `FibonacciAnyon(:τ)` with fusion rules ``τ ⊗ τ = 1 ⊕ τ``.
"""
struct FibonacciAnyon <: Sector
    isone::Bool
    function FibonacciAnyon(s::Symbol)
        s in (:I, :τ, :tau) || throw(ArgumentError("Unknown FibonacciAnyon $s."))
        new(s === :I)
    end
end

Base.iterate(::SectorValues{FibonacciAnyon}, i = 0) =
    i == 0 ? (FibonacciAnyon(:I), 1) : (i == 1 ? (FibonacciAnyon(:τ), 2) : nothing)

Base.IteratorSize(::Type{SectorValues{FibonacciAnyon}}) = HasLength()
Base.length(::SectorValues{FibonacciAnyon}) = 2

function Base.getindex(S::SectorValues{FibonacciAnyon}, i)
    if i == 1
        return FibonacciAnyon(:I)
    elseif i == 2
        return FibonacciAnyon(:τ)
    else
        throw(BoundsError(S, i))
    end
end
findindex(::SectorValues{FibonacciAnyon}, s::FibonacciAnyon) = 2 - s.isone

Base.hash(a::FibonacciAnyon, h::UInt) = hash(a.isone, h)
Base.convert(::Type{FibonacciAnyon}, s::Symbol) = FibonacciAnyon(s)

Base.one(::Type{FibonacciAnyon}) = FibonacciAnyon(:I)
Base.conj(s::FibonacciAnyon) = s
Base.isreal(::Type{FibonacciAnyon}) = false
Base.isless(a::FibonacciAnyon, b::FibonacciAnyon) = isless(!a.isone, !b.isone)

FusionStyle(::Type{FibonacciAnyon}) = SimpleFusion()

⊗(a::FibonacciAnyon, b::FibonacciAnyon) = FibonacciIterator(a, b)

struct FibonacciIterator
    a::FibonacciAnyon
    b::FibonacciAnyon
end
Base.IteratorSize(::Type{FibonacciIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{FibonacciIterator}) = Base.HasEltype()
Base.length(iter::FibonacciIterator) = (isone(iter.a) || isone(iter.b)) ? 1 : 2
Base.eltype(::Type{FibonacciIterator}) = FibonacciAnyon
function Base.iterate(iter::FibonacciIterator, state = 1)
    I = FibonacciAnyon(:I)
    τ = FibonacciAnyon(:τ)
    if state == 1
        iter.a == I && return (iter.b, 2)
        iter.b == I && return (iter.a, 2)
        return (I, 2)
    elseif state == 2
        (iter.a == iter.b == τ) && return (τ, 3)
        return nothing
    else
        return nothing
    end
end

Nsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon) =
    isone(a) + isone(b) + isone(c) != 2

function Fsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon,
                 d::FibonacciAnyon, e::FibonacciAnyon, f::FibonacciAnyon)
    Nsymbol(a, b, e) || return zero(_goldenratio)
    Nsymbol(e, c, d) || return zero(_goldenratio)
    Nsymbol(b, c, f) || return zero(_goldenratio)
    Nsymbol(a, f, d) || return zero(_goldenratio)

    I = FibonacciAnyon(:I)
    τ = FibonacciAnyon(:τ)
    if a == b == c == d == τ
        if e == f == I
            return +1/_goldenratio
        elseif e == f == τ
            return -1/_goldenratio
        else
            return +1/sqrt(_goldenratio)
        end
    else
        return one(_goldenratio)
    end
end

const _goldenratio = Float64(MathConstants.golden)

dim(a::FibonacciAnyon) = isone(a) ? one(_goldenratio) : _goldenratio

BraidingStyle(::Type{FibonacciAnyon}) = Anyonic()

function Rsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon)
    Nsymbol(a, b, c) || return 0*cis(0π/1)
    if isone(a) || isone(b)
        return cis(0π/1)
    else
        return isone(c) ? cis(4π/5) : cis(-3π/5)
    end
end

function Base.show(io::IO, a::FibonacciAnyon)
    s = isone(a) ? ":I" : ":τ"
    return get(io, :typeinfo, nothing) === FibonacciAnyon ?
        print(io, s) : print(io, "FibonacciAnyon(", s, ")")
end

# IsingAnyons
"""
    struct IsingAnyon <: Sector
    IsingAnyon(s::Symbol)

Represents the anyons (isomorphism classes of simple objects) of the Ising fusion category.
It can take three values, corresponding to the unit object `IsingAnyon(:I)` and the
non-trivial objects `IsingAnyon(:σ)` and `IsingAnyon(:ψ)`, with fusion rules
``ψ ⊗ ψ = 1``, ``σ ⊗ ψ = σ``, and ``σ ⊗ σ = 1 ⊕ ψ``.
"""
struct IsingAnyon <: Sector
    s::Symbol
    function IsingAnyon(s::Symbol)
        s == :sigma && (s = :σ)
        s == :psi && (s = :ψ)
        if !(s in (:I, :σ, :ψ))
            throw(ValueError("Unknown IsingAnyon $s."))
        end
        new(s)
    end
end

const all_isinganyons = (IsingAnyon(:I), IsingAnyon(:σ), IsingAnyon(:ψ))

Base.iterate(::SectorValues{IsingAnyon}, i = 1) = iterate(all_isinganyons, i)

Base.IteratorSize(::Type{SectorValues{IsingAnyon}}) = HasLength()
Base.length(::SectorValues{IsingAnyon}) = length(all_isinganyons)
Base.getindex(S::SectorValues{IsingAnyon}, i) = getindex(all_isinganyons, i)

function findindex(::SectorValues{IsingAnyon}, a::IsingAnyon)
    a == all_isinganyons[1] && return 1
    a == all_isinganyons[2] && return 2
    return 3
end

Base.hash(s::IsingAnyon, h::UInt) = hash(s.s, h)
Base.convert(::Type{IsingAnyon}, s::Symbol) = IsingAnyon(s)

Base.one(::Type{IsingAnyon}) = IsingAnyon(:I)
Base.conj(s::IsingAnyon) = s
Base.isreal(::Type{IsingAnyon}) = false

function Base.isless(a::IsingAnyon, b::IsingAnyon)
    vals = SectorValues{IsingAnyon}()
    return isless(findindex(vals, a), findindex(vals, b))
end

FusionStyle(::Type{IsingAnyon}) = SimpleFusion()

⊗(a::IsingAnyon, b::IsingAnyon) = IsingIterator(a, b)

struct IsingIterator
    a::IsingAnyon
    b::IsingAnyon
end

Base.IteratorSize(::Type{IsingIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{IsingIterator}) = Base.HasEltype()
Base.eltype(::Type{IsingIterator}) = IsingAnyon

function Base.length(iter::IsingIterator)
    σ = IsingAnyon(:σ)
    return (iter.a == σ && iter.b == σ) ? 2 : 1
end

function Base.iterate(iter::IsingIterator, state = 1)
    I, σ, ψ = all_isinganyons
    if state == 1
        iter.a == I && return (iter.b, 2)
        iter.b == I && return (iter.a, 2)
        iter.a == σ && iter.b == ψ && return (σ, 2)
        iter.a == ψ && iter.b == σ && return (σ, 2)
        return (I, 2)
    elseif state == 2
        (iter.a == iter.b == σ) && return (ψ, 3)
        return nothing
    else
        return nothing
    end
end

function Nsymbol(a::IsingAnyon, b::IsingAnyon, c::IsingAnyon)
    I, σ, ψ = all_isinganyons
    return ((a == I && b == c)
            || (b == I && a == c)
            || (c == I && a == b)
            || (a == σ && b == σ && c == ψ)
            || (a == σ && b == ψ && c == σ)
            || (a == ψ && b == σ && c == σ)
           )
end

function Fsymbol(a::IsingAnyon, b::IsingAnyon, c::IsingAnyon,
                 d::IsingAnyon, e::IsingAnyon, f::IsingAnyon)
    Nsymbol(a, b, e) || return 0.0
    Nsymbol(e, c, d) || return 0.0
    Nsymbol(b, c, f) || return 0.0
    Nsymbol(a, f, d) || return 0.0
    I, σ, ψ = all_isinganyons
    if a == b == c == d == σ
        if e == f == ψ
            return -1.0/sqrt(2.0)
        else
            return 1.0/sqrt(2.0)
        end
    end
    if e == f == σ
        if a == c == σ && b == d == ψ
            return -1.0
        elseif a == c == ψ && b == d == σ
            return -1.0
        end
    end
    return 1.0
end

dim(a::IsingAnyon) = a.s == :σ ? sqrt(2) : 1.0

BraidingStyle(::Type{IsingAnyon}) = Anyonic()

function Rsymbol(a::IsingAnyon, b::IsingAnyon, c::IsingAnyon)
    Nsymbol(a, b, c) || return complex(0.0)
    I, σ, ψ = all_isinganyons
    if c == I
        if b == a == σ
            return cis(-π/8)
        elseif b == a == ψ
            return complex(-1.0)
        end
    elseif c == σ && (a == σ && b == ψ || a == ψ && b == σ)
        return -1.0im
    elseif c == ψ && a == b == σ
        return cis(3π/8)
    end
    return complex(1.0)
end

function Base.show(io::IO, a::IsingAnyon)
    if get(io, :typeinfo, nothing) === IsingAnyon
        return print(io, ":$(a.s)")
    else
        return print(io, "IsingAnyon(:$(a.s))")
    end
end


text:
module Epochs

using LeapSeconds: offset_tai_utc

import Base: -, +, <, ==, isapprox, isless, show
import Dates
import Dates: format, parse

import ..AstroDates:
    Date,
    DateTime,
    Time,
    calendar,
    day,
    dayofyear,
    fractionofday,
    hour,
    j2000,
    julian,
    julian_twopart,
    millisecond,
    minute,
    month,
    second,
    secondinday,
    year,
    yearmonthday

export Epoch,
    CCSDS_EPOCH,
    FIFTIES_EPOCH,
    FUTURE_INFINITY,
    GALILEO_EPOCH,
    GPS_EPOCH,
    J2000_EPOCH,
    J2000_TO_JULIAN,
    J2000_TO_MJD,
    JULIAN_EPOCH,
    MODIFIED_JULIAN_EPOCH,
    PAST_INFINITY,
    UNIX_EPOCH,
    DateTime,
    Date,
    Time,
    day,
    dayofyear,
    fractionofday,
    hour,
    j2000,
    julian,
    julian_period,
    julian_twopart,
    millisecond,
    minute,
    modified_julian,
    month,
    now,
    second,
    secondinday,
    timescale,
    year,
    yearmonthday,
    -

using ..TimeScales
using ..AstroDates
using ..Periods

const J2000_TO_JULIAN = 2.451545e6days
const J2000_TO_MJD = 51544.5days

@inline function two_sum(a, b)
    hi = a + b
    a1 = hi - b
    b1 = hi - a1
    lo = (a - a1) + (b - b1)
    hi, lo
end

struct Epoch{S<:TimeScale, T} <: Dates.AbstractDateTime
    scale::S
    second::Int64
    fraction::T
    function Epoch{S}(second::Int64, fraction::T) where {S<:TimeScale, T}
        return new{S, T}(S(), second, fraction)
    end
end

Epoch{S,T}(ep::Epoch{S,T}) where {S,T} = ep

@inline function apply_offset(second::Int64, fraction, offset)
    sum, residual = two_sum(fraction, offset)
    if !isfinite(sum)
        fraction′ = sum
        second′ = sum < 0 ? typemin(Int64) : typemax(Int64)
    else
        int_secs = floor(Int64, sum)
        fraction′ = sum - int_secs + residual
        second′ = second + int_secs
    end
    return second′, fraction′
end

function Epoch{S}(ep::Epoch{S}, Δt) where {S<:TimeScale}
    second, fraction = apply_offset(ep.second, ep.fraction, Δt)
    Epoch{S}(second, fraction)
end

"""
    Epoch{S}(jd1::T, jd2::T=zero(T); origin=:j2000) where {S, T<:Period}

Construct an `Epoch` with time scale `S` from a Julian date
(optionally split into `jd1` and `jd2`). `origin` determines the
variant of Julian date that is used. Possible values are:

- `:j2000`: J2000 Julian date, starts at `2000-01-01T12:00`
- `:julian`: Julian date, starts at `-4712-01-01T12:00`
- `:modified_julian`: Modified Julian date, starts at `1858-11-17T00:00`

### Examples ###

```jldoctest
julia> Epoch{CoordinatedUniversalTime}(0.0days, 0.5days)
2000-01-02T00:00:00.000 UTC

julia> Epoch{CoordinatedUniversalTime}(2.451545e6days, origin=:julian)
2000-01-01T12:00:00.000 UTC
```
"""
function Epoch{S}(jd1::T, jd2::T=zero(T), args...; origin=:j2000) where {S, T<:Period}
    if jd2 > jd1
        jd1, jd2 = jd2, jd1
    end

    u = unit(jd1)

    if origin == :j2000
        # pass
    elseif origin == :julian
        jd1 -= u(J2000_TO_JULIAN)
    elseif origin == :modified_julian
        jd1 -= u(J2000_TO_MJD)
    else
        throw(ArgumentError("Unknown Julian epoch: $origin"))
    end

    jd1v = jd1 |> seconds |> value
    jd2v = jd2 |> seconds |> value

    sum, residual = two_sum(jd1v, jd2v)
    epoch = floor(Int64, sum)
    offset = (sum - epoch) + residual
    return Epoch{S}(epoch, offset)
end

"""
    julian_period([T,] ep::Epoch; origin=:j2000, scale=timescale(ep), unit=days)

Return the period since Julian Epoch `origin` within the time scale `scale` expressed in
`unit` for a given epoch `ep`. The result is a [`Period`](@ref) object by default.
If the type argument `T` is present, the result is converted to `T` instead.

### Example ###

```jldoctest
julia> ep = UTCEpoch(2018, 2, 6, 20, 45, 0.0)
2018-02-06T20:45:00.000 UTC

julia> julian_period(ep; scale=TAI)
6611.365011574074 days

julia> julian_period(ep; unit=years)
18.100929728496464 years

julia> julian_period(Float64, ep)
6611.364583333333
```
"""
function julian_period(ep::Epoch; origin=:j2000, scale=timescale(ep), unit=days)
    ep1 = Epoch(ep, scale)
    jd1 = unit(ep1.second * seconds)
    jd2 = unit(ep1.fraction * seconds)

    if origin == :j2000
        # pass
    elseif origin == :julian
        jd1 += unit(J2000_TO_JULIAN)
    elseif origin == :modified_julian
        jd1 += unit(J2000_TO_MJD)
    else
        throw(ArgumentError("Unknown Julian epoch: $origin"))
    end

    return jd2 + jd1
end

function julian_period(::Type{T}, ep::Epoch; kwargs...) where T
    jd = julian_period(ep; kwargs...)
    return T(value(jd))
end

"""
    j2000(ep)

Return the J2000 Julian Date for epoch `ep`.

### Example ###

```jldoctest
julia> j2000(UTCEpoch(2000, 1, 1, 12))
0.0 days
```
"""
j2000(ep::Epoch) = julian_period(ep)

"""
    julian(ep)

Return the Julian Date for epoch `ep`.

### Example ###

```jldoctest
julia> julian(UTCEpoch(2000, 1, 1, 12))
2.451545e6 days
```
"""
julian(ep::Epoch) = julian_period(ep; origin=:julian)

"""
    modified_julian(ep)

Return the Modified Julian Date for epoch `ep`.

### Example ###

```jldoctest
julia> modified_julian(UTCEpoch(2000, 1, 1, 12))
51544.5 days
```
"""
modified_julian(ep::Epoch) = julian_period(ep; origin=:modified_julian)

"""
    julian_twopart(ep)

Return the two-part Julian Date for epoch `ep`, which is a tuple consisting
of the Julian day number and the fraction of the day.

### Example ###

```jldoctest
julia> julian_twopart(UTCEpoch(2000, 1, 2))
(2.451545e6 days, 0.5 days)
```
"""
function julian_twopart(ep::Epoch)
    sec_in_days = ep.second / SECONDS_PER_DAY
    frac_in_days = ep.fraction / SECONDS_PER_DAY
    j2k1, j2k2 = divrem(sec_in_days, 1.0)
    jd1 = j2k1 * days + J2000_TO_JULIAN
    jd2 = (frac_in_days + j2k2) * days
    return jd1, jd2
end

include("offsets.jl")
include("accessors.jl")

show(io::IO, ep::Epoch) = print(io, DateTime(ep), " ", timescale(ep))

function Epoch{S}(date::Date, time::Time, args...) where S
    hr = hour(time)
    mn = minute(time)
    leap = getleap(S(), date)
    # We care only about discontinuities
    leap = ifelse(hr == 23 && mn == 59 && abs(leap) == 1.0, leap, 0.0)
    s, fraction = divrem(second(Float64, time) - leap, 1.0)
    daysec = Int64((j2000(date) - 0.5) * SECONDS_PER_DAY)
    hoursec = Int64(hour(time) * SECONDS_PER_HOUR)
    minutesec = Int64(minute(time) * SECONDS_PER_MINUTE)
    sec = Int64(s) + minutesec + hoursec + daysec
    return Epoch{S}(sec, fraction)
end

Dates.default_format(::Type{Epoch}) = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss ttt")

"""
    Epoch(str[, format])

Construct an `Epoch` from a string `str`. Optionally a `format` definition can
be passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
object or as a string. In addition to the character codes supported by `DateFormat` the character
code `D` is supported which is parsed as "day of year" (see the example below) and the character
code `t` which is parsed as the time scale.  The default format is `yyyy-mm-ddTHH:MM:SS.sss ttt`.

**Note:** Please be aware that this constructor requires that the time scale is part of `str`, e.g.
`2018-02-06T00:00 UTC`. Otherwise use an explicit constructor, e.g. `Epoch{UTC}`.

### Example ###

```jldoctest
julia> Epoch("2018-02-06T20:45:00.0 UTC")
2018-02-06T20:45:00.000 UTC

julia> Epoch("2018-037T00:00 UTC", "yyyy-DDDTHH:MM ttt")
2018-02-06T00:00:00.000 UTC
```
"""
Epoch(str::AbstractString, format::Dates.DateFormat=Dates.default_format(Epoch)) = parse(Epoch, str, format)

Epoch(str::AbstractString, format::AbstractString) = Epoch(str, Dates.DateFormat(format))

Dates.default_format(::Type{Epoch{S}}) where {S} = Dates.ISODateTimeFormat

"""
    Epoch{S}(str[, format]) where S

Construct an `Epoch` with time scale `S` from a string `str`. Optionally a `format` definition can
be passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
object or as a string. In addition to the character codes supported by `DateFormat` the code `D` can
be used which is parsed as "day of year" (see the example below).  The default format is
`yyyy-mm-ddTHH:MM:SS.sss`.

### Example ###

```jldoctest
julia> Epoch{CoordinatedUniversalTime}("2018-02-06T20:45:00.0")
2018-02-06T20:45:00.000 UTC

julia> Epoch{CoordinatedUniversalTime}("February 6, 2018", "U d, y")
2018-02-06T00:00:00.000 UTC

julia> Epoch{CoordinatedUniversalTime}("2018-037T00:00", "yyyy-DDDTHH:MM")
2018-02-06T00:00:00.000 UTC
```
"""
Epoch{S}(str::AbstractString,
         format::Dates.DateFormat=Dates.default_format(Epoch{S})) where {S} = parse(Epoch{S}, str, format)

Epoch{S}(str::AbstractString, format::AbstractString) where {S} = Epoch{S}(str, Dates.DateFormat(format))

Epoch{S}(d::Date) where {S} = Epoch{S}(d, AstroDates.H00)

Epoch{S}(dt::DateTime) where {S} = Epoch{S}(Date(dt), Time(dt))
Epoch{S}(dt::Dates.DateTime) where {S} = Epoch{S}(DateTime(dt))

"""
    now()

Get the current date and time as a `UTCEpoch`.
"""
now() = UTCEpoch(Dates.now())

"""
    Epoch{S}(year, month, day, hour=0, minute=0, second=0.0) where S

Construct an `Epoch` with time scale `S` from date and time components.

### Example ###

```jldoctest
julia> Epoch{CoordinatedUniversalTime}(2018, 2, 6, 20, 45, 0.0)
2018-02-06T20:45:00.000 UTC

julia> Epoch{CoordinatedUniversalTime}(2018, 2, 6)
2018-02-06T00:00:00.000 UTC
```
"""
function Epoch{S}(year::Int, month::Int, day::Int, hour::Int=0,
                  minute::Int=0, second::Float64=0.0, args...) where S
    Epoch{S}(Date(year, month, day), Time(hour, minute, second), args...)
end

function Epoch{S}(year::Int64, month::Int64, day::Int64, dayofyear::Int64,
                  hour::Int64, minute::Int64, second::Int64, milliseconds::Int64) where S
    if dayofyear != 0
        date = Date(year, dayofyear)
    else
        date = Date(year, month, day)
    end
    Epoch{S}(date, Time(hour, minute, second + 1e-3milliseconds))
end

function Epoch(year::Int64, month::Int64, day::Int64, dayofyear::Int64,
               hour::Int64, minute::Int64, second::Int64, milliseconds::Int64,
               scale::S) where S<:TimeScale
    if scale === TimeScales.NotATimeScale()
        throw(ArgumentError("Could not parse the provided string as an `Epoch`." *
                            " No time scale was provided."))
    end

    if dayofyear != 0
        date = Date(year, dayofyear)
    else
        date = Date(year, month, day)
    end
    Epoch{S}(date, Time(hour, minute, second + 1e-3milliseconds))
end

"""
    Epoch{S}(Δtai, ep::TAIEpoch) where S

Convert `ep`, a `TAIEpoch`, to an `Epoch` with time scale `S` by overriding
the offset between `S2` and `TAI` with `Δtai`.

### Examples ###

```jldoctest
julia> ep = TAIEpoch(2000,1,1)
2000-01-01T00:00:00.000 TAI

julia> TTEpoch(32.184, ep)
2000-01-01T00:00:32.184 TT
```
"""
function Epoch{S2}(offset, ep::Epoch{S1}) where {S1<:TimeScale, S2<:TimeScale}
    second, fraction = apply_offset(ep.second, ep.fraction, offset)
    Epoch{S2}(second, fraction)
end

"""
    Epoch{S2}(ep::Epoch{S1}) where {S1, S2}

Convert `ep`, an `Epoch` with time scale `S1`, to an `Epoch` with time
scale `S2`.

### Examples ###

```jldoctest
julia> ep = TTEpoch(2000,1,1)
2000-01-01T00:00:00.000 TT

julia> TAIEpoch(ep)
1999-12-31T23:59:27.816 TAI
```
"""
function Epoch{S2}(ep::Epoch{S1}) where {S1<:TimeScale, S2<:TimeScale}
    second, fraction = apply_offset(ep.second, ep.fraction, S1(), S2())
    Epoch{S2}(second, fraction)
end

"""
    Epoch(ep::Epoch{S1}, scale::S2) where {S1, S2}

Convert `ep`, an `Epoch` with time scale `S1`, to an `Epoch` with time
scale `S2`.

### Examples ###

```jldoctest
julia> ep = TTEpoch(2000,1,1)
2000-01-01T00:00:00.000 TT

julia> Epoch(ep, TAI)
1999-12-31T23:59:27.816 TAI
```
"""
function Epoch(ep::Epoch{S1}, scale::S2) where {S1<:TimeScale, S2<:TimeScale}
    second, fraction = apply_offset(ep.second, ep.fraction, S1(), S2())
    Epoch{S2}(second, fraction)
end

function Epoch{S2}(ep::Epoch{S1}, args...) where {S1<:TimeScale, S2<:TimeScale}
    offset = getoffset(S1(), S2(), ep.second, ep.fraction, args...)
    second, fraction = apply_offset(ep.second, ep.fraction, offset)
    Epoch{S2}(second, fraction)
end

Epoch{S}(ep::Epoch{S}) where {S<:TimeScale} = ep
Epoch(ep::Epoch{S}, ::S) where {S<:TimeScale} = ep

function isapprox(a::Epoch{S}, b::Epoch{S}; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps()) where S <: TimeScale
    a.second == b.second && isapprox(a.fraction, b.fraction; atol=atol, rtol=rtol)
end

function ==(a::Epoch, b::Epoch)
    a.second == b.second && a.fraction == b.fraction
end

<(ep1::Epoch, ep2::Epoch) = value(ep1 - ep2) < 0.0
isless(ep1::Epoch, ep2::Epoch) = isless(value(ep1 - ep2), 0.0)

+(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, value(seconds(p)))
-(ep::Epoch{S}, p::Period) where {S} = Epoch{S}(ep, -value(seconds(p)))

"""
    -(a::Epoch, b::Epoch)

Return the duration between epoch `a` and epoch `b`.

### Examples ###

```jldoctest
julia> UTCEpoch(2018, 2, 6, 20, 45, 20.0) - UTCEpoch(2018, 2, 6, 20, 45, 0.0)
20.0 seconds
```
"""
function -(a::Epoch{S}, b::Epoch{S}) where S<:TimeScale
    return ((a.second - b.second) + (a.fraction - b.fraction)) * seconds
end

# Generate aliases for all defined time scales so we can use
# e.g. `TTEpoch` instead of `Epoch{TT}`
for (scale, acronym) in zip(TimeScales.NAMES, TimeScales.ACRONYMS)
    epoch = Symbol(acronym, "Epoch")
    name = string(epoch)
    @eval begin
        const $epoch = Epoch{$scale}
        export $epoch

        """
            $($name)(str[, format])

        Construct a $($name) from a string `str`. Optionally a `format` definition can be
        passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
        object or as a string. In addition to the character codes supported by `DateFormat`
        the code `D` is supported which is parsed as "day of year" (see the example below).
        The default format is `yyyy-mm-ddTHH:MM:SS.sss`.

        ### Example ###

        ```jldoctest
        julia> $($name)("2018-02-06T20:45:00.0")
        2018-02-06T20:45:00.000 $($acronym)

        julia> $($name)("February 6, 2018", "U d, y")
        2018-02-06T00:00:00.000 $($acronym)

        julia> $($name)("2018-37T00:00", "yyyy-DDDTHH:MM")
        2018-02-06T00:00:00.000 $($acronym)
        ```
        """
        $epoch(::AbstractString)

        """
            $($name)(jd1::T, jd2::T=zero(T); origin=:j2000) where T<:Period

        Construct a $($name) from a Julian date (optionally split into
        `jd1` and `jd2`). `origin` determines the variant of Julian
        date that is used. Possible values are:

        - `:j2000`: J2000 Julian date, starts at `2000-01-01T12:00`
        - `:julian`: Julian date, starts at `-4712-01-01T12:00`
        - `:modified_julian`: Modified Julian date, starts at `1858-11-17T00:00`

        ### Examples ###

        ```jldoctest
        julia> $($name)(0.0days, 0.5days)
        2000-01-02T00:00:00.000 $($acronym)

        julia> $($name)(2.451545e6days, origin=:julian)
        2000-01-01T12:00:00.000 $($acronym)
        ```
        """
        $epoch(::Number, ::Number)

        """
            $($name)(year, month, day, hour=0, minute=0, second=0.0)

        Construct a $($name) from date and time components.

        ### Example ###

        ```jldoctest
        julia> $($name)(2018, 2, 6, 20, 45, 0.0)
        2018-02-06T20:45:00.000 $($acronym)

        julia> $($name)(2018, 2, 6)
        2018-02-06T00:00:00.000 $($acronym)
        ```
        """
        $epoch(::Int, ::Int, ::Int)
    end
end

include("range.jl")

const JULIAN_EPOCH = TTEpoch(AstroDates.JULIAN_EPOCH, AstroDates.H12)
const J2000_EPOCH = TTEpoch(AstroDates.J2000_EPOCH, AstroDates.H12)
const MODIFIED_JULIAN_EPOCH = TTEpoch(AstroDates.MODIFIED_JULIAN_EPOCH, AstroDates.H00)
const FIFTIES_EPOCH = TTEpoch(AstroDates.FIFTIES_EPOCH, AstroDates.H00)
const CCSDS_EPOCH = TTEpoch(AstroDates.CCSDS_EPOCH, AstroDates.H00)
const GALILEO_EPOCH = TTEpoch(AstroDates.GALILEO_EPOCH, AstroDates.H00)
const GPS_EPOCH = TTEpoch(AstroDates.GPS_EPOCH, AstroDates.H00)
const UNIX_EPOCH = TAIEpoch(AstroDates.UNIX_EPOCH, Time(0, 0, 10.0))

const PAST_INFINITY = TAIEpoch(UNIX_EPOCH, -Inf)
const FUTURE_INFINITY = TAIEpoch(UNIX_EPOCH, Inf)

const EPOCH_77 = TAIEpoch(1977, 1, 1)

function Dates.validargs(::Type{Epoch}, y::Int64, m::Int64, d::Int64,
                         h::Int64, mi::Int64, s::Int64, ms::Int64, ts::S) where S<:TimeScale
    err = Dates.validargs(Dates.DateTime, y, m, d, h, mi, s, ms)
    err !== nothing || return err
    return Dates.argerror()
end

abstract type DayOfYearToken end

@inline function Dates.tryparsenext(d::Dates.DatePart{'D'}, str, i, len, locale)
    next = Dates.tryparsenext_base10(str, i, len, 1, 3)
    next === nothing && return nothing
    val, i = next
    (val >= 1 && val <= 366) || throw(ArgumentError("Day number must be within 1 and 366."))
    return val, i
end

function Dates.format(io, d::Dates.DatePart{'D'}, ep)
    print(io, dayofyear(ep))
end

function Dates.format(io, d::Dates.DatePart{'t'}, ep)
    print(io, timescale(ep))
end

end

ext: jl


text:
using SafeTestsets
@safetestset "My f tests" begin include("myftests.jl") end


text:
function topologicalsort!(g::AbstractGraph)
    for v = 1:countvertices(g)
        setmark!(g, v, 0)
    end
    for v = 1:countvertices(g)
        if getmark(g, v) == 0
            _topologicalsort!(g, v)
        end
    end
end

function _topologicalsort!(g::AbstractGraph, v::Vertex)
    setmark!(g, v, 1)
    let w = first(g, v)
        while w <= countvertices(g)
            if getmark(g, w) == 0
                _topologicalsort!(g, w)
            end
            w = next(g, v, w)
        end
    end
    #postvisit(g, v)
    println(v)
end


text:
@symbol_func function q_DV(cur_reactor::AbstractReactor)
  cur_q = 18.31e-3

  cur_q *= 1 + cur_reactor.Q / 5.0

  cur_q *= 2 ^ 0.6

  cur_q /= ( 1.0 + cur_reactor.kappa_95 ^ 2 ) ^ 0.6

  cur_q *= P_H(cur_reactor)

  cur_q *= cur_reactor.I_P ^ 1.2

  cur_q /= cur_reactor.epsilon ^ 1.2

  cur_q /= cur_reactor.R_0 ^ 2.2

  cur_q
end


text:
n_draws = 100
grid = 0:10:500
wts = rand(10)
resampling = BinnedMeanWeightedResampling(grid, wts, n_draws)

@test resampling isa BinnedMeanWeightedResampling
@test typeof(resampling) <: AbstractBinnedSummarisedResampling


text:
module Colorsys

export rgb2yiq, yiq2rgb, rgb2hls, hls2rgb, rgb2hsv, hsv2rgb, yiq2hls, yiq2hsv, hls2yiq, hls2hsv, hsv2yiq, hsv2hls

include("CS.jl")

end # module


text:
using Random
using JuMP
using Ipopt
using MosekTools

function run()
    ind_others = [1,2]
    h = 4/2
    w = 2/2

    prob = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
    # prob = Model(with_optimizer(Mosek.Optimizer))
    ego_x,ego_y,ego_θ = 0, 0, 0

    oth_xs = zeros(length(ind_others))
    oth_ys = zeros(length(ind_others))
    oth_θs = zeros(length(ind_others))
    # for (i,ind) in enumerate(ind_others)
    #     oth_xs[i],oth_ys[i],oth_θs[i] = scene[ind].state.posG.x, scene[ind].state.posG.y, scene[ind].state.posG.θ
    # end

    # TEMP
    oth_xs = [0 6 6 0];
    oth_ys = [3 0 3 -3];
    oth_θs = [0 0 0 0];

    @variable(prob, ego_p); @variable(prob, ego_q) # oth_p: x-coord oth_q: y-coord
    @variable(prob, oth_p); @variable(prob, oth_q)
    @objective(prob, MOI.MIN_SENSE, (ego_p-oth_p)^2+(ego_q-oth_q)^2)
    # @objective(prob, Min, (ego_p-oth_p)^2+(ego_q-oth_q)^2)

    # constraints
    @constraint(prob, ((ego_p-ego_x)*cos(ego_θ) + (ego_q-ego_y)*sin(ego_θ))^2/h^2
                    + ((ego_p-ego_x)*sin(ego_θ) + (ego_q-ego_y)*cos(ego_θ))^2/w^2
                    <= 1.0
                )
    for (i,ind) in enumerate(ind_others)
        @constraint(prob, ((oth_p-oth_xs[i])*cos(oth_θs[i]) + (oth_q-oth_ys[i])*sin(oth_θs[i]))^2/h^2
                        + ((oth_p-oth_xs[i])*sin(oth_θs[i]) + (oth_q-oth_ys[i])*cos(oth_θs[i]))^2/w^2
                        <= 1.0
                    )
    end

    # @variable(prob, xᵢ′); @variable(prob, yᵢ′)
    # @constraint(prob, ((x′-x)*cos(θ) + (y′-y)*sin(θ))^2/(h^2) + ((x′-x)*sin(θ) + (y′-y)*cos(θ))^2/(w^2) <= 1.0)
    # @constraint(prob, ((xᵢ′-xᵢ)*cos(θᵢ) + (yᵢ′-yᵢ)*sin(θᵢ))^2/(h^2) + ((xᵢ′-xᵢ)*sin(θᵢ) + (yᵢ′-yᵢ)*cos(θᵢ))^2/(w^2) <= 1.0)
    optimize!(prob)

    return prob
end

prob = run()
println("solution:", objective_value(prob))


text:
# Overengineered set of strategy types and structs? Not gonna argue with that, but I do this for fun and sometimes I have a wierd idea of what fun is.
# Also, everything is the fault of size transparent vertices, especially SizeStack and those blasted Δnoutfactors!
"""
    AbstractConnectStrategy

Base type for strategies for how to (re)connect vertices when doing structural
mutation.
"""
abstract type AbstractConnectStrategy end
struct ConnectAll <: AbstractConnectStrategy end
struct ConnectNone <: AbstractConnectStrategy end

"""
    AbstractAlignSizeStrategy

Base type for strategies for how to align size (`nin`/`nout`) when doing structural mutation.

Note that all strategies are not guaranteed to work in all cases.

Default strategies should however be selected based on case so that things always work out.
"""
abstract type AbstractAlignSizeStrategy end

"""
    NoSizeChange <: AbstractAlignSizeStrategy
    NoSizeChange()

Don't do any size change.
"""
struct NoSizeChange <: AbstractAlignSizeStrategy end

"""
    ChangeNinOfOutputs <: AbstractAlignSizeStrategy
    ChangeNinOfOutputs(Δoutsize)

Just sets `nin` of each output to the provided value. Sometimes you just know the answer...
"""
struct ChangeNinOfOutputs <: AbstractAlignSizeStrategy
    Δoutsize
end

"""
    FailAlignSizeError <: AbstractAlignSizeStrategy
    FailAlignSizeError()

Throws an error.
"""
struct FailAlignSizeError <: AbstractAlignSizeStrategy end

"""
    FailAlignSizeWarn <: AbstractAlignSizeStrategy
    FailAlignSizeWarn()
    FailAlignSizeWarn(;andthen, msgfun)

Logs warning and then proceeds with the next action.
"""
struct FailAlignSizeWarn <: AbstractAlignSizeStrategy
    andthen
    msgfun
end
FailAlignSizeWarn(;andthen=FailAlignSizeRevert(), msgfun=(vin,vout) -> "Could not align sizes of $(vin) and $(vout)!") = FailAlignSizeWarn(andthen, msgfun)

"""
    FailAlignSizeRevert <: AbstractAlignSizeStrategy
    FailAlignSizeRevert()

Reverts new/removed edges (if any).

In other words, `create_edge!` and `remove_edge!` with `strategy = FailAlignSizeRevert` is a noop.

Note: Only works if input vertex is only input once to output vertex due to lazy coding.
For example, if `vout` has inputs `[v1,v2,v1]` and `remove_edge!(v1, vout, startegy=FailAlignSizeRevert)` is called, funtion will exit with graph in invalid state.
Same if `vout` has inputs `[v1,v2]` and `create_edge!(v1,vout)` is called.
"""
struct FailAlignSizeRevert <: AbstractAlignSizeStrategy end

"""
    AlignSizeBoth <: AbstractAlignSizeStrategy
    AlignSizeBoth()
    AlignSizeBoth(fallback)

Align sizes by changing both input and output considering any Δfactors.
Fallback to another strategy (default `FailAlignSizeError`) if size change is not possible.
"""
struct AlignSizeBoth <: AbstractAlignSizeStrategy
    fallback
end
AlignSizeBoth() = AlignSizeBoth(FailAlignSizeError())

"""
    DecreaseBigger <: AbstractAlignSizeStrategy
    DecreaseBigger()
    DecreaseBigger(fallback)

Try to align size by decreasing in the direction (in/out) which has the bigger size.
Fallback to another strategy (default `AlignSizeBoth`) if size change is not possible.
"""
struct DecreaseBigger <: AbstractAlignSizeStrategy
    fallback
end
DecreaseBigger() = DecreaseBigger(AlignSizeBoth())

"""
    IncreaseSmaller <: AbstractAlignSizeStrategy
    IncreaseSmaller()
    IncreaseSmaller(fallback)

Try to align size by increasing in the direction (in/out) which has the smaller size.
Fallback to another strategy (default `DecreaseBigger`) if size change is not possible.
"""
struct IncreaseSmaller <: AbstractAlignSizeStrategy
    fallback
end
IncreaseSmaller() = IncreaseSmaller(DecreaseBigger())

"""
    SelectOutputs <: AbstractAlignSizeStrategy
    SelectOutputs(;select=OutSelectExact(),align=IncreaseSmaller(), valuefun=v -> ones(nout_org(v)))

First align size using `align`, then select outputs (through `Δoutputs`) using `select` and `valuefun` if size alignment is successful.

Motivation is that selecting outputs is more efficient to do when original sizes are aligned.
"""
struct SelectOutputs <: AbstractAlignSizeStrategy
    selectstrategy::AbstractSelectionStrategy
    alignstrategy::AbstractAlignSizeStrategy
    valuefun::Function
end
SelectOutputs(;select=OutSelectExact(),align=IncreaseSmaller(), valuefun=v -> ones(nout_org(v))) = SelectOutputs(select, align, valuefun)

"""
    PostSelectOutputs <: AbstractAlignSizeStrategy
    PostSelectOutputs(;select=OutSelectExact(),align=IncreaseSmaller(), valuefun=v -> ones(nout_org(v)))

Post change alignment strategy which first aligns size using `align`, then select outputs (through `Δoutputs`) using `select` and `valuefun` if size alignment is successful. If `Δoutputs` is not successful, the `fallback` strategy will be invoked.

Motivation is basically convenience when creating new edges between vertices.
"""
struct PostSelectOutputs <: AbstractAlignSizeStrategy
    selectstrategy::AbstractSelectionStrategy
    alignstrategy::AbstractAlignSizeStrategy
    valuefun::Function
    fallback::AbstractAlignSizeStrategy
end
PostSelectOutputs(;select=OutSelectExact(),align=PostAlignJuMP(), valuefun=v -> ones(nout_org(v)), fallback=FailAlignSizeRevert()) = PostSelectOutputs(select, align, valuefun, fallback)

"""
    ApplyMutation <: AbstractAlignSizeStrategy
    ApplyMutation()
    ApplyMutation(strategy::AbstractAlignSizeStrategy)

First align size using `strategy`, then invoke `apply_mutation` if size alignment is successful.

Motivation is that selecting outputs is more efficient to do when original sizes are aligned.
"""
struct ApplyMutation <: AbstractAlignSizeStrategy
    strategy::AbstractAlignSizeStrategy
end
ApplyMutation() = ApplyMutation(SelectOutputs())

"""
    PostApplyMutation <: AbstractAlignSizeStrategy
    PostApplyMutation()
    PostApplyMutation(strategy::AbstractAlignSizeStrategy)

Post change alignment strategy which first aligns size using `strategy`, then invoke `apply_mutation` if size alignment is successful.

Motivation is basically convenience when creating new edges between vertices.
"""
struct PostApplyMutation <: AbstractAlignSizeStrategy
    strategy::AbstractAlignSizeStrategy
end
PostApplyMutation() = PostApplyMutation(PostSelectOutputs())

"""
    CheckNoSizeCycle <: AbstractAlignSizeStrategy
    CheckNoSizeCycle()
    CheckNoSizeCycle(;ifok, ifnok)

Check if a size change in one direction causes a change in the other direction and execute strategy `ifnok` (default `FailAlignSizeWarn`) if this is the case.
Motivation is that removing will result in the computation graph being in an invalid state as one of the vertices must fulfill the impossible criterion `nout(v) == nout(v) + a` where `a > 0`.

If no such cycle is detected, then proceed to execute strategy `ifok` (default `IncreaseSmaller`).

Will execute strategy `ifok` if vertex shall not to be removed.
"""
struct CheckNoSizeCycle <: AbstractAlignSizeStrategy
    ifok
    ifnok
end
CheckNoSizeCycle(;ifok=IncreaseSmaller(), ifnok=FailAlignSizeWarn(msgfun = (vin,vout) -> "Can not remove vertex $(vin)! Size cycle detected!")) = CheckNoSizeCycle(ifok, ifnok)

"""
    CheckCreateEdgeNoSizeCycle <: AbstractAlignSizeStrategy
    CheckCreateEdgeNoSizeCycle()
    CheckCreateEdgeNoSizeCycle(;ifok, ifnok)

Check if adding an edge creates the same type of size cycle that `CheckNoSizeCycle` checks for and execute `ifnok` (default `FailAlignSizeWarn`) if this is the case.
Motivation is that removing will result in the computation graph being in an invalid state as one of the vertices must fulfill the impossible criterion `nout(v) == nout(v) + a` where `a > 0`.

If no such cycle is detected, then proceed to execute strategy `ifok` (default `IncreaseSmaller`).

Will check both at `prealignsizes` (i.e before edge is added) and at `postalignsizes` (i.e after edge is added).
"""
struct CheckCreateEdgeNoSizeCycle <: AbstractAlignSizeStrategy
    ifok
    ifnok
end
CheckCreateEdgeNoSizeCycle(;ifok=IncreaseSmaller(), ifnok=FailAlignSizeWarn(msgfun = (vin,vout) -> "Can not add edge between $(vin) and $(vout)! Size cycle detected!")) = CheckCreateEdgeNoSizeCycle(ifok, ifnok)

"""
    CheckAligned <:AbstractAlignSizeStrategy
    CheckAligned()
    CheckAligned(ifnot)

Check if sizes are already aligned before making a change and return "go ahead" (`true`) if this is the case.
If not, proceed to execute another strategy (default `CheckNoSizeCycle`).
"""
struct CheckAligned <:AbstractAlignSizeStrategy
    ifnot
end
CheckAligned() = CheckAligned(CheckNoSizeCycle())

"""
    PostAlignJuMP <: AbstractAlignSizeStrategy
    PostAlignJuMP()
    PostAlignJuMP(s::AbstractAlignSizeStrategy)
    PostAlignJuMP(s::AbstractAlignSizeStrategy, fallback)

Align sizes using a `AbstractJuMPΔSizeStrategy`.

This is a post-align strategy, i.e it will be applied after a structural change has been made.
"""
struct PostAlignJuMP <: AbstractAlignSizeStrategy
    sizestrat::AbstractJuMPΔSizeStrategy
    fallback
end
PostAlignJuMP() = PostAlignJuMP(DefaultJuMPΔSizeStrategy())
PostAlignJuMP(s::AbstractJuMPΔSizeStrategy; fallback = FailAlignSizeError()) = PostAlignJuMP(AlignNinToNout(s, ΔSizeFailNoOp()), fallback)
PostAlignJuMP(s::AlignNinToNout; fallback=FailAlignSizeError()) = PostAlignJuMP(s, fallback)

"""
    RemoveStrategy
    RemoveStrategy()
    RemoveStrategy(rs::AbstractConnectStrategy)
    RemoveStrategy(as::AbstractAlignSizeStrategy)
    RemoveStrategy(rs::AbstractConnectStrategy, as::AbstractAlignSizeStrategy)

Strategy for removal of a vertex.

Consists of an `AbstractConnectStrategy` for how to treat inputs and outputs of
the removed vertex and an `AbstractAlignSizeStrategy` for how to align sizes of
inputs and outputs.
"""
struct RemoveStrategy
    reconnect::AbstractConnectStrategy
    align::AbstractAlignSizeStrategy
end
RemoveStrategy() = RemoveStrategy(ConnectAll(), CheckAligned())
RemoveStrategy(rs::AbstractConnectStrategy) = RemoveStrategy(rs, CheckAligned())
RemoveStrategy(as::AbstractAlignSizeStrategy) = RemoveStrategy(ConnectAll(), as)

"""
    remove!(v::MutationVertex, strategy=RemoveStrategy())

Removes `v` from the graph by removing it from its `inputs` and `outputs`.

It is possible to supply a strategy for how to 1) reconnect the inputs and outputs
of `v` and 2) align the input and output sizes of the `inputs` and `outputs` of `v`.

Default strategy is to first set `nin==nout` for `v` and then connect all its `inputs`
to all its `outputs`.
"""
function remove!(v::MutationVertex, strategy=RemoveStrategy())
    prealignsizes(strategy.align, v, vx -> vx == v) || return

    remove!(v, inputs, outputs, strategy.reconnect)
    remove!(v, outputs, inputs, strategy.reconnect)

    postalignsizes(strategy.align, v)
end

## Helper function to avoid code duplication. I don't expect this to be able to do
## anything useful unless f1=inputs and f2=outputs or vise versa.
function remove!(v::MutationVertex, f1::Function, f2::Function, s::AbstractConnectStrategy)
    for f1_v in f1(v)
        f1_f2_v = f2(f1_v)
        inds = findall(vx -> vx == v, f1_f2_v)
        deleteat!(f1_f2_v, inds)
        # Connect f2 of v to f2 of v1, e.g connect the inputs to v as inputs to v1 if f2 == inputs
        connect!(f1_v, f1_f2_v, inds, f2(v), s)
    end
end

"""
    connect!(v, to, inds, items, ::AbstractConnectStrategy)

Does connection of `items` to `to` at position `inds` depending on the given strategy.
"""
#Not sure broadcasting insert! like this is supposed to work, but it does...
connect!(v, to, inds, items, ::ConnectAll) = insert!.([to], inds, items)
function connect!(v, to, inds, items, ::ConnectNone)
    # Stupid generic function forces me to check this...
    # ... and no, it is not easy to understand what is going on here...
    if to == inputs(v) && !isempty(inds)
        rem_input!.(op.(v), inds...)
    end
end

tot_nin(v) = tot_nin(trait(v), v)
tot_nin(t::DecoratingTrait, v) = tot_nin(base(t), v)
tot_nin(::MutationTrait, v) = nin(v)[]
tot_nin(::SizeInvariant, v) = length(unique(nin(v))) == 1 ? unique(nin(v))[] : nothing
tot_nin(::SizeTransparent, v) = sum(nin(v))

# Boilerplate
prealignsizes(s::AbstractAlignSizeStrategy, v, will_rm::Function) = prealignsizes(s, v, v, will_rm)
prealignsizes(s::AbstractAlignSizeStrategy, vin, vout, will_rm) = true

# Failure cases
prealignsizes(::FailAlignSizeError, vin, vout, will_rm) = error("Could not align sizes of $(vin) and $(vout)!")
function prealignsizes(s::FailAlignSizeWarn, vin, vout, will_rm)
    @warn s.msgfun(vin, vout)
    return prealignsizes(s.andthen, vin, vout, will_rm)
end
prealignsizes(::FailAlignSizeRevert, vin, vout, will_rm) = false # No action needed?

# Actual actions
function prealignsizes(s::CheckAligned, vin, vout, will_rm)
    nout(vin) == tot_nin(vout) && return true
    return prealignsizes(s.ifnot, vin, vout, will_rm)
end

function prealignsizes(s::CheckNoSizeCycle, vin, vout, will_rm)
    if will_rm(vout) && vin==vout
        is_cyclic(ΔnoutSizeGraph(vin)) && return prealignsizes(CheckAligned(s.ifnok), vin, vout, will_rm)
    end
    return prealignsizes(s.ifok, vin, vout, will_rm)
end

function prealignsizes(s::CheckCreateEdgeNoSizeCycle, vin, vout, will_rm)
    sg = ΔnoutSizeGraph(vin)
    if vout in keys(sg.metaindex[:vertex])
        add_edge!(sg, sg[vin, :vertex], sg[vout, :vertex])
        is_cyclic(ΔnoutSizeGraph(vin)) && return prealignsizes(CheckAligned(s.ifnok), vin, vout, will_rm)
    end
    return prealignsizes(s.ifok, vin, vout, will_rm)
end

function prealignsizes(s::ChangeNinOfOutputs, vin, vout, will_rm)
    expected = nout(vin) + s.Δoutsize
    Δsize(ΔNout{Exact}(vin, s.Δoutsize, ΔSizeFailNoOp()), all_in_Δsize_graph(vin, Output()))
    return nout(vin) == expected
end

function prealignsizes(s::ApplyMutation, vin, vout, will_rm)
    if prealignsizes(s.strategy, vin, vout, will_rm)
        apply_mutation.(all_in_graph(vin))
        return true
    end
    return false
end

prealignsizes(s::PostApplyMutation, vin, vout, will_rm) = prealignsizes(s.strategy, vin, vout, will_rm)

function prealignsizes(s::SelectOutputs, vin, vout, will_rm)
    if prealignsizes(s.alignstrategy, vin, vout, will_rm)
        return Δoutputs(s.selectstrategy, vin, s.valuefun)
    end
    return false
end

prealignsizes(s::PostSelectOutputs, vin, vout, will_rm) = prealignsizes(s.alignstrategy, vin, vout, will_rm)

function prealignsizes(s::Union{IncreaseSmaller, DecreaseBigger}, vin, vout, will_rm)
    Δinsize = nout(vin) - tot_nin(vout)
    Δoutsize = -Δinsize

    strat = proceedwith(s, Δinsize) ? Δninstrat(vout, Δinsize) : Δnoutstrat(vin, Δoutsize)
    success = prealignsizes(strat, vin, vout)

    if !success
        return prealignsizes(s.fallback, vin, vout, will_rm)
    end
    return true
end
proceedwith(::DecreaseBigger, Δ::Integer) = Δ <= 0
proceedwith(::IncreaseSmaller, Δ::Integer) = Δ >= 0

Δninstrat(v, Δ) = Δninstrat(trait(v), v, Δ)
Δninstrat(t::DecoratingTrait, v, Δ) = Δninstrat(base(t), v, Δ)
Δninstrat(::Immutable, v, Δ) = ΔSizeFailNoOp()
Δninstrat(::MutationSizeTrait, v, Δ) = ΔNin{Exact}(v, [Δ], ΔSizeFailNoOp())
Δninstrat(::SizeTransparent, v, Δ) = ΔNout{Exact}(v, Δ, ΔSizeFailNoOp())

Δnoutstrat(v, Δ) = Δnoutstrat(trait(v), v, Δ)
Δnoutstrat(t::DecoratingTrait, v, Δ) = Δnoutstrat(base(t), v, Δ)
Δnoutstrat(::Immutable, v, Δ) = ΔSizeFailNoOp()
Δnoutstrat(::MutationSizeTrait, v, Δ) = ΔNout{Exact}(v, Δ, ΔSizeFailNoOp())

function prealignsizes(s::AlignSizeBoth, vin, vout, will_rm)

    strat = AlignNinToNoutVertices(vin, vout, 1:length(nin(vout)), AlignNinToNout(), ΔSizeFailNoOp())
    success = prealignsizes(strat, vin, vout)

    if !success
        return prealignsizes(s.fallback, vin, vout, will_rm)
    end
    return true
end

function prealignsizes(s::AbstractΔSizeStrategy, vin, vout)
    vin_all = all_in_Δsize_graph(vin, Output())
    vout_all = all_in_Δsize_graph(vout, Input())

    verts = union(vin_all, vout_all)
    success, nins, nouts = newsizes(s, verts)

    if success
        Δsize(nins, nouts, verts)
        return true
    end
    return false
end

# Boilerplate
postalignsizes(s::AbstractAlignSizeStrategy, v) = postalignsizes(s, v, v, missing)
postalignsizes(s::AbstractAlignSizeStrategy, vin, vout, pos) = true

# Failure cases
postalignsizes(::FailAlignSizeError, vin, vout, pos) = error("Could not align sizes of $(vin) and $(vout)!")
function postalignsizes(s::FailAlignSizeWarn, vin, vout, pos)
    @warn s.msgfun(vin, vout)
    return postalignsizes(s.andthen, vin, vout, pos)
end
function postalignsizes(s::FailAlignSizeRevert, vin, vout, pos)
    n = sum(inputs(vout) .== vin)
    # Can maybe be supported by comparing nin_org(vout) to nout_org(vin): If they match then vin was there before, else it shall be removed?
    @assert n <= 1 "Case when vin is input to vout multiple times not implemented!"

    if n == 1
        remove_edge!(vin, vout, strategy=NoSizeChange())
    else #if n == 0, but n > 1 not implemented
        create_edge!(vin, vout, pos=pos, strategy=NoSizeChange())
        reset_in!(op(vout)) # create_edge! sets nin_org to 0 for new edges
    end
    return false
end

function postalignsizes(s::PostSelectOutputs, vin, vout, pos)
    if postalignsizes(s.alignstrategy, vin, vout, pos)
        vin_all = all_in_Δsize_graph(vin, Output())
        vout_all = all_in_Δsize_graph(vout, Input())

        verts = union(vin_all, vout_all)
        isempty(verts) && return true

        success = Δoutputs(s.selectstrategy, verts, s.valuefun)

        if !success
            return postalignsizes(s.fallback, vin, vout, pos)
        end
        return success
    end
    return false
end

function postalignsizes(s::PostApplyMutation, vin, vout, pos)
    if postalignsizes(s.strategy, vin, vout, pos)
        apply_mutation.(all_in_graph(vin))
        return true
    end
    return false
end

function postalignsizes(s::CheckCreateEdgeNoSizeCycle, vin, vout, pos)
    sg = ΔnoutSizeGraph(vin)
    is_cyclic(ΔnoutSizeGraph(vin)) && return postalignsizes(s.ifnok, vin, vout, pos)
    return postalignsizes(s.ifok, vin, vout, pos)
end

# Ok, this one actually does something...
function postalignsizes(s::PostAlignJuMP, vin, vout, pos)
    vin_all = all_in_Δsize_graph(vin, Output())
    vout_all = all_in_Δsize_graph(vout, Input())

    verts = union(vin_all, vout_all)
    success, nins, nouts = newsizes(s.sizestrat, verts)
    if !success
        return postalignsizes(s.fallback, vin, vout, pos)
    end
    Δsize(nins, nouts, verts)
    return success
end

"""
    insert!(vin::AbstractVertex, factory::Function, outselect::Function=identity)

Replace `vin` as input to all outputs of `vin` with vertex produced by `factory`.

The function `factory` takes `vin` as input.

Example:

```text
Before:

    vin
   / | \\
  /  |  \\
 v₁ ... vₙ

After:

    vin
     |
   vnew₁
     ⋮
   vnewₖ
   / | \\
  /  |  \\
 v₁ ... vₙ
```
Note that the connection `vin` -> `vnew₁` as well as all connections `vnewₚ -> vnewₚ₊₁` is done by `factory` in the above example.

The function `outselect` can be used to select a subset of outputs to replace (default all).

"""
function Base.insert!(vin::AbstractVertex, factory::Function, outselect::Function=identity)

    prevouts = copy(outselect(outputs(vin)))
    deleteat!(outputs(vin), findall(vx -> vx in prevouts, outputs(vin)))

    vnew = factory(vin)
    for vout in prevouts
        inds = findall(vx -> vx == vin, inputs(vout))
        inputs(vout)[inds] .= vnew
        push!(outputs(vnew), vout)
    end
end

"""
    create_edge!(from::AbstractVertex, to::AbstractVertex; [pos], [strategy])

Create and edge from `from` to `to` at index `pos` in `inputs(to)`.

Sizes will be adjusted based on given `strategy`.
"""
function create_edge!(from::AbstractVertex, to::AbstractVertex; pos = length(inputs(to))+1, strategy = default_create_edge_strat(to))

    prealignsizes(strategy, from, to, v -> false) || return

    push!(outputs(from), to) # Order should never matter for outputs
    insert!(inputs(to), pos, from)

    add_input!(op(to), pos, nout(from))
    add_output!(op(to), trait(to))

    postalignsizes(strategy, from, to, pos)
end

default_create_edge_strat(v::AbstractVertex) = default_create_edge_strat(trait(v),v)
default_create_edge_strat(t::DecoratingTrait,v) = default_create_edge_strat(base(t),v)
default_create_edge_strat(::SizeStack,v) = CheckCreateEdgeNoSizeCycle(ifok=PostAlignJuMP())
default_create_edge_strat(::SizeInvariant,v) = CheckCreateEdgeNoSizeCycle(ifok=IncreaseSmaller())
default_create_edge_strat(::SizeAbsorb,v) = NoSizeChange()

function add_input!(::MutationOp, pos, size) end
add_input!(s::IoSize, pos, size) = insert!(s.nin, pos, size)
add_input!(s::IoIndices, pos, size) = insert!(s.in, pos, collect(1:size+1))
function add_input!(s::IoChange, pos, size)
    add_input!(s.size, pos, 0)
    add_input!(s.indices, pos, 0)
    insert!(s.inΔ, pos, size)
end

function add_output!(::MutationOp, t::MutationTrait) end
add_output!(s::MutationOp, t::DecoratingTrait) = add_output!(s, base(t))
add_output!(s::IoSize, ::SizeStack) = Δnout(s, sum(nin(s)) - nout(s))
add_output!(s::IoIndices, ::SizeStack) = Δnout(s, vcat(s.out, (length(s.out):length(s.out)+(sum(nin(s)) - nout(s)))))
add_output!(s::IoChange, t::SizeStack) = Δnout(s, sum(nin(s)) - nout(s))

"""
    remove_edge!(from::AbstractVertex, to::AbstractVertex; [nr], [strategy])

Remove edge from `from` to `to`.

If there are multiple edges from `from` to `to` then `nr` can be used to distinguish which one shall be removed.

Sizes will be adjusted based on given `strategy`.
"""
function remove_edge!(from::AbstractVertex, to::AbstractVertex; nr = 1, strategy = default_remove_edge_strat(to))

    prealignsizes(strategy, from, to, v -> false) || return

    in_ind  = findall(vx -> vx == from, inputs(to))[nr]
    out_ind = findall(vx -> vx == to, outputs(from))[nr]
    deleteat!(inputs(to), in_ind)
    deleteat!(outputs(from), out_ind)

    rem_input!(op(to), in_ind...)
    add_output!(op(to), trait(to))

    postalignsizes(strategy, from, to, in_ind)
end

default_remove_edge_strat(v::AbstractVertex) = default_remove_edge_strat(trait(v),v)
default_remove_edge_strat(t::DecoratingTrait,v) = default_remove_edge_strat(base(t),v)
default_remove_edge_strat(::SizeStack,v) = PostAlignJuMP()
default_remove_edge_strat(::SizeInvariant,v) = NoSizeChange()
default_remove_edge_strat(::SizeAbsorb,v) = NoSizeChange()

function rem_input!(::MutationOp, pos...) end
rem_input!(s::IoSize, pos...) = deleteat!(s.nin, collect(pos))
rem_input!(s::IoIndices, pos...) = deleteat!(s.in, collect(pos))
function rem_input!(s::IoChange, pos...)
    rem_input!(s.size, pos...)
    rem_input!(s.indices, pos...)
    deleteat!(s.inΔ, collect(pos))
end


text:
@def F_YR begin #(V, U, Ax, r, SA) # rear tire force

	if length(Ax) != 1
		Fyr = zeros(Float64, (length(Ax),1))
		for ii in eachindex(Ax)
			Fyr[ii] = ((PD2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)) + ((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))^2+EP^2)^(0.5))*0.001)+EP))*((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)) - ((PE2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PE1)*(1 - PE3*(((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1))/((((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)^2 + EP^2)^(0.5))))*((((PK1*sin(2*atan(PK2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)) + ((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))^2+EP^2)^(0.5))*0.001)+EP))*((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)) + ((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)))/(((PD2*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff))^2+EP^2)^(0.5))*0.001)+EP))*((atan((V[ii] - lb*R[ii])/(U[ii]+EP))) + PH2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff) + PH1)))))) + (PV2*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)^2 + PV1*(FzR0 + (Ax[ii] - V[ii]*R[ii])*dFzx_coeff)));
		end
	else
		Fyr = Float64
		Fyr = ((PD2*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*(FzR0 + (Ax - V*R)*dFzx_coeff))*sin(PC1*atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax - V*R)*dFzx_coeff))))/(((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)) + ((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)))/(((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff))^2+EP^2)^(0.5))*0.001)+EP))*((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*dFzx_coeff) + PH1)) - ((PE2*(FzR0 + (Ax - V*R)*dFzx_coeff) + PE1)*(1 - PE3*(((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*dFzx_coeff) + PH1))/((((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*dFzx_coeff) + PH1)^2 + EP^2)^(0.5))))*((((PK1*sin(2*atan(PK2*(FzR0 + (Ax - V*R)*dFzx_coeff))))/(((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)) + ((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)))/(((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff))^2+EP^2)^(0.5))*0.001)+EP))*((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*dFzx_coeff) + PH1)) - atan((((PK1*sin(2*atan(PK2*(FzR0 + (Ax - V*R)*dFzx_coeff))))/(((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)) + ((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)))/(((PD2*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PD1*PC1*(FzR0 + (Ax - V*R)*dFzx_coeff))^2+EP^2)^(0.5))*0.001)+EP))*((atan((V - lb*R)/(U+EP))) + PH2*(FzR0 + (Ax - V*R)*dFzx_coeff) + PH1)))))) + (PV2*(FzR0 + (Ax - V*R)*dFzx_coeff)^2 + PV1*(FzR0 + (Ax - V*R)*dFzx_coeff)));
	end
	Fyr
end


text:
"""
    _A_method(parameters::Parameters1PL, latent::Latent1D)
"""
function _A_method(parameters::Parameters1PL, latent::Latent1D)
    return _item_expected_information(parameters, latent)
end

"""
    _A_method(parameters::Parameters2PL, latent::Latent1D)
"""
function _A_method(parameters::Parameters2PL, latent::Latent1D)
    return LinearAlgebra.tr(_item_expected_information(parameters, latent))
end

"""
    A_method(item::AbstractItem, examinee::AbstractExaminee)

# Description
Computes the trace of the expected information matrix for an item and an examinee.

# Arguments
- **`item::AbstractItem`**: The item.
- **`examinee::AbstractExaminee`**: The examinee at which computing the information.

# Output
It returns a `Float64` scalar.
"""
function A_method(item::AbstractItem, examinee::AbstractExaminee)
    return _A_method(item.parameters, examinee.latent)
end


text:
"""
ResTrack(args...)  -> PyPlot Plot

Return a plot of resistivities logs in a single track.

 The next table show the list of variables allowed:
|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Depth|Mandatory| -- | Array{Number,1}| 1D Array of Depths|
|Res|Mandatory| -- | Array{Number,n}|2D Array of Resistivities Values sorted by colums|
|DepthFrom|Optional| minimum(Depth) | Number |Set the Depth upper limit to plot|
|DepthTo|Optional| maximum(Depth) | Number |  Set the Depth lower limit to plot|
|c1|Optional|c1=RGB(0.0,0.0,1.0)|c1=RGB(r,g,b)|Set the initial color for
 a Color gradient Resistivities curves|
|c1|Optional|c1=RGB(0.1,0.6,0.1)|c1=RGB(r,g,b)|Set the final color for
 a Color gradient Resistivities curves|
|Dtick|Optional|true|Dtick=Bool|Set the Depth ticks On/Off|

"""
function ResTrack(Depth;Res,DepthFrom=minimum(Depth),DepthTo=maximum(Depth),
        c1=RGB(0.0,0.0,1.0),c2=RGB(0.1,0.6,0.1),Dtick=false,ResRange=[0.2,20000])

n=size(Res,2)
pcolor=range(c1, stop=c2, length=n)
pcols = map(col -> (red(col), green(col), blue(col)), pcolor)
linestyles="-"
linewidths=collect(range(1, stop=1.5, length=n))
rsax=gca()
for i=1:n
    rsax.plot(Res[:,i],Depth,color=pcols[i],linestyle=linestyles,linewidth=linewidths[i])
    end
rsax.set_ylim([DepthTo,DepthFrom])
setp(rsax.get_yticklabels(),visible=Dtick)
rsax.set_xscale("log")
xlim(ResRange[1],ResRange[2])
ticks=round.(10 .^(range(log10(ResRange[1]), stop=log10(ResRange[2]), length=1+Int(log10(ResRange[2]/ResRange[1])))),digits=1)
rsax.set_xticks(ticks)
rsax.set_xticklabels(ticks)
xlabel("Resistivity")
rsax.grid(true)
rsax.grid(true, which="minor", axis="x",linewidth=0.5, linestyle="--")
rsax.xaxis.tick_top()
rsax.xaxis.set_label_position("top")
rsax.tick_params("both",labelsize=6)
rsax.xformatter = :auto
    rsax.set_yticks(range(DepthFrom, stop=DepthTo, length=51),minor=true)
    rsax.set_yticks(range(DepthFrom, stop=DepthTo, length=11))

    end


text:
# two-parameter tuning plot (with nested parameters)

using MLJ
using DataFrames, Statistics

Xraw = rand(300,3)
y = exp(Xraw[:,1] - Xraw[:,2] - 2Xraw[:,3] + 0.1*rand(300))
X = DataFrame(Xraw)

train, test = partition(eachindex(y), 0.70); # 70:30 split

knn_model = KNNRegressor(K=10)
knn = machine(knn_model, X, y)

ensemble_model = EnsembleModel(atom=knn_model, n=20)
ensemble = machine(ensemble_model, X, y)

B_range = range(ensemble_model, :bagging_fraction, lower= 0.5, upper=1.0, scale = :linear)
K_range = range(knn_model, :K, lower=1, upper=100, scale=:log10)
nested_ranges = (atom = (K = K_range,), bagging_fraction = B_range)

tuning = Grid(resolution=12)
resampling = Holdout(fraction_train=0.8)

tuned_ensemble_model = TunedModel(model=ensemble_model,
    tuning=tuning, resampling=resampling, nested_ranges=nested_ranges)

tuned_ensemble = machine(tuned_ensemble_model, X[train,:], y[train])
fit!(tuned_ensemble);

# data needed for plotting is here:
tuned_ensemble.report

# uncomment 4 lines to see plots
# using Plots
# pyplot()
# plot(tuned_ensemble)
# heatmap(tuned_ensemble)

true


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

if !isdefined(@__MODULE__, :DefsSubteamId)
    const DefsSubteamId = String
else
    @warn("Skipping redefinition of DefsSubteamId to String")
end


text:
const Diencephalon = BilateralStructure("diencephalon")

const Epithalamus = BilateralStructure("epithalamus")

const PinealBodyGland = BilateralStructure("pineal body gland")

const HabenularNuclei = BilateralStructure("habenular nuclei")

const StriaMedullaris = BilateralStructure("stria medullaris")

const TaeniaThalami = BilateralStructure("taenia thalami")

"""
    Thalamus
"""
const Thalamus = BilateralStructure("thalamus")
const ThalamusProper = Thalamus

"""
    AnteriorNuclearGroup
"""
const AnteriorNuclearGroup = BilateralStructure("anterior nuclear group")
const AnteroventralNucleus = BilateralStructure("anteroventral nucleus")
const AnterodorsalNucleus = BilateralStructure("anterodorsal nucleus")
const AnteromedialNucleus = BilateralStructure("anteromedial nucleus")

const MedialDorsalNucleus = BilateralStructure("medial dorsal nucleus")

const MedialNuclearGroup = BilateralStructure("medial nuclear group")
const ParvocellularPartOfThalamus = BilateralStructure("parvocellular part of thalamus")
const MagnocellularPartOfThalamus = BilateralStructure("magnocellular part of thalamus")

const MidlineNuclearGroup = BilateralStructure("midline nuclear group")
const ParatenialNucleus = BilateralStructure("paratenial nucleus")
const ParaventricularNucleusOfThalamus = BilateralStructure("pareventricular nucleus of thalamus")
const ReuniensNucleus = BilateralStructure("reuniens nucleusi")
const RhomboidalNucleus = BilateralStructure("rhomboidal nucleus")

const IntralaminarNuclearGroup = BilateralStructure("intralaminar nuclear group")

const CentromedianNucleus = BilateralStructure("centromedian nucleus")

const ParafascicularNucleus = BilateralStructure("parafascicular nucleus")

const ParacentralNucleus = BilateralStructure("paracentral nucleus")

const CentralLateralNucleus = BilateralStructure("central lateral nucleus")

const LateralNuclearGroup = BilateralStructure("lateral nuclear group")
const LateralDorsalNucleus = BilateralStructure("lateral dorsal nucleus")
const LateralPosteriorNucleus = BilateralStructure("lateral posterior nucleus")
const Pulvinar = BilateralStructure("pulvinar")

"""
    VentralNuclearGroup
"""
const VentralNuclearGroup = BilateralStructure("ventral nuclear group")
const VentralAnteriorNucleus = BilateralStructure("ventral anterior nucleus")
const VentralLateralNucleus = BilateralStructure("ventral lateral nucleus")
const VentralPosteriorLateralNucleus = BilateralStructure("ventral posterior lateral nucleus")
const VentralPosteriorMedialNucleus = BilateralStructure("ventral posterior medial nucleus")
const VentralPosteriorNucleus = BilateralStructure("ventral posterior nucleus")

"""
    Metathalamus
"""
const Metathalamus = BilateralStructure("metathalamus")
const MedialGeniculateBody = BilateralStructure("medial geniculate body")
const LateralGeniculateBody = BilateralStructure("lateral geniculate body")

const ThalamicReticularNucleus = BilateralStructure("thalamic reticular nucleus")


text:
struct Thing2
  other_thing::Thing1
end


text:
module Common

using Base.Iterators
using Unicode

export urinormalize

"""
    urinormalize(s::AbstractString)

Normalize `s` in a way expected to be suitable for URIs.
"""
function urinormalize(s::AbstractString)::String
    t = Unicode.normalize(s; compat=true, stripcc=true, rejectna=true,
                             stripignore=true, casefold=true)
    t = replace(t, r"\s" => "-")
    t = filter(c -> isletter(c) || isnumeric(c) || c == '-', t)
    if isempty(t)
        throw(ArgumentError(
            "Argument $(repr(s)) cannot be meaningfully URI-normalized."))
    end
    join(collect(take(t, 30)))
end

end


text:
0.0015020215673811625
0.00017680822734654275
0.00048637064841981195
-0.0006597479269052878
-0.019413611173469538
-0.042048427967227915
0.01708748187542781
0.0608144136683526
-0.06968149360629872
-0.09505504623756987
0.016626151803912966
-0.06490767715200972
0.04853521707350193
0.3123471343784102
-0.020250492628022258
-0.13315364999147442
0.2894887276848721
0.11039229876148088
-0.02574979051592033
-0.04543847658479455
-0.093534171223262
-0.07040279950792366
0.3687247182251427
-0.16206760719005292
-0.03730254719682338
-0.07141684250723922
-0.053982572774078884
-0.22491093718238434
-0.19093489443858308
-0.12230409564769004
0.03332591659058223
0.11451200993031567
-0.028250597885822647
0.17563864881873825
-0.08613987990076162
-0.07344723570235402
-0.07913030260496898
0.1420463815893184
-0.05787227576797076
0.16855632503151965
-0.007762606626547975
0.05442187707509019
0.10654738017237805
0.04705537808367204
-0.10828112731084161
-0.05286746608265549
-0.08487119065889706
0.04250851760058027
0.14287508117261785
0.29004243904501104
-0.0634003268210348
-0.35896795566489664

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
using BatchedRoutines, Test, CuArrays

@testset "Testing batched_scal! with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    A = rand(elty, 100)
    B = rand(elty, 10, 10, 100)

    dA, dB = CuArray(A), CuArray(B)
    batched_scal!(dA, dB)

    for k in 1:100
        B[:, :, k] .= A[k] * B[:, :, k]
    end

    test_B = Array(dB)

    @test B ≈ test_B
end

@testset "Testing batched_scal with $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    A = rand(elty, 100)
    B = rand(elty, 10, 10, 100)

    dA, dB = CuArray(A), CuArray(B)
    test_B = Array(batched_scal(A, B))

    for k in 1:100
        B[:, :, k] .= A[k] * B[:, :, k]
    end

    @test B ≈ test_B
end


text:
function xunderlying(x)
  2x
end


text:
#
# tests of test stuff!, exp .ExampleHelpers
#
using Test
using ScottishTaxBenefitModel
using .ModelHousehold
using .Definitions
using .ExampleHelpers

@testset "Adding/Removing people tests" begin

    hh = get_example( single_hh )
    head = get_head( hh )
    @test num_people( hh ) == 1
    @test num_adults( hh ) == 1
    @test num_children( hh ) == 0
    @test get_spouse( hh ) === nothing

    add_spouse!( hh, 30, Female )
    @test num_people( hh ) == 2
    @test num_adults( hh ) == 2
    @test num_children( hh ) == 0
    @test get_spouse( hh ) !== nothing

    spouse = get_spouse( hh )
    @test spouse.age == 30
    @test spouse.sex == Female

    bus = get_benefit_units( hh )
    @test size(bus)[1] == 1
    bu = bus[1]
    @test num_people( bu ) == 2
    @test num_adults( bu ) == 2
    @test num_children( bu ) == 0
    spb = get_spouse( bu )
    @test spb !== nothing

    @test spb.sex == spouse.sex
    @test spb.age == spouse.age
    @test spb.relationships[head.pid] == Spouse
    @test head.relationships[spouse.pid] == Spouse

    add_child!( hh, 3, Female )
    bus = get_benefit_units( hh )
    @test size(bus)[1] == 1
    bu = bus[1]
    @test num_people( bu ) == 3
    @test num_adults( bu ) == 2
    @test num_children( bu ) == 1
    for pid in bu.children
        @test spb.relationships[pid] == Parent
        @test head.relationships[pid] == Parent
    end
    ch1 = bu.people[bu.children[1]]
    @test ch1.relationships[spb.pid] == Son_or_daughter_incl_adopted
    @test ch1.relationships[head.pid] == Son_or_daughter_incl_adopted
    @test ch1.age == 3
    @test ch1.sex == Female

    add_child!( hh, 5, Male )
    bus = get_benefit_units( hh )
    bu = bus[1]
    @test num_people( bu ) == 4
    @test num_adults( bu ) == 2
    @test num_children( bu ) == 2

    # ch1 is ELDEST in bus, always..
    ch1 = bu.people[bu.children[1]]
    @test ch1.relationships[spb.pid] == Son_or_daughter_incl_adopted
    @test ch1.relationships[head.pid] == Son_or_daughter_incl_adopted
    @test ch1.age == 5
    @test ch1.sex == Male
    for pid in bu.children
        @test spb.relationships[pid] == Parent
        @test head.relationships[pid] == Parent
    end
    @test size(bu.children)[1] == 2

    ## TODO TEST REMOVING
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using MPI
using Test
using LinearAlgebra
using Random
using KernelAbstractions, StaticArrays

using ClimateMachine
using ClimateMachine.LinearSolvers
using ClimateMachine.ColumnwiseLUSolver:
    band_lu_kernel!, band_forward_kernel!, band_back_kernel!

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const device = ArrayType == Array ? CPU() : CUDA()

function band_to_full(B, p, q)
    _, n = size(B)

    A = similar(B, n, n) # assume square
    fill!(A, 0)
    for j in 1:n, i in max(1, j - q):min(j + p, n)
        A[i, j] = B[q + i - j + 1, j]
    end
    A
end

let
    Nq = 2
    nstate = 3
    nvertelem = 5
    nhorzelem = 4
    eband = 2

    FT = Float64
    m = n = Nq * nstate * nvertelem
    p = q = Nq * nstate * eband - 1

    Random.seed!(1234)
    AB = rand(FT, Nq, Nq, p + q + 1, n, nhorzelem)
    AB[:, :, q + 1, :, :] .+= 10 # Make A's diagonally dominate

    Random.seed!(5678)
    b = rand(FT, Nq, Nq, Nq, nstate, nvertelem, nhorzelem)
    x = similar(b)

    perm = (4, 3, 5, 1, 2, 6)
    bp = reshape(PermutedDimsArray(b, perm), n, Nq, Nq, nhorzelem)
    xp = reshape(PermutedDimsArray(x, perm), n, Nq, Nq, nhorzelem)

    d_F = ArrayType(AB)

    groupsize = (Nq, Nq)
    ndrange = (Nq, Nq, nhorzelem)

    event = Event(device)
    event = band_lu_kernel!(device, groupsize, ndrange)(
        d_F,
        Val(Nq),
        Val(Nq),
        Val(Nq),
        Val(nstate),
        Val(nvertelem),
        Val(nhorzelem),
        Val(eband),
        dependencies = (event,),
    )
    wait(device, event)

    F = Array(d_F)

    for h in 1:nhorzelem, j in 1:Nq, i in 1:Nq
        B = AB[i, j, :, :, h]
        G = band_to_full(B, p, q)
        GLU = lu!(G, Val(false))

        H = band_to_full(F[i, j, :, :, h], p, q)

        @test H ≈ G

        xp[:, i, j, h] .= GLU \ bp[:, i, j, h]
    end

    b = reshape(b, Nq * Nq * Nq, nstate, nvertelem * nhorzelem)
    x = reshape(x, Nq * Nq * Nq, nstate, nvertelem * nhorzelem)

    d_x = ArrayType(b)

    event = Event(device)
    event = band_forward_kernel!(device, groupsize, ndrange)(
        d_x,
        d_F,
        Val(Nq),
        Val(Nq),
        Val(nstate),
        Val(nvertelem),
        Val(nhorzelem),
        Val(eband),
        dependencies = (event,),
    )

    event = band_back_kernel!(device, groupsize, ndrange)(
        d_x,
        d_F,
        Val(Nq),
        Val(Nq),
        Val(nstate),
        Val(nvertelem),
        Val(nhorzelem),
        Val(eband),
        dependencies = (event,),
    )
    wait(device, event)

    @test x ≈ Array(d_x)
end


text:
export iauIr
"""
Initialize an r-matrix to the identity matrix.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Returned:
   r       double[3][3]    r-matrix

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauIr(double r[3][3])
function iauIr()
   # Allocate return values
   I  = zeros(Float64, 3, 3)

   ccall((:iauIr, libsofa_c), Cvoid,
        (Ptr{Cdouble},),
        I)

   return SMatrix{3,3}(I')
end


text:
using NCDatasets

sz = (40,10)
filename = tempname()
#filename = "/tmp/test-7.nc"
# The mode "c" stands for creating a new file (clobber)

NCDataset(filename,"c") do ds

    # define the dimension "lon" and "lat"
    defDim(ds,"lon",sz[1])
    defDim(ds,"lat",sz[2])

    # variables
    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64]
    #for T in [Float32]
        local data
        data = fill(T(123),sz)

        v = defVar(ds,"var-$T",T,("lon","lat");
                              shuffle = true,
                              chunksizes = [20,5],
                              deflatelevel = 9,
                              checksum = :nochecksum
                              )
        # check checksum method
        checksummethod = checksum(v)
        @test checksummethod == :nochecksum

        # change checksum method
        checksum(v,:fletcher32)
        checksummethod = checksum(v)
        @test checksummethod == :fletcher32

        # check chunking
        storage,chunksizes = chunking(v)
        @test storage == :chunked
        @test chunksizes[1] == 20

        # change chunking
        chunking(v,:chunked,[3,3])
        storage,chunksizes = chunking(v)
        @test storage == :chunked
        #@show chunksizes
        @test chunksizes[1] == 3

        # check compression
        isshuffled,isdeflated,deflate_level = deflate(v)
        @test isshuffled == true
        @test isdeflated == true
        @test deflate_level == 9

        # change compression
        deflate(v,false,true,4)
        isshuffled,isdeflated,deflate_level = deflate(v)
        @test isshuffled == false
        @test isdeflated == true
        @test_broken deflate_level == 4

        # write an array
        v[:,:] = data
        @test all(v[:,:] .== data)

        v = defVar(ds,"var2-$T",T,("lon","lat");
                              shuffle = true,
                              chunksizes = [20,5],
                              deflatelevel = 9,
                              checksum = :fletcher32
                              )
        checksummethod = checksum(v)
        @test checksummethod == :fletcher32

    end
end


text:
"""
    SimpleRandomWalkEnv

A simple random walk environment for tutorial.

# Example
```julia
julia> env = SimpleRandomWalkEnv()
SimpleRandomWalkEnv(7, 3, 3, [-1, 1])

julia> env(1)
(observation = 2, reward = 0.0, isdone = false)

julia> reset!(env)
(observation = 3, isdone = false)
```
"""
mutable struct SimpleRandomWalkEnv <: AbstractSyncEnvironment{DiscreteSpace,DiscreteSpace,1}
    N::Int
    start::Int
    state::Int
    actions::Vector{Int}
    SimpleRandomWalkEnv(;N=7, actions=[-1, 1], start=3) = new(N, start, start, actions)
end

function (env::SimpleRandomWalkEnv)(a::Int)
    env.state = min(max(env.state + env.actions[a], 1), env.N)
    (observation = env.state,
     reward      = env.state == env.N ? 3.0 : (env.state == 1 ? 1. : 0.),
     isdone      = env.state == env.N || env.state == 1)
end

function reset!(env::SimpleRandomWalkEnv)
    env.state = env.start
    (observation = env.state, isdone = false)
end

observe(env::SimpleRandomWalkEnv) = (observation = env.state, isdone = env.state == env.N || env.state == 1)
observationspace(env::SimpleRandomWalkEnv) = DiscreteSpace(env.N)
actionspace(env::SimpleRandomWalkEnv) = DiscreteSpace(length(env.actions))


text:
using HamiltonianABC
using Base.Test
using Parameters
using Distributions
using Cubature
using StatsBase

import HamiltonianABC: logdensity, simulate!

# use consistent random seed for unit tests
srand(UInt32[0x8c736cc0, 0x63556b2b, 0x808f767c, 0xc912a668])

include("utilities.jl")
include("test-mcmc.jl")

# toy model: exponential through normal
include("test-toyexponential.jl")

# toy model: G-and-K through normal mixture
include("EM_GaussianMixtureModel.jl")
include("g_and_k_quantile.jl")


text:
# Allow this to be set manually
function rasterpath()
    if haskey(ENV, "RASTERDATASOURCES_PATH") && isdir(ENV["RASTERDATASOURCES_PATH"])
        ENV["RASTERDATASOURCES_PATH"]
    else
        error("You must set `ENV[\"RASTERDATASOURCES_PATH\"]` to a path in your system")
    end
end

function remove_rasters()
    # May need an "are you sure"? - this could be a lot of GB of data to lose
    ispath(rasterpath()) && rm(rasterpath())
end

function remove_rasters(T::Type)
    ispath(rasterpath(T)) && rm(rasterpath(T))
end

function remove_rasters(::Type{TS}, ::Type{TD}) where {TS <: RasterDataSource, TD <: RasterDataSet}
    ispath(_raster_assets_folder(TS, TD)) && rm(_raster_assets_folder(TS, TD); recursive=false)
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# SolveBalances: Solves model equations from TSTART to TSTOP given parameters in data_dictionary.
# Type: GRN-JULIA
# Version: 1.0
#
# Input arguments:
# TSTART  - Time start
# TSTOP  - Time stop
# Ts - Time step
# data_dictionary  - Data dictionary instance (holds model parameters)
#
# Return arguments:
# TSIM - Simulation time vector
# X - Simulation state array (NTIME x NSPECIES)
# ----------------------------------------------------------------------------------- #
function SolveBalances(TSTART,TSTOP,Ts,data_dictionary)

  # Get required stuff from DataFile struct -
  TSIM = collect(TSTART:Ts:TSTOP);
  initial_condition_vector = data_dictionary["initial_condition_array"];

  # Call the ODE solver -
  fbalances(t,y) = Balances(t,y,data_dictionary);
  (t,y) = ode23s(fbalances,initial_condition_vector,TSIM;points=:specified);

  # Map -
  number_of_timesteps = length(t)
  number_of_states = length(initial_condition_vector)
  X = zeros(number_of_timesteps,number_of_states)
  for state_index = 1:number_of_states
    tmp = map(y->y[state_index],y)
    for time_index = 1:number_of_timesteps
      X[time_index,state_index] = tmp[time_index]
    end
  end

  # Check for negatives -
  idx_n = find(X.<0)
  X[idx_n] = 0.0

  # return time and state -
  return (t,X);
end