text:
"""
    mutable struct AbstractGlobal
        t::Any     # analyzed type
        iscd::Bool # whether this abstract global variable is declarared as constant or not
    end

Wraps a global variable whose type is analyzed by abtract interpretation.
`AbstractGlobal` object will be actually evaluated into the context module, and a later
  analysis may refer to its type or alter it on another assignment.

!!! note
    The type of the wrapped global variable will be propagated only when in a toplevel frame,
      and thus we don't care about the analysis cache invalidation on a refinement of the
      wrapped global variable, since JET doesn't cache the toplevel frame.
"""
mutable struct AbstractGlobal
    # analyzed type
    t::Any
    iscd::Bool

    function AbstractGlobal(@nospecialize(t),
                            iscd::Bool,
                            )
        return new(t,
                   iscd,
                   )
    end
end

@reportdef struct NoMethodErrorReport <: InferenceErrorReport
    @nospecialize(t::Union{Type,Vector{Type}})
end
# TODO count invalid union split case
get_msg(::Type{NoMethodErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, @nospecialize(t::Type)) =
    "no matching method found for call signature ($t)"
get_msg(::Type{NoMethodErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, ts::Vector{Type}) =
    "for $(length(ts)) of union split cases, no matching method found for call signatures ($(join(ts, ", "))))"

function is_empty_match(info::MethodMatchInfo)
    res = info.results
    isa(res, MethodLookupResult) || return false # when does this happen ?
    return isempty(res.matches)
end

# branching on https://github.com/JuliaLang/julia/pull/41020
@static if isdefined(CC, :MethodCallResult)

import .CC:
    bail_out_toplevel_call,
    bail_out_call,
    add_call_backedges!,
    MethodCallResult,
    ConstCallInfo,
    const_prop_entry_heuristic,
    const_prop_rettype_heuristic

# TODO:
# - report "too many method matched"
# - maybe "cound not identify method table for call" won't happen since we eagerly propagate bottom for e.g. undef var case, etc.

function CC.abstract_call_gf_by_type(analyzer::AbstractAnalyzer, @nospecialize(f),
                                     fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any}, @nospecialize(atype),
                                     sv::InferenceState, max_methods::Int = InferenceParams(analyzer).MAX_METHODS)
    ret = @invoke abstract_call_gf_by_type(analyzer::AbstractInterpreter, @nospecialize(f),
                                           fargs::Union{Nothing,Vector{Any}}, argtypes::Vector{Any}, @nospecialize(atype),
                                           sv::InferenceState, max_methods::Int)

    info = ret.info
    if isa(info, ConstCallInfo)
        info = info.call # unwrap to `MethodMatchInfo` or `UnionSplitInfo`
    end
    # report passes for no matching methods error
    if isa(info, UnionSplitInfo)
        ReportPass(analyzer)(NoMethodErrorReport, analyzer, sv, info, argtypes)
    elseif isa(info, MethodMatchInfo)
        ReportPass(analyzer)(NoMethodErrorReport, analyzer, sv, info, atype)
    end

    analyze_task_parallel_code!(analyzer, f, argtypes, sv)

    return ret
end

function (::SoundBasicPass)(::Type{NoMethodErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, info::UnionSplitInfo, argtypes::Vector{Any})
    # check each match for union-split signature
    split_argtypes = nothing
    ts = nothing

    for (i, matchinfo) in enumerate(info.matches)
        if is_empty_match(matchinfo)
            isnothing(split_argtypes) && (split_argtypes = switchtupleunion(argtypes))
            isnothing(ts) && (ts = Type[])
            sig_n = argtypes_to_type(split_argtypes[i])
            push!(ts, sig_n)
        end
    end

    if !isnothing(ts)
        add_new_report!(NoMethodErrorReport(analyzer, sv, ts), analyzer)
    end
end
function (::SoundBasicPass)(::Type{NoMethodErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, info::MethodMatchInfo, @nospecialize(atype))
    if is_empty_match(info)
        add_new_report!(NoMethodErrorReport(analyzer, sv, atype), analyzer)
    end
end

@doc """
    bail_out_toplevel_call(analyzer::AbstractAnalyzer, ...)

An overload for `abstract_call_gf_by_type(analyzer::AbstractAnalyzer, ...)`, which keeps
  inference on non-concrete call sites in a toplevel frame created by
  [`virtual_process`](@ref).
"""
function CC.bail_out_toplevel_call(analyzer::AbstractAnalyzer, @nospecialize(sig), sv)
    return isa(sv.linfo.def, Module) && !isdispatchtuple(sig) && !istoplevel(analyzer, sv)
end

@doc """
    bail_out_call(analyzer::AbstractAnalyzer, ...)

With this overload, `abstract_call_gf_by_type(analyzer::AbstractAnalyzer, ...)` doesn't bail out
  inference even after the current return type grows up to `Any` and collects as much error
  points as possible.
Of course this slows down inference performance, but hoopefully it stays to be "practical"
  speed since the number of matching methods are limited beforehand.
"""
CC.bail_out_call(analyzer::AbstractAnalyzer, @nospecialize(t), sv) = false

@doc """
    add_call_backedges!(analyzer::AbstractAnalyzer, ...)

An overload for `abstract_call_gf_by_type(analyzer::AbstractAnalyzer, ...)`, which always add
  backedges (even if a new method can't refine the return type grew up to `Any`).
This is because a new method definition always has a potential to change the JET analysis result.
"""
:(add_call_backedges!)
# branch on https://github.com/JuliaLang/julia/pull/41633
@static if isdefined(CC, :find_matching_methods)
import .CC:
    MethodMatches,
    UnionSplitMethodMatches
function CC.add_call_backedges!(analyzer::AbstractAnalyzer, @nospecialize(rettype), edges::Vector{MethodInstance},
                                matches::Union{MethodMatches,UnionSplitMethodMatches}, @nospecialize(atype),
                                sv::InferenceState)
    # NOTE a new method may refine analysis, so we always add backedges
    # # for `NativeInterpreter`, we don't add backedges when a new method couldn't refine (widen) this type
    # rettype === Any && return
    # for edge in edges
    #     add_backedge!(edge, sv)
    # end
    # also need an edge to the method table in case something gets
    # added that did not intersect with any existing method
    if isa(matches, MethodMatches)
        matches.fullmatch || add_mt_backedge!(matches.mt, atype, sv)
    else
        for (thisfullmatch, mt) in zip(matches.fullmatches, matches.mts)
            thisfullmatch || add_mt_backedge!(mt, atype, sv)
        end
    end
end
else # @static if isdefined(CC, :find_matching_methods)
function CC.add_call_backedges!(analyzer::AbstractAnalyzer, @nospecialize(rettype), edges::Vector{MethodInstance},
                                fullmatch::Vector{Bool}, mts::Vector{Core.MethodTable}, @nospecialize(atype),
                                sv::InferenceState)
    # NOTE a new method may refine analysis, so we always add backedges
    # if rettype === Any
    #     # for `NativeInterpreter`, we don't add backedges when a new method couldn't refine
    #     # (widen) this type
    #     return
    # end
    for edge in edges
        add_backedge!(edge, sv)
    end
    for (thisfullmatch, mt) in zip(fullmatch, mts)
        if !thisfullmatch
            # also need an edge to the method table in case something gets
            # added that did not intersect with any existing method
            add_mt_backedge!(mt, atype, sv)
        end
    end
end
end # @static if isdefined(CC, :find_matching_methods)

function CC.abstract_call_method_with_const_args(analyzer::AbstractAnalyzer, result::MethodCallResult,
                                                 @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch,
                                                 sv::InferenceState, va_override::Bool)
    const_result =
        @invoke abstract_call_method_with_const_args(analyzer::AbstractInterpreter, result::MethodCallResult,
                                                     @nospecialize(f), argtypes::Vector{Any}, match::MethodMatch,
                                                     sv::InferenceState, va_override::Bool)
    # update reports if constant prop' was successful
    # branch on https://github.com/JuliaLang/julia/pull/41697/
    @static if VERSION ≥ v"1.8.0-DEV.282"
        if !isnothing(const_result)
            # successful constant prop', we also need to update reports
            update_reports!(analyzer, sv)
        end
    else
        if !isnothing(getfield(const_result, 2))
            # successful constant prop', we also need to update reports
            update_reports!(analyzer, sv)
        end
    end
    return const_result
end

@doc """
    const_prop_entry_heuristic(analyzer::AbstractAnalyzer, result::MethodCallResult, sv::InferenceState)

This overload for `abstract_call_method_with_const_args(analyzer::AbstractAnalyzer, ...)` forces
constant prop' even if an inference result can't be improved anymore _with respect to the
return type_, e.g. when `result.rt` is already `Const`.
Especially, this overload implements an heuristic to force constant prop' when any error points
have been reported while the previous abstract method call without constant arguments.
The reason we want much more aggressive constant propagation by that heuristic is that it's
highly possible constant prop' can produce more accurate analysis result, by throwing away
false positive error reports by cutting off the unreachable control flow or detecting
must-reachable `throw` calls.
"""
function CC.const_prop_entry_heuristic(analyzer::AbstractAnalyzer, result::MethodCallResult, sv::InferenceState)
    edge = result.edge # edge associated with the previous non-constant analysis
    if !isnothing(edge)
        # if any error has been reported within the previous `abstract_call_method` (associated with `edge`),
        # force constant prop' and hope it can cut-off false positives
        any(is_from_same_frame(sv.linfo, edge), get_reports(analyzer)) && return true
    end
    return @invoke const_prop_entry_heuristic(analyzer::AbstractInterpreter, result::MethodCallResult, sv::InferenceState)
end

else # @static if isdefined(CC, :MethodCallResult)

include("legacy/abstractinterpretation")

end # @static if isdefined(CC, :MethodCallResult)

"""
    analyze_task_parallel_code!(analyzer::AbstractAnalyzer, @nospecialize(f), argtypes::Vector{Any}, sv::InferenceState)

Adds special cased analysis pass for task parallelism (xref: <https://github.com/aviatesk/JET.jl/issues/114>).
In Julia's task parallelism implementation, parallel code is represented as closure and it's
  wrapped in a `Task` object. `NativeInterpreter` doesn't run type inference nor optimization
  on the body of those closures when compiling code that creates parallel tasks, but JET will
  try to run additional analysis pass by recurring into the closures.

!!! note
    JET won't do anything other than doing JET analysis, e.g. won't annotate return type
    of wrapped code block in order to not confuse the original `AbstractInterpreter` routine
    track <https://github.com/JuliaLang/julia/pull/39773> for the changes in native abstract
    interpretation routine.
"""
function analyze_task_parallel_code!(analyzer::AbstractAnalyzer, @nospecialize(f), argtypes::Vector{Any}, sv::InferenceState)
    # TODO ideally JET should analyze a closure wrapped in a `Task` only when it's `schedule`d
    # but the `Task` construction may not happen in the same frame where it's `schedule`d
    # and so we may not be able to access to the closure at that point
    # as a compromise, JET now invokes the additional analysis on `Task` construction,
    # regardless of whether it's really `schedule`d or not
    if f === Task &&
       length(argtypes) ≥ 2 &&
       (v = argtypes[2]; v ⊑ Function)
        # if we encounter `Task(::Function)`, try to get its inner function and run analysis on it
        # the closure can be a nullary lambda that really doesn't depend on
        # the captured environment, and in that case we can retrieve it as
        # a function object, otherwise we will try to retrieve the type of the closure
        ft = (isa(v, Const) ? Core.Typeof(v.val) :
              isa(v, Core.PartialStruct) ? v.typ :
              isa(v, DataType) ? v :
              return)::Type
        analyze_additional_pass_by_type!(analyzer, Tuple{ft}, sv)
        return
    end
    return
end

# run additional interpretation with a new interpreter,
# and then append the reports to the original interpreter
function analyze_additional_pass_by_type!(analyzer::AbstractAnalyzer, @nospecialize(tt::Type{<:Tuple}), sv::InferenceState)
    newanalyzer = AbstractAnalyzer(analyzer)

    # in order to preserve the inference termination, we keep to use the current frame
    # and borrow the `AbstractInterpreter`'s cycle detection logic
    # XXX the additional analysis pass by `abstract_call_method` may involve various site-effects,
    # but what we're doing here is essentially equivalent to modifying the user code and inlining
    # the threaded code block as a usual code block, and thus the side-effects won't (hopefully)
    # confuse the abstract interpretation, which is supposed to terminate on any kind of code
    mm = get_single_method_match(tt, InferenceParams(newanalyzer).MAX_METHODS, get_world_counter(newanalyzer))
    result = abstract_call_method(newanalyzer, mm.method, mm.spec_types, mm.sparams, false, sv)

    rt = @static @isdefined(MethodCallResult) ? result.rt : first(result)

    # corresponding to the same logic in `analyze_frame!`
    if rt === Bottom
        if !isempty(get_uncaught_exceptions(newanalyzer))
            append!(get_reports(newanalyzer), get_uncaught_exceptions(newanalyzer))
        end
    end

    append!(get_reports(analyzer), get_reports(newanalyzer))
end

# works within inter-procedural context
function CC.abstract_call_method(analyzer::AbstractAnalyzer, method::Method, @nospecialize(sig), sparams::SimpleVector, hardlimit::Bool, sv::InferenceState)
    ret = @invoke abstract_call_method(analyzer::AbstractInterpreter, method::Method, sig, sparams::SimpleVector, hardlimit::Bool, sv::InferenceState)

    update_reports!(analyzer, sv)

    return ret
end

function update_reports!(analyzer::AbstractAnalyzer, sv::InferenceState)
    rs = get_to_be_updated(analyzer)
    if !isempty(rs)
        vf = get_virtual_frame(analyzer, sv)
        for r in rs
            pushfirst!(r.vst, vf)
        end
        empty!(rs)
    end
end

@static if isdefined(CC, :abstract_invoke)

import .CC:
    abstract_invoke,
    InvokeCallInfo,
    instanceof_tfunc

function CC.abstract_invoke(analyzer::AbstractAnalyzer, argtypes::Vector{Any}, sv::InferenceState)
    ret = @invoke abstract_invoke(analyzer::AbstractInterpreter, argtypes::Vector{Any}, sv::InferenceState)

    # if the `abstract_invoke(::AbstractInterpreter)` was successful, we need to update reports
    # since it's an inter-procedural inference that internally uses `typeinf_edge`
    info = ret.info
    if isa(info, InvokeCallInfo)
        update_reports!(analyzer, sv)
    end

    ReportPass(analyzer)(InvalidInvokeErrorReport, analyzer, sv, ret, argtypes)

    return ret
end

@reportdef struct InvalidInvokeErrorReport <: InferenceErrorReport
    argtypes::Vector{Any}
end
function get_msg(::Type{InvalidInvokeErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, argtypes::Vector{Any})
    fallback_msg = "invalid invoke" # mostly because of runtime unreachable

    ft = widenconst(argtype_by_index(argtypes, 2))
    ft === Bottom && return fallback_msg
    t = argtype_by_index(argtypes, 3)
    (types, isexact, isconcrete, istype) = instanceof_tfunc(t)
    if types === Bottom
        if isa(t, Const)
            type = typeof(t.val)
            return "argument type should be `Type`-object (given `$type`)"
        end
        return fallback_msg
    end

    argtype = argtypes_to_type(argtype_tail(argtypes, 4))
    nargtype = typeintersect(types, argtype)
    @assert nargtype === Bottom
    return "actual argument type (`$argtype`) doesn't intersect with specified argument type (`$types`)"
end

function (::SoundBasicPass)(::Type{InvalidInvokeErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, ret::CallMeta, argtypes::Vector{Any})
    if ret.rt === Bottom
        # here we report error that happens at the call of `invoke` itself.
        # if the error type (`Bottom`) is propagated from the `invoke`d call, the error has
        # already been reported within `typeinf_edge`, so ignore that case
        if !isa(ret.info, InvokeCallInfo)
            add_new_report!(InvalidInvokeErrorReport(analyzer, sv, argtypes), analyzer)
        end
    end
end

end # @static if isdefined(CC, :abstract_invoke)

function CC.abstract_eval_special_value(analyzer::AbstractAnalyzer, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    toplevel = istoplevel(analyzer, sv)
    if toplevel
        if isa(e, Slot) && is_global_slot(analyzer, e)
            if get_slottype((sv, get_currpc(sv)), e) === Bottom
                # if this abstract global variable is not initialized, form the global
                # reference and abstract intepret it; we may have abstract interpreted this
                # variable and it may have a type
                # if it's really not defined, the error will be generated later anyway
                e = GlobalRef(get_toplevelmod(analyzer), get_slotname(sv, e))
            end
        elseif isa(e, Symbol)
            # (already concretized) toplevel global symbols
            e = GlobalRef(get_toplevelmod(analyzer), e)
        end
    end

    ret = @invoke abstract_eval_special_value(analyzer::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    if isa(e, GlobalRef)
        mod, name = e.mod, e.name
        if isdefined(mod, name)
            if name === :Main
                # special case and propagate `Main` module as constant
                # XXX this was somewhat critical for accuracy and performance, but I'm not sure this still holds
                ret = Const(Main)
            elseif toplevel
                # here we will eagerly propagate the type of this global variable
                # of course the traced type might be difference from its type in actual execution
                # e.g. we don't track a global variable assignment wrapped in a function,
                # but it's highly possible this is a toplevel callsite and we need to take a
                # risk here, otherwise we can't enter the analysis !
                val = getfield(mod, name)
                ret = isa(val, AbstractGlobal) ? val.t : Const(val)
            else
                # TODO: add report pass to detect non-constant global variable here (for performance linting)
            end
        else
            # report pass for undefined global reference
            ReportPass(analyzer)(GlobalUndefVarErrorReport, analyzer, sv, mod, name)

            # `ret` at this point should be annotated as `Any` by `NativeInterpreter`, and we
            # just pass it as is to collect as much error points as possible within this frame
            # IDEA we can change it to `Bottom` to suppress any further abstract interpretation
            # with this variable, but at the same time we need to make sure we will invalidate
            # the cache for this frame on the future definition of this (currently) undefined binding
            # return Bottom
        end
    end

    return ret
end

@reportdef struct GlobalUndefVarErrorReport <: InferenceErrorReport
    mod::Module
    name::Symbol
end
get_msg(::Type{GlobalUndefVarErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol) =
    "variable $(mod).$(name) is not defined"

function (::SoundPass)(::Type{GlobalUndefVarErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol)
    add_new_report!(GlobalUndefVarErrorReport(analyzer, sv, mod, name), analyzer)
end

function (::BasicPass)(::Type{GlobalUndefVarErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol)
    is_corecompiler_undefglobal(mod, name) && return
    add_new_report!(GlobalUndefVarErrorReport(analyzer, sv, mod, name), analyzer)
end

"""
    is_corecompiler_undefglobal

Returns `true` if this global reference is undefined inside `Core.Compiler`, but the
  corresponding name exists in the `Base` module.
`Core.Compiler` reuses the minimum amount of `Base` code and there're some of missing
definitions, and `BasicPass` will exclude reports on those undefined names since they
usually don't matter and `Core.Compiler`'s basic functionality is battle-tested and validated
exhausively by its test suite and real-world usages
"""
is_corecompiler_undefglobal(mod::Module, name::Symbol) =
    return mod === CC ? isdefined(Base, name) :
           mod === CC.Sort ? isdefined(Base.Sort, name) :
           false

function CC.abstract_eval_value(analyzer::AbstractAnalyzer, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke abstract_eval_value(analyzer::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)

    # HACK if we encounter `_INACTIVE_EXCEPTION`, it means `ConcreteInterpreter` tried to
    # concretize an exception which was not actually thrown – yet the actual error hasn't
    # happened thanks to JuliaInterpreter's implementation detail, i.e. JuliaInterpreter
    # could retrieve `FrameData.last_exception`, which is initialized with
    # `_INACTIVE_EXCEPTION.instance` – but it's obviously not a sound approximation of an
    # actual execution and so here we will fix it to `Any`, since we don't analyze types of
    # exceptions in general
    if ret ⊑ _INACTIVE_EXCEPTION
        ret = Any
    end

    # report non-boolean condition error
    stmt = get_stmt((sv, get_currpc(sv)))
    if isa(stmt, GotoIfNot)
        t = widenconst(ret)
        if t !== Bottom
            ReportPass(analyzer)(NonBooleanCondErrorReport, analyzer, sv, t)
            # if this condition leads to an "non-boolean (t) used in boolean context" error,
            # we can turn it into Bottom and bail out early
            # TODO upstream this
            if typeintersect(Bool, t) !== Bool
                ret = Bottom
            end
        end
    end

    return ret
end

@reportdef struct NonBooleanCondErrorReport <: InferenceErrorReport
    @nospecialize(t::Union{Type,Vector{Type}})
end
get_msg(::Type{NonBooleanCondErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, @nospecialize(t::Type)) =
    "non-boolean ($t) used in boolean context"
get_msg(::Type{NonBooleanCondErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, ts::Vector{Type}) =
    "for $(length(ts)) of union split cases, non-boolean ($(join(ts, ", "))) used in boolean context"

function (::SoundPass)(::Type{NonBooleanCondErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, @nospecialize(t))
    if isa(t, Union)
        ts = Type[]
        for t in Base.uniontypes(t)
            if !(t ⊑ Bool)
                push!(ts, t)
            end
        end
        if !isempty(ts)
            add_new_report!(NonBooleanCondErrorReport(analyzer, sv, ts), analyzer)
        end
    else
        if !(t ⊑ Bool)
            add_new_report!(NonBooleanCondErrorReport(analyzer, sv, t), analyzer)
        end
    end
end

function (::BasicPass)(::Type{NonBooleanCondErrorReport}, analyzer::AbstractAnalyzer, sv::InferenceState, @nospecialize(t))
    if isa(t, Union)
        ts = Type[]
        for t in Base.uniontypes(t)
            if typeintersect(Bool, t) !== Bool
                # TODO move this to abstractinterpretation.jl
                if JETAnalysisParams(analyzer).strict_condition_check ||
                   !(t <: Function || # !(::Function)
                     t === Missing || # ==(::Missing, ::Any), ==(::Any, ::Missing), ...
                     false)
                    push!(ts, t)
                end
            end
        end
        if !isempty(ts)
            add_new_report!(NonBooleanCondErrorReport(analyzer, sv, ts), analyzer)
        end
    else
        if typeintersect(Bool, t) !== Bool
            add_new_report!(NonBooleanCondErrorReport(analyzer, sv, t), analyzer)
        end
    end
end

function CC.abstract_eval_statement(analyzer::AbstractAnalyzer, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    if istoplevel(analyzer, sv)
        if get_concretized(analyzer)[get_currpc(sv)]
            return Any # bail out if it has been interpreted by `ConcreteInterpreter`
        end
    end

    return @invoke abstract_eval_statement(analyzer::AbstractInterpreter, e, vtypes::VarTable, sv::InferenceState)
end

function CC.finish(me::InferenceState, analyzer::AbstractAnalyzer)
    @invoke finish(me::InferenceState, analyzer::AbstractInterpreter)

    if istoplevel(analyzer, me)
        # find assignments of abstract global variables, and assign types to them,
        # so that later analysis can refer to them

        stmts = me.src.code
        cfg = compute_basic_blocks(stmts)
        assigns = Dict{Int,Bool}() # slot id => is this deterministic
        for (pc, stmt) in enumerate(stmts)
            if @isexpr(stmt, :(=))
                lhs = first(stmt.args)
                if isa(lhs, Slot)
                    slot = slot_id(lhs)
                    if is_global_slot(analyzer, slot)
                        isnd = is_assignment_nondeterministic(cfg, pc)

                        # COMBAK this approach is really not true when there're multiple
                        # assignments in different basic blocks
                        if haskey(assigns, slot)
                            assigns[slot] |= isnd
                        else
                            assigns[slot] = isnd
                        end
                    end
                end
            end
        end

        if !isempty(assigns)
            slottypes = collect_slottypes(me)
            for (slot, isnd) in assigns
                slotname = get_global_slots(analyzer)[slot]
                typ = slottypes[slot]
                set_abstract_global!(analyzer, get_toplevelmod(analyzer), slotname, typ, isnd, me)
            end
        end
    end
end

# simple cfg analysis to check if the assignment at `pc` will happen non-deterministically
function is_assignment_nondeterministic(cfg::CFG, pc::Int)
    isnd = false

    blocks = cfg.blocks
    for (idx, block) in enumerate(blocks)
        if pc in rng(block)
            for block′ in blocks
                succs = block′.succs
                if idx in succs
                    isnd |= length(succs) > 1
                end
            end
        end
    end

    return isnd
end

# at this point all the types of SSA values are iterated to maximum fixed point,
# and we can compute types of slot as least upper bound of types of all the possible
# assignment of the slot (the type of assignment statement is available as SSA value type)
# the implementation is mostly same as `record_slot_assign!(sv::InferenceState)`, but
# we don't `widenconst` each SSA value type
function collect_slottypes(sv::InferenceState)
    states = sv.stmt_types
    ssavaluetypes = sv.src.ssavaluetypes::Vector{Any}
    stmts = sv.src.code::Vector{Any}
    slottypes = Any[Bottom for _ in sv.src.slottypes::Vector{Any}]
    for i = 1:length(stmts)
        stmt = stmts[i]
        state = states[i]
        # find all reachable assignments to locals
        if isa(state, VarTable) && @isexpr(stmt, :(=))
            lhs = first(stmt.args)
            if isa(lhs, Slot)
                vt = ssavaluetypes[i] # don't widen const
                if vt !== Bottom
                    id = slot_id(lhs)
                    otherTy = slottypes[id]
                    slottypes[id] = tmerge(otherTy, vt)
                end
            end
        end
    end
    return slottypes
end

function set_abstract_global!(analyzer::AbstractAnalyzer, mod::Module, name::Symbol, @nospecialize(t), isnd::Bool, sv::InferenceState)
    prev_agv = nothing
    prev_t = nothing
    iscd = is_constant_declared(name, sv)

    # check if this global variable is already assigned previously
    if isdefined(mod, name)
        val = getfield(mod, name)
        if isa(val, AbstractGlobal)
            prev_t = val.t
            if val.iscd && (prev_t′ = widenconst(prev_t)) !== (t′ = widenconst(t))
                warn_invalid_const_global!(name)
                ReportPass(analyzer)(InvalidConstantRedefinition, analyzer, sv, mod, name, prev_t′, t′)
                return
            end
            prev_agv = val
        else
            prev_t = Core.Typeof(val)
            if isconst(mod, name)
                invalid = prev_t !== (t′ = widenconst(t))
                if invalid || !isa(t, Const)
                    warn_invalid_const_global!(name)
                    if invalid
                        ReportPass(analyzer)(InvalidConstantRedefinition, analyzer, sv, mod, name, prev_t, t′)
                    end
                    return
                end
                # otherwise, we can just redefine this constant, and Julia will warn it
                ex = iscd ? :(const $name = $(QuoteNode(t.val))) : :($name = $(QuoteNode(t.val)))
                return Core.eval(mod, ex)
            end
        end
    end

    isnew = isnothing(prev_t)

    # if this constant declaration is invalid, just report it and bail out
    if iscd && !isnew
        warn_invalid_const_global!(name)
        ReportPass(analyzer)(InvalidConstantDeclaration, analyzer, sv, mod, name)
        return
    end

    # if this assignment happens non-deterministically, we need to take the previous type into account
    if isnd
        if !isnew # if this assignment is an initialization, we just need to use `t`
            t = tmerge(prev_t, t)
        end
    else
        # if this assignment happens deterministically, and the assigned value is known to be
        # constant statically, let's concretize it for good reasons;
        # we will be able to use it in concrete interpretation and so this allows to define
        # structs with type aliases, etc.
        if isa(t, Const)
            if iscd
                @assert isnew # means, this is a valid constant declaration
                return Core.eval(mod, :(const $name = $(QuoteNode(t.val))))
            else
                # we've checked `mod.name` wasn't declared as constant previously
                return Core.eval(mod, :($name = $(QuoteNode(t.val))))
            end
        end
    end

    # okay, we will define new abstract global variable from here on
    if isa(prev_agv, AbstractGlobal)
        return Core.eval(mod, :(let name = $name::$AbstractGlobal
            name.t = $t
            name
        end))
    else
        return Core.eval(mod, :($name = $(AbstractGlobal(t, iscd))))
    end
end

warn_invalid_const_global!(name::Symbol) = @warn """
JET.jl can't update the definition of this constant declared global variable: `$name`
This may fail, cause incorrect analysis, or produce unexpected errors.
"""

# IDEA we may want to hoist `InvalidConstXXX` errors into top-level errors

@reportdef struct InvalidConstantRedefinition <: InferenceErrorReport
    mod::Module
    name::Symbol
    @nospecialize(t′::Any)
    @nospecialize(t::Any)
end
get_msg(::Type{InvalidConstantRedefinition}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol, @nospecialize(t′::Any), @nospecialize(t::Any)) =
    "invalid redefinition of constant $(mod).$(name) (from $(t′) to $(t))"

function (::SoundBasicPass)(::Type{InvalidConstantRedefinition}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol, @nospecialize(prev_t), @nospecialize(t))
    add_new_report!(InvalidConstantRedefinition(analyzer, sv, mod, name, prev_t, t), analyzer)
end

@reportdef struct InvalidConstantDeclaration <: InferenceErrorReport
    mod::Module
    name::Symbol
end
get_msg(::Type{InvalidConstantDeclaration}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol) =
    "cannot declare a constant $(mod).$(name); it already has a value"

function (::SoundBasicPass)(::Type{InvalidConstantDeclaration}, analyzer::AbstractAnalyzer, sv::InferenceState, mod::Module, name::Symbol)
    add_new_report!(InvalidConstantDeclaration(analyzer, sv, mod, name), analyzer)
end

function is_constant_declared(name::Symbol, sv::InferenceState)
    return any(sv.src.code) do @nospecialize(x)
        if @isexpr(x, :const)
            arg = first(x.args)
            # `transform_global_symbols!` replaces all the global symbols in this toplevel frame with `Slot`s
            if isa(arg, Slot)
                return get_slotname(sv, arg) === name
            end
        end
        return false
    end
end


text: Dict{String, Any}()

ext: jl


text:
struct TransformationStructure{T <: Function,U <: Function}
    name::Union{String,Nothing}
    p_transform::T
    inv_p_transform::U
end

"""
returns TransformationStructure that flips the signs of negative parameters, and then does the transform p -> log(p).
"""
function logabs_transform(p0)
    name = "logabs_"
    is_positive = convert.(Int64, p0 .>= 0)
    is_positive[is_positive .== 0] .= -1
    pos(p) = p .* is_positive
    p_transform(p) = log.(p .* is_positive)
    inv_p_transform(logp) = exp.(logp) .* is_positive

    return TransformationStructure(name, p_transform, inv_p_transform)
end

"""
returns TransformationStructure that fixes parameters[indices]
"""
function fix_params(p0, indices)
    indices |> unique! |> sort!
    not_indices = setdiff(collect(1:length(p0)), indices)
    p0 = deepcopy(p0)
    function p_transform(p)
        return deleteat!(deepcopy(p), indices)
    end

    function inv_p_transform(p)
        out = [p[1] for el in p0]
        out[not_indices] .= p
        out[indices] .= p0[indices]
        return out
    end
    return TransformationStructure(nothing, p_transform, inv_p_transform)
end

"""
returns TransformationStructure. params[indices] -> biases.*params[indices]
"""
function bias_transform(p0, indices, biases)
    indices |> unique! |> sort!
    not_indices = setdiff(collect(1:length(p0)), indices)
    all_biases = ones(size(p0))
    all_biases[indices] = biases
    all_inv_biases = 1. ./ all_biases

    function p_transform(p)
        return  p .* all_biases
    end

    function inv_p_transform(p)
        return p .* all_inv_biases
    end

    return TransformationStructure(nothing, p_transform, inv_p_transform)
end

"""
returns TransformationStructure that fixes **all but** parameters[indices]
"""
function only_free_params(p0, indices)
    name = "only free indices are $indices"
    indices |> unique! |> sort!
    not_indices = setdiff(collect(1:length(p0)), indices)
    return fix_params(p0, not_indices)
end

"""
    transform_cost(cost, p0, tr::TransformationStructure; unames=nothing, pnames=nothing)
    return DiffCost(new_cost, new_cost2), newp0

Given a cost function C(p), makes a new differentiable cost function D(q), where q = tr(p) and D(q) = C(p)
"""
function transform_cost(cost, p0, tr::TransformationStructure; unames=nothing, pnames=nothing)
    newp0 = tr.p_transform(p0)
    jac = ForwardDiff.jacobian

    function new_cost(p)
        return p |> tr.inv_p_transform |> cost
    end

    function new_cost2(p, g)
        orig_p = tr.inv_p_transform(p)
        orig_grad  = deepcopy(orig_p)
        val = cost(orig_p, orig_grad)
        g[:] = jac(tr.inv_p_transform, p) * orig_grad
        return val
    end
    return DiffCost(new_cost, new_cost2), newp0
end

"""
    new_od = transform_ODESystem(od::ODESystem, tr::TransformationStructure)

    - reparameterises the parameters of an ODE system via the transformation tr.
    - within the ODE eqs, any instance of a parameter p is replaced with tr.inv_p_transform(p)
    - if there are default parameter values p0, they are changed to tr.p_transform(p0)
"""
function transform_ODESystem(od::ModelingToolkit.AbstractSystem, tr::TransformationStructure)
    t = ModelingToolkit.get_iv(od)
    unames = ModelingToolkit.get_states(od)
    eqs = ModelingToolkit.get_eqs(od)
    ps = ModelingToolkit.get_ps(od)
    new_ps = transform_names(ps, tr) # modified names under transformation
    new_p0 = tr.p_transform
    of = ODEFunction(od, eval_expression=false) # to solve world age issues
    rhs = similar(unames, Any)
    of(rhs, unames, tr.inv_p_transform(new_ps), t) # in place modification of rhs
    lhs = [el.lhs for el in  eqs]

    _defaults = ModelingToolkit.get_defaults(od)

    if length(_defaults) >  0
        p_collect = [el => _defaults[el] for el in ps] # in correct order
        if length(p_collect) > 0
            new_p_vals = tr.p_transform(last.(p_collect))
            new_p_dict = Dict(new_ps .=> new_p_vals)
            _defaults = merge(_defaults, new_p_dict)
        end
    end

    de = ODESystem(lhs .~ rhs, t, unames, new_ps,  defaults=_defaults, name=nameof(od))
    return de # (vars .=> last.(ic)), (new_ps .=> newp0)
end

"""
Reparameterises prob::ODEProblem via the transformation tr. so newprob.p = tr(p) is an equivalent ODEProblem to prob.p = p
"""

function transform_problem(prob::ODEProblem, tr::TransformationStructure; unames=nothing, pnames=nothing)
    println(pnames)
    sys = modelingtoolkitize(prob)
    eqs = ModelingToolkit.get_eqs(sys)
    pname_tr = ModelingToolkit.get_ps(sys) .=> pnames
    uname_tr = ModelingToolkit.get_states(sys) .=> unames
    neweqs = eqs
    if !(pnames === nothing)
        neweqs = [el.lhs ~ substitute(el.rhs, pname_tr) for el in neweqs]
    else
        pnames = ModelingToolkit.get_ps(sys)
    end

    if !(unames === nothing)
        neweqs = [substitute(el.lhs, uname_tr) ~ substitute(el.rhs, uname_tr) for el in neweqs]
    else
        unames = ModelingToolkit.get_states(sys)
    end
    named_sys = ODESystem(neweqs, ModelingToolkit.get_iv(sys), unames, pnames,  defaults=merge(Dict(unames .=> prob.u0), Dict(pnames .=> prob.p)), name=nameof(sys))
    newp0 = tr.p_transform(prob.p)
    t_sys = transform_ODESystem(named_sys, tr)
    return t_sys, (ModelingToolkit.get_states(t_sys) .=> prob.u0),
     (ModelingToolkit.get_ps(t_sys) .=> newp0)
end

"""
Transforms a set of parameter names via the name provided in tr
The output names are ModelingToolkit.Variable types
# Example

nv = [m, c, k]
tr.name = log
returns the variables: [log(m)], log(c), log(k)]
"""
function transform_names(nv, tr::TransformationStructure)
    if tr.name === nothing
        names = Symbol.(repr.(tr.p_transform(nv)))
    else
        names = Symbol.(tr.name .* repr.(nv))
    end
    new_vars = [Num(Variable(el)) for el in names]
    return new_vars
end

"""
searches for the parameter indices corresponding to names.
ps is an array of names
"""
function get_name_ids(ps, names::Array{String,1})
    # can't do a single findall as this doesn't preserve ordering
    all_names = repr.(ps)
    ids = [first(findall(x -> x == names[i], all_names)) for (i, el) in enumerate(names)]
    return ids
end

"""
searches for the parameter indices corresponding to names.
ps is an array of pairs: names .=> vals
"""
function get_name_ids(ps::Array{Pair{T,U},1}, names::Array{String,1}) where T where U
    return get_name_ids(first.(ps), names)
end


text:
let
    order = 2

    # Haverly1 standard pooling instance
    x14, x24, x35, x36, x45, x46, w4 = variables(["x_{14}", "x_{24}", "x_{35}", "x_{36}", "x_{45}", "x_{46}", "w_4"])

    c14, c24, c35, c36, c45, c46 = [6, 16, 1, -9, -5, -15]
    q1, q2, q3, q5, q6 = [ 3, 1, 2, 5//2, 3//2 ]
    ub1, ub2, ub3, ub4, ub5, ub6 = [ 300, 300, 300, 300, 100, 200 ] // 300

    flow_eq    = [ x14+x24-(x45+x46) ]
    blend_eq   = [ w4*(x45+x46)-(q1*x14+q2*x24) ]
    qual_bnd   = [ q5*(x45+x36)-(w4*x45+q3*x35), q6*(x46+x36)-(w4*x46+q3*x36) ]
    cap_bnd    = [ ub1-x14, ub2-x24, ub3-(x35+x36), ub4-(x45+x46), ub5-(x35+x45), ub6-(x36+x46)]
    flow_bnd   = [ x14,x24,x35,x36,x45,x46]

    obj = dot([c14, c24, c35, c36, c45, c46], [x14, x24, x35, x36, x45, x46])
    prob = momentprob(order, obj, [qual_bnd, cap_bnd, flow_bnd], [flow_eq,blend_eq])

    y, objval, solsta = solve_mosek(prob)
    @test norm(y[1:8]-[1,0,1/3,0,1/3,0,1/3,1]) < 1e-4
end


text:
using Documenter, GenomicFeatures

makedocs(
    format = :html,
    sitename = "GenomicFeatures.jl",
    pages = [
        "Home" => "index.md",
        "Intervals" => "intervals.md",
        "I/O" => [
            "BED" => "io/bed.md",
            "GFF3" => "io/gff3.md",
            "BigWig" => "io/bigwig.md",
            "BigBed" => "io/bigbed.md"
        ]
    ],
    authors = "Kenta Sato, D. C. Jones, Ben J. Ward, The BioJulia Organisation and other contributors."
)
deploydocs(
    repo = "github.com/BioJulia/GenomicFeatures.jl.git",
    julia = "0.6",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
theme(:wong)
gr(aspect_ratio=1,legend=:bottomright,axis=nothing,border=nothing,ticks=nothing,legendfontsize=14,xlims=[-1,7.3],
   markerstrokewidth=0,markersize=1)
r1D, w1D = gauss_quad(0,0,2)
Nfp = length(r1D) # number of points per face
e = ones(Nfp) # vector of all ones
z = zeros(Nfp) # vector of all zeros
rf = [r1D; -r1D; -e];
sf = [-e; r1D; -r1D];

r_tri,s_tri,_ = Basis2DTri.quad_nodes_tri(4)
plot([-1;-1],[-1;1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;1],[1;-1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;1],[-1;-1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;-1].+3/4,[-1;1].+2/4,seriescolor=:black,lw=3,ls=:dash,label = nothing)
plot!([-1;1].+3/4,[1;-1].+2/4,seriescolor=:black,lw=3,ls=:dash,label = nothing)
plot!([-1;1].+3/4,[-1;-1].+2/4,seriescolor=:black,lw=3,ls=:dash,label = nothing)
plot!([-1;-1].+3,[-1;1].+2,seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+3,[1;-1].+2,seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+3,[-1;-1].+2,seriescolor=:black,lw=3,label = nothing)
plot!([1;4],[-1,1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;2],[-1,1],seriescolor=:black,lw=3,ls=:dash,label = nothing)
plot!([-1;2],[1,3],seriescolor=:black,lw=3,label = nothing)
scatter!(r_tri.+3/4,s_tri.+2/4,markersize=5,markercolor=:deepskyblue,markerstrokecolor=:black,markerstrokewidth=2,label = nothing)
scatter!(rf.+3/4,sf.+2/4,markersize=5,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2,label = nothing)
scatter!(r_tri.+3,s_tri.+2,markersize=5,markercolor=:deepskyblue,markerstrokecolor=:black,markerstrokewidth=2,label = "Volume quadrature points")
scatter!(rf.+3,sf.+2,markersize=5,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=2,label = "Face quadrature points")
scatter!([1.5;2;2.5],[0;0.5*2/3;1*2/3],markercolor=:black,markersize=4,label = nothing)
savefig("reference_elem.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using BenchmarkTools
using JLD2

"""
```
print_all_benchmarks(trial, ref, trial_name; stats = [:min, :median, :mean, :max], factors = [:times, :gctimes, :memory, :allocs])
```

For printing the benchmark comparisons given by the BenchmarkTools package of a single trial.

### Arguments
- `trial::BenchmarkTools.Trial`: A trial of type BenchmarkTools.Trial, outputted from `@benchmark`.
- `ref::String`: The location of the previous trial, saved by the `write_ref_trial` method
- `trial_name::String`: The name of the trial being benchmarked, i.e. `forecast_one`, if the `forecast_one` function is the trial being benchmarked
"""
function print_all_benchmarks(trial::BenchmarkTools.Trial, ref::String,
                              trial_name::String;
                              stats::Vector{Symbol} = [:min, :median, :mean, :max],
                              factors::Vector{Symbol} = [:times, :gctimes, :memory, :allocs])
    print_all_benchmarks(trial_to_dict(trial), ref, trial_name, stats = stats, factors = factors)
end

# For a single trial comparison
function print_all_benchmarks(trial::Dict{Symbol, Union{Vector, Int}}, ref::String,
                              trial_name::String;
                              stats::Vector{Symbol} = [:min, :median, :mean, :max],
                              factors::Vector{Symbol} = [:times, :gctimes, :memory, :allocs])
    prev_trial = load(ref, trial_name)
    print_all_benchmarks(trial, prev_trial, trial_name, stats = stats, factors = factors)
end

"""
```
print_all_benchmarks(trial, prev_trial, trial_name; stats = [:min, :median, :mean, :max], factors = [:times, :gctimes, :memory, :allocs])
```

The underlying method for printing the benchmark comparisons given by the BenchmarkTools package of a single trial.

### Arguments
- `trial::Dict{Symbol, Union{Vector, Int}}`: A dictionary mapping the `factors` to their values for an individual trial.
- `prev_trial::Dict{Symbol, Union{Vector, Int}}`: The comparison factors that `trial` is being benchmarked against.
- `trial_name::String`: The name of the trial that is being benchmarked, i.e. `forecast_one`, if the `forecast_one` function is the trial being benchmarked.

### Keyword Arguments
- `stats::Vector{Symbol}`: The statistics that can be shown when printing the benchmarks, which include `[:min, :median, :mean, :max]`.
- `factors::Vector{Symbol}`: The different factors within the trial that can be shown when printing the benchmarks, which include `[:times, :gctimes, :memory, :allocs]`.
- `group::Bool`: Whether or not this benchmark was part of a group of trials being benchmarked, for the purposes of indentation of the stdout.
"""
function print_all_benchmarks(trial::Dict{Symbol, Union{Vector, Int}},
                              prev_trial::Dict{Symbol, Union{Vector, Int}},
                              trial_name::String;
                              stats::Vector{Symbol} = [:min, :median, :mean, :max],
                              factors::Vector{Symbol} = [:times, :gctimes, :memory, :allocs],
                              group::Bool = false)
    # Print formatting of the trial sections
    trial_delimiter = repeat("-", length(trial_name))
    trial_name = group ? "  "*trial_name : trial_name
    trial_delimiter = group ? "  "*trial_delimiter : trial_delimiter
    println(trial_delimiter*"\n")
    print_with_color(:bold, trial_name)
    print("\n")

    # Printing benchmarks for each factor within a single trial
    for f in factors
        times = (f == :times)
        if f in [:memory, :allocs]
            print_single_benchmark(trial, prev_trial, :identity, f, times, group)
            print("\n")
            continue
        end
        for s in stats
            print_single_benchmark(trial, prev_trial, s, f, times, group)
        end
        print("\n")
    end
end

"""
```
print_all_benchmarks(group, ref, group_name, trial_names; stats = [:min, :median, :mean, :max], factors = [:times, :gctimes, :memory, :allocs])
```

For printing the benchmark comparisons given by the BenchmarkTools package of a group of trials.

### Arguments
- `group::Dict{Symbol, Dict}`: A dictionary mapping the `trial_names` to their corresponding trials, which are of type, `Dict{Symbol, Union{Vector, Int}}`.
- `ref::String`: The location of the .jld file containing the previous group to be benchmarked against.
- `trial_names::Vector{Symbol}`: The names of the trials being benchmarked, i.e. `forecast_one`, if the `forecast_one` function is one of the trials being benchmarked.

### Keyword Arguments
- `stats::Vector{Symbol}`: The statistics that can be shown when printing the benchmarks, which include `[:min, :median, :mean, :max]`.
- `factors::Vector{Symbol}`: The different factors within the trial that can be shown when printing the benchmarks, which include `[:times, :gctimes, :memory, :allocs]`.
"""
function print_all_benchmarks(group::Dict{Symbol, Dict}, ref::String,
                              group_name::String,
                              trial_names::Vector{Symbol};
                              stats::Vector{Symbol} = [:min, :median, :mean, :max],
                              factors::Vector{Symbol} = [:times, :gctimes, :memory, :allocs])
    # Load in the previous group to benchmark against
    prev_group = load(ref, group_name)

    # Print the group_name section header
    print_with_color(:green, group_name, bold = true)
    print("\n")

    # For each trial within a group, print all of the benchmarks.
    for trial_name in trial_names
        print_all_benchmarks(group[trial_name], prev_group[trial_name], string(trial_name),
                             stats = stats, factors = factors, group = true)
    end
end

"""
```
print_single_benchmark(trial, prev_trial, stat, factor, times)
```

### Arguments
- `trial::Dict{Symbol, Union{Vector, Int}}`: A dictionary mapping the `factors` to their values for an individual trial.
- `prev_trial::Dict{Symbol, Union{Vector, Int}}`: The comparison factors that `trial` is being benchmarked against.
- `stats::Symbol`: The statistic that will be shown when printing the benchmarks, one of `[:min, :median, :mean, :max]`.
- `factors::Symbol`: The factor within the trial that will be shown when printing the benchmarks, one of `[:times, :gctimes, :memory, :allocs]`.
- `times::Bool`: Whether or not the factor chosen was the `:times` factor, for the purposes of printing more grammatically correct output.
- `group::Bool`: Whether or not this benchmark was part of a group of trials being benchmarked, for the purposes of indentation of the stdout.
"""
function print_single_benchmark(trial::Dict{Symbol, Union{Vector, Int}},
                                prev_trial::Dict{Symbol, Union{Vector, Int}},
                                stat::Symbol, factor::Symbol,
                                times::Bool = false,
                                group::Bool = false)
    title = stat == :identity ? factor : Symbol(stat, " ", factor)
    calc_stat(x) = stat_maps()[stat](x)
    print_single_benchmark(title, calc_stat(trial[factor]),
                           calc_stat(prev_trial[factor]), times, group)
end

function print_single_benchmark(title::Symbol, new::T, old::T, times::Bool = false,
                                group::Bool = false) where T<:Real
    bad_state = times ? "slower" : "larger"
    good_state = times ? "faster" : "smaller"

    diff = round((new/old - 1)*100, digits = 2)
    if diff > 0
        state = bad_state
    else
        state = good_state
    end
    color = calculate_color(diff, title)

    # Make sure the printing is indented properly to more aesthetically display group trial results
    if group
        print_with_color(color, "  $title is $(abs(diff))% $state\n")
    else
        print_with_color(color, "$title is $(abs(diff))% $state\n")
    end
end


text:
include(joinpath(dirname(@__FILE__), "crank_nicolson_CS", "calc_line.jl"))

function test_CN()
  @testset "---- Crank-Nicolson Convergence Tests, Complex Step Jacobian -----" begin
    start_dir = pwd()

    cd(dirname(@__FILE__))
    resize!(ARGS, 1)

    # =================== CN, CS tests ===================
    cd("./crank_nicolson_CS/")
    println("======",pwd())

    cd("./m1")
    println("======", pwd())
    fname = "input_vals1.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("../m2")
    println("======", pwd())
    fname = "input_vals1.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

#    cd("../m3")
#    println("======", pwd())
#    fname = "input_vals1.jl"
#    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("..")
    println("======", pwd())
    fname = "calc_line.jl"  #???
#    mesh, sbp, eqn, opts = solvePDE(fname)

    slope = calc_line()
    # println("slope = ", slope)

    data = readdlm("err_data.dat")
    err_vals = data[:, 2]
    #println("err_vals = ", err_vals)

    slope_val = 2.00
    slope_margin = 0.1

    @test  slope  > slope_val - slope_margin
    @test  slope  < slope_val + slope_margin

    err_val = 0.09095728504176116
    slope_fac = 1.25
    # println("err_vals[1] = ", err_vals[1])
    @test  err_vals[1]  > err_val/slope_fac
    @test  err_vals[1]  < err_val*slope_fac

    cd("../")
    # =================== CN, FD tests ===================
    cd("./crank_nicolson_FD/")

    cd("./m1")
    fname = "input_vals1.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("../m2")
    fname = "input_vals1.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

#    cd("../m3")
#    fname = "input_vals1.jl"
#    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("..")
    fname = "calc_line.jl"
#    mesh, sbp, eqn, opts = solvePDE(fname)  #???

    slope = calc_line()
    # println("slope = ", slope)

    data = readdlm("err_data.dat")
    err_vals = data[:, 2]
    #println("err_vals = ", err_vals)

    slope_val = 2.00
    slope_margin = 0.1

    @test  slope  > slope_val - slope_margin
    @test  slope  < slope_val + slope_margin

    err_val = 0.09095728504176116
    slope_fac = 1.25
    # println("err_vals[1] = ", err_vals[1])
    @test  err_vals[1]  > err_val/slope_fac
    @test  err_vals[1]  < err_val*slope_fac

    cd("../")
    # =================== CN, PETSc CS tests ===================
    cd("./crank_nicolson_PETSc_serial/")

    cd("./m1")
    fname = "input_vals1.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("../m2")
    fname = "input_vals1.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

#    cd("../m3")
#    fname = "input_vals1.jl"
#    mesh, sbp, eqn, opts = solvePDE(fname)

    cd("..")
    fname = "calc_line.jl"
#    mesh, sbp, eqn, opts = solvePDE(fname) #???

    slope = calc_line()
    # println("slope = ", slope)

    data = readdlm("err_data.dat")
    err_vals = data[:, 2]
    #println("err_vals = ", err_vals)

    slope_val = 2.00
    slope_margin = 0.1

    @test  slope  > slope_val - slope_margin
    @test  slope  < slope_val + slope_margin

    err_val = 0.09095728504176116
    slope_fac = 1.25
    # println("err_vals[1] = ", err_vals[1])
    @test  err_vals[1]  > err_val/slope_fac
    @test  err_vals[1]  < err_val*slope_fac

    cd(start_dir)
  end  # end facts block

end

add_func1!(AdvectionTests, test_CN, [TAG_SHORTTEST, TAG_CN])


text:
using Mass2, Base.Test
let
  # make a single noiseless pulse with known properties
  num_points = 1024
  rise_points = 50
  fall_points = 100
  quiescent_value = 1000
  arrival_point = div(num_points,4)
  amplitude = 10000
  pulse_float = two_exponential_pulses(num_points, rise_points, fall_points, quiescent_value, (arrival_point,), (amplitude,))
  pulse = to_type_and_white_noise(UInt16, 0, pulse_float)
  pulse_int = convert(Vector{Int}, pulse) # for functions where uint causes problems
  frame_time = 5e-6 # time per point, seconds

  # summarize it
  out = compute_summary((pulse,),arrival_point,frame_time)

  (pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv,
  peak_index, peak_value, min_value) = map(first,out)
  npre = arrival_point # number of pretrigger points
  npost = length(pulse)-npre # number of posttrigger points

  @test pretrig_mean ≈ mean(pulse[1:npre]) # average of points in pretrigger
  @test pretrig_rms ≈ std(pulse[1:npre]) # std deviatino of points in pretrigger
  @test pulse_average ≈ mean(pulse[npre+1:end]-quiescent_value)# pulse average deviation from pretrig_mean
  @test pulse_rms ≈ sqrt(mean((pulse_int[npre+1:end]-quiescent_value).^2)) # pulse rms relative to pretrig_rms
  @test rise_time ≈ frame_time*(find(pulse.>10000)[1]-find(pulse.>2000)[1]) # defined as 10%-90% rise time
  # @test postpeak_deriv<minimum(diff(pulse_int))/frame_time # for a noiseless pulse, the postpeak deriv should be negative
  @test peak_index ≈ indmax(pulse) # index of the maximum value
  @test peak_value+pretrig_mean ≈ maximum(pulse) #peak value is defined relative to pretrig_mean
  @test min_value ≈ minimum(pulse) #min_value is the minimum value, not defined relative to anything
end

let
  # make a piled up noiseless pulse with known properties, peak is on 2nd pulse
  num_points = 1024
  rise_points = 50
  fall_points = 100
  quiescent_value = 1000
  arrival_point = div(num_points,4)
  arrival_point2 = div(num_points,2)
  amplitude = 10000
  amplitude2 = 10000
  pulse_float = two_exponential_pulses(num_points, rise_points, fall_points, quiescent_value,
                (arrival_point,arrival_point2), (amplitude, amplitude2))
  pulse = to_type_and_white_noise(UInt16, 0, pulse_float)
  pulse_int = convert(Vector{Int}, pulse) # for functions where uint causes problems
  frame_time = 5e-6 # time per point, seconds

  # summarize it
  out = compute_summary((pulse,),arrival_point,frame_time)

  (pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv,
  peak_index, peak_value, min_value) = map(first,out)
  npre = arrival_point # number of pretrigger points
  npost = length(pulse)-npre # number of posttrigger points

  @test pretrig_mean ≈ mean(pulse[1:npre]) # average of points in pretrigger
  @test pretrig_rms ≈ std(pulse[1:npre]) # std deviatino of points in pretrigger
  @test pulse_average ≈ mean(pulse[npre+1:end]-quiescent_value)# pulse average deviation from pretrig_mean
  @test pulse_rms ≈ sqrt(mean((pulse_int[npre+1:end]-quiescent_value).^2)) # pulse rms relative to pretrig_rms
  @test rise_time > frame_time*(find(pulse.>10000)[1]-find(pulse.>2000)[1]) # defined as 10%-90% rise time, will be confused by pileup
  # @test postpeak_deriv<minimum(diff(pulse_int))/frame_time # for a noiseless pulse, the postpeak deriv should be negative
  @test peak_index ≈ indmax(pulse) # index of the maximum value
  @test peak_value+pretrig_mean ≈ maximum(pulse) #peak value is defined relative to pretrig_mean
  @test min_value ≈ minimum(pulse) #min_value is the minimum value, not defined relative to anything
end

let
  # make a piled up noiseless pulse with known properties, peak is on 1st pulse
  num_points = 1024
  rise_points = 50
  fall_points = 100
  quiescent_value = 1000
  arrival_point = div(num_points,4)
  arrival_point2 = div(num_points,2)
  amplitude = 10000
  amplitude2 = 5000
  pulse_float = two_exponential_pulses(num_points, rise_points, fall_points, quiescent_value,
                (arrival_point,arrival_point2), (amplitude, amplitude2))
  pulse = to_type_and_white_noise(UInt16, 0, pulse_float)
  pulse_int = convert(Vector{Int}, pulse) # for functions where uint causes problems
  frame_time = 5e-6 # time per point, seconds

  # summarize it
  out = compute_summary((pulse,),arrival_point,frame_time)

  (pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv,
  peak_index, peak_value, min_value) = map(first,out)
  npre = arrival_point # number of pretrigger points
  npost = length(pulse)-npre # number of posttrigger points

  @test pretrig_mean ≈ mean(pulse[1:npre]) # average of points in pretrigger
  @test pretrig_rms ≈ std(pulse[1:npre]) # std deviatino of points in pretrigger
  @test pulse_average ≈ mean(pulse[npre+1:end]-quiescent_value)# pulse average deviation from pretrig_mean
  @test pulse_rms ≈ sqrt(mean((pulse_int[npre+1:end]-quiescent_value).^2)) # pulse rms relative to pretrig_rms
  @test rise_time ≈ frame_time*(find(pulse.>10000)[1]-find(pulse.>2000)[1]) # defined as 10%-90% rise time
  @test postpeak_deriv ≈ maximum(diff(pulse_int)[peak_index:end]) # for a noiseless pulse, the postpeak deriv should be negative
  @test peak_index ≈ indmax(pulse) # index of the maximum value
  @test peak_value+pretrig_mean ≈ maximum(pulse) #peak value is defined relative to pretrig_mean
  @test min_value ≈ minimum(pulse) #min_value is the minimum value, not defined relative to anything
end

let
  # make a single noisy pulse with known properties
  num_points = 10240
  rise_points = 50
  fall_points = 100
  quiescent_value = 1000
  arrival_point = div(num_points,4)
  amplitude = 10000
  noise_std = 1
  pulse_float = two_exponential_pulses(num_points, rise_points, fall_points, quiescent_value, (arrival_point,), (amplitude,))
  pulse = to_type_and_white_noise(UInt16, noise_std, pulse_float)
  pulse_int = convert(Vector{Int}, pulse) # for functions where uint causes problems
  frame_time = 5e-6 # time per point, seconds

  # summarize it
  out = compute_summary((pulse,),arrival_point,frame_time)

  (pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv,
  peak_index, peak_value, min_value) = map(first,out)
  npre = arrival_point # number of pretrigger points
  npost = length(pulse)-npre # number of posttrigger points

  @test isapprox(pretrig_mean, mean(pulse[1:npre]) ;rtol=1e-3) # average of points in pretrigger
  @test isapprox(pretrig_rms, std(pulse[1:npre]) ;rtol=1e-3) # std deviatino of points in pretrigger
  @test isapprox(pulse_average, mean(pulse_int[npre+1:end]-quiescent_value) ;rtol=1e-3)# pulse average deviation from pretrig_mean
  @test isapprox(pulse_rms, sqrt(mean((pulse_int[npre+1:end]-quiescent_value).^2));rtol=1e-3) # pulse rms relative to pretrig_rms
  @test rise_time ≈ frame_time*(find(pulse.>10000)[1]-find(pulse.>2000)[1]) # defined as 10%-90% rise time
  @test postpeak_deriv ≈ maximum(diff(pulse_int)[peak_index:end]) # for a noiseless pulse, the postpeak deriv should be negative
  @test peak_index ≈ indmax(pulse) # index of the maximum value
  @test isapprox(peak_value+pretrig_mean, maximum(pulse) ;rtol=1e-3) #peak value is defined relative to pretrig_mean
  @test min_value ≈ minimum(pulse) #min_value is the minimum value, not defined relative to anything
end


text:
import StatsBase.countmap
using DataFrames
import Combinatorics.combinations

macro no_gc(ex)
    quote
        GC.enable(false)
        local val = $(esc(ex))
        while !GC.enable(true)
            nothing
        end
        GC.gc()
        val
    end
end

function get_unique_items(transactions::Array{Array{T,1},1}) where {T}
    unique_item_list = unique(vcat(transactions...))
    unique_items = Dict(zip(unique_item_list, range(1, stop=length(unique_item_list))))
end

function get_sorted_integer_transactions(transactions::Array{Array{T,1},1}, unique_items::Dict{T,Int}) where {T}
    sort.(map(x->map(y->unique_items[y],x), transactions))
end

function get_support_dict(transactions::Array{Array{Int,1},1}, item_set_size::Int, min_support::Int) where {T}
    if item_set_size == 1
        @no_gc item_sets = map(x->[x], vcat(transactions...))
    else
        @no_gc item_sets = vcat(collect.(combinations.(transactions,item_set_size))...)
    end
    support_dict = countmap(item_sets)
    return filter(x -> x[2] >= min_support, support_dict)
end

function get_support_dictionaries(transactions::Array{Array{Int,1},1}, max_length::Int, min_support::Int) where {T}
    support_dictionaries = Array{Dict{Array{Int,1},Int},1}()
    for item_set_size in range(1, stop=max_length)
        push!(support_dictionaries, get_support_dict(transactions,item_set_size, min_support))
    end
    return support_dictionaries
end

function frequent(transactions::Array{Array{T,1},1};max_length::Int,min_support::Int) where {T}
    unique_items = get_unique_items(transactions)
    integer_transactions = get_sorted_integer_transactions(transactions,unique_items)
    support_dictionaries = get_support_dictionaries(integer_transactions,max_length,min_support)

    unique_item_list = map(x->x[1],sort(collect(unique_items), by=x->x[2]))
    df=DataFrame(map(x->(unique_item_list[x[1]],x[2]), vcat(collect.(support_dictionaries)...)))
    rename!(df,[:set,:support])
end

function get_rules(support_dictionaries::Array{Dict{Array{Int,1},Int},1}, unique_items::Dict{T,Int}) where {T}
    rules = Array{Tuple{Array{T,1}, Array{T,1}, Int, Int, Float64},1}()
    unique_item_list_sorted = map(x->x[1], sort(collect(unique_items), by=x->x[2]))
    for item_set_size in range(2, stop=length(support_dictionaries))
        println("generating rules for item-set-size: $item_set_size  ")
        for (itemset, support_suggestion) in support_dictionaries[item_set_size]
            for subset in collect(combinations(itemset))[1:end-1]
                suggestion = setdiff(itemset, subset)
                support_subset = support_dictionaries[length(subset)][subset]
                confidence = support_suggestion / support_subset
                push!(rules, (unique_item_list_sorted[subset], unique_item_list_sorted[suggestion], support_subset, support_suggestion, confidence))
            end
        end
    end
    rules
end

function get_rules_df(support_dictionaries::Array{Dict{Array{Int,1},Int64},1}, min_confidence, unique_items)
    rules = get_rules(support_dictionaries, unique_items)
    rules_df = DataFrame(filter(x-> x[5] >= min_confidence, rules))
    col_names = [:subset, :suggestion, :support_subset, :support_suggestion, :confidence]
    rename!(rules_df, col_names)
end

function apriori(transactions::Array{Array{T,1},1}; max_length::Int, min_support::Int, min_confidence::Union{Int,Float64}) where {T}
    unique_items = get_unique_items(transactions)
    integer_transactions = get_sorted_integer_transactions(transactions, unique_items)
    support_dictionaries = get_support_dictionaries(integer_transactions, max_length, min_support)
    rules_df = get_rules_df(support_dictionaries, min_confidence, unique_items)
end


text:
using BinaryBuilder

name = "ncurses"
version = v"6.1"

sources = [
    "https://invisible-mirror.net/archives/ncurses/ncurses-6.1.tar.gz" =>
        "aa057eeeb4a14d470101eff4597d5833dcef5965331be3528c08d99cebaa0d17",
]

script = raw"""
# tic for host
cd ${WORKSPACE}/srcdir/ncurses-*
./configure CC="$CC_FOR_BUILD" CXX="$CXX_FOR_BUILD" RANLIB=$RANLIB_FOR_BUILD LD=$LD_FOR_BUILD CFLAGS="" LDFLAGS="" --without-shared --with-normal --disable-mixed-case
make -C include
make -C progs tic
cp progs/tic /usr/bin

make distclean
./configure --prefix=${prefix} --host=${target} --build=x86_64-linux-gnu --with-shared --without-progs --with-build-cc="$CC_FOR_BUILD" --with-termlib --without-progs --without-tack --without-tests --without-manpages --without-normal --disable-mixed-case
make -j${nproc}
make install

# avoid case-insensitive file system problems...
rm -f $prefix/share/terminfo/*/hp70092A
rm -f $prefix/share/terminfo/*/2621a
rm -f $prefix/share/terminfo/*/hp2621a
"""

products(prefix) = [
    LibraryProduct(prefix, "libncurses", :libncurses),
    LibraryProduct(prefix, "libtinfo", :libtinfo)
]

platforms = filter(x->!isa(x,Windows),supported_platforms())

dependencies = [
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
function GenerateExamples(example::Char)
    # Generate damping examples of dimension n=801, 1601, 2001
    # with example=='S','L'or 'H', for Small, Large and Homogeneous, respectively.

    if example=='S'
        n₁=400
    elseif example=='L'
        n₁=800
    else
        n₁=1000
    end

    # Dimension of the system
    n=2*n₁+1
    m₁=1
    m₂=1

    if example=='S'
        m₀=1200.0
    elseif example=='L'
        m₀=1800.0
    else
        m₀=2000.0
    end

    k₁=100
    k₂=150
    k₃=200

    first_part=Array{Int,1}(undef,n₁)
    second_part=Array{Int,1}(undef,n₁)

    if example=='S'
        for ii=1:400
            first_part[ii]=1000-4*ii
        end
        for ii=201:400
            first_part[ii]=3*ii-400
        end
        for ii=1:400
            second_part[ii]=500+ii
        end
    elseif example=='L'
        for ii=1:div(n₁,2)
            first_part[ii]=2000-4*ii
        end
        for ii=div(n₁,2)+1:n₁
            first_part[ii]=3*ii-800
        end
        for ii=1:n₁
            second_part[ii]=500+ii
        end
    else
        for ii=1:n₁
           first_part[ii]=1000
        end
        for ii=1:n₁
          second_part[ii]=1500
        end
    end

    M=Diagonal([first_part; second_part; m₀])

    K₁=SymTridiagonal(2*ones(n₁),-ones(n₁-1))
    α=0.02
    K₁=Matrix(K₁)

    K=zeros(2*n₁+1,2*n₁+1)
    K[1:n₁,1:n₁]=k₁*K₁
    K[n₁+1:2*n₁,n₁+1:2*n₁]=k₂*K₁

    K[1:n₁,end]=-[zeros(1,n₁-1) k₁]'
    K[n₁+1:2*n₁,end]=-[zeros(1,n₁-1) k₂]'

    K[end,1:n₁]=-[zeros(1,n₁-1) k₁]
    K[end,n₁+1:2*n₁]=-[zeros(1,n₁-1) k₂]
    K[2*n₁+1, 2*n₁+1]=k₁+k₂+k₃

    # Return stiffness and mass matrix
    return K, M
end


text:
using PyCall

# i = 1, 2, ...
up_index(i) = 2*(i-1) + 1
down_index(i) = 2*i

@enum PauliID pauli_I=0 pauli_X=1 pauli_Y=2 pauli_Z=3
pauli_id_lookup = Dict("I"=>pauli_I, "X"=>pauli_X, "Y"=>pauli_Y, "Z"=>pauli_Z)
pauli_id_str = ["I", "X", "Y", "Z"]

################################################################################
############################# QUANTUM  CIRCUIT #################################
################################################################################
include("gate.jl")

################################################################################
############################# QUANTUM STATE  ###################################
################################################################################
include("quantum_state.jl")

################################################################################
############################# QUANTUM  CIRCUIT #################################
################################################################################
include("circuit.jl")

################################################################################
################################# OPERATOR #####################################
################################################################################
abstract type SecondQuantOperator end

function Base.:(==)(op1::SecondQuantOperator, op2::SecondQuantOperator)
    op1.pyobj == op2.pyobj
end

function Base.:*(op1::SecondQuantOperator, op2::SecondQuantOperator)
    typeof(op1)(op1.pyobj * op2.pyobj)
end

function Base.:+(op1::SecondQuantOperator, op2::SecondQuantOperator)
    typeof(op1)(op1.pyobj + op2.pyobj)
end

function Base.:-(op1::SecondQuantOperator, op2::SecondQuantOperator)
    typeof(op1)(op1.pyobj - op2.pyobj)
end

function Base.:/(op::SecondQuantOperator, x::Number)
    typeof(op)(op.pyobj/x)
end

################################################################################
############################# FERMION OPERATOR #################################
################################################################################
# Wrap openfermion.ops.operators.FermionOperator
struct FermionOperator <: SecondQuantOperator
    pyobj::PyObject
end

function _parse_fermion_operator_str(str::String)
    regex = r"\d+\^*"
    res = Tuple{Int64,Int64}[]
    for x in eachmatch(regex, str)
        m = x.match
        op_type = m[end] == '^' ? 1 : 0
        idx = (op_type == 1 ? parse(Int64, m[1:end-1]) : parse(Int64, m))
        push!(res, (idx, op_type))
    end
    res
end

function FermionOperator(ops::Vector{Tuple{Int,Int}}, coeff::Number=1.0)
    for (idx, _) in ops
        if idx <= 0
            error("idx for fermion operator must be positive.")
        end
    end
    ops_py = [(idx-1, op_type) for (idx, op_type) in ops]
    FermionOperator(ofermion.ops.operators.FermionOperator(Tuple(ops_py), coeff))
end

function FermionOperator()
    FermionOperator("", 0.0)
end

function FermionOperator(op_str::String, coeff::Number=1.0)
    FermionOperator(_parse_fermion_operator_str(op_str), coeff)
end

function get_number_preserving_sparse_operator(ham::FermionOperator, n_qubit::Int, n_electron::Int)::PyObject
    ofermion.linalg.get_number_preserving_sparse_operator(ham.pyobj, n_qubit, n_electron)
end

function jordan_wigner(op::FermionOperator)
    OFQubitOperator(ofermion.transforms.jordan_wigner(op.pyobj))
end

################################################################################
############################# QUBIT OPERATOR ###################################
################################################################################
abstract type QubitOperator <: SecondQuantOperator end

# Wrap openfermion.ops.operators.qubit_operator.QubitOperator
struct OFQubitOperator <: QubitOperator
    pyobj
end

function _convert_qubitop_str_from_py_to_jl(str::String)
    regex = r"[XYZI]\d+"
    res = ""
    for x in eachmatch(regex, str)
        m = x.match
        pauli_op_char = m[1]
        idx = parse(Int64, m[2:end])
        if idx <= 0
            error("idx for qubit must be positive!")
        end
        res *= " $pauli_op_char$(idx-1)"
    end
    res
end

function OFQubitOperator(str::String, coeff::Number=1.0)
    return OFQubitOperator(ofermion.ops.operators.qubit_operator.QubitOperator(
        _convert_qubitop_str_from_py_to_jl(str), coeff))
end

function hermitian_conjugated(op::OFQubitOperator)
    OFQubitOperator(ofermion.utils.hermitian_conjugated(op.pyobj))
end

function get_n_qubit(op::OFQubitOperator)
    count_qubit_in_qubit_operator(op.pyobj)
end

function get_term_count(op::OFQubitOperator)
    length(op.pyobj.terms)
end

function terms_dict(op::OFQubitOperator)::Dict{Any,Any}
    # Change index convention
    dict_jl = Dict()
    for (k, v) in op.pyobj.terms
        k_new = Tuple(( (k_[1]+1, k_[2]) for k_ in k))
        dict_jl[k_new] = v
    end
    dict_jl
end

function is_hermitian(op::OFQubitOperator)
    ofermion.utils.operator_utils.is_hermitian(op.pyobj)
end

function convert_to_qulacs_op(op::OFQubitOperator, n_qubit::Int)
    convert_openfermion_op(n_qubit, op.pyobj)
end

function get_expectation_value(op::OFQubitOperator, state::QulacsQuantumState)
    if !is_hermitian(op)
        error("op must be hermite for get_expectation_value!")
    end
    convert_to_qulacs_op(op, get_n_qubit(state)).get_expectation_value(state.pyobj)
end

function get_transition_amplitude(op::OFQubitOperator, state_bra::QulacsQuantumState, state_ket::QulacsQuantumState)
    convert_to_qulacs_op(op, get_n_qubit(state_bra)).get_transition_amplitude(state_bra.pyobj, state_ket.pyobj)
end

function inner_product(state_bra::QulacsQuantumState, state_ket::QulacsQuantumState)
    qulacs.state.inner_product(state_bra.pyobj, state_ket.pyobj)
end

"""
Parse a tuple representing a Pauli string
  When x is ((1, "X"), (5, "Y")), returns [1, 5], [PauliID.X, PauliID.Y]
"""
function parse_pauli_str(x)::Tuple{Vector{Int64}, Vector{PauliID}}
    collect(e[1] for e in x), collect(pauli_id_lookup[e[2]] for e in x)
end

"""
Update a state using a circuit
"""
function update_quantum_state!(circuit::QulacsParametricQuantumCircuit, state::QulacsQuantumState)
    circuit.pyobj.update_quantum_state(state.pyobj)
end

function update_quantum_state!(circuit::QulacsQuantumCircuit, state::QulacsQuantumState)
    circuit.pyobj.update_quantum_state(state.pyobj)
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
export Trajectory,trajectory,
       kspaceNodes, echoTime, acqTimePerProfile, readoutTimes,
       numSamplingPerProfile, numProfiles, numSlices, findCenters, samplingFreq,
       scale

import Base: string

"""
struct describing a trajectory

# Fields
* `name::String`                  - name of the trajectory
* `nodes::Matrix{Float64}`        - sampling locations in k-space.
                                    (1.dim <-> dimensions of k-space, 2. dim <-> sampling points)
* `times::Vector{Float64}`        - sampling times in s
* `TE::Float64`                   - echo time in s
* `AQ::Float64``                  - readout duration in s (per profile)
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* `numSlices::Int64`              - number of slices (for 3d trajectories)
* `cartesian::Bool`               - true if sampling points lie on a cartesian grid
* `circular::Bool`                - true if kspace is covered in a circular domain
"""
mutable struct Trajectory
  name::String
  nodes::Matrix{Float64}
  times::Vector{Float64}
  TE::Float64
  AQ::Float64
  numProfiles::Int64
  numSamplingPerProfile::Int64
  numSlices::Int64
  cartesian::Bool
  circular::Bool
end

function Trajectory(nodes::AbstractMatrix{T}, numProfiles::Int64, numSamplingPerProfile;
                    times=nothing, TE::Float64=0.0, AQ::Float64=1.e-3, numSlices::Int64=1,
                    cartesian::Bool=false, circular::Bool=false) where T <: AbstractFloat
  if times != nothing
    ttimes = times
  else
    ttimes = readoutTimes(numProfiles,numSamplingPerProfile; TE=TE, AQ=AQ)
  end

  return Trajectory("Custom", Float64.(nodes), ttimes, TE, AQ, numProfiles, numSamplingPerProfile, numSlices, cartesian, circular)
end

Base.vec(tr::Trajectory) = [tr]
Base.vec(tr::Vector{Trajectory}) = tr
Base.size(tr::Trajectory) = size(kspaceNodes(tr))
Base.size(tr::Trajectory,i::Int) = size(kspaceNodes(tr),i)
function scale(tr::Trajectory, factor::Real)
  trC = deepcopy(tr)
  trC.nodes ./= factor
  return trC
end

include("2D/Cartesian2D.jl")
include("2D/Radial2D.jl")
include("2D/Spiral2D.jl")
include("2D/OneLine2D.jl")
include("2D/Spiral2DVariableDens.jl")
include("2D/Spiral2DDualDens.jl")
include("2D/Spiral2DPerturbed.jl")
include("2D/EPI.jl")

include("3D/Kooshball.jl")
include("3D/Cartesian3D.jl")
include("3D/StackofStars.jl")

#### Factory method ###

# This dispatches on the file extension and automatically
# generates the correct type
"""
    trajectory(trajName::AbstractString, numProfiles::Int, numSamplingPerProfile::Int; numSlices::Int64=1, TE::Float64=0.0, AQ::Float64=1.e-3, kargs...)

is a factory method to construct a trajectory from its `name`

# Arguments
* `name::String`                  - name of the trajectory
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`numSlices::Int64=1`)          - number of slices (for 3d trajectories)
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* `kargs...`                      - addional keyword arguments
"""
function trajectory(trajName::AbstractString, numProfiles::Int, numSamplingPerProfile::Int; numSlices::Int64=1, TE::Float64=0.0, AQ::Float64=1.e-3, kargs...)
  if trajName == "Spiral"
    tr = SpiralTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "Radial"
    tr = RadialTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "Cartesian"
    tr = CartesianTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "EPI"
    tr = EPITrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "OneLine"
    tr = OneLine2dTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "SpiralVarDens"
    tr = SpiralTrajectoryVarDens(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
elseif trajName == "SpiralPerturbed"
    tr = SpiralPerturbedTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "SpiralDualDens"
    tr = SpiralTrajectoryDualDens(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  elseif trajName == "Cartesian3D"
    tr = CartesianTrajectory3D(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, numSlices=numSlices, kargs...)
  elseif trajName == "StackOfStars"
    tr = StackOfStarsTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, numSlices=numSlices, kargs...)
  elseif trajName == "Kooshball"
    tr = KooshballTrajectory(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ, kargs...)
  else
    @error "The trajectory $trajName is not yet supported!"
  end

  # return Trajectory(trajName, nodes, TE, AQ, numProfiles, numSamplingPerProfile, numSlices, false, true)
  return tr
end

## some general methods that might be overloaded by the real implementations
""" `string(tr::Trajectory)` returns the name of a trajectory """
string(tr::Trajectory) = tr.name

""" `echoTime(tr::Trajectory)` returns the echo time of a trajectory """
echoTime(tr::Trajectory) = tr.TE

""" `acqTimePerProfile(tr::Trajectory)` returns the acquisition time per profile of a trajectory """
acqTimePerProfile(tr::Trajectory) = tr.AQ

""" `string(tr::Trajectory)` returns the name of a trajectory """
numProfiles(tr::Trajectory) = tr.numProfiles

""" `numSamplingPerProfile(tr::Trajectory)` returns the number of samples per profile of a trajectory """
numSamplingPerProfile(tr::Trajectory) = tr.numSamplingPerProfile

""" `numSlices(tr::Trajectory)` returns the number of slices of a trajectory """
numSlices(tr::Trajectory) = tr.numSlices

""" `isCircular(tr::Trajectory)` returns whether the trajectory has circular k-space coverage """
isCircular(tr::Trajectory) = tr.circular

""" `isCartesian(tr::Trajectory)` returns whether the trajectory nodes lie on a cartesian grid"""
isCartesian(tr::Trajectory) = tr.cartesian

# 2d or 3d encoding
""" `dims(tr::Trajectory)` returns the number of dimensions of a trajectory """
dims(tr::Trajectory) = size(tr.nodes,1)

""" `kspaceNodes(tr::Trajectory)` returns the kspace sampling points of a trajectory """
kspaceNodes(tr::Trajectory) = tr.nodes

""" `readoutTimes(tr::Trajectory)` returns the readoutTimes for the sampling points of a trajectory """
readoutTimes(tr::Trajectory) = tr.times

""" `profileNodes(tr::Trajectory,prof::Int,slice::Int)` returns kspace nodes for a given profile """
function profileNodes(tr::Trajectory,prof::Int,slice::Int)
  nodes = reshape(tr.nodes,:,tr.numSamplingPerProfile,tr.numProfiles, tr.numSlices)
  return nodes[:,:,prof,slice]
end

function readoutTimes(numProfiles, numSamplingPerProfile, numSlices=1; TE=0.0, AQ=1.e-3)
  times = zeros(numSamplingPerProfile, numProfiles, numSlices)
  for j = 1:numSlices
    for l = 1:numProfiles
      for k = 1:numSamplingPerProfile
        times[k,l,j] = TE + AQ*(k-1)/numSamplingPerProfile
      end
    end
  end
  return vec(times)
end

function zax(n::Int64,N::Int64)
  return (2 .*n - N -1.)/N
end

function xax(n::Int64,N::Int64)
  return cos(sqrt(N*pi)*asin(zax(n,N)))*sqrt(1.0-zax(n,N)^2.)
end

function yax(n::Int64,N::Int64)
  return sin(sqrt(N*pi)*asin(zax(n,N)))*sqrt(1.0-zax(n,N)^2.)
end

# function findCenters(tr::Abstract2DTrajectory)
#   knodes = kspaceNodes(tr)
#   temp = find(x->x==0,knodes[1,:])
#
#   zeroInd=Int64[]
#
#   for i in temp
#
#     if knodes[2,i] == 0
#       push!(zeroInd,i)
#     end
#
#   end
#
#   return zeroInd
#
# end
#
# function findCenters(tr::Trajectory)
#   knodes = kspaceNodes(tr)
#   temp = find(x->x==0,knodes[1,:])
#
#   zeroIndY=Int64[]
#   zeroInd = Int64[]
#
#   for i in temp
#     if knodes[2,i] == 0
#       push!(zeroIndY,i)
#     end
#   end
#
#   for i in zeroIndY
#     if knodes[3,i] == 0
#       push!(zeroInd,i)
#     end
#   end
#
#   return zeroInd
#
# end
#
# function findCenters(tr::StackOfStarsTrajectory)
#   knodes = kspaceNodes(tr)
#   temp = find(x->x==0,knodes[1,:])
#
#   zeroInd = Int64[]
#
#   for i in temp
#     if knodes[2,i] == 0
#       push!(zeroInd,i)
#     end
#   end
#
#   return zeroInd
#
# end
#
# function findCenterNeighbours(tr::CartesianTrajectory,range::Int64)
#   knodes = kspaceNodes(tr)
#   pRange = minimum(knodes[1,find(x->x>0,knodes[1,:])])
#   mRange = maximum(knodes[1,find(x->x<0,knodes[1,:])])
#
#   idx=Int64[]
#
#   zeroInd = findCenters(tr)
#
#   push!(idx,zeroInd[1])
#
#
#   for j=1:range
#
#     temp = find(x->x==0,knodes[1,:])
#
#     for i in temp
#       for k=range*mRange:pRange:range*pRange
#         if k!=0
#           if knodes[2,i] == k
#             push!(idx,i)
#           end
#         end
#       end
#     end
#
#     temp = find(x->x==j*pRange,knodes[1,:])
#
#     for i in temp
#       for k=range*mRange:pRange:range*pRange
#         if knodes[2,i] == k
#           push!(idx,i)
#         end
#       end
#     end
#
#     temp = find(x->x==j*mRange,knodes[1,:])
#
#     for i in temp
#       for k=range*mRange:pRange:range*pRange
#         if knodes[2,i] == k
#           push!(idx,i)
#         end
#       end
#     end
#
#   end
#
#   return idx
#
# end

function improvedTrajPlot(tr::Trajectory)
  numProfiles = numProfiles(tr)
  numSamplePerProfile = numSamplingPerProfile(tr)
  nodes = kspaceNodes(tr)
  figure()
  for i=1:numProfiles
     indices = 1 + (i-1)*numSamplePerProfile : i*numSamplePerProfile
     plot(nodes[1,indices],nodes[2,indices],color="blue" )
     plot(nodes[1,indices],nodes[2,indices],".",color="red", )
  end
end


text:
#=
This code is loosely based on
Strid and Walentin (2008) "Block Kalman filtering for large-scale DSGE models"
=#

mutable struct BlockKalmanFilter{S<:AbstractFloat}
    A1::Matrix{S}
    A2::Matrix{S}
    A3::Matrix{S}
    B1::Matrix{S}
    B2::Matrix{S}
    B3::Matrix{S}
    Cblock::Matrix{S}
    Rup::Matrix{S}
    Rlo::Matrix{S}
    C::Vector{S}
    Q::Matrix{S}
    Z1::Matrix{S}
    Z2::Matrix{S}
    Z3::Matrix{S}
    Z4::Matrix{S}
    D::Vector{S}
    E::Matrix{S}
    s1_t::Vector{S} # s1_{t|t-1} or s1_{t|t}, and so on
    s2_t::Vector{S}
    s3_t::Vector{S}
    s4_t::Vector{S}
    P11_t::Matrix{S} # P11_{t|t-1} or P11_{t|t}
    P12_t::Matrix{S}
    P13_t::Matrix{S}
    P14_t::Matrix{S}
    P22_t::Matrix{S}
    P23_t::Matrix{S}
    P24_t::Matrix{S}
    P33_t::Matrix{S}
    P34_t::Matrix{S}
    P44_t::Matrix{S}
    M::Matrix{S}
    Mtild::Matrix{S}
    block_dims::Vector{Int64}
    block_num::Int64
    true_block::Bool
    loglh_t::S     # P(y_t | y_{1:t})
end

"""
```
KalmanFilter(T, R, C, Q, Z, D, E, [s_0, P_0])
```

Outer constructor for the `KalmanFilter` type.
"""
function BlockKalmanFilter(Ttild::Matrix{S}, Rtild::Matrix{S}, Ctild::Vector{S}, Qtild::Matrix{S},
                      Ztild::Matrix{S}, Dtild::Vector{S}, Etild::Matrix{S}, M::Matrix{S}, Mtild::Matrix{S},
                           block_dims::Vector{Int64}, s_0tild::Vector{S} = Vector{S}(0),
                           P_0tild::Matrix{S} = Matrix{S}(0, 0); block_num::Int64 = 2, true_block::Bool = false) where {S<:AbstractFloat}
    if isempty(s_0tild) || isempty(P_0tild)
        s_0tild, P_0tild = init_stationary_states(Ttild, Rtild, Ctild, Qtild)
    end
    T = M' * Ttild * M
    R = M' * Rtild * Mtild
    C = M * Ctild
    Q = Mtild' * Qtild * Mtild
    Z = Ztild * M
    # D and E are left as is
    s_0 = M * s_0tild
    P_0 = M' * P_0tild * M
    if block_dims[2] == 0 & block_dims[3] == 0
        # 2-block case
        dim1 = block_dims[1]; dim2 = block_dims[4]
        A1 = T[1:dim1, 1:dim1]
        B1 = R[dim1 + 1:end, 1:dim1]
        Cblock = T[dim1 + 1:end, dim1 + 1:end]
        Z1 = Matrix{S}(0, 0)
        Z2 = Z[:, dim1 + 1:end]
        Rup = Matrix{S}(0, 0)
        Rlo = Matrix{S}(0, 0)

        # Check if dim1 is equal to number of exogenous shocks; if not, need to track more matrices
        if dim1 != size(R, 2)
            true_block = false
            Rup = R[1:dim1, dim1 + 1:end]
            Rlo = R[dim1 + 1:end, dim1 + 1:end]
        end
        # Check Z is in true block form; if not (e.g. unit root), need to track more matrices
        if Z[:, 1:dim1] != zeros(Z[:, 1:dim1])
            true_block = false
            Z1 = Z[:, 1:dim1]
        end
        s1_0 = s_0[1:dim1]
        s2_0 = s_0[dim1 + 1:end]
        P11_0 = P_0[1:dim1, 1:dim1]
        P12_0 = P_0[1:dim1, dim1 + 1:end]
        P22_0 = P_0[dim1 + 1:end, dim1 + 1:end]

        return BlockKalmanFilter(A1, Matrix{S}(0, 0), Matrix{S}(0, 0), B1, Matrix{S}(0, 0), Matrix{S}(0, 0),
                          Cblock, Rup, Rlo, C, Q, Z1, Z2, Matrix{S}(0, 0),
                          Matrix{S}(0, 0), D, E, s1_0, s2_0, Vector{S}(0), Vector{S}(0),
                          P11_0, P12_0, Matrix{S}(0, 0), Matrix{S}(0, 0), P22_0,
                          Matrix{S}(0, 0), Matrix{S}(0, 0), Matrix{S}(0, 0), Matrix{S}(0, 0),
                          Matrix{S}(0, 0), M, Mtild, block_dims, block_num, true_block, NaN)
    end
end

"""
```
kalman_filter(y, T, R, C, Q, Z, D, E, s_0 = Vector(), P_0 = Matrix();
    outputs = [:loglh, :pred, :filt], Nt0 = 0)

kalman_filter(regime_indices, y, Ts, Rs, Cs, Qs, Zs, Ds, Es,
    s_0 = Vector(), P_0 = Matrix(); outputs = [:loglh, :pred, :filt],
    Nt0 = 0)
```

This function implements the Kalman filter for the following state-space model:

```
s_{t+1} = C + T*s_t + R*ϵ_t    (transition equation)
y_t     = D + Z*s_t + u_t      (measurement equation)

ϵ_t ∼ N(0, Q)
u_t ∼ N(0, E)
Cov(ϵ_t, u_t) = 0
```

### Inputs

- `y`: `Ny` x `Nt` matrix containing data `y_1, ... , y_T`
- `s_0`: optional `Ns` x 1 initial state vector
- `P_0`: optional `Ns` x `Ns` initial state covariance matrix

**Method 1 only:**

- `T`: `Ns` x `Ns` state transition matrix
- `R`: `Ns` x `Ne` matrix in the transition equation mapping shocks to states
- `C`: `Ns` x 1 constant vector in the transition equation
- `Q`: `Ne` x `Ne` matrix of shock covariances
- `Z`: `Ny` x `Ns` matrix in the measurement equation mapping states to
  observables
- `D`: `Ny` x 1 constant vector in the measurement equation
- `E`: `Ny` x `Ny` matrix of measurement error covariances

**Method 2 only:**

- `regime_indices`: `Vector{Range{Int}}` of length `n_regimes`, where
  `regime_indices[i]` indicates the time periods `t` in regime `i`
- `Ts`: `Vector{Matrix{S}}` of `T` matrices for each regime
- `Rs`
- `Cs`
- `Qs`
- `Zs`
- `Ds`
- `Es`

where:

- `Nt`: number of time periods for which we have data
- `Ns`: number of states
- `Ne`: number of shocks
- `Ny`: number of observables

### Keyword Arguments

- `outputs`: some subset of `[:loglh, :pred, :filt]` specifying which outputs to
  compute and return. There will always be the same number of return values,
  but, for example, `s_pred` and `P_pred` will be returned as empty arrays if
  `:pred` is not in `outputs`
- `Nt0`: number of presample periods to omit from all return values

### Outputs

- `loglh`: length `Nt` vector of conditional log-likelihoods P(y_t | y_{1:t-1})
- `s_pred`: `Ns` x `Nt` matrix of one-step predicted state vectors s_{t|t-1}
- `P_pred`: `Ns` x `Ns` x `Nt` array of mean squared errors P_{t|t-1} of
  predicted state vectors
- `s_filt`: `Ns` x `Nt` matrix of filtered state vectors s_{t|t}
- `P_filt`: `Ns` x `Ns` x `Nt` matrix containing mean squared errors P_{t|t} of
  filtered state vectors
- `s_0`: `Ns` x 1 initial state vector. This may have been reassigned to the
  last presample state vector if `Nt0 > 0`
- `P_0`: `Ns` x `Ns` initial state covariance matrix. This may have been
  reassigned to the last presample state covariance if `Nt0 > 0`
- `s_T`: `Ns` x 1 final filtered state `s_{T|T}`
- `P_T`: `Ns` x `Ns` final filtered state covariance matrix `P_{T|T}`

### Notes

When `s_0` and `P_0` are omitted, they are computed using
`init_stationary_states`.
"""
function block_kalman_filter(y::Matrix{Float64}, Ttild::Matrix{Float64}, Rtild::Matrix{Float64},
                             Ctild::Vector{Float64}, Qtild::Matrix{Float64}, Ztild::Matrix{Float64},
                             Dtild::Vector{Float64}, Etild::Matrix{Float64},
                             M::Matrix{Float64}, Mtild::Matrix{Float64}, block_dims::Vector{Int64},
                             s_0tild::Vector{Float64} = Vector{Float64}(0),
                             P_0tild::Matrix{Float64} = Matrix{Float64}(0,0);
                             outputs::Vector{Symbol} = [:loglh, :pred, :filt],
                             Nt0::Int = 0, block_num::Int64 = 2, true_block::Bool = false)

    # Determine outputs
    return_loglh = :loglh in outputs
    return_pred  = :pred in outputs
    return_filt  = :filt in outputs

    # Dimensions
    Ns = size(T,1) # number of states
    Nt = size(y,2) # number of periods of data

    # Initialize Inputs and outputs, populate initial states
    s_0 = M * s_0tild
    P_0 = M' * P_0tild * M
    k = BlockKalmanFilter(Ttild, Rtild, Ctild, Qtild, Ztild, Dtild, Etild, M, Mtild,
                          block_dims, s_0tild, P_0tild; block_num = block_num, true_block = true_block)

    mynan  = convert(Float64, NaN)
    loglh  = return_loglh ? fill(mynan, Nt)         : Vector{Float64}(0)
    s_pred = return_pred  ? fill(mynan, Ns, Nt)     : Matrix{Float64}(0, 0)
    P_pred = return_pred  ? fill(mynan, Ns, Ns, Nt) : Array{Float64, 3}(0, 0, 0)
    s_filt = return_filt  ? fill(mynan, Ns, Nt)     : Matrix{Float64}(0, 0)
    P_filt = return_filt  ? fill(mynan, Ns, Ns, Nt) : Array{Float64, 3}(0, 0, 0)

    # Loop through periods t
    for t = 1:Nt
        # Forecast
        forecast!(k)
        if return_pred
            if k.block_num == 2
                s_pred[:,    t] = vcat(k.s1_t, k.s2_t)
                P_pred[:, :, t] = [k.P11_t k.P12_t; k.P12_t' k.P22_t]
            end
        end

        # Update and compute log-likelihood
        update!(k, y[:, t]; return_loglh = return_loglh)
        if return_filt
            if k.block_num == 2
                s_filt[:,    t] = vcat(k.s1_t, k.s2_t)
                P_filt[:, :, t] = [k.P11_t k.P12_t; k.P12_t' k.P22_t]
            end
        end
        if return_loglh
            loglh[t]        = k.loglh_t
        end

        # # Update s_0 and P_0 if Nt_0 > 0, presample periods
        # if t == Nt0
        #     s_0 = k.s_t
        #     P_0 = k.P_t
        # end
    end

    # Populate final states
    if k.block_num == 2
        s_T = vcat(k.s1_t, k.s2_t)
        P_T = [k.P11_t k.P12_t; k.P12_t' k.P22_t]
    end

    # # Remove presample periods
    # loglh, s_pred, P_pred, s_filt, P_filt =
    #     remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt; outputs = outputs)

    return loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T
end

# Computes the one-step-ahead states s_{t|t-1} and state covariances P_{t|t-1}
# and assign to 'k'
function forecast!(k::BlockKalmanFilter{Float64})

    if k.block_num == 2
        # Get block matrices that don't depend on true_block
        dim1 = k.block_dims[1] # dimension of exogenous AR(1) processes
        dim2 = k.block_dims[4] # dimension of endogenous states
        s1_filt = k.s1_t
        s2_filt = k.s2_t
        P11_filt = k.P11_t
        P12_filt = k.P12_t
        P22_filt = k.P22_t
        Ablock = k.A1
        Bblock = k.B1
        Cblock = k.Cblock
        Q1 = k.Q[1:dim1, 1:dim1]
        Q2 = k.Q[dim1 + 1:end, dim1 + 1:end]

        # Compute P_t, s_t
        if k.true_block
            L_t = (P12_filt * Cblock') .* (diag(Ablock) * ones(1, dim2))
            P11_t = P11_filt .* (diag(Ablock) * diag(Ablock)') + Q1
            P12_t = P11_t * Bblock' + L_t
            G = Bblock * (P12_t + L_t)
            P22_t = (G + G')/2 + Cblock * P22_filt * Cblock'

            s1_t = diag(Ablock) .* s1_filt
            s2_t = Bblock * k.s1_t + Cblock * s2_filt
        else
            # RQR' = |  Q1 + Rup * Q2 * Rup', Q1*Bblock' + Rup * Q2 * Rlo'        |
            #        |     ---           , Bblock * Q1 * Bblock' + Rlo * Q2 * Rlo'|
            L_t = (P12_filt * Cblock') .* (diag(Ablock) * ones(1, dim2))
            RupQ2 = k.Rup * Q2

            # Compute intermediate matrices to minimize additions
            P11_t = P11_filt .* (diag(Ablock) * diag(Ablock)') + Q1
            P12_t = P11_t * Bblock' + L_t
            G = Bblock * (P12_t + L_t)

            # Compute final matrices
            P22_t = (G + G')/2 + k.Rlo * Q2 * k.Rlo' + Cblock * P22_filt * Cblock'
            P12_t = P12_t + RupQ2 * k.Rlo'
            P11_t = P11_t + RupQ2 * k.Rup'

            s1_t = diag(Ablock) .* s1_filt
            s2_t = Bblock * s1_t + Cblock * s2_filt
        end

        # Save values
        k.s1_t = s1_t; k.s2_t = s2_t; k.P11_t = P11_t; k.P12_t = P12_t; k.P22_t = P22_t
    end
    return nothing
end

function update!(k::BlockKalmanFilter{Float64}, y_obs::Vector{Float64}; return_loglh::Bool = true)
    # Keep rows of measurement equation corresponding to non-NaN observables
    nonnan = .!isnan.(y_obs)
    y_obs = y_obs[nonnan]
    D = k.D[nonnan]
    E = k.E[nonnan, nonnan]
    Ny = length(y_obs)

    if k.block_num == 2
        # Get block matrices
        dim1 = k.block_dims[1] # dimension of exogenous AR(1) processes
        dim2 = k.block_dims[4] # dimension of endogenous states
        s1_pred  = k.s1_t
        s2_pred  = k.s2_t
        P11_pred = k.P11_t
        P12_pred = k.P12_t
        P22_pred = k.P22_t
        Z2 = k.Z2[nonnan, :]

        # Compute predicted y, measurement covariance matrix, error, and auxiliary matrices
        if k.true_block
            y_pred = Z2 * s2_pred + D
            V_pred = Z2 * P22_pred * Z2' + E
        else
            Z1 = k.Z1[nonnan, :]
            off_diag_mat = Z1 * P12_pred * Z2'
            y_pred = Z1 * s1_pred + Z2 * s2_pred + D
            V_pred = Z1 * P11_pred * Z1' + off_diag_mat + off_diag_mat' + Z2 * P22_pred * Z2' + E
        end
        V_pred_inv = inv(V_pred)
        dy = y_obs - y_pred # prediction error

        # Update state covariance matrix
        if k.true_block
            PZV_1 = P12_pred * Z2' * V_pred_inv
            PZV_2 = P22_pred * Z2' * V_pred_inv

            Z2_P22_pred = Z2 * P22_pred
            P11_t = P11_pred - PZV_1 * Z2 * P12_pred'
            P12_t = P12_pred - PZV_1 * Z2_P22_pred
            P22_t = P22_pred - PZV_2 * Z2_P22_pred
        else
            kalman_gain_1 = P11_pred * Z1' + P12_pred * Z2'
            kalman_gain_2 = P12_pred' * Z1' + P22_pred * Z2'
            PZV_1 = kalman_gain_1 * V_pred_inv
            PZV_2 = kalman_gain_2 * V_pred_inv

            P11_t = P11_pred - PZV_1 * kalman_gain_1'
            P12_t = P12_pred - PZV_1 * kalman_gain_2'
            P22_t = P22_pred - PZV_2 * kalman_gain_2'
        end

        # Update states
        s1_t = s1_pred + PZV_1 * dy
        s2_t = s2_pred + PZV_2 * dy

        # Save matrices
        k.s1_t = s1_t; k.s2_t = s2_t
        k.P11_t = P11_t; k.P12_t = P12_t; k.P22_t = P22_t

        if return_loglh
            k.loglh_t = -(Ny * log(2*π) + log(det(V_pred)) + dy'*V_pred_inv*dy)/2
        end
    end
end

# """
# ```
# remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt)
# ```

# Remove the first `Nt0` periods from all other input arguments and return.
# """
# function remove_presample!(Nt0::Int, loglh::Vector{S},
#                            s_pred::Matrix{S}, P_pred::Array{S, 3},
#                            s_filt::Matrix{S}, P_filt::Array{S, 3};
#                            outputs::Vector{Symbol} = [:loglh, :pred, :filt]) where {S<:AbstractFloat}
#     if Nt0 > 0
#         if :loglh in outputs
#             loglh  = loglh[(Nt0+1):end]
#         end
#         if :pred in outputs
#             s_pred = s_pred[:,    (Nt0+1):end]
#             P_pred = P_pred[:, :, (Nt0+1):end]
#         end
#         if :filt in outputs
#             s_filt = s_filt[:,    (Nt0+1):end]
#             P_filt = P_filt[:, :, (Nt0+1):end]
#         end
#     end
#     return loglh, s_pred, P_pred, s_filt, P_filt
# end