text:
"""
    ArbComplex(x)
    ArbComplex(x; [bits])
    ArbComplex(x; [digits [, base=10]])

    ArbComplex(re, im)
    ArbComplex(re, im; [bits])
    ArbComplex(re, im; [digits [, base=10]])

- `ArbComplex`(re, im) == `ArbComplex`(re, im, bits=workingprecision(ArbComplex))
- `ArbComplex`(re, im, digits=digits, base=2) == `ArbComplex`(re, im, bits=digits)
- `ArbComplex`(re, im, digits=digits, base=10) == `ArbComplex`(re, im, digits=digits)

`ArbComplexs` are arbitrary precision complex binary floating-point regions with some specific precision established at construction.

The precision of an `ArbComplex` must be >= 24 bits or >= 8 digits.

Internally, an `ArbComplex` represents as a pair of `ArbReals`, where the real part and the imaginary part are each
intervals given by midpoint and radius (half-width of the interval). This representation serves as a _ball_ over the complex numbers.
Calculations with `ArbComplex` generate results that are assured to contain the true mathematical result.
There is no a priori assurance on the width of the interval that results.
Good practice is to set the precision substantively greater than the precision
of the resultant width required, and to check the radius of results.

See also: [`ArbFloat`](@ref), [`ArbReal`](@ref), [`ball`](@ref), [`setball`](@ref), [`midpoint`](@ref), [`radius`](@ref)
""" ArbComplex

# ArbComplex structs hold complex balls given as ArbReal pairs

mutable struct ArbComplex{P}  <: Number  # P is the precision in bits
                         #      real midpoint
    real_mid_exp::Int    # fmpz         exponent of 2 (2^exp)
    real_mid_size::UInt  # mp_size_t    nwords and sign (lsb holds sign of significand)
    real_mid_d1::UInt    # significand  unsigned, immediate value or the initial span
    real_mid_d2::UInt    #   (d1, d2)   the final part indicating the significand, or 0
                         #      real radius
    real_rad_exp::Int    # fmpz       exponent of 2 (2^exp)
    real_rad_man::UInt   # mp_limb_t  radius, unsigned by definition
                         #      imaginary midpoint
    imag_mid_exp::Int    # fmpz         exponent of 2 (2^exp)
    imag_mid_size::UInt  # mp_size_t    nwords and sign (lsb holds sign of significand)
    imag_mid_d1::UInt    # significand  unsigned, immediate value or the initial span
    imag_mid_d2::UInt    #   (d1, d2)   the final part indicating the significand, or 0
                         #      imaginary radius
    imag_rad_exp::Int    # fmpz       exponent of 2 (2^exp)
    imag_rad_man::UInt   # mp_limb_t  radius, unsigned by definition

    function ArbComplex{P}() where {P}
        z = new{P}(0,0,0,0,0,0,0,0,0,0,0,0)
        ccall(@libarb(acb_init), Cvoid, (Ref{ArbComplex},), z)
        finalizer(clear_acb, z)
        return z
    end
end

ArbComplex{P}(x::T) where {P, T<:Number} = ArbComplex{P}(real(x), imag(x))

# for use within a struct, eg. ArbComplexMatrix
const PtrToArbComplex = Ptr{ArbComplex} # acb_ptr
const PtrToPtrToArbComplex = Ptr{Ptr{ArbComplex}} # acb_ptr*

clear_acb(x::ArbComplex{P}) where {P} = ccall(@libarb(acb_clear), Cvoid, (Ref{ArbComplex},), x)

ArbComplex{P}(x::ArbComplex{P}) where {P} = x
ArbComplex(x::ArbComplex{P}) where {P} = x

ArbComplex{P}(x::Missing) where {P} = missing
ArbComplex(x::Missing) = missing

function ArbComplex(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real}
    bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2"))
    digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10"))
    if base === 10
        bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x)
    elseif base === 2
        bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbComplex{bits}(x, zero(T))
end

function ArbComplex(x::T, y::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Real}
    bits > 0 && bits < MINIMUM_PRECISION_BASE2 && throw(DomainError("bit precision $bits < $MINIMUM_PRECISION_BASE2"))
    digits > 0 && digits < MINIMUM_PRECISION_BASE10 && throw(DomainError("digit precision $digits < $MINIMUM_PRECISION_BASE10"))
    if base === 10
        bits = digits > 0 ? bits4digits(digits)+extrabits() : (bits > 0 ? bits+extrabits() : DEFAULT_PRECISION.x)
    elseif base === 2
        bits = bits > 0 ? bits+extrabits() : (digits > 0 ? digits+extrabits() : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbComplex{bits}(x, y)
end

function ArbComplex(x::T1, y::T2; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T1<:Number, T2<:Number}
    ArbComplex(promote(x, y)...,)
end

ArbComplex(re::T, im::T) where T<:AbstractFloat = ArbComplex(ArbFloat(re), ArbFloat(im))

const Analytic = Cint(0) # prefer the non-analytic versions

@inline sign_bit(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size)
@inline sign_bits(x::ArbComplex{P}) where {P} = isodd(x.real_mid_size), isodd(x.imag_mid_size)
@inline sign_bit(x::ArbComplex{P}, ::Type{RealPart}) where {P} = isodd(x.real_mid_size)
@inline sign_bit(x::ArbComplex{P}, ::Type{ImagPart}) where {P} = isodd(x.imag_mid_size)

# ArbComplex{P}(sqrt(2.0))
function ArbComplex{P}(rea::Float64) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_d), Cvoid, (Ref{ArbComplex}, Float64), z, rea)
    return z
end

function ArbComplex{P}(rea::Cint) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_si), Cvoid, (Ref{ArbComplex}, Cint), z, rea)
    return z
end

function ArbComplex{P}(rea::ArbFloat{P}) where {P}
    x = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, rea)
    y = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), y, x)
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_arb), Cvoid, (Ref{ArbComplex}, Ref{ArbReal}), z, y)
    return z
end

function ArbComplex{P}(rea::ArbReal{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_arb), Cvoid, (Ref{ArbComplex}, Ref{ArbReal}), z, rea)
    return z
end

# ArbComplex{P}(sqrt(2.0), inv(17.0))
function ArbComplex{P}(re::Float64, im::Float64) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_d_d), Cvoid, (Ref{ArbComplex}, Float64, Float64), z, re, im)
    return z
end

function ArbComplex{P}(x::Int64, y::Int64) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_si_si), Cvoid, (Ref{ArbComplex}, Cint, Cint, Clong), z, x, y, P)
    return z
end

function copy(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}), z, x)
    return z
end

deepcopy(x::ArbComplex{P}) where {P} = copy(x)

function ArbComplex{P}(x::ArbComplex{Q}, roundingmode::RoundingMode) where {P,Q}
    z = ArbComplex{P}()
    rounding = match_rounding_mode(roundingmode)
    res = ccall(@libarb(acb_set_round), Cint,
                (Ref{ArbComplex}, Ref{ArbComplex}, Clong, Cint), z, x, P, rounding)
    return z
end

ArbComplex(x::ArbFloat{P}) where {P} = ArbComplex{P}(x)
ArbComplex(x::ArbReal{P}) where {P} = ArbComplex{P}(x)
ArbComplex(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ArbComplex{P}(x,y)
ArbComplex(x::ArbReal{P}, y::ArbReal{P}) where {P} = ArbComplex{P}(x,y)
ArbComplex(x::ArbFloat{P}, y::ArbReal{P}) where {P} = ArbComplex{P}(ArbReal{P}(x),y)
ArbComplex(x::ArbReal{P}, y::ArbFloat{P}) where {P} = ArbComplex{P}(x,ArbReal{P}(y))

Base.Complex(x::ArbFloat{P}) where {P} = ArbComplex(x)
Base.Complex(x::ArbReal{P}) where {P} = ArbComplex(x)
Base.Complex(re::ArbFloat{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, im)
Base.Complex(re::ArbReal{P}, im::ArbReal{P}) where {P} = ArbComplex(re, im)

Base.complex(re::ArbFloat{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, im)
Base.complex(re::ArbReal{P}, im::ArbReal{P}) where {P} = ArbComplex(re, im)

Base.Complex(re::ArbFloat{P}, im::ArbReal{P}) where {P} = ArbComplex(ArbReal{P}(re), im)
Base.Complex(re::ArbReal{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, ArbReal{P}(im))
Base.complex(re::ArbFloat{P}, im::ArbReal{P}) where {P} = ArbComplex(ArbReal{P}(re), im)
Base.complex(re::ArbReal{P}, im::ArbFloat{P}) where {P} = ArbComplex(re, ArbReal{P}(im))

Base.Complex{ArbFloat{P}}(x::ArbComplex{P}) where {P} = x
Base.Complex{ArbReal{P}}(x::ArbComplex{P}) where {P} = x

function ArbComplex{P}(x::T) where {P,T<:Integer}
    y = ArbReal{P}(x)
    z = ArbComplex{P}(y)
    return z
end

function ArbComplex{P}(x::T1, y::T2) where {P,T1<:Union{Integer,AbstractFloat},T2<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z  = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::Irrational{S1}, y::Irrational{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z  = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::Irrational{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::T, y::Irrational{S}) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::T, y::ArbFloat{P}) where {P,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

ArbFloat{P}(x::ArbComplex{P}, bits::Int) where P = ArbFloat(real(x), bits=bits)
ArbReal{P}(x::ArbComplex{P}, bits::Int) where P = ArbReal(real(x), bits=bits)
ArbComplex{P}(x::ArbReal{P}, bits::Int) where P = ArbComplex(x, bits=bits)
ArbComplex{P}(x::ArbFloat{P}, bits::Int) where P = ArbComplex(x, bits=bits)

function ArbComplex{P}(x::ArbFloat{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function ArbComplex{P}(x::T, y::ArbReal{P}) where {P,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    z = ArbComplex{P}(x1, y)
    return z
end
function ArbComplex{P}(x::ArbReal{P}, y::T) where {P,T<:Union{Integer,AbstractFloat}}
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x, y1)
    return z
end

function ArbComplex{P}(x::ArbFloat{P}, y::ArbReal{P}) where {P}
    x1 = ArbReal{P}(x)
    z = ArbComplex{P}(x1, y)
    return z
end
function ArbComplex{P}(x::ArbReal{P}, y::ArbFloat{P}) where {P}
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x, y1)
    return z
end

function ArbComplex{P}(x::T, y::ArbFloat{S}) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::T, y::ArbReal{S}) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbReal{S}, y::T) where {P,S,T<:Union{Integer,AbstractFloat}}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{S1}, y::ArbFloat{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbReal{S1}, y::ArbReal{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbFloat{S1}, y::ArbReal{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end
function ArbComplex{P}(x::ArbReal{S1}, y::ArbFloat{S2}) where {P,S1,S2}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    z = ArbComplex{P}(x1, y1)
    return z
end

function conj(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_conj), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}), z, x)
    return z
end

function midpoint(x::ArbComplex{P}) where {P}
    xreal = real(x)
    ximag = imag(x)
    zreal = ArbReal{P}()
    zimag = ArbReal{P}()
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zreal, xreal)
    ccall(@libarb(arb_get_mid_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zimag, ximag)
    return ArbComplex{P}(zreal, zimag)
end

function radius(x::ArbComplex{P}) where {P}
    xreal = real(x)
    ximag = imag(x)
    zreal = ArbReal{P}()
    zimag = ArbReal{P}()
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zreal, xreal)
    ccall(@libarb(arb_get_rad_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), zimag, ximag)
    return ArbComplex{P}(zreal, zimag)
end

trunc(x::ArbComplex{P}) where {P} = ArbComplex{P}(trunc(real(x)), trunc(imag(x)))
trunc(::Type{T}, x::ArbComplex{P}) where {P, T<:Integer} = T(trunc(real(x))), T(trunc(imag(x)))

floor(x::ArbComplex{P}) where {P} = ArbComplex{P}(floor(real(x)), floor(imag(x)))
floor(::Type{T}, x::ArbComplex{P}) where {P, T<:Integer} = T(floor(real(x))), T(floor(imag(x)))

ceil(x::ArbComplex{P}) where {P} = ArbComplex{P}(ceil(real(x)), ceil(imag(x)))
ceil(::Type{T}, x::ArbComplex{P}) where {P, T<:Integer} = T(ceil(real(x))), T(ceil(imag(x)))

function modf(x::ArbComplex{P}) where {P}
    repart = modf(real(x))
    impart = modf(imag(x))
    return (repart, impart)
end

fmod(repart::Tuple{ArbReal{P1}, ArbReal{P1}}, impart::Tuple{ArbReal{P2},ArbReal{P2}}) where {P1, P2} =
    ArbComplex(fmod(repart), fmod(impart))

# phase angle
function Base.angle(x::ArbComplex{P}) where {P}
    rea, ima = reim(x / hypot(reim(x)...,))

    y = hypot(rea - 1.0, ima)
    x = hypot(rea + 1.0, ima)

    a = 2 * atan(y , x)

    T = ArbFloat{P}
    !(signbit(a) || signbit(T(pi) - a)) ? a : (signbit(a) ? zero(T) : T(pi))
end

# needed for GenericLinearAlgebra

flipsign(x::ArbComplex{P}, y::T) where {P, T<:Base.IEEEFloat} =
    signbit(y) ? -x : x
flipsign(x::ArbComplex{P}, y::T) where {P, T<:Real} =
    signbit(y) ? -x : x
flipsign(x::ArbComplex{P}, y::T) where {P, T<:ArbFloat} =
    signbit(y) ? -x : x
flipsign(x::ArbComplex{P}, y::T) where {P, T<:ArbReal} =
    signbit(y) ? -x : x

copysign(x::ArbComplex{P}, y::T) where {P, T<:Base.IEEEFloat} =
    signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x)
copysign(x::ArbComplex{P}, y::T) where {P, T<:Real} =
    signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x)
copysign(x::ArbComplex{P}, y::T) where {P, T<:ArbFloat} =
    signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x)
copysign(x::ArbComplex{P}, y::T) where {P, T<:ArbReal} =
    signbit(y) ? (signbit(x) ? x : -x) : (signbit(x) ? -x : x)

# a type specific hash function helps the type to 'just work'
const hash_arbcomplex_lo = (UInt === UInt64) ? 0x76143ad985246e79 : 0x5b6a64dc
const hash_0_arbcomplex_lo = hash(zero(UInt), hash_arbcomplex_lo)
Base.hash(z::ArbComplex{P}, h::UInt) where {P} =
    hash(z.real_mid_d1 ⊻ z.real_rad_exp ⊻ z.imag_mid_d1 ⊻ z.imag_rad_exp,
           h ⊻ hash(z.real_mid_d2 ⊻ (~reinterpret(UInt,P)), hash_arbcomplex_lo)
               ⊻ hash_0_arbcomplex_lo)


text:
using Random: Random

struct Integrated{T,F}
    gen::F
end
function Integrated(m::Generator{T,F}) where {T,F}
    gen(rng) = unfold(Shuffle ∘ shrink, generate(rng, m))
    Integrated{Tree{T},typeof(gen)}(gen)
end
function Integrated(el::T) where T
    gen(_) = unfold(Shuffle ∘ shrink, el)
    Integrated{Tree{T},typeof(gen)}(gen)
end
Integrated(::Type{T}) where T = Integrated(Generator(T))
generate(rng, i::Integrated) = i.gen(rng)

function Base.iterate(g::Integrated, state=nothing)
    el = generate(default_rng(), g)
    el === nothing && return nothing
    (el, nothing)
end
Base.IteratorEltype(::Type{<:Integrated}) = Base.HasEltype()
Base.IteratorSize(::Type{<:Integrated}) = Base.SizeUnknown()
Base.eltype(::Type{Integrated{T,F}}) where {T,F} = T

################################################
# utility for working with integrated generators
################################################

const iBool = Integrated(Generator(Bool))
iWord(hi) = Integrated(mWord(hi))

freeze(i::Integrated{T,F}) where {T,F} = Generator{T}(i.gen)
# TODO: check if this should just produce the root continously instead
dontShrink(i::Integrated{T,F}) where {T,F} = Generator{T}((rng) -> root(i.gen(rng)))
dependent(g::Generator{T,F}) where {T,F} = Integrated{T,F}(g.gen)

function listAux(genLen::Integrated, genA::Integrated{T, F}) where {T, F}
    n = dontShrink(genLen)
    genF(rng) = [ generate(rng, freeze(genA)) for _ in 1:generate(rng, n) ]
    Generator{Vector{T}}(genF)
end
listAux(genLen, genA) = listAux(Integrated(genLen), genA)

# TODO: this should be possible to make more generic
function vector(genLen, genA::Integrated{T,F}) where {T,F}
    genF(rng) = interleave(generate(rng, listAux(genLen, genA)))
    gen = Generator{Vector{eltype(T)}}(genF)
    dependent(gen)
end

function tuple(genLen, genA::Integrated{T,F}) where {T,F}
    genF(rng) = interleave((generate(rng, listAux(genLen, genA))...,))
    gen = Generator{NTuple{N, eltype(T)} where N}(genF)
    dependent(gen)
end

function interleave(integrated...)
    T = typeof(map((root ∘ generate), integrated))
    genF(rng) = interleave(map(i -> generate(rng, i), integrated))
    gen = Generator{Tree{T}}(genF)
    dependent(gen)
end

function Base.filter(p, genA::Integrated{T,F}, trim=false) where {T,F}
    function genF(rng)
        element = iterate(filter(p, generate(rng, freeze(genA)), trim))
        element === nothing && return nothing
        return first(element)
    end
    gen = Generator{Union{Nothing, Tree{T}}}(genF)
    dependent(gen)
end

ext: jl


text:
const EAGER_INIT = Ref{Bool}(false)
const EAGER_THUNK_CHAN = Channel(typemax(Int))
const EAGER_THUNK_MAP = Dict{Int,Int}()
const EAGER_CONTEXT = Ref{Context}()

eager_context() = isassigned(EAGER_CONTEXT) ? EAGER_CONTEXT[] : nothing

function init_eager()
    EAGER_INIT[] && return
    EAGER_INIT[] = true
    if eager_context() === nothing
        EAGER_CONTEXT[] = Context([myid(),workers()...])
    end
    ctx = EAGER_CONTEXT[]
    @async try
        sopts = SchedulerOptions(;allow_errors=true)
        topts = ThunkOptions(;single=1)
        Dagger.compute(ctx, Dagger.delayed(eager_thunk;options=topts)(); options=sopts)
    catch err
        iob = IOContext(IOBuffer(), :color=>true)
        println(iob, "Error in eager scheduler:")
        Base.showerror(iob, err)
        Base.show_backtrace(iob, catch_backtrace())
        println(iob)
        seek(iob.io, 0)
        write(stderr, iob)
    finally
        EAGER_INIT[] = false
    end
end

"Adjusts the scheduler's cached pressure indicator for the specified worker by
the specified amount."
function adjust_pressure!(h::SchedulerHandle, proctype::Type, pressure)
    uid = Dagger.get_tls().sch_uid
    lock(ACTIVE_TASKS_LOCK) do
        ACTIVE_TASKS[uid][proctype][] += pressure
        notify(TASK_SYNC)
    end
    exec!(_adjust_pressure!, h, myid(), proctype, pressure)
end
function _adjust_pressure!(ctx, state, task, tid, (pid, proctype, pressure))
    state.worker_pressure[pid][proctype] += pressure
    nothing
end

"Allows a thunk to safely wait on another thunk, by temporarily reducing its
effective pressure to 0."
function thunk_yield(f)
    if Dagger.in_thunk()
        h = sch_handle()
        tls = Dagger.get_tls()
        proctype = typeof(tls.processor)
        util = tls.utilization
        adjust_pressure!(h, proctype, -util)
        try
            f()
        finally
            adjust_pressure!(h, proctype, util)
        end
    else
        f()
    end
end

function eager_thunk()
    h = sch_handle()
    util = Dagger.get_tls().utilization
    # Don't apply pressure from this thunk
    adjust_pressure!(h, Dagger.ThreadProc, -util)
    while isopen(EAGER_THUNK_CHAN)
        try
            future, uid, f, args, opts = take!(EAGER_THUNK_CHAN)
            args = map(x->x isa Dagger.EagerThunk ? ThunkID(EAGER_THUNK_MAP[x.uid]) : x, args)
            tid = add_thunk!(f, h, args...; opts...)
            register_future!(h, tid, future)
            EAGER_THUNK_MAP[uid] = tid.id
        catch err
            iob = IOContext(IOBuffer(), :color=>true)
            println(iob, "Error in eager listener:")
            Base.showerror(iob, err)
            Base.show_backtrace(iob, catch_backtrace())
            println(iob)
            seek(iob.io, 0)
            write(stderr, iob)
        end
    end
end


text:
# Normalization: remove trailing zeros in coefficient
function normalize(x::Decimal; rounded::Bool=false)
    p = 0
    if x.c != 0
        while x.c % 10^(p+1) == 0
            p += 1
        end
    end
    c = BigInt(x.c / 10^p)
    q = (c == 0 && x.s == 0) ? 0 : x.q + p
    if rounded
        Decimal(x.s, abs(c), q)
    else
        round(Decimal(x.s, abs(c), q), digits=DIGITS, normal=true)
    end
end


text:
# This file is a part of StruckVMEDevices.jl, licensed under the MIT License (MIT).

__precompile__(true)

module StruckVMEDevices

include("MemRegisters/MemRegisters.jl")
include("VMEGateways/VMEGateways.jl")
include("SIS3316Digitizers/SIS3316Digitizers.jl")

end # module


text:
#=

A call to Resynthesizer using large image.

For testing mallocs.

Usage:

>julia
julia> include("large.jl")
julia> @time test(image, mask)

OR

>julia --track-allocation=user
and the usual process for analyzing mallocs...
julia> include("large.jl")
julia> test(image, mask)
julia>using Profile
julia>Profile.clear_malloc_data()
julia> test(image, mask)
=#

# For "load" local  files
using FileIO
# The tested item
include("../../src/apps/inpaint.jl")

# mock test data
# global so it persists and out of the test
image = load("/home/bootch/git/resynthesizerJulia/test/data/in/lighthouse.png")
# large mask
mask = falses(size(image))
mask[50:350,300:400] .= true

function test(image, mask)
    inpaint(image, mask)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
module TestPromise

using Julio
using Test

function test_int()
    p = Julio.Promise{Int}()
    p[] = 111
    @test p[] == 111
end

function test_any()
    p = Julio.Promise()
    p[] = 111
    @test p[] == 111
end

function test_someany()
    p = Julio.Promise{Some{Any}}()
    p[] = Some{Any}(111)
    @test p[] === Some{Any}(111)
end

function test_union_someany_missing()
    p = Julio.Promise{Union{Some{Any},Missing}}()
    p[] = Some{Any}(111)
    @test p[] === Some{Any}(111)

    p = Julio.Promise{Union{Some{Any},Missing}}()
    p[] = missing
    @test p[] === missing
end

function test_maybefetch()
    @testset for (T, value) in [(Nothing, nothing), (Int, 111)]
        value::T
        test_maybefetch(T, value)
        test_maybefetch(Any, value)
    end
end

function test_maybefetch(T, value)
    p = Julio.Promise{T}()
    @test Julio.maybefetch(p) === nothing
    p[] = value
    @test Julio.maybefetch(p) === Some{T}(value)
end

end  # module


text:
module TimescaleDBs

using HTTP, JSON2
using Dates
using DataFrames

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils
import AutoMLPipeline.AbsTypes: fit!, transform!

export fit!,transform!
export TimescaleDB

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.
mutable struct TimescaleDB <: Transformer
  model
  args

  function TimescaleDB(args=Dict())
    default_args = Dict(
        :uri => "http://localhost:3000",
        :db => "dateval",
        :dbusername =>"postgres",
        :dbpassword => "mysecretpassword",
        :query => "select=date,value"
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(tdb::TimescaleDB, features::DataFrame=DataFrame(), labels::Vector=[])
  (features == DataFrame() && labels == [])  || error("features and labels should be empty because data are from http request")
  uri = tdb.args[:uri]; db = tdb.args[:db]; query = tdb.args[:query]
  (uri != "" && query != "" && db != "") || error("missing uri/query/db")
  tdb.model = tdb.args
end

function transform!(tdb::TimescaleDB, features::DataFrame=DataFrame())
  features == DataFrame()  || error("features should be empty because data are from http request")
  uri = tdb.args[:uri]; db = tdb.args[:db]; query = tdb.args[:query]
  (uri != "" && query != "" && db != "") || error("missing uri/query/db")
  payload = uri*"/"*db*"?"*query
  req = HTTP.get(payload)
  req.status == 200 || error(req.status)
  body = String(req.body)
  df = body |> JSON2.read |> DataFrame
  ncol(df) == 2 || error("data should have 2 columns")
  nrow(df) > 0 || return DataFrame()
  rename!(df,names(df)[1] => :Date, names(df)[2] => :Value)
  df.Date = DateTime.(df.Date)
  eltype(df.Value) <: Number || error("values in second column are not numeric")
  return df
end

end


text:
using Nuklear
using Nuklear.LibNuklear
using Nuklear.GLFWBackend
using GLFW, ModernGL
using Images, ImageMagick

const WINDOW_WIDTH = 800
const WINDOW_HEIGHT = 600

const MAX_VERTEX_BUFFER = 512 * 1024
const MAX_ELEMENT_BUFFER = 128 * 1024

@static if Sys.isapple()
    const VERSION_MAJOR = 3
    const VERSION_MINOR = 3
end

@static if Sys.isapple()
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, VERSION_MAJOR)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, VERSION_MINOR)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
else
    GLFW.DefaultWindowHints()
end

# set up GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
win = GLFW.CreateWindow(WINDOW_WIDTH, WINDOW_HEIGHT, "Demo")
@assert win != C_NULL
GLFW.MakeContextCurrent(win)
glViewport(0, 0, WINDOW_WIDTH, WINDOW_HEIGHT)

# init context
ctx = nk_glfw3_init(win, NK_GLFW3_INSTALL_CALLBACKS, MAX_VERTEX_BUFFER, MAX_ELEMENT_BUFFER)

# init font
nk_glfw3_font_stash_begin()
nk_glfw3_font_stash_end()

# creat texture for image drawing
img_width, img_height = 256, 256
tex_idx = nk_glfw3_create_texture(img_width, img_height)
image = rand(GLubyte, 4, img_width, img_height)

juliatan = load(joinpath(@__DIR__, "julia-tan.png")) |> transpose
juliatan_width, juliatan_height = size(juliatan)
juliatan_idx = nk_glfw3_create_texture(juliatan_width, juliatan_height,format=GL_RGBA, type=GL_FLOAT)
juliatan = channelview(juliatan)
juliatan = Array{Float32}(juliatan)

# states
const EASY = 0
const HARD = 1

op = EASY
property = Ref{Cint}(20)
bg = nk_colorf(0.10, 0.18, 0.24, 1.0)
while !GLFW.WindowShouldClose(win)
    global op; global bg; global ctx; global image; global juliatan;
    global tex_idx; global img_width; global img_height;

    GLFW.PollEvents()
    nk_glfw3_new_frame()

    if Bool(nk_begin(ctx, "Random Color", nk_rect(500, 0, 256, 256),
            NK_WINDOW_BORDER|NK_WINDOW_MOVABLE|NK_WINDOW_SCALABLE|
            NK_WINDOW_MINIMIZABLE|NK_WINDOW_TITLE))
        canvas = nk_window_get_canvas(ctx)
        region = nk_window_get_content_region(ctx)
        nk_glfw3_update_texture(tex_idx, image, img_width, img_height)
        img = Ref(create_nk_image(tex_idx))
        nk_draw_image(canvas, region, img, nk_rgba(255, 255, 255, 255))
    end
    nk_end(ctx)

    if Bool(nk_begin(ctx, "Julia-tan", nk_rect(0, 0, 500, 500),
            NK_WINDOW_BORDER|NK_WINDOW_MOVABLE|NK_WINDOW_SCALABLE|
            NK_WINDOW_MINIMIZABLE|NK_WINDOW_TITLE))
        canvas = nk_window_get_canvas(ctx)
        region = nk_window_get_content_region(ctx)
        nk_glfw3_update_texture(juliatan_idx, juliatan, juliatan_width, juliatan_height, format=GL_RGBA, type=GL_FLOAT)
        img = Ref(create_nk_image(juliatan_idx))
        nk_draw_image(canvas, region, img, nk_rgba(255, 255, 255, 255))
    end
    nk_end(ctx)

    # draw
    width, height = GLFW.GetWindowSize(win)
    glViewport(0, 0, width, height)
    glClear(GL_COLOR_BUFFER_BIT)
    glClearColor(bg.r, bg.g, bg.b, bg.a)
    nk_glfw3_render(NK_ANTI_ALIASING_ON, MAX_VERTEX_BUFFER, MAX_ELEMENT_BUFFER)
    GLFW.SwapBuffers(win)
end

nk_glfw3_shutdown()
GLFW.DestroyWindow(win)


text:
using StaticArrays
using BenchmarkTools

using RecursiveArrayTools

N = 5
K = 64*64

a,b = 1.0,1.0
Np = (N+1)*(N+2)÷2
Q = SVector{4}([randn(Np,K) for i = 1:4])
rhsQ = SVector{4}([randn(Np,K) for i = 1:4])
function inplace1!(Q,rhsQ,a,b)
    for i = 1:length(Q)
        @. Q[i] = a*Q[i] + b*rhsQ[i]
    end
end
inplace2!(Q,rhsQ,a,b) = Q = @. a*Q + b*rhsQ
inplace3!(Q,rhsQ,a,b) = @. Q = a*Q + b*rhsQ

# bcopy!(x,y) = x.=y

Q2 = VectorOfArray(Q)
rhsQ2 = VectorOfArray(rhsQ)

@btime inplace1!($Q,$rhsQ,$a,$b)   # 208.795 μs (0 allocations: 0 bytes)
@btime inplace2!.($Q,$rhsQ,$a,$b)  # 263.483 μs (8 allocations: 2.63 MiB)
@btime inplace3!.($Q,$rhsQ,$a,$b)  # 216.073 μs (0 allocations: 0 bytes)
@btime inplace3!($Q2,$rhsQ2,$a,$b) # 217.012 μs (0 allocations: 0 bytes)
@btime axpby!.($a, $Q, $b, $rhsQ)  # 113.659 μs (0 allocations: 0 bytes)


text:
using DiffusionMoments
using Test

@testset "DiffusionMoments.jl" begin
    # Write your tests here.
end

max_stars_repo_name: SosUts/DiffusionMoments

max_issues_repo_name: SosUts/DiffusionMoments

max_forks_repo_name: SosUts/DiffusionMoments


text:
module Juqbox

using LinearAlgebra
using Plots
#pyplot()
using Printf
using Random
using LaTeXStrings
using SparseArrays
using Ipopt
using FileIO
using FFTW
using DelimitedFiles

export splineparams, bspline2, gradbspline2
export bcparams, bcarrier2, gradbcarrier2!

export objparams, traceobjgrad, identify_guard_levels, identify_forbidden_levels, initial_cond
export plotunitary, plotspecified, evalctrl, plot_results, plot_energy, plot_final_unitary
export setup_ipopt_problem, Working_Arrays, estimate_Neumann!, setup_rotmatrices
export run_optimizer, plot_conv_hist, wmatsetup
export zero_start_end!, assign_thresholds, assign_thresholds_freq, assign_thresholds_ctrl_freq
export calculate_timestep, marginalize3
export save_pcof, read_pcof, juq2qis

# Julia versions prior to v"1.3.1" can't use LinearAlgebra's 5 argument mul!, routines
# included here for backwards compatability
if(VERSION < v"1.3.1")
    include("backwards_compat.jl")
end

include("bsplines.jl") # add all B-spline functionality to the Juqbox module

include("StormerVerlet.jl") # add in time-stepping functionality

# union type for different ctrl parameterizations
BsplineParams = Union{splineparams, bcparams}
# union type for different matrix representations
MyRealMatrix = Union{Array{Float64,2}, SparseMatrixCSC{Float64, Int64}}

include("evalobjgrad.jl")

include("plotstatectrl.jl")

include("plot-results.jl")

include("ipopt_interface.jl")

include("save_pcof.jl")

end # module


text:
using DFTK
using Plots

# This example plots the free-electron bands for a 1D system.
# TODO Convert into an example in the docs
#      show how one could add a potential to this and see the bands change.

lattice = diagm([5., 0, 0])
model   = Model(lattice; n_electrons=4, terms=[Kinetic()])
basis   = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1));

n_bands = 6
ρ0 = guess_density(basis)  # Just dummy, has no meaning in this model
p = plot_bandstructure(basis, ρ0, n_bands, kline_density=15)


text:
# websafe colors
# using Colors, ColorSchemes
# cs = ColorScheme([parse(RGB{Float64}, "#$(string.([r,g,b], base=16)...)") for r in 0x0:3:0xf for g in 0x0:3:0xf for b in 0x0:3:0xf])

loadcolorscheme(:websafe, [
        RGB{Float64}(0.0,0.0,0.0),
        RGB{Float64}(0.0,0.0,0.2),
        RGB{Float64}(0.0,0.0,0.4),
        RGB{Float64}(0.0,0.0,0.6),
        RGB{Float64}(0.0,0.0,0.8),
        RGB{Float64}(0.0,0.0,1.0),
        RGB{Float64}(0.0,0.2,0.0),
        RGB{Float64}(0.0,0.2,0.2),
        RGB{Float64}(0.0,0.2,0.4),
        RGB{Float64}(0.0,0.2,0.6),
        RGB{Float64}(0.0,0.2,0.8),
        RGB{Float64}(0.0,0.2,1.0),
        RGB{Float64}(0.0,0.4,0.0),
        RGB{Float64}(0.0,0.4,0.2),
        RGB{Float64}(0.0,0.4,0.4),
        RGB{Float64}(0.0,0.4,0.6),
        RGB{Float64}(0.0,0.4,0.8),
        RGB{Float64}(0.0,0.4,1.0),
        RGB{Float64}(0.0,0.6,0.0),
        RGB{Float64}(0.0,0.6,0.2),
        RGB{Float64}(0.0,0.6,0.4),
        RGB{Float64}(0.0,0.6,0.6),
        RGB{Float64}(0.0,0.6,0.8),
        RGB{Float64}(0.0,0.6,1.0),
        RGB{Float64}(0.0,0.8,0.0),
        RGB{Float64}(0.0,0.8,0.2),
        RGB{Float64}(0.0,0.8,0.4),
        RGB{Float64}(0.0,0.8,0.6),
        RGB{Float64}(0.0,0.8,0.8),
        RGB{Float64}(0.0,0.8,1.0),
        RGB{Float64}(0.0,1.0,0.0),
        RGB{Float64}(0.0,1.0,0.2),
        RGB{Float64}(0.0,1.0,0.4),
        RGB{Float64}(0.0,1.0,0.6),
        RGB{Float64}(0.0,1.0,0.8),
        RGB{Float64}(0.0,1.0,1.0),
        RGB{Float64}(0.2,0.0,0.0),
        RGB{Float64}(0.2,0.0,0.2),
        RGB{Float64}(0.2,0.0,0.4),
        RGB{Float64}(0.2,0.0,0.6),
        RGB{Float64}(0.2,0.0,0.8),
        RGB{Float64}(0.2,0.0,1.0),
        RGB{Float64}(0.2,0.2,0.0),
        RGB{Float64}(0.2,0.2,0.2),
        RGB{Float64}(0.2,0.2,0.4),
        RGB{Float64}(0.2,0.2,0.6),
        RGB{Float64}(0.2,0.2,0.8),
        RGB{Float64}(0.2,0.2,1.0),
        RGB{Float64}(0.2,0.4,0.0),
        RGB{Float64}(0.2,0.4,0.2),
        RGB{Float64}(0.2,0.4,0.4),
        RGB{Float64}(0.2,0.4,0.6),
        RGB{Float64}(0.2,0.4,0.8),
        RGB{Float64}(0.2,0.4,1.0),
        RGB{Float64}(0.2,0.6,0.0),
        RGB{Float64}(0.2,0.6,0.2),
        RGB{Float64}(0.2,0.6,0.4),
        RGB{Float64}(0.2,0.6,0.6),
        RGB{Float64}(0.2,0.6,0.8),
        RGB{Float64}(0.2,0.6,1.0),
        RGB{Float64}(0.2,0.8,0.0),
        RGB{Float64}(0.2,0.8,0.2),
        RGB{Float64}(0.2,0.8,0.4),
        RGB{Float64}(0.2,0.8,0.6),
        RGB{Float64}(0.2,0.8,0.8),
        RGB{Float64}(0.2,0.8,1.0),
        RGB{Float64}(0.2,1.0,0.0),
        RGB{Float64}(0.2,1.0,0.2),
        RGB{Float64}(0.2,1.0,0.4),
        RGB{Float64}(0.2,1.0,0.6),
        RGB{Float64}(0.2,1.0,0.8),
        RGB{Float64}(0.2,1.0,1.0),
        RGB{Float64}(0.4,0.0,0.0),
        RGB{Float64}(0.4,0.0,0.2),
        RGB{Float64}(0.4,0.0,0.4),
        RGB{Float64}(0.4,0.0,0.6),
        RGB{Float64}(0.4,0.0,0.8),
        RGB{Float64}(0.4,0.0,1.0),
        RGB{Float64}(0.4,0.2,0.0),
        RGB{Float64}(0.4,0.2,0.2),
        RGB{Float64}(0.4,0.2,0.4),
        RGB{Float64}(0.4,0.2,0.6),
        RGB{Float64}(0.4,0.2,0.8),
        RGB{Float64}(0.4,0.2,1.0),
        RGB{Float64}(0.4,0.4,0.0),
        RGB{Float64}(0.4,0.4,0.2),
        RGB{Float64}(0.4,0.4,0.4),
        RGB{Float64}(0.4,0.4,0.6),
        RGB{Float64}(0.4,0.4,0.8),
        RGB{Float64}(0.4,0.4,1.0),
        RGB{Float64}(0.4,0.6,0.0),
        RGB{Float64}(0.4,0.6,0.2),
        RGB{Float64}(0.4,0.6,0.4),
        RGB{Float64}(0.4,0.6,0.6),
        RGB{Float64}(0.4,0.6,0.8),
        RGB{Float64}(0.4,0.6,1.0),
        RGB{Float64}(0.4,0.8,0.0),
        RGB{Float64}(0.4,0.8,0.2),
        RGB{Float64}(0.4,0.8,0.4),
        RGB{Float64}(0.4,0.8,0.6),
        RGB{Float64}(0.4,0.8,0.8),
        RGB{Float64}(0.4,0.8,1.0),
        RGB{Float64}(0.4,1.0,0.0),
        RGB{Float64}(0.4,1.0,0.2),
        RGB{Float64}(0.4,1.0,0.4),
        RGB{Float64}(0.4,1.0,0.6),
        RGB{Float64}(0.4,1.0,0.8),
        RGB{Float64}(0.4,1.0,1.0),
        RGB{Float64}(0.6,0.0,0.0),
        RGB{Float64}(0.6,0.0,0.2),
        RGB{Float64}(0.6,0.0,0.4),
        RGB{Float64}(0.6,0.0,0.6),
        RGB{Float64}(0.6,0.0,0.8),
        RGB{Float64}(0.6,0.0,1.0),
        RGB{Float64}(0.6,0.2,0.0),
        RGB{Float64}(0.6,0.2,0.2),
        RGB{Float64}(0.6,0.2,0.4),
        RGB{Float64}(0.6,0.2,0.6),
        RGB{Float64}(0.6,0.2,0.8),
        RGB{Float64}(0.6,0.2,1.0),
        RGB{Float64}(0.6,0.4,0.0),
        RGB{Float64}(0.6,0.4,0.2),
        RGB{Float64}(0.6,0.4,0.4),
        RGB{Float64}(0.6,0.4,0.6),
        RGB{Float64}(0.6,0.4,0.8),
        RGB{Float64}(0.6,0.4,1.0),
        RGB{Float64}(0.6,0.6,0.0),
        RGB{Float64}(0.6,0.6,0.2),
        RGB{Float64}(0.6,0.6,0.4),
        RGB{Float64}(0.6,0.6,0.6),
        RGB{Float64}(0.6,0.6,0.8),
        RGB{Float64}(0.6,0.6,1.0),
        RGB{Float64}(0.6,0.8,0.0),
        RGB{Float64}(0.6,0.8,0.2),
        RGB{Float64}(0.6,0.8,0.4),
        RGB{Float64}(0.6,0.8,0.6),
        RGB{Float64}(0.6,0.8,0.8),
        RGB{Float64}(0.6,0.8,1.0),
        RGB{Float64}(0.6,1.0,0.0),
        RGB{Float64}(0.6,1.0,0.2),
        RGB{Float64}(0.6,1.0,0.4),
        RGB{Float64}(0.6,1.0,0.6),
        RGB{Float64}(0.6,1.0,0.8),
        RGB{Float64}(0.6,1.0,1.0),
        RGB{Float64}(0.8,0.0,0.0),
        RGB{Float64}(0.8,0.0,0.2),
        RGB{Float64}(0.8,0.0,0.4),
        RGB{Float64}(0.8,0.0,0.6),
        RGB{Float64}(0.8,0.0,0.8),
        RGB{Float64}(0.8,0.0,1.0),
        RGB{Float64}(0.8,0.2,0.0),
        RGB{Float64}(0.8,0.2,0.2),
        RGB{Float64}(0.8,0.2,0.4),
        RGB{Float64}(0.8,0.2,0.6),
        RGB{Float64}(0.8,0.2,0.8),
        RGB{Float64}(0.8,0.2,1.0),
        RGB{Float64}(0.8,0.4,0.0),
        RGB{Float64}(0.8,0.4,0.2),
        RGB{Float64}(0.8,0.4,0.4),
        RGB{Float64}(0.8,0.4,0.6),
        RGB{Float64}(0.8,0.4,0.8),
        RGB{Float64}(0.8,0.4,1.0),
        RGB{Float64}(0.8,0.6,0.0),
        RGB{Float64}(0.8,0.6,0.2),
        RGB{Float64}(0.8,0.6,0.4),
        RGB{Float64}(0.8,0.6,0.6),
        RGB{Float64}(0.8,0.6,0.8),
        RGB{Float64}(0.8,0.6,1.0),
        RGB{Float64}(0.8,0.8,0.0),
        RGB{Float64}(0.8,0.8,0.2),
        RGB{Float64}(0.8,0.8,0.4),
        RGB{Float64}(0.8,0.8,0.6),
        RGB{Float64}(0.8,0.8,0.8),
        RGB{Float64}(0.8,0.8,1.0),
        RGB{Float64}(0.8,1.0,0.0),
        RGB{Float64}(0.8,1.0,0.2),
        RGB{Float64}(0.8,1.0,0.4),
        RGB{Float64}(0.8,1.0,0.6),
        RGB{Float64}(0.8,1.0,0.8),
        RGB{Float64}(0.8,1.0,1.0),
        RGB{Float64}(1.0,0.0,0.0),
        RGB{Float64}(1.0,0.0,0.2),
        RGB{Float64}(1.0,0.0,0.4),
        RGB{Float64}(1.0,0.0,0.6),
        RGB{Float64}(1.0,0.0,0.8),
        RGB{Float64}(1.0,0.0,1.0),
        RGB{Float64}(1.0,0.2,0.0),
        RGB{Float64}(1.0,0.2,0.2),
        RGB{Float64}(1.0,0.2,0.4),
        RGB{Float64}(1.0,0.2,0.6),
        RGB{Float64}(1.0,0.2,0.8),
        RGB{Float64}(1.0,0.2,1.0),
        RGB{Float64}(1.0,0.4,0.0),
        RGB{Float64}(1.0,0.4,0.2),
        RGB{Float64}(1.0,0.4,0.4),
        RGB{Float64}(1.0,0.4,0.6),
        RGB{Float64}(1.0,0.4,0.8),
        RGB{Float64}(1.0,0.4,1.0),
        RGB{Float64}(1.0,0.6,0.0),
        RGB{Float64}(1.0,0.6,0.2),
        RGB{Float64}(1.0,0.6,0.4),
        RGB{Float64}(1.0,0.6,0.6),
        RGB{Float64}(1.0,0.6,0.8),
        RGB{Float64}(1.0,0.6,1.0),
        RGB{Float64}(1.0,0.8,0.0),
        RGB{Float64}(1.0,0.8,0.2),
        RGB{Float64}(1.0,0.8,0.4),
        RGB{Float64}(1.0,0.8,0.6),
        RGB{Float64}(1.0,0.8,0.8),
        RGB{Float64}(1.0,0.8,1.0),
        RGB{Float64}(1.0,1.0,0.0),
        RGB{Float64}(1.0,1.0,0.2),
        RGB{Float64}(1.0,1.0,0.4),
        RGB{Float64}(1.0,1.0,0.6),
        RGB{Float64}(1.0,1.0,0.8),
        RGB{Float64}(1.0,1.0,1.0),
    ], "general", "216 websafe colors")


text:
import DemeAssemblies
using DemeAssemblies: SystemConfig, datadir, serve

using PeaceVote
setnamespace(@__MODULE__)
uuid = PeaceVote.uuid("DemeAssemblies")
server = PeaceVote.Signer(uuid,"server") ### One needs to take care of accounts

MAINTAINER_ID = PeaceVote.Signer(uuid,"maintainer").id

deme = MAINTAINER_ID
config = SystemConfig(deme)
system = System(config,deme,server)

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
"""
    modular_exponentiation(base::Int64, exponent::Int64, modulus::Int64)

Compute the `residuum` of the base raised to the `exponent`, which is divided by the
`modulus`.

# Arguments
- `base::Int64`: base
- `exponent::Int64`: exponent
- `modulus ::Int64`: modulus

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> ClassicAlgorithmsCollections.modular_exponentiation(2, 3, 15)
8
```
"""
function modular_exponentiation(base::Int64, exponent::Int64, modulus::Int64)
    residuum = 1     # Initialize result

    # Update base if it is more
    # than or equal to modulus
    base = base % modulus

    if (base == 0)
        return 0
    end

    while (exponent > 0)

        # If exponent is odd, multiplexponent
        # base with result
        if ((exponent & 1) == 1)
            residuum = (residuum * base) % modulus
        end
        # exponent must be even now
        exponent = exponent >> 1      # exponent = exponent/2
        base = (base * base) % modulus
    end

    return residuum
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
#--- test linear solve method (Foster PhD)
using Plots, Revise, VortexDistributions

pv1 = PointVortex(0.1,0.3,1)
v1 = ScalarVortex(pv1)

x = LinRange(-10,10,256)
y = x

psi = Torus(one.(x.*y') |> complex,x,y)

vortex!(psi,v1)

psiv = psi.ψ

heatmap(x,y,angle.(psiv))

#--- define linear solve (Foster)

function vortex_detect2d(psi)
# Detect vortices in a field psi, on plaquettes of a 2D grid.
#
# First we localize vortices on the grid plaquettes using the phase winding
# algorithm.  After we've found which grid plaquettes contain vortices, we
# estimate the positions more accurately by solving the equations
#
#   real(psi(x,y)) = 0, imag(psi(x,y)) = 0
#
# In principle we could do this using the exact psi as computed via spectral
# interpolation (provided psi is defined in a spectral basis) but that would
# probably be relatively slow.  Instead, we make a linear approximation over
# the plaquette for the real and imaginary parts of psi, and find the position
# where these approximations are mutually zero.  This corresponds to the
# position of the vortex.
#
# The linear approximations are computed from the values of psi at the corners
# of the plaquette using least squares.
#
# Inputs:
#   psi - complex wavefunction in which to detect vortices.
#
# Outputs:
#   positions - Array of vortex positions with shape 2xN where N is the number
#               of vortices:  [ x1  x2  x3 ... xN ]
#                             [ y1  y2  y3 ... yN ]
#   winding - The vortex winding number of the positions array, +1 or -1.

vorticity = phase_winding2d(angle(psi));
[iy,ix] = find(vorticity);
positions = [ix,iy]';

if nargout > 1
    winding = vorticity(vorticity~=0);
end

# % Compute least squares estimator for the linear function coefficients on a
# % plaquette.  Suppose the real part of psi on a plaquette containing a vortex
# % is
# %
# %   [r1 r2;
# %    r3 r4]
# %
# % We want to fit a plane (a*x + b*y + c) to this such that
# %
# %   r1 = a*0 + b*0 + c*1
# %   r2 = a*1 + b*0 + c*1
# %   r3 = a*0 + b*1 + c*1
# %   r4 = a*1 + b*1 + c*1
# %
# % In matrix form, that would be L * [a;b;c] = [r1;r2;r3;r4], which is the
# % origin of the mysterious looking matrix below.
# %
L = [0 0 1;
     1 0 1;
     0 1 1;
     1 1 1];
lsEst = inv(L'*L)*L';
# % (For future reference, lsEst actually turns out to be the following)
# % lsEst = 1/4 * [-2,  2, -2,  2;
# %                -2, -2,  2,  2;
# %                 3,  1,  1, -1];
for i = 1:length(ix)
    # Extract relevant grid plaquette containing the vortex.
    p = psi(iy(i)+[0,1], ix(i)+[0,1]).'
    # Compute linear approximations of real & imag parts of p
    rcoeff = lsEst*real(p[:])
    icoeff = lsEst*imag(p[:])
    # Find position where these are mutually zero
    pos = [rcoeff(1:2), icoeff(1:2)]' \ -[rcoeff(3); icoeff(3)]
    positions[:,i] = positions[:,i] + pos;
    return positions, winding
end


text:
# adjoint-rotate.jl
# test adjoint consistency for rotate methods on very small case

using SPECTrecon: plan_rotate
using SPECTrecon: imrotate, imrotate_adj
using SPECTrecon: imrotate!, imrotate_adj!
using LinearMapsAA: LinearMapAA
using Test: @test, @testset

@testset "rotate" begin
    nx = 10
    plan = plan_rotate(nx)
    show(isinteractive() ? stdout : devnull, "text/plain", plan)
    show(isinteractive() ? stdout : devnull, "text/plain", plan[1])
    @test sizeof(plan) isa Int
    @test sizeof(plan[1]) isa Int
end

@testset "rotate3" begin
    for method in (:one, :two)
        nx = 20
        θ = π/6
        T = Float32
        plan = plan_rotate(nx; T, method)
        image3 = rand(T, nx, nx, 4)
        result1 = similar(image3)
        result2 = similar(result1)
        imrotate!(result1, image3, θ, plan)
        @test maximum(result1) ≤ 1.01
        imrotate!(result2, image3, θ, plan, :thread)
        @test result1 == result2

        imrotate_adj!(result1, image3, θ, plan)
        @test maximum(result1) ≤ 1.5
        imrotate_adj!(result2, image3, θ, plan, :thread)
        @test result1 == result2
    end
end

@testset "adjoint-rotate" begin
    Ntest = 7
    θ_list = (0:Ntest-1) / Ntest * 2π
    M = 16
    N = M
    T = Float32

    for method in (:one, :two), θ in θ_list
        idim = (M, N)
        odim = (M, N)
        x = randn(T, M, N)
        forw = x -> imrotate(x, θ; method)
        back = x -> imrotate_adj(x, θ; method)
        A = LinearMapAA(forw, back, (prod(odim), prod(idim)); T, odim, idim)
        @test Matrix(A') ≈ Matrix(A)'
    end
end


text:
module ML

using JuliaDB
using Dagger
using OnlineStats
using PooledArrays
using DataValues

import Dagger: ArrayOp, DArray, treereduce
import JuliaDB: Dataset, DDataset, nrows

import Base: merge

# Core schema types

schema(xs::AbstractArray) = nothing # catch-all
schema(xs::AbstractArray, ::Void) = nothing # catch-all
merge(::Void, ::Void) = nothing
width(::Void) = 0
featuremat!(A, ::Void, xs) = A

schema(xs::ArrayOp) = schema(compute(xs))
schema(xs::ArrayOp, T) = schema(compute(xs), T)

struct Continuous
  series::Series
end

function schema(xs, ::Type{Continuous})
  Continuous(fit!(Series(Variance()),xs))
end
function schema(xs::AbstractArray{<:Real})
    schema(xs, Continuous)
end
width(::Continuous) = 1
function merge(c1::Continuous, c2::Continuous)
    Continuous(merge(c1.series, c2.series))
end
Base.mean(c::Continuous)::Float64 = mean(c.series.stats[1])
Base.std(c::Continuous)::Float64 = std(c.series.stats[1])
function Base.show(io::IO, c::Continuous)
    write(io, "Continous(μ=$(mean(c)), σ=$(std(c)))")
end
Base.@propagate_inbounds function featuremat!(A, c::Continuous, xs,
                                         dropna=Val{false}())
    m = mean(c)
    s = std(c)
    for i in 1:length(xs)
        x =  xs[i]
        if dropna isa Val{true}
            if isnull(x)
                continue
            else
                A[i, 1] = (get(x) - m) / s
            end
        else
            A[i, 1] = (x - m) / s
        end
    end
    A
end

struct Categorical
  series::Series
end
function Categorical(xs::AbstractArray)
    Categorical(Series(CountMap(Dict(zip(xs, zeros(Int, length(xs)))))))
end

function schema(xs::AbstractArray, ::Type{Categorical})
  Categorical(fit!(Series(CountMap(eltype(xs))),xs))
end

function schema(xs::PooledArray)
    schema(xs, Categorical)
end

Base.keys(c::Categorical) = keys(c.series.stats[1])
width(c::Categorical) = length(keys(c))
merge(c1::Categorical, c2::Categorical) = Categorical(merge(c1.series, c2.series))
function Base.show(io::IO, c::Categorical)
    write(io, "Categorical($(collect(keys(c))))")
end
Base.@propagate_inbounds function featuremat!(A, c::Categorical, xs, dropna=Val{false}())
    ks = keys(c)
    labeldict = Dict{eltype(ks), Int}(zip(ks, 1:length(ks)))
    for i = 1:length(xs)
        if dropna isa Val{true} && isnull(xs[i])
            continue
        end
        A[i, labeldict[xs[i]]] = one(eltype(A))
    end
    A
end

# distributed schema calculation
function schema(xs::DArray)
    collect(treereduce(delayed(merge),
                       delayedmap(x -> schema(x), xs.chunks)))
end

for S in (Continuous, Categorical)
    @eval function schema(xs::DArray, T::Type{$S})
        collect(treereduce(delayed(merge),
                           delayedmap(x -> schema(x, T), xs.chunks)))
    end
end

struct Maybe{T}
  feature::T
end
function schema(xs::DataValueArray, T::Type)
    Maybe(schema(dropna(xs), T))
end
schema(xs::DataValueArray) = Maybe(schema(dropna(xs)))
width(c::Maybe) = width(c.feature) + 1
merge(m1::Maybe, m2::Maybe) = Maybe(merge(m1.feature, m2.feature))
nulls(xs) = Base.Generator(isnull, xs)
nulls(xs::DataValueArray) = xs.isnull
Base.@propagate_inbounds function featuremat!(A, c::Maybe, xs, dropna=Val{true}())
    copy!(A, CartesianRange((1:length(xs), 1:1)), reshape(nulls(xs), (length(xs), 1)), CartesianRange((1:length(xs), 1:1)))
    featuremat!(view(A, 1:length(xs), 2:size(A, 2)), c.feature, xs, Val{true}())
    A
end

# Schema inference

const Schema = Dict{Symbol,Any}

# vecTs: type of column vectors in each chunk
function schema(cols, names; hints=Dict())
    d = Schema()
    for (col, name) in zip(cols, names)
        if haskey(hints, name)
            d[name] = schema(col, hints[name])
        else
            d[name] = schema(col)
        end
    end
    d
end

function schema(t::Union{Dataset, DDataset}; hints=Dict())
    schema(collect(columns(t)), colnames(t), hints=hints)
end

width(sch::Schema) = sum(width(s) for s in values(sch))
function featuremat!(A, schemas::Schema, t::Dataset)
    j = 0
    for col in keys(schemas)
        schema = schemas[col]
        featuremat!(view(A, 1:length(t), j+1:j+width(schema)),
                   schema, column(t, col))

        j += width(schema)
    end
    A
end

splitschema(xs::Schema, ks...) =
    filter((k,v) -> k ∉ ks, xs),
    filter((k,v) -> k ∈ ks, xs)

function featuremat(sch, xs)
    featuremat!(zeros(Float32, length(xs), width(sch)), sch, xs)'
end
featuremat(t) = featuremat(schema(t), t)

function featuremat(s, t::DDataset)
    t = compute(t)
    w = width(s)
    h = length(t)
    lengths = get.(nrows.(t.domains))
    domains = Dagger.DomainBlocks((1,1), ([w], cumsum(lengths)))

    DArray(Float32,
           Dagger.ArrayDomain(w, sum(lengths)),
           domains,
           reshape(delayedmap(x->featuremat(s,x), t.chunks),
                   (1, length(t.chunks))))
end

end


text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

# create pde structure and mesh structure
pde, mesh = Preprocessing.initializeexasim();

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde.model = "ModelD";            # ModelC, ModelD, ModelW
include("pdemodel4.jl");          # include the PDE model file

# Set discretization parameters, physical parameters, and solver parameters
pde.porder = 4;          # polynomial degree
pde.torder = 2;          # time-stepping order of accuracy
pde.nstage = 2;          # time-stepping number of stages
pde.dt = 0.002*ones(5000);   # time step sizes
pde.soltime = collect(50:50:length(pde.dt)); # steps at which solution are collected
pde.visdt = 1; # visualization timestep size

gam = 1.0;                      # gravity
pde.physicsparam = [gam 0.0];
pde.tau = [1.0];          # DG stabilization parameter
pde.GMRESrestart=15;            # number of GMRES restarts
pde.linearsolvertol=1e-12;      # GMRES tolerance
pde.linearsolveriter=16;        # number of GMRES iterations
pde.precMatrixType=2;           # preconditioning type
pde.NLtol = 1e-12;              # Newton tolerance
pde.NLiter=2;                   # Newton iterations

# parameters for differential algebraic equations
pde.dae_alpha = 1e-4;
pde.dae_beta = 0.0;

# Choose computing platform and set number of processors
#pde.platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde.mpiprocs = 1;                # number of MPI processors

# create a linear mesh for a square domain
mesh.p, mesh.t = Mesh.SquareMesh(64,64,1); # a mesh of 8 by 8 quadrilaterals
mesh.p = (4*pi)*mesh.p .- 2*pi;
# expressions for disjoint boundaries
mesh.boundaryexpr = [p -> (p[2,:] .< -2*pi+1e-3), p -> (p[1,:] .> 2*pi-1e-3), p -> (p[2,:] .> 2*pi-1e-3), p -> (p[1,:] .< -2*pi+1e-3)];
mesh.boundarycondition = [1 1 1 1]; # Set boundary condition for each disjoint boundary
mesh.periodicexpr = [2 p->p[2,:] 4 p->p[2,:]; 1 p->p[1,:] 3 p->p[1,:]];

# call exasim to generate and run C++ code to solve the PDE model
sol, pde, mesh,~,~,~,~  = Postprocessing.exasim(pde,mesh);

# # search compilers and set options
# pde = Gencode.setcompilers(pde);

# # generate input files and store them in datain folder
# pde, mesh, master, dmd = Preprocessing.preprocessing(pde,mesh);

# # generate source codes and store them in app folder
# Gencode.gencode(pde);

# # compile source codes to build an executable file and store it in app folder
# compilerstr = Gencode.compilecode(pde);

# # run executable file to compute solution and store it in dataout folder
# runstr = Gencode.runcode(pde, 1);

# # get solution from output files in dataout folder
# sol = Postprocessing.fetchsolution(pde,master,dmd);

# visualize the numerical solution of the PDE model using Paraview
pde.visscalars = ["uy", 4, "vx", 5];  # list of scalar fields for visualization
pde.visvectors = ["velocity", [1, 2]]; # list of vector fields for visualization
Postprocessing.vis(sol,pde,mesh); # visualize the numerical solution
print("Done!");


text:
"""
SamplerState for relativistic Hamiltonian Monte Carlo

X. Lu, V. Perrone, L. Hasenclever, Y. W. Teh, Sebastian J. Vollmer relativistic Monte Carlo https://arxiv.org/abs/1609.04388

Fields:
    - x: the current state
    - p: the current momentum
    - stepsize
    - niters: number of iterations
    - c: speed of light
    - mass
"""
type RelHMCState <: SamplerState
    x::Array{Float64}
    p::Array{Float64} # included in sampler state to allow autoregressive updates.

    stepsize::Float64 #stepsize
    niters::Int64#L
    c::Float64
    mass
    function RelHMCState(x::Array{Float64},stepsize::Float64;niters=10,c = 1.0,mass=1.0,p=:none)
        if isa(mass,Number)
          mass = mass * ones(length(x))
        end
        if p == :none
            p = sample_rel_p(mass, c, length(x))
        end
        new(x,p,stepsize,niters,c,mass)
    end
end

"""
sample! performs one relativistic HMC update.
"""
function sample!(s::RelHMCState, llik, grad)
  nparams = length(s.x)
  mass = s.mass
  stepsize = s.stepsize
  niters = s.niters
  c = s.c

  #resample relativistic momentum
  s.p = sample_rel_p(mass, c, nparams)

 #pp = sqrt(mass).*randn(nparams)
  curx = s.x
  curp = s.p
  s.p += .5*stepsize * grad(s.x)
  for iter = 1:niters
    s.x += stepsize * s.p./(mass .*sqrt(s.p.^2 ./ (mass.*c).^2 + 1))
    s.p += (iter<niters ? stepsize : .5*stepsize) * grad(s.x) # two leapfrog steps rolled in one unless at the end.
  end

  #current kinetic energy
  cur_ke = sum(mass.*c.^2 .* sqrt(curp.^2 ./(mass.*c).^2 +1))[1]
  ke = sum(mass.*c.^2 .* sqrt(s.p.^2 ./(mass.*c).^2 +1))[1]

  logaccratio = llik(s.x) - llik(curx) - ke + cur_ke
  if 0.0 > logaccratio - log(rand())
      #reject
      s.x = curx
      s.p = curp
  else
      #accept
      #negate momentum for symmetric Metropolis proposal
      s.p = -s.p
  end
  return s
end

# logpdf of the relhmc momentum variable
function gen_rel_p_logpdf(m, c)
    function rel_p_logpdf(p)
        -m*c^2*sqrt(p^2/(m^2*c^2)+1)
    end

    rel_p_logpdf
end

function sample_rel_p(mass, c, nparams; bounds=[-Inf, Inf])

  if length(mass) == 1 && length(c) == 1
    p_logpdf = gen_rel_p_logpdf(mass[1], c[1])
    pp = ars(p_logpdf, -10.0, 10.0, bounds, nparams)
  else
    mass = length(mass) == 1 ? mass * ones(nparams) : mass
    c = length(c) == 1 ? c .* ones(nparams) : c

    pp = zeros(nparams)

    for i = 1:nparams
      p_logpdf = gen_rel_p_logpdf(mass[i],c[i])
      pp[i] = ars(p_logpdf, -10.0, 10.0, bounds, 1)[1]
    end

  end

  pp
end


text:
using DataFrames
using CSV
using Plots
#% %
file1 = open("rand_decay.dat")
file2 = open("exp_decay.dat")
df1=CSV.read(file1, header=false, delim=" ")
df2=CSV.read(file2, header=false, delim=" ")
#%%
plot(df2[1],df2[2], title="radioactive decay")
scatter!(df1[1],df1[2], markersize=2)
xaxis!("t", xlims=(1,Inf))
yaxis!("N", ylims=(1,Inf), yscale= :log10)

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
function build_lambda(y_size::Int64)
    λ = Array{Float64}(undef, y_size, y_size)
    if y_size == 30
        λ = [-0.83533612 0.014477463 0.79244709 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0015233196 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0011365784 -0.83341414 0.0 0.011937726 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0015233196 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.024511464 0.0 -0.067400493 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.014477463 0.024511464 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0;
             0.0011365784 0.0015233196 0.0 -0.83312459 0.010219256 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0015233196 0.0 0.0 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 -0.83284759 0.009430315 6.8793952e-6 9.4399775e-5 0.000518823 0.0015233196 0.0 0.0 0.0 0.79244709 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 9.4399749e-5 -0.82820667 0.0047018731 9.4399775e-5 0.000518823 0.0015233196 0.0 0.0 0.0 0.0 0.79244709 0.0 0.0 0.0 0.0011365783 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 9.4399749e-5 0.004701886 -0.82820668 9.4399775e-5 0.000518823 0.0015233196 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0 0.0 0.0011365783 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 0.0094303151 -0.83284759 0.000518823 0.0015233196 0.0 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0 0.0011365783 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 0.010219256 -0.83312459 0.0015233196 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0011365783 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.011937727 -0.83341414 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0;
             0.0 0.024511464 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.065478512 0.0 0.011937726 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0;
             0.024511464 0.0 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.83533612 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.014477463 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0011365784 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0015233196 0.0 -0.065188961 0.010219256 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0;
             0.0 0.0 0.0011365784 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0015233196 0.0 0.00051882295 -0.064911964 0.009430315 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0011365783 0.0 0.0;
             0.0 0.0 0.0011365784 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0015233196 0.0 0.00051882295 9.4399749e-5 -0.060271042 0.0047018731 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0011365783 0.0 0.0;
             0.0 0.0 0.0011365784 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0015233196 0.0 0.00051882295 9.4399749e-5 0.004701886 -0.060271055 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0011365783 0.0 0.0;
             0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 0.0094303151 -0.064911964 0.000518823 0.0 0.0015233196 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0011365783 0.0 0.0;
             0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 0.010219256 -0.065188961 0.0 0.0015233196 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0011365783 0.0 0.0;
             0.0011365784 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.014477463 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.83533612 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0 0.024511464;
             0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.011937727 0.0 -0.065478512 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365783 0.024511464 0.0;
             0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.83341414 0.011937726 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365784 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0015233196 -0.83312459 0.010219256 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365784 0.0 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0 0.0015233196 0.00051882295 -0.83284759 0.009430315 6.8793952e-6 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0 0.0011365784 0.0 0.0 0.79244709 0.0 0.0 0.0 0.0 0.0 0.0015233196 0.00051882295 9.4399749e-5 -0.82820667 0.0047018731 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0 0.0011365784 0.0 0.0 0.0 0.79244709 0.0 0.0 0.0 0.0 0.0015233196 0.00051882295 9.4399749e-5 0.004701886 -0.82820668 9.4399775e-5 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.79244709 0.0 0.0 0.0 0.0015233196 0.00051882295 9.4399749e-5 6.879389e-6 0.0094303151 -0.83284759 0.000518823 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0 0.0 0.0015233196 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 0.010219256 -0.83312459 0.0 0.0015233196 0.0011365783;
             0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0015233196 0.0 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.024511464 0.014477463 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.067400493 0.0 0.024511464;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.79244709 0.0015233196 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.011937727 0.0 -0.83341414 0.0011365783;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0011365784 0.0 0.0 0.0 0.0 0.0 0.0 0.024511464 0.0 0.0015233196 0.00051882295 9.4399749e-5 6.879389e-6 6.8793952e-6 9.4399775e-5 0.000518823 0.79244709 0.014477463 -0.83533612]

    elseif y_size == 33
        λ = [-0.85623069505 0.032502463 0.00078435813 0.79470408 0.0008046192 0.00057385784 0.0 0.00039669164 0.00057385784 0.0 0.0008046192 0.023382971 0.0 0.00078435813 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.00036481447 -0.85006657839 0.02675789 0.0 0.0008046192 0.00057385784 0.79470408 0.00039669164 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.023382971 0.0 0.0005540046 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.00036481447 0.0005540046 -0.8438983407899999 0.0 0.020840267 0.00057385784 0.0 0.00039669164 0.00057385784 0.79470408 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0 0.0005540046 0.0 0.023382971 0.0 0.00036481447 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.023382971 0.0 0.0 -0.08490958605 0.0 0.0 0.032502463 0.0 0.0 0.00078435813 0.0 0.023382971 0.0008046192 0.0 0.00057385784 0.0 0.00039669164 0.0 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0005540046 0.0 0.0 0.00036481447 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 -0.83805346908 0.014784895 0.0 0.00039669164 0.00057385784 0.0 0.0008046192 0.0 0.79470408 0.00078435813 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00036481447 0.023382971 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 -0.83346319354 0.0 0.0097866919 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.79470408 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.023382971 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.0 0.023382971 0.0 0.00036481447 0.0 0.0 -0.07874546938999999 0.0 0.0 0.02675789 0.0 0.0 0.0008046192 0.0 0.00057385784 0.023382971 0.00039669164 0.0 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0005540046 0.0 0.0 0.00036481447 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 0.00057385784 0.0 -0.8242503594799999 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.79470408 0.0005540046 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.023382971 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 0.00057385784 0.0 0.0097866919 -0.83346319354 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0 0.0005540046 0.79470408 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0 0.0 0.023382971 0.0 0.0 0.0 0.0 0.0;
             0.0 0.0 0.023382971 0.00036481447 0.0 0.0 0.0005540046 0.0 0.0 -0.07257723178999999 0.0 0.0 0.020840267 0.0 0.00057385784 0.0 0.00039669164 0.0 0.00057385784 0.023382971 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0005540046 0.0 0.0 0.00036481447 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 0.00057385784 0.0 0.00039669164 0.014784895 0.0 -0.83805346908 0.0 0.0 0.00078435813 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.79470408 0.00036481447 0.0 0.0 0.0 0.0 0.0 0.0 0.023382971 0.0 0.0 0.0 0.0;
             0.023382971 0.0 0.0 0.79470408 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.85623069505 0.0 0.0 0.0 0.032502463 0.0 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.00039669164 0.0 0.00057385784 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.00036481447 0.023382971 0.0 0.0005540046 0.0 0.0 0.00078435813 0.0 0.0 -0.06673236008 0.0 0.014784895 0.0 0.00039669164 0.0 0.00057385784 0.0 0.0008046192 0.0 0.023382971 0.00078435813 0.0 0.0 0.0005540046 0.0 0.0 0.00036481447 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 0.00057385784 0.0 0.00039669164 0.00057385784 0.0 0.020840267 0.0 0.0 -0.84389834079 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00036481447 0.0 0.79470408 0.0 0.0 0.0 0.0 0.0 0.0 0.023382971 0.0 0.0;
             0.0 0.0 0.0 0.00036481447 0.0 0.023382971 0.0005540046 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 -0.062142084539999996 0.0 0.0097866919 0.0 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.023382971 0.0 0.0005540046 0.0 0.0 0.00036481447 0.0 0.0 0.0;
             0.0 0.023382971 0.0 0.0 0.0 0.0 0.79470408 0.0 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 -0.85006657839 0.0 0.0 0.0 0.02675789 0.0 0.0 0.0008046192 0.0 0.00057385784 0.00039669164 0.0 0.00057385784 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0005540046 0.023382971 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.0 -0.052929250479999995 0.0 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.023382971 0.0005540046 0.0 0.0 0.00036481447 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 0.00057385784 0.0 0.00039669164 0.00057385784 0.0 0.0008046192 0.0 0.0 0.02675789 0.0 0.0 0.0 -0.85006657839 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0 0.79470408 0.0 0.0 0.0 0.0 0.023382971 0.0;
             0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0005540046 0.0 0.023382971 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.0 0.0097866919 0.0 -0.062142084539999996 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0005540046 0.023382971 0.0 0.00036481447 0.0 0.0 0.0;
             0.0 0.0 0.023382971 0.0 0.0 0.0 0.0 0.0 0.0 0.79470408 0.0 0.00036481447 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.0 -0.84389834079 0.0 0.0 0.020840267 0.0 0.00057385784 0.00039669164 0.0 0.00057385784 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0005540046 0.0 0.0 0.00078435813 0.023382971 0.0 0.0008046192 0.0 0.00057385784 0.0 0.00039669164 0.0 0.014784895 0.0 -0.06673236008 0.0 0.0 0.00078435813 0.0 0.0 0.0005540046 0.0 0.023382971 0.00036481447 0.0 0.0 0.0;
             0.00036481447 0.0005540046 0.00078435813 0.0 0.0008046192 0.00057385784 0.0 0.00039669164 0.00057385784 0.0 0.0008046192 0.0 0.0 0.00078435813 0.0 0.0 0.0 0.032502463 0.0 0.0 0.0 -0.85623069505 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.79470408 0.0 0.0 0.023382971;
             0.0 0.0 0.0 0.0 0.023382971 0.0 0.0 0.0 0.0 0.0 0.0 0.00036481447 0.79470408 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00078435813 0.0 0.0 -0.83805346908 0.0 0.014784895 0.00039669164 0.0 0.00057385784 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0005540046 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.023382971 0.00057385784 0.0 0.00039669164 0.0 0.00057385784 0.0 0.020840267 0.0 0.0 -0.07257723179 0.0 0.0 0.0005540046 0.0 0.0 0.00036481447 0.023382971 0.0 0.0;
             0.0 0.0 0.0 0.0 0.0 0.023382971 0.0 0.0 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.79470408 0.0005540046 0.0 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 -0.83346319354 0.0097866919 0.0 0.00057385784 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.023382971 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0005540046 0.79470408 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 -0.8242503594799999 0.0 0.00057385784 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0005540046 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.0 0.00039669164 0.023382971 0.00057385784 0.0 0.0008046192 0.0 0.0 0.02675789 0.0 0.0 -0.07874546938999999 0.0 0.0 0.00036481447 0.0 0.023382971 0.0;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.023382971 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.79470408 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.0097866919 0.0 -0.83346319354 0.0008046192 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.023382971 0.00036481447 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00078435813 0.79470408 0.0 0.0008046192 0.0 0.00057385784 0.00039669164 0.0 0.014784895 -0.83805346908 0.0 0.00078435813 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0005540046 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.0 0.00039669164 0.0 0.00057385784 0.0 0.0008046192 0.023382971 0.0 0.00078435813 0.0 0.0 0.032502463 0.0 0.0 -0.08490958605 0.0 0.0 0.023382971;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00036481447 0.0 0.023382971 0.0 0.0005540046 0.0 0.0 0.0 0.00078435813 0.0 0.0 0.0008046192 0.79470408 0.00057385784 0.00039669164 0.0 0.00057385784 0.020840267 0.0 -0.8438983407899999 0.0005540046 0.00036481447;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0005540046 0.0 0.023382971 0.0 0.00078435813 0.0 0.0 0.0008046192 0.0 0.00057385784 0.00039669164 0.79470408 0.00057385784 0.0008046192 0.0 0.02675789 -0.85006657839 0.00036481447;
             0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00036481447 0.0 0.0 0.0 0.0005540046 0.0 0.0 0.0 0.00078435813 0.0 0.023382971 0.0008046192 0.0 0.00057385784 0.00039669164 0.0 0.00057385784 0.0008046192 0.79470408 0.00078435813 0.032502463 -0.85623069505]
    end
    return λ
end


text:
using CRC

width = 8
poly = 0x31
init = 0x00
refin = false
refout = false
xorout = 0x00
check = 0x00

s = spec(width, poly, init, refin, refout, xorout, check)
d = b"xyz"

for t in (NoTables, Single, Multiple)
    println(t)
    c = crc(s, tables=t)
    c(d)
    @time c(d)
end


text:
fname = :TRMA97_abs
c = Combi(hessian_dense,PDataMA97,solve_modelTRDiagAbs,preprocessMA97,decreaseFact,Tparam())
include("Template.jl")


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/cudnn/secure/8.3.2/local_installers/11.5/cudnn-linux-x86_64-8.3.2.44_cuda11.5-archive.tar.xz",
                      "5500953c08c5e5d1dddcfda234f9efbddcdbe43a53b26dc0a82c723fa170c457")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/cudnn/secure/8.3.2/local_installers/11.5/cudnn-linux-ppc64le-8.3.2.44_cuda11.5-archive.tar.xz",
                      "0581bce48023a3ee71c3a819aaefcabe693eca18b61e2521dc5f8e6e71567b1b")],
    Platform("aarch64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/cudnn/secure/8.3.2/local_installers/11.5/cudnn-linux-sbsa-8.3.2.44_cuda11.5-archive.tar.xz",
                      "7eb8c96bfeec98e8aa7cea1e95633d2a9481fc99040eb0311d31bf137a7aa6ea")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/cudnn/secure/8.3.2/local_installers/11.5/cudnn-windows-x86_64-8.3.2.44_cuda11.5-archive.zip",
                      "9e36eef803f1cf9ab24846dc133a3014fdc548775ee29073e8466d415957a1c0")],
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module DarcyTests

using Test
using Gridap
import Gridap: ∇, divergence
using LinearAlgebra

u(x) = VectorValue(2*x[1],x[1]+x[2])

divergence(::typeof(u)) = (x) -> 3

p(x) = x[1]-x[2]

∇p(x) = VectorValue(1,-1)

∇(::typeof(p)) = ∇p

f(x) = u(x) + ∇p(x)

domain = (0,1,0,1)
partition = (4,4)
order = 1
model = CartesianDiscreteModel(domain,partition)

V = FESpace(
  reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64},
  conformity=:Hdiv, model=model, dirichlet_tags=[5,6])

Q = FESpace(
  reffe=:QLagrangian, order=order, valuetype=Float64,
  conformity=:L2, model=model)

U = TrialFESpace(V,u)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

neumanntags = [7,8]
btrian = BoundaryTriangulation(model,neumanntags)
degree = 2*(order+1)
bquad = CellQuadrature(btrian,degree)
nb = get_normal_vector(btrian)

function a(x,y)
  u, p = x
  v, q = y
  u⋅v - p*(∇⋅v) + q*(∇⋅u)
end

function l(y)
  v, q = y
  v⋅f + q*(∇⋅u)
end

function l_Γ(y)
  v, q = y
  -(v⋅nb)*p
end

t_Ω = AffineFETerm(a,l,trian,quad)
t_Γ = FESource(l_Γ,btrian,bquad)
op = AffineFEOperator(X,Y,t_Ω,t_Γ)
xh = solve(op)
uh, ph = xh

eu = u - uh
ep = p - ph

l2(v) = v⋅v
h1(v) = v*v + ∇(v)⋅∇(v)

eu_l2 = sum(integrate(l2(eu),trian,quad))
ep_l2 = sum(integrate(l2(ep),trian,quad))
ep_h1 = sum(integrate(h1(ep),trian,quad))

tol = 1.0e-9
@test eu_l2 < tol
@test ep_l2 < tol
@test ep_h1 < tol

end # module


text:
# Minimal sanity checks for MakieLayout
@testset "Layoutables constructors" begin
    scene, layout = layoutscene()
    ax = layout[1, 1] = Axis(scene)
    cb = layout[1, 2] = Colorbar(scene)
    gl2 = layout[2, :] = MakieLayout.GridLayout()
    bu = gl2[1, 1] = Button(scene)
    sl = gl2[1, 2] = Slider(scene)

    scat = scatter!(ax, rand(10))
    le = gl2[1, 3] = Legend(scene, [scat], ["scatter"])

    to = gl2[1, 4] = Toggle(scene)
    te = layout[0, :] = Label(scene, "A super title")
    me = layout[end + 1, :] = Menu(scene, options=["one", "two", "three"])
    tb = layout[end + 1, :] = Textbox(scene)
    @test true
end


text:
# implement julia_commander
# for the moment, only help is implemented

function filterranges(u::Vector{UnitRange{Int}})
    v = Vector{Int}()
    while length(u) > 1
        if all(m -> (u[1] ⊈ m), u[2:end])
            push!(v, u[1][end])
        end
        u = u[findall(m -> m ⊈ u[1], u[2:end]).+1]
    end
    if length(u) == 1
        push!(v,u[1][end])
    end
    return v
end

const STYLES = [
    r"```.+?```"s, r"`.+?`", r"~~.+?~~", r"_.+?_", r"__.+?__",
    r"\*.+?\*", r"\*\*.+?\*\*", r"\*\*\*.+?\*\*\*",
]

function split_message_fixed(text::AbstractString, chunk_limit::Int=2000; extraregex::Vector{Regex}=Vector{Regex}())
    length(text) <= chunk_limit && return String[text]
    chunks = String[]

    while !isempty(text)
        mranges = vcat(findall.(union(STYLES,extraregex),Ref(text))...)

        stop = maximum(filter(i -> length(text[1:i]) ≤ chunk_limit, filterranges(mranges)))

        # give up if first chunk cannot be broken down
        if stop == 0
            push!(chunks, strip(text))
            @warn "message could not be broken down into chunks smaller than the desired length $chunk_limit"
            return chunks
        end

        push!(chunks, strip(text[1:stop]))
        text = strip(text[stop+1:end])
    end

    return chunks
end

function parse_doc(doc::AbstractString)
    doc = replace(doc, r"\n\n\n+" => "\n\n")
    for m in eachmatch(r"(^|\n)(#+ |!!! )(.*)\n",doc)
        if m.captures[2] == "# "
            doc = replace(doc, m.match => m.captures[1]*"**"*m.captures[3]*"**\n"*"≡"^(length(m.captures[3])), count = 1)
        elseif m.captures[2] == "!!! "
            doc = replace(doc, m.match => m.captures[1]*"__"*m.captures[3]*"__\n", count = 1)
        else
            doc = replace(doc, m.match => m.captures[1]*"*"*m.captures[3]*"*\n"*"-"^(length(m.captures[3])), count = 1)
        end
    end
    doc = replace(doc, r"(```.+)\n" => "```julia\n")
    doc = replace(doc, "```\n" => "```")
    for m in eachmatch(r"\[([^ ]*)\]\(@ref\)",doc)
        doc = replace(doc, m.match => "`"*m.captures[1]*"`", count = 1)
    end
    return doc
end

function commander(c::Client, m::Message, ::Val{:julia_doc})
    # @info "julia_commander called"
    # @info "Message content" m.content m.author.username m.author.discriminator
    startswith(m.content, COMMAND_PREFIX * "j") || return
    regex = Regex(COMMAND_PREFIX * raw"j(\?| help| doc)? *(.*)$")
    matches = match(regex, m.content)
    if matches === nothing || matches.captures[1] in (" help", nothing)
        help_commander(c, m, Val(:julia_doc))
    elseif matches.captures[1] in ("?", " doc")
        handle_julia_help_commander(c, m, matches.captures[2])
    else
        reply(c, m, "Sorry, I don't understand the request; use `j help` for help")
    end
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:julia_doc})
    # @info "Sending help for message" m.id m.author
    reply(c, m, """
        For the moment, only "doc" is accept, for showing the docstring.

        How to use the `j` command:
        ```
        j help
        j? <name>
        j doc <name>
        ```
        `j help` returns this help
        `j? <name>` and `j doc <name>` return the documentation for `<name>`
        """)
    return nothing
end

function handle_julia_help_commander(c::Client, m::Message, name)
    # @info "julia_help_commander called"
    try
        doc = string(eval(:(Base.Docs.@doc $(Symbol(name)))))
        doc = parse_doc(doc)
        docs = split_message_fixed(doc, extraregex = [r"\n≡.+\n", r"n-.+\n"])
        for doc_chunck in docs
            reply(c, m, doc_chunck)
        end
    catch ex
        @show ex
        reply(c, m, "Sorry, it didn't work.")
    end
    return nothing
end


text:
using PMIx
using Test

using prrte_jll

if !haskey(ENV, "PRTE_LAUNCHED")
    @testset "Outside PRTE session" begin
        @test PMIx.initialized() == false
        @test_throws PMIx.PMIxException PMIx.init()
    end

    @testset "Examples" begin
        example = realpath(joinpath(@__DIR__, "..", "examples", "client.jl"))
        prrte_jll.prterun() do prterun
            cmd = `$prterun -np 2 $(Base.julia_cmd()) $example`
            @test success(pipeline(cmd, stdout=stdout, stderr=stderr))
        end
    end

    @info "relaunching under PRTE"
    current_file = @__FILE__ # bug in 1.5 can't be directly interpolated
    jlcmd = `$(Base.julia_cmd()) $(current_file)`

    prrte_jll.prterun() do prterun
        cmd = `$prterun -np 2 $jlcmd`
        @test success(pipeline(cmd, stdout=stdout, stderr=stderr))
    end
    exit()
end

# These tests now execute on multiple ranks
const myproc = PMIx.init()

@testset "Inside PRTE session" begin
    @test PMIx.initialized() == true
    @test myproc.rank ∈ (0, 1)
end

@testset "Values" begin
    optional = PMIx.Value(true, PMIx.API.PMIX_BOOL)
    @test convert(Bool, optional)

    optional = PMIx.Value(false, PMIx.API.PMIX_BOOL)
    @test !convert(Bool, optional)
end

@testset "Job size" begin
    proc = PMIx.Proc(myproc.nspace, PMIx.API.PMIX_RANK_WILDCARD)
    optional = PMIx.Value(true, PMIx.API.PMIX_BOOL)
    value = PMIx.get(proc, PMIx.API.PMIX_JOB_SIZE, [PMIx.Info(PMIx.API.PMIX_OPTIONAL, optional)])
    @test PMIx.get_number(value) == 2
end

PMIx.finalize()

max_issues_repo_name: JuliaParallel/PMIx

max_forks_repo_name: JuliaParallel/PMIx


text:
# ### m12.5t.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame);

@model ppl12_5(pulled_left, actor, block, condition, prosoc_left) = begin
    # Total num of y
    N = length(pulled_left)

    # Separate σ priors for each actor and block
    σ_actor ~ truncated(Cauchy(0, 1), 0, Inf)
    σ_block ~ truncated(Cauchy(0, 1), 0, Inf)

    # Number of unique actors in the data set
    N_actor = length(unique(actor)) #7
    N_block = length(unique(block))

    # Vector of actors (1,..,7) which we'll set priors on
    α_actor ~ filldist(Normal(0, σ_actor), N_actor)
    α_block ~ filldist(Normal(0, σ_block), N_block)

    # Prior for intercept, prosoc_left, and the interaction
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logitp = α .+ α_actor[actor] + α_block[block] .+
            (βp .+ βpC * condition) .* prosoc_left

    pulled_left .~ BinomialLogit.(1, logitp)

end

# Sample

m12_5t = ppl12_5(
    Vector{Int64}(df.pulled_left),
    Vector{Int64}(df.actor),
    Vector{Int64}(df.block),
    Vector{Int64}(df.condition),
    Vector{Int64}(df.prosoc_left)
);
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns12_5t = mapreduce(c -> sample(m12_5t, sampler, nsamples), chainscat, 1:nchains)

# Results rethinking

m125rethinking = "
             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
a_actor[1]  -1.19   0.98      -2.66       0.35  3451    1
a_actor[2]   4.14   1.59       1.87       6.38  5514    1
a_actor[3]  -1.49   0.99      -2.96       0.05  3493    1
a_actor[4]  -1.50   0.98      -3.00       0.01  3340    1
a_actor[5]  -1.19   0.99      -2.68       0.34  3447    1
a_actor[6]  -0.25   0.99      -1.79       1.25  3361    1
a_actor[7]   1.30   1.01      -0.21       2.89  3673    1
a_block[1]  -0.19   0.23      -0.56       0.10  3825    1
a_block[2]   0.04   0.19      -0.24       0.34  9346    1
a_block[3]   0.06   0.19      -0.22       0.37  9202    1
a_block[4]   0.00   0.18      -0.29       0.29 11314    1
a_block[5]  -0.03   0.19      -0.33       0.25 10596    1
a_block[6]   0.11   0.21      -0.16       0.45  6040    1
a            0.47   0.97      -1.01       1.99  3273    1
bp           0.83   0.26       0.40       1.24 13225    1
bpc         -0.15   0.30      -0.61       0.36  8492    1
sigma_actor  2.27   0.91       1.03       3.35  5677    1
sigma_block  0.23   0.18       0.01       0.44  2269    1
";

# End of m12.5t.jl


text:
import IterationControl

# # SAVE

struct Save{K}
    filename::String
    kwargs::K
end

# constructor:
function Save(filename="machine.jlso"; kwargs...)
    if :filename in keys(kwargs)
        error("`filename` is not a keyword argumnent of `Save`. "*
              "Usage: `Save(filename::String; kwargs...)`. ")
    end
    return Save(filename, kwargs)
end

IterationControl.@create_docs(Save,
             header="Save(filename=\"machine.jlso\"; kwargs...)",
             example="Save(\"run3/machine.jlso\", compression=:gzip)",
             body="Save the current state of the machine being iterated to "*
             "disk, using the provided `filename`, decorated with a number, "*
             "as in \"run3/machine_42.jlso\". The specified `kwargs` "*
             "are passed to the model-specific serializer "*
             "(JLSO for most Julia models).\n\n"*
             "For more on what is meant by \"the machine being iterated\", "*
             "see [`IteratedModel`](@ref).")

function IterationControl.update!(c::Save,
                                  ic_model,
                                  verbosity,
                                  n,
                                  state=(filenumber=0, ))
    filenumber = state.filenumber + 1
    root, suffix = splitext(c.filename)
    filename = string(root, filenumber, suffix)
    train_mach = IterationControl.expose(ic_model)
    verbosity > 0 && @info "Saving \"$filename\". "
    MLJSerialization.save(filename, train_mach, c.kwargs...)
    return (filenumber=filenumber, )
end


text:
@testset "Unscented Kalman filter" begin

    @testset "Time update" begin

        x = randn(6)
        PL = randn(6,6)
        P = PL'PL
        QL = randn(6,6)
        Q = QL'QL
        F = randn(6,6)
        f(x) = F * x

        tu = time_update(x, P, F, Q)
        tu_ukf = time_update(x, P, f, Q)
        @test get_covariance(tu_ukf) ≈ get_covariance(tu)
        @test get_state(tu_ukf) ≈ get_state(tu)

        f!(y, x) = mul!(y, F, x)
        tu_inter = UKFTUIntermediate(6)
        tu = time_update!(tu_inter, x, P, f!, Q)
        @test get_covariance(tu_ukf) ≈ get_covariance(tu)
        @test get_state(tu_ukf) ≈ get_state(tu)

    end

    @testset "Measurement update" begin

        x = randn(6)
        PL = randn(6,6)
        P = PL'PL
        RL = randn(3,3)
        R = RL'RL
        y = randn(3)
        H = randn(3,6)
        h(x) = H * x

        mu = measurement_update(x, P, y, H, R)
        mu_ukf = @inferred measurement_update(x, P, y, h, R)
        @test get_covariance(mu_ukf) ≈ get_covariance(mu)
        @test get_state(mu_ukf) ≈ get_state(mu)

        h!(y, x) = mul!(y, H, x)
        mu_inter = @inferred UKFMUIntermediate(6, 3)
        mu_ukf_inplace = @inferred measurement_update!(mu_inter, x, P, y, h!, R)
        @test get_covariance(mu_ukf_inplace) ≈ get_covariance(mu)
        @test get_state(mu_ukf_inplace) ≈ get_state(mu)
    end
end


text:
################################################################################
# set_scenario!
################################################################################

function set_scenario!(vis::Visualizer, sce::Scenario{T};
	roadway_vis_opts::RoadwayVisualizationOptions=RoadwayVisualizationOptions{T}(),
    player_vis_opts::PlayerVisualizationOptions=PlayerVisualizationOptions{T}(),
	color::Vector{Symbol}=Vector{Symbol}()) where {T}
    clean!(vis)
	set_roadway!(vis, sce.roadway.opts, vis_opts=roadway_vis_opts)
	set_player!(vis, sce.player, vis_opts=player_vis_opts, color=color)
    return nothing
end


text:
struct NoSuchAttr end

struct Attr{S} end

@inline getattr(x, f::Symbol) = getattr(x, Attr{f}())

@inline setattr!(x, f::Symbol, y) = setattr!(x, Attr{f}(), y)

@inline literal_getattr(x, ::Attr{F}) where F = getproperty(x, F)

@inline literal_setattr!(x, ::Attr{F}, y) where F = setproperty!(x, F, y)

@generated function trygetfield(x, ::Attr{F}) where F
    hasfield = F in fieldnames(x)
    quote
        $(Expr(:meta, :inline))
        $(hasfield ? :(Base.getfield(x, $(Meta.quot(F)))) : :(NoSuchAttr()))
    end
end

@generated function trygetfield(x, f::Symbol)
    fields = Meta.quot.(fieldnames(x))
    quote
        $(Expr(:meta, :inline))
        $([:(f === $a && return Base.getfield(x, $a)) for a in fields]...)
        NoSuchAttr()
    end
end

@generated function trysetfield!(x, ::Attr{F}, y) where F
    hasfield = F in fieldnames(x)
    quote
        $(Expr(:meta, :inline))
        $(hasfield ? :(Base.setfield!(x, $(Meta.quot(F)), y)) : :(NoSuchAttr()))
    end
end

@generated function trysetfield!(x, f::Symbol, y)
    fields = Meta.quot.(fieldnames(x))
    quote
        $(Expr(:meta, :inline))
        $([:(f === $a && return Base.setfield!(x, $a, y)) for a in fields]...)
        NoSuchAttr()
    end
end

@inline attrnames(::Type{<:NamedTuple{Names}}) where Names = Names

function attrnames(::Type{T}) where T
    fields = fieldnames(T)
    attrs = collect(Symbol, fields)
    attrset = Set(attrs)
    for m in methods(getattr, (T, Attr)).ms
        sig = m.sig
        while isa(sig, UnionAll)
            sig = sig.body
        end
        attr = sig.parameters[3].parameters[1]
        if !(attr isa TypeVar)
            attr::Symbol
            if !(attr in attrset)
                push!(attrset, attr)
                push!(attrs, attr)
            end
        end
    end

    if length(attrs) == length(fields)
        fields
    else
        sort!(view(attrs, (length(fields) + 1):length(attrs)))
        Tuple(attrs)
    end
end

@inline attrnames(x) = attrnames(typeof(x))

@inline function default_literal_getattr(x, f::Attr)
    res = trygetfield(x, f)
    if res ≡ NoSuchAttr()
        getattr(x, f)
    else
        res
    end
end

@inline function default_literal_setattr!(x, f::Attr, y)
    res = trysetfield!(x, f, y)
    if res ≡ NoSuchAttr()
        setattr!(x, f, y)
    else
        res
    end
end

@inline function default_getproperty(x, f::Symbol)
    res = trygetfield(x, f)
    if res ≡ NoSuchAttr()
        getattr(x, Attr{f}())
    else
        res
    end
end

@inline function default_setproperty!(x, f::Symbol, y)
    res = trysetfield!(x, f, y)
    if res ≡ NoSuchAttr()
        setattr!(x, Attr{f}(), y)
    else
        res
    end
end


text:
"""
    submodule CNumbers
Provides C-style numbers, whose conversion rule is different from julia's:
* Downcast to a smaller type is simply truncating, with no bounds check. (Although C11 6.3.1.3p3  permits throwing an InexactError, this may not be expected to C programmers)
* Integers promoted to int before performing arithmetic, possibly converted back later(C11 6.3 Conversions)
* Floating-point to integer conversion is plain `trunc`
Also it exports C-style lowercase type names like `float`, `int`, or `unsigned_int`.
"""

"""
    struct CNumber{T <: Number} <: Number
Wrapper of a C number, conforming to C-style integer/floating-point number conversion/arithmetic. Only intended for internal usage.
To get the wrapped value, call `value`.
"""
struct CNumber{T <: Number} <: Number
    v::T
end

"Get the wrapped value of a CNumber."
value(x::CNumber) = x.v

"Get the wrapped type of a CNumber."
basetype(::CNumber{T}) where T = T

"Integer rank is for comparing types. Defined on julia's primitive types."
function rank end

# Aliases

for (ind, (name, type)) in [
    (:bool, Bool),
    (:char, Cchar),
    (:unsigned_char, Cuchar),
    (:short, Cshort),
    (:unsigned_short, Cushort),
    (:int, Cint),
    (:unsigned_int, Cuint),
    (:long, Clong),
    (:unsigned_long, Culong),
    (:long_long, Clonglong),
    (:unsigned_long_long, Culonglong),
    (:float, Cfloat),
    (:double, Cdouble)] |> enumerate
    @eval begin
        const $name = CNumber{$type}
        # Avoid redefinition
        if !hasmethod(rank, Tuple{Type{$type}})
            rank(::Type{$type}) = $ind ÷ 2
        end
    end
end

# Convert to normal julia types (outer constructor is not defined)
Base.convert(::Type{T}, x::CNumber) where {T <: Number} = Base.convert(T, value(x))
# Convert between integer types by truncating
Base.convert(::Type{CNumber{T}}, x::CNumber{<:Integer}) where {T <: Integer} = CNumber{T}(value(x) % T)
# Convert to floating point number by direct julian conversion
Base.convert(::Type{CNumber{T}}, x::CNumber) where {T <: AbstractFloat} = CNumber{T}(Base.convert(T, value(x)))
# Convert from floating point number to integer by truncating
Base.convert(::Type{CNumber{T}}, x::CNumber{<:AbstractFloat}) where {T <: Integer} = CNumber{T}(trunc(value(x)))
# Convert to bool is !iszero
Base.convert(::Type{Bool}, x::CNumber) = !iszero(value(x))
Base.convert(::Type{bool}, x::CNumber) = bool(Base.convert(Bool, x))

"Integers with rank smaller than int are promoted to int"
promote_integer(::Type{T}) where T = rank(T) ≥ rank(Cint) ? T : Cint

signedness(::Type{<:Signed}) = true
signedness(::Type{<:Unsigned}) = false

function c_promote(::Type{T1}, ::Type{T2}) where {T1 <: AbstractFloat,T2}
    return promote_type(T1, T2)
end

function c_promote(::Type{T1}, ::Type{T2}) where {T1,T2 <: AbstractFloat}
    return promote_type(T1, T2)
end

function c_promote(::Type{T1}, ::Type{T2}) where {T1,T2}
    T3, T4 = promote_integer(T1), promote_integer(T2)
    if !(signedness(T3) ⊻ signedness(T4))
        return rank(T3) > rank(T4) ? T3 : T4
    elseif T3 <: Unsigned
        U, S = T3, T4
    else
        U, S = T4, T3
    end

    if rank(U) ≥ rank(S)
        return U
    elseif sizeof(U) < sizeof(S)
        return S
    else
        return U
    end
end

function Base.promote_rule(::Type{CNumber{T1}}, ::Type{CNumber{T2}}) where {T1,T2}
    return CNumber{c_promote(T1, T2)}
end

# Always favor CNumber over julian nubers
function Base.promote_rule(::Type{CNumber{T1}}, ::Type{T2}) where {T1, T2}
    Base.promote_rule(CNumber{T1}, CNumber{T2})
end

# Arithmetic

for op in [:+, :-, :*, :/, :%, :&, :|, :⊻, :<<, :>>] # XOR is ⊻ in julia
    @eval function Base.$op(a::CNumber{T}, b::CNumber{T}) where T
        CNumber{T}($op(value(a), value(b)))
    end
end

for op in [:+, :-, :~]
    @eval function Base.$op(a::CNumber{T}) where T
        CNumber{T}($op(value(a)))
    end
end

# Return Bool
for op in [:(==), :!=, :<, :>, :<=, :(>=)]
    @eval function Base.$op(a::CNumber, b::CNumber)
        a, b = Base.promote(a, b)
        Base.$op(value(a), value(b))
    end
end

function Base.:!(a::CNumber)
    iszero(a)
end

"@post_inc i is equivalent to i++"
macro post_inc(i)
	quote
		local v = $(esc(i))
		$(esc(i)) += 1
		v
	end
end

"@post_dec i is equivalent to i--"
macro post_dec(i)
	quote
		local v = $(esc(i))
		$(esc(i)) -= 1
		v
	end
end


text:
immutable Quantity{T <: Number, UNIT <: Unit} <: Number
    _::T
end
Quantity(n::Number, unit::AbstractUnit) = Quantity{typeof(n), typeof(unit)}(n)

""" Physical unit of the input quantity """
unit(q::Quantity) = unit(typeof(q))
unit{T <: Number, U <: Unit}(q::Type{Quantity{T, U}}) = U()
unit{A <: AbstractArray}(q::Type{A}) = unit(eltype(A))
unit(q::AbstractArray) = unit(eltype(q))

Base.convert(::Type{Number}, q::Quantity) = q._
Base.show(io::IO, q::Quantity) = print(io, q._, "⋅", unit(q))
Base.eltype{T <: Number, U <: Unit}(::Type{Quantity{T, U}}) = T
Base.eltype(q::Quantity) = eltype(typeof(q))

*(x::Number, unit::AbstractUnit) = Quantity{typeof(x), typeof(unit)}(x)
*(unit::AbstractUnit, x::Number) = Quantity{typeof(x), typeof(unit)}(x)
/(x::Number, unit::AbstractUnit) = Quantity{typeof(x), typeof(unit^-1)}(x)
/(unit::AbstractUnit, x::Number) =
    Quantity{typeof(one(x)/x), typeof(unit)}(one(x)/x)

dimensionality(q::Quantity) = dimensionality(unit(q))
dimensionality(q::AbstractArray) = dimensionality(unit(q))
dimensionality{T <: Number, U <: Unit}(::Type{Quantity{T, U}}) =
    dimensionality(U)
unit_system(q::Quantity) = unit_system(unit(q))
unit_system{T <: Number, U <: Unit}(::Type{Quantity{T, U}}) = unit_system(U)
unit_system{A <: AbstractArray}(::Type{A}) = unit_system(eltype(A))
unit_system(q::AbstractArray) = unit_system(eltype(q))

""" converts quantity to given unit """
function conversion(a::Quantity, b::AbstractUnit)
    dimensionality(a) == dimensionality(b) ||
        error("Dimensionality do not match")
    unit_system(unit(a)) == unit_system(b) && return a
    const slope, offset = conversion_factors(unit(a), b)
    Quantity{eltype(a), typeof(b)}(slope * a._ + offset)
end
""" converts quantity to another system """
function conversion(a::Quantity, system::Symbol)
    unit_system(a) == system && return a
    const U = Unit{system, dimensionality(a)}()
    const slope, offset = conversion_factors(unit(a), U)
    const value = a._ * slope + offset
    Quantity{typeof(value), typeof(U)}(value)
end
""" converts quantity to another system """
function conversion(u::Unit, system::Symbol)
    unit_system(u) == system && return u
    const U = Unit{system, dimensionality(u)}()
    const slope, offset = conversion_factors(u, U)
    const value = 1 * slope + offset
    Quantity{typeof(value), typeof(U)}(value)
end