Mathematical Functions

Mathematical functions and utilities.

This module contains a wide variety of mathematical functions including: - Basic matrix operations - Combinatorics (permutations, combinations) - Continuous optimization - Geometry primitives - Interpolation methods - Numerical integration - Polynomials - 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.
triaSqrt: 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

Uses the recurrence relation for numerical stability: J_{n+1}(x) / J_n(x) = 2n/x - 1/(J_n(x)/J_{n-1}(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.

Return type:

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.

Return type:

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.

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

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

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 high accuracy (~1e-14 relative) and ~10-50x speedup for batch computations compared to scipy.integrate.quad.

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, tol=1e-10, max_iter=100)[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.
tol (float, optional) – Tolerance for convergence. Default is 1e-10.
max_iter (int, optional) – Maximum number of iterations. Default is 100.

Returns:

b – Values such that Q_m(a, b) = q.

Return type:

ndarray

Notes

Uses Newton-Raphson iteration with the noncentral chi-squared distribution.

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.nuttall_q(a, b)[source]

Deprecated alias for rician_cdf().

The name was wrong: this computes 1 - Q_1(a, b), the Rician CDF, not the Nuttall Q function Q_{m,n}(a, b), which is a different integral (gh-20). The computation was always correct.

Deprecated since version Use: rician_cdf(). This alias will be removed in a future release.

Examples

>>> import warnings
>>> with warnings.catch_warnings():
...     warnings.simplefilter("ignore", DeprecationWarning)
...     round(float(nuttall_q(2, 2)), 6)
0.396499

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. nuttall_q remains as a deprecated alias.

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.