Drop Turing/DynamicPPL support

CHANGELOG.md

@@ -22,8 +22,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `ADTypes.AutoReactant()` can supply derivatives and sampler HVPs through
   Reactant.jl. It requires `using Reactant`, a traceable log-density, and
   `AutoReactant()` on both AD passes of an HVP (#37, #52).
-- The `DensityModel` constructors in `DynamicPPLExt` and `LogDensityProblemsExt`
-  now accept `hvp` and the batched derivative slots.
+- The `DensityModel` constructor in `LogDensityProblemsExt` now accepts `hvp`
+  and the batched derivative slots.
 
 ### Fixed
 
@@ -54,7 +54,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Removed
 
-- `DynamicPPLExt` no longer requires `ForwardDiff` as a triggering library to load.
+- Turing/DynamicPPL support (#74). `DynamicPPLExt` and the
+  `DensityModel(::DynamicPPL.Model)` constructor are gone, since DynamicPPL
+  models cannot run on the GPU. Use ParallelMCMC 0.2 for Turing models, or wrap
+  the log-density in a `LogDensityProblems` object or a plain `DensityModel`.
 
 ## [0.2.0] - 2026-06-29
 

Project.toml

@@ -15,7 +15,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-DynamicPPL = "366bfd00-2699-11ea-058f-f148b4cae6d8"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
@@ -25,7 +24,6 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
 CUDAExt = "CUDA"
-DynamicPPLExt = ["DynamicPPL", "FlexiChains", "LogDensityProblems"]
 EnzymeExt = "Enzyme"
 LogDensityProblemsExt = "LogDensityProblems"
 ReactantExt = ["Reactant", "Enzyme"]
@@ -36,7 +34,6 @@ AbstractMCMC = "5.10.0"
 CUDA = "5.11.0"
 CUDA_Runtime_jll = "0.21"
 DifferentiationInterface = "0.7.13"
-DynamicPPL = "0.41.6, 0.42"
 Enzyme = "0.13.146"
 FlexiChains = "0.6.6"
 ForwardDiff = "1"

README.md

@@ -41,7 +41,7 @@ The approach and its scaling tricks (stochastic Hutchinson Jacobian estimators,
 | [`MALASampler`](src/interface.jl) | Baseline — sequential MALA with a fixed step size |
 | [`AdaptiveMALASampler`](src/interface.jl) | Baseline — sequential MALA with dual-averaging step-size adaptation |
 
-All samplers implement the [AbstractMCMC](https://github.com/TuringLang/AbstractMCMC.jl) interface and return [`FlexiChains`](https://pysm.dev/FlexiChains.jl) objects, so they slot into existing Turing.jl / AbstractMCMC workflows.
+All samplers implement the [AbstractMCMC](https://github.com/TuringLang/AbstractMCMC.jl) interface and return [`FlexiChains`](https://pysm.dev/FlexiChains.jl) objects, so they slot into existing AbstractMCMC workflows.
 
 
 ### Quick start
@@ -66,7 +66,7 @@ sampler = ParallelMALASampler(0.1; T=64, jacobian=:stoch_diag,
 chain = sample(model, sampler, 500; chain_type=VNChain)
 ```
 
-See the [Getting Started guide](docs/src/10-getting-started.md) for worked examples, GPU usage, Turing.jl integration, and step-size tuning.
+See the [Getting Started guide](docs/src/10-getting-started.md) for worked examples, GPU usage, and step-size tuning.
 
 ## How to Cite
 

docs/src/10-getting-started.md

@@ -47,7 +47,7 @@ Hand-written and AD-derived slots can be mixed.
 | backend | backend | `SecondOrder(hvp, grad_logdensity)` on `logdensity` |
 | either | `SecondOrder(...)` | that pair on `logdensity`, gradient slot unused |
 
-An explicit `SecondOrder` bypasses the gradient slot. Use it for an AD-derived HVP on Turing or LogDensityProblems models.
+An explicit `SecondOrder` bypasses the gradient slot. Use it when the gradient can't be differentiated.
 
 See [GPU Execution](15-gpu.md) for backend compatibility.
 
@@ -103,59 +103,37 @@ The prefix-scan kernel runs as pure array broadcasts and is array-type-agnostic.
 2. Pass for example, `backend = ADTypes.AutoEnzyme()` to `ParallelMALASampler`.
 3. Pass a `CuVector` as `initial_params` to `sample`.
 
-GPU execution has its own caveats (Turing models are CPU-only, `AutoEnzyme()` currently needs the `pmcmc_matmul` / `pmcmc_dot` / `pmcmc_dotsum` wrappers, gradients should avoid scalar indexing).  See [GPU Execution](15-gpu.md) for the worked example and the full set of limitations.
+GPU execution has its own caveats (`AutoEnzyme()` currently needs the `pmcmc_matmul` / `pmcmc_dot` / `pmcmc_dotsum` wrappers, gradients should avoid scalar indexing).  See [GPU Execution](15-gpu.md) for the worked example and the full set of limitations.
 
 ---
 
-## Turing.jl integration
+## LogDensityProblems models
 
-ParallelMCMC.jl integrates with Turing.jl models through the `DynamicPPL` and `LogDensityProblems` extensions.
-
-### One-step convenience constructor
-
-Load `DynamicPPL` (part of Turing.jl) and a single-argument `DensityModel` constructor becomes available:
+With `LogDensityProblems` loaded, `DensityModel` accepts any object that implements the [LogDensityProblems](https://github.com/tpapp/LogDensityProblems.jl) interface with gradients:
 
 ```julia
-using Turing, ParallelMCMC, FlexiChains
+using LogDensityProblems, ParallelMCMC, FlexiChains
 
-@model function normal_model(y)
-    μ ~ Normal(0.0, 1.0)
-    y ~ Normal(μ, 0.5)
+# μ ~ Normal(0, 1), y ~ Normal(μ, 0.5)
+struct NormalPosterior
+    y::Float64
+end
+LogDensityProblems.capabilities(::Type{NormalPosterior}) =
+    LogDensityProblems.LogDensityOrder{1}()
+LogDensityProblems.dimension(::NormalPosterior) = 1
+LogDensityProblems.logdensity(ld::NormalPosterior, x) = -0.5 * x[1]^2 - 2 * (ld.y - x[1])^2
+function LogDensityProblems.logdensity_and_gradient(ld::NormalPosterior, x)
+    return LogDensityProblems.logdensity(ld, x), [-x[1] + 4 * (ld.y - x[1])]
 end
 
-model = DensityModel(normal_model(1.5))   # param_names=[:μ] extracted automatically
-
-chain = sample(model, ParallelMALASampler(0.1; T=64, backend=AutoEnzyme()), 500;
-               chain_type=VNChain)
+model = DensityModel(NormalPosterior(1.5); param_names=[:μ])
+chain = sample(model, AdaptiveMALASampler(0.3; n_warmup=500), 2_000;
+               chain_type=VNChain, discard_warmup=true)
 ```
 
-Much like Turing's own samplers, the resulting chain will always have parameters in the original (possibly constrained) space, even though the MCMC sampling itself is performed in unconstrained space.
-Furthermore, parameter names are automatically extracted from the Turing model (and will always be the same as those when using Turing's own samplers).
-
-Note that when sampling with a Turing model the returned chain will have `:logjoint`, `:logprior`, and `:loglikelihood` columns, since Turing models provide enough information to separate these contributions to the log-density.
-This is in contrast to sampling with a manually constructed `DensityModel`, which returns only a single `:logp` column.
+The object must implement `capabilities` (at least `LogDensityOrder{1}`), `dimension`, `logdensity`, and `logdensity_and_gradient`. `hvp` and the batched slots are forwarded to `DensityModel`.
 
-### Manual `LogDensityProblems` path
-
-For explicit control over the AD backend, construct the `LogDensityFunction`
-directly with DynamicPPL's `adtype` interface:
-
-```julia
-using Turing, LogDensityProblems, ADTypes
-using ParallelMCMC, FlexiChains
-
-ld = DynamicPPL.LogDensityFunction(
-    normal_model(1.5),
-    DynamicPPL.getlogjoint_internal,
-    DynamicPPL.LinkAll();
-    adtype=ADTypes.AutoEnzyme(),
-)
-
-model = DensityModel(ld; param_names=[:μ])
-```
-
-This also accepts any other `LogDensityProblems`-compatible object.
-As above, the returned chain will always contain parameters in the original space: the use of `LinkAll()` above only stipulates that MCMC sampling itself is to be performed in unconstrained space, and has no result on the form of the returned chain.
+A backend in `hvp` differentiates `logdensity_and_gradient`, so that function must be differentiable by the backend. If it is not (for example, it calls AD internally), pass an explicit `SecondOrder`, which differentiates `logdensity` twice instead. `AutoReactant()` is not supported for these models.
 
 ---
 
@@ -207,8 +185,6 @@ model = DensityModel(...; param_names=[:x1, :x2, :x3])
 model = DensityModel(...; param_names=[:x, (:y, (1, 2))])
 ```
 
-For Turing.jl models, parameter names are automatically derived from the model and do not need to be specified manually.
-
 ---
 
 ## Sequential MALA baselines

docs/src/15-gpu.md

@@ -22,13 +22,7 @@ A useful heuristic: if your CPU `grad_logdensity_batch` already saturates BLAS a
 
 ## Limitations
 
-### 1. DynamicPPL / Turing models are CPU-only
-
-The one-argument `DensityModel(::DynamicPPL.Model)` constructor (and the `LogDensityProblems` path that backs it) executes the model with `Vector{Float64}` parameters on CPU.  There is no path today that runs a `@model` against `CuArray` parameters.
-
-To run on GPU, hand-write `logdensity` / `grad_logdensity` (and the batched variants) directly over `CuArray`s and pass them to `DensityModel`.  This is shown in the example below.
-
-### 2. Enzyme on GPU currently needs `pmcmc_matmul` / `pmcmc_dot` / `pmcmc_dotsum`
+### 1. Enzyme on GPU currently needs `pmcmc_matmul` / `pmcmc_dot` / `pmcmc_dotsum`
 
 `ParallelMCMC` exports three thin wrappers:
 
@@ -58,9 +52,9 @@ When you do **not** need them:
 
 These wrappers are a workaround pending upstream Enzyme support for GPU matmul rules.  Once that lands, plain `A * B` on `CuArray`s will work and the wrappers will be deprecated.
 
-### 3. Some Enzyme GPU gradients hit a `gc-transition` abort. Staging them as single ops is a workaround
+### 2. Some Enzyme GPU gradients hit a `gc-transition` abort. Staging them as single ops is a workaround
 
-Certain `CuArray` gradient expressions abort during Enzyme compilation on GPU with the same `gc-transition` failure as limitation 2. In practice the fused form below aborts while the staged form compiles, so **writing the gradient as single-op stages is a reliable workaround**:
+Certain `CuArray` gradient expressions abort during Enzyme compilation on GPU with the same `gc-transition` failure as limitation 1. In practice the fused form below aborts while the staged form compiles, so **writing the gradient as single-op stages is a reliable workaround**:
 
 ```julia
 # ABORTS during Enzyme compile on GPU
@@ -253,7 +247,7 @@ sampler = ParallelMALASampler(0.005f0; T=16, backend=AutoReactant())
 
     **Reactant chooses the execution device independently of the input array.** With a CUDA client, `CuArray` inputs are staged through pinned host memory and results are copied back on the device. With the CPU client, both directions go through the host and the sampler warns when it prepares the model. Select a GPU client with `Reactant.set_default_backend` when available.
 
-- The log-density must be Reactant-traceable. DynamicPPL-built log-densities are not.
+- The log-density must be Reactant-traceable.
 - Across an HVP's two AD passes, use `AutoReactant()` for both or neither. A hand-written gradient may pair with it. `SecondOrder` cannot contain `AutoReactant()`.
 - Only the default `AutoReactant()` mode is supported.
 
@@ -268,7 +262,7 @@ sampler = ParallelMALASampler(0.005f0; T=16, backend=AutoReactant())
 
 - **The HVP has a clean closed form.**  Quadratic priors, Gaussian likelihoods, GLMs (logistic, Poisson, probit) — the second derivative is a known function of intermediate quantities you already compute in `gradlogp`.  A few extra lines and you skip the AD pipeline entirely.
 - **Performance matters and the AD compile is heavy.**  Enzyme and Mooncake both pay a one-shot compilation cost on the user's gradient.  For long-running chains this amortizes, but for many short runs the analytical HVP wins.
-- **You're hitting AD-backend-specific GPU restrictions.**  The [Enzyme limitations](#2-enzyme-on-gpu-currently-needs-pmcmc_matmul-pmcmc_dot-pmcmc_dotsum) above (`pmcmc_*` wrappers, staged broadcasts) only matter when the AD backend is invoked.  Supplying analytical HVP sidesteps them — your `gradlogp` and `hvp` can use plain `*`, `dot`, `sum`, and the sampler's `backend` can be omitted entirely because no AD is ever invoked.
+- **You're hitting AD-backend-specific GPU restrictions.**  The [Enzyme limitations](#1-enzyme-on-gpu-currently-needs-pmcmc_matmul-pmcmc_dot-pmcmc_dotsum) above (`pmcmc_*` wrappers, staged broadcasts) only matter when the AD backend is invoked.  Supplying analytical HVP sidesteps them — your `gradlogp` and `hvp` can use plain `*`, `dot`, `sum`, and the sampler's `backend` can be omitted entirely because no AD is ever invoked.
 - **You can reuse intermediates between gradient and HVP.**  When `hvp` shares $X\beta$, $\sigma(X\beta)$, or similar with the gradient computation, an analytical version can be both faster *and* shorter than what AD produces.
 
 ### Same example with analytical HVP

docs/src/95-reference.md

@@ -20,14 +20,9 @@ derivative-slot backend on `DensityModel` and as `ParallelMALASampler`'s
 ParallelMCMC.needs_host_staging
 ```
 
-## Extension constructors
+## Extension constructor
 
-`DensityModel` also has extension constructors for common probabilistic-programming interfaces:
-
-- `DensityModel(ld; param_names=nothing)` for `LogDensityProblems` models with gradients
-- `DensityModel(turing_model)` for `DynamicPPL` / Turing models when the relevant extension packages are loaded
-
-See [Getting Started](10-getting-started.md) for end-to-end examples of both.
+`DensityModel(ld; param_names=nothing)` wraps any `LogDensityProblems` object with gradients when `LogDensityProblems` is loaded. See [LogDensityProblems models](@ref) for an example.
 
 ## Model
 

docs/src/index.md

@@ -95,31 +95,8 @@ chain = sample(model, sampler, 2_000;
                progress=true)
 ```
 
-### Turing.jl integration
-
-When `DynamicPPL` (part of Turing.jl) is loaded, a one-argument `DensityModel` constructor is available that wraps a `@model` directly.
-Parameter names are automatically extracted, and values transformed back to the original model space:
-
-```julia
-using Turing, ParallelMCMC, FlexiChains
-
-@model function normal_model(y)
-    μ ~ Normal(0.0, 1.0)
-    y ~ Normal(μ, 0.5)
-end
-
-model   = DensityModel(normal_model(1.5))
-sampler = AdaptiveMALASampler(0.3; n_warmup=500)
-
-chain = sample(model, sampler, 2_000;
-               chain_type=VNChain,
-               discard_warmup=true)
-```
-
 See [Getting Started](10-getting-started.md) for worked examples and guidance on choosing samplers, and [Algorithm Details](20-algorithms.md) for the mathematics behind DEER.
 
-For Turing models the chain type used must be `VNChain` (not `SymChain`), as that is the natural parameter type for Turing models.
-
 ## Contributors
 
 ```@raw html