From 072e211e900450d54beed29b43a80f95573acc2c Mon Sep 17 00:00:00 2001
From: Dominik Demuth <dominik.demuth@physik.tu-darmstadt.de>
Date: Sun, 8 Mar 2026 15:02:47 +0100
Subject: [PATCH] Dynamic scheduling of threads

---
 src/sims.cpp            | 26 +++++++++++++++++++-------
 src/sims.h              |  4 ++--
 src/utils/functions.cpp |  7 +++++++
 3 files changed, 28 insertions(+), 9 deletions(-)

diff --git a/src/sims.cpp b/src/sims.cpp
index c1948f4..6ac53e4 100644
--- a/src/sims.cpp
+++ b/src/sims.cpp
@@ -53,9 +53,11 @@ void run_simulation(std::vector<Experiment *> experiments,
     auto &local_dynamics = *thread_dynamics[tid];
     auto &local_experiments = thread_experiments[tid];
 
-#pragma omp for schedule(static)
+    Trajectory traj;
+
+#pragma omp for schedule(dynamic, 16)
     for (int mol_i = 0; mol_i < num_walker; mol_i++) {
-      auto traj = make_trajectory(local_dynamics, tmax, local_rng);
+      make_trajectory(traj, local_dynamics, tmax, local_rng);
       const double init_omega = local_dynamics.getInitOmega();
 
       for (auto &exp : local_experiments) {
@@ -86,8 +88,21 @@ void run_simulation(std::vector<Experiment *> experiments,
   printEnd(start);
 }
 
-Trajectory make_trajectory(Dynamics &dynamics, const double t_max,
-                           std::mt19937_64 &rng) {
+void make_trajectory(Trajectory &traj, Dynamics &dynamics, const double t_max,
+                     std::mt19937_64 &rng) {
+  traj.time.clear();
+  traj.phase.clear();
+  traj.omega.clear();
+
+  // Reserve is only useful on first call per thread; after that
+  // clear() preserves capacity from the previous walker
+  if (traj.time.capacity() == 0) {
+    const size_t estimate = 1024;
+    traj.time.reserve(estimate);
+    traj.phase.reserve(estimate);
+    traj.omega.reserve(estimate);
+  }
+
   double t_passed = 0;
   double phase = 0;
 
@@ -95,7 +110,6 @@ Trajectory make_trajectory(Dynamics &dynamics, const double t_max,
 
   double omega = dynamics.getInitOmega();
 
-  Trajectory traj;
   traj.time.emplace_back(t_passed);
   traj.phase.emplace_back(phase);
   traj.omega.emplace_back(omega);
@@ -110,8 +124,6 @@ Trajectory make_trajectory(Dynamics &dynamics, const double t_max,
     traj.phase.emplace_back(phase);
     traj.omega.emplace_back(omega);
   }
-
-  return traj;
 }
 
 std::chrono::system_clock::time_point
diff --git a/src/sims.h b/src/sims.h
index 5285e3e..41e54c9 100644
--- a/src/sims.h
+++ b/src/sims.h
@@ -15,8 +15,8 @@ void run_simulation(std::vector<Experiment *> experiments,
                     const std::unordered_map<std::string, double> &optional,
                     Dynamics &dynamics, std::mt19937_64 &rng);
 
-Trajectory make_trajectory(Dynamics &dynamics, double t_max,
-                           std::mt19937_64 &rng);
+void make_trajectory(Trajectory &traj, Dynamics &dynamics, double t_max,
+                     std::mt19937_64 &rng);
 
 std::chrono::system_clock::time_point
 printStart(const std::unordered_map<std::string, double> &optional);
diff --git a/src/utils/functions.cpp b/src/utils/functions.cpp
index bfa1799..ae62913 100644
--- a/src/utils/functions.cpp
+++ b/src/utils/functions.cpp
@@ -1,3 +1,4 @@
+#include <algorithm>
 #include <chrono>
 #include <iostream>
 #include <vector>
@@ -5,6 +6,12 @@
 int nearest_index(const std::vector<double> &x_ref, const double x,
                   int start = 0) {
   const int last = static_cast<int>(x_ref.size()) - 2;
+  if (start == 0 && last > 32) {
+    // Binary search for cold starts on large trajectories
+    auto it = std::upper_bound(x_ref.begin(), x_ref.end(), x);
+    int idx = static_cast<int>(it - x_ref.begin()) - 1;
+    return std::clamp(idx, 0, last);
+  }
   while (start < last && x > x_ref[start + 1]) {
     start++;
   }