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

const IEEEFloat = Union{Float16, Float32, Float64}

## floating point traits ##

"""
    Inf16

Positive infinity of type [`Float16`](@ref).
"""
const Inf16 = bitcast(Float16, 0x7c00)
"""
    NaN16

A not-a-number value of type [`Float16`](@ref).
"""
const NaN16 = bitcast(Float16, 0x7e00)
"""
    Inf32

Positive infinity of type [`Float32`](@ref).
"""
const Inf32 = bitcast(Float32, 0x7f800000)
"""
    NaN32

A not-a-number value of type [`Float32`](@ref).
"""
const NaN32 = bitcast(Float32, 0x7fc00000)
const Inf64 = bitcast(Float64, 0x7ff0000000000000)
const NaN64 = bitcast(Float64, 0x7ff8000000000000)

const Inf = Inf64
"""
    Inf, Inf64

Positive infinity of type [`Float64`](@ref).
"""
Inf, Inf64

const NaN = NaN64
"""
    NaN, NaN64

A not-a-number value of type [`Float64`](@ref).
"""
NaN, NaN64

## conversions to floating-point ##
Float16(x::Integer) = convert(Float16, convert(Float32, x))
for t in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
    @eval promote_rule(::Type{Float16}, ::Type{$t}) = Float16
end
promote_rule(::Type{Float16}, ::Type{Bool}) = Float16

for t1 in (Float32, Float64)
    for st in (Int8, Int16, Int32, Int64)
        @eval begin
            (::Type{$t1})(x::($st)) = sitofp($t1, x)
            promote_rule(::Type{$t1}, ::Type{$st}) = $t1
        end
    end
    for ut in (Bool, UInt8, UInt16, UInt32, UInt64)
        @eval begin
            (::Type{$t1})(x::($ut)) = uitofp($t1, x)
            promote_rule(::Type{$t1}, ::Type{$ut}) = $t1
        end
    end
end
(::Type{T})(x::Float16) where {T<:Integer} = T(Float32(x))

Bool(x::Real) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

promote_rule(::Type{Float64}, ::Type{UInt128}) = Float64
promote_rule(::Type{Float64}, ::Type{Int128}) = Float64
promote_rule(::Type{Float32}, ::Type{UInt128}) = Float32
promote_rule(::Type{Float32}, ::Type{Int128}) = Float32

function Float64(x::UInt128)
    x == 0 && return 0.0
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 53
        y = ((x % UInt64) << (53-n)) & 0x000f_ffff_ffff_ffff
    else
        y = ((x >> (n-54)) % UInt64) & 0x001f_ffff_ffff_ffff # keep 1 extra bit
        y = (y+1)>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt64(trailing_zeros(x) == (n-54)) # fix last bit to round to even
    end
    d = ((n+1022) % UInt64) << 52
    reinterpret(Float64, d + y)
end

function Float64(x::Int128)
    x == 0 && return 0.0
    s = ((x >>> 64) % UInt64) & 0x8000_0000_0000_0000 # sign bit
    x = abs(x) % UInt128
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 53
        y = ((x % UInt64) << (53-n)) & 0x000f_ffff_ffff_ffff
    else
        y = ((x >> (n-54)) % UInt64) & 0x001f_ffff_ffff_ffff # keep 1 extra bit
        y = (y+1)>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt64(trailing_zeros(x) == (n-54)) # fix last bit to round to even
    end
    d = ((n+1022) % UInt64) << 52
    reinterpret(Float64, s | d + y)
end

function Float32(x::UInt128)
    x == 0 && return 0f0
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 24
        y = ((x % UInt32) << (24-n)) & 0x007f_ffff
    else
        y = ((x >> (n-25)) % UInt32) & 0x00ff_ffff # keep 1 extra bit
        y = (y+one(UInt32))>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt32(trailing_zeros(x) == (n-25)) # fix last bit to round to even
    end
    d = ((n+126) % UInt32) << 23
    reinterpret(Float32, d + y)
end

function Float32(x::Int128)
    x == 0 && return 0f0
    s = ((x >>> 96) % UInt32) & 0x8000_0000 # sign bit
    x = abs(x) % UInt128
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 24
        y = ((x % UInt32) << (24-n)) & 0x007f_ffff
    else
        y = ((x >> (n-25)) % UInt32) & 0x00ff_ffff # keep 1 extra bit
        y = (y+one(UInt32))>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt32(trailing_zeros(x) == (n-25)) # fix last bit to round to even
    end
    d = ((n+126) % UInt32) << 23
    reinterpret(Float32, s | d + y)
end

# Float32 -> Float16 algorithm from:
#   "Fast Half Float Conversion" by Jeroen van der Zijp
#   ftp://ftp.fox-toolkit.org/pub/fasthalffloatconversion.pdf
#
# With adjustments for round-to-nearest, ties to even.
#
let _basetable = Vector{UInt16}(undef, 512),
    _shifttable = Vector{UInt8}(undef, 512)
    for i = 0:255
        e = i - 127
        if e < -25  # Very small numbers map to zero
            _basetable[i|0x000+1] = 0x0000
            _basetable[i|0x100+1] = 0x8000
            _shifttable[i|0x000+1] = 25
            _shifttable[i|0x100+1] = 25
        elseif e < -14  # Small numbers map to denorms
            _basetable[i|0x000+1] = 0x0000
            _basetable[i|0x100+1] = 0x8000
            _shifttable[i|0x000+1] = -e-1
            _shifttable[i|0x100+1] = -e-1
        elseif e <= 15  # Normal numbers just lose precision
            _basetable[i|0x000+1] = ((e+15)<<10)
            _basetable[i|0x100+1] = ((e+15)<<10) | 0x8000
            _shifttable[i|0x000+1] = 13
            _shifttable[i|0x100+1] = 13
        elseif e < 128  # Large numbers map to Infinity
            _basetable[i|0x000+1] = 0x7C00
            _basetable[i|0x100+1] = 0xFC00
            _shifttable[i|0x000+1] = 24
            _shifttable[i|0x100+1] = 24
        else  # Infinity and NaN's stay Infinity and NaN's
            _basetable[i|0x000+1] = 0x7C00
            _basetable[i|0x100+1] = 0xFC00
            _shifttable[i|0x000+1] = 13
            _shifttable[i|0x100+1] = 13
        end
    end
    global const shifttable = (_shifttable...,)
    global const basetable = (_basetable...,)
end

function Float16(val::Float32)
    f = reinterpret(UInt32, val)
    if isnan(val)
        t = 0x8000 ⊻ (0x8000 & ((f >> 0x10) % UInt16))
        return reinterpret(Float16, t ⊻ ((f >> 0xd) % UInt16))
    end
    i = ((f & ~significand_mask(Float32)) >> significand_bits(Float32)) + 1
    @inbounds sh = shifttable[i]
    f &= significand_mask(Float32)
    # If `val` is subnormal, the tables are set up to force the
    # result to 0, so the significand has an implicit `1` in the
    # cases we care about.
    f |= significand_mask(Float32) + 0x1
    @inbounds h = (basetable[i] + (f >> sh) & significand_mask(Float16)) % UInt16
    # round
    # NOTE: we maybe should ignore NaNs here, but the payload is
    # getting truncated anyway so "rounding" it might not matter
    nextbit = (f >> (sh-1)) & 1
    if nextbit != 0 && (h & 0x7C00) != 0x7C00
        # Round halfway to even or check lower bits
        if h&1 == 1 || (f & ((1<<(sh-1))-1)) != 0
            h += UInt16(1)
        end
    end
    reinterpret(Float16, h)
end

function Float32(val::Float16)
    local ival::UInt32 = reinterpret(UInt16, val)
    local sign::UInt32 = (ival & 0x8000) >> 15
    local exp::UInt32  = (ival & 0x7c00) >> 10
    local sig::UInt32  = (ival & 0x3ff) >> 0
    local ret::UInt32

    if exp == 0
        if sig == 0
            sign = sign << 31
            ret = sign | exp | sig
        else
            n_bit = 1
            bit = 0x0200
            while (bit & sig) == 0
                n_bit = n_bit + 1
                bit = bit >> 1
            end
            sign = sign << 31
            exp = ((-14 - n_bit + 127) << 23) % UInt32
            sig = ((sig & (~bit)) << n_bit) << (23 - 10)
            ret = sign | exp | sig
        end
    elseif exp == 0x1f
        if sig == 0  # Inf
            if sign == 0
                ret = 0x7f800000
            else
                ret = 0xff800000
            end
        else  # NaN
            ret = 0x7fc00000 | (sign<<31) | (sig<<(23-10))
        end
    else
        sign = sign << 31
        exp  = ((exp - 15 + 127) << 23) % UInt32
        sig  = sig << (23 - 10)
        ret = sign | exp | sig
    end
    return reinterpret(Float32, ret)
end

#convert(::Type{Float16}, x::Float32) = fptrunc(Float16, x)
Float32(x::Float64) = fptrunc(Float32, x)
Float16(x::Float64) = Float16(Float32(x))

#convert(::Type{Float32}, x::Float16) = fpext(Float32, x)
Float64(x::Float32) = fpext(Float64, x)
Float64(x::Float16) = Float64(Float32(x))

AbstractFloat(x::Bool)    = Float64(x)
AbstractFloat(x::Int8)    = Float64(x)
AbstractFloat(x::Int16)   = Float64(x)
AbstractFloat(x::Int32)   = Float64(x)
AbstractFloat(x::Int64)   = Float64(x) # LOSSY
AbstractFloat(x::Int128)  = Float64(x) # LOSSY
AbstractFloat(x::UInt8)   = Float64(x)
AbstractFloat(x::UInt16)  = Float64(x)
AbstractFloat(x::UInt32)  = Float64(x)
AbstractFloat(x::UInt64)  = Float64(x) # LOSSY
AbstractFloat(x::UInt128) = Float64(x) # LOSSY

Bool(x::Float16) = x==0 ? false : x==1 ? true : throw(InexactError(:Bool, Bool, x))

"""
    float(x)

Convert a number or array to a floating point data type.
"""
float(x) = AbstractFloat(x)

"""
    float(T::Type)

Return an appropriate type to represent a value of type `T` as a floating point value.
Equivalent to `typeof(float(zero(T)))`.

# Examples
```jldoctest
julia> float(Complex{Int})
Complex{Float64}

julia> float(Int)
Float64
```
"""
float(::Type{T}) where {T<:Number} = typeof(float(zero(T)))
float(::Type{T}) where {T<:AbstractFloat} = T

"""
    unsafe_trunc(T, x)

Return the nearest integral value of type `T` whose absolute value is
less than or equal to `x`. If the value is not representable by `T`, an arbitrary value will
be returned.
"""
function unsafe_trunc end

for Ti in (Int8, Int16, Int32, Int64)
    @eval begin
        unsafe_trunc(::Type{$Ti}, x::Float16) = unsafe_trunc($Ti, Float32(x))
        unsafe_trunc(::Type{$Ti}, x::Float32) = fptosi($Ti, x)
        unsafe_trunc(::Type{$Ti}, x::Float64) = fptosi($Ti, x)
    end
end
for Ti in (UInt8, UInt16, UInt32, UInt64)
    @eval begin
        unsafe_trunc(::Type{$Ti}, x::Float16) = unsafe_trunc($Ti, Float32(x))
        unsafe_trunc(::Type{$Ti}, x::Float32) = fptoui($Ti, x)
        unsafe_trunc(::Type{$Ti}, x::Float64) = fptoui($Ti, x)
    end
end

function unsafe_trunc(::Type{UInt128}, x::Float64)
    xu = reinterpret(UInt64,x)
    k = Int(xu >> 52) & 0x07ff - 1075
    xu = (xu & 0x000f_ffff_ffff_ffff) | 0x0010_0000_0000_0000
    if k <= 0
        UInt128(xu >> -k)
    else
        UInt128(xu) << k
    end
end
function unsafe_trunc(::Type{Int128}, x::Float64)
    copysign(unsafe_trunc(UInt128,x) % Int128, x)
end

function unsafe_trunc(::Type{UInt128}, x::Float32)
    xu = reinterpret(UInt32,x)
    k = Int(xu >> 23) & 0x00ff - 150
    xu = (xu & 0x007f_ffff) | 0x0080_0000
    if k <= 0
        UInt128(xu >> -k)
    else
        UInt128(xu) << k
    end
end
function unsafe_trunc(::Type{Int128}, x::Float32)
    copysign(unsafe_trunc(UInt128,x) % Int128, x)
end

unsafe_trunc(::Type{UInt128}, x::Float16) = unsafe_trunc(UInt128, Float32(x))
unsafe_trunc(::Type{Int128}, x::Float16) = unsafe_trunc(Int128, Float32(x))

# matches convert methods
# also determines floor, ceil, round
trunc(::Type{Signed}, x::Float32) = trunc(Int,x)
trunc(::Type{Signed}, x::Float64) = trunc(Int,x)
trunc(::Type{Unsigned}, x::Float32) = trunc(UInt,x)
trunc(::Type{Unsigned}, x::Float64) = trunc(UInt,x)
trunc(::Type{Integer}, x::Float32) = trunc(Int,x)
trunc(::Type{Integer}, x::Float64) = trunc(Int,x)
trunc(::Type{T}, x::Float16) where {T<:Integer} = trunc(T, Float32(x))

# fallbacks
floor(::Type{T}, x::AbstractFloat) where {T<:Integer} = trunc(T,round(x, RoundDown))
floor(::Type{T}, x::Float16) where {T<:Integer} = floor(T, Float32(x))
ceil(::Type{T}, x::AbstractFloat) where {T<:Integer} = trunc(T,round(x, RoundUp))
ceil(::Type{T}, x::Float16) where {T<:Integer} = ceil(T, Float32(x))
round(::Type{T}, x::AbstractFloat) where {T<:Integer} = trunc(T,round(x, RoundNearest))
round(::Type{T}, x::Float16) where {T<:Integer} = round(T, Float32(x))

round(x::Float64, r::RoundingMode{:ToZero})  = trunc_llvm(x)
round(x::Float32, r::RoundingMode{:ToZero})  = trunc_llvm(x)
round(x::Float64, r::RoundingMode{:Down})    = floor_llvm(x)
round(x::Float32, r::RoundingMode{:Down})    = floor_llvm(x)
round(x::Float64, r::RoundingMode{:Up})      = ceil_llvm(x)
round(x::Float32, r::RoundingMode{:Up})      = ceil_llvm(x)
round(x::Float64, r::RoundingMode{:Nearest}) = rint_llvm(x)
round(x::Float32, r::RoundingMode{:Nearest}) = rint_llvm(x)

round(x::Float16, r::RoundingMode{:ToZero}) = Float16(round(Float32(x), r))
round(x::Float16, r::RoundingMode{:Down}) = Float16(round(Float32(x), r))
round(x::Float16, r::RoundingMode{:Up}) = Float16(round(Float32(x), r))
round(x::Float16, r::RoundingMode{:Nearest}) = Float16(round(Float32(x), r))

## floating point promotions ##
promote_rule(::Type{Float32}, ::Type{Float16}) = Float32
promote_rule(::Type{Float64}, ::Type{Float16}) = Float64
promote_rule(::Type{Float64}, ::Type{Float32}) = Float64

widen(::Type{Float16}) = Float32
widen(::Type{Float32}) = Float64

## floating point arithmetic ##
-(x::Float64) = neg_float(x)
-(x::Float32) = neg_float(x)
-(x::Float16) = reinterpret(Float16, reinterpret(UInt16, x) ⊻ 0x8000)

for op in (:+, :-, :*, :/, :\, :^)
    @eval ($op)(a::Float16, b::Float16) = Float16(($op)(Float32(a), Float32(b)))
end
+(x::Float32, y::Float32) = add_float(x, y)
+(x::Float64, y::Float64) = add_float(x, y)
-(x::Float32, y::Float32) = sub_float(x, y)
-(x::Float64, y::Float64) = sub_float(x, y)
*(x::Float32, y::Float32) = mul_float(x, y)
*(x::Float64, y::Float64) = mul_float(x, y)
/(x::Float32, y::Float32) = div_float(x, y)
/(x::Float64, y::Float64) = div_float(x, y)

muladd(x::Float32, y::Float32, z::Float32) = muladd_float(x, y, z)
muladd(x::Float64, y::Float64, z::Float64) = muladd_float(x, y, z)
function muladd(a::Float16, b::Float16, c::Float16)
    Float16(muladd(Float32(a), Float32(b), Float32(c)))
end

# TODO: faster floating point div?
# TODO: faster floating point fld?
# TODO: faster floating point mod?

for func in (:div,:fld,:cld,:rem,:mod)
    @eval begin
        $func(a::Float16,b::Float16) = Float16($func(Float32(a),Float32(b)))
    end
end

rem(x::Float32, y::Float32) = rem_float(x, y)
rem(x::Float64, y::Float64) = rem_float(x, y)

cld(x::T, y::T) where {T<:AbstractFloat} = -fld(-x,y)

function mod(x::T, y::T) where T<:AbstractFloat
    r = rem(x,y)
    if r == 0
        copysign(r,y)
    elseif (r > 0) ⊻ (y > 0)
        r+y
    else
        r
    end
end

## floating point comparisons ##
function ==(x::Float16, y::Float16)
    ix = reinterpret(UInt16,x)
    iy = reinterpret(UInt16,y)
    if (ix|iy)&0x7fff > 0x7c00 #isnan(x) || isnan(y)
        return false
    end
    if (ix|iy)&0x7fff == 0x0000
        return true
    end
    return ix == iy
end
==(x::Float32, y::Float32) = eq_float(x, y)
==(x::Float64, y::Float64) = eq_float(x, y)
!=(x::Float32, y::Float32) = ne_float(x, y)
!=(x::Float64, y::Float64) = ne_float(x, y)
<( x::Float32, y::Float32) = lt_float(x, y)
<( x::Float64, y::Float64) = lt_float(x, y)
<=(x::Float32, y::Float32) = le_float(x, y)
<=(x::Float64, y::Float64) = le_float(x, y)

isequal(x::Float32, y::Float32) = fpiseq(x, y)
isequal(x::Float64, y::Float64) = fpiseq(x, y)
isless( x::Float32, y::Float32) = fpislt(x, y)
isless( x::Float64, y::Float64) = fpislt(x, y)
for op in (:<, :<=, :isless)
    @eval ($op)(a::Float16, b::Float16) = ($op)(Float32(a), Float32(b))
end

# Exact Float (Tf) vs Integer (Ti) comparisons
# Assumes:
# - typemax(Ti) == 2^n-1
# - typemax(Ti) can't be exactly represented by Tf:
#   => Tf(typemax(Ti)) == 2^n or Inf
# - typemin(Ti) can be exactly represented by Tf
#
# 1. convert y::Ti to float fy::Tf
# 2. perform Tf comparison x vs fy
# 3. if x == fy, check if (1) resulted in rounding:
#  a. convert fy back to Ti and compare with original y
#  b. unsafe_convert undefined behaviour if fy == Tf(typemax(Ti))
#     (but consequently x == fy > y)
for Ti in (Int64,UInt64,Int128,UInt128)
    for Tf in (Float32,Float64)
        @eval begin
            function ==(x::$Tf, y::$Ti)
                fy = ($Tf)(y)
                (x == fy) & (fy != $(Tf(typemax(Ti)))) & (y == unsafe_trunc($Ti,fy))
            end
            ==(y::$Ti, x::$Tf) = x==y

            function <(x::$Ti, y::$Tf)
                fx = ($Tf)(x)
                (fx < y) | ((fx == y) & ((fx == $(Tf(typemax(Ti)))) | (x < unsafe_trunc($Ti,fx)) ))
            end
            function <=(x::$Ti, y::$Tf)
                fx = ($Tf)(x)
                (fx < y) | ((fx == y) & ((fx == $(Tf(typemax(Ti)))) | (x <= unsafe_trunc($Ti,fx)) ))
            end

            function <(x::$Tf, y::$Ti)
                fy = ($Tf)(y)
                (x < fy) | ((x == fy) & (fy < $(Tf(typemax(Ti)))) & (unsafe_trunc($Ti,fy) < y))
            end
            function <=(x::$Tf, y::$Ti)
                fy = ($Tf)(y)
                (x < fy) | ((x == fy) & (fy < $(Tf(typemax(Ti)))) & (unsafe_trunc($Ti,fy) <= y))
            end
        end
    end
end
for op in (:(==), :<, :<=)
    @eval begin
        ($op)(x::Float16, y::Union{Int128,UInt128,Int64,UInt64}) = ($op)(Float64(x), Float64(y))
        ($op)(x::Union{Int128,UInt128,Int64,UInt64}, y::Float16) = ($op)(Float64(x), Float64(y))

        ($op)(x::Union{Float16,Float32}, y::Union{Int32,UInt32}) = ($op)(Float64(x), Float64(y))
        ($op)(x::Union{Int32,UInt32}, y::Union{Float16,Float32}) = ($op)(Float64(x), Float64(y))

        ($op)(x::Float16, y::Union{Int16,UInt16}) = ($op)(Float32(x), Float32(y))
        ($op)(x::Union{Int16,UInt16}, y::Float16) = ($op)(Float32(x), Float32(y))
    end
end

abs(x::Float16) = reinterpret(Float16, reinterpret(UInt16, x) & 0x7fff)
abs(x::Float32) = abs_float(x)
abs(x::Float64) = abs_float(x)

"""
    isnan(f) -> Bool

Test whether a number value is a NaN, an indeterminate value which is neither an infinity
nor a finite number ("not a number").
"""
isnan(x::AbstractFloat) = x != x
isnan(x::Float16) = reinterpret(UInt16,x)&0x7fff > 0x7c00
isnan(x::Real) = false

"""
    isfinite(f) -> Bool

Test whether a number is finite.

# Examples
```jldoctest
julia> isfinite(5)
true

julia> isfinite(NaN32)
false
```
"""
isfinite(x::AbstractFloat) = x - x == 0
isfinite(x::Float16) = reinterpret(UInt16,x)&0x7c00 != 0x7c00
isfinite(x::Real) = decompose(x)[3] != 0
isfinite(x::Integer) = true

"""
    isinf(f) -> Bool

Test whether a number is infinite.
"""
isinf(x::Real) = !isnan(x) & !isfinite(x)

## hashing small, built-in numeric types ##

hx(a::UInt64, b::Float64, h::UInt) = hash_uint64((3a + reinterpret(UInt64,b)) - h)
const hx_NaN = hx(UInt64(0), NaN, UInt(0  ))

hash(x::UInt64,  h::UInt) = hx(x, Float64(x), h)
hash(x::Int64,   h::UInt) = hx(reinterpret(UInt64, abs(x)), Float64(x), h)
hash(x::Float64, h::UInt) = isnan(x) ? (hx_NaN ⊻ h) : hx(fptoui(UInt64, abs(x)), x, h)

hash(x::Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32}, h::UInt) = hash(Int64(x), h)
hash(x::Float32, h::UInt) = hash(Float64(x), h)

"""
    precision(num::AbstractFloat)

Get the precision of a floating point number, as defined by the effective number of bits in
the significand.
"""
function precision end

precision(::Type{Float16}) = 11
precision(::Type{Float32}) = 24
precision(::Type{Float64}) = 53
precision(::T) where {T<:AbstractFloat} = precision(T)

"""
    uabs(x::Integer)

Return the absolute value of `x`, possibly returning a different type should the
operation be susceptible to overflow. This typically arises when `x` is a two's complement
signed integer, so that `abs(typemin(x)) == typemin(x) < 0`, in which case the result of
`uabs(x)` will be an unsigned integer of the same size.
"""
uabs(x::Integer) = abs(x)
uabs(x::BitSigned) = unsigned(abs(x))

"""
    nextfloat(x::AbstractFloat, n::Integer)

The result of `n` iterative applications of `nextfloat` to `x` if `n >= 0`, or `-n`
applications of `prevfloat` if `n < 0`.
"""
function nextfloat(f::IEEEFloat, d::Integer)
    F = typeof(f)
    fumax = reinterpret(Unsigned, F(Inf))
    U = typeof(fumax)

    isnan(f) && return f
    fi = reinterpret(Signed, f)
    fneg = fi < 0
    fu = unsigned(fi & typemax(fi))

    dneg = d < 0
    da = uabs(d)
    if da > typemax(U)
        fneg = dneg
        fu = fumax
    else
        du = da % U
        if fneg ⊻ dneg
            if du > fu
                fu = min(fumax, du - fu)
                fneg = !fneg
            else
                fu = fu - du
            end
        else
            if fumax - fu < du
                fu = fumax
            else
                fu = fu + du
            end
        end
    end
    if fneg
        fu |= sign_mask(F)
    end
    reinterpret(F, fu)
end

"""
    nextfloat(x::AbstractFloat)

Return the smallest floating point number `y` of the same type as `x` such `x < y`. If no
such `y` exists (e.g. if `x` is `Inf` or `NaN`), then return `x`.
"""
nextfloat(x::AbstractFloat) = nextfloat(x,1)

"""
    prevfloat(x::AbstractFloat, n::Integer)

The result of `n` iterative applications of `prevfloat` to `x` if `n >= 0`, or `-n`
applications of `nextfloat` if `n < 0`.
"""
prevfloat(x::AbstractFloat, d::Integer) = nextfloat(x, -d)

"""
    prevfloat(x::AbstractFloat)

Return the largest floating point number `y` of the same type as `x` such `y < x`. If no
such `y` exists (e.g. if `x` is `-Inf` or `NaN`), then return `x`.
"""
prevfloat(x::AbstractFloat) = nextfloat(x,-1)

for Ti in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
    for Tf in (Float32, Float64)
        if Ti <: Unsigned || sizeof(Ti) < sizeof(Tf)
            # Here `Tf(typemin(Ti))-1` is exact, so we can compare the lower-bound
            # directly. `Tf(typemax(Ti))+1` is either always exactly representable, or
            # rounded to `Inf` (e.g. when `Ti==UInt128 && Tf==Float32`).
            @eval begin
                function trunc(::Type{$Ti},x::$Tf)
                    if $(Tf(typemin(Ti))-one(Tf)) < x < $(Tf(typemax(Ti))+one(Tf))
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError(:trunc, $Ti, x))
                    end
                end
                function (::Type{$Ti})(x::$Tf)
                    if ($(Tf(typemin(Ti))) <= x <= $(Tf(typemax(Ti)))) && (round(x, RoundToZero) == x)
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                    end
                end
            end
        else
            # Here `eps(Tf(typemin(Ti))) > 1`, so the only value which can be truncated to
            # `Tf(typemin(Ti)` is itself. Similarly, `Tf(typemax(Ti))` is inexact and will
            # be rounded up. This assumes that `Tf(typemin(Ti)) > -Inf`, which is true for
            # these types, but not for `Float16` or larger integer types.
            @eval begin
                function trunc(::Type{$Ti},x::$Tf)
                    if $(Tf(typemin(Ti))) <= x < $(Tf(typemax(Ti)))
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError(:trunc, $Ti, x))
                    end
                end
                function (::Type{$Ti})(x::$Tf)
                    if ($(Tf(typemin(Ti))) <= x < $(Tf(typemax(Ti)))) && (round(x, RoundToZero) == x)
                        return unsafe_trunc($Ti,x)
                    else
                        throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                    end
                end
            end
        end
    end
end

"""
    issubnormal(f) -> Bool

Test whether a floating point number is subnormal.
"""
function issubnormal(x::T) where {T<:IEEEFloat}
    y = reinterpret(Unsigned, x)
    (y & exponent_mask(T) == 0) & (y & significand_mask(T) != 0)
end

@eval begin
    typemin(::Type{Float16}) = $(bitcast(Float16, 0xfc00))
    typemax(::Type{Float16}) = $(Inf16)
    typemin(::Type{Float32}) = $(-Inf32)
    typemax(::Type{Float32}) = $(Inf32)
    typemin(::Type{Float64}) = $(-Inf64)
    typemax(::Type{Float64}) = $(Inf64)
    typemin(x::T) where {T<:Real} = typemin(T)
    typemax(x::T) where {T<:Real} = typemax(T)

    floatmin(::Type{Float16}) = $(bitcast(Float16, 0x0400))
    floatmin(::Type{Float32}) = $(bitcast(Float32, 0x00800000))
    floatmin(::Type{Float64}) = $(bitcast(Float64, 0x0010000000000000))
    floatmax(::Type{Float16}) = $(bitcast(Float16, 0x7bff))
    floatmax(::Type{Float32}) = $(bitcast(Float32, 0x7f7fffff))
    floatmax(::Type{Float64}) = $(bitcast(Float64, 0x7fefffffffffffff))

    eps(x::AbstractFloat) = isfinite(x) ? abs(x) >= floatmin(x) ? ldexp(eps(typeof(x)), exponent(x)) : nextfloat(zero(x)) : oftype(x, NaN)
    eps(::Type{Float16}) = $(bitcast(Float16, 0x1400))
    eps(::Type{Float32}) = $(bitcast(Float32, 0x34000000))
    eps(::Type{Float64}) = $(bitcast(Float64, 0x3cb0000000000000))
    eps() = eps(Float64)
end

"""
    floatmin(T)

The smallest in absolute value non-subnormal value representable by the given
floating-point DataType `T`.
"""
floatmin(x::T) where {T<:AbstractFloat} = floatmin(T)

"""
    floatmax(T)

The highest finite value representable by the given floating-point DataType `T`.

# Examples
```jldoctest
julia> floatmax(Float16)
Float16(6.55e4)

julia> floatmax(Float32)
3.4028235f38
```
"""
floatmax(x::T) where {T<:AbstractFloat} = floatmax(T)

floatmin() = floatmin(Float64)
floatmax() = floatmax(Float64)

"""
    eps(::Type{T}) where T<:AbstractFloat
    eps()

Return the *machine epsilon* of the floating point type `T` (`T = Float64` by
default). This is defined as the gap between 1 and the next largest value representable by
`typeof(one(T))`, and is equivalent to `eps(one(T))`.  (Since `eps(T)` is a
bound on the *relative error* of `T`, it is a "dimensionless" quantity like [`one`](@ref).)

# Examples
```jldoctest
julia> eps()
2.220446049250313e-16

julia> eps(Float32)
1.1920929f-7

julia> 1.0 + eps()
1.0000000000000002

julia> 1.0 + eps()/2
1.0
```
"""
eps(::Type{<:AbstractFloat})

"""
    eps(x::AbstractFloat)

Return the *unit in last place* (ulp) of `x`. This is the distance between consecutive
representable floating point values at `x`. In most cases, if the distance on either side
of `x` is different, then the larger of the two is taken, that is

    eps(x) == max(x-prevfloat(x), nextfloat(x)-x)

The exceptions to this rule are the smallest and largest finite values
(e.g. `nextfloat(-Inf)` and `prevfloat(Inf)` for [`Float64`](@ref)), which round to the
smaller of the values.

The rationale for this behavior is that `eps` bounds the floating point rounding
error. Under the default `RoundNearest` rounding mode, if ``y`` is a real number and ``x``
is the nearest floating point number to ``y``, then

```math
|y-x| \\leq \\operatorname{eps}(x)/2.
```

# Examples
```jldoctest
julia> eps(1.0)
2.220446049250313e-16

julia> eps(prevfloat(2.0))
2.220446049250313e-16

julia> eps(2.0)
4.440892098500626e-16

julia> x = prevfloat(Inf)      # largest finite Float64
1.7976931348623157e308

julia> x + eps(x)/2            # rounds up
Inf

julia> x + prevfloat(eps(x)/2) # rounds down
1.7976931348623157e308
```
"""
eps(::AbstractFloat)

## byte order swaps for arbitrary-endianness serialization/deserialization ##
bswap(x::IEEEFloat) = bswap_int(x)

# bit patterns
reinterpret(::Type{Unsigned}, x::Float64) = reinterpret(UInt64, x)
reinterpret(::Type{Unsigned}, x::Float32) = reinterpret(UInt32, x)
reinterpret(::Type{Signed}, x::Float64) = reinterpret(Int64, x)
reinterpret(::Type{Signed}, x::Float32) = reinterpret(Int32, x)

sign_mask(::Type{Float64}) =        0x8000_0000_0000_0000
exponent_mask(::Type{Float64}) =    0x7ff0_0000_0000_0000
exponent_one(::Type{Float64}) =     0x3ff0_0000_0000_0000
exponent_half(::Type{Float64}) =    0x3fe0_0000_0000_0000
significand_mask(::Type{Float64}) = 0x000f_ffff_ffff_ffff

sign_mask(::Type{Float32}) =        0x8000_0000
exponent_mask(::Type{Float32}) =    0x7f80_0000
exponent_one(::Type{Float32}) =     0x3f80_0000
exponent_half(::Type{Float32}) =    0x3f00_0000
significand_mask(::Type{Float32}) = 0x007f_ffff

sign_mask(::Type{Float16}) =        0x8000
exponent_mask(::Type{Float16}) =    0x7c00
exponent_one(::Type{Float16}) =     0x3c00
exponent_half(::Type{Float16}) =    0x3800
significand_mask(::Type{Float16}) = 0x03ff

for T in (Float16, Float32, Float64)
    @eval significand_bits(::Type{$T}) = $(trailing_ones(significand_mask(T)))
    @eval exponent_bits(::Type{$T}) = $(sizeof(T)*8 - significand_bits(T) - 1)
    @eval exponent_bias(::Type{$T}) = $(Int(exponent_one(T) >> significand_bits(T)))
    # maximum float exponent
    @eval exponent_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)) - exponent_bias(T))
    # maximum float exponent without bias
    @eval exponent_raw_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)))
end

# integer size of float
uinttype(::Type{Float64}) = UInt64
uinttype(::Type{Float32}) = UInt32
uinttype(::Type{Float16}) = UInt16

Base.iszero(x::Float16) = reinterpret(UInt16, x) & ~sign_mask(Float16) == 0x0000

## Array operations on floating point numbers ##

float(A::AbstractArray{<:AbstractFloat}) = A

function float(A::AbstractArray{T}) where T
    if !isconcretetype(T)
        error("`float` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(float(zero(T)))}, A)
end

float(r::StepRange) = float(r.start):float(r.step):float(last(r))
float(r::UnitRange) = float(r.start):float(last(r))
float(r::StepRangeLen{T}) where {T} =
    StepRangeLen{typeof(float(T(r.ref)))}(float(r.ref), float(r.step), length(r), r.offset)
function float(r::LinRange)
    LinRange(float(r.start), float(r.stop), length(r))
end

max_stars_repo_name: Dictino/julia

max_issues_repo_name: Dictino/julia

max_forks_repo_name: Dictino/julia


text:
@testset "334.increasing-triplet-subsequence.jl" begin
    @test increasing_triplet(collect(1:5)) == true
    @test increasing_triplet([5, 4, 3, 2, 1]) == false
    @test increasing_triplet([2, 2, 2, 1, 5, 0, 4, 6]) == true
end

ext: jl


text:
module Conductor

using ModelingToolkit,
      Unitful,
      Unitful.DefaultSymbols,
      InteractiveUtils,
      IfElse,
      Symbolics,
      SymbolicUtils,
      Setfield

import Symbolics:
    get_variables,
    Symbolic,
    value,
    tosymbol,
    VariableDefaultValue,
    wrap

import ModelingToolkit:
    toparam,
    isparameter,
    Equation,
    defaults,
    AbstractSystem,
    get_eqs,
    get_states,
    get_observed,
    get_defaults,
    get_ps,
    get_systems

import ModelingToolkit.SciMLBase: parameterless_type

import Unitful:
    Time,
    Voltage,
    Current,
    Molarity,
    ElectricalConductance

import SymbolicUtils: FnType
import Unitful: mV, mS, cm, µF, mF, µm, pA, nA, mA, µA, ms, mM, µM
import Base: show, display

export Gate, AlphaBetaRates, SteadyStateTau, IonChannel, PassiveChannel, SynapticChannel
export EquilibriumPotential, Equilibrium, Equilibria, MembranePotential, MembraneCurrent
export AuxConversion, D, Network
export Soma, Simulation, Concentration, IonConcentration
export @named
export Calcium, Sodium, Potassium, Chloride, Cation, Anion, Leak, Ion

const ℱ = Unitful.q*Unitful.Na # Faraday's constant
const t = let name = :t; only(@parameters $name) end
const D = Differential(t)

# Temporary fix until https://github.com/SciML/ModelingToolkit.jl/issues/1223
# gets resolved. Non-flattening form of extend
function _extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys))
    T = parameterless_type(basesys)
    ivs = independent_variables(basesys)
    if !(typeof(sys) <: T)
        if length(ivs) == 0
            sys = convert_system(T, sys)
        elseif length(ivs) == 1
            sys = convert_system(T, sys, ivs[1])
        else
            throw("Extending multivariate systems is not supported")
        end
    end

    eqs = union(get_eqs(basesys), get_eqs(sys))
    sts = union(get_states(basesys), get_states(sys))
    ps = union(get_ps(basesys), get_ps(sys))
    obs = union(get_observed(basesys), get_observed(sys))
    defs = merge(get_defaults(basesys), get_defaults(sys)) # prefer `sys`
    syss = union(get_systems(basesys), get_systems(sys))

    if length(ivs) == 0
        T(eqs, sts, ps, observed = obs, defaults = defs, name=name, systems = syss)
    elseif length(ivs) == 1
        T(eqs, ivs[1], sts, ps, observed = obs, defaults = defs, name = name, systems = syss)
    end
end

# Helper utils -- FIXME: MTK has these now
#hasdefault(x::Symbolic) = hasmetadata(x, VariableDefaultValue) ? true : false
#hasdefault(x::Num) = hasdefault(ModelingToolkit.value(x))
#hasdefault(x) = false
#
#getdefault(x::Symbolic) = hasdefault(x) ? getmetadata(x, VariableDefaultValue) : nothing
#getdefault(x::Num) = getdefault(ModelingToolkit.value(x))

# Basic symbols
function MembranePotential()
    name = :Vₘ
    return only(@variables $name(t))
end

@enum Location Outside Inside

# Custom Unitful.jl quantities
@derived_dimension SpecificConductance 𝐈^2*𝐋^-4*𝐌^-1*𝐓^3 # conductance per unit area
@derived_dimension SpecificCapacitance 𝐈^2*𝐋^-4*𝐌^-1*𝐓^4 # capacitance per unit area
@derived_dimension ConductancePerFarad 𝐓^-1 # S/F cancels out to 1/s; perhaps find a better abstract type?

# Metadata IDs
abstract type ConductorCurrentCtx end
abstract type ConductorEquilibriumCtx end
abstract type ConductorConcentrationCtx end
abstract type ConductorAggregatorCtx end

# Ion species
abstract type Ion end
abstract type Cation <: Ion end
abstract type Anion <: Ion end
struct Calcium  <: Ion end
struct Sodium <: Ion end
struct Potassium <: Ion end
struct Chloride <: Ion end

const Ca = Calcium
const Na = Sodium
const K = Potassium
const Cl = Chloride
const Mixed = Ion # non-specific ion
const Leak = Mixed

const PERIODIC_SYMBOL = IdDict(Na => :Na, K  => :K, Cl => :Cl, Ca => :Ca, Leak => :l)

# Concentrations of ions
abstract type AbstractConcentration end

struct IonConcentration{I<:Ion, L<:Location, V<:Union{Nothing, Num,Symbolic,Molarity}} <:AbstractConcentration
    ion::Type{I}
    val::V
    loc::L
end

# FIXME: handle default values better
function Concentration(::Type{I}, val = nothing, loc::Location = Inside, name::Symbol = PERIODIC_SYMBOL[I]) where {I <: Ion}
    sym = Symbol(name,(loc == Inside ? "ᵢ" : "ₒ"))
    # FIXME: Not necessarily a parameter when used as a primitive...but we should support
    # this. Use a flag?
    var = #=val isa Molarity ? only(@parameters $sym) :=# only(@variables $sym(t))
    return setmetadata(var,  ConductorConcentrationCtx, IonConcentration(I, val, loc))
end

isconcentration(x::Symbolic) = hasmetadata(x, ConductorConcentrationCtx)
isconcentration(x::Num) = isconcentration(value(x))
getconcentration(x::Symbolic) = isconcentration(x) ? getmetadata(x, ConductorConcentrationCtx) : nothing
getconcentration(x::Num) = getconcentration(value(x))

# Currents
struct MembraneCurrent{I<:Ion,V<:Union{Nothing,Num,Symbolic,Current}}
    ion::Type{I}
    val::V
end

# TODO: add aggregator as field of Membrane current struct
function MembraneCurrent{I}(val = nothing; name::Symbol = PERIODIC_SYMBOL[I], aggregate::Bool = false) where {I <: Ion}
    sym = Symbol("I", name)
    var = val isa Current ? only(@parameters $sym) : only(@variables $sym(t))
    var = setmetadata(var, ConductorCurrentCtx, MembraneCurrent(I, val))
    return setmetadata(var, ConductorAggregatorCtx, aggregate)
end

ismembranecurrent(x::Symbolic) = hasmetadata(x, ConductorCurrentCtx)
ismembranecurrent(x::Num) = ismembranecurrent(ModelingToolkit.value(x))
getmembranecurrent(x::Union{Num, Symbolic}) = ismembranecurrent(x) ? getmetadata(x, ConductorCurrentCtx) : nothing
iontype(x::Union{Num, Symbolic}) = getmembranecurrent(x).ion
isaggregator(x::Union{Num, Symbolic})  = getmetadata(x, ConductorAggregatorCtx)

# Equilibrium potential implicitly defines an ionic gradient
abstract type AbstractIonGradient end

struct EquilibriumPotential{I<:Ion,V<:Union{Num,Symbolic,Voltage}} <: AbstractIonGradient
    ion::Type{I}
    val::V
end

const Equilibrium{I} = EquilibriumPotential{I}

function EquilibriumPotential{I}(val, name::Symbol = PERIODIC_SYMBOL[I]) where {I <: Ion}
    sym = Symbol("E", name)
    var = val isa Voltage ? only(@parameters $sym) : only(@variables $sym(t))
    return setmetadata(var, ConductorEquilibriumCtx, EquilibriumPotential(I, val))
end

# Alternate constructor
function Equilibria(equil::Vector)
    out = Num[]
    for x in equil
        !(x.first <: Ion) && throw("Equilibrium potential must be associated with an ion type.")
        if typeof(x.second) <: Tuple
            tup = x.second
            typeof(tup[2]) !== Symbol && throw("Second tuple argument for $(x.first) must be a symbol.")
            push!(out, Equilibrium{x.first}(tup...))
        else
            push!(out, Equilibrium{x.first}(x.second))
        end
    end
    return out
end

# Gating variables (as an interface)
abstract type AbstractGatingVariable end

hassteadystate(x::AbstractGatingVariable) = hasfield(typeof(x), :ss) ? !(isnothing(x.ss)) : false
hasexponent(x::AbstractGatingVariable) = hasfield(typeof(x), :p) ? x.p !== one(Float64) : false
getsymbol(x::AbstractGatingVariable) = x.sym
getequation(x::AbstractGatingVariable) = x.df

struct Gate <: AbstractGatingVariable
    sym::Num # symbol/name (e.g. m, h)
    df::Equation # differential equation
    ss::Union{Nothing, Num} # optional steady-state expression for initialization
    p::Float64 # optional exponent (defaults to 1)
end

abstract type AbstractGateModel end
struct SteadyStateTau <: AbstractGateModel end
struct AlphaBetaRates <: AbstractGateModel end

function Gate(::Type{SteadyStateTau}, name::Symbol, ss::Num, tau::Num, p::Real)
    sym = only(@variables $name(t))
    df = D(sym) ~ (ss-sym)/tau # (m∞ - m)/τₘ
    return Gate(sym, df, ss, p)
end

function Gate(::Type{AlphaBetaRates}, name::Symbol, alpha::Num, beta::Num, p::Float64)
    sym = only(@variables $name(t))
    df = D(sym) ~ alpha * (1 - sym) - beta*sym # αₘ(1 - m) - βₘ*m
    ss = alpha/(alpha + beta) # αₘ/(αₘ + βₘ)
    return Gate(sym, df, ss, p)
end

Gate(t::Type{AlphaBetaRates}, name::Symbol, alpha::Num, beta::Num, p::Real) =
Gate(AlphaBetaRates, name, alpha, beta, Float64(p))

# TODO: find a nicer way to do this
function Gate(::Type{SteadyStateTau}; p = one(Float64), kwargs...)
    syms = keys(kwargs)
    if length(syms) !== 2
        throw("Invalid number of input equations.")
    elseif issetequal([:m∞, :τₘ], syms)
        Gate(SteadyStateTau, :m, kwargs[:m∞], kwargs[:τₘ], p)
    elseif issetequal([:h∞, :τₕ], syms)
        Gate(SteadyStateTau, :h, kwargs[:h∞], kwargs[:τₕ], p)
    else
        throw("invalid keywords")
    end
end

function Gate(::Type{AlphaBetaRates}; p = one(Float64), kwargs...)
    syms = keys(kwargs)
    if length(syms) !== 2
        throw("Invalid number of input equations.")
    elseif issetequal([:αₘ, :βₘ], syms)
        Gate(AlphaBetaRates, :m, kwargs[:αₘ], kwargs[:βₘ], p)
    elseif issetequal([:αₕ, :βₕ], syms)
        Gate(AlphaBetaRates, :h, kwargs[:αₕ], kwargs[:βₕ], p)
    elseif issetequal([:αₙ, :βₙ], syms)
        Gate(AlphaBetaRates, :n, kwargs[:αₙ], kwargs[:βₙ], p)
    else
        throw("invalid keywords")
    end
end

mutable struct AuxConversion
    params::Vector{Num}
    eqs::Vector{Equation}
end

# Conductance types (conductance as in "g")
abstract type AbstractConductance end

isbuilt(x::AbstractConductance) = x.sys !== nothing

struct IonChannel <: AbstractConductance
    gbar::SpecificConductance # scaling term - maximal conductance per unit area
    conducts::DataType # ion permeability
    inputs::Vector{Num} # cell states dependencies (input args to kinetics); we can infer this
    params::Vector{Num}
    kinetics::Vector{<:AbstractGatingVariable} # gating functions; none = passive channel
    sys::Union{ODESystem, Nothing} # symbolic system
end

# Return ODESystem pretty printing for our wrapper types
Base.show(io::IO, ::MIME"text/plain", x::IonChannel) = Base.display(isbuilt(x) ? x.sys : x)

function _conductance(gbar_val::T, gate_vars::Vector{<:AbstractGatingVariable};
                      passive::Bool = false, null_init::Bool = false, name::Symbol) where {T <: Real}
    inputs = Set{Num}()
    states = Set{Num}()
    eqs = Equation[]
    # retrieve all variables present in the RHS of kinetics equations
    # g = total conductance (e.g. g(m,h) ~ ̄gm³h)
    if passive
        params = @parameters g
        defaultmap = Pair[g => gbar_val]
    else
        gates = Set{Num}(getsymbol(x) for x in gate_vars)
        @variables g(t)
        push!(states, g)
        params = @parameters gbar
        defaultmap = Pair[gbar => gbar_val]

        for i in gate_vars
            syms = value.(get_variables(getequation(i)))
            for j in syms
                if j ∉ gates
                    isparameter(j) ? push!(params, j) : push!(inputs, j)
                end
            end
        end
        union!(states, gates, inputs)
        push!(eqs, g ~ gbar * prod(hasexponent(x) ? getsymbol(x)^x.p : getsymbol(x) for x in gate_vars))
        append!(eqs, getequation(x) for x in gate_vars)
        append!(defaultmap, getsymbol(x) => (hassteadystate(x) && !null_init) ? x.ss : 0.0 for x in gate_vars) # fallback to zero
    end
    system = ODESystem(eqs, t, states, params; defaults = defaultmap, name = name)
    return (collect(inputs), params, system)
end

# General purpose constructor
function IonChannel(conducts::Type{I},
                    gate_vars::Vector{<:AbstractGatingVariable},
                    max_g::SpecificConductance = 0mS/cm^2;
                    passive::Bool = false, name::Symbol) where {I <: Ion}
    # TODO: Generalize to other possible units (e.g. S/F)
    gbar_val = ustrip(Float64, mS/cm^2, max_g)
    (inputs, params, system) = _conductance(gbar_val, gate_vars, passive = passive, name = name)
    return IonChannel(max_g, conducts, inputs, params, gate_vars, system)
end

function (chan::IonChannel)(newgbar::SpecificConductance)
    newchan = deepcopy(chan)
    @set! newchan.gbar = newgbar
    gbar_val = ustrip(Float64, mS/cm^2, newgbar)
    gsym = length(newchan.kinetics) > 0 ?
           value(first(@parameters gbar)) :
           value(first(@parameters g))

    mapping = Dict([gsym => gbar_val])
    new_defaults = merge(defaults(newchan.sys), mapping)
    @set! newchan.sys.defaults = new_defaults
    return newchan
end

# Alias for ion channel with static conductance
function PassiveChannel(conducts::Type{I}, max_g::SpecificConductance = 0mS/cm^2;
                        name::Symbol = Base.gensym(:Leak)) where {I <: Ion}
    gate_vars = AbstractGatingVariable[]
    return IonChannel(conducts, gate_vars, max_g; name = name, passive = true)
end

struct SynapticChannel <: AbstractConductance
    gbar::ElectricalConductance
    conducts::DataType
    reversal::Num
    inputs::Vector{Num}
    params::Vector{Num}
    kinetics::Vector{<:AbstractGatingVariable}
    sys::Union{ODESystem, Nothing}
end

function SynapticChannel(conducts::Type{I}, gate_vars::Vector{<:AbstractGatingVariable},
                         reversal::Num, max_g::ElectricalConductance = 0mS;
                         passive::Bool = false, name::Symbol) where {I <: Ion}
    gbar_val = ustrip(Float64, mS, max_g)
    (inputs, params, system) = _conductance(gbar_val, gate_vars, passive = passive, null_init = true, name = name)
    return SynapticChannel(max_g, conducts, reversal, inputs, params, gate_vars, system)
end

function GapJunction(conducts::Type{I}, reversal::Num, max_g::ElectricalConductance = 0mS;
                     passive::Bool = false, name::Symbol) where {I <: Ion}
    SynapticChannel(conducts, AbstractGatingVariable[], reversal, max_g, passive = true, name)
end

function (chan::SynapticChannel)(newgbar::ElectricalConductance)
    newchan = deepcopy(chan)
    @set! newchan.gbar = newgbar
    gbar_val = ustrip(Float64, mS, newgbar)
    gsym = length(newchan.kinetics) > 0 ?
           value(first(@parameters gbar)) :
           value(first(@parameters g))

    mapping = Dict([gsym => gbar_val])
    new_defaults = merge(defaults(newchan.sys), mapping)
    @set! newchan.sys.defaults = new_defaults
    return newchan
end

abstract type Geometry end
abstract type Sphere <: Geometry end
abstract type Cylinder <: Geometry end

struct Compartment{G}
    cap::SpecificCapacitance
    chans::Vector{<:AbstractConductance}
    states::Vector
    params::Vector
    sys::ODESystem
end

function Compartment{Sphere}(channels::Vector{<:AbstractConductance},
                             gradients; name::Symbol, area::Float64 = 0.628e-3, #radius = 20µm,
                             capacitance::SpecificCapacitance = 1µF/cm^2,
                             V0::Voltage = -65mV,
                             holding::Current = 0nA,
                             stimulus::Union{Function,Nothing} = nothing,
                             aux::Union{Nothing, Vector{AuxConversion}} = nothing)

    Vₘ = MembranePotential()
    @variables Iapp(t) Isyn(t)
    params = @parameters cₘ aₘ
    grad_meta = getmetadata.(gradients, ConductorEquilibriumCtx)
    #r_val = ustrip(Float64, cm, radius) # FIXME: make it so we calculate area from dims as needed

    channel_systems = AbstractSystem[]
    eqs = Equation[] # equations must be a vector
    required_states = [] # states not produced or intrinsic (e.g. not currents or Vm)
    states = Any[Vₘ, Iapp, Isyn] # grow this as we discover/generate new states
    currents = []
    defaultmap = Pair[Iapp => ustrip(Float64, µA, holding),
                      aₘ => area,
                      Vₘ => ustrip(Float64, mV, V0),
                      Isyn => 0,
                      cₘ => ustrip(Float64, mF/cm^2, capacitance)]

    # By default, applied current is constant (just a bias/offset/holding current)
    # TODO: lift "stimulus" to a pass that happens at a higher level (i.e. after neuron
    # construction; with its own data type)
    if stimulus == nothing
        append!(eqs, [D(Iapp) ~ 0])
    else
        push!(eqs, Iapp ~ stimulus(t,Iapp))
    end

    # auxillary state transformations (e.g. net calcium current -> Ca concentration)
    if aux !== nothing
        for i in aux
            append!(params, i.params)
            # TODO: Probably can eliminate this block; MTK should find it automatically
            #for x in i.params
            #    hasdefault(x) && push!(defaultmap, x => getdefault(x))
            #end

            # gather all unique variables)
            inpvars = value.(vcat((get_variables(x.rhs) for x in i.eqs)...))
            unique!(inpvars)
            filter!(x -> !isparameter(x), inpvars) # exclude parameters
            append!(required_states, inpvars)

            # isolate states produced
            outvars = vcat((get_variables(x.lhs) for x in i.eqs)...)
            append!(states, outvars)
            append!(eqs, i.eqs)
            for j in outvars
                # FIXME: consider more consistent use of default variable ctx
                isconcentration(j) && push!(defaultmap, j => ustrip(Float64, µM, getconcentration(j).val))
            end
        end
    end

    # parse and build channel equations
    for chan in channels
        ion = chan.conducts
        sys = chan.sys
        push!(channel_systems, sys)

        # auto forward cell states to channels
        for inp in chan.inputs
            push!(required_states, inp)
            subinp = getproperty(sys, tosymbol(inp, escape=false))
            push!(eqs, inp ~ subinp)
            # Workaround for: https://github.com/SciML/ModelingToolkit.jl/issues/1013
            push!(defaultmap, subinp => inp)
        end

        # write the current equation state
        I = MembraneCurrent{ion}(name = nameof(sys), aggregate = false)
        push!(states, I)
        push!(currents, I)

        # for now, take the first reversal potential with matching ion type
        idx = findfirst(x -> x.ion == ion, grad_meta)
        Erev = gradients[idx]
        eq = [I ~ aₘ * sys.g * (Vₘ - Erev)]
        rhs = grad_meta[idx].val

        # check to see if reversal potential already defined
        if any(isequal(Erev, x) for x in states)
            append!(eqs, eq)
        else
            if typeof(rhs) <: Voltage
                push!(defaultmap, Erev => ustrip(Float64, mV, rhs))
                push!(params, Erev)
                append!(eqs, eq)
            else # symbolic/dynamic reversal potentials
                push!(eq, Erev ~ rhs)
                push!(states, Erev)
                rhs_vars = get_variables(rhs)
                filter!(x -> !isequal(x, value(Erev)), rhs_vars)
                rhs_ps = filter(x -> isparameter(x), rhs_vars)
                append!(eqs, eq)
                append!(required_states, rhs_vars)
            end
        end
    end

    required_states = unique(value.(required_states))
    states = Any[unique(value.(states))...]
    filter!(x -> !any(isequal(y, x) for y in states), required_states)

    if !isempty(required_states)
        newstateeqs = Equation[]
        for s in required_states
            # Handled based on metadata of each state (for now just one)
            if ismembranecurrent(s) && isaggregator(s)
                push!(newstateeqs, s ~ sum(filter(x -> iontype(x) == iontype(s), currents)))
                push!(states, s)
            end
            # ... other required state handlers ...
        end
        append!(eqs, newstateeqs)
    end

    # propagate default parameter values to channel systems
    vm_eq = D(Vₘ) ~ (Iapp - (+(currents..., Isyn)))/(aₘ*cₘ)
    push!(eqs, vm_eq)
    compartment = ODESystem(eqs, t, states, params; defaults = defaultmap, name = name)
    comp_sys = compose(compartment, channel_systems)
    return Compartment{Sphere}(capacitance, channels, states, params, comp_sys)
end

const Soma = Compartment{Sphere}

function Compartment{Cylinder}() end
# should also be able to parse "collections" of compartments" that have an adjacency list/matrix

# takes a topology, which for now is just an adjacency list; also list of neurons, but we
# should be able to just auto-detect all the neurons in the topology
function Network(neurons, topology; name = Base.gensym(:Network))

    all_neurons = Set(getproperty.(neurons, :sys))
    eqs = Equation[]
    params = Set()
    states = Set() # place holder
    defaultmap = Pair[]
    all_systems = AbstractSystem[]
    push!(all_systems, all_neurons...)
    post_neurons = Set()

    # what types of synapses do we have
    all_synapses = [synapse.second[2] for synapse in topology]
    all_synapse_types = [nameof(x.sys) for x in all_synapses]
    synapse_types = unique(all_synapse_types)

    # how many of each synapse type are there
    synapse_counts = Dict{Symbol,Int64}()

    for n in synapse_types
        c = count(x -> isequal(x, n), all_synapse_types)
        push!(synapse_counts, n => c)
    end

    # Extract reversal potentials
    reversals = unique([x.reversal for x in all_synapses])
    push!(params, reversals...)
    rev_meta = [getmetadata(x, ConductorEquilibriumCtx).val for x in reversals]

    for (i,j) in zip(reversals, rev_meta)
        push!(defaultmap, i => ustrip(Float64, mV, j))
    end

    voltage_fwds = Set()

    # create a unique synaptic current for each post-synaptic target
    for synapse in topology
        pre = synapse.first.sys # pre-synaptic neuron system
        post = synapse.second[1].sys # post-synaptic neuron system
        push!(post_neurons, post)
        Erev = synapse.second[2].reversal
        syntype = synapse.second[2].sys # synapse system
        synname = nameof(syntype)
        syn = @set syntype.name = Symbol(synname, synapse_counts[synname]) # each synapse is uniquely named
        synapse_counts[synname] -= 1
        push!(all_systems, syn)
        push!(voltage_fwds, syn.Vₘ ~ pre.Vₘ)

        if post.Isyn ∈ Set(vcat((get_variables(x.lhs) for x in eqs)...))
            idx = findfirst(x -> isequal([post.Isyn], get_variables(x.lhs)), eqs)
            eq = popat!(eqs, idx)
            eq = eq.lhs ~ eq.rhs + (syn.g * (post.Vₘ - Erev))
            push!(eqs, eq)
        else
            push!(eqs, post.Isyn ~ syn.g * (post.Vₘ - Erev))
        end
    end

    for nonpost in setdiff(all_neurons, post_neurons)
        push!(eqs, D(nonpost.Isyn) ~ 0)
    end

    append!(eqs, collect(voltage_fwds))
    network_system = ODESystem(eqs, t, states, params; defaults = defaultmap, name = name )
    return compose(network_system, all_systems)
end

function Simulation(network; time::Time, system = false)
    t_val = ustrip(Float64, ms, time)
    simplified = structural_simplify(network)
    if system
        return simplified
    else
        @info repr("text/plain", simplified)
        return ODAEProblem(simplified, [], (0., t_val), [])
    end
end

function Simulation(neuron::Soma; time::Time, system = false)
    # for a single neuron, we just need a thin layer to set synaptic current constant
    old_sys = neuron.sys
    Isyn = getproperty(old_sys, :Isyn, namespace=false)
    wrapper = ODESystem([D(Isyn) ~ 0]; name = Base.gensym(:wrapper))
    # Use non-flattening extend
    new_neuron_sys = _extend(wrapper, old_sys; name = nameof(old_sys))
    t_val = ustrip(Float64, ms, time)
    simplified = structural_simplify(new_neuron_sys)
    if system
        return simplified
    else
        @info repr("text/plain", simplified)
        return ODAEProblem(simplified, [], (0., t_val), [])
    end
end

end # module


text:
module MCMCDiagnosticTools

using AbstractFFTs: AbstractFFTs
using DataAPI: DataAPI
using Distributions: Distributions
using MLJModelInterface: MLJModelInterface
using SpecialFunctions: SpecialFunctions
using StatsBase: StatsBase
using Tables: Tables

using LinearAlgebra: LinearAlgebra
using Random: Random
using Statistics: Statistics

export discretediag
export ess_rhat, ESSMethod, FFTESSMethod, BDAESSMethod
export gelmandiag, gelmandiag_multivariate
export gewekediag
export heideldiag
export mcse
export rafterydiag
export rstar

include("discretediag.jl")
include("ess.jl")
include("gelmandiag.jl")
include("gewekediag.jl")
include("heideldiag.jl")
include("mcse.jl")
include("rafterydiag.jl")
include("rstar.jl")
end


text:
module SpringCollab2020FlagToggleWaterfall

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/FlagToggleWaterfall" FlagToggleWaterfall(x::Integer, y::Integer, flag::String="flag_toggle_waterfall", inverted::Bool=false)

const fillColor = Ahorn.XNAColors.LightBlue .* 0.3
const surfaceColor = Ahorn.XNAColors.LightBlue .* 0.8

const waterSegmentMatrix = [
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 0 0 0 0 0 0 1 1 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 1 0 0 0 0 0 0 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 0 0 0 0 0 0 1 1;
    1 1 0 0 0 0 0 0 1 1;
]

const waterSegment = Ahorn.matrixToSurface(
    waterSegmentMatrix,
    [
        fillColor,
        surfaceColor
    ]
)

function getHeight(entity::FlagToggleWaterfall, room::Maple.Room)
    waterEntities = filter(e -> e.name == "water" || e.name == "SpringCollab2020/FlagToggleWater", room.entities)
    waterRects = Ahorn.Rectangle[
        Ahorn.Rectangle(
            Int(get(e.data, "x", 0)),
            Int(get(e.data, "y", 0)),
            Int(get(e.data, "width", 8)),
            Int(get(e.data, "height", 8))
        ) for e in waterEntities
    ]

    width, height = room.size
    x, y = Int(get(entity.data, "x", 0)), Int(get(entity.data, "y", 0))
    tx, ty = floor(Int, x / 8) + 1, floor(Int, y / 8) + 1

    wantedHeight = 8 - y % 8
    while wantedHeight < height - y
        rect = Ahorn.Rectangle(x, y + wantedHeight, 8, 8)

        if any(Ahorn.checkCollision.(waterRects, Ref(rect)))
            break
        end

        if get(room.fgTiles.data, (ty + 1, tx), '0') != '0'
            break
        end

        wantedHeight += 8
        ty += 1
    end

    return wantedHeight
end

const placements = Ahorn.PlacementDict(
    "Flag Toggle Waterfall (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FlagToggleWaterfall
    )
)

function Ahorn.selection(entity::FlagToggleWaterfall, room::Maple.Room)
    x, y = Ahorn.position(entity)
    height = getHeight(entity, room)

    return Ahorn.Rectangle(x, y, size(waterSegmentMatrix, 2), height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleWaterfall, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = getHeight(entity, room)
    segmentHeight, segmentWidth = size(waterSegmentMatrix)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, segmentWidth, height)
    Ahorn.clip(ctx)

    for i in 0:segmentHeight:ceil(Int, height / segmentHeight) * segmentHeight
        Ahorn.drawImage(ctx, waterSegment, 0, i)
    end

    Ahorn.restore(ctx)
end

end


text:
# This file was generated, do not modify it. # hide
using MLJ, DataFrames, Statistics, PrettyPrinting

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.7);


text:
import Dates
using AutoHashEquals

const WEEKDAYS = Dict(map(reverse,enumerate(map(Symbol ∘ uppercase, Dates.ENGLISH.days_of_week_abbr))))
const NTH_PERIODS = ["1st", "2nd", "3rd", "4th", "5th"]
const NTH_LAST_PERIODS = ["Last", "2nd Last", "3rd Last", "4th Last", "5th Last"]
const PERIODS = merge(Dict(map(reverse,enumerate(NTH_PERIODS))), Dict(map(reverse,enumerate(NTH_LAST_PERIODS))))
const MONTHS = Dict(zip(map(Symbol ∘ uppercase, Dates.ENGLISH.months_abbr),range(1,stop=12)))

@auto_hash_equals struct FDOM <: RDate
    calendars::Union{Vector{String}, Nothing}
    FDOM() = new(nothing)
    FDOM(calendars) = new(calendars)
end
function apply(rdate::FDOM, d::Dates.Date, cal_mgr::CalendarManager)
    fdom = Dates.firstdayofmonth(d)
    if rdate.calendars !== nothing
        fdom = apply(HolidayRoundingNBD(), fdom, calendar(cal_mgr, rdate.calendars))
    end

    return fdom
end
multiply_roll(x::FDOM, ::Integer) = x
multiply_no_roll(x::FDOM, ::Integer) = x
Base.:-(x::FDOM) = x

Base.show(io::IO, rdate::FDOM) = print(io, rdate.calendars === nothing ? "FDOM" : "FDOM@$(join(rdate.calendars, "|"))")
register_grammar!(E"FDOM" > FDOM)
register_grammar!(E"FDOM@" + p"[a-zA-Z\\\\\\s\\|]+" > calendar_name -> FDOM(map(String, split(calendar_name, "|"))))

@auto_hash_equals struct LDOM <: RDate
    calendars::Union{Vector{String}, Nothing}
    LDOM() = new(nothing)
    LDOM(calendars) = new(calendars)
end
function apply(rdate::LDOM, d::Dates.Date, cal_mgr::CalendarManager)
    ldom = Dates.lastdayofmonth(d)
    if rdate.calendars !== nothing
        ldom = apply(HolidayRoundingPBD(), ldom, calendar(cal_mgr, rdate.calendars))
    end

    return ldom
end
multiply_roll(x::LDOM, ::Integer) = x
multiply_no_roll(x::LDOM, ::Integer) = x
Base.:-(x::LDOM) = x

Base.show(io::IO, rdate::LDOM) = print(io, rdate.calendars === nothing ? "LDOM" : "LDOM@$(join(rdate.calendars, "|"))")
register_grammar!(E"LDOM" > LDOM)
register_grammar!(E"LDOM@" + p"[a-zA-Z\\\\\\s\\|]+" > calendar_name -> LDOM(map(String, split(calendar_name, "|"))))

struct Easter <: RDate
    yearδ::Int64
end

function apply(rdate::Easter, date::Dates.Date, cal_mgr::CalendarManager)
    y = Dates.year(date) + rdate.yearδ
    a = rem(y, 19)
    b = div(y, 100)
    c = rem(y, 100)
    d = div(b, 4)
    e = rem(b, 4)
    f = div(b + 8, 25)
    g = div(b - f + 1, 3)
    h = rem(19*a + b - d - g + 15, 30)
    i = div(c, 4)
    k = rem(c, 4)
    l = rem(32 + 2*e + 2*i - h - k, 7)
    m = div(a + 11*h + 22*l, 451)
    n = div(h + l - 7*m + 114, 31)
    p = rem(h + l - 7*m + 114, 31)
    return Dates.Date(y, n, p + 1)
end
multiply_roll(x::Easter, count::Integer) = Easter(x.yearδ*count)
multiply_no_roll(x::Easter, count::Integer) = Easter(x.yearδ*count)

Base.:-(rdate::Easter) = Easter(-rdate.yearδ)
Base.show(io::IO, rdate::Easter) = print(io, "$(rdate.yearδ)E")
register_grammar!(PInt64() + E"E" > Easter)

struct Day <: RDate
    days::Int64
end

apply(x::Day, y::Dates.Date, cal_mgr::CalendarManager) = y + Dates.Day(x.days)
multiply_roll(x::Day, count::Integer) = Day(x.days*count)
multiply_no_roll(x::Day, count::Integer) = Day(x.days*count)
Base.:-(x::Day) = Day(-x.days)
Base.:+(x::Day, y::Day) = Day(x.days + y.days)

Base.show(io::IO, rdate::Day) = print(io, "$(rdate.days)d")
register_grammar!(PInt64() + E"d" > Day)

struct Week <: RDate
    weeks::Int64
end

apply(x::Week, y::Dates.Date, cal_mgr::CalendarManager) = y + Dates.Week(x.weeks)
multiply_roll(x::Week, count::Integer) = Week(x.weeks*count)
multiply_no_roll(x::Week, count::Integer) = Week(x.weeks*count)
Base.:-(x::Week) = Week(-x.weeks)
Base.:+(x::Week, y::Week) = Week(x.weeks + y.weeks)
Base.:+(x::Day, y::Week) = Day(x.days + 7*y.weeks)
Base.:+(x::Week, y::Day) = Day(7*x.weeks + y.days)

Base.show(io::IO, rdate::Week) = print(io, "$(rdate.weeks)w")
register_grammar!(PInt64() + E"w" > Week)

struct Month <: RDate
    months::Int64
    idc::InvalidDayConvention
    mic::MonthIncrementConvention

    Month(months::Int64) = new(months, InvalidDayLDOM(), MonthIncrementPDOM())
    Month(months::Int64, idc::InvalidDayConvention, mic::MonthIncrementConvention) = new(months, idc, mic)
end

function apply(rdate::Month, date::Dates.Date, cal_mgr::CalendarManager)
    y, m = Dates.yearmonth(date)
    ny = Dates.yearwrap(y, m, rdate.months)
    nm = Dates.monthwrap(m, rdate.months)
    ay, am, ad = adjust(rdate.mic, date, nm, ny, cal_mgr)
    ld = Dates.daysinmonth(ay, am)
    return ad <= ld ? Dates.Date(ay, am, ad) : adjust(rdate.idc, ad, am, ay)
end
multiply_no_roll(x::Month, count::Integer) = Month(x.months*count, x.idc, x.mic)
Base.:-(x::Month) = Month(-x.months)

Base.show(io::IO, rdate::Month) = (print(io, "$(rdate.months)m["), show(io, rdate.idc), print(io, ";"), show(io, rdate.mic), print(io,"]"))
register_grammar!(PInt64() + E"m" > Month)
register_grammar!(PInt64() + E"m[" + Alt(map(Pattern, collect(keys(INVALID_DAY_MAPPINGS)))...) + E";" + Alt(map(Pattern, collect(keys(MONTH_INCREMENT_MAPPINGS)))...) + E"]" > (d,idc,mic) -> Month(d, INVALID_DAY_MAPPINGS[idc], MONTH_INCREMENT_MAPPINGS[mic]))
# We also have the more complex PDOMEOM month increment which we handle separately.
register_grammar!(PInt64() + E"m[" + Alt(map(Pattern, collect(keys(INVALID_DAY_MAPPINGS)))...) + E";PDOMEOM@" + p"[a-zA-Z\\\\\\s\\|]+" + E"]" > (d,idc,calendar_name) -> Month(d, INVALID_DAY_MAPPINGS[idc], MonthIncrementPDOMEOM(map(String, split(calendar_name, "|")))))

struct Year <: RDate
    years::Int64
    idc::InvalidDayConvention
    mic::MonthIncrementConvention

    Year(years::Int64) = new(years, InvalidDayLDOM(), MonthIncrementPDOM())
    Year(years::Int64, idc::InvalidDayConvention, mic::MonthIncrementConvention) = new(years, idc, mic)
end

function apply(rdate::Year, date::Dates.Date, cal_mgr::CalendarManager)
    oy, m = Dates.yearmonth(date)
    ny = oy + rdate.years
    (ay, am, ad) = adjust(rdate.mic, date, m, ny, cal_mgr)
    ld = Dates.daysinmonth(ay, am)
    return ad <= ld ? Dates.Date(ay, am, ad) : adjust(rdate.idc, ad, am, ay)
end
multiply_no_roll(x::Year, count::Integer) = Year(x.years*count, x.idc, x.mic)
Base.:-(x::Year) = Year(-x.years)

Base.show(io::IO, rdate::Year) = (print(io, "$(rdate.years)y["), show(io, rdate.idc), print(io, ";"), show(io, rdate.mic), print(io,"]"))
register_grammar!(PInt64() + E"y" > Year)
register_grammar!(PInt64() + E"y[" + Alt(map(Pattern, collect(keys(INVALID_DAY_MAPPINGS)))...) + E";" + Alt(map(Pattern, collect(keys(MONTH_INCREMENT_MAPPINGS)))...) + E"]" > (d,idc,mic) -> Year(d, INVALID_DAY_MAPPINGS[idc], MONTH_INCREMENT_MAPPINGS[mic]))
# We also have the more complex PDOMEOM month increment which we handle separately.
register_grammar!(PInt64() + E"y[" + Alt(map(Pattern, collect(keys(INVALID_DAY_MAPPINGS)))...) + E";PDOMEOM@" + p"[a-zA-Z\\\\\\s\\|]+" + E"]" > (d,idc,calendar_name) -> Year(d, INVALID_DAY_MAPPINGS[idc], MonthIncrementPDOMEOM(map(String, split(calendar_name, "|")))))

struct DayMonth <: RDate
    day::Int64
    month::Int64

    DayMonth(day::Int64, month::Int64) = new(day, month)
end

apply(rdate::DayMonth, date::Dates.Date, cal_mgr::CalendarManager) = Dates.Date(Dates.year(date), rdate.month, rdate.day)
multiply_roll(x::DayMonth, ::Integer) = x
multiply_no_roll(x::DayMonth, ::Integer) = x
Base.:-(x::DayMonth) = x

Base.show(io::IO, rdate::DayMonth) = print(io, "$(rdate.day)$(uppercase(Dates.ENGLISH.months_abbr[rdate.month]))")
register_grammar!(PPosInt64() + month_short > (d,m) -> DayMonth(d,MONTHS[Symbol(m)]))

struct Date <: RDate
    date::Dates.Date

    Date(date::Dates.Date) = new(date)
    Date(y::Int64, m::Int64, d::Int64) = new(Dates.Date(y, m, d))
end

multiply_roll(x::Date, ::Integer) = x
multiply_no_roll(x::Date, ::Integer) = x
Base.:-(x::Date) = x

apply(rdate::Date, date::Dates.Date, cal_mgr::CalendarManager) = rdate.date
Base.show(io::IO, rdate::Date) = print(io, "$(Dates.day(rdate.date))$(uppercase(Dates.ENGLISH.months_abbr[Dates.month(rdate.date)]))$(Dates.year(rdate.date))")
register_grammar!(PPosInt64() + month_short + PPosInt64() > (d,m,y) -> Date(Dates.Date(y, MONTHS[Symbol(m)], d)))

struct NthWeekdays <: RDate
    dayofweek::Int64
    period::Int64

    NthWeekdays(dayofweek::Int64, period::Int64) = new(dayofweek, period)
end

function apply(rdate::NthWeekdays, date::Dates.Date, cal_mgr::CalendarManager)
    wd = Dates.dayofweek(date)
    wd1st = mod(wd - mod(Dates.day(date), 7), 7) + 1
    wd1stdiff = wd1st - rdate.dayofweek
    period = wd1stdiff > 0 ? rdate.period : rdate.period - 1
    days = 7*period - wd1stdiff + 1
    return Dates.Date(Dates.year(date), Dates.month(date), days)
end
multiply_roll(x::NthWeekdays, ::Integer) = x
multiply_no_roll(x::NthWeekdays, ::Integer) = x
Base.:-(x::NthWeekdays) = x

Base.show(io::IO, rdate::NthWeekdays) = print(io, "$(NTH_PERIODS[rdate.period]) $(uppercase(Dates.ENGLISH.days_of_week_abbr[rdate.dayofweek]))")

register_grammar!(Alt(map(Pattern, NTH_PERIODS)...) + space + weekday_short > (p,wd) -> NthWeekdays(WEEKDAYS[Symbol(wd)], PERIODS[p]))

struct NthLastWeekdays <: RDate
    dayofweek::Int64
    period::Int64

    NthLastWeekdays(dayofweek::Int64, period::Int64) = new(dayofweek, period)
end

function apply(rdate::NthLastWeekdays, date::Dates.Date, cal_mgr::CalendarManager)
    ldom = LDOM() + date
    ldom_dow = Dates.dayofweek(ldom)
    ldom_dow_diff = ldom_dow - rdate.dayofweek
    period = ldom_dow_diff >= 0 ? rdate.period - 1 : rdate.period
    days_to_sub = 7*period + ldom_dow_diff
    days = Dates.day(ldom) - days_to_sub
    return Dates.Date(Dates.year(date), Dates.month(date), days)
end
multiply_roll(x::NthLastWeekdays, ::Integer) = x
multiply_no_roll(x::NthLastWeekdays, ::Integer) = x
Base.:-(x::NthLastWeekdays) = x

Base.show(io::IO, rdate::NthLastWeekdays) = print(io, "$(NTH_LAST_PERIODS[rdate.period]) $(uppercase(Dates.ENGLISH.days_of_week_abbr[rdate.dayofweek]))")

register_grammar!(Alt(map(Pattern, NTH_LAST_PERIODS)...) + space + weekday_short > (p,wd) -> NthLastWeekdays(WEEKDAYS[Symbol(wd)], PERIODS[p]))

struct Weekdays <: RDate
    dayofweek::Int64
    count::Int64
    inclusive::Bool

    Weekdays(dayofweek::Int64, count::Int64) = new(dayofweek, count, false)
    Weekdays(dayofweek::Int64, count::Int64, inclusive::Bool) = new(dayofweek, count, true)
end

function apply(rdate::Weekdays, date::Dates.Date, cal_mgr::CalendarManager)
    dayδ = Dates.dayofweek(date) - rdate.dayofweek
    weekδ = rdate.count

    if rdate.count < 0 && (dayδ > 0 || (rdate.inclusive && dayδ == 0))
        weekδ += 1
    elseif rdate.count > 0 && (dayδ < 0 || (rdate.inclusive && dayδ == 0))
        weekδ -= 1
    end

    return date + Dates.Day(weekδ*7 - dayδ)
end

multiply_no_roll(x::Weekdays, count::Integer) = Weekdays(x.dayofweek, x.count * count, x.inclusive)
Base.:-(rdate::Weekdays) = Weekdays(rdate.dayofweek, -rdate.count)

Base.show(io::IO, rdate::Weekdays) = print(io, "$(rdate.count)$(uppercase(Dates.ENGLISH.days_of_week_abbr[rdate.dayofweek]))")
register_grammar!(PNonZeroInt64() + weekday_short > (i,wd) -> Weekdays(WEEKDAYS[Symbol(wd)], i))
register_grammar!(PNonZeroInt64() + weekday_short + E"!" > (i,wd) -> Weekdays(WEEKDAYS[Symbol(wd)], i, true))

@auto_hash_equals struct BizDays{S <: HolidayRoundingConvention} <: RDate
    days::Int64
    calendar_names::Vector{String}
    rounding::S

    # I have twisted and turned on what the right default is for holiday start dates.
    # Going with next business day (or previous business day if negated) for 0b and no rounding otherwise.
    # This means that 0b and 1b will give you the same result if you start from a holiday, but I'm comfortable with that.
    BizDays(days::Int64, calendar_names) = BizDays(days, calendar_names, days == 0 ? HolidayRoundingNBD() : HolidayRoundingNR())
    BizDays(days::Int64, calendar_names, rounding::S) where {S <: HolidayRoundingConvention} = new{S}(days, calendar_names, rounding)
end

function apply(rdate::BizDays, date::Dates.Date, cal_mgr::CalendarManager)
    cal = calendar(cal_mgr, rdate.calendar_names)
    date = apply(rdate.rounding, date, cal)

    count = rdate.days
    if rdate.days > 0
        rounding = HolidayRoundingNBD()
        while count > 0
            date = apply(rounding, date + Dates.Day(1), cal)
            count -= 1
        end
    elseif rdate.days < 0
        rounding = HolidayRoundingPBD()
        while count < 0
            date = apply(rounding, date - Dates.Day(1), cal)
            count += 1
        end
    end

    return date
end

function multiply_no_roll(x::BizDays, count::Integer)
    x = count < 0 ? -x : x
    return BizDays(x.count * abs(count), x.calendar_names, x.rounding)
end

Base.:-(x::BizDays{HolidayRoundingNBD}) = BizDays(-x.days, x.calendar_names, HolidayRoundingPBD())
Base.:-(x::BizDays{HolidayRoundingPBD}) = BizDays(-x.days, x.calendar_names, HolidayRoundingNBD())
Base.:-(x::BizDays{HolidayRoundingNBDSM}) = BizDays(-x.days, x.calendar_names, HolidayRoundingPBDSM())
Base.:-(x::BizDays{HolidayRoundingPBDSM}) = BizDays(-x.days, x.calendar_names, HolidayRoundingNBDSM())
Base.:-(x::BizDays{HolidayRoundingNR}) = BizDays(-x.days, x.calendar_names, HolidayRoundingNR())

register_grammar!(PPosZeroInt64() + E"b@" + p"[a-zA-Z\\\\\\s\\|]+" > (days,calendar_name) -> BizDays(days, map(String, split(calendar_name, "|"))))


text:
# In this script, we'll do some benchmarking to determine how to do the dynamics update. Should we do it inplace, using static arrays etc.?
using StaticArrays, BenchmarkTools, Statistics

N = 100 # Number of particles
d = 4   # Number of dimensions
x = [zeros(d) for _ in 1:N]
x2 = [zeros(d) for _ in 1:N]
xs = [@SVector(zeros(d)) for _ in 1:N]
xs2 = [@SVector(zeros(d)) for _ in 1:N]

f = x->x.^2 # Simple dynamics functions
g = x->x^2  # Simple dynamics functions

# below are four different ways of calling the dynamics
function f_inplace!(x2,x,f)
    for i = eachindex(x)
        x2[i] .= f(x[i])
    end
end

function f_inplace_b!(x2,x,f)
    for i = eachindex(x)
        x2[i] .= f.(x[i])
    end
end

function f_update!(x2,x,f)
    for i = eachindex(x)
        x2[i] = f(x[i])
    end
end

function f_update_b!(x2,x,f)
    for i = eachindex(x)
        x2[i] = f.(x[i])
    end
end;

# Now we do some benchmarking. We do `evals=1000` to simulate the performance for 1000 time steps, including potential gc time etc.
b1 = @benchmark f_inplace!($x2,$x,$f) evals=1000
b2 = @benchmark f_inplace_b!($x2,$x,$g) evals=1000
b3 = @benchmark f_update!($x2,$x,$f) evals=1000
b4 = @benchmark f_update!($xs2,$xs,$f) evals=1000
b5 = @benchmark f_update_b!($x2,$x,$g) evals=1000
b6 = @benchmark f_update_b!($xs2,$xs,$g) evals=1000

bs = [b1,b2,b3,b4,b5,b6]
display([ratio(mean(b1),mean(b2)) for b1 in bs, b2 in bs])
display([ratio(memory(b1),memory(b2)) for b1 in bs, b2 in bs])

# We now add the `::F) where F` to force specialization on `f`
function f_inplace!(x2,x,f::F) where F
    for i = eachindex(x)
        x2[i] .= f(x[i])
    end
end

function f_inplace_b!(x2,x,f::F) where F
    for i = eachindex(x)
        x2[i] .= f.(x[i])
    end
end

function f_update!(x2,x,f::F) where F
    for i = eachindex(x)
        x2[i] = f(x[i])
    end
end

function f_update_b!(x2,x,f::F) where F
    for i = eachindex(x)
        x2[i] = f.(x[i])
    end
end;

# Now we do some benchmarking. We do `evals=1000` to simulate the performance for 1000 time steps, including potential gc time etc.
b7 = @benchmark f_inplace!($x2,$x,$f) evals=1000
b8 = @benchmark f_inplace_b!($x2,$x,$g) evals=1000
b9 = @benchmark f_update!($x2,$x,$f) evals=1000
b10 = @benchmark f_update!($xs2,$xs,$f) evals=1000
b11 = @benchmark f_update_b!($x2,$x,$g) evals=1000
b12 = @benchmark f_update_b!($xs2,$xs,$g) evals=1000

bs2 = [b7, b8, b9, b10, b11, b12]
display(ratio.(mean.(bs), mean.(bs2)))
display(ratio.(memory.(bs), memory.(bs2)))
# NOTE: It was absolutely crucial to do the `f::F) where F`. Without this, b6 was much slower than b4, and b2 and b5 were really slow

display([ratio(mean(b1),mean(b2)) for b1 in bs2, b2 in bs2])
display([ratio(memory(b1),memory(b2)) for b1 in bs2, b2 in bs2])
# It seems like a dynamics function which can return a StaticVector (f_update!) is 7 times faster than inplace modification of a regular array (f_inplace_b!)

# Just for fun, let's also try inplace for MVectors
xm = [@MVector(zeros(d)) for _ in 1:N]
xm2 = [@MVector(zeros(d)) for _ in 1:N]
b13 = @benchmark f_inplace_b!($xm2,$xm,$g) evals=1000
bs3 = [b2,b4,b13]
display([ratio(mean(b1),mean(b2)) for b1 in bs3, b2 in bs3])
display([ratio(memory(b1),memory(b2)) for b1 in bs3, b2 in bs3])
# It seems we pay a factor of 2-3 for the MMatrix

# We now consider statically sized regular vectors
xv = [SizedVector{d}(zeros(d)) for _ in 1:N]
xv2 = [SizedVector{d}(zeros(d)) for _ in 1:N]
b14 = @benchmark f_inplace_b!($xv2,$xv,$g) evals=1000
bs4 = [b2,b4,b14]
display([ratio(mean(b1),mean(b2)) for b1 in bs4, b2 in bs4])
display([ratio(memory(b1),memory(b2)) for b1 in bs4, b2 in bs4])
# It seems we pay a factor of 5 for the SizedVector

#src Compile using literateweave("benchstatic.jl", doctype="github")


text:
cat_tuple_types(T1::Type{<:Tuple}, T2::Type{<:Tuple}) = Tuple{T1.types..., T2.types...}
cat_tuple_types(T1::Type{<:NamedTuple}, ::Type{Tuple{}}) = T1
cat_tuple_types(::Type{Tuple{}}, T2::Type{<:NamedTuple}) = T2
cat_tuple_types(T1::Type{<:NamedTuple}, T2::Type{<:NamedTuple}) =
    NamedTuple{(T1.names..., T2.names...), Tuple{T1.types..., T2.types...}}
cat_tuple_types(T1, T2, T3, T...) = cat_tuple_types(cat_tuple_types(T1, T2), T3, T...)


text:
module DOFBases

using Test
using Gridap
using Gridap.Helpers

export DOFBasis
export test_dof_basis
export numlocaldofs
export dof_type
export LagrangianDOFBasis

import Gridap: evaluate
import Gridap: evaluate!
import Base: length

# Interfaces

"""
Abstract DOF basis
"""
abstract type DOFBasis{D,T} end

"""
Returns the dimension of the DOF basis.
"""
length(::DOFBasis)::Int =  @abstractmethod

"""
Compute the DOF values of the given field. The results in written in place
in the last argument.
E is a subtype of Real (typically Float64). E can be computed as
E = dof_type(T)
"""
function evaluate!(
  ::DOFBasis{D,T},::Field{D,T},::AbstractVector{E}) where {D,T,E}
  @abstractmethod
end

"""
Compute the DOF values of the given basis. The results in written in place
in the last argument.
E is a subtype of Real (typically Float64). E can be computed as
E = dof_type(T). The result is a matrix. The rows correspond to
the Basis and the cols to the DOFBasis.
"""
function evaluate!(
  ::DOFBasis{D,T},::Basis{D,T},::AbstractMatrix{E}) where {D,T,E}
  @abstractmethod
end

"""
Returns the type representing a DOF value. Typically Float64,
but can be other types depending on the input T.
"""
function dof_type(::Type{T}) where T
  E = eltype(T)
  @notimplementedif !(E<:Real)
  E
end

"""
Returns the dimension of the DOF basis.
"""
numlocaldofs(b::DOFBasis) = length(b)

function evaluate(b::DOFBasis{D,T},f::Field{D,T}) where {D,T}
  E = dof_type(T)
  n = length(b)
  v = zeros(E,n)
  evaluate!(b,f,v)
  v
end

function evaluate(b::DOFBasis{D,T},f::Basis{D,T}) where {D,T}
  E = dof_type(T)
  n = length(b)
  m = length(f)
  @assert n == m
  v = zeros(E,n,m)
  evaluate!(b,f,v)
  v
end

# Testers

function test_dof_basis(
  dofbasis::DOFBasis{D,T},
  field::Field{D,T},
  basis::Basis{D,T},
  field_dofs::AbstractVector{E},
  basis_dofs::AbstractMatrix{E}) where {D,T,E}

  @test dof_type(T) == E
  @test length(dofbasis) == length(basis)
  @test numlocaldofs(dofbasis) == length(basis)

  a = evaluate(dofbasis,field)
  @test a ≈ field_dofs

  b = evaluate(dofbasis,basis)
  @test b ≈ basis_dofs

end

# Concrete implementations

"""
Lagrangian DOF basis, which consists on evaluating the polynomial basis
(prebasis) in a set of points (nodes)
"""
struct LagrangianDOFBasis{D,T} <: DOFBasis{D,T}
  nodes::Vector{Point{D,Float64}}
  dof_to_node::Vector{Int}
  dof_to_comp::Vector{Int}
  node_and_comp_to_dof::Matrix{Int}
  _cache_field::Vector{T}
  _cache_basis::Matrix{T}
end

function LagrangianDOFBasis(::Type{T}, nodes::Vector{Point{D,Float64}}) where {D,T}

  nnodes = length(nodes)
  ndofs = _num_dofs(T,nodes)
  dof_to_node, dof_to_comp, node_and_comp_to_dof =
    _generate_dof_layout(T,nnodes)
  cache_field = zeros(T,nnodes)
  cache_basis = zeros(T,ndofs,nnodes)

  LagrangianDOFBasis{D,T}(
    nodes,
    dof_to_node,
    dof_to_comp,
    node_and_comp_to_dof,
    cache_field,
    cache_basis)
end

function length(b::LagrangianDOFBasis{D,T}) where {D,T}
  _num_dofs(T,b.nodes)
end

function evaluate!(
  b::LagrangianDOFBasis{D,T},f::Field{D,T},dofs::AbstractVector{E}) where {D,T,E}
  cache = b._cache_field
  evaluate!(f,b.nodes,cache)
  for (node,v) in enumerate(cache)
    for (comp,vi) in enumerate(v)
      dof = b.node_and_comp_to_dof[node,comp]
      dofs[dof] = vi
    end
  end
end

function evaluate!(
  b::LagrangianDOFBasis{D,T},f::Basis{D,T},dofs::AbstractMatrix{E}) where {D,T,E}
  cache = b._cache_basis
  evaluate!(f,b.nodes,cache)
  for node in 1:size(cache,2)
    for i in 1:size(cache,1)
      v = cache[i,node]
      for (comp,vk) in enumerate(v)
        dof = b.node_and_comp_to_dof[node,comp]
        dofs[i,dof] = vk
      end
    end
  end
end

# Helpers

# Node major implementation
function _generate_dof_layout(::Type{T},nnodes::Integer) where T
  ncomps = _length(T)
  ndofs = ncomps*nnodes
  dof_to_comp = zeros(Int,ndofs)
  dof_to_node = zeros(Int,ndofs)
  node_and_comp_to_dof = zeros(Int,nnodes,ncomps)
  v = zero(T)
  for node in 1:nnodes
    for comp in 1:ncomps
      o = nnodes*(comp-1)
      dof = node+o
      dof_to_comp[dof] = comp
      dof_to_node[dof] = node
      node_and_comp_to_dof[node,comp] = dof
    end
  end
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

function _num_dofs(::Type{T}, nodes::Vector) where T
  ncomps = _length(T)
  nnodes = length(nodes)
  ncomps*nnodes
end

_length(::Type{<:Real}) = 1

_length(::Type{T}) where T = length(T)

end # module


text:
# code for working with expectations of nodes

# define an expectation method on Distributions
"Calculate the expected value of a Node x."
E(x) = x
E(x::Distribution) = mean(x)
V(x) = zero(x)
V(x::Distribution) = var(x)
C(x) = zero(x)
C(x::Distribution{Univariate}) = V(x)
C(x::AbstractMvNormal) = cov(x)
C(x::HMM) = cov(x)
C{D <: Distribution}(x::Vector{D}) = diagm(map(V, x))
H(x) = zero(x)
H(x::Distribution) = entropy(x)

# Define functions for nonrandom nodes.
# In each case, a specialized method is already defined for distributions.
Elog(x) = log(x)
Einv(x) = inv(x)
Eloggamma(x) = lgamma(x)
Elogmvgamma(p, x) = logmvgamma(p, x)
Elogbeta(x) = lbeta(x)
ElogB(x) = lB(x)
Elogdet(x) = logdet(x)
Elogdet{D <: Distribution{Univariate}}(x::Array{D}) = prod(Elog(x))
Elogdet(x::Distribution{Univariate}) = Elog(x)

# now make version of all these functions that work on nodes
# by working elementwise
macro make_mapped_version(fn)
    esc(:($fn{D <: Distribution}(n::Array{D}) = map($fn, n)))
end

to_map = [E, V, C, H, Einv, Elog, ElogB, Eloggamma, Elogmvgamma, Elogdet]
for fn in to_map
    @make_mapped_version fn
end


text:
abstract Interpolation
abstract Interpolation2D <: Interpolation

abstract DerivativeApprox
abstract BoundaryCondition

type Spline <: DerivativeApprox end
type Lagrange <: BoundaryCondition end

# general interpolation methods #
# update if value passed in
function update!(interp::Interpolation, idx::Int, val::Float64)
  interp.y_vals[idx] = val

  update!(interp, idx)

  return interp
end

# locate x
function locate(interp::Interpolation, val::Float64)
  if val < interp.x_vals[1]
    return 1
  elseif val >= interp.x_vals[end - 1]
    # return interp.x_vals[end] - interp.x_vals[1] - 2
    return length(interp.x_vals) - 1
  else
    # return findfirst(interp.x_vals .> val) - 1 # need to look at this
    return searchsortedlast(interp.x_vals, val)
  end
end


text:
using RingPolymerArrays: RingPolymerArrays

struct BCBwithTsit5 <: OrdinaryDiffEq.OrdinaryDiffEqAlgorithm end

OrdinaryDiffEq.isfsal(::BCBwithTsit5) = false

mutable struct BCBwithTsit5Cache{uType,rType,vType,rateType,uEltypeNoUnits,T} <: OrdinaryDiffEq.OrdinaryDiffEqMutableCache
    u::uType
    uprev::uType
    tmp::uType
    rtmp::rType
    vtmp::vType
    k::rateType
    cayley::Vector{Matrix{uEltypeNoUnits}}
    tsit5cache::T
end

function OrdinaryDiffEq.alg_cache(::BCBwithTsit5,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,inplace::Val{true}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
    σ = zero(DynamicsUtils.get_quantum_subsystem(u))
    σrate = zero(DynamicsUtils.get_quantum_subsystem(rate_prototype))
    σprev = zero(DynamicsUtils.get_quantum_subsystem(uprev))
    σprev2 = zero(DynamicsUtils.get_quantum_subsystem(uprev2))
    tsit5cache = OrdinaryDiffEq.alg_cache(OrdinaryDiffEq.Tsit5(), σ, σrate, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, σprev, σprev2, f, t, dt, reltol, p, calck, inplace)
    tmp = zero(u)
    rtmp = zero(DynamicsUtils.get_positions(u))
    vtmp = zero(DynamicsUtils.get_velocities(u))
    k = zero(DynamicsUtils.get_positions(rate_prototype))
    cayley = RingPolymers.cayley_propagator(p.beads, dt; half=false)
    BCBwithTsit5Cache(u, uprev, tmp, rtmp, vtmp, k, cayley, tsit5cache)
end

function OrdinaryDiffEq.initialize!(integrator, cache::BCBwithTsit5Cache)

    r = DynamicsUtils.get_positions(integrator.u)
    v = DynamicsUtils.get_velocities(integrator.u)
    σprev = DynamicsUtils.get_quantum_subsystem(integrator.u)
    Calculators.update_electronics!(integrator.p.calculator, r)
    if integrator.p.method isa DynamicsMethods.EhrenfestMethods.AbstractEhrenfest
        DynamicsMethods.EhrenfestMethods.acceleration!(cache.k, v, r, integrator.p, integrator.t, σprev)
    elseif integrator.p.method isa DynamicsMethods.SurfaceHoppingMethods.SurfaceHopping
        DynamicsMethods.SurfaceHoppingMethods.acceleration!(cache.k, v, r, integrator.p, integrator.t, integrator.p.method.state)
    end
end

@muladd function OrdinaryDiffEq.perform_step!(integrator, cache::BCBwithTsit5Cache, repeat_step=false)
    @unpack t, dt, uprev, u, p = integrator
    @unpack k, rtmp, vtmp, cayley = cache

    rprev = DynamicsUtils.get_positions(uprev)
    vprev = DynamicsUtils.get_velocities(uprev)
    σprev = DynamicsUtils.get_quantum_subsystem(uprev)

    copyto!(rtmp, rprev)

    step_B!(vtmp, vprev, dt/2, k)
    RingPolymerArrays.transform_to_normal_modes!(rtmp, p.beads.transformation)
    RingPolymerArrays.transform_to_normal_modes!(vtmp, p.beads.transformation)
    step_C!(vtmp, rtmp, cayley)
    RingPolymerArrays.transform_from_normal_modes!(rtmp, p.beads.transformation)
    RingPolymerArrays.transform_from_normal_modes!(vtmp, p.beads.transformation)

    Calculators.update_electronics!(p.calculator, rtmp)
    if p.method isa DynamicsMethods.EhrenfestMethods.AbstractEhrenfest
        DynamicsMethods.EhrenfestMethods.acceleration!(k, vtmp, rtmp, p, t, σprev)
    elseif p.method isa DynamicsMethods.SurfaceHoppingMethods.SurfaceHopping
        DynamicsMethods.SurfaceHoppingMethods.acceleration!(k, vtmp, rtmp, p, t, p.method.state)
    end
    step_B!(DynamicsUtils.get_velocities(u), vtmp, dt/2, k)

    copyto!(DynamicsUtils.get_positions(u), rtmp)

    # Tsit5 algorithm modified for only quantum subsystem
    @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76 = cache.tsit5cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,tmp = cache.tsit5cache

    DynamicsUtils.set_quantum_derivative!(k1, vtmp, σprev, p)
    a = dt*a21
    @.. tmp = σprev+a*k1
    DynamicsUtils.set_quantum_derivative!(k2, vtmp, tmp, p)
    @.. tmp = σprev+dt*(a31*k1+a32*k2)
    DynamicsUtils.set_quantum_derivative!(k3, vtmp, tmp, p)
    @.. tmp = σprev+dt*(a41*k1+a42*k2+a43*k3)
    DynamicsUtils.set_quantum_derivative!(k4, vtmp, tmp, p)
    @.. tmp = σprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
    DynamicsUtils.set_quantum_derivative!(k5, vtmp, tmp, p)
    @.. tmp = σprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
    DynamicsUtils.set_quantum_derivative!(k6, vtmp, tmp, p)
    @.. tmp = σprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
    copyto!(DynamicsUtils.get_quantum_subsystem(u), tmp)

    integrator.destats.nf += 6

end

DynamicsMethods.select_algorithm(::RingPolymerSimulation{<:DynamicsMethods.SurfaceHoppingMethods.SurfaceHopping}) = BCBwithTsit5()
DynamicsMethods.select_algorithm(::RingPolymerSimulation{<:DynamicsMethods.EhrenfestMethods.AbstractEhrenfest}) = BCBwithTsit5()


text:
__precompile__()
module FunctionalUtilities

import Base.==

include("FunctionalCore.jl")

include("BitmappedVectorTrie.jl")

include("PersistentVector.jl")

const PVector = PersistentVector
const pvec = PersistentVector

export PersistentVector, PVector, pvec,
       append, push, assoc, peek, pop

include("PersistentDict.jl")

const PDict = PersistentAutoDict

export PersistentArrayDict, PersistentHashDict, PersistentAutoDict, PDict,
       assoc, dissoc

include("PersistentSet.jl")

const PSet = PersistentSet

export PersistentSet, PSet,
       disj

include("PersistentList.jl")

const PList = PersistentList

export PersistentList, PList, EmptyList, cons, .., head, tail

include("PersistentQueue.jl")

const PQueue = PersistentQueue

export PersistentQueue, PQueue, queue, enq

include("ImmutableRef.jl")

export ImmutableRef, Promise, Delay, deliver, isrealized

include("LazyDict.jl")

const LDict = LazyDict

export LazyDict, LDict, islazy

export @Persistent

fromexpr(ex::Expr, ::Type{pvec}) = :(pvec($(esc(ex))))
fromexpr(ex::Expr, ::Type{pset}) = :(pset($(map(esc, ex.args[2:end])...)))
function fromexpr(ex::Expr, ::Type{phmap})
    kvtuples = [:($(esc(kv.args[end-1])), $(esc(kv.args[end])))
                for kv in ex.args[2:end]]
    :(phmap($(kvtuples...)))
end

using Base.Meta: isexpr
macro Persistent(ex)
    if isexpr(ex, [:vcat, :vect])
        fromexpr(ex, pvec)
    elseif isexpr(ex, :call) && ex.args[1] === :Set
        fromexpr(ex, pset)
    elseif isexpr(ex, :call) && ex.args[1] === :Dict
        fromexpr(ex, phmap)
    else
        error("Unsupported @Persistent syntax")
    end
end

################################################################################
# Definitions for freeze, deepfreeze_internal, thaw, and deepthaw_internal

# #TODO These are wrong; actively under development...

#deepfreeze_internal
deepfreeze_internal(x::Core.SimpleVector, stackdict::IdDict) = begin
    haskey(stackdict, x) && return stackdict[x]
    y = PVector(Any[deepfreeze_internal(x[i], stackdict) for i = 1:length(x)]...)
    stackdict[x] = y
    return y
end
deepfreeze_internal(x::Array, stackdict::IdDict) = begin
    haskey(stackdict, x) && return stackdict[x]
    _deepfreeze_array_t(x, eltype(x), stackdict)
end
_deepfreeze_array_t(@nospecialize(x), T, stackdict::IdDict) = begin
    isbitstype(T) && return (stackdict[x]=freeze(x))
    y = PVector(Any[deepfreeze_internal(x[i], stackdict) for i = 1:length(x)]...)

    dest = similar(x)
    stackdict[x] = dest
    for i = 1:(length(x)::Int)
        if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
            xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
            if !isbits(xi)
                xi = deepfreeze_internal(xi, stackdict)
            end
            ccall(:jl_arrayset, Cvoid, (Any, Any, Csize_t), dest, xi, i-1)
        end
    end
    return dest
end
function deepfreeze_internal(x::Union{Dict,IdDict}, stackdict::IdDict)
    if haskey(stackdict, x)
        return stackdict[x]::typeof(x)
    end

    if isbitstype(eltype(x))
        return (stackdict[x] = copy(x))
    end

    dest = empty(x)
    stackdict[x] = dest
    for (k, v) in x
        dest[deepfreeze_internal(k, stackdict)] = deepfreeze_internal(v, stackdict)
    end
    dest
end

end # module FunctionalUtilities

max_issues_repo_name: noahbenson/FunctionalUtilities

max_forks_repo_name: noahbenson/FunctionalUtilities


text:
function solve_lp()
    glpk_optimizer =
        optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)
    model = Model(optimizer_with_attributes(
        CS.Optimizer,
        "lp_optimizer" => glpk_optimizer,
        "logging" => [],
        "seed" => 1,
    ))

    # Variables
    @variable(model, inclusion[h = 1:3], Bin)
    @variable(model, 0 <= allocations[h = 1:3, a = 1:3] <= 1, Int)
    @variable(model, 0 <= days[h = 1:3, a = 1:3] <= 5, Int)

    # Constraints
    @constraint(
        model,
        must_include[h = 1:3],
        sum(allocations[h, a] for a in 1:3) <= inclusion[h]
    )
    # at least n
    @constraint(model, min_hospitals, sum(inclusion[h] for h in 1:3) >= 3)
    # every h must be allocated at most one a
    @constraint(model, must_visit[h = 1:3], sum(allocations[h, a] for a in 1:3) <= 1)
    # every allocated h must have fewer than 5 days of visits per week
    @constraint(
        model,
        max_visits[h = 1:3],
        sum(days[h, a] for a in 1:3) <= 5 * inclusion[h]
    )

    @objective(model, Max, sum(days[h, a] * 5 for h in 1:3, a in 1:3))
    optimize!(model)
end


text:
using Test

@testset "PProf.jl" begin
    include("PProf.jl")
end


text:
using LazyFym
using Transducers

using Random
using StructArrays

## Policy
abstract type Policy end
struct ExamplePolicy <: Policy
end
# command
command(policy::Policy, x, t) = -sum(x)

## simulation envs
struct IntegralReward <: Fym
end
struct RLEnv <: Fym
    env::LazyFym.InputAffineQuadraticCostEnv
    ∫r::IntegralReward
    policy::Policy
end
# dynamics
function ẋ(rlenv::RLEnv, x_nt, t; action=0.0)
    x = x_nt.env
    ẋ = LazyFym.ẋ(rlenv.env, x, t, action)
    r = LazyFym.r(rlenv.env, x, action)
    (; env = ẋ, ∫r = r)
end
# update (RL agent generates zero-order-hold (ZOH) action in this scenario)
function update(rlenv::RLEnv, ẋ, x_nt, t, Δt)
    action = command(rlenv.policy, x_nt.env, t)
    _datum = Dict(:x => x_nt.env, :t => t, :action => action)
    x_next_nt = ∫(rlenv, ẋ, x_nt, t, Δt; action=action)  # `action` will be a kwarg of ẋ
    _datum[:x_next] = x_next_nt.env
    _datum[:Δr] = x_next_nt.∫r - x_nt.∫r
    datum = (; zip(keys(_datum), values(_datum))...)
    return datum, x_next_nt
end
# initial condition
function initial_condition(rlenv::RLEnv)
    (; env=rand(2), ∫r = 0.0)
end
# LazyFym.initial_condition(env::LazyFym.InputAffineQuadraticCostEnv) = rand(2)  # NOTE: `LazyFym.initial_condition` will automatically guess the appropriate initial condition for nested envs.

function main()
    Random.seed!(1)
    rlenv = RLEnv(LazyFym.InputAffineQuadraticCostEnv(), IntegralReward(), ExamplePolicy())
    t0 = 0.0
    tf = 5.00
    Δt = 0.01
    ts = t0:Δt:tf
    x0 = initial_condition(rlenv)
    # x0 = LazyFym.initial_condition(rlenv)
    traj(x0) = Sim(rlenv, x0, ts, ẋ, update) |> collect |> StructArray
    data = traj(x0)
end