Benders Decomposition with callback

I am trying to implement Benders Decomposition for a graph problem, I have a master problem and subproblem and in each iteration, I have to solve the master problem change the RHS of the subproblem constraint and solve the subproblem and add lazy constraint until both problems converge, my problem that I can’t access the objective of the master problem with the following error

LoadError: MathOptInterface.OptimizeInProgress{MathOptInterface.TerminationStatus}: Cannot get the result as the the MOI.optimize! has not finished.
Moreover, to change the constraint I am deleting the constraint and creating new ones with the new values, there is any way to do the following in a more efficient way ??

using JuMP
using Random
import Gurobi
######################################
initial = 0
n = 32
s = 33
E = [(0, 1), (7, 9), (15, 27), (11, 1), (14, 9), (17, 16), (2, 33), (33, 17),
                  (1, 10), (7, 27), (15, 28), (19, 1), (15, 9), (18, 16), (33, 2), (18, 33),
                  (1, 11), (9, 14), (16, 17), (20, 1), (16, 9), (19, 16), (3, 33), (33, 18),
                  (1, 19), (9, 15), (16, 18), (3, 2), (17, 9), (20, 16), (33, 3), (19, 33),
                  (1, 20), (9, 16), (16, 19), (4, 2), (20, 9), (22, 16), (4, 33), (33, 19),
                  (2, 3), (9, 17), (16, 20), (5, 2), (21, 9), (25, 16), (33, 4), (20, 33),
                  (2, 4), (9, 20), (16, 22), (6, 2), (22, 9), (26, 16), (5, 33), (33, 20),
                  (2, 5), (9, 21), (16, 25), (7, 2), (23, 9), (27, 16), (33, 5), (21, 33),
                  (2, 6), (9, 22), (16, 26), (8, 2), (26, 9), (32, 16), (6, 33), (33, 21),
                  (2, 7), (9, 23), (16, 27), (9, 2), (27, 9), (23, 17), (33, 6), (22, 33),
                  (2, 8), (9, 26), (16, 32), (5, 3), (32, 9), (19, 18), (7, 33), (33, 22),
                  (2, 9), (9, 27), (17, 23), (9, 3), (11, 10), (20, 19), (33, 7), (23, 33),
                  (3, 5), (9, 32), (18, 19), (13, 3), (12, 10), (25, 19), (8, 33), (33, 23),
                  (3, 9), (10, 11), (19, 20), (14, 3), (13, 10), (26, 19), (33, 8), (24, 33),
                  (3, 13), (10, 12), (19, 25), (15, 3), (13, 12), (29, 19), (9, 33), (33, 24),
                  (3, 14), (10, 13), (19, 26), (16, 3), (26, 12), (30, 19), (33, 9), (25, 33),
                  (3, 15), (12, 13), (19, 29), (17, 3), (28, 12), (21, 20), (10, 33), (33, 25),
                  (3, 16), (12, 26), (19, 30), (9, 5), (19, 13), (28, 20), (33, 10), (26, 33),
                  (3, 17), (12, 28), (20, 21), (15, 5), (25, 13), (23, 22), (11, 33), (33, 26),
                  (5, 9), (13, 19), (20, 28), (23, 5), (26, 13), (24, 22), (33, 11), (27, 33),
                  (5, 15), (13, 25), (22, 23), (25, 5), (28, 13), (25, 22), (12, 33), (33, 27),
                  (5, 23), (13, 26), (22, 24), (29, 5), (29, 13), (26, 22), (33, 12), (28, 33),
                  (5, 25), (13, 28), (22, 25), (7, 6), (16, 14), (27, 22), (13, 33), (33, 28),
                  (5, 29), (13, 29), (22, 26), (9, 6), (16, 15), (27, 24), (33, 13), (29, 33),
                  (6, 7), (14, 16), (22, 27), (14, 6), (20, 15), (27, 26), (14, 33), (33, 29),
                  (6, 9), (15, 16), (24, 27), (21, 6), (21, 15), (30, 29), (33, 14), (30, 33),
                  (6, 14), (15, 20), (26, 27), (22, 6), (22, 15), (31, 29), (15, 33), (33, 30),
                  (6, 21), (15, 21), (29, 30), (26, 6), (24, 15), (0, 33), (33, 15), (31, 33),
                  (6, 22), (15, 22), (29, 31), (27, 6), (26, 15), (33, 0), (16, 33), (33, 31),
                  (6, 26), (15, 24), (1, 0), (9, 7), (27, 15), (1, 33), (33, 16), (32, 33),
                  (6, 27), (15, 26), (10, 1), (27, 7), (28, 15), (33, 1), (17, 33), (33, 32)]

E = sort(E)
###################################################
sub_model = direct_model(Gurobi.Optimizer())
set_optimizer_attribute(sub_model, "Threads", 1)
set_optimizer_attribute(sub_model, "presolve", 0)

@variable(sub_model, y[v=initial:n] >= 0)
@variable(sub_model,u[v=initial:s] >= 0)

x_H = Dict((i, j) => 0 for i in initial:s for j in initial:s if (i,j) in E)
x_H[(28, 15)] = 1
x_H[(15, s)] = 1
x_H[(s, 28)] = 1

for i = initial:n
    sub_c1=@constraint(sub_model, y[i] == sum(x_H[(i,j)] for j in initial:s if (i,j) in E))
    set_name(sub_c1, "sub_cons1_($(i))")
end
for (i, j) in E
    if j != s
            sub_c2=@constraint(sub_model,u[i] - u[j] <= n * (1 - x_H[(i, j)]) - 1)
            set_name(sub_c2, "sub_cons2_($(i),$(j))")
    end
end
sub_c3=@constraint(sub_model,u[s] == 1)
set_name(sub_c3, "sub_cons3")

@objective(sub_model, Min, -1*(sum(y[w] for w in initial:n)))
############################################
master_model = direct_model(Gurobi.Optimizer())
set_optimizer_attribute(master_model, "Threads", 1)
set_optimizer_attribute(master_model, "presolve", 0)

@variable(master_model, Z <= 0)
@variable(master_model,x[E],binary = true)

for (i, j) in E
    m_c1=@constraint(master_model,x[(i, j)] + x[(j, i)] <= 1)
    set_name(m_c1, "master_cons1_($(i),$(j))")
end

m_c2=@constraint(master_model,sum(x[(s, v)] for v in initial:n) == 1)
set_name(m_c2, "master_cons2")
m_c3=@constraint(master_model,sum(x[(v,s)] for v in initial:n) == 1)
set_name(m_c3, "master_cons3")

for (i, j) in E
    if i != s && j != s
        m_c4=@constraint(master_model,(sum(x[(i, g)] for g in initial:s
        if (g, i) in E && g != i && g != j)+sum(x[(j, g)] for g in initial:s
        if (g, j) in E && g != i && g != j)) <= 1)
          set_name(m_c4, "master_cons4_($(i),$(j))")
    end
end

for i in initial:n
    m_c5=@constraint(master_model,sum(x[(i, j)] for j in initial:s if (i, j) in E) <= 1)
    set_name(m_c5, "master_cons5_($(i))")

    m_c6=@constraint(master_model,sum(x[(j,i)] for j in initial:s if (j,i) in E) <= 1)
    set_name(m_c6, "master_cons6_($(i))")

    m_c7=@constraint(master_model,sum(x[(j, i)] for j in initial:s if (j, i) in E) ==
                    sum(x[(i, j)] for j in initial:s if (i, j) in E))
    set_name(m_c7, "master_con7_($(i))")
end
@objective(master_model, Min, Z)
################################################
gamma = Dict((1, v) => 0 for v in initial:n)
zeta = Dict((1,i, j) => 0 for i in initial:s for j in initial:n if (i,j) in E)
mu = Dict(1 => 0)
#################################################
optimize!(sub_model)
if termination_status(sub_model) == MOI.OPTIMAL
      k = 1
      for v in initial:n
           gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
       end
       for (i, j) in E
           if j != s
               zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
            end
      end
      mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))


      cut=@constraint(master_model, (sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
                                    +sum((n * (1 - x[(i, j)]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
                                    + mu[k]) <= Z)
      set_name(cut, "cut_($(k))")
else
       k = 1
       for v in initial:n
            gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
       end
       for (i, j) in E
            if j != s
                zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
             end
       end
       mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))

       cut=@constraint(master_model, (sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
                                     +sum((n * (1 - x[(i, j)]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
                                     + mu[k]) <= Z)
       set_name(cut, "cut_($(k))")

end
###################################################
function mycallback2(cb_data)
    model_status = callback_node_status(cb_data, master_model)
    if model_status == MOI.CALLBACK_NODE_STATUS_INTEGER
        Z_val = callback_value.(Ref(cb_data), Z)
        x_val = callback_value.(Ref(cb_data), x)
    end
    mp_solution =JuMP.objective_bound(master_model)
    status = getSubSolution(x_val,mp_solution)

    if status == 3
        add_feasibility_cut()
   elseif status == 2
        add_optimality_cut()
   else
       println("terminate")
       return
   end
end

function getSubSolution(x,mas_solution)

   for (a,b) in E
       if b != s
           delete(sub_model,constraint_by_name(sub_model,"sub_cons2_($(a),$(b))"))
           sub_c2=@constraint(sub_model,u[a] - u[b] <= n * (1 - x[(a, b)]) - 1)
           set_name(sub_c2, "sub_cons2_($(a),$(b))")
       end
   end
   println("change constraint")
   for a in initial:n
        temp = 0
        for b in initial:s
            if (a,b) in E
                temp+=x[(a,b)]
            end
        end
        delete(sub_model,constraint_by_name(sub_model,"sub_cons1_($(a))"))
        sub_c1=@constraint(sub_model, y[a] == temp)
        set_name(sub_c1, "sub_cons1_($(a))")
   end
   optimize!(sub_model)
   sub_solution = JuMP.objective_value(sub_model)
   if termination_status(sub_model) == MOI.OPTIMAL
       if abs(mas_solution) - abs(sub_solution) <= 0.0001
           println("True")
           return 1
       else
           return 2
       end
   elseif termination_status(sub_model) != MOI.OPTIMAL
          return 3
   else
       println("non of the above")
   end
end

function add_feasibility_cut()
    println("feas_cut")
    k = rand(Int, 1)
    for v in initial:n
         gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
     end
     for (i, j) in E
         if j != s
             zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
          end
    end
    mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))
    con = @build_constraint((sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
                      + sum((n * (1 - x[i, j]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
                      + mu[k]) <= 0)
    MOI.submit(cec_model, MOI.LazyConstraint(cb_data), con)

end

function add_optimality_cut()
    println("opt_cut")
    k = rand(Int, 1)
    for v in initial:n
         gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
     end
     for (i, j) in E
         if j != s
             zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
          end
    end
    mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))
    con = @build_constraint((sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
                      + sum((n * (1 - x[i, j]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
                      + mu[k]) <= Z)
    MOI.submit(cec_model, MOI.LazyConstraint(cb_data), con)

end

MOI.set(master_model, MOI.RawParameter("LazyConstraints"), 1)
MOI.set(master_model, MOI.LazyConstraintCallback(),mycallback2)
optimize!(master_model)