Metadata-Version: 2.5
Name: sft-wick
Version: 0.6.0
Summary: Wick's theorem contractions for statistical field theory perturbative calculations
Project-URL: Homepage, https://github.com/StatFieldTheory/sft-wick
Project-URL: Repository, https://github.com/StatFieldTheory/sft-wick
Project-URL: Documentation, https://sft-wick.readthedocs.io
Author-email: Zheng Zhang <zzhang0123@gmail.com>
License: BSD 3-Clause License
        
        Copyright (c) 2026-present, Zheng Zhang
        
        Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
        
        1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
        
        2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
        
        3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
        
        THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
License-File: LICENSE
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: BSD License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.10
Requires-Dist: joblib>=1.3
Requires-Dist: matplotlib>=3.7
Requires-Dist: networkx>=3.0
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: pyyaml>=6.0
Requires-Dist: scipy>=1.10
Requires-Dist: tabulate>=0.9
Provides-Extra: dev
Requires-Dist: mpmath>=1.3; extra == 'dev'
Requires-Dist: pytest-cov; extra == 'dev'
Requires-Dist: pytest>=7.0; extra == 'dev'
Provides-Extra: docs
Requires-Dist: furo>=2024.1; extra == 'docs'
Requires-Dist: ipykernel; extra == 'docs'
Requires-Dist: nbsphinx>=0.9; extra == 'docs'
Requires-Dist: sphinx-copybutton>=0.5; extra == 'docs'
Requires-Dist: sphinx>=7.0; extra == 'docs'
Provides-Extra: parallel
Requires-Dist: joblib>=1.3; extra == 'parallel'
Provides-Extra: progress
Requires-Dist: tqdm>=4.60; extra == 'progress'
Description-Content-Type: text/markdown

# sft-wick

[![PyPI](https://img.shields.io/pypi/v/sft-wick)](https://pypi.org/project/sft-wick/)
[![Python versions](https://img.shields.io/pypi/pyversions/sft-wick)](https://pypi.org/project/sft-wick/)
[![Documentation](https://img.shields.io/badge/docs-readthedocs-blue)](https://sft-wick.readthedocs.io)
[![arXiv](https://img.shields.io/badge/arXiv-2606.19480-b31b1b)](https://arxiv.org/abs/2606.19480)
[![DOI](https://zenodo.org/badge/1162069108.svg)](https://doi.org/10.5281/zenodo.20776358)
[![License: BSD-3-Clause](https://img.shields.io/badge/license-BSD--3--Clause-green)](LICENSE)

**Feynman-diagram expansion and evaluation for stochastic field theories.**

> 📖 **Documentation (API reference · user guide · theory background): <https://sft-wick.readthedocs.io>**

`sft-wick` automates perturbative calculations for **stochastic (partial) differential equations** in the Martin–Siggia–Rose (MSR) response-field formalism. The input is a Langevin-type field equation: a deterministic drift plus noise that may be non-Gaussian and spatially correlated. From it the package builds the interaction action, applies **Wick's theorem** to expand arbitrary field moments order by order, and writes every term using two two-point propagators, the correlation function *C* = ⟨φφ⟩ and the response (Green's) function *R* = ⟨φψ⟩. Diagrams are enumerated symbolically, rendered as Feynman graphs, and evaluated numerically, from a YAML config to theory-vs-simulation curves. For **self-consistent** (DMFT-style) problems, where the propagators define a self-energy that in turn defines the propagators, `solve_self_consistency` supplies the fixed-point iteration and its diagnostics.

## Installation

```bash
pip install sft-wick
```

For development (editable install with the test suite):

```bash
git clone https://github.com/StatFieldTheory/sft-wick.git
cd sft-wick
pip install -e ".[dev]"
```

Dependencies: `numpy`, `scipy`, `networkx`, `matplotlib`, `pandas`, `pyyaml`, `tabulate`, `joblib`. Optional: `pip install "sft-wick[progress]"` adds `tqdm` progress bars (plain stderr lines are printed without it). The `parallel` extra is kept for compatibility with older install commands. For development: `pytest`, `pytest-cov`.

The install also registers a CLI entry point:

```bash
sft-wick run config.yaml         # execute a full YAML-configured workflow
sft-wick run config.yaml --override sweep.seed=7 --dry-run
```

## Three-layer API: start with L2

`sft-wick` exposes three progressively higher-level entry points.
**The recommended entry point for any new analysis is L2: write a
YAML config, run it with the CLI, iterate.**  Drop down only when
the physics requires Python-level control.

| Layer | Entry point | Use when |
|---|---|---|
| **L2: YAML + CLI** ✨ | `sft-wick run config.yaml` | **Default choice.** Reproducible, shareable, diff-able; no Python code at the call site; runs identically on a laptop or a cluster; `--override` lets you scan parameters from the shell. |
| **L1: Python workflow** | `System`, `Expansion`, `Propagators`, `SweepResult` | You need to script custom pre/post-processing around the sweep, or compose multiple systems programmatically. |
| **L0: raw symbolic** | `compute_moment`, `Field`, `Vertex`, `Action`, `PropagatorCache`, `DiagramIntegrand` | You need fine-grained control over pairings, Itô flags, canonical forms, or custom simplifications. Research into the symbolic machinery. |

The Sphinx docs' "Workflow API" chapter (`docs/user_guide/workflow.rst`) covers L1/L2 end-to-end; the L0 reference below is the complete specification of the underlying symbolic machinery.

## Quick Start: L2 (config file)

One command writes a small config to the current directory and runs it
(a few seconds on a laptop, with a progress bar):

```bash
sft-wick quickstart
```

The file it writes is `examples/quickstart.yaml`:

```yaml
system:
  field: {name: phi, n_components: 2}
  linear: {type: diagonal, gamma: [1.0, 1.0]}
  vertices:
    - name: F
      coupling:                                    # bare F; MSR factor
        - [[0.0, 0.0], [0.0, 1.0]]                 # applied automatically
        - [[0.0, 0.5], [0.5, 0.0]]
  noise:
    kappa2:
      type: separable_translation
      temporal: {type: exponential, lam: 0.05, sigma_t: 0.3}
      spatial:  {type: exponential, sigma_x: 1.0}

expand:
  observable: ["phi_a(x)", "phi_b(y)"]
  orders: [0, 2]

propagators:
  t_max: 5.0
  n_grid_t: 60
  c_closed_form: auto      # built-in closed form for this kernel family

sweep:
  positions_grid: {x: [0.0], y: [0.0, 0.5, 1.5]}
  t_final_grid: [1.0, 5.0]
  component_pairs: [[0, 0], [1, 1]]
  n_samples: 4096
  seed: 42

output:
  - {type: table, format: markdown}
```

Run any config from the shell:

```bash
sft-wick run quickstart.yaml                   # full pipeline, progress bars
sft-wick run quickstart.yaml --override sweep.seed=7
sft-wick run quickstart.yaml --dry-run         # validate + cost estimate
sft-wick run quickstart.yaml --quiet           # no bars / banners
```

The next step up is `examples/demo1_config.yaml` (orders 0, 2 and 4;
four separations; the sweep behind the paper's Gaussian-noise figures),
then `examples/demo2_config.yaml` (non-local κ⁽³⁾ vertex with a dynamic
coupling), then `examples/demo3/config_FK.yaml` (filtered Poisson shot
noise, whose κ⁽³⁾ vertex is R-contracted exactly).

### How long will this take?

`sft-wick run CONFIG --dry-run` prints an estimate before you commit to a
run.  Measured wall-clock, serial, `pip install sft-wick` with no extras:

| Config | Apple M3 Ultra (one core) | Laptop proxy¹ | GitHub `ubuntu-latest` runner² | What it computes |
|---|---|---|---|---|
| `sft-wick quickstart` | 2 s | 10 s | 6 s | orders 0–2, 12 grid points, 4096 samples |
| `examples/demo1_config.yaml` | 50 s | 5.2 min | 1.5 min | orders 0–4 (71 diagrams), 16 grid points, 8192 samples |
| `examples/demo1/L2/config.yaml` (paper figures) | ~15 min (28 workers) | | | 672 grid points × 71 diagrams, `method: gauss_legendre` with `n_gauss: 24`, `sweep.n_jobs: -1` |
| `examples/demo2/L2/*.yaml` (κ⁽³⁾ figures) | 39 s (28 workers) | | | FF (QMC) + FK (R-contracted κ⁽³⁾, 1-D Gauss-Legendre), `sweep.n_jobs: -1` |
| `examples/demo3/config_FK.yaml` | 0.6 s | | | order-2 `Fκ³` channel (2 diagrams), 5 separations × 6 times |
| `examples/demo3/config_F3K.yaml` | 178 s (cold, `sweep.n_jobs: -1`) | | | the order-4 `F³κ³ + F³κ⁵` channels (36 diagrams, 3 times); 143 s of that is the one-off order-4 enumeration, cached afterwards |

¹ the same machine throttled with `taskpolicy -c background` and
`OMP_NUM_THREADS=1`; a 2023 MacBook Pro should be within a factor of two.
² the CI job `examples-time-gate` runs both on every push and fails above
600 s / 1200 s.

What drives the cost, in order: `expand.orders` (1 / 6 / 64 diagrams at
orders 0 / 2 / 4 for the cubic vertex), the number of grid points
(`sweep.positions_grid` × `t_final_grid` × `component_pairs`),
`sweep.n_samples` per diagram, and, only when no closed form applies,
`propagators.n_grid_t²` quadrature calls for the C table.  The exponential-
temporal kernel family gets a built-in closed form
(`propagators.c_closed_form: auto`), so its C table costs nothing.  That
covers the quick start, demo 1 and demo 2; demo 3 is the exception and
supplies its own (`c_closed_form_module`) for speed; its spatial envelope
has been expressible in YAML as `{type: custom}` since 0.6.0.
Kernels with neither are integrated by
Gauss-Legendre with a node count checked for convergence at the table's
extreme cells (`c_method: auto`).

## Quick Start: L1 (Python, for programmatic use)

The same workflow in Python, for use inside a larger script:

```python
import numpy as np
import sft_wick as sw

F = np.zeros((2, 2, 2))
F[0, 1, 1] = 1.0
F[1, 0, 1] = F[1, 1, 0] = 0.5

system = sw.System(
    field=sw.FieldSpec("phi", n_components=2),
    linear=sw.DiagonalA(gamma=[1.0, 1.0]),
    vertices=[sw.LocalVertex("F", coupling=F)],   # bare F
    noise=sw.GaussianNoise(kappa2=sw.SeparableTranslation(
        temporal=sw.ExponentialTemporal(lam=0.05, sigma_t=0.3),
        spatial=sw.ExponentialSpatial(sigma_x=1.0),
    )),
)

expansion = system.expand(("phi_a(x)", "phi_b(y)"), orders=[0, 2, 4])
props = system.propagators(t_max=15.0, n_grid_t=60)

sweep = expansion.sweep(
    props,
    positions_grid={"x": [0.0], "y": [0.0, 0.5, 1.0, 2.5]},
    t_final_grid=[1.0, 15.0],
    component_pairs=[(0, 0), (1, 1)],
)

print(sweep.totals())    # long-format pandas DataFrame
```

## Raw API (L0) Quick Start

```python
from sft_wick import Field, Vertex, Action, compute_moment

# Define scalar fields
phi = Field('phi', 'physical')
psi = Field('psi', 'response')

# Compute <psi(x) phi(x) phi(x) phi(x)>_{S_0}
obs = [psi('x'), phi('x'), phi('x'), phi('x')]
result = compute_moment(obs, Action(vertices=[]), order=0)
print(result.order(0).to_latex())
# Output: 3 R(x, x) C(x, x)
```

## Self-consistent solutions

Some problems are a **fixed point**: propagators define a self-energy, the
self-energy defines new propagators, repeat. `solve_self_consistency` runs
that loop, mixes, and reports what the loop did.

```python
from sft_wick import solve_self_consistency

def step(state):
    sigma = self_energy_from_diagrams(state)   # sft-wick computes this
    return dyson_solve(sigma)                  # you supply this

result = solve_self_consistency(initial, step, tol=1e-8, damping=0.3)
if not result:                       # bool(result) is result.converged
    raise RuntimeError(result.summary())
R, C = result.state
```

The **Dyson solve is not provided**: it is model-specific, and it is an
integral-equation solve rather than a diagram evaluation.

The result is not a bare state: a non-converged iteration looks like a
converged one if you print only the last state. It carries `converged`, the
full residual history, and a `reason`: `converged`, `diverged`,
`oscillating` (use damping) or `max_iter`. See
[the API page](https://sft-wick.readthedocs.io/en/latest/api/selfconsistency.html)
for the four ways a fixed-point loop can report a solution it never found,
and what this one does about each.

## Background

In the MSR formalism for SDEs, path-integral averages involve two types of fields:

- **Physical field** φ (phi): the field of interest
- **Response field** ψ (psi): the auxiliary conjugate field

The free two-point functions (propagators) are:

| Contraction | Propagator | Meaning |
|---|---|---|
| ⟨φ\_i(x) φ\_j(x')⟩\_{S₀} | C\_{ij}(x, x') | Correlation function |
| ⟨φ\_i(x) ψ\_j(x')⟩\_{S₀} | R\_{ij}(x, x') | Response (Green's) function |
| ⟨ψ\_i(x) ψ\_j(x')⟩\_{S₀} | 0 | Vanishes by construction |

Since the MSR partition function Z = 1, the perturbative expansion is simply:

```
⟨O⟩_S = Σ_{n=0}^{N} (-1)^n / n! ⟨O S_int^n⟩_{S_0}
```

Each term is evaluated via **Wick's theorem**: the expectation of a product of fields equals the sum over all complete pairings of the product of two-point functions.

## Usage Guide

### 1. Defining Fields

```python
from sft_wick import Field

# Scalar fields (single component)
phi = Field('phi', 'physical')
psi = Field('psi', 'response')

# Multi-component fields
phi = Field('phi', 'physical', n_components=3)
psi = Field('psi', 'response', n_components=3)
```

### 2. Creating Field Operators

Field operators are concrete instances with bound component indices and spatial arguments.

```python
# Scalar: phi(spatial_arg)
op = phi('x')          # φ(x)

# Multi-component: phi(component_index, spatial_arg)
op = phi('a', 'x')     # φ_a(x)
op = psi('b', 'y')     # ψ_b(y)
```

### 3. Defining Interaction Vertices

Vertices represent terms in the interaction action S\_int.

```python
from sft_wick import Vertex

# Local vertex: ∫ F_{ijk} φ_i(x) φ_j(x) ψ_k(x) dx
# All fields share the same spatial argument.
v1 = Vertex(fields=[phi, phi, psi], coupling='F')

# Non-local vertex: ∬ K_{ij}(x, x') ψ_i(x) ψ_j(x') dx dx'
# Each field gets its own spatial argument.
v2 = Vertex(fields=[psi, psi], coupling='K', local=False)
```

### 4. Computing Perturbative Expansions

```python
from sft_wick import Action, compute_moment

action = Action(vertices=[v1])
obs = [psi('a', 'x'), phi('b', 'x'), phi('c', 'x'), phi('d', 'x')]

result = compute_moment(obs, action, order=1)

# Access individual orders
print(result.order(0).to_latex())
print(result.order(1).to_latex())

# Full result
print(result.to_latex())
```

### 5. Feynman Diagrams

Each non-vanishing Wick contraction corresponds to a Feynman diagram:

- **Vertices** (■): interaction points from S\_int
- **External points** (●): observable field operators
- **C propagator** (blue solid line): correlation φ-φ
- **R propagator** (red dashed arrow): response φ-ψ

```python
# Draw all diagrams
result.draw_diagrams()

# Draw only diagrams at a specific order
result.draw_diagrams(order=1)

# Access diagram topology
for d_info in result.diagrams_by_order[1]:
    fd = d_info.to_feynman_diagram()
    print(fd.summary())
    print(f"  Loops: {fd.n_loops}, Connected: {fd.is_connected}")
```

### 6. LaTeX Formatting

```python
from sft_wick import LaTeXFormatter

# Default names
print(result.order(0).to_latex())
# C_{ab}(x, y)

# Custom propagator names
fmt = LaTeXFormatter(propagator_names={
    'C': 'G',
    'R': r'R^{\mathrm{ret}}'
})
print(fmt.format(result.order(0)))
# G_{ab}(x, y)

# LaTeX align environment for order-by-order display
print(fmt.format_aligned(result.order_terms))
```

### 7. Direct Wick Contraction

For low-level access without the perturbative machinery:

```python
from sft_wick import wick_contract, contract_pair

# Contract a product of fields
ops = [phi('a', 'x'), phi('b', 'y'), phi('c', 'z'), phi('d', 'w')]
expr, pairings = wick_contract(ops)
print(expr.to_latex())
# C_{ab}(x, y) C_{cd}(z, w) + C_{ac}(x, z) C_{bd}(y, w) + C_{ad}(x, w) C_{bc}(y, z)

# Contract a single pair
prop = contract_pair(phi('a', 'x'), psi('b', 'y'))
print(prop.to_latex())
# R_{ab}(x, y)
```

## Examples

### Zeroth-Order Moment

```python
phi = Field('phi', 'physical')
psi = Field('psi', 'response')

obs = [psi('x'), phi('x'), phi('x'), phi('x')]
result = compute_moment(obs, Action(vertices=[]), order=0)
print(result.order(0).to_latex())
# 3 R(x, x) C(x, x)
```

The three terms arise because ψ can pair with any of the three φ's (producing R), and the remaining two φ's pair together (producing C).

### First-Order Perturbation

```python
phi = Field('phi', 'physical')
psi = Field('psi', 'response')

v = Vertex(fields=[phi, psi], coupling='g')
action = Action(vertices=[v])

obs = [phi('x'), phi('y')]
result = compute_moment(obs, action, order=1)
print(result.order(1).to_latex())
# ∫ dy₀ (-g) [R(x, y₀) C(y, y₀) + R(y, y₀) C(x, y₀) + R(y₀, y₀) C(x, y)]
```

### Multi-Component Four-Point Function

```python
phi = Field('phi', 'physical', n_components=3)

obs = [phi('a', 'x'), phi('b', 'y'), phi('c', 'z'), phi('d', 'w')]
result = compute_moment(obs, Action(vertices=[]), order=0)
print(result.order(0).to_latex())
# C_{ab}(x, y) C_{cd}(z, w) + C_{ac}(x, z) C_{bd}(y, w) + C_{ad}(x, w) C_{bc}(y, z)
```

### Non-Local Interaction

```python
phi = Field('phi', 'physical', n_components=2)
psi = Field('psi', 'response', n_components=2)

v_nonlocal = Vertex(fields=[psi, psi], coupling='K', local=False)
action = Action(vertices=[v_nonlocal])

obs = [phi('a', 'x'), phi('b', 'y')]
result = compute_moment(obs, action, order=1)
```

## API Reference

### Core Functions

| Function | Description |
|---|---|
| `compute_moment(observable, action, order, ito=True, response_phase=True, collect_topology=True)` | Perturbative expansion of ⟨O⟩\_S up to given order |
| `compute_moment_numerical(observable, action, order, coupling_values, fixed_indices, ..., n_jobs=1)` | Fast numerical path using nauty canonical labeling for diagram grouping. Enables order-6 calculations. Requires `pynauty`. Parallelization uses `joblib` (`n_jobs=-1`). |
| `wick_contract(operators, ito=True)` | Apply Wick's theorem to a product of field operators |
| `contract_pair(op1, op2, ito=True)` | Contract two field operators into a propagator |
| `apply_response_phase(expr)` | Multiply each term by (−i)^n for n response propagators |
| `collect_by_diagram(expr)` | Group terms by Feynman diagram isomorphism, factor out propagators |
| `collect_by_topology(expr)` | Alias for `collect_by_diagram` (backward compat) |
| `integrate_moment(integrand, lambda_f, cache, ...)` | Integrate a single diagram's contribution (QMC, Gauss-Legendre or nquad) |
| `integrate_diagrams(diagram_terms, coupling_values, lambda_f, cache, ..., n_jobs=1)` | Batch-integrate a list of diagram terms, optionally in parallel (`n_jobs=-1`) |
| `simplify(expr)` | Simplify an expression (flatten, collect terms, eliminate zeros) |
| `reset_uid_counter()` | Reset field operator UID counter (for reproducible tests) |

### Classes

| Class | Description |
|---|---|
| `Field` | Field declaration (name, type, component count) |
| `FieldOperator` | Concrete field instance with bound index and position |
| `Vertex` | Interaction vertex template (local or non-local) |
| `VertexInstance` | Instantiated vertex with fresh internal indices |
| `Action` | Collection of vertices defining S\_int |
| `PerturbativeResult` | Result container with order-by-order expressions and diagrams |
| `FeynmanDiagram` | Graph representation of a diagram (networkx MultiGraph) |
| `DiagramRenderer` | Matplotlib-based diagram visualizer |
| `PropagatorModel` | User-provided R\_time and κ² callables for numerical evaluation |
| `PropagatorCache` | Caches C propagators (spline-interpolated or dblquad) |
| `DiagramIntegrand` | Ready-to-integrate object combining coupling coefficients and spatial structure |
| `LaTeXFormatter` | Configurable LaTeX output |
| `ImaginaryUnit` | The imaginary unit i, used in phase factors |

### Expression Types

| Type | Description | Example LaTeX |
|---|---|---|
| `Rational(num, den)` | Exact rational number | `\frac{1}{2}` |
| `Symbol(name, indices, spatial_args)` | Named tensor/coupling | `F_{ijk}` |
| `Propagator(kind, il, ir, sl, sr)` | Two-point function | `C_{ab}(x, y)` |
| `ImaginaryUnit()` | Imaginary unit | `\mathrm{i}` |
| `Sum(terms)` | Sum of expressions | `a + b + c` |
| `Product(factors)` | Product of expressions | `a b c` |
| `SumOverIndex(index, dim, body)` | Index summation | `\sum_{i=1}^{N} ...` |
| `IntegralOver(var, body)` | Spatial integration | `\int dx ...` |
| `KroneckerDelta(i, j)` | Component delta | `δ_{ij}` |
| `DiracDelta(x, y)` | Spatial delta | `δ(x - y)` |

## Conventions and Options

### Distinct external labels (breaking change in 0.4.0)

Each external operator must carry its own spatial label.  Since 0.4.0,
`("phi_a(x)", "phi_b(x)")` raises `ValueError` at `System.expand` (L1)
and `compute_moment` (L0), and since 0.4.2 at `compute_moment_numerical`
(L0), **at interacting orders**; order 0 is exempt and unchanged.
Upgrading from 0.3.x, give each external a distinct
label and put them at the same point through `positions`:

```python
expansion = system.expand(("phi_a(x)", "phi_b(y)"), orders=[0, 2])
expansion.evaluate(props, positions={"x": 0.0, "y": 0.0}, ...)
```

Coincident external *points* are supported; the shared *label* is not.
The spatial contraction is keyed by label, so two operators sharing one
were collapsed without their distinct component-index routings, giving a
silently wrong answer: measured on demo 2's system, the order-2 `F`
channel was low by a factor 2 while its `FK` channel was right.  The
0.4.0 entry in `CHANGELOG.md` says why the spelling is refused rather
than repaired with a multiplicity factor.

### Itô prescription (`ito=True`, default)

By default, the Itô discretisation convention Θ(0)=0 is applied:

- **Equal-point R vanishes**: R(x,x) = 0, which eliminates self-response contractions and intra-vertex tadpoles in local vertices.
- **Causal R-loops vanish**: Any closed loop of response propagators R(a,b)R(b,c)...R(z,a) = 0, since this would require a cyclic time ordering t\_a > t\_b > ... > t\_a, which is impossible for the retarded propagator.

Pass `ito=False` to keep these terms symbolic.  The numerical layer
evaluates every R at equal times as 0, i.e. it computes the Itô SDE.  On a
local vertex with one ψ leg that is exact: the equal-point term and the
Stratonovich functional Jacobian cancel, and sft-wick emits neither, so the
number is the `ito=True` one.  On a local vertex with two or more ψ legs (a
φ-dependent noise covariance) or between two external operators, the Itô and
Stratonovich values differ, and the numerical layer raises instead of
returning the Itô value.

A Stratonovich SDE is computed in its Itô form.  At L1,
`GaussianNoise(sigma2=MultiplicativeImpulse(g0, g1, interpretation=...))`
declares white noise with the amplitude `g(φ) = g0 + g1 φ`: `D0 = g0 g0ᵀ`
enters C, the rest of `D = g gᵀ` becomes local ψψφⁿ vertices (factor
`−i²/2! = ½`), and `interpretation='stratonovich'` adds the noise-induced
drift `½ Σ g ∂g` as a source (ψ) and a linear (ψφ) vertex.  See
`docs/user_guide/workflow.rst`.

### Response phase convention (`response_phase=True`, default)

The MSR convention ⟨φ(a) ψ(b)⟩ = −i R(a,b) is implemented by multiplying each term by (−i)^n, where n is the number of response propagators R in that term. This is applied after simplification so that like-term collection is unaffected.

Pass `response_phase=False` to get raw R propagators without the phase factor.

### Diagram-based term collection (`collect_topology=True`, default)

Terms whose Feynman diagrams are isomorphic are grouped together, under relabeling of dummy integration variables and accounting for C propagator symmetry C(x,y) = C(y,x). The algorithm computes a canonical graph form for each term by trying all permutations of internal spatial variables. Propagators are factored out with canonical component indices, and coupling coefficients are summed with appropriately permuted indices to produce expressions like (F\_{ijk} + F\_{ikj}) R C.

At second order and above, spatial-variable relabeling (e.g. y\_0 ↔ y\_1 for two copies of the same vertex) merges additional equivalent pairings.

Pass `collect_topology=False` to keep all pairings expanded individually.

## Design Notes

- **No SymPy dependency**: Uses a custom lightweight expression tree with `fractions.Fraction` for exact rational arithmetic.
- **Frozen dataclasses**: All expression types are immutable and hashable, safe for use in sets and dicts.
- **Unique operator IDs**: Each `FieldOperator` carries a unique integer ID, so that two copies of φ\_a(x) in the same product are properly distinguished during contraction.
- **Optimized contraction**: Two engines are available. The operator-level engine (`generate_valid_pairings`) skips ψ-ψ pairings at construction time. The spatial-level engine (`wick_contract_spatial`, used by default when `collect_topology=True`) enumerates spatial topologies instead of operator-level pairings, computing a multiplicity for each. This avoids the combinatorial explosion from component-index routing and is orders of magnitude faster at high perturbative orders.
- **Feynman diagrams**: Built on `networkx.MultiGraph` (supporting multiple edges between the same pair of nodes) with `matplotlib` rendering.

## Performance

### `compute_moment` (symbolic path)

The default `compute_moment` builds full symbolic expressions and groups diagrams via brute-force canonical form search (trying all k! permutations of internal spatial variables). Performance at each perturbative order:

| Order | Operators | Topologies | Diagrams | Time  |
|-------|-----------|------------|----------|-------|
| 2     | 8         | 12         | 6        | <0.01s |
| 4     | 14        | 1,416      | 64       | ~1.5s |
| 6     | 20        | 738,900    | 1,088    | infeasible (hours/OOM) |

### `compute_moment_numerical` (nauty path)

Replaces the O(k!) canonical form search with the nauty graph isomorphism algorithm (via `pynauty`), reducing diagram grouping from hours to seconds at order 6:

| Order | Topologies | Nauty grouping | Component routing | Total  |
|-------|------------|----------------|-------------------|--------|
| 4     | 1,416      | 0.02s          | 0.25s             | ~0.3s  |
| 6     | 738,900    | ~12s           | ~8 min            | ~10 min |

Pass `n_jobs=-1` to parallelize across CPU cores (`joblib` is installed by default).

### Known bottlenecks and future directions

At order 6, the dominant cost is **component routing** (`_enumerate_component_routings`), called once per spatial topology (738K calls). Potential improvements:

- **Cache routing for isomorphic topologies**: topologies in the same nauty canonical group are graph-isomorphic. If the nauty permutation can be applied at the operator level, routing need only run once per canonical group (~1K calls instead of ~738K). This requires mapping operator UIDs across isomorphic graphs.
- **Einsum-based coupling evaluation**: for constant coupling tensors, the coupling sum can be computed via `np.einsum` tensor contraction rather than symbolic expression evaluation, eliminating the combinatorial component-index enumeration entirely.
- **Vectorized QMC integration** ✅: `DiagramIntegrand.integrate_moment_qmc_vectorized()` replaces the Python for-loop over Sobol samples with batch propagator evaluation via `PropagatorCache.C_diagonal_batch()` and `R_time_batch()`. Achieves 18–22× speedup over the scalar `integrate_moment_qmc()` with identical results.
- **GPU acceleration**: the simulation (Euler–Maruyama) and QMC integration are embarrassingly parallel across realizations/samples and would benefit from JAX `vmap` or similar frameworks.

## Testing

```bash
pytest tests/ -v
```

Every test file, its test count (parametrised cases counted), what it
checks, the independent reference it is checked against and its
tolerance are listed in the generated validation catalogue,
`docs/verification/catalog.rst` (`python tools/gen_test_catalog.py`).
The core is organised into eight deductive phases:

1. Phase 1: Symbolic expansion (`test_deductive_expansion.py`)
2. Phase 2: Propagator numerics (`test_deductive_numerics.py::TestClosedFormC` etc.)
3. Phase 3: Full diagram evaluation
4. Phase 4: Alternative-path consistency (vectorised, parallel, nauty)
5. Phase 5: Spatial homogeneity modes (translation / rotation / general)
6. Phase 6: White-noise component
7. Phase 7: Dynamic non-local coupling + L1/L2 workflow round-trip
   (`test_workflow.py`, `test_workflow_config.py`)
8. Phase 8: Time-dependent linear operator (`test_diagonal_A_time_dependent.py`)

See `docs/verification/index.rst` for the per-phase test matrix, tolerances, and design rationale.

## Repository layout

| Path | Contents |
|------|----------|
| `src/sft_wick/` | Package source: diagram enumeration, propagators, numerical evaluation, drawing, and the `workflow/` high-level API + CLI |
| `examples/` | Worked examples: `demo1/` (Gaussian noise), `demo2/` (non-Gaussian, non-zero κ³), `demo3/` (filtered Poisson shot noise), `demo4/` (compound-Poisson noise asymmetric in points and components), `demo5/` (white noise on every integrator; multiplicative noise at L0 and L1, Itô and Stratonovich), `demo6/` (repeated and static non-local vertices, cubic plus quartic drift), `demo7/` (observables in space, angle and time), `demo8/` (time-dependent coefficients and non-exponential dynamics), `reference/` (`ito_moments.py`, the exact Itô moment hierarchy demos 4-8 are checked against; `hormander_moments.py`, the same hierarchy from the vector fields as written, for demo 5's multiplicative noise; `decaying_drift.py`, the reference behind `tests/test_local_callable_coupling.py`), and tutorial notebooks |
| `tests/` | pytest suite (eight deductive phases) |
| `docs/` | Sphinx documentation (ReadTheDocs source) |

## Worked examples (reproducible test runs)

Three of the eight demos are checked against a direct Langevin simulation,
each covering both symbolic diagram expansion and numerical evaluation.
All eight ship committed outputs; these three ship the cached simulation
as well:

```bash
# demo1: Gaussian driving noise
python examples/demo1/run_simulation.py            # writes sim_cache.npz (~50k realisations)
# then run examples/demo1/analysis.ipynb           # diagrams + theory-vs-simulation figures

# demo2: non-Gaussian noise (non-zero third cumulant kappa^3)
python examples/demo2/run_simulation.py --alpha 0.6   # non-Gaussian
python examples/demo2/run_simulation.py --alpha 0.0   # Gaussian control
# then run examples/demo2/analysis.ipynb           # kappa^3 cross-check + FK channel

# demo3: filtered Poisson (shot) noise, exact in the free-field limit
cd examples/demo3
python level_a.py                 # ~3 min   level A (F = 0): the exact check
python level_b.py                 # ~6 min   level B: Fkappa^3 + F^3kappa^3 + F^3kappa^5
python make_figures.py            #          figures + TikZ diagram sources
sft-wick run config_FK.yaml       # the same level-B physics through the L2 CLI
```

Demo 3's level A runs through the L1 Python API; both of demo 3's own
configs (`config_FK.yaml`, `config_F3K.yaml`) are level B.  An n-point
observable does go through YAML: `sweep.component_tuples` takes one index
per operator, as `examples/demo4/config_level_a.yaml` and
`config_level_a_4pt.yaml` do for a level-A triple and quadruple.  See
`examples/demo3/INTERPRETATION.md` for demo 3's validation ledger and
error budget.

The cached simulation outputs (`sim_cache.npz` for demos 1 and 2,
`level_a_results.npz` / `level_b_results.npz` for demo 3) and the
reference figures are committed, so a reviewer can re-run the scripts and
diff against the shipped outputs. See also
`examples/nonlocal_vertex_2pt.ipynb` for a non-local-vertex tutorial.

Demos 4 to 8 have an exact reference instead of a simulation: the moment
hierarchy of the Itô process at the observation points
(`examples/reference/ito_moments.py`, which imports nothing from
sft-wick), solved order by order in the couplings, so each package channel
is compared with one exact coefficient.  Demo 5's multiplicative-noise part
uses `examples/reference/hormander_moments.py`, which builds the same
hierarchy from the drift and the noise columns as written, in Hörmander
form for the Stratonovich reading, so the noise-induced drift is never
formed on the reference side.  A third module there,
`examples/reference/decaying_drift.py`, is a quadratic drift decaying in
time; it backs `tests/test_local_callable_coupling.py` rather than a demo.

```bash
# demo4: cumulants asymmetric in points and components (the leg-order defect's class)
cd examples/demo4 && python level_a.py && python level_b.py     # ~10 s
cd examples/demo4 && python poisson_level_b_order4.py           # ~3 min, the order-4 F^3kappa^3 channel
# demo5: white noise on every integrator; multiplicative noise at L0 and L1
cd examples/demo5 && python run.py && python multiplicative.py
cd examples/demo5 && python white_l1_multiplicative.py   # Itô and Stratonovich
# demo6: repeated and static non-local vertices, m = 5, quartic plus cubic drift
cd examples/demo6 && python vertex6_repeated.py && python vertex6_interacting.py
cd examples/demo6 && python vertex6_cubic.py && python vertex6_high_cumulants.py
cd examples/demo6 && python vertex6_gaussian_vertex.py
# demo7: two-time, angular and 3-D observables
cd examples/demo7 && python space7_run.py && python space7_shot3d.py
# demo8: a rate varying in time, a damped-oscillator R, custom kernels, sigma^2(t)
cd examples/demo8 && python time8_run.py
```

Each README (`examples/demo4/README.md` through
`examples/demo8/README.md`) lists the agreement per channel and
integrator, and what the same scripts report on the code before each fix
they cover.

## Documentation

Full documentation (API reference, user guide, theory background) is hosted at
<https://sft-wick.readthedocs.io>.

## Citation

If you use `sft-wick`, please cite the paper:

```bibtex
@misc{zhang2026sftwickformalismpackagefeynmandiagram,
      title={sft-wick: A formalism and package for Feynman-diagram expansion and evaluation in stochastic field theories},
      author={Zheng Zhang},
      year={2026},
      eprint={2606.19480},
      archivePrefix={arXiv},
      primaryClass={physics.comp-ph},
      url={https://arxiv.org/abs/2606.19480},
}
```

The software is additionally archived on Zenodo.  [DOI:10.5281/zenodo.20776358](https://doi.org/10.5281/zenodo.20776358) is the *concept* DOI: it covers all versions and always resolves to the most recent release.  To reference the *specific* version your results were produced with, cite that release's own version DOI, listed under "Versions" on the Zenodo record page.  For 0.5.0 that is [DOI:10.5281/zenodo.22715757](https://doi.org/10.5281/zenodo.22715757).

## License

`sft-wick` is released under the BSD 3-Clause License.  See [LICENSE](LICENSE).
