diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml
index 143ed12..fbe9ce1 100644
--- a/.github/workflows/main.yaml
+++ b/.github/workflows/main.yaml
@@ -44,8 +44,10 @@ jobs:
       - name: Install the project dependencies
         run: poetry install
 
-      - name: Run static analysis with ruff
-        run: poetry run ruff check .
+      - name: Run static analysis and linters
+        run: |
+          poetry run ruff check .
+          poetry run isort -c .
 
       - name: Run tests
         run: poetry run pytest -v --color=yes --cov=csaps --cov-report=term --cov-report=lcov:coverage.info
diff --git a/csaps/__init__.py b/csaps/__init__.py
index bd0291f..1e507f7 100644
--- a/csaps/__init__.py
+++ b/csaps/__init__.py
@@ -1,30 +1,13 @@
-# -*- coding: utf-8 -*-
-
 """
 Cubic spline approximation (smoothing)
-
 """
 
-from csaps._version import __version__  # noqa
-
-from csaps._base import (
-    ISplinePPForm,
-    ISmoothingSpline,
-)
-from csaps._sspumv import (
-    SplinePPForm,
-    CubicSmoothingSpline,
-)
-from csaps._sspndg import (
-    NdGridSplinePPForm,
-    NdGridCubicSmoothingSpline,
-)
-from csaps._types import (
-    UnivariateDataType,
-    MultivariateDataType,
-    NdGridDataType,
-)
-from csaps._shortcut import csaps, AutoSmoothingResult
+from csaps._base import ISmoothingSpline, ISplinePPForm
+from csaps._shortcut import AutoSmoothingResult, csaps
+from csaps._sspndg import NdGridCubicSmoothingSpline, NdGridSplinePPForm
+from csaps._sspumv import CubicSmoothingSpline, SplinePPForm
+from csaps._types import MultivariateDataType, NdGridDataType, UnivariateDataType
+from csaps._version import __version__
 
 __all__ = [
     # Shortcut
@@ -43,4 +26,5 @@
     'UnivariateDataType',
     'MultivariateDataType',
     'NdGridDataType',
+    '__version__',
 ]
diff --git a/csaps/_base.py b/csaps/_base.py
index 38c14b3..d5ab3dc 100644
--- a/csaps/_base.py
+++ b/csaps/_base.py
@@ -1,16 +1,13 @@
-# -*- coding: utf-8 -*-
-
 """
 The base classes and interfaces
-
 """
 
+from typing import Generic, Optional, Tuple
 import abc
-from typing import Generic, Tuple, Optional
 
 import numpy as np
 
-from ._types import TData, TProps, TSmooth, TXi, TNu, TExtrapolate, TSpline
+from ._types import TData, TExtrapolate, TNu, TProps, TSmooth, TSpline, TXi
 
 
 class ISplinePPForm(abc.ABC, Generic[TData, TProps]):
@@ -76,7 +73,7 @@ def ndim(self) -> int:
 
     @property
     @abc.abstractmethod
-    def shape(self) -> Tuple[int]:
+    def shape(self) -> Tuple[int, ...]:
         """Returns the source data shape
 
         Returns
@@ -105,9 +102,6 @@ def spline(self) -> TSpline:
         """
 
     @abc.abstractmethod
-    def __call__(self,
-                 xi: TXi,
-                 nu: Optional[TNu] = None,
-                 extrapolate: Optional[TExtrapolate] = None) -> np.ndarray:
+    def __call__(self, xi: TXi, nu: Optional[TNu] = None, extrapolate: Optional[TExtrapolate] = None) -> np.ndarray:
         """Evaluates spline on the data sites
         """
diff --git a/csaps/_reshape.py b/csaps/_reshape.py
index 98e520c..4479680 100644
--- a/csaps/_reshape.py
+++ b/csaps/_reshape.py
@@ -1,9 +1,7 @@
-# -*- coding: utf-8 -*-
-
+import typing as ty
 import functools
-import operator
 from itertools import chain
-import typing as ty
+import operator
 
 import numpy as np
 from numpy.lib.stride_tricks import as_strided
@@ -152,13 +150,12 @@ def umv_coeffs_to_flatten(arr: np.ndarray):
         strides = arr.strides[:-3:-1]
         arr_view = as_strided(arr, shape=shape, strides=strides)
     else:  # pragma: no cover
-        raise ValueError(
-            f"The array ndim must be 2 or 3, but given array has ndim={arr.ndim}.")
+        raise ValueError(f"The array ndim must be 2 or 3, but given array has ndim={arr.ndim}.")
 
     return arr_view
 
 
-def ndg_coeffs_to_canonical(arr: np.ndarray, pieces: ty.Tuple[int]) -> np.ndarray:
+def ndg_coeffs_to_canonical(arr: np.ndarray, pieces: ty.Tuple[int, ...]) -> np.ndarray:
     """Returns array canonical view for given n-d grid coeffs flatten array
 
     Creates n-d array canonical view with shape (k0, ..., kn, p0, ..., pn) for given
diff --git a/csaps/_shortcut.py b/csaps/_shortcut.py
index 31e66db..56233b5 100644
--- a/csaps/_shortcut.py
+++ b/csaps/_shortcut.py
@@ -1,17 +1,14 @@
-# -*- coding: utf-8 -*-
-
 """
-The module provised `csaps` shortcut function for smoothing data
-
+The module provided `csaps` shortcut function for smoothing data
 """
 
+from typing import NamedTuple, Optional, Sequence, Union, overload
 from collections import abc as c_abc
-from typing import Optional, Union, Sequence, NamedTuple, overload
 
 from ._base import ISmoothingSpline
+from ._sspndg import NdGridCubicSmoothingSpline, ndgrid_prepare_data_vectors
 from ._sspumv import CubicSmoothingSpline
-from ._sspndg import ndgrid_prepare_data_vectors, NdGridCubicSmoothingSpline
-from ._types import UnivariateDataType, MultivariateDataType, NdGridDataType
+from ._types import MultivariateDataType, NdGridDataType, UnivariateDataType
 
 
 class AutoSmoothingResult(NamedTuple):
@@ -28,78 +25,100 @@ class AutoSmoothingResult(NamedTuple):
 # csaps signatures
 #
 @overload
-def csaps(xdata: UnivariateDataType,
-          ydata: MultivariateDataType,
-          *,
-          weights: Optional[UnivariateDataType] = None,
-          smooth: Optional[float] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> ISmoothingSpline: ...
+def csaps(
+    xdata: UnivariateDataType,
+    ydata: MultivariateDataType,
+    *,
+    weights: Optional[UnivariateDataType] = None,
+    smooth: Optional[float] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> ISmoothingSpline:  # pragma: no cover
+    ...
 
 
 @overload
-def csaps(xdata: UnivariateDataType,
-          ydata: MultivariateDataType,
-          xidata: UnivariateDataType,
-          *,
-          weights: Optional[UnivariateDataType] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> AutoSmoothingResult: ...
+def csaps(
+    xdata: UnivariateDataType,
+    ydata: MultivariateDataType,
+    xidata: UnivariateDataType,
+    *,
+    weights: Optional[UnivariateDataType] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> AutoSmoothingResult:  # pragma: no cover
+    ...
 
 
 @overload
-def csaps(xdata: UnivariateDataType,
-          ydata: MultivariateDataType,
-          xidata: UnivariateDataType,
-          *,
-          smooth: float,
-          weights: Optional[UnivariateDataType] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> MultivariateDataType: ...
+def csaps(
+    xdata: UnivariateDataType,
+    ydata: MultivariateDataType,
+    xidata: UnivariateDataType,
+    *,
+    smooth: float,
+    weights: Optional[UnivariateDataType] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> MultivariateDataType:  # pragma: no cover
+    ...
 
 
 @overload
-def csaps(xdata: NdGridDataType,
-          ydata: MultivariateDataType,
-          *,
-          weights: Optional[NdGridDataType] = None,
-          smooth: Optional[Sequence[float]] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> ISmoothingSpline: ...
+def csaps(
+    xdata: NdGridDataType,
+    ydata: MultivariateDataType,
+    *,
+    weights: Optional[NdGridDataType] = None,
+    smooth: Optional[Sequence[float]] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> ISmoothingSpline:  # pragma: no cover
+    ...
 
 
 @overload
-def csaps(xdata: NdGridDataType,
-          ydata: MultivariateDataType,
-          xidata: NdGridDataType,
-          *,
-          weights: Optional[NdGridDataType] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> AutoSmoothingResult: ...
+def csaps(
+    xdata: NdGridDataType,
+    ydata: MultivariateDataType,
+    xidata: NdGridDataType,
+    *,
+    weights: Optional[NdGridDataType] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> AutoSmoothingResult:  # pragma: no cover
+    ...
 
 
 @overload
-def csaps(xdata: NdGridDataType,
-          ydata: MultivariateDataType,
-          xidata: NdGridDataType,
-          *,
-          smooth: Sequence[float],
-          weights: Optional[NdGridDataType] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> MultivariateDataType: ...
+def csaps(
+    xdata: NdGridDataType,
+    ydata: MultivariateDataType,
+    xidata: NdGridDataType,
+    *,
+    smooth: Sequence[float],
+    weights: Optional[NdGridDataType] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> MultivariateDataType:  # pragma: no cover
+    ...
+
+
 #
 # csaps signatures
 # **************************************
 
 
-def csaps(xdata: Union[UnivariateDataType, NdGridDataType],
-          ydata: MultivariateDataType,
-          xidata: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
-          *,
-          weights: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
-          smooth: Optional[Union[float, Sequence[float]]] = None,
-          axis: Optional[int] = None,
-          normalizedsmooth: bool = False) -> Union[MultivariateDataType, ISmoothingSpline, AutoSmoothingResult]:
+def csaps(
+    xdata: Union[UnivariateDataType, NdGridDataType],
+    ydata: MultivariateDataType,
+    xidata: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
+    *,
+    weights: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
+    smooth: Optional[Union[float, Sequence[float]]] = None,
+    axis: Optional[int] = None,
+    normalizedsmooth: bool = False,
+) -> Union[MultivariateDataType, ISmoothingSpline, AutoSmoothingResult]:
     """Smooths the univariate/multivariate/gridded data or computes the corresponding splines
 
     This function might be used as the main API for smoothing any data.
@@ -213,10 +232,22 @@ def csaps(xdata: Union[UnivariateDataType, NdGridDataType],
 
     if umv:
         axis = -1 if axis is None else axis
-        sp = CubicSmoothingSpline(xdata, ydata, weights=weights, smooth=smooth, axis=axis,
-                                  normalizedsmooth=normalizedsmooth)
+        sp = CubicSmoothingSpline(
+            xdata,
+            ydata,
+            weights=weights,
+            smooth=smooth,
+            axis=axis,
+            normalizedsmooth=normalizedsmooth,
+        )
     else:
-        sp = NdGridCubicSmoothingSpline(xdata, ydata, weights, smooth, normalizedsmooth=normalizedsmooth)
+        sp = NdGridCubicSmoothingSpline(
+            xdata,
+            ydata,
+            weights,
+            smooth,
+            normalizedsmooth=normalizedsmooth,
+        )
 
     if xidata is None:
         return sp
diff --git a/csaps/_sspndg.py b/csaps/_sspndg.py
index 8143cb3..32ecc8f 100644
--- a/csaps/_sspndg.py
+++ b/csaps/_sspndg.py
@@ -1,27 +1,24 @@
-# -*- coding: utf-8 -*-
-
 """
 ND-Gridded cubic smoothing spline implementation
-
 """
 
+from typing import Optional, Sequence, Tuple, Union
 import collections.abc as c_abc
 from numbers import Number
-from typing import Tuple, Sequence, Optional, Union
 
 import numpy as np
-from scipy.interpolate import PPoly, NdPPoly
+from scipy.interpolate import NdPPoly, PPoly
 
-from ._base import ISplinePPForm, ISmoothingSpline
-from ._types import UnivariateDataType, NdGridDataType
-from ._sspumv import CubicSmoothingSpline
+from ._base import ISmoothingSpline, ISplinePPForm
 from ._reshape import (
+    ndg_coeffs_to_canonical,
+    ndg_coeffs_to_flatten,
     prod,
     umv_coeffs_to_canonical,
     umv_coeffs_to_flatten,
-    ndg_coeffs_to_canonical,
-    ndg_coeffs_to_flatten,
 )
+from ._sspumv import CubicSmoothingSpline
+from ._types import NdGridDataType, UnivariateDataType
 
 
 def ndgrid_prepare_data_vectors(data, name, min_size: int = 2) -> Tuple[np.ndarray, ...]:
@@ -41,8 +38,7 @@ def ndgrid_prepare_data_vectors(data, name, min_size: int = 2) -> Tuple[np.ndarr
     return tuple(data)
 
 
-class NdGridSplinePPForm(ISplinePPForm[Tuple[np.ndarray, ...], Tuple[int, ...]],
-                         NdPPoly):
+class NdGridSplinePPForm(ISplinePPForm[Tuple[np.ndarray, ...], Tuple[int, ...]], NdPPoly):
     """N-D grid spline representation in PP-form
 
     N-D grid spline is represented in piecewise tensor product polynomial form.
@@ -80,10 +76,12 @@ def ndim(self) -> int:
     def shape(self) -> Tuple[int, ...]:
         return tuple(len(xi) for xi in self.x)
 
-    def __call__(self,
-                 x: Sequence[UnivariateDataType],
-                 nu: Optional[Tuple[int, ...]] = None,
-                 extrapolate: Optional[bool] = None) -> np.ndarray:
+    def __call__(
+        self,
+        x: Sequence[UnivariateDataType],
+        nu: Optional[Tuple[int, ...]] = None,
+        extrapolate: Optional[bool] = None,
+    ) -> np.ndarray:
         """Evaluate the spline for given data
 
         Parameters
@@ -111,11 +109,10 @@ def __call__(self,
         x = ndgrid_prepare_data_vectors(x, 'x', min_size=1)
 
         if len(x) != self.ndim:
-            raise ValueError(
-                f"'x' sequence must have length {self.ndim} according to 'breaks'")
+            raise ValueError(f"'x' sequence must have length {self.ndim} according to 'breaks'")
 
         if nu is None:
-            nu = (0,) * len(x)
+            nu = (0, ) * len(x)
 
         if extrapolate is None:
             extrapolate = True
@@ -157,13 +154,15 @@ def __repr__(self):  # pragma: no cover
         )
 
 
-class NdGridCubicSmoothingSpline(ISmoothingSpline[
-                                     NdGridSplinePPForm,
-                                     Tuple[float, ...],
-                                     NdGridDataType,
-                                     Tuple[int, ...],
-                                     bool,
-                                 ]):
+class NdGridCubicSmoothingSpline(
+    ISmoothingSpline[
+        NdGridSplinePPForm,
+        Tuple[float, ...],
+        NdGridDataType,
+        Tuple[int, ...],
+        bool,
+    ]
+):
     """N-D grid cubic smoothing spline
 
     Class implements N-D grid data smoothing (piecewise tensor product polynomial).
@@ -203,12 +202,14 @@ class NdGridCubicSmoothingSpline(ISmoothingSpline[
 
     __module__ = 'csaps'
 
-    def __init__(self,
-                 xdata: NdGridDataType,
-                 ydata: np.ndarray,
-                 weights: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
-                 smooth: Optional[Union[float, Sequence[Optional[float]]]] = None,
-                 normalizedsmooth: bool = False) -> None:
+    def __init__(
+        self,
+        xdata: NdGridDataType,
+        ydata: np.ndarray,
+        weights: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
+        smooth: Optional[Union[float, Sequence[Optional[float]]]] = None,
+        normalizedsmooth: bool = False,
+    ) -> None:
 
         x, y, w, s = self._prepare_data(xdata, ydata, weights, smooth)
         coeffs, smooth = self._make_spline(x, y, w, s, normalizedsmooth)
@@ -216,10 +217,12 @@ def __init__(self,
         self._spline = NdGridSplinePPForm.construct_fast(coeffs, x)
         self._smooth = smooth
 
-    def __call__(self,
-                 x: Union[NdGridDataType, Sequence[Number]],
-                 nu: Optional[Tuple[int, ...]] = None,
-                 extrapolate: Optional[bool] = None) -> np.ndarray:
+    def __call__(
+        self,
+        x: Union[NdGridDataType, Sequence[Number]],
+        nu: Optional[Tuple[int, ...]] = None,
+        extrapolate: Optional[bool] = None,
+    ) -> np.ndarray:
         """Evaluate the spline for given data
 
         Parameters
@@ -274,13 +277,11 @@ def _prepare_data(cls, xdata, ydata, weights, smooth):
         data_ndim = len(xdata)
 
         if ydata.ndim != data_ndim:
-            raise ValueError(
-                f"'ydata' must have dimension {data_ndim} according to 'xdata'")
+            raise ValueError(f"'ydata' must have dimension {data_ndim} according to 'xdata'")
 
         for axis, (yd, xs) in enumerate(zip(ydata.shape, map(len, xdata))):
             if yd != xs:
-                raise ValueError(
-                    f"'ydata' ({yd}) and xdata ({xs}) sizes mismatch for axis {axis}")
+                raise ValueError(f"'ydata' ({yd}) and xdata ({xs}) sizes mismatch for axis {axis}")
 
         if not weights:
             weights = [None] * data_ndim
@@ -288,14 +289,12 @@ def _prepare_data(cls, xdata, ydata, weights, smooth):
             weights = ndgrid_prepare_data_vectors(weights, 'weights')
 
         if len(weights) != data_ndim:
-            raise ValueError(
-                f"'weights' ({len(weights)}) and 'xdata' ({data_ndim}) dimensions mismatch")
+            raise ValueError(f"'weights' ({len(weights)}) and 'xdata' ({data_ndim}) dimensions mismatch")
 
         for axis, (w, x) in enumerate(zip(weights, xdata)):
             if w is not None:
                 if w.size != x.size:
-                    raise ValueError(
-                        f"'weights' ({w.size}) and 'xdata' ({x.size}) sizes mismatch for axis {axis}")
+                    raise ValueError(f"'weights' ({w.size}) and 'xdata' ({x.size}) sizes mismatch for axis {axis}")
 
         if smooth is None:
             smooth = [None] * data_ndim
@@ -308,7 +307,8 @@ def _prepare_data(cls, xdata, ydata, weights, smooth):
         if len(smooth) != data_ndim:
             raise ValueError(
                 'Number of smoothing parameter values must '
-                f'be equal number of dimensions ({data_ndim})')
+                f'be equal number of dimensions ({data_ndim})'
+            )
 
         return xdata, ydata, weights, smooth
 
@@ -318,8 +318,13 @@ def _make_spline(xdata, ydata, weights, smooth, normalizedsmooth):
 
         if ndim == 1:
             s = CubicSmoothingSpline(
-                xdata[0], ydata, weights=weights[0], smooth=smooth[0], normalizedsmooth=normalizedsmooth)
-            return s.spline.coeffs, (s.smooth,)
+                xdata[0],
+                ydata,
+                weights=weights[0],
+                smooth=smooth[0],
+                normalizedsmooth=normalizedsmooth,
+            )
+            return s.spline.coeffs, (s.smooth, )
 
         shape = ydata.shape
         coeffs = ydata
@@ -334,7 +339,12 @@ def _make_spline(xdata, ydata, weights, smooth, normalizedsmooth):
                 coeffs = coeffs.reshape(prod(coeffs.shape[:-1]), coeffs.shape[-1])
 
             s = CubicSmoothingSpline(
-                xdata[i], coeffs, weights=weights[i], smooth=smooth[i], normalizedsmooth=normalizedsmooth)
+                xdata[i],
+                coeffs,
+                weights=weights[i],
+                smooth=smooth[i],
+                normalizedsmooth=normalizedsmooth,
+            )
 
             smooths.append(s.smooth)
             coeffs = umv_coeffs_to_flatten(s.spline.coeffs)
diff --git a/csaps/_sspumv.py b/csaps/_sspumv.py
index d7f3568..ba486aa 100644
--- a/csaps/_sspumv.py
+++ b/csaps/_sspumv.py
@@ -1,22 +1,18 @@
-# -*- coding: utf-8 -*-
-
 """
 Univariate/multivariate cubic smoothing spline implementation
-
 """
 
+from typing import List, Optional, Tuple, Union
 import functools
-from typing import Optional, Union, Tuple, List
 
 import numpy as np
-
+from scipy.interpolate import PPoly
 import scipy.sparse as sp
 import scipy.sparse.linalg as la
-from scipy.interpolate import PPoly
 
-from ._base import ISplinePPForm, ISmoothingSpline
-from ._types import UnivariateDataType, MultivariateDataType
-from ._reshape import to_2d, prod
+from ._base import ISmoothingSpline, ISplinePPForm
+from ._reshape import prod, to_2d
+from ._types import MultivariateDataType, UnivariateDataType
 
 
 class SplinePPForm(ISplinePPForm[np.ndarray, int], PPoly):
@@ -61,7 +57,7 @@ def ndim(self) -> int:
         return prod(shape)
 
     @property
-    def shape(self) -> Tuple[int]:
+    def shape(self) -> Tuple[int, ...]:
         """Returns the source data shape
         """
         shape: List[int] = list(self.c.shape[2:])
@@ -82,13 +78,7 @@ def __repr__(self):  # pragma: no cover
         )
 
 
-class CubicSmoothingSpline(ISmoothingSpline[
-                               SplinePPForm,
-                               float,
-                               UnivariateDataType,
-                               int,
-                               Union[bool, str]
-                           ]):
+class CubicSmoothingSpline(ISmoothingSpline[SplinePPForm, float, UnivariateDataType, int, Union[bool, str]]):
     """Cubic smoothing spline
 
     The cubic spline implementation for univariate/multivariate data.
@@ -125,13 +115,15 @@ class CubicSmoothingSpline(ISmoothingSpline[
 
     __module__ = 'csaps'
 
-    def __init__(self,
-                 xdata: UnivariateDataType,
-                 ydata: MultivariateDataType,
-                 weights: Optional[UnivariateDataType] = None,
-                 smooth: Optional[float] = None,
-                 axis: int = -1,
-                 normalizedsmooth: bool = False):
+    def __init__(
+        self,
+        xdata: UnivariateDataType,
+        ydata: MultivariateDataType,
+        weights: Optional[UnivariateDataType] = None,
+        smooth: Optional[float] = None,
+        axis: int = -1,
+        normalizedsmooth: bool = False,
+    ):
 
         x, y, w, shape, axis = self._prepare_data(xdata, ydata, weights, axis)
         coeffs, smooth = self._make_spline(x, y, w, smooth, shape, normalizedsmooth)
@@ -140,10 +132,12 @@ def __init__(self,
         self._smooth = smooth
         self._spline = spline
 
-    def __call__(self,
-                 x: UnivariateDataType,
-                 nu: Optional[int] = None,
-                 extrapolate: Optional[Union[bool, str]] = None) -> np.ndarray:
+    def __call__(
+        self,
+        x: UnivariateDataType,
+        nu: Optional[int] = None,
+        extrapolate: Optional[Union[bool, str]] = None,
+    ) -> np.ndarray:
         """Evaluate the spline for given data
 
         Parameters
@@ -211,7 +205,8 @@ def _prepare_data(xdata, ydata, weights, axis):
         if ydata.shape[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"that is equal to 'xdata' size ({xdata.size})"
+            )
 
         # Rolling axis for using its shape while constructing coeffs array
         shape = np.rollaxis(ydata, axis).shape
@@ -251,9 +246,9 @@ def _normalize_smooth(x: np.ndarray, w: np.ndarray, smooth: Optional[float]):
 
         span = np.ptp(x)
 
-        eff_x = 1 + (span ** 2) / np.sum(np.diff(x) ** 2)
-        eff_w = np.sum(w) ** 2 / np.sum(w ** 2)
-        k = 80 * (span ** 3) * (x.size ** -2) * (eff_x ** -0.5) * (eff_w ** -0.5)
+        eff_x = 1 + (span**2) / np.sum(np.diff(x)**2)
+        eff_w = np.sum(w)**2 / np.sum(w**2)
+        k = 80 * (span**3) * (x.size**-2) * (eff_x**-0.5) * (eff_w**-0.5)
 
         s = 0.5 if smooth is None else smooth
         p = s / (s + (1 - s) * k)
@@ -266,8 +261,7 @@ def _make_spline(x, y, w, smooth, shape, normalizedsmooth):
         dx = np.diff(x)
 
         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(y, axis=1)
         dy_dx = dy / dx
@@ -291,8 +285,7 @@ def _make_spline(x, y, w, smooth, shape, normalizedsmooth):
         diags_qtw = np.vstack((dx_recip[:-1], -(dx_recip[1:] + dx_recip[:-1]), dx_recip[1:]))
         diags_sqrw_recip = 1. / np.sqrt(w)
 
-        qtw = (sp.diags(diags_qtw, [0, 1, 2], (pcount - 2, pcount)) @
-               sp.diags(diags_sqrw_recip, 0, (pcount, pcount)))
+        qtw = (sp.diags(diags_qtw, [0, 1, 2], (pcount - 2, pcount)) @ sp.diags(diags_sqrw_recip, 0, (pcount, pcount)))
         qtw = qtw @ qtw.T
 
         p = smooth
diff --git a/csaps/_types.py b/csaps/_types.py
index 2aec8fc..810ce13 100644
--- a/csaps/_types.py
+++ b/csaps/_types.py
@@ -1,19 +1,17 @@
-# -*- coding: utf-8 -*-
-
 """
 Type-hints and type vars
-
 """
 
+from typing import Sequence, Tuple, TypeVar, Union
 from collections import abc
-from typing import Union, Sequence, Tuple, TypeVar
 from numbers import Number
-import numpy as np
 
+from typing_extensions import TypeAlias
+import numpy as np
 
-UnivariateDataType = Union[np.ndarray, Sequence[Number]]
-MultivariateDataType = Union[np.ndarray, abc.Sequence]
-NdGridDataType = Sequence[UnivariateDataType]
+UnivariateDataType: TypeAlias = Union[np.ndarray, Sequence[Number]]
+MultivariateDataType: TypeAlias = Union[np.ndarray, abc.Sequence]
+NdGridDataType: TypeAlias = Sequence[UnivariateDataType]
 
 TData = TypeVar('TData', np.ndarray, Sequence[np.ndarray])
 TProps = TypeVar('TProps', int, Tuple[int, ...])
diff --git a/docs/conf.py b/docs/conf.py
index e985ae5..393e297 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -16,25 +16,21 @@
 ROOT_DIR = pathlib.Path(__file__).parent.parent
 sys.path.insert(0, str(ROOT_DIR))
 
-
 # -- Project information -----------------------------------------------------
 
 project = 'csaps'
-copyright = '2017-2020, Eugene Prilepin'
+copyright = '2017-2023, Eugene Prilepin'  # noqa
 author = 'Eugene Prilepin'
 
 
 def _get_version():
-    about = {}
-    ver_mod = ROOT_DIR / 'csaps' / '_version.py'
-    exec(ver_mod.read_text(), about)
-    return about['__version__']
+    from csaps import __version__
+    return __version__
 
 
 # The full version, including alpha/beta/rc tags
 release = _get_version()
 
-
 # -- General configuration ---------------------------------------------------
 
 # Add any Sphinx extension module names here, as strings. They can be
@@ -50,9 +46,7 @@ def _get_version():
     'm2r2',
 ]
 
-intersphinx_mapping = {
-    'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None)
-}
+intersphinx_mapping = {'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None)}
 
 # Extension settings
 plot_apply_rcparams = True
@@ -111,13 +105,10 @@ def univariate_data(n=25, seed=1234):
 html_theme_options = {
     'fixed_sidebar': 'true',
     'show_powered_by': 'false',
-
     'description': 'Cubic spline approximation (smoothing)',
-
     'github_user': 'espdev',
     'github_repo': 'csaps',
     'github_type': 'star',
-
     'extra_nav_links': {
         'GitHub repository': 'https://github.com/espdev/csaps',
         'PyPI': 'https://pypi.org/project/csaps',
diff --git a/poetry.lock b/poetry.lock
index 0c3bf5a..5c7697c 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -465,6 +465,20 @@ files = [
     {file = "iniconfig-2.0.0.tar.gz", hash = "sha256:2d91e135bf72d31a410b17c16da610a82cb55f6b0477d1a902134b24a455b8b3"},
 ]
 
+[[package]]
+name = "isort"
+version = "5.13.2"
+description = "A Python utility / library to sort Python imports."
+optional = false
+python-versions = ">=3.8.0"
+files = [
+    {file = "isort-5.13.2-py3-none-any.whl", hash = "sha256:8ca5e72a8d85860d5a3fa69b8745237f2939afe12dbf656afbcb47fe72d947a6"},
+    {file = "isort-5.13.2.tar.gz", hash = "sha256:48fdfcb9face5d58a4f6dde2e72a1fb8dcaf8ab26f95ab49fab84c2ddefb0109"},
+]
+
+[package.extras]
+colors = ["colorama (>=0.4.6)"]
+
 [[package]]
 name = "jinja2"
 version = "3.1.2"
@@ -988,6 +1002,21 @@ files = [
 docs = ["furo", "olefile", "sphinx (>=2.4)", "sphinx-copybutton", "sphinx-inline-tabs", "sphinx-removed-in", "sphinxext-opengraph"]
 tests = ["check-manifest", "coverage", "defusedxml", "markdown2", "olefile", "packaging", "pyroma", "pytest", "pytest-cov", "pytest-timeout"]
 
+[[package]]
+name = "platformdirs"
+version = "4.1.0"
+description = "A small Python package for determining appropriate platform-specific dirs, e.g. a \"user data dir\"."
+optional = false
+python-versions = ">=3.8"
+files = [
+    {file = "platformdirs-4.1.0-py3-none-any.whl", hash = "sha256:11c8f37bcca40db96d8144522d925583bdb7a31f7b0e37e3ed4318400a8e2380"},
+    {file = "platformdirs-4.1.0.tar.gz", hash = "sha256:906d548203468492d432bcb294d4bc2fff751bf84971fbb2c10918cc206ee420"},
+]
+
+[package.extras]
+docs = ["furo (>=2023.7.26)", "proselint (>=0.13)", "sphinx (>=7.1.1)", "sphinx-autodoc-typehints (>=1.24)"]
+test = ["appdirs (==1.4.4)", "covdefaults (>=2.3)", "pytest (>=7.4)", "pytest-cov (>=4.1)", "pytest-mock (>=3.11.1)"]
+
 [[package]]
 name = "pluggy"
 version = "1.3.0"
@@ -1438,6 +1467,22 @@ brotli = ["brotli (>=1.0.9)", "brotlicffi (>=0.8.0)"]
 socks = ["pysocks (>=1.5.6,!=1.5.7,<2.0)"]
 zstd = ["zstandard (>=0.18.0)"]
 
+[[package]]
+name = "yapf"
+version = "0.40.2"
+description = "A formatter for Python code"
+optional = false
+python-versions = ">=3.7"
+files = [
+    {file = "yapf-0.40.2-py3-none-any.whl", hash = "sha256:adc8b5dd02c0143108878c499284205adb258aad6db6634e5b869e7ee2bd548b"},
+    {file = "yapf-0.40.2.tar.gz", hash = "sha256:4dab8a5ed7134e26d57c1647c7483afb3f136878b579062b786c9ba16b94637b"},
+]
+
+[package.dependencies]
+importlib-metadata = ">=6.6.0"
+platformdirs = ">=3.5.1"
+tomli = ">=2.0.1"
+
 [[package]]
 name = "zipp"
 version = "3.17.0"
@@ -1456,4 +1501,4 @@ testing = ["big-O", "jaraco.functools", "jaraco.itertools", "more-itertools", "p
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.8"
-content-hash = "81677aa49cd174f202b0307e4a092938f181bdbe7c0271fa9748918f3fbe5a50"
+content-hash = "78eb0c54af951e30c3aa351860b953eb70270bd6eb800a7e79c4ca7e143e9d1d"
diff --git a/pyproject.toml b/pyproject.toml
index 4419824..16ad0a1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -44,6 +44,7 @@ packages = [
 
 [tool.poetry.dependencies]
 python = ">=3.8"
+typing-extensions = "*"
 numpy = [
     {version = "*", python = "<3.12"},
     {version = ">=1.26.2", python = ">=3.12"},
@@ -63,6 +64,8 @@ sphinx = "^7.1.2"
 numpydoc = "^1.6.0"
 matplotlib = "^3.7.4"
 ruff = "^0.1.8"
+isort = "^5.13.2"
+yapf = "^0.40.2"
 mypy = "^1.7.1"
 
 
@@ -82,15 +85,41 @@ select = [
     "T201",    # print found
     "T203",    # pprint found
 ]
-ignore = [
-    "A001",
-]
+ignore = []
 exclude = [
     ".ruff_cache",
     ".venv",
 ]
 
 
+[tool.isort]
+multi_line_output = 3
+include_trailing_comma = true
+force_sort_within_sections = true
+force_to_top = ["typing", "typing_extensions", "pytest"]
+line_length = 120
+skip = [
+    ".venv",
+]
+
+
+[tool.yapf]
+based_on_style = "pep8"
+column_limit = 120
+blank_line_before_module_docstring = true
+coalesce_brackets = false
+dedent_closing_brackets = true
+split_arguments_when_comma_terminated = true
+split_before_dict_set_generator = true
+allow_split_before_dict_value = false
+
+
+[tool.yapfignore]
+ignore_patterns = [
+    ".venv/**/*",
+]
+
+
 [tool.mypy]
 python_version = "3.8"
 
diff --git a/tests/conftest.py b/tests/conftest.py
index 436aed4..43cd7fe 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,5 +1,3 @@
-# -*- coding: utf-8 -*-
-
 from typing import NamedTuple, Tuple
 
 import pytest
@@ -23,6 +21,8 @@ def univariate_data() -> UnivariateData:
     """Univariate exponential noisy data sample
     """
 
+    # yapf: disable
+
     x = np.linspace(-5., 5., 25)
     y = np.array([
         0.015771474002402, 0.161329316958106, 0.133494845724251, 0.281006799675995, 0.343006057841707,
@@ -125,6 +125,8 @@ def univariate_data() -> UnivariateData:
         5.401286443327005, 5.415015799831229, 5.432185121115394, 5.453047420840733,
     ])
 
+    # yapf: enable
+
     integral = 5.453047420840733
     smooth = 0.992026535689226
 
@@ -147,66 +149,111 @@ def ndgrid_2d_data() -> NdGrid2dData:
     xy = (np.linspace(-3.0, 3.0, 5), np.linspace(-3.5, 3.5, 7))
 
     z = np.array(
-        [[-0.153527183790132, 0.360477327564227, -0.400800187993851, -0.661751768834967,
-          1.39827150034968, 1.044246054228617, 0.069681364921588],
-         [0.211217178485871, 0.592683030706752, 1.294599451385471, -4.924883983709012,
-          -2.641771353280953, 0.245330967159293, 0.171928943618129],
-         [1.012132959440344, 0.132792505223302, -3.970096642307903, 0.702129940268655,
-          4.729521910567126, 0.208213433055832, -0.40275495492284],
-         [0.35749708646856, 2.409904780478664, 0.892801916808247, 7.563804764350773,
-          2.510824654279176, 0.317303593544217, 0.393080231785911],
-         [0.000706314884567, 1.009080744382149, -0.45874273220015, -0.323494125914201,
-          -1.700362064179427, -1.394209767885332, -0.645566364768713]
-         ])
+        [
+            [
+                -0.153527183790132, 0.360477327564227, -0.400800187993851, -0.661751768834967, 1.39827150034968,
+                1.044246054228617, 0.069681364921588
+            ],
+            [
+                0.211217178485871, 0.592683030706752, 1.294599451385471, -4.924883983709012, -2.641771353280953,
+                0.245330967159293, 0.171928943618129
+            ],
+            [
+                1.012132959440344, 0.132792505223302, -3.970096642307903, 0.702129940268655, 4.729521910567126,
+                0.208213433055832, -0.40275495492284
+            ],
+            [
+                0.35749708646856, 2.409904780478664, 0.892801916808247, 7.563804764350773, 2.510824654279176,
+                0.317303593544217, 0.393080231785911
+            ],
+            [
+                0.000706314884567, 1.009080744382149, -0.45874273220015, -0.323494125914201, -1.700362064179427,
+                -1.394209767885332, -0.645566364768713
+            ]
+        ]
+    )
 
     zi = np.array(
-        [[-0.055377680470166, 0.195656616225213, -0.295030253111251, -0.830533929888634,
-          0.193176060095987, 0.770374649757329, 0.252865339751650],
-         [0.471994652733459, 0.293417006151304, -0.758106516247562, -1.960431309380293,
-          -0.781936045165379, 0.216341632490716, 0.180333235920312],
-         [0.875919224697303, 0.067344259041702, -0.735889940425535, 0.882313890047783,
-          2.056305063365266, 0.896850201038262, -0.190314083560006],
-         [0.606245376082951, 0.941947682137626, 1.225331206624579, 3.379540894700002,
-          2.581257432070516, 0.581783850872262, -0.187728390603794],
-         [0.183397630824828, 0.805748594104382, 0.503605325241657, 0.264260721868410,
-          -0.874052860773297, -1.188420383689933, -0.617919628357980]
-         ])
+        [
+            [
+                -0.055377680470166, 0.195656616225213, -0.295030253111251, -0.830533929888634, 0.193176060095987,
+                0.770374649757329, 0.252865339751650
+            ],
+            [
+                0.471994652733459, 0.293417006151304, -0.758106516247562, -1.960431309380293, -0.781936045165379,
+                0.216341632490716, 0.180333235920312
+            ],
+            [
+                0.875919224697303, 0.067344259041702, -0.735889940425535, 0.882313890047783, 2.056305063365266,
+                0.896850201038262, -0.190314083560006
+            ],
+            [
+                0.606245376082951, 0.941947682137626, 1.225331206624579, 3.379540894700002, 2.581257432070516,
+                0.581783850872262, -0.187728390603794
+            ],
+            [
+                0.183397630824828, 0.805748594104382, 0.503605325241657, 0.264260721868410, -0.874052860773297,
+                -1.188420383689933, -0.617919628357980
+            ]
+        ]
+    )
 
     zi_d1 = np.array(
-        [[-0.121812472223565, -0.157326929297246, -0.828962829944406, -0.640465943383180,
-          0.503884652683063, 0.662549465069741, 0.501482011524879],
-         [-0.470497745190257, -0.466533613272337, 0.403062598673842, 0.767903075058928,
-          -0.334118185393139, -0.545086449738045, -0.154825057274055],
-         [0.194607321700336, 0.192434581435788, 1.293210437573223, 0.562712202846881,
-          -1.370259457802591, -0.823157992133989, -0.176531696944242],
-         [0.607431089728563, 0.510938352696447, -0.237734672001804, -1.058160086418725,
-          -0.269494123744448, 0.632430110330214, 0.577443254979813],
-         [0.006464517899505, -0.138732754726110, -1.426802191857966, -1.065472485030174,
-          1.029409644646945, 1.158453840382574, 0.696817777775121]
-         ])
+        [
+            [
+                -0.121812472223565, -0.157326929297246, -0.828962829944406, -0.640465943383180, 0.503884652683063,
+                0.662549465069741, 0.501482011524879
+            ],
+            [
+                -0.470497745190257, -0.466533613272337, 0.403062598673842, 0.767903075058928, -0.334118185393139,
+                -0.545086449738045, -0.154825057274055
+            ],
+            [
+                0.194607321700336, 0.192434581435788, 1.293210437573223, 0.562712202846881, -1.370259457802591,
+                -0.823157992133989, -0.176531696944242
+            ],
+            [
+                0.607431089728563, 0.510938352696447, -0.237734672001804, -1.058160086418725, -0.269494123744448,
+                0.632430110330214, 0.577443254979813
+            ],
+            [
+                0.006464517899505, -0.138732754726110, -1.426802191857966, -1.065472485030174, 1.029409644646945,
+                1.158453840382574, 0.696817777775121
+            ]
+        ]
+    )
 
     zi_d2 = np.array(
-        [[0., 0., 0., 0., 0., 0., 0.],
-         [0., 0.090236774837945, 3.432579482518258, -3.029508420063717,
-          -2.105055823406707, 1.260180219448802, 0.],
-         [0., -0.104263911255016, -2.890141730806884, -0.016560671330926,
-          3.251809714350141, -0.674203298932788, 0.],
-         [0., -0.111324652785139, -1.121581432777390, 3.822735995622450,
-          -0.837099042473098, -0.929483902301833, 0.],
-         [0., 0., 0., 0., 0., 0., 0.]
-         ])
+        [
+            [0., 0., 0., 0., 0., 0., 0.],
+            [0., 0.090236774837945, 3.432579482518258, -3.029508420063717, -2.105055823406707, 1.260180219448802, 0.],
+            [0., -0.104263911255016, -2.890141730806884, -0.016560671330926, 3.251809714350141, -0.674203298932788, 0.],
+            [0., -0.111324652785139, -1.121581432777390, 3.822735995622450, -0.837099042473098, -0.929483902301833, 0.],
+            [0., 0., 0., 0., 0., 0., 0.]
+        ]
+    )
 
     zi_ad1 = np.array(
-        [[0., 0., 0., 0., 0., 0., 0.],
-         [0., 0.509395864624591, 0.476401125691113, -1.598603993938397,
-          -3.760699480840437, -3.750794105947393, -3.076469947082991],
-         [0., 1.318187655218603, 0.563946491827388, -3.255188072198998,
-          -5.625801866969895, -4.504932326883326, -3.364137445816383],
-         [0., 2.511070517473814, 2.161223166080562, 0.664448782633726,
-          3.261109214575526, 7.566477396430049, 9.182402413317154],
-         [0., 3.799550963870154, 5.127806810771760, 6.567161731071815,
-          12.269550101243563, 17.158134933110624, 18.133537415972746]
-         ])
+        [
+            [0., 0., 0., 0., 0., 0., 0.],
+            [
+                0., 0.509395864624591, 0.476401125691113, -1.598603993938397, -3.760699480840437, -3.750794105947393,
+                -3.076469947082991
+            ],
+            [
+                0., 1.318187655218603, 0.563946491827388, -3.255188072198998, -5.625801866969895, -4.504932326883326,
+                -3.364137445816383
+            ],
+            [
+                0., 2.511070517473814, 2.161223166080562, 0.664448782633726, 3.261109214575526, 7.566477396430049,
+                9.182402413317154
+            ],
+            [
+                0., 3.799550963870154, 5.127806810771760, 6.567161731071815, 12.269550101243563, 17.158134933110624,
+                18.133537415972746
+            ]
+        ]
+    )
 
     integral = 18.133537415972746
     smooth = (0.727272727272727, 0.850021862702230)
@@ -228,9 +275,10 @@ def surface_data() -> SurfaceData:
     xy = (np.linspace(-3.0, 3.0, 61), np.linspace(-3.5, 3.5, 51))
     i, j = np.meshgrid(*xy, indexing='ij')
 
-    z = (3 * (1 - j) ** 2. * np.exp(-(j ** 2) - (i + 1) ** 2)
-         - 10 * (j / 5 - j ** 3 - i ** 5) * np.exp(-j ** 2 - i ** 2)
-         - 1 / 3 * np.exp(-(j + 1) ** 2 - i ** 2))
+    z = (
+        3 * (1 - j)**2. * np.exp(-(j**2) - (i + 1)**2) - 10 * (j / 5 - j**3 - i**5) * np.exp(-j**2 - i**2) -
+        1 / 3 * np.exp(-(j + 1)**2 - i**2)
+    )
 
     z += (np.random.randn(*z.shape) * 0.75)
 
diff --git a/tests/test_analysis.py b/tests/test_analysis.py
index 7663b4d..c4c79cd 100644
--- a/tests/test_analysis.py
+++ b/tests/test_analysis.py
@@ -1,6 +1,5 @@
-# -*- coding: utf-8 -*-
-
 import pytest
+
 import csaps
 
 
diff --git a/tests/test_csaps.py b/tests/test_csaps.py
index 3fbb4df..4808056 100644
--- a/tests/test_csaps.py
+++ b/tests/test_csaps.py
@@ -1,9 +1,7 @@
-# -*- coding: utf-8 -*-
-
 import pytest
 import numpy as np
 
-from csaps import csaps, AutoSmoothingResult, CubicSmoothingSpline, NdGridCubicSmoothingSpline
+from csaps import AutoSmoothingResult, CubicSmoothingSpline, NdGridCubicSmoothingSpline, csaps
 
 
 @pytest.fixture(scope='module')
@@ -11,7 +9,8 @@ def curve():
     np.random.seed(12345)
 
     x = np.linspace(-5., 5., 25)
-    y = np.exp(-(x / 2.5) ** 2) + (np.random.rand(25) - 0.2) * 0.3
+    y = np.exp(-(x / 2.5)**2) + (np.random.rand(25) - 0.2) * 0.3
+
     return x, y
 
 
@@ -22,10 +21,12 @@ def surface():
     x = [np.linspace(-3, 3, 61), np.linspace(-3.5, 3.5, 51)]
     i, j = np.meshgrid(*x, indexing='ij')
 
-    y = (3 * (1 - j) ** 2. * np.exp(-(j ** 2) - (i + 1) ** 2)
-         - 10 * (j / 5 - j ** 3 - i ** 5) * np.exp(-j ** 2 - i ** 2)
-         - 1 / 3 * np.exp(-(j + 1) ** 2 - i ** 2))
+    y = (
+        3 * (1 - j)**2. * np.exp(-(j**2) - (i + 1)**2) - 10 * (j / 5 - j**3 - i**5) * np.exp(-j**2 - i**2) -
+        1 / 3 * np.exp(-(j + 1)**2 - i**2)
+    )
     y += np.random.randn(*y.shape) * 0.75
+
     return x, y
 
 
@@ -74,14 +75,16 @@ def test_shortcut_output(data, tolist):
     assert isinstance(sp, sp_cls)
 
 
-@pytest.mark.parametrize('smooth, cls', [
-    (0.85, np.ndarray),
-    ([0.85, 0.85], np.ndarray),
-    (None, AutoSmoothingResult),
-    ([None, 0.85], AutoSmoothingResult),
-    ([0.85, None], AutoSmoothingResult),
-    ([None, None], AutoSmoothingResult),
-])
+@pytest.mark.parametrize(
+    'smooth, cls', [
+        (0.85, np.ndarray),
+        ([0.85, 0.85], np.ndarray),
+        (None, AutoSmoothingResult),
+        ([None, 0.85], AutoSmoothingResult),
+        ([0.85, None], AutoSmoothingResult),
+        ([None, None], AutoSmoothingResult),
+    ]
+)
 def test_shortcut_ndgrid_smooth_output(surface, smooth, cls):
     x, y = surface
 
@@ -96,12 +99,12 @@ def test_normalized_smooth(data, smooth, scale):
     x, y, xi, *_ = data
 
     x2 = (
-        [scale * np.array(xx, dtype=np.float64) for xx in x]
-        if isinstance(x, list) else scale * np.array(x, dtype=np.float64)
+        [scale * np.array(xx, dtype=np.float64)
+         for xx in x] if isinstance(x, list) else scale * np.array(x, dtype=np.float64)
     )
     xi2 = (
-        [scale * np.array(xx, dtype=np.float64) for xx in xi]
-        if isinstance(x, list) else scale * np.array(xi, dtype=np.float64)
+        [scale * np.array(xx, dtype=np.float64)
+         for xx in xi] if isinstance(x, list) else scale * np.array(xi, dtype=np.float64)
     )
 
     smoothed_data_a = csaps(x, y, xi, smooth=smooth, normalizedsmooth=True)
diff --git a/tests/test_ndg.py b/tests/test_ndg.py
index c1975cf..0fab16c 100644
--- a/tests/test_ndg.py
+++ b/tests/test_ndg.py
@@ -1,35 +1,37 @@
-# -*- coding: utf-8 -*-
-
 import pytest
-
 import numpy as np
 from scipy.interpolate import NdPPoly
+
 import csaps
 
 
-@pytest.mark.parametrize('x,y,w,p', [
-    ([1, 2, 3], np.ones((10, 10)), None, None),
-    ([[1], [1]], np.ones((1, 1)), None, None),
-    ([[1, 2, 3], [1, 2, 3]], np.ones((4, 3)), None, None),
-    ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3, 3)), None, None),
-    ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), [1, 2, 3], None),
-    ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), [[1, 2, 3]], None),
-    ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), [[1, 2], [1, 2]], None),
-    ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), None, [0.5, 0.4, 0.2]),
-    (np.array([[1, 2, 3], [4, 5, 6]]), np.ones((3, 3)), None, None),
-    ([np.arange(6).reshape(2, 3), np.arange(6).reshape(2, 3)], np.ones((6, 6)), None, None),
-])
+@pytest.mark.parametrize(
+    'x,y,w,p', [
+        ([1, 2, 3], np.ones((10, 10)), None, None),
+        ([[1], [1]], np.ones((1, 1)), None, None),
+        ([[1, 2, 3], [1, 2, 3]], np.ones((4, 3)), None, None),
+        ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3, 3)), None, None),
+        ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), [1, 2, 3], None),
+        ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), [[1, 2, 3]], None),
+        ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), [[1, 2], [1, 2]], None),
+        ([[1, 2, 3], [1, 2, 3]], np.ones((3, 3)), None, [0.5, 0.4, 0.2]),
+        (np.array([[1, 2, 3], [4, 5, 6]]), np.ones((3, 3)), None, None),
+        ([np.arange(6).reshape(2, 3), np.arange(6).reshape(2, 3)], np.ones((6, 6)), None, None),
+    ]
+)
 def test_invalid_data(x, y, w, p):
     with pytest.raises((ValueError, TypeError)):
         csaps.NdGridCubicSmoothingSpline(x, y, w, p)
 
 
-@pytest.mark.parametrize('x, xi', [
-    (([1, 2, 3],), ([],)),
-    (([1, 2, 3], [1, 2, 3]), ([1, 2], [])),
-    (([1, 2, 3], [1, 2, 3], [1, 2, 3]), ([1, 2, 3], [1, 2, 3])),
-    (([1, 2, 3], [1, 2, 3]), ([1, 2, 3], [1, 2, 3], [1, 2, 3])),
-])
+@pytest.mark.parametrize(
+    'x, xi', [
+        (([1, 2, 3], ), ([], )),
+        (([1, 2, 3], [1, 2, 3]), ([1, 2], [])),
+        (([1, 2, 3], [1, 2, 3], [1, 2, 3]), ([1, 2, 3], [1, 2, 3])),
+        (([1, 2, 3], [1, 2, 3]), ([1, 2, 3], [1, 2, 3], [1, 2, 3])),
+    ]
+)
 def test_invalid_evaluate_data(x, xi):
     np.random.seed(1234)
 
@@ -40,18 +42,20 @@ def test_invalid_evaluate_data(x, xi):
         s(xi)
 
 
-@pytest.mark.parametrize('shape, coeffs_shape, order, pieces, ndim', [
-    ((2,), (2, 1), (2,), (1,), 1),
-    ((3,), (4, 2), (4,), (2,), 1),
-    ((4,), (4, 3), (4,), (3,), 1),
-    ((2, 4), (2, 4, 1, 3), (2, 4), (1, 3), 2),
-    ((3, 4), (4, 4, 2, 3), (4, 4), (2, 3), 2),
-    ((4, 4), (4, 4, 3, 3), (4, 4), (3, 3), 2),
-    ((2, 3, 4), (2, 4, 4, 1, 2, 3), (2, 4, 4), (1, 2, 3), 3),
-    ((3, 4, 5), (4, 4, 4, 2, 3, 4), (4, 4, 4), (2, 3, 4), 3),
-    ((2, 3, 2, 6), (2, 4, 2, 4, 1, 2, 1, 5), (2, 4, 2, 4), (1, 2, 1, 5), 4),
-    ((3, 4, 5, 6), (4, 4, 4, 4, 2, 3, 4, 5), (4, 4, 4, 4), (2, 3, 4, 5), 4),
-])
+@pytest.mark.parametrize(
+    'shape, coeffs_shape, order, pieces, ndim', [
+        ((2, ), (2, 1), (2, ), (1, ), 1),
+        ((3, ), (4, 2), (4, ), (2, ), 1),
+        ((4, ), (4, 3), (4, ), (3, ), 1),
+        ((2, 4), (2, 4, 1, 3), (2, 4), (1, 3), 2),
+        ((3, 4), (4, 4, 2, 3), (4, 4), (2, 3), 2),
+        ((4, 4), (4, 4, 3, 3), (4, 4), (3, 3), 2),
+        ((2, 3, 4), (2, 4, 4, 1, 2, 3), (2, 4, 4), (1, 2, 3), 3),
+        ((3, 4, 5), (4, 4, 4, 2, 3, 4), (4, 4, 4), (2, 3, 4), 3),
+        ((2, 3, 2, 6), (2, 4, 2, 4, 1, 2, 1, 5), (2, 4, 2, 4), (1, 2, 1, 5), 4),
+        ((3, 4, 5, 6), (4, 4, 4, 4, 2, 3, 4, 5), (4, 4, 4, 4), (2, 3, 4, 5), 4),
+    ]
+)
 def test_ndsplineppform(shape, coeffs_shape, order, pieces, ndim):
     x = tuple(np.arange(s) for s in shape)
     y = np.arange(float(np.prod(shape))).reshape(shape)
@@ -78,46 +82,44 @@ def test_surface(surface_data):
     assert noisy_s.shape == ydata.shape
 
 
-@pytest.mark.parametrize('shape, smooth', [
-    ((4,), 0.0),
-    ((4,), 0.5),
-    ((4,), 1.0),
-
-    ((4,), (0.0,)),
-    ((4,), (0.5,)),
-    ((4,), (1.0,)),
-
-    ((4, 5), 0.0),
-    ((4, 5), 0.5),
-    ((4, 5), 1.0),
-    ((4, 5), (0.0, 0.0)),
-    ((4, 5), (0.0, 0.5)),
-    ((4, 5), (0.5, 0.0)),
-    ((4, 5), (0.5, 0.7)),
-    ((4, 5), (1.0, 0.0)),
-    ((4, 5), (0.0, 1.0)),
-    ((4, 5), (1.0, 1.0)),
-
-    ((4, 5, 6), 0.0),
-    ((4, 5, 6), 0.5),
-    ((4, 5, 6), 1.0),
-    ((4, 5, 6), (0.0, 0.0, 0.0)),
-    ((4, 5, 6), (0.5, 0.0, 0.0)),
-    ((4, 5, 6), (0.5, 0.6, 0.0)),
-    ((4, 5, 6), (0.0, 0.5, 0.6)),
-    ((4, 5, 6), (0.4, 0.5, 0.6)),
-
-    ((4, 5, 6, 7), 0.0),
-    ((4, 5, 6, 7), 0.5),
-    ((4, 5, 6, 7), 1.0),
-    ((4, 5, 6, 7), (0.0, 0.0, 0.0, 0.0)),
-    ((4, 5, 6, 7), (0.5, 0.0, 0.0, 0.0)),
-    ((4, 5, 6, 7), (0.0, 0.5, 0.0, 0.0)),
-    ((4, 5, 6, 7), (0.5, 0.6, 0.0, 0.0)),
-    ((4, 5, 6, 7), (0.0, 0.5, 0.6, 0.0)),
-    ((4, 5, 6, 7), (0.0, 0.5, 0.6, 0.7)),
-    ((4, 5, 6, 7), (0.4, 0.5, 0.6, 0.7)),
-])
+@pytest.mark.parametrize(
+    'shape, smooth', [
+        ((4, ), 0.0),
+        ((4, ), 0.5),
+        ((4, ), 1.0),
+        ((4, ), (0.0, )),
+        ((4, ), (0.5, )),
+        ((4, ), (1.0, )),
+        ((4, 5), 0.0),
+        ((4, 5), 0.5),
+        ((4, 5), 1.0),
+        ((4, 5), (0.0, 0.0)),
+        ((4, 5), (0.0, 0.5)),
+        ((4, 5), (0.5, 0.0)),
+        ((4, 5), (0.5, 0.7)),
+        ((4, 5), (1.0, 0.0)),
+        ((4, 5), (0.0, 1.0)),
+        ((4, 5), (1.0, 1.0)),
+        ((4, 5, 6), 0.0),
+        ((4, 5, 6), 0.5),
+        ((4, 5, 6), 1.0),
+        ((4, 5, 6), (0.0, 0.0, 0.0)),
+        ((4, 5, 6), (0.5, 0.0, 0.0)),
+        ((4, 5, 6), (0.5, 0.6, 0.0)),
+        ((4, 5, 6), (0.0, 0.5, 0.6)),
+        ((4, 5, 6), (0.4, 0.5, 0.6)),
+        ((4, 5, 6, 7), 0.0),
+        ((4, 5, 6, 7), 0.5),
+        ((4, 5, 6, 7), 1.0),
+        ((4, 5, 6, 7), (0.0, 0.0, 0.0, 0.0)),
+        ((4, 5, 6, 7), (0.5, 0.0, 0.0, 0.0)),
+        ((4, 5, 6, 7), (0.0, 0.5, 0.0, 0.0)),
+        ((4, 5, 6, 7), (0.5, 0.6, 0.0, 0.0)),
+        ((4, 5, 6, 7), (0.0, 0.5, 0.6, 0.0)),
+        ((4, 5, 6, 7), (0.0, 0.5, 0.6, 0.7)),
+        ((4, 5, 6, 7), (0.4, 0.5, 0.6, 0.7)),
+    ]
+)
 def test_smooth_factor(shape, smooth):
     x = [np.arange(s) for s in shape]
     y = np.arange(0, np.prod(shape)).reshape(shape)
@@ -132,27 +134,25 @@ def test_smooth_factor(shape, smooth):
     assert sp.smooth == pytest.approx(expected_smooth)
 
 
-@pytest.mark.parametrize('shape', [
-    (2,),
-
-    (2, 3),
-    (2, 2),
-
-    (2, 3, 4),
-    (2, 2, 3),
-    (2, 2, 2),
-
-    (2, 3, 4, 5),
-    (2, 2, 3, 4),
-    (2, 2, 2, 3),
-    (2, 2, 2, 2),
-
-    (2, 3, 4, 5, 6),
-    (2, 2, 3, 4, 5),
-    (2, 2, 2, 3, 4),
-    (2, 2, 2, 2, 3),
-    (2, 2, 2, 2, 2),
-])
+@pytest.mark.parametrize(
+    'shape', [
+        (2, ),
+        (2, 3),
+        (2, 2),
+        (2, 3, 4),
+        (2, 2, 3),
+        (2, 2, 2),
+        (2, 3, 4, 5),
+        (2, 2, 3, 4),
+        (2, 2, 2, 3),
+        (2, 2, 2, 2),
+        (2, 3, 4, 5, 6),
+        (2, 2, 3, 4, 5),
+        (2, 2, 2, 3, 4),
+        (2, 2, 2, 2, 3),
+        (2, 2, 2, 2, 2),
+    ]
+)
 def test_nd_2pt_array(shape: tuple):
     xdata = [np.arange(s) for s in shape]
     ydata = np.arange(0, np.prod(shape)).reshape(shape)
@@ -164,18 +164,20 @@ def test_nd_2pt_array(shape: tuple):
     assert ydata_s == pytest.approx(ydata)
 
 
-@pytest.mark.parametrize('shape', [
-    (2,),
-    (3,),
-    (2, 3),
-    (3, 4),
-    (3, 2, 4),
-    (3, 4, 5),
-    (2, 4, 2, 6),
-    (3, 4, 5, 6),
-    (3, 2, 2, 6, 2),
-    (3, 4, 5, 6, 7),
-])
+@pytest.mark.parametrize(
+    'shape', [
+        (2, ),
+        (3, ),
+        (2, 3),
+        (3, 4),
+        (3, 2, 4),
+        (3, 4, 5),
+        (2, 4, 2, 6),
+        (3, 4, 5, 6),
+        (3, 2, 2, 6, 2),
+        (3, 4, 5, 6, 7),
+    ]
+)
 def test_nd_array(shape: tuple):
     xdata = [np.arange(s) for s in shape]
     ydata = np.arange(0, np.prod(shape)).reshape(shape)
diff --git a/tests/test_reshape.py b/tests/test_reshape.py
index b62e46e..5a1f71e 100644
--- a/tests/test_reshape.py
+++ b/tests/test_reshape.py
@@ -1,25 +1,25 @@
-# -*- coding: utf-8 -*-
-
 import pytest
 import numpy as np
 
 from csaps._reshape import (  # noqa
-    umv_coeffs_to_flatten,
-    umv_coeffs_to_canonical,
-    ndg_coeffs_to_flatten,
     ndg_coeffs_to_canonical,
+    ndg_coeffs_to_flatten,
+    umv_coeffs_to_canonical,
+    umv_coeffs_to_flatten,
 )
 
 
-@pytest.mark.parametrize('shape_canonical, shape_flatten, pieces', [
-    ((2, 1), (1, 2), 1),
-    ((3, 6), (1, 18), 6),
-    ((4, 3), (1, 12), 3),
-    ((4, 30), (1, 120), 30),
-    ((4, 5, 2), (2, 20), 5),
-    ((4, 6, 3), (3, 24), 6),
-    ((4, 120, 53), (53, 480), 120),
-])
+@pytest.mark.parametrize(
+    'shape_canonical, shape_flatten, pieces', [
+        ((2, 1), (1, 2), 1),
+        ((3, 6), (1, 18), 6),
+        ((4, 3), (1, 12), 3),
+        ((4, 30), (1, 120), 30),
+        ((4, 5, 2), (2, 20), 5),
+        ((4, 6, 3), (3, 24), 6),
+        ((4, 120, 53), (53, 480), 120),
+    ]
+)
 def test_umv_coeffs_reshape(shape_canonical: tuple, shape_flatten: tuple, pieces: int):
     np.random.seed(1234)
     arr_canonical_expected = np.random.randint(0, 99, size=shape_canonical)
@@ -31,35 +31,38 @@ def test_umv_coeffs_reshape(shape_canonical: tuple, shape_flatten: tuple, pieces
     np.testing.assert_array_equal(arr_canonical_actual, arr_canonical_expected)
 
 
-@pytest.mark.parametrize('shape_canonical, shape_flatten, pieces', [
-    # 1-d 2-ordered
-    ((2, 3), (2, 3), (3,)),
-    ((2, 4), (2, 4), (4,)),
-    ((2, 5), (2, 5), (5,)),
+@pytest.mark.parametrize(
+    'shape_canonical, shape_flatten, pieces',
+    [
+        # 1-d 2-ordered
+        ((2, 3), (2, 3), (3, )),
+        ((2, 4), (2, 4), (4, )),
+        ((2, 5), (2, 5), (5, )),
 
-    # 1-d 3-ordered
-    ((3, 3), (3, 3), (3,)),
-    ((3, 4), (3, 4), (4,)),
-    ((3, 5), (3, 5), (5,)),
+        # 1-d 3-ordered
+        ((3, 3), (3, 3), (3, )),
+        ((3, 4), (3, 4), (4, )),
+        ((3, 5), (3, 5), (5, )),
 
-    # 1-d 4-ordered
-    ((4, 3), (4, 3), (3,)),
-    ((4, 4), (4, 4), (4,)),
-    ((4, 5), (4, 5), (5,)),
+        # 1-d 4-ordered
+        ((4, 3), (4, 3), (3, )),
+        ((4, 4), (4, 4), (4, )),
+        ((4, 5), (4, 5), (5, )),
 
-    # 2-d {2,4}-ordered
-    ((2, 4, 3, 4), (6, 16), (3, 4)),
-    ((4, 2, 3, 3), (12, 6), (3, 3)),
-    ((4, 2, 4, 3), (16, 6), (4, 3)),
-    ((2, 4, 4, 4), (8, 16), (4, 4)),
+        # 2-d {2,4}-ordered
+        ((2, 4, 3, 4), (6, 16), (3, 4)),
+        ((4, 2, 3, 3), (12, 6), (3, 3)),
+        ((4, 2, 4, 3), (16, 6), (4, 3)),
+        ((2, 4, 4, 4), (8, 16), (4, 4)),
 
-    # 2-d {4,4}-ordered
-    ((4, 4, 3, 3), (12, 12), (3, 3)),
+        # 2-d {4,4}-ordered
+        ((4, 4, 3, 3), (12, 12), (3, 3)),
 
-    # 3-d {4,4,4}-ordered
-    ((4, 4, 4, 3, 3, 3), (12, 12, 12), (3, 3, 3)),
-    ((4, 4, 4, 3, 5, 7), (12, 20, 28), (3, 5, 7)),
-])
+        # 3-d {4,4,4}-ordered
+        ((4, 4, 4, 3, 3, 3), (12, 12, 12), (3, 3, 3)),
+        ((4, 4, 4, 3, 5, 7), (12, 20, 28), (3, 5, 7)),
+    ]
+)
 def test_ndg_coeffs_reshape(shape_canonical: tuple, shape_flatten: tuple, pieces: tuple):
     np.random.seed(1234)
     arr_canonical_expected = np.random.randint(0, 99, size=shape_canonical)
diff --git a/tests/test_umv.py b/tests/test_umv.py
index 3040a80..07346f6 100644
--- a/tests/test_umv.py
+++ b/tests/test_umv.py
@@ -1,114 +1,119 @@
-# -*- coding: utf-8 -*-
-
-from itertools import chain, product, permutations
+from itertools import chain, permutations, product
 
+import pytest
 import numpy as np
 from scipy.interpolate import CubicSpline
-import pytest
 
 import csaps
 
 
-@pytest.mark.parametrize('x,y,w', [
-    ([1], [2], None),
-    ([1, 2, 3], [1, 2], None),
-    ([(1, 2, 3), (1, 2, 3)], [1, 2, 3], None),
-    ([1, 2, 3], [1, 2, 3], [1, 1]),
-    ([1, 2, 3], [1, 2, 3], [1, 1, 1, 1]),
-    ([1, 2, 3], [1, 2, 3], [(1, 1, 1), (1, 1, 1)]),
-    ([1, 2, 3], [(1, 2, 3, 4), (1, 2, 3, 4)], None),
-    ([1, 2, 3], np.ones((2, 4, 5)), None),
-    ([1, 2, 3], np.ones((2, 4, 3)), np.ones((2, 4, 4))),
-    ([1, 2, 3], [(1, 2, 3), (1, 2, 3)], [(1, 1, 1, 1), (1, 1, 1, 1)]),
-    ([1, 2, 3], [(1, 2, 3), (1, 2, 3)], [(1, 1, 1), (1, 1, 1), (1, 1, 1)])
-])
+@pytest.mark.parametrize(
+    'x,y,w', [
+        ([1], [2], None), ([1, 2, 3], [1, 2], None), ([(1, 2, 3), (1, 2, 3)], [1, 2, 3], None),
+        ([1, 2, 3], [1, 2, 3], [1, 1]), ([1, 2, 3], [1, 2, 3], [1, 1, 1, 1]),
+        ([1, 2, 3], [1, 2, 3], [(1, 1, 1), (1, 1, 1)]), ([1, 2, 3], [(1, 2, 3, 4), (1, 2, 3, 4)], None),
+        ([1, 2, 3], np.ones((2, 4, 5)), None), ([1, 2, 3], np.ones((2, 4, 3)), np.ones((2, 4, 4))),
+        ([1, 2, 3], [(1, 2, 3), (1, 2, 3)], [(1, 1, 1, 1), (1, 1, 1, 1)]),
+        ([1, 2, 3], [(1, 2, 3), (1, 2, 3)], [(1, 1, 1), (1, 1, 1), (1, 1, 1)])
+    ]
+)
 def test_invalid_data(x, y, w):
     with pytest.raises(ValueError):
         csaps.CubicSmoothingSpline(x, y, weights=w)
 
 
-@pytest.mark.parametrize('y', [
-    # 1D (2, )
-    [2, 4],
-
-    # 2D (2, 2)
-    [(2, 4), (3, 5)],
-
-    # 2D (3, 2)
-    [(2, 4), (3, 5), (4, 6)],
-
-    # 3D (2, 2, 2)
-    [[(1, 2), (3, 4)],
-     [(5, 6), (7, 8)]],
-
-    # 1D (3, )
-    [2, 4, 6],
-
-    # 2D (2, 5)
-    [(1, 2, 3, 4, 5), (3, 4, 5, 6, 7)],
-
-    # 2D (3, 3)
-    [(2, 4, 6), (3, 5, 7), (4, 6, 8)],
-
-    # 3D (2, 2, 3)
-    [[(2, 4, 6), (3, 5, 7)],
-     [(2, 4, 6), (3, 5, 7)]],
-
-    # 1D (4, )
-    [2, 4, 6, 8],
-
-    # 2D (2, 4)
-    [(2, 4, 6, 8), (3, 5, 7, 9)],
-
-    # 2D (3, 4)
-    [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-
-    # 3D (2, 2, 4)
-    [[(2, 4, 6, 8), (3, 5, 7, 9)],
-     [(2, 4, 6, 8), (3, 5, 7, 9)]],
-
-    # 3D (2, 3, 4)
-    [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-     [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
-
-    # 3D (3, 2, 4)
-    [[(2, 4, 6, 8), (3, 5, 7, 9)],
-     [(2, 4, 6, 8), (3, 5, 7, 9)],
-     [(2, 4, 6, 8), (3, 5, 7, 9)]],
-
-    # 3D (3, 3, 4)
-    [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-     [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-     [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
-
-    # 4D (2, 2, 2, 4)
-    [[[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
-     [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]]],
-
-    # 4D (3, 2, 2, 4)
-    [[[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
-     [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
-     [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]]],
-
-    # 4D (3, 2, 3, 4)
-    [[[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
-     [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
-     [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]]],
-
-    # 4D (3, 3, 3, 4)
-    [[[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-      [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-      [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
-
-     [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-      [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-      [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
-
-     [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-      [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
-      [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]]],
-
-])
+@pytest.mark.parametrize(
+    'y',
+    [
+        # 1D (2, )
+        [2, 4],
+
+        # 2D (2, 2)
+        [(2, 4), (3, 5)],
+
+        # 2D (3, 2)
+        [(2, 4), (3, 5), (4, 6)],
+
+        # 3D (2, 2, 2)
+        [[(1, 2), (3, 4)], [(5, 6), (7, 8)]],
+
+        # 1D (3, )
+        [2, 4, 6],
+
+        # 2D (2, 5)
+        [(1, 2, 3, 4, 5), (3, 4, 5, 6, 7)],
+
+        # 2D (3, 3)
+        [(2, 4, 6), (3, 5, 7), (4, 6, 8)],
+
+        # 3D (2, 2, 3)
+        [[(2, 4, 6), (3, 5, 7)], [(2, 4, 6), (3, 5, 7)]],
+
+        # 1D (4, )
+        [2, 4, 6, 8],
+
+        # 2D (2, 4)
+        [(2, 4, 6, 8), (3, 5, 7, 9)],
+
+        # 2D (3, 4)
+        [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+
+        # 3D (2, 2, 4)
+        [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+
+        # 3D (2, 3, 4)
+        [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
+
+        # 3D (3, 2, 4)
+        [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+
+        # 3D (3, 3, 4)
+        [
+            [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+            [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+            [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+        ],
+
+        # 4D (2, 2, 2, 4)
+        [
+            [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+            [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+        ],
+
+        # 4D (3, 2, 2, 4)
+        [
+            [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+            [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+            [[(2, 4, 6, 8), (3, 5, 7, 9)], [(2, 4, 6, 8), (3, 5, 7, 9)]],
+        ],
+
+        # 4D (3, 2, 3, 4)
+        [
+            [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
+            [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
+            [[(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)], [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)]],
+        ],
+
+        # 4D (3, 3, 3, 4)
+        [
+            [
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+            ],
+            [
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+            ],
+            [
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+                [(2, 4, 6, 8), (3, 5, 7, 9), (4, 6, 8, 10)],
+            ],
+        ],
+    ]
+)
 def test_vectorize(y):
     x = np.arange(np.array(y).shape[-1])
 
@@ -116,21 +121,28 @@ def test_vectorize(y):
     np.testing.assert_allclose(ys, y)
 
 
-@pytest.mark.parametrize('shape, axis', chain(
-    *map(product, [
-        # shape
-        [(2,), (4,)],
-        permutations((2, 4), 2),
-        permutations((3, 4, 5), 3),
-        permutations((3, 4, 5, 6), 4)
-    ], [
-        # axis
-        range(-1, 1),
-        range(-2, 2),
-        range(-3, 3),
-        range(-4, 4),
-    ])
-))
+@pytest.mark.parametrize(
+    'shape, axis',
+    chain(
+        *map(
+            product,
+            [
+                # shape
+                [(2, ), (4, )],
+                permutations((2, 4), 2),
+                permutations((3, 4, 5), 3),
+                permutations((3, 4, 5, 6), 4),
+            ],
+            [
+                # axis
+                range(-1, 1),
+                range(-2, 2),
+                range(-3, 3),
+                range(-4, 4),
+            ]
+        )
+    )
+)
 def test_axis(shape, axis):
     y = np.arange(int(np.prod(shape))).reshape(shape)
     x = np.arange(np.array(y).shape[axis])
@@ -145,7 +157,7 @@ def test_axis(shape, axis):
     ndim = int(np.prod(shape)) // shape[axis]
     order = 2 if len(x) < 3 else 4
     pieces = len(x) - 1
-    coeffs_shape = (order, pieces) + shape[:axis] + shape[axis+1:]
+    coeffs_shape = (order, pieces) + shape[:axis] + shape[axis + 1:]
 
     assert ss.breaks == pytest.approx(x)
     assert ss.coeffs.shape == coeffs_shape
@@ -166,10 +178,7 @@ def test_zero_smooth():
 
     ys = sp(x)
 
-    assert ys == pytest.approx([2.440677966101695,
-                                3.355932203389830,
-                                5.186440677966102,
-                                7.016949152542373])
+    assert ys == pytest.approx([2.440677966101695, 3.355932203389830, 5.186440677966102, 7.016949152542373])
 
 
 def test_auto_smooth(univariate_data):
@@ -182,16 +191,18 @@ def test_auto_smooth(univariate_data):
     assert yi == pytest.approx(yi_expected)
 
 
-@pytest.mark.parametrize('x,y,xi,yid', [
-    ([1., 2.], [3., 4.], [1., 1.5, 2.], [3., 3.5, 4.]),
-    ([1., 2., 3.], [3., 4., 5.], [1., 1.5, 2., 2.5, 3.], [3., 3.5, 4., 4.5, 5.]),
-    ([1., 2., 4., 6.], [2., 4., 5., 7.], np.linspace(1., 6., 10), [
-        2.2579392157892, 3.0231172855707, 3.6937304019483,
-        4.21971044584031, 4.65026761247821, 5.04804510368134,
-        5.47288175793241, 5.94265482897362, 6.44293945952166,
-        6.95847986982311
-    ]),
-])
+@pytest.mark.parametrize(
+    'x,y,xi,yid', [
+        ([1., 2.], [3., 4.], [1., 1.5, 2.], [3., 3.5, 4.]),
+        ([1., 2., 3.], [3., 4., 5.], [1., 1.5, 2., 2.5, 3.], [3., 3.5, 4., 4.5, 5.]),
+        (
+            [1., 2., 4., 6.], [2., 4., 5., 7.], np.linspace(1., 6., 10), [
+                2.2579392157892, 3.0231172855707, 3.6937304019483, 4.21971044584031, 4.65026761247821, 5.04804510368134,
+                5.47288175793241, 5.94265482897362, 6.44293945952166, 6.95847986982311
+            ]
+        ),
+    ]
+)
 def test_npoints(x, y, xi, yid):
     sp = csaps.CubicSmoothingSpline(x, y)
     yi = sp(xi)
@@ -199,14 +210,16 @@ def test_npoints(x, y, xi, yid):
     np.testing.assert_allclose(yi, yid)
 
 
-@pytest.mark.parametrize('w,yid', [
-    ([0.5, 1, 0.7, 1.2], [
-        2.39572102230177, 3.13781163365086, 3.78568993197139,
-        4.28992448591238, 4.7009959256016, 5.08290363789967,
-        5.49673867759808, 5.9600748344541, 6.45698622142886,
-        6.97068522346297
-    ])
-])
+@pytest.mark.parametrize(
+    'w,yid', [
+        (
+            [0.5, 1, 0.7, 1.2], [
+                2.39572102230177, 3.13781163365086, 3.78568993197139, 4.28992448591238, 4.7009959256016,
+                5.08290363789967, 5.49673867759808, 5.9600748344541, 6.45698622142886, 6.97068522346297
+            ]
+        )
+    ]
+)
 def test_weighted(w, yid):
     x = [1., 2., 4., 6.]
     y = [2., 4., 5., 7.]