Replace the tree update with a lock-free version

cpp/src/branch_and_bound/branch_and_bound.hpp

@@ -157,7 +157,7 @@ class branch_and_bound_t {
   void set_initial_pseudocost(const pseudo_costs_t<i_t, f_t>& parent_pc,
                               const std::vector<i_t>& reduced_to_original);
 
-  f_t get_upper_bound() const { return upper_bound_.load(); }
+  f_t get_upper_bound() const { return upper_bound_.load(std::memory_order::relaxed); }
   bool has_solver_space_incumbent() const { return incumbent_.has_incumbent; }
 
   // Repair a low-quality solution from the heuristics.

cpp/src/branch_and_bound/mip_node.hpp

@@ -80,8 +80,9 @@ class mip_node_t {
       objective_estimate(std::numeric_limits<f_t>::infinity()),
       packed_vstatus(0)
   {
-    children[0] = nullptr;
-    children[1] = nullptr;
+    children[0]  = nullptr;
+    children[1]  = nullptr;
+    num_children = 0;
   }
 
   mip_node_t(f_t root_lower_bound, const std::vector<simplex::variable_status_t>& basis)
@@ -96,8 +97,9 @@ class mip_node_t {
       objective_estimate(std::numeric_limits<f_t>::infinity()),
       packed_vstatus(compress_vstatus(basis))
   {
-    children[0] = nullptr;
-    children[1] = nullptr;
+    children[0]  = nullptr;
+    children[1]  = nullptr;
+    num_children = 0;
   }
 
   mip_node_t(const simplex::lp_problem_t<i_t, f_t>& problem,
@@ -126,6 +128,7 @@ class mip_node_t {
                                                                     : problem.upper[branch_var];
     children[0]      = nullptr;
     children[1]      = nullptr;
+    num_children     = 0;
   }
 
   void get_variable_bounds(std::vector<f_t>& lower,
@@ -177,54 +180,24 @@ class mip_node_t {
     children[1] = std::move(up_child);
     // When we add children we no longer need to store our basis
     packed_vstatus = {};
+    num_children   = 2;
   }
 
-  bool is_inactive() const
-  {
-    if (inactive_status(status)) { return true; }
-    if ((children[0] != nullptr && inactive_status(children[0]->status)) &&
-        (children[1] != nullptr && inactive_status(children[1]->status))) {
-      return true;
-    }
-    if (children[0] == nullptr && inactive_status(children[1]->status)) { return true; }
-    if (children[1] == nullptr && inactive_status(children[0]->status)) { return true; }
-    return false;
-  }
-
-  void update_bound()
-  {
-    if (children[0] != nullptr && children[1] != nullptr) {
-      if (inactive_status(children[0]->status) && inactive_status(children[1]->status)) {
-        lower_bound = std::min(children[0]->lower_bound, children[1]->lower_bound);
-      }
-    }
-    if (children[0] != nullptr && children[1] == nullptr) {
-      if (inactive_status(children[0]->status)) { lower_bound = children[0]->lower_bound; }
-    }
-    if (children[1] != nullptr && children[0] == nullptr) {
-      if (inactive_status(children[1]->status)) { lower_bound = children[1]->lower_bound; }
-    }
-  }
-
-  // outputs a stack containing inactive nodes in the tree that can be freed
-  void set_status(node_status_t node_status, std::vector<mip_node_t*>& stack)
+  void set_status(node_status_t node_status)
   {
     status = node_status;
-    if (inactive_status(status)) {
-      update_bound();
-      stack.push_back(this);
-      // Propagate to parent
-      mip_node_t* parent_ptr = parent;
-      while (parent_ptr != nullptr) {
-        if (parent_ptr->is_inactive()) {
-          parent_ptr->status = node_status_t::FATHOMED;
-          parent_ptr->update_bound();
-          stack.push_back(parent_ptr);
-        } else {
-          break;
-        }
-        parent_ptr = parent_ptr->parent;
-      }
+    if (!inactive_status(node_status)) { return; }
+
+    mip_node_t* parent_ptr = parent;
+    while (parent_ptr != nullptr) {
+      i_t old = parent_ptr->num_children.fetch_sub(1, std::memory_order::acq_rel);
+      if (old != 1) { return; }  // The sibling is still active.
+      parent_ptr->status = node_status_t::FATHOMED;
+      parent_ptr->lower_bound =
+        std::min(parent_ptr->children[0]->lower_bound, parent_ptr->children[1]->lower_bound);
+      parent_ptr->children[0].reset();
+      parent_ptr->children[1].reset();
+      parent_ptr = parent_ptr->parent;
     }
   }
 
@@ -302,8 +275,9 @@ class mip_node_t {
   f_t fractional_val;
   i_t integer_infeasible;
 
-  mip_node_t<i_t, f_t>* parent;
+  mip_node_t* parent;
   std::unique_ptr<mip_node_t> children[2];
+  omp_atomic_t<i_t> num_children;
 
   std::vector<uint8_t> packed_vstatus;
 
@@ -344,31 +318,13 @@ class mip_node_t {
   }
 };
 
-template <typename i_t, typename f_t>
-void remove_fathomed_nodes(std::vector<mip_node_t<i_t, f_t>*>& stack)
-{
-  for (int i = 0; i < stack.size(); ++i) {
-    for (int child = 0; child < 2; ++child) {
-      if (stack[i]->children[child] != nullptr) { stack[i]->children[child].reset(); }
-    }
-  }
-}
-
 template <typename i_t, typename f_t>
 class search_tree_t {
  public:
   search_tree_t() : num_nodes(0) {}
 
   search_tree_t(mip_node_t<i_t, f_t>&& node) : root(std::move(node)), num_nodes(0) {}
 
-  void update(mip_node_t<i_t, f_t>* node_ptr, node_status_t status)
-  {
-    std::lock_guard<omp_mutex_t> lock(mutex);
-    std::vector<mip_node_t<i_t, f_t>*> stack;
-    node_ptr->set_status(status, stack);
-    remove_fathomed_nodes(stack);
-  }
-
   void branch(mip_node_t<i_t, f_t>* parent_node,
               const i_t branch_var,
               const f_t fractional_val,
@@ -377,7 +333,7 @@ class search_tree_t {
               const simplex::lp_problem_t<i_t, f_t>& original_lp,
               simplex::logger_t& log)
   {
-    i_t id = num_nodes.fetch_add(2);
+    i_t id = num_nodes.fetch_add(2, std::memory_order::relaxed);
 
     auto down_child = std::make_unique<mip_node_t<i_t, f_t>>(original_lp,
                                                              parent_node,
@@ -443,7 +399,6 @@ class search_tree_t {
   }
 
   mip_node_t<i_t, f_t> root;
-  omp_mutex_t mutex;
   omp_atomic_t<i_t> num_nodes;
 
   static constexpr bool write_graphviz = false;

cpp/src/branch_and_bound/node_queue.hpp

@@ -44,7 +44,7 @@ class heap_t {
     buffer.emplace_back(std::forward<Args>(args)...);
     std::push_heap(buffer.begin(), buffer.end(), comp);
     ++num_entries_;
-    assert(num_entries_.load() == buffer.size());
+    assert(num_entries_.load(std::memory_order::relaxed) == buffer.size());
   }
 
   T pop()
@@ -53,7 +53,7 @@ class heap_t {
     std::pop_heap(buffer.begin(), buffer.end(), comp);
     T node = std::move(buffer.back());
     buffer.pop_back();
-    assert(num_entries_.load() == buffer.size());
+    assert(num_entries_.load(std::memory_order::relaxed) == buffer.size());
     return node;
   }
 
@@ -154,7 +154,8 @@ class node_queue_t {
 
   f_t get_lower_bound()
   {
-    return best_first_heap_.empty() ? std::numeric_limits<f_t>::infinity() : lower_bound_.load();
+    return best_first_heap_.empty() ? std::numeric_limits<f_t>::infinity()
+                                    : lower_bound_.load(std::memory_order::relaxed);
   }
 
  private:

cpp/src/branch_and_bound/worker.hpp

@@ -301,11 +301,11 @@ class diving_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
   // Set this node inactive
   void set_inactive()
   {
-    if (!this->is_active.load()) { return; }
+    if (!this->is_active.load(std::memory_order::relaxed)) { return; }
     this->is_active = false;
 
     if (bfs_worker) {
-      assert(bfs_worker->active_diving_workers.load() > 0);
+      assert(bfs_worker->active_diving_workers.load(std::memory_order::relaxed) > 0);
       --bfs_worker->active_diving_workers;
     }
   }

cpp/src/branch_and_bound/worker_pool.hpp

@@ -53,7 +53,7 @@ class worker_pool_t {
       workers_[i]->symmetry_ptr = symmetry;
     }
 
-    is_initialized_ = true;
+    is_initialized_.store(true, std::memory_order::release);
   }
 
   WorkerType* pop_idle_worker()
@@ -65,7 +65,8 @@ class worker_pool_t {
       i_t idx = idle_workers_.front();
       idle_workers_.pop_front();
       num_idle_workers_--;
-      assert(idle_workers_.size() == static_cast<size_t>(num_idle_workers_.load()));
+      assert(idle_workers_.size() ==
+             static_cast<size_t>(num_idle_workers_.load(std::memory_order::relaxed)));
       assert(idx >= 0 && static_cast<size_t>(idx) < workers_.size());
       return workers_[idx].get();
     }
@@ -75,13 +76,14 @@ class worker_pool_t {
   {
     assert(worker != nullptr);
     worker->set_inactive();
-    assert(!worker->is_active.load());
+    assert(!worker->is_active.load(std::memory_order::relaxed));
 
-    if (!is_initialized_) return;
+    if (!is_initialized()) return;
 
     std::lock_guard lock(mutex_);
     assert(workers_[worker->worker_id].get() == worker);
-    assert(static_cast<size_t>(num_idle_workers_.load()) == idle_workers_.size());
+    assert(static_cast<size_t>(num_idle_workers_.load(std::memory_order::relaxed)) ==
+           idle_workers_.size());
     assert(idle_workers_.size() <= workers_.size());
 
     idle_workers_.push_back(worker->worker_id);
@@ -101,14 +103,14 @@ class worker_pool_t {
     return workers_[id].get();
   }
 
-  bool is_initialized() const { return is_initialized_; }
+  bool is_initialized() const { return is_initialized_.load(std::memory_order::acquire); }
 
   i_t num_idle() const { return num_idle_workers_; }
   i_t size() const { return workers_.size(); }
 
  private:
   std::vector<std::unique_ptr<WorkerType>> workers_;
-  bool is_initialized_ = false;
+  omp_atomic_t<bool> is_initialized_{false};
 
   omp_mutex_t mutex_;
   circular_deque_t<i_t> idle_workers_;

cpp/src/mip_heuristics/local_search/local_search.cu

@@ -111,13 +111,12 @@ void local_search_t<i_t, f_t>::start_cpufj_scratch_threads(population_t<i_t, f_t
         f_t obj, const std::vector<f_t>& h_vec, double /*work_units*/) {
         population.add_external_solution(h_vec, obj, solution_origin_t::CPUFJ);
         (void)problem_ptr;
-        if (obj < this->local_search_best_obj) {
+        if (obj < this->local_search_best_obj.fetch_min(obj, std::memory_order::relaxed)) {
           CUOPT_LOG_TRACE("******* New local search best obj %g, best overall %g",
                           problem_ptr->get_user_obj_from_solver_obj(obj),
                           problem_ptr->get_user_obj_from_solver_obj(
                             population.is_feasible() ? population.best_feasible().get_objective()
                                                      : std::numeric_limits<f_t>::max()));
-          this->local_search_best_obj = obj;
         }
       };
     counter++;
@@ -229,13 +228,12 @@ void local_search_t<i_t, f_t>::start_cpufj_lptopt_scratch_threads(
   scratch_cpu_fj_on_lp_opt->improvement_callback =
     [this, &population](f_t obj, const std::vector<f_t>& h_vec, double /*work_units*/) {
       population.add_external_solution(h_vec, obj, solution_origin_t::CPUFJ);
-      if (obj < this->local_search_best_obj) {
+      if (obj < this->local_search_best_obj.fetch_min(obj, std::memory_order::relaxed)) {
         CUOPT_LOG_DEBUG("******* New local search best obj %g, best overall %g",
                         context.problem_ptr->get_user_obj_from_solver_obj(obj),
                         context.problem_ptr->get_user_obj_from_solver_obj(
                           population.is_feasible() ? population.best_feasible().get_objective()
                                                    : std::numeric_limits<f_t>::max()));
-        this->local_search_best_obj = obj;
       }
     };
 

cpp/src/mip_heuristics/root_heuristics.hpp

@@ -305,8 +305,8 @@ struct root_heuristics_t {
     // On the first pass we use a thread to generate the clique table
     const i_t cut_generation     = cut_pass == 0 ? 2 : 1;
     const i_t persistent_workers = persistent_lanes_.size();
-    return std::max(i_t{0},
-                    max_workers_ - worker_count_->load() - persistent_workers - cut_generation);
+    const i_t worker_count       = worker_count_->load(std::memory_order::relaxed);
+    return std::max(i_t{0}, max_workers_ - worker_count - persistent_workers - cut_generation);
   }
 
   std::shared_ptr<cut_pass_heuristics_t<i_t, f_t>> create_new_cut_pass_heuristic(