text:
module DecFP

using DecFP_jll

import Printf, Random, SpecialFunctions

export Dec32, Dec64, Dec128, @d_str, @d32_str, @d64_str, @d128_str, exponent10, ldexp10, sigexp

const _buffer = Vector{Vector{UInt8}}()

import Base.promote_rule, Base.RefArray

const flags = [0x00000000]

# clear exception flags and return x
function nox(x)
    flags[Threads.threadid()] = 0
    return x
end

# Check exception flags in mask & throw, otherwise returning x;
# always clearing exceptions.  This is a macros so that
# the error message is only evaluated if an exception occurs.
macro xchk(x, exc, args...)
    mask=0x3f
    if !isempty(args) && Meta.isexpr(args[end], :(=)) && args[end].args[1] == :mask # mask=... keyword at end
        mask = esc(args[end].args[2])
        args = args[1:end-1]
    end
    quote
        ret = $(esc(x))
        if $exc === nothing
            flags[Threads.threadid()] = 0
        else
            f = flags[Threads.threadid()]
            flags[Threads.threadid()] = 0
            f & $mask != 0 && throw($exc($(map(esc,args)...)))
        end
        ret
    end
end

#############################################################################

@enum DecFPRoundingMode begin
    DecFPRoundNearest
    DecFPRoundDown
    DecFPRoundUp
    DecFPRoundToZero
    DecFPRoundFromZero
end

Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:Nearest})  = DecFPRoundNearest
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:Down})     = DecFPRoundDown
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:Up})       = DecFPRoundUp
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:ToZero})   = DecFPRoundToZero
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:FromZero}) = DecFPRoundFromZero

function Base.convert(::Type{RoundingMode}, r::DecFPRoundingMode)
    if r == DecFPRoundNearest
        return RoundNearest
    elseif r == DecFPRoundDown
        return RoundDown
    elseif r == DecFPRoundUp
        return RoundUp
    elseif r == DecFPRoundToZero
        return RoundToZero
    elseif r == DecFPRoundFromZero
        return RoundFromZero
    else
        throw(ArgumentError("invalid DecFP rounding mode code: $c"))
    end
end

const roundingmode = [DecFPRoundNearest]

# global vectors must be initialized at runtime (via __init__)
function __init__()
    resize!(_buffer, Threads.nthreads())
    resize!(flags, Threads.nthreads())
    resize!(roundingmode, Threads.nthreads())
    for i = 1:Threads.nthreads()
        _buffer[i] = fill(0x00, 1024)
        flags[i] = 0x00000000
        roundingmode[i] = DecFPRoundNearest
    end
end

# status flags from bid_functions.h:
const INVALID    = 0x01
const UNNORMAL   = 0x02
const DIVBYZERO  = 0x04
const OVERFLOW   = 0x08
const UNDERFLOW  = 0x10
const INEXACT    = 0x20

bidsym(w,s...) = string("__bid", w, "_", s...)

abstract type DecimalFloatingPoint <: AbstractFloat end

Base.Rounding.rounding_raw(::Type{T}) where {T<:DecimalFloatingPoint} =
    roundingmode[Threads.threadid()]
Base.Rounding.setrounding_raw(::Type{T}, r::DecFPRoundingMode) where {T<:DecimalFloatingPoint} =
    roundingmode[Threads.threadid()] = r

Base.Rounding.rounding(::Type{T}) where {T<:DecimalFloatingPoint} =
    convert(RoundingMode, Base.Rounding.rounding_raw(T))
Base.Rounding.setrounding(::Type{T}, r::RoundingMode) where {T<:DecimalFloatingPoint} =
    Base.Rounding.setrounding_raw(T, convert(DecFPRoundingMode, r))

for w in (32,64,128)
    BID = Symbol(string("Dec",w))
    Ti = Symbol(string("UInt",w))
    @eval struct $BID <: DecimalFloatingPoint
        x::$Ti
        $BID(x::Number) = convert($BID, x)
        Base.reinterpret(::Type{$BID}, x::$Ti) = new(x)
    end

    @eval function $BID(x::Real, mode::RoundingMode)
        setrounding($BID, mode) do
            convert($BID, x)
        end
    end

    @eval function $BID(sign::Integer, significand::Integer, exponent::Integer)
        -1 <= sign <= 1 || throw(DomainError(sign, "sign must be -1, 0, or +1"))
        significand == 0 && return flipsign(zero($BID), sign)
        sign == 0 && throw(DomainError(sign, "sign must be -1 or +1 for non-zero significand"))
        p = 9 * $w ÷ 32 - 2
        emax = 3 * 2^($w ÷ 16 + 3)
        bias = emax + p - 2
        bemax = emax * 2 - 1
        ebits = $w ÷ 16 + 6
        sbits = $w - ebits - 1
        sb = signbit(sign) ? one($Ti) << ($w - 1) : zero($Ti)
        e = exponent + bias
        s = $Ti(abs(significand))
        while s >= $Ti(maxintfloat($BID))
            q, r = divrem(s, $Ti(10))
            r != 0 && throw(InexactError(Symbol($BID), $BID, (sign, significand, exponent)))
            s = q
            e += 1
        end
        while e > bemax && s < $Ti(maxintfloat($BID)) / 10
            s *= $Ti(10)
            e -= 1
        end
        e > bemax && throw(InexactError(Symbol($BID), $BID, (sign, significand, exponent)))
        while e < 0
            q, r = divrem(s, $Ti(10))
            r != 0 && throw(InexactError(Symbol($BID), $BID, (sign, significand, exponent)))
            s = q
            e += 1
        end
        if s < one($Ti) << sbits
            return reinterpret($BID, sb | ($Ti(e) << sbits) | s)
        end
        return reinterpret($BID, sb | ($Ti(0x3) << ($w - 3)) | ($Ti(e) << (sbits - 2)) | (s & (typemax($Ti) >> (ebits + 3))))
    end

    @eval $BID(significand::Integer, exponent::Integer) = $BID(sign(significand), significand, exponent)

    @eval @doc """
        $($BID)(x::Union{Real, AbstractString} [, mode::RoundingMode])
        $($BID)([sign::Integer,] significand::Integer, exponent::Integer)

    Create a $($w)-bit IEEE 754-2008 decimal floating point number. The `mode` argument
    specifies the direction in which the result should be rounded if the conversion cannot
    be done exactly. If not provided, the `mode` is set by the current
    `rounding(DecFP.DecimalFloatingPoint)` mode.

    `$($BID)(x::Real)` is the same as `convert($($BID), x)`.

    `$($BID)(x::AbstractString)` is the same as `parse($($BID), x)`. This is provided for
    convenience since decimal literals are converted to `Float64` when parsed and may not
    produce what you expect.

    `$($BID)(sign, significand, exponent)` returns `sign * significand * 10^exponent`.
    If `sign` isn't passed, use the sign of `significand`.

    # Examples
    ```julia-repl
    julia> $($BID)(1)
    1.0

    julia> $($BID)(1.5)
    1.5

    julia> $($BID)("0.99999999999999999999999999999999999")
    1.0

    julia> $($BID)("0.99999999999999999999999999999999999", RoundDown)
    0.$("9"^(9*$w÷32-2))

    julia> $($BID)(-1, 123456, -4)
    -12.3456
    ```
    """ $BID

    # fix method ambiguities:
    @eval $BID(x::Rational{T}) where {T} = convert($BID, x)
end

# quickly check whether s begins with "±nan"
function isnanstr(s::AbstractString)
    st = iterate(s)
    c, i = '\0', 0
    while st !== nothing
        c, i = st
        isspace(c) || break
        st = iterate(s, i)
    end
    st === nothing && return false
    if c == '+' || c == '-'
        st = iterate(s, i)
        st === nothing && return false
        c, i = st
    end
    (c == 'n' || c == 'N') || return false
    st = iterate(s, i)
    if st !== nothing
        c, i = st
        if c == 'a' || c == 'A'
            st = iterate(s, i)
            if st !== nothing
                c, i = st
                if c == 'n' || c == 'N'
                    st = iterate(s, i)
                    if st === nothing
                        return true
                    end
                end
            end
        end
    end
    return false
end

"""
    sigexp(x::DecFP.DecimalFloatingPoint)

Return `(sign, significand, exponent)` such that `x` is equal to `sign * significand * 10^exponent`.
Throws `DomainError` for infinite or `NaN` arguments.

# Examples
```jldoctest
julia> sigexp(Dec64(1.25))
(1, 125, -2)
```
"""
sigexp(x::DecimalFloatingPoint)

"""
    exponent10(x::DecFP.DecimalFloatingPoint)

Get the exponent of the base 10 representation of a normalized floating-point number.

# Examples
```jldoctest
julia> exponent10(Dec64(123))
2
```
"""
exponent10(x::DecimalFloatingPoint)

"""
    ldexp10(x::DecFP.DecimalFloatingPoint, n::Integer)

Compute ``x * 10^n``.

# Examples
```jldoctest
julia> ldexp10(Dec64(15), 2)
1500.0
```
"""
ldexp10(x::DecFP.DecimalFloatingPoint, n::Integer)

for w in (32,64,128)
    BID = Symbol(string("Dec",w))
    Ti = eval(Symbol(string("UInt",w)))
    Tsi = eval(Symbol(string("Int",w)))
    T = eval(BID)

    # hack: we need an internal parsing function that doesn't check exceptions, since
    # flags isn't defined until __init__ runs.  Similarly for nextfloat/prevfloat
    @eval begin
        _parse(::Type{$BID}, s::AbstractString) =
            ccall(($(bidsym(w,"from_string")), libbid), $BID, (Ptr{UInt8},Cuint,Ref{Cuint}), s, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid()))
        _nextfloat(x::$BID) = ccall(($(bidsym(w,"nexttoward")), libbid), $BID, ($BID,Dec128,Ref{Cuint}), x, pinf128, RefArray(flags, Threads.threadid()))
        _prevfloat(x::$BID) = ccall(($(bidsym(w,"nexttoward")), libbid), $BID, ($BID,Dec128,Ref{Cuint}), x, minf128, RefArray(flags, Threads.threadid()))
        _sub(x::$BID, y::$BID) = ccall(($(bidsym(w,"sub")), libbid), $BID, ($BID,$BID,Cuint,Ref{Cuint}), x, y, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid()))
    end

    @eval begin
        function Base.parse(::Type{$BID}, s::AbstractString)
            x = _parse($BID, s)
            if isnan(x) && !isnanstr(s)
                throw(ArgumentError("invalid number format $s"))
            end
            return @xchk(x, nothing)
        end

        function Base.tryparse(::Type{$BID}, s::AbstractString)
            x = _parse($BID, s)
            if isnan(x) && !isnanstr(s)
                @xchk(x, nothing)
                return nothing
            end
            return @xchk(x, nothing)
        end

        $BID(x::AbstractString) = parse($BID, x)

        function $BID(x::AbstractString, mode::RoundingMode)
            setrounding($BID, mode) do
                parse($BID, x)
            end
        end

        function tostring(x::$BID)
            # fills global _buffer
            ccall(($(bidsym(w,"to_string")), libbid), Cvoid, (Ptr{UInt8},$BID,Ref{Cuint}), _buffer[Threads.threadid()], x, RefArray(flags, Threads.threadid()))
        end

        function Base.show(io::IO, x::$BID)
            isnan(x) && (print(io, "NaN"); return)
            isinf(x) && (print(io, signbit(x) ? "-Inf" : "Inf"); return)
            x == 0 && (print(io, signbit(x) ? "-0.0" : "0.0"); return)
            tostring(x)
            buffer = _buffer[Threads.threadid()]
            if buffer[1] == UInt8('-')
                print(io, '-')
            end
            normalized_exponent = exponent10(x)
            lastdigitindex = findfirst(isequal(UInt8('E')), buffer) - 1
            lastnonzeroindex = findlast(!isequal(UInt8('0')), view(buffer, 1:lastdigitindex))
            if -5 < normalized_exponent < 6
                # %f
                if normalized_exponent >= 0
                    if normalized_exponent >= lastnonzeroindex - 2
                        GC.@preserve buffer unsafe_write(io, pointer(buffer, 2), lastnonzeroindex - 1)
                        printzeros(io, normalized_exponent - lastnonzeroindex + 2)
                        print(io, ".0")
                    else
                        GC.@preserve buffer unsafe_write(io, pointer(buffer, 2), normalized_exponent + 1)
                        print(io, '.')
                        GC.@preserve buffer unsafe_write(io, pointer(buffer, normalized_exponent + 3), lastnonzeroindex - normalized_exponent - 2)
                    end
                else
                    print(io, "0.")
                    printzeros(io, -normalized_exponent - 1)
                    GC.@preserve buffer unsafe_write(io, pointer(buffer, 2), lastnonzeroindex - 1)
                end
            else
                # %e
                print(io, Char(buffer[2]), '.')
                if lastnonzeroindex == 2
                    print(io, '0')
                else
                    GC.@preserve buffer unsafe_write(io, pointer(buffer, 3), lastnonzeroindex - 2)
                end
                print(io, 'e')
                if normalized_exponent < 0
                    print(io, '-')
                    normalized_exponent = -normalized_exponent
                end
                b_lb = div(normalized_exponent, 10)
                b = 1
                while b <= b_lb
                    b *= 10
                end
                r = normalized_exponent
                while b > 0
                    q, r = divrem(r, b)
                    print(io, '0' + q)
                    b = div(b, 10)
                end
            end
            return
        end

        function Printf.Printf.fix_dec(x::$BID, n::Int, digits)
            if n > length(digits) - 1
                n = length(digits) - 1
            end
            rounded = round(ldexp10(x, n), RoundNearestTiesAway)
            if rounded == 0
                digits[1] = UInt8('0')
                return Int32(1), Int32(1), signbit(x)
            end
            tostring(rounded)
            buffer = _buffer[Threads.threadid()]
            trailing_zeros = 0
            i = 2
            while buffer[i] != UInt8('E')
                digits[i - 1] = buffer[i]
                if buffer[i] == UInt8('0')
                    trailing_zeros += 1
                else
                    trailing_zeros = 0
                end
                i += 1
            end
            ndigits = i - 2
            len = ndigits - trailing_zeros
            i += 1
            if buffer[i] == UInt8('+')
                expsign = +1
            elseif buffer[i] == UInt8('-')
                expsign = -1
            end
            exponent = 0
            i += 1
            while buffer[i] != 0x00
                exponent = exponent * 10 + buffer[i] - UInt8('0')
                i += 1
            end
            exponent *= expsign
            pt = ndigits + exponent - n
            neg = signbit(x)
            return Int32(len), Int32(pt), neg
        end

        function Printf.Printf.ini_dec(x::$BID, n::Int, digits)
            if n > length(digits) - 1
                n = length(digits) - 1
            end
            if x == 0
                for i = 1:n
                    digits[i] = UInt8('0')
                end
                return Int32(1), Int32(1), signbit(x)
            end
            normalized_exponent = exponent10(x)
            rounded = round(ldexp10(x, n - 1 - normalized_exponent), RoundNearestTiesAway)
            rounded_exponent = exponent10(rounded)
            tostring(rounded)
            buffer = _buffer[Threads.threadid()]
            i = 2
            while buffer[i] != UInt8('E')
                digits[i - 1] = buffer[i]
                i += 1
            end
            while i <= n + 1
                digits[i - 1] = UInt8('0')
                i += 1
            end
            pt = normalized_exponent + rounded_exponent - n + 2
            neg = signbit(x)
            return Int32(n), Int32(pt), neg
        end

        function sigexp(x::$BID)
            isnan(x) && throw(DomainError(x, "sigexp(x) is not defined for NaN."))
            isinf(x) && throw(DomainError(x, "sigexp(x) is only defined for finite x."))
            p = 9 * $w ÷ 32 - 2
            emax = 3 * 2^($w ÷ 16 + 3)
            bias = emax + p - 2
            ebits = $w ÷ 16 + 6
            sbits = $w - ebits - 1
            n = reinterpret($Ti, x)
            if n & $Ti(0x3) << ($w - 3) == $Ti(0x3) << ($w - 3)
                s = ((n & typemax($Ti) >> (ebits + 3)) | one($Ti) << sbits) % $Tsi
                e = ((n & typemax($Ti) << (sbits + 1) >> 3) >> (sbits - 2)) % Int
            else
                s = (n & typemax($Ti) >> (ebits + 1)) % $Tsi
                e = ((n & typemax($Ti) << (sbits + 1) >> 1) >> sbits) % Int
            end
            e -= bias
            return signbit(x) ? -1 : 1, s, e
        end

        Base.fma(x::$BID, y::$BID, z::$BID) = nox(ccall(($(bidsym(w,"fma")), libbid), $BID, ($BID,$BID,$BID,Cuint,Ref{Cuint}), x, y, z, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
        Base.muladd(x::$BID, y::$BID, z::$BID) = fma(x,y,z) # faster than x+y*z

        Base.one(::Union{Type{$BID},$BID}) = $(_parse(T, "1"))
        Base.zero(::Union{Type{$BID},$BID}) = $(_parse(T, "0"))

        Base.signbit(x::$BID) = $(zero(Ti)) != $(Ti(1) << (Ti(w - 1))) & x.x
        Base.sign(x::$BID) = ifelse(isnan(x) || iszero(x), x, ifelse(signbit(x), $(_parse(T, "-1")), $(_parse(T, "1"))))

        Base.nextfloat(x::$BID) = nox(_nextfloat(x))
        Base.prevfloat(x::$BID) = nox(_prevfloat(x))
        Base.eps(x::$BID) = ifelse(isfinite(x), @xchk(nextfloat(x) - x, OverflowError, "$($BID) value overflow", mask=OVERFLOW), $(_parse(T, "NaN")))

        # the meaning of the exponent is different than for binary FP: it is 10^n, not 2^n:
        exponent10(x::$BID) = nox(ccall(($(bidsym(w,"ilogb")), libbid), Cint, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())))
        ldexp10(x::$BID, n::Integer) = nox(ccall(($(bidsym(w,"ldexp")), libbid), $BID, ($BID,Cint,Cuint,Ref{Cuint}), x, n, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
    end

    for (f,c) in ((:isnan,"isNaN"), (:isinf,"isInf"), (:isfinite,"isFinite"), (:issubnormal,"isSubnormal"))
        @eval Base.$f(x::$BID) = ccall(($(bidsym(w,c)), libbid), Cint, ($BID,), x) != 0
    end

    for (f,c) in ((:+,"add"), (:-,"sub"), (:*,"mul"), (:/, "div"), (:hypot,"hypot"), (:atan,"atan2"), (:^,"pow"))
        @eval Base.$f(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,$BID,Cuint,Ref{Cuint}), x, y, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
    end
    @eval Base.copysign(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,"copySign")), libbid), $BID, ($BID,$BID,Ref{Cuint}), x, y, RefArray(flags, Threads.threadid())))

    for f in (:exp,:log,:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh,:atanh,:log1p,:expm1,:log10,:log2,:exp2,:exp10,:sqrt,:cbrt)
        @eval Base.$f(x::$BID) = @xchk(ccall(($(bidsym(w,f)), libbid), $BID, ($BID,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)
    end
    @eval Base.abs(x::$BID) = @xchk(ccall(($(bidsym(w,"abs")), libbid), $BID, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)

    @eval Base.rem(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,"fmod")), libbid), $BID, ($BID,$BID,Ref{Cuint}), x, y, RefArray(flags, Threads.threadid())))
    @eval Base.rem(x::$BID, y::$BID, ::RoundingMode{:Nearest}) = nox(ccall(($(bidsym(w,"rem")), libbid), $BID, ($BID,$BID,Ref{Cuint}), x, y, RefArray(flags, Threads.threadid())))

    @eval begin
        function Base.modf(x::$BID)
            ipart = Ref{$BID}()
            fpart = nox(ccall(($(bidsym(w,"modf")), libbid), $BID, ($BID,Ref{$BID},Ref{Cuint}), x, ipart, RefArray(flags, Threads.threadid())))
            fpart, ipart[]
        end
    end

    for (f,c) in ((:trunc,"round_integral_zero"), (:floor,"round_integral_negative"), (:ceil,"round_integral_positive"))
        @eval Base.$f(x::$BID) = @xchk(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)
    end
    @eval Base.:-(x::$BID) = @xchk(ccall(($(bidsym(w,"negate")), libbid), $BID, ($BID,), x), DomainError, x, mask=INVALID)
    @eval Base.round(x::$BID) = @xchk(ccall(($(bidsym(w,"nearbyint")), libbid), $BID, ($BID,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)

    @eval function SpecialFunctions.logabsgamma(x::$BID)
        isequal(modf(x)[1], -zero(x)) && return typemax(x), one(Int32)
        signgam = signbit(x) && mod(x, 2) > 1 ? -one(Int32) : one(Int32)
        y = @xchk(ccall(($(bidsym(w,:lgamma)), libbid), $BID, ($BID,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)
        return y, signgam
    end

    @eval SpecialFunctions.gamma(x::$BID) = @xchk(ccall(($(bidsym(w,:tgamma)), libbid), $BID, ($BID,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)

    @eval Random.rand(r::Random.AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{$BID}}) =  $BID(1, rand(r, zero($Ti):$Ti(maxintfloat($BID)) - one($Ti)), -(9 * $w ÷ 32 - 2))

    for (r,c) in ((RoundingMode{:Nearest},"round_integral_nearest_even"), (RoundingMode{:NearestTiesAway},"round_integral_nearest_away"), (RoundingMode{:ToZero},"round_integral_zero"), (RoundingMode{:Up},"round_integral_positive"), (RoundingMode{:Down},"round_integral_negative"))
        @eval Base.round(x::$BID, ::$r) = @xchk(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), DomainError, x, mask=INVALID)
    end

    for (f,c) in ((:(==),"quiet_equal"), (:>,"quiet_greater"), (:<,"quiet_less"), (:(>=), "quiet_greater_equal"), (:(<=), "quiet_less_equal"))
        @eval Base.$f(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,c)), libbid), Cint, ($BID,$BID,Ref{Cuint}), x, y, RefArray(flags, Threads.threadid())) != 0)
    end

    for Tf in (Float32,Float64)
        bT = string("binary",sizeof(Tf)*8)
        @eval begin
            Base.convert(::Type{$Tf}, x::$BID) = nox(ccall(($(bidsym(w,"to_",bT)), libbid), $Tf, ($BID,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
            Base.$(Symbol("$Tf"))(x::$BID) = convert($Tf, x)
            Base.convert(::Type{$BID}, x::$Tf) = nox(ccall(($(string("__",bT,"_to_","bid",w)), libbid), $BID, ($Tf,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
        end
    end

    for c in (:π, :e, :ℯ, :γ, :catalan, :φ)
        @eval begin
            Base.convert(::Type{$BID}, ::Irrational{$(QuoteNode(c))}) = $(_parse(T, setprecision(256) do
                                                                                string(BigFloat(getfield(MathConstants, c)))
                                                                                end))
        end
    end

    @eval promote_rule(::Type{$BID}, ::Type{Irrational}) = $BID

    for w′ in (32,64,128)
        BID′ = Symbol(string("Dec",w′))
        if w > w′
            @eval promote_rule(::Type{$BID}, ::Type{$BID′}) = $BID
        end
        if w > w′
            @eval Base.convert(::Type{$BID}, x::$BID′) = @xchk(ccall(($(string("__bid",w′,"_to_","bid",w)), libbid), $BID, ($BID′,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), nothing)
        elseif w < w′
            @eval Base.convert(::Type{$BID}, x::$BID′) = @xchk(ccall(($(string("__bid",w′,"_to_","bid",w)), libbid), $BID, ($BID′,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())), nothing)
        end

        # promote binary*decimal -> decimal, for consistency with other operations above
        # (there doesn't seem to be any clear standard for this)
        if w′ <= 64
            FP′ = Symbol(string("Float",w′))
            @eval promote_rule(::Type{$BID}, ::Type{$FP′}) = $(Symbol(string("Dec",max(w,w′))))
            for (i′, i′str) in (("Int$w′", "int$w′"), ("UInt$w′", "uint$w′"))
                Ti′ = eval(Symbol(i′))
                if w > w′
                    @eval Base.convert(::Type{$BID}, x::$Ti′) = nox(ccall(($(bidsym(w,"from_",i′str)), libbid), $BID, ($Ti′,Ref{Cuint}), x, RefArray(flags, Threads.threadid())))
                else
                    @eval Base.convert(::Type{$BID}, x::$Ti′) = nox(ccall(($(bidsym(w,"from_",i′str)), libbid), $BID, ($Ti′,Cuint,Ref{Cuint}), x, roundingmode[Threads.threadid()], RefArray(flags, Threads.threadid())))
                end
            end
        end
    end

    for w′ in (8,16,32,64)
        for (i′, i′str) in (("Int$w′", "int$w′"), ("UInt$w′", "uint$w′"))
            Ti′ = eval(Symbol(i′))
            @eval begin
                Base.trunc(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xint")), libbid), $Ti′, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), InexactError, :trunc, $BID, x, mask=INVALID | OVERFLOW)
                Base.floor(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xfloor")), libbid), $Ti′, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), InexactError, :floor, $BID, x, mask=INVALID | OVERFLOW)
                Base.ceil(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xceil")), libbid), $Ti′, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), InexactError, :ceil, $BID, x, mask=INVALID | OVERFLOW)
                Base.round(::Type{$Ti′}, x::$BID, ::RoundingMode{:NearestTiesAway}) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xrninta")), libbid), $Ti′, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), InexactError, :round, $BID, x, mask=INVALID | OVERFLOW)
                Base.convert(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xfloor")), libbid), $Ti′, ($BID,Ref{Cuint}), x, RefArray(flags, Threads.threadid())), InexactError, :convert, $BID, x)
                Base.$(Symbol("$Ti′"))(x::$BID) = convert($Ti′, x)
            end
        end
    end

    @eval Base.bswap(x::$BID) = reinterpret($BID, bswap(x.x))
    @eval Base.convert(::Type{Float16}, x::$BID) = convert(Float16, convert(Float32, x))
    @eval Base.Float16(x::$BID) = convert(Float16, x)
    @eval Base.reinterpret(::Type{$Ti}, x::$BID) = x.x
end # widths w

Base.round(x::DecimalFloatingPoint, ::RoundingMode{:FromZero}) = signbit(x) ? floor(x) : ceil(x)

for (f) in (:trunc, :floor, :ceil)
    @eval Base.$f(::Type{Signed}, x::DecimalFloatingPoint) = $f(Int, x)
    @eval Base.$f(::Type{Unsigned}, x::DecimalFloatingPoint) = $f(UInt, x)
    @eval Base.$f(::Type{Integer}, x::DecimalFloatingPoint) = $f(Int, x)

    @eval function Base.$f(::Type{I}, x::DecimalFloatingPoint) where {I<:Integer}
        x′ = $f(x)
        typemin(I) <= x′ <= typemax(I) || throw(InexactError(Symbol($f), I, x))
        _, s, e = sigexp(x′)
        return I(flipsign(s * I(10)^e, x))
    end
end

Base.convert(::Type{Signed}, x::DecimalFloatingPoint) = convert(Int, x)
Base.convert(::Type{Unsigned}, x::DecimalFloatingPoint) = convert(UInt, x)
Base.convert(::Type{Integer}, x::DecimalFloatingPoint) = convert(Int, x)

function Base.convert(::Type{I}, x::DecimalFloatingPoint) where {I<:Integer}
    x != trunc(x) && throw(InexactError(:convert, I, x))
    typemin(I) <= x <= typemax(I) || throw(InexactError(:convert, I, x))
    _, s, e = sigexp(x)
    return I(flipsign(s * I(10)^e, x))
end

Base.Int128(x::DecimalFloatingPoint) = convert(Int128, x)
Base.UInt128(x::DecimalFloatingPoint) = convert(UInt128, x)

Base.round(::Type{T}, x::DecimalFloatingPoint) where {T<:Integer} = convert(T, round(x))
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:Nearest}) where {T<:Integer} = convert(T, round(x, RoundNearest))
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:NearestTiesAway}) where {T<:Integer} = convert(T, round(x, RoundNearestTiesAway))
function Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:NearestTiesUp}) where {T<:Integer}
    y = floor(T, x)
    ifelse(x==y, y, copysign(floor(T, 2*x-y), x))
end
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:ToZero}) where {T<:Integer} = trunc(T, x)
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:FromZero}) where {T<:Integer} = (x>=0 ? ceil(T, x) : floor(T, x))
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:Up}) where {T<:Integer} = ceil(T, x)
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:Down}) where {T<:Integer} = floor(T, x)

# the complex-sqrt function in base doesn't work for use, because it requires base-2 ldexp
function Base.sqrt(z::Complex{T}) where {T<:DecimalFloatingPoint}
    x, y = reim(z)
    x==y==0 && return Complex(zero(x),y)
    ρ = sqrt((abs(x) + hypot(x,y)) * 0.5)
    ξ = ρ
    η = y
    if isfinite(η) η=(η/ρ)/2 end
    if x<0
        ξ = abs(η)
        η = copysign(ρ,y)
    end
    Complex(ξ,η)
end

# see issue #92 — the fallback power_by_squaring fails for p < 0, and this is more accurate:
Base.:^(x::DecimalFloatingPoint, p::Integer) = x^oftype(x, p)
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{0}) = one(x)
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{1}) = x
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{2}) = x*x
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{3}) = x*x*x

# decompose is needed for hash() and comparison with Rational
function Base.decompose(x::DecimalFloatingPoint)::Tuple{BigInt, Int, BigInt}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(signbit(x), -1, 1), 0, 0
    iszero(x) && return 0, 0, ifelse(signbit(x), -1, 1)
    sign, s, e = sigexp(x)
    if e >= 0
        if e <= 27
            return s * BigInt(Int64(5)^e), e, sign
        else
            return s * BigInt(5)^e, e, sign
        end
    else
        e2 = -e
        q, r = divrem(s, oftype(s, 5))
        while e2 > 0 && r == 0
            s = q
            q, r = divrem(s, oftype(s, 5))
            e2 -= 1
        end
        if e2 <= 27
            return copysign(s, x), e, Int64(5)^e2
        else
            return copysign(s, x), e, BigInt(5)^e2
        end
    end
end

function Base.:(==)(dec::DecimalFloatingPoint, rat::Rational)
    if isfinite(dec)
        rat.den == 1 && return dec == rat.num
        t = rat.den
        t >>= trailing_zeros(t)
        q, r = divrem(t, 5)
        while r == 0
            t = q
            q, r = divrem(t, 5)
        end
        return t == 1 && dec*rat.den == rat.num
    else
        return dec == rat.num/rat.den
    end
end

Base.:(==)(dec::T, flt::Union{Float16,Float32,Float64}) where {T<:DecimalFloatingPoint} = dec == T(flt, RoundUp) == T(flt, RoundDown)
Base.:>(dec::T, flt::Union{Float16,Float32,Float64}) where {T<:DecimalFloatingPoint} = dec > T(flt, RoundDown)
Base.:<(dec::T, flt::Union{Float16,Float32,Float64}) where {T<:DecimalFloatingPoint} = dec < T(flt, RoundUp)
Base.:(>=)(dec::T, flt::Union{Float16,Float32,Float64}) where {T<:DecimalFloatingPoint} = dec >= T(flt, RoundUp)
Base.:(<=)(dec::T, flt::Union{Float16,Float32,Float64}) where {T<:DecimalFloatingPoint} = dec <= T(flt, RoundDown)
Base.:(==)(flt::Union{Float16,Float32,Float64}, dec::T) where {T<:DecimalFloatingPoint} = dec == flt
Base.:>(flt::Union{Float16,Float32,Float64}, dec::T) where {T<:DecimalFloatingPoint} = dec < flt
Base.:<(flt::Union{Float16,Float32,Float64}, dec::T) where {T<:DecimalFloatingPoint} = dec > flt
Base.:(>=)(flt::Union{Float16,Float32,Float64}, dec::T) where {T<:DecimalFloatingPoint} = dec <= flt
Base.:(<=)(flt::Union{Float16,Float32,Float64}, dec::T) where {T<:DecimalFloatingPoint} = dec >= flt

# used for next/prevfloat:
const pinf128 = _parse(Dec128, "+Inf")
const minf128 = _parse(Dec128, "-Inf")

for T in (Dec32, Dec64, Dec128)
    @eval begin
        Base.eps(::Type{$T}) = $(_sub(_nextfloat(one(T)), one(T)))
        Base.typemax(::Type{$T}) = $(_parse(T, "+inf"))
        Base.typemin(::Type{$T}) = $(_parse(T, "-inf"))
        Base.floatmax(::Type{$T}) = $(_prevfloat(_parse(T, "+inf")))
    end
end

Base.floatmin(::Type{Dec32}) = reinterpret(Dec32, 0x03000001) # Dec32("1.0e-95")
Base.floatmin(::Type{Dec64}) = reinterpret(Dec64, 0x01e0000000000001) # Dec64("1.0e-383")
Base.floatmin(::Type{Dec128}) = reinterpret(Dec128, 0x00420000000000000000000000000001) # Dec128("1.0e-6143")

Base.maxintfloat(::Type{Dec32}) = reinterpret(Dec32, 0x36000001) # Dec32("1e7")
Base.maxintfloat(::Type{Dec64}) = reinterpret(Dec64, 0x33c0000000000001) # Dec64("1e16")
Base.maxintfloat(::Type{Dec128}) = reinterpret(Dec128, 0x30840000000000000000000000000001) # Dec128("1e34")

Base.convert(::Type{F}, x::Union{Int8,UInt8,Int16,UInt16}) where {F<:DecimalFloatingPoint} = F(Int32(x))
Base.convert(::Type{F}, x::Integer) where {F<:DecimalFloatingPoint} = F(string(x))
Base.convert(::Type{F}, x::Rational) where {F<:DecimalFloatingPoint} = F(x.num) / F(x.den)
Base.convert(::Type{F}, x::Float16) where {F<:DecimalFloatingPoint} = F(Float32(x))
promote_rule(::Type{F}, ::Type{Float16}) where {F<:DecimalFloatingPoint} = F
promote_rule(::Type{F}, ::Type{T}) where {F<:DecimalFloatingPoint,T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}} = F

# fallback floating-point conversion via strings
Base.convert(::Type{F}, x::AbstractFloat) where {F<:DecimalFloatingPoint} = F(string(x))
Base.convert(::Type{F}, x::F) where {F<:DecimalFloatingPoint} = x # don't call AbstractFloat method!

# so that mathconsts get promoted to Dec32, not Dec64, like Float32
promote_rule(::Type{Irrational{s}}, ::Type{F}) where {s,F<:DecimalFloatingPoint} = F
promote_rule(::Type{Irrational{s}}, T::Type{Complex{F}}) where {s,F<:DecimalFloatingPoint} = T

Base.widen(::Type{Dec32}) = Dec64
Base.widen(::Type{Dec64}) = Dec128

macro d_str(s, flags...) parse(Dec64, s) end
macro d32_str(s, flags...) parse(Dec32, s) end
macro d64_str(s, flags...) parse(Dec64, s) end
macro d128_str(s, flags...) parse(Dec128, s) end

# for zero-padding in printing routines above
function printzeros(io::IO, n::Int)
    for i = 1:n
        print(io, '0')
    end
end

end # module


text:
abstract type DelayDiffEqAlgorithm <: AbstractDDEAlgorithm end
abstract type AbstractMethodOfStepsAlgorithm{constrained} <: DelayDiffEqAlgorithm end

struct MethodOfSteps{algType,AType,RType,NType,constrained} <:
    AbstractMethodOfStepsAlgorithm{constrained}

    alg::algType
    fixedpoint_abstol::AType
    fixedpoint_reltol::RType
    fixedpoint_norm::NType
    max_fixedpoint_iters::Int
end

"""
    MethodOfSteps(alg; constrained::Bool=false, fixedpoint_abstol=nothing,
                  fixedpoint_reltol=nothing, fixedpoint_norm=nothing,
                  max_fixedpoint_iters::Int=10)

Construct an algorithm that solves delay differential equations by the method of steps,
where `alg` is an ODE algorithm from OrdinaryDiffEq.jl without lazy interpolation upon
which the calculation of steps is based.

If the algorithm is `constrained` only steps of size at most the minimal delay will be
taken. If it is unconstrained, fixed-point iteration is applied for step sizes that exceed
the minimal delay.

The absolute and relative tolerance of the fixed-point iterations can be set by
`fixedpoint_abstol` and `fixedpoint_reltol`, respectively, either as scalars or vectors.
Based on these tolerances error estimates are calculated during the fixed-point iterations
with a norm that may be specified as `fixedpoint_norm`. Fixed-point iterations are stopped
if the error estimate is less than 1 or after the maximal number `max_fixedpoint_iters` of
iteration steps.
"""
Base.@pure function MethodOfSteps(alg; constrained::Bool=false, fixedpoint_abstol=nothing,
                                  fixedpoint_reltol=nothing, fixedpoint_norm=nothing,
                                  max_fixedpoint_iters::Int=10)
    MethodOfSteps{typeof(alg),typeof(fixedpoint_abstol),typeof(fixedpoint_reltol),
                  typeof(fixedpoint_norm),constrained}(alg,fixedpoint_abstol,
                                                       fixedpoint_reltol,
                                                       fixedpoint_norm, max_fixedpoint_iters)
end

ext: jl


text:
##    Error messages for Unicode / UTF support

struct UnicodeError <: Exception
    errmsg::AbstractString   ##< A UTF_ERR_ message
    errpos::Int32            ##< Position of invalid character
    errchr::UInt32           ##< Invalid character
end

show(io::IO, exc::UnicodeError) = print(io, replace(replace(string("UnicodeError: ",exc.errmsg),
    "<<1>>" => string(exc.errpos)),
    "<<2>>" => string(exc.errchr, base=16)))


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

using Markdown
using InteractiveUtils

# ╔═╡ ad96c6b4-add2-11eb-1a72-8993033e1259
using Plots, LaTeXStrings, QuadGK, Roots, DifferentialEquations

# ╔═╡ 5168971b-2b60-494c-a8d8-cdfa8e1b8102
md"# 1D Bistable Example"

# ╔═╡ 261580d9-f1a5-442c-8e7f-3479f993f9b2
md"_Define separate production and degradation terms_"

# ╔═╡ 02d0fcf2-ca87-49ae-a8f3-f14310d0182a
begin
	n_max = 200.
	n_min = 1 # need to avoid 0 for log we'll evaluate later
	n = n_min:0.1:n_max
	γ = 1
	ks = 1e-4
	ν0 = 12.5
	ν1 = 200
	production(n) = (ν0+ ν1 * ks*n^2)/(1+ks*n^2)
	degradation(n) = γ*n
end

# ╔═╡ 13ddc727-e9bd-44f2-ba5e-c6c311b41723
begin

	plot(production, n, label=L"f(n)")
	p1 = plot!(degradation, n, label=L"g(n)", xlabel=L"n", title="Rate of Change", scale=:log10, xlimits=(15, n_max), ylimits=(15, 200), legend=:topleft)
# plot!(n->production(n)  - degradation(n), n, label=L"\frac{dn}{dt}",
# 		ylims=(-20, 150), xlims=(0, 150), legend=:topleft, title="Time evolution")
end

# ╔═╡ 868ba091-d0d2-42a2-bca7-f9e996b3870d
md"_Compute the fixed points to split up the 2 attractors into **Basins of Attractions**_"

# ╔═╡ 348718a0-8e19-4911-8998-5de9689bce1d
fixed_points = find_zeros(n-> production(n) - degradation(n), n_min, n_max)

# ╔═╡ 51d6af76-5ed3-41f8-9351-91e683c7bb57
md"_Integrate to get the quasi-potential_"

# ╔═╡ acb0adcf-fc19-45cf-bf0f-eb5532743b1a
begin
	f = production
	g = degradation

	ϕ(n) = quadgk(n -> -2*(f(n) -g(n))/(f(n) + g(n)), n_min, n)[1]
end

# ╔═╡ b4c511ea-fc96-475e-a5fc-1c9e282ca0e5
ϕ_vec = ϕ.(n)

# ╔═╡ 3373f26d-a6a4-4378-a8e1-c011c2e9f7a9
p2 = plot(n, ϕ_vec, ylabel=L"\phi(n)", xlabel=L"n", label="", title="Quasi-potential")

# ╔═╡ e147802a-2cc5-4610-821b-bf351ed62f9f
md"_Show the the probability due to the quasi-potential dominates compared to the other term in the equation_"

# ╔═╡ 178564ff-da53-4dd5-ba1f-b3e5db24b2ea
begin
	plot(n->log10.(f(n) + g(n)), n, label=L"f + g")
	plot!(n, log10.(ℯ.^(-ϕ_vec)), label=L"e^{-\phi}")
	p3 = plot!(xaxis=:log10, legend=:topleft)
end

# ╔═╡ 33cf16fa-679d-4447-aa59-b6d48e8fd15a
md"_Compute the normalized probability_"

# ╔═╡ a7e9be12-a1d4-4f1a-bbfc-d2642e5a9dd5
begin
	fn = f.(n)
	gn = g.(n)
	p_unnormal(n) =  ℯ^(-ϕ(n)) / (f(n)+g(n))
	norm = quadgk(p_unnormal, n_min, n_max)[1]
	p_normal(n) = 1/norm * p_unnormal(n)
	p4 = plot(p_normal,n,  title = "Normalized Probability", label="")
end

# ╔═╡ f5b8a081-222f-4dd6-ae4c-af7060aa4223
md"_Set up the jump process simulation_"

# ╔═╡ 467386e3-2ee8-46d4-9b40-831ae2778a57
random_count() = floor(Int, n_min + rand() * n_max)

# ╔═╡ 6ae4ed84-80c2-4c06-b69a-c81e6485b78b
begin
	affect1!(integrator) = integrator.u += 1
	affect2!(integrator) = integrator.u -= 1
	prodJump = ConstantRateJump((u, p, t) -> f(u), affect1!)
	degJump = ConstantRateJump((u, p, t) -> g(u), affect2!)
	u0 = random_count()
	t = (0.0, 5000.0)
	dprob = DiscreteProblem(u0, t)
	jprob = JumpProblem(dprob, Direct(), prodJump, degJump)
end

# ╔═╡ f32b2ab8-a941-463b-bb27-968c09e24777
begin
	sim_res = solve(jprob, SSAStepper())
	p5 = plot(sim_res, xlabel=L"t", label=L"n(t)", title="Simulation trace")
end

# ╔═╡ 50c911b1-57cf-4de4-9d65-054318c445fe
begin
	pf = plot(p1, p2, p4, p5, layout = grid(2,2), size=(750, 750))
	# savefig(pf, "ch3_bistable_prob.pdf")
	pf
end

# ╔═╡ 4547d20c-dab1-4df2-aa54-b8c52b1b5eeb
md"_Setup functions to do a longer time simulation and compare to integrated probability distribution_"

# ╔═╡ 456090be-e1c0-452e-b0ad-98fbc9f0576f
function calculate_time_fraction(sim_result, split)
	dts = diff(sim_result.t)
	dt_vals = collect(zip(dts, map(last, sim_result.u[1:end])))
	high = sum(map(first, filter(dtu -> last(dtu) > split, dt_vals)))
	low = sum(map(first, filter(dtu -> last(dtu) < split, dt_vals)))
	high, low
end

# ╔═╡ 23494ad1-82e3-4eee-86e3-695511de502e
function big_sim(num_sim = 25)
	function prob_func(prob, i, repeat)
        remake(prob, u0=random_count())
    end
	num_steps = 0
    ensemble_prob = EnsembleProblem(jprob, prob_func=prob_func)
	sim = solve(ensemble_prob, SSAStepper(), trajectories=num_sim)
	totals = [0., 0.]
	for j in 1:num_sim
		totals.+= calculate_time_fraction(sim[j], fixed_points[2])
		num_steps += length(sim[j])
	end
	totals, num_steps
end

# ╔═╡ 0e38a5c3-d970-4844-b560-fe4ed5452d85
totals, num_steps = big_sim()

# ╔═╡ 573a4914-b37e-4c2b-bd23-c10b58a6d4db
begin
	sim_frac = totals./sum(totals)
end

# ╔═╡ c62dd4a3-80ce-4b1f-b5a1-6142d30d0880
md"Time fractions from simulations is $(sim_frac)"

# ╔═╡ 86ee8326-51e1-453d-8d83-fd7021613c76
begin
	p_low = quadgk(p_normal, n_min, fixed_points[2])[1]
	p_high = quadgk(p_normal, fixed_points[2], n_max)[1]
	comp_frac = p_high, p_low
end

# ╔═╡ 5cabf565-b573-4446-b34d-8812986b25fe
md"Time fractions from integration of probability is $(comp_frac)"

# ╔═╡ Cell order:
# ╠═ad96c6b4-add2-11eb-1a72-8993033e1259
# ╟─5168971b-2b60-494c-a8d8-cdfa8e1b8102
# ╟─261580d9-f1a5-442c-8e7f-3479f993f9b2
# ╠═02d0fcf2-ca87-49ae-a8f3-f14310d0182a
# ╠═13ddc727-e9bd-44f2-ba5e-c6c311b41723
# ╟─868ba091-d0d2-42a2-bca7-f9e996b3870d
# ╠═348718a0-8e19-4911-8998-5de9689bce1d
# ╟─51d6af76-5ed3-41f8-9351-91e683c7bb57
# ╠═acb0adcf-fc19-45cf-bf0f-eb5532743b1a
# ╠═b4c511ea-fc96-475e-a5fc-1c9e282ca0e5
# ╠═3373f26d-a6a4-4378-a8e1-c011c2e9f7a9
# ╟─e147802a-2cc5-4610-821b-bf351ed62f9f
# ╠═178564ff-da53-4dd5-ba1f-b3e5db24b2ea
# ╟─33cf16fa-679d-4447-aa59-b6d48e8fd15a
# ╠═a7e9be12-a1d4-4f1a-bbfc-d2642e5a9dd5
# ╟─f5b8a081-222f-4dd6-ae4c-af7060aa4223
# ╠═467386e3-2ee8-46d4-9b40-831ae2778a57
# ╠═6ae4ed84-80c2-4c06-b69a-c81e6485b78b
# ╠═f32b2ab8-a941-463b-bb27-968c09e24777
# ╠═50c911b1-57cf-4de4-9d65-054318c445fe
# ╟─4547d20c-dab1-4df2-aa54-b8c52b1b5eeb
# ╠═456090be-e1c0-452e-b0ad-98fbc9f0576f
# ╠═23494ad1-82e3-4eee-86e3-695511de502e
# ╠═0e38a5c3-d970-4844-b560-fe4ed5452d85
# ╟─c62dd4a3-80ce-4b1f-b5a1-6142d30d0880
# ╟─5cabf565-b573-4446-b34d-8812986b25fe
# ╠═573a4914-b37e-4c2b-bd23-c10b58a6d4db
# ╠═86ee8326-51e1-453d-8d83-fd7021613c76


text:
# 2D LSRTM/FWI parallelization test with 2 sources and 2 vintages
#
# Ziyi Yin, ziyi.yin@gatech.edu
# August 2021

@testset "test reduction from distributed tasks" begin
    @timeit TIMEROUTPUT "Test reduction" begin
        function judiWeightsConst(i::Int)
            return judiWeights(i*ones(Float32,10*i,10*i))
        end
        f = Vector{Future}(undef, 6)
        for i = 1:6
            f[i] = @spawn judiWeightsConst(i)
        end
        using JUDI:reduce!
        w = reduce!(f)
        @info "test if stacking follows the correct order through reduction"
        for i = 1:6
            @test w.data[i] == i*ones(Float32,10*i,10*i)
        end

        # Test objective function
        f1 = @spawn (Ref{Float32}(1f0), randn(10))
        f2 = @spawn (Ref{Float32}(2f0), randn(10))
        JUDI.local_reduce!(f1, f2)
        res = fetch(f1)

        @test res[1][] == 3f0
        @test typeof(res[1]) == Base.RefValue{Float32}

        res = JUDI.as_vec(res, Val(false))
        @test res[1] == 3f0
        @test typeof(res[1]) == Float32
    end
end

### Model
model, model0, dm = setup_model(tti, viscoacoustic, 4)
q, srcGeometry, recGeometry, f0 = setup_geom(model; nsrc=4)
q1 = q[[1,4]]
q2 = q[[2,3]]
srcGeometry1 = subsample(srcGeometry,[1,4])
srcGeometry2 = subsample(srcGeometry,[2,3])
recGeometry1 = subsample(recGeometry,[1,4])
recGeometry2 = subsample(recGeometry,[2,3])

opt = Options(sum_padding=true, free_surface=fs, f0=f0, dt_comp=dt)
F1 = judiModeling(model, srcGeometry1, recGeometry1; options=opt)
F2 = judiModeling(model, srcGeometry2, recGeometry2; options=opt)

# Observed data
dobs1 = F1*q1
dobs2 = F2*q2

# Perturbations
dm1 = 2f0*circshift(dm, 10)
dm2 = 2f0*circshift(dm, 30)

@testset "Multi-experiment arg processing" begin
    @timeit TIMEROUTPUT "Process multi exp" begin
        for (nm, m) in zip([1, 2], [model0, [model0, model0]])
            for (nq, q) in zip([1, 2], [q1, [q1, q2]])
                for (nd, d) in zip([1, 2], [dobs1, [dobs1, dobs2]])
                    for dmloc in [dm1, dm1[:]]
                        for (ndm, dm) in zip([1, 2], [dmloc, [dmloc, dmloc]])
                            args = m, q, d, dm
                            @test JUDI.check_args(args...) == maximum((nm, nq, nd, ndm))
                            for (a, expected) in zip(args, (model0, q1, dobs1, dmloc))
                                @test JUDI.get_exp(a, 1) == expected
                            end
                        end
                    end
                end
            end
        end

        @test_throws ArgumentError JUDI.check_args([model0, model0], [dobs1, dobs2, dobs1])
        @test_throws ArgumentError JUDI.check_args([model0, model0], [dobs1, dobs2], [dm1, dm2, dm1])
        @test_throws ArgumentError JUDI.check_args(model0, [dobs1, dobs2], [dm1, dm2, dm1])
    end
end

@testset "FWI/LSRTM objective multi-level parallelization test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
    @timeit TIMEROUTPUT "Multi FWI/LSRTM" begin
        ftol = 1f-5

        Jm0, grad = lsrtm_objective([model0, model0], [q1, q2], [dobs1, dobs2], [dm1, dm2]; options=opt, nlind=true)
        Jm01, grad1 = lsrtm_objective(model0, q1, dobs1, dm1; options=opt, nlind=true)
        Jm02, grad2 = lsrtm_objective(model0, q2, dobs2, dm2; options=opt, nlind=true)
        @test isapprox(Jm01+Jm02, Jm0; rtol=ftol)
        @test isapprox(grad1, grad[1]; rtol=ftol)
        @test isapprox(grad2, grad[2]; rtol=ftol)

        _Jm0, _grad = fwi_objective([model0, model0], [q1, q2], [dobs1, dobs2]; options=opt)
        _Jm01, _grad1 = fwi_objective(model0, q1, dobs1; options=opt)
        _Jm02, _grad2 = fwi_objective(model0, q2, dobs2; options=opt)
        @test isapprox(_Jm01+_Jm02, _Jm0; rtol=ftol)
        @test isapprox(_grad1, _grad[1]; rtol=ftol)
        @test isapprox(_grad2, _grad[2]; rtol=ftol)
    end
end


text:
function makefigure1(dist::DataFrame;
  Fret::Symbol = DIST_RET,
  figurepath::String = FIGURE_PATH)

  local sdist::SubDataFrame
  local ssdist::SubDataFrame
  local aggstd::DataFrame

  generatepath(panelstring::String) = "$figurepath\\figure1$panelstring-$(REPLICATION_TYPE[]).pdf"
  generateagg(df::AbstractDataFrame) = aggregatebyday(x->std(x)*100, df, selectcols=[:day, Fret])
  #sdist::SubDataFrame = view(dist, completecases(dist[!, [Fret, :day]]), :)

  #display(allstdXday)

  #=panel22 = plot(allstdXday, x=:day, y=Fret, Guide.xlabel("Day"), Guide.ylabel("Std"),
    Geom.line, Geom.point, Guide.title("All Dividend Events"))
  draw(PDF(figurepath, 9inch, 7inch),panel22)
  allstdXday |> CSV.write("output\\allstdXday-$(REPLICATION_TYPE[]).csv")=#

  #=panel22 = plot(allstdXday[!,Fret], allstdXday[!,:day],
    xlabel="Event Day",
    marker=:dot,
    ylabel = "SD/day in \\%")
  savefig(panel22, figurepath)=#

  #=panel22 = allstdXday|> @vlplot(
    title="All Dividend Events",
    mark={:line,
      point={filled=false}},
    x={:day,
      title="Event Day",
      type=:quantitative},
    y={Fret,
      axis={title="SD in %",
        tickCount=5,
        grid=false},
      type=:quantitative,
      scale={zero=false},
      },
    )=#
  sdist = view(dist, (r->r.primary && (!ismissing(r[Fret]))).(eachrow(dist)), :)
  aggstd = generateagg(sdist)
  panel22 = makedividendeventplot(aggstd,
    Fret = Fret,
    title = "All Dividend Events",
    yaxisname = "SD in %")
  save(generatepath("panel22"), panel22)

  ssdist = view(sdist,sdist.yield .< 0.0075, :)
  aggstd = generateagg(ssdist)
  panel11 = makedividendeventplot(aggstd,
    Fret = Fret,
    title = "Low-Yield: <0.75% Dividend",
    yaxisname = "SD in %")
  save(generatepath("panel11"), panel11)

  ssdist = view(sdist,(f::Float64-> 0.0075≤f && f<.015).(sdist.yield), :)
  aggstd = generateagg(ssdist)
  panel12 = makedividendeventplot(aggstd,
    Fret = Fret,
    title = "Mid-Yield: 0.75 to 1.5% Dividend",
    yaxisname = "SD in %")
  save(generatepath("panel12"), panel12)

  ssdist = view(sdist,0.015 .≤ sdist.yield, :)
  aggstd = generateagg(ssdist)
  panel21 = makedividendeventplot(aggstd,
    Fret = Fret,
    title = "Mid-Yield: >1.5% Dividend",
    yaxisname = "SD in %")
  save(generatepath("panel21"), panel21)
end

function makedividendeventplot(df::AbstractDataFrame;title::String = error("title is required"),
  xaxisname::String = "Event Day",
  yaxisname::String = error("yaxisname is required"),
  Fret::Symbol = error("Fret is required"))

  panel = df|> @vlplot(
    title=title,
    mark={:line,
      point={filled=false}},
    x={:day,
      axis={title=xaxisname,
        grid=true,
        gridDash=[2,2]},
      type=:quantitative},
    y={Fret,
      axis={title=yaxisname,
        grid=true,
        gridDash=[2,2]},
      type=:quantitative,
      scale={zero=false},
      },
    )

  return panel
end

function aggregatebyday(F::Function, dist::AbstractDataFrame; Fgroup::Symbol = :day,
  selectcols::Vector{Symbol} = names(dist), suffix::String = "")

  sdist::SubDataFrame = view(dist, :, selectcols)
  aggdf::DataFrame = aggregate(groupby(sdist, Fgroup), F)
  names!(aggdf, (n::Symbol->Symbol(replace(string(n), "_function"=>suffix))).(names(aggdf)))

  return aggdf
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using GeneralizedGenerated
using JuliaVariables
using Test
using BenchmarkTools
using DataStructures

rmlines = NGG.rmlines

@testset "no kwargs" begin

@gg function f1(a)
    quote
        x -> a + x
    end |> rmlines
end

@test f1(1)(2) == 3

@gg function f2(a)
    quote
        x -> begin
            a += 2
            x + a
        end
    end |> rmlines
end

@test f2(1)(2) == 5

@gg function f3(a)
    quote
        k = 20
        x -> begin
            a += 2
            x + a + k
        end
    end
end

@test f3(1)(2) == 25

end

@testset "kwargs" begin

@gg function f4(a)
    quote
        function z(x, k=1)
            x + 20 + a + k
        end
    end
end

@test f4(10)(2) == 33
end

@testset "namedtuple" begin

@gg function f5(a)
    quote
        function z(x, k=1)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=1, q=10)
end

@testset "mk funcs" begin

f_ = mk_function(:((x, y) -> x + y))
@test f_(1, 2) == 3

f_ = mk_function(:(function (x, y) x + y end))
@test f_(1, 2) == 3

end

@testset "type encoding more datatypes" begin

@gg function f5(a)
    tp = (1, 2, 3)
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=(1, 2, 3), q=10)

@gg function f5(a)
    tp = (a1=1, a2=2, a3=3)
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=(a1=1, a2=2, a3=3), q=10)

@gg function f5(a)
    tp = "233"
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k="233", q=10)

@gg function f5(a)
    tp = list(1, 2, 3)
    quote
        function z(x; k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=list(1, 2, 3), q=10)
@test f5(10)(2; k=10) == (x=2, k=10, q=10)

end

@testset "runtime eval" begin

a = to_type(:(1 + 2))
@test :(1 + 2) == from_type(a)
@test string(from_type(a)) == string(:(1 + 2))

@test runtime_eval(1) == 1
@test mk_function(:(
    x -> x + 1
))(2) == 3

@test_throws Any mk_function(quote
    x -> x + 1
end)

@test runtime_eval(quote
    x -> x + 1
end)(1) == 2

end

@testset "self recursive" begin
    to_test = quote
        g(x, r=0) = x === 0 ? r : begin
            g = g # required for self recur
            g(x-1, r + x)
        end
        g(10)
    end |> runtime_eval

    g(x, r=0) = x === 0 ? r : g(x-1, r + x)
    expected = g(10)
    @test expected == to_test
end

@testset "self recursive" begin
    to_test = quote
        g(x, r=0) = x === 0 ? r : begin
            g = g # required for self recur
            g(x-1, r + x)
        end
        g(10)
    end |> runtime_eval

    g(x, r=0) = x === 0 ? r : g(x-1, r + x)
    expected = g(10)
    @test expected == to_test
end

@testset "support where clauses and return type annotations for @gg" begin
    @gg function foo(x::T) where T
        :(x, T)
    end
    @test foo(1) == (1, Int)
    @gg function bar(x::T) where T
        quote
            g = x + 20
            x = 10
            () -> begin
                x = g
                x
            end
        end
    end
    @test bar(2)() == 2 + 20

    @gg function foobar(x::T, y::A) where {T <: Number, A <: AbstractArray{T}}
        quote
            g = x + 20
            x = 10
            () -> begin
                x = g
                (A, x + y[1])
            end
        end
    end
    @test foobar(2, [3])() == (Vector{Int}, 2 + 20 + 3)
end

@testset "support default arguments" begin
    @gg function h(x, c)
        quote
            d = x + 10
            function g(x, y=c)
                x + y + d
            end
        end
    end
    @test h(1, 2)(3) == 16
end

module S
    run(y) = y + 1
end

struct K
    f1::Function
    f2::Function
end
@testset "specifying evaluation modules" begin
    @gg m function g(m::Module, y) :(run(y)) end
    @test g(S, 1) == 2

    @gg m function h(m, y)
        quote
        c = m.f1(y)
        () -> begin c = m.f2(c) end
        end
    end
    k = K(x -> x + 1, x -> x * 9)
    next = h(k, 1)
    @test next() == 18
    @test next() == 18 * 9
end

@testset "test free variables of let bindings" begin
    @gg function test_free_of_let()
        quote
            let x = 1
                f = () -> begin
                    x * 3
                end
                x = 2
                f
            end
        end
    end
    @test test_free_of_let()() == 6
end

@testset "show something" begin
    f1 = mk_function(:(x -> x + 1))
    f2 = mk_function(:((x :: Int = 2, ) -> x + 1))
    @test f1(1) == 2
    @test f2() == 3
    println(f1)
    println(f2)
end

@testset "omit func argname: #34" begin
   f1 = mk_function(:( (:: Int) -> 0 ))
   @test f1(1) == 0
   @test_throws MethodError f1("")
end
# # From Chris Rackauckas: https://github.com/JuliaLang/julia/pull/32737
# @inline @generated function _invokefrozen(f, ::Type{rt}, args...) where rt
#     tupargs = Expr(:tuple,(a==Nothing ? Int : a for a in args)...)
#     quote
#         _f = $(Expr(:cfunction, Base.CFunction, :f, rt, :((Core.svec)($((a==Nothing ? Int : a for a in args)...))), :(:ccall)))
#         return ccall(_f.ptr,rt,$tupargs,$((:(getindex(args,$i) === nothing ? 0 : getindex(args,$i)) for i in 1:length(args))...))
#     end
# end
# # @cscherrer's modification of `invokelatest` does better on kwargs
# export invokefrozen
# @inline function invokefrozen(f, rt, args...; kwargs...)
#     g(kwargs, args...) = f(args...; kwargs...)
#     kwargs = (;kwargs...)
#     _invokefrozen(g, rt, (;kwargs...), args...)
# end
# @inline function invokefrozen(f, rt, args...)
#     _invokefrozen(f, rt, args...)
# end

# @testset "GeneralizedGenerated.jl" begin
#     # Write your own tests here.

#     quote
#         a = 1
#         [a](b,) -> a + b
#     end |>  closure_conv_static |> println
#     @generated function f(x)
#         quote
#             [x](a) ->  x + a
#         end |> gg
#     end
#     @test f(1)(2) == 3

#     g_f = mk_function([:x, :y], [], :(x + y))
#     @test g_f(10, 2) == 12

#     eval(:(hpre(x,y) = x+y))

#     static_f(x,y) = x+y

#     g_f(1, 2)
#     hpre(1, 2)
#     static_f(1, 2)

#     @info :GG_approach
#     @btime $g_f(1, 2)

#     @info :eval
#     @btime $hpre(1, 2)

#     @btime $static_f(1, 2)

# end


text:
# NOTE: this assumes "libsnappy" is an available library dependency and in the load path.

# Native ccall versions, temporary for comparisons
function ccall_compress(input::Vector{UInt8})
    size_in = length(input)
    output = zeros(UInt8, ccall((:snappy_max_compressed_length, "libsnappy"), Csize_t, (Csize_t,), size_in))
    outputindex = Ref{Csize_t}(length(output))
    ccall((:snappy_compress, "libsnappy"),
        Cint,
        (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Ref{Csize_t}),
        input, size_in, output, outputindex)
    resize!(output, outputindex[])
    return output
end

function ccall_uncompress(input::Array{UInt8})
    size_in = length(input)
    uncompressed_len = Ref{Csize_t}(0)
    ccall((:snappy_uncompressed_length, "libsnappy"),
        Cint,
        (Ptr{UInt8}, Csize_t, Ref{Csize_t}),
        input, size_in, uncompressed_len)
    output = zeros(UInt8, uncompressed_len[])
    ccall((:snappy_uncompress, "libsnappy"),
        Cint,
        (Ptr{UInt8}, Csize_t, Ptr{UInt8}, Ref{Csize_t}),
        input, size_in, output, uncompressed_len)
    return output
end


text:
# Conversion of all complex variables to real ones in the Poly structures

"""
  realPart, imagPart = cplx2real(expo::Exponent)

  Convert a complex Exponent in complex variables into `realPart` and
  `imagPart` polynomials of twice as many variables, real and imag parts of
  `expo` variables. Done recursively with the `cplx2real_rec` function.
"""
function cplx2real(expo::Exponent)
  vars, inds = collect(keys(expo.expo)), collect(values(expo.expo))
  return cplx2real_rec(vars, inds, Polynomial()+1, Polynomial()+0, length(expo)+1, Degree(0,0))
end

"""
realPart, imagPart = cplx2real_rec(vars::Array{Variable}, degs::Array{Degree}, realPart::Polynomial, imagPart::Polynomial, cur_ind::Int, cur_deg::Degree)

Transform recursively the complex exponent represented by the `vars` and `degs`
arrays into its real and imag parts, functions of its imag and real part
variables.
`cur_ind` decreases to 0, `cur_deg` decreases to Degree(0,0) for each step of
cur_ind. Terminaison case is reached at 0, Degree(0,0).
Initial arrays `vars` and `degs` are read only.

### Arguments
- vars::Array{Variable}
- degs::Array{Degree}
- realPart::Polynomial
- imagPart::Polynomial
- cur_ind::Int
- cur_deg::Degree
"""
function cplx2real_rec(vars::Array{Variable}, degs::Array{Degree}, realPart::Polynomial, imagPart::Polynomial, cur_ind::Int, cur_deg::Degree)
  ## Final case:
  if cur_ind == 1 && cur_deg == Degree(0,0)
    return (realPart, imagPart)
  ## One less variable to deal with:
  elseif cur_deg == Degree(0,0)
    return cplx2real_rec(vars, degs, realPart, imagPart, cur_ind-1, degs[cur_ind-1])
  ## Recursion rule, decrease the current variable exponent until it reaches Degree(0,0):
  else
    var = vars[cur_ind]
    if iscomplex(var)
      var_real, var_imag = varname_cplx2real(var.name)
      var_R, var_I = Variable(var_real, Real), Variable(var_imag, Real)
      if cur_deg.explvar > 0
        return cplx2real_rec(vars, degs, var_R * realPart - var_I * imagPart, var_R * imagPart + var_I * realPart, cur_ind, Degree(cur_deg.explvar-1, cur_deg.conjvar))
      elseif cur_deg.conjvar > 0
        cur_deg.explvar == 0 || warn("cur_deg.explvar should be 0 (and not $(cur_deg.explvar)), set to this value")
        return cplx2real_rec(vars, degs, var_R * realPart + var_I * imagPart, var_R * imagPart - var_I * realPart, cur_ind, Degree(0, cur_deg.conjvar-1))
      end
    elseif isbool(var)
      return cplx2real_rec(vars, degs, var*realPart, var*imagPart, cur_ind, Degree(0,0))
    else
      return cplx2real_rec(vars, degs, Exponent(Dict(var=>cur_deg))*realPart, Exponent(Dict(var=>cur_deg))*imagPart, cur_ind, Degree(0,0))
    end
  end
end

"""
  realPart, imagPart = cplx2real(pol::Polynomial)

  Convert a complex polynomial in complex variables into `realPart` and
  `imagPart` polynomials of twice as many variables, real and imag parts of
  `pol` variables.
"""
function cplx2real(pol::Polynomial)
  realPart = Polynomial()
  imagPart = Polynomial()

  for (expo, λ) in pol
    realexpo, imagexpo = cplx2real(expo)

    realPart += realexpo*real(λ) - imagexpo*imag(λ)
    imagPart += imagexpo*real(λ) + realexpo*imag(λ)
  end
  return (realPart, imagPart)
end

function cplx2real(pt_C::Point)
  pt = Point()
  for (var, val) in pt_C
    if var.kind <: Complex
      var_real, var_imag = varname_cplx2real(var.name)
      if real(val) != 0
        pt[Variable(var_real, Real)] = real(val)
      end
      if imag(val) != 0
        pt[Variable(var_imag, Real)] = imag(val)
      end
    else
      if val != 0
        pt[var] = real(val)
      end
    end
  end
  return pt
end

function real2cplx(pt::Point)
    ptC = Point()
    for (var, val) in pt
        if ismatch(r"(_Re|_Im)", var.name)
            if ismatch(r"_Re", var.name)
                varname = var.name[1:(end-3)]
                ptC[Variable(varname, Complex)] = val + im*pt[Variable(varname*"_Im", Real)]
            end
        else
            ptC[var] = val
        end
    end
    return ptC
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# This is a Julia script for plotting the
# orbits that show oscillations (position vs. time)

# Define the floating point type used across the script
F = Float64

# Define the integer type used across the script
I = UInt64

# Define default values for optional arguments
POSTFIX = ""
H = 0.01

"Check if the value of the option is the last argument"
function check_last(i)
    if i + 1 == length(ARGS)
        println("The last argument should be the input directory.")
        exit(1)
    end
end

# Parse the options
for i in eachindex(ARGS)
    # Time step
    if ARGS[i] == "-h"
        check_last(i)
        try
            global H = parse(F, ARGS[i+1])
        catch
            println("Couldn't parse the value of the `-h` argument.")
            exit(1)
        end
    end
    # A postfix for the names of output files
    if ARGS[i] == "--postfix"
        check_last(i)
        try
            global POSTFIX = " ($(ARGS[i+1]))"
        catch
            println("Couldn't parse the value of the `--postfix` argument.")
            exit(1)
        end
    end
end

# Check for required arguments
if length(ARGS) == 0
    println("""
        Usage:
        julia --project=. scripts/oscillations.jl [-h <H>] [--postfix <POSTFIX>] <INPUT> """
    )
    exit(1)
end

# Multiply the time step
H *= π / 2

# Define the input directory
INPUT_DIR = ARGS[end]

println('\n', " "^4, "> Loading the packages...")

using LaTeXStrings
using Plots

# Use the GR backend for plots
gr()

# Change some of the default parameters for plots
default(fontfamily = "Computer Modern", dpi = 300, legend = :topright)

# Define the paths to output directories
CURRENT_DIR = @__DIR__
ROOT_DIR = basename(CURRENT_DIR) == "scripts" ? dirname(CURRENT_DIR) : CURRENT_DIR
PLOTS_DIR = joinpath(ROOT_DIR, "plots")
OUTPUT_DIR = joinpath(PLOTS_DIR, "oscillations")

# Make sure the needed directories exist
mkpath(OUTPUT_DIR)

# Define the paths to the binary files
z_path = joinpath(INPUT_DIR, "z.bin")

"Read binary files in the `bincode` format"
function read_bincode(path::AbstractString)::Tuple{I, Vector{F}}
    open(path, "r") do io
        n = read(io, I)
        a = reinterpret(F, read(io, sizeof(F) * n))
        n, a
    end
end

# Plot the phase space portrait if the corresponding data files exist
if isfile(z_path)
    # Read the data
    n, z = read_bincode(z_path)

    println(" "^4, "> Plotting the orbit...")

    # Plot the phase space portrait
    p = plot(
        [ i * H / (2 * π) for i in 0:(n-1) ],
        z;
        label = "",
        title = "Orbit",
        xlabel = L"t \; \mathrm{[2 \pi]}",
        ylabel = L"z",
        xminorticks = 5,
    )

    # Save the figure as PDF and PNG
    savefig(p, joinpath(OUTPUT_DIR, "Orbit$(POSTFIX).pdf"))
    savefig(p, joinpath(OUTPUT_DIR, "Orbit$(POSTFIX).png"))
end

println()


text:
__precompile__(true)

module SimplexSplitting

using LinearAlgebra
using Printf
using Distributions
using Parameters

include("../even_sampling.jl")
include("complementary.jl")
include("simplex_split.jl")
include("simplex_subdivision.jl")
include("simplex_subdivision_single.jl")
include("heaviside.jl")
include("embed.jl")
include("triangulate.jl")
include("refine_triangulation.jl")
include("refine_triangulation_with_images.jl")
include("invariantset.jl")
include("centroids_radii.jl")
include("simplexvolumes.jl")
include("refine_recursive.jl")
include("refine_recursive_withimages.jl")
include("refine_t.jl")
include("refine_variable_k.jl")
include("orientations.jl")

export tensordecomposition, simplex_split, simplicial_subdivision_single,
simplicial_subdivision, embed, Embedding, embedding, embedding_ex, triangulate, Triangulation, triang_from_embedding, embedding_example, refine_triangulation, simplex_volumes,
refine_triangulation_images, invariantset, centroids_radii2, refine_recursive, refine_recursive_images, refine_variable_k!, refine_variable_k_new, refine_variable_k_newnew, query_refinement,
refine_t!,
gaussian_embedding, centroids_radii2, example_triangulation,
gaussian_embedding_arr,
random_embedding,
orientations,
Simplex,
find_simplex, find_imsimplex,
get_simplices, get_imagesimplices,
maybeintersecting_simplices, maybeintersecting_imsimplices
end # module


text:
const REQUEST = 0
const RESPONSE = 1
const NOTIFICATION = 2

import Distributed.RemoteChannel

mutable struct NvimClient{S} <: NvimObject
    input::S # Input to nvim
    output::S

    channel_id::Int
    next_reqid::Int
    waiting::Dict{Int,RemoteChannel}
    reader::Task
    NvimClient(i, o, c, n, w, r) = new{typeof(i)}(i, o, c, n, w, r)
    NvimClient(i, o, c, n, w) = new{typeof(i)}(i, o, c, n, w)
end

function NvimClient(input, output, handler=DummyHandler())
    c = NvimClient(input, output, -1, 0, Dict{Int,RemoteChannel}())
    c.reader = @async readloop(c, handler)
    c.channel_id, metadata = send_request(c, "nvim_get_api_info", [])
    if symbolize(metadata) != _metadata
        println("warning: possibly incompatible api metadata")
    end
    c
end

# this is probably not most efficient in the common case (no contention)
# but it's the simplest way to assure task-safety of reading from the stream
function readloop(c::NvimClient, handler)
    while !eof(c.output)
        msg = unpack(c.output)
        kind = msg[1]
        if kind == RESPONSE
            serial = msg[2]
            ref = pop!(c.waiting, serial)
            put!(ref, (msg[3], msg[4]))
        elseif kind == NOTIFICATION
            name = String(msg[2])
            args = retconvert(Any, c, msg[3])
            try
                on_notify(handler, c, name, args)
            catch err
                logerr(err, catch_backtrace(), "handler", "notification", name, args)
            end
        elseif kind == REQUEST
            serial = msg[2]
            method = String(msg[3])
            args = retconvert(Any, c, msg[4])
            try
                on_request(handler, c, serial, method, args)
            catch err
                reply_error(c, serial, "Caught Exception in request handler: $err")
                logerr(err, catch_backtrace(), "handler", "request", method, args)
            end
        end
        flush(stderr)
    end
end

function logerr(e::Exception, bt, wher, what, name, args)
    println(stderr, "Caught Exception in $wher for $what \"$name\" with arguments $args:")
    showerror(stderr, e, bt)
    println(stderr, "")
    flush(stderr)
end

Base.wait(c::NvimClient) = wait(c.reader)

# we cannot use pack(stream, msg) as it's not synchronous
_send(c, msg) = (write(c.input, pack(msg)), flush(c.input))

function send_request(c::NvimClient, meth, args)
    reqid = c.next_reqid
    c.next_reqid += 1
    # TODO: are these things cheap to alloc or should they be reused
    rc = RemoteChannel()
    c.waiting[reqid] = rc
    meth = string(meth)

    _send(c, Any[REQUEST, reqid, meth, args])
    (err, res) = take!(rc) # blocking
    # TODO: make these recoverable
    if err !== nothing
        error(string(meth, ": ", String(err[2])))
    end
    res
end

function reply_error(c, serial, err)
    _send(c, Any[RESPONSE, serial, err, nothing])
end

function reply_result(c, serial, res)
    _send(c, Any[RESPONSE, serial, nothing, res])
end

# when overriding these, note that this runs in the reader task,
# use @async/@spawn when doing anything long-running or blocking,
# including method calls to nvim
struct DummyHandler; end
function on_notify(::DummyHandler, c, name, args)
    println("notification: $name $args")
end

function on_request(::DummyHandler, c, serial, name, args)
    println("WARNING: ignoring request $name $args, please override `on_request`")
    reply_error(c, serial, "Client cannot handle request, please override `on_request`")
end


text:
using LatexLH, Test

function symbol_info_test()
    @testset "SymbolInfo" begin
        si = SymbolInfo(:capShare, "\\alpha", "Capital share", "Technology");
        println(si);
        @test description(si) == "Capital share"
        @test group(si) == "Technology"
        ncStr = newcommand(si);
        println(ncStr)
    end
end

function symbol_table_test()
    @testset "SymbolTable" begin
        s = SymbolTable();
        @test isempty(s)
        @test length(s) == 0

        nGroups = 3;
        nSyms = 4;
        s = LatexLH.test_symbol_table(nGroups, nSyms);
        @test length(s) == nGroups * nSyms

        println(s);
        gl = LatexLH.group_list(s);
        @test length(gl) == nGroups

        fPath = joinpath(testDir, "preamble.tex");
        open(fPath, "w") do io
            write_preamble(io, s);
        end

        fPath = joinpath(testDir, "notation_tex.tex");
        open(fPath, "w") do io
            write_notation_tex(io, s);
        end

        erase!(s);
        @test isempty(s)

        s2 = SymbolTable(["capShare"  "\\alpha"  "Capital share"  "Technology";
            "tfp"  "A"  "Total factor productivity"  "Technology"]);
        @test length(s2) == 2
        @test has_symbol(s2, :capShare)
    end
end

function write_st_csv()
    fPath = joinpath(testDir, "st_csv.csv");
    open(fPath, "w") do io
        for ig = 1 : 3
            println(io, "G$ig,,,");
            for j = 1 : 2
                println(io, ",n$ig$j,l$ig$j,descr$ig$j");
            end
        end
    end
    return fPath
end

function load_st_test()
    fPath = write_st_csv();
    @test isfile(fPath)
    s = SymbolTable(fPath);
    @test has_symbol(s, :n11)
    @test has_symbol(s, :n32)
    si = s[:n21];
    @test group(si) == "G2"
    @test latex(si) == "l21"
end

@testset "SymbolTable" begin
    symbol_info_test();
    symbol_table_test();
    load_st_test();
end

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


text:
using YAXArrayBase
using Test
@testset "Arrays" begin
  include("arrays.jl")
end
@testset "Datasets" begin
include("datasets.jl")
end


text:
@testset "Ensembles" begin
  d = PointSet(rand(2,100))
  r = (z=[1:100 for i in 1:10],)
  s = Ensemble(d, r)
  @test domain(s) == d
  @test s[:z] == r[:z]
  for i in 1:10
    @test s[i] == georef((z=1:100,), d)
  end

  s1 = Ensemble(d, r)
  s2 = Ensemble(d, r)
  @test s1 == s2

  @test sprint(show, s) == "2D Ensemble"
  @test sprint(show, MIME"text/plain"(), s) == "2D Ensemble\n  domain: 100 PointSet{2,Float64}\n  variables: z\n  N° reals:  10"

  d = CartesianGrid(10,10)
  r = (z=[1:100 for i in 1:10],)
  s = Ensemble(d, r)
  @test domain(s) == d
  @test s[:z] == r[:z]
  for i in 1:10
    @test s[i] == georef((z=1:100,), d)
  end

  @test sprint(show, s) == "2D Ensemble"
  @test sprint(show, MIME"text/plain"(), s) == "2D Ensemble\n  domain: 10×10 CartesianGrid{2,Float64}\n  variables: z\n  N° reals:  10"

  if visualtests
    @test_reference "data/ensemble.png" plot(s,size=(800,300))
  end

  grid = CartesianGrid(3,3)
  reals = (z=[i*ones(nelements(grid)) for i in 1:3],)
  ensemble = Ensemble(grid, reals)

  # mean
  mean2D = mean(ensemble)
  @test mean2D.z == 2.0*ones(nelements(mean2D))
  @test domain(mean2D) == ensemble.domain

  # variance
  var2D = var(ensemble)
  @test var2D.z == 1.0*ones(nelements(var2D))
  @test domain(var2D) == ensemble.domain

  # quantile (scalar)
  p = 0.5
  quant2D = quantile(ensemble, p)
  @test quant2D.z == 2.0*ones(nelements(quant2D))
  @test domain(quant2D) == ensemble.domain

  # quantile (vector)
  ps = [0.0, 0.5, 1.0]
  quants2D = quantile(ensemble, ps)
  @test quants2D[2].z == quant2D.z
end


text:
# rbf kernel: exp(-gamma*|u-v|^2)

type Rbfk <: Kernel; n; s; w; v; x; y; dy; gamma;
    function Rbfk(nclass;o...)
        l = setparam!(new(nclass); o...)
        isdefined(l,:gamma) || (l.gamma = 1)
        return l
    end
end

# copy(l::Rbfk; o...)=Rbfk(s=copy(l.s), w=copy(l.w), v=copy(l.v), gamma=l.gamma)
setparam!(l::Rbfk; o...)=(for (n,v) in o; l.(n)=v; end; l)

function kernel(l::Rbfk)
    s2 = sum(l.s.^2, 1)
    x2 = sum(l.x.^2, 1)
    exp(-l.gamma * broadcast(+, x2, broadcast(+, s2', -2*(l.s' * l.x))))
end


text:
# This file includes stuff that both Hubbard models use

# conf === hsfield === discrete Hubbard Stratonovich field (Hirsch field)
const HubbardConf = Array{Int8, 2}
const HubbardDistribution = (Int8(-1), Int8(1))

abstract type HubbardModel <: Model end

"""
    HubbardModel([args..., ]; U[, kwargs...])

Creates an attractive or repulsive Hubbard model.
"""
function HubbardModel(args...; U, kwargs...)
    if U > 0.0
        HubbardModelRepulsive(args..., U = U; kwargs...)
    else
        HubbardModelAttractive(args..., U = -U; kwargs...)
    end
end

function choose_lattice(::Type{<: HubbardModel}, dims, L)
    if dims == 1
        return Chain(L)
    elseif dims == 2
        return SquareLattice(L)
    else
        return CubicLattice(dims, L)
    end
end

interaction_matrix_type(::Type{DQMC}, ::HubbardModel) = Diagonal{Float64, Vector{Float64}}

function init_interaction_matrix(m::HubbardModel)
    N = length(lattice(m))
    flv = nflavors(m)
    Diagonal(zeros(Float64, N*flv))
end

# deprecate this? maybe add nstates instead to avoid accessing model.flv?
@inline nflavors(m::HubbardModel) = m.flv
@inline lattice(m::HubbardModel) = m.l

# implement `DQMC` interface: mandatory
@inline function Base.rand(::Type{DQMC}, m::HubbardModel, nslices::Int)
    rand(HubbardDistribution, length(m.l), nslices)
end

# implement DQMC interface: optional
# Green's function is real for the repulsive Hubbard model.
@inline greenseltype(::Type{DQMC}, m::HubbardModel) = Float64

# See configurations.jl - compression of configurations
compress(mc::DQMC, ::HubbardModel, c) = BitArray(c .== 1)
compressed_conf_type(::Type{<: DQMC}, ::Type{<: HubbardModel}) = BitArray
function decompress(mc::DQMC{M, CB, CT}, ::HubbardModel, c) where {M, CB, CT}
    CT(2c .- 1)
end

include("HubbardModelAttractive.jl")
include("HubbardModelRepulsive.jl")


text:
mutable struct Object
l::Float64
μ::Float64
ϵ::Float64
m::Float64
R::Float64
M::Matrix{Float64}
M_inv::Matrix{Float64}

    function Object(l,μ,ϵ,m,R,M)
        new(l,μ,ϵ,m,R,M,inv(M))
    end
end

mutable struct Physics
    dt::Float64
    g::Vector{Float64}
    Δ::Float64
    T::Int

    function Physics(dt,g,Δ,T)
        new(dt,g,Δ,T)
    end
end

mutable struct Walker{T}
    M::T
    m::T
    I::T
    l::T
    c::T
    g::T
    γ::T
    use_controller::Bool
    tf::T
    Kp::T
    Θf::T
    Θ₂::T
    Θ₂̇ ::T
    t₀::T
end


text:
# Convert an AbstractNLPModel to a MathProgBase model that can be passed
# to a MathProgBase standard solver.
# Modeled after Miles Lubin's AmplMathProgInterface

using MathProgBase

export NLPModelEvaluator, NLPtoMPB

mutable struct NLPModelEvaluator{T <: AbstractNLPModel} <: SolverInterface.AbstractNLPEvaluator
  nlp :: T
end

MathProgBase.initialize(::NLPModelEvaluator, requested_features) = nothing

# SAutodetect?
MathProgBase.features_available(::NLPModelEvaluator) = [:Grad, :Jac, :JacVec, :HessVec, :Hess]
MathProgBase.features_available(::NLPModelEvaluator{NLPModels.ADNLPModel}) = [:Grad, :Jac, :Hess]

MathProgBase.eval_f(d::NLPModelEvaluator, x) = obj(d.nlp, x)

# use grad! ?
MathProgBase.eval_grad_f(d::NLPModelEvaluator, g, x) = copyto!(g, grad(d.nlp, x))

MathProgBase.eval_g(d::NLPModelEvaluator, g, x) = copyto!(g, cons(d.nlp, x))

function MathProgBase.jac_structure(d::NLPModelEvaluator)
  rows, cols, _ = jac_coord(d.nlp, [0.317i for i = 1:d.nlp.meta.nvar])
  return rows, cols
end

function MathProgBase.eval_jac_g(d::NLPModelEvaluator, J, x)
  _, _, vals = jac_coord(d.nlp, x)
  copyto!(J, vals)
end

# use jprod! ?
function MathProgBase.eval_jac_prod(d::NLPModelEvaluator, y, x, w)
  y = jprod(d.nlp, x, w)
end

# use jtprod! ?
function MathProgBase.eval_jac_prod_t(d::NLPModelEvaluator, y, x, w)
  y = jtprod(d.nlp, x, w)
end

function MathProgBase.hesslag_structure(d::NLPModelEvaluator)
  rows, cols, _ = hess_coord(d.nlp, [0.317i for i = 1:d.nlp.meta.nvar],
                             y=[0.618i for i = 1:d.nlp.meta.ncon])
  return rows, cols
end

function MathProgBase.eval_hesslag(d::NLPModelEvaluator, H, x, σ, μ)
  rows, cols, vals = hess_coord(d.nlp, x, y=μ, obj_weight=σ)
  copyto!(H, vals)
end

function MathProgBase.eval_hesslag_prod(d::NLPModelEvaluator, h, x, v, σ, μ)
  hprod!(d.nlp, x, v, h, y=μ, obj_weight=σ)
end

MathProgBase.isobjlinear(d::NLPModelEvaluator) = d.nlp.meta.nlo == 0

MathProgBase.isobjquadratic(d::NLPModelEvaluator) = false

MathProgBase.isconstrlinear(d::NLPModelEvaluator, i::Integer) = i in d.nlp.meta.lin

function loadNLPModel!(m::MathProgBase.AbstractMathProgModel, nlp::T) where T <: AbstractNLPModel
  MathProgBase.loadproblem!(m, nlp.meta.nvar, nlp.meta.ncon,
                            nlp.meta.lvar, nlp.meta.uvar,
                            nlp.meta.lcon, nlp.meta.ucon,
                            nlp.meta.minimize ? :Min : :Max,
                            NLPModelEvaluator(nlp))
  MathProgBase.setwarmstart!(m, nlp.meta.x0)

  nnlvar = max(nlp.meta.nlvc, nlp.meta.nlvo)
  narcvar = nlp.meta.nwv
  nlinvar = nlp.meta.nvar - (nnlvar + narcvar + nlp.meta.nbv + nlp.meta.niv)
  nbinvar = nlp.meta.nbv
  nintvar = nlp.meta.niv

  v = fill(:Cont, nlp.meta.nvar)
  # First populate Table 4
  varidx = 1
  for i = 1:(nlp.meta.nlvb - nlp.meta.nlvbi)
    varidx += 1
  end
  for i = 1:nlp.meta.nlvbi
    v[varidx] = :Int
    varidx += 1
  end
  for i = 1:(nlp.meta.nlvc - (nlp.meta.nlvb + nlp.meta.nlvci))
    varidx += 1
  end
  for i = 1:nlp.meta.nlvci
    v[varidx] = :Int
    varidx += 1
  end
  for i = 1:(nlp.meta.nlvo - (nlp.meta.nlvc + nlp.meta.nlvoi))
    varidx += 1
  end
  for i = 1:nlp.meta.nlvoi
    v[varidx] = :Int
    varidx += 1
  end
  # Now populate Table 3
  varidx += narcvar + nlinvar
  for i = 1:nbinvar
    v[varidx] = :Bin
    varidx += 1
  end
  for i = 1:nintvar
    v[varidx] = :Int
    varidx += 1
  end
  @assert varidx == nlp.meta.nvar + 1
  # Set variable types
  if any(vtype -> vtype == :Int || vtype == :Bin, v)
      MathProgBase.setvartype!(m, v)
  end
  return m
end

"""
    mp = NLPtoMPB(nlp, solver)

Return a `MathProgBase` model corresponding to an `AbstractNLPModel`.

#### Arguments
- `nlp::AbstractNLPModel`
- `solver::AbstractMathProgSolver` a solver instance, e.g., `IpoptSolver()`

Currently, all models are treated as nonlinear models.

#### Return values
The function returns a `MathProgBase` model `mpbmodel` such that it should
be possible to call

    MathProgBase.optimize!(mpbmodel)
"""
function NLPtoMPB(nlp :: T, solver :: MathProgBase.AbstractMathProgSolver) where T <: AbstractNLPModel
  model = MathProgBase.NonlinearModel(solver)
  return loadNLPModel!(model, nlp)
end

# MathProgNLPModels are straightforward
# Note that this method doesn't appear to be shadowed by the parametric NLPtoMPB
function NLPtoMPB(model :: Union{MathProgNLPModel, MathProgNLSModel},
                  solver :: MathProgBase.AbstractMathProgSolver)
  mpbmodel = MathProgBase.NonlinearModel(solver)
  MathProgBase.loadproblem!(mpbmodel,
                            model.meta.nvar, model.meta.ncon,
                            model.meta.lvar, model.meta.uvar,
                            model.meta.lcon, model.meta.ucon,
                            model.meta.minimize ? :Min : :Max,
                            NLPModelEvaluator(model))
  MathProgBase.setwarmstart!(mpbmodel, model.meta.x0)
  return mpbmodel
end


text:
# # Tutorial
#
# Here we will see what are the basic building blocks of SplitApplyPlot, and how to
# combine them to create complex plots based on tables or other formats.
#
# ## Basic building blocks
#
# The most important functions are `mapping`, and `visual`.
# `mapping` determines the mappings from data to plot. Its positional arguments correspond to
# the `x`, `y` or `z` axes of the plot, whereas the keyword arguments correspond to plot
# attributes that can vary continuously or discretely, such as `color` or `markersize`.
# Variables in `mapping` are split according to the categorical attributes in it, and then converted
# to plot attributes using a default palette.
# `visual` can be used to give data-independent visual information about the plot
# (plotting function or attributes).
#
# Finally, `data` determines what is the dataset to be used.
#
# ## Operations
#
# The outputs of `mapping`, `visual`, and `data` can be combined with `+` or `*`,
# to generate the specification of a complex plot.
#
# The operation `+` is used to create separate layers. `a + b` has as many layers as `la + lb`,
# where `la` and `lb` are the number of layers in `a` and `b` respectively.
#
# The operation `a * b` create `la * lb` layers, where `la` and `lb` are the number of layers
# in `a` and `b` respectively. Each layer of `a * b` contains the combined information of
# the corresponding layer in `a` and the corresponding layer in `b`. In simple cases,
# however, both `a` and `b` will only have one layer, and `a * b` simply combines the
# information.
#
# ## Working with tables

using RDatasets: dataset
using SplitApplyPlot, CairoMakie
mpg = dataset("ggplot2", "mpg");
cols = mapping(:Displ => "Displacement", :Hwy => "Highway miles");
grp = mapping(color = :Cyl => nonnumeric => "Cylinders");
scat = visual(Scatter)
pipeline = cols * scat
draw(data(mpg) * pipeline)
AbstractPlotting.save("scatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](scatter.svg)
#
# Now let's simply add `grp` to the pipeline to color according to `:Cyl`.

draw(data(mpg) * grp * pipeline)
AbstractPlotting.save("grouped_scatter.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](grouped_scatter.svg)
# Traces can be added together with `+`.

pipenew = cols * (scat + linear())
draw(data(mpg) * pipenew)
AbstractPlotting.save("linear.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](linear.svg)
# We can put grouping in the pipeline (we get a warning because of a degenerate group).

draw(data(mpg) * grp * pipenew)
AbstractPlotting.save("grouped_linear.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](grouped_linear.svg)
# This is a more complex example, where we split the scatter plot,
# but do the linear regression with all the data.

different_grouping = grp * scat + linear()
draw(data(mpg) * cols * different_grouping)
AbstractPlotting.save("semi_grouped.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](semi_grouped.svg)
#


text:
function igyp02b()
#
#  generalization of Igarash and Ypma
#
    igyp02(7)
end


text:
__precompile__()

module Watcher
include("code.jl")

###########################
# start the watchers

function __init__()
    try
        (dirs, filetypes, cmd, runnow) = parseargs()

        files = filter(x -> any(map(y -> endswith(x,y), filetypes)), walk().files)

        processes = Any[nothing]
        watchers = Dict()

        f(h) = (runcmd(cmd, processes, watchers, h); watchfiles(f, files, watchers))
        watchfiles(f, files, watchers)

        if runnow
            runcmd(cmd, processes)
        end

        while true
            sleep(1)
        end

    catch e
        if isa(e, InterruptException)
            exit()
        else
            rethrow(e)
        end
    end
end

end # module


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

# generic operations on associative collections

const secret_table_token = :__c782dbf1cf4d6a2e5e3865d7e95634f2e09b5902__

haskey(d::Associative, k) = in(k,keys(d))

function in(p::Pair, a::Associative, valcmp=(==))
    v = get(a,p[1],secret_table_token)
    if v !== secret_table_token
        valcmp(v, p[2]) && return true
    end
    return false
end

function in(p, a::Associative)
    error("""Associative collections only contain Pairs;
             Either look for e.g. A=>B instead, or use the `keys` or `values`
             function if you are looking for a key or value respectively.""")
end

function summary(t::Associative)
    n = length(t)
    return string(typeof(t), " with ", n, (n==1 ? " entry" : " entries"))
end

struct KeyIterator{T<:Associative}
    dict::T
end
struct ValueIterator{T<:Associative}
    dict::T
end

summary(iter::T) where {T<:Union{KeyIterator,ValueIterator}} =
    string(T.name, " for a ", summary(iter.dict))

show(io::IO, iter::Union{KeyIterator,ValueIterator}) = show(io, collect(iter))

length(v::Union{KeyIterator,ValueIterator}) = length(v.dict)
isempty(v::Union{KeyIterator,ValueIterator}) = isempty(v.dict)
_tt1(::Type{Pair{A,B}}) where {A,B} = A
_tt2(::Type{Pair{A,B}}) where {A,B} = B
eltype(::Type{KeyIterator{D}}) where {D} = _tt1(eltype(D))
eltype(::Type{ValueIterator{D}}) where {D} = _tt2(eltype(D))

start(v::Union{KeyIterator,ValueIterator}) = start(v.dict)
done(v::Union{KeyIterator,ValueIterator}, state) = done(v.dict, state)

function next(v::KeyIterator, state)
    n = next(v.dict, state)
    n[1][1], n[2]
end

function next(v::ValueIterator, state)
    n = next(v.dict, state)
    n[1][2], n[2]
end

in(k, v::KeyIterator) = get(v.dict, k, secret_table_token) !== secret_table_token

"""
    keys(a::Associative)

Return an iterator over all keys in a collection.
`collect(keys(a))` returns an array of keys.
Since the keys are stored internally in a hash table,
the order in which they are returned may vary.
But `keys(a)` and `values(a)` both iterate `a` and
return the elements in the same order.

```jldoctest
julia> a = Dict('a'=>2, 'b'=>3)
Dict{Char,Int64} with 2 entries:
  'b' => 3
  'a' => 2

julia> collect(keys(a))
2-element Array{Char,1}:
 'b'
 'a'
```
"""
keys(a::Associative) = KeyIterator(a)
eachindex(a::Associative) = KeyIterator(a)

"""
    values(a::Associative)

Return an iterator over all values in a collection.
`collect(values(a))` returns an array of values.
Since the values are stored internally in a hash table,
the order in which they are returned may vary.
But `keys(a)` and `values(a)` both iterate `a` and
return the elements in the same order.

```jldoctest
julia> a = Dict('a'=>2, 'b'=>3)
Dict{Char,Int64} with 2 entries:
  'b' => 3
  'a' => 2

julia> collect(values(a))
2-element Array{Int64,1}:
 3
 2
```
"""
values(a::Associative) = ValueIterator(a)

function copy(a::Associative)
    b = similar(a)
    for (k,v) in a
        b[k] = v
    end
    return b
end

"""
    merge!(d::Associative, others::Associative...)

Update collection with pairs from the other collections.
See also [`merge`](@ref).

```jldoctest
julia> d1 = Dict(1 => 2, 3 => 4);

julia> d2 = Dict(1 => 4, 4 => 5);

julia> merge!(d1, d2);

julia> d1
Dict{Int64,Int64} with 3 entries:
  4 => 5
  3 => 4
  1 => 4
```
"""
function merge!(d::Associative, others::Associative...)
    for other in others
        for (k,v) in other
            d[k] = v
        end
    end
    return d
end

"""
    merge!(combine, d::Associative, others::Associative...)

Update collection with pairs from the other collections.
Values with the same key will be combined using the
combiner function.

```jldoctest
julia> d1 = Dict(1 => 2, 3 => 4);

julia> d2 = Dict(1 => 4, 4 => 5);

julia> merge!(+, d1, d2);

julia> d1
Dict{Int64,Int64} with 3 entries:
  4 => 5
  3 => 4
  1 => 6

julia> merge!(-, d1, d1);

julia> d1
Dict{Int64,Int64} with 3 entries:
  4 => 0
  3 => 0
  1 => 0
```
"""
function merge!(combine::Function, d::Associative, others::Associative...)
    for other in others
        for (k,v) in other
            d[k] = haskey(d, k) ? combine(d[k], v) : v
        end
    end
    return d
end

# very similar to `merge!`, but accepts any iterable and extends code
# that would otherwise only use `copy!` with arrays.
function copy!(dest::Union{Associative,AbstractSet}, src)
    for x in src
        push!(dest, x)
    end
    return dest
end

"""
    keytype(type)

Get the key type of an associative collection type. Behaves similarly to [`eltype`](@ref).

```jldoctest
julia> keytype(Dict(Int32(1) => "foo"))
Int32
```
"""
keytype(::Type{Associative{K,V}}) where {K,V} = K
keytype(a::Associative) = keytype(typeof(a))
keytype(::Type{A}) where {A<:Associative} = keytype(supertype(A))

"""
    valtype(type)

Get the value type of an associative collection type. Behaves similarly to [`eltype`](@ref).

```jldoctest
julia> valtype(Dict(Int32(1) => "foo"))
String
```
"""
valtype(::Type{Associative{K,V}}) where {K,V} = V
valtype(::Type{A}) where {A<:Associative} = valtype(supertype(A))
valtype(a::Associative) = valtype(typeof(a))

"""
    merge(d::Associative, others::Associative...)

Construct a merged collection from the given collections. If necessary, the
types of the resulting collection will be promoted to accommodate the types of
the merged collections. If the same key is present in another collection, the
value for that key will be the value it has in the last collection listed.

```jldoctest
julia> a = Dict("foo" => 0.0, "bar" => 42.0)
Dict{String,Float64} with 2 entries:
  "bar" => 42.0
  "foo" => 0.0

julia> b = Dict("baz" => 17, "bar" => 4711)
Dict{String,Int64} with 2 entries:
  "bar" => 4711
  "baz" => 17

julia> merge(a, b)
Dict{String,Float64} with 3 entries:
  "bar" => 4711.0
  "baz" => 17.0
  "foo" => 0.0

julia> merge(b, a)
Dict{String,Float64} with 3 entries:
  "bar" => 42.0
  "baz" => 17.0
  "foo" => 0.0
```
"""
merge(d::Associative, others::Associative...) =
    merge!(emptymergedict(d, others...), d, others...)

"""
    merge(combine, d::Associative, others::Associative...)

Construct a merged collection from the given collections. If necessary, the
types of the resulting collection will be promoted to accommodate the types of
the merged collections. Values with the same key will be combined using the
combiner function.

```jldoctest
julia> a = Dict("foo" => 0.0, "bar" => 42.0)
Dict{String,Float64} with 2 entries:
  "bar" => 42.0
  "foo" => 0.0

julia> b = Dict("baz" => 17, "bar" => 4711)
Dict{String,Int64} with 2 entries:
  "bar" => 4711
  "baz" => 17

julia> merge(+, a, b)
Dict{String,Float64} with 3 entries:
  "bar" => 4753.0
  "baz" => 17.0
  "foo" => 0.0
```
"""
merge(combine::Function, d::Associative, others::Associative...) =
    merge!(combine, emptymergedict(d, others...), d, others...)

promoteK(K) = K
promoteV(V) = V
promoteK(K, d, ds...) = promoteK(promote_type(K, keytype(d)), ds...)
promoteV(V, d, ds...) = promoteV(promote_type(V, valtype(d)), ds...)
function emptymergedict(d::Associative, others::Associative...)
    K = promoteK(keytype(d), others...)
    V = promoteV(valtype(d), others...)
    Dict{K,V}()
end

function filter!(f, d::Associative)
    badkeys = Vector{keytype(d)}(0)
    for (k,v) in d
        # don't delete!(d, k) here, since associative types
        # may not support mutation during iteration
        f(k,v) || push!(badkeys, k)
    end
    for k in badkeys
        delete!(d, k)
    end
    return d
end
function filter(f, d::Associative)
    # don't just do filter!(f, copy(d)): avoid making a whole copy of d
    df = similar(d)
    for (k,v) in d
        if f(k,v)
            df[k] = v
        end
    end
    return df
end

eltype(::Type{Associative{K,V}}) where {K,V} = Pair{K,V}

function isequal(l::Associative, r::Associative)
    l === r && return true
    if isa(l,ObjectIdDict) != isa(r,ObjectIdDict)
        return false
    end
    if length(l) != length(r) return false end
    for pair in l
        if !in(pair, r, isequal)
            return false
        end
    end
    true
end

function ==(l::Associative, r::Associative)
    l === r && return true
    if isa(l,ObjectIdDict) != isa(r,ObjectIdDict)
        return false
    end
    if length(l) != length(r) return false end
    for pair in l
        if !in(pair, r, ==)
            return false
        end
    end
    true
end

const hasha_seed = UInt === UInt64 ? 0x6d35bb51952d5539 : 0x952d5539
function hash(a::Associative, h::UInt)
    hv = hasha_seed
    for (k,v) in a
        hv ⊻= hash(k, hash(v))
    end
    hash(hv, h)
end

function getindex(t::Associative, key)
    v = get(t, key, secret_table_token)
    if v === secret_table_token
        throw(KeyError(key))
    end
    return v
end

# t[k1,k2,ks...] is syntactic sugar for t[(k1,k2,ks...)].  (Note
# that we need to avoid dispatch loops if setindex!(t,v,k) is not defined.)
getindex(t::Associative, k1, k2, ks...) = getindex(t, tuple(k1,k2,ks...))
setindex!(t::Associative, v, k1, k2, ks...) = setindex!(t, v, tuple(k1,k2,ks...))

push!(t::Associative, p::Pair) = setindex!(t, p.second, p.first)
push!(t::Associative, p::Pair, q::Pair) = push!(push!(t, p), q)
push!(t::Associative, p::Pair, q::Pair, r::Pair...) = push!(push!(push!(t, p), q), r...)

# hashing objects by identity

"""
    ObjectIdDict([itr])

`ObjectIdDict()` constructs a hash table where the keys are (always)
object identities.  Unlike `Dict` it is not parameterized on its key
and value type and thus its `eltype` is always `Pair{Any,Any}`.

See [`Dict`](@ref) for further help.
"""
mutable struct ObjectIdDict <: Associative{Any,Any}
    ht::Vector{Any}
    ndel::Int
    ObjectIdDict() = new(Vector{Any}(32), 0)

    function ObjectIdDict(itr)
        d = ObjectIdDict()
        for (k,v) in itr; d[k] = v; end
        d
    end

    function ObjectIdDict(pairs::Pair...)
        d = ObjectIdDict()
        for (k,v) in pairs; d[k] = v; end
        d
    end

    ObjectIdDict(o::ObjectIdDict) = new(copy(o.ht))
end

similar(d::ObjectIdDict) = ObjectIdDict()

function rehash!(t::ObjectIdDict, newsz = length(t.ht))
    t.ht = ccall(:jl_idtable_rehash, Any, (Any, Csize_t), t.ht, newsz)
    t
end

function setindex!(t::ObjectIdDict, v::ANY, k::ANY)
    if t.ndel >= ((3*length(t.ht))>>2)
        rehash!(t, max(length(t.ht)>>1, 32))
        t.ndel = 0
    end
    t.ht = ccall(:jl_eqtable_put, Array{Any,1}, (Any, Any, Any), t.ht, k, v)
    return t
end

get(t::ObjectIdDict, key::ANY, default::ANY) =
    ccall(:jl_eqtable_get, Any, (Any, Any, Any), t.ht, key, default)

function pop!(t::ObjectIdDict, key::ANY, default::ANY)
    val = ccall(:jl_eqtable_pop, Any, (Any, Any, Any), t.ht, key, default)
    # TODO: this can underestimate `ndel`
    val === default || (t.ndel += 1)
    return val
end

function pop!(t::ObjectIdDict, key::ANY)
    val = pop!(t, key, secret_table_token)
    val !== secret_table_token ? val : throw(KeyError(key))
end

function delete!(t::ObjectIdDict, key::ANY)
    pop!(t, key, secret_table_token)
    t
end

function empty!(t::ObjectIdDict)
    resize!(t.ht, 32)
    ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), t.ht, 0, sizeof(t.ht))
    t.ndel = 0
    return t
end

_oidd_nextind(a, i) = reinterpret(Int,ccall(:jl_eqtable_nextind, Csize_t, (Any, Csize_t), a, i))

start(t::ObjectIdDict) = _oidd_nextind(t.ht, 0)
done(t::ObjectIdDict, i) = (i == -1)
next(t::ObjectIdDict, i) = (Pair{Any,Any}(t.ht[i+1],t.ht[i+2]), _oidd_nextind(t.ht, i+2))

function length(d::ObjectIdDict)
    n = 0
    for pair in d
        n+=1
    end
    n
end

copy(o::ObjectIdDict) = ObjectIdDict(o)

get!(o::ObjectIdDict, key, default) = (o[key] = get(o, key, default))

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia