text:
import Base: show

"""
    UserNSRunner

A `UserNSRunner` represents an "execution context", an object that bundles all
necessary information to run commands within the container that contains
our crossbuild environment.  Use `run()` to actually run commands within the
`UserNSRunner`, and [`runshell()`](@ref) as a quick way to get an interactive shell
within the crossbuild environment.
"""
mutable struct UserNSRunner <: Runner
    sandbox_cmd::Cmd
    env::Dict{String, String}
    platform::AbstractPlatform

    shards::Vector{CompilerShard}
    workspace_root::String
end

function UserNSRunner(workspace_root::String;
                      cwd = nothing,
                      workspaces::Vector = Pair[],
                      mappings::Vector = Pair[],
                      platform::AbstractPlatform = HostPlatform(),
                      extra_env=Dict{String, String}(),
                      verbose::Bool = false,
                      compiler_wrapper_path::String = mktempdir(),
                      src_name::AbstractString = "",
                      shards = nothing,
                      kwargs...)
    global use_ccache, runner_override

    # Check that our kernel is new enough to use this runner
    kernel_version_check()

    # Check to make sure we're not going to try and bindmount within an
    # encrypted directory, as that triggers kernel bugs
    check_encryption(workspace_root; verbose=verbose)

    # Construct environment variables we'll use from here on out
    platform = get_concrete_platform(platform; extract_kwargs(kwargs, (:preferred_gcc_version,:preferred_llvm_version,:compilers))...)
    envs = merge(platform_envs(platform, src_name; verbose=verbose), extra_env)

    # JIT out some compiler wrappers, add it to our mounts
    generate_compiler_wrappers!(platform; bin_path=compiler_wrapper_path, extract_kwargs(kwargs, (:compilers,:allow_unsafe_flags,:lock_microarchitecture))...)
    push!(workspaces, compiler_wrapper_path => "/opt/bin")

    # the workspace_root is always a workspace, and we always mount it first
    insert!(workspaces, 1, workspace_root => "/workspace")

    # If we're enabling ccache, then mount in a read-writeable volume at /root/.ccache
    if use_ccache
        if !isdir(ccache_dir())
            mkpath(ccache_dir())
        end
        push!(workspaces, ccache_dir() => "/root/.ccache")
    end

    if isnothing(shards)
        # Choose the shards we're going to mount
        shards = choose_shards(platform; extract_kwargs(kwargs, (:preferred_gcc_version,:preferred_llvm_version,:bootstrap_list,:compilers))...)
    end

    # Construct sandbox command to look at the location it'll be mounted under
    mpath = mount_path(shards[1], workspace_root)
    sandbox_cmd = `$(mpath)/sandbox`
    if verbose
        sandbox_cmd = `$sandbox_cmd --verbose`
    end
    sandbox_cmd = `$sandbox_cmd --rootfs $(mpath)`
    if cwd != nothing
        sandbox_cmd = `$sandbox_cmd --cd $cwd`
    end

    # Add in read-only mappings and read-write workspaces
    for (outside, inside) in workspaces
        sandbox_cmd = `$sandbox_cmd --workspace $outside:$inside`
    end

    # Mount in compiler shards (excluding the rootfs shard)
    for shard in shards[2:end]
        mpath = mount_path(shard, workspace_root)
        sandbox_cmd = `$sandbox_cmd --map $(mpath):$(map_target(shard))`
    end

    # Mount in externally-defined read-only mappings
    for (outside, inside) in mappings
        sandbox_cmd = `$sandbox_cmd --map $(outside):$(inside)`
    end

    # If runner_override is not yet set, let's probe to see if we can use
    # unprivileged containers, and if we can't, switch over to privileged.
    if runner_override == ""
        if !probe_unprivileged_containers()
            msg = strip("""
            Unable to run unprivileged containers on this system!
            This may be because your kernel does not support mounting overlay
            filesystems within user namespaces. To work around this, we will
            switch to using privileged containers. This requires the use of
            sudo. To choose this automatically, set the BINARYBUILDER_RUNNER
            environment variable to "privileged" before starting Julia.
            """)
            @warn(replace(msg, "\n" => " "))
            runner_override = "privileged"
        else
            runner_override = "userns"
        end
    end

    # Check to see if we need to run privileged containers.
    if runner_override == "privileged"
        # Next, prefer `sudo`, but allow fallback to `su`. Also, force-set
        # our environmental mappings with sudo, because it is typically
        # lost and forgotten.  :(
        if sudo_cmd() == `sudo`
            sudo_envs = vcat([["-E", "$k=$(envs[k])"] for k in keys(envs)]...)
            sandbox_cmd = `$(sudo_cmd()) $(Cmd(sudo_envs)) $(sandbox_cmd)`
        else
            sandbox_cmd = `$(sudo_cmd()) "$sandbox_cmd"`
        end
    end

    # Finally, return the UserNSRunner in all its glory
    return UserNSRunner(sandbox_cmd, envs, platform, shards, workspace_root)
end

function show(io::IO, x::UserNSRunner)
    p = x.platform
    # Displays as, e.g., Linux x86_64 (glibc) UserNSRunner
    write(io, "$(typeof(p).name.name)", " ", arch(p), " ",
          Sys.islinux(p) ? "($(p.libc)) " : "",
          "UserNSRunner")
end

prompted_userns_run_privileged = false
function warn_priviledged()
    global prompted_userns_run_privileged
    if runner_override == "privileged" && !prompted_userns_run_privileged
        @info("Running privileged container via `sudo`, may ask for your password:")
        prompted_userns_run_privileged = true
    end
end

function Base.run(ur::UserNSRunner, cmd, logger::IO = stdout; verbose::Bool = false, tee_stream=stdout)
    warn_priviledged()

    did_succeed = false
    try
        mount_shards(ur; verbose=verbose)
        sandbox_cmd = `$(ur.sandbox_cmd) -- $(cmd)`
        if cmd.ignorestatus
            sandbox_cmd = ignorestatus(sandbox_cmd)
        end
        oc = OutputCollector(setenv(sandbox_cmd, ur.env); verbose=verbose, tee_stream=tee_stream)
        did_succeed = wait(oc)

        # First print out the command, then the output
        println(logger, cmd)
        print(logger, merge(oc))
    finally
        unmount_shards(ur; verbose=verbose)
    end

    # Return whether we succeeded or not
    return did_succeed
end

function Base.read(ur::UserNSRunner, cmd; verbose=false)
    warn_priviledged()

    local oc

    did_succeed = false
    try
        mount_shards(ur; verbose=verbose)
        oc = OutputCollector(setenv(`$(ur.sandbox_cmd) -- $(cmd)`, ur.env))
        did_succeed = wait(oc)
    finally
        unmount_shards(ur; verbose=verbose)
    end

    if !did_succeed
        print(stderr, collect_stderr(oc))
        return nothing
    end

    return collect_stdout(oc)
end

const AnyRedirectable = Union{Base.AbstractCmd, Base.TTY, IOStream}
function run_interactive(ur::UserNSRunner, user_cmd::Cmd; stdin = nothing, stdout = nothing, stderr = nothing, verbose::Bool = false)
    warn_priviledged()

    cmd = setenv(`$(ur.sandbox_cmd) -- $(user_cmd.exec)`, ur.env)
    if user_cmd.ignorestatus
        cmd = ignorestatus(cmd)
    end
    @debug("About to run: $(cmd)")

    if stdin isa AnyRedirectable
        cmd = pipeline(cmd, stdin=stdin)
    end
    if stdout isa AnyRedirectable
        cmd = pipeline(cmd, stdout=stdout)
    end
    if stderr isa AnyRedirectable
        cmd = pipeline(cmd, stderr=stderr)
    end

    try
        mount_shards(ur; verbose=verbose)
        if stdout isa IOBuffer
            if !(stdin isa IOBuffer)
                stdin = devnull
            end
            process = open(cmd, "r", stdin)
            @async begin
                while !eof(process)
                    write(stdout, read(process))
                end
            end
            wait(process)
            return success(process)
        else
            return success(run(cmd))
        end
    finally
        unmount_shards(ur)
    end
end

"""
    uname()

On Linux systems, return the strings returned by the `uname()` function in libc
"""
function uname()
    # Get libc and handle to uname
    libcs = filter(x -> occursin("libc.so", x), dllist())
    if isempty(libcs)
        error("Could not find libc, unable to call uname()")
    end
    libc = dlopen(first(libcs))
    uname_hdl = dlsym(libc, :uname)

    # The uname struct can have wildly differing layouts; we take advantage
    # of the fact that it is just a bunch of NULL-terminated strings laid out
    # one after the other, and that it is (as best as I can tell) at maximum
    # around 1.5KB long.  We bump up to 2KB to be safe.
    uname_struct = zeros(UInt8, 2048)
    ccall(uname_hdl, Cint, (Ptr{UInt8},), uname_struct)

    # Parse out all the strings embedded within this struct
    strings = String[]
    idx = 1
    while idx < length(uname_struct)
        # Extract string
        new_string = unsafe_string(pointer(uname_struct, idx))
        push!(strings, new_string)
        idx += length(new_string) + 1

        # Skip trailing zeros
        while uname_struct[idx] == 0 && idx < length(uname_struct)
            idx += 1
        end
    end

    return strings
end

function kernel_version_check(;verbose::Bool = false)
    # If we're not on Linux, just say everything is okay.
    if !Sys.islinux()
        return
    end

    uname_strings = try
        uname()
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end

        @warn("Unable to run `uname()` to check version number; assuming kernel version >= 3.18")
        return
    end

    # Otherwise, get the strings, convert to VersionNumber
    kernel_version = nothing

    # Some distributions tack extra stuff onto the version number.  We walk backwards
    # from the end, searching for the longest string that we can extract a VersionNumber
    # out of.  We choose a minimum length of 5, as all kernel version numbers will be at
    # least `X.Y.Z`.
    for end_idx in length(uname_strings[3]):-1:5
        try
            kernel_version = VersionNumber(uname_strings[3][1:end_idx])
            break
        catch e
            if isa(e, InterruptException)
                rethrow(e)
            end
        end
    end

    # If we were unable to parse any part of the version number, then warn and exit.
    if kernel_version === nothing
        @warn("Unable to check version number; assuming kernel version >= 3.18")
        return
    end

    # Otherwise, we have a kernel version and if it's too old, we should freak out.
    if kernel_version < v"3.18"
        error("Kernel version too old: detected $(kernel_version), need at least 3.18!")
    end

    if verbose
        @info("Parsed kernel version \"$(kernel_version)\"")
    end
end

function probe_unprivileged_containers(;verbose::Bool=false)
    # Ensure we're not about to make fools of ourselves by trying to mount an
    # encrypted directory, which triggers kernel bugs.  :(
    check_encryption(tempdir())

    function test_sandbox(; verbose=verbose, workspace_tmpdir=false, map_shard=false)
        # Choose and prepare our shards
        shards = choose_shards(Platform("x86_64", "linux"; libc="musl"))
        root_shard = first(shards)
        other_shard = last(shards)

        mktempdir() do tmpdir
            mpath = mount(root_shard, tmpdir; verbose=verbose)

            workspace_flag = []
            if workspace_tmpdir
                mkdir(joinpath(tmpdir, "workspace"))
                touch(joinpath(tmpdir, "workspace", "foo"))
                workspace_flag = `--workspace $(tmpdir)/workspace:/workspace`
            end

            map_flag = []
            if map_shard
                shard_mpath = mount(other_shard, tmpdir; verbose=verbose)
                map_flag = `--map $(shard_mpath):$(map_target(other_shard))`
            end
            try
                cmd = `$(mpath)/sandbox --rootfs $(mpath) $(workspace_flag) $(map_flag) -- /bin/sh -c "echo hello julia"`
                oc = OutputCollector(cmd; tail_error=verbose)
                return wait(oc) && merge(oc) == "hello julia\n"
            finally
                if map_shard
                    unmount(other_shard, tmpdir)
                end
                unmount(root_shard, tmpdir)
            end
        end
    end

    if verbose
        @info("Probing for unprivileged container capability...")
    end

    flags = []
    if !test_sandbox(;verbose=verbose)
        if verbose
            @warn("Unable to run simple unprivileged container test")
        end
        return false
    else
        if verbose
            @info(" * Bare-bones test passed")
        end
    end

    if !test_sandbox(;verbose=verbose, workspace_tmpdir=true)
        if verbose
            @warn("Unable to mount workspaces!")
        end
        return false
    else
        if verbose
            @info(" * Workspacing test passed")
        end
    end

    if !test_sandbox(;verbose=verbose, map_shard=true)
        if verbose
            @warn("Unable to map in a compiler shard!")
        end
        return false
    else
        if verbose
            @info(" * Shard mapping test passed")
        end
    end

    return true
end

"""
    is_ecryptfs(path::AbstractString; verbose::Bool=false)

Checks to see if the given `path` (or any parent directory) is placed upon an
`ecryptfs` mount.  This is known not to work on current kernels, see this bug
for more details: https://bugzilla.kernel.org/show_bug.cgi?id=197603

This method returns whether it is encrypted or not, and what mountpoint it used
to make that decision.
"""
function is_ecryptfs(path::AbstractString; verbose::Bool=false)
    # Canonicalize `path` immediately, and if it's a directory, add a "/" so
    # as to be consistent with the rest of this function
    path = abspath(path)
    if isdir(path)
        path = abspath(path * "/")
    end

    if verbose
        @info("Checking to see if $path is encrypted...")
    end

    # Get a listing of the current mounts.  If we can't do this, just give up
    if !isfile("/proc/mounts")
        if verbose
            @info("Couldn't open /proc/mounts, returning...")
        end
        return false, path
    end
    mounts = String(read("/proc/mounts"))

    # Grab the fstype and the mountpoints
    mounts = [split(m)[2:3] for m in split(mounts, "\n") if !isempty(m)]

    # Canonicalize mountpoints now so as to dodge symlink difficulties
    mounts = [(abspath(m[1]*"/"), m[2]) for m in mounts]

    # Fast-path asking for a mountpoint directly (e.g. not a subdirectory)
    direct_path = [m[1] == path for m in mounts]
    local parent
    if any(direct_path)
        parent = mounts[findfirst(direct_path)]
    else
        # Find the longest prefix mount:
        parent_mounts = [m for m in mounts if startswith(path, m[1])]
        if isempty(parent_mounts)
            # This is weird; this means that we can't find any mountpoints that
            # hold the given path.  I've only ever seen this in `chroot`'ed scenarios.
            return false, path
        end
        parent = parent_mounts[argmax(map(m->length(m[1]), parent_mounts))]
    end

    # Return true if this mountpoint is an ecryptfs mount
    return parent[2] == "ecryptfs", parent[1]
end

function check_encryption(workspace_root::AbstractString;
                          verbose::Bool = false)
    # If we've explicitly allowed ecryptfs, just quit out immediately
    global allow_ecryptfs
    if allow_ecryptfs
        return
    end
    msg = []

    is_encrypted, mountpoint = is_ecryptfs(workspace_root; verbose=verbose)
    if is_encrypted
        push!(msg, replace(strip("""
        Will not launch a user namespace runner within $(workspace_root), it
        has been encrypted!  Change your working directory to one outside of
        $(mountpoint) and try again.
        """), "\n" => " "))
    end

    is_encrypted, mountpoint = is_ecryptfs(storage_dir(); verbose=verbose)
    if is_encrypted
        push!(msg, replace(strip("""
        Cannot mount rootfs at $(storage_dir()), it has been encrypted!  Change
        your rootfs cache directory to one outside of $(mountpoint) by setting
        the BINARYBUILDER_ROOTFS_DIR environment variable and try again.
        """), "\n" => " "))
    end

    if !isempty(msg)
        error(join(msg, "\n"))
    end
end


text:
dThresh = 0.5 #Detect Threshold (minimum acceptable confidence)
oThresh = 0.5 #Overlap Threshold (maximum acceptable IoU)
@info "Testing all models with detectThresh = $dThresh, overlapThresh = $oThresh"

testimages = ["dog-cycle-car_nonsquare","dog-cycle-car"]
pretrained_list = [
                    YOLO.v2_tiny_416_COCO,
                    # YOLO.v2_608_COCO,
                    YOLO.v3_tiny_416_COCO,
                    YOLO.v3_320_COCO,
                    YOLO.v3_416_COCO,
                    YOLO.v3_608_COCO,
                    # YOLO.v3_spp_608_COCO
                    ]

header = ["Model" "loaded?" "load time (s)" "ran?" "run time (s)" "objects detected"]
table = Array{Any}(undef, length(pretrained_list), 6)
for (k, pretrained) in pairs(pretrained_list)
    global table
    modelname = string(pretrained)
    table[k,:] .= [modelname, false, "-", "-", "-", "-"]
    @testset "Pretrained Model: $modelname" begin
        global table

        t_load = @elapsed begin
            yolomod = pretrained(silent=true)
        end
        table[k, 2] = true
        table[k, 3] = round(t_load, digits=3)
        @info "$modelname: Loaded in $(round(t_load, digits=2)) seconds."

        batch = emptybatch(yolomod)
        for (j, imagename) in pairs(testimages)

            @info """Testing image "$imagename" """
            IMG = load(joinpath(@__DIR__,"images","$imagename.png"))
            resultsdir = joinpath(@__DIR__,"results",imagename)
            !isdir(resultsdir) && mkdir(resultsdir)
            batch[:,:,:,1], padding = prepareImage(IMG, yolomod)

            val, t_run, bytes, gctime, m = @timed res = yolomod(batch, detectThresh=dThresh, overlapThresh=oThresh);
            @test size(res,2) > 0
            table[k, 4] = true
            table[k, 5] = round(t_run, digits=4)
            table[k, 6] = size(res, 2)
            @info "$modelname: Ran in $(round(t_run, digits=2)) seconds. (bytes $bytes, gctime $gctime)"

            imgBoxes = drawBoxes(IMG, yolomod, padding, res)
            resfile = joinpath(resultsdir,"$(modelname).jpg")
            save(resfile, imgBoxes)
            @info "$modelname: View result: $resfile"

        end
    end
    GC.gc()
end
pretty_table(table, header)
@info "Times approximate. For more accurate benchmarking run ObjectDetector.benchmark()"

@testset "Custom cfg's" begin
    @testset "Valid non-square dimensions (512x384)" begin
        IMG = load(joinpath(@__DIR__,"images","dog-cycle-car.png"))
        yolomod = YOLO.v3_COCO(silent=true, cfgchanges=[(:net, 1, :width, 512), (:net, 1, :height, 384)])
        batch = emptybatch(yolomod)
        batch[:,:,:,1], padding = prepareImage(IMG, yolomod)
        res = yolomod(batch, detectThresh=dThresh, overlapThresh=oThresh) #run once
        @test size(res,2) > 0
    end
    @testset "Invalid non-square dimensions" begin
        IMG = load(joinpath(@__DIR__,"images","dog-cycle-car.png"))
        # invalid height
        @test_throws AssertionError YOLO.v3_COCO(silent=false, w=512, h=383)
        # invalid width
        @test_throws AssertionError YOLO.v3_COCO(silent=false, w=511, h=384)
    end
end

ext: jl


text:
using DifferentiableStateSpaceModels, ModelingToolkit, HMCExamples
overwrite_model_cache = true
model_cache_location = joinpath(pkgdir(HMCExamples),"src/generated_models")

include(joinpath(pkgdir(HMCExamples),"deps/rbc.jl"))

# Save off the rbc_observables_benchmark example
H, mod_vals, _ = rbc()
model_name = "rbc"
@info "RBC parameters" mod_vals.p mod_vals.p_f
save_first_order_module(H; model_name = "$(model_name)_1", model_cache_location, overwrite_model_cache, mod_vals...)
save_second_order_module(H; model_name = "$(model_name)_2", model_cache_location, overwrite_model_cache, mod_vals...)

include(joinpath(pkgdir(HMCExamples),"deps/FVGQ.jl"))

# The FVGQ20 is in this repo
H, mod_vals, _ = FVGQ20() # in this repo,
model_name = "FVGQ20"
@info "FVGQ20 parameters" mod_vals.p mod_vals.p_f
save_first_order_module(H; model_name = "$(model_name)_1", model_cache_location, overwrite_model_cache, mod_vals...)
save_second_order_module(H; model_name = "$(model_name)_2", model_cache_location, overwrite_model_cache, mod_vals...)


text:
@testset "739.daily-temperatures.jl" begin
    @test daily_temperatures([73, 74, 75, 71, 69, 72, 76, 73]) == [1, 1, 4, 2, 1, 1, 0, 0]
end


text:
mutable struct BadAgent <: AbstractAgent
    useless::Int
    pos::Int
end
mutable struct BadAgentId <: AbstractAgent
    id::Float64
end
struct ImmutableAgent <: AbstractAgent
    id::Int
end
mutable struct DiscreteVelocity <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Int}
    diameter::Float64
end
mutable struct ParametricAgent{T<:Integer} <: AbstractAgent
    id::T
    pos::NTuple{2,T}
    weight::T
    info::String
end

@testset "Model construction" begin
    # Shouldn't use ImmutableAgent since it cannot be edited
    agent = ImmutableAgent(1)
    @test_logs (
        :warn,
        "AgentType should be mutable. Try adding the `mutable` keyword infront of `struct` in your agent definition.",
    ) ABM(agent)
    # Warning is suppressed if flag is set
    @test Agents.agenttype(ABM(agent; warn = false)) <: AbstractAgent
    # Cannot use BadAgent since it has no `id` field
    @test_throws ArgumentError ABM(BadAgent)
    agent = BadAgent(1, 1)
    @test_throws ArgumentError ABM(agent)
    # Cannot use BadAgent in a grid space context since `pos` has an invalid type
    @test_throws ArgumentError ABM(BadAgent, GridSpace((1, 1)))
    @test_throws ArgumentError ABM(agent, GridSpace((1, 1)))
    # Cannot use BadAgentId since `id` has an invalid type
    @test_throws ArgumentError ABM(BadAgentId)
    agent = BadAgentId(1.0)
    @test_throws ArgumentError ABM(agent)
    # Cannot use Agent0 in a grid space context since it has no `pos` field
    @test_throws ArgumentError ABM(Agent0, GridSpace((1, 1)))
    agent = Agent0(1)
    @test_throws ArgumentError ABM(agent, GridSpace((1, 1)))
    # Cannot use Agent3 in a graph space context since `pos` has an invalid type
    @test_throws ArgumentError ABM(Agent3, GraphSpace(Agents.Graph(1)))
    agent = Agent3(1, (1, 1), 5.3)
    @test_throws ArgumentError ABM(agent, GraphSpace(Agents.Graph(1)))
    # Cannot use Agent3 in a continuous space context since `pos` has an invalid type
    @test_throws ArgumentError ABM(Agent3, ContinuousSpace((1, 1), 0.1))
    @test_throws ArgumentError ABM(agent, ContinuousSpace((1, 1), 0.1))
    # Cannot use Agent4 in a continuous space context since it has no `vel` field
    @test_throws ArgumentError ABM(Agent4, ContinuousSpace((1, 1), 0.1))
    agent = Agent4(1, (1, 1), 5)
    @test_throws ArgumentError ABM(agent, ContinuousSpace((1, 1), 0.1))
    # Shouldn't use DiscreteVelocity in a continuous space context since `vel` has an invalid type
    @test_logs (
        :warn,
        "`vel` field in Agent struct should be of type `NTuple{<:AbstractFloat}` when using ContinuousSpace.",
    ) ABM(DiscreteVelocity, ContinuousSpace((1, 1), 0.1))
    agent = DiscreteVelocity(1, (1, 1), (2, 3), 2.4)
    @test_logs (
        :warn,
        "`vel` field in Agent struct should be of type `NTuple{<:AbstractFloat}` when using ContinuousSpace.",
    ) ABM(agent, ContinuousSpace((1, 1), 0.1))
    # Warning is suppressed if flag is set
    @test Agents.agenttype(ABM(agent, ContinuousSpace((1, 1), 0.1); warn = false)) <:
          AbstractAgent
    # Shouldn't use ParametricAgent since it is not a concrete type
    @test_logs (
        :warn,
        "AgentType is not concrete. If your agent is parametrically typed, you're probably seeing this warning because you gave `Agent` instead of `Agent{Float64}` (for example) to this function. You can also create an instance of your agent and pass it to this function. If you want to use `Union` types for mixed agent models, you can silence this warning.",
    ) ABM(ParametricAgent, GridSpace((1, 1)))
    # Warning is suppressed if flag is set
    @test Agents.agenttype(ABM(ParametricAgent, GridSpace((1, 1)); warn = false)) <:
          AbstractAgent
    # ParametricAgent{Int} is the correct way to use such an agent
    @test Agents.agenttype(ABM(ParametricAgent{Int}, GridSpace((1, 1)))) <: AbstractAgent
    #Type inferance using an instance can help users here
    agent = ParametricAgent(1, (1, 1), 5, "Info")
    @test Agents.agenttype(ABM(agent, GridSpace((1, 1)))) <: AbstractAgent
    #Mixed agents
    @test Agents.agenttype(ABM(Union{Agent0,Agent1}; warn = false)) <: AbstractAgent
    @test_logs (
        :warn,
        "AgentType is not concrete. If your agent is parametrically typed, you're probably seeing this warning because you gave `Agent` instead of `Agent{Float64}` (for example) to this function. You can also create an instance of your agent and pass it to this function. If you want to use `Union` types for mixed agent models, you can silence this warning.",
    ) ABM(Union{Agent0,Agent1})
    @test_throws ArgumentError ABM(Union{Agent0,BadAgent}; warn = false)
    @test_throws ArgumentError ABM(Agent6, GridSpace((50, 50)))
    @test_throws ErrorException Agents.notimplemented(ABM(Agent0))
    # Test @agent macro
    @agent A3 GridAgent{2} begin
        weight::Float64
    end
    @test A3 <: AbstractAgent
    @test fieldnames(A3) == fieldnames(Agent3)
    @test A3.types == Agent3.types
end


text:
using PersistentHomologyTransfer
using JLDEirene
using Base.Test

@test unittest() == []


text:
"""
Extend this type for different `Potential` choices
"""
abstract type Potential end

"""
Extend this type for different `InitialData` choices
"""
abstract type InitialData end

"""
Extend this type for different `Diagnostics` operations
"""
abstract type Diagnostics end

"""
Extend this type for different `GaugeCondition`s
"""
abstract type GaugeCondition end

Base.@kwdef struct ConstantGauge <: GaugeCondition
    # order of the FD operator for solving the AH equation
    fd_order :: Int = 2
end

Base.@kwdef struct Advect_xi{T} <: GaugeCondition
    xi_vx    :: T   = 0.0
    xi_vy    :: T   = 0.0
    # order of the FD operator for solving the AH equation
    fd_order :: Int = 2
end

Base.@kwdef struct ConstantAH{T} <: GaugeCondition
    u_AH     :: T   = 1.0
    kappa    :: T   = 1.0
    # order of the FD operator for solving the xi_t and AH equations
    fd_order :: Int = 2
end

"""
Parameters for the Apparent Horizon Finder
"""
Base.@kwdef struct AHF
    itmax     :: Int      = 20
    epsilon   :: Float64  = 1e-12
end

"""
Extend this type for different `EvolutionEquations`. Needs members `ahf` and
`gaugecondition`.
"""
abstract type EvolutionEquations end

Base.@kwdef struct EvolTest0{TG<:GaugeCondition} <: EvolutionEquations
    gaugecondition :: TG  = ConstantGauge()
    ahf            :: AHF = AHF()
end

Base.@kwdef struct AffineNull{T,TP<:Potential,TG<:GaugeCondition} <: EvolutionEquations
    phi0           :: T   = 0.0
    potential      :: TP  = ZeroPotential()
    gaugecondition :: TG  = ConstantAH()
    ahf            :: AHF = AHF()
end

"""
Parameters for the time evolution
"""
Base.@kwdef struct Integration{T,Tdt,S}
    dt              :: Tdt  = :auto
    tmax            :: T
    ODE_method      :: S    = AB3()
    adaptive        :: Bool = false
    # relative tolerance for adaptive integrators
    reltol          :: Float64 = 1e-6
    filter_poststep :: Bool = true
end

"""
Parameters for Input/Output
"""
Base.@kwdef struct InOut
    # negative values suppress output
    out_boundary_every          :: Int  = -1
    out_bulk_every              :: Int  = -1
    out_gauge_every             :: Int  = -1
    out_bulkconstrained_every   :: Int  = -1

    out_boundary_every_t        :: Float64  = -1.0
    out_bulk_every_t            :: Float64  = -1.0
    out_gauge_every_t           :: Float64  = -1.0
    out_bulkconstrained_every_t :: Float64  = -1.0

    checkpoint_every_walltime_hours :: Float64 = -1.0

    # stop and checkpoint upon reaching this walltime
    max_walltime       :: Float64 = 1.e20

    # trigger termination
    termination_from_file :: Bool    = true
    check_file_every      :: Int     = 10
    termination_file      :: String  = "TERMINATE"

    # name of script
    _parfile           :: String  = splitext(basename(Base.source_path()))[1]

    # use name of script by default
    out_dir            :: String  = _parfile
    checkpoint_dir     :: String  = _parfile

    recover            :: Symbol  = :auto
    recover_dir        :: String  = _parfile

    # be very careful with this option! it will remove the whole folder contents
    # if set to true! use only for fast debugging runs
    remove_existing    :: Bool  = false
end

abstract type AbstractVars{T} <: AbstractVector{T} end

struct BulkEvolved{T} <: AbstractVars{T}
    B1  :: Array{T,3}
    B2  :: Array{T,3}
    G   :: Array{T,3}
    phi :: Array{T,3}
end

struct BulkConstrained{T} <: AbstractVars{T}
    S    :: Array{T,3}
    Fx   :: Array{T,3}
    Fy   :: Array{T,3}
    B1d  :: Array{T,3}
    B2d  :: Array{T,3}
    Gd   :: Array{T,3}
    phid :: Array{T,3}
    Sd   :: Array{T,3}
    A    :: Array{T,3}
end

struct Bulk{T} <: AbstractVars{T}
    B1   :: Array{T,3}
    B2   :: Array{T,3}
    G    :: Array{T,3}
    phi  :: Array{T,3}
    S    :: Array{T,3}
    Fx   :: Array{T,3}
    Fy   :: Array{T,3}
    B1d  :: Array{T,3}
    B2d  :: Array{T,3}
    Gd   :: Array{T,3}
    phid :: Array{T,3}
    Sd   :: Array{T,3}
    A    :: Array{T,3}
end

struct BulkDeriv{T}
    Du_B1   :: Array{T,3}
    Du_B2   :: Array{T,3}
    Du_G    :: Array{T,3}
    Du_phi  :: Array{T,3}
    Du_S    :: Array{T,3}
    Du_Fx   :: Array{T,3}
    Du_Fy   :: Array{T,3}
    Du_Sd   :: Array{T,3}
    Du_B1d  :: Array{T,3}
    Du_B2d  :: Array{T,3}
    Du_Gd   :: Array{T,3}
    Du_A    :: Array{T,3}
    Duu_B1  :: Array{T,3}
    Duu_B2  :: Array{T,3}
    Duu_G   :: Array{T,3}
    Duu_phi :: Array{T,3}
    Duu_S   :: Array{T,3}
    Duu_Fx  :: Array{T,3}
    Duu_Fy  :: Array{T,3}
    Duu_A   :: Array{T,3}
end

struct Boundary{T} <: AbstractVars{T}
    a4  :: Array{T,3}
    fx2 :: Array{T,3}
    fy2 :: Array{T,3}
end

struct Gauge{T} <: AbstractVars{T}
    xi  :: Array{T,3}
end

@inline varlist(::BulkEvolved)     = (:B1, :B2, :G, :phi)
@inline varlist(::BulkConstrained) = (:S, :Fx, :Fy, :B1d, :B2d, :Gd, :phid, :Sd, :A)
@inline varlist(::Bulk)            = (:B1, :B2, :G, :phi, :S, :Fx, :Fy, :B1d, :B2d,
                                      :Gd, :phid, :Sd, :A)
@inline varlist(::Boundary)        = (:a4, :fx2, :fy2)
@inline varlist(::Gauge)           = (:xi,)

"""
    BulkEvolved{T}(undef, Nu, Nx, Ny)

Construct a container of uninitialized Arrays to hold all the bulk variables
that are evolved in time: B1, B2, G, phi
"""
function BulkEvolved{T}(::UndefInitializer, Nu::Int, Nx::Int, Ny::Int) where {T<:Real}
    B1  = Array{T}(undef, Nu, Nx, Ny)
    B2  = Array{T}(undef, Nu, Nx, Ny)
    G   = Array{T}(undef, Nu, Nx, Ny)
    phi = Array{T}(undef, Nu, Nx, Ny)
    BulkEvolved{T}(B1, B2, G, phi)
end

"""
    BulkConstrained{T}(undef, Nu, Nx, Ny)

Construct a container of uninitialized Arrays to hold all the bulk variables
that are constrained (not evolved in time): S, Fx, Fy, B1d, B2d, Gd, phid, Sd, A
"""
function BulkConstrained{T}(::UndefInitializer, Nu::Int, Nx::Int, Ny::Int) where {T<:Real}
    S    = Array{T}(undef, Nu, Nx, Ny)
    Fx   = Array{T}(undef, Nu, Nx, Ny)
    Fy   = Array{T}(undef, Nu, Nx, Ny)
    B1d  = Array{T}(undef, Nu, Nx, Ny)
    B2d  = Array{T}(undef, Nu, Nx, Ny)
    Gd   = Array{T}(undef, Nu, Nx, Ny)
    phid = Array{T}(undef, Nu, Nx, Ny)
    Sd   = Array{T}(undef, Nu, Nx, Ny)
    A    = Array{T}(undef, Nu, Nx, Ny)
    BulkConstrained{T}(S, Fx, Fy, B1d, B2d, Gd, phid, Sd, A)
end

"""
    Bulk(bulkevol::BulkEvolved, bulkconstrain::BulkConstrained)

Construct a container to hold all the bulk variables where the evolved variables
point to the given bulkevol struct and the constrained ones point to the bulkconstrain struct
"""
function Bulk(bulkevol::BulkEvolved{T}, bulkconstrain::BulkConstrained{T}) where {T}
    B1    = bulkevol.B1
    B2    = bulkevol.B2
    G     = bulkevol.G
    phi   = bulkevol.phi
    S     = bulkconstrain.S
    Fx    = bulkconstrain.Fx
    Fy    = bulkconstrain.Fy
    B1d   = bulkconstrain.B1d
    B2d   = bulkconstrain.B2d
    Gd    = bulkconstrain.Gd
    phid  = bulkconstrain.phid
    Sd    = bulkconstrain.Sd
    A     = bulkconstrain.A
    Bulk{T}(B1, B2, G, phi, S, Fx, Fy, B1d, B2d, Gd, phid, Sd, A)
end

function BulkDeriv{T}(::UndefInitializer, Nu::Int, Nx::Int, Ny::Int) where {T<:Real}
    Du_B1    = Array{T}(undef, Nu, Nx, Ny)
    Du_B2    = Array{T}(undef, Nu, Nx, Ny)
    Du_G     = Array{T}(undef, Nu, Nx, Ny)
    Du_phi   = Array{T}(undef, Nu, Nx, Ny)
    Du_S     = Array{T}(undef, Nu, Nx, Ny)
    Du_Fx    = Array{T}(undef, Nu, Nx, Ny)
    Du_Fy    = Array{T}(undef, Nu, Nx, Ny)
    Du_Sd    = Array{T}(undef, Nu, Nx, Ny)
    Du_B1d   = Array{T}(undef, Nu, Nx, Ny)
    Du_B2d   = Array{T}(undef, Nu, Nx, Ny)
    Du_Gd    = Array{T}(undef, Nu, Nx, Ny)
    Du_A     = Array{T}(undef, Nu, Nx, Ny)
    Duu_B1   = Array{T}(undef, Nu, Nx, Ny)
    Duu_B2   = Array{T}(undef, Nu, Nx, Ny)
    Duu_G    = Array{T}(undef, Nu, Nx, Ny)
    Duu_phi  = Array{T}(undef, Nu, Nx, Ny)
    Duu_S    = Array{T}(undef, Nu, Nx, Ny)
    Duu_Fx   = Array{T}(undef, Nu, Nx, Ny)
    Duu_Fy   = Array{T}(undef, Nu, Nx, Ny)
    Duu_A    = Array{T}(undef, Nu, Nx, Ny)
    BulkDeriv{T}(Du_B1, Du_B2, Du_G, Du_phi, Du_S, Du_Fx, Du_Fy, Du_Sd, Du_B1d,
                 Du_B2d, Du_Gd, Du_A,
                 Duu_B1, Duu_B2, Duu_G, Duu_phi, Duu_S, Duu_Fx, Duu_Fy, Duu_A)
end

"""
    Boundary{T}(undef, Nx, Ny)

Construct a container of uninitialized Arrays to hold all the boundary
variables: a4, fx2, fy2

These variables are automatically defined on a `(1,Nx,Ny)` grid, rather than
a `(Nx,Ny)` one, so that the same Dx and Dy differential operators defined for
the bulk quantities can also straightforwardly apply on them. Remember that the
axis along which the operator applies is defined on the operator itself. So, by
defining things this way, the Dx operator (which acts along the 2nd index) will
also do the correct thing when acting on a4, fx2, fy2.
"""
function Boundary{T}(::UndefInitializer, Nx::Int, Ny::Int) where {T<:Real}
    a4  = Array{T}(undef, 1, Nx, Ny)
    fx2 = Array{T}(undef, 1, Nx, Ny)
    fy2 = Array{T}(undef, 1, Nx, Ny)
    Boundary{T}(a4, fx2, fy2)
end

"""
    Gauge{T}(undef, Nx, Ny)

Construct a container with an uninitialized Array to hold the gauge variable xi

As in the Boundary struct, xi is automatically defined on a `(1,Nx,Ny)` grid,
rather than a `(Nx,Ny)` one, so that the same Dx and Dy differential operators
defined for the bulk quantities can also straightforwardly apply on it.
"""
function Gauge{T}(::UndefInitializer, Nx::Int, Ny::Int) where {T<:Real}
    xi  = Array{T}(undef, 1, Nx, Ny)
    Gauge{T}(xi)
end

@inline getB1(ff::BulkEvolved)       = ff.B1
@inline getB2(ff::BulkEvolved)       = ff.B2
@inline getG(ff::BulkEvolved)        = ff.G
@inline getphi(ff::BulkEvolved)      = ff.phi

@inline getS(ff::BulkConstrained)    = ff.S
@inline getFx(ff::BulkConstrained)   = ff.Fx
@inline getFy(ff::BulkConstrained)   = ff.Fy
@inline getB1d(ff::BulkConstrained)  = ff.B1d
@inline getB2d(ff::BulkConstrained)  = ff.B2d
@inline getGd(ff::BulkConstrained)   = ff.Gd
@inline getphid(ff::BulkConstrained) = ff.phid
@inline getSd(ff::BulkConstrained)   = ff.Sd
@inline getA(ff::BulkConstrained)    = ff.A

@inline getB1(ff::Bulk)              = ff.B1
@inline getB2(ff::Bulk)              = ff.B2
@inline getG(ff::Bulk)               = ff.G
@inline getphi(ff::Bulk)             = ff.phi
@inline getS(ff::Bulk)               = ff.S
@inline getFx(ff::Bulk)              = ff.Fx
@inline getFy(ff::Bulk)              = ff.Fy
@inline getB1d(ff::Bulk)             = ff.B1d
@inline getB2d(ff::Bulk)             = ff.B2d
@inline getGd(ff::Bulk)              = ff.Gd
@inline getphid(ff::Bulk)            = ff.phid
@inline getSd(ff::Bulk)              = ff.Sd
@inline getA(ff::Bulk)               = ff.A

@inline geta4(ff::Boundary)          = ff.a4
@inline getfx2(ff::Boundary)         = ff.fx2
@inline getfy2(ff::Boundary)         = ff.fy2

@inline getxi(ff::Gauge)             = ff.xi

@inline function Base.length(ff::AbstractVars)
    vars = varlist(ff)
    sum_l = 0
    for x in vars
        f = getproperty(ff,x)   # f will point to each variable in the ff struct
        sum_l += length(f)
    end
    sum_l
end
@inline Base.size(ff::AbstractVars) = (length(ff),)

# indexing. this is just a linear indexing through all the arrays. adapted from
# RecursiveArrayTools
@inline Base.firstindex(ff::AbstractVars) = 1
@inline Base.lastindex(ff::AbstractVars) = length(ff)

@inline function Base.getindex(ff::AbstractVars, i::Int)
    vars = varlist(ff)
    @inbounds for x in vars
        f  = getproperty(ff,x)
        i -= length(f)
        if i <= 0
            return f[length(f)+i]
        end
    end
end
@inline function Base.setindex!(ff::AbstractVars, v, i::Int)
    vars = varlist(ff)
    @inbounds for x in vars
        f  = getproperty(ff,x)
        i -= length(f)
        if i <= 0
            f[length(f)+i] = v
            break
        end
    end
end

Base.similar(ff::BulkEvolved{T}) where{T} =
    BulkEvolved{T}(similar(ff.B1), similar(ff.B2), similar(ff.G), similar(ff.phi))

Base.similar(ff::Boundary{T}) where{T} =
    Boundary{T}(similar(ff.a4), similar(ff.fx2), similar(ff.fy2))

Base.similar(ff::Gauge{T}) where{T} = Gauge{T}(similar(ff.xi))

"""
    unpack(ff::AbstractVars)

Return `Array` with all the different fields in the structure
"""
function unpack(ff::AbstractVars)
    vars = varlist(ff)
    [getproperty(ff,x) for x in vars]
end

"""
    BulkPartition{N,A} <: AbstractPartition{N,A}

Container to store (bulk) quantities that may be spread across different grid
partitions. This is to be thought as a `Tuple` (or an `Array`) of `Bulk` objects.
"""
struct BulkPartition{N,A} <: AbstractPartition{N,A}
    x :: NTuple{N,A}
end
BulkPartition(x...) = BulkPartition(tuple(x...))

function BulkEvolved(bulks::BulkPartition{N}) where{N}
    f = ntuple(i -> BulkEvolved(bulks[i]), N)
    BulkPartition(f)
end

function BulkConstrained(bulks::BulkPartition{N}) where{N}
    f = ntuple(i -> BulkConstrained(bulks[i]), N)
    BulkPartition(f)
end

struct EvolVars{T,N,A} <: FlattenedVector{T,N,A}
    x :: NTuple{N,A}
end
EvolVars(x::NTuple{N,A}) where {N,A} = EvolVars{eltype(A),N,A}(x)
EvolVars(x...) = EvolVars(tuple(x...))

function EvolVars(ff::AbstractVector{A}) where{A<:AbstractArray}
    x = Tuple(ff)
    EvolVars{eltype(A),length(x),A}(x)
end

"""
    EvolVars(boundary::Boundary, gauge::Gauge, bulkevols::NTuple)

Build a container to store all the evolved quantities as elements of an
`NTuple`. The idea is to treat them as a single column vector for the point of
view of the time evolution routine.
"""
function EvolVars(boundary::Boundary, gauge::Gauge,
                  bulkevols::BulkPartition{Nsys}) where {Nsys}
    f1 = unpack(boundary)
    f2 = unpack(gauge)

    # this technique allocates a bit of memory and is not fully type-stable, but
    # i think it's not a problem since we only call this function once
    f3 = [unpack(bulkevol) for bulkevol in bulkevols]
    EvolVars([f1; f2; f3...])
end

Base.similar(ff::EvolVars{T,N,S}) where {T,N,S} = EvolVars{T,N,S}(similar.(ff.x))

@inline getudomains(::EvolVars{T,N}) where {T,N} = div(N-4, 4)

@inline geta4(ff::EvolVars)   = ff.x[1]
@inline getfx2(ff::EvolVars)  = ff.x[2]
@inline getfy2(ff::EvolVars)  = ff.x[3]
@inline getxi(ff::EvolVars)   = ff.x[4]

function getB1(ff::EvolVars, i::Int)
    Nsys = getudomains(ff)
    @assert i > 0
    @assert i <= Nsys
    ff.x[5 + (i-1)*4]
end

function getB2(ff::EvolVars, i::Int)
    Nsys = getudomains(ff)
    @assert i > 0
    @assert i <= Nsys
    ff.x[6 + (i-1)*4]
end

function getG(ff::EvolVars, i::Int)
    Nsys = getudomains(ff)
    @assert i > 0
    @assert i <= Nsys
    ff.x[7 + (i-1)*4]
end

function getphi(ff::EvolVars, i::Int)
    Nsys = getudomains(ff)
    @assert i > 0
    @assert i <= Nsys
    ff.x[8 + (i-1)*4]
end

@inline getboundary(ff::EvolVars) = Boundary(geta4(ff), getfx2(ff), getfy2(ff))
@inline getgauge(ff::EvolVars)    = Gauge(getxi(ff))

@inline getbulkevolved(ff::EvolVars, i::Int) =
    BulkEvolved(getB1(ff,i), getB2(ff,i), getG(ff,i), getphi(ff,i))

function getbulkevolvedpartition(ff::EvolVars)
    Nsys = getudomains(ff)
    f = ntuple(i -> getbulkevolved(ff,i), Nsys)
    BulkPartition(f)
end

struct BulkHorizon{T}
    B1_uAH      :: Array{T,3}
    B2_uAH      :: Array{T,3}
    G_uAH       :: Array{T,3}
    phi_uAH     :: Array{T,3}
    S_uAH       :: Array{T,3}
    Fx_uAH      :: Array{T,3}
    Fy_uAH      :: Array{T,3}
    Sd_uAH      :: Array{T,3}
    B1d_uAH     :: Array{T,3}
    B2d_uAH     :: Array{T,3}
    Gd_uAH      :: Array{T,3}
    phid_uAH    :: Array{T,3}
    A_uAH       :: Array{T,3}

    Du_B1_uAH   :: Array{T,3}
    Du_B2_uAH   :: Array{T,3}
    Du_G_uAH    :: Array{T,3}
    Du_phi_uAH  :: Array{T,3}
    Du_S_uAH    :: Array{T,3}
    Du_Fx_uAH   :: Array{T,3}
    Du_Fy_uAH   :: Array{T,3}
    Du_Sd_uAH   :: Array{T,3}
    Du_B1d_uAH  :: Array{T,3}
    Du_B2d_uAH  :: Array{T,3}
    Du_Gd_uAH   :: Array{T,3}
    Du_A_uAH    :: Array{T,3}

    Duu_B1_uAH  :: Array{T,3}
    Duu_B2_uAH  :: Array{T,3}
    Duu_G_uAH   :: Array{T,3}
    Duu_S_uAH   :: Array{T,3}
    Duu_Fx_uAH  :: Array{T,3}
    Duu_Fy_uAH  :: Array{T,3}
    Duu_A_uAH   :: Array{T,3}
end
function BulkHorizon{T}(Nx::Int, Ny::Int) where {T<:Real}
    B1_uAH      = Array{T}(undef, 1, Nx, Ny)
    B2_uAH      = Array{T}(undef, 1, Nx, Ny)
    G_uAH       = Array{T}(undef, 1, Nx, Ny)
    phi_uAH     = Array{T}(undef, 1, Nx, Ny)
    S_uAH       = Array{T}(undef, 1, Nx, Ny)
    Fx_uAH      = Array{T}(undef, 1, Nx, Ny)
    Fy_uAH      = Array{T}(undef, 1, Nx, Ny)
    Sd_uAH      = Array{T}(undef, 1, Nx, Ny)
    B1d_uAH     = Array{T}(undef, 1, Nx, Ny)
    B2d_uAH     = Array{T}(undef, 1, Nx, Ny)
    Gd_uAH      = Array{T}(undef, 1, Nx, Ny)
    phid_uAH    = Array{T}(undef, 1, Nx, Ny)
    A_uAH       = Array{T}(undef, 1, Nx, Ny)

    Du_B1_uAH   = Array{T}(undef, 1, Nx, Ny)
    Du_B2_uAH   = Array{T}(undef, 1, Nx, Ny)
    Du_G_uAH    = Array{T}(undef, 1, Nx, Ny)
    Du_phi_uAH  = Array{T}(undef, 1, Nx, Ny)
    Du_S_uAH    = Array{T}(undef, 1, Nx, Ny)
    Du_Fx_uAH   = Array{T}(undef, 1, Nx, Ny)
    Du_Fy_uAH   = Array{T}(undef, 1, Nx, Ny)
    Du_Sd_uAH   = Array{T}(undef, 1, Nx, Ny)
    Du_B1d_uAH  = Array{T}(undef, 1, Nx, Ny)
    Du_B2d_uAH  = Array{T}(undef, 1, Nx, Ny)
    Du_Gd_uAH   = Array{T}(undef, 1, Nx, Ny)
    Du_A_uAH    = Array{T}(undef, 1, Nx, Ny)

    Duu_B1_uAH  = Array{T}(undef, 1, Nx, Ny)
    Duu_B2_uAH  = Array{T}(undef, 1, Nx, Ny)
    Duu_G_uAH   = Array{T}(undef, 1, Nx, Ny)
    Duu_S_uAH   = Array{T}(undef, 1, Nx, Ny)
    Duu_Fx_uAH  = Array{T}(undef, 1, Nx, Ny)
    Duu_Fy_uAH  = Array{T}(undef, 1, Nx, Ny)
    Duu_A_uAH   = Array{T}(undef, 1, Nx, Ny)

    BulkHorizon{T}(B1_uAH, B2_uAH, G_uAH, phi_uAH, S_uAH, Fx_uAH, Fy_uAH,
                   Sd_uAH, B1d_uAH, B2d_uAH, Gd_uAH, phid_uAH, A_uAH, Du_B1_uAH,
                   Du_B2_uAH, Du_G_uAH, Du_phi_uAH, Du_S_uAH, Du_Fx_uAH,
                   Du_Fy_uAH, Du_Sd_uAH, Du_B1d_uAH, Du_B2d_uAH, Du_Gd_uAH,
                   Du_A_uAH, Duu_B1_uAH, Duu_B2_uAH, Duu_G_uAH, Duu_S_uAH,
                   Duu_Fx_uAH, Duu_Fy_uAH, Duu_A_uAH)
end

struct HorizonCache{T1,T2,D}
    bulkhorizon :: BulkHorizon{T1}
    axx         :: Vector{T2}
    ayy         :: Vector{T2}
    axy         :: Vector{T2}
    bx          :: Vector{T2}
    by          :: Vector{T2}
    cc          :: Vector{T2}
    b_vec       :: Vector{T2}
    Dx_2D       :: D
    Dy_2D       :: D
    Dxx_2D      :: D
    Dyy_2D      :: D
    Dxy_2D      :: D
    _Dx_2D      :: D
    _Dy_2D      :: D
    _Dxx_2D     :: D
    _Dyy_2D     :: D
    _Dxy_2D     :: D
end


text:
"""
    Wrapper to SlepcInitializeNoPointers

# Implementation
I don't know if I am supposed to use PetscInt or not...
"""
function SlepcInitialize(args::Vector{String}, filename::String, help::String)
    args2 = ["julia"; args]
    nargs = Cint(length(args2))
    error = ccall( (:SlepcInitializeNoPointers, libslepc),
            PetscErrorCode,
            (Cint,
            Ptr{Ptr{UInt8}},
            Cstring,
            Cstring),
            nargs, args2, filename, help
    )
    @assert iszero(error)
end

SlepcInitialize(args::Vector{String}) = SlepcInitialize(args, "", "")
SlepcInitialize(args::String) = SlepcInitialize(convert(Vector{String}, split(args)), "", "")

"""
    Initialize SLEPc.

If `cmd_line_args == true`, then command line arguments passed to Julia are used as
arguments for SLEPc (leading to a call to `SlepcInitializeNoPointers`).

Otherwise, if `cmd_line_args == false`, initialize SLEPc without arguments (leading
to a call to `SlepcInitializeNoArguments`).
"""
function SlepcInitialize(cmd_line_args::Bool = false)
    if (cmd_line_args)
        SlepcInitialize(ARGS)
    else
        error = ccall((:SlepcInitializeNoArguments, libslepc), PetscErrorCode, ())
        @assert iszero(error)
    end
end

"""
    Wrapper to SlepcFinalize
"""
function SlepcFinalize()
    error = ccall( (:SlepcFinalize, libslepc), PetscErrorCode, ())
    @assert iszero(error)
end


text:
# the main visualization module for OpenGL abstraction rendering
module vis

using GLFW
using ModernGL
using Serialization

struct Context
  window::GLFW.Window
  program::GLuint
  width::Int
  height::Int
  attributes::Array{GLint}
  Tmat_location::GLint
end

include("vis_util.jl")
include("vis_drawing.jl")
include("vis_scene.jl")

function make_shader_program(vsh::String, fsh::String)
  createcontextinfo()
  version_str = get_glsl_version_string()

  vertexShader = createShader(version_str * vsh, GL_VERTEX_SHADER)
  fragmentShader = createShader(version_str * fsh, GL_FRAGMENT_SHADER)
  program = createShaderProgram(vertexShader, fragmentShader)

  return program
end

function make_context(width, height)
  # Create a window and its OpenGL context
  GLFW.WindowHint(GLFW.SAMPLES, 4)
  GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
  GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 3)
  GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, Cint(1)) # true, required by OS X
  GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
  window = GLFW.CreateWindow(width, height, "window.jl")

  # Make the window's context current
  GLFW.MakeContextCurrent(window)

  # Compile the Shader Program
  dir_path = @__DIR__
  vsh = read(dir_path * "/shaders/shader.vert", String)
  fsh = read(dir_path * "/shaders/shader.frag", String)
  program = make_shader_program(vsh, fsh)
  glUseProgram(program)

  positionAttribute = glGetAttribLocation(program, "position")
  colorAttribute = glGetAttribLocation(program, "color")
  usetexAttribute = glGetAttribLocation(program, "usetex")
  texcoordAttribute = glGetAttribLocation(program, "texcoord")
  attributes = [positionAttribute, colorAttribute, usetexAttribute,
                texcoordAttribute]
  Tmat_location = glGetUniformLocation(program, "Tmat")

  return Context(window, program, width, height, attributes, Tmat_location)
end

function print_error_if_not_empty()
  err = glErrorMessage()
  if err != ""
    println(err)
  end
end

function setup()
  # Load the Window -----------------------------------------------------------
  width = 1080
  height = round(Int, 1080 * 9 / 16)
  context = make_context(width, height)
  # ---------------------------------------------------------------------------

  # window scaling matrix -----------------------------------------------------
  W = scale_mat(context.height/context.width, 1)
  glUniformMatrix4fv(glGetUniformLocation(context.program, "Wmat"), 1,
                     GL_FALSE, W)
  # ---------------------------------------------------------------------------

  # Load Font Texture ---------------------------------------------------------
  textureBuffer = glGenTexture()
  glActiveTexture(GL_TEXTURE0)

  dir_path = @__DIR__
  fp = open(dir_path * "/font/font.bin", "r")
  pixels = deserialize(fp)
  close(fp)
  (w, h) = size(pixels)
  w = div(w, 4)

  glEnable(GL_BLEND)
  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)
  glBindTexture(GL_TEXTURE_2D, textureBuffer)
  glBindTexture(GL_TEXTURE_2D, textureBuffer)
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, w, h, 0, GL_RGBA, GL_FLOAT, pixels);

  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);

  glUniform1i(glGetUniformLocation(context.program, "tex"), 0)
  # ---------------------------------------------------------------------------

  return context
end

function test()
  context = setup()

  # Make a Bunch of Test Objects ----------------------------------------------
  # triangle
  position_data = GLfloat[-1.0, 1.0,
                          -1.0, -1.0,
                          -0.1, 0.0]
  position = RenderData(position_data, 2, GL_STATIC_DRAW)

  color_data = GLfloat[0.0, 0.0, 1.0,
                       1.0, 0.0, 0.0,
                       0.0, 1.0, 0.0]
  color = RenderData(color_data, 3, GL_STATIC_DRAW)
  usetex = RenderData(fill(GLfloat(0), 3), 1, GL_STATIC_DRAW)
  texcoord = RenderData(fill(GLfloat(0), 6), 2, GL_STATIC_DRAW)
  idx = GLuint[0, 1, 2]
  object2 = RenderObject(context, [position, color, usetex, texcoord], idx)

  # road
  th = range(0, stop=(2 * pi), length=30)
  r = 0.6
  x = fill(0.0, length(th))
  y = fill(0.0, length(th))
  for i in 1:length(th)
    x[i] = r * cos(th[i])
    y[i] = r * sin(th[i])
  end
  object3 = make_road(context, x, y, 0.1)

  # font square
  position_data = GLfloat[-1.0, 0.5,
                          -1.0, 1.0,
                          0.0, 1.0,
                          0.0, 0.5]
  position = RenderData(position_data, 2, GL_STATIC_DRAW)
  color = RenderData(repeat(GLfloat[1.0, 0.0, 0.0], 4), 3, GL_STATIC_DRAW)
  usetex = RenderData(fill(GLfloat(1.0), 4), 1, GL_STATIC_DRAW)
  texcoord_data = GLfloat[0.0, 1.0,
                          0.0, 0.0,
                          1.0, 0.0,
                          1.0, 1.0]
  texcoord = RenderData(texcoord_data, 2, GL_STATIC_DRAW)
  idx = GLuint[0, 1, 2,
               0, 2, 3]
  object4 = RenderObject(context, [position, color, usetex, texcoord], idx)

  # line
  color = [0.0, 0.0, 0.0]
  l1 = make_line(context, 0.0, 0.0, 0.5, -0.5, color)
  s1 = make_text(context, string(time_ns()), 0.0, 0.0)

  # vector
  vec1 = make_vector_xy(context, 0.2, -0.2, 0.5, 0.5)

  # car
  car1 = make_car(context)
  car1.T = scale_mat(0.01, 0.01)

  # Main Render Loop ----------------------------------------------------------
  k = 0
  text = ""

  println("Set run to false")
  t0 = time_ns() / 1e9
  window_status = true

  k = 0
  t = 0
  while window_status
    #println(1e3 * (time_ns() / 1e9 - t))
    t = time_ns() / 1e9
    points = GLfloat[0.0, 0.0, cos(t), sin(t)]
    update_buffer!(l1, points, context.attributes[1])
    if k == 60
      #text = string(rand(1:100))
      #update_text!(s1, text, 0.0, 0.0)
      car_lights!(car1)
      k = 0
    end
    update_text!(s1, string(time_ns()), 0.0, 0.0)
    car1.T = scale_mat(0.01, 0.01) * translate_mat(0, (t - t0) * 0.1)

    k += 1

    window_status = visualize(context, [object2, object3, object4, l1,
                                        s1, car1, vec1])
  end

  if window_status == true
    GLFW.DestroyWindow(context.window)
  end

  return
end

function visualize(context::Context, objects::Array{Any, 1})
  if GLFW.WindowShouldClose(context.window)
    GLFW.DestroyWindow(context.window)

    return false
  end

  glClearColor(0.3, 0.7, 0.0, 1.0)
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

  for obj in objects
    render(obj)
  end

  GLFW.SwapBuffers(context.window)
  GLFW.PollEvents()

  return true
end

#test()

end


text:
module LHEF

using EzXML

export parse_lhe, flatparticles

struct Event
    header
    particles
end

function Base.show(io::IO, evt::Event)
    pads = (3, 6, 6, 7, 7, 6, 6, 14, 14, 14, 14, 6, 8, 6)

    println(io, "  Event header: ", evt.header)
    parts = evt.particles
    println(io, "  Event particles:")
    print(io, "  ")
    ks = lpad.(keys(first(parts)), pads)
    print(io, join(ks, "| ")...)
    println(io)

    for p in parts
        ps = lpad.(values(p), pads)
        print(io, "  ")
        print(io, join(ps, ", ")...)
        println(io)
    end
end

function parse_event(event)
    lines = split(event, '\n'; keepempty=false)
    headerdata = split(lines[1], ' '; keepempty=false)
    header = (;
        nparticles=parse(Int16, headerdata[1]), # Number of particles
        pid=parse(Int16, headerdata[2]),        # Process type
        weight=parse(Float64, headerdata[3]),   # Event weight
        scale=parse(Float64, headerdata[4]),    # Scale
        aqed=parse(Float64, headerdata[5]),     # AQEDUP
        aqcd=parse(Float64, headerdata[6]),     # AQCDUP
    )
    particles = [
        begin
            fields = split(line, ' '; keepempty=false)
            p = (;
                idx=idx - 1, # zero-based to match `mother1`, `mother2`
                id=parse(Int32, fields[1]),
                status=parse(Int8, fields[2]),
                mother1=parse(Int16, fields[3]),
                mother2=parse(Int16, fields[4]),
                color1=parse(Int32, fields[5]),
                color2=parse(Int32, fields[6]),
                px=parse(Float64, fields[7]),
                py=parse(Float64, fields[8]),
                pz=parse(Float64, fields[9]),
                e=parse(Float64, fields[10]),
                m=parse(Float64, fields[11]),
                lifetime=parse(Float64, fields[12]),
                spin=parse(Float64, fields[13]),
            )
            p
        end for (idx, line) in enumerate(lines[2:(2 + header.nparticles - 1)])
    ]
    return Event(header, particles)
end

function parse_lhe(filename)
    reader = open(EzXML.StreamReader, filename)
    function f(item)
        return (item != nothing) &&
               (reader.name == "event") &&
               (reader.type == EzXML.READER_ELEMENT)
    end
    return (parse_event(reader.content) for _ in Iterators.filter(f, reader))
end

function flatparticles(filename)
    return vcat(
        [
            [(; eventnum=ievt, p...) for p in evt.particles] for
            (ievt, evt) in enumerate(parse_lhe(filename))
        ]...,
    )
end

end # module


text:
"""
    rayleigh(λs, nH_I, nHe_II, nH2)

Absorption coefficient from Rayleigh scattering by neutral H, He, and H2.  Formulations for H and He
are via [Colgan+ 2016](https://ui.adsabs.harvard.edu/abs/2016ApJ...817..116C/abstract).  Formulation
for H2 from
[Dalgarno and Williams 1962](https://ui.adsabs.harvard.edu/abs/1962ApJ...136..690D/abstract).

The Dalgarno and Williams H2 is applicable redward of 1300 Å.  Since Rayleigh scattering breaks down
when the particle size to wavelength ratio gets large, we that all frequencies passed to this
function be equivalent to 1300 Å or greater.

The formulations for H is adapted from
[Lee 2005](https://ui.adsabs.harvard.edu/abs/1962ApJ...136..690D/abstract), which states that it is
applicable redward of Lyman alpha. See Colgan 2016 for details on He.
"""
function rayleigh(νs::AbstractVector{<:Real}, nH_I, nHe_I, nH2)
    @assert c_cgs/maximum(νs) > 1.3e-5
    σth = 6.65246e-25 #Thompson scattering cross section [cm^2]

    #(ħω/ 2E_H)^2 in Colgan+ 2016.  The photon energy over 2Ryd
    E_2Ryd_2 = @. (hplanck_eV * νs / (2 * Rydberg_eV ))^2
    E_2Ryd_4 = E_2Ryd_2.^2
    E_2Ryd_6 = E_2Ryd_2.*E_2Ryd_4
    E_2Ryd_8 = E_2Ryd_4.^2

    #Colgan+ 2016 equation 6
    σH_σth = @. 20.24E_2Ryd_4 + 239.2E_2Ryd_6 + 2256E_2Ryd_8
    #Colgan+ 2016 equation 7
    σHe_σth = @. 1.913E_2Ryd_4 + 4.52E_2Ryd_6 + 7.90E_2Ryd_8

    αH_HE = @. (nH_I*σH_σth + nHe_I*σHe_σth) * σth

    #Dalgarno & Williams 1962 equation 3 (which assumes λ in Å)
    invλ2 = (νs ./ (1e8 * c_cgs)) .^2
    invλ4 = invλ2 .^ 2
    invλ6 = invλ2 .* invλ4
    invλ8 = invλ4 .^ 2
    αH2 = @. (8.14e-13*invλ4 + 1.28e-6*invλ6 + 1.61*invλ8) * nH2

    αH_HE .+ αH2
end

"""
    electron_scattering(nₑ)

Compute the linear absorption coefficient, α, from scattering off of free electrons. This has no
wavelength dependence. It assumes isotropic scattering.  (See, e.g. Gray p 160.)

# Arguments
- `nₑ::F`: number density of free electrons (in cgs)
"""
function electron_scattering(nₑ::F) where {F<:Real}
    8π/3 * (electron_charge_cgs^2/(electron_mass_cgs * c_cgs^2))^2 * nₑ
end


text:
using Pkg
Pkg.instantiate()
using JSON, Statistics, GitHub, Base64, SHA, Downloads, Dates
cd(@__DIR__)
include("benchmark-library.jl")

ctx = github_context()

url = upload_file(ctx, "test.png", "test3.png")

comment = """
## Compile Times compared to tagged

![]($(url))
"""

pr = GitHub.PullRequest(1691)

make_or_edit_comment(ctx, pr, comment)

# save("test.png", series(rand(4,10)))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Glibc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Glibc")
JLLWrappers.@generate_main_file("Glibc", UUID("452aa2e7-e185-58db-8ff9-d3c1fa4bc997"))
end  # module Glibc_jll


text:
using CUDAdrv, CUDAnative

struct CUDABackend{R<:CUDALibs} <: Backend{R}
    resource::R
end

function backend_init(resource::CUDALibs, domain, sim)
    CUDABackend(resource)
end

function backend_cleanup!(backend::CUDABackend, state, sim)
    nothing
end

function move_backend(A::Array{T,N}, backend::CUDABackend) where {T<:Union{Float32,Float64},N}
    CuArray(A)
end

function move_backend(A::Array{Int8,N}, backend::CUDABackend) where N
    move_backend(Float32.(A), backend)
end

# --------------------------------------------------------------------

function backend_update!(state::BoxState{2,T}, backend::CUDABackend, sim) where T
    h,w = size(state.current)
    threads = (32,32)
    yblocks = ceil(Int, h / threads[2])
    xblocks = ceil(Int, w / threads[1])
    # shmem = 2 * prod(threads.+2) * sizeof(T)
    @cuda blocks=(xblocks,yblocks) threads=threads cuda_kernel!(
        state.previous, state.current, state.q,
        sim.wave.λ,
        state.box.γ)
    state
end

function cuda_kernel!(Ψₜ₋₁::CuDeviceArray{T,2}, Ψₜ::CuDeviceArray{T,2}, q::CuDeviceArray, λ::T, γ::T) where T
    y = (blockIdx().y-1) * blockDim().y + threadIdx().y
    x = (blockIdx().x-1) * blockDim().x + threadIdx().x
    h, w = size(Ψₜ₋₁)
    Ψₜ₊₁ = Ψₜ₋₁ # store next state in previous
    eins = one(T)
    zwei = eins + eins
    λsq = λ*λ
    λhalf = λ/zwei

    @inbounds if 1 < y < h && 1 < x < w
        q₁₋ = q[y-1,x]
        q₁₊ = q[y+1,x]
        q₂₋ = q[y,x-1]
        q₂₊ = q[y,x+1]
        K = q₁₋ + q₁₊ + q₂₋ + q₂₊

        Q = q₁₋ * Ψₜ[y-1,x] + q₁₊ * Ψₜ[y+1,x] +
            q₂₋ * Ψₜ[y,x-1] + q₂₊ * Ψₜ[y,x+1]

        γ̄ = γ * ((eins-q₁₋) + (eins-q₁₊) + (eins-q₂₋) + (eins-q₂₊))

        Ψₜ₊₁[y,x] = eins/(eins+λhalf*γ̄) * ((zwei-λsq*K) * Ψₜ[y,x] + (λhalf*γ̄-eins) * Ψₜ₋₁[y,x] + λsq*Q)
    end

    return nothing
end

# --------------------------------------------------------------------

function backend_update!(state::BoxState{3,T}, backend::CUDABackend, sim) where T
    h,w,v = size(state.current)
    threads = (16,8,8)
    yblocks = ceil(Int, h / threads[2])
    xblocks = ceil(Int, w / threads[1])
    zblocks = ceil(Int, v / threads[3])
    #shmem = Int(2 * prod(threads.+2) * sizeof(T))
    @cuda blocks=(xblocks,yblocks,zblocks) threads=threads cuda_kernel!(
        state.previous, state.current, state.q,
        sim.wave.λ,
        state.box.γ)
    state
end

function cuda_kernel!(Ψₜ₋₁::CuDeviceArray{T,3}, Ψₜ::CuDeviceArray{T,3}, q::CuDeviceArray, λ::T, γ::T) where T
    y = (blockIdx().y-1) * blockDim().y + threadIdx().y
    x = (blockIdx().x-1) * blockDim().x + threadIdx().x
    z = (blockIdx().z-1) * blockDim().z + threadIdx().z
    h, w, v = size(Ψₜ₋₁)
    Ψₜ₊₁ = Ψₜ₋₁ # store next state in previous
    eins = one(T)
    zwei = eins + eins
    λsq = λ*λ
    λhalf = λ/zwei

    @inbounds if 1 < y < h && 1 < x < w && 1 < z < v
        q₁₋ = q[y-1,x,z]
        q₁₊ = q[y+1,x,z]
        q₂₋ = q[y,x-1,z]
        q₂₊ = q[y,x+1,z]
        q₃₋ = q[y,x,z-1]
        q₃₊ = q[y,x,z+1]
        K = q₁₋ + q₁₊ + q₂₋ + q₂₊ + q₃₋ + q₃₊

        Q = q₁₋ * Ψₜ[y-1,x,z] + q₁₊ * Ψₜ[y+1,x,z] +
            q₂₋ * Ψₜ[y,x-1,z] + q₂₊ * Ψₜ[y,x+1,z] +
            q₃₋ * Ψₜ[y,x,z-1] + q₃₊ * Ψₜ[y,x,z+1]

        γ̄ = γ * ((eins-q₁₋) + (eins-q₁₊) + (eins-q₂₋) + (eins-q₂₊)+ (eins-q₃₋) + (eins-q₃₊))

        Ψₜ₊₁[y,x,z] = eins/(eins+λhalf*γ̄) * ((zwei-λsq*K) * Ψₜ[y,x,z] + (λhalf*γ̄-eins) * Ψₜ₋₁[y,x,z] + λsq*Q)
    end

    return nothing
end


text:
@testset "molgraph" begin

    @testset "graphmol" begin
        mol = graphmol(SDFileAtom, SDFileBond)
        @test nodecount(mol) == 0
        edges = [(1, 2)]
        atoms = [SDFileAtom(),SDFileAtom(),SDFileAtom()]
        bonds = [SDFileBond()]
        mol = graphmol(edges, atoms, bonds)
        newbond = addedge!(mol, 2, 3, SDFileBond())
        @test newbond == 2
        # TODO: updateedge! for plaingraph
        # updateedge!(mol, newbond, 2)
        # @test edgecount(mol) == 2
    end

    @testset "querymol" begin
        q = querymol(SmartsAtom, SmartsBond)
        addnode!(q, SmartsAtom(:Any => true))
        addnode!(q, SmartsAtom(:Any => true))
        addedge!(q, 1, 2, SmartsBond(:Any => true))
        @test nodecount(q) == 2
        @test hasedge(q, 1, 2)
    end

end


text:
@testset "send_close.jl" begin
    @testset "FIN sent test" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        sender_close(sender)
        expect_state(sender, JLSponge.FIN_SENT)
        @test sender.bytes_in_flight == 1
        expect_seg(sender, fin=true, seqno=WrappingInt32(isn + 1))
    end

    @testset "FIN acked test" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        sender_close(sender)
        expect_state(sender, JLSponge.FIN_SENT)
        expect_seg(sender, fin=true, seqno=WrappingInt32(isn + 1))
        ack_received_test(sender, WrappingInt32(isn + 2))
        expect_state(sender, JLSponge.FIN_ACKED)
        @test sender.bytes_in_flight == 0
        expect_no_seg(sender)
    end
    @testset "FIN not acked test" begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        sender_close(sender)
        expect_state(sender, JLSponge.FIN_SENT)
        expect_seg(sender, fin=true, seqno=WrappingInt32(isn + 1))
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.FIN_SENT)
        @test sender.bytes_in_flight == 1
        expect_no_seg(sender)
    end

    @testset "FIN retx test" begin
        TIMEOUT_DFLT = UInt16(1000)
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(;fixed_isn=WrappingInt32(isn), retx_timeout=TIMEOUT_DFLT)
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        sender_close(sender)
        expect_state(sender, JLSponge.FIN_SENT)
        expect_seg(sender, fin=true, seqno=WrappingInt32(isn + 1))
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.FIN_SENT)
        @test sender.bytes_in_flight == 1
        expect_no_seg(sender)
        tick!(sender, TIMEOUT_DFLT - 1)
        expect_state(sender, JLSponge.FIN_SENT)
        @test sender.bytes_in_flight == 1
        expect_no_seg(sender)

        tick!(sender, 1)
        expect_state(sender, JLSponge.FIN_SENT)
        @test sender.bytes_in_flight == 1
        expect_seg(sender, fin=true, seqno=WrappingInt32(isn + 1))
        expect_no_seg(sender)
        tick!(sender, 1)
        expect_state(sender, JLSponge.FIN_SENT)
        @test sender.bytes_in_flight == 1
        expect_no_seg(sender)

        ack_received_test(sender, WrappingInt32(isn + 2))
        expect_state(sender, JLSponge.FIN_ACKED)
        @test sender.bytes_in_flight == 0
        expect_no_seg(sender)
        tick!(sender, 1)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
###################################
"""
    fgmaxgrids = fgmaxdata("simlation/path/_output"::AbstractString)
    fgmaxgrids = fgmaxdata("simlation/path/_output/fgmax_grids.data"::AbstractString)

Function: fgmax_grids.data reader
"""
function fgmaxdata(outdir::AbstractString)
    ## definition of filename
    fname = occursin("fgmax_grids.data", basename(outdir)) ? outdir : joinpath(outdir, "fgmax_grids.data")

    ## check
    isfile(fname) || error("File $fname is not found.")
    ## read all lines
    open(fname,"r") do f
        global txt = readlines(f)
    end

    ## parse parameters
    num_fgmax_val = parse(Int64, split(txt[occursin.("num_fgmax_val",txt)][1],r"\s+")[1])
    num_fgmax_grids = parse(Int64, split(txt[occursin.("num_fgmax_grids",txt)][1],r"\s+")[1])

    ## preallocate
    fg = Vector{VisClaw.FixedGrid}(undef, num_fgmax_grids)
    ## load
    global baseline = 10
    for i = 1:num_fgmax_grids
        ## basic parameters
        FGid = parse(Int64, split(strip(txt[baseline+1]), r"\s+")[1])
        point_style = parse(Int64, split(strip(txt[baseline+8]), r"\s+")[1])
        ## check point style
        (point_style > 4 || point_style < 0) && error("point_style $point_style is not supported yet.")

        ## 0
        if point_style == 0
            ## npts
            npts = parse(Int64, split(strip(txt[baseline+9]), r"\s+")[1])
            ## x, y
            if npts == 0
                fgmax_deffile = String(strip(txt[baseline+10])[2:end-1])
                open(fgmax_deffile,"r") do ff; global fgtxt = readlines(ff); end
                npts = parse(Int64, split(strip(fgtxt[1]), r"\s+")[1])
                x = zeros(Float64,npts)
                y = zeros(Float64,npts)
                for ip = 1:npts
                    x[ip], y[ip] = parse.(Float64, split(strip(fgtxt[1+ip]), r"\s+")[1:2])
                end
                baseline += 11
            else
                x = zeros(Float64,npts)
                y = zeros(Float64,npts)
                for ip = 1:npts
                    x[ip], y[ip] = parse.(Float64, split(strip(txt[baseline+9+ip]), r"\s+")[1:2])
                end
                baseline += 10 + npts
            end
            ## instance
            fg[i] = VisClaw.FixedGrid(FGid, point_style, num_fgmax_val, npts, x, y)

        ## 1
        elseif point_style == 1
            npts = parse(Int64, split(strip(txt[baseline+9]), r"\s+")[1])
            x1, y1 = parse.(Float64, split(strip(txt[baseline+10]), r"\s+")[1:2])
            x2, y2 = parse.(Float64, split(strip(txt[baseline+11]), r"\s+")[1:2])
            x = LinRange(x1, x2, npts)
            y = LinRange(y1, y2, npts)
            # instance
            fg[i] = VisClaw.FixedGrid(FGid, point_style, num_fgmax_val, npts, x, y)
            baseline += 12

        ## 2
        elseif point_style == 2
            nx, ny = parse.(Int64, split(strip(txt[baseline+9]), r"\s+")[1:2])
            x1, y1 = parse.(Float64, split(strip(txt[baseline+10]), r"\s+")[1:2])
            x2, y2 = parse.(Float64, split(strip(txt[baseline+11]), r"\s+")[1:2])
            # instance
            fg[i] = VisClaw.FixedGrid(FGid, point_style, num_fgmax_val, nx, ny, (x1,x2), (y1,y2))
            baseline += 12

        ## 3
        elseif point_style == 3
            error("point_style $point_style is not supported yet.")

        ## 4
        elseif point_style == 4
            fgmax_deffile = String(strip(txt[baseline+9])[2:end-1])
            topoflag = loadtopo(fgmax_deffile, 3)

            flag = findall(convert(BitArray, topoflag.elevation))
            X = repeat(topoflag.x', outer=(topoflag.nrows,1))
            Y = repeat(topoflag.y,  outer=(1,topoflag.ncols))
            x = vec(X[flag])
            y = vec(Y[flag])
            npts = length(x)

            ind = sortperm([flag[i][1] for i=1:npts], rev=true)
            flag = flag[ind]
            x = x[ind]
            y = y[ind]

            ## instance
            fg[i] = VisClaw.FixedGrid(FGid, point_style, num_fgmax_val,
                                      topoflag.ncols, topoflag.nrows,
                                      extrema(topoflag.x), extrema(topoflag.y),
                                      npts, x, y, flag)
            baseline += 10
        end
    end

    ## return
    return fg
end


text:
module Utils

export deleteitems!, get_n_words, hascoordinates, select_element, matrix2list, list2matrix, check_file, isnotemptyfile,

# eachline,

Format,

AbstractTest, TestType, TestOperation, Is, In, Not, capture, isobject, findobjects, collectobjects, collectcaptures, guess_type,

Commandline

include("generalutils.jl")
include("EachLineString.jl")
include("Read.jl")
include("Write.jl")
include("FindObjects.jl")

# COMMAND LINE: Scripts Module
# ============================

module Scripts

export parse_commandline,
runscript, run_single_script, script,
set_parallel, open_output, close_output,
readorparse

using ArgParse, GZip

include("scripts.jl")

end # Scripts

end # Utils


text:
module GapFindFill

using JuMP, Gurobi, CPLEX, GLPK
include("GFFJInterface.jl")

# Export the interface function -
export find_gaps, fill_gaps_min

end # module


text:
#### options.jl ####
# A framework for providing optional arguments to functions

# Harlan Harris & Timothy E. Holy, with contributions from Stefan
# Karpinski, Patrick O'Leary, and Jeff Bezanson
module OptionsMod
using Base
# can't get Base.ht_keyindex from dict.jl -- will pull it manually

import Base: convert, copy, show, getindex, setindex!

export Options,
	CheckNone, CheckWarn, CheckError,
	add_defaults!,
	ischeck, docheck, clearcheck,
	@options, @defaults, @check_used, @set_options

#### Options type ####
abstract OptionsChecking
type CheckNone <: OptionsChecking
end
type CheckWarn <: OptionsChecking
end
type CheckError <: OptionsChecking
end
type Options{T<:OptionsChecking}
    key2index::Dict{Symbol,Int}
    vals::Vector
    used::BitArray
    check_lock::BitArray
end
# Constructor: supply type followed by parameter/value pairs,
#   o = Options(CheckWarn,:a,5,:b,rand(3),...)
# Note, the macro @options makes construction easier
function Options{T<:OptionsChecking}(::Type{T},args...)
    if isodd(length(args))
        error("Parameter/value lists must come in pairs")
    end
    n = div(length(args),2)
    keys, index, vals = if n > 0
        (args[1:2:end], tuple(collect(1:n)...), Any[args[2:2:end]...])
    else
        ((), (), Array(Any, 0))
    end
    if VERSION < v"0.3-"
        ht = Dict{Symbol,Int}(keys,index)
    else
        ht = Dict{Symbol,Int}(zip(keys,index))
    end
    used = falses(n)
    check_lock = falses(n)
    Options{T}(ht,vals,used,check_lock)
end
# Constructor: supply type followed by list of assignment expressions, e.g.,
#   o = Options(CheckNone,:(a=5),:(b=rand(3)),...)
function Options{T<:OptionsChecking}(::Type{T},ex::Expr...)
    ht = Dict{Symbol,Int}()
    vals = Array(Any,0)
    n = length(ex)
    for i = 1:n
        if ex[i].head != :(=)
            error("All expressions must be assignments")
        end
        ht[ex[i].args[1]] = i
        push!(vals,ex[i].args[2])
    end
    used = falses(n)
    check_lock = falses(n)
    Options{T}(ht,vals,used,check_lock)
end
# If no type specified, CheckError is the default
Options(p1::Symbol,v1,args...) = Options(CheckError,p1,v1,args...)
Options(ex::Expr...) = Options(CheckError,ex...)

function copy{T<:OptionsChecking}(o::Options{T})
    Options{T}(copy(o.key2index),
            copy(o.vals),
            copy(o.used),
            copy(o.check_lock))
end
function convert{Tnew<:OptionsChecking,Told<:OptionsChecking}(::Type{Tnew},o::Options{Told})
    Options{Tnew}(o.key2index,o.vals,o.used,o.check_lock)
end
function show{T<:OptionsChecking}(io::IO, o::Options{T})
    # Put them in the same order specified by the user
    key = Array(ASCIIString,length(o.vals))
    for (k, v) = o.key2index
        key[v] = string(k)
    end
    for i = 1:length(key)
        print(io, "$(key[i]) = $(o.vals[i])")
        if i < length(key)
            print(io, ", ")
        end
    end
    print(io, " ($T)")
end

#### Functions ####
# Return an options setting
function getindex(o::Options,s::Symbol)
    index = Base.ht_keyindex(o.key2index,s)
    if index > 0
        index = o.key2index.vals[index]
        return o.vals[index]
    else
        return nothing
    end
end
# Re-set an options setting, or add a new one
function setindex!(o::Options,v,s::Symbol)
    index = Base.ht_keyindex(o.key2index,s)
    if index > 0
        index = o.key2index.vals[index]
        o.vals[index] = v
        o.used[index] = false
        o.check_lock[index] = false
    else
        o.key2index[s] = length(o.vals)+1
        push!(o.vals,v)
        push!(o.used,false)
        push!(o.check_lock,false)
    end
end
# Functions for "claiming" and checking usage of individual options
function ischeck(o::Options)
    ret = !o.check_lock   # items marked true are to be checked by caller
    o.check_lock[ret] = true
    return ret
end
function ischeck(o::Options{CheckNone})
    return falses(length(o.vals))
end
function ischeck(a)
    error("First argument must be an options type")
end
function docheck_common(o::Options,checkflag::BitArray)
    unused = checkflag & !o.used[1:length(checkflag)]
    msg = ""
    if any(unused)
        s = Array(ASCIIString,0)
        for (k, v) = o.key2index
            if v <= length(unused) && unused[v]
                push!(s,string(k))
            end
        end
        msg = string("The following option(s) were unused: ",s)
    end
    return unused, msg
end
function docheck(o::Options{CheckNone},checkflag::BitArray)
    clearcheck(o,checkflag)
end
function docheck(o::Options{CheckWarn},checkflag::BitArray)
    unused, msg = docheck_common(o,checkflag)
    if any(unused)
        println("Warning: ",msg)
    end
    clearcheck(o,checkflag)
end
function docheck(o::Options{CheckError},checkflag::BitArray)
    unused, msg = docheck_common(o,checkflag)
    clearcheck(o,checkflag)  # in case it's in a try/catch block...
    if any(unused)
        error(msg)
    end
end
# Reset status on handled options (in case o is reused later)
function clearcheck(o::Options,checkflag::BitArray)
    for i = 1:length(checkflag)
        if checkflag[i]
            o.used[i] = false
            o.check_lock[i] = false
        end
    end
end

#### Macros ####
# Macro to set the defaults. Usage:
#     @defaults opts a=3 b=a+2 c=f(a,b)
# After executing this macro, you can use a, b, and c as ordinary variables
macro defaults(opts,ex...)
    # Create a new variable storing the checkflag
    varname = string("_",string(opts),"_checkflag")
    exret = :($(esc(symbol(varname))) = ischeck($(esc(opts))))
    # Transform the tuple into a vector, so that
    # we can manipulate it
    ex = Any[ex...]
    # Check each argument in the assignment list
    i = 1
    while i <= length(ex)
        y = ex[i]
        if isa(y, Expr) && y.head == :block
            # Found a begin..end block: expand its contents in-place
            # and restart from the same position
            splice!(ex, i, y.args)
            continue
        elseif isa(y, LineNumberNode) || (isa(y,Expr) && y.head == :line)
            # A line number node, ignore
            i += 1
            continue
        elseif !isa(y,Expr) || !(y.head == :(=) || y.head == :(=>) || y.head == :(:=) || y.head == :kw)
            error("Arguments to @defaults following the options variable must be assignments, e.g., a=5 or a=>5")
        end
        y.head = :(=)

        sym = y.args[1]
        exret = quote
            $exret
            htindex = Base.ht_keyindex($(esc(opts)).key2index,$(Expr(:quote,sym)))
            if htindex > 0
                htindex = $(esc(opts)).key2index.vals[htindex]
                $(esc(sym)) = $(esc(opts)).vals[htindex]
                $(esc(opts)).used[htindex] = true
            else
                $(esc(y))
            end
        end
        i += 1
    end
    exret
end

# Macro to check whether options were used (optional, use only if you
# want to protect the user from typos or other nonsensical settings)
# Usage:
#    @check_used opts
macro check_used(opts)
    varname = string("_",string(opts),"_checkflag")
    :(docheck($(esc(opts)),$(esc(symbol(varname)))))
end

# Macro for setting options. Usage:
#    opts = @options a=5 b=7
#    opts = @options CheckWarn a=5 b=7
macro options(ex...)
    callargs = Any[:Options]
    # Transform the tuple into a vector, so that
    # we can manipulate it
    ex = Any[ex...]
    i = 1
    if length(ex) >= 1 && isa(ex[1], Symbol)
        push!(callargs, esc(ex[1]))
        i += 1
    end
    while i <= length(ex)
        y = ex[i]
        if isa(y, Expr) && y.head == :block
            # Found a begin..end block: expand its contents in-place
            # and restart from the same position
            splice!(ex, i, y.args)
            continue
        elseif isa(y, LineNumberNode) || (isa(y,Expr) && y.head == :line)
            # A line number node, ignore
            i += 1
            continue
        elseif !isa(y,Expr) || !(y.head == :(=) || y.head == :(=>) || y.head == :(:=) || y.head == :kw)
            error("Arguments to @options must be assignments, e.g., a=5 or a=>5")
        end
        push!(callargs, Expr(:quote, y.args[1]))
        push!(callargs, esc(y.args[2]))
        i += 1
    end
    Expr(:call, callargs...)
end

# Macro to reset or append additional options. Usage:
#    @set_options opts c=25 d=4
macro set_options(opts,ex...)
    exret = quote
    end
    # Transform the tuple into a vector, so that
    # we can manipulate it
    ex = Any[ex...]
    # Check each argument in the assignment list
    i = 1
    while i <= length(ex)
        y = ex[i]
        if isa(y, Expr) && y.head == :block
            # Found a begin..end block: expand its contents in-place
            # and restart from the same position
            splice!(ex, i, y.args)
            continue
        elseif isa(y, LineNumberNode) || (isa(y,Expr) && y.head == :line)
            # A line number node, ignore
            i += 1
            continue
        elseif !isa(y,Expr) || !(y.head == :(=) || y.head == :(=>) || y.head == :(:=) || y.head == :kw)
            error("Arguments to @set_options following the options variable must be assignments, e.g., a=5 or a=>5")
        end
        y.head = :(=)

        sym = y.args[1]
        val = y.args[2]
        exret = quote
            $exret
            $(esc(opts))[$(Expr(:quote,sym))] = $(esc(val))
        end
        i += 1
    end
    exret
end

end # module


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    promote_rule(Exo1Real{Int}, Float64) == Exo1Real{Float64}
end

@test begin
    promote_rule(Exo1Real{Int}, Exo1Real{Float64}) == Exo1Real{Float64}
end

@test begin
    -(1.0 - Exo1Real(2, 3)) == Exo1Real(1.0, 3.0) == Exo1Real(2, 3) - 1.0
end

@test begin
    1.0 + Exo1Real(2, 3) == Exo1Real(3.0, 3.0) == Exo1Real(2, 3) + 1.0
end

@test begin
    convert(Exo1Real{Float64}, 1) == Exo1Real(1.0, 0)
end

@test begin
    convert(Exo1Real{Float64}, Exo1Real(1,2)) == Exo1Real(1.0, 2.0)
end

@test begin
    convert(Exo1Real, 2) == Exo1Real(2, 0)
end

@test begin
    z = Exo1Real{Int}(1, 2)
    convert(Parabolic{Int}, z) == z
end

@test begin
    convert(Int, Exo1Real(3)) == 3
end

@test_throws InexactError begin
    convert(Int, Exo1Real(3,4))
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(Exo1Real(a))
end

@test begin
    a = 1
    !isreal(Exo1Real(a, a))
end

@test begin
    z = random(Exo1Real{Int})
    z == +(z)
end

@test begin
    a = 1
    real(Exo1Real(a)) == a
end

@test begin
    a = 1
    l = unreal(Exo1Real(0, a))
    r = [a]
    l == r
end

@test begin
    io = IOBuffer()
    show(io, Exo1Real(1, 2))
    l = String(take!(io))
    r = "[1: 1, a: 2]"
    l == r
end

@test begin
    x = random(Exo1Real{BigFloat})
    y = random(Exo1Real{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo1Real{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(Exo1Real{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    z = random(Exo1Real{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    z = random(Exo1Real{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(Exo1Real{BigInt})
    cloak(cloak(x)) == x
end

@test begin
    x = random(Exo1Real{BigInt})
    dagger(dagger(x)) == x
end

@test begin
    x = random(Exo1Real{BigInt})
    star(star(x)) == x
end

@test begin
    x = random(Exo1Real{Rational{BigInt}})
    selfstar(selfstar(x)) == selfstar(x)
end

@test begin
    x = random(Exo1Real{Rational{BigInt}})
    anti_selfstar(anti_selfstar(x)) == anti_selfstar(x)
end

@test begin
    x = random(Exo1Real{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(Exo1Real{Rational{BigInt}})
    y = random(Exo1Real{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(Exo1Real{Rational{BigInt}})
    y = random(Exo1Real{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(Exo1Real{Rational{BigInt}})
    y = random(Exo1Real{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    z = random(Exo1Real{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(Exo1Real{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    z = random(Exo1Real{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    z = random(Exo1Real{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    z = random(Exo1Real{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(Exo1Real{BigInt})
    quadrance(x) > big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(Exo1Real{BigInt})
    y = random(Exo1Real{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(Exo1Real(0, 0.0))
end

@test_throws ErrorException begin
    inv(Exo1Real(0, 1.0))
end

@test_throws ErrorException begin
    random(Exo1Real{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(Exo1Real{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo1Real{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(Exo1Real{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(Exo1Real{BigInt})
    z == z + zero(Exo1Real{BigInt})
end

@test begin
    z = random(Exo1Real{BigInt})
    z == z * one(Exo1Real{BigInt})
end

@test begin
    z = random(Exo1Real{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(Exo1Real{BigInt})
    z == z * one(z)
end

@test begin
    float(Exo1Real{Float64}) == Exo1Real{Float64}
end


text:
using Test

using Yao
using YaoExtensions

@testset "constructor" begin

    g = Sequence(
        kron(2, X, Y),
        kron(2, 1=>phase(0.1)),
    )

    @test g isa Sequence
    @test g.blocks == [kron(2, X, Y), kron(2, 1=>phase(0.1))]
end

@testset "apply" begin
    g = Sequence(
        kron(2, X, Y),
        kron(2, 1=>phase(0.1)),
    )

    reg = rand_state(2)
    @test statevec(apply!(copy(reg), g)) ≈ mat(chain(g...)) * reg.state
end

@testset "iteration" begin
    test_list = [X, Y, phase(0.1), rot(X, 0.0)]
    g = Sequence(test_list)

    for (src, tg) in zip(g, test_list)
        @test src == tg
    end

    for (src, tg) in zip(eachindex(g), 1:length(test_list))
        @test src == tg
    end
end

@testset "additional" begin
    g = Sequence(X, Y)
    push!(g, Z)
    @test g[3] == Z

    append!(g, [rot(X, 0.0), rot(Y, 0.0)])
    @test g[4] == rot(X, 0.0)
    @test g[5] == rot(Y, 0.0)

    prepend!(g, [phase(0.1)])
    @test g[1] == phase(0.1)
    @test g[2] == X
    @test g[end] == rot(Y, 0.0)
    gg = insert!(g, 4, Z)
    @test gg[4] == Z
end

@testset "traits" begin
    # TODO: check traits when primitive blocks' traits are all defined
    g = Sequence(X, Y)
    @test length(g) == 2
    @test eltype(g) == eltype(g.blocks)
end


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "EMBRAKE"
cases = ["BRKDC01", "BRKDC01-D", "BRKNC01", "BRKNC01-D", "BRKNP01", "BRKNP01-D"]
SUITE[model] = BenchmarkGroup()

include("embrake.jl")
validation = []
times = []

# ----------------------------------------
#  BRKDC01
# ----------------------------------------

prob_no_pv_no_jit = embrake_no_pv(ζ=0., Tsample=1e-4)
alg = GLGM06(δ=1e-7, max_order=1, static=true, dim=4, ngens=4)
sol_no_pv_no_jit = solve(prob_no_pv_no_jit, max_jumps=1001, alg=alg)

# verify that specification holds
property = ρ(eₓ, sol_no_pv_no_jit) < x0
push!(validation, Int(property))

t = max_time(sol_no_pv_no_jit)
push!(times, trunc(t, digits=4))
println("maximum time that x < x0 , case $(cases[1]) : $t")

# benchmark
SUITE[model][cases[1]] = @benchmarkable solve($prob_no_pv_no_jit, max_jumps=1001, alg=$alg)

# ----------------------------------------
#  BRKDC01 (discrete-time)
# ----------------------------------------

alg = GLGM06(δ=1e-7, max_order=1, static=true, dim=4, ngens=4, approx_model=NoBloating())
sol_no_pv_no_jit_discrete = solve(prob_no_pv_no_jit, max_jumps=1001, alg=alg)

# verify that specification holds
property = ρ(eₓ, sol_no_pv_no_jit_discrete) < x0
push!(validation, Int(property))

t = max_time(sol_no_pv_no_jit_discrete)
push!(times, trunc(t, digits=4))
println("maximum time that x < x0 , case $(cases[2]) : $t")

# benchmark
SUITE[model][cases[2]] = @benchmarkable solve($prob_no_pv_no_jit, max_jumps=1001, alg=$alg)

sol_no_pv_no_jit_discrete = nothing
GC.gc()

# ----------------------------------------
#  BRKNC01
# ----------------------------------------

prob_no_pv_with_jit = embrake_no_pv(ζ=[-1e-8, 1e-7], Tsample=1e-4)
alg = GLGM06(δ=1e-8, max_order=1, static=true, dim=4, ngens=4)
sol_no_pv_with_jit = solve(prob_no_pv_with_jit, max_jumps=1001, alg=alg)

# verify that specification holds
property = ρ(eₓ, sol_no_pv_with_jit) < x0
push!(validation, Int(property))

GC.gc()
t = max_time(sol_no_pv_with_jit)
push!(times, trunc(t, digits=4))
println("maximum time that x < x0 , case $(cases[3]) : $t")

# benchmark
SUITE[model][cases[3]] = @benchmarkable solve($prob_no_pv_with_jit, max_jumps=1001, alg=$alg)

sol_no_pv_with_jit = nothing
GC.gc()

# ----------------------------------------
#  BRKNC01 (discrete-time)
# ----------------------------------------

alg = GLGM06(δ=1e-8, max_order=1, static=true, dim=4, ngens=4, approx_model=NoBloating())
sol_no_pv_no_jit_discrete = solve(prob_no_pv_no_jit, max_jumps=1001, alg=alg)

# verify that specification holds
property = ρ(eₓ, sol_no_pv_no_jit_discrete) < x0
push!(validation, Int(property))

GC.gc()
t = max_time(sol_no_pv_no_jit_discrete)
push!(times, trunc(t, digits=4))
println("maximum time that x < x0 , case $(cases[4]) : $t")

# benchmark
SUITE[model][cases[4]] = @benchmarkable solve($prob_no_pv_no_jit, max_jumps=1001, alg=$alg)

sol_no_pv_no_jit_discrete = nothing
GC.gc()

# ----------------------------------------
#  BRKNP01
# ----------------------------------------

prob_pv_1_with_jit = embrake_pv_1(ζ=[-1e-8, 1e-7], Tsample=1e-4)
alg = ASB07(δ=1e-8, max_order=1, static=true, dim=4, ngens=4)
sol_pv_1_with_jit = solve(prob_pv_1_with_jit, max_jumps=1001, alg=alg)

# verify that specification holds
property = ρ(eₓ, sol_pv_1_with_jit) < x0
push!(validation, Int(property))

GC.gc()
t = max_time(sol_pv_1_with_jit)
push!(times, trunc(t, digits=4))
println("maximum time that x < x0 , case $(cases[5]) : $t")

# benchmark
SUITE[model][cases[5]] = @benchmarkable solve($prob_pv_1_with_jit, max_jumps=1001, alg=$alg)

sol_pv_1_with_jit = nothing
GC.gc()

# ----------------------------------------
#  BRKNP01 (discrete-time)
# ----------------------------------------

alg = ASB07(δ=1e-8, max_order=1, static=true, dim=4, ngens=4, approx_model=NoBloating(exp=:interval))
sol_pv_1_with_jit_discrete = solve(prob_pv_1_with_jit, max_jumps=1001, alg=alg)

# verify that specification holds
property = ρ(eₓ, sol_pv_1_with_jit_discrete) < x0
push!(validation, Int(property))

GC.gc()
t = max_time(sol_pv_1_with_jit_discrete)
push!(times, trunc(t, digits=4))
println("maximum time that x < x0 , case $(cases[6]) : $t")

# benchmark
SUITE[model][cases[6]] = @benchmarkable solve($prob_pv_1_with_jit, max_jumps=1001, alg=$alg)

sol_pv_1_with_jit_discrete = nothing
GC.gc()

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
    t = median(results[model][c]).time * 1e-9
    runtimes[c] = round(t, digits=4)
end

for (i, c) in enumerate(cases)
    print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c]), $(times[i])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

polys = Vector{VPolygon{Float64, Vector{Float64}}}()
for fp in sol_no_pv_no_jit
    for (j, X) in enumerate(fp)
        sfpos, sfneg = ρ(eₓ, X), ρ(-eₓ, X)
        dt = tspan(fp, j)
        ti, tf = inf(dt), sup(dt)
        p = VPolygon([[ti, sfpos], [ti, -sfneg], [tf, sfpos], [tf, -sfneg]])
        push!(polys, p)
    end
end

fig = Plots.plot()

plot!(fig, [X for X in polys[1:500:end]], color=:black, lw=1.0, linecolor=:black,
        tickfont=font(30, "Times"), guidefontsize=45,
        xlab=L"t",
        ylab=L"x",
        xtick=[0.025, 0.05, 0.075, 0.1], ytick=[0.0, 0.01, 0.02, 0.03, 0.04, 0.05],
        xlims=(0.0, 0.1), ylims=(0.0, 0.05),
        bottom_margin=6mm, left_margin=2mm, right_margin=12mm, top_margin=3mm,
        size=(1000, 1000))

savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-EMBrake.png"))
savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-EMBrake.pdf"))

sol_no_pv_no_jit = nothing
GC.gc()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CpuTemplate <: SwaggerModel

    function CpuTemplate(;)
        o = new()
        o
    end
end # type CpuTemplate

const _name_map_CpuTemplate = Dict{String,Symbol}([])
const _field_map_CpuTemplate = Dict{Symbol,String}([])
Swagger.name_map(::Type{ CpuTemplate }) = _name_map_CpuTemplate
Swagger.field_map(::Type{ CpuTemplate }) = _field_map_CpuTemplate

function check_required(o::CpuTemplate)
    true
end

function validate_field(o::CpuTemplate, name::Symbol, val)
end