Expose ring-generator utilities and surface non-ergodic Markov generators

jump_diffusion_engine/__init__.py

@@ -1,3 +1,5 @@
-from .engine import JumpDiffusionEngine
+from .engine import (JumpDiffusionEngine, ergodic_components,
+                     ring_generator, reduce_ring)
 
-__all__ = ["JumpDiffusionEngine"]
+__all__ = ["JumpDiffusionEngine", "ergodic_components",
+           "ring_generator", "reduce_ring"]

jump_diffusion_engine/engine.py

@@ -1,8 +1,97 @@
+import warnings
+from math import gcd as _gcd
+
 import numpy as np
 from scipy.optimize import root_scalar
 import matplotlib.pyplot as plt
 from typing import Callable, List, Optional, Tuple, Dict
 
+
+def ergodic_components(L: np.ndarray, rtol: float = 1e-9) -> int:
+    """Number of ergodic components of a column-generator matrix.
+
+    For a conservative generator (columns sum to zero), dim ker(L) equals
+    the number of closed communicating classes: each component carries its
+    own stationary distribution, so the zero eigenvalue has one copy per
+    component. Computed via SVD (robust to non-normality): the count of
+    singular values below ``rtol * s_max``.
+
+    Discrete check of the gcd law: for a ring of N sites with jump channels
+    {k_i}, dim ker(L) = gcd(k_1, ..., k_m, N). Anything downstream that
+    assumes a *unique* stationary distribution (stationary densities,
+    transfer times, spectral gaps) silently requires this to return 1.
+    """
+    s = np.linalg.svd(np.asarray(L, dtype=float), compute_uv=False)
+    if s.size == 0:
+        return 0
+    smax = s[0] if s[0] > 0 else 1.0
+    return int(np.sum(s <= rtol * smax))
+
+
+def ring_generator(N: int, channels: Dict[int, float]) -> np.ndarray:
+    """Conservative column-generator for jump channels on the ring Z_N.
+
+    ``channels`` maps a (signed) step dn to a rate. Steps are taken mod N;
+    dn ≡ 0 (mod N) is a no-op. Columns sum to zero by construction.
+    """
+    L = np.zeros((N, N))
+    for i in range(N):
+        for dn, r in channels.items():
+            j = (i + dn) % N
+            if j != i and r > 0:
+                L[j, i] += r
+                L[i, i] -= r
+    return L
+
+
+def reduce_ring(N: int, channels: Dict[int, float]) -> Dict:
+    """Exact decimation of a ring jump generator to its coprime core.
+
+    Since gcd(m·k, N) = m·gcd(k, N/m) for m | N, the cosets mod
+    g = GCD(all steps, N) are exactly invariant, and
+
+        L_N(channels)  ≅  I_g ⊗ L_{N/g}(channels/g)      (permutation-similar)
+
+    so the full spectrum is the core spectrum with multiplicity g, and
+    dim ker(L_N) = g. Returns the core parameters; analyses (spectra,
+    mixing rates, transfer times) can be run on the (N/g)-site core and
+    replicated g times at zero error. If no *active* channels remain after
+    applying the same acceptance rule as ``ring_generator`` (positive rate and
+    non-zero modular step), then no jumps occur and the reduction yields
+    ``g = N`` and ``N_core = 1``.
+    When multiple channels map to the same reduced step, their rates are
+    aggregated in ``channels_core``.
+
+    Returns
+    -------
+    dict with keys 'g' (sector count / multiplicity), 'N_core',
+    'channels_core', 'L_core'.
+    """
+    mod_n = lambda dn: ((dn % N) + N) % N
+    active = []
+    for dn, r in channels.items():
+        dn_mod = mod_n(dn)
+        if r > 0 and dn_mod != 0:
+            active.append((dn_mod, r))
+
+    if not active:
+        g = N
+        N_core = 1
+        ch_core = {}
+    else:
+        g = N
+        for dn_mod, _ in active:
+            g = _gcd(g, dn_mod)
+        N_core = N // g
+        ch_core = {}
+        for dn_mod, r in active:
+            dn_c = dn_mod // g
+            ch_core[dn_c] = ch_core.get(dn_c, 0.0) + r
+
+    return {'g': g, 'N_core': N_core, 'channels_core': ch_core,
+            'L_core': ring_generator(N_core, ch_core)}
+
+
 class JumpDiffusionEngine:
     """
     Universal Jump-Diffusion Simulator & Stochastic Stability Framework.
@@ -595,12 +684,22 @@ def markov_generator(self, lambda_val: float,
                 L[n - 1, n] = rate_down[n]  # downward neighbour receives this rate
             L[n, n] = -(rate_up[n] + rate_down[n])  # diagonal: exact negative sum
 
+        n_comp = ergodic_components(L)
+        if n_comp != 1:
+            warnings.warn(
+                f"markov_generator: generator has {n_comp} ergodic components; "
+                "stationary quantities and transfer times are ill-defined "
+                "(initial-condition dependent).",
+                RuntimeWarning
+            )
+
         return {
             'L': L,
             'x': x,
             'dx': float(dx),
             'rate_up': rate_up,
             'rate_down': rate_down,
+            'n_components': n_comp,
         }
 
     def stationary_density(self, lambda_val: float, x_range: Tuple[float, float] = (-10, 10), n_points: int = 2000):

tests/test_engine.py

@@ -6,7 +6,7 @@
 import numpy as np
 import pytest
 
-from jump_diffusion_engine import JumpDiffusionEngine
+from jump_diffusion_engine import JumpDiffusionEngine, reduce_ring, ring_generator
 
 
 # ---------------------------------------------------------------------------
@@ -500,6 +500,62 @@ def test_accepts_single_realization(self, engine, monkeypatch):
         assert fig is not None
 
 
+# ---------------------------------------------------------------------------
+# reduce_ring
+# ---------------------------------------------------------------------------
+
+class TestReduceRing:
+    def test_gcd_uses_only_active_channels(self):
+        out = reduce_ring(12, {2: 1.0, 3: -5.0, 6: 0.0})
+        assert out['g'] == 2
+        assert out['N_core'] == 6
+        assert out['channels_core'] == {1: 1.0}
+
+    def test_effectively_empty_channels_reduce_to_trivial_core(self):
+        out = reduce_ring(12, {1: 0.0, -1: -2.0, 12: 4.0})
+        assert out['g'] == 12
+        assert out['N_core'] == 1
+        assert out['channels_core'] == {}
+        np.testing.assert_allclose(out['L_core'], np.zeros((1, 1)))
+
+
+# ---------------------------------------------------------------------------
+# entropy production
+# ---------------------------------------------------------------------------
+
+class TestEntropyProduction:
+    def test_entropy_production_equals_current_times_affinity(self):
+        N = 7
+        rate_fwd = 1.7
+        rate_bwd = 0.4
+        L = ring_generator(N, {1: rate_fwd, -1: rate_bwd})
+        A = np.vstack([L, np.ones(N)])
+        b = np.zeros(N + 1)
+        b[-1] = 1.0
+        p = np.linalg.lstsq(A, b, rcond=None)[0]
+        np.testing.assert_allclose(L @ p, 0.0, rtol=0, atol=1e-12)
+        np.testing.assert_allclose(p.sum(), 1.0, rtol=0, atol=1e-12)
+        assert np.all(p > 0.0)
+
+        entropy_production = 0.0
+        bond_currents = []
+        affinity = 0.0
+        for i in range(N):
+            j = (i + 1) % N
+            f_ij = L[j, i] * p[i]
+            f_ji = L[i, j] * p[j]
+            assert f_ij > 0.0 and f_ji > 0.0
+            bond_current = f_ij - f_ji
+            bond_affinity = np.log(f_ij / f_ji)
+            bond_currents.append(bond_current)
+            affinity += bond_affinity
+            entropy_production += bond_current * bond_affinity
+
+        current = float(np.mean(bond_currents))
+        np.testing.assert_allclose(bond_currents, current, rtol=0, atol=1e-12)
+        np.testing.assert_allclose(entropy_production, current * affinity, rtol=0, atol=1e-12)
+
+
 # ---------------------------------------------------------------------------
 # markov_generator
 # ---------------------------------------------------------------------------