using JuMP, Gurobi, LinearAlgebra

# Problem data
v = [9, 3, 2, 6, 10, 7, 6, 6, 4, 5, 4, 1, 1, 
    6, 3, 2, 3, 2, 10, 7, 9, 8, 3, 10, 8, 5, 
    3, 8, 3, 6, 5, 7, 8, 6, 9, 7, 5, 5, 1, 5, 
    9, 5, 4, 5, 5, 3, 4, 8, 8, 10]
w = [4, 10, 9, 8, 8, 4, 7, 3, 1, 10, 5, 4, 2, 
    3, 9, 9, 9, 5, 8, 8, 4, 2, 6, 10, 5, 7, 7, 
    8, 4, 7, 7, 8, 4, 4, 10, 3, 9, 2, 9, 10, 3, 
    7, 3, 5, 5, 7, 10, 10, 5, 8]
C = 50
n = length(v)

knapsack_model = Model(Gurobi.Optimizer)

# I am turing these off just for illustration, 
# otherwise Gurobi ends up solving the problem during the presolve phase 😲

set_optimizer_attribute(knapsack_model, "Presolve", 0)

set_optimizer_attribute(knapsack_model, "Heuristics", 0)

@variable(knapsack_model, x[1:n] >= 0, Bin)

@objective(knapsack_model, Max, v' * x)

@constraint(knapsack_model, w' * x <= C)

optimize!(knapsack_model)

println(abs.(value.(x)))

objective_value(knapsack_model)

small_knapsack_model = Model(Gurobi.Optimizer)
@variable(small_knapsack_model , x, Bin)
@variable(small_knapsack_model , y, Bin)
@constraint(small_knapsack_model , x+2*y<=1.5)
@objective(small_knapsack_model , Max, x+y)
optimize!(small_knapsack_model )
println("X is: ", value(x), " Y is: ", value(y))

# Load the data
include("img//mip_data.jl");

println(f)

using JuMP, Gurobi

sfMipModel =  Model(Gurobi.Optimizer)

set_optimizer_attribute(sfMipModel, "MIPFocus", 3)
# If you are more interested in good quality feasible solutions, you can select MIPFocus=1.
# If you believe the solver is having no trouble finding the optimal solution, and wish to focus more
# attention on proving optimality, select MIPFocus=2.
# If the best objective bound is moving very slowly (or not at all), you may want to try MIPFocus=3 to focus on the bound.

 @variable(sfMipModel, x[1:n] >=0)

@variable(sfMipModel, y[1:p] >= 0, Bin)

@objective(sfMipModel, Min, sum(c[i] * x[i] for i in 1:n)+sum(d[i]*y[i] for i in 1:p))

for i in 1:m
    @constraint(sfMipModel, sum(A[i,j]*x[j] for j in 1:n)+ sum(B[i,j]*y[j] for j in 1:p) == f[i])
end

optimize!(sfMipModel)

x_sol = value.(x)    

y_sol = value.(y)

# data
n = 20
p = 5
m = 5
A = [-0.25149258982896966 0.525287826731042 -0.6656342865169662 -1.2960398772338568 0.9101252450751471 -0.07298248406603706 -1.9148704340623688 -0.5073995272298952 -1.6558511941286744 2.1006659624182586 -0.502359026079171 -0.1708406495941207 1.4161030569819415 0.33735905844690145 -0.10748103771152817 -1.4416340194991104 -0.9578960880328228 -0.1738984790022571 -1.0517643292393826 -0.8271987403702172; -0.8518984622275042 -1.374417162491773 0.35920232369365845 0.8437691247979793 -0.6281577605465152 0.5921264379488904 -0.12756078132890084 -0.44447204507515287 0.7293549470299907 0.715267254601397 -0.9155882952842251 0.35378971887340743 -0.8043683142766567 0.5324223554306367 0.2088297149344946 -0.5594219995662968 0.5632674521223371 -0.2302672262989066 0.7676987968759004 1.0910836780489503; 1.473815753336459 -0.6525879104965642 0.8519411939583414 -0.05756789248416348 -0.7125581880426727 0.4736057460760701 -2.231011632648033 0.5409490591899032 0.1583459153748212 -1.1509507228336586 1.3994115761547934 0.14111034486559706 -0.15123003344570538 -0.547672888665602 -1.9555427017941476 -1.8533179907378834 0.13798098584535193 0.2337797562353011 -0.8404137158789987 0.060446173222487264; 0.6396537082966063 0.856847789735652 1.4805326346158039 0.012324563610194187 1.3977381658753554 1.3342059688961259 -0.4784014905808825 -0.14890693050929457 0.8155899237370289 -0.5697405099766946 1.3446291535605297 0.6840341944716634 1.2484251582871324 1.2891201777341281 -0.2791170028219365 0.43803145770984425 0.4677466100676112 1.3639772554672074 0.8359137382225728 0.8723882764319004; 0.40517145458826054 1.169967487661965 1.7941045258417871 1.035188318563306 1.118495812473751 -0.04027314052794962 0.5438620238342292 -0.8181792084120169 0.08209406151847103 0.6469386439823567 -0.9467882212470308 -0.2866858366079981 0.8836746919594078 -2.0164799655225094 -0.42541330025958685 -0.09145359012880407 0.8920741578367347 -0.5093214870138416 -0.23118339134915017 -1.262095251923206]
D = [-2.0831200526934572 1.5615062392183003 0.18155475431473142 0.7599664330346061 0.2668125999785801 1.2573917809432806 0.04093911087719557 -0.21735884773163816 1.6340683454668137 1.2804063399285137 0.6443447646682994 -0.8939206481655969 -0.22507735354054656 -0.3094500241639488 0.2954219652735186 -0.3880055556581179 -0.9578091110001216 -0.922838957226308 -1.3571981131817064 -0.05064205091038371; 0.6324558903791979 -0.6652227624026216 -1.7276813422697317 -0.689479323056671 -0.7816950255371095 -1.802372714658206 0.14164169531740856 0.46527545010854515 1.1036334528775311 1.305688755430216 0.3797915504580763 1.7602878854059907 0.8155022785101097 -0.6315798865662005 -0.30581893855662595 -0.33305541977514697 -0.8854550960099433 0.03621076189783246 -0.2788693026530316 1.6348184993843053; -2.4242492532538096 0.5770923004878816 0.4961795092913091 -0.04867957395223246 0.6673192113922586 -0.8855274251005092 0.4568856956106312 -0.6097467636667095 -0.3639368866683512 0.5625062854293676 -0.5017933052998168 0.3836052143475206 1.8936021566371188 -0.209665904941543 -1.3278566925429307 0.6736435837748775 -0.5017805315844522 0.6452850279162458 0.7881379044649098 0.45931063639065656; -0.5372061672955294 -0.7867658024064041 0.09657107558669714 -1.0587831441218194 -1.2444111679260363 -0.6885359430850223 -0.4325022650433149 0.7912025990048928 0.1141092892049688 0.918388681602774 0.41101551577580714 1.3817408040581836 -0.030696867637192395 0.7523799268172908 -0.3083505028779359 0.27074896029609724 -0.5303233440211375 -1.1025104074431078 -2.3234729621969272 0.705781546034252; 0.6313935384296541 0.699436081582457 -0.5461167811714385 -0.5108885871403681 -1.59867913729006 -0.7350106637365103 -1.1033454910986003 0.898936477156721 -0.07365263955313238 -0.48330570412807555 1.0944353481674236 0.19009504425867302 -0.5106409146420512 0.647369019147565 -0.40169454919669595 -0.1575677053433486 1.4977052939251785 -0.39775746600800416 1.821753152147404 0.10823903833322045]
c = [-5.724274247670441; 0.872180569098502; -13.289441619771337; -1.284190729811427; -10.090248808669005; 19.226305603166242; -14.323991615118594; -5.8893709748795615; -9.801477473401599; -2.533420139004128; -4.017671140375999; -12.25747413009963; -6.5234635509559205; -18.204822439044932; 10.759517257360242; -2.678837137470244; 2.6260424150728934; -9.929728378562064; -1.08448828496412; -14.883711771646961]
b = [-7.3620594939851545; 1.2864238748369008; -5.191670678971517; 3.4603990262005064; -0.9930993374575348]
h = [-5.110534843409266; -0.12459704540647681; -2.8197063217437384; -1.9945660953801996; 3.5137561922861824];

# a proper code is as follows:

lazy_model = Model(Gurobi.Optimizer)

@variable(lazy_model, x[1:n], Bin)

for i in 1:m
    @constraint(lazy_model, sum(A[i,j]*x[j] for j in 1:n) <= b[i])
end

@objective(lazy_model, Min, sum(c[i] * x[i] for i in 1:n))


# lazy constraint is added through the following callback function
function my_callback_function(cb_data)
    # check the status of the solution at the branch-and-bound tree node 
    status = callback_node_status(cb_data, lazy_model)
    if status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
        # if we enter this block then it means that
        # `callback_value(cb_data, x)` is not integer (to some tolerance).
        # Our lazy constraint generator requires an
        # integer-feasible primal solution, so we can add a `return` here.
        return
    elseif status == MOI.CALLBACK_NODE_STATUS_INTEGER
        # `callback_value(cb_data, x)` is integer (to some tolerance).
        println("[😁] Current node has an integer solution")
    else
        @assert status == MOI.CALLBACK_NODE_STATUS_UNKNOWN
        # `callback_value(cb_data, x)` might be integer feasible or infeasible,
        # in such a case it is best if we call our lazy constraint
        # @assert ensures that if the node status is none of the three conditions above, 
        # then it will throw an error
    end
    x_val = callback_value.(cb_data, x) # load the value of the interim solution at the node 
    for i in 1:p
        if sum(D[i,j]*x_val[j] for j in 1:n) > h[i] + 1e-6 
            println("[😻] Submitting $(i)-th constraint d[i]'x <= h[i] as a lazy constraint")
            # 💀 note that we using x_val which is a known vector, not x which is the variable
            # also, added a small number 1e-6 to ensure numberical stability
            # now add the lazy constraint
            con = @build_constraint(sum(D[i,j]*x[j] for j in 1:n) <= h[i])
            MOI.submit(lazy_model, MOI.LazyConstraint(cb_data), con) #💀 need to submit the constraint to our model
            break
        end
    end
end

MOI.set(lazy_model, MOI.LazyConstraintCallback(), my_callback_function) 
# 💀 need to invoke this command so that JuMP knows that "my_callback_function" is used to generate lazy constraints

optimize!(lazy_model)

x_sol = value.(x)

println(x_sol[15:20])

sum(x_sol)

D*x_sol <= h

lazy_obj_val = objective_value(lazy_model)

# data
n = 20
m = 5
A = [-0.25149258982896966 0.525287826731042 -0.6656342865169662 -1.2960398772338568 0.9101252450751471 -0.07298248406603706 -1.9148704340623688 -0.5073995272298952 -1.6558511941286744 2.1006659624182586 -0.502359026079171 -0.1708406495941207 1.4161030569819415 0.33735905844690145 -0.10748103771152817 -1.4416340194991104 -0.9578960880328228 -0.1738984790022571 -1.0517643292393826 -0.8271987403702172; -0.8518984622275042 -1.374417162491773 0.35920232369365845 0.8437691247979793 -0.6281577605465152 0.5921264379488904 -0.12756078132890084 -0.44447204507515287 0.7293549470299907 0.715267254601397 -0.9155882952842251 0.35378971887340743 -0.8043683142766567 0.5324223554306367 0.2088297149344946 -0.5594219995662968 0.5632674521223371 -0.2302672262989066 0.7676987968759004 1.0910836780489503; 1.473815753336459 -0.6525879104965642 0.8519411939583414 -0.05756789248416348 -0.7125581880426727 0.4736057460760701 -2.231011632648033 0.5409490591899032 0.1583459153748212 -1.1509507228336586 1.3994115761547934 0.14111034486559706 -0.15123003344570538 -0.547672888665602 -1.9555427017941476 -1.8533179907378834 0.13798098584535193 0.2337797562353011 -0.8404137158789987 0.060446173222487264; 0.6396537082966063 0.856847789735652 1.4805326346158039 0.012324563610194187 1.3977381658753554 1.3342059688961259 -0.4784014905808825 -0.14890693050929457 0.8155899237370289 -0.5697405099766946 1.3446291535605297 0.6840341944716634 1.2484251582871324 1.2891201777341281 -0.2791170028219365 0.43803145770984425 0.4677466100676112 1.3639772554672074 0.8359137382225728 0.8723882764319004; 0.40517145458826054 1.169967487661965 1.7941045258417871 1.035188318563306 1.118495812473751 -0.04027314052794962 0.5438620238342292 -0.8181792084120169 0.08209406151847103 0.6469386439823567 -0.9467882212470308 -0.2866858366079981 0.8836746919594078 -2.0164799655225094 -0.42541330025958685 -0.09145359012880407 0.8920741578367347 -0.5093214870138416 -0.23118339134915017 -1.262095251923206]
c = [-5.724274247670441; 0.872180569098502; -13.289441619771337; -1.284190729811427; -10.090248808669005; 19.226305603166242; -14.323991615118594; -5.8893709748795615; -9.801477473401599; -2.533420139004128; -4.017671140375999; -12.25747413009963; -6.5234635509559205; -18.204822439044932; 10.759517257360242; -2.678837137470244; 2.6260424150728934; -9.929728378562064; -1.08448828496412; -14.883711771646961]
b = [-7.3620594939851545; 1.2864238748369008; -5.191670678971517; 3.4603990262005064; -0.9930993374575348]


# a proper code is as follows:

hrstc_model = Model(Gurobi.Optimizer)

set_optimizer_attribute(hrstc_model, "Presolve", 0)

set_optimizer_attribute(hrstc_model, "Heuristics", 0)

@variable(hrstc_model, x[1:n], Bin)

for i in 1:m
    @constraint(hrstc_model, sum(A[i,j]*x[j] for j in 1:n) <= b[i])
end

@objective(hrstc_model, Min, sum(c[i] * x[i] for i in 1:n))


# lazy constraint is added through the following callback function
function my_heuristic_callback_function(cb_data)
    # load the current values of the variables at a node of the branch-and-bound tree
    x_val = callback_value.(cb_data, x)
    if sum(x_val) <= 10 # this is the condition for applying the heuristic
        println("[🐼 ] Heuristic condition satisfied")
        println("[🐯 ] Status of the callback node is: $(callback_node_status(cb_data,hrstc_model))")
        # show the status of the callback node
        println("[🐰 ]The values of the variables at the callback node $(x_val)")
        # Load the heuristic solution values in a vertical vector pointwise
        heuristic_values = vcat(1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0)
        # Submit the heuristic solution for potentially improving the current solution
        status = MOI.submit(
             hrstc_model, MOI.HeuristicSolution(cb_data), x, heuristic_values
             )
      println("[🙀 ] Status of the submitted heuristic solution is: ", status) # The status shows if the submitted heuristic solution is accepted or not
  end
end

# IMPORTANT: 💀  This enables the heuristic
MOI.set(hrstc_model, MOI.HeuristicCallback(), my_heuristic_callback_function)

optimize!(hrstc_model)

x_hrstc_opt = value.(x)

println(x_hrstc_opt)

# a proper code is as follows:

x_warm_start = [1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0];

ws_model = Model(Gurobi.Optimizer)

set_optimizer_attribute(ws_model, "Presolve", 0)

set_optimizer_attribute(ws_model, "Heuristics", 0)

@variable(ws_model, x2[1:n], Bin)

for i in 1:m
    @constraint(ws_model, sum(A[i,j]*x2[j] for j in 1:n) <= b[i])
end

# warm-starting here:

set_start_value.(x2, x_warm_start)


@objective(ws_model, Min, sum(c[i] * x2[i] for i in 1:n))

optimize!(ws_model)

x_ws_opt = value.(x)

# Data, change it accordingly
# ---------------------------

m = 5
n = 10
A = randn(m,n)
b = randn(m)
M = 1
k = convert(Int64, round(m/3))

# Renaming a bunch of variables
S = A'*A
c = -2*A'*b
d = norm(b)^2

# Define the model
# ----------------

model = Model(with_optimizer(Gurobi.Optimizer)) # define name of the model, it could be anything, not necessarily "model"

# Variables
# ---------

@variable(model, x[1:n]) # define variable x 

@variable(model, y[1:n], Bin) # define the binary variable y

# Objective
# ---------

sense = MOI.MIN_SENSE # by this command, we are programatically defining a quadratic objective to be minimized 

@objective(model, sense, sum(S[i,j]*x[i]*x[j] for i in 1:n, j in 1:n)+ sum(c[i]*x[i] for i in 1:n) + d) # define the objective

# Constraints
# ------------

@constraint(model, con_lb[i=1:n], -M*y[i] <= x[i]) # lower bound constraint

@constraint(model, con_ub[i=1:n], x[i] <= M*y[i]) # upper bound constraint

@constraint(model, con_bd_sum, sum(y[i] for i in 1:n) <= k) # cardinality constraint in terms of y

# Run the optimizer
# -----------------

status=optimize!(model) # time to optimize!