text:
function jump_display()
    atol = config.atol

    # config.bound_hint = true

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4
    #       x +2y <= 4
    #       x     >= 0
    #           y >= 0
    #
    # sol: x = 2, y = 0
    # obj_upper = -8
    # obj_lower =  0

    atol = config.atol

    model = BilevelModel()

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    @constraints(Lower(model), begin
        c1, 2x+y <= 4
        c2, x+2y <= 4
        c3, x >= 0
        c4, y >= 0
    end)

    display(x)
    println()
    display(c2)
    println()
    display(model)
    println()
    display(Upper(model))
    println()
    display(Lower(model))
    println()

    xx  = JuMP.variable_by_name(model, "x")
    cc2 = JuMP.constraint_by_name(model, "c2")

    # set_optimizer(model, MOIU.Model{Float64})
    # display(model)
    # println()

end

function invalid_lower_objective()
    model = BilevelModel(CBC_BRIDGED, mode = BilevelJuMP.SOS1Mode())

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    optimize!(model)
    return
end

jump_01vec(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_01(optimizer, true, mode, config)
jump_01(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_01(optimizer, false, mode, config)
function _jump_01(optimizer, vectorized::Bool, mode, config)

    atol = config.atol

    # config.bound_hint = true

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4
    #       x +2y <= 4
    #       x     >= 0
    #           y >= 0
    #
    # sol: x = 2, y = 0
    # obj_upper = -8
    # obj_lower =  0

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -4x -3y)

    @objective(Lower(model), Min, y)

    if vectorized
        @constraints(Lower(model), begin
            c1, 2x+y <= 4
            c2, x+2y <= 4
            c3, x >= 0
            c4, y >= 0
        end)
    else
        @constraint(Lower(model), c1, 2x+y <= 4)
        @constraint(Lower(model), c2, x+2y <= 4)
        @constraint(Lower(model), c3, x >= 0)
        @constraint(Lower(model), c4, y >= 0)
    end

    if config.bound_hint
        for cref in [c1, c2, c3, c4]
            BilevelJuMP.set_dual_upper_bound(cref, 10)
            BilevelJuMP.set_dual_lower_bound(cref, -10)
        end
        BilevelJuMP.set_primal_lower_bound_hint(x, -1)
        BilevelJuMP.set_primal_lower_bound_hint(y, -1)
        BilevelJuMP.set_primal_upper_bound_hint(x, 5)
        BilevelJuMP.set_primal_upper_bound_hint(y, 5)
    end

    # solve twice to check robustness
    optimize!(model)
    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ -8 atol=atol

    @test value(x) ≈  2 atol=atol
    @test value(y) ≈  0 atol=atol

    # @test dual(c1) ≈ [0] atol=atol #NLP fail
    @test dual(c2) ≈ [0] atol=atol
    @test dual(c3) ≈ [0] atol=atol
    # @test dual(c4) ≈ [1] atol=atol #NLP fail

    # solve again to check robustness
    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ -8 atol=atol

    @test value(x) ≈  2 atol=atol
    @test value(y) ≈  0 atol=atol

    # @test dual(c1) ≈ [0] atol=atol #NLP fail
    @test dual(c2) ≈ [0] atol=atol
    @test dual(c3) ≈ [0] atol=atol
    # @test dual(c4) ≈ [1] atol=atol #NLP fail

    # display(x)
    # display(c2)
    # display(model)
    # display(Upper(model))
    # display(Lower(model))

end

function jump_02(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 6)
    @variable(Lower(model), y, start = 2)

    @objective(Upper(model), Min, x + 3y)

    @objective(Lower(model), Min, -y)

    @constraint(Lower(model), c1, x+y <= 8)
    @constraint(Lower(model), c2, x+4y >= 8)
    @constraint(Lower(model), c3, x+2y <= 13)
    @constraint(Lower(model), c4, x >= 1)
    @constraint(Lower(model), c5, x <= 6)
    @constraint(Lower(model), c6, y >= 0)

    if config.bound_hint
        for cref in [c1, c2, c3, c4, c5, c6]
            BilevelJuMP.set_dual_upper_bound(cref, 10)
            BilevelJuMP.set_dual_lower_bound(cref, -10)
        end
        BilevelJuMP.set_primal_upper_bound_hint(y, 10)
        BilevelJuMP.set_primal_lower_bound_hint(y, -1)
        BilevelJuMP.set_primal_lower_bound_hint(x, 0)
        BilevelJuMP.set_primal_upper_bound_hint(x, 7)
    end

    optimize!(model)

    primal_status(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED]

    @test objective_value(model) ≈ 12 atol=atol
    @test BilevelJuMP.lower_objective_value(model) ≈ -2 atol=atol

    @test value(x) ≈ 6 atol=atol
    @test value(y) ≈ 2 atol=atol

    @test dual(c1) ≈ [-1] atol=atol
    @test dual(c2) ≈ [0] atol=atol
    @test dual(c3) ≈ [0] atol=atol
    @test dual(c6) ≈ [0] atol=atol
    @test value(c1) ≈ 8 atol=atol
    @test value(c2) ≈ 14 atol=atol
    @test value(c3) ≈ 10 atol=atol
    @test value(c4) ≈ 6 atol=atol
    @test value(c5) ≈ 6 atol=atol
    @test value(c6) ≈ 2 atol=atol

end

#=
    From:
    Foundations of Bilevel Programming
    by: Stephan Dempe
    in: Nonconvex Optimization and Its Applications
=#

# obs: example 2 is from the book

# Cap 3.2, Pag 25
jump_03(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_03(optimizer, false, mode, config)
jump_03_vec(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_03(optimizer, true, mode, config)
function _jump_03(optimizer, vec::Bool, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x)
    @variable(Upper(model), y)
    if start
        JuMP.set_start_value(x, 3.5*8/15)
        JuMP.set_start_value(y, 8/15)
    end

    @objective(Upper(model), Min, 3x + y)
    @constraint(Upper(model), u1, x <= 5)
    @constraint(Upper(model), u2, y <= 8)
    @constraint(Upper(model), u3, y >= 0)

    @objective(Lower(model), Min, -x)

    if vec
        A = [
             1   1;
            -4  -1;
             2   1;
             2  -7;
        ]
        b = [8, -8, 13, 0]
        @constraint(Lower(model), l,  A*[x, y] .<= b)
        l1 = l[1]
        l3 = l[3]
    else
        @constraint(Lower(model), l1,  x +  y <= 8)
        @constraint(Lower(model), l2, 4x +  y >= 8)
        @constraint(Lower(model), l3, 2x +  y <= 13)
        @constraint(Lower(model), l4, 2x - 7y <= 0)
    end

    if config.bound_hint
        for c in [l1, l2, l3, l4]
            BilevelJuMP.set_dual_upper_bound(c, 15)
            BilevelJuMP.set_dual_lower_bound(c, -15)
        end
        BilevelJuMP.set_primal_lower_bound_hint(x, -10)
        BilevelJuMP.set_primal_upper_bound_hint(x, 6)
        BilevelJuMP.set_primal_lower_bound_hint(y, -1)
        BilevelJuMP.set_primal_upper_bound_hint(y, 9)
    end

    optimize!(model)

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test primal_status(Upper(model)) == MOI.FEASIBLE_POINT
    @test primal_status(Lower(model)) == MOI.FEASIBLE_POINT

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test objective_value(model) ≈ 3* (3.5*8/15) + (8/15) atol=atol
    @test BilevelJuMP.lower_objective_value(model) ≈ -3.5*8/15 atol=atol
    @test objective_value(Lower(model)) ≈ -3.5*8/15 atol=atol

    @test value(x) ≈ 3.5*8/15 atol=atol
    @test value(y) ≈ 8/15 atol=atol
    @test value(u1) ≈ 3.5*8/15 atol=atol
    @test value(l1) ≈ 4.5*8/15 atol=atol

    @test JuMP.dual_status(Lower(model)) == MOI.FEASIBLE_POINT

    @test dual(l1) ≈ [0] atol=atol
    # @test dual(l2) #≈ [0] atol=atol
    @test dual(l3) ≈ [0] atol=atol
    # @show dual(l4) #≈ [0] atol=atol

    if typeof(mode) <: BilevelJuMP.ProductMode
        @test JuMP.dual_status(Upper(model)) == MOI.FEASIBLE_POINT
        @test dual(u1) ≈ 0 atol=atol
    end

end
# change the bound on x to lower level
function jump_04(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x, start = 3.5*8/15)
    @variable(Upper(model), y, start = 8/15)

    @objective(Upper(model), Min, 3x + y)
    @constraint(Upper(model), y <= 8)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, -x)

    @constraint(Lower(model),  x +  y <= 8)
    @constraint(Lower(model), 4x +  y >= 8)
    @constraint(Lower(model), 2x +  y <= 13)
    @constraint(Lower(model), 2x - 7y <= 0)
    @constraint(Lower(model), x <= 5)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 3* (3.5*8/15) + (8/15) atol=1e-4

    @test value(x) ≈ 3.5*8/15 atol=1e-4
    @test value(y) ≈ 8/15 atol=1e-4

end

# Sec 3.3 , pag 30 -> product of x and y in lower level objective

# Sec 3.4.1 , pag 32
function jump_05(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x[i=1:2])
    @variable(Upper(model), y)

    @objective(Upper(model), Min, 2x[2] - y)
    @constraint(Upper(model), y <= 3)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, -x[1])

    @constraint(Lower(model),  100x[1] -x[2] <= 1)
    @constraint(Lower(model),  x[2] <= y)
    @constraint(Lower(model), x[1] >= 0)
    @constraint(Lower(model), x[2] >= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 0 atol=1e-5

    @test value(x[1]) ≈ 1/100 atol=atol
    @test value(x[2]) ≈ 0 atol=atol
    @test value(y) ≈ 0 atol=1e-5

end

# sec 3.5.2.2 pag 44 -> product

# sec 3.7 pag 59
function jump_3SAT(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    n = 7 # maximum literals
    clauses = [[1,2,3],[-1,-4,3],[7,-6,4],[5,6,7]]

    @variable(Lower(model), x[i=1:n])
    @variable(Upper(model), ya[i=1:n])
    @variable(Upper(model), yb[i=1:n])
    @variable(Upper(model), z)
    # @variable(Upper(model), z)

    @objective(Upper(model), Min, sum(x[i] for i in 1:n) - z)
    @constraint(Upper(model), z <= 1)
    @constraint(Upper(model), z >= 0)
    @constraint(Upper(model), [i=1:n], ya[i] >= 0)
    @constraint(Upper(model), [i=1:n], ya[i] <= 1)
    @constraint(Upper(model), [i=1:n], yb[i] >= 0)
    @constraint(Upper(model), [i=1:n], yb[i] <= 1)
    @constraint(Upper(model), [i=1:n], ya[i] + yb[i] == 1)
    for c in clauses
        @constraint(Upper(model),
            sum(i > 0 ? ya[i] : yb[-i] for i in c) >= z)
    end

    @objective(Lower(model), Min, -sum(x[i] for i in 1:n))

    @constraint(Lower(model), [i=1:n], x[i] >= 0)
    @constraint(Lower(model), [i=1:n], x[i] <= ya[i])
    @constraint(Lower(model), [i=1:n], x[i] <= yb[i])

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ -1 atol=atol

    # 3SAT is yese IFF obj = -1

end

# sec 5.1 pag 121 -> product

# sec 5.1 pag 127
function jump_quad_01_a(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x, start = 0)
    @variable(Upper(model), y, start = 0)

    @objective(Upper(model), Min, x^2 + y)
    @constraint(Upper(model), u1, -x -y <= 0)

    @objective(Lower(model), Min, x)
    @constraint(Lower(model), l1, x >= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 0 atol=atol
    @test value(x) ≈ 0 atol=atol
    @test value(y) ≈ 0 atol=atol

end
function jump_quad_01_b(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x, start = 0.5)
    @variable(Upper(model), y, start = -0.5)

    @objective(Upper(model), Min, x^2 + y)

    @objective(Lower(model), Min, x)
    @constraint(Lower(model), l1, -x -y <= 0)
    @constraint(Lower(model), l2, x >= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 0.5^2 - 0.5 atol=atol
    @test value(x) ≈ 0.5 atol=atol
    @test value(y) ≈ -0.5 atol=atol

end
function jump_quad_01_c(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x)
    @variable(Upper(model), y)

    @objective(Upper(model), Min, x^2 + y)
    @constraint(Upper(model), -x -y <= 0)
    @constraint(Upper(model), x >= 0)

    @objective(Lower(model), Min, x)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # lower level migh not be always feasible
    # KKT and global methods should work

end

# sec 5.1 pag 126 -> product

# sec 5.2 pag 130 -> quadracti second level! (good test for new KKT)
# sec 5.5 pag 149 -> quadracti
# sec 5.7.1 pag 159, 161, 169, 171 -> quadracti

# sec 6.1 pag 196 -> quadratic objective in second level
# sec 7.2 pag 226,227, 233 -> quad

# sec 8.1 pag 255 example from [279]
# only the second level is described
function jump_int_01(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x)
    @variable(Upper(model), y)

    @objective(Upper(model), Min, 0)
    # @constraint(Upper(model), -x -y <= 0)
    # @constraint(Upper(model), x >= 0)

    # include integrality

    @objective(Lower(model), Min, x)

    @constraint(Lower(model), x + y <= 2)
    @constraint(Lower(model), x - y <= 2)
    @constraint(Lower(model), -4x + 5y <= 10)
    @constraint(Lower(model), -4x - 5y <= 10)

    optimize!(model)

    primal_status(model)

    termination_status(model)

end

# sec 8.1 pag 257
function jump_int_02(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Lower(model), x)
    @variable(Upper(model), y)

    @objective(Upper(model), Min, -10x - y)

    # include integrality on Y

    @objective(Lower(model), Min, x)

    @constraint(Lower(model), 20x - 25y <= 30)
    @constraint(Lower(model), 2x + y <= 10)
    @constraint(Lower(model), -x + 2y <= 15)
    @constraint(Lower(model), 10x + 2y >= 15)

    # DONT include integrality on X (KKT wont work)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # continuous solution x = 1, y = 8
    # same for integer on upper
    # integer on both levels is 2,2

end

# sec 8.3 pag 265 -> products + integrality on lower

#=
    From:
    Practical Bilevel Optimization
    by: J. F. Bard
    in:
=#

# pag 197 ex 5.1.1
jump_06(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_06(optimizer, false, mode, config)
jump_06_sv(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_06(optimizer, true, mode, config)
function _jump_06(optimizer, sv::Bool, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    if sv
        @variable(Upper(model), x >= 0, start = 4)
        @variable(Lower(model), y >= 0, start = 4)
    else
        @variable(Upper(model), x, start = 4)
        @variable(Lower(model), y, start = 4)
    end

    @objective(Upper(model), Min, x - 4y)
    if !sv
        @constraint(Upper(model), x >= 0)
    end

    @objective(Lower(model), Min, y)

    @constraint(Lower(model), -x  - y <= -3)
    @constraint(Lower(model), -2x + y <= 0)
    @constraint(Lower(model), 2x  + y <= 12)
    @constraint(Lower(model), 3x + -2y <= 4) # signs are wrong in some versions of the book
    if !sv
        @constraint(Lower(model), y >= 0)
    end

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ -12 atol=atol

    @test value(x) ≈ 4 atol=atol
    @test value(y) ≈ 4 atol=atol

end

# pag 208 ex 5.3.1
function jump_07(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2], start = [0, 0.9][i])
    @variable(Lower(model), y[i=1:3], start = [0, 0.6, 0.4][i])

    @objective(Upper(model), Min,
        -8x[1] -4x[2] + 4y[1] - 40y[2] - 4y[3])
    @constraint(Upper(model), [i=1:2], x[i] >= 0)

    @objective(Lower(model), Min,
        x[1]+ 2x[2] + y[1] + y[2] + 2y[3])
    @constraint(Lower(model), [i=1:3], y[i] >= 0)

    @constraint(Lower(model), -y[1] + y[2] + y[3] <= 1)
    @constraint(Lower(model), 2x[1] - y[1] + 2y[2] - 0.5y[3] <= 1)
    @constraint(Lower(model), 2x[2] +2y[1] - y[2] - 0.5y[3] <= 1)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value.(x) ≈ [0, 0.9] atol=atol
    @test value.(y) ≈ [0, 0.6, 0.4] atol=atol

end

# pag 208 ex 5.3.1
function jump_08(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 8/9)
    @variable(Lower(model), y, start = 20/9)

    @objective(Upper(model), Min, x + y)
    @constraint(Upper(model), x >= 0)

    @objective(Lower(model), Min, -5x - y)
    @constraint(Lower(model), y >= 0)

    @constraint(Lower(model), -x - y/2 <= -2)
    @constraint(Lower(model), -x/4 + y <= 2)
    @constraint(Lower(model), x + y/2 <= 8)
    @constraint(Lower(model), x - 2y <= 4)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 8/9 atol=atol
    @test value(y) ≈ 20/9 atol=atol

end

# pag 271 ex 7.1.1 -> quadratic terms
# pag 281 ex 7.2.2 -> quadratic terms
# pag 302 ex 8.1.1 -> quadratic terms

# pag 304 ex 8.1.2
function jump_09a(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    # degenerate second level
    # its case that show that the KKT approach is OPTIMISTIC

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 0.0)
    @variable(Lower(model), y[i=1:2], start = [0, 1][i])

    @objective(Upper(model), Min, -10.1x + 10y[1] - y[2])
    @constraint(Upper(model), x >= 0)

    @constraint(Upper(model), y[1] >= 0)
    @constraint(Upper(model), y[2] >= 0)

    @objective(Lower(model), Min, -y[1] - y[2])
    @constraint(Lower(model), y[1] >= 0)
    @constraint(Lower(model), y[2] >= 0)

    @constraint(Lower(model), x - y[1] <= 1)
    @constraint(Lower(model), x + y[2] <= 1)
    @constraint(Lower(model), y[1] + y[2] <= 1)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 0 atol=atol
    @test value.(y) ≈ [0, 1] atol=atol

end
function jump_09b(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y[i=1:2])

    @objective(Upper(model), Min, -x + 10y[1] - 2y[2])
    @constraint(Upper(model), x >= 0)

    @objective(Lower(model), Min, -y[1] - y[2])
    @constraint(Lower(model), y[1] >= 0)
    @constraint(Lower(model), y[2] >= 0)

    @constraint(Lower(model), x - y[1] <= 1)
    @constraint(Lower(model), x + y[2] <= 1)
    @constraint(Lower(model), y[1] + y[2] <= 1)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 0 atol=atol
    @test value.(y) ≈ [0, 1] atol=atol

end

# pag 306 ex 8.1.3 -> quardratics
# pag 336 ex 8.5.1 -> quardratics
# pag 358 ex 8.7.1 -> quardratics

#=
    From:
    Bilevel programming problems
    by: S. Dempe et al.
    in: Energy systems
=#

# pag 4 ex 1.2 -> product

# pag 9 ex 1.4 TODO obtain solution
function jump_10(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    n = 10
    f = 100
    b_l = 0
    b_u = 50

    c_lw =
    [10 30 8 60 4 16 32 30 120 6]
    a_lw =
    [5 3 2 5 1 8 4 3 6 3]
    c_up =
    [10 15 -24 20 -40 80 -32 -60 -12 -60]

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), b)
    @variable(Lower(model), x[i=1:n])

    @objective(Upper(model), Min, sum(c_up[i]*x[i] for i in 1:n) + f*b)
    @constraint(Upper(model), b >= b_l)
    @constraint(Upper(model), b <= b_u)

    @objective(Lower(model), Min, sum(c_lw[i]*x[i] for i in 1:n) )
    @constraint(Lower(model), [i=1:n], x[i] >= 0)
    @constraint(Lower(model), [i=1:n], x[i] <= 1)
    @constraint(Lower(model), sum(a_lw[i]*x[i] for i in 1:n) >= b)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # @test value(x) ≈ 0
    # @test value(y) ≈ [1, 0]

end

# pag 21 ex 2.1
function jump_11a(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 8)
    @variable(Lower(model), y, start = 6)

    @objective(Upper(model), Min, -x - 2y)
    @constraint(Upper(model), 2x - 3y >= -12)
    @constraint(Upper(model), x + y <= 14)

    @objective(Lower(model), Min, -y)
    @constraint(Lower(model), -3x + y <= -3)
    @constraint(Lower(model), 3x + y <= 30)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 8 atol=atol
    @test value(y) ≈ 6 atol=atol

end
function jump_11b(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 6)
    @variable(Lower(model), y, start = 8)

    @objective(Upper(model), Min, -x - 2y)

    @objective(Lower(model), Min, -y)
    @constraint(Lower(model), 2x - 3y >= -12)
    @constraint(Lower(model), x + y <= 14)
    @constraint(Lower(model), -3x + y <= -3)
    @constraint(Lower(model), 3x + y <= 30)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 6 atol=atol
    @test value(y) ≈ 8 atol=atol

end

# pag 45 ex 3.3
function jump_12(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    for a in [0 0.1 0.2]
        atol = config.atol

        MOI.empty!(optimizer)
        model = BilevelModel(()->optimizer, mode = mode)

        @variable(Upper(model), x)
        @variable(Lower(model), y)
        @variable(Lower(model), z)

        @objective(Upper(model), Min, 0.5x - y + 3z)

        @objective(Lower(model), Min, -y - z)
        @constraint(Lower(model), x + y <= 1 + a)
        @constraint(Lower(model), -x + y <= 1)
        @constraint(Lower(model), z <= 1)
        @constraint(Lower(model), z >= 0)

        MOI.empty!(optimizer)
        @test MOI.is_empty(optimizer)

        optimize!(model)

        primal_status(model)

        termination_status(model)

        @test value(x) ≈ a/2 atol=atol
        @test value(y) ≈ 1 + a/2 atol=atol
        @test value(z) ≈ 1 atol=atol
    end
end

# pag 45 ex 3.3
function jump_13_quad(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 3/2)
    @variable(Lower(model), y, start = 1/2)

    @objective(Upper(model), Min, (x -1 )^2 + y^2)
    @constraint(Upper(model), x >= -2)
    @constraint(Upper(model), x <= +2)

    @objective(Lower(model), Min, -y)
    @constraint(Lower(model), x + y <= 2)
    @constraint(Lower(model), -x + y <= 2)
    @constraint(Lower(model), y >= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 3/2 atol=atol
    @test value(y) ≈ 1/2 atol=atol
end

# pag 290 ex 8.2
function jump_14(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Upper(model), y)
    @variable(Lower(model), z)

    @objective(Upper(model), Min, x - 2y + z)
    @constraint(Upper(model), x + y +z >= 15)
    @constraint(Upper(model), x <= 10)
    @constraint(Upper(model), y <= 10)
    @constraint(Upper(model), z <= 10)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), y >= 0)
    @constraint(Upper(model), z >= 0)

    @objective(Lower(model), Min, 2x-y+z)
    @constraint(Lower(model), x + y -z <= 5)
    @constraint(Lower(model), z >= 0)
    @constraint(Lower(model), z <= 10)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 0 atol=atol
    @test value(y) ≈ 10 atol=atol
    @test value(z) ≈ 5 atol=atol
end

# pag 300 ex 8.5.2
function jump_15a_INT(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)
    @variable(Lower(model), z)

    @objective(Upper(model), Min, x - 2y + z)
    @constraint(Upper(model), x <= 1)
    @constraint(Upper(model), x >= 0)
    # ADD INTEGRALITY (binary) ON X
    @constraint(Upper(model), y <= 100)
    @constraint(Upper(model), z >= 0)
    @constraint(Upper(model), z <= 100)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, - 60y - 8z)
    @constraint(Lower(model), 10x + 2y + 3z <= 225)
    @constraint(Lower(model),  5x + 3y + 0z <= 230)
    @constraint(Lower(model),  5x + 0y + 1z <= 85)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 75 atol=atol
    @test value(Z) ≈ 21+2/3 atol=atol
end
function jump_15b_INT(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Upper(model), y)
    @variable(Lower(model), z)

    @objective(Upper(model), Min, x - 2y + z)
    @constraint(Upper(model), x <= 1)
    @constraint(Upper(model), x >= 0)
    # ADD INTEGRALITY (binary) ON X
    @constraint(Upper(model), y <= 100)
    @constraint(Upper(model), z >= 0)
    @constraint(Upper(model), z <= 100)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, - 60y - 8z)
    @constraint(Lower(model), 10x + 2y + 3z <= 225)
    @constraint(Lower(model),  5x + 3y + 0z <= 230)
    @constraint(Lower(model),  5x + 0y + 1z <= 85)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 75 atol=atol
    @test value(Z) ≈ 21+2/3 atol=atol
end

#=
    Princeton Handbook of Test Problems in Local and Global Optimization
    Floudas c., Pardalos P., et al.
    HTP
=#

# 9.2.2
function jump_HTP_lin01(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y[i=1:2])
    if start
        JuMP.set_start_value(x, 5.0)
        JuMP.set_start_value(y[1], 4.0)
        JuMP.set_start_value(y[2], 2.0)
    end

    @objective(Upper(model), Min, -x - 3y[1] + 2y[2])
    @constraint(Upper(model), x <= 8)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), [i=1:2], y[i] <= 4)
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min, - y[1])
    @constraint(Lower(model), -y[1] <= 0)
    @constraint(Lower(model),  y[1] <= 4)
    @constraint(Lower(model), -2x + 2y[1] + 4y[2] <= 16)
    @constraint(Lower(model),  8x + 3y[1] - 2y[2] <= 48)
    @constraint(Lower(model), -2x +  y[1] - 3y[2] <= -12)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 5 atol=atol
    @test value.(y) ≈ [4, 2] atol=atol
end

# 9.2.3
function jump_HTP_lin02(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 4)
    @variable(Lower(model), y, start = 4)

    @objective(Upper(model), Min, -x - 3y)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, y)
    @constraint(Lower(model), c1, -y <= 0)
    @constraint(Lower(model), c2, -x + y <= 3)
    @constraint(Lower(model),   x + 2y <= 12)
    @constraint(Lower(model),  4x -  y <= 12)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 4 atol=atol
    @test value(y) ≈ 4 atol=atol
    @test dual(c1) ≈ [0] atol=atol
    @test dual(c2) ≈ [0] atol=atol
end

# 9.2.4 - parg 211
jump_HTP_lin03(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_HTP_lin03(optimizer, false, mode, config)
jump_HTP_lin03_vec(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_HTP_lin03(optimizer, true, mode, config)
function _jump_HTP_lin03(optimizer, vec::Bool, mode, config)
    # send y to upper level

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2], start = [0, 0.9][i])
    @variable(Lower(model), y[i=1:6], start = [0, 0.6, 0.4, 0, 0, 0][i])

    @objective(Upper(model), Min,
        4y[1] - 40y[2] -4y[3] -8x[1] -4x[2])
    @constraint(Upper(model), [i=1:2], x[i] >= 0)
    @constraint(Upper(model), [i=1:6], y[i] >= 0)

    @objective(Lower(model), Min,
        y[1] + y[2] + 2y[3] +x[1] +2x[2])

    @constraint(Lower(model), [i=1:6], y[i] >= 0)
    if vec
        H1 = [ -1  1  1   1  0  0;
               -1  2 -1/2 0  1  0;
                2 -1 -1/2 0  0  1
            ]
        H2 = [0 0;
              2 0;
              0 2
            ]
        b = [1, 1, 1]
        @constraint(Lower(model), H1*y + H2*x .== b)
    else
        @constraint(Lower(model), c1, -y[1] +  y[2] +       y[3] + y[4]         == 1)
        @constraint(Lower(model), c2, -y[1] + 2y[2] - (1/2)*y[3] + y[5] + 2x[1] == 1)
        @constraint(Lower(model), c3, 2y[1] -  y[2] - (1/2)*y[3] + y[6] + 2x[2] == 1)
    end

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value.(x) ≈ [0, 0.9] atol=atol
    @test value.(y) ≈ [0, 0.6, 0.4, 0, 0, 0] atol=atol
end

# 9.2.5
function jump_HTP_lin04(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, x - 4y)
    # @constraint(Upper(model), x >= 0)
    # @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, y)
    @constraint(Lower(model), -2x +  y <= 0)
    @constraint(Lower(model),  2x + 5y <= 108)
    @constraint(Lower(model),  2x - 3y <= -4)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 19 atol=atol
    @test value(y) ≈ 14 atol=atol
end

# 9.2.6
function jump_HTP_lin05(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 0)
    @variable(Lower(model), y[i=1:2], start = [0, 1][i])

    @objective(Upper(model), Min, -x + 10y[1] - y[2])
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min, -y[1] - y[2])
    @constraint(Lower(model),    x -  y[1] <= 1)
    @constraint(Lower(model),    x +  y[2] <= 1)
    @constraint(Lower(model), y[1] +  y[2] <= 1)
    @constraint(Lower(model),  -y[1] <= 0)
    @constraint(Lower(model),  -y[2] <= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # book solution looks incorrect
    # the actual solutions seems to be
    @test value(x) ≈ 0 atol=atol
    @test value.(y) ≈ [0, 1] atol=atol
end

# 9.2.7
function jump_HTP_lin06(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 16)
    @variable(Lower(model), y, start = 11)

    @objective(Upper(model), Min, -x - 3y)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, -x + 3y)
    @constraint(Lower(model),   - x - 2y <= -10)
    @constraint(Lower(model),     x - 2y <= 6) # sign of 2y was wrong on book
    @constraint(Lower(model),    2x -  y <= 21)
    @constraint(Lower(model),     x + 2y <= 38)
    @constraint(Lower(model),    -x + 2y <= 18)
    @constraint(Lower(model),         -y <= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # book solution looks incorrect
    # the actual solutions seems to be
    @test value(x) ≈ 16 atol=atol
    @test value(y) ≈ 11 atol=atol
end

# 9.2.8 - parg 216
function jump_HTP_lin07(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol
    rtol = config.rtol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2], start = [0, 0.9][i])
    @variable(Lower(model), y[i=1:3], start = [0, 0.6, 0.4][i])

    @objective(Upper(model), Min,
        -8x[1] -4x[2] + 4y[1] - 40y[2] + 4y[3])
        # the sign of 4y[3] seems off, but solution is the same
        # model from the book left as is
    @constraint(Upper(model), [i=1:2], x[i] >= 0)
    @constraint(Upper(model), [i=1:3], y[i] >= 0)

    @objective(Lower(model), Min,
        x[1] + 2x[2] + y[1] + y[2] + 2y[3])

    @constraint(Lower(model), [i=1:3], y[i] >= 0)

    @constraint(Lower(model),         -  y[1] +  y[2] +       y[3] <= 1)
    @constraint(Lower(model), + 2x[1] -  y[1] + 2y[2] - (1/2)*y[3] <= 1)
    @constraint(Lower(model), + 2x[2] + 2y[1] -  y[2] - (1/2)*y[3] <= 1)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value.(x) ≈ [0, 0.9] atol=atol rtol=rtol
    @test value.(y) ≈ [0, 0.6, 0.4] atol=atol rtol=rtol
end

# 9.2.9 - parg 217
function jump_HTP_lin08(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2])
    @variable(Lower(model), y[i=1:2])

    @objective(Upper(model), Min,
        -2x[1] + x[2] + 0.5*y[1])
    @constraint(Upper(model), [i=1:2], x[i] >= 0)
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min,
        x[1] + x[2] - 4y[1] + y[2])

    @constraint(Lower(model), [i=1:2], y[i] >= 0)

    @constraint(Lower(model), - 2x[1]         +  y[1] -  y[2] <= -2.5)
    @constraint(Lower(model), +  x[1] - 3x[2]         +  y[2] <= 2)
    @constraint(Lower(model), +  x[1] +  x[2]                 <= 2)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # solution frm the book is incorrect
    # using solution from: https://www.gams.com/latest/emplib_ml/libhtml/emplib_flds918.html
    @test value.(x) ≈ [2, 0] atol=atol
    @test value.(y) ≈ [1.5, 0] atol=atol
end

# 9.2.10
function jump_HTP_lin09(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, x + y)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), y >= 0)

    @objective(Lower(model), Min, -5x - y)
    @constraint(Lower(model),   - x - 0.5*y <= -2)
    @constraint(Lower(model), -0.25*x + y <= 2) # sign of 2y was wrong on book
    @constraint(Lower(model),     x + 0.5*y <= 8)
    @constraint(Lower(model),     x - 2y <= 2)
    @constraint(Lower(model),         -y <= 0)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 8/9 atol=atol#0.888_888_888_888
    @test value(y) ≈ 20/9 atol=atol# 2.222_222_222_222
end

# 9.2.11 - parg 219
function jump_HTP_lin10(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2])
    @variable(Lower(model), y[i=1:2])

    @objective(Upper(model), Min,
        -2x[1] + x[2] + 0.5*y[1])
    @constraint(Upper(model), [i=1:2], x[i] >= 0)
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    # this upper constraint is missing in the book but is correct in the GAMS file
    @constraint(Upper(model), +  x[1] +  x[2]  <= 2)

    @objective(Lower(model), Min,
        x[1] + x[2] - 4y[1] + y[2])

    @constraint(Lower(model), [i=1:2], y[i] >= 0)

    @constraint(Lower(model), - 2x[1]         +  y[1] -  y[2] <= -2.5)
    @constraint(Lower(model), +  x[1] - 3x[2]         +  y[2] <= 2)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    # solution frm the book is incorrect
    # using solution from: https://www.gams.com/latest/emplib_ml/libhtml/emplib_flds918.html
    @test value.(x) ≈ [2, 0] atol=atol
    @test value.(y) ≈ [1.5, 0] atol=atol
end

# TODO - add quadratic problems

# 9.3.2 - parg 221
function jump_HTP_quad01(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min,
        (x-5)^2 + (2y+1)^2)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), y >= 0) # only in lowrrin GAMS

    @objective(Lower(model), Min,
        (y-1)^2 -1.5*x*y)

    @constraint(Lower(model), -3x +    y <= -3)
    @constraint(Lower(model),   x - 0.5y <= 4)
    @constraint(Lower(model),   x +    y <= 7)
    # @constraint(Lower(model), y >= 0) # GAMS file has this constraint

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 0 atol=atol
end

# 9.3.3- parg 222
function jump_HTP_quad02(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min,
        x^2 + (y-10)^2)
    @constraint(Upper(model), - x + y <= 0)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), x <= 15)
    @constraint(Upper(model), y >= 0)
    @constraint(Upper(model), y <= 20)

    @objective(Lower(model), Min,
        (x + 2y - 30)^2)

    @constraint(Lower(model),   x +    y <= 20)
    @constraint(Lower(model),     -    y <= 0)
    @constraint(Lower(model),          y <= 20)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 10 atol=1e-4
    @test value(y) ≈ 10 atol=1e-4
end

# 9.3.4- parg 223
function jump_HTP_quad03(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2], start = 0)
    @variable(Lower(model), y[i=1:2], start = -10)

    @objective(Upper(model), Min,
        2x[1] + 2x[2] -3y[1] - 3y[2] -60)
    @constraint(Upper(model), x[1] + x[2] + y[1] -2y[2] -40 <= 0)
    @constraint(Upper(model), [i=1:2], x[i] >= 0)
    @constraint(Upper(model), [i=1:2], x[i] <= 50)
    @constraint(Upper(model), [i=1:2], y[i] >= -10)
    @constraint(Upper(model), [i=1:2], y[i] <= 20)

    @objective(Lower(model), Min,
        (-x[1] + y[1] + 40)^2 + (-x[2] + y[2] + 20)^2)
        # the boo does not contain the 40, it is a 20 there
        # however the solution does not match
        # the file https://www.gams.com/latest/emplib_ml/libhtml/emplib_flds923.html
        # has a 40

    @constraint(Lower(model),  [i=1:2],- x[i] + 2y[i] <= -10)
    @constraint(Lower(model), [i=1:2], y[i] >= -10)
    @constraint(Lower(model), [i=1:2], y[i] <= 20)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    @test objective_value(model) ≈ 0 atol=atol

    sol = vcat(value.(x), value.(y))
    @test sol ≈ [0 ; 0 ; -10; -10] || sol ≈ [0 ; 30; -10; 10] #atol=atol
    # gurobi found the second solution  which is actually feasible and
    # has the same objective value
end

# 9.3.5- parg 225
function jump_HTP_quad04(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 3)
    @variable(Lower(model), y[i=1:2], start = [1, 2][i])

    @objective(Upper(model), Min,
        0.5*((y[1]-2)^2+(y[2]-2)^2) )
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min,
        0.5*(y[1]^2) + y[2])

    @constraint(Lower(model), y[1] + y[2] == x)
    @constraint(Lower(model), [i=1:2], y[i] >= 0)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    @test objective_value(model) ≈ 0.5 atol=atol

    @test value(x) ≈ 3 atol=atol
    @test value.(y) ≈ [1 ; 2] atol=atol
end

# 9.3.6- parg 226
function jump_HTP_quad05(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 1)
    @variable(Lower(model), y, start = 3)

    @objective(Upper(model), Min,
        (x-3)^2+(y-2)^2 )
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), x <= 8)

    @objective(Lower(model), Min,
        (y-5)^2 )

    @constraint(Lower(model), -2x+y <= 1)
    @constraint(Lower(model),  x-2y <= 2)
    @constraint(Lower(model),  x+2y <= 14)
    @constraint(Lower(model), y >= 0)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    @test objective_value(model) ≈ 5 atol=1e-5

    @test value(x) ≈ 1 atol=1e-5
    @test value(y) ≈ 3 atol=1e-5
end

# 9.3.7- parg 227
function jump_HTP_quad06(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2], start = 0.5)
    @variable(Lower(model), y[i=1:2], start = 0.5)

    @objective(Upper(model), Min,
        x[1]^2 -2x[1] +x[2]^2 - 2x[2] +y[1]^2 +y[2]^2)
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min,
        (-x[1] + y[1])^2 + (-x[2] + y[2])^2)

    @constraint(Lower(model), [i=1:2], y[i] >= 0.5)
    @constraint(Lower(model), [i=1:2], y[i] <= 1.5)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    @test objective_value(model) ≈ -1 atol=atol

    sol = vcat(value.(x), value.(y))
    @test sol ≈ [0.5 ; 0.5; 0.5; 0.5] atol=atol

end
function jump_HTP_quad06b(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # TODO reviews the behaviour comment
    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:2], start = 0.5*sqrt(3))
    @variable(Lower(model), y[i=1:2], start = 0.5*sqrt(3))

    @objective(Upper(model), Min,
        x[1]^2 +x[2]^2 +y[1]^2 - 3y[1] +y[2]^2 - 3y[2])
    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min,
        (-x[1] + y[1])^2 + (-x[2] + y[2])^2)

    @constraint(Lower(model), [i=1:2], y[i] >= 0.5)
    @constraint(Lower(model), [i=1:2], y[i] <= 1.5)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    # @test objective_value(model) ≈ -1 atol=atol

    sol = vcat(value.(x), value.(y))
    @test sol ≈ [0.5 ; 0.5; 0.5; 0.5]*sqrt(3) atol=atol

end

# 9.3.8- parg 228
function jump_HTP_quad07(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min,
        (x-5)^2 + (2y+1)^2)
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), y >= 0) # only in lowrrin GAMS

    @objective(Lower(model), Min,
        (y-1)^2 -1.5*x*y)

    @constraint(Lower(model), -3x +    y <= -3)
    @constraint(Lower(model),   x - 0.5y <= 4)
    @constraint(Lower(model),   x +    y <= 7)
    @constraint(Lower(model), y >= 0) # only difference from problem 1

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 0 atol=atol
end

# 9.3.9 - parg 229
function jump_HTP_quad08(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # Q objective is not PSD

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 0.25)
    @variable(Lower(model), y, start = 0)

    @objective(Upper(model), Min,
        -(4x-3)*y+(2x+1) )
    @constraint(Upper(model), x >= 0)
    @constraint(Upper(model), x <= 1)
    @constraint(Upper(model), y >= 0)
    @constraint(Upper(model), y <= 1)

    @objective(Lower(model), Min,
        -(1-4x)*y -(2x+2) )

    @constraint(Lower(model), y >= 0)
    @constraint(Lower(model), y <= 1)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    @test objective_value(model) ≈ 1.5 atol=atol

    @test value(x) ≈ 0.25 atol=atol
    @test value(y) ≈ 0 atol=atol
end

# 9.3.10- parg 230
function jump_HTP_quad09(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 2)
    @variable(Lower(model), y[i=1:2], start = [6, 0][i])

    @objective(Upper(model), Min,
        x+y[2] )
    @constraint(Upper(model), [i=1:2], y[i] >= 0)
    @constraint(Upper(model), x >= 0)

    @objective(Lower(model), Min,
        2y[1]+x*y[2])

    @constraint(Lower(model), x + 4 <= y[1] + y[2])
    # this 4 is missing from the book
    # see https://www.gams.com/latest/emplib_ml/libhtml/emplib_flds929.html
    @constraint(Lower(model), [i=1:2], y[i] >= 0)

    optimize!(model)

    primal_status(model)
    termination_status(model)
    @test objective_value(model) ≈ 2 atol=atol

    @test value(x) ≈ 2 atol=atol
    @test value.(y) ≈ [6 ; 0] atol=atol
end

#=
    GAMS educational examples by Michael Ferris
=#

# from https://www.gams.com/latest/emplib_ml/libhtml/emplib_jointc1.html
function jump_jointc1(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    ###
    ### PART 1
    ###

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, x)
    @constraint(Upper(model), x >= 1)
    @constraint(Upper(model), y >= +x)
    @constraint(Upper(model), y >= -x)

    # this lower problem leads to x == -x (infeasible because x >= 1)
    @objective(Lower(model), Min, y)
    @constraint(Lower(model), y >= -x)

    optimize!(model)

    @test termination_status(model) in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED, MOI.LOCALLY_INFEASIBLE]

    ###
    ### PART 2
    ###

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 1)
    @variable(Lower(model), y, start = 1)

    @objective(Upper(model), Min, x)
    @constraint(Upper(model), x >= 1)
    @constraint(Upper(model), y >= +x)
    @constraint(Upper(model), y >= -x)

    # this lower problem leads to x == -x (infeasible because x >= 1)
    @objective(Lower(model), Min, y)
    @constraint(Lower(model), y >= -x)
    @constraint(Lower(model), y >= +x)

    optimize!(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED]

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 1 atol=atol
end

# from https://www.gams.com/latest/emplib_ml/libhtml/emplib_jointc2.html
function jump_jointc2(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    ###
    ### PART 1
    ###

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 1)
    @variable(Lower(model), y, start = -1)

    @objective(Upper(model), Max, x)

    @constraint(Upper(model), x >= 1)
    @constraint(Upper(model), x <= 2)
    @constraint(Upper(model), y <= 100)

    @constraint(Upper(model), y >= x - 2)
    @constraint(Upper(model), y >= -x)

    # this lower problem leads to x == -x (infeasible because x >= 1)
    @objective(Lower(model), Min, y)
    @constraint(Lower(model), y >= -x)

    optimize!(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test value(x) ≈ +1 atol=atol
    @test value(y) ≈ -1 atol=atol

    ###
    ### PART 2
    ###

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 2)
    @variable(Lower(model), y, start = 0)

    @objective(Upper(model), Max, x)

    @constraint(Upper(model), x >= 1)
    @constraint(Upper(model), x <= 2)
    @constraint(Upper(model), y <= 100)

    @constraint(Upper(model), y >= x - 2)
    @constraint(Upper(model), y >= -x)

    # this lower problem leads to x == -x (infeasible because x >= 1)
    @objective(Lower(model), Min, y)
    @constraint(Lower(model), y >= -x)
    @constraint(Lower(model), y >= x - 2)

    optimize!(model)

    @test termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

    @test value(x) ≈ 2 atol=atol
    @test value(y) ≈ 0 atol=atol
end

#=
    An Efficient Point Algorithm for Two-Stage Optimization Problem
    J.F. Bard
    Operations Research, 1983
=#
function jump_EffPointAlgo(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())
    # send y to upper level

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:6])
    @variable(Lower(model), y[i=1:3])

    @objective(Upper(model), Max,
        2x[1] - x[2] - x[3] + 2x[4] + x[5] -3.5x[6]
        - y[1] -1.5y[2] + 3y[3])
    @constraint(Upper(model), [i=1:6], x[i] >= 0)

    @objective(Lower(model), Max,
        2x[2] - x[5] + 3y[1] -y[2] -4y[3])
    @constraint(Lower(model), [i=1:3], y[i] >= 0)
    @constraint(Lower(model), -  x[1] +0.2x[2]                + x[5] +  2x[6] - 4y[1] + 2y[2] + y[3] <= 12)
    @constraint(Lower(model),    x[1]          + x[3] - 2x[4]                         - 4y[2] + y[3] <= 10)
    @constraint(Lower(model),   5x[1]                 +  x[4]        +3.2x[6] + 2y[1] + 2y[2]        <= 15)
    @constraint(Lower(model),         -  3x[2]        -  x[4] + x[5]          - 2y[1]                <= 12)
    @constraint(Lower(model), - 2x[1] -   x[2]                                        -  y[2] + y[3] <= -2)
    @constraint(Lower(model),                                                 -  y[1] - 2y[2] - y[3] <= -2)
    @constraint(Lower(model),         -  2x[2] - 3x[3]        - x[5]                                 <= -3)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value.(x) ≈ [0, 4, 0, 15, 9.2, 0] atol=atol
    @test value.(y) ≈ [0, 0, 2] atol=atol
end

#=
    Test Problem Construction for Linear Bilevel Programming Problems
    K. Moshirvaziri et al.
    Journal of Global Optimization, 1996
=#
function jump_SemiRand(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:4])
    @variable(Lower(model), y[i=1:2])

    @objective(Upper(model), Min,
        -4x[1] +8x[2] +x[3] -x[4] + 9y[1] -9y[2])
    @constraint(Upper(model), [i=1:4], x[i] >= 0)
    @constraint(Upper(model), -9x[1] + 3x[2] -8x[3] + 3x[4] +3y[1]        <= -1)
    @constraint(Upper(model), +4x[1] -10x[2] +3x[3] + 5x[4] +8y[1] +8y[2] <= 25)
    @constraint(Upper(model), +4x[1] - 2x[2] -2x[3] +10x[4] -5y[1] +8y[2] <= 21)
    @constraint(Upper(model), +9x[1] - 9x[2] +4x[3] - 3x[4] - y[1] -9y[2] <= -1)
    @constraint(Upper(model), -2x[1] - 2x[2] +8x[3] - 5x[4] +5y[1] +8y[2] <= 20)
    @constraint(Upper(model), +7x[1] + 2x[2] -5x[3] + 4x[4] -5y[1]        <= 11)

    @constraint(Upper(model), [i=1:2], y[i] >= 0)

    @objective(Lower(model), Min,
        -9y[1] + 9y[2])
    @constraint(Lower(model), [i=1:2], y[i] >= 0)
    @constraint(Lower(model), -6x[1] + x[2] + x[3] - 3x[4] -9y[1] -7y[2] <= -15)
    @constraint(Lower(model),        +4x[2] +5x[3] +10x[4]               <= 26)
    @constraint(Lower(model), -9x[1] +9x[2] -9x[3] + 5x[4] -5y[1] -4y[2] <= -5)
    @constraint(Lower(model), +5x[1] +3x[2] + x[3] + 9x[4] + y[1] +5y[2] <= 32)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value.(x) ≈ [1.57504, 0.83668, 0.188807, 2.17092] atol=1e-4
    @test value.(y) ≈ [1.88765, 0.0] atol=1e-4
end

#=
    Decomposition Techniques in Mathematical Programming
    A. Conejo et al.
    Springer
=#

# Chapter 7.2, pag 281
function jump_DTMP_01(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y)

    @objective(Upper(model), Min, -x + 4y)
    @constraint(Upper(model), y + 2x <= 8)

    @objective(Lower(model), Min, -x - y)
    @constraint(Lower(model), -y <= 0)
    @constraint(Lower(model), x + y <= 7)
    @constraint(Lower(model), -x <= 0)
    @constraint(Lower(model),  x <= 4)

    optimize!(model)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 6 atol=atol
end

#=
    modification to test variables used in a single level
=#

function jump_DTMP_01_mod1(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(UpperOnly(model), z)
    @variable(Lower(model), y)
    @variable(LowerOnly(model), w)

    @objective(Upper(model), Min, -x + 4y + z)
    @constraint(Upper(model), y + 2x + z <= 9)
    @constraint(Upper(model), z == 1)

    @objective(Lower(model), Min, -x - y + w)
    @constraint(Lower(model),  y >= 0)
    @constraint(Lower(model), x + y + w <= 8)
    @constraint(Lower(model),  x >= 0)
    @constraint(Lower(model),  x <= 4)
    @constraint(Lower(model),  w == 1)

    optimize!(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 6 atol=atol
    @test value(z) ≈ 1 atol=atol
    @test value(w) ≈ 1 atol=atol
end

function jump_DTMP_01_mod2_error(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(UpperOnly(model), z)
    @variable(Lower(model), y)
    @variable(LowerOnly(model), w)

    @test_throws ErrorException @objective(Upper(model), Min, -x + 4y + z + w)
    @constraint(Upper(model), y + 2x + z <= 9)
    @constraint(Upper(model), z == 1)

    @test_throws ErrorException @objective(Lower(model), Min, -x - y + z)
end

jump_DTMP_01_mod3vec(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_DTMP_01_mod3(optimizer, true, mode, config)
jump_DTMP_01_mod3(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) = _jump_DTMP_01_mod3(optimizer, false, mode, config)
function _jump_DTMP_01_mod3(optimizer, vectorized::Bool, mode, config)

    atol = config.atol

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    if vectorized
        @variables(Upper(model), begin
            x
            z
        end)
    else
        @variable(Upper(model), x)
        @variable(Upper(model), z)
    end
    if vectorized
        @variables(Lower(model), begin
            y
            w
        end)
    else
        @variable(Lower(model), y)
        @variable(Lower(model), w)
    end

    @objective(Upper(model), Min, -x + 4y + z)
    if vectorized
        @constraints(Upper(model), begin
            y + 2x + z <= 9
            z == 1
        end)
    else
        @constraint(Upper(model), y + 2x + z <= 9)
        @constraint(Upper(model), z == 1)
    end

    @objective(Lower(model), Min, -x - y + w)
    if vectorized
        @constraints(Lower(model), begin
            y >= 0
            x + y + w <= 8
            x >= 0
            x <= 4
            w == 1
        end)
    else
        @constraint(Lower(model),  y >= 0)
        @constraint(Lower(model), x + y + w <= 8)
        @constraint(Lower(model),  x >= 0)
        @constraint(Lower(model),  x <= 4)
        @constraint(Lower(model),  w == 1)
    end

    optimize!(model)

    @test value(x) ≈ 1 atol=atol
    @test value(y) ≈ 6 atol=atol
    @test value(z) ≈ 1 atol=atol
    @test value(w) ≈ 1 atol=atol
end

#=
    Conejo, A. J., Baringo, L., Kazempour, S. J., and Siddiqui, A. S. (2016).
    Investment in Electricity Generation and Transmission. Springer.
=#
function jump_conejo2016(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config(); bounds = false)

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x, start = 50)
    if bounds
        @variable(Lower(model), -1 <= y[i=1:3] <= 300, start = [50, 150, 0][i])
    else
        @variable(Lower(model), y[i=1:3], start = [50, 150, 0][i])
    end
    # 2 and 3 are lower only

    @constraint(Upper(model), lb0, x >= 0)
    @constraint(Upper(model), ub0, x <= 250)

    @objective(Lower(model), Min, 10y[1] + 12y[2] + 15y[3])
    @constraint(Lower(model), b, y[1] + y[2] + y[3] == 200)
    @constraint(Lower(model), ub1, y[1] <= x)
    @constraint(Lower(model), ub2, y[2] <= 150)
    @constraint(Lower(model), ub3, y[3] <= 100)
    @constraint(Lower(model), lb[i=1:3], y[i] >= 0)
    if bounds
        @variable(Upper(model), 0 <= lambda <= 20, DualOf(b), start = 15)
    else
        @variable(Upper(model), lambda, DualOf(b), start = 15)
    end

    @objective(Upper(model), Min, 40_000x + 8760*(10y[1]-lambda*y[1]))

    optimize!(model)

    primal_status(model)
    termination_status(model)

    @test objective_value(model) ≈ -190_000 atol=1e-1 rtol=1e-2
    @test value(x) ≈ 50 atol=1e-3 rtol=1e-2
    @test value.(y) ≈ [50, 150, 0] atol=1e-3 rtol=1e-2
    @test value(lambda) ≈ 15 atol=1e-3 rtol=1e-2
end

#=
    Bruno Fanzeres PhD thesis Robust Strategic Bidding in Auction-Based Markets.
=#
function jump_fanzeres2017(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), q1, start = 20)
    @variable(Lower(model), g[i=1:4], start = [20, 40, 40, 0][i])

    @constraint(Upper(model), q1 >= 0)
    @constraint(Upper(model), q1 <= 100)

    @objective(Lower(model), Min, 50g[2] + 100g[3] + 1000g[4])
    @constraint(Lower(model), b, g[1] + g[2] + g[3] + g[4] == 100)
    @constraint(Lower(model), g[1] <= q1)
    @constraint(Lower(model), g[2] <= 40)
    @constraint(Lower(model), g[3] <= 40)
    @constraint(Lower(model), g[4] <= 100)
    @constraint(Lower(model), lb[i=1:4], g[i] >= 0)

    @variable(Upper(model), lambda, DualOf(b), start = 1_000)

    @objective(Upper(model), Max, lambda*g[1])

    optimize!(model)

    primal_status(model)
    termination_status(model)

    @test objective_value(model) ≈ 20_000  atol=1e-1
    @test BilevelJuMP.lower_objective_value(model) ≈ 6_000  atol=1e-1
    @test value(q1) ≈ 20 atol=1e-3
    @test value.(g) ≈ [20, 40, 40, 0] atol=1e-3
    @test value(lambda) ≈ 1_000 atol=1e-3
end

function jump_16(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config())

    atol = config.atol
    start = config.start_value

    MOI.empty!(optimizer)
    m = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(m), 0 <= x <= 3 )
    @variable(Lower(m), 0 <= y[1:2])
    @variable(Lower(m), 0 <= z[1:2])

    a = JuMP.Containers.DenseAxisArray([y[i]+ z[i] for i = 1:2], 1:2)
    d = JuMP.Containers.SparseAxisArray(Dict((1, i) => y[i] for i in 1:2))

    c = y[1] + y[2]

    @constraint(Upper(m), x >=1)

    @constraint(Lower(m), y[1] >= 0.5)
    @constraint(Lower(m), y[2] >= 0.5)

    @objective(Upper(m), Min, x)

    @objective(Lower(m), Min, y[1]+z[2])

    optimize!(m)

    # TODO add tests
    _a = value.(a)
    _c = value.(c)
    _d = value.(d)

end

function jump_conic01(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config(); bounds = false)

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x[i=1:3])
    @variable(Lower(model), y[i=1:3])

    MOI.set(optimizer, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

    @constraint(Upper(model), soc_up, x in SecondOrderCone())
    @constraint(Lower(model), soc_lw, y in SecondOrderCone())
    if bounds
        BilevelJuMP.set_dual_upper_bound(soc_lw, +[5., 5., 5.])
        BilevelJuMP.set_dual_lower_bound(soc_lw, -[5., 5., 5.])
        # bounds defined in the upper level are not dualized
        for i in 1:3
            @constraint(Upper(model), y[i] in MOI.LessThan(+5.0))
            @constraint(Upper(model), y[i] in MOI.GreaterThan(-5.0))
        end
    end

    @objective(Upper(model), Min, x[1])
    @objective(Lower(model), Min, y[1])

    optimize!(model)

    primal_status(model)
    termination_status(model)

    @test objective_value(model) ≈ 0  atol=1e-1
    @test value.(x) ≈ [0, 0, 0] atol=1e-3
    @test value.(y) ≈ [0, 0, 0] atol=1e-3
end

#=
    The models of bilevel programming with lower level second-order cone programs
    Chi et al. Journal of Inequalities and Applications 2014, 2014:168
    https://journalofinequalitiesandapplications.springeropen.com/articles/10.1186/1029-242X-2014-168
    https://core.ac.uk/download/pdf/81261904.pdf
=#

function jump_conic02(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config(); bounds = false)

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y[i=1:2])

    @objective(Upper(model), Min, x + 3(y[1] -y[2]))
    if bounds
        @constraint(Upper(model), x in MOI.LessThan(+6.0))
        @constraint(Upper(model), x in MOI.GreaterThan(+2.0))
    else
        @constraint(Upper(model), x >= 2)
        @constraint(Upper(model), x <= 6)
    end

    @objective(Lower(model), Min, - (y[1] - y[2]))
    @constraint(Lower(model), lb_y_1, y[1] >= 0)
    @constraint(Lower(model), lb_y_2, y[2] >= 0)
    @constraint(Lower(model), con1, x +  (y[1] - y[2]) <=  8)
    @constraint(Lower(model), con2, x + 4(y[1] - y[2]) >=  8)
    @constraint(Lower(model), con3, x + 2(y[1] - y[2]) <= 12)
    @constraint(Lower(model), soc_lw, y in SecondOrderCone())

    if bounds
        BilevelJuMP.set_dual_upper_bound(soc_lw, +[5., 5.])
        BilevelJuMP.set_dual_lower_bound(soc_lw, -[5., 5.])
        # require lower bounds
        for con in [con1, con3]
            BilevelJuMP.set_dual_lower_bound(con, -15)
        end
        # require upper bounds
        for con in [lb_y_1, lb_y_2, con2]
            BilevelJuMP.set_dual_upper_bound(con, +15)
        end
        # bounds defined in the upper level are not dualized
        for i in 1:2
            @constraint(Upper(model), y[i] in MOI.LessThan(+5.0))
            @constraint(Upper(model), y[i] in MOI.GreaterThan(-5.0))
        end
    end

    optimize!(model)

    primal_status(model)
    termination_status(model)
    value.(y)

    @test objective_value(model) ≈ 12  atol=1e-1
    @test value(x) ≈ 6 atol=1e-3
    @test value(y[2]) >= 0 - 1e-3
    @test value(y[1]) - value(y[2]) ≈ 2 atol=1e-3
end
function jump_conic03(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config(); bounds = false)

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y[i=1:2])

    @objective(Upper(model), Min, x + 2(y[1] + y[2]))
    @constraint(Upper(model), y[1] + y[2] <= 3) # creates disconnected region
    if bounds
        @constraint(Upper(model), x in MOI.LessThan(+6.0))
        @constraint(Upper(model), x in MOI.GreaterThan(+0.0))
    else
        @constraint(Upper(model), x >= 0)
        @constraint(Upper(model), x <= 6)
    end

    @objective(Lower(model), Min, - y[1] - y[2])
    @constraint(Lower(model), con1, y[1] >= 0)
    @constraint(Lower(model), con2, y[2] <= 0)
    @constraint(Lower(model), con3, x +  (y[1] + y[2]) <= 8)
    @constraint(Lower(model), con4, x + 3(y[1] + y[2]) >= 8)
    @constraint(Lower(model), con5,-x +  (y[1] + y[2]) <= 0)
    @constraint(Lower(model), soc_lw, y in SecondOrderCone())

    if bounds
        BilevelJuMP.set_dual_upper_bound(soc_lw, +[5., 5.])
        BilevelJuMP.set_dual_lower_bound(soc_lw, -[5., 5.])
        # require lower bounds
        for con in [con2, con3, con5]
            BilevelJuMP.set_dual_lower_bound(con, -15)
        end
        # require upper bounds
        for con in [con1, con4]
            BilevelJuMP.set_dual_upper_bound(con, +15)
        end
        # bounds defined in the upper level are not dualized
        for i in 1:2
            @constraint(Upper(model), y[i] in MOI.LessThan(+5.0))
            @constraint(Upper(model), y[i] in MOI.GreaterThan(-5.0))
        end
    end

    optimize!(model)

    JuMP.raw_status(model)
    primal_status(model)
    termination_status(model)
    value.(y)

    @test objective_value(model) ≈ 6  atol=1e-1
    @test value(x) ≈ 2 atol=1e-3
    @test value(y[2]) <= 0 + 1e-3
    @test value(y[1]) + value(y[2]) ≈ 2 atol=1e-3
end
function jump_conic04(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config(); bounds = false)

    MOI.empty!(optimizer)
    model = BilevelModel(()->optimizer, mode = mode)

    @variable(Upper(model), x)
    @variable(Lower(model), y[i=1:3])

    @objective(Upper(model), Min, x + 3y[1])
    if bounds
        @constraint(Upper(model), x in MOI.GreaterThan(+2.0))
        @constraint(Upper(model), x in MOI.LessThan(+6.0))
    else
        @constraint(Upper(model), x >= 2)
        @constraint(Upper(model), x <= 6)
    end

    @objective(Lower(model), Min, - y[1])
    @constraint(Lower(model), a, x +  y[1] <=  8)
    @constraint(Lower(model), b, x + 4y[1] >=  8)
    @constraint(Lower(model), c, x + 2y[1] <= 12)
    # @constraint(Lower(model), d, y[1] >= 0)
    @constraint(Lower(model), soc_lw, y in SecondOrderCone())

    if bounds
        BilevelJuMP.set_dual_upper_bound(soc_lw, +[5., 5., 5.])
        BilevelJuMP.set_dual_lower_bound(soc_lw, -[5., 5., 5.])
        # require lower bounds
        for con in [a, c]
            BilevelJuMP.set_dual_lower_bound(con, -15)
        end
        # require upper bounds
        for con in [b]
            BilevelJuMP.set_dual_upper_bound(con, +15)
        end
        # bounds defined in the upper level are not dualized
        for i in 1:3
            @constraint(Upper(model), y[i] in MOI.LessThan(+5.0))
            @constraint(Upper(model), y[i] in MOI.GreaterThan(-5.0))
        end
    end

    optimize!(model)

    JuMP.raw_status(model)
    primal_status(model)
    termination_status(model)
    value.(y)

    @test objective_value(model) ≈ 12 atol=1e-1
    @test value(x) ≈ 6 atol=1e-3
    @test value(y[1]) ≈ 2 atol=1e-3
    @test sqrt(value(y[2])^2 + value(y[3])^2) <= 2 + 1e-3
end


text:
# ======================================
# Kronecker powers
# ======================================

using LinearAlgebra: checksquare

"""
    kron_pow(x::IA.Interval, pow::Int)

Given an interval x and an integer pow, compute `x^pow`.

### Input

- `x`   -- interval
- `pow` -- integer

### Output

An interval enclosure of `x^pow`.

### Examples

```julia
julia> [kron_pow(2 .. 3, i) for i in 1:3]
3-element Array{IntervalArithmetic.Interval{Float64},1}:
  [2, 3]
  [4, 9]
 [8, 27]
```
"""
function kron_pow(x::IA.Interval, pow::Int)
    if pow == 1
        return x
    end

    return x^pow
end

"""
    kron_pow(x::Interval, pow::Int)

Given an interval x and an integer pow, compute `x^pow`.

### Input

- `x`   -- interval
- `pow` -- integer

### Output

An interval enclosure of `x^pow` as a LazySets `Interval`.
"""
function kron_pow(x::Interval, pow::Int)
    if pow == 1
        return x
    end

    return Interval(kron_pow(x.dat, pow))
end

"""
    kron_pow(H::AbstractHyperrectangle, pow::Int)

Given hyperrectangular set `H` and an integer `pow`, compute the Kronecker power
`H^{⊗ pow}`.

### Input

- `H`   -- hyperrectangular set
- `pow` -- integer power

### Output

A hyperrectangle.

### Algorithm

We compute `H^{⊗ pow}` where `H` is a hyperrectangular set by working with `H`
as a product of intervals.

See also `kron_pow` which requires `DynamicPolynomials.jl`.
"""
function kron_pow(H::AbstractHyperrectangle, pow::Int, algorithm=nothing)
    if pow == 1
        return H
    end

    if isnothing(algorithm)
        if isdefined(@__MODULE__, :DynamicPolynomials)
            algorithm = "symbolic"
        else
            algorithm = "explicit"
        end
    end
    if algorithm == "explicit"
        return _kron_pow_explicit(H, pow)
    elseif algorithm == "symbolic"
        return _kron_pow_symbolic(H, pow)
    else
        throw(ArgumentError("algorithm \"$algorithm\" is not known"))
    end
end

function _kron_pow_explicit(H::AbstractHyperrectangle, pow::Int)
    x = [Interval(low(H, j), high(H, j)) for j in 1:dim(H)]
    r = reduce(kron, fill(x, pow))

    low_r = min.(r)
    high_r = max.(r)
    return Hyperrectangle(low=low_r, high=high_r)
end

"""
    kron_pow_stack(x::IA.Interval, pow::Int)

Return a hyperrectangle with the interval powers `[x, x^2, …, x^pow]`.

### Input

- `x`   -- interval
- `pow` -- integer power

### Output

A hyperrectangle such that the `i`-th dimension is the interval `x^i`.
"""
function kron_pow_stack(x::IA.Interval, pow::Int)
    return convert(Hyperrectangle, IntervalBox([kron_pow(x, i) for i in 1:pow]))
end

"""
    kron_pow_stack(x::Interval, pow::Int)

Return a hyperrectangle with the interval powers `[x, x^2, …, x^pow]`.

### Input

- `x`   -- interval
- `pow` -- integer power

### Output

A hyperrectangle such that the `i`-th dimension is the interval `x^i`.
"""
function kron_pow_stack(x::Interval, pow::Int)
    return convert(Hyperrectangle, CartesianProductArray([kron_pow(x, i) for i in 1:pow]))
end

"""
    kron_pow_stack(H::AbstractHyperrectangle, pow::Int)

Return the Cartesian product array ``H × H^{⊗2} × ⋯ × H^{⊗pow}`` where ``H`` is a
hyperrectangular set and ``H^{⊗ i}`` is the `i`-th Kronecker power of ``H``.

### Input

- `H`   -- hyperrectangular set
- `pow` -- integer power

### Output

A Cartesian product array of hyperrectangles.
"""
function kron_pow_stack(H::AbstractHyperrectangle, pow::Int)
    out = [kron_pow(H, i) for i in 1:pow]
    return CartesianProductArray(out)
end

# ------------------------------------------------------------
# Functionality that requires DynamicPolynomials.jl
# ------------------------------------------------------------

function load_kron_dynamicpolynomials()
return quote

function _kron_pow_symbolic(H::AbstractHyperrectangle{N}, pow::Int) where {N}
    n = dim(H)
    @polyvar x[1:n]
    B = convert(IntervalBox, H)
    dict = Dict((x[i] => i for i in 1:n)...)
    y = kron_pow(x, pow)

    out = Vector{IA.Interval{N}}(undef, length(y))
    for (i, p) in enumerate(y)
        aux = interval(1)
        for (xj, j) in powers(p)
            aux = aux * B[dict[xj]]^j
        end
        out[i] = aux
    end
    Bpow = IntervalBox(out)
    Hpow = convert(Hyperrectangle, Bpow)
    return Hpow
end

end end  # quote / load_kron_dynamicpolynomials()

# ------------------------------------------------------------
# Functionality that requires MultivariatePolynomials.jl
# ------------------------------------------------------------

function load_kron_multivariate()
return quote

"""
    kron_pow(x::Vector{<:AbstractVariable}, pow::Int)

Compute the higher order concrete Kronecker power: `x ⊗ x ⊗ ... ⊗ x`, `pow` times
for a vector of symbolic monomials.

### Input

- `x`   -- polynomial variable
- `pow` -- integer

### Output

Vector of multivariate monomial corresponding to `x^{⊗ pow}`.

### Examples

```julia
julia> using DynamicPolynomials

julia> @polyvar x[1:2]
(PolyVar{true}[x₁, x₂],)

julia> x
2-element Array{PolyVar{true},1}:
 x₁
 x₂

julia> kron_pow(x, 2)
4-element Array{Monomial{true},1}:
 x₁²
 x₁x₂
 x₁x₂
 x₂²
```
"""
function kron_pow(x::Vector{<:AbstractVariable}, pow::Int)
    @assert pow > 0 "expected positive power, got $pow"
    if pow == 1
        return x
    else
        return kron(x, kron_pow(x, pow-1))
    end
end

end end  # quote / load_kron_multivariate()

ext: jl


text:
using GraphUtilities, GenericTensorNetworks, Graphs
using Test

@testset "json solve" begin
    for property in ["SizeMax", "SizeMin",
            "SizeMax3", "SizeMin3", "CountingAll",
            "CountingMax", "CountingMin", "CountingMax3", "CountingMin3", "GraphPolynomial",
            "SingleConfigMax", "SingleConfigMin", "SingleConfigMax3", "SingleConfigMin3",
            "ConfigsMaxTree", "ConfigsMin",
            "ConfigsMaxTree3", "ConfigsMin3",
            "ConfigsAll", "ConfigsAllTree"
        ]
        res = GraphUtilities.application(Dict(
            "api"=>"solve",
            "api.solve"=>Dict(
                "graph"=>Dict("nv"=>10, "edges"=>[[e.src, e.dst] for e in edges(smallgraph(:petersen))]),
                "problem"=>"MaximalIS",
                "property"=>property,
                "openvertices"=>[],
                "fixedvertices"=>Dict(),
                "optimizer"=>Dict(
                    "method"=>"TreeSA",
                    "TreeSA"=>Dict(
                        "sc_target"=>20,
                        "sc_weight"=>1.0,
                        "ntrials"=>1,
                        "niters"=>5,
                        "openvertices"=>[],
                        "fixedvertices"=>Dict(),
                        "nslices"=>0,
                        "rw_weight"=>1.0,
                        "betas"=>collect(0.01:0.1:30)
                        )
                    ),
                "cudadevice"=>-1
            )
        ))
        @test res isa Union{Dict, Vector}
    end
    res = GraphUtilities.application(Dict(
        "api"=>"solve",
        "api.solve"=>Dict(
            "graph"=>Dict("nv"=>10, "edges"=>[[e.src, e.dst] for e in edges(smallgraph(:petersen))]),
            "problem"=>"MaximalIS",
            "property"=>"SizeMin",
        )
    ))
    @test res == Dict("size"=>3.0)
    @test_throws GraphUtilities.VerificationError GraphUtilities.application(Dict(
        "api"=>"solve",
        "api.solve"=>Dict(
        "graph"=>Dict("nv"=>10, "edges"=>[[e.src, e.dst] for e in edges(smallgraph(:petersen))]),
        "problem"=>"MaximalIS",
        "property"=>"SizeMi",)
    ))
end

@testset "json graph" begin
    d_graph = Dict(
        "api"=>"graph",
        "api.graph" => Dict(
            "type"=>"kings",
            "type.kings"=> Dict(
                "m"=>8,
                "n"=>8,
                "filling"=>0.8,
                "seed"=>2
            )
        )
    )
    res = GraphUtilities.application(d_graph)
    @test res isa Dict
end

@testset "json opteinsum" begin
    d_opteinsum = Dict("api"=>"opteinsum",
        "api.opteinsum" => Dict(
            "inputs"=>[[1,2], [2,3], [3,4], [5,4]],
            "output"=>[1,6],
            "method"=>"TreeSA",
            "sizes"=>[1=>2, 2=>2, 3=>2,
                4=>2, 5=>2]
        )
    )
    res = GraphUtilities.application(d_opteinsum)
    @test res isa Dict
    @test GraphUtilities.application(Dict("api"=>"help")) isa Dict
end


text:
using WriteVTK
using Test

function main()
    # Suppose that the mesh is made of 5 points:
    cells = [
        MeshCell(VTKCellTypes.VTK_TRIANGLE, [1, 4, 2]),
        MeshCell(VTKCellTypes.VTK_QUAD, [2, 4, 3, 5]),
    ]
    x = [0.95, 0.16, 0.07, 0.21, 0.60]
    y = [0.32, 0.55, 0.87, 0.12, 0.85]
    vtufile = "simulation/simulation_1.vtu"
    rm(vtufile, force = true)
    @time outfiles = pvtk_grid("simulation", x, y, cells; part = 1, nparts = 1) do pvtk # 2D
        pvtk["Pressure"] = x
        pvtk["Processor"] = [1,2]
        pvtk["Temperature", VTKPointData()] = y
        pvtk["Id", VTKCellData()] = [2, 1]
    end
    @test isfile(vtufile)
    @test vtufile ∈ outfiles
    println("Saved:  ", join(outfiles, "  "))
    @test WriteVTK._serial_filename(3, 100, "prefix", ".ext") == "prefix_003.ext"
    outfiles
end

main()


text:
function add_disjunction(m::Model,disj...;reformulation=:BMR,M=missing,eps=1e-6,name=missing)
    @assert m isa Model "A valid JuMP Model must be provided."
    @assert reformulation in [:BMR, :CHR] "Invalid reformulation method passed to keyword argument `:reformulation`. Valid options are :BMR (Big-M Reformulation) and :CHR (Convex-Hull Reormulation)."

    #create binary indicator variables for each disjunction
    disj_name = ismissing(name) ? Symbol("disj",gensym()) : name
    @assert isnothing(variable_by_name(m,string(disj_name))) "Name for disjunction binary cannot be the same as an existing variable."
    eval(:(@variable($m, $disj_name[i = 1:length($disj)], Bin)))
    #enforce exclussive OR
    eval(:(@constraint($m,sum($disj_name[i] for i = 1:length($disj)) == 1)))
    #apply reformulation
    if reformulation == :BMR
        reformulate(m, disj, disj_name, reformulation, M)
    elseif reformulation == :CHR
        reformulate(m, disj, disj_name, reformulation, eps)
    end
end


text:
using Coverage: Coverage
using CoverageBase
using Test

@testset "locate tests" begin
    @test isdir(CoverageBase.julia_top())
    tn = testnames()
    @test !isempty(tn)
    @test "goto" in tn
end

@testset "locate source" begin
    thisfile = pathof(CoverageBase)
    @test fixpath(thisfile) == thisfile
    @test fixpath(joinpath("something", "relative.jl")) == joinpath("base", "something", "relative.jl")
    files = [Coverage.FileCoverage(f[2], "", Coverage.CovCount[]) for f in Base._included_files]
    @test all(isabspath(c.filename) for c in files) || files
    fixpath!(files)
    push!(files, Coverage.FileCoverage("test/goto.jl", "", Coverage.CovCount[]))
    for c in files
        @test !isabspath(c.filename) || c
        @test c.filename == "test/goto.jl" || startswith(c.filename, "base") || startswith(c.filename, "stdlib") || c
        lpath = CoverageBase.fixabspath(c.filename)
        @test isabspath(lpath)
        @test isfile(lpath)
        @test fixpath(lpath) == c.filename
    end
    files[1].source = "fakesource"
    push!(files, Coverage.FileCoverage(thisfile, "", Coverage.CovCount[]))
    push!(files, Coverage.FileCoverage(thisfile * "_nonexistant", "", Coverage.CovCount[]))
    readsource!(files)
    @test files[1].source == "fakesource"
    @test pop!(files).source == ""
    @test pop!(files).source == read(thisfile, String)
    @test all(!isempty(c.source) for c in files) || files
end

@testset "run tests" begin
    @test runtests(["goto", "Test"])
    @test_warn "ERROR: failed process" @test !runtests(["fail_does_not_even_exist", "goto"])
end


text:
using PosteriorAnalysis
using Base.Test
using MacroTools

import PosteriorAnalysis: _transform_expr

######################################################################
# test transformations (internal)
######################################################################

"Test that expression is passed through."
macro test_passthrough(expr)
    (expr2, captured_names) = _transform_expr(expr)
    @test isequal(expr, expr2)
    @test isempty(captured_names)
end

test_template(f, expr) = f(_transform_expr(expr)...)

macro test_template(expr, template, variables...)
    evar = gensym(:expr)
    cvar = gensym(:captured_names)
    quote
        test_template($(Meta.quot(expr))) do $(evar), $(cvar)
            @match $(evar) begin
                $template => begin
                    $([:(@test get($cvar, $(Meta.quot(v)), nothing) == $v)
                       for v in variables]...)
                end
                _ => throw(error("Result $($evar) does not match template"))
            end
        end
    end
end

@testset "macros" begin

    @test_throws ErrorException _transform_expr(:(@v())) # no arguments
    @test_throws ErrorException _transform_expr(:(@v a b)) # too many arguments

    @test_passthrough(:(a+b))
    @test_passthrough(:(a.b))
    @test_passthrough(:(quote a+@v(b)+c end))

    @test_template @v(a) a_ a
    @test_template @v(a)+b a_+b a
    @test_template @v(a).b a_.b a
    @test_template a+@v(b)+@v(c) a+b_+c_ b c

end

######################################################################
# test operations
######################################################################

@testset "operations" begin

    a = 1:3
    b = 4:6
    pd_z = PosteriorDraws(length(a))
    pd_a = PosteriorDraws(:a => 1:3)
    pd_ab = PosteriorDraws([:b,:a], [b,a]) # constructor will enforce sorting

    @test size(pd_ab) == (2, length(a))
    @test keys(pd_ab) == [:a,:b]

    @test addvars(pd_a, :b => b) == pd_ab
    @test addvars(pd_z, :a => a, :b => b) == pd_ab
    @test vcat(pd_z, pd_ab) == pd_ab
    @test dropvars(pd_ab, :b) == pd_a

    @test all(pd_ab[:a] .== a)
    @test all(pd_ab[:a,2:3] .== a[2:3])

    pd_abc0 = addvars(pd_ab, :c => a+b)
    pd_abc1 = map(pd_ab, :c, +, :a, :b)
    pd_abc2 = @pdmap pd_ab c => @v(a) + @v(b)

    @test pd_abc0 == pd_abc1
    @test pd_abc0 == pd_abc2

    io = IOBuffer()
    @test (show(io, pd_ab); true)
    @test takebuf_string(io) ==
        """PosteriorDraws with 3 observations
            a => Int64
            b => Int64
        """
end

@testset "elttype test" begin

    f(a,b) = Float64(a+b)
    g(a,b) = [string(a),string(b)]

    pd1 = PosteriorDraws(:a => 1:3, :b => 11:13)

    @test pd1[:a] == collect(1:3)
    @test pd1[:b] == collect(11:13)
    @test eltype(pd1[:a]) ≡ eltype(pd1[:b]) ≡ Int

    pd2 = @pdmap pd1 c => f(@v(a), @v(b))
    pd3 = @pdmap pd2 d => g(@v(a), @v(b))

    @test eltype(pd2[:c]) ≡ eltype(pd3[:c]) ≡ Float64

    @test eltype(pd3[:d]) ≡ Array{String, 1}

end


text:
import RCall: rcopy,
              RObject,
              LangSxp,
              @rput,
              @rget,
              @R_str

"""
    fulldummy2dummy(lhs, grp, bar="|")

"""
function _fulldummy2dummy(lhs, grp, bar="|")
    terms = Vector{String}(undef, length(lhs.terms))

    for (idx, tt) in enumerate(lhs.terms)
        if tt isa CategoricalTerm &&
            tt.contrasts isa StatsModels.ContrastsMatrix{StatsModels.FullDummyCoding}
            terms[idx] = "dummy($(tt.sym), levels($(tt.sym)))"
        else
            terms[idx] = string(tt)
        end

    end

    "( " * foldl((x,y) -> "$x + $y", terms) * " $bar $grp)"

end

function convert_julia_to_r(f::StatsModels.FormulaTerm)::AbstractString
    rhs = Vector{String}(undef, length(f.rhs))

    for (idx, trm) in enumerate(f.rhs)

        if trm isa MixedModels.ZeroCorr
            grp = trm.term.rhs.sym
            terms = Vector{String}(undef, length(trm.term.lhs.terms))

            if LMER == "afex::lmer_alt"
                rhs[idx] = _fulldummy2dummy(trm.term.lhs, trm.term.rhs.sym, "||")

            else
                terms_assemble = true
                for (idx,tt) in enumerate(trm.term.lhs.terms)
                    if tt isa CategoricalTerm && size(tt.contrasts.matrix, 2) != 2
                        AFEX_INSTALLED ||
                                throw(ArgumentError("zerocorr for factors with more than 2 levels not supported without using afex::lmer_alt"))

                        @info "Swapping to afex::lmer_alt instead of lme4::g/lmer from here on out..."
                        global GLMER = global LMER = "afex::lmer_alt"
                        rhs[idx] = _fulldummy2dummy(trm.term.lhs, trm.term.rhs.sym, "||")
                        terms_assemble = false
                        break
                    elseif tt isa InterceptTerm
                        terms[idx] = "( 1 | $grp)"
                    elseif tt isa CategoricalTerm && tt.contrasts isa StatsModels.ContrastsMatrix{StatsModels.FullDummyCoding}
                        terms[idx] = "( 0 + dummy($(tt.sym), levels($(tt.sym))) | $grp)"
                    else
                        terms[idx] = "( 0 + $(tt.sym) | $grp)"
                    end

                end

                if terms_assemble
                    rhs[idx] = foldl((x,y) -> "$x + $y", terms)
                end
            end

        elseif trm isa MixedModels.RandomEffectsTerm
            rhs[idx] = _fulldummy2dummy(trm.lhs, trm.rhs.sym)
        else
            rhs[idx] = string(trm)
        end

    end

    formula = "$(f.lhs) ~ $(foldl((x,y) -> "$x + $y", string.(rhs)))"
    # fulldummy in FE
    formula = replace(formula, r"fulldummy\(([^)]*)\)" => s"dummy(\1,levels(\1))")

    formula = replace(formula, ":(log" => "(log")
    formula = replace(formula, ":(exp" => "(exp")
    # should we consider supporting division on the lhs?
    # what about nesting syntax on the rhs?
    # zscore won't work here because the operations within a formula are broadcast
    # and never see the whole array
    # formula = replace(formula, ":(zscore" => "(scale")

    occursin(":(", formula) && throw(ArgumentError("Formula contains a transformation."))

    formula = replace(formula, "&" => ":")
end

# NOTE: this is known to yet to work properly!
function convert_r_to_julia(f::RObject{LangSxp})::StatsModels.FormulaTerm
    _, lhs, rhs = rcopy(R"as.character($f)")
    cbind =  match(r"cbind\((.*),(.*)\)", lhs);
    if isnothing(cbind)
        # we do weird things because of the way / is translated
        vv = rcopy(Expr, f)
    else
        # we don't need to compute the weights from the cbind formula because we can
        # get them directly from the glmerMod object
        numerator = successes = first(cbind.captures);
        failures = last(cbind.captures);
        denominator = "( $successes + $failures )";
        lhs = "$(numerator) / $(denominator)";
        form = "$lhs ~ $rhs"
        # we do weird things because of the way / is translated
        vv = rcopy(Expr, R"""as.formula($form)""")

    end
    vv = eval(:@formula($vv))
    # @info vv
    # @info typeof(vv)
    return vv
end

convert_r_to_julia(f::Ptr{LangSxp})::StatsModels.FormulaTerm =
                convert_r_to_julia(RObject{LangSxp}(f))


text:
export MeterProvider

using Dates: time

const N_MAX_METRICS = 1_000

Base.@kwdef mutable struct Metric{A<:AbstractAggregation}
    name::String
    description::Union{String,Nothing}
    attribute_keys::Union{Tuple{Vararg{String}},Nothing}
    aggregation::A
    criteria::Criteria
end

function (metric::Metric)(ms)
    for m in ms
        metric(m)
    end
end

Base.iterate(m::Metric, args...) = iterate(m.aggregation, args...)
Base.getindex(m::Metric, k) = getindex(m.aggregation, k)
Base.length(m::Metric) = length(m.aggregation)

function (metric::Metric)(m::Measurement)
    if isnothing(metric.attribute_keys)
        filtered_attributes = StaticAttrs()
    else
        interested_keys = keys(m.attributes) ∩ metric.attribute_keys
        filtered_keys = setdiff(keys(m.attributes), interested_keys)
        m = Measurement(m.value, m.attributes[Tuple(interested_keys)])
        filtered_attributes = m.attributes[Tuple(filtered_keys)]
    end

    span_ctx = span_context()
    if isnothing(span_ctx)
        trace_id = INVALID_TRACE_ID
        span_id = INVALID_SPAN_ID
    else
        trace_id = span_ctx.trace_id
        span_id = span_ctx.span_id
    end

    exemplar = Exemplar(;
        value = m,
        time_unix_nano = UInt(time() * 10^9),
        filtered_attributes = filtered_attributes,
        trace_id = trace_id,
        span_id = span_id,
    )
    metric.aggregation(exemplar)
end

struct MeterProvider <: AbstractMeterProvider
    resource::Resource
    meters::IdDict{Meter,Nothing}
    instrument_associated_metric_names::IdDict{AbstractInstrument,Set{String}}
    async_instruments::IdDict{AbstractAsyncInstrument,Nothing}
    views::Vector{View}
    metrics::Dict{String,Metric}
    n_max_metrics::UInt
end

function MeterProvider(;
    resource = Resource(),
    views = View[],
    n_max_metrics = N_MAX_METRICS,
)
    if isempty(views)
        push!(views, View(; instrument_name = "*"))
    end

    MeterProvider(
        resource,
        Dict{String,Meter}(),
        IdDict{AbstractInstrument,Vector{String}}(),
        IdDict{AbstractAsyncInstrument,Nothing}(),
        views,
        Dict{String,Metric}(),
        n_max_metrics,
    )
end

function Base.push!(p::MeterProvider, m::Meter)
    p.meters[m] = nothing
    for ins in m.instruments
        push!(p, ins)
    end
end

function Base.push!(p::MeterProvider, ins::AbstractInstrument)
    if !haskey(p.meters, ins.meter)
        p.meters[ins.meter] = nothing
    end

    if !haskey(p.instrument_associated_metric_names, ins)
        drop_views = (v for v in p.views if v.aggregation === DROP)
        valid_views = (v for v in p.views if v.aggregation !== DROP)

        is_drop = false
        for v in drop_views
            if occursin(ins, v.criteria)
                is_drop = true
                break
            end
        end

        if is_drop
            @info "Metric won't be created since the view for the instrument [$(ins.name)] is configured to DROP."
        else
            for v in valid_views
                if occursin(ins, v.criteria)
                    if length(p.metrics) >= p.n_max_metrics
                        @warn "Maximum number of metrics [$(p.n_max_metrics)] reached. Instrument [$(ins.name)] related metrics are dropped!"
                    else
                        metric_name = something(v.name, ins.name)
                        if haskey(p.metrics, metric_name)
                            @info "Found a duplicate metric [$metric_name]. The original one will be reused."
                        else
                            p.metrics[metric_name] = Metric(
                                name = metric_name,
                                description = something(v.description, ins.description),
                                attribute_keys = v.attribute_keys,
                                aggregation = if isnothing(v.aggregation)
                                    default_aggregation(ins)
                                else
                                    v.aggregation
                                end,
                                criteria = v.criteria,
                            )
                        end
                        push!(
                            get!(p.instrument_associated_metric_names, ins, Set{String}()),
                            metric_name,
                        )

                        if ins isa AbstractAsyncInstrument
                            p.async_instruments[ins] = nothing
                        end
                    end
                end
            end
        end
    end
end

function Base.push!(p::MeterProvider, (instrument, measurement))
    if haskey(p.instrument_associated_metric_names, instrument)
        for metric_name in p.instrument_associated_metric_names[instrument]
            metric = p.metrics[metric_name]
            metric(measurement)
        end
    else
        @debug "Instrument [$(instrument.name)] is not registered in the meter provider."
    end
end


text:
import Base: exp, expm1, log, log10, log1p, sqrt, cbrt, exponent,
             significand, sin, sinpi, cos, cospi, tan, sec, cot, csc,
             sinh, cosh, tanh, coth, sech, csch,
             asin, acos, atan, acot, asec, acsc,
             asinh, acosh, atanh, acoth, asech, acsch, sinc, cosc,
             +, -, %, &, *

import Base: broadcast
import Base: Broadcast
import Base.Broadcast: Broadcasted, BroadcastStyle, combine_eltypes

"""
This is a way of suggesting that stage should call
stage_operand with the operation and other arguments
"""
struct PromotePartition{T,N} <: ArrayOp{T,N}
    data::AbstractArray{T,N}
end

size(p::PromotePartition) = size(domain(p.data))

struct BCast{B, T, Nd} <: ArrayOp{T, Nd}
    bcasted::B
end

BCast(b::Broadcasted) = BCast{typeof(b), combine_eltypes(b.f, b.args), length(axes(b))}(b)

size(x::BCast) = map(length, axes(x.bcasted))

function stage_operands(ctx, ::BCast, xs::ArrayOp...)
    map(x->cached_stage(ctx, x), xs)
end

function stage_operands(ctx, ::BCast, x::ArrayOp, y::PromotePartition)
    stg_x = cached_stage(ctx, x)
    y1 = Distribute(domain(stg_x), y.data)
    stg_x, cached_stage(ctx, y1)
end

function stage_operands(ctx, ::BCast, x::PromotePartition, y::ArrayOp)
    stg_y = cached_stage(ctx, y)
    x1 = Distribute(domain(stg_y), x.data)
    cached_stage(ctx, x1), stg_y
end

struct DaggerBroadcastStyle <: BroadcastStyle end

BroadcastStyle(::Type{<:ArrayOp}) = DaggerBroadcastStyle()
BroadcastStyle(::DaggerBroadcastStyle, ::BroadcastStyle) = DaggerBroadcastStyle()
BroadcastStyle(::BroadcastStyle, ::DaggerBroadcastStyle) = DaggerBroadcastStyle()

function Base.copy(b::Broadcast.Broadcasted{<:DaggerBroadcastStyle})
    BCast(b)
end

function stage(ctx, node::BCast)
    bc = Broadcast.flatten(node.bcasted)
    args = bc.args
    args1 = map(args) do x
        x isa ArrayOp ? cached_stage(ctx, x) : x
    end
    ds = map(x->x isa DArray ? domainchunks(x) : nothing, args1)
    sz = size(node)
    dss = filter(x->x !== nothing, collect(ds))
    cumlengths = ntuple(ndims(node)) do i
        idx = findfirst(d -> i <= length(d.cumlength), dss)
        if idx === nothing
            [sz[i]] # just one slice
        end
        dss[idx].cumlength[i]
    end

    args2 = map(args1) do arg
        if arg isa AbstractArray
            s = size(arg)
            splits = map(enumerate(s)) do dim
                i, n = dim
                if n == 1
                    return [1]
                else
                    cumlengths[i]
                end
            end |> Tuple
            dmn = DomainBlocks(ntuple(_->1, length(s)), splits)
            cached_stage(ctx, Distribute(dmn, arg)).chunks
        else
            arg
        end
    end
    blcks = DomainBlocks(map(_->1, size(node)), cumlengths)

    thunks = broadcast(delayed((args...)->broadcast(bc.f, args...); ),
                       args2...)
    DArray(eltype(node), domain(node), blcks, thunks)
end

export mappart, mapchunk

struct MapChunk{F, Ni, T, Nd} <: ArrayOp{T, Nd}
    f::F
    input::NTuple{Ni, ArrayOp{T,Nd}}
end

mapchunk(f, xs::ArrayOp...) = MapChunk(f, xs)
Base.@deprecate mappart(args...) mapchunk(args...)
function stage(ctx, node::MapChunk)
    inputs = map(x->cached_stage(ctx, x), node.input)
    thunks = map(map(chunks, inputs)...) do ps...
        Thunk(node.f, ps...)
    end

    DArray(Any, domain(inputs[1]), domainchunks(inputs[1]), thunks)
end