text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GTK4"
version = v"4.6.0"

# Collection of sources required to build GTK
sources = [
    # https://download.gnome.org/sources/gtk/
    ArchiveSource("https://download.gnome.org/sources/gtk/4.6/gtk-4.6.0.tar.xz",
                  "782d5951fbfd585fc9ec76c09d07e28e6014c72db001fb567fff217fb96e4d8c"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gtk*/

# We need to run some commands with a native Glib
apk add glib-dev

# This is awful, I know
ln -sf /usr/bin/glib-compile-resources ${bindir}/glib-compile-resources
ln -sf /usr/bin/glib-compile-schemas ${bindir}/glib-compile-schemas
ln -sf /usr/bin/gdk-pixbuf-pixdata ${bindir}/gdk-pixbuf-pixdata
# Remove gio-2.0 pkgconfig file so that it isn't picked up by post-install script.
rm ${prefix}/lib/pkgconfig/gio-2.0.pc

# llvm-ar seems to generate corrupted static archives:
#
#     [974/980] Linking target gtk/libgtk-4.1.dylib
#     ninja: job failed: [...]
#     ld: warning: ignoring file gtk/css/libgtk_css.a, building for macOS-x86_64 but attempting to link with file built for unknown-unsupported file format ( 0x21 0x3C 0x74 0x68 0x69 0x6E 0x3E 0x0A 0x2F 0x20 0x20 0x20 0x20 0x20 0x20 0x20 )
#     ld: warning: ignoring file gtk/libgtk.a, building for macOS-x86_64 but attempting to link with file built for unknown-unsupported file format ( 0x21 0x3C 0x74 0x68 0x69 0x6E 0x3E 0x0A 0x2F 0x20 0x20 0x20 0x20 0x20 0x20 0x20 )
#     ld: warning: ignoring file gsk/libgsk.a, building for macOS-x86_64 but attempting to link with file built for unknown-unsupported file format ( 0x21 0x3C 0x74 0x68 0x69 0x6E 0x3E 0x0A 0x2F 0x20 0x20 0x20 0x20 0x20 0x20 0x20 )
#     ld: warning: ignoring file gsk/libgsk_f16c.a, building for macOS-x86_64 but attempting to link with file built for unknown-unsupported file format ( 0x21 0x3C 0x74 0x68 0x69 0x6E 0x3E 0x0A 0x2F 0x20 0x20 0x20 0x20 0x20 0x20 0x20 )
#     Undefined symbols for architecture x86_64:
#       "_gtk_make_symbolic_pixbuf_from_data", referenced from:
#           _main in encodesymbolic.c.o
#     ld: symbol(s) not found for architecture x86_64
if [[ "${target}" == *apple* ]]; then
    sed -i "s?^ar = .*?ar = '/opt/${target}/bin/${target}-ar'?g" "${MESON_TARGET_TOOLCHAIN}"
fi

# https://gitlab.gnome.org/GNOME/gtk/-/merge_requests/4443
atomic_patch -p1 ../patches/0001-Include-gdk-private.h-to-fix-error-about-g_source_se.patch
# https://gitlab.gnome.org/GNOME/gtk/-/merge_requests/4445
atomic_patch -p1 ../patches/0001-gdkjpeg-include-stdlib.h-necessary-for-free.patch

FLAGS=()
if [[ "${target}" == *-apple-* ]]; then
    FLAGS+=(-Dx11-backend=false -Dwayland-backend=false)
elif [[ "${target}" == *-freebsd* ]]; then
    FLAGS+=(-Dwayland-backend=false)
elif [[ "${target}" == *-mingw* ]]; then
    # Need to tell we're targeting at least Windows 7 so that `GC_ALLGESTURES` is defined
    sed -ri "s/^c_args = \[(.*)\]/c_args = [\1, '-DWINVER=_WIN32_WINNT_WIN7']/" ${MESON_TARGET_TOOLCHAIN}
fi

mkdir build-gtk && cd build-gtk
meson .. \
    --buildtype=release \
    -Dmedia-gstreamer=disabled \
    -Dintrospection=disabled \
    -Ddemos=false \
    -Dbuild-examples=false \
    -Dbuild-tests=false \
    -Dgtk_doc=false \
    "${FLAGS[@]}" \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install

# post-install script is disabled when cross-compiling
glib-compile-schemas ${prefix}/share/glib-2.0/schemas

# Remove temporary links
rm ${bindir}/gdk-pixbuf-pixdata ${bindir}/glib-compile-{resources,schemas}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libgtk-4", :libgtk4),
]

x11_platforms = filter(p -> Sys.islinux(p) || Sys.isfreebsd(p), platforms)

# Dependencies that must be installed before this package can be built
dependencies = [
    # Need a native `sassc`
    HostBuildDependency("SassC_jll"),
    # Need a host Wayland for wayland-scanner
    HostBuildDependency("Wayland_jll"; platforms=x11_platforms),
    BuildDependency("Xorg_xorgproto_jll"; platforms=x11_platforms),
    Dependency("Glib_jll"; compat="2.68.3"),
    Dependency("Graphene_jll"; compat="1.10.6"),
    Dependency("Cairo_jll"),
    Dependency("Pango_jll"; compat="1.50.3"),
    Dependency("FriBidi_jll"),
    Dependency("FreeType2_jll"),
    Dependency("gdk_pixbuf_jll"),
    Dependency("Libepoxy_jll"),
    Dependency("HarfBuzz_jll"),
    Dependency("xkbcommon_jll"; platforms=x11_platforms),
    Dependency("iso_codes_jll"),
    Dependency("Wayland_jll"; platforms=x11_platforms),
    Dependency("Wayland_protocols_jll"; compat="1.23", platforms=x11_platforms),
    Dependency("Xorg_libXrandr_jll"; platforms=x11_platforms),
    Dependency("Xorg_libX11_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXrender_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXi_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXext_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXcursor_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXdamage_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXfixes_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXinerama_jll"; platforms=x11_platforms),
    Dependency("Fontconfig_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
    using PersistenceBarcodes
using Base.Test
using Plots; plotly()
using Suppressor

@testset "PersistenceBarcodes" begin
    @testset "PersistencePair" begin
        # Types
        pair1 = PersistencePair(0, 2)
        @test pair1 isa PersistencePair{Int, Void}
        pair2 = PersistencePair(0, Inf, [1])
        @test pair2 isa PersistencePair{Float64, Vector{Int}}

        # Getters
        @test birth(pair1) ≡ 0
        @test birth(pair2) ≡ 0.0
        @test death(pair1) ≡ 2
        @test death(pair2) ≡ Inf
        @test data(pair1)  ≡ nothing
        @test data(pair2) == [1]

        # Deconstruction
        b1, d1 = pair1
        @test b1 ≡ 0 && d1 ≡ 2
        b2, d2 = pair2
        @test b2 ≡ 0.0 && d2 ≡ Inf

        # Show
        @test @capture_out(print(pair1)) ==
            "PersistencePair{$Int, Void}(0, 2)"
        @test @capture_out(print(pair2)) ==
            "PersistencePair{Float64, Array{$Int,1}}(0.0, Inf; [1])"
        @test @capture_out(showcompact(pair1)) == "(0, 2)"
        @test @capture_out(showcompact(pair2)) == "(0.0, Inf; [1])"

        # Make sure the correct IO is used.
        @test @capture_err(print(STDERR, pair1)) ==
            "PersistencePair{$Int, Void}(0, 2)"
        @test @capture_err(print(STDERR, pair2)) ==
            "PersistencePair{Float64, Array{$Int,1}}(0.0, Inf; [1])"
        @test @capture_err(showcompact(STDERR, pair1)) == "(0, 2)"
        @test @capture_err(showcompact(STDERR, pair2)) == "(0.0, Inf; [1])"
    end

    @testset "PersistenceBarcode" begin
        # Construction
        barcode1 = PersistenceBarcode(PersistencePair.(Float64[0,0,0], [1,2,3], [1,2,3]))
        barcode2 = PersistenceBarcode([PersistencePair.(Float64[0,0,0], [1,2,3]),
                                       PersistencePair.(Float64[1,0],   [Inf,2])])
        barcode3 = PersistenceBarcode(PersistencePair.(Float64[0,0,0], [1,2,3]),
                                      PersistencePair.(Float64[1,0],   [Inf,2]))
        @test barcode2 == barcode3
        @test barcode1 isa PersistenceBarcode{Float64, Int}
        @test barcode2 isa PersistenceBarcode{Float64, Void}

        @test dim(barcode1) == 0
        @test dim(barcode2) == 1

        # Indexing
        @test barcode1[0] == PersistencePair.([0.0,0.0,0.0], [1.0,2.0,3.0], [1,2,3])
        @test barcode2[1] == PersistencePair.([1.0,0.0], [Inf,2.0])

        # Show
        bc1_str =
            """
        0-d PersistenceBarcode{Float64, $Int}:
         dim 0:
          (0.0, 1.0; 1)
          (0.0, 2.0; 2)
          (0.0, 3.0; 3)
        """
        @test @capture_out(println(barcode1)) == bc1_str
        @test @capture_err(println(STDERR, barcode1)) == bc1_str
        @test @capture_out(showcompact(barcode1)) ==
            "0-d PersistenceBarcode{Float64, $Int}"
        @test @capture_err(showcompact(STDERR, barcode1)) ==
            "0-d PersistenceBarcode{Float64, $Int}"

        bc2_str =
            """
        1-d PersistenceBarcode{Float64, Void}:
         dim 0:
          (0.0, 1.0)
          (0.0, 2.0)
          (0.0, 3.0)
         dim 1:
          (1.0, Inf)
          (0.0, 2.0)
        """
        @test @capture_out(println(barcode2)) == bc2_str
        @test @capture_err(println(STDERR, barcode2)) == bc2_str
        @test @capture_out(showcompact(barcode2)) ==
            "1-d PersistenceBarcode{Float64, Void}"
        @test @capture_err(showcompact(STDERR, barcode2)) ==
            "1-d PersistenceBarcode{Float64, Void}"

        # Iteration
        l = Int[]
        T = eltype(barcode2)
        for p in barcode2
            push!(l, length(p))
            @test p isa T
        end
        @test l == [3, 2]

        l = map(length, collect(barcode2))
        @test l == [3, 2]

        @test filter(p -> !isfinite(death(p)), barcode2) ==
            PersistenceBarcode(PersistencePair{Float64, Void}[],
                               [PersistencePair(1.0, Inf)])

        @test map(p -> PersistencePair(death(p),0), barcode2) ==
            PersistenceBarcode(PersistencePair.([1.0,2.0,3.0],[0,0,0]),
                               PersistencePair.([Inf,2.0], [0,0]))

        # map, filter
    end

    @testset "plot" begin
        barcode = PersistenceBarcode(PersistencePair.(Float64[0,0,0], [1,2,3]),
                                     PersistencePair.(Float64[1,0],   [Inf,2]))

        @test PersistenceBarcodes.getlastdeath(barcode, 0:1) ≡ 3.0
        @test PersistenceBarcodes.getlastdeath(barcode, 0)   ≡ 3.0
        @test PersistenceBarcodes.getlastdeath(barcode, 1)   ≡ 2.0

        # Barcode plot
        @test plot(barcode)                ≠ nothing
        @test plot(barcode, dims = 0:1)    ≠ nothing
        @test plot(barcode, dims = 1)      ≠ nothing
        @test plot(barcode, dims = 0)      ≠ nothing
        @test plot(barcode, dims = [1, 0]) ≠ nothing

        @test_throws ArgumentError plot(barcode, dims = [1, 0, -1])
        @test_throws ArgumentError plot(barcode, dims = [0, 2])
        @test_throws ArgumentError plot(barcode, dims = 2)
        @test_throws ArgumentError plot(barcode, dims = -1)

        # Diagram plot
        @test persistencediagram(barcode)                ≠ nothing
        @test persistencediagram(barcode, dims = 0:1)    ≠ nothing
        @test persistencediagram(barcode, dims = 1)      ≠ nothing
        @test persistencediagram(barcode, dims = 0)      ≠ nothing
        @test persistencediagram(barcode, dims = [1, 0]) ≠ nothing

        @test_throws ArgumentError persistencediagram(barcode, dims = [1, 0, -1])
        @test_throws ArgumentError persistencediagram(barcode, dims = [0, 2])
        @test_throws ArgumentError persistencediagram(barcode, dims = 2)
        @test_throws ArgumentError persistencediagram(barcode, dims = -1)
        @test_throws ArgumentError persistencediagram(barcode, [1,2,3])
        @test_throws ArgumentError persistencediagram([1,2,3])

        # Empty plots
        T = PersistencePair{Float64, Void}
        @test plot(PersistenceBarcode(T[], T[]))               ≠ nothing
        @test persistencediagram(PersistenceBarcode(T[], T[])) ≠ nothing
    end
end

ext: jl


text:
using Test
using acopf
using TimerOutputs

@testset "Ipopt CPU evaluation" begin
  arraytype = Array
  include("Ipopt.jl")
  case="data/case9"
  max_iter=100
  opfdata = acopf.opf_loaddata(case)
  Pg0, Qg0, Vm0, Va0 = acopf.initialPt_IPOPT(opfdata)
  obj_val = test(Pg0, Qg0, Vm0, Va0, timeroutput, opfdata, arraytype; max_iter = max_iter)
  @test obj_val ≈ 5.2966862e+03
end

@testset "Ipopt GPU evaluation" begin
  arraytype = CuArray
  include("Ipopt.jl")
  case="data/case9"
  max_iter=100
  opfdata = acopf.opf_loaddata(case)
  Pg0, Qg0, Vm0, Va0 = acopf.initialPt_IPOPT(opfdata)
  obj_val = test(Pg0, Qg0, Vm0, Va0, timeroutput, opfdata, arraytype; max_iter = max_iter)
  @test obj_val ≈ 5.2966862e+03
end

@testset "JuMP ACOPF Reference" begin
  include("jumpmodel.jl")
  case="data/case9"
  max_iter=100
  opfdata = acopf.opf_loaddata(case)
  Pg0, Qg0, Vm0, Va0 = acopf.initialPt_IPOPT(opfdata)
  opfmodel, Pg, Qg, Va, Vm = model(opfdata; max_iter = max_iter)
  opfmodel,status = solve(opfmodel,opfdata)
  if status==MOI.LOCALLY_SOLVED
    outputAll(opfmodel,opfdata, Pg, Qg, Va, Vm)
  end
  @test objective_value(opfmodel) ≈ 5.2966862e+03
end

@testset "JuMP ACOPF with HiOp" begin
  include("jumpmodel.jl")
  case="data/case9"
  max_iter=100
  opfdata = acopf.opf_loaddata(case)
  Pg0, Qg0, Vm0, Va0 = acopf.initialPt_IPOPT(opfdata)
  opfmodel, Pg, Qg, Va, Vm = model(opfdata; solver="Hiop")
  opfmodel,status = solve(opfmodel,opfdata)
  if status==MOI.LOCALLY_SOLVED
    outputAll(opfmodel,opfdata, Pg, Qg, Va, Vm)
  end
  @test objective_value(opfmodel) ≈ 5.2966862e+03
end


text:
module MoonCreature

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Moon Creature" => Ahorn.EntityPlacement(
        Maple.MoonCreature
    )
)

sprite = "scenery/moon_creatures/tiny05"

function Ahorn.selection(entity::Maple.MoonCreature)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.MoonCreature, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end


text:
using Remark, FileWatching
using ArgParse

function _parse_args()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--watch"
        action = :store_true
        help = "watch files and rebuilt on changes"
    end
    parse_args(s)
end

function main(watch::Bool)
    Remark.slideshow(
        @__DIR__;
        options = Dict("ratio" => "16:9"),
        title = "Meta-Learning and Differentiable Convex Optimization with Flux and Zygote",
    )
    @info "slides build done"
    while watch
        FileWatching.watch_folder(joinpath(@__DIR__, "src"))
        @info "Rebuilt"
        Remark.slideshow(
            @__DIR__;
            options = Dict("ratio" => "16:9"),
            title = "Meta-Learning and Differentiable Convex Optimization with Flux and Zygote",
        )
    end
end

if nameof(@__MODULE__) == :Main
    args = _parse_args()
    main(args["watch"])
end


text:
__precompile__()
module TraceCalls

using Requires
using MacroTools
using MacroTools: combinedef, combinearg, longdef1
using Base.Test: @inferred
using DataStructures: OrderedDict, OrderedSet
using Memoize
using Base: url
using Base.Core: MethodInstance
using Crayons: Crayon, inv
import Base: +, -

export @traceable, @trace, Trace, prune, FontColor, Bold,
    is_inferred, map_is_inferred, redgreen, greenred,
    compare_past_trace, filter_func, apply_macro, @stacktrace, measure, tree_size,
    is_mutating, REPR, filter_cutting, NoTraceable, trace_log, filter_lineage,
    bottom, top, highlight, @show_val_only_type, objects_in, signature, group_by,
    map_groups, trace_benchmark, @compilation_times, specializations, is_root,
    function_is

const traceable_definitions = OrderedDict()  # We use an OrderedDict because in case we
                                             # accidentally store the same definition
                                             # twice, at least the latter one takes
                                             # precedence.

################################################################################
# Utils

macro ignore_errors(error_rval, expr)
    quote
        try
            $(esc(expr))
        catch e
            $(esc(error_rval))
        end
    end
end

macro return_error(expr)
    quote
        try
            $(esc(expr))
        catch e
            return e
        end
    end
end

return_error_fn(f) = x->@return_error(f(x))

################################################################################

""" A `Trace` object represents a function call. It has fields `func, args, kwargs,
called, value`, with `func(args...; kwargs...) = value`. `called::Vector{Trace}` of the
immediate traceable function calls that happened during the execution of this trace.

There are no accessors for `Trace`; please reference its fields directly.
For instance, `filter(tr->isa(tr.args[1], Int), trace)` will select all calls whose
first argument is an `Int`. """
struct Trace
    func             # the function/callable called
    args::Tuple      # the positional arguments
    kwargs::Tuple    # the keyword arguments
    called::Vector{Trace}
    value       # This is the return value of the func(args...) call, but it's also where
                # the result of `map(f, ::Trace)` will be stored.
end
# Convenience constructors
Trace(tr::Trace, new_value) = Trace(tr.func, tr.args, tr.kwargs, tr.called, new_value)
Trace(tr::Trace, called::Vector{Trace}) =
    Trace(tr.func, tr.args, tr.kwargs, called, tr.value)

function Base.copy!(dest::Trace, src::Trace)
    dest.func = src.func
    dest.args = src.args
    dest.kwargs = src.kwargs
    dest.called = src.called
    dest.value = src.value
end

tree_size(tr::Trace) = 1 + mapreduce(tree_size, +, 0, tr.called)
Base.copy(tr::Trace) = Trace(tr.func, tr.args, tr.kwargs, tr.called, tr.value)
Base.push!(tr::Trace, sub_trace::Trace) = push!(tr.called, sub_trace)
Base.getindex(tr::Trace, i::Int, args...) = tr.called[i][args...]
Base.getindex(tr::Trace) = tr
Base.length(tr::Trace) = length(tr.called)
Base.endof(tr::Trace) = length(tr)
Base.url(tr::Trace) = @ignore_errors "" Base.url(which(tr))
Base.which(tr::Trace) = apply_macro(:@which, tr)
Base.code_llvm(tr::Trace) = apply_macro(:@code_llvm, tr)
Base.code_lowered(tr::Trace) = apply_macro(:@code_lowered, tr)
Base.code_native(tr::Trace) = apply_macro(:@code_native, tr)
Base.code_typed(tr::Trace) = apply_macro(:@code_typed, tr)
Base.code_warntype(tr::Trace) = apply_macro(:@code_warntype, tr)
Base.less(tr::Trace) = apply_macro(:@less, tr)
Base.edit(tr::Trace) = apply_macro(:@edit, tr)
Base.sum(f::Function, tr::Trace) = sum(f, collect(tr))
Base.sum(tr::Trace) = sum(t->t.value, tr)
Base.all(f::Function, tr::Trace) = all(f, collect(tr))
Base.all(tr::Trace) = all(t->value(t), collect(tr))
Base.maximum(tr::Trace) = maximum(sub.value for sub in collect(tr))
Base.round(tr::Trace, n::Int) = map(sub->round(sub.value, n), tr)
Base.signif(tr::Trace, n::Int) = map(sub->signif(sub.value, n), tr)
Base.normalize(tr::Trace, div=value(tr)) = map(apply_to_value_fn(x->x/div), tr)
Base.sort(tr::Trace; rev=false, by=tr->value(tr)) =
    Trace(tr, sort(tr.called; rev=rev, by=by))

narrow_typeof{T}(t::Type{T}) = Type{T}
narrow_typeof{T}(t::T) = T
""" `signature(tr::Trace) = (tr.func, map(typeof, tr.args)...)` - the arguments that this
call dispatches on. """
signature(tr::Trace) = (tr.func, map(narrow_typeof, tr.args)...)
arg_names(method::Method) = [Symbol(first(a))
                             for a in Base.arg_decl_parts(method)[2][2:end]
                             if first(a)!=""]
arg_names(tr::Trace) = arg_names(which(tr))  # positional arguments only
function Base.get(tr::Trace, arg::Symbol, default)
    if (i = findfirst(a->a==arg, arg_names(tr))) > 0
        tr.args[i]
    elseif (k = findfirst(p->p[1]==arg, tr.kwargs)) > 0
        tr.kwargs[k][2]
    else
        default
    end
end
struct Dummy end
Base.getindex(tr::Trace, arg::Symbol) =
    (r=Base.get(tr, arg, Dummy())) == Dummy() ? throw(KeyError(arg)) : r
Base.keys(tr::Trace) = [arg_names(tr)...; map(first, tr.kwargs)...]
Base.haskey(tr::Trace, key::Symbol) = key in arg_names(tr) || any(key==first(kwa)
                                                                  for kwa in tr.kwargs)
""" `objects_in(tr::Trace)` returns the call's arguments and return value in a vector.
Useful for filtering, eg. `filter(tr->any(obj isa Number && obj < 0 for obj in objects_in(tr)), trace)` will keep all calls that contain some negative number. """
objects_in(tr::Trace) = [tr.args..., map(last, tr.kwargs)..., value(tr)]

function_is(funs::Function...) = tr->tr.func in funs

# I've disabled iteration because it doesn't align with our desired `Base.map`'s
# behaviour, and it's kinda useless anyway.
# Base.start(tr::Trace) = 1
# Base.next(tr::Trace, i::Int) = (tr[i], i+1)
# Base.done(tr::Trace, i::Int) = i == length(tr)+1
(tr::Trace)() = tr.func(tr.args...; tr.kwargs...)

empty_trace() = nothing
const empty_trace_dummy = Trace(empty_trace, (), (), [], nothing)

struct Bottom
    delta::Int   # a positive number
end
""" `some_trace[bottom-2]` returns the subtrace that is the caller of the caller of the
deepest call in `some_trace`. """
const bottom = Bottom(0)
deepest_call(tr) = tr.called[findmax(map(depth, tr.called))[2]]
-(bot::Bottom, d::Int) = Bottom(bot.delta + d)
Base.getindex(tr::Trace, bot::Bottom, args...) =
    (depth(tr) <= bot.delta + 1) ? tr[args...] : deepest_call(tr)[bot, args...]
depth(tr::Trace) = isempty(tr.called) ? 1 : 1 + maximum(depth, tr.called)

struct Top
    delta::Int   # a positive number
end
""" `some_trace[top+2]` is equivalent to `some_trace[i, j]` where `i::Int` and `j::Int`
are the indexes of the deepest call trees in `some_trace` and `some_trace[i]`,
respectively. """
const top = Top(0)
+(top::Top, d::Int) = Top(top.delta + d)
Base.getindex(tr::Trace, top::Top, args...) =
    top.delta==0 ? tr[args...] : deepest_call(tr)[Top(top.delta-1), args...]

apply_macro(mac::Symbol, tr::Trace, mod::Module=Main) =
    eval(mod, Expr(:macrocall, mac,
                   (isempty(tr.kwargs) ?
                    # Needs special casing because @inferred chokes on kwargs-less
                    # funcalls otherwise.
                    :($(tr.func)($(map(QuoteNode, tr.args)...))) :
                    :($(tr.func)($(map(QuoteNode, tr.args)...); $(tr.kwargs...))))))
apply_macro(mac::Expr, tr::Trace, mod::Module=Main) =
    (mac.head==:macrocall ? apply_macro(only(mac.args), tr, mod) :
     error("Unable to call macro $mac"))

""" `map(f, tr::Trace)` recursively applies the function f to each `Trace` in `tr`,
and stores the result in `Trace.value` """
Base.map(f::Function, tr::Trace) = Trace(tr.func, tr.args, tr.kwargs,
                                         [map(f, c) for c in tr.called], f(tr))
""" `map_trace(f, tr)` is like `map`, but we expect `f(trace)` to return a new Trace """
map_trace(f::Function, tr::Trace) = f(Trace(tr.func, tr.args, tr.kwargs,
                                            [map_trace(f, c) for c in tr.called],
                                            tr.value))::Trace
apply_macro_fn(mac::Union{Symbol, Expr}) = sub->apply_macro(mac, sub)
Base.map(mac::Union{Symbol, Expr}, tr::Trace) = map(apply_macro_fn(mac), tr)
""" Similar to `map`, but we apply `filter_fun(f(tr))` to each call, and if it's false,
we filter out that call and all its descendents. Efficiently implemented. """
function map_filter(f::Function, filter_fun::Function, tr::Trace)
    # This was written for measure, but we turned out not to use it. Delete?
    res = f(tr)
    if filter_fun(res)
        Trace(tr.func, tr.args, tr.kwargs,
              filter(x->x!==empty_trace_dummy,
                     [map_filter(f, filter_fun, c) for c in tr.called]), res)
    else
        empty_trace_dummy
    end
end
handle_mac(f::Function) = f
handle_mac(mac::Union{Symbol, Expr}) = apply_macro_fn(mac)

################################################################################

struct REPR
    text
    html
    function REPR(x)
        s_text =
        s_html = IOBuffer(); try_show_val(s_html, MIME"text/html"(), x)
        new(get_io_output() do io; try_show_val(io, MIME"text/plain"(), x) end,
            get_io_output() do io; try_show_val(io, MIME"text/html"(), x) end)
    end
end
show_val(io::IO, ::MIME"text/plain", r::REPR) = write(io, r.text)
show_val(io::IO, ::MIME"text/html", r::REPR) = write(io, r.html)

################################################################################

""" `filter_cutting(f::Function, tr::Trace)` removes all subtraces (and their
callees) for which `f(tr)` is false. """
filter_cutting(f::Function, tr::Trace) =
    Trace(tr, Trace[filter_cutting(f, sub_tr) for sub_tr in tr.called if f(sub_tr)])

"""    filter_lineage(f::Function, tr::Trace; highlight=true, keep_descendents=true)

keeps all subtraces for which `f(::Trace)` is true of some of its descendents OR
ancestors. """
function filter_lineage(f::Function, tr::Trace; highlight=true, keep_descendents=true)
    if f(tr)
        # FIXME: we should probably go down and keep those descendents for which f(tr) is
        # true.
        res = keep_descendents ? tr : prune(tr)
    else
        called0 = Trace[filter_lineage(f, sub_tr; highlight=false,
                                       keep_descendents=keep_descendents)
                        for sub_tr in tr.called]
        called = filter(c->c!=empty_trace_dummy, called0)
        res = isempty(called) ? empty_trace_dummy : Trace(tr, called)
    end
    return highlight ? TraceCalls.highlight(f, res) : res
end

filter_descendents(f, tr) = # helper
    # Special casing because of #18852
    isempty(tr.called) ? Trace[] : Trace[t for sub in tr.called for t in filter_(f, sub)]
filter_(f, tr) =
    f(tr) ? [Trace(tr, filter_descendents(f, tr))] : filter_descendents(f, tr)
Base.filter(f::Function, tr::Trace) = Trace(tr, filter_descendents(f, tr))

""" `filter_func(functions::Vector, tr::Trace)` keeps only Trace objects whose function
is one of `functions` """
filter_func(functions::Vector, tr::Trace) = filter(tr->tr.func in functions, tr)
filter_func(func::Function, tr::Trace) = filter_func([func], tr)

function collect_!(trace_list, tr::Trace)
    push!(trace_list, tr)
    for sub in tr.called
        collect_!(trace_list, sub)
    end
end

""" `collect(tr::Trace)` returns a vector of all `Trace` objects within `tr`. """
function Base.collect(tr::Trace)
    trace_list = Trace[]
    collect_!(trace_list, tr)
    trace_list
end

""" `prune(tr::Trace, max_depth::Int=0, max_length::Int=1000000000)` prunes the Trace-tree
maximum tree depth, and maximum length (number of branches in each node).
(convenient to first explore a trace at a high-level).

`prune(tr)` prunes all of `tr`'s children (so all that remains is the function call). """
prune(tr::Trace, max_depth::Int=0, max_length::Int=1000000000) =
    Trace(tr, Trace[prune(sub_tr, max_depth-1, max_length)
                    for sub_tr in tr.called[1:min(length(tr), max_length)]
                    if max_depth > 0])

is_mutating(tr::Trace) = is_mutating(tr.func)
is_mutating(f::Function) = last(string(f)) == '!'
is_mutating(typ::Type) = typ

function is_inferred(tr::Trace)
    try
        apply_macro(:@inferred, tr, Base.Test)
        return true
    catch e
        if isa(e, ErrorException)&&contains(e.msg, "does not match inferred return type")
            return false
        else
            rethrow()
        end
    end
end
map_is_inferred(tr::Trace) = redgreen(map(is_inferred, tr))

number(x::Number) = x

################################################################################
# group_by

struct Group
    #key    # temporarily taken out until https://github.com/simonster/JLD2.jl/issues/37
            # is resolved
    traces::Vector{Trace}
end
Group(key, traces) = Group(traces)
value(grp::Group) = value(grp[1])

Base.getindex(gr::Group, i) = gr.traces[i]
Base.push!(gr::Group, tr::Trace) = push!(gr.traces, tr)
Base.length(gr::Group) = length(gr.traces)
Base.endof(gr::Group) = length(gr)
is_root(grp::Group) = is_root(grp[1])
(grp::Group)() = grp[1]()
function show_group(io, mime, gr)
    N = length(gr)
    write(io, "$N call", N>1?"s":"", " like ")
    show_call_(io, mime, gr[1])
end
# Necessary to split because otherwise it's ambiguous
Base.show(io::IO, mime::MIME"text/plain", gr::Group) = show_group(io, mime, gr)
function Base.show(io::IO, mime::MIME"text/html", gr::Group)
    write(io, "<pre>")
    show_group(io, mime, gr)
    write(io, "</pre>")
end
Base.show(io::IO, gr::Group) = show_group(io, MIME"text/plain"(), gr)

apply(mac_or_fun, group::Group) =
    # Only apply mac_or_fun to the first trace, which is assumed to be representative
    Group([map(mac_or_fun, group[1]); group[2:end]])

measure(mac_or_fun::Union{Expr, Function}, groups::Vector{Group}) =
    # Hmmm, not a great definition?
    map_groups(mac_or_fun, groups)

map_groups(mac_or_fun, groups::Vector{Group}) =
    Group[apply(return_error_fn(handle_mac(mac_or_fun)), grp) for grp in groups]
# map_groups(mac_or_fun, groups1::Vector{Group}, groups2::Vector{Group}) =
#     Group[apply(return_error_fn(handle_mac(mac_or_fun)), grp1, grp2)
#           for (grp1, grp2) in zip(groups1, groups2)]

if_not_error(fn::Function, group::Group) =
    value(group) isa Exception ? apply(tr->value(group), group) : fn()
if_not_error(fn::Function, group1::Group, group2::Group) =
    (value(group1) isa Exception || value(group2) isa Exception ?
     # Return the exception
     apply(tr->value(value(group1) isa Exception ? group1 : group2), group1) :
     fn())

macro delegate_to_group(funs...)
    esc(quote
        $([:($fun(grp::Group, args...; kwargs...) = apply($fun, grp, args...; kwargs...))
           for fun in funs]...)
        end)
end
is_inferred(grp::Group) = is_inferred(grp[1])

redgreen(grp::Group) = Group([redgreen(grp[1]); grp[2:end]])

function group_by(by::Function, trace::Trace)
    di = OrderedDict{Any, Group}()
    for tr in collect(trace)
        key = by(tr)
        group = get!(()->Group(key, []), di, key)
        push!(group, prune(tr))
    end
    collect(values(di))
end

################################################################################

include("scrub_stderr.jl")
include("code_update.jl")
include("tracing.jl")
include("show.jl")
include("specializations.jl")
include("profile.jl")
include("debug.jl")

""" `@traceable trace_log(args...; kwargs...) = nothing` is a dummy, always-traced
function you can call to add extra information to your traces. For instance,
`trace_log(:inside_foo; y=2, sqrt_of_x=sqrt(x))` """
@traceable trace_log(args...; kwargs...) = nothing

end # module


text:
# helpers for writing device functionality

# local method table for device functions
@static if isdefined(Base.Experimental, Symbol("@overlay"))
Base.Experimental.@MethodTable(method_table)
else
const method_table = nothing
end

# list of overrides (only for Julia 1.6)
const overrides = Expr[]

macro device_override(ex)
    code = quote
        $GPUCompiler.@override(CUDA.method_table, $ex)
    end
    if isdefined(Base.Experimental, Symbol("@overlay"))
        return esc(code)
    else
        push!(overrides, code)
        return
    end
end

macro device_function(ex)
    ex = macroexpand(__module__, ex)
    def = splitdef(ex)

    # generate a function that errors
    def[:body] = quote
        error("This function is not intended for use on the CPU")
    end

    esc(quote
        $(combinedef(def))
        @device_override $ex
    end)
end

macro device_functions(ex)
    ex = macroexpand(__module__, ex)

    # recursively prepend `@device_function` to all function definitions
    function rewrite(block)
        out = Expr(:block)
        for arg in block.args
            if Meta.isexpr(arg, :block)
                # descend in blocks
                push!(out.args, rewrite(arg))
            elseif Meta.isexpr(arg, [:function, :(=)])
                # rewrite function definitions
                push!(out.args, :(@device_function $arg))
            else
                # preserve all the rest
                push!(out.args, arg)
            end
        end
        out
    end

    esc(rewrite(ex))
end


text:
using RecipesBase
using ColorTypes
import PlotUtils: cgrad, palette, color_list
using LaTeXStrings

@userplot Trajectories

@recipe function f(h::Trajectories)
  if length(h.args) != 2
      error("`trajectories` should be given two arguments.  Got: $(typeof(h.args))")
  end
  traj_array, sys = h.args

  @series begin
    sys.bodies
  end

  for traj in traj_array
    @series begin
      linewidth := 2
      traj[1,:], traj[2,:]
    end
  end

  @series begin
    xlims --> (-1,3)
    ylims --> (-1,1)
    xguide := L"x"
    yguide := L"y"
    title := "Particle trajectories"
    aspect_ratio := 1
    size --> (800,400)
    ()
  end

end

@userplot FieldTrajectory

@recipe function f(h::FieldTrajectory;fieldlabel="Field",deriv=0)
  if length(h.args) != 3
      error("`fieldtrajectory` should be given three arguments.  Got: $(typeof(h.args))")
  end
  traj, field, sys = h.args

  xtraj = traj[1,:]
  ytraj = traj[2,:]

  body = sys.bodies[1]

  layout := (2,1)
  size --> (600,600)
  xlims := (-1,2)

  @series begin
    subplot := 1
    ylims --> (-0.5,0.5)
    xguide := L"x"
    yguide := L"y"
    aspect_ratio := 1
    body
  end

  @series begin
    subplot := 1
    aspect_ratio := 1
    ylims --> (-0.5,0.5)
    title := "Particle trajectory"
    xtraj, ytraj
  end

  yb = -20
  yt = 20
  Xr = [minimum(body.x), maximum(body.x), maximum(body.x), minimum(body.x), minimum(body.x)]
  Yr = [yb,yb,yt,yt,yb]

  @series begin
    subplot := 2
    fillrange := 0
    fillalpha --> 0.2
    linealpha --> 0.2
    fillcolor := :lightgray
    linecolor := :lightgray
    label := "Body location"
    Xr,Yr
  end

  if typeof(field) <: VectorGridData
    utraj,vtraj = field_along_trajectory(field,sys,traj,deriv=deriv)
    minuv = min(minimum(utraj),minimum(vtraj)) - 0.5
    maxuv = max(maximum(utraj),maximum(vtraj)) + 0.5

    @series begin
      subplot := 2
      xguide := L"x"
      label := L"$x$ %$fieldlabel"
      xtraj, utraj
    end

    @series begin
      subplot := 2
      xguide := L"x"
      label := L"$y$ %$fieldlabel"
      title := "$fieldlabel components along trajectory"
      legend := true
      ylims := (minuv,maxuv)
      xtraj, vtraj
    end
  elseif typeof(field) <: ScalarGridData
    straj = field_along_trajectory(field,sys,traj,deriv=deriv)

    @series begin
      subplot := 2
      xguide := L"x"
      label := "$fieldlabel"
      title := "$fieldlabel along trajectory"
      ylims := (minimum(straj)-0.5,maximum(straj)+0.5)
      legend := true
      xtraj, straj
    end
  end

end


text:
using Plots

x = collect(0:0.01:0.99)
y = 1 ./ (1 .- x)

plot(x, y, xlim=[0, 1], dpi=600, xlabel="Utilization", ylabel="Normalized Waiting Time", legend = false)
png("waitingtime.png")


text:
module Factories

using MacroTools: block
using MacroTools: splitarg, combinefield, combinearg
using MacroTools: splitstructdef, combinestructdef
using MacroTools: splitdef, combinedef
using Printf: @printf, @sprintf
using ..TypeTraits: efficientoperations

import MacroTools: rmlines

export FExpr, escape
export Escaped, UnEscaped, MixEscaped
export RawExpr, rawexpr
export Inference, @inference
export Argument, @argument
export Parameter, @parameter
export Field, @field
export Block, @block, @push!, rmlines!, @rmlines!, rmlines, @rmlines
export FunctionFactory, @functionfactory, addargs!, @addargs!, addkwargs!, @addkwargs!, addwhereparams!, @addwhereparams!, extendbody!, @extendbody!
export TypeFactory, @typefactory, addfields!, @addfields!, addconstructors!, @addconstructors!
export addparams!, @addparams!

@generated function maxfieldnamelength(::T) where T
    result = 1
    for name in fieldnames(T)
        result = max(result, name|>string|>length)
    end
    return :($result)
end

"Factory expression types, which is defined as `Union{Symbol, Expr}`."
const FExpr = Union{Symbol, Expr}

"Abstract escape mechanism."
abstract type EscapeMechanism end

"""
    Escaped(names::Symbol...)

Indicate that symbols of a factory should be escaped if they are in `names`.
"""
struct Escaped{N} <: EscapeMechanism
    names::NTuple{N, Symbol}
    Escaped(names::Symbol...) = new{length(names)}(names)
end

"""
    UnEscaped(names::Symbol...)

IIndicate that symbols of a factory should be escaped if they are not in `names`.
"""
struct UnEscaped{N} <: EscapeMechanism
    names::NTuple{N, Symbol}
    UnEscaped(names::Symbol...) = new{length(names)}(names)
end

"""
    MixEscaped(escaped::Escaped)
    MixEscaped(unescaped::UnEscaped)
    MixEscaped(escaped::Escaped, unescaped::UnEscaped)
    MixEscaped(unescaped::UnEscaped, escaped::Escaped)

Indicate that some parts of a factory use the `Escaped` mechanism while other parts use the `UnEscaped` mechanism.
"""
struct MixEscaped{N, M} <: EscapeMechanism
    escaped::Escaped{N}
    unescaped::UnEscaped{M}
end
MixEscaped() = MixEscaped(Escaped(), UnEscaped())
MixEscaped(escaped::Escaped) = MixEscaped(escaped, UnEscaped())
MixEscaped(unescaped::UnEscaped) = MixEscaped(Escaped(), unescaped)
MixEscaped(unescaped::UnEscaped, escaped::Escaped) = MixEscaped(escaped, unescaped)

"Raw expression without any variable escaped."
struct RawExpr <: EscapeMechanism end
"""
    rawexpr

Indicate that no variable in a factory should be escaped.
"""
const rawexpr = RawExpr()

"""
    escape(expr, ::RawExpr) -> Any
    escape(expr, ::Escaped) -> Any
    escape(expr, ::UnEscaped) -> Any
    escape(expr::Symbol, em::Escaped) -> FExpr
    escape(expr::Expr, em::Escaped) -> Expr
    escape(expr::Symbol, em::UnEscaped) -> FExpr
    escape(expr::Expr, em::UnEscaped) -> Expr

Escape the variables in the input expression.
"""
escape(expr, ::RawExpr) = expr
escape(expr, ::Escaped) = expr
escape(expr, ::UnEscaped) = expr
escape(expr::Symbol, em::Escaped) = (expr ∈ em.names) ? :($(esc(expr))) : expr
escape(expr::Expr, em::Escaped) = (expr.head == :escape) ? expr : Expr(expr.head, (escape(arg, em) for arg in expr.args)...)
escape(expr::Symbol, em::UnEscaped) = (expr ∉ em.names) ? :($(esc(expr))) : expr
function escape(expr::Expr, em::UnEscaped)
    if expr.head == :escape
        @assert length(expr.args) == 1 "escape error: internal error..."
        escape(expr.args[1], em)
    else
        Expr(expr.head, (escape(arg, em) for arg in expr.args)...)
    end
end

"""
    names(expr::Symbol) -> Vector{Symbol}
    names(expr::Expr) -> Vector{Symbol}

Get all the symbols in an expression.
"""
Base.names(expr::Symbol) = [expr]
function Base.names(expr::Expr)
    result = Symbol[]
    for arg in expr.args
        (isa(arg, Symbol) || isa(arg, Expr)) && append!(result, names(arg))
    end
    return result
end

"""
    AbstractFactory

Abstract type for all concrete factories.
"""
abstract type AbstractFactory end

"""
    ==(f1::F, f2::F) where F<:AbstractFactory -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(f1::F, f2::F) where {F<:AbstractFactory} = ==(efficientoperations, f1, f2)

"""
    isequal(f1::F, f2::F) where F<:AbstractFactory -> Bool

Overloaded equivalent operator.
"""
Base.isequal(f1::F, f2::F) where {F<:AbstractFactory} = isequal(efficientoperations, f1, f2)

"""
    show(io::IO, f::AbstractFactory)

Show a concrete `AbstractFactory`.
"""
function Base.show(io::IO, f::AbstractFactory)
    @printf io "%s(\n" nameof(typeof(f))
    for name in fieldnames(typeof(f))
        nstr = " "^2 * string(name) * ":" * " "^(maxfieldnamelength(f)+1-length(string(name)))
        value = getfield(f, name)
        if isa(value, Nothing)
            vstr = "nothing"
        elseif isa(value, AbstractFactory)
            vstr = string(value(rawexpr))
        elseif isa(value, Vector{<:AbstractFactory})
            eltypename = nameof(eltype(value))
            strings = Tuple(string(af(rawexpr)) for af in value)
            if ((length(strings) > 0) ? sum(map(length, strings)) : 0) < 90-maxfieldnamelength(f)
                vstr = @sprintf "%s[%s]" eltypename join(strings, ", ")
            else
                vstr = replace(@sprintf("%s[\n%s\n]", eltypename, join(strings, "\n")), "\n"=>"\n  ")
            end
        else
            vstr = string(value)
        end
        ('\n' ∈ vstr) && (vstr = replace(vstr, "\n"=>"\n"*" "^(maxfieldnamelength(f)+4)))
        @printf(io, "%s%s\n", nstr, vstr)
    end
    @printf io ")"
end

"""
    replace(f::AbstractFactory; kwargs...) -> typeof(f)

Return a copy of a concrete `AbstractFactory` with some of the field values replaced by the keyword arguments.
"""
Base.replace(f::AbstractFactory; kwargs...) = replace(efficientoperations, f; kwargs...)

"""
    Inference(head::Union{Symbol, Nothing}, name::Union{Symbol, Nothing}, params::Union{Inference, Vector{Inference}, Nothing})
    Inference(;
            head::Union{Symbol, Nothing}=nothing,
            name::Union{Symbol, Nothing}=nothing,
            params::Union{Inference, Vector{Inference}, Nothing}=nothing,
            )
    Inference(expr::FExpr)

The struct to describe a type inference.
"""
mutable struct Inference <: AbstractFactory
    head::Union{Symbol, Nothing}
    name::Union{Symbol, Nothing}
    params::Union{Inference, Vector{Inference}, Nothing}
    function Inference(head::Union{Symbol, Nothing}, name::Union{Symbol, Nothing}, params::Union{Inference, Vector{Inference}, Nothing})
        @assert (head === nothing) == (params === nothing) "Inference error: `head` and `params` must be or not be `nothing` at the same time."
        new(head, name, params)
    end
end
function Inference(; head::Union{Symbol, Nothing}=nothing, name::Union{Symbol, Nothing}=nothing, params::Union{Inference, Vector{Inference}, Nothing}=nothing)
    Inference(head, name, params)
end
function Inference(expr::FExpr)
    if isa(expr, Symbol)
        return Inference(nothing, expr, nothing)
    elseif expr.head == :(<:)
        @assert length(expr.args) == 1 "Inference error: wrong input expr."
        return Inference(expr.head, nothing, Inference(expr.args[1]))
    elseif expr.head == :curly
        return Inference(expr.head, expr.args[1], Inference.(@views(expr.args[2:end])))
    else
        error("Inference error: wrong input expr.")
    end
end

"""
    @inference expr::FExpr

Construct an `Inference` directly from a type inference.
"""
macro inference(expr::FExpr)
    expr = [expr]
    return :(Inference($expr...))
end

"""
    (i::Inference)(em::RawExpr) -> FExpr
    (i::Inference)(em::UnEscaped) -> FExpr
    (i::Inference)(em::MixEscaped) -> FExpr

Convert a `Inference` to the `Expr` representation of the type inference it describes.
"""
(i::Inference)(em::MixEscaped) = i(em.unescaped)
function (i::Inference)(em::Union{RawExpr, <:UnEscaped})
    if i.head === nothing
        return escape(i.name, em)
    elseif i.head == :(<:)
        return Expr(i.head, i.params(em))
    else
        return Expr(i.head, escape(i.name, em), [param(em) for param in i.params]...)
    end
end

"""
    Argument(name::Union{Symbol, Nothing}, type::Inference, slurp::Bool, default::Any)
    Argument(; name::Union{Symbol, Nothing}=nothing, type::Inference=Inference(:Any), slurp::Bool=false, default::Any=nothing)
    Argument(expr::FExpr)

The struct to describe a argument of a `function`.
"""
mutable struct Argument <: AbstractFactory
    name::Union{Symbol, Nothing}
    type::Inference
    slurp::Bool
    default::Any
end
Argument(; name::Union{Symbol, Nothing}=nothing, type::Inference=Inference(:Any), slurp::Bool=false, default::Any=nothing) = Argument(name, type, slurp, default)
Argument(expr::FExpr) = (cache = splitarg(expr); Argument(cache[1], Inference(cache[2]), cache[3], cache[4]))

"""
    @argument expr::FExpr

Construct an `Argument` directly from an argument statement.
"""
macro argument(expr::FExpr)
    expr = [expr]
    return :(Argument($expr...))
end

"""
    (a::Argument)(em::RawExpr) -> Expr
    (a::Argument)(em::MixEscaped) -> Expr

Convert an `Argument` to the `Expr` representation of the argument it describes.
"""
(a::Argument)(em::RawExpr) = combinearg(a.name, a.type(em), a.slurp, escape(a.default, em))
(a::Argument)(em::MixEscaped) = combinearg(a.name, a.type(em.unescaped), a.slurp, escape(a.default, em.escaped))

"""
    Parameter(head::Union{Symbol, Nothing}, name::Union{Symbol, Nothing}, type::Union{Inference, Nothing})
    Parameter(; head::Union{Symbol, Nothing}=nothing, name::Union{Symbol, Nothing}=nothing, type::Union{Inference, Nothing}=nothing)
    Parameter(expr::FExpr)

The struct to describe a parameter of a `function` or a `type`.
"""
mutable struct Parameter <: AbstractFactory
    head::Union{Symbol, Nothing}
    name::Union{Symbol, Nothing}
    type::Union{Inference, Nothing}
end
Parameter(; head::Union{Symbol, Nothing}=nothing, name::Union{Symbol, Nothing}=nothing, type::Union{Inference, Nothing}=nothing) = Parameter(head, name, type)
function Parameter(expr::FExpr)
    if isa(expr, Symbol)
        head, name, type = nothing, expr, nothing
    elseif expr.head == :(<:)
        head = :(<:)
        (name, type) = (length(expr.args) == 2) ? expr.args : (nothing, expr.args[1])
        type = Inference(type)
    elseif expr.head == :(::)
        @assert length(expr.args) == 1 "Parameter error: wrong input expr."
        head, name, type = :(::), nothing, Inference(expr.args[1])
    else
        error("Parameter error: wrong input expr.")
    end
    Parameter(head, name, type)
end

"""
    @parameter expr::FExpr

Construct a `Parameter` directly from an parameter statement.
"""
macro parameter(expr::FExpr)
    expr = [expr]
    return :(Parameter($expr...))
end

"""
    (p::Parameter)(em::RawExpr) -> FExpr
    (p::Parameter)(em::UnEscaped) -> FExpr
    (p::Parameter)(em::MixEscaped) -> FExpr

Convert a `Parameter` to the `Expr` representation of the parameter it describes.
"""
function (p::Parameter)(em::Union{RawExpr, <:UnEscaped, <:MixEscaped})
    (p.name === nothing) ?
        ((p.head == :(::)) ? p.type(em) : Expr(:(<:), p.type(em))) :
        ((p.type === nothing) ? p.name : Expr(:(<:), p.name, p.type(em)))
end

"""
    Field(name::Symbol, type::Inference)
    Field(; name::Symbol, type::FExpr=Inference(:Any))
    Field(expr::FExpr)

The struct to describe a field of a `struct`.
"""
mutable struct Field <: AbstractFactory
    name::Symbol
    type::Inference
end
Field(; name::Symbol, type::Inference=Inference(:Any)) = Field(name, type)
Field(expr::FExpr) = (cache = splitarg(expr); Field(cache[1], Inference(cache[2])))

"""
    @field expr::FExpr

Construct a `Field` directly from a field statement.
"""
macro field(expr::FExpr)
    expr = [expr]
    return :(Field($expr...))
end

"""
    (f::Field)(em::RawExpr) -> Expr
    (f::Field)(em::UnEscaped) -> Expr
    (f::Field)(em::MixEscaped) -> Expr

Convert a `Field` to the `Expr` representation of the field it describes.
"""
(f::Field)(em::Union{RawExpr, <:UnEscaped, <:MixEscaped}) = combinefield((f.name, f.type(em)))

"""
    Block(parts::FExpr...)

The struct to describe a `begin ... end` block.
"""
mutable struct Block <: AbstractFactory
    body::Vector{Any}
    Block(parts::FExpr...) = new(vcat((block(part).args for part in parts)...))
end

"""
    @block parts::FExpr...

Construct a `Block` directly from a `begin ... end` block definition.
"""
macro block(parts::FExpr...)
    return :(Block($parts...))
end

"""
    (b::Block)(em::RawExpr) -> Expr
    (b::Block)(em::Escaped) -> Expr
    (b::Block)(em::MixEscaped) -> Expr

Convert a `Block` to the `Expr` representation of the `begin ... end` block it describes.
"""
(b::Block)(em::MixEscaped) = b(em.escaped)
(b::Block)(em::Union{RawExpr, <:Escaped}) = Expr(:block, (escape(part, em) for part in b.body)...)

"""
    push!(b::Block, parts::FExpr...) -> Block
    push!(b::Block, parts::Block...) -> Block

Push other parts into the body of a block.
"""
Base.push!(b::Block) = b
Base.push!(b::Block, parts::FExpr...) = push!(b, Block(parts...))
function Base.push!(b::Block, parts::Block...)
    for part in parts
        append!(b.body, part.body)
    end
    return b
end

"""
    @push! b parts::FExpr...

Push other parts into the body of a block.
"""
macro push!(b, parts::FExpr...)
    return :(push!($(esc(b)), $parts...))
end

"""
    rmlines!(b::Block) -> Block

Remove line number nodes in the body of a block.
"""
rmlines!(b::Block) = (filter!(part->!isa(part, LineNumberNode), b.body); b)

"""
    @rmlines! b::Expr

Remove line number nodes in the body of a block.
"""
macro rmlines!(b::Expr)
    return :(rmlines!($(esc(b))))
end

"""
    rmlines(b::Block) -> Block

Return a copy of a block with the line number nodes removed.
"""
rmlines(b::Block) = Block(filter(part->!isa(part, LineNumberNode), b.body)...)

"""
    @rmlines b::Expr

Return a copy of a block with the line number nodes removed.
"""
macro rmlines(b::Expr)
    return :(rmlines($(esc(b))))
end

"""
    FunctionFactory(name::FExpr,
                    params::Vector{Inference},
                    args::Vector{Argument},
                    kwargs::Vector{Argument},
                    rtype::Inference,
                    whereparams::Vector{Parameter},
                    body::Block
                    )
    FunctionFactory(;   name::FExpr,
                        params::Vector{Inference}=Inference[],
                        args::Vector{Argument}=Argument[],
                        kwargs::Vector{Argument}=Argument[],
                        rtype::Inference=Inference(:Any),
                        whereparams::Vector{Parameter}=Parameter[],
                        body::Block=Block()
                        )
    FunctionFactory(expr::Expr)

The struct to describe a `function`.
"""
mutable struct FunctionFactory <: AbstractFactory
    name::FExpr
    params::Vector{Inference}
    args::Vector{Argument}
    kwargs::Vector{Argument}
    rtype::Inference
    whereparams::Vector{Parameter}
    body::Block
end
function FunctionFactory(;  name::FExpr,
                            params::Vector{Inference}=Inference[],
                            args::Vector{Argument}=Argument[],
                            kwargs::Vector{Argument}=Argument[],
                            rtype::Inference=Inference(:Any),
                            whereparams::Vector{Parameter}=Parameter[],
                            body::Block=Block()
                            )
    FunctionFactory(name, params, args, kwargs, rtype, whereparams, body)
end
function FunctionFactory(expr::Expr)
    dict = splitdef(expr)
    FunctionFactory(dict[:name],
                    Inference.(get(dict, :params, [])),
                    Argument.(dict[:args]),
                    Argument.(dict[:kwargs]),
                    Inference(get(dict, :rtype, :Any)),
                    Parameter.(collect(FExpr, dict[:whereparams])),
                    Block(dict[:body])
                    )
end

"""
    @functionfactory expr::FExpr

Construct a `FunctionFactory` directly from a function definition.
"""
macro functionfactory(expr::Expr)
    expr = [expr]
    return :(FunctionFactory($expr...))
end

"""
    (ff::FunctionFactory)(em::RawExpr) -> Expr
    (ff::FunctionFactory)(em::MixEscaped) -> Expr

Convert a `FunctionFactory` to the `Expr` representation of the `function` it describes.
"""
function (ff::FunctionFactory)(em::Union{RawExpr, <:MixEscaped})
    combinedef(
        Dict(:name          =>      escape(ff.name, isa(em, RawExpr) ? em : em.escaped),
            :params         =>      [param(em) for param in ff.params],
            :args           =>      [arg(em) for arg in ff.args],
            :kwargs         =>      [kwarg(em) for kwarg in ff.kwargs],
            :rtype          =>      ff.rtype(em),
            :whereparams    =>      [whereparam(em) for whereparam in ff.whereparams],
            :body           =>      ff.body(em)
            )
        )
end

"""
    addargs!(ff::FunctionFactory, args::Argument...) -> FunctionFactory
    addargs!(ff::FunctionFactory, args::FExpr...) -> FunctionFactory

Add a couple of positional arguments to a function factory.
"""
addargs!(ff::FunctionFactory) = ff
addargs!(ff::FunctionFactory, args::Argument...) = (push!(ff.args, args...); ff)
addargs!(ff::FunctionFactory, args::FExpr...) = addargs!(ff, Argument.(args)...)

"""
    @addargs! ff args::FExpr...

Add a couple of positional arguments to a function factory.
"""
macro addargs!(ff, args::FExpr...)
    return :(addargs!($(esc(ff)), $args...))
end

"""
    addkwargs!(ff::FunctionFactory, kwargs::Argument...) -> FunctionFactory
    addkwargs!(ff::FunctionFactory, kwargs::FExpr...) -> FunctionFactory

Add a couple of keyword arguments to a function factory.
"""
addkwargs!(ff::FunctionFactory) = ff
addkwargs!(ff::FunctionFactory, kwargs::Argument...) = (push!(ff.kwargs, kwargs...); ff)
addkwargs!(ff::FunctionFactory, kwargs::FExpr...) = addkwargs!(ff, Argument.(kwargs)...)

"""
    @addkwargs! ff kwargs::FExpr...

Add a couple of keyword arguments to a function factory.
"""
macro addkwargs!(ff, kwargs::FExpr...)
    return :(addkwargs!($(esc(ff)), $kwargs...))
end

"""
    extendbody!(ff::FunctionFactory, parts::FExpr...) -> FunctionFactory
    extendbody!(ff::FunctionFactory, parts::Block...) -> FunctionFactory

Extend the body of a function factory.
"""
extendbody!(ff::FunctionFactory) = ff
extendbody!(ff::FunctionFactory, parts::FExpr...) = extendbody!(ff, Block(parts...))
extendbody!(ff::FunctionFactory, parts::Block...) = (push!(ff.body, parts...); ff)

"""
    @extendbody! ff parts::FExpr...

Extend the body of a function factory.
"""
macro extendbody!(ff, parts::FExpr...)
    return :(extendbody!($(esc(ff)), $parts...))
end

"""
    addwhereparams!(ff::FunctionFactory, whereparams::Parameter...) -> FunctionFactory
    addwhereparams!(ff::FunctionFactory, whereparams::FExpr...) -> FunctionFactory

Add a couple of method where parameters to a function factory or a type factory.
"""
addwhereparams!(ff::FunctionFactory) = ff
addwhereparams!(ff::FunctionFactory, whereparams::Parameter...) = (push!(ff.whereparams, whereparams...); ff)
addwhereparams!(ff::FunctionFactory, whereparams::FExpr...) = addwhereparams!(ff, Parameter.(whereparams)...)

"""
    @addwhereparams! ff whereparams::FExpr...

Add a couple of method parameters to a function factory or a type factory.
"""
macro addwhereparams!(ff, whereparams::FExpr...)
    return :(addwhereparams!($(esc(ff)), $whereparams...))
end

"""
    TypeFactory(name::Symbol, mutable::Bool, params::Vector{Parameter}, supertype::Inference, fields::Vector{Field}, constructors::Vector{FunctionFactory})
    TypeFactory(;   name::Symbol,
                    mutable::Bool=false,
                    params::Vector{Parameter}=Parameter[],
                    supertype::Inference=Inference(:Any),
                    fields::Vector{Field}=Field[],
                    constructors::Vector{FunctionFactory}=FunctionFactory[]
                    )
    TypeFactory(expr::Expr)

The struct to describe a `struct`.
"""
mutable struct TypeFactory <: AbstractFactory
    name::Symbol
    mutable::Bool
    params::Vector{Parameter}
    supertype::Inference
    fields::Vector{Field}
    constructors::Vector{FunctionFactory}
end
function TypeFactory(;  name::Symbol,
                        mutable::Bool=false,
                        params::Vector{Parameter}=Parameter[],
                        supertype::Inference=Inference(:Any),
                        fields::Vector{Field}=Field[],
                        constructors::Vector{FunctionFactory}=FunctionFactory[]
                        )
    TypeFactory(name, mutable, params, supertype, fields, constructors)
end
function TypeFactory(expr::Expr)
    dict = splitstructdef(expr)
    TypeFactory(dict[:name],
                dict[:mutable],
                Parameter.(dict[:params]),
                Inference(dict[:supertype]),
                [Field(field[1], Inference(field[2])) for field in dict[:fields]],
                FunctionFactory.(dict[:constructors])
                )
end

"""
    @typefactory expr::Expr

Construct a `TypeFactory` directly from a type definition.
"""
macro typefactory(expr::Expr)
    expr = [expr]
    return :(TypeFactory($expr...))
end

"""
    (tf::TypeFactory)(em::RawExpr) -> Expr
    (tf::TypeFactory)(em::MixEscaped) -> Expr

Convert a `TypeFactory` to the `Expr` representation of the `struct` it describes.
"""
function (tf::TypeFactory)(em::Union{RawExpr, <:MixEscaped})
    combinestructdef(
        Dict(:name          =>      escape(tf.name, isa(em, RawExpr) ? em : em.escaped),
            :mutable        =>      tf.mutable,
            :params         =>      [param(em) for param in tf.params],
            :supertype      =>      tf.supertype(em),
            :fields         =>      [(field.name, field.type(em)) for field in tf.fields],
            :constructors   =>      [constructor(em) for constructor in tf.constructors]
            )
    )
end

"""
    addfields!(tf::TypeFactory, fields::Field...) -> TypeFactory
    addfields!(tf::TypeFactory, fields::FExpr...) -> TypeFactory

Add a couple of fields to a type factory.
"""
addfields!(tf::TypeFactory) = tf
addfields!(tf::TypeFactory, fields::Field...) = (push!(tf.fields, fields...); tf)
addfields!(tf::TypeFactory, fields::FExpr...) = addfields!(tf, Field.(fields)...)

"""
    @addfields! tf fields::FExpr...

Add a couple of fields to a type factory.
"""
macro addfields!(tf, fields::FExpr...)
    return :(addfields!($(esc(tf)), $fields...))
end

"""
    addconstructors!(tf::TypeFactory, constructors::FunctionFactory...) -> TypeFactory
    addconstructors!(tf::TypeFactory, constructors::Expr...) -> TypeFactory

Add a couple of constructors to a type factory.
"""
addconstructors!(tf::TypeFactory) = tf
addconstructors!(tf::TypeFactory, constructors::FunctionFactory...) = (push!(tf.constructors, constructors...); tf)
addconstructors!(tf::TypeFactory, constructors::Expr...) = addconstructors!(tf, FunctionFactory.(constructors)...)

"""
    @addconstructors! tf constructors::Expr...

Add a couple of constructors to a type factory.
"""
macro addconstructors!(tf, constructors::Expr...)
    return :(addconstructors!($(esc(tf)), $constructors...))
end

"""
    addparams!(f::FunctionFactory, params::Inference...) ->FunctionFactory
    addparams!(f::FunctionFactory, params::FExpr...) -> FunctionFactory
    addparams!(f::TypeFactory, params::Parameter...) -> TypeFactory
    addparams!(f::TypeFactory, params::FExpr...) -> TypeFactory

Add a couple of parameters to a function factory or a type factory.
"""
addparams!(f::FunctionFactory) = f
addparams!(f::FunctionFactory, params::Inference...) = (push!(f.params, params...); f)
addparams!(f::FunctionFactory, params::FExpr...) = addparams!(f, Inference.(params)...)
addparams!(f::TypeFactory) = f
addparams!(f::TypeFactory, params::Parameter...) = (push!(f.params, params...); f)
addparams!(f::TypeFactory, params::FExpr...) = addparams!(f, Parameter.(params)...)

"""
    @addparams! f params::FExpr...

Add a couple of method parameters to a function factory or a type factory.
"""
macro addparams!(f, params::FExpr...)
    return :(addparams!($(esc(f)), $params...))
end

end #module


text:
mutable struct BackoffSchedulerEntry
    match_limit::Int
    ban_length::Int
    times_banned::Int
    banned_until::Int
end

"""
A Rewrite Scheduler that implements exponential rule backoff.
For each rewrite, there exists a configurable initial match limit.
If a rewrite search yield more than this limit, then we ban this rule
for number of iterations, double its limit, and double the time it
will be banned next time.

This seems effective at preventing explosive rules like
associativity from taking an unfair amount of resources.
"""
mutable struct BackoffScheduler <: AbstractScheduler
    data::IdDict{AbstractRule, BackoffSchedulerEntry}
    G::EGraph
    theory::Vector{<:AbstractRule}
    curr_iter::Int
end

cansearch(s::BackoffScheduler, r::AbstractRule)::Bool = s.curr_iter > s.data[r].banned_until

function BackoffScheduler(g::EGraph, theory::Vector{<:AbstractRule})
    # BackoffScheduler(g, theory, 128, 4)
    BackoffScheduler(g, theory, 1000, 5)
end

function BackoffScheduler(G::EGraph, theory::Vector{<:AbstractRule}, match_limit::Int, ban_length::Int)
    gsize = length(G.uf)
    data = IdDict{AbstractRule, BackoffSchedulerEntry}()

    for rule ∈ theory
        data[rule] = BackoffSchedulerEntry(match_limit, ban_length, 0, 0)
    end

    return BackoffScheduler(data, G, theory, 1)
end

# can saturate if there's no banned rule
cansaturate(s::BackoffScheduler)::Bool = all(kv -> s.curr_iter > last(kv).banned_until, s.data)

function inform!(s::BackoffScheduler, rule::AbstractRule, n_matches)
    # println(s.data[rule])

    rd = s.data[rule]
    treshold = rd.match_limit << rd.times_banned
    if n_matches > treshold
        ban_length = rd.ban_length << rd.times_banned
        rd.times_banned += 1
        rd.banned_until = s.curr_iter + ban_length
        # @info "banning rule $rule until $(rd.banned_until)!"
        return false
    end
    return true
end

function setiter!(s::BackoffScheduler, curr_iter)
    s.curr_iter = curr_iter
end


text:
module Visual
using GLFW
using Makie
using GeometryTypes

struct PsychoScreen
    scene
    monitor::GLFW.Monitor
    distance::Float64
end

PsychoSreen(scene, distance::Float64) = PsychoScreen(scene, GLFW.GetPrimaryMonitor(), distance)

function px2deg(screen::PsychoScreen, px)
    widths = GLFW.GetMonitorPhysicalSize(screen.monitor)
    r = sqrt(widths[1]^2 + widths[2]^2)
    dg = r/screen.distance
    vm = GLFW.GetVideoMode(screen.monitor)
    dpx = sqrt(vm.height^2 + vm.width^2)
    px*dg/dpx
end

function get_screen(resolution, distance, monitor=GLFW.GetPrimaryMonitor())
    scene = Scene(resolution=resolution,raw=true, show_axis=false)
    campixel!(scene)
    PsychoScreen(scene, monitor, distance)
end

function draw_circle(screen::PsychoScreen, pos::Point2f0, radius::Float32, color="black")
    points = decompose(Point2f0, Sphere(pos,radius))
    draw_poly(screen, points, color)
    nothing
end

function draw_square(screen::PsychoScreen, pos::Point2f0, width::Float32, color="black")
    points = [pos, pos + Point2f0(width, 0.0), pos + Point2f0(width, width), pos + Point2f0(0.0, width), pos]
    draw_poly(screen, points, color)
    nothing
end

function draw_poly(screen, points::AbstractVector{Point2f0}, color)
    poly!(screen.scene, points, show_axis=false, raw=true, color=color)
end

#Makie.opengl_renderloop[] = (screen)-> begin
#    GLFW.PollEvents()
#    Makie.render_frame(screen)
#    GLFW.SwapBuffers(Makie.to_native(screen))
#    glFinish()
#end

end# module


text:
import PerronFrobenius: AbstractTriangulationInvariantMeasure
import CausalityToolsBase: RectangularBinning, CustomReconstruction
import StateSpaceReconstruction: Simplex, generate_interior_points
import StaticArrays: SVector
import Distances: Metric, Chebyshev

export transferentropy,
    BinningTransferEntropyEstimator,
    TransferEntropyEstimator,
        TransferOperatorGrid,
        VisitationFrequency,
        NearestNeighbourMI

"""
    TransferEntropyEstimator

An abstract type for transfer entropy estimators. This type has several concrete subtypes
that are accepted as inputs to the [`transferentropy`](@ref) methods.

- [`VisitationFrequency`](@ref)
- [`TransferOperatorGrid`](@ref)
- [`NearestNeighbourMI`](@ref)
"""
abstract type TransferEntropyEstimator end

function Base.show(io::IO, estimator::TransferEntropyEstimator)
    s = "$(typeof(estimator))($(estimator.b))"
    print(io, s)
end

"""
    BinningTransferEntropyEstimator <: TransferEntropyEstimator

An abstract type for transfer entropy estimators that works on a discretization
of the [reconstructed state space](@ref custom_delay_reconstruction). Has the following concrete subtypes

- [`VisitationFrequency`](@ref)
- [`TransferOperatorGrid`](@ref)

## Used by

Concrete subtypes are accepted as inputs by

- [`transferentropy`](@ref te_estimator_rectangular) (low-level method)
- [`VisitationFrequencyTest`](@ref)
- [`TransferOperatorGridTest`](@ref),
- [`ExactSimplexIntersectionTest`](@ref)
- [`ApproximateSimplexIntersectionTest`](@ref)
"""
abstract type BinningTransferEntropyEstimator end

"""
    NearestNeighbourMI(k1::Int = 2, k2::Int = 3, metric::Metric = Chebyshev, b::Number)

A transfer entropy estimator counting of nearest neighbours nearest neighbours
to estimate mutual information over an appropriate
[custom delay reconstruction](@ref custom_delay_reconstruction) (the method from [1], as implemented in [2]).

## Fields

- **`k1::Int = 2`**: The number of nearest neighbours for the highest-dimensional mutual
    information estimate. To minimize bias, choose ``k_1 < k_2`` if
    ``min(k_1, k_2) < 10`` (see fig. 16 in [1]). Beyond dimension 5, choosing
    ``k_1 = k_2`` results in fairly low bias, and a low number of nearest
    neighbours, say `k1 = k2 = 4`, will suffice.
- **`k2::Int = 3`**: The number of nearest neighbours for the lowest-dimensional mutual
    information estimate. To minimize bias, choose ``k_1 < k_2`` if
    if ``min(k_1, k_2) < 10`` (see fig. 16 in [1]). Beyond dimension 5, choosing
    ``k_1 = k_2`` results in fairly low bias, and a low number of nearest
    neighbours, say `k1 = k2 = 4`, will suffice.
- **`metric::Metric = Chebyshev()`**: The metric used for distance computations.
- **`b::Number = 2`**: The base of the logarithm, controlling the unit of the transfer
    entropy estimate (e.g. `b = 2` will give the transfer entropy in bits).

## Used by

This estimator is accepted as input by

- [`transferentropy`](@ref te_estimator_nn) (low level method),
- [`NearestNeighbourMITest`](@ref)

## References

1. Kraskov, Alexander, Harald Stögbauer, and Peter Grassberger. "Estimating
    mutual information." Physical review E 69.6 (2004): 066138.
2. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
@Base.kwdef struct NearestNeighbourMI <: TransferEntropyEstimator
    k1::Int = 2
    k2::Int = 3
    metric::Metric = Chebyshev()
    b::Number = 2
end

function Base.show(io::IO, estimator::NearestNeighbourMI)
    k1 = estimator.k1
    k2 = estimator.k2
    metric = estimator.metric
    b = estimator.b
    s = "$(typeof(estimator))(b=$(b), k1=$(k1), k2=$(k2), metric=$(typeof(metric)))"
    print(io, s)
end

"""
    TransferOperatorGrid(; b::Number = 2)

An transfer entropy estimator which computes transfer entropy over a
dicretization of an appropriate [delay reconstruction](@ref custom_delay_reconstruction),
using the logarithm to base `b`. Invariant probabilities over the
partition are computed using an approximation to the transfer (Perron-Frobenius)
operator over the grid [1], which explicitly gives the transition probabilities
between states.

## Fields

- **`b::Number = 2`**: The base of the logarithm, controlling the unit of the transfer
    entropy estimate (e.g. `b = 2` will give the transfer entropy in bits).

## Used by

This estimator is accepted as input by

- [`transferentropy`](@ref te_estimator_rectangular) (low-level method)
- [`TransferOperatorGridTest`](@ref)
- [`ExactSimplexIntersectionTest`](@ref)
- [`ApproximateSimplexIntersectionTest`](@ref)

## References

[1] Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
Base.@kwdef struct TransferOperatorGrid <: BinningTransferEntropyEstimator
    """ The base of the logarithm usen when computing transfer entropy. """
    b::Number = 2.0

    TransferOperatorGrid(b) = new(b)
end

"""
    VisitationFrequency(; b::Number = 2)

An transfer entropy estimator which computes transfer entropy over a
dicretization of an appropriate [delay reconstruction](@ref custom_delay_reconstruction), using the
logarithm to the base `b`. The invariant probabilities over the partition
are computed using an approximation to the transfer (Perron-Frobenius) operator
over the grid [1], which explicitly gives the transition probabilities between states.

## Used by

This estimator is accepted as input by

- [`transferentropy`](@ref te_estimator_rectangular) (low-level method)
- [`VisitationFrequencyTest`](@ref)
- [`ExactSimplexIntersectionTest`](@ref)
- [`ApproximateSimplexIntersectionTest`](@ref)

## References

[1] Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
Base.@kwdef struct VisitationFrequency <: BinningTransferEntropyEstimator
    """ The base of the logarithm usen when computing transfer entropy. """
    b::Number = 2.0
end

"""
    transferentropy(μ::AbstractTriangulationInvariantMeasure, vars::TEVars,
        binning_scheme::RectangularBinning;
        estimator::BinningTransferEntropyEstimator = VisitationFrequency(),
        n::Int = 10000) -> Float64

#### Transfer entropy using a precomputed invariant measure over a triangulated partition

Estimate transfer entropy from an invariant measure over a triangulation of
an appropriate [generalised delay reconstruction](@ref custom_delay_reconstruction).

The invariant measure has been precomputed either as

1. `μ = invariantmeasure(pts, TriangulationBinning(), ApproximateIntersection())`, or
2. `μ = invariantmeasure(pts, TriangulationBinning(), ExactIntersection())`

`μ` contains all the information needed to compute transfer entropy, and the estimate
of the measure will be slightly different depending on the method. Approximate simplex
intersections is faster, while exact simplex intersections are (very) slow.

Note: `pts` must be a vector of states, not a vector of
variables/(time series). Wrap your time series in a `Dataset`
first if the latter is the case.

#### Computing transfer entropy (triangulation -> rectangular partition)

Computing transfer entropy directly on a triangulation is not possible.
To compute marginals, we need a rectangular grid.

Therefore, first we create a cloud of points (approximately `n` points)
with known measure by sampling points uniformly from within the simplices
of the triangulation. Each point has measure proportional to the measure
of the simplex from which it was sampled.

Introducing multiple points as representatives for the partition elements
does not introduce any bias, because in computing the
invariant measure, we use no more information than what is encoded in the
dynamics of the original data points. But from the invariant measure,
we can get a practically infinite amount of points to estimate transfer
entropy from, by generating points as described above.

After sampling, a rectangular grid is superimposed on the new point cloud
according to `binning_scheme`.
Transfer entropy is then computed using a [`BinningTransferEntropyEstimator`].
frequency estimator on those (roughly `n`) points.

#### Common use case

The invariant measure, which encodes the dynamical information, is slow to compute over
the triangulation, but only needs to be computed once. After that, transfer entropy may
be estimated at multiple scales (grids) very quickly.

This method is therefore useful if you want to explore the sensitivity of transfer entropy
to the bin size in the final rectangular grid, when you have few observations in the time series.

### Example

```julia
# Assume these points are an appropriate delay embedding {(x(t), y(t), y(t+1))} and
# that we're measure transfer entropy from x -> y.
pts = invariantize([rand(3) for i = 1:30])

v = TEVars(Tf = [3], Tpp = [2], Spp = [1])

# Compute invariant measure over a triangulation using approximate
# simplex intersections. This is relatively slow.
μ = invariantmeasure(pts, TriangulationBinning(), ApproximateIntersection())

# Compute transfer entropy from the invariant measure over multiple
# bin sizes. This is fast, because the measure has been precomputed.
tes = map(ϵ -> transferentropy(μ, v, RectangularBinning(ϵ)), 2:50)
```
"""
function transferentropy(μ::AbstractTriangulationInvariantMeasure, vars::TEVars,
        binning_scheme::RectangularBinning;
        estimator::BinningTransferEntropyEstimator = VisitationFrequency(b = 2), n::Int = 20000)

    # Get the base of the logarithm
    b = estimator.b

    dim = length(μ.points[1])
    triang = μ.triangulation.simplexindices
    n_simplices = length(triang)

    simplices = [Simplex(μ.points[triang[i]]) for i = 1:n_simplices]

    # Find a number of points to fill each simplex with so that we
    # obey the measure over the simplices of the triangulation.
    # The total number of points will be roughly `n`, but slightly
    # higher because we need integer numbers of points and use `ceil`
    # for this.
    n_fillpts_persimplex = ceil.(Int, μ.measure.dist .* n)

    # Array to store the points filling the simplices
    fillpts = Vector{SVector{dim, Float64}}()
    sizehint!(fillpts, sum(n_fillpts_persimplex))

    for i = 1:n_simplices
        sᵢ = simplices[i]

        if n_fillpts_persimplex[i] > 0
            pts = generate_interior_points(sᵢ, n_fillpts_persimplex[i])
            append!(fillpts, [SVector{dim, Float64}(pt) for pt in pts])
        end
    end

    sizehint!(fillpts, length(fillpts))

    transferentropy(fillpts, vars, binning_scheme, estimator)
end

"""
    transferentropy(pts, vars::TEVars, ϵ::RectangularBinning,
        estimator::BinningTransferEntropyEstimator) -> Float64

Compute the transfer entropy for a set of `pts` over the state space
partition specified by `ϵ` (a [`RectangularBinning`](@ref) instance).

## Fields

- **`pts`**: An ordered set of `m`-dimensional points (`pts`) representing
    an appropriate [generalised embedding]((@ref custom_delay_reconstruction))
    of some data series. Must be vector of states, not a vector of variables/time series.
    Wrap your time series in a `DynamicalSystemsBase.Dataset` first if the latter is the case.
- **`vars::TEVars`**: A [`TEVars`](@ref) instance specifying how the `m` different
    variables of `pts` are to be mapped into the marginals required for transfer
    entropy computation.
- **`ϵ::RectangularBinning`**: A [`RectangularBinning`](@ref) instance that
    dictates how the point cloud (state space reconstruction) should be discretized.
- **`estimator::BinningTransferEntropyEstimator`**. There are different ways of
    computing the transfer entropy over a discretization of a point cloud.
    The `estimator` should be a valid [`TransferEntropyEstimator`](@ref)
    that works for rectangular partitions, for example `VisitationFrequency()`
    or `TransferOperatorGrid()`. The field `estimator.b` sets the base of the
    logarithm used for the computations (e.g `VisitationFrequency(b = 2)` computes
    the transfer entropy in bits using the [VisitationFrequency](@ref) estimator).

## Returns

A single number that is the transfer entropy for the discretization of `pts`
using the binning scheme `ϵ`, using the provided `estimator`.

## Long example

```
using TransferEntropy, DynamicalSystems
```

#### 1. Generate some example time series

Let's generate some random noise series and use those as our time series.

```julia
x = rand(100)
y = rand(100)
```

We need `pts` to be a vector of states. Therefore, collect the time series in a
`DynamicalSystems.Dataset` instance. This way, the states of the composite system
will be represented as a `Vector{SVector}`.

```julia
raw_timeseries = Dataset(x, y)
```

#### 2. Generalised embedding

*Note: If your data are already organised in a form of a
[generalised embedding](@ref custom_delay_reconstruction),
where columns of the dataset correspond to lagged variables of the time series,
you can skip to step 3.*

Say we want to compute transfer entropy from ``x`` to ``y``, and that we
require a 4-dimensional embedding. For that, we need to decide on a
generalised state space reconstruction of the time series. One possible choice
is

```math
E = \\{S_{pp}, T_{pp}, T_f \\}= \\{x_t, (y_t, y_{t-\\tau}), y_{t+\\eta} \\}
```

If so, we're computing the following TE

```math
TE_{x \\to y} =  \\int_E P(x_t, y_{t-\\tau} y_t, y_{t + \\eta}) \\log{\\left( \\dfrac{P(y_{t + \\eta} | (y_t, y_{t - \\tau}, x_t)}{P(y_{t + \\eta} | y_t, y_{t-\\tau})} \\right)}.
```

To create the embedding, we'll use the [`customembed`](@ref) function (check its
documentation for a detailed explanation on how it works). This is basically just
making lagged copies of the time series, and stacking them next to each other
as column vectors. The order in which we arrange the lagged time series is not
important per se, but we need to keep track of the ordering, because that
information is crucial to the transfer entropy estimator.

According to the reconstruction we decided on above, we need to put the lagged
time series for `y` (the target variable) in the first three columns. The lags
for those columns are `η, 0, -τ`, in that order. Next, we need to put the
time series for `x` (the source variable) in the fourth column (which is not
lagged).

```julia
τ = optimal_delay(y) # find the optimal embedding lag
η = 2 # prediction lag (your choice as an analyst)
embedding_pts = customembed(raw_timeseries, Positions(2, 2, 2, 1), Lags(η, 0, -τ, 0))
```

The combination of the `Positions` instance and the `Lags` instance gives us
the necessary information about which time series in `raw_timeseries` that
corresponds to columns of the embedded dataset, and which lags each of the
columns have.

#### 3. Instructions to the estimator

The transfer entropy estimators needs the following about the columns of
the generalised reconstruction of your time series (`embedding_pts` in
our case):

- Which columns correspond to the future of the target variable (``T_f```)?
- Which columns correspond to the present and past of the target variable (``T_{pp}```)?
- Which columns correspond to the present and past of the source variable (``S_{pp}```)?
- Which columns correspond to the present/past/future of any variables
    that we are to condition on (``C_{pp}``)?

This information is needed to ensure that marginals are properly assigned during transfer
entropy computation. The estimators accept this information in the form of a `TEVars`
instance, which can be constructed like so:

```julia
vars = TEVars(Tf = [1], Tpp = [2, 3], Spp = [4])
```

#### 4. Rectangular grid specification

Entropy is essentially a property of a collection of states. To meaningfully talk about
states for our generalised state space reconstruction, we will divide the coordinate
axes of the reconstruction into a rectangular grid. Each box in the grid will
be considered a state, and the probability of visitation is equally distributed
within the box.

In this example, we'll use a rectangular partition where the box sizes are
determined by splitting each coordinate axis into 6 equally spaced
intervals, spanning the range of the data.

```julia
binning = RectangularBinning(6)
```

#### 5. Compute transfer entropy

Now we're ready to compute transfer entropy. First, let's use the
[`VisitationFrequency`](@ref) estimator with logarithm to the base 2. This
gives the transfer entropy in units of bits.

```julia
estimator = VisitationFrequency(b = 2)
te_vf = transferentropy(embedding_pts, vars, binning, estimator) #, or
```

Okay, but what if we want to use another estimator and want the transfer
entropy in units of nats? Easy.

```julia
estimator = TransferOperatorGrid(b = Base.MathConstants.e)
transferentropy(embedding_pts, vars, binning, estimator)
```

Above, we computed transfer entropy for one particular choice of partition.
Transfer entropy is a function of the  partition, and care must be taken with
the choice of partition. Below is an example where we compute transfer entropy
over 15 different cubic grids spanning the range of the data, with differing box sizes
all having fixed edge lengths  (logarithmically spaced from 0.001 to 0.3).

```julia
# Define estimator
est = VisitationFrequency(b = 2)

# Define binning schemes based on different box sizes
edgelengths = 10 .^ range(log(10, 0.001), log10(0.3), length = 15)
bs = [RectangularBinning(ϵ) for ϵ in edgelengts]

tes = map(b -> transferentropy(embedding_pts, vars, b, est), bs)
```

`tes` now contains 15 different values of the transfer entropy, one for each of
the discretization schemes. For the smallest bin sizes, the transfer entropy
is close to or equal to zero, because there are not enough points distributed
among the bins (of which there are many when the box edge length is small).
"""
function transferentropy(pts, vars::TEVars, ϵ::RectangularBinning,
    estimator::BinningTransferEntropyEstimator) end

function transferentropy(pts, vars::TEVars, ϵ::RectangularBinning,
        estimator::VisitationFrequency)

    # Base of the logarithm
    b = estimator.b

    # Collect variables for the marginals
    C = vars.conditioned_presentpast
    XY = [vars.target_future;      vars.target_presentpast; C]
    YZ = [vars.target_presentpast; vars.source_presentpast; C]
    Y =  [vars.target_presentpast;                          C]

    # Find the bins visited by the joint system (and then get
    # the marginal visits from that, so we don't have to encode
    # bins multiple times).
    joint_bin_visits = joint_visits(pts, ϵ)

    # Compute visitation frequencies for nonempty bi
    p_Y = non0hist(marginal_visits(joint_bin_visits, Y))
    p_XY = non0hist(marginal_visits(joint_bin_visits, XY))
    p_YZ = non0hist(marginal_visits(joint_bin_visits, YZ))
    p_joint = non0hist(joint_bin_visits)

    te = StatsBase.entropy(p_YZ, b) +
            StatsBase.entropy(p_XY, b) -
            StatsBase.entropy(p_Y, b) -
            StatsBase.entropy(p_joint, b)
end

function transferentropy(pts, vars::TEVars, ϵ::RectangularBinning,
        estimator::TransferOperatorGrid)

    # Base of the logarithm
    b = estimator.b

    # Collect variables for the marginals
    C = vars.conditioned_presentpast
    XY = [vars.target_future;      vars.target_presentpast; C]
    YZ = [vars.target_presentpast; vars.source_presentpast; C]
    Y =  [vars.target_presentpast;                          C]

    # Calculate the invariant distribution over the bins.
    μ = invariantmeasure(pts, ϵ)

    # Find the unique visited bins, then find the subset of those bins
    # with nonzero measure.
    positive_measure_bins = unique(μ.encoded_points)[μ.measure.nonzero_inds]

    p_Y  = marginal(Y, Array(transpose(hcat(positive_measure_bins...,))), μ.measure)
    p_XY = marginal(XY, Array(transpose(hcat(positive_measure_bins...,))), μ.measure)
    p_YZ = marginal(YZ, Array(transpose(hcat(positive_measure_bins...,))), μ.measure)

    te = StatsBase.entropy(p_YZ, b) +
            StatsBase.entropy(p_XY, b) -
            StatsBase.entropy(p_Y, b) -
            StatsBase.entropy(μ.measure.dist[μ.measure.nonzero_inds], b)
end

# Allow CustomReconstructions to be used.
transferentropy(pts::CustomReconstruction, vars::TEVars, ϵ::RectangularBinning,
    estimator::VisitationFrequency) =
    transferentropy(pts.reconstructed_pts, vars, ϵ, estimator)

transferentropy(pts::CustomReconstruction, vars::TEVars, ϵ::RectangularBinning,
    estimator::TransferOperatorGrid) =
    transferentropy(pts.reconstructed_pts, vars, ϵ, estimator)


text:
function get_dipoles(xi,yi,n,nplus,Lx,Ly)

if isempty(xi) || isempty(yi)
  @info "Input position vector is empty"
  xd = []
  yd = []
  is_dipole = []
  num_found_in_dipoles = 0
  return xd,yd,is_dipole,num_found_in_dipoles
end

κi = [ones(nplus); -ones(n-nplus)]

# x and y distances
xij,yij = distances(xi,yi)

# incorporate periodicity
periodic_distances!(xij,yij,Lx,Ly)

# inter-vortex distance matrix, set self distance to much larger than box size
rij = hypot.(xij,yij)
rij += diagm(0=>1e6*hypot(Lx,Ly)*ones(n))

# get index of minimum separations
_, nn_index = findmin(rij,dims=2)

nn_mat = ones(n)*collect(1:n)' .|> Int
nn_vec = nn_mat[nn_index][:,1]

mutuals = (collect(1:n) .== nn_vec[nn_vec])

# find dipoles
is_dipole = zeros(n)
for i in eachindex(mutuals)
is_dipole[i] = (mutuals[i] == 1 && κi[i]*κi[nn_vec[i]] == -1.0)
end

dipoles_index = findall(is_dipole .> 0.0)
dipoles_index2 = nn_vec[dipoles_index]
num_found_in_dipoles = length(dipoles_index)

xd = xi[dipoles_index]
yd = yi[dipoles_index]
xd2 = xi[dipoles_index2]
yd2 = yi[dipoles_index2]
xd = [xd[1:Int(end/2)]; xd2[1:Int(end/2)]]
yd = [yd[1:Int(end/2)]; yd2[1:Int(end/2)]]

return xd,yd,is_dipole,num_found_in_dipoles
end


text:
module DetectAlpha
export find_peak,find_peaks,alphamodel,valid_peak

using CSV
using DataFrames
using RadiationSpectra
using StatsBase
using Plots

TEST = true
DEBUG = true

include("Utils.jl")
include("Isotope.jl")
include("Spectrum.jl")
include("Models.jl")

@static if @isdefined(DEBUG) || @isdefined(TEST) TEST_DATA_FOLDER = "/home/b/Data/Spectra" end

struct Peak
    range
    channel
    function Peak()
        new(StepRange(typemin(Int32),Int32(1),typemax(Int32)),typemin(Int32))
    end
end

valid_peak(pk1::Peak) = pk1.channel != typemin(Int32) &&
    first(pk1.range) != typemin(Int32)

"""
find a peak within a certain channel range of a AlphaSpectrum
"""
function find_peak(channelrange::StepRange,as::AlphaSpectrum)
    #get the channels from the range

    # @show last(channelrange)
    # @show length(as.channels)

    if last(channelrange) <= length(as.channels) #issubset(channels,as.channels)
        println("channelrange is within alpha spectrum channel arrays")
        v = view(as.channels,channelrange)
        @show v

        return Peak()
    else
        throw(BoundsError())
    end
end

"""
detect peaks in a alpha spectrum
"""
function find_peaks(as::AlphaSpectrum,model = alphamodel,params = [256.0,15.0,4.2,10000.0])
    # println("find an alpha peak in spectrum ",as)
    #params = [256.0,15.0,4.2,10000.0]
end

end #end module DetectAlpha


text:
module JuTools

using Reexport

# include("analysis.jl")
# @reexport using JuTools.Analysis

#include("GAP.jl")
#@reexport using JuTools.GAP

include("Zm_par.jl")
@reexport using JuTools.Zm_par

include("Zm_analysis.jl")
@reexport using JuTools.Zm_analysis

end # module


text:
Base.IteratorSize(::Type{T}) where {T <: SimpleSDMLayer} = Base.HasLength()
Base.IteratorEltype(::Type{T}) where {T <: SimpleSDMLayer} = Base.HasEltype()

function Base.getindex(layer::T, i::CartesianIndex{2}) where {T <: SimpleSDMLayer}
    return layer.grid[i]
end

function Base.iterate(layer::T) where {T <: SimpleSDMLayer}
    length(layer) == 0 && return nothing
    idx = findfirst(!isnothing, layer.grid)
    return (layer[idx], LinearIndices(layer.grid)[idx])
end

function Base.iterate(layer::T, state) where {T <: SimpleSDMLayer}
    state == prod(size(layer)) && return nothing
    idx = findnext(!isnothing, layer.grid, CartesianIndices(layer.grid)[state+1])
    isnothing(idx) && return nothing
    return (layer[idx], LinearIndices(layer.grid)[idx])
end

function Base.length(layer::T) where {T <: SimpleSDMLayer}
    return count(!isnothing, layer.grid)
end


text:
# These utilities are for working with the signature_def `Dict` that comes out of
# `signature`/`splitdef`

"""
    args_tuple_expr(signature_def::Dict{Symbol})
    args_tuple_expr(arg_exprs)

For `arg_exprs` being a list of positional argument expressions from a signature, of a form
such as `[:(x::Int), :(y::Float64), :(z::Vararg)]`, or being a whole signature_def `Dict`
containing a `signature_def[:args]` value of that form.

This returns a tuple expresion containing all of the args by name. It correctly handles
splatting for things that are `Vararg` typed, e.g for the prior example `:((x, y, z...))`

This is useful for modifiying the `signature_def[:body]`.
For example, one could printout all the arguments via

```julia
signature_def[:body] = quote
    args = \$(args_tuple_expr(signature_def))
    println("args = ",args)
    \$(signature_def[:body]) # insert old body
end
```

A more realistic use case is if you want to insert a call to another function
that accepts the same arguments as the original function.
"""
function args_tuple_expr end

args_tuple_expr(signature_def::Dict{Symbol}) = args_tuple_expr(signature_def[:args])

function args_tuple_expr(arg_exprs)
    ret = Expr(:tuple)
    ret.args = map(arg_exprs) do arg

        # remove splatting (will put back on at end)
        was_splatted = Meta.isexpr(arg, :(...), 1)
        if was_splatted
            arg = arg.args[1]
        end

        # handle presence or absence of type constraints
        if Meta.isexpr(arg, :(::), 2)
            arg_name, Texpr = arg.args
        elseif arg isa Symbol
            arg_name = arg
            Texpr = nothing
        else
            error("unexpected form of argument: $arg")
        end

        # Clean up types so we can recognise if it is `Vararg`
        # remove where clauses (because they interfere with recognizing Vararg)
        if Meta.isexpr(Texpr, :where)
            Texpr = Texpr.args[1]
        end
        # remove curlies from type constraints Needs to be after removing `where`
        # important because we want to make Vararg{T,N}` into just `Vararg`
        if Meta.isexpr(Texpr, :curly)
            Texpr = Texpr.args[1]
        end
        # now can detect if should be splatted because of using Vararg in some form
        was_splatted |= Texpr == :Vararg

        #Finally apply splatting if required.
        if was_splatted
            return :($arg_name...)
        else
            return arg_name
        end
    end
    return ret
end


text:
using IOUtils
using Test
using LinearAlgebra

@testset "Output redirection" begin
    @test redirect_output(io -> print(io, "Hello")) == "Hello"
    @test redirect_output() do io
        println(io, "Hello")
        println(io, "World")
    end == "Hello\nWorld\n"
    @test redirect_output() do io
        printstyled(io, "Hello World", color=:light_red)
    end == "\e[91mHello World\e[39m"
end

@testset "Print boxed" begin
    @testset "One line" begin
        @test redirect_output() do io
            print_boxed(io, ["Hello world"])
        end == "\e[90m\e[1m[ \e[22m\e[39m Hello world\n"

        @test redirect_output() do io
            print_boxed(io, ["Hello world"], color=:light_red)
        end == "\e[91m\e[1m[ \e[22m\e[39m Hello world\n"

        @test redirect_output() do io
            print_boxed(io, ["Hello world"], "Prefix", "Suffix", color=:light_red)
        end == "\e[91m\e[1m[ \e[22m\e[39m\e[91m\e[1mPrefix \e[22m\e[39m Hello world \e[90mSuffix\e[39m\n"

        @test redirect_output() do io
            print_boxed(io, ["Hello world"], chars="(╭│╰")
        end == "\e[90m\e[1m( \e[22m\e[39m Hello world\n"
    end

    @testset "Multiple lines" begin
        @test redirect_output() do io
            print_boxed(io, ["Hello", "World"], "Prefix", "Suffix", color=:light_red)
        end == """\e[91m\e[1m┌ \e[22m\e[39m\e[91m\e[1mPrefix \e[22m\e[39m Hello
\e[91m\e[1m└ \e[22m\e[39m World \e[90mSuffix\e[39m
"""
        @test redirect_output() do io
            print_boxed(io, ["Hello", "Beautiful", "World"],
                        "Prefix", "Suffix", color=:light_red)
        end == """\e[91m\e[1m┌ \e[22m\e[39m\e[91m\e[1mPrefix \e[22m\e[39m Hello
\e[91m\e[1m│ \e[22m\e[39m Beautiful
\e[91m\e[1m└ \e[22m\e[39m World \e[90mSuffix\e[39m
"""
    end

    @testset "Block interface" begin
        @test redirect_output() do io
            print_boxed(io,
                        "Prefix", "Suffix", color=:light_red) do io
                            println(io, "Hello")
                            println(io, "Beautiful")
                            println(io, "World")
                        end
        end == """\e[91m\e[1m┌ \e[22m\e[39m\e[91m\e[1mPrefix \e[22m\e[39m Hello
\e[91m\e[1m│ \e[22m\e[39m Beautiful
\e[91m\e[1m└ \e[22m\e[39m World \e[90mSuffix\e[39m
"""

        @test redirect_output() do io
            print_boxed(io,
                        "Prefix", "Suffix", color=:light_red, chars="(╭│╰") do io
                            println(io, "Hello")
                            println(io, "Beautiful")
                            println(io, "World")
                        end
        end == """\e[91m\e[1m╭ \e[22m\e[39m\e[91m\e[1mPrefix \e[22m\e[39m Hello
\e[91m\e[1m│ \e[22m\e[39m Beautiful
\e[91m\e[1m╰ \e[22m\e[39m World \e[90mSuffix\e[39m
"""
    end
end

@testset "Indentation" begin
    @test redirect_output() do io
        indent(io, 5) do io
            println(io, "Hello")
            println(io, "World")
        end
    end == """     \e[0mHello
     \e[0mWorld
"""
    @test redirect_output() do io
        indent(io, "Important information: ") do io
            println(io, "Hello")
            println(io, "World")
        end
    end == """\e[0mImportant information: \e[0mHello
                       \e[0mWorld
"""
    @test redirect_output() do io
        indent(io, "Important information: ", color=:light_red) do io
            println(io, "Hello")
            println(io, "World")
        end
    end == """\e[91mImportant information: \e[39m\e[91mHello\e[39m
                       \e[91mWorld\e[39m
"""
    @test redirect_output() do io
        indent(io, "Important iñformation: ", color=:light_red) do io
            println(io, "Hell̃o")
            println(io, "Wôrld")
        end
    end == """\e[91mImportant iñformation: \e[39m\e[91mHell̃o\e[39m
                       \e[91mWôrld\e[39m
"""
end

@testset "@display" begin
    filename = @__FILE__
    line = @__LINE__
    @test redirect_output() do io
        @eval @display "Hello world" $io
    end == "\e[32m\e[1m[ \e[22m\e[39m\e[32m\e[1mHello world = \e[22m\e[39m \"Hello world\" \e[90m@ $(filename):$(line+2)\e[39m\n"

    line = @__LINE__
    v = identity.(1.0:10)
    @test redirect_output() do io
        @eval @display identity.(1.0:10) $io
    end == """\e[32m\e[1m┌ \e[22m\e[39m\e[32m\e[1midentity.(1.0:10) = \e[22m\e[39m 10-element $(typeof(v)):
\e[32m\e[1m│ \e[22m\e[39m   1.0
\e[32m\e[1m│ \e[22m\e[39m   2.0
\e[32m\e[1m│ \e[22m\e[39m   3.0
\e[32m\e[1m│ \e[22m\e[39m   4.0
\e[32m\e[1m│ \e[22m\e[39m   5.0
\e[32m\e[1m│ \e[22m\e[39m   6.0
\e[32m\e[1m│ \e[22m\e[39m   7.0
\e[32m\e[1m│ \e[22m\e[39m   8.0
\e[32m\e[1m│ \e[22m\e[39m   9.0
\e[32m\e[1m│ \e[22m\e[39m  10.0
\e[32m\e[1m└ \e[22m\e[39m \e[90m@ $(filename):$(line+3)\e[39m
"""
end

@testset "Horizontal lines" begin
    @test redirect_output(io -> horizontal_line(io)) ==
        "\e[90m"*repeat("━", 80)*"\e[39m\n"
    @test redirect_output(io -> horizontal_line(io,char="-")) ==
        "\e[90m"*repeat("-", 80)*"\e[39m\n"
    @test redirect_output(io -> horizontal_line(io,color=:light_red)) ==
        "\e[91m"*repeat("━", 80)*"\e[39m\n"
    @test redirect_output(io -> horizontal_line(IOContext(io, :displaysize => (10,10)))) ==
        "\e[90m"*repeat("━", 10)*"\e[39m\n"
end

@testset "Display matrices" begin
    n = 10; o = ones(n);
    T = Tridiagonal(o[2:end], -2o, o[2:end]);
    @test redirect_output() do io
        display_matrix(io, T)
    end == """ 10×10 $(typeof(T)):
⎡ -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ ⎤
⎢  1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅ ⎢
⎢   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0⎢
⎣   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0⎦
"""

    M = reshape([1], 1,1)
    @test redirect_output() do io
        display_matrix(io, M)
    end == """ 1×1 $(typeof(M)):
[ 1]
"""
end

@testset "Columns" begin
    n = 6; o = ones(n);
    T = Tridiagonal(o[2:end], -2o, o[2:end]);

    @test redirect_output() do io
        io2 = IOContext(io, :displaysize => (24,93))
        columns(io2, 1//2, 2, 1//2) do ios
            display_matrix(ios[1], T)

            println(ios[2])
            println(ios[2], "-1")

            display_matrix(ios[3], T)
        end
    end == " 6×6 $(typeof(T)):     6×6 $(typeof(T)):   \n⎡ -2.0   1.0    ⋅     ⋅     ⋅     ⋅ ⎤        -1⎡ -2.0   1.0    ⋅     ⋅     ⋅     ⋅ ⎤         \n⎢  1.0  -2.0   1.0    ⋅     ⋅     ⋅ ⎢          ⎢  1.0  -2.0   1.0    ⋅     ⋅     ⋅ ⎢         \n⎢   ⋅    1.0  -2.0   1.0    ⋅     ⋅ ⎢          ⎢   ⋅    1.0  -2.0   1.0    ⋅     ⋅ ⎢         \n⎢   ⋅     ⋅    1.0  -2.0   1.0    ⋅ ⎢          ⎢   ⋅     ⋅    1.0  -2.0   1.0    ⋅ ⎢         \n⎢   ⋅     ⋅     ⋅    1.0  -2.0   1.0⎢          ⎢   ⋅     ⋅     ⋅    1.0  -2.0   1.0⎢         \n⎣   ⋅     ⋅     ⋅     ⋅    1.0  -2.0⎦          ⎣   ⋅     ⋅     ⋅     ⋅    1.0  -2.0⎦         \n"

    @test redirect_output() do io
        io2 = IOContext(io, :displaysize => (24,93))
        columns(io2, 1//2, 2, 1//2, trim=true) do ios
            display_matrix(ios[1], T)

            println(ios[2])
            println(ios[2], "-1")

            display_matrix(ios[3], T)
        end
    end == " 6×6 $(typeof(T)):   6×6 $(typeof(T)):\n⎡ -2.0   1.0    ⋅     ⋅     ⋅     ⋅ ⎤      -1⎡ -2.0   1.0    ⋅     ⋅     ⋅     ⋅ ⎤      \n⎢  1.0  -2.0   1.0    ⋅     ⋅     ⋅ ⎢        ⎢  1.0  -2.0   1.0    ⋅     ⋅     ⋅ ⎢      \n⎢   ⋅    1.0  -2.0   1.0    ⋅     ⋅ ⎢        ⎢   ⋅    1.0  -2.0   1.0    ⋅     ⋅ ⎢      \n⎢   ⋅     ⋅    1.0  -2.0   1.0    ⋅ ⎢        ⎢   ⋅     ⋅    1.0  -2.0   1.0    ⋅ ⎢      \n⎢   ⋅     ⋅     ⋅    1.0  -2.0   1.0⎢        ⎢   ⋅     ⋅     ⋅    1.0  -2.0   1.0⎢      \n⎣   ⋅     ⋅     ⋅     ⋅    1.0  -2.0⎦        ⎣   ⋅     ⋅     ⋅     ⋅    1.0  -2.0⎦      \n"

    @test redirect_output() do io
        io2 = IOContext(io, :displaysize => (10,4))
        columns(io2, 1//2, 1//2) do ios
            println(ios[1], "ℓ̃")
            println(ios[2], "l")
            println(ios[1], "l")
        end
    end == "ℓ̃ l \nl   \n"
end


text:
module Scalings

n_pox_pairs(N) = N*(N-1)/2
export n_pox_pairs

n_pox_dir_links(N) = N*(N-1)
export n_pox_dir_links

"""
Number of links for various network size, in dense regime
"""
denseLScal_DirBin(N, C)  = C * n_pox_dir_links(N)

denseLScal_DirBin(N)  = denseLScal_DirBin(N, 0.2)
export denseLScal_DirBin

"""
Number of links for various network size, in semiDense regime
"""
semiDenseLScal_DirBin(N, C)  = C * n_pox_dir_links(N)/sqrt(N)

semiDenseLScal_DirBin(N)  = semiDenseLScal_DirBin(N, 0.1)
export semiDenseLScal_DirBin

"""
Number of links for various network size, in sparse regime
"""
sparseLScal_DirBin(N, C) =  C * N

sparseLScal_DirBin(N) = sparseLScal_DirBin(N, 4)
export sparseLScal_DirBin

"""
Average number of reciprocal pairs in Erdos Reny Model
 - avgL is the expected number of links, related to the probability p of a single link being present by expL = (N^2-N) * p

"""
erdosRenyRecScal_DirBin(expL, N) = expL^2/(2*(N^2-N))
export erdosRenyRecScal_DirBin

end