sft_wick.workflow — High-Level Workflow API (L1 + L2)

The workflow subpackage is the user-facing layer. Users who need to “declare physics → expand → integrate → inspect” should reach for these types first; drop down to the raw API only if a concrete requirement demands it.

Top-level system object

System — the top-level user-facing spec for a stochastic-field problem. Lowers high-level specification objects to raw PropagatorModel / Field / Vertex / Action / etc.

class sft_wick.workflow.system.System(field, linear, noise, vertices=<factory>, nonlocal_vertices=<factory>, t_min=0.0, explicit_R=None)

Bases: object

High-level specification of a stochastic field-theory problem.

Combines all the physics inputs the user needs to describe in a single immutable object. Use expand() and propagators() to get to the computational layer.

Parameters:

• field (FieldSpec) – FieldSpec — the physical field (φ) spec; the response field (ψ) is derived automatically.

• linear (LinearOp) – LinearOp variant — defines R. Required unless explicit_R is supplied.

• noise (GaussianNoise) – GaussianNoise — κ² (+ optional σ²) defines C together with R.

• vertices (tuple[LocalVertex, ...]) – list of LocalVertex — F^(n) local interactions. May be empty (linear theory).

• nonlocal_vertices (tuple[NonLocalVertex, ...]) – list of NonLocalVertex — κ^(m) for m ≥ 3 non-Gaussian driving contributions. May be empty.

• t_min (float) – lower time bound for propagator integrations.

• explicit_R (Callable | None) – escape hatch — if set, overrides linear. The structured alternative is the ExplicitR LinearOp variant; both work, but ExplicitR is preferred for serialisation and YAML configs.

field: FieldSpec

linear: LinearOp

noise: GaussianNoise

vertices: tuple[LocalVertex, ...]

nonlocal_vertices: tuple[NonLocalVertex, ...]

t_min: float = 0.0

explicit_R: Callable | None = None

property n_components: int

property homogeneity: str

Inferred from the noise κ² spec ('translation' / …).

property iso_R: bool

build_propagator_model(diag_C=True)

Lower the spec to a PropagatorModel.

Parameters:

diag_C (bool) – when True (default), the numerical C propagator is represented as a diagonal vector (n, N). When False, the full (n, N, N) matrix is preserved – required for observables that probe cross-component C entries (e.g. lensing kappa-gamma cross-correlation). Off-diagonal support is only meaningful in combination with Propagators.build(c_closed_form_only=True) because the spline-table paths build only diagonal entries.

Return type:

PropagatorModel

build_fields()

Return (phi, psi) Field objects.

Return type:

tuple[Field, Field]

build_action()

Lower vertices to a raw Action.

Return type:

Action

build_coupling_values()

Dict passed to DiagramTerm.evaluate_coupling / build_integrand.

The MSR prefactors are applied here:

• Local F^(n): multiplied by -i (F_MSR = -i F).

• Non-local κ^(m): multiplied by -(i^m) / m! (K_MSR = -(i^m)/m! κ^(m)).

Users pass the bare physical tensors when constructing LocalVertex / NonLocalVertex; this method is the single point of truth for the MSR convention. For a callable non-local coupling the returned dict value is a wrapped callable that applies the factor at evaluation time.

Return type:

dict[str, Any]

expand(observable, orders, *, response_phase=True, ito=True, collect_topology=True, iso_R=None, diag_R=True, diag_C=True, iso_C=False, cache_path=None)

Run the perturbative expansion and return an Expansion object.

Parameters:

• observable (Iterable) – iterable of either FieldOperator objects or strings like "phi_a(x)". Strings are parsed as field_compIndex(spatialArg).

• orders (Iterable[int]) – iterable of perturbative orders to compute.

• response_phase (bool) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• ito (bool) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• collect_topology (bool) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• iso_R (bool | None) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• diag_R (bool) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• diag_C (bool) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• iso_C (bool) – forwarded to compute_moment(). When iso_R is None the value is inferred from self.linear.

• cache_path (Any) – directory (or file) for on-disk caching. None disables caching (prints a one-shot reminder).

Return type:

Expansion

propagators(t_max, n_grid_t=60, *, homogeneity=None, r_max=None, n_grid_r=None, n_grid_cos=None, x_max=None, n_grid_x=None, n_jobs=1, c_closed_form=None, cache_path=None, interp_method='linear', c_closed_form_only=False, c_closed_form_vectorized=False, c_method='dblquad', c_n_gauss=20, diag_C=True)

Build a Propagators object with a precomputed spatial table matching self.homogeneity.

By default (all grid args None) the cache enters lazy mode for the inferred homogeneity — recommended for moment calculations at a small fixed set of external positions.

Parameters:

• t_max (float) – upper time bound.

• n_grid_t (int) – number of time grid points.

• homogeneity (str | None) – override the inferred homogeneity ('translation' | 'rotation' | 'general').

• r_max (float | None) – translation full-grid parameters (both given ⇒ full 3-D spline; either None ⇒ lazy).

• n_grid_r (int | None) – translation full-grid parameters (both given ⇒ full 3-D spline; either None ⇒ lazy).

• n_grid_cos (int | None) – rotation full-grid parameter (given ⇒ full 3-D spline over cos θ).

• x_max (float | None) – general full-grid parameters (both given ⇒ full 4-D spline).

• n_grid_x (int | None) – general full-grid parameters (both given ⇒ full 4-D spline).

• n_jobs (int) – parallel workers for grid build (-1 = all cores, via joblib).

• c_closed_form (Callable | None) – optional callable (n1, t1, n2, t2) -> (N, N) returning the full C matrix directly. When provided, the spline-table builder bypasses scipy.integrate.dblquad — use this when the user knows C in closed form (e.g. OU kernel).

• cache_path (Any) – directory (or file) for on-disk caching of the constructed PropagatorCache.

• interp_method (str) – RegularGridInterpolator method for the full-grid C tables. 'linear' (default) is monotone and safe for steep tails; 'cubic' gives O(h⁴) on smooth grids. Forwarded to PropagatorCache.

• c_closed_form_only (bool) – when True (with c_closed_form set), skip every spline and route C lookups directly through the user’s c_fn – machine-precision agreement with the analytical form. Forwarded to Propagators.build().

• c_closed_form_vectorized (bool) – c_fn accepts batched arrays and returns (n, N, N) (only meaningful with c_closed_form_only=True).

• c_method (str) –

  Quadrature method for the inner C-propagator ∫ R κ² R integral when the cache builds its table.

  • 'dblquad' (default) – adaptive Gauss-Kronrod; robust on any κ² but slow (10-80 ms / cell).

  • 'gauss_legendre' – tensor-product GL with a diagonal-aware sub-region split. Recommended for any κ² that is piecewise analytic – the package’s exponential / Gaussian / Legendre kernels all qualify. 18-100× faster than dblquad at machine precision (c_n_gauss=20). Not the right tool for non-smooth κ² (discontinuous, oscillatory at high frequency, integrable singularities like 1/√t).

  Ignored when c_closed_form_only=True (the user’s C_fn is exact and bypasses both quadrature paths). See Workflow API (L1 + L2) “Choosing an integrator” for the full decision matrix.

• c_n_gauss (int) – Per-dim GL node count for c_method='gauss_legendre' (default 20 → machine precision on smooth kernels). Cost c_n_gauss² per sub-region.

• diag_C (bool) – when True (default), the numerical C propagator is represented as a diagonal vector (n, N) – the bit-identical behaviour for every pre-existing caller. When False, the full (n, N, N) matrix is preserved so observables can read off-diagonal entries C[a, b] with a != b (e.g. the lensing κ-γ₊ cross-correlation). Only meaningful with c_closed_form_only=True: spline-table paths build diagonal entries only.

Return type:

Propagators

Specification objects

Physical specification objects for the high-level System API.

Each object is a frozen dataclass that cleanly separates what the user wants to express (linear operator, interaction vertices, noise covariance) from how it is lowered to the raw package primitives (PropagatorModel, Field, Vertex, Action, PropagatorCache.precompute_C_table_*).

Design choices:

• Discriminated union via subclasses: e.g. LinearOp has DiagonalA, ExplicitR, … subclasses. Each subclass knows how to lower itself to a raw callable via build_R_callable().

• Escape hatches everywhere: any object can be bypassed by supplying a raw callable (ExplicitR, GeneralKappa2, etc.).

• Everything immutable: frozen dataclasses, safe to share between parallel workers and to hash for caching.

class sft_wick.workflow.specs.FieldSpec(name='phi', n_components=1)

Bases: object

Specification of the physical field φ of the theory.

The response field ψ is auto-generated under the hood and shares n_components; users should never need to reference it explicitly.

Parameters:

• name (str) – Physical field name (e.g. "phi"). Used only for symbolic rendering and as a convention.

• n_components (int) – Number of field components N. Use 1 for a scalar theory.

name: str = 'phi'

n_components: int = 1

property response_name: str

The conventional response-field name ("psi").

class sft_wick.workflow.specs.LinearOp

Bases: object

Base class — do not instantiate. See subclasses.

build_R_callable()

Return an R_time(t1, t2) -> float | (N, N) callable suitable for PropagatorModel.R_time.

Return type:

Callable

property is_iso_R: bool

True when R is scalar (no component structure).

class sft_wick.workflow.specs.DiagonalA(gamma, t_max_cache=100.0, n_grid_cache=200)

Bases: LinearOp

Diagonal linear operator A_{ab} = −γ_a δ_{ab}.

Two calling conventions for gamma:

Static (constant decay rates):

gamma is a length-N sequence of floats.

R_{aa}(t_1, t_2) = Θ(t_1 - t_2) · exp(-γ_a (t_1 - t_2))

Time-dependent (callable):

gamma is a callable γ(t) -> np.ndarray(shape=(N,)) returning the per-component instantaneous decay rate. In this case the wrapper pre-computes Γ_a(t) = ∫_0^t γ_a(τ) dτ on a time grid, caches it as a cubic spline, and evaluates R via:

R_{aa}(t_1, t_2) = Θ(t_1 - t_2) · exp(-(Γ_a(t_1) − Γ_a(t_2)))

The spline build cost is a one-time n_grid_cache calls to γ(t) + a cumulative trapezoidal integral. For the full-matrix (non-diagonal) time-dependent case — which requires a time-ordered matrix exponential — use ExplicitR with your own R callable.

Parameters:

• gamma (Any) – length-N sequence OR a callable γ(t) -> array(N).

• t_max_cache (float) – Upper bound of the Γ-spline grid, only used in the callable case. Queries beyond this bound extrapolate the spline (may be inaccurate — set this ≥ your maximum lambda_f).

• n_grid_cache (int) – Number of grid points for the Γ-spline build.

gamma: Any

t_max_cache: float = 100.0

n_grid_cache: int = 200

build_R_callable()

Return an R_time(t1, t2) -> float | (N, N) callable suitable for PropagatorModel.R_time.

Return type:

Callable

property is_iso_R: bool

True when R is scalar (no component structure).

class sft_wick.workflow.specs.ExplicitR(R_time, iso_R=True)

Bases: LinearOp

Escape hatch: user provides R directly, bypassing A.

Parameters:

• R_time (Callable) – (t1, t2) -> float | (N, N) callable. Must enforce causality (return 0 when t1 < t2).

• iso_R (bool) – True if the callable returns a scalar, False if matrix.

R_time: Callable

iso_R: bool = True

build_R_callable()

Return an R_time(t1, t2) -> float | (N, N) callable suitable for PropagatorModel.R_time.

Return type:

Callable

property is_iso_R: bool

True when R is scalar (no component structure).

class sft_wick.workflow.specs.ExponentialTemporal(lam, sigma_t)

Bases: object

κ²_t(Δt) = λ · exp(−|Δt|/σ_t) (OU kernel).

Parameters:

• lam (float)

• sigma_t (float)

lam: float

sigma_t: float

class sft_wick.workflow.specs.GaussianTemporal(lam, sigma_t)

Bases: object

κ²_t(Δt) = λ · exp(−Δt² / (2 σ_t²)) (Gaussian kernel).

Parameters:

• lam (float)

• sigma_t (float)

lam: float

sigma_t: float

class sft_wick.workflow.specs.ExponentialSpatial(sigma_x)

Bases: object

κ²_x(|Δx|) = exp(−|Δx|/σ_x) (exponential spatial envelope).

Parameters:

sigma_x (float)

sigma_x: float

class sft_wick.workflow.specs.GaussianSpatial(sigma_x)

Bases: object

κ²_x(|Δx|) = exp(−Δx² / (2 σ_x²)).

Parameters:

sigma_x (float)

sigma_x: float

class sft_wick.workflow.specs.LegendreAngular(coeffs)

Bases: object

Angular kernel on the sphere: κ²_x(cos θ) = Σ_ℓ C_ℓ P_ℓ(cos θ).

Parameters:

coeffs (Any) – [C_0, C_1, …, C_L] — ℓ-th entry is the Legendre coefficient for order ℓ.

coeffs: Any

class sft_wick.workflow.specs.CustomKernel(fn)

Bases: object

Escape hatch for a user-supplied 1-D kernel callable.

Parameters:

fn (Callable)

fn: Callable

class sft_wick.workflow.specs.Kappa2

Bases: object

Base class — do not instantiate. See subclasses.

build_callable(n_components)

Lower to a kappa2(n1, t1, n2, t2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

property homogeneity: str

One of 'translation' | 'rotation' | 'general'.

class sft_wick.workflow.specs.SeparableTranslation(temporal, spatial)

Bases: Kappa2

Translation-invariant + separable: κ²(n1, t1, n2, t2) = κ²_t(t1 − t2) · κ²_x(|n1 − n2|) · I_N.

Parameters:

• temporal (Callable) – callable accepting Δt, returning scalar.

• spatial (Callable) – callable accepting |Δx|, returning scalar.

temporal: Callable

spatial: Callable

build_callable(n_components)

Lower to a kappa2(n1, t1, n2, t2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

property homogeneity: str

One of 'translation' | 'rotation' | 'general'.

class sft_wick.workflow.specs.SeparableRotation(temporal, angular)

Bases: Kappa2

Rotation-invariant + separable: κ²(n1, t1, n2, t2) = κ²_t(t1 − t2) · κ²_Ω(cos θ) · I_N where cos θ = n1·n2/(|n1||n2|) and n1, n2 are direction vectors (typically on S²).

Parameters:

• temporal (Callable) – callable Δt → scalar.

• angular (Callable) – callable cos θ ∈ [−1, 1] → scalar.

temporal: Callable

angular: Callable

build_callable(n_components)

Lower to a kappa2(n1, t1, n2, t2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

property homogeneity: str

One of 'translation' | 'rotation' | 'general'.

class sft_wick.workflow.specs.GeneralKappa2(fn)

Bases: Kappa2

Escape hatch: user-supplied κ² callable without symmetry assumptions. Use this when κ² is not translation- or rotation-invariant and you want the full 4-D (t1, t2, x1, x2) spline build.

Parameters:

fn (Callable) – callable (n1, t1, n2, t2) -> (N, N).

fn: Callable

build_callable(n_components)

Lower to a kappa2(n1, t1, n2, t2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

property homogeneity: str

One of 'translation' | 'rotation' | 'general'.

class sft_wick.workflow.specs.Sigma2

Bases: object

Base — white-noise δ(t1 − t2)·σ²(t; n1, n2) component.

build_callable(n_components)

Lower to a sigma2(n1, t, n2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

class sft_wick.workflow.specs.ConstantImpulse(amplitude)

Bases: Sigma2

Time- and position-independent white-noise amplitude.

Parameters:

amplitude (Any) – scalar (isotropic) or (N, N) matrix.

amplitude: Any

build_callable(n_components)

Lower to a sigma2(n1, t, n2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

class sft_wick.workflow.specs.CustomImpulse(fn)

Bases: Sigma2

Escape hatch: arbitrary σ²(t; n1, n2) -> (N, N) callable.

Parameters:

fn (Callable)

fn: Callable

build_callable(n_components)

Lower to a sigma2(n1, t, n2) -> (N, N) callable.

Parameters:

n_components (int)

Return type:

Callable

class sft_wick.workflow.specs.GaussianNoise(kappa2, sigma2=None)

Bases: object

Gaussian driving: W_G = ½ ∫ψ ⊗ κ² ⊗ ψ with an optional white-noise impulse part.

Parameters:

• kappa2 (Kappa2) – the smooth (time-continuous) κ² cumulant.

• sigma2 (Sigma2 | None) – optional white-noise δ(t-t')·σ² part; None by default.

kappa2: Kappa2

sigma2: Sigma2 | None = None

class sft_wick.workflow.specs.LocalVertex(name, coupling)

Bases: object

A local MSR interaction vertex from the deterministic nonlinearity F^(n): one ψ leg and (n − 1) φ legs at a single spacetime point.

Parameters:

• name (str) – symbolic coupling name (e.g. "F"). Must be unique per system (used as a dict key in coupling_values).

• coupling (Any) – the bare F^(n) tensor — the coefficient as it appears in the deterministic equation of motion (dφ_a/dt = … + F^(n)_{a b_1 … b_{n−1}} φ_{b_1} … φ_{b_{n−1}} + …). Shape (N,)*n; the first axis is the ψ leg. The wrapper multiplies by the MSR factor -i internally (so demo1’s F_MSR = -1j * F_bare is automated).

Notes

Use msr_coupling to retrieve the MSR-factor-applied tensor — that is what is forwarded to the raw compute_moment / DiagramTerm.evaluate_coupling layer.

name: str

coupling: Any

property msr_coupling: ndarray

Bare F multiplied by the MSR factor -i.

class sft_wick.workflow.specs.NonLocalVertex(name, order, coupling, coupling_vectorized=False, equal_time=False, already_R_contracted=False)

Bases: object

A non-local MSR vertex from a higher driving-field cumulant κ^(m) (m ≥ 3): m ψ legs at m distinct spacetime points.

Parameters:

• name (str) – symbolic coupling name (e.g. "K").

• order (int) – m, the number of ψ legs (= rank of κ^(m)).

• coupling (Any) – the bare κ^(m) tensor — either a numeric tensor of shape (N,)*m (when κ^(m) is spacetime-independent) or a callable with signature fn(n_list, t_list) -> np.ndarray(shape=(N,)*m) where n_list / t_list are length-m sequences of the spatial / time coordinates. The wrapper multiplies by the MSR factor -(i^m) / m! internally (so demo2’s K = (i/6) * κ^(3) is automated).

• coupling_vectorized (bool) – only meaningful when coupling is a callable. False (default) signals the per-sample contract – the workflow calls fn with 1-D length-m n_list / t_list once per QMC sample. True signals the batched contract – the workflow calls fn with shape (m, n_samples) n_list / t_list once per integrand and expects an output of shape (n_samples,) + (N,)*m. Use the batched form when the user callable can amortise its cost across samples (numpy ufuncs, special functions, etc.); for cheap callables the per-sample form has lower overhead.

• equal_time (bool)

• already_R_contracted (bool)

Notes

The MSR factor values:

m	factor	simplified
1	−i	−i (mean drift)
2	−i²/2!	+½ (Gaussian kernel)
3	−i³/3!	+i/6 (demo2’s K)
4	−i⁴/4!	−1/24

Use msr_coupling to retrieve the MSR-factor-applied tensor or callable — that is what is forwarded to the raw coupling-values dict.

name: str

order: int

coupling: Any

coupling_vectorized: bool = False

equal_time: bool = False

already_R_contracted: bool = False

property msr_factor: complex

The −(i^m) / m! prefactor applied to the bare κ^(m).

property msr_coupling: Any

Bare κ^(m) multiplied by msr_factor.

If coupling is a callable, returns a wrapped callable that applies the factor at evaluation time — preserving the original signature.

Expansion

Expansion — diagram-level view of the perturbative expansion.

Everything a user wants to do between compute_moment and the final numeric result: inspect diagrams, classify by vertex composition, draw, render LaTeX, integrate point-by-point or as a sweep.

class sft_wick.workflow.expansion.Expansion(system, dts_by_order, orders, observable_repr, raw_result)

Bases: object

Result of System.expand(). Opaque to construct; use the accessors below.

Variables:

• system (Any) – the System this expansion was built from.

• dts_by_order (Mapping[int, list]) – {order: [DiagramTerm, ...]}.

• orders (tuple) – tuple of computed orders (sorted).

• observable_repr (tuple) – hashable description of the observable.

• raw_result (Any) – the underlying PerturbativeResult (advanced users).

Parameters:

• system (Any)

• dts_by_order (Mapping[int, list])

• orders (tuple)

• observable_repr (tuple)

• raw_result (Any)

system: Any

dts_by_order: Mapping[int, list]

orders: tuple

observable_repr: tuple

raw_result: Any

diagrams(order)

All DiagramTerm objects at this order.

Parameters:

order (int)

Return type:

list

summary()

Per-order count, plus a vertex_type histogram.

Return type:

dict[int, dict[str, int]]

by_vertex_type(order)

Group this order’s diagrams by vertex composition label.

Label format: the concatenation of the unique sorted coupling-symbol names appearing in the diagram’s coupling_sum. E.g.:

• A diagram using only the local F vertex → "F".

• Order-2 diagram mixing local F and non-local K (demo2 FK channel) → "FK".

• Pure-K diagram → "K".

The label is a set of vertex types, not a multiset — this matches demo2’s classification convention.

Parameters:

order (int)

Return type:

dict[str, list]

latex(order)

Concatenated LaTeX of every diagram at this order.

Parameters:

order (int)

Return type:

str

plot(order, i=0, **kwargs)

Render the i-th diagram at order via the package’s DiagramRenderer.

Parameters:

• order (int) – Perturbative order to look up in raw_result.

• i (int) – Index of the diagram within that order.

• **kwargs – Renderer-level kwargs (figsize, style, external_label_fn, vertex_label_fn, label_format) plus per-call draw kwargs (ax, title, external_labels, vertex_labels, positions, show_legend).

Returns:

The matplotlib Figure containing the rendered diagram.

evaluate(propagators, *, positions, t_final, component_pair=(0, 0), orders=None, vertex_types=None, integrate_over=None, method='qmc_vectorized', n_samples=8192, seed=42, n_jobs=1, n_gauss=8)

Integrate the expansion at a single (positions, t_final, component_pair) point. Returns a Result.

Parameters:

• propagators – Propagators from System.propagators().

• positions (dict[str, Any]) – {spatial_arg: x_value} mapping — e.g. {"x": 0.0, "y": 0.5}.

• t_final (float) – upper time bound for external-time integration (lambda_f).

• component_pair (tuple) – (a, b) component indices for the observable (e.g. (1, 1)). For scalar observables use (0, 0).

• orders (Iterable[int] | None) – subset of the expansion’s orders; None uses all.

• vertex_types (Iterable[str] | None) – subset of the vertex-composition labels (e.g. {"F"}, {"FK"}) to include — labels match by_vertex_type() keys. None ⇒ all. Useful for computing a single channel, or for skipping channels that require a bespoke integrator (e.g. non-local K whose coupling is spacetime-dependent; see demo2).

• integrate_over (Any) –

  Controls which external points have their time integrated.

  • None (default — physics observable): all externals held fixed at t_final. Matches the equal-time correlator ⟨φ(t_f) · φ(t_f)⟩ that is compared to MC data and demo notebooks.

  • "all": all externals integrated over [t_min, t_final] — the time-integrated moment ⟨∫φ(t)dt · ∫φ(t')dt'⟩. Natural e.g. for weak-lensing line-of-sight integrals.

  • Iterable of external-point names: mixed — those listed are integrated, others fixed. E.g. {"x"} for a source integrated along the line of sight × a detector field at t_final.

• method (str) –

  time-integrator selector. Recommended choice depends on the diagram’s number of internal time-integration variables d = len(time_integration_vars) and integrand smoothness:

  Method	Best for	Trade-off
  'qmc_vectorized' (default)	d >= 6 / non-smooth integrand	~ 1/sqrt(n_samples) bias
  'gauss_legendre'	d <= 5 smooth (the typical sft-wick case)	exponential convergence in n_gauss; cost n_gauss^d
  'nquad'	Adaptive 1-3D fallback	slow; raises NotImplementedError on dynamic-coupling
  'qmc' / 'qmc_scalar'	Compatibility / debugging	slow Python loop

  See Workflow API (L1 + L2) “Choosing an integrator” for the full decision matrix and worked examples.

• n_samples (int) – forwarded to the integrator (QMC only).

• seed (int | None) – forwarded to the integrator (QMC only).

• n_gauss (int) – nodes per dimension for method='gauss_legendre' (default 8 — exact for polynomials up to degree 15). Cost scales as n_gauss^d; bump to 12-20 for stiff integrands at large t_final.

• n_jobs (int)

sweep(propagators, *, positions_grid, t_final_grid, component_pairs=((0, 0),), orders=None, vertex_types=None, integrate_over=None, method='qmc_vectorized', n_samples=8192, seed=42, n_jobs=1, evaluate_n_jobs=1, n_gauss=8)

Cartesian-product sweep over positions, t_final, and component pairs.

Parameters:

• positions_grid (dict[str, list]) – {spatial_arg: [list of values]}. Each key’s list is swept independently; result is the full Cartesian product. E.g. {"x": [0.0], "y": [0.0, 0.5, 1.0]}.

• t_final_grid (list) – list of upper time bounds.

• component_pairs (Iterable[tuple]) – list of (a, b) component index tuples.

• vertex_types (Iterable[str] | None) – optional filter — same semantics as in evaluate(); only diagrams whose _vertex_type_label() lies in this set are integrated. None ⇒ all channels.

• n_jobs (int) – parallelise over Cartesian-product grid points (positions × t_final × component_pairs). 1 (default) preserves the original sequential behaviour — bit-identical when seed is fixed. -1 uses all CPU cores via joblib loky.

• evaluate_n_jobs (int) – parallelise over diagrams inside each grid point’s evaluate() call. Mutually exclusive with n_jobs > 1 (nested loky pools are not supported); the dispatcher raises if both are set. Use n_jobs > 1 when the sweep grid is large; use evaluate_n_jobs > 1 when each grid point has many diagrams (typical at orders >= 2).

• method (str) – integrator knobs – see evaluate() for the recommendation matrix and Workflow API (L1 + L2) “Choosing an integrator”. 'gauss_legendre' with n_gauss=8 is the right default for d ≤ 5 smooth integrands (exponential convergence, deterministic, no seed).

• n_samples (int) – integrator knobs – see evaluate() for the recommendation matrix and Workflow API (L1 + L2) “Choosing an integrator”. 'gauss_legendre' with n_gauss=8 is the right default for d ≤ 5 smooth integrands (exponential convergence, deterministic, no seed).

• seed (int | None) – integrator knobs – see evaluate() for the recommendation matrix and Workflow API (L1 + L2) “Choosing an integrator”. 'gauss_legendre' with n_gauss=8 is the right default for d ≤ 5 smooth integrands (exponential convergence, deterministic, no seed).

• n_gauss (int) – integrator knobs – see evaluate() for the recommendation matrix and Workflow API (L1 + L2) “Choosing an integrator”. 'gauss_legendre' with n_gauss=8 is the right default for d ≤ 5 smooth integrands (exponential convergence, deterministic, no seed).

• orders (Iterable[int] | None)

• integrate_over (Any)

Returns:

SweepResult with a pandas-friendly tidy table.

sft_wick.workflow.expansion.count_cross_group_c(dt)

Number of C propagators whose endpoints land in distinct direction groups (same helper as tests.test_deductive_numerics.TestSpatialAwareCache and examples.demo1.validate_phase5).

Return type:

int

Propagators

Propagators — thin wrapper around PropagatorCache that auto-dispatches precompute to the right homogeneity builder.

class sft_wick.workflow.propagators.Propagators(cache, homogeneity, is_lazy)

Bases: object

Holds a PropagatorCache preconfigured for a System. Opaque to the user — the only thing they do with this object is pass it to Expansion.evaluate() / Expansion.sweep().

Variables:

• cache (sft_wick.evaluate.PropagatorCache) – the underlying PropagatorCache.

• homogeneity (str) – resolved homogeneity string.

• is_lazy (bool) – whether the cache is in lazy-spline mode for its spatial dimension.

Parameters:

• cache (PropagatorCache)

• homogeneity (str)

• is_lazy (bool)

cache: PropagatorCache

homogeneity: str

is_lazy: bool

classmethod build(system, *, t_max, n_grid_t=60, homogeneity=None, r_max=None, n_grid_r=None, n_grid_cos=None, x_max=None, n_grid_x=None, n_jobs=1, c_closed_form=None, cache_path=None, interp_method='linear', c_closed_form_only=False, c_closed_form_vectorized=False, c_method='dblquad', c_n_gauss=20, diag_C=True)

Construct a Propagators for system. Called indirectly via System.propagators().

Parameters:

• c_closed_form (Callable | None) – optional fast path for C evaluation. If provided, it must be a callable (n1, t1, n2, t2) -> (N, N) returning the full C matrix at that spacetime-pair — the same signature as PropagatorCache._C_value_direct(). When set, the wrapper builds a PropagatorCache subclass that uses it instead of dblquad, collapsing the spline-table build time from minutes (typical for scipy.integrate.dblquad on fine grids) to milliseconds. Intended for kernels with known closed-form C (OU, separable exponentials).

• c_closed_form_only (bool) – when True (and c_closed_form is provided), skip the spline interpolator entirely. cache.C_at_batch then routes every lookup straight through the user’s c_fn – machine-precision agreement with the analytical C, no truncation error from grid spacing or spline order. Use this when the closed form is fast and the spline error would dominate over QMC noise (e.g. demo1’s OU kernel where sigma_t = 0.3 requires dt < 0.1 for sub-percent spline accuracy).

• c_closed_form_vectorized (bool) – only meaningful when c_closed_form_only=True. When True, the user’s c_fn must accept batched (t1, t2, x1, x2) arrays of shape (n,) and return a (n, N, N) tensor in a single call. When False, the cache falls back to a Python per-sample loop (slow; useful only for small point evaluations or when migrating an existing scalar c_fn).

• interp_method (str) – RegularGridInterpolator method used by full-grid C tables. 'linear' (default) is monotone and safe for steep cosmological tails; 'cubic' gives O(h⁴) accuracy on smooth, well-sampled grids. See PropagatorCache docstring for the full list of accepted methods and the linear-vs-cubic trade-off. Ignored under c_closed_form_only=True.

• c_method (str) –

  How the inner 2-D C-propagator integral ∫ R κ² R is evaluated when the cache builds its table.

  • 'dblquad' (default) – scipy.integrate.dblquad adaptive Gauss-Kronrod, robust on any κ² but slow (10-80 ms / call).

  • 'gauss_legendre' – tensor-product GL with a diagonal-aware sub-region split at λ1 = λ2. 18-100× faster on κ² that is piecewise analytic with at most a single |λ1−λ2| cusp on the diagonal (the standard exponential / Gaussian / OU family used in demo1, demo2, and the test suite). Returns near-machine-precision agreement with 'dblquad' at c_n_gauss=20.

  Ignored when c_closed_form_only=True (the C lookup bypasses the cache entirely in that mode).

• c_n_gauss (int) – Per-dimension GL node count for c_method='gauss_legendre' (default 20, enough for machine precision on smooth OU / Gaussian kernels). Cost scales as c_n_gauss² per sub-region.

• t_max (float)

• n_grid_t (int)

• homogeneity (str | None)

• r_max (float | None)

• n_grid_r (int | None)

• n_grid_cos (int | None)

• x_max (float | None)

• n_grid_x (int | None)

• n_jobs (int)

• cache_path (Any)

• diag_C (bool)

Return type:

Propagators

Result & SweepResult

Result and SweepResult — structured output from Expansion.evaluate() / Expansion.sweep().

These are what users actually look at. Accessors are designed around three common analyses:

1. Total: summed moment value.

2. By order: per-perturbative-order contribution. Reveals convergence.

3. By vertex type: demo2-style FF/FK/KK decomposition.

class sft_wick.workflow.result.Result(total, by_order, by_vertex_type, per_diagram, positions, t_final, component_pair, n_samples, seed)

Bases: object

Single-point evaluation result.

Variables:

• total (float) – summed moment value across all orders.

• by_order (dict) – {order: value}.

• by_vertex_type (dict) – {label: value} (e.g. {"F": 0.4, "FK": 0.01}).

• per_diagram (list) – list of dicts, one per diagram evaluated. Keys: order, diagram_idx, vertex_type, n_cross_C, value, error.

• positions (dict) – the {spatial_arg: x} at which the evaluation was performed.

• t_final (float) – the external-time upper bound used.

• component_pair (tuple) – observable’s component indices (e.g. (1, 1)).

• seed (n_samples,) – integrator settings (for reproducibility).

Parameters:

• total (float)

• by_order (dict)

• by_vertex_type (dict)

• per_diagram (list)

• positions (dict)

• t_final (float)

• component_pair (tuple)

• n_samples (int)

• seed (Any)

total: float

by_order: dict

by_vertex_type: dict

per_diagram: list

positions: dict

t_final: float

component_pair: tuple

n_samples: int

seed: Any

to_dict()

Flat dict suitable for serialisation or tabulation.

Return type:

dict

class sft_wick.workflow.result.SweepResult(rows, position_keys)

Bases: object

Grid-sweep result. Internally a tidy row store.

Each row is a dict with:

• the sweep coordinate fields (x, y, …, t_final, a, b),

• the diagram coordinate fields (order, diagram_idx, vertex_type, n_cross_C),

• the integrand output (value, error).

Typical usage:

sweep.to_dataframe().groupby(['r', 't_final']).value.sum()

Parameters:

• rows (list)

• position_keys (tuple)

rows: list

position_keys: tuple

to_dataframe()

Convert to a pandas.DataFrame. Requires pandas.

Columns include position keys, t_final, a, b, order, diagram_idx, vertex_type, n_cross_C, value, error.

totals()

Sum across diagrams at each sweep point, grouped by order and sweep coordinates. Returns a DataFrame with one row per (positions, t_final, a, b, order) and a value column.

Return type:

pd.DataFrame

by_vertex_type_totals()

Sum across diagrams grouped by (positions, t_final, a, b, vertex_type). Produces the demo2-style channel decomposition.

Return type:

pd.DataFrame

plot(*, x, y='value', hue='order', facet_col=None, aggregate='sum', **kwargs)

Quick line plot via matplotlib.

Parameters:

• x (str) – column to put on the x-axis (typically a position key).

• y (str) – value column (default 'value').

• hue (str | None) – column whose unique values become separate lines (default 'order').

• facet_col (str | None) – column whose unique values become subplot panels. None ⇒ single panel.

• aggregate (str) – 'sum' or 'none'. 'sum' aggregates y across all other (non-x/hue/facet) coordinates before plotting.

• kwargs – forwarded to matplotlib.axes.Axes.plot.

YAML + CLI (L2)

YAML-driven configuration for the workflow API.

Lets users declare an entire System / Expansion / Propagators / sweep pipeline in a single YAML file, with fields mapping 1:1 to the Python L1 API so the config is self-documenting.

Invoke via the CLI:

sft-wick run examples/demo1_config.yaml

Or programmatically:

from sft_wick.workflow.config import load_workflow_config, run_workflow

cfg = load_workflow_config("examples/demo1_config.yaml")
sweep, totals_df = run_workflow(cfg)

class sft_wick.workflow.config.WorkflowConfig(system, expand, propagators, sweep, output=<factory>)

Bases: object

Top-level parsed config.

Parameters:

• system (SystemConfig)

• expand (ExpandConfig)

• propagators (PropagatorsConfig)

• sweep (SweepConfig)

• output (list[OutputConfig])

system: SystemConfig

expand: ExpandConfig

propagators: PropagatorsConfig

sweep: SweepConfig

output: list[OutputConfig]

class sft_wick.workflow.config.SystemConfig(field_name: 'str', n_components: 'int', linear: 'dict', noise: 'dict', vertices: 'list', nonlocal_vertices: 'list' = <factory>, t_min: 'float' = 0.0, base_dir: 'Path | None' = None)

Bases: object

Parameters:

• field_name (str)

• n_components (int)

• linear (dict)

• noise (dict)

• vertices (list)

• nonlocal_vertices (list)

• t_min (float)

• base_dir (Path | None)

field_name: str

n_components: int

linear: dict

noise: dict

vertices: list

nonlocal_vertices: list

t_min: float = 0.0

base_dir: Path | None = None

class sft_wick.workflow.config.ExpandConfig(observable: 'tuple', orders: 'tuple', response_phase: 'bool' = True, ito: 'bool' = True, collect_topology: 'bool' = True, iso_R: 'Any' = None, diag_R: 'bool' = True, diag_C: 'bool' = True, iso_C: 'bool' = False, cache_path: 'Any' = None, n_jobs: 'int' = 1)

Bases: object

Parameters:

• observable (tuple)

• orders (tuple)

• response_phase (bool)

• ito (bool)

• collect_topology (bool)

• iso_R (Any)

• diag_R (bool)

• diag_C (bool)

• iso_C (bool)

• cache_path (Any)

• n_jobs (int)

observable: tuple

orders: tuple

response_phase: bool = True

ito: bool = True

collect_topology: bool = True

iso_R: Any = None

diag_R: bool = True

diag_C: bool = True

iso_C: bool = False

cache_path: Any = None

n_jobs: int = 1

class sft_wick.workflow.config.PropagatorsConfig(t_max: 'float', n_grid_t: 'int' = 60, dt: 'float | None' = None, homogeneity: 'Any' = None, r_max: 'Any' = None, n_grid_r: 'Any' = None, n_grid_cos: 'Any' = None, x_max: 'Any' = None, n_grid_x: 'Any' = None, n_jobs: 'int' = 1, c_closed_form_module: 'Any' = None, c_closed_form_attr: 'str' = 'C_fn', c_closed_form_only: 'bool' = False, c_closed_form_vectorized: 'bool' = False, cache_path: 'Any' = None, interp_method: 'str' = 'linear', c_method: 'str' = 'dblquad', c_n_gauss: 'int' = 20, diag_C: 'bool' = True)

Bases: object

Parameters:

• t_max (float)

• n_grid_t (int)

• dt (float | None)

• homogeneity (Any)

• r_max (Any)

• n_grid_r (Any)

• n_grid_cos (Any)

• x_max (Any)

• n_grid_x (Any)

• n_jobs (int)

• c_closed_form_module (Any)

• c_closed_form_attr (str)

• c_closed_form_only (bool)

• c_closed_form_vectorized (bool)

• cache_path (Any)

• interp_method (str)

• c_method (str)

• c_n_gauss (int)

• diag_C (bool)

t_max: float

n_grid_t: int = 60

dt: float | None = None

homogeneity: Any = None

r_max: Any = None

n_grid_r: Any = None

n_grid_cos: Any = None

x_max: Any = None

n_grid_x: Any = None

n_jobs: int = 1

c_closed_form_module: Any = None

c_closed_form_attr: str = 'C_fn'

c_closed_form_only: bool = False

c_closed_form_vectorized: bool = False

cache_path: Any = None

interp_method: str = 'linear'

c_method: str = 'dblquad'

c_n_gauss: int = 20

diag_C: bool = True

class sft_wick.workflow.config.SweepConfig(positions_grid: 'dict', t_final_grid: 'list', component_pairs: 'list', orders: 'Any' = None, vertex_types: 'Any' = None, integrate_over: 'Any' = None, method: 'str' = 'qmc_vectorized', n_samples: 'int' = 8192, seed: 'int' = 42, n_jobs: 'int' = 1, n_gauss: 'int' = 8)

Bases: object

Parameters:

• positions_grid (dict)

• t_final_grid (list)

• component_pairs (list)

• orders (Any)

• vertex_types (Any)

• integrate_over (Any)

• method (str)

• n_samples (int)

• seed (int)

• n_jobs (int)

• n_gauss (int)

positions_grid: dict

t_final_grid: list

component_pairs: list

orders: Any = None

vertex_types: Any = None

integrate_over: Any = None

method: str = 'qmc_vectorized'

n_samples: int = 8192

seed: int = 42

n_jobs: int = 1

n_gauss: int = 8

class sft_wick.workflow.config.OutputConfig(type: 'str', path: 'Any' = None, format: 'str' = 'markdown', x: 'Any' = None, y: 'str' = 'value', hue: 'Any' = 'order', facet_col: 'Any' = None)

Bases: object

Parameters:

• type (str)

• path (Any)

• format (str)

• x (Any)

• y (str)

• hue (Any)

• facet_col (Any)

type: str

path: Any = None

format: str = 'markdown'

x: Any = None

y: str = 'value'

hue: Any = 'order'

facet_col: Any = None

sft_wick.workflow.config.load_workflow_config(path, overrides=None)

Load and validate a workflow YAML config.

Parameters:

• path (str | Path) – path to a YAML file.

• overrides (dict | None) – optional {dotted.key: value} dict to patch the loaded config (e.g. {"sweep.seed": 7}).

Returns:

A WorkflowConfig ready to pass to run_workflow().

Return type:

WorkflowConfig

sft_wick.workflow.config.build_system(cfg)

Lower a parsed SystemConfig to a sft_wick.System instance.

Parameters:

cfg (SystemConfig)

sft_wick.workflow.config.run_workflow(cfg)

Execute the full pipeline — expand, build propagators, sweep, emit outputs.

Returns (sweep, totals_dataframe) for programmatic use.

Parameters:

cfg (WorkflowConfig)

sft-wick CLI — run a full workflow from a YAML config.

Usage:

sft-wick run CONFIG.yaml
    Execute expansion → propagators → sweep → emit outputs.

sft-wick run CONFIG.yaml --override sweep.seed=7 --override sweep.n_samples=4096
    Patch config fields without editing the file.

sft-wick run CONFIG.yaml --dry-run
    Parse + validate the config, print a resolved summary, exit 0.

sft_wick.workflow.cli.main(argv=None)

Parameters:

argv (list[str] | None)

Return type:

int