Native particle MCMC (SMC/PG), revised: dropping AdvancedPS

AGENTS.md

@@ -38,7 +38,7 @@ CI matrix: Julia stable + min, Ubuntu/Windows/macOS, 1 and 2 threads.
 
 Most complexity is in DynamicPPL. Turing.jl contains:
 
-  - **Sampler implementations** (`src/mcmc/`): HMC/NUTS/HMCDA (wrapping AdvancedHMC), MH (wrapping AdvancedMH), particle samplers SMC/PG/CSMC (wrapping AdvancedPS), ESS (wrapping EllipticalSliceSampling), SGLD/SGHMC, Emcee, and Gibbs.
+  - **Sampler implementations** (`src/mcmc/`): HMC/NUTS/HMCDA (wrapping AdvancedHMC), MH (wrapping AdvancedMH), particle samplers SMC/PG/CSMC (implemented natively in `particle_mcmc.jl` on top of Libtask coroutines), ESS (wrapping EllipticalSliceSampling), SGLD/SGHMC, Emcee, and Gibbs.
   - **External sampler interface** (`src/mcmc/external_sampler.jl`): The `externalsampler()` wrapper lets any `AbstractMCMC.AbstractSampler` that implements `step` for `LogDensityModel` work with Turing models. This is the easier path for new samplers — it only requires a dependency on AbstractMCMC and the LogDensityProblems.jl interface, with no Turing internals. The tradeoff is less power: you can only interact with the model as a black-box log-density function, just like using `LogDensityFunction` directly.
   - **Variational inference** (`src/variational/`): Wraps AdvancedVI algorithms.
   - **Mode estimation** (`src/optimisation/`): MAP and MLE via Optimization.jl.

HISTORY.md

@@ -1,3 +1,24 @@
+# 0.47.0
+
+## Breaking changes
+
+### Particle MCMC (SMC and PG)
+
+`SMC` and `PG` / `CSMC` have been reimplemented natively and no longer depend on AdvancedPS.
+
+Resampling schemes are now types rather than functions — `StratifiedResampler()`, `SystematicResampler()`, and `MultinomialResampler()` (in `Turing.Inference`), optionally wrapped in `ESSThresholdResampler(threshold, scheme)` to resample only when the effective sample size falls below `threshold * nparticles`; for example `SMC(Turing.Inference.SystematicResampler())`, `SMC(0.5)`, or `PG(10, Turing.Inference.MultinomialResampler(), 0.5)`.
+The old function-based API (`resample_systematic`, `AdvancedPS.ResampleWithESSThreshold`, …) is gone.
+
+The default scheme is now **stratified** rather than systematic: it stays consistent as the number of particles grows, which systematic does not.
+The selected scheme applies to unconditional sweeps; `PG` / `CSMC` draw the ancestors of a conditional sweep from the categorical over the weights, since a correct conditional version of stratified or systematic resampling is scheme-specific rather than "pin one draw and keep the rest".
+Exact draws may therefore differ from previous releases, but remain statistically consistent (the same target distribution).
+
+The rewrite also brings:
+
+  - **Reproducibility.** Internal seeds are derived through a counter-based (Philox) generator, so a fixed user seed gives the same draws on every Julia version and platform, and splitting one stream into many is better decorrelated. Previously, results could drift between Julia versions even under a `StableRNG` (https://github.com/TuringLang/Turing.jl/issues/2781).
+  - **Parallelism** at two independent levels. *Across chains*, SMC/PG work with AbstractMCMC's `MCMCThreads()` / `MCMCDistributed()` like any other sampler — each chain is an independent run. *Within a single sweep*, `SMC(; multithreaded=true)` / `PG(n; multithreaded=true)` spread that sweep's particles across threads. These are separate knobs: the ensemble does not parallelise a sweep, `multithreaded` does not parallelise chains, and they compose. Neither changes the results; start Julia with multiple threads (e.g. `julia -t auto`) for the thread-based paths to take effect.
+  - **Equal-weight draws & diagnostics.** `SMC` resamples once at the end of the sweep so the returned particles are an equal-weight sample — `mean(chain[...])` and other summaries need no weighting. `SMC`, `PG`, and `CSMC` chains all carry `log_normalizing_constant`; `SMC` chains additionally carry `ess_per_step`, the per-observation effective sample size across the sweep (a degeneracy diagnostic). For `SMC`, `exp(log_normalizing_constant)` is an unbiased estimator of the marginal likelihood `p(y)` under the usual particle-filter assumptions. For `PG` / `CSMC`, the exponentiated estimate is biased and must not be used for model comparison; see the `PG` docstring.
+
 # 0.46.1
 
 Fixed a bug, present since v0.41.0, that biased `PG` / `CSMC` posteriors, whether sampled on their own or as a Gibbs component.

Project.toml

@@ -1,6 +1,6 @@
 name = "Turing"
 uuid = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
-version = "0.46.1"
+version = "0.47.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -9,7 +9,6 @@ AbstractPPL = "7a57a42e-76ec-4ea3-a279-07e840d6d9cf"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d"
 AdvancedMH = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170"
-AdvancedPS = "576499cb-2369-40b2-a588-c64705576edc"
 AdvancedVI = "b5ca4192-6429-45e5-a2d9-87aec30a685c"
 BangBang = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
 Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
@@ -29,6 +28,7 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Random123 = "74087812-796a-5b5d-8853-05524746bad3"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -52,7 +52,6 @@ AbstractPPL = "0.15"
 Accessors = "0.1"
 AdvancedHMC = "0.8.3"
 AdvancedMH = "0.8.9"
-AdvancedPS = "0.7.2"
 AdvancedVI = "0.7"
 BangBang = "0.4.2"
 Bijectors = "0.15.17, 0.16"
@@ -74,6 +73,7 @@ OptimizationOptimJL = "0.1 - 0.4"
 OrderedCollections = "1, 2"
 Printf = "1"
 Random = "1"
+Random123 = "1.7.1"
 Reexport = "0.2, 1"
 SciMLBase = "2, 3"
 SpecialFunctions = "0.7.2, 0.8, 0.9, 0.10, 1, 2"

src/mcmc/Inference.jl

@@ -32,7 +32,6 @@ import AdvancedHMC
 const AHMC = AdvancedHMC
 import AdvancedMH
 const AMH = AdvancedMH
-import AdvancedPS
 import EllipticalSliceSampling
 import LogDensityProblems
 import Random