text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using CSV, StatsPlots
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.22"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data1 = CSV.read(datadir("machine1.csv"), header=false)[:,1]
	data2 = CSV.read(datadir("machine2.csv"), header=false)[:,1]
	data3 = CSV.read(datadir("machine3.csv"), header=false)[:,1]
end

# ╔═╡ 2ac4c7a6-250e-11eb-34b3-3b420e053d44
boxplot([data1,data2,data3], c=[:blue :red :green], label="",
	xticks=([1:1:3;],["1", "2", "3"]), xlabel="Machine type",
	ylabel="Pipe Diameter (mm)")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.22"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═2ac4c7a6-250e-11eb-34b3-3b420e053d44
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
export DubinsCar
export init_traj_nothing, init_traj_straightline

mutable struct DubinsCar <: DynamicsModel
	x_dim	# x, y, θ
	u_dim	# turning rate control

	v 		# speed
	k 		# curvature parameter

	x_min
	x_max
	u_min
	u_max
	clearance

	f::Vector
	A::Vector
	B
end

function DubinsCar()
	x_dim = 3
	u_dim = 1
	v = 2.
	k = 1.
	x_max = [100., 100., 2*π]
	x_min = -x_max
	u_max = 10.
	u_min = -u_max
	clearance = 0.01
	DubinsCar(x_dim, u_dim, v, k, x_min, x_max, u_min, u_max, clearance, [], [], [])
end

function SCPParam(model::DubinsCar, fixed_final_time::Bool)
  convergence_threshold = 1e-4
  SCPParam(fixed_final_time, convergence_threshold)
end

function SCPParam_GuSTO(model::DubinsCar)
  Δ0 = 10000.
  ω0 = 1.
  ω_max = 1.0e10
  ε = 1.0e-6
  ρ0 = 0.4
  ρ1 = 1.5
  β_succ = 2.
  β_fail = 0.5
  γ_fail = 5.

  SCPParam_GuSTO(Δ0, ω0, ω_max, ε, ρ0, ρ1, β_succ, β_fail, γ_fail)
end

function cost_true(traj, traj_prev::Trajectory, OAP::A) where A <: OptAlgorithmProblem{Car, DubinsCar, E} where E
	U, N = traj.U, OAP.N
  dtp = traj_prev.dt 	# TODO(ambyld): Change to current trajectory dt?

  Jm = 0

  # Forward Euler
  # Jm += sum(dt*U[k-1]^2 for k = 2:N)

  # Trapezoidal
  Jm += sum(1/2*dtp*(U[k-1]^2 + U[k]^2) for k = 2:N)
end

function cost_true_convexified(traj, traj_prev::Trajectory, OAP::A) where A <: OptAlgorithmProblem{Car, DubinsCar, E} where E
	cost_true(traj, traj_prev, OAP)
end

#############################
# Trajectory Initializations
#############################
function init_traj_nothing(TOP::TrajectoryOptimizationProblem{Car, DubinsCar, E}) where E
	model, x_init, goal_set = TOP.PD.model, TOP.PD.x_init, TOP.PD.goal_set
	x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess

	goal_final = get_first_goal_at_time(goal_set, tf_guess)
	x_goal = goal_final.params.point

	X = repmat(0.5(x_init + x_goal),1,N)
	U = zeros(u_dim,N)
	Trajectory(X, U, tf_guess)
end

function init_traj_straightline(TOP::TrajectoryOptimizationProblem{Car, DubinsCar, E}) where E
	model, x_init, goal_set = TOP.PD.model, TOP.PD.x_init, TOP.PD.goal_set
	x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess

	goal_final = get_first_goal_at_time(goal_set, tf_guess)
	x_goal = goal_final.params.point

	X = hcat(range(x_init, stop=x_goal, length=N)...)
	U = zeros(u_dim,N)
	Trajectory(X, U, tf_guess)
end

####################
# Constraint-adding
####################
function initialize_model_params!(SCPP::SCPProblem{Car, DubinsCar, E}, traj_prev::Trajectory) where E
	N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
	x_dim, u_dim = model.x_dim, model.u_dim
	Xp, Up = traj_prev.X, traj_prev.U

	model.f, model.A, model.B = [], [], B_dyn(Xp[:,1],robot,model)
	for k = 1:N
		push!(model.f, f_dyn(Xp[:,k],Up[:,k],robot,model))
		push!(model.A, A_dyn(Xp[:,k],robot,model))
	end
end

function update_model_params!(SCPP::SCPProblem{Car, DubinsCar, E}, traj_prev::Trajectory) where E
	N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
	Xp, Up, f, A = traj_prev.X, traj_prev.U, model.f, model.A

	for k = 1:N
		update_f!(f[k], Xp[:,k], Up[:,k], robot, model)
		update_A!(A[k], Xp[:,k], robot, model)
	end
end

macro constraint_abbrev_dubinscar(traj, traj_prev, SCPP)
	quote
		X, U, Tf = $(esc(traj)).X, $(esc(traj)).U, $(esc(traj)).Tf
    Xp, Up, Tfp, dtp = $(esc(traj_prev)).X, $(esc(traj_prev)).U, $(esc(traj_prev)).Tf, $(esc(traj_prev)).dt
    robot, model, WS, x_init, goal_set = $(esc(SCPP)).PD.robot, $(esc(SCPP)).PD.model, $(esc(SCPP)).WS, $(esc(SCPP)).PD.x_init, $(esc(SCPP)).PD.goal_set
    x_dim, u_dim, N, dh = model.x_dim, model.u_dim, $(esc(SCPP)).N, $(esc(SCPP)).dh
    X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh
	end
end

## Dynamics constraints
function dynamics_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Car, DubinsCar, E}, k::Int) where E
	# Where i is the state index, and k is the timestep index
	X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @constraint_abbrev_dubinscar(traj, traj_prev, SCPP)

  Xkp, Ukp, Xk = X[:,k-1], U[k-1], X[:,k]
  fpkp, Apkp, Bpkp, Xpkp, Upkp = get_f(k-1, model), get_A(k-1, model), get_B(k-1, model), Xp[:,k-1], Up[k-1]

  # Just Trapezoidal rule
  Uk = U[k]
  fpk, Apk, Bpk, Xpk, Upk = get_f(k, model), get_A(k, model), get_B(k, model), Xp[:,k], Up[k]
  return (Xkp-Xk) + 1/2*dtp.*(fpkp + Apkp*(Xkp-Xpkp) + Bpkp.*(Ukp-Upkp) +
                              fpk  + Apk *(Xk-Xpk)   + Bpk .*(Uk-Upk))
end

# Get current dynamics structures for a particular time step
get_f(k::Int, model::DubinsCar) = model.f[k]
get_A(k::Int, model::DubinsCar) = model.A[k]
get_B(k::Int, model::DubinsCar) = model.B

# Generate full dynamics structures for one time step
function f_dyn(x::Vector, u::Vector, robot::Robot, model::DubinsCar)
	x_dim = model.x_dim
	f = zeros(x_dim)
  update_f!(f, x, u, robot, model)
  return f
end

function update_f!(f, x::Vector, u::Vector, robot::Robot, model::DubinsCar)
	f[1] = model.v*cos(x[3])
	f[2] = model.v*sin(x[3])
	f[3] = model.k*u[1]
end

function A_dyn(x::Vector, robot::Robot, model::DubinsCar)
	x_dim = model.x_dim
	A = zeros(x_dim, x_dim)
  update_A!(A, x, robot, model)
  return A
end

function update_A!(A, x::Vector, robot::Robot, model::DubinsCar)
	A[1,3] = -model.v*sin(x[3])
	A[2,3] = model.v*cos(x[3])
end

function B_dyn(x::Vector, robot::Robot, model::DubinsCar)
  [0; 0; model.k]
end

## Constructing full list of constraint functions
function SCPConstraints(SCPP::SCPProblem{Car, DubinsCar, E}) where E
	model, goal_set = SCPP.PD.model, SCPP.PD.goal_set
	x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, SCPP.N, SCPP.tf_guess

	SCPC = SCPConstraints()

	## Dynamics constraints
  add_constraint_category!(SCPC.dynamics, dynamics_constraints, :array, 2:N)

	## Convex state equality constraints
	# Init and goal
	add_constraint_category!(SCPC.convex_state_eq, cse_init_constraints, :scalar, 0, 1:x_dim)
	for goal in values(inclusive(goal_set.goals, searchsortedfirst(goal_set.goals, tf_guess), searchsortedlast(goal_set.goals, tf_guess)))
		if typeof(goal.params) == PointGoal
			add_constraint_category!(SCPC.convex_state_boundary_condition_eq, csbce_goal_constraints, goal, :array)
		else
			add_constraint_category!(SCPC.convex_state_boundary_condition_ineq, csbci_goal_constraints, goal, :array)
		end
	end

	## Convex state inequality constraints
	# State bounds
	add_constraint_category!(SCPC.convex_state_ineq, csi_max_bound_constraints, :scalar, 1:N, 1:x_dim)
	add_constraint_category!(SCPC.convex_state_ineq, csi_min_bound_constraints, :scalar, 1:N, 1:x_dim)

	## Nonconvex state equality constraints
	nothing
	## Nonconvex state inequality constraints
	nothing
	## Nonconvex state equality constraints (convexified)
	nothing
	## Nonconvex state inequality constraints (convexified)
	nothing
	## Convex control equality constraints
	nothing

	## Convex control inequality constraints
	# Control bounds
	add_constraint_category!(SCPC.convex_control_ineq, cci_max_bound_constraints, :scalar, 1:N-1, 1)
	add_constraint_category!(SCPC.convex_control_ineq, cci_min_bound_constraints, :scalar, 1:N-1, 1)

	return SCPC
end

# TODO(ambyld): Replace f_dyn
function trust_region_ratio_gusto(traj, traj_prev::Trajectory, SCPP::SCPProblem{Car, DubinsCar, E}) where E
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @constraint_abbrev_dubinscar(traj, traj_prev, SCPP)
  fp, Ap = model.f, model.A
  num, den = 0, 0

  for k in 1:N-1
    linearized = fp[k] + Ap[k]*(X[:,k]-Xp[:,k])
    num += norm(f_dyn(X[:,k],U[:,k],robot,model) - linearized)
    den += norm(linearized)
  end

  return num/den
end

##################
# Shooting Method
##################
function get_dual_cvx(prob::Convex.Problem, SCPP::SCPProblem{Car, DubinsCar, E}, solver) where E
	if solver == "Mosek"
		return -MathProgBase.getdual(prob.model)[1:SCPP.PD.model.x_dim]
	else
		return []
	end
end

function get_dual_jump(SCPC::SCPConstraints, SCPP::SCPProblem{Car, DubinsCar, E}) where E
	init_constraint_category = SCPC.convex_state_eq[:cse_init_constraints][1]
  -JuMP.dual.([init_constraint_category.con_reference[0,(i,)] for i = init_constraint_category.ind_other[1]])
end

function shooting_ode!(xdot, x, SP::ShootingProblem{Car, DubinsCar, E}, t) where E
	X, p = x[1:3], x[4:6]
	U = get_control(X, p, SP)[1]
	xx, y, θ = X
  px, py, pθ = p
	robot, model = SP.PD.robot, SP.PD.model

	# State variables
	xdot[1] = model.v*cos(θ) 	 	# xdot
	xdot[2] = model.v*sin(θ)		# ydot
	xdot[3] = model.k*U 				# θdot

	# Dual variables
	xdot[4] = 0
	xdot[5] = 0
	xdot[6] = px*model.v*sin(θ) - py*model.v*cos(θ)
end

function get_control(x, p, SP::ShootingProblem{Car, DubinsCar, E}) where E
	model = SP.PD.model
	U = model.k/2*p[3,:]
end

#######
# JuMP
#######
function add_objective!(solver_model::Model, SCPV::SCPVariables, SCPP::SCPProblem{Astrobee3D{T}, AstrobeeSE3Manifold, E}) where {T,E}
  robot, model = SCPP.PD.robot, SCPP.PD.model
  x_dim, u_dim, N = model.x_dim, model.u_dim, SCPP.N

  U = SCPV.U
  N, dt = SCPP.N, SCPP.tf_guess/SCPP.N

  @NLobjective(solver_model, Min, sum(dt*U[i,k]^2 for i = 1:u_dim, k = 1:N-1))
end

ext: jl


text:
include("test_snr.jl")
include("test_radon.jl")


text:
# This is necessary or Memento will attempt to use `println` with one argument (unsupported by Syslog)
# This should be handled by optional dependency stuff in the future
"""
    emit(
        handler::DefaultHandler{F, O},
        rec::Record
    ) where {F<:Formatter, O<:Syslogs.Syslog}

Handles printing any records with any `Formatter` and a `Syslog` `IO` type.
"""
function Memento.emit(
    handler::DefaultHandler{F, O}, rec::Record
) where {F<:Formatter, O<:Syslogs.Syslog}

    println(handler.io, getlevel(rec), format(handler.fmt, rec))
    flush(handler.io)
end


text:
module BloomFilters
	export BloomFilter, add!

	if VERSION < v"1.0.0"
		import Base: contains
	else
		export contains
	end

	include("bloom-filter.jl")
end


text:
export Configuration

lazy_artifact(x) = @artifact_str(x)

const rootfs_lock = ReentrantLock()
const rootfs_cache = Dict()
function prepare_rootfs(distro="debian"; uid=1000, user="pkgeval", gid=1000, group="pkgeval", home="/home/$user")
    lock(rootfs_lock) do
        get!(rootfs_cache, (distro, uid, user, gid, group, home)) do
            base = lazy_artifact(distro)

            # a bare rootfs isn't usable out-of-the-box
            derived = mktempdir()
            cp(base, derived; force=true)

            # add a user and group
            chmod(joinpath(derived, "etc/passwd"), 0o644)
            open(joinpath(derived, "etc/passwd"), "a") do io
                println(io, "$user:x:$uid:$gid::$home:/bin/bash")
            end
            chmod(joinpath(derived, "etc/group"), 0o644)
            open(joinpath(derived, "etc/group"), "a") do io
                println(io, "$group:x:$gid:")
            end
            chmod(joinpath(derived, "etc/shadow"), 0o640)
            open(joinpath(derived, "etc/shadow"), "a") do io
                println(io, "$user:*:::::::")
            end

            # replace resolv.conf
            rm(joinpath(derived, "etc/resolv.conf"); force=true)
            write(joinpath(derived, "etc/resolv.conf"), read("/etc/resolv.conf"))

            return (path=derived, uid, user, gid, group, home)
        end
    end
end

"""
    run_sandboxed_julia(install::String, args=``; env=Dict(), mounts=Dict(),
                        wait=true, stdin=stdin, stdout=stdout, stderr=stderr,
                        install_dir="/opt/julia", kwargs...)

Run Julia inside of a sandbox, passing the given arguments `args` to it. The argument `wait`
determines if the process will be waited on. Streams can be connected using the `stdin`,
`stdout` and `sterr` arguments. Returns a `Process` object.

Further customization is possible using the `env` arg, to set environment variables, and the
`mounts` argument to mount additional directories. With `install_dir`, the directory where
Julia is installed can be chosen.
"""
function run_sandboxed_julia(install::String, args=``; wait=true,
                             mounts::Dict{String,String}=Dict{String,String}(),
                             kwargs...)
    config, cmd = runner_sandboxed_julia(install, args; kwargs...)

    # XXX: even when preferred_executor() returns UnprivilegedUserNamespacesExecutor,
    #      sometimes a stray sudo happens at run time? no idea how.
    exe_typ = UnprivilegedUserNamespacesExecutor
    exe = exe_typ()
    proc = Base.run(exe, config, cmd; wait)

    # TODO: introduce a --stats flag that has the sandbox trace and report on CPU, network, ... usage

    if wait
        cleanup(exe)
    else
        @async begin
            try
                Base.wait(proc)
                cleanup(exe)
            catch err
                @error "Unexpected error while cleaning up process" exception=(err, catch_backtrace())
            end
        end
    end

    return proc
end

# global Xvfb process for use by all containers
const xvfb_lock = ReentrantLock()
const xvfb_proc = Ref{Union{Base.Process,Nothing}}(nothing)

function runner_sandboxed_julia(install::String, args=``; install_dir="/opt/julia",
                                stdin=stdin, stdout=stdout, stderr=stderr,
                                env::Dict{String,String}=Dict{String,String}(),
                                mounts::Dict{String,String}=Dict{String,String}(),
                                xvfb::Bool=true, cpus::Vector{Int}=Int[])
    julia_path = installed_julia_dir(install)
    rootfs = prepare_rootfs()
    read_only_maps = Dict(
        "/"                                 => rootfs.path,
        install_dir                         => julia_path,
        "/usr/local/share/julia/registries" => registry_dir(),
    )

    artifacts_path = joinpath(storage_dir(), "artifacts")
    mkpath(artifacts_path)
    read_write_maps = merge(mounts, Dict(
        joinpath(rootfs.home, ".julia/artifacts")   => artifacts_path
    ))

    env = merge(env, Dict(
        # PkgEval detection
        "CI" => "true",
        "PKGEVAL" => "true",
        "JULIA_PKGEVAL" => "true",

        # use the provided registry
        # NOTE: putting a registry in a non-primary depot entry makes Pkg use it as-is,
        #       without needingb to set Pkg.UPDATED_REGISTRY_THIS_SESSION.
        "JULIA_DEPOT_PATH" => "::/usr/local/share/julia",

        # some essential env vars (since we don't run from a shell)
        "PATH" => "/usr/local/bin:/usr/local/sbin:/usr/bin:/usr/sbin:/bin:/sbin",
        "HOME" => rootfs.home,
    ))
    if haskey(ENV, "TERM")
        env["TERM"] = ENV["TERM"]
    end

    if xvfb
        lock(xvfb_lock) do
            if xvfb_proc[] === nothing || !process_running(xvfb_proc[])
                proc = Base.run(`Xvfb :1 -screen 0 1024x768x16`; wait=false)
                sleep(1)
                process_running(proc) || error("Could not start Xvfb")

                xvfb_proc[] === nothing && atexit() do
                    kill(xvfb_proc[])
                    wait(xvfb_proc[])
                end
                xvfb_proc[] = proc
            end
        end

        env["DISPLAY"] = ":1"
        read_write_maps["/tmp/.X11-unix"] = "/tmp/.X11-unix"
    end

    cmd = `$install_dir/bin/julia`

    # restrict resource usage
    if !isempty(cpus)
        cmd = `/usr/bin/taskset --cpu-list $(join(cpus, ',')) $cmd`
        env["JULIA_CPU_THREADS"] = string(length(cpus)) # JuliaLang/julia#35787
    end

    # NOTE: we use persist=true so that modifications to the rootfs are backed by
    #       actual storage on the host, and not just the (1G hard-coded) tmpfs,
    #       because some packages like to generate a lot of data during testing.

    config = SandboxConfig(read_only_maps, read_write_maps, env;
                           rootfs.uid, rootfs.gid, pwd=rootfs.home, persist=true,
                           stdin, stdout, stderr, verbose=isdebug(:sandbox))

    return config, `$cmd $args`
end

function process_children(pid)
    pids = Int[]
    for tid in readdir("/proc/$pid/task")
        children = read("/proc/$pid/task/$tid/children", String)
        append!(pids, parse.(Int, split(children)))
    end
    pids
end

function cpu_time(pid)
    stats = read("/proc/$pid/stat", String)
    m = match(r"^(\d+) \((.+)\) (.+)", stats)
    @assert m !== nothing
    fields = [[m.captures[1], m.captures[2]]; split(m.captures[3])]
    utime = parse(Int, fields[14])
    stime = parse(Int, fields[15])
    cutime = parse(Int, fields[16])
    cstime = parse(Int, fields[17])
    total_time = (utime + stime + cutime + cstime) / Sys.SC_CLK_TCK

    # cutime and cstime are only updated when the child exits,
    # so recursively scan all known children
    total_time += sum(cpu_time, process_children(pid); init=0.0)

    return total_time
end

"""
    run_sandboxed_test(install::String, pkg; do_depwarns=false,
                       log_limit=2^20, time_limit=60*60)

Run the unit tests for a single package `pkg` inside of a sandbox using a Julia installation
at `install`. If `do_depwarns` is `true`, deprecation warnings emitted while running the
package's tests will cause the tests to fail. Test will be forcibly interrupted after
`time_limit` seconds (defaults to 1h) or if the log becomes larger than `log_limit`
(defaults to 1MB).

A log for the tests is written to a version-specific directory in the PkgEval root
directory.

Refer to `run_sandboxed_julia`[@ref] for more possible `keyword arguments.
"""
function run_sandboxed_test(install::String, pkg; log_limit = 2^20 #= 1 MB =#,
                            time_limit = 60*60, do_depwarns=false,
                            kwargs...)
    @assert log_limit > 0

    # prepare for launching a container
    script = raw"""
        try
            using Dates
            print('#'^80, "\n# PkgEval set-up: $(now())\n#\n\n")

            using InteractiveUtils
            versioninfo()
            println()

            print("\n\n", '#'^80, "\n# Installation: $(now())\n#\n\n")

            using Pkg
            Pkg.add(ARGS[1])

            print("\n\n", '#'^80, "\n# Testing: $(now())\n#\n\n")

            Pkg.test(ARGS[1])

            println("\nPkgEval succeeded")
        catch err
            print("\nPkgEval failed: ")
            showerror(stdout, err)
            Base.show_backtrace(stdout, catch_backtrace())
            println()
        finally
            print("\n\n", '#'^80, "\n# PkgEval teardown: $(now())\n#\n\n")
        end"""
    cmd = do_depwarns ? `--depwarn=error` : ``
    cmd = `$cmd --eval 'eval(Meta.parse(read(stdin,String)))' $(pkg.name)`

    input = Pipe()
    output = Pipe()

    env = Dict(
        "JULIA_PKG_PRECOMPILE_AUTO" => "0",
        # package hacks
        "PYTHON" => "",
        "R_HOME" => "*"
    )
    if haskey(ENV, "JULIA_PKG_SERVER")
        env["JULIA_PKG_SERVER"] = ENV["JULIA_PKG_SERVER"]
    end

    proc = run_sandboxed_julia(install, cmd; env, wait=false,
                               stdout=output, stderr=output, stdin=input,
                               kwargs...)
    close(output.in)

    # pass the script over standard input to avoid exceeding max command line size,
    # and keep the process listing somewhat clean
    println(input, script)
    close(input)

    function stop()
        if process_running(proc)
            # FIXME: if we only kill proc, we sometimes only end up killing the sandbox.
            #        shouldn't the sandbox handle this, e.g., by creating a process group?
            function recursive_kill(proc, sig)
                parent_pid = getpid(proc)
                for pid in reverse([parent_pid; process_children(parent_pid)])
                    ccall(:uv_kill, Cint, (Cint, Cint), pid, Base.SIGKILL)
                end
                return
            end

            recursive_kill(proc, Base.SIGINT)
            terminator = Timer(5) do timer
                recursive_kill(proc, Base.SIGTERM)
            end
            killer = Timer(10) do timer
                recursive_kill(proc, Base.SIGKILL)
            end
            wait(proc)
            close(terminator)
            close(killer)
        end
        close(output)
    end

    status = nothing
    reason = missing

    # kill on timeout
    timeout_monitor = Timer(time_limit) do timer
        process_running(proc) || return
        status = :kill
        reason = :time_limit
        stop()
    end

    # kill on inactivity (less than 1 second of CPU usage every minute)
    previous_cpu_time = nothing
    inactivity_monitor = Timer(60; interval=30) do timer
        process_running(proc) || return
        pid = getpid(proc)
        current_cpu_time = cpu_time(pid)
        if current_cpu_time > 0 && previous_cpu_time !== nothing
            cpu_time_diff = current_cpu_time - previous_cpu_time
            if 0 <= cpu_time_diff < 1
                status = :kill
                reason = :inactivity
                stop()
            end
        end
        previous_cpu_time = current_cpu_time
    end

    # collect output
    log_monitor = @async begin
        io = IOBuffer()
        while !eof(output)
            print(io, readline(output; keep=true))

            # kill on too-large logs
            if io.size > log_limit
                process_running(proc) || break
                status = :kill
                reason = :log_limit
                stop()
                break
            end
        end
        return String(take!(io))
    end

    wait(proc)
    close(timeout_monitor)
    close(inactivity_monitor)
    log = fetch(log_monitor)

    if sizeof(log) > log_limit
        # even though the monitor above should have limited the log size,
        # a single line may still have exceeded the limit, so make sure we truncate.
        ind = prevind(log, log_limit)
        log = log[1:ind]
    end

    # pick up the installed package version from the log
    version_match = match(Regex("Installed $(pkg.name) .+ v(.+)"), log)
    version = if version_match !== nothing
        try
            VersionNumber(version_match.captures[1])
        catch
            @error "Could not parse installed package version number '$(version_match.captures[1])'"
            v"0"
        end
    else
        missing
    end

    # try to figure out the failure reason
    if status === nothing
        if occursin("PkgEval succeeded", log)
            status = :ok
        else
            status = :fail

            # figure out a more accurate failure reason from the log
            reason = if occursin("Unsatisfiable requirements detected for package", log)
                # NOTE: might be the package itself, or one of its dependencies
                :unsatisfiable
            elseif occursin("Package $(pkg.name) did not provide a `test/runtests.jl` file", log)
                :untestable
            elseif occursin("cannot open shared object file: No such file or directory", log)
                :binary_dependency
            elseif occursin(r"Package .+ does not have .+ in its dependencies", log)
                :missing_dependency
            elseif occursin(r"Package .+ not found in current path", log)
                :missing_package
            elseif occursin("Some tests did not pass", log) || occursin("Test Failed", log)
                :test_failures
            elseif occursin("ERROR: LoadError: syntax", log)
                :syntax
            elseif occursin("GC error (probable corruption)", log)
                :gc_corruption
            elseif occursin("signal (11): Segmentation fault", log)
                :segfault
            elseif occursin("signal (6): Abort", log)
                :abort
            elseif occursin("Unreachable reached", log)
                :unreachable
            else
                :unknown
            end
        end
    end

    return version, status, reason, log
end

Base.@kwdef struct Configuration
    julia::VersionNumber = Base.VERSION
    # TODO: depwarn, checkbounds, etc
    # TODO: also move buildflags here?
end

# behave as a scalar in broadcast expressions
Base.broadcastable(x::Configuration) = Ref(x)

function run(configs::Vector{Configuration}, pkgs::Vector;
             ninstances::Integer=Sys.CPU_THREADS, retries::Integer=2, kwargs...)
    # here we deal with managing execution: spawning workers, output, result I/O, etc

    # Julia installation
    instantiated_configs = Dict{Configuration,String}()
    for config in configs
        install = prepare_julia(config.julia)
        instantiated_configs[config] = install
    end

    jobs = vec(collect(Iterators.product(instantiated_configs, pkgs)))

    # use a random test order to (hopefully) get a more reasonable ETA
    shuffle!(jobs)

    njobs = length(jobs)
    ninstances = min(njobs, ninstances)
    running = Vector{Any}(nothing, ninstances)
    times = DateTime[now() for i = 1:ninstances]
    all_workers = Task[]

    done = false
    did_signal_workers = false
    function stop_work()
        if !done
            done = true
            if !did_signal_workers
                for (i, task) in enumerate(all_workers)
                    task == current_task() && continue
                    Base.istaskdone(task) && continue
                    try; schedule(task, InterruptException(); error=true); catch; end
                    running[i] = nothing
                end
                did_signal_workers = true
            end
        end
    end

    start = now()
    io = IOContext(IOBuffer(), :color=>true)
    p = Progress(njobs; barlen=50, color=:normal)
    function update_output()
        # known statuses
        o = count(==(:ok),      result[!, :status])
        s = count(==(:skip),    result[!, :status])
        f = count(==(:fail),    result[!, :status])
        k = count(==(:kill),    result[!, :status])
        # remaining
        x = njobs - nrow(result)

        function runtimestr(start)
            time = Dates.canonicalize(Dates.CompoundPeriod(now() - start))
            if isempty(time.periods) || first(time.periods) isa Millisecond
                "just started"
            else
                # be coarse in the name of brevity
                string(first(time.periods))
            end
        end

        if !isinteractive()
            println("$x combinations to test ($o succeeded, $f failed, $k killed, $s skipped, $(runtimestr(start)))")
            sleep(10)
        else
            print(io, "Success: ")
            printstyled(io, o; color = :green)
            print(io, "\tFailed: ")
            printstyled(io, f; color = Base.error_color())
            print(io, "\tKilled: ")
            printstyled(io, k; color = Base.error_color())
            print(io, "\tSkipped: ")
            printstyled(io, s; color = Base.warn_color())
            println(io, "\tRemaining: ", x)
            for i = 1:ninstances
                job = running[i]
                str = if job === nothing
                    " #$i: -------"
                else
                    config, pkg = job
                    " #$i: $(pkg.name) @ $(config.julia) ($(runtimestr(times[i])))"
                end
                if i%2 == 1 && i < ninstances
                    print(io, rpad(str, 50))
                else
                    println(io, str)
                end
            end
            print(String(take!(io.io)))
            p.tlast = 0
            update!(p, nrow(result))
            sleep(1)
            CSI = "\e["
            print(io, "$(CSI)$(ceil(Int, ninstances/2)+1)A$(CSI)1G$(CSI)0J")
        end
    end

    # NOTE: we expand the Configuration into separate columns
    result = DataFrame(julia = VersionNumber[],
                       name = String[],
                       uuid = UUID[],
                       version = Union{Missing,VersionNumber}[],
                       status = Symbol[],
                       reason = Union{Missing,Symbol}[],
                       duration = Float64[],
                       log = Union{Missing,String}[])

    # Printer
    @async begin
        try
            while (!isempty(jobs) || !all(==(nothing), running)) && !done
                update_output()
            end
            stop_work()
        catch e
            stop_work()
            isa(e, InterruptException) || rethrow(e)
        end
    end

    # Workers
    try @sync begin
        for i = 1:ninstances
            push!(all_workers, @async begin
                try
                    while !isempty(jobs) && !done
                        (config, install), pkg = pop!(jobs)
                        times[i] = now()
                        running[i] = (config, pkg)

                        # can we even test this package?
                        julia_supported = Dict{VersionNumber,Bool}()
                        ctx = Pkg.Types.Context()
                        pkg_version_info = Pkg.Operations.load_versions(ctx, pkg.path)
                        pkg_versions = sort!(collect(keys(pkg_version_info)))
                        pkg_compat =
                            Pkg.Operations.load_package_data(ctx, Pkg.Types.VersionSpec,
                                                             joinpath(pkg.path,
                                                                      "Compat.toml"),
                                                             pkg_versions)
                        for (pkg_version, bounds) in pkg_compat
                            if haskey(bounds, "julia")
                                julia_supported[pkg_version] =
                                    config.julia ∈ bounds["julia"]
                            end
                        end
                        if length(julia_supported) != length(pkg_version_info)
                            # not all versions have a bound for Julia,
                            # so we need to be conservative
                            supported = true
                        else
                            supported = any(values(julia_supported))
                        end
                        if !supported
                            push!(result, [config.julia,
                                           pkg.name, pkg.uuid, missing,
                                           :skip, :unsupported, 0, missing])
                            running[i] = nothing
                            continue
                        elseif pkg.name in skip_lists[pkg.registry]
                            push!(result, [config.julia,
                                           pkg.name, pkg.uuid, missing,
                                           :skip, :explicit, 0, missing])
                            running[i] = nothing
                            continue
                        elseif endswith(pkg.name, "_jll")
                            push!(result, [config.julia,
                                           pkg.name, pkg.uuid, missing,
                                           :skip, :jll, 0, missing])
                            running[i] = nothing
                            continue
                        end

                        # perform an initial run
                        pkg_version, status, reason, log =
                            run_sandboxed_test(install, pkg; cpus=[i-1], kwargs...)

                        # certain packages are known to have flaky tests; retry them
                        for j in 1:retries
                            if status == :fail && reason == :test_failures &&
                               pkg.name in retry_lists[pkg.registry]
                                times[i] = now()
                                pkg_version, status, reason, log =
                                    run_sandboxed_test(install, pkg; cpus=[i-1], kwargs...)
                            end
                        end

                        duration = (now()-times[i]) / Millisecond(1000)
                        push!(result, [config.julia,
                                       pkg.name, pkg.uuid, pkg_version,
                                       status, reason, duration, log])
                        running[i] = nothing
                    end
                catch e
                    stop_work()
                    isa(e, InterruptException) || rethrow(e)
                end
            end)
        end
    end
    catch e
        isa(e, InterruptException) || rethrow(e)
    finally
        stop_work()
        println()

        # clean-up
        for (config, install) in instantiated_configs
            rm(install; recursive=true)
        end
    end

    return result
end

"""
    run(configs::Vector{Configuration}=[Configuration()],
        pkg_names::Vector{String}=[]]; registry=General, update_registry=true, kwargs...)

Run all tests for all packages in the registry `registry`, or only for the packages as
identified by their name in `pkgnames`, using the configurations from `configs`.
The registry is first updated if `update_registry` is set to true.

Refer to `run_sandboxed_test`[@ref] and `run_sandboxed_julia`[@ref] for more possible
keyword arguments.
"""
function run(configs::Vector{Configuration}=[Configuration()],
             pkg_names::Vector{String}=String[];
             registry::String=DEFAULT_REGISTRY, update_registry::Bool=true, kwargs...)
    prepare_registry(registry; update=update_registry)
    pkgs = read_pkgs(pkg_names)

    run(configs, pkgs; kwargs...)
end

# for backwards compatibility
run(julia_versions::Vector{VersionNumber}, args...; kwargs...) =
    run([Configuration(julia=julia_version) for julia_version in julia_versions], args...;
        kwargs...)
prepare_runner() = return


text:
"""
    NonAct()
    (l::NonAct)(x) = x

Identity activation does nothing
"""
struct NonAct <: Activation end
@inline (l::NonAct)(x) = x

"""
    ReLU()
    (l::ReLU)(x) = max(0,x)

Rectified Linear Unit function.
"""
struct ReLU <: Activation end
@inline (l::ReLU)(x) = relu.(x)

"""
    Sigm()
    (l::Sigm)(x) = sigm(x)

Sigmoid function
"""
struct Sigm <: Activation end
@inline (l::Sigm)(x) = sigm.(x)

"""
    Tanh()
    (l::Tanh)(x) = tanh(x)

Tangent hyperbolic function
"""
struct Tanh <: Activation end
@inline (l::Tanh)(x) = tanh.(x)

"""
    ELU()
    (l::ELU)(x) = elu(x) -> Computes x < 0 ? exp(x) - 1 : x

Exponential Linear Unit nonlineariy.
"""
struct ELU <: Activation end
@inline (l::ELU)(x) = elu.(x)

"""
    LeakyReLU(α=0.2)
    (l::LeakyReLU)(x) -> Computes x < 0 ? α*x : x
"""
mutable struct LeakyReLU <: Activation
    α::AbstractFloat
    LeakyReLU(alpha::AbstractFloat=0.2) = new(alpha)
end
@inline (l::LeakyReLU)(x) = relu.(x) .+ l.α*min.(0,x) # relu.(x) .+ l.α*relu.(-x) ?

#FIXME: x*x*x -> x.^3 after Knet new version
mutable struct GeLU <: Activation
end
"""
    gelu(x) = 0.5x*(1 + tanh(√(2/π)*(x + 0.044715x^3)))
[Gaussian Error Linear Unit](https://arxiv.org/pdf/1606.08415.pdf)
activation function.
"""
function (::GeLU)(x)
    T = eltype(x)
    λ = T(√(2/π))
    α = T(0.044715)
    h = T(0.5)
    return h .* x .* (one(T) .+ tanh.(λ .* (x .+ α * x .* x .* x)))
end
"""
    Dropout(p=0)

Dropout Layer. `p` is the droput probability.
"""
mutable struct Dropout <: Activation
    p::Real
end
Dropout(;p=0) = Dropout(p)
@inline (l::Dropout)(x; enable=true) = enable ? dropout(x,l.p) : x

"""
    LogSoftMax(dims=:)
    (l::LogSoftMax)(x)

Treat entries in x as as unnormalized log probabilities and return normalized log probabilities.

dims is an optional argument, if not specified the normalization is over the whole x, otherwise the normalization is performed over the given dimensions. In
particular, if x is a matrix, dims=1 normalizes columns of x, dims=2 normalizes rows of x.
"""
struct LogSoftMax <: Activation
    dims::Union{Integer,Colon}
end
LogSoftMax(;dims=:) = LogSoftMax(dims)
@inline (l::LogSoftMax)(x) = logp(x;dims=l.dims)

"""
    SoftMax(dims=:)
    (l::SoftMax)(x)

Treat entries in x as as unnormalized scores and return softmax probabilities.

dims is an optional argument, if not specified the normalization is over the whole x, otherwise the normalization is performed over the given dimensions. In
particular, if x is a matrix, dims=1 normalizes columns of x, dims=2 normalizes rows of x.
"""
struct SoftMax <: Activation
    dims::Union{Integer,Colon}
end
SoftMax(;dims=:) = SoftMax(dims)
@inline (l::SoftMax)(x) = softmax(x;dims=l.dims)

"""
    LogSumExp(dims=:)
    (l::LogSumExp)(x)

  Compute log(sum(exp(x);dims)) in a numerically stable manner.

  dims is an optional argument, if not specified the summation is over the whole x, otherwise the summation is performed over the given dimensions. In particular if x
  is a matrix, dims=1 sums the columns of x, dims=2 sums the rows of x.
"""
struct LogSumExp <: Activation
    dims::Union{Integer,Colon}
end
LogSumExp(;dims=:) = LogSumExp(dims)
@inline (l::LogSumExp)(x) = logsumexp(x;dims=l.dims)

"""
    normalise(x; dims=1)
    Normalises x to mean 0 and standard deviation 1, across the dimensions given by dims. Defaults to normalising over columns.
"""
function normalise(x; dims=1)
  μ′ = mean(x, dims = dims)
  σ′ = std(x, dims = dims, mean = μ′, corrected=false)
  return (x .- μ′) ./ σ′
end


text:
for (T, TP) in [(Cheb, ChebParams), (Lin, LinParams), (Spline, SplineParams)]
    @eval _param(::$(T)) = $TP
    @eval _param(::Type{$(T)}) = $TP
end

# params of same type are equal if all their fields are equal
==(p1::T, p2::T) where {T<:BasisParams} =
    all(map(nm->getfield(p1, nm) == getfield(p2, nm), fieldnames(T)))::Bool

# Bases of different dimension can't be equal
==(::T1, ::T2) where {T1<:BasisParams,T2<:BasisParams} = false

# ---------- #
# Basis Type #
# ---------- #

struct Basis{N,TP<:Tuple}
    params::TP     # params to construct basis
end

Base.min(b::Basis) = min.(b.params)
Base.max(b::Basis) = max.(b.params)
Base.ndims(::Basis{N}) where {N} = N
Base.ndims(::Type{Basis{N,TP}}) where {N,TP} = N

_get_TP(::Basis{N,TP}) where {N,TP} = TP
_get_TP(::Type{Basis{N,TP}}) where {N,TP} = TP

function Base.show(io::IO, b::Basis{N}) where N
    m = """
    $N dimensional Basis on the hypercube formed by $(min(b)) × $(max(b)).
    Basis families are $(join(string.(family_name.(b.params)), " × "))
    """
    print(io, m)
end

Basis(params::BasisParams...) = _Basis(params)
Basis(params::Tuple) = _Basis(params)

# hack to make method above type stable
@generated function _Basis(params)
    N = length(params.parameters)
    quote
        Basis{$N,$params}(params)
    end
end

# combining basis -- fundef-esque method
Basis(bs::Basis...) = _Basis2(bs)

# Another hack to make the above type stable
@generated function _Basis2(bs)
    N = sum(ndims, bs.parameters)

    # tup_of_tups will be a tuple where each element is a tuple of types
    # we want to concatenate them and end up with with a single tuple of types
    # to do that we first put them in a vector, then splat that vector into a
    # Tuple. Note that we don't use `tuple` because that will create a tuple
    # obejct, not the `Tuple` type.
    tup_of_tups = map(_get_TP, bs.parameters)
    basis_types = []
    for x in tup_of_tups
        push!(basis_types, x.parameters...)
    end
    TP = Tuple{basis_types...}
    quote
        new_params = []
        for x in bs
            push!(new_params, x.params...)
        end
        Basis{$N,$TP}(tuple(new_params...))
    end
end

# fundefn type method
Basis(bt::BasisFamily, n::Int, a, b) = Basis(_param(bt)(n, a, b))

Basis(::Type{T}, n::Int, a, b) where {T<:BasisFamily} = Basis(T(), n, a, b)

Basis(bt::T, n::Vector, a::Vector, b::Vector) where {T<:BasisFamily} =
    Basis(map(_param(T), n, a, b)...)

# special method for Spline that adds `k` argument
Basis(::Spline, n::Int, a, b, k) = Basis(SplineParams(n, a, b, k))
Basis(::Spline, n::Vector, a::Vector, b::Vector, k::Vector=ones(Int, length(n))) =
    Basis(map(SplineParams, n, a, b, k)...)::Basis{length(n)}

# ----------------- #
# Basis API methods #
# ----------------- #

# separating Basis -- just re construct it from the nth set of params
function Base.getindex(basis::Basis{N}, n::Int) where N
    n < 0 || n > N && error("n must be between 1 and $N")
    Basis(basis.params[n])::Basis{1}
end

_all_sparse(b::Basis{N,TP}) where {N,TP} = all(issparse, TP.parameters)

# other AbstractArray like methods for Basis
Base.length(b::Basis) = prod(length, b.params)
Base.size(b::Basis, i::Int) = length(b[i])  # uses method on previous line
Base.size(b::Basis{N}) where {N} = map(length, b.params)

# Bases of different dimension can't be equal
==(::Basis{N}, ::Basis{M}) where {N,M} = false

# basis are equal if all fields of the basis are equal
==(b1::Basis{N}, b2::Basis{N}) where {N} =
    all(map(nm->getfield(b1, nm) == getfield(b2, nm), fieldnames(Basis)))::Bool

function nodes(b::Basis{1})
    x = nodes(b.params[1])
    (x, (x,))
end

function nodes(b::Basis)  # funnode method
    xcoord = nodes.(b.params)
    x = gridmake(xcoord...)
    return x, xcoord
end

@generated function bmat_type(::Type{TO}, bm::Basis{N,TP}, x=1.0) where {N,TP,TO}
    if N == 1
        out = bmat_type(TO, TP.parameters[1], x)
    else
        out = bmat_type(TO, TP.parameters[1], x)
        for this_TP in TP.parameters[2:end]
            this_out = bmat_type(TO, this_TP, x)
            if this_out != out
                out = AbstractMatrix{promote_type(eltype(out), eltype(this_out))}
            end
        end
    end
    return :($out)
end

bmat_type(b::Basis) = bmat_type(Nothing, b)
bmat_type(b::Basis, x) = bmat_type(Nothing, b, x)


text:
"""
    MSE
"""
function eval_mse_kernel!(eval::CuDeviceVector{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(y)
        @inbounds eval[i] = (p[1,i] - y[i])^2
    end
    return nothing
end

function eval_metric(::Val{:mse}, eval::AbstractVector{T}, p::AbstractMatrix{T}, y::AbstractVector{T}, α; MAX_THREADS=1024) where T <: AbstractFloat
    threads = min(MAX_THREADS, length(y))
    blocks = ceil(Int, length(y) / threads)
    @cuda blocks = blocks threads = threads eval_mse_kernel!(eval, p, y)
    CUDA.synchronize()
    return mean(eval)
end

"""
    Logloss
"""
function eval_logloss_kernel!(eval::CuDeviceVector{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(y)
        @inbounds pred = sigmoid(p[1,i])
        @inbounds eval[i] = -y[i] * log(pred) + (y[i] - 1) * log(1 - pred)
    end
    return nothing
end

function eval_metric(::Val{:logloss}, eval::AbstractVector{T}, p::AbstractMatrix{T}, y::AbstractVector{T}, α; MAX_THREADS=1024) where T <: AbstractFloat
    threads = min(MAX_THREADS, length(y))
    blocks = ceil(Int, length(y) / threads)
    @cuda blocks = blocks threads = threads eval_logloss_kernel!(eval, p, y)
    CUDA.synchronize()
    return mean(eval)
end

"""
    Gaussian
"""
function eval_gaussian_kernel!(eval::CuDeviceVector{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if i <= length(y)
        @inbounds eval[i] = p[2,i] + (y[i] - p[1,i])^2 / (2 * exp(2 * p[2,i]))
    end
    return nothing
end

function eval_metric(::Val{:gaussian}, eval::AbstractVector{T}, p::AbstractMatrix{T}, y::AbstractVector{T}, α; MAX_THREADS=1024) where T <: AbstractFloat
    threads = min(MAX_THREADS, length(y))
    blocks = ceil(Int, length(y) / threads)
    @cuda blocks = blocks threads = threads eval_gaussian_kernel!(eval, p, y)
    CUDA.synchronize()
    return mean(eval)
end

"""
    Poisson
        Unsupported factorial on CUDA at the moment
"""
# function eval_pois_kernel!(eval::CuDeviceVector{T}, p::CuDeviceMatrix{T}, y::CuDeviceVector{T}) where {T <: AbstractFloat}
#     i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
#     if i <= length(y)
#         @inbounds eval[i] = exp(p[1,i]) * (1 - y[i]) + log(factorial(y[i]))
#     end
#     return nothing
# end

# function eval_metric(::Val{:poisson}, eval::AbstractVector{T}, p::AbstractMatrix{T}, y::AbstractVector{T}, α; MAX_THREADS=1024) where T <: AbstractFloat
#     threads = min(MAX_THREADS, length(y))
#     blocks = ceil(Int, length(y) / threads)
#     @cuda blocks = blocks threads = threads eval_pois_kernel!(eval, p, y)
#     CUDA.synchronize()
#     return mean(eval)
# end

# function eval_metric(::Val{:mlogloss}, pred::AbstractMatrix{T}, y::Vector{S}, α=0.0) where {T <: AbstractFloat, S <: Integer}
#     eval = zero(T)
#     L = length(y)
#     K = size(pred,2)
#     # pred = pred - maximum.(pred)
#     @inbounds for i in 1:L
#         pred[i] = pred[i] .- maximum(pred[i])
#         soft_pred = exp.(pred[i]) / sum(exp.(pred[i]))
#         eval -= log(soft_pred[y[i]])
#     end
#     eval /= length(y)
#     return eval
# end

# function eval_metric(::Val{:quantile}, pred::AbstractMatrix{T}, y::AbstractVector{T}, α=0.0) where T <: AbstractFloat
#     eval = zero(T)
#     for i in 1:length(y)
#         eval += α * max(y[i] - pred[i,1], zero(T)) + (1-α) * max(pred[i,1] - y[i], zero(T))
#     end
#     eval /= length(y)
#     return eval
# end


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ 56a604d4-ec32-404f-91bc-86fa24b9bb47
using Pkg, DrWatson

# ╔═╡ 286f9602-6773-4fe0-ac75-36b19d8333f0
begin
	@quickactivate "NumericalMethodsforEngineers"
	using Symbolics, SymbolicUtils, Latexify, LaTeXStrings, ModelingToolkit
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ 987efa36-893d-45cc-88b5-234b80071ab2
using NamedTupleTools

# ╔═╡ 278c20ed-066d-436d-836a-10c34661c3eb
md"
!!! tip
	Output from `versioninfo()` is shown in terminal.
"

# ╔═╡ 8e5d2539-81c6-4396-81ff-a6908b582096
versioninfo()

# ╔═╡ f419c6fc-7b77-4ec6-80c0-8672b064c7f4
begin
	@variables t x y
	D = Differential(t)
end;

# ╔═╡ 303f5542-8e24-4994-a701-2bb664560699
md"
!!! note
	Helper functions such as `LATEXIFY()` and `L()` are defined in pluto-helpers.jl
"

# ╔═╡ 0eb5fb16-765e-4bfc-9300-e7a867037eb5
begin
	a = [x^2+y 0 2x;
		 0 0 2y;
		 y^2 + x 0 0]
	LATEXIFY("a", a)
end

# ╔═╡ 62588b63-213e-44b4-b3a0-eec4b1c930c3
L("a", a)

# ╔═╡ 2292c3af-6ca4-415f-b637-610dbbe42cab
begin
	z = t*t^4 + 3t^5
	D(z) # symbolic representation of derivative(t*t^4 + 15t^4, t)
	L(expand_derivatives(D(z))) # 10t^4
end

# ╔═╡ c899af71-1ae1-4b17-ba3a-d3ba989e49b6
L(z)

# ╔═╡ 5e617b9a-eedc-4a0a-82df-57dad2b91b3e
md"b0 = $(latexify(z))"

# ╔═╡ 340e4b31-4de7-40e3-816e-6e705df711de
begin
	b0 = Symbolics.jacobian([x + x*y, x^2 + y],[x, y])
	md"b0 = $(latexify(b0))"
end

# ╔═╡ 0d87482f-7a1a-495f-94fb-a6a625534801
begin
	b = [
		t^2 + t + t^2  2t + 4t
		x + y + y + 2t  x^2 - x^2 + y^2
	]
	b1 = simplify.(b)
	Text(sprint(show, "text/plain", b1))
end

# ╔═╡ ede004a7-e53c-4102-a4dd-6529dcd8f855
md"b1 = $(latexify(b1))"

# ╔═╡ a9376da1-b58a-4134-ab77-3db8ef42695e
begin
	b2 = simplify.(substitute.(b, (Dict(x => y^2),)))
	md"b2 = $(latexify(b2))"
end

# ╔═╡ 65237ef6-7c92-4687-96cb-766016f5ed73
md"
!!! note
	Element b3[2,1] == 16
"

# ╔═╡ c0df2952-f0f0-414a-92a3-3b1c6c622206
begin
	#2×2 Matrix{Num}:
	# 36.0  24.0
	# 16.0   9.0

	b3 = substitute.(b, (Dict(x => 2.0, y => 3.0, t => 4.0),))
	LATEXIFY("b3", b3)
end

# ╔═╡ e4998080-32ef-4095-bad0-92cc27bf5170
LATEXIFY(b3)

# ╔═╡ 9b7a5d1d-1a3e-40bc-86c2-6cf7cedc57a0
begin
	#2×2 Matrix{Num}:
	#   t + 2(t^2)   6t
	# 2t + 2y  y^2

	b4 = simplify.(substitute.(b, (Dict(x => 2y),)))
	Text(sprint(show, "text/plain", b4))
end

# ╔═╡ 412ef249-a7ac-42ef-b30e-4e77cddcd9a6
begin
	b5 = substitute.(b4, (Dict(x => 2.0, y => 3.0, t => 4.0),))
	Text(sprint(show, "text/plain", b5))
end

# ╔═╡ 1e48e58c-509b-42fd-b7e8-a943ff48d7f9
struct MyStruct
           tally::Int
           team::String
       end

# ╔═╡ 170bf340-b629-4cce-831f-570b1f7f4494
mystruct = MyStruct(5, "hometeam")

# ╔═╡ 7f3c0230-0619-498f-85de-9ab5982f72c7
mynamedtuple = ntfromstruct(mystruct)

# ╔═╡ 93ff9a77-8419-465a-8f7b-b565179364bf
ntstruct = structfromnt(MyStruct, mynamedtuple)

# ╔═╡ 2aeebb4c-379e-4876-9d1e-99e7bc8622f2
mystruct == ntstruct

# ╔═╡ Cell order:
# ╠═56a604d4-ec32-404f-91bc-86fa24b9bb47
# ╠═286f9602-6773-4fe0-ac75-36b19d8333f0
# ╟─278c20ed-066d-436d-836a-10c34661c3eb
# ╠═8e5d2539-81c6-4396-81ff-a6908b582096
# ╠═f419c6fc-7b77-4ec6-80c0-8672b064c7f4
# ╟─303f5542-8e24-4994-a701-2bb664560699
# ╠═0eb5fb16-765e-4bfc-9300-e7a867037eb5
# ╠═62588b63-213e-44b4-b3a0-eec4b1c930c3
# ╠═2292c3af-6ca4-415f-b637-610dbbe42cab
# ╠═c899af71-1ae1-4b17-ba3a-d3ba989e49b6
# ╠═5e617b9a-eedc-4a0a-82df-57dad2b91b3e
# ╠═340e4b31-4de7-40e3-816e-6e705df711de
# ╠═0d87482f-7a1a-495f-94fb-a6a625534801
# ╠═ede004a7-e53c-4102-a4dd-6529dcd8f855
# ╠═a9376da1-b58a-4134-ab77-3db8ef42695e
# ╠═e4998080-32ef-4095-bad0-92cc27bf5170
# ╟─65237ef6-7c92-4687-96cb-766016f5ed73
# ╠═c0df2952-f0f0-414a-92a3-3b1c6c622206
# ╠═9b7a5d1d-1a3e-40bc-86c2-6cf7cedc57a0
# ╠═412ef249-a7ac-42ef-b30e-4e77cddcd9a6
# ╠═987efa36-893d-45cc-88b5-234b80071ab2
# ╠═1e48e58c-509b-42fd-b7e8-a943ff48d7f9
# ╠═170bf340-b629-4cce-831f-570b1f7f4494
# ╠═7f3c0230-0619-498f-85de-9ab5982f72c7
# ╠═93ff9a77-8419-465a-8f7b-b565179364bf
# ╠═2aeebb4c-379e-4876-9d1e-99e7bc8622f2


text:
# ===============================
# Test minimal (default) options
# ===============================
using Reachability

# linear ODE: x' = Ax
A = [ 0.0509836  0.168159  0.95246   0.33644
      0.42377    0.67972   0.129232  0.126662
      0.518654   0.981313  0.489854  0.588326
      0.38318    0.616014  0.518412  0.778765]

# initial set
X0 = BallInf(ones(4), 0.1)

# default options (computes all variables)
s = solve(InitialValueProblem(LinearContinuousSystem(A), X0), :T=>0.1)

# two variables and custom partition
s = solve(ContinuousSystem(A, X0), :T=>0.1, :partition=>[1:2, 3:4], :vars=>[1,3])

# the default partition is used. uses Interval for set_type_init and set_type_iter
# but Hyperrectangle for set_type_proj
s = solve(ContinuousSystem(A, X0), :T=>0.1, :vars=>[1,3])

# ===============================
# Test projection
# ===============================
s = solve(ContinuousSystem(A, X0), :T=>0.1, :partition=>[1:2, 3:4],
          :vars=>[1,3], :project_reachset=>false);
ps = project(s);

# ===============
# Test check mode
# ===============

A = [ 0.0509836  0.168159  0.95246   0.33644
      0.42377    0.67972   0.129232  0.126662
      0.518654   0.981313  0.489854  0.588326
      0.38318    0.616014  0.518412  0.778765]

# check that x1 + x2 <= 2 doesn't hold
s = solve(ContinuousSystem(A, X0), :T=>0.1, :mode=>"check",
          :partition=>[1:2, 3:4], :vars=>[1,2,3],
          :property=>LinearConstraintProperty([1., 1., 0., 0.], 2.))
@test s.violation == 1

# check that x1 - x2 <= 2 holds
s = solve(ContinuousSystem(A, X0), :T=>0.1, :mode=>"check",
          :partition=>[1:2, 3:4], :vars=>[1,2,3],
          :property=>LinearConstraintProperty([1., -1., 0., 0.], 2.))
@test s.violation == -1

# ===============================
# Test reachability options
# ===============================
s = solve(ContinuousSystem(A, X0), :T=>0.1, :partition=>[1:2, 3:4],
          :vars=>[1,3], :lazy_sih=>true);
s = solve(ContinuousSystem(A, X0), :T=>0.1, :partition=>[1:2, 3:4],
          :vars=>[1,3], :lazy_sih=>false);

s = solve(ContinuousSystem(sparse(A), X0), :T=>0.1, :partition=>[1:2, 3:4],
          :vars=>[1,3], :lazy_expm=>true, :assume_sparse=>false);
s = solve(ContinuousSystem(sparse(A), X0), :T=>0.1, :partition=>[1:2, 3:4],
          :vars=>[1,3], :lazy_expm=>true, :assume_sparse=>true);

s = solve(ContinuousSystem(sparse(A), X0), :T=>0.1, :partition=>[[i] for i in 1:4],
          :set_type=>Interval, :vars=>[1,3], :lazy_expm=>true);

# ===============================
# System with an odd dimension
# ===============================
A = randn(5, 5); X0 = BallInf(ones(5), 0.1)
system = ContinuousSystem(A, X0)
options = Options(Dict(:T=>0.1, :partition=>[1:2, 3:4, [5]], :block_types=>Dict(HPolygon=>[1:2, 3:4], Interval=>[[5]]), :vars=>[1,3]))
s = solve(system, options)

# ===============================
# System with an odd dimension
# ===============================
A = [1 2 1.; 0 0. 1; -2 1 4]
X0 = BallInf(ones(3), 0.1)
system = ContinuousSystem(A, X0)
s = solve(system, :T=>0.1, :vars=>1:3, :partition=>[1:2, 3:3])


text:
# Numerically stable expit
@inline function expit(t::T) where T
    if t >= 0.0
        return one(T) / (one(T) + exp(-t))
    else
        e_t = exp(t)
        return e_t / (one(T) + e_t)
    end
end

# A numerically robust function to evaluate -log(1 + exp(-t))
# See: http://fa.bianp.net/blog/2019/evaluate_logistic/
@inline function log1pexp(t::T)::T where T
    if t < -33.3
        return t
    elseif t <= -18.0
        return t - exp(t)
    elseif t <= 37.0
        return -log1p(exp(-t))
    else
        return -exp(-t)
    end
end

"Compute Fenchel transform of f(x) = log(1 + exp(-x))."
function logloss(λ::T) where T <: Real
    # TODO: numerically robust version
    if λ == 0.0 || λ == -1.0
        return 0.0
    elseif -1.0 < λ < 0.0
        return (1.0 + λ) * log(1.0 + λ) - λ * log(-λ)
    else
        return Inf
    end
end

abstract type AbstractScaler end

"Scale dataset before optimization."
scale!(::AbstractScaler, X::AbstractArray{T, 2}) where T = nothing

"""
Scale each feature in dataset by standard deviation.

```math
μ_j = mean(x_j)
σ_j = std(x_j)
x_j^+ =  (x_j - μ_j) / σ_j

```

"""
struct NormalScaler <: AbstractScaler end

function scale!(::NormalScaler, X::Array{T, 2}) where T
    n, d = size(X)
    μ = mean(X, dims=1)
    σ = std(X, dims=1)

    @inbounds for i in 1:d
        X[:, i] .= (X[:, i] .- μ[i]) ./ σ[i]
    end
end

function scale!(::NormalScaler, X::AbstractSparseMatrix{T, Int}) where T
    n, d = size(X)
    σ = std(X, dims=1)
    @inbounds for i in 1:d
        X[:, i] ./= σ[i]
    end
end

"""
Format labels `y = [y_1, ..., y_n]` to ensure that `y_i ∈ { -1, 1}`
for all index i.

# Examples

```julia
    y = [-2.0, 1.0, -1.0, 2.0]
    format_label!(y)

```
"""
function format_label!(y::AbstractVector)
    set_reference = [-1.0, 1.0]
    elts = unique(y)
    if union(elts, set_reference) != set_reference
        intersets = intersect(elts, set_reference)
        replacements = setdiff(set_reference, intersets)
        # For all elements not in reference set
        count = 1
        for elt in setdiff(elts, set_reference)
            y[y .== elt] .= replacements[count]
            count += 1
        end
    end
    return nothing
end

# Return upper index (ordered by rows) of a pxp triangular matrix
function triul(p)
    index = Int[]
    for i in 1:p
        for j in i:p
            push!(index, i + (j - 1) * p)
        end
    end
    return index
end


text:
# Create example data
data = 2 .* randn(100) .- 1
# Draw the histogram with 19 bins
histogram(data, nbins=19)
# Horizontal histogram with log scale
histogram(data, horizontal=true, xlog=true)


text:
module EventsWrapper

using JavaCall

const Action = @jimport org.hipparchus.ode.events.Action
const BracketedRealFieldUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver
const BracketedUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedUnivariateSolver
const Class = @jimport java.lang.Class
const Enum = @jimport java.lang.Enum
const EventFilter = @jimport org.hipparchus.ode.events.EventFilter
const EventState = @jimport org.hipparchus.ode.events.EventState
const EventState_EventOccurrence = @jimport org.hipparchus.ode.events.EventState$EventOccurrence
const FieldEventState = @jimport org.hipparchus.ode.events.FieldEventState
const FieldEventState_EventOccurrence = @jimport org.hipparchus.ode.events.FieldEventState$EventOccurrence
const FieldODEEventHandler = @jimport org.hipparchus.ode.events.FieldODEEventHandler
const FieldODEState = @jimport org.hipparchus.ode.FieldODEState
const FieldODEStateAndDerivative = @jimport org.hipparchus.ode.FieldODEStateAndDerivative
const FieldODEStateInterpolator = @jimport org.hipparchus.ode.sampling.FieldODEStateInterpolator
const FilterType = @jimport org.hipparchus.ode.events.FilterType
const JString = @jimport java.lang.String
const ODEEventHandler = @jimport org.hipparchus.ode.events.ODEEventHandler
const ODEState = @jimport org.hipparchus.ode.ODEState
const ODEStateAndDerivative = @jimport org.hipparchus.ode.ODEStateAndDerivative
const ODEStateInterpolator = @jimport org.hipparchus.ode.sampling.ODEStateInterpolator
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const Transformer = @jimport org.hipparchus.ode.events.Transformer

export Action
export EventFilter
export EventState
export FieldEventState
export FieldODEEventHandler
export FilterType
export ODEEventHandler
export Transformer
export compare_to
export describe_constable
export do_event
export equals
export evaluate_step
export event_occurred
export g
export get_class
export get_convergence
export get_declaring_class
export get_event_handler
export get_event_time
export get_max_check_interval
export get_max_iteration_count
export hash_code
export init
export name
export notify
export notify_all
export ordinal
export reinitialize_begin
export reset_state
export to_string
export try_advance
export value_of
export values
export wait

include("action.jl")
include("event_filter.jl")
include("event_state.jl")
include("field_event_state.jl")
include("field_ode_event_handler.jl")
include("filter_type.jl")
include("ode_event_handler.jl")
include("transformer.jl")

end