text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Profiling support, main entry point is the [`@profile`](@ref) macro.
"""
module Profile

import Base.StackTraces: lookup, UNKNOWN, show_spec_linfo, StackFrame

# deprecated functions: use `getdict` instead
lookup(ip::UInt) = lookup(convert(Ptr{Cvoid}, ip))

export @profile

"""
    @profile

`@profile <expression>` runs your expression while taking periodic backtraces. These are
appended to an internal buffer of backtraces.
"""
macro profile(ex)
    return quote
        try
            status = start_timer()
            if status < 0
                error(error_codes[status])
            end
            $(esc(ex))
        finally
            stop_timer()
        end
    end
end

####
#### User-level functions
####

"""
    init(; n::Integer, delay::Real))

Configure the `delay` between backtraces (measured in seconds), and the number `n` of
instruction pointers that may be stored per thread. Each instruction pointer corresponds to a single
line of code; backtraces generally consist of a long list of instruction pointers. Current
settings can be obtained by calling this function with no arguments, and each can be set
independently using keywords or in the order `(n, delay)`.

!!! compat "Julia 1.8"
    As of Julia 1.8, this function allocates space for `n` instruction pointers per thread being profiled.
    Previously this was `n` total.
"""
function init(; n::Union{Nothing,Integer} = nothing, delay::Union{Nothing,Real} = nothing)
    n_cur = ccall(:jl_profile_maxlen_data, Csize_t, ())
    delay_cur = ccall(:jl_profile_delay_nsec, UInt64, ())/10^9
    if n === nothing && delay === nothing
        return Int(n_cur), delay_cur
    end
    nnew = (n === nothing) ? n_cur : n
    delaynew = (delay === nothing) ? delay_cur : delay
    init(nnew, delaynew)
end

function init(n::Integer, delay::Real)
    nthreads = Sys.iswindows() ? 1 : Threads.nthreads() # windows only profiles the main thread
    sample_size_bytes = sizeof(Ptr) # == Sys.WORD_SIZE / 8
    buffer_samples = n * nthreads
    buffer_size_bytes = buffer_samples * sample_size_bytes
    if buffer_size_bytes > 2^29 && Sys.WORD_SIZE == 32
        buffer_size_bytes_per_thread = floor(Int, 2^29 / nthreads)
        buffer_samples_per_thread = floor(Int, buffer_size_bytes_per_thread / sample_size_bytes)
        buffer_samples = buffer_samples_per_thread * nthreads
        buffer_size_bytes = buffer_samples * sample_size_bytes
        @warn "Requested profile buffer limited to 512MB (n = $buffer_samples_per_thread per thread) given that this system is 32-bit"
    end
    status = ccall(:jl_profile_init, Cint, (Csize_t, UInt64), buffer_samples, round(UInt64,10^9*delay))
    if status == -1
        error("could not allocate space for ", n, " instruction pointers per thread being profiled ($nthreads threads, $(Base.format_bytes(buffer_size_bytes)) total)")
    end
end

# init with default values
# Use a max size of 10M profile samples, and fire timer every 1ms
# (that should typically give around 100 seconds of record)
if Sys.iswindows() && Sys.WORD_SIZE == 32
    # The Win32 unwinder is 1000x slower than elsewhere (around 1ms/frame),
    # so we don't want to slow the program down by quite that much
    __init__() = init(1_000_000, 0.01)
else
    __init__() = init(10_000_000, 0.001)
end

"""
    clear()

Clear any existing backtraces from the internal buffer.
"""
clear() = ccall(:jl_profile_clear_data, Cvoid, ())

const LineInfoDict = Dict{UInt64, Vector{StackFrame}}
const LineInfoFlatDict = Dict{UInt64, StackFrame}

struct ProfileFormat
    maxdepth::Int
    mincount::Int
    noisefloor::Float64
    sortedby::Symbol
    combine::Bool
    C::Bool
    recur::Symbol
    function ProfileFormat(;
        C = false,
        combine = true,
        maxdepth::Int = typemax(Int),
        mincount::Int = 0,
        noisefloor = 0,
        sortedby::Symbol = :filefuncline,
        recur::Symbol = :off)
        return new(maxdepth, mincount, noisefloor, sortedby, combine, C, recur)
    end
end

"""
    print([io::IO = stdout,] [data::Vector]; kwargs...)

Prints profiling results to `io` (by default, `stdout`). If you do not
supply a `data` vector, the internal buffer of accumulated backtraces
will be used.

The keyword arguments can be any combination of:

 - `format` -- Determines whether backtraces are printed with (default, `:tree`) or without (`:flat`)
   indentation indicating tree structure.

 - `C` -- If `true`, backtraces from C and Fortran code are shown (normally they are excluded).

 - `combine` -- If `true` (default), instruction pointers are merged that correspond to the same line of code.

 - `maxdepth` -- Limits the depth higher than `maxdepth` in the `:tree` format.

 - `sortedby` -- Controls the order in `:flat` format. `:filefuncline` (default) sorts by the source
    line, `:count` sorts in order of number of collected samples, and `:overhead` sorts by the number of samples
    incurred by each function by itself.

 - `noisefloor` -- Limits frames that exceed the heuristic noise floor of the sample (only applies to format `:tree`).
    A suggested value to try for this is 2.0 (the default is 0). This parameter hides samples for which `n <= noisefloor * √N`,
    where `n` is the number of samples on this line, and `N` is the number of samples for the callee.

 - `mincount` -- Limits the printout to only those lines with at least `mincount` occurrences.

 - `recur` -- Controls the recursion handling in `:tree` format. `:off` (default) prints the tree as normal. `:flat` instead
    compresses any recursion (by ip), showing the approximate effect of converting any self-recursion into an iterator.
    `:flatc` does the same but also includes collapsing of C frames (may do odd things around `jl_apply`).
"""
function print(io::IO,
        data::Vector{<:Unsigned} = fetch(),
        lidict::Union{LineInfoDict, LineInfoFlatDict} = getdict(data)
        ;
        format = :tree,
        C = false,
        combine = true,
        maxdepth::Int = typemax(Int),
        mincount::Int = 0,
        noisefloor = 0,
        sortedby::Symbol = :filefuncline,
        recur::Symbol = :off)
    print(io, data, lidict, ProfileFormat(
            C = C,
            combine = combine,
            maxdepth = maxdepth,
            mincount = mincount,
            noisefloor = noisefloor,
            sortedby = sortedby,
            recur = recur),
        format)
end

function print(io::IO, data::Vector{<:Unsigned}, lidict::Union{LineInfoDict, LineInfoFlatDict}, fmt::ProfileFormat, format::Symbol)
    cols::Int = Base.displaysize(io)[2]
    data = convert(Vector{UInt64}, data)
    fmt.recur ∈ (:off, :flat, :flatc) || throw(ArgumentError("recur value not recognized"))
    if format === :tree
        tree(io, data, lidict, cols, fmt)
    elseif format === :flat
        fmt.recur === :off || throw(ArgumentError("format flat only implements recur=:off"))
        flat(io, data, lidict, cols, fmt)
    else
        throw(ArgumentError("output format $(repr(format)) not recognized"))
    end
end

"""
    print([io::IO = stdout,] data::Vector, lidict::LineInfoDict; kwargs...)

Prints profiling results to `io`. This variant is used to examine results exported by a
previous call to [`retrieve`](@ref). Supply the vector `data` of backtraces and
a dictionary `lidict` of line information.

See `Profile.print([io], data)` for an explanation of the valid keyword arguments.
"""
print(data::Vector{<:Unsigned} = fetch(), lidict::Union{LineInfoDict, LineInfoFlatDict} = getdict(data); kwargs...) =
    print(stdout, data, lidict; kwargs...)

"""
    retrieve() -> data, lidict

"Exports" profiling results in a portable format, returning the set of all backtraces
(`data`) and a dictionary that maps the (session-specific) instruction pointers in `data` to
`LineInfo` values that store the file name, function name, and line number. This function
allows you to save profiling results for future analysis.
"""
function retrieve()
    data = fetch()
    return (data, getdict(data))
end

function getdict(data::Vector{UInt})
    dict = LineInfoDict()
    return getdict!(dict, data)
end

function getdict!(dict::LineInfoDict, data::Vector{UInt})
    for ip in data
        # Lookup is expensive, so do it only once per ip.
        haskey(dict, UInt64(ip)) && continue
        st = lookup(convert(Ptr{Cvoid}, ip))
        # To correct line numbers for moving code, put it in the form expected by
        # Base.update_stackframes_callback[]
        stn = map(x->(x, 1), st)
        try Base.invokelatest(Base.update_stackframes_callback[], stn) catch end
        dict[UInt64(ip)] = map(first, stn)
    end
    return dict
end

"""
    flatten(btdata::Vector, lidict::LineInfoDict) -> (newdata::Vector{UInt64}, newdict::LineInfoFlatDict)

Produces "flattened" backtrace data. Individual instruction pointers
sometimes correspond to a multi-frame backtrace due to inlining; in
such cases, this function inserts fake instruction pointers for the
inlined calls, and returns a dictionary that is a 1-to-1 mapping
between instruction pointers and a single StackFrame.
"""
function flatten(data::Vector, lidict::LineInfoDict)
    # Makes fake instruction pointers, counting down from typemax(UInt)
    newip = typemax(UInt64) - 1
    taken = Set(keys(lidict))  # make sure we don't pick one that's already used
    newdict = Dict{UInt64,StackFrame}()
    newmap  = Dict{UInt64,Vector{UInt64}}()
    for (ip, trace) in lidict
        if length(trace) == 1
            newdict[ip] = trace[1]
        else
            newm = UInt64[]
            for sf in trace
                while newip ∈ taken && newip > 0
                    newip -= 1
                end
                newip == 0 && error("all possible instruction pointers used")
                push!(newm, newip)
                newdict[newip] = sf
                newip -= 1
            end
            newmap[ip] = newm
        end
    end
    newdata = UInt64[]
    for ip::UInt64 in data
        if haskey(newmap, ip)
            append!(newdata, newmap[ip])
        else
            push!(newdata, ip)
        end
    end
    return (newdata, newdict)
end

# Take a file-system path and try to form a concise representation of it
# based on the package ecosystem
function short_path(spath::Symbol, filenamecache::Dict{Symbol, String})
    return get!(filenamecache, spath) do
        path = string(spath)
        if isabspath(path)
            if ispath(path)
                # try to replace the file-system prefix with a short "@Module" one,
                # assuming that profile came from the current machine
                # (or at least has the same file-system layout)
                root = path
                while !isempty(root)
                    root, base = splitdir(root)
                    isempty(base) && break
                    @assert startswith(path, root)
                    for proj in Base.project_names
                        project_file = joinpath(root, proj)
                        if Base.isfile_casesensitive(project_file)
                            pkgid = Base.project_file_name_uuid(project_file, "")
                            isempty(pkgid.name) && return path # bad Project file
                            # return the joined the module name prefix and path suffix
                            path = path[nextind(path, sizeof(root)):end]
                            return string("@", pkgid.name, path)
                        end
                    end
                end
            end
            return path
        elseif isfile(joinpath(Sys.BINDIR::String, Base.DATAROOTDIR, "julia", "base", path))
            # do the same mechanic for Base (or Core/Compiler) files as above,
            # but they start from a relative path
            return joinpath("@Base", normpath(path))
        else
            # for non-existent relative paths (such as "REPL[1]"), just consider simplifying them
            return normpath(path) # drop leading "./"
        end
    end
end

"""
    callers(funcname, [data, lidict], [filename=<filename>], [linerange=<start:stop>]) -> Vector{Tuple{count, lineinfo}}

Given a previous profiling run, determine who called a particular function. Supplying the
filename (and optionally, range of line numbers over which the function is defined) allows
you to disambiguate an overloaded method. The returned value is a vector containing a count
of the number of calls and line information about the caller. One can optionally supply
backtrace `data` obtained from [`retrieve`](@ref); otherwise, the current internal
profile buffer is used.
"""
function callers end

function callers(funcname::String, bt::Vector, lidict::LineInfoFlatDict; filename = nothing, linerange = nothing)
    if filename === nothing && linerange === nothing
        return callersf(li -> String(li.func) == funcname,
            bt, lidict)
    end
    filename === nothing && throw(ArgumentError("if supplying linerange, you must also supply the filename"))
    filename = String(filename)
    if linerange === nothing
        return callersf(li -> String(li.func) == funcname && String(li.file) == filename,
            bt, lidict)
    else
        return callersf(li -> String(li.func) == funcname && String(li.file) == filename && in(li.line, linerange),
            bt, lidict)
    end
end

callers(funcname::String, bt::Vector, lidict::LineInfoDict; kwargs...) =
    callers(funcname, flatten(bt, lidict)...; kwargs...)
callers(funcname::String; kwargs...) = callers(funcname, retrieve()...; kwargs...)
callers(func::Function, bt::Vector, lidict::LineInfoFlatDict; kwargs...) =
    callers(string(func), bt, lidict; kwargs...)
callers(func::Function; kwargs...) = callers(string(func), retrieve()...; kwargs...)

##
## For --track-allocation
##
# Reset the malloc log. Used to avoid counting memory allocated during
# compilation.

"""
    clear_malloc_data()

Clears any stored memory allocation data when running julia with `--track-allocation`.
Execute the command(s) you want to test (to force JIT-compilation), then call
[`clear_malloc_data`](@ref). Then execute your command(s) again, quit
Julia, and examine the resulting `*.mem` files.
"""
clear_malloc_data() = ccall(:jl_clear_malloc_data, Cvoid, ())

# C wrappers
start_timer() = ccall(:jl_profile_start_timer, Cint, ())

stop_timer() = ccall(:jl_profile_stop_timer, Cvoid, ())

is_running() = ccall(:jl_profile_is_running, Cint, ())!=0

is_buffer_full() = ccall(:jl_profile_is_buffer_full, Cint, ())!=0

get_data_pointer() = convert(Ptr{UInt}, ccall(:jl_profile_get_data, Ptr{UInt8}, ()))

len_data() = convert(Int, ccall(:jl_profile_len_data, Csize_t, ()))

maxlen_data() = convert(Int, ccall(:jl_profile_maxlen_data, Csize_t, ()))

error_codes = Dict(
    -1=>"cannot specify signal action for profiling",
    -2=>"cannot create the timer for profiling",
    -3=>"cannot start the timer for profiling",
    -4=>"cannot unblock SIGUSR1")

"""
    fetch() -> data

Returns a copy of the buffer of profile backtraces. Note that the
values in `data` have meaning only on this machine in the current session, because it
depends on the exact memory addresses used in JIT-compiling. This function is primarily for
internal use; [`retrieve`](@ref) may be a better choice for most users.
"""
function fetch()
    maxlen = maxlen_data()
    len = len_data()
    if is_buffer_full()
        @warn """The profile data buffer is full; profiling probably terminated
                 before your program finished. To profile for longer runs, call
                 `Profile.init()` with a larger buffer and/or larger delay."""
    end
    data = Vector{UInt}(undef, len)
    GC.@preserve data unsafe_copyto!(pointer(data), get_data_pointer(), len)
    return data
end

## Print as a flat list
# Counts the number of times each line appears, at any nesting level and at the topmost level
# Merging multiple equivalent entries and recursive calls
function parse_flat(::Type{T}, data::Vector{UInt64}, lidict::Union{LineInfoDict, LineInfoFlatDict}, C::Bool) where {T}
    lilist = StackFrame[]
    n = Int[]
    m = Int[]
    lilist_idx = Dict{T, Int}()
    recursive = Set{T}()
    first = true
    totalshots = 0
    for ip in data
        if ip == 0
            totalshots += 1
            empty!(recursive)
            first = true
            continue
        end
        frames = lidict[ip]
        nframes = (frames isa Vector ? length(frames) : 1)
        for i = 1:nframes
            frame = (frames isa Vector ? frames[i] : frames)
            !C && frame.from_c && continue
            key = (T === UInt64 ? ip : frame)
            idx = get!(lilist_idx, key, length(lilist) + 1)
            if idx > length(lilist)
                push!(recursive, key)
                push!(lilist, frame)
                push!(n, 1)
                push!(m, 0)
            elseif !(key in recursive)
                push!(recursive, key)
                n[idx] += 1
            end
            if first
                m[idx] += 1
                first = false
            end
        end
    end
    @assert length(lilist) == length(n) == length(m) == length(lilist_idx)
    return (lilist, n, m, totalshots)
end

function flat(io::IO, data::Vector{UInt64}, lidict::Union{LineInfoDict, LineInfoFlatDict}, cols::Int, fmt::ProfileFormat)
    lilist, n, m, totalshots = parse_flat(fmt.combine ? StackFrame : UInt64, data, lidict, fmt.C)
    if isempty(lilist)
        warning_empty()
        return
    end
    if false # optional: drop the "non-interpretable" ones
        keep = map(frame -> frame != UNKNOWN && frame.line != 0, lilist)
        lilist = lilist[keep]
        n = n[keep]
        m = m[keep]
    end
    filenamemap = Dict{Symbol,String}()
    print_flat(io, lilist, n, m, cols, filenamemap, fmt)
    Base.println(io, "Total snapshots: ", totalshots)
    nothing
end

function print_flat(io::IO, lilist::Vector{StackFrame},
        n::Vector{Int}, m::Vector{Int},
        cols::Int, filenamemap::Dict{Symbol,String},
        fmt::ProfileFormat)
    if fmt.sortedby === :count
        p = sortperm(n)
    elseif fmt.sortedby === :overhead
        p = sortperm(m)
    else
        p = liperm(lilist)
    end
    lilist = lilist[p]
    n = n[p]
    m = m[p]
    filenames = String[short_path(li.file, filenamemap) for li in lilist]
    funcnames = String[string(li.func) for li in lilist]
    wcounts = max(6, ndigits(maximum(n)))
    wself = max(9, ndigits(maximum(m)))
    maxline = 1
    maxfile = 6
    maxfunc = 10
    for i in 1:length(lilist)
        li = lilist[i]
        maxline = max(maxline, li.line)
        maxfunc = max(maxfunc, length(funcnames[i]))
        maxfile = max(maxfile, length(filenames[i]))
    end
    wline = max(5, ndigits(maxline))
    ntext = max(20, cols - wcounts - wself - wline - 3)
    maxfunc += 25 # for type signatures
    if maxfile + maxfunc <= ntext
        wfile = maxfile
        wfunc = ntext - maxfunc # take the full width (for type sig)
    else
        wfile = 2*ntext÷5
        wfunc = 3*ntext÷5
    end
    println(io, lpad("Count", wcounts, " "), " ", lpad("Overhead", wself, " "), " ",
            rpad("File", wfile, " "), " ", lpad("Line", wline, " "), " Function")
    println(io, lpad("=====", wcounts, " "), " ", lpad("========", wself, " "), " ",
            rpad("====", wfile, " "), " ", lpad("====", wline, " "), " ========")
    for i = 1:length(n)
        n[i] < fmt.mincount && continue
        li = lilist[i]
        Base.print(io, lpad(string(n[i]), wcounts, " "), " ")
        Base.print(io, lpad(string(m[i]), wself, " "), " ")
        if li == UNKNOWN
            if !fmt.combine && li.pointer != 0
                Base.print(io, "@0x", string(li.pointer, base=16))
            else
                Base.print(io, "[any unknown stackframes]")
            end
        else
            file = filenames[i]
            isempty(file) && (file = "[unknown file]")
            Base.print(io, rpad(rtruncto(file, wfile), wfile, " "), " ")
            Base.print(io, lpad(li.line > 0 ? string(li.line) : "?", wline, " "), " ")
            fname = funcnames[i]
            if !li.from_c && li.linfo !== nothing
                fname = sprint(show_spec_linfo, li)
            end
            isempty(fname) && (fname = "[unknown function]")
            Base.print(io, ltruncto(fname, wfunc))
        end
        println(io)
    end
    nothing
end

## A tree representation

# Representation of a prefix trie of backtrace counts
mutable struct StackFrameTree{T} # where T <: Union{UInt64, StackFrame}
    # content fields:
    frame::StackFrame
    count::Int          # number of frames this appeared in
    overhead::Int       # number frames where this was the code being executed
    flat_count::Int     # number of times this frame was in the flattened representation (unlike count, this'll sum to 100% of parent)
    max_recur::Int      # maximum number of times this frame was the *top* of the recursion in the stack
    count_recur::Int    # sum of the number of times this frame was the *top* of the recursion in a stack (divide by count to get an average)
    down::Dict{T, StackFrameTree{T}}
    # construction workers:
    recur::Int
    builder_key::Vector{UInt64}
    builder_value::Vector{StackFrameTree{T}}
    up::StackFrameTree{T}
    StackFrameTree{T}() where {T} = new(UNKNOWN, 0, 0, 0, 0, 0, Dict{T, StackFrameTree{T}}(), 0, UInt64[], StackFrameTree{T}[])
end

const indent_s = "    ╎"^10
const indent_z = collect(eachindex(indent_s))
function indent(depth::Int)
    depth < 1 && return ""
    depth <= length(indent_z) && return indent_s[1:indent_z[depth]]
    div, rem = divrem(depth, length(indent_z))
    indent = indent_s^div
    rem != 0 && (indent *= SubString(indent_s, 1, indent_z[rem]))
    return indent
end

function tree_format(frames::Vector{<:StackFrameTree}, level::Int, cols::Int, maxes, filenamemap::Dict{Symbol,String}, showpointer::Bool)
    nindent = min(cols>>1, level)
    ndigoverhead = ndigits(maxes.overhead)
    ndigcounts = ndigits(maxes.count)
    ndigline = ndigits(maximum(frame.frame.line for frame in frames)) + 6
    ntext = max(30, cols - ndigoverhead - nindent - ndigcounts - ndigline - 6)
    widthfile = 2*ntext÷5 # min 12
    widthfunc = 3*ntext÷5 # min 18
    strs = Vector{String}(undef, length(frames))
    showextra = false
    if level > nindent
        nextra = level - nindent
        nindent -= ndigits(nextra) + 2
        showextra = true
    end
    for i = 1:length(frames)
        frame = frames[i]
        li = frame.frame
        stroverhead = lpad(frame.overhead > 0 ? string(frame.overhead) : "", ndigoverhead, " ")
        base = nindent == 0 ? "" : indent(nindent - 1) * " "
        if showextra
            base = string(base, "+", nextra, " ")
        end
        strcount = rpad(string(frame.count), ndigcounts, " ")
        if li != UNKNOWN
            if li.line == li.pointer
                strs[i] = string(stroverhead, "╎", base, strcount, " ",
                    "[unknown function] (pointer: 0x",
                    string(li.pointer, base = 16, pad = 2*sizeof(Ptr{Cvoid})),
                    ")")
            else
                if !li.from_c && li.linfo !== nothing
                    fname = sprint(show_spec_linfo, li)
                else
                    fname = string(li.func)
                end
                filename = short_path(li.file, filenamemap)
                if showpointer
                    fname = string(
                        "0x",
                        string(li.pointer, base = 16, pad = 2*sizeof(Ptr{Cvoid})),
                        " ",
                        fname)
                end
                strs[i] = string(stroverhead, "╎", base, strcount, " ",
                    rtruncto(filename, widthfile),
                    ":",
                    li.line == -1 ? "?" : string(li.line),
                    "; ",
                    ltruncto(fname, widthfunc))
            end
        else
            strs[i] = string(stroverhead, "╎", base, strcount, " [unknown stackframe]")
        end
    end
    return strs
end

# turn a list of backtraces into a tree (implicitly separated by NULL markers)
function tree!(root::StackFrameTree{T}, all::Vector{UInt64}, lidict::Union{LineInfoFlatDict, LineInfoDict}, C::Bool, recur::Symbol) where {T}
    parent = root
    tops = Vector{StackFrameTree{T}}()
    build = Vector{StackFrameTree{T}}()
    startframe = length(all)
    for i in startframe:-1:1
        ip = all[i]
        if ip == 0
            # sentinel value indicates the start of a new backtrace
            empty!(build)
            root.recur = 0
            if recur !== :off
                # We mark all visited nodes to so we'll only count those branches
                # once for each backtrace. Reset that now for the next backtrace.
                push!(tops, parent)
                for top in tops
                    while top.recur != 0
                        top.max_recur < top.recur && (top.max_recur = top.recur)
                        top.recur = 0
                        top = top.up
                    end
                end
                empty!(tops)
            end
            let this = parent
                while this !== root
                    this.flat_count += 1
                    this = this.up
                end
            end
            parent.overhead += 1
            parent = root
            root.count += 1
            startframe = i
        else
            pushfirst!(build, parent)
            if recur === :flat || recur === :flatc
                # Rewind the `parent` tree back, if this exact ip was already present *higher* in the current tree
                found = false
                for j in 1:(startframe - i)
                    if ip == all[i + j]
                        if recur === :flat # if not flattening C frames, check that now
                            frames = lidict[ip]
                            frame = (frames isa Vector ? frames[1] : frames)
                            frame.from_c && break # not flattening this frame
                        end
                        push!(tops, parent)
                        parent = build[j]
                        parent.recur += 1
                        parent.count_recur += 1
                        found = true
                        break
                    end
                end
                found && continue
            end
            builder_key = parent.builder_key
            builder_value = parent.builder_value
            fastkey = searchsortedfirst(builder_key, ip)
            if fastkey < length(builder_key) && builder_key[fastkey] === ip
                # jump forward to the end of the inlining chain
                # avoiding an extra (slow) lookup of `ip` in `lidict`
                # and an extra chain of them in `down`
                # note that we may even have this === parent (if we're ignoring this frame ip)
                this = builder_value[fastkey]
                let this = this
                    while this !== parent && (recur === :off || this.recur == 0)
                        this.count += 1
                        this.recur = 1
                        this = this.up
                    end
                end
                parent = this
                continue
            end
            frames = lidict[ip]
            nframes = (frames isa Vector ? length(frames) : 1)
            this = parent
            # add all the inlining frames
            for i = nframes:-1:1
                frame = (frames isa Vector ? frames[i] : frames)
                !C && frame.from_c && continue
                key = (T === UInt64 ? ip : frame)
                this = get!(StackFrameTree{T}, parent.down, key)
                if recur === :off || this.recur == 0
                    this.frame = frame
                    this.up = parent
                    this.count += 1
                    this.recur = 1
                end
                parent = this
            end
            # record where the end of this chain is for this ip
            insert!(builder_key, fastkey, ip)
            insert!(builder_value, fastkey, this)
        end
    end
    function cleanup!(node::StackFrameTree)
        stack = [node]
        while !isempty(stack)
            node = pop!(stack)
            node.recur = 0
            empty!(node.builder_key)
            empty!(node.builder_value)
            append!(stack, values(node.down))
        end
        nothing
    end
    cleanup!(root)
    return root
end

function maxstats(root::StackFrameTree)
    maxcount = Ref(0)
    maxflatcount = Ref(0)
    maxoverhead = Ref(0)
    maxmaxrecur = Ref(0)
    stack = [root]
    while !isempty(stack)
        node = pop!(stack)
        maxcount[] = max(maxcount[], node.count)
        maxoverhead[] = max(maxoverhead[], node.overhead)
        maxflatcount[] = max(maxflatcount[], node.flat_count)
        maxmaxrecur[] = max(maxmaxrecur[], node.max_recur)
        append!(stack, values(node.down))
    end
    return (count=maxcount[], count_flat=maxflatcount[], overhead=maxoverhead[], max_recur=maxmaxrecur[])
end

# Print the stack frame tree starting at a particular root. Uses a worklist to
# avoid stack overflows.
function print_tree(io::IO, bt::StackFrameTree{T}, cols::Int, fmt::ProfileFormat) where T
    maxes = maxstats(bt)
    filenamemap = Dict{Symbol,String}()
    worklist = [(bt, 0, 0, "")]
    println(io, "Overhead ╎ [+additional indent] Count File:Line; Function")
    println(io, "=========================================================")
    while !isempty(worklist)
        (bt, level, noisefloor, str) = popfirst!(worklist)
        isempty(str) || println(io, str)
        level > fmt.maxdepth && continue
        isempty(bt.down) && continue
        # Order the line information
        nexts = collect(values(bt.down))
        # Generate the string for each line
        strs = tree_format(nexts, level, cols, maxes, filenamemap, T === UInt64)
        # Recurse to the next level
        if fmt.sortedby === :count
            counts = collect(frame.count for frame in nexts)
            p = sortperm(counts)
        elseif fmt.sortedby === :overhead
            m = collect(frame.overhead for frame in nexts)
            p = sortperm(m)
        elseif fmt.sortedby === :flat_count
            m = collect(frame.flat_count for frame in nexts)
            p = sortperm(m)
        else
            lilist = collect(frame.frame for frame in nexts)
            p = liperm(lilist)
        end
        for i in reverse(p)
            down = nexts[i]
            count = down.count
            count < fmt.mincount && continue
            count < noisefloor && continue
            str = strs[i]
            noisefloor_down = fmt.noisefloor > 0 ? floor(Int, fmt.noisefloor * sqrt(count)) : 0
            pushfirst!(worklist, (down, level + 1, noisefloor_down, str))
        end
    end
end

function tree(io::IO, data::Vector{UInt64}, lidict::Union{LineInfoFlatDict, LineInfoDict}, cols::Int, fmt::ProfileFormat)
    if fmt.combine
        root = tree!(StackFrameTree{StackFrame}(), data, lidict, fmt.C, fmt.recur)
    else
        root = tree!(StackFrameTree{UInt64}(), data, lidict, fmt.C, fmt.recur)
    end
    if isempty(root.down)
        warning_empty()
        return
    end
    print_tree(io, root, cols, fmt)
    Base.println(io, "Total snapshots: ", root.count)
    nothing
end

function callersf(matchfunc::Function, bt::Vector, lidict::LineInfoFlatDict)
    counts = Dict{StackFrame, Int}()
    lastmatched = false
    for id in bt
        if id == 0
            lastmatched = false
            continue
        end
        li = lidict[id]
        if lastmatched
            if haskey(counts, li)
                counts[li] += 1
            else
                counts[li] = 1
            end
        end
        lastmatched = matchfunc(li)
    end
    k = collect(keys(counts))
    v = collect(values(counts))
    p = sortperm(v, rev=true)
    return [(v[i], k[i]) for i in p]
end

# Utilities
function rtruncto(str::String, w::Int)
    if length(str) <= w
        return str
    else
        return string("...", str[prevind(str, end, w-4):end])
    end
end
function ltruncto(str::String, w::Int)
    if length(str) <= w
        return str
    else
        return string(str[1:nextind(str, 1, w-4)], "...")
    end
end

truncto(str::Symbol, w::Int) = truncto(string(str), w)

# Order alphabetically (file, function) and then by line number
function liperm(lilist::Vector{StackFrame})
    function lt(a::StackFrame, b::StackFrame)
        a == UNKNOWN && return false
        b == UNKNOWN && return true
        fcmp = cmp(a.file, b.file)
        fcmp < 0 && return true
        fcmp > 0 && return false
        fcmp = cmp(a.func, b.func)
        fcmp < 0 && return true
        fcmp > 0 && return false
        fcmp = cmp(a.line, b.line)
        fcmp < 0 && return true
        return false
    end
    return sortperm(lilist, lt = lt)
end

warning_empty() = @warn """
            There were no samples collected. Run your program longer (perhaps by
            running it multiple times), or adjust the delay between samples with
            `Profile.init()`."""

end # module

max_stars_repo_name: tkf/julia

max_issues_repo_name: tkf/julia

max_forks_repo_name: tkf/julia


text:
# MArray.jl

mutable struct Variables{N,T} <: TupleVector{N,T}
    v::NTuple{N,T}
    Variables{N,T}(x::NTuple{N,T}) where {N,T} = new{N,T}(x)
    Variables{N,T}(x::NTuple{N,Any}) where {N,T} = new{N,T}(convert_ntuple(T, x))
    Variables{N,T}(::UndefInitializer) where {N,T} = new{N,T}()
end

@inline Variables(a::TupleVector{N}) where N = Variables{N}(Tuple(a))
@generated function (::Type{Variables{N,T}})(x::Tuple) where {N,T}
    return quote
        $(Expr(:meta, :inline))
        Variables{N,T}(x)
    end
end
@generated function (::Type{Variables{N}})(x::T) where {N,T<:Tuple}
    return quote
        $(Expr(:meta, :inline))
        Variables{N,promote_tuple_eltype(T)}(x)
    end
end

@propagate_inbounds function Base.getindex(v::Variables, i::Int)
    @boundscheck checkbounds(v,i)
    T = eltype(v)
    if isbitstype(T)
        return GC.@preserve v unsafe_load(Base.unsafe_convert(Ptr{T}, pointer_from_objref(v)), i)
    end
    v.v[i]
end
@propagate_inbounds function Base.setindex!(v::Variables, val, i::Int)
    @boundscheck checkbounds(v,i)
    T = eltype(v)
    if isbitstype(T)
        GC.@preserve v unsafe_store!(Base.unsafe_convert(Ptr{T}, pointer_from_objref(v)), convert(T, val), i)
    else
        # This one is unsafe (#27)
        # unsafe_store!(Base.unsafe_convert(Ptr{Ptr{Nothing}}, pointer_from_objref(v.data)), pointer_from_objref(val), i)
        error("setindex!() with non-isbitstype eltype is not supported by TupleVectors. Consider using FixedVector.")
    end
    return val
end

@inline Base.Tuple(v::Variables) = v.v
Base.dataids(ma::Variables) = (UInt(pointer(ma)),)

@inline function Base.unsafe_convert(::Type{Ptr{T}}, a::Variables{N,T}) where {N,T}
    Base.unsafe_convert(Ptr{T}, pointer_from_objref(a))
end

function Base.promote_rule(::Type{<:Variables{N,T}}, ::Type{<:Variables{N,U}}) where {N,T,U}
    Variables{N,promote_type(T,U)}
end

# MVector.jl

@inline Variables(x::NTuple{N,Any}) where N = Variables{N}(x)
@inline Variables{N}(x::NTuple{N,T}) where {N,T} = Variables{N,T}(x)
@inline Variables{N}(x::NTuple{N,Any}) where N = Variables{N, promote_tuple_eltype(typeof(x))}(x)

# Some more advanced constructor-like functions
@inline Base.zeros(::Type{Variables{N}}) where N = zeros(Variables{N,Float64})
@inline Base.ones(::Type{Variables{N}}) where N = ones(Variables{N,Float64})

ext: jl


text:
@doc """ Give an Array of pics file paths, return data
""" ->
function read_pics(fls::Array{UTF8String,1})
    imgs =  pmap(load, fls)
    dt = @parallel (vcat) for img in imgs
        data(reinterpret(UInt8, img))[:]'
    end

    convert(Array{Float32,2},dt)
end

@doc """ get train data
""" ->
function train_data()
    pics = readdir(trainresize_dir)
    pics_path = map(x->joinpath(trainresize_dir,x), pics)

    pics = map(y->convert(ASCIIString,y), map(x->x[1:end-4], pics))

    sample_label_dict = JLD.load(sample_label_fl,"sample_label_dict")
    labels = map(x->sample_label_dict[x], pics)
    read_pics(pics_path),labels
end

@doc """ get test data
""" ->
function test_data()
    pics = readdir(testresize_dir)
    pics_path = map(x->joinpath(testresize_dir,x), pics)

    pics,read_pics(pics_path)
end

@doc """ return sample_names=>labels dict
""" ->
function readlabel()
    dt = readcsv(train_lb,ASCIIString)
    dt =  dt[2:end,:]
    classes = sort(unique(dt[:,2]))
    sample_names = dt[:,1]
    num_sample = length(sample_names)
    labels = Array{Int64,1}(num_sample)
    for i=1:num_sample
        lbs = searchsorted(classes,dt[i,2])
        @assert length(lbs) == 1
        labels[i] = lbs[1]-1
    end
    sample_label_dict = Dict(zip(sample_names,labels))
    JLD.save(sample_label_fl,"sample_label_dict",sample_label_dict)

    label_sample_dict = Dict(map(reverse, collect(sample_label_dict)))
    JLD.save(label_sample_fl, "label_sample_dict", label_sample_dict)

    nothing
end


text:
export vectosymm, symmtovec, dynp_init, dynp_vect2dict, dynp_dict2vect

vectosymm(U::AbstractVector) =
    [U[1] U[2] U[3]
     U[2] U[4] U[5]
     U[3] U[5] U[6]]

symmtovec(S::AbstractMatrix) =
    [S[1,1], S[1,2], S[1,3], S[2,2], S[2,3], S[3,3]]

abstract type DynParmType{T} end

flatter(::Type{DynParmType{:scalar}}) = x->[x]
flatter(::Type{DynParmType{:vec3}}) = identity
flatter(::Type{DynParmType{:symm3}}) = symmtovec

expander(::Type{DynParmType{:scalar}}) = v->v[1]
expander(::Type{DynParmType{:vec3}}) = identity
expander(::Type{DynParmType{:symm3}}) = vectosymm

dynptype(::Type{DynParmType{:scalar}}, T) = T
dynptype(::Type{DynParmType{:vec3}}, T) = Vector{T}
dynptype(::Type{DynParmType{:symm3}}, T) = Matrix{T}

dynplen(::Type{DynParmType{:scalar}}) = 1
dynplen(::Type{DynParmType{:vec3}}) = 3
dynplen(::Type{DynParmType{:symm3}}) = 6

function dynp_init(;args...)
    ndof = length(args[1][1])
    vtype = eltype(eltype(args[1][1]))

    format = Tuple{Symbol, Symbol}[]
    v = Array{vtype, 1}[[] for _ in 1:ndof]

    for (s, (x, typ)) in pairs(args)
        flat = flatter(DynParmType{typ})
        for i in 1:ndof append!(v[i], flat(x[i])) end
        push!(format, (s, typ))
    end

    vcat(v...), format
end

function dynp_vect2dict(v::AbstractVector{T}, ndof::Int, format::AbstractVector{Tuple{Symbol, Symbol}}) where T
    ndynpdof = div(length(v), ndof)
    if !isinteger(ndynpdof)
        throw("Lenght of dynamic parameter vector is not multiple of DOF number")
    end

    d = Dict{Symbol, Any}()

    for (s, typ) in format
        d[s] = dynptype(DynParmType{typ}, T)[]
    end

    offset = 0
    for (s, typ) in format
        len = dynplen(DynParmType{typ})
        expand = expander(DynParmType{typ})
        for i in 1:ndof
            idx = ndynpdof*(i-1)+offset
            range = idx+1:idx+len
            push!(d[s], expand(v[range]))
        end
        offset += len
    end

    d
end

function dynp_dict2vect(d, format::AbstractVector{Tuple{Symbol, Symbol}})
    ndof = length(d[format[1][1]])
    vtype = eltype(eltype(d[format[1][1]]))

    v = Array{vtype, 1}[[] for _ in 1:ndof]

    for (s, typ) in format
        x = d[s]
        flat = flatter(DynParmType{typ})
        for i in 1:ndof append!(v[i], flat(x[i])) end
    end

    vcat(v...)
end


text:
const rev_ops = Dict(
                    :+     => :plus_rev,
                    :*     => :mul_rev,
                    :^     => :power_rev,
                    :-     => :minus_rev,
                    :sqrt  => :sqrt_rev
                    )

function plus_rev(a::Interval, b::Interval, c::Interval)  # a = b + c
    # a = a ∩ (b + c)
    b = b ∩ (a - c)
    c = c ∩ (a - b)

    return a, b, c
end

plus_rev(a,b,c) = plus_rev(promote(a,b,c)...)

function minus_rev(a::Interval, b::Interval, c::Interval)  # a = b - c
    # a = a ∩ (b - c)
    b = b ∩ (a + c)
    c = c ∩ (b - a)

    return a, b, c
end

minus_rev(a,b,c) = minus_rev(promote(a,b,c)...)

function mul_rev(a::Interval, b::Interval, c::Interval)  # a = b * c
    # a = a ∩ (b * c)
    b = b ∩ (a / c)
    c = c ∩ (a / b)

    return a, b, c
end

mul_rev(a,b,c) = mul_rev(promote(a,b,c)...)

Base.iseven(x::Interval) = isinteger(x) && iseven(round(Int, x.lo))

function power_rev(a::Interval, b::Interval, c::Interval)  # a = b^c,  log(a) = c.log(b),  b = a^(1/c)

    # special if c is an even integer: include the possibility of the negative root

    if c == 2  # a = b^2
        b1 = b ∩ √a
        b2 = b ∩ (-√a)

        b = hull(b1, b2)

    elseif iseven(c)
        b1 = b ∩ ( a^(inv(c) ))
        b2 = b ∩ ( -( a^(inv(c)) ) )

        b = hull(b1, b2)

    else

        b = b ∩ ( a^(inv(c) ))
    end

    # a = a ∩ (b ^ c)
    c = c ∩ (log(a) / log(b))

    return a, b, c
end

power_rev(a,b,c) = power_rev(promote(a,b,c)...)

function sqrt_rev(a::Interval, b::Interval)  # a = sqrt(b)
    # a1 = a ∩ √b
    # a2 = a ∩ (-(√b))
    # a = hull(a1, a2)

    b = b ∩ (a^2)

    return a, b
end

sqrt_rev(a,b) = sqrt_rev(promote(a,b)...)

# IEEE-1788 style

function sqr_rev(c, x)   # c = x^2;  refine x
    x1 = sqrt(c) ∩ x
    x2 = -(sqrt(c)) ∩ x

    return hull(x1, x2)
end

sqr_rev(c) = sqr_rev(c, -∞..∞)

"""
∘_rev1(b, c, x) is the subset of x such that x ∘ b is defined and in c
∘_rev2(a, c, x) is the subset of x such that a ∘ x is defined and in c

If these agree (∘ is commutative) then call it ∘_rev(b, c, x)
"""

function mul_rev_new(b, c, x)   # c = b*x
    return x ∩ (c / b)
end

function pow_rev1(b, c, x)   # c = x^b
    return x ∩ c^(1/b)
end

function pow_rev2(a, c, x)   # c = a^x
    return x ∩ (log(c) / lob(a))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Conversion/Promotion

for T in
    (:DayDate, :WeekDate, :MonthDate, :QuarterDate, :SemesterDate, :YearDate, :UndatedDate)
    @eval begin
        $T(dt::Dates.TimeType) = convert($T, dt)
    end
end


text:
isbanded(A::AbstractTriangular) = isbanded(parent(A))
bandwidths(A::Union{UpperTriangular,UnitUpperTriangular}) =
    (min(0,bandwidth(parent(A),1)), bandwidth(parent(A),2))
bandwidths(A::Union{LowerTriangular,UnitLowerTriangular}) =
    (bandwidth(parent(A),1), min(0,bandwidth(parent(A),2)))

triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:BandedColumns} = Tri{ML}()
triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:BandedRows} = Tri{ML}()
triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:ConjLayout{<:BandedRows}} = Tri{ML}()

function tribandeddata(::TriangularLayout{'U'}, A)
    B = triangulardata(A)
    u = bandwidth(B,2)
    D = bandeddata(B)
    view(D, 1:u+1, :)
end

function tribandeddata(::TriangularLayout{'L'}, A)
    B = triangulardata(A)
    l,u = bandwidths(B)
    D = bandeddata(B)
    view(D, u+1:l+u+1, :)
end

tribandeddata(A) = tribandeddata(MemoryLayout(typeof(A)), A)

# Mul

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'U',UNIT,<:BandedColumnMajor},
                                                <:AbstractStridedLayout}) where UNIT
    A,x = M.A,M.B
    tbmv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'L',UNIT,<:BandedColumnMajor},
                                                <:AbstractStridedLayout}) where UNIT
    A,x = M.A,M.B
    tbmv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                                <:AbstractStridedLayout}) where {UPLO,UNIT}
    A,x = M.A,M.B
    tbmv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                                <:AbstractStridedLayout}) where {UPLO,UNIT}
    A,x = M.A,M.B
    tbmv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
end

# Ldiv
for UNIT in ('N', 'U')
    @eval begin
        @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{'U',$UNIT,<:BandedColumnMajor},
                                        <:AbstractStridedLayout})
            A,x = M.A,M.B
            tbsv!('U', 'N', $UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), x)
        end

        @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{'L',$UNIT,<:BandedColumnMajor},
                                                        <:AbstractStridedLayout})
            A,x = M.A,M.B
            tbsv!('L', 'N', $UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), x)
        end
    end
    for UPLO in ('U', 'L')
        @eval begin
            @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{$UPLO,$UNIT,BandedRowMajor},
                                                        <:AbstractStridedLayout})
                A,x = M.A,M.B
                tbsv!($UPLO, 'T', $UNIT, transpose(tribandeddata(A)), x)
            end

            @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{$UPLO,$UNIT,ConjLayout{BandedRowMajor}},
                                                        <:AbstractStridedLayout})
                A,x = M.A,M.B
                tbsv!($UPLO, 'C', $UNIT, tribandeddata(A)', x)
            end
        end
    end
end


text:
import Base: iterate
using Printf
export gmres, gmres!

struct ArnoldiDecomp{T, matT}
    A::matT
    V::Matrix{T} # Orthonormal basis vectors
    H::Matrix{T} # Hessenberg matrix
end

ArnoldiDecomp(A::matT, order::Int, T::Type) where {matT} = ArnoldiDecomp{T, matT}(
    A,
    zeros(T, size(A, 1), order + 1),
    zeros(T, order + 1, order)
)

mutable struct Residual{numT, resT}
    current::resT # Current, absolute, preconditioned residual
    accumulator::resT # Used to compute the residual on the go
    nullvec::Vector{numT} # Vector in the null space of H to compute residuals
    β::resT # the initial residual
end

Residual(order, T::Type) = Residual{T, real(T)}(
    one(real(T)),
    one(real(T)),
    ones(T, order + 1),
    one(real(T))
)

mutable struct GMRESIterable{preclT, precrT, solT, rhsT, vecT, arnoldiT <: ArnoldiDecomp, residualT <: Residual, resT <: Real, orthmethT}
    Pl::preclT
    Pr::precrT
    x::solT
    b::rhsT
    Ax::vecT # Some room to work in.

    arnoldi::arnoldiT
    residual::residualT

    mv_products::Int
    restart::Int
    k::Int
    maxiter::Int
    tol::resT
    β::resT

    orth_meth::orthmethT
end

converged(g::GMRESIterable) = g.residual.current ≤ g.tol

start(::GMRESIterable) = 0

done(g::GMRESIterable, iteration::Int) = iteration ≥ g.maxiter || converged(g)

function iterate(g::GMRESIterable, iteration::Int=start(g))
    # Check for termination first
    if done(g, iteration)
        return nothing
    end

    # Arnoldi step: expand
    expand!(g.arnoldi, g.Pl, g.Pr, g.k, g.Ax)
    g.mv_products += 1

    # Orthogonalize V[:, k + 1] w.r.t. V[:, 1 : k]
    g.arnoldi.H[g.k + 1, g.k] = orthogonalize_and_normalize!(
        view(g.arnoldi.V, :, 1 : g.k),
        view(g.arnoldi.V, :, g.k + 1),
        view(g.arnoldi.H, 1 : g.k, g.k),
        g.orth_meth
    )

    # Implicitly computes the residual
    update_residual!(g.residual, g.arnoldi, g.k)

    g.k += 1

    # Computation of x only at the end of the iterations
    # and at restart.
    if g.k == g.restart + 1 || done(g, iteration + 1)

        # Solve the projected problem Hy = β * e1 in the least-squares sense
        rhs = solve_least_squares!(g.arnoldi, g.β, g.k)

        # And improve the solution x ← x + Pr \ (V * y)
        update_solution!(g.x, view(rhs, 1 : g.k - 1), g.arnoldi, g.Pr, g.k, g.Ax)

        g.k = 1

        # Restart when not done.
        if !done(g, iteration)

            # Set the first basis vector
            g.β = init!(g.arnoldi, g.x, g.b, g.Pl, g.Ax)

            # And initialize the residual
            init_residual!(g.residual, g.β)

            g.mv_products += 1
        end
    end

    g.residual.current, iteration + 1
end

function gmres_iterable!(x, A, b;
                         Pl = Identity(),
                         Pr = Identity(),
                         abstol::Real = zero(real(eltype(b))),
                         reltol::Real = sqrt(eps(real(eltype(b)))),
                         restart::Int = min(20, size(A, 2)),
                         maxiter::Int = size(A, 2),
                         initially_zero::Bool = false,
                         orth_meth::OrthogonalizationMethod = ModifiedGramSchmidt())
    T = eltype(x)

    # Approximate solution
    arnoldi = ArnoldiDecomp(A, restart, T)
    residual = Residual(restart, T)
    mv_products = initially_zero ? 1 : 0

    # Workspace vector to reduce the # allocs.
    Ax = similar(x)
    residual.current = init!(arnoldi, x, b, Pl, Ax, initially_zero = initially_zero)
    init_residual!(residual, residual.current)

    tolerance = max(reltol * residual.current, abstol)

    GMRESIterable(Pl, Pr, x, b, Ax,
        arnoldi, residual,
        mv_products, restart, 1, maxiter, tolerance, residual.current,
        orth_meth
    )
end

"""
    gmres(A, b; kwargs...) -> x, [history]

Same as [`gmres!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
gmres(A, b; kwargs...) = gmres!(zerox(A, b), A, b; initially_zero = true, kwargs...)

"""
    gmres!(x, A, b; kwargs...) -> x, [history]

Solves the problem ``Ax = b`` with restarted GMRES.

# Arguments

- `x`: Initial guess, will be updated in-place;
- `A`: linear operator;
- `b`: right-hand side.

## Keywords

- `initially_zero::Bool`: If `true` assumes that `iszero(x)` so that one
  matrix-vector product can be saved when computing the initial
  residual vector;
- `abstol::Real = zero(real(eltype(b)))`,
  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
  tolerance for the stopping condition
  `|r_k| ≤ max(reltol * |r_0|, abstol)`, where `r_k = A * x_k - b`
- `restart::Int = min(20, size(A, 2))`: restarts GMRES after specified number of iterations;
- `maxiter::Int = size(A, 2)`: maximum number of inner iterations of GMRES;
- `Pl`: left preconditioner;
- `Pr`: right preconditioner;
- `log::Bool`: keep track of the residual norm in each iteration;
- `verbose::Bool`: print convergence information during the iterations.
- `orth_meth::OrthogonalizationMethod = ModifiedGramSchmidt()`: orthogonalization method (ModifiedGramSchmidt(), ClassicalGramSchmidt(), DGKS())

# Return values

**if `log` is `false`**

- `x`: approximate solution.

**if `log` is `true`**

- `x`: approximate solution;
- `history`: convergence history.
"""
function gmres!(x, A, b;
                Pl = Identity(),
                Pr = Identity(),
                abstol::Real = zero(real(eltype(b))),
                reltol::Real = sqrt(eps(real(eltype(b)))),
                restart::Int = min(20, size(A, 2)),
                maxiter::Int = size(A, 2),
                log::Bool = false,
                initially_zero::Bool = false,
                verbose::Bool = false,
                orth_meth::OrthogonalizationMethod = ModifiedGramSchmidt())
    history = ConvergenceHistory(partial = !log, restart = restart)
    history[:abstol] = abstol
    history[:reltol] = reltol
    log && reserve!(history, :resnorm, maxiter)

    iterable = gmres_iterable!(x, A, b; Pl = Pl, Pr = Pr,
                               abstol = abstol, reltol = reltol, maxiter = maxiter,
                               restart = restart, initially_zero = initially_zero,
                               orth_meth = orth_meth)

    verbose && @printf("=== gmres ===\n%4s\t%4s\t%7s\n","rest","iter","resnorm")

    for (iteration, residual) = enumerate(iterable)
        if log
            nextiter!(history)
            history.mvps = iterable.mv_products
            push!(history, :resnorm, residual)
        end

        verbose && @printf("%3d\t%3d\t%1.2e\n", 1 + div(iteration - 1, restart), 1 + mod(iteration - 1, restart), residual)
    end

    verbose && println()
    setconv(history, converged(iterable))
    log && shrink!(history)

    log ? (x, history) : x
end

function update_residual!(r::Residual, arnoldi::ArnoldiDecomp, k::Int)
    if iszero(arnoldi.H[k + 1, k])
        r.current = zero(r.current)
    else
        # Cheaply computes the current residual
        r.nullvec[k + 1] = -conj(dot(view(r.nullvec, 1 : k), view(arnoldi.H, 1 : k, k)) / arnoldi.H[k + 1, k])
        r.accumulator += abs2(r.nullvec[k + 1])
        r.current = r.β / √r.accumulator
    end
end

function init!(arnoldi::ArnoldiDecomp{T}, x, b, Pl, Ax; initially_zero::Bool = false) where {T}
    # Initialize the Krylov subspace with the initial residual vector
    # This basically does V[1] = Pl \ (b - A * x) and then normalize

    first_col = view(arnoldi.V, :, 1)

    copyto!(first_col, b)

    # Potentially save one MV product
    if !initially_zero
        mul!(Ax, arnoldi.A, x)
        first_col .-= Ax
    end

    ldiv!(Pl, first_col)

    # Normalize
    β = norm(first_col)
    first_col .*= inv(β)
    β
end

@inline function init_residual!(r::Residual{numT, resT}, β) where {numT,resT}
    r.accumulator = one(resT)
    r.β = β
end

function solve_least_squares!(arnoldi::ArnoldiDecomp{T}, β, k::Int) where {T}
    # Compute the least-squares solution to Hy = β e1 via Given's rotations
    rhs = zeros(T, k)
    rhs[1] = β

    H = FastHessenberg(view(arnoldi.H, 1 : k, 1 : k - 1))
    ldiv!(H, rhs)

    rhs
end

function update_solution!(x, y, arnoldi::ArnoldiDecomp{T}, Pr::Identity, k::Int, Ax) where {T}
    # Update x ← x + V * y
    mul!(x, view(arnoldi.V, :, 1 : k - 1), y, one(T), one(T))
end

function update_solution!(x, y, arnoldi::ArnoldiDecomp{T}, Pr, k::Int, Ax) where {T}
    # Computing x ← x + Pr \ (V * y) and use Ax as a work space
    mul!(Ax, view(arnoldi.V, :, 1 : k - 1), y)
    ldiv!(Pr, Ax)
    x .+= Ax
end

function expand!(arnoldi::ArnoldiDecomp, Pl::Identity, Pr::Identity, k::Int, Ax)
    # Simply expands by A * v without allocating
    mul!(view(arnoldi.V, :, k + 1), arnoldi.A, view(arnoldi.V, :, k))
end

function expand!(arnoldi::ArnoldiDecomp, Pl, Pr::Identity, k::Int, Ax)
    # Expands by Pl \ (A * v) without allocating
    nextV = view(arnoldi.V, :, k + 1)
    mul!(nextV, arnoldi.A, view(arnoldi.V, :, k))
    ldiv!(Pl, nextV)
end

function expand!(arnoldi::ArnoldiDecomp, Pl, Pr, k::Int, Ax)
    # Expands by Pl \ (A * (Pr \ v)). Avoids allocation by using Ax.
    nextV = view(arnoldi.V, :, k + 1)
    ldiv!(nextV, Pr, view(arnoldi.V, :, k))
    mul!(Ax, arnoldi.A, nextV)
    copyto!(nextV,  Ax)
    ldiv!(Pl, nextV)
end


text:
#!/usr/bin/env julia
"""Lempel-Ziv complexity for a binary sequence, in naive Julia code.

- How to use it? From Julia, it's easy:

julia> using LempelZiv
julia> s = "1001111011000010"
julia> LempelZiv.lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
6

- Run this .jl file with an argument "test" to run a small benchmark (100 strings of size 10000).
- Note: there is also a Python version, if you need.

- MIT Licensed, (C) 2017 Lilian Besson (Naereen)
  https://GitHub.com/Naereen/LempelZiv.jl
"""

__author__ = "Lilian Besson (Naereen)"
__version__ = "0.0.1"

module LempelZiv
export lempel_ziv_complexity

"""
    lempel_ziv_complexity(binary_sequence)

Returns the Lempel-Ziv complexity for a binary sequence, in naive Julia code.

- How to use it? From Julia, it's easy:

```julia
julia> using LempelZiv
julia> s = "1001111011000010"
julia> LempelZiv.lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
6
```

- [MIT Licensed](https://lbesson.mit-license.org)
  (C) 2017 Lilian Besson (Naereen)
  See [Naereen/LempelZiv.jl](https://GitHub.com/Naereen/LempelZiv.jl)
"""
function lempel_ziv_complexity(binary_sequence)
    u, v, w = 0, 1, 1
    v_max = 1
    size = length(binary_sequence)
    complexity = 1
    while true
        if binary_sequence[u + v] == binary_sequence[w + v]
            v += 1
            if w + v >= size
                complexity += 1
                break
            end
        else
            if v > v_max
                v_max = v
            end
            u += 1
            if u == w
                complexity += 1
                w += v_max
                if w > size
                    break
                else
                    u = 0
                    v = 1
                    v_max = 1
                end
            else
                v = 1
            end
        end
    end
    return complexity
end

end

if "test" in ARGS
    # import LempelZiv
    s = "1001111011000010"
    # LempelZiv.lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
    println("For s = ", s, " its Lempel-Ziv Complexity is = ", LempelZiv.lempel_ziv_complexity(s))

    M = 100;
    N = 5000;
    for _ in 1:M
        s = join(rand(0:1, rand(N:10*N)));
        println("For a random string s of size = ", length(s), " its Lempel-Ziv Complexity is = ", LempelZiv.lempel_ziv_complexity(s))
        println(@time LempelZiv.lempel_ziv_complexity(s))
    end
end


text:
"""
    TriComplex{T <: Real} = Complexification{BiComplex{T}}

A tri-complex number is a complexification of bi-complex numbers.
"""
const TriComplex{T <: Real} = Complexification{BiComplex{T}}

function TriComplex(z::BiComplex{T}) where T <: Real
    TriComplex{T}(z, zero(BiComplex{T}))
end

function TriComplex(z::Cmplex{T}) where T <: Real
    TriComplex{T}(
        BiComplex{T}(z, zero(Cmplex{T})),
        zero(BiComplex{T})
    )
end

function TriComplex(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        BiComplex{T}(
            Cmplex{T}(f, g),
            Cmplex{T}(h, j)
        )
    )
end

function TriComplex(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        BiComplex{T}(
            Cmplex{T}(f, g),
            Cmplex{T}(h, zero(T))
        )
    )
end

function TriComplex(a::T, b::T, c::T, d::T, f::T, g::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        BiComplex{T}(
            Cmplex{T}(f, g),
            zero(Cmplex{T})
        )
    )
end

function TriComplex(a::T, b::T, c::T, d::T, f::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        BiComplex{T}(
            Cmplex{T}(f, zero(T)),
            zero(Cmplex{T})
        )
    )
end

function TriComplex(a::T, b::T, c::T, d::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, d)
        ),
        zero(BiComplex{T})
    )
end

function TriComplex(a::T, b::T, c::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            Cmplex{T}(c, zero(T))
        ),
        zero(BiComplex{T})
    )
end

function TriComplex(a::T, b::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, b),
            zero(Cmplex{T})
        ),
        zero(BiComplex{T})
    )
end

function TriComplex(a::T) where T <: Real
    TriComplex{T}(
        BiComplex{T}(
            Cmplex{T}(a, zero(T)),
            zero(Cmplex{T})
        ),
        zero(BiComplex{T})
    )
end

function TriComplex(a::Real, b::Real)
    TriComplex(promote(a, b)...)
end

function TriComplex(a::Real, b::Real, c::Real)
    TriComplex(promote(a, b, c)...)
end

function TriComplex(a::Real, b::Real, c::Real, d::Real)
    TriComplex(promote(a, b, c, d)...)
end

function TriComplex(a::Real, b::Real, c::Real, d::Real, f::Real)
    TriComplex(promote(a, b, c, d, f)...)
end

function TriComplex(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    TriComplex(promote(a, b, c, d, f, g)...)
end

function TriComplex(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    TriComplex(promote(a, b, c, d, f, g, h)...)
end

function TriComplex(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    TriComplex(promote(a, b, c, d, f, g, h, j)...)
end

"""
    unreal(z::TriComplex)

The unreal part of a tri-complex number is a 7-vector.
"""
function unreal(z::TriComplex)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::TriComplex)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", H: ")
    print(io, z.l.l.r)
    print(io, ", G: ")
    print(io, z.l.r.l)
    print(io, ", GH: ")
    print(io, z.l.r.r)
    print(io, ", F: ")
    print(io, z.r.l.l)
    print(io, ", FH: ")
    print(io, z.r.l.r)
    print(io, ", FG: ")
    print(io, z.r.r.l)
    print(io, ", FGH: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function random(::Type{TriComplex{T}}) where T <: Real
    TriComplex{T}(random(BiComplex{T}), random(BiComplex{T}))
end

function barH(z::TriComplex{T}) where T <: Real
    TriComplex{T}(barH(z.l), barH(z.r))
end

function barG(z::TriComplex{T}) where T <: Real
    TriComplex{T}(barG(z.l), barG(z.r))
end

function barF(z::TriComplex{T}) where T <: Real
    TriComplex{T}(z.l, -(z.r))
end

function quadrance(z::TriComplex)
    (z * barF(z)).l
end

function norm(z::TriComplex)
    abs2(quadrance(quadrance(z)))
end

function iszerodivisor(z::TriComplex)
    iszerodivisor(quadrance(z))
end

function inv(z::TriComplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    q = quadrance(z)
    qq = quadrance(q)
    barF(z) * (barG(q) * barH(qq)) / abs2(qq)
end


text:
import Pkg

if dirname(Pkg.project().path) != @__DIR__
    Pkg.activate(@__DIR__)
end

using Test
using LennardJones
using LinearAlgebra

@testset "Test" begin

@test LennardJones.lj_potential_uij(LennardJones.R_MIN) ≈ -1.0

r₁ = [0.916914, 0.111325, 0.594877]
r₂ = [0.454572, 0.979098, 0.988925]
dr = r₂ - r₁
r = norm(dr)
ϵ = 0.8
σ = 0.9
dist = (a, b) -> a - b

u = LennardJones.lj_potential_uij(r, ϵ = ϵ, σ = σ)
@test u ≈ -0.7509440391688087
@test LennardJones.lj_potential_uij(dr, ϵ = ϵ, σ = σ) ≈ u
@test LennardJones.lj_potential_uij(r₁, r₂, ϵ = ϵ, σ = σ, dist = dist) ≈ u

f = LennardJones.lj_potential_fij(dr, ϵ = ϵ, σ = σ)
@test f ≈ [0.7369705345615358, -1.3832252568186918, -0.6281102846008022]
@test LennardJones.lj_potential_fij(r₁, r₂, ϵ = ϵ, σ = σ, dist = dist) ≈ -f

cutoff = 2.5
u₂ = LennardJones.lj_potential_uij_cutoff(dr, ϵ = ϵ, σ = σ, r_cutoff = cutoff)
@test u₂ ≈ -0.7439934985138913
@test LennardJones.lj_potential_uij_cutoff(r₁, r₂, ϵ = ϵ, σ = σ, r_cutoff = cutoff, dist = dist) ≈ u₂

f₂ = LennardJones.lj_potential_fij_cutoff(dr, ϵ = ϵ, σ = σ, r_cutoff = cutoff)
@test f₂ ≈ f
@test LennardJones.lj_potential_fij_cutoff(r₁, r₂, ϵ = ϵ, σ = σ, r_cutoff = cutoff, dist = dist) ≈ -f₂

r₃ = [0.669256, 0.300412, 0.671123]
dr₂ = r₃ - r₁
u₃ = LennardJones.lj_potential_wca_uij(dr₂, ϵ = ϵ, σ = σ)
@test u₃ ≈ 759679.8217248165
@test LennardJones.lj_potential_wca_uij(r₁, r₃, ϵ = ϵ, σ = σ, dist = dist) ≈ u₃

f₃ = LennardJones.lj_potential_wca_fij(dr₂, ϵ = ϵ, σ = σ)
@test f₃ ≈ [-2.1962740885915317e7, 1.6768563042159226e7, 6.761627492701628e6]
@test LennardJones.lj_potential_wca_fij(r₁, r₃, ϵ = ϵ, σ = σ, dist = dist) ≈ -f₃

end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sso_oidc
using AWS.Compat
using AWS.UUIDs

"""
    CreateToken()

Creates and returns an access token for the authorized client. The access token issued will be used to fetch short-term credentials for the assigned roles in the AWS account.

# Required Parameters
- `clientId`: The unique identifier string for each client. This value should come from the persisted result of the RegisterClient API.
- `clientSecret`: A secret string generated for the client. This value should come from the persisted result of the RegisterClient API.
- `deviceCode`: Used only when calling this API for the device code grant type. This short-term code is used to identify this authentication attempt. This should come from an in-memory reference to the result of the StartDeviceAuthorization API.
- `grantType`: Supports grant types for authorization code, refresh token, and device code request.

# Optional Parameters
- `code`: The authorization code received from the authorization service. This parameter is required to perform an authorization grant request to get access to a token.
- `redirectUri`: The location of the application that will receive the authorization code. Users authorize the service to send the request to this location.
- `refreshToken`: The token used to obtain an access token in the event that the access token is invalid or expired. This token is not issued by the service.
- `scope`: The list of scopes that is defined by the client. Upon authorization, this list is used to restrict permissions when granting an access token.
"""
create_token(clientId, clientSecret, deviceCode, grantType; aws_config::AWSConfig=global_aws_config()) = sso_oidc("POST", "/token", Dict{String, Any}("clientId"=>clientId, "clientSecret"=>clientSecret, "deviceCode"=>deviceCode, "grantType"=>grantType); aws_config=aws_config)
create_token(clientId, clientSecret, deviceCode, grantType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sso_oidc("POST", "/token", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientId"=>clientId, "clientSecret"=>clientSecret, "deviceCode"=>deviceCode, "grantType"=>grantType), args)); aws_config=aws_config)

"""
    RegisterClient()

Registers a client with AWS SSO. This allows clients to initiate device authorization. The output should be persisted for reuse through many authentication requests.

# Required Parameters
- `clientName`: The friendly name of the client.
- `clientType`: The type of client. The service supports only public as a client type. Anything other than public will be rejected by the service.

# Optional Parameters
- `scopes`: The list of scopes that are defined by the client. Upon authorization, this list is used to restrict permissions when granting an access token.
"""
register_client(clientName, clientType; aws_config::AWSConfig=global_aws_config()) = sso_oidc("POST", "/client/register", Dict{String, Any}("clientName"=>clientName, "clientType"=>clientType); aws_config=aws_config)
register_client(clientName, clientType, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sso_oidc("POST", "/client/register", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientName"=>clientName, "clientType"=>clientType), args)); aws_config=aws_config)

"""
    StartDeviceAuthorization()

Initiates device authorization by requesting a pair of verification codes from the authorization service.

# Required Parameters
- `clientId`: The unique identifier string for the client that is registered with AWS SSO. This value should come from the persisted result of the RegisterClient API operation.
- `clientSecret`: A secret string that is generated for the client. This value should come from the persisted result of the RegisterClient API operation.
- `startUrl`: The URL for the AWS SSO user portal. For more information, see Using the User Portal in the AWS Single Sign-On User Guide.

"""
start_device_authorization(clientId, clientSecret, startUrl; aws_config::AWSConfig=global_aws_config()) = sso_oidc("POST", "/device_authorization", Dict{String, Any}("clientId"=>clientId, "clientSecret"=>clientSecret, "startUrl"=>startUrl); aws_config=aws_config)
start_device_authorization(clientId, clientSecret, startUrl, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sso_oidc("POST", "/device_authorization", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientId"=>clientId, "clientSecret"=>clientSecret, "startUrl"=>startUrl), args)); aws_config=aws_config)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# Defines functions used to parse the input optimization problem into
# a solvable form including routines used to classify input problems as a
# LP, SOCP, MILP, MISOCP, and convex problem types.
#############################################################################

"""
$(TYPEDSIGNATURES)

Converts `MOI.MAX_SENSE` objective to equivalent `MOI.MIN_SENSE` objective
`max(f) = -min(-f)`.
"""
function convert_to_min!(m::Optimizer)

    m._working_problem._optimization_sense = MOI.MIN_SENSE
    if m._input_problem._optimization_sense === MOI.MAX_SENSE

        obj_type = m._input_problem._objective_type
        if obj_type === SINGLE_VARIABLE
            m._working_problem._objective_type = SCALAR_AFFINE
            m._working_problem._objective_saf = MOIU.operate(-, Float64, m._working_problem._objective_sv)
            m._working_problem._objective_saf_parsed = AffineFunctionIneq(m._working_problem._objective_saf, LT_ZERO)

        elseif obj_type === SCALAR_AFFINE
            m._working_problem._objective_saf = MOIU.operate(-, Float64, m._working_problem._objective_saf)
            m._working_problem._objective_saf_parsed = AffineFunctionIneq(m._working_problem._objective_saf, LT_ZERO)

        elseif obj_type === SCALAR_QUADRATIC
            sqf = m._working_problem._objective_sqf.sqf
            m._working_problem._objective_sqf.sqf = MOIU.operate(-, Float64, sqf)

        elseif obj_type === NONLINEAR

            # updates tape for nlp_data block (used by local optimizer)
            nd = m._working_problem._nlp_data.evaluator.m.nlp_data.nlobj.nd
            pushfirst!(nd, NodeData(JuMP._Derivatives.CALLUNIVAR, 2, -1))
            nd[2] = NodeData(nd[2].nodetype, nd[2].index, 1)
            for i = 3:length(nd)
                @inbounds nd[i] = NodeData(nd[i].nodetype, nd[i].index, nd[i].parent + 1)
            end

            # updates tape used by evaluator for the nonlinear objective (used by the relaxed optimizer)
            nd = m._working_problem._objective_nl.expr.nd
            pushfirst!(nd, NodeData(JuMP._Derivatives.CALLUNIVAR, 2, -1))
            nd[2] = NodeData(nd[2].nodetype, nd[2].index, 1)
            for i = 3:length(nd)
                @inbounds nd[i] = NodeData(nd[i].nodetype, nd[i].index, nd[i].parent + 1)
            end
            I, J, V = findnz(m._working_problem._objective_nl.expr.adj)
            I .+= 1
            J .+= 1
            pushfirst!(I, 2)
            pushfirst!(J, 1)
            pushfirst!(V, true)
            m._working_problem._objective_nl.expr.adj = sparse(I, J, V)

            set_val = copy(m._working_problem._objective_nl.expr.setstorage[1])
            pushfirst!(m._working_problem._objective_nl.expr.setstorage, set_val)
            pushfirst!(m._working_problem._objective_nl.expr.numberstorage, 0.0)
            pushfirst!(m._working_problem._objective_nl.expr.isnumber, false)
        end
    end

    return nothing
end

"""

Performs an epigraph reformulation assuming the working_problem is a minimization problem.
"""
function reform_epigraph!(m::Optimizer)

    if m._parameters.presolve_epigraph_flag
        #=
        # add epigraph variable
        obj_variable_index = MOI.add_variable(m)

        # converts ax + b objective to ax - y <= -b constraint with y objective
        obj_type = m._working_problem._objective_type
        if obj_type === SCALAR_AFFINE

            # update unparsed expression
            objective_saf = m._working_problem._objective_saf
            push!(objective_saf, SAT(-1.0, obj_variable_index))
            obj_ci = MOI.add_constraint(m, saf, LT(-objective_saf.constant))

            # update parsed expression (needed for interval bounds)

        # converts ax + b objective to ax - y <= -b constraint with y objective
        elseif obj_type === SCALAR_QUADRATIC

            # update parsed expression
            objective_sqf = m._working_problem._objective_sqf
            obj_ci = MOI.add_constraint(m, saf, LT())

        elseif obj_type === NONLINEAR

            # updated parsed expressions
            objective_nl = m._working_problem._objective_nl

        end

        MOI.set(m, MOI.ObjectiveFunction{SV}(), SV(obj_variable_index))
        =#
    end

    return nothing
end

function check_set_is_fixed(v::VariableInfo)
    v.is_fixed && return true
    v.is_fixed = x.lower_bound === x.upper_bound
    return v.is_fixed
end

"""
$(TYPEDSIGNATURES)

Detects any variables set to a fixed value by equality or inequality constraints
and populates the `_fixed_variable` storage array.
"""
function label_fixed_variables!(m::Optimizer)
    map!(x -> check_set_is_fixed(x), m._fixed_variable, m._working_problem._variable_info)
end

"""
$(TYPEDSIGNATURES)

Detects any variables participating in nonconvex terms and populates the
`_branch_variables` storage array.
"""
function label_branch_variables!(m::Optimizer)

    m._user_branch_variables = !isempty(m._parameters.branch_variable)
    if m._user_branch_variables
        append!(m._branch_variables, m._parameters.branch_variable)
    else

        append!(m._branch_variables, fill(false, m._working_problem._variable_count))

        # adds nonlinear terms in quadratic constraints
        sqf_leq = m._working_problem._sqf_leq
        for i = 1:m._working_problem._sqf_leq_count
            quad_ineq = @inbounds sqf_leq[i]
            for term in quad_ineq.func.quadratic_terms
                variable_index_1 = term.variable_index_1.value
                variable_index_2 = term.variable_index_2.value
                @inbounds m._branch_variables[variable_index_1] = true
                @inbounds m._branch_variables[variable_index_2] = true
            end
        end

        sqf_eq = m._working_problem._sqf_eq
        for i = 1:m._working_problem._sqf_eq_count
            quad_eq = @inbounds sqf_eq[i]
            for term in quad_eq.func.quadratic_terms
                variable_index_1 = term.variable_index_1.value
                variable_index_2 = term.variable_index_2.value
                @inbounds m._branch_variables[variable_index_1] = true
                @inbounds m._branch_variables[variable_index_2] = true
            end
        end

        obj_type = m._working_problem._objective_type
        if obj_type === SCALAR_QUADRATIC
            for term in m._working_problem._objective_sqf.func.quadratic_terms
                variable_index_1 = term.variable_index_1.value
                variable_index_2 = term.variable_index_2.value
                @inbounds m._branch_variables[variable_index_1] = true
                @inbounds m._branch_variables[variable_index_2] = true
            end
        end

        # label nonlinear branch variables (assumes affine terms have been extracted)
        nl_constr = m._working_problem._nonlinear_constr
        for i = 1:m._working_problem._nonlinear_count
            nl_constr_eq = @inbounds nl_constr[i]
            grad_sparsity = nl_constr_eq.expr.grad_sparsity
            for indx in grad_sparsity
                @inbounds m._branch_variables[indx] = true
            end
        end

        if obj_type === NONLINEAR
            grad_sparsity = m._working_problem._objective_nl.expr.grad_sparsity
            for indx in grad_sparsity
                @inbounds m._branch_variables[indx] = true
            end
        end
    end

    # add a map of branch/node index to variables in the continuous solution
    for i = 1:m._working_problem._variable_count
        if m._working_problem._variable_info[i].is_fixed
            m._branch_variables[i] = false
            continue
        end
        if m._branch_variables[i]
            push!(m._branch_to_sol_map, i)
        end
    end

    # creates reverse map
    m._sol_to_branch_map = zeros(m._working_problem._variable_count)
    for i = 1:length(m._branch_to_sol_map)
        j = m._branch_to_sol_map[i]
        m._sol_to_branch_map[j] = i
    end

    # adds branch solution to branch map to evaluator
    m._working_problem._relaxed_evaluator.node_to_variable_map = m._branch_to_sol_map
    m._working_problem._relaxed_evaluator.variable_to_node_map = m._sol_to_branch_map
    m._working_problem._relaxed_evaluator.node_count = length(m._branch_to_sol_map)

    return nothing
end

add_nonlinear_functions!(m::Optimizer) = add_nonlinear_functions!(m, m._input_problem._nlp_data.evaluator)

add_nonlinear_functions!(m::Optimizer, evaluator::Nothing) = nothing
add_nonlinear_functions!(m::Optimizer, evaluator::EmptyNLPEvaluator) = nothing
function add_nonlinear_functions!(m::Optimizer, evaluator::JuMP.NLPEvaluator)

    nlp_data = m._input_problem._nlp_data
    MOI.initialize(evaluator, Symbol[:Grad, :ExprGraph])

    # set nlp data structure
    m._working_problem._nlp_data = nlp_data

    # add subexpressions (assumes they are already ordered by JuMP)
    # creates a dictionary that lists the subexpression sparsity
    # by search each node for variables dict[2] = [2,3] indicates
    # that subexpression 2 depends on variables 2 and 3
    # this is referenced when subexpressions are called by other
    # subexpressions or functions to determine overall sparsity
    # the sparsity of a function is the collection of indices
    # in all participating subexpressions and the function itself
    # it is necessary to define this as such to enable reverse
    # McCormick constraint propagation
    relax_evaluator = m._working_problem._relaxed_evaluator
    has_subexpressions = length(evaluator.m.nlp_data.nlexpr) > 0
    dict_sparsity = Dict{Int64,Vector{Int64}}()
    if has_subexpressions

        println("########################")
        println(length(evaluator.subexpressions))
        println("########################")

        for i = 1:length(evaluator.subexpressions)
            try
                println("########################")
                println(i)
                println(evaluator.subexpressions[i])
                println("########################")

                subexpr = evaluator.subexpressions[i]
                push!(relax_evaluator.subexpressions, NonlinearExpression!(subexpr, dict_sparsity, i,
                                                                          evaluator.subexpression_linearity,
                                                                          m._parameters.relax_tag))
            catch
                @warn "fail to get evaluator.subexpressions[i]"
            end

        end
    end

    # scrubs udf functions using Cassette to remove odd data structures...
    # alternatively convert udfs to JuMP scripts...
    m._parameters.presolve_scrubber_flag && Script.scrub!(m._working_problem._nlp_data)
    if m._parameters.presolve_to_JuMP_flag
        Script.udf_loader!(m)
    end

    parameter_values = copy(evaluator.parameter_values)

    # add nonlinear objective
    if evaluator.has_nlobj
        m._working_problem._objective_nl = BufferedNonlinearFunction(evaluator.objective, MOI.NLPBoundsPair(-Inf, Inf),
                                                                     dict_sparsity, evaluator.subexpression_linearity,
                                                                     m._parameters.relax_tag)
    end

    # add nonlinear constraints
    constraint_bounds = m._working_problem._nlp_data.constraint_bounds
    for i = 1:length(evaluator.constraints)
        constraint = evaluator.constraints[i]
        bnds = constraint_bounds[i]
        push!(m._working_problem._nonlinear_constr, BufferedNonlinearFunction(constraint, bnds, dict_sparsity,
                                                                              evaluator.subexpression_linearity,
                                                                              m._parameters.relax_tag))
    end

    m._input_problem._nonlinear_count = length(m._working_problem._nonlinear_constr)
    m._working_problem._nonlinear_count = length(m._working_problem._nonlinear_constr)

    return nothing
end

function add_nonlinear_evaluator!(m::Optimizer)
    evaluator = m._input_problem._nlp_data.evaluator
    add_nonlinear_evaluator!(m, evaluator)
    return nothing
end

add_nonlinear_evaluator!(m::Optimizer, evaluator::Nothing) = nothing
add_nonlinear_evaluator!(m::Optimizer, evaluator::EmptyNLPEvaluator) = nothing
function add_nonlinear_evaluator!(m::Optimizer, evaluator::JuMP.NLPEvaluator)
    m._working_problem._relaxed_evaluator = Evaluator()

    relax_evaluator = m._working_problem._relaxed_evaluator
    relax_evaluator.variable_count = length(m._working_problem._variable_info)
    relax_evaluator.user_operators = evaluator.m.nlp_data.user_operators

    relax_evaluator.lower_variable_bounds = zeros(relax_evaluator.variable_count)
    relax_evaluator.upper_variable_bounds = zeros(relax_evaluator.variable_count)
    relax_evaluator.x                     = zeros(relax_evaluator.variable_count)
    relax_evaluator.num_mv_buffer         = zeros(relax_evaluator.variable_count)
    relax_evaluator.treat_x_as_number     = fill(false, relax_evaluator.variable_count)
    relax_evaluator.ctx                   = GuardCtx(metadata = GuardTracker(m._parameters.domain_violation_ϵ,
                                                                             m._parameters.domain_violation_guard_on))
    relax_evaluator.subgrad_tol           = m._parameters.subgrad_tol

    m._nonlinear_evaluator_created = true

    return nothing
end

function add_subexpression_buffers!(m::Optimizer)
    relax_evaluator = m._working_problem._relaxed_evaluator
    relax_evaluator.subexpressions_eval = fill(false, length(relax_evaluator.subexpressions))

    return nothing
end

"""
Translates input problem to working problem. Routines and checks and optional manipulation is left to the presolve stage.
"""
function initial_parse!(m::Optimizer)

    # reset initial time and solution statistics
    m._time_left = m._parameters.time_limit

    # add variables to working model
    ip = m._input_problem
    append!(m._working_problem._variable_info, ip._variable_info)
    m._working_problem._variable_count = ip._variable_count

    # add linear constraints to the working problem
    linear_leq = ip._linear_leq_constraints
    for i = 1:ip._linear_leq_count
        linear_func, leq_set = @inbounds linear_leq[i]
        push!(m._working_problem._saf_leq, AffineFunctionIneq(linear_func, leq_set))
        m._working_problem._saf_leq_count += 1
    end

    linear_geq = ip._linear_geq_constraints
    for i = 1:ip._linear_geq_count
        linear_func, geq_set = @inbounds linear_geq[i]
        push!(m._working_problem._saf_leq, AffineFunctionIneq(linear_func, geq_set))
        m._working_problem._saf_leq_count += 1
    end

    linear_eq = ip._linear_eq_constraints
    for i = 1:ip._linear_eq_count
        linear_func, eq_set = @inbounds linear_eq[i]
        push!(m._working_problem._saf_eq, AffineFunctionEq(linear_func, eq_set))
        m._working_problem._saf_eq_count += 1
    end

    # add quadratic constraints to the working problem
    quad_leq = ip._quadratic_leq_constraints
    for i = 1:ip._quadratic_leq_count
        quad_func, leq_set = @inbounds quad_leq[i]
        push!(m._working_problem._sqf_leq, BufferedQuadraticIneq(quad_func, leq_set))
        m._working_problem._sqf_leq_count += 1
    end

    quad_geq = ip._quadratic_geq_constraints
    for i = 1:ip._quadratic_geq_count
        quad_func, geq_set = @inbounds quad_geq[i]
        push!(m._working_problem._sqf_leq, BufferedQuadraticIneq(quad_func, geq_set))
        m._working_problem._sqf_leq_count += 1
    end

    quad_eq = ip._quadratic_eq_constraints
    for i = 1:ip._quadratic_eq_count
        quad_func, eq_set = @inbounds quad_eq[i]
        push!(m._working_problem._sqf_eq, BufferedQuadraticEq(quad_func, eq_set))
        m._working_problem._sqf_eq_count += 1
    end

    # add conic constraints to the working problem
    soc_vec = m._input_problem._conic_second_order
    for i = 1:ip._conic_second_order_count
        soc_func, soc_set = @inbounds soc_vec[i]
        first_variable_loc = soc_func.variables[1].value
        prior_lbnd = m._working_problem._variable_info[first_variable_loc].lower_bound
        m._working_problem._variable_info[first_variable_loc].lower_bound = max(prior_lbnd, 0.0)
        push!(m._working_problem._conic_second_order, BufferedSOC(soc_func, soc_set))
        m._working_problem._conic_second_order_count += 1
    end

    # set objective function
    m._working_problem._objective_type = ip._objective_type
    m._working_problem._objective_sv = ip._objective_sv
    m._working_problem._objective_saf = ip._objective_saf
    m._working_problem._objective_saf_parsed = AffineFunctionIneq(ip._objective_saf, LT_ZERO)
    m._working_problem._objective_sqf = BufferedQuadraticIneq(ip._objective_sqf, LT_ZERO)

    # add nonlinear constraints
    # the nonlinear evaluator loads with populated subexpressions which are then used
    # to asssess the linearity of subexpressions
    add_nonlinear_evaluator!(m)
    add_nonlinear_functions!(m)
    add_subexpression_buffers!(m)

    # converts a maximum problem to a minimum problem (internally) if necessary
    # this is placed after adding nonlinear functions as this prior routine
    # copies the nlp_block from the input_problem to the working problem
    convert_to_min!(m)
    reform_epigraph!(m)

    # labels the variable info and the _fixed_variable vector for each fixed variable
    label_fixed_variables!(m)

    # labels variables to branch on
    label_branch_variables!(m)

    # updates run and parse times
    new_time = time() - m._start_time
    m._parse_time = new_time
    m._run_time = new_time

    return nothing
end

### Routines for parsing the full nonconvex problem
"""
[FUTURE FEATURE] Reformulates quadratic terms in SOC constraints if possible.
For <= or >=, the quadratic term is deleted if an SOCP is detected. For ==,
the SOC check is done for each >= and <=, the convex constraint is reformulated
to a SOC, the concave constraint is keep as a quadratic.
"""
function parse_classify_quadratic!(m::Optimizer)
    #=
    for (id, cinfo) in m._quadratic_constraint
        is_soc, add_concave, cfunc, cset, qfunc, qset = check_convexity(cinfo.func, cinfo.set)
        if is_soc
            MOI.add_constraint(m, cfunc, cset)
            deleteat!(m._quadratic_constraint, id)
            if add_concave
                MOI.add_constraint(m, qfunc, qset)
            end
        end
    end
    =#
    return nothing
end

"""
[FUTURE FEATURE] Parses provably convex nonlinear functions into a convex
function buffer
"""
function parse_classify_nlp(m)
    return nothing
end

"""
Classifies the problem type
"""
function parse_classify_problem!(m::Optimizer)

    ip = m._input_problem
    integer_variable_number = count(is_integer.(ip._variable_info))

    nl_expr_number = ip._objective_type === NONLINEAR ? 1 : 0
    nl_expr_number += ip._nonlinear_count
    cone_constraint_number = ip._conic_second_order_count
    quad_constraint_number = ip._quadratic_leq_count + ip._quadratic_geq_count + ip._quadratic_eq_count

    linear_or_sv_objective = (ip._objective_type === SINGLE_VARIABLE || ip._objective_type === SCALAR_AFFINE)
    relaxed_supports_soc = false
    #TODO: relaxed_supports_soc = MOI.supports_constraint(m.relaxed_optimizer, VECOFVAR, SOC)

    if integer_variable_number === 0

        if cone_constraint_number === 0 && quad_constraint_number === 0 &&
            nl_expr_number === 0 && linear_or_sv_objective
            m._working_problem._problem_type = LP
            #println("LP")
        elseif quad_constraint_number === 0 && relaxed_supports_soc &&
               nl_expr_number === 0 && linear_or_sv_objective
            m._working_problem._problem_type = SOCP
            #println("SOCP")
        else
            #parse_classify_quadratic!(m)
            #if iszero(m._input_nonlinear_constraint_number)
            #    if isempty(m._quadratic_constraint)
            #        m._problem_type = SOCP
            #    end
            #else
            #    # Check if DIFF_CVX, NS_CVX, DIFF_NCVX, OR NS_NCVX
            #    m._problem_type = parse_classify_nlp(m)
            #end
            m._working_problem._problem_type = MINCVX
            #println("MINCVX")
        end
    else
        #=
        if cone_constraint_number === 0 && quad_constraint_number === 0 && linear_or_sv_objective
        elseif quad_constraint_number === 0 && relaxed_supports_soc && linear_or_sv_objective
            m._working_problem._problem_type = MISOCP
        else
            #parse_classify_quadratic!(m)
            #=
            if iszero(m._nonlinear_constraint_number)
                if iszero(m._quadratic_constraint_number)
                    m._problem_type = MISOCP
                end
            else
                # Performs parsing
                _ = parse_classify_nlp(m)
            end
            =#
            m._problem_type = MINCVX
        end
        =#
    end

    return nothing
end

"""

Basic parsing for global solutions (no extensive manipulation)
"""
function parse_global!(t::ExtensionType, m::Optimizer)
    return nothing
end


text:
"""
The goal of the module is to attach physical coordinaes to array elements, so `a[i,j]` could be interpreted as ``f(x_i, y_j)`` for some function ``f``.
"""
module SampledDomains

abstract type AbstractDomain end

import Base.:*
import Base.size

export CartesianDomain2D, make_centered_domain2D, dualRange, dualDomain

# TODO rewrite as parametric type
struct CartesianDomain2D<:AbstractDomain
    xrange::AbstractRange
    yrange::AbstractRange
end

getranges(dom::AbstractDomain) = [getfield(dom,n) for n in fieldnames(typeof(dom))]

size(dom::AbstractDomain) = tuple(length.(reverse(getranges(dom)))...) # we need reverse to follow column major rule
# size(dom::AbstractDomain, d) = length(getfield(dom,d))
size(dom::AbstractDomain, d) = size(dom)[d]

*(dom::CartesianDomain2D, s::Real) = CartesianDomain2D(dom.xrange * s, dom.yrange*s)
*( s::Real, dom::CartesianDomain2D) = CartesianDomain2D(dom.xrange * s, dom.yrange*s)

function make_centered_domain2D(xlength, ylength, pixelsize)
    xrange = ((1:xlength) .- (1 + xlength) / 2) .* pixelsize
    yrange = ((1:ylength) .- (1 + ylength) / 2) .* pixelsize
    CartesianDomain2D(xrange, yrange)
end

"""
    dualRange(xrange::AbstractRange, q::Int =1)

Construct range in Fourier-transform-dual domain with upsampling factor `q`. If sampling in `x` has step size Δx, the dual
domain is sampled in interval (-q/2Δx, q/2Δx).
"""
function dualRange(xrange::AbstractRange, q::Int =1)
    len = length(xrange) *q
    st = step(xrange) /q
    return UnitRange(-floor(Int, len / 2), ceil(Int, len / 2) - 1) / (st * len)

end

function dualDomain(dom::CartesianDomain2D, q::Int = 1)
    kxrange = dualRange(dom.xrange, q)
    kyrange = dualRange(dom.yrange, q)
    return CartesianDomain2D(kxrange,kyrange)
end

"""
    SampledDomain(vals::Array, dom::AbstractDomain)
    SampledDomain(f::Function, dom::AbstractDomain)

Contains two fileds: `vals` and `dom` representing sampled values of function `f` on `dom`ain.
"""
struct SampledDomain
    vals::Array
    dom::AbstractDomain
    SampledDomain(vals::Array, dom::AbstractDomain) = size(vals) != size(dom) ? ErrorException("Different size of the domain and values array") : new(vals, dom)
end

SampledDomain(f::Function, dom::AbstractDomain) = SampledDomain(map(x->f(x...), collect(Iterators.product(reverse(getranges(dom))...))), dom)

end


text:
using Catlab.WiringDiagrams, Catlab.Graphs

function to_parts(xs::Vector{Int}; distinguish::Int=0)::String
  join([let xi=(findall(==(i), xs)), x=join(xi, " "), b=i==distinguish;
          b ? ""*x*"" : x end
        for i in 1:max0(xs)]," | ")
end

function view_cdict(cd::CDict; distinguish::Pair{Symbol, Int}=(:X=>0))::String
  ds, di = distinguish
  join(["$k:$(to_parts(v; distinguish=k==ds ? di : 0))" for (k,v) in pairs(cd)],"
 ") * " "
end

"""Adds tree to wiring diagram, returns box id of root"""
function add_tree!(wd::WiringDiagram, t::Tree)::Int
  n = length(t.children)
  bx= Box(view_cdict(t.coloring), [Int],repeat([Int], n))
  root = add_box!(wd, bx)
  for (i,(k,v)) in enumerate(t.children)
    child = add_tree!(wd, v)
    add_wire!(wd, Wire(k, (root, i), (child, 1)))
  end
  return root
end

function view_tree(t::Tree)::WiringDiagram
  wd = WiringDiagram([],[])
  add_tree!(wd, t)
  return wd
end

struct VizState
  wd::WiringDiagram
  tree::Tree
  msg::String
  boxids::Dict{VPSI,Int}
  colors::Dict{Int,String}
end

VizState(wd::WiringDiagram,tree::Tree,msg::String,boxids::Dict{VPSI,Int}) = VizState(wd,tree,msg,boxids,Dict{Int,String}())

function to_vizstate(hs::Vector{History}; sat::Bool=false)::Vector{VizState}
  vs = VizState[]
  add_history!(vs, hs; sat)
  return vs
end

function add_history!(vs::Vector{VizState}, hs::Vector{History}; sat::Bool=false)
  if isempty(hs)
    return nothing
  end
  h = hs[1]
  v = isempty(vs) ? nothing : deepcopy(vs[end])

  if (h.action == "start")
    isnothing(v) || error()
  elseif h.action == "start_iter"
    # unpack history
    split_seq, color_seq, new_ind, splt_data = h.val
    split_tab, split_color, split_inds = splt_data
    # Load data from previous vizstate
    if isnothing(v)
      wd, vt, bi, c = WiringDiagram([],[]), Tree(), Dict{VPSI,Int}(), Dict{Int,String}()
    else
      wd, vt, bi, c = v.wd, v.tree, v.boxids, Dict([
        b=>(clr == "lightblue" ? "bisque" : clr) for (b,clr) in v.colors])
    end

    # Add new box
    bi[split_seq] = boxid = add_box!(wd, Box(view_cdict(color_seq[1]),[Int],[Int]))
    c[boxid] = "lightblue"

    if !isempty(split_seq)
      add_wire!(wd, Wire(split_seq[end][2], (bi[split_seq[1:end-1]], 1), (boxid, 1)))
    end

    # Update tree
    vt[split_seq].coloring = color_seq[1]
    vt[split_seq].saturated = color_seq[end]
    vt[split_seq].indicator = new_ind
    for i in split_inds
      vt[split_seq].children[split_tab=>i] = Tree()
    end

    push!(vs, VizState(deepcopy(wd), vt, "Add node", bi, c))
      # Add states for each color saturation step
    if sat
      for (i, color) in enumerate(color_seq[2:end])
        set_subpart!(wd.diagram, boxid, :value,  view_cdict(color))
        msg = "Color saturate $i/$(length(color_seq)-1)"
        push!(vs, VizState(deepcopy(wd), vt, msg, bi, deepcopy(c)))
      end
    end
    # Add state for splitting
    vt.saturated, vt.indicator = color_seq[end], new_ind
    set_subpart!(wd.diagram, boxid, :value,  view_cdict(
      color_seq[end]; distinguish=(split_tab => split_color)))
    col = deepcopy(c)
    if split_tab == :_nothing # we have a leaf
      msg = "Found leaf"
      col[boxid] = "green"
    else
      msg = "Split on $splt_data"
    end
    push!(vs, VizState(deepcopy(wd), vt, msg, bi, col))
  elseif h.action == "auto_prune"
    push!(vs, VizState(v.wd,v.tree,"AUTO-PRUNE $(h.val)", v.boxids,v.colors))
  elseif h.action == "order_prune"
    push!(vs, VizState(v.wd,v.tree,"ORDER-PRUNE $(h.val)", v.boxids,v.colors))
  elseif h.action == "orbit_prune"
    push!(vs, VizState(v.wd,v.tree,"ORBIT-PRUNE $(h.val)", v.boxids,v.colors))
  elseif h.action == "flag_skip"
    push!(vs, VizState(v.wd,v.tree,"FLAGGED FOR SKIPPING", v.boxids,v.colors))
  elseif h.action == "return"
    # for everything below h.val (a split_seq) remove coloring
    n = length(h.val)
    newcolor = Dict{Int,String}()
    for (b,clr) in v.colors
      pth=only([seq for (seq,bxid) in v.boxids if bxid==b])
      if pth == h.val
        newcolor[b]="lightblue"
      elseif common(pth, h.val) < n
        newcolor[b]=clr
      end
    end
    v = VizState(v.wd,v.tree, "Return", v.boxids, newcolor)
    push!(vs, v)
  else
    error((h.action, v))
  end
  add_history!(vs, hs[2:end])
end

function process_history(hs::Vector{History})::Vector{Pair{WiringDiagram, String}}
  @assert hs[1].action == "start"
  res, wd = Pair{WiringDiagram, String}[], WiringDiagram([],[])
  add_box!(wd, Box(view_cdict(t.coloring),[],[]))
  for h in hs[2:end]
    push!(res, add_history(wd, h))
  end
  return res
end