Mathematical Functions

Mathematical functions and utilities.

This module contains a wide variety of mathematical functions including: - Basic matrix operations - Combinatorics (permutations, combinations) - Geometry primitives - Interpolation methods - Numerical integration - Signal processing - Special functions - Statistics and distributions

Basic Matrix Operations

Basic matrix operations and constructions.

This module provides: - Matrix decompositions (Cholesky, SVD-based, QR) - Special matrix constructions (Vandermonde, Toeplitz, Hankel, etc.) - Matrix vectorization operations (vec, unvec, Kronecker products)

Decompositions

Matrix decomposition utilities.

This module provides matrix decomposition functions that wrap numpy/scipy with consistent APIs matching the MATLAB TrackerComponentLibrary conventions.

pytcl.mathematical_functions.basic_matrix.decompositions.chol_semi_def(A, upper=False, tol=1e-10)[source]

Compute Cholesky decomposition of a positive semi-definite matrix.

For positive semi-definite matrices that may be singular or near-singular, this function uses eigenvalue decomposition with thresholding to produce a valid Cholesky-like factor.

Parameters:

A (array_like) – Symmetric positive semi-definite matrix of shape (n, n).
upper (bool, optional) – If True, return upper triangular factor R such that A ≈ R.T @ R. If False (default), return lower triangular factor L such that A ≈ L @ L.T.
tol (float, optional) – Eigenvalues below tol * max(eigenvalues) are treated as zero. Default is 1e-10.

Returns:

L_or_R – Lower (or upper if upper=True) triangular Cholesky factor. Shape is (n, n).

Return type:

ndarray

Examples

>>> A = np.array([[4, 2], [2, 1]])  # Singular but positive semi-definite
>>> L = chol_semi_def(A)
>>> np.allclose(L @ L.T, A)
True

See also

numpy.linalg.cholesky: Standard Cholesky for positive definite matrices.
scipy.linalg.cholesky: Standard Cholesky with more options.

pytcl.mathematical_functions.basic_matrix.decompositions.tria(A)[source]

Compute lower triangular square root factor of a symmetric matrix.

Given a symmetric positive semi-definite matrix A, returns a lower triangular matrix S such that A = S @ S.T. This is useful for square-root Kalman filtering implementations.

Parameters:: A (array_like) – Symmetric positive semi-definite matrix of shape (n, n).
Returns:: S – Lower triangular matrix of shape (n, n) such that A ≈ S @ S.T.
Return type:: ndarray

Notes

This function is equivalent to the lower Cholesky factor for positive definite matrices. For semi-definite matrices, it uses the eigenvalue-based approach from chol_semi_def.

See also

chol_semi_def: More general function with tolerance control.
tria_sqrt: Square root of concatenated matrices for filter updates.

pytcl.mathematical_functions.basic_matrix.decompositions.tria_sqrt(A, B=None)[source]

Compute triangular square root of [A, B] @ [A, B].T.

This is commonly used in square-root Kalman filter implementations where we need to compute the square root of a sum of outer products.

Parameters:

A (array_like) – Matrix of shape (n, m).
B (array_like, optional) – Matrix of shape (n, p). If None, computes sqrt of A @ A.T.

Returns:

S – Lower triangular matrix of shape (n, n) such that S @ S.T = A @ A.T + B @ B.T (or just A @ A.T if B is None).

Return type:

ndarray

Notes

Uses QR decomposition for numerical stability: [A, B].T = Q @ R implies [A, B] @ [A, B].T = R.T @ R

Examples

>>> A = np.random.randn(3, 4)
>>> B = np.random.randn(3, 2)
>>> S = tria_sqrt(A, B)
>>> expected = A @ A.T + B @ B.T
>>> np.allclose(S @ S.T, expected)
True

pytcl.mathematical_functions.basic_matrix.decompositions.pinv_truncated(A, tol=None, rank=None)[source]

Compute truncated pseudo-inverse using SVD.

Computes the Moore-Penrose pseudo-inverse with explicit control over which singular values are retained.

Parameters:

A (array_like) – Input matrix of shape (m, n).
tol (float, optional) – Singular values below tol * max(singular values) are set to zero. Default is max(m, n) * eps * max(singular values).
rank (int, optional) – If provided, only the largest rank singular values are used. Overrides tol if specified.

Returns:

A_pinv – Pseudo-inverse of A with shape (n, m).

Return type:

ndarray

Examples

>>> A = np.array([[1, 2], [3, 4], [5, 6]])
>>> A_pinv = pinv_truncated(A)
>>> np.allclose(A @ A_pinv @ A, A)
True

See also

numpy.linalg.pinv: Standard pseudo-inverse.
scipy.linalg.pinv: Scipy version with rcond parameter.

pytcl.mathematical_functions.basic_matrix.decompositions.matrix_sqrt(A, method='schur')[source]

Compute the principal matrix square root.

Finds matrix S such that S @ S = A. This is different from the Cholesky factor which satisfies L @ L.T = A.

Parameters:

A (array_like) – Square matrix of shape (n, n).
method ({'schur', 'eigenvalue', 'denman_beavers'}, optional) – Algorithm to use: - ‘schur’: Uses Schur decomposition (default, most stable). - ‘eigenvalue’: Uses eigenvalue decomposition (faster for normal matrices). - ‘denman_beavers’: Iterative method (good for ill-conditioned cases).

Returns:

S – Principal square root matrix of shape (n, n).

Return type:

ndarray

Notes

The principal square root has eigenvalues with positive real parts.

Examples

>>> A = np.array([[4, 0], [0, 9]])
>>> S = matrix_sqrt(A)
>>> np.allclose(S @ S, A)
True
>>> S
array([[2., 0.],
       [0., 3.]])

pytcl.mathematical_functions.basic_matrix.decompositions.rank_revealing_qr(A, tol=None)[source]

Compute rank-revealing QR decomposition with column pivoting.

Computes A[:, P] = Q @ R where P is a permutation that reveals the numerical rank of A through the diagonal of R.

Parameters:

A (array_like) – Input matrix of shape (m, n).
tol (float, optional) – Tolerance for determining numerical rank. Diagonal elements of R below tol * |R[0,0]| indicate rank deficiency. Default is max(m, n) * eps * |R[0,0]|.

Returns:

Q (ndarray) – Orthogonal matrix of shape (m, k) where k = min(m, n).
R (ndarray) – Upper triangular matrix of shape (k, n).
P (ndarray) – Permutation indices such that A[:, P] = Q @ R.
rank (int) – Numerical rank determined by tolerance.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[int64]], int]

Examples

>>> A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])  # Rank 2
>>> Q, R, P, rank = rank_revealing_qr(A)
>>> rank
2

pytcl.mathematical_functions.basic_matrix.decompositions.null_space(A, tol=None)[source]

Compute orthonormal basis for the null space of A.

Parameters:

A (array_like) – Input matrix of shape (m, n).
tol (float, optional) – Singular values below tol are considered zero. Default is max(m, n) * eps * max(singular values).

Returns:

N – Orthonormal basis for null(A) with shape (n, k) where k is the dimension of the null space.

Return type:

ndarray

Examples

>>> A = np.array([[1, 2, 3], [4, 5, 6]])
>>> N = null_space(A)
>>> np.allclose(A @ N, 0, atol=1e-10)
True

pytcl.mathematical_functions.basic_matrix.decompositions.range_space(A, tol=None)[source]

Compute orthonormal basis for the range (column space) of A.

Parameters:

A (array_like) – Input matrix of shape (m, n).
tol (float, optional) – Singular values below tol are considered zero. Default is max(m, n) * eps * max(singular values).

Returns:

R – Orthonormal basis for range(A) with shape (m, r) where r is the rank.

Return type:

ndarray

Examples

>>> A = np.array([[1, 2], [3, 6], [5, 10]])  # Rank 1
>>> R = range_space(A)
>>> R.shape
(3, 1)

Special Matrices

Special matrix constructions.

This module provides functions for constructing special matrices commonly used in numerical algorithms and signal processing.

pytcl.mathematical_functions.basic_matrix.special_matrices.vandermonde(x, n=None, increasing=False)[source]

Construct a Vandermonde matrix.

The Vandermonde matrix has columns that are powers of the input vector. By default (decreasing order): V[i,j] = x[i]^(n-1-j) With increasing=True: V[i,j] = x[i]^j

Parameters:

x (array_like) – Input vector of length m.
n (int, optional) – Number of columns. Default is len(x).
increasing (bool, optional) – If True, powers increase left to right. Default is False.

Returns:

V – Vandermonde matrix of shape (m, n).

Return type:

ndarray

Examples

>>> vandermonde([1, 2, 3], 3)
array([[1., 1., 1.],
       [4., 2., 1.],
       [9., 3., 1.]])

>>> vandermonde([1, 2, 3], 3, increasing=True)
array([[1., 1., 1.],
       [1., 2., 4.],
       [1., 3., 9.]])

pytcl.mathematical_functions.basic_matrix.special_matrices.toeplitz(c, r=None)[source]

Construct a Toeplitz matrix.

A Toeplitz matrix has constant diagonals. It is fully specified by its first column and first row.

Parameters:

c (array_like) – First column of the matrix.
r (array_like, optional) – First row of the matrix. If None, r = conjugate(c) is assumed (Hermitian Toeplitz). Note: r[0] is ignored; c[0] is used.

Returns:

T – Toeplitz matrix.

Return type:

ndarray

Examples

>>> toeplitz([1, 2, 3], [1, 4, 5])
array([[1., 4., 5.],
       [2., 1., 4.],
       [3., 2., 1.]])

See also

scipy.linalg.toeplitz: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.hankel(c, r=None)[source]

Construct a Hankel matrix.

A Hankel matrix has constant anti-diagonals. It is fully specified by its first column and last row.

Parameters:

c (array_like) – First column of the matrix.
r (array_like, optional) – Last row of the matrix. If None, zeros are used except for c[-1]. Note: r[0] should equal c[-1].

Returns:

H – Hankel matrix.

Return type:

ndarray

Examples

>>> hankel([1, 2, 3], [3, 4, 5])
array([[1., 2., 3.],
       [2., 3., 4.],
       [3., 4., 5.]])

See also

scipy.linalg.hankel: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.circulant(c)[source]

Construct a circulant matrix.

A circulant matrix is a special Toeplitz matrix where each row is a cyclic shift of the row above it.

Parameters:: c (array_like) – First column of the matrix.
Returns:: C – Circulant matrix of shape (n, n) where n = len(c).
Return type:: ndarray

Examples

>>> circulant([1, 2, 3])
array([[1., 3., 2.],
       [2., 1., 3.],
       [3., 2., 1.]])

Notes

Circulant matrices are diagonalized by the DFT matrix.

See also

scipy.linalg.circulant: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.block_diag(*arrs)[source]

Create a block diagonal matrix from provided arrays.

Parameters:: *arrs (sequence of array_like) – Input arrays. Each array becomes a block on the diagonal.
Returns:: D – Block diagonal matrix.
Return type:: ndarray

Examples

>>> A = np.array([[1, 2], [3, 4]])
>>> B = np.array([[5, 6, 7]])
>>> block_diag(A, B)
array([[1., 2., 0., 0., 0.],
       [3., 4., 0., 0., 0.],
       [0., 0., 5., 6., 7.]])

See also

scipy.linalg.block_diag: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.companion(c)[source]

Create a companion matrix.

The companion matrix is used for polynomial root finding. For a monic polynomial p(x) = x^n + c_{n-1}*x^{n-1} + … + c_1*x + c_0, the eigenvalues of the companion matrix are the roots of p(x).

Parameters:: c (array_like) – Coefficients of the polynomial (excluding leading 1), in order [c_{n-1}, c_{n-2}, …, c_1, c_0] or [c_0, c_1, …, c_{n-1}] depending on convention used.
Returns:: C – Companion matrix of shape (n, n) where n = len(c).
Return type:: ndarray

Examples

>>> # Polynomial: x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)
>>> c = [6, -11, 6]  # Coefficients (negated, reversed)
>>> C = companion(c)

See also

scipy.linalg.companion: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.hilbert(n)[source]

Create a Hilbert matrix.

The Hilbert matrix H has entries H[i,j] = 1 / (i + j + 1). This matrix is notoriously ill-conditioned.

Parameters:: n (int) – Size of the matrix.
Returns:: H – Hilbert matrix of shape (n, n).
Return type:: ndarray

Examples

>>> hilbert(3)
array([[1.        , 0.5       , 0.33333333],
       [0.5       , 0.33333333, 0.25      ],
       [0.33333333, 0.25      , 0.2       ]])

See also

scipy.linalg.hilbert: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.invhilbert(n)[source]

Compute the inverse of the Hilbert matrix.

Uses an exact formula to compute the inverse, which is known to have integer entries.

Parameters:: n (int) – Size of the matrix.
Returns:: H_inv – Inverse of the n x n Hilbert matrix.
Return type:: ndarray

See also

scipy.linalg.invhilbert: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.hadamard(n)[source]

Construct a Hadamard matrix.

A Hadamard matrix H satisfies H @ H.T = n * I, where all entries are +1 or -1.

Parameters:: n (int) – Size of the matrix. Must be a power of 2, or 1, 2.
Returns:: H – Hadamard matrix of shape (n, n).
Return type:: ndarray
Raises:: ValueError – If n is not a power of 2.

Examples

>>> hadamard(4)
array([[ 1.,  1.,  1.,  1.],
       [ 1., -1.,  1., -1.],
       [ 1.,  1., -1., -1.],
       [ 1., -1., -1.,  1.]])

See also

scipy.linalg.hadamard: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.dft_matrix(n, normalized=False)[source]

Construct the DFT (Discrete Fourier Transform) matrix.

The DFT matrix F has entries F[j,k] = exp(-2*pi*i*j*k/n).

Parameters:

n (int) – Size of the matrix.
normalized (bool, optional) – If True, return unitary DFT matrix (scaled by 1/sqrt(n)). Default is False.

Returns:

F – DFT matrix of shape (n, n), complex-valued.

Return type:

ndarray

Examples

>>> F = dft_matrix(4)
>>> x = np.array([1, 2, 3, 4])
>>> np.allclose(F @ x, np.fft.fft(x))
True

See also

scipy.linalg.dft: Equivalent scipy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.kron(a, b)[source]

Compute the Kronecker product of two arrays.

Parameters:

a (array_like) – First input array.
b (array_like) – Second input array.

Returns:

K – Kronecker product of a and b.

Return type:

ndarray

Examples

>>> a = np.array([[1, 2], [3, 4]])
>>> b = np.array([[1, 0], [0, 1]])
>>> kron(a, b)
array([[1., 0., 2., 0.],
       [0., 1., 0., 2.],
       [3., 0., 4., 0.],
       [0., 3., 0., 4.]])

See also

numpy.kron: Equivalent numpy function.

pytcl.mathematical_functions.basic_matrix.special_matrices.vec(A)[source]

Vectorize a matrix by stacking its columns.

This is the standard vec operation from matrix calculus.

Parameters:: A (array_like) – Input matrix of shape (m, n).
Returns:: v – Column vector of shape (m*n,) containing columns of A stacked.
Return type:: ndarray

Examples

>>> A = np.array([[1, 2], [3, 4]])
>>> vec(A)
array([1., 3., 2., 4.])

See also

unvec: Inverse operation.

pytcl.mathematical_functions.basic_matrix.special_matrices.unvec(v, m, n)[source]

Reshape a vector back to a matrix (inverse of vec).

Parameters:

v (array_like) – Input vector of length m*n.
m (int) – Number of rows in output matrix.
n (int) – Number of columns in output matrix.

Returns:

A – Matrix of shape (m, n).

Return type:

ndarray

Examples

>>> v = np.array([1, 3, 2, 4])
>>> unvec(v, 2, 2)
array([[1., 2.],
       [3., 4.]])

See also

vec: Forward operation.

pytcl.mathematical_functions.basic_matrix.special_matrices.commutation_matrix(m, n)[source]

Construct the commutation matrix K_{m,n}.

The commutation matrix satisfies K @ vec(A) = vec(A.T) for any m x n matrix A.

Parameters:

m (int) – Number of rows of the matrix to be transposed.
n (int) – Number of columns of the matrix to be transposed.

Returns:

K – Commutation matrix of shape (m*n, m*n).

Return type:

ndarray

Examples

>>> K = commutation_matrix(2, 3)
>>> A = np.array([[1, 2, 3], [4, 5, 6]])
>>> np.allclose(K @ vec(A), vec(A.T))
True

pytcl.mathematical_functions.basic_matrix.special_matrices.duplication_matrix(n)[source]

Construct the duplication matrix D_n.

For a symmetric n x n matrix A, D_n @ vech(A) = vec(A), where vech is the half-vectorization operator.

Parameters:: n (int) – Size of the symmetric matrix.
Returns:: D – Duplication matrix of shape (n*n, n*(n+1)/2).
Return type:: ndarray

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.basic_matrix import duplication_matrix, vec
>>> # Create duplication matrix for 2x2 symmetric matrices
>>> D = duplication_matrix(2)
>>> D.shape
(4, 3)
>>> # For symmetric matrix A = [[1, 2], [2, 3]], half-vec has 3 elements
>>> A = np.array([[1.0, 2.0], [2.0, 3.0]])
>>> vech_A = A[np.tril_indices(2)]  # Half-vectorization [1, 2, 3]
>>> # Duplication matrix should reconstruct full vectorization
>>> vec_A = D @ vech_A
>>> bool(np.allclose(vec_A, vec(A)))
True

pytcl.mathematical_functions.basic_matrix.special_matrices.elimination_matrix(n)[source]

Construct the elimination matrix L_n.

For any n x n matrix A, L_n @ vec(A) = vech(A), where vech is the half-vectorization operator that extracts the lower triangle.

Parameters:: n (int) – Size of the matrix.
Returns:: L – Elimination matrix of shape (n*(n+1)/2, n*n).
Return type:: ndarray

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.basic_matrix import elimination_matrix, vec
>>> # Create elimination matrix for 2x2 matrices
>>> L = elimination_matrix(2)
>>> L.shape
(3, 4)
>>> # For matrix A, extracts unique elements: [A[0,0], A[1,0], A[1,1]]
>>> A = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> # Elimination extracts lower-triangular elements
>>> vech_A = L @ vec(A)
>>> bool(np.allclose(vech_A, [1.0, 3.0, 4.0]))
True

Special Functions

Special mathematical functions.

This module provides special functions commonly used in mathematical physics, signal processing, and statistical applications: - Bessel functions (cylindrical and spherical) - Gamma and beta functions - Error functions - Elliptic integrals - Marcum Q function (radar detection) - Hypergeometric functions - Lambert W function - Debye functions (thermodynamics)

Bessel Functions

Bessel functions and related special functions.

This module provides Bessel functions commonly used in signal processing, antenna theory, and scattering problems in tracking applications.

pytcl.mathematical_functions.special_functions.bessel.besselj(n, x)[source]

Bessel function of the first kind.

Computes J_n(x), the Bessel function of the first kind of order n.

Parameters:

n (int, float, or array_like) – Order of the Bessel function.
x (array_like) – Argument of the Bessel function.

Returns:

J – Values of J_n(x).

Return type:

ndarray

Examples

>>> float(besselj(0, 0))
1.0
>>> besselj(1, np.array([0, 1, 2]))
array([0.        , 0.44005059, 0.57672481])

See also

scipy.special.jv: Bessel function of first kind of real order.

pytcl.mathematical_functions.special_functions.bessel.bessely(n, x)[source]

Bessel function of the second kind (Neumann function).

Computes Y_n(x), the Bessel function of the second kind of order n.

Parameters:

n (int, float, or array_like) – Order of the Bessel function.
x (array_like) – Argument of the Bessel function. Must be positive.

Returns:

Y – Values of Y_n(x).

Return type:

ndarray

Notes

Y_n(x) is singular at x = 0.

Examples

>>> round(float(bessely(0, 1)), 6)
0.088257

See also

scipy.special.yv: Bessel function of second kind of real order.

pytcl.mathematical_functions.special_functions.bessel.besseli(n, x)[source]

Modified Bessel function of the first kind.

Computes I_n(x), the modified Bessel function of the first kind.

Parameters:

n (int, float, or array_like) – Order of the Bessel function.
x (array_like) – Argument of the Bessel function.

Returns:

I – Values of I_n(x).

Return type:

ndarray

Examples

>>> float(besseli(0, 0))
1.0

See also

scipy.special.iv: Modified Bessel function of first kind.

pytcl.mathematical_functions.special_functions.bessel.besselk(n, x)[source]

Modified Bessel function of the second kind.

Computes K_n(x), the modified Bessel function of the second kind.

Parameters:

n (int, float, or array_like) – Order of the Bessel function.
x (array_like) – Argument of the Bessel function. Must be positive.

Returns:

K – Values of K_n(x).

Return type:

ndarray

Notes

K_n(x) is singular at x = 0.

Examples

>>> round(float(besselk(0, 1)), 6)
0.421024
>>> round(float(besselk(1, 2)), 6)
0.139866

See also

scipy.special.kv: Modified Bessel function of second kind.

pytcl.mathematical_functions.special_functions.bessel.besselh(n, k, x)[source]

Hankel function (Bessel function of the third kind).

Computes H^(k)_n(x), the Hankel function of the first (k=1) or second (k=2) kind.

Parameters:

n (int, float, or array_like) – Order of the Hankel function.
k (int) – Kind of Hankel function. Must be 1 or 2.
x (array_like) – Argument of the Hankel function.

Returns:

H – Complex values of H^(k)_n(x).

Return type:

ndarray

Notes

H^(1)_n(x) = J_n(x) + i*Y_n(x) H^(2)_n(x) = J_n(x) - i*Y_n(x)

Examples

>>> h = besselh(0, 1, 1)  # H^(1)_0(1)
>>> round(float(h.real), 6)
0.765198
>>> round(float(h.imag), 6)
0.088257

See also

scipy.special.hankel1: Hankel function of first kind.
scipy.special.hankel2: Hankel function of second kind.

pytcl.mathematical_functions.special_functions.bessel.spherical_jn(n, x, derivative=False)[source]

Spherical Bessel function of the first kind.

Computes j_n(x), the spherical Bessel function of the first kind.

Parameters:

n (int) – Order of the function (non-negative).
x (array_like) – Argument of the function.
derivative (bool, optional) – If True, return the derivative j_n’(x) instead. Default is False.

Returns:

j – Values of j_n(x) or j_n’(x).

Return type:

ndarray

Notes

j_n(x) = sqrt(pi / (2*x)) * J_{n+1/2}(x)

Examples

>>> round(float(spherical_jn(0, 1)), 6)  # sin(1)/1
0.841471
>>> round(float(spherical_jn(0, 1, derivative=True)), 6)  # Derivative
-0.301169

See also

scipy.special.spherical_jn: Spherical Bessel function of first kind.

pytcl.mathematical_functions.special_functions.bessel.spherical_yn(n, x, derivative=False)[source]

Spherical Bessel function of the second kind.

Computes y_n(x), the spherical Bessel function of the second kind.

Parameters:

n (int) – Order of the function (non-negative).
x (array_like) – Argument of the function. Must be positive.
derivative (bool, optional) – If True, return the derivative y_n’(x) instead. Default is False.

Returns:

y – Values of y_n(x) or y_n’(x).

Return type:

ndarray

Examples

>>> round(float(spherical_yn(0, 1)), 6)  # -cos(1)/1
-0.540302

See also

scipy.special.spherical_yn: Spherical Bessel function of second kind.

pytcl.mathematical_functions.special_functions.bessel.spherical_in(n, x, derivative=False)[source]

Modified spherical Bessel function of the first kind.

Computes i_n(x), the modified spherical Bessel function of the first kind.

Parameters:

n (int) – Order of the function (non-negative).
x (array_like) – Argument of the function.
derivative (bool, optional) – If True, return the derivative i_n’(x) instead. Default is False.

Returns:

i – Values of i_n(x) or i_n’(x).

Return type:

ndarray

Examples

>>> round(float(spherical_in(0, 1)), 6)  # sinh(1)/1
1.175201

See also

scipy.special.spherical_in: Modified spherical Bessel function of first kind.

pytcl.mathematical_functions.special_functions.bessel.spherical_kn(n, x, derivative=False)[source]

Modified spherical Bessel function of the second kind.

Computes k_n(x), the modified spherical Bessel function of the second kind.

Parameters:

n (int) – Order of the function (non-negative).
x (array_like) – Argument of the function. Must be positive.
derivative (bool, optional) – If True, return the derivative k_n’(x) instead. Default is False.

Returns:

k – Values of k_n(x) or k_n’(x).

Return type:

ndarray

Examples

>>> round(float(spherical_kn(0, 1)), 6)  # (pi/2) * exp(-1)
0.577864

See also

scipy.special.spherical_kn: Modified spherical Bessel function of second kind.

pytcl.mathematical_functions.special_functions.bessel.airy(x)[source]

Airy functions and their derivatives.

Computes Ai(x), Ai’(x), Bi(x), Bi’(x).

Parameters:

x (array_like) – Argument of the Airy functions.

Returns:

Ai (ndarray) – Airy function Ai(x).
Aip (ndarray) – Derivative of Airy function Ai’(x).
Bi (ndarray) – Airy function Bi(x).
Bip (ndarray) – Derivative of Airy function Bi’(x).

Return type:

tuple[ndarray[Any, Any], ndarray[Any, Any], ndarray[Any, Any], ndarray[Any, Any]]

Examples

>>> Ai, Aip, Bi, Bip = airy(0)
>>> round(float(Ai), 6)
0.355028
>>> round(float(Bi), 6)
0.614927

See also

scipy.special.airy: Airy functions.

pytcl.mathematical_functions.special_functions.bessel.bessel_ratio(n, x, kind='j')[source]

Ratio of Bessel functions J_{n+1}(x) / J_n(x) or I_{n+1}(x) / I_n(x).

Parameters:

n (int or float) – Order of the Bessel function in the denominator.
x (array_like) – Argument of the Bessel function.
kind (str, optional) – Type of Bessel function: ‘j’ for J_n, ‘i’ for I_n. Default is ‘j’.

Returns:

ratio – Values of J_{n+1}(x) / J_n(x) or I_{n+1}(x) / I_n(x).

Return type:

ndarray

Notes

Evaluated with the modified Lentz method (Thompson & Barnett, J. Comput. Phys. 64:490, 1986; Numerical Recipes 3rd ed., Sec. 6.5) applied to the continued fraction that follows from the three-term recurrence:

J_{n+1}(x)/J_n(x) = 1/(2(n+1)/x - 1/(2(n+2)/x - 1/(2(n+3)/x - …))) I_{n+1}(x)/I_n(x) = 1/(2(n+1)/x + 1/(2(n+2)/x + 1/(2(n+3)/x + …)))

J_{n+k}(x) and I_{n+k}(x) are the minimal solutions of their recurrences as k grows, so by Pincherle’s theorem the fractions converge to the exact ratio without ever forming the numerator and denominator separately – the ratio therefore stays finite where a direct quotient of float64 Bessel values underflows to 0/0 (e.g. n = 170, x = 1). Measured against a 50-digit mpmath reference over the grid n in {0, 1, 5, 20, 80, 170, 400} x x in {0.5, 1, 10, 50, 100}, the worst relative error is 9.6e-15 for kind ‘j’ and 1.5e-15 for kind ‘i’ (tests/validation/test_special_functions_audit.py).

The ‘j’ fraction needs roughly max(n, abs(x)) terms to converge; where it misses the iteration cap of 10000 (abs(x) in the thousands and beyond, small n) the direct quotient jv(n+1, x)/jv(n, x) is used instead – machine-accurate there since neither value underflows at large abs(x). Measured worst relative error 7.9e-14 at x in {9000, 15000, 30000}, n in {0, 5}; the all-positive ‘i’ fraction converges within the cap (1.3e-15 measured at x = 30000).

Near a zero of J_n(x) the ratio is well-defined and the fraction converges to it, but roundoff in evaluating the fraction is amplified like the reciprocal of abs(J_n(x)): measured agreement with mpmath at n = 0 loosens from 5e-15 at x = 2.404 to 1e-8 at x = 2.4048255576, i.e. 1e-10 from the first zero of J_0 at x ~= 2.40482555769577, and reaches O(1) relative error (measured 1.4) at the float64 neighbor of the zero. At x = 0 the limit value 0 is returned for every order (the two-sided limit of the ratio); the ratio is nan for non-finite x.

Examples

>>> round(float(bessel_ratio(0, 1)), 6)  # J_1(1) / J_0(1)
0.575081

pytcl.mathematical_functions.special_functions.bessel.bessel_deriv(n, x, kind='j')[source]

Derivative of Bessel function d/dx[B_n(x)].

Parameters:

n (int or float) – Order of the Bessel function.
x (array_like) – Argument of the Bessel function.
kind (str, optional) – Type of Bessel function: ‘j’, ‘y’, ‘i’, or ‘k’. Default is ‘j’.

Returns:

deriv – Values of dB_n(x)/dx.

Return type:

ndarray

Notes

Uses the identity: dJ_n/dx = (J_{n-1}(x) - J_{n+1}(x)) / 2 dY_n/dx = (Y_{n-1}(x) - Y_{n+1}(x)) / 2 dI_n/dx = (I_{n-1}(x) + I_{n+1}(x)) / 2 dK_n/dx = -(K_{n-1}(x) + K_{n+1}(x)) / 2

Examples

>>> round(float(bessel_deriv(0, 1, kind='j')), 6)  # -J_1(1)
-0.440051

pytcl.mathematical_functions.special_functions.bessel.struve_h(n, x)[source]

Struve function H_n(x).

The Struve function is defined by the integral:

H_n(x) = (2/sqrt(pi)) * (x/2)^n * integral from 0 to pi/2 of
         sin(x*cos(t)) * sin^(2n)(t) dt

Parameters:

n (int or float) – Order of the Struve function.
x (array_like) – Argument of the function.

Returns:

H – Values of H_n(x).

Return type:

ndarray

Notes

Related to Bessel functions through: H_0(x) is the particular solution of y’’ + y’/x + y = 2/(pi*x)

Examples

>>> round(float(struve_h(0, 1)), 6)
0.568657

pytcl.mathematical_functions.special_functions.bessel.struve_l(n, x)[source]

Modified Struve function L_n(x).

The modified Struve function is related to the Struve function by: L_n(x) = -i * exp(-i*n*pi/2) * H_n(i*x)

Parameters:

n (int or float) – Order of the modified Struve function.
x (array_like) – Argument of the function.

Returns:

L – Values of L_n(x).

Return type:

ndarray

Examples

>>> round(float(struve_l(0, 1)), 6)
0.710243

pytcl.mathematical_functions.special_functions.bessel.bessel_zeros(n, nt, kind='j')[source]

Zeros of Bessel functions.

Computes the first nt zeros of J_n(x), Y_n(x), or their derivatives.

Parameters:

n (int) – Order of the Bessel function.
nt (int) – Number of zeros to compute.
kind (str, optional) – Type: ‘j’ for J_n zeros, ‘y’ for Y_n zeros, ‘jp’ for J_n’ zeros, ‘yp’ for Y_n’ zeros. Default is ‘j’.

Returns:

zeros – Array of zeros.

Return type:

ndarray

Examples

>>> bessel_zeros(0, 3, kind='j')  # First 3 zeros of J_0
array([2.40482556, 5.52007811, 8.65372791])

pytcl.mathematical_functions.special_functions.bessel.kelvin(x)[source]

Kelvin functions ber, bei, ker, kei.

Kelvin functions are the real and imaginary parts of the Bessel functions with argument x*exp(3*pi*i/4).

Parameters:

x (array_like) – Argument of the Kelvin functions.

Returns:

ber (ndarray) – Kelvin function ber(x).
bei (ndarray) – Kelvin function bei(x).
ker (ndarray) – Kelvin function ker(x).
kei (ndarray) – Kelvin function kei(x).

Return type:

tuple[ndarray[Any, Any], ndarray[Any, Any], ndarray[Any, Any], ndarray[Any, Any]]

Notes

ber(x) + i*bei(x) = J_0(x * exp(3*pi*i/4)) ker(x) + i*kei(x) = K_0(x * exp(pi*i/4))

Examples

>>> ber, bei, ker, kei = kelvin(1)
>>> round(float(ber), 6)
0.984382

Gamma Functions

Gamma and related functions.

This module provides gamma functions, factorials, and related special functions used in statistics and probability calculations.

pytcl.mathematical_functions.special_functions.gamma_functions.gamma(x)[source]

Gamma function.

Computes Γ(x) = ∫_0^∞ t^(x-1) * e^(-t) dt.

Parameters:: x (array_like) – Argument of the gamma function.
Returns:: Γ – Values of Γ(x).
Return type:: ndarray

Notes

For positive integers, Γ(n) = (n-1)!

Examples

>>> float(gamma(5))  # 4! = 24
24.0
>>> round(float(gamma(0.5)), 6)  # sqrt(pi)
1.772454

See also

scipy.special.gamma: Gamma function.

pytcl.mathematical_functions.special_functions.gamma_functions.gammaln(x)[source]

Natural logarithm of the absolute value of the gamma function.

Computes ln|Γ(x)|. This is more numerically stable than computing log(gamma(x)) for large x.

Parameters:: x (array_like) – Argument of the function.
Returns:: lng – Values of ln|Γ(x)|.
Return type:: ndarray

Examples

>>> round(float(gammaln(100)), 6)  # log(99!)
359.134205

See also

scipy.special.gammaln: Log of gamma function.

pytcl.mathematical_functions.special_functions.gamma_functions.gammainc(a, x)[source]

Regularized lower incomplete gamma function.

Computes P(a, x) = γ(a, x) / Γ(a), where γ(a, x) = ∫_0^x t^(a-1) * e^(-t) dt.

Parameters:

a (array_like) – Parameter of the function (must be positive).
x (array_like) – Upper limit of integration (must be non-negative).

Returns:

P – Values of the regularized lower incomplete gamma function.

Return type:

ndarray

Notes

This is the CDF of the gamma distribution.

Examples

>>> round(float(gammainc(1, 1)), 6)  # 1 - exp(-1)
0.632121

See also

scipy.special.gammainc: Regularized lower incomplete gamma function.
gammaincc: Upper incomplete gamma (complement).

pytcl.mathematical_functions.special_functions.gamma_functions.gammaincc(a, x)[source]

Regularized upper incomplete gamma function.

Computes Q(a, x) = Γ(a, x) / Γ(a) = 1 - P(a, x).

Parameters:

a (array_like) – Parameter of the function (must be positive).
x (array_like) – Lower limit of integration (must be non-negative).

Returns:

Q – Values of the regularized upper incomplete gamma function.

Return type:

ndarray

Examples

>>> round(float(gammaincc(1, 1)), 6)  # exp(-1)
0.367879

See also

scipy.special.gammaincc: Regularized upper incomplete gamma function.

pytcl.mathematical_functions.special_functions.gamma_functions.gammaincinv(a, y)[source]

Inverse of the regularized lower incomplete gamma function.

Finds x such that P(a, x) = y.

Parameters:

a (array_like) – Parameter of the function.
y (array_like) – Target probability (between 0 and 1).

Returns:

x – Values where P(a, x) = y.

Return type:

ndarray

Examples

>>> round(float(gammaincinv(1, 0.5)), 6)  # Median of exponential distribution
0.693147

See also

scipy.special.gammaincinv: Inverse of lower incomplete gamma.

pytcl.mathematical_functions.special_functions.gamma_functions.digamma(x)[source]

Digamma (psi) function.

Computes ψ(x) = d/dx ln(Γ(x)) = Γ’(x) / Γ(x).

Parameters:: x (array_like) – Argument of the function.
Returns:: ψ – Values of the digamma function.
Return type:: ndarray

Examples

>>> round(float(digamma(1)), 6)  # -γ (negative Euler-Mascheroni constant)
-0.577216

See also

scipy.special.digamma: Digamma function.
polygamma: Higher derivatives.

pytcl.mathematical_functions.special_functions.gamma_functions.polygamma(n, x)[source]

Polygamma function.

Computes ψ^(n)(x) = d^(n+1)/dx^(n+1) ln(Γ(x)).

Parameters:

n (int) – Order of the derivative (n=0 gives digamma, n=1 gives trigamma, etc.).
x (array_like) – Argument of the function.

Returns:

ψn – Values of the n-th polygamma function.

Return type:

ndarray

Examples

>>> round(float(polygamma(0, 1)), 6)  # Digamma at 1 = -γ
-0.577216
>>> round(float(polygamma(1, 1)), 6)  # Trigamma at 1 = π²/6
1.644934

See also

scipy.special.polygamma: Polygamma function.

pytcl.mathematical_functions.special_functions.gamma_functions.beta(a, b)[source]

Beta function.

Computes B(a, b) = Γ(a) * Γ(b) / Γ(a + b).

Parameters:

a (array_like) – First parameter.
b (array_like) – Second parameter.

Returns:

B – Values of the beta function.

Return type:

ndarray

Examples

>>> float(beta(1, 1))
1.0
>>> round(float(beta(0.5, 0.5)), 6)  # pi
3.141593

See also

scipy.special.beta: Beta function.

pytcl.mathematical_functions.special_functions.gamma_functions.betaln(a, b)[source]

Natural logarithm of the beta function.

Computes ln(B(a, b)) = ln(Γ(a)) + ln(Γ(b)) - ln(Γ(a + b)).

Parameters:

a (array_like) – First parameter.
b (array_like) – Second parameter.

Returns:

lnB – Values of ln(B(a, b)).

Return type:

ndarray

Examples

>>> import numpy as np
>>> round(float(betaln(100, 100)), 6)  # More stable than log(beta(100, 100))
-139.665259

See also

scipy.special.betaln: Log of beta function.

pytcl.mathematical_functions.special_functions.gamma_functions.betainc(a, b, x)[source]

Regularized incomplete beta function.

Computes I_x(a, b) = B(x; a, b) / B(a, b), where B(x; a, b) = ∫_0^x t^(a-1) * (1-t)^(b-1) dt.

Parameters:

a (array_like) – First parameter (must be positive).
b (array_like) – Second parameter (must be positive).
x (array_like) – Upper limit of integration (between 0 and 1).

Returns:

I – Values of the regularized incomplete beta function.

Return type:

ndarray

Notes

This is the CDF of the beta distribution.

Examples

>>> round(float(betainc(1, 1, 0.5)), 6)  # Uniform distribution CDF at 0.5
0.5

See also

scipy.special.betainc: Regularized incomplete beta function.

pytcl.mathematical_functions.special_functions.gamma_functions.betaincinv(a, b, y)[source]

Inverse of the regularized incomplete beta function.

Finds x such that I_x(a, b) = y.

Parameters:

a (array_like) – First parameter.
b (array_like) – Second parameter.
y (array_like) – Target probability (between 0 and 1).

Returns:

x – Values where I_x(a, b) = y.

Return type:

ndarray

Examples

>>> round(float(betaincinv(1, 1, 0.5)), 6)  # Median of uniform distribution
0.5

See also

scipy.special.betaincinv: Inverse of incomplete beta function.

pytcl.mathematical_functions.special_functions.gamma_functions.factorial(n, exact=False)[source]

Factorial function.

Computes n! = n * (n-1) * … * 2 * 1.

Parameters:

n (array_like) – Input values (non-negative integers).
exact (bool, optional) – If True, compute exact integer factorial (may overflow for large n). If False (default), use gamma function approximation.

Returns:

nfact – Values of n!.

Return type:

ndarray

Examples

>>> float(factorial(5))
120.0
>>> factorial(np.array([1, 2, 3, 4, 5]))
array([  1.,   2.,   6.,  24., 120.])

See also

scipy.special.factorial: Factorial function.

pytcl.mathematical_functions.special_functions.gamma_functions.factorial2(n, exact=False)[source]

Double factorial.

Computes n!! = n * (n-2) * (n-4) * … * (2 or 1).

Parameters:

n (array_like) – Input values (non-negative integers).
exact (bool, optional) – If True, compute exact integer result.

Returns:

nfact2 – Values of n!!.

Return type:

ndarray

Examples

>>> round(float(factorial2(5)), 6)  # 5 * 3 * 1 = 15
15.0
>>> round(float(factorial2(6)), 6)  # 6 * 4 * 2 = 48
48.0

See also

scipy.special.factorial2: Double factorial.

pytcl.mathematical_functions.special_functions.gamma_functions.comb(n, k, exact=False, repetition=False)[source]

Binomial coefficient (combinations).

Computes C(n, k) = n! / (k! * (n-k)!).

Parameters:

n (array_like) – Number of elements to choose from.
k (array_like) – Number of elements to choose.
exact (bool, optional) – If True, compute exact integer result.
repetition (bool, optional) – If True, compute combinations with repetition.

Returns:

C – Values of C(n, k).

Return type:

ndarray

Examples

>>> float(comb(5, 2))
10.0
>>> float(comb(10, 3))
120.0

See also

scipy.special.comb: Combinations.

pytcl.mathematical_functions.special_functions.gamma_functions.perm(n, k, exact=False)[source]

Permutation coefficient.

Computes P(n, k) = n! / (n-k)!.

Parameters:

n (array_like) – Number of elements to arrange.
k (array_like) – Number of elements in arrangement.
exact (bool, optional) – If True, compute exact integer result.

Returns:

P – Values of P(n, k).

Return type:

ndarray

Examples

>>> float(perm(5, 2))
20.0

See also

scipy.special.perm: Permutations.

Error Functions

Error functions and related special functions.

This module provides error functions and their variants, commonly used in probability theory and statistical analysis.

pytcl.mathematical_functions.special_functions.error_functions.erf(x)[source]

Error function.

Computes erf(x) = (2/√π) * ∫_0^x e^(-t²) dt.

Parameters:: x (array_like) – Argument of the error function.
Returns:: y – Values of erf(x).
Return type:: ndarray

Notes

erf(0) = 0
erf(∞) = 1
erf(-x) = -erf(x)

The error function is related to the normal distribution CDF by: Φ(x) = (1 + erf(x/√2)) / 2

Examples

>>> float(erf(0))
0.0
>>> round(float(erf(1)), 6)
0.842701

See also

scipy.special.erf: Error function.
erfc: Complementary error function.

pytcl.mathematical_functions.special_functions.error_functions.erfc(x)[source]

Complementary error function.

Computes erfc(x) = 1 - erf(x) = (2/√π) * ∫_x^∞ e^(-t²) dt.

Parameters:: x (array_like) – Argument of the function.
Returns:: y – Values of erfc(x).
Return type:: ndarray

Notes

This function is more accurate than computing 1 - erf(x) for large x.

Examples

>>> float(erfc(0))
1.0
>>> round(float(erfc(3)), 10)  # Very small
2.20905e-05

See also

scipy.special.erfc: Complementary error function.

pytcl.mathematical_functions.special_functions.error_functions.erfcx(x)[source]

Scaled complementary error function.

Computes erfcx(x) = exp(x²) * erfc(x).

Parameters:: x (array_like) – Argument of the function.
Returns:: y – Values of erfcx(x).
Return type:: ndarray

Notes

This function is useful when erfc(x) underflows but the scaled version remains representable.

Examples

>>> float(erfcx(0))
1.0
>>> round(float(erfcx(10)), 6)  # Remains finite even when erfc(10) underflows
0.056141

See also

scipy.special.erfcx: Scaled complementary error function.

pytcl.mathematical_functions.special_functions.error_functions.erfi(x)[source]

Imaginary error function.

Computes erfi(x) = -i * erf(i*x) = (2/√π) * ∫_0^x e^(t²) dt.

Parameters:: x (array_like) – Argument of the function.
Returns:: y – Values of erfi(x).
Return type:: ndarray

Examples

>>> float(erfi(0))
0.0
>>> round(float(erfi(1)), 6)
1.650426

See also

scipy.special.erfi: Imaginary error function.

pytcl.mathematical_functions.special_functions.error_functions.erfinv(y)[source]

Inverse error function.

Finds x such that erf(x) = y.

Parameters:: y (array_like) – Values in the range (-1, 1).
Returns:: x – Inverse error function values.
Return type:: ndarray

Examples

>>> float(erfinv(0))
0.0
>>> round(float(erf(erfinv(0.5))), 6)
0.5

See also

scipy.special.erfinv: Inverse error function.

pytcl.mathematical_functions.special_functions.error_functions.erfcinv(y)[source]

Inverse complementary error function.

Finds x such that erfc(x) = y.

Parameters:: y (array_like) – Values in the range (0, 2).
Returns:: x – Inverse complementary error function values.
Return type:: ndarray

Examples

>>> abs(float(erfcinv(1)))
0.0
>>> round(float(erfc(erfcinv(0.5))), 6)
0.5

See also

scipy.special.erfcinv: Inverse complementary error function.

pytcl.mathematical_functions.special_functions.error_functions.dawsn(x)[source]

Dawson’s integral.

Computes F(x) = exp(-x²) * ∫_0^x exp(t²) dt.

Parameters:: x (array_like) – Argument of Dawson’s integral.
Returns:: F – Values of Dawson’s integral.
Return type:: ndarray

Notes

Dawson’s integral is related to the imaginary error function by: F(x) = (√π/2) * exp(-x²) * erfi(x)

Examples

>>> float(dawsn(0))
0.0
>>> round(float(dawsn(1)), 6)
0.53808

See also

scipy.special.dawsn: Dawson’s integral.

pytcl.mathematical_functions.special_functions.error_functions.fresnel(x)[source]

Fresnel integrals.

Computes the Fresnel sine and cosine integrals: S(x) = ∫_0^x sin(π*t²/2) dt C(x) = ∫_0^x cos(π*t²/2) dt

Parameters:

x (array_like) – Argument of the Fresnel integrals.

Returns:

S (ndarray) – Fresnel sine integral.
C (ndarray) – Fresnel cosine integral.

Return type:

tuple[ndarray[Any, Any], ndarray[Any, Any]]

Examples

>>> S, C = fresnel(1)
>>> round(float(S), 6)
0.438259
>>> round(float(C), 6)
0.779893

See also

scipy.special.fresnel: Fresnel integrals.

pytcl.mathematical_functions.special_functions.error_functions.wofz(z)[source]

Faddeeva function.

Computes w(z) = exp(-z²) * erfc(-i*z).

Parameters:: z (array_like) – Argument (can be complex).
Returns:: w – Complex Faddeeva function values.
Return type:: ndarray

Notes

This function is useful in spectral line modeling and plasma physics.

Examples

>>> w = wofz(0)
>>> float(w.real)
1.0
>>> float(w.imag)
0.0

See also

scipy.special.wofz: Faddeeva function.

pytcl.mathematical_functions.special_functions.error_functions.voigt_profile(x, sigma, gamma)[source]

Voigt profile.

The Voigt profile is a convolution of a Gaussian and Lorentzian profile, commonly used in spectroscopy and line shape analysis.

Parameters:

x (array_like) – Position parameter.
sigma (float) – Standard deviation of the Gaussian component.
gamma (float) – Half-width at half-maximum of the Lorentzian component.

Returns:

V – Voigt profile values (normalized to unit area).

Return type:

ndarray

Examples

>>> round(float(voigt_profile(0, 1, 0)), 6)  # Pure Gaussian at x=0
0.398942

See also

scipy.special.voigt_profile: Voigt profile.

Elliptic Functions

Elliptic integrals and functions.

This module provides elliptic integrals used in various physical applications including orbits, pendulums, and electromagnetic calculations.

pytcl.mathematical_functions.special_functions.elliptic.ellipk(m)[source]

Complete elliptic integral of the first kind.

Computes K(m) = ∫_0^(π/2) (1 - m*sin²(θ))^(-1/2) dθ.

Parameters:: m (array_like) – Parameter m (not the modulus k). Note: m = k². Must be in [0, 1).
Returns:: K – Values of the complete elliptic integral of the first kind.
Return type:: ndarray

Notes

As m → 1, K(m) → ∞.

Examples

>>> round(float(ellipk(0)), 6)  # K(0) = π/2
1.570796
>>> round(float(ellipk(0.5)), 6)
1.854075

See also

scipy.special.ellipk: Complete elliptic integral of first kind.

pytcl.mathematical_functions.special_functions.elliptic.ellipkm1(p)[source]

Complete elliptic integral of the first kind around m = 1.

Computes K(1 - p) for small p, more accurate than ellipk(1 - p).

Parameters:: p (array_like) – Parameter p = 1 - m.
Returns:: K – Values of K(1 - p).
Return type:: ndarray

Examples

>>> round(float(ellipkm1(0.1)), 6)  # K(0.9)
2.578092

See also

scipy.special.ellipkm1: Elliptic integral near m = 1.

pytcl.mathematical_functions.special_functions.elliptic.ellipe(m)[source]

Complete elliptic integral of the second kind.

Computes E(m) = ∫_0^(π/2) (1 - m*sin²(θ))^(1/2) dθ.

Parameters:: m (array_like) – Parameter m (not the modulus k). Note: m = k². Must be in [0, 1].
Returns:: E – Values of the complete elliptic integral of the second kind.
Return type:: ndarray

Examples

>>> round(float(ellipe(0)), 6)  # E(0) = π/2
1.570796
>>> float(ellipe(1))  # E(1) = 1
1.0

See also

scipy.special.ellipe: Complete elliptic integral of second kind.

pytcl.mathematical_functions.special_functions.elliptic.ellipeinc(phi, m)[source]

Incomplete elliptic integral of the second kind.

Computes E(φ, m) = ∫_0^φ (1 - m*sin²(θ))^(1/2) dθ.

Parameters:

phi (array_like) – Amplitude (in radians).
m (array_like) – Parameter m = k².

Returns:

E – Values of the incomplete elliptic integral of the second kind.

Return type:

ndarray

Examples

>>> import numpy as np
>>> round(float(ellipeinc(np.pi/2, 0)), 6)  # Same as ellipe(0) = π/2
1.570796

See also

scipy.special.ellipeinc: Incomplete elliptic integral of second kind.

pytcl.mathematical_functions.special_functions.elliptic.ellipkinc(phi, m)[source]

Incomplete elliptic integral of the first kind.

Computes F(φ, m) = ∫_0^φ (1 - m*sin²(θ))^(-1/2) dθ.

Parameters:

phi (array_like) – Amplitude (in radians).
m (array_like) – Parameter m = k².

Returns:

F – Values of the incomplete elliptic integral of the first kind.

Return type:

ndarray

Examples

>>> import numpy as np
>>> round(float(ellipkinc(np.pi/2, 0)), 6)  # Same as ellipk(0) = π/2
1.570796

See also

scipy.special.ellipkinc: Incomplete elliptic integral of first kind.

pytcl.mathematical_functions.special_functions.elliptic.elliprd(x, y, z)[source]

Carlson symmetric elliptic integral R_D.

Computes the symmetric elliptic integral: R_D(x, y, z) = (3/2) ∫_0^∞ [(t+x)(t+y)]^(-1/2) (t+z)^(-3/2) dt

Parameters:

x (array_like) – First argument (non-negative).
y (array_like) – Second argument (non-negative).
z (array_like) – Third argument (positive).

Returns:

R_D – Values of the Carlson R_D integral.

Return type:

ndarray

Examples

>>> round(float(elliprd(1, 2, 3)), 6)
0.29046

See also

scipy.special.elliprd: Carlson R_D integral.

pytcl.mathematical_functions.special_functions.elliptic.elliprf(x, y, z)[source]

Carlson symmetric elliptic integral R_F.

Computes the symmetric elliptic integral: R_F(x, y, z) = (1/2) ∫_0^∞ [(t+x)(t+y)(t+z)]^(-1/2) dt

Parameters:

x (array_like) – First argument (non-negative).
y (array_like) – Second argument (non-negative).
z (array_like) – Third argument (non-negative). At most one of x, y, z can be zero.

Returns:

R_F – Values of the Carlson R_F integral.

Return type:

ndarray

Notes

The complete elliptic integral of the first kind is: K(m) = R_F(0, 1-m, 1)

Examples

>>> float(elliprf(1, 1, 1))  # R_F(a, a, a) = 1/sqrt(a)
1.0

See also

scipy.special.elliprf: Carlson R_F integral.

pytcl.mathematical_functions.special_functions.elliptic.elliprg(x, y, z)[source]

Carlson symmetric elliptic integral R_G.

Computes the symmetric elliptic integral R_G(x, y, z).

Parameters:

x (array_like) – First argument (non-negative).
y (array_like) – Second argument (non-negative).
z (array_like) – Third argument (non-negative).

Returns:

R_G – Values of the Carlson R_G integral.

Return type:

ndarray

Notes

The complete elliptic integral of the second kind is: E(m) = 2 * R_G(0, 1-m, 1)

Examples

>>> float(elliprg(1, 1, 1))  # R_G(a, a, a) = sqrt(a)
1.0

See also

scipy.special.elliprg: Carlson R_G integral.

pytcl.mathematical_functions.special_functions.elliptic.elliprj(x, y, z, p)[source]

Carlson symmetric elliptic integral R_J.

Computes the symmetric elliptic integral: R_J(x, y, z, p) = (3/2) ∫_0^∞ [(t+x)(t+y)(t+z)]^(-1/2) (t+p)^(-1) dt

Parameters:

x (array_like) – First argument (non-negative).
y (array_like) – Second argument (non-negative).
z (array_like) – Third argument (non-negative).
p (array_like) – Fourth argument (non-zero).

Returns:

R_J – Values of the Carlson R_J integral.

Return type:

ndarray

Notes

The complete elliptic integral of the third kind can be computed using R_J.

Examples

>>> round(float(elliprj(1, 2, 3, 4)), 6)
0.239848

See also

scipy.special.elliprj: Carlson R_J integral.

pytcl.mathematical_functions.special_functions.elliptic.elliprc(x, y)[source]

Carlson degenerate elliptic integral R_C.

Computes R_C(x, y) = R_F(x, y, y).

Parameters:

x (array_like) – First argument (non-negative).
y (array_like) – Second argument (non-zero).

Returns:

R_C – Values of the Carlson R_C integral.

Return type:

ndarray

Notes

R_C(x, y) = arctanh(sqrt((x-y)/x)) / sqrt(x-y) for x > y
R_C(x, y) = arctan(sqrt((y-x)/x)) / sqrt(y-x) for x < y

Examples

>>> float(elliprc(1, 1))  # R_C(a, a) = 1/sqrt(a)
1.0

See also

scipy.special.elliprc: Carlson R_C integral.

Statistics

Statistics and probability distributions.

This module provides: - Probability distribution classes with consistent APIs - Descriptive statistics (mean, variance, correlation) - Robust estimators (MAD, IQR) - Filter consistency metrics (NEES, NIS)

Distributions

Probability distributions.

This module provides probability distribution classes with consistent APIs for PDF, CDF, sampling, and moment calculations. These wrap scipy.stats distributions with additional functionality useful for tracking applications.

class pytcl.mathematical_functions.statistics.distributions.Distribution[source]

Bases: ABC

Abstract base class for probability distributions.

All distribution classes inherit from this and provide consistent methods for probability calculations.

abstractmethod pdf(x)[source]

Probability density function.

abstractmethod logpdf(x)[source]

Log of probability density function.

abstractmethod cdf(x)[source]

Cumulative distribution function.

abstractmethod ppf(q)[source]

Percent point function (inverse of CDF).

abstractmethod sample(size=None)[source]

Generate random samples.

abstractmethod mean()[source]

Distribution mean.

abstractmethod var()[source]

Distribution variance.

std()[source]

Distribution standard deviation.

class pytcl.mathematical_functions.statistics.distributions.Gaussian(mean=0.0, var=1.0)[source]

Bases: Distribution

Univariate Gaussian (Normal) distribution.

Parameters:

mean (float) – Mean of the distribution.
var (float) – Variance of the distribution.

Examples

>>> g = Gaussian(mean=0, var=1)
>>> round(float(g.pdf(0)), 6)
0.398942
>>> round(float(g.cdf(0)), 6)
0.5

__init__(mean=0.0, var=1.0)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.MultivariateGaussian(mean, cov)[source]

Bases: Distribution

Multivariate Gaussian (Normal) distribution.

Parameters:

mean (array_like) – Mean vector of shape (n,).
cov (array_like) – Covariance matrix of shape (n, n).

Examples

>>> mg = MultivariateGaussian(mean=[0, 0], cov=[[1, 0], [0, 1]])
>>> round(float(mg.pdf([0, 0])), 6)
0.159155

__init__(mean, cov)[source]

property dim: int: Dimension of the distribution.

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Return diagonal of covariance (marginal variances).

cov()[source]

Return full covariance matrix.

mahalanobis(x)[source]

Compute Mahalanobis distance from the mean.

Parameters:: x (array_like) – Point(s) to compute distance for.
Returns:: d – Mahalanobis distance(s).
Return type:: ndarray

class pytcl.mathematical_functions.statistics.distributions.Uniform(low=0.0, high=1.0)[source]

Bases: Distribution

Continuous uniform distribution.

Parameters:

low (float) – Lower bound of the distribution.
high (float) – Upper bound of the distribution.

__init__(low=0.0, high=1.0)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.Exponential(rate=1.0)[source]

Bases: Distribution

Exponential distribution.

Parameters:: rate (float) – Rate parameter (λ). Mean is 1/λ.

__init__(rate=1.0)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.Gamma(shape, rate=None, scale=None)[source]

Bases: Distribution

Gamma distribution.

Parameters:

shape (float) – Shape parameter (k or α).
rate (float, optional) – Rate parameter (β = 1/θ). Default is 1.
scale (float, optional) – Scale parameter (θ = 1/β). Alternative to rate.

Notes

Either rate or scale should be specified, not both.

__init__(shape, rate=None, scale=None)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.ChiSquared(df)[source]

Bases: Distribution

Chi-squared distribution.

Parameters:: df (int) – Degrees of freedom.

__init__(df)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.StudentT(df, loc=0.0, scale=1.0)[source]

Bases: Distribution

Student’s t-distribution.

Parameters:

df (float) – Degrees of freedom.
loc (float, optional) – Location parameter (default 0).
scale (float, optional) – Scale parameter (default 1).

__init__(df, loc=0.0, scale=1.0)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.Beta(a, b)[source]

Bases: Distribution

Beta distribution.

Parameters:

a (float) – First shape parameter (α > 0).
b (float) – Second shape parameter (β > 0).

__init__(a, b)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.Poisson(rate)[source]

Bases: Distribution

Poisson distribution (discrete).

Parameters:: rate (float) – Rate parameter (λ), also the mean.

__init__(rate)[source]

pdf(x)[source]

Probability mass function (PMF).

logpdf(x)[source]

Log of probability mass function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.VonMises(mu=0.0, kappa=1.0)[source]

Bases: Distribution

Von Mises distribution (circular normal).

Useful for angular/directional data in tracking applications.

Parameters:

mu (float) – Mean direction (in radians).
kappa (float) – Concentration parameter.

__init__(mu=0.0, kappa=1.0)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

class pytcl.mathematical_functions.statistics.distributions.Wishart(df, scale)[source]

Bases: Distribution

Wishart distribution (matrix-valued).

The Wishart distribution is used for covariance matrix estimation in multivariate statistics.

Parameters:

df (float) – Degrees of freedom.
scale (array_like) – Scale matrix (positive definite).

__init__(df, scale)[source]

pdf(x)[source]

Probability density function.

logpdf(x)[source]

Log of probability density function.

cdf(x)[source]

Cumulative distribution function.

ppf(q)[source]

Percent point function (inverse of CDF).

sample(size=None)[source]

Generate random samples.

mean()[source]

Distribution mean.

var()[source]

Distribution variance.

Estimators

Statistical estimators and descriptive statistics.

This module provides functions for computing sample statistics, robust estimators, and related quantities used in tracking applications.

pytcl.mathematical_functions.statistics.estimators.weighted_mean(x, weights, axis=None)[source]

Compute weighted mean.

Parameters:

x (array_like) – Input data.
weights (array_like) – Weights for each data point.
axis (int, optional) – Axis along which to compute. Default is None (all elements).

Returns:

mean – Weighted mean. A scalar when axis is None, otherwise an array.

Return type:

float or ndarray

Examples

>>> weighted_mean([1, 2, 3], [1, 1, 2])
2.25

pytcl.mathematical_functions.statistics.estimators.weighted_var(x, weights, ddof=0, axis=None)[source]

Compute weighted variance.

Parameters:

x (array_like) – Input data.
weights (array_like) – Weights for each data point.
ddof (int, optional) – Delta degrees of freedom. Default is 0 (population variance).
axis (int, optional) – Axis along which to compute.

Returns:

var – Weighted variance.

Return type:

ndarray

Examples

>>> x = [1, 2, 3]
>>> weights = [1, 1, 2]
>>> float(weighted_var(x, weights))
0.6875

pytcl.mathematical_functions.statistics.estimators.weighted_cov(x, weights, ddof=0)[source]

Compute weighted covariance matrix.

Parameters:

x (array_like) – Data matrix of shape (n_samples, n_features).
weights (array_like) – Weights of shape (n_samples,).
ddof (int, optional) – Delta degrees of freedom. Default is 0.

Returns:

cov – Weighted covariance matrix of shape (n_features, n_features).

Return type:

ndarray

Examples

>>> x = [[1, 2], [2, 3], [3, 4]]
>>> weights = [1, 1, 1]
>>> cov = weighted_cov(x, weights)
>>> cov.shape
(2, 2)

pytcl.mathematical_functions.statistics.estimators.sample_mean(x, axis=None)[source]

Compute sample mean.

Parameters:

x (array_like) – Input data.
axis (int, optional) – Axis along which to compute.

Returns:

mean – Sample mean.

Return type:

ndarray

pytcl.mathematical_functions.statistics.estimators.sample_var(x, ddof=1, axis=None)[source]

Compute sample variance.

Parameters:

x (array_like) – Input data.
ddof (int, optional) – Delta degrees of freedom. Default is 1 (unbiased estimator).
axis (int, optional) – Axis along which to compute.

Returns:

var – Sample variance.

Return type:

ndarray

Examples

>>> x = [1, 2, 3, 4, 5]
>>> sample_var(x)
2.5

pytcl.mathematical_functions.statistics.estimators.sample_cov(x, y=None, ddof=1)[source]

Compute sample covariance matrix.

Parameters:

x (array_like) – Data matrix of shape (n_samples, n_features) or 1D array.
y (array_like, optional) – Second variable for cross-covariance.
ddof (int, optional) – Delta degrees of freedom. Default is 1.

Returns:

cov – Covariance matrix, or a scalar variance when x is 1D and y is None.

Return type:

float or ndarray

Examples

>>> x = [[1, 2], [2, 3], [3, 4]]
>>> cov = sample_cov(x)
>>> cov.shape
(2, 2)

pytcl.mathematical_functions.statistics.estimators.sample_corr(x)[source]

Compute sample correlation matrix.

Parameters:: x (array_like) – Data matrix of shape (n_samples, n_features).
Returns:: corr – Correlation matrix of shape (n_features, n_features).
Return type:: ndarray

Examples

>>> x = [[1, 2], [2, 3], [3, 4]]
>>> corr = sample_corr(x)
>>> corr.shape
(2, 2)
>>> corr[0, 0]  # Correlation of feature 1 with itself
1.0

pytcl.mathematical_functions.statistics.estimators.median(x, axis=None)[source]

Compute median.

Parameters:

x (array_like) – Input data.
axis (int, optional) – Axis along which to compute.

Returns:

med – Median value(s).

Return type:

ndarray

Examples

>>> median([1, 2, 3, 4, 5])
3.0
>>> median([1, 2, 3, 4])
2.5

pytcl.mathematical_functions.statistics.estimators.mad(x, axis=None, scale=1.4826)[source]

Median Absolute Deviation (MAD).

A robust measure of statistical dispersion.

Parameters:

x (array_like) – Input data.
axis (int, optional) – Axis along which to compute.
scale (float, optional) – Scale factor for consistency with standard deviation for normal distributions. Default is 1.4826.

Returns:

mad – MAD value(s).

Return type:

ndarray

Examples

>>> mad([1, 2, 3, 4, 5])
1.4826

Notes

For normally distributed data, scale * MAD approximates the standard deviation.

pytcl.mathematical_functions.statistics.estimators.iqr(x, axis=None)[source]

Interquartile range (IQR).

Parameters:

x (array_like) – Input data.
axis (int, optional) – Axis along which to compute.

Returns:

iqr – Interquartile range (Q3 - Q1).

Return type:

ndarray

Examples

>>> float(iqr([1, 2, 3, 4, 5, 6, 7, 8, 9]))
4.0

pytcl.mathematical_functions.statistics.estimators.skewness(x, axis=None, bias=True)[source]

Compute sample skewness.

Parameters:

x (array_like) – Input data.
axis (int, optional) – Axis along which to compute.
bias (bool, optional) – If False, apply bias correction. Default is True.

Returns:

skew – Skewness value(s).

Return type:

ndarray

Examples

>>> float(skewness([1, 2, 3, 4, 5]))
0.0

pytcl.mathematical_functions.statistics.estimators.kurtosis(x, axis=None, fisher=True, bias=True)[source]

Compute sample kurtosis.

Parameters:

x (array_like) – Input data.
axis (int, optional) – Axis along which to compute.
fisher (bool, optional) – If True, return excess kurtosis (Fisher definition). If False, return Pearson kurtosis. Default is True.
bias (bool, optional) – If False, apply bias correction. Default is True.

Returns:

kurt – Kurtosis value(s).

Return type:

ndarray

Examples

>>> float(kurtosis([1, 2, 3, 4, 5]))
-1.3

pytcl.mathematical_functions.statistics.estimators.moment(x, order, axis=None, central=True)[source]

Compute sample moment.

Parameters:

x (array_like) – Input data.
order (int) – Order of the moment.
axis (int, optional) – Axis along which to compute.
central (bool, optional) – If True, compute central moment. Default is True.

Returns:

m – Moment value(s).

Return type:

ndarray

Examples

>>> float(moment([1, 2, 3, 4, 5], order=2))
2.0
>>> float(moment([1, 2, 3, 4, 5], order=2, central=False))
11.0

pytcl.mathematical_functions.statistics.estimators.nees(error, covariance)[source]

Normalized Estimation Error Squared (NEES).

A consistency metric for estimators. For a consistent estimator, NEES should be chi-squared distributed with n degrees of freedom.

Parameters:

error (array_like) – Estimation error vector(s) of shape (n,) or (m, n).
covariance (array_like) – Covariance matrix of shape (n, n).

Returns:

nees – NEES value(s). Scalar if error is 1D, array if error is 2D.

Return type:

ndarray

Examples

>>> error = np.array([1.0, 0.5])
>>> cov = np.array([[1, 0], [0, 1]])
>>> nees(error, cov)
1.25

pytcl.mathematical_functions.statistics.estimators.nis(innovation, innovation_covariance)[source]

Normalized Innovation Squared (NIS).

A filter consistency metric based on measurement innovations. For a consistent filter, NIS should be chi-squared distributed.

Parameters:

innovation (array_like) – Innovation (measurement residual) vector(s).
innovation_covariance (array_like) – Innovation covariance matrix.

Returns:

nis – NIS value(s).

Return type:

ndarray

Notes

This is equivalent to NEES applied to innovations.

Interpolation

Interpolation methods.

This module provides: - 1D interpolation (linear, spline, PCHIP, Akima) - 2D/3D interpolation on regular grids - RBF interpolation for scattered data - Spherical interpolation

pytcl.mathematical_functions.interpolation.interp1d(x, y, kind='linear', fill_value=nan, bounds_error=False)[source]

Create a 1D interpolation function.

Parameters:

x (array_like) – Sample points (must be monotonically increasing).
y (array_like) – Sample values.
kind (str, optional) – Interpolation method: - ‘linear’: Linear interpolation (default) - ‘nearest’: Nearest neighbor - ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’: Spline of order 0, 1, 2, 3 - ‘previous’, ‘next’: Previous/next value
fill_value (float, tuple, or 'extrapolate', optional) – Value for points outside data range. Default is NaN. Use ‘extrapolate’ to extrapolate beyond bounds.
bounds_error (bool, optional) – If True, raise error for out-of-bounds. Default is False.

Returns:

f – Interpolation function that takes x values and returns y values.

Return type:

callable

Examples

>>> x = np.array([0, 1, 2, 3])
>>> y = np.array([0, 1, 4, 9])
>>> f = interp1d(x, y, kind='quadratic')
>>> f(1.5)
array(2.25)

See also

scipy.interpolate.interp1d: Underlying implementation.

pytcl.mathematical_functions.interpolation.linear_interp(x, xp, fp, left=None, right=None)[source]

One-dimensional linear interpolation.

Parameters:

x (array_like) – X-coordinates at which to evaluate.
xp (array_like) – X-coordinates of data points (must be increasing).
fp (array_like) – Y-coordinates of data points.
left (float, optional) – Value for x < xp[0]. Default is fp[0].
right (float, optional) – Value for x > xp[-1]. Default is fp[-1].

Returns:

y – Interpolated values.

Return type:

ndarray

Examples

>>> float(linear_interp(2.5, [1, 2, 3], [1, 4, 9]))
6.5

See also

numpy.interp: Underlying implementation.

pytcl.mathematical_functions.interpolation.cubic_spline(x, y, bc_type='not-a-knot')[source]

Create a cubic spline interpolation.

Parameters:

x (array_like) – Sample points (must be strictly increasing).
y (array_like) – Sample values.
bc_type (str, optional) – Boundary condition type: - ‘not-a-knot’: Default, uses continuity conditions. - ‘clamped’: First derivatives at endpoints are zero. - ‘natural’: Second derivatives at endpoints are zero. - ‘periodic’: Periodic boundary conditions.

Returns:

cs – Cubic spline object. Call cs(x_new) to interpolate.

Return type:

CubicSpline

Examples

>>> x = np.linspace(0, 2*np.pi, 10)
>>> y = np.sin(x)
>>> cs = cubic_spline(x, y)
>>> round(float(cs(np.pi/2)), 6)
0.999912

See also

scipy.interpolate.CubicSpline: Underlying implementation.

pytcl.mathematical_functions.interpolation.pchip(x, y)[source]

Piecewise Cubic Hermite Interpolating Polynomial (PCHIP).

PCHIP preserves monotonicity and avoids overshooting, making it suitable for data that should not have spurious oscillations.

Parameters:

x (array_like) – Sample points (must be strictly increasing).
y (array_like) – Sample values.

Returns:

p – PCHIP interpolator object.

Return type:

PchipInterpolator

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.interpolation import pchip
>>> # Monotonic data: stock prices (should not overshoot)
>>> x = np.array([0, 1, 2, 3, 4])
>>> y = np.array([10, 12, 11, 15, 18])  # Non-monotonic but generally increasing
>>> p = pchip(x, y)
>>> # Evaluate at intermediate points
>>> x_new = np.array([0.5, 1.5, 2.5])
>>> y_new = p(x_new)
>>> # PCHIP preserves bounds: should stay within observed values
>>> np.all((y_new >= y.min()) & (y_new <= y.max()))
True

Notes

Unlike cubic splines, PCHIP will not overshoot if the data is monotonic, making it more suitable for physical quantities that must stay positive or bounded.

See also

scipy.interpolate.PchipInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.akima(x, y)[source]

Akima interpolation.

Akima interpolation is a smooth interpolation method that avoids excessive oscillation compared to cubic splines.

Parameters:

x (array_like) – Sample points (must be strictly increasing).
y (array_like) – Sample values.

Returns:

a – Akima interpolator object.

Return type:

Akima1DInterpolator

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.interpolation import akima
>>> # Noisy data where smoothness without oscillation is desired
>>> x = np.array([0, 1, 2, 3, 4, 5])
>>> y = np.array([0, 1, 1.5, 1.2, 2.0, 2.5])  # Noisy measurements
>>> a = akima(x, y)
>>> # Evaluate at intermediate points
>>> x_new = np.array([0.5, 1.5, 3.5])
>>> y_new = a(x_new)
>>> # Akima should produce smooth, non-oscillating results
>>> y_new.shape
(3,)

See also

scipy.interpolate.Akima1DInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.interp2d(x, y, z, kind='linear')[source]

Create a 2D interpolation function on a regular grid.

Parameters:

x (array_like) – Grid coordinates along first axis.
y (array_like) – Grid coordinates along second axis.
z (array_like) – Values on the grid of shape (len(x), len(y)).
kind (str, optional) – Interpolation method: ‘linear’, ‘cubic’, or ‘quintic’. Default is ‘linear’.

Returns:

f – Interpolation function. Call f((xi, yi)) to interpolate.

Return type:

RegularGridInterpolator

Examples

>>> x = np.linspace(0, 4, 5)
>>> y = np.linspace(0, 4, 5)
>>> z = np.outer(x, y)  # z = x * y
>>> f = interp2d(x, y, z)
>>> f([[2.5, 2.5]])
array([6.25])

See also

scipy.interpolate.RegularGridInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.interp3d(x, y, z, values, kind='linear')[source]

Create a 3D interpolation function on a regular grid.

Parameters:

x (array_like) – Grid coordinates along first axis.
y (array_like) – Grid coordinates along second axis.
z (array_like) – Grid coordinates along third axis.
values (array_like) – Values on the grid of shape (len(x), len(y), len(z)).
kind (str, optional) – Interpolation method: ‘linear’ or ‘nearest’. Default is ‘linear’.

Returns:

f – Interpolation function.

Return type:

RegularGridInterpolator

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.interpolation import interp3d
>>> # Create a 3D grid: temperature field
>>> x = np.array([0, 1, 2])
>>> y = np.array([0, 1, 2])
>>> z = np.array([0, 1])
>>> # Temperature values at grid points (3x3x2)
>>> values = np.arange(18).reshape((3, 3, 2), order='C').astype(float)
>>> f = interp3d(x, y, z, values, kind='linear')
>>> # Interpolate at intermediate points
>>> pts = np.array([[0.5, 0.5, 0.5], [1.5, 1.5, 0.5]])
>>> result = f(pts)
>>> result.shape
(2,)

See also

scipy.interpolate.RegularGridInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.rbf_interpolate(points, values, kernel='thin_plate_spline', smoothing=0.0)[source]

Radial Basis Function (RBF) interpolation.

RBF interpolation works with scattered (non-grid) data in any dimension.

Parameters:

points (array_like) – Data point coordinates of shape (n_samples, n_dims).
values (array_like) – Values at data points of shape (n_samples,) or (n_samples, n_values).
kernel (str, optional) – RBF kernel function. Default is ‘thin_plate_spline’.
smoothing (float, optional) – Smoothing parameter. 0 means exact interpolation. Default is 0.

Returns:

rbf – RBF interpolation object.

Return type:

RBFInterpolator

Examples

>>> points = np.array([[0, 0], [1, 0], [0, 1], [1, 1]])
>>> values = np.array([0, 1, 1, 2])
>>> rbf = rbf_interpolate(points, values)
>>> rbf([[0.5, 0.5]])
array([1.])

See also

scipy.interpolate.RBFInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.barycentric(x, y)[source]

Barycentric polynomial interpolation.

This is a numerically stable method for polynomial interpolation.

Parameters:

x (array_like) – Sample points.
y (array_like) – Sample values.

Returns:

p – Barycentric interpolator object.

Return type:

BarycentricInterpolator

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.interpolation import barycentric
>>> # Interpolate a polynomial at Chebyshev nodes (most stable)
>>> n = 5
>>> x = np.cos(np.linspace(0, np.pi, n))  # Chebyshev nodes
>>> y = x**2 + 2*x + 1  # Quadratic function
>>> poly = barycentric(x, y)
>>> # Evaluate interpolant at new points
>>> x_new = np.linspace(-0.8, 0.8, 3)
>>> y_interp = poly(x_new)
>>> # Should match the original function well
>>> y_exact = x_new**2 + 2*x_new + 1
>>> np.allclose(y_interp, y_exact, atol=0.01)
True

See also

scipy.interpolate.BarycentricInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.krogh(x, y)[source]

Krogh interpolation.

Polynomial interpolation using divided differences.

Parameters:

x (array_like) – Sample points.
y (array_like) – Sample values (can include derivatives at points).

Returns:

k – Krogh interpolator object.

Return type:

KroghInterpolator

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.interpolation import krogh
>>> # Hermite interpolation: repeat a sample point to specify its derivative
>>> x = np.array([0.0, 0.0, 1.0, 1.0])
>>> y = np.array([1.0, 0.0, 2.0, 1.0])  # f(0)=1, f'(0)=0, f(1)=2, f'(1)=1
>>> k = krogh(x, y)
>>> # Interpolant passes through the specified function values
>>> float(k(0.0)), float(k(1.0))
(1.0, 2.0)

See also

scipy.interpolate.KroghInterpolator: Underlying implementation.

pytcl.mathematical_functions.interpolation.spherical_interp(lat, lon, values)[source]

Interpolation on a spherical surface.

Converts lat/lon to 3D Cartesian coordinates and uses RBF interpolation.

Parameters:

lat (array_like) – Latitude in radians of shape (n_samples,).
lon (array_like) – Longitude in radians of shape (n_samples,).
values (array_like) – Values at sample points.

Returns:

interp – Interpolation function. Call with 3D Cartesian coordinates.

Return type:

RBFInterpolator

Notes

To interpolate at new lat/lon points:

1. Convert lat/lon to Cartesian: x=cos(lat)*cos(lon),
   y=cos(lat)*sin(lon), z=sin(lat)
2. Call interp([[x, y, z]])

Numerical Integration

Numerical integration (quadrature) methods.

This module provides: - Gaussian quadrature rules (Legendre, Hermite, Laguerre, Chebyshev) - Adaptive integration functions - Multi-dimensional cubature rules for filtering (CKF, UKF)

pytcl.mathematical_functions.numerical_integration.gauss_legendre(n)[source]

Gauss-Legendre quadrature points and weights.

For integrating f(x) over [-1, 1]: ∫_{-1}^{1} f(x) dx ≈ Σ w_i * f(x_i)

Parameters:

n (int) – Number of quadrature points.

Returns:

x (ndarray) – Quadrature points of shape (n,).
w (ndarray) – Quadrature weights of shape (n,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x, w = gauss_legendre(5)
>>> # Integrate x^2 from -1 to 1 (exact = 2/3)
>>> round(float(np.sum(w * x**2)), 6)
0.666667

See also

numpy.polynomial.legendre.leggauss: Equivalent function.

pytcl.mathematical_functions.numerical_integration.gauss_hermite(n)[source]

Gauss-Hermite quadrature points and weights.

For integrating f(x) * exp(-x^2) over (-∞, ∞): ∫_{-∞}^{∞} f(x) * exp(-x²) dx ≈ Σ w_i * f(x_i)

Useful for expectations over Gaussian distributions.

Parameters:

n (int) – Number of quadrature points.

Returns:

x (ndarray) – Quadrature points of shape (n,).
w (ndarray) – Quadrature weights of shape (n,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x, w = gauss_hermite(5)
>>> # Compute E[X^2] for X ~ N(0, 1) (exact = 1)
>>> result = np.sum(w * (np.sqrt(2) * x)**2) / np.sqrt(np.pi)
>>> abs(result - 1.0) < 1e-10
True

Notes

For computing E[f(X)] where X ~ N(μ, σ²): E[f(X)] = (1/√π) * Σ w_i * f(μ + √2 * σ * x_i)

See also

numpy.polynomial.hermite.hermgauss: Equivalent function.

pytcl.mathematical_functions.numerical_integration.gauss_laguerre(n)[source]

Gauss-Laguerre quadrature points and weights.

For integrating f(x) * exp(-x) over [0, ∞): ∫_0^∞ f(x) * exp(-x) dx ≈ Σ w_i * f(x_i)

Parameters:

n (int) – Number of quadrature points.

Returns:

x (ndarray) – Quadrature points of shape (n,).
w (ndarray) – Quadrature weights of shape (n,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x, w = gauss_laguerre(5)
>>> # Integrate x * exp(-x) from 0 to inf (exact = 1)
>>> round(float(np.sum(w * x)), 10)
1.0

See also

numpy.polynomial.laguerre.laggauss: Equivalent function.

pytcl.mathematical_functions.numerical_integration.gauss_chebyshev(n, kind=1)[source]

Gauss-Chebyshev quadrature points and weights.

For kind=1, integrates f(x) / sqrt(1-x²) over [-1, 1]. For kind=2, integrates f(x) * sqrt(1-x²) over [-1, 1].

Parameters:

n (int) – Number of quadrature points.
kind ({1, 2}, optional) – Type of Chebyshev polynomial. Default is 1.

Returns:

x (ndarray) – Quadrature points of shape (n,).
w (ndarray) – Quadrature weights of shape (n,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x, w = gauss_chebyshev(5, kind=1)
>>> x.shape
(5,)

See also

numpy.polynomial.chebyshev.chebgauss: Type 1 Chebyshev.

pytcl.mathematical_functions.numerical_integration.quad(f, a, b, **kwargs)[source]

Adaptive quadrature integration.

Computes ∫_a^b f(x) dx using adaptive Gaussian quadrature.

Parameters:

f (callable) – Function to integrate.
a (float) – Lower limit.
b (float) – Upper limit.
**kwargs (Any) – Additional arguments passed to scipy.integrate.quad.

Returns:

result (float) – Estimated integral value.
error (float) – Estimate of the absolute error.

Return type:

Tuple[float, float]

Examples

>>> result, error = quad(lambda x: x**2, 0, 1)
>>> round(result, 6)
0.333333

See also

scipy.integrate.quad: Underlying implementation.

pytcl.mathematical_functions.numerical_integration.dblquad(f, a, b, gfun, hfun, **kwargs)[source]

Double integration.

Computes ∫_a^b ∫_{g(x)}^{h(x)} f(y, x) dy dx.

Parameters:

f (callable) – Function f(y, x) to integrate.
a (float) – Lower limit of x.
b (float) – Upper limit of x.
gfun (callable) – Lower limit of y as function of x.
hfun (callable) – Upper limit of y as function of x.
**kwargs (Any) – Additional arguments passed to scipy.integrate.dblquad.

Returns:

result (float) – Estimated integral value.
error (float) – Estimate of the absolute error.

Return type:

Tuple[float, float]

Examples

>>> # Integrate x*y over unit square
>>> result, error = dblquad(lambda y, x: x*y, 0, 1, lambda x: 0, lambda x: 1)
>>> round(result, 10)  # Should be 0.25
0.25

See also

scipy.integrate.dblquad: Underlying implementation.

pytcl.mathematical_functions.numerical_integration.tplquad(f, a, b, gfun, hfun, qfun, rfun, **kwargs)[source]

Triple integration.

Computes ∫_a^b ∫_{g(x)}^{h(x)} ∫_{q(x,y)}^{r(x,y)} f(z, y, x) dz dy dx.

Parameters:

f (callable) – Function f(z, y, x) to integrate.
a (float) – Lower limit of x.
b (float) – Upper limit of x.
gfun (callable) – Lower limit of y as function of x.
hfun (callable) – Upper limit of y as function of x.
qfun (callable) – Lower limit of z as function of x, y.
rfun (callable) – Upper limit of z as function of x, y.
**kwargs (Any) – Additional arguments passed to scipy.integrate.tplquad.

Returns:

result (float) – Estimated integral value.
error (float) – Estimate of the absolute error.

Return type:

Tuple[float, float]

Examples

>>> # Integrate x*y*z over unit cube
>>> result, error = tplquad(
...     lambda z, y, x: x*y*z,
...     0, 1,
...     lambda x: 0, lambda x: 1,
...     lambda x, y: 0, lambda x, y: 1
... )
>>> abs(result - 0.125) < 1e-6  # 1/8
True

See also

scipy.integrate.tplquad: Underlying implementation.

pytcl.mathematical_functions.numerical_integration.fixed_quad(f, a, b, n=5)[source]

Fixed-order Gaussian quadrature.

Computes ∫_a^b f(x) dx using n-point Gauss-Legendre quadrature.

Parameters:

f (callable) – Function to integrate. Should accept and return arrays.
a (float) – Lower limit.
b (float) – Upper limit.
n (int, optional) – Number of quadrature points. Default is 5.

Returns:

result (float) – Estimated integral value.
None – Placeholder for compatibility (no error estimate).

Return type:

tuple[float, None]

Examples

>>> result, _ = fixed_quad(lambda x: x**2, 0, 1, n=5)
>>> round(result, 6)
0.333333

See also

scipy.integrate.fixed_quad: Underlying implementation.

pytcl.mathematical_functions.numerical_integration.romberg(f, a, b, tol=1e-08, max_steps=20)[source]

Romberg integration.

Uses Richardson extrapolation to accelerate the trapezoidal rule.

Parameters:

f (callable) – Function to integrate.
a (float) – Lower limit.
b (float) – Upper limit.
tol (float, optional) – Desired tolerance. Default is 1e-8.
max_steps (int, optional) – Maximum number of extrapolation steps. Default is 20.

Returns:

result – Estimated integral value.

Return type:

float

Examples

>>> # Integrate x^2 from 0 to 1 (exact = 1/3)
>>> result = romberg(lambda x: x**2, 0, 1)
>>> abs(result - 1/3) < 1e-8
True

Notes

This is a native implementation that does not depend on scipy.integrate.romberg, which was deprecated in scipy 1.12 and removed in scipy 1.15.

pytcl.mathematical_functions.numerical_integration.simpson(y, x=None, dx=1.0)[source]

Simpson’s rule integration from samples.

Parameters:

y (array_like) – Array of function values.
x (array_like, optional) – Sample points. If None, uses uniform spacing dx.
dx (float, optional) – Spacing between samples if x is None. Default is 1.

Returns:

result – Estimated integral.

Return type:

float

Examples

>>> import numpy as np
>>> x = np.linspace(0, np.pi, 101)
>>> y = np.sin(x)
>>> result = simpson(y, x)  # Should be ~2.0
>>> abs(result - 2.0) < 1e-5
True

See also

scipy.integrate.simpson: Underlying implementation.

pytcl.mathematical_functions.numerical_integration.trapezoid(y, x=None, dx=1.0)[source]

Trapezoidal rule integration from samples.

Parameters:

y (array_like) – Array of function values.
x (array_like, optional) – Sample points. If None, uses uniform spacing dx.
dx (float, optional) – Spacing between samples if x is None. Default is 1.

Returns:

result – Estimated integral.

Return type:

float

Examples

>>> import numpy as np
>>> x = np.linspace(0, 1, 11)
>>> y = x**2  # Integrate x^2 from 0 to 1
>>> result = trapezoid(y, x)
>>> abs(result - 1/3) < 0.01  # Approximation
True

See also

scipy.integrate.trapezoid: Underlying implementation.

pytcl.mathematical_functions.numerical_integration.cubature_gauss_hermite(n_dim, n_points_per_dim)[source]

Tensor product Gauss-Hermite cubature rule.

Creates a multi-dimensional quadrature rule for integrating over a multivariate Gaussian distribution.

Parameters:

n_dim (int) – Number of dimensions.
n_points_per_dim (int) – Number of quadrature points per dimension.

Returns:

points (ndarray) – Cubature points of shape (n_points_per_dim^n_dim, n_dim).
weights (ndarray) – Cubature weights of shape (n_points_per_dim^n_dim,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

The number of points grows exponentially with dimension. For high dimensions, consider using sparse grid methods.

Examples

>>> points, weights = cubature_gauss_hermite(2, 3)
>>> points.shape
(9, 2)

pytcl.mathematical_functions.numerical_integration.spherical_cubature(n_dim)[source]

Spherical cubature rule for Gaussian integrals.

A 2n-point cubature rule that is exact for polynomials up to degree 3. This is the rule used in the Cubature Kalman Filter (CKF).

Parameters:

n_dim (int) – Number of dimensions.

Returns:

points (ndarray) – Cubature points of shape (2*n_dim, n_dim).
weights (ndarray) – Cubature weights of shape (2*n_dim,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

Points are at ±√n along each axis, scaled for use with standard normal distributions.

For computing E[f(X)] where X ~ N(μ, P): - Transform points: x_i = μ + chol(P) @ points[i] - E[f(X)] ≈ Σ weights[i] * f(x_i)

References

Arasaratnam & Haykin, “Cubature Kalman Filters”, IEEE TAC, 2009.

Examples

>>> points, weights = spherical_cubature(3)
>>> points.shape  # 2*n = 6 points in 3D
(6, 3)
>>> weights.shape
(6,)
>>> round(float(np.sum(weights)), 6)  # Weights sum to 1
1.0

pytcl.mathematical_functions.numerical_integration.unscented_transform_points(n_dim, alpha=0.001, beta=2.0, kappa=None)[source]

Generate sigma points and weights for unscented transform.

Parameters:

n_dim (int) – Number of dimensions.
alpha (float, optional) – Spread of sigma points. Default is 1e-3.
beta (float, optional) – Prior knowledge parameter (2 is optimal for Gaussian). Default is 2.
kappa (float, optional) – Secondary scaling parameter. Default is 3 - n_dim.

Returns:

sigma_points (ndarray) – Relative sigma point positions of shape (2*n_dim + 1, n_dim). Center point is at index 0, followed by ±directions.
wm (ndarray) – Weights for computing mean, shape (2*n_dim + 1,).
wc (ndarray) – Weights for computing covariance, shape (2*n_dim + 1,).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

For a random variable X ~ N(μ, P), the sigma points are: - χ_0 = μ - χ_i = μ + (√((n+λ)P))_i for i = 1..n - χ_{n+i} = μ - (√((n+λ)P))_i for i = 1..n

where (√A)_i is the i-th column of the matrix square root.

References

Julier & Uhlmann, “Unscented Filtering and Nonlinear Estimation”, Proc. IEEE, 2004.

See also

pytcl.mathematical_functions.numerical_integration.cubature_points.second_order_cubature_points: A structurally different (n+2)-point rule from the same paper’s Appendix III (spherical simplex sigma points), exact through degree 2 only (this function’s 2n+1 points are exact through degree 3). Not an alternate parameterization of this rule – see its Notes section for the reconciliation.

Examples

>>> sigma_points, wm, wc = unscented_transform_points(3)
>>> sigma_points.shape  # 2*n+1 = 7 points in 3D
(7, 3)
>>> wm.shape
(7,)
>>> np.abs(np.sum(wm) - 1.0) < 1e-10  # Mean weights sum to 1
True

pytcl.mathematical_functions.numerical_integration.second_order_cubature_points(n, w0=0.3333333333333333, alpha=1.0)[source]

Degree-2 cubature points for the standard normal N(0, I): the scaled unscented transformation’s minimal (n+2)-point spherical-simplex rule.

Counterpart of the MATLAB TCL’s secondOrderCubPoints: the scaled unscented transformation of Julier (2002), itself a scaled generalization of the spherical simplex sigma points of Appendix III of Julier and Uhlmann (2004) (alpha=1 reproduces those unscaled). A center point at the origin plus n+1 simplex vertices, built recursively one dimension at a time and then rescaled so the point set has unit sample covariance.

This is a genuinely different construction from both cubature rules already in pytcl, not an alternate parameterization of either – see Notes.

Parameters:

n (int) – Dimension, n >= 1.
w0 (float, optional) – Weight of the center point before alpha-scaling, 0 < w0 < 1. Default 1/3 (MATLAB’s default). The rescaling below divides by sqrt((1/w0 - 1) / (n + 1)), which -> 0 as w0 -> 1, so values close to the open upper bound blow the non-center points up without limit and degrade the sample covariance accordingly (e.g. w0 = 1 - 1e-12 at n=3 places points near 1.7e6 and leaves the sample covariance off from the identity by ~1.1e-4) – the same kind of near-boundary blow-up student_t_cubature_points() documents for dof near its own open lower bound. The 0 < w0 < 1 check does not, and is not meant to, guard against this.
alpha (float, optional) – Positive spread factor for the non-center points. Default 1.0 (unscaled spherical simplex points). Values of alpha and w0 with alpha**2 < 1 - w0 drive the center weight negative – see Notes.

Returns:

points (ndarray) – Shape (n + 2, n).
weights (ndarray) – Shape (n + 2,), summing to 1. The center weight w0 / alpha**2 + (1 - 1 / alpha**2) can be negative (see Notes); this is inherent to the scaled construction, not an error. Covariances assembled from these points must not use a sqrt-of-weights factorization.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

Relation to pytcl’s other sigma/cubature-point generators. All three trace to the same Julier & Uhlmann unscented-transform lineage but are structurally distinct rules, not reparameterizations of one another:

unscented_transform_points() is the classic symmetric sigma-point set: 2n+1 points (center plus an antipodal pair per axis), exact through degree 3 for N(0, I) because the antipodal pairs cancel every odd-order term.
ckf_spherical_cubature_points() is the CKF’s 2n-point spherical-radial rule (no center point, all weights equal), also exact through degree 3.
This function is the (n+2)-point spherical simplex set: no antipodal symmetry (only the first coordinate axis has a mirrored pair; every later dimension’s simplex vertices are one-sided), so it is exact through degree 2 only – matching the mean and covariance exactly but not third moments in general. E.g. at n=3, w0=1/3, alpha=1, the rule gives E[x2^3] = sqrt(2) against a true value of 0, even though E[x0^3] = 0 by the one axis that does have a mirrored pair.

A user choosing between these should pick unscented_transform_points or ckf_spherical_cubature_points (degree-3 exact, standard choices for the UKF/CKF) unless the point budget must be as small as possible: this rule’s n+2 points is the fewest of the three, at the cost of the third-moment accuracy the other two get for free from symmetry.

Negative center weight. Unlike the unscaled (alpha=1) spherical simplex points, for which 0 < w0 < 1 keeps every weight positive, the scaled center weight w0 / alpha**2 + (1 - 1 / alpha**2) goes negative whenever alpha**2 < 1 - w0 (e.g. w0=1/3, alpha=0.5 gives a center weight of -5/3). This is a real, disclosed property of Julier’s scaled construction – corrected here per Equation 15 of Julier (2002) (MATLAB’s own header notes that Equation 24 of the same paper has a typo the code does not follow), not suppressed or clamped.

randomize (MATLAB’s optional post hoc random-orthonormal-rotation parameter, used in Straka et al. (2012) and Dunik et al. (2011) to avoid repeated-orientation artifacts in tracking) is not exposed; callers who want it can rotate the returned points themselves.

Examples

>>> pts, w = second_order_cubature_points(3)
>>> pts.shape
(5, 3)
>>> round(float(w.sum()), 12)
1.0
>>> np.allclose(np.sum(w[:, None] * pts, axis=0), 0.0, atol=1e-12)  # E[x] = 0
True
>>> resid = pts - np.sum(w[:, None] * pts, axis=0)
>>> np.allclose(resid.T @ (w[:, None] * resid), np.eye(3), atol=1e-12)  # Cov = I
True

References

S. J. Julier, “The scaled unscented transformation,” in Proc. American Control Conference, Anchorage, AK, 8-10 May 2002, pp. 4555-4559.

S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proceedings of the IEEE, vol. 92, no. 3, pp. 401-422, Mar. 2004.

O. Straka, D. Dunik, and M. Simandl, “Randomized unscented Kalman filter in tracking,” in Proc. 15th Int. Conf. on Information Fusion, Singapore, 2012, pp. 503-510.

J. Dunik, O. Straka, and M. Simandl, “The development of a randomised unscented Kalman filter,” in Proc. 18th World Congress, IFAC, Milan, Italy, 2011, pp. 8-13.

pytcl.mathematical_functions.numerical_integration.fifth_order_cubature_points(n)[source]

Degree-5 cubature points for the standard normal N(0, I).

The 2n^2 + 1 point fully-symmetric rule E_n^{r^2} 5-3 of Stroud (1971), the counterpart of the MATLAB TCL’s fifthOrderCubPoints. Exactly integrates every polynomial of total degree <= 5 against N(0, I).

Parameters:

n (int) – Dimension, n >= 1.

Returns:

points (ndarray) – Shape (2*n*n + 1, n).
weights (ndarray) – Shape (2*n*n + 1,), summing to 1. For n > 4 the axis-point weight (4 - n)/(2 (n+2)^2) is negative; this is inherent to the rule, not an error. Covariances assembled from these points must not use a sqrt-of-weights factorization.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = fifth_order_cubature_points(3)
>>> pts.shape
(19, 3)
>>> round(float(w.sum()), 12)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 4)), 12)  # E[x^4] = 3
3.0

pytcl.mathematical_functions.numerical_integration.fourteenth_order_cubature_points(n, beta=0.0)[source]

Degree-14 cubature points for the standard normal N(0, I), n = 3 only.

The counterpart of the MATLAB TCL’s fourteenthOrderCubPoints: lifts the 72-point degree-14 spherical-surface rule (_fourteenth_order_unit_sphere_points_3d(), Stroud’s U3 14-1 (1971)) to N(0, I) times |x|^beta via sphere_surface_to_gauss_points() (the same adapter spherical_radial_points() uses), rather than duplicating the radial weight machinery. Exactly integrates every polynomial of total degree <= 14 against N(0, I).

Unlike fifth_order_cubature_points(), seventh_order_cubature_points(), and spherical_radial_points(), this rule has no documented n-dimensional generalization in the source – MATLAB’s fourteenthOrderCubPoints and fourteenthOrderSpherSurfCubPoints both hardcode if(numDim~=3) error('Only 3D points are supported'); end. So n = 3 here is not a lower bound, it is the only supported value.

Parameters:

n (int) – Dimension; only n = 3 is supported (matches the MATLAB source’s restriction).
beta (float, optional) – Exponent of |x| in the weighting function, beta > -n. Default 0.0 (plain N(0, I)).

Returns:

points (ndarray) – Shape (288, 3) – 72 surface points times 4 radial nodes.
weights (ndarray) – Shape (288,). Sums to 1 when beta=0; otherwise to 2**(beta / 2) * gamma((n + beta) / 2) / gamma(n / 2) (see sphere_surface_to_gauss_points()).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = fourteenth_order_cubature_points(3)
>>> pts.shape
(288, 3)
>>> round(float(w.sum()), 9)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 14)), 3)  # E[x^14] = 135135
135135.0

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Formula U3 14-1, p. 302.

pytcl.mathematical_functions.numerical_integration.genz_keister_points(n, m, algorithm=0, eps_val=None)[source]

Genz-Keister nested cubature points for the standard normal N(0, I).

Counterpart of the MATLAB TCL’s GenzKeisterPoints: a fully-symmetric interpolatory rule (Genz (1986), extended by Genz and Keister (1996)) built from a single 1-D sequence of NESTED generator magnitudes lambda_0=0, lambda_1=sqrt(3), lambda_2, ..., lambda_M – each algorithm’s full table is a strict superset of the previous one’s prefix, which is what makes the resulting n-D point sets nest across increasing m too (see Notes). This is the enabling primitive for Smolyak sparse grids: sparse-grid construction reuses function evaluations across levels only when the levels’ point sets nest, which plain (non-nested) Gauss-Hermite point sets do not do.

Parameters:

n (int) – Dimension, n >= 1.
m (int) – Controls both the number of generator magnitudes available and the rule’s exactness degree – see Notes for exactly how. Must satisfy 1 <= m <= 17 for algorithm=0, 1 <= m <= 15 for algorithm=1 (len(lambda table) - 1 for each; MATLAB’s own documented bounds).
algorithm (int, optional) –
Which tabulated generator to use (MATLAB’s algorithm selector):
- 0 (default): \(Q_P\), built from the increment vector nu = [3, 5, 8] – the first column of Table 3.4 in Genz and Keister (1996).
- 1: \(\hat{Q}_P\), built from nu = [4, 10] – the second column of Table 3.4 in Genz and Keister (1996).
MATLAB also accepts a vector of custom nu increments (via the internal getGenzKeisterGenerators/computeAValues routines) to derive a new generator from scratch. This port does not implement that path – MATLAB’s own docstring calls it “generally a bad idea” and warns of “a severe loss of precision” for the higher-order lambda terms it produces; the two tabulated generators above were computed by Genz and Keister in extended precision specifically to avoid that, and are what every practical use of this rule wants. algorithm values other than 0 or 1 raise ValueError.
eps_val (float, optional) – Points whose |weight| is at or below this threshold are dropped (MATLAB’s pruning of “numerically zero” weights, which is what keeps the point count far below the naive \((2m+1)^n\) bound – see Notes). Default None uses eps(1) in double precision (numpy.finfo(numpy.float64).eps), matching MATLAB’s default.

Returns:

points (ndarray) – Shape (num_points, n).
weights (ndarray) – Shape (num_points,), summing to 1. Genz-Keister rules commonly produce negative weights (visible directly in the worked examples below); these are inherent to the construction and are never suppressed, only the numerically-zero ones are dropped per eps_val.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

What `m` means. Write a “level” as a partition \(p = (p_1, ..., p_n)\) of nonnegative integers with \(p_1 + ... + p_n \le m\); each \(p_i\) selects \(\lambda_{p_i}\) for dimension \(i\). The rule sums, over every such partition and every one of its signed permutations, a shared per-partition weight (MATLAB’s computeW) – so m is a shared budget for how many generator “levels” the n dimensions may jointly spend, not a per-dimension point count. This module verified two independent, useful consequences of that budget structure numerically (both cross-checked against a literal, line-for-line translation of MATLAB’s computeW/innerTerms4Weights state machine before being trusted, and against an unoptimized brute-force reconstruction that evaluates every raw n-tuple in {0,...,m}^n instead of one representative per symmetry orbit):

General guarantee, any n, up to an n-DEPENDENT ceiling well below each algorithm’s maximum m (verified empirically for n = 1..8 – see the boxed number below for why this document refuses to state it unconditionally, and do not assume it extrapolates past n=8 or past this floor without checking): the rule is exact through total polynomial degree \(2m + 1\) (equivalently, since odd-total-degree monomials are already exact at every m via antipodal symmetry, through the largest even degree \(2m\)), and this bound is sharp there (a mixed-exponent monomial of degree \(2m+2\) fails) – e.g. at n=2, algorithm=0, m=4: E[x1^8 x2^2] rule-integrates to 153.38 against a true value of 105, even though every single-axis monomial through degree 16 is still exact there (next point).

This does NOT hold all the way up to each algorithm’s true maximum m – an earlier version of this note claimed the guarantee held for every m except the algorithm’s own maximum (17 for algorithm 0, 15 for algorithm 1), which is false: the breakdown starts before the maximum once n >= 2, and gets worse, not better, as n grows further. A second attempt narrowed this to a floor of m <= 15 (algorithm 0) / m <= 13 (algorithm 1), verified only for n <= 4 – ALSO wrong, because it stopped measuring exactly where a monomial needing five simultaneously nonzero exponents first becomes expressible: no n <= 4 sweep can even construct (6, 6, 6, 6, 6), which is where algorithm 0’s m = 15 actually breaks. Measured largest m for which every degree-\(2m\) monomial (mixed-exponent included – the worst case is usually mixed, not single-axis, and can require more simultaneously nonzero exponents than a low-n sweep can express at all; see the worked examples below) is exact to a relative error at or below the ~1e-12 roundoff-noise floor, by n:

algorithm

n=1

n=2

n=3

n=4

n=5

n=6

n=7

n=8

0 (max 17)

16

15

15

15

14

14

14

14 (b)

1 (max 15)

14

15 (a)

13

13

13

13

13

13 (c)

(a) No violation was found anywhere in algorithm 1’s valid range at n=2 (worst measured relative error 2.1e-13, at m=15, its true maximum) – unlike every other column, this one is not known to break before the algorithm’s own m ceiling, it was simply never observed to in the swept range.

(b) m=14 is confirmed safe through n=8 (worst measured 6.5e-14). m=15 was confirmed broken at n=5, 6, 7 (3.07e-2 each, remarkably flat) but was NOT itself tested at n=8 – listed as 14 because that is the largest confirmed-safe m, not because m=15 was separately confirmed to break at n=8 too (though nothing measured suggests it would recover).

(c) Algorithm 1’s error AT its own floor m=13 is climbing as n grows – 1.7e-13 at n=5 to 1.1e-12 at n=8 – i.e. it is approaching the same ~1e-12 threshold used to call every other cell in this table “safe”. This floor is not comfortably safe indefinitely; re-verify before using it past n=8.

A single conservative number usable without checking n first: the guarantee above is verified for m <= 14 (algorithm 0), m <= 13 (algorithm 1), for every n in 1..8 – and ONLY for n in 1..8. Do not use this number, or any other number in this docstring, for n > 8 without remeasuring: a measurement ceiling silently becoming an implied claim boundary is exactly how the (false) m <= 15 floor above happened, twice, in successive revisions of this note. Above the stated floor, consult the table (or remeasure for your own n) rather than assume the generic bound holds – e.g. at n=3, algorithm=0, m=16 (one step above n=3’s floor of 15), the worst degree-32 monomial – E[x2^32] (0, 32, 0) – is off by a relative 3.2e-2; at n=3, algorithm=1, m=14 (one step above n=3’s floor of 13 there), the worst degree-28 monomial is off by a relative 4.5e-2 – notably higher than the 3.1e-3 a single-axis-only probe (28, 0, 0) would suggest, because the true worst case there, (0, 6, 22), is mixed-exponent. Most strikingly, at n=5, algorithm=0, m=15 (one step above n=5’s floor of 14), the worst monomial is (6, 6, 6, 6, 6) – true value 759375, rule value 736043.78, a relative 3.07e-2 – while the single-axis probe (30, 0, 0, 0, 0) at that SAME rule reads a clean 3.8e-12: five simultaneously nonzero exponents were required to see the failure at all. See “Precision at the top of the range” below for why this degrades smoothly rather than cutting off, and why it starts sooner as n grows.

Bonus at specific “milestone” m values, single-axis moments only: MATLAB’s tabulated \(\lambda\)/\(a\) values were tuned by Genz and Keister so that at the m where each nu stage completes, single-axis moments (equivalently, the n=1 rule itself) are exact to a degree well beyond \(2m+1\) – this is what “the first column of Table 3.4” in Genz and Keister (1996) actually tabulates (this port does not have that table’s literal text, only the MATLAB source’s transcribed \(\lambda\)/\(a\) constants, so the milestone degrees below were independently determined by direct computation against the closed-form N(0,1) moments, not copied from the paper):

algorithm	m	points (n=1)	degree	nu stage
0 (\(Q_P\))	1	3	5	base (0, +-sqrt(3))
0	3	7	7	(no bonus – not a stage boundary)
0	4	9	15	nu[0]=3 complete
0	8	17	17	(no bonus)
0	9	19	29	nu[1]=5 complete
0	15	31	31	(no bonus)
0	16	33	33	(no bonus)
0	17	33	see below	nu[2]=8 “complete” – precision breakdown
1 (\(\hat{Q}_P\))	1	3	5	base
1	3	7	7	(no bonus)
1	4	9	9	(no bonus)
1	5	11	19	nu[0]=4 complete
1	10	21	21	(no bonus)
1	14	29	29	(no bonus)
1	15	29	see below	nu[1]=10 “complete” – precision breakdown

(m values between milestones, and m=2/m=1 or m=14/m=… duplicates, reuse the previous milestone’s point set post-pruning – see the nesting property below.) Mixed-exponent moments do not inherit this bonus at any m, milestone or not (point 1 above).

Precision at the top of the range (disclosed, not a porting bug). The rows marked “see below” in the milestone table above are real and reproducible, not an artifact of this port, but they are not a “degree” claim at all: at m = 17 (algorithm 0) and m = 15 (algorithm 1) – the last stage boundary, which is also the maximum m MATLAB’s own docstring allows – there is no accuracy cliff to report a single degree for, for the n=1 marginal those milestone rows describe. Every other tested m has relative error staying flat (~1e-16) up to some degree and then cutting off sharply; at the maximum m of each algorithm the relative error instead grows smoothly with degree from the start, so “the degree this rule is exact to” becomes purely a function of whatever tolerance is used to define “exact”. Measured directly for algorithm 0 at m=17 (n=1 marginal, single-axis moments): relative error is 2.2e-16 at degree 0 and grows to 5.6e-7 by degree 24 with no step anywhere in between; the resulting degree at a chosen relative tolerance is 35 @ 1e-4, 25 @ 1e-6, 19 @ 1e-8, 13 @ 1e-10, 9 @ 1e-12, 5 @ 1e-14 – no tolerance reproduces a “15” here; there is no such number, at any tolerance.

The SAME smooth-growth-not-cliff pattern, with the SAME underlying cause, is what makes item 1’s n-dependent floor above sit below each algorithm’s true maximum m once n >= 2: e.g. at n=3, algorithm=0, m=16 (n=3’s floor there is 15, one below), the single-axis marginal E[x2^d] measured at even d from 0 to 32 grows from 4.9e-15 to 3.2e-2 with no step in between either – it is the identical phenomenon as the n=1/m=17 case above, just triggered one m earlier. This traces to the published double-precision \(\lambda\)/\(a\) constants themselves: achieving the top stage’s designed accuracy requires near-exact cancellation among terms spanning roughly 15 to 16 orders of magnitude (a ranges up to 2.92e15 at algorithm 0’s last entry), which is at the edge of what float64 can represent faithfully even with exact-precision inputs; at n=1 only that single generator sequence is ever combined against itself, but at higher n the shared per-partition weight table (computeW) combines it across combinatorially more partitions, so the same precision loss surfaces at a lower m (this explains the direction of the n-dependence empirically measured in item 1, not a claim independently re-derived from the constants). Do not rely on exactness AT ALL above the n-dependent floor established in item 1 above (m <= 14 for algorithm 0, m <= 13 for algorithm 1, and ONLY as verified there for n <= 8) – not even up to the generic \(2m+1\) bound that holds at or below that floor, within that same n range (see the counterexamples there: a 3.2e-2 relative error at n=3, algorithm=0, m=16, degree 32; a 4.5e-2 relative error at n=3, algorithm=1, m=14, degree 28; and, requiring a mixed monomial no n<=4 sweep could express, a 3.07e-2 relative error at n=5, algorithm=0, m=15, degree 30, at (6,6,6,6,6), versus a clean 3.8e-12 at the single-axis probe (30,0,0,0,0) on that same rule).

Nesting. For every consecutive pair m-1, m except the last one (m = 17 for algorithm 0, m = 15 for algorithm 1), the point set at m contains every point of m - 1 to floating tolerance – verified explicitly for n = 1, 2, 3 across the full valid range. At the excluded top pair, the newly available generator magnitude does not simply add to the previous point set: recomputing the shared weight table with the larger m shifts an existing extremal point’s weight across the eps_val pruning threshold, dropping it as a new, unrelated extremal point becomes available. This is the same precision effect described above, viewed through its effect on the point set rather than on exactness degree, and follows directly from it (both are consequences of the last stage’s constants being computed by Genz and Keister in extended precision, but published here, and by MATLAB, at double precision).

Point count. MATLAB’s docstring describes the point-count formula “given in the text after Equation 1” of Genz and Keister (1996); this port does not preallocate (MATLAB does, to size xi/w before filling them) and instead grows the point list directly, so that formula was not needed here.

randomize (MATLAB’s optional post hoc random-orthonormal-rotation parameter, see Straka et al. (2012), Dunik et al. (2011)) is not exposed; callers who want it can rotate the returned points themselves.

Examples

>>> pts, w = genz_keister_points(1, 1)
>>> sorted(pts.ravel().round(10).tolist())
[-1.7320508076, 0.0, 1.7320508076]
>>> round(float(w.sum()), 12)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 4)), 9)  # E[x^4] = 3, exact at m=1
3.0

>>> pts9, w9 = genz_keister_points(1, 4)  # milestone m=4: bonus degree 15
>>> pts9.shape
(9, 1)
>>> round(float(np.sum(w9 * pts9[:, 0] ** 14)), 3)  # E[x^14] = 135135
135135.0

>>> _, w19 = genz_keister_points(1, 9)  # milestone m=9: bonus degree 29
>>> bool((w19 < 0).any())  # Genz-Keister rules commonly have negative weights
True

>>> pts1, w1 = genz_keister_points(1, 1)
>>> pts3, w3 = genz_keister_points(1, 3)  # nesting: m=1's points subset of m=3's
>>> all(np.any(np.isclose(pts3, p, atol=1e-12)) for p in pts1)
True

References

A. Genz and B. D. Keister, “Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight,” Journal of Computational and Applied Mathematics, vol. 71, no. 2, pp. 299-309, Jul. 1996.

A. Genz, “Fully symmetric interpolatory rules for multiple integrals,” SIAM Journal on Numerical Analysis, vol. 23, no. 6, pp. 1273-1283, Dec. 1986.

O. Straka, D. Dunik, and M. Simandl, “Randomized unscented Kalman filter in tracking,” in Proc. 15th Int. Conf. on Information Fusion, Singapore, 2012, pp. 503-510.

J. Dunik, O. Straka, and M. Simandl, “The development of a randomised unscented Kalman filter,” in Proc. 18th World Congress, IFAC, Milan, Italy, 2011, pp. 8-13.

pytcl.mathematical_functions.numerical_integration.seventh_order_cubature_points(n, algorithm=None)[source]

Degree-7 cubature points for the standard normal N(0, I).

The counterpart of the MATLAB TCL’s seventhOrderCubPoints, full algorithm surface. algorithm=None reproduces MATLAB’s default selection logic: n == 1 -> algorithm 9, n == 2 -> algorithm 2, otherwise -> algorithm 0. Each algorithm exactly integrates every polynomial of total degree <= 7 against N(0, I) for its own valid n; that has been verified by the test suite exactly for the (algorithm, n) pairs in the table below – no wider claim is made.

algorithm rule (Stroud 1971) valid n points ——— ————————— ————– —————- 0 E_n^{r^2} 7-3, p. 319 n = 3..6 (note) 2*(2^n + 2n^2) 1 E_n^{r^2} 7-1, p. 318 3, 4, 6, 7 2^n + 2n^2 + 1 2 E_2^{r^2} 7-1, p. 324 2 12 3 E_2^{r^2} 7-2, p. 324 2 17 (see note) 4 / 5 E_3^{r^2} 7-1, p. 327 3 27 6 / 7 E_3^{r^2} 7-2, p. 328 3 33 8 E_4^{r^2} 7-1, p. 329 4 49 (see note) 9 quadraturePoints1D(4) 1 4

Note on algorithm 0: it predates this table (it is the pre-existing default for n > 2) and its code, unlike every other algorithm here, does not restrict n beyond n >= 3 – it will run for any such n. But the test suite exercises it only at n = 3..6, so per this docstring’s own rule (“verified… exactly for the (algorithm, n) pairs in the table below – no wider claim is made”) the exactness CLAIM is bounded to n = 3..6, matching every other row; n > 6 is accepted by the code but is an unverified extrapolation, not a documented guarantee.

Deviations from MATLAB (algorithms 3 and 8). MATLAB’s comments document scale corrections for these two algorithms (a missing 4/3 factor for algorithm 3’s r and s; a sqrt(4/5) factor for algorithm 8’s r, s, and t). Both were checked here with exact symbolic arithmetic against every monomial of total degree <= 7 and neither actually achieves degree-7 exactness – MATLAB’s own documented fixes are themselves incorrect (or fix a different bug than the one that matters), independent of what the underlying formula’s own literal coefficients are. Algorithm 8’s real defect is a plain transcription typo (the book’s t = 3 + sqrt(3) is missing an outer square root; with t = sqrt(3 + sqrt(3)) and Stroud’s other coefficients unscaled, no sqrt(4/5) correction is needed and the rule is exact). Algorithm 3’s 16-point, no-origin layout is provably incapable of degree-7 exactness for any choice of its two radii and three weights (5 free parameters short by one against 6 independent moment constraints, confirmed by direct solve and by least-squares residual minimization); adding a 17th point at the origin supplies the missing degree of freedom and, with Stroud’s original (unscaled) r and s, yields an exact rule. See each private constructor’s docstring (_e2_7_2(), _e4_7_1()) for the full derivation. Both corrected rules still integrate exactly through degree 7 and fail at degree 8 (confirmed sharp); neither reproduces MATLAB’s seventhOrderCubPoints numeric output for that algorithm.

Parameters:

n (int) – Dimension. Valid values depend on algorithm – see the table above; an unsupported (algorithm, n) pair raises ValueError.
algorithm (int, optional) – Which of the 10 algorithms (0-9) above to use. Default None reproduces MATLAB’s default selection (see above).

Returns:

points (ndarray) – Shape (num_points, n); num_points per the table above.
weights (ndarray) – Shape (num_points,), summing to 1. For algorithm 0, the axis shell’s surface weight (8 - n)/(n(n+2)(n+4)) is negative for n > 8; this is inherent to the rule, not an error. Algorithm 3’s corrected 17-point rule (see the deviation note above) likewise has a negative weight, D = -2/3, on its origin point – also inherent, not an error. In both cases, covariances assembled from these points must not use a sqrt-of-weights factorization.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = seventh_order_cubature_points(3)
>>> pts.shape
(52, 3)
>>> round(float(w.sum()), 12)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 6)), 6)  # E[x^6] = 15
15.0

>>> pts, w = seventh_order_cubature_points(2, algorithm=3)
>>> pts.shape
(17, 2)
>>> round(float(w.sum()), 12)
1.0

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Formulas E_n^{r^2} 7-1 through 7-3 and E_2^{r^2}/E_3^{r^2}/E_4^{r^2} 7-1/7-2, pp. 318-329. Algorithm 0’s formula as printed there contains a typo; the corrected form used here follows Stroud’s original papers, Stroud (1967) and Stroud (1968), summarized in Crouse (2014).

A. H. Stroud, “Some seventh degree integration formulas for symmetric regions,” SIAM Journal on Numerical Analysis, vol. 4, no. 1, pp. 37-44, Mar. 1967.

A. H. Stroud, “Some seventh degree integration formulas for the surface of an n-sphere,” Numerische Mathematik, vol. 11, no. 3, pp. 273-276, Mar. 1968.

D. F. Crouse, “Basic tracking using nonlinear 3D monostatic and bistatic measurements,” IEEE Aerospace and Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.

pytcl.mathematical_functions.numerical_integration.smolyak_points(n, level, algorithm=0)[source]

Smolyak sparse-grid cubature points over nested Genz-Keister sequences.

ORIGINAL DESIGN: unlike the rest of this module, this function has no MATLAB TCL counterpart – the MATLAB library provides the nested Genz-Keister 1-D sequences (genz_keister_points) but never the Smolyak combination over them. It is the standard sparse-grid construction from the literature (Smolyak (1963); Genz and Keister (1996) build their sequences expressly for it), implemented and verified here from scratch; every exactness claim below is bounded to the measured grid it was verified on.

The classic combination formula: for multi-indices \(k = (k_1, ..., k_n)\) with \(L - n + 1 \le |k| \le L\) (\(L\) = level, components >= 0), sum the tensor products of the 1-D rules at levels \(k_i\) with coefficient \((-1)^{L - |k|} \binom{n-1}{L - |k|}\). The 1-D rule at level q is the Genz-Keister rule at the m where the algorithm’s q-th nu stage completes (the “milestone” m values of genz_keister_points’ bonus table; see _SMOLYAK_GK_M for the mapping and the full rationale):

algorithm	level	GK m	1-D points	1-D degree
0	0	0	1	1
0	1	1	3	5
0	2	4	9	15
0	3	9	19	29
1	0	0	1	1
1	1	1	3	5
1	2	5	11	19

Each algorithm’s ladder stops BELOW its table’s top milestone (m=17 / m=15): there the published double-precision constants produce no exactness cliff at any tolerance (see “Precision at the top of the range” in genz_keister_points()’ docstring), so no Smolyak level is built on them – hence the level caps of 3 (algorithm 0) and 2 (algorithm 1).

Because the 1-D point sets nest exactly (each level’s nodes are a strict superset of the previous level’s – the same table constants, so exact float equality), points repeated across the combination’s tensor grids are merged into single points with summed weights, which is what keeps the point count far below the tensor product’s.

Parameters:

n (int) – Dimension, n >= 1.
level (int) – Smolyak accuracy level, 0 <= level <= 3 for algorithm=0, 0 <= level <= 2 for algorithm=1 (the milestone ladders above). level=0 is the single point at the origin.
algorithm (int, optional) – Which tabulated Genz-Keister generator to build on – 0 (default, \(Q_P\), nu=[3,5,8]) or 1 (\(\hat{Q}_P\), nu=[4,10]); see genz_keister_points().

Returns:

points (ndarray) – Shape (num_points, n), rows sorted lexicographically. The grids nest across levels: every point of level appears (exactly, as floats) in level + 1’s grid.
weights (ndarray) – Shape (num_points,), summing to 1. Commonly contains negative values – already at n=4, level=1 the origin’s weight is exactly -3 + 4*(2/3) = -1/3. This is a real, disclosed property of the Smolyak combination (its coefficients alternate in sign) compounded by the Genz-Keister rules’ own negative weights, not suppressed or clamped. Covariances assembled from these points must not use a sqrt-of-weights factorization.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

Measured exactness (the claim, and its exact boundary). The standard result for nested 1-D sequences whose level-q member is exact through degree >= 2q+1 gives total-degree exactness \(2 \cdot \mathrm{level} + 1\); the GK milestone degrees (1, 5, 15, 29) dominate 2q+1 at every rung, so that floor applies here. Because those milestone degrees grow much faster than 2q+1, low dimensions do better. Rather than assert the generic floor and hope, the actual total-degree exactness was measured per (algorithm, n, level) cell – every all-even-exponent monomial scanned per total degree at relative tolerance 1e-12 against the closed-form N(0, I) moments, odd exponents exact by antipodal symmetry (re-verified at looser tolerance in the test suite’s through-degree sweep):

algorithm 0	level 0	level 1	level 2	level 3
n=1	1	5	15	29
n=2	1	3	7	11
n=3	1	3	5	9
n=4 .. 8	1	3	5	7

algorithm 1	level 0	level 1	level 2
n=1	1	5	19
n=2	1	3	7
n=3 .. 6	1	3	5

Every cell meets or exceeds 2*level+1, and every cell is sharp: the next even degree fails, with a measured relative error between 2.1e-4 (algorithm 0, n=1, level 3, on E[x^30]) and 1.0 (the (2,2,…) mixed monomials), against worst in-range noise of 2.7e-14. These tables are claims ONLY for the (algorithm, n, level) cells they list – measured for n <= 8 (algorithm 0) and n <= 6 (algorithm 1). For larger n the generic 2*level+1 floor is the standard theoretical result but was NOT measured here; re-verify before relying on it (see tests/unit/test_cubature_points.py::TestSmolyakPoints for the measurement reproduced as tests).

Point count vs. the tensor grid. At level 2 (degree >= 5) in n = 8, this grid has 177 points; the tensor product of the 3-point degree-5 Gauss-Hermite rule has 3^8 = 6561. At level 3 (degree >= 7), 1377 points against the 4-point rule’s 4^8 = 65536. The sparse count grows polynomially in n at fixed level, the tensor count exponentially.

Examples

>>> pts, w = smolyak_points(2, 1)
>>> pts.shape  # cross of the 3-point GK rule, origin merged: 5, not 9
(5, 2)
>>> pts.round(10).tolist()[:2]
[[-1.7320508076, 0.0], [0.0, -1.7320508076]]
>>> round(float(w.sum()), 12)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 2)), 12)  # E[x^2] = 1
1.0

>>> pts, w = smolyak_points(8, 2)  # degree 5 in n=8: 177 points, not 3^8
>>> pts.shape
(177, 8)
>>> bool((w < 0).any())  # negative weights are inherent, see above
True

References

S. A. Smolyak, “Quadrature and interpolation formulas for tensor products of certain classes of functions,” Doklady Akademii Nauk SSSR, vol. 148, no. 5, pp. 1042-1045, 1963.

A. Genz and B. D. Keister, “Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight,” Journal of Computational and Applied Mathematics, vol. 71, no. 2, pp. 299-309, Jul. 1996.

F. Heiss and V. Winschel, “Likelihood approximation by numerical integration on sparse grids,” Journal of Econometrics, vol. 144, no. 1, pp. 62-80, May 2008.

pytcl.mathematical_functions.numerical_integration.sphere_surface_to_gauss_points(surface_points, surface_weights, degree, beta=0.0)[source]

Lift a spherical-surface cubature rule to a Gaussian(-times-|x|^beta) rule.

Counterpart of the MATLAB TCL’s spherSurfPoints2GaussPoints: given cubature points/weights for the uniform measure on the unit sphere S^(n-1) (weights summing to 1, e.g. from _sphere_surface_points()), produces points/weights for the weighting function w(x) = N(x; 0, I) * |x|^beta (beta = 0 is the plain N(0, I) density).

The MATLAB source builds the radial rule from a 1-D quadrature for |x|^c1 * exp(-x^2) (quadraturePoints1D, algorithm 9, a three-term-recursion construction restricted to integer c1), then rescales by x -> sqrt(2) x. This port instead reuses the generalized Gauss-Laguerre substitution t = r^2/2 already used by spherical_radial_points(): with alpha = (n + beta)/2 - 1, the quadrature nodes/weights from scipy.special.roots_genlaguerre match the same target radial moments integral_0^inf r^(n-1+beta) exp(-r^2/2) r^(2k) dr for every k needed up to degree, so it reproduces the identical family of rules while additionally allowing non-integer beta (MATLAB’s three-term-recursion route cannot). Both routes implement the same “spherical shell plus a 1-D |x|^beta * exp(-x^2)-type formula” construction described in Chapter 2.8 of Stroud (1971), cited by the MATLAB source. Randomization (MATLAB’s randomize flag, a random orthonormal rotation applied post hoc to reduce repeated-orientation artifacts in tracking – see Straka et al. (2012), Dunik et al. (2011)) is not exposed; callers who want it can rotate the returned points themselves.

Parameters:

surface_points (array_like) – Points on the unit sphere S^(n-1), shape (num_surface_points, n). Precondition (not validated): every row has unit norm – a rule built from points off the sphere silently integrates a different weighting than the one this function documents.
surface_weights (array_like) – Weights for the uniform measure on S^(n-1), shape (num_surface_points,), summing to 1 (checked: must be 1-D and its length must match surface_points’s row count, or ValueError is raised). Precondition (not validated): the weights actually sum to 1.
degree (int) – The polynomial degree the surface rule (and thus this rule) is exact through, an integer >= 1. Unlike spherical_radial_points(), EVEN values are accepted here – this is the lower-level adapter fourteenth_order_cubature_points() itself calls with degree=14 – so only integer-ness is checked, not oddness. Precondition (not validated): this must match the degree the supplied surface_points/surface_weights rule was actually built for – an inconsistent value is accepted silently and just mislabels the resulting rule’s true degree.
beta (float, optional) – Exponent of |x| in the weighting function, beta > -n. Default 0.0 (plain N(0, I)).

Returns:

points (ndarray, shape (num_radii * num_surface_points, n))
weights (ndarray, shape (num_radii * num_surface_points,). Sums to) – 2**(beta / 2) * gamma((n + beta) / 2) / gamma(n / 2), the beta-th absolute moment of the chi_n distribution – 1.0 when beta=0, not 1 in general.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> surf_pts, surf_w = _sphere_surface_points(3, 5)
>>> pts, w = sphere_surface_to_gauss_points(surf_pts, surf_w, 5)
>>> round(float(w.sum()), 12)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 4)), 9)  # E[x^4] = 3
3.0

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Ch. 2.8.

O. Straka, D. Dunik, and M. Simandl, “Randomized unscented Kalman filter in tracking,” in Proc. 15th Int. Conf. on Information Fusion, Singapore, 2012, pp. 503-510.

J. Dunik, O. Straka, and M. Simandl, “The development of a randomised unscented Kalman filter,” in Proc. 18th World Congress, IFAC, Milan, Italy, 2011, pp. 8-13.

pytcl.mathematical_functions.numerical_integration.spherical_radial_points(n, degree, beta=0.0)[source]

Arbitrary-odd-degree spherical-radial cubature points for N(0, I) times |x|^beta.

Product of a generalized Gauss-Laguerre radial rule (exact for all required even powers of r) with a dimension-recursive surface rule on the unit sphere. Generalizes the 3rd-degree spherical-radial rule of the CKF to any odd degree. beta selects the sphere_surface_to_gauss_points() weight family (|x|^beta times the Gaussian); the default beta=0.0 is the plain N(0, I) case and its code path is byte-for-byte the original (pre-beta-parameter) implementation, so existing callers are unaffected.

The point count grows roughly as (degree/2)^(n-1) from the surface rule; for the common degrees 5 and 7 prefer fifth_order_cubature_points() and seventh_order_cubature_points(), which grow polynomially in n.

Parameters:

n (int) – Dimension, n >= 1.
degree (int) – Odd polynomial degree >= 3 the rule integrates exactly.
beta (float, optional) – Exponent of |x| in the weighting function, beta > -n. Default 0.0 (plain N(0, I); weights then sum to 1).

Returns:

points (ndarray, shape (num_points, n))
weights (ndarray, shape (num_points,). Summing to 1 when beta=0;) – otherwise to 2**(beta / 2) * gamma((n + beta) / 2) / gamma(n / 2) (see sphere_surface_to_gauss_points()).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

Deviation from MATLAB at n=1. MATLAB’s arbOrderGaussCubPoints (the beta != 0 routine this function’s beta path delegates to, via sphere_surface_to_gauss_points()) hard-errors at numDim==1 (error('numDim must be >1.')); this port accepts n=1 for both the beta=0 and beta != 0 paths and returns a correct rule instead – S^0 = {-1, +1} is a valid, if degenerate, sphere, and the construction handles it without special-casing. This was already true of the pre-existing beta=0 path; the beta path added alongside this note now exercises the same n=1 case (via sphere_surface_to_gauss_points(), which never had MATLAB’s numDim>1 restriction to begin with), so it inherits the same disclosed divergence.

Examples

>>> pts, w = spherical_radial_points(2, 5)
>>> round(float(w.sum()), 12)
1.0
>>> round(float(np.sum(w * pts[:, 0] ** 4)), 10)  # E[x^4] = 3
3.0

pytcl.mathematical_functions.numerical_integration.student_t_cubature_points(n, dof)[source]

Third-order cubature points for the standard multivariate Student-t.

Counterpart of the MATLAB TCL’s thirdOrderStudentTCubPoints: unit points for the multivariate Student-t distribution with zero mean, identity scale matrix, and dof degrees of freedom – the Student-t analogue of ckf_spherical_cubature_points()’s 2n-point Gaussian rule. To use Student-t process or measurement noise in a cubature filter, swap these points/weights in wherever the Gaussian cubature points would otherwise go (e.g. via ckf_predict’s points/weights arguments), then affinely map by the target mean and scale-matrix square root exactly as the Gaussian points are (see transform_cubature_points()).

MATLAB’s mu (mean) and SR (lower-triangular scale-matrix square root) arguments are folded into that same affine-map step used everywhere else in this module rather than taken here directly: with SR = I MATLAB’s xi(:,curPoint) = nuConst * SR(:,curDim) reduces to nuConst times the standard basis vectors, which is exactly what this function returns before the caller (or transform_cubature_points) maps them to a specific mean/scale.

Parameters:

n (int) – Dimension, n >= 1.
dof (float) – Degrees of freedom, dof > 2. The constraint comes directly from nuConst = sqrt(dof * n / (dof - 2)) in the source: below dof = 2 the multivariate Student-t’s covariance dof / (dof - 2) * Sigma itself does not exist (the term the rule is built to reproduce, see Notes), so there is no covariance left for a third-order rule to match.

Returns:

points (ndarray) – Shape (2n, n): pairs (+nuConst * e_i, -nuConst * e_i) for each axis i, interleaved in that order (MATLAB’s curPoint loop – not grouped into all-positive-then-all-negative blocks the way ckf_spherical_cubature_points() is).
weights (ndarray) – Shape (2n,), each 1 / (2n), summing to 1.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

Why third order, not higher. The rule matches the distribution’s mean (0, by antipodal symmetry – every odd-total-degree monomial vanishes for the same reason it does in ckf_spherical_cubature_points()) and its per-axis second moment, E[x_i^2] = dof / (dof - 2), via nuConst**2 = dof * n / (dof - 2) (see the worked example below). It is not exact at degree 4: e.g. the rule’s E[x_i^4] reduces to nuConst**4 / n, which generally disagrees with the Student-t distribution’s true fourth moment 3 * dof**2 / ((dof-2)*(dof-4)) (itself only defined for dof > 4) – see tests/unit/test_cubature_points.py::TestStudentT for the verified numeric gap. As dof -> inf the Student-t distribution converges to N(0, I) and nuConst -> sqrt(n), recovering ckf_spherical_cubature_points exactly.

Examples

>>> pts, w = student_t_cubature_points(3, 6.0)
>>> pts.shape
(6, 3)
>>> round(float(w.sum()), 12)
1.0
>>> nu_const_sq = 6.0 * 3 / (6.0 - 2)
>>> round(float(np.sum(w * pts[:, 0] ** 2)), 10) == round(nu_const_sq / 3, 10)
True

References

Y. Huang, Y. Zhang, N. Li, S. M. Naqvi, and J. Chambers, “A robust Student’s t based cubature filter,” in Proc. 19th Int. Conf. on Information Fusion, Heidelberg, Germany, 5-8 Jul. 2016.

pytcl.mathematical_functions.numerical_integration.transform_cubature_points(points, weights, mean, sqrt_cov)[source]

Affinely map unit cubature points to a given mean and covariance.

Parameters:

points (array_like) – Unit points for N(0, I), shape (num_points, n).
weights (array_like) – Weights, shape (num_points,).
mean (array_like) – Target mean, shape (n,).
sqrt_cov (array_like) – Square root of the target covariance (lower-triangular Cholesky factor S with S @ S.T = P), shape (n, n).

Returns:

points (ndarray) – Transformed points mean + points @ sqrt_cov.T.
weights (ndarray) – Unchanged weights (copied).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> unit = np.array([[1.0], [-1.0]])
>>> w = np.array([0.5, 0.5])
>>> pts, wts = transform_cubature_points(unit, w, [10.0], [[2.0]])
>>> pts.ravel().tolist()
[12.0, 8.0]

pytcl.mathematical_functions.numerical_integration.cubature_point_moments(points, weights, func, mean=None, cov=None)[source]

Propagate a distribution’s first two moments through func by cubature integration.

Counterpart of the MATLAB TCL’s calcCubPointMoments: given a set of cubature points/weights, transform them through func and return the resulting mean and covariance. This is the public, filter-independent form of the transform-then-propagate pattern ckf_predict/ ckf_update (pytcl.dynamic_estimation.kalman.unscented) apply internally to a specific dynamics or measurement function.

MATLAB’s z/S (mean and lower-triangular covariance square root) are mandatory positional arguments there – calcCubPointMoments always affinely maps its xi unit points by them before applying h. This port makes that step optional via mean/cov: when both are given, points are treated as unit points for N(0, I) (or the corresponding unit rule for whatever distribution points/ weights target) and are mapped through transform_cubature_points() first, exactly as MATLAB’s transformCubPoints does; when both are omitted, points are passed to func unchanged, e.g. because they are already points of the target distribution (as when reusing a filter’s own cubature points, already scaled by its state covariance).

MATLAB’s optional innovTrans (custom difference function, e.g. for circular quantities) and meanFun (custom weighted-average function, e.g. for angular means) are not exposed: func’s output is always averaged and differenced with plain arithmetic here. Callers needing a non-Euclidean mean/innovation should wrap func accordingly or average/difference the returned points by hand.

Parameters:

points (array_like) – Cubature points, shape (num_points, n). Unit points for N(0, I) if mean/cov are given; already-in-distribution points otherwise.
weights (array_like) – Weights matching points, shape (num_points,), normally summing to 1. May contain negative values (e.g. higher-degree rules); never suppressed.
func (callable) – Maps a single point, shape (n,), to a transformed point, shape (m,).
mean (array_like, optional) – Target mean, shape (n,). Must be given together with cov.
cov (array_like, optional) – Target covariance, shape (n, n). Must be given together with mean. Need only be positive semi-definite – factored with chol_semi_def(), matching MATLAB’s own docstring recommendation (cholSemiDef(R,'lower')) and the same fallback ckf_predict/ckf_update already use for a filter’s own (possibly rank-deficient) state covariance – rather than a raw Cholesky that would raise on a singular or near-singular cov.

Returns:

mean_out (ndarray) – Weighted mean of func applied to the (possibly transformed) points, shape (m,).
cov_out (ndarray) – Weighted covariance of func applied to the (possibly transformed) points, shape (m, m).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = second_order_cubature_points(2)
>>> mean = np.array([1.0, -2.0])
>>> cov = np.diag([0.5, 2.0])
>>> A = np.array([[2.0, 0.0], [1.0, 1.0]])
>>> b = np.array([0.5, 0.5])
>>> mu, P = cubature_point_moments(pts, w, lambda x: A @ x + b, mean, cov)
>>> np.allclose(mu, A @ mean + b, atol=1e-10)
True
>>> np.allclose(P, A @ cov @ A.T, atol=1e-10)
True

See also

transform_cubature_points: The affine unit-point mapping applied when mean/cov are given.

pytcl.mathematical_functions.numerical_integration.gaussian_lcd_samples(n, num_points, *, force_cov_match=True, rng=None, max_iter=1000)[source]

Localized cumulative distribution (LCD) cubature points for N(0, I).

Port of the MATLAB Tracker Component Library’s GaussianLCDSamples.m (commit 593ce51) entry point, per the Gaussian-LCD port-feasibility design spec (local-only, untracked): the CvM objective/gradient are a faithful transcription (this module’s private functions); the optimizer is scipy.optimize.minimize(method= "L-BFGS-B", jac=True) in place of MATLAB’s MEX-only liblbfgs call, since liblbfgs is a vendored third-party C library with no MATLAB source of its own to port fidelity against (spec Section 6).

Points are generated as 2*(num_points // 2) symmetric pairs +-s_i obtained by minimizing the modified Cramer-von Mises distance between the points’ localized CDF and a standard normal’s, plus one fixed point at the origin when num_points is odd (MATLAB lines 212-218). All weights are uniform, 1 / num_points (MATLAB line 222: w=1/numSamples), summing to 1 – the Gaussian-weight convention this module’s sibling cubature_points functions use (region-measure weighting, where weights sum to a region’s volume rather than 1, is a separate contract used only by region_cubature.py; see that module’s docstring).

Rotation-invariance caveat (read before comparing outputs across runs or against MATLAB). For n >= 2 the CvM cost is provably invariant under any global orthogonal transform applied to every point simultaneously (design spec Section 3): O(n) is a continuous symmetry group of the objective, so a minimizer sits on a flat manifold of equally-optimal solutions, not an isolated point. Two calls that converge correctly – including one call of this function versus a MATLAB GaussianLCDSamples run given the same starting matrix – generically land on different points of that manifold: same CvM cost, different raw coordinates, and this is expected, not a bug. Only n == 1 has a discrete symmetry group (O(1) = {+1, -1}) where raw coordinates (up to sign/permutation) are a meaningful comparison.

Seeding is NOT cross-compatible with MATLAB. rng (default: a fresh numpy.random.default_rng(), PCG64-backed per this repo’s convention) mirrors MATLAB’s sInit=randn(sDims) initialization scheme – drawing a (n, num_points // 2) matrix of standard normal variates via rng.standard_normal – but NumPy’s Generator and MATLAB’s randn use different underlying bit generators and uniform-to-Gaussian transforms. No seed value reproduces the same sInit matrix in both ecosystems; giving this function and a MATLAB call “the same seed” only means “the same kind of random start,” not “the same numbers.” A fixed rng does give bit-identical output across repeated Python calls (checked in the test suite).

L-BFGS-B option mapping (honest, not exact, parity with liblbfgs). MATLAB’s quasiNewtonLBFGS defaults (all left at default by GaussianLCDSamples.m) come from a MEX wrapper around Naoaki Okazaki’s liblbfgs (a C port of Nocedal’s L-BFGS with a More-Thuente line search). scipy’s L-BFGS-B uses a different Fortran implementation (Zhu/Byrd/Lu/Nocedal) with its own line search, so no option mapping can reproduce liblbfgs’s exact step sequence; the mapping below matches intent, not mechanics:

numCorr=6 (history size) -> maxcor=6 (scipy’s history-size option shares the same meaning and default value).
max_iterations=1000 -> maxiter=max_iter (this function’s parameter, default 1000, matching MATLAB’s default exactly) and max_linesearch=20 -> maxfun=max_iter*20. liblbfgs’s max_iterations bounds only outer iterations; scipy’s L-BFGS-B additionally caps total function evaluations via maxfun (default 15000, independent of maxiter), which would silently cut optimization short before liblbfgs’s own worst-case evaluation budget (max_iterations * max_linesearch = 20000 at the defaults) is reached, so maxfun is raised to match that budget rather than left at scipy’s unrelated default.
epsilon=1e-6 (stop when the Euclidean gradient norm drops below epsilon * max(1, ||x||)) -> gtol=1e-6 (the minimize option name; the underlying Fortran variable is called pgtol, which is also the keyword the legacy, non-minimize scipy.optimize.fmin_l_bfgs_b entry point uses – minimize renames it to gtol). This is a genuinely different criterion, not just a differently-named equivalent: scipy’s gtol stops on the max-absolute-component of the (here unconstrained, so unprojected) gradient, with no scaling by the parameter norm. Both use the same numeric threshold as a reasonable value-level match; they are not the same stopping rule.
delta=0 / past=0 (liblbfgs’s relative-function-decrease test explicitly disabled) -> ftol=0.0. This one is a faithful mapping: scipy computes factr = ftol / eps internally and the underlying Fortran L-BFGS-B code treats factr=0 as “suppress this termination test” (its own documented convention), mirroring liblbfgs’s disabled delta-test exactly.
wolfe=0.9 (curvature/Wolfe condition in the More-Thuente line search) and ftol=1e-6 / xtol=1e-16 / min_step / max_step (liblbfgs’s own internal line-search parameters) have no exposed equivalent in scipy’s L-BFGS-B minimize interface – its internal line search (dcsrch) hardcodes its own strong-Wolfe curvature parameter (conventionally 0.9, coincidentally the same value liblbfgs defaults to, but not user-settable through this wrapper) and is not configurable from Python. Not mapped; noted here so a future reader does not assume silent parity.

Parameters:

n (int) – Dimensionality of the cubature points, n >= 1.
num_points (int) – Total number of points to generate, num_points >= 2. For a non-singular sample covariance when force_cov_match is True, MATLAB’s own documented requirement is num_points >= 2*n (see Raises below).
force_cov_match (bool, optional) – When True (the default – MATLAB’s default is the same True only when num_points >= 2*n; this port always defaults True and raises instead of silently falling back to False), whiten the optimized points post-hoc by Cholesky factorization so their sample covariance is exactly the identity (to float64 rounding), correcting the base algorithm’s tendency to underestimate the diagonal of the covariance (MATLAB lines 224-240). When False, the raw optimized (mirrored) points are returned with whatever sample covariance the CvM optimum happens to produce.
rng (numpy.random.Generator, optional) – Generator used to draw the MATLAB-randn-equivalent initialization matrix (see the seeding caveat above). Default None constructs a fresh numpy.random.default_rng().
max_iter (int, optional) – Maximum L-BFGS-B outer iterations, default 1000 (MATLAB’s max_iterations default, see the option-mapping notes above).

Returns:

points (ndarray) – Shape (num_points, n).
weights (ndarray) – Shape (num_points,), every entry exactly 1 / num_points (uniform weighting; MATLAB line 222), summing to 1.

Raises:

ValueError – If n < 1 or num_points < 2.
SingularMatrixError – If force_cov_match is True and the optimized points’ sample covariance is singular (MATLAB’s exitCode=-2 case, lines 230-235) – generically occurs when num_points < 2*n, since a symmetric point set spanning fewer than n free directions cannot have a full-rank covariance.
ConvergenceError – If L-BFGS-B does not report success within max_iter iterations.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

Measured on this campaign’s validation grid {(1,5),(2,10),(2,20),(3,15),(4,20)} (macOS/Apple Silicon, 2026-08-18): every case converges (result.success True) with the optimized objective strictly below its value at the random initialization, and, under force_cov_match (default), the returned points’ sample mean matches 0 and sample covariance matches the identity to float64 rounding (~1e-14 or tighter absolute, exact test tolerances recorded in tests/unit/test_lcd_samples.py). No claim is made about convergence quality outside this grid.

Examples

>>> pts, w = gaussian_lcd_samples(2, 10, rng=np.random.default_rng(0))
>>> pts.shape
(10, 2)
>>> round(float(w.sum()), 12)
1.0
>>> bool(np.allclose(pts.mean(axis=0), 0.0, atol=1e-10))
True

pytcl.mathematical_functions.numerical_integration.cube_cubature_points(n, degree, algorithm=None)[source]

Cubature points for the n-dimensional cube [-1, 1]^n.

Counterpart of the MATLAB TCL’s Cube_Space top-level, general- dimension files (see this module’s docstring for the exact per-degree algorithm coverage and the two corrected upstream defects). Every degree here is exact through that total polynomial degree – verified against the closed-form cube-monomial oracle in tests/unit/test_region_cubature.py for the (n, degree, algorithm) grid its test classes sweep; no wider claim is made (per the claims-inherit-measurement-range convention).

Weight convention (region measure, not probability). weights sum to the cube’s volume 2**n, NOT to 1 – this module targets the plain Lebesgue measure on [-1, 1]^n, unlike cubature_points’s Gaussian-weight rules. A caller wanting a probability-normalized rule divides by the volume themselves: weights / weights.sum().

Parameters:

n (int) – Dimension, n >= 1 (algorithm/degree combinations below may require more, e.g. degree 9 requires n >= 4).
degree (int) – Polynomial degree the rule is exact through. One of 1, 2, 3, 5, 7, 9 – the degrees MATLAB’s Cube_Space top-level files provide. Degrees 4, 6, 8 have no file in this directory (MATLAB provides only odd degrees plus degree 2, matching Gaussian-quadrature convention of even-order-exact rules living at the next odd degree up).
algorithm (int, optional) –
Which MATLAB algorithm variant to use; see each degree’s section in the module docstring for the ported subset (default+general-n algorithms only – fixed-n=2/n=3 literature variants are deferred, no current pytcl consumer). Default None reproduces MATLAB’s own default selection for that degree/n:
- degree 1: algorithm 1 (2^n points; algorithm 0 is 1 point).
- degree 2: algorithm 0 (n+1 points; algorithm 1 is 2n+1 points).
- degree 3: algorithm 0 (2n points, Table-I-corrected Cn 3-1).
- degree 5: algorithm 0 (2n^2+1 points, Cn 5-2).
- degree 7: algorithm 0 if n == 2, else algorithm 5 if n == 3 (no other n is supported – no general-n degree-7 cube formula exists in MATLAB).
- degree 9: algorithm 0 (n >= 4 required). Algorithm 1 is a genuinely DIFFERENT rule (MATLAB’s own comment calls it “variant 2” of the same Cn 9-1 formula, not a relabeling of algorithm 0) – both are degree-9 exact; see this module’s docstring for how they differ.

Returns:

points (ndarray) – Shape (num_points, n).
weights (ndarray) – Shape (num_points,), summing to 2**n (the cube’s volume), not 1. Commonly contains negative entries at larger n – inherent to these Stroud formulas, not suppressed. Measured examples: degree 2 algorithm 1 (every n tested, n >= 2); degree 3 algorithm 1 (n >= 4); degree 5’s DEFAULT algorithm 0 (n >= 3) and algorithms 1, 3; degree 9’s DEFAULT algorithm (n >= 4, i.e. every n this function accepts for that degree). Covariances or other quantities assembled from these points must not use a sqrt-of-weights factorization.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = cube_cubature_points(3, 3)
>>> pts.shape
(6, 3)
>>> round(float(w.sum()), 12)
8.0
>>> round(float(np.sum(w * pts[:, 0] ** 2)), 9)  # integral of x^2 over [-1,1]^3
2.666666667

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Formulas Cn 1-1/1-2, Cn 2-1/2-2, Cn 3-1/3-3/3-4/ 3-5/3-6, Cn 5-2 through 5-9, C2 7-1, C3 7-2, Cn 9-1, pp. 229-266.

R. Cools, “An encyclopedia of cubature formulas,” Journal of Complexity, vol. 19, no. 3, pp. 445-453, Jun. 2003.

pytcl.mathematical_functions.numerical_integration.simplex_cubature_points(n, degree, algorithm=None)[source]

Cubature points for the standard n-simplex {x >= 0, sum(x) <= 1}.

Counterpart of the MATLAB TCL’s Simplex top-level, general-dimension files (see this module’s docstring for the exact per-degree algorithm coverage and the corrected upstream defect). Every degree here is exact through that total polynomial degree – verified against the closed-form Dirichlet-integral simplex-monomial oracle in tests/unit/test_region_cubature.py for the (n, degree, algorithm) grid its test classes sweep; no wider claim is made (per the claims-inherit-measurement-range convention).

Weight convention (region measure, not probability). weights sum to the simplex’s volume 1 / n!, NOT to 1 – this module targets the plain Lebesgue measure on the simplex, unlike cubature_points’s Gaussian-weight rules. A caller wanting a probability-normalized rule divides by the volume themselves: weights / weights.sum().

Parameters:

n (int) – Dimension, n >= 2 (the design spec’s tested/captured dimension range starts at n=2; algorithm/degree combinations below may require more, e.g. degree 4 requires n >= 3). n=1 (the degenerate 1-D “simplex”, the interval [0, 1]) is out of scope for this port even where a formula would evaluate without error there.
degree (int) – Polynomial degree the rule is exact through. One of 2, 3, 4, 5 – the degrees MATLAB’s Simplex top-level files provide.
algorithm (int, optional) –
Which MATLAB algorithm variant to use; see each degree’s section in the module docstring for the ported subset. Default None reproduces MATLAB’s own default selection for that degree/n:
- degree 2: algorithm 0 (the only variant; MATLAB’s secondOrderSimplexCubPoints.m takes no algorithm argument).
- degree 3: algorithm 0 ((n+2) points, T_n 3-1). Algorithms 1-9 are the other general-n variants (each with its own dimension restriction, some hardened beyond MATLAB’s own checks – see module docstring); algorithms 10 (T_2 3-1, n=2 only) and 11 (T_5 3-1, n=5 only) are the fixed-dimension literature variants, not ported (module docstring).
- degree 4: algorithm 0 (the only variant; MATLAB’s fourthOrderSimplexCubPoints.m takes no algorithm argument), n >= 3 required.
- degree 5: algorithm 0 (T_n 5-2, n >= 4) if n >= 4; algorithm 1 (T_2 5-1, 7 points) if n == 2; algorithm 2 (T_3 5-1, 15 points) if n == 3 – matching MATLAB’s own per-n default selection.

Returns:

points (ndarray) – Shape (num_points, n).
weights (ndarray) – Shape (num_points,), summing to 1 / n! (the simplex’s volume), not 1. Contains negative entries for several algorithms at some dimensions – e.g. degree 3 algorithm 0’s B weight is negative for every n this function accepts – inherent to these Stroud formulas, not suppressed.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = simplex_cubature_points(2, 2)
>>> pts.shape
(6, 2)
>>> round(float(w.sum()), 12)
0.5
>>> round(float(np.sum(w * pts[:, 0])), 6)  # integral of x over the 2-simplex
0.166667

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Formulas T_n 2-2, T_n 3-1/3-2/3-3/3-4/3-5/3-7/ 3-8/3-9/3-10/3-11, T_n 4-1, T_n 5-2, T_2 5-1, T_3 5-1, pp. 307-315.

R. Cools, “An encyclopedia of cubature formulas,” Journal of Complexity, vol. 19, no. 3, pp. 445-453, Jun. 2003.

pytcl.mathematical_functions.numerical_integration.ball_cubature_points(n, degree, algorithm=None, alpha=0.0)[source]

Cubature points for the unit n-ball {x : |x| <= 1}, weight |x|**alpha.

Counterpart of the MATLAB TCL’s Sphere top-level, general-dimension files (see this module’s docstring for the exact per-degree algorithm coverage, the two excluded Sphere files with reasons, and the corrected/confirmed issues specific to this region – a documentation-only alpha-sign defect and a wrong-formula defect in degree 7 algorithm 2). Every odd degree >= 9 dispatches to a from-scratch general-order, general-alpha construction (folding in arbOrderSpherCubPoints.m, see module docstring); every other degree is exact through that total polynomial degree – verified against the closed-form ball-monomial oracle in tests/unit/test_region_cubature.py for the (n, degree, algorithm, alpha) grid its test classes sweep; no wider claim is made (per the claims-inherit-measurement-range convention).

Weight convention (region measure, not probability). weights sum to the |x|**alpha-weighted unit-ball measure 2/(n+alpha) * pi**(n/2) / gamma(n/2) (the plain unit-ball volume pi**(n/2) / gamma(n/2+1) when alpha == 0), NOT to 1 – this module targets the plain (alpha-weighted) Lebesgue measure on the unit ball, unlike cubature_points’s Gaussian-weight rules. A caller wanting a probability-normalized rule divides by the volume themselves: weights / weights.sum(). NOTE the sign of alpha here is CORRECTED relative to MATLAB’s own docstrings (which describe the weight as |x|**(-alpha)) – see the module docstring’s “confirmed MATLAB documentation defect” note; this module’s alpha matches what the MATLAB CODE actually computes, not what its comments claim.

Parameters:

n (int) – Dimension, n >= 2.
degree (int) – Polynomial degree the rule is exact through. One of 2, 3, 5, 7, or any ODD degree >= 9 (9, 11, 13, ...): 2/3/5/7 are the top-level, named-formula degrees MATLAB’s Sphere directory provides; degree 9 and above dispatch to the general-order arbOrderSpherCubPoints.m port (module docstring) at order = (degree + 1) // 2 – an EVEN degree >= 9 (e.g. 10) has no MATLAB formula in this directory at all and raises ValueError (ninthOrderSpherCubPoints.m/eleventhOrderSpherCubPoints.m are n == 2-only named files, superseded here by the general-n arbOrderSpherCubPoints path at the same degrees – see module docstring).
algorithm (int, optional) –
Which MATLAB algorithm variant to use; see each degree’s section in the module docstring for the ported subset. Default None reproduces MATLAB’s own default selection for that degree:
- degree 2: algorithm 0 (the only variant; MATLAB’s secondOrderSpherCubPoints.m takes no algorithm argument, and does not expose alpha at all – weight fixed at 1).
- degree 3: algorithm 0 (Sn 3-1, 2n points, alpha-capable). Algorithms 1 (2^n points, alpha=0), 2/3 (S2 3-1/3-2, n=2, alpha=0), 4 (S3 3-1, n=3, alpha=0) are the other ported variants.
- degree 5: algorithm 0 (Sn 5-2, 2n^2+1 points, alpha-capable). Algorithms 1 (Sn 5-3, alpha=0), 2 (Sn 5-4, alpha-capable), 3 (Sn 5-5, alpha=0), 4 (Sn 5-6, alpha=0) are general-n; 5 (S2 5-1, n=2, alpha-capable), 6 (S2 5-2, n=2, alpha=0), 7 (S3 5-1, n=3, alpha-capable), 8 (S3 5-2, n=3, alpha=0), 9 (S4 5-1, n=4, alpha=0) are the fixed-dimension variants.
- degree 7: algorithm 0 (Sn 7-2, general n>=3) always, matching MATLAB’s own unconditional default (seventhOrderSpherCubPoints.m does not auto-select per n the way fifthOrderSimplexCubPoints.m does elsewhere in this module – calling this with n == 2 and no explicit algorithm raises ValueError from algorithm 0’s own n >= 3 guard, matching real MATLAB’s behavior exactly; pass algorithm=2 explicitly for the n=2 case). No algorithm here exposes alpha at all. Algorithm 0 depends on a private port of seventhOrderSpherSurfCubPoints.m algorithm 0 (see _seventh_order_sphere_surface_alg0()); algorithms 1 (S2 7-1, n=2), 2 (S2 7-2, n=2, corrected – see module docstring), 3 (S3 7-2, n=3), 4 (S3 7-3, n=3), 5 (S4 7-2, n=4) are the fixed-dimension variants.
- degree >= 9 (odd): algorithm must be None or 0 (MATLAB’s arbOrderSpherCubPoints.m takes no algorithm argument); alpha-capable for any real alpha > -n.
alpha (float, optional) – Exponent of the radial weighting function |x|**alpha, alpha > -n. Default 0.0 (plain Lebesgue ball measure). Not every algorithm supports alpha != 0 – see the per-degree notes above; an unsupported nonzero alpha raises ValueError.

Returns:

points (ndarray) – Shape (num_points, n).
weights (ndarray) – Shape (num_points,), summing to the |x|**alpha-weighted ball measure (see above), not 1.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = ball_cubature_points(3, 3)
>>> pts.shape
(6, 3)
>>> round(float(w.sum()), 9)  # unit 3-ball volume, 4*pi/3
4.188790205
>>> round(float(np.sum(w * pts[:, 0] ** 2)), 9)  # integral of x^2, 4*pi/15
0.837758041

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Formulas Sn 2-1, Sn 3-1/3-2, S2 3-1/3-2, S3 3-1, Sn 5-2 through 5-6, S2 5-1/5-2, S3 5-1/5-2, S4 5-1, Sn 7-2, S2 7-1/7-2, S3 7-2/7-3, S4 7-2, pp. 267-292.

A. H. Stroud, “Some seventh degree integration formulas for the surface of an n-sphere,” Numerische Mathematik, vol. 11, no. 3, pp. 273-276, Mar. 1968.

R. Cools, “An encyclopedia of cubature formulas,” Journal of Complexity, vol. 19, no. 3, pp. 445-453, Jun. 2003.

pytcl.mathematical_functions.numerical_integration.spherical_surface_cubature_points(n, degree, algorithm=None)[source]

Cubature points for the unit sphere surface S^(n-1) = {x : |x| = 1}.

Counterpart of the MATLAB TCL’s Spherical_Surface top-level files (see this module’s docstring for the full per-degree algorithm coverage). Degrees 1, 3, 5, and 7 are direct ports of the matching named-formula MATLAB file; degree 14 (n=3 only) and every odd degree >= 9 REUSE existing private helpers from cubature_points rather than re-deriving those constructions (design spec Section 4, rows 179-180 of its inventory): fourteenthOrderSpherSurfCubPoints.m wraps _fourteenth_order_unit_sphere_points_3d(), and every general-n, general-odd-degree >= 9 case (superseding arbOrderSpherSurfCubPoints.m and, at n=2, arbOrder2DSpherSurfCubPoints.m) wraps _sphere_surface_points() – both private helpers normalize to sum(w) == 1, so this function rescales by the closed-form surface area (below) rather than transcribing a second, functionally-equivalent-but-differently-pointed construction. This is the one place region_cubature.py depends on cubature_points.py (one-directional; never the reverse – module docstring).

Every degree here is exact through that total polynomial degree – verified against the closed-form surface-monomial oracle (sphere_surface_monomial_integral) in tests/unit/test_region_cubature.py for the (n, degree, algorithm) grid its test classes sweep; no wider claim is made (per the claims-inherit-measurement-range convention). This also includes the two reused-construction cases: the wrapped _sphere_surface_points general path is checked at low order against the same oracle to confirm it agrees with the closed-form moments despite using a different point set than MATLAB’s own arbOrderSpherSurfCubPoints construction would (design spec’s stated purpose for that capture case).

Weight convention (region measure, not probability). weights sum to the unit sphere’s surface area 2*pi**(n/2) / gamma(n/2), NOT to 1 – this module targets the plain (uniform) surface measure on S^(n-1), unlike cubature_points’s Gaussian-weight rules (whose surface-adjacent helpers, e.g. sphere_surface_to_gauss_points(), normalize to 1 because their consumers are Kalman-family filters computing E[f(X)]). A caller wanting a probability-normalized rule divides by the surface area themselves: weights / weights.sum().

Parameters:

n (int) – Dimension, n >= 1 (individual algorithms below may require more, e.g. degree 14 requires n == 3).
degree (int) – Polynomial degree the rule is exact through. One of 1, 3, 5, 7, 14, or any ODD degree >= 9 (9, 11, 13, …): 1/3/5/7 are the top-level, named-formula degrees MATLAB’s Spherical_Surface directory provides for general n; 14 is the fixed-n=3 Stroud U3 14-1 construction (reused, see above); odd degrees >= 9 dispatch to the general-n, general-order reused construction – this also SUPERSEDES MATLAB’s ninthOrderSpherSurfCubPoints.m and eleventhOrderSpherSurfCubPoints.m (both n=3-only, Tier 3 per the design spec’s Section 8 – not individually ported here, same rationale as ball_cubature_points()’s analogous exclusion of ninthOrderSpherCubPoints.m/eleventhOrderSpherCubPoints.m). An EVEN degree other than 14 (e.g. 10) has no MATLAB formula in this directory at all and raises ValueError.
algorithm (int, optional) –
Which MATLAB algorithm variant to use; see each degree’s section in the module docstring for the ported subset. Default None reproduces MATLAB’s own default selection for that degree:
- degree 1: algorithm 0 (the only variant; MATLAB’s firstOrderSpherSurfCubPoints.m takes no algorithm argument).
- degree 3: algorithm 0 (Un 3-1, 2n points). Algorithms 1 (Un 3-2, 2^n points), 2 (Mysovskikh, 2*(n+1) points), 3 (U3 3-1, n=3, 12 points) are the other ported variants.
- degree 5: algorithm 0 (Un 5-1, 2n^2 points, negative weights for n>4). Algorithms 1-4 are general-n (Un 5-2 through Un 5-4, plus Mysovskikh); 5-8 are the n=3 fixed-dimension variants (U3 5-1 through a 30-point formula MATLAB’s own docstring mislabels – see module docstring’s fifth confirmed finding). Algorithm index 9 does not exist (MATLAB raises “Unknown algorithm specified”).
- degree 7: algorithm 0 (Formula I, general n>=3) always, matching MATLAB’s own unconditional default (mirroring ball_cubature_points()’s degree-7 dispatch note – calling this with n < 3 and no explicit algorithm raises ValueError from algorithm 0’s own guard). Algorithms 1-2 are general-n (n>=4); 3-4 are general-n (n>=3, n>=2 respectively) and are VERIFIED to compute the identical rule as algorithm 0 (see module docstring); 5 is general-n (n>=2); 6-8 are the fixed-n=3/ n=3/n=4 variants.
- degree 14: algorithm 0 (the only variant; MATLAB’s fourteenthOrderSpherSurfCubPoints.m takes no algorithm argument, and only n=3 is supported).
- degree >= 9 (odd): algorithm must be None or 0 (MATLAB’s arbOrderSpherSurfCubPoints.m takes no algorithm parameter).

Returns:

points (ndarray) – Shape (num_points, n).
weights (ndarray) – Shape (num_points,), summing to the sphere surface area (see above), not 1.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> pts, w = spherical_surface_cubature_points(3, 3)
>>> pts.shape
(6, 3)
>>> round(float(w.sum()), 9)  # surface area of S^2, 4*pi
12.566370614
>>> round(float(np.sum(w * pts[:, 0] ** 2)), 9)  # integral of x^2, 4*pi/3
4.188790205

References

A. H. Stroud, “Approximate Calculation of Multiple Integrals,” Prentice-Hall, 1971, Formulas Un 3-1/3-2, U3 3-1, Un 5-1 through 5-4, U3 5-1/5-2/5-3/5-5, Un 7-1/7-2, U3 7-1/7-2, U4 7-1, U3 14-1, pp. 292-302.

A. H. Stroud, “Some seventh degree integration formulas for the surface of an n-sphere,” Numerische Mathematik, vol. 11, no. 3, pp. 273-276, Mar. 1968.

I. P. Mysovskikh, “The approximation of multiple integrals by using interpolatory cubature formulae,” in Quantitative Approximation, R. A. DeVore and K. Scherer, eds., Academic Press, 1980, pp. 217-243.

R. Cools, “An encyclopedia of cubature formulas,” Journal of Complexity, vol. 19, no. 3, pp. 445-453, Jun. 2003.

Geometry

Geometric primitives and calculations.

This module provides: - Point-in-polygon tests - Convex hull computation - Line and plane intersections - Triangle operations - Bounding box computation

pytcl.mathematical_functions.geometry.point_in_polygon(point, polygon)[source]

Test if a point is inside a polygon.

Uses the ray casting algorithm.

Parameters:

point (array_like) – Point coordinates (x, y).
polygon (array_like) – Polygon vertices of shape (n, 2), ordered.

Returns:

inside – True if point is inside the polygon.

Return type:

bool

Examples

>>> polygon = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
>>> point_in_polygon([0.5, 0.5], polygon)
True
>>> point_in_polygon([2, 2], polygon)
False

pytcl.mathematical_functions.geometry.points_in_polygon(points, polygon)[source]

Test if multiple points are inside a polygon.

Parameters:

points (array_like) – Point coordinates of shape (n, 2).
polygon (array_like) – Polygon vertices of shape (m, 2).

Returns:

inside – Boolean array of shape (n,).

Return type:

ndarray

pytcl.mathematical_functions.geometry.convex_hull(points)[source]

Compute the convex hull of a set of points.

Parameters:

points (array_like) – Point coordinates of shape (n, d).

Returns:

vertices (ndarray) – Vertices of the convex hull.
indices (ndarray) – Indices into points of the hull vertices.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[int64]]]

Examples

>>> points = np.array([[0, 0], [1, 0], [0, 1], [0.5, 0.5]])
>>> vertices, indices = convex_hull(points)
>>> len(indices)
3

pytcl.mathematical_functions.geometry.convex_hull_area(points)[source]

Compute the area (or volume) of the convex hull.

Parameters:: points (array_like) – Point coordinates of shape (n, d).
Returns:: area – Area (2D) or volume (3D) of the convex hull.
Return type:: float

Examples

>>> points = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
>>> area = convex_hull_area(points)
>>> area
1.0

pytcl.mathematical_functions.geometry.polygon_area(vertices)[source]

Compute the area of a polygon using the shoelace formula.

Parameters:: vertices (array_like) – Polygon vertices of shape (n, 2), ordered.
Returns:: area – Signed area (positive if counterclockwise).
Return type:: float

Examples

>>> polygon_area([[0, 0], [1, 0], [1, 1], [0, 1]])
1.0

pytcl.mathematical_functions.geometry.polygon_centroid(vertices)[source]

Compute the centroid of a polygon.

Parameters:: vertices (array_like) – Polygon vertices of shape (n, 2), ordered.
Returns:: centroid – Centroid coordinates (x, y).
Return type:: ndarray

Examples

>>> polygon = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])
>>> centroid = polygon_centroid(polygon)
>>> np.allclose(centroid, [0.5, 0.5])
True

pytcl.mathematical_functions.geometry.line_intersection(p1, p2, p3, p4)[source]

Find the intersection point of two line segments.

Parameters:

p1 (array_like) – Endpoints of first line segment.
p2 (array_like) – Endpoints of first line segment.
p3 (array_like) – Endpoints of second line segment.
p4 (array_like) – Endpoints of second line segment.

Returns:

intersection – Intersection point, or None if segments don’t intersect.

Return type:

ndarray or None

Examples

>>> line_intersection([0, 0], [1, 1], [0, 1], [1, 0])
array([0.5, 0.5])

pytcl.mathematical_functions.geometry.line_plane_intersection(line_point, line_dir, plane_point, plane_normal)[source]

Find the intersection of a line and a plane.

Parameters:

line_point (array_like) – A point on the line.
line_dir (array_like) – Direction vector of the line.
plane_point (array_like) – A point on the plane.
plane_normal (array_like) – Normal vector of the plane.

Returns:

intersection – Intersection point, or None if line is parallel to plane.

Return type:

ndarray or None

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.geometry import line_plane_intersection
>>> # Line: passes through origin with direction (0, 0, 1) [vertical]
>>> line_point = np.array([0.0, 0.0, 0.0])
>>> line_dir = np.array([0.0, 0.0, 1.0])
>>> # Plane: z = 5, normal is (0, 0, 1)
>>> plane_point = np.array([0.0, 0.0, 5.0])
>>> plane_normal = np.array([0.0, 0.0, 1.0])
>>> intersection = line_plane_intersection(line_point, line_dir, plane_point, plane_normal)
>>> np.allclose(intersection, [0, 0, 5])
True
>>> # Parallel case: line and plane parallel, no intersection
>>> line_dir_parallel = np.array([1.0, 0.0, 0.0])
>>> intersection = line_plane_intersection(line_point, line_dir_parallel, plane_point, plane_normal)
>>> intersection is None
True

pytcl.mathematical_functions.geometry.point_to_line_distance(point, line_p1, line_p2)[source]

Compute the distance from a point to a line.

Parameters:

point (array_like) – Point coordinates.
line_p1 (array_like) – Two points defining the line.
line_p2 (array_like) – Two points defining the line.

Returns:

distance – Perpendicular distance from point to line.

Return type:

float

Examples

>>> point_to_line_distance([0, 1], [0, 0], [1, 0])
1.0

pytcl.mathematical_functions.geometry.point_to_line_segment_distance(point, seg_p1, seg_p2)[source]

Compute the distance from a point to a line segment.

Parameters:

point (array_like) – Point coordinates.
seg_p1 (array_like) – Endpoints of the line segment.
seg_p2 (array_like) – Endpoints of the line segment.

Returns:

distance – Distance from point to nearest point on segment.

Return type:

float

pytcl.mathematical_functions.geometry.triangle_area(p1, p2, p3)[source]

Compute the area of a triangle.

Parameters:

p1 (array_like) – Vertices of the triangle.
p2 (array_like) – Vertices of the triangle.
p3 (array_like) – Vertices of the triangle.

Returns:

area – Area of the triangle.

Return type:

float

Examples

>>> triangle_area([0, 0], [1, 0], [0, 1])
0.5

pytcl.mathematical_functions.geometry.barycentric_coordinates(point, p1, p2, p3)[source]

Compute barycentric coordinates of a point in a triangle.

Parameters:

point (array_like) – Point coordinates.
p1 (array_like) – Triangle vertices.
p2 (array_like) – Triangle vertices.
p3 (array_like) – Triangle vertices.

Returns:

coords – Barycentric coordinates (λ1, λ2, λ3) where point = λ1*p1 + λ2*p2 + λ3*p3.

Return type:

ndarray

Notes

If all coordinates are in [0, 1], the point is inside the triangle.

pytcl.mathematical_functions.geometry.delaunay_triangulation(points)[source]

Compute Delaunay triangulation.

Parameters:: points (array_like) – Point coordinates of shape (n, 2) or (n, 3).
Returns:: tri – Delaunay triangulation object. - tri.simplices: Indices of triangle vertices - tri.neighbors: Indices of neighboring triangles
Return type:: Delaunay

pytcl.mathematical_functions.geometry.bounding_box(points)[source]

Compute axis-aligned bounding box.

Parameters:

points (array_like) – Point coordinates of shape (n, d).

Returns:

min_corner (ndarray) – Minimum coordinates.
max_corner (ndarray) – Maximum coordinates.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> points = np.array([[0, 1], [2, 3], [1, 2]])
>>> min_c, max_c = bounding_box(points)
>>> min_c
array([0., 1.])
>>> max_c
array([2., 3.])

pytcl.mathematical_functions.geometry.minimum_bounding_circle(points, rng=None)[source]

Compute minimum enclosing circle (2D).

Parameters:

points (array_like) – Point coordinates of shape (n, 2).
rng (int or numpy.random.Generator, optional) – Seed or generator for the shuffle Welzl’s algorithm depends on for its expected linear running time. Default draws from fresh entropy, so repeated calls on the same points may return centers differing at the level of floating-point tie-breaking. Pass a seed for a reproducible pipeline.

Returns:

center (ndarray) – Center of the enclosing circle.
radius (float) – Radius of the enclosing circle.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], float]

Notes

Welzl’s algorithm, expected O(n).

Two robustness problems are fixed here relative to earlier versions (gh-26). The shuffle used the global np.random state, so results depended on unrelated code having drawn from it and could not be reproduced; and the implementation recursed once per point, so a few thousand points raised RecursionError. The formulation below is the standard three-loop incremental one, which recurses not at all.

The circle itself is unique, so seeding changes only which of several equally-valid representations is returned when points are co-circular – never the radius beyond floating-point noise.

Examples

>>> points = np.array([[0.0, 0.0], [2.0, 0.0], [1.0, 1.0]])
>>> center, radius = minimum_bounding_circle(points, rng=0)
>>> bool(np.isclose(radius, 1.0))
True

pytcl.mathematical_functions.geometry.oriented_bounding_box(points)[source]

Compute minimum-area oriented bounding box (2D).

Parameters:

points (array_like) – Point coordinates of shape (n, 2).

Returns:

center (ndarray) – Center of the bounding box.
extents (ndarray) – Half-widths along each principal direction.
angle (float) – Rotation angle in radians.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], float]

Combinatorics

Combinatorics utilities.

This module provides: - Permutation and combination generation - Permutation ranking/unranking - Integer partitions - Combinatorial numbers (Stirling, Bell, Catalan)

pytcl.mathematical_functions.combinatorics.factorial(n)[source]

Compute factorial of n.

Parameters:: n (int) – Non-negative integer.
Returns:: n! – Factorial of n.
Return type:: int

Examples

>>> factorial(5)
120

pytcl.mathematical_functions.combinatorics.n_choose_k(n, k)[source]

Compute binomial coefficient C(n, k).

Parameters:

n (int) – Total number of items.
k (int) – Number of items to choose.

Returns:

C(n, k) – Number of ways to choose k items from n.

Return type:

int

Examples

>>> n_choose_k(5, 2)
10

pytcl.mathematical_functions.combinatorics.n_permute_k(n, k)[source]

Compute number of k-permutations of n items.

Parameters:

n (int) – Total number of items.
k (int) – Number of items in each permutation.

Returns:

P(n, k) – Number of k-permutations: n! / (n-k)!

Return type:

int

Examples

>>> n_permute_k(5, 2)
20

pytcl.mathematical_functions.combinatorics.permutations(items, k=None)[source]

Generate all k-permutations of items.

Parameters:

items (array_like) – Items to permute.
k (int, optional) – Length of permutations. Default is len(items).

Yields:

perm (tuple) – Each k-permutation of items.

Examples

>>> list(permutations([1, 2, 3], 2))
[(1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2)]

pytcl.mathematical_functions.combinatorics.combinations(items, k)[source]

Generate all k-combinations of items.

Parameters:

items (array_like) – Items to combine.
k (int) – Size of each combination.

Yields:

comb (tuple) – Each k-combination of items.

Examples

>>> list(combinations([1, 2, 3, 4], 2))
[(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)]

pytcl.mathematical_functions.combinatorics.combinations_with_replacement(items, k)[source]

Generate all k-combinations with replacement.

Parameters:

items (array_like) – Items to combine.
k (int) – Size of each combination.

Yields:

comb (tuple) – Each k-combination with replacement.

Examples

>>> list(combinations_with_replacement([1, 2], 2))
[(1, 1), (1, 2), (2, 2)]

pytcl.mathematical_functions.combinatorics.permutation_rank(perm)[source]

Compute the lexicographic rank of a permutation.

The rank is the zero-based index of the permutation in the lexicographically sorted list of all permutations.

Parameters:: perm (array_like) – Permutation of integers 0, 1, …, n-1.
Returns:: rank – Lexicographic rank (0-indexed).
Return type:: int

Examples

>>> permutation_rank([0, 1, 2])  # First permutation
0
>>> permutation_rank([2, 1, 0])  # Last permutation
5

pytcl.mathematical_functions.combinatorics.permutation_unrank(rank, n)[source]

Compute the permutation with a given lexicographic rank.

Parameters:

rank (int) – Lexicographic rank (0-indexed).
n (int) – Length of the permutation.

Returns:

perm – Permutation of [0, 1, …, n-1] with the given rank.

Return type:

list

Examples

>>> permutation_unrank(0, 3)
[0, 1, 2]
>>> permutation_unrank(5, 3)
[2, 1, 0]

pytcl.mathematical_functions.combinatorics.next_permutation(perm)[source]

Generate the next permutation in lexicographic order.

Parameters:: perm (array_like) – Current permutation.
Returns:: next_perm – Next permutation, or None if perm is the last permutation.
Return type:: list or None

Examples

>>> next_permutation([1, 2, 3])
[1, 3, 2]
>>> print(next_permutation([3, 2, 1]))  # Last permutation
None

pytcl.mathematical_functions.combinatorics.partition_count(n, k=None)[source]

Count the number of integer partitions of n.

A partition of n is a way of writing n as a sum of positive integers, where order doesn’t matter.

Parameters:

n (int) – Number to partition.
k (int, optional) – If specified, count only partitions with exactly k parts.

Returns:

count – Number of partitions.

Return type:

int

Examples

>>> partition_count(5)  # 5 = 5 = 4+1 = 3+2 = 3+1+1 = 2+2+1 = 2+1+1+1 = 1+1+1+1+1
7
>>> partition_count(5, 2)  # 5 = 4+1 = 3+2
2

pytcl.mathematical_functions.combinatorics.partitions(n, k=None)[source]

Generate all integer partitions of n.

Parameters:

n (int) – Number to partition.
k (int, optional) – If specified, generate only partitions with exactly k parts.

Yields:

partition (tuple) – Each partition as a tuple of integers in descending order.

Examples

>>> list(partitions(4))
[(4,), (3, 1), (2, 2), (2, 1, 1), (1, 1, 1, 1)]

pytcl.mathematical_functions.combinatorics.multinomial_coefficient(*args)[source]

Compute multinomial coefficient.

multinomial(n1, n2, …, nk) = (n1 + n2 + … + nk)! / (n1! * n2! * … * nk!)

Parameters:: *args (int) – Non-negative integers.
Returns:: coeff – Multinomial coefficient.
Return type:: int

Examples

>>> multinomial_coefficient(2, 3, 1)  # 6! / (2! * 3! * 1!)
60

pytcl.mathematical_functions.combinatorics.stirling_second(n, k)[source]

Stirling number of the second kind.

S(n, k) is the number of ways to partition n elements into k non-empty subsets.

Parameters:

n (int) – Number of elements.
k (int) – Number of subsets.

Returns:

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

Return type:

int

Examples

>>> stirling_second(4, 2)  # {{1,2,3},{4}}, {{1,2,4},{3}}, ...
7

pytcl.mathematical_functions.combinatorics.bell_number(n)[source]

Bell number B_n.

B_n is the number of ways to partition a set of n elements.

Parameters:: n (int) – Number of elements.
Returns:: B_n – n-th Bell number.
Return type:: int

Examples

>>> bell_number(4)
15

pytcl.mathematical_functions.combinatorics.catalan_number(n)[source]

Catalan number C_n.

Catalan numbers count many combinatorial structures including: - Valid parenthesizations - Full binary trees with n+1 leaves - Triangulations of a polygon with n+2 sides

Parameters:: n (int) – Non-negative integer.
Returns:: C_n – n-th Catalan number.
Return type:: int

Examples

>>> catalan_number(5)
42

pytcl.mathematical_functions.combinatorics.derangements_count(n)[source]

Count the number of derangements.

A derangement is a permutation with no fixed points.

Parameters:: n (int) – Number of elements.
Returns:: D_n – Number of derangements.
Return type:: int

Examples

>>> derangements_count(4)  # {2,1,4,3}, {2,3,4,1}, ...
9

pytcl.mathematical_functions.combinatorics.subfactorial(n)[source]

Subfactorial (number of derangements).

Alias for derangements_count.

Parameters:: n (int) – Number of elements.
Returns:: !n – Subfactorial of n.
Return type:: int

Polynomials

Simultaneous multivariate polynomial root finding via the Macaulay null-space method.

Simultaneous multivariate polynomial root finding.

Port of the MATLAB TCL polyRootsMultiDim.m and the helpers it draws from (multiDimPolyMat2Terms, rankTComposition, unrankTComposition, nullspace), implementing the affine null-space Macaulay-matrix method of [1] (Algorithm 3): the root finding problem becomes a generalized eigenvalue problem on the null space of a degree-augmented Macaulay matrix.

References

class pytcl.mathematical_functions.polynomials.multivariate.PolyRootsResult(roots, exit_code)[source]

Bases: NamedTuple

Result of poly_roots_multi_dim().

Variables:

roots (ndarray) – (n, num_sol) matrix of the affine roots found (complex in general), or an empty (n, 0) array when exit_code is nonzero.
exit_code (int) – 0 on success; 1 if the maximum number of degree increases elapsed; 2 if a finite-precision error made the Macaulay nullity change after stabilizing or decrease with degree.

roots: ndarray[tuple[Any, ...], dtype[complexfloating]]: Alias for field number 0

exit_code: int: Alias for field number 1

pytcl.mathematical_functions.polynomials.multivariate.poly_roots_multi_dim(poly_coeff_mats, max_deg_increases=None, use_motzkin_null=False)[source]

Roots of a system of simultaneous multivariate polynomials.

Only the affine roots are found (generally the only ones desired), not the roots at infinity. Due to finite-precision effects and the combinatorial growth of the Macaulay matrix, the method is best suited to systems of at most 3 variables and degree at most 3; sparse systems fare much better than dense ones.

Parameters:

poly_coeff_mats (sequence of array_like) – n coefficient hypermatrices, one per polynomial in n variables. coeffs[a1, a2, ..., an] is the coefficient of x1**a1 * x2**a2 * ... * xn**an (note: zero-based exponents, the reverse of MATLAB’s 1-based indices with the same layout).
max_deg_increases (int, optional) – Maximum number of degree increases of the Macaulay matrix. Too small a value makes the solve fail with exit code 1. Default: 10 * n.
use_motzkin_null (bool, optional) – Use the Motzkin null-space algorithm of [1] instead of the SVD. Generally less numerically stable; provided to allow stepping through the reference values in [1]. Default False.

Returns:

result – The roots as an (n, num_sol) complex matrix and the exit code.

Return type:

PolyRootsResult

Examples

The two-variable system from Section 2.1 of [1], whose four roots are all real: (4, -5), (1, 0), (3, -2) and (0, -1).

>>> import numpy as np
>>> p = np.zeros((3, 3))
>>> p[0, 0], p[2, 0], p[1, 1], p[0, 2] = -4.0, -1.0, 2.0, 1.0
>>> p[1, 0], p[0, 1] = 5.0, -3.0
>>> q = np.zeros((3, 3))
>>> q[0, 0], q[2, 0], q[1, 1], q[0, 2] = -1.0, 1.0, 2.0, 1.0
>>> roots, exit_code = poly_roots_multi_dim([p, q])
>>> exit_code
0
>>> sorted(np.round(roots.real.T, 6).tolist())
[[0.0, -1.0], [1.0, 0.0], [3.0, -2.0], [4.0, -5.0]]

Notes

Port of polyRootsMultiDim.m, implementing Algorithm 3 of [1]. The shift function g(x) is the arbitrary choice sum_i i * x_i made by the original. Monomials are tracked with composition ranking/unranking; the Macaulay matrix’s sparsity is not exploited, as in the original.

Debye

Debye functions.

Debye functions appear in solid-state physics for computing thermodynamic properties of solids (heat capacity, entropy).

Performance

This module uses Numba JIT compilation with rapidly convergent series expansions (Abramowitz & Stegun 27.1.1-27.1.3), providing ~10-50x speedup for batch computations compared to scipy.integrate.quad.

Accuracy

Typical relative error is ~1e-16 (machine precision), measured against a 30-40-digit mpmath oracle for n in {1, 4, 6, 8, 9, 10}, x in {0.5, 0.99, 1.0, 1.01, 1.9, 1.99, 1.999, 2.0, 10.0}. The worst points on that grid are ~2e-11 just below the x0=2.0 boundary (n=9/10 near x=1.99-1.999, where the small-x series is nearing the edge of its useful precision) and ~4e-12 exactly at the boundary (n=9/10 at x=2.0, where the large-x branch’s own alternating-sum cancellation happens to be smallest). The small-x/large-x branch switch is at x0=2.0: before this was tightened from x0=1.0, the same grid measured up to 8.6e-9 at x=1.0-1.01 for n>=8 (catastrophic cancellation in the large-x branch right at the old boundary).

pytcl.mathematical_functions.special_functions.debye.debye(n, x)[source]

Debye function D_n(x).

The Debye function of order n is defined as: D_n(x) = (n/x^n) * integral from 0 to x of t^n / (exp(t) - 1) dt

Parameters:

n (int) – Order of the Debye function (positive integer).
x (array_like) – Argument of the function, x >= 0.

Returns:

D – Values of D_n(x).

Return type:

ndarray

Notes

Special cases: - D_n(0) = 1 - D_n(inf) = n! * zeta(n+1) / x^n -> 0

The Debye function D_3(x) appears in the heat capacity of solids at low temperatures.

This implementation uses Numba JIT compilation for performance, achieving ~10-50x speedup compared to scipy.integrate.quad for batch computations.

Examples

>>> float(debye(3, 0)[0])  # D_3(0) = 1
1.0
>>> round(float(debye(3, 1)[0]), 6)
0.674416
>>> round(float(debye(3, 10)[0]), 6)
0.019296

References

Debye, P. (1912). “Zur Theorie der spezifischen Wärmen”. Annalen der Physik, 344(14), 789-839.

pytcl.mathematical_functions.special_functions.debye.debye_1(x)[source]

First-order Debye function D_1(x).

Parameters:: x (array_like) – Argument of the function, x >= 0.
Returns:: D – Values of D_1(x).
Return type:: ndarray

Notes

D_1(x) = (1/x) * integral from 0 to x of t / (exp(t) - 1) dt

pytcl.mathematical_functions.special_functions.debye.debye_2(x)[source]

Second-order Debye function D_2(x).

Parameters:: x (array_like) – Argument of the function, x >= 0.
Returns:: D – Values of D_2(x).
Return type:: ndarray

Notes

D_2(x) = (2/x^2) * integral from 0 to x of t^2 / (exp(t) - 1) dt

pytcl.mathematical_functions.special_functions.debye.debye_3(x)[source]

Third-order Debye function D_3(x).

This is the most commonly used Debye function, appearing in the heat capacity of solids.

Parameters:: x (array_like) – Argument of the function, x >= 0.
Returns:: D – Values of D_3(x).
Return type:: ndarray

Notes

D_3(x) = (3/x^3) * integral from 0 to x of t^3 / (exp(t) - 1) dt

The heat capacity of a solid in the Debye model is: C_V = 9 * N * k_B * (T/Θ_D)^3 * D_3(Θ_D/T)

where Θ_D is the Debye temperature.

pytcl.mathematical_functions.special_functions.debye.debye_4(x)[source]

Fourth-order Debye function D_4(x).

Parameters:: x (array_like) – Argument of the function, x >= 0.
Returns:: D – Values of D_4(x).
Return type:: ndarray

Notes

D_4(x) = (4/x^4) * integral from 0 to x of t^4 / (exp(t) - 1) dt

This appears in computing the entropy of solids.

pytcl.mathematical_functions.special_functions.debye.debye_heat_capacity(temperature, debye_temperature)[source]

Debye model heat capacity (normalized).

Computes C_V / (3*N*k_B) using the Debye model.

Parameters:

temperature (array_like) – Temperature in Kelvin.
debye_temperature (float) – Debye temperature Θ_D in Kelvin.

Returns:

cv_normalized – Normalized heat capacity C_V / (3*N*k_B). Multiply by 3*N*k_B for actual heat capacity.

Return type:

ndarray

Notes

The Debye model heat capacity is: C_V / (3*N*k_B) = 4*D_3(x) - 3*x/(e^x - 1), with x = Θ_D/T

Limits: - High T (T >> Θ_D): C_V -> 3*N*k_B (classical) - Low T (T << Θ_D): C_V ~ (4*π^4/5) * (T/Θ_D)^3 (quantum)

Examples

>>> # Aluminum at room temperature (Θ_D ≈ 428 K)
>>> cv = debye_heat_capacity(300, 428)  # ~0.91

pytcl.mathematical_functions.special_functions.debye.debye_entropy(temperature, debye_temperature)[source]

Debye model entropy (normalized).

Computes S / (3*N*k_B) using the Debye model.

Parameters:

temperature (array_like) – Temperature in Kelvin.
debye_temperature (float) – Debye temperature Θ_D in Kelvin.

Returns:

s_normalized – Normalized entropy S / (3*N*k_B).

Return type:

ndarray

Notes

The entropy in the Debye model is: S / (3*N*k_B) = (4/3)*D_3(Θ_D/T) - ln(1 - exp(-Θ_D/T))

Hypergeometric

Hypergeometric functions.

This module provides hypergeometric functions commonly used in mathematical physics, probability theory, and special function evaluation.

Performance

The generalized hypergeometric function uses Numba JIT compilation for the series summation loop, providing significant speedup for the general case (p > 2 or q > 1).

pytcl.mathematical_functions.special_functions.hypergeometric.hyp0f1(b, z)[source]

Confluent hypergeometric limit function 0F1(b; z).

The function 0F1(b; z) is defined by the series: 0F1(b; z) = sum_{k=0}^inf z^k / ((b)_k * k!)

where (b)_k is the Pochhammer symbol (rising factorial).

Parameters:

b (array_like) – Numerator parameter. Must not be a non-positive integer.
z (array_like) – Argument of the function.

Returns:

F – Values of 0F1(b; z).

Return type:

ndarray

Notes

Related to Bessel functions: J_n(x) = (x/2)^n / Gamma(n+1) * 0F1(n+1; -x^2/4) I_n(x) = (x/2)^n / Gamma(n+1) * 0F1(n+1; x^2/4)

Examples

>>> float(hyp0f1(1, 0))  # 0F1(1; 0) = 1
1.0
>>> round(float(hyp0f1(1, 1)), 6)
2.279585

References

NIST Digital Library of Mathematical Functions, Chapter 16.

pytcl.mathematical_functions.special_functions.hypergeometric.hyp1f1(a, b, z)[source]

Confluent hypergeometric function 1F1(a; b; z) (Kummer’s function M).

The function 1F1(a; b; z) is defined by the series: 1F1(a; b; z) = sum_{k=0}^inf (a)_k * z^k / ((b)_k * k!)

Parameters:

a (array_like) – Numerator parameter.
b (array_like) – Denominator parameter. Must not be a non-positive integer.
z (array_like) – Argument of the function.

Returns:

F – Values of 1F1(a; b; z).

Return type:

ndarray

Notes

Also known as Kummer’s function M(a, b, z).

Special cases: - 1F1(0; b; z) = 1 - 1F1(a; a; z) = exp(z) - 1F1(1; 2; 2z) = sinh(z) * exp(z) / z

Related to incomplete gamma: gammainc(a, z) = z^a * exp(-z) * 1F1(1; 1+a; z) / (a * Gamma(a))

Examples

>>> round(float(hyp1f1(1, 1, 1)), 6)  # exp(1)
2.718282
>>> round(float(hyp1f1(0.5, 1.5, -1)), 6)  # erf(1) * sqrt(pi) / 2
0.746824

References

Abramowitz & Stegun, “Handbook of Mathematical Functions”, Ch. 13.

pytcl.mathematical_functions.special_functions.hypergeometric.hyp2f1(a, b, c, z)[source]

Gauss hypergeometric function 2F1(a, b; c; z).

The function 2F1(a, b; c; z) is defined by the series:

2F1(a, b; c; z) = sum_{k=0}^inf (a)_k * (b)_k * z^k / ((c)_k * k!)

converging for |z| < 1.

Parameters:

a (array_like) – First numerator parameter.
b (array_like) – Second numerator parameter.
c (array_like) – Denominator parameter. Must not be a non-positive integer.
z (array_like) – Argument of the function. For |z| >= 1, analytic continuation is used.

Returns:

F – Values of 2F1(a, b; c; z).

Return type:

ndarray

Notes

Many elementary and special functions are special cases: - (1-z)^(-a) = 2F1(a, b; b; z) - log(1+z)/z = 2F1(1, 1; 2; -z) - arcsin(z)/z = 2F1(1/2, 1/2; 3/2; z^2) - Complete elliptic integrals K(k) and E(k)

Examples

>>> round(float(hyp2f1(1, 1, 2, 0.5)), 6)  # -log(1-0.5)/0.5 = 2*log(2)
1.386294
>>> round(float(hyp2f1(0.5, 0.5, 1.5, 0.25)), 6)  # arcsin(0.5)/0.5 = pi/3
1.047198

References

NIST DLMF, Chapter 15.

pytcl.mathematical_functions.special_functions.hypergeometric.hyperu(a, b, z)[source]

Confluent hypergeometric function U(a, b, z) (Tricomi function).

The function U(a, b, z) is defined as:

U(a, b, z) = Gamma(1-b)/Gamma(a-b+1) * 1F1(a; b; z)
             + Gamma(b-1)/Gamma(a) * z^(1-b) * 1F1(a-b+1; 2-b; z)

Parameters:

a (array_like) – First parameter.
b (array_like) – Second parameter.
z (array_like) – Argument of the function (must be positive for real result).

Returns:

U – Values of U(a, b, z).

Return type:

ndarray

Notes

Also known as Tricomi’s function or Kummer’s function of the second kind.

Asymptotic behavior for large z: U(a, b, z) ~ z^(-a) as z -> inf

Examples

>>> round(float(hyperu(1, 1, 1)), 6)
0.596347

pytcl.mathematical_functions.special_functions.hypergeometric.hyp1f1_regularized(a, b, z)[source]

Regularized confluent hypergeometric function 1F1(a; b; z) / Gamma(b).

This is useful when b may be near a non-positive integer.

Parameters:

a (array_like) – Numerator parameter.
b (array_like) – Denominator parameter.
z (array_like) – Argument of the function.

Returns:

F – Values of 1F1(a; b; z) / Gamma(b).

Return type:

ndarray

Examples

>>> import numpy as np
>>> from pytcl.mathematical_functions.special_functions import hyp1f1_regularized
>>> # Regularized form avoids overflow for problematic b values
>>> a, b, z = 0.5, 1.5, 1.0
>>> f_reg = hyp1f1_regularized(a, b, z)
>>> # Should give finite, non-overflowing result
>>> np.isfinite(f_reg)
True
>>> # Compare to regular hypergeometric computation
>>> from pytcl.mathematical_functions.special_functions import hyp1f1
>>> import scipy.special as sp
>>> f_normal = hyp1f1(a, b, z) / sp.gamma(b)
>>> np.allclose(f_reg, f_normal)
True

Notes

This function remains finite even when b is a non-positive integer, unlike the standard 1F1.

pytcl.mathematical_functions.special_functions.hypergeometric.pochhammer(a, n)[source]

Pochhammer symbol (rising factorial) (a)_n.

The Pochhammer symbol is defined as: (a)_n = a * (a+1) * (a+2) * … * (a+n-1) = Gamma(a+n) / Gamma(a)

Parameters:

a (array_like) – Base value.
n (array_like) – Number of terms (can be non-integer for generalization).

Returns:

p – Values of (a)_n.

Return type:

ndarray

Notes

Special cases: - (a)_0 = 1 - (1)_n = n! - (a)_1 = a

Examples

>>> float(pochhammer(1, 5))  # 5!
120.0
>>> float(pochhammer(3, 4))  # 3*4*5*6
360.0

pytcl.mathematical_functions.special_functions.hypergeometric.falling_factorial(a, n)[source]

Falling factorial (a)_n (Pochhammer symbol variant).

The falling factorial is defined as: (a)_n = a * (a-1) * (a-2) * … * (a-n+1)

Parameters:

a (array_like) – Base value.
n (array_like) – Number of terms.

Returns:

f – Values of the falling factorial.

Return type:

ndarray

Notes

Related to rising factorial: (a)_n (falling) = (-1)^n * (-a)_n (rising)

Examples

>>> float(falling_factorial(5, 3))  # 5*4*3
60.0

pytcl.mathematical_functions.special_functions.hypergeometric.generalized_hypergeometric(a, b, z, max_terms=500, tol=1e-15)[source]

Generalized hypergeometric function pFq(a; b; z).

Computes the generalized hypergeometric function with p numerator and q denominator parameters.

Parameters:

a (array_like) – Numerator parameters (1D array of length p).
b (array_like) – Denominator parameters (1D array of length q).
z (array_like) – Argument of the function.
max_terms (int, optional) – Maximum number of series terms. Default is 500.
tol (float, optional) – Tolerance for series convergence. Default is 1e-15.

Returns:

F – Values of pFq(a; b; z).

Return type:

ndarray

Notes

The series converges for:

- p <= q: all z
- p = q + 1: |z| < 1
- p > q + 1: diverges except for polynomial cases

Uses Numba JIT compilation for the general case (p > 2 or q > 1), providing 5-10x speedup over pure Python loops.

Examples

>>> round(float(generalized_hypergeometric([1], [2], 1)), 6)  # 1F1(1; 2; 1) = e - 1
1.718282

Lambert W

Lambert W function and related functions.

The Lambert W function appears in diverse applications including delay differential equations, combinatorics, and physics.

pytcl.mathematical_functions.special_functions.lambert_w.lambert_w(z, k=0, tol=1e-10)[source]

Lambert W function W_k(z).

The Lambert W function is defined as the inverse of f(w) = w * exp(w), satisfying W(z) * exp(W(z)) = z.

Parameters:

z (array_like) – Argument of the function. Can be complex.
k (int, optional) – Branch index. Default is 0 (principal branch). - k = 0: Principal branch, real for z >= -1/e - k = -1: Lower branch, real for -1/e <= z < 0 - Other k: Complex branches
tol (float, optional) – Tolerance for convergence (used in edge cases). Default is 1e-10.

Returns:

W – Values of W_k(z).

Return type:

ndarray

Notes

The principal branch W_0(z) satisfies: - W_0(0) = 0 - W_0(e) = 1 - W_0(-1/e) = -1

The function has a branch point at z = -1/e ≈ -0.3679.

Examples

>>> float(lambert_w(0).real)  # W(0) = 0
0.0
>>> float(lambert_w(np.e).real)  # W(e) = 1
1.0
>>> float(lambert_w(-np.exp(-1)).real)  # W(-1/e) = -1
-1.0

References

Corless, R.M., et al. (1996). “On the Lambert W Function”. Advances in Computational Mathematics, 5, 329-359.

pytcl.mathematical_functions.special_functions.lambert_w.lambert_w_real(x, branch=0)[source]

Real-valued Lambert W function.

Returns only the real part of the Lambert W function for real inputs.

Parameters:

x (array_like) – Real argument. For branch 0: x >= -1/e. For branch -1: -1/e <= x < 0.
branch (int, optional) – Branch index: 0 (principal) or -1 (lower). Default is 0.

Returns:

W – Real values of W(x).

Return type:

ndarray

Raises:

ValueError – If x is outside the valid range for real-valued output.

Examples

>>> round(float(lambert_w_real(1)), 6)
0.567143
>>> round(float(lambert_w_real(-0.2, branch=-1)), 6)
-2.542641

pytcl.mathematical_functions.special_functions.lambert_w.omega_constant()[source]

Omega constant (principal value of W(1)).

The omega constant Ω is the unique real solution to Ω * exp(Ω) = 1, satisfying Ω = W_0(1).

Returns:: omega – Ω ≈ 0.5671432904097838729999686622…
Return type:: float

Notes

The omega constant appears in: - Growth of the iterated logarithm - Stirling’s approximation refinements - Analysis of tree structures

Examples

>>> omega = omega_constant()
>>> round(omega * np.exp(omega), 12)  # Should equal 1
1.0

pytcl.mathematical_functions.special_functions.lambert_w.wright_omega(z)[source]

Wright omega function ω(z).

The Wright omega function is defined as ω(z) = W_k(e^z) for the appropriate branch k.

Parameters:: z (array_like) – Argument of the function. Can be complex.
Returns:: omega – Values of the Wright omega function.
Return type:: ndarray

Notes

The Wright omega function satisfies: ω(z) + log(ω(z)) = z

It is entire (analytic everywhere) unlike the Lambert W function.

Examples

>>> round(float(wright_omega(0).real), 6)  # Omega constant
0.567143

References

Wright, E.M. (1959). “Solution of the equation z*exp(z) = a”. Bull. Amer. Math. Soc., 65, 89-93.

pytcl.mathematical_functions.special_functions.lambert_w.solve_exponential_equation(a, b, c)[source]

Solve a*x*exp(b*x) = c using Lambert W.

Finds x such that a*x*exp(b*x) = c.

Parameters:

a (array_like) – Coefficient of x.
b (array_like) – Coefficient in the exponential.
c (array_like) – Right-hand side constant.

Returns:

x – Solution(s) to the equation.

Return type:

ndarray

Notes

The solution is: x = W(b*c/a) / b

Examples

>>> x = solve_exponential_equation(1, 1, np.e)  # x*exp(x) = e
>>> round(float(x.real), 6)  # Should be 1
1.0

pytcl.mathematical_functions.special_functions.lambert_w.time_delay_equation(a, tau)[source]

Solve characteristic equation for first-order delay system.

Finds s such that s + a*exp(-s*tau) = 0, which appears in delay differential equations.

Parameters:

a (array_like) – Coefficient in the characteristic equation.
tau (array_like) – Time delay.

Returns:

s – Root(s) of the characteristic equation.

Return type:

ndarray

Notes

The solution is: s = W(a*tau)/tau

This is the dominant eigenvalue for the delay system: dx/dt = -a * x(t - tau)

Examples

>>> s = time_delay_equation(1, 1)  # s + exp(-s) = 0
>>> bool(abs(s + np.exp(-s)) < 1e-10)  # Should be approximately 0
True

Marcum Q

Marcum Q function and related functions.

The Marcum Q function is crucial in radar and communications for analyzing detection probabilities and signal statistics.

pytcl.mathematical_functions.special_functions.marcum_q.marcum_q(a, b, m=1)[source]

Generalized Marcum Q function Q_m(a, b).

The Marcum Q function is the complementary cumulative distribution function of the noncentral chi-squared distribution and appears in radar detection theory.

Parameters:

a (array_like) – First argument (non-centrality parameter), a >= 0.
b (array_like) – Second argument (threshold), b >= 0.
m (int, optional) – Order of the Marcum Q function (positive integer). Default is 1.

Returns:

Q – Values of Q_m(a, b).

Return type:

ndarray

Notes

For m = 1, this is the standard Marcum Q function: Q_1(a, b) = integral from b to inf of x * exp(-(x^2 + a^2)/2) * I_0(a*x) dx

The function is related to the noncentral chi-squared distribution: Q_m(a, b) = P(X > b^2) where X ~ chi^2(2m, a^2)

Special cases: - Q_m(0, b) = 1 - gammainc(m, b^2/2) = gammaincc(m, b^2/2) - Q_m(a, 0) = 1

Examples

>>> float(marcum_q(0, 0))  # Q_1(0, 0) = 1
1.0
>>> round(float(marcum_q(3, 4)), 6)  # Standard Marcum Q
0.196512

References

Marcum, J.I. (1950). “Table of Q Functions”.
Shnidman, D.A. (1989). “The Calculation of the Probability of Detection and the Generalized Marcum Q-Function”. IEEE Trans. on Information Theory, 35(2), 389-400.

pytcl.mathematical_functions.special_functions.marcum_q.marcum_q1(a, b)[source]

Standard Marcum Q function Q_1(a, b).

Convenience function for the first-order Marcum Q function.

Parameters:

a (array_like) – First argument (non-centrality parameter), a >= 0.
b (array_like) – Second argument (threshold), b >= 0.

Returns:

Q – Values of Q_1(a, b).

Return type:

ndarray

Examples

>>> round(float(marcum_q1(2, 2)), 6)
0.603501

See also

marcum_q: Generalized Marcum Q function.

pytcl.mathematical_functions.special_functions.marcum_q.log_marcum_q(a, b, m=1)[source]

Natural logarithm of the Marcum Q function.

Computes log(Q_m(a, b)) with better numerical precision for small values of Q.

Parameters:

a (array_like) – First argument (non-centrality parameter), a >= 0.
b (array_like) – Second argument (threshold), b >= 0.
m (int, optional) – Order of the Marcum Q function. Default is 1.

Returns:

log_Q – Values of log(Q_m(a, b)).

Return type:

ndarray

Notes

For small Q values (large b), this function provides better numerical accuracy than computing log(marcum_q(a, b)).

Examples

>>> round(float(log_marcum_q(1, 5)), 6)  # log(Q_1(1, 5))
-9.506564

pytcl.mathematical_functions.special_functions.marcum_q.marcum_q_inv(a, q, m=1)[source]

Inverse Marcum Q function.

Finds b such that Q_m(a, b) = q.

Parameters:

a (array_like) – First argument (non-centrality parameter), a >= 0.
q (array_like) – Target probability value, 0 < q < 1.
m (int, optional) – Order of the Marcum Q function. Default is 1.

Returns:

b – Values such that Q_m(a, b) = q.

Return type:

ndarray

Notes

Closed form via the noncentral chi-squared inverse survival function: b = sqrt(ncx2.isf(q, 2m, a**2)). No iteration is involved – the previous tol/max_iter parameters were accepted and never read, alongside a Notes claim of Newton-Raphson iteration that no code performed.

Examples

>>> b = marcum_q_inv(3, 0.5)  # Find b where Q_1(3, b) = 0.5
>>> round(float(marcum_q(3, b)), 6)  # Verify
0.5

pytcl.mathematical_functions.special_functions.marcum_q.rician_cdf(a, b)[source]

Rician cumulative distribution function, 1 - Q_1(a, b).

Parameters:

a (array_like) – Non-centrality parameter, a >= 0.
b (array_like) – Threshold, b >= 0.

Returns:

P – Values of 1 - Q_1(a, b).

Return type:

ndarray

Notes

This is the probability P(X <= b^2) for X ~ chi^2(2, a^2).

Formerly exported as nuttall_q, which was a misnomer: the Nuttall Q function Q_{m,n}(a, b) is a different integral, a generalization of the Marcum Q with an extra power of the integration variable. This routine computes neither – it is the complementary Marcum Q, which is exactly the Rician CDF, and it always did so correctly (gh-20). Only the name was wrong. The deprecated nuttall_q alias was removed in v2.8.0.

Examples

>>> round(float(rician_cdf(2, 2)), 6)  # 1 - Q_1(2, 2)
0.396499

See also

marcum_q: Marcum Q function.

pytcl.mathematical_functions.special_functions.marcum_q.swerling_detection_probability(snr, pfa, n_pulses=1, swerling_case=1)[source]

Detection probability for Swerling target models.

Computes probability of detection for different Swerling cases using the Marcum Q function.

Parameters:

snr (array_like) – Signal-to-noise ratio (linear, not dB).
pfa (float) – Probability of false alarm (0 < pfa < 1).
n_pulses (int, optional) – Number of integrated pulses. Default is 1.
swerling_case (int, optional) – Swerling case (0, 1, 2, 3, or 4). Default is 1. - 0: Non-fluctuating (Marcum) - 1: Slow fluctuation, Rayleigh PDF - 2: Fast fluctuation, Rayleigh PDF - 3: Slow fluctuation, one dominant + Rayleigh - 4: Fast fluctuation, one dominant + Rayleigh

Returns:

Pd – Probability of detection.

Return type:

ndarray

Notes

The detection threshold T is set from the false alarm probability via pfa = gammaincc(n, T/2) (square-law detector, n integrated pulses).

For Swerling 0 (non-fluctuating):: P_d = Q_n(sqrt(2*n*SNR), sqrt(T))

Swerling 1 and 2 use the exact closed forms for chi-squared (2 DOF) target fluctuation with scan-to-scan (1) or pulse-to-pulse (2) decorrelation. Swerling 3 uses the DiFranco-Rubin closed form for chi-squared (4 DOF) scan-to-scan fluctuation, and Swerling 4 the exact finite-sum for pulse-to-pulse chi-squared (4 DOF) fluctuation.

Examples

>>> pd = swerling_detection_probability(10, 1e-6, n_pulses=10, swerling_case=0)
>>> pd > 0.9  # High probability of detection with 10 dB SNR
True

References

Swerling, P. (1960). “Probability of Detection for Fluctuating Targets”. IRE Trans. on Information Theory, IT-6, 269-308.

algorithm	n=1	n=2	n=3	n=4	n=5	n=6	n=7	n=8
0 (max 17)	16	15	15	15	14	14	14	14 (b)
1 (max 15)	14	15 (a)	13	13	13	13	13	13 (c)