PG silently gives incorrect results with `Threads.@threads`

[penelopeysm]

Minimal working example

using Turing, Random
@model function h(y)
    x ~ MvNormal(zeros(length(y)), I)
    Threads.@threads for i in eachindex(y)
        y[i] ~ Normal(x[i])
    end
end
chn = sample(Xoshiro(468), h([1.0, 2.0, 3.0]), PG(10), MCMCThreads(), 2000, 4)
describe(chn)

Description

x should roughly be [0.5, 1.0, 1.5]. That is indeed what you get with any other sampler, but PG gives this:

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ess_per_sec
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64       Float64

        x[1]    0.0010    0.9965    0.0123   6530.6140   6203.9280    1.0004       23.5784
        x[2]   -0.0169    1.0029    0.0125   6395.6826   5863.2289    1.0002       23.0912
        x[3]    0.0109    1.0102    0.0125   6480.7591   6038.9386    1.0000       23.3984

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

        x[1]   -1.9788   -0.6553   -0.0075    0.6790    1.9822
        x[2]   -1.9955   -0.6864   -0.0164    0.6486    1.9556
        x[3]   -1.9920   -0.6598    0.0108    0.6944    1.9834

It is known that SMC / PG don't work with multithreaded models. The reason for this is obvious, we don't know which observation will happen first, so in general the weighting and resampling will be all over the place.

Turing.jl/src/mcmc/particle_mcmc.jl

Lines 386 to 390 in ac93c30

	function DynamicPPL.use_threadsafe_eval(
	::DynamicPPL.SamplingContext{<:Sampler{<:Union{PG,SMC}}}, ::AbstractVarInfo
	)
	return false
	end

However, this is not respected. Or, to be precise, setting use_threadsafe_eval = false does not actually forbid multithreaded models from being run: it merely prevents DynamicPPL from using ThreadSafeVarInfo. So the model will still run but with ordinary VarInfo, which is arguably even worse.

Julia version info

versioninfo()

Julia Version 1.11.6
Commit 9615af0f269 (2025-07-09 12:58 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 10 × Apple M1 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 4 default, 0 interactive, 2 GC (on 8 virtual cores)

Manifest

Turing v0.40.1

[penelopeysm]

TuringLang/DynamicPPL.jl#1151 and #2715 will fix this

[penelopeysm]

Fixed by the above ^

Now all samplers will give incorrect results with Threads.@threads if threadsafe evaluation is not enabled. If threadsafe evaluation is enabled, PG will error instead of silently giving a wrong result.