text:
# Floating point arithmetic functions

for op ∈ (
    :(+), #:(-),
    :abs, :floor, :ceil, :round,
    :inv,# :log10, :log2, :exp2, :sin, :cos,
    :sqrt, :trunc
    # :exp, :log
    )
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @inline $rename(s1::FloatingTypes) = @fastmath $op(s1)

        @vectordef $rename function Base.$op(v1) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1))
        end
    end
end
# @inline vexp10(s1::FloatingTypes) = exp10(s1)
# @vectordef vexp10 function Base.exp10(v1) where {W,T <: FloatingTypes}
    # vpow(vbroadcast(Vec{W,T}, 10), extract_data(v1))
# end
# @inline vexp10(v1::AbstractSIMDVector{W,T}) where {W,T<:FloatingTypes} = vpow(vbroadcast(Vec{W,T}, 10), v1)
@inline vsign(s1::FloatingTypes) = sign(s1)
@vectordef vsign function Base.sign(v1) where {W,T<:FloatingTypes}
    vifelse(
        extract_data(v1) == vbroadcast(Vec{W,T}, zero(T)),
        vbroadcast(Vec{W,T}, zero(T)),
        copysign(vbroadcast(Vec{W,T}, one(T)), extract_data(v1))
    )
end
# @inline function vsign(v1::SVec{W,T}) where {W,T<:FloatingTypes}
#     SVec(vsign(extract_data(v1)))
# end
# @inline function Base.sign(v1::SVec{W,T}) where {W,T<:FloatingTypes}
#     SVec(vsign(extract_data(v1)))
# end

for op ∈ (
        :(+), :(-), :(*), :(/), :(%),# :(^),
        :copysign#, :max, :min
    )
    rename = VECTOR_SYMBOLS[op]
    # exact / explicit version
    erename = Symbol(:e, rename)
    @eval begin
        @inline $rename(s1::FloatingTypes, s2::FloatingTypes) = @fastmath $op(s1, s2)

        @vectordef $rename function Base.$op(v1, v2) where {W,T <: FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
        end

        @inline $erename(s1::FloatingTypes, s2::FloatingTypes) = $op(s1, s2)

        @evectordef $erename function Base.$op(v1, v2) where {W,T <: FloatingTypes}
            llvmwrap_notfast(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
        end

        @inline function Base.$op(v1::SVec{W,T1}, v2::SVec{W,T2}) where {W,T1,T2}
            T = promote_type(T1, T2)
            $op(vconvert(SVec{W,T}, v1), vconvert(SVec{W,T}, v2))
        end
        @inline function $rename(v1::V1, v2::V2) where {W,T1,T2,V1<:AbstractSIMDVector{W,T1},V2<:AbstractSIMDVector{W,T2}}
            V = promote_vtype(V1, V2)
            $rename(vconvert(V, v1), vconvert(V, v2))
        end
        # @inline function $rename(v1::V1, v2::V2) where {V1,V2}
            # V = vpromote(V1, V2)
            # $rename(vconvert(V, v1), vconvert(V, v2))
        # end
        # @inline $rename(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:FloatingTypes} =
        #     llvmwrap(Val{$(QuoteNode(op))}(), v1, v2)
        # @inline $rename(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T<:FloatingTypes} =
        #     SVec(llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2)))
        # @inline Base.$op(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T<:FloatingTypes} =
        #     SVec(llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2)))
    end
end
@inline vfdiv(v1::Vec{W,Int32}, v2::Vec{W,Int32}) where {W} = vfdiv(vconvert(Vec{W,Float32}, v1), vconvert(Vec{W,Float32}, v2))
@inline vfdiv(v1::Vec{W,Int64}, v2::Vec{W,Int64}) where {W} = vfdiv(vconvert(Vec{W,Float64}, v1), vconvert(Vec{W,Float64}, v2))
@inline function Base.:(/)(v1::SVec{W,I},v2::SVec{W,I}) where {W,I<:Integer}
    T = sizeequivalentfloat(I)
    SVec(vfdiv(vconvert(Vec{W,T}, v1), vconvert(Vec{W,T}, v2)))
end
# @vpromote vfdiv 2
@inline vmax(x::Number, y::Number) = Base.FastMath.max_fast(x,y)
@vectordef vmax function Base.max(v1, v2) where {W,T<:FloatingTypes}
    vifelse(vgreater(extract_data(v1), extract_data(v2)), extract_data(v1), extract_data(v2))
end
@inline vmin(x::Number, y::Number) = Base.FastMath.min_fast(x,y)
@vectordef vmin function Base.min(v1, v2) where {W,T<:FloatingTypes}
    vifelse(vless(extract_data(v1), extract_data(v2)), extract_data(v1), extract_data(v2))
end

# let op = :(*)
#     rename = VECTOR_SYMBOLS[op]
#     @eval begin
#         @inline $rename(s1::FloatingTypes, s2::FloatingTypes) = $op(s1, s2)
#         # @inline $rename(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:FloatingTypes} =
#         #     VecProduct(v1, v2)
#         @inline $rename(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:FloatingTypes} =
#             VecProduct(extract_data(v1), extract_data(v2))
#         @inline $rename(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T<:FloatingTypes} =
#             SVecProduct(extract_data(v1), extract_data(v2))
#         @inline Base.$op(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T<:FloatingTypes} =
#             SVecProduct(extract_data(v1), extract_data(v2))
#     end
# end

# @inline vpow(s1::FloatingTypes, x2::Integer) = s1^x2
# @vectordef vpow function Base.:^(v1, x2::Integer) where {W,T<:FloatingTypes}
    # llvmwrap(Val{:powi}, extract_data(v1), Int(x2))
# end
# @inline function vpow(v1::Vec{W,T}, x2::Integer) where {W,T<:FloatingTypes}
#     llvmwrap(Val{:powi}, v1, Int(x2))
# end
# @inline function vpow(v1::SVec{W,T}, x2::Integer) where {W,T<:FloatingTypes}
#     SVec(llvmwrap(Val{:powi}, extract_data(v1), Int(x2)))
# end
# @inline function Base.:^(v1::SVec{W,T}, x2::Integer) where {W,T<:FloatingTypes}
#     SVec(llvmwrap(Val{:powi}, extract_data(v1), Int(x2)))
# end

@inline vflipsign(s1::FloatingTypes, s2::FloatingTypes) = flipsign(s1, s2)

@vectordef vflipsign function Base.flipsign(v1, v2) where {W,T<:FloatingTypes}
    vifelse(vsignbit(extract_data(v2)), vsub(extract_data(v1)), extract_data(v1))
end

@inline vcopysign(v1::Vec{W,T}, v2::Vec{W,U}) where {W,T,U<:Unsigned} = vcopysign(v1, vreinterpret(Vec{W,T}, v2))

# @inline vflipsign(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:FloatingTypes} =
#     vifelse(vsignbit(v2), -v1, v1)
# @inline vflipsign(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T<:FloatingTypes} =
#     SVec(vifelse(vsignbit(v2), -v1, v1))
# @inline Base.flipsign(v1::SVec{W,T}, v2::SVec{W,T}) where {W,T<:FloatingTypes} =
#     SVec(vifelse(vsignbit(v2), -v1, v1))

for op ∈ (:fma, :muladd)
# let op = :fma
    rename = VECTOR_SYMBOLS[op]
    # if op == :muladd
        # rename = Symbol(:e,rename)
    # end
    @eval begin

        @vectordef $rename function Base.$op(v1, v2, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(s1::T, v2, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), vbroadcast(Vec{W,T}, s1), extract_data(v2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(v1, s2::T, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), vbroadcast(Vec{W,T}, s2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(v1, v2, s3::T) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2), vbroadcast(Vec{W,T}, s3))
        end
        @vectordef $rename function Base.$op(s1::T, s2::T, v3) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), vbroadcast(Vec{W,T}, s1), vbroadcast(Vec{W,T}, s2), extract_data(v3))
        end
        @vectordef $rename function Base.$op(s1::T, v2, s3::T) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), vbroadcast(Vec{W,T}, s1), extract_data(v2), vbroadcast(Vec{W,T}, s3))
        end
        @vectordef $rename function Base.$op(v1, s2::T, s3::T) where {W,T<:FloatingTypes}
            llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), vbroadcast(Vec{W,T}, s2), vbroadcast(Vec{W,T}, s3))
        end
        @vpromote $rename 3
        # scalar default already set in integer_arithmetic.jl
        # @inline function $rename(v1::Vec{W,T},
        #         v2::Vec{W,T}, v3::Vec{W,T}) where {W,T<:FloatingTypes}
        #     llvmwrap(Val{$(QuoteNode(op))}(), v1, v2, v3)
        # end
        # @inline function $rename(v1::AbstractSIMDVector{W,T},
        #         v2::AbstractSIMDVector{W,T}, v3::AbstractSIMDVector{W,T}) where {W,T<:FloatingTypes}
        #     SVec(llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2), extract_data(v3)))
        # end
        # @inline function Base.$op(v1::SVec{W,T},
        #         v2::SVec{W,T}, v3::SVec{W,T}) where {W,T<:FloatingTypes}
        #     SVec(llvmwrap(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2), extract_data(v3)))
        # end
    end
end

# Type promotion

# Promote scalars of all ScalarTypes to vectors of FloatingTypes, leaving the
# vector type unchanged

for op ∈ (
        :(==), :(!=), :(<), :(<=), :(>), :(>=),
        :+, :-, :*, :/,# :^,
        :copysign, :flipsign, :max, :min, :%
    )
    rename = VECTOR_SYMBOLS[op]
    @eval begin
        @vectordef $rename function Base.$op(s1::ScalarTypes, v2) where {W,T<:FloatingTypes}
            $rename(vbroadcast(Vec{W,T}, s1), extract_data(v2))
        end
        @vectordef $rename function Base.$op(v1, s2::ScalarTypes) where {W,T<:FloatingTypes}
            $rename(extract_data(v1), vbroadcast(Vec{W,T}, s2))
        end

        @inline function Base.$op(s1::T2, v2::SVec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(vbroadcast(Vec{W,T2}, s1), vconvert(Vec{W,T2}, extract_data(v2))))
        end
        @inline function Base.$op(v1::SVec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(vconvert(Vec{W,T2}, extract_data(v1)), vbroadcast(Vec{W,T2}, s2)))
        end
        @inline function $rename(s1::T2, v2::SVec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(vbroadcast(Vec{W,T2}, s1), extract_data(v2)))
        end
        @inline function $rename(v1::SVec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            SVec($rename(extract_data(v1), vbroadcast(Vec{W,T2}, s2)))
        end
        @inline function $rename(s1::T2, v2::Vec{W,T}) where {W,T<:Integer,T2<:FloatingTypes}
            $rename(vbroadcast(Vec{W,T2}, s1), v2)
        end
        @inline function $rename(v1::Vec{W,T}, s2::T2) where {W,T<:Integer,T2<:FloatingTypes}
            $rename(v1, vbroadcast(Vec{W,T2}, s2))
        end
        @vpromote $rename 2
    end
end
for op ∈ (
        :(+), :(-), :(*), :(/), :(%),# :(^),
        :copysign#, :max, :min
    )
    # exact / explicit version
    erename = Symbol(:e, VECTOR_SYMBOLS[op])
    @eval begin

        @evectordef $erename function Base.$op(v1, s2::ScalarTypes) where {W,T <: FloatingTypes}
            $erename(extract_data(v1), vbroadcast(Vec{W,T}, s2))
        end
        @evectordef $erename function Base.$op(s1::ScalarTypes, v2) where {W,T <: FloatingTypes}
            $erename(vbroadcast(Vec{W,T}, s1), extract_data(v2))
        end
    end
end

# Reduction operations

# TODO: map, mapreduce

function getneutral(op::Symbol, ::Type{T}) where {T}
    if op == :&
        return ~zero(T)
    elseif op == :|
        return zero(T)
    elseif op == :max
        return typemin(T)
    elseif op == :min
        return typemax(T)
    elseif op == :+
        return zero(T)
    elseif op == :*
        return one(T)
    end
    throw("Op $op not recognized.")
end

# We cannot pass in the neutral element via Val{}; if we try, Julia refuses to
# inline this function, which is then disastrous for performance
@generated function llvmwrapreduce(::Val{Op}, v::Vec{W,T}) where {Op,W,T}
    @assert isa(Op, Symbol)
    z = getneutral(Op, T)
    typ = llvmtype(T)
    decls = String[]
    instrs = String[]
    n = W
    nam = "%0"
    nold,n = n, VectorizationBase.nextpow2(n)
    if n > nold
        namold,nam = nam,"%vec_$n"
        append!(instrs,
            extendvector(namold, nold, typ, n, n-nold,
                llvmtypedconst(T,z), nam))
    end
    while n > 1
        nold,n = n, div(n, 2)
        namold,nam = nam,"%vec_$n"
        vtyp = "<$n x $typ>"
        ins = llvmins(Op, n, T)
        append!(instrs, subvector(namold, nold, typ, "$(nam)_1", n, 0))
        append!(instrs, subvector(namold, nold, typ, "$(nam)_2", n, n))
        if ins[1] == '@'
            push!(decls, "declare $vtyp $ins($vtyp, $vtyp)")
            push!(instrs,
                "$nam = call $vtyp $ins($vtyp $(nam)_1, $vtyp $(nam)_2)")
        else
            push!(instrs, "$nam = $ins $vtyp $(nam)_1, $(nam)_2")
        end
    end
    push!(instrs, "%res = extractelement <$n x $typ> $nam, i32 0")
    push!(instrs, "ret $typ %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            $T, Tuple{Vec{$W,$T}}, v
        )
    end
end

@static if Base.libllvm_version >= v"8"
    @generated function vsum(v::Vec{W,T}) where {W,T<:FloatingTypes}
        decls = String[]
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fadd.f$(bits).v$(W)f$(bits)"
        push!(decls, "declare $(typ) $(ins)($(typ), $(vtyp))")
        push!(instrs, "%res = call fast $typ $ins($typ 0.0, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function vprod(v::Vec{W,T}) where {W,T<:FloatingTypes}
        decls = String[]
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fmul.f$(bits).v$(W)f$(bits)"
        push!(decls, "declare $(typ) $(ins)($(typ), $(vtyp))")
        push!(instrs, "%res = call fast $typ $ins($typ 1.0, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function reduced_add(v::Vec{W,T}, s::T) where {W,T<:FloatingTypes}
        decls = String[]
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fadd.f$(bits).v$(W)f$(bits)"
        push!(decls, "declare $(typ) $(ins)($(typ), $(vtyp))")
        push!(instrs, "%res = call fast $typ $ins($typ %1, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T},$T}, v, s
            )
        end
    end
    @generated function reduced_prod(v::Vec{W,T}, s::T) where {W,T<:FloatingTypes}
        decls = String[]
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.v2.fmul.f$(bits).v$(W)f$(bits)"
        push!(decls, "declare $(typ) $(ins)($(typ), $(vtyp))")
        push!(instrs, "%res = call fast $typ $ins($typ %1, $vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T},$T}, v, s
            )
        end
    end
    @generated function vsum(v::Vec{W,T}) where {W,T<:IntegerTypes}
        decls = String[]
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.add.v$(W)i$(bits)"
        push!(decls, "declare $(typ) $(ins)($(vtyp))")
        push!(instrs, "%res = call $typ $ins($vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function vprod(v::Vec{W,T}) where {W,T<:IntegerTypes}
        decls = String[]
        instrs = String[]
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        bits = 8sizeof(T)
        ins = "@llvm.experimental.vector.reduce.mul.v$(W)i$(bits)"
        push!(decls, "declare $(typ) $(ins)($(vtyp))")
        push!(instrs, "%res = call $typ $ins($vtyp %0)")
        push!(instrs, "ret $typ %res")
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall(
                $((join(decls, "\n"), join(instrs, "\n"))),
                $T, Tuple{Vec{$W,$T}}, v
            )
        end
    end
    @generated function vsub(v::Vec{W,T}) where {W,T<:FloatingTypes}
        typ = llvmtype(T)
        vtyp = "<$W x $typ>"
        instrs = "%res = fneg fast $vtyp %0\nret $vtyp %res"
        quote
            $(Expr(:meta, :inline))
            Base.llvmcall( $instrs, Vec{$W,$T}, Tuple{Vec{$W,$T}}, v )
        end
    end
    Base.:(-)(v::SVec{W,T}) where {W,T} = SVec(vsub(extract_data(v)))
    vsub(v::SVec{W,T}) where {W,T} = SVec(vsub(extract_data(v)))
else
    # @generated function vsub(v::Vec{W,T}) where {W,T<:FloatingTypes}
    #     typ = llvmtype(T)
    #     vtyp = "<$W x $typ>"
    #     instrs = "%res = fsub fast $vtyp zeroinitializer, %0\nret $vtyp %res"
    #     quote
    #         $(Expr(:meta, :inline))
    #         Base.llvmcall( $instrs, Vec{$W,$T}, Tuple{Vec{$W,$T}}, v )
    #     end
    # end
    # Base.:(-)(v::SVec{W,T}) where {W,T} = SVec(vsub(extract_data(v)))
    # vsub(v::SVec{W,T}) where {W,T} = SVec(vsub(extract_data(v)))
    rename = VECTOR_SYMBOLS[:-]
    @eval begin
        @vectordef $rename function Base.:(-)(v1) where {W,T<:FloatingTypes}
            llvmwrap(Val{:(-)}(), extract_data(v1))
        end
    end
end
vsub(x::FloatingTypes) = Base.FastMath.sub_fast(x)

for (name, rename, op) ∈ ((:(Base.all),:vall,:&), (:(Base.any),:vany,:|),
                                    (:(Base.maximum), :vmaximum, :max), (:(Base.minimum), :vminimum, :min),
                                    (:(Base.sum),:vsum,:+), (:(Base.prod),:vprod,:*))
    @eval begin
        @inline $rename(v::AbstractSIMDVector{W,T}) where {W,T} = llvmwrapreduce(Val{$(QuoteNode(op))}(), extract_data(v))
        @inline $name(v::SVec{W,T}) where {W,T} = llvmwrapreduce(Val{$(QuoteNode(op))}(), extract_data(v))
    end
end

# @inline vall(v::Vec{W,T}) where {W,T<:IntegerTypes} = llvmwrapreduce(Val{:&}, v)
# @inline vany(v::Vec{W,T}) where {W,T<:IntegerTypes} = llvmwrapreduce(Val{:|}, v)
# @inline vmaximum(v::Vec{W,T}) where {W,T<:FloatingTypes} =
#     llvmwrapreduce(Val{:max}, v)
# @inline vminimum(v::Vec{W,T}) where {W,T<:FloatingTypes} =
#     llvmwrapreduce(Val{:min}, v)
# @inline vprod(v::Vec{W,T}) where {W,T} = llvmwrapreduce(Val{:*}, v)
# @inline vsum(v::Vec{W,T}) where {W,T} = llvmwrapreduce(Val{:+}, v)

@generated function vreduce(::Val{Op}, v::Vec{W,T}) where {Op,W,T}
    @assert isa(Op, Symbol)
    z = getneutral(Op, T)
    stmts = String[]
    n = W
    push!(stmts, :($(Symbol(:v,n)) = v))
    nold,n = n,nextpow2(n)
    if n > nold
        push!(stmts,
            :($(Symbol(:v,n)) = $(Expr(:tuple,
                [:($(Symbol(:v,nold))[$i]) for i in 1:nold]...,
                [:(VE($z)) for i in nold+1:n]...))))
    end
    while n > 1
        nold,n = n, div(n, 2)
        push!(stmts,
            :($(Symbol(:v,n,"lo")) = $(Expr(:tuple,
                [:($(Symbol(:v,nold))[$i]) for i in 1:n]...,))))
        push!(stmts,
            :($(Symbol(:v,n,"hi")) = $(Expr(:tuple,
                [:($(Symbol(:v,nold))[$i]) for i in n+1:nold]...))))
        push!(stmts,
            :($(Symbol(:v,n)) =
                $Op($(Symbol(:v,n,"lo")), $(Symbol(:v,n,"hi")))))
    end
    push!(stmts, :(v1[1].value))
    Expr(:block, Expr(:meta, :inline), stmts...)
end

@vectordef vmaximum function Base.maximum(v) where {W,T<:IntegerTypes}
    vreduce(Val{:max}, extract_data(v))
end
@vectordef vminimum function Base.minimum(v) where {W,T<:IntegerTypes}
    vreduce(Val{:min}, extract_data(v))
end

# @inline vmaximum(v::Vec{W,T}) where {W,T<:IntegerTypes} = vreduce(Val{:max}, v)
# @inline vminimum(v::Vec{W,T}) where {W,T<:IntegerTypes} = vreduce(Val{:min}, v)

# TODO: Handle cases with vectors of different lengths correctly!
@inline vmul(x::Number, y::Number) = Base.FastMath.mul_fast(x, y)
@inline vadd(x::Number, y::Number) = Base.FastMath.add_fast(x, y)
@inline vmul(x,y,z...) = vmul(x,vmul(y,z...))
@inline vadd(x,y,z...) = vadd(x,vadd(y,z...))
@inline vadd(a,b,c,d) = evadd(vadd(a,b), vadd(c,d))
@inline vadd(a,b,c,d,e,f) = evadd(evadd(vadd(a,b), vadd(c,d)),vad(e,f))
@inline vadd(a,b,c,d,e,f,g,h) = vadd(evadd(vadd(a,b), vadd(c,d)),evadd(vadd(e,f),vadd(g,h)))
@inline Base.:(+)(a::SVec, b, c, d) = evadd(vadd(a,b), vadd(c,d))
@inline vmuladd(a::Number, b::Number, c::Number) = muladd(a, b, c)
# These intrinsics are not FastMath.

@generated function vfmadd(v1::Vec{W,T}, v2::Vec{W,T}, v3::Vec{W,T}) where {W,T<:FloatingTypes}
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    # ins = "@llvm.fma.v$(W)f$(8*sizeof(T))"
    ins = "@llvm.fmuladd.v$(W)f$(8*sizeof(T))"
    decls = "declare $vtyp $ins($vtyp, $vtyp, $vtyp)"
    instrs = "%res = call $vtyp $ins($vtyp %0, $vtyp %1, $vtyp %2)\nret $vtyp %res"
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $((decls,instrs)),
            Vec{$W,$T}, Tuple{Vec{$W,$T}, Vec{$W,$T}, Vec{$W,$T}},
            v1, v2, v3
        )
    end
end
vfmadd(a::Number, b::Number, c::Number) = muladd(a, b, c)
vfnmadd(a::Number, b::Number, c::Number) = muladd(-a, b, c)
vfmsub(a::Number, b::Number, c::Number) = muladd(a, b, -c)
vfnmsub(a::Number, b::Number, c::Number) = -muladd(a, b, c)
@inline vfmadd(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T <: Integer} = vmuladd(a, b, c)
@inline vfnmadd(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T} = vfmadd(vsub(a), b, c)
@inline vfmsub(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T} = vfmadd(a, b, vsub(c))
@inline vfnmsub(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T} = vsub(vfmadd(a, b, c))
@vpromote vfmadd 3
@vpromote vfnmadd 3
@vpromote vfmsub 3
@vpromote vfnmsub 3

@generated function vfmadd_fast(v1::Vec{W,T}, v2::Vec{W,T}, v3::Vec{W,T}) where {W,T<:FloatingTypes}
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    # ins = "@llvm.fmuladd.v$(W)f$(8*sizeof(T))"
    # decls = "declare $vtyp $ins($vtyp, $vtyp, $vtyp)"
    # instrs = "%res = call fast $vtyp $ins($vtyp %0, $vtyp %1, $vtyp %2)\nret $vtyp %res"
    # I don't really want reassoc on the fmul part, so this is what I'm going with.
    instrs = """
    %prod = fmul nnan ninf nsz arcp contract $vtyp %0, %1
    %res = fadd fast $vtyp %prod, %2
    ret $vtyp %res
    """
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall(
            $instrs, Vec{$W,$T}, Tuple{Vec{$W,$T}, Vec{$W,$T}, Vec{$W,$T}}, v1, v2, v3
        )
    end
end
vfmadd_fast(a::Number, b::Number, c::Number) = Base.FastMath.add_fast(Base.FastMath.mul_fast(a, b), c)
vfnmadd_fast(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(c, Base.FastMath.mul_fast(a, b))
vfmsub_fast(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.mul_fast(a, b), c)
vfnmsub_fast(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.add_fast(Base.FastMath.mul_fast(a, b), c))
@inline vfmadd_fast(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T <: Integer} = vmuladd(a, b, c)
@inline vfnmadd_fast(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T} = vfmadd_fast(vsub(a), b, c)
@inline vfmsub_fast(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T} = vfmadd_fast(a, b, vsub(c))
@inline vfnmsub_fast(a::Vec{W,T}, b::Vec{W,T}, c::Vec{W,T}) where {W,T} = vsub(vfmadd_fast(a, b, c))
@vpromote vfmadd_fast 3
@vpromote vfnmadd_fast 3
@vpromote vfmsub_fast 3
@vpromote vfnmsub_fast 3

# Lowers to same split mul-add llvm as the
# @inline vfmadd(a, b, c) = vadd(vmul( a, b), c)
# definition, so I wont bother implementing
# vfnmadd, vfmsub, and vfnmsub
# in this manner.

# @inline Base.:*(a::IntegerTypes, b::SVec{W,T}) where {W,T} = SVec{W,T}(a) * b
# @inline Base.:*(a::T, b::SVec{W,<:IntegerTypes}) where {W,T<:FloatingTypes} = SVec(vmul(vbroadcast(Vec{W,T}, a), vconvert(Vec{W,T}, extract_data(b))))

# const Vec{W,T} = NTuple{W,Core.VecElement{T}}
@inline vfmadd231(a::Vec, b, c) = vfmadd(a, b, c)
@inline vfnmadd231(a::Vec, b, c) = vfnmadd(a, b, c)
@inline vfmsub231(a::Vec, b, c) = vfmsub(a, b, c)
@inline vfnmsub231(a::Vec, b, c) = vfnmsub(a, b, c)
@inline vfmadd231(a::Number, b::Number, c::Number) = Base.FastMath.add_fast(Base.FastMath.mul_fast(a, b), c)
@inline vfnmadd231(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(c, Base.FastMath.mul_fast(a, b))
@inline vfmsub231(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.mul_fast(a, b), c)
@inline vfnmsub231(a::Number, b::Number, c::Number) = Base.FastMath.sub_fast(Base.FastMath.sub_fast(c), Base.FastMath.mul_fast(a, b))
@inline function vifelse(f::typeof(vfmadd231), m::Mask{W}, a, b, c) where {W}
    vifelse(m, vfmadd_fast(a, b, c), c)
end
@inline function vifelse(f::typeof(vfnmadd231), m::Mask{W}, a, b, c) where {W}
    vifelse(m, vfnmadd_fast(a, b, c), c)
end
@inline function vifelse(f::typeof(vfmsub231), m::Mask{W}, a, b, c) where {W}
    vifelse(m, vfmsub_fast(a, b, c), c)
end
@inline function vifelse(f::typeof(vfnmsub231), m::Mask{W}, a, b, c) where {W}
    vifelse(m, vfnmsub_fast(a, b, c), c)
end

@vpromote vfmadd231 3
@vpromote vfnmadd231 3
@vpromote vfmsub231 3
@vpromote vfnmsub231 3

if VectorizationBase.FMA3
    for T ∈ [Float32,Float64]
        W = 16 ÷ sizeof(T)
        local suffix = T == Float32 ? "ps" : "pd"
        typ = llvmtype(T)
        while W <= VectorizationBase.REGISTER_SIZE ÷ sizeof(T)
            vfmadd_str = """%res = call <$W x $(typ)> asm "vfmadd231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            vfnmadd_str = """%res = call <$W x $(typ)> asm "vfnmadd231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            vfmsub_str = """%res = call <$W x $(typ)> asm "vfmsub231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            vfnmsub_str = """%res = call <$W x $(typ)> asm "vfnmsub231$(suffix) \$3, \$2, \$1", "=v,0,v,v"(<$W x $(typ)> %2, <$W x $(typ)> %1, <$W x $(typ)> %0)
                ret <$W x $(typ)> %res"""
            @eval begin
                @inline function vfmadd231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    Base.llvmcall($vfmadd_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
                @inline function vfnmadd231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    Base.llvmcall($vfnmadd_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
                @inline function vfmsub231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    Base.llvmcall($vfmsub_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
                @inline function vfnmsub231(a::Vec{$W,$T}, b::Vec{$W,$T}, c::Vec{$W,$T})
                    Base.llvmcall($vfnmsub_str, Vec{$W,$T}, Tuple{Vec{$W,$T},Vec{$W,$T},Vec{$W,$T}}, a, b, c)
                end
            end
            W += W
        end
    end
end

# @generated function rsqrt_fast(x::NTuple{16,Core.VecElement{Float32}})
#     if VectorizationBase.REGISTER_SIZE == 64
#         return quote
#             $(Expr(:meta,:inline))
#             Base.llvmcall(
#             """ %rs = call <16 x float> asm "vrsqrt14ps \$1, \$0", "=x,x"(<16 x float> %0)
#                 ret <16 x float> %rs""",
#             NTuple{16,Core.VecElement{Float32}}, Tuple{NTuple{16,Core.VecElement{Float32}}}, x)
#         end
#     else
#         return quote
#             vinv(vsqrt(x))
#         end
#     end
# end
# @inline function rsqrt(x::NTuple{16,Core.VecElement{Float32}})
#     r = rsqrt_fast(x)
#     # Performs a Newton step to increase accuracy.
#     # ns = vmuladd(vmul(-0.5f0, x), vmul(r, r), 1.5f0)
#     # vmul(r, ns)
#     ns = vfma(vmul(r,r), x, -3.0f0)
#     vmul(vmul(-0.5f0, r), ns)
# end
@inline rsqrt_fast(v) = vinv(vsqrt(v))
@inline rsqrt_fast(x::SVec) = SVec(rsqrt_fast(extract_data(x)))
@inline rsqrt(x::SVec) = SVec(rsqrt(extract_data(x)))
@inline rsqrt(x) = vinv(vsqrt(x))
@inline vinv(x::IntegerTypes) = vinv(float(x))
@inline vinv(v::Vec{W,I}) where {W,I<:Union{UInt64,Int64}} = vinv(pirate_convert(Vec{W,Float64}, v))
@inline vinv(v::Vec{W,I}) where {W,I<:Union{UInt32,Int32}} = vinv(pirate_convert(Vec{W,Float32}, v))

# for accumulating vector results of different sizes.
@generated function vadd(v1::Vec{W1,T}, v2::Vec{W2,T}) where {W1,W2,T}
    @assert ispow2(W1)
    @assert ispow2(W2)
    W3 = W1
    W4 = W2
    v1e = :v1
    while W3 < W2
        W3 <<= 1
        v1e = Expr(:call, :zeropad, v1e)
    end
    v2e = :v2
    while W4 < W1
        W4 <<= 1
        v2e = Expr(:call, :zeropad, v2e)
    end
    Expr(:block, Expr(:meta, :inline), Expr(:call, :vadd, v1e, v2e))
end
@inline Base.:(+)(x::AbstractStructVec{W1,T}, y::AbstractStructVec{W2,T}) where {W1,W2,T} = SVec(vadd(extract_data(x), extract_data(y)))

@inline Base.abs2(v::SVec) = vmul(v,v)
@inline vabs2(v) = vmul(v,v)
@inline vsum(s::FloatingTypes) = s
@inline vprod(s::FloatingTypes) = s

@inline reduced_add(v::AbstractSIMDVector{W,T}, s::T) where {W,T} = Base.FastMath.add_fast(s, vsum(v))
@inline reduced_add(s::T, v::AbstractSIMDVector{W,T}) where {W,T} = Base.FastMath.add_fast(s, vsum(v))
# @inline reduced_add(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T} = vadd(v1, v2)
@inline reduced_prod(v::AbstractSIMDVector{W,T}, s::T) where {W,T} = Base.FastMath.mul_fast(s, vprod(v))
@inline reduced_prod(s::T, v::AbstractSIMDVector{W,T}) where {W,T} = Base.FastMath.mul_fast(s, vprod(v))
# @inline reduced_prod(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T} = vmul(v1, v2)

@inline reduced_add(v::SVec{W,T}, s::T) where {W,T} = reduced_add(extract_data(v), s)
@inline reduced_prod(v::SVec{W,T}, s::T) where {W,T} = reduced_prod(extract_data(v), s)
@inline reduced_add(v::AbstractSIMDVector{W,T1}, s::T2) where {W,T1,T2} = Base.FastMath.add_fast(s, convert(T2,vsum(v)))
@inline reduced_add(s::T2, v::AbstractSIMDVector{W,T1}) where {W,T1,T2} = Base.FastMath.add_fast(s, convert(T2,vsum(v)))
@inline reduced_add(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1,T2} = vadd(v1, v2)
@inline reduced_add(v1::T, v2::T) where {T<:Number} = Base.FastMath.add_fast(v1,v2)
@inline reduced_prod(v::AbstractSIMDVector{W,T1}, s::T2) where {W,T1,T2} = Base.FastMath.mul_fast(s, convert(T2,vprod(v)))
@inline reduced_prod(s::T2, v::AbstractSIMDVector{W,T1}) where {W,T1,T2} = Base.FastMath.mul_fast(s, convert(T2,vprod(v)))
@inline reduced_prod(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1,T2} = vmul(v1, v2)
@inline reduced_prod(v1::T, v2::T) where {T<:Number} = Base.FastMath.mul_fast(v1,v2)

@inline reduce_to_add(v::AbstractSIMDVector{W,T}, ::T) where {W,T} = vsum(v)
@inline reduce_to_prod(v::AbstractSIMDVector{W,T}, ::T) where {W,T} = vprod(v)
@inline reduce_to_add(v::AbstractSIMDVector{W,T1}, ::T2) where {W,T1,T2} = convert(T2,vsum(v))
@inline reduce_to_add(v::T, ::T) where {T<:Number} = v
@inline reduce_to_prod(v::AbstractSIMDVector{W,T1}, ::T2) where {W,T1,T2} = convert(T2,vprod(v))
@inline reduce_to_add(v::AbstractSIMDVector{W,T}, ::AbstractSIMDVector{W}) where {W,T} = v
@inline reduce_to_prod(v::AbstractSIMDVector{W,T}, ::AbstractSIMDVector{W}) where {W,T} = v
@inline reduce_to_prod(s::T, ::T) where {T<:Number} = v

# @inline reduced_all(u1::Unsigned, u1::Unsigned) where {W,T} = vall(v) & s
# @inline reduced_any(u::Unsigned, b::Bool) where {W,T} = b || vany(v)
# @inline reduced_all(v::AbstractSIMDVector, s::T) where {W,T} = vall(v) & s
# @inline reduced_any(v::AbstractSIMDVector, s::T) where {W,T} = s || vany(v)
@inline reduced_max(v::AbstractSIMDVector{W,T}, s::T) where {W,T} = max(vmaximum(v), s)
@inline reduced_min(v::AbstractSIMDVector{W,T}, s::T) where {W,T} = min(vminimum(v), s)
@inline reduced_max(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T} = vmax(v1, v2)
@inline reduced_min(v1::AbstractSIMDVector{W,T}, v2::AbstractSIMDVector{W,T}) where {W,T} = vmin(v1, v2)
@inline reduce_to_max(v::AbstractSIMDVector{W,T}, s::T) where {W,T} = vmaximum(v)
@inline reduce_to_min(v::AbstractSIMDVector{W,T}, s::T) where {W,T} = vminimum(v)
@inline reduce_to_max(v::AbstractSIMDVector{W,T}, ::AbstractSIMDVector{W,T}) where {W,T} = v
@inline reduce_to_min(v::AbstractSIMDVector{W,T}, ::AbstractSIMDVector{W,T}) where {W,T} = v

@inline reduced_max(v::AbstractSIMDVector{W,T1}, s::T2) where {W,T1,T2} = max(convert(T2,vmaximum(v)), s)
@inline reduced_min(v::AbstractSIMDVector{W,T1}, s::T2) where {W,T1,T2} = min(convert(T2,vminimum(v)), s)
@inline reduced_max(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1,T2} = vmax(vconvert(Vec{W,T2},v1), v2)
@inline reduced_min(v1::AbstractSIMDVector{W,T1}, v2::AbstractSIMDVector{W,T2}) where {W,T1,T2} = vmin(vconvert(Vec{W,T2},v1), v2)
@inline reduce_to_max(v::AbstractSIMDVector{W,T1}, s::T2) where {W,T1,T2} = convert(T2,vmaximum(v))
@inline reduce_to_min(v::AbstractSIMDVector{W,T1}, s::T2) where {W,T1,T2} = convert(T2,vminimum(v))
@inline reduce_to_max(v::AbstractSIMDVector{W,T1}, ::AbstractSIMDVector{W,T2}) where {W,T1,T2} = vconvert(Vec{W,T2}, v)
@inline reduce_to_min(v::AbstractSIMDVector{W,T1}, ::AbstractSIMDVector{W,T2}) where {W,T1,T2} = vconvert(Vec{W,T2}, v)

@inline vnmul(x,y) = vsub(vmul(x, y))
@inline vnsub(x,y) = vsub(y, x)
@inline vmul2(x) = vadd(x,x)
@inline vmul3(x::T) where {T <: Number} = Base.FastMath.mul_fast(T(3), x)
@inline vmul3(v::V) where {W, T, V <: AbstractSIMDVector{W,T}} = vmul(vbroadcast(V, T(3)), v)
@inline vadd1(x::T) where {T <: Number} = Base.FastMath.add_fast(x, one(x))
@inline vadd1(v::V) where {W, T, V <: AbstractSIMDVector{W,T}} = vadd(vbroadcast(V, one(T)), v)

@generated function addscalar(v::Vec{W,T}, s::T) where {W, T <: Integer}
    typ = "i$(8sizeof(T))"
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp zeroinitializer, $typ %1, i32 0")
    push!(instrs, "%v = add $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@generated function addscalar(v::Vec{W,T}, s::T) where {W, T <: Union{Float16,Float32,Float64}}
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    instrs = String[]
    push!(instrs, "%ie = insertelement $vtyp zeroinitializer, $typ %1, i32 0")
    push!(instrs, "%v = fadd fast $vtyp %0, %ie")
    push!(instrs, "ret $vtyp %v")
    quote
        $(Expr(:meta,:inline))
        Base.llvmcall( $(join(instrs,"\n")), NTuple{$W,Core.VecElement{$T}}, Tuple{NTuple{$W,Core.VecElement{$T}},$T}, v, s )
    end
end
@inline addscalar(v::SVec, s) = addscalar(extract_data(v), s)
@inline addscalar(s::T, v::Union{Vec{W,T},SVec{W,T}}) where {W,T} = addscalar(v, s)
@inline addscalar(a, b) = vadd(a, b)


text:
using Triangle
export Triangulation, RandomTriangulation
"""
`RandomTriangulation(n)` creates a random `n`-vertex planar graph
all of whose bounded faces are triangles.
"""
function RandomTriangulation(n::Int)::SimpleGraph{Int}
    pts = rand(n,2)
    G = Triangulation(pts)
    embed(G,:tutte)
    name(G,"Random Triangulation")
    return G
end

"""
`Triangle(point_list)` creates a graph that is the Delauney triangulation
of the `n`-by-`2` matrix `point_list` whose rows give the coordinates
of the points.
"""
function Triangulation(pts::Array{Float64,2})::SimpleGraph
    n,c = size(pts)
    G = IntGraph(n)
    idx = collect(1:n)
    TT = basic_triangulation(pts,idx)
    xy = Dict{Int,Array{Float64,1}}()
    for v = 1:n
        xy[v] = pts[v,:]
    end

    for T in TT
        a,b,c = T
        add!(G,a,b)
        add!(G,a,c)
        add!(G,b,c)
    end
    embed(G,xy)
    embed_rot(G)
    name(G,"Triangulation")
    return G
end

ext: jl


text:
function water_map(r::MersenneTwister)
    w = 5; h = 5
    map = Dict{String, Any}("height" => h, "width"  => w)
    sx, sy = rand(r, [(1, 1), (1, w), (h, 1), (h, w)])
    map["starts"] = [Dict("x"=>sx, "y"=>sy)]
    map
end

function water_meta()
    Dict{String, Any}("max_step" => 100)
end

function cov_terminate(e::Episode)
    terminate(e) || all(map(c->object(c) == "empty", e.grid.cells))
end

function coverage_fitness(cont_f::Function; seed::Int64=0)
    r = MersenneTwister(seed)
    map = water_map(r)
    g = Grid(map)
    for c in g.cells
        c.obj = type_to_float("food")
        c.color = color_to_float("green")
    end
    e = Episode(g; reward=food_reward, meta=water_meta())
    total_reward = run!(e, cont_f; terminate=cov_terminate)
    steps = terminate(e) ? e.meta["max_step"] : e.meta["step"]
    (total_reward + (1.0 - steps / e.meta["max_step"])) / (length(e.grid.cells) + 1)
end

function food_exit_map(seed::Int64=0)
    r = MersenneTwister(seed)
    map = water_map(r)
    ex, ey = rand(r, 2:4), rand(r, 2:4)
    map["exits"] = [Dict("x"=>ex, "y"=>ey)]
    map["objects"] = [Dict("type"=>"food", "color"=>"green", "x"=>ex, "y"=>ey)]
    map
end

function water_search_fitness(cont_f::Function; seed::Int64=0)
    e = Episode(Grid(food_exit_map(seed)); reward=food_reward, meta=water_meta())
    run!(e, cont_f)
    steps = terminate(e) ? e.meta["max_step"] : e.meta["step"]
    1 - steps / e.meta["max_step"]
end

function water_memorize_fitness(cont_f::Function; seed::Int64=0)
    map_function() = food_exit_map(seed)
    learn_fitness(cont_f, map_function, water_meta)
end


text:
"""
    function atp_ranking(
        players::Vector{Player},
        ranking_scheme_to_update::String
    )

Update the player's rankings using ATP Ranking system.

Arguments
=========
- `players`:
    an array of players in the simulation.
- `ranking_scheme_to_update`:
    the name of the ranking scheme to update for each Player.
"""

function atp_ranking(
    players::Vector{Player},
    ranking_scheme_to_update::String
)
    ordered_ranking = []
    for player in players
        point_total = 0
        top6_results = Vector{Int64}([])
        for (tournament_name, summary) in player.doubles_results
            if summary.tournament_type == "Grandslam" ||
               summary.tournament_type == "ATP1000"
                point_total +=
                    compute_points(summary.result, summary.tournament_type)
            else
                push!(
                    top6_results,
                    compute_points(summary.result, summary.tournament_type),
                )
            end
        end
        sort!(top6_results, rev = true)
        point_total += sum(top6_results[1:min(length(top6_results), 6)])
        player.rankings[ranking_scheme_to_update] = point_total
    end
    return nothing
end

"""
    function compute_points(
        result::String,
        tournament_type::String
    )::Int64

Retrieve the number of points corresponding to a player's performance in a tournament.

Arguments
=========
- `result`:
    the final round that a player finished in a tournament.
- `tournament_type`:
    the type of the tournament. Ex: Grandslam
"""

function compute_points(
    result::String,
    tournament_type::String
)::Int64
    if tournament_type == "Grandslam"
        return ATPGrandSlamDoublesPoints[result]
    elseif tournament_type == "ATPNitto"
        return ATPNittoFinalsPoints[result]
    elseif tournament_type == "ATP1000"
        return ATP1000DoublesPoints[result]
    elseif tournament_type == "ATP500"
        return ATP500DoublesPoints[result]
    elseif tournament_type == "ATP250"
        return ATP250DoublesPoints[result]
    elseif tournament_type == "ATPChallenger125"
        return ATPChallengerTour125DoublesPoints[result]
    elseif tournament_type == "ATPChallenger110"
        return ATPChallengerTour110DoublesPoints[result]
    elseif tournament_type == "ATPChallenger100"
        return ATPChallengerTour100DoublesPoints[result]
    elseif tournament_type == "ATPChallenger90"
        return ATPChallengerTour90DoublesPoints[result]
    elseif tournament_type == "ATPChallenger80"
        return ATPChallengerTour80DoublesPoints[result]
    elseif tournament_type == "ITFWorldTennisTour25k"
        return ITFWorldTennisTour25kPoints[result]
    elseif tournament_type == "ITFWorldTennisTour15k"
        return ITFWorldTennisTour15kPoints[result]
    end
    return throw(DomainError(tournament_type, "invalid tournament type"))
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
"""
    Grid2DTopology{M <: Mesh2D} <: AbstractTopology

An unstructured 2D conformal mesh of elements.
"""
struct Grid2DTopology{M} <: AbstractTopology
    mesh::M
end

domain(topology::Grid2DTopology{M}) where {M <: AbstractMesh} =
    topology.mesh.domain

function nlocalelems(topology::Grid2DTopology)
    return topology.mesh.nelems
end

eachslabindex(topology::Grid2DTopology) = 1:nlocalelems(topology)

function vertex_coordinates(topology::Grid2DTopology, elem::Integer)
    @assert 1 <= elem <= nlocalelems(topology)
    coords = topology.mesh.coordinates
    verts = view(topology.mesh.elem_verts, elem, :)
    dim = size(coords, 2)
    CT = Topologies.coordinate_type(topology)
    if dim == 2
        return (
            CT(coords[verts[1], 1], coords[verts[1], 2]),
            CT(coords[verts[2], 1], coords[verts[2], 2]),
            CT(coords[verts[3], 1], coords[verts[3], 2]),
            CT(coords[verts[4], 1], coords[verts[4], 2]),
        )
    else
        return (
            CT(coords[verts[1], 1], coords[verts[1], 2], coords[verts[1], 3]),
            CT(coords[verts[2], 1], coords[verts[2], 2], coords[verts[2], 3]),
            CT(coords[verts[3], 1], coords[verts[3], 2], coords[verts[3], 3]),
            CT(coords[verts[4], 1], coords[verts[4], 2], coords[verts[4], 3]),
        )
    end
end

function opposing_face(topology::Grid2DTopology, elem::Integer, face::Integer)
    @assert 1 <= elem <= nlocalelems(topology)
    @assert 1 <= face <= 4

    face_neighbors = topology.mesh.face_neighbors
    fc = topology.mesh.elem_faces[elem, face]

    if face_neighbors[fc, 1] == elem
        opelem = face_neighbors[fc, 3]
        opface = face_neighbors[fc, 4]
    elseif face_neighbors[fc, 3] == elem
        opelem = face_neighbors[fc, 1]
        opface = face_neighbors[fc, 2]
    else
        error("opposing_face: fatal error in connectivity information")
    end
    if opelem == 0
        opface = face
    end
    reversed = face_neighbors[5] == -1
    return opelem, opface, reversed
end

# InteriorFaceIterator
function Base.length(fiter::InteriorFaceIterator{T}) where {T <: Grid2DTopology}
    topology = fiter.topology
    mesh = topology.mesh
    loc = findfirst(mesh.boundary_tags .== 0) # interior tag is 0
    if isnothing(loc)
        return 0
    else
        return (
            mesh.face_boundary_offset[loc + 1] - mesh.face_boundary_offset[loc]
        )
    end
end

function Base.iterate(
    fiter::InteriorFaceIterator{T},
    faceno = 1,
) where {T <: Grid2DTopology}
    mesh = fiter.topology.mesh
    face_neighbors = mesh.face_neighbors
    face_boundary = mesh.face_boundary
    face_boundary_offset = mesh.face_boundary_offset
    boundary_tags = mesh.boundary_tags

    if boundary_tags[1] ≠ 0
        return nothing
    end
    st, en = face_boundary_offset[1], face_boundary_offset[2]
    nfcs = en - st
    if faceno > nfcs
        return nothing
    end
    face = face_boundary[faceno]
    reversed = face_neighbors[face, 5] == -1
    return tuple(face_neighbors[face, 1:4]..., reversed), faceno + 1
end

# BoundaryFaceIterator
function boundary_names(topology::Grid2DTopology)
    boundary_tag_names = topology.mesh.boundary_tag_names
    if boundary_tag_names[1] == :interior
        return tuple(boundary_tag_names[2:end]...)
    else
        return boundary_tag_names
    end
end

function boundary_tag(topology::Grid2DTopology, name::Symbol)
    boundary_tag_names = topology.mesh.boundary_tag_names
    boundary_tags = topology.mesh.boundary_tags
    location = findfirst(boundary_tag_names .== name)
    if isnothing(location)
        error("Invalid boundary name")
    else
        return boundary_tags[location]
    end
end

function boundaries(topology::Grid2DTopology)
    boundary_tags = topology.mesh.boundary_tags
    if boundary_tags[1] == 0
        return tuple(boundary_tags[2:end]...)
    else
        return tuple(boundary_tags...)
    end
end

function Base.length(
    bfiter::BoundaryFaceIterator{T},
) where {T <: Grid2DTopology}
    boundary = bfiter.boundary
    topology = bfiter.topology
    face_boundary_offset = topology.mesh.face_boundary_offset
    boundary_tags = topology.mesh.boundary_tags
    location = findfirst(boundary_tags .== boundary)
    if isnothing(location)
        return 0
    end
    return face_boundary_offset[location + 1] - face_boundary_offset[location]
end

function Base.iterate(
    bfiter::BoundaryFaceIterator{T},
    faceno = 1,
) where {T <: Grid2DTopology}
    boundary = bfiter.boundary
    topology = bfiter.topology
    mesh = topology.mesh
    face_neighbors = mesh.face_neighbors
    face_boundary = mesh.face_boundary
    face_boundary_offset = mesh.face_boundary_offset
    boundary_tags = mesh.boundary_tags
    location = findfirst(boundary_tags .== boundary)

    if isnothing(location)
        return nothing
    end
    st, en = face_boundary_offset[location:(location + 1)]
    nbndry = en - st
    if faceno < 1 || faceno > nbndry
        return nothing
    end
    face = face_boundary[faceno + st - 1]
    if face_neighbors[face, 1] == 0
        return tuple(face_neighbors[face, 3:4]...), faceno + 1
    else
        return tuple(face_neighbors[face, 1:2]...), faceno + 1
    end
end

# VertexIterator
function Base.length(viter::VertexIterator{T}) where {T <: Grid2DTopology}
    return length(viter.topology.mesh.unique_verts)
end

function Base.iterate(
    viter::VertexIterator{T},
    uvertno = 1,
) where {T <: Grid2DTopology}
    topology = viter.topology
    nuverts = length(viter.topology.mesh.unique_verts)
    if uvertno > 0 && uvertno ≤ nuverts
        return Vertex(topology, uvertno), uvertno + 1
    else
        return nothing
    end
end

# Vertex
function Base.length(vertex::Vertex{T}) where {T <: Grid2DTopology}
    mesh = vertex.topology.mesh
    vtno = vertex.num
    offset = mesh.uverts_offset
    return div(offset[vtno + 1] - offset[vtno], 2)
end

function Base.iterate(vertex::Vertex{T}, iterno = 1) where {T <: Grid2DTopology}
    # iterator of (element, vertnum) that share global vertex
    mesh = vertex.topology.mesh
    vtno = vertex.num
    offset = mesh.uverts_offset
    uverts_conn = mesh.uverts_conn
    st, en = offset[vtno], offset[vtno + 1]
    niter = div(en - st, 2)

    if iterno < 1 || iterno > niter
        return nothing
    end
    loc = (iterno - 1) * 2 + st
    return tuple(uverts_conn[loc:(loc + 1)]...), iterno + 1
end


text:
"""
    get_mode(A, z, pseudo, verbosity)

Compute the (pseudo)mode of `A` associated with `z`.
"""
function get_mode(A, z, pseudo, verbosity)
    if pseudo
        n = size(A, 2)
        q = randn(n) + im*randn(n)
        normalize!(q)
        niter = 20

        SS, q = psmode_inv_lanczos(A, q, z, 1e-15, niter)
    else
        SS, q = oneeigcond(A,z, verbosity)
    end
    return SS, q
end

"""
    get_mtxpowersnorm(A, nmax)

Compute ``|A^k|`` for k up to `nmax`.
"""
function get_mtxpowersnorm(A, nmax)
    transient = Vector{Float64}(undef, nmax+1)
    transient[1] = 1.0
    (min_tp, max_tp) = (Inf, -Inf)
    mat = copy(A)
    i = 2
    transient[i] = norm(mat)
    stop = false
    while !stop && i < nmax + 1
        mat = A * mat
        i += 1
        transient[i] = norm(mat)
        min_tp = min(min_tp, transient[i])
        max_tp = max(max_tp, transient[i])
        if (max_tp > 1e130) || ((min_tp < 1e-130) && (min_tp != 0.0))
            @warn("stopping: bounds going out of range")
            stop = true
        end
        if min_tp == 0.0
            @warn("stopping: exactly zero matrix")
            stop = true
        end
    end
    return 0:(i-1), transient[1:i]
end

"""
    get_mtxexpnorm(A, dt, nmax)

Compute ``|e^{A*k*dt}|`` for k up to `nmax`.
"""
function get_mtxexpnorm(A, dt, nmax)
    eAdt = exp(dt*A)
    ns, transient = get_mtxpowersnorm(eAdt, nmax)
    return dt*ns, transient
end


text:
module ParallelRun

using Base.Threads
using Base.Dates

import FITSIO
import JLD
import WCS
using DataStructures

import ..Celeste: detect_sources, BoundingBox, Config
import ..Log
using ..Model
import ..PSF
import ..Coordinates: angular_separation, match_coordinates
import ..MCMC

using ..DeterministicVI
using ..DeterministicVI.ConstraintTransforms: ConstraintBatch, DEFAULT_CHUNK
using ..DeterministicVI.ElboMaximize: ElboConfig, maximize!

include("partition.jl")

"""
Return a Cyclades partitioning or an equal partitioning of the target
sources.
"""
function partition_box(npartitions::Int, target_sources::Vector{Int},
                       neighbor_map::Dict{Int,Vector{Int}}, ea_vec;
                       cyclades_partition=true,
                       batch_size=7000)
    if cyclades_partition
        #return partition_cyclades(npartitions, target_sources,
        #                          neighbor_map,
        #                          batch_size=batch_size)
        #return partition_cyclades_dynamic(target_sources,
        #                                  neighbor_map,
        #                                  batch_size=batch_size)
	return partition_cyclades_dynamic_auto_batchsize(target_sources, neighbor_map, ea_vec)
    else
        return partition_equally(npartitions, length(target_sources))
    end
end

function estimate_time(patches)
    sum(sum(p.active_pixel_bitmap) for p in patches)
end

function load_balance_across_threads(n_threads, times)
    ts = [0 for i=1:n_threads]
    for t in times
        minimum,index_min = findmin(ts)
        ts[index_min] += t
    end
    return ts
end

"""
Partitions sources via the cyclades algorithm. Finds the batch size which returns most balanced CC distribution.
- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources
"""
function partition_cyclades_dynamic_auto_batchsize(target_sources::Vector{Int}, neighbor_map::Dict{Int,Vector{Int}}, ea_vec)
    n_threads = nthreads()

    # Sample batch sizes at intervals
    n_to_sample = 100
    stepsize = max(1, trunc(Int, length(target_sources) / n_to_sample))

    # Best load imbalance of the batch sizes to test
    best_score = Inf
    best_result = Inf
    best_batch_size = -Inf

    for batch_size_to_use = 1 : stepsize : length(target_sources)+1
        ccs = partition_cyclades_dynamic(target_sources, neighbor_map, batch_size=batch_size_to_use)
        score = 0
        for batch in ccs
            # Find average load imbalance within the batch as a percentage
	    times = [sum([estimate_time(ea_vec[source_index].patches) for source_index in component]) for component in batch]
            times = load_balance_across_threads(n_threads, times)
            estimated_imbalance = mean(maximum(times) - times)
	    score += estimated_imbalance
        end
        if score <= best_score
            best_result = ccs
            best_score = score
            best_batch_size = batch_size_to_use
        end
    end
    for batch in best_result
        sizes = [length(x) for x in batch]
    end
    best_result
end

"""
Pre-allocate elbo args variational parameters; configurations need to
be persisted across calls to `maximize!()` so that location constraints
don't shift from their initial position.
"""
function setup_vecs(n_sources::Int, target_sources::Vector{Int},
                    catalog::Vector{CatalogEntry})
    ts_vp = Dict{Int64, Array{Float64}}()
    for ts in target_sources
        cat = catalog[ts]
        ts_vp[ts] = generic_init_source(cat.pos)
    end

    ea_vec = Vector{ElboArgs}(n_sources)
    vp_vec = Vector{VariationalParams{Float64}}(n_sources)
    cfg_vec = Vector{ElboConfig{DEFAULT_CHUNK,Float64}}(n_sources)

    return ea_vec, vp_vec, cfg_vec, ts_vp
end

"""
Uses multiple threads on one node to fit the Celeste
model over numerous iterations.

catalog - the catalog of light sources
target_sources - light sources to optimize
neighbor_map - light_source id -> neighbor light_source ids

cyclades_partition - use the cyclades algorithm to partition into non conflicting batches for updates.
batch_size - size of a single batch of sources for updates
within_batch_shuffling - whether or not to process sources within a batch randomly
joint_inference_terminate - whether to terminate once sources seem to be stable

Returns:

- Vector of OptimizedSource results
"""
function one_node_joint_infer(catalog, patches, target_sources, neighbor_map,
                              images;
                              cyclades_partition::Bool=true,
                              batch_size::Int=7000,
                              within_batch_shuffling::Bool=true,
                              config=Config())
    # Seed random number generator to ensure the same results per run.
    srand(42)

    n_threads = nthreads()

    # Partition the sources
    n_sources = length(target_sources)
    Log.info("Optimizing $(n_sources) sources")

    #Log.info(batched_connected_components)

    Log.info("Done assigning sources to threads for processing")

    ea_vec, vp_vec, cfg_vec, target_source_variational_params =
            setup_vecs(n_sources, target_sources, catalog)

    # Initialize elboargs in parallel
    thread_initialize_sources_assignment::Vector{Vector{Vector{Int64}}} = partition_equally(n_threads, n_sources)

    initialize_elboargs_sources!(config, ea_vec, vp_vec, cfg_vec, thread_initialize_sources_assignment,
                                 catalog, patches, target_sources, neighbor_map, images,
                                 target_source_variational_params)

    #thread_sources_assignment = partition_box(n_threads, target_sources,
    #                                  neighbor_map;
    #                                  cyclades_partition=cyclades_partition,
    #                                  batch_size=batch_size)
    batched_connected_components = partition_box(n_threads, target_sources,
                                                 neighbor_map, ea_vec;
                                                 cyclades_partition=cyclades_partition,
                                                 batch_size=batch_size)

    # Process sources in parallel
    #process_sources!(images, ea_vec, vp_vec, cfg_vec,
    #                 thread_sources_assignment,
    #                 n_iters, within_batch_shuffling)
    niters = config.num_joint_vi_iters
    process_sources_dynamic!(images, ea_vec, vp_vec, cfg_vec,
                             batched_connected_components,
                             niters, within_batch_shuffling)

    # Return add results to vector
    results = OptimizedSource[]

    for i = 1:n_sources
        entry = catalog[target_sources[i]]
        is_sky_bad = bad_sky(entry, images)
        result = OptimizedSource(entry.pos[1], entry.pos[2], vp_vec[i][1],
                                 is_sky_bad)
        push!(results, result)
    end

    show_pixels_processed()

    results
end

function initialize_elboargs_sources!(config::Config, ea_vec, vp_vec, cfg_vec,
                                      thread_initialize_sources_assignment,
                                      catalog, patches, target_sources, neighbor_map, images,
                                      target_source_variational_params;
                                      termination_callback=nothing)
    Threads.@threads for i in 1:nthreads()
        try
            for batch in 1:length(thread_initialize_sources_assignment[i])
                for source_index in thread_initialize_sources_assignment[i][batch]
                    init_elboargs(config, source_index, catalog, patches, target_sources,
                                  neighbor_map, images, ea_vec, vp_vec, cfg_vec,
                                  target_source_variational_params;
                                  termination_callback=termination_callback)
                end
            end
        catch ex
            Log.exception(ex)
            rethrow()
        end
    end
end

"""
Initialize elbo args for the specified target source.
"""
function init_elboargs(config::Config,
                       ts::Int,
                       catalog::Vector{CatalogEntry},
                       patches::Matrix{ImagePatch},
                       target_sources::Vector{Int},
                       neighbor_map::Dict{Int,Vector{Int}},
                       images::Vector{<:Image},
                       ea_vec::Vector{ElboArgs},
                       vp_vec::Vector{VariationalParams{Float64}},
                       cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                       ts_vp::Dict{Int64,Array{Float64}};
                       termination_callback=nothing)
    try
        entry_id = target_sources[ts]
        entry = catalog[entry_id]
        neighbor_ids = neighbor_map[entry_id]
        neighbors = catalog[neighbor_ids]
        cat_local = vcat([entry], neighbors)
        ids_local = vcat([entry_id], neighbor_ids)

        # Limit patches to just the active source and its neighbors.
        patches = patches[ids_local, :]

        # Load vp with shared target source params, and also vp
        # that doesn't share target source params
        vp = Vector{Float64}[haskey(ts_vp, x) ?
                             ts_vp[x] :
                             catalog_init_source(catalog[x])
                             for x in ids_local]
        ea = ElboArgs(images, patches, [1])

        ea_vec[ts] = ea
        vp_vec[ts] = vp
        cfg_vec[ts] = ElboConfig(ea, vp;
                termination_callback=termination_callback)
    catch exc
        if is_production_run || nthreads() > 1
            Log.exception(exc)
        else
            rethrow()
        end
    end
end

function process_sources!(images::Vector{<:Model.Image},
                          ea_vec::Vector{ElboArgs},
                          vp_vec::Vector{VariationalParams{Float64}},
                          cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                          thread_sources_assignment::Vector{Vector{Vector{Int64}}},
                          n_iters::Int,
                          within_batch_shuffling::Bool)
    n_threads::Int = nthreads()
    n_batches::Int = length(thread_sources_assignment[1])
    for iter in 1:n_iters
        # Process every batch of every iteration. We do the batches on the outside
        # Since there is an implicit barrier after the inner threaded for loop below.
        # We want this barrier because there may be conflict _between_ Cyclades batches.
        for batch in 1:n_batches
            process_sources_elapsed_times::Vector{Float64} = Vector{Float64}(n_threads)
            # Process every batch of every iteration with n_threads
            Threads.@threads for i in 1:n_threads
                try
                    tic()
                    process_sources_kernel!(ea_vec, vp_vec, cfg_vec,
                                            thread_sources_assignment[i::Int][batch::Int],
                                            within_batch_shuffling)
                    process_sources_elapsed_times[i::Int] = toq()
                catch exc
                    Log.exception(exc)
                    rethrow()
                end
            end
            Log.info("Batch $(batch) - $(process_sources_elapsed_times)")
        end
    end
end

function process_sources_dynamic!(images::Vector{<:Model.Image},
                                  ea_vec::Vector{ElboArgs},
                                  vp_vec::Vector{VariationalParams{Float64}},
                                  cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                                  thread_sources_assignment::Vector{Vector{Vector{Int64}}},
                                  n_iters::Int,
                                  within_batch_shuffling::Bool)
    Log.info("Processing with dynamic connected components load balancing")

    n_threads::Int = nthreads()
    n_batches::Int = length(thread_sources_assignment)

    l = SpinLock()

    total_idle_time = 0
    total_sum_of_thread_times = 0

    for iter in 1:n_iters

        # Process every batch of every iteration. We do the batches on the outside
        # Since there is an implicit barrier after the inner threaded for loop below.
        # We want this barrier because there may be conflict _between_ Cyclades batches.
        for batch in 1:n_batches

            connected_components_index = 1
            process_sources_elapsed_times::Vector{Float64} = Vector{Float64}(n_threads)

            # Process every batch of every iteration with n_threads
            Threads.@threads for i in 1:n_threads
                try
                    tic()
                    while true
                        cc_index = -1

                        #ccall(:jl_,Void,(Any,), "this is thread number $(Threads.threadid()) $(cc_index)")

                        lock(l)
                        cc_index = connected_components_index
                        connected_components_index += 1
                        unlock(l)

                        if cc_index > length(thread_sources_assignment[batch::Int])
                            break
                        end

                        #ccall(:jl_,Void,(Any,), "this is thread number $(Threads.threadid()) $(cc_index) popped")

                        process_sources_kernel!(ea_vec, vp_vec, cfg_vec,
                                                thread_sources_assignment[batch::Int][cc_index],
                                                within_batch_shuffling)
                    end
                    process_sources_elapsed_times[i::Int] = toq()
                catch exc
                    Log.exception(exc)
                    rethrow()
                end
            end
            Log.info("Batch $(batch) - $(process_sources_elapsed_times)")
	    avg_thread_idle_time = mean(maximum(process_sources_elapsed_times) - process_sources_elapsed_times)
	    maximum_thread_time = maximum(process_sources_elapsed_times)
            idle_percent = 100.0 * (avg_thread_idle_time / maximum_thread_time)
            Log.info("Batch $(batch) avg threads idle: $(round(Int, idle_percent))% ($(avg_thread_idle_time) / $(maximum_thread_time))")
	    total_idle_time += sum(maximum(process_sources_elapsed_times) - process_sources_elapsed_times)
            total_sum_of_thread_times += sum(process_sources_elapsed_times)
        end
    end
    Log.info("Total idle time: $(round(Int, total_idle_time)), Total sum of threads times: $(round(Int, total_sum_of_thread_times))")
end

# Process partition of sources. Multiple threads call this function in parallel.
function process_sources_kernel!(ea_vec::Vector{ElboArgs},
                                 vp_vec::Vector{VariationalParams{Float64}},
                                 cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                                 source_assignment::Vector{Int64},
                                 within_batch_shuffling::Bool)
    try
        # Shuffle the source assignments within each batch of each process.
        # This is disabled by default because it ruins the deterministic outcome
        # required by the test cases.
        if within_batch_shuffling
            shuffle!(source_assignment)
        end

	#tic()
        for i in source_assignment
            maximize!(ea_vec[i], vp_vec[i], cfg_vec[i])
        end
	#ccall(:jl_,Void,(Any,), "this is thread number $(Threads.threadid()) took $(toq()) seconds to process batch with $(length(source_assignment)) sources")
    catch ex
        if is_production_run || nthreads() > 1
            Log.exception(ex)
        else
            rethrow(ex)
        end
    end
end

function get_pixels_processed()
    n_active = 0
    n_inactive = 0
    for elbo_vars in DeterministicVI.ElboMaximize.ELBO_VARS_POOL
        n_active += elbo_vars.active_pixel_counter[]
        n_inactive += elbo_vars.inactive_pixel_counter[]
        elbo_vars.active_pixel_counter[] = 0
        elbo_vars.inactive_pixel_counter[] = 0
    end
    return (n_active, n_inactive)
end

# Count and display the number of active and inactive pixels processed.
function show_pixels_processed()
    n_active, n_inactive = get_pixels_processed()
    Log.info("$(Time(now())): (active,inactive) pixels processed: \($n_active,$n_inactive\)")
end

# In production mode, rather the development mode, always catch exceptions
const is_production_run = haskey(ENV, "CELESTE_PROD") && ENV["CELESTE_PROD"] != ""

# ------
# optimization

# optimization result container
struct OptimizedSource
    init_ra::Float64
    init_dec::Float64
    vs::Vector{Float64}
    is_sky_bad::Bool
end

# The sky intensity does not appear to match the true background intensity
# for some light sources. This function flags light sources whose infered
# brightness may be inaccurate because the background intensity estimates
# are off.
function bad_sky(ce, images)
    # The 'i' band is a pretty bright one, so I used it here.
    img_index = findfirst(img -> img.b == 4, images)
    if img_index < 0
        return false
    end

    img = images[img_index]

    # Get sky at location of the object
    pixel_center = WCS.world_to_pix(img.wcs, ce.pos)
    h = max(1, min(round(Int, pixel_center[1]), size(img.pixels, 1)))
    w = max(1, min(round(Int, pixel_center[2]), size(img.pixels, 2)))

    claimed_sky = img.sky[h, w] * img.nelec_per_nmgy[h]

    # Determine background in a 50-pixel radius box around the object
    box = Model.box_around_point(img.wcs, ce.pos, 50.0)
    h_range, w_range = Model.clamp_box(box, (img.H, img.W))
    observed_sky = median(filter(!isnan, img.pixels[h_range, w_range]))

    # A 5 photon-per-pixel disparity can really add up if the light sources
    # covers a lot of pixels.
    return (claimed_sky + 5) < observed_sky
end

"""
Optimize the `s`th element in `catalog`.
Used only for one_node_single_infer, not one_node_joint_infer.
"""
function process_source(config::Config,
                        s::Int,
                        catalog::Vector{CatalogEntry},
                        patches::Matrix{ImagePatch},
                        neighbor_ids::Vector{Int},
                        images::Vector{<:Image})
    neighbors = catalog[neighbor_ids]
    if length(neighbors) > 100
        msg = string("object at RA, Dec = $(entry.pos) has an excessive",
                     "number ($(length(neighbors))) of neighbors")
        Log.warn(msg)
    end

    entry = catalog[s]
    cat_local = vcat([entry], neighbors)

    # Limit patches to just the active source and its neighbors.
    idxs = vcat([s], neighbor_ids)
    patches = patches[idxs, :]

    vp = DeterministicVI.init_sources([1], cat_local)
    ea = DeterministicVI.ElboArgs(images, patches, [1])

    tic()
    f_evals, max_f, max_x, nm_result = DeterministicVI.ElboMaximize.maximize!(ea, vp)
    Log.info("#$(s) at ($(entry.pos[1]), $(entry.pos[2])): $(toq()) secs")

    vs_opt = vp[1]
    is_sky_bad = bad_sky(entry, images)
    return OptimizedSource(entry.pos[1], entry.pos[2], vs_opt, is_sky_bad)
end

"""
Run MCMC to process a source.  Returns
"""
function process_source_mcmc(config::Config,
                             s::Int,
                             catalog::Vector{CatalogEntry},
                             patches::Matrix{ImagePatch},
                             neighbor_ids::Vector{Int},
                             images::Vector{<:Image};
                             use_ais::Bool=true)
    # subselect source, select active source and neighbor set
    entry = catalog[s]
    neighbors = catalog[neighbor_ids]
    if length(neighbors) > 100
        msg = string("objid $(entry.objid) [ra: $(entry.pos)] has an excessive",
                     "number ($(length(neighbors))) of neighbors")
        Log.warn(msg)
    end

    # Limit patches to just the active source and its neighbors.
    idxs = vcat([s], neighbor_ids)
    patches = patches[idxs, :]

    # create smaller images for the MCMC sampler to use
    patch_images = [MCMC.patch_to_image(patches[1, i], images[i])
                    for i in 1:length(images)]

    # render a background image on the active source (first in list)
    background_images = [MCMC.render_patch_nmgy(patch_images[i], patches[1, i], neighbors)
                         for i in 1:length(patch_images)]

    # run mcmc sampler on this image/patch/background initialized at entry
    if use_ais
        mcmc_results = MCMC.run_ais(entry, patch_images, patches, background_images;
            num_temperatures=config.num_ais_temperatures,
            num_samples=config.num_ais_samples)
    else
        mcmc_results = MCMC.run_mcmc(entry, patch_images, patches, background_images)
    end

    # summary
    return mcmc_results
end

function one_node_single_infer(catalog::Vector{CatalogEntry},
                               patches::Matrix{ImagePatch},
                               target_sources::Vector{Int},
                               neighbor_map::Dict{Int,Vector{Int}},
                               images::Vector{<:Image};
                               config=Config(),
                               do_vi=true)
    curr_source = 1
    last_source = length(target_sources)
    sources_lock = SpinLock()
    results_lock = SpinLock()
    if do_vi
        results = OptimizedSource[]
    else
        results = []
    end

    # iterate over sources
    function process_sources()
        tid = threadid()

        while true
            lock(sources_lock)
            ts = curr_source
            curr_source += 1
            unlock(sources_lock)
            if ts > last_source
                break
            end

            s = target_sources[ts]
            try
                if do_vi
                    result = process_source(config, s, catalog, patches,
                                            neighbor_map[s], images)
                else
                    result = process_source_mcmc(config, s, catalog, patches,
                                                 neighbor_map[s], images)
                end

                lock(results_lock)
                push!(results, result)
                unlock(results_lock)
            catch ex
                if is_production_run || nthreads() > 1
                    Log.exception(ex)
                else
                    rethrow(ex)
                end
            end
        end
    end

    if nthreads() == 1
        process_sources()
    else
        ccall(:jl_threading_run, Void, (Any,), Core.svec(process_sources))
        Log.LEVEL[] >= Log.INFO && ccall(:jl_threading_profile, Void, ())
    end

    return results
end

function _infer_box(images::Vector{<:Image}, catalog::Vector{CatalogEntry},
                    patches::Matrix{ImagePatch}, box::BoundingBox;
                    method=:joint_vi, config=Config())
    # Get indices of entries in the RA/Dec range of interest.
    # (Some images can have regions that are outside the box, so not
    # all sources are necessarily in the box.)
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_ids = find(entry_in_range, catalog)

    # find objects with patches overlapping
    neighbor_map = Dict(id=>Model.find_neighbors(patches, id)
                        for id in target_ids)

    if method == :joint_vi
        results = one_node_joint_infer(catalog, patches, target_ids,
                                       neighbor_map, images; config=config)
    elseif method == :single_vi
        results = one_node_single_infer(catalog, patches, target_ids,
                                        neighbor_map, images; do_vi=true,
                                        config=config)
    elseif method == :mcmc
        results = one_node_single_infer(catalog, patches, target_ids,
                                        neighbor_map, images; do_vi=false,
                                        config=config)
    else
        error("unknown method: $method")
    end
    return results
end

"""
    infer_box(images, box; options...)

Detect objects in images and run inference on all sources
within the bounding box. If `catalog` is given, run inference on objects in
`catalog` that fall within the bounding box.

Objects outside the bounding box (either detected or given)
may be used as "neighbors".
"""
function infer_box(images::Vector{<:Image}, box::BoundingBox;
                   catalog::Union{Void, Vector{CatalogEntry}}=nothing,
                   method=:joint_vi, config=Config())
    tic()
    Log.info("processing box $(box.ramin), $(box.ramax), $(box.decmin), ",
             "$(box.decmax) with $(nthreads()) threads")

    if catalog === nothing
        cat, patches = detect_sources(images)
    else
        cat = catalog
        patches = Model.get_sky_patches(images, cat)
    end

    results = _infer_box(images, cat, patches, box;
                         method=method, config=config)

    Log.info("infer_box for $box took $(toq()) seconds")

    return results
end

end


text:
module Semirings

# SimpleTraits is only used for a secondary hierarchy of semirings (commutative, idempotent, etc)

using SimpleTraits

# imports from base

import Base: convert, eltype, promote_rule
import Base: +, *, ==, -, /, inv, conj, ≈
import Base: zero, one, iszero, isone

# exporting types, new functions, and new symbols

export AbstractSemiringElement
export ⊕, ⊙, star, mulinv, addinv, val
export issamesemiring, semiring

# traits
export IsCommutative, iscommutative, IsIdempotent, isidempotent, IsStar, isstar

# export implementations

export MaxPlusSemiringElement, MinPlusSemiringElement, RealSemiringElement, BooleanSemiringElement

# includes

include("interface.jl")
include("traits.jl")
include("core.jl")
include("implementations.jl")

end # module


text:
# Contexts are execution states for the core LLVM IR system.

export Context, dispose, GlobalContext

@checked struct Context
    ref::API.LLVMContextRef
end
reftype(::Type{Context}) = API.LLVMContextRef

function Context()
    ctx = Context(API.LLVMContextCreate())
    _install_handlers(ctx)
    ctx
end

dispose(ctx::Context) = API.LLVMContextDispose(ref(ctx))

function Context(f::Core.Function)
    ctx = Context()
    try
        f(ctx)
    finally
        dispose(ctx)
    end
end

GlobalContext() = Context(API.LLVMGetGlobalContext())

## wrapper exception type

export LLVMException

struct LLVMException <: Exception
    info::String
end

function Base.showerror(io::IO, err::LLVMException)
    @printf(io, "LLVM error: %s", err.info)
end

## diagnostics

export DiagnosticInfo, severity, message

@checked struct DiagnosticInfo
    ref::API.LLVMDiagnosticInfoRef
end
reftype(::Type{DiagnosticInfo}) = API.LLVMDiagnosticInfoRef

severity(di::DiagnosticInfo) = API.LLVMGetDiagInfoSeverity(ref(di))
message(di::DiagnosticInfo) = unsafe_string(API.LLVMGetDiagInfoDescription(ref(di)))

## handlers

function handle_diagnostic(diag_ref::API.LLVMDiagnosticInfoRef, args::Ptr{Cvoid})
    di = DiagnosticInfo(diag_ref)
    @assert args == C_NULL

    sev = severity(di)
    msg = message(di)

    if sev == API.LLVMDSError
        # NOTE: it might be more true to the API to just report an error here,
        #       and require callers to verify whether operations succeeded,
        #       but throwing here fails faster with better backtraces.
        throw(LLVMException(msg))
    elseif sev == API.LLVMDSWarning
        @warn msg
    elseif sev == API.LLVMDSRemark || sev == API.LLVMDSNote
        @debug msg
    else
        error("unknown diagnostic severity level $sev")
    end

    return nothing
end

function yield_callback(ctx_ref::API.LLVMContextRef, args::Ptr{Cvoid})
    ctx = Context(ctx_ref)
    @assert args == C_NULL

    # TODO: is this allowed? can we yield out of an active `ccall`?
    yield()
end

function _install_handlers(ctx::Context)
    # set yield callback
    callback = @cfunction(yield_callback, Cvoid, (Context, Ptr{Cvoid}))
    # NOTE: disabled until proven safe
    #API.LLVMContextSetYieldCallback(ref(ctx), callback, C_NULL)

    # set diagnostic callback
    handler = @cfunction(handle_diagnostic, Cvoid, (API.LLVMDiagnosticInfoRef, Ptr{Cvoid}))
    API.LLVMContextSetDiagnosticHandler(ref(ctx), handler, C_NULL)

    return nothing
end

function handle_error(reason::Cstring)
    throw(LLVMException(unsafe_string(reason)))
end

function _install_handlers()
    handler = @cfunction(handle_error, Cvoid, (Cstring,))

    # NOTE: LLVM doesn't support re-installing the error handler, which happens when we
    #       `reload("LLVM")`, so reset it instead. This isn't correct, as the installed
    #       handler might not have been ours. Ideally we'd have module finalizers...
    API.LLVMResetFatalErrorHandler()
    API.LLVMInstallFatalErrorHandler(handler)
end


text:
_default_basesize(n::Integer) = max(1, cld(n, Threads.nthreads()))
default_basesize(xs) = _default_basesize(length(xs))
const DEFAULT_SMALLSIZE = 64

Base.@kwdef struct ParallelStableMSDRadixSortAlg{Alg,BaseSize} <: RadixSortAlgorithm
    smallsort::Alg = Base.Sort.DEFAULT_STABLE
    smallsize::Int = DEFAULT_SMALLSIZE
    basesize::BaseSize = nothing  # lazily determined
end

Base.@kwdef struct StableMSDRadixSortAlg{Alg} <: RadixSortAlgorithm
    smallsort::Alg = Base.Sort.DEFAULT_STABLE
    smallsize::Int = DEFAULT_SMALLSIZE
end

const ParallelStableMSDRadixSort = ParallelStableMSDRadixSortAlg()
const StableMSDRadixSort = StableMSDRadixSortAlg()

const ParallelMSDRadixSort =
    ParallelStableMSDRadixSortAlg(smallsort = Base.Sort.DEFAULT_UNSTABLE)
const MSDRadixSort = StableMSDRadixSortAlg(smallsort = Base.Sort.DEFAULT_UNSTABLE)

_compose(::typeof(identity), ::typeof(identity)) = identity
_compose(::typeof(identity), f) = f
_compose(f, ::typeof(identity)) = f
_compose(::typeof(-), ::typeof(-)) = identity  # OK?
_compose(f, g) = f ∘ g

# TODO: reversing order should be done after the projection
_ord_by(o) = error("unsupported ordering: $o")
_ord_by(::ForwardOrdering) = identity
_ord_by(::ReverseOrdering{ForwardOrdering}) = -  # TODO; fix it for strings
_ord_by(order::Order.By) = _compose(_ord_by(order.order), order.by)

function Base.sort!(
    v::AbstractVector,
    lo::Integer,
    hi::Integer,
    a::ParallelStableMSDRadixSortAlg,
    o::Ordering,
)
    xs = view(v, lo:hi)
    _stablemsd!(
        xs,
        copy(xs),
        _ord_by(o),
        static(1),
        something(a.basesize, _default_basesize(hi - lo + 1)),
        a.smallsize,
        a.smallsort,
        Val(true),  # xs_mutable
    )
    return xs
end

function Base.sort!(
    v::AbstractVector,
    lo::Integer,
    hi::Integer,
    a::StableMSDRadixSortAlg,
    o::Ordering,
)
    xs = view(v, lo:hi)
    _stablemsd_seq!(
        xs,
        copy(xs),
        _ord_by(o),
        static(1),
        a.smallsize,
        a.smallsort,
        Val(true),
    )
    return xs
end

# TODO: NaN ordering correctness
Base.Sort.Float.fpsort!(v::AbstractVector, a::RadixSortAlgorithm, o::Ordering) =
    sort!(v, firstindex(v), lastindex(v), a, o)

# https://cs.stackexchange.com/a/96420
# https://www.youtube.com/watch?v=zqs87a_7zxw

function _stablemsd!(
    ys,
    xs,
    by::BY,
    ibyte,
    basesize,
    smallsize,
    smallsort,
    xs_mutable::Union{Val{true},Val{false}} = Val(false),
) where {BY}
    @assert firstindex(xs) == 1
    @assert firstindex(ys) == 1

    will_be_allpadded(eltype(xs), by, ibyte) && return ys

    if length(xs) <= basesize
        return _stablemsd_seq!(ys, xs, by, ibyte, smallsize, smallsort, xs_mutable)
    end

    # Split the input into at most `nthreads` parts since there's no recursion
    # in `_countmsd!` to benefit from constructive cache sharing:
    chunksize = max(basesize, cld(length(xs), Threads.nthreads()))

    chunks = Iterators.partition(xs, chunksize)
    allcounts = alloc!(COUNTS_POOL, length(chunks))
    allpaddeds = Vector{Bool}(undef, length(chunks))
    let allcounts = allcounts
        @sync for (i, xs) in enumerate(chunks)
            @spawn begin
                fill!(allcounts[i], 0)
                (_, allpaddeds[i]) = _countmsd!(allcounts[i], xs, by, ibyte)
            end
        end
    end
    alloffsets = allcounts
    @DBG allcounts = map(copy, alloffsets)
    all(allpaddeds) && return ys

    # TODO: parallel prefix
    # exclusive scan
    foldl(Iterators.product(alloffsets, 1:256); init = 0) do acc, (offsets, i)
        @inbounds (acc, offsets[i]) = (acc + offsets[i], acc)
        return acc
    end

    @sync for (idx, exclusive_offsets) in
              zip(Iterators.partition(1:length(xs), chunksize), alloffsets)
        @spawn _scatter!(ys, xs, idx, exclusive_offsets, by, ibyte)
    end
    @DBG @check vec(reduce(hcat, alloffsets)') == cumsum(vec(reduce(hcat, allcounts)'))

    free!(COUNTS_POOL, @view alloffsets[1:end-1])
    _spawn_foreach_remaining_subrange(alloffsets[end], basesize) do idx
        ys_chunk = view(ys, idx)
        if length(idx) <= smallsize
            sort!(ys_chunk, smallsort, Base.ord(isless, by, nothing, Base.Forward))
        else
            xs_chunk = if xs_mutable === Val(true)
                view(xs, idx)
            else
                similar(ys_chunk)
            end
            _stablemsd!(
                ys_chunk,
                copyto!(xs_chunk, ys_chunk),
                by,
                static_if_leq(@stat(ibyte + 1), static(8)),
                basesize,
                smallsize,
                smallsort,
                Val(true),
            )
        end
        return
    end
    free!(COUNTS_POOL, @view alloffsets[end:end])

    return ys
end

function _stablemsd_seq!(
    ys,
    xs,
    by::BY,
    ibyte,
    smallsize,
    smallsort::SORT,
    xs_mutable::Union{Val{true},Val{false}} = Val(false),
    bufs = [zeros(Int, 256)],
    ibuf = 1,
) where {BY,SORT}

    offsets, allpadded = _countmsd!(fill!(bufs[ibuf], 0), xs, by, ibyte)
    @DBG counts = copy(offsets)
    allpadded && return ys

    # exclusive scan
    (acc, offsets[1]) = (offsets[1], 0)
    for i in 2:length(offsets)
        @inbounds (acc, offsets[i]) = (acc + offsets[i], acc)
    end

    _scatter!(ys, xs, eachindex(xs), offsets, by, ibyte)
    @DBG @check offsets == cumsum(counts)

    # Skip next round if know we don't need it by the type information:
    will_be_allpadded(eltype(xs), by, @stat(ibyte + 1)) && return ys

    if ibuf + 1 > length(bufs)
        push!(bufs, zeros(Int, 256))
    end
    _foreach_remaining_subrange(offsets) do idx
        ys_chunk = view(ys, idx)
        if length(idx) <= smallsize
            sort!(ys_chunk, smallsort, Base.ord(isless, by, nothing, Base.Forward))
        else
            if xs_mutable === Val(true)
                xs_chunk = copyto!(view(xs, idx), ys_chunk)
            else
                xs_chunk = copy(ys_chunk)
            end
            _stablemsd_seq!(
                ys_chunk,
                xs_chunk,
                by,
                static_if_leq(@stat(ibyte + 1), static(8)),
                smallsize,
                smallsort,
                Val(true),
                bufs,
                ibuf + 1,
            )
        end
        return
    end

    return ys
end

@inline function will_be_allpadded(::Type{ElType}, by::F, ibyte) where {ElType,F}
    if ibyte == 1
        return will_be_allpadded(ElType, by, static(1))
    elseif ibyte == 2
        return will_be_allpadded(ElType, by, static(2))
    elseif ibyte == 3
        return will_be_allpadded(ElType, by, static(3))
    elseif ibyte == 4
        return will_be_allpadded(ElType, by, static(4))
    elseif ibyte == 5
        return will_be_allpadded(ElType, by, static(5))
    elseif ibyte == 6
        return will_be_allpadded(ElType, by, static(6))
    elseif ibyte == 7
        return will_be_allpadded(ElType, by, static(7))
    elseif ibyte == 8
        return will_be_allpadded(ElType, by, static(8))
    elseif ibyte == 9
        return will_be_allpadded(ElType, by, static(9))
    end
    return false
end
@inline function will_be_allpadded(
    ::Type{ElType},
    by::F,
    ibyte::StaticInteger,
) where {ElType,F}
    T = Core.Compiler.return_type(Tuple{ElType}) do x
        msd_at(by(x), ibyte)
    end
    return T === Nothing
end

zerocounts() = zeros(MVector{256,Int})
const COUNTS_POOL = Pool(zerocounts)
module _CountsPoolInitializer
function __init__ end
end
_CountsPoolInitializer.__init__() = initpool!(COUNTS_POOL)

"""
    _countmsd!(counts, xs, by, ibyte) -> (counts, allpadded::Bool)
"""
function _countmsd!(counts, xs, by, ibyte)
    allpadded = true
    i = firstindex(xs)
    n = lastindex(xs) - 8
    @inbounds while i <= n
        tryprefetch(xs[i+4], ibyte)
        tryprefetch(xs[i+5], ibyte)
        tryprefetch(xs[i+6], ibyte)
        tryprefetch(xs[i+7], ibyte)
        x1 = xs[i]
        x2 = xs[i+1]
        x3 = xs[i+2]
        x4 = xs[i+3]
        k1 = msd_at(by(x1), ibyte)
        k2 = msd_at(by(x2), ibyte)
        k3 = msd_at(by(x3), ibyte)
        k4 = msd_at(by(x4), ibyte)
        if k1 !== nothing
            allpadded = false
        end
        if k2 !== nothing
            allpadded = false
        end
        if k3 !== nothing
            allpadded = false
        end
        if k4 !== nothing
            allpadded = false
        end
        counts[something(k1, 1)] += 1
        counts[something(k2, 1)] += 1
        counts[something(k3, 1)] += 1
        counts[something(k4, 1)] += 1
        i += 4
    end
    while i <= lastindex(xs)
        k = msd_at(by(@inbounds xs[i]), ibyte)
        if k !== nothing
            allpadded = false
        end
        @inbounds counts[something(k, 1)] += 1
        i += 1
    end
    return counts, allpadded
end

function _scatter!(ys, xs, idx, exclusive_offsets, by, ibyte)
    offsets = exclusive_offsets  # exclusive scan of `counts`
    i = first(idx)
    n = last(idx) - 8
    @inbounds while i <= n
        tryprefetch(xs[i+4], ibyte)
        tryprefetch(xs[i+5], ibyte)
        tryprefetch(xs[i+6], ibyte)
        tryprefetch(xs[i+7], ibyte)
        x1 = xs[i]
        x2 = xs[i+1]
        x3 = xs[i+2]
        x4 = xs[i+3]
        k1 = something(msd_at(by(x1), ibyte), 1)
        k2 = something(msd_at(by(x2), ibyte), 1)
        k3 = something(msd_at(by(x3), ibyte), 1)
        k4 = something(msd_at(by(x4), ibyte), 1)
        j1 = offsets[k1] += 1
        j2 = offsets[k2] += 1
        j3 = offsets[k3] += 1
        j4 = offsets[k4] += 1
        ys[j1] = x1
        ys[j2] = x2
        ys[j3] = x3
        ys[j4] = x4
        i += 4
    end
    while i <= last(idx)
        j = offsets[something(msd_at(by(xs[i]), ibyte), 1)] += 1
        ys[j] = xs[i]
        i += 1
    end
    # `offsets` here is now inclusive scan of `counts`
end

function _foreach_remaining_subrange(
    f,
    inclusive_offsets,
    indices::UnitRange = firstindex(inclusive_offsets):lastindex(inclusive_offsets),
)
    prev = prev = indices[1] == 1 ? 0 : inclusive_offsets[indices[1]-1]
    for x in indices
        i = @inbounds inclusive_offsets[x]
        if !(prev == i || prev + 1 == i)
            f(prev+1:i)
        end
        prev = i
    end
end

function _spawn_foreach_remaining_subrange(
    f,
    inclusive_offsets,
    basesize::Integer,
    indices::UnitRange = firstindex(inclusive_offsets):lastindex(inclusive_offsets),
)
    @inline function maybe_nonempty_subrange(subindices::UnitRange)
        isempty(subindices) && return nothing
        local prev = subindices[1] == 1 ? 0 : inclusive_offsets[subindices[1]-1]
        local curr = inclusive_offsets[subindices[end]]
        if prev == curr || prev + 1 == curr
            return nothing
        else
            return prev+1:curr
        end
    end

    let subrange = maybe_nonempty_subrange(indices)
        subrange === nothing && return
        if length(subrange) <= basesize || length(indices) < 2
            _foreach_remaining_subrange(f, inclusive_offsets, indices)
            return
        end
    end

    m = (last(indices) - first(indices) + 1) ÷ 2 + first(indices)
    left = first(indices):m-1
    right = m:last(indices)
    leftrange = maybe_nonempty_subrange(left)
    rightrange = maybe_nonempty_subrange(right)
    if leftrange === nothing
        if rightrange === nothing
        else
            _spawn_foreach_remaining_subrange(f, inclusive_offsets, basesize, right)
        end
    else
        if rightrange === nothing
            _spawn_foreach_remaining_subrange(f, inclusive_offsets, basesize, left)
        else
            @sync begin
                @spawn _spawn_foreach_remaining_subrange(
                    f,
                    inclusive_offsets,
                    basesize,
                    right,
                )
                _spawn_foreach_remaining_subrange(f, inclusive_offsets, basesize, left)
            end
        end
    end
    return
end

@inline msd_at(x, ibyte::Integer) = sizeof(x) < ibyte ? nothing : _msd_at(x, ibyte) + 1

@inline _msd_at(x::Unsigned, ibyte::Integer) =
    if ibyte == sizeof(x)
        Int(x % UInt8)
    else
        Int((x >> (8 * (sizeof(x) - ibyte))) % UInt8)
    end

# Taken from SortingAlgorithms.jl:
@inline _msd_at(x::Signed, ibyte::Integer) = _msd_at(unsigned(xor(x, typemin(x))), ibyte)
@inline function _msd_at(x::Float32, ibyte::Integer)
    y = reinterpret(Int32, x)
    return _msd_at(reinterpret(UInt32, ifelse(y < 0, ~y, xor(y, typemin(Int32)))), ibyte)
end
@inline function _msd_at(x::Float64, ibyte::Integer)
    y = reinterpret(Int64, x)
    return _msd_at(reinterpret(UInt64, ifelse(y < 0, ~y, xor(y, typemin(Int64)))), ibyte)
end
# TODO: normalize NaN?

@inline _msd_at(x::AbstractChar, ibyte::Integer) = _msd_at(UInt32(x), ibyte)

@inline function msd_at(x::AbstractString, ibyte::Integer)
    n, r = divrem(ibyte - 1, sizeof(eltype(x)))
    if ibyte isa StaticInteger
        c = _get(x, static(n + 1))
    else
        c = _get(x, n + 1)
    end
    c === nothing && return nothing
    return msd_at(c, r + 1)
end

@inline function _get(x::AbstractString, i::Integer)
    @DBG @check i > 0
    y = iterate(x)
    while true
        y === nothing && return nothing
        i -= 1
        i == 0 && return y[1]
        y = iterate(x, y[2])
    end
end

@inline function _get(x::AbstractString, i::StaticInteger)
    @DBG @check i > 1
    y0 = @inbounds iterate(x)
    y0 === nothing && return nothing
    i === static(1) && return y0[1]
    return _afoldl(y0[2], ntuple(identity, static(i - 1))...) do state, j
        j::StaticInteger
        y = @inbounds iterate(x, state)
        y === nothing && return Reduced(nothing)
        j === static(i - 1) && return Reduced(y[1])
        y[2]
    end
end

@inline function msd_at(x::AbstractArray, ibyte::Integer)
    n, r = divrem(ibyte - 1, sizeof(eltype(x)))
    checkbounds(Bool, x, n + 1) || return nothing
    return msd_at((@inbounds x[n+1]), r + 1)
end

msdsort!(
    xs;
    basesize = default_basesize(xs),
    by = identity,
    smallsize = DEFAULT_SMALLSIZE,
    smallsort = Base.Sort.DEFAULT_UNSTABLE,
    _copy = copy(xs),
) = _stablemsd!(xs, _copy, by, static(1), basesize, smallsize, smallsort, Val(true))

msdsort(
    xs;
    basesize = default_basesize(xs),
    by = identity,
    smallsize = DEFAULT_SMALLSIZE,
    smallsort = Base.Sort.DEFAULT_UNSTABLE,
) = _stablemsd!(copy(xs), xs, by, static(1), basesize, smallsize, smallsort)

msdsort_seq!(
    xs;
    by = identity,
    smallsize = DEFAULT_SMALLSIZE,
    smallsort = Base.Sort.DEFAULT_UNSTABLE,
    _copy = copy(xs),
) = _stablemsd_seq!(xs, _copy, by, static(1), smallsize, smallsort, Val(true))

msdsort_seq(
    xs;
    by = identity,
    smallsize = DEFAULT_SMALLSIZE,
    smallsort = Base.Sort.DEFAULT_UNSTABLE,
) = _stablemsd_seq!(copy(xs), xs, by, static(1), smallsize, smallsort)


text:
# d for integers:
@printf("%d\n", 1e5) #> 100000
x = 7.35679
# f = float format, rounded if needed:
@printf("x = %0.3f\n", x) #> 7.357
aa = 1.5231071779744345
bb = 33.976886930000695
@printf("%.2f %.2f\n", aa, bb) #> 1.52 33.98
# or to create another string:
str = @sprintf("%0.3f", x)
show(str) #> "7.357"
println()
# e = scientific format with e:
@printf("%0.6e\n", x) #> 7.356790e+00
# c = for characters:
@printf("output: %c\n", 'α') #> output: α
# s for strings:
@printf("%s\n", "I like Julia")
# right justify:
@printf("%50s\n", "text right justified!")
# p for pointers:
@printf("a pointer: %p\n", 1e10) #> a pointer: 0x00000002540be400

max_stars_repo_path: Module 1/Chapter02/formatting.jl

max_issues_repo_path: Module 1/Chapter02/formatting.jl

max_forks_repo_path: Module 1/Chapter02/formatting.jl


text:
@pyimport openfermion
@pyimport openfermionpyscf

function get_molecule_openfermion()
    basis = "sto-3g"
    multiplicity = 1
    run_scf = 1
    run_mp2 = 1
    run_cisd = 0
    run_ccsd = 0
    run_fci = 1
    bond_length = 0.7414
    geometry = [("H", (0., 0., 0.)), ("H", (0., 0., bond_length))]
    molecule = openfermion.chem.MolecularData(
        geometry, basis, multiplicity)

    molecule = openfermionpyscf.run_pyscf(molecule)
                         # run_scf=run_scf,
                         # run_mp2=run_mp2,
                         # run_cisd=run_cisd,
                         # run_ccsd=run_ccsd,
                         # run_fci=run_fci)
    return molecule
end

mol_of = get_molecule_openfermion();


text:
using Documenter
using PlasmaProperties

push!(LOAD_PATH,"../src/")
makedocs(sitename="PlasmaProperties.jl",
         pages = [
            "Index" => "index.md",
         ],
         format = Documenter.HTML(prettyurls = false)
)
deploydocs(
    repo = "github.com/ep2lab/PlasmaProperties.jl.git",
    devbranch = "main"
)


text:
# Write header

for q in 0:100000
    for a in 0:100000
        x = q * 0.001
        y = a * 0.001
        row = "$(x)         ,$(y)  "  # DO
        print("Progress: $(row) \r")
        flush(stdout)

    end
end


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

using ValueShapes
using Test

@testset "const_value_shape" begin
    @inferred(size(ConstValueShape(42))) == ()
    @inferred(eltype(ConstValueShape(42))) == Int
    @inferred(totalndof(ConstValueShape(42))) == 0

    @inferred(size(ConstValueShape(rand(2,3)))) == (2,3)
    @inferred(ValueShapes.default_unshaped_eltype(ConstValueShape(rand(Float32,2,3)))) == Float32
    @inferred(totalndof(ConstValueShape(rand(2,3)))) == 0

    @test @inferred(ConstValueShape([1 4; 3 2])(undef)) == [1 4; 3 2]
    @test @inferred(ConstValueShape([1 4; 3 2])(Int[])) == [1 4; 3 2]

    shape = ConstValueShape([1 4; 3 2])

    @test typeof(@inferred(Vector{Int32}(undef, shape))) == Vector{Int32}
    @test size(@inferred(Vector{Int32}(undef, shape))) == (0,)
end


text:
using Random, SnpArrays
using LinearAlgebra,  BenchmarkTools
Random.seed!(1234)
import TraitSimulation: snparray_simulation
n = 10
p = 2
d = 2
m = 2

function generateSPDmatrix(n)
	A = rand(n)
	m = 0.5 * (A * A')
	PDmat = m + (n * Diagonal(ones(n)))
end

function generateRandomVCM(n::Int64, p::Int64, d::Int64, m::Int64)
	# n-by-p design matrix
	X = randn(n, p)

	# p-by-d mean component regression coefficient for each trait
	B = hcat(ones(p, 1), rand(p))

	V = ntuple(x -> zeros(n, n), m)
	for i = 1:m-1
	  copy!(V[i], generateSPDmatrix(n))
	end
	copy!(V[end], Diagonal(ones(n))) # last covarianec matrix is identity

	# a tuple of m d-by-d variance component parameters
	Σ = ntuple(x -> zeros(d, d), m)
	for i in 1:m
	  copy!(Σ[i], generateSPDmatrix(d))
	end
	return(X, B, Σ, V)
end

X, B, Σ, V = generateRandomVCM(n, p, d, m)
# specify the variance compoents of the model
varcomp = @vc Σ[1] ⊗ V[1] + Σ[2] ⊗ V[2]
test_vcm1 = VCMTrait(X, B, varcomp)

model_info_vcm = @capture_out show(test_vcm1)
vc_info = @capture_out show(varcomp)

@test model_info_vcm != nothing
@test vc_info != nothing

@test nsamplesize(test_vcm1) == n
@test nvc(test_vcm1) == m
@test ntraits(test_vcm1) == d
@test neffects(test_vcm1) == p

y_vcm = simulate(test_vcm1)

@test size(y_vcm) == size(X*B)
number_independent_simulations = 5
@test length(simulate(test_vcm1, number_independent_simulations)) == number_independent_simulations

@test test_vcm1 isa TraitSimulation.AbstractTraitModel

#testing for types
minor_alle_frequency  = 0.2
nsnps = p
@test snparray_simulation([minor_alle_frequency], nsnps) isa SnpArrays.SnpArray

variance_formula2  = @vc [minor_alle_frequency][:,:] ⊗ V[1] + [minor_alle_frequency][:,:] ⊗ V[1]
trait2 = VCMTrait(["x2 + 3x1", "2 +2x1"], DataFrame(X, :auto), variance_formula2)
sigma, v = vcobjtuple(variance_formula2)
trait2_equivalent = VCMTrait(["x2 + 3x1", "2 +2x1"], DataFrame(X, :auto), [sigma...], [v...])

@test trait2_equivalent.vc[1].V == trait2.vc[1].V

# simulate the SnpArray
G = snparray_simulation([0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5], n)
γ = rand(7, 2)

vcmOBJ =  VCMTrait(X, B, G, γ, varcomp)
vcmOBJ_equivalent =  VCMTrait(X, B, G, γ, [Σ...], [V...])

@test vcmOBJ.vc[1].V == V[1]
@test vcmOBJ_equivalent.vc[1].V == V[1]
# check data was copied correctly
@test vcmOBJ.X == X

vcmOBJ2 = VCMTrait(X, B, varcomp)

@test eltype(simulate(vcmOBJ2)) == Float64

@test length(simulate(vcmOBJ2, number_independent_simulations))  == number_independent_simulations

vcmOBJ2_equivalent = VCMTrait(X, B, [Σ...], [V...])
@test vcmOBJ2.vc[1].V == vcmOBJ2_equivalent.vc[1].V
@test vcmOBJ2.G != nothing
@test vcmOBJ2.μ == X*B

# check if the preallocated memory that will loop over for the matrix normal  (to sum the variance components) before we sum the mean matrix
@test size(vcmOBJ2.Z) == size(vcmOBJ2.μ)

y_alternative = zeros(size(vcmOBJ2.μ))
function testing_allocations(vcmOBJ2, y_alternative)
    TraitSimulation.simulate!(y_alternative, vcmOBJ2)
end
bmrk  = @benchmark testing_allocations(vcmOBJ2, y_alternative)

@test bmrk.memory == 0
@test bmrk.allocs ==  0


text:
module STL

using CxxInterface
using STL_jll
using TestAbstractTypes

# Run all self-tests when the package is loaded
const selftests = Any[]
function __init__()
    # If the library name is just a placeholder
    if libSTL ≠ "libSTL"
        for test in selftests
            test()
        end
    end
    return empty!(selftests)
end

eval(cxxsetup())

include("Defs.jl")
include("Types.jl")

include("StdMap.jl")
include("StdSharedPtr.jl")
include("StdString.jl")
include("StdVector.jl")

end