Loading...
Searching...
No Matches
Optimizers_BranchAndBound.h
1//
2// Created by henri on 13/03/23.
3//
4
5#ifndef IDOL_OPTIMIZERS_BRANCHANDBOUND_H
6#define IDOL_OPTIMIZERS_BRANCHANDBOUND_H
7
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"
20
21#include <memory>
22#include <cassert>
23#include <fstream>
24
25namespace idol::Optimizers {
26 template<class NodeInfoT> class BranchAndBound;
27}
28
29template<class NodeInfoT>
30class idol::Optimizers::BranchAndBound : public Algorithm {
31 using TreeNode = Node<NodeInfoT>;
32 using SetOfActiveNodes = NodeSet<Node<NodeInfoT>>;
33
34 std::unique_ptr<OptimizerFactory> m_relaxation_optimizer_factory;
35
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;
39
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;
43 std::unique_ptr<NodeInfoT> m_root_node_info;
44
45 std::unique_ptr<AbstractBranchAndBoundCallbackI<NodeInfoT>> m_callback;
46
47 bool m_has_integer_objective = false;
48 std::vector<unsigned int> m_steps = { std::numeric_limits<unsigned int>::max(), 0, 0 };
49 unsigned int m_n_created_nodes = 0;
50 unsigned int m_n_solved_nodes = 0;
51 unsigned int m_n_active_nodes = 0;
52 double m_root_node_best_bound = -Inf;
53 double m_root_node_best_obj = +Inf;
54
55 std::optional<TreeNode> m_incumbent;
56protected:
57 void build() override;
58 void hook_before_optimize() override;
59 void hook_optimize() override;
60 void hook_after_optimize() override;
61 void add(const Var &t_var) override;
62 void add(const Ctr &t_ctr) override;
63 void add(const QCtr &t_ctr) override;
64 void remove(const Var &t_var) override;
65 void remove(const Ctr &t_ctr) override;
66 void remove(const QCtr &t_ctr) override;
67 void update() override;
68 void write(const std::string &t_name) override;
69
70 void detect_integer_objective();
71 void create_relaxations();
72 Node<NodeInfoT> create_root_node();
73 void explore(TreeNode& t_node, unsigned int t_relaxation_id, SetOfActiveNodes& t_active_nodes, unsigned int t_step);
74 void analyze(const TreeNode& t_node, unsigned int t_relaxation_id, bool* t_explore_children_flag, bool* t_reoptimize_flag);
75 Node<NodeInfoT> select_node_for_branching(SetOfActiveNodes& t_active_nodes);
76 std::vector<TreeNode> create_child_nodes(const TreeNode& t_node);
77 void backtrack(SetOfActiveNodes& t_actives_nodes, SetOfActiveNodes& t_sub_active_nodes);
78 void set_as_incumbent(const TreeNode& t_node);
79 [[nodiscard]] bool gap_is_closed() const;
80 void prune_nodes_by_bound(SetOfActiveNodes& t_active_nodes);
81 void update_lower_bound(const SetOfActiveNodes& t_active_nodes);
82 bool is_valid(const TreeNode& t_node) const;
83
84 void log_node_after_solve(const Node<NodeInfoT>& t_node);
85 void log_after_termination();
86
87 SideEffectRegistry call_callbacks(CallbackEvent t_event, const TreeNode& t_node, unsigned int t_relaxation_id);
88 void update_obj_sense() override;
89 void update_obj() override;
90 void update_rhs() override;
91 void update_obj_constant() override;
92 void update_mat_coeff(const Ctr &t_ctr, const Var &t_var) override;
93 void update_ctr_type(const Ctr &t_ctr) override;
94 void update_ctr_rhs(const Ctr &t_ctr) override;
95 void update_var_type(const Var &t_var) override;
96 void update_var_lb(const Var &t_var) override;
97 void update_var_ub(const Var &t_var) override;
98
99 void update_var_obj(const Var &t_var) override;
100 void set_status(SolutionStatus t_status) override;
101 void set_reason(SolutionReason t_reason) override;
102 void set_best_bound(double t_value) override;
103 void set_best_obj(double t_value) override;
104public:
105 explicit BranchAndBound(const Model& t_model,
106 const OptimizerFactory& t_node_optimizer,
107 const BranchingRuleFactory<NodeInfoT>& t_branching_rule_factory,
108 const NodeSelectionRuleFactory<NodeInfoT>& t_node_selection_rule_factory,
109 AbstractBranchAndBoundCallbackI<NodeInfoT>* t_callback,
110 const Logs::BranchAndBound::Factory<NodeInfoT>& t_logger_factory);
111
112 [[nodiscard]] std::string name() const override { return "Branch-and-Bound"; }
113
114 void solve(TreeNode& t_node, unsigned int t_relaxation_id) const;
115
116 virtual void set_subtree_depth(unsigned int t_depth) { m_steps.at(1) = t_depth; }
117
118 [[nodiscard]] unsigned int subtree_depth() const { return m_steps.at(1); }
119
120 virtual void add_callback(BranchAndBoundCallback<NodeInfoT>* t_cb);
121
122 void submit_heuristic_solution(NodeInfoT* t_info);
123
124 void submit_lower_bound(double t_lower_bound);
125
126 [[nodiscard]] unsigned int n_created_nodes() const { return m_n_created_nodes; }
127
128 [[nodiscard]] unsigned int n_solved_nodes() const { return m_n_solved_nodes; }
129
130 [[nodiscard]] unsigned int n_active_nodes() const { return m_n_active_nodes; }
131
132 [[nodiscard]] const Model& relaxation() const { return *m_relaxations.front(); }
133
134 [[nodiscard]] Model& relaxation() { return *m_relaxations.front(); }
135
136 [[nodiscard]] double root_node_best_bound() const { return m_root_node_best_bound; }
137
138 [[nodiscard]] double root_node_best_obj() const { return m_root_node_best_obj; }
139
140 BranchingRule<NodeInfoT>& branching_rule() { return *m_branching_rule; }
141
142 const BranchingRule<NodeInfoT>& branching_rule() const { return *m_branching_rule; }
143
144 [[nodiscard]] bool has_incumbent() const { return m_incumbent.has_value(); }
145
146 const TreeNode& incumbent() const { return m_incumbent.value(); }
147
148 [[nodiscard]] double get_var_primal(const Var &t_var) const override;
149
150 [[nodiscard]] double get_var_reduced_cost(const Var &t_var) const override;
151
152 [[nodiscard]] double get_var_ray(const Var &t_var) const override;
153
154 [[nodiscard]] double get_ctr_dual(const Ctr &t_ctr) const override;
155
156 [[nodiscard]] double get_ctr_farkas(const Ctr &t_ctr) const override;
157
158 [[nodiscard]] unsigned int get_n_solutions() const override;
159
160 [[nodiscard]] unsigned int get_solution_index() const override;
161
162 void set_solution_index(unsigned int t_index) override;
163
164 void terminate() override;
165
166 void set_relaxation_optimizer_factory(const OptimizerFactory& t_factory);
167
168 void set_root_node_info(const NodeInfoT& t_node_info);
169
170 const OptimizerFactory& get_relaxation_optimizer_factory() const;
171};
172
173template<class NodeInfoT>
174void idol::Optimizers::BranchAndBound<NodeInfoT>::remove(const idol::QCtr &t_ctr) {
175 for (auto& relaxation : m_relaxations) {
176 relaxation->remove(t_ctr);
177 }
178}
179
180template<class NodeInfoT>
181void idol::Optimizers::BranchAndBound<NodeInfoT>::add(const idol::QCtr &t_ctr) {
182 for (auto& relaxation : m_relaxations) {
183 relaxation->add(t_ctr);
184 }
185}
186
187template<class NodeInfoT>
188bool idol::Optimizers::BranchAndBound<NodeInfoT>::is_valid(const TreeNode &t_node) const {
189 bool result;
190#pragma omp critical
191 result = m_branching_rule->is_valid(t_node);
192 return result;
193}
194
195template<class NodeInfoT>
196void idol::Optimizers::BranchAndBound<NodeInfoT>::terminate() {
197#pragma omp critical
198 Optimizer::terminate();
199}
200
201template <class NodeInfoT> void idol::Optimizers::BranchAndBound<NodeInfoT>::set_relaxation_optimizer_factory(
202 const OptimizerFactory& t_factory) {
203
204 m_relaxation_optimizer_factory.reset(t_factory.clone());
205
206}
207
208template <class NodeInfoT>
209void idol::Optimizers::BranchAndBound<NodeInfoT>::set_root_node_info(const NodeInfoT& t_node_info) {
210 m_root_node_info.reset(t_node_info.clone());
211}
212
213template <class NodeInfoT> const idol::OptimizerFactory& idol::Optimizers::BranchAndBound<NodeInfoT>::
214get_relaxation_optimizer_factory() const {
215 return *m_relaxation_optimizer_factory;
216}
217
218template<class NodeInfoT>
219void idol::Optimizers::BranchAndBound<NodeInfoT>::detect_integer_objective() {
220 const auto& src_model = parent();
221 const auto& objective = src_model.get_obj_expr();
222
223 if (!is_integer(objective.affine().constant(), Tolerance::Integer)) {
224 m_has_integer_objective = false;
225 return;
226 }
227
228 for (const auto& [var, val] : objective.affine().linear()) {
229 if (src_model.get_var_type(var) == Continuous || !is_integer(val, Tolerance::Integer)) {
230 m_has_integer_objective = false;
231 return;
232 }
233 }
234
235 m_has_integer_objective = true;
236}
237
238template<class NodeInfoT>
239double idol::Optimizers::BranchAndBound<NodeInfoT>::get_var_reduced_cost(const idol::Var &t_var) const {
240 if (m_n_solved_nodes > 1) {
241 throw Exception("Reduced cost not available.");
242 }
243 return m_relaxations.front()->get_var_reduced_cost(t_var);
244}
245
246template<class NodeInfoT>
247void idol::Optimizers::BranchAndBound<NodeInfoT>::log_after_termination() {
248
249 if (!get_param_logs()) {
250 return;
251 }
252
253 if (m_callback) {
254 m_callback->log_after_termination();
255 }
256
257 m_logger->log_after_termination();
258
259}
260
261template<class NodeInfoT>
262void idol::Optimizers::BranchAndBound<NodeInfoT>::log_node_after_solve(const idol::Node<NodeInfoT> &t_node) {
263 m_logger->log_node_after_solve(t_node);
264}
265
266template<class NodeInfoT>
267void idol::Optimizers::BranchAndBound<NodeInfoT>::set_solution_index(unsigned int t_index) {
268 if (t_index != 0) {
269 throw Exception("Solution index out of bounds.");
270 }
271}
272
273template<class NodeInfoT>
274unsigned int idol::Optimizers::BranchAndBound<NodeInfoT>::get_solution_index() const {
275 return 0;
276}
277
278template<class NodeInfoT>
279unsigned int idol::Optimizers::BranchAndBound<NodeInfoT>::get_n_solutions() const {
280 return m_incumbent.has_value();
281}
282
283template<class NodeInfoT>
284void idol::Optimizers::BranchAndBound<NodeInfoT>::update_obj() {
285 for (auto& relaxation : m_relaxations) {
286 relaxation->set_obj_expr(parent().get_obj_expr());
287 }
288}
289
290template<class NodeInfoT>
291void idol::Optimizers::BranchAndBound<NodeInfoT>::update_rhs() {
292 for (auto& relaxation : m_relaxations) {
293 relaxation->set_rhs_expr(parent().get_rhs_expr());
294 }
295}
296
297template<class NodeInfoT>
298void idol::Optimizers::BranchAndBound<NodeInfoT>::update_obj_constant() {
299 for (auto& relaxation : m_relaxations) {
300 relaxation->set_obj_const(parent().get_obj_expr().affine().constant());
301 }
302}
303
304template<class NodeInfoT>
305void idol::Optimizers::BranchAndBound<NodeInfoT>::update_mat_coeff(const Ctr &t_ctr, const Var &t_var) {
306 for (auto& relaxation : m_relaxations) {
307 relaxation->set_mat_coeff(t_ctr, t_var, parent().get_mat_coeff(t_ctr, t_var));
308 }
309}
310
311template<class NodeInfoT>
312void idol::Optimizers::BranchAndBound<NodeInfoT>::update_ctr_type(const Ctr &t_ctr) {
313 for (auto& relaxation : m_relaxations) {
314 relaxation->set_ctr_type(t_ctr, parent().get_ctr_type(t_ctr));
315 }
316}
317
318template<class NodeInfoT>
319void idol::Optimizers::BranchAndBound<NodeInfoT>::update_ctr_rhs(const Ctr &t_ctr) {
320 for (auto& relaxation : m_relaxations) {
321 relaxation->set_ctr_rhs(t_ctr, parent().get_ctr_rhs(t_ctr));
322 }
323}
324
325template<class NodeInfoT>
326void idol::Optimizers::BranchAndBound<NodeInfoT>::update_var_type(const Var &t_var) {
327 for (auto& relaxation : m_relaxations) {
328 relaxation->set_var_type(t_var, parent().get_var_type(t_var));
329 }
330}
331
332template<class NodeInfoT>
333void idol::Optimizers::BranchAndBound<NodeInfoT>::update_var_lb(const Var &t_var) {
334 for (auto& relaxation : m_relaxations) {
335 relaxation->set_var_lb(t_var, parent().get_var_lb(t_var));
336 }
337}
338
339template<class NodeInfoT>
340void idol::Optimizers::BranchAndBound<NodeInfoT>::update_var_ub(const Var &t_var) {
341 for (auto& relaxation : m_relaxations) {
342 relaxation->set_var_ub(t_var, parent().get_var_ub(t_var));
343 }
344}
345
346template<class NodeInfoT>
347void idol::Optimizers::BranchAndBound<NodeInfoT>::update_var_obj(const Var &t_var) {
348 for (auto& relaxation : m_relaxations) {
349 relaxation->set_var_obj(t_var, parent().get_var_obj(t_var));
350 }
351}
352
353template<class NodeInfoT>
354void idol::Optimizers::BranchAndBound<NodeInfoT>::update_obj_sense() {
355 throw Exception("Not implemented");
356}
357
358template<class NodeInfoT>
359double idol::Optimizers::BranchAndBound<NodeInfoT>::get_ctr_farkas(const Ctr &t_ctr) const {
360 if (m_n_solved_nodes > 1) {
361 throw Exception("Accessing Farkas certificate for MIP is not available after the root node.");
362 }
363 return m_relaxations.front()->get_ctr_farkas(t_ctr);
364}
365
366template<class NodeInfoT>
367double idol::Optimizers::BranchAndBound<NodeInfoT>::get_ctr_dual(const Ctr &t_ctr) const {
368 if (m_n_solved_nodes > 1) {
369 throw Exception("Accessing dual values for MIP is not available after the root node.");
370 }
371 return m_relaxations.front()->get_ctr_dual(t_ctr);
372}
373
374template<class NodeInfoT>
375double idol::Optimizers::BranchAndBound<NodeInfoT>::get_var_ray(const Var &t_var) const {
376 if (m_n_solved_nodes > 1) {
377 throw Exception("Ray not available.");
378 }
379 return m_relaxations.front()->get_var_ray(t_var);
380}
381
382template<class NodeInfoT>
383double idol::Optimizers::BranchAndBound<NodeInfoT>::get_var_primal(const Var &t_var) const {
384
385 if (!m_incumbent.has_value()){
386 throw Exception("Trying to access primal values while no incumbent was found.");
387 }
388
389 return m_incumbent->info().primal_solution().get(t_var);
390}
391
392template<class NodeInfoT>
393void idol::Optimizers::BranchAndBound<NodeInfoT>::submit_lower_bound(double t_lower_bound) {
394 if (t_lower_bound > get_best_obj()) {
395 set_status(Fail);
396 set_reason(NotSpecified);
397 terminate();
398 return;
399 }
400 if (t_lower_bound > get_best_bound()) {
401 set_best_bound(t_lower_bound);
402 // New best LB
403 }
404}
405
406template<class NodeInfoT>
407void idol::Optimizers::BranchAndBound<NodeInfoT>::submit_heuristic_solution(NodeInfoT* t_info) {
408
409 auto t_node = Node<NodeInfoT>::create_detached_node(t_info);
410
411 if (t_node.info().objective_value() < get_best_obj()) {
412
413 //if (m_branching_rule->is_valid(t_node)) {
414 set_as_incumbent(t_node);
415 // New incumbent by submission
416 //} else {
417 // idol_Log(Trace, "Ignoring submitted heuristic solution, solution is not valid.");
418 //}
419
420 } else {
421 // idol_Log(Trace, "Ignoring submitted heuristic solution, objective value is " << t_node.info().objective_value() << " while best obj is " << get_best_obj() << '.');
422 }
423
424}
425
426template<class NodeInfoT>
427idol::SideEffectRegistry
428idol::Optimizers::BranchAndBound<NodeInfoT>::call_callbacks(CallbackEvent t_event, const Optimizers::BranchAndBound<NodeInfoT>::TreeNode &t_node, unsigned int t_relaxation_id) {
429
430 if (!m_callback) {
431 return {};
432 }
433
434 return m_callback->operator()(this, t_event, t_node, m_relaxations[t_relaxation_id].get());
435}
436
437template<class NodeInfoT>
438void idol::Optimizers::BranchAndBound<NodeInfoT>::add_callback(BranchAndBoundCallback<NodeInfoT> *t_cb) {
439 m_callback->add_callback(t_cb);
440}
441
442template<class NodeInfoT>
443void idol::Optimizers::BranchAndBound<NodeInfoT>::build() {
444
445}
446
447template<class NodeInfoT>
448idol::Optimizers::BranchAndBound<NodeInfoT>::BranchAndBound(const Model &t_model,
449 const OptimizerFactory& t_node_optimizer,
450 const BranchingRuleFactory<NodeInfoT>& t_branching_rule_factory,
451 const NodeSelectionRuleFactory<NodeInfoT>& t_node_selection_rule_factory,
452 AbstractBranchAndBoundCallbackI<NodeInfoT>* t_callback,
453 const Logs::BranchAndBound::Factory<NodeInfoT>& t_logger_factory)
454 : Algorithm(t_model),
455 m_relaxation_optimizer_factory(t_node_optimizer.clone()),
456 m_branching_rule(t_branching_rule_factory(*this)),
457 m_node_selection_rule(t_node_selection_rule_factory(*this)),
458 m_callback(t_callback),
459 m_logger(t_logger_factory.operator()(*this)) {
460
461 create_relaxations();
462
463}
464
465template<class NodeInfoT>
466void idol::Optimizers::BranchAndBound<NodeInfoT>::create_relaxations() {
467
468 const auto &original_model = parent();
469
470 m_relaxations.clear();
471 m_relaxations.reserve(m_n_threads);
472
473 m_node_updators.clear();
474 m_node_updators.reserve(m_n_threads);
475
476 for (unsigned int i = 0 ; i < m_n_threads ; ++i) {
477 auto* relaxation = original_model.clone();
478 relaxation->use(*m_relaxation_optimizer_factory);
479 m_relaxations.emplace_back(relaxation);
480 m_node_updators.emplace_back(dynamic_cast<NodeUpdator<NodeInfoT>*>(NodeInfoT::create_updator(parent(), *relaxation)));
481 }
482
483}
484
485template<class NodeInfoT>
486idol::Node<NodeInfoT> idol::Optimizers::BranchAndBound<NodeInfoT>::create_root_node() {
487
488 NodeInfoT* node_info = nullptr;
489 if (m_root_node_info) {
490 node_info = m_root_node_info->clone();
491 } else if constexpr (std::is_default_constructible_v<NodeInfoT>) {
492 node_info = new NodeInfoT();
493 } else {
494 throw Exception("The given node info class has no default constructor and no root node info has been given.");
495 }
496 auto root_node = Node<NodeInfoT>::create_root_node(node_info);
497 assert(root_node.id() == 0);
498 ++m_n_created_nodes;
499
500 return root_node;
501
502}
503
504template<class NodeInfoT>
505void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_before_optimize() {
506 Algorithm::hook_before_optimize();
507
508 // Reset solution
509 set_best_bound(std::max(-Inf, get_param_best_obj_stop()));
510 set_best_obj(std::min(+Inf, get_param_best_bound_stop()));
511 m_incumbent.reset();
512
513 detect_integer_objective();
514
515 m_n_created_nodes = 0;
516 m_n_solved_nodes = 0;
517 m_n_active_nodes = 0;
518
519 m_root_node_best_bound = -Inf;
520 m_root_node_best_obj = Inf;
521
522 m_branching_rule->initialize();
523 for (auto& node_updator : m_node_updators) {
524 node_updator->initialize();
525 }
526 m_callback->initialize(this);
527
528 m_logger->initialize();
529}
530
531template<class NodeInfoT>
532void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_optimize() {
533
534 auto root_node = create_root_node();
535
536 SetOfActiveNodes active_nodes;
537
538 explore(root_node, 0, active_nodes, 0);
539
540 if (active_nodes.empty()) {
541 set_best_bound(get_best_obj());
542 }
543
544 m_node_updators[0]->clear();
545
546 if (get_status() == Fail) {
547 return;
548 }
549
550 if (!m_incumbent.has_value()) {
551
552 if (is_pos_inf(get_best_obj())) {
553 set_status(Infeasible);
554 return;
555 }
556
557 if (is_neg_inf(get_best_obj())) {
558 set_status(Unbounded);
559 return;
560 }
561
562 }
563
564 set_status(Feasible);
565
566 if (gap_is_closed()) {
567 set_status(Optimal);
568 set_reason(Proved);
569 return;
570 }
571
572}
573
574template<class NodeInfoT>
575void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_after_optimize() {
576
577 Optimizer::hook_after_optimize();
578
579 log_after_termination();
580
581}
582
583template<class NodeInfoT>
584void idol::Optimizers::BranchAndBound<NodeInfoT>::explore(
585 TreeNode &t_node,
586 unsigned int t_relaxation_id,
587 SetOfActiveNodes & t_active_nodes,
588 unsigned int t_step) {
589
590 if (is_terminated() || gap_is_closed()) { return; }
591
592 bool reoptimize_flag;
593 bool explore_children_flag;
594
595 do {
596
597 solve(t_node, t_relaxation_id);
598
599 analyze(t_node, t_relaxation_id, &explore_children_flag, &reoptimize_flag);
600
601 log_node_after_solve(t_node);
602
603 } while (reoptimize_flag);
604
605#pragma omp atomic
606 ++m_n_solved_nodes;
607
608 if (is_terminated() || gap_is_closed()) { return; }
609
610 if (!explore_children_flag) { return; }
611
612 t_active_nodes.emplace(t_node);
613
614 const unsigned int max_levels = m_steps.at(t_step);
615
616 for (unsigned int level = 0 ; level < max_levels ; ++level) {
617
618 prune_nodes_by_bound(t_active_nodes);
619
620 if (t_step == 0) {
621 update_lower_bound(t_active_nodes);
622 m_n_active_nodes = t_active_nodes.size();
623 }
624
625 if (t_active_nodes.empty()) { break; }
626
627 if (is_terminated() || gap_is_closed()) { break; }
628
629 auto selected_node = select_node_for_branching(t_active_nodes);
630
631 auto children = create_child_nodes(selected_node);
632
633 const unsigned int n_children = children.size();
634
635 std::vector<SetOfActiveNodes> active_nodes(n_children);
636
637 if (m_n_threads > 1 && t_step == 0) {
638
639 std::vector<std::vector<unsigned int>> to_explore(m_n_threads);
640
641 for (unsigned int q = 0 ; q < n_children ; ++q) {
642 const unsigned int relaxation_id = q % m_n_threads;
643 to_explore[relaxation_id].emplace_back(q);
644 }
645
646#pragma omp parallel for num_threads(m_n_threads)
647 for (unsigned int k = 0 ; k < m_n_threads ; ++k) {
648 for (unsigned int j = 0, n_tasks = to_explore[k].size() ; j < n_tasks ; ++j) {
649 const unsigned int q = to_explore[k][j];
650 explore(children[q], k, active_nodes[q], t_step + 1);
651 }
652 }
653
654 } else {
655
656 for (unsigned int q = 0 ; q < n_children ; ++q) {
657 explore(children[q], t_relaxation_id, active_nodes[q], t_step + 1);
658 }
659
660 }
661
662 for (unsigned int q = 0 ; q < n_children ; ++q) {
663 backtrack(t_active_nodes, active_nodes[q]);
664 }
665
666 }
667
668}
669
670template<class NodeInfoT>
671void idol::Optimizers::BranchAndBound<NodeInfoT>::solve(TreeNode& t_node,
672 unsigned int t_relaxation_id) const {
673
674 auto& node_updator = *m_node_updators[t_relaxation_id];
675 auto& relaxation = *m_relaxations[t_relaxation_id];
676
677 relaxation.optimizer().set_param_best_bound_stop(std::min(get_best_obj(), get_param_best_bound_stop()));
678 relaxation.optimizer().set_param_time_limit(get_remaining_time());
679
680 node_updator.prepare(t_node);
681
682 // TODO call callback "before node"
683
684 /*
685 for (const auto& var : parent().vars()) {
686 const double lb = relaxation.get_var_lb(var);
687 const double ub = relaxation.get_var_ub(var);
688 if (lb > ub + Tolerance::Integer) {
689 std::cerr << "Inconsistent bounds for variable " << var << ": " << lb << " > " << ub << std::endl;
690 }
691 }
692 */
693
694 relaxation.optimize();
695
696 t_node.info().save(parent(), relaxation);
697
698 m_branching_rule->on_node_solved(t_node);
699
700}
701
702template<class NodeInfoT>
703void idol::Optimizers::BranchAndBound<NodeInfoT>::analyze(const BranchAndBound::TreeNode &t_node, unsigned int t_relaxation_id, bool* t_explore_children_flag, bool* t_reoptimize_flag) {
704
705 *t_explore_children_flag = false;
706 *t_reoptimize_flag = false;
707
708 if (get_best_bound() > get_best_obj()) {
709 set_status(Fail);
710 set_reason(NotSpecified);
711 terminate();
712 return;
713 }
714
715 const auto& status = t_node.info().status();
716 const auto& reason = t_node.info().reason();
717
718 if (get_remaining_time() == 0 || reason == TimeLimit) {
719 set_reason(TimeLimit);
720 terminate();
721 return;
722 }
723
724 if (status == Unbounded) {
725 set_reason(Proved);
726 set_best_obj(-Inf);
727 terminate();
728 return;
729 }
730
731 if (status == Infeasible || status == InfOrUnbnd) {
732
733 if (t_node.level() == 0) {
734 set_status(status);
735 }
736
737 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
738 return;
739 }
740
741 if (status == Feasible && reason == ObjLimit) {
742 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
743 return;
744 }
745
746 if (status == Fail || status == Loaded) {
747
748 set_status(Fail);
749 set_reason(NotSpecified);
750 terminate();
751 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
752 return;
753 }
754
755 if (t_node.level() == 0) {
756 m_root_node_best_bound = t_node.info().objective_value();
757 set_best_bound(std::max(m_root_node_best_bound, get_best_bound()));
758 }
759
760 if (t_node.info().objective_value() < get_best_obj()) {
761
762 if (is_valid(t_node)) {
763
764 auto side_effects = call_callbacks(IncumbentSolution, t_node, t_relaxation_id);
765
766 if (side_effects.n_added_user_cuts > 0 || side_effects.n_added_lazy_cuts > 0) {
767 *t_reoptimize_flag = true;
768 return;
769 }
770
771 set_as_incumbent(t_node);
772 return;
773
774 }
775
776 } else {
777
778 call_callbacks(PrunedSolution, t_node, t_relaxation_id);
779 return;
780
781 }
782
783 auto side_effects = call_callbacks(InvalidSolution, t_node, t_relaxation_id);
784
785 if (t_node.level() == 0) {
786 m_root_node_best_obj = get_best_obj();
787 }
788
789 if (side_effects.n_added_lazy_cuts > 0 || side_effects.n_added_user_cuts > 0 || side_effects.n_added_local_variable_branching > 0) {
790 *t_reoptimize_flag = true;
791 return;
792 }
793
794 *t_explore_children_flag = true;
795
796}
797
798template<class NodeInfoT>
799void idol::Optimizers::BranchAndBound<NodeInfoT>::set_as_incumbent(const BranchAndBound::TreeNode &t_node) {
800#pragma omp critical
801 m_incumbent = t_node;
802 set_best_obj(t_node.info().objective_value());
803 set_status(Feasible);
804}
805
806template<class NodeInfoT>
807void idol::Optimizers::BranchAndBound<NodeInfoT>::update_lower_bound(const BranchAndBound::SetOfActiveNodes &t_active_nodes) {
808
809 if (t_active_nodes.empty()) { return; }
810
811 auto& lowest_node = *t_active_nodes.by_objective_value().begin();
812 const double raw_lower_bound = lowest_node.info().objective_value();
813 const double lower_bound = m_has_integer_objective && !is_integer(raw_lower_bound, Tolerance::Integer) ? std::ceil(raw_lower_bound) : raw_lower_bound;
814 if (lower_bound > get_best_bound()) {
815 set_best_bound(lower_bound);
816 }
817
818}
819
820template<class NodeInfoT>
821void idol::Optimizers::BranchAndBound<NodeInfoT>::prune_nodes_by_bound(BranchAndBound::SetOfActiveNodes& t_active_nodes) {
822
823 const double upper_bound = get_best_obj();
824
825 auto it = t_active_nodes.by_objective_value().begin();
826 auto end = t_active_nodes.by_objective_value().end();
827
828 while (it != end) {
829
830 const auto& node = *it;
831 const double raw_lower_bound = node.info().objective_value();
832 const double lower_bound = m_has_integer_objective && !is_integer(raw_lower_bound, Tolerance::Integer) ? std::ceil(raw_lower_bound) : raw_lower_bound;
833
834 if (lower_bound >= upper_bound - get_tol_mip_absolute_gap()) {
835 it = t_active_nodes.erase(it);
836 end = t_active_nodes.by_objective_value().end();
837 } else {
838 break;
839 }
840
841 }
842
843}
844
845template<class NodeInfoT>
846idol::Node<NodeInfoT>
847idol::Optimizers::BranchAndBound<NodeInfoT>::select_node_for_branching(BranchAndBound::SetOfActiveNodes &t_active_nodes) {
848 auto iterator = m_node_selection_rule->operator()(t_active_nodes);
849 auto result = *iterator;
850 t_active_nodes.erase(iterator);
851 return result;
852}
853
854template<class NodeInfoT>
855void idol::Optimizers::BranchAndBound<NodeInfoT>::write(const std::string &t_name) {
856 m_relaxations.front()->write(t_name);
857}
858
859template<class NodeInfoT>
860std::vector<idol::Node<NodeInfoT>> idol::Optimizers::BranchAndBound<NodeInfoT>::create_child_nodes(const BranchAndBound::TreeNode &t_node) {
861
862 auto children_info = m_branching_rule->create_child_nodes(t_node);
863
864 std::vector<Node<NodeInfoT>> result;
865 result.reserve(children_info.size());
866
867 for (auto* ptr_to_info : children_info) {
868 result.emplace_back(ptr_to_info, m_n_created_nodes++, t_node);
869 }
870
871 m_branching_rule->on_nodes_have_been_created();
872
873 return result;
874}
875
876template<class NodeInfoT>
877bool idol::Optimizers::BranchAndBound<NodeInfoT>::gap_is_closed() const {
878 return get_relative_gap() <= get_tol_mip_relative_gap()
879 || get_absolute_gap() <= get_tol_mip_absolute_gap();
880}
881
882template<class NodeInfoT>
883void idol::Optimizers::BranchAndBound<NodeInfoT>::backtrack(BranchAndBound::SetOfActiveNodes &t_actives_nodes,
884 SetOfActiveNodes &t_sub_active_nodes) {
885
886 t_actives_nodes.merge(std::move(t_sub_active_nodes));
887
888}
889
890template<class NodeInfoT>
891void idol::Optimizers::BranchAndBound<NodeInfoT>::update() {
892 for (auto& relaxation : m_relaxations) {
893 relaxation->update();
894 }
895}
896
897template<class NodeInfoT>
898void idol::Optimizers::BranchAndBound<NodeInfoT>::remove(const Ctr &t_ctr) {
899 for (auto& relaxation : m_relaxations) {
900 relaxation->remove(t_ctr);
901 }
902}
903
904template<class NodeInfoT>
905void idol::Optimizers::BranchAndBound<NodeInfoT>::remove(const Var &t_var) {
906 for (auto& relaxation : m_relaxations) {
907 relaxation->remove(t_var);
908 }
909}
910
911template<class NodeInfoT>
912void idol::Optimizers::BranchAndBound<NodeInfoT>::add(const Ctr &t_ctr) {
913 for (auto& relaxation : m_relaxations) {
914 relaxation->add(t_ctr);
915 }
916}
917
918template<class NodeInfoT>
919void idol::Optimizers::BranchAndBound<NodeInfoT>::add(const Var &t_var) {
920 for (auto& relaxation : m_relaxations) {
921 relaxation->add(t_var);
922 }
923}
924
925template<class NodeInfoT>
926void idol::Optimizers::BranchAndBound<NodeInfoT>::set_best_obj(double t_value) {
927#pragma omp critical
928 Algorithm::set_best_obj(t_value);
929 assert(get_best_bound() <= get_best_obj() + get_tol_mip_absolute_gap());
930}
931
932template<class NodeInfoT>
933void idol::Optimizers::BranchAndBound<NodeInfoT>::set_best_bound(double t_value) {
934#pragma omp critical
935 Algorithm::set_best_bound(t_value);
936 assert(get_best_bound() <= get_best_obj() + get_tol_mip_absolute_gap());
937}
938
939template<class NodeInfoT>
940void idol::Optimizers::BranchAndBound<NodeInfoT>::set_reason(idol::SolutionReason t_reason) {
941#pragma omp critical
942 Algorithm::set_reason(t_reason);
943}
944
945template<class NodeInfoT>
946void idol::Optimizers::BranchAndBound<NodeInfoT>::set_status(idol::SolutionStatus t_status) {
947#pragma omp critical
948 Algorithm::set_status(t_status);
949}
950
951#endif //IDOL_OPTIMIZERS_BRANCHANDBOUND_H
idol::Tolerance::Integer
static double Integer
Definition numericals.h:73
Generated by doxygen 1.16.1