Correct sampler lifecycle and Bayesian posterior evaluation

.github/workflows/ci.yml

@@ -0,0 +1,14 @@
+name: CI
+on: [push, pull_request]
+permissions:
+  contents: read
+jobs:
+  regression:
+    strategy:
+      matrix:
+        os: [ubuntu-latest, macos-latest]
+    runs-on: ${{ matrix.os }}
+    steps:
+      - uses: actions/checkout@v4
+      - run: c++ -std=c++17 -O2 -Wall -Wextra -Wpedantic -pthread -Ibest tests/regression.cpp -o regression && ./regression
+      - run: c++ -std=c++17 -O1 -g -fsanitize=address,undefined -fno-omit-frame-pointer -pthread -Ibest tests/regression.cpp -o sanitized && ./sanitized

.gitignore

@@ -1,67 +1,4 @@
-# Compiled Object files
-*.slo
-*.lo
-*.o
-*.obj
-
-# Precompiled Headers
-*.gch
-*.pch
-
-# Compiled Dynamic libraries
-*.so
-*.dylib
-*.dll
-
-# Fortran module files
-*.mod
-*.smod
-
-# Compiled Static libraries
-*.lai
-*.la
-*.a
-*.lib
-
-# Executables
-*.exe
-*.out
-*.app
-
-# OSX / XCode
+build*/
+install/
+bazel-*
 .DS_Store
-*.pbxuser
-xcuserdata/
-
-# Visual Studio user-specific files
-*.suo
-*.user
-*.userosscache
-*.sln.docstates
-*.vcxproj.filters
-
-# Visual Studio cache/options directory
-.vs/
-
-# Visual C++ cache files
-ipch/
-*.aps
-*.ncb
-*.opendb
-*.opensdf
-*.sdf
-*.cachefile
-*.VC.db
-*.VC.VC.opendb
-
-# Build results
-[Dd]ebug/
-[Dd]ebugPublic/
-[Rr]elease/
-[Rr]eleases/
-x64/
-x86/
-bld/
-[Bb]in/
-[Oo]bj/
-[Ll]og/

CHANGELOG.md

@@ -0,0 +1,13 @@
+# Changelog
+
+## Unreleased
+
+- Replace speculative parallel proposals with sequential adaptive Metropolis-within-Gibbs. Remove shared random-generator races and proposal-selection bias. Remove the thread-count argument from `AMWG::Init`; callers now pass only the starting state and log-posterior callback.
+- `Sample(n)` records exactly n completed sweeps; `Burn(n)` advances without retaining samples. Adaptation tracks sweeps independently of chain storage. Reinitialization resets the seeded generator and adaptation. Sampling no longer prints progress.
+- Validate initialization, observations, batch sizes and statistical parameters. Callback exceptions propagate on the calling thread; a partially completed sweep may have advanced the state, but is not added to the chain.
+- Enforce both sigma bounds and nu >= 1, evaluate the model in log space, retain double precision, and reuse Student-t normalization across observations.
+- Own model observations to prevent dangling references. BEST cannot be copied or moved because its callback refers to the model.
+- Define the empirical 95% interval as the shortest window containing ceil(0.95*n) observations, choosing the first window on ties. Reject empty/nonfinite inputs.
+- Compute population standard deviation without a temporary allocation. Add numerical, lifecycle, seeded-distribution and independent-chain regression checks.
+
+These corrections change seeded sequences, inference results and invalid-input behavior. Old results should be recomputed. Reproducibility is within the same implementation and standard library; C++ random-distribution algorithms are not portable bit-for-bit.

best/amwg.hpp

@@ -1,100 +1,79 @@
 #pragma once
-
+#include <algorithm>
 #include <array>
-#include <atomic>
-#include <condition_variable>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
 #include <functional>
+#include <limits>
 #include <random>
-#include <thread>
+#include <stdexcept>
 #include <type_traits>
+#include <utility>
 #include <vector>
 
 static const size_t kDefaultBatchSize = 50;
 
-/**
- * Implementation of an Adaptive Metropolis within Gibbs sampler.
- */
+// Sequential adaptive Metropolis-within-Gibbs; one proposal per coordinate.
+// Instances are not safe for concurrent mutation. Independent chains may run
+// on separate threads with independent seeds and posterior state.
 template<typename RealType, size_t NumParams>
 class AMWG {
+  static_assert(std::is_floating_point<RealType>::value, "floating point required");
+  static_assert(NumParams > 0, "at least one parameter required");
 public:
-
   using ParamArray = std::array<RealType, NumParams>;
   using PosteriorFunc = std::function<RealType(const ParamArray&)>;
-
-  AMWG(size_t batchSize = kDefaultBatchSize, uint32_t seed = std::mt19937::default_seed)
-    : rng_(seed)
-    , rand_(0.0, 1.0)
-    , currentPosteriorDensity_(0.0)
-    , batchSize_(batchSize)
-    , batchCount_(0)
-    , logSD_()
-    , acceptanceCount_()
-    , running_(false)
-    , epoch_(0)
-    , completedCount_(0)
-    , curParam_(0)
-  {};
-
-  ~AMWG() {
-    {
-      std::lock_guard<std::mutex> lock(mutex_);
-      running_ = false;
-      epoch_ = SIZE_MAX;
-    }
-
-    cv_.notify_all();
-
-    for (auto&& t : threads_) {
-      if (t.joinable()) t.join();
-    }
+  explicit AMWG(size_t batchSize = kDefaultBatchSize,
+                uint32_t seed = std::mt19937::default_seed)
+    : rng_(seed), seed_(seed), batchSize_(batchSize) {
+    if (!batchSize) throw std::invalid_argument("batch size must be positive");
   }
-
-  void Init(std::array<RealType, NumParams> startValues, PosteriorFunc posteriorFunc, uint32_t threads = std::thread::hardware_concurrency()) {
-    state_ = startValues;
-    posteriorFunc_ = posteriorFunc;
+  void Init(ParamArray start, PosteriorFunc posterior) {
+    if (!posterior) throw std::invalid_argument("posterior is required");
+    for (auto x : start)
+      if (!std::isfinite(x)) throw std::invalid_argument("nonfinite initial state");
+    RealType density = posterior(start);
+    if (!std::isfinite(density)) throw std::invalid_argument("initial log density must be finite");
+    posterior_ = std::move(posterior);
+    state_ = start;
+    density_ = density;
     chain_.clear();
-    currentPosteriorDensity_ = posteriorFunc_(state_);
-
-    running_ = true;
-
-    threads_.reserve(threads);
-    threadStates_.resize(threads);
-
-    for (uint32_t threadId = 0; threadId < threads; threadId++) {
-      threads_.emplace_back(std::thread(&AMWG::threadWorker, this, threadId));
-    }
-  };
-
-  size_t NextSample();
-  void Sample(size_t n);
-  void Burn(size_t n);
-  RealType posterior_density() { return currentPosteriorDensity_; };
-  std::vector<ParamArray>& chain() { return chain_; };
-
+    rng_.seed(seed_);
+    logSD_.fill(0);
+    accepted_.fill(0);
+    batchCount_ = withinBatch_ = 0;
+    initialized_ = true;
+  }
+  size_t NextSample() {
+    requireInit();
+    chain_.reserve(chain_.size() + 1);
+    step();
+    chain_.push_back(state_);
+    return NumParams;
+  }
+  void Sample(size_t n) {
+    requireInit();
+    if (n > chain_.max_size() - chain_.size()) throw std::length_error("chain too large");
+    chain_.reserve(chain_.size() + n);
+    for (size_t i=0; i<n; ++i) { step(); chain_.push_back(state_); }
+  }
+  void Burn(size_t n) { requireInit(); for (size_t i=0; i<n; ++i) step(); }
+  RealType posterior_density() const { requireInit(); return density_; }
+  const ParamArray& state() const { requireInit(); return state_; }
+  std::vector<ParamArray>& chain() { return chain_; }
+  const std::vector<ParamArray>& chain() const { return chain_; }
 private:
-
-  void threadWorker(uint32_t threadId);
-  void createProposal(uint32_t threadId);
-
-  PosteriorFunc posteriorFunc_;
+  void requireInit() const { if (!initialized_) throw std::logic_error("Init must be called first"); }
+  void step();
   std::mt19937 rng_;
-  std::uniform_real_distribution<RealType> rand_;
+  uint32_t seed_;
+  size_t batchSize_, batchCount_=0, withinBatch_=0;
+  bool initialized_=false;
+  PosteriorFunc posterior_;
+  ParamArray state_{}, logSD_{};
+  std::array<size_t, NumParams> accepted_{};
+  RealType density_=0;
   std::vector<ParamArray> chain_;
-  RealType currentPosteriorDensity_;
-  size_t batchSize_;
-  size_t batchCount_;
-  ParamArray state_;
-  ParamArray logSD_;
-  std::array<size_t, NumParams> acceptanceCount_;
-
-  std::atomic_bool running_;
-  std::atomic<size_t> epoch_;
-  std::atomic<size_t> completedCount_;
-  std::atomic<size_t> curParam_;
-  std::vector<std::thread> threads_;
-  std::vector<std::tuple<RealType, RealType, bool>> threadStates_;
-  std::mutex mutex_;
-  std::condition_variable cv_;
 };
-
 #include "amwg.inl"