Bspline

docs/developers_notes/basis_module.md

@@ -13,7 +13,11 @@ Abstract Class Basis
 │
 ├─ Abstract Subclass SplineBasis
 │   │
-│   └─ Concrete Subclass MSplineBasis
+│   ├─ Concrete Subclass MSplineBasis
+│   │
+│   ├─ Concrete Subclass BSplineBasis
+│   │
+│   └─ Concrete Subclass CyclicBSplineBasis
 │
 ├─ Abstract Subclass RaisedCosineBasis
 │   │

src/neurostatslib/basis.py

@@ -9,11 +9,14 @@
 import numpy as np
 import scipy.linalg
 from numpy.typing import NDArray
+from scipy.interpolate import splev
 
 from neurostatslib.utils import row_wise_kron
 
 __all__ = [
     "MSplineBasis",
+    "BSplineBasis",
+    "CyclicBSplineBasis",
     "RaisedCosineBasisLinear",
     "RaisedCosineBasisLog",
     "OrthExponentialBasis",
@@ -91,12 +94,6 @@ def evaluate(self, *xi: NDArray) -> NDArray:
         -------
         :
             The generated basis functions.
-
-        Raises
-        ------
-        ValueError
-            If the time point number is inconsistent between inputs or if the number of inputs doesn't match what
-            the Basis object requires.
         """
         # checks on input and outputs
         self._check_samples_consistency(*xi)
@@ -465,6 +462,13 @@ def _generate_knots(
         )
         return self.knot_locs
 
+    @staticmethod
+    def _check_samples_non_empty(sample_pts):
+        if sample_pts.shape[0] == 0:
+            raise ValueError(
+                "Empty sample array provided. At least one sample is required for evaluation!"
+            )
+
 
 class MSplineBasis(SplineBasis):
     """M-spline 1-dimensional basis functions.
@@ -533,6 +537,218 @@ def _check_n_basis_min(self) -> None:
             )
 
 
+class BSplineBasis(SplineBasis):
+    """
+    B-spline 1-dimensional basis functions.
+
+    Parameters
+    ----------
+    n_basis_funcs :
+        Number of basis functions.
+    order :
+        Order of the splines used in basis functions. Must lie within [1, n_basis_funcs].
+        The B-splines have (order-2) continuous derivatives at each interior knot.
+        The higher this number, the smoother the basis representation will be.
+
+    Attributes
+    ----------
+    order :
+        Spline order.
+
+    Methods
+    -------
+    _evaluate(x_tuple)
+       Evaluate the basis function at the samples x_tuple[0]. x_tuple must be of length 1 in order to pass the checks
+       of super().evaluate
+
+    References
+    ----------
+    ..[2] Prautzsch, H., Boehm, W., Paluszny, M. (2002). B-spline representation. In: Bézier and B-Spline Techniques.
+    Mathematics and Visualization. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-04919-8_5
+
+    """
+
+    def __init__(self, n_basis_funcs: int, order: int = 2):
+        super().__init__(n_basis_funcs, order=order)
+
+    def _evaluate(self, sample_pts: NDArray) -> NDArray:
+        """
+        Evaluate the B-spline basis functions with given sample points.
+
+        Parameters
+        ----------
+        sample_pts :
+            The sample points at which the B-spline is evaluated.
+
+        Returns
+        -------
+        NDArray
+            The basis function evaluated at the samples (Time points x number of basis)
+
+        Raises
+        ------
+        AssertionError
+            If the sample points are not within the B-spline knots range unless `outer_ok=True`.
+
+        Notes
+        -----
+        This method evaluates the B-spline basis functions at the given sample points.
+        It requires the knots to be defined through the `_generate_knots` method.
+
+        The evaluation is performed by looping over each element and using `splev`
+        from SciPy to compute the basis values.
+        """
+        super()._check_samples_non_empty(sample_pts)
+
+        # add knots
+        self._generate_knots(sample_pts, 0.0, 1.0)
+
+        basis_eval = bspline(
+            sample_pts, self.knot_locs, order=self.order, der=0, outer_ok=False
+        )
+
+        return basis_eval
+
+    def _check_n_basis_min(self) -> None:
+        """Check that the user required enough basis elements.
+
+        Check that BSplineBasis has at least as many basis as the order of the spline.
+
+        Raises
+        ------
+        ValueError
+            If an insufficient number of basis element is requested for the basis type
+        """
+        if self.n_basis_funcs < self.order:
+            raise ValueError(
+                f"{self.__class__.__name__} `order` parameter cannot be larger "
+                "than `n_basis_funcs` parameter."
+            )
+
+
+class CyclicBSplineBasis(SplineBasis):
+    """
+    B-spline 1-dimensional basis functions for cyclic splines.
+
+    Parameters
+    ----------
+    n_basis_funcs :
+        Number of basis functions.
+    order :
+        Order of the splines used in basis functions. Must lie within [1, n_basis_funcs].
+        The B-splines have (order-2) continuous derivatives at each interior knot.
+        The higher this number, the smoother the basis representation will be.
+
+    Attributes
+    ----------
+    n_basis_funcs : int
+        Number of basis functions.
+    order : int
+        Order of the splines used in basis functions.
+
+    Methods
+    -------
+    _evaluate(sample_pts, der)
+        Evaluate the B-spline basis functions with given sample points.
+    """
+
+    def __init__(self, n_basis_funcs: int, order: int = 2):
+        super().__init__(n_basis_funcs, order=order)
+        if self.order < 2:
+            raise ValueError(
+                f"Order >= 2 required for cyclic B-spline, "
+                f"order {self.order} specified instead!"
+            )
+
+    def _evaluate(self, sample_pts: NDArray) -> NDArray:
+        """
+        Evaluate the B-spline basis functions with given sample points.
+
+        Parameters
+        ----------
+        sample_pts :
+            The sample points at which the B-spline is evaluated. Must be a tuple of length 1.
+
+        Returns
+        -------
+        NDArray
+            The basis function evaluated at the samples (Time points x number of basis)
+
+        Raises
+        ------
+        AssertionError
+            If the sample points are not within the B-spline knots range unless `outer_ok=True`.
+
+        Notes
+        -----
+        This method evaluates the B-spline basis functions at the given sample points.
+        It requires the knots to be defined through the `_generate_knots` method.
+
+        The evaluation is performed by looping over each element and using `splev` from
+        SciPy to compute the basis values.
+        """
+        super()._check_samples_non_empty(sample_pts)
+
+        self._generate_knots(sample_pts, 0.0, 1.0, is_cyclic=True)
+
+        # for cyclic, do not repeat knots
+        self.knot_locs = np.unique(self.knot_locs)
+
+        knots_orig = self.knot_locs.copy()
+
+        nk = knots_orig.shape[0]
+
+        # make sure knots are sorted
+        knots_orig.sort()
+
+        # extend knots
+        xc = knots_orig[nk - 2 * self.order + 1]
+        knots = np.hstack(
+            (
+                self.knot_locs[0]
+                - self.knot_locs[-1]
+                + self.knot_locs[nk - self.order : nk - 1],
+                self.knot_locs,
+            )
+        )
+        ind = sample_pts > xc
+
+        # temporarily set the extended knots as attribute
+        self.knot_locs = knots
+
+        basis_eval = bspline(
+            sample_pts, self.knot_locs, order=self.order, der=0, outer_ok=True
+        )
+        sample_pts[ind] = sample_pts[ind] - knots.max() + knots_orig[0]
+
+        if np.sum(ind):
+            basis_eval[ind] = basis_eval[ind] + bspline(
+                sample_pts[ind], self.knot_locs, order=self.order, outer_ok=True, der=0
+            )
+        # restore points
+        sample_pts[ind] = sample_pts[ind] + knots.max() - knots_orig[0]
+        # restore the original knots
+        self.knot_locs = knots_orig
+
+        return basis_eval
+
+    def _check_n_basis_min(self) -> None:
+        """Check that the user required enough basis elements.
+
+        Check that Cuclic-BSplineBasis has at least as many basis as the order of the spline +2
+        and at least 2*order - 2 basis.
+
+        Raises
+        ------
+        ValueError
+            If an insufficient number of basis element is requested for the basis type
+        """
+        if self.n_basis_funcs < max(self.order * 2 - 2, self.order + 2):
+            raise ValueError(
+                f"Insufficient basis elements for {self.__class__.__name__} instantiation."
+            )
+
+
 class RaisedCosineBasis(Basis, abc.ABC):
     def __init__(self, n_basis_funcs: int) -> None:
         super().__init__(n_basis_funcs)
@@ -870,3 +1086,83 @@ def mspline(x: NDArray, k: int, i: int, T: NDArray):
             )
             / ((k - 1) * (T[i + k] - T[i]))
         )
+
+
+def bspline(
+    sample_pts: NDArray,
+    knots: NDArray,
+    order: int = 4,
+    der: int = 0,
+    outer_ok: bool = False,
+):
+    """
+    Calculate and return the evaluation of B-spline basis.
+
+    This function evaluates B-spline basis for given sample points. It checks for
+    out of range points and optionally handles them. It also handles the NaNs if present.
+
+    Parameters
+    ----------
+    sample_pts :
+        An array containing sample points for which B-spline basis needs to be evaluated.
+    knots :
+        An array containing knots for the B-spline basis. The knots are sorted in ascending order.
+    order :
+        The order of the B-spline basis. Default is 4.
+    der :
+        The derivative of the B-spline basis to be evaluated. Default is 0.
+    outer_ok :
+        If True, allows for evaluation at points outside the range of knots.
+        Default is False, in which case an assertion error is raised when
+        points outside the knots range are encountered.
+
+    Returns
+    -------
+    basis_eval :
+        An array containing the evaluation of B-spline basis for the given sample points.
+        Shape (n_samples, n_basis_funcs).
+
+    Raises
+    ------
+    AssertionError
+        If `outer_ok` is False and the sample points lie outside the B-spline knots range.
+
+    Notes
+    -----
+    The function uses splev function from scipy.interpolate library for the basis evaluation.
+    """
+    knots.sort()
+    nk = knots.shape[0]
+
+    # check for out of range points (in cyclic b-spline need_outer must be set to False)
+    need_outer = any(sample_pts < knots[order - 1]) or any(
+        sample_pts > knots[nk - order]
+    )
+    assert (
+        not need_outer
+    ) | outer_ok, 'sample points must lie within the B-spline knots range unless "outer_ok==True".'
+
+    # select knots that are within the knots range (this takes care of eventual NaNs)
+    in_sample = (sample_pts >= knots[0]) & (sample_pts <= knots[-1])
+
+    if need_outer:
+        reps = order - 1
+        knots = np.hstack((np.ones(reps) * knots[0], knots, np.ones(reps) * knots[-1]))
+        nk = knots.shape[0]
+    else:
+        reps = 0
+
+    # number of basis elements
+    n_basis = nk - order
+
+    # initialize the basis element container
+    basis_eval = np.zeros((n_basis - 2 * reps, sample_pts.shape[0]))
+
+    # loop one element at the time and evaluate the basis using splev
+    id_basis = np.eye(n_basis, nk, dtype=np.int8)
+    for i in range(reps, len(knots) - order - reps):
+        basis_eval[i - reps, in_sample] = splev(
+            sample_pts[in_sample], (knots, id_basis[i], order - 1), der=der
+        )
+
+    return basis_eval.T