# Data for the problem
#---------------------
c1=[-1;-4]
c2=[-2; -3]
dimX=length(c1)
dimU=length(c2)
b=[-2;-3]
A1=[
    1 -3;
   -1 -3
   ]
A2=[
    1 -2;
   -1 -1
   ]
M=1000

# Loading the necessary packages
#-------------------------------
using JuMP 
using GLPK
using MathOptInterface
const MOI = MathOptInterface

# Master Problem Description
# --------------------------

masterProblemModel = Model(optimizer=GLPK.GLPKOptimizerMIP()) 

# Variable Definition 
# ----------------------------------------------------------------
@variable(masterProblemModel, 0<= x[1:dimX]<= 1e6  , Int) 
@variable(masterProblemModel, t<=1e6)

# Objective Setting
# -----------------
@objective(masterProblemModel, Max, t)
iC=1 # iC is the iteration counter 

print(masterProblemModel)

# Trying the entire benders decomposition in one step
while(true)
    println("\n-----------------------")
    println("Iteration number = ", iC)
    println("-----------------------\n")
    println("The current master problem is")
    print(masterProblemModel)
     
    JuMP.optimize(masterProblemModel)
    t_status = JuMP.terminationstatus(masterProblemModel)# == MOI.Success
    p_status = JuMP.primalstatus(masterProblemModel)# == MOI.FeasiblePoint
    
    if p_status == MOI.InfeasiblePoint
        println("The problem is infeasible :-(")
        break
    end

    if t_status == MOI.InfeasibleOrUnbounded
        fmCurrent = M
        xCurrent=M*ones(dimX)
    end


    if p_status == MOI.FeasiblePoint
        fmCurrent = JuMP.resultvalue(t)
        xCurrent=Float64[]
            for i in 1:dimX
            push!(xCurrent, JuMP.resultvalue(x[i]))
        end
    end

    println("Status of the master problem is", t_status, 
            "\nwith fmCurrent = ", fmCurrent, 
            "\nxCurrent = ", xCurrent)


    subProblemModel = Model(optimizer=GLPK.GLPKOptimizerLP())

    cSub=b-A1*xCurrent

    @variable(subProblemModel, u[1:dimU]>=0)


    @constraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)],sum(A2[i,j]*u[i] for i in 1:size(A2,1))>=c2[j])
    # The second argument of @constraint macro, constrRefSubProblem[j=1:size(A2,2)] means that the j-th constraint is
    # referenced by constrRefSubProblem[j]. 

    
    @objective(subProblemModel, Min, dot(c1, xCurrent) + sum(cSub[i]*u[i] for i in 1:dimU))

    print("\nThe current subproblem model is \n", subProblemModel)

    JuMP.optimize(subProblemModel)  
    t_status_sub = JuMP.terminationstatus(subProblemModel)# == MOI.Success
    p_status_sub = JuMP.primalstatus(subProblemModel)# == MOI.FeasiblePoint

    fsxCurrent = JuMP.objectivevalue(subProblemModel) 

    uCurrent = Float64[]

    for i in 1:dimU
        push!(uCurrent, JuMP.resultvalue(u[i]))
    end

    γ=dot(b,uCurrent)

    println("Status of the subproblem is ", t_status_sub, 
        "\nwith fsxCurrent= ", fsxCurrent, 
        "\nand fmCurrent= ", fmCurrent) 
    
    if p_status_sub == MOI.FeasiblePoint &&  fsxCurrent == fmCurrent # we are done
        println("\n################################################")
        println("Optimal solution of the original problem found")
        println("The optimal objective value t is ", fmCurrent)
        println("The optimal x is ", xCurrent)
                println("The optimal v is ", JuMP.resultdual.(constrRefSubProblem))
        println("################################################\n")
        break
    end  
    
    if p_status_sub == MOI.FeasiblePoint && fsxCurrent < fmCurrent 
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv= A1'*uCurrent - c1
        @constraint(masterProblemModel, t+sum(cv[i]*x[i] for i in 1:dimX) <= γ )
        println("t + ", cv, "ᵀ x <= ", γ)
    end 
    
    if t_status_sub == MOI.InfeasibleOrUnbounded
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1'*uCurrent
        @constraint(masterProblemModel, sum(ce[i]*x[i] for i in 1:dimX) <= γ)
        println(ce, "ᵀ x <= ", γ)
    end
   
    iC=iC+1
                    sleep(0.5)
    
end



































                    # Loading the necessary packages
                    #-------------------------------
using JuMP 
#using CPLEX
using Gurobi

                    # Master Problem Description
                    # --------------------------

# Model name
 
#masterProblemModel = Model(solver = CplexSolver())
masterProblemModel = Model(optimizer = GurobiOptimizer(Heuristics=0, Cuts = 0, OutputFlag = 0)) # If we want to add Benders lazy constraints 
# in Gurobi, then we have to turn off Gurobi's own Cuts and Heuristics in the master problem

# Variable Definition (Only CplexSolver() works properly for these)
# ----------------------------------------------------------------
#@variable(masterProblemModel,  x[1:dimX] >=0 , Int) 
#@variable(masterProblemModel, t)

# ***************ALTERNATIVE VARIABLE DEFINITION FOR GUROBI************
#If we replace the two lines above with the follwoing:
@variable(masterProblemModel,  0<= x[1:dimX] <= 1e6 , Int)
@variable(masterProblemModel, t <= 1e6)
# then all the solvers give the expected solution
#**********************************************************************

# Objective Setting
# -----------------
@objective(masterProblemModel, Max, t)

print(masterProblemModel)

stringOfBenderCuts=AbstractString[] # this is an array of strings which will contain all the
# Benders cuts added to be displayed later

# There are no constraints when we start, so we will add all the constraints in the
# form of Benders cuts lazily

function addBendersCut(cb)
    #***************************************************************************
    # First we store the master problem solution in conventional data structures
    println("----------------------------")
    println("ITERATION NUMBER = ", length(stringOfBenderCuts)+1)
    println("---------------------------\n")
    
    fmCurrent = getvalue(t)
    xCurrent=Float64[]
    for i in 1:dimX
        push!(xCurrent, JuMP.resultvalue(x[i]))
    end
    
    # Display the current solution of the master problem
    println("MASTERPROBLEM INFORMATION")
    println("-------------------------")
    println("The master problem that was solved was:")
    print(masterProblemModel)
    println("with ", length(stringOfBenderCuts), " added lazy constraints")
    println(stringOfBenderCuts)
    println("Current Value of x is: ", xCurrent)
    println("Current objective value of master problem, fmCurrent is: ", fmCurrent)
    println("\n")
    
    #************************************************************************
    
    # ========================================================================
    #                         Now we solve the subproblem

    #  subProblemModel=Model(solver=CplexSolver())

    subProblemModel = Model(optimizer=GurobiOptimizer(Presolve=0, OutputFlag = 0))
    
    cSub=b-A1*xCurrent
    
    @variable(subProblemModel, u[1:dimU]>=0)
    

    @constraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)], sum(A2[i,j]*u[i] for i in 1:size(A2,1))>=c2[j])

    
    @objective(subProblemModel, Min, dot(c1, xCurrent) + sum(cSub[i]*u[i] for i in 1:dimU))
    
    println("The subproblem is being solved")
    
            statusSubProblem = JuMP.optimize(subProblemModel) 

    # We store the results achieved from the subproblem in conventional data structures    
    
    fsxCurrent = JuMP.objectivevalue(subProblemModel) 

    uCurrent = Float64[]
    for i in 1:dimU
                push!(uCurrent, JuMP.resultvalue(u[i]))
    end
    
    # Display the solution corresponding to the subproblem
    
    println("SUBPROBLEM INFORMATION")
    println("----------------------")
    println("The subproblem that was solved was: ")
    print(subProblemModel)
    println("Current status of the subproblem is ", statusSubProblem)
    println("Current Value of u is: ", uCurrent) # JuMP will return an extreme ray
    # automatically (if the solver supports it), so we do not need to change the syntax
    println("Current Value of fs(xCurrent) is: ", fsxCurrent)
    println("\n")
    
    # ==========================================================================
    
    
    
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Now we check the status of the algorithm and add Benders cut when necessary
    γ=dot(b,uCurrent)
      

    
    if statusSubProblem == :Optimal &&  fsxCurrent==fmCurrent # we are done
        println("OPTIMAL SOLUTION OF THE ORIGINAL PROBLEM FOUND :-)")
        println("The optimal objective value t is ", fmCurrent)
        println("The optimal x is ", xCurrent)
                println("The optimal v is ", JuMP.resultdual(constrRefSubProblem))
        println("\n")
        return
    end 

    println("-------------------ADDING LAZY CONSTRAINT----------------")  
        if statusSubProblem == :Optimal && fsxCurrent < fmCurrent 
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv= A1'*uCurrent - c1
        @lazyconstraint(cb, t+sum(cv[i]*x[i] for i in 1:dimX) <= γ)
        println("t + ", cv, "ᵀ x <= ", γ)
        push!(stringOfBenderCuts, string("t+", cv, "'x <=", γ))
    end
    
    if statusSubProblem == :Unbounded 
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1'*uCurrent
        @lazyconstraint(cb, sum(ce[i]*x[i] for i in 1:dimX) <= γ)
        println(ce, "x <= ", γ)
        push!(stringOfBenderCuts, string(ce, "ᵀ x <= ", γ))
    end
    println("\n") 
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        
end

addlazycallback(masterProblemModel, addBendersCut) # Telling the solver to use the 
# callback function

JuMP.optimize(masterProblemModel)