Mathematical Utilities

class KernelInterpolator:
    """
    Kernel-based interpolation of a 1-D function using convolution operations.
    
    Parameters:
    - x: independent variable values
    - y: dependent variable values to interpolate
    - window: width of the window in samples on each side (default: 3)
    - kernel: kernel function to use (default: kernel_lanczos)
    - kernel_kwargs: arguments for the kernel function
    - padding_kwargs: arguments for numpy.pad when extending data
    - dtype: data type for internal conversions
    """
    def __init__(self, x: ArrayLike, y: ArrayLike, window: float = 3, kernel: Callable = None, kernel_kwargs: dict = None, padding_kwargs: dict = None, dtype: Type = None): ...
    
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
    
    @property
    def x(self) -> NDArrayFloat: ...
    @property 
    def y(self) -> NDArrayFloat: ...
    @property
    def window(self) -> float: ...
    @property
    def kernel(self) -> Callable: ...
    @property
    def kernel_kwargs(self) -> dict: ...
    @property
    def padding_kwargs(self) -> dict: ...

class LinearInterpolator:
    """
    Linear interpolation of a 1-D function using numpy.interp.
    
    Parameters:
    - x: independent variable values
    - y: dependent variable values to interpolate
    - dtype: data type for internal conversions
    """
    def __init__(self, x: ArrayLike, y: ArrayLike, dtype: Type = None): ...
    
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
    
    @property
    def x(self) -> NDArrayFloat: ...
    @property
    def y(self) -> NDArrayFloat: ...

class SpragueInterpolator:
    """
    Fifth-order polynomial interpolation using Sprague (1880) method.
    Recommended by CIE for uniformly spaced independent variables.
    
    Parameters:
    - x: independent variable values
    - y: dependent variable values to interpolate (minimum 6 points required)
    - dtype: data type for internal conversions
    """
    def __init__(self, x: ArrayLike, y: ArrayLike, dtype: Type = None): ...
    
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
    
    @property
    def x(self) -> NDArrayFloat: ...
    @property
    def y(self) -> NDArrayFloat: ...

class CubicSplineInterpolator:
    """
    Cubic spline interpolation using scipy.interpolate.interp1d.
    
    Wrapper around scipy's cubic spline interpolation.
    """
    def __init__(self, *args, **kwargs): ...

class PchipInterpolator:
    """
    Piecewise Cubic Hermite Interpolating Polynomial interpolation.
    
    Wrapper around scipy.interpolate.PchipInterpolator.
    """
    def __init__(self, x: ArrayLike, y: ArrayLike, *args, **kwargs): ...
    
    @property
    def y(self) -> NDArrayFloat: ...

class NearestNeighbourInterpolator(KernelInterpolator):
    """
    Nearest-neighbour interpolation using kernel-based approach.
    
    Specialized KernelInterpolator using nearest-neighbour kernel.
    """
    def __init__(self, *args, **kwargs): ...

class NullInterpolator:
    """
    Null interpolation returning existing values within tolerances.
    
    Parameters:
    - x: independent variable values
    - y: dependent variable values
    - absolute_tolerance: absolute tolerance for matching
    - relative_tolerance: relative tolerance for matching
    - default: default value for points outside tolerances
    - dtype: data type for internal conversions
    """
    def __init__(self, x: ArrayLike, y: ArrayLike, absolute_tolerance: float = 1e-10, relative_tolerance: float = 1e-5, default: float = np.nan, dtype: Type = None): ...
    
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
    
    @property
    def x(self) -> NDArrayFloat: ...
    @property
    def y(self) -> NDArrayFloat: ...
    @property
    def absolute_tolerance(self) -> float: ...
    @property
    def relative_tolerance(self) -> float: ...
    @property
    def default(self) -> float: ...

def kernel_nearest_neighbour(x: ArrayLike) -> NDArrayFloat:
    """
    Nearest-neighbour kernel returning 1 for |x| < 0.5, 0 otherwise.
    
    Parameters:
    - x: samples at which to evaluate the kernel
    
    Returns:
    Kernel values evaluated at given samples
    """

def kernel_linear(x: ArrayLike) -> NDArrayFloat:
    """
    Linear kernel returning (1 - |x|) for |x| < 1, 0 otherwise.
    
    Parameters:
    - x: samples at which to evaluate the kernel
    
    Returns:
    Kernel values evaluated at given samples
    """

def kernel_sinc(x: ArrayLike, a: float = 3) -> NDArrayFloat:
    """
    Sinc kernel for high-quality interpolation.
    
    Parameters:
    - x: samples at which to evaluate the kernel
    - a: kernel size parameter (must be >= 1)
    
    Returns:
    Kernel values evaluated at given samples
    """

def kernel_lanczos(x: ArrayLike, a: float = 3) -> NDArrayFloat:
    """
    Lanczos kernel for high-quality resampling.
    
    Parameters:
    - x: samples at which to evaluate the kernel
    - a: kernel size parameter (must be >= 1)
    
    Returns:
    Kernel values evaluated at given samples
    """

def kernel_cardinal_spline(x: ArrayLike, a: float = 0.5, b: float = 0.0) -> NDArrayFloat:
    """
    Cardinal spline kernel with configurable parameters.
    
    Notable parameterizations:
    - Catmull-Rom: (a=0.5, b=0)
    - Cubic B-Spline: (a=0, b=1) 
    - Mitchell-Netravalli: (a=1/3, b=1/3)
    
    Parameters:
    - x: samples at which to evaluate the kernel
    - a: control parameter a
    - b: control parameter b
    
    Returns:
    Kernel values evaluated at given samples
    """

def lagrange_coefficients(r: float, n: int = 4) -> NDArrayFloat:
    """
    Compute Lagrange coefficients for interpolation.
    
    Parameters:
    - r: point to get coefficients at
    - n: degree of the Lagrange coefficients
    
    Returns:
    Lagrange coefficients array
    """

def table_interpolation_trilinear(V_xyz: ArrayLike, table: ArrayLike) -> NDArrayFloat:
    """
    Trilinear interpolation using 4-dimensional lookup table.
    
    Parameters:
    - V_xyz: values to interpolate
    - table: 4-dimensional interpolation table (NxNxNx3)
    
    Returns:
    Interpolated values
    """

def table_interpolation_tetrahedral(V_xyz: ArrayLike, table: ArrayLike) -> NDArrayFloat:
    """
    Tetrahedral interpolation using 4-dimensional lookup table.
    
    Parameters:
    - V_xyz: values to interpolate  
    - table: 4-dimensional interpolation table (NxNxNx3)
    
    Returns:
    Interpolated values
    """

def table_interpolation(V_xyz: ArrayLike, table: ArrayLike, method: str = "Trilinear") -> NDArrayFloat:
    """
    Generic table interpolation with selectable method.
    
    Parameters:
    - V_xyz: values to interpolate
    - table: 4-dimensional interpolation table (NxNxNx3) 
    - method: interpolation method ("Trilinear" or "Tetrahedral")
    
    Returns:
    Interpolated values
    """

# Available table interpolation methods
TABLE_INTERPOLATION_METHODS = {
    "Trilinear": table_interpolation_trilinear,
    "Tetrahedral": table_interpolation_tetrahedral,
}

class Extrapolator:
    """
    Extrapolate 1-D functions beyond their interpolation range.
    
    Supports linear and constant extrapolation methods with customizable boundary values.
    
    Parameters:
    - interpolator: interpolator object with x and y properties
    - method: extrapolation method ("Linear" or "Constant")
    - left: value to return for x < xi[0] (overrides method)
    - right: value to return for x > xi[-1] (overrides method)
    - dtype: data type for internal conversions
    """
    def __init__(self, interpolator: ProtocolInterpolator = None, method: str = "Linear", left: Real = None, right: Real = None, dtype: Type = None): ...
    
    def __call__(self, x: ArrayLike) -> NDArrayFloat: ...
    
    @property
    def interpolator(self) -> ProtocolInterpolator: ...
    @property
    def method(self) -> str: ...
    @property
    def left(self) -> Real: ...
    @property
    def right(self) -> Real: ...

def sdiv(a: ArrayLike, b: ArrayLike) -> NDArrayFloat:
    """
    Safe division handling zero-division gracefully with configurable modes.
    
    Available modes (configured via sdiv_mode):
    - "Numpy": standard numpy zero-division handling
    - "Ignore": zero-division occurs silently
    - "Warning": zero-division with warnings
    - "Raise": raise exception on zero-division
    - "Ignore Zero Conversion": silent with NaN/inf converted to zero
    - "Warning Zero Conversion": warning with NaN/inf converted to zero
    - "Ignore Limit Conversion": silent with NaN/inf converted to limits
    - "Warning Limit Conversion": warning with NaN/inf converted to limits
    
    Parameters:
    - a: numerator array
    - b: denominator array
    
    Returns:
    Safely divided result
    """

def get_sdiv_mode() -> str:
    """Get current safe division mode."""

def set_sdiv_mode(mode: str) -> None:
    """Set safe division mode."""

class sdiv_mode:
    """
    Context manager for temporarily setting safe division mode.
    
    Parameters:
    - mode: temporary safe division mode
    """
    def __init__(self, mode: str = None): ...
    def __enter__(self) -> "sdiv_mode": ...
    def __exit__(self, *args): ...

def spow(a: ArrayLike, p: ArrayLike) -> NDArrayFloat:
    """
    Safe power function: sign(a) * |a|^p
    
    Avoids NaN generation when base is negative and exponent is fractional.
    
    Parameters:
    - a: base array
    - p: power/exponent array
    
    Returns:
    Safely computed power result
    """

def is_spow_enabled() -> bool:
    """Check if safe power function is enabled."""

def set_spow_enable(enable: bool) -> None:
    """Enable or disable safe power function."""

class spow_enable:
    """
    Context manager for temporarily setting safe power function state.
    
    Parameters:
    - enable: whether to enable safe power function
    """
    def __init__(self, enable: bool): ...
    def __enter__(self) -> "spow_enable": ...
    def __exit__(self, *args): ...

def normalise_vector(a: ArrayLike) -> NDArrayFloat:
    """
    Normalize vector to unit length using L2 norm.
    
    Parameters:
    - a: vector to normalize
    
    Returns:
    Normalized vector
    """

def normalise_maximum(a: ArrayLike, axis: int = None, factor: float = 1, clip: bool = True) -> NDArrayFloat:
    """
    Normalize array values by maximum value with optional clipping.
    
    Parameters:
    - a: array to normalize
    - axis: normalization axis
    - factor: normalization factor
    - clip: whether to clip values to domain [0, factor]
    
    Returns:
    Maximum-normalized array
    """

def vecmul(m: ArrayLike, v: ArrayLike) -> NDArrayFloat:
    """
    Batched matrix-vector multiplication with broadcasting.
    
    Equivalent to np.einsum('...ij,...j->...i', m, v).
    
    Parameters:
    - m: matrix array 
    - v: vector array
    
    Returns:
    Matrix-vector multiplication result
    """

def is_identity(a: ArrayLike) -> bool:
    """
    Test whether array is an identity matrix.
    
    Parameters:
    - a: array to test
    
    Returns:
    True if array is identity matrix
    """

def eigen_decomposition(a: ArrayLike, eigen_w_v_count: int = None, descending_order: bool = True, covariance_matrix: bool = False) -> tuple:
    """
    Compute eigenvalues and eigenvectors with optional ordering.
    
    Parameters:
    - a: array to decompose
    - eigen_w_v_count: number of eigenvalues/eigenvectors to return
    - descending_order: whether to return in descending eigenvalue order
    - covariance_matrix: whether to use covariance matrix A^T * A
    
    Returns:
    Tuple of (eigenvalues, eigenvectors)
    """

def euclidean_distance(a: ArrayLike, b: ArrayLike) -> NDArrayFloat:
    """
    Euclidean (L2) distance between point arrays.
    
    For 2D: sqrt((xa - xb)² + (ya - yb)²)
    
    Parameters:
    - a: first point array
    - b: second point array
    
    Returns:
    Euclidean distance
    """

def manhattan_distance(a: ArrayLike, b: ArrayLike) -> NDArrayFloat:
    """
    Manhattan (L1, City-Block) distance between point arrays.
    
    For 2D: |xa - xb| + |ya - yb|
    
    Parameters:
    - a: first point array
    - b: second point array
    
    Returns:
    Manhattan distance
    """

def linear_conversion(a: ArrayLike, old_range: ArrayLike, new_range: ArrayLike) -> NDArrayFloat:
    """
    Linear conversion between value ranges.
    
    Parameters:
    - a: array to convert
    - old_range: [old_min, old_max] current range
    - new_range: [new_min, new_max] target range
    
    Returns:
    Linearly converted values
    """

def linstep_function(x: ArrayLike, a: ArrayLike = 0, b: ArrayLike = 1, clip: bool = False) -> NDArrayFloat:
    """
    Linear interpolation between two values: (1-x)*a + x*b
    
    Parameters:
    - x: interpolation parameter
    - a: start value  
    - b: end value
    - clip: whether to clip result to [a, b]
    
    Returns:
    Linearly interpolated result
    """

# Alias for linstep_function
lerp = linstep_function

def smoothstep_function(x: ArrayLike, a: ArrayLike = 0, b: ArrayLike = 1, clip: bool = False) -> NDArrayFloat:
    """
    Smooth sigmoid-like interpolation function.
    
    Uses polynomial: t² * (3 - 2t) where t = (x-a)/(b-a)
    
    Parameters:
    - x: input array
    - a: low input domain limit (left edge)
    - b: high input domain limit (right edge)
    - clip: whether to scale and clip input to [a, b]
    
    Returns:
    Smoothly interpolated result
    """

# Alias for smoothstep_function  
smooth = smoothstep_function

def cartesian_to_spherical(a: ArrayLike) -> NDArrayFloat:
    """
    Transform Cartesian coordinates (x,y,z) to spherical (ρ,θ,φ).
    
    Parameters:
    - a: Cartesian coordinates [x, y, z]
    
    Returns:
    Spherical coordinates [ρ, θ, φ] where:
    - ρ: radial distance [0, +∞)
    - θ: inclination [0, π] radians
    - φ: azimuth [-π, π] radians
    """

def spherical_to_cartesian(a: ArrayLike) -> NDArrayFloat:
    """
    Transform spherical coordinates (ρ,θ,φ) to Cartesian (x,y,z).
    
    Parameters:
    - a: Spherical coordinates [ρ, θ, φ]
    
    Returns:
    Cartesian coordinates [x, y, z]
    """

def cartesian_to_polar(a: ArrayLike) -> NDArrayFloat:
    """
    Transform Cartesian coordinates (x,y) to polar (ρ,φ).
    
    Parameters:
    - a: Cartesian coordinates [x, y]
    
    Returns:
    Polar coordinates [ρ, φ] where:
    - ρ: radial coordinate [0, +∞)
    - φ: angular coordinate [-π, π] radians
    """

def polar_to_cartesian(a: ArrayLike) -> NDArrayFloat:
    """
    Transform polar coordinates (ρ,φ) to Cartesian (x,y).
    
    Parameters:
    - a: Polar coordinates [ρ, φ]
    
    Returns:
    Cartesian coordinates [x, y]
    """

def cartesian_to_cylindrical(a: ArrayLike) -> NDArrayFloat:
    """
    Transform Cartesian coordinates (x,y,z) to cylindrical (ρ,φ,z).
    
    Parameters:
    - a: Cartesian coordinates [x, y, z]
    
    Returns:
    Cylindrical coordinates [ρ, φ, z] where:
    - ρ: radial distance [0, +∞) 
    - φ: azimuth [-π, π] radians
    - z: height [0, +∞)
    """

def cylindrical_to_cartesian(a: ArrayLike) -> NDArrayFloat:
    """
    Transform cylindrical coordinates (ρ,φ,z) to Cartesian (x,y,z).
    
    Parameters:
    - a: Cylindrical coordinates [ρ, φ, z]
    
    Returns:
    Cartesian coordinates [x, y, z]
    """

def least_square_mapping_MoorePenrose(y: ArrayLike, x: ArrayLike) -> NDArrayFloat:
    """
    Least-squares mapping using Moore-Penrose inverse.
    
    Computes optimal linear mapping from y to x variables.
    
    Parameters:
    - y: dependent variable (known values)
    - x: independent variable (target values)
    
    Returns:
    Least-squares mapping matrix
    """

def random_triplet_generator(size: int, limits: ArrayLike = [[0,1], [0,1], [0,1]], random_state: np.random.RandomState = None) -> NDArrayFloat:
    """
    Generate random triplets within specified limits.
    
    Parameters:
    - size: number of triplets to generate
    - limits: [[x_min, x_max], [y_min, y_max], [z_min, z_max]] limits for each axis
    - random_state: random number generator state for reproducibility
    
    Returns:
    Array of random triplets shape (size, 3)
    """

# Global random state for reproducible results
RANDOM_STATE = np.random.RandomState()

# CIE Constants
CONSTANT_K_M: float  # CIE luminous efficacy constant (683 lm/W)
CONSTANT_KP_M: float  # CIE scotopic luminous efficacy constant (1700 lm/W)

# CODATA Physical Constants  
CONSTANT_AVOGADRO: float      # Avogadro constant (6.02214076e23 mol⁻¹)
CONSTANT_BOLTZMANN: float     # Boltzmann constant (1.380649e-23 J/K)
CONSTANT_LIGHT_SPEED: float   # Speed of light in vacuum (299792458 m/s)
CONSTANT_PLANCK: float        # Planck constant (6.62607015e-34 J⋅s)

# Numerical Constants
EPSILON: float                    # Machine epsilon for floating point comparisons
DTYPE_INT_DEFAULT: type          # Default integer data type (np.int64)  
DTYPE_FLOAT_DEFAULT: type        # Default float data type (np.float64)
TOLERANCE_ABSOLUTE_DEFAULT: float # Default absolute tolerance (1e-10)
TOLERANCE_RELATIVE_DEFAULT: float # Default relative tolerance (1e-5)
TOLERANCE_ABSOLUTE_TESTS: float  # Absolute tolerance for tests (1e-7)
TOLERANCE_RELATIVE_TESTS: float  # Relative tolerance for tests (1e-7)

# Utility Constants
THRESHOLD_INTEGER: float         # Threshold for integer detection
PATTERN_FLOATING_POINT_NUMBER: str # Regex pattern for floating point numbers

import numpy as np
from colour.algebra import LinearInterpolator, SpragueInterpolator

# Linear interpolation
x = np.array([1, 2, 3, 4, 5])
y = np.array([1, 4, 9, 16, 25])  # x²
interp = LinearInterpolator(x, y)

# Interpolate at intermediate points
result = interp([1.5, 2.5, 3.5])
print(result)  # [2.5, 6.5, 12.5]

# High-quality Sprague interpolation
sprague = SpragueInterpolator(x, y)
result = sprague([1.5, 2.5, 3.5])
print(result)  # More accurate cubic interpolation

from colour.algebra import sdiv, spow, sdiv_mode

# Safe division handling zero denominators
a = np.array([1, 2, 3])
b = np.array([2, 0, 4])

# Default mode converts inf/nan to zero
result = sdiv(a, b)
print(result)  # [0.5, 0.0, 0.75]

# Change division behavior
with sdiv_mode("Warning"):
    result = sdiv(a, b)  # Will warn about zero division

# Safe power for negative bases
result = spow(-2, 0.5)  # Returns -√2 instead of NaN
print(result)  # -1.4142135...

from colour.algebra import normalise_vector, euclidean_distance, vecmul

# Normalize vectors
vector = np.array([3, 4, 5])
normalized = normalise_vector(vector)
print(normalized)  # [0.424, 0.566, 0.707] (unit length)

# Calculate distances
point_a = np.array([1, 2, 3])
point_b = np.array([4, 6, 8])
distance = euclidean_distance(point_a, point_b)
print(distance)  # 7.071...

# Matrix-vector multiplication with broadcasting
matrix = np.random.random((100, 3, 3))  # 100 3x3 matrices
vectors = np.random.random((100, 3))     # 100 3D vectors
result = vecmul(matrix, vectors)         # 100 transformed vectors
print(result.shape)  # (100, 3)

from colour.algebra import (
    cartesian_to_spherical, spherical_to_cartesian,
    cartesian_to_polar, polar_to_cartesian
)

# 3D coordinate transformations
cartesian = np.array([3, 1, 6])
spherical = cartesian_to_spherical(cartesian)
print(spherical)  # [6.782, 0.485, 0.322] (ρ, θ, φ)

# Convert back
back_to_cartesian = spherical_to_cartesian(spherical)
print(back_to_cartesian)  # [3, 1, 6] (approximately)

# 2D coordinate transformations  
cartesian_2d = np.array([3, 4])
polar = cartesian_to_polar(cartesian_2d)
print(polar)  # [5.0, 0.927] (ρ, φ)

from colour.algebra import LinearInterpolator, Extrapolator

# Create base interpolator
x = np.array([2, 3, 4])
y = np.array([4, 9, 16])
interp = LinearInterpolator(x, y)

# Linear extrapolation
extrap = Extrapolator(interp, method="Linear")
result = extrap([0, 1, 5, 6])  # Extrapolate beyond [2, 4] range
print(result)  # [-11, -4, 25, 32] (linear extrapolation)

# Constant extrapolation with custom boundaries
extrap_const = Extrapolator(interp, method="Constant", left=0, right=100)
result = extrap_const([0, 1, 5, 6])
print(result)  # [0, 0, 100, 100] (constant boundaries)

from colour.algebra import table_interpolation
import colour

# Load a 3D LUT for color transformation
lut_path = "path/to/color_lut.cube"
lut = colour.read_LUT(lut_path)

# Interpolate RGB values using the LUT
rgb_values = np.array([
    [0.5, 0.3, 0.8],
    [0.2, 0.7, 0.1],
    [0.9, 0.1, 0.4]
])

# Trilinear interpolation (default)
result_trilinear = table_interpolation(rgb_values, lut.table, method="Trilinear")

# Tetrahedral interpolation (more accurate)
result_tetrahedral = table_interpolation(rgb_values, lut.table, method="Tetrahedral")

print("Trilinear:", result_trilinear)
print("Tetrahedral:", result_tetrahedral)