sft_wick.evaluate — Numerical Evaluation Pipeline

Numerical evaluation pipeline for DiagramTerm objects.

Provides the 4-step workflow: 1. Coupling coefficients — handled by DiagramTerm.evaluate_coupling() (existing) 2. Spatial structure analysis — analyze_spatial() 3. Propagator evaluation — PropagatorModel + PropagatorCache 4. Contraction & integration — DiagramIntegrand

class sft_wick.evaluate.SpatialStructure(direction_groups, direction_map, time_orderings, r_propagators, c_propagators, time_integration_vars, direction_integration_vars, external_points, equal_time_aliases=(), r_absorbed_pairs=())[source]

Bases: object

Result of analyzing R-propagator connectivity and time orderings.

R propagators of the form δ(n−n') Θ(t−t') R_time(t,t') identify directions along R-chains and impose causal time orderings.

Parameters:
direction_groups: tuple[frozenset[str], ...]

Groups of spatial points connected by R propagators (share direction).

direction_map: dict[str, str]

spatial_point → representative direction variable name (e.g. 'n_x').

time_orderings: tuple[tuple[str, str], ...]

(earlier_point, later_point) from R causality. R(a, b) with Θ(t_a − t_b) means t_b ≤ t_a.

Type:

Time ordering pairs

r_propagators: tuple[tuple[str, str], ...]

R propagators as (spatial_left, spatial_right) pairs.

c_propagators: tuple[tuple[str, str, str | None, str | None], ...]

(spatial_left, spatial_right, index_left, index_right).

Type:

C propagators

time_integration_vars: tuple[str, ...]

Time integration variables, topologically sorted (innermost first). For an equal_time non-local vertex, only the canonical representative leg appears here; the other (m-1) legs of that vertex are removed and their time values are read through the equal_time_aliases map below.

direction_integration_vars: tuple[str, ...]

Surviving direction integration variables (one per R-component that contains only integration points — typically empty).

external_points: tuple[str, ...]

External (non-integration) spatial points.

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

Maps non-representative internal spatial labels to their canonical time representative, propagated from DiagramTerm.equal_time_aliases. The integration loops look up every leg’s time via equal_time_aliases.get(label, label) so aliased legs share a single integration variable while keeping independent spatial labels. Empty tuple ⇒ no aliasing (original cross-spacetime cumulant behaviour).

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

(partner_label, leg_label) pairs identifying R-propagators whose factor has been absorbed into an upstream NonLocalVertex(already_R_contracted=True) callable. Propagated straight from DiagramTerm.r_absorbed_pairs. The integrand R-product loop iterates over this set to skip the matching propagators (so the factor becomes 1 rather than R(t, t) = 0 under Itô after the leg time aliases onto the partner’s). The leg’s time / direction collapse onto the partner’s via the accompanying equal_time_aliases entries.

sft_wick.evaluate.analyze_spatial(dt)[source]

Analyze a DiagramTerm’s propagator topology.

Determines direction identification groups from R-propagator connectivity, time ordering constraints from R causality, and which integration variables survive after δ-function elimination.

Parameters:

dt (DiagramTerm)

Return type:

SpatialStructure

class sft_wick.evaluate.PropagatorModel(R_time, kappa2, n_components, iso_R=True, diag_C=False, t_min=0.0, sigma2=None)[source]

Bases: object

Physical propagator functions for numerical evaluation.

The R propagator has the form:

R_{ab}(n,t; n',t') = δ(n−n') Θ(t−t') R_time(t,t') δ_{ab}   [iso_R]
R_{ab}(n,t; n',t') = δ(n−n') Θ(t−t') R_time_ab(t,t')        [general]

The C propagator is computed from R and the 2-point cumulant κ^(2):

C_{ab}(n₁,t₁; n₂,t₂) = ∫∫ R_time(t₁,λ') κ_{ab}(n₁,λ'; n₂,λ'')
                              R_time(t₂,λ'') dλ' dλ''
Variables:

  • R_time (Callable) – Callable (t_left, t_right) → scalar when iso_R=True, or (t_left, t_right) → (N, N) array otherwise.

  • kappa2 (Callable) – Callable (n1, t1, n2, t2) → (N, N) array. The 2-point cumulant.

  • n_components (int) – Number of field components.

  • iso_R (bool) – Whether R is isotropic (proportional to identity).

  • diag_C (bool) – Whether C is diagonal in component indices.

  • t_min (float) – Lower bound for the λ integrals in C computation.

Parameters:
R_time: Callable
kappa2: Callable
n_components: int
iso_R: bool = True
diag_C: bool = False
t_min: float = 0.0
sigma2: Callable | None = None
class sft_wick.evaluate.PropagatorCache(model, quad_opts=None, homogeneity='translation', interp_method='linear', c_method='dblquad', n_gauss=20)[source]

Bases: object

Evaluates and caches propagator values.

Computes C via double integration of R · κ · R and caches results to avoid redundant quadrature.

Parameters:

  • model (PropagatorModel) – The physical propagator model.

  • quad_opts (dict | None) – Options passed to scipy.integrate.dblquad (e.g. {'epsabs': 1e-8, 'epsrel': 1e-8}).

  • homogeneity (str)

  • interp_method (str)

  • c_method (str)

  • n_gauss (int)

R_time(t_left, t_right)[source]

Evaluate R_time(t_left, t_right).

Returns scalar if model.iso_R=True, else (N, N) array. The Θ function is NOT enforced here — the integration domain handles causality.

Parameters:
Return type:

float | ndarray

R_product(r_pairs, times)[source]

Product of R_time over all R propagator pairs.

Only valid when model.iso_R=True (R is scalar).

Parameters:
Return type:

float

C_value(n1, t1, n2, t2)[source]

Compute the full C matrix C_{ab}(n1, t1; n2, t2).

If precompute_C_table() has been called and t1, t2 are within the table range, uses fast spline interpolation. Otherwise falls back to scipy.integrate.dblquad:

C_{ab} = ∫_{t_min}^{t1} dλ' ∫_{t_min}^{t2} dλ''
         R_time(t1,λ') κ_{ab}(n1,λ'; n2,λ'') R_time(t2,λ'')
Returns:

(N, N) array.

Parameters:
Return type:

ndarray

C_diagonal(n, t1, n_prime=None, t2=None)[source]

Return diagonal [C_{00}, C_{11}, …] as a 1D array.

If n_prime and t2 are None, uses equal-point C(n,t; n,t). Uses spline table when available for fast lookup.

Parameters:
Return type:

ndarray

clear_cache()[source]

Clear the C value cache and spline table.

Return type:

None

precompute_C_table(t_max, n_grid=100, direction=0)[source]

Pre-compute C on a grid and build spline interpolators.

For iso_R + diag_C with direction-independent kappa2, C_{aa}(t1, t2) depends only on (t1, t2). This method evaluates C on an n_grid × n_grid grid via dblquad, then builds a RectBivariateSpline for each diagonal component.

After calling this, C_value() and C_diagonal() use fast spline interpolation instead of dblquad.

Parameters:

  • t_max (float) – Upper bound of the time grid.

  • n_grid (int) – Number of grid points per axis (default 100).

  • direction (Any) – Direction value to use for kappa2 evaluation.

Return type:

None

property has_table: bool

Whether a pre-computed C table is available.

C_diagonal_batch(t1, t2)[source]

Evaluate diagonal C at arrays of time pairs.

Requires precompute_C_table() to have been called.

Parameters:

  • t1 (ndarray) – Array of shape (n,) — first time coordinates.

  • t2 (ndarray) – Array of shape (n,) — second time coordinates.

Returns:

Array of shape (n, N) where N is the number of field components. result[s, a] is C_{aa}(t1[s], t2[s]).

Return type:

ndarray

R_time_batch(t1, t2)[source]

Evaluate R_time at arrays of time pairs (iso_R only).

Parameters:

  • t1 (ndarray) – Array of shape (n,) — left times.

  • t2 (ndarray) – Array of shape (n,) — right times.

Returns:

Array of shape (n,). R(t1, t2) = Θ(t1−t2) × R_time(t1, t2).

Return type:

ndarray

precompute_C_table_translation(t_max, n_grid_t=60, r_max=None, n_grid_r=None, n_jobs=1, c_method='dblquad', n_gauss=20)[source]

Enable spatial support under the translation-invariance assumption:

C(t1, t2; x1, x2) depends only on |x1 − x2|.

Two operating modes depending on whether r_max and n_grid_r are supplied:

Lazy mode (default, recommended for fixed-point moments)

r_max=None and n_grid_r=None. No r-grid is pre-computed. When C_at_batch() encounters a new r = |x1−x2| value it builds and caches a 2-D (t1, t2) spline on demand. For a moment calculation at fixed external points (x_1, …, x_n) there are at most O(n²) distinct r-values, so lazy mode saves an enormous amount of work vs a densely-sampled r-grid — particularly when the underlying _C_value_direct is scipy.integrate.dblquad (expensive).

Full-grid mode (right for sweeping r curves)

Both r_max and n_grid_r provided. A 3-D (t1, t2, r) spline is built up-front over the full range, giving O(1) evaluation per query — faster once n_distinct_r >> n_grid_r. Appropriate when the user will plot C or a derived observable as a function of r over a continuous range.

Parameters:

  • t_max (float) – upper bound of the (t1, t2) time grid.

  • n_grid_t (int) – grid size along each time axis.

  • r_max (float | None) – upper bound of the r = |x1-x2| grid (full-grid mode). None ⇒ lazy mode.

  • n_grid_r (int | None) – grid size along r (full-grid mode). None ⇒ lazy mode.

  • n_jobs (int) – parallel workers for grid-point _C_value_direct evaluations. 1 serial, -1 all cores via joblib (loky backend). Applies in both full-grid mode (fans out across the 3-D (t1, t2, r) grid) and lazy mode (fans out across each on-demand 2-D (t1, t2) build; each new r pays a ~1 s worker startup cost).

  • c_method (str) –

    Quadrature method passed to _C_value_direct() for every grid-point evaluation (full-grid mode only – lazy mode will use the cache instance’s c_method default).

    • 'dblquad' (default): backwards-compatible adaptive quadrature.

    • 'gauss_legendre': tensor-product GL with diagonal-aware sub-region splitting. 10-1000× faster on smooth κ² with a single |λ1−λ2| cusp. See _C_value_direct_gl().

  • n_gauss (int) – Per-dimension GL node count when c_method='gauss_legendre'. Default 20 is enough for machine precision on demo1/demo2 OU kernels.

Return type:

None

Requires self.homogeneity == 'translation'.

Side effect: sets either _c_translation_splines (full) or _lazy_translation (lazy).

precompute_C_table_rotation(t_max, n_grid_t=60, n_grid_cos=None, n_jobs=1, c_method='dblquad', n_gauss=20)[source]

Enable spatial support under the rotation-invariance assumption:

C(t1, t2; x1, x2) depends only on x1 · x2.

Appropriate when x is a unit direction vector (e.g. on the sphere). The effective spatial parameter is cos θ ∈ [−1, 1].

Two modes, analogous to precompute_C_table_translation():

  • Lazy (n_grid_cos=None): 2-D (t1, t2) splines built on demand per distinct cos θ value. Best when only a few cos θ values arise (e.g. a fixed-set moment calculation).

  • Full-grid (n_grid_cos integer): 3-D (t1, t2, cos) spline over the whole cos ∈ [−1, 1] range. Best for angular-correlation sweeps.

Under the hood, the reference evaluation at a specific cos θ uses two representative unit vectors e_1 = (1, 0, …) and e_2 chosen so that e_1 · e_2 = cos θ (a 2-D rotation of e_1).

Parameters:

  • n_jobs (int) – parallel workers for grid-point _C_value_direct evaluations; applies in both full-grid (fans out across the 3-D (t1, t2, cos) grid) and lazy mode (fans out across each on-demand 2-D build). 1 serial, -1 all cores via joblib.

  • c_method (str) – 'dblquad' (default) or 'gauss_legendre' – forwarded to _C_value_direct() for every full-grid point. See precompute_C_table_translation() for details.

  • n_gauss (int) – Per-dim GL node count for c_method='gauss_legendre'; default 20.

  • t_max (float)

  • n_grid_t (int)

  • n_grid_cos (int | None)

Return type:

None

Requires self.homogeneity == 'rotation'.

precompute_C_table_general(t_max, n_grid_t=40, x_max=None, n_grid_x=None, n_jobs=1, c_method='dblquad', n_gauss=20)[source]

Pre-compute C with no spatial symmetry assumed.

Two modes:

  • Lazy (x_max=None or n_grid_x=None): 2-D (t1, t2) splines built on demand per distinct (x1, x2) pair. Best for fixed-point moments where only a handful of pairs occur.

  • Full-grid: 4-D (t1, t2, x1, x2) interpolator over the pre-allocated grid. Build cost is n_grid_t² · n_grid_x² independent _C_value_direct calls — strongly consider n_jobs=-1.

Parameters:

  • n_jobs (int) – parallel workers for _C_value_direct evaluations; applies in both full-grid (fans out across the 4-D grid) and lazy mode (fans out across each on-demand 2-D build). 1 serial, -1 all cores via joblib.

  • c_method (str) – 'dblquad' (default) or 'gauss_legendre' – forwarded to _C_value_direct() for every full-grid point. See precompute_C_table_translation() for details.

  • n_gauss (int) – Per-dim GL node count for c_method='gauss_legendre'; default 20.

  • t_max (float)

  • n_grid_t (int)

  • x_max (float | None)

  • n_grid_x (int | None)

Return type:

None

Requires self.homogeneity == 'general'.

C_at_batch(t1, t2, x1, x2)[source]

Batch-evaluate C at arbitrary (t1, t2, x1, x2) arrays.

Dispatches on homogeneity:

  • 'translation': C depends only on |x1 − x2|. Uses the 3-D (t1, t2, r) full-grid spline if built via precompute_C_table_translation(r_max=..., n_grid_r=...); otherwise uses a lazy per-r 2-D spline cache; otherwise (legacy 2-D spline only) falls back to C_diagonal_batch() — x is ignored in that case.

  • 'rotation': C depends only on x1 · x2 / (|x1| |x2|). Uses the 3-D (t1, t2, cos) spline or per-cos lazy cache.

  • 'general': 4-D (t1, t2, x1, x2) spline, or a lazy per-pair 2-D cache.

Parameters:

  • t1 (ndarray) – Time arrays of shape (n,).

  • t2 (ndarray) – Time arrays of shape (n,).

  • x1 (ndarray | float) – Spatial coordinates. Scalar or array of shape (n,) (broadcast if scalar). # TODO(d-dim): support (n, d).

  • x2 (ndarray | float) – Spatial coordinates. Scalar or array of shape (n,) (broadcast if scalar). # TODO(d-dim): support (n, d).

Returns:

Array of shape (n, N) — diagonal C values per component.

Return type:

ndarray

class sft_wick.evaluate.DynamicCouplingPromise(diagram_term, static_values, dynamic_values, spatial_args_by_name, fixed_indices)[source]

Bases: object

Defers coupling-tensor materialisation to QMC-sample time.

When any entry in coupling_values is a callable (e.g. a spacetime-dependent non-local vertex such as demo2’s κ^{(3)}(x₁,t₁; x₂,t₂; x₃,t₃)), the usual pre-QMC DiagramTerm.evaluate_coupling() path doesn’t apply — the callable’s output changes with each sample’s QMC time coordinates. This class packages everything needed for the per-sample path:

  • the list of static (ndarray) symbol values,

  • the list of dynamic (callable) symbol values,

  • for each dynamic symbol, the ψ-leg spatial-label tuple extracted from the DiagramTerm’s coupling_sum so we can look up each leg’s time and position per sample.

The per-sample evaluator evaluate_at() materialises the dynamic tensors using the sample’s (times, positions) and then delegates to the static-path DiagramTerm.evaluate_coupling() with a fully-numeric coupling_values dict — so the contraction code stays the same.

Parameters:

  • diagram_term (Any)

  • static_values (dict)

  • dynamic_values (dict)

  • spatial_args_by_name (dict)

  • fixed_indices (dict)

diagram_term: Any

The parent DiagramTerm.

static_values: dict

Static coupling values — already materialised arrays.

dynamic_values: dict

Dynamic coupling values — mapping name -> callable(n_list, t_list) returning an ndarray.

spatial_args_by_name: dict

Per-dynamic-symbol tuple of ψ-leg spatial labels, as they appear in the diagram’s coupling_sum.

fixed_indices: dict

Component-index pins forwarded from build_integrand.

evaluate_at(times, positions)[source]

Materialise the dynamic callables at this sample’s (times, positions) and return the per-sample coupling_array.

Parameters:

  • times (dict) – {spatial_label: time_value} for all spatial labels referenced by any dynamic symbol’s legs.

  • positions (dict) – {spatial_label: position_value} same.

Return type:

ndarray

evaluate_at_batch(label_t, label_x, n_samples)[source]

Vectorised batch evaluator – returns a (n_samples,) complex array of the contracted scalar coupling per sample.

label_t and label_x map every spatial label appearing in any dynamic symbol’s legs to a (n_samples,) time array / position respectively. Each label_x entry may be:

  • a scalar (the historical / 1-D translation case) — produces a per-leg broadcast of shape (n_samples,);

  • a (d,) vector (e.g. a 3-D unit vector under homogeneity='rotation') — produces a per-leg broadcast of shape (n_samples, d).

The user callable then sees per-leg slices of shape (m,) or (m, d) respectively (or, in the vectorised contract below, (m, n_samples) / (m, n_samples, d)). All legs of one symbol must have the same position shape; mixing scalar and vector legs raises from np.stack.

Per-symbol fast path: when the user marks a callable vectorized=True (e.g. via NonLocalVertex(coupling_vectorized=True) or by setting fn.vectorized = True), the wrapped callable receives (m_legs, n_samples) arrays — or (m_legs, n_samples, d) for d-dim positions — in a single call and returns a tensor of shape (n_samples,) + (N,)*order. Otherwise we fall back to the n_samples per-sample calls of evaluate_at().

Either way, we still loop in Python over n_samples to contract each sample’s tensor down to a scalar via DiagramTerm.evaluate_coupling(). Vectorising the contraction itself across a sample axis is the deferred prop-indexed-dynamic refactor and is not done here.

Raises NotImplementedError if any sample’s contracted coupling is not scalar (the prop-indexed case locked by DC1 in tests/test_dynamic_coupling.py).

Parameters:
Return type:

ndarray

class sft_wick.evaluate.DiagramIntegrand(diagram_term, spatial, coupling_array, fixed_indices=<factory>, dynamic_coupling=None)[source]

Bases: object

A diagram’s contribution decomposed for numerical integration.

Combines coupling coefficients (Step 1) with spatial analysis (Step 2) for efficient numerical evaluation.

The integrand has the factored form (when iso_R=True):

R_product(times) × Σ_{prop_idx} coeff[idx] × Π_C C_values[idx]
Parameters:
diagram_term: DiagramTerm

The underlying DiagramTerm.

spatial: SpatialStructure

Spatial analysis result.

coupling_array: np.ndarray

Coupling coefficient array from evaluate_coupling(). When dynamic_coupling is set, this is a placeholder zero array and the real coupling value is computed per-sample in the integrator.

fixed_indices: dict[str, int]

Fixed component indices (e.g. {'a': 0, 'b': 1}). Used by the QMC/GL integrators to resolve C-propagator component indices that are not summation variables.

dynamic_coupling: Any = None

Optional DynamicCouplingPromise carrying per-sample coupling materialisation logic for spacetime-dependent (callable) coupling values. None ⇒ the fully-static fast path is used (all coupling values were ndarrays).

evaluate(times, directions, cache)[source]

Evaluate the integrand at specific time + direction coordinates.

Parameters:

  • times (dict[str, float]) – {spatial_point: time_value} for ALL spatial points. Aliased (non-representative) legs of equal-time non-local vertices may be absent — their values are filled in from the canonical representatives inside this method.

  • directions (dict[str, Any]) – {direction_var: value} for independent direction variables (one per R-connected component, keyed by the representative name from spatial.direction_map).

  • cache (PropagatorCache) – A PropagatorCache for propagator evaluation.

Returns:

Scalar value of the integrand (complex if coupling is complex).

Return type:

complex

make_scipy_integrand(external_times, external_directions, cache)[source]

Return a callable for scipy.integrate.nquad.

The returned function accepts the integration time variables as positional arguments in the order of spatial.time_integration_vars and returns the integrand value. External times and directions are baked in.

Parameters:
Returns:

Callable f(*time_args) → float.

Return type:

Callable

integration_bounds(external_times, t_min=0.0)[source]

Return integration bounds for scipy.integrate.nquad.

Time ordering constraints from R causality translate to variable-dependent bounds. For nquad, bounds can be callables f(*earlier_args) → (lo, hi).

The variables are ordered as spatial.time_integration_vars.

Returns:

List of bounds, one per integration variable. Each is either (lo, hi) or a callable for dependent bounds.

Parameters:
Return type:

list

to_latex()[source]

Render the decomposed integrand with explicit coordinates.

Shows direction identifications, time-ordered integrals, factored R product, and the component sum over C propagators.

Return type:

str

integrate_moment_qmc(lambda_f, cache, t_min=0.0, direction=0, n_samples=16384, seed=None, positions=None, integrate_over=None)[source]

Integrate over all time variables using Quasi-Monte Carlo.

Scalar-loop sibling of integrate_moment_qmc_vectorized(); see that method’s docstring for the integrate_over kwarg semantics.

Parameters:
Return type:

tuple[float, float]

integrate_moment_qmc_vectorized(lambda_f, cache, t_min=0.0, direction=0, n_samples=16384, seed=None, positions=None, integrate_over=None)[source]

Vectorized QMC integration — no Python loop over samples.

Same algorithm as integrate_moment_qmc() but evaluates all Sobol samples simultaneously using batch propagator lookups. Requires precompute_C_table for the batch C evaluation.

For the iso_R + diag_C case with scalar coupling (no propagator indices), the integrand at each sample is:

coupling × prod_R R(t_l, t_r) × prod_C C_diag(t_l, t_r)

where all R and C values are looked up in batch.

Parameters:

  • positions (dict[str, float] | None) – Optional mapping {point_name: spatial_coord} (e.g. {'x': 0.0, 'y': 0.5}). When the cache has a spatial table built (any homogeneity kind), the x-coordinate of each propagator endpoint is derived from its direction group’s external-point position and flows through cache.C_at_batch. Ignored when no spatial table is built (legacy behaviour — C_diagonal_batch has no x-dependence regardless).

  • integrate_over (Any) –

    Controls which external points have their time integrated over [t_min, lambda_f].

    • None (default, physics observable): all externals held fixed at ``lambda_f`` — this gives the equal-time correlator ⟨φ(t_f) φ(t_f)⟩ that demo notebooks and MC data compare against.

    • "all": all externals integrated — the time-integrated moment ⟨∫₀^{t_f}φ(t) dt · ∫₀^{t_f}φ(t') dt'⟩. Useful e.g. for weak-lensing line-of-sight integrals of the Sachs field.

    • Iterable of external point names (subset): those listed are integrated, the rest fixed at lambda_f. Enables mixed observables such as an integrated source × a detector field.

    Internal vertex times are always integrated; their causal parent-map uses the (possibly fixed) external times as upper bounds.

  • lambda_f (float)

  • cache (PropagatorCache)

  • t_min (float)

  • direction (Any)

  • n_samples (int)

  • seed (int | None)

Returns:

(estimate, std_error)

Return type:

tuple[float, float]

integrate_moment_gauss_legendre(lambda_f, cache, t_min=0.0, direction=0, n_gauss=8, positions=None, integrate_over=None)[source]

Tensor-product Gauss-Legendre quadrature on the causal simplex.

Mirrors integrate_moment_qmc_vectorized() node-for-node – the SAME causal-simplex mapping (parents → upper bounds → Jacobians) and the SAME vectorised batch path through PropagatorCache.R_time_batch() and PropagatorCache.C_at_batch() – but replaces the Sobol n_samples quasi-random nodes with the deterministic n_gauss^d tensor product of 1-D Gauss-Legendre nodes (mapped from [-1, 1] to [0, 1]).

For smooth integrands (a finite product of exponentials – the typical R/C/κ kernel structure) this gives exponential convergence in n_gauss, vastly outperforming Sobol QMC at modest dimensionality. In particular, demo2’s FK channel (4D smooth integrand on a causal simplex of area t_f^2 / 2) is dominated at large t_f by a narrow peak band of area ~ σ_t/γ near the upper-right corner; Sobol QMC under-resolves it unless n_samples is enormous, while n_gauss=8 (4096 nodes for d=4) recovers the notebook’s hand-derived value to ~5 sig figs.

Cost trade-off. Tensor-product GL scales as n_gauss^d – fine for d ≤ 5 (demo2 FK), fast at d ≤ 4, painful at d ≥ 7. For high-d integrands prefer method='qmc_vectorized' with a large n_samples instead.

Parameters:
Returns:

(estimate, 0.0) – GL is deterministic so there is no statistical error to report. The 0.0 mirrors the return shape of the QMC variants for downstream code.

Return type:

tuple[float, float]

integrate_moment_nquad(lambda_f, cache, t_min=0.0, direction=0, positions=None, integrate_over=None)[source]

Integrate over time variables using nested adaptive quadrature.

See integrate_moment_qmc_vectorized() for the integrate_over kwarg semantics (external-time partition into integrated vs fixed-at-lambda_f).

Parameters:
Return type:

tuple[float, float]

sft_wick.evaluate.integrate_moment(integrand, lambda_f, cache, t_min=0.0, direction=0, method='qmc', n_samples=16384, seed=None, positions=None, integrate_over=None, n_gauss=8)[source]

Integrate a diagram’s contribution over all time variables.

Convenience wrapper around DiagramIntegrand.integrate_moment_qmc(), DiagramIntegrand.integrate_moment_qmc_vectorized(), and DiagramIntegrand.integrate_moment_nquad().

Dispatch logic. With method='qmc' (default) the function auto-selects the fastest QMC path compatible with the supplied cache:

  • If the cache supports batched C evaluation (either PropagatorCache.precompute_C_table has been called, or a custom cache implements C_diagonal_batch natively) and cache.model.iso_R is true → integrate_moment_qmc_vectorized() (~200× faster than the scalar loop on typical workloads).

  • Otherwise → integrate_moment_qmc() (scalar Python loop).

Users who explicitly want the scalar loop (e.g. for debugging or with a non-batch-capable custom cache) can pass method='qmc_scalar'.

For best performance, call PropagatorCache.precompute_C_table() before integrating so the vectorised path is selected automatically.

Parameters:

  • integrand (DiagramIntegrand) – A DiagramIntegrand built from a DiagramTerm.

  • lambda_f (float) – Upper integration limit for external times.

  • cache (PropagatorCache) – A PropagatorCache (or compatible custom cache).

  • t_min (float) – Lower time bound (default 0).

  • direction (Any) – Direction value for all spatial points.

  • method (str) – 'qmc' (default, auto-selects vectorised vs scalar), 'qmc_scalar' (force scalar loop), 'qmc_vectorized' (force vectorised; raises if cache doesn’t support batch), or 'nquad' (nested adaptive quadrature).

  • n_samples (int) – Number of Sobol samples (QMC only, should be a power of 2).

  • seed (int | None) – Random seed for reproducibility (QMC only).

  • positions (dict[str, float] | None)

  • integrate_over (Any)

  • n_gauss (int)

Returns:

(estimate, error) tuple.

Return type:

tuple[float, float]

sft_wick.evaluate.integrate_diagrams(diagram_terms, coupling_values, lambda_f, cache, t_min=0.0, direction=0, method='qmc', n_samples=16384, seed=None, fixed_indices=None, n_jobs=1, positions=None, integrate_over=None, n_gauss=8)[source]

Integrate a batch of diagram terms, optionally in parallel.

Builds integrands and evaluates all diagrams, returning the total and per-diagram results. When n_jobs != 1, diagrams are evaluated in parallel using joblib.

Note

n_jobs defaults to 1 (sequential), matching the L1 Expansion.evaluate / Expansion.sweep defaults. Pass -1 to use all CPU cores. Sequential and parallel paths are bit-identical when the QMC seed is fixed.

Parameters:

  • diagram_terms (list) – List of DiagramTerm.

  • coupling_values (dict) – Coupling tensor arrays (e.g. {'F': F_code}).

  • lambda_f (float) – Upper integration limit.

  • cache (PropagatorCache) – A PropagatorCache.

  • t_min (float) – Lower time bound (default 0).

  • direction (Any) – Direction value for all spatial points.

  • method (str) – 'qmc' or 'nquad'.

  • n_samples (int) – Number of Sobol samples (QMC only).

  • seed (int | None) – Random seed (QMC only).

  • fixed_indices (dict[str, int] | None) – Optional pinned index values for build_integrand().

  • n_jobs (int) – Number of parallel jobs. 1 = sequential (default, safe under nested-callers), -1 = all CPU cores. A built-in guard falls back to sequential for len(diagram_terms) <= 2 to avoid joblib’s ~1 s startup overhead on trivial batches. Requires joblib when n_jobs != 1.

  • positions (dict[str, Any] | None)

  • integrate_over (Any)

  • n_gauss (int)

Returns:

(total, details) where total is the scalar sum and details is a list of (estimate, error) per diagram.

Return type:

tuple[float, list[tuple[float, float]]]

sft_wick.evaluate.integrate_two_point_qmc(integrands, t_f, positions, cache, t_min=0.0, n_samples=16384, seed=None)[source]

Vectorised QMC integration for two-point correlation functions.

Evaluates the sum of diagram integrands at external time t_f with spatial points at the positions given by positions. This is the standard entry point for computing correlators such as ⟨φ(x, t) φ(y, t)⟩ where x and y may differ.

Unlike integrate_moment(), which assigns a single direction value to all spatial points, this function supports per-point spatial positions and applies the appropriate spatial factor to each correlation propagator automatically.

The spatial factor for a C-propagator connecting points at positions n₁ and n₂ is computed as:

κ²(n₁, t_ref; n₂, t_ref) / κ²(0, t_ref; 0, t_ref)

which is exact for separable kernels (where the factor is time-independent) and a good approximation otherwise.

Requires PropagatorCache.precompute_C_table() to have been called for the batch C evaluation fast path.

Parameters:

  • integrands (list[DiagramIntegrand]) – List of DiagramIntegrand objects (one per non-vanishing Feynman diagram).

  • t_f (float) – External (observation) time — all external points are set to this time.

  • positions (dict[str, float]) – {point_name: position} mapping each spatial point name (e.g. "x", "y") to a scalar spatial coordinate. Points not in the dict default to 0.

  • cache (PropagatorCache) – A PropagatorCache with pre-computed C table.

  • t_min (float) – Lower time bound (default 0).

  • n_samples (int) – Number of Sobol samples (should be a power of 2).

  • seed (int | None) – Random seed for the Sobol sequence.

Returns:

(estimate, std_error) for the summed two-point function.

Return type:

tuple[float, float]