text:
"""
    ToplevelErrorReport

An interface type of error reports that JET collects while top-level concrete interpration.
All `ToplevelErrorReport` should have the following fields:
- `file::String`: the path to the file containing the interpretation context
- `line::Int`: the line number in the file containing the interpretation context

See also: [`virtual_process`](@ref), [`ConcreteInterpreter`](@ref)
"""
abstract type ToplevelErrorReport end

# `ToplevelErrorReport` interface
function Base.getproperty(er::ToplevelErrorReport, sym::Symbol)
    return if sym === :file
        getfield(er, sym)::String
    elseif sym === :line
        getfield(er, sym)::Int
    else
        getfield(er, sym) # fallback
    end
end

"""
Configurations for top-level analysis.
These configurations will be active for all the top-level entries explained in the
[top-level analysis entry points](@ref jetanalysis-toplevel-entry) section.

---
- `context::Bool = Main` \\
  The module context in which the top-level execution will be simulated.

  This configuration can be useful when you just want to analyze a submodule, without
  starting entire analysis from the root module.
  For example, we can analyze `Base.Math` like below:
  ```julia-repl
  julia> report_file(JET.fullbasepath("math.jl");
                     context = Base,                  # `Base.Math`'s root module
                     analyze_from_definitions = true, # there're only definitions in `Base`
                     )
  ```

  Note that this module context will be virtualized by default so that JET can repeat analysis
  in the same session without having "invalid redefinition of constant ..." error etc.
  In other word, JET virtualize the module context of `context` and make sure the original
  module context isn't polluted by JET.
---
- `target_defined_modules::Bool = false` \\
  If `true`, automatically set the [`target_modules`](@ref result-config) configuration so that
  JET filters out errors that are reported within modules that JET doesn't analyze directly.
---
- `analyze_from_definitions::Bool = false` \\
  If `true`, JET will start analysis using signatures of top-level definitions (e.g. method signatures),
  after the top-level interpretation has been done (unless no serious top-level error has
  happened, like errors involved within a macro expansion).

  This can be handy when you want to analyze a package, which usually contains only definitions
  but not their usages (i.e. top-level callsites).
  With this option, JET can enter analysis just with method or type definitions, and we don't
  need to pass a file that uses the target package.

  !!! warning
      This feature is very experimental at this point, and you may face lots of false positive
      errors, especially when trying to analyze a big package with lots of dependencies.
      If a file that contains top-level callsites (e.g. `test/runtests.jl`) is available,
      JET analysis using the file is generally preferred, since analysis entered from
      concrete call sites will produce more accurate results than analysis entered from
      (maybe not concrete-typed) method signatures.

  Also see: [`report_file`](@ref), [`report_and_watch_file`](@ref)
---
- `concretization_patterns::Vector{<:Any} = Expr[]` \\
  Specifies a customized top-level code concretization strategy.

  When analyzing a top-level code, JET first splits the entire code into appropriate units
  of code (i.e. code blocks), and then iterate a virtual top-level code execution process
  on each code block in order to simulate Julia's sequential top-level code execution.
  In virtual code execution, JET will selectively interpret "top-level definitions" (like a function definition)
  just like Julia's top-level code execution, while it tries to avoid executing any other
  parts of code like function calls and leaves them to succeeding static analysis by
  abstract interpretation.

  However, currently, JET doesn't track the "inter-code-block" level code dependencies, and
  so the selective interpretation of top-level definitions can fail if it needs an access to
  global variables defined in other code blocks that are not actually interpreted (i.e. "concretized")
  but just left for abstract interpreation (i.e. "abstracted").

  For example, the issue happens when your macro accesses to a global variable during its expansion, e.g.:
  > test/fixtures/concretization_patterns.jl
  $(let
      text = read(normpath(@__DIR__, "..", "..", "test", "fixtures", "concretization_patterns.jl"), String)
      lines = split(text, '\n')
      pushfirst!(lines, "```julia"); push!(lines, "```")
      join(lines, "\n  ")
  end)

  To circumvent this issue, JET offers the `concretization_patterns::Vector{<:Any}` configuration,
  which allows us to customize JET's top-level code concretization strategy.
  `concretization_patterns` specifies the _patterns of code_ that should be concretized.
  To put in other word, when JET sees a code that matches any of code patterns specified by
  an user, JET will try to interpret and concretize the code, regardless of whether JET's
  code selection logic decides to concretize it or not.

  JET uses [MacroTools.jl's expression pattern match](https://fluxml.ai/MacroTools.jl/stable/pattern-matching/),
  and we can specify whatever code pattern expected by `MacroTools.@capture` macro.
  For example, in order to solve the issue explained above, we can have:
  ```julia
  concretization_patterns = [:(const GLOBAL_CODE_STORE = Dict())]
  ```
  Then `GLOBAL_CODE_STORE` will just be concretized and so any top-level error won't happen
  at the macro expansion.

  Although configuring `concretization_patterns` properly could be really tricky, we can
  effectively debug JET's top-level code concretization plan using [`toplevel_logger`](@ref JETLogger)
  configuration with the logging level above than `$DEBUG_LOGGER_LEVEL` ("debug") level.
  With the `toplevel_logger` configuration, we can see:
  - which code is matched with `concretization_patterns` and forcibly concretized
  - which code is selected to be concretized or not by JET's code selection logic:
    where `t`-annotated statements are concretized while `f`-annotated statements are abstracted
    and left abstract interpretation
  ```julia-repl
  julia> report_file("test/fixtures/concretization_patterns.jl";
                     concretization_patterns = [:(const GLOBAL_CODE_STORE = Dict())],
                     toplevel_logger = IOContext(stdout, :JET_LOGGER_LEVEL => 1))
  ```
  ```
  [toplevel-debug] virtualized the context of Main (took 0.003 sec)
  [toplevel-debug] entered into test/fixtures/concretization_patterns.jl
  [toplevel-debug] concretization pattern `const GLOBAL_CODE_STORE = Dict()` matched `const GLOBAL_CODE_STORE = Dict()` at test/fixtures/concretization_patterns.jl:2
  [toplevel-debug] concretization plan at test/fixtures/concretization_patterns.jl:4:
  1 f 1 ─      \$(Expr(:thunk, CodeInfo(
      @ none within `top-level scope`
  1 ─     return \$(Expr(:method, Symbol("@with_code_record")))
  )))
  2 t │        \$(Expr(:method, Symbol("@with_code_record")))
  3 t │   %3 = Core.Typeof(var"@with_code_record")
  4 t │   %4 = Core.svec(%3, Core.LineNumberNode, Core.Module, Core.Any)
  5 t │   %5 = Core.svec()
  6 t │   %6 = Core.svec(%4, %5, \$(QuoteNode(:(#= test/fixtures/concretization_patterns.jl:4 =#))))
  7 t │        \$(Expr(:method, Symbol("@with_code_record"), :(%6), CodeInfo(
      @ test/fixtures/concretization_patterns.jl:5 within `none`
  1 ─      \$(Expr(:meta, :nospecialize, :(a)))
  │        Base.setindex!(GLOBAL_CODE_STORE, a, __source__)
  │   @ test/fixtures/concretization_patterns.jl:6 within `none`
  │   %3 = esc(a)
  └──      return %3
  )))
  8 f └──      return var"@with_code_record"
  [toplevel-debug] concretization plan at test/fixtures/concretization_patterns.jl:11:
  1 f 1 ─      \$(Expr(:thunk, CodeInfo(
      @ none within `top-level scope`
  1 ─     return \$(Expr(:method, :foo))
  )))
  2 t │        \$(Expr(:method, :foo))
  3 t │   %3 = Core.Typeof(foo)
  4 t │   %4 = Core.svec(%3, Core.Any)
  5 t │   %5 = Core.svec()
  6 t │   %6 = Core.svec(%4, %5, \$(QuoteNode(:(#= test/fixtures/concretization_patterns.jl:11 =#))))
  7 t │        \$(Expr(:method, :foo, :(%6), CodeInfo(
      @ test/fixtures/concretization_patterns.jl:11 within `none`
  1 ─ %1 = identity(a)
  └──      return %1
  )))
  8 f └──      return foo
  [toplevel-debug] concretization plan at test/fixtures/concretization_patterns.jl:13:
  1 f 1 ─ %1 = foo(10)
  2 f └──      return %1
  [toplevel-debug]  exited from test/fixtures/concretization_patterns.jl (took 0.032 sec)
  ```

  Also see: [JET's logging configurations](@ref logging-config), [`virtual_process`](@ref).
---
- `virtualize::Bool = true` \\
  When `true`, JET will virtualize the given root module context.

  This configuration is supposed to be used only for testing or debugging.
  See [`virtualize_module_context`](@ref) for the internal.
---
"""
struct ToplevelConfig{CP<:Any}
    context::Module
    analyze_from_definitions::Bool
    concretization_patterns::Vector{CP}
    virtualize::Bool
    @jetconfigurable function ToplevelConfig(; context::Module                    = Main,
                                               analyze_from_definitions::Bool     = false,
                                               concretization_patterns::Vector{T} = Expr[],
                                               virtualize::Bool                   = true,
                                               ) where {T<:Any}
        if typeintersect(T, Expr) === Bottom
            CP = tmerge(T, Expr)
            concretization_patterns = convert(Vector{CP}, concretization_patterns)
        else
            CP = T
        end
        push!(concretization_patterns,
              # `@enum` macro is fairly complex and especially the `let insts = (Any[ $(esc(typename))(v) for v in $values ]...,)`
              # (adapted from https://github.com/JuliaLang/julia/blob/e5d7ef01b06f44cb75c871e54c81eb92eceae738/base/Enums.jl#L198)
              # part requires the statement selection logic to choose statements that only
              # pushes elements (`v`) into a slot representing an array (`insts`),
              # which is very hard to be generalized;
              # here we add them as pre-defined concretization patterns and make sure
              # false postive top-level errors won't happen by the macro expansion
              :(@enum(args__)), :(Base.@enum(args__)),
              )
        concretization_patterns = CP[striplines(normalise(x)) for x in concretization_patterns]
        return new{CP}(context,
                       analyze_from_definitions,
                       concretization_patterns,
                       virtualize,
                       )
    end
end

const Actual2Virtual = Pair{Module,Module}

"""
    res::VirtualProcessResult

- `res.included_files::Set{String}`: files that have been analyzed
- `res.defined_modules::Set{Module}`: module contexts created while this top-level analysis
- `res.toplevel_error_reports::Vector{ToplevelErrorReport}`: toplevel errors found during the
    text parsing or partial (actual) interpretation; these reports are "critical" and should
    have precedence over `inference_error_reports`
- `res.inference_error_reports::Vector{InferenceErrorReport}`: possible error reports found
    by `AbstractAnalyzer`
- `res.toplevel_signatures`: signatures of methods defined within the analyzed files
- `res.actual2virtual::$Actual2Virtual`: keeps actual and virtual module
"""
struct VirtualProcessResult
    included_files::Set{String}
    defined_modules::Set{Module}
    toplevel_error_reports::Vector{ToplevelErrorReport}
    inference_error_reports::Vector{InferenceErrorReport}
    toplevel_signatures::Vector{Type}
    actual2virtual::Union{Actual2Virtual,Nothing}
end

function VirtualProcessResult(actual2virtual, context)
    return VirtualProcessResult(Set{String}(),
                                Set{Module}((context,)),
                                ToplevelErrorReport[],
                                InferenceErrorReport[],
                                Type[],
                                actual2virtual,
                                )
end

"""
    virtual_process(s::AbstractString,
                    filename::AbstractString,
                    analyzer::AbstractAnalyzer,
                    config::ToplevelConfig,
                    ) -> res::VirtualProcessResult

Simulates Julia's toplevel execution and collects error points, and finally returns $(@doc VirtualProcessResult)

This function first parses `s::AbstractString` into `toplevelex::Expr` and then iterate the
following steps on each code block (`blk`) of `toplevelex`:
1. if `blk` is a `:module` expression, recusively enters analysis into an newly defined
   virtual module
2. `lower`s `blk` into `:thunk` expression `lwr` (macros are also expanded in this step)
3. if the context module is virtualized, replaces self-references of the original context
   module with virtualized one: see `fix_self_references`
4. `ConcreteInterpreter` partially interprets some statements in `lwr` that should not be
   abstracted away (e.g. a `:method` definition); see also [`partially_interpret!`](@ref)
5. finally, `AbstractAnalyzer` analyzes the remaining statements by abstract interpretation

!!! warning
    In order to process the toplevel code sequentially as Julia runtime does, `virtual_process`
    splits the entire code, and then iterate a simulation process on each code block.
    With this approach, we can't track the inter-code-block level dependencies, and so a
    partial interpretation of toplevle definitions will fail if it needs an access to global
    variables defined in other code blocks that are not interpreted but just abstracted.
    We can circumvent this issue using JET's `concretization_patterns` configuration, which
    allows us to customize JET's concretization strategy.
    See [`ToplevelConfig`](@ref) for more details.
"""
function virtual_process(x::Union{AbstractString,Expr},
                         filename::AbstractString,
                         analyzer::AbstractAnalyzer,
                         config::ToplevelConfig,
                         )
    if config.virtualize
        actual  = config.context

        start = time()
        virtual = virtualize_module_context(actual)
        with_toplevel_logger(analyzer) do @nospecialize(io)
            sec = round(time() - start; digits = 3)
            println(io, "virtualized the context of $actual (took $sec sec)")
        end

        actual2virtual = Actual2Virtual(actual, virtual)
        context = virtual
    else
        actual2virtual = nothing
        context = config.context
    end

    res = _virtual_process!(x,
                            filename,
                            analyzer,
                            config,
                            context,
                            VirtualProcessResult(actual2virtual, context),
                            )

    # analyze collected signatures unless critical error happened
    if config.analyze_from_definitions && isempty(res.toplevel_error_reports)
        analyze_from_definitions!(analyzer, res)
    end

    # TODO want to do in `analyze_from_definitions!` ?
    unique!(aggregation_policy(analyzer), res.inference_error_reports)

    return res
end

"""
    virtualize_module_context(actual::Module)

HACK: Returns a module where the context of `actual` is virtualized.

The virtualization will be done by 2 steps below:
1. loads the module context of `actual` into a sandbox module, and export the whole context from there
2. then uses names exported from the sandbox

This way, JET's runtime simulation in the virtual module context will be able to define
a name that is already defined in `actual` without causing
"cannot assign a value to variable ... from module ..." error, etc.
It allows JET to virtualize the context of already-existing module other than `Main`.

!!! warning "TODO"
    Currently this function relies on `Base.names`, and thus it can't restore the `using`ed
    names.
"""
function virtualize_module_context(actual::Module)
    modpath = split_module_path(actual)

    if length(modpath) ≥ 2
        if modpath[1] === :Main && modpath[2] === :anonymous
            error(ArgumentError("can't virtualize an anonymous module"))
        end
    end

    virtual = gen_virtual_module(actual)
    sandbox = gen_virtual_module(actual; name = :JETSandboxModule)

    uex = Expr(:(:), Expr(:., :., :., modpath...))
    usage = Expr(:using, uex)
    exprt = Expr(:export)
    unames = uex.args
    enames = exprt.args
    for n in names(actual; all = true, imported = true)
        isdefined(sandbox, n) && continue
        isdefined(actual, n) && push!(unames, Expr(:., n)) # an exported name can be undefined, and `using` of it will throw otherwise

        push!(enames, n)
    end
    Core.eval(sandbox, usage)
    Core.eval(sandbox, exprt)

    usage = Expr(:using, Expr(:., :., :., split_module_path(sandbox)...))
    Core.eval(virtual, usage)

    return virtual
end

const VIRTUAL_MODULE_NAME = :JETVirtualModule

gen_virtual_module(root = Main; name = VIRTUAL_MODULE_NAME) =
    Core.eval(root, :(module $(gensym(name)) end))::Module

# NOTE when `@generated` function has been defined, signatures of both its entry and
# generator should have been collected, and we will just analyze them separately
# if code generation has failed given the entry method signature, the overload of
# `InferenceState(..., ::AbstractAnalyzer)` will collect `GeneratorErrorReport`
function analyze_from_definitions!(analyzer::AbstractAnalyzer, res::VirtualProcessResult)
    n = length(res.toplevel_signatures)
    succeeded = 0
    for (i, tt) in enumerate(res.toplevel_signatures)
        mms = _methods_by_ftype(tt, -1, get_world_counter())
        if isa(mms, Vector{Any})
            filter!(mm::MethodMatch->mm.spec_types===tt, mms)
            if length(mms) == 1
                succeeded += 1
                with_toplevel_logger(analyzer; pre=clearline) do @nospecialize(io)
                    (i == n ? println : print)(io, "analyzing from top-level definitions ... $succeeded/$n")
                end
                analyzer = AbstractAnalyzer(analyzer, _CONCRETIZED, _TOPLEVELMOD)
                analyzer, result = analyze_method!(
                    analyzer, (first(mms)::MethodMatch).method;
                    # JETAnalyzer{BasicPass}: don't report errors unless this frame is concrete
                    set_entry = false)
                append!(res.inference_error_reports, get_reports(analyzer, result))
                continue
            end
        end
        # something went wrong
        with_toplevel_logger(analyzer, ≥(DEBUG_LOGGER_LEVEL); pre=clearline) do @nospecialize(io)
            println(io, "couldn't find a single method matching the signature `$tt`")
        end
    end
    return nothing
end
clearline(io) = print(io, '\r')

function _virtual_process!(s::AbstractString,
                           filename::AbstractString,
                           analyzer::AbstractAnalyzer,
                           config::ToplevelConfig,
                           context::Module,
                           res::VirtualProcessResult,
                           )
    start = time()

    with_toplevel_logger(analyzer) do @nospecialize(io)
        println(io, "entered into $filename")
    end

    push!(res.included_files, filename)

    toplevelex = parse_input_line(s; filename)

    if isexpr(toplevelex, (:error, :incomplete))
        # if there's any syntax error, try to identify all the syntax error location
        append!(res.toplevel_error_reports, collect_syntax_errors(s, filename))
    elseif isnothing(toplevelex)
        # just return if there is nothing to analyze
    else
        res = _virtual_process!(toplevelex, filename, analyzer, config, context, res)
    end

    with_toplevel_logger(analyzer) do @nospecialize(io)
        sec = round(time() - start; digits = 3)
        println(io, " exited from $filename (took $sec sec)")
    end

    return res
end

function _virtual_process!(toplevelex::Expr,
                           filename::AbstractString,
                           analyzer::AbstractAnalyzer,
                           config::ToplevelConfig,
                           context::Module,
                           res::VirtualProcessResult,
                           force_concretize::Bool = false,
                           )
    @assert isexpr(toplevelex, :toplevel)

    local lnn::LineNumberNode = LineNumberNode(0, filename)

    function macroexpand_err_handler(@nospecialize(err), st)
        report = is_missing_concretization(err) ?
                 MissingConcretization(err, st, filename, lnn.line) :
                 ActualErrorWrapped(err, st, filename, lnn.line)
        push!(res.toplevel_error_reports, report)
        return nothing
    end
    # `scrub_offset = 4` corresponds to `with_err_handling` -> `f` -> `macroexpand` -> kwfunc (`macroexpand`)
    macroexpand_with_err_handling(mod, x) = with_err_handling(macroexpand_err_handler, #= scrub_offset =# 4) do
        # XXX we want to non-recursive, sequential partial macro expansion here, which allows
        # us to collect more fine-grained error reports within macro expansions
        # but it can lead to invalid macro hygiene escaping because of https://github.com/JuliaLang/julia/issues/20241
        return macroexpand(mod, x; recursive = true #= but want to use `false` here =#)
    end
    function eval_err_handler(@nospecialize(err), st)
        report = is_missing_concretization(err) ?
                 MissingConcretization(err, st, filename, lnn.line) :
                 ActualErrorWrapped(err, st, filename, lnn.line)
        push!(res.toplevel_error_reports, report)
        return nothing
    end
    # `scrub_offset = 3` corresponds to `with_err_handling` -> `f` -> `eval`
    eval_with_err_handling(mod, x) = with_err_handling(eval_err_handler, #= scrub_offset =# 3) do
        return Core.eval(mod, x)
    end
    function lower_err_handler(@nospecialize(err), st)
        report = is_missing_concretization(err) ?
                 MissingConcretization(err, st, filename, lnn.line) :
                 ActualErrorWrapped(err, st, filename, lnn.line)
        push!(res.toplevel_error_reports, report)
        return nothing
    end
    # `scrub_offset = 3` corresponds to `with_err_handling` -> `f` -> `lower`
    lower_with_err_handling(mod, x) = with_err_handling(lower_err_handler, #= scrub_offset =# 3) do
        lwr = lower(mod, x)

        # here we should capture syntax errors found during lowering
        if isexpr(lwr, :error)
            msg = first(lwr.args)
            push!(res.toplevel_error_reports, SyntaxErrorReport("syntax: $msg", filename, lnn.line))
            return nothing
        end

        return lwr
    end

    # transform, and then analyze sequentially
    # IDEA the following code has some of duplicated work with `JuliaInterpreter.ExprSpliter` and we may want to factor them out
    exs = push_vex_stack!(VExpr[], toplevelex, force_concretize)
    while !isempty(exs)
        (; x, force_concretize) = pop!(exs)

        # with_toplevel_logger(analyzer, ≥(DEBUG_LOGGER_LEVEL)) do @nospecialize(io)
        #     println(io, "analyzing ", x)
        # end

        # update line info
        if isa(x, LineNumberNode)
            lnn = x
            continue
        end

        # apply user-specified concretization strategy, which is configured as expression
        # pattern match on surface level AST code representation; if any of the specified
        # patterns matches `x`, JET just concretizes everything involved with it
        # since patterns are expected to work on surface level AST, we should configure it
        # here before macro expansion and lowering
        if !force_concretize
            for pat in config.concretization_patterns
                if @capture(x, $pat)
                    with_toplevel_logger(analyzer, ≥(DEBUG_LOGGER_LEVEL)) do @nospecialize(io)
                        line, file = lnn.line, lnn.file
                        x′ = striplines(normalise(x))
                        println(io, "concretization pattern `$pat` matched `$x′` at $file:$line")
                    end
                    force_concretize = true
                    break
                end
            end
        end

        # we will end up lowering `x` later, but special case `macrocall`s and expand it here
        # this is because macros can arbitrarily generate `:toplevel` and `:module` expressions
        if isexpr(x, :macrocall)
            newx = macroexpand_with_err_handling(context, x)

            # if any error happened during macro expansion, bail out now and continue
            isnothing(newx) && continue

            # special case and flatten the resulting expression expanded from `@doc` macro
            # the macro expands to a block expression and so it makes it difficult to specify
            # concretization pattern correctly since `@doc` macro is attached implicitly
            if first(x.args) === GlobalRef(Core, Symbol("@doc"))
                # `@doc` macro usually produces :block expression, but may also produce :toplevel
                # one when attached to a module expression
                @assert isexpr(newx, :block) || isexpr(newx, :toplevel)
                push_vex_stack!(exs, newx::Expr, force_concretize)
            else
                push!(exs, VExpr(newx, force_concretize))
            end

            continue
        end

        # flatten container expression
        if isexpr(x, :toplevel)
            push_vex_stack!(exs, x, force_concretize)
            continue
        end

        # handle `:module` definition and module usage;
        # should happen here because modules need to be loaded sequentially while
        # "toplevel definitions" inside of the loaded modules shouldn't be evaluated in a
        # context of `context` module

        if isexpr(x, :module)
            newblk = x.args[3]
            @assert isexpr(newblk, :block)
            newtoplevelex = Expr(:toplevel, newblk.args...)

            x.args[3] = Expr(:block) # empty module's code body
            newcontext = eval_with_err_handling(context, x)

            isnothing(newcontext) && continue # error happened, e.g. duplicated naming

            newmod = newcontext::Module
            push!(res.defined_modules, newmod)
            _virtual_process!(newtoplevelex, filename, analyzer, config, newmod, res, force_concretize)

            continue
        end

        # can't wrap `:global` declaration into a block
        if isexpr(x, :global)
            eval_with_err_handling(context, x)
            continue
        end

        blk = Expr(:block, lnn, x) # attach current line number info
        lwr = lower_with_err_handling(context, blk)

        isnothing(lwr) && continue # error happened during lowering
        isexpr(lwr, :thunk) || continue # literal

        src = first((lwr::Expr).args)::CodeInfo

        fix_self_references!(res.actual2virtual, src)

        interp = ConcreteInterpreter(filename,
                                     lnn,
                                     eval_with_err_handling,
                                     context,
                                     analyzer,
                                     config,
                                     res,
                                     )
        if force_concretize
            JuliaInterpreter.finish!(interp, Frame(context, src), true)
            continue
        end
        concretized = partially_interpret!(interp, context, src)

        # bail out if nothing to analyze (just a performance optimization)
        all(concretized) && continue

        analyzer = AbstractAnalyzer(analyzer, concretized, context)

        _, result = analyze_toplevel!(analyzer, src)

        append!(res.inference_error_reports, get_reports(analyzer, result)) # collect error reports
    end

    return res
end

struct VExpr
    x
    force_concretize::Bool
    VExpr(@nospecialize(x), force_concretize::Bool) = new(x, force_concretize)
end

function push_vex_stack!(exs::Vector{VExpr}, newex::Expr, force_concretize::Bool)
    nargs = length(newex.args)
    for i in 0:(nargs-1)
        push!(exs, VExpr(newex.args[nargs-i], force_concretize))
    end
    return exs
end

function split_module_path(m::Module)
    ret = Symbol[]
    split_module_path!(m, ret)
    return ret
end
function split_module_path!(m::Module, ret)
    mp = parentmodule(m)
    if m === Main || m === Base || m === Core || mp === m
        push!(ret, nameof(m))
        return
    end
    split_module_path!(mp, ret)
    push!(ret, nameof(m))
end

# if virtualized, replace self references of `actualmod` with `virtualmod` (as is)
fix_self_references!(::Nothing, x) = return nothing
function fix_self_references!((actualmod, virtualmod)::Actual2Virtual, x)
    actualmodsym   = Symbol(actualmod)
    virtualmodsyms = split_module_path(virtualmod)

    # we can't use `Module` object in module usage expression, so we need to replace it with
    # its symbolic representation
    function fix_modpath(modpath)
        ret = Symbol[]
        for s in modpath
            if s === virtualmod # postwalk, so we compare to `virtualmod`
                push!(ret, virtualmodsyms...)
            else
                push!(ret, s)
            end
        end
        return ret
    end

    fix_module_usage(head, modpath) =
        Expr(head, Expr(:., fix_modpath(modpath)...))
    fix_module_usage_specific(head, modpath, names) =
        Expr(head, Expr(:(:), Expr(:., fix_modpath(modpath)...), Expr(:., names...)))

    function fix_simple_module_usage(usage)
        return @capture(usage, import modpath__) ? fix_module_usage(:import, modpath) :
               @capture(usage, using modpath__) ? fix_module_usage(:using, modpath) :
               @capture(usage, import modpath__: names__) ? fix_module_usage_specific(:import, modpath, names) :
               @capture(usage, using modpath__: names__) ? fix_module_usage_specific(:using, modpath, names) :
               usage # :export
    end

    function fix_module_usage(usage)
        head = usage.head
        ret = Expr(head)
        for simple in to_simple_module_usages(usage)
            fixed = fix_simple_module_usage(simple)
            @assert fixed.head === head && length(fixed.args) == 1
            push!(ret.args, first(fixed.args))
        end
        return ret
    end

    return postwalk_and_transform!(x) do x, scope
        if ismoduleusage(x)
            return fix_module_usage(x)
        elseif x === actualmodsym
            return virtualmod
        end
        return x
    end
end

postwalk_and_transform!(f, x, scope = Symbol[]) =
    walk_and_transform!(x, x->postwalk_and_transform!(f, x, scope), f, scope)
prewalk_and_transform!(f, x, scope = Symbol[]) =
    walk_and_transform!(f(x, scope), x->prewalk_and_transform!(f, x, scope), (x,scope)->x, scope)

walk_and_transform!(@nospecialize(x), inner, outer, scope) = outer(x, scope)
function walk_and_transform!(node::GotoIfNot, inner, outer, scope)
    push!(scope, :gotoifnot)
    cond = inner(walk_and_transform!(node.cond, inner, outer, scope))
    dest = inner(walk_and_transform!(node.dest, inner, outer, scope))
    pop!(scope)
    return outer(GotoIfNot(cond, dest), scope)
end
function walk_and_transform!(ex::Expr, inner, outer, scope)
    push!(scope, ex.head)
    for (i, x) in enumerate(ex.args)
        ex.args[i] = inner(walk_and_transform!(x, inner, outer, scope))
    end
    pop!(scope)
    return outer(ex, scope)
end
function walk_and_transform!(src::CodeInfo, inner, outer, scope)
    for (i, x) in enumerate(src.code)
        src.code[i] = inner(walk_and_transform!(x, inner, outer, scope))
    end
    return outer(src, scope)
end

# wraps an error that might happen because of inappropriate top-level code abstraction
struct MissingConcretization <: ToplevelErrorReport
    err
    st::Base.StackTraces.StackTrace
    file::String
    line::Int
end

# wraps general errors from actual execution
struct ActualErrorWrapped <: ToplevelErrorReport
    err
    st::Base.StackTraces.StackTrace
    file::String
    line::Int

    function ActualErrorWrapped(@nospecialize(err), st, file, line)
        if isa(err, ErrorException) && startswith(err.msg, "syntax: ")
            # forward syntax error
            return SyntaxErrorReport(err.msg, file, line)
        end
        return new(err, st, file, line)
    end
end

"""
    ConcreteInterpreter

The trait to inject code into JuliaInterpreter's interpretation process; JET.jl overloads:
- `JuliaInterpreter.step_expr!` to add error report pass for module usage expressions and
  support package analysis
- `JuliaInterpreter.evaluate_call_recurse!` to special case `include` calls
- `JuliaInterpreter.handle_err` to wrap an error happened during interpretation into
  `ActualErrorWrapped`
"""
struct ConcreteInterpreter{Analyzer<:AbstractAnalyzer,CP}
    filename::String
    lnn::LineNumberNode
    eval_with_err_handling::Function
    context::Module
    analyzer::Analyzer
    config::ToplevelConfig{CP}
    res::VirtualProcessResult
end

"""
    partially_interpret!(interp::ConcreteInterpreter, mod::Module, src::CodeInfo)

Partially interprets statements in `src` using JuliaInterpreter.jl:
- concretizes "toplevel definitions", i.e. `:method`, `:struct_type`, `:abstract_type` and
  `:primitive_type` expressions and their dependencies
- concretizes user-specified toplevel code (see [`ToplevelConfig`](@ref))
- directly evaluates module usage expressions and report error of invalid module usages
  (TODO: enter into the loaded module and keep JET analysis)
- special-cases `include` calls so that top-level analysis recursively enters the included file
"""
function partially_interpret!(interp::ConcreteInterpreter, mod::Module, src::CodeInfo)
    concretize = select_statements(src)

    with_toplevel_logger(interp.analyzer, ≥(DEBUG_LOGGER_LEVEL)) do @nospecialize(io)
        line, file = interp.lnn.line, interp.lnn.file
        println(io, "concretization plan at $file:$line:")
        print_with_code(io, src, concretize)
    end

    # NOTE if `JuliaInterpreter.optimize!` may modify `src`, `src` and `concretize` can be inconsistent
    # here we create `JuliaInterpreter.Frame` by ourselves disabling the optimization (#277)
    frame = Frame(mod, src; optimize=false)
    @assert length(frame.framecode.src.code) == length(concretize)
    selective_eval_fromstart!(interp, frame, concretize, #= istoplevel =# true)

    return concretize
end

# select statements that should be concretized, and actually interpreted rather than abstracted
function select_statements(src::CodeInfo)
    stmts = src.code
    edges = CodeEdges(src)

    concretize = fill(false, length(stmts))

    select_direct_requirement!(concretize, stmts, edges)

    select_dependencies!(concretize, src, edges)

    return concretize
end

function select_direct_requirement!(concretize, stmts, edges)
    for (i, stmt) in enumerate(stmts)
        if begin
                ismethod(stmt)      || # don't abstract away method definitions
                istypedef(stmt)     || # don't abstract away type definitions
                ismoduleusage(stmt)    # module usages are handled by `ConcreteInterpreter`
            end
            concretize[i] = true
            continue
        end

        if isexpr(stmt, :(=))
            lhs, rhs = stmt.args
            stmt = rhs
        end
        if isexpr(stmt, :call)
            if is_known_call(stmt, :include)
                # `include` calls are special cased
                select_stmt_requirement!(concretize, i, edges)
                continue
            elseif is_known_call(stmt, :eval)
                # analysis of `eval` calls are difficult, let's give up it and just evaluate
                # toplevel `eval` calls; they may contain toplevel definitions
                select_stmt_requirement!(concretize, i, edges)
                continue
            end
        end
    end
end

# adapted from
# - https://github.com/timholy/Revise.jl/blob/266ed68d7dd3bea67c39f96513cda30bbcd7d441/src/lowered.jl#L53
# - https://github.com/timholy/Revise.jl/blob/266ed68d7dd3bea67c39f96513cda30bbcd7d441/src/lowered.jl#L87-L88
function is_known_call(stmt::Expr, func::Symbol)
    f = stmt.args[1]
    isa(f, Symbol) && f === func && return true
    isa(f, GlobalRef) && f.name === :eval && return true
    callee_matches(f, Base, :getproperty) && is_quotenode_egal(stmt.args[end], func) && return true
    callee_matches(f, Core, :getproperty) && is_quotenode_egal(stmt.args[end], func) && return true
    callee_matches(f, Core.Compiler, :getproperty) && is_quotenode_egal(stmt.args[end], func) && return true
    return false
end

# statement `idx` may be the equivalent of `f = known_call`, so require each
# stmt that calls `known_call` via `f(expr)`
function select_stmt_requirement!(concretize, idx, edges)
    concretize[edges.succs[idx]] .= true
    concretize[idx] = true
end

# TODO implement these prototypes and support a following pattern:
# N = 10
# let tpl = ([i for i in 1:N]...,)
#     @eval gettpl() = $tpl # `tpl` here should be fully concretized
# end
function find_iterblocks(src) end
function add_iterblocks!(cocretize, src, edges, iterblocks) end

# implementation of https://github.com/aviatesk/JET.jl/issues/196
function select_dependencies!(concretize, src, edges)
    debug = false
    debug && println("initially selected:", findall(concretize))

    # find statement sets that come from iteration/iterator protocol
    # TODO iterblocks = find_iterblocks(src)

    # We'll mostly use generic graph traversal to discover all the lines we need,
    # but structs are in a bit of a different category (especially on Julia 1.5+).
    # It's easiest to discover these at the beginning.
    typedefs = find_typedefs(src)

    changed = true
    while changed
        changed = false

        # track SSA predecessors of initial requirements
        changed |= add_ssa_preds!(concretize, src, edges, ())

        # add some domain-specific information
        # TODO changed |= add_iterblocks!(concretized, src, edges, iterblocks)
        changed |= add_typedefs!(concretize, src, edges, typedefs, ())
    end
    debug && (println("after initial requirements discovery:"); print_with_code(stdout::IO, src, concretize))

    # find a loop region and check if any of the requirements discovered so far is involved
    # with it, and if require everything involved with the loop in order to properly
    # concretize the requirement — "require everything involved with the loop" means we will
    # care about "strongly connected components" rather than "cycles" of a directed graph, and
    # strongly connected components detection runs in linear time ``O(|V|+|E|)`` as opposed to
    # cycle dection (, whose worst time complexity is exponential with the number of vertices)
    # and thus the analysis here should terminate in reasonable time even with a fairly
    # complex control flow graph
    cfg = compute_basic_blocks(src.code)
    loops = filter!(>(1)∘length, strongly_connected_components(cfg)) # filter
    debug && @show loops

    critical_blocks = BitSet()
    for (i, block) in enumerate(cfg.blocks)
        if any(view(concretize, block.stmts))
            for loop in loops
                if i in loop
                    push!(critical_blocks, i)
                    for j in loop
                        push!(critical_blocks, j)
                    end
                end
            end
        end
    end
    for loop in loops
        if any(in(critical_blocks), loop)
            for j in loop
                push!(critical_blocks, j)
            end
            # push!(critical_blocks, minimum(loop) - 1)
        end
    end
    debug && @show critical_blocks

    norequire = BitSet()
    for (i, block) in enumerate(cfg.blocks)
        if i ∉ critical_blocks
            pushall!(norequire, rng(block))
        end
    end
    debug && @show norequire

    changed = true
    while changed
        changed = false

        # track SSA predecessors and control flows of the critical blocks
        changed |= add_ssa_preds!(concretize, src, edges, norequire)
        changed |= add_control_flow!(concretize, cfg, norequire)
    end
    debug && (println("after all requirements discovery:"); print_with_code(stdout::IO, src, concretize))
end

function JuliaInterpreter.step_expr!(interp::ConcreteInterpreter, frame::Frame, @nospecialize(node), istoplevel::Bool)
    @assert istoplevel "JET.ConcreteInterpreter can only work for top-level code"

    # TODO:
    # - support package analysis
    # - add report pass (report usage of undefined name, etc.)
    if ismoduleusage(node)
        for ex in to_simple_module_usages(node::Expr)
            # NOTE: usages of abstract global variables also work here, since they are supposed
            # to be actually evaluated into `interp.context` (as `AbstractGlobal`
            # object) at this point

            interp.eval_with_err_handling(interp.context, ex)
        end
        return frame.pc += 1
    end

    res = @invoke JuliaInterpreter.step_expr!(interp, frame, node, true::Bool)

    interp.config.analyze_from_definitions && collect_toplevel_signature!(interp, frame, node)

    return res
end

function collect_toplevel_signature!(interp::ConcreteInterpreter, frame::Frame, @nospecialize(node))
    if isexpr(node, :method, 3)
        sigs = node.args[2]
        atype_params, sparams, _ = @lookup(moduleof(frame), frame, sigs)::SimpleVector
        # t = atype_params[1]
        # if isdefined(t, :name)
        #     # XXX ignore constructor methods, just because it can lead to false positives ...
        #     t.name === CC._TYPE_NAME && return
        # end
        atype = Tuple{(atype_params::SimpleVector)...}
        push!(interp.res.toplevel_signatures, atype)
    end
end

ismoduleusage(@nospecialize(x)) = isexpr(x, (:import, :using, :export))

# assuming `ismoduleusage(x)` holds
function to_simple_module_usages(x::Expr)
    if length(x.args) != 1
        # using A, B, export a, b
        return Expr.(x.head, x.args)
    else
        arg = first(x.args)
        if isa(arg, Symbol)
            # export a
            return [x]
        else
            # import Pkg as P
            if arg.head === :as
                arg = first(arg.args)
            end
            if arg.head === :.
                # using A
                return [x]
            else
                # using A: sym1, sym2, ...
                @assert isexpr(arg, :(:))
                a, as... = arg.args
                return Expr.(x.head, Expr.(arg.head, Ref(a), as))::Vector{Expr}
            end
        end
    end
end

# adapted from https://github.com/JuliaDebug/JuliaInterpreter.jl/blob/2f5f80034bc287a60fe77c4e3b5a49a087e38f8b/src/interpret.jl#L188-L199
# works almost same as `JuliaInterpreter.evaluate_call_compiled!`, but also added special
# cases for JET analysis
function JuliaInterpreter.evaluate_call_recurse!(interp::ConcreteInterpreter, frame::Frame, call_expr::Expr; enter_generated::Bool=false)
    # @assert !enter_generated
    pc = frame.pc
    ret = bypass_builtins(interp, frame, call_expr, pc)
    isa(ret, Some{Any}) && return ret.value
    ret = @invokelatest maybe_evaluate_builtin(frame, call_expr, false)
    isa(ret, Some{Any}) && return ret.value
    args = collect_args(interp, frame, call_expr)
    f = first(args)
    if isinclude(f)
        return handle_include(interp, args)
    else
        popfirst!(args)  # now it's really just `args`
        # `_virtual_process!` iteratively interpret toplevel expressions but it doesn't hit toplevel
        # we may want to make `_virtual_process!` hit the toplevel on each interation rather than
        # using `invokelatest` here, but assuming concretized calls are supposed only to be
        # used for other toplevel definitions and as such not so computational heavy,
        # I'd like to go with this simplest way
        return @invokelatest f(args...)
    end
end

isinclude(@nospecialize f) = isa(f, Function) && nameof(f) === :include

function handle_include(interp, args)
    filename = interp.filename
    res = interp.res
    lnn = interp.lnn
    context = interp.context

    function handle_actual_method_error!(args)
        err = MethodError(args[1], args[2:end])
        report = ActualErrorWrapped(err, [], filename, lnn.line)
        push!(res.toplevel_error_reports, report)
    end

    nargs = length(args)
    if nargs == 2
        fname = args[2]
    elseif nargs == 3
        x = args[2]
        fname = args[3]
        if isa(x, Module)
            context = x
        elseif isa(x, Function)
            @warn "JET is unable to analyze `include(mapexpr::Function, filename::String)` call currently."
        else
            handle_actual_method_error!(args)
            return nothing
        end
    else
        handle_actual_method_error!(args)
        return nothing
    end
    if !isa(fname, AbstractString)
        handle_actual_method_error!(args)
        return nothing
    end

    include_file = normpath(dirname(filename), fname)
    # handle recursive `include`s
    if include_file in res.included_files
        report = RecursiveIncludeErrorReport(include_file, res.included_files, filename, lnn.line)
        push!(res.toplevel_error_reports, report)
        return nothing
    end

    function read_err_handler(@nospecialize(err), st)
        report = ActualErrorWrapped(err, st, filename, lnn.line)
        push!(res.toplevel_error_reports, report)
        return nothing
    end
    # `scrub_offset = 3` corresponds to `with_err_handling` -> `f`
    include_text = with_err_handling(read_err_handler, #= scrub_offset =# 2) do
        read(include_file, String)
    end
    isnothing(include_text) && return nothing # typically no file error

    _virtual_process!(include_text::String, include_file, interp.analyzer, interp.config, context, interp.res)

    # TODO: actually, here we need to try to get the lastly analyzed result of the `_virtual_process!` call above
    return nothing
end

struct RecursiveIncludeErrorReport <: ToplevelErrorReport
    duplicated_file::String
    files::Set{String}
    file::String
    line::Int
end

const JET_VIRTUALPROCESS_FILE = Symbol(@__FILE__)
const JULIAINTERPRETER_BUILTINS_FILE = let
    jlfile = pathof(JuliaInterpreter)::String
    Symbol(normpath(jlfile, "..", "builtins.jl"))
end

# handle errors from toplevel user code
function JuliaInterpreter.handle_err(interp::ConcreteInterpreter, frame, err)
    # catch stack trace
    bt = catch_backtrace()
    st = stacktrace(bt)

    # scrub the original stacktrace so that it only contains frames from user code
    i = 0
    for (j, frame) in enumerate(st)
        # if errors happen in `JuliaInterpreter.maybe_evaluate_builtin`, we just discard all
        # the stacktrace assuming they are enough self-explanatory (corresponding to the last logic below)
        if frame.file === JULIAINTERPRETER_BUILTINS_FILE && frame.func === :maybe_evaluate_builtin
            break # keep `i = 0`
        end

        # find an error frame that happened at `@invokelatest f(fargs...)` in the overload
        # `JuliaInterpreter.evaluate_call_recurse!(interp::ConcreteInterpreter, frame::Frame, call_expr::Expr; enter_generated::Bool=false)`
        if frame.file === JET_VIRTUALPROCESS_FILE && frame.func === :evaluate_call_recurse!
            i = j - 4 # offset: `evaluate_call_recurse` -> kwfunc (`evaluate_call_recurse`) -> `invokelatest` -> kwfunc (`invokelatest`)
            break
        end

        # other general errors may happen at `collect_args`, etc.
        # we don't show any stacktrace for those errors (by keeping the original `i = 0`)
        # since they are hopefully enough self-explanatory (XXX is it really so ?) and even
        # Julia's base error handler only shows something like `[1] top-level scope` in these cases
        continue
    end
    st = st[1:i]

    report = is_missing_concretization(err) ?
             MissingConcretization(err, st, interp.filename, interp.lnn.line) :
             ActualErrorWrapped(err, st, interp.filename, interp.lnn.line)
    push!(interp.res.toplevel_error_reports, report)

    return nothing # stop further interpretation
end

function with_err_handling(f, err_handler, scrub_offset)
    return try
        f()
    catch err
        bt = catch_backtrace()
        st = stacktrace(bt)

        # scrub the original stacktrace so that it only contains frames from user code
        i = findfirst(st) do frame
            frame.file === JET_VIRTUALPROCESS_FILE && frame.func === :with_err_handling
        end
        @assert i !== nothing
        st = st[1:(i - scrub_offset)]

        err_handler(err, st)
    end
end

let s = string(nameof(AbstractGlobal))
    global function is_missing_concretization(@nospecialize(err))
        io = IOBuffer()
        showerror(io, err)
        occursin(s, String(take!(io)))
    end
end

function collect_syntax_errors(s, filename)
    reports = SyntaxErrorReport[]
    index = line = 1
    while begin
            ex, nextindex = @static VERSION ≥ v"1.8.0-DEV.1370" ?
                _parse_string(s, filename, line, index, :statement) :
                _parse_string(s, filename, index, :statement)
            !isnothing(ex)
        end
        line += count(==('\n'), s[index:nextindex-1])
        report = isexpr(ex, :error) ? SyntaxErrorReport(string("syntax: ", first(ex.args)), filename, line) :
                 isexpr(ex, :incomplete) ? SyntaxErrorReport(first(ex.args), filename, line) :
                 nothing
        isnothing(report) || push!(reports, report)
        index = nextindex
    end
    return reports
end

struct SyntaxErrorReport <: ToplevelErrorReport
    err::ErrorException
    file::String
    line::Int
    SyntaxErrorReport(msg::AbstractString, file, line) = new(ErrorException(msg), file, line)
end


text:
using Documenter, Backtest

makedocs(
    sitename="Backtest Documentation",
    modules=[Backtest]
)

deploydocs(
    repo = "github.com/maxholloway/Backtest.jl.git",
)

ext: jl

max_stars_repo_name: maxholloway/backtest

max_issues_repo_name: maxholloway/backtest

max_forks_repo_name: maxholloway/backtest


text:
using .DifferentialEquations
const DE = DifferentialEquations

import .DifferentialEquations: controls

export trajectory,
       trajectories,
       controls,
       inputs,
       solution

struct SimulationSolution{TT, CT, IT}
    trajectory::TT  # trajectory pieces for each control cycle
    controls::CT  # control inputs for each control cycle
    inputs::IT  # nondeterministic inputs for each control cycle
end

trajectory(sol::SimulationSolution) = sol.trajectory
controls(sol::SimulationSolution) = sol.controls
inputs(sol::SimulationSolution) = sol.inputs

struct EnsembleSimulationSolution{TT, CT, IT}
    solutions::Vector{SimulationSolution{TT, CT, IT}}
end

# constructor from a bulk input
function EnsembleSimulationSolution(simulations, controls, inputs)
    n = length(simulations)  # number of pieces
    m = length(simulations[1])  # number of trajectories
    @assert n == length(controls) == length(inputs) "incompatible lengths"
    @assert all(m == length(piece) for piece in simulations)

    simulations_new = @inbounds [[simulations[i][j] for i in 1:n] for j in 1:m]
    controls_new = @inbounds [[controls[i][j] for i in 1:n] for j in 1:m]
    inputs_new = @inbounds [[(isassigned(inputs, i) ? inputs[i][j] : nothing)
        for i in 1:n] for j in 1:m]
    solutions = @inbounds [SimulationSolution(simulations_new[j],
        controls_new[j], inputs_new[j]) for j in 1:m]
    return EnsembleSimulationSolution(solutions)
end

Base.length(ess::EnsembleSimulationSolution) = length(ess.solutions)
solution(ess::EnsembleSimulationSolution, i) = ess.solutions[i]
trajectory(ess::EnsembleSimulationSolution, i) = trajectory(solution(ess, i))
trajectories(ess::EnsembleSimulationSolution) = trajectory.(ess.solutions)
controls(ess::EnsembleSimulationSolution, i) = controls(solution(ess, i))
controls(ess::EnsembleSimulationSolution) = controls.(ess.solutions)
inputs(ess::EnsembleSimulationSolution, i) = inputs(solution(ess, i))
inputs(ess::EnsembleSimulationSolution) = inputs.(ess.solutions)

# simulation of multiple trajectories for an ODE system and a time span
# currently we can't use this method from RA because the sampling should be made from outside the function
function _solve_ensemble(ivp, X0_samples, tspan;
                         trajectories_alg=DE.Tsit5(),
                         ensemble_alg=DE.EnsembleThreads(),
                         inplace=true,
                         kwargs...)
    if inplace
        field = ReachabilityAnalysis.inplace_field!(ivp)
    else
        field = ReachabilityAnalysis.outofplace_field(ivp)
    end

    _prob_func(prob, i, repeat) = remake(prob, u0=X0_samples[i])
    ensemble_prob = EnsembleProblem(ODEProblem(field, first(X0_samples), tspan),
                                    prob_func=_prob_func)
    return DE.solve(ensemble_prob, trajectories_alg, ensemble_alg;
                    trajectories=length(X0_samples))
end


text:
using FastGaussQuadrature, Compat
using Compat.Test, Compat.LinearAlgebra, Compat.Random, SpecialFunctions

if VERSION < v"0.7-"
    const floatmax = realmax
end

println("Chebyshev tests")
include("test_gausschebyshev.jl")
println("Legendre tests")
include("test_gausslegendre.jl")
println("Jacobi tests")
include("test_gaussjacobi.jl")
println("Radau tests")
include("test_gaussradau.jl")
println("Lobatto tests")
include("test_gausslobatto.jl")
println("Bessel roots tests")
include("test_besselroots.jl")
println("Laguerre tests")
include("test_gausslaguerre.jl")
println("Hermite tests")
include("test_gausshermite.jl")


text:
module FawziFangOperator;

include("HarmonicBases.jl")
include("SphericalQuadrature.jl")
using LinearAlgebra, DynamicPolynomials, FastGaussQuadrature, PolynomialBases, SpecialFunctions, .HarmonicBases, .SphericalQuadrature
using Convex, SCS

function risingfactorial(a,b)
    res = 1
    if b > 0
        res = a*risingfactorial(a+1,b-1)
    elseif b < 0
        throw(DomainError(b,"b must be positive"))
    end
    return res
end

Acoefficient(n,i,j) = 2.0^((n-2)/2 + 1) * binomial(i+n-2-1,i) * binomial(j+n-2-1,j)

bcoefficient(n,i,j,r) = sum([ sum([ (r+(n-2)/2)*(-1.0)^(p+r) * (binomial(i,p-m)*binomial(j,m)*(n-2+j)^m * risingfactorial(n-2+n-j,p-m)*risingfactorial(p,r) * risingfactorial((n-2)/2 + 1/2, p) )/(risingfactorial((n-2)/2 + 1/2, m)*risingfactorial((n-2)/2 + 1/2, p-m)*risingfactorial(n-2,r+p+1)) for m in 0:1:p]) for p in r:1:n])

function eigenvalues(n,m,k)
    e = Convex.Variable(m+1)
end

function evaluation(m,p,vars)
    desc = HarmonicBases.harmonicdecomposition(p,vars)
    n = length(vars)
    k = floor(Int,maxdegree(p)/2)
    if 2*k != maxdegree(p)
        throw(DomainError(p,"polynomial p must have even degree"))
    end
    if 2*k != mindegree(p)
        throw(DomainError(p,"polynomial p must be homogeneous"))
    end
    eigvals = eigenvalues(n,m,2*k)
    norm = sum(vars .* vars)
    return sum( eigvals[floor(Int,u[1]/2)+1]*norm^(floor(Int,u[1]/2))*u[2] for u in desc )
end

function inverseeval(m,p,vars)
    desc = HarmonicBases.harmonicdecomposition(p,vars)
    n = length(vars)
    k = floor(Int,maxdegree(p)/2)
    if 2*k != maxdegree(p)
        throw(DomainError(p,"polynomial p be of even degree"))
    end
    if 2*k != mindegree(p)
        throw(DomainError(p,"polynomial p must be homogeneous"))
    end
    eigvals = eigenvalues(n,m,k)
    norm = sum(vars .* vars)
    return sum( (1/eigvals[floor(Int,u[1]/2)+1])*norm^(floor(Int,u[1]/2))*u[2] for u in desc )
end

end


text:
using Plots
include(joinpath(@__DIR__, "quadrotor_broken_propeller.jl"))
include(joinpath(pwd(), "models/visualize.jl"))

x̄, ū = unpack(z̄, prob)
x, u = unpack(z, prob)

@show sum([ū[t][end] for t = 1:T-1])
t_nominal = [0.0, [sum([ū[i][end] for i = 1:t]) for t = 1:T-1]...]
@show sum([u[t][end] for t = 1:T-1])
t_dpo = [0.0, [sum([u[i][end] for i = 1:t]) for t = 1:T-1]...]

plot(t_nominal, hcat(x̄...)[1:3,:]', linetype=:steppost, width = 2.0)
plot(t_nominal[1:end-1], hcat(ū...)', linetype=:steppost, width = 2.0)

plot(t_dpo, hcat(x...)[1:3,:]', linetype=:steppost, width = 2.0)
plot(t_dpo[1:end-1], hcat(u...)', linetype=:steppost, width = 2.0)

vis = Visualizer()
render(vis)
visualize!(vis, model, x̄, Δt = ū[1][end])

# # Plots results
using Plots
# # Position trajectory
x_nom_pos = [X_nom[t][1] for t = 1:T]
y_nom_pos = [X_nom[t][2] for t = 1:T]
# pts = Plots.partialcircle(0,2π,100,r)
# cx,cy = Plots.unzip(pts)
# cx1 = [_cx + xc1 for _cx in cx]
# cy1 = [_cy + yc1 for _cy in cy]
# cx2 = [_cx + xc2 for _cx in cx]
# cy2 = [_cy + yc2 for _cy in cy]
# cx3 = [_cx + xc3 for _cx in cx]
# cy3 = [_cy + yc3 for _cy in cy]
# cx4 = [_cx + xc4 for _cx in cx]
# cy4 = [_cy + yc4 for _cy in cy]
# # cx5 = [_cx + xc5 for _cx in cx]
# # cy5 = [_cy + yc5 for _cy in cy]
#
# plt = plot(Shape(cx1,cy1),color=:red,label="",linecolor=:red)
# # plt = plot!(Shape(cx2,cy2),color=:red,label="",linecolor=:red)
# # plt = plot!(Shape(cx3,cy3),color=:red,label="",linecolor=:red)
# # plt = plot!(Shape(cx4,cy4),color=:red,label="",linecolor=:red)
# # # plt = plot(Shape(cx5,cy5),color=:red,label="",linecolor=:red)
#
# for i = 1:N
#     x_sample_pos = [X_sample[i][t][1] for t = 1:T]
#     y_sample_pos = [X_sample[i][t][2] for t = 1:T]
#     plt = plot!(x_sample_pos,y_sample_pos,aspect_ratio=:equal,
#         color=:cyan,label= i != 1 ? "" : "sample")
# end
# display(plt)
plt = plot()
plt = scatter!(x_nom_pos,y_nom_pos,aspect_ratio=:equal,xlabel="x",ylabel="y",width=4.0,label="TO",color=:purple,legend=:topleft)
x_sample_pos = [X_nom_sample[t][1] for t = 1:T]
y_sample_pos = [X_nom_sample[t][2] for t = 1:T]
plt = plot!(x_sample_pos,y_sample_pos,aspect_ratio=:equal,width=4.0,label="DPO",color=:orange,legend=:bottomright)
#
# savefig(plt,joinpath(@__DIR__,"results/quadrotor_trajectory.png"))
#
# Control
plt = plot(t_nominal[1:T-1],Array(hcat(U_nom...))',color=:purple,width=2.0,
    title="quad",xlabel="time (s)",ylabel="control",
    legend=:bottom,linetype=:steppost)
plt = plot!(t_sample[1:T-1],Array(hcat(U_nom_sample...))',color=:orange,
    width=2.0,linetype=:steppost)
savefig(plt,joinpath(@__DIR__,"results/quadrotor_control.png"))

# # Samples
X_sample
# State samples
plt1 = plot(title="Sample states",legend=:bottom,xlabel="time (s)");
for i = 1:N
    # t_sample = zeros(T)
    # for t = 2:T
    #     t_sample[t] = t_sample[t-1] + H_nom_sample[t-1]
    # end
    plt1 = plot!(hcat(X_sample[i]...)',label="");
end
plt1 = plot!(hcat(X_nom_sample...)',color=:red,width=2.0,
    label=["nominal" "" ""])
display(plt1)
# savefig(plt1,joinpath(@__DIR__,"results/quadrotor_sample_states.png"))
#
# # Control samples
# plt2 = plot(title="Sample controls",xlabel="time (s)",legend=:bottom);
# for i = 1:N
#     t_sample = zeros(T)
#     for t = 2:T
#         t_sample[t] = t_sample[t-1] + H_nom_sample[t-1]
#     end
#     plt2 = plot!(t_sample[1:end-1],hcat(U_sample[i]...)',label="",
#         linetype=:steppost);
# end
# plt2 = plot!(t_sample[1:end-1],hcat(U_nom_sample...)',color=:red,width=2.0,
#     label=["nominal" ""],linetype=:steppost)
# display(plt2)
# savefig(plt2,joinpath(@__DIR__,"results/quadrotor_sample_controls.png"))

# plot states
plt_x = plot(t_nom,hcat(X_nom...)[1:model.nx,:]',
	legend=:topright,color=:red,
    label="",width=2.0,xlabel="time (s)",
    title="Quadrotor",ylabel="state")
plt_x = plot!(t_sim_nom,hcat(z_tvlqr1...)[1:model.nx,:]',color=:black,label="",
    width=1.0)
plt_x = plot!(t_sim_nom,hcat(z_tvlqr2...)[1:model.nx,:]',color=:black,label="",
    width=1.0)
plt_x = plot!(t_sim_nom,hcat(z_tvlqr3...)[1:model.nx,:]',color=:black,label="",
    width=1.0)
plt_x = plot!(t_sim_nom,hcat(z_tvlqr4...)[1:model.nx,:]',color=:black,label="",
    width=1.0)

plt_x = plot(t_nom_sample,hcat(X_nom_sample...)[1:model.nx,:]',
	legend=:topright,color=:red,
    label="",width=2.0,xlabel="time (s)",
    title="Quadrotor",ylabel="state")
plt_x = plot!(t_sim_nom_sample,hcat(z_sample1...)[1:model.nx,:]',color=:black,label="",
    width=1.0)
plt_x = plot!(t_sim_nom_sample,hcat(z_sample2...)[1:model.nx,:]',color=:black,label="",
	width=1.0)
plt_x = plot!(t_sim_nom_sample,hcat(z_sample3...)[1:model.nx,:]',color=:black,label="",
    width=1.0)
plt_x = plot!(t_sim_nom_sample,hcat(z_sample4...)[1:model.nx,:]',color=:black,label="",
    width=1.0)

# plot COM
plot_traj = plot(hcat(X_nom...)[1,:],hcat(X_nom...)[2,:],
	legend=:topright,color=:red,
    label="",width=2.0,xlabel="y",ylabel="z",
    title="Quadrotor")
plot_traj = plot!(hcat(z_tvlqr1...)[1,:],hcat(z_tvlqr1...)[2,:],
	color=:black,
    label="",width=1.0)
plot_traj = plot!(hcat(z_tvlqr2...)[1,:],hcat(z_tvlqr2...)[2,:],
	color=:black,
	label="",width=1.0)
plot_traj = plot!(hcat(z_tvlqr3...)[1,:],hcat(z_tvlqr3...)[2,:],
	color=:black,
    label="",width=1.0)
plot_traj = plot!(hcat(z_tvlqr4...)[1,:],hcat(z_tvlqr4...)[2,:],
	color=:black,
    label="",width=1.0)

plot_traj = plot(hcat(X_nom_sample...)[1,:],hcat(X_nom_sample...)[2,:],
	legend=:topright,color=:red,
    label="",width=2.0,xlabel="y",ylabel="z",
    title="Quadrotor")
plot_traj = plot!(hcat(z_sample1...)[1,:],hcat(z_sample1...)[2,:],
	color=:black,
    label="",width=1.0)
plot_traj = plot!(hcat(z_sample2...)[1,:],hcat(z_sample2...)[2,:],
	color=:black,
    label="",width=1.0)
plot_traj = plot!(hcat(z_sample3...)[1,:],hcat(z_sample3...)[2,:],
	color=:black,
    label="",width=1.0)
plot_traj = plot!(hcat(z_sample4...)[1,:],hcat(z_sample4...)[2,:],
	color=:black,
    label="",width=1.0)

plot(t_nom[1:end-1],hcat(U_nom...)',
	linetype=:steppost,width=2.0)
plot!(t_sim_nom[1:end-1],hcat(u_tvlqr1...)',
	linetype=:steppost,color=:black,width=1.0)

plot(t_nom[1:end-1],hcat(U_nom_sample...)',
	linetype=:steppost,width=2.0)
plot!(t_sim_nom[1:end-1],hcat(u_sample1...)',
	linetype=:steppost,color=:black,width=1.0)

# state tracking
(Jx_tvlqr1 + Jx_tvlqr2 + Jx_tvlqr3 + Jx_tvlqr4)/4.0
(Jx_sample1 + Jx_sample2 + Jx_sample3 + Jx_sample4)/4.0

# control tracking
(Ju_tvlqr1 + Ju_tvlqr2 + Ju_tvlqr3 + Ju_tvlqr4)/4.0
(Ju_sample1 + Ju_sample2 + Ju_sample3 + Ju_sample4)/4.0

# objective value
(J_tvlqr1 + J_tvlqr2 + J_tvlqr3 + J_tvlqr4)/4.0
(J_sample1 + J_sample2 + J_sample3 + J_sample4)/4.0

using PGFPlots
const PGF = PGFPlots

# TO trajectory
p_u1_nom = PGF.Plots.Linear(t_nominal[1:end],hcat(ū..., ū[end])[1,:],
    mark="none",style="const plot, color=cyan, line width=2pt, solid")
p_u2_nom = PGF.Plots.Linear(t_nominal[1:end],hcat(ū..., ū[end])[2,:],
    mark="none",style="const plot, color=cyan, line width=2pt, solid")
p_u3_nom = PGF.Plots.Linear(t_nominal[1:end],hcat(ū..., ū[end])[3,:],
    mark="none",style="const plot, color=cyan, line width=2pt, solid")
p_u4_nom = PGF.Plots.Linear(t_nominal[1:end],hcat(ū..., ū[end])[4,:],
    mark="none",style="const plot, color=cyan, line width=2pt, solid")

p_u1_dpo = PGF.Plots.Linear(t_dpo[1:end],hcat(u..., u[end])[1,:],
    mark="none",style="const plot, color=orange, line width=2pt, solid")
p_u2_dpo = PGF.Plots.Linear(t_dpo[1:end],hcat(u..., u[end])[2,:],
    mark="none",style="const plot, color=orange, line width=2pt, solid")
p_u3_dpo = PGF.Plots.Linear(t_dpo[1:end],hcat(u..., u[end])[3,:],
    mark="none",style="const plot, color=orange, line width=2pt, solid")
p_u4_dpo = PGF.Plots.Linear(t_dpo[1:end],hcat(u..., u[end])[4,:],
    mark="none",style="const plot, color=orange, line width=2pt, solid")

a1 = Axis([p_u1_nom;p_u1_dpo
    ],
    xmin=0, ymin=0,
    hideAxis=false,
	ylabel="u1",
	xlabel="time",
	)
a2 = Axis([p_u2_nom;p_u2_dpo
    ],
    xmin=0, ymin=0,
    hideAxis=false,
	ylabel="u2",
	xlabel="time",
	)
a3 = Axis([p_u3_nom;p_u3_dpo
    ],
    xmin=0, ymin=0,
    hideAxis=false,
	ylabel="u3",
	xlabel="time",
	)
a4 = Axis([p_u4_nom;p_u4_dpo
    ],
    xmin=0, ymin=0,
    hideAxis=false,
	ylabel="u4",
	xlabel="time",
	)

# Save to tikz format
dir = joinpath(pwd(), "examples/direct_policy_optimization/figures")
PGF.save(joinpath(dir,"quad_prop_u1.tikz"), a1, include_preamble=false)
PGF.save(joinpath(dir,"quad_prop_u2.tikz"), a2, include_preamble=false)
PGF.save(joinpath(dir,"quad_prop_u3.tikz"), a3, include_preamble=false)
PGF.save(joinpath(dir,"quad_prop_u4.tikz"), a4, include_preamble=false)

function pad_trajectory(x, shift, T_shift)
	[[x[1] + shift for i = 1:T_shift]..., [_x + shift for _x in x]..., [x[end] + shift for i = 1:T_shift]...]
end

# # visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
open(vis)
shift = zero(x[1])
shift[1] = -1.5
shift[2] = -1.5
_shift = shift[1:3]
settransform!(vis["/Cameras/default"],
	compose(Translation(0.0, 0.0, 3.0),LinearMap(RotY(0.0*-pi/2.5))))
# visualize!(vis,model,[x̄[t] + shift for t = 1:T],Δt=ū[1][end])
# visualize!(vis,model,[x[t] + shift for t = 1:T],Δt=u[1][end])

pts_nom = collect(eachcol(hcat([z[1:3] + _shift for z in x̄]...)))
material_nom = LineBasicMaterial(color = colorant"cyan", linewidth = 10.0)
setobject!(vis["com_traj_nom"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_nom"],false)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift for z in z_lqr1]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_lqr1"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_lqr1"],false)
visualize!(vis,model,pad_trajectory(z_lqr1, shift, 100),Δt=dt_sim_nom)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift for z in z_lqr2]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_lqr2"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_lqr2"],false)
visualize!(vis,model,pad_trajectory(z_lqr2, shift, 100),Δt=dt_sim_nom)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift for z in z_lqr3]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_lqr3"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_lqr3"],false)
visualize!(vis,model,pad_trajectory(z_lqr3, shift, 100),Δt=dt_sim_nom)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift for z in z_lqr4]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_lqr4"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_lqr4"],false)
visualize!(vis,model,pad_trajectory(z_lqr4[27], shift, 100),Δt=dt_sim_nom)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift for z in x]...)))
material_nom = LineBasicMaterial(color = colorant"orange", linewidth = 10.0)
setobject!(vis["com_traj_dpo"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_dpo"],true)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift .+ 1.0e-3 for z in z_dpo1]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_dpo1"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_dpo1"],false)
visualize!(vis,model,pad_trajectory(z_dpo1, shift, 100),Δt=dt_sim_dpo)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift .+ 1.0e-3 for z in z_dpo2]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_dpo2"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_dpo2"],false)
visualize!(vis,model,pad_trajectory(z_dpo2, shift, 100),Δt=dt_sim_dpo)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift .+ 1.0e-3 for z in z_dpo3]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_dpo3"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_dpo3"],false)
visualize!(vis,model,pad_trajectory(z_dpo3, shift, 100),Δt=dt_sim_dpo)

pts_nom = collect(eachcol(hcat([z[1:3] + _shift .+ 1.0e-3 for z in z_dpo4]...)))
material_nom = LineBasicMaterial(color = colorant"gray", linewidth = 5.0)
setobject!(vis["com_traj_dpo4"], Object(PointCloud(pts_nom), material_nom, "Line"))
setvisible!(vis["com_traj_dpo4"],true)
visualize!(vis,model,pad_trajectory(z_dpo4, shift, 100),Δt=dt_sim_dpo)

# for (k,z) in enumerate([z_tvlqr1,z_tvlqr2,z_tvlqr3,z_tvlqr4])
# 	i = k
# 	q_to = z
# 	for t = 1:3:T_sim
# 		setobject!(vis["traj_to$t$i"], Sphere(Point3f0(0),
# 			convert(Float32,0.025)),
# 			MeshPhongMaterial(color=RGBA(128.0/255.0,128.0/255.0,128.0/255.0,1.0)))
# 		settransform!(vis["traj_to$t$i"], Translation((q_to[t][1],q_to[t][2],q_to[t][3])))
# 		setvisible!(vis["traj_to$t$i"],false)
# 	end
# end
# q_to_nom = X_nom
# for t = 1:T
# 	setobject!(vis["traj_to_nom$t"], Sphere(Point3f0(0),
# 		convert(Float32,0.05)),
# 		MeshPhongMaterial(color=RGBA(0.0,255.0/255.0,255.0/255.0,1.0)))
# 	settransform!(vis["traj_to_nom$t"], Translation((q_to_nom[t][1],q_to_nom[t][2],q_to_nom[t][3])))
# 	setvisible!(vis["traj_to_nom$t"],false)
# end
#
# for (k,z) in enumerate([z_sample1,z_sample2,z_sample3,z_sample4])
# 	i = k
# 	q_dpo = z
# 	for t = 1:3:T_sim
# 		setobject!(vis["traj_dpo$t$i"], Sphere(Point3f0(0),
# 			convert(Float32,0.025)),
# 			MeshPhongMaterial(color=RGBA(128.0/255.0,128.0/255.0,128.0/255.0,1.0)))
# 		settransform!(vis["traj_dpo$t$i"], Translation((q_dpo[t][1],q_dpo[t][2],q_dpo[t][3])))
# 		setvisible!(vis["traj_dpo$t$i"],true)
# 	end
# end
#
# q_dpo_nom = X_nom_sample
# for t = 1:T
# 	setobject!(vis["traj_dpo_nom$t"], Sphere(Point3f0(0),
# 		convert(Float32,0.05)),
# 		MeshPhongMaterial(color=RGBA(255.0/255.0,127.0/255.0,0.0,1.0)))
# 	settransform!(vis["traj_dpo_nom$t"], Translation((q_dpo_nom[t][1],q_dpo_nom[t][2],q_dpo_nom[t][3])))
# 	setvisible!(vis["traj_dpo_nom$t"],true)
# end
#
obj_path = joinpath(pwd(),
   "models/quadrotor/drone.obj")
 mtl_path = joinpath(pwd(),
   "models/quadrotor/drone.mtl")
ctm = ModifiedMeshFileObject(obj_path,mtl_path,scale=1.0)
setobject!(vis["drone2"],ctm)
settransform!(vis["drone2"], compose(Translation(x[10][1:3]),LinearMap(MRP(x[10][4:6]...)*RotX(pi/2.0))))

# open(vis)


text:
module testMuSigmaFunctions

using Base.Test
using DataFrames
using TimeData
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running mu-sigma function tests\n")

####################
## load test data ##
####################

include(joinpath(homedir(),
                 "research/julia/AssetMgmt/test/createTestData.jl"))
## scacapE6Names and mus20011202 and covMatr20011202

#############################
## test efficient frontier ##
#############################

## create mu-sigma model
mod = AssetMgmt.SampleMoments(mus20011202,
                              covMatr20011202,
                              scacapE6Names)

## get moments of efficient portfolios directly
nPfs = 50
effMus, effSigmas = AssetMgmt.getEffFrontier(mod, nPoints = nPfs)

## no asset may be better in both first and second moment
nAss = size(mod.mu, 1)
for ii=1:nPfs
    for jj=1:nAss
        betterMu = effMus[ii] < mod.mu[jj]
        betterVola = effSigmas[ii] > sqrt(mod.sigma[jj,jj])
        @test !(betterMu & betterVola)
    end
end

## get efficient portfolio weights
xEff = AssetMgmt.getEffPortfolios(mod, nPoints = nPfs)

## calculate moments for efficient portfolios
effMus2 = zeros(nPfs)
effSigmas2 = zeros(nPfs)
for ii=1:nPfs
    effMus2[ii] = AssetMgmt.getPMean(xEff[ii, :][:], mod.mu)
    effSigmas2[ii] =
        sqrt(AssetMgmt.getPVar(xEff[ii, :][:], mod.sigma))
end

## test equality of both approaches
@test_approx_eq_eps effMus effMus2 1e-14
@test_approx_eq_eps effSigmas effSigmas2 1e-14

## test that efficient portfolios are weights
@test_approx_eq_eps sum(xEff, 2) ones(Float64, nPfs) 1e-12

end


text:
using TSVD: tsvd

"""
    PCA(ncomps; options...)
    PCA(;ncomps=nothing, options...)

Use principal components analysis (PCA) to form a low-rank orthonormal basis of the input. Uses deterministic singular-value decomposition (SVD) to decompose data.

If `ncomps` is `nothing`, the basis will not be truncated (i.e. `ncomps` is equal to the number of frames). `ncomps` can be set to `:noise` or `:pratio` to automatically choose the number of components using the residual frame noise or principal ratio, respectively. For more information, see the extended help.

# References
1. [Soummer, Pueyo, and Larkin (2012)](https://ui.adsabs.harvard.edu/abs/2012ApJ...755L..28S) "Detection and Characterization of Exoplanets and Disks Using Projections on Karhunen-Loève Eigenimages"

# Extended help

## Optimizing `ncomps`

There are a few ways to optimize `ncomps` using the input data. Additional options for the optimization are listed below
1. `ncomps=:noise` - residual noise optimization
2. `ncomps=:pratio` - principal ratio optimization

### Residual noise optimization

This technique progressively increases `ncomps` at each step measuring the pixel-to-pixel noise (standard deviation) in the residual data. Iteration will stop when the noise is not improving beyond a threshold. This is suited for data with similar statistical characteristics, such as an annulus more so than a full-frame cube.
* `collapse=false` - if true, the temporal median of the residual data will be used for measuring the noise.
* `noise_error=1e-3` - the threshold for the minimal noise improvement looking back 2 iterations

### Principal ratio optimization

This technique chooses the number of components required to explain some ratio of the total variance in the data. This is known as the *principal ratio* or the *explained variance ratio*. The explained variance is measured by transforming the singular values of the SVD decomposition (`Λ = @. S^2 / (n - 1)`).
* `pratio=0.9` - the target principal ratio (between 0 and 1)
"""
struct PCA <: ADIAlgorithm
    ncomps
    opts
end
PCA(ncomps; options...) = PCA(ncomps, options)
PCA(; ncomps=nothing, options...) = PCA(ncomps, options)

function fit(alg::PCA, data::AbstractMatrix; ref=data, kwargs...)
    # get number of components (using dispatch for symbolic args)
    k = get_ncomps(alg.ncomps, ref; alg.opts...)
    # fit SVD to get principal subspace of reference
    decomp = svd(collect(ref))
    # Get the principal components (principal subspace) and weights
    P = decomp.Vt[1:k, :]
    weights = data * P'
    return LinearDesign(P, weights)
end

# get ncomps using given value or num frames, whichever is smaller
get_ncomps(n::Int, data; kwargs...) = min(n, size(data, 1))
get_ncomps(::Nothing, data; kwargs...) = size(data, 1)

# get ncomps using automatic methods
function get_ncomps(s::Symbol, data; kwargs...)
    if s === :noise
        noise_decay_ncomps(data; kwargs...)
    elseif s === :pratio
        pratio_ncomps(data; kwargs...)
    else
        error("Invalid `ncomps`. Did you mean :noise or :pratio?")
    end
end

function noise_decay_ncomps(data; collapse=false, noise_error=1e-3)
    if collapse
        μ = mean(data; dims=1)
        σ2 = var(data; dims=1, mean=μ)
        X = @. (data - μ) / σ2
    else
        X = data .- mean(data, dims=1)
    end
    P = svd(X).Vt
    tmpr = similar(data)
    τ1 = τ2 = 0
    @progress name="Optimizing ncomps using residual noise" for ncomp in axes(data, 1)
        Pv = @view P[1:ncomp, :]
        tmpr .= X * (I - Pv'Pv)
        # calculate noise (standard deviation) optionally collapsing
        noise = collapse ? std(median(tmpr, dims=1)) : std(tmpr)
        # test if we've reached the noise decay tolerance
        if ncomp > firstindex(data, 1) + 2
            px_noise_decay = τ2 - noise
            @debug noise_decay=px_noise_decay noise=noise
            if px_noise_decay < noise_error
                @info "noise threshold reached with $ncomp components"
                return ncomp
            end
        end
        # update recursion variables
        τ2, τ1 = τ1, noise
    end
    return lastindex(data, 1)
end

function pratio_ncomps(data; pratio=0.9)
    @debug "Choosing ncomps required to explain $(pratio*100)% of data's temporal variance"
    X = data .- mean(data, dims=1)
    Λ = svd!(X).S
    n = length(Λ)
    exp_var = @. Λ^2 / (n - 1)
    ratio_cumsum = cumsum(exp_var ./ sum(exp_var))
    k = last(searchsorted(ratio_cumsum, pratio))
    @info "$(pratio*100)% of variance explained with $k components"
    return k
end

"""
    TPCA(;ncomps=nothing, options...)
    TPCA(ncomps; options...)

Perform principal components analysis (PCA) using truncated SVD (TSVD; provided by TSVD.jl) instead of deterministic SVD. This is often faster than [`PCA`](@ref) but is non-deterministic, so the results may be different. See the [`PCA`](@ref) docstring for more information about the options.

# See Also
* [`PCA`](@ref), [`TSVD.tsvd`](https://github.com/JuliaLinearAlgebra/TSVD.jl)
"""
struct TPCA <: ADIAlgorithm
    ncomps
    opts
end
TPCA(ncomps; options...) = TPCA(ncomps, options)
TPCA(; ncomps=nothing, options...) = TPCA(ncomps, options)

function fit(alg::TPCA, data::AbstractMatrix; ref=data, kwargs...)
    # get number of components (using dispatch for symbolic args)
    # TODO using automatic methods is not valid for TPCA using different SVD methods
    k = get_ncomps(alg.ncomps, ref; alg.opts...)
    # fit SVD to get principal subspace of reference
    U, Σ, V = tsvd(collect(ref), k)
    # Get the principal components (principal subspace) and weights
    P = V'
    weights = data * P'
    return LinearDesign(P, weights)
end


text:
using CSV

df = CSV.File(joinpath(@__DIR__, "benchmarks.csv"))

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

passes = [:forward, :backward]

plt = data(df) * mapping(
    :machine => sorter("dense", "convolution", "recurrent"),
    passes .=> (t -> (t / 1000)) .=> "runtime (μs)",
    color = dims(1) => renamer(passes) => "pass",
    dodge = dims(1) => renamer(passes) => "pass",
    col = :device,
    row = :problemsize => renamer("small" => "small size", "medium" => "medium size")
)

fig = Figure(resolution=(600, 800))

ag = draw!(
    fig, plt * visual(BarPlot, fillto = 1, dodge_gap = 0.),
    axis=(xticklabelrotation=π/6, yscale=log10, yticks= LogTicks(LinearTicks(6)))
)
legend!(fig[0, 1:2], ag, orientation=:horizontal)

display(fig)


text:
@require TimeSeries begin

    """
    Convert a TimeArray into an Array{Timestamp}

    ```julia
    Array{Timestamp}(a::TimeArray)
    ```

    ### Arguments

    * `a::TimeArray`

    ### Returns

    * `::Array{Timestamp}` : an array of Timestamp from a TimeArray

    ### Examples

    ```julia
    using MarketData
    Apple = Array{Timestamp}(AAPL)
    ```
    """

function convert{T}(::Type{TimestampArray}, A::TimeSeries.TimeArray{T,1})
    stamps = Array(Timestamp, length(A))
    for i in 1:length(A)
        stamps[i] = Timestamp(A.timestamp[i], A.values[i])
    end
    stamps
end

# convert{T}(at::Array{Timestamp{T},1}, ta::TimeSeries.TimeArray{T,1}) = convert(at, ta)
#Timestamp{T,N}(ta::TimeSeries.TimeArray{T,N}) = convert{T}(Array{Timestamp}, A::TimeSeries.TimeArray{T})

end # of begin statement


text:
export MultiCutGenerator, AvgCutGenerator

abstract type AbstractCutGenerator end
abstract type AbstractOptimalityCutGenerator <: AbstractCutGenerator end

struct FeasibilityCutGenerator <: AbstractCutGenerator end
function gencut(::FeasibilityCutGenerator, sp::SOI.AbstractStochasticProgram, parent, pool::SOI.AbstractSolutionPool, to::TimerOutput, ztol)
    for tr in SOI.get(sp, SOI.OutTransitions(), parent)
        if SOI.hassolution(pool, tr)
            trsol = SOI.getsolution(pool, tr)
            if SOI.getstatus(trsol) == :Infeasible
                @timeit to SOI.FCUTS_KEY SOI.addcut!(sp, tr, FeasibilityCut(SOI.feasibility_cut(trsol)...))
            end
        end
    end
end
function applycut(::FeasibilityCutGenerator, sp::SOI.AbstractStochasticProgram, node)
    for tr in SOI.get(sp, SOI.OutTransitions(), node)
        SOI.applycuts!(sp, tr, FeasibilityCut)
    end
end

# No Optimality Cut Generator
struct NoOptimalityCutGenerator <: AbstractOptimalityCutGenerator
end
nθ(::NoOptimalityCutGenerator, proba) = 0
needallsolutions(::NoOptimalityCutGenerator) = false
function gencut(::NoOptimalityCutGenerator, sp::SOI.AbstractStochasticProgram, parent, pool::SOI.AbstractSolutionPool, to::TimerOutput, ztol)
end
function applycut(::NoOptimalityCutGenerator, sp::SOI.AbstractStochasticProgram, node)
end

# Multi Cut Generator
struct MultiCutGenerator <: AbstractOptimalityCutGenerator
end
nθ(::MultiCutGenerator, proba) = length(proba)
needallsolutions(::MultiCutGenerator) = false
function gencut(::MultiCutGenerator, sp::SOI.AbstractStochasticProgram, parent, pool::SOI.AbstractSolutionPool, to::TimerOutput, ztol)
    for tr in SOI.get(sp, SOI.OutTransitions(), parent)
        if SOI.hassolution(pool, tr)
            trsol = SOI.getsolution(pool, tr)
            status = SOI.getstatus(trsol)
            if status != :Unbounded
                a, β = SOI.optimality_cut(trsol)
                @assert status == :Optimal
                if tr.childT !== nothing
                    aT = tr.childT' * a
                else
                    aT = a
                end
                sol = SOI.getsolution(pool)
                x = SOI.getnodevalue(sol)
                θ = SOI.getbellmanvalue(sp, tr, sol)
                if SOI.getstatus(sol) == :Unbounded || _lt(θ, β - dot(aT, x), ztol)
                    @timeit to SOI.OCUTS_KEY SOI.addcut!(sp, tr, MultiOptimalityCut(a, β))
                end
            end
        end
    end
end
function applycut(::MultiCutGenerator, sp::SOI.AbstractStochasticProgram, node)
    for tr in SOI.get(sp, SOI.OutTransitions(), node)
        SOI.applycuts!(sp, tr, MultiOptimalityCut)
    end
end

# Average Cut Generator
struct AvgCutGenerator <: AbstractOptimalityCutGenerator
end
nθ(::AvgCutGenerator, proba) = 1
needallsolutions(::AvgCutGenerator) = true
function gencut(::AvgCutGenerator, sp::SOI.AbstractStochasticProgram, parent, pool::SOI.AbstractSolutionPool, to::TimerOutput, ztol)
    # We need all transitions to be solved, feasible and bounded to generate an averaged cut
    (!SOI.allfeasible(pool) || !SOI.allbounded(pool)) && return
    avga = zeros(SOI.get(sp, SOI.Dimension(), parent))
    avgβ = 0.
    for tr in SOI.get(sp, SOI.OutTransitions(), parent)
        SOI.hassolution(pool, tr) || error("Average Cut Generator needs a solution for each transition")
        trsol = SOI.getsolution(pool, tr)
        status = SOI.getstatus(trsol)
        @assert status != :Unbounded
        a, β = SOI.optimality_cut(trsol)
        @assert status == :Optimal
        if tr.childT !== nothing
            aT = tr.childT' * a
        else
            aT = a
        end
        proba = SOI.get(sp, SOI.Probability(), tr)
        avga += proba * aT
        avgβ += proba * β
    end
    sol = SOI.getsolution(pool)
    x = SOI.getnodevalue(sol)
    θ = SOI.getbellmanvalue(sp, parent, sol)
    if SOI.getstatus(sol) == :Unbounded || _lt(θ, avgβ - dot(avga, x), ztol)
        @timeit to SOI.OCUTS_KEY SOI.addcut!(sp, parent, AveragedOptimalityCut(avga, avgβ))
    end
end
function applycut(::AvgCutGenerator, sp::SOI.AbstractStochasticProgram, node)
    SOI.applycuts!(sp, node, AveragedOptimalityCut)
end


text:
function elements(rv::Vector, mu::FloatingPoint)
    r, v = rv[1:3], rv[4:6]
    rm = norm(r)
    vm = norm(v)
    h = cross(r, v)
    hm = norm(h)
    k = [0, 0, 1]
    n = cross(k, h)
    nm = norm(n)
    xi = vm^2/2 - mu/rm
    ec = ((vm^2 - mu/rm)*r - v*dot(r, v))/mu
    ecc = norm(ec)
    if ecc != 1
        sma = -mu/(2*xi)
        p = sma*(1-ecc^2)
    else
        p = hm^2/mu
        sma = p
    end
    inc = acos(h[3]/hm)
    node = acos(n[1]/nm)
    peri = acos(dot(n, ec)/(ecc*nm))
    ano = acos(dot(ec, r)/(ecc*rm))
    if n[2] < 0
        node = 2*pi - node
    end
    if ec[3] < 0
        peri = 2*pi - peri
    end
    if dot(r, v) < 0
        ano = 2*pi - ano
    end
    return [sma, ecc, inc, node, peri, ano]
end

function elements(rv::Matrix, mu::FloatingPoint)
    m, n = size(rv)
    if m != 6 && n != 6
        error("'rv' must be a 6xN or Nx6 matrix.")
    elseif m == 6
        rv = rv'
        m, n = n, m
    end
    ele = zeros(m,n)
    for i = 1:m
        ele[i,:] = elements(vec(rv[i,:]), mu)
    end
    return ele
end

function elements(s::State)
    return elements(s.rv, planets[s.body]["mu"])
end

function elements(s::State, deg::Bool)
    ele = elements(s)
    if deg
        ele[3:end] = ele[3:end]*180/pi
        return ele
    else
        return ele
    end
end

function cartesian(ele::Vector, mu::FloatingPoint)
    sma, ecc, inc, lan, per, ano = ele
    u = per + ano
    if ecc == 1
        p = sma
    else
        p = sma*(1 - ecc^2)
    end
    r = p/(1 + ecc*cos(ano))
    x = r*(cos(lan)*cos(u) - sin(lan)*cos(inc)*sin(u))
    y = r*(sin(lan)*cos(u) + cos(lan)*cos(inc)*sin(u))
    z = r*sin(inc)*sin(u)
    vr = sqrt(mu/p)*ecc*sin(ano)
    vf = sqrt(mu*p)/r
    vx = ((vr*(cos(lan)*cos(u) - sin(lan)*cos(inc)*sin(u))
        - vf*(cos(lan)*sin(u) + sin(lan)*cos(u)*cos(inc))))
    vy = ((vr*(sin(lan)*cos(u) + cos(lan)*cos(inc)*sin(u))
        - vf*(sin(lan)*sin(u) - cos(lan)*cos(u)*cos(inc))))
    vz = vr*sin(inc)*sin(u) + vf*cos(u)*sin(inc)
    return [x, y, z, vx, vy, vz]
end

function cartesian(ele::Matrix, mu::FloatingPoint)
    m, n = size(ele)
    if m != 6 && n != 6
        error("'ele' must be a 6xN or Nx6 matrix.")
    elseif m == 6
        ele = ele'
        m, n = n, m
    end
    rv = zeros(m,n)
    for i = 1:m
        rv[i,:] = cartesian(vec(ele[i,:]), mu)
    end
    return rv
end


text:
module Debug
import Base.*
export @debug, debug_hook, debug_eval, Scope, getdefs, code_debug

const doublecolon = symbol("::")

quot(ex) = expr(:quote, ex)

is_expr(ex, head::Symbol) = (isa(ex, Expr) && (ex.head == head))
is_expr(ex, head::Symbol, n::Int) = is_expr(ex, head) && length(ex.args) == n

is_symbol(ex::Symbol)      = true
is_symbol(ex::SymbolNode ) = true
is_symbol(ex)              = false

get_symbol(ex::Symbol)     = ex
get_symbol(ex::SymbolNode) = ex.name

is_linenumber(ex::LineNumberNode) = true
is_linenumber(ex::Expr)           = is(ex.head, :line)
is_linenumber(ex)                 = false

get_linenumber(ex::Expr)           = ex.args[1]
get_linenumber(ex::LineNumberNode) = ex.line

# ---- macro_expand -----------------------------------------------------------

# expand() seems to do a number of things, including macro expansion.
# macro_expand() should confine it to macros.
function macro_expand(ex::Expr)
    if ex.head === :macrocall;                 return expand(ex); end
    if ex.head === :quote || ex.head === :top; return ex;         end
    expr(ex.head, {macro_expand(arg) for arg in ex.args})
end
macro_expand(ex) = ex

# ---- getdefs: gather defined symbols by scope -------------------------------

type DefinedSyms
    scopes::ObjectIdDict
    syms::Set{Symbol}
end

enter(c::DefinedSyms, ex) = enter(c.scopes, ex)
function enter(scopes::ObjectIdDict, ex)
    scopes[ex] = syms = Set{Symbol}()
    DefinedSyms(scopes, syms)
end

const comprehensions = [:comprehension, symbol("cell-comprehension")]

getdefs(c::DefinedSyms, exs::Vector) = (for ex in exs; getdefs(c, ex); end)
function getdefs(c::DefinedSyms, ex::Expr)
    head = ex.head
    args = ex.args

    if contains([:line, :quote, :top, :type], head); return; end
    if head === :macrocall; return; end  # don't go in there for now
    if head === :let
        body = args[1]
        c_outer, c_inner = c, enter(c, body)
        getdefs(c_inner, body)
        for arg in args[2:end]
            if is_expr(arg, :(=))
                lhs, rhs = arg.args[1], arg.args[2]
                getdefs(c_outer, rhs)
            else
                lhs = arg
            end
            getdefs_lhs(c_inner, lhs)
            c.scopes[lhs] = c_inner.syms
        end
        return
    end
    if contains(comprehensions, head)
        c = enter(c, args[1])
        getdefs(c, args)
        for arg in args[2:end]
            c.scopes[arg] = c.syms
        end
    end
    if head === :try
        try_body, catch_arg, catch_body = args[1], args[2], args[3]

        getdefs(enter(c, try_body), try_body)

        c = enter(c, catch_body)
        getdefs_lhs(c, catch_arg)
        getdefs(c, catch_body)
        c.scopes[catch_arg] = c.syms
        return
    end

    if contains([:function, :for, :while, :(->)], head)
        c = enter(c, ex)
    end
    if head === :(=) || head === :(->)
        getdefs_lhs(c, args[1])
        getdefs(c, args[2])
    elseif contains([:local, :global], head)
        getdefs_lhs(c, args)
    elseif head === :function
        getdefs_lhs(c, args[1])
        getdefs(c, args[2])
    else
        getdefs(c, args)
    end
end
getdefs(c::DefinedSyms, ex) = nothing

getdefs_lhs(c::DefinedSyms, exs::Vector) = (for e in exs; getdefs_lhs(c,e);end)
function getdefs_lhs(c::DefinedSyms, ex::Expr)
    head = ex.head
    args = ex.args
    nargs = length(args)

    if head === doublecolon && nargs == 2
        getdefs_lhs(c, args[1])
        getdefs(c, args[2])
    elseif contains([:call, :tuple], head)
        getdefs_lhs(c, args)
    elseif head === :ref
        getdefs(c, args)
    else
        error("getdefs_lhs: don't know how to handle head=$head")
    end
end
getdefs_lhs(c::DefinedSyms, ex::Symbol)     = (add(c.syms, ex); nothing)
getdefs_lhs(c::DefinedSyms, ex::SymbolNode) = getdefs(c, ex.name)
getdefs_lhs(c::DefinedSyms, ex)             = nothing

function getdefs(ex::Expr)
    scopes = ObjectIdDict()
    getdefs(enter(scopes, ex), ex)
    scopes
end

# ---- Scope: runtime symbol table with getters and setters -------------------

abstract Scope

type NoScope <: Scope; end
has(scope::NoScope, sym::Symbol) = false

type LocalScope <: Scope
    parent::Scope
    syms::Dict
end

function has(scope::LocalScope, sym::Symbol)
    has(scope.syms, sym) || has(scope.parent, sym)
end
function get_entry(scope::LocalScope, sym::Symbol)
    has(scope.syms, sym) ? scope.syms[sym] : get_entry(scope.parent, sym)
end

get_getter(scope::LocalScope, sym::Symbol) = get_entry(scope, sym)[1]
get_setter(scope::LocalScope, sym::Symbol) = get_entry(scope, sym)[2]
ref(   scope::LocalScope,     sym::Symbol) = get_getter(scope, sym)()
assign(scope::LocalScope, x,  sym::Symbol) = get_setter(scope, sym)(x)

function code_getset(sym::Symbol)
    val = gensym(string(sym))
    :( ()->($sym), ($val)->(($sym)=($val)) )
end
function code_scope(scope::Symbol, parent, syms)
    pairs = {expr(:(=>), quot(sym), code_getset(sym)) for sym in syms}
    :(local ($scope) = ($quot(LocalScope))(($parent), ($expr(:cell1d, pairs))))
end

# ---- @debug: instrument an AST with debug code ------------------------------

type DbgShared
    line::Integer
    file
    scopes::ObjectIdDict
    DbgShared(scopes::ObjectIdDict) = new(-1, nothing, scopes)
end

type CodeDebug
    shared::DbgShared
    scope_sym
    syms::Set{Symbol}
end
CodeDebug(shared::DbgShared, scope) = CodeDebug(shared, scope, Set{Symbol}())

function enter(c::CodeDebug, ex)
    CodeDebug(c.shared, c.scope_sym, union(c.syms, c.shared.scopes[ex]))
end

code_debug(c::CodeDebug, exs::Vector) = {code_debug(c, ex) for ex in exs}
function code_debug(c::CodeDebug, ex::Expr)
    if contains([:line, :quote, :top, :type], ex.head); return ex; end
    if ex.head === :macrocall; return ex; end  # don't go in there for now

    scopes = c.shared.scopes
    if has(scopes, ex); c = enter(c, ex); end

    if ex.head === :block || ex.head === :body
        args = {}
        if !isempty(c.syms)
            # emit Scope object
            scope_sym = gensym("scope")
            push(args, code_scope(scope_sym, c.scope_sym, c.syms))
            del_all(c.syms)
            c.scope_sym = scope_sym
        end
        for arg in ex.args
            if is_linenumber(arg)
                c.shared.line = get_linenumber(arg)
                if is_expr(arg, :line, 2)
                    c.shared.file = arg.args[2]
                end
                push(args, arg)
                continue
            end
            line, file, scope = c.shared.line, c.shared.file, c.scope_sym
            push(args, :( debug_hook(($quot(line)), ($quot(file)), ($scope)) ))
            push(args, code_debug(c, arg))
        end
        return expr(ex.head, args)
    end

    return expr(ex.head, code_debug(c, ex.args))
end
code_debug(c::CodeDebug, ex) = ex

function code_debug(ex)
#    ex = macro_expand(ex)
    println("ex =\n", ex)
    code_debug(CodeDebug(DbgShared(getdefs(ex)), quot(NoScope())), ex)
end

debug_hook(line::Int, scope::Scope) = nothing

macro debug(ex)
    globalvar = esc(gensym("globalvar"))
    quote
        ($globalvar) = false
        try
            global ($globalvar)
            ($globalvar) = true
        end
        if !($globalvar)
            error("@debug: must be applied in global scope!")
        end
        ($esc(code_debug(ex)))
    end
end

# ---- debug_eval: eval an expression inside a Scope (by substitution) --------

const updating_ops = {
 :+= => :+,   :-= => :-,  :*= => :*,  :/= => :/,  ://= => ://, :.//= => :.//,
:.*= => :.*, :./= => :./, :\= => :\, :.\= => :.\,  :^= => :^,   :.^= => :.^,
 :%= => :%,   :|= => :|,  :&= => :&,  :$= => :$,  :<<= => :<<,  :>>= => :>>,
 :>>>= => :>>>}

function resym(s::LocalScope, ex::Expr)
    head, args = ex.head, ex.args
    if head === :(=)
        if is_symbol(args[1])
            lhs = get_symbol(args[1])
            if !has(s, lhs) error("Cannot assign to ($lhs): not in scope.") end
            setter = get_setter(s, lhs)
            rhs = resym(s, args[2])
            expr(:call, quot(setter), rhs)
        elseif is_expr(args[1], :tuple)
            tup = esc(gensym("tuple"))
            resym(s, expr(:block, :( ($tup) = ($args[2])                ),
                                 {:( ($dest) = ($tupleref)(($tup),($k)) ) for
                                     (dest,k)=enumerate(args[1].args)}...))
        else
            expr(head, {resym(s, arg) for arg in args})
        end
    elseif has(updating_ops, head) # Translate updating ops, e g x+=1 ==> x=x+1
        op = updating_ops[head]
        resym(s, :( ($args[1]) = ($op)(($args[1]), ($args[2])) ))
    elseif head === :quote || head === :top; ex
    elseif head === :escape; ex.args[1]  # bypasses substitution
    else                     expr(head, {resym(s, arg) for arg in args})
    end
end
resym(s::LocalScope, node::SymbolNode) = resym(s, node.name)
function resym(s::LocalScope, ex::Symbol)
    has(s, ex) ? :(($quot(get_getter(s, ex)))()) : ex
end
resym(s::Scope, ex) = ex

debug_eval(scope::Scope, ex) = eval(resym(scope, ex))
# function debug_eval(scope::Scope, ex)
#     ex = resym(scope, ex)
#     println("ex = ", ex)
#     eval(ex)
# end

end  # module


text:
module Accounts

export EMAIL, WORK # Domain
export create # API
export add_to_file, read_from_file # Infrastructure

include("Domain.jl"); using .Domain
include("API.jl"); using .API
include("Infrastructure.jl"); using .Infrastructure

end