text:
function jump_01(optimizer, mode = BilevelJuMP.SOS1Mode())

    # 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

    model = BilevelModel()

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

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

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

    @constraint(Lower(model), 2x+y <= 4)
    @constraint(Lower(model), x+2y <= 4)
    @constraint(Lower(model), x >= 0)
    @constraint(Lower(model), y >= 0)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ -8

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

end

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

    model = BilevelModel()

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

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

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

    @constraint(Lower(model), x+y <= 8)
    @constraint(Lower(model), x+4y >= 8)
    @constraint(Lower(model), x+2y <= 13)
    @constraint(Lower(model), x >= 1)
    @constraint(Lower(model), x <= 6)
    @constraint(Lower(model), y >= 0)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 12

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

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
function jump_03(optimizer, mode = BilevelJuMP.SOS1Mode())

    model = BilevelModel()

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

    @objective(Upper(model), Min, 3x + y)
    @constraint(Upper(model), x <= 5)
    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 3* (3.5*8/15) + (8/15)

    @test value(x) ≈ 3.5*8/15
    @test value(y) ≈ 8/15

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 3* (3.5*8/15) + (8/15)

    @test value(x) ≈ 3.5*8/15
    @test value(y) ≈ 8/15

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())

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ 0

    @test value(x[1]) ≈ 1/100
    @test value(x[2]) ≈ 0
    @test value(y) ≈ 0

end

# sec 3.5.2.2 pag 44 -> product

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

    model = BilevelModel()

    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])

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ -1

    # 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())

    model = BilevelModel()

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

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

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

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

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

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

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

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    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())

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

end

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    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
function jump_06(optimizer, mode = BilevelJuMP.SOS1Mode())

    model = BilevelModel()

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

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

    @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
    @constraint(Lower(model), y >= 0)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test objective_value(model) ≈ -12

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

end

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

end

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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())

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

    model = BilevelModel()

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

    @objective(Upper(model), Min, -10.1x + 10y[1] - y[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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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())

    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]

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    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())

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

end

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

    for a in [0 0.1 0.2]
        model = BilevelModel()

        @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, optimizer, mode)

        primal_status(model)

        termination_status(model)

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

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 1
    @test value(y) ≈ 75
    @test value(Z) ≈ 21+2/3
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())
    # send y to upper level

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

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

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

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

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test value(x) ≈ 4
    @test value(y) ≈ 4
end

# 9.2.4 - parg 211
function jump_HTP_lin03(optimizer, mode = BilevelJuMP.SOS1Mode())
    # send y to upper level

    model = BilevelModel()

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

    @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)

    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]

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

    # TODO fix broadcast
    # @constraint(Lower(model), H1*y + H2*x .== b)
    @constraint(Lower(model), [i=1:6], y[i] >= 0)

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

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

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

    @objective(Upper(model), Min,
        -8x[1] -4x[2] + 4y[1] - 40y[2] + 4y[3])
        # the sign of 4y[3] seems off, but slutionis 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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    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]
    @test value.(y) ≈ [1.5, 0]
end

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    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]
    @test value.(y) ≈ [1.5, 0]
end

# TODO - add quadratic problems

#=
    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())

    ###
    ### PART 1
    ###

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

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

    ###
    ### PART 2
    ###

    model = BilevelModel()

    @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)
    @constraint(Lower(model), y >= +x)

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

    optimize!(model, optimizer, mode)

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

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

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

    ###
    ### PART 1
    ###

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

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

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

    ###
    ### PART 2
    ###

    model = BilevelModel()

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

    @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)

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

    optimize!(model, optimizer, mode)

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

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

    @test value.(x) ≈ [0, 4, 0, 15, 9.2, 0]
    @test value.(y) ≈ [0, 0, 2]
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())

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    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())

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

    primal_status(model)

    termination_status(model)

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

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

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

    model = BilevelModel()

    @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)

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

    optimize!(model, optimizer, mode)

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

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

    model = BilevelModel()

    @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


text:
module day18

include("inputs.jl")
ex1 = inputs.ex1_d18()
ex2 = inputs.ex2_d18()
ex3 = inputs.ex3_d18()
ex4 = inputs.ex4_d18()
ex5 = inputs.ex5_d18()
ex6 = inputs.ex6_d18()
input = inputs.input_d18()

# Part 1
"""Use the power of Julia to define a new binary operator ⨦ (plussim)
this operator will perform multiplcation, but be parsed as having
the same precendence as addition."""

⨦(a,b) = a * b

proc1(in) = eval(Meta.parse(replace(in, "*"=>"⨦")))

@assert proc1(ex1) == 71
@assert proc1(ex2) == 51
@assert proc1(ex3) == 26
@assert proc1(ex4) == 437
@assert proc1(ex5) == 12240
@assert proc1(ex6) == 13632

part1(input) = sum([proc1(line) for line in input])
@info "Day 18, Part 1 answer: $(part1(input))"

# Part 2
⨱(a,b) = a + b

proc2(in) = eval(Meta.parse(replace(replace(in, "*"=>"⨦"), "+"=>"⨱")))

@assert proc2(ex1) == 231
@assert proc2(ex2) == 51
@assert proc2(ex3) == 46
@assert proc2(ex4) == 1445
@assert proc2(ex5) == 669060
@assert proc2(ex6) == 23340

part2(input) = sum([proc2(line) for line in input])
@info "Day 18, Part 2 answer: $(part2(input))"

end

ext: jl


text:
"""
    characterize(x::AbstractInteractor)

Characterize the state of a single interactor. Calls the `characterize`
method on the replicator carried by the interactor.
"""
function characterize(x::AbstractInteractor)
  characterize(x.replicator)
end


text:
# This source file is though to put that code required in order to
# customize ConformingFESpaces.jl to H(div)-conforming global FE Spaces built
# out of RaviartThomas FEs. In particular, this customization is in the
# definition of the shape functions (get_cell_shapefuns) and the DoFs
# (get_cell_dof_basis) of the **global** FE space, which requires a sign flip
# for those sitting on facets of the slave cell of the facet.

# Two key ingredients in the implementation of this type of ReferenceFE are the
# get_cell_shapefuns(model,cell_reffes,::Conformity) and
# get_cell_dof_basis(mode,cell_reffes,::Conformity) overloads.
# These are written such that, for each cell K, they return the shape functions
# and dof values in the *global* RT space. For a dof owned by a face which is shared by
# two cells, there is a master and a slave cell. The slave cell first computes the
# shape functions and dof values using local-to-cell data structures, but then flips the
# sign of both in order to get their corresponding counterparts in the **global**
# RT space. As as result we have the following:

# * When we interpolate a function into the global FE space, and we perform the cell-wise
#   DoF values to global DoF values gather operation, we can either extract the global DoF value
#   from the master or slave cell without worrying about the sign.
# * When we evaluate a global FE function, and we perform the global DoF values to
#   cell-wise DoF values scatter operation, we don't have to worry about the sign either.
#   On the slave cell, we will have both the sign of the DoF value, and the sign of the
#   shape function corresponding to the global DoF.
# * We do NOT have to use the signed determinant, but its absolute value, in the Piola Map.

struct TransformRTDofBasis{Dc,Dp} <: Map end ;

function get_cell_dof_basis(model::DiscreteModel,
                            cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}},
                            ::DivConformity)
    sign_flip = get_sign_flip(model, cell_reffe)
    cell_map  = get_cell_map(Triangulation(model))
    phi       = cell_map[1]
    Jt        = lazy_map(Broadcasting(∇),cell_map)
    x         = lazy_map(get_nodes,lazy_map(get_dof_basis,cell_reffe))
    Jtx       = lazy_map(evaluate,Jt,x)
    reffe     = cell_reffe[1]
    Dc        = num_dims(reffe)
    et        = return_type(get_prebasis(reffe))
    pt        = Point{Dc,et}
    Dp        = first(size(return_type(phi,zero(pt))))
    k         = TransformRTDofBasis{Dc,Dp}()
    lazy_map(k,cell_reffe,Jtx,sign_flip)
end

function get_cell_shapefuns(model::DiscreteModel,
                            cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}},
                            ::DivConformity)
    sign_flip = get_sign_flip(model, cell_reffe)
    cell_reffe_shapefuns=lazy_map(get_shapefuns,cell_reffe)
    k=ContraVariantPiolaMap()
    lazy_map(k,
             cell_reffe_shapefuns,
             get_cell_map(Triangulation(model)),
             lazy_map(Broadcasting(constant_field), sign_flip))
end

struct SignFlipMap{T} <: Map
  model::T
end

function return_cache(k::SignFlipMap,reffe,cell_id)
  model = k.model
  D = num_cell_dims(model)
  gtopo = get_grid_topology(model)

  # Extract composition among cells and facets
  cell_wise_facets_ids = get_faces(gtopo, D, D - 1)
  cache_cell_wise_facets_ids = array_cache(cell_wise_facets_ids)

  # Extract cells around facets
  cells_around_facets = get_faces(gtopo, D - 1, D)
  cache_cells_around_facets = array_cache(cells_around_facets)

  (cell_wise_facets_ids,
   cache_cell_wise_facets_ids,
   cells_around_facets,
   cache_cells_around_facets,
   CachedVector(Bool))

end

function evaluate!(cache,k::SignFlipMap,reffe,cell_id)
  model = k.model

  cell_wise_facets_ids,
  cache_cell_wise_facets_ids,
  cells_around_facets,
  cache_cells_around_facets,
  sign_flip_cached = cache

  setsize!(sign_flip_cached, (num_dofs(reffe),))
  sign_flip = sign_flip_cached.array
  sign_flip .= false

  D = num_dims(reffe)
  face_own_dofs = get_face_own_dofs(reffe)
  facet_lid = get_offsets(get_polytope(reffe))[D] + 1
  cell_facets_ids = getindex!(cache_cell_wise_facets_ids,
                              cell_wise_facets_ids,
                              cell_id)
  for facet_gid in cell_facets_ids
      facet_cells_around = getindex!(cache_cells_around_facets,
           cells_around_facets,
           facet_gid)
      is_slave = (findfirst((x) -> (x == cell_id), facet_cells_around) == 2)
      if is_slave
          for dof in face_own_dofs[facet_lid]
              sign_flip[dof] = true
          end
      end
      facet_lid = facet_lid + 1
  end
  sign_flip
end

function get_sign_flip(model::DiscreteModel,
                       cell_reffe::AbstractArray{<:GenericRefFE{RaviartThomas}})
    lazy_map(SignFlipMap(model),
            cell_reffe,
            get_cell_to_bgcell(model))
end

function return_cache(::TransformRTDofBasis{Dc,Dp},
                      reffe::GenericRefFE{RaviartThomas},
                      Jtx,
                      ::AbstractVector{Bool}) where {Dc,Dp}
  p = get_polytope(reffe)
  prebasis = get_prebasis(reffe)
  order = get_order(prebasis)
  et = return_type(prebasis)
  dofs = get_dof_basis(reffe)
  nodes, nf_nodes, nf_moments =  get_nodes(dofs),
                                 get_face_nodes_dofs(dofs),
                                 get_face_moments(dofs)
  db = MomentBasedDofBasis(nodes,nf_moments,nf_nodes)
  face_moments = [ similar(i,VectorValue{Dp,et})  for i in nf_moments ]

  cache = (db.nodes, db.face_nodes, nf_moments, face_moments)
  cache
end

function evaluate!(cache,
                   ::TransformRTDofBasis,
                   reffe::GenericRefFE{RaviartThomas},
                   Jt_q,
                   sign_flip::AbstractVector{Bool})
  nodes, nf_nodes, nf_moments, face_moments = cache
  face_own_dofs=get_face_own_dofs(reffe)
  for face in 1:length(face_moments)
    nf_moments_face   = nf_moments[face]
    face_moments_face = face_moments[face]
    if length(nf_moments_face) > 0
      sign = (-1)^sign_flip[face_own_dofs[face][1]]
      num_qpoints, num_moments = size(nf_moments_face)
      for i in 1:num_qpoints
        Jt_q_i = Jt_q[nf_nodes[face][i]]
        change = sign * meas(Jt_q_i) * pinvJt(Jt_q_i)
        for j in 1:num_moments
          face_moments_face[i,j] = change ⋅ nf_moments_face[i,j]
        end
      end
    end
  end
  MomentBasedDofBasis(nodes,face_moments,nf_nodes)
end

# Support for DIV operator
function DIV(f::LazyArray{<:Fill{T}}) where T
  df=DIV(f.args[1])
  k=f.maps.value
  lazy_map(k,df)
end
function DIV(f::LazyArray{<:Fill{Broadcasting{Operation{ContraVariantPiolaMap}}}})
  ϕrgₖ       = f.args[1]
  fsign_flip = f.args[4]
  div_ϕrgₖ = lazy_map(Broadcasting(divergence),ϕrgₖ)
  fsign_flip=lazy_map(Broadcasting(Operation(x->(-1)^x)), fsign_flip)
  lazy_map(Broadcasting(Operation(*)),fsign_flip,div_ϕrgₖ)
end
function DIV(a::LazyArray{<:Fill{typeof(linear_combination)}})
  i_to_basis = DIV(a.args[2])
  i_to_values = a.args[1]
  lazy_map(linear_combination,i_to_values,i_to_basis)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b36fb690-3447-11eb-3629-edbddae4ebd1
input = readlines("./input.txt");

# ╔═╡ 4ab2f908-3449-11eb-299d-53d5a0a3044e
"""
	Converts Input to a Int matrix
"""
function parseInput(input::Array{String,1})
	rows = length(input)
	cols = length(input[1])
	boolMat = Matrix{Int}(undef,rows,cols)
	for row ∈ 1:rows
		for col ∈ 1:cols
			boolMat[row,col] = input[row][col] == '#' ? 1 : 0
		end
	end
	return boolMat
end

# ╔═╡ e1af2152-3835-11eb-1940-4744bf549a6f

begin
	"""
 	Gets the tree count for all given slopes
	"""
	function GetTreeCountForSlope((right,down)::Tuple{Int,Int})
		function(boolMat::Matrix{Int})
			r,c = size(boolMat)
			x,y = 1,1
			count = 0
			while x < r
				count += boolMat[x,y]
				x=(x+down)
				y=mod1(right+y,c)
			end
			return count
		end
	end

	function GetTreeCountForSlope(slopes::Array{Tuple{Int64,Int64},1})
		function(boolMat::Matrix{Int})
			counts = Array{Int}(undef,size(slopes))
			for i ∈ 1:size(slopes)[1]
				(right,down) = slopes[i]
				r,c = size(boolMat)
				x,y = 1,1
				count = 0
				while x < r
					count += boolMat[x,y]
					x=(x+down)
					y=mod1(right+y,c)
				end
				counts[i] = count
			end
			return counts
		end
	end

end

# ╔═╡ 7b8df328-3838-11eb-269d-21098ee6f5d5
parsedInput = input |> parseInput;

# ╔═╡ d839bbcc-3836-11eb-2fb0-a50603b97bcd
Part1 = parsedInput |> GetTreeCountForSlope((3,1))

# ╔═╡ 86e320c0-383a-11eb-3160-e3da943b134d
Part2 = reduce(*,parsedInput |> GetTreeCountForSlope([(1,1),(3,1),(5,1),(7,1),(1,2)]) )

# ╔═╡ Cell order:
# ╠═b36fb690-3447-11eb-3629-edbddae4ebd1
# ╠═4ab2f908-3449-11eb-299d-53d5a0a3044e
# ╠═e1af2152-3835-11eb-1940-4744bf549a6f
# ╠═7b8df328-3838-11eb-269d-21098ee6f5d5
# ╠═d839bbcc-3836-11eb-2fb0-a50603b97bcd
# ╠═86e320c0-383a-11eb-3160-e3da943b134d


text:
using GithubMarkdown
using Test

cd(@__DIR__) do
@testset "GithubMarkdown.jl" begin
    @testset "in-memory" begin
        ioin = IOBuffer("foo")
        ioout = IOBuffer()
        expectedout = "
foo
\n"
        rendergfm(ioout, ioin)
        @test String(take!(ioout)) == expectedout
    end

    @testset "html" begin
        io = IOBuffer()
        rendergfm(io, "test1.md"; format = "html")
        @test String(take!(io)) == read("test1.html", String)

        rendergfm("testout.html", joinpath(@__DIR__, "test1.md"); format = "html")
        @test read("testout.html", String) == read("test1.html", String)
    end

    @testset "latex" begin
        io = IOBuffer()
        rendergfm(io, "test1.md"; format = "latex")
        @test String(take!(io)) == read("test1.tex", String)
    end

    @testset "html input" begin
        io = IOBuffer()
        rendergfm(io, "malicious.md", extensions = [])
        @test String(take!(io)) == "\n
asd
\n"
        rendergfm(io, "malicious.md")
        @test String(take!(io)) == "<script>console.log(\"this is bad.\")</script>\n
asd
\n"
        rendergfm(io, "malicious.md", removehtml = true)
        @test String(take!(io)) == "\n\n"
    end

    @testset "error handling" begin
        @test_throws ErrorException rendergfm("dontwriteme", "idontexist.md")
        @test !isfile("dontwriteme")
        @test_throws ArgumentError rendergfm("dontwriteme", "test1.md"; format = "nahhh")
        @test !isfile("dontwriteme")
        @test_throws ArgumentError rendergfm("dontwriteme", "test1.md"; extensions = ["nope"])
        @test !isfile("dontwriteme")
    end

    @testset "string inputs" begin
        @test strip(rendergfm("a **b** c")) == "
a b c
"
    end
end
end


text:
const t=2.1

@inline initial_condition_generator(rng,n,model) = (X=[0.0],φ=[1])

const tau = fixed_time(t)

const x = range(0,10,length=10_001)


text:
##This is like the toffset_fits function but doesn't evaluate all toffsets:
##does a recursive binary-ish search.  If you want to eliminate a higher or lower fraction of points per-iteration you can adjust elim_frac higher or lower.
##Also differsf from toffset_fits in that it only returns the tfm and mm of the best fit
#function bsearch_toffsets(target, testimgs, toffsets; elim_frac = 0.4, sigmas=(1.0,1.0))
#    @show ntoffsets = length(toffsets)
#    ntrials = convert(Int, size(testimgs,3) / ntoffsets)
#    @show ntrials
#    firstmov = squeeze(mean(Float64.(view(testimgs, :, :, 1:ntrials)), 3),3)
#    firsttfm, firstmm = align2d(target, firstmov; sigmas=sigmas)
#    lastmov = squeeze(mean(Float64.(view(testimgs, :, :, (ntoffsets-1)*ntrials+1:ntoffsets*ntrials)), 3),3)
#    lasttfm, lastmm = align2d(target, lastmov; sigmas=sigmas)
#	if firstmm <= lastmm
#		toffsets_new = toffsets[1: round(Int, ntoffsets*(1-elim_frac))]
#		@show ntoffsets_new =length(toffsets_new)
#		n_elim = ntoffsets - ntoffsets_new
#        @show toffsets_new
#		if ntoffsets_new == 1 || toffsets == toffsets_new
#            return first(toffsets), firsttfm, firstmm
#		else
#    		testimgs = Float64.(view(testimgs, :, :, 1:ntoffsets_new*ntrials))
#            print("here1\n")
#			return bsearch_toffsets(target, testimgs, toffsets_new; elim_frac = elim_frac, sigmas = sigmas)
#		end
#	else
#        toffsets_new = toffsets[(round(Int, ntoffsets*elim_frac)+1):ntoffsets]
#		@show ntoffsets_new = length(toffsets_new)
#		n_elim = ntoffsets - ntoffsets_new
#        @show toffsets_new
#		if ntoffsets_new == 1 || toffsets == toffsets_new
#            return last(toffsets), lasttfm, lastmm
#		else
#            @show n_elim
#            @show ntoffsets
#            @show size(testimgs)
#    		testimgs = Float64.(view(testimgs, :, :, n_elim*ntrials+1:ntoffsets*ntrials))
#            print("here2\n")
#            @show size(testimgs)
#            @show length(toffsets_new)
#			return bsearch_toffsets(target, testimgs, toffsets_new; elim_frac = elim_frac, sigmas = sigmas)
#		end
#	end
#end
#


text:
"""
Recursive function to print the type tree for a given `type`, including its
supertype and all its subtypes.
"""
function show_type_tree(type, indentation_level=0)
    if type === Any
        throw(ArgumentError("Any is subtype and supertype of Any"))
    end
    if indentation_level == 0
        println(supertype(type))
        indentation_level += 1
    end
    println("    " ^ indentation_level, type)
    for subtype in subtypes(type)
            show_type_tree(subtype, indentation_level + 1)
    end
end

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
include("1.jl")
function mark_krossxxx(r::Robot)
    for side in (Nord,Sud,Ost,West)
        putmarkholexxx(r,side)
        backingxxx(r, inverse(side))
    end
    putmarker!(r)
end

function putmarkholexxx(r::Robot,side::HorizonSide)
    while (!isborder(r,side) && !isborder(r,HorizonSide(mod(Int(side)+1, 4))))
        move!(r,side)
        move!(r,HorizonSide(mod(Int(side)+1, 4)))
        putmarker!(r)
    end
end

function backingxxx(r::Robot,side::HorizonSide)
    while ismarker(r)
        move!(r,side)
        move!(r,HorizonSide(mod(Int(side)+1, 4)))
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
export agent_step!, model_step!, households, municipalities

households(model::ABM) = Iterators.filter(a -> isa(a, Household), allagents(model))
households(municipality::Municipality, model::ABM) = Iterators.filter(
    a -> isa(a, Household) && a.municipality == municipality.id,
    allagents(model),
)
municipalities(model::ABM) = Iterators.filter(a -> isa(a, Municipality), allagents(model))

"""
    active_interventions(municipality, year)

Returns all active interventions in the planner. Due to the `year = -1` -> always active
convention we must merge the current year's plan with the `-1` key (if extant).
"""
active_interventions(m::Municipality, year::Int) = vcat(
    get(m.interventions, -1, Intervention[]),
    get(m.interventions, year, Intervention[]),
)

function agent_step!(house::Household, model) # yearly
    # Once a year households may update their oss
    if !house.oss && house.information
        #Agent has an understanding of compliance but has not yet upgraded.
        if rand(model.rng) < house.compliance
            house.oss = true
            municipality = model[house.municipality]
            if municipality.houseowner_type isa Social
                # Tell neighbours to increase their compliance
                neighbors = nearby_ids(house, model, municipality.neighbor_distance)
                # We must filter out bordering municipalities first
                for nid in Iterators.filter(
                    n ->
                        model[n] isa Household &&
                            model[n].municipality == house.municipality,
                    neighbors,
                )
                    model[nid].compliance = min(model[nid].compliance * 1.5, 0.99)
                end
            end
        else
            house.implementation_lag += 1
        end
    end
end

function agent_step!(municipality::Municipality, model)
    # Update pollution level: nutrient inflow by affectors
    municipality.affectors = count(a -> !a.oss, households(municipality, model))
    # Municipality decides on sewage treatment rule
    if municipality.regulate
        # Check whether update of rules is needed.
        # When information from monitoring indicates need to act/rule, the legislation status is updated
        if municipality.information && !municipality.legislation
            municipality.legislation = true
            #TODO: Not useful in a multi municipality context
            #model.outcomes_year_when_informed = model.year - 1 # Netlogo comment: remember this as t1 for later calculation of implementation efficiency
        end

        #TODO: This check now essentially supersedes `municipality.regulate`
        any(
            i -> i isa WastewaterTreatment,
            active_interventions(municipality, model.year),
        ) && water_treatment!(model, municipality)
    end
end

function model_step!(model)
    # Run the decision optimiser only if this is not a test model.
    # This check avoids an recursion based stack overflow.
    haskey(model.properties, :test) || set_policy!(model)

    household_log!(model)
    # Update bream stock, pike stock, vegetation and nutrients (daily)
    OrdinaryDiffEq.step!(model.lake, 365.0, true)
    # Nutrients affect lake dynamics
    nutrient_load!(model, model.nutrient_series)
    OrdinaryDiffEq.u_modified!(model.lake, true)

    # Intervention which necessitate lake-wide changes
    aggregate_regulate!(model)
    model.year += 1

    # Monitoring (in this order to prevent immediate action based on new monitoring data)
    monitor!(model)

    # remember the year when the desired state is restored
    # look only for this year after degradation and regulation of system has started
    #TODO: This makes no sense in the context of multiple municipalities
    # if any(m->m.legislation, municipalities(model)) &&
    #     model.outcomes_year_when_desired_pike_is_back == 0
    #     restoration_log!(model)
    # end
end

function set_policy!(model::ABM)
    apply_knowledge!(model)
    if model.policy.start == model.year || (
        model.year > model.policy.start &&
        mod(model.year - model.policy.start, model.policy.every) == 0
    )
        make_decision!(model)
    end
end

function apply_knowledge!(model::ABM)
    for knowledge in model.knowledge
        apply_knowledge!(model, knowledge)
    end
end

function apply_knowledge!(model, knowledge::VegetationImbalance)
    # We know that it's appropriate to plant in deep turbid states, but
    # also that bream flips before everything else. This causes a vegetation
    # imbalance that puts the lake into a dramatically high vegetation state
    # that cannot be recovered by the optimiser. In addition, planting is no
    # longer the cheapest or most efficient option in this case, so we swap
    # all planned Planting interventions to Trawling.

    knowledge.year_target_reached > 0 && return nothing

    #TODO: This assumes `VegetationImbalance` only applies to Turbid->Clear
    # which is not explicitly true.
    if clear_state_region(model, model.year)
        # We have reached an appropriate state, stop interventions
        knowledge.year_target_reached = model.year
        for municipality in municipalities(model)
            for plans in
                values(filter(i -> i.first >= model.year, municipality.interventions))
                for (idx, p) in enumerate(plans)
                    deleteat!(plans, idx)
                end
            end
        end
    else
        B, P, V = model.lake.u
        if B < knowledge.bream_density_flip
            for municipality in municipalities(model)
                for plans in
                    values(filter(i -> i.first >= model.year, municipality.interventions))
                    for (idx, p) in enumerate(plans)
                        if p isa Planting
                            replace = Trawling(p.rate, p.cost)
                            deleteat!(plans, idx)
                            push!(plans, replace)
                        end
                    end
                end
            end
        end
    end
end

function household_log!(model::ABM)
    # checking the success of upgrade
    if model.outcomes_year_when_informed > 0 &&
       model.outcomes_year_of_full_upgrade <= model.year
        model.outcomes_upgraded_households_sum += count(a -> a.oss, households(model)) #NOTE: This is some cumulative trickery in Netlogo, and is not a true count of upgraded households.
    end
    # remember the year when all households were upgraded
    if model.outcomes_year_of_full_upgrade == 0 && all(a -> a.oss, households(model))
        model.outcomes_year_of_full_upgrade = model.year
    end
end

function monitor!(model::ABM)
    loss = model.pike_expectation - model.lake.u[2] #pike
    if loss > 0 && model.outcomes_year_when_pike_became_critical == 0
        model.outcomes_year_when_pike_became_critical = model.year
    end
    if model.outcomes_year_when_nutrients_became_critical == 0 &&
       model.lake.p.nutrients > model.critical_nutrients
        model.outcomes_year_when_nutrients_became_critical = model.year
    end
    for municipality in municipalities(model)
        if !municipality.information && threshold_monitor(loss, municipality, model)
            municipality.information = true
            if municipality.agents_uniform
                # Immediate legislation when lake state is found to be undesirable
                municipality.legislation = true
                #TODO: Not useful in a multi municipality context
                #TODO: Should be brought into a municipality based governance strategy
                #model.outcomes_year_when_informed = model.year
            end #Otherwise, it takes a year to act.
        end
    end
end

function threshold_monitor(loss::Float64, municipality::Municipality, model::ABM)
    if municipality.threshold_variable isa Nutrients
        return model.lake.p.nutrients > model.critical_nutrients
    elseif municipality.threshold_variable isa Pike
        trigger = municipality.respond_direct ? 0 : rand(model.rng)
        return trigger < pike_loss_perception(loss)
    end
    false
end

function pike_loss_perception(loss::Float64)
    #TODO: This is the partial defiancy of pike loss, look it up, seems pretty arbitrary currently.
    loss > 0 ? loss^2 / (loss^2 + 0.75^2) : 0
end

#TODO: This makes no sense in the context of multiple municipalities
#function restoration_log!(model::ABM)
#    if model.threshold_variable isa Nutrients
#        if model.lake.p.nutrients < model.critical_nutrients
#            model.outcomes_year_when_desired_level_is_back = model.year
#        end
#    elseif model.threshold_variable isa Pike
#        if model.lake.pike > model.pike_expectation
#            model.outcomes_year_when_desired_pike_is_back = model.year
#        end
#    end
#end

function aggregate_regulate!(model::ABM)
    # Wastewater is internal to the municipality and therefore does not need to update rates here.
    schedule = Iterators.flatten(
        active_interventions(m, model.year) for m in municipalities(model)
    )

    planting = Iterators.filter(i -> i isa Planting, schedule)
    model.lake.p.pv = 0.0
    for project in planting
        model.lake.p.pv += project.rate
    end

    trawling = Iterators.filter(i -> i isa Trawling, schedule)
    model.lake.p.tb = 0.0
    for project in trawling
        model.lake.p.tb += project.rate
    end

    angling = Iterators.filter(i -> i isa Angling, schedule)
    model.lake.p.mp = model.init_pike_mortality
    for project in angling
        model.lake.p.mp += project.rate
    end

    OrdinaryDiffEq.u_modified!(model.lake, true)
end

function water_treatment!(model::ABM, municipality::Municipality)
    #Inform households, eventually enforcing.
    if municipality.legislation && municipality.affectors > 0
        if count(a -> a.oss, households(municipality, model)) == 0
            #TODO: This assumes perfect information transfer, and only happens once.
            # Delays will happen here and we should structure this as an information campaign in the future
            for house in Iterators.filter(
                a -> !a.oss && !a.information,
                households(municipality, model),
            )
                house.information = true
                house.implementation_lag = 0
            end
        end
        # Enforcement after 5 years
        if municipality.houseowner_type isa Enforced
            for house in Iterators.filter(
                a -> a.implementation_lag > 4,
                households(municipality, model),
            )
                house.compliance = min(house.compliance * 1.5, 0.99)
            end
        end
    end
end


text:
# Utils for both simulation and measuring of nets.
module SpinGlass

using DataFrames

export Net, metropolis!, magnetization, energy, χsg, measure

#
#
# Constructor methods
#
#

type Net
    Jup     :: Array{Int8,3}
    Jright  :: Array{Int8,3}
    Jfront  :: Array{Int8,3}
    σ       :: Array{Int8,3}
    L       :: Int64
end

function Net(L::Int64; kind=:spinglass)
    if kind == :spinglass
    Net(  rand([-1,1],L,L,L) # Jup
        , rand([-1,1],L,L,L) # Jright
        , rand([-1,1],L,L,L) # Jfront
        , rand([-1,1],L,L,L) # σ
        , L)                 # L
    elseif kind == :ferromagnet
    Net(  ones(Int8,L,L,L)   # Jup
        , ones(Int8,L,L,L)   # Jright
        , ones(Int8,L,L,L)   # Jfront
        , rand([-1,1],L,L,L) # σ
        , L)                 # L
    elseif kind == :antiferromagnet
    Net(  -ones(Int8,L,L,L)  # Jup
        , -ones(Int8,L,L,L)  # Jright
        , -ones(Int8,L,L,L)  # Jfront
        , rand([-1,1],L,L,L) # σ
        , L)                 # L
    else
        error("Unsuported initializer kind \"$kind\"")
    end
end

#
#
# Evolution
#
#

# Contribution to the hamiltonian of a spin.
function H(sample,x,y,z)
    L = sample.L
    left   = sample.Jright[mod1(x-1,L),y,z] * sample.σ[mod1(x-1,L), y, z]
    right  = sample.Jright[x,y,z]           * sample.σ[mod1(x+1,L), y, z]
    up     =    sample.Jup[x,y,z]           * sample.σ[x, mod1(y+1,L), z]
    down   =    sample.Jup[x,mod1(y-1,L),z] * sample.σ[x, mod1(y-1,L), z]
    front  = sample.Jfront[x,y,z]           * sample.σ[x, y, mod1(z+1,L)]
    behind = sample.Jfront[x,y,mod1(z-1,L)] * sample.σ[x, y, mod1(z-1,L)]

    neighbours_term = left + right + up + down + behind + front

    return -sample.σ[x,y,z] * neighbours_term
end

function metropolis!(sample,β)
    L = sample.L
    for z in 1:L, y in 1:L, x in 1:L
        # Get ΔE and probability of the tentative flip.
        ΔE = -2*H(sample,x,y,z)
        P = exp(-β*ΔE)
        if P > rand()
            sample.σ[x,y,z] *= -1
        end
    end
end

#
#
# Observables
#
#

magnetization(sample::Net)     = mean(sample.σ)
magnetization(spins::BitArray) = 2*(mean(spins)-0.5)

function energy(sample)
    L = sample.L
    acc = 0
    for x in 1:L, y in 1:L, z in 1:L
        acc += H(sample,x,y,z)
    end
    return acc / (L*L*L) / 3 / 2
end

function χsg(sample)
    L = sample.L
    SiSj = 0
    Si = 0
    for x in 1:L, y in 1:L, z in 1:L
        right  = sample.σ[mod1(x+1,L), y, z]
        up     = sample.σ[x, mod1(y+1,L), z]
        front  = sample.σ[x, y, mod1(z+1,L)]
        SiSj  += sample.σ[x,y,z] * (right + up + front)
        Si    += sample.σ[x,y,z]
    end
    SiSj /= 3*L*L*L
    Si   /= L*L*L
    return (SiSj-Si*Si)^2
end

#
#
# Measuring utilities.
#
#

# +-+-+--+- → BitArray
line_to_spins(line) = [x=='+'?true:false for x in line] |> BitArray
# Bitarray to spin list
function bits_to_spins(b,L)
    # From bits to +1 and -1
    ints = [x ? Int8(+1) : Int8(-1) for x in b]
    # Reshape, return
    return reshape(ints,(L,L,L))
end

# Get fields from a file created with annealer.c
function get_data(file)
    data = readlines(file)

    # Get the global
    Lline = split(data[1])
    L = parse(Int64,Lline[2])

    Jup_line = split(data[9])
    Jup = line_to_spins(Jup_line[2])

    Jright_line = split(data[10])
    Jright = line_to_spins(Jright_line[2])

    Jfront_line = split(data[11])
    Jfront = line_to_spins(Jfront_line[2])

    # Iterate over the rest of the lines, build a list of spins, β and
    # mc step.
    mcs   = Int64[]
    βs    = Float64[]
    spins = BitArray[]

    for line in data[12:end]
        #if is a comment, do nothing.
        if line[1] != '#'
            # Read mc, β, spin list
            fields = split(line,['\t','\n'])
            push!(mcs   ,parse(Int64  ,fields[1]) )
            push!(βs    ,parse(Float64,fields[2]) )
            push!(spins ,line_to_spins(fields[3]) )
        end
    end

    # Return everything
    return Jup,Jright,Jfront,mcs,βs,spins
end

# Parse the fields of a file created with annealer.c in a DataFrame,
# measuring the observables.
# Observables is a string that can contain E,C,χ,M.
function measure(file;T0=10,observables="")
    Jup,Jright,Jfront,mcs,βs,spins = get_data(file)
    Jhash = hash([Jup Jright Jfront])
    N = length(spins[1])
    L = floor(Int64,cbrt(N))
    Nβ = 1 + (βs |> diff |> countnz) # Number of β changes + 1
    # Build a tentative net with the J's. Spins will be changed on the
    # fly.
    net = Net( bits_to_spins(Jup,L)
               ,bits_to_spins(Jright,L)
               ,bits_to_spins(Jfront,L)
               ,zeros(Int8,L,L,L)
               ,L    )
    # Fill this.
    df = DataFrame()
    for i in 1:length(spins)
        σ   = spins[i]
        # Fill the net spins
        S = bits_to_spins(σ,L)
        net.σ = S
        # Compute observables if asked for
        M       = NaN
        E       = NaN
        χsglass = NaN
        if 'M' ∈ observables
            M       = magnetization(σ)
        end
        if 'E' ∈ observables
            E       = energy(net)
        end
        if 'χ' ∈ observables || 'X' ∈ observables
            χsglass = χsg(net)
        end
        # For the correlation, we need to compare spins[i] and
        # spins[i+T0]. But i is the measurement number, and we need
        # the measurement number at that temperature, nothing more.
        # That's easy, because the number of measurements per β is
        # just length(spins)/Nβ. Use a modulo operator.
        C = NaN
        if 'C' ∈ observables
            Nmeas_per_β = length(spins)/Nβ
            if mod1(i,Nmeas_per_β) > Nmeas_per_β-T0
                C = NaN
            else
                # Correlation is just the scalar product. Because they are
                # bit arrays with elements ∈(0,1), they must translated to
                # spins in ∈(-1,1).
                thisσ = S
                nextσ = bits_to_spins(spins[i+T0], L)
                # Now, they are cubes of spins. We can't use dot product
                # on that kind of tensor; use a workaround:
                C = map(*,thisσ,nextσ) |> mean
            end
        end
        μdf = DataFrame( mc  = mcs[i]
                        ,β   = βs[i]
                        ,J   = Jhash
                        ,M   = M
                        ,E   = E
                        ,χsg = χsglass
                        ,C   = C
                        ,C2  = βs[i]*(1-C)  )
        rename!(μdf,:C2,Symbol("β(1-C)"))
        df = vcat(df,μdf)
    end
    return df
end

end#module

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
using Base.Iterators: partition
using Flux
using Flux.Optimise: update!
using Flux: logitbinarycrossentropy
using Images
using MLDatasets
using Statistics
using Parameters: @with_kw
using Printf
using Random

@with_kw struct HyperParams
    batch_size::Int = 128
    latent_dim::Int = 100
    epochs::Int = 20
    verbose_freq::Int = 1000
    output_x::Int = 6
    output_y::Int = 6
    lr_dscr::Float64 = 0.0002
    lr_gen::Float64 = 0.0002
end

function create_output_image(gen, fixed_noise, hparams)
    @eval Flux.istraining() = false
    fake_images = @. cpu(gen(fixed_noise))
    @eval Flux.istraining() = true
    image_array = permutedims(dropdims(reduce(vcat, reduce.(hcat, partition(fake_images, hparams.output_y))); dims=(3, 4)), (2, 1))
    image_array = @. Gray(image_array + 1f0) / 2f0
    return image_array
end

function Discriminator()
    return Chain(
            Conv((4, 4), 1 => 64; stride = 2, pad = 1),
            x->leakyrelu.(x, 0.2f0),
            Dropout(0.25),
            Conv((4, 4), 64 => 128; stride = 2, pad = 1),
            x->leakyrelu.(x, 0.2f0),
            Dropout(0.25),
            x->reshape(x, 7 * 7 * 128, :),
            Dense(7 * 7 * 128, 1))
end

function Generator()
    return Chain(
            Dense(hparams.latent_dim, 7 * 7 * 256),
            BatchNorm(7 * 7 * 256, relu),
            x->reshape(x, 7, 7, 256, :),
            ConvTranspose((5, 5), 256 => 128; stride = 1, pad = 2),
            BatchNorm(128, relu),
            ConvTranspose((4, 4), 128 => 64; stride = 2, pad = 1),
            BatchNorm(64, relu),
            ConvTranspose((4, 4), 64 => 1, tanh; stride = 2, pad = 1),
            )
end

# Loss functions
function discriminator_loss(real_output, fake_output)
    real_loss = mean(logitbinarycrossentropy.(real_output, 1f0))
    fake_loss = mean(logitbinarycrossentropy.(fake_output, 0f0))
    return real_loss + fake_loss
end

generator_loss(fake_output) = mean(logitbinarycrossentropy.(fake_output, 1f0))

function train_discriminator!(gen, dscr, x, opt_dscr, hparams)
    noise = randn!(similar(x, (hparams.latent_dim, hparams.batch_size)))
    fake_input = gen(noise)
    ps = Flux.params(dscr)
    # Taking gradient
    loss, back = Flux.pullback(ps) do
        discriminator_loss(dscr(x), dscr(fake_input))
    end
    grad = back(1f0)
    update!(opt_dscr, ps, grad)
    return loss
end

function train_generator!(gen, dscr, x, opt_gen, hparams)
    noise = randn!(similar(x, (hparams.latent_dim, hparams.batch_size)))
    ps = Flux.params(gen)
    # Taking gradient
    loss, back = Flux.pullback(ps) do
        generator_loss(dscr(gen(noise)))
    end
    grad = back(1f0)
    update!(opt_gen, ps, grad)
    return loss
end

function train(; kws...)
    # Model Parameters
    hparams = HyperParams(; kws...)

    # Load MNIST dataset
    images, _ = MLDatasets.MNIST.traindata(Float32)
    # Normalize to [-1, 1]
    image_tensor = reshape(@.(2f0 * images - 1f0), 28, 28, 1, :)
    # Partition into batches
    data = [image_tensor[:, :, :, r] |> gpu for r in partition(1:60000, hparams.batch_size)]

    fixed_noise = [randn(hparams.latent_dim, 1) |> gpu for _=1:hparams.output_x*hparams.output_y]

    # Discriminator
    dscr = Discriminator() |> gpu

    # Generator
    gen =  Generator() |> gpu

    # Optimizers
    opt_dscr = ADAM(hparams.lr_dscr)
    opt_gen = ADAM(hparams.lr_gen)

    # Training
    train_steps = 0
    for ep in 1:hparams.epochs
        @info "Epoch $ep"
        for x in data
            # Update discriminator and generator
            loss_dscr = train_discriminator!(gen, dscr, x, opt_dscr, hparams)
            loss_gen = train_generator!(gen, dscr, x, opt_gen, hparams)

            if train_steps % hparams.verbose_freq == 0
                @info("Train step $(train_steps), Discriminator loss = $(loss_dscr), Generator loss = $(loss_gen)")
                # Save generated fake image
                output_image = create_output_image(gen, fixed_noise, hparams)
                save(@sprintf("output/dcgan_steps_%06d.png", train_steps), output_image)
            end
            train_steps += 1
        end
    end

    output_image = create_output_image(gen, fixed_noise, hparams)
    save(@sprintf("output/dcgan_steps_%06d.png", train_steps), output_image)
end

cd(@__DIR__)
train()

max_stars_repo_name: thomascherickal/model-zoo

max_issues_repo_name: thomascherickal/model-zoo

max_forks_repo_name: thomascherickal/model-zoo


text:
"""
    A System is a set of variables and equations to be solved simultaneously.

    The parameterization is for the type used for Variable values.
    This controls whether or not automatic diffentiation is used.

    If T = Float64: Then Variable values are just POD.
                    In this case `computeQpJacobian` statements will be used to form the Jacobian

    If T = Dual{X, Float64}: Then Variable values are automatically differentiated types
                             In this case the Jacobians of Residuals will be automatically
                             computed using the partials of the Dual numbers
"""
type System{T}
    " The Mesh this System will use "
    mesh::Mesh

    " All of the Variables belonging to this System "
    variables::Array{Variable{T}}

    " All of the Kernels belonging to this System "
    kernels::Array{Kernel}

    " All of the BoundaryConditions belonging to this System "
    bcs::Array{BoundaryCondition}

    " The total number of degrees of freedom "
    n_dofs::Int64

    " Number of local degrees of freedom "
    local_n_dofs::Int64

    " The first DOF on this processor "
    first_local_dof::Int64

    " The last DOF on this processor "
    last_local_dof::Int64

    " A list of DoFs that need to be ghosted to this processor.  These are DoFs which are off-processor, but we need their values for Assembly. "
    ghosted_dofs::Set{Int64}

    " Whether or not initialize!() has been called for this System "
    initialized::Bool

    " The Quadrature Rule "
    q_rule::QuadratureRule

    " The Function Space "
    func_space::FunctionSpace

    " The Current FEValues "
    fe_values::FECellValues

    " Number of local dofs connected to each local dof.  Must call `generateSparsity!()` to fill "
    local_dof_sparsity::Array{Int32}

    " Number of non-local dofs connected to each local dof. Must call `generateSparsity!()` to fill "
    off_proc_dof_sparsity::Array{Int32}

    " A Mesh must be provided "
    System(mesh::Mesh) = (x = new(mesh,
                                  Array{Variable}(0),
                                  Array{Kernel}(0),
                                  Array{BoundaryCondition}(0),
                                  0,
                                  0,
                                  0,
                                  0,
                                  Set{Int64}(),
                                  false,
                                  QuadratureRule{2,RefCube}(:legendre, 2),
                                  Lagrange{2, RefCube, 1}()) ;
                          x.fe_values = FECellValues(x.q_rule, x.func_space);
                          return x)
end

" Add a Variable named name to the System "
function addVariable!{T}(sys::System{T}, name::String)
    n_vars = length(sys.variables)

    var = Variable{T}(n_vars+1, name)

    push!(sys.variables, var)

    return var
end

" Add a Kernel to the System "
function addKernel!(sys::System, kernel::Kernel)
    push!(sys.kernels, kernel)
end

" Add a BC to the System "
function addBC!(sys::System, bc::BoundaryCondition)
    push!(sys.bcs, bc)
end

" Distribute the DoFs to the DofObjects "
function distributeDofs(sys::System)
    startLog(main_perf_log, "distributeDofs()")

    # These will be "node major"
    n_vars = length(sys.variables)

    # Find the number of dofs on each processor
    dofs_per_proc = Array{Int64}(MPI.Comm_size(MPI.COMM_WORLD))
    fill!(dofs_per_proc, 0)

    # Here's the way this works... what we want is the dof IDs for each processor to be contiguous
    # In addition, we want the dof IDs to be increasing with processor_id
    # To do that, we're going to find the total number of DoFs on each processor and communicate
    # that into a vector on every processor.  Next, we'll compute the cumsum of that vector.
    # That cumsum will give us the starting DoF ID for each processor (if we add 1 to it)
    # with one small change: the first processor needs to start at 1...
    # and the final entry in the cumsum is useless
    # So we'll do a little vector manipulation to fix that up.
    #
    # At that point we can loop over the elements and assign dof IDs for ALL dof objects
    # (Yes, even for the ones not owned by this processor)
    # We'll keep track of the current DoF ID for each processor as we go.
    # Every processor will be doing the same thing simultaneously so all processors will
    # end up with exactly the same information.
    #
    # You could attempt to parallelize this... but the gains are minimal unless you want
    # to go to a true distributed mesh.  The reason why is that you would need to add
    # a communication step... but you would still have to assign them to every node on every
    # processor anyway.  At least by doing it this way there is no communication.

    # For now we're just going to say that the number of dofs is the number of
    # local nodes * number of variables.  This would be much more complicated to do if we
    # had element degrees of freedom or higher order shape functions
    sys.local_n_dofs = n_vars * length(sys.mesh.local_nodes)

    dofs_per_proc = MPI.Allgather([sys.local_n_dofs], 1, MPI.COMM_WORLD)

    proc_current_dof = cumsum(dofs_per_proc)

    proc_current_dof .+= 1

    # The first processor is going to start at 1
    unshift!(proc_current_dof, 1)

    # Remove the last entry
    pop!(proc_current_dof)

    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    # Set our local information
    sys.first_local_dof = proc_current_dof[proc_id + 1] # Don't forget that MPI is 0-based!
    sys.last_local_dof = sys.first_local_dof + sys.local_n_dofs - 1

    # Go over the elements and set the DoFs for the nodes (if they haven't already been set)
    # We go over elements for a couple of reasons... but the main one for now is to match
    # the way libMesh does it to make it easy to compare Jacobian/Residual matrices!
    for elem in sys.mesh.elements
        for node in elem.nodes
            node_proc_id = node.processor_id
            if length(node.dofs) == 0 # Only set DoFs if they haven't been set before
                current_dof = proc_current_dof[node.processor_id + 1] # Don't forget about 0-based indexing!

                node.dofs = [x for x in current_dof:((n_vars+current_dof)-1)]

                proc_current_dof[node.processor_id + 1] += n_vars
            end
        end
    end

    # Save off the total number of DoFs distributed
    sys.n_dofs = sum(dofs_per_proc)

    stopLog(main_perf_log, "distributeDofs()")
end

"""
    Find all of the off-processor DoFs that we will need for Assembly
"""
function findGhostedDofs!(sys::System)
    # Basic plan:
    # 1. Loop over local elements
    # 2. Interrogate nodes connected to those elems...
    # 3. If they are off-processor then ghost their dofs

    # Note!  This doesn't create any sort of "halo"...
    # it just barely gets the info we need for Lagrange FE Assembly
    # A lot more could be done here

    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    for elem in sys.mesh.local_elements
        for node in elem.nodes
            if node.processor_id != proc_id
                dofs = node.dofs
                for dof in node.dofs
                    push!(sys.ghosted_dofs, dof)
                end
            end
        end
    end
end

"""
    PRIVATE: Helper function
"""
function _addDofs!(sys::System, original_node::Node, other_node::Node, local_dofs::Array, off_proc_dofs::Array)
    for dof in original_node.dofs
        local_dof = (dof - sys.first_local_dof) + 1

        for other_dof in other_node.dofs
            if sys.first_local_dof <= other_dof && other_dof <= sys.last_local_dof
                local_dofs[local_dof] += 1
            else
                off_proc_dofs[local_dof] += 1
            end
        end
    end
end

"""
    Generate Sparsity information
"""
function generateSparsity!(sys::System)
    startLog(main_perf_log, "generateSparsity!()")

    # Basic plan:
    # 1. Loop over local nodes
    # 2. Add up local dofs on each node as local nonzeros
    # 3. Use node_to_elem_map to get all elems connected to the node
    # 4. Go over nodes connected to those elements and add up the number of local and off-processor dofs

    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    local_dofs = Array{Int32}(sys.local_n_dofs)
    fill!(local_dofs, 0)
    off_proc_dofs = Array{Int32}(sys.local_n_dofs)
    fill!(off_proc_dofs, 0)

    for node in sys.mesh.local_nodes
        # Couple all of the dofs on this node to eachother
        _addDofs!(sys, node, node, local_dofs, off_proc_dofs)

        connected_elems = sys.mesh.node_to_elem_map[node.id]

        visited_nodes = Set{Int64}()

        push!(visited_nodes, node.id)

        for elem in connected_elems
            for other_node in elem.nodes
                if !(other_node.id in visited_nodes)
                    _addDofs!(sys, node, other_node, local_dofs, off_proc_dofs)
                    push!(visited_nodes, other_node.id)
                end
            end
        end
    end

    sys.local_dof_sparsity = local_dofs
    sys.off_proc_dof_sparsity = off_proc_dofs

    stopLog(main_perf_log, "generateSparsity!()")
end

"""
    Should be called after all Variables have been added to the System and before
    a `Solver` is created.

    The main purpose of this function is to distribute the DoFs across the mesh.
"""
function initialize!(sys::System)
    startLog(main_perf_log, "initialize!(System)")

    distributeDofs(sys)

    # Set a flag that this System is initialized
    sys.initialized = true

    stopLog(main_perf_log, "initialize!(System)")
end

"""
    Grab all of the DoF indices for the current element for a given variable
"""
function connectedDofIndices(elem::Element, var::Variable)
    [ node.dofs[var.id] for node in elem.nodes ]
end

"""
    Create the vector of dof_values from the dof_indices and solution vector

    Specialization for Float64
"""
function dofValues!(dof_values::Array, sys::System{Float64}, solution::AbstractArray, var::Variable, dof_indices::Array)
    resize!(dof_values, length(dof_indices))
    dof_values[:] = solution[dof_indices]
end

"""
    Create the vector of dof_values from the dof_indices and solution vector

    Specialization for Dual
"""
function dofValues!{N,T}(dof_values::Array, sys::System{Dual{N,T}}, solution::AbstractArray, var::Variable, dof_indices::Array)
    resize!(dof_values, length(dof_indices))

    values = solution[dof_indices]

    n_vars = length(sys.variables)
    n_dofs_per_var = length(dof_indices)
    total_n_dofs = n_vars * n_dofs_per_var

    @assert total_n_dofs <= N

    partial_offset = ((var.id - 1) * n_dofs_per_var)

    for i in 1:n_dofs_per_var
        dof_values[i] = dualVariable(values[i], partial_offset + i, N)
    end
end

"""
    Reinitialize all of the data and objects in the system for the current Element
"""
function reinit!{T}(sys::System{T}, elem::Element, solution::AbstractArray)
    # Grab all of the coordinattes for the current element
    coords = [node.coords for node in elem.nodes]

    # Reinitialize the shape functions
    JuAFEM.reinit!(sys.fe_values, coords)

    # Reinitialize the variable values
    for var in sys.variables
        # Need to grab the dof_indices for this variable on this element
        dof_indices = connectedDofIndices(elem, var)

        # Now pull out those pieces of the solution vector
        dof_values = Array{T}(0)
        dofValues!(dof_values, sys, solution, var, dof_indices)

        # Recompute the Variable values
        reinit!(var, sys.fe_values, dof_indices, dof_values)
    end
end

" Helper function for getting a Float64 dof_value at a Node "
function dofValue(solution::AbstractArray, var::Variable{Float64}, dof_index::Int64)
    return solution[dof_index]
end

" Helper function for getting a Dual dof_value at a Node "
function dofValue{N,T}(solution::AbstractArray, var::Variable{Dual{N,T}}, dof_index::Int64)
    return dualVariable(solution[dof_index], var.id, N)
end

"""
    Reinitialize all of the data and objects in the system for the current Node
"""
function reinit!{T}(sys::System{T}, node::Node, solution::AbstractArray)
    # Reinitialize the variable values
    for var in sys.variables
        dof_index = node.dofs[var.id]
        dof_value = dofValue(solution, var, dof_index)

        # Recompute the Variable values
        reinit!(var, dof_index, dof_value)
    end
end


text:
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __julia-version__ = 1.7.2
# __author__        = "Max Kannenberg"
# __copyright__     = "2020-2022"
# __license__       = "ISC"

## create a moving section and its containing characteristic sections with secured braking, accelerating and cruising behavior
function determineCharacteristics(path::Path, train::Train, settings::Settings)
    movingSection = createMovingSection(path, train.v_limit, train.length)
    movingSection = secureBrakingBehavior!(movingSection, train.a_braking)
    movingSection = secureAcceleratingBehavior!(movingSection, settings, train)
    #movingSection = secureCruisingBehavior!(movingSection, settings, train)

    return movingSection
end #function determineCharacteristics

## define the intersection velocities between the characterisitc sections to secure braking behavior
function secureBrakingBehavior!(movingSection::Dict, a_braking::Real)
    # this function limits the entry and exit velocity of the characteristic sections to secure that the train stops at the moving sections end
        CSs = movingSection[:characteristicSections]

        csId = length(CSs)
        followingCSv_entry = 0.0     # the exit velocity of the last characteristic section is 0.0 m/s
        while csId >= 1
            CS = CSs[csId]

            CS[:v_exit] = min(CS[:v_limit], followingCSv_entry)

            v_entryMax = calcBrakingStartVelocity(CS[:v_exit], a_braking, CS[:length])

            CS[:v_entry] = min(CS[:v_limit], v_entryMax)
            CS[:v_peak] = CS[:v_entry]

            # reset the characteristic section (CS), delete behavior sections (BS) that were used during the preperation for setting v_entry, v_peak and v_exit
            CS[:behaviorSections] = Dict()
            CS[:E] = 0.0
            CS[:t] = 0.0

            followingCSv_entry = CS[:v_entry]
            csId = csId - 1
        end #while
    return movingSection
end #function secureBrakingBehavior!

## define the intersection velocities between the characterisitc sections to secure accelerating behavior
function secureAcceleratingBehavior!(movingSection::Dict, settings::Settings, train::Train)
    # this function limits the entry and exit velocity of the characteristic sections in case that the train accelerates in every section and cruises aterwards
    CSs = movingSection[:characteristicSections]

    CSs[1][:v_entry] = 0.0     # the entry velocity of the first characteristic section is 0.0 m/s
    startingPoint = createDataPoint()
    startingPoint[:i] = 1

    previousCSv_exit = CSs[1][:v_entry]
    for CS in CSs
        CS[:v_entry] = min(CS[:v_entry], previousCSv_exit)
        startingPoint[:s] = CS[:s_entry]
        startingPoint[:v] = CS[:v_entry]
        calculateForces!(startingPoint, CSs, CS[:id], "accelerating", train, settings.massModel) # traction effort and resisting forces (in N)
        acceleratingCourse::Vector{Dict} = [startingPoint]    # List of data points

        if CS[:v_entry] < CS[:v_peak]
            # conditions for entering the accelerating phase
            stateFlags = Dict(:endOfCSReached => false,
                              :brakingStartReached => false,
                              :tractionDeficit => false,
                              :resistingForceNegative => false,
                              :previousSpeedLimitReached => false,
                              :speedLimitReached => false,
                              :error => false,
                              :usedForDefiningCharacteristics => true)      # because usedForDefiningCharacteristics == true the braking distance will be ignored during securing the accelerating phase
            v_peak = CS[:v_entry]
            (CS, acceleratingCourse, stateFlags) = addBreakFreeSection!(CS, acceleratingCourse, stateFlags, settings, train, CSs)
            while !stateFlags[:speedLimitReached] && !stateFlags[:endOfCSReached]
                if !stateFlags[:tractionDeficit]
                    if !stateFlags[:previousSpeedLimitReached]
                    (CS, acceleratingCourse, stateFlags) = addAcceleratingSection!(CS, acceleratingCourse, stateFlags, settings, train, CSs)        # this function changes the acceleratingCourse

                    elseif stateFlags[:previousSpeedLimitReached]
                        (CS, acceleratingCourse, stateFlags) = addClearingSection!(CS, acceleratingCourse, stateFlags, settings, train, CSs)        # this function is needed in case the train is not allowed to accelerate because of a previous speed limit
                    end
                else
                    if settings.massModel == :mass_point || acceleratingCourse[end][:s] > CS[:s_entry] + train.length
                        break
                    else
                        (CS, acceleratingCourse, stateFlags) = addDiminishingSection!(CS, acceleratingCourse, stateFlags, settings, train, CSs)        # this function is needed in case the resisitng forces are higher than the maximum possible tractive effort
                    end
                end
                v_peak = max(v_peak, acceleratingCourse[end][:v])
            end

#            CS[:v_peak] = max(CS[:v_entry], acceleratingCourse[end][:v])
            CS[:v_peak] = v_peak
            CS[:v_exit] = min(CS[:v_exit], CS[:v_peak], acceleratingCourse[end][:v])
        else #CS[:v_entry] == CS[:v_peak]
            # v_exit stays the same
        end #if

        previousCSv_exit = CS[:v_exit]

        # reset the characteristic section (CS), delete behavior sections (BS) that were used during the preperation for setting v_entry, v_peak and v_exit
        CS[:behaviorSections] = Dict()
        CS[:E] = 0.0
        CS[:t] = 0.0
    end #for

    return movingSection
end #function secureAcceleratingBehavior!

#=
## define the intersection velocities between the characterisitc sections to secure cruising behavior
function secureCruisingBehavior!(movingSection::Dict, settings::Settings, train::Train)
    # limit the exit velocity of the characteristic sections in case that the train cruises in every section at v_peak
    CSs = movingSection[:characteristicSections]

    startingPoint = createDataPoint()
    startingPoint[:i] = 1

    previousCSv_exit = CSs[1][:v_entry]

    for CS in CSs
        # conditions for entering the cruising phase
        stateFlags = Dict(:endOfCSReached => false,
                          :brakingStartReached => false,
                          :tractionDeficit => false,
                          :resistingForceNegative => false,
                          :previousSpeedLimitReached => false,
                          :speedLimitReached => false,
                          :error => false,
                          :usedForDefiningCharacteristics => true)

        CS[:v_entry] = min(CS[:v_entry], previousCSv_exit)

        startingPoint[:s] = CS[:s_entry]
        startingPoint[:v] = CS[:v_peak]
        cruisingCourse::Vector{Dict} = [startingPoint]    # List of data points

        while !stateFlags[:endOfCSReached] #&& s_cruising > 0.0
            if !stateFlags[:tractionDeficit]
                s_cruising = CS[:s_exit] - cruisingCourse[end][:s]
                if !stateFlags[:resistingForceNegative]# cruisingCourse[end][:F_R] >= 0
                    (CS, cruisingCourse, stateFlags) = addCruisingSection!(CS, cruisingCourse, stateFlags, s_cruising, settings, train, CSs, "cruising")        # this function changes the cruisingCourse
                else
                    (CS, cruisingCourse, stateFlags) = addCruisingSection!(CS, cruisingCourse, stateFlags, s_cruising, settings, train, CSs, "downhillBraking")
                end
            else
                if settings.massModel == :mass_point || cruisingCourse[end][:s] > CS[:s_entry] + train.length
                    break
                else
                    (CS, cruisingCourse, stateFlags) = addDiminishingSection!(CS, cruisingCourse, stateFlags, settings, train, CSs)        # this function is needed in case the resisitng forces are higher than the maximum possible tractive effort
                end
            end
        end

        CS[:v_exit] = min(CS[:v_exit], cruisingCourse[end][:v])

        previousCSv_exit = CS[:v_exit]

        # reset the characteristic section (CS), delete behavior sections (BS) that were used during the preperation for setting v_entry, v_peak and v_exit
        CS[:behaviorSections] = Dict()
        CS[:E] = 0.0
        CS[:t] = 0.0
    end #for

    return movingSection
end #function secureCruisingBehavior!
=#


text:
@everywhere using JuMP
@everywhere using Ipopt
@everywhere using MathOptInterface
@everywhere using Optim
@everywhere using LinearAlgebra
@everywhere using DelimitedFiles
@everywhere using Distances

@everywhere function evalSingleSample(exposures, allSignatures, genome)

    maxMutations = sum(genome);

    numberOfSignatures = count(!iszero, exposures);

    x0 = maxMutations * rand(numberOfSignatures);
    x0 = x0 / sum(x0) * maxMutations;
    A = ones(1, numberOfSignatures);
    max_channels = size(genome, 1);

    # select only the signatures with exposure
    subSignatures = allSignatures[:,findall(!iszero, exposures)];

    # JuMP model instance
    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0, "sb" => "yes"));

    # set variables with bounds
    @variable(model, 0 <= x[i=1:numberOfSignatures] <= maxMutations, start=x0[i]);

    # set objective
    @NLobjective(model, Min,  sum( (sum(subSignatures[i,j]*x[j] for j=1:numberOfSignatures) - genome[i])^2 for i=1:max_channels));

    # set additional constraints
    @constraint(model, sum(A[i]*x[i] for i=1:numberOfSignatures) == maxMutations);

    # solve it
    JuMP.optimize!(model);
    x = JuMP.value.(x);

    x = round.(x);
    maximum_x_index = argmax(x);

    if sum(x) != maxMutations
         x[maximum_x_index] = x[maximum_x_index] + maxMutations - sum(x);
    end

    exposuresSample = zeros(size(allSignatures, 2));
    exposuresSample[findall(!iszero, exposures)] = x ;
    recon = allSignatures*exposuresSample;
    accr = 1-cosine_dist(recon, genome);
    error_fro = norm(recon-genome);
    norm_genome_fro = norm(genome);
    frob_rel_div = error_fro/norm_genome_fro;
    norm_one_dif = norm(recon - genome, 1);

    return accr, frob_rel_div, norm_one_dif, exposuresSample
end

max_stars_repo_name: ypradat/sigprofilerjulia

max_issues_repo_name: ypradat/sigprofilerjulia

max_forks_repo_name: ypradat/sigprofilerjulia