text:
### Constraints for Thermal Generation without commitment variables ####

"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitmentconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64; args...) where {T <: PowerSystems.ThermalGen, D <: StandardThermalCommitment, S <: AbstractDCPowerModel}

    on_th = m[:on_th]
    start_th = m[:start_th]
    stop_th = m[:stop_th]

    name_index = m[:on_th].axes[1]
    time_index = m[:on_th].axes[2]

    # set args default values
    if :initialstatus in keys(args)
        initialstatus = args[:initialstatus]
    else
        initialstatus = Dict(zip(name_index,[devices[ix].tech.activepower > 0.0 ? 1 : 0  for (ix,name) in enumerate(name_index) if name == devices[ix].name ]))
    end

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    commitment_th  = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_periods), name_index, time_index)

    for (ix,name) in enumerate(name_index)
        if name == devices[ix].name
            t1 = time_index[1]
            commitment_th[name,t1] = @constraint(m, on_th[name,t1] == initialstatus[name] + start_th[name,t1] - stop_th[name,t1])
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    for t in time_index[2:end], (ix,name) in enumerate(name_index)
        if name == devices[ix].name
            commitment_th[name,t] = @constraint(m, on_th[name,t] == on_th[name,t-1] + start_th[name,t] - stop_th[name,t])
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    JuMP.register_object(m, :commitment_th, commitment_th)

    return m
end

function timeconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64; args...) where {T <: PowerSystems.ThermalGen, D <: StandardThermalCommitment, S <: AbstractDCPowerModel}

    devices = [d for d in devices if !isa(d.tech.timelimits, Nothing)]

    if !isempty(devices)
        # TODO: Change loop orders, loop over time first and then over the device names

        on_th = m[:on_th]
        start_th = m[:start_th]
        stop_th = m[:stop_th]
        time_index = m[:on_th].axes[2]
        name_index = [d.name for d in devices]

        # set args default values
        if :initialonduration in keys(args)
            initialonduration = args[:initialonduration]
        else
            initialonduration = Dict(zip(name_index,ones(Float64,length(devices))*9999))
        end

        if :initialoffduration in keys(args)
            initialoffduration = args[:initialoffduration]
        else
            initialoffduration = Dict(zip(name_index,ones(Float64,length(devices))*9999))
        end

        (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

        mindown_th = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_periods), name_index, time_index)
        minup_th = JuMP.JuMPArray(Array{ConstraintRef}(undef, length(name_index), time_periods), name_index, time_index)

        for t in time_index, (ix,name) in enumerate(name_index)
            if name==devices[ix].name
                if t - devices[ix].tech.timelimits.up >= 1
                    tst = devices[ix].tech.timelimits.up
                else
                    tst = max(0, devices[ix].tech.timelimits.up - initialonduration[name])
                end
                if t - devices[ix].tech.timelimits.down >= 1
                    tsd = devices[ix].tech.timelimits.down
                else
                    tsd = max(0, devices[ix].tech.timelimits.down - initialoffduration[name])
                end

                minup_th[name,t] = @constraint(m,sum([start_th[name,i] for i in ((t - tst - 1) :t) if i > 0 ]) <= on_th[name,t])
                mindown_th[name,t] = @constraint(m,sum([stop_th[name,i] for i in ((t - tsd - 1) :t) if i > 0]) <= (1 - on_th[name,t]))

            else
                @error "Bus name in Array and variable do not match"

            end
        end

        JuMP.register_object(m, :minup_th, minup_th)
        JuMP.register_object(m, :mindown_th, mindown_th)

    else
        @warn "There are no generators with Min-up -down limits data in the system"

    end

    return m
end

#TODO: Add the Knueven Model

ext: jl


text:
module BenchRips

using BenchmarkTools
using Ripserer

include(joinpath(@__DIR__, "utils.jl"))

full_benchmarks = [
    (file="o3_1024.pts", dim_max=3, threshold=1.8),
    (file="fract-r.dist", dim_max=2, threshold=nothing),
    (file="dragon2000.pts", dim_max=1, threshold=nothing),
]
sparse_benchmarks = [(file="alpha_4_sphere_2000.spdist", dim_max=4)]

suite = BenchmarkGroup()
suite["sparse"] = BenchmarkGroup()
suite["dense"] = BenchmarkGroup()

for b in full_benchmarks
    for sparse in (true, false)
        name = b.file
        data = load_data(joinpath(@__DIR__, "data", b.file))
        suite[sparse ? "sparse" : "dense"][name] = @benchmarkable ripserer(
            $data; dim_max=$(b.dim_max), threshold=$(b.threshold), sparse=$sparse
        ) samples = 2
    end
end

for b in sparse_benchmarks
    name = b.file
    data = load_data(joinpath(@__DIR__, "data", b.file))
    suite["sparse"][name] = @benchmarkable ripserer($data; dim_max=$(b.dim_max)) samples = 2
end

end
BenchRips.suite


text:
"""
`ClimbingImage(∇V, τ, ∇V_climb,integrate!, dist, Δt)` - Set up the Climbing Image Method

### Fields
* `∇V`   - In place gradient of the potential
* `τ` - Approximate tangent vector of the unstable direction
* `∇V_climb   `- In place gradient of the potential with reflector
* `integrate!` - In place integrator
* `Δt`    - Time step
"""
struct ClimbingImage{TGV, TT, TGVC, TI, TF<:AbstractFloat} <: ClimbingImageMethod
    ∇V::TGV
    τ::TT
    ∇V_climb::TGVC
    integrate!::TI
    Δt::TF
end

function ClimbingImage(∇V::TGV, τ::TT, integrate!::TI, Δt::TF) where{TGV, TT, TI, TF}

    function ∇V_climb(x)
        gradV = ∇V(x);
        gradV .-= 2 * (gradV⋅τ) * τ;
        return gradV
    end

    return ClimbingImage(∇V, τ, ∇V_climb, integrate!, Δt);
end

"""
`climbing_image`: Run the climbing image method...
### Fields
* `u` - Initial guess for saddle
* `C` - Climbing image data structure
### Optional Fields
"""
function climbing_image(u₀, C::TC; options=SaddleOptions()) where {TC <: ClimbingImage}
    u = copy(u₀);
    gradV = copy(u);

    u_trajectory = typeof(u)[copy(u)];

    err_est = 0.;

    for n in 1:options.nmax
        C.integrate!(u, C.∇V_climb, C.Δt);
        gradV = C.∇V(u);
        err_est =norm(gradV);
        if(options.verbose && mod(n, options.print_iters)==0)
            @printf("[%d]: error = %g\n", n, err_est);
        end

        if(options.save_trajectory)
            push!(u_trajectory, copy(u));
        end
        if (err_est< options.tol)
            if(!options.save_trajectory)
                push!(u_trajectory, copy(u));
            end
            break
        end
    end

    if(err_est >= options.tol)
        @printf("ERROR: Did not converge after %d iterartions", options.nmax)
    end
    return u_trajectory
end

"""
`climbing_image!`: Run the climbing image method...
### Fields
* `u` - Initial guess for saddle
* `C` - Climbing image data structure
### Optional Fields
"""
function climbing_image!(u, C::TC; options=SaddleOptions()) where {TC <: ClimbingImage}
    gradV = copy(u);

    err_est = 0.;

    for n in 1:options.nmax
        C.integrate!(u, C.∇V_climb, C.Δt);
        gradV = C.∇V(u);
        err_est =norm(gradV);
        if(options.verbose && mod(n, options.print_iters)==0)
            @printf("[%d]: error = %g\n", n, err_est);
        end
        if (err_est< options.tol)
            break
        end
    end

    if(err_est >= options.tol)
        @printf("ERROR: Did not converge after %d iterartions", options.nmax)
    end
    u
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/lower_problem.jl
# Functions which define how relaxations of subproblems are constructed,
# when domain reductions algorithms are run, and the lower bounding (relaxed
# problem solution subroutines).
#############################################################################

"""
$(FUNCTIONNAME)

Adds linear objective cut constraint to the `x.relaxed_optimizer`.
"""
function objective_cut!(m::GlobalOptimizer, check_safe::Bool)
    f = m._working_problem._objective_saf
    u = m._global_upper_bound
    if  u < Inf
        b = f.constant
        f.constant = 0.0
        if check_safe && is_safe_cut!(m, f)
            s = LT(u - b + _constraint_tol(m))
            m._affine_objective_cut_ci = MOI.add_constraint(_relaxed_optimizer(m), f, s)
        end
        f.constant = b
        m._new_eval_objective = false
    end
    return
end

"""
    RelaxResultStatus

Status code used internally to determine how to interpret the results from the
solution of a relaxed problem.
"""
@enum(RelaxResultStatus, RRS_OPTIMAL, RRS_DUAL_FEASIBLE, RRS_INFEASIBLE, RRS_INVALID)

"""
$(SIGNATURES)

Takes an `MOI.TerminationStatusCode` and a `MOI.ResultStatusCode` and returns
the tuple `(valid_result::Bool, feasible::Bool)`. The value `valid_result` is
`true` if the pair of codes prove that either the subproblem solution was solved
to global optimality or the subproblem solution is infeasible. The value of
`feasible` is true if the problem is feasible and false if the problem is infeasible.
"""
function relaxed_problem_status(t::MOI.TerminationStatusCode,
                                p::MOI.ResultStatusCode,
                                d::MOI.ResultStatusCode)

    if (t == MOI.OPTIMAL) && (p == MOI.FEASIBLE_POINT)
        return RRS_OPTIMAL
    elseif t == MOI.INFEASIBLE
        if (p == MOI.INFEASIBILITY_CERTIFICATE) ||
           (p == MOI.NO_SOLUTION) || (p == MOI.UNKNOWN_RESULT_STATUS)
            return RRS_INFEASIBLE
        end
    elseif (t == MOI.INFEASIBLE_OR_UNBOUNDED && p == MOI.NO_SOLUTION)
        return RRS_INFEASIBLE
    end
    (d == MOI.FEASIBLE_POINT) && return RRS_DUAL_FEASIBLE
    return RRS_INVALID
end

"""
$(SIGNATURES)

Updates the relaxed constraint by setting the constraint set of `v == x*`` ,
`xL_i <= x_i`, and `x_i <= xU_i` for each such constraint added to the relaxed
optimizer. Resets integral valued constraints to either `EqualTo` or `Interval`
constraints.
"""
function update_relaxed_problem_box!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    d = _relaxed_optimizer(m)
    for i = 1:_variable_num(BranchVar(), m)
        l = _lower_bound(BranchVar(), m, i)
        u = _upper_bound(BranchVar(), m, i)
        v = VI(_bvi(m, i))
        if l == u
            ci_vi_et = MOI.add_constraint(d, v, ET(l))
            push!(m._relaxed_variable_et, (ci_vi_et,i))
        else
            ci_vi_lt = MOI.add_constraint(d, v, LT(u))
            ci_vi_gt = MOI.add_constraint(d, v, GT(l))
            m._node_to_sv_leq_ci[i] = ci_vi_lt
            m._node_to_sv_geq_ci[i] = ci_vi_gt
            push!(m._relaxed_variable_lt, (ci_vi_lt,i))
            push!(m._relaxed_variable_gt, (ci_vi_gt,i))
        end
    end
    return
end

const SOLUTION_EPS = 0.05
function store_lower_solution!(m::GlobalOptimizer{R,S,Q}, d::T) where {R,S,Q<:ExtensionType,T}
    for i = 1:_variable_num(FullVar(), m)
        l = _lower_bound(FullVar(), m, i)
        u = _upper_bound(FullVar(), m, i)
        ladj = l + SOLUTION_EPS*(u - l)
        uadj = u - SOLUTION_EPS*(u - l)
        x = MOI.get(d, MOI.VariablePrimal(), m._relaxed_variable_index[i])
        (x < ladj) && (x = ladj)
        (x > uadj) && (x = uadj)
        m._lower_solution[i] = x
    end
    return
end

function reset_relaxation!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    d = _relaxed_optimizer(m)
    m._cut_iterations = 1
    m._obbt_performed_flag = false
    m._working_problem._relaxed_evaluator.is_first_eval = true
    fill!(m._working_problem._relaxed_evaluator.subexpressions_eval, false)

    m._new_eval_objective = true
    m._new_eval_constraint = true

    # delete added affine constraints
    foreach(c -> MOI.delete(d, c), m._affine_relax_ci)
    empty!(m._affine_relax_ci)

    # delete variable
    foreach(c -> MOI.delete(d, c[1]), m._relaxed_variable_et)
    foreach(c -> MOI.delete(d, c[1]), m._relaxed_variable_lt)
    foreach(c -> MOI.delete(d, c[1]), m._relaxed_variable_gt)
    foreach(c -> MOI.delete(d, c), m._relaxed_variable_integer)
    empty!(m._relaxed_variable_et)
    empty!(m._relaxed_variable_lt)
    empty!(m._relaxed_variable_gt)
    empty!(m._relaxed_variable_integer)

    # delete objective cut
    !isnothing(m._affine_objective_cut_ci) && MOI.delete(d, m._affine_objective_cut_ci)
    return
end

"""
$(FUNCTIONNAME)

"""
function set_first_relax_point!(m::GlobalOptimizer)
    if m._cut_iterations == 1
        m._working_problem._relaxed_evaluator.is_first_eval = true
        m._new_eval_constraint = true
        m._new_eval_objective = true
        for i = 1:_variable_num(FullVar(), m)
            l = _lower_bound(FullVar(), m, i)
            u = _upper_bound(FullVar(), m, i)
            if isfinite(l) && isfinite(u)
                x = 0.5*(l + u)
            elseif isfinite(l)
                x = min(0.0, u)
            elseif isfinite(u)
                x = max(0.0, l)
            else
                x = 0.0
            end
            _set_lower_solution!(FullVar(), m, x, i)
        end
    end
    return
end

"""
$(TYPEDSIGNATURES)

A routine that adds relaxations for all nonlinear constraints and quadratic constraints
corresponding to the current node to the relaxed problem. This adds an objective cut
(if specified by `objective_cut_on`) and then sets the `_new_eval_constraint` flag
to false indicating that an initial evaluation of the constraints has occurred. If
the `objective_cut_on` flag is `true` then the `_new_eval_objective` flag is also
set to `false` indicating that the objective expression was evaluated.
"""
function relax_all_constraints!(t::ExtensionType, m::GlobalOptimizer, k::Int)
    check_safe = (k == 1) ? false : m._parameters.cut_safe_on
    wp = m._working_problem
    wp._relaxed_evaluator.is_first_eval = m._new_eval_constraint
    foreach(f -> relax!(m, f, k, check_safe), wp._sqf_leq)
    foreach(f -> relax!(m, f, k, check_safe), wp._sqf_eq)
    valid_relax_flag = true
    num_feasible_relax_flag = true
    if valid_relax_flag
        for nl in wp._nonlinear_constr
            valid_cut, feas_cut = relax!(m, nl, k, check_safe)
            valid_relax_flag &= valid_cut
            num_feasible_relax_flag &= feas_cut
        end
    end
    m._new_eval_constraint = false
    (k == 1) && objective_cut!(m, check_safe)
    return valid_relax_flag, num_feasible_relax_flag
end
relax_constraints!(t::ExtensionType, m::GlobalOptimizer, k::Int) = relax_all_constraints!(t, m, k)
relax_constraints!(m::GlobalOptimizer{R,S,Q}, k::Int) where {R,S,Q<:ExtensionType} = relax_constraints!(_ext(m), m, k)

function relax_problem!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    wp = m._working_problem
    if m._cut_iterations == 1
        reset_relaxation!(m)
        if m._nonlinear_evaluator_created
            set_node!(wp._relaxed_evaluator, m._current_node)
            set_reference_point!(m)
            fill!(wp._relaxed_evaluator.subexpressions_eval, false)
        end
        wp._relaxed_evaluator.is_post = m._parameters.subgrad_tighten
        _set_has_value!(wp._objective, false)
        wp._relaxed_evaluator.interval_intersect = false
        update_relaxed_problem_box!(m)
        set_first_relax_point!(m)
    else
        set_reference_point!(m)
    end
    valid_relax_flag, num_feasible_relax_flag = relax_constraints!(m, m._cut_iterations)
    MOI.set(_relaxed_optimizer(m), MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return valid_relax_flag, num_feasible_relax_flag
end

"""
$(SIGNATURES)

Retrieves the lower and upper duals for variable bounds from the
`relaxed_optimizer` and sets the appropriate values in the
`_lower_lvd` and `_lower_uvd` storage fields.
"""
function set_dual!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    d = _relaxed_optimizer(m)
    if MOI.get(d, MOI.DualStatus()) == MOI.FEASIBLE_POINT
        for (c, i) in m._relaxed_variable_lt
            m._lower_uvd[i] = MOI.get(d, MOI.ConstraintDual(), c)
        end
        for (c, i) in m._relaxed_variable_gt
            m._lower_lvd[i] = MOI.get(d, MOI.ConstraintDual(), c)
        end
    else
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
    end
    return
end

"""
"""
function interval_objective_bound! end
interval_objective_bound!(m::GlobalOptimizer, f::Nothing, is_first_eval) = nothing
function interval_objective_bound!(m::GlobalOptimizer, f::AffineFunctionIneq, is_first_eval)
    m._working_problem._relaxed_evaluator.is_first_eval = is_first_eval
    fL, fU = bound_objective(m)
    if fL > m._lower_objective_value
        m._lower_objective_value = fL
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
        m._cut_add_flag = false
    end
end
function interval_objective_bound!(m::GlobalOptimizer, f, is_first_eval)
    m._working_problem._relaxed_evaluator.is_first_eval = is_first_eval
    if is_first_eval
        m._working_problem._relaxed_evaluator.pass_number = 1
    end
    fL, fU = bound_objective(m)
    fv = _is_input_min(m) ? fL : -fU
    if fv > m._lower_objective_value
        m._lower_objective_value = fv
        fill!(m._lower_lvd, 0.0)
        fill!(m._lower_uvd, 0.0)
        m._cut_add_flag = false
    end
    return
end
interval_objective_bound!(m::GlobalOptimizer, is_first_eval = true) = interval_objective_bound!(m, m._working_problem._objective, is_first_eval)

"""
$(TYPEDSIGNATURES)

Runs contractor methods prior to solving lower bounding problem. By default linear and quadratic
contractor methods followed by interval constraint propagation then optimization-based bound
tightening for a specified number of iterations while the subproblem at current node `n` has
not been proven infeasible.
"""
function preprocess!(t::ExtensionType, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}

    feasible_flag = true
    reset_relaxation!(m)
    if _fbbt_lp_depth(m) >= _iteration_count(m)
        load_fbbt_buffer!(m)
        for _ = 1:_fbbt_lp_repetitions(m)
            ns = NodeBB(_current_node(m))
            for f in m._working_problem._saf_leq
                !(feasible_flag = feasible_flag && fbbt!(m, f)) && break
            end
            !feasible_flag && break
            for f in m._working_problem._saf_eq
                !(feasible_flag = feasible_flag && fbbt!(m, f)) && break
            end
            (same_box(ns,_current_node(m),0.0) || !feasible_flag) && break
        end
        unpack_fbbt_buffer!(m)
    end

    # done after cp to prevent using cp specific flags in cut generation
    set_first_relax_point!(m)
    # nonlinear CP can detect infeasibility and bound objective even if
    # the relaxation is ill-posed, so one is always used to mitigate numerical issues
    cp_reps = _cp_depth(m) >= _iteration_count(m) ? _cp_repetitions(m) : 1
    for _ = 1:_cp_repetitions(m)
        ns = NodeBB(_current_node(m))
        feasible_flag = feasible_flag && set_constraint_propagation_fbbt!(m)
        (same_box(ns,_current_node(m),0.0) || !feasible_flag) && break
    end

    if _obbt_depth(m) >= _iteration_count(m)
        for k = 1:_obbt_repetitions(m)
            ns = NodeBB(_current_node(m))
            feasible_flag = feasible_flag && obbt!(m)
            m._obbt_performed_flag = true
            (same_box(ns,_current_node(m),0.0) || !feasible_flag) && break
        end
    end
    m._preprocess_feasibility = feasible_flag
    return
end
preprocess!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = preprocess!(_ext(m), m)

"""
$(TYPEDSIGNATURES)

Returns `true` if a cut should be added and computes a new reference point to add the
cut at. By default, checks that `cut_max_iterations` are not exceeded and that the
improvement in the objective value associated with the previous cut is greater than
both an absolute tolerance `cut_ϵ_abs` and a relative tolerance `cut_ϵ_rel`. Returns
`false` otherwise.
"""
function cut_condition(t::ExtensionType, m::GlobalOptimizer)
    obj_old = m._last_cut_objective
    obj_new = m._lower_objective_value
    flag = m._cut_iterations < _cut_max_iterations(m)
    flag &= obj_new - obj_old > _cut_ϵ_abs(m)
    flag &= obj_new - obj_old > _cut_ϵ_rel(m)*abs(obj_new)
    return flag
end
cut_condition(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = cut_condition(_ext(m), m)

"""
$(SIGNATURES)

Returns `true` that the subproblem at the current node `n` has participating integer variables
that have not been fixed to constant valued as the branch-and-bound algorithm progresses. Returns
`false` otherwise.
"""
is_integer_subproblem(m) = !continuous(_current_node(m))

"""
$(TYPEDSIGNATURES)

Constructs a relaxation of the MINLP on node `y` and solves it using the default EAGO
relaxation scheme. By default, EAGO applies Kelley's algorithm (from Kelley Jr., J.E.:
The cutting-plane method for solving convex programs. J. Soc. Ind. Appl. Math. 8(4),
703 to 712 (1960)) while `cut_condition(m)` returns `true` then activates the integrality
constraints of the relaxed problems and solves the resulting MILP relaxation. results
are stored to the `_lower_solution`, `_lower_termination_status`, `_lower_primal_status`,
`_lower_dual_status`, `_lower_objective_value`, and `_lower_feasibility`. Further, lower
and upper variable duals are stored `_lower_lvd` and `_lower_uvd`, respectively, for use
in duality based bound tightening. If relaxation-based bounds are weaker or cutting-planes
are numerically poorly ill-posed, then interval bounds are used instead. If the problem is
dual feasible but the primal status is ambiguous the dual objective value is used for the
lower bound to avoid numerical issues.
"""
function lower_problem!(t::ExtensionType, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}

    num_feasible_relax_flag = true

    d = _relaxed_optimizer(m)
    m._last_cut_objective = typemin(Float64)
    m._lower_objective_value = typemin(Float64)

    t_status = MOI.OPTIMIZE_NOT_CALLED
    p_status = MOI.OTHER_RESULT_STATUS
    d_status = MOI.OTHER_RESULT_STATUS
    status = RRS_INVALID

    set_first_relax_point!(m)
    MOI.set(d, MOI.ObjectiveFunction{SAF}(), m._working_problem._objective_saf)

    while true
        valid_prob, feas_flag = relax_problem!(m)
        if !feas_flag
            num_feasible_relax_flag = false
            break
        end
        m._last_cut_objective = m._lower_objective_value
        MOI.optimize!(d)
        t_status = MOI.get(d, MOI.TerminationStatus())
        p_status = MOI.get(d, MOI.PrimalStatus())
        d_status = MOI.get(d, MOI.DualStatus())
        status = relaxed_problem_status(t_status, p_status, d_status)
        if status != RRS_OPTIMAL
            break
        end
        m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
        if cut_condition(m)
            store_lower_solution!(m, d)
            m._cut_iterations += 1
        else
            break
        end
    end
    if !num_feasible_relax_flag
        status = RRS_INFEASIBLE
    end

    # activate integrality conditions for MIP & solve MIP subproblem
    if is_integer_subproblem(m) && (status !== RRS_INFEASIBLE)
        m._last_cut_objective = m._lower_objective_value
        for i = 1:_variable_num(BranchVar(), m)
            l = _lower_bound(BranchVar(), m, i)
            u = _upper_bound(BranchVar(), m, i)
            if is_integer(BranchVar(), m, i) && (l != u)
                c_integer = MOI.add_constraint(d, VI(_bvi(m, i)), MOI.Integer())
                push!(m._relaxed_variable_integer, c_integer)
            end
        end
        MOI.optimize!(d)
        t_status = MOI.get(d, MOI.TerminationStatus())
        p_status = MOI.get(d, MOI.PrimalStatus())
        d_status = MOI.get(d, MOI.DualStatus())
        status = relaxed_problem_status(t_status, p_status, d_status)
        if status == RRS_OPTIMAL
            m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
        end
    end

    # check status -- if not feasible/infeasible then fallback to interval bounds
    if status == RRS_OPTIMAL
        m._lower_objective_value = MOI.get(d, MOI.ObjectiveValue())
    end

    m._lower_termination_status = t_status
    m._lower_primal_status = p_status
    m._lower_dual_status = d_status
    status = relaxed_problem_status(t_status, p_status, d_status)
    if !num_feasible_relax_flag
        status = RRS_INFEASIBLE
    end
    if status == RRS_INFEASIBLE
        m._lower_feasibility  = false
        m._lower_objective_value = -Inf
        return
    end

    # set dual values
    set_dual!(m)
    m._lower_feasibility = true
    store_lower_solution!(m, d)
    if status == RRS_DUAL_FEASIBLE
        m._lower_objective_value = MOI.get(d, MOI.DualObjectiveValue())
    end
    interval_objective_bound!(m, true)
    return
end
lower_problem!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = lower_problem!(_ext(m), m)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
workspace()
#using EAGOSemiInfinite
#using EAGOGlobalSolver
using IntervalArithmetic
using EAGOBranchBound
#using EAGOSmoothMcCormickGrad
#using StaticArrays

function FDI_invApprox(x::Interval)
  Ap = 0.00005
  g = 9.8
  low = -(x^2)/(2*g*Ap^2)
  high = (x^2)/(2*g*Ap^2)
  if x.lo>0
    return high
  elseif x.hi<0
    return low
  else
    return Interval(0,max(low.hi,high.hi))
  end
end

function FDI_invApprox(x::Float64)
  Ap = 0.00005
  g = 9.8
  if x>=0
    return (x^2)/(2*g*Ap^2)
  else x<0
    return -(x^2)/(2*g*Ap^2)
  end
end

function xbnds(u,p)

         # parameters
         Ap = 0.00005
         g = 9.8
         z = []
         # bound z via McCormick operators
         z1 = (((u[1]/p[3])/Ap)^2)/(2.0*g)
         z3 = (1.0/(2.0*g))*((u[1]+u[2])/(p[4]*Ap))^2
         z2 = (-((u[2]-p[4]*Ap*sqrt(2*g*z3))/(p[5])/Ap)^2)/(2.0*g)
         z4 = (((u[1]*p[1]/p[3])/Ap)^2)/(2.0*g)
         z6 = (1.0/(2.0*g))*((p[1]*u[1]+u[2])/(p[4]*Ap))^2
         z5 = (-((u[2]-p[4]*Ap*sqrt(2*g*z6))/(p[5])/Ap)^2)/(2.0*g)
         z7 = (((u[1]/p[3])/Ap)^2)/(2.0*g)
         z9 = (1.0/(2.0*g))*((u[1]+u[2])/(p[4]*(Ap+3.14159*p[2]^2)))^2
         z8 = FDI_invApprox((u[2]-p[4]*(Ap+3.14159*p[2]^2)*sqrt(2.0*g*z9))/(p[5]))
         push!(z,z1,z2,z3,z4,z5,z6,z7,z8,z9)

         # rotate coordinates back to X & intersects

         x = []
         x1 = min(max(z1-z2+z3,0.0),0.75)
         x2 = min(max(z3,0.0),0.75)
         x3 = min(max(z3-z2,0.0),0.75)
         x4 = min(max(z4-z5+z6,0.0),0.75)
         x5 = min(max(z6,0.0),0.75)
         x6 = min(max(z6-z5,0.0),0.75)
         x7 = min(max(z7-z8+z9,0.0),0.75)
         x8 = min(max(z9,0.0),0.75)
         x9 = min(max(z9-z8,0.0),0.75)
         push!(x,x1,x2,x3,x4,x5,x6,x7,x8,x9)

         return x
end

# sets up SIP constraint with imbedded x(p) calculation
function SIP_g(un,p)
  u = un[1:2]
  xbnd = xbnds(u,p)
  r_xbnds = [zero(xbnd[1]) for i=1:3,j=1:3]
  r_xbnds[:,1] = xbnd[1:3]
  r_xbnds[:,2] = xbnd[4:6]
  r_xbnds[:,3] = xbnd[7:9]
  out = un[3]
  for i=1:3
    for j=1:2
      for k=(j+1):3
        out = out - (r_xbnds[i,j]-r_xbnds[i,k])^2
      end
    end
  end
  return out
end

# sets up SIP objective
function SIP_f(un)
  return -un[3]
end
# sets up SIP
ex_SIP_X = [(1E-5)..(1E-4),(1E-5)..(1E-4),(-10.0)..(10.0)]
ex_SIP_P = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)]

#Intv = SIP_g(ex_SIP_X,ex_SIP_P)
#Trial1 = [1E-5,1E-5,0.0]
#Intv1 = SIP_g(Trial1,ex_SIP_P)
#Trial2 = [1E-5,1E-4,0.0]
#Intv2 = SIP_g(Trial2,ex_SIP_P)
#Trial3 = [1E-4,1E-5,0.0]
#Intv3 = SIP_g(Trial3,ex_SIP_P)
#Trial4 = [1E-4,1E-4,0.0]
#Intv4 = SIP_g(Trial4,ex_SIP_P)
#Trial5 = [(1E-5+1E-4)/2,(1E-5+1E-4)/2,0.0]
#Intv5 = SIP_g(Trial5,ex_SIP_P)
#Trial6 = [(1E-5+1E-4)/3,(1E-5+1E-4)/3,0.0]
#Intv6 = SIP_g(Trial6,ex_SIP_P)
#Trial7 = [2*(1E-5+1E-4)/3,2*(1E-5+1E-4)/3,0.0]
#Intv7 = SIP_g(Trial7,ex_SIP_P)
#Trial8 = [(1E-5+1E-4)/2,(1E-5+1E-4)/3,0.0]
#Intv8 = SIP_g(Trial8,ex_SIP_P)
#Trial9 = [(1E-5+1E-4)/2,2*(1E-5+1E-4)/3,0.0]
#Intv9 = SIP_g(Trial9,ex_SIP_P)
#Trial10 = [(1E-5+1E-4)/3,(1E-5+1E-4)/2,0.0]
#Intv10 = SIP_g(Trial10,ex_SIP_P)
Trial11 = [2*(1E-5+1E-4)/3,(1E-5+1E-4)/2,0.0]
Intv11 = SIP_g(Trial11,ex_SIP_P)

ex_SIP_P1a = [0.66,(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)]
ex_SIP_P2a = [(0.54)..(0.66),0.005,(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)]
ex_SIP_P3a = [(0.54)..(0.66),(0.0005)..(0.005),1.15,(0.65)..(0.95),(0.85)..(1.15)]
ex_SIP_P4a = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.95,(0.85)..(1.15)]
ex_SIP_P5a = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),1.15]
Intv11a1 = SIP_g(Trial11,ex_SIP_P1a)
Intv11b1 = SIP_g(Trial11,ex_SIP_P2a)
Intv11c1 = SIP_g(Trial11,ex_SIP_P3a)
Intv11d1 = SIP_g(Trial11,ex_SIP_P4a)
Intv11e1 = SIP_g(Trial11,ex_SIP_P5a)

ex_SIP_P1b = [0.54,(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)]
ex_SIP_P2b = [(0.54)..(0.66),0.0005,(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)]
ex_SIP_P3b = [(0.54)..(0.66),(0.0005)..(0.005),0.85,(0.65)..(0.95),(0.85)..(1.15)]
ex_SIP_P4b = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.65,(0.85)..(1.15)]
ex_SIP_P5b = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),0.85]
Intv11a2 = SIP_g(Trial11,ex_SIP_P1b)
Intv11b2 = SIP_g(Trial11,ex_SIP_P2b)
Intv11c2 = SIP_g(Trial11,ex_SIP_P3b)
Intv11d2 = SIP_g(Trial11,ex_SIP_P4b)
Intv11e2 = SIP_g(Trial11,ex_SIP_P5b)
ex_SIP_P4ca = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.70,(0.85)..(1.15)]
ex_SIP_P4c = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.75,(0.85)..(1.15)]
ex_SIP_P4da = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.80,(0.85)..(1.15)]
ex_SIP_P4d = [(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),0.85,(0.85)..(1.15)]
Intv11d2a = SIP_g(Trial11,ex_SIP_P4ca)
Intv11d3 = SIP_g(Trial11,ex_SIP_P4c)
Intv11d3a = SIP_g(Trial11,ex_SIP_P4da)
Intv11d4 = SIP_g(Trial11,ex_SIP_P4d)

disc1 = [0.6,0.0005,1.0,0.65,1.00]
ex_SIP_Xt = [(1E-5)..(1E-4),(1E-5)..(1E-4),(0.0)..(0.0)]
Xtest = SIP_g(ex_SIP_Xt,disc1)
m = BnBModel([(1E-5)..(1E-4),(1E-5)..(1E-4),(-10.0)..(10.0)])
s = BnBSolver()
set_to_default!(s)
s.BnB_atol = 1E-2

# solve optimization problem via interval extension
function test_lower(X,k,pos,opt_inner,UBD)
  FInt = -X[3]
  disc1 = [0.6,0.005,1.0,0.65,1.00]
  GInt = SIP_g(X,disc1)
  feas = true
  for i=1:length(GInt)
    if (GInt[i].lo>0.0)
      feas = false
      break
    end
  end
  val = FInt.lo
  pnt = mid.(X)
  return pnt, val, feas, X, [FInt]
end

function test_upper(X,k,pos,opt_inner,temp)
  feas = true
  val = temp[1].hi
  pnt = mid.(X)
  return pnt, val, feas, X
end

s.Lower_Prob = test_lower
s.Upper_Prob = test_upper
solveBnB!(s,m)

mi = BnBModel([(0.54)..(0.66),(0.0005)..(0.005),(0.85)..(1.15),(0.65)..(0.95),(0.85)..(1.15)])
si = BnBSolver()
set_to_default!(si)
si.BnB_atol = 5E-2

function inner_lower(X,k,pos,opt_inner,UBD)
  disc1 = [5.5E-5,5.5E-5,4.99]
  FInt = -SIP_g(disc1,X)
  feas = true
  val = FInt.lo
  pnt = mid.(X)
  return pnt, val, feas, X, [FInt]
end

function inner_upper(X,k,pos,opt_inner,temp)
  feas = true
  val = temp[1].hi
  pnt = mid.(X)
  return pnt, val, feas, X
end

si.Lower_Prob = inner_lower
si.Upper_Prob = inner_upper
#solveBnB!(si,mi)

m1 = BnBModel([(1E-5)..(1E-4),(1E-5)..(1E-4),(-10.0)..(10.0)])
s1 = BnBSolver()
set_to_default!(s1)
s1.BnB_atol = 1E-2

# solve optimization problem via interval extension
function test_lower(X,k,pos,opt_inner,UBD)
  FInt = -X[3]
  disc1 = [0.6,0.005,1.0,0.65,1.00]
  GInt = SIP_g(X,disc1) +0.5
  feas = true
  for i=1:length(GInt)
    if (GInt[i].lo>0.0)
      feas = false
      break
    end
  end
  val = FInt.lo
  pnt = mid.(X)
  return pnt, val, feas, X, [FInt]
end

function test_upper(X,k,pos,opt_inner,temp)
  feas = true
  val = temp[1].hi
  pnt = mid.(X)
  return pnt, val, feas, X
end

s1.Lower_Prob = test_lower
s1.Upper_Prob = test_upper
solveBnB!(s1,m1)

#=
SIPopt = SIP_opts()
set_to_default!(SIPopt)
SIPopt.Verbosity = "Full" #"Normal"

# solves sample LBP
vals = 1
count = 1
p_sto = zeros(vals^5,5)
for i=1:vals
  for j=1:vals
    for k=1:vals
      for p=1:vals
        for q=1:vals
          t1 = ex_SIP_P[1].lo + (ex_SIP_P[1].hi - ex_SIP_P[1].lo)*(i-1)/(vals-1)
          t2 = ex_SIP_P[2].lo + (ex_SIP_P[2].hi - ex_SIP_P[2].lo)*(j-1)/(vals-1)
          t3 = ex_SIP_P[3].lo + (ex_SIP_P[3].hi - ex_SIP_P[3].lo)*(k-1)/(vals-1)
          t4 = ex_SIP_P[4].lo + (ex_SIP_P[4].hi - ex_SIP_P[4].lo)*(p-1)/(vals-1)
          t5 = ex_SIP_P[5].lo + (ex_SIP_P[5].hi - ex_SIP_P[5].lo)*(q-1)/(vals-1)
          p_sto[count,:] = [t1,t2,t3,t4,t5]
          count += 1
        end
      end
    end
  end
end

bnd_z2 = ex_SIP_X[2]-ex_SIP_P[4]*0.00005*sqrt((ex_SIP_P[1]*ex_SIP_X[1]+ex_SIP_X[2])/(ex_SIP_P[5]*(ex_SIP_P[4]*0.00005)^2))
bnd_z5 = ex_SIP_X[2]-ex_SIP_P[4]*0.00005*sqrt((ex_SIP_P[1]*ex_SIP_X[1]+ex_SIP_X[2])/(ex_SIP_P[5]*(ex_SIP_P[4]*0.00005)^2))
bnd_z8 = (ex_SIP_X[2]-ex_SIP_P[4]*(0.00005+3.14159*ex_SIP_P[2]^2)*sqrt((((ex_SIP_X[1]+ex_SIP_X[2])/(ex_SIP_P[4]*(0.00005+3.14159*ex_SIP_P[2]^2)))^2)))/(ex_SIP_P[5])

SIPopt.tol = 1E-1
SIPopt.eps_g0 = 1.0
SIPopt.P_LBD = [p_sto[i,:] for i=1:vals^5]
SIPopt.P_UBD = [p_sto[i,:] for i=1:vals^5]
SIPopt.LBP_Opt.DAG_depth = -1
SIPopt.LBP_Opt.f = SIP_f
SIPopt.LBP_Opt.g = SIP_g
SIPopt.LBP_Opt.X = ex_SIP_X
SIPopt.LBP_Opt.g = x -> BndProb_reform(x,[],SIP_g,[],0.0)
#SIPopt.LBP_Opt.LBD_func_relax = "Diff2-MV-OFF"
#SIPopt.LBP_Opt.LBD_problem_relax = "NLP2"
#SIPopt.LBP_Opt.LBD_problem_solver = "Ipopt"
SIPopt.LBP_Opt.LBD_func_relax = "Interval"
SIPopt.LBP_Opt.LBD_problem_relax = "Interval"
SIPopt.LBP_Opt.LBD_problem_solver = "Interval"
SIPopt.LBP_Opt.UBD_func_relax = "Interval"
SIPopt.LBP_Opt.UBD_problem_relax = "Interval"
SIPopt.LBP_Opt.UBD_problem_solver = "Interval"
SIPopt.LBP_Opt.BnB_options.Verbosity = "Normal" #"Full"
SIPopt.LBP_Opt.BnB_options.BnB_tol = 1E-2
#LBDg, xbar, feas, tLBP1,tLBP2 = EAGO_Global_Explicit(SIPopt.LBP_Opt)
#time_lower = @time EAGO_Global_Explicit(SIPopt.LBP_Opt)
#time_lower = @time EAGO_Global_Explicit(SIPopt.LBP_Opt)

xbar = mid.(ex_SIP_X)
SIPopt.LLP_Opt.DAG_depth = -1
SIPopt.LLP_Opt.f = p -> -SIP_g(xbar,p)
SIPopt.LLP_Opt.X = ex_SIP_P
#SIPopt.LLP_Opt.LBD_func_relax = "Diff2-MV-OFF"
#SIPopt.LLP_Opt.LBD_problem_relax = "NLP2"
#SIPopt.LLP_Opt.LBD_problem_solver = "Ipopt"
SIPopt.LLP_Opt.LBD_func_relax = "Interval"
SIPopt.LLP_Opt.LBD_problem_relax = "Interval"
SIPopt.LLP_Opt.LBD_problem_solver = "Interval"
SIPopt.LLP_Opt.UBD_func_relax = "Interval"
SIPopt.LLP_Opt.UBD_problem_relax = "Interval"
SIPopt.LLP_Opt.UBD_problem_solver = "Interval"
SIPopt.LLP_Opt.BnB_options.Verbosity = "Normal"
SIPopt.LLP_Opt.BnB_options.BnB_tol = 1E-2
#INNg1, pbar,feas,tLLP1,tLLP2 = EAGO_Global_Explicit(SIPopt.LLP_Opt)
#time_inner = @time EAGO_Global_Explicit(SIPopt.LLP_Opt)
#time_inner = @time EAGO_Global_Explicit(SIPopt.LLP_Opt)

SIPopt.UBP_Opt.DAG_depth = -1
SIPopt.UBP_Opt.f = SIP_f
SIPopt.UBP_Opt.g = SIP_g
SIPopt.UBP_Opt.X = ex_SIP_X
SIPopt.UBP_Opt.g = x -> BndProb_reform(x,[],SIP_g,[],1.0)
#SIPopt.UBP_Opt.LBD_func_relax = "Diff2-MV-OFF"
#SIPopt.UBP_Opt.LBD_problem_relax = "NLP2"
#SIPopt.UBP_Opt.LBD_problem_solver = "Ipopt"
SIPopt.UBP_Opt.LBD_func_relax = "Interval"
SIPopt.UBP_Opt.LBD_problem_relax = "Interval"
SIPopt.UBP_Opt.LBD_problem_solver = "Interval"
SIPopt.UBP_Opt.UBD_func_relax = "Interval"
SIPopt.UBP_Opt.UBD_problem_relax = "Interval"
SIPopt.UBP_Opt.UBD_problem_solver = "Interval"
SIPopt.UBP_Opt.BnB_options.Verbosity = "Normal"
SIPopt.UBP_Opt.BnB_options.BnB_tol = 1E-2
#UBD_temp, xbar, feas,tUBP1,tUBP2 = EAGO_Global_Explicit(SIPopt.UBP_Opt)
#time_upper = @time EAGO_Global_Explicit(SIPopt.UBP_Opt)
#time_upper = @time EAGO_Global_Explicit(SIPopt.UBP_Opt)
#pbar = [0.66, 0.00275, 1.13595, 0.95, 1.15]
#ubar = [5.5e-05, 0.000100001] × [5.5e-05, 0.000100001] × [9.9781, 9.97815]

xbnd_range1 = xbnds(ex_SIP_X[1:2],ex_SIP_P)
SIP_range = SIP_g(ex_SIP_X,ex_SIP_P)
ex_out = Explicit_SIP_Solve(SIP_f,[],SIP_g,ex_SIP_X,ex_SIP_P,SIPopt)
=#


text:
module TestConstraintZeroOne

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_ZeroOne()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.ZeroOne{Float64}(mock)

    bridge_type = MOI.Bridges.Constraint.ZeroOneBridge{Float64}
    @test MOI.supports_constraint(bridge_type, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.Bridges.Constraint.concrete_bridge_type(
        bridge_type,
        MOI.VariableIndex,
        MOI.ZeroOne,
    ) == bridge_type
    MOI.Test.test_basic_VariableIndex_ZeroOne(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [1, 0, 0, 1, 1]),
    )
    MOI.Test.test_linear_integer_knapsack(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(),
        ),
    )
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), ci) == 1
    _test_delete_bridge(
        bridged_mock,
        ci,
        5,
        (
            (MOI.VariableIndex, MOI.Integer, 0),
            (MOI.VariableIndex, MOI.Interval{Float64}, 0),
        ),
        num_bridged = 5,
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 20.0)
            MOI.Utilities.mock_optimize!(mock, [4, 5, 1])
        end,
    )
    MOI.Test.test_linear_integer_integration(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [1.0; zeros(10)]),
    )
    MOI.Test.test_linear_integer_solve_twice(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(),
        ),
    )
    s = """
    variables: x, y
    y == 1.0
    x in ZeroOne()
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    sb = """
    variables: x, y
    y == 1.0
    x in Integer()
    x in Interval(0.0,1.0)
    minobjective: x
    """
    modelb = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(modelb, sb)
    MOI.empty!(bridged_mock)
    @test MOI.is_empty(bridged_mock)
    MOI.Utilities.loadfromstring!(bridged_mock, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        ["x", "y"],
        String[],
        [("y", MOI.EqualTo{Float64}(1.0)), ("x", MOI.ZeroOne())],
    )
    MOI.Test.util_test_models_equal(
        mock,
        modelb,
        ["x", "y"],
        String[],
        [
            ("y", MOI.EqualTo{Float64}(1.0)),
            ("x", MOI.Integer()),
            ("x", MOI.Interval{Float64}(0.0, 1.0)),
        ],
    )
    return
end

end  # module

TestConstraintZeroOne.runtests()


text:
using MusicModels
using NoteSequences
using Test

include("musicencoders.jl")


text:
"""
    _calculate_interval_inner_counts(order::Dict{Int,String},
                                          intervals::Dict{String,<:Dates.TimePeriod},
                                          step_resolution::Dates.TimePeriod)

Calculates how many times a stage is executed for every interval of the previous stage
"""
function _calculate_interval_inner_counts(
    order::Dict{Int, String},
    intervals::Dict{String, Tuple{Dates.TimePeriod, <:FeedForwardChronology}},
    step_resolution::Dates.TimePeriod,
)
    reverse_order = sort(collect(keys(order)), rev = true)
    interval_run_counts = Dict{Int, Int}()
    for k in reverse_order[1:(end - 1)]
        stage_name = order[k]
        previous_stage_name = order[k - 1]
        stage_interval = intervals[stage_name][1]
        previous_stage_interval = intervals[previous_stage_name][1]
        if Dates.Millisecond(previous_stage_interval % stage_interval) !=
           Dates.Millisecond(0)
            throw(IS.ConflictingInputsError("The interval configuration provided results in a fractional number of executions of stage $stage_name"))
        end
        interval_run_counts[k] = previous_stage_interval / stage_interval
        @debug "Stage $k is executed $(interval_run_counts[k]) time within each interval of Stage $(k-1)"
    end
    stage_name = order[1]
    stage_interval = intervals[stage_name][1]
    interval_run_counts[1] = 1
    return interval_run_counts
end

""" Function calculates the total number of stage executions in the simulation and allocates the appropiate vector"""
function _allocate_execution_order(interval_run_counts::Dict{Int, Int})
    total_size_of_vector = 0
    for (k, counts) in interval_run_counts
        mult = 1
        for i in 1:k
            mult *= interval_run_counts[i]
        end
        total_size_of_vector += mult
    end
    return -1 * ones(Int, total_size_of_vector)
end

function _fill_execution_order(
    execution_order::Vector{Int},
    interval_run_counts::Dict{Int, Int},
)
    function _fill_stage(index::Int, stage::Int)
        if stage < last_stage
            next_stage = stage + 1
            for i in 1:interval_run_counts[next_stage]
                index = _fill_stage(index, next_stage)
            end
        end
        execution_order[index] = stage
        index -= 1
    end

    index = length(execution_order)
    stages = sort!(collect(keys(interval_run_counts)))
    last_stage = stages[end]
    _fill_stage(index, stages[1])
    return
end

function _get_execution_order_vector(
    order::Dict{Int, String},
    intervals::Dict{String, Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}},
    step_resolution::Dates.TimePeriod,
)
    length(order) == 1 && return [1]
    interval_run_counts =
        _calculate_interval_inner_counts(order, intervals, step_resolution)
    execution_order_vector = _allocate_execution_order(interval_run_counts)
    _fill_execution_order(execution_order_vector, interval_run_counts)
    @assert isempty(findall(x -> x == -1, execution_order_vector))
    return execution_order_vector
end

function _check_stage_order(order::Dict{Int, String})
    length(order) == 1 && return
    if !mapreduce(x -> x == 1, *, diff(sort!(collect(keys(order)))))
        throw(IS.InvalidValue("Keys in the order dictionary aren't specified as consecutive integrers 1 -> N"))
    end
    return
end

function _check_all_inputs_present(
    order::Dict{Int, String},
    intervals::Dict{String, <:Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}},
    horizons::Dict{String, Int},
)
    for stage_name in values(order)
        if !(stage_name in keys(horizons))
            throw(IS.ConflictingInputsError("Horizon not defined for stage $(stage_name)"))
        end
        if !(stage_name in keys(intervals))
            throw(IS.ConflictingInputsError("Interval not defined for stage $(stage_name)"))
        end
    end
    return
end

function _check_feedforward(
    feedforward::Dict{Tuple{String, Symbol, Symbol}, <:AbstractAffectFeedForward},
    feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology},
)
    isempty(feedforward) && return
    for stage_key in keys(feedforward)
        @debug stage_key
        if !mapreduce(x -> x.second == stage_key[1], |, keys(feedforward_chronologies))
            throw(ArgumentError("No valid Chronology has been defined for the feedforward added to $(stage_key[1])"))
        end
    end
    return
end

function _check_chronology_consistency(
    order::Dict{Int, String},
    feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology},
    ini_cond_chronology::InitialConditionChronology,
)
    if isempty(feedforward_chronologies)
        @warn("No Feedforward Chronologies have been defined. This configuration assummes that there is no information passing between stages")
    end
    if length(order) == 1
        if isa(ini_cond_chronology, InterStageChronology)
            @warn("Single stage detected, the default Initial Condition Chronology is IntraStageChronology(), other values will be ignored.")
        end
    end
    #TODO: Add more consistency checks
    return
end

function _check_step_interval_consistency(
    step_resolution::Dates.TimePeriod,
    first_stage_interval::Dates.TimePeriod,
)
    if step_resolution != first_stage_interval
        throw(IS.ConflictingInputsError("step_resolution != first_stage_interval : this specification is not currently supported, adjust step_resolution or intervals arguments."))
    end
    return
end

function _check_cache_defination(cache::Dict{<:Tuple, <:AbstractCache})
    for (stage_names, c) in cache
        if typeof(c) == TimeStatusChange && length(stage_names) > 1
            error("TimeStatusChange cache currently only supports single stage. Please consider changing your cache definations")
        end
    end
    return
end

# TODO: Add DocString
@doc raw"""
    SimulationSequence(horizons::Dict{String, Int}
                        step_resolution::Dates.TimePeriod
                        intervals::Dict{String, <:Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}}
                        order::Dict{Int, String}
                        feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology}
                        feedforward::Dict{Tuple{String, Symbol, Symbol}, <:AbstractAffectFeedForward}
                        ini_cond_chronology::Dict{String, <:FeedForwardChronology}
                        cache::Dict{String, AbstractCache}
                        )
"""
mutable struct SimulationSequence
    horizons::Dict{String, Int}
    step_resolution::Dates.TimePeriod
    # The string here is the name of the stage
    intervals::Dict{String, Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}}
    order::Dict{Int, String}
    feedforward_chronologies::Dict{Pair{String, String}, <:FeedForwardChronology}
    feedforward::Dict{Tuple{String, Symbol, Symbol}, <:AbstractAffectFeedForward}
    ini_cond_chronology::InitialConditionChronology
    cache::Dict{Tuple, AbstractCache}
    execution_order::Vector{Int}
    current_execution_index::Int

    function SimulationSequence(;
        horizons::Dict{String, Int},
        step_resolution::Dates.TimePeriod,
        intervals::Dict{String, <:Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}},
        order::Dict{Int, String},
        feedforward_chronologies = Dict{Pair{String, String}, FeedForwardChronology}(),
        feedforward = Dict{Tuple{String, Symbol, Symbol}, AbstractAffectFeedForward}(),
        ini_cond_chronology = InterStageChronology(),
        cache = Dict{Tuple, AbstractCache}(),
    )
        _check_stage_order(order)
        _check_all_inputs_present(order, intervals, horizons)
        _intervals = Dict{String, Tuple{<:Dates.TimePeriod, <:FeedForwardChronology}}()
        for (k, v) in intervals
            _intervals[k] = (IS.time_period_conversion(intervals[k][1]), intervals[k][2])
        end
        step_resolution = IS.time_period_conversion(step_resolution)
        _check_step_interval_consistency(step_resolution, intervals[order[1]][1])
        _check_feedforward(feedforward, feedforward_chronologies)
        _check_chronology_consistency(order, feedforward_chronologies, ini_cond_chronology)
        _check_cache_defination(cache)
        if length(order) == 1
            ini_cond_chronology = IntraStageChronology()
        end
        new(
            horizons,
            step_resolution,
            _intervals,
            order,
            feedforward_chronologies,
            feedforward,
            ini_cond_chronology,
            cache,
            _get_execution_order_vector(order, _intervals, step_resolution),
            0,
        )
    end
end

function get_stage_horizon(s::SimulationSequence, stage::String)
    horizon = get(s.horizons, stage, nothing)
    isnothing(horizon) &&
        throw(ArgumentError("Stage $(stage.internal.number) not present in the simulation"))
    return horizon
end

get_stage_interval(s::SimulationSequence, stage::String) = s.intervals[stage][1]

function get_stage_name(s::SimulationSequence, stage::Stage)
    name = get(s.order, get_number(stage), nothing)
    isnothing(name) &&
        throw(ArgumentError("Stage $(stage.internal.number) not present in the simulation"))
    return name
end

get_step_resolution(s::SimulationSequence) = s.step_resolution

function get_stage_interval_chronology(s::SimulationSequence, stage::String)
    return s.intervals[stage][2]
end


text: hist(randn(100); col="powderblue")


text:
#### DEFINE A PEDESTRIAN MODEL AND A PEDESTRIAN FLOW #####

const PED_WIDTH = 1.0
const PED_LENGTH = 1.0
const PEDESTRIAN_DEF = VehicleDef(AgentClass.PEDESTRIAN, PED_LENGTH, PED_WIDTH)

"""
Action type for the pedestrian model
"""
struct ConstantSpeedDawdling
    v::Float64 # instantaneous change in velocity
    lat::Float64 # lateral displacement
end

"""
    Simple model, pedestrian walking at constant speed with noise
    and dawdling in the y direction
"""
@with_kw mutable struct ConstantPedestrian <: DriverModel{ConstantSpeedDawdling}
    update_freq::Int64 = 1
    tick::Int64 = 0.
    a::ConstantSpeedDawdling = ConstantSpeedDawdling(0., 0.)
    dawdling_amp::Float64 = 0.
    v_noise::Float64 = 0.5
    v_desired::Float64 = 1.0
    dt::Float64 = 0.1
end

AutomotiveDrivingModels.get_name(model::ConstantPedestrian) = "ConstantPedestrian"
Base.rand(rng::AbstractRNG, model::ConstantPedestrian) = model.a
Distributions.pdf(model::ConstantPedestrian, a::ConstantSpeedDawdling) = 1/3*pdf(Normal(0., model.dawdling_amp))
AutomotiveDrivingModels.reset_hidden_state!(model::ConstantPedestrian) = model

function AutomotiveDrivingModels.observe!(model::ConstantPedestrian,
                                          scene::EntityFrame{VehicleState, VehicleDef, Int64},
                                          roadway::Roadway,
                                          egoid::Int)
    if model.tick % model.update_freq == 0
        model.a = ConstantSpeedDawdling(model.v_desired + rand([-model.v_noise, 0, model.v_noise]),
                                        model.dawdling_amp*randn())
    end
    model.a = ConstantSpeedDawdling(model.a.v, model.dawdling_amp*randn())
    model.tick += 1
end

function AutomotiveDrivingModels.propagate(veh::VehicleState, action::ConstantSpeedDawdling, roadway::Roadway, ΔT::Float64)
    lane_tag_orig = veh.posF.roadind.tag

    # update value in Frenet coordinates
    v = veh.v
    lat = action.lat
    ϕ = veh.posF.ϕ
    s = veh.posF.s
    Δs = v*cos(ϕ)*ΔT + lat*sin(ϕ)
    t = veh.posF.t
    Δt = v*sin(ϕ)*ΔT + lat*cos(ϕ)
    lane = get_lane(roadway, veh)
    # cannot go offroad
    if t + Δt > lane.width/2
        Δt = 0.
    elseif t - Δt < -lane.width/2
        Δt = 0.
    end

    speed = action.v
    # println(s + Δs)
    s_ = clamp(s + Δs, 0., get_end(lane))
    # println(s_)

    #XXX what exactly is happening here, ask Tim ?
    # roadind = move_along(veh.posF.roadind, roadway, Δs)
    # footpoint = roadway[roadind]
    # posG = convert(VecE2, footpoint.pos) + polar(t + Δt, footpoint.pos.θ + π/2)

    # posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ)

    # state = VehicleState(posG, roadway, action.v)

    # roadproj = proj(state.posG, roadway[lane_tag_orig], roadway, move_along_curves=false)
    # retval = VehicleState(Frenet(roadproj, roadway), roadway, state.v)
    # return retval
    return VehicleState(Frenet(lane, s_, t + Δt, ϕ), roadway, speed)
end

function AutomotiveDrivingModels.propagate(veh::Vehicle, action::ConstantSpeedDawdling, roadway::Roadway, ΔT::Float64)
    return propagate(veh.state, action, roadway, ΔT)
end

## Discrete distribution
# P(a | s)

function action_space(model::ConstantPedestrian)
    return (ConstantSpeedDawdling(0., 0.), ConstantSpeedDawdling(1., 0.), ConstantSpeedDawdling(2., 0.))
end

function get_distribution(model::ConstantPedestrian)
    actions = action_space(model)
    probs = (1/3, 1/3, 1/3)
    return actions, probs
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Catlab.Present
using Catlab.Graphs
using JSON
using AutoHashEquals
using DataStructures: DefaultDict

using Catlab.Graphs.BasicGraphs: TheoryGraph
using Catlab.Theories
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.CSetDataStructures: struct_acset

import Base: isempty

include(joinpath(@__DIR__, "../../CSetAutomorphisms.jl/src/CSetAutomorphisms.jl"))

"""Edges and vertices labeled by symbols"""
@present TheoryLabeledGraph <: TheoryGraph begin
  Label::AttrType
  vlabel::Attr(V,Label)
  elabel::Attr(E,Label)
end;

@acset_type LabeledGraph_(TheoryLabeledGraph, index=[:src,:tgt])
const LabeledGraph = LabeledGraph_{Symbol}

"""Forget about the labels"""
function to_graph(lg::LabeledGraph_)::Graph
  G = Graph(nparts(lg, :V))
  add_edges!(G, lg[:src], lg[:tgt])
  return G
end

@auto_hash_equals struct Cone
  d::LabeledGraph
  apex::Symbol
  legs::Vector{Pair{Int, Symbol}}
  function Cone(d::LabeledGraph, apex::Symbol, legs::Vector{Pair{Int, Symbol}})
    l1, _ = zip(legs...) # l2 might have duplicates, e.g. monomorphism cone
    length(Set(l1)) == length(legs) || error("nonunique legs $legs")
    return new(d, apex, legs)
  end
end

@auto_hash_equals struct FLS
  name::Symbol
  schema::LabeledGraph
  cones::Vector{Cone}
  eqs::Vector{Tuple{Symbol, Vector{Symbol}, Vector{Symbol}}}
  cset::Type
  crel::Type
  function FLS(name::Symbol, schema::LabeledGraph, cones::Vector{Cone},
               eqs::Vector)
    # Check eqs are well-formed
    for (_, p, q) in eqs
      pe, qe = [[only(incident(schema, edge, :elabel)) for edge in x] for x in [p,q]]
      isempty(qe) || schema[:src][pe[1]] == schema[:src][qe[1]] || error("path eq don't share start point $p \n$q")
      isempty(qe) || schema[:src][pe[end]] == schema[:src][qe[end]] || error("path eq don't share end point $p \n$q")
      !isempty(qe) || schema[:src][pe[1]] == schema[:tgt][pe[end]] || error("path eq has self loop but p doesn't have same start/end $p \n$q")
      all([schema[:tgt][p1]==schema[:src][p2] for (p1, p2) in zip(pe, pe[2:end])]) || error("head/tail mismatch in p $p \n$q")
      all([schema[:tgt][q1]==schema[:src][q2] for (q1, q2) in zip(qe, qe[2:end])]) || error("head/tail mismatch in q $p \n$q")
    end
    # Check cones are well-formed
    for c in cones
      vert = only(incident(schema, c.apex, :vlabel))
      for (v, l) in c.legs
        edge = only(incident(schema, l, :elabel))
        schema[:src][edge] == vert || error("Leg does not come from apex $c")
        schema[:vlabel][schema[:tgt][edge]] == c.d[:vlabel][v] || error("Leg $l -> $v does not go to correct vertex $c")
        is_homomorphic(c.d, schema) || error("Cone diagram does not map into schema $c")
      end
    end
    return new(name, schema, cones, eqs, grph_to_cset(name, schema),
                grph_to_crel(name, schema))
  end
end

@auto_hash_equals struct ChaseStepData
  cones::DefaultDict{Symbol, Vector{Int}}
  tgds::DefaultDict{Symbol, Vector{Pair{Int, Int}}}
  path_eqs::DefaultDict{Symbol, Vector{Pair{Int, Int}}}
  fun_eqs::DefaultDict{Symbol, Vector{Vector{Int}}}
  cone_eqs::DefaultDict{Symbol, Vector{Vector{Int}}}
  function ChaseStepData()
    return new(DefaultDict{Symbol, Vector{Int}}(Vector{Int}),
               DefaultDict{Symbol, Vector{Pair{Int,Int}}}(Vector{Pair{Int,Int}}),
               DefaultDict{Symbol, Vector{Pair{Int,Int}}}(Vector{Pair{Int,Int}}),
               DefaultDict{Symbol, Vector{Vector{Int}}}(Vector{Vector{Int}}),
               DefaultDict{Symbol, Vector{Vector{Int}}}(Vector{Vector{Int}}))
  end
end

isempty(c::ChaseStepData) = all(map(isempty, [
  c.cones, c.tgds, c.path_eqs, c.fun_eqs, c.cone_eqs]))

csd_to_dict(csd::ChaseStepData) = Dict([
  "cones"=>Dict([string(k)=>v for (k,v) in collect(csd.cones)]),
  "tgds" =>Dict([string(k)=>v for (k,v) in csd.tgds]),
  "path_eqs"=>Dict([string(k)=>v for (k,v) in csd.path_eqs]),
  "fun_eqs"=>Dict([string(k)=>collect(v) for (k,v) in csd.fun_eqs]),
  "cone_eqs"=>Dict([string(k)=>collect(v) for (k,v) in csd.cone_eqs])])

src(F::FLS, e::Symbol) = F.schema[:vlabel][F.schema[:src][only(incident(F.schema, e, :elabel))]]
tgt(F::FLS, e::Symbol) = F.schema[:vlabel][F.schema[:tgt][only(incident(F.schema, e, :elabel))]]

cone_to_dict(c::Cone) = Dict([
  "d"=>generate_json_acset(c.d),
  "apex"=>string(c.apex),"legs"=>c.legs])
dict_to_cone(d::Dict)::Cone = Cone(
  parse_json_acset(LabeledGraph,d["d"]), Symbol(d["apex"]),
  [parse(Int, k)=>Symbol(v) for (k, v) in map(only, d["legs"])])
"""TO DO: add cone and eq info to the hash...prob requires CSet for FLS"""
Base.hash(F::FLS) = canonical_hash(to_graph(F.schema))
to_json(F::FLS) = JSON.json(Dict([
  :name=>F.name, :schema=>generate_json_acset(F.schema),
  :cones => [cone_to_dict(c) for c in F.cones],
  :eqs => [Dict([:name=>n,:p=>p,:q=>q]) for (n,p,q) in F.eqs]]))
function fls_from_json(s::String)::FLS
  p = JSON.parse(s)
  return FLS(Symbol(p["name"]), parse_json_acset(LabeledGraph, p["schema"]),
             [dict_to_cone(d) for d in p["cones"]],
             [(Symbol(pq["name"]), map(Symbol, pq["p"]),
               map(Symbol, pq["q"])) for pq in p["eqs"]])
end

add_srctgt(x::Symbol) = Symbol("src_$(x)") => Symbol("tgt_$(x)")

"""This should only be called inside constructing FLS ... move inside?"""
function grph_to_crel(name::Symbol,fls::LabeledGraph)::Type
  name_ = Symbol("rel_$name")
  pres = Presentation(FreeSchema)
  nv = length(fls[:vlabel])
  alledge = vcat([add_srctgt(e) for e in fls[:elabel]]...)
  xobs = [Ob(FreeSchema, s) for s in [fls[:vlabel]...,fls[:elabel]...]]
  for x in xobs
    add_generator!(pres, x)
  end
  for (i,(e, src, tgt)) in enumerate(zip(fls[:elabel],fls[:src], fls[:tgt]))
    s, t = add_srctgt(e)
    add_generator!(pres, Hom(s, xobs[nv+i], xobs[src]))
    add_generator!(pres, Hom(t, xobs[nv+i], xobs[tgt]))
  end
  expr = struct_acset(name_, StructACSet, pres, index=alledge)
  eval(expr)
  return eval(name_)
end

"""This should only be called inside constructing FLS ... move inside?"""
function grph_to_cset(name::Symbol, fls::LabeledGraph)::Type
  pres = Presentation(FreeSchema)
  xobs = [Ob(FreeSchema, s) for s in fls[:vlabel]]
  for x in xobs
    add_generator!(pres, x)
  end
  for (e, src, tgt) in zip(fls[:elabel], fls[:src], fls[:tgt])
    add_generator!(pres, Hom(e, xobs[src], xobs[tgt]))
  end
  expr = struct_acset(name, StructACSet, pres, index=fls[:elabel])
  eval(expr)
  return eval(name)
end

"""
Convert a functional C-Rel to a C-Set. Elements that are
not mapped by a relation are given a target value of 0. If this happens at all,
an output bool will be false
"""
function crel_to_cset(F::FLS, J::StructACSet)::Pair{StructACSet, Bool}
  res = F.cset() # grph_to_cset(F.name, F.schema)
  for o in F.schema[:vlabel]
    add_parts!(res, o, nparts(J, o))
  end
  partial = false
  for m in F.schema[:elabel]
    msrc, mtgt = add_srctgt(m)
    length(J[msrc]) == length(Set(J[msrc])) || error("nonfunctional $J")
    partial |= length(J[msrc]) != nparts(J, src(F, m))
    for (domval, codomval) in zip(J[msrc], J[mtgt])
      set_subpart!(res, domval, m, codomval)
    end
  end
  return res => partial
end

"""
Initialize a relation on a schema from either a model or a dict of cardinalities
for each object.
"""
function initrel(F::FLS,
                 I::Union{Nothing, Dict{Symbol, Int}, StructACSet}=nothing,
                 )::StructACSet
  if !(I isa StructACSet)
    dic = deepcopy(I)
    I = F.cset() # grph_to_cset(F.name, F.schema)
    for (k, v) in (dic === nothing ? [] : collect(dic))
      add_parts!(I, k, v)
    end
  end
  J = F.crel() # grph_to_crel(F.name, F.schema)
  # Initialize data in J from I
  for o in F.schema[:vlabel]
    add_parts!(J, o, nparts(I, o))
  end
  for d in F.schema[:elabel]
    hs, ht = add_srctgt(d)
    for (i, v) in filter(x->x[2]!=0, collect(enumerate(I[d])))
      n = add_part!(J, d)
    set_subpart!(J, n, hs, i)
    set_subpart!(J, n, ht, v)
    end
  end
  return J
end

function realobs(F::FLS)::Set{Symbol}
  return setdiff(Set(F.schema[:vlabel]), [c.apex for c in F.cones])
end
function relsize(F::FLS, I::StructACSet)::Int
  return sum([nparts(I, x) for x in realobs(F)])
end

function get_eq(F::FLS,name::Symbol)::Pair{Vector{Symbol}, Vector{Symbol}}
  return only([p=>q for (n,p,q) in F.eqs if n==name])
end

"""
Query that returns all instances of the base pattern. External variables
are labeled by the legs of the cone.
"""
function cone_query(c::Cone; verbose::Bool=false)::StructACSet
  vars = [Symbol("x$i") for i in nparts(c.d, :V)]
  typs = ["$x(_id=x$i)" for (i, x) in enumerate(c.d[:vlabel])]
  bodstr = vcat(["begin"], typs)
  for (e, s, t) in zip(c.d[:elabel], c.d[:src], c.d[:tgt])
    push!(bodstr, "$e(src_$e=x$s, tgt_$e=x$t)")
  end
  push!(bodstr, "end")
  exstr = "($(join(["$(v)_$i=x$k" for vs in values(vars) for (i, (k,v)) in enumerate(c.legs)],",") ))"
  ctxstr = "($(join(vcat(["x$i::$x" for (i, x) in enumerate(c.d[:vlabel])],),",")))"
  ex  = Meta.parse(exstr)
  ctx = Meta.parse(ctxstr)
  hed = Expr(:where, ex, ctx)
  bod = Meta.parse(join(bodstr, "\n"))
  if verbose
    println("ex $exstr\n ctx $ctxstr\n bod $(join(bodstr, "\n"))")
  end
  res = parse_relation_diagram(hed, bod)
  return res
end


text:
# # Solving the Huggett model

using Test #jl
using Aiyagari
using Optim, QuantEcon, Parameters

PKG_HOME = joinpath(dirname(pathof(Aiyagari)), "..")

# ## Setting up the state space

a_grid = LinRange(-2, 5, 20)

endo = EndogenousStateSpace((a=a_grid,))

z_grid = [0.5; 1.0; 1.5]
z_prob = [0.7 0.15 0.15;
          0.2 0.6 0.2;
          0.15 0.15 0.7]
z_MC = MarkovChain(z_prob, z_grid, :z)

exo = ExogenousStateSpace([z_MC]);

# ## Define the household's problem in the Huggett model

u(c) = c > 0 ? log(c) : 10000 * c - 10000 * one(c)

function get_c(choices, states, aggregate_state)
  @unpack a_next = choices
  @unpack a, z = states
  @unpack r = aggregate_state

  c = a + z - a_next/(1+r)
end

function objective(choices, states, aggregate_state, 𝔼V, params)
  @unpack a_next = choices
  @unpack β = params

  c = get_c(choices, states, aggregate_state)

  u(c) + β * 𝔼V(a_next)
end

function Aiyagari.get_optimum(states, agg_state, 𝔼V, params, endo, hh::Consumer)
  a_grid = endo.grids.a

  a_min, a_max = extrema(a_grid)

  res = Optim.optimize(a_next -> - objective((a_next=a_next,), states, agg_state, 𝔼V, params), a_min, a_max)

  a_next = Optim.minimizer(res)
  val    = - Optim.minimum(res)
  conv   = Optim.converged(res)

  pol = a_next
  pol_full = (a_next=a_next, c=get_c((a_next=a_next, ), states, agg_state))

  (pol=pol, pol_full=pol_full, val=val, conv=conv)
end

mutable struct HuggettAS{T1,T2} <: AggregateState
  r::T1
  dist::T2 # the distribution over idiosynchratic states
end

function HuggettAS(r, a_grid, z_MC)
  dist_proto = zeros((length(a_grid), length(z_MC.state_values)))
  HuggettAS(r, dist_proto)
end;

# ## Initialize the aggregate state and solve the model

agg_state = HuggettAS(0.05, a_grid, z_MC)
param = (β = 0.9, )

#using BenchmarkTools
@unpack val, policy, policies_full = solve_bellman(endo, exo, agg_state, param, Consumer(), rtol=√eps())
#22 ms 176 itr

using DelimitedFiles #jl
using DelimitedFiles #src
#writedlm(joinpath(PKG_HOME, "test/matrices", "huggett_value.txt"), val) #src
value_test = readdlm(joinpath(PKG_HOME, "test/matrices", "huggett_value.txt")) #jl

@test all(val .≈ value_test) #jl

a_min, a_max = a_grid[[1;end]]
@test all(a_min .< policy .< a_max) #jl

using Plots #md
plot(a_grid, val) #md
#- #md
plot(a_grid, policies_full.a_next) #md
#- #md
plot(a_grid, policies_full.c) #md

#using BenchmarkTools #src

# let value = value, policy=policy, z_mc = z_MC #src
#   lin_ind = LinearIndices(size(value)) #src
#  #src
#   ngp_exo = length(z_mc.state_values) #src
#   n = length(policy) #src
#   len = n * ngp_exo * 2 #src
#  #src
#   I = zeros(Int, len) #src
#   J = zeros(Int, len) #src
#   V = zeros(len) #src
#  #src
#   @btime Aiyagari.controlled_markov_chain!($I, $J, $V, $lin_ind, $z_mc, #src $a_grid, $policy)
# end #src

dist = stationary_distribution(endo, exo, policy)
#926 μs

#writedlm("test/matrices/huggett_dist.txt", dist) #src
dist_test = readdlm(joinpath(PKG_HOME, "test/matrices", "huggett_dist.txt")) #jl
@test maximum(abs, dist .- dist_test) < 1e-14 #jl

# ## Looking at the equilibrium

function excess_demand(r)
  agg_state = HuggettAS(r, a_grid, z_MC)
  @unpack val, policy, policies_full = solve_bellman(endo, exo, agg_state, param, Consumer(), rtol=√eps())

  dist = stationary_distribution(endo, exo, policy)

  sum(dist .* policies_full.a_next)
end

r_vec = 0.07:0.005:0.11 #md

plot(r_vec, excess_demand.(r_vec)) #md

# ## Finite Horizon: Solving the OLG version of the model

function Aiyagari.get_optimum_last(states, agg_state, params, endo, hh::Consumer)
  a_grid = endo.grids.a

  a_next = 0.0
  val    = objective((a_next=a_next,), states, agg_state, x -> 0, params)
  conv   = true

  pol = a_next
  pol_full = (a_next=a_next, c=get_c((a_next=a_next, ), states, agg_state))

  (pol=pol, pol_full=pol_full, val=val, conv=conv)
end

age_grid = 25:65

@unpack val, policy, policies_full = solve_bellman(endo, exo, agg_state, param, Consumer(), age_grid)

using Plots, Colors #md
#md
plt = plot(legend=false) #md
 ran = 1:length(age_grid)- 9 #md
 col = range(colorant"deepskyblue", stop=colorant"navyblue", length=length(ran)) #md
 map(ran) do i #md
   plot!(plt, a_grid, val[:,:,i], color=col[i]) #md
 end #md
 plt #md


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior.Gherkin.Experimental: BadExpectationParseResult

@testset "Combinators          " begin
    @testset "Line" begin
        @testset "Match Foo; Foo; OK" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            p = Line("Foo")
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "Match Foo; Bar; Not OK" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            p = Line("Bar")
            result = p(input)

            # Assert
            @test result isa BadParseResult{String}
            @test result.expected == "Bar"
            @test result.actual == "Foo"
        end

        @testset "Match Foo, then Bar; Foo Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                Bar
            """)

            # Act
            foo = Line("Foo")
            bar = Line("Bar")
            result1 = foo(input)
            result2 = bar(result1.newinput)

            # Assert
            @test result1 isa OKParseResult{String}
            @test result1.value == "Foo"
            @test result2 isa OKParseResult{String}
            @test result2.value == "Bar"
        end

        @testset "Match Foo, then Bar; Foo Baz; Not OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                Baz
            """)

            # Act
            foo = Line("Foo")
            bar = Line("Bar")
            result1 = foo(input)
            result2 = bar(result1.newinput)

            # Assert
            @test result1 isa OKParseResult{String}
            @test result1.value == "Foo"
            @test result2 isa BadParseResult{String}
            @test result2.expected == "Bar"
            @test result2.actual == "Baz"
        end

        @testset "Match Foo; Bar; Not OK, state is unchanged" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            p = Line("Bar")
            result = p(input)

            # Assert
            @test result isa BadParseResult{String}
            @test result.newinput == input
        end

        @testset "Match Foo; No more input; Unexpected EOF" begin
            # Arrange
            input = ParserInput("")

            # Act
            p = Line("Foo")
            result = p(input)

            # Assert
            @test result isa BadUnexpectedEOFParseResult{String}
            @test result.newinput == input
        end
    end

    @testset "Optionally" begin
        @testset "Optionally Foo; Foo; OK" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            parser = Optionally{String}(Line("Foo"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Union{Nothing, String}}
            @test result.value == "Foo"
        end

        @testset "Optionally Foo; Bar; OK with nothing" begin
            # Arrange
            input = ParserInput("Bar")

            # Act
            parser = Optionally{String}(Line("Foo"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Union{Nothing, String}}
            @test result.value === nothing
        end

        @testset "Optionally Bar; Bar; OK" begin
            # Arrange
            input = ParserInput("Bar")

            # Act
            parser = Optionally{String}(Line("Bar"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Union{Nothing, String}}
            @test result.value == "Bar"
        end

        @testset "Optionally Foo then Bar; Foo Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                Bar
            """)

            # Act
            parser1 = Optionally{String}(Line("Foo"))
            result1 = parser1(input)

            parser2 = Line("Bar")
            result2 = parser2(result1.newinput)

            # Assert
            @test result1 isa OKParseResult{Union{Nothing, String}}
            @test result1.value == "Foo"
            @test result2 isa OKParseResult{String}
            @test result2.value == "Bar"
        end

        @testset "Optionally Foo then Bar; Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Bar
            """)

            # Act
            parser1 = Optionally{String}(Line("Foo"))
            result1 = parser1(input)

            parser2 = Line("Bar")
            result2 = parser2(result1.newinput)

            # Assert
            @test result1 isa OKParseResult{Union{Nothing, String}}
            @test result1.value === nothing
            @test result2 isa OKParseResult{String}
            @test result2.value == "Bar"
        end
    end

    @testset "Or" begin
        @testset "Foo or Bar; Foo; OK" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            p = Line("Foo") | Line("Bar")
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "Foo or Bar; Bar; OK" begin
            # Arrange
            input = ParserInput("Bar")

            # Act
            p = Line("Foo") | Line("Bar")
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Bar"
        end

        @testset "Foo or Bar; Baz; Not OK" begin
            # Arrange
            input = ParserInput("Baz")

            # Act
            p = Line("Foo") | Line("Bar")
            result = p(input)

            # Assert
            @test result isa BadParseResult{String}
        end
    end

    @testset "Transformer" begin
        @testset "Transform to Int; 1; OK" begin
            # Arrange
            input = ParserInput("1")

            # Act
            digit = Line("1")
            p = Transformer{String, Int}(digit, x -> parse(Int, x))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Int}
            @test result.value == 1
        end

        @testset "Transform to Int; 2; OK" begin
            # Arrange
            input = ParserInput("2")

            # Act
            digit = Line("1") | Line("2")
            p = Transformer{String, Int}(digit, x -> parse(Int, x))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Int}
            @test result.value == 2
        end
    end

    @testset "Sequence" begin
        @testset "Sequence Line Baz; Baz; OK" begin
            # Arrange
            input = ParserInput("Baz")

            # Act
            p = Sequence{String}(Line("Baz"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Baz"]
        end

        @testset "Sequence Line Baz then Quux; Baz Quux; OK" begin
            # Arrange
            input = ParserInput("""
                Baz
                Quux
            """)

            # Act
            p = Sequence{String}(Line("Baz"), Line("Quux"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Baz", "Quux"]
        end

        @testset "Sequence Line Baz then Quux; Baz Bar; Not OK" begin
            # Arrange
            input = ParserInput("""
                Baz
                Bar
            """)

            # Act
            p = Sequence{String}(Line("Baz"), Line("Quux"))
            result = p(input)

            # Assert
            @test result isa BadParseResult{Vector{String}}
        end

        @testset "Sequence Line Baz then Quux; Foo Quux; Not OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                Quux
            """)

            # Act
            p = Sequence{String}(Line("Baz"), Line("Quux"))
            result = p(input)

            # Assert
            @test result isa BadParseResult{Vector{String}}
        end

        @testset "Sequence Ints 1 then 2; 1 then 2; OK" begin
            # Arrange
            input = ParserInput("""
                1
                2
            """)

            # Act
            digits = Line("1") | Line("2")
            intparser = Transformer{String, Int}(digits, x -> parse(Int, x))
            p = Sequence{Int}(intparser, intparser)
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{Int}}
            @test result.value == [1, 2]
        end
    end

    @testset "Joined" begin
        @testset "Join sequence of Baz or Quux; Baz Quux; OK" begin
            # Arrange
            input = ParserInput("""
                Baz
                Quux
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Joined(Sequence{String}(s, s))
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Baz\nQuux"
        end

        @testset "Join sequence of Baz or Quux; Quux Baz; OK" begin
            # Arrange
            input = ParserInput("""
                Quux
                Baz
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Joined(Sequence{String}(s, s))
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Quux\nBaz"
        end

        @testset "Join sequence of Baz or Quux; Foo Bar; Not OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                Bar
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Joined(Sequence{String}(s, s))
            result = p(input)

            # Assert
            @test result isa BadParseResult{String}
        end
    end

    @testset "Repeating" begin
        @testset "Repeating Baz or Quux; Bar; OK and empty" begin
            # Arrange
            input = ParserInput("""
                Bar
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Repeating{String}(s)
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == []
        end

        @testset "Repeating Baz or Quux; Baz Quux Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Baz
                Quux
                Bar
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Repeating{String}(s)
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Baz", "Quux"]
        end

        @testset "Repeating Baz or Quux followed by Bar; Baz Quux Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Baz
                Quux
                Bar
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Sequence{String}(Joined(Repeating{String}(s)), Line("Bar"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Baz\nQuux", "Bar"]
        end

        @testset "Repeating digits; 3 2 1; OK" begin
            # Arrange
            input = ParserInput("""
                3
                2
                1
                Not a digit
            """)

            # Act
            digits = Line("1") | Line("2") | Line("3")
            intparser = Transformer{String, Int}(digits, x -> parse(Int, x))
            p = Repeating{Int}(intparser)
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{Int}}
            @test result.value == [3, 2, 1]
        end

        @testset "Repeating Baz or Quux, at least 1; Bar; Not OK" begin
            # Arrange
            input = ParserInput("""
                Bar
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Repeating{String}(s, atleast=1)
            result = p(input)

            # Assert
            @test result isa BadParseResult{Vector{String}}
        end

        @testset "Repeating Baz or Quux, at least 1; Baz Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Baz
                Bar
            """)

            # Act
            s = Line("Baz") | Line("Quux")
            p = Repeating{String}(s, atleast=1)
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Baz"]
        end
    end

    @testset "LineIfNot" begin
        @testset "LineIfNot Baz; Baz; Not OK" begin
            # Arrange
            input = ParserInput("""
                Baz
            """)

            # Act
            p = LineIfNot(Line("Baz"))
            result = p(input)

            # Assert
            @test result isa BadParseResult{String}
            @test result.unexpected == "Baz"
        end

        @testset "LineIfNot Baz; Foo Baz; OK" begin
            # Arrange
            input = ParserInput("""
                Foo
            """)

            # Act
            p = LineIfNot(Line("Baz"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "LineIfNot Baz then Baz; Bar Baz; OK" begin
            # Arrange
            input = ParserInput("""
                Bar
                Baz
            """)

            # Act
            p = Sequence{String}(LineIfNot(Line("Baz")), Line("Baz"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Bar", "Baz"]
        end

        @testset "LineIfNot Baz; EOF; Not OK" begin
            # Arrange
            input = ParserInput("")

            # Act
            p = LineIfNot(Line("Baz"))
            result = p(input)

            # Assert
            @test result isa BadUnexpectedEOFParseResult{String}
            @test result.newinput == input
        end
    end

    @testset "StartsWith" begin
        @testset "Foo; Foo; OK" begin
            # Arrange
            input = ParserInput(
                """
                Foo
                """
            )

            # Act
            parser = StartsWith("Foo")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "Foo; Bar; Not OK" begin
            # Arrange
            input = ParserInput(
                """
                Bar
                """
            )

            # Act
            parser = StartsWith("Foo")
            result = parser(input)

            # Assert
            @test result isa BadParseResult{String}
        end

        @testset "Foo; Foo Bar; OK" begin
            # Arrange
            input = ParserInput(
                """
                Foo Bar
                """
            )

            # Act
            parser = StartsWith("Foo")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo Bar"
        end

        @testset "Foo the Quux; Foo Bar, Quux; OK" begin
            # Arrange
            input = ParserInput(
                """
                Foo Bar
                Quux
                """
            )

            # Act
            parser = Sequence{String}(StartsWith("Foo"), Line("Quux"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Foo Bar", "Quux"]
        end

        @testset "Foo; EOF; Not OK" begin
            # Arrange
            input = ParserInput("")

            # Act
            parser = StartsWith("Foo")
            result = parser(input)

            # Assert
            @test result isa BadUnexpectedEOFParseResult{String}
            @test result.newinput == input
        end
    end

    @testset "Whitespace and comments" begin
        @testset "Match Foo; Blank line, then Foo; OK" begin
            # Arrange
            input = ParserInput("""

                Foo
            """)

            # Act
            p = Line("Foo")
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "Match Foo then Bar; Blank line, then Foo, Bar; OK" begin
            # Arrange
            input = ParserInput("""

                Foo
                Bar
            """)

            # Act
            p = Sequence{String}(Line("Foo"), Line("Bar"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value[1] == "Foo"
            @test result.value[2] == "Bar"
        end

        @testset "Match Foo then Bar; Blank line between, then Foo, Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Foo

                Bar
            """)

            # Act
            p = Sequence{String}(Line("Foo"), Line("Bar"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value[1] == "Foo"
            @test result.value[2] == "Bar"
        end

        @testset "Match Foo then Bar; 3 blank line before and between, then Foo, Bar; OK" begin
            # Arrange
            input = ParserInput("""

                Foo

                Bar
            """)

            # Act
            p = Sequence{String}(Line("Foo"), Line("Bar"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value[1] == "Foo"
            @test result.value[2] == "Bar"
        end

        @testset "Match Foo; Comment, then Foo; OK" begin
            # Arrange
            input = ParserInput("""
                # Some comment
                Foo
            """)

            # Act
            p = Line("Foo")
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "Match Foo then Bar; Comment between, then Foo, Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                # Skip this comment
                Bar
            """)

            # Act
            p = Sequence{String}(Line("Foo"), Line("Bar"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value[1] == "Foo"
            @test result.value[2] == "Bar"
        end
    end

    @testset "EOF" begin
        @testset "No non-blank lines left; OK" begin
            # Arrange
            input = ParserInput("")

            # Act
            parser = EOFParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Nothing}
        end

        @testset "Line Foo; Not OK" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            parser = EOFParser()
            result = parser(input)

            # Assert
            @test result isa BadExpectedEOFParseResult{Nothing}
        end

        @testset "Foo, then EOF; OK" begin
            # Arrange
            input = ParserInput("Foo")

            # Act
            parser = Sequence{Union{Nothing, String}}(Line("Foo"), EOFParser())
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{Union{Nothing, String}}}
            @test result.value[1] == "Foo"
            @test result.value[2] === nothing
        end

        @testset "Foo, then Bar; Not OK" begin
            # Arrange
            input = ParserInput("""
                Foo
                Bar
            """)

            # Act
            parser = Sequence{Union{Nothing, String}}(Line("Foo"), EOFParser())
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Vector{Union{Nothing, String}}}
        end

        @testset "Feature, then EOF; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario: Some scenario

                    Scenario: Other scenario
            """)

            # Act
            parser = Sequence{Union{Nothing, Feature}}(FeatureParser(), EOFParser())
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{Union{Nothing, Feature}}}
            @test result.value[1] isa Feature
            @test result.value[2] === nothing
        end

        @testset "Feature, then unallowed new feature; Not OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario: Some scenario

                    Scenario: Other scenario

                    Feature: Not allowed here
            """)

            # Act
            parser = Sequence{Union{Nothing, Feature}}(FeatureParser(), EOFParser())
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Vector{Union{Nothing, Feature}}}
        end
    end
end