sft_wick.perturbation — Perturbative Expansion Driver

Perturbative expansion driver.

Computes <O>_S = sum_{n=0}^{N} (-1)^n / n! * <O * S_int^n>_{S_0} using Wick’s theorem for each order.

class sft_wick.perturbation.PerturbativeResult(order_terms, total, diagrams_by_order, diagram_terms_by_order=None)

Bases: object

Container for the result of a perturbative calculation.

Stores the symbolic expression at each perturbative order together with the corresponding Feynman diagrams.

Variables:

• order_terms – Mapping from perturbative order n to the simplified expression for that order.

• total – Sum of all non-zero order contributions.

• diagrams_by_order – Mapping from perturbative order n to the list of DiagramInfo records at that order.

Parameters:

• order_terms (dict[int, Expr])

• total (Expr)

• diagrams_by_order (dict[int, list[DiagramInfo]])

• diagram_terms_by_order (dict[int, list[DiagramTerm]] | None)

order(n)

Get the contribution at a specific perturbative order.

Parameters:

n (int) – The perturbative order (0, 1, 2, …).

Returns:

The simplified expression at order n, or ZERO if that order has no contribution.

Return type:

Expr

diagram_terms(order)

Structured diagram contributions for numerical evaluation.

Each DiagramTerm carries the propagators, coupling coefficient, prefactor, and index structure needed to evaluate a single Feynman diagram numerically.

Populated when collect_topology=True (the default).

Parameters:

order (int) – The perturbative order.

Returns:

List of DiagramTerm at that order, or empty list.

Return type:

list[DiagramTerm]

to_latex()

Generate LaTeX for the full result, order by order.

When diagram terms are available (the default collect_topology=True path), each order is rendered as a sum of fully-wrapped diagram contributions — with summation indices, integration variables, response phase, and prefactor all explicit.

Falls back to the raw symbolic expression when diagram terms are not populated.

Returns:

A multi-line string with one O(n): <latex> line per non-zero order.

Return type:

str

draw_diagrams(order=None, **kwargs)

Draw Feynman diagrams using matplotlib.

Topologically identical diagrams are drawn only once, with a ×N multiplicity label when N > 1.

Parameters:

• order (int | None) – If given, draw only diagrams at this perturbative order. Otherwise draw all diagrams.

• **kwargs – Forwarded to DiagramRenderer (e.g. figsize). Grid-level options accepted by draw_all() (ncols, suptitle, suptitle_kwargs, subtitle_fn, shared_legend, wspace, hspace, show, external_labels) are forwarded there instead.

Returns:

The matplotlib Figure containing the rendered diagrams, or None when there are no diagrams. Returning the figure makes the method display reliably as the final expression in notebooks and lets scripts call fig.savefig(...).

Return type:

Figure | None

class sft_wick.perturbation.DiagramInfo(observable_ops, vertex_instances, pairing, coefficient, order)

Bases: object

Lightweight record of a Feynman diagram for deferred construction.

The actual FeynmanDiagram graph is built lazily via to_feynman_diagram() to avoid up-front cost when many diagrams are generated.

Variables:

• observable_ops – Field operators forming the observable.

• vertex_instances – Instantiated vertices contributing to this diagram.

• pairing – The Wick contraction pairing (tuple of index pairs).

• coefficient – Rational prefactor for this diagram.

• order – The perturbative order at which this diagram appears.

Parameters:

• observable_ops (list[FieldOperator])

• vertex_instances (list[VertexInstance])

• pairing (Pairing)

• coefficient (Rational)

• order (int)

to_feynman_diagram()

Construct a FeynmanDiagram from this record.

Returns:

A fully-constructed FeynmanDiagram graph with external nodes, vertices, and propagator edges.

class sft_wick.perturbation.DiagramTerm(propagators, coupling_sum, rational_prefactor, integration_vars, summation_indices, n_response, equal_time_aliases=(), r_absorbed_pairs=())

Bases: object

A single Feynman diagram’s contribution, structured for numerical evaluation.

The full contribution is:

rational_prefactor × response_phase_factor × coupling_sum × ∏ propagators

summed over summation_indices and integrated over integration_vars.

Variables:

• propagators (tuple[sft_wick.expressions.Propagator, ...]) – Tuple of propagators forming the diagram.

• coupling_sum (sft_wick.expressions.Expr) – Symbolic coupling expression (sum of permuted couplings), without the rational prefactor.

• rational_prefactor (sft_wick.expressions.Rational) – The (-1)^n / n! × multinomial coefficient.

• integration_vars (tuple[str, ...]) – Spatial variables to integrate over.

• summation_indices (tuple[tuple[str, int], ...]) – (index_name, dimension) pairs for component index summations.

• n_response (int) – Number of R propagators (determines the response phase).

Parameters:

• propagators (tuple[Propagator, ...])

• coupling_sum (Expr)

• rational_prefactor (Rational)

• integration_vars (tuple[str, ...])

• summation_indices (tuple[tuple[str, int], ...])

• n_response (int)

• equal_time_aliases (tuple[tuple[str, str], ...])

• r_absorbed_pairs (tuple[tuple[str, str], ...])

propagators: tuple[Propagator, ...]

coupling_sum: Expr

rational_prefactor: Rational

integration_vars: tuple[str, ...]

summation_indices: tuple[tuple[str, int], ...]

n_response: int

equal_time_aliases: tuple[tuple[str, str], ...] = ()

r_absorbed_pairs: tuple[tuple[str, str], ...] = ()

property propagator_indices: tuple[tuple[str, int], ...]

Summation indices that appear in at least one propagator.

property coupling_only_indices: tuple[tuple[str, int], ...]

Summation indices that appear only in the coupling sum, not in any propagator.

spatial_topology()

Return (kind, spatial_left, spatial_right) for each propagator.

Return type:

list[tuple[str, str, str]]

response_phase_factor()

Return (-i)^n_response as a complex number.

Return type:

complex

evaluate_coupling(coupling_values, fixed_indices=None)

Substitute numeric coupling tensor values and return an array.

The output is indexed by propagator indices only. For each combination of propagator index values the coupling sum is evaluated by summing over coupling-only indices (indices that appear in the coupling expression but in no propagator). This matches the natural contraction order: fix the propagator legs, then sum the vertex factors internally.

Parameters:

• coupling_values (dict) – {name: array} mapping coupling names to NumPy arrays. For a rank-3 coupling F, the array shape should be (N, N, N) where N is the number of field components.

• fixed_indices (dict[str, int] | None) – Optional {index_name: int_value} for indices pinned by external constraints (e.g. observable component indices like {'1': 0}).

Returns:

NumPy array with shape (dim_i0, dim_i1, ...) for the propagator indices, with rational_prefactor and the MSR response phase (-i)^n_response already applied. Returns a 0-d array when there are no propagator indices. The result is complex whenever n_response % 4 in {1, 3} or the coupling tensors are complex-valued.

Return type:

ndarray

evaluate_coupling_batched(coupling_values_batched, n_samples, fixed_indices=None)

Vectorised counterpart of evaluate_coupling().

Same contract as evaluate_coupling(), but every per-sample Symbol value carries a leading sample axis of length n_samples. Returns an array whose shape is (n_samples,) + prop_shape (or just (n_samples,) when there are no propagator indices).

Parameters:

• coupling_values_batched (dict) – {name: ndarray} – arrays may be either (n_samples, *kappa_shape) (dynamic / per-sample) or (*kappa_shape,) (static, broadcast across the sample axis).

• n_samples (int) – Length of the sample axis.

• fixed_indices (dict[str, int] | None) – Optional {index_name: int_value} for indices pinned by external constraints (e.g. observable component indices like {'a': 0}).

Returns:

Complex / float array. When this DiagramTerm has propagator indices, the returned shape is (n_samples,) + prop_shape. Otherwise it is (n_samples,).

Raises:

NotImplementedError – If the symbolic coupling_sum contains a node type the batched evaluator cannot handle. Callers MUST catch this and fall back to a per-sample loop over evaluate_coupling().

Return type:

ndarray

apply_diagonal(*, diag_R=False, diag_C=False, iso_R=False, iso_C=False)

Return a new term with diagonal (and optionally isotropic) propagator constraints applied.

Eliminates summation indices that are pinned by diagonal propagators and substitutes index equalities into the coupling expression. The delta constraint collapses the double sum into a single sum, so no extra factor is applied.

With iso_R=True (implies diag_R=True), the remaining equal component indices are stripped from R propagators entirely, expressing the assumption that all diagonal R entries are the same scalar R.

Parameters:

• diag_R (bool) – Enforce diagonal response propagators.

• diag_C (bool) – Enforce diagonal correlation propagators.

• iso_R (bool) – Enforce isotropic (index-free) response propagators.

• iso_C (bool) – Enforce isotropic (index-free) correlation propagators.

Returns:

A new DiagramTerm with reduced summation indices.

Return type:

DiagramTerm

to_latex()

Full LaTeX representation including summations, integrals, and phase.

Renders the complete contribution:

\sum_{i_0} ... \int dy_0 ... coeff × (coupling) × propagators

where coeff combines the rational prefactor and the MSR phase (-i)^{n_R} into a single coefficient.

Return type:

str

to_latex_evaluated(coupling_values, fixed_indices=None)

LaTeX with numerically evaluated coupling coefficients.

Like to_latex(), but replaces the symbolic coupling sum with the numerical coefficient array obtained from evaluate_coupling(). Summation indices that survive in the propagators are shown as explicit sums over the numerical coefficients; coupling-only indices are already contracted.

Parameters:

• coupling_values (dict) – {name: numpy_array} for coupling tensors.

• fixed_indices (dict[str, int] | None) – Optional pinned observable indices.

Returns:

LaTeX string with numerical coefficients.

Return type:

str

analyze_spatial()

Analyze R-propagator connectivity and time orderings.

Returns a SpatialStructure describing direction identification groups, causal time orderings, and surviving integration variables.

Return type:

SpatialStructure

build_integrand(coupling_values, fixed_indices=None)

Build a DiagramIntegrand for numerical evaluation.

Combines coupling coefficient evaluation (Step 1) with spatial structure analysis (Step 2) into a ready-to-evaluate object.

Parameters:

• coupling_values (dict) –

  dict mapping coupling-symbol name → value. Each value is either:

  • a numeric numpy.ndarray of shape (N,)*rank (e.g. a spacetime-independent local-vertex coefficient F, or a placeholder for a non-local tensor whose spacetime dependence is trivial),

  • a callable fn(n_list, t_list) → ndarray for spacetime-dependent couplings (e.g. demo2’s κ^{(3)}(x₁,t₁; x₂,t₂; x₃,t₃)). Here n_list and t_list are length-order sequences of the vertex’s ψ-leg positions and times.

  When any value is callable the integrand enters a per-QMC-sample evaluation path (see DynamicCouplingPromise); all-ndarray values keep the fast vectorised path.

• fixed_indices (dict[str, int] | None) – Optional pinned index values for observable component labels (e.g. {'a': 1, 'b': 1}).

Returns:

A DiagramIntegrand that can evaluate the integrand at specific time/direction coordinates.

Return type:

DiagramIntegrand

sft_wick.perturbation.compute_moment(observable, action, order, ito=True, response_phase=True, collect_topology=True, diag_R=False, diag_C=False, iso_R=False, iso_C=False)

Compute the perturbative expansion of an observable.

Evaluates

\[\langle \mathcal{O} \rangle_S = \sum_{n=0}^{N} \frac{(-1)^n}{n!}\, \langle \mathcal{O}\, S_{\mathrm{int}}^{\,n} \rangle_{S_0}\]

up to the requested perturbative order. In the MSR formalism the partition function \(Z = 1\), so there is no denominator.

Parameters:

• observable (list[FieldOperator]) – List of field operators defining the observable \(\mathcal{O}\).

• action (Action) – The interaction action \(S_{\mathrm{int}}\).

• order (int) – Maximum perturbative order N to compute.

• ito (bool) – If True, apply the Itô prescription \(\Theta(0)=0\): the response propagator vanishes at equal spatial points, \(R(x,x)=0\), and causal R-loops are eliminated.

• response_phase (bool) – If True, multiply each term by \((-\mathrm{i})^n\) where n is the number of response propagators in that term, implementing the convention \(\langle\phi\,\psi\rangle = -\mathrm{i}\,R\).

• collect_topology (bool) – If True, group terms that share the same propagator spatial structure and factor out the propagators, summing the coupling coefficients with appropriately permuted indices.

• diag_R (bool) – If True, enforce diagonal response propagators \(R_{ij} = \delta_{ij} R\), eliminating one summation index per R propagator.

• diag_C (bool) – If True, enforce diagonal correlation propagators \(C_{ij} = \delta_{ij} C\), eliminating one summation index per C propagator.

• iso_R (bool) – If True (implies diag_R), further assume all diagonal R entries are equal, \(R_{ii} = R\). The component index is dropped from R propagators entirely.

• iso_C (bool) – If True (implies diag_C), the same for C.

Returns:

A PerturbativeResult containing order-by-order expressions, a combined total, and Feynman diagram information.

Return type:

PerturbativeResult

sft_wick.perturbation.compute_moment_numerical(observable, action, order, coupling_values, fixed_indices=None, ito=True, diag_R=False, diag_C=False, iso_R=False, iso_C=False, n_jobs=1)

Compute diagram terms using nauty graph isomorphism.

A faster alternative to compute_moment() that replaces the \(k!\) brute-force canonical form search with the nauty graph isomorphism algorithm (via pynauty), enabling order-6 calculations that were previously infeasible.

Optionally parallelizes the nauty canonicalization and component routing steps using joblib (install with pip install sft-wick[parallel]).

Parameters:

• observable (list[FieldOperator]) – Field operators defining the observable.

• action (Action) – The interaction action.

• order (int) – Maximum perturbative order.

• coupling_values (dict) – {name: numpy_array} mapping coupling names to numeric tensors.

• fixed_indices (dict[str, int] | None) – {index_label: int_value} for observable component indices (e.g. {'a': 1, 'b': 1}).

• ito (bool) – Apply Itô prescription (default True).

• diag_R (bool) – Enforce diagonal R propagators.

• diag_C (bool) – Enforce diagonal C propagators.

• iso_R (bool) – Isotropic R (implies diag_R).

• iso_C (bool) – Isotropic C (implies diag_C).

• n_jobs (int) – Number of parallel workers for nauty canonicalization and component routing. 1 = sequential (default), -1 = use all CPU cores. Requires joblib.

Returns:

{order_n: [DiagramTerm, ...]} for each order with non-zero contributions. Use DiagramTerm.build_integrand() for numerical evaluation.

Return type:

dict[int, list[DiagramTerm]]