text:
# Functional implementation of https://egraphs-good.github.io/
# https://dl.acm.org/doi/10.1145/3434304

using DataStructures
using TermInterface

"""
Abstract type representing an [`EGraph`](@ref) analysis,
attaching values from a join semi-lattice domain to
an EGraph
"""
const ClassMem = Dict{EClassId,EClass}
const HashCons = Dict{AbstractENode,EClassId}
const Analyses = Set{Type{<:AbstractAnalysis}}
const SymbolCache = Dict{Any, Set{EClassId}}
const TermTypes = Dict{Tuple{Any, EClassId}, Type}

"""
A concrete type representing an [`EGraph`].
See the [egg paper](https://dl.acm.org/doi/pdf/10.1145/3434304)
for implementation details
"""
mutable struct EGraph
    """stores the equality relations over e-class ids"""
    # uf::IntDisjointSets{EClassId}
    uf::IntDisjointSet
    """map from eclass id to eclasses"""
    classes::ClassMem
    memo::HashCons             # memo
    """worklist for ammortized upwards merging"""
    dirty::Vector{EClassId}
    root::EClassId
    """A vector of analyses associated to the EGraph"""
    analyses::Analyses
    # """
    # a cache mapping function symbols to e-classes that
    # contain e-nodes with that function symbol.
    # """
    # symcache::SymbolCache
    default_termtype::Type
    termtypes::TermTypes
    numclasses::Int
    numnodes::Int
    # number of rules that have been applied
    # age::Int
end

function EGraph()
    EGraph(
        IntDisjointSet{EClassId}(),
        # IntDisjointSets{EClassId}(0),
        ClassMem(),
        HashCons(),
        # ParentMem(),
        EClassId[],
        -1,
        Analyses(),
        # SymbolCache(),
        Expr,
        TermTypes(),
        0,
        0,
        # 0
    )
end

function EGraph(e; keepmeta=false)
    g = EGraph()
    if keepmeta
        push!(g.analyses, MetadataAnalysis)
    end

    rootclass, rootnode = addexpr!(g, e; keepmeta=keepmeta)
    g.root = rootclass.id
    g
end

function settermtype!(g::EGraph, f, ar, T)
    g.termtypes[(f,ar)] = T
end

function settermtype!(g::EGraph, T)
    g.default_termtype = T
end

function gettermtype(g::EGraph, f, ar)
    if haskey(g.termtypes, (f,ar))
        g.termtypes[(f,ar)]
    else
        g.default_termtype
    end
end

"""
Returns the canonical e-class id for a given e-class.
"""
# function find(g::EGraph, a::EClassId)::EClassId
#     find_root_if_normal(g.uf, a)
# end
function find(g::EGraph, a::EClassId)::EClassId
    find_root(g.uf, a)
end
find(g::EGraph, a::EClass)::EClassId = find(g, a.id)

function geteclass(g::EGraph, a::EClassId)::EClass
    id = find(g, a)
    ec = g.classes[id]
    # @show ec.id id a
    # @assert ec.id == id
    # ec.id = id
    ec
end
# geteclass(g::EGraph, a::EClass)::EClassId = geteclass()
Base.getindex(g::EGraph, i::EClassId) = geteclass(g, i)

### Definition 2.3: canonicalization
iscanonical(g::EGraph, n::ENodeTerm) = n == canonicalize(g, n)
iscanonical(g::EGraph, n::ENodeLiteral) = true
iscanonical(g::EGraph, e::EClass) = find(g, e.id) == e.id

canonicalize(g::EGraph, n::ENodeLiteral) = n

function canonicalize(g::EGraph, n::ENodeTerm{T}) where {T}
    if arity(n) > 0
        new_args = map(x -> find(g, x), arguments(n))
        return ENodeTerm{T}(exprhead(n), operation(n), new_args)
    end
    return n
end

function canonicalize!(g::EGraph, n::ENodeTerm)
    args = arguments(n)
    for i ∈ 1:arity(n)
        args[i] = find(g, args[i])
    end
    n.hash[] = UInt(0)
    return n
end

canonicalize!(g::EGraph, n::ENodeLiteral) = n

function canonicalize!(g::EGraph, e::EClass)
    e.id = find(g, e.id)
end

function lookup(g::EGraph, n::AbstractENode)
    cc = canonicalize(g, n)
    if !haskey(g.memo, cc)
        return nothing
    end
    return find(g, g.memo[cc])
end

"""
Inserts an e-node in an [`EGraph`](@ref)
"""
function add!(g::EGraph, n::AbstractENode)::EClass
    @debug("adding ", n)

    n = canonicalize(g, n)
    if haskey(g.memo, n)
        eclass = geteclass(g, g.memo[n])
        return eclass
    end
    @debug(n, " not found in memo")

    id = push!(g.uf) # create new singleton eclass

    if n isa ENodeTerm
        for c_id ∈ arguments(n)
            addparent!(g.classes[c_id], n, id)
        end
    end

    g.memo[n] = id

    classdata = EClass(id, AbstractENode[n], Pair{AbstractENode, EClassId}[])
    g.classes[id] = classdata
    g.numclasses += 1

    for an ∈ g.analyses
        if !islazy(an) && an !== MetadataAnalysis
            setdata!(classdata, an, make(an, g, n))
            modify!(an, g, id)
        end
    end
    return classdata
end

"""
Extend this function on your types to do preliminary
preprocessing of a symbolic term before adding it to
an EGraph. Most common preprocessing techniques are binarization
of n-ary terms and metadata stripping.
"""
function preprocess(e::Expr)
    cleanast(e)
end
preprocess(x) = x

"""
Recursively traverse an type satisfying the `TermInterface` and insert terms into an
[`EGraph`](@ref). If `e` has no children (has an arity of 0) then directly
insert the literal into the [`EGraph`](@ref).
"""
addexpr!(g::EGraph, se::EClass; keepmeta=false) = (se, se[1])

function addexpr!(g::EGraph, se; keepmeta=false)::Tuple{EClass, AbstractENode}
    # println("========== $e ===========")
    e = preprocess(se)
    T = typeof(e)
    node = nothing

    if istree(T)
        exhead = exprhead(e)
        op = operation(e)
        args = arguments(e)

        n = length(args)

        class_ids = EClassId[
            first(addexpr!(g, child; keepmeta=keepmeta)).id
        for child in args]

        node = ENodeTerm{typeof(e)}(exhead, op, class_ids)
    else
        # constant enode
        node = ENodeLiteral(e)
    end

    ec = add!(g, node)
    if keepmeta
        # TODO check if eclass already has metadata?
        meta = TermInterface.metadata(e)
        setdata!(ec, MetadataAnalysis, meta)
    end
    return (ec, node)
end

"""
Given an [`EGraph`](@ref) and two e-class ids, set
the two e-classes as equal.
"""
function Base.merge!(g::EGraph, a::EClassId, b::EClassId)::EClassId
    id_a = find(g, a)
    id_b = find(g, b)

    id_a == id_b && return id_a
    to = union!(g.uf, id_a, id_b)

    @debug "merging" id_a id_b

    from = (to == id_a) ? id_b : id_a

    push!(g.dirty, to)

    from_class = g.classes[from]
    to_class = g.classes[to]
    to_class.id = to

    # I (was) the troublesome line!
    g.classes[to] = union!(to_class, from_class)
    delete!(g.classes, from)
    g.numclasses -= 1

    return to
end

function in_same_class(g::EGraph, a, b)
    find(g, a) == find(g, b)
end

# TODO new rebuilding from egg
"""
This function restores invariants and executes
upwards merging in an [`EGraph`](@ref). See
the [egg paper](https://dl.acm.org/doi/pdf/10.1145/3434304)
for more details.
"""
function rebuild!(g::EGraph)
    # normalize!(g.uf)

    while !isempty(g.dirty)
        # todo = unique([find(egraph, id) for id ∈ egraph.dirty])
        todo = unique(g.dirty)
        empty!(g.dirty)
        for x ∈ todo
            repair!(g, x)
        end
    end

    if g.root != -1
        g.root = find(g, g.root)
    end

    normalize!(g.uf)

    # for i ∈ 1:length(egraph.uf)
    #     find_root!(egraph.uf, i)
    # end
    # INVARIANTS ASSERTIONS
    # for (id, c) ∈  egraph.classes
    #     # ecdata.nodes = map(n -> canonicalize(egraph.uf, n), ecdata.nodes)
    #     println(id, "=>", c.id)
    #     @assert(id == c.id)
    #     # for an ∈ egraph.analyses
    #     #     if haskey(an, id)
    #     #         @assert an[id] == mapreduce(x -> make(an, x), (x, y) -> join(an, x, y), c.nodes)
    #     #     end
    #     # end

    #     for n ∈ c
    #         println(n)
    #         println("canon = ", canonicalize(egraph, n))
    #         hr = egraph.memo[canonicalize(egraph, n)]
    #         println(hr)
    #         @assert hr == find(egraph, id)
    #     end
    # end
    # display(egraph.classes); println()
    # @show egraph.dirty

end

function repair!(g::EGraph, id::EClassId)
    id = find(g, id)
    ecdata = geteclass(g, id)
    ecdata.id = id
    @debug "repairing " id

    # for (p_enode, p_eclass) ∈ ecdata.parents
    #     clean_enode!(g, p_enode, find(g, p_eclass))
    # end

    new_parents = (length(ecdata.parents) > 30 ? OrderedDict : LittleDict){AbstractENode,EClassId}()

    for (p_enode, p_eclass) ∈ ecdata.parents
        p_enode = canonicalize!(g, p_enode)
        # deduplicate parents
        if haskey(new_parents, p_enode)
            @debug "merging classes" p_eclass (new_parents[p_enode])
            merge!(g, p_eclass, new_parents[p_enode])
        end
        n_id = find(g, p_eclass)
        g.memo[p_enode] = n_id
        new_parents[p_enode] = n_id
    end

    ecdata.parents = collect(new_parents)
    @debug "updated parents " id g.parents[id]

    # ecdata.nodes = map(n -> canonicalize(g.uf, n), ecdata.nodes)

    # Analysis invariant maintenance
    for an ∈ g.analyses
        hasdata(ecdata, an) && modify!(an, g, id)
        # modify!(an, id)
        # id = find(g, id)
        for (p_enode, p_id) ∈ ecdata.parents
            # p_eclass = find(g, p_eclass)
            p_eclass = geteclass(g, p_id)
            if !islazy(an) && !hasdata(p_eclass, an)
                setdata!(p_eclass, an, make(an, g, p_enode))
            end
            if hasdata(p_eclass, an)
                p_data = getdata(p_eclass, an)

                new_data = join(an, p_data, make(an, g, p_enode))
                if new_data != p_data
                    setdata!(p_eclass, an, new_data)
                    push!(g.dirty, p_id)
                end
            end
        end
    end

    unique!(ecdata.nodes)

    # ecdata.nodes = map(n -> canonicalize(g.uf, n), ecdata.nodes)

end

"""
Recursive function that traverses an [`EGraph`](@ref) and
returns a vector of all reachable e-classes from a given e-class id.
"""
function reachable(g::EGraph, id::EClassId)
    id = find(g, id)
    hist = EClassId[id]
    todo = EClassId[id]

    function reachable_node(xn::ENodeTerm)
        x = canonicalize(g, xn)
        for c_id in arguments(x)
            if c_id ∉ hist
                push!(hist, c_id)
                push!(todo, c_id)
            end
        end
    end
    function reachable_node(x::ENodeLiteral)
    end

    while !isempty(todo)
        curr = find(g, pop!(todo))
        for n ∈ g.classes[curr]
            reachable_node(n)
        end
    end

    return hist
end


text:
function generate_proposal_position(	x::Float64,
									y::Float64,
									z::Float64,
									Lx::Float64,
									Ly::Float64,
									Lz::Float64,
									sigma_translation::Float64)

	x_star::Float64 = x + sigma_translation * randn()
	y_star::Float64 = y + sigma_translation * randn()
	z_star::Float64 = z + sigma_translation * randn()

	if x_star < 0.0
		x_star += Lx
	elseif x_star > Lx
		x_star -= Lx
	end

	if y_star < 0.0
		y_star += Ly
	elseif y_star > Ly
		y_star -= Ly
	end

	if z_star < 0.0
		z_star += Lz
	elseif z_star > Lz
		z_star -= Lz
	end

	return (x_star, y_star, z_star)
end

ext: jl

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
# run this to make a precompiled system image with the packages used
# in Econometrics pre-compiled, and with some commonly used functions
# also pre-compiled.
#
# once it's done, start Julia from the Econometrics directory with
# julia --proj -JEconometricsSysimage.so"
#
# delete the .so file to return to normal.
#
# NOTE: this system image does not work with VS Code.

using Pkg, PackageCompiler
Pkg.activate(".")
PackageCompiler.create_sysimage(:Econometrics; sysimage_path="JuliaSysimage.so",
                                                precompile_execution_file="warmup.jl")

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using DynamicalSystems
using Plots
using Random

# Simulate an AR process using DynamicalSystems.jl
X0 = [0.0, 0.1]
p = [0.1, 0.1]
T = 25
@inline @inbounds function AR(X, p, t)
    r = randn()
    dX1 = r + p[1].*X[1] + X[2]
    dX2 = p[2].*X[1]
    return SVector{2}(dX1, dX2)
end
Random.seed!(32)
ds = DiscreteDynamicalSystem(AR, X0, p)

# ----------- Look more closely at the integrator-----------
# function cut_step!(integ, N)
#     #println(integ.u[1])
#     for i in 1:N
#         integ.u = integ.f(integ.u, integ.p, integ.t)
#         integ.t += 1
#     end
#     return
# end

# function cut_trajectory(ds, t, u = ds.u0; dt = 1, Ttr = 0) where {S}
#     dt = round(Int, dt)
#     ti = ds.t0
#     tvec = ti:dt:t+ti
#     L = length(tvec)
#     integ = integrator(ds, u)
#     Ttr != 0 && step!(integ, Ttr)
#     data = zeros(L)
#     data[1] = integ.u[1]
#     for i in 2:L
#         #print(integ)
#         cut_step!(integ, dt)
#         data[i] = integ.u[1]
#     end
#     return data
# end
#data = cut_trajectory(ds, 25; dt=1)
#----------------------------------------
data = trajectory(ds, T; dt=1)

plot(data[:, 1], seriestype=:path)

# Check that this EOM is in the correct form by reproducing:
function f(X, p)
    dX1 = randn() + p[1].*X[1] + X[2]
    dX2 = p[2].*X[1]
    return [dX1, dX2]
end

function integrateAR(X0, p, T)
    X = zeros(T+length(X0), 2)
    X[1, :] = X0
    for t = 1:T
        X[t+length(X0), :] = f(X[t, :], p)
    end
    return X
end
Random.seed!(32)
t = integrateAR(X0, p, T)
plot!(t[:, 1], seriestype=:path)


text:
#=
    Regression tests verify that the sw_driver program successfully ran for each datasize
=#

# include source file of the driver program
include("../src/Csw.jl")

using Test, TOML, .program

input_file_name = joinpath(@__DIR__, "data/inputs/inputs.toml")

configfile = TOML.parsefile(input_file_name)

# Iterate over the data sets in the config file
for (dataset, dataIODict) in configfile
    @testset "$dataset Regression " begin
        main(dataIODict)
        @test true
    end
end

return true


text:
"""
`poro_elasticity_model(;
    n = 980,
    rho = 1e-3,
    alpha = 0.79,
    M = 1/7.80e3,
    kappanu = 633.33,
    eta = 1e-4,
    force_download = false
  )`

This function returns a port-Hamiltonian model of linear poroelasticity in a
bounded Lipschitz domain as described in
    Altmann, Mehrmann, Unger: Port-Hamiltonian Formulations of Poroelastic
    Network Models
# Arguments
- `n`: System dimension (can only be either: 320, 980, or 1805). Default = 980.
- `rho`: density. Default = `1e-3`.
- `alpha`: Biot-Willis fluid-solid coupling coefficient. Default = 0.79.
- `bm`: Biot-Modulus. Default = `1/7.8e3`.
- `kappanu`: Quotient kappa/Nu, where kappa denotes the permeability and nu denotes the fluid viscosity. Default = 633.33.
- `eta`: artificial damping coefficient. Default = `1e-4`.
"""
function poro_elasticity_model(;
    n = 980,
    rho = 1e-3,
    alpha = 0.79,
    bm = 1/7.80e3,
    kappanu = 633.33,
    eta = 1e-4,
    force_download = false
  )
  Y, D, M, K, Bp, Bf, A = load_poro_raw_data(
    n=n,
    force_download = force_download
  )
  Y = rho*sparse(Y)
  D = alpha*sparse(D)
  M = 1/bm .* sparse(M)
  K = kappanu*sparse(K)
  A = sparse(A)
  Bp = Bp'
  Bf = Bf'
  n = size(A, 1);
  m = size(M, 1);
  E = [Y spzeros(n,n+m); spzeros(n,n) A spzeros(n,m); spzeros(m,n+n) M];
  J = [spzeros(n,n) -A D';A spzeros(n,n+m); -D spzeros(m,n+m)];
  R = [spzeros(n,2*n+m); spzeros(n,2*n+m); spzeros(m,2*n) K] + eta*I
  B = [zeros(n,1); Bf; Bp];
  return E, J, R, B
end

function load_poro_raw_data(;
    n = 980,
    force_download = false
  )
  filename = "poro-n$n.mat"
  url = "https://zenodo.org/record/5702554/files/poro-n$n.mat?download=1"
  if n == 980
    md5_hash = hex(0x2961a189be7049ffe2d476b18cb1f678)
  elseif n == 320
    md5_hash = hex(0x97afe8c34f0e9a56bbe86d0a51b7b626)
  elseif n == 1805
    md5_hash = hex(0xc61f6687da9cd26cbf2d880d7d3a9ac9)
  else
    throw(ArgumentError("Model size is either 320, 980, or 1805"))
  end
  download_system_data_if_required(
    filename,
    url,
    md5_hash,
    force_download = force_download
  )
  dd = loadMAT(get_filepath(filename))
  return dd["Y"], dd["D"], dd["M"], dd["K"], dd["Bp"], dd["Bf"], dd["A"]
end

export poro_elasticity_model


text:
# entry
# =====

Base.show(io::IO, res::JETToplevelResult) = print_reports(io, res)
function print_reports(io::IO, res::JETToplevelResult)
    return print_reports(io,
                         get_reports(res),
                         gen_postprocess(res.res.actual2virtual);
                         res.jetconfigs...)
end

Base.show(io::IO, res::JETCallResult) = print_reports(io, res)
function print_reports(io::IO, res::JETCallResult)
    return print_reports(io,
                         get_reports(res);
                         res.jetconfigs...)
end

# entry for test
print_reports(args...; jetconfigs...) =
    return print_reports(stdout::IO, args...; jetconfigs...)

# configuration
# =============

"""
Configurations for report printing.
The configurations below will be active whenever `show`ing [JET's analysis result](@ref analysis-result) within REPL.

---
- `annotate_types::Bool = false` \\
  When set to `true`, annotates types when printing analyzed call stack.

  > Examples:
  * with `annotate_types = false` (default):
    ```julia-repl
    julia> @report_call sum("julia")
    ═════ 2 possible errors found ═════
    ┌ @ reduce.jl:530 Base.#sum#241(Base.pairs(Core.NamedTuple()), #self#, a)
    │┌ @ reduce.jl:530 Base.sum(Base.identity, a)
    ││┌ @ reduce.jl:503 Base.#sum#240(Base.pairs(Core.NamedTuple()), #self#, f, a)
    │││┌ @ reduce.jl:503 Base.mapreduce(f, Base.add_sum, a)
    ││││┌ @ reduce.jl:289 Base.#mapreduce#237(Base.pairs(Core.NamedTuple()), #self#, f, op, itr)
    │││││┌ @ reduce.jl:289 Base.mapfoldl(f, op, itr)
    ││││││┌ @ reduce.jl:162 Base.#mapfoldl#233(Base._InitialValue(), #self#, f, op, itr)
    │││││││┌ @ reduce.jl:162 Base.mapfoldl_impl(f, op, init, itr)
    ││││││││┌ @ reduce.jl:44 Base.foldl_impl(op′, nt, itr′)
    │││││││││┌ @ reduce.jl:48 Base._foldl_impl(op, nt, itr)
    ││││││││││┌ @ reduce.jl:62 op(v, Base.getindex(y, 1))
    │││││││││││┌ @ reduce.jl:81 Base.getproperty(op, :rf)(acc, x)
    ││││││││││││┌ @ reduce.jl:24 Base.+(x, y)
    │││││││││││││ no matching method found for call signature: Base.+(x::Char, y::Char)
    ││││││││││││└────────────────
    │││││││││┌ @ reduce.jl:49 Base.reduce_empty_iter(op, itr)
    ││││││││││┌ @ reduce.jl:356 Base.reduce_empty_iter(op, itr, Base.IteratorEltype(itr))
    │││││││││││┌ @ reduce.jl:357 Base.reduce_empty(op, Base.eltype(itr))
    ││││││││││││┌ @ reduce.jl:330 Base.reduce_empty(Base.getproperty(op, :rf), _)
    │││││││││││││┌ @ reduce.jl:322 Base.reduce_empty(Base.+, _)
    ││││││││││││││┌ @ reduce.jl:313 Base.zero(_)
    │││││││││││││││ no matching method found for call signature: Base.zero(_::Type{Char})
    ││││││││││││││└─────────────────
    Char
    ```
  * with `annotate_types = true`
    ```julia-repl
    julia> @report_call annotate_types = true sum("julia")
    ═════ 2 possible errors found ═════
    ┌ @ reduce.jl:530 Base.#sum#241(Base.pairs(Core.NamedTuple()::NamedTuple{(), Tuple{}})::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, #self#::typeof(sum), a::String)
    │┌ @ reduce.jl:530 Base.sum(Base.identity, a::String)
    ││┌ @ reduce.jl:503 Base.#sum#240(Base.pairs(Core.NamedTuple()::NamedTuple{(), Tuple{}})::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, #self#::typeof(sum), f::typeof(identity), a::String)
    │││┌ @ reduce.jl:503 Base.mapreduce(f::typeof(identity), Base.add_sum, a::String)
    ││││┌ @ reduce.jl:289 Base.#mapreduce#237(Base.pairs(Core.NamedTuple()::NamedTuple{(), Tuple{}})::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, #self#::typeof(mapreduce), f::typeof(identity), op::typeof(Base.add_sum), itr::String)
    │││││┌ @ reduce.jl:289 Base.mapfoldl(f::typeof(identity), op::typeof(Base.add_sum), itr::String)
    ││││││┌ @ reduce.jl:162 Base.#mapfoldl#233(Base._InitialValue()::Base._InitialValue, #self#::typeof(mapfoldl), f::typeof(identity), op::typeof(Base.add_sum), itr::String)
    │││││││┌ @ reduce.jl:162 Base.mapfoldl_impl(f::typeof(identity), op::typeof(Base.add_sum), init::Base._InitialValue, itr::String)
    ││││││││┌ @ reduce.jl:44 Base.foldl_impl(op′::Base.BottomRF{typeof(Base.add_sum)}, nt::Base._InitialValue, itr′::String)
    │││││││││┌ @ reduce.jl:49 Base.reduce_empty_iter(op::Base.BottomRF{typeof(Base.add_sum)}, itr::String)
    ││││││││││┌ @ reduce.jl:356 Base.reduce_empty_iter(op::Base.BottomRF{typeof(Base.add_sum)}, itr::String, Base.IteratorEltype(itr::String)::Base.HasEltype)
    │││││││││││┌ @ reduce.jl:357 Base.reduce_empty(op::Base.BottomRF{typeof(Base.add_sum)}, Base.eltype(itr::String)::Type{Char})
    ││││││││││││┌ @ reduce.jl:330 Base.reduce_empty(Base.getproperty(op::Base.BottomRF{typeof(Base.add_sum)}, :rf::Symbol)::typeof(Base.add_sum), _::Type{Char})
    │││││││││││││┌ @ reduce.jl:322 Base.reduce_empty(Base.+, _::Type{Char})
    ││││││││││││││┌ @ reduce.jl:313 Base.zero(_::Type{Char})
    │││││││││││││││ no matching method found for call signature: Base.zero(_::Type{Char})
    ││││││││││││││└─────────────────
    │││││││││┌ @ reduce.jl:48 Base._foldl_impl(op::Base.BottomRF{typeof(Base.add_sum)}, nt::Base._InitialValue, itr::String)
    ││││││││││┌ @ reduce.jl:62 op::Base.BottomRF{typeof(Base.add_sum)}(v::Char, Base.getindex(y::Tuple{Char, Int64}, 1)::Char)
    │││││││││││┌ @ reduce.jl:81 Base.getproperty(op::Base.BottomRF{typeof(Base.add_sum)}, :rf::Symbol)::typeof(Base.add_sum)(acc::Char, x::Char)
    ││││││││││││┌ @ reduce.jl:24 Base.+(x::Char, y::Char)
    │││││││││││││ no matching method found for call signature: Base.+(x::Char, y::Char)
    ││││││││││││└────────────────
    Char
    ```

  !!! note
      JET always annotates types when printing the error point, e.g. in the example above,
      the error points below are always type-annotated regardless of this configuration:
      - `no matching method found for call signature: Base.zero(_::Type{Char})`
      - `no matching method found for call signature: Base.+(x::Char, y::Char)`
---
- `fullpath::Bool = false` \\
  Controls whether or not expand a file path to full path when printing analyzed call stack.
  Note that paths of Julia's `Base` files will also be expanded when set to `true`.
---
- `print_toplevel_success::Bool = false` \\
  If `true`, prints a message when there is no toplevel errors found.
---
- `print_inference_success::Bool = true` \\
  If `true`, print a message when there is no errors found in abstract interpretation based analysis pass.
---
"""
struct PrintConfig
    print_toplevel_success::Bool
    print_inference_success::Bool
    annotate_types::Bool
    fullpath::Bool
    @jetconfigurable PrintConfig(; print_toplevel_success::Bool  = false,
                                   print_inference_success::Bool = true,
                                   annotate_types::Bool          = false,
                                   fullpath::Bool                = false,
                                   ) =
        return new(print_toplevel_success,
                   print_inference_success,
                   annotate_types,
                   fullpath,
                   )
end

# utility
# =======

const ERROR_COLOR = :light_red
const NOERROR_COLOR = :light_green
# TODO other nicer color scheme ?
const RAIL_COLORS = ( # Julia color + yellow
    :green,
    :magenta,
    :blue,
    :yellow,
)
const N_RAILS = length(RAIL_COLORS)
const LEFT_ROOF  = "═════ "
const RIGHT_ROOF = " ═════"
const HEADER_COLOR = :reverse
const ERROR_SIG_COLOR = :bold
const TYPE_ANNOTATION_COLOR = :light_cyan
const HINT_COLOR = :light_green

pluralize(n::Integer, one::AbstractString, more::AbstractString = string(one, 's')) =
    string(n, ' ', isone(n) ? one : more)

printlnstyled(args...; kwarg...) = printstyled(args..., '\n'; kwarg...)

function print_rails(io, depth)
    for i = 1:depth
        color = RAIL_COLORS[i%N_RAILS+1]
        printstyled(io, '│'; color)
    end
end

function with_bufferring(f, args...)
    buf = IOBuffer()
    io = IOContext(buf, args...)
    f(io)
    return String(take!(buf))
end

# toplevel
# ========

function print_reports(io::IO,
                       reports::Vector{ToplevelErrorReport},
                       @nospecialize(postprocess = identity);
                       jetconfigs...)
    config = PrintConfig(; jetconfigs...)

    n = length(reports)
    if n == 0
        if config.print_toplevel_success
            printlnstyled(io, "No toplevel errors detected"; color = NOERROR_COLOR)
        end
        return 0
    end

    arg = :color => get(io, :color, false)
    with_bufferring(arg) do io
        s = string(pluralize(n, "toplevel error"), " found")
        printlnstyled(io, LEFT_ROOF, s, RIGHT_ROOF; color = HEADER_COLOR)

        color = ERROR_COLOR

        rail = with_bufferring(arg) do io
            printstyled(io, "│ "; color)
        end

        for report in reports
            s = string("┌ @ ", (config.fullpath ? tofullpath : identity)(string(report.file)), ':', report.line, ' ')
            printlnstyled(io, s; color)

            errlines = with_bufferring(arg) do io
                print_report(io, report)
            end |> strip
            println(io, join(string.(rail, split(errlines, '\n')), '\n'))

            s = string('└', '─'^(length(s)-1))
            printlnstyled(io, s; color)
        end
    end |> postprocess |> (x->print(io::IO,x))

    return n
end

# don't show stacktrace for syntax errors
print_report(io, report::SyntaxErrorReport) = showerror(io, report.err)
function print_report(io, report::RecursiveIncludeErrorReport)
    printstyled(io, "ERROR: "; bold = true, color = ERROR_COLOR)
    println(io, "recursive `include` call detected:")
    println(io, " ⚈ duplicated file: ", report.duplicated_file)
    println(io, " ⚈  included files: ", join(report.files, ' '))
end
# TODO: add context information, i.e. during macroexpansion, defining something
print_report(io, report::ActualErrorWrapped) = showerror(io, report.err, report.st)
function print_report(io, report::MissingConcretization)
    printstyled(io, "HINT: "; bold = true, color = HINT_COLOR)
    printlnstyled(io, """
    the following error happened mostly because of the missing concretization of global variables,
    and this could be fixed with the `concretization_patterns` configuration.
    Check https://aviatesk.github.io/JET.jl/dev/config/#JET.ToplevelConfig for the details.
    ---"""; color = HINT_COLOR)

    showerror(io, report.err, report.st)
end

# inference
# =========

function print_reports(io::IO,
                       reports::Vector{InferenceErrorReport},
                       @nospecialize(postprocess = identity);
                       jetconfigs...)
    config = PrintConfig(; jetconfigs...)

    n = length(reports)

    if n == 0
        if config.print_inference_success
            printlnstyled(io, "No errors detected"; color = NOERROR_COLOR)
        end
        return 0
    end

    with_bufferring(:color => get(io, :color, false)) do io
        s = string(pluralize(length(reports), "possible error"), " found")
        printlnstyled(io, LEFT_ROOF, s, RIGHT_ROOF; color = HEADER_COLOR)

        # don't duplicated virtual stack frames for reports from the same toplevel frame
        toplevel_linfo_hash = hash(:dummy)
        wrote_linfos = Set{UInt64}()
        for report in reports
            new_toplevel_linfo_hash = hash(first(report.vst))
            if toplevel_linfo_hash != new_toplevel_linfo_hash
                toplevel_linfo_hash = new_toplevel_linfo_hash
                wrote_linfos = Set{UInt64}()
            end
            print_report(io, report, config, wrote_linfos)
        end
    end |> postprocess |> (x->print(io::IO,x))

    return n
end

# traverse abstract call stack, print frames
function print_report(io, report::InferenceErrorReport, config, wrote_linfos, depth = 1)
    if length(report.vst) == depth # error here
        return print_error_frame(io, report, config, depth)
    end

    frame = report.vst[depth]

    # cache current frame info
    linfo_hash = hash(frame)
    should_print = linfo_hash ∉ wrote_linfos
    push!(wrote_linfos, linfo_hash)

    # print current frame and go into deeper
    should_print && print_frame(io, frame, config, depth, false)
    print_report(io, report, config, wrote_linfos, depth + 1)
end

function print_error_frame(io, report, config, depth)
    frame = report.vst[depth]

    len = print_frame(io, frame, config, depth, true)
    print_rails(io, depth-1)
    print_error_report(io, report)

    print_rails(io, depth-1)
    printlnstyled(io, '└', '─'^len; color = ERROR_COLOR)
end

function print_frame(io, frame, config, depth, is_err)
    print_rails(io, depth-1)

    color = is_err ? ERROR_COLOR : RAIL_COLORS[(depth)%N_RAILS+1]
    s = string("┌ @ ", (config.fullpath ? tofullpath : identity)(string(frame.file)), ':', frame.line)
    printstyled(io, s, ' '; color)
    print_signature(io, frame.sig, config)

    return length(s) # the length of frame info string
end

function print_signature(io, sig::Signature, config; kwargs...)
    for a in sig
        _print_signature(io, a, config; kwargs...)
    end
    println(io)
end
function _print_signature(io, @nospecialize(x), config; kwargs...)
    if isa(x, AbstractChar)
        printstyled(io, x; kwargs...)
    elseif isa(x, AbstractString)
        printstyled(io, x; kwargs...)
    else
        if config.annotate_types
            printstyled(io, "::", string(x); color = TYPE_ANNOTATION_COLOR, kwargs...)
        end
    end
end

# default error report printer
function print_error_report(io, report::InferenceErrorReport)
    printstyled(io, "│ ", report.msg, ": "; color = ERROR_COLOR)
    print_signature(io, report.sig,
                    (; annotate_types = true); # always annotate types for errored signatures
                    bold = true,
                    )
end


text:
#### Types and constructors ####

if !isdefined(:StoredArray)
    const StoredArray = AbstractArray
end

# Plain arrays can be treated as images. Other types will have
# metadata associated, make yours a child of one of the following:
abstract AbstractImage{T,N} <: StoredArray{T,N}         # image with metadata
abstract AbstractImageDirect{T,N} <: AbstractImage{T,N}   # each pixel has own value/color
abstract AbstractImageIndexed{T,N} <: AbstractImage{T,N}  # indexed images (i.e., lookup table)

# converts keyword argument to a dictionary
function kwargs2dict(kwargs)
    d = (ASCIIString=>Any)[]
    for (k, v) in kwargs
        d[string(k)] = v
    end
    return d
end

# Direct image (e.g., grayscale, RGB)
type Image{T,N,A<:StoredArray} <: AbstractImageDirect{T,N}
    data::A
    properties::Dict
end
Image(data::StoredArray, props::Dict) = Image{eltype(data),ndims(data),typeof(data)}(data,props)
Image(data::StoredArray; kwargs...) = Image(data, kwargs2dict(kwargs))

# Indexed image (colormap)
type ImageCmap{T,N,A<:StoredArray,C<:AbstractArray} <: AbstractImageIndexed{T,N}
    data::A
    cmap::C
    properties::Dict
end
ImageCmap(data::StoredArray, cmap::AbstractArray, props::Dict) = ImageCmap{eltype(cmap),ndims(data),typeof(data),typeof(cmap)}(data, cmap, props)
ImageCmap(data::StoredArray, cmap::AbstractArray; kwargs...) = ImageCmap(data, cmap, kwargs2dict(kwargs))

# Convenience constructors
grayim{T}(A::StoredArray{T,2}) = Image(A; colorspace="Gray", spatialorder=["x","y"])
grayim{T}(A::StoredArray{T,3}) = Image(A; colorspace="Gray", spatialorder=["x","y","z"])

function colorim{T}(A::StoredArray{T,3})
    if size(A, 1) == 4 || size(A, 3) == 4
        error("The array looks like a 4-channel color image. Please specify the colorspace explicitly (e.g. \"ARGB\" or \"RGBA\".)")
    end

    colorim(A, "RGB")
end

function colorim{T}(A::StoredArray{T,3}, colorspace)
    if 3 <= size(A, 1) <= 4 && 3 <= size(A, 3) <= 4
        error("Both first and last dimensions are of size 3 or 4; impossible to guess which is for color. Use the Image constructor directly.")
    elseif 3 <= size(A, 1) <= 4  # Image as returned by imread for regular 2D RGB images
        Image(A; colorspace=colorspace, colordim=1, spatialorder=["x","y"])
    elseif 3 <= size(A, 3) <= 4  # "Matlab"-style image, as returned by converT(Array, im).
        Image(A; colorspace=colorspace, colordim=3, spatialorder=["y","x"])
    else
        error("Neither the first nor the last dimension is of size 3. This doesn't look like an RGB image.")
    end
end

# Dispatch-based scaling/clipping/type conversion
abstract ScaleInfo{T}

# An array type for colorized overlays of grayscale images
type Overlay{AT<:(AbstractArray...),N,SIT<:(ScaleInfo...)} <: StoredArray{RGB,N}
    channels::AT   # this holds the grayscale arrays
    colors::Vector{RGB}
    scalei::SIT
    visible::BitVector

    function Overlay(channels::(AbstractArray...), colors, scalei::(ScaleInfo...), visible::BitVector)
        nc = length(channels)
        for i = 1:nc
            if ndims(channels[i]) != N
                error("All arrays must have the same dimensionality")
            end
        end
        sz = size(channels[1])
        for i = 2:nc
            if size(channels[i]) != sz
                error("All arrays must have the same size")
            end
        end
        if length(colors) != nc || length(scalei) != nc || length(visible) != nc
            error("All input must have the same length")
        end
        new(channels, [convert(RGB, c) for c in colors], scalei, visible)
    end
end
Overlay(channels::(AbstractArray...), colors, scalei::(ScaleInfo...), visible = trues(length(channels))) =
    Overlay{typeof(channels),ndims(channels[1]),typeof(scalei)}(channels,colors,scalei,convert(BitVector, visible))

function Overlay(channels::(AbstractArray...), colors, clim = ntuple(length(channels), i->limits(channels[i])))
    n = length(channels)
    for i = 1:n
        if length(clim[i]) != 2
            error("clim must be a 2-vector")
        end
    end
    scalei = ntuple(n, i->scaleminmax(Float64, channels[i], clim[i][1], clim[i][2]))
    Overlay{typeof(channels),ndims(channels[1]),typeof(scalei)}(channels, colors, scalei, trues(n))
end

# Returns the overlay as an image, if possible
function OverlayImage(channels::(AbstractArray...), colors::(ColorValue...), args...)
    ovr = Overlay(channels, colors, args...)
    for i = 1:length(channels)
        if isa(channels[i], AbstractImage)
            prop = copy(properties(channels[i]))
            haskey(prop, "colorspace") && delete!(prop, "colorspace")
            haskey(prop, "limits") && delete!(prop, "limits")
            return Image(ovr, prop)
        end
    end
    colorim(ovr)
end

#### Core operations ####

eltype{T}(img::AbstractImage{T}) = T

size(img::AbstractImage) = size(img.data)
size(img::AbstractImage, i::Integer) = size(img.data, i)
size(img::AbstractImage, dimname::String) = size(img.data, dimindex(img, dimname))

ndims(img::AbstractImage) = ndims(img.data)

strides(img::AbstractImage) = strides(img.data)

copy(img::AbstractImage) = deepcopy(img)

# copy, replacing the data
copy(img::AbstractArray, data::AbstractArray) = data

copy(img::AbstractImageDirect, data::StoredArray) = Image(data, copy(img.properties))

copy(img::AbstractImageIndexed, data::StoredArray) = ImageCmap(data, copy(img.cmap), copy(img.properties))

copy(img::AbstractImageDirect, _data::AbstractImageDirect) = copy(img, data(_data))

# Provide new data but reuse the properties & cmap
share(img::AbstractArray, data::AbstractArray) = data

share(img::AbstractImageDirect, data::StoredArray) = Image(data, img.properties)

share(img::AbstractImageIndexed, data::StoredArray) = ImageCmap(data, img.cmap, img.properties)

# similar
similar(img::AbstractImageDirect) = Image(similar(img.data), copy(img.properties))
similar(img::AbstractImageDirect, ::NTuple{0}) = Image(similar(img.data), copy(img.properties))

similar(img::AbstractImageDirect, dims::Dims) = Image(similar(img.data, dims), copy(img.properties))

similar{T}(img::AbstractImageDirect, ::Type{T}) = Image(similar(img.data, T), copy(img.properties))

similar{T}(img::AbstractImageDirect, ::Type{T}, dims::Dims) = Image(similar(img.data, T, dims), copy(img.properties))

similar(img::AbstractImageIndexed) = ImageCmap(similar(img.data), copy(img.cmap), copy(img.properties))
similar(img::AbstractImageIndexed, ::NTuple{0}) = ImageCmap(similar(img.data), copy(img.cmap), copy(img.properties))

similar(img::AbstractImageIndexed, dims::Dims) = ImageCmap(similar(img.data, dims), copy(img.cmap), copy(img.properties))

similar{T}(img::AbstractImageIndexed, ::Type{T}) = ImageCmap(similar(img.data, T), copy(img.cmap), copy(img.properties))

similar{T}(img::AbstractImageIndexed, ::Type{T}, dims::Dims) = ImageCmap(similar(img.data, T, dims), copy(img.cmap), copy(img.properties))

# copy properties
function copy!(imgdest::AbstractImage, imgsrc::AbstractImage, prop1::ASCIIString, props::ASCIIString...)
    imgdest[prop1] = imgsrc[prop1]
    for p in props
        imgdest[p] = imgsrc[p]
    end
    imgdest
end

# convert
convert{T}(::Type{Image{T}}, img::Image{T}) = img
convert(::Type{Image}, A::AbstractArray) = Image(A, properties(A))
# Convert an indexed image (cmap) to a direct image
function convert{II<:AbstractImageIndexed}(::Type{Image}, img::II)
    if size(img.cmap, 2) == 1
        data = reshape(img.cmap[img.data[:]], size(img.data))
        prop = img.properties
        return Image(data, prop)
    else
        newsz = tuple(size(img.data)...,size(img.cmap,2))
        data = reshape(img.cmap[img.data[:],:], newsz)
        prop = copy(img.properties)
        prop["colordim"] = length(newsz)
        return Image(data, prop)
    end
end
# Convert an Image to an array. We convert the image into the canonical storage order convention for arrays. We restrict this to 2d images because for plain arrays this convention exists only for 2d.
# In other cases---or if you don't want the storage order altered---just grab the .data field and perform whatever manipulations you need directly.
function convert{T}(::Type{Array{T}}, img::AbstractImage)
    assert2d(img)
    dat = convert(Array{T}, data(img))
    # put in canonical storage order
    p = spatialpermutation(spatialorder(Matrix), img)
    p = coords_spatial(img)[p]
    cd = colordim(img)
    if cd > 0
        push!(p, cd)
    end
    if issorted(p)
        return dat
    else
        return permutedims(dat, p)
    end
end
convert(::Type{Array}, img::AbstractImage) = convert(Array{eltype(img)}, img)
# See scaling for conversion of Arrays to Images

# Indexing. In addition to conventional array indexing, support syntax like
#    img["x", 100:400, "t", 32]
# where anything not mentioned by name is taken to include the whole range

# setindex!
setindex!(img::AbstractImage, X, i::Real) = setindex!(img.data, X, i)
setindex!{T<:Real}(img::AbstractImage, X, I::Union(Real,AbstractArray{T})) = setindex!(img.data, X, I)
setindex!{T<:Real}(img::AbstractImage, X, I::Union(Real,AbstractArray{T}), J::Union(Real,AbstractArray{T})) = setindex!(img.data, X, I, J)
setindex!{T<:Real}(img::AbstractImage, X, I::Union(Real,AbstractArray{T}), J::Union(Real,AbstractArray{T}),
                   K::Union(Real,AbstractArray{T})) = setindex!(img.data, X, I, J, K)
setindex!{T<:Real}(img::AbstractImage, X, I::Union(Real,AbstractArray{T}), J::Union(Real,AbstractArray{T}),
                   K::Union(Real,AbstractArray{T}), L::Union(Real,AbstractArray{T})) = setindex!(img.data, X, I, J, K, L)
setindex!{T<:Real}(img::AbstractImage, X, I::Union(Real,AbstractArray{T})...) = setindex!(img.data, X, I...)
setindex!{T<:Real}(img::AbstractImage, X, dimname::String, ind::Union(Real,AbstractArray{T}), nameind...) = setindex!(img.data, X, coords(img, dimname, ind, nameind...)...)

# Adding a new property via setindex!
setindex!(img::AbstractImage, X, propname::String) = setindex!(img.properties, X, propname)

# Delete a property!
delete!(img::AbstractImage, propname::String) = delete!(img.properties, propname)

# getindex, sub, and slice return a value or AbstractArray, not an Image
getindex(img::AbstractImage, i::Real) = getindex(img.data, i)
getindex(img::AbstractImage, I::Union(Real,AbstractVector)) = getindex(img.data, I)
getindex(img::AbstractImage, I::Union(Real,AbstractVector), J::Union(Real,AbstractVector)) = getindex(img.data, I, J)
getindex(img::AbstractImage, I::Union(Real,AbstractVector), J::Union(Real,AbstractVector),
            K::Union(Real,AbstractVector)) = getindex(img.data, I, J, K)
getindex(img::AbstractImage, I::Union(Real,AbstractVector), J::Union(Real,AbstractVector),
            K::Union(Real,AbstractVector), L::Union(Real,AbstractVector)) = getindex(img.data, I, J, K, L)
getindex(img::AbstractImage, I::Union(Real,AbstractVector)...) = getindex(img.data, I...)

# getindex{T<:Real}(img::AbstractImage, dimname::String, ind::Union(Real,AbstractArray{T}), nameind...) = getindex(img.data, coords(img, dimname, ind, nameind...)...)
getindex(img::AbstractImage, dimname::ASCIIString, ind, nameind...) = getindex(img.data, coords(img, dimname, ind, nameind...)...)

getindex(img::AbstractImage, propname::ASCIIString) = getindex(img.properties, propname)

sub(img::AbstractImage, I::RangeIndex...) = sub(img.data, I...)

sub(img::AbstractImage, dimname::ASCIIString, ind::RangeIndex, nameind::RangeIndex...) = sub(img.data, coords(img, dimname, ind, nameind...)...)

slice(img::AbstractImage, I::RangeIndex...) = slice(img.data, I...)

slice(img::AbstractImage, dimname::ASCIIString, ind::RangeIndex, nameind::RangeIndex...) = slice(img.data, coords(img, dimname, ind, nameind...)...)

# getindexim, subim, and sliceim return an Image. The first two share properties, the last requires a copy.
function getindexim{T<:Real}(img::AbstractImage, I::Union(Real,AbstractArray{T})...)
    ret = copy(img, data(img)[I...])
    cd = colordim(img)
    nd = ndims(ret)
    if cd > nd
        ret["colordim"] = 0
        ret["colorspace"] = "Unknown"
    end
    td = timedim(img)
    if td > nd
        ret["timedim"] = 0
    end
    sp = spatialproperties(img)
    for p in sp
        val = ret[p]
        if length(val) > nd
            ret[p] = val[1:nd]
        end
    end
    ret
end

getindexim(img::AbstractImage, dimname::ASCIIString, ind::Union(Real,AbstractArray), nameind...) = getindexim(img, coords(img, dimname, ind, nameind...)...)

subim(img::AbstractImage, I::RangeIndex...) = share(img, sub(img.data, I...))

subim(img::AbstractImage, dimname::ASCIIString, ind::RangeIndex, nameind...) = subim(img, coords(img, dimname, ind, nameind...)...)

function sliceim(img::AbstractImage, I::RangeIndex...)
    dimmap = Array(Int, ndims(img))
    n = 0
    for j = 1:ndims(img)
        if !isa(I[j], Int); n += 1; end;
        dimmap[j] = n
    end
    ret = copy(img, slice(img.data, I...))
    cd = colordim(img)
    if cd > 0
        if isa(I[cd], Int)
            ret.properties["colordim"] = 0
            ret.properties["colorspace"] = "Unknown"
        else
            ret.properties["colordim"] = dimmap[cd]
            if I[cd] != 1:size(img, cd)
                ret.properties["colorspace"] = "channels"
            end
        end
    end
    td = timedim(img)
    if td > 0
        ret.properties["timedim"] = isa(I[td], Int) ? 0 : dimmap[td]
    end
    sp = spatialproperties(img)
    if !isempty(sp)
        c = coords_spatial(img)
        keep = Bool[map(x -> isa(x, AbstractVector), I[c])...]
        if !all(keep)
            for pname in sp
                p = img.properties[pname]
                if isa(p, Vector)
                    ret.properties[pname] = p[keep]
                elseif isa(p, Matrix)
                    ret.properties[pname] = p[keep, keep]
                else
                    error("Do not know how to handle property ", pname)
                end
            end
        end
    end
    ret
end

sliceim(img::AbstractImage, dimname::String, ind::RangeIndex, nameind...) = subim(img, coords(img, dimname, ind, nameind...)...)

sliceim(img::AbstractImage, dimname::String, ind::RangeIndex, nameind...) = sliceim(img, coords(img, dimname, ind, nameind...)...)

subim(img::AbstractImage, I::AbstractVector...) = error("Indexes must be integers or ranges")
sliceim(img::AbstractImage, I::AbstractVector...) = error("Indexes must be integers or ranges")

# Support colon indexes
getindexim(img::AbstractImage, I...) = getindexim(img, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(img,i)) : I[i])...)
subim(img::AbstractImage, I...) = subim(img, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(img,i)) : I[i])...)
sliceim(img::AbstractImage, I...) = sliceim(img, ntuple(length(I), i-> isa(I[i], Colon) ? (1:size(img,i)) : I[i])...)

# We'll frequently want to pull out different 2d slices from the same image, so here's a type and set of functions making that easier.
# We deliberately do not require the user to specify the full list of new slicing/ranging parameters, as often we'll want to change some aspects (e.g., z-slice) but not others (e.g., color coordinates)
type SliceData
    slicedims::(Int...,)
    slicestrides::(Int...,)
    rangedims::(Int...,)

    function SliceData(A::AbstractArray, slicedims::Int...)
        keep = trues(ndims(A))
        for i = 1:length(slicedims)
            keep[slicedims[i]] = false
        end
        s = strides(A)
        new(slicedims, ntuple(length(slicedims), i->s[slicedims[i]]), tuple((1:ndims(A))[keep]...))
    end
end

SliceData(img::AbstractImage, dimname::String, dimnames::String...) = SliceData(img, dimindexes(img, dimname, dimnames...)...)

function _slice(A::AbstractArray, sd::SliceData, I::Int...)
    if length(I) != length(sd.slicedims)
        throw(BoundsError())
    end
    indexes = RangeIndex[1:size(A, i) for i = 1:ndims(A)]
    for i = 1:length(I)
        indexes[sd.slicedims[i]] = I[i]
    end
    indexes
end

function slice(A::AbstractArray, sd::SliceData, I::Int...)
    indexes = _slice(A, sd, I...)
    slice(A, indexes...)
end

function sliceim(img::AbstractImage, sd::SliceData, I::Int...)
    indexes = _slice(img, sd, I...)
    sliceim(img, indexes...)
end

function first_index(A::SubArray, sd::SliceData)
    newfirst = 1
    for i = 1:length(sd.slicedims)
        newfirst += (A.indexes[sd.slicedims[i]]-1)*sd.slicestrides[i]
    end
    for i = 1:length(sd.rangedims)
        newfirst += (A.indexes[sd.rangedims[i]][1]-1)*A.strides[i]
    end
    newfirst
end

function reslice!(A::SubArray, sd::SliceData, I::Int...)
    indexes = RangeIndex[A.indexes...]
    for i = 1:length(I)
        indexes[sd.slicedims[i]] = I[i]
    end
    A.indexes = tuple(indexes...)
    A.first_index = first_index(A, sd)
    A
end

function reslice!(img::AbstractImage, sd::SliceData, I::Int...)
    reslice!(img.data, sd, I...)
    img
end

function rerange!(A::SubArray, sd::SliceData, I::(RangeIndex...,))
    indexes = RangeIndex[A.indexes...]
    for i = 1:length(I)
        indexes[sd.rangedims[i]] = I[i]
    end
    A.indexes = tuple(indexes...)
    A.first_index = first_index(A, sd)
    A
end

function rerange!(img::AbstractImage, sd::SliceData, I::(RangeIndex...,))
    rerange!(img.data, sd, I...)
    img
end

const emptyset = Set()
function showim(io::IO, img::AbstractImageDirect)
    IT = typeof(img)
    print(io, colorspace(img), " ", IT.name, " with:\n  data: ", summary(img.data), "\n  properties:")
    showdictlines(io, img.properties, get(img, "suppress", emptyset))
end
function showim(io::IO, img::AbstractImageIndexed)
    IT = typeof(img)
    print(io, colorspace(img), " ", IT.name, " with:\n  data: ", summary(img.data), "\n  cmap: ", summary(img.cmap), "\n  properties:")
    showdictlines(io, img.properties, get(img, "suppress", emptyset))
end
show(io::IO, img::AbstractImageDirect) = showim(io, img)
writemime(io::IO, ::MIME"text/plain", img::AbstractImageDirect) = showim(io, img)
show(io::IO, img::AbstractImageIndexed) = showim(io, img)
writemime(io::IO, ::MIME"text/plain", img::AbstractImageIndexed) = showim(io, img)

data(img::AbstractArray) = img
data(img::AbstractImage) = img.data

minimum(img::AbstractImageDirect) = minimum(img.data)
maximum(img::AbstractImageDirect) = maximum(img.data)
# min/max deliberately not defined for AbstractImageIndexed

# Overlays
eltype{T}(o::Overlay{T}) = T
ndims{T,N}(o::Overlay{T,N}) = N
ndims{T,N}(::Type{Overlay{T,N}}) = N

size(o::Overlay) = size(o.channels[1])
size(o::Overlay, i::Integer) = size(o.channels[1], i)

show(io::IO, o::Overlay) = print(io, summary(o), " with colors ", o.colors)

function squeeze(img::AbstractImage, dims)
    imgret = copy(img, squeeze(data(img), dims))
    td = timedim(img)
    if td > 0 && in(td, dims)
        imgret["timedim"] = 0
    end
    cd = colordim(img)
    if cd > 0 && in(cd, dims)
        imgret["colordim"] = 0
    end
    imgret
end

#### Properties ####

# Generic programming with images uses properties to obtain information. The strategy is to define a particular property name, and then write an accessor function of the same name. The accessor function provides default behavior for plain arrays and when the property is not defined. Alternatively, use get(img, "propname", default) or haskey(img, "propname") to define your own default behavior.

# You can define whatever properties you want. Here is a list of properties used
# in some algorithms:
#   colorspace: "RGB", "ARGB", "Gray", "Binary", "RGB24", "Lab", "HSV", etc.
#   colordim: the array dimension used to store color information, or 0 if there
#     is no dimension corresponding to color
#   timedim: the array dimension used for time (i.e., sequence), or 0 for single images
#   limits: (minvalue,maxvalue) for this type of image (e.g., (0,255) for Uint8
#     images, even if pixels do not reach these values)
#   pixelspacing: the spacing between adjacent pixels along spatial dimensions
#   spacedirections: the direction of each array axis in physical space (a vector-of-vectors, one per dimension)
#   spatialorder: a string naming each spatial dimension, in the storage order of
#     the data array. Names can be arbitrary, but the choices "x" and "y" have special
#     meaning (horizontal and vertical, respectively, irrespective of storage order).
#     If supplied, you must have one entry per spatial dimension.

properties(img::AbstractArray) = [
    "colorspace" => colorspace(img),
    "colordim" => colordim(img),
    "timedim" => timedim(img),
    "limits" => limits(img),
    "pixelspacing" => pixelspacing(img),
    "spatialorder" => spatialorder(img)]
properties(img::AbstractImage) = img.properties

haskey(a::AbstractArray, k::String) = false
haskey(img::AbstractImage, k::String) = haskey(img.properties, k)

get(img::AbstractArray, k::String, default) = default
get(img::AbstractImage, k::String, default) = get(img.properties, k, default)

# So that defaults don't have to be evaluated unless they are needed, we also define a @get macro (thanks Toivo Hennington):
macro get(img, k, default)
    quote
        img, k = $(esc(img)), $(esc(k))
        local val
        if !isa(img, AbstractImage)
            val = $(esc(default))
        else
            index = Base.ht_keyindex(img.properties, k)
            val = (index > 0) ? img.properties.vals[index] : $(esc(default))
        end
        val
    end
end

# Using plain arrays, we have to make all sorts of guesses about colorspace and storage order. This can be a big problem for three-dimensional images, image sequences, cameras with more than 16-bits, etc. In such cases use an AbstractImage type.

# Here are the two most important assumptions (see also colorspace below):
defaultarraycolordim = 3
# defaults for plain arrays ("vertical-major")
const yx = ["y", "x"]
# order used in Cairo & most image file formats (with color as the very first dimension)
const xy = ["x", "y"]
spatialorder(::Type{Matrix}) = yx
spatialorder(img::AbstractArray) = (sdims(img) == 2) ? spatialorder(Matrix) : error("Wrong number of spatial dimensions for plain Array, use an AbstractImage type")

isdirect(img::AbstractArray) = true
isdirect(img::AbstractImageDirect) = true
isdirect(img::AbstractImageIndexed) = false

colorspace{C<:ColorValue}(img::AbstractMatrix{C}) = string(C)
colorspace{C<:ColorValue}(img::AbstractArray{C,3}) = string(C)
colorspace{C<:ColorValue}(img::AbstractImage{C}) = string(C)
colorspace(img::AbstractMatrix{Bool}) = "Binary"
colorspace(img::AbstractArray{Bool}) = "Binary"
colorspace(img::AbstractArray{Bool,3}) = "Binary"
colorspace{T<:Union(Int32,Uint32)}(img::AbstractMatrix{T}) = "RGB24"
colorspace(img::AbstractMatrix) = "Gray"
colorspace{T}(img::AbstractArray{T,3}) = (size(img, defaultarraycolordim) == 3) ? "RGB" : error("Cannot infer colorspace of Array, use an AbstractImage type")
# colorspace{T}(img::AbstractArray{T,3}) = (size(img, defaultarraycolordim) == 3) ? "RGB" : 0
colorspace(img::AbstractImage{Bool}) = "Binary"
colorspace{T,N,A<:AbstractArray,C<:ColorValue}(img::ImageCmap{T,N,A,Array{C,1}}) = string(C)
colorspace(img::AbstractImageIndexed) = @get img "colorspace" csinfer(eltype(img.cmap))
colorspace{T}(img::AbstractImageIndexed{T,2}) = @get img "colorspace" csinfer(eltype(img.cmap))
csinfer{C<:ColorValue}(::Type{C}) = string(C)
csinfer(C) = "Unknown"
colorspace(img::AbstractImage) = get(img.properties, "colorspace", "Unknown")

colordim{C<:ColorValue}(img::AbstractVector{C}) = 0
colordim{C<:ColorValue}(img::AbstractMatrix{C}) = 0
colordim{C<:ColorValue}(img::AbstractArray{C,3}) = 0
colordim{C<:ColorValue}(img::AbstractImage{C}) = 0
colordim(img::AbstractMatrix) = 0
colordim{T}(img::AbstractImageDirect{T,3}) = get(img, "colordim", 0)
colordim{T}(img::AbstractArray{T,3}) = (size(img, defaultarraycolordim) == 3) ? 3 : 0
colordim(img::AbstractImageDirect) = get(img, "colordim", 0)
colordim(img::AbstractImageIndexed) = 0

timedim(img) = get(img, "timedim", 0)

limits(img::AbstractArray{Bool}) = 0,1
limits{T<:Integer}(img::AbstractArray{T}) = typemin(T), typemax(T)
limits{T<:FloatingPoint}(img::AbstractArray{T}) = zero(T), one(T)
limits(img::AbstractImage{Bool}) = 0,1
limits{T<:ColorValue}(img::AbstractArray{T}) = 0,1
limits{T}(img::AbstractImageDirect{T}) = get(img, "limits", (typemin(T), typemax(T)))
limits(img::AbstractImageIndexed) = @get img "limits" (minimum(img.cmap), maximum(img.cmap))

pixelspacing{T}(img::AbstractArray{T,3}) = (size(img, defaultarraycolordim) == 3) ? [1.0,1.0] : error("Cannot infer pixelspacing of Array, use an AbstractImage type")
pixelspacing(img::AbstractMatrix) = [1.0,1.0]
pixelspacing{T}(img::AbstractImage{T}) = @get img "pixelspacing" _pixelspacing(img)
function _pixelspacing(img::AbstractImage)
    if haskey(img, "spacedirections")
        sd = img["spacedirections"]
        return [maximum(abs(sd[i])) for i = 1:length(sd)]
    end
    ones(sdims(img))
end

spacedirections(img::AbstractArray) = @get img "spacedirections" _spacedirections(img)
function _spacedirections(img::AbstractArray)
    ps = pixelspacing(img)
    T = eltype(ps)
    nd = length(ps)
    Vector{T}[(tmp = zeros(T, nd); tmp[i] = ps[i]; tmp) for i = 1:nd]
end

spatialorder(img::AbstractImage) = @get img "spatialorder" _spatialorder(img)
_spatialorder(img::AbstractImage) = (sdims(img) == 2) ? spatialorder(Matrix) : error("Cannot guess default spatial order for ", sdims(img), "-dimensional images")

# This is mostly for user information---in code it's generally better to use spatialorder, colordim, and timedim directly
function storageorder(img::AbstractArray)
    so = Array(ASCIIString, ndims(img))
    so[coords_spatial(img)] = spatialorder(img)
    cd = colordim(img)
    if cd != 0
        so[cd] = "color"
    end
    td = timedim(img)
    if td != 0
        so[td] = "t"
    end
    so
end

# number of spatial dimensions in the image
sdims(img) = ndims(img) - (colordim(img) != 0) - (timedim(img) != 0)

# number of time slices
function nimages(img)
    sd = timedim(img)
    if sd > 0
        return size(img, sd)
    else
        return 1
    end
end

# number of array elements used for each pixel/voxel
function ncolorelem(img)
    cd = colordim(img)
    return cd > 0 ? size(img, cd) : 1
end

# indices of spatial coordinates
function coords_spatial(img)
    nd = ndims(img)
    cd = colordim(img)
    td = timedim(img)
    if cd > nd || td > nd
        error("Properties are inconsistent with the array dimensionality")
    end
    ind = [1:nd]
    if cd > td
        splice!(ind, cd)
        if td > 0
            splice!(ind, td)
        end
    elseif td > cd
        splice!(ind, td)
        if cd > 0
            splice!(ind, cd)
        end
    end
    ind
end

# size of the spatial grid
function size_spatial(img)
    sz = size(img)
    sz[coords_spatial(img)]
end

#### Utilities for writing "simple algorithms" safely ####
# If you don't feel like supporting multiple representations, call these

# Two-dimensional images
function assert2d(img::AbstractArray)
    if sdims(img) != 2
        error("Only two-dimensional images are supported")
    end
    if timedim(img) != 0
        error("Image sequences are not supported")
    end
end

# "Scalar color", either grayscale, RGB24, or an immutable type
function assert_scalar_color(img::AbstractArray)
    if colordim(img) != 0
        error("Only 'scalar color' is supported")
    end
end

# Check that the time dimension, if present, is last
function assert_timedim_last(img::AbstractArray)
    if 0 < timedim(img) < ndims(img)
        error("Time dimension is not last")
    end
end

# Spatial storage order
isyfirst(img::AbstractArray) = spatialorder(img)[1] == "y"
function assert_yfirst(img)
    if !isyfirst(img)
        error("Image must have y as its first dimension")
    end
end
isxfirst(img::AbstractArray) = spatialorder(img)[1] == "x"
function assert_xfirst(img::AbstractArray)
    if !isxfirst(img)
        error("Image must have x as its first dimension")
    end
end

#### Permutations over dimensions ####

# width and height, translating "x" and "y" spatialorder into horizontal and vertical, respectively
function widthheight(img::AbstractArray, p)
    c = coords_spatial(img)
    size(img, c[p[1]]), size(img, c[p[2]])
end
widthheight(img::AbstractArray) = widthheight(img, spatialpermutation(xy, img))

width(img::AbstractArray) = widthheight(img)[1]
height(img::AbstractArray) = widthheight(img)[2]

# Calculate the permutation needed to put the spatial dimensions into a specified order
spatialpermutation(to, img::AbstractArray) = default_permutation(to, spatialorder(img))
function spatialpermutation(to, img::AbstractImage)
    so = spatialorder(img)
    if so != nothing
        return default_permutation(to, so)
    else
        if sdims(img) != 2
            error("Cannot guess default spatialorder when there are more than 2 spatial dimensions")
        end
        return default_permutation(to, yx)
    end
end

# Permute the dimensions of an image, also permuting the relevant properties. If you have non-default properties that are vectors or matrices relative to spatial dimensions, include their names in the list of spatialprops.
permutedims(img::AbstractImage, p::(), spatialprops::Vector = spatialproperties(img)) = img

function permutedims(img::AbstractImage, p::Union(Vector{Int}, (Int...)), spatialprops::Vector = spatialproperties(img))
    if length(p) != ndims(img)
        error("The permutation must have length equal to the number of dimensions")
    end
    if issorted(p) && length(p) == ndims(img)
        return img   # should we return a copy?
    end
    ip = invperm(to_vector(p))
    cd = colordim(img)
    sd = timedim(img)
    ret = copy(img, permutedims(img.data, p))
    if cd > 0
        ret.properties["colordim"] = ip[cd]
        p = setdiff(p, cd)
    end
    if sd > 0
        ret.properties["timedim"] = ip[sd]
        p = setdiff(p, sd)
    end
    if !isempty(spatialprops)
        ip = sortperm(p)
        for prop in spatialprops
            a = img.properties[prop]
            if isa(a, AbstractVector)
                ret.properties[prop] = a[ip]
            elseif isa(a, AbstractMatrix) && size(a,1) == size(a,2)
                ret.properties[prop] = a[ip,ip]
            else
                error("Do not know how to handle property ", prop)
            end
        end
    end
    ret
end

permutedims{S<:String}(img::AbstractImage, pstr::Union(Vector{S}, (S...)), spatialprops::Vector = spatialproperties(img)) = permutedims(img, dimindexes(img, pstr...), spatialprops)

# Define the transpose of a 2d image
function ctranspose(img::AbstractImage)
    assert2d(img)
    s = coords_spatial(img)
    p = [1:ndims(img)]
    p[s] = s[2:-1:1]
    permutedims(img, p)
end

# Default list of spatial properties possessed by an image
function spatialproperties(img::AbstractImage)
    if haskey(img, "spatialproperties")
        return img.properties["spatialproperties"]
    end
    spatialprops = ASCIIString[]
    if haskey(img, "spatialorder")
        push!(spatialprops, "spatialorder")
    end
    if haskey(img, "pixelspacing")
        push!(spatialprops, "pixelspacing")
    end
    if haskey(img, "spacedirections")
        push!(spatialprops, "spacedirections")
    end
    spatialprops
end
spatialproperties(img::AbstractVector) = ASCIIString[]  # these are not mutable

#### Low-level utilities ####
function permutation(to, from)
    n = length(to)
    nf = length(from)
    d = Dict(tuple(from...), tuple([1:length(from)]...))
    ind = Array(Int, max(n, nf))
    for i = 1:n
        ind[i] = get(d, to[i], 0)
    end
    ind[n+1:nf] = n+1:nf
    ind
end

function default_permutation(to, from)
    p = permutation(to, from)
    pzero = p .== 0
    if any(pzero)
        p[pzero] = setdiff(1:length(to), p)
    end
    p
end

function showdictlines(io::IO, dict::Dict, suppress::Set)
    for (k, v) in dict
        if k == "suppress"
            continue
        end
        if !in(k, suppress)
            print(io, "\n    ", k, ": ")
            printdictval(io, v)
        else
            print(io, "\n    ", k, ": <suppressed>")
        end
    end
end

printdictval(io::IO, v) = print(io, v)
function printdictval(io::IO, v::Vector)
    for i = 1:length(v)
        print(io, " ", v[i])
    end
end

# Support indexing via
#    img["t", 32, "x", 100:400]
# where anything not mentioned by name is assumed to include the whole range
function coords(img::AbstractImage, dimname::ASCIIString, ind, nameind...)
    c = Any[1:d for d in size(img)]
    so = spatialorder(img)
    c[require_dimindex(img, dimname, so)] = ind
    for i = 1:2:length(nameind)
        c[require_dimindex(img, nameind[i], so)] = nameind[i+1]
    end
    tuple(c...)
end

# Use keyword arguments
# e.g. coord(x=1:100, y=1:50)
function coords(img::AbstractImage; kwargs...)
    c = Any[1:d for d in size(img)]
    so = spatialorder(img)
    for (k, v) in kwargs
        c[require_dimindex(img, string(k), so)] = v
    end
    tuple(c...)
end

function dimindex(img::AbstractImage, dimname::ASCIIString, so = spatialorder(img))
    n::Int = 0
    if dimname == "color"
        n = colordim(img)
    elseif dimname == "t"
        n = timedim(img)
    else
        cd = colordim(img)
        td = timedim(img)
        j = 1
        tn = 1
        while j <= length(so)
            while tn == cd || tn == td
                tn += 1
            end
            if dimname == so[j]
                n = tn
                break
            end
            tn += 1
            j += 1
        end
    end
    n
end

require_dimindex(img::AbstractImage, dimname, so) = (di = dimindex(img, dimname, so); di > 0 || error("No dimension called ", dimname); di)

dimindexes(img::AbstractImage, dimnames::String...) = Int[dimindex(img, nam, spatialorder(img)) for nam in dimnames]

to_vector(v::AbstractVector) = v
to_vector(v::Tuple) = [v...]


text:
export BatchSampleMsg,
    BatchDataMsg,
    FetchParamMsg,
    LoadParamMsg,
    InsertTrajectoryMsg,
    Trainer,
    TrajectoryManager,
    Worker,
    WorkerProxy,
    Orchestrator,
    InsertSampleLoadRateLimiter

using Flux

#####
# Messages
#####

struct BatchSampleMsg <: AbstractMessage
    from::RemoteChannel
end

struct BatchDataMsg{D} <: AbstractMessage
    data::D
end

Base.@kwdef struct FetchParamMsg <: AbstractMessage
    from::RemoteChannel = self()
end

struct LoadParamMsg{P} <: AbstractMessage
    data::P
end

struct InsertTrajectoryMsg{B} <: AbstractMessage
    data::B
end

#####
# Trainer
#####

"""
    Trainer(;policy, sealer=deepcopy)

A wrapper around `AbstractPolicy`, the `sealer` is used to create an immutable
object from the inner parameters of `policy` when received a `FetchParamMsg`.
"""
Base.@kwdef struct Trainer{P,S}
    policy::P
    sealer::S = deepcopy
end

Trainer(p) = Trainer(;policy=p)

function (trainer::Trainer)(msg::BatchDataMsg)
    update!(trainer.policy, msg.data)
end

function (trainer::Trainer)(msg::FetchParamMsg)
    ps = trainer.sealer(params(trainer.policy))
    put!(msg.from, LoadParamMsg(ps))
end

#####
# TrajectoryManager
#####

Base.@kwdef mutable struct TrajectoryManager{T,S,I}
    trajectory::T
    sampler::S
    inserter::I
end

(t::TrajectoryManager)(msg::InsertTrajectoryMsg) = push!(t.trajectory, msg.data, t.inserter)

function (t::TrajectoryManager)(msg::BatchSampleMsg)
    s = sample(t.trajectory, t.sampler)  # !!! sample must ensure no sharing data
    put!(msg.from, BatchDataMsg(s))
end

#####
# Worker
#####

"""
    Worker(()->ex::Experiment)
"""
mutable struct Worker
    init::Any
    experiment::Experiment
    task::Task
    Worker(f) = new(f)
end

function (w::Worker)(msg::StartMsg)
    w.experiment = w.init(msg.args...;msg.kwargs...)
    w.task = Threads.@spawn run(w.experiment)
end

function (w::Worker)(msg::StopMsg)
    msg(w.experiment.stop_condition)
    wait(w.task)
end

(w::Worker)(msg::LoadParamMsg) = msg(w.experiment.hook)

#####
# WorkerProxy
#####

mutable struct WorkerProxy
    is_fetch_msg_sent::Ref{Bool}
    workers::Vector{<:RemoteChannel}
    target::RemoteChannel
    WorkerProxy(workers) = new(Ref(false), workers)
end

function (wp::WorkerProxy)(msg::StartMsg)
    wp.target = msg.args[1]
    for w in wp.workers
        put!(w, StartMsg(self()))
    end
end

(wp::WorkerProxy)(msg::InsertTrajectoryMsg) = put!(wp.target, msg)

function (wp::WorkerProxy)(::FetchParamMsg)
    if !wp.is_fetch_msg_sent[]
        put!(wp.target, FetchParamMsg(self()))
        wp.is_fetch_msg_sent[] = true
    end
end

function (wp::WorkerProxy)(msg::LoadParamMsg)
    for w in wp.workers
        put!(w, msg)
    end
    wp.is_fetch_msg_sent[] = false
end

#####
# Orchestrator
#####

Base.@kwdef mutable struct InsertSampleLoadRateLimiter
    min_insert_before_sampling::Int = 100
    n_sample::Int = 0
    n_insert::Int = 0
    n_load::Int = 0
    sample_insert_ratio::Int = 1
    sample_load_ratio::Int = 1
end

Base.@kwdef struct Orchestrator
    trainer::RemoteChannel
    trajectory_proxy::RemoteChannel
    worker::RemoteChannel{Channel{Any}}
    limiter::InsertSampleLoadRateLimiter = InsertSampleLoadRateLimiter()
end

function (orc::Orchestrator)(msg::StartMsg)
    put!(orc.worker, StartMsg(self()))
end

function (orc::Orchestrator)(msg::InsertTrajectoryMsg)
    L = orc.limiter
    put!(orc.trajectory_proxy, msg)
    L.n_insert += 1
    if L.n_insert > L.min_size_to_sample
        for i in 1:L.sample_insert_ratio
            put!(orc.trajectory_proxy, BatchSampleMsg(orc.trainer))
            L.n_sample += 1
            if L.n_sample == (L.n_load + 1) * L.sample_load_ratio
                put!(orc.trajectory_proxy, ProxyMsg(to=orc.trainer, msg=FetchParamMsg(orc.worker)))
                L.n_load += 1
            end
        end
    end
end


text:
function EstimatedTroposphericModel(arg0::MappingFunction, arg1::jdouble)
    return EstimatedTroposphericModel((MappingFunction, jdouble), arg0, arg1)
end

function EstimatedTroposphericModel(arg0::jdouble, arg1::jdouble, arg2::MappingFunction, arg3::jdouble)
    return EstimatedTroposphericModel((jdouble, jdouble, MappingFunction, jdouble), arg0, arg1, arg2, arg3)
end

function compute_zenith_delay(obj::EstimatedTroposphericModel, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::FieldAbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2)
end

function compute_zenith_delay(obj::EstimatedTroposphericModel, arg0::jdouble, arg1::Vector{jdouble}, arg2::AbsoluteDate)
    return jcall(obj, "computeZenithDelay", Vector{jdouble}, (jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_parameters_drivers(obj::EstimatedTroposphericModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function mapping_factors(obj::EstimatedTroposphericModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{RealFieldElement}, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function mapping_factors(obj::EstimatedTroposphericModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "mappingFactors", Vector{jdouble}, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function path_delay(obj::EstimatedTroposphericModel, arg0::RealFieldElement, arg1::RealFieldElement, arg2::Vector{RealFieldElement}, arg3::FieldAbsoluteDate)
    return jcall(obj, "pathDelay", RealFieldElement, (RealFieldElement, RealFieldElement, Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function path_delay(obj::EstimatedTroposphericModel, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::AbsoluteDate)
    return jcall(obj, "pathDelay", jdouble, (jdouble, jdouble, Vector{jdouble}, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using SAMTools
using AvailablePotentialEnergyFramework
using DelimitedFiles

file_control = "/global/homes/a/aramreye/SAM_model/SAM6.10.6/exp/f5e-4_2km_1000km_control/init_snd/snd300"
file_ensemble1 = "/global/homes/a/aramreye/SAM_model/SAM6.10.6/exp/f5e-4_2km_1000km_control_ensemble1/init_snd/snd300"
file_ensemble2 = "/global/homes/a/aramreye/SAM_model/SAM6.10.6/exp/f5e-4_2km_1000km_control_ensemble2/init_snd/snd300"

profile_control = readdlm(file_control;header=true,skipstart=1)[1][1:56,2:4]
profile_ensemble1= readdlm(file_ensemble1;header=true,skipstart=1)[1][1:56,2:4]
profile_ensemble2 = readdlm(file_ensemble2;header=true,skipstart=1)[1][1:56,2:4]

profile_enseble1_modified = copy(profile_ensemble1)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(chain::Chain)
  seriestype --> :path
  seriescolor --> :black
  label --> "chain"

  xs = coordinates.(vertices(chain))

  Tuple.(xs)
end


text:
using CausalTrees
using Plots
import Statistics

# Functions to generate data ##################################################
"""
Generate RCT like data, with heterogeneity in the treatment effect
n: number of samples (total)
pnoise: number of noise variable
"""
function gen_trial1(n::Int,pnoise::Int=0)
    # covariates
    x = randn((n,3+pnoise))
    # treatment group
    w = rand([0,1],n)
    # The treatment effect: E(Y1-Y0|X)
    tau =  2 .+ 4.35 .* x[:,1] + 0.8 .* x[:,1].^2 + 0.3 .* x[:,1].^3 + 1.4 .* x[:,2] .* x[:,3]
    # The potential outcomes
    y0 = 1 .* x[:,1] + x[:,2] + randn(n)
    y1 = 1 .* x[:,1] + x[:,2] + tau + randn(n)
    truth = y1 .- y0
    # observed outcome
    y = (1 .- w) .* y0 + w .* y1
    y = reshape(y,n)
    (x,w,y,tau)
end

"""
logistic
"""
logistic(α) = exp(α) / (1 + exp(α))

"""
Generate RCT like data, with heterogeneity in the treatment effect
n: number of samples (total)
pnoise: number of noise variable
"""
function gen_trial2(n::Int,pnoise::Int=0)
    # covariates
    x = randn((n,3+pnoise))
    # treatment group
    w = rand([0,1],n)
    # part of treatment effect
    α = 2 .* x[:,1] + 0.5 .* x[:,2] + 0.7 .* w
    τ = 0.4 .* x[:,1]
    y0 = logistic.(α)
    y1 = logistic.(α) .+ logistic.(τ + 0.1 .* randn(n))
    truth = y1 .- y0
    # observed outcome
    ypr = (1 .- w) .* y0 + w .* y1
    y = zeros(n)
    rr = rand(n)
    y[ypr .> rr] .= 1.0
    y = reshape(y,n)
    (x,w,y,ypr)
end

# TEST 1 ######################################################################
# generate data
x,w,y,tau = gen_trial1(200)
# plot
plotly()
scatter(x[:,1],y,colour=w .+ 1,label="σ",legend=:bottomright)
# test split finding algorithm: split #1
@time split_i,thres,loss = CausalTrees.evaluate_possible_splits(x[:,1],w,y,10)
#@time pos_thres = CausalTrees.compute_possible_splits(x[:,1],10)
vline!([thres], label = L"split #1")
# test split finding algorithm: split #2
msk = x[:,1] .< thres
@time split_i,thres,loss = CausalTrees.evaluate_possible_splits(x[msk,1],w[msk],y[msk],10)
vline!([thres], label = "split #2")

# TEST 2 #######################################################################
# generate data
x,w,y,tau = gen_trial1(200)
# plot
scatter(x[:,1],tau,color= w,label= "")
# fit tree and predict
@time t1 = causal_tree(x,w,y,mse_tau,min_leaf_size=10)
yhat = predict(t1,x)
Plots.scatter!(x[:,1],yhat,label="",markershape=:hline,markercolor=:red)
# plot the tree
@time CausalTrees.plot(t1,leaf_size=5,split_xsize=2,split_ysize=1)
# print the results
CausalTrees.print(t1)

# TEST 3 #######################################################################
# generate data
x,w,y,tau = gen_trial1(200)
# plot
scatter(x[:,1],tau,color= w,label="")
# fit tree and predict
@time t1 = causal_tree(x,w,y,mse_tau,min_leaf_size=20)
@time t2 = causal_tree(x,w,y,mse_tau,min_leaf_size=20,honesty=true)
yhat1 = predict(t1,x)
yhat2 = predict(t2,x)
scatter!(x[:,1],yhat1,label="not honest")
scatter!(x[:,1],yhat2,label="honest")

# TEST 4 #######################################################################
# generate data
x,w,y,tau = gen_trial1(400)
# plot
scatter(x[:,1],tau,colour=w,label="",legend=:bottomright)
hline!([0],label="null effect")
# fit random forest and predict
@time rf1 = causal_forest(x,w,y,mse_tau,n_trees=500,row_subsamp_p=0.8,
    min_leaf_size=20)
yhat = predict(rf1,x)
scatter!(x[:,1],yhat,label="causal forest")

# TEST 4 #######################################################################
import Random
Random.seed!(1234)
# generate data (as above with larger dataset)
x,w,y,tau = gen_trial1(1000)
# plot
scatter(x[:,1],tau,colour=w,label="",legend=:bottomright,markersize=2)
hline!([0],label="null effect")
# fit random forest and predict
@time t1 = causal_tree(x,w,y,mse_tau,min_leaf_size=10,max_depth=5)
@time rf1 = causal_forest(x,w,y,mse_tau,n_trees=100,row_subsamp_p=0.5,min_leaf_size=20,max_depth=10)
tauhat_rf = predict(rf1,x)
tauhat_t = predict(t1,x)
scatter!(x[:,1],tauhat_rf,label="causal forest",colour=:lightblue)
scatter!(x[:,1],tauhat_t,label="causal tree",colour=:pink)
ate_ = CausalTrees.ate(y,w)
hline!([ate_],label="ate")

# TEST 5 #######################################################################
# binary data
n = 5000
# covariates
x = randn((n,5))
# treatment group
w = rand([0,1],n)
# part of treatment effect
α = 0.6 .* x[:,1] + 0.7 .* x[:,2]
τ = 0.4 .+ 0.4 .* x[:,1]
y0 = logistic.(α + 0.01 .* randn(n))
y1 = logistic.(α + τ + 0.01 .* randn(n))
scatter(x[:,1],y1 .- y0)
ypr = (1 .- w) .* y0 + w .* y1
y = zeros(n)
rr = rand(n)
y[ypr .> rr] .= 1.0
y = reshape(y,n)

println("True: ",Statistics.mean(y1 .- y0))

println("Est: ",Statistics.mean(y[w .== 1]) - Statistics.mean(y[w .== 0]))
# you need to add w to x for this to work!!!!!!

# plots
#histogram(tau)
scatter(x[:,1],y1 .- y0,colour=w,label="",legend=:bottomright,markersize=5)
# fit random forest and predict
@time t1 = causal_tree(x,w,y,mse_tau,min_leaf_size=40,max_depth=3,honesty=true)
tauhat_t = predict(t1,x)
scatter!(x[:,1],tauhat_t,label="causal tree",colour=:pink)
ate_ = CausalTrees.ate(y,w)
println(ate_," ", Statistics.mean(tau))
hline!([ate_],label="ate")


text:
@testset "transform" begin

    @testset "basic" begin
        @test transform(2, 3) == Fronts.ϕ(2,3)
    end

    @testset "ODE types" begin
        eq = DiffusionEquation(identity)

        odefun = @inferred transform(eq)

        @test isa(odefun, ODEFunction)

        v = @inferred odefun((@SVector [1.0, 0.0]), NullParameters(), 0.0)
        @test isa(v, SVector)
        @test v == [0.0, 0.0]

        prob = CauchyProblem(eq, b=1, d_dϕb=0)

        @test isa(transform(prob), ODEProblem)
    end

end


text:
module PWscf

using REPL.Terminals: TTYTerminal
using REPL.TerminalMenus: RadioMenu, request

using Crayons.Box: RED_FG, GREEN_FG
using QuantumESPRESSOBase.Inputs: entryname, qestring
using QuantumESPRESSOBase.Inputs.PWscf:
    ControlNamelist,
    SystemNamelist,
    ElectronsNamelist,
    IonsNamelist,
    CellNamelist,
    AtomicSpecies,
    AtomicSpeciesCard,
    AtomicPosition,
    AtomicPositionsCard,
    KPointsCard,
    CellParametersCard,
    GammaPoint, MonkhorstPackGrid, KPointsCard,
    PWInput

using ..Inputs: namelist_builder, card_builder, input_builder, setfield_helper

import QuantumESPRESSOHelpers.Inputs

function Inputs.namelist_builder(
    terminal::TTYTerminal,
    ::Type{T},
) where {T<:ControlNamelist}
    calculations = pairs(("scf", "nscf", "bands", "relax", "md", "vc-relax", "vc-md"))
    restart_modes = pairs(("from_scratch", "restart"))
    calculation = calculations[request(
        terminal,
        GREEN_FG("What exact calculation do you want to run?") |> string,
        RadioMenu(collect(values(calculations))),
    )]
    restart_mode = restart_modes[request(
        terminal,
        GREEN_FG("Starting from scratch?") |> string,
        RadioMenu(["yes", "no"]),
    )]
    print(terminal, GREEN_FG("Convergence threshold on total energy (a.u): ") |> string)
    etot_conv_thr = parse(Float64, readline(terminal))
    print(terminal, GREEN_FG("Convergence threshold on forces (a.u): ") |> string)
    forc_conv_thr = parse(Float64, readline(terminal))
    control = T(
        calculation = calculation,
        restart_mode = restart_mode,
        etot_conv_thr = etot_conv_thr,
        forc_conv_thr = forc_conv_thr,
    )
    return setfield_helper(terminal, control)
end # function namelist_builder
function Inputs.namelist_builder(
    terminal::TTYTerminal,
    ::Type{T},
) where {T<:SystemNamelist}
    print(terminal, GREEN_FG("Please input the Bravais lattice index `ibrav`: ") |> string)
    ibrav = parse(Int, readline(terminal))
    print(
        terminal,
        GREEN_FG("Please input a `celldm` 1-6 (separated by spaces): ") |> string,
    )
    celldm = map(x -> parse(Float64, x), split(readline(terminal), " ", keepempty = false))
    print(
        terminal,
        GREEN_FG("Please input the number of atoms in the unit cell `nat`: ") |> string,
    )
    nat = parse(Int, readline(terminal))
    print(
        terminal,
        GREEN_FG("Please input the number of types of atoms in the unit cell `ntyp`: ") |>
        string,
    )
    ntyp = parse(Int, readline(terminal))
    print(
        terminal,
        GREEN_FG("Please input the kinetic energy cutoff (Ry) for wavefunctions `ecutwfc`: ") |>
        string,
    )
    ecutwfc = parse(Float64, readline(terminal))
    print(
        terminal,
        GREEN_FG("Please input the Kinetic energy cutoff (Ry) for charge density `ecutrho`: ") |>
        string,
    )
    ecutrho = parse(Float64, readline(terminal))
    system = T(
        ibrav = ibrav,
        celldm = celldm,
        nat = nat,
        ntyp = ntyp,
        ecutwfc = ecutwfc,
        ecutrho = ecutrho,
    )
    return setfield_helper(terminal, system)
end # function namelist_builder
function Inputs.namelist_builder(
    terminal::TTYTerminal,
    ::Type{T},
) where {T<:ElectronsNamelist}
    print(
        terminal,
        GREEN_FG("Please input the convergence threshold for selfconsistency `conv_thr`: ") |>
        string,
    )
    conv_thr = parse(Float64, readline(terminal))
    diagonalizations = pairs(("david", "cg", "cg-serial", "david-serial"))
    diagonalization = diagonalizations[request(
        terminal,
        GREEN_FG("Please input the diagonalization method `diagonalization`: ") |> string,
        RadioMenu(collect(values(diagonalizations))),
    )]
    electrons = T(conv_thr = conv_thr, diagonalization = diagonalization)
    return setfield_helper(terminal, electrons)
end # function namelist_builder
function Inputs.namelist_builder(
    terminal::TTYTerminal,
    ::Type{T},
) where {T<:IonsNamelist}
    ion_dynamics_pool =
        pairs(("none", "bfgs", "damp", "verlet", "langevin", "langevin-smc", "beeman"))
    ion_dynamics = ion_dynamics_pool[request(
        terminal,
        GREEN_FG("Please input the type of ionic dynamics `ion_dynamics`: ") |> string,
        RadioMenu(collect(values(ion_dynamics_pool))),
    )]
    ion_temperature_pool = (
        "rescaling",
        "rescale-v",
        "rescale-T",
        "reduce-T",
        "berendsen",
        "andersen",
        "initial",
        "not_controlled",
    )
    ion_temperature = ion_temperature_pool[request(
        terminal,
        GREEN_FG("Please input the ions temperature `ion_temperature`: ") |> string,
        RadioMenu(collect(values(ion_temperature_pool))),
    )]
    ions = T(ion_dynamics = ion_dynamics, ion_temperature = ion_temperature)
    return setfield_helper(terminal, ions)
end # function namelist_builder
function Inputs.namelist_builder(
    terminal::TTYTerminal,
    ::Type{T},
) where {T<:CellNamelist}
    cell_dynamics_pool = pairs(("none", "sd", "damp-pr", "damp-w", "bfgs", "pr", "w"))
    cell_dynamics = cell_dynamics_pool[request(
        terminal,
        GREEN_FG("Please input the type of dynamics for the cell `cell_dynamics`: ") |>
        string,
        RadioMenu(collect(values(cell_dynamics_pool))),
    )]
    print(
        terminal,
        GREEN_FG("Please input the target pressure [KBar] in a variable-cell md or relaxation run `press`: ") |>
        string,
    )
    press = parse(Float64, readline(terminal))
    print(
        terminal,
        GREEN_FG("Please input the fictitious cell mass [amu] for variable-cell simulations `wmass`: ") |>
        string,
    )
    wmass = parse(Float64, readline(terminal))
    print(
        terminal,
        GREEN_FG("Please input the Convergence threshold on the pressure for variable cell `press_conv_thr`: ") |>
        string,
    )
    press_conv_thr = parse(Float64, readline(terminal))
    cell = T(
        cell_dynamics = cell_dynamics,
        press = press,
        wmass = wmass,
        press_conv_thr = press_conv_thr,
    )
    return setfield_helper(terminal, cell)
end # function namelist_builder

function Inputs.card_builder(terminal::TTYTerminal, ::Type{T}) where {T<:KPointsCard}
    kpt_style = request(
        terminal,
        GREEN_FG("What k-point style do you want?") |> string,
        RadioMenu(["gamma", "automatic"]),
    )
    return if kpt_style == 1
        KPointsCard("gamma", GammaPoint())
    else  # "automatic"
        print(
            terminal,
            GREEN_FG("What 3-element k-point grid do you want (separated by spaces): ") |>
            string,
        )
        grid = map(x -> parse(Int, x), split(readline(terminal), " ", keepempty = false))
        print(
            terminal,
            GREEN_FG("What 3-element k-point offsets do you want (separated by spaces): ") |>
            string,
        )
        offsets = map(x -> parse(Int, x), split(readline(terminal), " ", keepempty = false))
        return KPointsCard("automatic", MonkhorstPackGrid(grid, offsets))
    end
end # function card_builder

function Inputs.input_builder(terminal::TTYTerminal, ::Type{T}) where {T<:PWInput}
    fields = Dict{Symbol,Any}()
    for S in (ControlNamelist, SystemNamelist, ElectronsNamelist)
        haserror = true
        while haserror
            try
                push!(fields, entryname(S) => namelist_builder(terminal, S))
                haserror = false
            catch e
                isa(e, InterruptException) && rethrow(e)
                println(terminal, RED_FG("Something wrong happens, try again!") |> string)
            end
        end
    end
    if fields[:control].calculation ∈ ("relax", "md", "vc-relax", "vc-md")
        haserror = true
        while haserror
            try
                push!(
                    fields,
                    entryname(IonsNamelist) => namelist_builder(terminal, IonsNamelist),
                )
                haserror = false
            catch e
                isa(e, InterruptException) && rethrow(e)
                println(terminal, RED_FG("Something wrong happens, try again!") |> string)
            end
        end
    else
        IonsNamelist()
    end
    if fields[:control].calculation ∈ ("vc-relax", "vc-md")
        haserror = true
        while haserror
            try
                push!(
                    fields,
                    entryname(CellNamelist) => namelist_builder(terminal, CellNamelist),
                )
                haserror = false
            catch e
                isa(e, InterruptException) && rethrow(e)
                println(terminal, RED_FG("Something wrong happens, try again!") |> string)
            end
        end
    else
        CellNamelist()
    end
    haserror = true
    while haserror
        try
            push!(fields, entryname(KPointsCard) => card_builder(terminal, KPointsCard))
            haserror = false
        catch e
            isa(e, InterruptException) && rethrow(e)
            println(terminal, RED_FG("Something wrong happens, try again!") |> string)
        end
    end
    push!(fields, entryname(AtomicSpeciesCard) => AtomicSpeciesCard(AtomicSpecies[]))
    push!(
        fields,
        entryname(AtomicPositionsCard) => AtomicPositionsCard("alat", AtomicPosition[]),
    )
    result = T(; fields...)
    saveresult = pairs((true, false))[request(
        terminal,
        GREEN_FG("Do you want to save the generated input to file?") |> string,
        RadioMenu(["yes", "no"]),
    )]
    if saveresult
        print(terminal, GREEN_FG("Input file name: ") |> string)
        write(readline(terminal) |> chomp, qestring(result))
    end
    return result
end # function input_builder

end # module PWscf


text:
# Assign the value 10 to the variable `x`
x = 10

# Doing math with the value of `x`
x + 1

# Reassign the value of `x`
x = 1 + 1

# Assign a value of a new **type** to `x`
x = "Hello, Julia World!"

# Flexible variable naming
x = 1.0

y = -3

Z = "My string"

customary_phrase = "Hello, Julia World!"

DeclarationOfIndependenceSelfEvidentTruths = "Wir halten diese Wahrheiten für selbstverständlich, dass alle Menschen gleich geschaffen sind, dass sie von ihrem Schöpfer mit bestimmen unveräußerlichen Rechten ausgestattet sind, darunter Leben, Freiheit and das Streben nach Glück."

# Unicode names (in UTF-8 encoding) are allowed
δ = 0.00001

selbstverständlich = "self-evident"

# Although one can redefine built-in constants and functions if necessary, it is **not** recommended to avoid
# potential confusion
#
# pi = 3
#
# pi  # prints "3"
#
# sqrt = 4  # changes the `sqrt` function

# If a built-in constant or function is already in use, Julia will report an error if you try to redefine it.
pi

pi = 3

sqrt(100)

sqrt = 4

# Allowed variable names
#
# Variables start with many Unicode code points in specific categories:
#
# - Lu (uppercase letter)
# - Ll (lowercase letter)
# - Lt
# - Lm
# - Lo
# - Nl (letters (numbers))
#
# - Sc / So currency and other symbols
# - A few other letter-like characters; for example, a subset of the Sm math symbols
#
# Subsequent characters include ! and digits (0-9 and other characters in categories Nd/No), as well as other
# Unicode code points:
# - Diacritics and other modifying marks (categories Mn/Mc/Me/Sk)
# - Some punctuation connectors (category Pc)
# - Primes
# - A few other characters

# Operators like `+` are also valid identifiers, but are parsed specially.

# The only explicitly disallowed names for variables are the names of the built-in `keywords`:

# else = false
# try = "No"

# Stylistic conventions for variable names
#
# - Names of variables are in lower case
# - Word separation can be indicated by underscore ('_'), but the use of underscores is **discouraged** unless
#   the name would be hard to read otherwise
# - Names of `Type`s and `Module`s begin with a capital letter and word separation is shown with upper camel case
#   instead of underscores (PascalCase)
# - Names of `functions` and `macros` are in lower case **without** underscores
# - Functions that mutate their arguments have names that end in '!'.
#
# For more information, see the Style Guide.


text:
@testset "SimpleEdgeIter" begin
    ga = @inferred(SimpleGraph(10, 20; seed=1))
    gb = @inferred(SimpleGraph(10, 20; seed=1))
    @test sprint(show, edges(ga)) == "SimpleEdgeIter 20"
    # note: we don't get the first iterator state,
    #since iterate returns the state after taking the first value
    @test sprint(show, iterate(edges(ga))[2]) == "SimpleEdgeIterState [1, 2]"

    @test length(collect(edges(Graph(0, 0)))) == 0

    @test @inferred(edges(ga)) == edges(gb)
    @test @inferred(edges(ga)) == collect(Edge, edges(gb))
    @test edges(ga) != collect(Edge, Base.Iterators.take(edges(gb), 5))

    @test collect(Edge, edges(gb)) == edges(ga)
    @test Set{Edge}(collect(Edge, edges(gb))) == edges(ga)
    @test @inferred(edges(ga)) == Set{Edge}(collect(Edge, edges(gb)))

    @test @inferred(eltype(edges(ga))) == eltype(typeof(edges(ga))) == edgetype(ga)
    @test eltype(collect(edges(ga))) == edgetype(ga)

    ga = @inferred(SimpleDiGraph(10, 20; seed=1))
    gb = @inferred(SimpleDiGraph(10, 20; seed=1))
    @test @inferred(edges(ga)) == edges(gb)
    @test @inferred(edges(ga)) == collect(SimpleEdge, edges(gb))
    @test collect(SimpleEdge, edges(gb)) == edges(ga)
    @test Set{Edge}(collect(SimpleEdge, edges(gb))) == edges(ga)
    @test @inferred(edges(ga)) == Set{SimpleEdge}(collect(SimpleEdge, edges(gb)))
    @test eltype(collect(edges(ga))) == edgetype(ga)

    ga = SimpleGraph(10)
    add_edge!(ga, 3, 2)
    add_edge!(ga, 3, 10)
    add_edge!(ga, 5, 10)
    add_edge!(ga, 10, 3)

    e1 = Edge(3, 10)
    e2 = (3, 10)
    @test e1 ∈ edges(ga)
    @test e2 ∈ edges(ga)
    @test (3, 9) ∉ edges(ga)

    for u in 1:20, v in 1:20
      b = has_edge(ga, u, v)
      @test b == @inferred (u, v) ∈ edges(ga)
      @test b == @inferred (u => v) ∈ edges(ga)
      @test b == @inferred Edge(u, v) ∈ edges(ga)
    end

    eit = edges(ga)
    # @inferred not valid for new interface anymore (return type is a Union)
    es = iterate(eit)[2]

    @test es.s == 3
    @test es.di == 2

    @test [e for e in eit] == [Edge(2, 3), Edge(3, 10), Edge(5, 10)]

    gb = copy(ga)
    add_vertex!(gb)
    @test edges(ga) == edges(gb)
    @test edges(gb) == edges(ga)

    ga = SimpleDiGraph(10)
    add_edge!(ga, 3, 2)
    add_edge!(ga, 3, 10)
    add_edge!(ga, 5, 10)
    add_edge!(ga, 10, 3)

    for u in 1:20, v in 1:20
      b = has_edge(ga, u, v)
      @test b == @inferred (u, v) ∈ edges(ga)
      @test b == @inferred (u => v) ∈ edges(ga)
      @test b == @inferred Edge(u, v) ∈ edges(ga)
    end

    eit = @inferred(edges(ga))
    es = iterate(eit)[2]

    @test es.s == 3
    @test es.di == 2

    @test [e for e in eit] == [
      SimpleEdge(3, 2), SimpleEdge(3, 10),
      SimpleEdge(5, 10), SimpleEdge(10, 3)
    ]

    gb = copy(ga)
    add_vertex!(gb)
    @test edges(ga) == edges(gb)
    @test edges(gb) == edges(ga)
end


text:
using ZernikePolynomials
using Test

@testset "ZernikePolynomials.jl" begin

    @testset "Sequential index conversion" begin
        @test mn2Noll(0,0) == 1
        @test mn2Noll(-1,3) == 7
        @test isnan(mn2Noll(-1,0))
        @test mn2Noll(3,5) == 18

        @test mn2OSA(0,0) == 0
        @test mn2OSA(-1,3) == 7
        @test isnan(mn2OSA(-1,0))
        @test mn2OSA(3,5) == 19

        @test all([mn2OSA(OSA2mn(i)...) for i in 0:30] .== collect(0:30))
        @test all([mn2Noll(Noll2mn(i)...) for i in 1:31] .== collect(1:31))
    end

    @testset "Generation of Zernike polynomials" begin

        Z = Zernike(1,1,coord=:cartesian)
        @test Z isa Function

        Z = Zernike(1,1,coord=:polar)
        @test Z isa Function

        x = LinRange(-1,1,64)
        ϕi(i) = evaluateZernike(x,i,1.)
        r = 1:10
        N = sum(evaluateZernike(x,0,1.))
        X = round.([sum(ϕi(i) .* ϕi(j))/N for i in r, j in r])
        for i in 1:10, j in 1:10
            if i == j
                @test X[i,j] ≈ 1.0
            else
                @test X[i,j] == 0.0
            end
        end

        # coefficient estimation
        ϕ = evaluateZernike(64, [2, 3], [0.5, 0.3], index=:OSA)
        @test Zernikecoefficients(ϕ, [1, 2, 3]) ≈ [0.0, 0.5,0.3]
    end
end


text:
module Outliernicers

using Dates
using DataFrames
using Statistics
using StatsBase: iqr, quantile, sample
using TSML.ValDateFilters

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils
import AutoMLPipeline.AbsTypes: fit!, transform!

export fit!,transform!
export Outliernicer

"""
    Outliernicer(Dict(
       :dateinterval => Dates.Hour(1),
       :nnsize => 1,
       :missdirection => :symmetric
    ))

Detects outliers below or above (q25-iqr,q75+iqr)
and calls DateValNNer to replace them with nearest neighbors.

Example:

    fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
    csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
    valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
    valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
    stfier = Statifier(Dict(:processmissing=>true))
    mono = Monotonicer(Dict())
    outliernicer = Outliernicer(Dict(:dateinterval=>Dates.Hour(1)))

    mpipeline = @pipeline csvfilter |> valgator |> mono |> valnner |> outliernicer |> stfier
    results = fit_transform!(mpipeline)

Implements: `fit!`, `transform!`
"""
mutable struct Outliernicer <: Transformer
  model
  args

  function Outliernicer(args=Dict())
    default_args = Dict(
        :dateinterval => Dates.Hour(1),
        :nnsize => 1,
        :missdirection => :symmetric
    )
    new(nothing, mergedict(default_args, args))
  end
end

"""
    fit!(st::Outliernicer, features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}

Check that `features` are two-colum data.
"""
function fit!(st::Outliernicer, features::DataFrame, labels::Vector=[])
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  st.model = st.args
end

"""
    transform!(st::Outliernicer, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Locate outliers based on IQR factor and calls DateValNNer to replace them with nearest neighbors.
"""
function transform!(st::Outliernicer, features::DataFrame)
  features != DataFrame() || return DataFrame()
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  sum(names(features) .== (:Date,:Value))  == 2 || error("wrong column names")
  mfeatures=deepcopy(features)
  rvals = mfeatures.Value
  # compute the outlier range
  # setup to store both missing and numbers
  mvals = Array{Union{Missing,eltype(rvals)},1}(missing,length(rvals))
  mvals .= rvals
  crvals = skipmissing(rvals) # stat of non-missing
  miqr = iqr(crvals)
  q25,q75 = quantile(crvals,[0.25,0.75])
  lower=q25-miqr; upper=q75+miqr
  missindx = findall(x -> !ismissing(x) && (x > upper || x < lower),rvals)
  mvals[missindx] .= missing
  mfeatures.Value = mvals
  # use ValNNer to replace missings
  valnner = DateValNNer(st.args)
  fit!(valnner,mfeatures)
  resdf = transform!(valnner,mfeatures)
  resdf.Value = collect(skipmissing(resdf.Value))
  resdf
end

end


text:
import Term: cleantext, chars, textlen, split_lines

function same_widths(text::String)::Bool
    widths = textlen.(split_lines(text))
    return length(unique(widths)) == 1
end

function check_widths(text, width)
    for line in split_lines(text)
        @test textlen(line) <= width
    end
end

"""
Extensively test a panel making sure it has
the right size and Measure
"""
function testpanel(p, w, h)
    # check all lines have the same length
    _p = string(p)

    dw = displaysize(stdout)[2]
    if isnothing(w) || w > dw
        return
    else
        widths = textwidth.(cleantext.(split(_p, "\n")))
    end

    # println(p, p.measure, widths)
    @test length(unique(widths)) == 1

    # check it has the right measure
    if !isnothing(w)
        @test p.measure.w == w
        @test textlen(cleantext(p.segments[1].text)) == w
        @test length(chars(cleantext(p.segments[1].text))) == w
    end

    if !isnothing(h)
        @test p.measure.h == h
        @test length(p.segments) == h
    end
end

nlines(x) = length(split(x, "\n"))
lw(x) = max(length.(split(x, "\n"))...)

"""
Include but with a timer
"""

macro timeit_include(path::AbstractString)
    return :(@timeit TIMEROUTPUT $path include($path))
end

# ------------------------------- test no throw ------------------------------ #

struct NoException <: Exception end

macro test_nothrow(ex)
    esc(:(@test_throws NoException ($(ex); throw(NoException()))))
end