From e390826d4d0b4d08a9c82eec95de44dd62aeccaa Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 18 Jul 2026 15:39:40 +0000 Subject: [PATCH 1/6] feat: add jump_operator method, 9 rigorous tests, and README explainer --- README.md | 71 ++++++++++ jump_diffusion_engine/engine.py | 51 +++++++ tests/test_engine.py | 243 ++++++++++++++++++++++++++++++++ 3 files changed, 365 insertions(+) diff --git a/README.md b/README.md index 399a862..6a18b38 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,77 @@ A simulation framework for systems with a Source (Λ(t)), a Medium (Δ), and a nonlinear Sink (f(Δ)) subject to continuous diffusion and discrete jump noise. +## The Jump Operator 𝒥 + +Every Markov generator decomposes uniquely (**Lévy–Khintchine**) as: + +``` +ℒ = b ∂ₓ + ½σ² ∂ₓ² + 𝒥 + ──────────────── ─ + local nonlocal +``` + +There is no fourth term. Any process that moves is drift, diffusion, or jump — and 𝒥 is the only term that evaluates *f at the destination* rather than derivatives at the origin. + +### Definition + +``` +𝒥f(x) = λ(x) ∫ [f(x+z) − f(x)] ν(dz|x) +``` + +| Symbol | Role | This engine | +|:---------|:-----|:------------| +| `λ(x)` | **Rate** — *when* jumps occur | `jump_rate` parameter | +| `ν(dz\|x)` | **Kernel** — *where* to land | `jump_size_dist` parameter (default: N(0,1)) | + +Use `engine.jump_operator(f, x)` to evaluate 𝒥 numerically for any test function. + +### Confirmed properties + +The implementation is verified by the `TestJumpOperator` suite: + +| Property | Result | +|:---------|:-------| +| **Annihilates constants** | `f = c → 𝒥f = 0` (exactly) | +| **Zero without jumps** | `jump_rate = 0 → 𝒥f = 0` (exactly) | +| **Linear in `f`** | `𝒥(αf + βg) = α𝒥f + β𝒥g` | +| **Quadratic identity** | `f(x) = x²`, kernel N(0,1): `𝒥f(x) = λ` (independent of `x`) | +| **Symmetric-kernel identity** | `f(x) = x`, kernel N(0,1): `𝒥f(x) = 0` | +| **Scales with rate** | Doubling `λ` doubles `𝒥f` | +| **Respects custom kernel** | `f(x)=x`, Exp(1) kernel: `𝒥f(x) = λ·E[Z] = λ` | + +### The flat-obstruction property + +For the smooth flat function `f(x) = e^{−1/x}` (x > 0), `f(0) = 0`: + +``` +b f′(0) + ½σ² f″(0) → 0 (all derivatives vanish at the flat point) +𝒥f(0⁺) > 0 (jump kernel has positive weight on x > 0) +``` + +Local operators are blind to a flat wall. **𝒥 is not.** This is why nonlocal dynamics can cross obstructions that purely diffusive noise cannot. + +### The crossing condition + +In a bistable potential the basin edge sits at distance `J_c = Δ_edge − Δ*` from the equilibrium. + +``` +J < J_c → P(escape) ≈ 0 (jump cannot clear the wall) +J > J_c → P(escape) ≫ 0 (jump clears in a single step) +``` + +This is a **cliff, not a slope**. The `test_crossing_condition_cliff_not_slope` test confirms empirically: `P_below < 0.15`, `P_above > 0.50`, gap `> 0.35`. + +### SDE form + +``` +dΔₜ = −Γ′(Δ*)(Δₜ − Δ*) dt + σ dWₜ + J dNₜ(λ) +``` + +- `dt` holds the trajectory in the bowl. +- `dW` explores continuously. +- **`dN` arrives** — and when `J > J_c`, it crosses. + ## Core Equation `dΔ = [Λ(t) − f(Δ)] dt + σ dW + J dN` diff --git a/jump_diffusion_engine/engine.py b/jump_diffusion_engine/engine.py index c17a850..dfb8245 100644 --- a/jump_diffusion_engine/engine.py +++ b/jump_diffusion_engine/engine.py @@ -466,6 +466,57 @@ def seat_and_release(self, t_max: float, x0: float = 0.0, 'contained_after_release': contained, } + def jump_operator(self, f: Callable[[float], float], x: float, + n_samples: int = 50_000, + seed: Optional[int] = 0) -> float: + """ + Evaluate the jump operator at a single point via Monte Carlo integration. + + Definition (Lévy–Khintchine nonlocal part): + + 𝒥f(x) = λ(x) ∫ [f(x+z) − f(x)] ν(dz|x) + + where λ(x) is the jump rate and ν(dz|x) is the jump kernel (where to + land). For the default Gaussian kernel ν = N(0,1) this reduces to: + + 𝒥f(x) = jump_rate · E_Z[f(x+Z) − f(x)], Z ~ N(0,1) + + Key properties (each covered by a test): + • Annihilates constants: f = c → 𝒥f = 0 + • Zero without jumps: jump_rate = 0 → 𝒥f = 0 + • Linear in f: 𝒥(αf + βg) = α𝒥f + β𝒥g + • Quadratic identity: f(x)=x², kernel N(0,1) → 𝒥f(x) = λ + • Flat obstruction (non-zero at smooth flat points): + f(x) = e^{-1/x} (x>0), f(0)=0 — all derivatives vanish at 0, + yet 𝒥f(0⁺) > 0. Local operators (drift, diffusion) go to zero; + the nonlocal operator does not. + + Parameters + ---------- + f : test function callable, f: float → float + x : evaluation point + n_samples : MC sample count (higher → lower variance) + seed : integer seed for the MC RNG (None → use self.rng) + + Returns + ------- + float : Monte Carlo estimate of 𝒥f(x) + """ + lam = self.jump_rate(x) if callable(self.jump_rate) else float(self.jump_rate) + if lam == 0.0: + return 0.0 + + rng = np.random.default_rng(seed) if seed is not None else self.rng + + if self.jump_size_dist is None: + z = rng.normal(0.0, 1.0, size=n_samples) + else: + z = np.array([self.jump_size_dist() for _ in range(n_samples)]) + + fx = f(x) + increments = np.vectorize(f)(x + z) - fx + return float(lam * np.mean(increments)) + def stationary_density(self, lambda_val: float, x_range: Tuple[float, float] = (-10, 10), n_points: int = 2000): x = np.linspace(x_range[0], x_range[1], n_points) V = self.potential(x, lambda_val) diff --git a/tests/test_engine.py b/tests/test_engine.py index e977a8f..f880359 100644 --- a/tests/test_engine.py +++ b/tests/test_engine.py @@ -230,6 +230,249 @@ def test_energy_not_recorded_when_skipped(self, engine): assert 'energy' not in results[0] +# --------------------------------------------------------------------------- +# jump_operator +# --------------------------------------------------------------------------- + +class TestJumpOperator: + """ + Rigorous tests for 𝒥f(x) = λ ∫ [f(x+z) − f(x)] ν(dz|x). + + All MC tests use seed=0 and n_samples=100_000 for tight tolerances. + """ + + # --- helpers shared across tests --- + + def _engine_with_jumps(self, jump_rate=0.5): + """Engine with Gaussian jump kernel N(0,1), fixed seed.""" + return JumpDiffusionEngine( + lambda_func=lambda t: 0.5, + sigma=0.3, + jump_rate=jump_rate, + dt=0.01, + seed=42, + ) + + N = 100_000 # MC samples — enough for 3-decimal accuracy on smooth f + + # 1. Annihilates constants ------------------------------------------------- + + def test_constant_function_gives_zero(self): + eng = self._engine_with_jumps() + f = lambda x: 7.0 + for x in [-2.0, 0.0, 1.5, 3.0]: + result = eng.jump_operator(f, x, n_samples=self.N, seed=0) + assert abs(result) < 1e-10, ( + f"Jf should be 0 for constant f, got {result} at x={x}" + ) + + # 2. Zero when jump_rate = 0 ----------------------------------------------- + + def test_zero_jump_rate_gives_zero(self): + eng = self._engine_with_jumps(jump_rate=0.0) + f = lambda x: x**2 + np.sin(x) + for x in [-1.0, 0.0, 2.5]: + result = eng.jump_operator(f, x, n_samples=self.N, seed=0) + assert result == 0.0, ( + f"Jf must be exactly 0 when jump_rate=0, got {result}" + ) + + # 3. Linearity: 𝒥(αf + βg) = α𝒥f + β𝒥g ------------------------------------ + + def test_linearity_in_f(self): + eng = self._engine_with_jumps() + f = lambda x: x**2 + g = lambda x: np.cos(x) + alpha, beta = 3.0, -2.0 + h = lambda x: alpha * f(x) + beta * g(x) + + for x in [0.5, 1.0, 2.0]: + Jf = eng.jump_operator(f, x, n_samples=self.N, seed=0) + Jg = eng.jump_operator(g, x, n_samples=self.N, seed=0) + Jh = eng.jump_operator(h, x, n_samples=self.N, seed=0) + expected = alpha * Jf + beta * Jg + assert abs(Jh - expected) < 0.02, ( + f"Linearity failed at x={x}: J(αf+βg)={Jh:.6f}, " + f"αJf+βJg={expected:.6f}" + ) + + # 4. Quadratic identity (analytical): f=x², N(0,1) kernel → 𝒥f(x) = λ ----- + + def test_quadratic_identity(self): + """ + For f(x)=x² and Z~N(0,1): + 𝒥f(x) = λ · E[(x+Z)² − x²] + = λ · E[2xZ + Z²] + = λ · (2x · 0 + 1) = λ + Result is independent of x. + """ + lam = 0.5 + eng = self._engine_with_jumps(jump_rate=lam) + f = lambda x: x**2 + for x in [-3.0, 0.0, 1.0, 4.0]: + result = eng.jump_operator(f, x, n_samples=self.N, seed=0) + assert abs(result - lam) < 0.02, ( + f"Quadratic identity: expected {lam}, got {result:.5f} at x={x}" + ) + + # 5. Linear f has zero jump operator under symmetric kernel ---------------- + + def test_linear_function_zero_under_symmetric_kernel(self): + """ + For f(x)=x and Z~N(0,1) (symmetric, mean 0): + 𝒥f(x) = λ · E[(x+Z) − x] = λ · E[Z] = 0 + """ + eng = self._engine_with_jumps() + f = lambda x: x + for x in [-1.0, 0.0, 2.0]: + result = eng.jump_operator(f, x, n_samples=self.N, seed=0) + assert abs(result) < 0.02, ( + f"Linear f under symmetric kernel: expected 0, got {result:.5f} at x={x}" + ) + + # 6. Flat obstruction: 𝒥f(0⁺) ≠ 0 even though f'(0)=0 (all orders) -------- + + def test_flat_obstruction_nonzero(self): + """ + f(x) = e^{-1/x} (x>0), 0 otherwise — a flat obstruction. + + All derivatives vanish at x=0 → drift and diffusion terms go to zero. + Yet 𝒥f(0⁺) > 0 because the jump kernel puts weight on x>0 where f>0. + This is the fundamental non-locality of 𝒥. + """ + eng = self._engine_with_jumps(jump_rate=1.0) + + def flat_f(x): + return float(np.exp(-1.0 / x)) if x > 1e-10 else 0.0 + + # At x→0⁺ with Gaussian jumps, about half the destinations are positive. + # E[f(0+Z)] = E[e^{-1/Z} · 1_{Z>0}] > 0 → 𝒥f(0⁺) > 0. + result = eng.jump_operator(flat_f, 0.0, n_samples=self.N, seed=0) + assert result > 0.0, ( + f"Flat obstruction: 𝒥f(0⁺) must be > 0, got {result:.6f}" + ) + + # Confirm that local terms (drift b·f' and diffusion ½σ²f'') are 0 at x=0. + # All derivatives of e^{-1/x} vanish as x→0⁺. + h = 1e-6 + f_prime_approx = (flat_f(h) - flat_f(0.0)) / h + assert abs(f_prime_approx) < 1e-100, ( + f"f'(0⁺) should be 0 (flat point), got {f_prime_approx}" + ) + + # 7. Crossing condition: cliff, not slope ---------------------------------- + + def test_crossing_condition_cliff_not_slope(self): + """ + Demonstrates the crossing condition J > J_c = Δ_edge − Δ*. + + Uses a fixed-size jump distribution: all jumps equal J (deterministic). + Below J_c: escape_probability stays low (≤ 0.15). + Above J_c: escape_probability jumps sharply (≥ 0.50). + + The cliff shape — not a gradual slope — reflects that the nonlocal + operator either reaches over the basin wall or it does not. + + Uses a bistable parameter regime (k=0.5, g=2.0, K=1.0) so the basin + walls are well-defined unstable fixed points. + """ + bistable_kwargs = dict(k=0.5, g=2.0, K=1.0) + eng_base = JumpDiffusionEngine( + lambda_func=lambda t: 0.5, + sigma=0.05, + jump_rate=0.2, + dt=0.01, + seed=42, + **bistable_kwargs, + ) + boundary = eng_base.identify_boundary(lambda_val=0.5) + x_star = boundary['x_star'] + half_width = boundary['half_width'] + + if half_width is None or x_star is None: + pytest.skip("No well-defined basin found for crossing test") + + J_c = half_width # minimum jump to clear the basin edge + + results = {} + for scale in [0.4, 1.8]: + J = scale * J_c + eng = JumpDiffusionEngine( + lambda_func=lambda t: 0.5, + sigma=0.05, + jump_rate=0.2, + jump_size_dist=lambda _J=J: _J, # deterministic jump size + dt=0.01, + seed=42, + **bistable_kwargs, + ) + results[scale] = eng.escape_probability( + threshold=J_c * 0.9, + t_max=20.0, + x0=x_star, + x_star=x_star, + n_trials=80, + ) + + p_below = results[0.4] + p_above = results[1.8] + assert p_below < 0.15, ( + f"Crossing condition (below J_c): expected P < 0.15, got {p_below:.3f}" + ) + assert p_above > 0.50, ( + f"Crossing condition (above J_c): expected P > 0.50, got {p_above:.3f}" + ) + assert p_above > p_below + 0.35, ( + f"Cliff: gap P_above−P_below should be > 0.35, got {p_above - p_below:.3f}" + ) + + # 8. Scaling with jump_rate ------------------------------------------------ + + def test_scales_linearly_with_jump_rate(self): + """𝒥f(x) is proportional to λ (rate doubles → operator doubles).""" + f = lambda x: x**3 + x = 1.5 + lam1, lam2 = 0.5, 1.0 + eng1 = self._engine_with_jumps(jump_rate=lam1) + eng2 = self._engine_with_jumps(jump_rate=lam2) + J1 = eng1.jump_operator(f, x, n_samples=self.N, seed=0) + J2 = eng2.jump_operator(f, x, n_samples=self.N, seed=0) + # Both use the same seed → same z draws → ratio should equal lam2/lam1 + ratio = J2 / J1 if abs(J1) > 1e-12 else None + if ratio is not None: + assert abs(ratio - (lam2 / lam1)) < 0.05, ( + f"Expected ratio {lam2/lam1}, got {ratio:.4f}" + ) + + # 9. Custom kernel is respected ------------------------------------------- + + def test_custom_kernel_respected(self): + """With a custom jump_size_dist, 𝒥 uses that distribution.""" + # Positive-only jumps: Z ~ Exp(1). For f(x)=x: + # 𝒥f(x) = λ · E[Z] = λ · 1 = λ + lam = 0.6 + rng_inner = np.random.default_rng(99) + + def exp_kernel(): + return rng_inner.exponential(1.0) + + eng = JumpDiffusionEngine( + lambda_func=lambda t: 0.5, + sigma=0.3, + jump_rate=lam, + jump_size_dist=exp_kernel, + dt=0.01, + seed=42, + ) + f = lambda x: x + x = 2.0 + result = eng.jump_operator(f, x, n_samples=self.N, seed=None) + # E[Z] = 1 for Exp(1) → 𝒥f(x)=λ + assert abs(result - lam) < 0.05, ( + f"Custom Exp(1) kernel: expected 𝒥f ≈ {lam}, got {result:.4f}" + ) + + # --------------------------------------------------------------------------- # plot_trajectories (non-display smoke test) # --------------------------------------------------------------------------- From e749739142c330a2e1a3896512688b423ea8bf52 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 18 Jul 2026 15:41:36 +0000 Subject: [PATCH 2/6] docs: move jump operator explainer under Core Equation, trim to essentials --- README.md | 77 +++++++------------------------------------- tests/test_engine.py | 6 +++- 2 files changed, 16 insertions(+), 67 deletions(-) diff --git a/README.md b/README.md index 6a18b38..a581904 100644 --- a/README.md +++ b/README.md @@ -8,84 +8,29 @@ A simulation framework for systems with a Source (Λ(t)), a Medium (Δ), and a nonlinear Sink (f(Δ)) subject to continuous diffusion and discrete jump noise. -## The Jump Operator 𝒥 +## Core Equation -Every Markov generator decomposes uniquely (**Lévy–Khintchine**) as: +`dΔ = [Λ(t) − f(Δ)] dt + σ dW + J dN` -``` -ℒ = b ∂ₓ + ½σ² ∂ₓ² + 𝒥 - ──────────────── ─ - local nonlocal -``` +- `Λ(t) = ε₀ + A·sin(ωt)` — the source (constant or time-varying) +- `f(Δ) = kΔ + gΔ²/(K²+Δ²)` — the nonlinear sink (linear + saturating) +- `Δ* : Λ = f(Δ), f′(Δ*) > 0` — a stable equilibrium (basin centre) -There is no fourth term. Any process that moves is drift, diffusion, or jump — and 𝒥 is the only term that evaluates *f at the destination* rather than derivatives at the origin. +### The Jump Operator 𝒥 -### Definition +Every Markov generator decomposes uniquely (**Lévy–Khintchine**) as `ℒ = b∂ₓ + ½σ²∂ₓ² + 𝒥` — drift, diffusion, jump, no fourth term. 𝒥 is the only term that evaluates *f at the destination* rather than derivatives at the origin. ``` 𝒥f(x) = λ(x) ∫ [f(x+z) − f(x)] ν(dz|x) ``` -| Symbol | Role | This engine | -|:---------|:-----|:------------| -| `λ(x)` | **Rate** — *when* jumps occur | `jump_rate` parameter | -| `ν(dz\|x)` | **Kernel** — *where* to land | `jump_size_dist` parameter (default: N(0,1)) | - -Use `engine.jump_operator(f, x)` to evaluate 𝒥 numerically for any test function. - -### Confirmed properties - -The implementation is verified by the `TestJumpOperator` suite: - -| Property | Result | -|:---------|:-------| -| **Annihilates constants** | `f = c → 𝒥f = 0` (exactly) | -| **Zero without jumps** | `jump_rate = 0 → 𝒥f = 0` (exactly) | -| **Linear in `f`** | `𝒥(αf + βg) = α𝒥f + β𝒥g` | -| **Quadratic identity** | `f(x) = x²`, kernel N(0,1): `𝒥f(x) = λ` (independent of `x`) | -| **Symmetric-kernel identity** | `f(x) = x`, kernel N(0,1): `𝒥f(x) = 0` | -| **Scales with rate** | Doubling `λ` doubles `𝒥f` | -| **Respects custom kernel** | `f(x)=x`, Exp(1) kernel: `𝒥f(x) = λ·E[Z] = λ` | - -### The flat-obstruction property - -For the smooth flat function `f(x) = e^{−1/x}` (x > 0), `f(0) = 0`: - -``` -b f′(0) + ½σ² f″(0) → 0 (all derivatives vanish at the flat point) -𝒥f(0⁺) > 0 (jump kernel has positive weight on x > 0) -``` - -Local operators are blind to a flat wall. **𝒥 is not.** This is why nonlocal dynamics can cross obstructions that purely diffusive noise cannot. +`λ(x)` is the jump rate (`jump_rate`); `ν(dz|x)` is the jump kernel (`jump_size_dist`, default N(0,1)). Use `engine.jump_operator(f, x)` to evaluate 𝒥 numerically. -### The crossing condition +**Confirmed properties** (verified by `TestJumpOperator`): annihilates constants · linear in `f` · quadratic identity `𝒥(x²) = λ` for N(0,1) kernel · zero under symmetric kernel for linear `f` · scales with rate · respects custom kernel. -In a bistable potential the basin edge sits at distance `J_c = Δ_edge − Δ*` from the equilibrium. +**Flat obstruction** — at a flat point `f(x) = e^{−1/x}` all derivatives vanish, so drift and diffusion go to zero. `𝒥f(0⁺) > 0` regardless. Local operators are blind to a flat wall; `𝒥` is not. -``` -J < J_c → P(escape) ≈ 0 (jump cannot clear the wall) -J > J_c → P(escape) ≫ 0 (jump clears in a single step) -``` - -This is a **cliff, not a slope**. The `test_crossing_condition_cliff_not_slope` test confirms empirically: `P_below < 0.15`, `P_above > 0.50`, gap `> 0.35`. - -### SDE form - -``` -dΔₜ = −Γ′(Δ*)(Δₜ − Δ*) dt + σ dWₜ + J dNₜ(λ) -``` - -- `dt` holds the trajectory in the bowl. -- `dW` explores continuously. -- **`dN` arrives** — and when `J > J_c`, it crosses. - -## Core Equation - -`dΔ = [Λ(t) − f(Δ)] dt + σ dW + J dN` - -- `Λ(t) = ε₀ + A·sin(ωt)` — the source (constant or time-varying) -- `f(Δ) = kΔ + gΔ²/(K²+Δ²)` — the nonlinear sink (linear + saturating) -- `Δ* : Λ = f(Δ), f′(Δ*) > 0` — a stable equilibrium (basin centre) +**Crossing condition** — basin half-width `J_c = Δ_edge − Δ*`. Below `J_c`: `P(escape) ≈ 0`. Above: `P(escape) ≫ 0`. A cliff, not a slope. Use the `jump_diffusion_engine/` package to **analyse** stochastic systems and **steer** trajectories toward stable basins: diff --git a/tests/test_engine.py b/tests/test_engine.py index f880359..4958419 100644 --- a/tests/test_engine.py +++ b/tests/test_engine.py @@ -397,11 +397,15 @@ def test_crossing_condition_cliff_not_slope(self): results = {} for scale in [0.4, 1.8]: J = scale * J_c + + def _make_jump(size): + return lambda: size + eng = JumpDiffusionEngine( lambda_func=lambda t: 0.5, sigma=0.05, jump_rate=0.2, - jump_size_dist=lambda _J=J: _J, # deterministic jump size + jump_size_dist=_make_jump(J), # deterministic jump size dt=0.01, seed=42, **bistable_kwargs, From 975cce162301e7d659e74b0f68be1fd99e50a73e Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 18 Jul 2026 16:17:43 +0000 Subject: [PATCH 3/6] feat: add markov_generator with safer column-generator approach and tests --- jump_diffusion_engine/engine.py | 86 +++++++++++++++++++++++++++++++++ tests/test_engine.py | 77 +++++++++++++++++++++++++++++ 2 files changed, 163 insertions(+) diff --git a/jump_diffusion_engine/engine.py b/jump_diffusion_engine/engine.py index dfb8245..bd76b26 100644 --- a/jump_diffusion_engine/engine.py +++ b/jump_diffusion_engine/engine.py @@ -517,6 +517,92 @@ def jump_operator(self, f: Callable[[float], float], x: float, increments = np.vectorize(f)(x + z) - fx return float(lam * np.mean(increments)) + def markov_generator(self, lambda_val: float, + x_range: Tuple[float, float] = (-10, 10), + n_points: int = 200) -> Dict: + """ + Build the Markov generator matrix L for the discretized SDE. + + Discretises the state space into ``n_points`` grid nodes and constructs + the (n_points × n_points) transition-rate matrix using the upwind + Fokker–Planck finite-difference scheme. The diagonal is set via the + *safer column-generator approach*:: + + L[n, n] = -(rate_up + rate_down) + + so that each column sums to exactly zero by construction, rather than + being derived as a residual (which can accumulate floating-point error). + + Off-diagonal entries (non-zero only for nearest neighbours): + + * ``L[n+1, n] = rate_up[n]`` — upward transition from state n + * ``L[n-1, n] = rate_down[n]`` — downward transition from state n + + Transition rates are computed from the drift ``μ(x) = Λ − f(x)`` and + diffusion coefficient ``σ`` via the standard upwind decomposition: + + .. math:: + + r_\\mathrm{up}[n] = \\max(\\mu(x_n),\\,0)/dx + \\sigma^2/(2\\,dx^2) + + r_\\mathrm{down}[n] = \\max(-\\mu(x_n),\\,0)/dx + \\sigma^2/(2\\,dx^2) + + Boundary nodes (n=0 and n=N−1) are treated as absorbing: their + off-grid rates are dropped and the diagonal is adjusted accordingly. + + Parameters + ---------- + lambda_val : float + Constant forcing Λ used to evaluate the drift μ = Λ − f(Δ). + x_range : (float, float) + Extent of the discretised state space. + n_points : int + Number of grid nodes (matrix size n_points × n_points). + + Returns + ------- + dict with keys: + + * ``'L'`` — (n_points, n_points) ndarray, the generator matrix + * ``'x'`` — (n_points,) ndarray, grid node positions + * ``'dx'`` — float, uniform grid spacing + * ``'rate_up'`` — (n_points,) ndarray, upward rates at each node + * ``'rate_down'`` — (n_points,) ndarray, downward rates at each node + """ + x = np.linspace(x_range[0], x_range[1], n_points) + dx = x[1] - x[0] + + # Drift at each grid node: μ(x) = Λ − f(x) + mu = np.array([lambda_val - self.f_func(xi) for xi in x]) + + # Upwind rates (non-negative by construction) + rate_up = np.maximum(mu, 0.0) / dx + self.sigma**2 / (2.0 * dx**2) + rate_down = np.maximum(-mu, 0.0) / dx + self.sigma**2 / (2.0 * dx**2) + + # Boundary nodes: no off-grid transitions + rate_up[-1] = 0.0 + rate_down[0] = 0.0 + + # --- Safer column-generator approach --- + # Fill super- and sub-diagonal first, then set diagonal explicitly so + # each column sums to exactly zero without relying on residual arithmetic. + L = np.zeros((n_points, n_points)) + + for n in range(n_points): + if n < n_points - 1: + L[n + 1, n] = rate_up[n] # upward neighbour receives this rate + if n > 0: + 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 + + return { + 'L': L, + 'x': x, + 'dx': float(dx), + 'rate_up': rate_up, + 'rate_down': rate_down, + } + def stationary_density(self, lambda_val: float, x_range: Tuple[float, float] = (-10, 10), n_points: int = 2000): x = np.linspace(x_range[0], x_range[1], n_points) V = self.potential(x, lambda_val) diff --git a/tests/test_engine.py b/tests/test_engine.py index 4958419..60a9f35 100644 --- a/tests/test_engine.py +++ b/tests/test_engine.py @@ -498,3 +498,80 @@ def test_accepts_single_realization(self, engine, monkeypatch): record_energy=False) fig = engine.plot_trajectories(results) assert fig is not None + + +# --------------------------------------------------------------------------- +# markov_generator +# --------------------------------------------------------------------------- + +class TestMarkovGenerator: + """Tests for the safer column-generator approach to building L.""" + + N = 50 # small grid for fast tests + + def test_returns_expected_keys(self, engine): + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + for key in ('L', 'x', 'dx', 'rate_up', 'rate_down'): + assert key in out, f"Missing key: {key}" + + def test_matrix_shape(self, engine): + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + L = out['L'] + assert L.shape == (self.N, self.N), f"Expected ({self.N},{self.N}), got {L.shape}" + + def test_column_sums_are_zero(self, engine): + """Safer column-generator: each column must sum to exactly zero.""" + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + col_sums = out['L'].sum(axis=0) + np.testing.assert_allclose( + col_sums, 0.0, atol=1e-12, + err_msg="Column sums of generator matrix are not zero" + ) + + def test_diagonal_equals_negative_rate_sum(self, engine): + """L[n, n] == -(rate_up[n] + rate_down[n]) for every n.""" + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + L = out['L'] + diag = np.diag(L) + expected = -(out['rate_up'] + out['rate_down']) + np.testing.assert_allclose( + diag, expected, atol=1e-14, + err_msg="Diagonal entries do not equal -(rate_up + rate_down)" + ) + + def test_off_diagonal_non_negative(self, engine): + """All off-diagonal entries must be ≥ 0 (valid transition rates).""" + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + L = out['L'] + off_diag = L - np.diag(np.diag(L)) + assert np.all(off_diag >= -1e-14), "Off-diagonal entries contain negative values" + + def test_tridiagonal_structure(self, engine): + """Generator matrix for a diffusion should be tridiagonal.""" + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + L = out['L'] + # All entries more than one step from the diagonal must be zero. + for i in range(self.N): + for j in range(self.N): + if abs(i - j) > 1: + assert L[i, j] == 0.0, ( + f"Non-zero entry at L[{i},{j}]={L[i,j]} (not tridiagonal)" + ) + + def test_boundary_rates_are_zero(self, engine): + """Absorbing boundaries: rate_up[-1] and rate_down[0] must be zero.""" + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + assert out['rate_up'][-1] == 0.0, "rate_up at last node should be zero" + assert out['rate_down'][0] == 0.0, "rate_down at first node should be zero" + + def test_grid_length_matches_n_points(self, engine): + n = 30 + out = engine.markov_generator(lambda_val=0.5, n_points=n) + assert len(out['x']) == n + assert len(out['rate_up']) == n + assert len(out['rate_down']) == n + + def test_dx_matches_grid(self, engine): + out = engine.markov_generator(lambda_val=0.5, n_points=self.N) + expected_dx = (out['x'][-1] - out['x'][0]) / (self.N - 1) + assert abs(out['dx'] - expected_dx) < 1e-12 From 0509bc17667dfe02496d1d603c7ea820d625a7eb Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 18 Jul 2026 16:18:28 +0000 Subject: [PATCH 4/6] fix: address code review feedback (spelling, assert style) --- jump_diffusion_engine/engine.py | 2 +- tests/test_engine.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/jump_diffusion_engine/engine.py b/jump_diffusion_engine/engine.py index bd76b26..7b43a62 100644 --- a/jump_diffusion_engine/engine.py +++ b/jump_diffusion_engine/engine.py @@ -523,7 +523,7 @@ def markov_generator(self, lambda_val: float, """ Build the Markov generator matrix L for the discretized SDE. - Discretises the state space into ``n_points`` grid nodes and constructs + Discretizes the state space into ``n_points`` grid nodes and constructs the (n_points × n_points) transition-rate matrix using the upwind Fokker–Planck finite-difference scheme. The diagonal is set via the *safer column-generator approach*:: diff --git a/tests/test_engine.py b/tests/test_engine.py index 60a9f35..7b98779 100644 --- a/tests/test_engine.py +++ b/tests/test_engine.py @@ -574,4 +574,4 @@ def test_grid_length_matches_n_points(self, engine): def test_dx_matches_grid(self, engine): out = engine.markov_generator(lambda_val=0.5, n_points=self.N) expected_dx = (out['x'][-1] - out['x'][0]) / (self.N - 1) - assert abs(out['dx'] - expected_dx) < 1e-12 + np.testing.assert_allclose(out['dx'], expected_dx, atol=1e-12) From 13c36a6939870415c06b1703afbf3e4490846de8 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 18 Jul 2026 16:26:57 +0000 Subject: [PATCH 5/6] examples: add master_equation.py demonstrating Pauli master equation via markov_generator --- examples/master_equation.py | 158 ++++++++++++++++++++++++++++++++++++ 1 file changed, 158 insertions(+) create mode 100644 examples/master_equation.py diff --git a/examples/master_equation.py b/examples/master_equation.py new file mode 100644 index 0000000..3cffed5 --- /dev/null +++ b/examples/master_equation.py @@ -0,0 +1,158 @@ +""" +master_equation.py — Pauli master equation via markov_generator + +Demonstrates that the generator matrix L built by JumpDiffusionEngine.markov_generator +is the exact matrix encoding of the Pauli (Kolmogorov forward) master equation: + + dp_n/dt = Σ_{m≠n} [ W_{n←m} p_m − W_{m←n} p_n ] + ^^^^^^^^^^^^^ ^^^^^^^^^^^^^ + in out + +In matrix form this is simply: + + dp/dt = L · p + +where L[n, m] = W_{n←m} (off-diagonal, ≥ 0) and + L[n, n] = −Σ_{m≠n} W_{m←n} = −(rate_up[n] + rate_down[n]) + +The diagonal is set by the *safer column-generator approach*: + + L[n, n] = −(rate_up[n] + rate_down[n]) + +so each column of L sums to exactly zero, which is the matrix-level statement +of probability conservation: d/dt Σ_n p_n = 0. + +This script: + 1. Builds L for the default engine (constant Λ, nonlinear sink f). + 2. Starts from a narrow Gaussian initial distribution p(0). + 3. Integrates dp/dt = L·p forward in time using scipy.integrate.solve_ivp. + 4. Plots the evolving distribution and verifies that probability is conserved + at every output time. + +Run: + python3 master_equation.py +""" + +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp + +from jump_diffusion_engine import JumpDiffusionEngine + + +# --------------------------------------------------------------------------- +# 1. Set up the engine +# --------------------------------------------------------------------------- +LAMBDA_VAL = 0.5 # constant forcing Λ +SIGMA = 0.4 # diffusion strength +K = 0.8 # linear restoring rate +G = 0.5 # nonlinear gain +K_SAT = 2.0 # saturation scale + +eng = JumpDiffusionEngine( + lambda_func=lambda t: LAMBDA_VAL, + sigma=SIGMA, + jump_rate=0.0, # pure diffusion + drift for clarity + dt=0.01, + seed=42, + k=K, g=G, K=K_SAT, +) + +# --------------------------------------------------------------------------- +# 2. Build the Markov generator matrix +# --------------------------------------------------------------------------- +X_RANGE = (-8.0, 12.0) +N = 300 # grid nodes (matrix is N × N) + +out = eng.markov_generator(lambda_val=LAMBDA_VAL, x_range=X_RANGE, n_points=N) +L = out['L'] # (N, N) generator — columns sum to zero +x = out['x'] # grid positions +dx = out['dx'] + +# Verify the column-sum property (should be ~machine epsilon) +max_col_err = np.abs(L.sum(axis=0)).max() +print(f"Max |column sum| of L: {max_col_err:.2e} (should be < 1e-12)") + +# --------------------------------------------------------------------------- +# 3. Initial distribution — narrow Gaussian centred at x = 2.0 +# --------------------------------------------------------------------------- +x0_dist = 2.0 # initial mean +sig0 = 0.4 # initial width + +p0 = np.exp(-0.5 * ((x - x0_dist) / sig0) ** 2) +_trapz = np.trapezoid if hasattr(np, 'trapezoid') else np.trapz +p0 /= _trapz(p0, x) # normalise to a proper density + +# --------------------------------------------------------------------------- +# 4. Integrate dp/dt = L · p forward in time +# --------------------------------------------------------------------------- +T_MAX = 4.0 +T_REPORT = [0.0, 0.5, 1.0, 2.0, 4.0] # snapshot times + +def rhs(t, p): + """dp/dt = L p""" + return L @ p + +sol = solve_ivp( + rhs, + t_span=(0.0, T_MAX), + y0=p0, + method='RK45', + t_eval=T_REPORT, + rtol=1e-6, + atol=1e-9, +) + +assert sol.success, f"ODE solver failed: {sol.message}" + +# --------------------------------------------------------------------------- +# 5. Verify probability conservation at every snapshot +# --------------------------------------------------------------------------- +print("\nProbability conservation check:") +print(f"{'Time':>8s} {'∫p dx':>12s} {'max p':>10s}") +for i, t in enumerate(sol.t): + p_t = sol.y[:, i] + norm = _trapz(p_t, x) + print(f"{t:8.2f} {norm:12.6f} {p_t.max():10.6f}") + +# --------------------------------------------------------------------------- +# 6. Plot the evolving distribution +# --------------------------------------------------------------------------- +fig, axes = plt.subplots(1, 2, figsize=(12, 4)) + +# Left: probability distributions at each snapshot +ax = axes[0] +colours = plt.cm.viridis(np.linspace(0.1, 0.9, len(sol.t))) +for i, t in enumerate(sol.t): + ax.plot(x, sol.y[:, i], color=colours[i], label=f"t = {t:.1f}") + +# Mark the stable fixed point +fps = eng.find_fixed_points(LAMBDA_VAL) +stable = [fp for fp in fps if fp['stable']] +if stable: + xs = stable[0]['x_star'] + ax.axvline(xs, color='crimson', linestyle='--', linewidth=1.2, + label=f"Δ* = {xs:.2f}") + +ax.set_xlabel("State Δ") +ax.set_ylabel("Probability density p(Δ, t)") +ax.set_title("Master equation: dp/dt = L · p") +ax.legend(fontsize=8) +ax.set_xlim(X_RANGE) +ax.set_ylim(bottom=0) + +# Right: probability norm over time (conservation check) +norms = [_trapz(sol.y[:, i], x) for i in range(len(sol.t))] +ax2 = axes[1] +ax2.plot(sol.t, norms, 'o-', color='steelblue') +ax2.axhline(1.0, color='gray', linestyle='--', linewidth=0.8) +ax2.set_xlabel("Time t") +ax2.set_ylabel("∫ p(Δ, t) dΔ") +ax2.set_title("Probability conservation (should stay = 1)") +ax2.set_ylim(0.98, 1.02) + +plt.tight_layout() +fig.savefig("/tmp/master_equation.png", dpi=120) +print("\nFigure saved to /tmp/master_equation.png") From 8b4423f877ae40b6e34c9975dc7f69e6d5e532e6 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 18 Jul 2026 16:33:32 +0000 Subject: [PATCH 6/6] examples: add fock_ladder_generator.py (column generator, pair processes, parity, detailed balance) --- examples/fock_ladder_generator.py | 167 ++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 examples/fock_ladder_generator.py diff --git a/examples/fock_ladder_generator.py b/examples/fock_ladder_generator.py new file mode 100644 index 0000000..846096f --- /dev/null +++ b/examples/fock_ladder_generator.py @@ -0,0 +1,167 @@ +""" +fock_ladder_generator.py — Column generator for a Fock-space ladder + +Demonstrates five properties of the Markov generator matrix L for a +truncated quantum harmonic oscillator / photon-number ladder: + + dp/dt = L p, with the column-generator convention 1^T L = 0 + +The single-step transitions are standard bosonic ladder rates: + W_{n-1 <- n} = gd * n (emission / stimulated decay) + W_{n+1 <- n} = gu * (n+1) (absorption / stimulated creation) + +Optional pair processes add dn=±2 jumps: + W_{n-2 <- n} = pair_dn * n*(n-1) (two-photon loss) + W_{n+2 <- n} = pair_up * (n+1)*(n+2) (two-photon gain) + +Demonstrations +-------------- +1. COLUMN GENERATOR, dn=±1 only + - Columns sum to zero (1^T L = 0) + - Matrix is tridiagonal (bandwidth ≤ 1) + - Stationary distribution is geometric: p_n ∝ r^n, r = gu/gd + +2. BOUNDARY TERMS — truncation matters + - Lower wall n=0: emission rate gd*0 = 0 → natural reflecting wall + - Upper wall n=N-1: absorption would leave space, so it is dropped → reflecting + - n_bar converges to r/(1-r) as N → ∞; finite-N bias is truncation, not physics + +3. PAIR PROCESSES: dn=±2 widens the kernel support + - Bandwidth becomes 2; 1^T L = 0 still holds + +4. PARITY: pure dn=±2 splits Fock space into even/odd sectors + - Off-diagonal even↔odd block is exactly zero: parity is conserved + - Squeezing / two-photon processes cannot change photon-number parity + +5. DETAILED BALANCE with both dn=±1 and dn=±2 present + - Net current on each bond J_{a→b} = L[b,a]*p[a] − L[a,b]*p[b] + - Non-zero currents signal irreversibility when competing processes coexist + +Run: + python3 fock_ladder_generator.py +""" + +import numpy as np + + +def L_gen(N, gd, gu, pair_up=0.0, pair_dn=0.0): + """ + Build the Markov generator matrix for a truncated Fock-space ladder. + + Convention: L[n, m] = W_{n←m} (off-diagonal, ≥ 0), + L[m, m] = −Σ_{n≠m} W_{n←m} (diagonal, ≤ 0). + + This is the *column-generator* form: dp/dt = L p, 1^T L = 0. + + Parameters + ---------- + N : int — number of Fock states {0, 1, …, N-1} + gd : float — single-photon loss rate (emission) + gu : float — single-photon gain rate (absorption) + pair_up : float — two-photon gain rate (dn = +2) + pair_dn : float — two-photon loss rate (dn = −2) + + Returns + ------- + L : (N, N) ndarray + """ + L = np.zeros((N, N)) + + def add(dest, src, rate): + """Add a transition src→dest with given rate (ignored if out of bounds).""" + if 0 <= dest < N and rate > 0: + L[dest, src] += rate + L[src, src] -= rate + + for n in range(N): + add(n - 1, n, gd * n) # emission: dn = −1 + add(n + 1, n, gu * (n + 1)) # absorption: dn = +1 + add(n - 2, n, pair_dn * n * (n - 1)) # pair loss: dn = −2 + add(n + 2, n, pair_up * (n + 1) * (n + 2)) # pair gain: dn = +2 + + return L + + +def stationary(L): + """Return the normalised stationary distribution (eigenvector for eigenvalue 0).""" + w, v = np.linalg.eig(L) + i = np.argmin(np.abs(w)) + p = np.real(v[:, i]) + p /= p.sum() + return p + + +def net_current(p, L, a, b): + """Net probability current on bond a↔b: J = L[b,a]*p[a] − L[a,b]*p[b].""" + return L[b, a] * p[a] - L[a, b] * p[b] + + +# --------------------------------------------------------------------------- +# 1. COLUMN GENERATOR, dn=±1 only +# --------------------------------------------------------------------------- +N = 8 +gd, gu = 1.0, 0.3 + +L = L_gen(N, gd, gu) +print("1. COLUMN GENERATOR, dn=+-1 only") +print(" 1^T L = ", np.round(L.sum(0), 12)) + +idx = np.abs(np.subtract.outer(np.arange(N), np.arange(N))) +bw_ok = np.all(idx[np.abs(L) > 1e-12] <= 1) +print(f" bandwidth: nonzero only on |n-m|<=1 -> {bw_ok}") + +p = stationary(L) +r = gu / gd +th = r ** np.arange(N) +th /= th.sum() +print(f" stationary matches geometric r^n: {np.allclose(p, th)}") + +# --------------------------------------------------------------------------- +# 2. BOUNDARY TERMS — truncation matters +# --------------------------------------------------------------------------- +print("\n2. BOUNDARY TERMS — the truncation matters") +print(" at n=0: emission rate gd*0 = 0 (no leak below vacuum) -> reflecting OK") +print(f" at n=N-1: absorption gu*N = {gu*N:.2f} would leave the space; we DROP it.") +print(" dropping = reflecting wall. This is why n_bar was slightly below r/(1-r):") + +nbar = np.sum(np.arange(N) * p) +print(f" n_bar(N={N}) = {nbar:.6f} vs r/(1-r) = {r/(1-r):.6f}") + +for NN in (8, 16, 32, 64): + LL = L_gen(NN, gd, gu) + pp = stationary(LL) + print(f" N={NN:>3}: n_bar={np.sum(np.arange(NN)*pp):.9f}") + +print(f" -> converges to {r/(1-r):.9f} as the wall recedes. Truncation, not physics.") + +# --------------------------------------------------------------------------- +# 3. PAIR PROCESSES: dn=±2 widens the support of W +# --------------------------------------------------------------------------- +print("\n3. PAIR PROCESSES: dn=+-2 widens the support of W") +L2 = L_gen(N, gd, gu, pair_up=0.02, pair_dn=0.05) +bw = np.max(idx[np.abs(L2) > 1e-12]) +print(f" bandwidth now = {bw} (was 1). Kernel support = {{+-1, +-2}}") +print(" 1^T L = ", np.round(L2.sum(0), 12)) + +# --------------------------------------------------------------------------- +# 4. PARITY: pure dn=±2 splits Fock space into even/odd sectors +# --------------------------------------------------------------------------- +print("\n4. PARITY: pure dn=+-2 splits Fock space into even/odd sectors") +Lp = L_gen(N, 0, 0, pair_up=0.02, pair_dn=0.05) +ev = np.arange(0, N, 2) +od = np.arange(1, N, 2) +cross = np.abs(Lp[np.ix_(od, ev)]).max() +print(f" coupling even->odd block: {cross:.3e}") +print(" -> exactly zero: parity is CONSERVED. Two disconnected ladders.") +print(" Squeezing/two-photon processes cannot change photon-number parity.") + +# --------------------------------------------------------------------------- +# 5. DETAILED BALANCE with dn=±1 and dn=±2 present +# --------------------------------------------------------------------------- +print("\n5. DETAILED BALANCE with dn=+-2 present") +p2 = stationary(L2) +print(" net current on each bond (0 = reversible):") +for n in range(4): + j1 = net_current(p2, L2, n, n + 1) + j2 = net_current(p2, L2, n, n + 2) + print(f" {n}->{n+1}: {j1:+.3e} {n}->{n+2}: {j2:+.3e}")