@@ -57,6 +57,7 @@ include("kernel.jl")
 include("implicit.jl")
 include("energy_network.jl")
 include("answers.jl")
+include("product.jl")
 include("factor.jl")
 include("implicit_factor.jl")
 include("proxdm.jl")
@@ -88,6 +89,7 @@ export ImplicitProx, implicit_solution, implicit_factor_pullback
 # --- Proximal diffusion models (ProxDM) ----------------------------------------
 export EnergyNetwork, energy, implicit_energy, denoising_gradient
 export implicit_roots, implicit_laplace, density_lambda
+export AbstractImplicitRelation, ProductRelation
 export AbstractProximalPredictor, proximal, ProxNetwork, MixtureProx
 export proxdm_sample, prox_infer, proximal_matching_loss
 
 
@@ -7,7 +7,7 @@
 >
 > Theory (CT-ML wiki): [Bayesian Inversion](https://mathstruct.org/CategoryTheory-ML-Wiki/Bayesian-Inversion) · [Statistical Game](https://mathstruct.org/CategoryTheory-ML-Wiki/Statistical-Game) · [Lens](https://mathstruct.org/CategoryTheory-ML-Wiki/Lens)
 
-## The chain, in ten files
+## The chain, in eleven files
 
 | file | note | supplies |
 |---|---|---|
@@ -19,6 +19,7 @@
 | `kernel.jl` | [[kernel]] | kernel density estimates as relations: batch KDE with bandwidth selection, an online forgetting KDE |
 | `implicit.jl` | [[implicit]] | deterministic implicit inference and its adjoint backward pass |
 | `energy_network.jl` | [[energy_network]] | energy-parametrised predictors: $\varepsilon = \sigma_t\nabla_x E_\theta$, a conservative score, an energy for the relation |
+| `product.jl` | [[product]] | products of relations on the same variables: intersections, with the adjoint into every factor |
 | `implicit_factor.jl` | [[implicit_factor]] | `DiffusionFactor` with `ImplicitProx`: report and per-channel pullback |
 | `proxdm.jl` | [[proxdm]] | proximal diffusion models: prox interface, exact oracle, sampler, proximal inference |
 | `ext/` | [[backends]] | the network's derivatives through any AD backend (DifferentiationInterface, Reactant) |
 
@@ -27,7 +27,7 @@ predictor is an [`EnergyNetwork`](@ref) (lowest first, i.e. the most plausible a
 and otherwise in the order found. There is no guarantee that every stable answer is found;
 more starts and a larger `spread` find more.
 """
-function implicit_roots(m::ImplicitDiffusion, z₀, ρ, ps, st; nstarts::Integer = 16, spread = 1.0,
+function implicit_roots(m::AbstractImplicitRelation, z₀, ρ, ps, st; nstarts::Integer = 16, spread = 1.0,
                         starts = nothing, rng::AbstractRNG = Random.Xoshiro(0), unique_tol = 1e-3, kwargs...)
     free = .!isinf.(ρ)
     inits = if starts === nothing
@@ -44,13 +44,16 @@ function implicit_roots(m::ImplicitDiffusion, z₀, ρ, ps, st; nstarts::Integer
         push!(roots, sol)
     end
     roots = [r for r in roots]                                   # concrete element type
-    if m.predictor isa NoisePredictor{<:EnergyNetwork} && length(roots) > 1
+    if _has_energy(m) && length(roots) > 1
         energies = [_query_energy(m, r.z, z₀, ρ, ps, st) for r in roots]
         roots = roots[sortperm(energies)]
     end
     return (roots, st)
 end
 
+_has_energy(m) = false
+_has_energy(::ImplicitDiffusion{<:NoisePredictor{<:EnergyNetwork}}) = true
+
 function _query_energy(m, z, z₀, ρ, ps, st)
     U, _ = implicit_energy(m, z, ps, st)
     soft = .!isinf.(ρ)
@@ -93,7 +96,7 @@ covariance is then calibrated to the density *smoothed at the field's noise leve
 the data distribution by those levels). With another λ it is a relative measure of how sharply
 the answer is determined. Throws unless `sol` converged to a stable answer.
 """
-function implicit_laplace(m::ImplicitDiffusion, sol::ImplicitSolution, ρ, ps, st)
+function implicit_laplace(m::AbstractImplicitRelation, sol::ImplicitSolution, ρ, ps, st)
     (sol.converged && sol.stable) ||
         throw(ArgumentError("the Laplace approximation needs a converged, stable answer (a local minimum)"))
     free = .!isinf.(ρ)
 
@@ -74,13 +74,24 @@ is the denoiser (`Deterministic Relaxation.md`).
 """
 noisefree_nodes(n::Integer, t::Real) = FieldNodes([float(t)], zeros(n, 1), [1.0])
 
+"""
+    AbstractImplicitRelation
+
+A learned relation that can be queried by root-finding: anything with a field
+[`prior_field`](@ref) and its Jacobian [`prior_jacobian`](@ref). [`ImplicitDiffusion`](@ref) is one
+diffusion model; [`ProductRelation`](@ref) pools several. The query functions
+([`implicit_infer`](@ref), [`implicit_roots`](@ref), [`implicit_laplace`](@ref),
+[`implicit_pullback`](@ref)) accept any of them.
+"""
+abstract type AbstractImplicitRelation end
+
 """
     ImplicitDiffusion(predictor, nodes; λ = 1.0)
 
 The relation ``R_\\theta = \\{z : r(z) = 0\\}`` defined by a noise predictor and a fixed node
 set, with RED-Diff's weighting ``\\lambda_t = \\lambda\\sigma_t/\\alpha_t``.
 """
-struct ImplicitDiffusion{P<:NoisePredictor,T<:Real}
+struct ImplicitDiffusion{P<:NoisePredictor,T<:Real} <: AbstractImplicitRelation
     predictor::P
     nodes::FieldNodes{T}
     λ::T
@@ -130,7 +141,7 @@ end
 ``r = g_\\theta(z) + \\rho^2\\odot(z - z_0)`` on the free coordinates (finite `ρ`), zero on the
 hard-clamped ones (``\\rho_i = \\infty``), where the constraint is ``z_i = z_{0,i}`` instead.
 """
-function implicit_residual(m::ImplicitDiffusion, z, z₀, ρ, ps, st)
+function implicit_residual(m::AbstractImplicitRelation, z, z₀, ρ, ps, st)
     g, st = prior_field(m, z, ps, st)
     free = .!isinf.(ρ)
     r = zero(g)
@@ -177,7 +188,7 @@ is only trusted in locally convex regions so that it is not drawn to saddles and
 `z₀` holds the clamp values on hard coordinates and the soft targets elsewhere; `ρ` is the
 precision vector, `Inf` for a hard clamp (cf. [`precision_vector`](@ref)).
 """
-function implicit_infer(m::ImplicitDiffusion, z₀, ρ, ps, st;
+function implicit_infer(m::AbstractImplicitRelation, z₀, ρ, ps, st;
                         z_init = z₀, tol = 1e-9, maxiters::Integer = 200, step = 0.05)
     free = .!isinf.(ρ)
     z = float.(copy(z_init))
@@ -244,7 +255,7 @@ the Jacobians are evaluated at ``z^\\star`` afresh. Throws if `sol` did not conv
 implicit function theorem says nothing about a point that is not a root
 (`Backpropagation through Implicit Inference.md` §6).
 """
-function implicit_pullback(m::ImplicitDiffusion, sol::ImplicitSolution, z₀, ρ, z̄, ps, st)
+function implicit_pullback(m::AbstractImplicitRelation, sol::ImplicitSolution, z₀, ρ, z̄, ps, st)
     sol.converged || throw(ArgumentError(
         "implicit_pullback at a non-converged state (residual $(sol.residual)): the implicit " *
         "function theorem does not apply. Re-solve, or use an unrolled/one-step surrogate."))
@@ -256,12 +267,7 @@ function implicit_pullback(m::ImplicitDiffusion, sol::ImplicitSolution, z₀, ρ
     λF = -(Jff' \ z̄[free])
     λ = zero(float.(z)); λ[free] .= λF
 
-    ps̄ = _zero(ps)
-    for k in (_hasparams(ps) ? eachindex(m.nodes.t) : ())
-        t = m.nodes.t[k]
-        g = epsilon_vjp_params(m.predictor, _node_input(m, z, k), t, ps, st, λ)
-        ps̄ = _axpy(m.nodes.w[k] * _weight(m, t), g, ps̄)
-    end
+    ps̄ = _params_cotangent(m, z, λ, ps, st)
 
     z̄₀ = zero(float.(z))
     z̄₀[free] .= .-(ρ[free] .^ 2) .* λF
@@ -271,6 +277,17 @@ function implicit_pullback(m::ImplicitDiffusion, sol::ImplicitSolution, z₀, ρ
     return (z₀ = z̄₀, ρ = ρ̄, ps = ps̄)
 end
 
+# (∂_θ g)ᵀ λ: the parameter part of the adjoint, one VJP per node
+function _params_cotangent(m::ImplicitDiffusion, z, λ, ps, st)
+    ps̄ = _zero(ps)
+    for k in (_hasparams(ps) ? eachindex(m.nodes.t) : ())
+        t = m.nodes.t[k]
+        g = epsilon_vjp_params(m.predictor, _node_input(m, z, k), t, ps, st, λ)
+        ps̄ = _axpy(m.nodes.w[k] * _weight(m, t), g, ps̄)
+    end
+    return ps̄
+end
+
 # --- As a factor's inversion ----------------------------------------------------
 
 """
 
@@ -0,0 +1,75 @@
+# ---------------------------------------------------------------------------
+# Several learned relations on the same variables, pooled: the product of their densities.
+#
+# Each implicit learner's field g_i is (for an ideal model) the gradient of a smoothed negative
+# log-density, so the weighted sum Σ β_i g_i is the field of the product Π p_i^{β_i}: its stable
+# roots are where all relations hold at once, e.g. the intersection points of a learned circle
+# and a learned line. Every query function works on the product unchanged.
+#
+# See `product.md` and `Composing Diffusion Factors` in the vault.
+# ---------------------------------------------------------------------------
+
+"""
+    ProductRelation(relations...; weights = ones)
+
+Several implicit relations on the **same** joint space, pooled: the field is the weighted sum
+``\\sum_i \\beta_i\\,g_i`` of the relations' fields, i.e. the field of the product of their densities
+``\\prod_i p_i^{\\beta_i}``. Its stable roots are the configurations that satisfy all relations
+together, e.g. the intersection points of a learned circle and a learned line. The weights
+``\\beta_i`` set how much each relation is trusted (a tempering exponent).
+
+Parameters and states are **tuples**, one entry per relation, in order: `ps = (ps₁, ps₂)`. Every
+query function accepts a product: [`implicit_infer`](@ref), [`implicit_roots`](@ref),
+[`implicit_laplace`](@ref), and [`implicit_pullback`](@ref) (which returns a tuple of parameter
+cotangents). Products may be nested.
+
+Pooling smoothed densities is exact only at zero noise: at the field's noise levels the
+product of the smoothed densities is not the smoothed product, so intersections are found
+with the same kind of smoothing bias as single relations (`product.md` §3).
+"""
+struct ProductRelation{M<:Tuple,W<:AbstractVector} <: AbstractImplicitRelation
+    relations::M
+    weights::W
+end
+function ProductRelation(relations::AbstractImplicitRelation...; weights = ones(length(relations)))
+    isempty(relations) && throw(ArgumentError("a product needs at least one relation"))
+    length(weights) == length(relations) || throw(ArgumentError("need one weight per relation"))
+    all(>(0), weights) || throw(ArgumentError("weights must be positive"))
+    return ProductRelation(relations, float.(collect(weights)))
+end
+
+_check_tuple(m::ProductRelation, x, what) = (x isa Tuple && length(x) == length(m.relations)) ||
+    throw(ArgumentError("a ProductRelation of $(length(m.relations)) relations needs `$what` as a tuple with one entry per relation"))
+
+function prior_field(m::ProductRelation, z, ps, st)
+    _check_tuple(m, ps, "ps"); _check_tuple(m, st, "st")
+    g = zero(float.(z))
+    sts = map(eachindex(m.relations)) do i
+        gi, sti = prior_field(m.relations[i], z, ps[i], st[i])
+        g = g .+ m.weights[i] .* gi
+        sti
+    end
+    return (g, Tuple(sts))
+end
+
+function prior_jacobian(m::ProductRelation, z, ps, st)
+    _check_tuple(m, ps, "ps"); _check_tuple(m, st, "st")
+    return sum(m.weights[i] .* prior_jacobian(m.relations[i], z, ps[i], st[i]) for i in eachindex(m.relations))
+end
+
+_params_cotangent(m::ProductRelation, z, λ, ps, st) =
+    Tuple(_params_cotangent(m.relations[i], z, m.weights[i] .* λ, ps[i], st[i]) for i in eachindex(m.relations))
+
+# a product of energy-parametrised relations has the weighted sum of their energies
+_has_energy(m::ProductRelation) = all(_has_energy, m.relations)
+
+function implicit_energy(m::ProductRelation, z, ps, st)
+    _has_energy(m) || throw(ArgumentError("implicit_energy needs every relation of the product to be energy-parametrised"))
+    U = zero(float(eltype(z)))
+    sts = map(eachindex(m.relations)) do i
+        Ui, sti = implicit_energy(m.relations[i], z, ps[i], st[i])
+        U += m.weights[i] * Ui
+        sti
+    end
+    return (U, Tuple(sts))
+end
@@ -0,0 +1,88 @@
+#implementation
+
+> Products of relations. Several implicit learners on the **same** variables, pooled: the field
+> of the product is the weighted sum of their fields, so its stable roots are where all the
+> relations hold at once, e.g. the two points where a learned circle meets a learned line.
+> `ProductRelation` is one more `AbstractImplicitRelation`, so every query function works on it
+> unchanged, and gradients flow back into each factor.
+
+> Sources: code: `product.jl`, `implicit.jl` (`AbstractImplicitRelation`, `_params_cotangent`), `answers.jl`; Hinton, *Training Products of Experts by Minimizing Contrastive Divergence*, Neural Computation 2002; Du et al., *Reduce, Reuse, Recycle: Compositional Generation with Energy-Based Diffusion Models and MCMC*, ICML 2023 (adding scores, and why it is not exact at $t > 0$)
+>
+> Bibliography: [[Bibliography#^hinton2002poe|Hinton 2002]] · [[Bibliography#^du2023reduce|Du et al. 2023]]
+>
+> Theory: [[Composing Diffusion Factors]] · [[Implicit Diffusion Learners]] · [[Open Problems in Implicit Diffusion Learning]]
+
+## 1. The construction
+
+An implicit learner's field $g(z) = \sum_k w_k\lambda_k\,(\varepsilon_\theta(\alpha_k z + \sigma_k\varepsilon_k, t_k) - \varepsilon_k)$
+is, for an ideal predictor, the gradient of a smoothed negative log-density. A weighted sum of
+such fields is therefore the field of the product of the densities:
+
+$$g_\Pi(z) = \sum_i \beta_i\, g_i(z) \quad\longleftrightarrow\quad p_\Pi \propto \prod_i p_i^{\beta_i}.$$
+
+`ProductRelation(m₁, m₂, …; weights = β)` implements exactly the three functions the solvers
+use:
+
+| function | product |
+|---|---|
+| `prior_field` | $\sum_i \beta_i g_i$ |
+| `prior_jacobian` | $\sum_i \beta_i J_i$ |
+| `_params_cotangent` | a tuple: factor $i$ gets the adjoint $\beta_i\lambda$ |
+
+Parameters and states are tuples, one entry per factor: `ps = (ps₁, ps₂)`. Passing a single
+parameter set raises an `ArgumentError`. Products nest, because a product is itself an
+`AbstractImplicitRelation`.
+
+## 2. What works on it unchanged
+
+- `implicit_infer`: one answer, with clamps (`ρ = Inf`) and soft anchors as usual;
+- `implicit_roots`: all answers by multi-start, e.g. the intersection points;
+- `implicit_laplace`: a covariance per answer;
+- `implicit_pullback`: the IFT adjoint. The cotangent of the parameters is a tuple, so a loss
+  on the intersection trains every factor;
+- `implicit_energy`, and ordering roots by energy, when every factor is an `EnergyNetwork`
+  (the energy of the product is $\sum_i\beta_i U_i$).
+
+The refactor that made this possible is small. `implicit.jl` now types its functions with
+the abstract supertype, and moves the parameter loop of the pullback into `_params_cotangent`,
+the one place where a product differs.
+
+## 3. Accuracy: the smoothing bias
+
+Pooling is exact only at zero noise. At the field's noise levels the product of the smoothed
+densities is not the smoothed product, $q_t * (p_1p_2) \ne (q_t * p_1)(q_t * p_2)$
+([[Open Problems in Implicit Diffusion Learning]] T9). For root finding this produces the same kind of
+bias that a single relation already has (its ridge sits slightly on the inside of a bend), not a
+failure. Closed-form mixtures (`GaussianMixtureEps`, $s = 0.05$), the default `field_nodes`:
+
+| product | answers found | exact | largest distance to exact |
+|---|---|---|---|
+| circle ∩ line $y = x$ | 2 | 2 | 0.022 |
+| circle ∩ ellipse (1.5, 0.6) | 4 | 4 | 0.029 |
+
+The adjoint is exact for the field that is actually solved. Clamping $x = 0.95$ on circle ∩
+ellipse, $\partial y/\partial x$ is −0.65669290204 by the adjoint and −0.65669290206 by finite
+differences. A parameter of the ellipse gives 0.0220369733 by both.
+
+## 4. Limits
+
+- **Same variables only.** All factors see the whole joint space. Factors on *different* but
+  overlapping variable sets (a product with a separator, $p_Ap_B/p_S$) are not covered.
+- **Tangent intersections** are ill-conditioned: the Jacobian of the sum is nearly singular,
+  so the Laplace covariance blows up along the common tangent. This is the right answer, but
+  slow to converge.
+- **Empty intersections.** If the relations do not meet, the stable roots are compromises
+  between them, weighted by β. They are still returned; nothing flags them yet.
+- **No correction steps.** The MCMC corrections of Du et al. are for sampling and are not built
+  (I8).
+
+## 5. Tests
+
+`test/runtests.jl`, "products of relations: intersections":
+
+- the field equals the weighted sum;
+- a non-tuple `ps` is rejected;
+- two and four roots within 0.05 of the exact intersections;
+- the adjoint for the clamp and for one ellipse parameter matches central differences
+  (rtol 1e-6);
+- the parameter cotangent is a tuple.
@@ -667,6 +667,41 @@ end
     @test 1.1 < answer(remember) < 1.35          # averages the whole history
 end
 
+@testset "products of relations: intersections" begin
+    curve(f, n) = reduce(hcat, [collect(f(θ)) for θ in range(0, 2π; length = n + 1)[1:n]])
+    nodes = field_nodes(Xoshiro(1), 2; samples = 8)
+    function model(c)
+        p = NoisePredictor(GaussianMixtureEps(SCHED, c; s = 0.05), SCHED)
+        return ImplicitDiffusion(p, nodes), LuxCore.setup(Xoshiro(0), p)...
+    end
+    m1, p1, s1 = model(curve(θ -> (cos(θ), sin(θ)), 64))                       # unit circle
+    m2, p2, s2 = model(reduce(hcat, [[u, u] for u in range(-1.6, 1.6; length = 65)]))   # line y = x
+    m3, p3, s3 = model(curve(θ -> (1.5cos(θ), 0.6sin(θ)), 96))                 # ellipse
+    # the field of a product is the weighted sum of the fields
+    prod = ProductRelation(m1, m2; weights = [1.0, 2.0])
+    z = [0.3, -0.2]
+    @test first(prior_field(prod, z, (p1, p2), (s1, s2))) ≈
+          first(prior_field(m1, z, p1, s1)) .+ 2 .* first(prior_field(m2, z, p2, s2))
+    @test_throws ArgumentError prior_field(prod, z, p1, s1)                    # ps must be a tuple
+    # circle ∩ line: two answers, circle ∩ ellipse: four, near the exact intersections
+    for (other, ps, st, exact) in ((m2, p2, s2, [[0.7071, 0.7071], [-0.7071, -0.7071]]),
+                                   (m3, p3, s3, [[sx * 0.873, sy * 0.488] for sx in (1, -1) for sy in (1, -1)]))
+        roots, _ = implicit_roots(ProductRelation(m1, other), [0.0, 0.0], [0.0, 0.0], (p1, ps), (s1, st);
+                                  nstarts = 40, spread = 1.5, maxiters = 400)
+        @test length(roots) == length(exact)
+        @test all(r -> minimum(norm(r.z .- e) for e in exact) < 0.05, roots)
+    end
+    # the adjoint through a product: w.r.t. the clamped input and a parameter of the second factor
+    pe = ProductRelation(m1, m3); ρ = [Inf, 0.0]; h = 1e-6
+    sol, _ = implicit_infer(pe, [0.95, 0.4], ρ, (p1, p3), (s1, s3); tol = 1e-12)
+    b = implicit_pullback(pe, sol, [0.95, 0.4], ρ, [0.0, 1.0], (p1, p3), (s1, s3))
+    y(x, μ) = first(implicit_infer(pe, [x, 0.4], ρ, (p1, (μ = μ,)), (s1, s3); tol = 1e-12)).z[2]
+    @test b.z₀[1] ≈ (y(0.95 + h, p3.μ) - y(0.95 - h, p3.μ)) / 2h rtol = 1e-6
+    μ₊ = copy(p3.μ); μ₊[2, 10] += h; μ₋ = copy(p3.μ); μ₋[2, 10] -= h
+    @test b.ps[2].μ[2, 10] ≈ (y(0.95, μ₊) - y(0.95, μ₋)) / 2h rtol = 1e-6
+    @test b.ps isa Tuple && length(b.ps) == 2
+end
+
 @testset "all answers of a query, and their Laplace uncertainty" begin
     # Gaussian data: the curvature is RED-Diff's κ, so the Laplace covariance has a closed form
     v0 = 0.5