sft-wick

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

pip install sft-wick

For development (editable install with the test suite):

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:

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.

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):

sft-wick quickstart

The file it writes is examples/quickstart.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:

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:

¹ 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:

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

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.

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 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:

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

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.

# 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.

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

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 φ-ψ

# 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

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:

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

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

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

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

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

Classes

Expression Types

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:

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:

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:

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

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

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:

# 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.

# 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:

@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 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.

License

sft-wick is released under the BSD 3-Clause License. See LICENSE.