mixle.utils.special module

Defines the log-pseudo-determinant, polgamma, trigamma, and digamma inverse functions.

This module is the canonical home for mixle’s hand-rolled stable-math helpers (log_erfcx, trigamma, digammainv, log1mexp, logsubexp, logsumexp, softmax, softmax_rows, valid_integer). Import from here rather than re-implementing privately so the numerically-careful versions stay in one place.

log_erfcx(x, out=None)

Stable natural log of the scaled complementary error function log(erfcx(x)).

erfcx(x) = exp(x**2) * erfc(x) is the scaled complementary error function, the workhorse of the exponentially-modified-Gaussian tail. A naive log(erfcx(x)) blows up at both ends: for large positive x it underflows (erfcx -> 0 so log -> -inf) and for large negative x it overflows (the exp(x**2) factor in erfcx -> inf so log -> inf). Three branches keep it finite and accurate everywhere:

• large positive x: the asymptotic series erfcx(x) ~ 1/(x*sqrt(pi)) * (1 - 1/(2 x^2) + 3/(4 x^4) - ...), in log space;

• large negative x: log(erfcx(x)) = x**2 + log(erfc(x)) with erfc(x) -> 2 finite, so no overflow;

• moderate x: the direct log(erfcx(x)).

The branches match log(erfcx) to machine precision on their overlaps. This is what keeps the EMG tail from underflowing/overflowing.

Parameters:

• x (ndarray | float) – Array-like or scalar argument.

• out (ndarray | None) – Optional output array (only used for the array path).

Returns:

log(erfcx(x)) as a numpy array (array input) or float (scalar input).

Return type:

ndarray | float

stirling2(n, k)

Stirling number of the second kind S(n, k).

Counts the number of ways to partition n labeled objects into k non-empty unlabeled subsets. Computed with the standard recurrence S(n, k) = k*S(n-1, k) + S(n-1, k-1) using exact integer arithmetic.

Parameters:

• n (int) – Number of objects (n >= 0).

• k (int) – Number of subsets (k >= 0).

Returns:

Integer value of S(n, k); 0 when k > n or when exactly one of n, k is 0.

Return type:

int

logpdet(x_mat)

Computes the log-pseudo-determinant for a symmetric dense matrix.

Parameters:

x_mat (np.ndarray) – 2-d Numpy array representing a matrix.

Returns:

float, log-pseudo-determinant.

Return type:

float

trigamma(y, out=None)

Trigamma function.

Parameters:

• y (Array-like) – An array-like or float/int.

• out (np.ndarray)

Returns:

Numpy array of trigamma function evaluated at y.

Return type:

ndarray | float

digammainv(y)

Inverse digamma function evaluated on y.

Parameters:

y (Union[np.ndarray, float]) – Numpy array of values to be evaluated or single value.

Returns:

Numpy array if y is numpy array else float.

Return type:

ndarray | float

log1mexp(x)

Return log(1 - exp(x)) for x <= 0, stable across the whole range.

Uses the two-regime split (Mächler, “Accurately Computing log(1 - exp(-|a|))”): log(-expm1(x)) when exp(x) is small and log1p(-exp(x)) when it is close to 1, so 1 - exp(x) is never formed by a catastrophically cancelling subtraction. Returns -inf for x >= 0 (where 1 - exp(x) <= 0).

Parameters:

x (float) – A non-positive log-probability log p with p in [0, 1].

Returns:

log(1 - exp(x)); -inf when x >= 0.

Return type:

float

logsubexp(log_hi, log_lo)

Return log(exp(log_hi) - exp(log_lo)) for log_hi >= log_lo, computed stably.

Evaluates log_hi + log1mexp(log_lo - log_hi) so a far-tail difference whose two operands are individually indistinguishable from 0 (or 1) in probability space still returns a finite large-negative log-mass instead of log(0) = -inf. Returns -inf when the difference is non-positive (log_hi <= log_lo).

Parameters:

• log_hi (float) – Log of the larger operand.

• log_lo (float) – Log of the smaller operand.

Returns:

log(exp(log_hi) - exp(log_lo)); -inf if log_hi <= log_lo.

Return type:

float

logsumexp(a, axis=None)

Stable log(sum(exp(a))) via the max-shift trick.

Thin wrapper over scipy.special.logsumexp() providing mixle’s canonical scalar/array fallback. The axis=None (full reduction) result is returned as a Python float to match the private re-implementations this replaces; a reduced array is returned otherwise. Empty input reduces to -inf and a non-finite running max propagates (+inf stays +inf).

Parameters:

• a (Any) – Array-like of log-values.

• axis (int | None) – Axis to reduce over, or None for a full reduction.

Returns:

float when axis is None; a numpy array otherwise.

Return type:

Any

softmax(log_scores, axis=-1)

Numerically stable softmax of log_scores along axis, with an all--inf guard.

Subtracts the per-slice maximum before exponentiating. A slice that is entirely -inf (no finite log-score) has no defined softmax and would otherwise yield nan; it is filled with a uniform distribution 1 / n over that axis instead.

Parameters:

• log_scores (ndarray) – Array of log-scores.

• axis (int) – Axis along which to normalize (default the last axis).

Returns:

An array the same shape as log_scores whose axis slices each sum to 1.

Return type:

ndarray

softmax_rows(log_scores)

Row-wise (axis=1) softmax of a (B, K) log-score matrix, with an all--inf guard.

Convenience wrapper for softmax() with axis=1: a row that is entirely -inf is replaced by the uniform distribution 1 / K instead of yielding nan.

Parameters:

log_scores (ndarray) – A (B, K) matrix of log-scores.

Returns:

A (B, K) matrix whose rows each sum to 1.

Return type:

ndarray

valid_integer(x, *, nonneg=False)

Return whether x is a finite (optionally non-negative) integer value.

Coerces x to float and checks it is finite and integer-valued (floor(x) == x). Any coercion failure returns False. With nonneg=True the value must additionally be >= 0.

Parameters:

• x (Any) – Candidate value.

• nonneg (bool) – If True, require x >= 0 (e.g. a count); if False, allow negatives.

Returns:

True if x is a valid (non-negative, if requested) integer; False otherwise.

Return type:

bool