5#ifndef IDOL_OPTIMIZERS_BRANCHANDBOUND_H
6#define IDOL_OPTIMIZERS_BRANCHANDBOUND_H
8#include "idol/mixed-integer/optimizers/branch-and-bound/nodes/NodeSet.h"
9#include "idol/general/optimizers/Algorithm.h"
10#include "idol/mixed-integer/optimizers/branch-and-bound/branching-rules/factories/BranchingRuleFactory.h"
11#include "idol/mixed-integer/optimizers/branch-and-bound/node-selection-rules/factories/NodeSelectionRuleFactory.h"
12#include "idol/mixed-integer/optimizers/branch-and-bound/branching-rules/impls/BranchingRule.h"
13#include "idol/mixed-integer/optimizers/branch-and-bound/node-selection-rules/impls/NodeSelectionRule.h"
14#include "idol/mixed-integer/optimizers/branch-and-bound/nodes/Node.h"
15#include "idol/mixed-integer/optimizers/branch-and-bound/nodes/NodeUpdator.h"
16#include "idol/mixed-integer/optimizers/branch-and-bound/callbacks/AbstractBranchAndBoundCallbackI.h"
17#include "idol/general/optimizers/OptimizerFactory.h"
18#include "idol/mixed-integer/modeling/models/Model.h"
19#include "idol/mixed-integer/optimizers/branch-and-bound/logs/Factory.h"
25namespace idol::Optimizers {
26 template<
class NodeInfoT>
class BranchAndBound;
29template<
class NodeInfoT>
34 std::unique_ptr<OptimizerFactory> m_relaxation_optimizer_factory;
36 unsigned int m_n_threads = 1;
37 std::vector<std::unique_ptr<Model>> m_relaxations;
38 std::vector<std::unique_ptr<NodeUpdator<NodeInfoT>>> m_node_updators;
40 std::unique_ptr<BranchingRule<NodeInfoT>> m_branching_rule;
41 std::unique_ptr<NodeSelectionRule<NodeInfoT>> m_node_selection_rule;
42 std::unique_ptr<typename Logs::BranchAndBound::Factory<NodeInfoT>::Strategy> m_logger;
44 std::unique_ptr<AbstractBranchAndBoundCallbackI<NodeInfoT>> m_callback;
46 bool m_has_integer_objective =
false;
47 std::vector<unsigned int> m_steps = { std::numeric_limits<unsigned int>::max(), 0, 0 };
48 unsigned int m_n_created_nodes = 0;
49 unsigned int m_n_solved_nodes = 0;
50 unsigned int m_n_active_nodes = 0;
51 double m_root_node_best_bound = -Inf;
52 double m_root_node_best_obj = +Inf;
54 std::optional<TreeNode> m_incumbent;
56 void build()
override;
57 void hook_before_optimize()
override;
58 void hook_optimize()
override;
59 void hook_after_optimize()
override;
60 void add(
const Var &t_var)
override;
61 void add(
const Ctr &t_ctr)
override;
62 void add(
const QCtr &t_ctr)
override;
63 void remove(
const Var &t_var)
override;
64 void remove(
const Ctr &t_ctr)
override;
65 void remove(
const QCtr &t_ctr)
override;
66 void update()
override;
67 void write(
const std::string &t_name)
override;
69 void detect_integer_objective();
70 void create_relaxations();
73 void analyze(
const TreeNode& t_node,
unsigned int t_relaxation_id,
bool* t_explore_children_flag,
bool* t_reoptimize_flag);
75 std::vector<TreeNode> create_child_nodes(
const TreeNode& t_node);
77 void set_as_incumbent(
const TreeNode& t_node);
78 [[nodiscard]]
bool gap_is_closed()
const;
81 bool is_valid(
const TreeNode& t_node)
const;
84 void log_after_termination();
87 void update_obj_sense()
override;
88 void update_obj()
override;
89 void update_rhs()
override;
90 void update_obj_constant()
override;
91 void update_mat_coeff(
const Ctr &t_ctr,
const Var &t_var)
override;
92 void update_ctr_type(
const Ctr &t_ctr)
override;
93 void update_ctr_rhs(
const Ctr &t_ctr)
override;
94 void update_var_type(
const Var &t_var)
override;
95 void update_var_lb(
const Var &t_var)
override;
96 void update_var_ub(
const Var &t_var)
override;
98 void update_var_obj(
const Var &t_var)
override;
99 void set_status(SolutionStatus t_status)
override;
100 void set_reason(SolutionReason t_reason)
override;
101 void set_best_bound(
double t_value)
override;
102 void set_best_obj(
double t_value)
override;
111 [[nodiscard]] std::string name()
const override {
return "Branch-and-Bound"; }
113 void solve(
TreeNode& t_node,
unsigned int t_relaxation_id)
const;
115 virtual void set_subtree_depth(
unsigned int t_depth) { m_steps.at(1) = t_depth; }
117 [[nodiscard]]
unsigned int subtree_depth()
const {
return m_steps.at(1); }
121 void submit_heuristic_solution(NodeInfoT* t_info);
123 void submit_lower_bound(
double t_lower_bound);
125 [[nodiscard]]
unsigned int n_created_nodes()
const {
return m_n_created_nodes; }
127 [[nodiscard]]
unsigned int n_solved_nodes()
const {
return m_n_solved_nodes; }
129 [[nodiscard]]
unsigned int n_active_nodes()
const {
return m_n_active_nodes; }
131 [[nodiscard]]
const Model& relaxation()
const {
return *m_relaxations.front(); }
133 [[nodiscard]]
Model& relaxation() {
return *m_relaxations.front(); }
135 [[nodiscard]]
double root_node_best_bound()
const {
return m_root_node_best_bound; }
137 [[nodiscard]]
double root_node_best_obj()
const {
return m_root_node_best_obj; }
143 [[nodiscard]]
bool has_incumbent()
const {
return m_incumbent.has_value(); }
145 const TreeNode& incumbent()
const {
return m_incumbent.value(); }
147 [[nodiscard]]
double get_var_primal(
const Var &t_var)
const override;
149 [[nodiscard]]
double get_var_reduced_cost(
const Var &t_var)
const override;
151 [[nodiscard]]
double get_var_ray(
const Var &t_var)
const override;
153 [[nodiscard]]
double get_ctr_dual(
const Ctr &t_ctr)
const override;
155 [[nodiscard]]
double get_ctr_farkas(
const Ctr &t_ctr)
const override;
157 [[nodiscard]]
unsigned int get_n_solutions()
const override;
159 [[nodiscard]]
unsigned int get_solution_index()
const override;
161 void set_solution_index(
unsigned int t_index)
override;
163 void terminate()
override;
166template<
class NodeInfoT>
168 for (
auto& relaxation : m_relaxations) {
169 relaxation->remove(t_ctr);
173template<
class NodeInfoT>
175 for (
auto& relaxation : m_relaxations) {
176 relaxation->add(t_ctr);
180template<
class NodeInfoT>
184 result = m_branching_rule->is_valid(t_node);
188template<
class NodeInfoT>
191 Optimizer::terminate();
194template<
class NodeInfoT>
196 const auto& src_model = parent();
197 const auto& objective = src_model.get_obj_expr();
199 if (!is_integer(objective.affine().constant(), Tolerance::Integer)) {
200 m_has_integer_objective =
false;
204 for (
const auto& [var, val] : objective.affine().linear()) {
205 if (src_model.get_var_type(var) == Continuous || !is_integer(val, Tolerance::Integer)) {
206 m_has_integer_objective =
false;
211 m_has_integer_objective =
true;
214template<
class NodeInfoT>
216 if (m_n_solved_nodes > 1) {
217 throw Exception(
"Reduced cost not available.");
219 return m_relaxations.front()->get_var_reduced_cost(t_var);
222template<
class NodeInfoT>
225 if (!get_param_logs()) {
230 m_callback->log_after_termination();
233 m_logger->log_after_termination();
237template<
class NodeInfoT>
239 m_logger->log_node_after_solve(t_node);
242template<
class NodeInfoT>
245 throw Exception(
"Solution index out of bounds.");
249template<
class NodeInfoT>
254template<
class NodeInfoT>
256 return m_incumbent.has_value();
259template<
class NodeInfoT>
261 for (
auto& relaxation : m_relaxations) {
262 relaxation->set_obj_expr(parent().get_obj_expr());
266template<
class NodeInfoT>
268 for (
auto& relaxation : m_relaxations) {
269 relaxation->set_rhs_expr(parent().get_rhs_expr());
273template<
class NodeInfoT>
275 for (
auto& relaxation : m_relaxations) {
276 relaxation->set_obj_const(parent().get_obj_expr().affine().constant());
280template<
class NodeInfoT>
282 for (
auto& relaxation : m_relaxations) {
283 relaxation->set_mat_coeff(t_ctr, t_var, parent().get_mat_coeff(t_ctr, t_var));
287template<
class NodeInfoT>
289 for (
auto& relaxation : m_relaxations) {
290 relaxation->set_ctr_type(t_ctr, parent().get_ctr_type(t_ctr));
294template<
class NodeInfoT>
296 for (
auto& relaxation : m_relaxations) {
297 relaxation->set_ctr_rhs(t_ctr, parent().get_ctr_rhs(t_ctr));
301template<
class NodeInfoT>
303 for (
auto& relaxation : m_relaxations) {
304 relaxation->set_var_type(t_var, parent().get_var_type(t_var));
308template<
class NodeInfoT>
310 for (
auto& relaxation : m_relaxations) {
311 relaxation->set_var_lb(t_var, parent().get_var_lb(t_var));
315template<
class NodeInfoT>
317 for (
auto& relaxation : m_relaxations) {
318 relaxation->set_var_ub(t_var, parent().get_var_ub(t_var));
322template<
class NodeInfoT>
324 for (
auto& relaxation : m_relaxations) {
325 relaxation->set_var_obj(t_var, parent().get_var_obj(t_var));
329template<
class NodeInfoT>
331 throw Exception(
"Not implemented");
334template<
class NodeInfoT>
336 if (m_n_solved_nodes > 1) {
337 throw Exception(
"Accessing Farkas certificate for MIP is not available after the root node.");
339 return m_relaxations.front()->get_ctr_farkas(t_ctr);
342template<
class NodeInfoT>
344 if (m_n_solved_nodes > 1) {
345 throw Exception(
"Accessing dual values for MIP is not available after the root node.");
347 return m_relaxations.front()->get_ctr_dual(t_ctr);
350template<
class NodeInfoT>
352 if (m_n_solved_nodes > 1) {
353 throw Exception(
"Ray not available.");
355 return m_relaxations.front()->get_var_ray(t_var);
358template<
class NodeInfoT>
361 if (!m_incumbent.has_value()){
362 throw Exception(
"Trying to access primal values while no incumbent was found.");
365 return m_incumbent->info().primal_solution().get(t_var);
368template<
class NodeInfoT>
370 if (t_lower_bound > get_best_obj()) {
372 set_reason(NotSpecified);
376 if (t_lower_bound > get_best_bound()) {
377 set_best_bound(t_lower_bound);
382template<
class NodeInfoT>
385 auto t_node = Node<NodeInfoT>::create_detached_node(t_info);
387 if (t_node.info().objective_value() < get_best_obj()) {
390 set_as_incumbent(t_node);
402template<
class NodeInfoT>
410 return m_callback->operator()(
this, t_event, t_node, m_relaxations[t_relaxation_id].get());
413template<
class NodeInfoT>
415 m_callback->add_callback(t_cb);
418template<
class NodeInfoT>
423template<
class NodeInfoT>
425 const OptimizerFactory& t_node_optimizer,
426 const BranchingRuleFactory<NodeInfoT>& t_branching_rule_factory,
427 const NodeSelectionRuleFactory<NodeInfoT>& t_node_selection_rule_factory,
428 AbstractBranchAndBoundCallbackI<NodeInfoT>* t_callback,
429 const Logs::BranchAndBound::Factory<NodeInfoT>& t_logger_factory)
430 : Algorithm(t_model),
431 m_relaxation_optimizer_factory(t_node_optimizer.clone()),
432 m_branching_rule(t_branching_rule_factory(*this)),
433 m_node_selection_rule(t_node_selection_rule_factory(*this)),
434 m_callback(t_callback),
435 m_logger(t_logger_factory.operator()(*this)) {
437 create_relaxations();
441template<
class NodeInfoT>
444 const auto &original_model = parent();
446 m_relaxations.clear();
447 m_relaxations.reserve(m_n_threads);
449 m_node_updators.clear();
450 m_node_updators.reserve(m_n_threads);
452 for (
unsigned int i = 0 ; i < m_n_threads ; ++i) {
453 auto* relaxation = original_model.clone();
454 relaxation->use(*m_relaxation_optimizer_factory);
455 m_relaxations.emplace_back(relaxation);
456 m_node_updators.emplace_back(
dynamic_cast<NodeUpdator<NodeInfoT>*
>(NodeInfoT::create_updator(parent(), *relaxation)));
461template<
class NodeInfoT>
464 auto root_node = Node<NodeInfoT>::create_root_node();
465 assert(root_node.id() == 0);
472template<
class NodeInfoT>
474 Algorithm::hook_before_optimize();
477 set_best_bound(std::max(-Inf, get_param_best_obj_stop()));
478 set_best_obj(std::min(+Inf, get_param_best_bound_stop()));
481 detect_integer_objective();
483 m_n_created_nodes = 0;
484 m_n_solved_nodes = 0;
485 m_n_active_nodes = 0;
487 m_root_node_best_bound = -Inf;
488 m_root_node_best_obj = Inf;
490 m_branching_rule->initialize();
491 for (
auto& node_updator : m_node_updators) {
492 node_updator->initialize();
494 m_callback->initialize(
this);
496 m_logger->initialize();
499template<
class NodeInfoT>
502 auto root_node = create_root_node();
504 SetOfActiveNodes active_nodes;
506 explore(root_node, 0, active_nodes, 0);
508 if (active_nodes.empty()) {
509 set_best_bound(get_best_obj());
512 m_node_updators[0]->clear();
514 if (get_status() == Fail) {
518 if (!m_incumbent.has_value()) {
520 if (is_pos_inf(get_best_obj())) {
521 set_status(Infeasible);
525 if (is_neg_inf(get_best_obj())) {
526 set_status(Unbounded);
532 set_status(Feasible);
534 if (gap_is_closed()) {
542template<
class NodeInfoT>
545 Optimizer::hook_after_optimize();
547 log_after_termination();
551template<
class NodeInfoT>
554 unsigned int t_relaxation_id,
555 SetOfActiveNodes & t_active_nodes,
556 unsigned int t_step) {
558 if (is_terminated() || gap_is_closed()) {
return; }
560 bool reoptimize_flag;
561 bool explore_children_flag;
565 solve(t_node, t_relaxation_id);
567 analyze(t_node, t_relaxation_id, &explore_children_flag, &reoptimize_flag);
569 log_node_after_solve(t_node);
571 }
while (reoptimize_flag);
576 if (is_terminated() || gap_is_closed()) {
return; }
578 if (!explore_children_flag) {
return; }
580 t_active_nodes.emplace(t_node);
582 const unsigned int max_levels = m_steps.at(t_step);
584 for (
unsigned int level = 0 ; level < max_levels ; ++level) {
586 prune_nodes_by_bound(t_active_nodes);
589 update_lower_bound(t_active_nodes);
590 m_n_active_nodes = t_active_nodes.size();
593 if (t_active_nodes.empty()) {
break; }
595 if (is_terminated() || gap_is_closed()) {
break; }
597 auto selected_node = select_node_for_branching(t_active_nodes);
599 auto children = create_child_nodes(selected_node);
601 const unsigned int n_children = children.size();
603 std::vector<SetOfActiveNodes> active_nodes(n_children);
605 if (m_n_threads > 1 && t_step == 0) {
607 std::vector<std::vector<unsigned int>> to_explore(m_n_threads);
609 for (
unsigned int q = 0 ; q < n_children ; ++q) {
610 const unsigned int relaxation_id = q % m_n_threads;
611 to_explore[relaxation_id].emplace_back(q);
614#pragma omp parallel for num_threads(m_n_threads)
615 for (
unsigned int k = 0 ; k < m_n_threads ; ++k) {
616 for (
unsigned int j = 0, n_tasks = to_explore[k].size() ; j < n_tasks ; ++j) {
617 const unsigned int q = to_explore[k][j];
618 explore(children[q], k, active_nodes[q], t_step + 1);
624 for (
unsigned int q = 0 ; q < n_children ; ++q) {
625 explore(children[q], t_relaxation_id, active_nodes[q], t_step + 1);
630 for (
unsigned int q = 0 ; q < n_children ; ++q) {
631 backtrack(t_active_nodes, active_nodes[q]);
638template<
class NodeInfoT>
640 unsigned int t_relaxation_id)
const {
642 auto& node_updator = *m_node_updators[t_relaxation_id];
643 auto& relaxation = *m_relaxations[t_relaxation_id];
645 relaxation.optimizer().set_param_best_bound_stop(std::min(get_best_obj(), get_param_best_bound_stop()));
646 relaxation.optimizer().set_param_time_limit(get_remaining_time());
648 node_updator.prepare(t_node);
660 relaxation.optimize();
662 t_node.info().save(parent(), relaxation);
664 m_branching_rule->on_node_solved(t_node);
668template<
class NodeInfoT>
671 *t_explore_children_flag =
false;
672 *t_reoptimize_flag =
false;
674 if (get_best_bound() > get_best_obj()) {
676 set_reason(NotSpecified);
681 const auto& status = t_node.info().status();
682 const auto& reason = t_node.info().reason();
684 if (get_remaining_time() == 0 || reason == TimeLimit) {
685 set_reason(TimeLimit);
690 if (status == Unbounded) {
697 if (status == Infeasible || status == InfOrUnbnd) {
699 if (t_node.level() == 0) {
703 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
707 if (status == Feasible && reason == ObjLimit) {
708 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
712 if (status == Fail || status == Loaded) {
715 set_reason(NotSpecified);
717 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
721 if (t_node.level() == 0) {
722 m_root_node_best_bound = t_node.info().objective_value();
725 if (t_node.info().objective_value() < get_best_obj()) {
727 if (is_valid(t_node)) {
729 auto side_effects = call_callbacks(IncumbentSolution, t_node, t_relaxation_id);
731 if (side_effects.n_added_user_cuts > 0 || side_effects.n_added_lazy_cuts > 0) {
732 *t_reoptimize_flag =
true;
736 set_as_incumbent(t_node);
743 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
748 auto side_effects = call_callbacks(InvalidSolution, t_node, t_relaxation_id);
750 if (t_node.level() == 0) {
751 m_root_node_best_obj = get_best_obj();
754 if (side_effects.n_added_lazy_cuts > 0 || side_effects.n_added_user_cuts > 0 || side_effects.n_added_local_variable_branching > 0) {
755 *t_reoptimize_flag =
true;
759 *t_explore_children_flag =
true;
763template<
class NodeInfoT>
766 m_incumbent = t_node;
767 set_best_obj(t_node.info().objective_value());
768 set_status(Feasible);
771template<
class NodeInfoT>
774 if (t_active_nodes.empty()) {
return; }
776 auto& lowest_node = *t_active_nodes.by_objective_value().begin();
777 const double raw_lower_bound = lowest_node.info().objective_value();
778 const double lower_bound = m_has_integer_objective && !is_integer(raw_lower_bound,
Tolerance::Integer) ? std::ceil(raw_lower_bound) : raw_lower_bound;
779 if (lower_bound > get_best_bound()) {
780 set_best_bound(lower_bound);
785template<
class NodeInfoT>
788 const double upper_bound = get_best_obj();
790 auto it = t_active_nodes.by_objective_value().begin();
791 auto end = t_active_nodes.by_objective_value().end();
795 const auto& node = *it;
796 const double raw_lower_bound = node.info().objective_value();
797 const double lower_bound = m_has_integer_objective && !is_integer(raw_lower_bound,
Tolerance::Integer) ? std::ceil(raw_lower_bound) : raw_lower_bound;
799 if (lower_bound >= upper_bound - get_tol_mip_absolute_gap()) {
800 it = t_active_nodes.erase(it);
801 end = t_active_nodes.by_objective_value().end();
810template<
class NodeInfoT>
813 auto iterator = m_node_selection_rule->operator()(t_active_nodes);
814 auto result = *iterator;
815 t_active_nodes.erase(iterator);
819template<
class NodeInfoT>
821 m_relaxations.front()->write(t_name);
824template<
class NodeInfoT>
827 auto children_info = m_branching_rule->create_child_nodes(t_node);
829 std::vector<Node<NodeInfoT>> result;
830 result.reserve(children_info.size());
832 for (
auto* ptr_to_info : children_info) {
833 result.emplace_back(ptr_to_info, m_n_created_nodes++, t_node);
836 m_branching_rule->on_nodes_have_been_created();
841template<
class NodeInfoT>
843 return get_relative_gap() <= get_tol_mip_relative_gap()
844 || get_absolute_gap() <= get_tol_mip_absolute_gap();
847template<
class NodeInfoT>
849 SetOfActiveNodes &t_sub_active_nodes) {
851 t_actives_nodes.merge(std::move(t_sub_active_nodes));
855template<
class NodeInfoT>
857 for (
auto& relaxation : m_relaxations) {
858 relaxation->update();
862template<
class NodeInfoT>
864 for (
auto& relaxation : m_relaxations) {
865 relaxation->remove(t_ctr);
869template<
class NodeInfoT>
871 for (
auto& relaxation : m_relaxations) {
872 relaxation->remove(t_var);
876template<
class NodeInfoT>
878 for (
auto& relaxation : m_relaxations) {
879 relaxation->add(t_ctr);
883template<
class NodeInfoT>
885 for (
auto& relaxation : m_relaxations) {
886 relaxation->add(t_var);
890template<
class NodeInfoT>
893 Algorithm::set_best_obj(t_value);
896template<
class NodeInfoT>
899 Algorithm::set_best_bound(t_value);
902template<
class NodeInfoT>
905 Algorithm::set_reason(t_reason);
908template<
class NodeInfoT>
911 Algorithm::set_status(t_status);