core.py — Computational Backend

The core module contains all numerical computation for ToFUL. The frontend (app.py) calls the public functions at the bottom of this module; it imports nothing from SciPy or NumPy directly.

Architecture 

toful_parser.py
    └── normalise_for_eval()    ← clean user input
    └── build_safe_dict()       ← eval namespace

core.py
    ├── _build_safe_dict()      ← canonical eval namespace
    ├── _cached_compile()       ← AST compiled once per expression
    ├── _eval_scalar()          ← scalar evaluation
    ├── _eval_array()           ← vectorised NumPy evaluation
    │
    ├── Convergence cascade (discrete)
    │   ├── _wynn_epsilon()
    │   ├── _aitken_delta2()
    │   ├── _cohen_villegas_zagier()
    │   └── _sum_series_with_acceleration()
    │
    ├── Quadrature (continuous)
    │   ├── _gauss_laguerre_integrate()
    │   ├── _gauss_hermite_integrate()
    │   ├── _mpmath_integrate()
    │   └── _quadrature_integrate()   ← dispatch
    │
    ├── Symbolic path
    │   └── _try_sympy_moment()
    │
    ├── Validation
    │   ├── validate_drv_probabilities()   ← public
    │   └── validate_crv_pdf()             ← public
    │
    ├── Moment computation
    │   ├── _compute_drv_moment()
    │   ├── _compute_crv_moment()
    │   └── compute_moments_parallel()     ← public
    │
    ├── Statistical summary
    │   └── compute_statistical_summary()  ← public
    │
    └── Range parsing
        ├── parse_range_input()            ← public
        └── parse_continuous_bound()       ← public

Data Classes 

class core.SeriesPattern(kind, params, description='')[source]

Detected pattern of a discrete sequence.

Parameters:

kind (str)
params (Dict)
description (str)

kind: str

params: Dict

description: str = ''

class core.SeriesAnalysis(value, converged, terms_used, method, info, partial_sums=<factory>, term_values=<factory>)[source]

Diagnostic output from a single discrete series computation.

Parameters:

value (float)
converged (bool)
terms_used (int)
method (str)
info (str)
partial_sums (List[float])
term_values (List[float])

value: float

converged: bool

terms_used: int

method: str

info: str

partial_sums: List[float]

term_values: List[float]

class core.ValidationResult(is_valid, message, integral_or_sum, analysis=<factory>)[source]

Returned by the probability validators.

Parameters:

is_valid (bool)
message (str)
integral_or_sum (float)
analysis (Dict)

is_valid: bool

message: str

integral_or_sum: float

analysis: Dict

class core.MomentResult(order, reference, value, method, converged, terms_or_nodes, info)[source]

Result for a single moment computation.

Parameters:

order (int)
reference (float)
value (float)
method (str)
converged (bool)
terms_or_nodes (int)
info (str)

order: int

reference: float

value: float

method: str

converged: bool

terms_or_nodes: int

info: str

class core.StatisticalSummary(mean, variance=None, std_dev=None, skewness=None, kurtosis=None, excess_kurtosis=None)[source]

Derived statistical measures computed from central moments.

Parameters:

mean (float)
variance (float | None)
std_dev (float | None)
skewness (float | None)
kurtosis (float | None)
excess_kurtosis (float | None)

mean: float

variance: float | None = None

std_dev: float | None = None

skewness: float | None = None

kurtosis: float | None = None

excess_kurtosis: float | None = None

Public API 

Validation 

core.validate_drv_probabilities(expr_str, x_values, is_infinite=False, tol=0.01, extra=None)[source]

Validate that a discrete PMF sums to 1 and has no negative values.

Parameters:

expr_str (str) – Eval-normalised PMF expression.
x_values (array-like) – Support values (finite window or extended series for infinite).
is_infinite (bool) – Whether the support is effectively unbounded beyond x_values.
tol (float) – Tolerance on abs(sum - 1.0). Loosened for non-converged infinite series due to truncation error.
extra (dict, optional)

Return type:

ValidationResult

core.validate_crv_pdf(expr_str, bounds, precision=8, extra=None)[source]

Validate that a continuous PDF integrates to 1 and is non-negative.

Quadrature rule selected based on the integration domain: - Finite [a, b] : scipy.integrate.quad (Gauss-Kronrod, adaptive) - Semi-infinite [0, inf): Gauss-Laguerre nodes (exact for exp weight) - Doubly-infinite : Gauss-Hermite nodes (exact for Gaussian weight) - High precision (>12dp): mpmath tanh-sinh quadrature

Parameters:

expr_str (str) – Eval-normalised PDF expression.
bounds ((float, float)) – Integration limits.
precision (int) – Target decimal places; triggers mpmath above _MPMATH_THRESHOLD.
extra (dict, optional)

Return type:

ValidationResult

Moment Computation 

core.compute_moments_parallel(expr_str, support_or_bounds, max_order, a, is_crv, is_infinite=False, precision=8, tol=1e-12, extra=None, max_workers=4)[source]

Compute moments of orders 1 through max_order in parallel.

Moment computations for different orders are independent. All orders are submitted to a ThreadPoolExecutor; results are collected as they complete. For max_order <= 2 the threading overhead is avoided and computation is sequential.

Parameters:

expr_str (str) – Eval-normalised PMF/PDF expression.
support_or_bounds (np.ndarray or (float, float)) – For DRV: array of support values. For CRV: (lower, upper) tuple.
max_order (int) – Compute moments mu_1 through mu_{max_order}.
a (float) – Reference point (0=raw, mu=central, other=custom).
is_crv (bool) – True for continuous RV, False for discrete.
is_infinite (bool) – For DRV: whether support is unbounded beyond support array.
precision (int) – Target decimal places (affects quadrature tolerance and mpmath).
tol (float) – Convergence tolerance for discrete series.
extra (dict, optional) – Extra variable bindings.
max_workers (int) – Thread pool size.

Return type:

dict mapping int -> MomentResult

Statistical Summary 

core.compute_statistical_summary(mean, central_moments)[source]

Derive variance, standard deviation, skewness, and kurtosis from central moments (moments computed about the mean).

Standard formulae:: variance = mu_2 std_dev = sqrt(abs(mu_2)) skewness = mu_3 / sigma^3 kurtosis = mu_4 / sigma^4 excess_kurtosis = kurtosis - 3

Parameters:

mean (float) – First raw moment (E[X]).
central_moments (dict) – Mapping order -> MomentResult for central moments (a = mean).

Return type:

StatisticalSummary

Range Parsing 

core.parse_range_input(range_input)[source]

Parse a discrete range string into a NumPy array of support values.

Accepts:: “1, 2, 3, 4, 5” -> finite support “0, 1, 2, 3, …” -> infinite arithmetic series (extended to _MAX_SERIES_TERMS) “1, 2, 4, 8, …” -> infinite geometric series

Parameters:: range_input (str) – Range string, already cleaned by toful_parser.normalise_range_input().
Return type:: (values_array, is_infinite, pattern_description)
Raises:: ValueError – If the string cannot be parsed as a comma-separated list of numbers.

core.parse_continuous_bound(bound_str)[source]

Parse a single bound string to a float, accepting infinity variants.

Recognised infinity tokens: inf, infinity, +inf, +infinity, and their Unicode variants (inf_sign + oo).

Parameters:: bound_str (str) – E.g. “0”, “1.5”, “inf”, “-inf”, “∞”, “-∞”.
Return type:: float

Utilities 

core.to_subscript(n)[source]

Convert a non-negative integer to its Unicode subscript string.

Examples

>>> to_subscript(12)
'12'  # rendered as subscript digits

Parameters:: n (int)
Return type:: str

core.detect_series_pattern(values)[source]

Detect whether seed values form an arithmetic or geometric sequence.

Parameters:: values (list of float) – Seed values from the user range input.
Return type:: SeriesPattern

core.generate_extended_series(pattern, max_terms=500)[source]

Generate a NumPy array of the first max_terms values matching the pattern.

Parameters:

pattern (SeriesPattern)
max_terms (int)

Return type:

np.ndarray

Convergence Algorithm Reference 

The following internal functions implement the convergence acceleration cascade. They are not part of the public API but are documented here for contributors and advanced users.

_wynn_epsilon(partial_sums)

Wynn ε-algorithm. Builds a recursive epsilon-table from the list of partial sums and returns the best limit estimate from the even-column diagonal.

Returns (estimate, converged, info_string).

_aitken_delta2(s0, s1, s2)

Aitken Δ² three-point extrapolation. Returns the extrapolated limit from three consecutive partial sums.

_cohen_villegas_zagier(expr_str, n_terms)

Cohen-Villegas-Zagier algorithm for alternating series. Computes near-machine-precision estimate in exactly n_terms evaluations.

Returns (estimate, True, info_string).

_sum_series_with_acceleration(expr_str, x_values, r, a, tol)

Main discrete summation engine. Evaluates Σ (x-a)^r · f(x) over x_values and applies the full convergence cascade.

Returns a SeriesAnalysis instance.

Quadrature Reference 

_gauss_laguerre_integrate(func, n=64): Gauss-Laguerre quadrature for \(\int_0^\infty f(x)\,dx\). Uses 64 nodes. Returns (estimate, error_estimate).
_gauss_hermite_integrate(func, n=64): Gauss-Hermite quadrature for \(\int_{-\infty}^{\infty} f(x)\,dx\). Uses 64 nodes. Returns (estimate, error_estimate).
_mpmath_integrate(func, lower, upper, precision): mpmath tanh-sinh (double-exponential) quadrature. Active when Calc precision ≥ 13. Returns (estimate, error_estimate, "mpmath-tanh-sinh").
_quadrature_integrate(func, lower, upper, precision, ...): Dispatch function. Selects the best quadrature method based on domain and precision, cross-checks results, returns (value, error, method_name).

Configuration Constants 

_MPMATH_THRESHOLD = 12    # Calc precision above which mpmath activates
_MAX_SERIES_TERMS = 500   # Hard cap on infinite-series truncation
_CONVERGENCE_TOL  = 1e-12 # Default convergence tolerance
_GAUSS_NODES      = 64    # Gauss-Laguerre / Hermite quadrature nodes

Legacy Compatibility 

The following classes and function from the original codebase are retained as thin shims over the new API. New code should use the public functions above directly.

class core.InfiniteSeriesHandler[source]: Legacy shim — use detect_series_pattern() and generate_extended_series().

class core.EnhancedProbabilityValidator[source]: Legacy shim — use validate_drv_probabilities() / validate_crv_pdf().

class core.EnhancedMomentCalculator[source]: Legacy shim — use compute_moments_parallel() or _compute_drv_moment().

core.rth_moment(pmf, support, r, c=0, tol=1e-12, max_iter=1000000)[source]: Legacy function — kept for backward compatibility.