text:
# top-level bridge
# ================

"""
    mutable struct AbstractGlobal
        t::Any     # analyzed type
        iscd::Bool # is this abstract global variable 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 or alter its type on future load and store operations.

!!! 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
    t::Any     # analyzed type
    iscd::Bool # is this abstract global variable declarared as constant or not

    function AbstractGlobal(@nospecialize(t),
                            iscd::Bool,
                            )
        return new(t,
                   iscd,
                   )
    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).
"""
CC.bail_out_toplevel_call(analyzer::AbstractAnalyzer, @nospecialize(sig), sv) = false

function CC.abstract_eval_special_value(analyzer::AbstractAnalyzer, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    if istoplevel(sv)
        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
        end
    end

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

    if istoplevel(sv)
        if isa(e, GlobalRef)
            mod, name = e.mod, e.name
            if isdefined(mod, name)
                # here we 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 global variable assignments that happen in a function,
                # but it's highly possible this is a toplevel callsite and we take a risk here,
                # otherwise we can't enter the analysis !
                val = getfield(mod, name)
                ret = isa(val, AbstractGlobal) ? val.t : Const(val)
            end
        end
    end

    return ret
end

function CC.abstract_eval_value(analyzer::AbstractAnalyzer, @nospecialize(e), vtypes::VarTable, sv::InferenceState)
    ret = @invoke CC.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

    return ret
end

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

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

function CC.builtin_tfunction(analyzer::AbstractAnalyzer, @nospecialize(f), argtypes::Array{Any,1},
                              sv::InferenceState) # `AbstractAnalyzer` isn't overloaded on `return_type`
    ret = @invoke CC.builtin_tfunction(analyzer::AbstractInterpreter, f, argtypes::Array{Any,1},
                                       sv::Union{InferenceState,Nothing})

    if f === getfield && 2 ≤ length(argtypes) ≤ 3
        obj, fld = argtypes
        if isa(fld, Const)
            name = fld.val
            if isa(name, Symbol)
                if isa(obj, Const)
                    mod = obj.val
                    if isa(mod, Module)
                        if isdefined(mod, name)
                            if istoplevel_globalref(analyzer, sv)
                                # when accessing to a global variable in a module
                                # concretized by `analyzer`, eagerly propagate its type
                                # NOTE logic here should be synced with that of `abstract_eval_special_value(::AbstractAnalyzer, ::Any, ::VarTable, ::InferenceState)`
                                val = getfield(mod, name)
                                return isa(val, AbstractGlobal) ? val.t : Const(val)
                            end
                        end
                    end
                end
            end
        end
    elseif f === fieldtype
        # the valid widest possible return type of `fieldtype_tfunc` is `Union{Type,TypeVar}`
        # because fields of unwrapped `DataType`s can legally be `TypeVar`s,
        # but this will cause lots of false positive `NoMethodErrorReport`s for inference
        # with accessing to abstract fields since most methods don't expect `TypeVar`
        # (e.g. `@report_call readuntil(stdin, 'c')`)
        # JET.jl further widens this case to `Any` and give up further analysis rather than
        # trying hard to do sound and noisy analysis
        # xref: https://github.com/JuliaLang/julia/pull/38148
        if ret === Union{Type, TypeVar}
            return Any
        end
    end

    return ret
end

# check if this frame is for `getproperty(::Module, ::Symbol)`, which accesses to a global
# variable traced by `analyzer`
function istoplevel_globalref(analyzer::AbstractAnalyzer, sv::InferenceState)
    def = sv.linfo.def
    def.name === :getproperty || return false
    def.sig === Tuple{typeof(getproperty), Module, Symbol} || return false
    parent = sv.parent
    return !isnothing(parent) && istoplevel(parent)
end

# inter-procedural
# ================

function collect_callee_reports!(analyzer::AbstractAnalyzer, sv::InferenceState)
    reports = get_caller_cache(analyzer)
    if !isempty(reports)
        vf = get_virtual_frame(sv)
        for report in reports
            pushfirst!(report.vst, vf)
            add_new_report!(sv.result, report)
        end
        empty!(reports)
    end
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 CC.abstract_call_method(analyzer::AbstractInterpreter, method::Method, sig, sparams::SimpleVector, hardlimit::Bool, sv::InferenceState)

    collect_callee_reports!(analyzer, sv)

    @static if VERSION < v"1.8.0-DEV.510"
        # manually take in https://github.com/JuliaLang/julia/pull/42195
        if method === ISEQUAL_ANY_ANY && ret.rt === Union{Bool,Missing}
            ret = MethodCallResult(Bool, ret.edgecycle, ret.edgelimited, ret.edge)
        end
    end

    return ret
end

@static if VERSION < v"1.8.0-DEV.510"
    # manually take in https://github.com/JuliaLang/julia/pull/42195
    const ISEQUAL_ANY_ANY = let
        ms = methods(isequal)
        i = findfirst(m->m.sig===Tuple{typeof(isequal),Any,Any}, ms)::Int
        ms[i]
    end
end

function CC.abstract_call_method_with_const_args(analyzer::AbstractAnalyzer, result::MethodCallResult,
                                                 @nospecialize(f), arginfo::(@static IS_AFTER_42529 ? ArgInfo : Argtypes), match::MethodMatch,
                                                 sv::InferenceState, va_override::Bool)
    set_cacher!(analyzer, :abstract_call_method_with_const_args => sv.result)

    const_result =
        @invoke CC.abstract_call_method_with_const_args(analyzer::AbstractInterpreter, result::MethodCallResult,
                                                        @nospecialize(f), arginfo::(@static IS_AFTER_42529 ? ArgInfo : Argtypes), match::MethodMatch,
                                                        sv::InferenceState, va_override::Bool)

    # we should make sure we reset the cacher because at this point we may have not hit
    # `CC.cache_lookup(linfo::MethodInstance, given_argtypes::Argtypes, cache::JETLocalCache)`
    set_cacher!(analyzer, nothing)

    if const_result !== nothing
        # successful constant prop', we also need to update reports
        collect_callee_reports!(analyzer, sv)
    end

    return const_result
end

@static if IS_AFTER_42529
function CC.abstract_call(analyzer::AbstractAnalyzer, arginfo::ArgInfo,
                          sv::InferenceState, max_methods::Int = InferenceParams(analyzer).MAX_METHODS)
    ret = @invoke CC.abstract_call(analyzer::AbstractInterpreter, arginfo::ArgInfo,
                                   sv::InferenceState, max_methods::Int)
    analyze_task_parallel_code!(analyzer, arginfo.argtypes, sv)
    return ret
end
else # @static if IS_AFTER_42529
function CC.abstract_call(analyzer::AbstractAnalyzer, fargs::Union{Nothing,Vector{Any}}, argtypes::Argtypes,
                          sv::InferenceState, max_methods::Int = InferenceParams(analyzer).MAX_METHODS)
    ret = @invoke CC.abstract_call(analyzer::AbstractInterpreter, fargs::Union{Nothing,Vector{Any}}, argtypes::Argtypes,
                                   sv::InferenceState, max_methods::Int)
    analyze_task_parallel_code!(analyzer, argtypes, sv)
    return ret
end
end # @static if IS_AFTER_42529

"""
    analyze_task_parallel_code!(analyzer::AbstractAnalyzer, argtypes::Argtypes, sv::InferenceState)

Adds special cased analysis pass for task parallelism.
In Julia's task parallelism implementation, parallel code is represented as closure and it's
wrapped in a `Task` object. `NativeInterpreter` doesn't infer nor optimize the bodies of
those closures when compiling code that creates parallel tasks, but JET will try to run
additional analysis pass by recurring into the closures.

See also: <https://github.com/aviatesk/JET.jl/issues/114>

!!! 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, argtypes::Argtypes, sv::InferenceState)
    f = singleton_type(argtypes[1])

    # TODO we 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, here we invoke 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 nothing
end

# run additional interpretation with a new analyzer
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))
    abstract_call_method(newanalyzer, mm.method, mm.spec_types, mm.sparams, false, sv)

    return nothing
end

# `return_type_tfunc` internally uses `abstract_call` to model `Core.Compiler.return_type`
# and here we should NOT catch error reports detected within the simulated call
# because it is really not any abstraction of actual execution
function CC.return_type_tfunc(analyzer::AbstractAnalyzer, argtypes::Argtypes, sv::InferenceState)
    # stash and discard the result from the simulated call, and keep the original result (`result0`)
    result = sv.result
    result0 = result.src::JETResult
    set_result!(result)
    ret = @invoke return_type_tfunc(analyzer::AbstractInterpreter, argtypes::Argtypes, sv::InferenceState)
    set_result!(sv.result, result0)
    return ret
end

# cache
# =====

# global
# ------

"""
    JET_CACHE::$(typeof(JET_CACHE))

Keeps `src::CodeInstance` cache associated with `mi::MethodInstace` that represents the
analysis result on `mi` performed by [`analyzer::AbstractAnalyzer`](@ref AbstractAnalyzer),
where [`src.inferred::JETCachedResult`](@ref JETCachedResult) caches JET's analysis result.
This cache is separated by the identities of `AbstractAnalyzer`s, which are hash keys
computed by [`get_cache_key(analyzer::AbstractAnalyzer)`](@ref get_cache_key).

`JET_CACHE` is completely separated from the `NativeInterpreter`'s global cache, so that
JET's analysis never interacts with actual code execution.
"""
const JET_CACHE = IdDict{UInt, IdDict{MethodInstance,CodeInstance}}()

# just used for interactive developments
__clear_caches!() = empty!(JET_CACHE)

function CC.code_cache(analyzer::AbstractAnalyzer)
    cache  = JETGlobalCache(analyzer)
    worlds = WorldRange(get_world_counter(analyzer))
    return WorldView(cache, worlds)
end

struct JETGlobalCache{Analyzer<:AbstractAnalyzer}
    analyzer::Analyzer
end

# cache existence for this `analyzer` is ensured on its construction
jet_cache(analyzer::AbstractAnalyzer)       = JET_CACHE[get_cache_key(analyzer)]
jet_cache(wvc::WorldView{<:JETGlobalCache}) = jet_cache(wvc.cache.analyzer)

CC.haskey(wvc::WorldView{<:JETGlobalCache}, mi::MethodInstance) = haskey(jet_cache(wvc), mi)

function CC.typeinf_edge(analyzer::AbstractAnalyzer, method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
    # enable the report cache restoration at `code = get(code_cache(interp), mi, nothing)`
    set_cacher!(analyzer, :typeinf_edge => caller.result)
    return @invoke typeinf_edge(analyzer::AbstractInterpreter, method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
end

function CC.get(wvc::WorldView{<:JETGlobalCache}, mi::MethodInstance, default)
    codeinf = get(jet_cache(wvc), mi, default) # will ignore native code cache for a `MethodInstance` that is not analyzed by JET yet

    analyzer = wvc.cache.analyzer

    # XXX this relies on a very dirty analyzer state manipulation, the reason for this is
    # that this method (and `code_cache(::AbstractAnalyzer)`) can be called from multiple
    # contexts including edge inference, constant prop' heuristics and inlining, where we
    # want to use report cache only in edge inference, but we can't tell which context is
    # the caller of this specific method call here and thus can't tell whether we should
    # enable report cache reconstruction without the information
    # XXX move this logic into `typeinf_edge` ?
    cacher = get_cacher(analyzer)
    if isa(cacher, Pair{Symbol,InferenceResult})
        setter, caller = cacher
        if setter === :typeinf_edge
            if isa(codeinf, CodeInstance)
                # cache hit, now we need to append cached reports associated with this `MethodInstance`
                for cached in get_cached_reports(codeinf.inferred::JETCachedResult)
                    restored = add_cached_report!(caller, cached)
                    @static if JET_DEV_MODE
                        actual, expected = first(restored.vst).linfo, mi
                        @assert actual === expected "invalid global cache restoration, expected $expected but got $actual"
                    end
                    add_caller_cache!(analyzer, restored) # should be updated in `abstract_call` (after exiting `typeinf_edge`)
                end
            end
            set_cacher!(analyzer, nothing)
        end
    end

    return codeinf
end

function CC.getindex(wvc::WorldView{<:JETGlobalCache}, mi::MethodInstance)
    r = CC.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function CC.transform_result_for_cache(analyzer::AbstractAnalyzer, linfo::MethodInstance,
                                       valid_worlds::WorldRange, @nospecialize(inferred_result))
    jetresult = inferred_result::JETResult
    cache = InferenceErrorReportCache[]
    for report in get_reports(jetresult)
        @static if JET_DEV_MODE
            actual, expected = first(report.vst).linfo, linfo
            @assert actual === expected "invalid global caching detected, expected $expected but got $actual"
        end
        cache_report!(cache, report)
    end
    return JETCachedResult(cache, get_source(jetresult))
end

function CC.setindex!(wvc::WorldView{<:JETGlobalCache}, ci::CodeInstance, mi::MethodInstance)
    setindex!(jet_cache(wvc), ci, mi)
    return nothing
end

function add_jet_callback!(linfo)
    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_jet_cache!]
    else
        callbacks = linfo.callbacks::Vector{Any}
        if !any(function (@nospecialize(cb),)
                    cb === invalidate_jet_cache!
                end,
                callbacks)
            push!(callbacks, invalidate_jet_cache!)
        end
    end
    return nothing
end

function invalidate_jet_cache!(replaced, max_world, depth = 0)
    for cache in values(JET_CACHE)
        delete!(cache, replaced)
    end

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !any(cache->haskey(cache, mi), values(JET_CACHE))
                continue # otherwise fall into infinite loop
            end
            invalidate_jet_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

# local
# -----

struct JETLocalCache{Analyzer<:AbstractAnalyzer}
    analyzer::Analyzer
    cache::Vector{InferenceResult}
end

CC.get_inference_cache(analyzer::AbstractAnalyzer) = JETLocalCache(analyzer, get_inference_cache(get_native(analyzer)))

function CC.cache_lookup(linfo::MethodInstance, given_argtypes::Argtypes, cache::JETLocalCache)
    # XXX the very dirty analyzer state observation again
    # this method should only be called from the single context i.e. `abstract_call_method_with_const_args`,
    # and so we should reset the cacher immediately we reach here
    analyzer = cache.analyzer
    setter, caller = get_cacher(analyzer)::Pair{Symbol,InferenceResult}
    @assert setter === :abstract_call_method_with_const_args
    set_cacher!(analyzer, nothing)

    inf_result = cache_lookup(linfo, given_argtypes, cache.cache)

    isa(inf_result, InferenceResult) || return inf_result

    # constant prop' hits a cycle (recur into same non-constant analysis), we just bail out
    isa(inf_result.result, InferenceState) && return inf_result

    # cache hit, try to restore local report caches

    # corresponds to the throw-away logic in `_typeinf(analyzer::AbstractAnalyzer, frame::InferenceState)`
    filter!(!is_from_same_frame(caller.linfo, linfo), get_reports(caller))

    for cached in get_cached_reports(inf_result)
        restored = add_cached_report!(caller, cached)
        @static if JET_DEV_MODE
            actual, expected = first(restored.vst).linfo, linfo
            @assert actual === expected "invalid local cache restoration, expected $expected but got $actual"
        end
        add_caller_cache!(analyzer, restored) # should be updated in `abstract_call_method_with_const_args`
    end

    return inf_result
end

CC.push!(cache::JETLocalCache, inf_result::InferenceResult) = CC.push!(cache.cache, inf_result)

# main driver
# ===========

# in this overload we will work on some meta/debug information management per inference frame
function CC.typeinf(analyzer::AbstractAnalyzer, frame::InferenceState)
    (; linfo, parent, result) = frame
    isentry = isnothing(parent)

    #= logging stage1 start =#
    local sec::Float64
    logger_activated = isa(JETLogger(analyzer).inference_logger, IO)
    depth = get_depth(analyzer)
    if logger_activated
        sec = time()
        with_inference_logger(analyzer, ==(DEBUG_LOGGER_LEVEL)) do @nospecialize(io)
            print_rails(io, depth)
            printstyled(io, "┌ @ "; color = RAIL_COLORS[(depth+1)%N_RAILS+1])
            print(io, linfo)
            file, line = get_file_line(linfo)
            print(io, ' ', file, ':', line)
            println(io)
            set_depth!(analyzer, get_depth(analyzer) + 1) # manipulate this only in debug mode
        end
    end
    #= logging stage1 end =#

    # some methods like `getproperty` can't propagate accurate types without actual values,
    # and constant prop' plays a somewhat critical role in those cases by overwriteing the
    # previous non-constant inference result (under the current design constant prop' always
    # happens after inference with non-constant abstract elements)
    # JET also needs that in order to reduce false positive reports, and here we will
    # throw-away previously-collected error reports that are "lineage" of this frame,
    # when it is being re-inferred with constants
    # NOTE `frame.linfo` is the exactly same object as that of the previous non-constant inference
    # IDEA we may still want to keep some "serious" error reports like `GlobalUndefVarErrorReport`
    # even when constant prop' reveals it never happ∫ens given the current constant arguments
    if is_constant_propagated(frame) && !isentry
        filter!(!is_from_same_frame(parent.linfo, linfo), get_reports(parent.result))
    end

    @assert isa(result.src, JETResult)

    ret = @invoke typeinf(analyzer::AbstractInterpreter, frame::InferenceState)

    ret && @assert isa(result.src, isentry ? JETResult : JETCachedResult)

    #= logging stage2 start =#
    if logger_activated
        elapsed = round(time() - sec; digits = 3)
        with_inference_logger(analyzer, ==(INFO_LOGGER_LEVEL)) do @nospecialize(io)
            println(io, "inference on $linfo finished in $elapsed sec")
        end
        with_inference_logger(analyzer, ==(DEBUG_LOGGER_LEVEL)) do @nospecialize(io)
            print_rails(io, depth)
            printstyled(io, "└─→ "; color = RAIL_COLORS[(depth+1)%N_RAILS+1])
            printstyled(io, frame.bestguess; color = TYPE_ANNOTATION_COLOR)
            println(io, " (", join(filter(!isnothing, (
                             linfo,
                             ret ? nothing : "in cycle",
                             string(length((isentry ? get_reports : get_cached_reports)(result)), " reports"),
                             string(elapsed, " sec"),
                             )), ", "),
                         ')')
            set_depth!(analyzer, get_depth(analyzer) - 1) # manipulate this only in debug mode
        end
    end
    #= logging stage2 end =#

    return ret
end

"""
    is_from_same_frame(parent_linfo::MethodInstance, current_linfo::MethodInstance) ->
        (report::InferenceErrorReport) -> Bool

Returns a function that checks if a given `InferenceErrorReport` is generated from `current_linfo`.
It also checks `current_linfo` is a "lineage" of `parent_linfo` (i.e. entered from it).

This function is supposed to be used to filter out reports collected from analysis on `current_linfo`
without using constants when entering into the constant analysis. As such, this function
assumes that when a report should be filtered out, the first elment of its virtual stack
frame `st` is for `parent_linfo` and the second element of that is for `current_linfo`.

Example: Assume `linfo2` will produce a report for some reason.
```
entry
└─ linfo1
   ├─ linfo2 (report1: linfo2)
   ├─ linfo3 (report1: linfo1->linfo2, report2: linfo3->linfo2)
   │  └─ linfo2 (report1: linfo1->linfo2, report2: linfo2)
   └─ linfo3′ (report1: linfo1->linfo2, ~~report2: linfo1->linfo3->linfo2~~)
```
In the example analysis above, `report2` will be filtered out on re-entering into `linfo3′`
(i.e. we're analyzing `linfo3` with constants argument), because
`is_from_same_frame(linfo1, linfo3)(report2)` returns `true`.
Note that `report1` is still kept there because of the lineage check, i.e.
`is_from_same_frame(linfo1, linfo3)(report1)` returns `false`.
"""
function is_from_same_frame(parent_linfo::MethodInstance,
                            current_linfo::MethodInstance,
                            )
    function (report::InferenceErrorReport)
        @inbounds begin
            vst = report.vst
            length(vst) > 1 || return false
            vst[1].linfo === parent_linfo || return false
            return vst[2].linfo === current_linfo
        end
    end
end

# in this overload we can work on `frame.src::CodeInfo` (and also `frame::InferenceState`)
# where type inference (and also optimization if applied) already ran on
function CC._typeinf(analyzer::AbstractAnalyzer, frame::InferenceState)
    CC.typeinf_nocycle(analyzer, frame) || return false # frame is now part of a higher cycle
    # with no active ip's, frame is done
    frames = frame.callers_in_cycle
    isempty(frames) && push!(frames, frame)
    valid_worlds = WorldRange()
    for caller in frames
        @assert !(caller.dont_work_on_me)
        caller.dont_work_on_me = true
        # might might not fully intersect these earlier, so do that now
        valid_worlds = CC.intersect(caller.valid_worlds, valid_worlds)
    end
    for caller in frames
        caller.valid_worlds = valid_worlds
        CC.finish(caller, analyzer)
        # finalize and record the linfo result
        caller.inferred = true
    end
    # NOTE we don't discard `InferenceState`s here so that some analyzers can use them in `finish!`
    # # collect results for the new expanded frame
    # results = Tuple{InferenceResult, Vector{Any}, Bool}[
    #         ( frames[i].result,
    #           frames[i].stmt_edges[1]::Vector{Any},
    #           frames[i].cached )
    #     for i in 1:length(frames) ]
    # empty!(frames)
    for frame in frames
        caller = frame.result
        opt = get_source(caller)
        if opt isa OptimizationState # implies `may_optimize(analyzer) === true`
            result_type = caller.result
            @assert !(result_type isa LimitedAccuracy)
            CC.optimize(analyzer, opt, OptimizationParams(analyzer), result_type)
            # # COMBAK we may want to enable inlining ?
            # if opt.const_api
            #     # XXX: The work in ir_to_codeinf! is essentially wasted. The only reason
            #     # we're doing it is so that code_llvm can return the code
            #     # for the `return ...::Const` (which never runs anyway). We should do this
            #     # as a post processing step instead.
            #     CC.ir_to_codeinf!(opt)
            #     if result_type isa Const
            #         caller.src = result_type
            #     else
            #         @assert CC.isconstType(result_type)
            #         caller.src = Const(result_type.parameters[1])
            #     end
            # end
            caller.valid_worlds = CC.getindex((opt.inlining.et::CC.EdgeTracker).valid_worlds)
        end
    end

    for frame in frames
        caller = frame.result
        edges = frame.stmt_edges[1]::Vector{Any}
        cached = frame.cached
        valid_worlds = caller.valid_worlds
        if CC.last(valid_worlds) >= get_world_counter()
            # if we aren't cached, we don't need this edge
            # but our caller might, so let's just make it anyways
            CC.store_backedges(caller, edges)
        end
        CC.finish!(analyzer, frame)

        # XXX this is a dirty fix for performance problem, we need more "proper" fix
        # https://github.com/aviatesk/JET.jl/issues/75
        unique!(aggregation_policy(analyzer), get_reports(caller))

        # global cache management
        # part 1: transform collected reports to `JETCachedResult` and put it into `CodeInstance.inferred`
        if cached && !istoplevel(frame)
            CC.cache_result!(analyzer, caller)
        end
        # part 2: register invalidation callback for JET cache
        add_jet_callback!(caller.linfo)

        reports = get_reports(caller)
        if frame.parent !== nothing
            # inter-procedural handling: get back to the caller what we got from these results
            add_caller_cache!(analyzer, reports)

            # local cache management
            # TODO there are duplicated work here and `transform_result_for_cache`
            cache = InferenceErrorReportCache[]
            for report in reports
                cache_report!(cache, report)
            end
            set_cached_result!(caller, cache)
        end
    end

    return true
end

function CC.finish(me::InferenceState, analyzer::AbstractAnalyzer)
    # @invoke CC.finish(me::InferenceState, analyzer::AbstractInterpreter)
    # prepare to run optimization passes on fulltree
    s_edges = me.stmt_edges[1]
    if s_edges === nothing
        s_edges = me.stmt_edges[1] = []
    end
    for edges in me.stmt_edges
        edges === nothing && continue
        edges === s_edges && continue
        append!(s_edges, edges)
        empty!(edges)
    end
    if me.src.edges !== nothing
        edges = me.src.edges
        if isa(edges, Vector{Any})
            append!(s_edges, edges)
        else
            # this pass should never happen in ordinal code, but some external
            # code generator like IRTools.jl may produce this case
            # see https://github.com/aviatesk/JET.jl/issues/271
            append!(s_edges, edges::Vector{MethodInstance})
        end
        me.src.edges = nothing
    end
    # inspect whether our inference had a limited result accuracy,
    # else it may be suitable to cache
    me.bestguess = CC.cycle_fix_limited(me.bestguess, me)
    limited_ret = me.bestguess isa LimitedAccuracy
    limited_src = false
    if !limited_ret
        gt = me.src.ssavaluetypes::Vector{Any}
        for j = 1:length(gt)
            gt[j] = gtj = CC.cycle_fix_limited(gt[j], me)
            if gtj isa LimitedAccuracy && me.parent !== nothing
                limited_src = true
                break
            end
        end
    end
    if limited_ret
        # a parent may be cached still, but not this intermediate work:
        # we can throw everything else away now
        set_source!(me.result, nothing)
        me.cached = false
        me.src.inlineable = false
        unlock_mi_inference(analyzer, me.linfo)
    elseif limited_src
        # a type result will be cached still, but not this intermediate work:
        # we can throw everything else away now
        set_source!(me.result, nothing)
        me.src.inlineable = false
    else
        # annotate fulltree with type information,
        # either because we are the outermost code, or we might use this later
        doopt = (me.cached || me.parent !== nothing)
        CC.type_annotate!(me, doopt)
        if doopt && may_optimize(analyzer)
            set_source!(me.result, OptimizationState(me, OptimizationParams(analyzer), analyzer))
        else
            set_source!(me.result, me.src::CodeInfo) # stash a convenience copy of the code (e.g. for reflection)
        end
    end
    me.result.valid_worlds = me.valid_worlds
    me.result.result = me.bestguess

    if istoplevel(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′) # ignored by default
                    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) # ignored by default
        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′
    @nospecialize t
end
get_msg(::Type{InvalidConstantRedefinition}, sv::InferenceState, mod::Module, name::Symbol, @nospecialize(t′::Any), @nospecialize(t::Any)) =
    "invalid redefinition of constant $(mod).$(name) (from $(t′) to $(t))"

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

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_abstract_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

function CC.finish!(analyzer::AbstractAnalyzer, frame::InferenceState)
    caller = frame.result

    # transform optimized `IRCode` to optimized `CodeInfo`
    src = get_source(caller)
    if src isa OptimizationState # implies `may_optimize(analyzer) === true`
        opt = src
        @assert opt.ir !== nothing # `_typeinf(::AbstractAnalyzer, ::InferenceState)` disabled `const_api`

        src = CC.ir_to_codeinf!(opt)
        set_source!(caller, src)
    end

    return src
end


text:
using LinearAlgebra, DataFrames,Optim, ForwardDiff, BenchmarkTools,Distributions,
    Expectations, QuantEcon, Statistics, GLM
using Distributions: invsqrt2π, log2π, sqrt2, invsqrt2
using DCDC
using Test
using Distributed, Suppressor
using Plots

begin "Dynamic Decision Process"
    σ₀ = 2;
    β = 0.8;
    nM = 50;
    nT = 5;
    ddc=DynamicDecisionProcess(2,0.8);
    computeEquilibrium(ddc);
    checked = [check_ee(ddc) for i = 1:100];
    plot(ddc.ValueFn.xdata,ddc.ValueFn.y);
    plot!(ddc.PolicyFn.xdata,ddc.PolicyFn.y);
    data = simulate_ddc(nM,nT,ddc);
end

# The moment condition is
# u'(a_t) = β d_Trans(a_t,x_t) u'(a_t+1)

reshaped_data = convert_data(data);

using Flux #Get derivatives
using Flux.Tracker
using Flux.Tracker: update!

θ₀ = ones(2)  #θ = (β,σ)
Transition = ddc.trans;
dtrans = (x,c) -> Tracker.gradient(Transition,x,c);
dtransx = (x,c) -> Tracker.gradient(Transition,x,c)[1];
R = dtransx.(reshaped_data.action_t,reshaped_data.state_t);
θ = θ₀

lb = [0,0];
ub = [1,100];
function _moment(θ)
    (β₀,σ₀ ) = θ₀;
    util = Utility(σ₀);
    dutil = (x) -> Tracker.gradient(util,x)[1].data;
    dudc  = dutil.(reshaped_data.action_t);
    dudc1 = dutil.(reshaped_data.action_t_1);
    η = dudc - β₀ .* dudc1 .* R;
    return(η'*η)
end

function _moment2(σ₀)
    dutil = (x) -> Tracker.gradient(Utility(σ₀);,x)[1].data;
    dudc  = dutil.(reshaped_data.action_t);
    dudc1 = dutil.(reshaped_data.action_t_1);
    η = dudc - 0.8 .* dudc1 .* R;
    return(η'*η)
end

optimize(_moment,zeros(2),BFGS(),autodiff=:forward)

ext: jl


text:
export Prostate

"""
    Prostate

The dataset is from the following reference:

> Dinesh Singh, Phillip G. Febbo, Kenneth Ross, Donald G. Jackson, Judith Manola,
 Christine Ladd, Pablo Tamayo, Andrew A. Renshaw, Anthony V. D’Amico, Jerome P. Richie,
 Eric S. Lander, Massimo Loda, Philip W. Kantoff, Todd R. Golub, and William R. Sellers.
 Gene expression correlates of clinical prostate cancer behavior.
 Cancer cell, 1(2): 203–209, 2002.

See the following monograph for further illustrations
of empirical Bayes methods on this dataset:

> Bradley Efron. Large-scale inference: Empirical Bayes methods for estimation, testing,
and prediction. Cambridge University Press, 2012.
"""
module Prostate

using CSV
using ..Empirikos: MultinomialSummary, StandardNormalSample

const URL = "https://web.stanford.edu/~hastie/CASI_files/DATA/prostz.txt"

function load_table()
    CSV.File(download(URL), header=false)
end

function ebayes_samples()
    tbl = load_table()
    Zs = StandardNormalSample.(tbl.Column1)
end

end


text:
using StaticArrays

module MyPkgTests

using LinearAlgebra
using Test
using StaticArrays
using JLD2
import CollisionDetection

include("test_core.jl")
include("test_bloated.jl")
include("test_searcheq.jl")

end # module


text:
using ROOTFramework

bindings = TTreeBindings()
idx = bindings[:idx] = Ref(zero(Int32))
s = bindings[:s] = Ref(zero(Float64))
a = bindings[:a] = zeros(Float64, 2)

open(TFile, "out.root", "recreate") do tfile
    ttree = create_ttree!(tfile, "data", "Data")
    output = TTreeOutput(ttree, bindings)

    @time for entry in 1:1000000
        idx[] = entry
        s[] = 2.0 + 2.0 * randn()

        a[1] = 12.0 + 2.0 * randn()
        a[2] = 22.0 + 2.0 * randn()

        push!(output)
    end
end


text:
import NbodyGradient: set_state!, zero_out!, amatrix

@testset "Transit Parameters" begin
    N = 3
    t0 = 7257.93115525 - 7300.0 - 0.5 # Initialize IC before first transit
    h = 0.04
    itime = 10.0 # Time integrator will run
    tmax = itime

    # Setup initial conditions:
    elements = readdlm("elements.txt", ',')[1:N,:]
    elements[2:end,3] .-= 7300.0
    elements[:,end] .= 0.0 # Set Ωs to 0
    elements[2,1] *= 10.0
    elements[3,1] *= 10.0
    ic = ElementsIC(t0, N, elements)

    function calc_times(h)
        intr = Integrator(ahl21!, h, 0.0, tmax)
        s = State(ic)
        ttbv = TransitParameters(itime, ic)
        intr(s, ttbv)
        return ttbv
    end

    # Calculate transit times
    tts = calc_times(h)

    mask = zeros(Bool, size(tts.dtbvdq0));
    for jq = 1:N
        for iq = 1:7
            if iq == 7; ivary = 1; else; ivary = iq + 1; end  # Shift mass variation to end
            for i = 2:N
                for k = 1:tts.count[i]
                    for itbv = 1:3
                        # Ignore inclination & longitude of nodes variations:
                        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
                            mask[itbv,i,k,iq,jq] = true
                        end
                    end
                end
            end
        end
    end

    function calc_finite_diff(h, t0, itime, tmax, elements)
    # Now do finite difference with big float
        dq0 = big(1e-10)
        ic_big = ElementsIC(big(t0), N, big.(elements))
        elements_big = copy(ic_big.elements)
        s_big = State(ic_big)
        ttp = TransitParameters(big(itime), ic_big);
        ttm = TransitParameters(big(itime), ic_big);
        dtde_num = zeros(BigFloat, size(ttp.dtbvdelements));
        intr_big = Integrator(big(h), zero(BigFloat), big(tmax))

        for jq in 1:N
            for iq in 1:7
                zero_out!(ttp); zero_out!(ttm);
                ic_big.elements .= elements_big
                ic_big.m .= elements_big[:,1]
                if iq == 7; ivary = 1; else; ivary = iq + 1; end
                ic_big.elements[jq,ivary] += dq0
                if ivary == 1; ic_big.m[jq] += dq0; amatrix(ic_big); end # Masses don't update with elements array
                sp = State(ic_big);
                intr_big(sp, ttp; grad=false)
                ic_big.elements[jq,ivary] -= 2dq0
                if ivary == 1; ic_big.m[jq] -= 2dq0; amatrix(ic_big); end
                sm = State(ic_big)
                intr_big(sm, ttm; grad = false)
                for i in 2:N
                    for k in 1:ttm.count[i]
                        for itbv=1:3
                            # Compute double-sided derivative for more accuracy:
                            dtde_num[itbv,i,k,iq,jq] = (ttp.ttbv[itbv,i,k] - ttm.ttbv[itbv,i,k]) / (2dq0)
                        end
                    end
                end
            end
        end
        return dtde_num
    end

    dtde_num = calc_finite_diff(h, t0, itime, tmax, elements)
    @test isapprox(asinh.(tts.dtbvdelements[mask]), asinh.(dtde_num[mask]);norm=maxabs)
    @test isapprox(asinh.(tts.dtbvdelements), asinh.(dtde_num);norm=maxabs)
end


text:
# This file is part of IntegerTriangles.jl.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module TrianglesExamples

using Nemo, TrianglesBase

export LahNumbers, LahTriangle, LahTransform, SchröderBTriangle, SchröderLTriangle, MotzkinTriangle
export Catalan, CatalanTriangle, CatalanTransform, CatalanBallot, ExtCatalanTriangle
export BernoulliPolynomial, PascalTriangle, SchroederBigTriangle, AitkenTriangle
export EulerianTriangle, EulerianTriangle2, NarayanaTriangle, NarayanaTransform
export EulerTriangle, EulerTanTriangle, EulerSecTriangle, UniTriangle
export EulerianTransform, MotzkinTransform, SchroederBigTransform, FubiniTriangle
export Laguerre, LaguerreTriangle, LaguerreTransform, TrinomialTriangle, TrinomialTransform
export JacobsthalTriangle, JacobsthalTransform, FibonacciTriangle, FibonacciTransform
export StirlingSetTriangle, StirlingCycleTriangle, FallingFactTriangle, RisingFactTriangle
export StirlingSetTransform, StirlingCycleTransform, FubiniTriangle, RencontresTriangle
export DArcaisTriangle, WorpitzkyTriangle, EulerianSO2Triangle, SwingTriangle, DelannoyTriangle
export Bessel1Triangle, Bessel1Transform, HermiteTriangle, HermiteTransform
export I132393, I048993, I271703, I094587, I008279, I225478, T132393, T048993
export T094587, T008279, T225478, T271703, FineTriangle, TTreeTriangle
export A046802Triangle, A046802Transform, A343237Triangle
export A000166, A038048, V000111, Fine
export TRIANGLES

#=
function __init__()
    global CacheA000166 = Dict{Int,ℤInt}(0 => ZZ(1))
    global CacheA038048 = Dict{Int,ℤInt}(0 => ZZ(0))
    global CacheFine = Dict{Int,ℤInt}(0 => ZZ(1), 1 => ZZ(1), 2 => ZZ(0))
    global CacheAndré = Dict{Tuple{Int,Int},fmpz}()

    global CacheBeta = Dict{Tuple{Int,Int},fmpz}()
    global CacheLah = Dict{Int,ℤSeq}([0 => [ZZ(1)]])
    global CacheFubini = Dict{Int,ℤSeq}([0 => [ZZ(1)]])
    global CacheDelannoy = Dict{Tuple{Int,Int},ℤInt}([(0,0) => ZZ(1), (1,0) => ZZ(1), (1,1) => ZZ(1)])
    global CacheAitken = Dict{Int,ℤSeq}([0 => [ZZ(1)]])
    global CacheBallot = Dict{Tuple{Int,Int},ℤInt}()
    global CacheDArcais = Dict{Tuple{Int,Int},ℤInt}()
    global CacheWorpitzky = Dict{Tuple{Int,Int},ℤInt}()
    global CacheEulerianSO2 = Dict{Tuple{Int,Int},ℤInt}((0,0) => ZZ(1))
    global CacheSwing = Dict{Tuple{Int,Int},ℤInt}((0,0) => ZZ(1))
    global CacheEulerian = Dict{Tuple{Int,Int},ℤInt}()
    global CacheEulerianClassic = Dict{Tuple{Int,Int},fmpz}()
    global CacheNarayana = Dict{Tuple{Int,Int},fmpz}()
    global CacheLaguerre = Dict{Int,ℤSeq}(0 => [ZZ(1)])
end
=#

const CacheA000166 = Dict{Int,ℤInt}(0 => ZZ(1))
const CacheA038048 = Dict{Int,ℤInt}(0 => ZZ(0))
const CacheFine = Dict{Int,ℤInt}(0 => ZZ(1), 1 => ZZ(1), 2 => ZZ(0))
const CacheAndré = Dict{Tuple{Int,Int},fmpz}()

const CacheBeta = Dict{Tuple{Int,Int},fmpz}()
const CacheLah = Dict{Int,ℤSeq}([0 => [ZZ(1)]])
const CacheFubini = Dict{Int,ℤSeq}([0 => [ZZ(1)]])
const CacheDelannoy = Dict{Tuple{Int,Int},ℤInt}([(0,0) => ZZ(1), (1,0) => ZZ(1), (1,1) => ZZ(1)])
const CacheAitken = Dict{Int,ℤSeq}([0 => [ZZ(1)]])
const CacheBallot = Dict{Tuple{Int,Int},ℤInt}()
const CacheDArcais = Dict{Tuple{Int,Int},ℤInt}()
const CacheWorpitzky = Dict{Tuple{Int,Int},ℤInt}()
const CacheEulerianSO2 = Dict{Tuple{Int,Int},ℤInt}((0,0) => ZZ(1))
const CacheSwing = Dict{Tuple{Int,Int},ℤInt}((0,0) => ZZ(1))
const CacheEulerian = Dict{Tuple{Int,Int},ℤInt}()
const CacheEulerianClassic = Dict{Tuple{Int,Int},fmpz}()
const CacheNarayana = Dict{Tuple{Int,Int},fmpz}()
const CacheLaguerre = Dict{Int,ℤSeq}(0 => [ZZ(1)])

"""
Recurrences and iterations for some triangles.
The examples given are:
```
Aitken, Bessel1, Binomial, Catalan, DArcais,

Delannoy, Eulerian, EulerianSO2, EulerSec,

EulerTan, Euler, FallingFact, Fibonacci,

Fine, Fubini, Hermite, Laguerre, Lah, Motzkin,

Narayana, Rencontres, RisingFact, SchröderB,

SchröderL, StirlingCycle, StirlingSet, Trinomial,

TTree, Uni, Worpitzky
```
"""
const ModuleTrianglesExamples = ""

# ------------------------------------------------

function A000166(n::Int)
    n <= 1 && return ZZ(1 - n)
    haskey(CacheA000166, n) && return CacheA000166[n]
    a = (n - 1) * (A000166(n - 1) + A000166(n - 2))
    return CacheA000166[n] = a
end

function A038048(n::Int)
    haskey(CacheA038048, n) && return CacheA038048[n]
	s = Factorial(n - 1) * divisor_sigma(n, 1)
    return CacheA038048[n] = s
end

function Fine(n::Int)
    haskey(CacheFine, n) && return CacheFine[n]
    s = div((7*n - 12)*Fine(n-1) + (4*n - 6)*Fine(n-2), 2*n)
    CacheFine[n] = s
end

function André(m::Int, n::Int)
    haskey(CacheAndré, (m, n)) && return CacheAndré[(m, n)]
    n ≤ 0 && return ZZ(1)
    r = range(0, step = m, stop = n - 1)
    S = sum(binomial(ZZ(n), ZZ(k)) * André(m, k) for k ∈ r)
    return CacheAndré[(m, n)] = n % m == 0 ? -S : S
end

"""
Return the up-down numbers (2-alternating permutations).
"""
V000111(n::Int) = abs(André(2, n))

function PrimeDivisors(n)
    n < 2 && return ℤInt[]
    sort!([p for (p, e) ∈ factor(ZZ(n))])
end

# ------------------------------------------------

Euler(n, k) = Binomial(n, k) * V000111(n - k)
EulerTriangle(dim) = [[Euler(n, k) for k = 0:n] for n = 0:dim - 1]
Euler(A::ℤSeq) = LinMap(Euler, A, length(A))
EulerTransform(A::ℤSeq) = Euler.(Telescope(A))

# ------------------------------------------------

"""
https://oeis.org/wiki/User:Peter_Luschny/SwissKnifePolynomials
"""
function BetaPoly(n)

    function v(n, k)
        if haskey(CacheBeta, (n, k))
            CacheBeta[(n, k)]
        else
            CacheBeta[(n, k)] = binomial(n, k)*subst(BetaPoly(k), 0)
        end
    end

    R, z = ZPolyRing("z")
    n == 0 && return R(1)
    R(sum(isodd(k) ? 0 : v(n, k)*(z - 1)^(n - k - 1) for k in 0:n-1))
end

SECH(n) = isodd(n)  ? 0 : subst(BetaPoly(n), 0)
TANH(n) = iseven(n) ? 0 : subst(BetaPoly(n), 0)

function SecPoly(n)
    R, z = ZPolyRing("z")
    R(sum(Binomial(n, k)*SECH(k)*z^(n-k) for k in 0:n))
end

function TanPoly(n)
    R, z = ZPolyRing("z")
    n == 0 && return R(0)
    R(sum(Binomial(n, k)*TANH(n - k)*z^k for k in 0:n-1))
end

function SwissKnifePoly(n, p, q)
    R, z = ZPolyRing("z")
    P = p == 0 ? 0 : SecPoly(n)
    Q = q == 0 ? 0 : TanPoly(n)
    p*P + q*Q
end

EulerTanTriangle(dim) = [Coefficients(TanPoly(n), n) for n in 0:dim-1]
EulerSecTriangle(dim) = Coefficients([SecPoly(n) for n in 0:dim-1])

# ------------------------------------------------

FineTriangle(dim::Int) = RiordanSquare(dim, Fine)
TernaryTree(n::Int) = div(Binomial(ZZ(3*n), n), ZZ(2*n + 1))
TTreeTriangle(dim::Int) =  RiordanSquare(dim, TernaryTree)

# ------------------------------------------------

"""
(SIGNATURES)
"""
function LahNumbers(n::Int64)
    haskey(CacheLah, n) && return CacheLah[n]
    prevrow = LahNumbers(n - 1)
    row = ZSeq(n + 1)
    row[1] = 0
    row[n + 1] = 1
    for k = 2:n
        row[k] = prevrow[k - 1] + prevrow[k] * (n + k - 2)
    end
    CacheLah[n] = row
end

LahNumbers(n, k) = LahNumbers(n)[k + 1]

"""
(SIGNATURES)
"""
function LahTriangle(size)
    length(CacheLah) < size && LahNumbers(size)
    [CacheLah[n] for n = 0:size - 1]
end

LahNumbers(A::ℤSeq)   = LinMap(LahNumbers, A, length(A))
LahTransform(A::ℤSeq) = LahNumbers.(Telescope(A))

# ------------------------------------------------

"""
(SIGNATURES)
"""
function FubiniNumbers(n::Int64)
    haskey(CacheFubini, n) && return CacheFubini[n]
    prevrow = FubiniNumbers(n - 1)
    row = ZSeq(n + 1)
    row[1] = 0
    row[2] = 1
    row[n + 1] = n * prevrow[n]
    for k = 3:n
        row[k] = (k - 1)*(prevrow[k - 1] + prevrow[k])
    end
    CacheFubini[n] = row
end

FubiniNumbers(n, k) = FubiniNumbers(n)[k + 1]

"""
(SIGNATURES)
"""
function FubiniTriangle(size)
    length(CacheFubini) < size && FubiniNumbers(size)
    [CacheFubini[n] for n = 0:size - 1]
end

FubiniNumbers(A::ℤSeq) = LinMap(FubiniNumbers, A, length(A))
FubiniTransform(A::ℤSeq) = FubiniNumbers.(Telescope(A))

# ------------------------------------------------

"""
(SIGNATURES)
"""
function Delannoy(n::Int64, k::Int64)
    haskey(CacheDelannoy, (n, k)) && return CacheDelannoy[(n, k)]
    if n < 0 || k < 0 || k > n
        return CacheDelannoy[(n,k)] = 0
    end
    CacheDelannoy[(n,k)] = Delannoy(n-1, k-1) + Delannoy(n-2, k-1) + Delannoy(n-1, k)
end

Delannoy(n) = [Delannoy(n, k) for k ∈ 0:n]
DelannoyTriangle(dim) = [Delannoy(n) for n = 0:dim - 1]
Delannoy(A::ℤSeq) = LinMap(Delannoy, A, length(A))
DelannoyTransform(A::ℤSeq) = Delannoy.(Telescope(A))

# ------------------------------------------------

"""
(SIGNATURES)
"""
function AitkenNumbers(n::Int64)
    haskey(CacheAitken, n) && return CacheAitken[n]
    prevrow = AitkenNumbers(n - 1)
    row = ZSeq(n + 1)
    row[1] = prevrow[n]
    for k = 2:n + 1
        row[k] = row[k-1] + prevrow[k - 1]
    end
    CacheAitken[n] = row
end

AitkenNumbers(n, k) = AitkenNumbers(n)[k + 1]

"""
(SIGNATURES)
"""
function AitkenTriangle(size)
    length(CacheAitken) < size && AitkenNumbers(size)
    [CacheAitken[n] for n = 0:size - 1]
end

AitkenNumbers(A::ℤSeq) = LinMap(AitkenNumbers, A, length(A))
AitkenTransform(A::ℤSeq) = AitkenNumbers.(Telescope(A))

# ------------------------------------------------
#row[k] = ( get(prevrow, k - 1, 0)
#+ get(prevrow, k, 0) * (2 * k - 1)
#+ get(prevrow, k + 1, 0) * k^2 )

#A011971 := proc(n, k) option remember; if n=0 and k=0 then 1 elif k=0
#then A011971(n-1, n-1) else A011971(n, k-1) + A011971(n-1, k-1); fi: end;

G271703(x, t) = exp(t * divexact(x, 1 - t))
T271703(dim) = EgfExpansionCoeff(dim, G271703)

egfBernoulli(x, t) = divexact(t * exp(x * t), exp(t) - 1)
BernoulliPolynomial(n) = EgfExpansionPoly(n, egfBernoulli)

G278075(x, t) = divexact(1, 1 - x * (1 - exp(-t)))
FubiniPolynomial(n) = EgfExpansionPoly(n, G278075)

function MotzkinTriangle2(dim::Int)
    T = ZTri(dim)
    for n = 1:dim
        T[n][n] = 1
    end
    for n = 2:dim
        u(k) = k == 0 || k == n ? 0 : T[n - 1][k]
        for k = 1:n - 1
            T[n][k] = u(k - 1) + u(k) + u(k + 1)
        end
    end
    T
end

# ------------------------------------------------

"""
(SIGNATURES)
"""
MotzkinTriangle(dim) = OrthoPoly(dim, n -> 1, n -> 1)

"""
(SIGNATURES)
"""
Motzkin(n) = MotzkinTriangle(n + 1)[n + 1] # TODO!!

Motzkin(n, k) = Motzkin(n)[k + 1]
Motzkin(A::ℤSeq) = LinMap(Motzkin, A, length(A))
MotzkinTransform(A::ℤSeq) = Motzkin.(Telescope(A))

# ------------------------------------------------

"""
(SIGNATURES)
"""
CatalanTriangle(dim) = OrthoPoly(dim, n -> 0, n -> 1)
Catalan(n) = CatalanTriangle(n + 1)[n + 1] # TODO!!
Catalan(n::Int, k::Int) = Catalan(n)[k + 1]
Catalan(A::ℤSeq) = LinMap(Catalan, A, length(A))
CatalanTransform(A::ℤSeq) = Catalan.(Telescope(A))

function CatalanBallot(n::Int, k::Int)
    haskey(CacheBallot, (n, k)) && return CacheBallot[(n, k)]
    (k > n || k < 0) && return ZZ(0)
    n == k && return ZZ(1)
    CacheBallot[(n, k)] = (
            CatalanBallot(n - 1, k - 1) +
            CatalanBallot(n - 1, k + 1) +
            (iseven(n - k) ? 0 : CatalanBallot(n - 1, k))
        )
end

CatalanBallot(n) = [CatalanBallot(n, k) for k ∈ 0:n]
ExtCatalanTriangle(dim) = [[CatalanBallot(n, k) for k = 0:n] for n = 0:dim - 1]

# ------------------------------------------------

# CatalanConvolution A106566
# T(n, k) = binomial(2n-k-1, n-k)*k/n for 0 <= k <= n
# with n > 0; T(0, 0) = 1; T(0, k) = 0 if k > 0.
# The following are all versions of (essentially) the same
# Catalan triangle: A009766, A030237, A033184, A059365,
# A099039, A106566, A130020, A047072.

# ------------------------------------------------

Rencontres(n, k) = A000166(n - k)*Binomial(n, k)
RencontresTriangle(dim) = [[Rencontres(n, k) for k = 0:n] for n = 0:dim - 1]

#A008290(n) = A008290Triangle(n + 1)[n + 1]
#A008290(n, k) = A008290(n)[k + 1]
Rencontres(A::ℤSeq) = LinMap(Rencontres, A, length(A))
RencontresTransform(A::ℤSeq) = Rencontres.(Telescope(A))

# ------------------------------------------------

function DArcais(n, k)
    haskey(CacheDArcais, (n, k)) && return CacheDArcais[(n, k)]
    k == 0 && return ZZ(0^n)
    S = sum(Binomial(n-1, j-1)*A038048(j)*DArcais(n-j, k-1) for j in 0:n-k+1)
    CacheDArcais[(n, k)] = S
end

DArcaisTriangle(dim) = [[DArcais(n, k) for k = 0:n] for n = 0:dim - 1]
DArcais(A::ℤSeq) = LinMap(DArcais, A, length(A))
DArcaisTransform(A::ℤSeq) = DArcais.(Telescope(A))

# ------------------------------------------------

function Worpitzky(n, k)
    haskey(CacheWorpitzky, (n, k)) && return CacheWorpitzky[(n, k)]
    S = sum(Eulerian(n, j)*Binomial(n-j, n-k) for j in 0:n)
    CacheWorpitzky[(n, k)] = S
end

WorpitzkyTriangle(dim) = [[Worpitzky(n, k) for k = 0:n] for n = 0:dim - 1]
Worpitzky(A::ℤSeq) = LinMap(Worpitzky, A, length(A))
WorpitzkyTransform(A::ℤSeq) = Worpitzky.(Telescope(A))

# ------------------------------------------------

"""
Return the second order Eulerian number E2(n,k).
"""
function EulerianSO2(n, k)
    n < 0 && return ZZ(0)
    haskey(CacheEulerianSO2, (n, k)) && return CacheEulerianSO2[(n, k)]
    S = EulerianSO2(n-1, k)*k + EulerianSO2(n-1, k-1)*(2*n - k)
    CacheEulerianSO2[(n, k)] = S
end

EulerianSO2Triangle(dim) = [[EulerianSO2(n, k) for k = 0:n] for n = 0:dim - 1]
EulerianSO2(A::ℤSeq) = LinMap(EulerianSO2, A, length(A))
EulerianSO2Transform(A::ℤSeq) = EulerianSO2.(Telescope(A))

# ------------------------------------------------

function Swing(n, k)
    haskey(CacheSwing, (n, k)) && return CacheSwing[(n, k)]
    p = ZZ(2)^n
    k == 0 && return p
    S = numerator(p*prod(QQ(2,j)^((-1)^(j )) for j in 1:k))
    CacheSwing[(n, k)] = S
end

SwingTriangle(dim) = [[Swing(n, k) for k = 0:n] for n = 0:dim - 1]
Swing(A::ℤSeq) = LinMap(Swing, A, length(A))
SwingTransform(A::ℤSeq) = Swing.(Telescope(A))

# ------------------------------------------------

Uni(n, k) = 1

UniTriangle(dim) = [[ZZ(1) for k = 0:n] for n = 0:dim - 1]
Uni(A::ℤSeq) = LinMap(Uni, A, length(A))
UniTransform(A::ℤSeq) = Uni.(Telescope(A))

# ------------------------------------------------

function Bessel1(n, k)
    (k == 0) && return ZZ(k^n)
    div(Factorial(2*n - k - 1),((-2)^(n - k)*Factorial(k - 1)*Factorial(n - k)))
end

Bessel1Triangle(dim) = [[Bessel1(n, k) for k = 0:n] for n = 0:dim - 1]
Bessel1(A::ℤSeq) = LinMap(Bessel1, A, length(A))
Bessel1Transform(A::ℤSeq) = Bessel1.(Telescope(A))

# ------------------------------------------------

function Hermite(n, k)
    isodd(n - k) && return ZZ(0)
    m = div(n - k, 2)
    div(Factorial(n), Factorial(k)*Factorial(m)*2^m)
end

HermiteTriangle(dim) = [[Hermite(n, k) for k = 0:n] for n = 0:dim - 1]
Hermite(A::ℤSeq) = LinMap(Hermite, A, length(A))
HermiteTransform(A::ℤSeq) = Hermite.(Telescope(A))

# ------------------------------------------------

function Trinomial(n, k)
    sum(Binomial(n, j) * Binomial(j, k-j) for j in 0:n)
end

TrinomialTriangle(dim) = [[Trinomial(n, k) for k = 0:n] for n = 0:dim - 1]
Trinomial(A::ℤSeq) = LinMap(Trinomial, A, length(A))
TrinomialTransform(A::ℤSeq) = Trinomial.(Telescope(A))

# ------------------------------------------------

bs(n) = iszero(n) ? 0 : isodd(n) ? 2 : 1
SchröderBTriangle(dim) = DelehamΔ(dim, bs, n -> 0^n)

ls(n) = isodd(n) ? 2 : 1
SchröderLTriangle(dim) = DelehamΔ(dim, ls, n -> 0^n)

"""
Recurrence for A132393, StirlingCycle numbers.
"""
function R132393(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    (n - 1) * prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A048993, StirlingSet numbers.
"""
function R048993(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    k * prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A271703, Lah numbers.
"""
function R271703(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    (k - 1 + n) * prevrow(k) + prevrow(k - 1)
end
"""
Recurrence for A094587, (Rising factorials).
"""
function R094587(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    (n - k) * prevrow(k) + prevrow(k - 1)
end

"""
Recurrence for A008279. Number of permutations of n things k at a time.
(Falling factorials)
"""
function R008279(n::Int, k::Int, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    prevrow(k) + k * prevrow(k - 1)
end

"""
xIterates over the first n rows of `A132393`.
Triangle of unsigned Stirling numbers of the first kind.
"""
I132393(n) = RecTriangle(n, R132393)
T132393(dim) = ZTri(I132393(dim))

"""
(SIGNATURES)
"""
StirlingCycleTriangle(dim) = T132393(dim)
StirlingCycle(n) = StirlingCycleTriangle(n + 1)[n + 1]
StirlingCycle(n, k) = StirlingCycle(n)[k + 1]
StirlingCycle(A::ℤSeq) = LinMap(StirlingCycle, A, length(A))
StirlingCycleTransform(A::ℤSeq) = StirlingCycle.(Telescope(A))

"""
Iterates over the first n rows of `A048993`.
Triangle of Stirling numbers of 2nd kind.
"""
I048993(n) = RecTriangle(n, R048993)
T048993(dim) = ZTri(I048993(dim))

"""
(SIGNATURES)
"""
StirlingSetTriangle(dim) = T048993(dim)
StirlingSet(n) = StirlingSetTriangle(n + 1)[n + 1]
StirlingSet(n, k) = StirlingSet(n)[k + 1]
StirlingSet(A::ℤSeq) = LinMap(StirlingSet, A, length(A))
StirlingSetTransform(A::ℤSeq) = StirlingSet.(Telescope(A))

# """
# Iterates over the first n rows of `A271703`.
# """
# Alternative definition already given above via EgfExpansion.
# I271703(n) = RecTriangle(n, R271703)
# T271703(dim) = ZTri(I271703(dim))
"""
Iterates over the first n rows of `A094587`.
"""
I094587(n) = RecTriangle(n, R094587)
T094587(dim) = ZTri(I094587(dim))
# T094587(dim) = Reverse(ZTri(I008279(dim)))

"""
Iterates over the first n rows of `A008279`.
"""
I008279(n) = RecTriangle(n, R008279)
T008279(dim) = ZTri(I008279(dim))
# Alias
FallingFactTriangle(dim) = ZTri(I008279(dim))
PermCoeffsTriangle(dim) = ZTri(I008279(dim))

function R225478(n, k, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    4 * prevrow(k - 1) + (4 * n - 1) * prevrow(k)
end
I225478(dim) = RecTriangle(dim, R225478)
T225478(dim) = ZTri(I225478(dim))

function R007318(n, k, prevrow::Function)
    (k == 0 && n == 0) && return ZZ(1)
    prevrow(k - 1) + prevrow(k)
end
I007318(dim) = RecTriangle(dim, R007318)
PascalTriangle(dim) = ZTri(I007318(dim))

function SchroederBig(n, k)
    div(
        (k + 1) * sum(Binomial(n + 1, k + j + 1) *
            Binomial(n + j, j) for j = 0:(n - k)),
        n + 1
    )
end

function SchroederBigTriangle(dim)
    [[SchroederBig(n, k) for k = 0:n] for n = 0:dim - 1]
end

SchroederBig(n) = SchroederBigTriangle(n + 1)[n + 1]
SchroederBig(A::ℤSeq) = LinMap(SchroederBig, A, length(A))
SchroederBigTransform(A::ℤSeq) = SchroederBig.(Telescope(A))

FJ(n, k) = n <= 2 ? 1 : FJ(n - 1, k) + k * FJ(n - 2, k)

"""
(SIGNATURES)
"""
JacobsthalTriangle(dim) = RiordanSquare(dim, n -> FJ(n, 2))
Jacobsthal(n) = JacobsthalTriangle(n + 1)[n + 1]
Jacobsthal(n, k) = Jacobsthal(n)[k + 1]
Jacobsthal(A::ℤSeq) = LinMap(Jacobsthal, A, length(A))
JacobsthalTransform(A::ℤSeq) = Jacobsthal.(Telescope(A))

"""
(SIGNATURES)
"""
FibonacciTriangle(dim) = RiordanSquare(dim, n -> FJ(n, 1))
Fibonacci(n) = FibonacciTriangle(n + 1)[n + 1]
Fibonacci(n, k) = Fibonacci(n)[k + 1]
Fibonacci(A::ℤSeq) = LinMap(Fibonacci, A, length(A))
FibonacciTransform(A::ℤSeq) = Fibonacci.(Telescope(A))

function EulerianNumbers(n, k)
    haskey(CacheEulerian, (n, k)) && return CacheEulerian[(n, k)]
    CacheEulerian[(n, k)] = if (k == 0) && (n >= 0)
        ZZ(1)
    elseif (k < 0) || (k > n)
        ZZ(0)
    else
        (n - k) * EulerianNumbers(n - 1, k - 1) +
        (k + 1) * EulerianNumbers(n - 1, k)
    end
end

"""
(SIGNATURES)
"""
EulerianTriangle(dim) = [[EulerianNumbers(n, k) for k = 0:n] for n = 0:dim - 1]
Eulerian(n) = EulerianTriangle(n + 1)[n + 1]
Eulerian(n, k) = Eulerian(n)[k + 1]
Eulerian(A::ℤSeq) = LinMap(Eulerian, A, length(A))
EulerianTransform(A::ℤSeq) = Eulerian.(Telescope(A))

# ------------------------------------------------

function A046802(n, k)
    sum(Binomial(n, j) * Eulerian(j, j-k) for j in k:n)
end

A046802Triangle(dim) = [[A046802(n, k) for k = 0:n] for n = 0:dim - 1]
A046802(A::ℤSeq) = LinMap(A046802, A, length(A))
A046802Transform(A::ℤSeq) = A046802.(Telescope(A))

# ------------------------------------------------

function EulerianNumbersClassic(n, k)
    haskey(CacheEulerianClassic, (n, k)) && return CacheEulerianClassic[(n, k)]
    CacheEulerianClassic[(n, k)] = if (k == n)
        ZZ(1)
    elseif (k <= 0) || (k > n)
        ZZ(0)
    else
        (n - k + 1) * EulerianNumbersClassic(n - 1, k - 1) +
        (k) * EulerianNumbersClassic(n - 1, k)
    end
end

EulerianTriangleClassic(dim) = [[EulerianNumbersClassic(n, k) for k = 0:n] for n = 0:dim - 1]
EulerianClassic(n) = EulerianTriangleClassic(n + 1)[n + 1]
EulerianClassic(n, k) = EulerianClassic(n)[k + 1]
EulerianClassic(A::ℤSeq) = LinMap(EulerianClassic, A, length(A))
EulerianTransform2(A::ℤSeq) = EulerianClassic.(Telescope(A))

function NarayanaNumbers(n::Int, k::Int)
    haskey(CacheNarayana, (n, k)) && return CacheNarayana[(n, k)]
    CacheNarayana[(n, k)] = if (k == n)
        ZZ(1)
    elseif k > n || k <= 0
        ZZ(0)
    else div(
            Binomial(ZZ(n), n - k) * Binomial(ZZ(n - 1), n - k),
            ZZ(n - k + 1)
         )
    end
end

"""
(SIGNATURES)
"""
NarayanaTriangle(dim) = [[NarayanaNumbers(n, k) for k = 0:n] for n = 0:dim]
Narayana(n) = NarayanaTriangle(n + 1)[n + 1]
Narayana(n, k) = Narayana(n)[k + 1]
Narayana(A::ℤSeq) = LinMap(Narayana, A, length(A))
NarayanaTransform(A::ℤSeq) = Narayana.(Telescope(A))

# ------------------------------------------------

A343237(n, k) = ZZ(k + 1)^(n - k + 1) - ZZ(k)^(n - k + 1)

A343237Triangle(dim) = [[A343237(n, k) for k = 0:n] for n = 0:dim - 1]
A343237(A::ℤSeq) = LinMap(A343237, A, length(A))
A343237Transform(A::ℤSeq) = A343237.(Telescope(A))

# ------------------------------------------------

"""
(SIGNATURES)
"""
RisingFactTriangle(dim) = [[rising_factorial(ZZ(n), ZZ(k)) for k = 0:n] for n = 0:dim]
RisingFact(n) = RisingFactTriangle(n + 1)[n + 1]
RisingFact(n, k) = RisingFact(n)[k + 1]
RisingFact(A::ℤSeq) = LinMap(RisingFact, A, length(A))
RisingFactTransform(A::ℤSeq) = RisingFact.(Telescope(A))

function Laguerre(n::Int)
    haskey(CacheLaguerre, n) && return CacheLaguerre[n]
    prevrow = Laguerre(n - 1)
    row = ZSeq(n + 1)
    row[n + 1] = ZZ(1)
    for k ∈ 1:n
        row[k] = (  get(prevrow, k - 1, 0)
                  + get(prevrow, k, 0) * (2 * k - 1)
                  + get(prevrow, k + 1, 0) * k^2 )
    end
    # println(row)
    CacheLaguerre[n] = row
end

Laguerre(n::Int, k::Int) = Laguerre(n)[k + 1]
Laguerre(A::ℤSeq) = LinMap(Laguerre, A, length(A))

"""
(SIGNATURES)
"""
function LaguerreTriangle(size::Int)
    length(CacheLaguerre) < size && Laguerre(size)
    [CacheLaguerre[n] for n in 0:size - 1]
end

LaguerreTransform(A::ℤSeq) = Laguerre.(Telescope(A))

function transforms(trans)
    len = 18

    println("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
    println(trans)
    println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
    println()
    Unos = [ZZ(1) for n = 0:len]
    print("UNOS ")
    Println(Unos)
    trans(Unos) |> Println

    println()
    AltUnos = [(-1)^n * ZZ(1) for n = 0:len]
    println("ALTS ")
    Println(AltUnos)
    trans(AltUnos) |> Println

    println()
    Squares = [ZZ(n^2) for n = 0:len]
    println("SQRS ")
    Println(Squares)
    trans(Squares) |> Println

    println()
    Nat = [ZZ(n) for n = 0:len]
    println("N0TS ")
    Println(Nat)
    trans(Nat) |> Println

    println()
    Nut = [ZZ(n) for n = 1:len]
    println("NATS ")
    Println(Nut)
    trans(Nut) |> Println
end

const TRIANGLES = Function[
    AitkenTriangle,
    Bessel1Triangle,
    BinomialTriangle,
    CatalanTriangle,
    DArcaisTriangle,
    DelannoyTriangle,
    EulerianTriangle,
    EulerianSO2Triangle,
    EulerSecTriangle,
    EulerTanTriangle,
    EulerTriangle,
    FallingFactTriangle,
    FibonacciTriangle,
    FineTriangle,
    FubiniTriangle,
    HermiteTriangle,
    LaguerreTriangle,
    LahTriangle,
    MotzkinTriangle,
    NarayanaTriangle,
    RencontresTriangle,
    RisingFactTriangle,
    SchröderBTriangle,
    SchröderLTriangle,
    StirlingCycleTriangle,
    StirlingSetTriangle,
    TrinomialTriangle,
    TTreeTriangle,
    UniTriangle,
    WorpitzkyTriangle
]

# START-TEST-########################################################
using Test

function test()
    @testset "Laguerre" begin
        A = ℤInt[0, 1, 6, 39, 292, 2505, 24306, 263431]

        @test Laguerre(0) == ℤInt[1]
        @test Laguerre(3) == ℤInt[6, 18, 9, 1]
        @test Laguerre(5) == ℤInt[120, 600, 600, 200, 25, 1]

        @test Laguerre(6, 1) == ZZ(4320)
        @test [Laguerre(7, k) for k in 0:7] == Laguerre(7)

        @test Laguerre(ℤInt[0, 1, 2, 3, 4, 5]) == ZZ(2505)
        @test Laguerre(ℤInt[0, 1, 4, 9, 16, 25]) == ZZ(5225)

        @test LaguerreTransform([ZZ(n) for n in 0:7]) == A

        # Test types
        @test typeof(Laguerre(3, 3))       === fmpz                   # ℤInt
        @test typeof(Laguerre(5))          === Array{fmpz,1}          # ℤSeq
        @test typeof(Laguerre(A))          === fmpz                   # ℤInt
        @test typeof(LaguerreTriangle(5))  === Array{Array{fmpz,1},1} # ℤTri
        @test typeof(LaguerreTransform(A)) === Array{fmpz,1}          # ℤSeq
    end
end

function demo()
    # isa(PascalTriangle(6), ℤTri) |> println

    FubiniPolynomial(6) |> println
    BernoulliPolynomial(6) |> println
    SchroederBigTriangle(8) |> println

    transforms(BinomialTransform)
    transforms(LahTransform)
    transforms(MotzkinTransform)
    transforms(SchroederBigTransform)

    Println.(JacobsthalTriangle(8))
    transforms(JacobsthalTransform)

    Println.(FibonacciTriangle(8))
    transforms(FibonacciTransform)

    Println.(StirlingSetTriangle(8))
    transforms(StirlingSetTransform)

    Println.(StirlingCycleTriangle(8))
    transforms(StirlingCycleTransform)

    Println.(LaguerreTriangle(8))
    transforms(LaguerreTransform)

    Println.(EulerianTriangle(8))
    #Println.(EulerianTriangle2(8))
    transforms(EulerianTransform)

   # Println.(Inverse(EulerianTriangle2(8)))

    Println.(NarayanaTriangle(8))
    transforms(NarayanaTransform)

    CatalanBallot.(0:6) |> println

    Println.(Reverse(DiagonalTriangle(FallingFactTriangle(10))))
    Println.(Reverse(FallingFactTriangle(10)))

    Println.(RisingFactTriangle(10))
    Println.(Inverse(RisingFactTriangle(10)))
    Println.(Reverse(RisingFactTriangle(10)))
    Println.(Inverse(Reverse(RisingFactTriangle(10))))

    for n in 0:7 BetaPoly(n) |> println end
    for n in 0:7 SecPoly(n) |> println end
    for n in 0:7 TanPoly(n) |> println end

end

function perf()
    # display(@benchmark LaguerreTriangle(1000))
    # display(@benchmark PolyArray(LaguerreTriangle(100)))
    # @time LaguerreTriangle(1000)
    # @time PolyArray(LaguerreTriangle(100))

    T = LaguerreTriangle(8)
    Println.(PolyTriangle(T))
    Println.(PolyArray(T))
    Println.(Inverse(PolyTriangle(T)))
    Println.(FubiniTriangle(8))
    Println.(DArcaisTriangle(8))
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
# Process an enwiki*.txt.bz2 file downloaded from Wikipedia to produce the data
# needed for Normalized Web/Wikipedia Distance (NWD).

const DefaultPathToEnwikiFile = "~/tmp/english_wikipedia"
const DefaultEnwikiFileName = "enwiki-20150205-pages-articles-multistream-index.txt.bz2"

const EnwikiFile = if length(ARGS) >= 1
  ARGS[1]
else
  DefaultPathToEnwikiFile * "/" * DefaultEnwikiFileName
end

println("enwiki file: ", EnwikiFile)


text:
A = [[+4.  +5.  -2.]
     [+7.  -1.  +2.]
     [+3.  +1.  +4.]]
b = [-14., +42., +28.]
x = [+4., -4., +5.]

lu = LUSolverLAPACK(A, b)

factorize!(lu)
solve!(lu)

@test lu.b == x


text:
function budget_constraint(ind_vec::NTuple{N,<:Integer}, mdl::Model, add_args = nothing) where N

    R,Rᴸ = mdl.Prices[1:2]
    ω, τ = mdl.AddPrices
    v = get_vars(ind_vec, mdl.GridDef)
    #l₁, l₂, d₁, d₂
    # function c = ... . If negative it means that this selection is invalid
    #TODO:fix below (lambdas are state (emp/unemp) variables)
    Δₗ = v.l₂ - (1 - δᶠ) * v.l₁ * ( (Rᴸ) / 𝜋 )
    dᴿ = v.d₂ - v.l₂ * δᶠ * δˡ ≥ 0

    𝔠 = v.d₁ * ( ρ * R / 𝜋) + ( 1 - τ ) * ω * (v.λᵢ==1) + ϐ * (v.λᵢ==2) + Δₗ -
    δᶠ * v.l₁ * ( (Rᴸ) / 𝜋 ) - v.d₂

    # C_ deposit requirement:
    c_min = max(0.0, ( ( 1 - τ ) * ω * (v.λᵢ==1) + ϐ * (v.λᵢ==2) + v.d₁ * ( ρ * R / 𝜋) -
                    - δᶠ * v.l₁ * ( (Rᴸ) / 𝜋 ) ) )
    c = min(c_min,𝔠) + ( 𝔠 - min(c_min,𝔠) ) * ξ

    tozero_l = (ind_vec[2] - 1) == ind_vec[4]
    # tozero_d = (ind_vec[1] - 1) == ind_vec[3]
    return(c,Δₗ≥0,dᴿ,tozero_l)

end


text:
# everything related to a DG semidiscretization in 1D,
# currently limited to Lobatto-Legendre nodes

# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::TreeMesh{1}, equations,
                      dg::DG, RealT, uEltype)
  # Get cells for which an element needs to be created (i.e. all leaf cells)
  leaf_cell_ids = local_leaf_cells(mesh.tree)

  elements = init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)

  interfaces = init_interfaces(leaf_cell_ids, mesh, elements)

  boundaries = init_boundaries(leaf_cell_ids, mesh, elements)

  cache = (; elements, interfaces, boundaries)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
  cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)

  return cache
end

# The methods below are specialized on the volume integral type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::TreeMesh{1}, equations,
                      volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype)
  create_cache(mesh, have_nonconservative_terms(equations), equations, volume_integral, dg, uEltype)
end

function create_cache(mesh::TreeMesh{1}, nonconservative_terms::Val{false}, equations,
                      ::VolumeIntegralFluxDifferencing, dg, uEltype)
  NamedTuple()
end

# TODO: MHD in 1D
# function create_cache(mesh::TreeMesh{1}, nonconservative_terms::Val{true}, equations,
#                       ::VolumeIntegralFluxDifferencing, dg, uEltype)
# end

function create_cache(mesh::TreeMesh{1}, equations,
                      volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
  element_ids_dg   = Int[]
  element_ids_dgfv = Int[]

  cache = create_cache(mesh, equations,
                       VolumeIntegralFluxDifferencing(volume_integral.volume_flux_dg),
                       dg, uEltype)

  A2dp1_x = Array{uEltype, 2}
  fstar1_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg)+1) for _ in 1:Threads.nthreads()]

  return (; cache..., element_ids_dg, element_ids_dgfv, fstar1_threaded)
end

function create_cache(mesh::TreeMesh{1}, equations,
                      volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype)

  A2dp1_x = Array{uEltype, 2}
  fstar1_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg)+1) for _ in 1:Threads.nthreads()]

  return (; fstar1_threaded)
end

# The methods below are specialized on the mortar type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::TreeMesh{1}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype)
  NamedTuple()
end

# TODO: Taal discuss/refactor timer, allowing users to pass a custom timer?

function rhs!(du, u, t,
              mesh::TreeMesh{1}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Reset du
  @timeit_debug timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @timeit_debug timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Prolong solution to interfaces
  @timeit_debug timer() "prolong2interfaces" prolong2interfaces!(
    cache, u, mesh, equations, dg)

  # Calculate interface fluxes
  @timeit_debug timer() "interface flux" calc_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg, cache)

  # Prolong solution to boundaries
  @timeit_debug timer() "prolong2boundaries" prolong2boundaries!(
    cache, u, mesh, equations, dg)

  # Calculate boundary fluxes
  @timeit_debug timer() "boundary flux" calc_boundary_flux!(
    cache, t, boundary_conditions, mesh, equations, dg)

  # Calculate surface integrals
  @timeit_debug timer() "surface integral" calc_surface_integral!(
    du, mesh, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
  @timeit_debug timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @timeit_debug timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  return nothing
end

function calc_volume_integral!(du, u,
                               mesh::Union{TreeMesh{1}, CurvedMesh{1}},
                               nonconservative_terms::Val{false}, equations,
                               volume_integral::VolumeIntegralWeakForm,
                               dg::DGSEM, cache)
  @unpack derivative_dhat = dg.basis

  @threaded for element in eachelement(dg, cache)
    for i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, element)

      flux1 = flux(u_node, 1, equations)
      for ii in eachnode(dg)
        integral_contribution = derivative_dhat[ii, i] * flux1
        add_to_node_vars!(du, integral_contribution, equations, dg, ii, element)
      end
    end
  end

  return nothing
end

function calc_volume_integral!(du, u,
                               mesh::TreeMesh{1},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralFluxDifferencing,
                               dg::DGSEM, cache)
  @threaded for element in eachelement(dg, cache)
    split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
  end
end

@inline function split_form_kernel!(du::AbstractArray{<:Any,3}, u,
                                    nonconservative_terms::Val{false}, equations,
                                    volume_flux, dg::DGSEM, cache,
                                    element, alpha=true)
  # true * [some floating point value] == [exactly the same floating point value]
  # This can (hopefully) be optimized away due to constant propagation.
  @unpack derivative_split = dg.basis

  # Calculate volume integral in one element
  for i in eachnode(dg)
    u_node = get_node_vars(u, equations, dg, i, element)

    # x direction
    # use consistency of the volume flux to make this evaluation cheaper
    flux1 = flux(u_node, 1, equations)
    integral_contribution = alpha * derivative_split[i, i] * flux1
    add_to_node_vars!(du, integral_contribution, equations, dg, i, element)
    # use symmetry of the volume flux for the remaining terms
    for ii in (i+1):nnodes(dg)
      u_node_ii = get_node_vars(u, equations, dg, ii, element)
      flux1 = volume_flux(u_node, u_node_ii, 1, equations)
      integral_contribution = alpha * derivative_split[i, ii] * flux1
      add_to_node_vars!(du, integral_contribution, equations, dg, i,  element)
      integral_contribution = alpha * derivative_split[ii, i] * flux1
      add_to_node_vars!(du, integral_contribution, equations, dg, ii, element)
    end
  end
end

# TODO: Taal dimension agnostic
function calc_volume_integral!(du, u,
                               mesh::TreeMesh{1},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralShockCapturingHG,
                               dg::DGSEM, cache)
  @unpack element_ids_dg, element_ids_dgfv = cache
  @unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral

  # Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α
  alpha = @timeit_debug timer() "blending factors" indicator(u, equations, dg, cache)

  # Determine element ids for DG-only and blended DG-FV volume integral
  pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)

  # Loop over pure DG elements
  @timeit_debug timer() "pure DG" @threaded for element in element_ids_dg
    split_form_kernel!(du, u, nonconservative_terms, equations, volume_flux_dg, dg, cache, element)
  end

  # Loop over blended DG-FV elements
  @timeit_debug timer() "blended DG-FV" @threaded for element in element_ids_dgfv
    alpha_element = alpha[element]

    # Calculate DG volume integral contribution
    split_form_kernel!(du, u, nonconservative_terms, equations, volume_flux_dg, dg, cache, element, 1 - alpha_element)

    # Calculate FV volume integral contribution
    fv_kernel!(du, u, equations, volume_flux_fv, dg, cache, element, alpha_element)
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_volume_integral!(du, u,
                               mesh::TreeMesh{1},
                               nonconservative_terms, equations,
                               volume_integral::VolumeIntegralPureLGLFiniteVolume,
                               dg::DGSEM, cache)
  @unpack volume_flux_fv = volume_integral

  # Calculate LGL FV volume integral
  @threaded for element in eachelement(dg, cache)
    fv_kernel!(du, u, equations, volume_flux_fv, dg, cache, element, true)
  end

  return nothing
end

@inline function fv_kernel!(du::AbstractArray{<:Any,3}, u::AbstractArray{<:Any,3},
                            equations, volume_flux_fv, dg::DGSEM, cache, element, alpha=true)
  @unpack fstar1_threaded = cache
  @unpack inverse_weights = dg.basis

  # Calculate FV two-point fluxes
  fstar1 = fstar1_threaded[Threads.threadid()]
  calcflux_fv!(fstar1, u, equations, volume_flux_fv, dg, element)

  # Calculate FV volume integral contribution
  for i in eachnode(dg)
    for v in eachvariable(equations)
      du[v, i, element] += ( alpha *
                             (inverse_weights[i] * (fstar1[v, i+1] - fstar1[v, i])) )

    end
  end

  return nothing
end

@inline function calcflux_fv!(fstar1, u::AbstractArray{<:Any,3},
                              equations, volume_flux_fv, dg::DGSEM, element)

  fstar1[:, 1,           ] .= zero(eltype(fstar1))
  fstar1[:, nnodes(dg)+1,] .= zero(eltype(fstar1))

  for i in 2:nnodes(dg)
    u_ll = get_node_vars(u, equations, dg, i-1, element)
    u_rr = get_node_vars(u, equations, dg, i,   element)
    flux = volume_flux_fv(u_ll, u_rr, 1, equations) # orientation 1: x direction
    set_node_vars!(fstar1, flux, equations, dg, i)
  end

  return nothing
end

function prolong2interfaces!(cache, u,
                             mesh::TreeMesh{1}, equations, dg::DG)
  @unpack interfaces = cache

  @threaded for interface in eachinterface(dg, cache)
    left_element  = interfaces.neighbor_ids[1, interface]
    right_element = interfaces.neighbor_ids[2, interface]

    # interface in x-direction
    for v in eachvariable(equations)
      interfaces.u[1, v, interface] = u[v, nnodes(dg), left_element]
      interfaces.u[2, v, interface] = u[v,          1, right_element]
    end
  end

  return nothing
end

function calc_interface_flux!(surface_flux_values,
                              mesh::TreeMesh{1},
                              nonconservative_terms::Val{false}, equations,
                              dg::DG, cache)
  @unpack surface_flux = dg
  @unpack u, neighbor_ids, orientations = cache.interfaces

  @threaded for interface in eachinterface(dg, cache)
    # Get neighboring elements
    left_id  = neighbor_ids[1, interface]
    right_id = neighbor_ids[2, interface]

    # Determine interface direction with respect to elements:
    # orientation = 1: left -> 2, right -> 1
    left_direction  = 2 * orientations[interface]
    right_direction = 2 * orientations[interface] - 1

    # Call pointwise Riemann solver
    u_ll, u_rr = get_surface_node_vars(u, equations, dg, interface)
    flux = surface_flux(u_ll, u_rr, orientations[interface], equations)

    # Copy flux to left and right element storage
    for v in eachvariable(equations)
      surface_flux_values[v, left_direction,  left_id]  = flux[v]
      surface_flux_values[v, right_direction, right_id] = flux[v]
    end
  end
end

# TODO: MHD in 1D
# function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1},
#                               nonconservative_terms::Val{true}, equations,
#                               dg::DG, cache)
# end

function prolong2boundaries!(cache, u,
                             mesh::TreeMesh{1}, equations, dg::DG)
  @unpack boundaries = cache
  @unpack orientations, neighbor_sides = boundaries

  @threaded for boundary in eachboundary(dg, cache)
    element = boundaries.neighbor_ids[boundary]

    # boundary in x-direction
    if neighbor_sides[boundary] == 1
      # element in -x direction of boundary
      for v in eachvariable(equations)
        boundaries.u[1, v, boundary] = u[v, nnodes(dg), element]
      end
    else # Element in +x direction of boundary
      for v in eachvariable(equations)
        boundaries.u[2, v, boundary] = u[v, 1,          element]
      end
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic,
                             mesh::TreeMesh{1}, equations, dg::DG)
  @assert isempty(eachboundary(dg, cache))
end

# TODO: Taal dimension agnostic
function calc_boundary_flux!(cache, t, boundary_condition,
                             mesh::TreeMesh{1}, equations, dg::DG)
  @unpack surface_flux_values = cache.elements
  @unpack n_boundaries_per_direction = cache.boundaries

  # Calculate indices
  lasts = accumulate(+, n_boundaries_per_direction)
  firsts = lasts - n_boundaries_per_direction .+ 1

  # Calc boundary fluxes in each direction
  for direction in eachindex(firsts)
    calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_condition,
                                     equations, dg, cache,
                                     direction, firsts[direction], lasts[direction])
  end
end

function calc_boundary_flux!(cache, t, boundary_conditions::Union{NamedTuple,Tuple},
                             mesh::TreeMesh{1}, equations, dg::DG)
  @unpack surface_flux_values = cache.elements
  @unpack n_boundaries_per_direction = cache.boundaries

  # Calculate indices
  lasts = accumulate(+, n_boundaries_per_direction)
  firsts = lasts - n_boundaries_per_direction .+ 1

  # Calc boundary fluxes in each direction
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1],
                                   equations, dg, cache, 1, firsts[1], lasts[1])
  calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2],
                                   equations, dg, cache, 2, firsts[2], lasts[2])
end

function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,3}, t,
                                          boundary_condition, equations, dg::DG, cache,
                                          direction, first_boundary, last_boundary)
  @unpack surface_flux = dg
  @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries

  @threaded for boundary in first_boundary:last_boundary
    # Get neighboring element
    neighbor = neighbor_ids[boundary]

    # Get boundary flux
    u_ll, u_rr = get_surface_node_vars(u, equations, dg, boundary)
    if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right
      u_inner = u_ll
    else # Element is on the right, boundary on the left
      u_inner = u_rr
    end
    x = get_node_coords(node_coordinates, equations, dg, boundary)
    flux = boundary_condition(u_inner, orientations[boundary], direction, x, t, surface_flux,
                              equations)

    # Copy flux to left and right element storage
    for v in eachvariable(equations)
      surface_flux_values[v, direction, neighbor] = flux[v]
    end
  end

  return nothing
end

function calc_surface_integral!(du, mesh::Union{TreeMesh{1}, CurvedMesh{1}},
                                equations, dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  @threaded for element in eachelement(dg, cache)
    for v in eachvariable(equations)
      # surface at -x
      du[v, 1,          element] -= surface_flux_values[v, 1, element] * boundary_interpolation[1,          1]
      # surface at +x
      du[v, nnodes(dg), element] += surface_flux_values[v, 2, element] * boundary_interpolation[nnodes(dg), 2]
    end
  end

  return nothing
end

function apply_jacobian!(du, mesh::Union{TreeMesh{1}, CurvedMesh{1}},
                         equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    factor = -cache.elements.inverse_jacobian[element]

    for i in eachnode(dg)
      for v in eachvariable(equations)
        du[v, i, element] *= factor
      end
    end
  end

  return nothing
end

# TODO: Taal dimension agnostic
function calc_sources!(du, u, t, source_terms::Nothing,
                       equations::AbstractEquations{1}, dg::DG, cache)
  return nothing
end

function calc_sources!(du, u, t, source_terms,
                       equations::AbstractEquations{1}, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for i in eachnode(dg)
      u_local = get_node_vars(u, equations, dg, i, element)
      x_local = get_node_coords(cache.elements.node_coordinates, equations, dg, i, element)
      du_local = source_terms(u_local, x_local, t, equations)
      add_to_node_vars!(du, du_local, equations, dg, i, element)
    end
  end

  return nothing
end


text:
export optimal_cell_size

export StructuredGrid, YeeGrid

export xnodes, ynodes, znodes,nodes
export cell_centroids, cell_sizes,cell_volumes

export Primal , Dual

include("./Grid/UnifromGridGeneration.jl")
include("./Grid/StructuredGrid.jl")
include("./Grid/YeeGrid.jl")


text:
using Test
using ContinuousObservationToyPOMDPs
using ParticleFilters
using QMDP
using POMDPs
using POMDPSimulators
using POMDPTesting

sld = SimpleLightDark()
@test isterminal(sld, sld.radius+1)
p = solve(LDHSolver(), sld)
filter = SIRParticleFilter(sld, 1000)
for (s, b, a, r, sp, o) in stepthrough(sld, p, filter, "s,b,a,r,sp,o", max_steps=10)
    @show (s, a, r, sp, o)
    @show mean(b)
end

qp = solve(QMDPSolver(), sld, verbose=true)
for (s, b, a, r, sp, o) in stepthrough(sld, qp, "s,b,a,r,sp,o", max_steps=10)
    @show (s, a, r, sp, o)
    @show b
end

m = COTigerPOMDP()
qp = solve(QMDPSolver(), m)
for (s, b, a, r, sp, o, t) in stepthrough(m, qp, "s,b,a,r,sp,o,t", max_steps=100)
    @show (t=t, s=s, a=a, r=r, sp=sp, o=o)
    @show collect(s=>pdf(b, s) for s in support(b))
end

m = TimedCOTigerPOMDP()
qp = solve(QMDPSolver(), m)
steps = collect(stepthrough(m, qp, "s,b,a,r,sp,o,t", max_steps=100))
@test length(steps) <= ContinuousObservationToyPOMDPs.horizon(m)
for (s, b, a, r, sp, o, t) in stepthrough(m, qp, "s,b,a,r,sp,o,t", max_steps=100)
    @show (t=t, s=s, a=a, r=r, sp=sp, o=o)
end

m = TimedDOTigerPOMDP()
qp = solve(QMDPSolver(), m)
steps = collect(stepthrough(m, qp, "s,b,a,r,sp,o,t", max_steps=100))
@test length(steps) <= ContinuousObservationToyPOMDPs.horizon(m)
for (s, b, a, r, sp, o, t) in stepthrough(m, qp, "s,b,a,r,sp,o,t", max_steps=100)
    @show (t=t, s=s, a=a, r=r, sp=sp, o=o)
end

trans_prob_consistency_check(COTigerPOMDP())
probability_check(DOTigerPOMDP())
# these don't work for problems with terminal states like this :'(
trans_prob_consistency_check(TimedCOTigerPOMDP())
probability_check(TimedDOTigerPOMDP())


text:
using RiskAdjustedLinearizations, Test
include(joinpath(dirname(@__FILE__), "..", "..", "examples", "crw", "crw.jl"))

# Solve model
m_crw = CoeurdacierReyWinant()
m = crw(m_crw)
try
    solve!(m, m.z, m.y, m.Ψ; algorithm = :homotopy, step = .5, verbose = :none)
catch e
    local sssout = JLD2.jldopen(joinpath(dirname(@__FILE__), "..", "..", "test", "reference/crw_sss.jld2"), "r")
    update!(ral, sssout["z_rss"], sssout["y_rss"], sssout["Psi_rss"])
end

# Verify impulse responses with a zero shock is the same as simulate with no shocks
horizon = 100
@testset "Calculate impulse responses for Coeurdacier et al. (2011)" begin

    # No shocks and start from steady state
    state1, jump1 = simulate(m, horizon)
    state2, jump2 = impulse_responses(m, horizon, 1, 0.; deviations = false)

    @test state1 ≈ state2
    @test jump1 ≈ jump2
    @test state1 ≈ repeat(m.z, 1, horizon)
    @test jump1 ≈ repeat(m.y, 1, horizon)

    # No shocks but perturb away from steady state
    state1, jump1 = simulate(m, horizon, 1.01 * m.z)
    state2, jump2 = impulse_responses(m, horizon, 1, 0., 1.01 * m.z; deviations = false)

    @test state1 ≈ state2
    @test jump1 ≈ jump2
    @test !(state1[:, 2] ≈ m.z)
    @test !(jump1[:, 2] ≈ m.y)

    # Now with shocks, from steady state
    shocks = zeros(2, horizon)
    shocks[1] = -3.

    state1, jump1 = impulse_responses(m, horizon, 1, -3.; deviations = false)
    state2, jump2 = impulse_responses(m, horizon, 1, -3., m.z; deviations = false)
    state3, jump3 = simulate(m, horizon, shocks)
    state4, jump4 = impulse_responses(m, horizon, 1, -3.; deviations = true)
    state5, jump5 = impulse_responses(m, horizon, 1, -3., m.z; deviations = true)

    @test state1 ≈ state2
    @test state1 ≈ state3
    @test state1 ≈ state4 .+ m.z
    @test state1 ≈ state5 .+ m.z
    @test jump1 ≈ jump2
    @test jump1 ≈ jump3
    @test jump1 ≈ jump4 .+ m.y
    @test jump1 ≈ jump5 .+ m.y
end