text:
"""
    get_split_pvals(com, ::Val{:Auto}, var::Variable)

Splits the possible values into two by using either :smallest or :biggest value and the rest.
It depends on whether it's a satisfiability or optimization problem and whether the variable has a positive
or negative coefficient + minimization or maximization
Return lb, leq, geq, ub => the bounds for the lower part and the bounds for the upper part
"""
function get_split_pvals(com, ::Val{:Auto}, var::Variable)
    @assert var.min != var.max
    if isa(com.objective, LinearCombinationObjective)
        linear_comb = com.objective.lc
        for i in 1:length(linear_comb.indices)
            if linear_comb.indices[i] == var.idx
                coeff = linear_comb.coeffs[i]
                factor = com.sense == MOI.MIN_SENSE ? -1 : 1
                if coeff * factor > 0
                    return get_split_pvals(com, Val(:Biggest), var)
                else
                    return get_split_pvals(com, Val(:Smallest), var)
                end
            end
        end
    end
    # fallback for satisfiability or not in objective
    return get_split_pvals(com, rand(CS_RNG, [Val(:Smallest), Val(:Biggest)]), var)
end

"""
    get_split_pvals(com, ::Val{:InHalf}, var::Variable)

Splits the possible values into two by obtaining the mean value.
Return lb, leq, geq, ub => the bounds for the lower part and the bounds for the upper part
"""
function get_split_pvals(com, ::Val{:InHalf}, var::Variable)
    pvals = values(var)
    mean_val = mean(pvals)
    leq = typemin(Int)
    geq = typemax(Int)
    lb = typemax(Int)
    ub = typemin(Int)
    for pval in pvals
        if pval <= mean_val && pval > leq
            leq = pval
        elseif pval > mean_val && pval < geq
            geq = pval
        end
        if pval < lb
            lb = pval
        end
        if pval > ub
            ub = pval
        end
    end
    return lb, leq, geq, ub
end

"""
    get_split_pvals(com, ::Val{:Random}, var::Variable)

Splits the possible values into two by using a random value for the left branch
**Attention:** The right branch is not useable atm.
"""
function get_split_pvals(com, ::Val{:Random}, var::Variable)
    @assert var.min != var.max
    val = rand(CS_RNG, values(var))
    return val, val, typemin(typeof(val)), typemin(typeof(val))
end

"""
    get_split_pvals(com, ::Val{:Smallest}, var::Variable)

Splits the possible values into two by using the smallest value and the rest.
Return lb, leq, geq, ub => the bounds for the lower part and the bounds for the upper part
"""
function get_split_pvals(com, ::Val{:Smallest}, var::Variable)
    @assert var.min != var.max
    right_lb = partialsort(values(var), 2)
    return var.min, var.min, right_lb, var.max
end

"""
    get_split_pvals(com, ::Val{:Biggest}, var::Variable)

Splits the possible values into two by using the biggest value and the rest.
Return lb, leq, geq, ub => the bounds for the lower part and the bounds for the upper part
"""
function get_split_pvals(com, ::Val{:Biggest}, var::Variable)
    @assert var.min != var.max
    left_ub = partialsort(values(var), 2; rev = true)
    return var.min, left_ub, var.max, var.max
end

"""
    get_next_branch_variable(com::CS.CoM)

Get the next weak index for backtracking. This will be the next branching variable.
Return whether there is an unfixed variable and a best index
"""
function get_next_branch_variable(com::CS.CoM)
    branch_var = get_next_branch_variable(com, com.branch_strategy)
    if !branch_var.is_feasible || branch_var.is_solution
        close_node!(com, com.c_backtrack_idx)
    end
    return branch_var
end

function get_next_branch_variable(com::CS.CoM, ::Val{:IMPS})
    lowest_num_pvals = typemax(Int)
    biggest_inf = -1
    best_vidx = -1
    biggest_dependent = typemax(Int)
    is_in_objective = false
    is_solution = true

    for vidx in 1:length(com.search_space)
        if !isfixed(com.search_space[vidx])
            num_pvals = nvalues(com.search_space[vidx])
            inf = com.bt_infeasible[vidx]
            if !is_in_objective && com.var_in_obj[vidx]
                is_in_objective = true
                lowest_num_pvals = num_pvals
                biggest_inf = inf
                best_vidx = vidx
                is_solution = false
                continue
            end
            if !is_in_objective || com.var_in_obj[vidx]
                if inf >= biggest_inf
                    if inf > biggest_inf || num_pvals < lowest_num_pvals
                        lowest_num_pvals = num_pvals
                        biggest_inf = inf
                        best_vidx = vidx
                        is_solution = false
                    end
                end
            end
        end
    end
    return BranchVarObj(true, is_solution, best_vidx)
end

function get_next_branch_variable(com::CS.CoM, ::Val{:Random})
    variables = com.search_space
    is_free = com.activity_vars.is_free
    if com.c_backtrack_idx == 1
        com.activity_vars.is_free = zeros(Bool, length(variables))
        is_free = com.activity_vars.is_free
    else
        is_free .= false
    end

    for variable in variables
        if !isfixed(variable)
            is_free[variable.idx] = true
        end
    end
    is_solution = !any(is_free)
    vidx = -1
    if !is_solution
        vidx = sample(CS_RNG, 1:length(variables), Weights(is_free))
    end
    return BranchVarObj(true, is_solution, vidx)
end

"""
    still_probing(n, μ, variance, max_deviation)

Determine whether probing should be continued.
This depends on
- the number of probes (at least 2)
- the variance and standard deviation
- whether the 95% confidence interval is in µ ± max_deviation * µ
- max_deviation should be options.activity.max_confidence_deviation / 100
"""
function still_probing(n, μ, variance, max_deviation)
    n < 2 && return true
    for i in 1:length(μ)
        σ = sqrt(variance[i])
        if tdistcdf(n - 1, 0.05) * (σ / sqrt(n)) > max_deviation * μ[i]
            range = tdistcdf(n - 1, 0.05) * (σ / sqrt(n))
            return true
        end
    end
    return false
end

"""
    probe_until(com::CS.CoM)

Start probing from current node:
- Create paths in a depth first search way and use a random branch variable strategy
- Update activity
This creates new `BacktrackObj` inside `backtrack_vec` which get overwritten in each probe
"""
function probe_until(com::CS.CoM)
    mean_activities = [var.activity for var in com.search_space]
    variance_activities = zeros(length(com.search_space))
    saved_traverse_strategy = com.traverse_strategy
    saved_traverse_strategy_option = com.options.traverse_strategy
    saved_branch_strategy = com.branch_strategy
    saved_branch_split = com.branch_split

    # options for probing
    com.options.traverse_strategy = :DFS
    com.traverse_strategy = Val(:DFS)
    com.branch_split = Val(:Random)
    com.branch_strategy = Val(:Random)
    global_feasible = true
    changed_traverse_strategy!(com, saved_traverse_strategy)

    backtrack_obj = BacktrackObj(com)
    backtrack_obj.idx = length(com.backtrack_vec) + 1
    backtrack_obj.parent_idx = 1
    backtrack_obj.depth = 1
    com.c_step_nr += 1
    set_current_step_nr!(backtrack_obj, com)

    addBacktrackObj2Backtrack_vec!(com.backtrack_vec, backtrack_obj, com)
    update_log_node!(com, 2)
    # the option itself is in percentage
    max_deviation = com.options.activity.max_confidence_deviation / 100

    n = 1
    while n <= com.options.activity.max_probes &&
              still_probing(n, mean_activities, variance_activities, max_deviation) &&
              global_feasible
        n += 1
        root_feasible, feasible, activities, last_backtrack_id = probe(com)
        for i in 1:length(com.search_space)
            new_mean = mean_activities[i] + (activities[i] - mean_activities[i]) / n
            # update variance: https://math.stackexchange.com/questions/102978/incremental-computation-of-standard-deviation
            if n != 1
                variance_activities[i] =
                    (
                        (n - 2) * variance_activities[i] +
                        (n - 1) * (mean_activities[i] - new_mean)^2 +
                        (activities[i] - new_mean)^2
                    ) / (n - 1)
            end
            mean_activities[i] = new_mean
        end
        # if root infeasible we can remove that setting completely
        if !root_feasible
            backtrack_obj = com.backtrack_vec[last_backtrack_id]
            vidx = backtrack_obj.vidx
            vidx == 0 && break # all variables are fixed already
            # use the variable from com not ccom to remove it from the actual model
            variable = com.search_space[vidx]
            lb = backtrack_obj.lb
            ub = backtrack_obj.ub
            push!(com.root_infeasible_vars, VarAndVal(vidx, lb, ub))
        end
    end
    for i in 1:length(com.search_space)
        com.search_space[i].activity = mean_activities[i]
    end
    com.branch_split = saved_branch_split
    com.traverse_strategy = saved_traverse_strategy
    com.options.traverse_strategy = saved_traverse_strategy_option
    com.branch_strategy = saved_branch_strategy
    changed_traverse_strategy!(com, Val(:DFS))

    # checkout root node
    checkout_from_to!(com, com.c_backtrack_idx, 1)
    # prune the last step as checkout_from_to! excludes the to part
    restore_prune!(com, 1)

    return global_feasible
end

"""
    probe(com::CS.CoM)

Probe from node root node
- Follow a DFS path by random selection of a branching variable
- Update activity
Return if feasible and the created backtrack ids along the way
"""
function probe(com::CS.CoM)
    activities = zeros(length(com.search_space))

    first_parent_idx = 2
    options = com.options

    status, last_backtrack_id = backtrack!(
        com,
        options.max_bt_steps;
        sorting = options.backtrack_sorting,
        log_table = false,
        first_parent_idx = first_parent_idx,
        single_path = true,
        compute_bounds = false,
        check_bounds = false,
        cb_finished_pruning = (com1) -> update_probe_activity!(activities, com1),
    )

    feasible = com.backtrack_vec[last_backtrack_id].is_feasible
    root_feasible = true
    if last_backtrack_id == first_parent_idx ||
       com.backtrack_vec[last_backtrack_id].parent_idx == first_parent_idx
        root_feasible = feasible
    end

    # checkout root node
    checkout_from_to!(com, last_backtrack_id, first_parent_idx)
    # prune the last step as checkout_from_to! excludes the to part
    restore_prune!(com, first_parent_idx)
    return root_feasible, feasible, activities, last_backtrack_id
end

function update_activity!(com)
    c_backtrack_idx = com.c_backtrack_idx
    c_step_nr = com.c_step_nr
    backtrack_obj = com.backtrack_vec[c_backtrack_idx]
    branch_vidx = backtrack_obj.vidx
    γ = com.options.activity.decay
    for variable in com.search_space
        variable.idx == branch_vidx && continue
        if nvalues(variable) > 1
            variable.activity *= γ
        end
        if num_changes(variable, c_step_nr) > 0
            variable.activity += 1
        end
    end
end

function update_probe_activity!(activities, com)
    c_backtrack_idx = com.c_backtrack_idx
    c_step_nr = com.c_step_nr
    backtrack_obj = com.backtrack_vec[c_backtrack_idx]
    branch_vidx = backtrack_obj.vidx
    for variable in com.search_space
        variable.idx == branch_vidx && continue
        if num_changes(variable, c_step_nr) > 0
            activities[variable.idx] += 1
        end
    end
end

function get_next_branch_variable(com::CS.CoM, ::Val{:ABS})
    update_activity!(com)

    if com.in_probing_phase
        com.in_probing_phase = false
        global_feasible = probe_until(com)
        if !global_feasible
            return BranchVarObj(false, false, -1)
        end
    end

    # activity based search
    is_in_objective = false
    highest_activity = -1.0
    best_vidx = -1
    is_solution = true

    for variable in com.search_space
        if !isfixed(variable)
            vidx = variable.idx
            num_pvals = nvalues(variable)
            activity_ratio = variable.activity / num_pvals
            if !is_in_objective && com.var_in_obj[vidx]
                highest_activity = activity_ratio
                is_in_objective = true
                is_solution = false
                best_vidx = vidx
                continue
            end
            if !is_in_objective || com.var_in_obj[vidx]
                if activity_ratio >= highest_activity
                    highest_activity = activity_ratio
                    best_vidx = vidx
                    is_solution = false
                end
            end
        end
    end
    return BranchVarObj(true, is_solution, best_vidx)
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Flux, StochasticDiffEq, LinearAlgebra, Statistics, Test, NeuralNetDiffEq
using Random
Random.seed!(0)

d = 100 # number of dimensions
x0 = fill(0.0f0,d)
tspan = (0.0f0, 1.0f0)
m = 100 # number of trajectories (batch size)
λ = 1.0f0
#
g(X) = log(0.5f0 + 0.5f0*sum(X.^2))
f(X,u,σᵀ∇u,p,t) = -λ*sum(σᵀ∇u.^2)
μ(X,p,t) = zero(X)  #Vector d x 1 λ
σ(X,p,t) = Diagonal(sqrt(2.0f0)*ones(Float32,d)) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ, σ, x0, tspan)

hls = 10 + d #hidden layer size
opt = Flux.ADAM(0.03)  #optimizer
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
#
@time ans = solve(prob, pdealg, verbose=true, maxiters=500, trajectories=m,
                            alg=LambaEM(), pabstol = 1f-2, reltol=1e-4, abstol=1e-4)

T = tspan[2]
MC = 10^5
W() = randn(d,1)
u_analytical(x, t) = -(1/λ)*log(mean(exp(-λ*g(x .+ sqrt(2.0)*abs.(T-t).*W())) for _ = 1:MC))
analytical_ans = u_analytical(x0, tspan[1])

error_l2 = sqrt((ans - analytical_ans)^2/ans^2)

println("Hamilton Jacobi Bellman Equation")
# println("numerical = ", ans)
# println("analytical = " , analytical_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 0.2

ext: jl


text:
using Test
using LogisticRegression: propagate!

w = [1.0 2.0]
b = [2.0]

x1 = [1.0, 3.0]
x2 = [2.0, 4.0]
x3 = [-1.0, -3.2]

y1 = 1.0
y2 = 0.0
y3 = 1.0

J = zeros(1)
dw = zeros(1, 2)
db = zeros(1)
z = zeros(1)
a = zeros(1)
dz = zeros(1)
propagate!(J, dw, db, w, b, [x1, x2, x3], [y1, y2, y3], z, a, dz)

@test isapprox(first(J), 5.801545319394553)
@test isapprox(dw, [0.99845601 2.39507239])
@test isapprox(db, [0.00145557813678])


text:
using Dates
using Test
using CliOptions

include("testutils.jl")

@testset "Option()" begin
    @testset "ctor" begin
        @test_throws MethodError Option()
        @test_throws MethodError Option("-f", "--foo", "--bar")
        @test_throws ArgumentError Option("")
        @test_throws ArgumentError Option("a")
        @test_throws ArgumentError Option("-")
        @test_throws ArgumentError Option("--")
        @test Option("-a").names == ("-a",)
        @test Option("-a", "-b").names == ("-a", "-b")

        @test Option(String, "-a").T == String
        @test Option(DateTime, "-a").T == DateTime
        @test Option(UInt32, "-a").T == UInt32
    end

    @testset "ctor; duplicates, $(v[1])" for v in [
        (["-f", "--foo"], (nothing, nothing)),
        (["-f", "-f"], (ArgumentError, "-f")),
    ]
        names, expected = v
        if expected[1] isa Type
            tr = @test_throws expected[1] Option(names...)
            if tr isa Test.Pass
                buf = IOBuffer()
                showerror(buf, tr.value)
                msg = String(take!(buf))
                @test msg == ("ArgumentError: Duplicate names for an Option found: " *
                              expected[2])
            end
        else
            @test Option(names...) !== nothing
        end
    end

    @testset "show(x); $(join(v[1],','))" for v in [
        (["-a"], "Option(:a)"),
        (["-a", "--foo-bar"], "Option(:a,:foo_bar)"),
    ]
        names, expected_repr = v
        option = Option(names...)
        @test repr(option) == expected_repr
    end

    @testset "show(io, x)" begin
        let option = Option("-a")
            buf = IOBuffer()
            redirect_stdout(buf) do
                show(option)
            end
            @test String(take!(buf)) == "Option(:a)"
        end
    end

    @testset "consume!(); $(v[1])" for v in [
        ([""], (0, nothing)),
        (["-a"], (0, nothing)),
        (["-d"], CliOptionError),
        (["-d", "3"], (2, "3")),
    ]
        args, expected = v
        option = Option("-d", "--depth")
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        if expected isa Type
            @test_throws expected CliOptions.consume!(d, option, args, ctx)
        else
            num_consumed = CliOptions.consume!(d, option, args, ctx)
            @test num_consumed == expected[1]
            if expected[2] !== nothing
                @test d["d"] == expected[2]
                @test d["depth"] == expected[2]
            end
        end
    end

    @testset "consume!(); type, $(v[1])" for v in [
        # ctor
        ("constructible", Date, ["-a", "2006-01-02"], (2, Date(2006, 1, 2))),

        # Base.parse
        ("parsable, -1", UInt8, ["-a", "-1"], CliOptionError),
        ("parsable, 0", UInt8, ["-a", "0"], (2, UInt8(0))),
        ("parsable, 255", UInt8, ["-a", "255"], (2, UInt8(255))),
        ("parsable, 256", UInt8, ["-a", "256"], CliOptionError),

        # N/A
        ("unacceptable", AbstractFloat, ["-a", "42"], CliOptionError),
    ]
        _, T, args, expected = v
        option = Option(T, "-a")
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        if expected isa Type
            @test_throws expected CliOptions.consume!(d, option, args, ctx)
        else
            num_consumed = CliOptions.consume!(d, option, args, ctx)
            @test num_consumed == expected[1]
            @test d["a"] == expected[2]
        end
    end

    @testset "consume!(); requirement, $(v[1])" for v in [
        ("[foo, bar], foo", ["-a", "foo"],
            String, ["foo", "bar"], (2, "foo")),
        ("[foo, bar], bar", ["-a", "bar"],
            String, ["foo", "bar"], (2, "bar")),
        ("[foo, bar], qux", ["-a", "qux"],
            String, ["foo", "bar"], (CliOptionError, "must be one of")),
        ("(foo, bar), foo", ["-a", "foo"],
            String, ("foo", "bar"), (2, "foo")),
        ("(foo, bar), bar", ["-a", "bar"],
            String, ("foo", "bar"), (2, "bar")),
        ("(foo, bar), qux", ["-a", "qux"],
            String, ("foo", "bar"), (CliOptionError, "must be one of")),
        ("/qu+x/, quux/", ["-a", "quux"],
            String, Regex("qu+x"), (2, "quux")),
        ("/qu+x/, qux", ["-a", "qux"],
            String, Regex("qu+x"), (2, "qux")),
        ("/qu+x/, qx", ["-a", "qx"],
            String, Regex("qu+x"), (CliOptionError, "must match for")),
        ("String -> Bool, foo", ["-a", "foo"],
            String, s -> startswith(s, "foo"), (2, "foo")),
        ("String -> Bool, 6", ["-a", "6"],
            Int, n -> iseven(n), (2, 6)),
        ("String -> Bool, 7", ["-a", "7"],
            Int, n -> iseven(n), (CliOptionError, "validation failed")),
        ("String -> String, foo", ["-a", "foo"],
            String, s -> startswith(s, "foo") ? "" : "It's not foo", (2, "foo")),
        ("String -> String, 6", ["-a", "6"],
            Int, n -> iseven(n) ? "" : "must be even", (2, 6)),
        ("String -> String, 7", ["-a", "7"],
            Int, n -> iseven(n) ? "" : "must be even", (CliOptionError, "must be even")),
    ]
        _, args, T, requirement, expected = v
        option = Option(T, "-a"; requirement = requirement)
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        if expected[1] isa Type
            try
                CliOptions.consume!(d, option, args, ctx)
                @test false  # Exception must be thrown
            catch ex
                @test ex isa expected[1]
                @test occursin(args[1], ex.msg)
                @test occursin(expected[2], ex.msg)
            end
        else
            num_consumed = CliOptions.consume!(d, option, args, ctx)
            @test num_consumed == expected[1]
            @test d["a"] == expected[2]
        end
    end

    @testset "consume!(); until, $(v[1])" for v in [
        ("nothing, String", String, nothing, (2, "1"))
        ("nothing, Int", Int, nothing, (2, 1))
        ("4, String", String, "4", (5, ["1", "2", "3"]))
        ("4, Int", Int, "4", (5, [1, 2, 3]))
        ("[4, 3], String", String, ["4", "3"], (4, ["1", "2"]))
        ("[4, 3], Int", Int, ["4", "3"], (4, [1, 2]))
        ("(4, 3), String", String, ("4", "3"), (4, ["1", "2"]))
        ("(4, 3), Int", Int, ("4", "3"), (4, [1, 2]))
    ]
        _, T, until, expected = v
        args = split("-a 1 2 3 4")
        option = Option(T, "-a"; until = until)
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        if expected[1] isa Type
            try
                CliOptions.consume!(d, option, args, ctx)
                @test false  # Exception must be thrown
            catch ex
                @test ex isa expected[1]
                @test occursin(args[1], ex.msg)
                @test occursin(expected[2], ex.msg)
            end
        else
            num_consumed = CliOptions.consume!(d, option, args, ctx)
            @test num_consumed == expected[1]
            @test d["a"] == expected[2]
        end
    end

    @testset "check_usage_count(); $(v[1])" for v in [
        ("required, 0", missing, 0, CliOptionError),
        ("required, 1", missing, 1, nothing),
        ("omittable, 0", "foo", 0, nothing),
    ]
        _, default, count, expected = v
        option = Option("-n", default = default)
        ctx = CliOptions.ParseContext()
        ctx.usage_count[option] = count
        if expected isa Type
            @test_throws expected CliOptions.check_usage_count(option, ctx)
        else
            @test true  # No exception was thrown
        end
    end
end


text:
# utilities for pipe diameters

"compute equivalent diameter from serial pipe segments"
function serialmerge(diams::Vector{Diameter}, lengths::Vector{Float64})
    @assert length(diams) == length(lengths)
    @assert sum(lengths) ≈ 1.0
    y = sum(l / d.value^5 for (d,l) in zip(diams, lengths))
    y^(-1/5)
end

"compute lengths of two segments for equivalent diameter"
function pipesplit(diams::Vector, equiv::Float64)
    ds = [d.value for d in diams]
    dmin, dmax = extrema(ds)
    @assert isapproxin(equiv, dmin, dmax, atol=1e-8)
    y, ys = equiv^(-5), ds.^(-5)
    @assert issorted(ys, rev=true) # decreasing

    l = fill(0.0, length(diams))

    right = findfirst(yi -> y ≥ yi, ys) # decreasing
    if right == 1 # on left boundary
        l[1] = 1.0
        return l
    end

    left = right - 1
    λ = (y - ys[left])/(ys[right] - ys[left])
    l[left] = 1 - λ
    l[right] = λ
    l
end


text:
# Overall structure
#----------------------------

export parse_commandline_counterfactuals, print_parsed_args_counterfactuals

function arg_settings_counterfactuals()
    s = ArgParseSettings()

    @add_arg_table! s begin

        "--jld2"
            arg_type = String
            default = "E:/projects/ShaleDrillingResults/2019-12/16481962/estimation-results-16481962-no-rigs-data_all_leases.jld2"

        "--noPar"
            action = :store_true

        "--dateStop"
            arg_type = Float64
            default = 2016.75

        "--techYearZero"
            arg_type = Int
            default = 2007

        "--rFileDir"
            arg_type = String
            default = "./"
    end
    return s
end

parse_commandline_counterfactuals() = parse_args(arg_settings_counterfactuals())

function print_parsed_args_counterfactuals(x::Dict)

    kys = keys(x)
    padlen = maximum(length.(kys))
    padk(k) = rpad(k,padlen)

    println("--------------------------------------")
    println("ARGUMENTS PASSED")
    println("--------------------------------------")
    for (k,v) in x
        println("  $(padk(k))  =>  $v")
    end
    println("--------------------------------------")
    print("\n\n")
    return nothing
end


text:
rng = StableRNG(666899)

@testset "built-in regressor measures" begin
    y    = [1, 2, 3, 4]
    yhat = [4, 3, 2, 1]
    w = [1, 2, 4, 3]
    @test isapprox(mae(yhat, y), 2)
    @test isapprox(mae(yhat, y, w), (1*3 + 2*1 + 4*1 + 3*3)/4)
    @test isapprox(rms(yhat, y), sqrt(5))
    @test isapprox(rms(yhat, y, w), sqrt((1*3^2 + 2*1^2 + 4*1^2 + 3*3^2)/4))
    @test isapprox(mean(l1(yhat, y)), 2)
    @test isapprox(mean(l1(yhat, y, w)), mae(yhat, y, w))
    @test isapprox(mean(l2(yhat, y)), 5)
    @test isapprox(mean(l2(yhat, y, w)), rms(yhat, y, w)^2)
    @test isapprox(mean(log_cosh(yhat, y)), 1.3715546675)

    yhat = y .+ 1
    @test isapprox(rmsl(yhat, y),
                   sqrt((log(1/2)^2 + log(2/3)^2 + log(3/4)^2 + log(4/5)^2)/4))
    @test isapprox(rmslp1(yhat, y),
                   sqrt((log(2/3)^2 + log(3/4)^2 + log(4/5)^2 + log(5/6)^2)/4))
    @test isapprox(rmsp(yhat, y), sqrt((1 + 1/4 + 1/9 + 1/16)/4))
    @test isapprox(mape(yhat, y), (1/1 + 1/2 + 1/3 + 1/4)/4)
end

@testset "MLJBase.value" begin
    yhat = randn(rng,5)
    X = (weight=randn(rng,5), x1 = randn(rng,5))
    y = randn(rng,5)
    w = randn(rng,5)

    @test MLJBase.value(mae, yhat, nothing, y, nothing) ≈ mae(yhat, y)
    @test MLJBase.value(mae, yhat, nothing, y, w) ≈ mae(yhat, y, w)

    spooky(yhat, y) = abs.(yhat - y) |> mean
    @test MLJBase.value(spooky, yhat, nothing, y, nothing) ≈ mae(yhat, y)

    cool(yhat, y, w) = abs.(yhat - y) .* w |> mean
    MLJBase.supports_weights(::Type{typeof(cool)}) = true
    @test MLJBase.value(cool, yhat, nothing, y, w) ≈ mae(yhat, y, w)

    funky(yhat, X, y) = X.weight .* abs.(yhat - y) |> mean
    MLJBase.is_feature_dependent(::Type{typeof(funky)}) = true
    @test MLJBase.value(funky, yhat, X, y, nothing) ≈ mae(yhat, y, X.weight)

    weird(yhat, X, y, w) = w .* X.weight .* abs.(yhat - y) |> mean
    MLJBase.is_feature_dependent(::Type{typeof(weird)}) = true
    MLJBase.supports_weights(::Type{typeof(weird)}) = true
    @test MLJBase.value(weird, yhat, X, y, w) ≈ mae(yhat, y, X.weight .* w)
end


text:
# Observations:
#
# By and large, we really smoke AMOS in speed. Especially for large arg.
# Assembling a reasonably large matrix (1024 x 1024) gets done in about half the
# time. Which may seem like nothing, but a factor of two never hurts...

using BenchmarkTools, SpecialFunctions, StaticArrays
include("../besk.jl")
include("shared.jl")

const NUs = (0.1, 0.25, 0.5, 0.75, 0.999, 1.0, 1.001, 1.25, 1.5, 1.75, 2.0,
             2.05, 2.75, 3.1, 3.8, 4.3, 4.9)

const TINY_X   = 1e-4
const SMALL_X  = 0.25
const MID_X    = 7.5
const LARGE_X  = 20.0
const Xs = (1e-4, 0.25, 5.0, 7.5, 8.5, 14.0, 17.0, 29.0, 31.0, 50.0)

const PTS  = [rand(3).*10.0 for _ in 1:1024]
const PRMS = @SVector [1.0, 1.0, 1.25]

print("\n\n")
println("##################")
println("HEAT ONE: pointwise timings.")
println("##################")
print("\n\n")
for (v, x) in Iterators.product(NUs, Xs)
  println("(v,x) = ($v, $x):")
  print("Timing for AMOS:")
  @btime SpecialFunctions.besselk($v, $x)
  print("Timing for us:")
  try
    @btime _besselk($v, $x)
  catch
    println("FAILURE/ERROR OUT.")
  end
  print("\n\n")
end

print("\n\n")
println("##################")
println("HEAT TWO: kernel matrix assembly.")
println("##################")
print("\n\n")
println("Timings for AMOS:")
@btime assemble_matrix(matern_amos, $PTS, $PRMS)
println("Timings for US:")
@btime assemble_matrix(matern_us, $PTS, $PRMS)


text:
# ----------------------------------------------------------
# ---------------  Apply operator to arrays  ---------------
# ----------------------------------------------------------

# In case of applying a single FunctionOperator to the data, everything is quite straightforward:
# We make sure that operator and the data (and maybe the buffer) works together, and then apply
# the forw or backw field of FunctionOperator based on the value of adjoint field.

Base.:*(FO::FunctionOperator, A::AbstractArray) = begin
    assertType(FO, A)
    FunctionOperators_global_settings.auto_reshape ? (A = reshape(A, FO.inDims)) : assertMultDim(FO, A)
    info("Allocation of buffer0, size: $(FO.outDims)")
    result = FO.adjoint ?
        FO.ismutating ? FO.backw(Array{eltype(FO)}(undef, FO.outDims), A) : FO.backw(A) :
        FO.ismutating ? FO.forw(Array{eltype(FO)}(undef, FO.outDims), A) : FO.forw(A)
    FunctionOperators_global_settings.auto_reshape ? reshape(result, FO.outDims) : result
end

LinearAlgebra.mul!(buffer::AbstractArray, FO::FunctionOperator, A::AbstractArray) = begin
    assertType(FO, A)
    assertTypeBuffer(FO, buffer)
    orig_buffer = buffer
    if FunctionOperators_global_settings.auto_reshape
        A = reshape(A, FO.inDims)
        buffer = reshape(buffer, FO.outDims)
    else
        assertMultDim(FO, A)
        assertMultDimBuffer(FO, buffer)
    end
    result = FO.adjoint ?
        (checkTwoInputs(FO.backw) ? FO.backw(buffer, A) : buffer .= FO.backw(A)) :
        (checkTwoInputs(FO.forw) ? FO.forw(buffer, A) : buffer .= FO.forw(A))
    (buffer !== result) && (buffer .= result)
    orig_buffer # buffer might have been reshaped, so we return the un-reshaped version
end

# When FunctionOperatorComposite is applied to data, the process changes only slightly:
# After checking if plan is already calculated (if not, then we allocate the buffer, and then calculate the plan), then we can apply the aggrataged function to the data.

Base.:*(FO::FunctionOperatorComposite, A::AbstractArray) = begin
    assertType(FO, A)
    FunctionOperators_global_settings.auto_reshape ? (A = reshape(A, FO.inDims)) : assertMultDim(FO, A)
    storage = Vector{Buffer}(undef, 0)
    buffer0 = Buffer(Array{eltype(FO)}(undef, FO.outDims), "buffer0", 1, true)
    push!(storage, buffer0)
    if FO.plan_function == noplan
        FO.plan_function, output, FO.plan_string = getPlan(FO, buffer0, FO.adjoint, "x", storage)
        info("Plan calculated: $(output.name) .= "*FO.plan_string)
        @assert output == buffer0 "Implementation error: Output of computation is written to $(output.name) instead of buffer1"
    end
    result = FO.plan_function(buffer0.buffer, A)
    FunctionOperators_global_settings.auto_reshape ? reshape(result, FO.outDims) : result
end

LinearAlgebra.mul!(buffer::AbstractArray, FO::FunctionOperatorComposite, A::AbstractArray) = begin
    assertType(FO, A)
    assertType(FO, buffer)
    orig_buffer = buffer
    if FunctionOperators_global_settings.auto_reshape
        A = reshape(A, FO.inDims)
        buffer = reshape(buffer, FO.outDims)
    else
        assertMultDim(FO, A)
        assertMultDimBuffer(FO, buffer)
    end
    if FO.plan_function == noplan
        storage = Array{Buffer}(undef, 0)
        buffer0 = Buffer(buffer, "buffer0", 1, true)
        push!(storage, buffer0)
        info("buffer0 = <previously allocated>")
        FO.plan_function, output, FO.plan_string = getPlan(FO, buffer0, FO.adjoint, "x", storage)
        info(("Plan calculated: $(output.name) .= "*FO.plan_string))
        @assert output.buffer === buffer "Implementation error: Output of computation is written to $(output.name) instead of buffer1"
    end
    result = FO.plan_function(buffer, A)
    (buffer !== result) && (buffer .= result)
    orig_buffer # buffer might have been reshaped, so we return the un-reshaped version
end


text:
import Base: getindex, iterate

"""
    damas(env::Environment, b [,f<:AbstractArray, ω::Real]; maxiter=10)
Performs exactly `maxiter` DAMAS iterations for all frequency bins in `env.fn` or frequencies in `f`. `f` must contain values that match exact a (sub)set of the values in `env.fn`. Successive Over Relaxation (SOR) can be set with relaxation parameter `ω`. Default is `ω=1` corresponding to no relaxation.
"""
function damas(env::Environment, b::FA, f::T=env.fn, ω::Real=1.0; maxiter::Int=10) where {T <: AbstractArray, FA <: FreqArray}
    env.fn == f ? (f_inds=1:env.Nf) : (f_inds = findall(x->x in f, env.fn))
    x = zeros(size(b,1),length(f_inds))
    for (i,fn) in enumerate(f_inds)
        @views damas!(x[:,i], env.steeringvec.arr[:,:,fn], b[:,fn], ω; maxiter=maxiter)
    end
    return FreqArray(x,env.fn[f_inds])
end

function damas!(x, steer, b, ω::Real=1.0; maxiter::Int=10)
    iterable = DAMASSORIterable(x, steer, similar(x), similar(x), b, ω, maxiter)
    for _ = iterable end
    x
end

mutable struct DAMASSORIterable{solT,steerT,vecT,rhsT,numT}
    x::solT
    steer::steerT
    tmp1::vecT
    tmp2::vecT
    b::rhsT
    ω::numT
    maxiter::Int
end

start(::DAMASSORIterable) = 1
done(it::DAMASSORIterable, iteration::Int) = iteration > it.maxiter
function iterate(s::DAMASSORIterable, iteration::Int=start(s))
    if done(s, iteration) return nothing end

    n = size(s.b, 1)

    for col = 1 : n
        AeroAcoustics.psf_col!(s.tmp2,s.steer,col)
        @simd for row = 1 : col - 1
            @inbounds s.tmp1[row] = max(0.0,s.tmp1[row] - s.tmp2[row] * s.x[col])
        end

        s.tmp1[col] = s.b[col]
    end

    for col = 1 : n
        AeroAcoustics.psf_col!(s.tmp2,s.steer,col)
        @inbounds s.x[col] += s.ω * (s.tmp1[col] / s.tmp2[col] - s.x[col])
        @simd for row = col + 1 : n
            @inbounds s.tmp1[row] = max(0.0,s.tmp1[row] - s.tmp2[row] * s.x[col])
        end
    end

    nothing, iteration + 1
end


text:
type RollingWindow{T} <: AbstractVector{T}
  cb::CircularBuffer{T}
  lags::IntIterable
end

function RollingWindow{T}(::Type{T}, lags::IntIterable)
  minlag = minimum(lags)
  @assert minlag >= 0
  maxlag = maximum(lags)
  @assert maxlag >= minlag
  RollingWindow{T}(CircularBuffer(T, maxlag+1), lags)
end

function Base.print(io::IO, window::RollingWindow)
  println(io, "RollingWindow{n=", length(window), (isfull(window) ? "" : " cap=$(capacity(window))"), " cb:")
  for item in window
    println(io, "  ", item)
  end
  print(io, "}")
end
Base.show(io::IO, window::RollingWindow) = print(io, window)

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

function bufferIndex(window::RollingWindow, i::Int)
  if i < 1 || i > length(window)
    error("RollingWindow out of range. window=$window i=$i")
  end
  length(window.cb) - window.lags[i]
end

Base.getindex(window::RollingWindow, i::Int) = window.cb[bufferIndex(window, i)]
# Base.getindex{T}(window::RollingWindow{T}, itr::IntIterable) = T[window[i] for i in itr]
function Base.setindex!{T}(window::RollingWindow, data::T, i::Int)
  window.cb[bufferIndex(window, i)] = data
  nothing
end

# # iterate with the syntax "for x in list"
# Base.start{T}(window::RollingWindow{T}) = 1
# Base.done{T}(window::RollingWindow{T}, state::Int) = state > length(window)
# Base.next{T}(window::RollingWindow{T}, state::Int) = (window.cb[bufferIndex(window, state)], state+1)

# note: length is 0 until the underlying buffer is full
Base.length(window::RollingWindow) = (isfull(window) ? capacity(window) : 0)
Base.size(cb::CircularBuffer) = (length(cb),)
# Base.isempty(window::RollingWindow) = length(window) == 0
# Base.endof(window::RollingWindow) = length(window)
# Base.first(window::RollingWindow) = (isempty(window) ? error("Rolling window is empty. can't take first.") : window[1])
# Base.last(window::RollingWindow) = (isempty(window) ? error("Rolling window is empty. can't take last.") : window[length(window)])

capacity(window::RollingWindow) = length(window.lags)
isfull(window::RollingWindow) = isfull(window.cb)
toarray{T}(window::RollingWindow{T}) = window[1:length(window)]

Base.push!{T}(window::RollingWindow{T}, data::T) = push!(window.cb, data)
Base.append!{T}(window::RollingWindow{T}, datavec::AbstractVector{T}) = append!(window.cb, datavec)


text:
"Distance between two values"
function d end

@inline d(x::Real, y::Real) = (xy = (x - y); xy * xy)
# @inline d(x::Vector{<:Real}, y::Vector{<:Real}) = norm(x - y)
@inline d(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = sum(d.(x,y))
@inline d(x::NTuple{N, <: Real}, y::NTuple{N, <:Real}) where N = sum(d.(x,y))
@inline d(x::AbstractArray{<:Real}, y::AbstractArray{<:Real}) = norm(x[:] - y[:])

"Distance from x to [a, b]"
function bound_loss(x, a, b)
  # @pre b >= a
  if x < a
    a - x
  elseif x > b
    x - b
  else
    zero(x)
  end
end


text:
# Structured matrices

#####
##### `Diagonal`
#####

function rrule(::Type{<:Diagonal}, d::AbstractVector)
    function Diagonal_pullback(ȳ::AbstractMatrix)
        return (NO_FIELDS, diag(ȳ))
    end
    function Diagonal_pullback(ȳ::Composite)
        # TODO: Assert about the primal type in the Composite, It should be Diagonal
        # infact it should be exactly the type of `Diagonal(d)`
        # but right now Zygote loses primal type information so we can't use it.
        # See https://github.com/FluxML/Zygote.jl/issues/603
        return (NO_FIELDS, ȳ.diag)
    end
    return Diagonal(d), Diagonal_pullback
end

function rrule(::typeof(diag), A::AbstractMatrix)
    function diag_pullback(ȳ)
        return (NO_FIELDS, @thunk(Diagonal(ȳ)))
    end
    return diag(A), diag_pullback
end
if VERSION ≥ v"1.3"
    function rrule(::typeof(diag), A::AbstractMatrix, k::Integer)
        function diag_pullback(ȳ)
            return (NO_FIELDS, @thunk(diagm(size(A)..., k => ȳ)), DoesNotExist())
        end
        return diag(A, k), diag_pullback
    end

    function rrule(::typeof(diagm), m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)
        function diagm_pullback(ȳ)
            return (NO_FIELDS, DoesNotExist(), DoesNotExist(), _diagm_back.(kv, Ref(ȳ))...)
        end
        return diagm(m, n, kv...), diagm_pullback
    end
end
function rrule(::typeof(diagm), kv::Pair{<:Integer,<:AbstractVector}...)
    function diagm_pullback(ȳ)
        return (NO_FIELDS, _diagm_back.(kv, Ref(ȳ))...)
    end
    return diagm(kv...), diagm_pullback
end

function _diagm_back(p, ȳ)
    return Thunk() do
        k, v = p
        d = diag(ȳ, k)[1:length(v)] # handle if diagonal was smaller than matrix
        return Composite{typeof(p)}(second = d)
    end
end

function rrule(::typeof(*), D::Diagonal{<:Real}, V::AbstractVector{<:Real})
    function times_pullback(Ȳ)
        return (NO_FIELDS, @thunk(Diagonal(Ȳ .* V)), @thunk(D * Ȳ))
    end
    return D * V, times_pullback
end

#####
##### `Symmetric`
#####

function rrule(::Type{<:Symmetric}, A::AbstractMatrix)
    function Symmetric_pullback(ȳ)
        return (NO_FIELDS, @thunk(_symmetric_back(ȳ)))
    end
    return Symmetric(A), Symmetric_pullback
end

_symmetric_back(ΔΩ) = UpperTriangular(ΔΩ) + LowerTriangular(ΔΩ)' - Diagonal(ΔΩ)
_symmetric_back(ΔΩ::Union{Diagonal,UpperTriangular}) = ΔΩ

#####
##### `Adjoint`
#####

# ✖✖✖TODO: Deal with complex-valued arrays as well
function rrule(::Type{<:Adjoint}, A::AbstractMatrix{<:Real})
    function Adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(adjoint(ȳ)))
    end
    return Adjoint(A), Adjoint_pullback
end

function rrule(::Type{<:Adjoint}, A::AbstractVector{<:Real})
    function Adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(vec(adjoint(ȳ))))
    end
    return Adjoint(A), Adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractMatrix{<:Real})
    function adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(adjoint(ȳ)))
    end
    return adjoint(A), adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractVector{<:Real})
    function adjoint_pullback(ȳ)
        return (NO_FIELDS, @thunk(vec(adjoint(ȳ))))
    end
    return adjoint(A), adjoint_pullback
end

#####
##### `Transpose`
#####

function rrule(::Type{<:Transpose}, A::AbstractMatrix)
    function Transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk transpose(ȳ))
    end
    return Transpose(A), Transpose_pullback
end

function rrule(::Type{<:Transpose}, A::AbstractVector)
    function Transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk vec(transpose(ȳ)))
    end
    return Transpose(A), Transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractMatrix)
    function transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk transpose(ȳ))
    end
    return transpose(A), transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractVector)
    function transpose_pullback(ȳ)
        return (NO_FIELDS, @thunk vec(transpose(ȳ)))
    end
    return transpose(A), transpose_pullback
end

#####
##### Triangular matrices
#####

function rrule(::Type{<:UpperTriangular}, A::AbstractMatrix)
    function UpperTriangular_pullback(ȳ)
        return (NO_FIELDS, @thunk Matrix(ȳ))
    end
    return UpperTriangular(A), UpperTriangular_pullback
end

function rrule(::Type{<:LowerTriangular}, A::AbstractMatrix)
    function LowerTriangular_pullback(ȳ)
        return (NO_FIELDS, @thunk Matrix(ȳ))
    end
    return LowerTriangular(A), LowerTriangular_pullback
end

function rrule(::typeof(triu), A::AbstractMatrix, k::Integer)
    function triu_pullback(ȳ)
        return (NO_FIELDS, @thunk(triu(ȳ, k)), DoesNotExist())
    end
    return triu(A, k), triu_pullback
end
function rrule(::typeof(triu), A::AbstractMatrix)
    function triu_pullback(ȳ)
        return (NO_FIELDS, @thunk triu(ȳ))
    end
    return triu(A), triu_pullback
end

function rrule(::typeof(tril), A::AbstractMatrix, k::Integer)
    function tril_pullback(ȳ)
        return (NO_FIELDS, @thunk(tril(ȳ, k)), DoesNotExist())
    end
    return tril(A, k), tril_pullback
end
function rrule(::typeof(tril), A::AbstractMatrix)
    function tril_pullback(ȳ)
        return (NO_FIELDS, @thunk tril(ȳ))
    end
    return tril(A), tril_pullback
end


text:
"""
     uniform_partition(n::Int, block_size::Int)

Compute a uniform block partition of the given size.

### Input

- `n`          -- number of dimensions of the partition
- `block_size` -- size of each block

### Output

A vector of ranges, `Vector{UnitRange{Int}}`, such that the size of each block
is the same, if possible.

### Examples

If the number of dimensions `n` is 2, we have two options: either two blocks
of size `1` or one block of size `2`:

```jldoctest partition
julia> LazySets.Approximations.uniform_partition(2, 1)
2-element Vector{UnitRange{Int64}}:
 1:1
 2:2

julia> LazySets.Approximations.uniform_partition(2, 2)
1-element Vector{UnitRange{Int64}}:
 1:2
```

If the block size argument is not compatible with (i.e. does not divide) `n`, the
output is filled with one block of the size needed to reach `n`:

```jldoctest partition
julia> LazySets.Approximations.uniform_partition(3, 1)
3-element Vector{UnitRange{Int64}}:
 1:1
 2:2
 3:3

julia> LazySets.Approximations.uniform_partition(3, 2)
2-element Vector{UnitRange{Int64}}:
 1:2
 3:3

julia> LazySets.Approximations.uniform_partition(10, 6)
2-element Vector{UnitRange{Int64}}:
 1:6
 7:10
```
"""
@inline function uniform_partition(n::Int, block_size::Int)
    m = div(n, block_size)
    r = n % block_size
    res = Vector{UnitRange{Int}}(undef, r > 0 ? m + 1 : m)
    k = 1
    @inbounds for i in 1:m
        l = k + block_size - 1
        res[i] = k:l
        k = l + 1
    end
    if r > 0
        res[m+1] = k:n
    end
    return res
end

"""
    decompose(S::LazySet{N},
              partition::AbstractVector{<:AbstractVector{Int}},
              block_options
             ) where {N}

Decompose a high-dimensional set into a Cartesian product of overapproximations
of the projections over the specified subspaces.

### Input

- `S`             -- set
- `partition`     -- vector of blocks (i.e., of vectors of integers) (see the
                     Notes below)
- `block_options` -- mapping from block indices in `partition` to a
                     corresponding overapproximation option; we only require
                     access via `[⋅]` (but see also the Notes below)

### Output

A `CartesianProductArray` containing the low-dimensional approximated
projections.

### Algorithm

For each block a specific `project` method is called, dispatching on the
corresponding overapproximation option.

### Notes

The argument `partition` requires some discussion.
Typically, the list of blocks should form a partition of the set
``\\{1, \\dots, n\\}`` represented as a list of consecutive blocks, where ``n``
is the ambient dimension of set `S`.

However, technically there is no problem if the blocks are not consecutive,
blocks are missing, blocks occur more than once, or blocks are overlapping.
This function will, however, stick to the order of blocks, so the resulting set
must be interpreted with care in such cases.
One use case is the need of a projection consisting of several blocks.

For convenience, the argument `block_options` can also be given as a single
option instead of a mapping, which is then interpreted as the option for all
blocks.

### Examples

This function supports different options: one can specify the target set, the
degree of accuracy, and template directions.
These options are exemplified below, where we use the following example.

```jldoctest decompose_examples
julia> S = Ball2(zeros(4), 1.0);  # set to be decomposed (4D 2-norm unit ball)

julia> P2d = [1:2, 3:4];  # a partition with two blocks of size two

julia> P1d = [[1], [2], [3], [4]];  # a partition with four blocks of size one
```

#### Different set types

We can decompose using polygons in constraint representation:

```jldoctest decompose_examples
julia> all(ai isa HPolygon for ai in array(decompose(S, P2d, HPolygon)))
true
```

For decomposition into 1D subspaces, we can use `Interval`:

```jldoctest decompose_examples
julia> all(ai isa Interval for ai in array(decompose(S, P1d, Interval)))
true
```

However, if you need to specify different set types for different blocks, the
interface presented so far does not apply.
See the paragraph *Advanced input for different block approximations* below for
how to do that.

#### Refining the decomposition I: ``ε``-close approximation

The ``ε`` option can be used to refine a decomposition, i.e., obtain a more
accurate result.
We use the [Iterative refinement](@ref) algorithm from the `Approximations`
module.

To illustrate this, consider again the set `S` from above.
We decompose into two 2D polygons.
Using smaller ``ε`` implies a better precision, thus more constraints in each 2D
decomposition.
In the following example, we look at the number of constraints in the first
block.

```jldoctest decompose_examples
julia> d(ε, bi) = array(decompose(S, P2d, (HPolygon => ε)))[bi]
d (generic function with 1 method)

julia> [length(constraints_list(d(ε, 1))) for ε in [Inf, 0.1, 0.01]]
3-element Vector{Int64}:
  4
  8
 32
```

#### Refining the decomposition II: template polyhedra

Another way to refine a decomposition is by using template polyhedra.
The idea is to specify a set of template directions and then to compute on each
block the polytopic overapproximation obtained by evaluating the support
function of the given input set over the template directions.

For example, octagonal 2D approximations of the set `S` are obtained with:

```jldoctest decompose_examples
julia> B = decompose(S, P2d, OctDirections);

julia> length(B.array) == 2 && all(dim(bi) == 2 for bi in B.array)
true
```

See [Template directions](@ref) for the available template directions.
Note that, in contrast to the polygonal ``ε``-close approximation from above,
this method can be applied to blocks of any size.

```jldoctest decompose_examples
julia> B = decompose(S, [1:4], OctDirections);

julia> length(B.array) == 1 && dim(B.array[1]) == 4
true
```

#### Advanced input for different block approximations

Instead of defining the approximation option uniformly for each block, we can
define different approximations for different blocks.
The third argument has to be a mapping from block index (in the partition) to
the corresponding approximation option.

For example:

```jldoctest decompose_examples
julia> res = array(decompose(S, P2d, Dict(1 => Hyperrectangle, 2 => 0.1)));

julia> typeof(res[1]), typeof(res[2])
(Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}, HPolygon{Float64, Vector{Float64}})
```
"""
function decompose(S::LazySet{N},
                   partition::AbstractVector{<:AbstractVector{Int}},
                   block_options
                  ) where {N}
    n = dim(S)
    result = Vector{LazySet{N}}(undef, length(partition))

    @inbounds for (i, block) in enumerate(partition)
        result[i] = project(S, block, block_options[i], n)
    end
    return CartesianProductArray(result)
end

# convenience method with uniform block options
function decompose(S::LazySet{N},
                   partition::AbstractVector{<:AbstractVector{Int}},
                   block_options::Union{Type{<:LazySet},
                                        Pair{<:UnionAll, <:Real},
                                        Real,
                                        Type{<:AbstractDirections},
                                        Nothing
                                       }
                  ) where {N}
    n = dim(S)
    result = Vector{LazySet{N}}(undef, length(partition))

    @inbounds for (i, block) in enumerate(partition)
        result[i] = project(S, block, block_options, n)
    end
    return CartesianProductArray(result)
end

# convenience method with uniform block size
function decompose(S::LazySet, block_options; block_size::Int=1)
    partition = uniform_partition(dim(S), block_size)
    return decompose(S, partition, block_options)
end


text:
struct PoissonProcess{T <: Real}
    intensity :: Union{T, Function}
    maxtime :: T
end

struct HomogPoissonProcess{T <: Real}
    intensity :: T
    maxtime :: T
end

struct HawkesProcess{T <: Real}
    α :: T
    β :: T
    μ :: T
    maxtime
end

struct EventStreamProcess{T <: Real, L}
    other_events :: Vector{Tuple{T, Int}}
    labels :: Vector{L}
    maxtime :: T
    basis :: BSplineBasis
    memory ::T
    λ_0 ::T
    other_kernels :: Array{Spline{BSplineBasis{Vector{T}}, Vector{T}}, 1}
    self_kernel :: Spline
end

function other_intensity(p :: EventStreamProcess, t)
    label_order = Dict(l => i for (i, l) in enumerate(p.labels))
    could_influence = searchsorted(p.other_events, prevfloat(t), by=(tup) -> memorylengths_away(tup[1], t, p.memory))
    influences = [(t - e[1], label_order[e[2]]) for e in p.other_events[could_influence]]
    total_intensity = sum([p.other_kernels[infl[2]](infl[1]) for infl in influences])
    return total_intensity
end

function rand(p :: EventStreamProcess)
    #relevant_spike_counts = zeros(length(p.labels) + 1)
    points = eltype(p.maxtime)[]
    kern_ubs = [max(0, maximum(spl.coeffs)) for spl in p.other_kernels]
    self_ub = max(0, maximum(p.self_kernel.coeffs))
    s=0.0
    earliestmemevnt_idx = 1
    rej_count=0
    while s < p.maxtime
        flush(stdout)
        # compute intensity UB
        intensity_ub = p.λ_0
        intensity_ub_contribs = 0.0
        next_event_time = p.maxtime
        for (e, n) in view(p.other_events, earliestmemevnt_idx:length(p.other_events))
            # increment by maximum possible impact of preceeding events
            if e > s
                next_event_time = e
                break
            end
            intensity_ub_contribs += kern_ubs[n]
        end
        for selfp in Iterators.reverse(points)
            s - selfp < p.memory || break
            intensity_ub_contribs += self_ub
        end
        intensity_ub *= exp(intensity_ub_contribs)
        print("s: $s, ub:$intensity_ub   \r")
        #println("s: $s, earliestmem: $earliestmemevnt_idx")
        u, d = rand(Float64, 2)
        w = -log(u)/intensity_ub # Exponential with mean 1/intensity_ub
        s = s + w
        s <= p.maxtime || break
        # If no sample until next event (which changes UB), then no need to reject/adopt
        if s > next_event_time
            s = next_event_time
            #update index of earliest memory
            for (e, n) in view(p.other_events, earliestmemevnt_idx:length(p.other_events))
                if s- e > p.memory
                    earliestmemevnt_idx +=1
                else
                    rej_count +=1
                    break
                end
            end
        else
            other_infl = 0.0
            for (e, n) in view(p.other_events, earliestmemevnt_idx:length(p.other_events))
                if s - e > p.memory
                    earliestmemevnt_idx += 1
                    continue
                end
                e <= s || break
                other_infl += p.other_kernels[n](s-e)
            end
            self_infl = 0.0
            for selfp in Iterators.reverse(points)
                s - selfp < p.memory || break
                self_infl += p.self_kernel(s - selfp)
            end
            intensity = p.λ_0 * exp(other_infl + self_infl)
            #println("s:$s  Intensity: $intensity, UB: $intensity_ub")
            if intensity > intensity_ub +0.01
                @warn "Intensity Upper bound is violated!"
                return points
            end
            if d <= intensity/intensity_ub
                push!(points, s)
            else
                rej_count +=1
            end
        end
    end
    println("$(length(points)/rej_count)")
    return points
end

function rand(p :: EventStreamProcess, intensity_ub::Real)
    # Hopefully can remove specification of intensity_ub, as difficult to specify
    # TODO: switch to Orgata, use fact that upper-bound for each spline is max(coefs)
    points = eltype(p.maxtime)[]
    s = 0
    while s < p.maxtime
        u, d = rand(Float64, 2) # Uniform(0, 1)
        w = -log(u)/intensity_ub # Exponential with mean 1/intensity_ub
        s = s + w # Proposed point
        s <=p.maxtime || break
        influential_diffs = eltype(p.maxtime)[]
        for x in Iterators.reverse(points)
            s - x < p.memory || break
            push!(influential_diffs, s - x)
        end
        int_value = other_intensity(p, s) + sum(p.self_kernel.(influential_diffs))
        int_value = p.λ_0 * exp.(int_value)
        if int_value > intensity_ub
            @warn "Intensity Upper Bound incorrectly specified: λ($s) = $int_value > $intensity_ub"
        end
        if d <= int_value/intensity_ub
            # accept with prop int_value/intensity_ub
            push!(points, s)
        end
        #ow reject
    end
    return points
end

function rand(p :: HawkesProcess)
    # Ogata Iterative Thinning
    points = eltype(p.α)[]
    s = 0
    while s < p.maxtime
        local_ub = p.μ + sum(p.α .* exp.(-1*p.β*(s .- points))) # Local UB on conditional intensity function, assuming no more points
        u, d = rand(Float64, 2)
        w = -log(u)/local_ub
        s= s + w
        s <=p.maxtime || break
        int_value = p.μ + sum(p.α .* exp.(-1*p.β*(s .- points)))
        if int_value > local_ub
            @warn "Intensity Upper Bound incorrectly specified: λ($s) = $int_value > $local_ub"
        end
        if d <= int_value/local_ub
            push!(points, s)
        end
    end
    return points
end

function rand(p :: HomogPoissonProcess)
    points = eltype(p.intensity)[]
    s=0
    while s<p.maxtime
        u = rand(Float64)
        w = -log(u)/p.intensity # Exponential wiht mean 1/intensity
        s = s+w
        s <= p.maxtime || break
        push!(points, s)
    end
    return points
end

function rand(p :: PoissonProcess, intensity_ub )
    # Shedler-Lewis Thinning
    points = eltype(intensity_ub)[]
    s = 0
    while s<p.maxtime
        u, d = rand(Float64, 2) # Uniform(0, 1)
        w = -log(u)/intensity_ub # Exponential with mean 1/intensity_ub
        s = s + w # Proposed point
        s <=p.maxtime || break
        int_value = p.intensity(s) :: Float64
        if int_value > intensity_ub
            @warn "Intensity Upper Bound incorrectly specified: λ($s) = $int_value > $intensity_ub"
        end
        if d <= int_value/intensity_ub
            # accept with prop int_value/intensity_ub
            push!(points, s)
        end
        #ow reject
    end
    return points
end


text:
const imaginary_unit_str =  split(string(Complex(0,0)), "0")[end]

const ARB_STR_FLAGS                = UInt64(0)

const ARB_STR_MORE                 = UInt64(1)        # extra, maybe incorrect digs
const ARB_STR_NO_RADIUS            = UInt64(2)        # omit radius
const ARB_STR_RADIUS               = UInt64(4)
const ARB_STR_CONDENSE             = UInt64(16)
const ARB_STR_CONDENSE_36_Digits   = UInt64(50*16)
const ARB_STR_CONDENSE_72_Digits   = UInt64(72*16)

const NO_FLAGS  = UInt64(0)
const NO_RADIUS = ARB_STR_NO_RADIUS

@inline function digit_precision(bitprecision::Int)
    bitprecision = evincedbits(bitprecision)
    return maximin_digits(bitprecision)
end

function trimzeros(str::String)
    if occursin('e', str)
        str1, str2 = String.(split(str, 'e'))
        if occursin('.', str1)
            str1a, str1b = String.(split(str1, '.'))
            str1b = trimallzeros(str1b)
            if str1b === ""
                str1b = "0"
            end
            str1 = join((str1a, str1b), '.')
        end
        str = join((str1, str2), 'e')
    elseif occursin('.', str)
        str1a, str1b = String.(split(str, '.'))
        str1b = trimallzeros(str1b)
        if str1b === ""
            str1b = "0"
        end
        str = join((str1a, str1b), '.')
    end
    return str
end

function trimallzeros(str::String)
    m = n = length(str)
    (n === 0 || str[n] !== '0') && return str
    while n>1 && str[n] === '0'
        n -= 1
    end
    m === n ? str : str[1:n]
end

function string(x::Mag, maxdigits::Int = maximin_digits(30), flags::UInt = NO_FLAGS)
    y = ArbFloat{precision(x)}()
    ccall(@libarb(arf_set_mag), Cvoid, (Ref{ArbFloat}, Ref{Mag}), y, x)
    z = ArbReal{precision(x)}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, y)
    unsafestr = ccall(@libarb(arb_get_str), Cstring,
                      (Ref{ArbReal}, Clong, Culong), z, maxdigits, flags)
    str = deepcopy( unsafe_string(pointer(unsafestr)) )
    ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr)
    return str
end

function string(x::ArbFloat{P}; midpoint::Bool=false) where {P}
    prec = midpoint ? digits4bits(P) : digit_precision(P)
    flags = ARB_STR_NO_RADIUS
    return arbstring(x, prec, flags=flags)
end

stringall(x::ArbFloat{P}) where {P} = return arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=NO_RADIUS)

function arbstring(x::ArbFloat{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_RADIUS) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, x)
    unsafestr = ccall(@libarb(arb_get_str), Cstring,
                      (Ref{ArbReal}, Clong, Culong), z, maxdigits, flags)
    str = deepcopy( unsafe_string(pointer(unsafestr)) )
    str = trimzeros(str)
    ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr)
    return str
end

function string(x::ArbReal{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    prec = midpoint ? digits4bits(P) : digit_precision(P)
    flags = radius ? ARB_STR_RADIUS : ARB_STR_NO_RADIUS
    return arbstring(x, prec, flags=flags)
end

stringall(x::ArbReal{P}; radius::Bool=false) where {P} =
    radius ? arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_RADIUS) :
             arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_NO_RADIUS)

function arbstring(x::ArbReal{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P}
    unsafestr = ccall(@libarb(arb_get_str), Cstring,
                      (Ref{ArbReal}, Clong, Culong), x, maxdigits, flags)
    str = deepcopy( unsafe_string(pointer(unsafestr)) )
    str = trimzeros(str)
    ccall(@libflint(flint_free), Cvoid, (Cstring,), unsafestr)
    return str
end

function string(x::ArbComplex{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    prec = midpoint ? digits4bits(P) : digit_precision(P)
    flags = radius ? ARB_STR_RADIUS : ARB_STR_NO_RADIUS
    return arbstring(x, prec, flags=flags)
end

stringall(x::ArbComplex{P}; radius::Bool=false) where {P} =
    radius ? arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_RADIUS) :
             arbstring(x, floor(Int,log10(2)*precision(x))+10; flags=ARB_STR_NO_RADIUS)

function arbstring(x::ArbComplex{P}, maxdigits::Int=digit_precision(P); flags::UInt = NO_FLAGS) where {P}
    # rea, ima = real(x), imag(x)
    rea = ArbReal{P}()
    ima = ArbReal{P}()
    ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), rea, x)
    ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), ima, x)

    ima_isneg = signbit(ima)
    ima_abs = ima_isneg ? -ima : ima
    connection = ima_isneg ? " - " : " + "
    rea_str = arbstring(rea, maxdigits, flags=flags)
    ima_str = arbstring(abs(ima), maxdigits, flags=flags)
    str = string(rea_str, connection, ima_str, imaginary_unit_str)
    return str
end


text:
export NeedPiecewiseCorrection

# Computes samplex of real or complex Haar matrices of size nxn
#
# For beta=1,2,4, generates random orthogonal, unitary and symplectic matrices
# of uniform Haar measure.
# These matrices are distributed with uniform Haar measure over the
# classical orthogonal, unitary and symplectic groups O(N), U(N) and
# Sp(N)~USp(2N) respectively.
#
# The last parameter specifies whether or not the piecewise correction
# is applied to ensure that it truly of Haar measure
# This addresses an inconsistency in the Householder reflections as
# implemented in most versions of LAPACK
# Method 0: No correction
# Method 1: Multiply rows by uniform random phases
# Method 2: Multiply rows by phases of diag(R)
# References:
#    Edelman and Rao, 2005
#    Mezzadri, 2006, math-ph/0609050
#TODO implement O(n^2) method
function rand(W::Haar, n::Integer, doCorrection::Integer)
    beta = W.beta
    M=rand(Ginibre(beta,n))
    q,r=qr(M)
    if doCorrection==0
        q
    elseif doCorrection==1
        if beta==1
            L = sign(rand(n).-0.5)
        elseif beta==2
            L = exp(im*rand(n)*2pi)
        elseif beta==4
            L = exp(im*rand(2n)*2pi)
        else
            error(string("beta = ",beta, " not implemented."))
        end
        q*diagm(L)
    elseif doCorrection==2
        if beta==1
            L=sign(diag(r))
        elseif (beta==2 || beta==4)
            L=diag(r)
            L=L./abs(L)
        else
            error(string("beta = ",beta, " not implemented."))
        end
        q*diagm(L)
    end
end

#By default, always do piecewise correction
#For most applications where you use the HaarMatrix as a similarity transform
#it doesn't matter, but better safe than sorry... let the user choose else
rand(W::Haar,n::Integer) = rand(W,n, 1)

#A utility method to check if the piecewise correction is needed
#This checks the R part of the QR factorization; if correctly done,
#the diagonals are all chi variables so are non-negative
function NeedPiecewiseCorrection()
    n=20
    R=qr(randn(Ginibre(2,n)))[2]
    return any([x<0 for x in diag(R)])
end

#TODO maybe, someday
#Haar measure on U(N) in terms of local coordinates
#Zyczkowski and Kus, Random unitary matrices, J. Phys. A: Math. Gen. 27,
#4235–4245 (1994).


text:
using GBIF
using Test

global anyerrors = false

tests = Dict{String,String}(
  "single occurrence functions" => "occurrence.jl",
  "multiple occurrences" => "occurrences.jl",
  "paging" => "paging.jl",
  "iteration" => "iteration.jl",
  "methods" => "methods.jl",
  "filtering" => "filter.jl",
  "species retrieval" => "taxon.jl"
)

for (name,test) in tests
  try
    include(test)
    println("\033[1m\033[32m✓\033[0m\t$(name)")
  catch e
    global anyerrors = true
    println("\033[1m\033[31m×\033[0m\t$(name)")
    println("\033[1m\033[38m→\033[0m\ttest/$(test)")
    showerror(stdout, e, backtrace())
    println()
    break
  end
end

if anyerrors
  throw("Tests failed")
end


text:
# mpiexec.exe -n 4 julia .\mpisum.jl
using ADCME

mpi_init()
r = mpi_rank()
a = constant(Float64.(Array(1:10) * r))
b = mpi_sum(a)

L = sum(b)
g = gradients(L, a)
sess = Session(); init(sess)
v, G = run(sess, [b,g])

@info r, G
if r==0
    ref = zeros(10)
    for k = 0:mpi_size()-1
        global ref
        ref += Array(1:10)*k
    end
    @info v, ref
end

mpi_finalize()


text:
using CoolPkg
using Test

@testset "CoolPkg.jl" begin
    @test add2(1,2) == 3
end


text:
"""
Represents everything needed to run and record a single simulation, including model, initial conditions, and metadata.

A vector of `Sim` objects can be executed with [`run`](@ref) or [`run_parallel`](@ref).

## Keyword Arguments
- `rng::AbstractRNG=Random.GLOBAL_RNG`
- `max_steps::Int=typemax(Int)`
- `simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps)`
- `metadata::NamedTuple a named tuple (or dictionary) of metadata for the sim that will be recorded, e.g. `(solver_iterations=500,)`.
"""
abstract type Sim end

struct POMDPSim <: Sim
    simulator::Simulator
    pomdp::POMDP
    policy::Policy
    updater::Updater
    initial_belief::Any
    initialstate::Any
    metadata::NamedTuple
end

problem(sim::POMDPSim) = sim.pomdp

struct MDPSim <: Sim
    simulator::Simulator
    mdp::MDP
    policy::Policy
    initialstate::Any
    metadata::NamedTuple
end

problem(sim::MDPSim) = sim.mdp

"""
    Sim(p::POMDP, policy::Policy, metadata=(note="a note",))
    Sim(p::POMDP, policy::Policy[, updater[, initial_belief[, initialstate]]]; kwargs...)

Create a `Sim` object that represents a POMDP simulation.
"""
function Sim(pomdp::POMDP,
             policy::Policy,
             up=updater(policy),
             initial_belief=initialstate_distribution(pomdp),
             initialstate=nothing;
             rng::AbstractRNG=Random.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata = NamedTuple()
            )

    if initialstate == nothing && statetype(pomdp) != Nothing
        is = rand(rng, initial_belief)
    else
        is = initialstate
    end
    return POMDPSim(simulator, pomdp, policy, up, initial_belief, is, merge(NamedTuple(), metadata))
end

"""
    Sim(p::MDP, policy::Policy, metadata=(note="a note",))
    Sim(p::MDP, policy::Policy[, initialstate]; kwargs...)

Create a `Sim` object that represents a MDP simulation.
"""
function Sim(mdp::MDP,
             policy::Policy,
             initialstate=nothing;
             rng::AbstractRNG=Random.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata = NamedTuple()
            )

    if initialstate == nothing && statetype(mdp) != Nothing
        is = POMDPs.initialstate(mdp, rng)
    else
        is = initialstate
    end
    return MDPSim(simulator, mdp, policy, is, merge(NamedTuple(), metadata))
end

POMDPs.simulate(s::POMDPSim) = simulate(s.simulator, s.pomdp, s.policy, s.updater, s.initial_belief, s.initialstate)
POMDPs.simulate(s::MDPSim) = simulate(s.simulator, s.mdp, s.policy, s.initialstate)

default_process(s::Sim, r::Float64) = (reward=r,)
default_process(s::Sim, hist::SimHistory) = default_process(s, discounted_reward(hist))

run_parallel(queue::AbstractVector; kwargs...) = run_parallel(default_process, queue; kwargs...)

"""
    run_parallel(queue::Vector{Sim})
    run_parallel(f::Function, queue::Vector{Sim})

Run `Sim` objects in `queue` in parallel and return results as a `DataFrame`.

By default, the `DataFrame` will contain the reward for each simulation and the metadata provided to the sim.

# Arguments
- `queue`: List of `Sim` objects to be executed
- `f`: Function to process the results of each simulation
This function should take two arguments, (1) the `Sim` that was executed and (2) the result of the simulation, by default a `SimHistory`. It should return a named tuple that will appear in the dataframe. See Examples below.

## Keyword Arguments
- `progress`: a `ProgressMeter.Progress` for showing progress through the simulations; `progress=false` will suppress the progress meter

# Examples

```julia
run_parallel(queue) do sim, hist
    return (n_steps=n_steps(hist), reward=discounted_reward(hist))
end
```
will return a dataframe with with the number of steps and the reward in it.
"""
function run_parallel(process::Function, queue::AbstractVector;
                      progress=Progress(length(queue), desc="Simulating..."),
                      proc_warn=true)

    #=
    frame_lines = pmap(progress, queue) do sim
        result = simulate(sim)
        return process(sim, result)
    end
    =#

    np = nprocs()
    if np == 1 && proc_warn
        @warn("""
             run_parallel(...) was started with only 1 process, so simulations will be run in serial.

             To supress this warning, use run_parallel(..., proc_warn=false).

             To use multiple processes, use addprocs() or the -p option (e.g. julia -p 4).
             """)
    end
    n = length(queue)
    i = 1
    prog = 0
    # based on the simple implementation of pmap here: https://docs.julialang.org/en/latest/manual/parallel-computing
    frame_lines = Vector{Any}(missing, n)
    nextidx() = (idx=i; i+=1; idx)
    prog_lock = ReentrantLock()
    @sync begin
        for p in 1:np
            if np == 1 || p != myid()
                @async begin
                    while true
                        idx = nextidx()
                        if idx > n
                            break
                        end
                        frame_lines[idx] = remotecall_fetch(p, queue[idx]) do sim
                            result = simulate(sim)
                            output = process(sim, result)
                            return merge(sim.metadata, output)
                        end
                        if progress isa Progress
                            lock(prog_lock)
                            update!(progress, prog+=1)
                            unlock(prog_lock)
                        end
                    end
                end
            end
        end
    end
    if progress isa Progress
        lock(prog_lock)
        finish!(progress)
        unlock(prog_lock)
    end

    return create_dataframe(frame_lines)
end

Base.run(queue::AbstractVector) = run(default_process, queue)

"""
    run(queue::Vector{Sim})
    run(f::Function, queue::Vector{Sim})

Run the `Sim` objects in `queue` on a single process and return the results as a dataframe.

See `run_parallel` for more information.
"""
function Base.run(process::Function, queue::AbstractVector; show_progress=true)
    lines = []
    if show_progress
        @showprogress for sim in queue
            result = simulate(sim)
            output = process(sim, result)
            line = merge(sim.metadata, output)
            push!(lines, line)
        end
    else
        for sim in queue
            result = simulate(sim)
            output = process(sim, result)
            line = merge(sim.metadata, output)
            push!(lines, line)
        end
    end
    return create_dataframe(lines)
end

function create_dataframe(lines::Vector)
    master = Dict{Symbol, AbstractVector}()
    for line in lines
        push_line!(master, line)
    end
    return DataFrame(master)
end

function push_line!(d::Dict{Symbol, AbstractVector}, line::NamedTuple)
    if isempty(d)
        len = 0
    else
        len = length(first(values(d)))
    end
    for (key, val) in pairs(line)
        T = typeof(val)
        if !haskey(d, key)
            d[key] = Vector{Union{T,Missing}}(missing, len)
        end
        data = d[key]
        if !isa(val,eltype(data))
            d[key] = convert(Array{promote_type(typeof(val), eltype(data)),1}, data)
        end
        push!(d[key], val)
    end
    for da in values(d)
        if length(da) < len + 1
            push!(da, missing)
        end
    end
    return d
end


text:
using DynamicGrids: normalise, minval, maxval, ismasked, rgb

"""
    RegionOutput(init; tspan, objective, [mask], [extent])

A minimal low-memory output that stores the inhabited regions for
each timestep, as required by the [`RegionObjective`](@ref).
"""
mutable struct RegionOutput{T,F<:AbstractVector{T},E,O} <: Output{T}
    frames::F
    running::Bool
    extent::E
    objective::O
end
RegionOutput(; frames, running, extent, objective, kwargs...) = begin
    predictions = [BitArray(zeros(Bool, size(objective.occurance)))]
    running = false
    RegionOutput(predictions, running, extent, objective)
end

objective(o::RegionOutput) = o.objective

# RegionObjective and supporting types/methods

DynamicGrids.storeframe!(o::RegionOutput, data::DynamicGrids.SimData) = begin
    f = DynamicGrids.frameindex(o, data)
    obj = objective(o)
    step = stepfromframe(obj, f)
    pred = predictions(obj, o)
    for I in CartesianIndices(first(DynamicGrids.init(data)))
        storeprediction!(first(data), o, step, pred, I)
    end
end

DynamicGrids.initgrids!(o::RegionOutput, init::AbstractArray) = begin
    obj = objective(o)
    step = stepfromframe(obj, 1)
    pred = predictions(obj, o)
    pred .= false
    for I in CartesianIndices(init)
        storeprediction!(init, o, step, pred, I)
    end
end
DynamicGrids.initgrids!(o::RegionOutput, init::NamedTuple) =
    DynamicGrids.initgrids!(o, first(init))

"""
Set region presence status in non-zero blocks
"""
@inline storeprediction!(grid, o::RegionOutput, step, predictions, I) = begin
    obj = objective(o)
    grid[I] > obj.detectionthreshold || return
    region = obj.regionlookup[I]
    region > zero(region) || return
    predictions[region, step] = true
    return
end

"""
    ColorRegionFit(objective, truescheme, falsescheme, falsepositivecolor,
                   truenegativecolor, falsenegativecolor, maskcolor)

An image procesor for visualising the match between predictions and observed
regional occupancy. A live version of the region fitting optimiser.

## Arguments:

- `objective` : A [`RegionObjective`](@ref) object
- `truescheme` : ColorSchemes.jl scheme or `Greyscale()` for true positive fit, etc.
- `falsescheme` : ColorSchemes.jl scheme or `Greyscale()` for true positive fit, etc.
- `truezerocolor` : `Color`, `Real` beteeen 0.0 and 1.0 or 3 tuple of `Real`.
- `falsezerocolor` : `Color`, `Real` beteeen 0.0 and 1.0 or 3 tuple of `Real`.
- `maskcolor` : `Color` or `Real` between 0.0 and 1.0 to use when a cell is masked.
"""
struct ColorRegionFit{O,T,F,TZ,FZ,M} <: GridProcessor
    objective::O
    truescheme::T
    falsescheme::F
    truezerocolor::TZ
    falsezerocolor::FZ
    maskcolor::M
end

DynamicGrids.grid2image(p::ColorRegionFit, o::ImageOutput, ruleset::Ruleset, grids::NamedTuple, t) = begin
    step = stepfromframe(p.objective, t)
    grid = first(grids)
    img = fill(RGB24(0), size(grid))
    obj = p.objective
    min, max = minval(o), maxval(o)
    for i in CartesianIndices(img)
        region = p.objective.regionlookup[i]
        img[i] = if !(p.maskcolor isa Nothing) && ismasked(mask(o), i)
            rgb(p.maskcolor)
        elseif region > zero(region)
            x = grid[i]
            normed = normalise(x, min, max)
            if p.objective.occurance[region, step]
                if !(p.truezerocolor isa Nothing) && normed == zero(normed)
                    rgb(p.falsezerocolor)
                else
                    rgb(p.truescheme, normed)
                end
            else
                if !(p.falsezerocolor isa Nothing) && normed == zero(normed)
                    rgb(p.truezerocolor)
                elseif x > obj.detectionthreshold
                    rgb(p.falsescheme, normed)
                else
                    rgb(p.truescheme, normed)
                end
            end
        else
            rgb(p.maskcolor)
        end
    end
    img
end

stepfromframe(objective::RegionObjective, t) =
    stepfromframe(objective.framesperstep, objective.start, t)
stepfromframe(framesperstep, start, t) =
    (t - 2one(t) + start) ÷ framesperstep + one(t)


text:
# Implicit solvent models
# Based on the OpenMM source code

export
    AbstractGBSA,
    ImplicitSolventOBC,
    born_radii_and_grad

"""
Generalized Born (GB) implicit solvent models augmented with the
hydrophobic solvent accessible surface area (SA) term.
Custom GBSA methods should sub-type this type.
"""
abstract type AbstractGBSA end

"""
    ImplicitSolventOBC(atoms, atoms_data, bonds)

Onufriev-Bashford-Case GBSA model.
Should be used along with a Coulomb or CoulombReactionField interaction.
The keyword argument `use_OBC2` determines whether to use parameter set
I (`false`, the default) or II (`true`).
"""
struct ImplicitSolventOBC{T, D, V, S, F, I} <: AbstractGBSA
    offset_radii::V
    scaled_offset_radii::V
    solvent_dielectric::T
    solute_dielectric::T
    offset::D
    cutoff::D
    use_ACE::Bool
    α::T
    β::T
    γ::T
    probe_radius::D
    sa_factor::S
    pre_factor::F
    is::I
    js::I
end

# Default solvent dielectric is 78.5 for consistency with AMBER
# Elsewhere it is 78.3
function ImplicitSolventOBC(atoms::AbstractArray{Atom{T, M, D, E}},
                            atoms_data,
                            bonds;
                            solvent_dielectric=78.5,
                            solute_dielectric=1.0,
                            offset=0.009u"nm",
                            cutoff=0.0u"nm",
                            probe_radius=0.14u"nm",
                            sa_factor=28.3919551u"kJ * mol^-1 * nm^-2",
                            use_ACE=true,
                            use_OBC2=false) where {T, M, D, E}
    default_radius = 0.15u"nm"
    element_to_radius = Dict(
        "N"  => 0.155u"nm",
        "O"  => 0.15u"nm" ,
        "F"  => 0.15u"nm" ,
        "Si" => 0.21u"nm" ,
        "P"  => 0.185u"nm",
        "S"  => 0.18u"nm" ,
        "Cl" => 0.17u"nm" ,
        "C"  => 0.17u"nm" ,
    )
    default_screen = 0.8
    element_to_screen = Dict(
        "H" => 0.85,
        "C" => 0.72,
        "N" => 0.79,
        "O" => 0.85,
        "F" => 0.88,
        "P" => 0.86,
        "S" => 0.96,
    )

    # Find atoms bonded to nitrogen
    n_atoms = length(atoms)
    atoms_bonded_to_N = falses(n_atoms)
    for (i, j) in zip(bonds.is, bonds.js)
        if atoms_data[i].element == "N"
            atoms_bonded_to_N[j] = true
        end
        if atoms_data[j].element == "N"
            atoms_bonded_to_N[i] = true
        end
    end

    offset_radii = D[]
    scaled_offset_radii = D[]
    for (at_data, bonded_to_N) in zip(atoms_data, atoms_bonded_to_N)
        if at_data.element in ("H", "D")
            radius = bonded_to_N ? 0.13u"nm" : 0.12u"nm"
        else
            radius = get(element_to_radius, at_data.element, default_radius)
        end
        offset_radius = radius - offset
        screen = get(element_to_screen, at_data.element, default_screen)
        if dimension(D) == u"𝐋"
            push!(offset_radii, offset_radius)
            push!(scaled_offset_radii, screen * offset_radius)
        else
            push!(offset_radii, ustrip(offset_radius))
            push!(scaled_offset_radii, ustrip(screen * offset_radius))
        end
    end

    if use_OBC2
        # GBOBCII parameters
        α, β, γ = T(1.0), T(0.8), T(4.85)
    else
        # GBOBCI parameters
        α, β, γ = T(0.8), T(0.0), T(2.909125)
    end

    inds_i = hcat([collect(1:n_atoms) for i in 1:n_atoms]...)
    inds_j = permutedims(inds_i, (2, 1))

    if !iszero(solute_dielectric) && !iszero(solvent_dielectric)
        pre_factor = T(-coulombconst) * (1/T(solute_dielectric) - 1/T(solvent_dielectric))
    else
        pre_factor = zero(T(coulombconst))
    end

    if isa(atoms, CuArray)
        or = CuArray(offset_radii) # cu converts to Float32
        sor = CuArray(scaled_offset_radii)
        is, js = cu(inds_i), cu(inds_j)
    else
        or = offset_radii
        sor = scaled_offset_radii
        is, js = inds_i, inds_j
    end

    if dimension(D) == u"𝐋"
        return ImplicitSolventOBC{T, D, typeof(or), typeof(T(sa_factor)), typeof(pre_factor), typeof(is)}(
                    or, sor, solvent_dielectric, solute_dielectric, offset, cutoff,
                    use_ACE, α, β, γ, probe_radius, T(sa_factor), pre_factor, is, js)
    else
        return ImplicitSolventOBC{T, T, typeof(or), T, T, typeof(is)}(
                    or, sor, solvent_dielectric, solute_dielectric, ustrip(offset), ustrip(cutoff),
                    use_ACE, α, β, γ, ustrip(probe_radius), ustrip(sa_factor), ustrip(pre_factor), is, js)
    end
end

function born_radii_loop(coord_i, coord_j, ori, srj, cutoff, box_size)
    I = zero(coord_i[1] / unit(cutoff)^2)
    r = norm(vector(coord_i, coord_j, box_size))
    if iszero(r) || (!iszero(cutoff) && r > cutoff)
        return I
    end
    U = r + srj
    if ori < U
        D_ij = abs(r - srj)
        L = max(ori, D_ij)
        I += (1/L - 1/U + (r - (srj^2)/r)*(1/(U^2) - 1/(L^2))/4 + log(L/U)/(2*r)) / 2
        if ori < (srj - r)
            I += 2 * (1/ori - 1/L)
        end
    end
    return I
end

"""
    born_radii_and_grad(inter, coords, box_size)

Calculate Born radii and gradients of Born radii with respect to atomic distance.
Custom GBSA methods should implement this function.
"""
function born_radii_and_grad(inter::ImplicitSolventOBC, coords, box_size)
    coords_i = @view coords[inter.is]
    coords_j = @view coords[inter.js]
    oris = @view inter.offset_radii[inter.is]
    srjs = @view inter.scaled_offset_radii[inter.js]
    Is = dropdims(sum(born_radii_loop.(coords_i, coords_j, oris, srjs, (inter.cutoff,),
                                        (box_size,)); dims=2); dims=2)

    ori = inter.offset_radii
    radii = ori .+ inter.offset
    ψs = Is .* ori
    ψs2 = ψs .^ 2
    α, β, γ = inter.α, inter.β, inter.γ
    tanh_sums = tanh.(α .* ψs .- β .* ψs2 .+ γ .* ψs2 .* ψs)
    Bs = 1 ./ (1 ./ ori .- tanh_sums ./ radii)
    grad_terms = ori .* (α .- 2 .* β .* ψs .+ 3 .* γ .* ψs2)
    B_grads = (1 .- tanh_sums .^ 2) .* grad_terms ./ radii

    return Bs, B_grads
end

# Store the results of the ij broadcasts during force calculation
struct ForceLoopResult1{T, V}
    bi::T
    bj::T
    fi::V
    fj::V
end

struct ForceLoopResult2{V}
    fi::V
    fj::V
end

function gb_force_loop_1(coord_i, coord_j, i, j, charge_i, charge_j,
                            Bi, Bj, cutoff, pre_factor, box_size)
    if j < i
        zero_force = zero(pre_factor ./ coord_i .^ 2)
        return ForceLoopResult1(zero_force[1], zero_force[1], zero_force, zero_force)
    end
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)
    if !iszero(cutoff) && r2 > cutoff^2
        zero_force = zero(pre_factor ./ coord_i .^ 2)
        return ForceLoopResult1(zero_force[1], zero_force[1], zero_force, zero_force)
    end
    alpha2_ij = Bi * Bj
    D = r2 / (4 * alpha2_ij)
    exp_term = exp(-D)
    denominator2 = r2 + alpha2_ij * exp_term
    denominator = sqrt(denominator2)
    Gpol = (pre_factor * charge_i * charge_j) / denominator
    dGpol_dr = -Gpol * (1 - exp_term/4) / denominator2
    dGpol_dalpha2_ij = -Gpol * exp_term * (1 + D) / (2 * denominator2)
    change_born_force_i = dGpol_dalpha2_ij * Bj
    if i != j
        change_born_force_j = dGpol_dalpha2_ij * Bi
        fdr = dr * dGpol_dr
        change_fs_i =  fdr
        change_fs_j = -fdr
        return ForceLoopResult1(change_born_force_i, change_born_force_j,
                                change_fs_i, change_fs_j)
    else
        zero_force = zero(pre_factor ./ coord_i .^ 2)
        return ForceLoopResult1(change_born_force_i, zero_force[1], zero_force, zero_force)
    end
end

function gb_force_loop_2(coord_i, coord_j, bi, ori, srj, cutoff, box_size)
    dr = vector(coord_i, coord_j, box_size)
    r = norm(dr)
    if iszero(r) || (!iszero(cutoff) && r > cutoff)
        zero_force = zero(bi ./ coord_i .^ 2)
        return ForceLoopResult2(zero_force, zero_force)
    end
    rsrj = r + srj
    if ori < rsrj
        D = abs(r - srj)
        L = inv(max(ori, D))
        U = inv(rsrj)
        rinv = inv(r)
        r2inv = rinv^2
        t3 = (1 + (srj^2)*r2inv)*(L^2 - U^2)/8 + log(U/L)*r2inv/4
        de = bi * t3 * rinv
        fdr = dr * de
        return ForceLoopResult2(-fdr, fdr)
    else
        zero_force = zero(bi ./ coord_i .^ 2)
        return ForceLoopResult2(zero_force, zero_force)
    end
end

function forces(inter::AbstractGBSA, sys, neighbors=nothing)
    coords, atoms, box_size = sys.coords, sys.atoms, sys.box_size
    Bs, B_grads = born_radii_and_grad(inter, coords, box_size)

    if inter.use_ACE
        radii = inter.offset_radii .+ inter.offset
        sa_terms = inter.sa_factor .* (radii .+ inter.probe_radius) .^ 2 .* (radii ./ Bs) .^ 6
        born_forces = (-6 .* sa_terms ./ Bs) .* (Bs .> zero(inter.offset))
    else
        born_forces = zeros(typeof(inter.sa_factor * inter.cutoff), length(sys))
    end

    coords_i = @view coords[inter.is]
    coords_j = @view coords[inter.js]
    charges = charge.(atoms)
    charges_i = @view charges[inter.is]
    charges_j = @view charges[inter.js]
    Bsi = @view Bs[inter.is]
    Bsj = @view Bs[inter.js]
    loop_res_1 = gb_force_loop_1.(coords_i, coords_j, inter.is, inter.js, charges_i, charges_j,
                                    Bsi, Bsj, (inter.cutoff,), (inter.pre_factor,), (box_size,))
    born_forces = born_forces .+ dropdims(sum(lr -> lr.bi, loop_res_1; dims=2); dims=2)
    born_forces = born_forces .+ dropdims(sum(lr -> lr.bj, loop_res_1; dims=1); dims=1)
    fs = dropdims(sum(lr -> lr.fi, loop_res_1; dims=2); dims=2) .+ dropdims(sum(lr -> lr.fj, loop_res_1; dims=1); dims=1)

    born_forces_2 = born_forces .* (Bs .^ 2) .* B_grads

    bis = @view born_forces_2[inter.is]
    oris = @view inter.offset_radii[inter.is]
    srjs = @view inter.scaled_offset_radii[inter.js]
    loop_res_2 = gb_force_loop_2.(coords_i, coords_j, bis, oris, srjs, (inter.cutoff,), (box_size,))

    return fs .+ dropdims(sum(lr -> lr.fi, loop_res_2; dims=2); dims=2) .+ dropdims(sum(lr -> lr.fj, loop_res_2; dims=1); dims=1)
end

function gb_energy_loop(coord_i, coord_j, i, j, charge_i, charge_j, Bi, Bj, ori, cutoff,
                        pre_factor, offset, probe_radius, sa_factor, use_ACE, box_size)
    if i == j
        E = pre_factor * (charge_i^2) / (2*Bi)
        if use_ACE && (Bi > zero(offset))
            radius_i = ori + offset
            E += sa_factor * (radius_i + probe_radius)^2 * (radius_i / Bi)^6
        end
        return E
    elseif j > i
        r2 = sum(abs2, vector(coord_i, coord_j, box_size))
        if !iszero(cutoff) && r2 > cutoff^2
            return zero(pre_factor / offset)
        end
        f = sqrt(r2 + Bi*Bj*exp(-r2/(4*Bi*Bj)))
        if iszero(cutoff)
            f_cutoff = 1/f
        else
            f_cutoff = (1/f - 1/cutoff)
        end
        return pre_factor * charge_i * charge_j * f_cutoff
    else
        return zero(pre_factor / offset)
    end
end

function potential_energy(inter::AbstractGBSA, sys, neighbors=nothing)
    coords, atoms, box_size = sys.coords, sys.atoms, sys.box_size
    Bs, B_grads = born_radii_and_grad(inter, coords, box_size)

    coords_i = @view coords[inter.is]
    coords_j = @view coords[inter.js]
    charges = charge.(atoms)
    charges_i = @view charges[inter.is]
    charges_j = @view charges[inter.js]
    Bsi = @view Bs[inter.is]
    Bsj = @view Bs[inter.js]
    oris = @view inter.offset_radii[inter.is]
    return sum(gb_energy_loop.(coords_i, coords_j, inter.is, inter.js, charges_i, charges_j,
                                Bsi, Bsj, oris, (inter.cutoff,), (inter.pre_factor,),
                                (inter.offset,), (inter.probe_radius,), (inter.sa_factor,),
                                (inter.use_ACE,), (box_size,)))
end


text:
using Smg2s
using SparseArrays
using MatrixMarket
using LinearAlgebra

n = 10
diag_l = -5
diag_u = -3
nbone = 2

function f(idx::Integer, n = n)
    return 10 * cos(((idx-1) / 2) * 2 * π / n) + 5
end

spec = zeros(ComplexF64, n)
Spectrum!(spec, f, n)

genMat = nonsym(nbone, n, diag_l, diag_u, spec)

MatrixMarket.mmwrite("genMat.mtx", genMat)


text:
using DynamicGrids, Aqua, SafeTestsets

if VERSION >= v"1.5.0"
    # Amibiguities are not owned by DynamicGrids
    # Aqua.test_ambiguities([DynamicGrids, Base, Core])
    Aqua.test_unbound_args(DynamicGrids)
    Aqua.test_undefined_exports(DynamicGrids)
    Aqua.test_project_extras(DynamicGrids)
    # Aqua.test_stale_deps(DynamicGrids)
    Aqua.test_deps_compat(DynamicGrids)
    Aqua.test_project_toml_formatting(DynamicGrids)
end

@time @safetestset "generated" begin include("generated.jl") end
@time @safetestset "rules" begin include("rules.jl") end
@time @safetestset "neighborhoods" begin include("neighborhoods.jl") end
@time @safetestset "simulationdata" begin include("simulationdata.jl") end
@time @safetestset "utils" begin include("utils.jl") end
@time @safetestset "chain" begin include("chain.jl") end
@time @safetestset "condition" begin include("condition.jl") end
@time @safetestset "outputs" begin include("outputs.jl") end
@time @safetestset "transformed" begin include("transformed.jl") end
@time @safetestset "integration" begin include("integration.jl") end
@time @safetestset "objectgrids" begin include("objectgrids.jl") end
@time @safetestset "parametersources" begin include("parametersources.jl") end
@time @safetestset "show" begin include("show.jl") end
@time @safetestset "textconfig" begin include("textconfig.jl") end
# ImageMagick breaks in windows travis for some reason
if !Sys.iswindows()
    @time @safetestset "image" begin include("image.jl") end
end