fix: give each solver its own seed instead of a process-wide counter

cpp/include/cuopt/routing/solver_settings.hpp

@@ -64,12 +64,27 @@ class solver_settings_t {
    */
   void dump_best_results(const std::string& file_path, i_t interval);
 
+  /**
+   * @brief Set the random seed used by the routing solver.
+   *
+   * Controls the initial seed for random number generation. Use -1 to derive the seed
+   * from the problem, which is the default and reproduces a given problem run to run.
+   *
+   * @param[in] seed The seed, or -1 to derive it from the problem
+   */
+  void set_seed(i_t seed);
+
   /**
    * @brief Return set solving time
    * @return Solving time set in seconds
    */
   f_t get_time_limit() const noexcept;
 
+  /**
+   * @brief Return the random seed, or -1 if it is derived from the problem
+   */
+  i_t get_seed() const noexcept;
+
   /**
    * @brief Return true if verbose mode is enabled
    */
@@ -93,6 +108,7 @@ class solver_settings_t {
   i_t dump_interval_{std::numeric_limits<i_t>::max()};
   bool dump_best_results_{false};
   std::string best_result_file_name_;
+  i_t seed_{-1};
 };
 
 }  // namespace CUOPT_EXPORT routing

cpp/src/CMakeLists.txt

@@ -3,8 +3,7 @@
 # SPDX-License-Identifier: Apache-2.0
 # cmake-format: on
 
-set(UTIL_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/utilities/seed_generator.cu
-                   ${CMAKE_CURRENT_SOURCE_DIR}/utilities/logger.cpp
+set(UTIL_SRC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/utilities/logger.cpp
                    ${CMAKE_CURRENT_SOURCE_DIR}/utilities/version_info.cpp
                    ${CMAKE_CURRENT_SOURCE_DIR}/utilities/timestamp_utils.cpp
                    ${CMAKE_CURRENT_SOURCE_DIR}/utilities/work_unit_scheduler.cpp)

cpp/src/mip_heuristics/diversity/diversity_manager.cu

@@ -77,10 +77,10 @@ diversity_manager_t<i_t, f_t>::diversity_manager_t(mip_solver_context_t<i_t, f_t
                             context.problem_ptr->handle_ptr),
     sub_mip_recombiner(
       context, population, context.problem_ptr->n_variables, context.problem_ptr->handle_ptr),
-    rng(cuopt::seed_generator::get_seed()),
+    rng(context.problem_ptr->seed_gen.get_seed()),
     stats(context.stats),
-    mab_recombiner(0, cuopt::seed_generator::get_seed(), recombiner_alpha, "recombiner"),
-    mab_ls(mab_ls_config_t<i_t, f_t>::n_of_arms, cuopt::seed_generator::get_seed(), ls_alpha, "ls"),
+    mab_recombiner(0, context.problem_ptr->seed_gen.get_seed(), recombiner_alpha, "recombiner"),
+    mab_ls(mab_ls_config_t<i_t, f_t>::n_of_arms, context.problem_ptr->seed_gen.get_seed(), ls_alpha, "ls"),
     ls_hash_map(*context.problem_ptr)
 {
   int max_config             = -1;

cpp/src/mip_heuristics/diversity/population.cu

@@ -42,7 +42,7 @@ population_t<i_t, f_t>::population_t(std::string const& name_,
     max_solutions(max_solutions_),
     infeasibility_importance(infeasibility_weight_),
     weights(0, context.problem_ptr->handle_ptr),
-    rng(cuopt::seed_generator::get_seed()),
+    rng(context.problem_ptr->seed_gen.get_seed()),
     early_exit_primal_generation(false),
     population_hash_map(*problem_ptr),
     timer(0)

cpp/src/mip_heuristics/diversity/recombiners/bound_prop_recombiner.cuh

@@ -26,7 +26,7 @@ class bound_prop_recombiner_t : public recombiner_t<i_t, f_t> {
                           const raft::handle_t* handle_ptr)
     : recombiner_t<i_t, f_t>(context, n_vars, handle_ptr),
       constraint_prop(constraint_prop_),
-      rng(cuopt::seed_generator::get_seed()),
+      rng(context.problem_ptr->seed_gen.get_seed()),
       vars_to_fix(n_vars, handle_ptr->get_stream())
   {
   }
@@ -65,7 +65,7 @@ class bound_prop_recombiner_t : public recombiner_t<i_t, f_t> {
        offspring_view,
        int_tol,
        probing_values = probing_values.data(),
-       seed           = cuopt::seed_generator::get_seed()] __device__(i_t idx) {
+       seed           = this->context.problem_ptr->seed_gen.get_seed()] __device__(i_t idx) {
         f_t guiding_val = guiding_view.assignment[idx];
         f_t other_val   = other_view.assignment[idx];
         cuopt_assert(guiding_view.problem.check_variable_within_bounds(idx, guiding_val), "");
@@ -151,7 +151,7 @@ class bound_prop_recombiner_t : public recombiner_t<i_t, f_t> {
     if (n_different_vars > (i_t)bp_recombiner_config_t::max_n_of_vars_from_other) {
       fixed_from_guiding = n_vars_from_other - bp_recombiner_config_t::max_n_of_vars_from_other;
       n_vars_from_other  = bp_recombiner_config_t::max_n_of_vars_from_other;
-      thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()};
+      thrust::default_random_engine g{(unsigned int)this->context.problem_ptr->seed_gen.get_seed()};
       thrust::shuffle(a.handle_ptr->get_thrust_policy(),
                       this->remaining_indices.data(),
                       this->remaining_indices.data() + n_different_vars,

cpp/src/mip_heuristics/diversity/recombiners/fp_recombiner.cuh

@@ -53,7 +53,7 @@ class fp_recombiner_t : public recombiner_t<i_t, f_t> {
     i_t n_vars_from_other = n_different_vars;
     if (n_vars_from_other > (i_t)fp_recombiner_config_t::max_n_of_vars_from_other) {
       n_vars_from_other = fp_recombiner_config_t::max_n_of_vars_from_other;
-      thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()};
+      thrust::default_random_engine g{(unsigned int)this->context.problem_ptr->seed_gen.get_seed()};
       thrust::shuffle(a.handle_ptr->get_thrust_policy(),
                       this->remaining_indices.data(),
                       this->remaining_indices.data() + n_different_vars,

cpp/src/mip_heuristics/diversity/recombiners/line_segment_recombiner.cuh

@@ -40,7 +40,7 @@ class line_segment_recombiner_t : public recombiner_t<i_t, f_t> {
     i_t n_vars_from_other = remaining_variables;
     if (n_vars_from_other > (i_t)ls_recombiner_config_t::max_n_of_vars_from_other) {
       n_vars_from_other = ls_recombiner_config_t::max_n_of_vars_from_other;
-      thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()};
+      thrust::default_random_engine g{(unsigned int)this->context.problem_ptr->seed_gen.get_seed()};
       thrust::shuffle(guiding_solution.handle_ptr->get_thrust_policy(),
                       this->remaining_indices.data(),
                       this->remaining_indices.data() + remaining_variables,

cpp/src/mip_heuristics/diversity/recombiners/recombiner.cuh

@@ -119,7 +119,7 @@ class recombiner_t {
                     objective_indices.size());
     if (objective_indices.size() > 0 &&
         objective_indices_in_subproblem.size() < 0.4 * remaining_variables) {
-      std::default_random_engine rng_host(cuopt::seed_generator::get_seed());
+      std::default_random_engine rng_host(context.problem_ptr->seed_gen.get_seed());
       std::vector<i_t> objective_indices_not_in_subproblem;
       std::set_difference(objective_indices.begin(),
                           objective_indices.end(),

cpp/src/mip_heuristics/diversity/recombiners/sub_mip.cuh

@@ -57,7 +57,7 @@ class sub_mip_recombiner_t : public recombiner_t<i_t, f_t> {
     i_t n_vars_from_other = n_different_vars;
     if (n_vars_from_other > (i_t)sub_mip_recombiner_config_t::max_n_of_vars_from_other) {
       n_vars_from_other = sub_mip_recombiner_config_t::max_n_of_vars_from_other;
-      thrust::default_random_engine g{(unsigned int)cuopt::seed_generator::get_seed()};
+      thrust::default_random_engine g{(unsigned int)this->context.problem_ptr->seed_gen.get_seed()};
       thrust::shuffle(a.handle_ptr->get_thrust_policy(),
                       this->remaining_indices.data(),
                       this->remaining_indices.data() + n_different_vars,

cpp/src/mip_heuristics/feasibility_jump/feasibility_jump.cu

@@ -134,7 +134,7 @@ void fj_t<i_t, f_t>::reset_weights(const rmm::cuda_stream_view& climber_stream,
 template <typename i_t, typename f_t>
 void fj_t<i_t, f_t>::randomize_weights(const raft::handle_t* handle_ptr)
 {
-  std::mt19937 rng(cuopt::seed_generator::get_seed());
+  std::mt19937 rng(context.problem_ptr->seed_gen.get_seed());
   constexpr f_t min_weight = 10.;
   constexpr f_t max_weight = 30.;
   // generate a range of weights between 10. and 30.
@@ -671,7 +671,7 @@ void fj_t<i_t, f_t>::run_step_device(const rmm::cuda_stream_view& climber_stream
 
   auto& data    = *climbers[climber_idx];
   auto v        = data.view();
-  settings.seed = cuopt::seed_generator::get_seed();
+  settings.seed = context.problem_ptr->seed_gen.get_seed();
   // ensure an updated copy of the settings is used device-side
   raft::copy(v.settings, &settings, 1, climber_stream);
 

cpp/src/mip_heuristics/feasibility_jump/fj_cpu.cu

@@ -1492,7 +1492,7 @@ static std::unique_ptr<fj_cpu_climber_t<i_t, f_t>> init_fj_cpu_from_host_lp(
   fj_settings.iteration_limit        = std::numeric_limits<int>::max();
   fj_settings.update_weights         = true;
   fj_settings.feasibility_run        = false;
-  fj_settings.seed                   = seed >= 0 ? seed : cuopt::seed_generator::get_seed();
+  fj_settings.seed                   = seed >= 0 ? seed : settings.random_seed;
 
   auto fj_cpu      = std::make_unique<fj_cpu_climber_t<i_t, f_t>>(preemption_flag);
   fj_cpu->view     = typename fj_t<i_t, f_t>::climber_data_t::view_t{};
@@ -1595,13 +1595,13 @@ std::unique_ptr<fj_cpu_climber_t<i_t, f_t>> fj_t<i_t, f_t>::create_cpu_climber(
   init_fj_cpu(*fj_cpu, solution, left_weights, right_weights, objective_weight);
   fj_cpu->settings = settings;
   if (randomize_params) {
-    auto rng                 = std::mt19937(cuopt::seed_generator::get_seed());
+    auto rng                 = std::mt19937(solution.problem_ptr->seed_gen.get_seed());
     fj_cpu->mtm_viol_samples = std::uniform_int_distribution<i_t>(15, 50)(rng);
     fj_cpu->mtm_sat_samples  = std::uniform_int_distribution<i_t>(10, 30)(rng);
     fj_cpu->nnz_samples      = std::uniform_int_distribution<i_t>(2000, 15000)(rng);
     fj_cpu->perturb_interval = std::uniform_int_distribution<i_t>(50, 500)(rng);
   }
-  fj_cpu->settings.seed = cuopt::seed_generator::get_seed();
+  fj_cpu->settings.seed = solution.problem_ptr->seed_gen.get_seed();
   return fj_cpu;  // move
 }
 
@@ -1786,7 +1786,7 @@ std::unique_ptr<fj_cpu_climber_t<i_t, f_t>> init_fj_cpu_standalone(
   std::vector<f_t> default_weights(problem.n_constraints, 1.0);
   init_fj_cpu(*fj_cpu, solution, default_weights, default_weights, 0.0);
   fj_cpu->settings      = settings;
-  fj_cpu->settings.seed = cuopt::seed_generator::get_seed();
+  fj_cpu->settings.seed = solution.problem_ptr->seed_gen.get_seed();
 
   return fj_cpu;
 }

cpp/src/mip_heuristics/feasibility_jump/fj_cpu_worker.cuh

@@ -35,9 +35,9 @@ struct fj_cpu_worker_t {
   ~fj_cpu_worker_t() { stop(); }
 
   // `seed` selects the FJ RNG seed: pass a non-negative value for a deterministic seed,
-  // or -1 to draw from the global cuopt::seed_generator (the historical behavior).
+  // or -1 to fall back to the simplex settings' random_seed.
   // In deterministic mode the caller MUST pass an explicit seed, otherwise the underlying
-  // seed_generator::get_seed() racing with concurrent callers breaks reproducibility.
+  // shared seed source racing with concurrent callers breaks reproducibility.
   void create_worker(const simplex::lp_problem_t<i_t, f_t>& problem,
                      const std::vector<simplex::variable_type_t>& variable_types,
                      const std::vector<f_t>& seed_assignment,

cpp/src/mip_heuristics/local_search/feasibility_pump/feasibility_pump.cu

@@ -51,7 +51,7 @@ feasibility_pump_t<i_t, f_t>::feasibility_pump_t(
     orig_variable_types(context.problem_ptr->n_variables,
                         context.problem_ptr->handle_ptr->get_stream()),
     lp_optimal_solution(lp_optimal_solution_),
-    rng(cuopt::seed_generator::get_seed()),
+    rng(context.problem_ptr->seed_gen.get_seed()),
     timer(20.)
 {
 }

cpp/src/mip_heuristics/local_search/local_search.cu

@@ -39,7 +39,7 @@ local_search_t<i_t, f_t>::local_search_t(mip_solver_context_t<i_t, f_t>& context
        constraint_prop,
        line_segment_search,
        lp_optimal_solution_),
-    rng(cuopt::seed_generator::get_seed()),
+    rng(context.problem_ptr->seed_gen.get_seed()),
     problem_with_objective_cut(*context.problem_ptr, context.problem_ptr->handle_ptr)
 {
   const int n_cpufj = context.settings.heuristic_params.num_cpufj_threads;

cpp/src/mip_heuristics/local_search/rounding/bounds_repair.cu

@@ -31,7 +31,7 @@ bounds_repair_t<i_t, f_t>::bounds_repair_t(const problem_t<i_t, f_t>& pb,
     violated_constraints(0, pb.handle_ptr->get_stream()),
     violated_cstr_map(0, pb.handle_ptr->get_stream()),
     total_vio(pb.handle_ptr->get_stream()),
-    gen(cuopt::seed_generator::get_seed()),
+    gen(pb.seed_gen.get_seed()),
     cycle_vector(MAX_CYCLE_SEQUENCE, -1)
 {
 }

cpp/src/mip_heuristics/local_search/rounding/constraint_prop.cu

@@ -41,7 +41,7 @@ constraint_prop_t<i_t, f_t>::constraint_prop_t(mip_solver_context_t<i_t, f_t>& c
     ub_restore(context.problem_ptr->n_variables, context.problem_ptr->handle_ptr->get_stream()),
     assignment_restore(context.problem_ptr->n_variables,
                        context.problem_ptr->handle_ptr->get_stream()),
-    rng(cuopt::seed_generator::get_seed(), 0, 0)
+    rng(context.problem_ptr->seed_gen.get_seed(), 0, 0)
 {
 }
 
@@ -604,7 +604,7 @@ thrust::pair<f_t, f_t> constraint_prop_t<i_t, f_t>::generate_double_probing_pair
   if (probing_config.has_value()) {
     // for now get the first one
     auto [from_first, from_second] = probing_config.value().get().probing_values[unset_var_idx];
-    std::mt19937 rng(cuopt::seed_generator::get_seed());
+    std::mt19937 rng(context.problem_ptr->seed_gen.get_seed());
     std::uniform_real_distribution<float> dist(0.0f, 1.0f);
     f_t random_value  = dist(rng);
     f_t average_value = (from_first + from_second) / 2;
@@ -848,7 +848,7 @@ bool constraint_prop_t<i_t, f_t>::find_integer(
 {
   using crit_t             = termination_criterion_t;
   auto& unset_integer_vars = unset_vars;
-  std::mt19937 rng(cuopt::seed_generator::get_seed());
+  std::mt19937 rng(context.problem_ptr->seed_gen.get_seed());
   lb_restore.resize(sol.problem_ptr->n_variables, sol.handle_ptr->get_stream());
   ub_restore.resize(sol.problem_ptr->n_variables, sol.handle_ptr->get_stream());
   assignment_restore.resize(sol.problem_ptr->n_variables, sol.handle_ptr->get_stream());
@@ -883,7 +883,7 @@ bool constraint_prop_t<i_t, f_t>::find_integer(
         sol.handle_ptr->get_thrust_policy(),
         unset_integer_vars.begin(),
         unset_integer_vars.begin() + n_to_round,
-        [sol = sol.view(), seed = cuopt::seed_generator::get_seed()] __device__(i_t var_idx) {
+        [sol = sol.view(), seed = context.problem_ptr->seed_gen.get_seed()] __device__(i_t var_idx) {
           raft::random::PCGenerator rng(seed, var_idx, 0);
           auto var_bnd            = sol.problem.variable_bounds[var_idx];
           sol.assignment[var_idx] = round_nearest(sol.assignment[var_idx],

cpp/src/mip_heuristics/local_search/rounding/lb_bounds_repair.cu

@@ -27,7 +27,7 @@ lb_bounds_repair_t<i_t, f_t>::lb_bounds_repair_t(const raft::handle_t* handle_pt
     violated_constraints(0, handle_ptr->get_stream()),
     violated_cstr_map(0, handle_ptr->get_stream()),
     total_vio(handle_ptr->get_stream()),
-    gen(cuopt::seed_generator::get_seed()),
+    gen(problem.seed_gen.get_seed()),
     cycle_vector(MAX_CYCLE_SEQUENCE, -1)
 {
 }

cpp/src/mip_heuristics/local_search/rounding/lb_constraint_prop.cu

@@ -33,7 +33,7 @@ lb_constraint_prop_t<i_t, f_t>::lb_constraint_prop_t(mip_solver_context_t<i_t, f
                    context.problem_ptr->handle_ptr->get_stream()),
     assignment_restore(context.problem_ptr->n_variables,
                        context.problem_ptr->handle_ptr->get_stream()),
-    rng(cuopt::seed_generator::get_seed(), 0, 0)
+    rng(context.problem_ptr->seed_gen.get_seed(), 0, 0)
 {
 }
 
@@ -765,7 +765,7 @@ bool lb_constraint_prop_t<i_t, f_t>::find_integer(
 
   using crit_t             = termination_criterion_t;
   auto& unset_integer_vars = unset_vars;
-  std::mt19937 rng(cuopt::seed_generator::get_seed());
+  std::mt19937 rng(context.problem_ptr->seed_gen.get_seed());
 
   bounds_restore.resize(2 * orig_sol.problem_ptr->n_variables, orig_sol.handle_ptr->get_stream());
   assignment_restore.resize(orig_sol.problem_ptr->n_variables, orig_sol.handle_ptr->get_stream());

cpp/src/mip_heuristics/local_search/rounding/simple_rounding.cu

@@ -112,7 +112,7 @@ void invoke_round_nearest(solution_t<i_t, f_t>& solution)
 
   i_t n_blocks = (solution.problem_ptr->n_integer_vars + TPB - 1) / TPB;
   nearest_rounding_kernel<i_t, f_t><<<n_blocks, TPB, 0, solution.handle_ptr->get_stream()>>>(
-    solution.view(), cuopt::seed_generator::get_seed());
+    solution.view(), solution.problem_ptr->seed_gen.get_seed());
   RAFT_CHECK_CUDA(solution.handle_ptr->get_stream());
 }
 
@@ -127,23 +127,23 @@ void invoke_random_round_nearest(solution_t<i_t, f_t>& solution, i_t n_target_ra
                   solution.problem_ptr->n_integer_vars);
   rmm::device_scalar<i_t> n_randomly_rounded(0, solution.handle_ptr->get_stream());
   random_nearest_rounding_kernel<i_t, f_t><<<n_blocks, TPB, 0, solution.handle_ptr->get_stream()>>>(
-    solution.view(), cuopt::seed_generator::get_seed(), n_randomly_rounded.data());
+    solution.view(), solution.problem_ptr->seed_gen.get_seed(), n_randomly_rounded.data());
   i_t h_n_random_rounds = n_randomly_rounded.value(solution.handle_ptr->get_stream());
   CUOPT_LOG_TRACE("Randomly rounded integers %d", h_n_random_rounds);
   i_t additional_roundings_needed = n_target_random_rounds - h_n_random_rounds;
   if (additional_roundings_needed > 0) {
     // TODO sort the remaining integers with fractionality and round them randomly
     rmm::device_uvector<i_t> shuffled_indices(solution.problem_ptr->integer_indices,
                                               solution.handle_ptr->get_stream());
-    thrust::default_random_engine rng(cuopt::seed_generator::get_seed());
+    thrust::default_random_engine rng(solution.problem_ptr->seed_gen.get_seed());
     // from the remaining integers, populate randomly.
     thrust::shuffle(solution.handle_ptr->get_thrust_policy(),
                     shuffled_indices.begin(),
                     shuffled_indices.end(),
                     rng);
     random_rounding_kernel<i_t, f_t>
       <<<1, 1, 0, solution.handle_ptr->get_stream()>>>(solution.view(),
-                                                       cuopt::seed_generator::get_seed(),
+                                                       solution.problem_ptr->seed_gen.get_seed(),
                                                        shuffled_indices.data(),
                                                        n_randomly_rounded.data(),
                                                        additional_roundings_needed);

cpp/src/mip_heuristics/problem/problem.cuh

@@ -14,6 +14,7 @@
 #include "presolve_data.cuh"
 
 #include <mip_heuristics/logger.hpp>
+#include <utilities/seed_generator.cuh>
 #include <mip_heuristics/relaxed_lp/lp_state.cuh>
 
 #include <cuopt/mathematical_optimization/mip/solver_settings.hpp>
@@ -327,6 +328,8 @@ class problem_t {
   /** name of the objective (only a single objective is currently allowed) */
   std::string objective_name;
   f_t objective_offset;
+  // Seed source for this problem, seeded from the solver settings.
+  seed_generator_t seed_gen;
   bool is_scaled_{false};
   bool preprocess_called{false};
   bool objective_is_integral{false};

cpp/src/mip_heuristics/solution/solution.cu

@@ -229,7 +229,7 @@ template <typename i_t, typename f_t>
 void solution_t<i_t, f_t>::assign_random_within_bounds(f_t ratio_of_vars_to_random_assign,
                                                        bool only_integers)
 {
-  std::mt19937 rng(cuopt::seed_generator::get_seed());
+  std::mt19937 rng(problem_ptr->seed_gen.get_seed());
   auto stream                   = handle_ptr->get_stream();
   std::vector<f_t> h_assignment = host_copy(assignment, stream);
   std::uniform_real_distribution<f_t> unif_prob(0, 1);

cpp/src/mip_heuristics/solve.cu

@@ -370,8 +370,6 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
 
     print_version_info();
 
-    // Initialize seed generator if a specific seed is requested
-    if (settings.seed >= 0) { cuopt::seed_generator::set_seed(settings.seed); }
 
     raft::common::nvtx::range fun_scope("Running solver");
     auto timer = timer_t(time_limit);
@@ -434,6 +432,7 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
     bool has_symmetry = false;
     if (settings.symmetry != 0) {
       mip::problem_t<i_t, f_t> problem(op_problem);
+        if (settings.seed >= 0) { problem.seed_gen.set_seed(settings.seed); }
       simplex_solver_settings_t<i_t, f_t> simplex_settings;
       simplex_settings.set_log(true);
       simplex_settings.time_limit = settings.time_limit;
@@ -453,6 +452,7 @@ mip_solution_t<i_t, f_t> solve_mip_helper(optimization_problem_t<i_t, f_t>& op_p
     std::optional<mip::third_party_presolve_device_result_t<i_t, f_t>> presolve_result_opt;
     mip::problem_t<i_t, f_t> problem(
       op_problem, settings.get_tolerances(), settings.determinism_mode == CUOPT_MODE_DETERMINISTIC);
+      if (settings.seed >= 0) { problem.seed_gen.set_seed(settings.seed); }
 
     auto run_presolve              = settings.presolver != presolver_t::None;
     bool has_set_solution_callback = false;