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,31 @@
+# 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 only; `PG` / `CSMC` draw a conditional sweep's ancestors from the categorical over the weights.
+Exact draws may therefore differ from previous releases, but remain statistically consistent (the same target distribution).
+
+Chain statistics have changed: `chain[:logevidence]` is now `chain[:log_normalizing_constant]`, and SMC's per-particle `weight` is gone, since the returned particles are now equal-weight.
+
+`SMC` no longer has a sampler state, because it runs one sweep rather than an MCMC loop: `save_state` and `initial_state` are now warned about and ignored, so `loadstate` has nothing to return for an `SMC` chain. `PG` / `CSMC` are unaffected.
+
+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. Previously, results could drift between Julia versions even under a `StableRNG` (https://github.com/TuringLang/Turing.jl/issues/2781).
+  - Parallelism within a sweep. `SMC(; multithreaded=true)` / `PG(n; multithreaded=true)` spread that sweep's particles across threads without changing the results; start Julia with multiple threads (e.g. `julia -t auto`) for this to take effect. It is independent of `MCMCThreads()` / `MCMCDistributed()`, which parallelise whole chains and work with SMC/PG as with any other sampler.
+  - Equal-weight draws. `SMC` resamples once at the end of the sweep, so `mean(chain[...])` and other summaries need no weighting.
+  - A degeneracy diagnostic. `SMC` chains carry `ess_per_step`, the effective sample size after each filtering step; one entry per likelihood term, so an `@addlogprob!` adds one alongside the observations. MCMCChains exposes the entries as `ess_per_step[1]`, `ess_per_step[2]`, and so on.
+  - 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` it 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"

ext/TuringMCMCChainsExt.jl

@@ -2,10 +2,10 @@ module TuringMCMCChainsExt
 
 using Turing
 using Turing: AbstractMCMC, DynamicPPL
-using Turing.Inference: HMC, NUTS, HMCDA, Emcee, EmceeState, _get_n_walkers
+using Turing.Inference: HMC, NUTS, HMCDA, Emcee, EmceeState, SMC, _get_n_walkers
 using MCMCChains: MCMCChains
 
-import Turing.Inference: post_sample_hook
+import Turing.Inference: bundle_smc_samples, post_sample_hook
 
 """
     loadstate(chain::MCMCChains.Chains)
@@ -44,6 +44,28 @@ function AbstractMCMC.bundle_samples(
     return AbstractMCMC.chainscat(chains...)
 end
 
+function bundle_smc_samples(
+    transitions::Vector{<:DynamicPPL.ParamsWithStats},
+    model::DynamicPPL.Model,
+    sampler::SMC,
+    state,
+    ::Type{MCMCChains.Chains};
+    kwargs...,
+)
+    # MCMCChains drops non-scalar statistics, so give each filtering step its own internal.
+    ess_per_step = first(transitions).stats.ess_per_step
+    ess_names = ntuple(i -> Symbol("ess_per_step[$i]"), length(ess_per_step))
+    ess_stats = NamedTuple{ess_names}(Tuple(ess_per_step))
+    scalar_transitions = map(transitions) do transition
+        names = filter(!=(:ess_per_step), keys(transition.stats))
+        stats = merge(NamedTuple{names}(transition.stats), ess_stats)
+        DynamicPPL.ParamsWithStats(transition.params, stats)
+    end
+    return AbstractMCMC.bundle_samples(
+        scalar_transitions, model, sampler, state, MCMCChains.Chains; kwargs...
+    )
+end
+
 """
     post_sample_hook(chain::MCMCChains.Chains, sampler::Union{HMC,NUTS,HMCDA}; kwargs...)
 

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