perf(kernel): persistent worker pool instead of threads per batch (6% faster)

.cspell/custom-dictionary-workspace.txt

@@ -25,6 +25,7 @@ armhf
 armv
 ASHP
 asyncio
+atfork
 authoriser
 autodocstring
 autoflake
@@ -380,6 +381,7 @@ preseeded
 prevs
 protobuf
 psum
+pthread
 pvbat
 pvenergytotal
 pvlib
@@ -537,10 +539,14 @@ venv
 Victron
 visibilitychange
 waf
+waitpid
 Wallbox
 weblink
 welink
 Werror
+WEXITSTATUS
+WIFEXITED
+WNOHANG
 workmode
 writeonly
 wrongsha

apps/predbat/plan.py

@@ -35,6 +35,26 @@
 PLAN_PASS_WINDOW_BUDGET = 8
 
 
+def resolve_batch_threads(threads, cpu_count_value):
+    """Map the threads setting onto how many kernel lanes one batch may use.
+
+    'auto' takes the core count and is deliberately not capped. On a fast machine the curve is very
+    flat and peaks slightly below the core count - measured on the 20-scenario benchmark, best of 3:
+    serial 26.33s, 4 threads 24.89s, 6 threads 24.71s, 8 threads 24.91s, 16 threads 25.04s - so a cap
+    looks attractive. But re-running with each job made eight times dearer, which is how a machine
+    where the kernel dominates behaves, the curve stops turning over entirely: 48.92s serial, 32.03s
+    at 4, 29.98s at 6, 29.85s at 8, 28.94s at 16.
+
+    That makes the risk asymmetric. Capping at 4 costs 0.7% on the fast machine but 10.7% on the
+    kernel-heavy one, while not capping costs 1.3% at worst. The worst case for a low cap is far
+    worse than the worst case for none, so 'auto' is left alone and anyone who wants fewer lanes sets
+    threads: explicitly.
+    """
+    if threads == "auto":
+        return max(cpu_count_value, 1)
+    return max(int(threads), 1)
+
+
 def slots_around(target_slots, slot_lengths):
     """
     Return a list of slot lengths around the target slots
@@ -1306,11 +1326,7 @@ def calculate_plan(self, recompute=True, debug_mode=False, publish=True):
         # The kernel spreads one batched fan-out across threads with the GIL released for the whole
         # call, so these are real cores - unlike a Python ThreadPool, which peaked at 1.15x on two
         # threads and then degraded below serial (perf/threadpool-prototype).
-        threads = self.get_arg("threads", "auto")
-        if threads == "auto":
-            self.prediction.batch_threads = cpu_count()
-        else:
-            self.prediction.batch_threads = max(int(threads), 1)
+        self.prediction.batch_threads = resolve_batch_threads(self.get_arg("threads", "auto"), cpu_count())
         self.log("Prediction batch using {} kernel thread(s)".format(self.prediction.batch_threads))
         kernel_message, kernel_is_warning = kernel_status_summary(self.prediction)
         self.log("{}Prediction kernel: {}".format("Warn: " if kernel_is_warning else "", kernel_message))

apps/predbat/prediction_kernel.cpp

@@ -16,6 +16,8 @@
 // Build: bash apps/predbat/build_kernel.sh (g++/clang, C++17, no dependencies)
 
 #include <algorithm>
+#include <atomic>
+#include <condition_variable>
 #include <cstdint>
 #include <limits>
 #include <cstdio>
@@ -28,6 +30,10 @@
 #include <thread>
 #include <vector>
 
+#if defined(__unix__) || defined(__APPLE__)
+#include <pthread.h>
+#endif
+
 #define PK_ABI_VERSION 3
 #define PK_PARITY_REVISION 5
 #define PK_MAX_CARS 4
@@ -299,12 +305,47 @@ struct PkBatchResult {
 // each scenario would otherwise heap-allocate ~70KB, and 1200 of those per batch put every worker
 // thread in contention on the allocator. Measured 4 threads at 0.6x of serial before this, because
 // the simulation is cheap enough that malloc, not arithmetic, was the bottleneck.
+//
+// See thread_scratch() below: one of these now serves every batch a thread runs, not one per call.
+
+// Counts PkScratch constructions for the life of the library, for pk_scratch_construct_count(). The
+// suite asserts this stops growing with the number of batches rather than trusting the reuse.
+std::atomic<int64_t> g_scratch_constructions{0};
+
+// Counts worker threads started for the life of the library, for pk_worker_thread_count(). Thread
+// creation measured 7.5us each and a benchmark plan used to start 201,545 of them - 1.5s of a 26s
+// plan, which is what the pool exists to stop paying. The suite asserts this stops growing.
+std::atomic<int64_t> g_worker_threads_started{0};
+
 struct PkScratch {
     std::vector<int32_t> charge_member, export_member;
     std::vector<int32_t> clipped_start, clipped_end;
     std::vector<double> clipped_limit;
+
+    PkScratch()
+    {
+        g_scratch_constructions.fetch_add(1, std::memory_order_relaxed);
+    }
 };
 
+// One scratch per thread, reused for that thread's lifetime.
+//
+// Every field is fully rewritten before it is read on each run - build_window_membership assigns
+// exactly n_steps entries and the clip helpers clear-then-push - which is the same invariant that
+// already let one scratch serve a whole stride within a batch. Extending it across calls means a
+// thread grows these buffers once instead of once per call: a benchmark plan issues ~25,000 batch
+// calls, and at a median of 6 jobs each the ~70KB allocation was the same order as the simulation it
+// wrapped.
+//
+// Thread-local rather than per-context because nothing here survives a simulation, so it never needs
+// to belong to a particular context. A thread that later runs a context with a different n_steps
+// simply has the vectors reassigned to the new length.
+static PkScratch &thread_scratch()
+{
+    static thread_local PkScratch scratch;
+    return scratch;
+}
+
 // Deep-copied context storage so Python-side buffers can be freed after create
 struct ContextStore {
     std::vector<double> rate_import, rate_export, alert_keep;
@@ -1276,8 +1317,7 @@ int32_t pk_run(int64_t handle, const PkScenario *s, PkResult *out)
     if (!c) {
         return 1;
     }
-    PkScratch scratch;
-    return pk_run_one(c, s, out, -1, -1, nullptr, nullptr, scratch);  // c is the store
+    return pk_run_one(c, s, out, -1, -1, nullptr, nullptr, thread_scratch());  // c is the store
 }
 
 // Run n_jobs scenarios against one context in a single call.
@@ -1312,6 +1352,170 @@ static void run_batch_job(const ContextStore *c, const PkBatchJob &job, PkBatchR
     out.status = pk_run_one(c, &scenario, &out.result, job.soc_range_start_step, job.soc_range_end_step, &out.soc_range_min, &out.soc_range_max, scratch);
 }
 
+
+// A pool of worker threads that outlives every batch and every context.
+//
+// pk_run_batch used to build and join a fresh std::vector<std::thread> per call. Measured over one
+// benchmark plan that was 201,545 thread creations costing 1.5s - against 3.1s of simulation, so the
+// creation cost was cancelling most of what the parallelism won and threading measured 1.6% overall.
+// Parking workers on a condition variable replaces each creation with a wake.
+//
+// Workers park until `generation` changes, then take lane `index` of a stride partition. The calling
+// thread takes lane 0 itself, so a six-job batch wakes five workers rather than six and the caller
+// does useful work instead of blocking. Nobody holds the mutex while simulating; run_batch_job only
+// reads the shared const context and writes its own result slot, which is what made the previous
+// per-call threading safe and is unchanged here.
+class PkThreadPool {
+  public:
+    // Start workers until at least `want` exist. Grows, never shrinks. Returns how many there are.
+    int32_t reserve_workers(int32_t want)
+    {
+        std::lock_guard<std::mutex> lock(mutex);
+        while (static_cast<int32_t>(workers.size()) < want) {
+            const int32_t lane = static_cast<int32_t>(workers.size()) + 1;  // lane 0 is the caller
+            std::unique_ptr<Worker> worker(new Worker());
+            try {
+                worker->thread = std::thread(&PkThreadPool::worker_loop, this, lane, worker.get());
+            } catch (const std::system_error &) {
+                // Out of thread resources - keep what we have and let the caller use fewer lanes.
+                break;
+            }
+            workers.push_back(std::move(worker));
+            g_worker_threads_started.fetch_add(1, std::memory_order_relaxed);
+        }
+        return static_cast<int32_t>(workers.size());
+    }
+
+    // Run jobs[0..n_jobs) across `use` lanes, the calling thread taking lane 0. `use` must be at
+    // least 1 and no more than worker_count() + 1.
+    void run(const ContextStore *c, const PkBatchJob *jobs, int32_t n_jobs, PkBatchResult *results, int32_t use)
+    {
+        {
+            std::lock_guard<std::mutex> lock(mutex);
+            cur_ctx = c;
+            cur_jobs = jobs;
+            cur_results = results;
+            cur_n_jobs = n_jobs;
+            cur_use = use;
+            outstanding = use - 1;
+            // Each worker has its own flag and its own condition variable, so exactly the lanes this
+            // batch needs are woken. Broadcasting instead would wake every worker in the pool - with
+            // a median batch of six jobs against sixteen workers that was ~250,000 wakeups per plan
+            // spent entirely on threads that immediately went back to sleep.
+            for (int32_t lane = 1; lane < use; lane++) {
+                workers[lane - 1]->has_work = true;
+            }
+        }
+        for (int32_t lane = 1; lane < use; lane++) {
+            workers[lane - 1]->ready.notify_one();
+        }
+
+        PkScratch &scratch = thread_scratch();
+        for (int32_t i = 0; i < n_jobs; i += use) {
+            run_batch_job(c, jobs[i], results[i], scratch);
+        }
+
+        std::unique_lock<std::mutex> lock(mutex);
+        work_done.wait(lock, [this] { return outstanding == 0; });
+    }
+
+  private:
+    // One parked worker. The flag is written under the pool mutex before its condition variable is
+    // signalled, so a wake cannot be lost even if the worker has not re-parked yet.
+    struct Worker {
+        std::condition_variable ready;
+        bool has_work = false;
+        std::thread thread;
+    };
+
+    void worker_loop(int32_t lane, Worker *self)
+    {
+        std::unique_lock<std::mutex> lock(mutex);
+        for (;;) {
+            self->ready.wait(lock, [self] { return self->has_work; });
+            self->has_work = false;
+            const ContextStore *c = cur_ctx;
+            const PkBatchJob *jobs = cur_jobs;
+            PkBatchResult *results = cur_results;
+            const int32_t n_jobs = cur_n_jobs;
+            const int32_t use = cur_use;
+            lock.unlock();
+
+            PkScratch &scratch = thread_scratch();
+            for (int32_t i = lane; i < n_jobs; i += use) {
+                run_batch_job(c, jobs[i], results[i], scratch);
+            }
+
+            lock.lock();
+            if (--outstanding == 0) {
+                work_done.notify_one();
+            }
+        }
+    }
+
+    std::mutex mutex;
+    std::condition_variable work_done;  // the caller waits here
+    const ContextStore *cur_ctx = nullptr;
+    const PkBatchJob *cur_jobs = nullptr;
+    PkBatchResult *cur_results = nullptr;
+    int32_t cur_n_jobs = 0;
+    int32_t cur_use = 1;
+    int32_t outstanding = 0;
+    std::vector<std::unique_ptr<Worker>> workers;
+};
+
+// The pool is created on first threaded use and deliberately never destroyed: a static destructor
+// running while CPython tears the process down risks a deadlock for no benefit, and the OS reclaims
+// the threads at exit.
+std::atomic<PkThreadPool *> g_pool{nullptr};
+std::mutex g_pool_create_mutex;
+// Held for a whole batch. Python can enter pk_run_batch from two threads at once because ctypes
+// releases the GIL, and one shared pool cannot serve two batches concurrently - its published job
+// state would race. Nothing in predbat does this today; serialising is the safe answer if it ever does.
+std::mutex g_dispatch_mutex;
+
+#if defined(__unix__) || defined(__APPLE__)
+// Threads do not survive fork(), so a child would inherit a pool whose workers do not exist and hang
+// waiting for them. The prepare/parent pair keep the dispatch mutex consistent across the fork, and
+// the child drops the pool so its next batch builds a fresh one. Nothing here allocates or takes a
+// new lock, which is what the child handler is allowed to do; the dead pool leaks in the child.
+void pk_atfork_prepare()
+{
+    g_dispatch_mutex.lock();
+}
+
+void pk_atfork_parent()
+{
+    g_dispatch_mutex.unlock();
+}
+
+void pk_atfork_child()
+{
+    g_dispatch_mutex.unlock();
+    g_pool.store(nullptr, std::memory_order_relaxed);
+}
+#endif
+
+// Fetch the pool, creating it on first use.
+static PkThreadPool *pool_singleton()
+{
+    PkThreadPool *pool = g_pool.load(std::memory_order_acquire);
+    if (pool) {
+        return pool;
+    }
+    std::lock_guard<std::mutex> lock(g_pool_create_mutex);
+    pool = g_pool.load(std::memory_order_relaxed);
+    if (!pool) {
+#if defined(__unix__) || defined(__APPLE__)
+        static std::once_flag atfork_once;
+        std::call_once(atfork_once, [] { pthread_atfork(pk_atfork_prepare, pk_atfork_parent, pk_atfork_child); });
+#endif
+        pool = new PkThreadPool();
+        g_pool.store(pool, std::memory_order_release);
+    }
+    return pool;
+}
+
 // Run n_jobs scenarios against one context in a single call.
 //
 // The context lookup and the Python/C boundary are paid once for a whole fan-out rather than once
@@ -1333,30 +1537,50 @@ int32_t pk_run_batch(int64_t handle, const PkBatchJob *jobs, int32_t n_jobs, PkB
     // Scenarios are independent - each reads the shared const context and writes only its own result
     // slot - so they are split across threads by stride.
     if (n_threads > 1 && n_jobs > 1) {
-        const int32_t use = n_threads < n_jobs ? n_threads : n_jobs;
-        std::vector<std::thread> pool;
-        pool.reserve(use);
-        for (int32_t t = 0; t < use; t++) {
-            pool.emplace_back([c, jobs, results, n_jobs, use, t]() {
-                PkScratch scratch;
-                for (int32_t i = t; i < n_jobs; i += use) {
-                    run_batch_job(c, jobs[i], results[i], scratch);
-                }
-            });
+        std::lock_guard<std::mutex> dispatch(g_dispatch_mutex);
+        int32_t use = n_threads < n_jobs ? n_threads : n_jobs;
+        PkThreadPool *pool = pool_singleton();
+        // Clamp to the lanes that actually exist: reserve_workers stops early if the system refuses a
+        // thread, and a batch must never wait on a worker that was never started.
+        const int32_t lanes = pool->reserve_workers(use - 1) + 1;
+        if (use > lanes) {
+            use = lanes;
         }
-        for (auto &th : pool) {
-            th.join();
+        if (use > 1) {
+            pool->run(c, jobs, n_jobs, results, use);
+            return 0;
         }
-        return 0;
+        // No workers available - fall through and run the batch inline.
     }
-    PkScratch scratch;
     for (int32_t i = 0; i < n_jobs; i++) {
-        run_batch_job(c, jobs[i], results[i], scratch);
+        run_batch_job(c, jobs[i], results[i], thread_scratch());
     }
     return 0;
 }
 
 
+// Test hook: how many PkScratch instances have been constructed since the library was loaded.
+//
+// The scratch buffers are the only heap allocation a simulation makes, and the kernel's whole
+// batching win rests on not paying them per call - so the suite asserts this counter stops growing
+// with the number of batches rather than trusting the code to be reusing them.
+int64_t pk_scratch_construct_count(void)
+{
+    return g_scratch_constructions.load(std::memory_order_relaxed);
+}
+
+
+// Test hook: how many worker threads have been started since the library was loaded.
+//
+// Batches used to build and join a fresh set of threads on every call. The suite asserts this count
+// stops growing once the pool is warm, because a thread creation costs several microseconds against
+// a median batch of six jobs - it was cancelling the parallelism it was there to provide.
+int64_t pk_worker_thread_count(void)
+{
+    return g_worker_threads_started.load(std::memory_order_relaxed);
+}
+
+
 // Test hook: sweep SoC densely and confirm the precomputed bucket boundaries give exactly the same
 // percent as the round_py path they replace. Returns the number of disagreements (0 = equivalent).
 int32_t pk_verify_soc_percent_table(double soc_max, int32_t samples)