text:
# Generic function wrappers

# Functions taking one argument
@inline function llvmwrap(::Val{Op}, v1::_Vec{_W,T1}) where {Op,_W,T1}
    llvmwrap(Val{Op}(), v1, T1)
end
@generated function llvmwrap(::Val{Op}, v1::_Vec{_W,T1}, ::Type{R}) where {Op,_W,T1,R}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    typr = llvmtype(R)
    vtypr = "<$W x $typr>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    flags = [""]
    if Op in FASTOPS# && T1 <: FloatingTypes
        push!(flags, fastflags(promote_type(T1,R)))
    end
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1)")
        push!(instrs, "%res = call" * join(flags, " ") * " $vtypr $ins($vtyp1 %0)")
    else
        if Op === :~
            @assert T1 <: IntegerTypes
            otherval = -1
        elseif Op === :inv
            @assert T1 <: FloatingTypes
            otherval = 1.0
        else
            otherval = 0
        end
        otherarg = llvmconst(W, T1, otherval)
        push!(instrs, "%res = $ins" * join(flags, " ") * " $vtyp1 $otherarg, %0")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$R}, Tuple{Vec{$W,$T1}}, v1
        )
    end
end

# Functions taking one Bool argument
@generated function llvmwrap(::Val{Op}, v1::_Vec{_W,Bool}, ::Type{Bool}) where {Op,_W}
    W = _W + 1
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$W x $btyp>"
    ins = llvmins(Op, W, Bool)
    decls = String[]
    instrs = String[]
    push!(instrs, "%arg1 = trunc $vbtyp %0 to <$W x i1>")
    otherarg = llvmconst(W, Bool, true)
    push!(instrs, "%res = $ins <$W x i1> $otherarg, %arg1")
    push!(instrs, "%resb = zext <$W x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,Bool}, Tuple{Vec{$W,Bool}}, v1)
    end
end

# Functions taking two arguments
@generated function llvmwrap(
    ::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T2}, ::Type{R} = T1
) where {Op,_W,T1,T2,R}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$W x $typ2>"
    typr = llvmtype(R)
    vtypr = "<$W x $typr>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    flags = [""]
    T = promote_type(T1,T2,R)
    if Op in FASTOPS# && T <: FloatingTypes
        push!(flags, fastflags(T))
    end
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2)")
        push!(instrs, "%res = call" * join(flags, " ") * " $vtypr $ins($vtyp1 %0, $vtyp2 %1)")
    else
        push!(instrs, "%res = $ins" * join(flags, " ") * " $vtyp1 %0, %1")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$R}, Tuple{Vec{$W,$T1}, Vec{$W,$T2}},
            v1, v2)
    end
end
@generated function llvmwrap_notfast(
    ::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T2}, ::Type{R} = T1
) where {Op,_W,T1,T2,R}
    @assert isa(Op, Symbol)
    W = _W + 1
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$W x $typ2>"
    typr = llvmtype(R)
    vtypr = "<$W x $typr>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2)")
        push!(instrs, "%res = call $vtypr $ins($vtyp1 %0, $vtyp2 %1)")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$R}, Tuple{Vec{$W,$T1}, Vec{$W,$T2}},
            v1, v2)
    end
end

# Functions taking two arguments, returning Bool
@generated function llvmwrap(
    ::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T2}, ::Type{Bool}
) where {Op,_W,T1,T2}
    @assert isa(Op, Symbol)
    W = _W + 1
    btyp = llvmtype(Bool)
    vbtyp = "<$W x $btyp>"
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$W x $typ2>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    # if false && W == 1
        # append!(instrs, array2vector("%arg1", W, typ1, "%0", "%arg1arr"))
        # append!(instrs, array2vector("%arg2", W, typ2, "%1", "%arg2arr"))
        # push!(instrs, "%cond = $ins $vtyp1 %arg1, %arg2")
        # push!(instrs, "%res = zext <$W x i1> %cond to $vbtyp")
        # append!(instrs, vector2array("%resarr", W, btyp, "%res"))
        # push!(instrs, "ret $abtyp %resarr")
    # else
    push!(instrs, "%res = $ins $vtyp1 %0, %1")
    push!(instrs, "%resb = zext <$W x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    # end
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,Bool}, Tuple{Vec{$W,$T1}, Vec{$W,$T2}},
            v1, v2)
    end
end
# Functions taking two arguments, returning a bitmask
@generated function llvmwrap_bitmask(::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T1}) where {Op,_W,T1}
    @assert isa(Op, Symbol)
    W = _W + 1
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    # if T1 <: FloatingTypes && T2 <: FloatingTypes
        # push!(instrs, "%res = $ins reassoc $vtyp1 %0, %1")
    # else
    push!(instrs, "%res = $ins $vtyp1 %0, %1")
    maskbits = max(8, VectorizationBase.nextpow2(W))
    bitcastname = maskbits == W ? "resu" : "resutrunc"
    push!(instrs, "%$(bitcastname) = bitcast <$W x i1> %res to i$W")
    if maskbits != W
        push!(instrs, "%resu = zext i$W %$(bitcastname) to i$maskbits")
    end
    push!(instrs, "ret i$maskbits %resu")
    julia_mask_type = VectorizationBase.mask_type(maskbits)
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            $julia_mask_type, Tuple{Vec{$W,$T1}, Vec{$W,$T1}},
            v1, v2)
    end
end
@generated function llvmwrap_bitmask_notfast(::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T1}) where {Op,_W,T1}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    ins = if Op == :(==)
        "fcmp oeq"
    elseif Op == :(!=)
        "fcmp une"
    elseif Op == :(>)
        "fcmp ogt"
    elseif Op == :(>=)
        "fcmp oge"
    elseif Op == :(<)
        "fcmp olt"
    elseif Op == :(<=)
        "fcmp ole"
    else
        throw("Op $Op not recognized.")
    end
    decls = String[]
    instrs = String[]
    # if T1 <: FloatingTypes && T2 <: FloatingTypes
        # push!(instrs, "%res = $ins reassoc $vtyp1 %0, %1")
    # else
    push!(instrs, "%res = $ins $vtyp1 %0, %1")
    # end
    # push!(instrs, "%resb = zext <$W x i1> %res to $vbtyp")
    maskbits = max(8, VectorizationBase.nextpow2(W))
    bitcastname = maskbits == W ? "resu" : "resutrunc"
    push!(instrs, "%$(bitcastname) = bitcast <$W x i1> %res to i$W")
    if maskbits != W
        push!(instrs, "%resu = zext i$W %$(bitcastname) to i$maskbits")
    end
    push!(instrs, "ret i$maskbits %resu")
    julia_mask_type = VectorizationBase.mask_type(maskbits)
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            $julia_mask_type, Tuple{Vec{$W,$T1}, Vec{$W,$T1}},
            v1, v2)
    end
end

# Functions taking two Bool arguments, returning Bool
@generated function llvmwrap(::Val{Op}, v1::_Vec{_W,Bool}, v2::_Vec{_W,Bool}, ::Type{Bool} = Bool) where {Op,_W}
    W = _W + 1
    @assert isa(Op, Symbol)
    btyp = llvmtype(Bool)
    vbtyp = "<$W x $btyp>"
    ins = llvmins(Op, W, Bool)
    decls = String[]
    instrs = String[]
    push!(instrs, "%arg1 = trunc $vbtyp %0 to <$W x i1>")
    push!(instrs, "%arg2 = trunc $vbtyp %1 to <$W x i1>")
    push!(instrs, "%res = $ins <$W x i1> %arg1, %arg2")
    push!(instrs, "%resb = zext <$W x i1> %res to $vbtyp")
    push!(instrs, "ret $vbtyp %resb")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,Bool}, Tuple{Vec{$W,Bool}, Vec{$W,Bool}},
            v1, v2)
    end
end

# Functions taking three arguments
@generated function llvmwrap(
    ::Val{Op}, s1::T, s2::T, s3::T
) where {Op, T <: Number}
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    ins = llvmins(Op, -1, T)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $typ $ins($typ, $typ, $typ)")
        push!(instrs, "%res = call $typ $ins($typ %0, $typ %1, $typ %2)")
    else
        push!(instrs, "%res = $ins $typ %0, %1, %2")
    end
    push!(instrs, "ret $typ %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            $T, Tuple{$T,$T,$T}, s1, s2, s3
        )
    end
end
@generated function llvmwrap(
    ::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T2}, v3::_Vec{_W,T3}, ::Type{R}
) where {Op,_W,T1,T2,T3,R}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$W x $typ2>"
    typ3 = llvmtype(T3)
    vtyp3 = "<$W x $typ3>"
    typr = llvmtype(R)
    vtypr = "<$W x $typr>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2, $vtyp3)")
        push!(instrs,
            "%res = call $vtypr $ins($vtyp1 %0, $vtyp2 %1, $vtyp3 %2)")
    else
        push!(instrs, "%res = $ins $vtyp1 %0, %1, %2")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$R},
            Tuple{Vec{$W,$T1}, Vec{$W,$T2}, Vec{$W,$T3}},
            v1, v2, v3)
    end
end

@inline llvmwrap(::Val{Op}, v1::_Vec{_W,R}, v2::_Vec{_W}, v3::_Vec{_W}) where {Op,_W,R} = llvmwrap(Val{Op}(), v1, v2, v3, R)
@inline llvmwrap_fast(::Val{Op}, v1::_Vec{_W,R}, v2::_Vec{_W}, v3::_Vec{_W}) where {Op,_W,R} = llvmwrap_fast(Val{Op}(), v1, v2, v3, R)
@generated function llvmwrap_fast(
    ::Val{Op}, v1::_Vec{_W,T1}, v2::_Vec{_W,T2}, v3::_Vec{_W,T3}, ::Type{R}
) where {Op,_W,T1,T2,T3,R}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ1 = llvmtype(T1)
    vtyp1 = "<$W x $typ1>"
    typ2 = llvmtype(T2)
    vtyp2 = "<$W x $typ2>"
    typ3 = llvmtype(T3)
    vtyp3 = "<$W x $typ3>"
    typr = llvmtype(R)
    vtypr = "<$W x $typr>"
    ins = llvmins(Op, W, T1)
    decls = String[]
    instrs = String[]
    if ins[1] == '@'
        push!(decls, "declare $vtypr $ins($vtyp1, $vtyp2, $vtyp3)")
        push!(instrs,
            "%res = call ninf nsz arcp contract afn reassoc $vtypr $ins($vtyp1 %0, $vtyp2 %1, $vtyp3 %2)")
    else
        push!(instrs, "%res = $ins ninf nsz arcp contract afn reassoc $vtyp1 %0, %1, %2")
    end
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$R},
            Tuple{Vec{$W,$T1}, Vec{$W,$T2}, Vec{$W,$T3}},
            v1, v2, v3)
    end
end

@generated function llvmwrapshift(
    ::Val{Op}, v1::_Vec{_W,T}, ::Val{I}
) where {Op,_W,T,I}
    W = _W + 1
    @assert isa(Op, Symbol)
    if I >= 0
        op = Op
        i = I
    else
        if Op === :>> || Op === :>>>
            op = :<<
        else
            @assert Op === :<<
            if T <: Unsigned
                op = :>>>
            else
                op = :>>
            end
        end
        i = -I
    end
    @assert op in (:<<, :>>, :>>>)
    @assert i >= 0
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    ins = llvmins(op, W, T)
    decls = String[]
    instrs = String[]
    nbits = 8*sizeof(T)
    if (op === :>> && T <: IntTypes) || i < nbits
        count = llvmconst(W, T, min(nbits-1, i))
        push!(instrs, "%res = $ins $vtyp %0, $count")
        push!(instrs, "ret $vtyp %res")
    else
        push!(instrs, "return $vtyp zeroinitializer")
    end
    quote
        $(Expr(:meta, :inline))
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$T}, Tuple{Vec{$W,$T}}, v1)
    end
end

@generated function llvmwrapshift(
    ::Val{Op}, v1::_Vec{_W,T}, x2::Unsigned
) where {Op,_W,T}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    ins = llvmins(Op, W, T)
    decls = String[]
    instrs = String[]
    append!(instrs, scalar2vector("%count", W, typ, "%1"))
    nbits = 8*sizeof(T)
    push!(instrs, "%res = $ins $vtyp %0, %count")
    # push!(instrs, "%inbounds = icmp ult $typ %1, $nbits")
    # if Op === :>> && T <: IntTypes
    #     nbits1 = llvmconst(W, T, 8*sizeof(T)-1)
    #     push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
    #     push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp %limit")
    # else
    #     push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
    # end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        # Wote that this function might be called with out-of-bounds
        # values for x2, assuming that the results are then ignored
        llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$T}, Tuple{Vec{$W,$T}, $T}, v1, x2 % $T)
    end
end

@generated function llvmwrapshift(
    ::Val{Op}, v1::_Vec{_W,T}, x2::Integer
) where {Op,_W,T}
    W = _W + 1
    if Op === :>> || Op === :>>>
        NegOp = :<<
    else
        @assert Op === :<<
        if T <: Unsigned
            NegOp = :>>>
        else
            NegOp = :>>
        end
    end
    ValOp = Val{Op}()
    ValNegOp = Val{NegOp}()
    quote
        $(Expr(:meta, :inline))
        x2 < 0 ? llvmwrapshift($ValNegOp, v1, unsigned(-x2)) : llvmwrapshift($ValOp, v1, unsigned(x2))
    end
end

@generated function llvmwrapshift(
    ::Val{Op},
    v1::_Vec{_W,T},
    v2::_Vec{_W,U}
) where {Op,_W,T,U<:UIntTypes}
    W = _W + 1
    @assert isa(Op, Symbol)
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    ins = llvmins(Op, W, T)
    decls = String[]
    instrs = String[]
    push!(instrs, "%res = $ins $vtyp %0, %1")
    # push!(instrs, "%tmp = $ins $vtyp %0, %1")
    # nbits = llvmconst(W, T, 8*sizeof(T))
    # push!(instrs, "%inbounds = icmp ult $vtyp %1, $nbits")
    # if Op === :>> && T <: IntTypes
    #     nbits1 = llvmconst(W, T, 8*sizeof(T)-1)
    #     push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
    #     push!(instrs,
    #         "%res = select <$W x i1> %inbounds, $vtyp %tmp, $vtyp %limit")
    # else
    #     push!(instrs,
    #         "%res = select <$W x i1> %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
    # end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        llvmcall(
            $((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$T}, Tuple{Vec{$W,$T}, Vec{$W,$T}},
            v1, extract_data(SVec(v2) % $T)
        )
    end
end

@generated function llvmwrapshift(
    ::Val{Op},
    v1::_Vec{_W,T},
    v2::_Vec{_W,U}
) where {Op,_W,T,U<:IntegerTypes}
    W = _W + 1
    if Op === :>> || Op === :>>>
        NegOp = :<<
    else
        @assert Op === :<<
        if T <: Unsigned
            NegOp = :>>>
        else
            NegOp = :>>
        end
    end
    ValOp = Val{Op}()
    ValNegOp = Val{NegOp}()
    quote
        $(Expr(:meta, :inline))
        vifelse(
            vgreater_or_equal(v2, 0),
            llvmwrapshift($ValOp, v1, extract_data(SVec(v2) % unsigned(U))),
            llvmwrapshift($ValNegOp, v1, extract_data(SVec(vsub(v2)) % unsigned(U)))
        )
    end
end

# @generated function llvmwrapshift(
#     ::Val{Op}, v1::_Vec{_W,T}, ::Val{I}
# ) where {Op,W,T,I}
#     @assert isa(Op, Symbol)
#     if I >= 0
#         op = Op
#         i = I
#     else
#         if Op === :>> || Op === :>>>
#             op = :<<
#         else
#             @assert Op === :<<
#             if T <: Unsigned
#                 op = :>>>
#             else
#                 op = :>>
#             end
#         end
#         i = -I
#     end
#     @assert op in (:<<, :>>, :>>>)
#     @assert i >= 0
#     typ = llvmtype(T)
#     vtyp = "<$W x $typ>"
#     ins = llvmins(op, W, T)
#     decls = String[]
#     instrs = String[]
#     nbits = 8*sizeof(T)
#     if (op === :>> && T <: IntTypes) || i < nbits
#         count = llvmconst(W, T, min(nbits-1, i))
#         push!(instrs, "%res = $ins $vtyp %0, $count")
#         push!(instrs, "ret $vtyp %res")
#     else
#         push!(instrs, "return $vtyp zeroinitializer")
#     end
#     quote
#         $(Expr(:meta, :inline))
#         llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{W,T}, Tuple{Vec{W,T}}, v1)
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op}, v1::_Vec{_W,T}, x2::Unsigned
# ) where {Op,W,T}
#     @assert isa(Op, Symbol)
#     typ = llvmtype(T)
#     vtyp = "<$W x $typ>"
#     ins = llvmins(Op, W, T)
#     decls = String[]
#     instrs = String[]
#     append!(instrs, scalar2vector("%count", W, typ, "%1"))
#     nbits = 8*sizeof(T)
#     push!(instrs, "%tmp = $ins $vtyp %0, %count")
#     push!(instrs, "%inbounds = icmp ult $typ %1, $nbits")
#     if Op === :>> && T <: IntTypes
#         nbits1 = llvmconst(W, T, 8*sizeof(T)-1)
#         push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
#         push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp %limit")
#     else
#         push!(instrs, "%res = select i1 %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
#     end
#     push!(instrs, "ret $vtyp %res")
#     quote
#         $(Expr(:meta, :inline))
#         # Note that this function might be called with out-of-bounds
#         # values for x2, assuming that the results are then ignored
#         llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{W,T}, Tuple{Vec{W,T}, T}, v1, x2 % T)
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op}, v1::_Vec{_W,T}, x2::Integer
# ) where {Op,W,T}
#     if Op === :>> || Op === :>>>
#         NegOp = :<<
#     else
#         @assert Op === :<<
#         if T <: Unsigned
#             NegOp = :>>>
#         else
#             NegOp = :>>
#         end
#     end
#     ValOp = Val{Op}()
#     ValNegOp = Val{NegOp}()
#     quote
#         $(Expr(:meta, :inline))
#         x2 < 0 ? llvmwrapshift($ValNegOp, v1, unsigned(-x2)) : llvmwrapshift($ValOp, v1, unsigned(x2))
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op},
#     v1::_Vec{_W,T},
#     v2::_Vec{_W,U}
# ) where {Op,W,T,U<:UIntTypes}
#     @assert isa(Op, Symbol)
#     typ = llvmtype(T)
#     vtyp = "<$W x $typ>"
#     ins = llvmins(Op, W, T)
#     decls = String[]
#     instrs = String[]
#     push!(instrs, "%res = $ins $vtyp %0, %1")
#     # push!(instrs, "%tmp = $ins $vtyp %0, %1")
#     # nbits = llvmconst(W, T, 8*sizeof(T))
#     # push!(instrs, "%inbounds = icmp ult $vtyp %1, $nbits")
#     # if Op === :>> && T <: IntTypes
#     #     nbits1 = llvmconst(W, T, 8*sizeof(T)-1)
#     #     push!(instrs, "%limit = $ins $vtyp %0, $nbits1")
#     #     push!(instrs,
#     #         "%res = select <$W x i1> %inbounds, $vtyp %tmp, $vtyp %limit")
#     # else
#     #     push!(instrs,
#     #         "%res = select <$W x i1> %inbounds, $vtyp %tmp, $vtyp zeroinitializer")
#     # end
#     push!(instrs, "ret $vtyp %res")
#     quote
#         $(Expr(:meta, :inline))
#         llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
#             Vec{W,T}, Tuple{Vec{W,T}, Vec{W,T}},
#             v1, vrem(v2, Vec{W,T}))
#     end
# end

# @generated function llvmwrapshift(
#     ::Val{Op},
#     v1::_Vec{_W,T},
#     v2::_Vec{_W,U}
# ) where {Op,W,T,U<:IntegerTypes}
#     if Op === :>> || Op === :>>>
#         NegOp = :<<
#     else
#         @assert Op === :<<
#         if T <: Unsigned
#             NegOp = :>>>
#         else
#             NegOp = :>>
#         end
#     end
#     ValOp = Val{Op}()
#     ValNegOp = Val{NegOp}()
#     quote
#         $(Expr(:meta, :inline))
#         vifelse(
#             vgreater_or_equal(v2, 0),
#             llvmwrapshift($ValOp, v1, vrem(v2, Vec{W,unsigned(U)})),
#             llvmwrapshift($ValNegOp, v1, vrem(vsub(v2), Vec{W,unsigned(U)}))
#         )
#     end
# end


text:
module ColorForwardDiff

using ColorTypes
using ForwardDiff
using ForwardDiff: Dual, Partials, value, partials

_partials(N, p::Dual) = partials(p)
_partials(N, k::Real) = ntuple(i->zero(k), N)

ForwardDiff.can_dual(::Type{<:Gray}) = true
ForwardDiff.can_dual(::Type{<:RGB}) = true

function dual_rgb(Tag, T, N, r, g, b)
    Dual{Tag, RGB{T}, N}(
        RGB{T}(value(r), value(g), value(b)),
        Partials(Tuple(RGB{T}.(_partials(N, r), _partials(N, g), _partials(N, b))))
    )
end

ColorTypes.RGB{T}(r::Dual{Tag,T,N}, g::Real, b::Real) where {Tag, T<:Real, N} =
    dual_rgb(Tag, T, N, r, g, b)
ColorTypes.RGB{T}(r::T, g::Dual{Tag, T, N}, b::Real) where {Tag, T<:Real, N} =
    dual_rgb(Tag, T, N, r, g, b)
ColorTypes.RGB{T}(r::T, g::T, b::Dual{Tag, T, N}) where {Tag, T<:Real, N} =
    dual_rgb(Tag, T, N, r, g, b)

ColorTypes.RGB(r::Dual{Tag,T,N}, g::Real, b::Real) where {Tag, T<:Real, N} =
    dual_rgb(Tag, T, N, r, g, b)
ColorTypes.RGB(r::T, g::Dual{Tag,T,N}, b::Real) where {Tag, T<:Real, N} =
    dual_rgb(Tag, T, N, r, g, b)
ColorTypes.RGB(r::T, g::T, b::Dual{Tag,T,N}) where {Tag, T<:Real, N} =
    dual_rgb(Tag, T, N, r, g, b)

function dual_gray(Tag, T, N, v)
    Dual{Tag, Gray{T}, N}(
        Gray{T}(value(d)),
        Partials(Tuple(Gray{T}.(partials(d))))
    )
end

ColorTypes.Gray{T}(v::Dual{Tag, T, N}) where {Tag, T, N} = dual_gray(Tag, T, N, v)

ColorTypes.Gray(v::Dual{Tag, T, N}) where {Tag, T, N} = dual_gray(Tag, T, N, v)

for f in (:red, :green, :blue)
    @eval begin
        ColorTypes.$f(c::Dual{Tag, C, N}) where {T, Tag, C<:AbstractRGB{T}, N} =
            Dual{Tag, T, N}($f(value(c)), Partials(Tuple($f.(partials(c)))))
    end
end

ColorTypes.gray(c::Dual{Tag, C, N}) where {T, Tag, C<:AbstractGray{T}, N} =
    Dual{Tag, T, N}(gray(value(c)), Partials(Tuple(gray.(partials(c)))))

end # module

ext: jl


text:
module Fonts
using FileType

TypeWoff  = FileType.Types("woff", MIME("application/font-woff"))
TypeWoff2 = FileType.Types("woff2", MIME("application/font-woff"))
TypeTtf   = FileType.Types("ttf", MIME("application/font-sfnt"))
TypeOtf   = FileType.Types("otf", MIME("application/font-sfnt"))

function Woff(buf)::Bool
	return length(buf) > 7 &&
		buf[1] == 0x77 && buf[2] == 0x4F &&
		buf[3] == 0x46 && buf[4] == 0x46 &&
		buf[5] == 0x00 && buf[6] == 0x01 &&
		buf[7] == 0x00 && buf[8] == 0x00
end

function Woff2(buf)::Bool
	return length(buf) > 7 &&
		buf[1] == 0x77 && buf[2] == 0x4F &&
		buf[3] == 0x46 && buf[4] == 0x32 &&
		buf[5] == 0x00 && buf[6] == 0x01 &&
		buf[7] == 0x00 && buf[8] == 0x00
end

function Ttf(buf)::Bool
	return length(buf) > 4 &&
		buf[1] == 0x00 && buf[2] == 0x01 &&
		buf[3] == 0x00 && buf[4] == 0x00 &&
		buf[5] == 0x00
end

function Otf(buf)::Bool
	return length(buf) > 4 &&
		buf[1] == 0x4F && buf[2] == 0x54 &&
		buf[3] == 0x54 && buf[4] == 0x4F &&
		buf[5] == 0x00
end

Font = Dict(
	TypeWoff =>  Woff,
	TypeWoff2=>  Woff2,
	TypeTtf  =>  Ttf,
	TypeOtf  =>  Otf,
)

end


text:
export BoundedEuclideanStateSpace, UnitHypercube

### Metric Typedefs/Evaluation
evaluate(M::Euclidean, v::SVector, w::SVector) = norm(w - v)
colwise{S<:AbstractVector}(d::Euclidean, v::AbstractVector, W::Vector{S}) = colwise(d, v, statevec2mat(W))
colwise{S<:AbstractVector}(d::Euclidean, W::Vector{S}, v::AbstractVector) = colwise(d, v, statevec2mat(W))

### Metric Space Instantiation
# TODO: implement show method with dynamics and cost function
BoundedEuclideanStateSpace(lo::AbstractVector, hi::AbstractVector) =
    BoundedStateSpace(SVector(lo), SVector(hi), Euclidean(), Identity())
UnitHypercube{T<:AbstractFloat}(d::Int, ::Type{T} = Float64) = BoundedEuclideanStateSpace(zeros(T,d), ones(T,d))

helper_data_structures{S}(V::Vector{S}, M::Euclidean) = (TreeDistanceDS(KDTree(V, M; reorder=false)),
                                                         EmptyControlCache())

### Steering
propagate(M::Euclidean, v::AbstractVector, u::StepControl) = v + u.t*typeof(v)(u.u)
steering_control(M::Euclidean, v::AbstractVector, w::AbstractVector) = StepControl(evaluate(M, v, w), normalize(w - v))
collision_waypoints(::Euclidean, v::AbstractVector, w::AbstractVector) = (v, w)
# TODO: specialized methods that don't check ends


text:
using DataFrames, CSV

function interpolate(x_in::AbstractArray{T1}, y_in::AbstractArray{T2}, x::T1) where {T1<:Number, T2<:Number}
    @assert x_in[1] <= x <= x_in[2]
    ((x_in[2]-x)*y_in[1]+(x-x_in[1])*y_in[2]) / (x_in[2]-x_in[1])
end

function interp_quantile_list(x::AbstractArray{T,1},q::Real;
            quantile_list=1:length(x) ) where T<:Number
   @assert issorted(quantile_list)
   @assert length(x) == length(quantile_list)
   idx_hi = searchsortedfirst(quantile_list,q)
   idx_lo = idx_hi-1
   if idx_hi>length(x)
     return x[end]
   elseif idx_lo<1
       return x[1]
   elseif q==quantile_list[idx_hi]
     return x[idx_hi]
   end
    interpolate(view(quantile_list,idx_lo:idx_hi),view(x,idx_lo:idx_hi), q)
end

mcmc_filename = "dr25_koi_mcmc_quant.csv"
koi_filename = "q1_q17_dr25_koi_original.csv"
mcmc_df = CSV.read(mcmc_filename)
tmp_koi_cat = readlines(koi_filename)
tmp_ind = 1
num_skip = 1
while tmp_koi_cat[tmp_ind][1] == '#'
    num_skip += 1
    tmp_ind += 1
end
df = CSV.read(koi_filename, allowmissing=:all, header=num_skip, rows_for_type_detect = 6000)
df = join(df, mcmc_df, on=:kepoi_name, kind=:left)

df[:koi_depth_cat] = map(x -> x, df[:koi_depth])
df[:koi_depth_err1_cat] = map(x -> x, df[:koi_depth_err1])
df[:koi_depth_err2_cat] = map(x -> x, df[:koi_depth_err2])
df[:koi_duration_cat] = map(x -> x, df[:koi_duration])
df[:koi_duration_err1_cat] = map(x -> x, df[:koi_duration_err1])
df[:koi_duration_err2_cat] = map(x -> x, df[:koi_duration_err2])

num_quantiles = 99
quantile_list = range(1/(num_quantiles+1), 1/(num_quantiles+1), num_quantiles)
q_hi = 0.5+0.34134
q_mid = 0.5
q_lo = 0.5-0.34134

for j in 1:length(df[:kepoi_name])
    if !ismissing(df[j, :depth_mean]) && df[j, :depth_mean] >= 0
        depth_quantiles = Array{Float64, 1}()
        duration_quantiles = Array{Float64, 1}()
        for i in 1:num_quantiles
            push!(depth_quantiles, df[j, Symbol("depth_q" * string(i))]*1e6)
            push!(duration_quantiles, df[j, Symbol("duration_q" * string(i))]*24)
        end
        df[j, :koi_depth] = interp_quantile_list(depth_quantiles,q_mid,quantile_list=quantile_list)
        df[j, :koi_depth_err1] =  interp_quantile_list(depth_quantiles,q_hi,quantile_list=quantile_list) - interp_quantile_list(depth_quantiles,q_mid,quantile_list=quantile_list)
        df[j, :koi_depth_err2] =  interp_quantile_list(depth_quantiles,q_lo,quantile_list=quantile_list) - interp_quantile_list(depth_quantiles,q_mid,quantile_list=quantile_list)
        df[j, :koi_duration] =  interp_quantile_list(duration_quantiles,q_mid,quantile_list=quantile_list)
        df[j, :koi_duration_err1] = interp_quantile_list(duration_quantiles,q_hi,quantile_list=quantile_list) - interp_quantile_list(duration_quantiles,q_mid,quantile_list=quantile_list)
        df[j, :koi_duration_err2] = interp_quantile_list(duration_quantiles,q_lo,quantile_list=quantile_list)- interp_quantile_list(duration_quantiles,q_mid,quantile_list=quantile_list)
    end
end

delete!(df, [:depth_mean, :depth_std, :duration_mean, :duration_std, :radr_gt1_cnt, :b_gt1_cnt])
for i in 1:num_quantiles
    delete!(df, Symbol("depth_q" * string(i)))
    delete!(df, Symbol("duration_q" * string(i)))
end

CSV.write("q1_q17_dr25_koi.csv", df)


text:
"""

ДАНО: Робот - в произвольной клетке ограниченного прямоугольной рамкой поля без внутренних перегородок и маркеров.

РЕЗУЛЬТАТ: Робот - в исходном положении в центре косого креста (в форме X) из маркеров.

(Подсказка: решение будет подобно решению задачи 1, если направление премещения Робота задавать кортежами пары значений типа HorizonSide)

"""

using HorizonSideRobots
originRobot = Robot(animate = true, "13example.sit")
#show!(originRobot)
include("../AbstractType.jl")
include("../functions.jl")

function get_line_of_cross!(robot::Robot, sides)
    path = []
    while ! ( isborder(robot,sides[1]) || isborder(robot, sides[2]) )
        for side in sides
            if (!isborder(robot,side))
                move!(robot, side)
                push!(path,Int(side))
            else
                continue
            end
        end
        putmarker!(robot)
    end
    get_back!(robot, path)
end

function master13!(robot::Robot)::Nothing

    for sides in [ [West,Sud], [West,Nord], [Ost,Sud], [Ost,Nord] ]
        get_line_of_cross!(robot, sides)
    end

    putmarker!(robot)
end

master13!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
#include("../src/ReconstructionEvaluations.jl")
using ReconstructionEvaluations
using Base.Test

include("test_count_table.jl")
include("test_overlap.jl")
include("test_geometry.jl")
include("test_syn_recovery.jl")
include("test_clustering.jl")
include("eval/runtests.jl")


text:
import Random

function generate_deck()
    cards::Vector{Card} = Card[]
    for i in 1:4
        for j in 2:14
            push!(cards, Card(Ranks.Rank(j), Suits.Suit(i)))
        end
    end
    return Deck(cards, 1)
end

function predeal_shuffle!(deck::Deck)
    Random.shuffle!(deck.cards)
end

function reshuffle!(deck::Deck)
    Random.shuffle!(deck.cards)
    deck.top = 1
end

function deal!(deck::Deck, players::Player...)
    i = 1
    while i <= length(players)
        players[i].hole.card1 = _pop_deck!(deck)
        i += 1
    end
    i = 1
    while i <= length(players)
        players[i].hole.card2 = _pop_deck!(deck)
        i += 1
    end
end

function flop!(deck::Deck, board::Board)
    _burn!(deck)
    board.flop1 = _pop_deck!(deck)
    board.flop2 = _pop_deck!(deck)
    board.flop3 = _pop_deck!(deck)
end

function turn!(deck::Deck, board::Board)
    _burn!(deck)
    board.turn = _pop_deck!(deck)
end

function river!(deck::Deck, board::Board)
    _burn!(deck)
    board.river = _pop_deck!(deck)
end

function _burn!(deck::Deck)
    _pop_deck!(deck)
end

function _pop_deck!(deck::Deck)
    card = deck.cards[deck.top]
    deck.top += 1
    return card
end

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
"""
Function creates an __Vector{String}__ of all files matching pattern. Made to imitate list.files() from R. Prints whole path.
# Additional Arguments:
- `path::String`: by default current working directory `pwd()`
- `pattern::Regex`: pattern for files to search for
- `sub::Bool`: Should subdirectories be searched to? By default `true`.
# Examples
```julia-repl
julia> lf(;pattern = Regex(".jl\\\$"))
```
"""
function lf(;path::String=pwd(), pattern::Regex=nothing, sub::Bool=true)

    files = Vector{String}()

    if sub

        for cpath = walkdir(path)

            x = joinpath.(cpath[1], cpath[3])

            if !isempty(x)
                idx = occursin.(pattern, cpath[3])
                append!(files, x[idx] )
            end

        end

    else

        files = filter(x->occursin(pattern,x), first(walkdir(path))[3])

    end

    return(files)

end


text:
"""
Translate mouse or modifier state into the corresponding XCB value.
"""
function state_xcb end

state_xcb(s::MouseState) = .&(UInt[s.left, s.middle, s.right, s.scroll_up, s.scroll_down, s.any] .* UInt[XCB_BUTTON_MASK_1, XCB_BUTTON_MASK_2, XCB_BUTTON_MASK_3, XCB_BUTTON_MASK_4, XCB_BUTTON_MASK_5, XCB_BUTTON_MASK_ANY]...)
state_xcb(s::KeyModifierState) = sum(2 .^ (0:3) .* Int[s.shift, s.ctrl, s.alt, s.super])
state_xcb(::Any) = 0

state_xcb(e::MouseEvent) = state_xcb(e.state)
state_xcb(e::KeyEvent) = state_xcb(e.modifiers)

"""
Translate an action into its corresponding XCB value.
"""
function response_type_xcb end

response_type_xcb(::Type{ButtonPressed}) = XCB_BUTTON_PRESS
response_type_xcb(::Type{ButtonReleased}) = XCB_BUTTON_RELEASE
response_type_xcb(::Type{KeyPressed}) = XCB_KEY_PRESS
response_type_xcb(::Type{KeyReleased}) = XCB_KEY_RELEASE
response_type_xcb(::Type{PointerMoves}) = XCB_MOTION_NOTIFY
response_type_xcb(::Type{PointerEntersWindow}) = XCB_ENTER_NOTIFY
response_type_xcb(::Type{PointerLeavesWindow}) = XCB_LEAVE_NOTIFY
response_type_xcb(::Type{Expose}) = XCB_EXPOSE
response_type_xcb(::Type{Resize}) = XCB_CONFIGURE_NOTIFY

event_type_xcb(::Type{<:Union{ButtonPressed, ButtonReleased}}) = xcb_button_press_event_t
event_type_xcb(::Type{<:Union{KeyPressed, KeyReleased}}) = xcb_key_press_event_t
event_type_xcb(::Type{PointerMoves}) = xcb_motion_notify_event_t
event_type_xcb(::Type{<:Union{PointerEntersWindow, PointerLeavesWindow}}) = xcb_enter_notify_event_t
event_type_xcb(::Type{ResizeEvent}) = xcb_configure_notify_event_t
event_type_xcb(::Type{ExposeEvent}) = xcb_expose_event_t

button_xcb(::ButtonLeft) = XCB_BUTTON_INDEX_1
button_xcb(::ButtonMiddle) = XCB_BUTTON_INDEX_2
button_xcb(::ButtonRight) = XCB_BUTTON_INDEX_3
button_xcb(::ButtonScrollUp) = XCB_BUTTON_INDEX_4
button_xcb(::ButtonScrollDown) = XCB_BUTTON_INDEX_5

detail_xcb(wm::XWindowManager, e::MouseEvent) = button_xcb(e.button)
detail_xcb(wm::XWindowManager, e::PointerMovesEvent) = XCB_MOTION_NORMAL
detail_xcb(wm::XWindowManager, e::PointerEntersWindowEvent) = XCB_ENTER_NOTIFY
detail_xcb(wm::XWindowManager, e::PointerLeavesWindowEvent) = XCB_LEAVE_NOTIFY
detail_xcb(wm::XWindowManager, e::EventData) = 0
detail_xcb(wm::XWindowManager, e::EventDetails) = detail_xcb(wm, e.data)
detail_xcb(wm::XWindowManager, e::EventDetails{<:KeyEvent}) = keycode(wm, e)

event_xcb(wm::XWindowManager, e::EventDetails) = event_type_xcb(action(e))(
    response_type_xcb(action(e)),
    detail_xcb(wm, e),
    0,
    e.time,
    e.win.parent_id,
    e.win.id,
    0,
    0,
    0,
    e.location...,
    state_xcb(e.data),
    true,
    false
)

event_xcb(wm::XWindowManager, e::EventDetails{<:ExposeEvent}) = event_type_xcb(action(e))(
    response_type_xcb(action(e)),
    0,
    0,
    e.win.id,
    e.location...,
    extent(e.win)...,
    0,
    0
)

event_xcb(wm::XWindowManager, e::EventDetails{<:ResizeEvent}) = event_type_xcb(action(e))(
    response_type_xcb(action(e)),
    0,
    0,
    e.win.id,
    e.win.id,
    0,
    e.location...,
    extent(e.win)...,
    0,
    0,
    0
)

send_event(wm::XWindowManager, e::EventDetails) = send_event(e.win, event_xcb(wm, e))

function send_event(win::XCBWindow, event)
    ref = Ref(event)
    GC.@preserve ref begin
        event_ptr = Base.reinterpret(Cstring, Base.unsafe_convert(Ptr{typeof(event)}, ref))
        @flush @check :error xcb_send_event(win.conn, false, win.id, 0, event_ptr)
    end
end

hex(x) = "0x$(string(x, base=16))"

function send(wm::XWindowManager, win::XCBWindow; location=(0, 0))
    (event; location=location) -> send_event(wm, EventDetails(event, location, floor(time()), win))
end


text:
function augment_states(m::SmetsWoutersOrig{T}, TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T}) where T<:AbstractFloat
    endo = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    exo = m.exogenous_shocks

    n_endo = n_states(m)
    n_exo = n_shocks_exogenous(m)
    @assert (n_endo, n_endo) == size(TTT)
    @assert (n_endo, n_exo) == size(RRR)
    @assert n_endo == length(CCC)

    # Initialize augmented matrices
    n_addl_states = length(endo_addl)
    n_addl_eqs = n_addl_states
    TTT_aug = zeros(n_endo + n_addl_eqs, n_endo + n_addl_states)
    TTT_aug[1:n_endo, 1:n_endo] = TTT
    RRR_aug = [RRR; zeros(n_addl_eqs, n_exo)]
    CCC_aug = [CCC; zeros(n_addl_eqs)]

    ### TTT modifications

    # Track Lags
    TTT_aug[endo_addl[:y_t1], endo[:y_t]] = 1.
    TTT_aug[endo_addl[:c_t1], endo[:c_t]] = 1.
    TTT_aug[endo_addl[:i_t1], endo[:i_t]] = 1.
    TTT_aug[endo_addl[:w_t1], endo[:w_t]] = 1.

    return TTT_aug, RRR_aug, CCC_aug
end


text:
module StanMamba

# package code goes here
using Mamba, Statistics

import CmdStan: convert_a3d
import Mamba: AbstractChains, Chains

include("utilities/convert_a3d.jl")

export
  convert_a3d

end # module


text:
"""
# Module FortranReader

FortranReader provides basic commands to read the 'unformatted' files written
by Fortran programs.

Although not technically portable, most Fortran programs write these files in
a predictable way, in 'records'.  Each record contains either one or several
variables, marked before and after with the length of the variable(s) in bytes.
This marker is written as a 4-byte integer.

function `read_record(::IO, ::DataType, rank::Int, dims)` returns the value in
the next record from the stream and data type given.
"""
module FortranReader

export
    read_record,
    skip_record,
    write_record

# Standard Fortran number sizes
const Real4 = Float32
const Real8 = Float64
const Integer4 = Int32
const Integer8 = Int64
const Complex4 = Complex{Float32}
const Complex8 = Complex{Float64}
const Character = String

const len4 = sizeof(Integer4)
const len8 = sizeof(Integer8)

const permissible_types = (Real4, Real8, Integer4, Integer8, Complex4, Complex8, Character)

"""
    read_record(f::IO, T::DataType, dims...) -> v::T

Return the next record in the stream `f`, in the format of `T`.  This means,
`read_record` assumes that the type of `v` is `T` (e.g., `Float64`).
`dims...` holds the dimensions; the default (no argument) assumes a scalar record.
"""
function read_record(f::IOStream, T::DataType, dims...)
    all(Bool[typeof(dim) <: Integer for dim in dims]) || error("`dims` must be integers")
    d, len = read_raw(f)
    if T <: String
        return String(d)
    end
    Tlen = sizeof(T)
    rank = length(dims)
    if rank == 0
        len == Tlen || error("Record at $(position(f) - len4 - len) is not correct " *
            "length for scalar of type $T")
        reinterpret(T, d)[1]
    elseif rank == 1
        n = len ÷ Tlen
        n != dims[1] && warn("Length of record does not match dimension supplied " *
            "(Requested $(dims[1]); read size $n)")
        reinterpret(T, d)
    elseif rank >= 2
        len ÷ Tlen == prod(dims) || error("Requested dimensions do not match size "*
            "of record.  (Requested $dims (size $(prod(dims))); record size $(len÷Tlen))")
        reshape(reinterpret(T, d), dims)
    else
        error("Arrays must have positive rank")
    end
end

"""
    skip_record(f::IOStream) -> n::Int

Skip over the next record, returning the length of the record.
"""
function skip_record(f::IOStream)
    nstart = reinterpret(Integer4, read(f, len4))[1]
    pos = position(f)
    skip(f, nstart)
    nend = reinterpret(Integer4, read(f, len4))[1]
    nstart == nend || error("Error skipping over record at byte $pos " *
        "in stream $(f.name)")
    nstart
end

"""
    read_raw(f) -> d::Array{UInt8,1}, n::Int

Return the raw data from stream `f`, checked that the start and end records match,
and the number of bytes `n` read.
"""
function read_raw(f::IOStream)
    nstart = reinterpret(Integer4, read(f, len4))[1]
    pos = position(f)
    d = read(f, nstart)
    nend = reinterpret(Integer4, read(f, len4))[1]
    nstart == nend || error("Error reading raw data at position $pos from stream "
        * "$(f.name): record length markers do not match ($nstart != $nend)")
    d, nstart
end

"""
    write_record(f::IOStream, T::DataType, x...) -> n::Int

Write a record to the stream `f` containing `x`, which will be converted to have
element type `T`.  `T` must therefore correspond to one of the Fortran types
`Float{32,64}`, `Int{32,64}`, `Complex{64,128}` or `String`.  You can also use the
equivalent Fortran names which are exported by the module:  `Integer{4,8}`, `Real{4,8}`,
`Complex{4,8}` and `Character`.

Return the **total** number of bytes written to the stream `f`, including the start-
and end-markers of the record.
"""
function write_record(f::IOStream, T::DataType, xs...)
    T in permissible_types || error("Type of record $T is not supported.  " *
        "Choose one of $permissible_types")
    bytes_written = 0
    for x in xs
        n = T == Character ? Integer4(length(x)) : Integer4(sizeof(eltype(T))*length(x))
        write(f, n)
        if T == Character
            write(f, x)
        elseif ndims(x) > 0
            write(f, convert(Array{T}, x))
        else
            write(f, convert(T, x))
        end
        write(f, n)
        bytes_written += 2*len4 + n
    end
    bytes_written
end

end # module


text:
A = [ 1. 2 3; 4 5 6; 7 8 9 ]

B = [12 -51 4; 6 167 -68; -4 24 -41]

function check_DV(A, D, V)
    I = eye(size(A)[1])
    for i = 1:size(A)[1]
        check = (A - D[i] * I) * V[:,i]
        println("For the ", i, "th eigen value and vector, (A - lambda I)^k * v = ", check, "\n\twhich has a L2 norm = ", norm(check))
        assert(norm(check) <= 1e-14)
    end
end

D_A, V_A = eig(A)

check_DV(A, D_A, V_A)

D_B, V_B = eig(B)

check_DV(B, D_B, V_B)

"Power iteration method on A, to find eigenpair with largest eigen value."
function poweriteration(A, maxiter=20, userandomstart=true)
    n = size(A)[1]
    if userandomstart
        X = randn((n, 1))
    else
        X = ones((n, 1))
    end
    for i = 1:maxiter
        next_X = A * X
        X = next_X / norm(next_X)
    end
    # lambda = (X^* ⋅ A ⋅ X) / (X^* ⋅ X)
    lambda = (transpose(conj(X)) * (A * X)) / (transpose(conj(X)) * X)
    lambda, X
end

poweriteration(A)

D_A, V_A = eig(A)

(D_A[1], -1 * V_A[:,1])

poweriteration

println(poweriteration(A, 1))
println(poweriteration(A, 2))
println(poweriteration(A, 5))
println(poweriteration(A, 10))

"Inner product of two vectors v, w. Can be vertical (n,1) or horizontal (1,n) vectors."
function inner(v, w)
    assert(size(v) == size(w))
    nm = size(v)
    if length(nm) == 1
        return transpose(conj(v)) * w
    elseif nm[1] == 1
        return conj(v) * transpose(w)
    end
end

inner([1 1 1], [2 3 4])

inner([1; 1; 1], [2; 3; 4])

"projection(a, e): projection operator of vector a onto base vector e. Uses inner(e, a) and inner(e, e)."
function projection(a, e)
    return (inner(e, a) / inner(e, e)) * e
end

"gramschmidt(A): Gram-Schmidt orthogonalization operator, returns us, es (unormalized matrix, base matrix)."
function gramschmidt(A, verbose=false)
    assert(size(A)[1] == size(A)[2])
    n = size(A)[1]
    if verbose
        println("n = ", n)
    end
    us = zeros((n, n))
    es = zeros((n, n))
    for i = 1:n
        if verbose
            println("i = ", i)
        end
        us[:,i] = A[:,i]
        for j = 1:i-1
            if verbose
                println("\tj = ", j)
                println("\tus[:,i] = ", us[:,i])
            end
            us[:,i] -= projection(A[:,i], us[:,j])
        end
        if verbose
            println("us[:,i] = ", us[:,i])
        end
        es[:,i] = us[:,i] / norm(us[:,i])
        if verbose
            println("es[:,i] = ", es[:,i])
        end
    end
    return us, es
end

A

us, es = gramschmidt(A)

us

es

inner(es[:,1], A[:,1]) * es[:,1]

inner(es[:,1], A[:,2]) * es[:,1] + inner(es[:,2], A[:,2]) * es[:,2]

rank(A)

inner(es[:,1], A[:,3]) * es[:,1] + inner(es[:,2], A[:,3]) * es[:,2] + inner(es[:,3], A[:,3]) * es[:,3]

[ norm(es[:,i]) for i = 1:size(es)[1] ]

(es' * es)[1:2,1:2]

us, es = gramschmidt(B)

[ norm(es[:,i]) for i = 1:size(es)[1] ]

es * es'

V = [3 2; 1 2]

us, es = gramschmidt(V)

es' * es

"QR decomposition, returns (Q, R): Q is orthogonal and R is upper-triangular, such that A = Q R."
function QR(A)
    assert(size(A)[1] == size(A)[2])
    n = size(A)[1]
    us, es = gramschmidt(A)
    Q = copy(es)
    R = zeros((n, n))
    for i = 1:n
        for j = i:n
            R[i,j] = inner(es[:,i], A[:,j])
        end
    end
    return Q, R
end

Q, R = QR(A)

A

Q

R

Q * R

Q2, R2 = QR(B)

B

Q2

R2

Q2 * R2

B == Q2 * R2

assert(norm(B - (Q2 * R2)) <= 1e-13)

"Apply the QR method for maxiter steps. Should return a triangular matrix similar to A (same eigenvalues)."
function QR_method(A, maxiter=50)
    Ak = A
    for k = 1:maxiter
        Qk, Rk = QR(Ak)
        Ak = Rk * Qk
    end
    return Ak
end

Ak = QR_method(A)

"Truncate to zero values under the diagonal (have to be smaller than tolerance)"
function truncate_to_zero_below_diag(A, tolerance=1e-12)
    assert(size(A)[1] == size(A)[2])
    n = size(A)[1]
    for j = 1:n
        for i = j+1:n
            assert(norm(A[i,j]) <= tolerance)
            A[i,j] = 0
        end
    end
    return A
end

truncate_to_zero_below_diag(Ak)

Ak = QR_method(B)

truncate_to_zero_below_diag(Ak)

"Extract the diagonal from the QR method."
function QR_eigenvalues(A, maxiter=50)
    return diag(QR_method(A, maxiter))
end

QR_eigenvalues(A)

eigvals(A)

QR_eigenvalues(B)

eigvals(B)

"Inverse iteration method to find the eigen vector corresponding to a given eigen value."
function inverseIteration(A, val, maxiter=20, userandomstart=true)
    mu_I = val * eye(A)
    inv_A_minus_muI = inv(A - mu_I)
    n = size(A)[1]
    if userandomstart
        X = randn((n, 1))
    else
        X = ones((n, 1))
    end
    for i = 1:maxiter
        next_X = inv_A_minus_muI * X
        X = next_X / norm(next_X)
    end
    X
end

"Approximation of D, V: D contains the eigen values (vector) and V the eigen vectors (column based)."
function QR_eigen(A, maxiter=20, userandomstart=true)
    n = size(A)[1]
    D = QR_eigenvalues(A, maxiter)
    V = zeros((n,n))
    for i = 1:n
        V[:,i] = inverseIteration(A, D[i], maxiter, userandomstart)
    end
    return D, V
end

D2, V2 = QR_eigen(B)

D2s, V2s = eig(B)

max_issues_repo_name: operade/notebooks

max_forks_repo_name: operade/notebooks


text:
# Converts a String to Languages.Language (using STR_TO_LANG)
convert(::Type{L}, lang::S) where {L<:Languages.Language, S<:AbstractString} = begin
    TypeLang = get(STR_TO_LANG, strip(lowercase(lang)), Languages.English)
    return TypeLang()
end

# Converts Languages.Language to String (using LANG_TO_STR)
convert(::Type{S}, lang::Type{L}) where {L<:Languages.Language, S<:AbstractString} =
    get(LANG_TO_STR, lang, string(L))

convert(::Type{S}, lang::L) where {L<:Languages.Language, S<:AbstractString} =
    convert(S, L)

# Overload ismatch to work matching any value within a vector
occursin(r::Regex, strings::T) where T<:AbstractArray{<:AbstractString} =
    any(occursin(r, si) for si in sv);

# Overload lowercase function to work with vectors of strings
lowercase(v::T) where T<:AbstractArray{S} where S<:AbstractString =
    Base.lowercase.(v)

"""
    detect_language(text [; default=DEFAULT_LANGUAGE])

Detects the language of a piece of `text`. Returns a language of
type `Languages.Language`. If the text is empty of the confidence
is low, return the `default` language.
"""
function detect_language(text::AbstractString; default=DEFAULT_LANGUAGE)
    isempty(text) && return default
    detector = LanguageDetector()
    l, _, c = detector(text)  # returns (language, script, confidence)
    if c < 0.15
        return default()
    else
        return l
    end
end

"""
    summarize(sentences [;ns=1, flags=DEFAULT_SUMMARIZATION_STRIP_FLAGS])

Build a summary of the text's `sentences`. The resulting summary will be
a `ns` sentence document; each sentence is pre-procesed using the
`flags` option.
"""
function summarize(sentences::Vector{S};
                   ns::Int=1,
                   flags::UInt32=DEFAULT_SUMMARIZATION_STRIP_FLAGS
                  ) where S<:AbstractString
    # Get document term matrix
    s = prepare.(sentences, flags)
    filter!(d->occursin(r"[a-zA-Z0-9]", d), s)
    t = dtm(DocumentTermMatrix{Float32}(s))
    tf_idf!(t)
    # Page rank
    α = 0.85  # damping factor
    n = 100  # number of iterations
    ϵ = 1.0e-6  # convergence threhshold
    G = Graph(t' * t)
    try
        p = pagerank(G, α, n, ϵ)
        # Sort sentences and return
        text_summary = sentences[sort(sortperm(p, rev=true)[1:min(ns, length(p))])]
        return text_summary
    catch
        @warn "Summarization failed during TextRank. No summarization done."
        return sentences
    end
end

"""
Post-processes a string to fit a certain length, adding … if necessary
at the end of its choped represenation.
"""
function chop_to_length(input, len)
     input = replace(input, "\n" => "")
     idxs = collect(eachindex(input))
     _idx = findlast(Base.:<=(len), idxs)
     if _idx == nothing
         _len=0
     else
         _len = idxs[findlast(Base.:<=(len), idxs)]
     end
     length(input) > len ? input[1:_len] * "…"  : input
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

@testset "test integration point" begin
    ip = IP(1, 1.0, sqrt(1.0/3.0)*[-1.0, -1.0])
    strain = [1.0 2.0; 3.0 4.0]
    update!(ip, "strain", 0.0 => strain)
    @test isapprox(ip("strain", 0.0), strain)
    @test isapprox(ip("strain"), strain)
end


text:
"""
Inplace sde solver.

Typical call is like
solve!(EulerMaruyama!(), X, x0, W, P)
or
solve!(StratonovichHeun!(), X, x0, W, P)

Here
- x0 is the starting point,
- W the driving Wiener process,
- P the ContinuousTimeProcess

Requires
- drift function _b!((i,t¯), y, tmp1, P)
- diffusion function σ!(t¯, y, dw, tmp2, P)

The result is written into X
"""
function solve!(solver::EulerMaruyama!, Y::VSamplePath, u::T, W, P::ProcessOrCoefficients) where {T}
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    tt = Y.tt
    tt[:] = W.tt
    yy = Y.yy
    y::T = copy(u)

    size(Y.yy) != (length(y), N) && error("Starting point has wrong length.")
    #assert(size(W.yy) == (length(y), N))
    tmp1 = copy(y)
    tmp2 = copy(y)
    dw = W.yy[.., 1]
    dw2 = copy(dw)
    for i in 1:N-1
        t¯ = tt[i]
        dt = tt[i+1] - t¯
        for k in eachindex(tmp1)
            @inbounds yy[k, i] = y[k]
        end
        for k in eachindex(dw)
            @inbounds dw[k] = W.yy[k, i+1] - W.yy[k, i]
        end
        _b!((i,t¯), y, tmp1, P)
        σ!(t¯, y, dw, tmp2, P)
        for k in eachindex(y)
            @inbounds y[k] = y[k] + tmp1[k]*dt + tmp2[k]
        end
    end
    yy[.., N] = endpoint(y, P)
    Y
end

function Bridge.solve!(solver::EulerMaruyama!, Y::SamplePath, u, W::SamplePath, P::Bridge.ProcessOrCoefficients)
    N = length(W)

    tt = W.tt
    yy = Y.yy
    copyto!(yy[1], u)

    i = 1
    dw = W.yy[i+1] - W.yy[i]
    t¯ = tt[i]
    dt = tt[i+1] - t¯

    tmp1 = Bridge._b!((i,t¯), u, 0*u, P)
    tmp2 = Bridge.σ(t¯, u, dw, P)

    y = u + tmp1*dt + tmp2

    for i in 2:N-1
        t¯ = tt[i]
        dt = tt[i+1] - t¯
        copyto!(yy[i], y)
        if dw isa Number
            dw = W.yy[i+1] - W.yy[i]
        else
            for k in eachindex(dw)
                dw[k] = W.yy[i+1][k] - W.yy[i][k]
            end
        end

        Bridge._b!((i,t¯), y, tmp1, P)
        Bridge.σ!(t¯, y, dw, tmp2, P)

        for k in eachindex(y)
            y[k] = y[k] + tmp1[k]*dt + tmp2[k]
        end
    end
    copyto!(yy[end], Bridge.endpoint(y, P))
    SamplePath(tt, yy)
end


text:
using BenchmarkTools
using DataFrames
using Statistics
using Transducers
using VegaLite

resultpath = joinpath(@__DIR__, "result.json")
result, = BenchmarkTools.load(resultpath)

df_raw =
    BenchmarkTools.leaves(result) |>
    Map() do ((nworks, ex), trial)
        (nworks = parse(Int, nworks), executor = ex, trial = trial)
    end |>
    DataFrame
#-

begin
    df_tmp = select(df_raw, Not(:trial))
    df_tmp[!, :minimum] = map(trial -> minimum(trial).time, df_raw.trial)
    df_tmp[!, :median] = map(trial -> median(trial).time, df_raw.trial)
    df_tmp[!, :memory] = map(trial -> trial.memory, df_raw.trial)
    df_stats = stack(
        df_tmp,
        [:minimum, :median],
        variable_name = :time_stat,
        value_name = :time_ns,
    )
end
#-

df = combine(groupby(df_stats, [:nworks, :time_stat])) do group
    d = Dict(zip(group.executor, group.time_ns))
    (speedup = d["ThreadedEx"] / d["WorkStealingEx"],)
end
#-

plt1 = @vlplot(
    layer = [
        {
            mark = {type = :line, point = true},
            encoding = {
                x = {field = :nworks},
                y = {field = :speedup, axis = {title = "Speedup (T_default / T_WS)"}},
                color = {field = :time_stat},
            },
        },
        {
            mark = :rule,
            encoding = {y = {datum = 1}},
        },
    ],
    data = df,
    width = 400,
    height = 200,
)
#-

plt2 = @vlplot(
    mark = {type = :line, point = true},
    x = :nworks,
    y = {field = :time_ns, axis = {title = "Time [ns]"}},
    color = {field = :time_stat},
    column = :executor,
    data = df_stats,
)
#-


text:
using DynamicHMCModels, StatsFuns

ProjDir = @__DIR__

delim = ';'
df = CSV.read(joinpath(ProjDir, "..", "..", "data", "chimpanzees.csv"), DataFrame; delim)

Base.@kwdef struct Chimpanzees_02
    "Number of actors"
    N_actors::Int
    pulled_left::Vector{Int}
    prosoc_left::Vector{Int}
    condition::Vector{Int}
    actor::Vector{Int}
end

function make_transformation(model::Chimpanzees_02)
    as((a = as(Vector, model.N_actors), bp = asℝ, bpC = asℝ))
end

model = Chimpanzees_02(; N_actors = maximum(df.actor), pulled_left = df.pulled_left,
                    prosoc_left = df.prosoc_left, condition = df.condition,
                    actor = df.actor)

function (model::Chimpanzees_02)(θ)
    @unpack pulled_left, prosoc_left, condition, actor = model
    @unpack a, bp, bpC = θ
    ℓ_likelihood = mapreduce(+, actor, condition, prosoc_left,
       pulled_left) do actor, condition, prosoc_left, pulled_left
           p = logistic(a[actor] + (bp + bpC * condition) * prosoc_left)
           logpdf(Bernoulli(p), pulled_left)
       end
    P = Normal(0, 10)
    ℓ_prior = logpdf(P, bpC) + logpdf(P, bp) + sum(a -> logpdf(P, a), a)
    ℓ_prior + ℓ_likelihood
end

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P)
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

p = as_particles(posterior)
display(p)

# Result rethinking

rethinking = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
        Mean        SD       Naive SE       MCSE      ESS
a.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000
a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000
a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000
a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000
a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000
a.6  0.21599365 0.26307574 0.0041595927 0.0045153523 1000
a.7  1.81090866 0.39318577 0.0062168129 0.0071483527 1000
 bp  0.83979926 0.26284676 0.0041559722 0.0059795826 1000
bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000
";


text:
import Base.mean
export mean

function mean(x::AffineExpr)
  return sum(x) / (x.size[1] * x.size[2])
end

function mean(x::AffineExpr, dim::Int64)
  return sum(x, dim) / x.size[dim]
end


text:
@testset "SOM training helper functions" begin

    @testset "gridRectangular" begin
        grid = GigaSOM.gridRectangular(5, 5)
        @test size(grid) == (25, 2)
    end

    @testset "Kernels" begin
        @test isapprox(
            gaussianKernel(Vector{Float64}(1:7), 2.0),
            [0.1760326, 0.1209853, 0.0647587, 0.0269954, 0.0087641, 0.0022159, 0.0004363],
            atol = 1e-4,
        )
        @test isapprox(
            bubbleKernel(Vector{Float64}(1:6), 5.0),
            [0.979796, 0.916515, 0.8, 0.6, 0, 0],
            atol = 0.001,
        )
        @test isapprox(
            thresholdKernel(Vector{Float64}(1:10), 5.0),
            [1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
            atol = 1e-3,
        )
        @test isapprox(
            thresholdKernel(Vector{Float64}(1:10), 10.1),
            [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0],
            atol = 1e-3,
        )
    end

    @testset "Radii-generating functions" begin
        radii1 = expRadius().(5.0, 0.5, Vector(1:10), 10)
        radii2 = expRadius(-10.0).(10.0, 0.1, Vector(1:20), 20)
        radii3 = linearRadius.(50.0, 0.001, Vector(1:30), 30)
        @test isapprox(radii1[1], 5.0)
        @test isapprox(radii1[10], 0.5)
        @test isapprox(radii2[1], 10.0)
        @test isapprox(radii2[20], 0.1)
        @test isapprox(radii3[1], 50.0)
        @test isapprox(radii3[30], 0.001)
        @test all(isapprox.(radii1[1:9] ./ radii1[2:10], radii1[1] / radii1[2]))
        #note: radius2 is adjusted and thus not really exponential
        @test all(isapprox.(radii3[1:29] .- radii3[2:30], radii3[1] - radii3[2]))
    end

    @testset "distMatrix" begin
        g = GigaSOM.gridRectangular(2, 2)
        dm = GigaSOM.distMatrix(Euclidean())(g)
        @test size(dm) == (4, 4)
        @test all([dm[i, i] == 0 for i = 1:4])
    end
end


text:
## =============================================================================
# Hidden markov model with zero inflated likelihood

struct ZINegBinomialProblem{Y,N,I,P} <: SurveyProblem
    y::Y
    n::N
    indices::I
    prior::P
end

function ZINegBinomialProblem(y,n,indices; c=nothing)
    prior = ZINegBinomialPrior1(n; c)
    ZINegBinomialProblem(y,n,indices, prior)
end

struct ZINegBinomialPrior1{P,L,F,N}
    prior_p::P
    prior_μ::L
    prior_ϕ::F
    n::N
end

function ZINegBinomialPrior1(n; c=nothing)
    # if isnothing(c)
    prior_p = Beta()
    prior_μ = filldist(truncated(GammaMeanCv(15,.8), 0, 1000), n)
    prior_ϕ = filldist(truncated(Normal(0, 20), 0, Inf), n)
    return ZINegBinomialPrior1( prior_p, prior_μ, prior_ϕ, n)
    # elseif c == :family
    #     prior_p  = truncated(Normal(0.51,0.05), 0, 1)
    #     prior_μ = truncated(Normal(3.5,1), 0, Inf)
    #     prior_ϕ = filldist(truncated(Normal(0, 20), 0, Inf), n)
    #     return ZINegBinomialPrior1( prior_p, prior_μ, prior_ϕ, n)
    # elseif c == :colleagues
    #     prior_p  = truncated(Normal(0.33,0.05), 0, 1)
    #     prior_μ = truncated(Normal(5,1.5), 0, Inf)
    #     prior_ϕ = filldist(truncated(Normal(0, 20), 0, Inf), n)
    #     return ZINegBinomialPrior1( prior_p, prior_μ, prior_ϕ, n)
    # elseif c == :friends
    #     prior_p  = truncated(Normal(0.46,0.05), 0, 1)
    #     prior_μ = truncated(Normal(4.5,2), 0, Inf)
    #     prior_ϕ = filldist(truncated(Normal(0, 20), 0, Inf), n)
    #     return ZINegBinomialPrior1( prior_p, prior_μ, prior_ϕ, n)
    # elseif c == :strangers
    #     prior_p  = truncated(Normal(0.32,0.05), 0, 1)
    #     prior_μ = truncated(Normal(20,8), 0, Inf)
    #     prior_ϕ = filldist(truncated(Normal(0, 20), 0, Inf), n)
    #     return ZINegBinomialPrior1( prior_p, prior_μ, prior_ϕ, n)
    # else
    #     error( "wrong choice" )
    # end
end

function Base.rand(priors::ZINegBinomialPrior1)
    @unpack prior_p, prior_μ, prior_ϕ, n = priors
    p  = rand(prior_p)
    μ  = rand(prior_μ)
    ϕ  = rand(prior_ϕ)
    return (; p, μ, ϕ)
end

function Distributions.logpdf(prior::ZINegBinomialPrior1, θ)
    @unpack p, μ, ϕ = θ
    @unpack prior_p, prior_μ, prior_ϕ = prior
    ℓ  = logpdf(prior_p, p)
    ℓ += logpdf(prior_μ, μ)
    ℓ += logpdf(prior_ϕ, ϕ)
    return ℓ
end

# likelihood_dist(p, μ, ϕ) =
#     ZeroInflatedDistribution.(
#         Ref(Bernoulli(p)),
#         NegativeBinomial2.(μ, Ref(ϕ)
#     )
# )
likelihood_dist(p, μ, ϕ) =
    ZeroInflatedDistribution.(
        Ref( Bernoulli(p) ),
        NegativeBinomial2.(μ, ϕ .+ 0.1)
)

function (problem::ZINegBinomialProblem)(θ)
    @unpack p, μ, ϕ  = θ               # extract the parameters
    @unpack y, indices, prior = problem
    # !(ϕ > 0) && (@error "not ϕ>0"; return -Inf)
    # !(0 <= p <= 1) && (@error "not 0 <= p <= 1"; return -Inf)

    ℓ  = logpdf(prior, θ)
    # !all(0 .<= μ ) && (@error "not all(0 .<= μ )"; return -Inf)

    d = likelihood_dist(p, μ, ϕ)
    ℓ += logpdf(arraydist(d[indices]),y)
    return ℓ
end

function problem_transformation(p::ZINegBinomialProblem)
    as((p = as𝕀, μ = as(Vector, as(Real, 0, 1000), p.n), ϕ = as(Vector, asℝ₊, p.n)))
end

function ZINegBinomialProblem(data, col)
    y = data[:,col]
    indices = date2index( data.date )
    n = last(indices)
    return ZINegBinomialProblem(y, n, indices; c=col)
end

function StatsBase.predict(problem::ZINegBinomialProblem, θ)
    @unpack p, μ, ϕ = θ               # extract the parameters
    @unpack y, indices = problem

    d = likelihood_dist(p, μ, ϕ)
    y = rand(arraydist(d[indices]))
    return y
end

generate_quantities(problem::ZINegBinomialProblem, θ) = θ.μ

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
#=
  Euler #24 in Julia.

  Euler 24:
  """
  A permutation is an ordered arrangement of objects. For example, 3124 is one
  possible permutation of the digits 1, 2, 3 and 4. If all of the permutations are
  listed numerically or alphabetically, we call it lexicographic order. The
  lexicographic permutations of 0, 1 and 2 are:

      012   021   102   120   201   210

  What is the millionth lexicographic permutation of the digits
  0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

#
function euler24a()
    a = Array(0:9)
    c = 1
    while c < 1000000
        a = next_permutation(a)
        c += 1
    end

    return parse(Int, join(a,""))
end

# Too slow: 15.24007616s
function euler24b()
    return parse(Int,join(all_permutations2(Array(0:9))[1000000]))
end

# Inspired by a solution on the 'net
# 0.05707951s
function euler24c()
    n = 999999
    p = 10
    eli = (1:p).|>(i->i % 10)
    answer = []
    for i in 1:p-1
        f = factorial(p-i)
        d = floor(Int,n / f)
        n %= f
        push!(answer,eli[d])
        # Note the indexOf()...
        # eli.splice(eli.indexOf(eli[d-1]),1);
        eli = filter(e->e != eli[d], eli)
    end
    return parse(Int,join(vcat(answer,eli)))
end

# all_permutations is much faster than all_permutations2
# 0.49170683s
function euler24d()
    return parse(Int,join(all_permutations(Array(0:9))[1000000]))
end

run_euler(euler24a)
# run_euler(euler24b) # too slow
# run_euler(euler24c)
# run_euler(euler24d)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
#functions for calling the matrix vector products
#i.e. function A(y_in,Nx,Ny,h) computes the product A*u and stores it in y

#Right now only_out[for  the second-order accurate SBP operators for constant coefficients from
#Mattsson and Nordstrom, 2004.

#Written for the 2D domain (x,y) \in (a, b) \times (c, d),
#stacking the grid function in the vertical direction,
#so that, e.g., \partial u / \partial x \approx D \kron I and
#\partial u / \partial y_out[\approx I \kron D
#the faces of the 2D domain are label bottom to top, left to right, i.e.
#side 1 is y_out[= c, side 2 is y_out[= d
#side 3 is x = a, side 4 is x = b

#D2x and D2y_out[compute approximations to
#\partial^2 u / \partial x^2 and
#\partial^2 u / \partial y^2, respectively

#FACEtoVOL(face,y_in_face,Nx,Ny) maps the face value y_in_face to a full length vector at the
#nodes corresponding to face

#BxSx \approx the traction on faces 3 and 4
#BySy_out[\approx the traction on faces 1 and 2

#BxSx_tran and BySy_tran are the transposes of BxSx and BySy

function D2x_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)
	idx = 1:Ny
	y_out[idx] = (y_in[idx] - 2 .* y_in[Ny.+ idx] + y_in[2*Ny .+ idx]) ./ h^2

	idx1 = Ny+1:N-Ny
	y_out[idx1] = (y_in[idx1 .- Ny] - 2 .* y_in[idx1] + y_in[idx1 .+ Ny]) ./ h^2

	idx2 = N-Ny+1:N
	y_out[idx2] = (y_in[idx2 .- 2*Ny] -2 .* y_in[idx2 .- Ny] + y_in[idx2]) ./ h^2

	return y_out
end

function D2y_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)
	idx = 1:Ny:N-Ny+1
	y_out[idx] = (y_in[idx] - 2 .* y_in[idx .+ 1] + y_in[idx .+ 2]) ./ h^2

	idx1 = Ny:Ny:N
	y_out[idx1] = (y_in[idx1 .- 2] - 2 .* y_in[idx1 .- 1] + y_in[idx1]) ./ h^2

	for j = 1:Nx
		idx = 2+(j-1)*Ny:j*Ny-1
		y_out[idx] = (y_in[idx .- 1] - 2 .* y_in[idx] + y_in[idx .+ 1]) ./ h^2
	end

	return y_out

end

function Dx_p4(y_in, Nx, Ny, h)
	d  = [1/12 -2/3 0 2/3 -1/12]
    bd = [-24/17  59/34  -4/17  -3/34  0     0;
           -1/2    0      1/2    0     0     0;
            4/43 -59/86   0     59/86 -4/43  0;
            3/98   0    -59/98   0    32/49 -4/49]
	N = Nx*Ny
	y_out = zeros(N)

	for row_num = 1:size(bd)[1]
		idx = (row_num - 1)*Ny+1:(row_num)*Ny
		y_out[idx] = (bd[row_num,1] * y_in[idx .+ (1 - row_num)*Ny] + bd[row_num,2] * y_in[idx .+ Ny*(2-row_num)] + bd[row_num,3] * y_in[idx .+ Ny*(3-row_num)] + bd[row_num,4] * y_in[idx .+ Ny*(4-row_num)] + bd[row_num,5] * y_in[idx .+ Ny*(5-row_num)] + bd[row_num,6] * y_in[idx .+ Ny*(6-row_num)]) ./ h
	end

	idx1 = (size(bd)[1])*Ny+1:N-(size(bd)[1])*Ny
	y_out[idx1] = (d[1]*y_in[idx1 .- 2*Ny] + d[2]*y_in[idx1 .- Ny] + d[3]*y_in[idx1] + d[4]*y_in[idx1 .+ Ny] + d[5] * y_in[idx1 .+ 2*Ny]) ./ (h)

	for row_num = 1:size(bd)[1]
		idx2 = N-(row_num)*Ny+1:N .- (row_num - 1)*Ny
		y_out[idx2] = -(bd[row_num,1] * y_in[idx2 .+ (row_num - 1)*Ny] + bd[row_num,2] * y_in[idx2 .- Ny*(2-row_num)] + bd[row_num,3] * y_in[idx2 .- Ny*(3-row_num)] + bd[row_num,4]*y_in[idx2 .- Ny*(4-row_num)] + bd[row_num,5] * y_in[idx2 .- Ny*(5-row_num)] + bd[row_num,6] * y_in[idx2 .- Ny*(6-row_num)]) ./ h
	end

	return y_out
end

function Dy_p4(y_in, Nx, Ny, h)
	d  = [1/12 -2/3 0 2/3 -1/12]
    bd = [-24/17  59/34  -4/17  -3/34  0     0;
           -1/2    0      1/2    0     0     0;
            4/43 -59/86   0     59/86 -4/43  0;
            3/98   0    -59/98   0    32/49 -4/49]
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = (y_in[idx .+ 1] - y_in[idx]) ./ h

	idx = Ny:Ny:N
	y_out[idx] = (y_in[idx] - y_in[idx .- 1]) ./h

	for j = 1:Nx
		idx = 2+(j-1)*Ny:j*Ny-1
		y_out[idx] = (y_in[idx .+ 1] - y_in[idx .- 1]) ./ (2*h)
	end

	return y_out
end

function Hxinv_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	idx = Ny+1:N-Ny
	y_out[idx] = (1*y_in[idx]) .* (1/h)

	idx = N-Ny+1:N
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	return y_out
end

function Hyinv_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	idx = Ny:Ny:N
	y_out[idx] = (2*y_in[idx]) .* (1/h)

	for i = 1:Nx
		idx = 2+(i-1).*Ny:i*Ny-1
		y_out[idx] = y_in[idx] .* (1/h)
	end

	return y_out

end

function Hx_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
    y_out = zeros(N)

    idx = 1:Ny
	y_out[idx] = h .* (1/2)*y_in[idx]

    idx = Ny+1:N-Ny
    y_out[idx] = h .* 1*y_in[idx]

    idx = N-Ny+1:N
	y_out[idx] = h .* (1/2)*y_in[idx]

    return y_out

end

function Hy_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
    y_out = zeros(N)

    idx = 1:Ny:N-Ny+1
	y_out[idx] = h .* (1/2)*y_in[idx]

    idx = Ny:Ny:N
	y_out[idx] = h .* (1/2)*y_in[idx]

    for i = 1:Nx
        idx = 2+(i-1).*Ny:i*Ny-1
        y_out[idx] = h .* y_in[idx]
    end

    return y_out

end

function FACEtoVOL(y_in_face, face, Nx, Ny)
	N = Nx*Ny
	y_out = zeros(N)

	if face == 1
		idx = 1:Ny:N-Ny+1
	elseif face == 2
		idx = Ny:Ny:N
	elseif face == 3
		idx = 1:Ny
	elseif face == 4
		idx = N-Ny+1:N
	else
	end

	y_out[idx] = y_in_face

	return y_out

end

function VOLtoFACE(y_in, face, Nx, Ny)
	N = Nx*Ny
    y_out = zeros(N)

    if face == 1
        idx = 1:Ny:N-Ny+1
    elseif face == 2
        idx = Ny:Ny:N
    elseif face == 3
        idx = 1:Ny
    elseif face == 4
        idx = N-Ny+1:N
    else
    end

	y_out[idx] = y_in[idx]
    return y
end

function Bx(Nx,Ny)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = -1 .* ones(Ny)

	idx = N-Ny+1:N
	y_out[idx] = 1 .* ones(Ny)
	return y_out
end

function By(Nx,Ny)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = -1 .* ones(Ny)

	idx = Ny:Ny:N
	y_out[idx] = 1 .* ones(Ny)
	return y_out
end

function BxSx_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny
	y_out[idx] = (1/h) .* (1.5 .* y_in[idx] - 2 .* y_in[idx .+ Ny] + 0.5 .* y_in[idx .+ 2*Ny])
	y_out[N-Ny.+ idx] = (1/h) .* (0.5 .* y_in[N-3*Ny .+ idx] - 2 .* y_in[N-2*Ny .+ idx] + 1.5 .* y_in[N-Ny .+ idx])

	return y_out

end

function BySy_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx = 1:Ny:N-Ny+1
	y_out[idx] = (1/h) .* (1.5 .* y_in[idx] - 2 .* y_in[idx .+ 1] + 0.5 .* y_in[idx .+ 2])

	idx = Ny:Ny:N
	y_out[idx] = (1/h) .* (0.5 .* y_in[idx .- 2] - 2 .* y_in[idx .- 1] + 1.5 .* y_in[idx])

	return y_out
end

function BxSx_tran_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx1 = 1:Ny
	y_out[idx1] += (1.5 .* y_in[idx1]) .* (1/h)
	idx = Ny+1:2*Ny
	y_out[idx] += (-2 .* y_in[idx1]) .* (1/h)
	idx  = 2*Ny+1:3*Ny
	y_out[idx] += (0.5 .* y_in[idx1]) .* (1/h)

	idxN = N-Ny+1:N
	y_out[idxN] += (1.5 .* y_in[idxN]) .* (1/h)
	idx = N-2*Ny+1:N-Ny
	y_out[idx] += (-2 .* y_in[idxN]) .* (1/h)
	idx = N-3*Ny+1:N-2*Ny
	y_out[idx] += (0.5 .* y_in[idxN]) .* (1/h)

	return y_out
end

function BySy_tran_p4(y_in, Nx, Ny, h)
	N = Nx*Ny
	y_out = zeros(N)

	idx1 = 1:Ny:N-Ny+1
	y_out[idx1] += (1.5 .* y_in[idx1]) .* (1/h)
	idx = 2:Ny:N-Ny+2
	y_out[idx] += (-2 .* y_in[idx1]) .* (1/h)
	idx = 3:Ny:N-Ny+3
	y_out[idx] += (0.5 .* y_in[idx1]) .* (1/h)

	idxN = Ny:Ny:N
	y_out[idxN] += (1.5 .* y_in[idxN]) .* (1/h)
	idx = Ny-1:Ny:N-1
	y_out[idx] += (-2 .* y_in[idxN]) .* (1/h)
	idx = Ny-2:Ny:N-2
	y_out[idx] += (0.5 .* y_in[idxN]) .* (1/h)

	return y_out
end


text:
using DSGE, StateSpaceRoutines
using QuantEcon: solve_discrete_lyapunov
using DataFrames

# Setup the model and data

m = AnSchorfheide()
df = readtable("us.txt", header = false, separator = ' ')
data = convert(Matrix{Float64}, df)'

params = [2.09, 0.98, 2.25, 0.65, 0.34, 3.16, 0.51, 0.81, 0.98, 0.93, 0.19, 0.65, 0.24,
          0.115985, 0.294166, 0.447587]
update!(m, params)

# Solution to a Linear DSGE Model w/ IID Gaussian Errors

system  = compute_system(m)
TTT = system[:TTT]
RRR = system[:RRR]
QQ  = system[:QQ]
Φ, Ψ, F_ϵ, F_u = compute_system_function(system)

# Tuning of the tempered particle filter algorithm

tuning = Dict(:r_star => 2.,
              :xtol => 1e-3, :resampling_method => :systematic,
              :n_particles => 1000, :n_presample_periods => 0,
              :allout => true, :parallel => false)

# Generation of the initial state draws

n_states = n_states_augmented(m)
s0 = zeros(n_states)
P0 = solve_discrete_lyapunov(TTT, RRR*QQ*RRR')
U, E, V = svd(P0)
s_init = s0 .+ U*diagm(sqrt.(E))*randn(n_states, tuning[:n_particles])
tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; tuning...)


text:
julia> using Statistics; mean([1,2,3])
2.0
julia> mean(1:10)
5.5
julia> mean([])
ERROR: mean of empty collection undefined: []

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#This Julia file has been generated by rosetta script of MatlabCompat library from an m-file
#The code generated need further corrections by you. Execute it line by line and correct the errors.
#
#
# Copyright © 2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich
#
#  MIT license.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including without
# limitation the rights to use, copy, modify, merge, publish, distribute,
# sublicense, and/or sell copies of the Software, and to permit persons to whom
# the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using MatlabCompat
using MatlabCompat.ImageTools.Morph

# Morphological stuff
# -------------------------------------------------
se = strel("disk", 1);
# Flat STREL object containing 5 neighbors.
#
# Neighborhood:
#      0     1     0
#      1     1     1
#      0     1     0

# imopen is an erosion followed by a dilation
imgErodedDilated = imopen(I, se);
imshow(imgErodedDilated);
# -------------------------------------------------


text:
# struct PlanarVector{T} <: AbstractPlanarVector{T}
#     u::T
#     v::T
# end

"""
    struct CartesianVector{T} <: AbstractCoordinateVector{T, Cartesian}

Describes a three-dimensional vector in Cartesian coordinates.

## Fields
* `x`: x-coordinate (in m).
* `y`: y-coordinate (in m).
* `z`: z-coordinate (in m).

See also [`CylindricalVector`](@ref).
"""
struct CartesianVector{T} <: AbstractCoordinateVector{T, Cartesian}
    x::T
    y::T
    z::T
end

zero(VT::Type{<:AbstractCoordinateVector{T}}) where {T} = VT(zero(T),zero(T),zero(T))

# @inline rotate(pt::CartesianPoint{T}, r::RotMatrix{3,T,TT}) where {T, TT} = r.mat * pt
# @inline rotate(pt::CylindricalPoint{T}, r::RotMatrix{3,T,TT}) where {T, TT} = CylindricalPoint(rotate(CartesianPoint(pt), r))
# @inline rotate!(vpt::Vector{<:AbstractCoordinatePoint{T}}, r::RotMatrix{3,T,TT}) where {T, TT} = begin for i in eachindex(vpt) vpt[i] = rotate(vpt[i], r) end; vpt end
# @inline rotate!(vvpt::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, r::RotMatrix{3,T,TT}) where {T, TT} = begin for i in eachindex(vvpt) rotate!(vvpt[i], r) end; vvpt end

# @inline scale(pt::CartesianPoint{T}, v::SVector{3,T}) where {T} = pt .* v
# @inline scale(pt::CylindricalPoint{T}, v::SVector{3,T}) where {T} = CylindricalPoint(scale(CartesianPoint(pt),v))
# @inline scale!(vpt::Vector{<:AbstractCoordinatePoint{T}}, v::SVector{3,T}) where {T} = begin for i in eachindex(vpt) vpt[i] = scale(vpt[i], v) end; vpt end
# @inline scale!(vvpt::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, v::SVector{3,T}) where {T} = begin for i in eachindex(vvpt) scale!(vvpt[i], v) end; vvpt end

# @inline translate(pt::CartesianPoint{T}, v::CartesianVector{T}) where {T} = pt + v
# @inline translate(pt::CylindricalPoint{T}, v::CartesianVector{T}) where {T} = CylindricalPoint(translate(CartesianPoint(pt), v))
# @inline translate!(vpt::Vector{<:AbstractCoordinatePoint{T}}, v::CartesianVector{T}) where {T} = begin for i in eachindex(vpt) vpt[i] = translate(vpt[i], v) end; vpt end
# @inline translate!(vvpt::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, v::CartesianVector{T}) where {T} =  begin for i in eachindex(vvpt) translate!(vvpt[i], v) end; vvpt end

"""
    struct CylindricalVector{T} <: AbstractCoordinateVector{T, Cylindrical}

Describes a three-dimensional vector in cylindrical coordinates.

## Fields
* `r`: Radius (in m).
* `φ`: Polar angle (in rad).
* `z`: `z`-coordinate (in m).

!!! note
    `φ == 0` corresponds to the `x`-axis in the Cartesian coordinate system.

See also [`CartesianVector`](@ref).
"""
struct CylindricalVector{T} <: AbstractCoordinateVector{T, Cylindrical}
    r::T
    φ::T
    z::T
end


text:
export det
export minor
export cofactor

sub(M::Array{Float64}, i, j) = M[1:end .!= i, 1:end .!= j]
det(a::Array{Float64}) = begin
    if size(a) == (2, 2)
        a[1, 1] * a[2, 2] - a[1, 2] * a[2, 1]
    else
        sum([(-1)^(i+1) * a[i, 1] * det(sub(a, i, 1)) for i in 1:size(a, 1)])
        # performance tip, look for the best row/column choice for expansion.
    end
end
minor(M::Array{Float64}, i, j) = det(sub(M, i, j))
cofactor(M::Array{Float64}, i, j) = minor(M, i, j) * (-1)^(i+j)

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
# Write function to automate plotting scale-location plot
# - Plots standard scale-location plot for assessing linear model fits

function plot_scale_location(;model, df, y)

    # Uses an if-else statement to either plot Standardised Pearson residuals
    # if the model is fit with a Gaussian/Normal distribution, or
    # Standardised deviance residuals if the model is fit with an
    # exponential family (e.g. Poisson) or any non-Gaussian data distribution.

    ###################################################
    # - Setup and defaults
    ###################################################

    # Set Gadfly theme
    default_theme = Gadfly.Theme(
        panel_fill = nothing,
        highlight_width = 0Gadfly.pt,
        point_size = 0.8Gadfly.mm,
        key_position = :inside,
        #grid_line_width = 0,
        panel_stroke = Colors.colorant"black",
        background_color = Colors.colorant"white")
    Gadfly.push_theme(default_theme)

    ###################################################
    # - Plot either Pearson or deviance residuals
    ###################################################

    # Extract statistical distribution specified in model
    distType = typeof(model).parameters[1].parameters[1].parameters[2]
    distType = string.(distType)

    # Calculate and plot Standardised Pearson Residuals if model was
    # fit with a Gaussian distribution
    if distType == "Normal{Float64}"

        ###################################################
        # - Calculate raw Pearson residuals
        ###################################################

        # Extract fitted/predicted values from model object
        pred = GLM.predict(model)

        # Extract vector of response values
        resp = df[!, y]

        # Calculate raw residuals
        PearsResids = resp .- pred

        # Calculate standard deviation of observed value
        sdResp = Statistics.std(resp)

        # Calculate length (n) of response vector
        nResp = Base.length(resp)

        # Calculate standard error of response
        #seResp = sdResp ./ Base.sqrt.(nResp)

        # Calculate Pearson residuals
        #PearsResids = (resids ./ seResp) ./ 10

        ###################################################
        # - Calculate hat matrix
        ###################################################

        # Extract model formula
        mod_formula = StatsModels.formula(model)

        # Calculate hat matrix
        reg = LinRegOutliers.createRegressionSetting(mod_formula, df)
        hat = LinRegOutliers.hatmatrix(reg)

        ###################################################
        # - Standardise Pearson residuals
        ###################################################

        # Extract deviance and residual df
        modDeviance = GLM.deviance(model)
        modResidualDf = GLM.dof_residual(model)

        # Calculate 's'
        s = Base.sqrt.(modDeviance ./ modResidualDf)

        # Standardise residuals by hat matrix
        # - Same as R: resid(mod1, type="deviance")/sqrt(1 - hatvalues(mod1))
        # - Necessary when dispersion parameter not equal to 1
        dHat = LinearAlgebra.Diagonal(hat)
        StdPearsResids = PearsResids ./ (s .* (Base.sqrt.(1 .- dHat)) )
        StdPearsResids = StdPearsResids[:,1]

        # Calculate square root of Std. Pearson residuals
        sqrPearResids = sqrt.(complex(StdPearsResids))
        sqrPearResids = abs.(sqrPearResids)

        # Define quantiles
        qx = Distributions.quantile.(Distributions.Normal(),
                                        Base.range(0.5,
                                                  stop = (nResp .- 0.5),
                                                  length = (nResp)) ./ (nResp .+ 1))

        ###################################################
        # - Plot scale-location plot
        ###################################################

        p = Gadfly.plot(
            # Add points layer
            x = pred,
            y = sqrPearResids,
            Gadfly.layer(Gadfly.Geom.smooth(method = :loess,
                                    smoothing = 0.9),
                                    color=[Colors.colorant"red"]),
            Gadfly.layer(Gadfly.Geom.point),
            # Change plot aesthetics
            Gadfly.Scale.x_continuous(format=:plain),
            Gadfly.Guide.title("Scale-Location"),
            Gadfly.Guide.xlabel("Fitted values",
                                orientation=:horizontal),
            Gadfly.Guide.ylabel("√ Std. Pearson residuals"))

        # Return the plot object
        return(p)

    # Calculate and plot Standardised deviance residuals if model was
    # fit with a non-Gaussian distribution
    else

        ###################################################
        # - Calculate raw deviance residuals
        ###################################################

        # Use if-elseif-else statements to calculate deviance residuals depending
        # on which error distribution was specified during model specification

        # Extract fitted/predicted values from model object
        pred = GLM.predict(model)

        # Calculate raw deviance residuals
        devResids = ModelCheck.calc_deviance_resids(model = model,
                                                    df = df,
                                                    y = y)

        # Calculate length (n) of residuals vector
        nResids = Base.length(devResids)

        ###################################################
        # - Calculate hat matrix
        ###################################################

        # Extract model formula
        mod_formula = StatsModels.formula(model)

        # Calculate hat matrix
        reg = LinRegOutliers.createRegressionSetting(mod_formula, df)
        hat = LinRegOutliers.hatmatrix(reg)

        ###################################################
        # - Standardise deviance residuals
        ###################################################

        # Extract deviance and residual df
        modDeviance = GLM.deviance(model)
        modResidualDf = GLM.dof_residual(model)

        # Calculate 's'
        s = Base.sqrt.(modDeviance ./ modResidualDf)

        # Standard residuals by hat matrix
        # - Same as R: resid(mod1, type="deviance")/sqrt(1 - hatvalues(mod1))
        # - Necessary when dispersion parameter not equal to 1
        dHat = LinearAlgebra.Diagonal(hat)
        StdDevResids = devResids ./ (s * (Base.sqrt.(1 .- dHat)) )
        StdDevResids = StdDevResids[:,1]

        # Calculate square root of Std. deviance residuals
        sqrDevResids = sqrt.(complex(StdDevResids))
        sqrDevResids = abs.(sqrDevResids)

        # Define quantiles
        qx = Distributions.quantile.(Distributions.Normal(),
                                        Base.range(0.5,
                                        stop = (nResids .- 0.5),
                                        length = (nResids))
                                    ./ (nResids .+ 1))

         ###################################################
        # - Plot scale-location plot
        ###################################################

        p = Gadfly.plot(
            # Add points layer
            x = pred,
            y = sqrDevResids,
            Gadfly.layer(Gadfly.Geom.smooth(method = :loess,
                                    smoothing = 0.9),
                                    color=[Colors.colorant"red"]),
            Gadfly.layer(Gadfly.Geom.point),
            # Change plot aesthetics
            Gadfly.Scale.x_continuous(format=:plain),
            Gadfly.Guide.title("Scale-Location"),
            Gadfly.Guide.xlabel("Fitted values",
                                orientation=:horizontal),
            Gadfly.Guide.ylabel("√ Std. deviance residuals"))

        # Return the plot object
        return(p)

    end # End the if-else statement for Pearson vs Deviance residual plot
end     # End the entire function


text:
F(t) = hcat(1.0)
G(t) = hcat(1.0)
Σ(t) = hcat(9.0)
Tau(t) = hcat(1.0)
μ0 = vcat(0.0)
Tau0 = hcat(100.0)

Y = [
    5.22896384735402,
    2.059848221254017,
    11.77696919584194,
    5.721925811242492,
    2.647894413351219,
    8.833040488128724,
    10.112474414603064,
    6.208953089526956,
    7.483218024927158,
    1.1450746725213778,
    2.3827281539488214,
    2.5154588929528154,
    7.6938119335857795,
    11.267170019825972,
    11.203900809640977,
    4.641402249798912,
    6.393560809006704,
    5.434277785253692,
    7.754107613915341,
    2.8018519282190573,
    8.844915959029999,
    6.179788030963477,
    8.986422350276268,
    8.129873884597485,
    10.464689898594969,
    7.43846896751153,
    6.195006607054435,
    7.684383628629019,
    6.432819127902091,
    9.784040586681458,
    3.9386281953008275,
    1.0524045623884053,
    5.2287370014479535,
    3.3793567409786296,
    9.109449328047496,
    2.4107833351064665,
    3.1442786924940815,
    8.694980615700167,
    9.34957815303634,
    1.9899759308134595,
    1.7676330607869204,
    6.164419273710474,
    7.378003270378701,
    3.730353213883039,
    8.15997706282842,
    3.912665621137947,
    2.8222087161027987,
    7.89108559137954,
    2.016432756238303,
    -0.7283326731242079,
    5.34191426139306,
    3.7422889485477566,
    6.879539434261266,
    1.697715484841611,
    5.7492808736295835,
    3.4683312150848638,
    7.897465524526299,
    9.358577329317843,
    9.01915967340802,
    5.098410508423983,
    10.161486865442141,
    12.905026331072838,
    11.193444598591357,
    5.284455585283422,
    5.403785641006497,
    2.984985087125228,
    -0.1108114821742765,
    3.8443055205297396,
    1.1335781598255243,
    7.313228382179326,
    1.1291030218984508,
    0.8188795644543472,
    13.044240581728545,
    4.551403331575773,
    5.182374381380807,
    4.064021944452551,
    13.385904690258663,
    2.7383489638189253,
    4.272321073129483,
    1.936449829003382,
    7.579915913135302,
    6.891186573209264,
    5.696714256363617,
    2.9401339639559136,
    4.045792031292042,
    -2.69331630340043,
    4.253305831682493,
    -6.440609695250053,
    1.2224659946193035,
    3.6232261847393623,
    5.212788618742227,
    -0.5932923334934714,
    6.036279124314067,
    6.862343503023183,
    -2.9726591794689665,
    4.290653945351619,
    -2.420184182764869,
    0.6612926367123761,
    0.6312700637171165,
    -0.041609972691161345,
    -0.24839219509891972,
]

Y = hcat(Y)
push_tuple = SMC.push_state_forward(vcat(0.0), G(1), Tau(1), μ0, Tau0)
@test push_tuple[1] ≈ vcat(0.0)
@test push_tuple[2] ≈ Tau(1) + Tau0

ll_tuple = SMC.calc_loglik_t(hcat(0.0), F(1), Σ(1), μ0, Tau0, Y[1, :])
@test size(ll_tuple[1]) == (1, 1)
@test size(ll_tuple[2]) == (1, 1)
@test isa(ll_tuple[3], Float64)

update_tuple = SMC.update_state(F(1), μ0, Tau0, ll_tuple[1], ll_tuple[2])
@test size(update_tuple[1]) == (1, 1)
@test size(update_tuple[2]) == (1, 1)
m, V, loglik_t = SMC.iterate_kalman_filter_mv(
    vcat(0.0),
    F(1),
    vcat(0.0),
    G(1),
    Σ(1),
    Tau(1),
    μ0,
    Tau0,
    @view(Y[1, :])
)
loglik, ms, vs =
    SMC.kalman_filter_mv((t) -> hcat(0.0), F, (t) -> hcat(0.0), G, Σ, Tau, μ0, Tau0, Y)

@assert size(ms) == (1, 101)
@assert size(vs) == (1, 1, 101)
@test loglik ≈ -275.019326809115

@inferred SMC.kalman_filter_mv(
    (t) -> hcat(0.0),
    F,
    (t) -> hcat(0.0),
    G,
    Σ,
    Tau,
    μ0,
    Tau0,
    Y,
)

model_nomutate = SMC.KalmanModel(
    mutate = false,
    F0 = (t) -> hcat(0.0),
    F1 = F,
    G0 = (t) -> hcat(0.0),
    G1 = G,
    Σ = Σ,
    Tau = Tau,
    μ0 = μ0,
    Tau0 = Tau0,
)

@inferred model_nomutate(Y)
ll_nomutate, state_nomutate = model_nomutate(Y)
@test ll_nomutate ≈ -275.019326809115
@test size(state_nomutate.ms) == (1, 101)
@test size(state_nomutate.Vs) == (1, 1, 101)
@test all(state_nomutate.ms ≈ ms)
@test all(state_nomutate.Vs ≈ vs)

# Test is a 3-dimensional concatenation
# Should be 3x the log-likelihood above

using LinearAlgebra
F0(t) = [0.0, 0.0, 0.0]
G0(t) = [0.0, 0.0, 0.0]
F(t) = diagm([1.0, 1.0, 1.0])
G(t) = diagm(([1.0, 1.0, 1.0]))
Σ(t) = 9.0 * diagm(([1.0, 1.0, 1.0]))
Tau(t) = diagm([1.0, 1.0, 1.0])
μ0 = vcat(0.0, 0.0, 0.0)
Tau0 = 100.0 * diagm([1.0, 1.0, 1.0])

y = [
    5.22896384735402,
    2.059848221254017,
    11.77696919584194,
    5.721925811242492,
    2.647894413351219,
    8.833040488128724,
    10.112474414603064,
    6.208953089526956,
    7.483218024927158,
    1.1450746725213778,
    2.3827281539488214,
    2.5154588929528154,
    7.6938119335857795,
    11.267170019825972,
    11.203900809640977,
    4.641402249798912,
    6.393560809006704,
    5.434277785253692,
    7.754107613915341,
    2.8018519282190573,
    8.844915959029999,
    6.179788030963477,
    8.986422350276268,
    8.129873884597485,
    10.464689898594969,
    7.43846896751153,
    6.195006607054435,
    7.684383628629019,
    6.432819127902091,
    9.784040586681458,
    3.9386281953008275,
    1.0524045623884053,
    5.2287370014479535,
    3.3793567409786296,
    9.109449328047496,
    2.4107833351064665,
    3.1442786924940815,
    8.694980615700167,
    9.34957815303634,
    1.9899759308134595,
    1.7676330607869204,
    6.164419273710474,
    7.378003270378701,
    3.730353213883039,
    8.15997706282842,
    3.912665621137947,
    2.8222087161027987,
    7.89108559137954,
    2.016432756238303,
    -0.7283326731242079,
    5.34191426139306,
    3.7422889485477566,
    6.879539434261266,
    1.697715484841611,
    5.7492808736295835,
    3.4683312150848638,
    7.897465524526299,
    9.358577329317843,
    9.01915967340802,
    5.098410508423983,
    10.161486865442141,
    12.905026331072838,
    11.193444598591357,
    5.284455585283422,
    5.403785641006497,
    2.984985087125228,
    -0.1108114821742765,
    3.8443055205297396,
    1.1335781598255243,
    7.313228382179326,
    1.1291030218984508,
    0.8188795644543472,
    13.044240581728545,
    4.551403331575773,
    5.182374381380807,
    4.064021944452551,
    13.385904690258663,
    2.7383489638189253,
    4.272321073129483,
    1.936449829003382,
    7.579915913135302,
    6.891186573209264,
    5.696714256363617,
    2.9401339639559136,
    4.045792031292042,
    -2.69331630340043,
    4.253305831682493,
    -6.440609695250053,
    1.2224659946193035,
    3.6232261847393623,
    5.212788618742227,
    -0.5932923334934714,
    6.036279124314067,
    6.862343503023183,
    -2.9726591794689665,
    4.290653945351619,
    -2.420184182764869,
    0.6612926367123761,
    0.6312700637171165,
    -0.041609972691161345,
    -0.24839219509891972,
];

Y = hcat(y, y, y)
push_tuple = SMC.push_state_forward(G0(1), G(1), Tau(1), μ0, Tau0)
@test push_tuple[1] ≈ vcat(0.0, 0.0, 0.0)
@test push_tuple[2] ≈ Tau(1) + Tau0

ll_tuple = SMC.calc_loglik_t(F0(1), F(1), Σ(1), μ0, Tau0, Y[1, :])
@test size(ll_tuple[1]) == (3,)
@test size(ll_tuple[2]) == (3, 3)
@test isa(ll_tuple[3], Float64)

update_tuple = SMC.update_state(F(1), μ0, Tau0, ll_tuple[1], ll_tuple[2])
@test size(update_tuple[1]) == (3,)
@test size(update_tuple[2]) == (3, 3)

loglik, _, _ = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y)
@test loglik ≈ -275.019326809115 * 3

Tau(t) = [
    1.0 0.5 0.0
    0.5 1.0 0.0
    0.0 0.0 1.0
]

# @test kalman_filter_loglik_mv(vcat(0.0, 0.0, 0.0), F, vcat(0.0, 0.0, 0.0), G, Σ, Tau, μ0, Tau0, Y) ≈ -818.388291668541

## Test filtered distribution
loglik, ms, Vs = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y)
@test loglik ≈ -818.388291668541
@test size(ms) == (3, 101)
@test size(Vs) == (3, 3, 101)

Y_missing = Matrix{Union{Missing,Float64}}(deepcopy(Y));
breaks = [2, 5, 18];
Y_missing[breaks, :] .= missing;
loglik, ms, Vs = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y_missing)

## The missing likelihood should be the same as if we estimated before/after each break
push_forward(m, v, t) = SMC.push_state_forward(G0(t), G(t), Tau(t), m, v)
ll_pt1, ms, Vs = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y[1:1, :])
m_next, V_next = push_forward(push_forward(ms[:, end], Vs[:, :, end], 2)..., 3)

ll_pt2, ms, Vs = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, m_next, V_next, Y[3:4, :])
m_next, V_next = push_forward(push_forward(ms[:, end], Vs[:, :, end], 5)..., 6)

ll_pt3, ms, Vs = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, m_next, V_next, Y[6:17, :])
m_next, V_next = push_forward(push_forward(ms[:, end], Vs[:, :, end], 18)..., 19)

ll_pt4, ms, Vs = SMC.kalman_filter_mv(F0, F, G0, G, Σ, Tau, m_next, V_next, Y[19:end, :])

@test (ll_pt1 + ll_pt2 + ll_pt3 + ll_pt4) ≈ loglik

## Smoother
ms_smoothed, Vs_smoothed = SMC.kalman_smoother_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y_missing)
@test size(ms_smoothed) == (3, 101)
@test size(Vs_smoothed) == (3, 3, 101)
@inferred SMC.kalman_smoother_mv(F0, F, G0, G, Σ, Tau, μ0, Tau0, Y_missing)

model_nomutate = SMC.KalmanModel(
    mutate = false,
    F0 = F0,
    G0 = G0,
    F1 = F,
    G1 = G,
    Σ = Σ,
    Tau = Tau,
    μ0 = μ0,
    Tau0 = Tau0,
)

state_nomutate = SMC.smooth(model_nomutate, Y_missing);
@test all(ms_smoothed ≈ state_nomutate.ms_smoothed)
@test all(Vs_smoothed ≈ state_nomutate.Vs_smoothed)
@inferred SMC.smooth(model_nomutate, Y_missing);

## Draw posterior pass_threshold
using Random: MersenneTwister
Xs_smoothed =
    @inferred SMC.draw_posterior_path(MersenneTwister(23), ms, Vs, G0, G, Σ, Tau, μ0, Tau0)


text:
module FunctionEnvelope

using BSON
using DataFrames
using NamedArrays
using SchumakerSpline
using Plots
using Distributions
using Random

include("1_structs.jl")
export make_function_dataframe
include("2_evaluate_functions.jl")
export evaluate_function, evaluate_function_and_add_to_dd, get_function_value
include("3_crawler.jl")
export initial_evaluations, maximum_possible_value, get_func, top_function_at_point, crawler
export get_intercept_with_constant, get_intercept, get_upper_envelope
include("4_envelope.jl")
export Envelope, evaluate, evaluate_integral, sample_envelope_given_func, sample_envelope
include("5_plotting.jl")
export plot
end