text:
__precompile__()
module ASTInterpreter

export enter, Environment, @enter
import Base: LineEdit, REPL

using AbstractTrees
using JuliaParser
using JuliaParser: Lexer, Tokens
using Base.Meta
import JuliaParser.Tokens: SourceNode, SourceRange, SourceExpr
using JuliaParser.Tokens: √
import JuliaParser.Diagnostics: diag, AbstractDiagnostic, display_diagnostic

import AbstractTrees: children, printnode

pretty_repr(x) = (string(x), :default)
pretty_repr(x::LineNumberNode) = (string("line: ", x.line), :red)
function pretty_repr(x::Expr)
    if x.head == :body
        ("Body", :blue)
    elseif x.head == :return
        ("Return", :green)
    elseif x.head == :call
        ("Call", :yellow)
    else
        return (repr(x), :default)
    end
end

function printnode(io::IO, x::Union{Expr, LineNumberNode}, color = :default)
    x = pretty_repr(x)
    if color == :default
        color = x[2]
    end
    if color == :default
        print(io, x[1])
    else
        print_with_color(color, io, x[1])
    end
end

# EvaluationTree

immutable Environment
    locals::Vector{Nullable{Any}}
    ssavalues::Vector{Any}
    sparams::Vector{Any}
    # A vector from names to the slotnumber of that name
    # for which a reference was last encountered.
    last_reference::Dict{Symbol, Int}
end
Environment() = Environment(Vector{Nullable{Any}}(), Any[], Any[], Dict{Symbol, Int}())

Base.copy(e::Environment) = Environment(copy(e.locals), copy(e.ssavalues), copy(e.sparams), copy(e.last_reference))

type Interpreter
    stack::Vector{Any}
    env::Environment
    linfo::LambdaInfo
    evaluating_staged::Bool
    ast::Any
    it::Any
    cur_state::Any
    next_expr::Any
    shadowtree::Any
    code::AbstractString
    loctree::Any
    retval::Any
    exception_frames::Vector{Int}
    did_wrappercall::Bool
    prehook::AbstractString
end

function make_shadowtree(tree)
    resulttree = copy(tree)
    annotations = AbstractTrees.make_annotations(resulttree, resulttree, false) do node, parent, parent_ev
        if parent_ev
            return parent_ev
        end
        unevaluated = isa(node, Expr) || isa(node, GlobalRef) || isa(node, Symbol) || isa(node,Slot) ||
            isa(node,SSAValue) || isa(node, GotoNode) || isa(node, QuoteNode)
        if isa(node, Expr) && (node.head == :meta || node.head == :boundscheck ||
            node.head == :inbounds || node.head == :line)
            unevaluated = false
        end
        if (isa(node, SSAValue) || isa(node, Symbol) || isa(node, Slot) || isa(node,GlobalRef)) && isexpr(parent,:(=)) && parent.args[1] == node
            unevaluated = false
        end
        if isexpr(parent, :static_typeof)
            unevaluated = false
        end
        !unevaluated
    end
    shadowtree = AbstractTrees.ShadowTree(Tree(resulttree), Tree(annotations))
    it = indenumerate(PostOrderDFS(resulttree))
    shadowtree, it
end

using AbstractTrees: getnode

function enter(linfo, tree::Expr, env, stack = Any[];
        loctree = nothing, code = "", evaluating_staged = false)
    shadowtree, it = make_shadowtree(tree)

    interp = Interpreter(stack, env, linfo, evaluating_staged, tree, it,
        start(it), nothing, shadowtree, code, loctree, nothing, Vector{Int}(),
        false, "")
    push!(stack, interp)
    ind, node = next_expr!(interp)

    while getnode(interp.shadowtree,ind).shadow.x.val
        ind, node = next_expr!(interp)
    end

    interp
end
function enter(meth::Union{Method, TypeMapEntry}, env::Environment, stack = Any[]; kwargs...)
    isa(meth, TypeMapEntry) && (meth = meth.func)
    linfo = meth.lambda_template
    code = Base.uncompressed_ast(linfo)
    tree = Expr(:body); tree.args = code
    enter(linfo, tree, env, stack; kwargs...)
end
function enter(linfo::LambdaInfo, env::Environment, stack = Any[]; kwargs...)
    if linfo.inferred
        return enter(linfo.def, env, stack; kwargs...)
    end
    code = Base.uncompressed_ast(linfo)
    tree = Expr(:body); tree.args = code
    enter(linfo, tree, env, stack; kwargs...)
end
enter(f::Function, env) = enter(first(methods(f)), env)

function print_shadowtree(shadowtree, highlight = nothing, inds = nothing)
    from = nothing
    to = nothing
    if inds !== nothing
        from = first(inds)
        to = last(inds)
    end
    AbstractTrees.print_tree(STDOUT, shadowtree, withinds = true, from = from, to = to) do io, annotatednode, inds
        node, evaluated = annotatednode.tree.x, annotatednode.shadow.x.val
        str, _ = pretty_repr(node)
        print_with_color((inds == highlight) ? :yellow : (evaluated ? :green : :red), io, str)
    end
end

function node_repr(x)
    isa(x,GlobalRef) ? string(x.mod,'.',x.name) : string(x)
end

abstract ReplacementLoc

immutable SimpleReplacementLoc <: ReplacementLoc
    replacing::SourceRange
    before::AbstractString
    after::AbstractString
end
Tokens.normalize(x::ReplacementLoc) = x.replacing
Tokens.merge(x::ReplacementLoc,y::ReplacementLoc) =
    Tokens.merge(Tokens.normalize(x),Tokens.normalize(y))
Tokens.merge(x,y::ReplacementLoc) = Tokens.merge(x,Tokens.normalize(y))
Tokens.merge(x::ReplacementLoc,y::Void) = Tokens.merge(Tokens.normalize(x),y)
Tokens.merge(x::ReplacementLoc,y::Lexer.Token) = Tokens.merge(Tokens.normalize(x),y)
Tokens.merge(x::ReplacementLoc,y) = Tokens.merge(Tokens.normalize(x),y)

# SequencingReplacementLoc
type SRLoc <: ReplacementLoc
    replacing::SourceRange
    sequence::Vector{Any}
    lastidx::Int
end

function sequence!(s::SRLoc, ind)
    sidx = findnext(x->x==0,s.sequence,s.lastidx+1)
    s.sequence[sidx] = ind
    s.lastidx = sidx
    s
end

immutable Suppressed{T}
    item::T
end
Base.show(io::IO, x::Suppressed) = print(io, "<suppressed ", x.item, '>')

immutable Coloring
    x
    color::Symbol
end
Base.print(io::IO, c::Coloring) = print_with_color(c.color, string(c.x))
Base.show_unquoted(io::IO, c::Coloring, x, y) =
    print_with_color(c.color, sprint(Base.show_unquoted, c.x, x, y))

function annotate_highlights!(x, highlights)
    wrapcolor = nothing
    for highlight in highlights
        if isempty(highlight[1])
            wrapcolor = highlight[2]
        else
            Tree(x)[highlight[1]] = Coloring(Tree(x)[highlight[1]], highlight[2])
        end
    end
    wrapcolor == nothing ? x : Coloring(x, wrapcolor)
end

is_loc_meta(expr, kind) = isexpr(expr, :meta) && length(expr.args) >= 1 && expr.args[1] === kind

function determine_line_and_file(interp, highlight)
    file = get_linfo(interp).def.file
    line = get_linfo(interp).def.line
    foundline = false
    extra_locs = Any[]
    # Find a line number node previous to this expression
    if highlight !== nothing && !isempty(highlight)
        exprtree = interp.shadowtree.tree.x
        i = highlight[1]
        while i >= 1
            expr = exprtree.args[i]
            if !foundline && isa(expr, LineNumberNode)
                line = expr.line
                foundline = true
            elseif !foundline && isexpr(expr, :line)
                line = expr.args[1]
                foundline = true
            elseif foundline && is_loc_meta(expr, :push_loc)
                file = expr.args[2]
                extra_locs = determine_line_and_file(interp, [i-1])
                break
            elseif is_loc_meta(expr, :pop_loc)
                npops = 1
                while npops >= 1
                    i -= 1
                    expr = exprtree.args[i]
                    is_loc_meta(expr, :pop_loc) && (npops += 1)
                    is_loc_meta(expr, :push_loc) && (npops -= 1)
                end
            end
            i -= 1
        end
    end
    [extra_locs; (file, line)]
end

print_sourcecode(linfo::LambdaInfo, code, line) =
    print_sourcecode(code, line, linfo.def.file, linfo.def.line; file = SourceFile(code))
function print_sourcecode(code, line, deffile, defline; file = SourceFile(code))
    startoffset, stopoffset = compute_source_offsets(code, file.offsets[line],
        max(1,file == deffile ? defline : line-1), line+3; file=file)

    if startoffset == -1
        print_with_color(:bold, STDOUT, "Line out of file range (bad debug info?)")
        return
    end

    # Compute necessary data for line numbering
    startline = compute_line(file, startoffset)
    stopline = compute_line(file, stopoffset)
    current_line = line
    stoplinelength = length(string(stopline))

    code = split(code[(startoffset:stopoffset)+1],'\n')
    lineno = startline

    if !isempty(code) && isempty(code[end])
        pop!(code)
    end

    print_with_color(:bold, STDOUT,
        string("In ", deffile,":",defline, "\n"))

    for textline in code
        print_with_color(lineno == current_line ? :yellow : :bold,
            string(lineno, " "^(stoplinelength-length(lineno)+1)))
        println(textline)
        lineno += 1
    end
    println()
end

"""
Determine the offsets in the source code to print, based on the offset of the
currently highlighted part of the code, and the start and stop line of the
entire function.
"""
function compute_source_offsets(code, offset, startline, stopline; file = SourceFile(interp.code))
    offsetline = compute_line(file, offset)
    if offsetline - 3 > length(file.offsets) || startline > length(file.offsets)
        return -1, -1
    end
    startoffset = max(file.offsets[max(offsetline-3,1)], file.offsets[startline])
    stopoffset = endof(code)-1
    if offsetline + 3 < endof(file.offsets)
        stopoffset = min(stopoffset, file.offsets[offsetline + 3]-1)
    end
    if stopline + 1 < endof(file.offsets)
        stopoffset = min(stopoffset, file.offsets[stopline + 1]-1)
    end
    startoffset, stopoffset
end

global fancy_mode = false

print_status(state, _::Void, args...) = nothing
function print_status(state, interp::Interpreter, highlight = idx_stack(interp); fancy = fancy_mode)
    if !fancy && !isempty(interp.code)
        fls = determine_line_and_file(interp, highlight)
        for (file,line) in fls[end:-1:2]
            # See if we can get the code for this
            dfile = Base.find_source_file(string(file))
            file = dfile == nothing ? file : dfile
            contents = readfileorhist(file)
            file = "macro expansion from $file"
            if contents == nothing
                print_with_color(:bold, STDOUT,
                    string("In ", file,":",line, "\n"))
                continue
            end
            print_sourcecode(contents, line, file, max(1,line-3))
        end
        print_sourcecode(interp.linfo, interp.code, fls[1][2])
        ex = interp.shadowtree[highlight].tree.x
        # print slots with their names
        wrap = !isa(ex,Expr)
        ex = wrap ? Expr(:block, ex) : copy(ex)
        treemap!(function (x)
                     if isa(x, Slot)
                         name = interp.linfo.slotnames[x.id]
                         return sym_visible(name) ? name : x
                     elseif isa(x, AbstractArray)
                         return length(x) <= 10 ? x : Suppressed(summary(x))
                     elseif isa(x, LambdaInfo)
                         return Suppressed("generated thunk")
                     else
                         return x
                     end
                 end,
                 PreOrderDFS(ex,node->isa(node, Expr)))
        nargs = length(ex.args)
        # Suppress arguments that are too long
        for (i, arg) in enumerate(ex.args)
            nbytes = length(repr(arg))
            if nbytes > max(40, div(200,nargs))
                ex.args[i] = Suppressed("$nbytes bytes of output")
            end
        end
        if wrap; ex = ex.args[1]; end
        println("About to run: ", ex)
    elseif interp.loctree === nothing
        print_shadowtree(interp.shadowtree, highlight)
    else
        union = Lexer.normalize(reduce(⤄,PostOrderDFS(interp.loctree)))
        file = SourceFile(interp.code)
        # Compute the start and stop line of the current function
        startline = compute_line(file, union.offset)
        stopline = compute_line(file, union.offset+union.length)
        if highlight != nothing
            x = AbstractTrees.getindexhighest(Tree(interp.loctree),highlight)
            highlighstart = x[1]
            while Lexer.normalize(Tree(interp.loctree)[highlighstart]) == SourceRange()
                highlighstart = highlighstart[1:(end-1)]
                !isempty(highlighstart) || break
            end
            locnode = Tree(interp.loctree)[highlighstart]
            offset = Lexer.normalize(locnode).offset

            if offset != -1 % UInt32
                startoffset, stopoffset =
                    compute_source_offsets(interp.code, offset, startline,
                        stopline, file = file)

                ###### New Algorithm
                idxstartwiwth(idxs, start) = length(idxs) >= length(start) && idxs[1:length(start)] == start
                ishighlight(x) = idxstartwiwth(x,highlight)

                # Figure out all indices we want to print
                indices = Set{Vector{Int}}()
                srlocs = Set{SRLoc}()
                locs = Vector{Any}()
                highlighting = Vector{Any}()
                for (ind, isevaluated) in indenumerate(PostOrderDFS(interp.shadowtree.shadow))
                    if (ishighlight(ind) && length(ind)==length(highlight)) ||
                       (!ishighlight(ind) && isevaluated.val)
                        orgind = copy(ind)
                        while !isempty(ind) && try
                                Lexer.normalize(Tree(interp.loctree)[ind])
                            catch e
                                isa(e, BoundsError) || rethrow(e)
                                SourceRange()
                            end == SourceRange()
                            ind = ind[1:end-1]
                        end
                        locnode = Tree(interp.loctree)[ind]
                        if Lexer.normalize(locnode) == SourceRange() || ind == []
                            continue
                        end
                        off = Lexer.normalize(locnode).offset
                        if off > stopoffset || off < startoffset
                            continue
                        end
                        push!(highlighting,(ind, orgind, ishighlight(orgind) ? :yellow : :green))
                        if isa(locnode.loc, SRLoc)
                            if !ishighlight(orgind)
                                continue
                            end
                            r = locnode.loc
                            r in srlocs && continue
                            push!(srlocs, r)
                        else
                            r = ind
                            r in indices && continue
                            push!(indices, r)
                        end
                        push!(locs,(Lexer.normalize(locnode),r))
                    end
                end

                # Collect indices that the srlocs care about as well
                for srloc in srlocs
                    for e in srloc.sequence
                        if isa(e,Array)
                            push!(indices, e)
                        end
                    end
                end

                # Turn indices into an array. We will use indices into this
                # array to associate them with the highlight information below
                indices = sort(collect(indices), lt = lexless)

                # Remove any indices whose prefix is on the list
                prev = [-1]
                newinds = Any[]
                for ind in indices
                    if !idxstartwiwth(ind, prev)
                        push!(newinds, ind)
                        prev = ind
                    end
                end
                indices = newinds

                ## Now go back and record highlight information
                highlighinfo = Any[Any[] for _ in 1:length(indices)]
                for x in highlighting
                    ind = x[1]
                    while true
                        r = searchsorted(indices, ind, lt = lexless)
                        if length(r) != 0
                            push!(highlighinfo[first(r)],(x[2][(length(ind)+1):end],x[3]))
                            break
                        end
                        if ind == []
                            break
                        end
                        ind = ind[1:end-1]
                    end
                end

                sort!(locs, by = x->Lexer.normalize(x[1]).offset)
                append!(locs,map(
                    Lexer.sortedcomplement(SourceRange(startoffset,stopoffset-startoffset+1,0),
                    filter(x->Lexer.normalize(x).offset >= startoffset, map(x->x[1],locs)))) do x
                    off = x.offset
                    code = interp.code[off+1:(off+x.length)]
                    x, code
                end)
                sort!(locs, by = x->Lexer.normalize(x[1]).offset)

                # Compute necessary data for line numbering
                startline = compute_line(file, startoffset)
                stopline = compute_line(file, stopoffset)
                current_line = startline
                stoplinelength = length(string(stopline))

                print(startline, " "^(stoplinelength-length(string(startline))+1))

                function print_piece(piece, implicit = true)
                    if isa(piece,Array)
                        r = searchsorted(indices, piece, lt = lexless)
                        val = Tree(interp.shadowtree.tree)[piece]
                        if length(r) == 0
                            return
                        end
                        highlights = highlighinfo[first(r)]
                        print(annotate_highlights!(deepcopy(val), highlights))
                    elseif isa(piece,Coloring)
                        print(piece)
                    else
                        lines = split(piece, '\n', keep = true)
                        print_with_color(:white, lines[1])
                        for line in lines[2:end]
                            if !implicit
                                current_line += 1
                                lineno = string(current_line)
                            else
                                lineno = string('↳')
                            end
                            print('\n',lineno, " "^(stoplinelength-length(lineno)+1),line)
                        end
                    end
                end

                for loc in locs
                    if isa(loc[2],SRLoc)
                        for x in loc[2].sequence
                            print_piece(x, true)
                        end
                    else
                        print_piece(loc[2], false)
                    end
                end

                println()
                return
            end
        end
        println(join(map(bytestring,file[startline:stopline]),""))
    end
end

function next_expr!(interp)
    x, next_it = next(interp.it, interp.cur_state)
    interp.cur_state = next_it
    interp.next_expr = get(x[1]), x[2]
end

function find_label_index(tree, label)
    for (i, node) in enumerate(tree.args)
        if isa(node, LabelNode) && node.label == label
            return i
        end
    end
    error("Label not found")
end

function goto!(interp, target)
    # Reset shadowtree
    interp.shadowtree, interp.it = make_shadowtree(interp.ast)
    lind = find_label_index(interp.ast, target)

    next_state  = typeof(interp.cur_state)(
        AbstractTrees.ImplicitNodeStack(Any[Tree(interp.ast)[lind]],
            AbstractTrees.ImplicitIndexStack([lind]))
    )

    # next_expr! below will move past the label node
    interp.cur_state = next(interp.it, next_state)[2]
    next_expr!(interp)
    return done!(interp)
end

const fancy_backtraces = true

function _step_expr(interp)
    ind, node = interp.next_expr
    local ret
    try
        if isa(node, Expr) && node.head == :return
            interp.retval = node.args[1]
            return false
        end
        if isa(node, Slot) || isa(node,SSAValue)
            # Check if we're the LHS of an assignment
            if ind.idx_stack.stack[end] == 1 && interp.shadowtree.tree[ind.idx_stack.stack[1:end-1]].head == :(=)
                ret = node
            elseif isa(node,SSAValue)
                ret = interp.env.ssavalues[node.id+1]
            else
                nslots = length(interp.env.locals)
                id = node.id
                interp.env.last_reference[interp.linfo.slotnames[id]] = id
                if id > nslots
                    # hack: static parameters placed at end of locals for interactive exprs
                    ret = interp.env.sparams[id - nslots]
                else
                    val = interp.env.locals[id]
                    if isnull(val)
                        sym = interp.linfo.slotnames[id]
                        error("local variable $sym not defined")
                    end
                    ret = get(val)
                end
            end
        elseif isa(node, Expr)
            if node.head == :(=)
                lhs = node.args[1]
                rhs = node.args[2]
                if isa(lhs, SSAValue)
                    interp.env.ssavalues[lhs.id+1] = rhs
                elseif isa(lhs, Slot)
                    interp.env.locals[lhs.id] = Nullable{Any}(rhs)
                    interp.env.last_reference[interp.linfo.slotnames[lhs.id]] =
                        lhs.id
                elseif isa(lhs, GlobalRef)
                    eval(lhs.mod,:($(lhs.name) = $(QuoteNode(rhs))))
                end
                # Special case hack for readability.
                # ret = rhs
                ret = node
            elseif node.head == :&
                ret = node
            elseif node.head == :gotoifnot
                ret = node
                if !isa(node.args[1], Bool)
                    throw(TypeError(interp.linfo.def.name, "if", Bool, node.args[1]))
                end
                if !node.args[1]
                    return goto!(interp, node.args[2])
                end
            elseif node.head == :call
                # Don't go through eval since this may have unqouted, symbols and
                # exprs
                f = to_function(node.args[1])
                if isa(f, Core.IntrinsicFunction)
                    # Special handling to quote any literal symbols that may still
                    # be in here, so we can pass it into eval
                    for i in 1:length(node.args)
                        if isa(node.args[i], Union{Symbol,SSAValue,Slot,GlobalRef})
                            node.args[i] = QuoteNode(node.args[i])
                        end
                    end
                    ret = eval(node)
                elseif isa(f, LambdaInfo)
                    ret = finish!(enter_call_expr(interp, node))
                else
                    # Don't go through eval since this may have unqouted, symbols and
                    # exprs
                    ret = f(node.args[2:end]...)
                end
            elseif node.head == :static_typeof
                ret = Any
            elseif node.head == :type_goto
                ret = nothing
            elseif node.head == :enter
                push!(interp.exception_frames, node.args[1])
                ret = node
            elseif node.head == :leave
                for _ = 1:node.args[1]
                    pop!(interp.exception_frames)
                end
                ret = node
            elseif node.head == :static_parameter
                ret = interp.env.sparams[node.args[1]]
            else
                ret = eval(node)
            end
        elseif isa(node, GotoNode)
            return goto!(interp, node.label)
        elseif isa(node, QuoteNode)
            ret = node.value
        else
            ret = eval(node)
        end
    catch err
        (!fancy_backtraces || interp.loctree == nothing) && rethrow(err)
        buf = IOBuffer()
        Base.showerror(buf, err)
        D = diag(SourceRange(),takebuf_string(buf))
        diag(D, Tree(interp.loctree)[ind.idx_stack.stack].loc, "while running this expression", :note)
        rethrow(D)
    end
    _evaluated!(interp, ret)
    true
end
step_expr(interp) = (r = _step_expr(interp); done!(interp); r)

function next_statement!(interp)
    ind, node = interp.next_expr
    move_past = ind[1]
    while step_expr(interp)
        ind, node = interp.next_expr
        if ind[1] != move_past
            return true
        end
    end
    return false
end

function next_until!(f,interp)
    ind, node = interp.next_expr
    while step_expr(interp)
        ind, node = interp.next_expr
        f(node) && return true
    end
    return false
end
next_call!(interp) = next_until!(node->isexpr(node,:call)||isexpr(node,:return), interp)

function changed_line!(expr, line, fls)
    if length(fls) == 1 && isa(expr, LineNumberNode)
        return expr.line != line
    elseif length(fls) == 1 && isa(expr, Expr) && isexpr(expr, :line)
        return expr.args[1] != line
    else
        if is_loc_meta(expr, :pop_loc)
            pop!(fls)
        elseif is_loc_meta(expr, :push_loc)
            push!(fls,(expr.args[2],0))
        end
        return false
    end
end

isgotonode(node) = isa(node, GotoNode) || isexpr(node, :gotoifnot)

"""
Determine whether we are calling a function for which the current function
is a wrapper (either because of optional arguments or becaue of keyword arguments).
"""
function iswrappercall(interp, expr)
    !isexpr(expr, :call) && return false
    r = determine_method_for_expr(interp, expr; enter_generated = false)
    if r !== nothing
        linfo, method, args, _ = r
        ours, theirs = interp.linfo.def, method
        # Check if this a method of the same function that shares a definition line/file.
        # If so, we're likely in an automatically generated wrapper.
        if ours.sig.parameters[1] == theirs.sig.parameters[1] &&
            ours.line == theirs.line && ours.file == theirs.file
            return true
        end
    end
    return false
end

function next_line!(interp; state = nothing)
    didchangeline = false
    fls = determine_line_and_file(interp, idx_stack(interp))
    line = fls[1][2]
    first = true
    while !didchangeline
        ind, node = interp.next_expr
        # Skip evaluated values (e.g. constants)
        while interp.shadowtree.shadow[ind].val
            didchangeline = changed_line!(node, line, fls)
            didchangeline && break
            ind, node = next_expr!(interp)
        end
        didchangeline && break
        ind, node = interp.next_expr
        # If this is a return node, interrupt execution. This is the same
        # special case as in `s`.
        (!first && isexpr(node, :return)) && return true
        first = false
        # If this is a goto node, step it and reevaluate
        if isgotonode(node)
            _step_expr(interp) || return false
            fls = determine_line_and_file(interp, idx_stack(interp))
            didchangeline = line != fls[1][2]
        elseif iswrappercall(interp, node)
            interp.did_wrappercall = true
            interp = enter_call_expr(interp, node)
        else
            _step_expr(interp) || return false
        end
    end
    if state !== nothing && interp != state.interp
        state.interp = state.top_interp = interp
    end
    done!(interp)
    # Ok, we stepped to the next line. Now step through to the next call
    call_or_assignment(node) = isexpr(node,:call) || isexpr(node,:(=)) || isexpr(node, :return)
    call_or_assignment(interp.next_expr[2]) ||
        next_until!(call_or_assignment, interp)
end

function advance_to_line(interp, line)
    while true
        at_line = determine_line_and_file(interp, idx_stack(interp))[1][2]
        at_line == line && break
        next_line!(interp) || break
    end
end

function _evaluated!(interp, ret, wasstaged = false)
    if wasstaged
        # If this is the result of a staged function, we replace the argument
        # the call rather than the call itself
        ind, node = interp.next_expr
        @assert isexpr(node, :call)
        interp.shadowtree[[ind.idx_stack.stack; 1]] = (ret, AnnotationNode{Any}(true,AnnotationNode{Any}[]))
    else
        ind, node = interp.next_expr
        interp.shadowtree[ind.idx_stack.stack] = (ret, AnnotationNode{Any}(true,AnnotationNode{Any}[]))
    end
end
evaluated!(interp, ret, wasstaged = false) = (_evaluated!(interp, ret, wasstaged); done!(interp))

"""
Advance to the next evaluatable statement
"""
function done!(interp)
    ind, node = interp.next_expr
    # Skip evaluated values (e.g. constants)
    while interp.shadowtree.shadow[ind.idx_stack.stack].val
        ind, node = next_expr!(interp)
    end
    return true
end

function to_function(x)
    if isa(x, Function) || isa(x, Core.IntrinsicFunction)
        x
    elseif isa(x, GlobalRef)
        eval(x)
    else
        x
    end
end

_Typeof(x) = isa(x,Type) ? Type{x} : typeof(x)

#TODO: The method should probably just start at it's first definition line
function is_function_def(ex)
    (isa(ex,Expr) && ex.head == :(=) && isexpr(ex.args[1],:call)) ||
    isexpr(ex,:function)
end

function collectcalls(file, parsedexpr, parsedloc, complete = true)
    thecalls = Any[]
    theassignments = Any[]
    forlocs = Any[]
    active_fors = Any[]
    for (ind,node) in indenumerate(PostOrderDFS(parsedexpr))
        parent = Tree(parsedexpr)[ind[1:end-1]]
        if isexpr(node, :call) || isa(node, Tuple)
            push!(thecalls, (Tree(parsedloc)[ind],node))
        end
        # If our parent is a for and we're the first child, insert, the entry
        # iteration lowering now.
        if isexpr(parent, :for) && ind[end] == 1
            forrange = Lexer.normalize(Tree(parsedloc)[ind[1:end-1]])
            argrange = Lexer.normalize(Tree(parsedloc)[ind])
            active_reprange = SourceRange(forrange.offset,argrange.offset+argrange.length-forrange.offset,0)
            l = compute_line(file,active_reprange.offset)
            padding = (active_reprange.offset - file.offsets[l])+1
            # Lowering for x in y
            loc = SRLoc(active_reprange,
                [0, # A = y
                string("\n"," "^padding),
                0, # B = start(A)
                string("\n"," "^padding,"while "),
                0,  # !done(A,B)
                string("\n"," "^(padding+4)),
                0,  # C = next(A,B)
                string("\n"," "^(padding+4)),
                0,  # x = C.1
                string("\n"," "^(padding+4)),
                0], # B = C.2
                0
            )

            push!(active_fors,(loc,active_reprange))

            push!(thecalls,(SourceRange(),nothing))         # Start
            push!(thecalls,(SourceRange(),nothing))         # done
            push!(thecalls,(loc,          nothing))         # !
            push!(thecalls,(SourceRange(),nothing))         # next
            push!(thecalls,(SourceRange(), nothing))        # getfield
            push!(thecalls,(SourceRange(), nothing))        # getfield

            push!(theassignments, (loc, nothing))                  # A = y
            push!(theassignments, (loc, nothing))                  # B = start(A)
            push!(theassignments, (loc, nothing))                  # ssavalue() = next()
            push!(theassignments, (loc, nothing))        # ssavalue() = ans.1
            push!(theassignments, (loc, nothing))        # ssavalue() = ans.2
        elseif isexpr(node, :(=))
            push!(theassignments, (Tree(parsedloc)[ind],node))
        elseif isexpr(node, :(&&))
            push!(theassignments, (SourceRange(), nothing))
            push!(theassignments, (SourceRange(), nothing))
        end
    end
    thecalls = complete ? thecalls[2:end] : thecalls, theassignments, forlocs
end

function expression_mismatch(loweredast, parsedexpr, thecalls, theassignments, forlocs)
    if !fancy_mode
        return nothing
    end
    if loweredast == nothing
        println("Failed to lower expression")
        return
    end
    println("Failed to match parsed and lowered ASTs. This is a bug (or rather missing coverage).")
    println("Falling back to unsugared mode.")
    println("I attempted the following matching:")
    display(Tree(loweredast))
    display(Tree(parsedexpr))
    treemap(PostOrderDFS(loweredast)) do ind, node, childlocs
        if isexpr(node, :call)
            if isempty(thecalls)
                println("Failed to match $node")
                return
            end
            candidate = shift!(thecalls)
            println("Matching call $node with $candidate")
        elseif isexpr(node, :(=))
            if isempty(theassignments)
                println("Failed to match $node")
                return
            end
            println("Matching assignment $node with $(shift!(theassignments))")
        end
        nothing
    end
end

immutable MatchingError
end

function process_loctree(res, contents, linfo, complete = true)
    parsedexpr = Lexer.¬(res)
    parsedloc = Lexer.√(res)
    loweredast = nothing
    local thecalls, theassignments, forlocs
    loctree = try
        res = lower!(res)
        parsedexpr = Lexer.¬(res)
        parsedloc = Lexer.√(res)
        stmts = Base.uncompressed_ast(linfo)
        loweredast = Expr(:body); loweredast.args = stmts
        thecalls, theassignments, forlocs = collectcalls(SourceFile(contents), parsedexpr, parsedloc, complete)
        treemap(PostOrderDFS(loweredast)) do ind, node, childlocs
            if isexpr(node, :call)
                call = node.args[1]
                isempty(thecalls) && throw(MatchingError())
                candidate = shift!(thecalls)
                if isa(candidate[1], SRLoc)
                    ASTInterpreter.sequence!(candidate[1], ind)
                    SourceNode(candidate[1],childlocs)
                else
                    candidate[1]
                end
            elseif isexpr(node, :(=))
                isempty(theassignments) && throw(MatchingError())
                x = shift!(theassignments)
                x = x[1]
                if isa(x, SRLoc)
                    ASTInterpreter.sequence!(x, ind)
                    SourceNode(x,childlocs)
                else
                    x
                end
            else
                SourceNode(SourceRange(),childlocs)
            end
        end
    catch err
        if isa(err, MatchingError)
            expression_mismatch(loweredast, parsedexpr, collectcalls(SourceFile(contents), parsedexpr, parsedloc, complete)...)
        elseif fancy_mode
            rethrow(err)
        end
        nothing
    end

    function postprocess!(loctree, forlocs)
        for (forloc,forloc2) in forlocs
            # Add a location to the parent of !done
            negind = forloc.sequence[5][1:(end-1)]
            newloc1 = deepcopy(forloc)
            newloc2 = deepcopy(forloc)
            # TODO: This really needs to have a better way
            newloc1.sequence[4] = ASTInterpreter.Coloring(newloc1.sequence[4],:yellow)
            Tree(loctree)[negind] = SourceNode(newloc1,Tree(loctree)[negind].children)
        end
    end
    if loctree != nothing
        postprocess!(loctree, forlocs)
        # Make sure we have the whole bounds of the function
        loctree = SourceNode(Tokens.normalize(reduce(⤄,PostOrderDFS(parsedloc))),loctree.children)
    end

    loctree, contents
end

const SEARCH_PATH = [joinpath(JULIA_HOME,"../share/julia/base/"),
    joinpath(JULIA_HOME,"../include/")]
function readfileorhist(file)
    if startswith(string(file),"REPL[")
        hist_idx = parse(Int,string(file)[6:end-1])
        isdefined(Base, :active_repl) || return nothing, ""
        hp = Base.active_repl.interface.modes[1].hist
        return hp.history[hp.start_idx+hist_idx]
    else
        for path in SEARCH_PATH
            fullpath = joinpath(path,string(file))
            if isfile(fullpath)
                return open(readstring, fullpath)
            end
        end
    end
    return nothing
end

function reparse_meth(meth)
    if isa(meth, LambdaInfo)
        meth = meth.def
    elseif isa(meth, TypeMapEntry)
        meth = meth.func
    end
    file, line = Base.find_source_file(string(meth.file)), meth.line
    contents = readfileorhist(file == nothing ? meth.file : file)
    contents == nothing && return (nothing, "")
    buf = IOBuffer(contents)
    for _ in line:-1:2
        readuntil(buf,'\n')
    end
    ts = Lexer.TokenStream{Lexer.SourceLocToken}(buf)
    local res, ts
    try
        res = Parser.parse(ts)
        if !is_function_def(res.expr)
            # Retry parsing the line before
            seekstart(buf)
            for _ in (line-1):-1:2
                readuntil(buf,'\n')
            end
            ts = Lexer.TokenStream{Lexer.SourceLocToken}(buf)
            res = Parser.parse(ts)
        end
    catch err
        if !fancy_mode
            return nothing, contents
        end
        rethrow(err)
    end
    process_loctree(res, contents, meth)
end

function prepare_locals(linfo, argvals = ())
    linfo = linfo.def.lambda_template
    # Construct the environment from the arguments
    argnames = linfo.slotnames[1:linfo.nargs]
    locals = Array(Nullable{Any}, length(linfo.slotflags))
    ng = isa(linfo.ssavaluetypes, Int) ? linfo.ssavaluetypes : length(linfo.ssavaluetypes)
    ssavalues = Array(Any, ng)
    sparams = Array(Any, length(linfo.sparam_syms))
    for i = 1:linfo.nargs
        if linfo.isva && i == length(argnames)
            locals[i] = length(argvals) >= i ? tuple(argvals[i:end]...) : Nullable{Any}(())
            break
        end
        locals[i] = length(argvals) >= i ? Nullable{Any}(argvals[i]) : Nullable{Any}()
    end
    # add local variables initially undefined
    for i = (linfo.nargs+1):length(linfo.slotnames)
        locals[i] = Nullable{Any}()
    end
    Environment(locals, ssavalues, sparams, Dict{Symbol,Int}())
end

function determine_method_for_expr(interp, expr; enter_generated = false)
    f = to_function(expr.args[1])
    allargs = expr.args
    # Extract keyword args
    local kwargs = Expr(:parameters)
    if length(allargs) > 1 && isexpr(allargs[2], :parameters)
        kwargs = splice!(allargs, 2)
    end
    if !isempty(kwargs.args)
        of = f
        f = Core.kwfunc(f)
        allargs = [f,reduce(vcat,Any[[ex.args[1];ex.args[2]] for ex in kwargs.args]),of,
            allargs[2:end]...]
    elseif is(f,Core._apply)
        f = to_function(allargs[2])
        allargs = Base.append_any((allargs[2],), allargs[3:end]...)
    end
    # Can happen for thunks created by generated functions
    if !isa(f, Core.Builtin) && !isa(f, Core.IntrinsicFunction)
        args = allargs[2:end]
        argtypes = Tuple{map(_Typeof,args)...}
        if isa(f, LambdaInfo)
            linfo = f
            method = linfo.def
            enter_generated = false
            lenv = linfo.sparam_vals
        else
            method = try
                which(f, argtypes)
            catch err
                @show typeof(f)
                println(f)
                println(argtypes)
                rethrow(err)
            end
            argtypes = Tuple{_Typeof(f), argtypes.parameters...}
            args = allargs
            sig, tvars = method.sig, method.tvars
            isa(method, TypeMapEntry) && (method = method.func)
            linfo = method.lambda_template
            # Get static parameters
            (ti, lenv) = ccall(:jl_match_method, Any, (Any, Any, Any),
                               argtypes, sig, tvars)::SimpleVector
            if method.isstaged && !enter_generated
                # If we're stepping into a staged function, we need to use
                # the specialization, rather than stepping thorugh the
                # unspecialized method.
                linfo = Core.Inference.specialize_method(method, argtypes, lenv, false)
            else
                if method.isstaged
                    args = map(_Typeof, args)
                end
            end
        end
        return linfo, method, args, lenv
    end
    nothing
end

function enter_call_expr(interp, expr; enter_generated = false)
    r = determine_method_for_expr(interp, expr; enter_generated = enter_generated)
    if r !== nothing
        linfo, method, args, lenv = r
        loctree, code = nothing, ""
        if !(method.isstaged && !enter_generated)
            loctree, code = reparse_meth(method)
        end
        env = prepare_locals(linfo, args)
        # Add static parameters to environment
        for i = 1:length(lenv)
            env.sparams[i] = lenv[i]
        end
        newinterp = enter(method.isstaged && !enter_generated ? linfo : method,
            env, interp != nothing ? interp.stack : Any[],
            loctree = loctree, code = code,
            evaluating_staged = method.isstaged && enter_generated)
        return newinterp
    end
    nothing
end

function print_linfo_desc(io::IO, linfo, specslottypes = false)
    slottypes = linfo.slottypes
    linfo = linfo.def.lambda_template
    argnames = linfo.slotnames[2:linfo.nargs]
    spectypes = specslottypes && (slottypes != nothing) ?
        slottypes[2:linfo.nargs] : Any[Any for i=1:length(argnames)]
    print(io, linfo.def.name,'(')
    first = true
    for (argname, argT) in zip(argnames, spectypes)
        first || print(io, ", ")
        first = false
        print(io, argname)
        !is(argT, Any) && print(io, "::", argT)
    end
    print(io, ") at ",linfo.def.file,":",linfo.def.line)
end

function sym_visible(name)
    sn = string(name)
    if startswith(sn, "#")
        lasthash = rsearchindex(sn, "#")
        if lasthash == 1  # mangled names have 2 '#'s in them,
            return false  # hidden names have 1.
        end
    end
    if sn == "#temp#" || name == :__temp__
        return false
    end
    return true
end

function print_var(io::IO, name, val::Nullable, undef_callback)
    if sym_visible(name)
        print("  | ")
        if isnull(val)
            print(io, name, " = ")
            undef_callback(io,name)
        else
            val = get(val)
            T = typeof(val)
            val = repr(val)
            if length(val) > 150
                val = Suppressed("$(length(val)) bytes of output")
            end
            println(io, name, "::", T, " = ", val)
        end
    end
end

function print_locals(io::IO, linfo, env::Environment, undef_callback)
    linfo = linfo.def.lambda_template
    for i = 1:length(env.locals)
        print_var(io, linfo.slotnames[i], env.locals[i], undef_callback)
    end
    for i = 1:length(env.sparams)
        print_var(io, linfo.sparam_syms[i], Nullable{Any}(env.sparams[i]), undef_callback)
    end
end

function print_backtrace(state, interp)
    num = 1
    for frame in reverse(interp.stack)
        print_frame(state, STDOUT, num, frame)
        num += 1
    end
end

function print_frame(_, io::IO, num, interp::Interpreter)
    print(io, "[$num] ")
    print_linfo_desc(io, interp.linfo)
    println(io)
    print_locals(io, interp.linfo, interp.env,
        (io,name)->println(io, "<undefined>"))
end
print_backtrace(_::Void) = nothing

include(joinpath(dirname(@__FILE__),"..","..","JuliaParser","src","interactiveutil.jl"))

get_env_for_eval(interp::Interpreter) = interp.env
get_linfo(interp::Interpreter) = interp.linfo

function execute_command(interp, command)
    print_with_color(:red,"\nUnknown command!\n");
end

execute_command(state, interp, cmd1, command) = (execute_command(interp, command); return false)

function finish_until!(top_interp, interp)
    while true
        if top_interp == interp
            return interp
        end
        finish!(top_interp)
        oldinterp = pop!(top_interp.stack)
        top_interp = oldinterp.stack[end]
        evaluated!(top_interp, oldinterp.retval)
    end
    interp
end

can_step(_) = false
can_step(interp::Interpreter) = true

loc_if_available(interp, ind) = try
    Tree(interp.loctree)[ind].loc
catch err
    isa(err, BoundsError) || rethrow(err)
    SourceRange()
end

function process_exception!(interp::Interpreter, D::AbstractDiagnostic, istop)
    if !isempty(interp.exception_frames)
        target = interp.exception_frames[end]
        goto!(interp, target)
        ind, _ = interp.next_expr
        diag(D, loc_if_available(interp, ind), "caught here",:note)
        return true
    else
        ind, _ = interp.next_expr
        !istop && diag(D, loc_if_available(interp, ind), "while evaluating expression",:note)
        return false
    end
end

function eval_in_interp(interp, body, slbody = nothing, code = "")
    selfsym = Symbol("#self#")  # avoid having 2 arguments called `#self#`
    unusedsym = Symbol("#unused#")
    env = get_env_for_eval(interp)
    linfo = get_linfo(interp)
    lnames = Any[linfo.slotnames[2:end]..., linfo.sparam_syms...]
    map!(x->(x===selfsym || !sym_visible(x) ? unusedsym : x), lnames)
    lnames = unique(lnames)
    f = Expr(:->,Expr(:tuple,lnames...), body)
    lam = linfo.def.module.eval(f)
    linfo = first(methods(lam)).lambda_template
    # New interpreter is on detached stack
    loctree = nothing
    if slbody != nothing
        loctree, code = process_loctree(slbody, code, linfo, false)
    end
    # Construct a new environment with the arguments inserted
    eval_env = ASTInterpreter.prepare_locals(linfo)
    for varname in lnames
        lidx = findfirst(linfo.slotnames, varname)
        lidx == 0 && continue
        if haskey(env.last_reference, varname)
            eval_env.locals[lidx] = env.locals[env.last_reference[varname]]
        else
            oldidx = findfirst(get_linfo(interp).slotnames, varname)
            if oldidx == 0
                sparamidx = findfirst(get_linfo(interp).sparam_syms, varname)
                sparamidx == 0 && continue
                eval_env.locals[lidx] = env.sparams[sparamidx]
            else
                eval_env.locals[lidx] = env.locals[oldidx]
            end
        end
    end
    stmts = Base.uncompressed_ast(linfo)
    body = Expr(:body); body.args = stmts
    einterp = enter(linfo,body,eval_env,Any[], loctree = loctree, code = code)
    ok, val = try
        true, finish!(einterp)
    catch err
        false, err
    end
    # Here we would reapply the variables
    ok, val
end

type InterpreterState
    top_interp
    interp
    level
    s
    julia_prompt
    main_mode
end

function make_linfo(method, ret)
    func = method.lambda_template
    argnames = [func.slotnames[i] for i = 1:func.nargs]
    lambda = Expr(:lambda, argnames, Expr(Symbol("scope-block"), Expr(:block, ret)))
    if !isempty(func.sparam_syms)
        lambda = Expr(Symbol("with-static-parameters"), lambda, func.sparam_syms...)
    end
    linfo = eval(method.module, lambda)
    linfo.def = method
    linfo
end

# Command Implementation
function done_stepping!(state, interp; to_next_call = false)
    s = state.s
    stack = state.interp.stack
    this_idx = findfirst(stack, interp)
    if this_idx == 0
        LineEdit.transition(s, :abort)
        interp = nothing
    else
        oldinterp = state.interp
        state.top_interp = state.interp = this_idx == 1 ? nothing : stack[this_idx-1]
        resize!(stack, this_idx-1)
        if !isa(state.interp, Interpreter)
            LineEdit.transition(s, :abort)
            return nothing
        end
        ret = oldinterp.retval
        if oldinterp.evaluating_staged
            ret = make_linfo(get_linfo(oldinterp).def, oldinterp.retval)
        end
        evaluated!(state.interp, ret, oldinterp.evaluating_staged)
        # For wrappercall, also pop the wrapper
        if to_next_call && state.interp.did_wrappercall
            finish!(state.interp)
            done_stepping!(state, state.interp; to_next_call = to_next_call)
            return interp
        end
        to_next_call &&
          (isexpr(state.interp.next_expr[2], :call) || next_call!(state.interp))
    end
    interp
end

function execute_command(state, interp::Interpreter, ::Val{:?}, cmd)
    display(
            Base.@md_str """
    Basic Commands:\\
    - `n` steps to the next line\\
    - `s` steps into the next call\\
    - `finish` runs to the end of the function\\
    - `bt` shows a simple backtrace\\
    - ``` `stuff ``` runs `stuff` in the current frame's context\\
    - `fr v` will show all variables in the current frame\\
    - `f n` where `n` is an integer, will go to the `n`-th frame.

    Advanced commands:\\
    - `nc` steps to the next call\\
    - `ns` steps to the next statement\\
    - `se` does one expression step\\
    - `si` does the same but steps into a call if a call is the next expression\\
    - `sg` steps into a generated function\\
    - `shadow` shows the internal representation of the expression tree\\
       (for debugger debugging only)\\
    - `loc` shows the column data for the current top frame,\\
        in the same format as JuliaParsers's testshell.\\
    """)
    return false
end

function execute_command(state, interp::Interpreter, ::Val{:finish}, cmd)
    finish!(state.interp)
    done_stepping!(state, state.interp; to_next_call = true)
    return true
end

function execute_command(state, interp::Interpreter, ::Val{:edit}, cmd)
    file, line = determine_line_and_file(interp, idx_stack(interp))[end]
    dfile = Base.find_source_file(string(file))
    file = dfile == nothing ? file : dfile
    Base.edit(file, line)
    return false
end

function execute_command(state, interp, ::Val{:bt}, cmd)
    print_backtrace(state, state.top_interp)
    println()
    return false
end

function execute_command(state, interp::Interpreter, ::Val{:shadow}, cmd)
    print_shadowtree(interp.shadowtree, idx_stack(interp))
    println()
    return false
end

function execute_command(state, interp, ::Val{:linfo}, cmd)
    eval(Main,:(linfo = $(get_linfo(interp))))
    LineEdit.transition(state.s, :abort)
    return false
end

function execute_command(state, interp::Interpreter, ::Val{:ind}, cmd)
    println("About to execute index ", idx_stack(interp))
    return false
end

function execute_command(state, interp::Interpreter, ::Val{:loc}, cmd)
    if interp.loctree == nothing
        print_with_color(:red, STDERR, "No loctree available\n")
        return true
    end
    w = create_widget(state.interp.loctree, state.interp.code)
    TerminalUI.print_snapshot(TerminalUI.InlineDialog(w,
        Base.Terminals.TTYTerminal("xterm", STDIN, STDOUT, STDERR)
        ))
    return false
end

function execute_command(state, interp, ::Val{:up}, cmd)
    new_stack_idx = length(state.top_interp.stack)-state.level
    if new_stack_idx == 0
        print_with_color(:red, STDERR, "Already at the top of the stack\n")
        return false
    end
    state.level += 1
    state.interp = state.top_interp.stack[new_stack_idx]
    return true
end

function execute_command(state, interp, ::Val{:down}, cmd)
    new_stack_idx = length(state.top_interp.stack)-(state.level-2)
    if new_stack_idx > length(state.top_interp.stack)
        print_with_color(:red, STDERR, "Already at the bottom of the stack\n")
        return false
    end
    state.level -= 1
    state.interp = state.top_interp.stack[new_stack_idx]
    return true
end

function execute_command(state, interp, ::Union{Val{:f},Val{:fr}}, command)
    subcmds = split(command,' ')[2:end]
    if isempty(subcmds) || subcmds[1] == "v"
        print_frame(state, STDOUT, state.level, state.interp)
        return false
    else
        new_level = parse(Int, subcmds[1])
        new_stack_idx = length(state.top_interp.stack)-(new_level-1)
        if new_stack_idx > length(state.top_interp.stack) || new_stack_idx < 1
            print_with_color(:red, STDERR, "Not a valid frame index\n")
            return false
        end
        state.level = new_level
        state.interp = state.top_interp.stack[new_stack_idx]
    end
    return true
end

function execute_command(state, interp::Interpreter, cmd::Union{Val{:s},Val{:si},Val{:sg}}, command)
    first = true
    while true
        expr = state.interp.next_expr[2]
        if isa(expr, Expr)
            if expr.head == :call && !isa(expr.args[1],Core.IntrinsicFunction)
                x = enter_call_expr(state.top_interp, expr, enter_generated = command == "sg")
                if x !== nothing
                    state.interp = state.top_interp = x
                    (cmd == Val{:s}() || cmd == Val{:sg}()) && next_call!(x)
                    return true
                end
            elseif !first && isexpr(expr, :return)
                # As a special case, do not step through a return
                # statement, unless the user was already there when they
                # hit `s`
                return true
            end
        end
        first = false
        command == "si" && break
        if !step_expr(state.interp)
            done_stepping!(state, state.top_interp; to_next_call = true)
            return true
        end
    end
    execute_command(state, interp, Val{:se}(), "se")
end

function execute_command(state, interp::Interpreter, ::Union{Val{:ns},Val{:nc},Val{:n},Val{:se}}, command)
    (state.top_interp != interp) && (state.top_interp = finish_until!(state.top_interp, interp))
    state.level = 1
    if command == "ns" ? !next_statement!(state.interp) :
       command == "nc" ? !next_call!(state.interp) :
       command == "n" ? !next_line!(state.interp; state = state) :
        !step_expr(state.interp) #= command == "se" =#
        done_stepping!(state, state.interp; to_next_call = command == "n")
        return true
    end
    return true
end

function eval_code(state, command)
    local theerr = nothing
    res = try
        ts = Lexer.TokenStream{Lexer.SourceLocToken}(command)
        ts.filename = "REPL"
        res = Main.JuliaParser.Parser.parse(ts)
    catch err
        theerr = err
    end
    if theerr == nothing
        body = Lexer.¬(res)
        ok, result = eval_in_interp(state.interp, body, res, command)
        ok && return (ok, result)
        theerr = result
    end
    if !isa(theerr, AbstractDiagnostic)
        REPL.display_error(STDERR, theerr, Base.catch_backtrace())
    else
        display_diagnostic(STDERR, command, theerr)
    end
    REPL.println(STDERR); REPL.println(STDERR)
    return false, nothing
end
eval_code(state, buf::IOBuffer) = eval_code(state, takebuf_string(buf))

function language_specific_prompt(state, stack::Interpreter)
    state.julia_prompt
end

function run_pre_hooks(state, interp::Interpreter)
    if !isempty(interp.prehook)
        ok, res = eval_code(state, interp.prehook)
        ok && res !== nothing && (display(res); println())
    end
end
run_pre_hooks(state, interp) = nothing

function execute_command(state, interp::Interpreter, ::Val{:prehook}, cmd)
    idx = findfirst(cmd, ' ')
    cmd = idx == 0 ? "" : cmd[idx+1:end]
    interp.prehook = cmd
    return false
end

function RunDebugREPL(top_interp)
    promptname(level, name) = "$level|$name > "

    repl = Base.active_repl
    state = InterpreterState(top_interp, top_interp, 1, nothing, nothing, nothing)

    # Setup debug panel
    panel = LineEdit.Prompt(promptname(state.level, "debug");
        prompt_prefix="\e[38;5;166m",
        prompt_suffix=Base.text_colors[:white],
        on_enter = s->true)

    # For now use the regular REPL completion provider
    replc = Base.REPL.REPLCompletionProvider(repl)

    # Set up the main Julia prompt
    julia_prompt = LineEdit.Prompt(promptname(state.level, "julia");
        # Copy colors from the prompt object
        prompt_prefix = repl.prompt_color,
        prompt_suffix = (repl.envcolors ? Base.input_color : repl.input_color),
        complete = replc,
        on_enter = Base.REPL.return_callback)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:debug => panel,
        :julia => julia_prompt))
    julia_prompt.hist = panel.hist
    Base.REPL.history_reset_state(panel.hist)

    search_prompt, skeymap = Base.LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = Base.REPL.LatexCompletions()

    state.interp = top_interp
    state.julia_prompt = julia_prompt
    state.main_mode = panel

    panel.on_done = (s,buf,ok)->begin
        state.s = s
        line = takebuf_string(buf)
        old_level = state.level
        if !ok || strip(line) == "q"
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if isempty(strip(line)) && length(panel.hist.history) > 0
            command = panel.hist.history[end]
        else
            command = strip(line)
        end
        do_print_status = true
        cmd1 = split(command,' ')[1]
        do_print_status = try
            execute_command(state, state.interp, Val{Symbol(cmd1)}(), command)
        catch err
            isa(err, AbstractDiagnostic) || rethrow(err)
            caught = false
            for interp_idx in length(state.top_interp.stack):-1:1
                if process_exception!(state.top_interp.stack[interp_idx], err, interp_idx == length(top_interp.stack))
                    interp = state.top_interp = state.top_interp.stack[interp_idx]
                    resize!(state.top_interp.stack, interp_idx)
                    caught = true
                    break
                end
            end
            !caught && rethrow(err)
            display_diagnostic(STDERR, state.interp.code, err)
            println(STDERR)
            LineEdit.reset_state(s)
            return true
        end
        if old_level != state.level
            panel.prompt = promptname(state.level,"debug")
            julia_prompt.prompt = promptname(state.level,"julia")
        end
        LineEdit.reset_state(s)
        if do_print_status
            print_status(state, state.interp)
            run_pre_hooks(state, state.interp)
        end
        return true
    end

    const all_commands = ("q", "s", "si", "finish", "bt", "loc", "ind", "shadow",
        "up", "down", "ns", "nc", "n", "se")

    julia_prompt.on_done = (s,buf,ok)->begin
        if !ok
            LineEdit.transition(s, :abort)
            return false
        end
        xbuf = copy(buf)
        ok, result = eval_code(state, buf)
        if ok
            display(result)
            println(); println()
        else
            command = takebuf_string(buf)
            # Convenience hack. We'll see if this is more useful or annoying
            for c in all_commands
                !startswith(command, c) && continue
                LineEdit.transition(s, panel)
                LineEdit.state(s, panel).input_buffer = xbuf
                break
            end
        end
        LineEdit.reset_state(s)
        return true
    end

    const repl_switch = Dict{Any,Any}(
        '`' => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                prompt = language_specific_prompt(state, state.interp)
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, prompt) do
                    LineEdit.state(s, prompt).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s,key)
            end
        end
    )

    b = Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap([repl_switch;b])
    julia_prompt.keymap_dict = LineEdit.keymap([Base.REPL.mode_keymap(panel);b])

    # Skip evaluated values (e.g. constants)
    done!(state.interp)
    print_status(state, state.interp)
    Base.REPL.run_interface(repl.t, LineEdit.ModalInterface([panel,julia_prompt,search_prompt]))
end

idx_stack(interp::Interpreter) = interp.next_expr[1].idx_stack.stack

"""
Single step the current function to exit, optionally printing the tree at every
step.

Returns the return value of the interpreted function.
"""
function finish!(interp::Interpreter; print_step::Bool = false, recursive = false)
    while true
        if recursive
            expr = interp.next_expr[2]
            if isa(expr, Expr)
                if expr.head == :call
                    newinterp = enter_call_expr(interp, expr)
                    if newinterp !== nothing
                        finish!(newinterp; print_step = print_step, recursive = true)
                        evaluated!(interp, newinterp.retval)
                        print_step && print_shadowtree(interp.shadowtree, idx_stack(interp))
                        continue
                    end
                end
            end
        end
        if !step_expr(interp)
            break
        end
        print_step && print_shadowtree(interp.shadowtree, idx_stack(interp))
    end
    interp.retval
end

macro enter(arg)
    arg = ASTInterpreter.lower!(arg)
    @assert isa(arg, Expr) && arg.head == :call
    quote
        theargs = $(esc(Expr(:tuple,map(
            x->isexpr(x,:parameters)?QuoteNode(x):x,arg.args)...)))
        interp = ASTInterpreter.enter_call_expr(nothing,Expr(:call,theargs...))
        ASTInterpreter.next_call!(interp)
        ASTInterpreter.RunDebugREPL(interp)
    end
end

include("lowering.jl")
include("precompile.jl")
_precompile_()

end


text:
"""
    robust_pmap(f::Function, xs::Any, num_retries::Int=3)

A pmap function with retries for common network errors.
"""
function robust_pmap(f::Function, args...; num_retries::Int=3)

    # This function returns true if we should retry and fails if we should error out.
    function retry_check(delay_state, err)
        # Below each condition to retry is listed along with an explanation about why
        # retrying should/might work.
        should_retry = (
            # Worker death is normally stocastic, if not then doesn't matter how many
            # retries as it will rapidly kill all workers
            err isa ProcessExitedException ||
            # If we are in the middle of fetching data and the process is killed we could
            # get an ArgumentError saying that the stream was closed or unusable.
            # So same as above.
            err isa ArgumentError && occursin("stream is closed or unusable", err.msg) ||
            # In general IO errors can be transient and related to network blips
            err isa Base.IOError
        )
        if should_retry
            info(LOGGER, "Retrying computation that failed due to a $(typeof(err)): $err")
        else
            warn(LOGGER, "Non-retryable $(typeof(err)) occurred: $err")
        end
        return should_retry
    end

    return pmap(f, CachingPool(workers()), args...;
        retry_check=retry_check,
        retry_delays=ExponentialBackOff(n=num_retries)
    )
end

ext: jl


text:
using Images, ImageMagick

function write_image(pixels::Array{HSV, 2}, output_file::String)
    save(output_file, pixels)
end

function f(x::Complex{Float64})
    r = sqrt(real(x)^2 + imag(x)^2)
    theta = atan(imag(x), real(x))
    complex_val = r*exp(theta*im)
    return complex_val
end

function f2(x::Complex{Float64})
    r = sqrt(real(x)^2 + imag(x)^2)
    theta = atan(imag(x), real(x))
    complex_val = (r*exp(theta*im))^2
    return complex_val
end

function poly(x::Complex{Float64})
    r = sqrt(real(x)^2 + imag(x)^2)
    theta = atan(imag(x), real(x))
    complex_val = r*exp(theta*im)
    x = complex_val
    #return (x^2 - 1)*((x - 2 - 1im)^2) / (x^2 + 2 + 2im)
    return x^3 - 1
end

function inverse(x::Complex{Float64})
    r = sqrt(real(x)^2 + imag(x)^2)
    theta = atan(imag(x), real(x))
    complex_val = r*exp(theta*im)
    if abs(complex_val) < 0.005
        complex_val = 0.004 + 0.003im
    end
    return 1/complex_val
end

function draw_grid(pixel_color::HSV, z::Complex{Float64}, threshold::Float64)

    shade_real = (abs(sin(real(z)*pi))^threshold)
    shade_imag = (abs(sin(imag(z)*pi))^threshold)

    return HSV(pixel_color.h,
               pixel_color.s,
               pixel_color.v*shade_real*shade_imag)
end

function contour(pixel_color::HSV, z::Complex{Float64})
    factor = 1-((abs(z) % 1)*0.5)
    pixel_color = HSV(pixel_color.h, pixel_color.s, pixel_color.v*factor)
    return pixel_color
end

function simple_domain(res::Int64, xmax::Float64, alpha::Float64,
                       output_file::String, f, threshold)
    a = Array{HSV,2}(undef, res, res)
    for i = 1:res
        for j = 1:res
            x = -xmax + 2*i*xmax/res + (-xmax + 2*j*xmax/res)im
            z = f(x)
            magnitude = abs(z)
            argument = angle(z) + pi

            a[i,j] = HSV(argument/(2*pi)*360, 1-alpha^(magnitude), 1)
            a[i,j] = draw_grid(a[i,j], z, threshold)
            a[i,j] = contour(a[i,j], z)
        end
    end
    write_image(a, output_file)
end

function write_domain(res::Int64, xmax::Float64, real_output::String,
                      imaginary_output::String, f)
    a = Array{Complex,2}(undef, res, res)

    imaginary_f = open(imaginary_output, "w")
    real_f = open(real_output, "w")

    for i = 1:res
        for j = 1:res
            x = -xmax + 2*i*xmax/res + (-xmax + 2*j*xmax/res)im
            z = f(x)
            if j != res
                write(imaginary_f, string(imag(z))*"\t")
                write(real_f, string(real(z))*"\t")
            else
                write(imaginary_f, string(imag(z)))
                write(real_f, string(real(z)))
            end
        end
        write(imaginary_f, "\n")
        write(real_f, "\n")
    end

    close(imaginary_f)
    close(real_f)
end

max_stars_repo_name: Amynka/simuleios

max_issues_repo_name: Amynka/simuleios

max_forks_repo_name: Amynka/simuleios


text:
module IntervalStatistic

include("Mean.jl")
include("Variance.jl")
include("Check.jl")

function mean(values :: Vector, method :: Mean.BaseEstimator)
    Mean.estimate(values, method)
end

function var(values :: Vector, method :: Variance.BaseEstimator)
    Variance.estimate(values, method)
end

function isDistribution(values :: Vector, method :: Check.BaseCheck)
    Check.isDistribution(values, method)
end

end


text:
module SIR

using RecordedArrays
using Random
using Test
using EvolutionaryModelingTools

# parameters
clock = ContinuousClock(300.0)
β = 0.001
γ = 0.01

# state
S = recorded(DynamicEntry, clock, 100)
I = recorded(DynamicEntry, clock, 1)
R = recorded(DynamicEntry, clock, 0)

# simulate with this package
@reaction infection begin
    β * S * I
    (S[] -= 1; I[] += 1)
end

@reaction recovery begin
    γ * I
    (I[] -= 1; R[] += 1)
end

ps = (; β, γ, S, I, R)
rs = (infection, recovery)

Random.seed!(1)
init!(clock)

ps_updated = gillespie(clock, ps, rs)[1]

# simulate manually
Random.seed!(1)
init!(clock)

for _ in clock
    infection_rate = β * S * I
    recovery_rate = γ * I

    summed = infection_rate + recovery_rate

    τ = -log(rand()) / summed # compute time interval

    increase!(clock, τ) # update current time

    # sample a reaction and adjust population size
    rn = rand() * summed
    if rn < infection_rate
        S[] -= 1
        I[] += 1
    else
        I[] -= 1
        R[] += 1
    end
end

for sym in (:S, :I, :R)
    @eval @test state($sym) == state(ps_updated.$sym)
    @eval @test getts(getentries($sym)) == getts(getentries(ps_updated.$sym))
    @eval @test getvs(getentries($sym)) == getvs(getentries(ps_updated.$sym))
end

end


text:
using CalibrationErrors
using Distances
using Distributions
using StatsBase

using LinearAlgebra
using Random
using Statistics
using Test

using CalibrationErrors: unsafe_skce_eval, unsafe_ucme_eval

Random.seed!(1234)

@testset "CalibrationErrors" begin
    @testset "binning" begin
        @testset "generic" begin
            include("binning/generic.jl")
        end
        @testset "uniform" begin
            include("binning/uniform.jl")
        end
        @testset "median variance" begin
            include("binning/medianvariance.jl")
        end
    end

    @testset "ECE" begin
        include("ece.jl")
    end

    @testset "Kernels" begin
        include("kernels.jl")
    end

    @testset "SKCE" begin
        @testset "generic" begin
            include("skce/generic.jl")
        end
        @testset "biased" begin
            include("skce/biased.jl")
        end
        @testset "unbiased" begin
            include("skce/unbiased.jl")
        end
    end

    @testset "UCME" begin
        include("ucme.jl")
    end

    @testset "deprecations" begin
        include("deprecated.jl")
    end
end


text:
################################################################################
#
#   FUNCTIONS FOR PLOTTING OF BRILLOUIN ZONES
#
################################################################################

# overwriting the plot function of PyPlot (2D)
function PyPlot.plot(
            bz :: BZ
            ;
            BZ_corners :: Bool = true,
            BZ_color :: Vector{<:Integer} = [10,10,10],
            BZ_corner_size :: Real = 1,
            BZ_linewidth :: Real = 1,
            kwargs...
        ) where {L,N,P<:AbstractReciprocalPoint{2},B<:AbstractBond{L,N},R<:AbstractReciprocalUnitcell{P,B},BZ<:AbstractBrillouinZone{R}}

    # scatter all corners
    if BZ_corners
        # scatter all corners
        scatter(
            [c[1] for c in corners(bz)],
            [c[2] for c in corners(bz)],
            color=BZ_color./255
        )
    end
    # draw the surrounding faces
    for f in faces(bz)
        # obtain lists of x y and z values
        xvals = [corners(bz)[i][1] for i in f]
        yvals = [corners(bz)[i][2] for i in f]
        # push the first element into the lists to close the loops
        push!(xvals, xvals[1])
        push!(yvals, yvals[1])
        # plot the BZ face
        plot(
            xvals, yvals,
            color = BZ_color./255,
            linewidth = BZ_linewidth,
            linestyle = "-",
        )
    end
end

# overwriting the plot function of PyPlot (3D)
function PyPlot.plot(
            bz :: BZ
            ;
            BZ_corners :: Bool = true,
            BZ_color :: Vector{<:Integer} = [10,10,10],
            BZ_corner_size :: Real = 1,
            BZ_linewidth :: Real = 1,
            kwargs...
        ) where {L,N,P<:AbstractReciprocalPoint{3},B<:AbstractBond{L,N},R<:AbstractReciprocalUnitcell{P,B},BZ<:AbstractBrillouinZone{R}}

    # scatter all corners
    if BZ_corners
        # scatter all corners
        scatter3D(
            [c[1] for c in corners(bz)],
            [c[2] for c in corners(bz)],
            [c[3] for c in corners(bz)],
            color=BZ_color./255
        )
    end
    # draw the surrounding faces
    for f in faces(bz)
        # obtain lists of x y and z values
        xvals = [corners(bz)[i][1] for i in f]
        yvals = [corners(bz)[i][2] for i in f]
        zvals = [corners(bz)[i][3] for i in f]
        # push the first element into the lists to close the loops
        push!(xvals, xvals[1])
        push!(yvals, yvals[1])
        push!(zvals, zvals[1])
        # plot the BZ face
        plot3D(
            xvals, yvals, zvals,
            color = BZ_color./255,
            linewidth = BZ_linewidth,
            linestyle = "-",
        )
    end
end

# export the plot function
export plot

# Explicit fancy function to plot a Brillouin zone
function plotBrillouinZone(
            bz :: BZ
            ;
            figsize :: Tuple{<:Real, <:Real} = (6,6),
            coordinate_system :: Bool = false,
            kwargs...
        ) where {D,L,N,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},R<:AbstractReciprocalUnitcell{P,B},BZ<:AbstractBrillouinZone{R}}

    # create a new figure
    fig = PyPlot.figure(figsize = figsize)

    # plot the lattice
    plot(bz; kwargs...)

    # get the current axis
    ax = fig.gca()
    # turn off everthing but the plot
    if !coordinate_system
        ax.set_axis_off()
    end
    # equal aspect ratio
    ax.set_aspect("equal")
    # tighten layout
    tight_layout()
    # show
    show()
end

# export the plot function
export plotBrillouinZone


text:
"""
    SimpleFlap <: NoStateModel

Linear, steady-state, two-dimensional control surface model with parameters
``c_{n,\\delta}``, ``c_{a,\\delta}``, and ``c_{m,\\delta}``.
"""
struct SimpleFlap <: NoStateModel end

"""
    SimpleFlap()

Initialize an object of type [`SimpleFlap`](@ref)
"""
SimpleFlap()

# --- Traits --- #

number_of_parameters(::Type{SimpleFlap}) = 3

# --- Convenience Methods --- #

function set_parameters!(p, model::SimpleFlap; cnd, cad, cmd)

    p[1] = cnd
    p[2] = cad
    p[3] = cmd

    return p
end

function separate_parameters(model::SimpleFlap, p)

    return (cnd=p[1], cad=p[2], cmd=p[3])
end

# --- Internal Methods --- #

function simpleflap_loads(b, U, ρ, cnδ, caδ, cmδ, δ)
    L = ρ*U^2*b*cnδ*δ
    D = ρ*U^2*b*caδ*δ
    M = 2*ρ*U^2*b^2*cmδ*δ
    return L, D, M
end


text:
using SatelliteImages
using Documenter

DocMeta.setdocmeta!(SatelliteImages, :DocTestSetup, :(using SatelliteImages); recursive=true)

makedocs(;
    modules=[SatelliteImages],
    authors="Ayush Patnaik, Ajay Shah, Anshul Tayal, Susan Thomas",
    repo="https://github.com/ayushpatnaikgit/SatelliteImages.jl/blob/{commit}{path}#{line}",
    sitename="SatelliteImages.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ayushpatnaikgit.github.io/SatelliteImages.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ayushpatnaikgit/SatelliteImages.jl",
    devbranch="main",
)


text:
include("testcalls.jl")
include("testoparser.jl")
include("testgparser.jl")
include("testbparser.jl")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using FEMBasis
using Test

function test_elements(elements::Symbol...)
    for i in 1:length(elements)
        element = getfield(FEMBasis, elements[i])
        N = zeros(size(element, 2))
        X = get_reference_element_coordinates(element)
        for i=1:length(element)
            eval_basis!(element, N, X[i])
            N_expected = zeros(length(element))
            N_expected[i] = 1.0
            @test isapprox(N, N_expected)
        end
    end
end

@testset "Continuous Lagrange elements" begin
    @testset "Segments" begin test_elements(:Seg2, :Seg3) end
    @testset "Quadrangles" begin test_elements(:Quad4, :Quad8, :Quad9) end
    @testset "Triangles" begin test_elements(:Tri3, :Tri6, :Tri7) end
    @testset "Tetrahedrons" begin test_elements(:Tet4, :Tet10) end
    @testset "Hexahedrons" begin test_elements(:Hex8, :Hex20, :Hex27) end
    @testset "Wedges" begin test_elements(:Wedge6, :Wedge15) end
    @testset "Pyramids" begin test_elements(:Pyr5, ) end
    # fixme: partition of unity is not satisfied
    #@testset "Pyramids" begin test_elements(:Pyr5A, ) end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Xdmf.jl/blob/master/LICENSE

module Xdmf
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

module Blake3Hash

ENDIAN_BOM != 0x04030201 && throw(
    "Sorry this package currently only works on little endian machines."
)

include("ref.jl")
include("sa_impl/_module.jl")
include("simd_128/_module.jl")

using .Simd128
const Blake3Ctx = Simd128.Blake3Ctx
const update!   = Simd128.update!
const digest    = Simd128.digest

export Blake3Ctx, update!, digest

end # module


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

using MacroTools

"""
    @rec [index=]range expression

Computes recursive operations on time series. The first argument is the range
and the second argument is an expression to be evaluated over that range.

The expression is meant to be an assignment, but it doesn't have to be.

The the range can specify an optional indexing variable (as in a for loop). If
not given, the variable is assumed to be `t`.

### Examples
```julia-repl
julia> s = TSeries(1U)
Empty TSeries{Unit} starting 5U

julia> s[1U] = s[2U] = 1; s
2-element TSeries{Unit} with range 1U:2U:
      1U : 1.0
      2U : 1.0

julia> @rec t=3U:10U s[t] = s[t-1] + s[t-2]

julia> s
10-element TSeries{Unit} with range 1U:10U:
      1U : 1.0
      2U : 1.0
      3U : 2.0
      4U : 3.0
      5U : 5.0
      6U : 8.0
      7U : 13.0
      8U : 21.0
      9U : 34.0
     10U : 55.0
```
"""
macro rec(arg_rng, arg_eqn)
    ind = nothing
    rng = nothing
    matched = @capture(arg_rng, ind_ = rng_) || @capture(arg_rng, ind_ in rng_) || @capture(arg_rng, ind_ ∈ rng_)
    if !matched
        rng = arg_rng
        @capture(arg_eqn, var_[ind_] = rhs_)
    end
    if !isa(ind, Symbol)
        ind = :t
    end
    return esc(quote
        for $(ind) in $(rng)
            $(arg_eqn)
        end
    end)
end


text:
NIScope_Measurements_Dict = Dict()
# Base on IVIScope measurement enum
NIScope_Measurements_Dict["RiseTime"]  		   = 0
NIScope_Measurements_Dict["FallTime"]  	 	   = 1
NIScope_Measurements_Dict["Frequency"] 		   = 2
NIScope_Measurements_Dict["Period"]     	   = 3
NIScope_Measurements_Dict["VoltageRMS"] 	   = 4
NIScope_Measurements_Dict["VoltagePeakToPeak"] = 5
NIScope_Measurements_Dict["VoltageMax"]  	   = 6
NIScope_Measurements_Dict["VoltageMin"]  	   = 7
NIScope_Measurements_Dict["VoltageHigh"]  	   = 8
NIScope_Measurements_Dict["VoltageLow"]  	   = 9
NIScope_Measurements_Dict["VoltageAverage"]    = 10
NIScope_Measurements_Dict["WidthNeg"] 		   = 11
NIScope_Measurements_Dict["WidthPos"] 		   = 12
NIScope_Measurements_Dict["DutyCycleNeg"] 	   = 13
NIScope_Measurements_Dict["DutyCyclePos"]      = 14
NIScope_Measurements_Dict["Amplitude"] 		   = 15
NIScope_Measurements_Dict["VoltageCycleRMS"]   = 16
NIScope_Measurements_Dict["VoltageCycleAverage"] = 17
NIScope_Measurements_Dict["OverShoot"] 		   = 18
NIScope_Measurements_Dict["Preshoot"]  		   = 19
# Set Measurement
function set_meas(obj::T;ch::String="1",fct::String="max",par::Int=1,vw::Int=1) where {T<:NIScope}
	error("not implemented yet")
end
# Get Measurement
function get_meas(obj::T;par::Int=1) where {T<:NIScope}
	error("not implemented yet")
end

# Read measurement
# Initiates an acquisition, waits for it to complete, and performs the specified waveform measurement for a single channel and record or for multiple channels and records.
#"""
#meas = niScope_ReadMeasurement(obj::T,ch="1",fct=7) where {T<:NIScope}
#"""
#function niScope_ReadMeasurement(obj::T;ch="",fct=7) where {T<:NIScope}
# niScope_ReadMeasurement (ViSession vi, ViConstString chlist, ViReal64 timeout, ViInt32 scalarMeasFunction, ViReal64* result);
#	chlist = [UInt8.(collect("0")); UInt8(0)] # terminate with NULL char
#	timeout = 0.0
#	result = Ref(0.0)
#	niScope_ReadMeasurement(obj.handle, chlist, timeout, Int32(fct),result)
#	result[]
#end


text:
"""copy a FEFunction on a FE Space"""
clone_fe_function(space,f)=FEFunction(space,copy(get_free_dof_values(f)))

"""
 Given an arbitrary number of FEFunction arguments,
 returns a tuple with their corresponding free DoF Values
"""
function Gridap.get_free_dof_values(functions...)
  map(get_free_dof_values,functions)
end

function setup_mixed_spaces(model, order)
  reffe_rt  = ReferenceFE(raviart_thomas, Float64, order)
  V = FESpace(model, reffe_rt ; conformity=:HDiv)
  U = TrialFESpace(V)
  reffe_lgn = ReferenceFE(lagrangian, Float64, order)
  Q = FESpace(model, reffe_lgn; conformity=:L2)
  P = TrialFESpace(Q)
  reffe_lgn = ReferenceFE(lagrangian, Float64, order+1)
  S = FESpace(model, reffe_lgn; conformity=:H1)
  R = TrialFESpace(S)

  R, S, U, V, P, Q
end

function setup_and_factorize_mass_matrices(dΩ, R, S, U, V, P, Q)
  amm(a,b) = ∫(a⋅b)dΩ
  H1MM = assemble_matrix(amm, R, S)
  RTMM = assemble_matrix(amm, U, V)
  L2MM = assemble_matrix(amm, P, Q)
  H1MMchol = lu(H1MM)
  RTMMchol = lu(RTMM)
  L2MMchol = lu(L2MM)

  H1MM, RTMM, L2MM, H1MMchol, RTMMchol, L2MMchol
end

function new_vtk_step(Ω,file,_cellfields)
  n = size(_cellfields)[1]
  createvtk(Ω,
            file,
            cellfields=_cellfields,
            nsubcells=n)
end


text:
module MPIReco
  using Reexport
  @reexport using ImageUtils
  using IniFile
  @reexport using MPIFiles
  @reexport using RegularizedLeastSquares
  @reexport using Images

  using AxisArrays
  using ProgressMeter
  using LinearAlgebra
  using SparseArrays
  @reexport using Unitful
  using Distributed
  using DistributedArrays
  # using TensorDecompositions
  import LinearOperators
  using IniFile
  import LinearAlgebra: ldiv!, \
  using DSP

  include("Utils.jl")
  include("MultiContrast.jl")
  include("RecoParameters.jl")
  include("SystemMatrixCenter.jl")
  include("SystemMatrix.jl")
  include("Weighting.jl")
  include("Reconstruction.jl")
  include("MultiPatch.jl")
  include("SystemMatrixRecovery.jl")
  include("MotionCompensation/MotionCompensation.jl")
end # module


text:
module Day3

using Advent: read_charmatrix, save_input

function slope(mat, (right, down))
    xs = 1:down:size(mat, 1)
    ys = range(1, step=right, length=length(xs))
    I = CartesianIndex.(xs, mod1.(ys, size(mat, 2)))
    count(c -> mat[c] == '#', I)
end

function solve_a(filename)
    mat = read_charmatrix(filename)
    slope(mat, (3, 1))
end

function solve_b(filename)
    mat = read_charmatrix(filename)
    slopes = ((1, 1), (3, 1), (5, 1), (7, 1), (1, 2))
    prod(Base.Fix1(slope, mat), slopes)
end

function solve_big(filename = "day3_big.in")
    mat = read_charmatrix(filename)
    right = (2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 32, 36, 48, 54, 64)
    down = (1, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 47)
    prod(Base.Fix1(slope, mat), Iterators.product(right, down))
end

solve(filename = save_input(2020, 3)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

#=
This file contains some dictionary and functions to build
the dictionaries `_dict_unary_calls` and `_dict_binary_calls`
which allow to call the internal mutating functions. This
may be used to improve memory usage, e.g., to construct the
jet-coefficients used to integrate ODEs.
=#

"""
`_InternalMutFuncs`

Contains parameters and expressions that allow a simple
programatic construction for calling the internal mutating
functions.
"""
struct _InternalMutFuncs
    namef    :: Symbol   # internal name of the function
    argsf    :: NTuple   # arguments
    defexpr  :: Expr     # defining expr
    auxexpr  :: Expr     # auxiliary expr
end

# Constructor
function _InternalMutFuncs( namef::Vector )
    if length(namef) == 3
    	return _InternalMutFuncs( namef[1], namef[2], namef[3], Expr(:nothing) )
    else
        return _InternalMutFuncs( namef...)
    end
end

"""
`_dict_binary_ops`

`Dict{Symbol, Array{Any,1}}` with the information to
construct the `_InternalMutFuncs` related to binary
operations.

The keys correspond to the function symbols.

The arguments of the array are the function name (e.g. `add!`), a tuple
with the function arguments, and an `Expr` with the calling pattern. The
convention for the arguments of the functions and the calling pattern
is to use `:_res` for the (mutated) result, `:_arg1` and `_arg2`
for the required arguments, and `:_k` for the computed order
of `:_res`.

"""
const _dict_binary_ops = Dict(
    :+ => [:add!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 + _arg2)],
    :- => [:subst!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 - _arg2)],
    :* => [:mul!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 * _arg2)],
    :/ => [:div!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 / _arg2)],
    :^ => [:pow!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 ^ _arg2)],
);

"""
`_dict_binary_ops`

`Dict{Symbol, Array{Any,1}}` with the information to
construct the `_InternalMutFuncs` related to unary
operations.

The keys correspond to the function symbols.

The arguments of the array are the function name (e.g. `add!`), a tuple
with the function arguments, and an `Expr` with the calling pattern. The
convention for the arguments of the functions and the calling pattern
is to use `:_res` for the (mutated) result, `:_arg1`, for the required
argument, possibly `:_aux` when there is an auxiliary expression
needed, and `:_k` for the computed order of `:_res`. When an auxiliary
expression is required, and `Expr` defining its calling pattern is
added as the last entry of the vector.

"""
const _dict_unary_ops = Dict(
    :+ => [:add!,   (:_res, :_arg1, :_k), :(_res = + _arg1)],
    :- => [:subst!, (:_res, :_arg1, :_k), :(_res = - _arg1)],
    :sqr =>  [:sqr!, (:_res, :_arg1, :_k), :(_res = sqr(_arg1))],
    :sqrt => [:sqrt!, (:_res, :_arg1, :_k), :(_res = sqrt(_arg1))],
    :exp =>  [:exp!, (:_res, :_arg1, :_k), :(_res = exp(_arg1))],
    :log =>  [:log!, (:_res, :_arg1, :_k), :(_res = log(_arg1))],
    :identity => [:identity!, (:_res, :_arg1, :_k), :(_res = identity(_arg1))],
    :zero => [:zero!, (:_res, :_arg1, :_k), :(_res = zero(_arg1))],
    :one => [:one!, (:_res, :_arg1, :_k), :(_res = one(_arg1))],
    :abs => [:abs!, (:_res, :_arg1, :_k), :(_res = abs(_arg1))],
    :abs2 => [:abs2!, (:_res, :_arg1, :_k), :(_res = abs2(_arg1))],
    :deg2rad => [:deg2rad!, (:_res, :_arg1, :_k), :(_res = deg2rad(_arg1))],
    :rad2deg => [:rad2deg!, (:_res, :_arg1, :_k), :(_res = rad2deg(_arg1))],
    #
    :sin =>  [:sincos!, (:_res, :_aux, :_arg1, :_k), :(_res = sin(_arg1)),
        :(_aux = cos(_arg1))],
    :cos => [:sincos!, (:_aux, :_res, :_arg1, :_k), :(_res = cos(_arg1)),
        :(_aux = sin(_arg1))],
    :tan => [:tan!, (:_res, :_arg1, :_aux, :_k), :(_res = tan(_arg1)),
        :(_aux = tan(_arg1)^2)],
    :asin => [:asin!, (:_res, :_arg1, :_aux, :_k), :(_res = asin(_arg1)),
        :(_aux = sqrt(1 - _arg1^2))],
    :acos => [:acos!, (:_res, :_arg1, :_aux, :_k), :(_res = acos(_arg1)),
        :(_aux = sqrt(1 - _arg1^2))],
    :atan => [:atan!, (:_res, :_arg1, :_aux, :_k), :(_res = atan(_arg1)),
        :(_aux = 1 + _arg1^2)],
    :sinh => [:sinhcosh!, (:_res, :_aux, :_arg1, :_k), :(_res = sinh(_arg1)),
        :(_aux = cosh(_arg1))],
    :cosh => [:sinhcosh!, (:_aux, :_res, :_arg1, :_k), :(_res = cosh(_arg1)),
        :(_aux = sinh(_arg1))],
    :tanh => [:tanh!, (:_res, :_arg1, :_aux, :_k), :(_res = tanh(_arg1)),
        :(_aux = tanh(_arg1)^2)],
);

"""
```
_internalmutfunc_call( fn :: _InternalMutFuncs )
```

Creates the appropriate call to the internal mutating
function defined by the `_InternalMutFuncs` object.
This is used to construct [`_dict_unary_calls`](@ref)
and [`_dict_binary_calls`](@ref).
The call contains the prefix `TaylorSeries.`.
"""
_internalmutfunc_call( fn :: _InternalMutFuncs ) = (
    Expr( :call, Meta.parse("TaylorSeries.$(fn.namef)"), fn.argsf... ), fn.defexpr, fn.auxexpr )

"""
`_dict_unary_calls::Dict{Symbol, NTuple{2,Expr}}`

Dictionary with the expressions that define the
internal unary functions and the auxiliary functions,
whenever they exist. The keys correspond to those
functions, passed as symbols, with the defined
internal mutating functions.

Evaluating the entries generates expressions that represent
the actual calls to the internal mutating functions.
"""
const _dict_unary_calls = Dict{Symbol, NTuple{3,Expr}}()

#Populates the constant vector `_dict_unary_calls`.
for kk in keys(_dict_unary_ops)
    res = _internalmutfunc_call( _InternalMutFuncs(_dict_unary_ops[kk]) )
    push!(_dict_unary_calls, kk => res )
end

"""
`_dict_binary_calls::Dict{Symbol, NTuple{2,Expr}}`

Dictionary with the expressions that define the
internal binary functions and the auxiliary functions,
whenever they exist. The keys correspond to those
functions, passed as symbols, with the defined
internal mutating functions.

Evaluating the entries generates symbols that represent
the actual calls to the internal mutating functions.
"""
const _dict_binary_calls = Dict{Symbol, NTuple{3,Expr}}()
#Populates the constant vector `_dict_binary_calls`.
for kk in keys(_dict_binary_ops)
    res = _internalmutfunc_call( _InternalMutFuncs(_dict_binary_ops[kk]) )
    push!(_dict_binary_calls, kk => res )
end


text:
using ADBCVUMPS
using ADBCVUMPS:num_grad
using BCVUMPS:model_tensor,qrpos,lqpos,Ising,Ising22
using BCVUMPS:leftorth,leftenv,rightorth,rightenv,FLmap,ACenv,Cenv,ACCtoALAR,obs2x2FL,obs2x2FR
using Test
using Zygote
using ChainRulesTestUtils
using ChainRulesCore
using Random
using OMEinsum
using LinearAlgebra
using Random

@testset "matrix autodiff" begin
    a = randn(10, 10)
    @test Zygote.gradient(norm, a)[1] ≈ num_grad(norm, a)

    function foo1(x)
        sum(Float64[x 2x; 3x x])
    end
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1)

    # example to solve differential of array of array
    # use `[]` list then reshape
    A = Array{Array,2}(undef, 2, 2)
    for j = 1:2,i = 1:2
        A[i,j] = rand(2,2)
    end
    function foo2(x)
        # B[i,j] = A[i,j].*x   # mistake
        B = reshape([A[i].*x for i=1:4],2,2)
        return sum(sum(B))
    end
    @test Zygote.gradient(foo2, 1)[1] ≈ num_grad(foo2, 1)
end

@testset "$(Ni)x$(Nj) model_tensor" for Ni in [1,2,3], Nj in [1,2,3]
    Random.seed!(100)
    function foo1(β)
        M = model_tensor(Ising(Ni, Nj), β)
        return norm(norm(M))
    end
    @test isapprox(Zygote.gradient(foo1, 0.4)[1], num_grad(foo1, 0.4), atol=1e-8)

    function foo2(β)
        M = model_tensor(Ising22(Ni*Nj*0.1), β)
        return norm(norm(M))
    end
    @test isapprox(Zygote.gradient(foo2, 1)[1], num_grad(foo2, 1), atol=1e-8)
end

@testset "QR factorization" begin
    M = rand(10,10)
    function foo5(x)
        A = M.*x
        Q, R = qrpos(A)
        return norm(Q) + norm(R)
    end
    @test isapprox(Zygote.gradient(foo5, 1)[1], num_grad(foo5, 1), atol = 1e-5)
end

@testset "LQ factorization" begin
    M = rand(10,10)
    function foo6(x)
        A = M .*x
        L, Q = lqpos(A)
        return norm(L) + norm(Q)
    end
    @test isapprox(Zygote.gradient(foo6, 1)[1], num_grad(foo6, 1), atol = 1e-5)
end

@testset "$(Ni)x$(Nj) leftenv and rightenv" for Ni in [1,2,3], Nj in [1,2,3]
    Random.seed!(50)
    D, d = 3, 2
    A = Array{Array,2}(undef, Ni, Nj)
    S = Array{Array,2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = rand(D, d, D)
        S[i,j] = rand(D, d, D, D, d, D)
    end

    AL, = leftorth(A)
    _, AR, = rightorth(A)

    function foo3(β)
        M = model_tensor(Ising(Ni, Nj), β)
        λL, FL = leftenv(AL, M)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γcη,ηcγαaβ,βaα -> "(FL[i,j], S[i,j], FL[i,j])[] / ein"γcη,ηcγ -> "(FL[i,j], FL[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo3, 1)[1], num_grad(foo3, 1), atol=1e-8)

    function foo4(β)
        M = model_tensor(Ising(Ni, Nj), β)
        λR, FR = rightenv(AR, M)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γcη,ηcγαaβ,βaα -> "(FR[i,j], S[i,j], FR[i,j])[] / ein"γcη,ηcγ -> "(FR[i,j], FR[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo4, 1)[1], num_grad(foo4, 1), atol=1e-8)
end

@testset "$(Ni)x$(Nj) ACenv and Cenv" for Ni in [1,2,3], Nj in [1,2,3]
    Random.seed!(50)
    D, d = 3, 2
    A = Array{Array,2}(undef, Ni, Nj)
    S1 = Array{Array,2}(undef, Ni, Nj)
    S2 = Array{Array,2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = rand(D, d, D)
        S1[i,j] = rand(D, d, D, D, d, D)
        S2[i,j] = rand(D, D, D, D)
    end

    AL, L, _ = leftorth(A)
    R, AR, _ = rightorth(A)
    M = model_tensor(Ising(Ni, Nj), 1.0)
    _, FL = leftenv(AL, M)
    _, FR = rightenv(AR, M)

    C = Array{Array,2}(undef, Ni, Nj)
    AC = Array{Array,2}(undef, Ni, Nj)
    for j in 1:Nj,i in 1:Ni
        jr = j + 1 - (j + 1 > Nj) * Nj
        C[i,j] = L[i,j] * R[i,jr]
        AC[i,j] = ein"asc,cb -> asb"(AL[i,j], C[i,j])
    end

    function foo1(β)
        M = model_tensor(Ising(Ni, Nj), β)
        λAC, AC = ACenv(AC, FL, M, FR)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γcη,ηcγαaβ,βaα -> "(AC[i,j], S1[i,j], AC[i,j])[] / ein"γcη,γcη -> "(AC[i,j], AC[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo1, 1)[1], num_grad(foo1, 1), atol=1e-8)

    function foo2(β)
        M = model_tensor(Ising(Ni, Nj), β)
        λL, FL = leftenv(AL, M)
        λR, FR = rightenv(AR, M)
        λC, C = Cenv(C, FL, FR)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γη,ηγαβ,βα -> "(C[i,j],S2[i,j],C[i,j])[]/ein"γη,γη -> "(C[i,j],C[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo2, 1)[1], num_grad(foo2, 1), atol=1e-8)
end

@testset "$(Ni)x$(Nj) ACCtoALAR" for Ni in [1,2,3], Nj in [1,2,3]
    Random.seed!(10)
    D, d = 3, 2
    A = Array{Array,2}(undef, Ni, Nj)
    S1 = Array{Array,2}(undef, Ni, Nj)
    S2 = Array{Array,2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = rand(D, d, D)
        S1[i,j] = rand(D, d, D, D, d, D)
        S2[i,j] = rand(D, D, D, D)
    end

    AL, L, _ = leftorth(A)
    R, AR, _ = rightorth(A)
    M = model_tensor(Ising(Ni, Nj), 1.0)
    _, FL = leftenv(AL, M)
    _, FR = rightenv(AR, M)

    C = Array{Array,2}(undef, Ni, Nj)
    for j in 1:Nj,i in 1:Ni
        jr = j + 1 - (j + 1 > Nj) * Nj
        C[i,j] = L[i,j] * R[i,jr]
    end

    function foo1(β)
        M = model_tensor(Ising(Ni, Nj), β)
        AL, C, AR = ACCtoALAR(AL, C, AR, M, FL, FR)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γcη,ηcγαaβ,βaα -> "(AL[i,j], S1[i,j], AL[i,j])[] / ein"γcη,γcη -> "(AL[i,j], AL[i,j])[]
            s += ein"γη,ηγαβ,βα -> "(C[i,j],S2[i,j],C[i,j])[]/ein"γη,γη -> "(C[i,j],C[i,j])[]
            s += ein"γcη,ηcγαaβ,βaα -> "(AR[i,j], S1[i,j], AR[i,j])[] / ein"γcη,γcη -> "(AR[i,j], AR[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo1, 1)[1], num_grad(foo1, 1), atol=1e-4)
end

@testset "obs2x2 leftenv and rightenv" for Ni in [2], Nj in [2]
    Random.seed!(50)
    D, d = 3, 2
    A = Array{Array,2}(undef, Ni, Nj)
    S = Array{Array,2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        A[i,j] = rand(D, d, D)
        S[i,j] = rand(D, d, D, D, d, D)
    end

    AL, = leftorth(A)
    _, AR, = rightorth(A)

    function foo3(β)
        M = model_tensor(Ising(Ni, Nj), β)
        λL, FL = obs2x2FL(AL, M)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γcη,ηcγαaβ,βaα -> "(FL[i,j], S[i,j], FL[i,j])[] / ein"γcη,ηcγ -> "(FL[i,j], FL[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo3, 1)[1], num_grad(foo3, 1), atol=1e-8)

    function foo4(β)
        M = model_tensor(Ising(Ni, Nj), β)
        λR, FR = obs2x2FL(AR, M)
        s = 0
        for j in 1:Nj, i in 1:Ni
            s += ein"γcη,ηcγαaβ,βaα -> "(FR[i,j], S[i,j], FR[i,j])[] / ein"γcη,ηcγ -> "(FR[i,j], FR[i,j])[]
        end
        return s
    end
    @test isapprox(Zygote.gradient(foo4, 1)[1], num_grad(foo4, 1), atol=1e-8)
end


text:
using ProbNumDiffEq
using Test

du0 = [0.0]
u0 = [2.0]
tspan = (0.0, 6.3)
p = [1e1]

function vanderpol!(ddu, du, u, p, t)
    μ = p[1]
    @. ddu = μ * ((1 - u^2) * du - u)
end
prob_iip = SecondOrderODEProblem(vanderpol!, du0, u0, tspan, p)

function vanderpol(du, u, p, t)
    μ = p[1]
    ddu = μ .* ((1 .- u .^ 2) .* du .- u)
    return ddu
end
prob_oop = SecondOrderODEProblem(vanderpol, du0, u0, tspan, p)

appxsol = solve(prob_iip, Tsit5(), abstol=1e-7, reltol=1e-7)

@testset "IIP" begin
    for alg in (EK0(), EK1())
        @testset "$alg" begin
            @test solve(prob_iip, alg) isa ProbNumDiffEq.ProbODESolution
            @test solve(prob_iip, alg).u[end] ≈ appxsol.u[end] rtol = 1e-3
        end
    end
end

@testset "OOP" begin
    for alg in (EK0(), EK1())
        @testset "$alg" begin
            @test solve(prob_oop, alg) isa ProbNumDiffEq.ProbODESolution
            @test solve(prob_oop, alg).u[end] ≈ appxsol.u[end] rtol = 1e-3
        end
    end
end

@testset "RungeKuttaInit for SecondOrderODEProblems" begin
    @test_broken solve(prob_iip, EK1(initialization=RungeKuttaInit())) isa
                 ProbNumDiffEq.ProbODESolution
end

@testset "Fixed Diffusion" begin
    @test solve(prob_iip, EK0(diffusionmodel=:fixed)) isa ProbNumDiffEq.ProbODESolution
end