text:
# Type-dependent optimization flags
function fastflags(@nospecialize(T))
    if T <: IntTypes
        s = "nsw"
    elseif T <: UIntTypes
        s = "nuw"
    else#if T <: FloatingTypes
        s = "fast"
        # s = "nnan ninf nsz arcp contract reassoc"
        # s = "nnan ninf nsz arcp contract reassoc"
    end
    return s
end

function suffix(N::Int, @nospecialize(T))
    if T <: Ptr
        if T <: Ptr{<:IntegerTypes}
            t = "p0i"
        else#if T <: Ptr{<:FloatingTypes}
            t = "p0f"
        end
    elseif T <: IntegerTypes
        t = "i"
    else#if T <: FloatingTypes
        t = "f"
    end
    "v$(N)$(t)$(8sizeof(T))"
end

# Type-dependent LLVM constants
function llvmconst(T, val)
    iszero(val) && return "zeroinitializer"
    typ = llvmtype(T)
    "$typ $val"
end
function llvmconst(::Type{Bool}, val)
    Bool(val) || return "zeroinitializer"
    typ = "i1"
    "$typ $(Int(val))"
end
function llvmconst(N::Integer, T, val)
    iszero(val) && return "zeroinitializer"
    typ = llvmtype(T)
    "<" * join(["$typ $val" for i in 1:N], ", ") * ">"
end
function llvmconst(N::Integer, ::Type{Bool}, val)
    Bool(val) || return "zeroinitializer"
    typ = "i1"
    "<" * join(["$typ $(Int(val))" for i in 1:N], ", ") * ">"
end
function llvmtypedconst(T, val)
    typ = llvmtype(T)
    iszero(val) && return "$typ zeroinitializer"
    "$typ $val"
end
function llvmtypedconst(::Type{Bool}, val)
    typ = "i1"
    Bool(val) || return "$typ zeroinitializer"
    "$typ $(Int(val))"
end

# Type-dependent LLVM intrinsics
const LLVM_INS_Integer = Dict{Symbol,String}()
const LLVM_INS_Int = Dict{Symbol,String}()
const LLVM_INS_UInt = Dict{Symbol,String}()
const LLVM_INS = Dict{Symbol,String}()
const LLVM_OVERLOADED_INS = Dict{Symbol,String}()

LLVM_INS_Integer[:+] = "add"
LLVM_INS_Integer[:-] = "sub"
LLVM_INS_Integer[:*] = "mul"
LLVM_INS_Int[:÷] = "sdiv"
LLVM_INS_Int[:div] = "sdiv"
LLVM_INS_Int[:rem] = "srem"
LLVM_INS_Int[:%] = "srem"
LLVM_INS_UInt[:÷] = "udiv"
LLVM_INS_UInt[:div] = "udiv"
LLVM_INS_UInt[:rem] = "urem"
LLVM_INS_UInt[:%] = "urem"

LLVM_INS_Integer[:~] = "xor"
LLVM_INS_Integer[:&] = "and"
LLVM_INS_Integer[:|] = "or"
LLVM_INS_Integer[:⊻] = "xor"

LLVM_INS_Integer[:<<] = "shl"
LLVM_INS_Integer[:>>>] = "lshr"
LLVM_INS_UInt[:>>] = "lshr"
LLVM_INS_Int[:>>] = "ashr"

LLVM_INS_Integer[:(==)] = "icmp eq"
LLVM_INS_Integer[:(!=)] = "icmp ne"
LLVM_INS_Int[:(>)] = "icmp sgt"
LLVM_INS_Int[:(>=)] = "icmp sge"
LLVM_INS_Int[:(<)] = "icmp slt"
LLVM_INS_Int[:(<=)] = "icmp sle"
LLVM_INS_UInt[:(>)] = "icmp ugt"
LLVM_INS_UInt[:(>=)] = "icmp uge"
LLVM_INS_UInt[:(<)] = "icmp ult"
LLVM_INS_UInt[:(<=)] = "icmp ule"

LLVM_INS[:vifelse] = "select"

LLVM_INS[:+] = "fadd"
LLVM_INS[:-] = "fsub"
LLVM_INS[:*] = "fmul"
LLVM_INS[:/] = "fdiv"
LLVM_INS[:inv] = "fdiv"
LLVM_INS[:rem] = "frem"

LLVM_INS[:(==)] = "fcmp fast oeq"
LLVM_INS[:(!=)] = "fcmp fast une"
LLVM_INS[:(>)] = "fcmp fast ogt"
LLVM_INS[:(>=)] = "fcmp fast oge"
LLVM_INS[:(<)] = "fcmp fast olt"
LLVM_INS[:(<=)] = "fcmp fast ole"

LLVM_OVERLOADED_INS[:^] = "@llvm.pow."
LLVM_OVERLOADED_INS[:abs] = "@llvm.fabs."
LLVM_OVERLOADED_INS[:ceil] = "@llvm.ceil."
LLVM_OVERLOADED_INS[:copysign] = "@llvm.copysign."
LLVM_OVERLOADED_INS[:cos] = "@llvm.cos."
LLVM_OVERLOADED_INS[:exp] = "@llvm.exp."
LLVM_OVERLOADED_INS[:exp2] = "@llvm.exp2."
LLVM_OVERLOADED_INS[:floor] = "@llvm.floor."
LLVM_OVERLOADED_INS[:fma] = "@llvm.fma."
LLVM_OVERLOADED_INS[:log] = "@llvm.log."
LLVM_OVERLOADED_INS[:log10] = "@llvm.log10."
LLVM_OVERLOADED_INS[:log2] = "@llvm.log2."
LLVM_OVERLOADED_INS[:max] = "@llvm.maxnum."
LLVM_OVERLOADED_INS[:min] = "@llvm.minnum."
LLVM_OVERLOADED_INS[:muladd] = "@llvm.fmuladd."
LLVM_OVERLOADED_INS[:powi] = "@llvm.powi."
LLVM_OVERLOADED_INS[:round] = "@llvm.rint."
LLVM_OVERLOADED_INS[:sin] = "@llvm.sin."
LLVM_OVERLOADED_INS[:sqrt] =  "@llvm.sqrt."
LLVM_OVERLOADED_INS[:trunc] = "@llvm.trunc."

function llvmins(func::Symbol, N::Int, T)::String
    ins = get(LLVM_OVERLOADED_INS, func, nothing)
    ins === nothing || return ins * suffix(N, T)
    if T <: IntegerTypes
        d = T <: IntTypes ? LLVM_INS_Int : LLVM_INS_UInt
        return get(d, func) do
            LLVM_INS_Integer[func]
        end
    else
        return LLVM_INS[func]
    end
end

# Convert between LLVM scalars, vectors, and arrays

function scalar2vector(vec, siz, typ, sca)
    tempvec = vec * "_temporary"
    String[
        "$tempvec = insertelement <$siz x $typ> undef, $typ $sca, i32 0",
        "$vec = shufflevector <$siz x $typ> $tempvec, <$siz x $typ> undef, <$siz x i32> zeroinitializer"
    ]
end

function array2vector(vec, siz, typ, arr, tmp="$(arr)_av")
    instrs = String[]
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_iter$i"
    for i in 0:siz-1
        push!(instrs, "$(tmp)_elem$i = extractvalue [$siz x $typ] $arr, $i")
        push!(instrs,
            "$(accum(vec,i)) = " *
                "insertelement <$siz x $typ> $(accum(vec,i-1)), " *
                "$typ $(tmp)_elem$i, i32 $i")
    end
    instrs
end

function vector2array(arr, siz, typ, vec, tmp="$(vec)_va")
    instrs = String[]
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_iter$i"
    for i in 0:siz-1
        push!(instrs,
            "$(tmp)_elem$i = extractelement <$siz x $typ> $vec, i32 $i")
        push!(instrs,
            "$(accum(arr,i)) = "*
                "insertvalue [$siz x $typ] $(accum(arr,i-1)), " *
                "$typ $(tmp)_elem$i, $i")
    end
    instrs
end

# TODO: change argument order
function subvector(vec, siz, typ, rvec, rsiz, roff, tmp="$(rvec)_sv")
    instrs = String[]
    accum(nam, i) = i<0 ? "undef" : i==rsiz-1 ? nam : "$(nam)_iter$i"
    @assert 0 <= roff
    @assert roff + rsiz <= siz
    for i in 0:rsiz-1
        push!(instrs,
            "$(tmp)_elem$i = extractelement <$siz x $typ> $vec, i32 $(roff+i)")
        push!(instrs,
            "$(accum(rvec,i)) = " *
                "insertelement <$rsiz x $typ> $(accum(rvec,i-1)), " *
                "$typ $(tmp)_elem$i, i32 $i")
    end
    instrs
end

function extendvector(vec, siz, typ, voff, vsiz, val, rvec, tmp="$(rvec)_ev")
    instrs = String[]
    accum(nam, i) = i<0 ? "undef" : i==siz+vsiz-1 ? nam : "$(nam)_iter$i"
    rsiz = siz + vsiz
    for i in 0:siz-1
        push!(instrs,
            "$(tmp)_elem$i = extractelement <$siz x $typ> $vec, i32 $i")
        push!(instrs,
            "$(accum(rvec,i)) = " *
                "insertelement <$rsiz x $typ> $(accum(rvec,i-1)), " *
                "$typ $(tmp)_elem$i, i32 $i")
    end
    for i in siz:siz+vsiz-1
        push!(instrs,
            "$(accum(rvec,i)) = " *
                "insertelement <$rsiz x $typ> $(accum(rvec,i-1)), $val, i32 $i")
    end
    instrs
end

# Type conversion

@generated function vreinterpret(::Type{Vec{N,R}},
        v1::Vec{N1,T1}) where {N,R,N1,T1}
    if N*sizeof(R) != N1*sizeof(T1)
        throw("N*sizeof(R) == N1*sizeof(T1) is not true; Trying to reinterpret to a size of $N * $(sizeof(R)) from a size of $N1 * $(sizeof(T1))")
    end
    # @assert N*sizeof(R) == N1*sizeof(T1)
    typ1 = llvmtype(T1)
    vtyp1 = "<$N1 x $typ1>"
    typr = R <: Ptr ? llvmtype(Int) : llvmtype(R)
    vtypr = "<$N x $typr>"
    decls = String[]
    instrs = String[]
    push!(instrs, "%res = bitcast $vtyp1 %0 to $vtypr")
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$R}, Tuple{Vec{$N1,$T1}}, v1)
    end
end
@generated function vreinterpret(::Type{Vec{N,R}},
        v1::Vec{N1,UInt128}) where {N,R,N1}
    T1 = UInt128
    @assert N*sizeof(R) == N1*sizeof(T1)
    typ1 = llvmtype(T1)
    vtyp1 = "[$N1 x $typ1]"
    typr = R <: Ptr ? llvmtype(Int) : llvmtype(R)
    vtypr = "<$N x $typr>"
    decls = String[]
    instrs = String[]
    push!(instrs, "%res = bitcast $vtyp1 %0 to $vtypr")
    push!(instrs, "ret $vtypr %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$N,$R}, Tuple{Vec{$N1,$T1}}, v1)
    end
end
@inline function Base.reinterpret(::Type{SVec{N,R}}, v1::SVec{N1,T1}) where {N,R,N1,T1}
    SVec(vreinterpret(Vec{N,R}, extract_data(v1)))
end
@inline function Base.reinterpret(::Type{SVec{N,R}}, v1::SVec{N1,UInt128}) where {N,R,N1}
    SVec(vreinterpret(Vec{N,R}, extract_data(v1)))
end
@inline vrem(v::Vec{W,I}, ::Type{Vec{W,U}}) where {W,I<:Signed,U<:Unsigned} = vrem(vreinterpret(Vec{W,U}, v), Vec{W,U})
@inline vrem(v::Vec{W,U}, ::Type{Vec{W,I}}) where {W,I<:Signed,U<:Unsigned} = vrem(vreinterpret(Vec{W,I}, v), Vec{W,I})

# Doesn't seem to do anything, eg
# assume(N < 10)
# for n ∈ 1:10
#
# end
# will still often be unrolled by factors greater than 10. =(
# @inline function assume(b::Bool)
#     decls = "declare void @llvm.assume(i1 %cond)"
#     instrs = """
#     %b = trunc i8 %0 to i1
#     call void @llvm.assume(i1 %b)
#     ret void
# """
#     Base.llvmcall((decls, instrs), Nothing, Tuple{Bool}, b)
# end

const FASTOPS = Set((:+, :-, :*, :/, :log, :log2, :log10, :exp, :exp2, :exp10, :sqrt, :pow, :sin, :cos))#, :inv, :muladd, :fma

const VECTOR_SYMBOLS = Dict{Symbol,Symbol}(
    :(==) => :visequal,
    :(!=) => :vnot_equal,
    :(<) => :vless,
    :(<=) => :vless_or_equal,
    :(>) => :vgreater,
    :(>=) => :vgreater_or_equal,
    :(<<) => :vleft_bitshift,
    :(>>) => :vright_bitshift,
    :(>>>) => :vuright_bitshift,
    :(&) => :vand,
    :(|) => :vor,
    :(⊻) => :vxor,
    :(+) => :vadd,
    :(-) => :vsub,
    :(*) => :vmul,
    :(/) => :vfdiv,
    :(÷) => :vdiv,
    :(%) => :vrem,
    :div => :vdiv,
    :rem => :vrem,
    :(~) => :vbitwise_not,
    :(!) => :vnot,
    :(^) => :vpow,
    :abs => :vabs,
    :floor => :vfloor,
    :ceil => :vceil,
    :round => :vround,
    :sin => :vsin,
    :cos => :vcos,
    :exp => :vexp,
    :exp2 => :vexp2,
    :exp10 => :vexp10,
    :inv => :vinv,
    :log => :vlog,
    :log10 => :vlog10,
    :log2 => :vlog2,
    :sqrt => :vsqrt,
    :trunc => :vtrunc,
    :sign => :vsign,
    :copysign => :vcopysign,
    :flipsign => :vflipsign,
    :max => :vmax,
    :min => :vmin,
    :fma => :vfma,
    :muladd => :vmuladd,
    :all => :vall,
    :any => :vany,
    :maximum => :vmaximum,
    :minimum => :vminimum,
    :prod => :vprod,
    :sum => :vsum,
    :reduce => :vreduce,
    :isfinite => :visfinite,
    :isinf => :visinf,
    :isnan => :visnan,
    :issubnormal => :vissubnormal,
    :fmadd => :vfmadd,
    :fmsub => :vfmsub,
    :fnmadd => :vfnmadd,
    :fnmsub => :vfnmsub
)


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# DFS iterator test
fst = WFST(Dict(1=>1))
add_arc!(fst, 1, 2, 1, 1, 1.0)
add_arc!(fst, 1, 5, 1, 1, 1.0)
add_arc!(fst, 1, 3, 1, 1, 1.0)
add_arc!(fst, 2, 6, 1, 1, 1.0)
add_arc!(fst, 2, 4, 1, 1, 1.0)
add_arc!(fst, 5, 4, 1, 1, 1.0)
add_arc!(fst, 3, 4, 1, 1, 1.0)
add_arc!(fst, 3, 7, 1, 1, 1.0)
initial!(fst,1)
println(fst)

println("DFS un-folded")
println("Starting from state 1")
dfs = FiniteStateTransducers.DFS(fst,1)
dfs_it = Int[]
completed_seq = Int[]
for (p,s,n,d,e,a) in dfs
  println("$p,$s,$n")
  if d == true
    println("visiting first time $n (Gray)")
    push!(dfs_it,n)
  else
    if e
      println("completed state $s (Black)")
      push!(completed_seq,s)
    else
      println("node $n already visited")
    end
  end
end
@test completed_seq == [6,4,2,5,7,3,1]
@test dfs_it == [2,6,4,5,3,7]

L = Int
fst = WFST(Dict(1=>1))
add_arc!(fst, 1, 2, 1, 1, 1.0)
add_arc!(fst, 2, 3, 1, 1, 1.0)
add_arc!(fst, 3, 1, 1, 1, 1.0)
add_arc!(fst, 3, 4, 1, 1, 1.0)
add_arc!(fst, 4, 5, 1, 1, 1.0)
add_arc!(fst, 5, 6, 1, 1, 1.0)
add_arc!(fst, 6, 7, 1, 1, 1.0)
add_arc!(fst, 7, 8, 1, 1, 1.0)
add_arc!(fst, 8, 5, 1, 1, 1.0)
add_arc!(fst, 2, 9, 1, 1, 1.0)
initial!(fst,1)
final!(fst,8)
println(fst)

println(fst)
scc, c, v = get_scc(fst)
@test scc == [[1, 2, 3], [9], [4], [5, 6, 7, 8]]

# is_visited test
L = Int
ϵ = get_eps(L)
fst = WFST(Dict(i=>i for i=1:4))
add_arc!(fst, 1 => 2, 1=>1,rand() )
add_arc!(fst, 1 => 3, 1=>1, rand())
add_arc!(fst, 2 => 4, 1=>ϵ, rand())
add_arc!(fst, 2 => 3, ϵ=>ϵ, rand())
add_arc!(fst, 2 => 3, 2=>2, rand())
add_arc!(fst, 3 => 4, 3=>3, rand())
add_arc!(fst, 3 => 3, 4=>4, rand())
#print(fst)
@test is_acceptor(fst) == false
## test rm_state!
@test count_eps(get_ilabel,fst) == 1
@test count_eps(get_olabel,fst) == 2
rm_state!(fst,2)
#print(fst)
@test is_acceptor(fst) == true
@test count_eps(get_ilabel,fst) == 0
@test count_eps(get_olabel,fst) == 0
@test size(fst) == (3,3)

# test DFS against recursive
I = O = Int32
fst = WFST(Dict(i=>i for i =1:5),Dict(100*i=>i for i =1:5))
w = randn(5)
add_arc!(fst, 1, 1, 1, 100*1, w[1])
add_arc!(fst, 1, 2, 1, 100*1, w[1])
add_arc!(fst, 1, 3, 2, 100*2, w[2])
add_arc!(fst, 1, 4, 3, 100*3, w[3])
add_arc!(fst, 2, 4, 4, 100*4, w[4])
add_arc!(fst, 3, 4, 5, 100*5, w[5])
initial!(fst, 1)
final!(fst, 4)
# fst have all nodes connected
v = FiniteStateTransducers.recursive_dfs(fst)
@test all(v)
fst2 = deepcopy(fst)
add_arc!(fst2, 8, 10, 5, 100*5, w[5])
add_arc!(fst2, 8, 9, 5, 100*5, w[5])
add_arc!(fst2, 8, 8, 5, 100*5, w[5])
# fst2 5 < state <= 10 are not connected with previous ones
add_arc!(fst2, 3, 11, 5, 100*5, w[5])
add_arc!(fst2, 11, 12, 5, 100*5, w[5])
add_arc!(fst2, 12, 12, 5, 100*5, w[5])
final!(fst2,12,w[1])
# fst2 state >= 11 are connected with first nodes
v = FiniteStateTransducers.recursive_dfs(fst2)
@test all(v[1:4])
@test all(.!v[5:10])
@test all(v[11:end])

# testing lazy version
v = is_visited(fst)
@test all(v)
v = is_visited(fst2)
@test all(v[1:4])
@test all(.!v[5:10])
@test all(v[11:end])
v = is_visited(fst2,8)
@test all(v[8:10])
@test all(.!v[1:7])
@test all(.!v[11:end])

# connect tests
#println(fst2)
## testing connect
fst3 = connect(fst2)
@test all(is_visited(fst3))
@test size(fst3) == (6,9)
@test isfinal(fst3,6)
@test isfinal(fst3,4)
@test isinitial(fst3,1)
@test get(get_weight(fst3,6)) ≈ w[1]

ext: jl


text:
facts("Inferring Factor structure") do

    context("All indices outer") do
        a[i, j] ~ Normal(rand(5, 5), rand(5, 5))
        b[j, k] ~ Gamma(rand(5, 3), rand(5, 3))

        nodes = Node[a, b]
        fi = get_structure(nodes...)

        # order should not be important, so compare Sets
        @fact Set(fi.indices) --> Set([:i, :j, :k])
        @fact Set(zip(fi.indices, fi.maxvals)) --> Set([(:i, 5), (:j, 5), (:k, 3)])

        # check definition of inds_in_factor
        @fact Set(fi.indices[fi.inds_in_factor[:a]]) --> Set([:i, :j])
        @fact Set(fi.indices[fi.inds_in_factor[:b]]) --> Set([:k, :j])

        # check definition of inds_in_node
        @fact Set(a.outerinds[fi.inds_in_node[:a]]) --> Set([:i, :j])
        @fact Set(b.outerinds[fi.inds_in_node[:b]]) --> Set([:j, :k])
    end

    context("Some indices inner") do
        dims = (5, 3)
        m = [rand(dims[1]) for x in 1:dims[2]]
        VV = [diagm(rand(dims[1])) for x in 1:dims[2]]

        d[i, j] ~ MvNormal(m, VV)
        a[i, k] ~ Normal(rand(5, 4), rand(5, 4))

        nodes = Node[a, d]
        fi = get_structure(nodes...)

        # order should not be important, so compare Sets
        @fact Set(fi.indices) --> Set([:j, :k])
        @fact Set(zip(fi.indices, fi.maxvals)) --> Set([(:j, 3), (:k, 4)])

        # check definition of inds_in_factor
        @fact Set(fi.indices[fi.inds_in_factor[:a]]) --> Set([:k])
        @fact Set(fi.indices[fi.inds_in_factor[:d]]) --> Set([:j])

        # check definition of inds_in_node
        @fact Set(a.outerinds[fi.inds_in_node[:a]]) --> Set([:k])
        @fact Set(d.outerinds[fi.inds_in_node[:d]]) --> Set([:j])
    end

    context("Dimension mismatch throws error") do
        a[i, j] ~ Normal(rand(5, 5), rand(5, 5))
        b[j, k] ~ Gamma(rand(4, 3), rand(4, 3))

        nodes = Node[a, b]

        @fact_throws ErrorException get_structure(nodes...)
    end

    context("Inner indices optional") do
        dims = (5, 3)
        m = [rand(dims[1]) for x in 1:dims[2]]
        VV = [diagm(rand(dims[1])) for x in 1:dims[2]]
        LL = rand(3)

        d[j] ~ MvNormal(m, VV)
        L[j] ~ Const(LL)

        nodes = Node[L, d]
        fi = get_structure(nodes...)
        inds = fi.indices
        maxvals = fi.maxvals

        @fact Set(inds) --> Set([:j])
    end

    context("Inner indices arbitrary") do
        dims = (5, 3)
        m = [rand(dims[1]) for x in 1:dims[2]]
        VV = [diagm(rand(dims[1])) for x in 1:dims[2]]
        LL = rand(5, 5, 3)

        d[a, j] ~ MvNormal(m, VV)
        L[a, a, j] ~ Const(LL)

        nodes = Node[L, d]
        fi = get_structure(nodes...)
        inds = fi.indices
        maxvals = fi.maxvals

        @fact Set(inds) --> Set([:j])
    end

    context("Some inner indices forced") do
        a[i, j] ~ Normal(rand(5, 5), rand(5, 5))
        b[j, k] ~ Gamma(rand(5, 3), rand(5, 3))

        nodes = Node[a, b]
        fi = get_structure([:j], nodes...)

        # order should not be important, so compare Sets
        # index j should be missing
        @fact Set(fi.indices) --> Set([:i, :k])
    end

end

facts("Basic factor construction") do

    context("Simple univariate nodes") do
        dims = (10, 2)

        a[j] ~ Normal(rand(dims[2]), ones(dims[2]))
        b ~ Gamma(1, 1)
        y[i, j] ~ Normal(rand(dims), ones(dims))

        f = @factor LogNormalFactor y a b
        inds = f.inds.indices
        maxvals = f.inds.maxvals

        @fact Set(inds) --> Set([:i, :j, :scalar])

        # internals use names defined for factor type, not node names
        @fact project_inds(:x, f, inds) --> [:i, :j]
        @fact project_inds(:μ, f, inds) --> [:j]
        @fact project_inds(:τ, f, inds) --> [:scalar]
        @fact isa(project(:x, f, maxvals), Normal) --> true
        @fact isa(project(:μ, f, maxvals), Normal) --> true
        @fact isa(project(:τ, f, maxvals), Gamma) --> true
        @fact value(f) --> isfinite
    end

    context("Multivariate nodes in factor") do
        dims = (5, 3)
        m = [rand(dims[1]) for x in 1:dims[2]]
        VV = [diagm(rand(dims[1])) for x in 1:dims[2]]

        y[i, j] ~ MvNormal(m, VV)
        μ[i] ~ Normal(zeros(dims[1]), ones(dims[1]))
        Λ[j] ~ Const([diagm(rand(dims[1])) for x in 1:dims[2]])

        f = @factor LogMvNormalCanonFactor y μ Λ
        inds = f.inds.indices
        maxvals = f.inds.maxvals

        @fact Set(inds) --> Set([:j])
        @fact project_inds(:x, f, inds) --> [:j]
        @fact project_inds(:μ, f, inds) --> [Colon()]
        @fact project_inds(:Λ, f, inds) --> [:j]
        @fact isa(project(:x, f, maxvals), MvNormal) --> true
        @fact isa(project(:μ, f, maxvals), Vector{Normal}) --> true
        @fact isa(project(:Λ, f, maxvals), Matrix{Float64}) --> true
        @fact value(f) --> isfinite
    end

    context("pmodel macro") do
        dims = (10, 2)

        a[j] ~ Normal(rand(dims[2]), ones(dims[2]))
        b ~ Gamma(1, 1)
        y[i, j] ~ Normal(rand(dims), ones(dims))

        f1 = @factor LogNormalFactor y a b
        @pmodel begin
            y ~ Normal(a, b)
        end
        f2 = pmodel_factors[1]

        @fact isa(pmodel_factors, Array{Factor, 1}) --> true
        @fact isa(f2, typeof(f1)) --> true
        for f in fieldnames(f1.inds)
            @fact getfield(f2.inds, f) --> getfield(f1.inds, f)
        end
        @fact f2.namemap --> f1.namemap
    end

end

facts("Unrolling and rerolling parameters") do
    context("Distributions") do
        dists = [Gamma(1.1, 1.), Normal(),
                 MvNormalCanon(ones(5), eye(5)),
                 Wishart(5, eye(2)),]

        for d in dists
            par_sizes = get_par_sizes(d)
            npars = mapreduce(prod, +, par_sizes)
            x = unroll_pars(d)
            dd = reroll_pars(d, par_sizes, x)
            @fact isa(x, Vector) --> true
            @fact length(x) --> npars
            @fact dd --> d
        end
    end
    context("Nodes") do
        p, d = 5, 7
        x[i] ~ Gamma(1.1 * ones(p), ones(p))
        y[i] ~ MvNormalCanon([ones(d) for p in 1:p], [eye(d) for p in 1:p])

        node_list = [x, y]

        for n in node_list
            # get parameter sizes for single distribution
            par_sizes = get_par_sizes(n.data[1])
            npars = mapreduce(prod, +, par_sizes)

            # make a copy of the original node
            n_bak = deepcopy(n)

            # unroll parameters
            v = unroll_pars(n)
            @fact length(v) --> prod(size(n)) * npars
            @fact isa(v, Vector) --> true
            update_pars!(n, v)
            @fact n --> n_bak
        end
    end
end


text:
@testset "Operators" begin
    @testset "intersect" begin
        g = rand_graph(10, 20, graph_type=GRAPH_T)
        @test intersect(g, g).num_edges == 20
    end
end


text:
struct UnionFind
    par::Vector{Int}
    size::Vector{Int}
    UnionFind(n) = new(collect(1:n), ones(Int, n))
end

root!(uf, x) = uf.par[x] == x ? x : uf.par[x] = root!(uf, uf.par[x])

same(uf, x, y) = root!(uf, x) == root!(uf, y)

function unite!(uf, x, y)
    x = root!(uf, x)
    y = root!(uf, y)
    x == y && return false

    if uf.size[x] > uf.size[y]
        uf.par[y] = x
        uf.size[x] += uf.size[y]
    else
        uf.par[x] = y
        uf.size[y] += uf.size[x]
    end
    return true
end
size!(uf, x) = uf.size[root!(uf, x)]

struct Edge
    s::Int
    t::Int
    c::Int
end

import Base: isless
isless(e₁::Edge, e₂::Edge) = isless(e₁.c, e₂.c) # For sorting

function solve(n, m)
    edges = [Edge((parse.(Int, split(readline())))...) for _ = 1:m]

    sort!(edges)
    uf = UnionFind(n)
    cnt_edges = 0
    ans = 0
    for e ∈ edges
        if !same(uf, e.s, e.t)
            unite!(uf, e.s, e.t)
            cnt_edges += 1
        end

        if cnt_edges == (n + 1) ÷ 2
            ans = e.c
            break
        end
    end

    println(ans)
end

function main()
    while true
        n, m = parse.(Int, split(readline()))
        if n == 0 && m == 0
            break
        else
            solve(n, m)
        end
    end
end

main()


text:
export tukey
"""
    tukey(N::Int, r::AbstractArray) -> AbstractArray

Tukey window is a rectangular window with the first and
last r/2 percent of the samples equal to parts of a cosine.
"""
function tukey(N::Int, r::AbstractFloat=0.5)
    n = (0:N-1)/(N-1)
    w = zeros(N)
    for i = 1:N
        x = n[i]
        if 0 <= x <= r/2
            w[i] = 0.5 + 0.5*cos(2*pi/r*(x-r/2))
        elseif r/2 <= x <= 1-r/2
            w[i] = 1.0
        else
            w[i] = 0.5 + 0.5*cos(2*pi/r*(x-1+r/2))
        end
    end
    return w
end

export gauss
"""
    gauss(N::Int, a::AbstractFloat=0.5)

a is proportional to the standard deviation
"""
function gauss(N::Int, a::AbstractFloat=0.4)
    n = 0:N-1
    m = (N-1)/2
    return @. exp(-((n - m)/m/a)^2/2)
end


text:
using LinearAlgebra
using JuMP
using Cbc

model = Model(Cbc.Optimizer)
set_optimizer_attribute(model, "cuts","off")
# set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ON)
# set_optimizer_attribute(model, "maxNodes", 5)
month_num = 8
month = 1:month_num
# paramerters
q = [5000] # product fix cost
p = [100] # cost per product
d = diagm([400, 400, 800, 800, 1200, 1200, 1200, 1200]) # demand of product by month
s_init = [200,0,0,0, 0,0,0,0] # initial stock
h = [5] # stock cost per product month

@variable(model, x[i in month, t in month] >= 0)
@variable(model, s[i in month, t in month] >= 0)
@variable(model, y[i in month], Bin)

# 库存平衡
@constraint(model ,dem_sat1[t in month], s_init[t] + x[1,t] == d[1, t] + s[1, t])
@constraint(model, dem_sat[i in month[2:end], t in month], s[i - 1, t] + x[i, t] == d[i, t] + s[i, t] )

# 最大产量约束
@constraint(model, vub1[i in month], x[i, 1] <= (d[1,1] - s_init[1]) * y[i])
@constraint(model, vub[i in month, t in month[2:end]], x[i, t] <= d[t,t] * y[i])
# 目标函数
@expression(model ,cost, sum(p[1] * x[i,t] + h[1] * s[i, t] for i in month, t in month[i:end]) + sum(q[1] * y[i] for i in month))
@objective(model, Min, cost)
for t in month
    fix(s[t,t], 0;force=true)
end
optimize!(model)

# 显示结果
cost_val = value.(cost)
x_val = value.(x)
s_val = value.(s)
@show cost_val;
@show x_val;
@show s_val;

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
using Sockets
using Serialization
using REPL
using Logging

mutable struct ServerSideSession
    display_properties::Dict
    in_module::Module
end

function send_header(io, ser_version=Serialization.ser_version)
    write(io, PROTOCOL_MAGIC, PROTOCOL_VERSION)
    write(io, UInt32(ser_version))
    flush(io)
end

# Like `sprint()`, but uses IOContext properties `ctx_properties`
#
# This is used to stringify results before sending to the client. This is
# beneficial because:
#   * It allows us to show the user-defined types which exist only on the
#     remote server
#   * It allows us to limit the amount of data transmitted (eg, large arrays
#     are truncated when the :limit IOContext property is set)
function sprint_ctx(f, session)
    io = IOBuffer()
    ctx = IOContext(io, :module=>session.in_module, session.display_properties...)
    f(ctx)
    String(take!(io))
end

# Server-side expression processing
function preprocess_expression!(ex, new_stdout)
    if ex isa Expr
        if ex.head ∉ (:quote, :inert)
            map!(e->preprocess_expression!(e, new_stdout), ex.args, ex.args)
        end
    elseif ex === STDOUT_PLACEHOLDER
        ex = new_stdout
    end
    return ex
end

function eval_message(session, messageid, messagebody)
    try
        if messageid in (:eval, :eval_and_get)
            result = nothing
            resultstr = sprint_ctx(session) do io
                with_logger(ConsoleLogger(io)) do
                    expr = preprocess_expression!(messagebody, io)
                    result = Base.eval(session.in_module, expr)
                    if messageid === :eval && !isnothing(result)
                        # We require invokelatest here in case the user
                        # modifies any method tables after starting the session,
                        # which change methods of `show`
                        Base.invokelatest(show, io, MIME"text/plain"(), result)
                    end
                end
            end
            return messageid === :eval ?
                (:eval_result, resultstr) :
                (:eval_and_get_result, (result, resultstr))
        elseif messageid === :help
            resultstr = sprint_ctx(session) do io
                md = Main.eval(REPL.helpmode(io, messagebody))
                Base.invokelatest(show, io, MIME"text/plain"(), md)
            end
            return (:help_result, resultstr)
        elseif messageid === :display_properties
            session.display_properties = messagebody::Dict
            return nothing
        elseif messageid === :in_module
            mod = Main.eval(messagebody)::Module
            session.in_module = mod
            resultstr = "Evaluating commands in module $mod"
            return (:in_module, resultstr)
        elseif messageid === :repl_completion
            # See REPL.jl complete_line(c::REPLCompletionProvider, s::PromptState)
            partial, full = messagebody
            ret, range, should_complete = REPL.completions(full, lastindex(partial),
                                                           session.in_module)
            result = (unique!(map(REPL.completion_text, ret)),
                      partial[range], should_complete)
            return (:completion_result, result)
        else
            return (:error, "Unknown message id: $messageid")
        end
    catch _
        resultstr = sprint_ctx(session) do io
            Base.invokelatest(Base.display_error, io, Base.catch_stack())
        end
        return (:error, resultstr)
    end
end

function evaluate_requests(request_chan, response_chan)
    session = ServerSideSession(Dict(), Main)

    while true
        try
            request = take!(request_chan)
            result = eval_message(session, request...)
            if !isnothing(result)
                put!(response_chan, result)
            end
        catch exc
            if exc isa InvalidStateException && !isopen(request_chan)
                break
            elseif exc isa InterruptException
                # Ignore any interrupts which are sent while we're not
                # evaluating a command.
                continue
            else
                rethrow()
            end
        end
    end
end

function deserialize_requests(socket, repl_backend, request_chan, response_chan)
    while isopen(socket)
        request = nothing
        try
            request = deserialize(socket)
        catch exc
            resultstr = sprint() do io
                if exc isa UndefVarError
                    Base.display_error(io, exc, nothing)

                    print(io, """
                        This can happen when you try to pass a custom type
                        from the client which doesn't exist on the server""")
                else
                    Base.display_error(io, Base.catch_stack())
                    println(io, """
                        Unexpected error deserializing RemoteREPL message""")
                end
            end
            put!(response_chan, (:error, resultstr))
            continue
        end
        @debug "Client command" request
        if request isa Tuple && length(request) == 2 && request[1] isa Symbol
            messageid, messagebody = request
            # Handle flow control messages in the RemoteREPL frontend, here.
            if messageid === :exit
                break
            elseif messageid == :interrupt
                # Soft interrupt - this will only work if the
                # `repl_backend` task yields.  It won't work if
                # `repl_backend` is executing a compute-heavy task
                # which doesn't yield to the scheduler.
                schedule(repl_backend, InterruptException(); error=true)
            else
                # All other messages are handled by the evaluator
                put!(request_chan, request)
            end
        else
            put!(response_chan, (:error, "Invalid message of type: $(typeof(request))"))
        end
    end
end

function serialize_responses(socket, response_chan)
    try
        while true
            response = take!(response_chan)
            serialize(socket, response)
        end
    catch exc
        if isopen(response_chan) && isopen(socket)
            rethrow()
        end
    end
end

# Serve a remote REPL session to a single client over `socket`.
function serve_repl_session(socket)
    send_header(socket)
    @sync begin
        request_chan = Channel(1)
        response_chan = Channel(1)

        repl_backend = @async try
            evaluate_requests(request_chan, response_chan)
        catch exc
            @error "RemoteREPL backend crashed" exception=exc,catch_backtrace()
        finally
            close(response_chan)
        end

        @async try
            serialize_responses(socket, response_chan)
        catch exc
            @error "RemoteREPL responder crashed" exception=exc,catch_backtrace()
        finally
            close(socket)
        end

        try
            deserialize_requests(socket, repl_backend, request_chan, response_chan)
        catch exc
            @error "RemoteREPL frontend crashed" exception=exc,catch_backtrace()
            rethrow()
        finally
            close(socket)
            close(request_chan)
        end
    end
end

"""
    serve_repl([address=Sockets.localhost,] port=$DEFAULT_PORT)
    serve_repl(server)

Start a REPL server listening on interface `address` and `port`. In normal
operation `serve_repl()` serves REPL clients indefinitely (ie., it does not
return), so you will generally want to launch it using `@async serve_repl()` to
do other useful work at the same time.

If you want to be able to stop the server you can pass an already-listening
`server` object (the result of `Sockets.listen()`). The server can then be
cancelled from another task using `close(server)` as necessary to control the
server lifetime.

## Security

`serve_repl()` uses an *unauthenticated, unecrypted protocol* so it should not
be used on open networks or multi-user machines where other users aren't
trusted. For open networks, use the default `address=Sockets.localhost` and the
automatic ssh tunnel support provided by the client-side `connect_repl()`.
"""
function serve_repl(address=Sockets.localhost, port::Integer=DEFAULT_PORT)
    server = listen(address, port)
    try
        serve_repl(server)
    finally
        close(server)
    end
end
serve_repl(port::Integer) = serve_repl(Sockets.localhost, port)

function serve_repl(server::Base.IOServer)
    open_sockets = Set()
    @sync try
        while isopen(server)
            socket = accept(server)
            push!(open_sockets, socket)
            peer=getpeername(socket)
            @async try
                serve_repl_session(socket)
            catch exc
                if !(exc isa EOFError && !isopen(socket))
                    @warn "Something went wrong evaluating client command" #=
                        =# exception=exc,catch_backtrace()
                end
            finally
                @info "REPL client exited" peer
                close(socket)
                pop!(open_sockets, socket)
            end
            @info "REPL client opened a connection" peer
        end
    catch exc
        if exc isa Base.IOError && !isopen(server)
            # Ok - server was closed
            return
        end
        @error "Unexpected server failure" isopen(server) exception=exc,catch_backtrace()
        rethrow()
    finally
        for socket in open_sockets
            close(socket)
        end
    end
end


text:
# ------------------------------------------------------------------------------------------
# # Factorizaciones y otras diversiones
# Autor: Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)
# (con edición de Jane Herriman=
#
# ## Sinopsis
#  - Factorizaciones
#  - Estructuras de matrices especiales
#  - Álgebra lineal genérica
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Empezamos con un sistema lineal de la forma
#
# `Ax = b`
# ------------------------------------------------------------------------------------------

A = randn(3,3)

x = fill(1.0, (3))
b = A*x

# ------------------------------------------------------------------------------------------
# ### Factorización
# La función `\` esconde cómo regularmente se resuelve el problema.
#
# Dependiendo de las dimensiones de `A`, distintos métodos son elegidos para resolver el
# problema.
#
# ```
# Ax = b
# ```
#
# Un paso intermedio en la solución es el cálculo de la factorización de la matriz `A`.
#
# Básicamente, una factorización de `A` es una manera de expresar `A` como el producto de
# matrices triangulares, unitarias, y de permutación.
#
# Julia guarda estas factorizaciones usando un tipo abstracto llamado `Factorization` y
# varios subtipos.
#
# Un objeto `Factorization` entonces debería ser pensado como una represenatación de la
# matriz `A`.
#
# #### LU
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Cuando `A` es cuadrada, el sistema linear es resuelto factorizando `A` vía
#
# ```
# PA=LU
# ```
#
# donde `P` es una matriz de permutación, `L` es un triangular inferior unitaria y `U` es
# triangular superior.
#
# Julia permita calcular la facorización LU y define un tipo de factorización compuesta para
# guardarlo.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Podemos hacer una factorización LU sobre `A` vía `lu(A)` ó `lufact(A)`.
#
# La función `lu` regresa matrices `l` y `u` y un vector de permutación `p`.
# ------------------------------------------------------------------------------------------

l,u,p = lu(A)

# ------------------------------------------------------------------------------------------
# El pivoteo está prendido por default, o sea que no podemos assumir que A == LU.
# Vamos a comprobar esto viendo la norma de `LU - A`:
# ------------------------------------------------------------------------------------------

norm(l*u - A)

# ------------------------------------------------------------------------------------------
# Esto nos muestra que queremos tomar en cuenta el pivoteo!
#
# Podemos pensar en `A[p,:]` como la sintaxis para `PA`, o como el producto de una matriz de
# permutación con `A`
# ------------------------------------------------------------------------------------------

norm(l*u - A[p,:])

# ------------------------------------------------------------------------------------------
# Por otro lado, podemos apagar el pivoteo usando `Val{false}` en Julia 0.6 ó `Val(false)`
# en versiones más modernas.
# ------------------------------------------------------------------------------------------

l,u,p = lu(A, Val{false})

# ------------------------------------------------------------------------------------------
# Cuando apagamos el pivoteo, `LU = A`
# ------------------------------------------------------------------------------------------

norm(l*u - A)

# ------------------------------------------------------------------------------------------
# Una segunda manera de hacer la factorización es con `lufact`.
# ------------------------------------------------------------------------------------------

Alu = lufact(A)

# ------------------------------------------------------------------------------------------
# Distintas partes de la factorización las puedes accesar con índices especiales
# ------------------------------------------------------------------------------------------

Alu[:P]

Alu[:L]

Alu[:U]

# ------------------------------------------------------------------------------------------
# Podemos calcular la solución de $Ax=b$ del objecto de factorización
# ------------------------------------------------------------------------------------------

# PA = LU
# A = P'LU
# P'LUx = b
# LUx = Pb
# Ux = L\Pb
# x = U\L\Pb
Alu[:U]\(Alu[:L]\(Alu[:P]b))

# ------------------------------------------------------------------------------------------
# *Más importantemente,* podemos despachar sobre el tipo de `LU` y simplemente resolver el
# problema por medio de
# ------------------------------------------------------------------------------------------

Alu\b

# ------------------------------------------------------------------------------------------
# Esto puede ser útil si el mismo lado izquierdo es usado para lados derechos.
#
# La factorización también se puede usar para calcular el determinante pues
#
# $\det(A)=\det(PLU)=\det(P)\det(U)=\pm \prod u_{ii}$
#
# porque $U$ es triangular y el signo está definido por $\det(P)$.
# ------------------------------------------------------------------------------------------

det(A)

det(Alu)

# ------------------------------------------------------------------------------------------
# #### QR
# Cuando `A` es alta,
# ------------------------------------------------------------------------------------------

Atall = randn(3, 2)

# ------------------------------------------------------------------------------------------
# Julia calcula la solución de mínimos cuadrados $\hat{x}$ que minimiza $\|Ax-b\|_2$.
#
# Esto se puede hacer factorizando
#
# ```
# A=QR
# ```
#
# donde $Q$ es unitaria/ortogonal y
#
# $R=\left(\begin{smallmatrix}R_0\\0\end{smallmatrix}\right)$ y  $R_0$ es triangular
# superior.
#
# con la factorización QR la norma mínima se puede expresar como
#
# \begin{equation*}
# \|Ax-b\|=\|QRx-b\|=\|Q(Rx-Q'b)\|=\|Rx-Q'b\|=\left\|\begin{pmatrix}R_0x-Q_0'b\\Q_1'b\end{pm
# atrix}\right\|=\|R_0x-Q_0'b\|+\|Q_1'b\|
# \end{equation*}
#
# Y entonces el problema se puede reducir a resolver el problema cuadrado $R_0x=Q_0'b$ para
# $x$.
#
# Podemos factorizar QR sobre `Atall` vía
# ------------------------------------------------------------------------------------------

Aqr = qrfact(Atall)

# ------------------------------------------------------------------------------------------
# Otra característica de la factorización QR es que los tipos `Q` para guardar las matrices
# unitarias $Q$. Se pueden extraer de tipos `QR` con los índices
# ------------------------------------------------------------------------------------------

Aqr[:Q]

# ------------------------------------------------------------------------------------------
# Similarmente, la matriz superior triangular $R$ se puede extraer con el índice
# ------------------------------------------------------------------------------------------

Aqr[:R]

# ------------------------------------------------------------------------------------------
# En este caso se guarda R como una matriz 2x2 en vez de 3x2 porque el último renglón de R
# está lleno de 0's.

#
#
# Aunque la matriz `Aqr[:Q]` se imprime como $3\times 3$ en el objeto de factorización, en
# la práctica puede representar la versión delgada también. Así
# ------------------------------------------------------------------------------------------

Aqr[:Q]*ones(2)

# ------------------------------------------------------------------------------------------
# funciona y representa a $3 x 2$ matrix por un vector de 2-elementos.
#
# Similarmente,
# ------------------------------------------------------------------------------------------

Aqr[:Q]*ones(3)

# ------------------------------------------------------------------------------------------
# funciona representando la matriz $3x3$ y un vector de 3 elementos.
#
# Sin embargo, esto no significa que podemos multiplicar a `Q` por vectores de longitued
# arbitraria.
# ------------------------------------------------------------------------------------------

Aqr[:Q]*ones(4)

# ------------------------------------------------------------------------------------------
# La matriz tiene representación interna compacta, entonces indexar sólo hace sentido si uno
# sabe cómo la factorización guarda datos.
# ------------------------------------------------------------------------------------------

Aqr[:Q][1]

# ------------------------------------------------------------------------------------------
# El objeto QRCompactWY `\` tiene un método para QR y entonces el problema de los mínimos
# cuadrados es resuelto con
# ------------------------------------------------------------------------------------------

Aqr\b

# ------------------------------------------------------------------------------------------
# Y si en vez escribimos
# ------------------------------------------------------------------------------------------

Atall\b

# ------------------------------------------------------------------------------------------
# En vez de factorizar con QR a `Atall` primero, Julia va a defaultear factorizar *con*
# pivoteo.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Este default a pivotear la factorización QR le permite a Julia resolver problemas
# deficientes de rango.
#
# Podemos explícitamente escoger pivotear durante la factorización QR (de una matriz
# singular, por ejemplo) con `Val{true}`.
# ------------------------------------------------------------------------------------------

v = randn(3)
# Tomar el producto exterio de un vector consigo mismo nos da una matriz singular
singmatriz = v * v'

Aqrp = qrfact(singmatriz,Val{true})

# ------------------------------------------------------------------------------------------
# Notamos que el tipo que resulta del objecto de Factorization es distinto que antes.
#
# `\` también tiene un método de `QRPivoted` y el problema con rango deficiente es entonces
# calculado como
# ------------------------------------------------------------------------------------------

Aqrp\b

# ------------------------------------------------------------------------------------------
# #### Eigendescompisición y los SVDs (Valores de descomposición Singular)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Los resultados de eigendescomposición y de la descomposición singular de valores se
# guardan en los tipos`Factorization`. Esto también incluye la factorización de Hessenberg y
# de Schur
#
# La eigendescomposición puede ser calculada como
# ------------------------------------------------------------------------------------------

Asym = A + A'
AsymEig = eigfact(Asym)

# ------------------------------------------------------------------------------------------
# Los valores y vectores se pueden recoger del tipo Eigen con un índice especial
# ------------------------------------------------------------------------------------------

AsymEig[:values]

AsymEig[:vectors]

# ------------------------------------------------------------------------------------------
# Una vez más, como la descomposición se guarda en un tipo, podemos despatchar sobre esos
# tipos y explotar un método especializado para cada factorización, e.g. que
# $A^{-1}=(V\Lambda V^{-1})^{-1}=V\Lambda^{-1}V^{-1}$.
# ------------------------------------------------------------------------------------------

inv(AsymEig)*Asym

# ------------------------------------------------------------------------------------------
# Julia también tiene una función `eig` que regresa una tupla con los valores y vectores
# ------------------------------------------------------------------------------------------

eig(Asym)

# ------------------------------------------------------------------------------------------
# No recomendamos esta versión.
#
# La función `svdfact` calcula la descomposición singular de valores
# ------------------------------------------------------------------------------------------

Asvd = svdfact(A[:,1:2])

# ------------------------------------------------------------------------------------------
# y de nuevo `\` tiene un método para el tipo que permite los mínimos cuadrados por SVD
# ------------------------------------------------------------------------------------------

Asvd\b

# ------------------------------------------------------------------------------------------
# Existen funciones especiales para proporcionar los valores sólamente: `eigvals` and
# `svdvals`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Estructuras especiales de matrices
#
# La estructura de matrices es muy importante en el álgebra linear. Ésta estructura se le
# puede hacer explícita a Julia por medio de los tipos compuestos. Ejemplos: `Diagonal`,
# `Triangular`, `Symmetric`, `Hermitian`, `Tridiagonal` y `SymTridiagonal`. Se han escrito
# métodos especializados para cada tipo de matriz especial para aprovechar su estructura.
# Siguen algunos ejemplos:
# ------------------------------------------------------------------------------------------

A

# ------------------------------------------------------------------------------------------
# Creando una matriz Diagonal
# ------------------------------------------------------------------------------------------

Diagonal(diag(A))

Diagonal(A)

# ------------------------------------------------------------------------------------------
# Creando una matriz triangular inferior
# ------------------------------------------------------------------------------------------

LowerTriangular(tril(A))

LowerTriangular(A)

# ------------------------------------------------------------------------------------------
# Creando una matriz simétrica
# ------------------------------------------------------------------------------------------

Symmetric(Asym)

SymTridiagonal(diag(Asym),diag(Asym,1))

# ------------------------------------------------------------------------------------------
# Cuando se sabe que una matriz es e.g. triangular o simétrica Julia puede que resuelva el
# problema más rápido convirtiendo a una matriz especial.
#
# Para algunos procedimientos, Julia checa si el input de matriz es triangular o simétrica y
# lo convierte a tal estructura si lo detecta.
#
# Debería notarse que `Symmetric`, `Hermitian` y `Triangular` no copian la matriz original.
#
# #### Eigenproblema simétrico
# La capacidad de Julia para poder detectar si una matriz es simétrica/Hermitian puede
# influenciar muchísimo sobre qué tan rápido se puede resolver un problema de eigenvalor.
# ------------------------------------------------------------------------------------------

n = 1000;
A = randn(n,n);

# ------------------------------------------------------------------------------------------
# Usamos `A` para genera una matriz simétrica `Asym`
# ------------------------------------------------------------------------------------------

Asym = A + A';

# ------------------------------------------------------------------------------------------
# Ahora creemos una matriz Asym para simular una matriz simétrica con errores de punto
# flotante
# ------------------------------------------------------------------------------------------

Asym_noisy = copy(Asym); Asym_noisy[1,2] += 5eps();

# ------------------------------------------------------------------------------------------
# ¿Puede Julia determinar que ambas `Asym` y `Asym_noisy` son matrices simétricas?
# ------------------------------------------------------------------------------------------

println("Is Asym symmetric? ", issymmetric(Asym))
println("Is Asym_noisy symmetric? ", issymmetric(Asym_noisy))

# ------------------------------------------------------------------------------------------
# Ahora veamos como el ruido de `Asym_noisy` impacta el tiempo en llevar a cabo una
# eigendescomposición
# ------------------------------------------------------------------------------------------

@time eigvals(Asym);

@time eigvals(Asym_noisy);

# ------------------------------------------------------------------------------------------
# Por suerte, le podeemos proveer información explícita sobre la estructura de la matriz a
# Julia
# En este ejemplo, usamos la palabra clave `Symmetric`
# ------------------------------------------------------------------------------------------

@time eigvals(Symmetric(Asym_noisy));

# ------------------------------------------------------------------------------------------
# Y así nuestros cálculos son mucho más eficientes :)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Un gran problema
# Usar matrices tridiagonales permite trabajar con problemas potencialmente muy grandes. El
# siguiente problema no seria posible resolverlo en una laptop si la matriz se tuviera que
# guardar como tipo `Matrix`.
# ------------------------------------------------------------------------------------------

n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)

# ------------------------------------------------------------------------------------------
# ### Álgebra lineal numérica
# La manera usual de agregar soporte para álgebra lineal numérica es haciendo un wrapper
# para subrutinas de BLAS y LAPACK. Para matrices con elementos `Float32`, `Float64`,
# `Complex{Float32}` ó `Complex{Float64}` esto es lo que hace Julia. Desde hace rato, Julia
# ha tnido soport para la multiplicación genérica de tipos. Así, cuando uno multiplica
# matrices de enteros, obtiene una matriz de enteros.
# ------------------------------------------------------------------------------------------

rand(1:10,3,3)*rand(1:10,3,3)

# ------------------------------------------------------------------------------------------
# Recientemente, más métodos de álgebra lineal se han añadido a Julia y ahora soporta
# factorizaciones generales de tipo `LU` y `QR`. Métodos generales para eigenvalores y SVD
# han sido escritos más recientemente en paquetería externa.
#
# En general, la factorización `LU` puede ser calculada cuando los elementos de la matriz se
# cierraon sobre los operadores `+`, `-`, `*` y `\`. Por supuesto, la matriz también deben
# tener rango completo. El método general de factorización `LU` en Julia aplica pivoteo y
# por lo tanto debe poder soportar `<` y `abs`. Por lo tanto es posible resolver sistemas de
# ecuaciones de e.g. números racionales como en los ejemplos que siguen.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Para usar números racionales, usamos un doble //:
# ------------------------------------------------------------------------------------------

1//2

# ------------------------------------------------------------------------------------------
# #### Ejemplo 1: Sistemas racionales lineales de ecuaciones
# Julia cuenta con números racionales ya instalados. El siguiente ejemplo consta de un
# sistema lineal de ecuaciones resulto sin promover a tipos de punto flotanto. Puede haber
# un error de overflow fácilmente al trabajar con racionales, así que usamos `BigInt`s.
# ------------------------------------------------------------------------------------------

Ar = convert(Matrix{Rational{BigInt}}, rand(1:10,3,3))/10

x = fill(1, (3))
b = Ar*x

Ar\b

lufact(Ar)

# ------------------------------------------------------------------------------------------
# #### Ejemplo 2: Matriz racional de eigenestructura
#
# El siguiente ejemplo muestra como la artimética de matriz racional puede ser usada para
# calcular una matriz dados los eigenvalores y eigenvectores racionales. Yo he encontrado
# ésto útil para mostrar ejemplos de sistemas dinámicos lineales.
# ------------------------------------------------------------------------------------------

λ1,λ2,λ3 = 1//1,1//2,1//4
v1,v2,v3 = [1,0,0],[1,1,0],[1,1,1]
V,Λ = [v1 v2 v3], Diagonal([λ1,λ2,λ3])
A = V*Λ/V

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# 11.1 ¿Cuáles son los eigenvalores de la Matriz A
#
# ```
# A =
# [
#  140   97   74  168  131
#   97  106   89  131   36
#   74   89  152  144   71
#  168  131  144   54  142
#  131   36   71  142   36
# ]
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 11.2 Crea una matriz diagonal de los eigenvalores de A
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 11.3 Realiza un factorización de Hessenberg sobre la matriz A. Verifica que `A = QHQ'`.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/intro-to-julia-ES/11.Factorizaciones y otras diversiones.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module DecFP

import Printf, SpecialFunctions

export Dec32, Dec64, Dec128, @d_str, @d32_str, @d64_str, @d128_str, exponent10, ldexp10

# Load libbid from our deps.jl
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("DecFP not installed properly, run Pkg.build(\"DecFP\"), restart Julia and try again")
end
include(depsjl_path)

const _buffer = fill(0x00, 1024)

import Base.promote_rule

#############################################################################
# exception handling via global flags
# (todo: recompile library with GLOBAL_FLAGS=0 for thread-safety)

const flags = Ref{Ptr{Cuint}}() # set to __bid_IDEC_glbflags in __init__

# clear exception flags and return x
function nox(x)
    unsafe_store!(flags[], 0)
    return x
end

# Check exception flags in mask & throw, otherwise returning x;
# always clearing exceptions.  This is a macros so that
# the error message is only evaluated if an exception occurs.
macro xchk(x, exc, args...)
    mask=0x3f
    if !isempty(args) && Meta.isexpr(args[end], :(=)) && args[end].args[1] == :mask # mask=... keyword at end
        mask = esc(args[end].args[2])
        args = args[1:end-1]
    end
    quote
        ret = $(esc(x))
        if $exc === nothing
            unsafe_store!(flags[], 0)
        else
            f = unsafe_load(flags[])
            unsafe_store!(flags[], 0)
            f & $mask != 0 && throw($exc($(map(esc,args)...)))
        end
        ret
    end
end

#############################################################################

@enum DecFPRoundingMode begin
    DecFPRoundNearest
    DecFPRoundDown
    DecFPRoundUp
    DecFPRoundToZero
    DecFPRoundFromZero
end

Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:Nearest})  = DecFPRoundNearest
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:Down})     = DecFPRoundDown
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:Up})       = DecFPRoundUp
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:ToZero})   = DecFPRoundToZero
Base.convert(::Type{DecFPRoundingMode}, ::RoundingMode{:FromZero}) = DecFPRoundFromZero

function Base.convert(::Type{RoundingMode}, r::DecFPRoundingMode)
    if r == DecFPRoundNearest
        return RoundNearest
    elseif r == DecFPRRoundDown
        return RoundDown
    elseif r == DecFPRRoundUp
        return RoundUp
    elseif r == DecFPRRoundToZero
        return RoundToZero
    elseif r == DecFPRRoundFromZero
        return RoundFromZero
    else
        throw(ArgumentError("invalid DecFP rounding mode code: $c"))
    end
end

const ROUNDING_PTR = Ref{Ptr{DecFPRoundingMode}}()

# global pointers and dicts must be initialized at runtime (via __init__)
function __init__()
    global ROUNDING_PTR[] = cglobal((:__bid_IDEC_glbround, libbid), DecFPRoundingMode) # rounding mode
    global flags[] = cglobal((:__bid_IDEC_glbflags, libbid), Cuint) # exception status
    unsafe_store!(flags[], 0)
end

# status flags from bid_functions.h:
const INVALID    = 0x01
const UNNORMAL   = 0x02
const DIVBYZERO  = 0x04
const OVERFLOW   = 0x08
const UNDERFLOW  = 0x10
const INEXACT    = 0x20

bidsym(w,s...) = string("__bid", w, "_", s...)

abstract type DecimalFloatingPoint <: AbstractFloat end

Base.Rounding.rounding_raw(::Type{T}) where {T<:DecimalFloatingPoint} =
    unsafe_load(ROUNDING_PTR[])
Base.Rounding.setrounding_raw(::Type{T}, r::DecFPRoundingMode) where {T<:DecimalFloatingPoint} =
    unsafe_store!(ROUNDING_PTR[],r)

Base.Rounding.rounding(::Type{T}) where {T<:DecimalFloatingPoint} =
    convert(RoundingMode, Base.Rounding.rounding_raw(T))
Base.Rounding.setrounding(::Type{T}, r::RoundingMode) where {T<:DecimalFloatingPoint} =
    Base.Rounding.setrounding_raw(T, convert(DecFPRoundingMode, r))

for w in (32,64,128)
    BID = Symbol(string("Dec",w))
    Ti = Symbol(string("UInt",w))
    @eval struct $BID <: DecimalFloatingPoint
        x::$Ti
        $BID(x::Number) = convert($BID, x)
        Base.reinterpret(::Type{$BID}, x::$Ti) = new(x)
    end
    # fix method ambiguities:
    @eval $BID(x::Rational{T}) where {T} = convert($BID, x)
end

# quickly check whether s begins with "±nan"
function isnanstr(s::AbstractString)
    st = iterate(s)
    c, i = '\0', 0
    while st !== nothing
        c, i = st
        isspace(c) || break
        st = iterate(s, i)
    end
    st === nothing && return false
    if c == '+' || c == '-'
        st = iterate(s, i)
        st === nothing && return false
        c, i = st
    end
    (c == 'n' || c == 'N') || return false
    st = iterate(s, i)
    if st !== nothing
        c, i = st
        if c == 'a' || c == 'A'
            st = iterate(s, i)
            if st !== nothing
                c, i = st
                if c == 'n' || c == 'N'
                    st = iterate(s, i)
                    if st === nothing
                        return true
                    end
                end
            end
        end
    end
    return false
end

"""
    exponent10(x::DecFP.DecimalFloatingPoint)

Get the exponent of the base 10 representation of a normalized floating-point number.

# Examples
```jldoctest
julia> exponent10(Dec64(123))
2
```
"""
exponent10(x::DecimalFloatingPoint)

"""
    ldexp10(x::DecFP.DecimalFloatingPoint, n::Integer)

Compute ``x * 10^n``.

# Examples
```jldoctest
julia> ldexp10(Dec64(15), 2)
1500.0
```
"""
ldexp10(x::DecFP.DecimalFloatingPoint, n::Integer)

# for compatibility with julia#29885
if isdefined(Base.Grisu, :getbuf)
    getdigitsbuf() = Base.Grisu.getbuf()
elseif isdefined(Base.Grisu, :DIGITSs)
    getdigitsbuf() = Base.Grisu.DIGITSs[Threads.threadid()]
else
    getdigitsbuf() = Base.Grisu.DIGITS
end

for w in (32,64,128)
    BID = Symbol(string("Dec",w))
    Ti = eval(Symbol(string("UInt",w)))
    T = eval(BID)

    # hack: we need an internal parsing function that doesn't check exceptions, since
    # flags isn't defined until __init__ runs.  Similarly for nextfloat/prevfloat
    @eval begin
        _parse(::Type{$BID}, s::AbstractString) =
            ccall(($(bidsym(w,"from_string")), libbid), $BID, (Ptr{UInt8},), s)
        _nextfloat(x::$BID) = ccall(($(bidsym(w,"nexttoward")), libbid), $BID, ($BID,Dec128), x, pinf128)
        _prevfloat(x::$BID) = ccall(($(bidsym(w,"nexttoward")), libbid), $BID, ($BID,Dec128), x, minf128)
        _sub(x::$BID, y::$BID) = ccall(($(bidsym(w,"sub")), libbid), $BID, ($BID,$BID), x, y)
    end

    @eval begin
        function Base.parse(::Type{$BID}, s::AbstractString)
            x = _parse($BID, s)
            if isnan(x) && !isnanstr(s)
                throw(ArgumentError("invalid number format $s"))
            end
            return @xchk(x, nothing)
        end

        $BID(x::AbstractString) = parse($BID, x)

        function tostring(x::$BID)
            # fills global _buffer
            ccall(($(bidsym(w,"to_string")), libbid), Cvoid, (Ptr{UInt8}, $BID), _buffer, x)
        end

        function Base.show(io::IO, x::$BID)
            isnan(x) && (print(io, "NaN"); return)
            isinf(x) && (print(io, signbit(x) ? "-Inf" : "Inf"); return)
            x == 0 && (print(io, signbit(x) ? "-0.0" : "0.0"); return)
            tostring(x)
            if _buffer[1] == UInt8('-')
                print(io, '-')
            end
            normalized_exponent = exponent10(x)
            lastdigitindex = findfirst(isequal(UInt8('E')), _buffer) - 1
            lastnonzeroindex = findlast(!isequal(UInt8('0')), view(_buffer, 1:lastdigitindex))
            if -5 < normalized_exponent < 6
                # %f
                if normalized_exponent >= 0
                    if normalized_exponent >= lastnonzeroindex - 2
                        GC.@preserve _buffer unsafe_write(io, pointer(_buffer, 2), lastnonzeroindex - 1)
                        printzeros(io, normalized_exponent - lastnonzeroindex + 2)
                        print(io, ".0")
                    else
                        GC.@preserve _buffer unsafe_write(io, pointer(_buffer, 2), normalized_exponent + 1)
                        print(io, '.')
                        GC.@preserve _buffer unsafe_write(io, pointer(_buffer, normalized_exponent + 3), lastnonzeroindex - normalized_exponent - 2)
                    end
                else
                    print(io, "0.")
                    printzeros(io, -normalized_exponent - 1)
                    GC.@preserve _buffer unsafe_write(io, pointer(_buffer, 2), lastnonzeroindex - 1)
                end
            else
                # %e
                print(io, Char(_buffer[2]), '.')
                if lastnonzeroindex == 2
                    print(io, '0')
                else
                    GC.@preserve _buffer unsafe_write(io, pointer(_buffer, 3), lastnonzeroindex - 2)
                end
                print(io, 'e')
                if normalized_exponent < 0
                    print(io, '-')
                    normalized_exponent = -normalized_exponent
                end
                b_lb = div(normalized_exponent, 10)
                b = 1
                while b <= b_lb
                    b *= 10
                end
                r = normalized_exponent
                while b > 0
                    q, r = divrem(r, b)
                    print(io, '0' + q)
                    b = div(b, 10)
                end
            end
            return
        end

        function Printf.Printf.fix_dec(x::$BID, n::Int, digits)
            if n > length(digits) - 1
                n = length(digits) - 1
            end
            rounded = round(ldexp10(x, n), RoundNearestTiesAway)
            if rounded == 0
                digits[1] = UInt8('0')
                return Int32(1), Int32(1), signbit(x)
            end
            tostring(rounded)
            trailing_zeros = 0
            i = 2
            while _buffer[i] != UInt8('E')
                digits[i - 1] = _buffer[i]
                if _buffer[i] == UInt8('0')
                    trailing_zeros += 1
                else
                    trailing_zeros = 0
                end
                i += 1
            end
            ndigits = i - 2
            len = ndigits - trailing_zeros
            i += 1
            if _buffer[i] == UInt8('+')
                expsign = +1
            elseif _buffer[i] == UInt8('-')
                expsign = -1
            end
            exponent = 0
            i += 1
            while _buffer[i] != 0x00
                exponent = exponent * 10 + _buffer[i] - UInt8('0')
                i += 1
            end
            exponent *= expsign
            pt = ndigits + exponent - n
            neg = signbit(x)
            return Int32(len), Int32(pt), neg
        end

        function Printf.Printf.ini_dec(x::$BID, n::Int, digits)
            if n > length(digits) - 1
                n = length(digits) - 1
            end
            if x == 0
                for i = 1:n
                    digits[i] = UInt8('0')
                end
                return Int32(1), Int32(1), signbit(x)
            end
            normalized_exponent = exponent10(x)
            rounded = round(ldexp10(x, n - 1 - normalized_exponent), RoundNearestTiesAway)
            rounded_exponent = exponent10(rounded)
            tostring(rounded)
            i = 2
            while _buffer[i] != UInt8('E')
                digits[i - 1] = _buffer[i]
                i += 1
            end
            while i <= n + 1
                digits[i - 1] = UInt8('0')
                i += 1
            end
            pt = normalized_exponent + rounded_exponent - n + 2
            neg = signbit(x)
            return Int32(n), Int32(pt), neg
        end

        # compatibility with julia#30373
        Printf.Printf.fix_dec(x::$BID, n::Int) = Printf.Printf.fix_dec(x, n, getdigitsbuf())
        Printf.Printf.ini_dec(x::$BID, n::Int) = Printf.Printf.ini_dec(x, n, getdigitsbuf())

        Base.fma(x::$BID, y::$BID, z::$BID) = nox(ccall(($(bidsym(w,"fma")), libbid), $BID, ($BID,$BID,$BID), x, y, z))
        Base.muladd(x::$BID, y::$BID, z::$BID) = fma(x,y,z) # faster than x+y*z

        Base.one(::Union{Type{$BID},$BID}) = $(_parse(T, "1"))
        Base.zero(::Union{Type{$BID},$BID}) = $(_parse(T, "0"))

        Base.signbit(x::$BID) = $(zero(Ti)) != $(Ti(1) << (Ti(w - 1))) & x.x
        Base.sign(x::$BID) = ifelse(signbit(x), $(_parse(T, "-1")), $(_parse(T, "1")))

        Base.nextfloat(x::$BID) = nox(_nextfloat(x))
        Base.prevfloat(x::$BID) = nox(_prevfloat(x))
        Base.eps(x::$BID) = ifelse(isfinite(x), @xchk(nextfloat(x) - x, OverflowError, "$($BID) value overflow", mask=OVERFLOW), $(_parse(T, "NaN")))

        # the meaning of the exponent is different than for binary FP: it is 10^n, not 2^n:
        exponent10(x::$BID) = nox(ccall(($(bidsym(w,"ilogb")), libbid), Cint, ($BID,), x))
        ldexp10(x::$BID, n::Integer) = nox(ccall(($(bidsym(w,"ldexp")), libbid), $BID, ($BID,Cint), x, n))
    end

    for (f,c) in ((:isnan,"isNaN"), (:isinf,"isInf"), (:isfinite,"isFinite"), (:issubnormal,"isSubnormal"))
        @eval Base.$f(x::$BID) = ccall(($(bidsym(w,c)), libbid), Cint, ($BID,), x) != 0
    end

    for (f,c) in ((:+,"add"), (:-,"sub"), (:*,"mul"), (:/, "div"), (:hypot,"hypot"), (:atan,"atan2"), (:^,"pow"), (:copysign,"copySign"))
        @eval Base.$f(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,$BID), x, y))
    end

    for f in (:exp,:log,:sin,:cos,:tan,:asin,:acos,:atan,:sinh,:cosh,:tanh,:asinh,:acosh,:atanh,:log1p,:expm1,:log10,:log2,:exp2,:exp10,:sqrt,:cbrt,:abs)
        @eval Base.$f(x::$BID) = @xchk(ccall(($(bidsym(w,f)), libbid), $BID, ($BID,), x), DomainError, x, mask=INVALID)
    end

    for (f,c) in ((:-,"negate"), (:trunc,"round_integral_zero"), (:floor,"round_integral_negative"), (:ceil,"round_integral_positive"), (:round,"nearbyint"))
        @eval Base.$f(x::$BID) = @xchk(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,), x), DomainError, x, mask=INVALID)
    end

    @eval SpecialFunctions.lgamma(x::$BID) = @xchk(ccall(($(bidsym(w,:lgamma)), libbid), $BID, ($BID,), x), DomainError, x, mask=INVALID)
    @eval SpecialFunctions.gamma(x::$BID) = @xchk(ccall(($(bidsym(w,:tgamma)), libbid), $BID, ($BID,), x), DomainError, x, mask=INVALID)

    for (r,c) in ((RoundingMode{:Nearest},"round_integral_nearest_even"), (RoundingMode{:NearestTiesAway},"round_integral_nearest_away"), (RoundingMode{:ToZero},"round_integral_zero"), (RoundingMode{:Up},"round_integral_positive"), (RoundingMode{:Down},"round_integral_negative"))
        @eval Base.round(x::$BID, ::$r) = @xchk(ccall(($(bidsym(w,c)), libbid), $BID, ($BID,), x), DomainError, x, mask=INVALID)
    end

    for (f,c) in ((:(==),"quiet_equal"), (:>,"quiet_greater"), (:<,"quiet_less"), (:(>=), "quiet_greater_equal"), (:(<=), "quiet_less_equal"))
        @eval Base.$f(x::$BID, y::$BID) = nox(ccall(($(bidsym(w,c)), libbid), Cint, ($BID,$BID), x, y) != 0)
    end

    for Tf in (Float32,Float64)
        bT = string("binary",sizeof(Tf)*8)
        @eval begin
            Base.convert(::Type{$Tf}, x::$BID) = nox(ccall(($(bidsym(w,"to_",bT)), libbid), $Tf, ($BID,), x))
            Base.$(Symbol("$Tf"))(x::$BID) = convert($Tf, x)
            Base.convert(::Type{$BID}, x::$Tf) = nox(ccall(($(string("__",bT,"_to_","bid",w)), libbid), $BID, ($Tf,), x))
        end
    end

    for c in (:π, :e, :ℯ, :γ, :catalan, :φ)
        @eval begin
            Base.convert(::Type{$BID}, ::Irrational{$(QuoteNode(c))}) = $(_parse(T, setprecision(256) do
                                                                                string(BigFloat(getfield(MathConstants, c)))
                                                                                end))
        end
    end

    @eval promote_rule(::Type{$BID}, ::Type{Irrational}) = $BID

    for w′ in (32,64,128)
        BID′ = Symbol(string("Dec",w′))
        if w > w′
            @eval promote_rule(::Type{$BID}, ::Type{$BID′}) = $BID
        end
        if w != w′
            @eval Base.convert(::Type{$BID}, x::$BID′) = @xchk(ccall(($(string("__bid",w′,"_to_","bid",w)), libbid), $BID, ($BID′,), x), nothing)
        end

        # promote binary*decimal -> decimal, for consistency with other operations above
        # (there doesn't seem to be any clear standard for this)
        if w′ <= 64
            FP′ = Symbol(string("Float",w′))
            @eval promote_rule(::Type{$BID}, ::Type{$FP′}) = $(Symbol(string("Dec",max(w,w′))))
            for (i′, i′str) in (("Int$w′", "int$w′"), ("UInt$w′", "uint$w′"))
                Ti′ = eval(Symbol(i′))
                @eval begin
                    Base.convert(::Type{$BID}, x::$Ti′) = nox(ccall(($(bidsym(w,"from_",i′str)), libbid), $BID, ($Ti′,), x))
                end
            end
        end
    end

    for w′ in (8,16,32,64)
        for (i′, i′str) in (("Int$w′", "int$w′"), ("UInt$w′", "uint$w′"))
            Ti′ = eval(Symbol(i′))
            @eval begin
                Base.trunc(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xint")), libbid), $Ti′, ($BID,), x), InexactError, :trunc, $BID, x, mask=INVALID | OVERFLOW)
                Base.floor(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xfloor")), libbid), $Ti′, ($BID,), x), InexactError, :floor, $BID, x, mask=INVALID | OVERFLOW)
                Base.ceil(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xceil")), libbid), $Ti′, ($BID,), x), InexactError, :ceil, $BID, x, mask=INVALID | OVERFLOW)
                Base.round(::Type{$Ti′}, x::$BID, ::RoundingMode{:NearestTiesAway}) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xrninta")), libbid), $Ti′, ($BID,), x), InexactError, :round, $BID, x, mask=INVALID | OVERFLOW)
                Base.convert(::Type{$Ti′}, x::$BID) = @xchk(ccall(($(bidsym(w,"to_",i′str,"_xfloor")), libbid), $Ti′, ($BID,), x), InexactError, :convert, $BID, x)
                Base.$(Symbol("$Ti′"))(x::$BID) = convert($Ti′, x)
            end
        end
    end

    @eval Base.bswap(x::$BID) = reinterpret($BID, bswap(x.x))
    @eval Base.convert(::Type{Float16}, x::$BID) = convert(Float16, convert(Float32, x))
    @eval Base.Float16(x::$BID) = convert(Float16, x)
    @eval Base.reinterpret(::Type{$Ti}, x::$BID) = x.x
end # widths w

Base.round(x::DecimalFloatingPoint, ::RoundingMode{:FromZero}) = signbit(x) ? floor(x) : ceil(x)

Base.trunc(::Type{Integer}, x::DecimalFloatingPoint) = trunc(Int, x)
Base.floor(::Type{Integer}, x::DecimalFloatingPoint) = floor(Int, x)
Base.ceil(::Type{Integer}, x::DecimalFloatingPoint) = ceil(Int, x)
Base.round(::Type{Integer}, x::DecimalFloatingPoint) = round(Int, x)
Base.round(::Type{Integer}, x::DecimalFloatingPoint, ::RoundingMode{:NearestTiesAway}) = round(Int, x, RoundNearestTiesAway)
Base.convert(::Type{Integer}, x::DecimalFloatingPoint) = convert(Int, x)

Base.round(::Type{T}, x::DecimalFloatingPoint) where {T<:Integer} = convert(T, round(x))
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:Nearest}) where {T<:Integer} = convert(T, round(x, RoundNearest))
function Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:NearestTiesUp}) where {T<:Integer}
    y = floor(T, x)
    ifelse(x==y, y, copysign(floor(T, 2*x-y), x))
end
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:ToZero}) where {T<:Integer} = trunc(T, x)
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:FromZero}) where {T<:Integer} = (x>=0 ? ceil(T, x) : floor(T, x))
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:Up}) where {T<:Integer} = ceil(T, x)
Base.round(::Type{T}, x::DecimalFloatingPoint, ::RoundingMode{:Down}) where {T<:Integer} = floor(T, x)

# the complex-sqrt function in base doesn't work for use, because it requires base-2 ldexp
function Base.sqrt(z::Complex{T}) where {T<:DecimalFloatingPoint}
    x, y = reim(z)
    x==y==0 && return Complex(zero(x),y)
    ρ = sqrt((abs(x) + hypot(x,y)) * 0.5)
    ξ = ρ
    η = y
    if isfinite(η) η=(η/ρ)/2 end
    if x<0
        ξ = abs(η)
        η = copysign(ρ,y)
    end
    Complex(ξ,η)
end

# see issue #92 — the fallback power_by_squaring fails for p < 0, and this is more accurate:
Base.:^(x::DecimalFloatingPoint, p::Integer) = x^oftype(x, p)
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{0}) = one(x)
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{1}) = x
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{2}) = x*x
@inline Base.literal_pow(::typeof(^), x::DecimalFloatingPoint, ::Val{3}) = x*x*x

# used for next/prevfloat:
const pinf128 = _parse(Dec128, "+Inf")
const minf128 = _parse(Dec128, "-Inf")

for T in (Dec32, Dec64, Dec128)
    @eval begin
        Base.eps(::Type{$T}) = $(_sub(_nextfloat(one(T)), one(T)))
        Base.typemax(::Type{$T}) = $(_parse(T, "+inf"))
        Base.typemin(::Type{$T}) = $(_parse(T, "-inf"))
        Base.floatmax(::Type{$T}) = $(_prevfloat(_parse(T, "+inf")))
        Base.floatmin(::Type{$T}) = $(_nextfloat(zero(T)))
    end
end

Base.maxintfloat(::Type{Dec32}) = reinterpret(Dec32, 0x36000001) # Dec32("1e7")
Base.maxintfloat(::Type{Dec64}) = reinterpret(Dec64, 0x33c0000000000001) # Dec64("1e16")
Base.maxintfloat(::Type{Dec128}) = reinterpret(Dec128, 0x30840000000000000000000000000001) # Dec128("1e34")

Base.convert(T::Type{F}, x::Union{Int8,UInt8,Int16,UInt16}) where {F<:DecimalFloatingPoint} = F(Int32(x))
Base.convert(T::Type{F}, x::Integer) where {F<:DecimalFloatingPoint} = F(Int64(x))
Base.convert(T::Type{F}, x::Unsigned) where {F<:DecimalFloatingPoint} = F(UInt64(x))
Base.convert(T::Type{F}, x::Rational) where {F<:DecimalFloatingPoint} = F(x.num) / F(x.den)
Base.convert(T::Type{F}, x::Float16) where {F<:DecimalFloatingPoint} = F(Float32(x))
promote_rule(::Type{F}, ::Type{Float16}) where {F<:DecimalFloatingPoint} = F
promote_rule(::Type{F}, ::Type{T}) where {F<:DecimalFloatingPoint,T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64}} = F

# so that mathconsts get promoted to Dec32, not Dec64, like Float32
promote_rule(::Type{Irrational{s}}, ::Type{F}) where {s,F<:DecimalFloatingPoint} = F
promote_rule(::Type{Irrational{s}}, T::Type{Complex{F}}) where {s,F<:DecimalFloatingPoint} = T

macro d_str(s, flags...) parse(Dec64, s) end
macro d32_str(s, flags...) parse(Dec32, s) end
macro d64_str(s, flags...) parse(Dec64, s) end
macro d128_str(s, flags...) parse(Dec128, s) end

# for zero-padding in printing routines above
function printzeros(io::IO, n::Int)
    for i = 1:n
        print(io, '0')
    end
end

end # module


text:
using FractalTools

# An example of one dimensional ifs
A1 = reshape([1/2],1,1)
b1 = [-1/8]
A2 = reshape([1/2],1,1)
b2 = [1/8]

w1 = Transformation(A1,b1)
w2 = Transformation(A2,b2)

w = [w1, w2]
ifs = IFS(w)

take!(ifs.generator)


text:
using DelaySSAToolkit
using Catalyst
# using Plots
################################################################################
# the following example illustrates how to deal with a oscillatory system
# d: X-> 0; k1/(1+Y^2): 0-> X; [trigger X-> Y after τ time;] k2*Y/(1+Y):  Y-> 0;

rn = @reaction_network begin
    1/(1+Y^2), 0 --> X
    1/(1+Y),   Y --> 0
end

jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)

u0 = [0,0]
de_chan0 = [[]]
tf = 400.
tspan = (0,tf)
τ = 20.
dprob = DiscreteProblem(jumpsys, u0, tspan)
# jumpsys.dep_graph

delay_trigger_affect! = function (integrator, rng)
  append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(1=>delay_trigger_affect!)
delay_complete = Dict(1=>[2=>1, 1=>-1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
# alg = DelayRejection()
# alg = DelayDirect()
# alg = DelayMNRM()
alg = DelayDirectCR()
djprob = DelayJumpProblem(jumpsys, dprob, alg, delayjumpset, de_chan0)
sol = solve(djprob, SSAStepper(), seed = 12345)
using Plots
plot(sol.t,[sol.u[i][1] for i in eachindex(sol.u)],alpha=0.3, label="X",color="green",linewidth = 2, legend = :top, ylabel = "# of individuals", xlabel = "Time", fmt=:svg)

sol = solve(djprob, SSAStepper(), seed = 1234)
plot!(sol.t,[sol.u[i][1] for i in eachindex(sol.u)], label="X",linewidth = 2, line=:dash,color="green")

Sample_size = Int(1e4)
ens_prob = EnsembleProblem(djprob)
ens = solve(ens_prob, SSAStepper(), EnsembleThreads(), trajectories = Sample_size, saveat = .1)
using StatsBase
tsm(t) = timepoint_mean(ens,t)
mean_X(t) = tsm(t)[1]
mean_Y(t) = tsm(t)[2]

timestamps=0:1.0:tf
plot(timestamps,mean_X.(timestamps),linewidth=3,line=:dash,label="X",xlabel="time",ylabel="Mean Value")
plot!(timestamps,mean_Y.(timestamps),linewidth=3,line=:dash, legend = :topright,label="Y")

# plot(mean_X.(timestamps),mean_Y.(timestamps))
# might be very slow, can try with small end point
# @gif for i in 1:200
#     plot(mean_X,mean_Y,0,i, linewidth=range(0, 5, length = 200),seriesalpha=range(0, 1, length = 200),xlim=(0,21),ylim=(0,7),label=false,xlabel="X",ylabel="Y")
# end


text:
# libminknap methods
macro pk_min_ccall(func, args...)
    args = map(esc, args)
    f = "$func"
    quote
        ccall(($f, PisingerKnapsack._jl_libminknap), $(args...))
    end
end

function minknap(p::Vector{Cint}, w::Vector{Cint}, capacity::Cint)
    nbitems = length(p)
    solution = fill(Cint(0), nbitems)
    obj =
        @pk_min_ccall minknap Clong (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint) nbitems p w solution capacity
    return (obj, solution)
end

function minknap(p::Vector, w::Vector, capacity)
    return minknap(convert(Vector{Cint}, p), convert(Vector{Cint}, w), Cint(capacity))
end

# libbouknap methods
macro pk_bou_ccall(func, args...)
    args = map(esc, args)
    f = "$func"
    quote
        ccall(($f, PisingerKnapsack._jl_libbouknap), $(args...))
    end
end

function bouknap(p::Vector{Cint}, w::Vector{Cint}, ub::Vector{Cint}, capacity::Cint)
    nbitems = length(p)
    solution = fill(Cint(0), nbitems)
    obj =
        @pk_bou_ccall bouknap Clong (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint) nbitems p w ub solution capacity
    return (obj, solution)
end

function bouknap(p::Vector, w::Vector, ub::Vector, capacity)
    return bouknap(
        convert(Vector{Cint}, p),
        convert(Vector{Cint}, w),
        convert(Vector{Cint}, ub),
        Cint(capacity),
    )
end

# libmulknap methods
macro pk_mul_ccall(func, args...)
    args = map(esc, args)
    f = "$func"
    quote
        ccall(($f, PisingerKnapsack._jl_libmulknap), $(args...))
    end
end

function mulknap(p::Vector{Cint}, w::Vector{Cint}, capacities::Vector{Cint})
    nbitems = length(p)
    nbknapsacks = length(capacities)
    solution = fill(Cint(0), nbitems)
    obj =
        @pk_mul_ccall mulknap Clong (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}) nbitems nbknapsacks p w solution capacities
    return (obj, solution)
end

function mulknap(p::Vector, w::Vector, capacities::Vector)
    return mulknap(
        convert(Vector{Cint}, p),
        convert(Vector{Cint}, w),
        convert(Vector{Cint}, capacities),
    )
end

# FYI : mcknap function generates random instances so we cannot use as provided by Pisinger

# macro pk_minmc_ccall(func, args...)
#     args = map(esc, args)
#     f = "$func"
#     quote
#         ccall(($f, PisingerKnapsack._jl_libminmcknap), $(args...))
#     end
# end

# libminmcknap functions
# function minmcknap(p::Vector{Cint}, w::Vector{Cint}, capacity::Cint, items_per_class::Vector{Cint})
#     nbitems = length(p)
#     nbclasses = length(items_per_class)
#     solution = fill(Cint(0), nbitems)
#     obj = @pk_minmc_ccall minmcknap Clong (Cint, Cint, Ptr{Cint}, Ptr{Cint},
#         Ptr{Cint}, Ptr{Cint}, Cint) nbitems nbclasses items_per_class p w solution capacity
#     return (obj, solution)
# end

# function minmcknap(p::Vector, w::Vector, capacity, items_per_class::Vector)
#     return minmcknap(convert(Vector{Cint}, p), convert(Vector{Cint}, w), Cint(capacity), convert(Vector{Cint}, items_per_class))
# end


text:
using StatsBase
using Yao.AD
using Yao.AD: Rotor

############# manipulate quantum differential node ############
Yao.AD.generator(c::CPhaseGate{N}) where N = ControlBlock{N}(c.ctrl_locs, c.ctrl_config, Z, c.locs)

function get_diffblocks(circuit::AbstractBlock)
    diffblocks = collect_blocks(Union{RotationGate, CPhaseGate}, circuit)
    length(diffblocks) == nparameters(circuit) || throw(ArgumentError("Input circuit contains non-differentiable/unsupported parameters."))
    return diffblocks
end

#### interface #####
export numdiff, faithful_grad

@inline function _perturb(func, gate::AbstractBlock, δ::Real)
    dispatch!(-, gate, (δ,))
    r1 = func()
    dispatch!(+, gate, (2δ,))
    r2 = func()
    dispatch!(-, gate, (δ,))
    r1, r2
end

@inline function _expect(op::AbstractBlock, reg::ArrayReg; nshots=nothing)
    if nshots === nothing
        expect(op, reg)
    else
        mean(measure(op, reg; nshots=nshots))
    end
end

"""
    numdiff(loss, circuit::AbstractBlock, δ::Real=1e-2) => Vector

Numeric differentiation a loss over a circuit, the loss take the circuit as input.
"""
@inline function numdiff(loss, circuit::AbstractBlock, δ::Real=1e-2)
    map(get_diffblocks(circuit)) do diffblock
        r1, r2 = _perturb(()->loss(circuit), diffblock, δ)
        (r2 - r1)/2δ
    end
end

"""
    faithful_grad(x, pair::Pair{<:ArrayReg, <:AbstractBlock}; nshots=nothing) -> Vector

Differentiate `x` over all parameters. `x` can be an `AbstractBlock`, `StatFunctional` or `MMD`.
"""
@inline function faithful_grad(op::AbstractBlock, pair::Pair{<:ArrayReg, <:AbstractBlock}; nshots=nothing)
    map(get_diffblocks(pair.second)) do diffblock
        r1, r2 = _perturb(()->_expect(op, copy(pair.first) |> pair.second; nshots=nshots) |> real, diffblock, π/2)
        (r2 - r1)/2
    end
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
using InfiniteLinearAlgebra, InfiniteArrays, BandedMatrices, LazyArrays, FillArrays, LinearAlgebra, Test
import InfiniteLinearAlgebra: LowerHessenbergQ, tail_de, toeptail, InfToeplitz, PertToeplitz
import BandedMatrices: _BandedMatrix

@testset "Inf QL" begin
    @testset "Toeplitz QLHessenberg" begin
        @testset "Tridiagonal Toeplitz" begin
            for (Z,A,B) in ((2.0,5.1,0.5), (2.0,2.2,0.5), (2.0,-2.2,0.5), (2.0,-5.1,0.5),
                        (0.5,5.1,2.0),  (0.5,-5.1,2.0))
                n = 100_000; T = Tridiagonal(Fill(Z,∞), Fill(A,∞), Fill(B,∞)); Qn,Ln = ql(BandedMatrix(T)[1:n,1:n]);
                Q,L = ql(T)
                @test L[1:10,1:10] ≈ Ln[1:10,1:10]
                @test Q[1:10,1:10] ≈ Qn[1:10,1:10]
                @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
            end

            for (Z,A,B)  in ((2,5.0im,0.5),
                        (2,2.1+0.1im,0.5),
                        (2,2.1-0.1im,0.5),
                        (2,1.5+0.1im,0.5),
                        (2,1.5+0.0im,0.5),
                        (2,1.5-0.0im,0.5),
                        (2,0.0+0.0im,0.5),
                        (2,-0.1+0.1im,0.5),
                        (2,-1.1+0.1im,0.5),
                        (2,-0.1-0.1im,0.5))
                T = Tridiagonal(Fill(ComplexF64(Z),∞), Fill(ComplexF64(A),∞), Fill(ComplexF64(B),∞))
                Q,L = ql(T)
                @test Q[1:10,1:12] * L[1:12,1:10] ≈ T[1:10,1:10]
            end

            for (Z,A,B) in ((0.5,2.1,2.0),(0.5,-1.1,2.0))
                @test_throws DomainError ql(Tridiagonal(Fill(Z,∞), Fill(A,∞), Fill(B,∞)))
            end
        end

        @testset "Hessenberg Toeplitz" begin
            a = [1,2,3,0.5]
            T = _BandedMatrix(reverse(a) * Ones(1,∞), ℵ₀, 2, 1)
            F = ql(T)
            @test F.Q[1:10,1:11]*F.L[1:11,1:10] ≈ T[1:10,1:10]

            a = [1,2,3+im,0.5]
            T = _BandedMatrix(reverse(a) * Ones{eltype(a)}(1,∞), ℵ₀, 2, 1)
            Q,L = ql(T)
            @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
            @test T isa InfToeplitz

            T = BandedMatrix(-2 => Fill(1,∞), 0 => Fill(0.5+eps()im,∞), 1 => Fill(0.25,∞))
            Q,L = ql(T)
            @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]

            for T in (BandedMatrix(-2 => Fill(1,∞), 0 => Fill(0.5,∞), 1 => Fill(0.25,∞)),
                        BandedMatrix(-2 => Fill(1/4,∞), 1 => Fill(1,∞))-im*I)
                Q,L = ql(T)
                @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
            end

            a =    [ -2.531640004434771-0.0im , 0.36995310821558014+2.5612894011525276im, -0.22944284364953327+0.39386202384951985im, -0.2700241133710857 + 0.8984628598798804im, 4.930380657631324e-32 + 0.553001215633963im ]
            T = _BandedMatrix(a * Ones{ComplexF64}(1,∞), ℵ₀, 3, 1)
            Q,L = ql(T)
            @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
        end
    end

    @testset "PertTriToeplitz QL" begin
        A = Tridiagonal(Vcat(Float64[], Fill(2.0,∞)),
                        Vcat(Float64[2.0], Fill(0.0,∞)),
                        Vcat(Float64[], Fill(0.5,∞)))
        for λ in (-2.1-0.01im,-2.1+0.01im,-2.1+eps()im,-2.1-eps()im,-2.1+0.0im,-2.1-0.0im,-1.0+im,-3.0+im,-3.0-im)
            Q, L = ql(A - λ*I)
            @test Q[1:10,1:12]*L[1:12,1:10] ≈ A[1:10,1:10] - λ*I
        end
    end

    @testset "Pert Hessenberg Toeplitz" begin
        a = [1,2,5,0.5]
        Random.seed!(0)
        A = _BandedMatrix(Hcat(randn(4,2), reverse(a) * Ones(1,∞)), ℵ₀, 2, 1)
        @test A isa PertToeplitz
        @test BandedMatrix(A, (3,1))[1:10,1:10] == A[1:10,1:10]
        Q,L = ql(A)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ A[1:10,1:10]

        a = [0.1,1,2,3,0.5]
        A = _BandedMatrix(Hcat([0.5 0.5; -1 3; 2 2; 1 1; 0.1 0.1], reverse(a) * Ones(1,∞)), ℵ₀, 3, 1)
        @test A isa PertToeplitz
        @test BandedMatrix(A, (3,1))[1:10,1:10] == A[1:10,1:10]

        B = BandedMatrix(A, (3,1))
        @test B[1:10,1:10] == A[1:10,1:10]
        Q,L = ql(A)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ A[1:10,1:10]

        A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),
                        0 => Vcat([0.0], Fill(1/2,∞)),
                        1 => Vcat([1/4], Fill(1/4,∞)))
        Q,L = ql(A)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ A[1:10,1:10]
        a = [-0.1,0.2,0.3]
        A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),  0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
        λ = 0.5+0.1im

        B = BandedMatrix(A-λ*I, (3,1))
        T = toeptail(B)
        Q,L = ql(T)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]

        Q,L = ql(A-λ*I)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ (A-λ*I)[1:10,1:10]

        a = [-0.1,0.2,0.3]
        A = BandedMatrix(-2 => Vcat([1], Fill(1,∞)),  0 => Vcat(a, Fill(0,∞)), 1 => Vcat([1/4], Fill(1/4,∞)))
        for λ in (0.1+0im,0.1-0im, 3.0, 3.0+1im, 3.0-im, -0.1+0im, -0.1-0im)
            Q, L = ql(A-λ*I)
            @test Q[1:10,1:11]*L[1:11,1:10] ≈ (A-λ*I)[1:10,1:10]
        end
    end

    @testset "Pert faux-periodic QL" begin
        a = [0.5794879759059747 + 0.0im,0.538107104952824 - 0.951620830938543im,-0.19352887774167749 - 0.3738926065520737im,0.4314153362874331,0.0]
        T = _BandedMatrix(a*Ones{ComplexF64}(1,∞), ℵ₀, 3,1)
        Q,L = ql(T)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ T[1:10,1:10]
        Qn,Ln = ql(T[1:1001,1:1001])
        @test Qn[1:10,1:10] * diagm(0=>[-1; (-1).^(1:9)]) ≈ Q[1:10,1:10]
        @test diagm(0=>[-1; (-1).^(1:9)]) * Ln[1:10,1:10] ≈ L[1:10,1:10]

        f =  [0.0+0.0im        0.522787+0.0im     ; 0.59647-1.05483im    0.538107-0.951621im; -0.193529-0.373893im  -0.193529-0.373893im;
                0.431415+0.0im        0.431415+0.0im; 0.0+0.0im             0.0+0.0im]
        A = _BandedMatrix(Hcat(f, a*Ones{ComplexF64}(1,∞)), ℵ₀, 3,1)
        Q,L = ql(A)
        @test Q[1:10,1:11]*L[1:11,1:10] ≈ A[1:10,1:10]
        Qn,Ln = ql(A[1:1000,1:1000])
        @test Qn[1:10,1:10] * diagm(0 => [Ones(5); -(-1).^(1:5)]) ≈ Q[1:10,1:10]
        @test diagm(0 => [Ones(5); -(-1).^(1:5)]) * Ln[1:10,1:10] ≈ L[1:10,1:10]
    end

    @testset "solve with QL" begin
        A = BandedMatrix(-1 => Fill(2,∞), 0 => Fill(5,∞), 1 => Fill(0.5,∞))
        @test (qr(A)\Vcat(1.0,Zeros(∞)))[1:1000] ≈ (ql(A)\Vcat(1.0,Zeros(∞)))[1:1000]

        J = BandedMatrix(0 => Vcat([1.0], Fill(0.0,∞)), 1 => Vcat(Float64[],Fill(0.5,∞)), -1 => Vcat(Float64[],Fill(0.5,∞)))
        z = 3.5
        A = J - z*I
        F = ql(A)
        Q,L = F
        b = [1; zeros(∞)]
        @test (Q'b)[1] ≈ 0.9892996329463546
        @test L[1] == L[1,1]

        @test L[1:5,1:5] isa BandedMatrix
        @test (L*(L \ b))[1:10] ≈ [1; zeros(9)]
        u = F \ b
        @test (A*u)[1:10] ≈ [1; zeros(9)]
    end

    @testset "Derivation" begin
        A = BandedMatrix(-1 => Fill(2,∞), 0 => Fill(5,∞), 1 => Fill(0.5,∞))
        a = reverse(A.data.args[1])
        d,e = tail_de(a)
        X = [transpose(a); 0 d e]
        Q = LowerHessenbergQ(Fill(ql!(X).Q,∞))
        L = _BandedMatrix(Hcat([e; X[2,2]; X[2,1]], X[2,end:-1:1] * Ones{Float64}(1,∞)), ℵ₀, 2, 0)
        Qn,Ln = ql(A[1:1000,1:1000])
        @test Q[1:10,1:10] ≈ Qn[1:10,1:10]
        @test Q'A isa MulMatrix
        @test Array((Q'A)[1:10,1:10]) ≈ [(Q'A)[k,j] for k=1:10,j=1:10]
        @test (Q'A)[1:10,1:10] ≈ Ln[1:10,1:10] ≈ L[1:10,1:10]

        A = BandedMatrix(-1 => Fill(2,∞), 0 => Fill(5+im,∞), 1 => Fill(0.5,∞))
        a = reverse(A.data.args[1])
        d,e = tail_de(a)
        X = [transpose(a); 0 d e]
        q = ql!(X).Q
        Q = LowerHessenbergQ(Fill(q,∞))
        L = _BandedMatrix(Hcat([e; X[2,2]; X[2,1]], X[2,end:-1:1] * Ones{Float64}(1,∞)), ℵ₀, 2, 0)
        Qn,Ln = ql(A[1:1000,1:1000])
        @test Q[1:10,1:10] ≈ Qn[1:10,1:10] * diagm(0 => [1; -Ones(9)] )
        @test (Q'A)[1:10,1:10] ≈ diagm(0 => [1; -Ones(9)] ) * Ln[1:10,1:10] ≈ L[1:10,1:10]
    end
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using ScoreDrivenERGM
using LightGraphs
using Distributions

# define test matrix
N = 100
A = rand(Bool,N,N)

g = SimpleDiGraph(A)


text:
using PGFPlotsX
@pgf a = Axis(
    {
        ybar,
        enlargelimits = 0.05,
        legend_style =
        {
            at = Coordinate(0.5, -0.15),
            anchor = "north",
            legend_columns = -1
        },
        ylabel = "Storage in megabyte",
        symbolic_x_coords=["MeshCore D", "MeshCore D/U", "MeshCore 32", "MeshCore 64", "MDS-RED", "MOAB", "MDS", "GRUMMP", "STK"],
        xtick = "data",
        xticklabel_style={
        rotate=45,
        anchor="east"
        },
        nodes_near_coords,
        nodes_near_coords_align={vertical},
    },
    Plot(Coordinates([("MeshCore D", 1.99), ("MeshCore D/U", 3.7), ("MeshCore 32", 15.3), ("MeshCore 64",  27.47), ("MDS-RED", 6.4), ("MOAB", 10), ("MDS", 23.4), ("GRUMMP", 39.8), ("STK", 67.6)]))
)

display(a)

# The MDS system of
# Figure 9 of~\cite{Ibanez2016}: MOAB 10 MB per 100K tets
# We need 3*8*100000/6/1024/1024+4*8*100000/1024/1024=3.5 MB (non--modifiable MOAB
# 10 MB, implementations support adaptive meshes PUMI 23 MB, GRUMMP 40 MB,  STK 67
# MB).


text:
using Test, Documenter, COVID_Tracking_Clean

@testset "Dry Run" begin
    @test COVID_Tracking_Clean.main(false)
end

@testset "Documentation" begin
    using Documenter, COVID_Tracking_Clean

    DocMeta.setdocmeta!(COVID_Tracking_Clean,
                        :DocTestSetup,
                        :(using COVID_Tracking_Clean;),
                        recursive = true)
    makedocs(sitename = "COVID_Tracking_Clean",
             modules = [COVID_Tracking_Clean],
             pages = [
                 "Home" => "index.md",
                 "API" => "api.md"
             ],
             source = joinpath("..", "docs", "src"),
             build = joinpath("..", "docs", "build"),
             )
    @test true
end