text:
struct Vec{N, T<:VecTypes}
    data::LVec{N, T}
end

# Constructors
@inline Vec(v::NTuple{N, T}) where {N, T<:VecTypes} = Vec(VE.(v))
@inline Vec(v::Vararg{T, N}) where {N, T<:VecTypes} = Vec(v)
@inline Vec(v::Vec) = v
# Numbers defines this and it is needed in power_by_squaring...
Base.copy(v::Vec) = v

# No throwing versions of convert
@inline _unsafe_convert(::Type{T}, v) where {T <: IntegerTypes} = v % T
@inline _unsafe_convert(::Type{T}, v) where {T <: VecTypes} = convert(T, v)
@inline constantvector(v::T1, ::Type{Vec{N, T2}}) where {N, T1, T2} =
    Vec(Intrinsics.constantvector(_unsafe_convert(T2, v), Intrinsics.LVec{N, T2}))

@inline Vec{N, T}(v::Vec{N, T}) where {N, T<:VecTypes} = v
@inline Vec{N, T}(v::Vec{N, T}) where {N, T<:FloatingTypes} = v
@inline Vec{N, T1}(v::T2) where {N, T1<:VecTypes, T2<:VecTypes} = constantvector(v, Vec{N, T1})
@inline Vec{N, T1}(v::Vec{N, T2}) where {N, T1<:Union{IntegerTypes, Ptr}, T2<:Union{IntegerTypes, Ptr}} =
    convert(Vec{N, T1}, v)

@inline Base.convert(::Type{Vec{N,T}}, v::Vec{N,T}) where {N,T} = v
@inline function Base.convert(::Type{Vec{N, T1}}, v::Vec{N, T2}) where {T1, T2, N}
    if T1 <: Union{IntegerTypes, Ptr}
        if T2 <: Union{IntegerTypes, Ptr, Bool}
            if sizeof(T1) < sizeof(T2)
                return Vec(Intrinsics.trunc(Intrinsics.LVec{N, T1}, v.data))
            elseif sizeof(T1) == sizeof(T2)
                return Vec(Intrinsics.bitcast(Intrinsics.LVec{N, T1}, v.data))
            else
                if T2 <: UIntTypes
                    return Vec(Intrinsics.zext(Intrinsics.LVec{N, T1}, v.data))
                else
                    return Vec(Intrinsics.sext(Intrinsics.LVec{N, T1}, v.data))
                end
            end
        elseif T2 <: FloatingTypes
            if T1 <: UIntTypes
                return Vec(Intrinsics.fptoui(Intrinsics.LVec{N, T1}, v.data))
            elseif T1 <: IntTypes
                return Vec(Intrinsics.fptosi(Intrinsics.LVec{N, T1}, v.data))
            end
        end
    end
    if T1 <: FloatingTypes
        if T2 <: UIntTypes
            return Vec(Intrinsics.uitofp(Intrinsics.LVec{N, T1}, v.data))
        elseif T2 <: IntTypes
            return Vec(Intrinsics.sitofp(Intrinsics.LVec{N, T1}, v.data))
        elseif T2 <: FloatingTypes
            if sizeof(T1) < sizeof(T2)
                return Vec(Intrinsics.fptrunc(Intrinsics.LVec{N, T1}, v.data))
            else
                return Vec(Intrinsics.fpext(Intrinsics.LVec{N, T1}, v.data))
            end
        end
    end
    _unreachable()
end
@noinline _unreachable() = error("unreachable")

Base.Tuple(v::Vec) = map(i -> i.value, v.data)
Base.NTuple{N, T}(v::Vec{N}) where {T, N} = map(i -> convert(T, i.value), v.data)

Base.eltype(::Type{Vec{N,T}}) where {N,T} = T
Base.ndims( ::Type{Vec{N,T}}) where {N,T} = 1
Base.length(::Type{Vec{N,T}}) where {N,T} = N
Base.size(  ::Type{Vec{N,T}}) where {N,T} = (N,)
Base.size(  ::Type{Vec{N,T}}, n::Integer) where {N,T} = n > N ? 1 : (N,)[n]

Base.eltype(V::Vec) = eltype(typeof(V))
Base.ndims(V::Vec) = ndims(typeof(V))
Base.length(V::Vec) = length(typeof(V))
Base.size(V::Vec) = size(typeof(V))
Base.size(V::Vec, n::Integer) = size(typeof(V), n)

if VERSION <= v"1.4.0-rc1.0"
    function Base.show(io::IO, v::Vec{N,T}) where {N,T}
        print(io, "<$N x $T>[")
        join(io, [x.value for x in v.data], ", ")
        print(io, "]")
    end
else
    # This crashes on pre 1.4-rc2
    function Base.show(io::IO, v::Vec{N,T}) where {N,T}
        io = IOContext(io, :typeinfo => eltype(v))
        print(io, "<$N x $T>[")
        join(io, [sprint(show, x.value; context=io) for x in v.data], ", ")
        print(io, "]")
    end
end

@inline Base.checkbounds(v::Vec, i::IntegerTypes) =
(i < 1 || i > length(v.data)) && Base.throw_boundserror(v, i)

function Base.getindex(v::Vec, i::IntegerTypes)
    @boundscheck checkbounds(v, i)
    return Intrinsics.extractelement(v.data, i-1)
end

@inline function Base.setindex(v::Vec{N,T}, x, i::IntegerTypes) where {N,T}
    @boundscheck checkbounds(v, i)
    Vec(Intrinsics.insertelement(v.data, _unsafe_convert(T, x), i-1))
end

Base.zero(::Type{Vec{N,T}}) where {N, T} = Vec{N,T}(zero(T))
Base.zero(::Vec{N,T}) where {N, T} = zero(Vec{N, T})
Base.one(::Type{Vec{N,T}}) where {N, T} = Vec{N, T}(one(T))
Base.one(::Vec{N,T}) where {N, T} = one(Vec{N, T})

Base.reinterpret(::Type{Vec{N, T}}, v::Vec) where {T, N} = Vec(Intrinsics.bitcast(Intrinsics.LVec{N, T}, v.data))
Base.reinterpret(::Type{Vec{N, T}}, v::ScalarTypes) where {T, N} = Vec(Intrinsics.bitcast(Intrinsics.LVec{N, T}, v))
Base.reinterpret(::Type{T}, v::Vec) where {T} = Intrinsics.bitcast(T, v.data)

const FASTMATH = Intrinsics.FastMathFlags(Intrinsics.FastMath.fast)

###################
# Unary operators #
###################

const UNARY_OPS = [
    (:sqrt           , FloatingTypes , Intrinsics.sqrt)       ,
    (:sin            , FloatingTypes , Intrinsics.sin)        ,
    (:trunc          , FloatingTypes , Intrinsics.trunc)      ,
    (:cos            , FloatingTypes , Intrinsics.cos)        ,
    (:exp            , FloatingTypes , Intrinsics.exp)        ,
    (:exp2           , FloatingTypes , Intrinsics.exp2)       ,
    (:log            , FloatingTypes , Intrinsics.log)        ,
    (:log10          , FloatingTypes , Intrinsics.log10)      ,
    (:log2           , FloatingTypes , Intrinsics.log2)       ,
    (:abs            , FloatingTypes , Intrinsics.fabs)       ,
    (:floor          , FloatingTypes , Intrinsics.floor)      ,
    (:ceil           , FloatingTypes , Intrinsics.ceil)       ,
    # (:rint         , FloatingTypes , Intrinsics)            ,
    # (:nearbyint    , FloatingTypes , Intrinsics)            ,
    (:round          , FloatingTypes , Intrinsics.round)      ,

    (:bswap          , IntegerTypes  , Intrinsics.bswap)      ,
    (:count_ones     , IntegerTypes  , Intrinsics.ctpop)      ,
    (:leading_zeros  , IntegerTypes  , Intrinsics.ctlz)       ,
    (:trailing_zeros , IntegerTypes  , Intrinsics.cttz)       ,
]

if isdefined(Base, :bitreverse)
    push!(UNARY_OPS,
        (:bitreverse   , IntegerTypes  , Intrinsics.bitreverse)
    )
end

for (op, constraint, llvmop) in UNARY_OPS
    @eval @inline (Base.$op)(x::Vec{<:Any, <:$constraint}) =
        Vec($(llvmop)(x.data))
end

Base.:+(v::Vec{<:Any, <:ScalarTypes}) = v
Base.:-(v::Vec{<:Any, <:IntegerTypes}) = zero(v) - v
Base.:-(v::Vec{<:Any, <:FloatingTypes}) = Vec(Intrinsics.fneg(v.data))
Base.FastMath.sub_fast(v::Vec{<:Any, <:FloatingTypes}) = Vec(Intrinsics.fneg(v.data, FASTMATH))
Base.:~(v::Vec{N, T}) where {N, T<:IntegerTypes} = Vec(Intrinsics.xor(v.data, Vec{N, T}(-1).data))
Base.:~(v::Vec{N, Bool}) where {N} = Vec(Intrinsics.xor(v.data, Vec{N, Bool}(true).data))
Base.abs(v::Vec{N, T}) where {N, T} = Vec(vifelse(v < zero(T), -v, v))
Base.:!(v1::Vec{N,Bool}) where {N} = ~v1
Base.inv(v::Vec{N, T}) where {N, T<:FloatingTypes} = one(T) / v

_unsigned(::Type{Float32}) = UInt32
_unsigned(::Type{Float64}) = UInt64
function Base.issubnormal(x::Vec{N, T}) where {N, T<:FloatingTypes}
    y = reinterpret(Vec{N, _unsigned(T)}, x)
    (y & Base.exponent_mask(T) == 0) & (y & Base.significand_mask(T) != 0)
end

@inline Base.signbit(x::Vec{N, <:IntegerTypes}) where {N} = x < 0

@inline Base.leading_ones(x::Vec{<:Any, <:IntegerTypes})  = leading_zeros(~(x))
@inline Base.trailing_ones(x::Vec{<:Any, <:IntegerTypes}) = trailing_zeros(~(x))
@inline Base.count_zeros(x::Vec{<:Any, <:IntegerTypes}) = count_ones(~(x))

@inline Base.isnan(v::Vec{<:Any, <:FloatingTypes}) = v != v
@inline Base.isfinite(v::Vec{<:Any, <:FloatingTypes}) = v - v == zero(v)
@inline Base.isinf(v::Vec{<:Any, <:FloatingTypes}) = !isnan(v) & !isfinite(v)
@inline Base.sign(v1::Vec{N,T}) where {N,T} =
    vifelse(v1 == zero(Vec{N,T}), zero(Vec{N,T}),
            vifelse(v1 < zero(Vec{N,T}), -one(Vec{N,T}), one(Vec{N,T})))

@inline Base.isnan(v::Vec{N, <:IntegerTypes}) where {N} = zero(Vec{N,Bool})
@inline Base.isfinite(v::Vec{N, <:IntegerTypes}) where {N} = one(Vec{N, Bool})
@inline Base.isinf(v::Vec{N, <:IntegerTypes}) where {N} = zero(Vec{N, Bool})

####################
# Binary operators #
####################

const BINARY_OPS = [
    (:(Base.:+)        , IntegerTypes  , Intrinsics.add)
    (:(Base.:-)        , IntegerTypes  , Intrinsics.sub)
    (:(Base.:*)        , IntegerTypes  , Intrinsics.mul)
    (:(Base.div)       , UIntTypes     , Intrinsics.udiv)
    (:(Base.div)       , IntTypes      , Intrinsics.sdiv)
    (:(Base.rem)       , UIntTypes     , Intrinsics.urem)
    (:(Base.rem)       , IntTypes      , Intrinsics.srem)

    (:(add_saturate) , IntTypes  , Intrinsics.sadd_sat)
    (:(add_saturate) , UIntTypes , Intrinsics.uadd_sat)
    (:(sub_saturate) , IntTypes  , Intrinsics.ssub_sat)
    (:(sub_saturate) , UIntTypes , Intrinsics.usub_sat)

    (:(Base.:+)        , FloatingTypes , Intrinsics.fadd)
    (:(Base.:-)        , FloatingTypes , Intrinsics.fsub)
    (:(Base.:*)        , FloatingTypes , Intrinsics.fmul)
    (:(Base.:^)        , FloatingTypes , Intrinsics.pow)
    (:(Base.:/)        , FloatingTypes , Intrinsics.fdiv)
    (:(Base.rem)       , FloatingTypes , Intrinsics.frem)
    (:(Base.min)       , FloatingTypes , Intrinsics.minnum)
    (:(Base.max)       , FloatingTypes , Intrinsics.maxnum)
    (:(Base.copysign)  , FloatingTypes , Intrinsics.copysign)
    (:(Base.:~)        , BIntegerTypes , Intrinsics.xor)
    (:(Base.:&)        , BIntegerTypes , Intrinsics.and)
    (:(Base.:|)        , BIntegerTypes , Intrinsics.or)
    (:(Base.:⊻)        , BIntegerTypes , Intrinsics.xor)

    (:(Base.:(==))   , BIntegerTypes  , Intrinsics.icmp_eq)
    (:(Base.:!=)     , BIntegerTypes  , Intrinsics.icmp_ne)
    (:(Base.:>)      , BIntTypes      , Intrinsics.icmp_sgt)
    (:(Base.:>=)     , BIntTypes      , Intrinsics.icmp_sge)
    (:(Base.:<)      , BIntTypes      , Intrinsics.icmp_slt)
    (:(Base.:<=)     , BIntTypes      , Intrinsics.icmp_sle)
    (:(Base.:>)      , UIntTypes      , Intrinsics.icmp_ugt)
    (:(Base.:>=)     , UIntTypes      , Intrinsics.icmp_uge)
    (:(Base.:<)      , UIntTypes      , Intrinsics.icmp_ult)
    (:(Base.:<=)     , UIntTypes      , Intrinsics.icmp_ule)

    (:(Base.:(==))   , FloatingTypes , Intrinsics.fcmp_oeq)
    (:(Base.:!=)     , FloatingTypes , Intrinsics.fcmp_une)
    (:(Base.:>)      , FloatingTypes , Intrinsics.fcmp_ogt)
    (:(Base.:>=)     , FloatingTypes , Intrinsics.fcmp_oge)
    (:(Base.:<)      , FloatingTypes , Intrinsics.fcmp_olt)
    (:(Base.:<=)     , FloatingTypes , Intrinsics.fcmp_ole)
]

function get_fastmath_function(op)
    if op isa Expr && op.head == Symbol(".") && op.args[1] == :Base &&
        op.args[2].value in keys(Base.FastMath.fast_op)
        return :(Base.FastMath.$(Base.FastMath.fast_op[op.args[2].value]))
    end
    return nothing
end

for (op, constraint, llvmop) in BINARY_OPS
    @eval @inline function $op(x::Vec{N, T}, y::Vec{N, T}) where {N, T <: $constraint}
        Vec($(llvmop)(x.data, y.data))
    end

    # Add a fast math version if applicable
    if (fast_op = get_fastmath_function(op)) !== nothing
        @eval @inline function $(fast_op)(x::Vec{N, T}, y::Vec{N, T}) where {N, T <: $constraint}
            Vec($(llvmop)(x.data, y.data, FASTMATH))
        end
    end
end

# overflow
const OVERFLOW_INTRINSICS = [
    (:(Base.Checked.add_with_overflow) , IntTypes  , Intrinsics.sadd_with_overflow)
    (:(Base.Checked.add_with_overflow) , UIntTypes , Intrinsics.uadd_with_overflow)
    (:(Base.Checked.sub_with_overflow) , IntTypes  , Intrinsics.ssub_with_overflow)
    (:(Base.Checked.sub_with_overflow) , UIntTypes , Intrinsics.usub_with_overflow)
    (:(Base.Checked.mul_with_overflow) , IntTypes  , Intrinsics.smul_with_overflow)
    (:(Base.Checked.mul_with_overflow) , UIntTypes , Intrinsics.umul_with_overflow)
]
for (op, constraint, llvmop) in OVERFLOW_INTRINSICS
    @eval @inline function $op(x::Vec{N, T}, y::Vec{N, T}) where {N, T <: $constraint}
        val, overflows = $(llvmop)(x.data, y.data)
        return Vec(val), Vec(overflows)
    end
end

# max min
@inline Base.max(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
    Vec(vifelse(v1 >= v2, v1, v2))
@inline Base.min(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntegerTypes} =
    Vec(vifelse(v1 >= v2, v2, v1))

# Pow
@inline Base.:^(x::Vec{N,T}, y::IntegerTypes) where {N,T<:FloatingTypes} =
    Vec(Intrinsics.powi(x.data, y))
# Do what Base does for HWNumber:
@inline Base.literal_pow(::typeof(^), x::Vec, ::Val{0}) = one(typeof(x))
@inline Base.literal_pow(::typeof(^), x::Vec, ::Val{1}) = x
@inline Base.literal_pow(::typeof(^), x::Vec, ::Val{2}) = x*x
@inline Base.literal_pow(::typeof(^), x::Vec, ::Val{3}) = x*x*x

# Sign
@inline Base.flipsign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T} =
    vifelse(signbit(v2), -v1, v1)
@inline Base.copysign(v1::Vec{N,T}, v2::Vec{N,T}) where {N,T<:IntTypes} =
    vifelse(signbit(v2), -abs(v1), abs(v1))
_signed(::Type{Float32}) = Int32
_signed(::Type{Float64}) = Int64
@inline Base.signbit(x::Vec{N, T}) where {N, T <:FloatingTypes} =
    signbit(reinterpret(Vec{N, _signed(T)}, x))

# Pointer arithmetic
for op in (:+, :-)
    @eval begin
        # Cast pointer to Int and back
        @inline Base.$op(x::Vec{N,Ptr{T}}, y::Vec{N,Ptr{T}}) where {N,T} =
            convert(Vec{N, Ptr{T}}, ($(op)(convert(Vec{N, Int}, x), convert(Vec{N, Int}, y))))
        @inline Base.$op(x::Vec{N,Ptr{T}}, y::Union{IntegerTypes}) where {N,T} = $(op)(x, Vec{N,Ptr{T}}(y))
        @inline Base.$op(x::IntegerTypes, y::Union{Vec{N,Ptr{T}}}) where {N,T} = $(op)(y, x)

        @inline Base.$op(x::Vec{N,<:IntegerTypes}, y::Ptr{T}) where {N,T} = $(op)(Vec{N,Ptr{T}}(x), Vec{N,Ptr{T}}(y))
        @inline Base.$op(x::Ptr{T}, y::Vec{N,<:IntegerTypes}) where {N,T} = $(op)(y, x)
    end
end

# Bitshifts
# See https://github.com/JuliaLang/julia/blob/7426625b5c07b0d93110293246089a259a0a677d/src/intrinsics.cpp#L1179-L1196
# Shifting with a value larger than the number of bits in the type is undefined behavior
# so set to zero in those cases.
@inline function shl_int(x::Vec{N, T1}, y::Vec{N, T2}) where {N, T1<:IntegerTypes, T2<:IntegerTypes}
    vifelse(y > sizeof(T1) * 8,
        zero(Vec{N, T1}),
        Vec(Intrinsics.shl(x.data, convert(Vec{N,T1}, y).data)))
end

@inline function lshr_int(x::Vec{N, T1}, y::Vec{N, T2}) where {N, T1<:IntegerTypes, T2<:IntegerTypes}
    vifelse(y > sizeof(T1) * 8,
        zero(Vec{N, T1}),
        Vec(Intrinsics.lshr(x.data, convert(Vec{N,T1}, y).data)))
end

@inline function ashr_int(x::Vec{N, T1}, y::Vec{N, T2}) where {N, T1<:IntegerTypes, T2<:IntegerTypes}
    vifelse(y > sizeof(T1) * 8,
            Vec(Intrinsics.ashr(x.data, Vec{N,T1}(sizeof(T1)*8-1).data)),
            Vec(Intrinsics.ashr(x.data, Vec{N,T1}(y).data)))
end

# See https://github.com/JuliaLang/julia/blob/a211abcdfacc05cb93c15774a59ce8961c16dac4/base/int.jl#L422-L435
@inline Base.:>>(x::Vec{N, <:IntTypes}, y::Vec{N, <:UIntTypes}) where {N} =
    ashr_int(x, y)
@inline Base.:>>(x::Vec{N, T1}, y::Vec{N, T2}) where {N, T1<:UIntTypes, T2<:UIntTypes} =
    lshr_int(x, y)
@inline Base.:<<(x::Vec{N, T1}, y::Vec{N, T2}) where {N, T1<:IntegerTypes, T2<:UIntTypes} =
    shl_int(x, y)
@inline Base.:>>>(x::Vec{N, T1}, y::Vec{N, T2}) where {N, T1<:IntegerTypes, T2<:UIntTypes} =
    lshr_int(x, y)

@inline unsigned(v::Vec{<:Any, <:UIntTypes}) = v
@inline unsigned(v::Vec{N, Int32}) where {N} = convert(Vec{N, UInt32}, v)
@inline unsigned(v::Vec{N, Int64}) where {N} = convert(Vec{N, UInt64}, v)

@inline Base.:>>(x::Vec{N, T1}, y::Vec{N, Int}) where {N, T1<:IntegerTypes} =
    vifelse(0 <= y, x >> unsigned(y), x << unsigned(-y))
@inline Base.:<<(x::Vec{N, T1}, y::Vec{N, Int}) where {N, T1<:IntegerTypes} =
    vifelse(0 <= y, x << unsigned(y), x >> unsigned(-y))
@inline Base.:>>>(x::Vec{N, T1}, y::Vec{N, Int}) where {N, T1<:IntegerTypes} =
    vifelse(0 <= y, x >>> unsigned(y), x << unsigned(-y))

for v in (:<<, :>>, :>>>)
    @eval begin
        @inline Base.$v(x::Vec{N,T}, y::ScalarTypes) where {N, T} = $v(x, Vec{N,T}(y))
        @inline Base.$v(x::Vec{N,T}, y::T2) where {N, T<:IntegerTypes, T2<:UIntTypes} = $v(x, Vec{N,T2}(y))
        @inline Base.$v(x::ScalarTypes, y::Vec{N,T}) where {N, T} = $v(Vec{N,T}(x), y)
        @inline Base.$v(x::Vec{N,T1}, y::Vec{N,T2}) where {N, T1<:IntegerTypes, T2<:IntegerTypes} =
            $v(x, convert(Vec{N, Int}, y))
    end
end

# Vectorize binary functions
for (op, constraint) in [BINARY_OPS;
        (:(Base.flipsign) , ScalarTypes)
        (:(Base.copysign) , ScalarTypes)
        (:(Base.signbit)  , ScalarTypes)
        (:(Base.min)      , IntegerTypes)
        (:(Base.max)      , IntegerTypes)
        (:(Base.:<<)      , IntegerTypes)
        (:(Base.:>>)      , IntegerTypes)
        (:(Base.:>>>)     , IntegerTypes)
        (:(Base.Checked.add_with_overflow) , IntTypes)
        (:(Base.Checked.add_with_overflow) , UIntTypes)
        (:(Base.Checked.sub_with_overflow) , IntTypes)
        (:(Base.Checked.sub_with_overflow) , UIntTypes)
        (:(Base.Checked.mul_with_overflow) , IntTypes)
        (:(Base.Checked.mul_with_overflow) , UIntTypes)
    ]
    ops = [op]
    if (fast_op = get_fastmath_function(op)) !== nothing
        push!(ops, fast_op)
    end
    for op in ops
        @eval @inline function $op(x::T2, y::Vec{N, T}) where {N, T2<:ScalarTypes, T <: $constraint}
            $op(Vec{N, T}(x), y)
        end
        @eval @inline function $op(x::Vec{N, T}, y::T2) where {N, T2 <:ScalarTypes, T <: $constraint}
            $op(x, Vec{N, T}(y))
        end
    end
end

#####################
# Ternary operators #
#####################

@inline vifelse(v::Bool, v1::Vec{N, T}, v2::Vec{N, T}) where {N, T} = ifelse(v, v1, v2)
@inline vifelse(v::Bool, v1::Vec{N, T}, v2::ScalarTypes) where {N, T} = ifelse(v, v1, Vec{N,T}(v2))
@inline vifelse(v::Bool, v1::ScalarTypes, v2::Vec{N, T}) where {N, T} = ifelse(v, Vec{N,T}(v1), v2)

@inline vifelse(v::Bool, v1::T, v2::T) where {T} = ifelse(v, v1, v2)
@inline vifelse(v::Vec{N, Bool}, v1::Vec{N, T}, v2::Vec{N, T}) where {N, T} =
    Vec(Intrinsics.select(v.data, v1.data, v2.data))
@inline vifelse(v::Vec{N, Bool}, v1::T2, v2::Vec{N, T}) where {N, T, T2 <:ScalarTypes} = vifelse(v, Vec{N, T}(v1), v2)
@inline vifelse(v::Vec{N, Bool}, v1::Vec{N, T}, v2::T2) where {N, T, T2 <:ScalarTypes} = vifelse(v, v1, Vec{N, T}(v2))

# fma, muladd and vectorization of these
for (op, llvmop) in [(:fma, Intrinsics.fma), (:muladd, Intrinsics.fmuladd)]
    @eval begin
        @inline Base.$op(a::Vec{N, T}, b::Vec{N, T}, c::Vec{N, T}) where {N,T<:FloatingTypes} =
            Vec($llvmop(a.data, b.data, c.data))
        @inline Base.$op(s1::ScalarTypes, v2::Vec{N,T}, v3::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), v2, v3)
        @inline Base.$op(v1::Vec{N,T}, s2::ScalarTypes, v3::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(v1, Vec{N,T}(s2), v3)
        @inline Base.$op(s1::ScalarTypes, s2::ScalarTypes, v3::Vec{N,T}) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), Vec{N,T}(s2), v3)
        @inline Base.$op(v1::Vec{N,T}, v2::Vec{N,T}, s3::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(v1, v2, Vec{N,T}(s3))
        @inline Base.$op(s1::ScalarTypes, v2::Vec{N,T}, s3::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(Vec{N,T}(s1), v2, Vec{N,T}(s3))
        @inline Base.$op(v1::Vec{N,T}, s2::ScalarTypes, s3::ScalarTypes) where {N,T<:FloatingTypes} =
            $op(v1, Vec{N,T}(s2), Vec{N,T}(s3))
    end
end

##############
# Reductions #
##############
const HORZ_REDUCTION_OPS = [
    (&   , Union{IntegerTypes, Bool}  , Intrinsics.reduce_and)
    (|   , Union{IntegerTypes, Bool}  , Intrinsics.reduce_or)
    (max , IntTypes      , Intrinsics.reduce_smax)
    (max , UIntTypes     , Intrinsics.reduce_umax)
    (max , FloatingTypes , Intrinsics.reduce_fmax)
    (min , IntTypes      , Intrinsics.reduce_smin)
    (min , UIntTypes     , Intrinsics.reduce_umin)
    (min , FloatingTypes , Intrinsics.reduce_fmin)
    (+   , Union{IntegerTypes, Bool}  , Intrinsics.reduce_add)
    (*   , IntegerTypes  , Intrinsics.reduce_mul)
    (+   , FloatingTypes , Intrinsics.reduce_fadd)
    (*   , FloatingTypes , Intrinsics.reduce_fmul)
]

for (op, constraint, llvmop) in HORZ_REDUCTION_OPS
    @eval @inline Base.reduce(::typeof($op), x::Vec{<:Any, <:$constraint}) =
        $(llvmop)(x.data)
end
Base.reduce(F::Any, v::Vec) = error("reduction not defined for SIMD.Vec on $F")

@inline Base.all(v::Vec{<:Any,Bool}) = reduce(&, v)
@inline Base.any(v::Vec{<:Any,Bool}) = reduce(|, v)
@inline Base.maximum(v::Vec) = reduce(max, v)
@inline Base.minimum(v::Vec) = reduce(min, v)
@inline Base.prod(v::Vec) = reduce(*, v)
@inline Base.sum(v::Vec) = reduce(+, v)

############
# Shuffles #
############

@inline function shufflevector(x::Vec{N, T}, ::Val{I}) where {N, T, I}
    Vec(Intrinsics.shufflevector(x.data, Val(I)))
end
@inline function shufflevector(x::Vec{N, T}, y::Vec{N, T}, ::Val{I}) where {N, T, I}
    Vec(Intrinsics.shufflevector(x.data, y.data, Val(I)))
end


text:
#   Unit "cv.jl" of the PosDefManifoldML Package for Julia language

#   MIT License
#   Copyright (c) 2019-2021
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home
#   Saloni Jain, Indian Institute of Technology, Kharagpur, India

# ? CONTENTS :
#   This unit implements cross-validation procedures for estimating
#   accuracy of all machine learning models.

"""
```julia
struct CVacc
    cvType    :: String
    scoring   :: Union{String, Nothing}
    modelType :: Union{String, Nothing}
    cnfs      :: Union{Vector{Matrix{T}}, Nothing} where T<:Real
    avgCnf    :: Union{Matrix{T}, Nothing} where T<:Real
    accs      :: Union{Vector{T}, Nothing} where T<:Real
    avgAcc    :: Union{Real, Nothing}
    stdAcc    :: Union{Real, Nothing}
end
```

A call to [`cvAcc`](@ref) results in an instance of this structure.
Fields:

`.cvTpe` is the type of cross-validation technique, given as a string
(e.g., "10-kfold")

`.scoring` is the type of accuracy that is computed, given as a string.
This has been passed as argument to [`cvAcc`](@ref).
Currently *accuracy* and *balanced accuracy* are supported.

`.modelType` is type of the machine learning used for performing the
cross-validation, given as a string.

`.cnfs` is a vector of matrices holding the *confusion matrices*
obtained at each fold of the cross-validation.

`.avgCnf` is the *average confusion matrix* across the folds of the
cross-validation.

`.accs` is a vector of real numbers holding the *accuracies* obtained
at each fold of the cross-validation.

`.avgAcc` is the *average accuracy* across the folds of the
cross-validation.

`.stdAcc` is the *standard deviation of the accuracy* across the folds of the
cross-validation.

"""
struct CVacc
    cvType    :: String
    scoring   :: Union{String, Nothing}
    modelType :: Union{String, Nothing}
    cnfs      :: Union{Vector{Matrix{T}}, Nothing} where T<:Real
    avgCnf    :: Union{Matrix{T}, Nothing} where T<:Real
    accs      :: Union{Vector{T}, Nothing} where T<:Real
    avgAcc    :: Union{Real, Nothing}
    stdAcc    :: Union{Real, Nothing}
end

"""
```julia
CVacc(s::String) =
     CVacc(s, nothing, nothing, nothing,
           nothing, nothing, nothing, nothing)
```

Construct an instance of the CVacc structure giving only the `.cvtype`
field. All other fields are filled with `nothing`. This is useful to construct
manually cvAcc objects.
"""
CVacc(s::String)=CVacc(s, nothing, nothing, nothing, nothing, nothing, nothing, nothing)

"""
```julia
function cvAcc(model    :: MLmodel,
               𝐏Tr     :: ℍVector,
               yTr      :: IntVector;
           nFolds       :: Int      = min(10, length(yTr)÷3),
           scoring      :: Symbol   = :b,
           shuffle      :: Bool     = false,
           verbose      :: Bool     = true,
           outModels    :: Bool     = false,
           ⏩           :: Bool     = true,
           fitArgs...)
```
Cross-validation accuracy for a machine learning `model`:
given an ℍVector `𝐏Tr` holding *k* Hermitian matrices,
an [IntVector](@ref) `yTr` holding the *k* labels for these matrices and
the number of folds `nFolds`,
return a [`CVacc`](@ref) structure.

**optional keyword arguments**

`nFolds` by default is set to the minimum between 10 and the number
of observation ÷ 3 (integer division).

If `scoring`=:b (default) the **balanced accuracy** is computed.
Any other value will make the function returning the regular **accuracy**.
Balanced accuracy is to be preferred for unbalanced classes.
For balanced classes the balanced accuracy reduces to the
regular accuracy, therefore there is no point in using regular accuracy
if not to avoid a few unnecessary computations when the class are balanced.

For the meaning of the `shuffle` argument (false by default),
see function [`cvSetup`](@ref), to which this argument is passed.

If `verbose` is true (default), information is printed in the REPL.
This option is included to allow repeated calls to this function
without crowding the REPL.

If `outModels` is true return a 2-tuple holding a [`CVacc`](@ref) structure
and a `nFolds`-vector of the model fitted for each fold,
otherwise (default), return only a [`CVacc`](@ref) structure.

If `⏩` the folds and some other computations are multi-threaded.
It is true by default. Set it to false if there are problems in running
this function.

`fitArgs` are optional keyword arguments that are passed to the
[`fit`](@ref) function called for each fold of the cross-validation.
For each machine learning model, all optional keyword arguments of
their fit method are elegible to be passed here, however,
the arguments listed in the following table for each model should not be passed.
Note that if they are passed, they will be disabled:

| MDM/MDMF |   ENLR    |    SVM    |
|:--------:|:---------:|:---------:|
| `verbose`| `verbose` | `verbose` |
|  `⏩`   | `⏩`      | `⏩`      |
|          | `meanISR` | `meanISR` |
|          | `meanInit`| `meanInit`|
|          | `fitType` |
|          | `offsets` |
|          | `lambda`  |
|          | `folds`   |

Also, if you pass a `w` (weights for tangent space projection) argument,
do not pass a vector of weights, just pass a symbol, e.g., `w=:b`
for balancing weights.

**See**: [notation & nomenclature](@ref), [the ℍVector type](@ref).

**See also**: [`fit`](@ref), [`predict`](@ref).

**Examples**
```julia
using PosDefManifoldML, PosDefManifold

# generate some data
PTr, PTe, yTr, yTe=gen2ClassData(10, 30, 40, 60, 80)

# perform 10-fold cross-validation using the minimum distance to mean classifier
cv=cvAcc(MDM(Fisher), PTr, yTr)

# ...using the lasso logistic regression classifier
cv=cvAcc(ENLR(Fisher), PTr, yTr)

# ...using the support-vector machine classifier
cv=cvAcc(SVM(Fisher), PTr, yTr)

# ...With a Polynomial kernel of order 3 (default)
cv=cvAcc(SVM(Fisher), PTr, yTr; kernel=kernel.Polynomial)

# perform 8-fold cross-validation instead
# (and see that you can go pretty fast if your PC has 8 threads)
cv=cvAcc(SVM(Fisher), PTr, yTr; nFolds=8)

# ...balance the weights for tangent space projection
cv=cvAcc(ENLR(Fisher), PTr, yTr; nFolds=8, w=:b)

# perform another cross-validation shuffling the folds
cv=cvAcc(ENLR(Fisher), PTr, yTr; shuffle=true, nFolds=8, w=:b)

```
"""
function cvAcc(model    :: MLmodel,
               𝐏Tr     :: ℍVector,
               yTr      :: IntVector;
           nFolds       :: Int      = min(10, length(yTr)÷3),
           scoring      :: Symbol   = :b,
           shuffle      :: Bool     = false,
           verbose      :: Bool     = true,
           outModels    :: Bool     = false,
           ⏩           :: Bool     = true,
           fitArgs...)

    ⌚ = now()
    verbose && println(greyFont, "\nPerforming $(nFolds)-fold cross-validation...")

    z  = length(unique(yTr))            # number of classes
    𝐐  = [ℍ[] for i=1:z]               # data arranged by class
    for j=1:length(𝐏Tr) push!(𝐐[yTr[j]], 𝐏Tr[j]) end

    # pre-allocated memory
    𝐐Tr = [ℍ[] for f=1:nFolds]                 # training data in 1 vector per folds
    zTr = [Int64[] for f=1:nFolds]              # training labels in 1 vector per fold
    𝐐Te = [[ℍ[] for i=1:z] for f=1:nFolds]     # testing data arranged by classes per fold
    CM  = [zeros(Float64, z, z) for f=1:nFolds] # confusion matrices per fold
    s   = Vector{Float64}(undef, nFolds)        # accuracy scores per fold
    pl  = [[Int[] for i=1:z] for f=1:nFolds]    # predicted labels per fold
    indTr = [[[]] for i=1:z]                    # indeces for training sets per fold
    indTe = [[[]] for i=1:z]                    # indeces for test sets per fold
    ℳ=Vector{MLmodel}(undef, nFolds)            # ML models

    # get indeces for all CVs (separated for each class)
    for i=1:z indTr[i], indTe[i] = cvSetup(length(𝐐[i]), nFolds; shuffle=shuffle) end

    fitArgs✔=()
    # make sure the user doesn't pass arguments that skrew up the cv
    if model isa ENLRmodel
        fitArgs✔=_rmArgs((:meanISR, :meanInit, :fitType, :verbose, :⏩,
                         :offsets, :lambda, :folds); fitArgs...)
        # overwrite the `alpha` value in `model` if the user has passed keyword `alpha`
        if (a=_getArgValue(:alpha; fitArgs...)) ≠ nothing model.alpha=a end
    end

    if model isa SVMmodel
        # make sure the user doesn't pass arguments that skrew up the cv
        fitArgs✔=_rmArgs((:meanISR, :meanInit, :verbose, :⏩); fitArgs...)
        # overwrite the `svmType` and `kernel` values in `model` if the user has passed then as kwargs
        if (a=_getArgValue(:svmType; fitArgs...)) ≠ nothing model.svmType=a end
        if (a=_getArgValue(:kernel;  fitArgs...)) ≠ nothing model.kernel =a end
    end

    # get means initializations. Be careful here with the behavior for other models
    # NB for the MDM model the initialization is a vector of means for each class,
    # for the TSmodels it is the mean of such means.
    # This is a quick approximation since the initialization is not critical,
    # but it hastens the computation time since itera. alg. require less iters.
    #=
    if      model.metric in (Fisher, logdet0)
                M0=means(logEuclidean, 𝐐; ⏩=⏩)
                if model isa TSmodel M0=mean(logEuclidean, M0; ⏩=⏩) end
    elseif  model.metric == Wasserstein
                M0=ℍVector([generalizedMean(𝐐[i], 0.5; ⏩=⏩) for i=1:length(𝐐)])
                if model isa Tsmodel M0=generalizedMean(M0, 0.5; ⏩=⏩) end
    else    M0=nothing;
    end
    =#

    # perform cv
    function fold(f::Int)
        @static if VERSION >= v"1.3" print(defaultFont, rand(dice), " ") end # print a random dice in the REPL

        # get testing data for current fold
        for i=1:z 𝐐Te[f][i] = [𝐐[i][j] for j ∈ indTe[i][f]] end

        # get training labels for current fold
        for i=1:z, j ∈ indTr[i][f] push!(zTr[f], Int64(i)) end

        # get training data for current fold
        for i=1:z, j ∈ indTr[i][f] push!(𝐐Tr[f], 𝐐[i][j]) end

        # fit machine learning model
        ℳ[f]=fit(model, 𝐐Tr[f], zTr[f];
                  #meanInit=M0,
                  verbose=false,
                  ⏩=false,
                  fitArgs✔...)

        # predict labels for current fold
        for i=1:z pl[f][i]=predict(ℳ[f], 𝐐Te[f][i], :l; verbose=false, ⏩=false) end

        # compute confusion matrix for current fold
        for i=1:z, s=1:length(pl[f][i]) CM[f][i, pl[f][i][s]]+=1. end

        # compute balanced accuracy or accuracy for current CV
        sumCM=sum(CM[f])
        scoring == :b ? s[f] = 𝚺(CM[f][i, i]/𝚺(CM[f][i, :]) for i=1:z) / z :
                        s[f] = 𝚺(CM[f][i, i] for i=1:z)/ sumCM

        CM[f]/=sumCM # confusion matrices in proportions
    end

    # This actually runs the cross-validation
    ⏩ ? (@threads for f=1:nFolds fold(f) end) : (for f=1:nFolds fold(f) end)
    verbose && println(greyFont, "\nDone in ", defaultFont, now()-⌚)

    # compute mean and sd (balanced) accuracy
    avg=mean(s);
    std=stdm(s, avg);
    scoStr = scoring == :b ? "balanced accuracy" : "accuracy"

    # create cv struct
    cv=CVacc("$nFolds-fold", scoStr, _modelStr(model), CM, mean(CM), s, avg, std)
    return outModels ? (cv, ℳ) : cv
end

"""
```julia
function cvSetup(k       :: Int,
                 nCV     :: Int;
                 shuffle :: Bool = false)
```

Given `k` elements and a parameter `nCV`, a nCV-fold cross-validation
is obtained defining `nCV` permutations of *k* elements
in *nTest=k÷nCV* (integer division) elements for the test and
*k-nTest* elements for the training,
in such a way that each element is represented in only one permutation.

Said differently, given a length `k` and the number of desired cross-validations
`nCV`, this function generates indices from the sequence of natural numbers
*1,..,k* to obtain all nCV-fold cross-validation sets.
Specifically, it generates `nCV` vectors of indices for generating test sets
and `nCV` vectors of indices for geerating training sets.

If optional keyword argument `shuffle` is true,
the sequence of natural numbers *1,..,k* is shuffled before
running the function, thus in this case two successive runs of this function
will give different cross-validation sets, hence different accuracy scores.
By default `shuffle` is false, so as to allow exactly the same result
in successive runs.
Note that no random initialization for the shuffling is provided, so as to
allow the replication of the same random sequences starting again
the random generation from scratch.

This function is used in [`cvAcc`](@ref). It constitutes the fundamental
basis to implement customized cross-validation procedures.

Return the 2-tuple with:

- A vector of `nCV` vectors holding the indices for the training sets,
- A vector of `nCV` vectors holding the indices for the corresponding test sets.

**Examples**
```julia
using PosDefManifoldML, PosDefManifold

cvSetup(10, 2)
# return:
# (Array{Int64,1}[[6, 7, 8, 9, 10], [1, 2, 3, 4, 5]],
#  Array{Int64,1}[[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]])

cvSetup(10, 2, shuffle=true)
# return:
# (Array{Int64,1}[[5, 4, 6, 1, 9], [3, 7, 8, 2, 10]],
#  Array{Int64,1}[[3, 7, 8, 2, 10], [5, 4, 6, 1, 9]])

cvSetup(10, 3)
# return:
# (Array{Int64,1}[[4, 5, 6, 7, 8, 9, 10], [1, 2, 3, 7, 8, 9, 10], [1, 2, 3, 4, 5, 6]],
#  Array{Int64,1}[[1, 2, 3], [4, 5, 6], [7, 8, 9, 10]])

```

"""
function cvSetup(k       :: Int,
                 nCV     :: Int;
                 shuffle :: Bool = false)

    if nCV == 1 @error 📌*", cvSetup function: The number of cross-validation must be bigger than one" end
    nTest = k÷nCV # nTrain = k-nTest
    #rng = MersenneTwister(1900)
    shuffle ? a=shuffle!( Vector(1:k)) : a=Vector(1:k)
    indTrain = [IntVector(undef, 0) for i=1:nCV]
    indTest  = [IntVector(undef, 0) for i=1:nCV]

    # vectors of indices for test and training sets
    j=1
    for i=1:nCV-1
        indTest[i]=a[j:j+nTest-1]
        for g=j+nTest:length(a) push!(indTrain[i], a[g]) end
        for l=i+1:nCV, g=j:j+nTest-1 push!(indTrain[l], a[g]) end
        j+=nTest
    end
    indTest[nCV]=a[j:end]
    return indTrain, indTest
end

# ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, cv::CVacc)
                            println(io, titleFont, "\n◕ Cross-Validation Accuracy")
                            println(io, separatorFont, "⭒  ⭒    ⭒       ⭒         ⭒", defaultFont)
                            println(io, separatorFont, ".cvType   :", defaultFont," $(cv.cvType)")
    cv.scoring   ≠ nothing && println(io, separatorFont, ".scoring  :", defaultFont," $(cv.scoring)")
    cv.modelType ≠ nothing && println(io, separatorFont, ".modelType:", defaultFont," $(cv.modelType)")
    cv.cnfs      ≠ nothing && println(io, separatorFont, ".cnfs      ", defaultFont,"(confusion mat. per fold)")
    cv.avgCnf    ≠ nothing && println(io, separatorFont, ".avgCnf    ", defaultFont,"(average confusion mat. )")
    cv.accs      ≠ nothing && println(io, separatorFont, ".accs      ", defaultFont,"(accuracies per fold    )")
    cv.avgAcc    ≠ nothing && println(io, separatorFont, ".avgAcc   :", defaultFont," $(round(cv.avgAcc; digits=3)) (average accuracy)")
    cv.stdAcc    ≠ nothing && println(io, separatorFont, ".stdAcc   :", defaultFont," $(round(cv.stdAcc; digits=3)) (st. dev accuracy)")
end

ext: jl


text:
@testset "PPW" begin
    k = 2
    U = [0, 0, 0, 0, 1]
    result = Dirichlet_sampling(k, U, 20)
    @test result[4] ≈ 0 atol = 1e-12
    @test result[6](0.5,π/4) ≈ 0 atol = 1e-12
end

@testset "EPW" begin
    k = 1
    U = [0, 0, 0, 0, 1]
    @testset "Uniform" begin
        result = Dirichlet_sampling(k, U, 40; smpl_type=uniform_sampling)
        @test result[4] ≈ 0 atol = 1e-8
        @test result[6](0.5,π/4) ≈ 0 atol = 1e-12
    end
    @testset "Sobol" begin
        result = Dirichlet_sampling(k, U, 40; smpl_type=sobol_sampling)
        @test result[4] ≈ 0 atol = 1e-8
        @test result[6](0.5,π/4) ≈ 0 atol = 1e-12
    end
    @testset "Random" begin
        result = Dirichlet_sampling(k, U, 40; smpl_type=random_sampling)
        @test result[4] ≈ 0 atol = 1e-7
        @test result[6](0.5,π/4) ≈ 0 atol = 1e-11
    end
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
module LossFunctions

export mse, cross_entropy

using Statistics

"""
    mse(x, y)

Computes the mean-squared error between `x` and `y`, returning a scalar.

For the multi-sample case, takes the error along columns (dimensions), and
averages along rows (count).

- `x`: Either (1) a vector containing a single prediction, or (2) a matrix where
  rows index predictions in a batch and columns indicate dimensions of the
  prediction.
- `y`: Either (1) a vector containing a single label, or (2) a matrix where rows
  index labels in a batch and columns indicate dimensions of the label.

`x` and `y` should have the same dimensionality.
"""
mse(x, y) = mean(sum((x .- y) .^2, dims=2), dims=1)[1]

"""
    cross_entropy(x, y)

Computes the cross-entropy error between `x` and `y`, returning a scalar.

For the multi-sample case, takes the error along columns (dimensions), and
averages along rows (count).

- `x`: Either (1) a vector containing a single prediction, or (2) a matrix where
  rows index predictions in a batch and columns indicate dimensions of the
  prediction.
- `y`: Either (1) a vector containing a single label, or (2) a matrix where rows
  index labels in a batch and columns indicate dimensions of the label.

`x` and `y` should have the same dimensionality.
"""
cross_entropy(x, y) = mean((-1) * sum(y .* log.(x), dims=2), dims=1)[1]

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
export AnyScope, UnionScope, NodeScope, ProcessScope, ExactScope

abstract type AbstractScope end

"Default scope that is unconstrained."
struct AnyScope <: AbstractScope end

"Union of two or more scopes."
struct UnionScope <: AbstractScope
    scopes::NTuple{N,<:AbstractScope} where N
end
UnionScope(scopes...) = UnionScope((scopes...,))
UnionScope(scopes::Vector{<:AbstractScope}) = UnionScope((scopes...,))
UnionScope(s::AbstractScope) = UnionScope((s,))
UnionScope() = throw(ArgumentError("Cannot construct empty UnionScope"))

"Scoped to the same physical node."
struct NodeScope <: AbstractScope
    uuid::UUID
end
NodeScope() = NodeScope(SYSTEM_UUIDS[myid()])

"Scoped to the same OS process."
struct ProcessScope <: AbstractScope
    parent::NodeScope
    wid::Int
end
function ProcessScope(wid)
    if wid == myid()
        ProcessScope(NodeScope(), wid)
    else
        ProcessScope(NodeScope(system_uuid(wid)), wid)
    end
end
ProcessScope() = ProcessScope(myid)

"Scoped to a specific processor."
struct ExactScope <: AbstractScope
    parent::ProcessScope
    processor::Processor
end
ExactScope(proc) = ExactScope(ProcessScope(get_parent(proc).pid), proc)

"Indicates that the applied scopes `x` and `y` are incompatible."
struct InvalidScope <: AbstractScope
    x::AbstractScope
    y::AbstractScope
end

# Comparisons and constraint checking

constrain(x, y) = x < y ? constrain(y, x) : throw(MethodError(constrain, x, y))

Base.isless(::AnyScope, ::AnyScope) = false
Base.isless(::AnyScope, x) = false
Base.isless(x, ::AnyScope) = true
constrain(::AnyScope, ::AnyScope) = AnyScope()
constrain(::AnyScope, y) = y

Base.isless(::UnionScope, ::UnionScope) = false
Base.isless(::UnionScope, ::AnyScope) = true
function constrain(x::UnionScope, y::UnionScope)
    zs = Vector{AbstractScope}()
    for xs in x.scopes
        for ys in y.scopes
            scope = constrain(xs, ys)
            scope isa InvalidScope && continue
            push!(zs, scope)
        end
    end
    isempty(zs) && return InvalidScope(x, y)
    UnionScope(zs)
end
constrain(x::UnionScope, y) = constrain(x, UnionScope((y,)))

Base.isless(::NodeScope, ::NodeScope) = false
Base.isless(::NodeScope, ::AnyScope) = true
constrain(x::NodeScope, y::NodeScope) =
    x == y ? y : InvalidScope(x, y)

Base.isless(::ProcessScope, ::ProcessScope) = false
Base.isless(::ProcessScope, ::NodeScope) = true
Base.isless(::ProcessScope, ::AnyScope) = true
constrain(x::ProcessScope, y::ProcessScope) =
    x == y ? y : InvalidScope(x, y)
constrain(x::NodeScope, y::ProcessScope) =
    x == y.parent ? y : InvalidScope(x, y)

Base.isless(::ExactScope, ::ExactScope) = false
Base.isless(::ExactScope, ::ProcessScope) = true
Base.isless(::ExactScope, ::NodeScope) = true
Base.isless(::ExactScope, ::AnyScope) = true
constrain(x::ExactScope, y::ExactScope) =
    x == y ? y : InvalidScope(x, y)
constrain(x::ProcessScope, y::ExactScope) =
    x == y.parent ? y : InvalidScope(x, y)
constrain(x::NodeScope, y::ExactScope) =
    x == y.parent.parent ? y : InvalidScope(x, y)


text:
println("========loading and compiling========")
include("try_gpu_accel.jl")

DEFAULT_N = 8
DEFAULT_TF = 32
println("\n========start testing========")

println("\n--------scaling n on cpu--------")
i=2
while i <= 10
    println("n=$(i); tf=$(DEFAULT_TF)")
    for t in 1:3
        r = nothing; GC.gc(true); CUDA.reclaim();
        @time r = anneal_spin_glass(i, DEFAULT_TF);
    end
    if i >= 8
        global i+=1;
    else
        global i+=2;
    end
end

println("\n--------scaling tf on cpu--------")
tf=4
while tf<=128
    println("n=$(DEFAULT_N); tf=$(tf)")
    for t in 1:3
        r = nothing; GC.gc(true); CUDA.reclaim();
        @time r = anneal_spin_glass(DEFAULT_N, tf);
    end
    global tf*=2;
end

println("\n--------scaling n on gpu--------")
i=2
while i<=10
    println("n=$(i); tf=$(DEFAULT_TF)")
    for t in 1:3
        r = nothing; GC.gc(true); CUDA.reclaim();
        @time r = anneal_spin_glass_gpu(i, DEFAULT_TF);
    end
    if i >= 8
        global i+=1;
    else
        global i+=2;
    end
end

println("\n--------scaling tf on cpu--------")
tf=4
while tf<=128
    println("n=$(DEFAULT_N); tf=$(tf)")
    for t in 1:3
        r = nothing; GC.gc(true); CUDA.reclaim();
        @time r = anneal_spin_glass(DEFAULT_N, tf);
    end
    global tf*=2;
end

max_stars_repo_name: naezzell/accelqat

max_issues_repo_name: naezzell/accelqat

max_forks_repo_name: naezzell/accelqat


text:
abstract type CircuitParams{T<:Real} <: AbstractParams end

"""
    Impedance{T<:Real} <: CircuitParams{T}
"""
mutable struct Impedance{T<:Real} <: CircuitParams{T}
    data::Complex{T}
end
Impedance(z::T) where {T<:Real} = Impedance(complex(z))

Impedance(zd::AbstractVector{Complex{T}}) where {T<:Real} =
    [Impedance(zd_) for zd_ in zd]
Impedance(zd::AbstractVector{T}) where {T<:Real} = Impedance(Complex{T}.(zd))

"""
    Admittance{T<:Real} <: CircuitParams{T}
"""
mutable struct Admittance{T<:Real} <: CircuitParams{Real}
    data::Complex{T}
end
Admittance(y::T) where {T<:Real} = Admittance(complex(y))

Admittance(yd::AbstractVector{Complex{T}}) where {T<:Real} =
    [Admittance(yd_) for yd_ in yd]
Admittance(yd::AbstractVector{T}) where {T<:Real} = Admittance(Complex{T}.(yd))

function show(io::IO, params::CircuitParams)
    write(io, "$(typeof(params)): $(params.data)")
end


text:
#=
P07 (**) Flatten a nested list structure.
Transform a list, possibly holding lists as elements into a `flat' list by replacing each list with its elements (recursively).

Example:
* (my-flatten '(a (b (c d) e)))
(A B C D E)

Hint: Use the predefined functions list and append.
=#

function flatten(coll, flattened)
    if isempty(coll)
        flattened
    else
        first, rest = coll[1], coll[2:end]
        if typeof(first) == Symbol
            flatten(rest, vcat(flattened, first))
        else
            vcat(flatten(first, flattened), flatten(rest, []))
        end
    end
end

function my_flatten(coll)
    flatten(coll, [])
end

@assert my_flatten((:a, (:b, (:c, :d), :e))) == [:a, :b, :c, :d, :e]
@assert my_flatten((:a, (:b, (:c,), (:d, :e, :f)), :g)) == [:a, :b, :c, :d, :e, :f, :g]
@assert my_flatten(()) == []

println("Tests passed: JL-07.jl")


text:
using Base: Float64
## Q2
days = 365
days_float = convert(Float64, days)

@assert days == 365
@assert days_float == 365.0

## Q3

convert(Int64, "1")
parse(Int64, "1")

max_stars_repo_path: Exercise Solutions/1-1 Variables and Printing.jl


text:
module LibHello
  using CxxWrap
  using libhello_jll
　const libname = libhello_jll.libhello_path
  @wrapmodule(libname)

  function __init__()
    @initcxx
  end

　export greet, get_array

end


text:
_get_iterable(global_fixtures, fx::AbstractGlobalFixture) = global_fixtures[fx]
_get_iterable(global_fixtures, fx::ConstantFixture) = setup(fx)
function _get_iterable(global_fixtures, fx::LocalFixture)
    # Workaround for Julia issue #29193.
    # If it is fixed, the last line should work for the empty `params` case as well.
    params = parameters(fx)
    if isempty(params)
        [()]
    else
        rowmajor_product([_get_iterable(global_fixtures, param) for param in params]...)
    end
end

get_iterable(global_fixtures) = fx -> _get_iterable(global_fixtures, fx)

struct DelayedTeardownValue
    lval :: LabeledValue
    rff :: Union{Nothing, RunningFixtureFactory}
    subvalues :: Array{DelayedTeardownValue, 1}
end

function setup(fx::LocalFixture, lvals)
    to_release = DelayedTeardownValue[]
    processed_args = []
    for (p, lval) in zip(parameters(fx), lvals)
        if typeof(p) == LocalFixture
            dval = setup(p, lval)
            push!(to_release, dval)
            lval = dval.lval
        end
        push!(processed_args, lval.value)
    end
    lval, rff = setup(fx.ff, processed_args)
    DelayedTeardownValue(lval, rff, to_release)
end

function release(val::DelayedTeardownValue)
    if !(val.rff === nothing)
        for v in val.subvalues
            release(v)
        end
        teardown(val.rff)
    end
end

unwrap_value(val::DelayedTeardownValue) = unwrap_value(val.lval)
unwrap_label(val::DelayedTeardownValue) = unwrap_label(val.lval)

instantiate(fx::LocalFixture, lval) = setup(fx, lval)
instantiate(fx, lval) = DelayedTeardownValue(lval, nothing, DelayedTeardownValue[])

function instantiate_global(global_fixtures, fx::GlobalFixture)

    for_teardown = RunningFixtureFactory[]

    iterables = Array{LabeledValue, 1}[_get_iterable(global_fixtures, p) for p in parameters(fx)]

    all_lvals = LabeledValue[]
    for lvals in rowmajor_product(iterables...)
        args = map(unwrap_value, lvals)
        lval, rff = setup(fx, args)

        push!(all_lvals, lval)
        if instant_teardown(rff)
            teardown(rff)
        else
            push!(for_teardown, rff)
        end
    end
    all_lvals, for_teardown
end

function run_testcases(run_options, tcs, doctest)

    global_fixtures = Dict{AbstractGlobalFixture, Array{LabeledValue, 1}}()
    gi = get_iterable(global_fixtures)
    for_teardown = DefaultDict{Int, Array{RunningFixtureFactory, 1}}(() -> RunningFixtureFactory[])

    progress = ProgressReporter(run_options[:verbosity], doctest)

    progress_start!(progress)

    capture_output = run_options[:capture_output] :: Bool
    fails_num = 0
    max_fails = run_options[:max_fails] :: Int
    max_fails_reached = false
    color_scheme = run_options[:report_color_scheme] :: ReportColorScheme

    for (i, entry) in enumerate(tcs)

        tcinfo, tc = entry

        for fx in dependencies(tc)
            if !haskey(global_fixtures, fx)
                if isa(fx, RunOptionsFixture)
                    global_fixtures[fx] = [LabeledValue(run_options, "run_options")]
                else
                    lvals, ftd = instantiate_global(global_fixtures, fx)
                    global_fixtures[fx] = lvals

                    if length(ftd) > 0
                        last_usage_idx = findlast(tcs) do entry
                            _, tcc = entry
                            fx in dependencies(tcc)
                        end
                        append!(for_teardown[last_usage_idx], ftd)
                    end
                end
            end
        end

        fixture_iterables = map(gi, parameters(tc))
        iterable_permutations = rowmajor_product(fixture_iterables...)

        for lvals in iterable_permutations
            dvals = map(instantiate, parameters(tc), lvals)
            args = map(unwrap_value, dvals)
            labels = map(unwrap_label, dvals)
            progress_start_testcase!(progress, tcinfo, labels)
            outcome = run_testcase(tc, args, capture_output)
            map(release, dvals)
            progress_finish_testcase!(progress, tcinfo, labels, outcome, color_scheme)

            if is_failed(outcome)
                fails_num += 1
            end

            if max_fails > 0 && fails_num == max_fails
                max_fails_reached = true
                break
            end
        end

        if haskey(for_teardown, i)
            map(teardown, for_teardown[i])
            delete!(for_teardown, i)
        end

        if max_fails_reached
            # Call all remaining teardowns
            for (i, rffs) in for_teardown
                map(teardown, rffs)
            end
            break
        end
    end

    progress_finish!(progress, color_scheme)

    fails_num == 0
end

function is_testcase_included(e_paths, i_paths, e_tags, i_tags, tcinfo::TestcaseInfo)
    full_tag = path_string(tcinfo)
    (
        (e_paths === nothing || !occursin(e_paths, full_tag))
        && (i_paths === nothing || occursin(i_paths, full_tag))
        && (isempty(e_tags) || isempty(intersect(e_tags, tcinfo.tags)))
        && (isempty(i_tags) || !isempty(intersect(i_tags, tcinfo.tags)))
        )
end

function filter_testcases(run_options, tcs)
    e_paths = run_options[:exclude]
    i_paths = run_options[:include_only]
    e_tags = Set(run_options[:exclude_tags])
    i_tags = Set(run_options[:include_only_tags])
    filter(p -> is_testcase_included(e_paths, i_paths, e_tags, i_tags, p[1]), tcs)
end

function runtests_internal(run_options, tcs, doctest=false)
    if run_options[:verbosity] > 0
        println("Collecting testcases...")
    end
    all_testcases = get_testcases(tcs)
    testcases = filter_testcases(run_options, all_testcases)

    if length(testcases) == 0 && run_options[:verbosity] > 0
        println("All $(length(all_testcases)) testcases were filtered out, nothing to run")
        return 0
    end

    if run_options[:verbosity] > 0
        println(
            "Using $(length(testcases)) out of $(length(all_testcases)) testcase definitions...")
        println("=" ^ 80)
    end

    all_successful = run_testcases(run_options, testcases, doctest)
    Int(!all_successful)
end

struct TestcaseAccumID
end

const TESTCASE_ACCUM_ID = TestcaseAccumID()

struct DoctestsFlagID
end

const DOCTESTS_FLAG_ID = DoctestsFlagID()

"""
Enable doctest mode for the reporting.
Replaces all variable parts of the reports (timings, module versions etc) with placeholders.
Warning: non-pure, sets a global flag.
"""
function jute_doctest()
    task_local_storage(DOCTESTS_FLAG_ID, true)
end

"""
    runtests(; options=nothing, ignore_commandline=false)

Run the test suite.

This function has several side effects:

* it parses the command-line arguments, using them to build the dictionary of run options
  (see [Run options](@ref run_options_manual) in the manual for the list);
* it picks up and includes the test files, selected according to the options.

`options` must be a dictionary with the keys corresponding to some of the options
from the above list. If `options` is given, command-line arguments are not parsed.

If `ignore_commandline` is `true`, command-line arguments passed to the program will not be used.
This option can be helpful if one wants to be sure that the run options used are exactly the ones
specified in the call.

Returns `0` if there are no failed tests, `1` otherwise.
"""
function runtests(tcs=nothing; options=nothing, ignore_commandline=false)
    args = ignore_commandline ? nothing : ARGS
    run_options = build_run_options(args=args, options=options)

    if tcs === nothing

        if haskey(task_local_storage(), TESTCASE_ACCUM_ID)
            tcs = task_local_storage(TESTCASE_ACCUM_ID)
            delete!(task_local_storage(), TESTCASE_ACCUM_ID)
        else
            runtests_dir = get_runtests_dir()
            test_files = find_test_files(runtests_dir, run_options[:test_file_postfix])

            if run_options[:verbosity] > 0
                println("Loading test files...")
            end

            tcs = task_local_storage(TESTCASE_ACCUM_ID, Any[]) do
                include_test_files!(
                    test_files, run_options[:dont_add_runtests_path] ? nothing : runtests_dir)
                task_local_storage(TESTCASE_ACCUM_ID)
            end
        end
    end

    doctest = haskey(task_local_storage(), DOCTESTS_FLAG_ID) && task_local_storage(DOCTESTS_FLAG_ID)

    res = runtests_internal(run_options, tcs, doctest)

    if doctest
        nothing
    else
        res
    end
end


text:
using Highlights.Tokens, Highlights.Themes

import Parameters: @with_kw

"""
    Theme

Stores colors for different semantically relevant items, used to
style outputs to terminal.
"""
@with_kw struct Theme
    docstring::String = "#a7de97"  # green
    string::String = "#A5D6A7"  # green
    type::String = "#d880e7"  # purple
    emphasis::String = "blue  bold"
    emphasis_light::String = "#bfe0fd"  # pale green
    code::String = "#ffd77a"  # light blue
    multiline_code::String = "#ffd77a"  # light blue
    symbol::String = "#ffc94d"  # orange light
    expression::String = "#F48FB1"  # pink light
    number::String = "#90CAF9"  # blue
    operator::String = "#fc6262 bold" # red
    func::String = "#fff126" # yellow

    # loggin levels
    info::String = "#90CAF9"  # blue
    debug::String = "#90CAF9"  # blue
    warn::String = "#ff681c"  # orange
    error::String = "bold red"  # orange
end

theme = Theme() # default theme

# ------------------------------ Highlighters.jl ----------------------------- #

"""
Custom hilighting theme for Highlighters.jl
https://juliadocs.github.io/Highlights.jl/stable/man/theme/
"""
abstract type CodeTheme <: AbstractTheme end

@theme CodeTheme Dict(
    :style => S"",
    :tokens => Dict(
        # TEXT    => S"fg: e6e8e6",
        # yellow
        NAME_FUNCTION => S"fg: FFF59D; bold",
        NAME_OTHER => S"fg: FFF59D; bold",

        # red
        KEYWORD => S"fg: fc6262; bold",
        OPERATOR => S"fg: fc6262; bold",
        PUNCTUATION => S"fg: fc7474",

        # green
        STRING => S"fg: A5D6A7",
        COMMENT => S"fg: C5E1A5; italic",
        STRING_DOC => S"fg: D4E157",

        # blue
        NUMBER => S"fg: 90CAF9",
    ),
)


text:
module TBhelpers

using JuLIP
using LinearAlgebra
using ACEtb.SlaterKoster

export set_SK_orbitals, find_row_in_matrix, get_latvec_ids, get_R, get_tbcells, wrap_shift, get_sparse_indexing, get_neighbours, get_translation_cells

function set_SK_orbitals(ntypes,nshells,orblist)
   orb_keys = Dict(0 => "s", 1 => "p", 2 => "d", 3 => "f")
   # Defining orbital keys for SK package
   k = 0
   skokeys = []
   for i = 1:ntypes
      okeys = []
      for j = 1:nshells[i]
         k += 1
         push!(okeys, SKOrbital(orb_keys[orblist[k]]))
      end
      push!(skokeys, okeys)
   end

   # generate bonds automatically
   SKH_list = []
   for ki in skokeys
       push!(SKH_list, SKH([i for i in ki]))
   end
   return SKH_list
end

function find_row_in_matrix(row, mat; round_digits=nothing)
    if round_digits != nothing
        row = round.(row, digits=round_digits)
    end
    rtn = []
    for i=1:size(mat,1)
        if round_digits != nothing
            matrow = round.(mat[i,:], digits=round_digits)
        else
            matrow = mat[i,:]
        end
        if row == matrow
            push!(rtn,i)
        end
    end
    return rtn
end

function get_latvec_ids(x,y,z)
    # Lattice vectors
    inds = [CartesianIndices((x,y,z))...]
    indx = [convert(Array{Int64,1}, [d[1],d[2],d[3]]) for d in inds]
    R_cN = hcat(indx...)
    N_c = Array{Int64,2}([x y z]')
    R_cN .+= N_c .÷ 2
    R_cN .%= N_c
    R_cN .-= N_c .÷ 2
    return R_cN
end

function get_latvec_ids(v::AbstractArray)
    # Lattice vectors
    inds = [CartesianIndices(tuple(v...))...]
    indx = [convert(Array{Int64,1}, [d[1],d[2],d[3]]) for d in inds]
    R_cN = hcat(indx...)
    N_c = Array{Int64,2}(v')
    R_cN .+= N_c .÷ 2
    R_cN .%= N_c
    R_cN .-= N_c .÷ 2
    return R_cN'
end

function get_latvec_ids(vmin::AbstractArray, vmax::AbstractArray; origin=true)
    if origin
       R_cN = [[0 0 0]]
    else
       R_cN = []
    end
    for i = vmin[1]:vmax[1]
       for j = vmin[1]:vmax[1]
          for k = vmin[1]:vmax[1]
              if ( origin && i==0 && j==0 && k==0 )
                 continue
              end
              push!(R_cN, [i j k])
          end
       end
    end
    return R_cN
end

function get_latvec_ids(v::AbstractArray)
    @assert length(v) < 3
    # Lattice vectors
    inds = [CartesianIndices(tuple(v))...]
    indx = [convert(Array{Int64,1}, [d[1],d[2],d[3]]) for d in inds]
    R_cN = hcat(indx...)
    N_c = Array{Int64,2}(v')
    R_cN .+= N_c .÷ 2
    R_cN .%= N_c
    R_cN .-= N_c .÷ 2
    return R_cN
end

function get_R(atoms, i, j; cell_vec=nothing)
    # i is always at unitcell (0,0,0)
    # @show atoms
    # @show i,j
    pos = positions(atoms)
    # @show pos
    if cell_vec != nothing
        c = cell(atoms)
        return pos[j] + (c * cell_vec) - pos[i]
    else
        return pos[j] - pos[i]
    end
end

function get_translation_cells(cell, cutoff; pos=1, neg=1)
    vmin = [-1 -1 -1]
    vmax = [1 1 1]
    invCell = inv(cell)

    for d = 1:3
       icell = norm(invCell[d,:])
       l = floor(cutoff * icell)
       vmin[d] = -(neg + l)
       vmax[d] = pos + l
    end
    return get_latvec_ids(vmin, vmax)
end

function get_tbcells(atoms, invcell, mesh, ia)
    tbcells = zeros(Int64, length(atoms), 3)
    for ja = 1:length(atoms)
        Rij = get_R(atoms, ia, ja)
        tbcells[ja,:] = convert(Array{Int64,1}, round.(invcell * Rij))
    end
    return tbcells
end

function wrap_shift(shift, mesh)
    m = mesh .÷ 2
    for a = 1:length(shift)
        if shift[a] < - m[a]
            shift[a] += mesh[a]
        elseif a > m[a]
            shift[a] -= mesh[a]
        end
    end
    return shift
end

function get_neighbours(atoms, tbcells, cutoff)

    nNeighs = [ 0 for i=1:length(atoms)]
    iNeighs = [[] for i=1:length(atoms)]
    img2CentCell = []
    iCellVec = []
    coords = []
    species = []
    cell = atoms.cell
    sqcut = cutoff^2

    nAllAtom = 0
    # Outer loop: All atoms in all translation cells
    for c = 1:length(tbcells)
        oldAtom_j = 0
        for j = 1:length(atoms)
            vec = tbcells[c] * cell
            neigh_pos = atoms.X[j][:] + vec[:]
            n = 0
            # Inner loop: All atoms in central cell
            for i = 1:length(atoms)
                dist = norm(neigh_pos .- atoms.X[i])
                if (dist > sqcut)
                    continue
                end
                if j != oldAtom_j
                   nAllAtom += 1
                   push!(coords, neigh_pos)
                   push!(img2CentCell, j)
                   push!(species, atoms.Z[j])
                   push!(iCellVec, c)
                   oldAtom_j = j
                end

                # Do not add self
                if (c == 1) && (i == j)
                   continue
                end

                nNeighs[i] += 1
                push!(iNeighs[i],nAllAtom)
            end
        end
    end
    return nNeighs, iNeighs, img2CentCell, iCellVec, coords, species
end

function get_sparse_indexing(atoms, nNeighs, iNeighs, img2CentCell, norbs)

    nAtoms = size(iNeighs,1)
    mNeighs = size(iNeighs,2)
    iPairs= [[] for i=1:nAtoms+1]
    all_types = [ atoms.Z[i].z for i=1:length(atoms)]
    types = unique!(sort!(all_types))
    a2t = [findfirst(isequal(atoms.Z[i].z), types) for i=1:length(atoms) ]

    ind = 1
    for i = 1:nAtoms
      isp = a2t[i]

      # Add onsite block starting index
      push!(iPairs[i],ind)
      ind += norbs[isp] * norbs[isp]

      # Add off-site cell blocks starting indexes for each neighbour
      for inei = 1:nNeighs[i]
        jsp = a2t[img2CentCell[iNeighs[i][inei]]]
        push!(iPairs[i],ind)
        ind += norbs[isp] * norbs[jsp]
      end

    end
    return iPairs
end

end # End of module


text:
struct Configuration{D,T,M,U}
  coord::Vector{NTuple{D,T}}
  tcode::Vector{Int}
  label::Vector{U}
end

function Configuration(lattice::Lattice{D,T,M,U}) where {D,T,M,U}
  coord, tcode = __simple_copy(lattice)
  label = lattice.types

  Configuration{D,T,M,U}(coord, tcode, label)
end

##### make sure SamplePath works for Configuration
SamplePath(xs::AbstractVector{T}, ts, dims::NTuple{N}) where {T <: Configuration, N} = SamplePath{T, N, typeof(xs), typeof(ts)}(xs, ts)

SamplePath(xs::AbstractVector{T}, ts::AbstractVector) where T <: Configuration = SamplePath(xs, ts, (1, 1))

##### plotting recipes

@recipe function f(config::Configuration{D,T,VonNeumann}) where {D,T}
  # unpack information
  coord = config.coord
  tcode = config.tcode
  label = config.label

  # set plot settings
  seriestype := :scatter
  aspect_ratio --> 1
  markerstrokewidth --> 0
  markershape --> :square
  grid --> nothing

  for t in label
    idx = findall(isequal(last(t)), tcode)

    @series begin
      label --> first(t)
      coord[idx]
    end
  end
end

@recipe function f(config::Configuration{2,T,Hexagonal}) where T
  # unpack information
  coord = config.coord
  tcode = config.tcode
  label = config.label

  # set plot settings
  seriestype := :scatter
  aspect_ratio --> 1
  markerstrokewidth --> 0
  markershape --> :hexagon
  grid --> nothing

  A = cos(pi / 3)
  B = sin(pi / 3)

  for t in label
    idx = findall(isequal(last(t)), tcode)
    new_coord = Vector{NTuple{2,Float64}}(undef, length(idx))

    @series begin
      label --> first(t)

      for i in eachindex(idx)
        c = coord[idx[i]]

        x = c[1] + A * c[2]
        y = B * c[2]

        new_coord[i] = (x, y)
      end

      new_coord
    end
  end
end

# display the initial and final configurations
@recipe function f(xw::SamplePath{T}) where T <: Configuration
  cfg_A = xw.u[1]
  cfg_B = xw.u[end]

  layout := 2

  # initial configuration
  @series begin
    subplot := 1
    cfg_A
  end

  # final configuration
  @series begin
    subplot := 2
    cfg_B
  end
end


text:
Base.@kwdef mutable struct Node
    P = nothing
    N = nothing
    E = nothing
    F = nothing
end

function Node(Point, Pos)
  Node(; N = Pos, P = Point')
end


text:
@auto_hash_equals struct DropObs <: Imputor
    vardim::Int
    context::AbstractContext
end

"""
    DropObs(; vardim=2, context=Context)

Removes missing observations from the `AbstractArray` or `Tables.table` provided.

# Keyword Arguments
* `vardim=2::Int`: Specify the dimension for variables in matrix input data
* `context::AbstractContext=Context()`: A context which keeps track of missing data
  summary information

# Example
```jldoctest
julia> using Impute: DropObs, Context, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Array{Union{Missing, Float64},2}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, DropObs(; vardim=1, context=Context(; limit=1.0)))
2×3 Array{Union{Missing, Float64},2}:
 1.0  2.0  5.0
 1.1  2.2  5.5
```
"""
DropObs(; vardim=2, context=Context()) = DropObs(vardim, context)

function impute!(data::AbstractVector, imp::DropObs)
    imp.context() do c
        filter!(x -> !ismissing(c, x), data)
    end
end

function impute!(data::AbstractMatrix, imp::DropObs)
    imp.context() do c
        return filterobs(imp, data) do obs
            !ismissing(c, obs)
        end
    end
end

# Deleting elements from subarrays doesn't work so we need to collect that data into
# a separate array.
impute!(data::SubArray, imp::DropObs) = impute!(collect(data), imp::DropObs)

function impute!(table, imp::DropObs)
    imp.context() do c
        @assert istable(table)
        rows = Tables.rows(table)

        # Unfortunately, we'll need to construct a new table
        # since Tables.rows is just an iterator
        filtered = Iterators.filter(rows) do r
            !any(x -> ismissing(c, x), propertyvalues(r))
        end

        table = materializer(table)(filtered)
        return table
    end
end

@auto_hash_equals struct DropVars <: Imputor
    vardim::Int
    context::AbstractContext
end

"""
    DropVars(; vardim=2, context=Context())

Finds variables with too many missing values in a `AbstractMatrix` or `Tables.table` and
removes them from the input data.

# Keyword Arguments
* `vardim=2::Int`: Specify the dimension for variables in matrix input data
* `context::AbstractContext`: A context which keeps track of missing data
  summary information

# Examples
```jldoctest
julia> using Impute: DropVars, Context, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Array{Union{Missing, Float64},2}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, DropVars(; vardim=1, context=Context(; limit=0.2)))
1×5 Array{Union{Missing, Float64},2}:
 1.1  2.2  3.3  missing  5.5
```
"""
DropVars(; vardim=2, context=Context()) = DropVars(vardim, context)

function impute!(data::AbstractMatrix, imp::DropVars)
    return filtervars(imp, data) do var
        try
            imp.context() do c
                for x in var
                    ismissing(c, x)
                end
            end
            return true
        catch e
            if isa(e, ImputeError)
                return false
            else
                rethrow(e)
            end
        end
    end
end

function impute!(table, imp::DropVars)
    istable(table) || throw(MethodError(impute!, (table, imp)))
    cols = Tables.columns(table)

    cnames = Iterators.filter(propertynames(cols)) do cname
        try
            imp.context() do c
                col = getproperty(cols, cname)
                for i in eachindex(col)
                    ismissing(c, col[i])
                end
            end
            return true
        catch e
            if isa(e, ImputeError)
                return false
            else
                rethrow(e)
            end
        end
    end

    selected = Tables.select(table, cnames...)
    table = materializer(table)(selected)
    return table
end


text:
"""
    UncertainEmpiricalVectorValue

An empirical value represented by a distribution estimated from actual data.

## Fields
- **`distribution`** The distribution describing the value.
- **`values`**: The values from which `distribution` is estimated.
"""
struct UncertainEmpiricalVectorValue{D <: Distribution, T} <: AbstractEmpiricalValue
    distribution::D
    values::AbstractVector{AbstractVector{T}}
end


text:
module FESOM
using UUIDs
using Random
using Statistics
using LinearAlgebra

export fesom

function quality_of_fit(isf_m::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1})
    mean(broadcast((x,y,z) -> ((x - y)/z)^2,isf,isf_m,isf_error))
end

function quality_of_fit!(isf_m::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1})
    mean(broadcast!((x,y,z) -> ((x - y)/z)^2,isf_m,isf,isf_m,isf_error))
end

function set_isf_term(imaginary_time::Array{Float64,1},frequency_bins::Array{Float64,1},beta::Float64)
    b = beta;
    isf_term = Array{Float64,2}(undef,size(imaginary_time,1),size(frequency_bins,1));
    for i in 1:size(imaginary_time,1)
        t = imaginary_time[i];
        for j in 1:size(frequency_bins,1)
            f = frequency_bins[j];
            isf_term[i,j] = (exp(-t*f) + exp(-(b - t)*f));
            #isf_term[i,j] = (exp(-t*f) + exp(-(b - t)*f))/(1 + exp(-b*f));
        end
    end
    isf_term
end

function set_isf_trapezoidal!( isf_m::Array{Float64,1}, isf_m2::Array{Float64,1},
                               isf_term::Array{Float64,2}, isf_term2::Array{Float64,2},
                               dsf::Array{Float64,1}
                               )
    mul!(isf_m,isf_term,dsf);
    mul!(isf_m2,isf_term2,dsf);
    isf_m .+= isf_m2;
    nothing
end

function fesom( dsf::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1},
                frequency::Array{Float64,1},imaginary_time::Array{Float64,1};
                temperature::Float64 = 1.2,
                number_of_iterations::Int64 = 10000,
                stop_minimum_fitness::Float64 = 1.0,
                seed::Int64 = 1,
                track_stats::Bool = false,
                number_of_blas_threads::Int64 = 0)
    rng = MersenneTwister(seed);
    fesom(rng,dsf,isf,isf_error,frequency,imaginary_time,
          temperature = temperature,
          number_of_iterations = number_of_iterations,
          stop_minimum_fitness = stop_minimum_fitness,
          number_of_blas_threads = number_of_blas_threads)
end

function fesom( rng::MersenneTwister,
                dsf::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1},
                frequency::Array{Float64,1},imaginary_time::Array{Float64,1};
                temperature::Float64 = 1.2,
                number_of_iterations::Int64=10000,
                stop_minimum_fitness::Float64 = 1.0,
                track_stats::Bool = false,
                number_of_blas_threads::Int64 = 0)
    #_bts = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
    #println("Number of BLAS threads: $(_bts)");
    if number_of_blas_threads > 0
        LinearAlgebra.BLAS.set_num_threads(number_of_blas_threads)
    end

    beta = 1/temperature;
    moment0 = isf[1];
    random_vector = Array{Float64,1}(undef,size(dsf,1));
    dsf_new = Array{Float64,1}(undef,size(dsf,1));

    isf_term = set_isf_term(imaginary_time,frequency,beta);
    isf_term2 = copy(isf_term);

    df = frequency[2:end] .- frequency[1:size(frequency,1) - 1];
    dfrequency1 = zeros(size(frequency,1));
    dfrequency2 = zeros(size(frequency,1));
    for i in 1:(size(frequency,1) - 1)
        dfrequency1[i] = df[i]/2
        dfrequency2[i+1] = df[i]/2
    end
    isf_term .*= dfrequency1';
    isf_term2 .*= dfrequency2';

    isf_m = Array{Float64,1}(undef,size(isf,1));
    isf_m2 = Array{Float64,1}(undef,size(isf,1));

    set_isf_trapezoidal!(isf_m, isf_m2, isf_term, isf_term2, dsf );

    χ2 = quality_of_fit!(isf_m,isf,isf_error);
    χ2_new = 0.0;
    normalization = 1.0;
    normalization_term1 = (1 .+ exp.(-beta .* frequency)) .* dfrequency1;
    normalization_term2 = (1 .+ exp.(-beta .* frequency)) .* dfrequency2;
    normalization_term = normalization_term1 .+ normalization_term2;

    if track_stats
        minimum_fitness = zeros(number_of_iterations);
    end
    nsteps = 0;
    for i in 1:number_of_iterations
        randn!(rng,random_vector);
        broadcast!((x,y) -> abs(x*(1 + y)),dsf_new, dsf, random_vector);

        #normalization: store intermediate results in random_vector
        normalization = dot(dsf_new,normalization_term);
        broadcast!(*,dsf_new, dsf_new, moment0/normalization);

        set_isf_trapezoidal!(isf_m, isf_m2, isf_term, isf_term2, dsf_new );

        χ2_new = quality_of_fit!(isf_m,isf,isf_error);
        if χ2_new < χ2
            dsf_new, dsf = dsf, dsf_new;
            χ2_new, χ2 = χ2, χ2_new;
        end

        nsteps +=1 ;

        if track_stats
            minimum_fitness[i] = χ2;
        end

        if χ2 < stop_minimum_fitness
            break
        end
    end
    println(χ2)

    if track_stats
        return rng,dsf,χ2,minimum_fitness[1:nsteps]
    else
        return rng,dsf,χ2,zeros(2)
    end

end
end


text:
include("A2C.jl")
include("ppo.jl")
include("A2CGAE.jl")
include("ddpg.jl")


text:
# Shaded line plots

function PlotShaded(label_id, values::AbstractArray{T}, count::Integer, y_ref::Real = 0., xscale::Real = 1., x0::Real = 0., offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotShaded(label_id, values, count, y_ref, xscale, x0, offset, stride)
end

function PlotShaded(label_id, x::AbstractArray{T}, y::AbstractArray{T}, count::Integer, y_ref::Real = 0., offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, stride)
end

function PlotShaded(label_id, x::AbstractArray{T}, y1::AbstractArray{T}, y2::AbstractArray{T}, count::Integer, offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, stride)
end

function PlotShaded(label_id, values::AbstractArray{T}, count::Integer, y_ref::Real = 0., xscale::Real = 1., x0::Real = 0., offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotShaded(label_id, Float64.(values), count, y_ref, xscale, x0, offset, stride)
end

function PlotShaded(label_id, x::AbstractArray{T}, y::AbstractArray{T}, count::Integer, y_ref::Real = 0., offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotShaded(label_id, Float64.(x), Float64.(y), count, y_ref, offset, stride)
end

function PlotShaded(label_id, x::AbstractArray{T}, y1::AbstractArray{T}, y2::AbstractArray{T}, count::Integer, offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotShaded(label_id, Float64.(x), Float64.(y1), Float64.(y2), count, offset, stride)
end

function PlotShaded(x::AbstractArray{T,1}, y1::AbstractArray{T,1}, y2::AbstractArray{T,1};
                               label_id::String = "", count::Integer = min(length(x), length(y2),
                               length(y2)), offset::Integer = 0, stride::Integer = 1) where {T <: ImPlotData}

    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, stride * sizeof(T))
end

function PlotShaded(x::AbstractArray{T1}, y1::AbstractArray{T2}, y2::AbstractArray{T3};
                    kwargs...) where {T1<:Real, T2<:Real, T3<:Real}
    PlotShaded(promote(x,y1,y2)...; kwargs...)
end

function PlotShaded(x::AbstractArray{T}, y::AbstractArray{T}, y_ref::T;
                    count::Integer = min(length(x),length(y)), offset::Integer = 0,
                    stride::Integer = 1, label_id::String = "") where {T <: ImPlotData}

    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, stride * sizeof(T))
end

function PlotShaded(x::AbstractArray{T1}, y::AbstractArray{T2}, y_ref::T3;
                    kwargs...) where {T1<:Real, T2<:Real, T3<:Real}
    (x, y) = promote(x, y)
    y_ref = eltype(x)(y_ref)
    PlotShaded(x, y, y_ref; kwargs...)
end

function PlotShaded(x::UnitRange{<:Integer}, y1::AbstractArray{T}, y2::AbstractArray{T}; label_id::String="") where {T<:ImPlotData}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")
    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, sizeof(T))
end

function PlotShaded(x::UnitRange{<:Integer}, y::AbstractArray{T}, y_ref::T; label_id::String="") where {T <: ImPlotData}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")
    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, sizeof(T))
end

function PlotShaded(x::StepRange, y1::AbstractArray{T}, y2::AbstractArray{T}; label_id::String="") where {T <: ImPlotData}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start -1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")
    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, stride * sizeof(T))
end

function PlotShaded(x::StepRange, y::AbstractArray{T}, y_ref::T; label_id::String = "") where {T <: ImPlotData}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start -1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")

    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, stride * sizeof(T))
end

function PlotShaded(x::OrdinalRange, y1::AbstractArray{T1}, y2::AbstractArray{T2};
                    kwargs...) where {T1<:Real, T2<:Real}
    PlotShaded(x, promote(y1,y2)...; kwargs...)
end

function PlotShaded(x::OrdinalRange, y::AbstractArray{T1}, y_ref::T2;
                    kwargs...) where {T1<:Real, T2<:Real}
    y_ref = eltype(y)(y_ref)
    PlotShaded(x, y, y_ref; kwargs...)
end

# xfield, yfield should be propertynames of eltype(structvec)
function PlotShaded(
    structvec::Vector{T}, xfield::Symbol, yfield::Symbol, y_ref::Real;
    count::Integer = length(structvec), offset::Integer = 0,
    stride::Integer = 1, label_id::String = ""
) where T

    Tx = fieldtype(T, xfield)
    Ty = fieldtype(T, yfield)
    x_offset = fieldoffset(T, Base.fieldindex(T, xfield))
    y_offset = fieldoffset(T, Base.fieldindex(T, yfield))
    x_ptr = (pointer(structvec, 1) + x_offset) |> Ptr{Tx}
    y_ptr = (pointer(structvec, 1) + y_offset) |> Ptr{Ty}

    if !T.mutable
        # this is somewhat illegal and is used only to pass a pointer through AbstractArray argument into ccall
        x = unsafe_wrap(Vector{Tx}, x_ptr, size(structvec); own = false)
        y = unsafe_wrap(Vector{Ty}, y_ptr, size(structvec); own = false)
        stride = stride * sizeof(T)
    else # two new vectors every 1/60 second...
        x = Vector{Tx}(undef, length(structvec))
        y = Vector{Ty}(undef, length(structvec))
        for (i, val) in enumerate(structvec)
            x[i] = getproperty(val, xfield)
            y[i] = getproperty(val, yfield)
        end
        if Tx !== Ty
            x, y = promote(x, y)
        end
        stride = stride * sizeof(eltype(x))
    end

    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, stride)
end

PlotShadedG(label_id, getter1, data1, getter2, data2, count, offset=0) =
LibCImPlot.PlotShadedG(label_id, getter1, data1, getter2, data2, count, offset)


text:
abstract type AbstractProgressiveHedging end

StochasticPrograms.num_subproblems(ph::AbstractProgressiveHedging) = StochasticPrograms.num_subproblems(ph.structure)

# Initialization #
# ======================================================================== #
function initialize!(ph::AbstractProgressiveHedging, penaltyterm::AbstractPenaltyterm)
    # Initialize progress meter
    ph.progress.thresh = sqrt(ph.parameters.ϵ₁ ^ 2 + ph.parameters.ϵ₂ ^ 2)
    # Initialize subproblems
    initialize_subproblems!(ph, scenarioproblems(ph.structure), penaltyterm)
    # Initialize penalty parameter (if applicable)
    ph.data.δ₁ = 1.0
    initialize_penalty!(ph)
    # Finish initialization
    finish_initilization!(ph, penalty(ph))
end
# ======================================================================== #

# Functions #
# ======================================================================== #
function decision(ph::AbstractProgressiveHedging)
    return ph.ξ
end

function decision(ph::AbstractProgressiveHedging, index::MOI.VariableIndex)
    i = something(findfirst(i -> i == index, ph.decisions.undecided), 0)
    if iszero(i)
        throw(MOI.InvalidIndex(index))
    end
    return decision(ph)[i]
end

function current_objective_value(ph::AbstractProgressiveHedging, Qs::AbstractVector)
    return sum(Qs)
end
current_objective_value(ph) = current_objective_value(ph, ph.subobjectives)

function get_objective_value(ph::AbstractProgressiveHedging)
    if !isempty(ph.Q_history)
        return ph.Q_history[end]
    else
        return calculate_objective_value(ph)
    end
end

function objective_value(ph::AbstractProgressiveHedging)
    return ph.data.Q
end

function log!(ph::AbstractProgressiveHedging; optimal = false, status = nothing)
    @unpack Q, δ₁, δ₂, iterations = ph.data
    @unpack ϵ₁, ϵ₂, keep, offset, indent = ph.parameters
    # Early termination log
    if status != nothing && ph.parameters.log
        ph.progress.thresh = Inf
        ph.progress.printed = true
        val = if status == MOI.INFEASIBLE
            Inf
        elseif status == MOI.DUAL_INFEASIBLE
            -Inf
        else
            0.0
        end
        ProgressMeter.update!(ph.progress, val,
                              showvalues = [
                                  ("$(indentstr(indent))Objective", val),
                                  ("$(indentstr(indent))Early termination", status),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
        return nothing
    end
    push!(ph.Q_history, Q)
    push!(ph.primal_gaps, δ₁)
    push!(ph.dual_gaps, δ₂)
    ph.data.iterations += 1
    if ph.parameters.log
        δ = sqrt(δ₁ ^ 2 + δ₂ ^ 2)
        if δ <= ph.progress.thresh && !(δ₁ <= ϵ₁ && δ₂ <= ϵ₂)
            δ = 1.01 * ph.progress.thresh
        end
        δ = optimal ? 0.0 : δ
        ProgressMeter.update!(ph.progress, δ,
                              showvalues = [
                                  ("$(indentstr(indent))Objective",Q),
                                  ("$(indentstr(indent))Primal gap", δ₁),
                                  ("$(indentstr(indent))Dual gap", δ₂),
                                  ("$(indentstr(indent))Iterations", iterations)
                              ], keep = keep, offset = offset)
    end
end

function indentstr(n::Integer)
    return repeat(" ", n)
end

function check_optimality(ph::AbstractProgressiveHedging)
    @unpack ϵ₁, ϵ₂ = ph.parameters
    @unpack δ₁, δ₂ = ph.data
    return δ₁ <= ϵ₁ && δ₂ <= ϵ₂
end
# ======================================================================== #
function show(io::IO, ph::AbstractProgressiveHedging)
    println(io, typeof(ph).name.name)
    println(io, "State:")
    show(io, ph.data)
    println(io, "Parameters:")
    show(io, ph.parameters)
end

function show(io::IO, ::MIME"text/plain", ph::AbstractProgressiveHedging)
    show(io, ph)
end


text:
# This file only forward to Pkg's command but under --project by default
module PkgCmd

using Pkg
using ..InstallCmd
using Comonicon.Tools: prompt, cmd_error
using ..Options

function withproject(command, glob, action_msg, compile_min=true, ion::Options.Ion=Options.read())
    script = "using Pkg;"
    if !glob
        msg = "cannot $action_msg in global environment, use -g, --glob to $action_msg to global environment"
        script *= "startswith(Base.active_project(), DEPOT_PATH[1]) && (print(\"$msg\"); exit(1););"
        # script *= "println(ENV);"
        # TODO: use Pkg API directly later
        script *= "if isdefined(Pkg.REPLMode, :PRINTED_REPL_WARNING); Pkg.REPLMode.PRINTED_REPL_WARNING[] = true; end;" # supress REPL warning
    end

    script *= command
    exename = Options.active_julia_bin(ion)

    options = []
    if compile_min
        push!(options, "--compile=min")
    end
    push!(options, "-g1", "--color=yes", "--startup-file=no", "-e", script)

    # PackageCompiler will set depot path to binary folder
    # we will want to use global default, unless user specified.
    if glob
        withenv("JULIA_PROJECT"=>nothing, "JULIA_LOAD_PATH"=>nothing, "JULIA_DEPOT_PATH"=>nothing) do
            run(ignorestatus(Cmd([exename, options...])))
        end
    else
        withenv("JULIA_PROJECT"=>"@.", "JULIA_LOAD_PATH"=>nothing, "JULIA_DEPOT_PATH"=>nothing) do
            run(ignorestatus(Cmd([exename, options...])))
        end
    end
    return
end

end

"""
add package/project to the closest project.

# Arguments

- `urls`: package names or urls to add.

# Flags

- `-g, --glob`: add package to global shared environment.

"""
@cast function add(urls...; glob::Bool=false)
    isempty(urls) && cmd_error("expect a package name or url")
    packages = join(urls, " ")

    PkgCmd.withproject(
        "pkg\"add $packages\"",
        glob,
        "install a package",
    )
end

"""
Make a package available for development. If pkg is an existing local path, that path will be recorded in the manifest and used.
Otherwise, a full git clone of pkg is made. Unlike the `dev/develop` command in Julia REPL pkg mode, `ion` will clone the package
to the dev folder of the current project. You can specify `--shared` flag to use shared `dev` folder under `~/.julia/dev`
(specified by `JULIA_PKG_DEVDIR`).

# Arguments

- `url`: URL or local path to the package.

# Flags

- `-s, --shared`: controls whether to use the shared develop folder.

"""
@cast function dev(url; shared::Bool=false, glob::Bool=false)
    shared_flag = shared ? "--shared" : "--local"
    cmd = Cmd(["-e", "using Pkg; pkg\"dev $shared_flag $url\" "])
    withenv("JULIA_PROJECT"=>"@.") do
        run(`$(Base.julia_cmd()) $cmd`)
    end
end

"""
Update a package. If no posistional argument is given, update all packages in current project.

# Arguments

- `pkg`: package name.

# Options

- `--level=<int>`: control how much packages are allowed to be upgraded, major: `3`, minor: `2`, patch: `1`, fixed: `0`.
"""
@cast function update(pkgs::String...; level::Int=3, glob::Bool=false)
    if isempty(pkgs)
        cmd = "Pkg.update(;level=UpgradeLevel($level))"
    else
        pkg_str = join(map(x->"\"$pkgs\""), ", ")
        cmd = "Pkg.update([$pkg_str])"
    end

    PkgCmd.withproject(cmd, glob, "update dependencies")
end

"""
build package/project/environment

# Arguments

- `pkg`: package name.

# Flags

- `-v, --verbose`: print verbose log
- `-g, --glob`: enable to build in global shared environment
"""
@cast function build(pkg=""; verbose::Bool=false, glob::Bool=false)
    if isempty(pkg)
        cmd = "pkg\"build\""
    else
        cmd = "pkg\"build $pkg\""
    end

    PkgCmd.withproject(cmd, glob, "build", false)
end

"""
test package/project

# Arguments

- `pkg`: package name.

# Flags

- `-g, --glob`: enable to test in global shared environment
"""
@cast function test(pkg=""; glob::Bool=false)
    if isempty(pkg)
        cmd = "Pkg.test()"
    else
        cmd = "Pkg.test(\"$pkg\")"
    end

    PkgCmd.withproject(cmd, glob, "test", false)
end

"""
show current environment status

# Arguments

- `pkg`: package name.

# Flags

- `-d, --diff`: show diff to last git commit
- `-g, --glob`: enable to show status in global shared environment
"""
@cast function status(pkg=""; diff::Bool=false, glob::Bool=false)
    if isempty(pkg)
        cmd = "Pkg.status(;diff=$diff)"
    else
        cmd = "Pkg.status(\"$pkg\"; diff=$diff)"
    end

    PkgCmd.withproject(cmd, glob, "show status")
end

"""
Update the current manifest. It will update manifest
with potential changes to the dependency graph from
packages that are tracking a path.

# Flags
- `-g, --glob`: enable to resolve in global shared environment
"""
@cast function resolve(; glob::Bool=false)
    PkgCmd.withproject("Pkg.resolve()", glob, "resolve dependencies")
end

"""
Remove a package from the current project. If the mode of pkg is PKGMODE_MANIFEST also
remove it from the manifest including all recursive dependencies of pkg.

# Arguments

- `pkg`: package name to remove.

# Flags

- `-g, --glob`: enable to remove package in global shared environment
"""
@cast function rm(pkgs...; glob::Bool=false)
    isempty(pkgs) && cmd_error("expect package name")
    pkgs_str = join(map(x->"\"$x\"", pkgs), ", ")
    PkgCmd.withproject("Pkg.rm([$pkgs_str])", glob, "rm package")
end

"""
Undoes the latest change to the active project. Only states in the current session are
stored, up to a maximum of 50 states.

# Flags

- `-g, --glob`: enable to execute in global shared environment
"""
@cast undo(; glob::Bool=false) = PkgCmd.withproject("Pkg.undo()", glob, "undo")

"""
Redoes the changes from the latest undo.

# Flags

- `-g, --glob`: enable to execute in global shared environment
"""
@cast redo(; glob::Bool=false) = PkgCmd.withproject("Pkg.redo()", glob, "redo")

"""
If a Manifest.toml file exists in the active project, download all the packages
declared in that manifest. Otherwise, resolve a set of feasible packages from the
Project.toml files and install them. If no Project.toml exist
in the current active project, create one with all the dependencies in the manifest
and instantiate the resulting project.

# Flags

- `-v, --verbose`: prints the build output to stdout/stderr instead of redirecting to the build.log file.
"""
@cast function instantiate(;verbose::Bool=false)
    PkgCmd.withproject("Pkg.instantiate(;verbose=$verbose)", false, "instantiate")
end

"""
If pkg is pinned, remove the pin. If pkg is tracking a path, e.g. after Pkg.develop,
go back to tracking registered versions.

# Flags

- `-g, --glob`: enable to execute in global shared environment

"""
@cast function free(pkg; glob::Bool=false)
    PkgCmd.withproject("Pkg.free(\"$pkg\")", glob, "free $pkg")
end

# """
# Precompile all the dependencies of the project by running import on all of them in a
# new process.
# """
# @cast function precompile()
#     PkgCmd.withproject("Pkg.precompile(\"$pkg\")", false, "precompile project")
# end

"""
registry tools
"""
module Registry

using Comonicon
using Pkg
using ..PkgCmd

"""
add a registry

# Arguments

- `url`: URL to the registry, or use name "General" to add the default general registry.
"""
@cast function add(url::String)
    # NOTE: registry operations are global only
    PkgCmd.withproject("Pkg.Registry.add(\"$url\")", true, "add registry $url")
end

"""
rm a registry

# Arguments

- `name`: registry name.
"""
@cast function rm(name::String)
    PkgCmd.withproject("Pkg.Registry.rm(\"$name\")", true, "rm registry $name")
end

"""
update a registry

# Arguments

- `name`: registry name.
"""
@cast function update(name::String="")
    if isempty(name)
        PkgCmd.withproject("Pkg.Registry.update()", true, "update registry $name")
    else
        PkgCmd.withproject("Pkg.Registry.update(\"$name\")", true, "update registry $name")
    end
end

end

@cast Registry

# @cast function register()
# end


text:
using FeedbackParticleFilters, Random, Statistics, LinearAlgebra, Distributions, FillArrays

println("Testing ppFPF.jl:")

@testset "ppFPF.jl" begin
    # helper function for keeping particles between 0 and 2*pi
    function mod2pi!(x::AbstractArray)
        for i in eachindex(x)
            x[i] = mod2pi(x[i])
        end
    end

    print("  inner constructor for ppFPFState")
    ens = UnweightedParticleEnsemble(2π .* randn(1,10))
    function h(x)
        20*exp.(10*cos.(x[1] .- 2*pi*collect(0:3)/4))/exp(10.)
    end
    h(x::AbstractMatrix) = mapslices(h, x, dims=1)
    eq_dt  = PoissonEquation(x -> -sum(h(x)), ens)
    eq_dN  = PoissonEquation(x -> log.(h(x)), ens)
    state = ppFPFState(ens, eq_dN, eq_dt)
    print("-")
    @test state.ensemble == ens
    print("-")
    @test state.eq_dt == eq_dt
    print("-")
    @test state.eq_dN == eq_dN
    println("DONE")

    print("  inner constructor for ppFPF")
    st_mod = BrownianMotionCircle()
    ob_mod = CountingObservationModel{Float64, Int, typeof(h)}(1, 4, h)
    filt_prob = FilteringProblem(st_mod, ob_mod)
    g_method = DifferentialRKHSMethodS1(0.1, 0.01)
    f_method = DeterministicFlow(100, g_method, mod2pi!)
    fpf = ppFPF(st_mod, ob_mod, g_method, f_method, 10)
    print("-")
    @test fpf.state_model == st_mod
    print("-")
    @test fpf.obs_model == ob_mod
    print("-")
    @test fpf.gain_method == g_method
    print("-")
    @test fpf.flow_method == f_method
    print("-")
    @test fpf.N == 10
    println("DONE")

    print("  method initial_condition")
    print("-")
    @test initial_condition(fpf) == Uniform(0,2π)
    println("DONE")

    print("  method no_of_particles")
    print("-")
    @test no_of_particles(fpf) == 10
    print("-")
    @test no_of_particles(state) == 10
    println("DONE")

    print("  method state_model")
    print("-")
    @test state_model(fpf) == st_mod
    println("DONE")

    print("  method obs_model")
    print("-")
    @test obs_model(fpf) == ob_mod
    println("DONE")

    print("  method gain_estimation_method")
    print("-")
    @test gain_estimation_method(fpf) == g_method
    println("DONE")

    print("  method gain_estimation_method")
    print("-")
    @test flow_method(fpf) == f_method
    println("DONE")

    print("  method for Statistics.mean")
    print("-")
    @test Statistics.mean(state) ≈ mean(ens)
    println("DONE")

    print("  method for Statistics.cov")
    print("-")
    @test Statistics.cov(state) ≈ cov(ens)
    println("DONE")

    print("  method for Statistics.var")
    print("-")
    @test Statistics.var(state) ≈ var(ens)
    println("DONE")

    print("  outer constructors for ppFPFState")
    state2 = ppFPFState(fpf)
    print("-")
    @test no_of_particles(state2) == 10
    print("-")
    @test all(0 .<= state2.ensemble.positions .<= 2π)
    filt_prob = FilteringProblem(st_mod, ob_mod)
    state3 = ppFPFState(filt_prob, 20)
    print("-")
    @test no_of_particles(state3) == 20
    println("DONE")

    print("  outer constructors for ppFPF")
    fpf2 = ppFPF(filt_prob, g_method, f_method, 100)
    print("-")
    @test no_of_particles(fpf2) == 100
    println("DONE")

    print("  method initialize")
    state4 = initialize(fpf)
    print("-")
    @test state4 isa ppFPFState
    print("-")
    @test all(0 .<= state4.ensemble.positions .<= 2π)
    println("DONE")

    print("  method assimilate!")
    oldstate = deepcopy(state)
    assimilate!([0, 1, 0, 1], state, fpf, 0.01)
    print("-")
    @test state.eq_dt.gain != oldstate.eq_dt.gain
    solve!(oldstate.eq_dt, g_method)
    @test state.eq_dt.gain ≈ oldstate.eq_dt.gain

    print("-")
    @test state.eq_dt.positions ≈ state.ensemble.positions
    print("-")
    @test state.ensemble != oldstate.ensemble
    println("DONE")

    print("  method propagate!")
    oldstate = deepcopy(state)
    propagate!(state, fpf, 0.01)
    print("-")
    @test all(state.ensemble.positions .!= oldstate.ensemble.positions)
    println("DONE")

    print("  method update!")
    fpf = ppFPF(st_mod, ob_mod, g_method, f_method, 100)
    state1 = initialize(fpf)
    state2 = deepcopy(state1)
    Random.seed!(0)
    update!(state1, fpf, [1, 1, 0, 0], 0.01)
    Random.seed!(0)
    propagate!(state2, fpf, 0.01)
    assimilate!([1, 1, 0, 0], state2, fpf, 0.01)
    print("-")
    @test state1.ensemble.positions ≈ state2.ensemble.positions
    print("-")
    @test state1.eq_dt.gain ≈ state2.eq_dt.gain
    println("DONE")
end; #ppFPF.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters