BnBTree.jl

include("table_log.jl")
include("bb_inits_and_defaults.jl")
include("bb_strategies.jl")
include("bb_user_limits.jl")
include("bb_type_correct.jl")
include("bb_integral_or_branch.jl")
include("bb_gains.jl")

"""
    get_branching_disc_idx!(m, step_obj, opts, disc2var_idx, gains, counter::Int64=1)

Get the index of a variable to branch on.
"""
function get_branching_disc_idx!(
    m,
    step_obj,
    opts,
    disc2var_idx,
    gains,
    counter::Int64 = 1,
)
    start = time()
    node = step_obj.node
    idx = 0
    strong_restarts = 0
    branch_strat = opts.branch_strategy
    status = :Normal
    if branch_strat == :MostInfeasible
        idx = branch_mostinfeasible(m, node, disc2var_idx)
        step_obj.branch_strategy = :MostInfeasible
    elseif branch_strat == :PseudoCost || branch_strat == :StrongPseudoCost
        if counter == 1 && branch_strat == :PseudoCost
            idx = branch_mostinfeasible(m, node, disc2var_idx)
            step_obj.branch_strategy = :MostInfeasible
        elseif counter <= opts.strong_branching_nsteps &&
               branch_strat == :StrongPseudoCost
            status, idx, strong_restarts =
                branch_strong!(m, opts, disc2var_idx, step_obj, counter)
            step_obj.branch_strategy = :Strong
        else
            idx = branch_pseudo(
                m,
                node,
                disc2var_idx,
                gains,
                opts.gain_mu,
                opts.atol,
            )
            step_obj.branch_strategy = :Pseudo
        end
    elseif branch_strat == :Reliability
        status, idx, strong_restarts =
            branch_reliable!(m, opts, step_obj, disc2var_idx, gains, counter)
        step_obj.branch_strategy = :Reliability
    end
    step_obj.state = status
    step_obj.var_idx = idx
    step_obj.nrestarts = strong_restarts
    step_obj.idx_time = time() - start
    return
end

"""
    process_node!(m, step_obj, cnode, disc2var_idx, temp)

Solve a child node `cnode` by relaxation.
Set the state and best_bound property.
Push integrals and new branch nodes to the step object
Return state
"""
function process_node!(m, step_obj, cnode, disc2var_idx, temp; restarts = 0)
    # set bounds
    for i in 1:m.num_var
        set_bounds(m.model, m.x[i], cnode.l_var[i], cnode.u_var[i])
    end
    if restarts > 0
        println("Doing one restart")
        restart_values = generate_random_restart(m)
        MOI.set(m.model, MOI.VariablePrimalStart(), m.x, restart_values)
    else
        MOI.set(m.model, MOI.VariablePrimalStart(), m.x, step_obj.node.solution)
    end

    old_mu_init =
        set_subsolver_option!(m.model, "Ipopt", "mu_init", 0.1 => 1e-5)

    MOI.optimize!(m.model)
    status = MOI.get(m.model, MOI.TerminationStatus())

    # reset mu_init
    set_subsolver_option!(m.model, "Ipopt", "mu_init", old_mu_init)

    objval = NaN
    cnode.solution = fill(NaN, m.num_var)
    if state_is_optimal(status; allow_almost = m.options.allow_almost_solved)
        objval = MOI.get(m.model, MOI.ObjectiveValue())
        cnode.solution = MOI.get(m.model, MOI.VariablePrimal(), m.x)
    end

    cnode.relaxation_state = status
    if !state_is_optimal(
        status;
        allow_almost = m.options.allow_almost_solved,
    ) && !state_is_infeasible(status)
        cnode.state = :Error
    elseif state_is_optimal(
        status;
        allow_almost = m.options.allow_almost_solved,
    )
        cnode.best_bound = objval
        set_cnode_state!(cnode, m, step_obj, disc2var_idx)
        if only_almost_solved(status) && (
            !(cnode.state == :Integral) ||
            m.options.allow_almost_solved_integral
        )
            @warn "Only almost solved"
        end

        # if almost_solved_integral is not allowed but it is integral => restart and hope for the best
        if cnode.state == :Integral &&
           only_almost_solved(status) &&
           !m.options.allow_almost_solved_integral
            if restarts >= 1
                cnode.state = :Almost_Solved
                return cnode.state
            else
                return process_node!(
                    m,
                    step_obj,
                    cnode,
                    disc2var_idx,
                    temp;
                    restarts = restarts + 1,
                )
            end
        end
        if !temp
            push_integral_or_branch!(step_obj, cnode)
        end
    else
        cnode.state = :Infeasible
    end
    return cnode.state
end

"""
    branch!(m, opts, step_obj, counter, disc2var_idx; temp=false)

Branch a node by using x[idx] <= floor(x[idx]) and x[idx] >= ceil(x[idx])
Solve both nodes and set current node state to done.
"""
function branch!(m, opts, step_obj, counter, disc2var_idx; temp = false)
    ps = opts.log_levels
    node = step_obj.node
    vidx = step_obj.var_idx

    start = time()
    # it might be already branched on
    if !temp && node.state != :Branch
        for cnode in [step_obj.l_nd, step_obj.r_nd]
            if cnode.state == :Branch || cnode.state == :Integral
                push_integral_or_branch!(step_obj, cnode)
            end
        end
        return step_obj.l_nd, step_obj.r_nd
    end

    l_nd_u_var = copy(node.u_var)
    r_nd_l_var = copy(node.l_var)
    l_nd_u_var[vidx] = floor(node.solution[vidx])
    r_nd_l_var[vidx] = ceil(node.solution[vidx])

    if opts.debug
        path_l = copy(node.path)
        path_r = copy(node.path)
        push!(path_l, node.hash)
        push!(path_r, node.hash)
    else
        path_l = []
        path_r = []
    end
    l_nd = new_left_node(node, l_nd_u_var; path = path_l)
    r_nd = new_right_node(node, r_nd_l_var; path = path_r)

    # save that this node branches on this particular variable
    node.var_idx = vidx

    check_print(ps, [:All, :FuncCall]) && println("branch")

    if !temp
        step_obj.l_nd = l_nd
        step_obj.r_nd = r_nd
        node.state = :Done
    end

    start_process = time()
    l_state = process_node!(m, step_obj, l_nd, disc2var_idx, temp)
    r_state = process_node!(m, step_obj, r_nd, disc2var_idx, temp)
    node_time = time() - start_process

    if !temp
        step_obj.node_branch_time += node_time
    end

    branch_strat = opts.branch_strategy

    if check_print(ps, [:All])
        println("State of left node: ", l_state)
        println("State of right node: ", r_state)
        println("l sol: ", l_nd.solution)
        println("r sol: ", r_nd.solution)
    end

    if !temp
        step_obj.branch_time += time() - start
    end

    return l_nd, r_nd
end

"""
    update_incumbent!(tree::BnBTreeObj, node::BnBNode)

Gets called if new integral solution was found.
Check whether it's a new incumbent and update if necessary
"""
function update_incumbent!(tree::BnBTreeObj, node::BnBNode)
    ps = tree.options.log_levels
    check_print(ps, [:All, :FuncCall]) && println("update_incumbent")

    factor = tree.obj_fac
    if !isdefined(tree, :incumbent) ||
       factor * node.best_bound > factor * tree.incumbent.objval
        objval = node.best_bound
        solution = copy(node.solution)
        status = node.relaxation_state
        @assert(
            state_is_optimal(status; allow_almost = false) || (
                only_almost_solved(status) &&
                tree.options.allow_almost_solved_integral
            )
        )
        tree.incumbent = Incumbent(objval, solution, only_almost_solved(status))
        if !tree.options.all_solutions
            bound!(tree)
        end
        return true
    end

    return false
end

"""
    bound!(tree::BnBTreeObj)
"""
function bound!(tree::BnBTreeObj)
    function isbetter(n)
        return f * n.best_bound > f * incumbent_val
    end
    incumbent_val = tree.incumbent.objval
    f = tree.obj_fac
    return filter!(isbetter, tree.branch_nodes)
end

"""
    add_incumbent_constr(m, incumbent)

Add a constraint >=/<= incumbent
"""
function add_incumbent_constr(m, incumbent)
    return add_obj_constraint(m, incumbent.objval)
end

"""
    add_obj_epsilon_constr(tree)

Add a constraint obj ≦ (1+ϵ)*LB or obj ≧ (1-ϵ)*UB
"""
function add_obj_epsilon_constr(tree)
    # add constr for objval
    if tree.options.obj_epsilon > 0
        ϵ = tree.options.obj_epsilon
        if tree.m.obj_sense == :Min
            rhs = (1 + ϵ) * tree.m.relaxation_objval
        else
            rhs = (1 - ϵ) * tree.m.relaxation_objval
        end
        add_obj_constraint(tree.m, rhs)
        tree.m.ncuts += 1
    end
end

"""
    get_next_branch_node!(tree)

Get the next branch node (sorted by best bound)
Return true,step_obj if there is a branch node and
false, nothing otherwise
"""
function get_next_branch_node!(tree)
    if length(tree.branch_nodes) == 0
        return false, false, nothing
    end
    bvalue, nidx =
        findmax([tree.obj_fac * n.best_bound for n in tree.branch_nodes])

    trav_strat = tree.options.traverse_strategy
    if trav_strat == :BFS
        # use the one with highest depth
        _value, nidx = findmax([
            if isapprox(tree.obj_fac * n.best_bound, bvalue)
                n.level
            else
                -Inf
            end for n in tree.branch_nodes
        ])
    end

    if trav_strat == :DFS ||
       (trav_strat == :DBFS && !isdefined(tree, :incumbent))
        _value, nidx = findmax([n.level for n in tree.branch_nodes])
    end

    tree.best_bound = tree.obj_fac * bvalue
    bbreak = isbreak_mip_gap(tree)
    if bbreak
        return false, true, nothing
    end

    node = tree.branch_nodes[nidx]
    deleteat!(tree.branch_nodes, nidx)

    return true, false, new_default_step_obj(tree.m, node)
end

"""
    one_branch_step!(m1, incumbent, opts, step_obj, disc2var_idx, gains, counter)

Get a branch variable using the specified strategy and branch on the node in step_obj
using that variable. Return the new updated step_obj
"""
function one_branch_step!(
    m1,
    incumbent,
    opts,
    step_obj,
    disc2var_idx,
    gains,
    counter,
)
    if m1 === nothing
        global m
        global is_newincumbent
        if opts.incumbent_constr && incumbent !== nothing && is_newincumbent
            is_newincumbent = false
            add_incumbent_constr(m, incumbent)
        end
    else
        m = m1
    end

    node = step_obj.node
    step_obj.counter = counter

    # get branch variable
    node_idx_start = time()
    get_branching_disc_idx!(m, step_obj, opts, disc2var_idx, gains, counter)
    step_obj.node_idx_time = time() - node_idx_start
    # if no variable got selected might be true that all variables are already type correct
    # node.solution can be updated if one child is infeasible and the other optimal (can result in discrete)
    if step_obj.var_idx == 0 && are_type_correct(
        step_obj.node.solution,
        m.var_type,
        disc2var_idx,
        opts.atol,
    )
        push!(step_obj.integral, node)
    else
        if step_obj.state != :Infeasible && step_obj.state != :PartlyInfeasible
            @assert step_obj.var_idx != 0
            branch!(m, opts, step_obj, counter, disc2var_idx)
        end
    end
    return step_obj
end

"""
    upd_time_obj!(time_obj, step_obj)

Add step_obj times to time_obj
"""
function upd_time_obj!(time_obj, step_obj)
    time_obj.solve_leaves_get_idx += step_obj.node_idx_time
    time_obj.solve_leaves_branch += step_obj.node_branch_time
    time_obj.branch += step_obj.branch_time
    time_obj.get_idx += step_obj.idx_time
    return time_obj.upd_gains += step_obj.upd_gains_time
end

"""
    upd_tree_obj!(tree, step_obj, time_obj)

Update the tree obj like new incumbent or new branch nodes using the step_obj
Return true if it's the end of the algorithm (checking different break rules)
"""
function upd_tree_obj!(tree, step_obj, time_obj)
    node = step_obj.node
    still_running = true

    if step_obj.node.level + 1 > tree.m.nlevels
        tree.m.nlevels = step_obj.node.level + 1
    end

    if step_obj.state == :Infeasible
        # if there is no incumbent yet
        if !isdefined(tree, :incumbent)
            # this can be locally infeasible or globally infeasible depending on the solver
            tree.limit = :Infeasible
        end # it will terminate and use the current solution as optimal (might want to rerun as an option)
        still_running = false
    end

    if still_running
        upd_gains_step!(tree, step_obj)
    end

    bbreak = upd_integral_branch!(tree, step_obj)
    if bbreak
        still_running = false
    end

    upd_time_obj!(time_obj, step_obj)
    return !still_running
end

"""
    solve_sequential(tree,
        last_table_arr,
        time_bnb_solve_start,
        fields,
        field_chars,
        time_obj)

Run branch and bound on a single processor
"""
function solve_sequential(
    tree,
    last_table_arr,
    time_bnb_solve_start,
    fields,
    field_chars,
    time_obj,
)
    m = tree.m
    opts = tree.options
    disc2var_idx = tree.disc2var_idx
    counter = 0
    ps = tree.options.log_levels

    dictTree = Dict{Any,Any}()
    while true
        # the _ is only needed for parallel
        exists, _, step_obj = get_next_branch_node!(tree)
        !exists && break
        isbreak_after_step!(tree) && break
        counter += 1
        if isdefined(tree, :incumbent)
            step_obj = one_branch_step!(
                m,
                tree.incumbent,
                opts,
                step_obj,
                disc2var_idx,
                tree.obj_gain,
                counter,
            )
        else
            step_obj = one_branch_step!(
                m,
                nothing,
                opts,
                step_obj,
                disc2var_idx,
                tree.obj_gain,
                counter,
            )
        end
        m.nnodes += 2 # two nodes explored per branch
        node = step_obj.node

        bbreak = upd_tree_obj!(tree, step_obj, time_obj)
        tree.options.debug &&
            (dictTree = push_step2treeDict!(dictTree, step_obj))

        if check_print(ps, [:Table])
            last_table_arr = print_table(
                1,
                tree,
                node,
                step_obj,
                time_bnb_solve_start,
                fields,
                field_chars;
                last_arr = last_table_arr,
            )
        end

        if bbreak
            break
        end
    end

    tree.options.debug && (tree.m.debugDict[:tree] = dictTree)
    return counter
end

function sendto(p::Int; args...)
    for (nm, val) in args
        remotecall_fetch(Core.eval, p, Juniper, Expr(:(=), nm, val))
    end
end

function dummysolve()
    global m
    return optimize!(m.model)
end

"""
    pmap(f, tree, last_table_arr, time_bnb_solve_start,
        fields, field_chars, time_obj)

Run the solving steps on several processors
"""
function pmap(
    f,
    tree,
    last_table_arr,
    time_bnb_solve_start,
    fields,
    field_chars,
    time_obj,
)
    np = nprocs()  # determine the number of processes available
    if tree.options.processors + 1 < np
        np = tree.options.processors + 1
    end

    # function to produce the next work item from the queue.
    # in this case it's just an index.
    ps = tree.options.log_levels
    still_running = true
    run_counter = 0
    counter = 0
    for p in 2:np
        seed = i -> Random.seed!(JUNIPER_RNG, i)
        remotecall_fetch(seed, p, 1)
    end
    @sync begin
        for p in 2:np
            @async begin
                sendto(p, m = tree.m)
                sendto(p, is_newincumbent = false)
            end
        end
    end

    for p in 3:np
        remotecall(dummysolve, p)
    end

    proc_counter = zeros(np)

    dictTree = Dict{Any,Any}()

    branch_strat = tree.options.branch_strategy
    opts = tree.options
    @sync begin
        for p in 1:np
            if p != myid() || np == 1
                @async begin
                    while true
                        exists, bbreak, step_obj = get_next_branch_node!(tree)
                        if bbreak
                            still_running = false
                            break
                        end

                        while !exists && still_running
                            sleep(0.1)
                            exists, bbreak, step_obj =
                                get_next_branch_node!(tree)
                            exists && break
                            if bbreak
                                still_running = false
                                break
                            end
                        end
                        !still_running && break

                        if isbreak_after_step!(tree)
                            still_running = false
                            break
                        end

                        run_counter += 1
                        counter += 1
                        if isdefined(tree, :incumbent)
                            step_obj = remotecall_fetch(
                                f,
                                p,
                                nothing,
                                tree.incumbent,
                                tree.options,
                                step_obj,
                                tree.disc2var_idx,
                                tree.obj_gain,
                                counter,
                            )
                        else
                            step_obj = remotecall_fetch(
                                f,
                                p,
                                nothing,
                                nothing,
                                tree.options,
                                step_obj,
                                tree.disc2var_idx,
                                tree.obj_gain,
                                counter,
                            )
                        end
                        tree.m.nnodes += 2 # two nodes explored per branch
                        run_counter -= 1

                        !still_running && break

                        bbreak = upd_tree_obj!(tree, step_obj, time_obj)
                        tree.options.debug &&
                            (dictTree = push_step2treeDict!(dictTree, step_obj))

                        if run_counter == 0 && length(tree.branch_nodes) == 0
                            still_running = false
                        end

                        bbreak && (still_running = false)

                        if check_print(ps, [:Table])
                            if length(tree.branch_nodes) > 0
                                bvalue, nidx = findmax([
                                    tree.obj_fac * n.best_bound for
                                    n in tree.branch_nodes
                                ])
                                tree.best_bound = tree.obj_fac * bvalue
                            end
                            last_table_arr = print_table(
                                p,
                                tree,
                                step_obj.node,
                                step_obj,
                                time_bnb_solve_start,
                                fields,
                                field_chars;
                                last_arr = last_table_arr,
                            )
                        end

                        !still_running && break
                    end
                end
            end
        end
    end
    tree.options.debug && (tree.m.debugDict[:tree] = dictTree)
    return counter
end

"""
    solvemip(tree::BnBTreeObj)

Solve the MIP part of a problem given by BnBTreeObj using branch and bound.
 - Identify the node to branch on
 - Get variable to branch on
 - Solve subproblems
"""
function solvemip(tree::BnBTreeObj)
    time_obj = init_time_obj()
    time_bnb_solve_start = time()

    ps = tree.options.log_levels

    # check if already integral
    if are_type_correct(
        tree.m.relaxation_solution,
        tree.m.var_type,
        tree.disc2var_idx,
        tree.options.atol,
    )
        tree.nsolutions = 1
        objval = MOI.get(tree.m.model, MOI.ObjectiveValue())
        sol = MOI.get(tree.m.model, MOI.VariablePrimal(), tree.m.x)
        bbound = try
            MOI.get(tree.m.model, MOI.ObjectiveBound())
        catch
            MOI.get(tree.m.model, MOI.ObjectiveValue())
        end
        tree.incumbent =
            Incumbent(objval, sol, only_almost_solved(tree.m.status))
        return tree.incumbent
    end

    # check if incumbent and if mip gap already fulfilled (i.e found by fp)
    if isbreak_mip_gap(tree)
        return tree.incumbent
    end

    last_table_arr = []
    fields = []
    field_chars = []
    # Print table init
    if check_print(ps, [:Table])
        fields, field_chars = get_table_config(tree.options)
        print_table_header(fields, field_chars)
    end

    check_print(ps, [:All, :FuncCall]) && println("Solve Tree")

    branch_strat = tree.options.branch_strategy

    add_obj_epsilon_constr(tree)

    # use pmap if more then one processor
    if tree.options.processors > 1 || tree.options.force_parallel
        counter = pmap(
            Juniper.one_branch_step!,
            tree,
            last_table_arr,
            time_bnb_solve_start,
            fields,
            field_chars,
            time_obj,
        )
    else
        counter = solve_sequential(
            tree,
            last_table_arr,
            time_bnb_solve_start,
            fields,
            field_chars,
            time_obj,
        )
    end

    if !isdefined(tree, :incumbent) && tree.limit == :None
        # infeasible
        tree.limit = :Infeasible
    end

    if tree.options.obj_epsilon != 0 && tree.limit == :Infeasible
        @warn "Maybe only infeasible because of obj_epsilon."
    end

    if !isnan(tree.options.best_obj_stop) && isdefined(tree, :incumbent)
        inc_val = tree.incumbent.objval
        bos = tree.options.best_obj_stop
        sense = tree.m.obj_sense
        if (sense == :Min && inc_val > bos) || (sense == :Max && inc_val < bos)
            @warn "best_obj_gap couldn't be reached."
        end
    end

    if length(tree.branch_nodes) > 0
        bvalue, nidx =
            findmax([tree.obj_fac * n.best_bound for n in tree.branch_nodes])
        tree.best_bound = tree.obj_fac * bvalue
    else
        if isdefined(tree, :incumbent) && !isnan(tree.incumbent.objval)
            tree.best_bound = tree.incumbent.objval
        end
    end

    tree.m.nbranches = counter

    time_bnb_solve = time() - time_bnb_solve_start
    check_print(ps, [:Table]) && println("")
    check_print(ps, [:All, :Info]) && println("#branches: ", counter)

    tree.options.debug && debug_set_tree_obj_gain!(tree)
    return check_print(ps, [:All, :Timing]) &&
           print_final_timing(time_bnb_solve, time_obj)
end