text:
using TrajectoryOptimization
using HDF5
include("N_plots.jl")

# Set up the problem
model, obj0 = Dynamics.cartpole_analytical
n,m = model.n, model.m

x0 = zeros(n)
xf = [0.5;pi;0;0]
tf = 2.0
obj = LQRObjective(obj0.cost.Q, obj0.cost.R, obj0.cost.Qf, tf, x0, xf)
u_bnd = 20
x_bnd = [0.6,Inf,Inf,Inf]

obj_c = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd)
obj_min = update_objective(obj_c,tf=:min)
dt = 0.1

# Params
N = 51
method = :hermite_simpson

# Initial Trajectory
U0 = ones(1,N)
X0 = line_trajectory(obj.x0,obj.xf,N)

X0_rollout = copy(X0)
solver = Solver(model,obj_c,N=N)
rollout!(X0_rollout,U0,solver)

# Dircol functions
# eval_f, eval_g, eval_grad_f, eval_jac_g = gen_usrfun_ipopt(solver,method)

function comparison_plot(results,sol_ipopt; kwargs...)
    n = size(results.X,1)
    if n == 2 # pendulum
        state_labels = ["pos" "vel"]
    else
        state_labels = ["pos" "angle"]
    end
    time_dircol = range(0,stop=obj.tf,length=size(sol_ipopt.X,2))
    time_ilqr = range(0,stop=obj.tf,length=size(results.X,1))

    colors = [:blue :red]

    p_u = plot(time_ilqr[1:end-1],to_array(results.U)',width=2,label="iLQR",color=colors)
    plot!(time_dircol,sol_ipopt.U',width=2, label="IPOPT",color=colors,style=:dashdot,ylabel="control")

    p_x = plot(time_ilqr,to_array(results.X)[1:2,:]',width=2,label=state_labels,color=colors)
    plot!(time_dircol,sol_ipopt.X[1:2,:]',width=2, label="",color=colors,style=:dashdot,ylabel="state")

    p = plot(p_x,p_u,layout=(2,1),xlabel="time (s)"; kwargs...)

    return p
end

function convergence_plot(stat_i,stat_d)
    plot(stat_i["cost"],width=2,label="iLQR")
    plot!(stat_d["cost"], ylim=[0,10], ylabel="Cost",xlabel="iterations",width=2,label="dircol")
end

colors_X = [:red :blue :orange :green]

#####################################
#          UNCONSTRAINED            #
#####################################

# Solver Options
opts = SolverOptions()
opts.cost_tolerance = 1e-6
opts.outer_loop_update_type = :default
opts.constraint_decrease_ratio = 0.75

# iLQR
solver = Solver(model,obj,N=N,opts=opts,integration=:rk3)
solver.opts.verbose = false
res_i, stat_i = solve(solver,U0)
plot(to_array(res_i.X)')
plot(to_array(res_i.U)')
max_violation(res_i)
cost(solver, res_i)
stat_i["runtime"]
stat_i["iterations"]
var_i = DircolVars(to_array(res_i.X),to_array(res_i.U))

# DIRCOL
res_d, stat_d = solve_dircol(solver, X0_rollout, U0; method=method)
plot(res_d.X')
plot(res_d.U')
stat_d["c_max"][end]
stat_d["cost"][end]
stat_d["runtime"]
stat_d["iterations"]

comparison_plot(res_i,res_d)
convergence_plot(stat_i,stat_d)

time_per_iter = stat_i["runtime"]/stat_i["iterations"]
time_per_iter = stat_d["runtime"]/stat_d["iterations"]

#####################################
#           CONSTRAINED             #
#####################################

# Solver Options
opts = SolverOptions()
opts.cost_tolerance = 1e-8
opts.cost_tolerance_intermediate = 1e-2
opts.constraint_tolerance = 1e-8
opts.outer_loop_update_type = :feedback
opts.constraint_decrease_ratio = .25
opts.square_root = true
opts.penalty_initial = 1
opts.penalty_scaling = 50

# iLQR
obj_c2 = update_objective(obj_c,tf=2)
solver = Solver(model,obj_c2,N=N,opts=opts,integration=:rk3)
res_i, stat_i = solve(solver,U0)
stat_i["iterations"]
comparison_plot(res_i,res_d)
plot(to_array(res_i.X)')
plot(to_array(res_i.U)')
stat_i["c_max"][end]
stat_i["runtime"]
constraint_plot(solver,U0)

# DIRCOL
res_d, stat_d = solve_dircol(solver, X0_rollout, U0; method=:hermite_simpson)
plot(res_d.X')
plot(res_d.U')
stat_d["c_max"][end]
stat_d["cost"][end]
stat_d["runtime"]
stat_d["iterations"]

convergence_plot(stat_i,stat_d)

cost(solver,res_i)
cost(solver,res_d)

time_per_iter = stat_i["runtime"]/stat_i["iterations"]
time_per_iter = stat_d["runtime"]/stat_d["iterations"]

#####################################
#        INFEASIBLE START           #
#####################################

# Solver Options
opts = SolverOptions()
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 5e-1
opts.constraint_tolerance = 1e-4
opts.outer_loop_update_type = :feedback
opts.constraint_decrease_ratio
opts.resolve_feasible = false
opts.square_root = true  # this makes a big difference!
opts.R_infeasible = 1
opts.penalty_initial = 1
opts.penalty_scaling = 10
opts.use_nesterov = false
opts.al_type = :default

# iLQR
solver = Solver(model,obj_c,N=N,opts=opts,integration=:rk3)
res_i, stat_i = solve(solver,X0,U0)
plot(to_array(res_i.X)')
plot(to_array(res_i.U)')
stat_i["c_max"][end]
stat_i["cost"][end]
stat_i["runtime"]
stat_i["iterations"]
var_i = DircolVars(res_i)

# DIRCOL
res_d, stat_d = solve_dircol(solver, X0, U0; method=:hermite_simpson)
plot(res_d.X')
plot(res_d.U')
stat_d["c_max"][end]
stat_d["cost"][end]
stat_d["runtime"]
stat_d["iterations"]
stat_d["cost"]

comparison_plot(res_i,res_d)
convergence_plot(stat_i,stat_d)

cost(solver,res_i)
cost(solver,res_d)

time_per_iter = stat_i["runtime"]/stat_i["iterations"]
time_per_iter = stat_d["runtime"]/stat_d["iterations"]

#####################################
#          MINIMUM TIME             #
#####################################

X0_warm = to_array(res_i.X)
U0_warm = to_array(res_i.U)

# Solver Options
opts = SolverOptions()
opts.max_dt = 0.7
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 0.01
opts.R_minimum_time = 100
opts.penalty_initial_minimum_time_equality = 100.
opts.penalty_scaling_minimum_time_equality = 10

opts.outer_loop_update_type = :default
opts.constraint_decrease_ratio
opts.resolve_feasible = false
opts.square_root = true  # this makes a big difference!
opts.penalty_initial = 1
opts.penalty_scaling = 10
opts.use_nesterov = false
opts.al_type = :default
opts.live_plotting = true
opts.resolve_feasible = false
opts.R_infeasible = 1
opts.minimum_time_tf_estimate = 1.5

# iLQR
solver_mintime = Solver(model,obj_min,N=N,opts=opts,integration=:rk3)
solver_mintime.opts.verbose = true
res_i, stat_i = solve(solver_mintime,X0_warm,U0_warm)
plot(to_array(res_i.X)')
plot(to_array(res_i.U)')
max_violation(res_i)
_cost(solver, res_i)
stat_i["runtime"]
stat_i["iterations"]
total_time(solver, res_i)
plot(stat_i["c_max"],yscale=:log10)

# DIRCOL
solver.opts.verbose = false
res_d, stat_d = solve_dircol(solver, X0_rollout, U0; method=:hermite_simpson)
stat_d["info"]
plot(res_d.X')
plot(res_d.U')
stat_d["c_max"][end]
stat_d["cost"][end]
stat_d["runtime"]
stat_d["iterations"]
stat_d["cost"]
total_time(solver, res_d)

comparison_plot(res_i,res_d)
convergence_plot(stat_i,stat_d)

eval_f(var_i.Z)
eval_f(res_d.Z)

stat_i["cost"][30]
stat_d["cost"][30]

time_per_iter = stat_i["runtime"]/stat_i["iterations"]
time_per_iter = stat_d["runtime"]/stat_d["iterations"]

dt_i = [res_i.U[k][end]^2 for k = 1:N]
dt_d = res_d.U[end,:]

plot(dt_i)
plot!(dt_d)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#                                                                              #
#                       KNOT POINT COMPARISONS                                 #
#                                                                              #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

Ns = [21,41,51,81,101,201,401,501,801,1001]
obj.tf ./ (Ns.-1)
dt_truth = 1e-3

#####################################
#          UNCONSTRAINED            #
#####################################
group = "cartpole/unconstrained"
run_step_size_comparison(model, obj, U0, group, Ns, integrations=[:midpoint,:rk3,:rk3_foh,:rk4],dt_truth=1e-3,opts=opts)
plot_stat("iterations",group)

#####################################
#            CONSTRAINED            #
#####################################
Ns = [21,41,51,81,101,201]
group = "cartpole/constrained"
plot_stat("error",group)

#####################################
#            INFEASIBLE             #
#####################################
solver_truth, res_truth,  = run_dircol_truth(model, obj_c, 1e-3, X0, U0, "cartpole/infeasible")
time_truth = get_time(solver_truth)
plot(res_truth.U')

err_mid, eterm_mid, stats_mid = run_Ns(model, obj_c, Ns, :midpoint, infeasible=true)
err_rk3, eterm_rk3, stats_rk3 = run_Ns(model, obj_c, Ns, :rk3, infeasible=true)
err_foh, eterm_foh, stats_foh = run_Ns(model, obj_c, Ns, :rk3_foh, infeasible=true)
err_rk4, eterm_rk4, stats_rk4 = run_Ns(model, obj_c, Ns, :rk4, infeasible=true)

save_data("cartpole/infeasible")


text:
# Read all entries in PCREO_DATABASE_DIRECTORY, compute canonical
# representatives of their nearest-neighbor graphs, and store the resulting
# graphs in PCREO_GRAPH_DIRECTORY in graph6 format. Does not recompute graphs
# that already exist in PCREO_GRAPH_DIRECTORY.

using Base.Threads: @threads, threadid, nthreads
using ProgressMeter
using UUIDs: UUID

push!(LOAD_PATH, @__DIR__)
using PCREO

function get_or_compute_canonical_graph(entrypath::AbstractString,
                                        graphpath::AbstractString)
    if isfile(graphpath)
        return readchomp(graphpath)
    else
        result = canonical_graph6(entrypath)
        write(graphpath, result, '\n')
        return result
    end
end

function main()

    println("Using $(nthreads()) threads.")
    flush(stdout)

    @assert isdir(ENV["PCREO_DATABASE_DIRECTORY"])
    @assert isdir(ENV["PCREO_GRAPH_DIRECTORY"])

    for num_dir in lsdir(ENV["PCREO_DATABASE_DIRECTORY"]; join=true)

        println("Reading PCREO database directory $num_dir...")
        flush(stdout)

        entry_paths = [
            entry_path
            for entry_dir in lsdir(num_dir; join=true)
            for entry_path in lsdir(entry_dir; join=true)
        ]

        n = length(entry_paths)
        canonical_graphs = Vector{String}(undef, n)
        uuids = Vector{UUID}(undef, n)

        println("Loading canonical nearest-neighbor graphs...")
        flush(stdout)

        p = Progress(n; dt=0.05, desc=basename(num_dir),
                        output=stdout, showspeed=true)

        @threads for i = 1 : n
            @inbounds entry = entry_paths[i]
            path = splitpath(entry)
            @inbounds uuids[i] = UUID(path[end-1])
            @inbounds name = path[end]
            @assert endswith(name, ".csv")
            @inbounds canonical_graphs[i] = get_or_compute_canonical_graph(
                entry,
                joinpath(ENV["PCREO_GRAPH_DIRECTORY"], name[1:end-3] * "g6")
            )
            next!(p)
        end

        flush(stdout)
        println("Verifying uniqueness...")
        flush(stdout)

        graph_dict = Dict{UUID,Set{String}}()
        for (graph, uuid) in zip(canonical_graphs, uuids)
            if haskey(graph_dict, uuid)
                push!(graph_dict[uuid], graph)
            else
                new_set = Set{String}()
                push!(new_set, graph)
                graph_dict[uuid] = new_set
            end
        end

        for (uuid, graph_set) in graph_dict
            @assert length(graph_set) == 1
        end

        @assert allunique(
            graph
            for (uuid, graph_set) in graph_dict
            for graph in graph_set
        )
    end
end

main()

ext: jl


text:
struct Rule
    predecessor::String

    successor_1::String
    successor_2::String

    reverse_input::Bool
    reverse_output::Bool

    inhibitory_1_2::Bool
    excitatory_1_2::Bool

    inhibitory_2_1::Bool
    excitatory_2_1::Bool

    start_index::UInt32
    end_index::UInt32
end

Base.length(rule::Rule) = length(rule.predecessor)
Base.isequal(rule1::Rule, rule2::Rule) = hash(rule1) == hash(rule2)
function Base.hash(rule1::Rule)
    hash((rule1.predecessor, rule1.successor_1, rule1.successor_2, rule1.reverse_input, rule1.reverse_output, rule1.inhibitory_1_2, rule1.excitatory_1_2, rule1.inhibitory_2_1, rule1.excitatory_2_1))
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module Environment
using  ...PatternStructure.Trees
using  ...PatternStructure.Special
using  ...PatternStructure.Checks
using  ...PatternStructure.SlurpTypes
using  ...Analyzer.SlurpOptimizations
using  ...Helper
export Variables, SlurpRanges, match_variable!, unmatch_variable!, extract!, add_slurp_range!

const Variables   = Dict{Symbol, Any}
const SlurpRanges = Dict{PatternNode, Vector{Range}}

function match_variable!(v::Variables, name::Symbol, value)
  haskey(v, name) ?
    (v[name] == value) :
    (v[name] =  value; true)
end

function extract!(v::Variables, tree::PatternNode, args, ai, mi, ranges)
  for b in bindings(tree)
    v[b] = Any[]
  end
  extract_slurp!(tree.head, v, tree, args, ai:mi, ranges)

end

function extract_values!(v::Variables, tree::PatternNode, value, ranges)
  args  = tree.step(value)
  index = 1

  for child in tree.children
    if is_slurp(child)
      range = shift!(ranges[child])
      index = last(range)+1

      for b in bindings(child)
        add_binding_iteration!(child.depth, b,v)
      end

      extract_slurp!(child.head, v, child, args, range, ranges)

    else
      extract_values!(v, child, args[index], ranges)
      index += 1

    end
  end
end

function extract_values!(v::Variables, tree::PatternGate, value, ranges)
  if isa(tree.check, Binding)
     binding = get_binding(v, tree.check.name, tree.depth)
     push!(binding, value)
  end

  extract_values!(v, tree.child, value, ranges)
end

function extract_values!(v::Variables, tree::PatternLeaf, value, ranges)
end

function extract_slurp!(h, v, tree, args, range, ranges)
  child = Looping(tree.children)
  for i in range
    extract_values!(v, current(child), args[i], ranges)
    next!(child)
  end
end

function extract_slurp!(h::SimpleSlurp, v, tree, args, range, ranges)
  binding = get_binding(v, tree.children[1].check.name, tree.depth)
  append!(binding, args[range])
end

function get_binding(v, name, depth)
  binding = v[name]

  while depth > 1
    depth  -= 1
    binding = binding[end]
  end

  return binding
end

function add_binding_iteration!(depth, b, v)
  binding = v[b]
  while depth > 2
    depth  -= 1
    binding = binding[end]
  end
  push!(binding, Any[])
end

function add_slurp_range!(ranges, slurp, range)
  if !haskey(ranges, slurp)
    ranges[slurp] = []
  end
  push!(ranges[slurp], range)
end

function unmatch_variable!(v::Variables, names)
  for name in names
  	if v.count[name] >  1
  	   v.count[name] -= 1
  	else
  	   delete!(v.count,  name)
  	   delete!(v.values, name)
  	end
  end
end

end


text:
module MiraiBots

using HTTP, JSON3, StructTypes, Dates

const Optional{T} = Union{T,Nothing}

# TODO: really make them different
const FriendId = Int
const GroupId = Int
const GroupOrFriendId = Union{FriendId,GroupId}
const MessageId = Int
const TimeStamp = Int
const DurationSeconds = Second
const EventId = Int
const DeviceKind = Int

include("MessageChains.jl")
include("EventsAndMessages.jl")
include("Commands.jl")

"""
    ProtocolAdapter

Adapter for different protocols that Mirai HTTP API provides.

- `loop`: start sending and receiving. The function blocks, so you would like to call it as `@async loop(...)`.
- `send`: post a message, wait and return the result. Because the adapter need to be initialized, do not send until the first message is received.
- `receive`: wait and retrieve a message (serialized when possible) or exception. The first message is always session info.
- `receive_or_throw`: like `receive`, but throws exception if an exception is received.
- `get_output_channel`: alternative to `receive`, you can manipulate the underlying channel directly.
- `close`: attempt to shutdown the adapter.
"""
abstract type ProtocolAdapter end

include("Broadcaster.jl")

get_output_channel(adp::ProtocolAdapter) = adp.output_channel
receive(adp::ProtocolAdapter) = take!(get_output_channel(adp))

function receive_or_throw(adp::ProtocolAdapter, show::Bool = true)
    msg = receive(adp)
    if msg isa ExceptionAndBacktrace
        show && @error "Error occurred receiving message" exception = (msg.exception, msg.backtrace)
        throw(msg.exception)
    else
        msg
    end
end

const sync_id_counter = Ref(0)
get_sync_id() = string(sync_id_counter[] += 1)

"A command object suitable for both adapters."
struct GeneralCommand{T}
    command::Symbol
    subCommand::Optional{Symbol}
    content::T
    method::Commands.CommandMethod
end

function make_command(cmd::Commands.AbstractCommand)
    GeneralCommand(Commands.command(cmd), Commands.subcommand(cmd), cmd, Commands.method(cmd))
end

struct ExceptionAndBacktrace
    exception::Exception
    backtrace::Union{Base.return_types(catch_backtrace, ())...}
end

const OutputChannel = Channel{Union{JSON3.Object,EventOrMessage,ExceptionAndBacktrace}}
const DEFAULT_BUFFER_SIZE = 16

function with_catch_backtrace(f, output_channel)
    try
        f()
    catch e
        put!(output_channel, ExceptionAndBacktrace(e, catch_backtrace()))
    end
end

struct AdapterConectionFailed <: Exception
    code::Int
    reason::String
end

include("WebSocketAdapter.jl")
include("HTTPAdapterBase.jl")
include("HTTPAdapter.jl")
include("HTTPCompatAdapter.jl")

function try_convert(T, x)
    try
        StructTypes.constructfrom(T, x)
    catch e
        @debug "Error deserializing data" exception = (e, catch_backtrace())
        x
    end
end

struct RESTfulRequestFailed <: Exception
    cmd::Commands.AbstractCommand
    code::Int
    msg::String
end

function send(adp::ProtocolAdapter, cmd::Commands.AbstractCommand, response_type_func = Commands.response_type;
    session_key_position::SessionKeyPosition = SESSION_KEY_IN_HEADERS, log::Bool = true)
    log && @info "Send $cmd"
    resp = send(adp, make_command(cmd); session_key_position)
    ret = try_convert(response_type_func(cmd), resp)
    if response_type_func(cmd) <: Union{Commands.RESTful,Commands.RESTfulErrorMessage} &&
            hasproperty(ret, :code) && ret.code != 0
        msg = hasproperty(ret, :msg) ? ret.msg : ""
        throw(RESTfulRequestFailed(cmd, ret.code, msg))
    else
        log && @info "Receive $ret"
        ret
    end
end

"""
    @send cmd(args...)

`@send cmd(args...)` is short for `send(bot, MiraiBots.Commands.cmd(args...)`.
"""
macro send(cmd)
    send_imp(cmd)
end

function send_imp(cmd::Expr)
    @assert cmd.head == :call && cmd.args[1] isa Symbol
    esc(Expr(:call, send, :bot,
        Expr(:call, Expr(:., MiraiBots.Commands, Meta.quot(cmd.args[1])), cmd.args[2:end]...)))
end

end


text:
using LinearAlgebra: norm_sqr

export Sphere

struct Sphere{T, M} <: Hittable
    center::Point3{T}
    radius::T
    material::M
end

function hit(s::Sphere, r::Ray, t_min, t_max)::Union{HitRecord, Nothing}
    oc = r.orig - s.center
    a = norm_sqr(r.dir)
    half_b = oc ⋅ r.dir
    c = norm_sqr(oc) - s.radius^2
    discriminant = half_b^2 - a*c
    if discriminant < 0
        return nothing
    end
    sqrtd = sqrt(discriminant)

    # Find the nearest root that lies in the acceptable range.
    root = (-half_b - sqrtd) / a
    if root < t_min || t_max < root
        root = (-half_b + sqrtd) / a
        if root < t_min || t_max < root
            return nothing
        end
    end

    t = root
    p = at(r, t)
    outward_normal = (p - s.center) / s.radius

    return HitRecord(r, p, outward_normal, s.material, t)
end


text:
module AirSeaFluxes

export simpleflux, bulkformulae
include("simpleflux.jl")
include("bulkformulae.jl")

end # module


text:
# ---
# title: 1351. Count Negative Numbers in a Sorted Matrix
# id: problem1351
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/count-negative-numbers-in-a-sorted-matrix/description/>
# hidden: true
# ---
#
# Given a `m * n` matrix `grid` which is sorted in non-increasing order both
# row-wise and column-wise.
#
# Return the number of **negative** numbers in `grid`.
#
#
#
# **Example 1:**
#
#
#
#     Input: grid = [[4,3,2,-1],[3,2,1,-1],[1,1,-1,-2],[-1,-1,-2,-3]]
#     Output: 8
#     Explanation: There are 8 negatives number in the matrix.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[3,2],[1,0]]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[1,-1],[-1,-1]]
#     Output: 3
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[-1]]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m, n <= 100`
#   * `-100 <= grid[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using GenerativeModels
using POMDPs
using POMDPToolbox
using ProfileView
using Multilane
using MCTS
using RobustMCTS

nb_lanes = 4
pp = PhysicalParam(nb_lanes,lane_length=100.) #2.=>col_length=8\n",
_discount = 1.
nb_cars = 10
rmodel = NoCrashRewardModel()
dmodel = NoCrashIDMMOBILModel(nb_cars, pp)
mdp = NoCrashMDP(dmodel, rmodel, _discount);
rng = MersenneTwister(5)

# dpws = DPWSolver(depth=20,
#                  n_iterations=500,
#                  exploration_constant=100.0,
#                  k_state=4.0,
#                  alpha_state=1/8,
#                  rollout_solver=SimpleSolver())
#
# policy = solve(dpws, mdp)

# policy = RandomPolicy(mdp, rng=rng)

rsolver = RobustMCTSSolver(
    depth=20,
    c = 100.0,
    n_iterations=500,
    k_state=4.0,
    alpha_state=1/8,
    k_nature=4.0,
    alpha_nature=1/8,
    rollout_solver=SimpleSolver(),
    rollout_nature=StochasticBehaviorNoCrashMDP(mdp))

solver = RobustMLSolver(rsolver)
policy = solve(solver, mdp)

sim = RolloutSimulator(max_steps=10)

@time simulate(sim, mdp, policy, initial_state(mdp, rng))
@time for i in 1
    simulate(sim, mdp, policy, initial_state(mdp, rng))
end

Profile.clear()
@profile for i in 1
    simulate(sim, mdp, policy, initial_state(mdp, rng))
end

ProfileView.view()


text:
@testset "Interface" begin
    @test PCA(5) == PCA(ncomps=5)
    @test isnothing(PCA().ncomps)
    @test PCA(5).pratio == 1
end

@testset "Decomposition" begin
    data = 4 .* randn(rng, 30, 101, 101) .+ 10
    angles = sort!(90randn(rng, 30)) |> normalize_par_angles

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        A, w = @inferred decompose(PCA(ncomps=N), data, angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (30, N)
    end
    @test_throws ErrorException decompose(PCA(40), data, angles)
    A, w = decompose(PCA(30; pratio=0.5), data, angles)
    @test size(A, 1) < 30
    @test size(w, 2) < 30

    # default is to use whole cube
    S = reconstruct(PCA(), data, angles)
    @test size(S) == (30, 101, 101)
end

@testset "ADI Trivial" begin
    cube = ones(10, 100, 100)
    angles = zeros(10)

    reduced_5 = PCA(5)(cube, angles)
    reduced_10 = PCA(10)(cube, angles)

    @test size(reduced_5) == size(reduced_10) == (100, 100)
    @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_5)
    @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_10)
end

@testset "RDI Trivial" begin
    data = randn(rng, 30, 101, 101)
    angles = sort!(90rand(rng, 30)) |> normalize_par_angles

    S = data .- reconstruct(PCA(1), data, angles, zeros(10, 101, 101))
    @test S ≈ data rtol=2e-2
end


text:
const colours = Dict(0=>:black,1=>:red,2=>:green,3=>:blue,4=>:cyan,5=>:magenta,6=>:yellow,7=>:orange)

## HierarchicalDomain

@recipe function f(H::HierarchicalDomain)
    H1,H2 = partition(H)
    @series begin H1 end
    @series begin H2 end
end


text:
using Cropbox
using CropRootBox
using Test

# using DataFrames
# using Statistics
# import Gadfly

root_maize = (
    :RootArchitecture => :maxB => 5,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 1 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.1 ± 0.01,
        la = 18.0 ± 1.8,
        ln = 0.6 ± 0.06,
        lmax = 89.7 ± 7.4,
        r = 6.0 ± 0.6,
        Δx = 0.5,
        σ = 10,
        θ = 80 ± 8,
        N = 1.5,
        a = 0.04 ± 0.004,
        color = CropRootBox.RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.2 ± 0.04,
        la = 0.4 ± 0.04,
        ln = 0.4 ± 0.03,
        lmax = 0.6 ± 1.6,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 15,
        N = 1,
        a = 0.03 ± 0.003,
        color = CropRootBox.RGBA(0, 1, 0, 1),
    ),
    :SecondOrderLateralRoot => (;
        lb = 0,
        la = 0.4 ± 0.02,
        ln = 0,
        lmax = 0.4,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 10,
        N = 2,
        a = 0.02 ± 0.002,
        color = CropRootBox.RGBA(0, 0, 1, 1),
    )
)
root_switchgrass_N = (
    :RootArchitecture => :maxB => 24,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 1 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.41 ± 0.26,
        la = 0.63 ± 0.50,
        ln = 0.27 ± 0.07,
        lmax = 33.92 ± 22.60,
        r = 1 ± 0.1,
        Δx = 0.5,
        σ = 9,
        θ = 60 ± 6,
        N = 1.5,
        a = (0.62 ± 0.41)u"mm",
        color = CropRootBox.RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.63 ± 0.45,
        la = 1.12 ± 1.42,
        ln = 0.23 ± 0.11,
        lmax = 7.03 ± 6.84,
        r = 0.21 ± 0.02,
        Δx = 0.1,
        σ = 18,
        θ = 60 ± 6,
        N = 1,
        a = (0.22 ± 0.07)u"mm",
        color = CropRootBox.RGBA(0, 1, 0, 1),
    ),
    :SecondOrderLateralRoot => (;
        lb = 0.45 ± 0.64,
        la = 0.71 ± 0.64,
        ln = 0.14 ± 0.10,
        lmax = 2.77 ± 2.68,
        r = 0.08 ± 0.01,
        Δx = 1,
        σ = 20,
        θ = 60 ± 6,
        N = 2,
        a = (0.19 ± 0.10)u"mm",
        color = CropRootBox.RGBA(0, 0, 1, 1),
    )
)
root_switchgrass_W = (
    :RootArchitecture => :maxB => 22,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 1 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.40 ± 0.45,
        la = 3.37 ± 2.89,
        ln = 0.31 ± 0.06,
        lmax = 42.33 ± 30.54,
        r = 1 ± 0.1,
        Δx = 0.5,
        σ = 9,
        θ = 60 ± 6,
        N = 1.5,
        a = (0.71 ± 0.65)u"mm",
        color = CropRootBox.RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.59 ± 0.48,
        la = 0.74 ± 0.87,
        ln = 0.09 ± 0.09,
        lmax = 1.79 ± 1.12,
        r = 0.04 ± 0.01,
        Δx = 0.1,
        σ = 18,
        θ = 60 ± 6,
        N = 1,
        a = (0.19 ± 0.08)u"mm",
        color = CropRootBox.RGBA(0, 1, 0, 1),
    ),
    :SecondOrderLateralRoot => (;
        lb = 0,
        la = 0.07 ± 0.04,
        ln = 0,
        lmax = 0.07 ± 0.04,
        r = 0.002 ± 0.001,
        Δx = 1,
        σ = 20,
        θ = 60 ± 6,
        N = 2,
        a = (0.19 ± 0.08)u"mm", # 0.68 ± 0.77
        color = CropRootBox.RGBA(0, 0, 1, 1),
    )
)
root_switchgrass_N2 = (
    :RootArchitecture => :maxB => 24,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 0 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.67 ± 0.25,
        la = 0.58 ± 0.96,
        ln = 0.25 ± 0.06,
        lmax = 11.85 ± 12.63,
        r = 1 ± 0.1,
        Δx = 0.5,
        σ = 9,
        θ = 60 ± 6,
        N = 1.5,
        a = (0.78 ± 0.27)u"mm",
        color = CropRootBox.RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.20 ± 0.15,
        la = 0.94 ± 1.12,
        ln = 0.12 ± 0.06,
        lmax = 8.04 ± 7.57,
        r = 0.68 ± 0.07,
        Δx = 0.1,
        σ = 18,
        θ = 60 ± 6,
        N = 1,
        a = (0.35 ± 0.38)u"mm",
        color = CropRootBox.RGBA(0, 1, 0, 1),
    ),
)
root_switchgrass_W2 = (
    :RootArchitecture => :maxB => 22,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 0 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.08 ± 0.07,
        la = 4.32 ± 1.86,
        ln = 0.21 ± 0.05,
        lmax = 22.80 ± 11.94,
        r = 1 ± 0.1,
        Δx = 0.5,
        σ = 9,
        θ = 60 ± 6,
        N = 1.5,
        a = (0.40 ± 0.15)u"mm",
        color = CropRootBox.RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.41 ± 0.35,
        la = 1.39 ± 1.00,
        ln = 0.10 ± 0.13,
        lmax = 2.37 ± 1.41,
        r = 0.10 ± 0.01,
        Δx = 0.1,
        σ = 18,
        θ = 60 ± 6,
        N = 1,
        a = (0.21 ± 0.06)u"mm",
        color = CropRootBox.RGBA(0, 1, 0, 1),
    ),
)

root_switchgrass_P = (root_switchgrass_W,
    :Clock => (;
        step = 1u"hr",
    ),
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 1 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        Δx = 0.5,
    ),
    :FirstOrderLateralRoot => (;
        Δx = 0.1,
    ),
    :SecondOrderLateralRoot => (;
        Δx = 1,
    )
)
root_switchgrass_KH2PO4 = (
    root_switchgrass_P,
    :PrimaryRoot => (;
        r = 1.28 ± 0.29,
    ),
    :FirstOrderLateralRoot => (;
        lmax = 1.63 ± 0.61,
        #ln = 0.163 ± 0.061,
        r = 0.31 ± 0.11,
        θ = 84.27 ± 3.03,
    ),
    :SecondOrderLateralRoot => (;
        lmax = 2.30 ± 0.65,
        r = 0.48 ± 0.15,
        θ = 41.87 ± 6.64,
    )
)
root_switchgrass_AlPO4 = (
    root_switchgrass_P,
    :PrimaryRoot => (;
        r = 1.10 ± 0.32,
    ),
    :FirstOrderLateralRoot => (;
        lmax = 1.56 ± 0.55,
        #ln = 0.156 ± 0.055,
        r = 0.33 ± 0.11,
        θ = 83.53 ± 3.42,
    ),
    :SecondOrderLateralRoot => (;
        lmax = 2.29 ± 0.60,
        r = 0.48 ± 0.13,
        θ = 34.87 ± 5.17,
    )
)
root_switchgrass_C6H17NaO24P6 = (
    root_switchgrass_P,
    :PrimaryRoot => (;
        r = 1.31 ± 0.35,
    ),
    :FirstOrderLateralRoot => (;
        lmax = 2.80 ± 0.61,
        #ln = 0.28 ± 0.061,
        r = 0.57 ± 0.09,
        θ = 78.07 ± 5.30,
    ),
    :SecondOrderLateralRoot => (;
        lmax = 3.05 ± 0.50,
        r = 0.65 ± 0.17,
        θ = 43.73 ± 6.71,
    )
)

container_pot = :Pot => (;
    r1 = 10,
    r2 = 6,
    height = 30,
)
container_rhizobox = :Rhizobox => (;
    l = 16u"inch",
    w = 10.5u"inch",
    h = 42u"inch",
)
soilcore = :SoilCore => (;
    d = 5,
    l = 20,
    x0 = 3,
    y0 = 3,
)

@testset "root" begin
    b = instance(CropRootBox.Pot, config=container_pot)
    s = instance(CropRootBox.RootArchitecture; config=root_maize, options=(; box=b), seed=0)
    r = simulate!(s, stop=100u"d")
    @test r.time[end] == 100u"d"
    # using GLMakie
    # scn = CropRootBox.render(s)
    # GLMakie.save("root_maize.png", scn)
    CropRootBox.writevtk(tempname(), s)
    # CropRootBox.writepvd(tempname(), CropRootBox.RootArchitecture, config=root_maize, stop=50)
    CropRootBox.writestl(tempname(), s)
end

# @testset "switchgrass" begin
#     C = Dict(
#         :KH2PO4 => root_switchgrass_KH2PO4,
#         :AlPO4 => root_switchgrass_AlPO4,
#         :C6H17NaO24P6 => root_switchgrass_C6H17NaO24P6,
#     )
#     b = instance(CropRootBox.Rhizobox, config=container_rhizobox)
#     P = [1, 4, 8, 12, 16, 20, 24]u"wk"
#     R = []
#     for i in 1:10, c in (:KH2PO4, :AlPO4, :C6H17NaO24P6)
#         n = "$c-$i"
#         r = simulate(CropRootBox.RootArchitecture; config=C[c], options=(; box=b), seed=i, stop=P[end]) do D, s
#             t = s.context.clock.time' |> u"wk"
#             if t in P
#                 p = deunitfy(t, u"wk") |> Int
#                 CropRootBox.writevtk("$n-w$p", s)
#                 G = gather!(s, CropRootBox.BaseRoot; callback=CropRootBox.gatherbaseroot!)
#                 D[1][:time] = t
#                 D[1][:treatment] = c
#                 D[1][:repetition] = i
#                 D[1][:length] = !isempty(G) ? sum([s.length' for s in G]) : 0.0u"cm"
#                 D[1][:volume] = !isempty(G) ? sum([s.length' * s.radius'^2 for s in G]) : 0.0u"mm^3"
#                 D[1][:count] = length(G)
#             else
#                 empty!(D)
#             end
#         end
#         push!(R, r)
#     end
#     df = vcat(R...)
#     combine(groupby(df, [:treatment, :time]), :length => mean, :length => std) |> display
#     Gadfly.plot(df, x=:time, y=:length, color=:treatment,
#         Gadfly.Geom.boxplot,
#         Gadfly.Scale.x_discrete,
#         Gadfly.Guide.xlabel("Time"),
#         Gadfly.Guide.ylabel("Total Root Length Per Plant"),
#         Gadfly.Guide.colorkey(title="", pos=[0.05*Gadfly.w, -0.4*Gadfly.h]),
#         Gadfly.Theme(boxplot_spacing=7*Gadfly.mm)
#     ) |> Gadfly.SVG("switchgrass_P.svg")
# end

# using DataStructures: OrderedDict
# @testset "switchgrass rhizobox" begin
#     C = OrderedDict(
#         :W => root_switchgrass_W,
#         :N => root_switchgrass_N,
#     )
#     R = []
#     for (k, c) in C
#         n = "switchgrass_$k"
#         b = instance(CropRootBox.Rhizobox, config=container_rhizobox)
#         s = instance(CropRootBox.RootArchitecture; config=c, options=(; box=b), seed=0)
#         r = simulate!(s, stop=1000) do D, s
#             G = gather!(s, CropRootBox.BaseRoot; callback=CropRootBox.gatherbaseroot!)
#             D[1][:length] = !isempty(G) ? sum([s.length' for s in G]) : 0.0u"cm"
#             D[1][:volume] = !isempty(G) ? sum([s.length' * s.radius'^2 for s in G]) : 0.0u"mm^3"
#             D[1][:count] = length(G)
#         end
#         push!(R, r)
#         plot(r, :time, :length, kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("$n-length.pdf")
#         plot(r, :time, :volume, yunit=u"mm^3", kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("$n-volume.pdf")
#         plot(r, :time, :count, kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("$n-count.pdf")
#         CropRootBox.writevtk(n, s)
#     end
#     save(x, y, f; kw...) = begin
#         K = collect(keys(C))
#         p = plot(R[1], x, y; name=string(K[1]), title=string(y), kind=:line, backend=:Gadfly, kw...)
#         for i in 2:length(R)
#             p = plot!(p, R[i], x, y; name=string(K[i]), kind=:line, backend=:Gadfly, kw...)
#         end
#         p[] |> Cropbox.Gadfly.PDF(f)
#     end
#     save(:time, :length, "switchgrass-length.pdf")
#     save(:time, :volume, "switchgrass-volume.pdf"; yunit=u"mm^3")
#     save(:time, :count, "switchgrass-count.pdf")
# end

# @testset "switchgrass layer" begin
#     C = Dict(
#         :W => root_switchgrass_W,
#         :N => root_switchgrass_N,
#     )
#     for (k, c) in C
#         n = "switchgrass_$k"
#         b = instance(CropRootBox.Rhizobox, config=container_rhizobox)
#         L = [instance(CropRootBox.SoilLayer, config=:SoilLayer => (; d, t=1)) for d in 0:1:10]
#         s = instance(CropRootBox.RootArchitecture; config=c, options=(; box=b), seed=0)
#         r = simulate!(s, stop=300) do D, s
#             G = gather!(s, CropRootBox.BaseRoot; callback=CropRootBox.gatherbaseroot!)
#             for (i, l) in enumerate(L)
#                 V = [s.length' for s in G if s.ii(l)]
#                 D[1][Symbol("L$(i-1)")] = !isempty(V) ? sum(V) : 0.0u"cm"
#             end
#         end
#         plot(r, :time, [Symbol("L$(i-1)") for i in 1:length(L)], kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("L$n.pdf")
#         CropRootBox.writevtk(n, s)
#     end
# end

# @testset "maize layer" begin
#     n = "maize"
#     b = instance(CropRootBox.Pot, config=container_pot)
#     t = 5
#     L = [instance(CropRootBox.SoilLayer, config=:SoilLayer => (; d, t)) for d in 0:t:30-t]
#     s = instance(CropRootBox.RootArchitecture; config=root_maize, options=(; box=b), seed=0)
#     r = simulate!(s, stop=500) do D, s
#         G = gather!(s, CropRootBox.BaseRoot; callback=CropRootBox.gatherbaseroot!)
#         for (i, l) in enumerate(L)
#             ll = [s.length' for s in G if s.ii(l)]
#             D[1][Symbol("L$(i-1)")] = !isempty(ll) ? sum(ll) : 0.0u"cm"
#             vl = [s.length' * s.radius'^2 for s in G if s.ii(l)]
#             D[1][Symbol("V$(i-1)")] = !isempty(vl) ? sum(vl) : 0.0u"cm^3"
#             D[1][Symbol("C$(i-1)")] = length(vl)
#         end
#     end
#     plot(r, :time, [Symbol("L$(i-1)") for i in 1:length(L)], legend="soil depth", names=["$((i-1)*t) - $(i*t) cm" for i in 1:length(L)], title="total length", kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("L$n.pdf")
#     plot(r, :time, [Symbol("V$(i-1)") for i in 1:length(L)], legend="soil depth", names=["$((i-1)*t) - $(i*t) cm" for i in 1:length(L)], title="total volume", yunit=u"mm^3", kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("V$n.pdf")
#     plot(r, :time, [Symbol("C$(i-1)") for i in 1:length(L)], legend="soil depth", names=["$((i-1)*t) - $(i*t) cm" for i in 1:length(L)], title="count", kind=:line, backend=:Gadfly)[] |> Cropbox.Gadfly.PDF("C$n.pdf")
#     CropRootBox.writevtk(n, s)
# end


text:
using DesignAutomation

@testset "Hooke & Jeeves" begin
    fn(x) = (x[1]-5)^2+2*(x[2]-5)^2
    x0 = [0., 0.]
    res = hooke_jeeves(fn, x0)

    @test res.minimizer == [5.000000000000004, 5.000000000000004]
    @test res.minimum == 5.916456789157589e-29
    @test res.iterations == 37
    @test res.converged == true

    fn(x) = (2*x[2]-3*x[1]^2)^4 + (5-2*x[1])^2
    x0 = [0., 0.]
    res = hooke_jeeves(fn, x0)

    @test res.minimizer == [2.4999999999999996, 9.37500000000003]
    @test res.minimum == 7.888609052210118e-31
    @test res.iterations == 99
    @test res.converged == true

    x0 = [0., 0.]
    res = hooke_jeeves(fn, x0; r=.5, h0=.5, max_iter=24)

    @test res.minimizer == [1.5, 3.25]
    @test res.minimum == 4.00390625
    @test res.iterations == 24
    @test res.converged == false

    x0 = [2., 9.]
    res = hooke_jeeves(fn, x0; r=.5, h0=.5, max_iter=24)

    @test res.minimizer == [2.5, 9.25]
    @test res.minimum == 0.00390625
    @test res.iterations == 24
    @test res.converged == false
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
# attribute
const EMPTY_ATTR = Dict{String, String}()
AttrType = Union{Dict{String, String}, Nothing}

abstract type CompoExpr end
abstract type SExpr <: CompoExpr end
abstract type CmdExpr <: SExpr end
abstract type GrpExpr <: SExpr end

struct ID <: SExpr
    val::String
    attr::AttrType
end
ID(x::String) = ID(x, nothing)

struct NUM <: SExpr
    val::String
    attr::AttrType
end
NUM(x::String) = NUM(x, nothing)

struct OPID <: SExpr
    val::Char
    attr::AttrType
end
OPID(x::String, attr) = OPID(x[1], attr)
OPID(x::String) = OPID(x, nothing)
OPID(x::Char) = OPID(x, nothing)

struct BIGOP <: SExpr
    val::Char
    size::Int
    attr::AttrType
end
BIGOP(x::Char) = BIGOP(x, 3, nothing)

struct SPACE <: SExpr
    val::String
    attr::AttrType
end
SPACE(x::String) = SPACE(x, nothing)

struct Group <: SExpr
    val::Vector{CompoExpr}
    attr::AttrType
end
Group(x::Vector{CompoExpr}) = Group(x, nothing)
function Group(::Int, args...)
    v = CompoExpr[]
    for i in args; push!(v, i); end
    return Group(v)
end

struct Environment <: SExpr
    tag::String
    val::Vector{CompoExpr}
    attr::AttrType
end
Environment(x::String, v::Vector{CompoExpr}) = Environment(x, v, nothing)

struct Fenced <: SExpr
    ldelim::Char
    rdelim::Char
    val::Vector{CompoExpr}
    attr::AttrType
end
Fenced(l::Char, r::Char, v::Vector{CompoExpr}) = Fenced(l, r, v, nothing)

# TeX built-in operators are converted into specific expression types
struct Newline <: CmdExpr
end

struct Ampersand <: CmdExpr
end

struct Sqrt <: CmdExpr
    base::SExpr
    attr::AttrType
end
Sqrt(b::SExpr) = Sqrt(b, nothing)

struct Frac <: CmdExpr
    num::SExpr
    den::SExpr
    attr::AttrType
end
Frac(n::SExpr, d::SExpr) = Frac(n, d, nothing)

struct Binom <: CmdExpr
    n::SExpr
    k::SExpr
    attr::AttrType
end
Binom(n::SExpr, k::SExpr) = Binom(n, k, nothing)

# composite expressions
struct Sub <: CompoExpr
    base::SExpr
    sub::SExpr
    attr::AttrType
end
Sub(b::SExpr, s::SExpr) = Sub(b, s, nothing)

struct Sup <: CompoExpr
    base::SExpr
    sup::SExpr
    attr::AttrType
end
Sup(b::SExpr, s::SExpr) = Sup(b, s, nothing)

struct SubSup <: CompoExpr
    base::SExpr
    sub::SExpr
    sup::SExpr
    attr::AttrType
end
SubSup(b::SExpr, sb::SExpr, Sp::SExpr) = SubSup(b, sb, sp, nothing)

struct Under <: CompoExpr
    base::SExpr
    under::SExpr
    attr::AttrType
end
Under(b::SExpr, u::SExpr) = Under(b, u, nothing)

struct Over <: CompoExpr
    base::SExpr
    over::SExpr
    attr::AttrType
end
Over(b::SExpr, o::SExpr) = Over(b, o, nothing)

struct UnderOver <: CompoExpr
    base::SExpr
    under::SExpr
    over::SExpr
    attr::AttrType
end
UnderOver(b::SExpr, u::SExpr, o::SExpr) = UnderOver(b, u, o, nothing)

# return a list of expressions
parse!(ts::TokenStream) = parse_expr_group!(ts, nothing)

# This function mainly deals with the sub/sup (`_` and `^`) operators.
# An expression that can be decomposed into simpler expressions connected
# by`_` and/or `^` is called a composite expression.
#
# We assume the input `ts` is not empty.
function parse_composite_expr!(ts::TokenStream)
    sexpr = parse_single_expr!(ts)

    # parse_single_expr!() could return nothing, example: {}
    if sexpr === nothing
        if peek(ts).val == "_"
            error("subscript _ behind an empty expression {}")
        elseif peek(ts).val == "^"
            error("superscript ^ behind an empty expression {}")
        elseif peek(ts).type == TK_FRACOVER
            error("\\over behind an empty expression {}")
        else
            return nothing
        end
    end
    peek(ts) === nothing && return sexpr
    # process composite expression. look-ahead
    if peek(ts).val == "_"
        # subscript
        get_token!(ts)
        peek(ts) == nothing && error("Illegal! expression ends with _")
        sub = parse_single_expr!(ts)
        if peek(ts) !== nothing && peek(ts).val == "^"
            get_token!(ts)
            peek(ts) == nothing && error("Illegal! expression ends with ^")
            sup = parse_single_expr!(ts)
            expr = make_sub_sup(sexpr, sub, sup)
        else
            expr = make_sub(sexpr, sub)
        end
    elseif peek(ts).val == "^"
        # superscript
        get_token!(ts)
        peek(ts) == nothing && error("Illegal! expression ends with ^")
        sup = parse_single_expr!(ts)
        if peek(ts) !== nothing && peek(ts).val == "_"
            get_token!(ts)
            peek(ts) == nothing && error("Illegal! expression ends with _")
            sub = parse_single_expr!(ts)
            expr = make_sub_sup(sexpr, sub, sup)
        else
            expr = make_sup(sexpr, sup)
        end
    elseif peek(ts).type == TK_FRACOVER
        # TODO: I should implement \over operator here (used by Maxima).
        # just implement the simplified version: Assume it is the infix
        # operator of two single-expressions
        #
        get_token!(ts)
        peek(ts) == nothing && error("Illegal! expression ends with \\over")
        den = parse_single_expr!(ts)
        expr = Frac(sexpr, den)
    else
        expr = sexpr
    end
    return expr
end

const token_type_to_expr_type = Dict(
# TODO: TK_LIM, TK_INT, TK_BIGOP should be processed specially, not
# just be treated as OPIDs.
    TK_ID     => ID,
    TK_NUM    => NUM,
    TK_OPID   => OPID,
    TK_PAREN  => OPID,
    TK_LIM    => ID,
    TK_INT    => OPID,
    TK_BIGOP  => OPID,
    TK_SPACE  => SPACE
    )

function parse_single_expr!(ts::TokenStream)
    t = get_token!(ts)

    if t.type in keys(token_type_to_expr_type)
        return token_type_to_expr_type[t.type](t.val)
    elseif t.type == TK_LBRACE
        v = parse_expr_group!(ts, TK_RBRACE)
        isempty(v) && return nothing
        # TODO: maybe I could do this optimization, will look at it later
        #length(v) == 1 && return v[1]
        return Group(v)
    elseif t.type == TK_SQRT
        peek(ts) === nothing && return error("unfinished \\sqrt command")
        b = parse_single_expr!(ts)
        b === nothing && return error("unfinished \\sqrt command")
        return Sqrt(b)
    elseif t.type == TK_FRAC
        #TODO: simplify this branch and the next branch
        peek(ts) === nothing && return error("unfinished \\frac command")
        num =  parse_single_expr!(ts)
        num === nothing && return error("unfinished \\frac command")
        peek(ts) === nothing && return error("unfinished \\frac command")
        den =  parse_single_expr!(ts)
        den === nothing && return error("unfinished \\frac command")
        return Frac(num, den)
    elseif t.type == TK_BINOM
        peek(ts) === nothing && return error("unfinished \\binom command")
        n =  parse_single_expr!(ts)
        n === nothing && return error("unfinished \\binom command")
        peek(ts) === nothing && return error("unfinished \\binom command")
        k =  parse_single_expr!(ts)
        k === nothing && return error("unfinished \\binom command")
        return Binom(n, k)
    elseif t.type == TK_AMP
        return Ampersand()
    elseif t.type == TK_NEWLINE
        return Newline()
    elseif t.type == TK_LEFT
        ldelimiter = get_token!(ts)
        v = parse_expr_group!(ts, TK_RIGHT)
        rdelimiter = get_token!(ts)
        return Fenced(ldelimiter.val[1], rdelimiter.val[1], v)
    elseif t.type == TK_BEGIN
        # t must be followed by a TK_ENV_NAME and optionally a TK_ENV_ATTR
        t1 = get_token!(ts)
        t1.type != TK_ENV_NAME && error("Expect a name after \\begin")
        # TODO: add support for TK_ENV_ATTR
        v = parse_expr_group!(ts, TK_END)
        # collect the envionment name right behind the the \end command
        peek(ts) === nothing && error("Expect an \\end tag, but stream ends")
        t2 = get_token!(ts)
        t2.type != TK_ENV_NAME && error("Expect a name after \\end")
        # verify the two tags following \begin and \end match
        t1.val != t2.val && error("\\begin{tag} and \\end{tag} mismatch!")
        return Environment(t1.val, v)
    elseif t.type == TK_ENV  #\bmatrix{}, ...
        # environments like \matrix{}, \pmatrix{}, ...
        tag = t.val
        t = get_token!(ts)
        t.type != TK_LBRACE && error("{ must follows a \\" + tag + " environment")
        return Environment(tag, parse_expr_group!(ts, TK_RBRACE))
    elseif t.type == TK_STYLE
        # Our handling of styling commands (\mathbb, \mathcal, \mathfrak, ...)
        # is simple. We assume there is a single group expression {...}
        # following the styling command, and inside the group expression, there
        # are only letters, each to be translated into a UCS codepoint that
        # represents a letter-variant, e.g, \mathbb{AC} should be interpreteded
        # as translations: A -> 𝔸 , C -> ℂ . The translation result is a
        # mathml group (<mrow>) that contains single letter IDs (<mi>).
        #
        # When we encounter a TK_STYLE, we translate all letters in the
        # proceeding {}. We do not do translation in scanning because there
        # may be a few optmizations that can be done here, so better keep
        # the scanning process simple.
        letter_variant = t.val   # variant can be mathbb, mathcal, ...
        t = get_token!(ts)
        t.type != TK_LBRACE &&
            error("\\", letter_variant, " command must be followed by {")
        v = CompoExpr[]
        while peek(ts) !== nothing
            t = get_token!(ts)
            t.type == TK_RBRACE && break
            # FIXME: what if t.val contains more than one element?
            !isletter(t.val[1]) && error("Incorrect input: ", t.val)
            d = LETTER_TRANSLATION[letter_variant]
            push!(v, ID(d[t.val]))
        end
        !(t.type == TK_RBRACE) && error("Incorrect \\", t.val, " command syntax")
        return (length(v) == 1 ? v[1] : Group(v))
    elseif t.type == TK_IGSTYLE
        #TODO: for now we ignore in-group styles (\bf, \it, ...)
    elseif t.type == TK_FUN
        return ID(t.val)
    else
        # Unknown tokens are treated as identifiers
        return ID(t.val, nothing)
    end
end

# Note that this function assumes the beginning token (e.g., \left, {, \begin) of
# the expression group in the input token stream `ts` is already parsed.
function parse_expr_group!(ts::TokenStream, end_tk_type::Union{Int, Nothing})
    v = CompoExpr[]

    while peek(ts) !== nothing && peek(ts).type != end_tk_type
        node = parse_composite_expr!(ts)
        node != nothing && push!(v, node)
    end

    ENDTYPE = Dict(
        TK_END    => "\\end",
        TK_RBRACE => "}",
        TK_RIGHT  => "\\right"
    )
    if peek(ts) == nothing && end_tk_type != nothing
        error("right delimiter ", ENDTYPE[end_tk_type], " is missing!")
    end

    if peek(ts) != nothing && peek(ts).type == end_tk_type
        # remove the ending token (\right, }, \end) from the token stream
        get_token!(ts)
    end
    return v
end

function make_sub(sexpr, sub)
    return Sub(sexpr, sub, nothing)
end

function make_sup(sexpr, sup)
    return Sup(sexpr, sup, nothing)
end

function make_sub_sup(sexpr, sub, sup)
    return SubSup(sexpr, sub, sup, nothing)
end

######################################################################################
# A simple nodes printer: `print_expr()`.
# This function prints a vector of `CompoExpr`s or one CompoExpr to screen.
const tabsize = 3
function print_expr(nodes::Vector{CompoExpr}, indent::Int)
    for n in nodes
        print_expr(n, indent)
    end
end
print_expr(nodes::Vector{CompoExpr}) = print_expr(nodes, 0)

function print_expr(n::Sub, indent::Int)
    println(" "^indent*"<Sub>: ")
    print_2childs(n.base, n.sub, indent + tabsize)
end
print_expr(n::Sub) = print_expr(n, 0)

function print_expr(n::Sup, indent::Int)
    println(" "^indent*"<Sup>: ")
    print_2childs(n.base, n.sup, indent + tabsize)
end
print_expr(n::Sup) = print_expr(n, 0)

function print_expr(n::SubSup, indent::Int)
    println(" "^indent*"<SubSup>: ")
    print_3childs(n.base, n.sub, n.sup, indent + tabsize)
end
print_expr(n::SubSup) = print_expr(n, 0)

function print_expr(n::Under, indent::Int)
    println(" "^indent*"<Under>: ")
    print_2childs(n.base, n.under, indent + tabsize)
end
print_expr(n::Under) = print_expr(n, 0)

function print_expr(n::Over, indent::Int)
    println(" "^indent*"<Over>: ")
    print_2childs(n.base, n.over, indent + tabsize)
end
print_expr(n::Over) = print_expr(n, 0)

function print_expr(n::UnderOver, indent::Int)
    println(" "^indent*"<UnderOver>: ")
    print_3childs(n.base, n.under, n.over, indent + tabsize)
end
print_expr(n::UnderOver) = print_expr(n, 0)

function print_expr(n::Sqrt, i::Int)
    println(" "^i*"<Sqrt>: ")
    print_expr(n.base, i + 2)
end
print_expr(n::Sqrt) = print_expr(n, 0)

function print_expr(n::Frac, i::Int)
    println(" "^i*"<Frac>: ")
    print_2childs(n.num, n.den, i + tabsize)
end
print_expr(n::Frac) = print_expr(n, 0)

function print_expr(n::Binom, i::Int)
    println(" "^i*"<Binom>: ")
    print_2childs(n.n, n.k, i + tabsize)
end
print_expr(n::Binom) = print_expr(n, 0)

function print_expr(n::Group, i::Int)
    println(" "^i*"<Group>: ")
    print_expr(n.val, i + tabsize)
end
print_expr(n::Group) = print_expr(n, 0)

function print_expr(n::Environment, i::Int)
    println(" "^i*"<Environment>: " * n.tag)
    print_expr(n.val, i + tabsize)
end
print_expr(n::Environment) = print_expr(n, 0)

function print_expr(n::Fenced, i::Int)
    println(" "^i * "<Fenced>: " * n.ldelim * " , " * n.rdelim)
    print_expr(n.val, i + tabsize)
end
print_expr(n::Fenced) = print_expr(n, 0)

print_expr(n::ID, i::Int) = println(" "^i*"<ID>: ", n.val)
print_expr(n::NUM, i::Int) = println(" "^i*"<NUM>: ", n.val)
print_expr(n::OPID, i::Int) = println(" "^i*"<OPID>: ", n.val)
print_expr(n::SPACE, i::Int) = println(" "^i*"<SPACE>: ", n.val)
print_expr(n::Ampersand, i::Int) = println(" "^i*"<Ampersand>")
print_expr(n::Newline, i::Int) = println(" "^i*"<Newline>")
print_expr(n, i::Int) = println(" "^i*"fallback: "*repr(n))

print_expr(n::ID) = print_expr(n, 0)
print_expr(n::NUM) = print_expr(n, 0)
print_expr(n::OPID) = print_expr(n, 0)
print_expr(n::SPACE) = print_expr(n, 0)
print_expr(n::Ampersand) = print_expr(n, 0)
print_expr(n::Newline) = print_expr(n, 0)

function print_2childs(c1, c2, indent)
    print_expr(c1, indent)
    print_expr(c2, indent)
end

function print_3childs(c1, c2, c3, indent)
    print_expr(c1, indent)
    print_expr(c2, indent)
    print_expr(c3, indent)
end


text:
using MosaicViews
using Images, TestImages

function calculate_averages(img::AbstractArray{T, N}, H𝚽::AbstractArray{T, M}) where {T<:Real, N, M}
    H𝚽ⁱ = @. 1. - H𝚽
    ∫H𝚽 = sum(H𝚽)
    ∫H𝚽ⁱ = sum(H𝚽ⁱ)
    if ndims(img) == 2
        ∫u₀H𝚽 = sum(img .* H𝚽)
        ∫u₀H𝚽ⁱ = sum(img .* H𝚽ⁱ)
    elseif ndims(img) == 3
        ∫u₀H𝚽 = sum(img .* H𝚽, dims=(1, 2))
        ∫u₀H𝚽ⁱ = sum(img .* H𝚽ⁱ, dims=(1, 2))
    end
    if ∫H𝚽 != 0
        c₁ = ∫u₀H𝚽 / ∫H𝚽
    end
    if ∫H𝚽ⁱ != 0
        c₂ = ∫u₀H𝚽ⁱ / ∫H𝚽ⁱ
    end

    return c₁, c₂
end

function difference_from_average_term(img::AbstractArray{T, N}, H𝚽::AbstractArray{T, M}, λ₁::Float64, λ₂::Float64) where {T<:Real, N, M}
    c₁, c₂ = calculate_averages(img, H𝚽)

    if ndims(img) == 2
        return @. -λ₁ * (img - c₁)^2 + λ₂ * (img - c₂)^2
    elseif ndims(img) == 3
        return -λ₁ .* sum((img .- c₁).^2, dims=3) .+ λ₂ .* sum((img .- c₂).^2, dims=3)
    end
end

function Hₕ(x::AbstractArray{T,N}, h::Float64=1.0) where {T<:Real, N}
    return @. 1. / 2. * (1. + 2. / pi * atan(x / h))
end

function δₕ(x::AbstractArray{T,N}, h::Float64=1.0) where {T<:Real, N}
    return @. h / (h^2 + x^2)
end

function initial_level_set(shape::Tuple)
    x₀ = reshape(collect(0:shape[begin]-1), shape[begin], 1)
    y₀ = reshape(collect(0:shape[begin+1]-1), 1, shape[begin+1])
    𝚽₀ = @. sin(pi / 5 * x₀) * sin(pi / 5 * y₀)
end

function calculate_variation(img::AbstractArray{T, N}, 𝚽ⁿ::AbstractArray{T, M}, μ::Float64, λ₁::Float64, λ₂::Float64, Δt::Float64) where {T<:Real, N, M}
    ϵ = 1e-16
    𝚽⁺ = padarray(𝚽ⁿ, Pad(1, 1))
    r = axes(𝚽⁺)
    @views 𝚽ᵢ₊ = 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 2:last(r[2])] .- 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1]
    @views 𝚽ᵢ₋ = 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1] .- 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]):last(r[2]) - 2]
    @views 𝚽ᵢ  = (𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 2:last(r[2])] .- 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]):last(r[2]) - 2]) / 2.

    @views 𝚽ⱼ₊ = 𝚽⁺[first(r[1]) + 2:last(r[1]), first(r[2]) + 1:last(r[2]) - 1] .- 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1]
    @views 𝚽ⱼ₋ = 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1] .- 𝚽⁺[first(r[1]):last(r[1]) - 2, first(r[2]) + 1:last(r[2]) - 1]
    @views 𝚽ⱼ  = (𝚽⁺[first(r[1]) + 2:last(r[1]), first(r[2]) + 1:last(r[2]) - 1] .- 𝚽⁺[first(r[1]):last(r[1]) - 2, first(r[2]) + 1:last(r[2]) - 1]) / 2.

    C₁ = @. 1. / sqrt(ϵ + 𝚽ᵢ₊^2 + 𝚽ⱼ^2)
    C₂ = @. 1. / sqrt(ϵ + 𝚽ᵢ₋^2 + 𝚽ⱼ^2)
    C₃ = @. 1. / sqrt(ϵ + 𝚽ᵢ^2 + 𝚽ⱼ₊^2)
    C₄ = @. 1. / sqrt(ϵ + 𝚽ᵢ^2 + 𝚽ⱼ₋^2)

    @views K = 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 2:last(r[2])] .* C₁ .+
               𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]):last(r[2]) - 2] .* C₂ .+
               𝚽⁺[first(r[1]) + 2:last(r[1]), first(r[2]) + 1:last(r[2]) - 1] .* C₃ .+
               𝚽⁺[first(r[1]):last(r[1]) - 2, first(r[2]) + 1:last(r[2]) - 1] .* C₄
    H𝚽 = @. 1. * (𝚽ⁿ > 0)

    diff = difference_from_average_term(img, H𝚽, λ₁, λ₂)
    𝚽ⁿ⁺¹ = (𝚽ⁿ .+ Δt .* δₕ(𝚽ⁿ) .* (μ * K .+ diff)) ./ (1 .+ μ .* Δt .* δₕ(𝚽ⁿ) .* (C₁ .+ C₂ .+ C₃ .+ C₄))
end

function calculate_reinitial(𝚽::AbstractArray{T, M}, Δt::Float64) where {T<:Real, M}
    ϵ = 1e-8
    𝚽⁺ = padarray(𝚽, Pad(1, 1))
    r = axes(𝚽⁺)
    @views a = @. 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1] - 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]):last(r[2]) - 2]
    @views b = @. 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 2:last(r[2])] - 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1]
    @views c = @. 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1] - 𝚽⁺[first(r[1]):last(r[1]) - 2, first(r[2]) + 1:last(r[2]) - 1]
    @views d = @. 𝚽⁺[first(r[1]) + 2:last(r[1]), first(r[2]) + 1:last(r[2]) - 1] - 𝚽⁺[first(r[1]) + 1:last(r[1]) - 1, first(r[2]) + 1:last(r[2]) - 1]

    a⁺ = max.(a, 0)
    a⁻ = min.(a, 0)
    b⁺ = max.(b, 0)
    b⁻ = min.(b, 0)
    c⁺ = max.(c, 0)
    c⁻ = min.(c, 0)
    d⁺ = max.(d, 0)
    d⁻ = min.(d, 0)

    G = zeros(size(𝚽))
    index⁺ = 𝚽 .> 0
    index⁻ = 𝚽 .< 0
    @. G = (sqrt(max(a⁺^2, b⁻^2) + max(c⁺^2, d⁻^2)) - 1) * index⁺ + (sqrt(max(a⁻^2, b⁺^2) + max(c⁻^2, d⁺^2)) - 1) * index⁻
    sign𝚽 = @. 𝚽 / sqrt(𝚽^2 + ϵ)
    𝚿 = @. 𝚽 - Δt * sign𝚽 * G
end

function reinitialize(𝚽::AbstractArray{T, M}, Δt::Float64, max_reiter::Int64=5) where {T<:Real, M}
    iter = 0
    while iter < max_reiter
        𝚽 .= calculate_reinitial(𝚽, Δt)
        iter += 1
    end

    return 𝚽
end

function chan_vese(img::AbstractArray{T,N};
                   μ::Float64=0.25,
                   λ₁::Float64=1.0,
                   λ₂::Float64=1.0,
                   tol::Float64=1e-3,
                   max_iter::Int64=500,
                   Δt::Float64=0.5,
                   reinitial_flag::Bool=false) where {T<:Real, N}
    iter = 0
    D = ndims(img)
    if D == 3
        img = PermutedDimsArray(img, (2, 3, 1))
    end
    𝚽ⁿ = initial_level_set(size(img))
    δ = tol + 1
    img .= img .- minimum(img)

    if maximum(img) != 0
        img .= img ./ maximum(img)
    end

    while (δ > tol) & (iter < max_iter)
        𝚽ⁿ⁺¹ = calculate_variation(img, 𝚽ⁿ, μ, λ₁, λ₂, Δt)
        δ = sqrt(meanfinite((𝚽ⁿ⁺¹ .- 𝚽ⁿ).^2, (1, 2))[1])
        if reinitial_flag
            𝚽ⁿ .= reinitialize(𝚽ⁿ⁺¹, Δt)
        else
            if D == 3
                r = axes(𝚽ⁿ⁺¹)
                @views 𝚽ⁿ .= 𝚽ⁿ⁺¹[:, :, first(r[3])]
            else
                @views 𝚽ⁿ .= 𝚽ⁿ⁺¹
            end
        end

        iter += 1
    end

    return 𝚽ⁿ, iter
end

img = float64.(channelview(testimage("cameraman")))
# img = float64.(channelview(testimage("lena_color_512")))
𝚽, iter_num = chan_vese(img, μ=0.25, λ₁=1.0, λ₂=1.0, tol=1e-3, max_iter=200, Δt=0.5, reinitial_flag=false)

segmentation = 𝚽 .> 0
print(iter_num)
𝚽 .= 𝚽 .- minimum(𝚽)

if maximum(𝚽) != 0
    𝚽 .= 𝚽 ./ maximum(𝚽)
end

output = mosaicview(img, segmentation, 𝚽; nrow=1, ncol=3, rowmajor=true)
save("demo.png", output)


text:
using JLD, Plots, Statistics

pyplot() # plotting backend

include("plot_defs.jl")

# plot the data for the simple periodic example

dirname = @__DIR__
fname = dirname * "/results/simple_periodic_pg_out.jld"
file = jldopen(fname,"r")

errs = read(file["errs"])
params = read(file["params"])

sigmas = params[1,:]

mederrs = dropdims(median(errs,dims=2),dims=2)

mederrs_l2 = mederrs[1,:]
mederrs_huber = mederrs[2,:]
mederrs_exact = mederrs[3,:]

ph = plot(sigmas, mederrs_exact, label="Exact DMD")
plot!(sigmas, mederrs_l2, label="Optimized DMD")
plot!(sigmas, mederrs_huber, label="Robust DMD (Huber)")
xaxis!("background noise",:log10)
yaxis!("error in recovered eigenvalues",:log10)
title!("Simple Periodic Example")

#  save figure
mkpath(dirname * "/figures")
fname = dirname * "/figures/simple_periodic_fig.pdf"
savefig(ph,fname)


text:
#
# Graph Functionality
#

function _show_conns(io, m, comp_name, which::Symbol)
    datumtype = which == :incoming ? "parameters" : "variables"
    println(io, "   $which $datumtype:")

    conns = get_connections(m, comp_name, which)

    if length(conns) == 0
        println(io, "      none")
    else
        for conn in conns
            if which == :incoming
                println(io, "      - $(conn.src_comp_path).$(conn.dst_par_name)")
            else
                println(io, "      - $(conn.dst_comp_path).$(conn.src_var_name)")
            end
        end
    end
end

show_conns(m::Model) = show_conns(stdout, m)

function show_conns(io::IO, m::Model)
    println(io, "Model component connections:")

    for (i, comp_name) in enumerate(compkeys(m.md))
        comp_def = compdef(m.md, comp_name)
        println(io, "$i. $(comp_def.comp_id) as :$(nameof(comp_def))")
        _show_conns(io, m, comp_name, :incoming)
        _show_conns(io, m, comp_name, :outgoing)
    end
end

function _filter_connections(conns::Vector{InternalParameterConnection}, comp_path::ComponentPath, which::Symbol)
    if which == :all
        f = obj -> (obj.src_comp_path == comp_path || obj.dst_comp_path == comp_path)
    elseif which == :incoming
        f = obj -> obj.dst_comp_path == comp_path
    elseif which == :outgoing
        f = obj -> obj.src_comp_path == comp_path
    else
        error("Invalid parameter for the 'which' argument; must be 'all' or 'incoming' or 'outgoing'.")
    end

    return collect(Iterators.filter(f, conns))
end

get_connections(m::Model, ci::LeafComponentInstance, which::Symbol) = get_connections(m, pathof(ci), which)

function get_connections(m::Model, cci::CompositeComponentInstance, which::Symbol)
    conns = []
    for ci in components(cci)
        append!(conns, get_connections(m, pathof(ci), which))
    end
    return conns
end

function get_connections(m::Model, comp_path::ComponentPath, which::Symbol)
    md = modeldef(m)
    return _filter_connections(internal_param_conns(md), comp_path, which)
end

function get_connections(m::Model, comp_name::Symbol, which::Symbol)
    comp = compdef(m, comp_name)
    get_connections(m, comp.comp_path, which)
end

function get_connections(mi::ModelInstance, comp_path::ComponentPath, which::Symbol)
    md = modeldef(mi)
    return _filter_connections(internal_param_conns(md), comp_path, which)
end


text:
"""
Julia package designed to build statistical models with associated systematic
uncertainties and provide inference into particular parameters of interest.
Both Bayesian and Frequentist analysis are supported.
"""
module Batman

using NLopt
using DataFrames
using DataStructures
using SpecialFunctions
using PyPlot
using Printf

# Logging
import Logging

function bat_format(level, _module, group, id, file, line)
  color, prefix, suffix = Logging.default_metafmt(level, _module, group, id, file, line)
  suffix = level == Logging.Debug ? "" : suffix
  return color, prefix, suffix
end

function bat_log(;kwargs...)
  kwargs = Dict(kwargs)
  debugLevel = haskey(ENV, "JULIA_DEBUG") ? Logging.Debug : Logging.Info
  debugLevel = get(kwargs, :level, debugLevel)
  Logging.ConsoleLogger( stdout, debugLevel; meta_formatter=bat_format )
end

Logging.global_logger( bat_log() )

# Primary definitions
include("predefined_functions.jl")
include("variable.jl")
include("likelihood.jl")
include("model.jl")
include("profile.jl")
include("plotter.jl")
include("show.jl")
# Kit definitions
include("modelkits/base.jl")
# Submodules
include("datastructures/DataStructures.jl")
include("toolbelt/Toolbelt.jl")
# Finally, macros
include("macros.jl")

# Batman exports based on a blacklist; where functions
# which begin with "_" are not exported.
const _EXCLUDE_SYMBOLS=[Symbol(@__MODULE__), :eval, :include]

for sym in names(@__MODULE__, all=true)
  sym_string = string(sym)
  if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_")
    continue
  end
  if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
                                  Base.isidentifier(sym_string[2:end])))
    continue
  end
  @eval export $sym
end

end # End module Batman


text:
# # Lasso, Ridge and Elastic Net Regressions
#
# This notebook presents a simple implementation of Lasso and elastic net regressions.

#-

# ## Load Packages and Extra Functions

using DelimitedFiles, LinearAlgebra, Statistics, Plots, Convex, SCS

import MathOptInterface
const MOI = MathOptInterface

# # Loading Data
#
# We use the diabetes data from Efron et al, downloaded from https://web.stanford.edu/~hastie/StatLearnSparsity_files/DATA/diabetes.html and then converted from a tab to a comma delimited file.
#
# All data series are standardised (see below) to have zero means and unit standard deviation, which improves the numerical stability. (Efron et al do not standardise the scale of the response variable.)

(x,header) = readdlm("aux_files/diabetes.csv",',',header=true)
#display(header)
#display(x)

x = (x .- mean(x,dims=1))./std(x,dims=1)          #standardise

(Y,X) = (x[:,end],x[:,1:end-1]);                  #to get traditional names
xNames = header[1:end-1];

# # Lasso, Ridge and Elastic Net Regressions
#
# (a)  The regression is $Y = Xb + u$,
# where $Y$ and $u$ are $T \times 1$, $X$ is $T \times K$, and $b$ is the $K$-vector of regression coefficients.
#
# (b) We want to minimize $(Y-Xb)'(Y-Xb) + \gamma \sum |b_i| + \lambda \sum b_i^2$.
#
# (c) We can equally well minimise $b'Qb - 2c'b + \gamma \sum |b_i| + \lambda \sum b_i^2$,
# where $Q = X'X$ and $c=X'Y$
#
# (d) Lasso: $\gamma>0,\lambda=0$; Ridge: $\gamma=0,\lambda>0$; elastic net: $\gamma>0,\lambda>0$.

"""
    LassoEN(Y,X,γ,λ)

Do Lasso (set γ>0,λ=0), ridge (set γ=0,λ>0) or elastic net regression (set γ>0,λ>0).

## Input
- `Y::Vector`:     T-vector with the response (dependent) variable
- `Y::VecOrMat`:   TxK matrix of covariates (regressors)
- `γ::Number`:     penalty on sum(abs.(b))
- `λ::Number`:     penalty on sum(b.^2)

"""
function LassoEN(Y,X,γ,λ=0.0)

    K = size(X,2)

    b_ls = X\Y                    #LS estimate of weights, no restrictions

    Q  = X'X
    c  = X'Y                      #c'b = Y'X*b

    b  = Variable(K)              #define variables to optimize over
    L1 = quadform(b,Q)            #b'Q*b
    L2 = dot(c,b)                 #c'b
    L3 = norm(b,1)                #sum(|b|)
    L4 = sumsquares(b)            #sum(b^2)

    Sol = minimize(L1-2*L2+γ*L3+λ*L4)      #u'u + γ*sum(|b|) + λsum(b^2), where u = Y-Xb
    solve!(Sol,()->SCS.Optimizer(verbose = false))
    Sol.status == MOI.OPTIMAL ? b_i = vec(evaluate(b)) : b_i = NaN

    return b_i, b_ls

end

# The next cell makes a Lasso regression for a single value of γ.

K = size(X,2)
γ = 100

(b,b_ls) = LassoEN(Y,X,γ)

println("OLS and Lasso coeffs (with γ=$γ)")
display([["" "OLS" "Lasso"];xNames b_ls b])

# # Redo the Lasso Regression with Different Gamma Values
#
#
# We now loop over $\gamma$ values.
#
# Remark: it would be quicker to put this loop inside the `LassoEN()` function so as to not recreate `L1`-`L4`.

nγ = 101
γM = range(0; stop=600, length=nγ)             #different γ values

bLasso   = fill(NaN,size(X,2),nγ)       #results for γM[i] are in bLasso[:,i]
for i = 1:nγ
    bLasso[:,i], = LassoEN(Y,X,γM[i])
end

#-

plot(log.(γM),bLasso',
     title = "Lasso regression coefficients",
     xlabel = "log(γ)",
     label = permutedims(xNames),
     size = (600,400))

# # Ridge Regression
#
# We use the same function to do a ridge regression. Alternatively, do `b = inv(X'X + λ*I)*X'Y`.

nλ = 101
λM = range(0; stop=3000, length=nλ)

bRidge  = fill(NaN,size(X,2),nλ)
for i = 1:nλ
    bRidge[:,i], = LassoEN(Y,X,0,λM[i])
end

#-

plot(log.(λM),bRidge',
     title = "Ridge regression coefficients",
     xlabel = "log(λ)",
     label = permutedims(xNames),
     size = (600,400))

# # Elastic Net Regression

λ = 200
println("redo the Lasso regression, but with λ=$λ: an elastic net regression")

bEN   = fill(NaN,size(X,2),nγ)
for i = 1:nγ
    bEN[:,i], = LassoEN(Y,X,γM[i],λ)
end

#-

plot(log.(γM),bEN',
     title = "Elastic Net regression coefficients",
     xlabel = "log(γ)",
     label = permutedims(xNames),
     size = (600,400))
#-


text:
using QPSReader
using Logging
using Printf

const RMPDIR = joinpath(@__DIR__, "rmp")
const RMPS = readdir(RMPDIR)

CG = Dict{Tuple{Int,Int}, Int}()
res = Dict{Tuple{Int,Int}, Dict}()

# Go through each RMP file to find number of CG iterations
for frmp in RMPS
    # File name has form DER_<T>_<R>_<iter>.mps
    fname, T, R, n = split(frmp, r"_|\.")[1:4]
    T = parse(Int, T)
    R = parse(Int, R)
    niter = parse(Int, n)

    # Check if such file already exists
    n_ = get(CG, (T, R), 0)
    CG[T, R] = max(niter, n_)
end

function parse_cg_log(flog)
    tmp = NaN
    tsp = NaN
    open(flog, "r") do f
        for ln in eachline(f)
            if occursin("Total time / MP", ln)
                tmp = parse(Float64, split(ln)[end][1:end-1])
                continue
            end
            if occursin("Total time / SP", ln)
                tsp = parse(Float64, split(ln)[end][1:end-1])
                continue
            end
        end
    end
    return tmp, tsp
end

# Go through the files again to access number of columns and linking constraints
for ((T, R), niter) in CG
    r_ = Dict()
    r_["cgiter"] = niter

    # Read
    frmp = "DER_$(T)_$(R)_$(niter).mps"
    qps = with_logger(Logging.NullLogger()) do
        readqps(joinpath(RMPDIR, frmp))
    end

    # Problem stats:
    #   * linking constraints
    #   * linking variables
    #   * Number of sub-problems
    #   * Density of lower blocks
    M = qps.ncon
    N = qps.nvar
    m0 = M - R

    nz = 0
    blockidx = zeros(Int, N)
    for (i, j, v) in zip(qps.arows, qps.acols, qps.avals)
        if 1 <= i <= m0
            nz += 1
        else
            @assert blockidx[j] == 0
            blockidx[j] = i - m0
        end
    end
    n0 = sum(blockidx .== 0)

    r_["m0"] = m0
    r_["n0"] = n0
    r_["R"] = R
    r_["N"] = N
    r_["density"] = nz / (m0 * N)

    # Parse log to get MP and SP total time
    flog = "cg_$(T)_$(R).log"
    tmp, tsp = parse_cg_log(joinpath(@__DIR__, "log", flog))

    r_["tmp"] = tmp
    r_["tsp"] = tsp

    @info frmp r_

    res[T, R] = r_
end

# Sort keys
KEYS = sort(collect(keys(res)))

# Print the table
for (T, R) in KEYS
    r = res[T, R]
    @printf "DER-%d & %5d & %3d & %5.1f & %5.1f & %3d & %6d & %.1f \$\\%%\$ \\\\\n" T R r["cgiter"] r["tmp"] r["tsp"] r["m0"]  r["N"] 100*r["density"]
end


text:
@testset "get_vertexval(g, v, ::Symbol)" begin

    g1 = ValGraph{Int8}(2, vertexval_types=(a=Int, b=Bool), vertexval_init=v -> (v, iseven(v)))
    g2 = ValOutDiGraph{Int16}(2, vertexval_types=(a=Int, b=Bool), vertexval_init=v -> (v, iseven(v)))
    g3 = ValDiGraph{UInt8}(2, vertexval_types=(a=Int, b=Bool), vertexval_init=v -> (v, iseven(v)))

    @test get_vertexval(g1, 1, :a) == 1
    @test get_vertexval(g1, Int8(2), :b) == true

    @test get_vertexval(g2, UInt8(1), :b) == false
    @test get_vertexval(g2, 2, :a) == 2

    @test get_vertexval(g3, 1, :b) == false
    @test get_vertexval(g3, 2, :a) == 2
end


text:
using Test
using DataTypesBasic

# Const
# =====

@test promote_type(Const{Int}, Const{Number}) == Const{Number}
@test promote_type(Const{Int}, Const) == Const
@test promote_type(Const, Const) == Const

@test Base.promote_typejoin(Const{Number}, Const{Number}) == Const{Number}
@test Base.promote_typejoin(Const{Int}, Const{Number}) == Const
@test Base.promote_typejoin(Const{Number}, Const{Int}) == Const
@test Base.promote_typejoin(Const, Const{Int}) == Const
@test Base.promote_typejoin(Const, Const) == Const

# Identity
# ========

@test promote_type(Identity{Int}, Identity{Number}) == Identity{Number}
@test promote_type(Identity{Int}, Identity) == Identity
@test promote_type(Identity, Identity) == Identity

@test Base.promote_typejoin(Identity{Number}, Identity{Number}) == Identity{Number}
@test Base.promote_typejoin(Identity{Int}, Identity{Number}) == Identity
@test Base.promote_typejoin(Identity{Number}, Identity{Int}) == Identity
@test Base.promote_typejoin(Identity{Int}, Identity) == Identity
@test Base.promote_typejoin(Identity, Identity{Int}) == Identity
@test Base.promote_typejoin(Identity, Identity) == Identity

# Option
# ======

@test promote_type(Identity{Int}, Const{Nothing}) == Union{Identity{Int}, Const{Nothing}}
@test promote_type(Const{Nothing}, Identity{<:Number}) == Union{Identity{<:Number}, Const{Nothing}}

@test Dict(:a => Identity(1), :b => Const(nothing)) isa Dict{Symbol, Option{Int}}

@test Base.promote_typejoin(Const{Nothing}, Identity{Number}) == Option{Number}
@test Base.promote_typejoin(Identity{Number}, Const{Nothing}) == Option{Number}
@test Base.promote_typejoin(Option{Number}, Identity{Int}) == Option
@test Base.promote_typejoin(Identity{Int}, Option{Number}) == Option
@test Base.promote_typejoin(Const{Nothing}, Option) == Option
@test Base.promote_typejoin(Option, Const{Nothing}) == Option
@test Base.promote_typejoin(Const{Nothing}, Identity{Number}) == Option{Number}
@test Base.promote_typejoin(Identity{Number}, Const{Nothing}) == Option{Number}
@test Base.promote_typejoin(Const{Nothing}, Identity) == Option
@test Base.promote_typejoin(Identity, Const{Nothing}) == Option

# Either
# ======

# Const & Identity
# ----------------

@test promote_type(Identity{String}, Const{Symbol}) == Union{Identity{String}, Const{Symbol}}
@test promote_type(Identity{AbstractString}, Const{Number}) == Union{Identity{AbstractString}, Const{Number}}

@test Base.promote_typejoin(Identity{Int}, Const{String}) == Either{String, Int}
@test Base.promote_typejoin(Const{String}, Identity{Int}) == Either{String, Int}
@test Base.promote_typejoin(Identity{Int}, Const) == Either{<:Any, Int}
@test Base.promote_typejoin(Const, Identity{Int}) == Either{<:Any, Int}
@test Base.promote_typejoin(Identity, Const{String}) == Either{String}
@test Base.promote_typejoin(Const{String}, Identity) == Either{String}
@test Base.promote_typejoin(Const, Identity) == Either
@test Base.promote_typejoin(Identity, Const) == Either

# Const & Either
# --------------

@test promote_type(Const{Int}, Either{Number}) == Either{Number}
@test promote_type(Const{Number}, Either{Int}) == Either{Number}
@test promote_type(Const{Number}, Either{Int}) == Either{Number}
@test promote_type(Const{Number}, Either{Int, String}) == Either{Number, String}
@test promote_type(Const{Int}, Either{Number, String}) == Either{Number, String}
@test promote_type(Const{Int}, Either{Number, String}) == Either{Number, String}
@test promote_type(Const{Int}, Either{<:Any, String}) == Either{<:Any, String}
@test promote_type(Const{Int}, Either) == Either

@test promote_type(Either{AbstractString}, Const{String}) == Either{AbstractString}
@test promote_type(Either{AbstractString}, Const) == Either
@test promote_type(Either{Integer}, Const{Number}) == Either{Number}
@test promote_type(Either{Number}, Const{Integer}) == Either{Number}
@test promote_type(Either{Number, String}, Const{Integer}) == Either{Number, String}

@test Base.promote_typejoin(Const, Either{Number, String}) == Either{<:Any, String}
@test Base.promote_typejoin(Either{Number, String}, Const) == Either{<:Any, String}

@test Base.promote_typejoin(Const{Number}, Either{Number, String}) == Either{Number, String}
@test Base.promote_typejoin(Either{Number, String}, Const{Number}) == Either{Number, String}

@test Base.promote_typejoin(Const{Int}, Either{Number, String}) == Either{<:Any, String}
@test Base.promote_typejoin(Either{Number, String}, Const{Int}) == Either{<:Any, String}

@test Base.promote_typejoin(Const, Either{Number}) == Either
@test Base.promote_typejoin(Either{Number}, Const) == Either
@test Base.promote_typejoin(Const{Int}, Either{Number}) == Either
@test Base.promote_typejoin(Either{Number}, Const{Int}) == Either

@test Base.promote_typejoin(Const, Either) == Either
@test Base.promote_typejoin(Either, Const) == Either

# Identity & Either
# -----------------

@test promote_type(Either{AbstractString}, Identity{Int}) == Either{AbstractString}
@test promote_type(Either{AbstractString}, Identity) == Either{AbstractString}
@test promote_type(Identity{Int}, Either{<:Any, Number}) == Either{<:Any, Number}
@test promote_type(Identity{Number}, Either{<:Any, Int}) == Either{<:Any, Number}
@test promote_type(Identity{Number}, Either{<:Any, Int}) == Either{<:Any, Number}
@test promote_type(Identity{Number}, Either{String, Int}) == Either{String, Number}
@test promote_type(Identity{Int}, Either{String, Number}) == Either{String, Number}
@test promote_type(Identity{Int}, Either{String, Number}) == Either{String, Number}
@test promote_type(Identity{Int}, Either{String, <:Any}) == Either{String, <:Any}
@test promote_type(Identity{Int}, Either) == Either

@test promote_type(Either{String, Integer}, Identity{Number}) == Either{String, Number}
@test promote_type(Either{String, Number}, Identity{Integer}) == Either{String, Number}
@test promote_type(Either{String, Integer}, Identity) == Either{String}
@test promote_type(Either{<:Any, Integer}, Identity) == Either
@test promote_type(Either{<:Any, Integer}, Identity{Number}) == Either{<:Any, Number}

@test Base.promote_typejoin(Identity, Either{String, Int}) == Either{String, <:Any}
@test Base.promote_typejoin(Either{String, Int}, Identity) == Either{String, <:Any}
@test Base.promote_typejoin(Identity{Int}, Either{String, Int}) == Either{String, Int}
@test Base.promote_typejoin(Either{String, Int}, Identity{Int}) == Either{String, Int}
@test Base.promote_typejoin(Identity{Number}, Either{String, Int}) == Either{String, <:Any}
@test Base.promote_typejoin(Either{String, Int}, Identity{Number}) == Either{String, <:Any}

@test Base.promote_typejoin(Identity, Either{<:Any, Int}) == Either
@test Base.promote_typejoin(Either{<:Any, Int}, Identity) == Either
@test Base.promote_typejoin(Identity{Int}, Either{<:Any, Int}) == Either{<:Any, Int}
@test Base.promote_typejoin(Either{<:Any, Int}, Identity{Int}) == Either{<:Any, Int}
@test Base.promote_typejoin(Identity{Number}, Either{<:Any, Int}) == Either
@test Base.promote_typejoin(Either{<:Any, Int}, Identity{Number}) == Either

@test Base.promote_typejoin(Identity, Either{Int, <:Any}) == Either{Int, <:Any}
@test Base.promote_typejoin(Either{Int, <:Any}, Identity) == Either{Int, <:Any}

@test Base.promote_typejoin(Identity, Either) == Either
@test Base.promote_typejoin(Either, Identity) == Either

# Either & Either
# ---------------

@test promote_type(Either{String, Integer}, Either{String, Integer}) == Either{String, Integer}
@test promote_type(Either{Integer, AbstractString}, Either{Integer, String}) == Either{Integer, AbstractString}
@test promote_type(Either{Integer, AbstractString}, Either{Number, String}) == Either{Number, AbstractString}
@test promote_type(Either{AbstractString, Integer}, Either{String, Integer}) == Either{AbstractString, Integer}
@test promote_type(Either{AbstractString, Integer}, Either{String, Number}) == Either{AbstractString, Number}
@test promote_type(Either{Number, String}, Either{Integer, AbstractString}) == Either{Number, AbstractString}

@test promote_type(Either{Integer, AbstractString}, Either{Integer, Bool}) == Either{Integer, <:Any}
@test promote_type(Either{Integer, AbstractString}, Either{Number, Bool}) == Either{Number, <:Any}
@test promote_type(Either{Integer, String}, Either{String, String}) == Either{<:Any, String}
@test promote_type(Either{Integer, AbstractString}, Either{String, String}) == Either{<:Any, AbstractString}

@test promote_type(Either{<:Any, Integer}, Either{<:Any, Integer}) == Either{<:Any, Integer}
@test promote_type(Either{<:Any, Integer}, Either{<:Any, Number}) == Either{<:Any, Number}
@test promote_type(Either{String, <:Any}, Either{String, <:Any}) == Either{String, <:Any}
@test promote_type(Either{String, <:Any}, Either{AbstractString, <:Any}) == Either{AbstractString, <:Any}
@test promote_type(Either, Either) == Either

@test Base.promote_typejoin(Either{Int, String}, Either{Int, String}) == Either{Int, String}

@test Base.promote_typejoin(Either{Int, String}, Either{Number, String}) == Either{<:Any, String}
@test Base.promote_typejoin(Either{Number, String}, Either{Int, String}) == Either{<:Any, String}
@test Base.promote_typejoin(Either{Number, String}, Either{Int, AbstractString}) == Either

@test Base.promote_typejoin(Either{String, Int}, Either{String, Number}) == Either{String, <:Any}
@test Base.promote_typejoin(Either{String, Number}, Either{String, Int}) == Either{String, <:Any}
@test Base.promote_typejoin(Either{String, Number}, Either{AbstractString, Int}) == Either

@test Base.promote_typejoin(Either{Int, <:Any}, Either{Int, <:Any}) == Either{Int, <:Any}
@test Base.promote_typejoin(Either{Int, <:Any}, Either{Number, <:Any}) == Either
@test Base.promote_typejoin(Either{Number, <:Any}, Either{Int, <:Any}) == Either

@test Base.promote_typejoin(Either{<:Any, Int}, Either{<:Any, Int}) == Either{<:Any, Int}
@test Base.promote_typejoin(Either{<:Any, Int}, Either{<:Any, Number}) == Either
@test Base.promote_typejoin(Either{<:Any, Number}, Either{<:Any, Int}) == Either

@test Base.promote_typejoin(Either{<:Any, Int}, Either) == Either
@test Base.promote_typejoin(Either, Either{<:Any, Int}) == Either
@test Base.promote_typejoin(Either{Int, <:Any}, Either) == Either
@test Base.promote_typejoin(Either, Either{Int, <:Any}) == Either
@test Base.promote_typejoin(Either, Either) == Either


text:
import HDF5

const HDF5_TS_ROOT_PATH = "time_series"

"""
Stores all time series data in an HDF5 file.

The file used is assumed to be temporary and will be automatically deleted when there are
no more references to the storage object.
"""
mutable struct Hdf5TimeSeriesStorage <: TimeSeriesStorage
    file_path::String
    read_only::Bool
end

"""
Constructs Hdf5TimeSeriesStorage by creating a temp file.
"""
function Hdf5TimeSeriesStorage()
    return Hdf5TimeSeriesStorage(true)
end

"""
Constructs Hdf5TimeSeriesStorage.

# Arguments
- `create_file::Bool`: create new file
- `filename=nothing`: if nothing, create a temp file, else use this name.
- `directory=nothing`: if set and filename is nothing, create a temp file in this
   directory. Use tempdir() if not set. This should be set if the time series data is larger
   than the tmp filesystem can hold.
- `read_only = false`: If true, don't allow changes to the file. Allows simultaneous read
   access.
"""
function Hdf5TimeSeriesStorage(
    create_file::Bool;
    filename = nothing,
    directory = nothing,
    read_only = false,
)
    if create_file
        if isnothing(filename)
            if isnothing(directory)
                directory = tempdir()
            end
            filename, io = mktemp(directory)
            close(io)
        end

        storage = Hdf5TimeSeriesStorage(filename, read_only)
        _make_file(storage)
    else
        storage = Hdf5TimeSeriesStorage(filename, read_only)
    end

    @debug "Constructed new Hdf5TimeSeriesStorage" storage.file_path read_only

    return storage
end

"""
Constructs Hdf5TimeSeriesStorage from an existing file.
"""
function from_file(
    ::Type{Hdf5TimeSeriesStorage},
    filename::AbstractString;
    read_only = false,
)
    if !isfile(filename)
        error("time series storage $filename does not exist")
    end

    if read_only
        file_path = abspath(filename)
    else
        file_path, io = mktemp()
        close(io)
        cp(filename, file_path; force = true)
    end
    storage = Hdf5TimeSeriesStorage(false; filename = file_path, read_only = read_only)
    @info "Loaded time series from storage file existing=$filename new=$(storage.file_path)"
    return storage
end

get_file_path(storage::Hdf5TimeSeriesStorage) = storage.file_path

function add_time_series!(
    storage::Hdf5TimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    label::AbstractString,
    ts::TimeSeriesData,
    columns = nothing,
)
    check_read_only(storage)
    uuid = string(get_uuid(ts))
    component_label = make_component_label(component_uuid, label)

    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        if !HDF5.exists(root, uuid)
            HDF5.g_create(root, uuid)
            path = root[uuid]
            @debug "Create new time series entry." uuid component_uuid label
            path["data"] = TimeSeries.values(ts.data)
            timestamps = [Dates.datetime2epochms(x) for x in TimeSeries.timestamp(ts.data)]
            path["timestamps"] = timestamps
            # Storing the UUID as an integer would take less space, but HDF5 library says
            # arrays of 128-bit integers aren't supported.
            path["components"] = [component_label]
            if !isnothing(columns)
                HDF5.attrs(path)["columns"] = string.(columns)
            end
        else
            path = root[uuid]
            if !isnothing(columns)
                if !HDF5.exists(HDF5.attrs(path), "columns")
                    throw(ArgumentError("columns are specified but existing array does not have columns"))
                end
                existing = Symbol.(HDF5.read(HDF5.attrs(path)["columns"]))
                if columns != existing
                    throw(ArgumentError("columns do not match $columns $existing"))
                end
            end
            @debug "Add reference to existing time series entry." uuid component_uuid label
            _append_item!(path, "components", component_label)
        end
    end
end

function add_time_series_reference!(
    storage::Hdf5TimeSeriesStorage,
    component_uuid::UUIDs.UUID,
    label::AbstractString,
    ts_uuid::UUIDs.UUID,
)
    check_read_only(storage)
    uuid = string(ts_uuid)
    component_label = make_component_label(component_uuid, label)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        path = root[uuid]
        _append_item!(path, "components", component_label)
        @debug "Add reference to existing time series entry." uuid component_uuid label
    end
end

function iterate_time_series(storage::Hdf5TimeSeriesStorage)
    Channel() do channel
        HDF5.h5open(storage.file_path, "r") do file
            root = _get_root(storage, file)
            for uuid_group in root
                uuid_path = HDF5.name(uuid_group)
                range = findlast("/", uuid_path)
                uuid_str = uuid_path[(range.start + 1):end]
                uuid = UUIDs.UUID(uuid_str)
                internal = InfrastructureSystemsInternal(uuid)
                time_series = TimeSeriesData(get_time_series(storage, uuid), internal)
                for item in HDF5.read(uuid_group["components"])
                    component, label = deserialize_component_label(item)
                    put!(channel, (component, label, time_series))
                end
            end
        end
    end
end

function remove_time_series!(
    storage::Hdf5TimeSeriesStorage,
    uuid::UUIDs.UUID,
    component_uuid::UUIDs.UUID,
    label::AbstractString,
)
    check_read_only(storage)
    HDF5.h5open(storage.file_path, "r+") do file
        root = _get_root(storage, file)
        path = _get_time_series_path(root, uuid)
        if _remove_item!(path, "components", make_component_label(component_uuid, label))
            @debug "$path has no more references; delete it."
            HDF5.o_delete(path)
        end
    end
end

function get_time_series(
    storage::Hdf5TimeSeriesStorage,
    uuid::UUIDs.UUID;
    index = 0,
    len = 0,
)::TimeSeries.TimeArray
    return HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        path = _get_time_series_path(root, uuid)

        if index == 0
            data = HDF5.read(path["data"])
            timestamps = HDF5.read(path["timestamps"])
        else
            @assert len != 0
            end_index = index + len - 1
            # HDF5.readmmap could be faster in many cases than this. However, experiments
            # resulted in various crashes if we tried to close the file before references
            # to the array data were garbage collected. May need to consult with the
            # Julia HDF5 library maintainers about that.
            sz_tuple = size(path["data"])
            if length(sz_tuple) == 1
                data = path["data"][index:end_index]
            else
                data = path["data"][index:end_index, 1:sz_tuple[2]]
            end
            timestamps = path["timestamps"][index:end_index]
        end

        # If the user set column names, return them. Otherwise, let TimeArray pick them.
        dtimestamps = [Dates.epochms2datetime(x) for x in timestamps]
        if HDF5.exists(HDF5.attrs(path), "columns")
            columns = Symbol.(HDF5.read(HDF5.attrs(path)["columns"]))
            return TimeSeries.TimeArray(dtimestamps, data, columns)
        end

        return TimeSeries.TimeArray(dtimestamps, data)
    end
end

function clear_time_series!(storage::Hdf5TimeSeriesStorage)
    check_read_only(storage)
    # Re-create the file. HDF5 will not actually free up the deleted space until h5repack
    # is run on the file.
    _make_file(storage)
    @info "Cleared all time series."
end

function get_num_time_series(storage::Hdf5TimeSeriesStorage)
    num = 0

    HDF5.h5open(storage.file_path, "r") do file
        root = _get_root(storage, file)
        for component in root
            num += 1
        end
    end

    return num
end

function _make_file(storage::Hdf5TimeSeriesStorage)
    HDF5.h5open(storage.file_path, "w") do file
        HDF5.g_create(file, HDF5_TS_ROOT_PATH)
    end
end

_get_root(storage::Hdf5TimeSeriesStorage, file) = file[HDF5_TS_ROOT_PATH]

function _get_time_series_path(root::HDF5.HDF5Group, uuid::UUIDs.UUID)
    uuid_str = string(uuid)
    if !HDF5.exists(root, uuid_str)
        throw(ArgumentError("UUID $uuid_str does not exist"))
    end

    return root[uuid_str]
end

function _append_item!(path::HDF5.HDF5Group, name::AbstractString, value::AbstractString)
    handle = HDF5.o_open(path, name)
    values = HDF5.read(handle)
    HDF5.close(handle)
    push!(values, value)

    ret = HDF5.o_delete(path, name)
    @assert ret == 0

    path[name] = values
    @debug "Appended $value to $name" values
end

"""
Removes value from the dataset called name.
Returns true if the array is empty afterwards.
"""
function _remove_item!(path::HDF5.HDF5Group, name::AbstractString, value::AbstractString)
    handle = HDF5.o_open(path, name)
    values = HDF5.read(handle)
    HDF5.close(handle)

    orig_len = length(values)
    filter!(x -> x != value, values)
    if length(values) != orig_len - 1
        throw(ArgumentError("$value wasn't stored in $name"))
    end

    ret = HDF5.o_delete(path, name)
    @assert ret == 0

    if isempty(values)
        is_empty = true
    else
        path[name] = values
        is_empty = false
    end
    @debug "Removed $value from $name" values

    return is_empty
end

function check_read_only(storage::Hdf5TimeSeriesStorage)
    if storage.read_only
        error("Operation not permitted; this time series file is read-only")
    end
end

function compare_values(x::Hdf5TimeSeriesStorage, y::Hdf5TimeSeriesStorage)::Bool
    data_x = sort!(collect(iterate_time_series(x)), by = z -> z[1])
    data_y = sort!(collect(iterate_time_series(y)), by = z -> z[1])
    if length(data_x) != length(data_y)
        @error "lengths don't match" length(data_x) length(data_y)
        return false
    end

    for ((uuid_x, label_x, ts_x), (uuid_y, label_y, ts_y)) in zip(data_x, data_y)
        if uuid_x != uuid_y
            @error "component UUIDs don't match" uuid_x uuid_y
            return false
        end
        if label_x != label_y
            @error "labels don't match" label_x label_y
            return false
        end
        if TimeSeries.timestamp(ts_x.data) != TimeSeries.timestamp(ts_y.data)
            @error "timestamps don't match" ts_x.data ts_y.data
            return false
        end
        if TimeSeries.values(ts_x.data) != TimeSeries.values(ts_y.data)
            @error "values don't match" ts_x.data ts_y.data
            return false
        end
    end

    return true
end


text:
"""
    fpcd!(rbm, data)

Trains the RBM on data using Persistent Contrastive divergence, with fast weights.
See http://dl.acm.org/citation.cfm?id=1553374.1553506.
"""
function fpcd!(rbm::RBM, data::AbstractArray;
    batchsize::Int = 1,
    epochs::Int = 1,
    optim = ADAM(),
    wts = nothing,
    steps::Int = 1,
    optimfast = ADAM(), # optimizer algorithm for fast parameters
    decayfast::Real = 19/20  # weight decay of fast parameters
)
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    @assert isnothing(wts) || _nobs(data) == _nobs(wts)
    stats = suffstats(rbm, data; wts)
    vm = sample_from_inputs(visible(rbm), falses(size(visible(rbm))..., batchsize))
    rbmfast = deepcopy(rbm) # store fast parameters
    # (Actually, the parameters of rbmfast are the sums θ_regular + θ_fast)
    for epoch in 1:epochs
        batches = minibatches(data, wts; batchsize = batchsize)
        Δt = @elapsed for (vd, wd) in batches
            vm = sample_v_from_v(rbmfast, vm; steps = steps)
            ∂ = ∂contrastive_divergence(rbm, vd, vm; wd, stats)
            update!(rbm, update!(∂, rbm, optim))
            # update fast parameters
            update!(rbmfast, update!(∂, rbmfast, optimfast))
            decayfast!(rbmfast, rbm; decay=decayfast)
        end
    end
    return rbm
end

#= Decays parameters of `rbmfast` towards those of `rbm`.
In other words, writing the parametrs of `rbmfast` as
ω_regular + ω_fast, where ω_regular are the prameters of
`rbm`, then here we decay the ω_fast part towards zero. =#
function decayfast!(rbmfast::M, rbm::M; decay::Real) where {M<:RBM}
    decayfast!(visible(rbmfast), visible(rbm); decay)
    decayfast!(hidden(rbmfast), hidden(rbm); decay)
    decayfast!(weights(rbmfast), weights(rbm); decay)
end

function decayfast!(fast::L, regular::L; decay::Real) where {L<:Union{Binary,Spin,Potts}}
    @assert size(fast) == size(regular)
    decayfast!(fast.θ, regular.θ; decay)
end

function decayfast!(fast::L, regular::L; decay::Real) where {L<:Union{Gaussian,ReLU}}
    @assert size(fast) == size(regular)
    decayfast!(fast.θ, regular.θ; decay)
    decayfast!(fast.γ, regular.γ; decay)
end

function decayfast!(fast::L, regular::L; decay::Real) where {L<:dReLU}
    @assert size(fast) == size(regular)
    decayfast!(fast.θp, regular.θp; decay)
    decayfast!(fast.θn, regular.θn; decay)
    decayfast!(fast.γp, regular.γp; decay)
    decayfast!(fast.γn, regular.γn; decay)
end

function decayfast!(ωfast::AbstractArray, ωregular::AbstractArray; decay::Real)
    @assert size(ωfast) == size(ωregular)
    ωfast .= decay .* ωfast + (1 - decay) .* ωregular
end


text:
using Bolt
using PyPlot
using Unitful, UnitfulAstro, NaturallyUnitful

par = CosmoParams()

# integrate saha to some transition redshift, we choose 1587.4
z_transition = 1587.4
x_transition = z2x(z_transition)
saha_z_grid = 1800:-10:z_transition
peebles_z_grid = z_transition:-10:1

early_time_Xₑ = Bolt.saha_Xₑ(par)
late_time_Xₑ = Bolt.peebles_Xₑ(
    par, early_time_Xₑ(x_transition), x_transition, 1.0)

clf()
plot(saha_z_grid, [early_time_Xₑ(z2x(z)) for z in saha_z_grid], "-", label="Saha")
plot(peebles_z_grid, [late_time_Xₑ(z2x(z)) for z in peebles_z_grid], "-", label="Peebles")
ylabel(raw"$X_e$")
xlim(1800,0)
ylim(2e-4, 2)
yscale("log")
xlabel(raw"$z$")
legend()
gcf()


text:
# Runs the StableDQMC.jl benchmarks.
# A branch name or commit id can be supplied as an argument.

using Pkg
Pkg.activate(@__DIR__)

using PkgBenchmark

branch_or_commit = nothing
if length(ARGS) > 0
    branch_or_commit = ARGS[1]
end

out = isnothing(branch_or_commit) ? "master" : string(branch_or_commit)

config = BenchmarkConfig(id = out,
                         juliacmd = `julia -O3 --project=.`,
                         env = Dict("JULIA_NUM_THREADS" => 1,
                                    "OMP_NUM_THREADS" => Sys.CPU_THREADS,
                                    "MKL_NUM_THREADS" => Sys.CPU_THREADS))

r = benchmarkpkg("StableDQMC", config;
             script = "benchmark/benchmarks.jl",
             resultfile = string(out, ".bench"))

export_markdown(string(out, ".md"), r)


text:
# ---
# title: 599. Minimum Index Sum of Two Lists
# id: problem599
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/minimum-index-sum-of-two-lists/description/>
# hidden: true
# ---
#
# Suppose Andy and Doris want to choose a restaurant for dinner, and they both
# have a list of favorite restaurants represented by strings.
#
# You need to help them find out their **common interest** with the **least list
# index sum**. If there is a choice tie between answers, output all of them with
# no order requirement. You could assume there always exists an answer.
#
#
#
# **Example 1:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["Piatti","The Grill at Torrey Pines","Hungry Hunter Steakhouse","Shogun"]
#     Output: ["Shogun"]
#     Explanation: The only restaurant they both like is "Shogun".
#
#
# **Example 2:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["KFC","Shogun","Burger King"]
#     Output: ["Shogun"]
#     Explanation: The restaurant they both like and have the least index sum is "Shogun" with index sum 1 (0+1).
#
#
# **Example 3:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["KFC","Burger King","Tapioca Express","Shogun"]
#     Output: ["KFC","Burger King","Tapioca Express","Shogun"]
#
#
# **Example 4:**
#
#
#
#     Input: list1 = ["Shogun","Tapioca Express","Burger King","KFC"], list2 = ["KNN","KFC","Burger King","Tapioca Express","Shogun"]
#     Output: ["KFC","Burger King","Tapioca Express","Shogun"]
#
#
# **Example 5:**
#
#
#
#     Input: list1 = ["KFC"], list2 = ["KFC"]
#     Output: ["KFC"]
#
#
#
#
# **Constraints:**
#
#   * `1 <= list1.length, list2.length <= 1000`
#   * `1 <= list1[i].length, list2[i].length <= 30`
#   * `list1[i]` and `list2[i]` consist of spaces `' '` and English letters.
#   * All the stings of `list1` are **unique**.
#   * All the stings of `list2` are **unique**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    Born-Mayer Potential
"""
mutable struct BornMayer <: Potential
    A::Float64
    ρ::Float64
end

"""
    BornMayer(params::Dict)

Creates a BM potential.
"""
function BornMayer(params::Dict)
    A = params["A"]
    ρ = params["ρ"]
    return BornMayer(A, ρ)
end

"""
    potential_energy(p::BornMayer, r::Position, args...)

Calculates Born-Mayer potential energy.
"""
function potential_energy(p::BornMayer, r::Position, args...)
    return p.A * exp(-norm(r) / p.ρ)
end