Merge bdd9488 into dd19cf4

CHANGELOG.md

@@ -1,5 +1,11 @@
 # Changelog
 
+## v0.10.0
+
+* Significant performance improvements for make/evaluate splines
+* Change format for storing spline coefficients (reshape coeffs array) to improve performance
+* Add shape property to `SplinePPForm`/`NdGridSplinePPForm` and axis property to `SplinePPForm`
+
 ## v0.9.0
 
 * Drop support of Python 3.5

csaps/_base.py

@@ -41,24 +41,24 @@ def coeffs(self) -> np.ndarray:
 
     @property
     @abc.abstractmethod
-    def pieces(self) -> TProps:
-        """Returns a spline pieces data
+    def order(self) -> TProps:
+        """Returns a spline order
 
         Returns
         -------
-        pieces : ty.Union[int, ty.Tuple[int, ...]]
-            Returns a spline pieces data
+        order : ty.Union[int, ty.Tuple[int, ...]]
+            Returns a spline order
         """
 
     @property
     @abc.abstractmethod
-    def order(self) -> TProps:
-        """Returns a spline order
+    def pieces(self) -> TProps:
+        """Returns a spline pieces data
 
         Returns
         -------
-        order : ty.Union[int, ty.Tuple[int, ...]]
-            Returns a spline order
+        pieces : ty.Union[int, ty.Tuple[int, ...]]
+            Returns a spline pieces data
         """
 
     @property
@@ -91,7 +91,7 @@ def __repr__(self):  # pragma: no cover
         return (
             f'{type(self).__name__}\n'
             f'  breaks: {self.breaks}\n'
-            f'  coeffs: {self.coeffs.shape} shape\n'
+            f'  coeffs: shape {self.coeffs.shape}\n'
             f'  pieces: {self.pieces}\n'
             f'  order: {self.order}\n'
             f'  ndim: {self.ndim}\n'

csaps/_sspndg.py

@@ -40,16 +40,21 @@ class NdGridSplinePPForm(SplinePPFormBase[ty.Sequence[np.ndarray], ty.Tuple[int,
     breaks : np.ndarray
         Breaks values 1-D array
     coeffs : np.ndarray
-        Spline coefficients 2-D array
+        Spline coefficients N-D array
     """
 
     def __init__(self, breaks: ty.Sequence[np.ndarray], coeffs: np.ndarray) -> None:
         self._breaks = breaks
         self._coeffs = coeffs
         self._pieces = tuple(x.size - 1 for x in breaks)
-        self._order = tuple(s // p for s, p in zip(coeffs.shape[1:], self._pieces))
         self._ndim = len(breaks)
 
+        if self._ndim > 1:
+            self._order = tuple(s // p for s, p in zip(coeffs.shape, self._pieces))
+        else:
+            # the corner case for univariate spline that is represented as 1d-grid
+            self._order = (coeffs.shape[1] // self._pieces[0], )
+
     @property
     def breaks(self) -> ty.Sequence[np.ndarray]:
         return self._breaks
@@ -58,36 +63,58 @@ def breaks(self) -> ty.Sequence[np.ndarray]:
     def coeffs(self) -> np.ndarray:
         return self._coeffs
 
-    @property
-    def pieces(self) -> ty.Tuple[int, ...]:
-        return self._pieces
-
     @property
     def order(self) -> ty.Tuple[int, ...]:
         return self._order
 
+    @property
+    def pieces(self) -> ty.Tuple[int, ...]:
+        return self._pieces
+
     @property
     def ndim(self) -> int:
         return self._ndim
 
+    @property
+    def shape(self) -> ty.Tuple[int, ...]:
+        """Returns the original data shape
+        """
+        return tuple(x.size for x in self.breaks)
+
     def evaluate(self, xi: ty.Sequence[np.ndarray]) -> np.ndarray:
-        yi = self.coeffs.copy()
-        sizey = list(yi.shape)
-        nsize = tuple(x.size for x in xi)
+        shape = tuple(x.size for x in xi)
+
+        coeffs = self.coeffs
+        coeffs_shape = list(coeffs.shape)
+
+        d = self.ndim - 1
+        permuted_axes = (d, *range(d))
 
         for i in reversed(range(self.ndim)):
-            ndim = int(np.prod(sizey[:self.ndim]))
-            coeffs = yi.reshape((ndim * self.pieces[i], self.order[i]), order='F')
+            xii = xi[i]
+            ndim = int(np.prod(coeffs_shape[:d]))
 
-            spp = SplinePPForm(self.breaks[i], coeffs, ndim=ndim, shape=(ndim, xi[i].size))
-            yi = spp.evaluate(xi[i])
+            if self.ndim > 2:
+                coeffs = coeffs.reshape((ndim, self.pieces[i] * self.order[i]))
 
-            yi = yi.reshape((*sizey[:self.ndim], nsize[i]), order='F')
-            axes = (0, self.ndim, *range(1, self.ndim))
-            yi = yi.transpose(axes)
-            sizey = list(yi.shape)
+            spp = SplinePPForm(
+                breaks=self.breaks[i],
+                coeffs=coeffs,
+                pieces=self.pieces[i],
+                order=self.order[i],
+                shape=(ndim, xii.size)
+            )
 
-        return yi.reshape(nsize, order='F')
+            coeffs = spp.evaluate(xii)
+
+            if self.ndim > 2:
+                coeffs = coeffs.reshape((*coeffs_shape[:d], shape[i]))
+
+            if self.ndim > 1:
+                coeffs = coeffs.transpose(permuted_axes)
+                coeffs_shape = list(coeffs.shape)
+
+        return coeffs.reshape(shape)
 
 
 class NdGridCubicSmoothingSpline(ISmoothingSpline[NdGridSplinePPForm, ty.Tuple[float, ...], NdGridDataType]):
@@ -209,25 +236,29 @@ def __call__(self, xi: NdGridDataType) -> np.ndarray:
         return self._spline.evaluate(xi)
 
     def _make_spline(self, smooth: ty.List[ty.Optional[float]]) -> ty.Tuple[NdGridSplinePPForm, ty.Tuple[float, ...]]:
-        sizey = [1] + list(self._ydata.shape)
-        ydata = self._ydata.reshape(sizey, order='F').copy()
-        _smooth = []
+        coeffs = self._ydata
+        coeffs_shape = list(coeffs.shape)
 
-        # Perform coordinatewise smoothing spline computing
+        smooths = []
+        permute_axes = (self._ndim - 1, *range(self._ndim - 1))
+
+        # computing coordinatewise smoothing spline
         for i in reversed(range(self._ndim)):
-            shape_i = (np.prod(sizey[:-1]), sizey[-1])
-            ydata_i = ydata.reshape(shape_i, order='F')
+            if self._ndim > 2:
+                coeffs = coeffs.reshape(np.prod(coeffs.shape[:-1]), coeffs.shape[-1])
 
             s = CubicSmoothingSpline(
-                self._xdata[i], ydata_i, weights=self._weights[i], smooth=smooth[i])
+                self._xdata[i], coeffs, weights=self._weights[i], smooth=smooth[i])
+
+            smooths.append(s.smooth)
+            coeffs = s.spline.coeffs
 
-            _smooth.append(s.smooth)
-            sizey[-1] = s.spline.pieces * s.spline.order
-            ydata = s.spline.coeffs.reshape(sizey, order='F')
+            if self._ndim > 2:
+                coeffs_shape[-1] = s.spline.pieces * s.spline.order
+                coeffs = coeffs.reshape(coeffs_shape)
 
             if self._ndim > 1:
-                axes = (0, self._ndim, *range(1, self._ndim))
-                ydata = ydata.transpose(axes)
-                sizey = list(ydata.shape)
+                coeffs = coeffs.transpose(permute_axes)
+                coeffs_shape = list(coeffs.shape)
 
-        return NdGridSplinePPForm(self._xdata, ydata), tuple(_smooth)
+        return NdGridSplinePPForm(breaks=self._xdata, coeffs=coeffs), tuple(smooths)

csaps/_sspumv.py

@@ -38,15 +38,15 @@ class SplinePPForm(SplinePPFormBase[np.ndarray, int]):
         Axis along which values are assumed to be varying
     """
 
-    def __init__(self, breaks: np.ndarray, coeffs: np.ndarray, ndim: int,
+    def __init__(self, breaks: np.ndarray, coeffs: np.ndarray, order: int, pieces: int,
                  shape: ty.Sequence[int], axis: int = -1) -> None:
         self._breaks = breaks
         self._coeffs = coeffs
-        self._pieces = np.prod(coeffs.shape[:-1]) // ndim
-        self._order = coeffs.shape[-1]
-        self._ndim = ndim
+        self._order = order
+        self._pieces = pieces
+        self._ndim = coeffs.shape[0]
 
-        self._shape = shape
+        self._shape = tuple(shape)
         self._axis = axis
 
     @property
@@ -57,21 +57,33 @@ def breaks(self) -> np.ndarray:
     def coeffs(self) -> np.ndarray:
         return self._coeffs
 
-    @property
-    def pieces(self) -> int:
-        return self._pieces
-
     @property
     def order(self) -> int:
         return self._order
 
+    @property
+    def pieces(self) -> int:
+        return self._pieces
+
     @property
     def ndim(self) -> int:
         return self._ndim
 
+    @property
+    def shape(self) -> ty.Tuple[int, ...]:
+        """Returns the original data shape
+        """
+        return self._shape
+
+    @property
+    def axis(self) -> int:
+        """Returns the data axis along the spline will be evaluated
+        """
+        return self._axis
+
     def evaluate(self, xi: np.ndarray) -> np.ndarray:
-        shape = list(self._shape)
-        shape[self._axis] = xi.size
+        shape = list(self.shape)
+        shape[self.axis] = xi.size
 
         # For each data site, compute its break interval
         mesh = self.breaks[1:-1]
@@ -82,34 +94,21 @@ def evaluate(self, xi: np.ndarray) -> np.ndarray:
         nanx = np.flatnonzero(index == 0)
         index = np.fmin(index, mesh.size + 1)
         index[nanx] = 1
+        index -= 1
 
         # Go to local coordinates
-        xi = xi - self.breaks[index - 1]
-        d = self.ndim
-        lx = len(xi)
-
-        if d > 1:
-            xi_shape = (1, d * lx)
-            xi_ndm = np.array(xi, ndmin=2)
-            xi = np.reshape(np.repeat(xi_ndm, d, axis=0), xi_shape, order='F')
-
-            index_rep = (np.repeat(np.array(1 + d * index, ndmin=2), d, axis=0)
-                         + np.repeat(np.array(np.r_[-d:0], ndmin=2).T, lx, axis=1))
-            index = np.reshape(index_rep, (d * lx, 1), order='F')
-
-        index -= 1
+        xi = xi - self.breaks[index]
 
         # Apply nested multiplication
-        values = self._coeffs[index, 0].T
+        values = self.coeffs[:, index]
 
-        for i in range(1, self._coeffs.shape[1]):
-            values = xi * values + self._coeffs[index, i].T
-
-        values = values.reshape((d, lx), order='F').squeeze()
+        for i in range(1, self.order):
+            index += self.pieces
+            values = xi * values + self.coeffs[:, index]
 
         if values.shape != shape:
             # Reshape values 2-D NxM array to N-D array with original shape
-            values = from_2d(values, shape, self._axis)
+            values = from_2d(values, shape, self.axis)
 
         return values
 
@@ -165,7 +164,7 @@ def __call__(self, xi: UnivariateDataType) -> np.ndarray:
         xi = ty.cast(np.ndarray, np.asarray(xi, dtype=np.float64))
 
         if xi.ndim > 1:  # pragma: no cover
-            raise ValueError('"xi" data must be a 1-d array.')
+            raise ValueError("'xi' data must be a 1-d array.")
 
         return self._spline.evaluate(xi)
 
@@ -197,16 +196,16 @@ def _prepare_data(xdata, ydata, weights, axis):
         ydata = np.asarray(ydata, dtype=np.float64)
 
         if xdata.ndim > 1:
-            raise ValueError('xdata must be a vector')
+            raise ValueError("'xdata' must be a vector")
         if xdata.size < 2:
-            raise ValueError('xdata must contain at least 2 data points.')
+            raise ValueError("'xdata' must contain at least 2 data points.")
 
         yshape = list(ydata.shape)
 
         if yshape[axis] != xdata.size:
             raise ValueError(
-                f'"ydata" data must be a 1-D or N-D array with shape[{axis}] '
-                f'that is equal to "xdata" size ({xdata.size})')
+                f"'ydata' data must be a 1-D or N-D array with shape[{axis}] "
+                f"that is equal to 'xdata' size ({xdata.size})")
 
         # Reshape ydata N-D array to 2-D NxM array where N is the data
         # dimension and M is the number of data points.
@@ -240,7 +239,8 @@ def _make_spline(self, smooth: ty.Optional[float]) -> ty.Tuple[SplinePPForm, flo
         dx = np.diff(self._xdata)
 
         if not all(dx > 0):  # pragma: no cover
-            raise ValueError('Items of xdata vector must satisfy the condition: x1 < x2 < ... < xN')
+            raise ValueError(
+                "Items of 'xdata' vector must satisfy the condition: x1 < x2 < ... < xN")
 
         dy = np.diff(self._ydata, axis=1)
         dy_dx = dy / dx
@@ -277,33 +277,39 @@ def _make_spline(self, smooth: ty.Optional[float]) -> ty.Tuple[SplinePPForm, flo
                 u = u.T
 
             dx = dx[:, np.newaxis]
-            d_pad = np.zeros((1, self._ydim))
 
-            d1 = np.diff(np.vstack((d_pad, u, d_pad)), axis=0) / dx
-            d2 = np.diff(np.vstack((d_pad, d1, d_pad)), axis=0)
+            pad_width = [(1, 1), (0, 0)]
+            d1 = np.diff(np.pad(u, pad_width), axis=0) / dx
+            d2 = np.diff(np.pad(d1, pad_width), axis=0)
 
             yi = self._ydata.T - ((6. * (1. - p)) * w) @ d2
-            c3 = np.vstack((d_pad, p * u, d_pad))
+            c3 = np.pad(p * u, pad_width)
             c2 = np.diff(yi, axis=0) / dx - dx * (2. * c3[:-1, :] + c3[1:, :])
 
-            coeffs = np.hstack((
-                (np.diff(c3, axis=0) / dx).T,
-                3. * c3[:-1, :].T,
-                c2.T,
-                yi[:-1, :].T
-            ))
+            coeffs = np.vstack((
+                np.diff(c3, axis=0) / dx,
+                3. * c3[:-1, :],
+                c2,
+                yi[:-1, :],
+            )).T
 
-            cf_shape = ((pcount - 1) * self._ydim, 4)
-            coeffs = coeffs.reshape(cf_shape, order='F')
+            order = 4
+            pieces = coeffs.shape[1] // order
         else:
-            p = 1.
+            # The corner case for the data with 2 points.
             yi = self._ydata[:, 0][:, np.newaxis]
-            coeffs = np.array(np.hstack((dy_dx, yi)), ndmin=2)
+            coeffs = np.hstack((dy_dx, yi))
+
+            order = 2
+            pieces = 1
+
+            p = 1.
 
         spline = SplinePPForm(
             breaks=self._xdata,
             coeffs=coeffs,
-            ndim=self._ydim,
+            order=order,
+            pieces=pieces,
             shape=self._shape,
             axis=self._axis
         )

csaps/_version.py

@@ -1,3 +1,3 @@
 # -*- coding: utf-8 -*-
 
-__version__ = '0.9.0'
+__version__ = '0.10.0'