Mathematical Functions<a class="headerlink" href="#module-pytcl.mathematical_functions" title="Link to this heading">

Examples

>>> result, error = quad(lambda x: x**2, 0, 1)
>>> result
0.33333333333333337

See also

scipy.integrate.quad: Underlying implementation.

pytcl.mathematical_functions.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,)
>>> np.sum(weights)  # Weights sum to 1
1.0

pytcl.mathematical_functions.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

pytcl.mathematical_functions.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.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

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

See also

numpy.interp: Underlying implementation.

pytcl.mathematical_functions.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)
>>> cs(np.pi/2)
array(0.99999...)

See also

scipy.interpolate.CubicSpline: Underlying implementation.

pytcl.mathematical_functions.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.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.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.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:

Examples

>>> n_choose_k(5, 2)
10

pytcl.mathematical_functions.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.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.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.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.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.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.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.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.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.butter_design(order, cutoff, fs, btype='low', output='sos')[source]

Design a Butterworth digital filter.

The Butterworth filter has maximally flat frequency response in the passband. It is often used as a “standard” filter when sharp transitions are not required.

Parameters:

order (int) – Filter order.
cutoff (float or tuple) – Cutoff frequency in Hz. For bandpass/bandstop, provide (low, high).
fs (float) – Sampling frequency in Hz.
btype ({'low', 'high', 'band', 'bandstop'}, optional) – Filter type. Default is ‘low’.
output ({'sos', 'ba'}, optional) – Output format. ‘sos’ is recommended for stability. Default is ‘sos’.

Returns:

coeffs – Filter coefficients (b, a, sos).

Return type:

FilterCoefficients

Examples

>>> fs = 1000  # 1 kHz
>>> coeffs = butter_design(4, 100, fs, btype='low')
>>> coeffs.sos.shape[0]  # Number of second-order sections
2

Notes

The Butterworth filter has no ripple in the passband or stopband. The -3 dB point occurs at the cutoff frequency.

pytcl.mathematical_functions.cfar_ca(signal, guard_cells, ref_cells, pfa=1e-06, alpha=None)[source]

Cell-Averaging CFAR detector.

CA-CFAR estimates the noise level by averaging the cells in the reference window (excluding guard cells around the cell under test).

Parameters:

signal (array_like) – Input signal (typically power or magnitude).
guard_cells (int) – Number of guard cells on each side of the cell under test.
ref_cells (int) – Number of reference cells on each side.
pfa (float, optional) – Probability of false alarm. Default is 1e-6.
alpha (float, optional) – Threshold multiplier. If provided, overrides pfa calculation.

Returns:

result – Named tuple with detections, threshold, indices, and noise estimate.

Return type:

CFARResult

Examples

>>> import numpy as np
>>> np.random.seed(42)
>>> # Noise with a few targets
>>> signal = np.random.exponential(1.0, 1000)
>>> signal[250] = 50  # Target 1
>>> signal[500] = 100  # Target 2
>>> signal[750] = 30  # Target 3
>>> result = cfar_ca(signal, guard_cells=2, ref_cells=16, pfa=1e-4)
>>> 250 in result.detection_indices
True

Notes

The CA-CFAR is optimal for homogeneous noise (noise power constant across all cells). It suffers in heterogeneous environments and near closely-spaced targets.

pytcl.mathematical_functions.matched_filter(signal, template, normalize=True, mode='same')[source]

Apply matched filtering in the time domain.

The matched filter maximizes the output signal-to-noise ratio (SNR) when the input contains a known signal plus white Gaussian noise.

Parameters:

signal (array_like) – Input signal (may contain noise).
template (array_like) – Template signal to match (the known waveform).
normalize (bool, optional) – If True, normalize output by template energy. Default is True.
mode ({'full', 'same', 'valid'}, optional) – Convolution mode. Default is ‘same’.

Returns:

result – Named tuple with filter output, peak location, and SNR gain.

Return type:

MatchedFilterResult

Examples

>>> import numpy as np
>>> # Create a signal with a pulse
>>> template = np.array([1, 1, 1, 1, 1])
>>> signal = np.zeros(100)
>>> signal[50:55] = template
>>> result = matched_filter(signal, template)
>>> 50 <= result.peak_index <= 54
True

Notes

The matched filter is the time-reversed, conjugated version of the template convolved with the signal. For real signals, this is equivalent to cross-correlation.

The theoretical SNR gain of a matched filter is equal to the number of samples in the template (for unit-energy template in unit-variance noise).

pytcl.mathematical_functions.fft(x, n=None, axis=-1, norm=None)[source]

Compute the one-dimensional discrete Fourier Transform.

Parameters:

x (array_like) – Input array, can be complex.
n (int, optional) – Length of the transformed axis of the output. If n is smaller than the length of the input, the input is cropped. If it is larger, the input is padded with zeros.
axis (int, optional) – Axis over which to compute the FFT. Default is -1.
norm ({None, "ortho", "forward", "backward"}, optional) – Normalization mode. Default is None (backward).

Returns:

out – The truncated or zero-padded input, transformed along the axis.

Return type:

ndarray

Examples

>>> import numpy as np
>>> x = np.array([1.0, 2.0, 1.0, -1.0])
>>> X = fft(x)
>>> np.allclose(X, [3.+0.j, 0.-2.j, 1.+0.j, 0.+2.j])
True

pytcl.mathematical_functions.ifft(X, n=None, axis=-1, norm=None)[source]

Compute the one-dimensional inverse discrete Fourier Transform.

Parameters:

X (array_like) – Input array, can be complex.
n (int, optional) – Length of the transformed axis of the output.
axis (int, optional) – Axis over which to compute the inverse FFT. Default is -1.
norm ({None, "ortho", "forward", "backward"}, optional) – Normalization mode. Default is None (backward).

Returns:

out – The truncated or zero-padded input, transformed along the axis.

Return type:

ndarray

Examples

>>> import numpy as np
>>> X = np.array([3.+0.j, 0.-2.j, 1.+0.j, 0.+2.j])
>>> x = ifft(X)
>>> np.allclose(x, [1.0, 2.0, 1.0, -1.0])
True

pytcl.mathematical_functions.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend=False, return_onesided=True, boundary='zeros', padded=True)[source]

Compute the Short-Time Fourier Transform.

Parameters:

x (array_like) – Input time-domain signal.
fs (float, optional) – Sampling frequency in Hz. Default is 1.0.
window (str, tuple, or array_like, optional) – Window function. Default is ‘hann’.
nperseg (int, optional) – Length of each segment. Default is 256.
noverlap (int, optional) – Number of points to overlap between segments. Default is nperseg // 2.
nfft (int, optional) – Length of the FFT used. Default is nperseg.
detrend (str or bool, optional) – Detrending: ‘constant’, ‘linear’, or False. Default is False.
return_onesided (bool, optional) – If True, return only non-negative frequencies for real input. Default is True.
boundary (str or None, optional) – Boundary extension: ‘zeros’, ‘even’, ‘odd’, or None. Default is ‘zeros’.
padded (bool, optional) – Whether to pad the signal. Default is True.

Returns:

result – Named tuple with frequencies, times, and STFT matrix.

Return type:

STFTResult

Examples

>>> import numpy as np
>>> fs = 1000
>>> t = np.arange(0, 1, 1/fs)
>>> x = np.sin(2 * np.pi * 50 * t)  # 50 Hz sine
>>> result = stft(x, fs=fs, nperseg=128)
>>> result.Zxx.shape  # (n_freq, n_time)
(65, 16)

Notes

The STFT provides a time-frequency representation where: - Time resolution = nperseg / fs - Frequency resolution = fs / nfft

There is a trade-off between time and frequency resolution (uncertainty principle): better time resolution requires shorter segments, which reduces frequency resolution, and vice versa.

pytcl.mathematical_functions.spectrogram(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend='constant', scaling='density', mode='psd')[source]

Compute a spectrogram (power spectral density over time).

Parameters:

x (array_like) – Input time-domain signal.
fs (float, optional) – Sampling frequency in Hz. Default is 1.0.
window (str, tuple, or array_like, optional) – Window function. Default is ‘hann’.
nperseg (int, optional) – Length of each segment. Default is 256.
noverlap (int, optional) – Overlap between segments. Default is nperseg // 8.
nfft (int, optional) – FFT length. Default is nperseg.
detrend (str or bool, optional) – Detrending: ‘constant’, ‘linear’, or False. Default is ‘constant’.
scaling ({'density', 'spectrum'}, optional) – ‘density’ for PSD (V^2/Hz), ‘spectrum’ for power (V^2). Default is ‘density’.
mode ({'psd', 'complex', 'magnitude', 'angle', 'phase'}, optional) – Return type. Default is ‘psd’.

Returns:

result – Named tuple with frequencies, times, and power spectrogram.

Return type:

Spectrogram

Examples

>>> import numpy as np
>>> fs = 1000
>>> t = np.arange(0, 2, 1/fs)
>>> # Chirp from 50 to 200 Hz
>>> x = np.sin(2 * np.pi * (50 + 75*t) * t)
>>> result = spectrogram(x, fs=fs, nperseg=128)
>>> result.power.shape  # (n_freq, n_time)
(65, 31)

Notes

The spectrogram is computed by taking the magnitude squared of the STFT. It shows how the spectral content of the signal evolves over time.

pytcl.mathematical_functions.power_spectrum(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', scaling='density')[source]

Estimate power spectral density using Welch’s method.

Parameters:

x (array_like) – Time series data.
fs (float, optional) – Sampling frequency in Hz. Default is 1.0.
window (str, optional) – Window function to use. Default is ‘hann’.
nperseg (int, optional) – Length of each segment. Default is 256.
noverlap (int, optional) – Number of points to overlap between segments. Default is nperseg//2.
nfft (int, optional) – Length of the FFT used. Default is nperseg.
detrend (str or bool, optional) – Detrending method: ‘constant’, ‘linear’, or False. Default is ‘constant’.
scaling ({'density', 'spectrum'}, optional) – ‘density’ for power spectral density (V^2/Hz), ‘spectrum’ for power spectrum (V^2). Default is ‘density’.

Returns:

result – Named tuple with frequencies and power spectral density.

Return type:

PowerSpectrum

Examples

>>> import numpy as np
>>> fs = 1000  # 1 kHz
>>> t = np.arange(0, 1, 1/fs)
>>> x = np.sin(2 * np.pi * 100 * t)  # 100 Hz sine
>>> result = power_spectrum(x, fs=fs)
>>> peak_freq = result.frequencies[np.argmax(result.psd)]
>>> abs(peak_freq - 100) < 5  # Peak near 100 Hz
True

Notes

Uses Welch’s method which averages modified periodograms of overlapping segments to reduce variance of the spectral estimate.

pytcl.mathematical_functions.cwt(signal, scales, wavelet='morlet', fs=1.0, method='fft')[source]

Compute the Continuous Wavelet Transform.

Parameters:

signal (array_like) – Input signal.
scales (array_like) – Scale values to use.
wavelet (str or callable, optional) – Wavelet to use. Options: ‘morlet’, ‘ricker’, ‘gaussian1’, ‘gaussian2’, or a callable. Default is ‘morlet’.
fs (float, optional) – Sampling frequency in Hz. Default is 1.0.
method ({'fft', 'conv'}, optional) – Computation method. ‘fft’ is faster for long signals. Default is ‘fft’.

Returns:

result – Named tuple with coefficients, scales, and frequencies.

Return type:

CWTResult

Examples

>>> import numpy as np
>>> fs = 1000
>>> t = np.arange(0, 1, 1/fs)
>>> x = np.sin(2 * np.pi * 50 * t)
>>> scales = np.arange(1, 128)
>>> result = cwt(x, scales, wavelet='morlet', fs=fs)
>>> result.coefficients.shape
(127, 1000)

Notes

The CWT is computed as:: W(a, b) = integral s(t) * (1/sqrt(a)) * psi*((t-b)/a) dt

where a is the scale, b is the translation, and psi is the wavelet.

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)

pytcl.mathematical_functions.basic_matrix.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.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.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.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.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.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.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.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)

pytcl.mathematical_functions.basic_matrix.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.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.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.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.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.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.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.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.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.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.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.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.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.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.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 import duplication_matrix, vech, 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 = vech(A)
>>> # Duplication matrix should reconstruct full vectorization
>>> vec_A = D @ vech_A
>>> np.allclose(vec_A, vec(A))
True

pytcl.mathematical_functions.basic_matrix.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 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)
>>> np.allclose(vech_A, [1.0, 3.0, 4.0])
True

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 import duplication_matrix, vech, 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 = vech(A)
>>> # Duplication matrix should reconstruct full vectorization
>>> vec_A = D @ vech_A
>>> 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 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)
>>> 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)

pytcl.mathematical_functions.special_functions.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

>>> 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.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

>>> bessely(0, 1)
0.088...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> besseli(0, 0)
1.0

See also

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

pytcl.mathematical_functions.special_functions.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

>>> besselk(0, 1)
0.421...
>>> besselk(1, 2)
0.139...

See also

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

pytcl.mathematical_functions.special_functions.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)
>>> h.real
0.765...
>>> h.imag
0.088...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> spherical_jn(0, 1)  # sin(1)/1
0.841...
>>> spherical_jn(0, 1, derivative=True)  # Derivative
0.301...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> spherical_yn(0, 1)  # -cos(1)/1
-0.540...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> spherical_in(0, 1)  # sinh(1)/1
1.175...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> spherical_kn(0, 1)  # (pi/2) * exp(-1)
0.578...

See also

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

pytcl.mathematical_functions.special_functions.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)
>>> Ai
0.355...
>>> Bi
0.614...

See also

scipy.special.airy: Airy functions.

pytcl.mathematical_functions.special_functions.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

>>> bessel_ratio(0, 1)  # J_1(1) / J_0(1)
0.5767...

pytcl.mathematical_functions.special_functions.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

>>> bessel_deriv(0, 1, kind='j')  # -J_1(1)
-0.4400...

pytcl.mathematical_functions.special_functions.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.404..., 5.520..., 8.653...])

pytcl.mathematical_functions.special_functions.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

>>> struve_h(0, 1)
0.5688...

pytcl.mathematical_functions.special_functions.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

>>> struve_l(0, 1)
0.710...

pytcl.mathematical_functions.special_functions.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)
>>> ber
0.984...

pytcl.mathematical_functions.special_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

>>> gamma(5)  # 4! = 24
24.0
>>> gamma(0.5)  # sqrt(pi)
1.7724538509055159

See also

scipy.special.gamma: Gamma function.

pytcl.mathematical_functions.special_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

>>> gammaln(100)  # log(99!)
359.1342053695754

See also

scipy.special.gammaln: Log of gamma function.

pytcl.mathematical_functions.special_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

>>> gammainc(1, 1)  # 1 - exp(-1)
0.632...

See also

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

pytcl.mathematical_functions.special_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

>>> gammaincc(1, 1)  # exp(-1)
0.367...

See also

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

pytcl.mathematical_functions.special_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

>>> gammaincinv(1, 0.5)  # Median of exponential distribution
0.693...

See also

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

pytcl.mathematical_functions.special_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

>>> digamma(1)  # -γ (negative Euler-Mascheroni constant)
-0.577...

See also

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

pytcl.mathematical_functions.special_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

>>> polygamma(0, 1)  # Digamma at 1 = -γ
-0.577...
>>> polygamma(1, 1)  # Trigamma at 1 = π²/6
1.644...

See also

scipy.special.polygamma: Polygamma function.

pytcl.mathematical_functions.special_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

>>> beta(1, 1)
1.0
>>> beta(0.5, 0.5)  # pi
3.141592653589793

See also

scipy.special.beta: Beta function.

pytcl.mathematical_functions.special_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
>>> betaln(100, 100)  # More stable than log(beta(100, 100))
-137.74...

See also

scipy.special.betaln: Log of beta function.

pytcl.mathematical_functions.special_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

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

See also

scipy.special.betainc: Regularized incomplete beta function.

pytcl.mathematical_functions.special_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

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

See also

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

pytcl.mathematical_functions.special_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

>>> 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.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

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

See also

scipy.special.factorial2: Double factorial.

pytcl.mathematical_functions.special_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

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

See also

scipy.special.comb: Combinations.

pytcl.mathematical_functions.special_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

>>> perm(5, 2)
20.0

See also

scipy.special.perm: Permutations.

pytcl.mathematical_functions.special_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

>>> erf(0)
0.0
>>> erf(1)
0.8427007929497149

See also

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

pytcl.mathematical_functions.special_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

>>> erfc(0)
1.0
>>> erfc(3)  # Very small
2.2090496998585438e-05

See also

scipy.special.erfc: Complementary error function.

pytcl.mathematical_functions.special_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

>>> erfcx(0)
1.0
>>> erfcx(10)  # Remains finite even when erfc(10) underflows
0.056...

See also

scipy.special.erfcx: Scaled complementary error function.

pytcl.mathematical_functions.special_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

>>> erfi(0)
0.0
>>> erfi(1)
1.650...

See also

scipy.special.erfi: Imaginary error function.

pytcl.mathematical_functions.special_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

>>> erfinv(0)
0.0
>>> erf(erfinv(0.5))
0.5

See also

scipy.special.erfinv: Inverse error function.

pytcl.mathematical_functions.special_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

>>> erfcinv(1)
0.0
>>> erfc(erfcinv(0.5))
0.5

See also

scipy.special.erfcinv: Inverse complementary error function.

pytcl.mathematical_functions.special_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

>>> dawsn(0)
0.0
>>> dawsn(1)
0.538...

See also

scipy.special.dawsn: Dawson’s integral.

pytcl.mathematical_functions.special_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)
>>> S
0.438...
>>> C
0.779...

See also

scipy.special.fresnel: Fresnel integrals.

pytcl.mathematical_functions.special_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)
>>> w.real
1.0
>>> w.imag
0.0

See also

scipy.special.wofz: Faddeeva function.

pytcl.mathematical_functions.special_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

>>> voigt_profile(0, 1, 0)  # Pure Gaussian at x=0
0.398...

See also

scipy.special.voigt_profile: Voigt profile.

pytcl.mathematical_functions.special_functions.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

>>> ellipk(0)  # K(0) = π/2
1.5707963267948966
>>> ellipk(0.5)
1.8540746773013719

See also

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

pytcl.mathematical_functions.special_functions.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

>>> ellipkm1(0.1)  # K(0.9)
2.578...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> ellipe(0)  # E(0) = π/2
1.5707963267948966
>>> ellipe(1)  # E(1) = 1
1.0

See also

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

pytcl.mathematical_functions.special_functions.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
>>> ellipeinc(np.pi/2, 0)  # Same as ellipe(0) = π/2
1.5707...

See also

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

pytcl.mathematical_functions.special_functions.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
>>> ellipkinc(np.pi/2, 0)  # Same as ellipk(0) = π/2
1.5707...

See also

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

pytcl.mathematical_functions.special_functions.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

>>> elliprd(1, 2, 3)
0.297...

See also

scipy.special.elliprd: Carlson R_D integral.

pytcl.mathematical_functions.special_functions.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

>>> 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.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

>>> 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.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

>>> elliprj(1, 2, 3, 4)
0.213...

See also

scipy.special.elliprj: Carlson R_J integral.

pytcl.mathematical_functions.special_functions.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) = arctan(sqrt((x-y)/y)) / sqrt(x-y) for x > y
R_C(x, y) = arctanh(sqrt((y-x)/y)) / sqrt(y-x) for x < y

Examples

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

See also

scipy.special.elliprc: Carlson R_C integral.

pytcl.mathematical_functions.special_functions.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

>>> marcum_q(0, 0)  # Q_1(0, 0) = 1
1.0
>>> marcum_q(3, 4)  # Standard Marcum Q
0.17789...

References

pytcl.mathematical_functions.special_functions.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

>>> marcum_q1(2, 2)
0.735...

See also

marcum_q: Generalized Marcum Q function.

pytcl.mathematical_functions.special_functions.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

>>> log_marcum_q(1, 5)  # log(Q_1(1, 5))
-10.96...

pytcl.mathematical_functions.special_functions.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
>>> marcum_q(3, b)  # Verify
0.5

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

Nuttall Q function (complementary Marcum Q).

Computes 1 - Q_1(a, b), which is the CDF of the Rician distribution.

Parameters:

a (array_like) – First argument (non-centrality parameter), a >= 0.
b (array_like) – Second argument (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).

Examples

>>> nuttall_q(2, 2)  # 1 - Q_1(2, 2)
0.264...

See also

marcum_q: Marcum Q function.

pytcl.mathematical_functions.special_functions.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

For Swerling 0 (non-fluctuating):: P_d = Q_n(sqrt(2*n*SNR), sqrt(threshold))
For Swerling 1:: P_d = exp(-threshold / (2 + 2*n*SNR)) * (1 + n*SNR/…)

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

pytcl.mathematical_functions.special_functions.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

>>> lambert_w(0)  # W(0) = 0
0.0
>>> lambert_w(np.e)  # W(e) = 1
1.0
>>> lambert_w(-np.exp(-1))  # W(-1/e) = -1
-1.0

References

pytcl.mathematical_functions.special_functions.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

>>> lambert_w_real(1)
0.5671432904097838
>>> lambert_w_real(-0.2, branch=-1)
-2.5426413577735264

pytcl.mathematical_functions.special_functions.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()
>>> omega * np.exp(omega)  # Should equal 1
1.0

pytcl.mathematical_functions.special_functions.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

>>> wright_omega(0)  # Omega constant
(0.5671...+0j)

References

pytcl.mathematical_functions.special_functions.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
>>> x  # Should be 1
1.0

pytcl.mathematical_functions.special_functions.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
>>> s + np.exp(-s)  # Should be approximately 0
~0

pytcl.mathematical_functions.special_functions.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

>>> debye(3, 0)  # D_3(0) = 1
1.0
>>> debye(3, 1)
0.674...
>>> debye(3, 10)
0.0192...

References

pytcl.mathematical_functions.special_functions.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_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_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_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_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) = 3 * (T/Θ_D)^3 * D_3(Θ_D/T)

Limits: - High T (T >> Θ_D): C_V -> 3*N*k_B (classical) - Low T (T << Θ_D): C_V ~ (T/Θ_D)^3 (quantum)

Examples

>>> # Aluminum at room temperature (Θ_D ≈ 428 K)
>>> cv = debye_heat_capacity(300, 428)  # ~0.95

pytcl.mathematical_functions.special_functions.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*D_3(Θ_D/T) - 3*ln(1 - exp(-Θ_D/T))

pytcl.mathematical_functions.special_functions.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

>>> hyp0f1(1, 0)  # 0F1(1; 0) = 1
1.0
>>> hyp0f1(1, 1)
2.279585...

References

pytcl.mathematical_functions.special_functions.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

>>> hyp1f1(1, 1, 1)  # exp(1)
2.718281828...
>>> hyp1f1(0.5, 1.5, -1)  # Related to erf
0.842700...

References

pytcl.mathematical_functions.special_functions.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

>>> hyp2f1(1, 1, 2, 0.5)  # -log(1-0.5)/0.5 = log(2)
1.386294...
>>> hyp2f1(0.5, 0.5, 1.5, 0.25)  # Related to arcsin
1.072379...

References

pytcl.mathematical_functions.special_functions.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

>>> hyperu(1, 1, 1)
0.596347...

pytcl.mathematical_functions.special_functions.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.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

>>> pochhammer(1, 5)  # 5!
120.0
>>> pochhammer(3, 4)  # 3*4*5*6
360.0

pytcl.mathematical_functions.special_functions.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

>>> falling_factorial(5, 3)  # 5*4*3
60.0

pytcl.mathematical_functions.special_functions.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

>>> generalized_hypergeometric([1], [2], 1)  # 1F1(1; 2; 1) ~ 1.718...
1.718281...

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

>>> 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

>>> bessely(0, 1)
0.088...

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

>>> 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

>>> besselk(0, 1)
0.421...
>>> besselk(1, 2)
0.139...

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)
>>> h.real
0.765...
>>> h.imag
0.088...

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

>>> spherical_jn(0, 1)  # sin(1)/1
0.841...
>>> spherical_jn(0, 1, derivative=True)  # Derivative
0.301...

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

>>> spherical_yn(0, 1)  # -cos(1)/1
-0.540...

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

>>> spherical_in(0, 1)  # sinh(1)/1
1.175...

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

>>> spherical_kn(0, 1)  # (pi/2) * exp(-1)
0.578...

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)
>>> Ai
0.355...
>>> Bi
0.614...

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

>>> bessel_ratio(0, 1)  # J_1(1) / J_0(1)
0.5767...

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

>>> bessel_deriv(0, 1, kind='j')  # -J_1(1)
-0.4400...

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

>>> struve_h(0, 1)
0.5688...

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

>>> struve_l(0, 1)
0.710...

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.404..., 5.520..., 8.653...])

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)
>>> ber
0.984...

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

>>> gamma(5)  # 4! = 24
24.0
>>> gamma(0.5)  # sqrt(pi)
1.7724538509055159

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

>>> gammaln(100)  # log(99!)
359.1342053695754

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

>>> gammainc(1, 1)  # 1 - exp(-1)
0.632...

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

>>> gammaincc(1, 1)  # exp(-1)
0.367...

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

>>> gammaincinv(1, 0.5)  # Median of exponential distribution
0.693...

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

>>> digamma(1)  # -γ (negative Euler-Mascheroni constant)
-0.577...

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

>>> polygamma(0, 1)  # Digamma at 1 = -γ
-0.577...
>>> polygamma(1, 1)  # Trigamma at 1 = π²/6
1.644...

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

>>> beta(1, 1)
1.0
>>> beta(0.5, 0.5)  # pi
3.141592653589793

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
>>> betaln(100, 100)  # More stable than log(beta(100, 100))
-137.74...

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

>>> betainc(1, 1, 0.5)  # 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

>>> betaincinv(1, 1, 0.5)  # 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

>>> 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

>>> factorial2(5)  # 5 * 3 * 1 = 15
15.0
>>> factorial2(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

>>> comb(5, 2)
10.0
>>> 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

>>> 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

>>> erf(0)
0.0
>>> erf(1)
0.8427007929497149

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

>>> erfc(0)
1.0
>>> erfc(3)  # Very small
2.2090496998585438e-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

>>> erfcx(0)
1.0
>>> erfcx(10)  # Remains finite even when erfc(10) underflows
0.056...

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

>>> erfi(0)
0.0
>>> erfi(1)
1.650...

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

>>> erfinv(0)
0.0
>>> erf(erfinv(0.5))
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

>>> erfcinv(1)
0.0
>>> erfc(erfcinv(0.5))
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

>>> dawsn(0)
0.0
>>> dawsn(1)
0.538...

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)
>>> S
0.438...
>>> C
0.779...

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)
>>> w.real
1.0
>>> 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

>>> voigt_profile(0, 1, 0)  # Pure Gaussian at x=0
0.398...

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

>>> ellipk(0)  # K(0) = π/2
1.5707963267948966
>>> ellipk(0.5)
1.8540746773013719

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

>>> ellipkm1(0.1)  # K(0.9)
2.578...

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

>>> ellipe(0)  # E(0) = π/2
1.5707963267948966
>>> 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
>>> ellipeinc(np.pi/2, 0)  # Same as ellipe(0) = π/2
1.5707...

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
>>> ellipkinc(np.pi/2, 0)  # Same as ellipk(0) = π/2
1.5707...

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

>>> elliprd(1, 2, 3)
0.297...

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

>>> 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

>>> 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

>>> elliprj(1, 2, 3, 4)
0.213...

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) = arctan(sqrt((x-y)/y)) / sqrt(x-y) for x > y
R_C(x, y) = arctanh(sqrt((y-x)/y)) / sqrt(y-x) for x < y

Examples

>>> 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)

class pytcl.mathematical_functions.statistics.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.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)
>>> g.pdf(0)
0.3989422804014327
>>> g.cdf(0)
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.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]])
>>> mg.pdf([0, 0])
0.15915494309189535

__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.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.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.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.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.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.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.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.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.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.

pytcl.mathematical_functions.statistics.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.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]
>>> weighted_var(x, weights)
0.5625

pytcl.mathematical_functions.statistics.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.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.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.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.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.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.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.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

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

pytcl.mathematical_functions.statistics.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

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

pytcl.mathematical_functions.statistics.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

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

pytcl.mathematical_functions.statistics.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

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

pytcl.mathematical_functions.statistics.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.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.

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)
>>> g.pdf(0)
0.3989422804014327
>>> g.cdf(0)
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]])
>>> mg.pdf([0, 0])
0.15915494309189535

__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]
>>> weighted_var(x, weights)
0.5625

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

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

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

>>> 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

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

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

>>> moment([1, 2, 3, 4, 5], order=2)
2.0
>>> 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

>>> 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)
>>> cs(np.pi/2)
array(0.99999...)

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 with function values and derivatives
>>> x = np.array([0, 1, 2])
>>> # For Hermite interpolation, y rows are: function values, then derivatives
>>> y = np.array([
...     [1, 2, 3],      # Function values at x=0, 1, 2
...     [0, 1, 2],      # Derivatives at x=0, 1, 2
... ])
>>> k = krogh(x, y)
>>> # Evaluate interpolant at new points
>>> x_new = np.array([0.5, 1.5])
>>> y_interp = k(x_new)
>>> # Interpolant passes through original points and matches derivatives
>>> np.allclose(k(x), y[0])
True

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)

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)
>>> np.sum(w * x**2)
0.6666666666666666

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)
>>> np.sum(w * x)
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:

Examples

>>> result, error = quad(lambda x: x**2, 0, 1)
>>> result
0.33333333333333337

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:

Examples

>>> # Integrate x*y over unit square
>>> result, error = dblquad(lambda y, x: x*y, 0, 1, lambda x: 0, lambda x: 1)
>>> result  # 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:

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)
>>> result
0.3333333333333333

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:

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:

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:

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,)
>>> np.sum(weights)  # 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:

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:

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:

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)[source]

Compute minimum enclosing circle (2D).

Parameters:

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

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

Uses Welzl’s algorithm with expected O(n) time complexity.

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:

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:

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]
>>> 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:

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: