text:
function ModelingToolkit.structural_simplify(sys::ODESystem)
    sys = initialize_system_structure(alias_elimination(sys))
    # check_consistency(sys)
    if sys isa ODESystem
        sys = dae_index_lowering(sys)
    end
    sys = tearing(sys)
    fullstates = [map(eq->eq.lhs, observed(sys)); states(sys)]
    ModelingToolkit.@set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullstates)
    return sys
end

function ModelingToolkit.initialize_system_structure(sys::ODESystem)
    sys = flatten(sys)

    iv = independent_variable(sys)
    eqs = copy(equations(sys))
    neqs = length(eqs)
    algeqs = trues(neqs)
    dervaridxs = OrderedSet{Int}()
    var2idx = Dict{Any,Int}()
    symbolic_incidence = []
    fullvars = []
    var_counter = Ref(0)
    addvar! = let fullvars=fullvars, var_counter=var_counter
        var -> begin
            get!(var2idx, var) do
                push!(fullvars, var)
                var_counter[] += 1
            end
        end
    end

    vars = OrderedSet()
    for (i, eq′) in enumerate(eqs)
        if ModelingToolkit._iszero(eq′.lhs)
            eq = eq′
        else
            eq = 0 ~ eq′.rhs - eq′.lhs
        end
        ModelingToolkit.vars!(vars, eq.rhs)
        isalgeq = true
        statevars = []
        for var in vars
            hasmetadata(var, VariableOutput) && continue
            isequal(var, iv) && continue
            if ModelingToolkit.isparameter(var) || (ModelingToolkit.istree(var) && ModelingToolkit.isparameter(ModelingToolkit.operation(var)))
                continue
            end
            varidx = addvar!(var)
            push!(statevars, var)

            dvar = var
            idx = varidx
            while ModelingToolkit.isdifferential(dvar)
                if !(idx in dervaridxs)
                    push!(dervaridxs, idx)
                end
                isalgeq = false
                dvar = ModelingToolkit.arguments(dvar)[1]
                idx = addvar!(dvar)
            end
        end
        push!(symbolic_incidence, copy(statevars))
        empty!(statevars)
        empty!(vars)
        algeqs[i] = isalgeq
        if isalgeq
            eqs[i] = eq
        end
    end

    # sort `fullvars` such that the mass matrix is as diagonal as possible.
    dervaridxs = collect(dervaridxs)
    sorted_fullvars = OrderedSet(fullvars[dervaridxs])
    for dervaridx in dervaridxs
        dervar = fullvars[dervaridx]
        diffvar = ModelingToolkit.arguments(dervar)[1]
        if !(diffvar in sorted_fullvars)
            push!(sorted_fullvars, diffvar)
        end
    end
    for v in fullvars
        if !(v in sorted_fullvars)
            push!(sorted_fullvars, v)
        end
    end
    fullvars = collect(sorted_fullvars)
    var2idx = Dict(fullvars .=> eachindex(fullvars))
    dervaridxs = 1:length(dervaridxs)

    nvars = length(fullvars)
    diffvars = []
    vartype = fill(ModelingToolkit.SystemStructures.DIFFERENTIAL_VARIABLE, nvars)
    varassoc = zeros(Int, nvars)
    inv_varassoc = zeros(Int, nvars)
    for dervaridx in dervaridxs
        vartype[dervaridx] = ModelingToolkit.SystemStructures.DERIVATIVE_VARIABLE
        dervar = fullvars[dervaridx]
        diffvar = ModelingToolkit.arguments(dervar)[1]
        diffvaridx = var2idx[diffvar]
        push!(diffvars, diffvar)
        varassoc[diffvaridx] = dervaridx
        inv_varassoc[dervaridx] = diffvaridx
    end

    algvars = setdiff(states(sys), diffvars)
    for algvar in algvars
        # hasmetadata(algvar, VariableOutput) && continue
        # it could be that a variable appeared in the states, but never appeared
        # in the equations.
        algvaridx = get(var2idx, algvar, 0)

        algvaridx == 0 && throw(ModelingToolkit.SystemStructures.InvalidSystemException("The system is missing "
            * "an equation for $algvar."
        ))
        vartype[algvaridx] = ModelingToolkit.SystemStructures.ALGEBRAIC_VARIABLE
    end

    graph = BipartiteGraph(neqs, nvars, Val(false))
    for (ie, vars) in enumerate(symbolic_incidence), v in vars
        jv = var2idx[v]
        LightGraphs.add_edge!(graph, ie, jv)
    end

    ModelingToolkit.@set! sys.eqs = eqs
    ModelingToolkit.@set! sys.structure = SystemStructure(
        fullvars = fullvars,
        vartype = vartype,
        varassoc = varassoc,
        inv_varassoc = inv_varassoc,
        varmask = iszero.(varassoc),
        algeqs = algeqs,
        graph = graph,
        solvable_graph = BipartiteGraph(ModelingToolkit.nsrcs(graph), ModelingToolkit.ndsts(graph), Val(false)),
        assign = Int[],
        inv_assign = Int[],
        scc = Vector{Int}[],
        partitions = ModelingToolkit.SystemStructures.SystemPartition[],
    )
    return sys
end

@variables v(t)=0.0 I(t)=0.0 out1(t)=0.0 out2(t)=0.0 dout1(t) dout2(t)
eqs = [
    D(v) ~ 0.1 + I
    I ~ 0.35*v*v
    D(out1) ~ dout1
    dout1 ~ 0.1 + I
    D(out2) ~ dout2
    out2 ~ I
]
test = ODESystem(eqs,t;name=:test)
sys = structural_simplify(test)
defs = ModelingToolkit.get_defaults(sys)
prob = ODAEProblem(sys,defs,(0.0,1.0))
sol = solve(prob,Tsit5())
plot(sol)
plot(sol.t, sol[1,:])
plot(sol, vars=[out1])
plot(sol, vars=[I])
plot(sol, vars=[out2])
plot(sol, vars=[dout1])

eqs = [D(v) ~ 0.1]
observed = [o1.x ~ v
            I ~ 0.1v
            o2.x ~ I]

=>

eqs = [D(v) ~ dv
       D(o1.x) ~ dv
       D(o2.x) ~ 0.1*dv]
observed = [I ~ 0.1*v
            dv ~ 0.1]

# outpins = [o1.x, o2.x]
# diffeqs = filter(e -> ModelingToolkit.isdiffeq(e), eqs)
# diffeqdict = Dict([e.lhs => e.rhs for e in diffeqs])
function move_observed_state_output(oeq)
    # create variable dv(iv...)
    # push!(observed, dv ~ diffeqdict[D(oeq.rhs)])
    # push!(eqs, D(oeq.lhs) ~ dv
    # remove oeq from observed
end
function move_observed_intermediate_output(oeq)
    # create variable 'du1x'
    # for eq in observed
    #   eq.lhs != oeq.rhs && continue
    #   rhs = eq.rhs
    #   diff_terms = []
    #   new_rhs = substutude
    #   new_oeq = du1x ~ new_rhs
    # end
    # push!(observed, ds1x ~ diffeqdict[D(oeq.rhs)])
    # push!(eqs, D(oeq.lhs) ~ du1x
    # remove oeq from observed
end
function create_states_for_outpins()
    for eq in observed
        # eq.lhs ∉ outpins && continue
        # if eq ∈ diffeqs
        #   move_observed_state_output()
        # else
        #   move_observed_intermediate_output()
        # end
    end
end


text:
output_clause(v) = join(v, " ") * " 0"

"""
Create string for DIMACS output.
Each clause is terminated by 0.

TODO: If clauses are too long they should be split up into separate lines.
"""
function dimacs_output(p::SATProblem)

    header = """
    c DIMACS created by sat.jl
    p cnf $(p.num_variables) $(length(p.clauses))
    """

    body = join(output_clause.(p.clauses), "\n")

    return header * body * "\n"
end

parse_line(s) = parse.(Int, split(s)[2:end])

# exit codes from CryptoMiniSat:
const SAT_status = Dict(10 => :sat,
                         15 => :unknown,
                         20 => :unsat)

"Use SAT solver to solve"
function call_solver(; solver="cryptominisat5", input="_input.cnf", output="_out.txt")

    # if solver == "cryptominisat5" && options == ""
    #     options = "--verb=0"  # remove verbose output
    # end

    # run the solver:
    cmd = `$solver $input`
    # @show cmd
    pipe = pipeline(cmd, stdout=output)
    out = run(pipe, wait=false);
    wait(out)

    status = SAT_status[out.exitcode]

    status ∈ (:unsat, :unknown) && return status, Int[]

    output = filter(line -> !startswith(line, "c"), readlines(output))

    # TODO: check for correct initial letter

    results = reduce(vcat, parse_line.(output[2:end]))
    pop!(results)  # remove final 0

    return status, results
end

function solve(p::SATProblem)
    s = dimacs_output(p)

    input = "_input.cnf"
    write(input, s)

    return @time call_solver(input=input)
end

ext: jl


text:
@enum(LintCodes,
MissingRef,
IncorrectCallArgs,
IncorrectIterSpec,
NothingEquality,
NothingNotEq,
ConstIfCondition,
EqInIfConditional,
PointlessOR,
PointlessAND,
UnusedBinding,
InvalidTypeDeclaration,
UnusedTypeParameter,
IncludeLoop,
MissingFile,
InvalidModuleName,
TypePiracy,
CannotDeclareConst,
InvalidRedefofConst)

const LintCodeDescriptions = Dict{LintCodes,String}(IncorrectCallArgs => "Possible method call error.",
    IncorrectIterSpec => "A loop iterator has been used that will likely error.",
    NothingEquality => "Compare against `nothing` using `===`",
    NothingNotEq => "Compare against `nothing` using `!==`",
    ConstIfCondition => "A boolean literal has been used as the conditional of an if statement - it will either always or never run.",
    EqInIfConditional => "Unbracketed assignment in if conditional statements is not allowed, did you mean to use ==?",
    PointlessOR => "The first argument of a `||` call is a boolean literal.",
    PointlessAND => "The first argument of a `&&` call is `false`.",
    UnusedBinding => "The variable name has been bound but not used.",
    InvalidTypeDeclaration => "A non-DataType has been used in a type declaration statement.",
    UnusedTypeParameter => "A DataType parameter has been specified but not used.",
    IncludeLoop => "Loop detected, this file has already been included.",
    MissingFile => "The included file can not be found.",
    InvalidModuleName => "Module name matches that of its parent.",
    TypePiracy => "An imported function has been extended without using module defined typed arguments.",
    CannotDeclareConst => "Cannot declare constant; it already has a value.",
    InvalidRedefofConst => "Invalid redefinition of constant.")

haserror(m::Meta) = m.error !== nothing
haserror(x::EXPR) = hasmeta(x) && haserror(x.meta)
errorof(x::EXPR) = hasmeta(x) ? x.meta.error : nothing
function seterror!(x::EXPR, e)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.error = e
end

function _typeof(x, state)
    if x isa EXPR
        if typof(x) in (CSTParser.Abstract, CSTParser.Primitive, CSTParser.Struct, CSTParser.Mutable)
            return CoreTypes.DataType
        elseif typof(x) in (CSTParser.ModuleH, CSTParser.BareModule)
            return CoreTypes.Module
        elseif CSTParser.defines_function(x)
            return CoreTypes.Function
        end
    elseif x isa SymbolServer.DataTypeStore
        return CoreTypes.DataType
    elseif x isa SymbolServer.FunctionStore
        return CoreTypes.Function
    end
end

# Call
function struct_nargs(x::EXPR)
    # struct defs wrapped in macros are likely to have some arbirtary additional constructors, so lets allow anything
    parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.MacroCall && return 0, typemax(Int), Symbol[], true
    minargs, maxargs, kws, kwsplat = 0, 0, Symbol[], false
    args = typof(x) === CSTParser.Mutable ? x[4] : x[3]
    length(args) == 0 && return 0, typemax(Int), kws, kwsplat
    inner_constructor = findfirst(a->CSTParser.defines_function(a), args.args)
    if inner_constructor !== nothing
        return func_nargs(args[inner_constructor])
    else
        minargs = maxargs = length(args)
    end
    return minargs, maxargs, kws, kwsplat
end

function func_nargs(x::EXPR)
    minargs, maxargs, kws, kwsplat = 0, 0, Symbol[], false
    sig = CSTParser.rem_where_decl(CSTParser.get_sig(x))
    for i = 2:length(sig)
        arg = sig[i]
        if typof(arg) === CSTParser.MacroCall && length(arg) > 1 &&
            _expr_assert(arg[1], MacroName, 2) && valof(arg[1][2]) == "nospecialize"
            if length(arg) == 2
                arg = arg[2]
            end
        end
        if typof(arg) === CSTParser.PUNCTUATION
            # skip
        elseif typof(arg) === CSTParser.Parameters
            for j = 1:length(arg)
                arg1 = arg[j]
                if typof(arg1) === CSTParser.Kw
                    push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg1[1]))))
                elseif typof(arg1) === CSTParser.BinaryOpCall && kindof(arg1[2]) === CSTParser.Tokens.DDDOT
                    kwsplat = true
                end
            end
        elseif typof(arg) === CSTParser.Kw
            if typof(arg[1]) === UnaryOpCall && kindof(arg[1][2]) === CSTParser.Tokens.DDDOT
                maxargs = typemax(Int)
            else
                maxargs !== typemax(Int) && (maxargs += 1)
            end
        elseif (typof(arg) === UnaryOpCall && kindof(arg[2]) === CSTParser.Tokens.DDDOT) ||
            (_binary_assert(arg, CSTParser.Tokens.DECLARATION) &&
            ((typof(arg[3]) === CSTParser.IDENTIFIER && valof(arg[3]) == "Vararg") ||
            (typof(arg[3]) === CSTParser.Curly && typof(arg[3][1]) === CSTParser.IDENTIFIER && valof(arg[3][1]) == "Vararg")))
            maxargs = typemax(Int)
        else
            minargs += 1
            maxargs !== typemax(Int) && (maxargs += 1)
        end
    end

    return minargs, maxargs, kws, kwsplat
end

function func_nargs(m::SymbolServer.MethodStore)
    minargs, maxargs, kws, kwsplat = 0, 0, Symbol[], false

    for arg in m.sig
        if CoreTypes.isva(last(arg))
            maxargs = typemax(Int)
        else
            minargs += 1
            maxargs !== typemax(Int) && (maxargs += 1)
        end
    end
    for kw in m.kws
        if endswith(String(kw), "...")
            kwsplat = true
        else
            push!(kws, kw)
        end
    end
    return minargs, maxargs, kws, kwsplat
end

function call_nargs(x::EXPR)
    minargs, maxargs, kws = 0, 0, Symbol[]
    if length(x) > 0
        for i = 2:length(x)
            arg = x[i]
            if typof(arg) === CSTParser.PUNCTUATION
                # skip
            elseif typof(x[i]) === CSTParser.Parameters
                for j = 1:length(x[i])
                    arg = x[i][j]
                    if typof(arg) === CSTParser.Kw
                        push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg[1]))))
                    end
                end
            elseif typof(arg) === CSTParser.Kw
                push!(kws, Symbol(CSTParser.str_value(CSTParser.get_arg_name(arg[1]))))
            elseif typof(arg) === CSTParser.UnaryOpCall && kindof(arg[2]) === CSTParser.Tokens.DDDOT
                maxargs = typemax(Int)
            else
                minargs += 1
                maxargs !== typemax(Int) && (maxargs += 1)
            end
        end
    else
        @info string("call_nargs: ", Expr(x))
    end

    return minargs, maxargs, kws
end

# compare_f_call(m_counts, call_counts) = true # fallback method

function compare_f_call(m_counts::Tuple{Int,Int,Array{Symbol},Bool}, call_counts::Tuple{Int,Int,Array{Symbol}})
    if call_counts[2] == typemax(Int)
        call_counts[1] <= call_counts[2] < m_counts[1] && return false
    else
        !(m_counts[1] <= call_counts[1] <= call_counts[2] <= m_counts[2]) && return false
    end
    if !m_counts[4] # no splatted kw in method sig
        length(call_counts[3]) > length(m_counts[3]) && return false # call has more kws than method accepts
        !all(kw in m_counts[3] for kw in call_counts[3]) && return false # call supplies a kw that isn't defined in the method
    else # splatted kw in method so accept any kw in call
    end
    return true
end

function check_call(x, server)
    if typof(x) === Call
        parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.Do && return # TODO: add number of args specified in do block.
        length(x) == 0 && return
        if isidentifier(first(x)) && hasref(first(x))
            func_ref = refof(first(x))
        elseif _binary_assert(x[1], CSTParser.Tokens.DOT) && typof(x[1]) === CSTParser.Quotenode && length(x[1][3]) > 0 && isidentifier(x[1][3][1]) && hasref(first(x)[3][1])
            func_ref = refof(first(last(first(x))))
        else
            return
        end
        if func_ref isa SymbolServer.FunctionStore || func_ref isa SymbolServer.DataTypeStore
            call_counts = call_nargs(x)
            function ff(m)
                m_counts = func_nargs(m)
                return compare_f_call(m_counts, call_counts)
            end
            tls = retrieve_toplevel_scope(x)
            tls == nothing && return
            iterate_over_ss_methods(func_ref, tls, server, ff) && return # returns if ff(m) == true for any methods
            seterror!(x, IncorrectCallArgs)
        elseif func_ref isa Binding && (func_ref.type === CoreTypes.Function || func_ref.type === CoreTypes.DataType)
            call_counts = call_nargs(x)
            function ff1(m)
                m_counts = func_nargs(m)
                return compare_f_call(m_counts, call_counts)
            end
            b = func_ref
            seen = Binding[]
            while b.next isa Binding && b.next.type == CoreTypes.Function && !(b in seen)
                push!(seen, b)
                b = b.next
            end
            seen = Binding[]
            while true
                b in seen && break
                if !(b isa Binding) # Needs to be cleaned up
                    if b isa SymbolServer.FunctionStore || b isa SymbolServer.DataTypeStore
                        tls = retrieve_toplevel_scope(x)
                        tls == nothing && return
                        iterate_over_ss_methods(b, tls, server, ff1) && return
                        break
                    else
                        return
                    end
                elseif b.type == CoreTypes.Function && b.val isa EXPR
                    m_counts = func_nargs(b.val)
                elseif b.type == CoreTypes.DataType && b.val isa EXPR && CSTParser.defines_struct(b.val)
                    m_counts = struct_nargs(b.val)
                elseif b.val isa SymbolServer.FunctionStore || b.val isa SymbolServer.DataTypeStore
                    tls = retrieve_toplevel_scope(x)
                    tls == nothing && return
                    iterate_over_ss_methods(b.val, tls, server, ff1) && return
                    break
                else
                    break
                end
                if compare_f_call(m_counts, call_counts)
                    return
                elseif CSTParser.rem_where_decl(CSTParser.get_sig(b.val)) == x
                    return
                end
                b.prev === nothing && break
                b in seen && break
                push!(seen, b)
                b = b.prev
            end
            seterror!(x, IncorrectCallArgs)
        end
    end
end

function check_loop_iter(x::EXPR, server)
    if typof(x) === CSTParser.For
        if length(x) > 1 && typof(x[2]) === CSTParser.BinaryOpCall && refof(x[2]) === nothing
            rng = x[2][3]
            if typof(rng) === CSTParser.LITERAL && kindof(rng) == CSTParser.Tokens.FLOAT || kindof(rng) == CSTParser.Tokens.INTEGER
                seterror!(x[2], IncorrectIterSpec)
            elseif typof(rng) === CSTParser.Call && refof(rng[1]) === getsymbolserver(server)[:Base][:length]
                seterror!(x[2], IncorrectIterSpec)
            end
        end
    elseif typof(x) === CSTParser.Generator
        for i = 3:length(x)
            if typof(x[i]) === CSTParser.BinaryOpCall && refof(x[i]) === nothing
                rng = x[i][3]
                if typof(rng) === CSTParser.LITERAL && kindof(rng) == CSTParser.Tokens.FLOAT || kindof(rng) == CSTParser.Tokens.INTEGER
                    seterror!(x[i], IncorrectIterSpec)
                elseif typof(rng) === CSTParser.Call && valof(rng[1]) == "length" && refof(rng[1]) === getsymbolserver(server)[:Base][:length]
                    seterror!(x[i], IncorrectIterSpec)
                end
            end
        end
    end
end

function check_nothing_equality(x::EXPR, server)
    if typof(x) === CSTParser.BinaryOpCall
        if kindof(x[2]) === CSTParser.Tokens.EQEQ && valof(x[3]) == "nothing" && refof(x[3]) === getsymbolserver(server)[:Core][:nothing]
            seterror!(x[2], NothingEquality)
        elseif kindof(x[2]) === CSTParser.Tokens.NOT_EQ && valof(x[3]) == "nothing" && refof(x[3]) === getsymbolserver(server)[:Core][:nothing]
            seterror!(x[2], NothingNotEq)
        end
    end
end

function _get_top_binding(x::EXPR, name::String)
    if scopeof(x) isa Scope
        return _get_top_binding(scopeof(x), name)
    elseif parentof(x) isa EXPR
        return _get_top_binding(parentof(x), name)
    else
        return nothing
    end
end
# Given the name of a function binding lookup the last declared method (i.e. the one
# stored in the scope)
_get_last_method(x, b::Binding, server) = b

function _get_last_method(x::EXPR, b::Binding, server)
    if scopeof(x) isa Scope
        if scopehasbinding(scopeof(x), b.name)
            if scopeof(x).names[b.name] isa Binding && scopeof(x).names[b.name].type == b.type
                return scopeof(x).names[b.name]
            end
        end
    elseif parentof(x) isa EXPR
        return _get_last_method(parentof(x), b.name, server)
    end
    return b
end

function _get_top_binding(s::Scope, name::String)
    if scopehasbinding(s, name)
        return s.names[name]
    elseif parentof(s) isa Scope
        return _get_top_binding(parentof(s), name)
    else
        return nothing
    end
end

function _get_global_scope(s::Scope)
    if !s.ismodule && parentof(s) isa Scope && parentof(s) != s
        return _get_global_scope(parentof(s))
    else
        return s
    end
end

function check_if_conds(x::EXPR)
    if typof(x) === CSTParser.If && length(x) > 2
        cond = typof(first(x)) == CSTParser.KEYWORD ? x[2] : x[1]
        if typof(cond) === CSTParser.LITERAL && (kindof(cond) === CSTParser.Tokens.TRUE || kindof(cond) === CSTParser.Tokens.FALSE)
            seterror!(cond, ConstIfCondition)
        elseif _binary_assert(cond, CSTParser.Tokens.EQ)
            seterror!(cond, EqInIfConditional)
        end
    end
end

function check_lazy(x::EXPR)
    if typof(x) === CSTParser.BinaryOpCall
        if kindof(x[2]) === CSTParser.Tokens.LAZY_OR
            if (typof(x[1]) === CSTParser.LITERAL && (kindof(x[1]) === CSTParser.Tokens.TRUE || kindof(x[1]) === CSTParser.Tokens.FALSE))
                seterror!(x, PointlessOR)
            end
        elseif kindof(x[2]) === CSTParser.Tokens.LAZY_AND
            if (typof(x[1]) === CSTParser.LITERAL && kindof(x[1]) === CSTParser.Tokens.FALSE) || (typof(x[3]) === CSTParser.LITERAL && kindof(x[3]) === CSTParser.Tokens.FALSE)
                seterror!(x, PointlessAND)
            end
        end
    end
end

function check_is_callable(x::EXPR, server)
    if typof(x) === CSTParser.Call
        func = first(x)

        if hasref(func)
        end
    end
end

function check_datatype_decl(x::EXPR, server)
    # Only call in function signatures
    if typof(x) === CSTParser.BinaryOpCall && length(x) === 3 && kindof(x[2]) === CSTParser.Tokens.DECLARATION &&
        parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.Call
        if hasref(last(x))
            dt = refof(last(x))
            if dt isa SymbolServer.DataTypeStore || dt isa Binding && dt.val isa SymbolServer.DataTypeStore
            elseif dt isa Binding && dt.type !== nothing
                safety_trip = 0
                while dt.type == CoreTypes.Function
                    safety_trip += 1
                    dt.prev === nothing && break
                    dt = dt.prev
                    if safety_trip > 50
                        @warn string(Expr(dt.name))
                        return
                    end
                end
                if dt.type !== nothing && dt.type !== CoreTypes.DataType
                    seterror!(x, InvalidTypeDeclaration)
                end
            end
        elseif typof(last(x)) === CSTParser.LITERAL
            seterror!(x, InvalidTypeDeclaration)
        end
    end
end

function check_modulename(x::EXPR)
    if (typof(x) === CSTParser.ModuleH || typof(x) === CSTParser.BareModule) && # x is a module
        scopeof(x) isa Scope && parentof(scopeof(x)) isa Scope && # it has a scope and a parent scope
        (typof(parentof(scopeof(x)).expr) === CSTParser.ModuleH ||
        typof(parentof(scopeof(x)).expr) === CSTParser.BareModule) && # the parent scope is a module
        valof(CSTParser.get_name(x)) == valof(CSTParser.get_name(parentof(scopeof(x)).expr)) # their names match
        seterror!(CSTParser.get_name(x), InvalidModuleName)
    end
end

mutable struct LintOptions
    call::Bool
    iter::Bool
    nothingcomp::Bool
    constif::Bool
    lazy::Bool
    datadecl::Bool
    typeparam::Bool
    modname::Bool
    pirates::Bool
end
LintOptions() = LintOptions(true, true, true, true, true, false, true, true, true)

function check_all(x::EXPR, opts::LintOptions, server)
    # Do checks
    opts.call && check_call(x, server)
    opts.iter && check_loop_iter(x, server)
    opts.nothingcomp && check_nothing_equality(x, server)
    opts.constif && check_if_conds(x)
    opts.lazy && check_lazy(x)
    # check_is_callable(x, server)
    opts.datadecl && check_datatype_decl(x, server)
    opts.typeparam && check_typeparams(x)
    opts.modname && check_modulename(x)
    opts.pirates && check_for_pirates(x)
    check_const_decl(x)
    check_const_redef(x)
    for i in 1:length(x)
        check_all(x[i], opts, server)
    end
end

function collect_hints(x::EXPR, missing = true, isquoted = false, errs = Tuple{Int,EXPR}[], pos = 0)
    if quoted(x)
        isquoted = true
    elseif isquoted && unquoted(x)
        isquoted = false
    end
    if typof(x) === CSTParser.ErrorToken
        # collect parse errors
        push!(errs, (pos, x))
    elseif !isquoted
        if missing && CSTParser.isidentifier(x) && !hasref(x) &&
            !(valof(x) == "var" && parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.NONSTDIDENTIFIER) &&
            !((valof(x) == "stdcall" || valof(x) == "cdecl" || valof(x) == "fastcall" || valof(x) == "thiscall" || valof(x) == "llvmcall") && is_in_fexpr(x, x->typof(x) === CSTParser.Call && isidentifier(x[1]) && valof(x[1]) == "ccall"))
            push!(errs, (pos, x))
        elseif haserror(x) && errorof(x) isa StaticLint.LintCodes
            # collect lint hints
            push!(errs, (pos, x))
        end
    end

    for i in 1:length(x)
        collect_hints(x[i], missing, isquoted, errs, pos)
        pos += x[i].fullspan
    end

    errs
end

function check_typeparams(x::EXPR)
    if typof(x) === CSTParser.WhereOpCall
        for i = 3:length(x)
            if hasbinding(x[i])
                if !(bindingof(x[i]).refs isa Vector)
                    seterror!(x[i], UnusedTypeParameter)
                elseif length(bindingof(x[i]).refs) == 1
                    # there should (will?) always be at least one reference in the where declaration
                    seterror!(x[i], UnusedTypeParameter)
                end
            end
        end
    end
end

function check_for_pirates(x::EXPR)
    if CSTParser.defines_function(x) && hasbinding(x) && overwrites_imported_function(bindingof(x))
        sig = CSTParser.rem_where_decl(CSTParser.get_sig(x))
        if typof(sig) == CSTParser.Call
            for i = 2:length(sig)
                if hasbinding(sig[i]) && bindingof(sig[i]).type isa Binding
                    return
                elseif refers_to_nonimported_type(sig[i])
                    return
                end
            end
            seterror!(x, TypePiracy)
        end
    end
end

function refers_to_nonimported_type(arg::EXPR)
    if hasref(arg) && refof(arg) isa Binding
        return true
    elseif typof(arg) === CSTParser.UnaryOpCall && length(arg) == 2 && (kindof(arg[1]) === CSTParser.Tokens.DECLARATION || kindof(arg[1]) === CSTParser.Tokens.ISSUBTYPE)
        return refers_to_nonimported_type(arg[2])
    elseif _binary_assert(arg, CSTParser.Tokens.DECLARATION)
        return refers_to_nonimported_type(arg[3])
    elseif typof(arg) === CSTParser.Curly
        for i = 1:length(arg)
            if refers_to_nonimported_type(arg[i])
                return true
            end
        end
        return false
    end
    return false
end

overwrites_imported_function(b, visited_bindings = Binding[]) = false
function overwrites_imported_function(b::Binding, visited_bindings = Binding[])
    if b in visited_bindings
        return false
    else
        push!(visited_bindings, b)
    end

    if b.prev isa Binding
        if b.prev.val isa EXPR
            # overwrites a within source bindig so lets check that
            return overwrites_imported_function(b.prev, visited_bindings)
        elseif (b.prev.val isa SymbolServer.FunctionStore || b.prev.val isa SymbolServer.DataTypeStore) && parentof(b.prev.name) isa EXPR && typof(parentof(b.prev.name)) === CSTParser.Import
            # explicitly imported, e.g. import ModuleName: somefunction
            return true
        end
    elseif b.prev isa SymbolServer.FunctionStore &&
        parentof(b.name) isa EXPR && typof(parentof(b.name)) === CSTParser.Quotenode && is_getfield(parentof(parentof(b.name)))
        fullname = parentof(parentof(b.name))
        # overwrites imported function by declaring full name, e.g. ModuleName.FunctionName
        if CSTParser.isidentifier(fullname[1]) && refof(fullname[1]) isa SymbolServer.ModuleStore
            return true
        elseif is_getfield(fullname[1]) && typof(fullname[1][3]) === CSTParser.Quotenode && refof(fullname[1][3][1]) isa SymbolServer.ModuleStore
            return true
        end
    end
    return false
end

function check_const_decl(x::EXPR)
    if CSTParser.defines_datatype(x) && hasbinding(x) && bindingof(x).prev !== nothing
        seterror!(x, CannotDeclareConst)
    end
end

function check_const_redef(x::EXPR)
    if hasbinding(x) && bindingof(x) isa Binding && bindingof(x).prev isa Binding && bindingof(x).prev.type === CoreTypes.DataType && bindingof(x).type !== CoreTypes.Function
        seterror!(x, InvalidRedefofConst)
    end
end


text:
#__precompile__()
module CircuitSimulation
import DifferentialEquations
import QuadGK
import LsqFit
import Plots
import ForwardDiff

Plots.pyplot()

# Package utilities
include("types.jl")

# FFT Utilities
include("ifft_times.jl")
include("fft_freqs.jl")
include("fft.jl")
include("plot.jl")

# Time-domain response utilities
include("td_simulation.jl")

# Network parameter utilities
include("convert_params.jl")

# Transmission line frequency-domain utilities
include("tline_diffeq_solve.jl")
include("abcd_tline.jl")
include("tline_z.jl")

# Transmission line fitting utilities
include("liu_fit.jl")
include("lossy_fit.jl")

end # module


text:
"""
    make_lindblad_operator(HK::Array{Float64,2},HS::Array{Float64,2},L_list,N::Int64=0)

# Argument:
- The real part and the imaginary part of the Hamiltonian: ``H = H_K-i H_S``
- A list of Lindblad terms ``L_k``'s
- N is number of states

# Output:
the operator for the vectorized system

- LK: real matrix corresponding to the real part of the Hamiltonian operator
- LS: real matrix corresponding to the imaginary part of the Hamiltonian operator
- LD: real matrix corresponding to the Lindblad operator
- Final system: ``\\rho_t = -i(L_K+i L_S)\\rho+L_D \\rho``
"""
function make_lindblad_operator(HK::Array{Float64,2},HS::Array{Float64,2},L_list,N::Int64=0)
    if(N==0)
        N = size(HK)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LK = (kron(It,HK) - kron(transpose(HK),It))
    LS = (kron(It,HS) - kron(transpose(HS),It))
    # Disspation part
    LD = zeros(Float64,N^2,N^2)
    for k = 1:length(L_list)
        LD += kron(L_list[k],L_list[k])-
              0.5*(kron(It,transpose(L_list[k])*L_list[k] )+kron(transpose(transpose(L_list[k])*L_list[k]),It))
    end
    return LK,LS,LD
end

function make_lindblad_operator(HK,HS,L_list,N::Int64=0)
    if(N==0)
        N = size(HK)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LK = (kron(It,HK) - kron(transpose(HK),It))
    LS = (kron(It,HS) - kron(transpose(HS),It))
    # Disspation part
    LD = kron(It,It)#zeros(Float64,N^2,N^2)
    for k = 1:length(L_list)
        LD += kron(L_list[k],L_list[k])-
              0.5*(kron(It,transpose(L_list[k])*L_list[k] )+kron(transpose(transpose(L_list[k])*L_list[k]),It))
    end
    return LK,LS,LD
end

"""
    make_lindblad_operator_complex(H::Array{ComplexF64,2},L_list,N::Int64=0)

# Argument:
- the Hamiltonain H and a list of Lindblad terms ``L_k``'s
- N is number of states

# Output:
the operator for the vectorized system

- LH: complex matrix corresponding to the Hamiltonian operator
- LD: complex matrix corresponding to the Lindblad operator
"""
function make_lindblad_operator_complex(H::Array{ComplexF64,2},L_list,N::Int64=0)
    if(N==0)
        N = size(H)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LH = -1.0im*(kron(It,H) - kron(transpose(H),It))
    # Disspation part
    LD = zeros(Float64,N^2,N^2)
    for k = 1:length(L_list)
        LD += kron(L_list[k],L_list[k])-
              0.5*(kron(It,transpose(L_list[k])*L_list[k] )+kron(transpose(transpose(L_list[k])*L_list[k]),It))
    end
    return LH,LD
end

"""
    make_Hamiltonian_operator(HK::Array{Float64,2},HS::Array{Float64,2},N::Int64=0)

# Argument:
- The real part and the imaginary part of the Hamiltonian: ``H = H_K-i H_S``
- N is number of states

# Output:
the operator for the vectorized system

- LK: real matrix corresponding to the real part of the Hamiltonian operator
- LS: real matrix corresponding to the imaginary part of the Hamiltonian operator
"""
function make_hamiltonian_operator(HK::Array{Float64,2},HS::Array{Float64,2},N::Int64=0)
    if(N==0)
        N = size(HK)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LK = (kron(It,HK) - kron(transpose(HK),It))
    LS = (kron(It,HS) - kron(transpose(HS),It))
    return LK,LS
end

function make_hamiltonian_operator(HK,HS,N::Int64=0)
    if(N==0)
        N = size(HK)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LK = (kron(It,HK) - kron(transpose(HK),It))
    LS = (kron(It,HS) - kron(transpose(HS),It))
    return LK,LS
end

"""
    make_Hamiltonian_operator(HK::Array{Float64,2},N::Int64=0)

# Argument:
- The real part and the imaginary part of the Hamiltonian: ``H = H_K``
- N is number of states

# Output:
the operator for the vectorized system

- LK: real matrix corresponding to the real part of the Hamiltonian operator
"""
function make_hamiltonian_operator(HK::Array{Float64,2},N::Int64=0)
    if(N==0)
        N = size(HK)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LK = (kron(It,HK) - kron(transpose(HK),It))
    return LK
end

function make_hamiltonian_operator(HK,N::Int64=0)
    if(N==0)
        N = size(HK)[1]
    end
    # define an identity operator
    It = Array{Float64, 2}(I, N, N)
    # Hamiltonian part
    LK = (kron(It,HK) - kron(transpose(HK),It))
    return LK
end


text:
function triangulate2D(m::CAGD.Model)
    return Lar.triangulate2D(convert(Lar.Points, m.G'), [m.T[1], m.T[2]])
end


text:
function local_metric_jacobian(M::MetricManifold, x)
    n = size(x, 1)
    ∂g = reshape(ForwardDiff.jacobian(x -> local_metric(M, x), x), n, n, n)
    return ∂g
end

function christoffel_symbols_second_jacobian(M::MetricManifold, x)
    n = size(x, 1)
    ∂Γ = reshape(
        ForwardDiff.jacobian(x -> christoffel_symbols_second(M, x), x),
        n, n, n, n
    )
    return ∂Γ
end


text:
# Autogenerated wrapper script for gsasl_jll for x86_64-w64-mingw32
export gsasl, libgsasl

JLLWrappers.@generate_wrapper_header("gsasl")
JLLWrappers.@declare_executable_product(gsasl)
JLLWrappers.@declare_library_product(libgsasl, "libgsasl-18.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        gsasl,
        "bin\\gsasl.exe",
    )

    JLLWrappers.@init_library_product(
        libgsasl,
        "bin\\libgsasl-18.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Autogenerated wrapper script for assimp_jll for aarch64-linux-musl-cxx03
export libassimp

using Zlib_jll
JLLWrappers.@generate_wrapper_header("assimp")
JLLWrappers.@declare_library_product(libassimp, "libassimp.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libassimp,
        "lib/libassimp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Copyright 2022 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    pressure(A :: Float64, A0 :: Float64, beta :: Float64, Pext :: Float64)

Compute pressure at a specific node by means of the linear-elastic constitutive
equation.
"""
function pressure(A :: Float64, A0 :: Float64, beta :: Float64, Pext :: Float64)
    return Pext + beta*(sqrt(A/A0) - 1.0)
end

# version with pre-computed `sqrt(A/A0)`
function pressure(s_A_over_A0 :: Float64, beta :: Float64, Pext :: Float64)
    return Pext + beta*(s_A_over_A0 - 1.0)
end

"""
    waveSpeed(A :: Float64, gamma :: Float64)

Compute pulse wave velocity at the given node.
"""
function waveSpeed(A :: Float64, gamma :: Float64)
    return sqrt(3*gamma*sqrt(A)*0.5)
end

function waveSpeedSA(sA :: Float64, gamma :: Float64)
    return sqrt(1.5*gamma*sA)
end

"""
    calculateDeltaT(vessels :: Array{Vessel,1}, Ccfl :: Float64)

Compute global minimum Δt through CFL condition. For each vessel the Δt is computed as

Delta t = C_{CFL} frac{Delta x}{S_{max}}

where Smax is the maximum between the forward and the backward characteristics, and Ccfl is the Courant-Friedrichs-Lewy condition defined by the user.
"""
function calculateDeltaT(vessels, Ccfl :: Float64)
    dt = 1.0
    for v in vessels
        Smax = 0.0
        for j in 1:v.M
            lambda = abs(v.u[j] + v.c[j])
            Smax < lambda && (Smax = lambda)
        end
        vessel_dt = v.dx*Ccfl/Smax
        dt > vessel_dt && (dt = vessel_dt)
    end
    return dt
end

"""
    solveModel(vessel :: Array){Vessel,1}, edges :: Array{Int,2}, blood :: Blood,
               dt :: Float64, current_time :: Float64)

Run the solver on each vessel one-by-one as listed in the .yml file.
"""
function solveModel(vessels :: Array{Vessel,1}, edges :: Array{Int,2}, blood :: Blood,
                    dt :: Float64, current_time :: Float64)
    @inbounds for j in 1:size(edges)[1]
        @inbounds i = edges[j,1]
        @inbounds v = vessels[i]

        solveVessel(v, blood, dt, current_time)
        solveOutlet(j, v, blood, vessels, edges, dt)
    end
end

"""
    solveVessel(vessel :: Vessel, blood :: Blood, dt :: Float64, current_time :: Float64)

Apply inlet boundary condition (if any) and run MUSCL solver along the vessel.
"""
function solveVessel(vessel :: Vessel, blood :: Blood, dt :: Float64,
                     current_time :: Float64)
    if vessel.inlet
        setInletBC(current_time, dt, vessel)
    end

    muscl(vessel, dt, blood)
end

"""
    solveOutlet(j :: Int, vessel :: Vessel, blood :: Blood, vessels :: Array{Vessel,1},
                edges :: Array{Int,2})

Join vessel in case of junction or apply outlet boundary condition otherwise.
"""
function solveOutlet(j :: Int, vessel :: Vessel, blood :: Blood, vessels :: Array{Vessel,1},
                     edges :: Array{Int,2}, dt :: Float64)
    @inbounds i = edges[j,1]
    @inbounds t = edges[j,3]

    if vessel.outlet != "none"
        setOutletBC(dt, vessel)

    elseif size(findall(edges[:,2] .== t))[1] == 2
        d1_i = findall(edges[:,2] .== t)[1]
        d2_i = findall(edges[:,2] .== t)[2]
        joinVessels(blood, vessel, vessels[d1_i], vessels[d2_i])

    elseif size(findall(edges[:,3] .== t))[1] == 1
        d_i = findall(edges[:,2] .== t)[1]
        joinVessels(blood, vessel, vessels[d_i])

    elseif size(findall(edges[:,3] .== t))[1] == 2
        p1_i = findall(edges[:,3] .== t)[1]
        p2_i = findall(edges[:,3] .== t)[2]

        if maximum([p1_i, p2_i]) == i
            p2_i = minimum([p1_i, p2_i])
            d = findall(edges[:,2] .== t)[1]
            solveAnastomosis(vessel, vessels[p2_i], vessels[d])
        end
    end
end

"""
    muscl(v :: Vessel, dt :: Float64, b :: Blood)

Run MUSCL solver along the vessel.
"""
function muscl(v :: Vessel, dt :: Float64, b :: Blood)
    @inbounds v.vA[1] = v.U00A
    @inbounds v.vA[end] = v.UM1A

    @inbounds v.vQ[1] = v.U00Q
    @inbounds v.vQ[end] = v.UM1Q

    @fastmath @inbounds @simd for i = 2:v.M+1
        v.vA[i] = v.A[i-1]
        v.vQ[i] = v.Q[i-1]
    end

    v.slopesA = computeLimiter(v, v.vA, v.invDx, v.dU, v.slopesA)
    v.slopesQ = computeLimiter(v, v.vQ, v.invDx, v.dU, v.slopesQ)

    @fastmath @inbounds @simd for i = 1:v.M+2
        slopeA_halfDx = v.slopesA[i]*v.halfDx
        v.Al[i] = v.vA[i] + slopeA_halfDx
        v.Ar[i] = v.vA[i] - slopeA_halfDx

        slopeQ_halfDX = v.slopesQ[i]*v.halfDx
        v.Ql[i] = v.vQ[i] + slopeQ_halfDX
        v.Qr[i] = v.vQ[i] - slopeQ_halfDX
    end

    v.Fl = computeFlux(v, v.Al, v.Ql, v.Fl)
    v.Fr = computeFlux(v, v.Ar, v.Qr, v.Fr)

    @fastmath dxDt = v.dx/dt
    @fastmath invDxDt = 1.0/dxDt

    @fastmath @inbounds @simd for i = 1:v.M+1
        v.flux[1,i] = 0.5*(v.Fr[1,i+1] + v.Fl[1,i] - dxDt*(v.Ar[i+1] - v.Al[i]))
        v.flux[2,i] = 0.5*(v.Fr[2,i+1] + v.Fl[2,i] - dxDt*(v.Qr[i+1] - v.Ql[i]))
    end

    @fastmath @inbounds @simd for i = 2:v.M+1
        v.uStar[1,i] = v.vA[i] + invDxDt*(v.flux[1,i-1] - v.flux[1,i])
        v.uStar[2,i] = v.vQ[i] + invDxDt*(v.flux[2,i-1] - v.flux[2,i])
    end

    @inbounds v.uStar[1,1] = v.uStar[1,2]
    @inbounds v.uStar[2,1] = v.uStar[2,2]
    @inbounds v.uStar[1,end] = v.uStar[1,end-1]
    @inbounds v.uStar[2,end] = v.uStar[2,end-1]

    v.slopesA = computeLimiter(v, v.uStar, 1, v.invDx, v.dU, v.slopesA)
    v.slopesQ = computeLimiter(v, v.uStar, 2, v.invDx, v.dU, v.slopesQ)

    @fastmath @inbounds @simd for i = 1:v.M+2
        v.Al[i] = v.uStar[1,i] + v.slopesA[i]*v.halfDx
        v.Ar[i] = v.uStar[1,i] - v.slopesA[i]*v.halfDx

        v.Ql[i] = v.uStar[2,i] + v.slopesQ[i]*v.halfDx
        v.Qr[i] = v.uStar[2,i] - v.slopesQ[i]*v.halfDx
    end

    v.Fl = computeFlux(v, v.Al, v.Ql, v.Fl)
    v.Fr = computeFlux(v, v.Ar, v.Qr, v.Fr)

    @fastmath @inbounds @simd for i = 1:v.M+1
        v.flux[1,i] = 0.5*(v.Fr[1,i+1] + v.Fl[1,i] - dxDt*(v.Ar[i+1] - v.Al[i]))
        v.flux[2,i] = 0.5*(v.Fr[2,i+1] + v.Fl[2,i] - dxDt*(v.Qr[i+1] - v.Ql[i]))
    end

    @fastmath @inbounds @simd for i = 2:v.M+1
        v.A[i-1] = 0.5*(v.A[i-1] + v.uStar[1,i] + invDxDt*(v.flux[1,i-1] - v.flux[1,i]))
        v.Q[i-1] = 0.5*(v.Q[i-1] + v.uStar[2,i] + invDxDt*(v.flux[2,i-1] - v.flux[2,i]))
    end

    #source term
    @fastmath @inbounds @simd for i = 1:v.M
        s_A = sqrt(v.A[i])
        Si = - v.viscT*v.Q[i]/v.A[i] - v.wallE[i]*(s_A - v.s_A0[i])*v.A[i]
        v.Q[i] += dt*Si

        v.P[i] = pressure(s_A*v.s_inv_A0[i], v.beta[i], v.Pext)
        v.c[i] = waveSpeedSA(s_A, v.gamma[i])

    end

    #parabolic system (viscoelastic part)
    if v.wallVa[1] != 0.0
        Td = 1.0/dt + v.wallVb
        Tlu = -v.wallVa
        T = Tridiagonal(Tlu[1:end-1], Td, Tlu[2:end])

        d = (1.0/dt - v.wallVb).*v.Q
        d[1] += v.wallVa[2]*v.Q[2]
        for i = 2:v.M-1
            d[i] += v.wallVa[i+1]*v.Q[i+1] + v.wallVa[i-1]*v.Q[i-1]
        end
        d[end] += v.wallVa[end-1]*v.Q[end-1]

        v.Q = T\d
    end

    @fastmath @inbounds @simd for i = 1:v.M
        v.u[i] = v.Q[i]/v.A[i]
    end
end

"""
    computeFlux(v :: Vessel, A :: Array{Float64,1}, Q :: Array{Float64,1},
                Flux :: Array{Float64,2})
"""
function computeFlux(v :: Vessel, A :: Array{Float64,1}, Q :: Array{Float64,1},
                     Flux :: Array{Float64,2})
    @fastmath @inbounds @simd for i in 1:v.M+2
        Flux[1,i] = Q[i]
        Flux[2,i] = Q[i]*Q[i]/A[i] + v.gamma_ghost[i]*A[i]*sqrt(A[i])
    end

    return Flux
end

"""
    maxMod(a :: Float64, b :: Float64)
"""
function maxMod(a :: Float64, b :: Float64)
    if a > b
        return a
    else
        return b
    end
end

"""
    minMod(a :: Float64, b :: Float64)
"""
function minMod(a :: Float64, b :: Float64)
    if a <= 0.0 || b <= 0.0
        return 0.0
    elseif a < b
        return a
    else
        return b
    end
end

"""
    superBee(v :: Vessel, dU :: Array{Float64,2}, slopes :: Array{Float64,1})
"""
function superBee(v :: Vessel, dU :: Array{Float64,2}, slopes :: Array{Float64,1})
    @fastmath @inbounds @simd for i in 1:v.M+2
        s1 = minMod(dU[1,i], 2*dU[2,i])
        s2 = minMod(2*dU[1,i], dU[2,i])
        slopes[i] = maxMod(s1, s2)
    end

    return slopes
end

"""
    computeLimiter(v :: Vessel, U :: Array{Float64,1}, invDx :: Float64,
                   dU :: Array{Float64,2}, slopes :: Array{Float64,1})
"""
function computeLimiter(v :: Vessel, U :: Array{Float64,1}, invDx :: Float64,
                        dU :: Array{Float64,2}, slopes :: Array{Float64,1})
    @fastmath @inbounds @simd for i = 2:v.M+2
        dU[1,i]   = (U[i] - U[i-1])*v.invDx
        dU[2,i-1] = dU[1,i]
    end

    return superBee(v, dU, slopes)
end

function computeLimiter(v :: Vessel, U :: Array{Float64,2}, idx :: Int,
                        invDx :: Float64, dU :: Array{Float64,2},
                        slopes :: Array{Float64,1})
    @inbounds U = U[idx,:]
    @fastmath @inbounds @simd for i = 2:v.M+2
        dU[1,i]   = (U[i] - U[i-1])*invDx
        dU[2,i-1] = dU[1,i]
    end

    return superBee(v, dU, slopes)
end


text:
using Test, Random, TensorNetworkEvolve, Graphs, Yao
using OMEinsumContractionOrders: TreeSA
using LinearAlgebra

@testset "initial state" begin
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end
    peps = zero_vidalpeps(ComplexF64, g, 2, optimizer=TreeSA(nslices=3))
    @test virtualbonds(peps) == [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
    @test vec(peps) ≈ (x=zeros(ComplexF64, 1<<5); x[1]=1; x)
    @test newlabel(peps, 2) == 11+2
    @test newlabel(peps, 4) == 11+4
    @test getphysicallabel(peps, 2) == 2
    @test length(getvlabel(peps, 2)) == 4
    apply_onsite!(peps, 1, [0 1; 1 0])
    @test vec(peps) ≈ (x=zeros(ComplexF64, 1<<5); x[2]=1; x)
    apply_onbond!(peps, 1, 2, reshape(Matrix(ConstGate.CNOT), 2, 2, 2, 2))
    @test vec(peps) ≈ (x=zeros(ComplexF64, 1<<5); x[4]=1; x)

    peps = rand_vidalpeps(ComplexF64, g, 2)
    normalize!(peps)
    @test norm(peps) ≈ 1
end

@testset "random gate" begin
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end
    peps = rand_vidalpeps(ComplexF64, g, 2; Dmax=4)
    reg = Yao.ArrayReg(vec(peps))

    # unary
    m = rand_unitary(2)
    apply_onsite!(peps, 1, m)
    reg |> put(5, (1,)=>matblock(m))
    @test vec(peps) ≈ statevec(reg)

    m = rand_unitary(4)
    reg |> put(5, (4,2)=>matblock(m))
    apply_onbond!(peps, 4, 2, reshape(m,2,2,2,2))
    @test vec(peps) ≈ statevec(reg)
end

@testset "conj, variables, and load" begin
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end
    peps = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    @test statevec(conj(peps)) ≈ conj(statevec(peps))
    P = deepcopy(peps)
    vars = variables(peps)
    randn!(peps)
    @test !(statevec(P) ≈ statevec(peps))
    load_variables!(peps, vars)
    @test statevec(P) ≈ statevec(peps)
    peps2 = randn!(peps)
    peps2 = load_variables(peps2, vars)
    @test statevec(P) ≈ statevec(peps2)
end

@testset "inner product and normalize" begin
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end

    p1 = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    @test inner_product(p1, p1)[] ≈ statevec(p1)' * statevec(p1)

    p1 = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    p2 = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    @test inner_product(p1, p2)[] ≈ statevec(p1)' * statevec(p2)

    p1 = rand_vidalpeps(ComplexF64, g, 2; Dmax=4)
    p2 = rand_vidalpeps(ComplexF64, g, 2; Dmax=4)
    @test inner_product(p1, p2)[] ≈ statevec(p1)' * statevec(p2)

    peps = rand_vidalpeps(ComplexF64, g, 2; Dmax=4)
    @test norm(vec(peps)) ≈ norm(peps)
    normalize!(peps)
    @test norm(peps) ≈ 1
end

@testset "expectation values" begin
    function rand_hamiltonian(g::SimpleGraph)
        blocks = []
        nbit = nv(g)
        for edge in edges(g)
            i,j = edge.src, edge.dst
            push!(blocks, put(nbit,(i,j)=>rand()*kron(X,X)) + rand()*kron(nbit, i=>Y,j=>Y) + control(nbit, i,j=>rand()*Z))
        end
        for i in vertices(g)
            push!(blocks, put(nbit,(i,)=>rand()*X))
        end
        return +(blocks...)
    end
    g = SimpleGraph(5)
    for (i,j) in [(1,2), (1,3), (2,4), (2,5), (3,4), (3,5)]
        add_edge!(g, i, j)
    end

    h = rand_hamiltonian(g)
    p1 = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    p2 = rand_simplepeps(ComplexF64, g, 2; Dmax=4)
    @test expect(h, p1, p2)[] ≈ statevec(p1)' * mat(h) * statevec(p2)
    @test expect(h, p1, p1)[] ≈ statevec(p1)' * mat(h) * statevec(p1)
end


text:
export poisson_trunc_point, poisson_cum_rtp

function poisson_trunc_point(λ, tol)
    k = 0.0
    expo = -λ
    term = exp(expo)
    term_sum = term
    log_λ = log(λ)
    probs = Vector{Float64}()

    while term_sum < tol / 2.0
        k += 1.0
        expo += log_λ - log(k)
        term = exp(expo)
        term_sum += term
    end

    ltp = Int(k)
    push!(probs, term)

    while 1.0 - term_sum > tol / 2.0
        k += 1.0
        expo += log_λ - log(k)
        term = exp(expo)
        term_sum += term
        push!(probs, term)
    end

    rtp = Int(k)
    @assert length(probs) == (rtp - ltp + 1)

    return ltp, rtp, probs
end

function poisson_cum_rtp(λ, t, tol)
    k = 0.0
    log_λ = log(λ)
    expo = -λ
    term = exp(expo)
    term_sum = term
    probs = Vector{Float64}()

    push!(probs, term)

    while tol < t * (1.0 - term_sum)
        k += 1.0
        expo += log_λ - log(k)
        term = exp(expo)
        term_sum += term
        push!(probs, term)
    end
    rtp = Int(k)
    @assert rtp + 1 == length(probs)
    return rtp, probs
end


text:
@testset "kernelpdmat" begin
    rng = MersenneTwister(123456)
    A = rand(rng, 10, 5)
    k = SqExponentialKernel()
    for obsdim in [1,2]
        @test all(Matrix(kernelpdmat(k,A,obsdim=obsdim)) .≈ Matrix(PDMat(kernelmatrix(k,A,obsdim=obsdim))))
        # @test_throws ErrorException kernelpdmat(k,ones(100,100),obsdim=obsdim)
    end
end