text:
module XInts

export XInt

using Base.GMP: Limb, BITS_PER_LIMB, SLimbMax, ULimbMax, ClongMax, CulongMax, CdoubleMax
import Base.GMP: isneg, ispos
import Base.GMP.MPZ
import Base.MPFR
using Base.GC: @preserve
using Base: BitInteger
import Base: +, -, *, ^, &, |, ==, /, ~, <<, >>, >>>, <, <=,
             string, widen, hastypemax, tryparse_internal,
             unsafe_trunc, trunc, mod, rem, iseven, isodd, gcd, lcm, xor, div, fld, cld,
             invmod, count_ones, trailing_zeros, trailing_ones, cmp, isqrt,
             flipsign, powermod, gcdx, promote_rule, factorial, binomial,
             digits!, ndigits0zpb, signbit, sign, iszero, isone,
             AbstractFloat, BigFloat, float, _prevpow2, _nextpow2,
             hash, hash_integer, sum, divrem

mutable struct Wrap
    b::BigInt
    d::Ptr{Limb}

    function Wrap(isoutput::Bool=false)
        b = BigInt()
        @assert b.alloc == 1
        @assert b.size == 0
        wrap = new(b, b.d)
        if !isoutput
            # init! was called
            finalizer(wrap) do wrap
                wrap.b.d = wrap.d
                wrap.b.alloc = 1
                wrap.b.size = 0
            end
        end
        wrap
    end
end

if BITS_PER_LIMB == 32
    const SLimb = Int32
    const SLimbW = Int64
    const LimbW = UInt64

    @inline _divLimb(n) = n >>> 5
    @inline _modLimb(n) = n & 31
    @inline _divByteLimb(n) = n >>> 2
    @inline _modByteLimb(n) = n & 3
elseif BITS_PER_LIMB == 64
    const SLimb = Int64
    const SLimbW = Int128
    const LimbW = UInt128

    @inline _divLimb(n) = n >>> 6
    @inline _modLimb(n) = n & 63
    @inline _divByteLimb(n) = n >>> 3
    @inline _modByteLimb(n) = n & 7
else
    error()
end

const slimbmin = typemin(SLimb)
const slimbmax = typemax(SLimb)
# when there is only one limb in a vector, this limb must always be >= limb1min
# reciprocally, a direct integer must have its absolute value < limb1min
const limb1min = slimbmin % Limb # typically 0x8000000000000000

const Short = Union{Limb,SLimb}
const ShortMax = Union{ULimbMax,SLimbMax}

abstract type XSigned <: Signed end

# T must be a "bit integer", assumed for efficiency, and also sizeof
struct Static{T<:Integer} <: XSigned
    x::T
end

# we could re-use Static for this, but Static <: Number already has getindex defined,
# so we clearly separate the two concepts
struct SVec{U<:Unsigned}
    u::U
end

is_short(s::Static) = sizeof(s.x) < sizeof(Limb) ? true : fits(s.x)
short(s::Static) = s.x % SLimb

vec(s::Static{<:Signed}) = SVec(unsigned(abs(s.x)))
vec(s::Static{<:Unsigned}) = SVec(s.x)

function lenvec(s::Static)
    # @assert !is_short(s)
    v = vec(s)
    lz = leading_zeros(v.u) >>> 3 # count leading full zero bytes
    usedbytes = sizeof(v.u)-lz
    _divByteLimb(usedbytes - 1) + 1, v
end

len(s::Static) = lenvec(s)[1]
vec(s::Static) = lenvec(s)[2]

iszero(s::Static) = iszero(s.x)
ispos(s::Static) = s.x > 0
signbit(s::Static) = s.x < 0
_sign(s::Static) = ispos(s) ? SLimb(1) : SLimb(-1)

function Base.getindex(s::SVec, i::Integer)
    # getindex is called only when sizeof(x) >= sizeof(Limb)
    if sizeof(s.u) == sizeof(Limb)
        # we always assume @inbounds
        s.u % Limb
    else
        (s.u >> ((i-1)*BITS_PER_LIMB)) % Limb
    end
end

Base.show(io::IO, s::Static) =
    if is_short(s)
        print(io, "Static($(s.x))")
    else
        sl, sv = lenvec(s)
        print(io, "Static(", [sv[i] for i=1:sl], ")")
    end

struct XInt <: XSigned
    x::SLimb # immediate integer or sign+length
    v::Union{Nothing,Vector{Limb}}

    @inline XInt(x::SLimb) =
        x === typemin(SLimb) ?
            new(-one(SLimb), [limb1min]) :
            new(x, nothing)

    global _XInt

    # unsafe version which doesn't check for typemin(SLimb)
    _XInt(x::SLimb, ::Nothing=nothing) = new(x, nothing)

    function _XInt(x::SLimb, v::Vector{Limb}, reduce::Bool=false)
        if reduce
            # we still assume length(v) >= abs(x); reduce here means to not store
            # a too small integer in a vector representation
            xl = abs(x)
            if xl <= 1
                if xl == 0
                    return new(SLimb(0), nothing)
                else # xl == 1
                    limb = @inbounds v[1]
                    limb < limb1min &&
                        return new(flipsign(limb % SLimb, x), nothing)
                end
            end
        end
        new(x, v)
    end
end

XInt(x::XInt) = x
XInt(x::Static) = XInt(x.x)

# reducing version of XInt(::XInt)
_XInt(x::XInt) =
    if x.v === nothing
        x
    else
        _XInt(x.x, vec(x), true)
    end

_XInt(u::Limb) = _XInt(u % SLimb) # no check done that u fits

is_short(x::XInt) = x.v === nothing
is_short(x::XInt, y::XInt) = x.v === nothing === y.v

short(x::XInt) = x.x
len(x::XInt) = abs(x.x)

include("MPN.jl")
include("mutants.jl")

MPN.ptr(x::XInt) = pointer(x.v::Vector{Limb})
MPN.len(x::XInt) = abs(x.x) % MPN.Size

function _XInt(z::BigInt, extra::SLimb=0)
    # @assert extra >= 0
    len = abs(z.size)
    if len == 0
        XInt(zero(SLimb))
    elseif len == 1 && (x = unsafe_load(z.d)) < limb1min
        XInt(flipsign(x % SLimb, z.size))
    else
        x = _XInt(z.size % SLimb, Vector{Limb}(undef, len+extra))
        @preserve x z unsafe_copyto!(pointer(x.v), z.d, len)
        x
    end
end

XInt(z::BigInt) = _XInt(z)

fits(z::Limb) = !Core.is_top_bit_set(z) # z < limb1min
fits(z::SLimb) = z !== slimbmin
fits(z::Signed) = slimbmin < z <= slimbmax
fits(z::Unsigned) = z <= slimbmax

XInt(z::Limb) = XInt!(nothing, z)

# only SLimb and Limb in this Union need to be validated, but are taken care of by
# other more specific constructors
XInt(z::Union{ULimbMax,SLimbMax}) = _XInt(z % SLimb)

XInt(z::Union{SLimbW,LimbW}) = XInt!(nothing, z)
_XInt(z::Union{SLimbW,LimbW}, extra::SLimb=0) = XInt!(nothing, z, extra)

# TODO: this can be made quite faster for "bit integers": for some reason, allocating
# limbs according to sizeof(z) (and breaking the loop when iszero(zz))
# is more efficient than via ndigits (although ndigits seems to be fast)
# NOTE: also, it allocates more than expected for e.g. Int256, and this disappears
# if the function is re-evaled (then we get 1 allocation as expected)
function _XInt(z::Integer, extra::SLimb=0)::XInt
    slimbmin < z <= slimbmax &&
        return _XInt((z % SLimb)::SLimb)
    nd = ndigits(z, base=2)
    xl = _divLimb(nd)
    xl += !iszero(_modLimb(nd))
    limbs = _vec(xl+extra)
    zz = abs(z)
    for i=1:xl
        @inbounds limbs[i] = zz % Limb
        zz >>>= BITS_PER_LIMB
    end
    _XInt(flipsign(xl, z), limbs)
end

@inline XInt(z::Integer) = _XInt(z)

function XInt(x::Float64)
    isinteger(x) || throw(InexactError(:XInt, XInt, x))
    unsafe_trunc(XInt, x)
end

# no Union, to disambiguate
XInt(x::Float16) = XInt(Float64(x))
XInt(x::Float32) = XInt(Float64(x))

function init!(z::Wrap, x::XInt)
    b = z.b
    if is_short(x)
        b.alloc = 1
        b.size = sign(x.x)
        b.d = z.d
        unsafe_store!(z.d, abs(x.x) % Limb)
    else
        len = x.x
        b.alloc = abs(len)
        b.size = len
        b.d = pointer(x.v)
    end
    z.b::BigInt
end

struct Wraps
    # read-write temporary storage (output arguments)
    U::Wrap
    V::Wrap
    W::Wrap
    # read-only shallow wrappers of BigInt objects ("const" arguments in libgmp)
    X::Wrap
    Y::Wrap
    Z::Wrap

    Wraps() = new(Wrap(true), Wrap(true), Wrap(true), Wrap(), Wrap(), Wrap())
end

const tWraps = Wraps[]

mutable struct Buffer
    vec::Vector{Limb}
    # TODO: determine whether `used` is really needed; it's to try to work
    # around the problem mentioned in Julia#17015, that affects also usage
    # of the @bigint macro, related to GC
    used::Bool

    Buffer() = new(sizehint!(Limb[], 20), false)
end

const tBuffer = Buffer[]

function tmpvec(n::Integer)
    buf = tBuffer[Threads.threadid()]
    if buf.used
        # if a previous call failed to be followed by freetmpvec,
        # or maybe finalizers triggered by GC call tmpvec while the
        # resource is already in use
        _vec(n)
    else
        v = vec!(buf.vec, n)
        buf.used = true # after possible allocation above, which might fail
        v
    end
end

function freetmpvec(v::Vector)
    buf = tBuffer[Threads.threadid()]
    if v === buf.vec
        @assert buf.used
        buf.used = false
    end
    nothing
end

function __init__()
    copy!(tWraps, [Wraps() for _=1:Threads.nthreads()])
    copy!(tBuffer, [Buffer() for _=1:Threads.nthreads()])
end

macro bigint(args...)
    @assert length(args) >= 2
    out = args[1]
    xs = esc.(args[2:end-1])
    body = args[end]

    lets = :(let
                 $(esc(body))
             end)
    vars = lets.args[1].args

    @gensym mpzs
    push!(vars, :($(esc(mpzs)) = tWraps[Threads.threadid()]))
    if !(out isa Expr)
        out = Expr(:tuple, out)
    end
    for (i, outi) in enumerate(out.args)
        push!(vars, :($(esc(outi)) = getfield($(esc(mpzs)), $i).b))
    end
    for (i, xi) in enumerate(xs)
        push!(vars, :($xi = init!(getfield($(esc(mpzs)), $i+3), $xi)))
    end
    Expr(:gc_preserve, lets, xs...)
end

function BigInt(x::XInt)
    if is_short(x)
        BigInt(x.x)
    else
        xl, xv = lenvec(x)
        nbits = xl*8*sizeof(Limb)
        b = VERSION >= v"1.3" ? BigInt(; nbits=nbits) :
                                MPZ.realloc2(nbits)
        @preserve b xv unsafe_copyto!(b.d, pointer(xv), xl)
        b.size = x.x
        b
    end
end

Base.show(io::IO, x::XInt) = print(io, string(x))

function string(n::XInt; base::Integer = 10, pad::Integer = 1)
    is_short(n) && return string(n.x; base=base, pad=pad)
    base < 0 && return Base._base(Int(base), n, pad, (base>0) & isneg(n))
    2 <= base <= 62 || throw(ArgumentError("base must be 2 ≤ base ≤ 62, got $base"))
    nd1 = ndigits(n, base=base)
    nd  = max(nd1, pad)
    sv  = Base.StringVector(nd + isneg(n))
    xl, xv = lenvec(n)
    if !ispow2(base)
        # MPN.get_str might clobber the input
        xv = copy(xv)
    end
    str_size = @preserve sv Int(MPN.get_str(pointer(sv) + nd - nd1 + isneg(n), base, xv, xl))
    @assert str_size == nd1
    digits = base <= 36 ? Base.base36digits : Base.base62digits
    @inbounds for i = (1:nd-nd1) .+ isneg(n)
        sv[i] = '0' % UInt8
    end
    @inbounds for j = nd-nd1+isneg(n)+1:lastindex(sv)
        sv[j] = digits[1+sv[j]]
    end
    isneg(n) && (sv[1] = '-' % UInt8)
    String(sv)
end

widen(::Type{XInt}) = XInt

hastypemax(::Type{XInt}) = false

promote_rule(::Type{XInt}, ::Type{<:Integer}) = XInt
promote_rule(::Type{XInt}, ::Type{BigInt}) = XInt
promote_rule(::Type{BigInt}, ::Type{XInt}) = XInt
promote_rule(::Type{XInt}, ::Type{<:AbstractFloat}) = BigFloat

tryparse_internal(::Type{XInt}, s::AbstractString, startpos::Int, endpos::Int,
                  base_::Integer, raise::Bool) =
                      XInt(tryparse_internal(BigInt, s, startpos, endpos, base_, raise))

unsafe_trunc(::Type{XInt}, x::Union{Float32,Float64}) = XInt(unsafe_trunc(BigInt, x))

function trunc(::Type{XInt}, x::Union{Float32,Float64})
    isfinite(x) || throw(InexactError(:trunc, XInt, x))
    unsafe_trunc(XInt, x)
end

rem(x::XInt, ::Type{Bool}) =
    is_short(x) ? x.x % Bool :
                  (@inbounds x.v[1]) % Bool

rem(x::XInt, ::Type{T}) where T<:Union{SLimbMax,ULimbMax} =
    is_short(x) ? x.x % T : flipsign((@inbounds x.v[1] % T), x.x)

# TODO: extend to "AbstractBitInteger", when available in Base
function rem(x::XInt, ::Type{T}) where T<:Union{Base.BitUnsigned,Base.BitSigned}
    is_short(x) && return x.x % T
    u = zero(T)
    # TODO: this is missing an optimization that BigInt has; it would be faster here
    # to have a chain of if/elseif based on the value of the min below
    for l = 1:min(abs(x.x), cld(sizeof(T), sizeof(Limb)))
        u += (@inbounds x.v[l] % T) << ((sizeof(Limb)<<3)*(l-1))
    end
    flipsign(u, x.x)
end

rem(x::Integer, ::Type{XInt}) = XInt(x)

isodd(x::XInt) = is_short(x) ? isodd(x.x) : isodd(@inbounds x.v[1])
iseven(x::XInt) = !isodd(x)

function (::Type{T})(x::XInt) where T<:Base.BitUnsigned
    if is_short(x)
        convert(T, x.x)
    elseif sizeof(T) < sizeof(Limb)
        convert(T, convert(Limb, x))
    else
        0 <= x.x <= cld(sizeof(T), sizeof(Limb)) || throw(InexactError(nameof(T), T, x))
        x % T
    end
end

ispos(x::XInt) = x.x > 0
ispos(x::XInt, y::XInt) = x.x > 0 && y.x > 0
isneg(x::XInt) = x.x < 0
signbit(x::XInt) = x.x < 0
sign(x::XInt) = _XInt(sign(x.x))
_sign(x::XInt) = sign(x.x)
iszero(x::XInt) = iszero(x.x)
isone(x::XInt) = is_short(x) ? isone(x.x) : isone(@inbounds x.v[1])
# TODO: remove is_short check when we guarantee that short integers are always stored in x.x

function (::Type{T})(x::XInt) where T<:Base.BitSigned
    is_short(x) && return T(x.x)
    n = abs(x.x)
    if sizeof(T) < sizeof(Limb)
        # @assert SLimb == typeof(Signed(one(Limb)))
        convert(T, convert(SLimb, x))
    else
        0 <= n <= cld(sizeof(T), sizeof(Limb)) || throw(InexactError(nameof(T), T, x))
        y = x % T
        ispos(x) ⊻ (y > 0) && throw(InexactError(nameof(T), T, x)) # catch overflow
        y
    end
end

function Float64(x::XInt, ::RoundingMode{:ToZero})
    if is_short(x)
        Float64(x.x, RoundToZero)
    else
        MPN.get_d(x, len(x), x.x)
    end
end

function (::Type{T})(x::XInt, ::RoundingMode{:ToZero}) where T<:Union{Float16,Float32}
    T(Float64(x, RoundToZero), RoundToZero)
end

function (::Type{T})(n::XInt, ::RoundingMode{:Down}) where T<:CdoubleMax
    x = T(n, RoundToZero)
    x > n ? prevfloat(x) : x
end

function (::Type{T})(n::XInt, ::RoundingMode{:Up}) where T<:CdoubleMax
    x = T(n, RoundToZero)
    x < n ? nextfloat(x) : x
end

for T = [Float16, Float32, Float64]
    @eval (::Type{$T})(x::XInt, r::RoundingMode{:Nearest}=RoundNearest) =
        is_short(x) ? $T(x.x, r) : @bigint () x $T(x, r)
end

AbstractFloat(x::XInt) = BigFloat(x)

BigFloat(x::XInt, r::MPFR.MPFRRoundingMode=MPFR.ROUNDING_MODE[];
         precision::Integer=MPFR.DEFAULT_PRECISION[]) =
             is_short(x) ? BigFloat(x.x, r; precision=precision) :
                           @bigint () x BigFloat(x, r; precision=precision)

float(::Type{XInt}) = BigFloat

# Binary ops
for (fJ, fC) in ((:mod, :fdiv_r), (:rem, :tdiv_r),
                 (:gcd, :gcd), (:lcm, :lcm))
    fC! = Symbol(fC, :!)
    @eval begin
        ($fJ)(x::XInt, y::XInt) =
            is_short(x, y) ? XInt(($fJ)(widen(x.x), widen(y.x))) :
                             @bigint z x y XInt(MPZ.$fC!(z, x, y))
    end
end

# TODO: 3+ args specializations for some ops, like in gmp.jl
# (already done for +)

# we put @inline, as it seems these methods don't by themselves,
# maybe because add! is already inline and too big
@inline +(x::XInt, y::XInt)       = add!(nothing, x, y)
@inline +(x::XInt, y::BitInteger) = add!(nothing, x, y)
@inline +(x::XInt, y::Integer)    = add!(_XInt(y, 1), x)
@inline +(x::Integer, y::XInt)    = y + x

@inline -(x::XInt, y::XInt)       = sub!(nothing, x, y)
@inline -(x::BitInteger, y::XInt) = sub!(nothing, x, y)
@inline -(x::Integer, y::XInt)    = sub!(_XInt(x, 1), y)
@inline -(x::XInt, y::Integer)    = neg!(y-x)

@inline *(x::XInt, y::XInt)     = mul!(nothing, x, y)
@inline *(x::XInt, c::ShortMax) = mul!(nothing, x, c)
@inline *(x::XInt, y::Integer)  = mul!(_XInt(y, safe_len(x)), x)
@inline *(x::Integer, y::XInt)  = y * x

@inline (&)(x::XInt, y::XInt)       = and!(nothing, x, y)
@inline (&)(x::XInt, y::BitInteger) = and!(nothing, x, y)
@inline (&)(x::XInt, y::Integer)    = and!(XInt(y), x)
@inline (&)(x::Integer, y::XInt)    = y & x

@inline (|)(x::XInt, y::XInt)       = ior!(nothing, x, y)
@inline (|)(x::XInt, y::BitInteger) = ior!(nothing, x, y)
@inline (|)(x::XInt, y::Integer)    = ior!(XInt(y), x)
@inline (|)(x::Integer, y::XInt)    = y | x

@inline xor(x::XInt, y::XInt)       = xor!(nothing, x, y)
@inline xor(x::XInt, y::BitInteger) = xor!(nothing, x, y)
@inline xor(x::XInt, y::Integer)    = xor!(XInt(y), x)
@inline xor(x::Integer, y::XInt)    = y ⊻ x

# it doesn't seem worthwhile to define 4+ args specific + methods,
# the vararg one is better
@inline +(x::XInt, y::XInt, z::XInt) = add!((x+y), z)
+(x::XInt, y::XInt, z::XInt, r::XInt...) = sum(tuple(x, y, z, r...))

function sum(arr::Union{AbstractArray{XInt},Tuple{Vararg{XInt}}})
    sp = _XInt(0) # > 0
    sn = _XInt(0) # < 0
    s0 = zero(Int128)

    for x = arr
        if is_short(x)
            s0 += short(x)
        elseif ispos(x)
            sp = addbig_positive!(sp, x)
        else
            sn = addbig_positive!(sn, x)
        end
    end
    if !iszero(sp)
        if !iszero(sn)
            s = addbig!(sp, sp, neg!(sn))
            add!(s, s0)
        else
            add!(sp, s0)
        end
    elseif !iszero(sn)
        add!(neg!(sn), s0)
    else
        XInt(s0)
    end
end

for (r, f!) in ((RoundToZero, :tdiv_q!),
               (RoundDown, :fdiv_q!),
               (RoundUp, :cdiv_q!))
    @eval div(x::XInt, y::XInt, ::typeof($r)) =
        is_short(x, y) ? XInt(div(widen(x.x), widen(y.x), $r)) :
                         @bigint z x y XInt(MPZ.$f!(z, x, y))
end

# For compat only. Remove in 2.0.
div(x::XInt, y::XInt) = div(x, y, RoundToZero)
fld(x::XInt, y::XInt) = div(x, y, RoundDown)
cld(x::XInt, y::XInt) = div(x, y, RoundUp)

divrem(x::XInt, y::XInt) =
    is_short(x, y) ? XInt.(divrem(x.x, y.x)) :
                     @bigint (u, v) x y XInt.(MPZ.tdiv_qr!(u, v, x, y))

/(x::XInt, y::XInt) = float(x)/float(y)

invmod(x::XInt, y::XInt) =
    is_short(x, y) ? XInt(invmod(widen(x.x), widen(y.x))) :
                     @bigint () x y XInt(invmod(x, y))

/(x::XInt, y::Union{ClongMax,CulongMax}) = float(x)/y
/(x::Union{ClongMax,CulongMax}, y::XInt) = x/float(y)

# unary ops
(-)(x::XInt) = neg!(x)
(~)(x::XInt) = com!(nothing, x)

<<(x::XInt, c::UInt) = lshift!(nothing, x, c)
>>(x::XInt, c::UInt) = rshift!(nothing, x, c)
>>>(x::XInt, c::UInt) = x >> c

trailing_zeros(x::XInt) =
    if is_short(x)
        iszero(x.x) && throw(DomainError(x)) # InexactError for BigInt
        trailing_zeros(x.x)
    else
        @bigint () x MPZ.scan1(x, 0)
    end

trailing_ones(x::XInt) =
    if is_short(x)
        x.x == -1 && throw(DomainError(x))
        trailing_ones(x.x)
    else
        @bigint () x MPZ.scan0(x, 0)
    end

count_ones(x::XInt) =
    if is_short(x)
        x.x < 0 && throw(DomainError(x))
        count_ones(x.x)
    else
        @bigint () x MPZ.popcount(x)
    end

cmp(x::XInt, y::XInt) = is_short(x, y) ? cmp(x.x, y.x) : @bigint () x y cmp(x, y)
cmp(x::XInt, y::Integer) = is_short(x) ? cmp(x.x, y) : @bigint () x cmp(x, y)
cmp(x::XInt, y::BigInt)  = is_short(x) ? cmp(x.x, y) : @bigint () x cmp(x, y) # disambiguation
cmp(y::Integer, x::XInt) = -cmp(x, y)
cmp(y::BigInt, x::XInt)  = -cmp(x, y) # disambiguation
cmp(x::XInt, y::CdoubleMax) = is_short(x) ? cmp(x.x, y) : @bigint () x cmp(x, y)
cmp(y::CdoubleMax, x::XInt) = -cmp(x, y)

function ==(x::XInt, y::XInt)
    if is_short(x, y)
        x.x == y.x
    else
        @bigint () x y x == y
    end
end

==(x::XInt, y::Integer) =
    if is_short(x)
        x.x == y
    else
        @bigint () x x == y
    end

==(x::Integer, y::XInt) = y == x

# disambiguation
==(x::XInt, y::BigInt) = invoke(==, Tuple{XInt, Integer}, x, y)
==(x::BigInt, y::XInt) = y == x

==(x::XInt, f::CdoubleMax) = isnan(f) ? false : cmp(x, f) == 0
==(f::CdoubleMax, x::XInt) = isnan(f) ? false : cmp(x, f) == 0

<=(x::XInt, y::XInt) = cmp(x,y) <= 0
<=(x::XInt, i::Integer) = cmp(x,i) <= 0
<=(x::XInt, i::BigInt)  = cmp(x,i) <= 0 # disambiguation
<=(i::Integer, x::XInt) = cmp(x,i) >= 0
<=(i::BigInt, x::XInt)  = cmp(x,i) >= 0 # disambiguation
<=(x::XInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) <= 0
<=(f::CdoubleMax, x::XInt) = isnan(f) ? false : cmp(x,f) >= 0

<(x::XInt, y::XInt) = cmp(x,y) < 0
<(x::XInt, i::Integer) = cmp(x,i) < 0
<(x::XInt, i::BigInt)  = cmp(x,i) < 0 # disambiguation
<(i::Integer, x::XInt) = cmp(x,i) > 0
<(i::BigInt, x::XInt)  = cmp(x,i) > 0 # disambiguation

<(x::XInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) < 0
<(f::CdoubleMax, x::XInt) = isnan(f) ? false : cmp(x,f) > 0

isqrt(x::XInt) = is_short(x) ? XInt(isqrt(x.x)) : @bigint z x XInt(MPZ.sqrt!(z, x))

^(x::XInt, y::Culong) = @bigint z x XInt(MPZ.pow_ui!(z, x, y))

function xint_pow(x::XInt, y::Integer)
    if y<0; throw(DomainError(y, "`y` cannot be negative.")); end
    @noinline throw1(y) =
        throw(OverflowError("exponent $y is too large and computation will overflow"))
    if x== 1; return x; end
    if x==-1; return isodd(y) ? x : -x; end
    if y>typemax(Culong)
       x==0 && return x

       #At this point, x is not 1, 0 or -1 and it is not possible to use
       #gmpz_pow_ui to compute the answer. Note that the magnitude of the
       #answer is:
       #- at least 2^(2^32-1) ≈ 10^(1.3e9) (if Culong === UInt32).
       #- at least 2^(2^64-1) ≈ 10^(5.5e18) (if Culong === UInt64).
       #
       #Assume that the answer will definitely overflow.

       throw1(y)
    end
    return x^convert(Culong, y)
end

^(x::XInt   , y::XInt   ) = xint_pow(x, y)
^(x::XInt   , y::BigInt ) = xint_pow(x, y)
^(x::BigInt , y::XInt   ) = Base.GMP.bigint_pow(x, y)
^(x::XInt   , y::Bool   ) = y ? x : one(x)
^(x::XInt   , y::Integer) = xint_pow(x, y)
^(x::Integer, y::XInt   ) = xint_pow(XInt(x), y)
^(x::Bool   , y::XInt   ) = Base.power_by_squaring(x, y)

powermod(x::XInt, p::XInt, m::XInt) =
    is_short(x, p) && is_short(m) ? XInt(powermod(x.x, p.x, m.x)) :
                                    @bigint () x p m XInt(powermod(x, p, m))

powermod(x::Integer, p::Integer, m::XInt) = powermod(XInt(x), XInt(p), m)

function gcdx(a::XInt, b::XInt)
    is_short(a, b) &&
        return XInt.(gcdx(widen(a.x), widen(b.x)))
        # widen above necessary at least for (a, b) == (typemin(Int), 0)
    if iszero(b) # shortcut this to ensure consistent results with gcdx(a,b)
        return a < 0 ? (-a,-one(XInt),b) : (a,one(XInt),b)
    end
    g, s, t = XInt.(@bigint (x, y, z) a b MPZ.gcdext!(x, y, z, a, b))
    if t == 0
        # work around a difference in some versions of GMP
        if a == b
            return g, t, s
        elseif abs(a)==abs(b)
            return g, t, -s
        end
    end
    g, s, t
end

flipsign(x::XInt, y::Integer) = signbit(y) ? -x : x
flipsign(x::XInt, y::XInt)  = signbit(y) ? -x : x

factorial(x::XInt) = signbit(x) ? throw(DomainError(x)) :
                                  @bigint z x XInt(MPZ.fac_ui!(z, x))

binomial(x::XInt, k::UInt) = @bigint z x XInt(MPZ.bin_ui!(z, x, k))
binomial(x::XInt, k::Integer) = k < 0 ? XInt(0) : binomial(x, UInt(k))

function digits!(a::AbstractVector{T}, n::XInt; base::Integer = 10) where {T<:Integer}
    if 2 ≤ base ≤ 62
        s = codeunits(string(n; base=base))
        i, j = firstindex(a)-1, length(s)+1
        lasti = min(lastindex(a), firstindex(a) + length(s)-1 - isneg(n))
        while i < lasti
            # base ≤ 36: 0-9, plus a-z for 10-35
            # base > 36: 0-9, plus A-Z for 10-35 and a-z for 36..61
            x = s[j -= 1]
            a[i += 1] = base ≤ 36 ? (x>0x39 ? x-0x57 : x-0x30) : (x>0x39 ? (x>0x60 ? x-0x3d : x-0x37) : x-0x30)
        end
        lasti = lastindex(a)
        while i < lasti
            a[i+=1] = zero(T)
        end
        isneg(n) ? map!(-,a,a) : a
    else
        invoke(digits!, Tuple{typeof(a), Integer}, a, n; base=base) # slow generic fallback
    end
end

ndigits0zpb(x::XInt, b::Integer) = is_short(x) ? ndigits0zpb(x.x, b) :
                                                 @bigint () x ndigits0zpb(x, b)

function _prevpow2(x::XInt)
    if is_short(x)
        XInt(_prevpow2(x.x))
    else
        len = abs(x.x)
        @assert length(x.v) >= len
        high = @inbounds x.v[len]
        @assert !iszero(high) # like for BigInt/mpz_t
        shift = BITS_PER_LIMB - leading_zeros(high) - 1
        v = fill(zero(Limb), len)
        @inbounds v[len] = one(Limb) << shift
        _XInt(x.x, v)
    end
end

function _nextpow2(x::XInt)
    if is_short(x)
        y = _nextpow2(x.x)
        if sign(x) != sign(y)
            @assert y == typemin(Int)
            _XInt(1, Limb[typemin(Int) % Limb])
        else
            XInt(y)
        end
    else
        len = abs(x.x)
        @assert length(x.v) >= len
        popcount = GC.@preserve x MPZ.mpn_popcount(pointer(x.v), len)
        if popcount <= 1
            @assert popcount == 1
            return x
        end
        high = @inbounds x.v[len]
        @assert !iszero(high) # like for BigInt/mpz_t
        shift = (BITS_PER_LIMB - leading_zeros(high)) & (BITS_PER_LIMB-1)
        newlen = len + iszero(shift)
        v = fill(zero(Limb), newlen)
        @inbounds v[newlen] = one(Limb) << shift
        _XInt(flipsign(newlen, x.x), v)
    end
end

Base.checked_abs(x::XInt) = abs(x)
Base.checked_neg(x::XInt) = -x
Base.checked_add(a::XInt, b::XInt) = a + b
Base.checked_sub(a::XInt, b::XInt) = a - b
Base.checked_mul(a::XInt, b::XInt) = a * b
Base.checked_div(a::XInt, b::XInt) = div(a, b)
Base.checked_rem(a::XInt, b::XInt) = rem(a, b)
Base.checked_fld(a::XInt, b::XInt) = fld(a, b)
Base.checked_mod(a::XInt, b::XInt) = mod(a, b)
Base.checked_cld(a::XInt, b::XInt) = cld(a, b)
Base.add_with_overflow(a::XInt, b::XInt) = a + b, false
Base.sub_with_overflow(a::XInt, b::XInt) = a - b, false
Base.mul_with_overflow(a::XInt, b::XInt) = a * b, false

_copy(x::XInt) = _XInt(x.x, _copy(x.v))
_copy(v::Vector) = copy(v)
_copy(::Nothing) = nothing

Base.deepcopy_internal(x::XInt, d::IdDict) = get!(() -> _copy(x), d, x)

if Limb === UInt
    using Base: hash_uint

    function hash_integer(x::XInt, h::UInt)
        is_short(x) && return hash_integer(x.x, h)
        sz = x.x
        xv = x.v
        b = @inbounds xv[1]
        h ⊻= hash_uint(ifelse(sz < 0, -b, b) ⊻ h)
        for k = 2:abs(sz)
            h ⊻= hash_uint(@inbounds xv[k] ⊻ h)
        end
        h
    end

    function hash(x::XInt, h::UInt)
        GC.@preserve x begin
            is_short(x) && return hash(x.x, h)
            sz = x.x
            xv = x.v
            if sz == 1
                return hash(@inbounds(xv[1]), h)
            elseif sz == -1
                limb = @inbounds xv[1]
                limb <= typemin(Int) % UInt && return hash(-(limb % Int), h)
            end
            pow = trailing_zeros(x)
            nd = Base.ndigits0z(x, 2)
            idx = _divLimb(pow) + 1
            shift = _modLimb(pow) % UInt
            upshift = BITS_PER_LIMB - shift
            asz = abs(sz)
            if shift == 0
                limb = @inbounds xv[idx]
            else
                limb1 = @inbounds xv[idx]
                limb2 = idx < asz ? @inbounds(xv[idx+1]) : UInt(0)
                limb = limb2 << upshift | limb1 >> shift
            end
            if nd <= 1024 && nd - pow <= 53
                # TODO: add test for this branch
                return hash(ldexp(flipsign(Float64(limb), sz), pow), h)
            end
            h = hash_integer(1, h)
            h = hash_integer(pow, h)
            h ⊻= hash_uint(flipsign(limb, sz) ⊻ h)
            for idx = idx+1:asz
                if shift == 0
                    limb = @inbounds xv[idx]
                else
                    limb1 = limb2
                    if idx == asz
                        limb = limb1 >> shift
                        limb == 0 && break # don't hash leading zeros
                    else
                        limb2 = @inbounds xv[idx+1]
                        limb = limb2 << upshift | limb1 >> shift
                    end
                end
                h ⊻= hash_uint(limb ⊻ h)
            end
            return h
        end
    end
end

## Random

import Random: Sampler, rand, rand!
using Random: Random, Repetition, AbstractRNG

abstract type SamplerXInt <: Sampler{XInt} end

struct SamplerXIntBig{SP<:Sampler{Limb}} <: SamplerXInt
    a::XInt           # first
    m::XInt           # range length - 1 (generate z ∈ [0, m])
    nlimbsmax::SLimb  # max number of limbs for z+a
    highsp::SP        # sampler for the highest limb
end

safe_len(x::XInt) = is_short(x) ? SLimb(!iszero(x)) : abs(x.x)

@noinline function SamplerXIntBig(::Type{RNG}, r1::XInt, r2::XInt) where {RNG<:AbstractRNG}
    m = r2-r1
    isneg(m) && throw(ArgumentError("range must be non-empty"))
    if is_short(m)
        mvl = m.x % Limb
        nlimbsmax = 0
    else
        nlimbs, mv = lenvec(m)
        if nlimbs == 1
            mvl = @inbounds mv[1]
            nlimbsmax = 0
        else
            nlimbsmax = max(nlimbs, safe_len(r2), safe_len(r1))
            nlimbsmax += ispos(r1) # add! will request one more limb just in case
            mvl = @inbounds mv[nlimbs]
        end
    end
    highsp = Sampler(RNG, zero(Limb):mvl)
    SamplerXIntBig(r1, m, nlimbsmax, highsp)
end

struct SamplerXIntShort{SP<:Sampler{SLimb}} <: SamplerXInt
    sp::SP
end

@inline function Sampler(::Type{RNG}, r::AbstractUnitRange{XInt}, N::Repetition
                 ) where {RNG<:AbstractRNG}
    r1, r2 = first(r), last(r)
    if is_short(r1, r2)
        SamplerXIntShort(Sampler(RNG, r1.x:r2.x, N))
    else
        SamplerXIntBig(RNG, r1, r2)
    end
end

rand(rng::AbstractRNG, sp::SamplerXInt) = rand!(rng, _XInt(0), sp)

function rand!(rng::AbstractRNG, x::XInt, sp::SamplerXIntBig)
    iszero(sp.nlimbsmax) &&
        return add!(x, sp.a, _XInt(rand(rng, sp.highsp)))
    nlimbs, mv = lenvec(sp.m)
    limbs = vec!(x, sp.nlimbsmax)
    mvl = @inbounds mv[nlimbs]
    while true
        @inbounds rand!(rng, @view(limbs[1:nlimbs-1]))
        @inbounds h = limbs[nlimbs] = rand(rng, sp.highsp)
        h < mvl && break
        MPN.cmp(limbs, sp.m, nlimbs) <= 0 && break
    end
    len = normalize(limbs, nlimbs)
    add!(_XInt(len, limbs, false), _XInt(len, limbs, true), sp.a)
end

rand!(rng::AbstractRNG, ::XInt, sp::SamplerXIntShort) = _XInt(rand(rng, sp.sp))

end # module


text:
using AdventOfCode2021
using Printf
using Test

@testset "AdventOfCode2021" begin
    for day in CURRENT_DAY:-1:1 # get current results first
        testset_name = @sprintf("Day %02d", day)
        @testset "$testset_name" begin
            D = Day(day)
            @test test1(D) == testresult1(D)
            @printf("Solution of Day %2d – Part 1: %d\n", day, solve1(D))
            @test test2(D) == testresult2(D)
            @printf("Solution of Day %2d – Part 2: %d\n", day, solve2(D))
        end
    end
    # println("Solution of Day $(CURRENT_DAY) – Part 1: $(solve1(CURRENT_DAY))")
    # println("Solution of Day $(CURRENT_DAY) – Part 2: $(solve2(CURRENT_DAY))")
end

ext: jl


text:
using PythontexProfiler
using Test

@testset "PythontexProfiler.jl" begin

profile("test.pytxcode")

end

max_stars_repo_name: sisl/PythontexProfiler


text:
#!/usr/bin/env julia

using ArgParse
using BioStructures

function parse_commandline()
    settings = ArgParseSettings(description="""
    This script reads the `hhr` file returned by `hhsearch` using the pdb70
    database and download the needed PDB files into the indicated folder.
    """,)

    @add_arg_table settings begin
        "--input", "-i"
            help = "path to the `hhr` file returned by `hhsearch`."
            required = true
        "--pdb", "-p"
            help = """
            path where the pdb files in mmCIF format are going to be downloaded.
            """
            default = "."
    end

    return parse_args(settings)
end

"""
Return the list of pdb codes in the `hhr` file returned by `hhsearch`.

For example, a `hhr` file containing the following list in the header:

```

 No Hit                             Prob E-value P-value  Score    SS Cols Query HMM  Template HMM
  1 2FE0_A small myristoylated pro  19.8      73  0.0012   23.4   0.0   16    1-16      6-21  (136)
  2 4MCB_B tRNA (guanine-N(1)-)-me  10.2   2E+02  0.0031   23.8   0.0   10   60-69     52-61  (246)

```

Is going to return:

```
10-element Array{String,1}:
 "2FE0"
 "4MCB"
```
"""
function get_pdb_list(hhr_file)
    open(hhr_file, "r") do file
        parse = false
        pdb_list = String[]
        for line in eachline(file)
            if occursin("No Hit", line)
                parse = true
                continue
            end
            if parse
                beginning = match(r"^\s*[0-9]+\s(\S+)_\S+\s", line)
                if beginning === nothing
                    return pdb_list
                end
                push!(pdb_list, beginning.captures[1])
            end
        end
        pdb_list
    end
end

"""
Download each pdb in pdb_list into the output_path in mmCIF format.
It will retry each pdb until 5 times if errors occur.
If the filesystem is case-sensitive, a lowercase symbolic link to
the uppercase file is created.
"""
function download_pdbs(pdb_list, output_path)
    for pdb in pdb_list
        path_upper = retry(downloadpdb, delays=ExponentialBackOff(n=5))(
            pdb, file_format=MMCIF, pdb_dir=output_path
        )
        path_lower = joinpath(
            (endswith(x, ".cif") ? lowercase(x) : x for x in splitpath(path_upper))...)
        if !isfile(path_lower)
            symlink(path_upper, path_lower)
        end
    end
end

function main()
    parsed_args = parse_commandline()

    pdb_list = get_pdb_list(parsed_args["input"])

    output_path = abspath(parsed_args["pdb"])
    mkpath(output_path)

    download_pdbs(pdb_list, output_path)
end

main()

max_stars_repo_name: elolaine/PhyloSofS

max_issues_repo_name: elolaine/Phylosofs

max_forks_repo_name: PhyloSofS-Team/PhyloSofS


text:
const Initialized = Union{Type{Val{true}}, Type{Val{false}}}

const Pointers = Union{Ptr{Cvoid}, IndirectPointer}

struct OnDiskRepresentation{Offsets,JLTypes,H5Types} end
odr_sizeof(::Nothing) = 0
@Base.pure odr_sizeof(x::DataType) = Int(x.size)

struct UnknownType{T}
    name::T
    parameters::Vector{Any}

    UnknownType{T}(name) where T = new(name)
    UnknownType{T}(name, parameters) where T = new(name, parameters)
end
UnknownType(name) = UnknownType{typeof(name)}(name)
UnknownType(name, parameters) = UnknownType{typeof(name)}(name, parameters)

struct Vlen{T}
    size::UInt32
    id::GlobalHeapID
end
odr_sizeof(::Type{T}) where {T<:Vlen} = 4 + jlsizeof(GlobalHeapID)

# Look up the corresponding committed datatype for a given type
macro lookup_committed(f, T)
    quote
        cdt = get($(esc(f)).jlh5type, $(esc(T)), nothing)
        cdt !== nothing && return cdt::CommittedDatatype
    end
end


text:
module KernelsTests

using Test
using UnstructuredGrids.Kernels

vv = [[1,2,3],[2,3],[5,8],Int[],[1,2,4]]

_data, _ptrs = generate_data_and_ptrs(vv)

data = [1, 2, 3, 2, 3, 5, 8, 1, 2, 4]
ptrs = [1, 4, 6, 8, 8, 11]

@test data == _data
@test ptrs == _ptrs

a = [9,2,1,2,4,7,4]
b = [1,9,2,1,2,4,7]

rewind_ptrs!(a)
@test a == b

a = [3,2,4,2]
b = [1,3,7,9]

length_to_ptrs!(a)
@test a == b

include("Mock2D.jl")

_vertex_to_cells_data, _vertex_to_cells_ptrs = generate_dual_connections(
  cell_to_vertices_data,cell_to_vertices_ptrs)

@test _vertex_to_cells_data == vertex_to_cells_data
@test _vertex_to_cells_ptrs == vertex_to_cells_ptrs

_cell_to_faces_data, _cell_to_faces_ptrs = generate_cell_to_faces(
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_cells_data,
  vertex_to_cells_ptrs)

@test cell_to_faces_data == _cell_to_faces_data
@test cell_to_faces_ptrs == _cell_to_faces_ptrs

_face_to_ftype = generate_face_to_ftype(
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_ftype)

@test face_to_ftype == _face_to_ftype

_face_to_vertices_data, _face_to_vertices_ptrs = generate_face_to_vertices(
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  cell_to_ctype,
  ctype_to_lface_to_lvertices)

@test face_to_vertices_data == _face_to_vertices_data
@test face_to_vertices_ptrs == _face_to_vertices_ptrs

pa = [1,3,5,7,9]
pb = [1,3,5,7]

pc = append_ptrs(pa,pb)

@test pc == [1, 3, 5, 7, 9, 11, 13, 15]

vertex_to_faces_data, vertex_to_faces_ptrs = generate_dual_connections(
  face_to_vertices_data,face_to_vertices_ptrs)

_cell_to_faces_data, _cell_to_faces_ptrs = find_cell_to_faces(
  cell_to_vertices_data,
  cell_to_vertices_ptrs,
  ctype_to_lface_to_lvertices,
  cell_to_ctype,
  vertex_to_faces_data,
  vertex_to_faces_ptrs)

@test cell_to_faces_data == _cell_to_faces_data
@test cell_to_faces_ptrs == _cell_to_faces_ptrs

face_to_cells_data, face_to_cells_ptrs = generate_dual_connections(
  cell_to_faces_data, cell_to_faces_ptrs)

_face_to_isboundary = generate_facet_to_isboundary(face_to_cells_ptrs)

@test face_to_isboundary == _face_to_isboundary

_vertex_to_isboundary = generate_face_to_isboundary(
  face_to_isboundary,
  vertex_to_faces_data,
  vertex_to_faces_ptrs)

@test vertex_to_isboundary == _vertex_to_isboundary

gface_to_vertices = [[1,2],[2,5],[4,5]]

gface_to_vertices_data, gface_to_vertices_ptrs =
  generate_data_and_ptrs(gface_to_vertices)

gface_to_face = find_gface_to_face(
  face_to_vertices_data,
  face_to_vertices_ptrs,
  vertex_to_faces_data,
  vertex_to_faces_ptrs,
  gface_to_vertices_data,
  gface_to_vertices_ptrs)

@test gface_to_face == [1,4,2]

end # module KernelsTests


text:
using DataFrames

function generate_dataframes(config, exp1, exp2, exptest; generate_uniform = false)
    max_trajectories = 100
    Ntrajectories = min(config["trainsamples"], max_trajectories)
    Ntrajectories_test = min(config["testsamples"], max_trajectories)
    Ntrajectorysamples = Int(ceil(config["trainsamples"]/max_trajectories))
    scaling = Int(0.01/config["Δtsim"])
    samplerange = 1:(2*Int(1/config["Δtsim"]) - scaling)
    traindf = DataFrame(sold = Vector{Vector{State}}(), scurr = Vector{Vector{State}}())
    generate_uniform && (traindf_uniform = DataFrame(sold = Vector{Vector{State}}(), scurr = Vector{Vector{State}}()))
    threadlock = ReentrantLock()  # Push to df not atomic
    Threads.@threads for _ in 1:div(Ntrajectories, 2)
        storage = _run_experiment!(exp1)
        _pushsamples!(storage, traindf, Ntrajectorysamples, samplerange, [0, scaling], scaling, threadlock)
        generate_uniform && _pushuniformsamples!(storage, traindf_uniform, Ntrajectorysamples, samplerange, [0, scaling], threadlock)
    end
    Threads.@threads for _ in 1:div(Ntrajectories, 2)
        storage = _run_experiment!(exp2)
        _pushsamples!(storage, traindf, Ntrajectorysamples, samplerange, [0, scaling], scaling, threadlock)
        generate_uniform && _pushuniformsamples!(storage, traindf_uniform, Ntrajectorysamples, samplerange, [0, scaling], threadlock)
    end
    testdf = DataFrame(sold = Vector{Vector{State}}(), sfuture = Vector{Vector{State}}())
    generate_uniform && (testdf_uniform = DataFrame(sold = Vector{Vector{State}}(), sfuture = Vector{Vector{State}}()))
    samplerange = 1:(2*Int(1/config["Δtsim"]) - scaling*(config["simsteps"]+1))
    Threads.@threads for _ in 1:Ntrajectories_test
        storage = _run_experiment!(exptest)
        _pushsamples!(storage, testdf, Int(ceil(config["testsamples"]/Ntrajectories_test)), samplerange, [0, scaling*(config["simsteps"]+1)], scaling, threadlock)
        generate_uniform && _pushuniformsamples!(storage, testdf_uniform, Int(ceil(config["testsamples"]/Ntrajectories_test)), samplerange, [0, scaling*(config["simsteps"]+1)], threadlock)
    end
    generate_uniform && return traindf, testdf, traindf_uniform, testdf_uniform
    return traindf, testdf
end

function _run_experiment!(experiment; maxruns = 10)
    for _ in 1:maxruns  # Retry experiment if simulation fails
        try
            return experiment()
        catch e
            @warn "Experiment failed, retrying..."
            display(e)
            continue
        end
    end
    throw(ErrorException("Experiment failed to execute $maxruns times"))
end

function _pushsamples!(storage, df, nsamples, samplerange, indexoffset, scaling, threadlock)
    @assert (4scaling+1)*nsamples < length(samplerange)
    indiceset = Set()
    for _ in 1:nsamples
        j = 0
        while true
            j = rand(samplerange)  # End of storage - required steps
            !any([j in ind-2scaling:ind+2scaling for ind in indiceset]) && break  # Sample j outside of existing indices
        end
        lock(threadlock)
        try
            push!(df, [getStates(storage, j+offset) for offset in indexoffset])
            push!(indiceset, j)
        finally
            unlock(threadlock)
        end
    end
end

function _pushuniformsamples!(storage, df, nsamples, samplerange, indexoffset, threadlock)
    @assert length(samplerange) > nsamples
    indiceset = Set()
    indices = nsamples > 1 ? [samplerange[i*div(length(samplerange),nsamples)] for i in 1:nsamples] : div(length(samplerange),2)
    for j in indices
        lock(threadlock)
        try
            push!(df, [getStates(storage, j+offset) for offset in indexoffset])
            push!(indiceset, j)
        finally
            unlock(threadlock)
        end
    end
end


text:
# THESE 4 CANNOT BE CHANGED BY GMT
const global GMT_SESSION_NORMAL   = 0   # Typical mode to GMT_Create_Session
const global GMT_SESSION_NOEXIT   = 1   # Call return and not exit when error
const global GMT_SESSION_EXTERNAL = 2   # Called by an external API
const global GMT_SESSION_COLMAJOR = 4   # External API uses column-major formats. [Row-major format]
const global GMT_SESSION_NOGDALCLOSE = 64   # External API tells GMT to not call GDALDestroyDriverManager()

const global GMT_SESSION_BITFLAGS = GMT_SESSION_NOEXIT + GMT_SESSION_EXTERNAL + GMT_SESSION_NOGDALCLOSE + GMT_SESSION_COLMAJOR
API = GMT_Create_Session("GMT", 2, GMT_SESSION_BITFLAGS)

enu = GMT_Get_Enum(API, "GMT_CHAR");	const global GMT_CHAR  = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_UCHAR");	const global GMT_UCHAR  = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_SHORT");	const global GMT_SHORT  = (enu != -99999) ? enu : 2
enu = GMT_Get_Enum(API, "GMT_USHORT");	const global GMT_USHORT  = (enu != -99999) ? enu : 3
enu = GMT_Get_Enum(API, "GMT_INT");		const global GMT_INT  = (enu != -99999) ? enu : 4
enu = GMT_Get_Enum(API, "GMT_UINT");	const global GMT_UINT  = (enu != -99999) ? enu : 5
enu = GMT_Get_Enum(API, "GMT_LONG");	const global GMT_LONG  = (enu != -99999) ? enu : 6
enu = GMT_Get_Enum(API, "GMT_ULONG");	const global GMT_ULONG  = (enu != -99999) ? enu : 7
enu = GMT_Get_Enum(API, "GMT_FLOAT");	const global GMT_FLOAT  = (enu != -99999) ? enu : 8
enu = GMT_Get_Enum(API, "GMT_DOUBLE");	const global GMT_DOUBLE  = (enu != -99999) ? enu : 9
enu = GMT_Get_Enum(API, "GMT_TEXT");	const global GMT_TEXT  = (enu != -99999) ? enu : 10

enu = GMT_Get_Enum(API, "GMT_OPT_INFILE");	const global GMT_OPT_INFILE = (enu != -99999) ? enu : 60
enu = GMT_Get_Enum(API, "GMT_IS_OUTPUT");	const global GMT_IS_OUTPUT  = (enu != -99999) ? enu : 1024
enu = GMT_Get_Enum(API, "GMT_VIA_MATRIX");	const global GMT_VIA_MATRIX = (enu != -99999) ? enu : 256

# GMT_enum_container
enu = GMT_Get_Enum(API, "GMT_CONTAINER_AND_DATA");	const global GMT_CONTAINER_AND_DATA  = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_CONTAINER_ONLY");	const global GMT_CONTAINER_ONLY  = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_DATA_ONLY");		const global GMT_DATA_ONLY  = (enu != -99999) ? enu : 2
enu = GMT_Get_Enum(API, "GMT_WITH_STRINGS");	const global GMT_WITH_STRINGS  = (enu != -99999) ? enu : 32
enu = GMT_Get_Enum(API, "GMT_NO_STRINGS");		const global GMT_NO_STRINGS  = (enu != -99999) ? enu : 0

# GMT_enum_read
enu = GMT_Get_Enum(API, "GMT_READ_DATA");	const global GMT_READ_DATA  = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_READ_TEXT");	const global GMT_READ_TEXT  = (enu != -99999) ? enu : 2
enu = GMT_Get_Enum(API, "GMT_READ_MIXED");	const global GMT_READ_MIXED = (enu != -99999) ? enu : 3

# begin enum GMT_enum_family
enu = GMT_Get_Enum(API, "GMT_IS_DATASET");	const global GMT_IS_DATASET = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_IS_GRID");		const global GMT_IS_GRID  = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_IS_IMAGE");	const global GMT_IS_IMAGE = (enu != -99999) ? enu : 2
enu = GMT_Get_Enum(API, "GMT_IS_PALETTE");		const global GMT_IS_PALETTE = (enu != -99999) ? enu : 3
enu = GMT_Get_Enum(API, "GMT_IS_POSTSCRIPT");	const global GMT_IS_POSTSCRIPT = (enu != -99999) ? enu : 4
enu = GMT_Get_Enum(API, "GMT_IS_MATRIX");		const global GMT_IS_MATRIX = (enu != -99999) ? enu : 5
enu = GMT_Get_Enum(API, "GMT_IS_VECTOR");		const global GMT_IS_VECTOR = (enu != -99999) ? enu : 6
enu = GMT_Get_Enum(API, "GMT_IS_CUBE");		const global GMT_IS_CUBE  = (enu != -99999) ? enu : 7

enu = GMT_Get_Enum(API, "GMT_COMMENT_IS_TEXT");		const global GMT_COMMENT_IS_TEXT = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_IMAGE_ALPHA_LAYER");	const global GMT_IMAGE_ALPHA_LAYER = (enu != -99999) ? enu : 8192

# begin enum GMT_module_enum
enu = GMT_Get_Enum(API, "GMT_MODULE_EXIST");	const global GMT_MODULE_EXIST = (enu != -99999) ? enu : -3
enu = GMT_Get_Enum(API, "GMT_MODULE_OPT");	const global GMT_MODULE_OPT = (enu != -99999) ? enu : -1
enu = GMT_Get_Enum(API, "GMT_MODULE_CMD");	const global GMT_MODULE_CMD = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_SYNOPSIS");	const global GMT_SYNOPSIS = (enu != -99999) ? enu : 1

# begin enum GMT_io_enum
enu = GMT_Get_Enum(API, "GMT_IN");	const global GMT_IN = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_OUT");	const global GMT_OUT = (enu != -99999) ? enu : 1

enu = GMT_Get_Enum(API, "GMT_ALLOC_EXTERNALLY");	const global GMT_ALLOC_EXTERNALLY = (enu != -99999) ? enu : 0
# begin enum GMT_enum_dimindex
enu = GMT_Get_Enum(API, "GMT_SEG");	const global GMT_SEG = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_ROW");	const global GMT_ROW = (enu != -99999) ? enu : 2
enu = GMT_Get_Enum(API, "GMT_COL");	const global GMT_COL = (enu != -99999) ? enu : 3

#const global GMT_GRID_ALL = 0		# This one is a #define

# begin enum GMT_enum_fmt
enu = GMT_Get_Enum(API, "GMT_IS_COL_FORMAT");	const global GMT_IS_COL_FORMAT = (enu != -99999) ? enu : 2

# begin enum GMT_enum_geometry
enu = GMT_Get_Enum(API, "GMT_IS_PLP");	const global GMT_IS_PLP = (enu != -99999) ? enu : 7
enu = GMT_Get_Enum(API, "GMT_IS_SURFACE");	const global GMT_IS_SURFACE = (enu != -99999) ? enu : 8
enu = GMT_Get_Enum(API, "GMT_IS_VOLUME");	const global GMT_IS_VOLUME  = (enu != -99999) ? enu : 9
enu = GMT_Get_Enum(API, "GMT_IS_NONE");	const global GMT_IS_NONE = (enu != -99999) ? enu : 16

# begin enum GMT_enum_color
enu = GMT_Get_Enum(API, "GMT_RGB");		const global GMT_RGB  = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_CMYK");	const global GMT_CMYK = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_HSV");		const global GMT_HSV = (enu != -99999) ? enu : 2

# GMT_enum_cptflags
enu = GMT_Get_Enum(API, "GMT_CPT_HINGED");	const global GMT_CPT_HINGED = (enu != -99999) ? enu : 4

# GMT_enum_CPT
enu = GMT_Get_Enum(API, "GMT_IS_PALETTE_KEY");		const global GMT_IS_PALETTE_KEY = (enu != -99999) ? enu : 1024
enu = GMT_Get_Enum(API, "GMT_IS_PALETTE_LABEL");	const global GMT_IS_PALETTE_LABEL = (enu != -99999) ? enu : 2048

# GMT_enum_workflowmode
enu = GMT_Get_Enum(API, "GMT_USE_WORKFLOW");		const global GMT_USE_WORKFLOW = (enu != -99999) ? enu : 0
enu = GMT_Get_Enum(API, "GMT_BEGIN_WORKFLOW");		const global GMT_BEGIN_WORKFLOW = (enu != -99999) ? enu : 1
enu = GMT_Get_Enum(API, "GMT_END_WORKFLOW");		const global GMT_END_WORKFLOW = (enu != -99999) ? enu : 2

GMT_Destroy_Session(API);


text:
# This file was generated, do not modify it. # hide
@load KNNClassifier pkg="NearestNeighbors"
@load MultinomialClassifier pkg="MLJLinearModels";

@pipeline KnnPipe(std=Standardizer(), clf=KNNClassifier()) is_probabilistic=true
@pipeline MnPipe(std=Standardizer(), clf=MultinomialClassifier()) is_probabilistic=true


text:
export
    ecl_to_equ,
    equ_to_ecl,
    equ_to_horiz,
    equ_to_horiz1,
    geocentric_to_topocentric

"""
    ecl_to_equ(long, lat, obliquity)

Convert ecliptic to equatorial coordinates.

[Meeus-1998: equations 13.3, 13.4]

Parameters:
    long : ecliptic longitude in radians
    lat : ecliptic latitude in radians
    obliquity : obliquity of the ecliptic in radians

Returns:
    Right ascension in radians (2 pi is 24 hours)
    Declination in radians
"""
function ecl_to_equ(long, lat, obliquity)
    cose = cos(obliquity)
    sine = sin(obliquity)
    sinl = sin(long)
    x = cos(long)
    y = sinl * cose - tan(lat) * sine
    alpha1 = atan(y, x)
    ra = mod2pi(alpha1)
    # convert to hours
    declin = asin(sin(lat) * cose + cos(lat) * sine * sinl)
    return (ra, declin)
end

"""
    equ_to_horiz(H, decl)

Convert equatorial to horizontal coordinates.

[Meeus-1998: equations 13.5, 13.6]

Note that azimuth is measured westward starting from the south.

This is not a good formula for using near the poles.

Parameters:
    H : hour angle in radians
    decl : declination in radians

Returns:
    azimuth in radians
    altitude in radians
"""
function equ_to_horiz(H, decl)
    global latitude
    cosH = cos(H)
    sinLat = sin(latitude)
    cosLat = cos(latitude)
    alt = asin((sinLat * sin(decl)) + (cosLat * cos(decl) * cosH))

    # atan(y,x)... not the other way round
    azi = atan((-1 * cos(decl)) * cosLat * sin(H), sin(decl) - (sinLat * sin(alt)))

    # Duffett-Smith says if azi positive, true azimuth is (360 - azi)
    if sin(H) > 0
    azi = mod2pi(azi)
    end
    return (azi, alt)
end

"""
    equ_to_ecl(ra, declin, obliquity)

Convert equatorial to ecliptic coordinates.

[Meeus-1998: equations 13.1, 13.2]

Parameters:
    ra : right accension radians
    declin : declination in radians
    obliquity : obliquity of the ecliptic in radians

Returns:
    ecliptic longitude in radians
    ecliptic latitude in radians
"""
function equ_to_ecl(ra, declin, obliquity)
    cose = cos(obliquity)
    sine = sin(obliquity)
    sina = sin(ra)
    long = mod2pi(atan(sina * cose + tan(declin) * sine, cos(ra)))
    lat = mod2pi(asin(sin(declin) * cose - cos(declin) * sine * sina))
    return (long, lat)
end

"""
    geocentric_to_topocentric(phi, H, L, ra, decl, d, jd)

Convert geocentric to topocentric coordinates.

[Meeus-1998: equations 40.2, 40.3]

Parameters:
    phi: geocentric radius in radians
    H : observer's altitude (above sea level) in meters
    L ; observer longitude (in radians)
    ra : body right ascension
    decl : body declination in radians
    d : body distance (in AU)
    jd : Julian date

Returns:
    topocentric right ascension in radians
    topocentric declination in radians
"""
function geocentric_to_topocentric(phi, H, L, ra, decl, d, jd)
     sinparallax = sin(deg2rad(8.794/3600))/d
     ha    = mod2pi(apparent_sidereal_time_greenwich(jd) - L - deg2rad(ra* 15))
     p1, psinp1, pcosp1 = geographical_to_geocentric_lat(phi, H)
     tmp   = cos(decl) - pcosp1 * sinparallax * cos(ha)
     da    = atan(-pcosp1 * sinparallax * sin(ha), tmp)
     tdecl = atan((sin(decl) - psinp1 * sinparallax) * cos(da), tmp)
    return (ra+da, tdecl)
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

precompile(Tuple{typeof(.init_FMU),                         Ptr{Cchar} })

precompile(Tuple{typeof(.fmi2GetTypesPlatform)})
precompile(Tuple{typeof(.fmi2GetVersion})
precompile(Tuple{typeof(.fmi2Instantiate),                  fmi2String, fmi2Type, fmi2String, fmi2String, Ptr{fmi2CallbackFunctions}, fmi2Boolean, fmi2Boolean})
precompile(Tuple{typeof(.fmi2FreeInstance),                 fmi2Component})
precompile(Tuple{typeof(.fmi2SetDebugLogging),              fmi2Component, fmi2Boolean, Csize_t, Ptr{fmi2String}})
precompile(Tuple{typeof(.fmi2SetupExperiment),              fmi2Component, fmi2Boolean, fmi2Real, fmi2Real, fmi2Boolean, fmi2Real})
precompile(Tuple{typeof(.fmi2EnterInitializationMode),      fmi2Component})
precompile(Tuple{typeof(.fmi2ExitInitializationMode),       fmi2Component})
precompile(Tuple{typeof(.fmi2Terminate),                    fmi2Component})
precompile(Tuple{typeof(.fmi2Reset),                        fmi2Component})
precompile(Tuple{typeof(.fmi2GetReal),                      fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Real}})
precompile(Tuple{typeof(.fmi2GetInteger),                   fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Integer}})
precompile(Tuple{typeof(.fmi2GetBoolean),                   fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Boolean}})
precompile(Tuple{typeof(.fmi2GetString),                    fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2String}})
precompile(Tuple{typeof(.fmi2SetReal),                      fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Real}})
precompile(Tuple{typeof(.fmi2SetInteger),                   fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Integer}})
precompile(Tuple{typeof(.fmi2SetBoolean),                   fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2Boolean}})
precompile(Tuple{typeof(.fmi2SetString),                    fmi2Component, Ptr{fmi2ValueReference}, Csize_t, Ptr{fmi2String}})
precompile(Tuple{typeof(.fmi2SetTime),                      fmi2Component, fmi2Real})
precompile(Tuple{typeof(.fmi2SetContinuousStates),          fmi2Component, Ptr{fmi2Real}, Csize_t})
precompile(Tuple{typeof(.fmi2EnterEventMode),               fmi2Component})
precompile(Tuple{typeof(.fmi2NewDiscreteStates),            fmi2Component, Ptr{fmi2EventInfo}})
precompile(Tuple{typeof(.fmi2EnterContinuousTimeMode),      fmi2Component})
precompile(Tuple{typeof(.fmi2CompletedIntegratorStep),      fmi2Component, fmi2Boolean, Ptr{fmi2Boolean}, Ptr{fmi2Boolean}})
precompile(Tuple{typeof(.fmi2GetDerivatives),               fmi2Component, Ptr{fmi2Real}, Csize_t})
precompile(Tuple{typeof(.fmi2GetEventIndicators),           fmi2Component, Ptr{fmi2Real}, Csize_t})
precompile(Tuple{typeof(.fmi2GetContinuousStates),          fmi2Component, Ptr{fmi2Real}, Csize_t})
precompile(Tuple{typeof(.fmi2GetNominalsOfContinuousStates),fmi2Component, Ptr{fmi2Real}, Csize_t})


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

function node1(url, msg)
    sock = nng_socket(0)
    rv = nng_push0_open(sock)
    check_err("nng_push0_open", rv)
    dp = nng_dialer(0)
    rv = nng_dial(sock, url, dp, 0)
    check_err("nng_dial", rv)
    println("NODE1: SENDING ", msg)
    rv = nng_send(sock, msg, 0)
    sleep(1)
    nng_close(sock)
end

node1(url, "HelloWorld")


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

mutable struct ApplicationGatewayRewriteRuleCondition <: SwaggerModel
    variable::Any # spec type: Union{ Nothing, String } # spec name: variable
    pattern::Any # spec type: Union{ Nothing, String } # spec name: pattern
    ignoreCase::Any # spec type: Union{ Nothing, Bool } # spec name: ignoreCase
    negate::Any # spec type: Union{ Nothing, Bool } # spec name: negate

    function ApplicationGatewayRewriteRuleCondition(;variable=nothing, pattern=nothing, ignoreCase=nothing, negate=nothing)
        o = new()
        validate_property(ApplicationGatewayRewriteRuleCondition, Symbol("variable"), variable)
        setfield!(o, Symbol("variable"), variable)
        validate_property(ApplicationGatewayRewriteRuleCondition, Symbol("pattern"), pattern)
        setfield!(o, Symbol("pattern"), pattern)
        validate_property(ApplicationGatewayRewriteRuleCondition, Symbol("ignoreCase"), ignoreCase)
        setfield!(o, Symbol("ignoreCase"), ignoreCase)
        validate_property(ApplicationGatewayRewriteRuleCondition, Symbol("negate"), negate)
        setfield!(o, Symbol("negate"), negate)
        o
    end
end # type ApplicationGatewayRewriteRuleCondition

const _property_map_ApplicationGatewayRewriteRuleCondition = Dict{Symbol,Symbol}(Symbol("variable")=>Symbol("variable"), Symbol("pattern")=>Symbol("pattern"), Symbol("ignoreCase")=>Symbol("ignoreCase"), Symbol("negate")=>Symbol("negate"))
const _property_types_ApplicationGatewayRewriteRuleCondition = Dict{Symbol,String}(Symbol("variable")=>"String", Symbol("pattern")=>"String", Symbol("ignoreCase")=>"Bool", Symbol("negate")=>"Bool")
Base.propertynames(::Type{ ApplicationGatewayRewriteRuleCondition }) = collect(keys(_property_map_ApplicationGatewayRewriteRuleCondition))
Swagger.property_type(::Type{ ApplicationGatewayRewriteRuleCondition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayRewriteRuleCondition[name]))}
Swagger.field_name(::Type{ ApplicationGatewayRewriteRuleCondition }, property_name::Symbol) =  _property_map_ApplicationGatewayRewriteRuleCondition[property_name]

function check_required(o::ApplicationGatewayRewriteRuleCondition)
    true
end

function validate_property(::Type{ ApplicationGatewayRewriteRuleCondition }, name::Symbol, val)
end


text:
using Test

@testset "FSimBase" begin
    include("minimal.jl")
    include("custom_example.jl")
end


text:
module IntervalWavelets

import OffsetArrays
import LinearAlgebra
import Wavelets

using RecipesBase

export
    DyadicRational,
    Filter,
    InteriorFilter,
    BoundaryFilter,
    Sides,
    LEFT,
    RIGHT,

    all_dyadic_rationals,
    boundary_filters,
    boundary_scaling_functions,
    coefficients,
    filters,
    interior_filter,
    interior_scaling_function,
    numerator,
    reduce,
    resolution,
    reconstruct,
    side,
    support,
    support_boundaries,
    vanishing_moments,
    interval_scaling_function_basis

const sqrt2 = sqrt(2)

include("DyadicRationals.jl")
include("Enums.jl")
include("Filters.jl")
include("Filters/Interior.jl")
include("Filters/Left.jl")
include("Filters/Right.jl")
include("InteriorScalingFunction.jl")
include("BoundaryScalingFunction.jl")
include("Basis.jl")
include("Misc.jl")

end # module


text:
"""
    MMAOptions

A struct that stores all the options of the MMA algorithms. Th following are the fields of `MMAOptions`:
 - `maxiter`: the maximum number of inner iterations. For `MMA87`, there is 1 inner iteration per outer iteration.
 - `outer_maxiter`: the maximum number of outer iterations.
 - `maxinner`: the maximum number of inner iterations per outer iteration of [`MMA02`](@ref). Not applicable for [`MMA87`](@ref).
 - `tol`: a tolerance struct of type [`Tolerance`](@ref).
 - `s_init`: defined in the original [`MMA02`](@ref) paper.
 - `s_incr`: defined in the original [`MMA02`](@ref) paper.
 - `s_decr`: defined in the original [`MMA02`](@ref) paper.
 - `store_trace`: if true, a trace will be stored.
 - `dual_options`: the options passed to the dual optimizer from [`Optim.jl`](https://github.com/JuliaNLSolvers/Optim.jl).
"""
@with_kw mutable struct MMAOptions{T, Ttol <: Tolerance, TSubOptions <: Optim.Options}
    maxiter::Int = 1000
    outer_maxiter::Int = 10^8
    maxinner::Int = 10
    tol::Ttol = Tolerance()
    s_init::T = 0.5
    s_incr::T = 1.2
    s_decr::T = 0.7
    store_trace::Bool = false
    show_trace::Bool = false
    auto_scale::Bool = false
    keep_best::Bool = false
    dual_options::TSubOptions = Optim.Options(
        allow_outer_f_increases = true,
        allow_f_increases = true,
        iterations = 1000,
        outer_iterations=1000,
    )
end


text:
using RoME

using Test

@testset "test tree based autoinit on small bimodal point2 trilateration..." begin

N=100
fg = initfg()

addVariable!(fg, :x0, Point2, N=N)
addFactor!(fg, [:x0], PriorPoint2(MvNormal([100.0;0], Matrix{Float64}(LinearAlgebra.I, 2,2))), autoinit=false)

addVariable!(fg, :x1, Point2, N=N)
addFactor!(fg, [:x1], PriorPoint2(MvNormal([0.0;100.0], Matrix{Float64}(LinearAlgebra.I, 2,2))), autoinit=false)

addVariable!(fg, :l1, Point2, N=N)
addFactor!(fg, [:x0;:l1], Point2Point2Range(Normal(100.0, 1.0)), autoinit=false)
addFactor!(fg, [:x1;:l1], Point2Point2Range(Normal(100.0, 1.0)), autoinit=false)

# tree = wipeBuildNewTree!(fg)
# drawTree(tree, filepath="/tmp/caesar/bt.pdf", show=true)

# eo = getEliminationOrder(fg, ordering=:qr)
# eo = [1;3;5]
# tree = buildTreeFromOrdering!(fg,eo)

tree, smt, hist = solveTree!(fg)

# cliq = tree.cliques[2]
# clst = cliqInitSolveUpByStateMachine!(fg, tree, cliq, drawtree=false, limititers=15 )

# cliq = tree.cliques[1]
# clst = cliqInitSolveUpByStateMachine!(fg, tree, cliq, drawtree=false, limititers=15 )

# using RoMEPlotting
#
# plotKDE(fg, [:x0;:x1]) |> SVG("/tmp/test.svg") || @async run(`eog /tmp/test.svg`)
# plotKDE(fg, :l1) |> SVG("/tmp/test.svg") || @async run(`eog /tmp/test.svg`)

@test 10 < sum( 80 .< getVal(fg, :l1)[1,:] .< 120 )
@test 10 < sum( -20 .< getVal(fg, :l1)[2,:] .< 20 )

@test 10 < sum( -20 .< getVal(fg, :l1)[1,:] .< 20 )
@test 10 < sum( 80 .< getVal(fg, :l1)[2,:] .< 120 )

@test 80 < sum( 80 .< getVal(fg, :x0)[1,:] .< 120 )
@test 80 < sum( -20 .< getVal(fg, :x0)[2,:] .< 20 )

@test 80 < sum( -20 .< getVal(fg, :x1)[1,:] .< 20 )
@test 80 < sum( 80 .< getVal(fg, :x1)[2,:] .< 120 )

# ett = ExploreTreeType(fg, tree, tree.cliques[1], nothing, NBPMessage[])
# downMsgPassingIterative!(ett,N=100, dbg=false, drawpdf=false);

@test 10 < sum( 80 .< getVal(fg, :l1)[1,:] .< 120 )
@test 10 < sum( -20 .< getVal(fg, :l1)[2,:] .< 20 )

@test 10 < sum( -20 .< getVal(fg, :l1)[1,:] .< 20 )
@test 10 < sum( 80 .< getVal(fg, :l1)[2,:] .< 120 )

@test 80 < sum( 80 .< getVal(fg, :x0)[1,:] .< 120 )
@test 80 < sum( -20 .< getVal(fg, :x0)[2,:] .< 20 )

@test 80 < sum( -20 .< getVal(fg, :x1)[1,:] .< 20 )
@test 80 < sum( 80 .< getVal(fg, :x1)[2,:] .< 120 )

end


text:
# Practical Machine learning
# Deep learning - Sparse autoencoder example
# Chapter 11

module SparseAutoencoder

function sigmoid(z)
        return one(z[1]) / (one(z[1]) + exp(-z))
end

function autoencodercost!(theta::Array{Float64,1}, visiblesize, hiddensize, data::Array{Float64,2}, grad::Vector{Float64};
                          lambda=0.0, sparsityparameter=0.01, beta=0.0)
        W1 = reshape(theta[1:hiddensize*visiblesize],hiddensize,visiblesize)
        W2 = reshape(theta[hiddensize*visiblesize+1:2*hiddensize*visiblesize],visiblesize,hiddensize)
        b1 = theta[2*hiddensize*visiblesize+1:2*hiddensize*visiblesize+hiddensize]
        b2 = theta[2*hiddensize*visiblesize+hiddensize+1:end]
        m = size(data,2)

        # Forward pass
        a1 = data
        z2 = broadcast(+,W1*a1,b1)
        a2 = sigmoid(z2)
        z3 = broadcast(+,W2*a2,b2)
        a3 = sigmoid(z3)
        output = a3

        avghiddenactivity = mean(a2,2)
        sparsitycost = sum(sparsityparameter*log(sparsityparameter/avghiddenactivity) + (1 - sparsityparameter)*log((1-sparsityparameter)/(1-avghiddenactivity)))
        cost = sum((output - data).^2)/(2*m) + lambda/2*sum(theta[1:2*hiddensize*visiblesize].^2) + beta*sparsitycost

        if length(grad) > 0
            # Back-propagation for sparsified objective
            a3prime = a3.*(one(a3[1])-a3)
            a2prime = a2.*(one(a2[1])-a2)
            sparsityterm = beta*(-sparsityparameter/avghiddenactivity + (1-sparsityparameter)/(1-avghiddenactivity))

            d3 = -(data - a3).* a3prime
            d2 = broadcast(+,W2'*d3,sparsityterm).*a2prime

            grad[1:visiblesize*hiddensize] = (d2*a1'/m + lambda*W1)[:] #dW1
            grad[visiblesize*hiddensize+1:2*visiblesize*hiddensize] = (d3*a2'/m + lambda*W2)[:] #dW2
            grad[2*visiblesize*hiddensize+1:2*hiddensize*visiblesize+hiddensize] = (sum(d2,2)/m)[:]
            grad[2*hiddensize*visiblesize+hiddensize+1:end] = (sum(d3,2)/m)[:]
        end
        return cost

end

function numericalgradient(f,x::Array{Float64,1})
        epsilon = 1e-4
        g = zeros(length(x))
        for i =1:length(x)
                x[i] += epsilon
                fplus = f(x)
                x[i] -= 2*epsilon
                fminus = f(x)
                g[i] = (fplus - fminus)/(2*epsilon)
                x[i] += epsilon
         end
         return g
end

function testGradient()
        vars = matread("IMAGES.mat")
        traindata = Array(Float64,SAMPLE_WIDTH*SAMPLE_HEIGHT,100)
        sampleimages!(vars["IMAGES"],traindata,100)

        hidden = 10
        visible = 64
        lambda = 1e-4
        beta = 3.0
        sparsityparameter=0.01
        theta = 0.1*randn(2*hidden*visible+hidden+visible)
        function justcost(x)
                c,g = autoencodercost(x,visible,hidden,traindata,lambda=lambda,sparsityparameter=sparsityparameter,beta=beta)
                return c
        end
        cost,grad = autoencodercost(theta,visible,hidden,traindata,lambda=lambda,sparsityparameter=sparsityparameter,beta=beta)
        numgrad = numericalgradient(justcost,theta)
        println([grad numgrad])
        diff = norm(numgrad - grad) / norm(numgrad + grad)^2
        return diff
end

function initializeparameters(hiddensize,visiblesize)
        r = sqrt(6) / sqrt(hiddensize+visiblesize+1)
        W1 = rand(hiddensize,visiblesize)*2*r - r
        W2 = rand(visiblesize,hiddensize)*2*r - r
        b1 = zeros(hiddensize)
        b2 = zeros(visiblesize)
        return [W1[:],W2[:],b1,b2]
end

using Plotly
export displaynetwork
function displaynetwork(A,filename)
        m,n = size(A)
        sz = int(sqrt(m))
        A -= mean(A)
        layout = [
            "autosize" => false,
            "width" => 500,
            "height"=> 500
        ]

        gridsize = int(ceil(sqrt(n)))
        buffer = 1
        griddata = ones(gridsize*(sz+1)+1,gridsize*(sz+1)+1)
        index = 1
        for i = 1:gridsize
                for j = 1:gridsize
                        if index > n
                                continue
                        end
                        columnlimit = maximum(abs(A[:,index]))
                        griddata[buffer+(i-1)*(sz+buffer)+(1:sz),buffer+(j-1)*(sz+buffer)+(1:sz)] = reshape(A[:,index],sz,sz)/columnlimit
                        index += 1
                end
        end

        Plotly.signin("kjchavez", "dd16t6j7li")
        data = [
          [
            "z" => griddata,
            "colorscale" => "Greys",
            "type" => "heatmap"
          ]
        ]
        response = Plotly.plot(data, ["layout" => layout, "filename" => filename, "fileopt" => "overwrite"])
        plot_url = response["url"]
end

using NLopt
export autoencode
function autoencode(data,hiddensize,visiblesize;lambda=1e-4,beta=3.0,rho=0.01,maxiter=600)
        theta::Vector{Float64} = initializeparameters(hiddensize,visiblesize)
        count = 0
        function objective(x,grad)
                count += 1
                println("iteration $count")
                return autoencodercost!(x,visiblesize,hiddensize,data, grad, lambda=lambda,sparsityparameter=rho,beta=beta)
        end
        opt = Opt(:LD_LBFGS,length(theta))
        min_objective!(opt,objective)
        maxeval!(opt,maxiter)
        minf,minx,ret = optimize(opt,theta)
        W1 = reshape(minx[1:hiddensize*visiblesize],hiddensize,visiblesize)
        W2 = reshape(minx[hiddensize*visiblesize+1:2*hiddensize*visiblesize],visiblesize,hiddensize)
        b1 = minx[2*hiddensize*visiblesize+1:2*hiddensize*visiblesize+hiddensize]
        b2 = minx[2*hiddensize*visiblesize+hiddensize+1:end]

        return minf,W1,W2,b1,b2
end

using MAT
const SAMPLE_WIDTH = 8
const SAMPLE_HEIGHT = 8
function sampleimages(images::Array{Float64,3},numsamples)
    width, height = size(images[:,:,1])
    array::Array{Float64,2} = zeros(SAMPLE_WIDTH*SAMPLE_HEIGHT,numsamples)
    for index=1:numsamples
        image_index = rand(1:size(images,3))
        x = rand(1:width-SAMPLE_WIDTH+1)
        y = rand(1:height-SAMPLE_HEIGHT+1)
        sample = images[x:x+SAMPLE_WIDTH-1,y:y+SAMPLE_HEIGHT-1,image_index]
        array[:,index] = reshape(sample,SAMPLE_WIDTH*SAMPLE_HEIGHT)
        array[:,index] -= mean(array[:,index]) #subtract mean
    end

    # rescale images to fit in range 0.1 to 0.9
    stddev = std(array)
    array = max(min(array,3*stddev),-3*stddev) / (3*stddev)
    array = (array + 1.0) * 0.4 + 0.1
    return array
end

function main()
        vars = matread("IMAGES.mat")
        traindata = sampleimages(vars["IMAGES"],10000)

        hiddensize = 25
        visiblesize = 64

        # tunable parameters
        lambda = 1e-4
        beta = 3.0
        rho = 0.01

        f,W1,W2,b1,b2 = autoencode(traindata,hiddensize,visiblesize)
        displaynetwork(W1',"autoencoder-tutorial")
end

end


text:
#######################################
#
# AVR file
#
#######################################

"""
Read AVR (.avr) file

#### Input
* `fname`: Name or path for the AVR file

#### Output
* `data`: Array of data read from AVR file. Each column represents a channel, and each row represents a point.
* `chanNames`: Channel Names

"""
function read_avr(fname::AbstractString)
    @info("Reading AVR file: $fname")

    # Open file
    file = open(fname, "r")

    # Header line
    header_exp =
        r"Npts= (\d*)\s+TSB= ([-+]?[0-9]*\.?[0-9]+)\s+DI= ([-+]?[0-9]*\.?[0-9]+)\s+SB= ([-+]?[0-9]*\.?[0-9]+)\s+SC= ([-+]?[0-9]*\.?[0-9]+)\s+Nchan= (\d*)"
    m = match(header_exp, readline(file))
    npts = parse(Int, ascii(m.captures[1])) # Number of points
    tsb = m.captures[2] # Omit?
    di = m.captures[3] # Omit?
    sb = m.captures[4] # Omit?
    sc = m.captures[5] # Omit?
    nchan = parse(Int, ascii(m.captures[6])) # Number of channels

    # Channel line
    names_exp = r"(\w+)"
    chanNames = collect((m.match for m in eachmatch(names_exp, readline(file))))

    # Data
    data = Array{Float64}(undef, (npts, nchan))
    for c = 1:nchan
        d = collect((m.match for m in eachmatch(r"([-+]?[0-9]*\.?[0-9]+)", readline(file))))
        for n = 1:npts
            data[n, c] = parse(Float64, ascii(d[n]))
        end
    end

    # Close file
    close(file)

    return data, chanNames
end

"""
Write AVR file
"""
function write_avr(fname::AbstractString, data::Array, chanNames::Array, fs::Number)

    @info("Saving avr to $fname")

    fs = float(fs)

    open(fname, "w") do fid
        @printf(
            fid,
            "Npts= %d   TSB= %2.6f DI= %2.6f SB= %2.3f SC= %3.1f Nchan= %d\n",
            size(data, 1),
            1000 / fs,
            1000 / fs,
            1.0,
            200.0,
            size(data, 2)
        )

        @printf(fid, "%s", chanNames[1])
        for c = 2:length(chanNames)
            @printf(fid, " ")
            @printf(fid, "%s", chanNames[c])
        end
        @printf(fid, "\n")

        for c = 1:size(data, 2)
            @printf(fid, "%2.6f", data[1, c])
            for p = 2:size(data, 1)
                @printf(fid, " ")
                @printf(fid, "%2.6f", data[p, c])
            end
            @printf(fid, " ")
            @printf(fid, "\n")
        end
    end
end


text:
#=
Annealed importance sampling (AIS) to estimate the partition function (and hence
the log-likelihood).
A nice explanation of AIS in general can be found in Goodfellow et al Deep Learning book.
Salakhutdinov et al (10.1145/1390156.1390266, http://www.cs.utoronto.ca/~rsalakhu/papers/bm.pdf)
discusses AIS for RBMs specifically.

AIS tends to understimate the log of the partition function (in probability).
In contrast, Reverse AIS estimator (RAISE) can be used to obtain a stochastic upper bound.
See http://proceedings.mlr.press/v38/burda15.html.
Combining the two we can "sandwiches" the true value to have an idea if the Monte Carlo
chains have converged.

Addendum: I think Burda's paper has a typo. The correct expression for the weights In
reverse AIS (which I use here) can be found in Upadhya et al 2015, Equation 10
(https://link.springer.com/chapter/10.1007/978-3-319-26535-3_62).

Bonus: A discussion of estimating partition function in RBMs, comparing several algorithms:

https://www.sciencedirect.com/science/article/pii/S0004370219301948

For a variant or RAISE: https://arxiv.org/abs/1511.02543
=#

"""
    ais(rbm0, rbm1, v0, βs)

Provided `v0` is an equilibrated sample from `rbm0`, returns `F` such that `mean(exp.(F))` is
an unbiased estimator of `Z1/Z0`, the ratio of partition functions of `rbm1` and `rbm0`.

!!! tip Use [`logmeanexp`](@ref)
    `logmeanexp(F)`, using the function `logmeanexp`[@ref] provided in this package,
    tends to give a better approximation of `log(Z1) - log(Z0)` than `mean(F)`.
"""
function ais(rbm0::RBM, rbm1::RBM, v::AbstractArray, βs::AbstractVector)
    @assert issorted(βs) && 0 == first(βs) ≤ last(βs) == 1
    F = free_energy(rbm0, v)
    for β in βs
        if iszero(β) || isone(β)
            continue
        else
            rbm = anneal(rbm0, rbm1; β)
            F -= free_energy(rbm, v)
            v = sample_v_from_v(rbm, v)
            F += free_energy(rbm, v)
        end
    end
    F -= free_energy(rbm1, v)
    return F
end

function ais(rbm0::RBM, rbm1::RBM, v0::AbstractArray; nbetas::Int=2)
    βs = range(0, 1, nbetas)
    return ais(rbm0, rbm1, v0, βs)
end

"""
    aise(rbm, [βs]; [nbetas], init=rbm.visible, nsamples=1)

AIS estimator of the log-partition function of `rbm`. It is recommended to fit `init` to
the single-site statistics of `rbm` (or the data).

!!! tip Use large `nbetas`
    For more accurate estimates, use larger `nbetas`. It is usually better to have
    large `nbetas` and small `nsamples`, rather than large `nsamples` and small `nbetas`.
"""
function aise(rbm::RBM, βs::AbstractVector{<:Real}; init::AbstractLayer=rbm.visible, nsamples::Int=1)
    rbm0 = anneal_zero(init, rbm)
    v0 = sample_from_inputs(init, Falses(size(init)..., nsamples))
    F = ais(rbm0, rbm, v0, βs)
    return F .+ log_partition_zero_weight(rbm0)
end

"""
    raise(rbm::RBM, βs; v, init)

Reverse AIS estimator of the log-partition function of `rbm`.
While `aise` tends to understimate the log of the partition function, `raise` tends to
overstimate it. `v` must be an equilibrated sample from `rbm`.

!!! tip Use [`logmeanexp`](@ref)
    If `F = raise(...)`, then `-logmeanexp(-F)`, using the function `logmeanexp`[@ref]
    provided in this package, tends to give a better approximation of `log(Z)` than `mean(F)`.

!!! tip Sandwiching the log-partition function
    If `Rf = aise(...)`, `Rr = raise(...)` are the AIS and reverse AIS estimators, we have the
    stochastic bounds `logmeanexp(Rf) ≤ log(Z) ≤ -logmeanexp(-Rr)`.
"""
function raise(rbm::RBM, βs::AbstractVector; v::AbstractArray, init::AbstractLayer=rbm.visible)
    rbm0 = anneal_zero(init, rbm)
    F = ais(rbm, rbm0, v, βs)
    return log_partition_zero_weight(rbm0) .- F
end

aise(rbm::RBM; nbetas::Int=10000, kw...) = aise(rbm, range(0, 1, nbetas); kw...)
raise(rbm::RBM; nbetas::Int=10000, kw...) = raise(rbm, range(0, 1, nbetas); kw...)

"""
    anneal(rbm0, rbm1; β)

Returns an RBM that interpolates between `rbm0` and `rbm1`.
Denoting by `E0(v, h)` and `E1(v, h)` the energies assigned by `rbm0` and `rbm1`,
respectively, the returned RBM assigns energies given by:

    E(v,h) = (1 - β) * E0(v) + β * E1(v, h)
"""
function anneal(rbm0::RBM, rbm1::RBM; β::Real)
    vis = anneal(rbm0.visible, rbm1.visible; β)
    hid = anneal(rbm0.hidden, rbm1.hidden; β)
    w = (1 - β) * rbm0.w + β * rbm1.w
    return RBM(vis, hid, w)
end

anneal(init::Binary, final::Binary; β::Real) = Binary(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ)
)
anneal(init::Spin, final::Spin; β::Real) = Spin(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ)
)
anneal(init::Potts, final::Potts; β::Real) = Potts(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ)
)
anneal(init::Gaussian, final::Gaussian; β::Real) = Gaussian(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ),
    oftype(final.γ, (1 - β) * init.γ + β * final.γ)
)
anneal(init::ReLU, final::ReLU; β::Real) = ReLU(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ),
    oftype(final.γ, (1 - β) * init.γ + β * final.γ)
)
anneal(init::dReLU, final::dReLU; β::Real) = dReLU(
    oftype(final.θp, (1 - β) * init.θp + β * final.θp),
    oftype(final.θn, (1 - β) * init.θn + β * final.θn),
    oftype(final.γp, (1 - β) * init.γp + β * final.γp),
    oftype(final.γn, (1 - β) * init.γn + β * final.γn)
)
anneal(init::pReLU, final::pReLU; β::Real) = pReLU(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ),
    oftype(final.γ, (1 - β) * init.γ + β * final.γ),
    oftype(final.Δ, (1 - β) * init.Δ + β * final.Δ),
    oftype(final.η, (1 - β) * init.η + β * final.η)
)
anneal(init::xReLU, final::xReLU; β::Real) = xReLU(
    oftype(final.θ, (1 - β) * init.θ + β * final.θ),
    oftype(final.γ, (1 - β) * init.γ + β * final.γ),
    oftype(final.Δ, (1 - β) * init.Δ + β * final.Δ),
    oftype(final.ξ, (1 - β) * init.ξ + β * final.ξ)
)

anneal_zero(init::AbstractLayer, rbm1::RBM) = RBM(init, anneal_zero(rbm1.hidden), Zeros(rbm1.w))

anneal_zero(l::Binary) = Binary(zero(l.θ))
anneal_zero(l::Spin) = Spin(zero(l.θ))
anneal_zero(l::Potts) = Potts(zero(l.θ))
anneal_zero(l::Gaussian) = Gaussian(zero(l.θ), l.γ)
anneal_zero(l::ReLU) = ReLU(zero(l.θ), l.γ)
anneal_zero(l::dReLU) = dReLU(zero(l.θp), zero(l.θn), l.γp, l.γn)
anneal_zero(l::pReLU) = pReLU(zero(l.θ), l.γ, zero(l.Δ), l.η)
anneal_zero(l::xReLU) = xReLU(zero(l.θ), l.γ, zero(l.Δ), l.ξ)

"""
    log_partition_zero_weight(rbm)

Log-partition function of a zero-weight version of `rbm`.
"""
log_partition_zero_weight(rbm::RBM) = -free_energy(rbm.visible) - free_energy(rbm.hidden)

"""
    logmeanexp(A; dims=:)

Computes `log.(mean(exp.(A); dims))`, in a numerically stable way.
"""
function logmeanexp(A::AbstractArray; dims=:)
    R = logsumexp(A; dims)
    N = length(A) ÷ length(R)
    return R .- log(N)
end

"""
    logvarexp(A; dims=:)

Computes `log.(var(exp.(A); dims))`, in a numerically stable way.
"""
function logvarexp(
    A::AbstractArray; dims=:, corrected::Bool=true, logmean=logmeanexp(A; dims)
)
    R = logsumexp(2logsubexp.(A, logmean); dims)
    N = length(A) ÷ length(R)
	if corrected
		return R .- log(N - 1)
    else
        return R .- log(N)
    end
end

"""
    logstdexp(A; dims=:)

Computes `log.(std(exp.(A); dims))`, in a numerically stable way.
"""
function logstdexp(
    A::AbstractArray; dims=:, corrected::Bool=true, logmean=logmeanexp(A; dims)
)
    return logvarexp(A; dims, corrected, logmean) / 2
end


text:
struct EmbeddedMap{T, As <: LinearMap, Rs <: AbstractVector{Int},
                    Cs <: AbstractVector{Int}} <: LinearMap{T}
    lmap::As
    dims::Dims{2}
    rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1)
    cols::Cs # typically j1:j2 with 1 <= j1 <= j2 <= size(map,2)

    function EmbeddedMap{T}(map::As, dims::Dims{2}, rows::Rs, cols::Cs) where {T,
                    As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}}
        check_index(rows, size(map, 1), dims[1])
        check_index(cols, size(map, 2), dims[2])
        return new{T,As,Rs,Cs}(map, dims, rows, cols)
    end
end

EmbeddedMap(map::LinearMap{T}, dims::Dims{2}; offset::Dims{2}) where {T} =
    EmbeddedMap{T}(map, dims, offset[1] .+ (1:size(map, 1)), offset[2] .+ (1:size(map, 2)))
EmbeddedMap(map::LinearMap, dims::Dims{2}, rows::AbstractVector{Int}, cols::AbstractVector{Int}) =
    EmbeddedMap{eltype(map)}(map, dims, rows, cols)

@static if VERSION >= v"1.8-"
    Base.reverse(A::LinearMap; dims=:) = _reverse(A, dims)
    function _reverse(A, dims::Integer)
        if dims == 1
            return EmbeddedMap(A, size(A), reverse(axes(A, 1)), axes(A, 2))
        elseif dims == 2
            return EmbeddedMap(A, size(A), axes(A, 1), reverse(axes(A, 2)))
        else
            throw(ArgumentError("invalid dims argument to reverse, should be 1 or 2, got $dims"))
        end
    end
    _reverse(A, ::Colon) = EmbeddedMap(A, size(A), map(reverse, axes(A))...)
    _reverse(A, dims::NTuple{1,Integer}) = _reverse(A, first(dims))
    function _reverse(A, dims::NTuple{M,Integer}) where {M}
        dimrev = ntuple(k -> k in dims, 2)
        if 2 < M || M != sum(dimrev)
            throw(ArgumentError("invalid dimensions $dims in reverse!"))
        end
        ax = ntuple(k -> dimrev[k] ? reverse(axes(A, k)) : axes(A, k), 2)
        return EmbeddedMap(A, size(A), ax...)
    end
end

function check_index(index::AbstractVector{Int}, dimA::Int, dimB::Int)
    length(index) != dimA && throw(ArgumentError("invalid length of index vector"))
    minimum(index) <= 0 && throw(ArgumentError("minimal index is below 1"))
    maximum(index) > dimB && throw(ArgumentError(
        "maximal index $(maximum(index)) exceeds dimension $dimB"
        ))
    # _isvalidstep(index) || throw(ArgumentError("non-monotone index set"))
    nothing
end

# _isvalidstep(index::AbstractRange) = step(index) > 0
# _isvalidstep(index::AbstractVector) = all(diff(index) .> 0)

Base.size(A::EmbeddedMap) = A.dims

# sufficient but not necessary conditions
LinearAlgebra.issymmetric(A::EmbeddedMap) =
    issymmetric(A.lmap) && (A.dims[1] == A.dims[2]) && (A.rows == A.cols)
LinearAlgebra.ishermitian(A::EmbeddedMap) =
    ishermitian(A.lmap) && (A.dims[1] == A.dims[2]) && (A.rows == A.cols)

Base.:(==)(A::EmbeddedMap, B::EmbeddedMap) =
    (eltype(A) == eltype(B)) && (A.lmap == B.lmap) &&
    (A.dims == B.dims) && (A.rows == B.rows) && (A.cols == B.cols)

LinearAlgebra.adjoint(A::EmbeddedMap) = EmbeddedMap(adjoint(A.lmap), reverse(A.dims), A.cols, A.rows)
LinearAlgebra.transpose(A::EmbeddedMap) = EmbeddedMap(transpose(A.lmap), reverse(A.dims), A.cols, A.rows)

for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
    @eval function _unsafe_mul!(y::$Out, A::EmbeddedMap, x::$In)
        fill!(y, zero(eltype(y)))
        _unsafe_mul!(selectdim(y, 1, A.rows), A.lmap, selectdim(x, 1, A.cols))
        return y
    end
    @eval function _unsafe_mul!(y::$Out, A::EmbeddedMap, x::$In, alpha::Number, beta::Number)
        LinearAlgebra._rmul_or_fill!(y, beta)
        _unsafe_mul!(selectdim(y, 1, A.rows), A.lmap, selectdim(x, 1, A.cols), alpha, !iszero(beta))
        return y
    end
end

function _unsafe_mul!(Y::AbstractMatrix, A::EmbeddedMap, x::Number)
    fill!(Y, zero(eltype(Y)))
    _unsafe_mul!(view(Y, A.rows, A.cols), A.lmap, x)
    return Y
end
function _unsafe_mul!(Y::AbstractMatrix, A::EmbeddedMap, x::Number, α::Number, β::Number)
    LinearAlgebra._rmul_or_fill!(Y, β)
    _unsafe_mul!(view(Y, A.rows, A.cols), A.lmap, x, α, !iszero(β))
    return Y
end


text:
"""
  calculate_auc(ya, yb)

  # "Calcualtes AUC related estimates"

  This function takes two DataArray arguments `ya` and `yb` and calculates variance of predicted AUC,
  logit of predicted AUC, and variance of logit of predicted AUC responses passed.
"""
calculate_auc(ya = DataFrames.DataArray([2,3,4,3]), yb = DataFrames.DataArray([3,2,1,3,4,3]))

function calculate_auc(; ya::DataArrays.DataArray = nothing, yb::DataArrays.DataArray = nothing)
  m = length(ya)
  p = length(yb)
  I = zeros(m, p)
    for i in range(1,m)
        for j in range(1,p)
            if ya[i] > yb[j]
              I[i,j] = 1
            elseif ya[i] == yb[j]
              I[i,j] = 0.5
            else
               I[i,j] = 0
            end
        end
    end
    finv(x::Float64) = return(-log((1/x)-1))
    auchat = mean(I)
    finvhat = finv(auchat)
    vya = mean(I,2)
    vyb = mean(I,1)
    svarya = var(vya)
    svaryb = var(vyb)
    vhat_auchat = (svarya/m) + (svaryb/p)
    v_finv_auchat = vhat_auchat/((auchat^2)*(1-auchat)^2)
    logitauchat = log(auchat/(1-auchat))
    var_logitauchat = vhat_auchat /((auchat^2)*(1-auchat)^2)
    return((var_logitauchat, logitauchat))
end


text:
#################### Slice Sampler ####################

#################### Types and Constructors ####################

const SliceForm = Union{Univariate, Multivariate}

mutable struct SliceTune{F<:SliceForm} <: SamplerTune
  logf::Union{Function, Missing}
  width::Union{Float64, Vector{Float64}}

  SliceTune{F}() where {F<:SliceForm} = new{F}()

  SliceTune{F}(x::Vector, width) where {F<:SliceForm} =
    SliceTune{F}(x, width, missing)

  SliceTune{F}(x::Vector, width::Real, logf::Union{Function, Missing}) where
    {F<:SliceForm} = new{F}(logf, Float64(width))

  SliceTune{F}(x::Vector, width::Vector, logf::Union{Function, Missing}) where
    {F<:SliceForm} = new{F}(logf, convert(Vector{Float64}, width))
end

const SliceUnivariate = SamplerVariate{SliceTune{Univariate}}
const SliceMultivariate = SamplerVariate{SliceTune{Multivariate}}

validate(v::SamplerVariate{SliceTune{F}}) where {F<:SliceForm} =
  validate(v, v.tune.width)

validate(v::SamplerVariate{SliceTune{F}}, width::Float64) where {F<:SliceForm} = v

function validate(v::SamplerVariate{SliceTune{F}}, width::Vector) where {F<:SliceForm}
  n = length(v)
  length(width) == n ||
    throw(ArgumentError("length(width) differs from variate length $n"))
  v
end

#################### Sampler Constructor ####################
"""
Slice(params::ElementOrVector{Symbol},
                width::ElementOrVector{T},
                ::Type{F}=Multivariate;
                transform::Bool=false) where {T<:Real, F<:SliceForm}

Construct a `Sampler` object for Slice sampling. Parameters are assumed to be
continuous, but may be constrained or unconstrained.

Returns a `Sampler{SliceTune{Univariate}}` or `Sampler{SliceTune{Multivariate}}`
type object if sampling univariately or multivariately, respectively.

  * `params`: stochastic node(s) to be updated with the sampler.

  * `width`: scaling value or vector of the same length as the combined elements of nodes `params`, defining initial widths of a hyperrectangle from which to simulate values.

  * `F` : sampler type. Options are
      * `:Univariate` : sequential univariate sampling of parameters.
      * `:Multivariate` : joint multivariate sampling.

  * `transform`: whether to sample parameters on the link-transformed scale (unconstrained parameter space). If `true`, then constrained parameters are mapped to unconstrained space according to transformations defined by the Stochastic `unlist()` function, and `width` is interpreted as being relative to the unconstrained parameter space. Otherwise, sampling is relative to the untransformed space.
"""
function Slice(params::ElementOrVector{Symbol},
                width::ElementOrVector{T},
                ::Type{F}=Multivariate;
                transform::Bool=false) where {T<:Real, F<:SliceForm}
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block, transform)
    v = SamplerVariate(block, width)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, SliceTune{F}())
end

#################### Sampling Functions ####################

sample!(v::Union{SliceUnivariate, SliceMultivariate}) = sample!(v, v.tune.logf)
"""
    sample!(v::Union{SliceUnivariate, SliceMultivariate}, logf::Function)

Draw one sample from a target distribution using the Slice univariate or
multivariate sampler. Parameters are assumed to be continuous, but may be
constrained or unconstrained.

Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::Union{SliceUnivariate, SliceMultivariate}, logf::Function)
    typeof(v.value[1]) <:GeneralNode ? sample_node!(v, logf) : sample_number!(v, logf)
end

function sample_node!(v::SliceUnivariate, logf::Function)
  tree = v.value[1]

  logf0 = logf(tree)
  blv = get_branchlength_vector(tree)

  n = length(blv)
  lower = blv - v.tune.width .* rand(n)
  lower[lower .< 0.0] .= 0.0
  upper = lower .+ v.tune.width

  for i in 1:n
    p0 = logf0 + log(rand())

    x = blv[i]
    blv[i] = rand(Uniform(lower[i], upper[i]))
    while true
      set_branchlength_vector!(tree, blv)
      logf0 = logf(tree)
      logf0 < p0 || break
      value = blv[i]
      if value < x
        lower[i] = value
      else
        upper[i] = value
      end
      blv[i] = rand(Uniform(lower[i], upper[i]))
    end
  end

  v
end

function sample_node!(v::SliceMultivariate, logf::Function)
  tree = v.value[1]

  p0 = logf(tree) + log(rand())
  blv = get_branchlength_vector(tree)
  org = deepcopy(blv)

  n = length(blv)
  lower = blv - v.tune.width .* rand(n)
  lower[lower .< 0.0] .= 0.0
  upper = lower .+ v.tune.width

  blv = v.tune.width .* rand(n) + lower
  set_branchlength_vector!(tree, blv)
  while logf(tree) < p0
    for i in 1:n
      value = blv[i]
      if value < org[i]
        lower[i] = value
      else
        upper[i] = value
      end
      blv[i] = rand(Uniform(lower[i], upper[i]))
    end
    set_branchlength_vector!(tree, blv)
  end

  v.value[1] = tree
  v
end

function sample_number!(v::SliceUnivariate, logf::Function)
  logf0 = logf(v.value)

  n = length(v)
  lower = v - v.tune.width .* rand(n)
  upper = lower .+ v.tune.width

  for i in 1:n
    p0 = logf0 + log(rand())

    x = v[i]
    v[i] = rand(Uniform(lower[i], upper[i]))
    while true
      logf0 = logf(v.value)
      logf0 < p0 || break
      value = v[i]
      if value < x
        lower[i] = value
      else
        upper[i] = value
      end
      v[i] = rand(Uniform(lower[i], upper[i]))
    end
  end

  v
end

function sample_number!(v::SliceMultivariate, logf::Function)
  p0 = logf(v.value) + log(rand())

  n = length(v)
  lower = v - v.tune.width .* rand(n)
  upper = lower .+ v.tune.width

  x = v.tune.width .* rand(n) + lower
  while logf(x) < p0
    for i in 1:n
      value = x[i]
      if value < v[i]
        lower[i] = value
      else
        upper[i] = value
      end
      x[i] = rand(Uniform(lower[i], upper[i]))
    end
  end
  v[:] = x

  v
end


text:
"""
    simulate linear feedback policy
"""
function simulate_linear_feedback(
        model,
        K,
        x̄, ū,
        x_ref, u_ref,
        Q, R,
        T_sim, Δt,
        x1,
        w;
        ul = -Inf * ones(length(ū[1])),
        uu = Inf * ones(length(ū[1])))

    T = length(x̄)
    times = [(t - 1) * Δt for t = 1:T-1]
    tf = Δt * (T-1)
    t_sim = range(0, stop = tf, length = T_sim)
    t_ctrl = range(0, stop = tf, length = T)
    dt_sim = tf / (T_sim - 1)

    x_sim = [x1]
    u_sim = []

    J = 0.0
    Jx = 0.0
    Ju = 0.0

    A_state = Array(hcat(x̄...))
    A_ctrl = Array(hcat(ū...))

    for tt = 1:T_sim-1
        t = t_sim[tt]
        k = searchsortedlast(times, t)
        x = x_sim[end]

        x_cubic = Array(zero(x̄[1]))

        for i = 1:length(x̄[1])
            interp_cubic = CubicSplineInterpolation(t_ctrl, A_state[i,:])
            x_cubic[i] = interp_cubic(t)
        end

        u = ū[k] + K[k] * (x - x_cubic)

        # clip controls
        u = max.(u, ul)
        u = min.(u, uu)

        push!(x_sim, fd(model, x, u, w[k], dt_sim, tt))
        push!(u_sim, u)

        J += (x_sim[end] - x_ref[k])' * Q[k + 1] * (x_sim[end] - x_ref[k])
        J += (u_sim[end] - u_ref[k])' * R[k] * (u_sim[end] - u_ref[k])
        Jx += (x_sim[end] - x_ref[k])' * Q[k + 1] * (x_sim[end] - x_ref[k])
        Ju += (u_sim[end] - u_ref[k])' * R[k] * (u_sim[end] - u_ref[k])
    end

    return x_sim, u_sim, J / (T_sim - 1), Jx / (T_sim - 1), Ju / (T_sim - 1)
end

function simulate_policy(
        model,
        θ,
        x_ref, u_ref,
        Q, R,
        T_sim, Δt,
        x1,
        w;
        ul = -Inf * ones(length(ū[1])),
        uu = Inf * ones(length(ū[1])))

    T = length(x_ref)
    times = [(t - 1) * Δt for t = 1:T-1]
    tf = Δt * (T-1)
    t_sim = range(0, stop = tf, length = T_sim)
    t_ctrl = range(0, stop = tf, length = T)
    dt_sim = tf / (T_sim - 1)

    x_sim = [x1]
    u_sim = []

    J = 0.0
    Jx = 0.0
    Ju = 0.0

    # A_state = Array(hcat(x̄...))
    # A_ctrl = Array(hcat(ū...))

    for tt = 1:T_sim-1
        t = t_sim[tt]
        k = searchsortedlast(times, t + 1.0e-16)
        x = x_sim[end]

        # println("t $t")
        # println("k $k")
        # println("")

        # x_cubic = Array(zero(x̄[1]))
        #
        # for i = 1:length(x̄[1])
        #     interp_cubic = CubicSplineInterpolation(t_ctrl, A_state[i,:])
        #     x_cubic[i] = interp_cubic(t)
        # end

        # # clip controls
        # u = max.(u, ul)
        # u = min.(u, uu)
        push!(x_sim, fd(model, x, [zeros(model.m); θ[k]], w[k], dt_sim, tt))
        push!(u_sim, policy(θ[k], x, nothing, model.n, model.m))

        J += (x_sim[end] - x_ref[k])' * Q[k + 1] * (x_sim[end] - x_ref[k])
        J += (u_sim[end] - u_ref[k])' * R[k] * (u_sim[end] - u_ref[k])
        Jx += (x_sim[end] - x_ref[k])' * Q[k + 1] * (x_sim[end] - x_ref[k])
        Ju += (u_sim[end] - u_ref[k])' * R[k] * (u_sim[end] - u_ref[k])
    end

    return x_sim, u_sim, J / (T_sim - 1), Jx / (T_sim - 1), Ju / (T_sim - 1)
end