diff --git a/src/xspline/bspl.py b/src/xspline/bspl.py index 8b16156..617a971 100644 --- a/src/xspline/bspl.py +++ b/src/xspline/bspl.py @@ -7,12 +7,24 @@ def cache_bspl(function: RawFunction) -> RawFunction: + """Cache implementation for bspline basis functions, to avoid repetitively + evaluate functions. + + Parameters + ---------- + function + Raw value, derivative and definite integral functions. + + Returns + ------- + describe + Cached version of the raw functions. + + """ cache = {} def wrapper_function(*args, **kwargs) -> NDArray: - key = tuple( - tuple(x.ravel()) if isinstance(x, np.ndarray) else x for x in args - ) + key = tuple(tuple(x.ravel()) if isinstance(x, np.ndarray) else x for x in args) if key in cache: return cache[key] result = function(*args, **kwargs) @@ -28,6 +40,22 @@ def cache_clear(): @cache_bspl def bspl_val(params: BsplParams, x: NDArray) -> NDArray: + """Value of the bspline function. + + Parameters + ---------- + params + Bspline function parameters as a tuple including, knots, degree and the + index of the spline basis. + x + Data points. + + Returns + ------- + describe + Function value of the bspline function. + + """ # knots, degree, and index t, k, i = params @@ -44,10 +72,10 @@ def bspl_val(params: BsplParams, x: NDArray) -> NDArray: if t[ii[0]] != t[ii[2]]: n0 = bspl_val((t, k - 1, i), x) - val0 = (x - t[ii[0]])*n0/(t[ii[2]] - t[ii[0]]) + val0 = (x - t[ii[0]]) * n0 / (t[ii[2]] - t[ii[0]]) if t[ii[1]] != t[ii[3]]: n1 = bspl_val((t, k - 1, i + 1), x) - val1 = (t[ii[3]] - x)*n1/(t[ii[3]] - t[ii[1]]) + val1 = (t[ii[3]] - x) * n1 / (t[ii[3]] - t[ii[1]]) val = val0 + val1 @@ -56,6 +84,22 @@ def bspl_val(params: BsplParams, x: NDArray) -> NDArray: @cache_bspl def bspl_der(params: BsplParams, x: NDArray, order: int) -> NDArray: + """Derivative of the bspline function. + + Parameters + ---------- + params + Bspline function parameters as a tuple including, knots, degree and the + index of the spline basis. + x + Data points. + + Returns + ------- + describe + Derivative of the bspline function. + + """ # knots, degree, and index t, k, i = params @@ -72,10 +116,10 @@ def bspl_der(params: BsplParams, x: NDArray, order: int) -> NDArray: if t[ii[0]] != t[ii[2]]: n0 = bspl_der((t, k - 1, i), x, order - 1) - val0 = k*n0/(t[ii[2]] - t[ii[0]]) + val0 = k * n0 / (t[ii[2]] - t[ii[0]]) if t[ii[1]] != t[ii[3]]: n1 = bspl_der((t, k - 1, i + 1), x, order - 1) - val1 = k*n1/(t[ii[3]] - t[ii[1]]) + val1 = k * n1 / (t[ii[3]] - t[ii[1]]) val = val0 - val1 @@ -84,6 +128,22 @@ def bspl_der(params: BsplParams, x: NDArray, order: int) -> NDArray: @cache_bspl def bspl_int(params: BsplParams, x: NDArray, order: int) -> NDArray: + """Definite integral of the bspline function. + + Parameters + ---------- + params + Bspline function parameters as a tuple including, knots, degree and the + index of the spline basis. + x + Data points. + + Returns + ------- + describe + Definite integral of the bspline function. + + """ # knots, degree, and index t, k, i = params @@ -101,14 +161,14 @@ def bspl_int(params: BsplParams, x: NDArray, order: int) -> NDArray: if t[ii[0]] != t[ii[2]]: val0 = ( - (x - t[ii[0]])*bspl_int((t, k - 1, i), x, order) + - order*bspl_int((t, k - 1, i), x, order - 1) - )/(t[ii[2]] - t[ii[0]]) + (x - t[ii[0]]) * bspl_int((t, k - 1, i), x, order) + + order * bspl_int((t, k - 1, i), x, order - 1) + ) / (t[ii[2]] - t[ii[0]]) if t[ii[1]] != t[ii[3]]: val1 = ( - (t[ii[3]] - x)*bspl_int((t, k - 1, i + 1), x, order) - - order*bspl_int((t, k - 1, i + 1), x, order - 1) - )/(t[ii[3]] - t[ii[1]]) + (t[ii[3]] - x) * bspl_int((t, k - 1, i + 1), x, order) + - order * bspl_int((t, k - 1, i + 1), x, order - 1) + ) / (t[ii[3]] - t[ii[1]]) val = val0 + val1 @@ -116,6 +176,10 @@ def bspl_int(params: BsplParams, x: NDArray, order: int) -> NDArray: def clear_bspl_cache() -> None: + """Clear all cache of the value, derivative and definite integral for + bspline function. + + """ bspl_val.cache_clear() bspl_der.cache_clear() bspl_int.cache_clear() @@ -149,6 +213,19 @@ def __init__(self, params: BsplParams) -> None: def get_bspl_funs(knots: tuple[float, ...], degree: int) -> tuple[Bspl]: - return tuple( - Bspl((knots, degree, i)) for i in range(-degree, len(knots) - 1) - ) + """Create the bspline basis functions give knots and degree. + + Parameters + ---------- + knots + Bspline knots. + degree + Bspline degree. + + Returns + ------- + describe + A full set of bspline functions. + + """ + return tuple(Bspl((knots, degree, i)) for i in range(-degree, len(knots) - 1)) diff --git a/src/xspline/indi.py b/src/xspline/indi.py index 10c3e58..c391b7f 100644 --- a/src/xspline/indi.py +++ b/src/xspline/indi.py @@ -1,13 +1,28 @@ from math import factorial import numpy as np -from numpy.typing import NDArray -from xspline.typing import IndiParams +from xspline.typing import IndiParams, NDArray from xspline.xfunction import BundleXFunction def indi_val(params: IndiParams, x: NDArray) -> NDArray: + """Indicator value function, + + Parameters + ---------- + params + Indicator parameters as a tuple consists of lower and upper bound of + the interval corresponding to the indicator function. + x + Data points. + + Returns + ------- + describe + Indicator function value. + + """ # lower and upper bounds lb, ub = params @@ -22,21 +37,54 @@ def indi_val(params: IndiParams, x: NDArray) -> NDArray: def indi_der(params: IndiParams, x: NDArray, order: int) -> NDArray: + """Indicator derivative function. Since indicator function is a piecewise + constant function, its derivative will always be zero. + + Parameters + ---------- + params + Indicator parameters as a tuple consists of lower and upper bound of + the interval corresponding to the indicator function. + x + Data points. + + Returns + ------- + describe + Indicator deviative value. + + """ return np.zeros(x.size, dtype=x.dtype) def indi_int(params: IndiParams, x: NDArray, order: int) -> NDArray: + """Indicator definite integral function. It is a piecewise polynomial + function. + + Parameters + ---------- + params + Indicator parameters as a tuple consists of lower and upper bound of + the interval corresponding to the indicator function. + x + Data points. + + Returns + ------- + describe + Indicator definite integral value. + + """ # lower and upper bounds lb, ub = params val = np.zeros(x.size, dtype=x.dtype) ind0 = (x >= lb[0]) & (x <= ub[0]) ind1 = x > ub[0] - val[ind0] = (x[ind0] - lb[0])**(-order)/factorial(-order) + val[ind0] = (x[ind0] - lb[0]) ** (-order) / factorial(-order) for i in range(-order): - val[ind1] += ( - ((ub[0] - lb[0])**(-order - i)/factorial(-order - i)) * - ((x[ind1] - ub[0])**i/factorial(i)) + val[ind1] += ((ub[0] - lb[0]) ** (-order - i) / factorial(-order - i)) * ( + (x[ind1] - ub[0]) ** i / factorial(i) ) return val diff --git a/src/xspline/poly.py b/src/xspline/poly.py index 7881939..ec250d4 100644 --- a/src/xspline/poly.py +++ b/src/xspline/poly.py @@ -1,21 +1,70 @@ from math import factorial import numpy as np -from numpy.typing import NDArray -from xspline.typing import PolyParams +from xspline.typing import PolyParams, NDArray from xspline.xfunction import BundleXFunction, XFunction def poly_val(params: PolyParams, x: NDArray) -> NDArray: + """Polynomial value function. + + Parameters + ---------- + params + Polynomial coefficients. + x + Data points. + + Returns + ------- + describe + Polynomial function values. + + """ return np.polyval(params, x) def poly_der(params: PolyParams, x: NDArray, order: int) -> NDArray: + """Polynomial derivative function. + + Parameters + ---------- + params + Polynomial coefficients. + x + Data points. + order + Order of differentiation. + + Returns + ------- + describe + Polynomial derivative values. + + """ return np.polyval(np.polyder(params, order), x) def poly_int(params: PolyParams, x: NDArray, order: int) -> NDArray: + """Polynomial definite integral function. + + Parameters + ---------- + params + Polynomial coefficients. + x + Data points. + order + Order of integration. Here we use negative integer. + + Returns + ------- + describe + Polynomial definite integral values. + + """ + return np.polyval(np.polyint(params, -order), x) @@ -46,16 +95,50 @@ def __init__(self, params: PolyParams) -> None: def get_poly_params(fun: XFunction, x: float, degree: int) -> tuple[float, ...]: + """Solve polynomial (taylor) coefficients provided the ``XFunction``. + + Parameters + ---------- + fun + Provided ``XFunction`` to be approximated. + x + The point where we want to approximate ``XFunction`` by the polynomial. + degree + Degree of the approximation polynomial. + + Returns + ------- + describe + The approximation polynomial coefficients. + + """ if degree < 0: return (0.0,) rhs = np.array([fun(x, order=i) for i in range(degree, -1, -1)]) mat = np.zeros((degree + 1, degree + 1)) for i in range(degree + 1): for j in range(i + 1): - mat[i, j] = factorial(degree - j) / factorial(i - j) * x**(i - j) + mat[i, j] = factorial(degree - j) / factorial(i - j) * x ** (i - j) return tuple(np.linalg.solve(mat, rhs)) def get_poly_fun(fun: XFunction, x: float, degree: int) -> Poly: + """Get the approximation polynomial function. + + Parameters + ---------- + fun + Provided ``XFunction`` to be approximated. + x + The point where we want to approximate ``XFunction`` by the polynomial. + degree + Degree of the approximation polynomial. + + Returns + ------- + describe + Instance of the ``Poly`` class to approximate provided ``XFunction``. + + """ params = get_poly_params(fun, x, degree) return Poly(params) diff --git a/src/xspline/xfunction.py b/src/xspline/xfunction.py index cabc3a1..aed0205 100644 --- a/src/xspline/xfunction.py +++ b/src/xspline/xfunction.py @@ -1,34 +1,58 @@ +from __future__ import annotations from functools import partial from math import factorial from operator import attrgetter -from typing import Optional import numpy as np -from numpy.typing import NDArray -from xspline.typing import (BoundaryPoint, Callable, RawDFunction, - RawIFunction, RawVFunction) +from xspline.typing import ( + BoundaryPoint, + Callable, + RawDFunction, + RawIFunction, + RawVFunction, + NDArray, +) class XFunction: + """Function interface that provide easy access to function value, + derivatives and definite integrals. - def __init__(self, fun: Optional[Callable] = None) -> None: + There are two different ways to use this class, either use this class to + wrap around a function that aligns with the ``XFunction`` function call + interface. The added benefit is ``XFunction`` will automatic check and parse + the inputs. The other way is to inherit this class and implement a ``fun`` + member function for the function call. + + Parameters + ---------- + fun + Optional function implementation. If this is ``None``, it will require + user to implement ``fun`` class member function. + + """ + + def __init__(self, fun: Callable | None = None) -> None: if not hasattr(self, "fun"): self.fun = fun def _check_args(self, x: NDArray, order: int) -> tuple[NDArray, int, bool]: x, order = np.asarray(x, dtype=float), int(order) if (x.ndim not in [0, 1, 2]) or (x.ndim == 2 and len(x) != 2): - raise ValueError("please provide a scalar, an 1d array, or a 2d " - "array with two rows") + raise ValueError( + "please provide a scalar, an 1d array, or a 2d " "array with two rows" + ) # reshape array isscalar = x.ndim == 0 if isscalar: x = x.ravel() if order >= 0 and x.ndim == 2: - raise ValueError("please provide an 1d array for function value " - "defivative computation") + raise ValueError( + "please provide an 1d array for function value " + "defivative computation" + ) if order < 0 and x.ndim == 1: x = np.vstack([np.repeat(x.min(), x.size), x]) @@ -44,7 +68,7 @@ def __call__(self, x: NDArray, order: int = 0) -> NDArray: Parameters ---------- x - Data points where the function is evaluated. If `order < 0` and `x` + Data points where the function is evaluated. If `order < 0` and `x` is a 2d array with two rows, the rows will be treated as the starting and ending points for definite interval. If `order < 0` and `x` is a 1d array, function will use the smallest number in `x` as @@ -85,7 +109,19 @@ def __call__(self, x: NDArray, order: int = 0) -> NDArray: result = result[0] return result - def append(self, other: "XFunction", sep: BoundaryPoint) -> "XFunction": + def append(self, other: XFunction, sep: BoundaryPoint) -> XFunction: + """Splice with another instance of ``XFunction`` to create a new + ``XFunction``. + + Parameters + ---------- + other + Another ``XFunction`` after the current function. + sep + The boundary point to separate two functions, before is the current + function and after the the ``other`` function. + + """ def fun(x: NDArray, order: int = 0) -> NDArray: left = x <= sep[0] if sep[1] else x < sep[0] @@ -114,18 +150,53 @@ def fun(x: NDArray, order: int = 0) -> NDArray: class BundleXFunction(XFunction): - - def __init__(self, - params: tuple, - val_fun: RawVFunction, - der_fun: RawDFunction, - int_fun: RawIFunction) -> None: + """This is one implementation of the ``XFunction``, it takes the value, + derivative and definite integral function and bundle them together as a + ``XFunction``. + + Parameters + ---------- + params + This is the parameters that is needed for the value, derivatives and + the definitely integral function. + val_fun + Value function. + der_fun + Derviative function. + int_fun + Defintie integral function. + + """ + + def __init__( + self, + params: tuple, + val_fun: RawVFunction, + der_fun: RawDFunction, + int_fun: RawIFunction, + ) -> None: self.params = params self.val_fun = partial(val_fun, params) self.der_fun = partial(der_fun, params) self.int_fun = partial(int_fun, params) def fun(self, x: NDArray, order: int = 0) -> NDArray: + """Function implementation, aligns with the ``XFunction`` function call + interface. + + Parameters + ---------- + x + Data points + order + Order of differentiation/integration. + + Returns + ------- + describe + Function value, derivatives or definite integrals. + + """ if order == 0: return self.val_fun(x) if order > 0: @@ -138,40 +209,93 @@ def fun(self, x: NDArray, order: int = 0) -> NDArray: class BasisXFunction(XFunction): + """This is one implementation of ``XFunction`` by taking in a set of + instances of ``XFunction`` as basis functions. And the linear combination + coefficients to provide function value, derivative and definite integral. + + Parameters + ---------- + basis_funs + A set of instances of ``XFunction`` as basis functions. + coefs + Coefficients for the linearly combine the basis functions. + + """ coefs = property(attrgetter("_coefs")) - def __init__(self, - basis_funs: tuple[XFunction, ...], - coefs: Optional[NDArray] = None) -> None: + def __init__( + self, basis_funs: tuple[XFunction, ...], coefs: NDArray | None = None + ) -> None: if not all(isinstance(fun, XFunction) for fun in basis_funs): - raise TypeError("basis functions must all be instances of " - "`XFunction`") + raise TypeError("basis functions must all be instances of " "`XFunction`") self.basis_funs = tuple(basis_funs) self.coefs = coefs @coefs.setter - def coefs(self, coefs: Optional[NDArray]) -> None: + def coefs(self, coefs: NDArray | None) -> None: if coefs is not None: coefs = np.asarray(coefs, dtype=float).ravel() if coefs.size != len(self): - raise ValueError("number of coeffcients does not match number " - "of basis functions") + raise ValueError( + "number of coeffcients does not match number " "of basis functions" + ) self._coefs = coefs - def get_design_mat(self, x: NDArray, - order: int = 0, - check_args: bool = True) -> NDArray: + def get_design_mat( + self, x: NDArray, order: int = 0, check_args: bool = True + ) -> NDArray: + """Provide design matrix from the set of basis functions. + + Parameters + ---------- + x + Data points + order + Order of differentiation/integration. + check_args + If ``True`` it will check and parse the arguments. + + Returns + ------- + describe + Design matrix with dimention number of data points by number of + basis functions. + + """ if check_args: x, order, _ = self._check_args(x, order) return np.vstack([fun.fun(x, order) for fun in self.basis_funs]).T def fun(self, x: NDArray, order: int = 0) -> NDArray: + """Function implementation, aligns with the ``XFunction`` function call + interface. + + Parameters + ---------- + x + Data points + order + Order of differentiation/integration. + + Returns + ------- + describe + Function value, derivatives or definite integrals. + + Raises + ------ + ValueError + Raised when the ``coefs`` is not provided. + + """ if self.coefs is None: - raise ValueError("please provide the coefficients for the basis " - "functions") + raise ValueError( + "please provide the coefficients for the basis " "functions" + ) design_mat = self.get_design_mat(x, order=order, check_args=False) return design_mat.dot(self.coefs) def __len__(self) -> int: + """Number of basis functions.""" return len(self.basis_funs) diff --git a/src/xspline/xspl.py b/src/xspline/xspl.py index 991d48d..4888e61 100644 --- a/src/xspline/xspl.py +++ b/src/xspline/xspl.py @@ -1,20 +1,37 @@ from typing import Optional -from numpy.typing import NDArray - from xspline.bspl import clear_bspl_cache, get_bspl_funs from xspline.poly import get_poly_fun from xspline.xfunction import BasisXFunction +from xspline.typing import NDArray class XSpline(BasisXFunction): + """Main class for xspline functions. + + Parameters + ---------- + knots + Knots of the spline. + degree + Degree of the spline. + ldegree + Left extrapolation polynomial degree. + rdegree + Right extrapolation polynomial degree. + coefs + The coefficients for linear combining the spline basis. + + """ - def __init__(self, - knots: tuple[float, ...], - degree: int, - ldegree: Optional[int] = None, - rdegree: Optional[int] = None, - coefs: Optional[NDArray] = None) -> None: + def __init__( + self, + knots: tuple[float, ...], + degree: int, + ldegree: Optional[int] = None, + rdegree: Optional[int] = None, + coefs: Optional[NDArray] = None, + ) -> None: # validate inputs knots, degree = tuple(sorted(map(float, knots))), int(degree) if len(set(knots)) < 2: @@ -37,8 +54,26 @@ def __init__(self, self.ldegree, self.rdegree = ldegree, rdegree super().__init__(funs, coefs=coefs) - def get_design_mat(self, x: NDArray, - order: int = 0, check_args: bool = True) -> NDArray: + def get_design_mat( + self, x: NDArray, order: int = 0, check_args: bool = True + ) -> NDArray: + """Create design matrix from spline basis functions. + + Parameters + ---------- + x + Data points. + order + Order of differentiation/integration. + check_args + If ``True``, it will automatically check and parse the arguments. + + Returns + ------- + describe + Design matrix from spline basis functions. + + """ design_mat = super().get_design_mat(x, order, check_args) clear_bspl_cache() return design_mat