diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 058303f65..61e753ea9 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -15,7 +15,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v5
       with:
-        python-version: "3.12"
+        python-version: "3.13"
 
     - name: install tests dependencies
       run: python -m pip install tox
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 75b82499a..1fa0d4420 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -16,7 +16,7 @@ jobs:
       - name: setup Python
         uses: actions/setup-python@v5
         with:
-          python-version: "3.12"
+          python-version: "3.13"
 
       - name: install tests dependencies
         run: python -m pip install tox
diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml
index 1108289fd..5ae3af316 100644
--- a/.github/workflows/lint.yml
+++ b/.github/workflows/lint.yml
@@ -14,7 +14,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v5
       with:
-        python-version: "3.12"
+        python-version: "3.13"
 
     - name: install tests dependencies
       run: python -m pip install tox
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index d9ebe3475..fe3609798 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -12,7 +12,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-22.04, macos-14, windows-2022]
-        python-version: ["3.9", "3.12"]
+        python-version: ["3.0", "3.13"]
 
     steps:
     - uses: actions/checkout@v4
diff --git a/.gitignore b/.gitignore
index 226933088..1f74b6f04 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
+*.coverage*
 *.pyc
 *.ipynb_checkpoints*
 __pycache__
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 71882236f..b4c827c42 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -9,7 +9,7 @@ version: 2
 build:
   os: ubuntu-22.04
   tools:
-    python: "3.12"
+    python: "3.13"
   
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
diff --git a/pyproject.toml b/pyproject.toml
index db7056b63..374d06293 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -19,7 +19,7 @@ authors = [
     {name = "Michele Ceriotti"}
 ]
 readme = "README.rst"
-requires-python = ">=3.9"
+requires-python = ">=3.10"
 license = {text = "BSD-3-Clause"}
 classifiers = [
     "Development Status :: 4 - Beta",
@@ -75,9 +75,6 @@ include = [
 [tool.coverage.xml]
 output = 'tests/coverage.xml'
 
-[tool.pytest.ini_options]
-testpaths = "tests"
-
 [tool.isort]
 skip = "__init__.py"
 profile = "black"
diff --git a/src/skmatter/_selection.py b/src/skmatter/_selection.py
index 50e55cff3..224c020a7 100644
--- a/src/skmatter/_selection.py
+++ b/src/skmatter/_selection.py
@@ -561,8 +561,7 @@ def score(self, X, y=None):
 
         Parameters
         ----------
-        X : numpy.ndarray of shape [n_samples, n_features]
-            The input samples.
+        X : ignored
         y : ignored
 
         Returns
@@ -570,8 +569,7 @@ def score(self, X, y=None):
         score : numpy.ndarray of (n_to_select_from_)
             :math:`\pi` importance for the given samples or features
         """
-        X, y = validate_data(self, X, y, reset=False)
-
+        validate_data(self, X, y, reset=False)  # present for API consistency
         return self.pi_
 
     def _init_greedy_search(self, X, y, n_to_select):
@@ -746,8 +744,7 @@ def score(self, X, y=None):
         score : numpy.ndarray of (n_to_select_from_)
             :math:`\pi` importance for the given samples or features
         """
-        X, y = validate_data(self, X, y, reset=False)
-
+        validate_data(self, X, y, reset=False)  # present for API consistency
         return self.pi_
 
     def _init_greedy_search(self, X, y, n_to_select):
@@ -941,8 +938,7 @@ def score(self, X, y=None):
         -------
         hausdorff : Hausdorff distances
         """
-        X, y = validate_data(self, X, y, reset=False)
-
+        validate_data(self, X, y, reset=False)
         return self.hausdorff_
 
     def get_distance(self):
@@ -1079,15 +1075,16 @@ def __init__(
         )
 
     def fit(self, X, y=None, warm_start=False):
-
         if self.mixing == 1.0:
             raise ValueError(
-                "Mixing = 1.0 corresponds to traditional FPS."
-                "Please use the FPS class."
+                "Mixing = 1.0 corresponds to traditional FPS. Please use the FPS class."
             )
 
         return super().fit(X, y)
 
+    # docstring is inherited and set from the base class
+    fit.__doc__ = GreedySelector.fit.__doc__
+
     def score(self, X, y=None):
         """Returns the Hausdorff distances of all samples to previous selections.
 
@@ -1104,8 +1101,7 @@ def score(self, X, y=None):
         -------
         hausdorff : Hausdorff distances
         """
-        X, y = validate_data(self, X, y, reset=False)
-
+        validate_data(self, X, y, reset=False)
         return self.hausdorff_
 
     def get_distance(self):
diff --git a/src/skmatter/decomposition/_kernel_pcovr.py b/src/skmatter/decomposition/_kernel_pcovr.py
index 65fe39a3e..825a0cf92 100644
--- a/src/skmatter/decomposition/_kernel_pcovr.py
+++ b/src/skmatter/decomposition/_kernel_pcovr.py
@@ -474,7 +474,7 @@ def score(self, X, y):
         """
         check_is_fitted(self, ["pkt_", "X_fit_"])
 
-        X, y = validate_data(self, X, y, reset=False)
+        X = validate_data(self, X, reset=False)
 
         K_NN = self._get_kernel(self.X_fit_, self.X_fit_)
         K_VN = self._get_kernel(X, self.X_fit_)
diff --git a/src/skmatter/decomposition/_pcovr.py b/src/skmatter/decomposition/_pcovr.py
index 8cdd24680..bc094a720 100644
--- a/src/skmatter/decomposition/_pcovr.py
+++ b/src/skmatter/decomposition/_pcovr.py
@@ -10,7 +10,7 @@
 from sklearn.decomposition._pca import _infer_dimension
 from sklearn.linear_model import LinearRegression, Ridge, RidgeCV
 from sklearn.linear_model._base import LinearModel
-from sklearn.utils import check_random_state
+from sklearn.utils import check_array, check_random_state
 from sklearn.utils._arpack import _init_arpack_v0
 from sklearn.utils.extmath import randomized_svd, stable_cumsum, svd_flip
 from sklearn.utils.validation import check_is_fitted, validate_data
@@ -585,7 +585,7 @@ def predict(self, X=None, T=None):
             X = validate_data(self, X, reset=False)
             return X @ self.pxy_
         else:
-            T = validate_data(self, T, reset=False)
+            T = check_array(T)
             return T @ self.pty_
 
     def transform(self, X=None):
diff --git a/src/skmatter/utils/_pcovr_utils.py b/src/skmatter/utils/_pcovr_utils.py
index 15286e341..8852a6386 100644
--- a/src/skmatter/utils/_pcovr_utils.py
+++ b/src/skmatter/utils/_pcovr_utils.py
@@ -9,18 +9,30 @@
 
 
 def check_lr_fit(regressor, X, y):
-    r"""
+    """
     Checks that a (linear) regressor is fitted, and if not,
-    fits it with the provided data
-
-    :param regressor: sklearn-style regressor
-    :type regressor: object
-    :param X: feature matrix with which to fit the regressor
-        if it is not already fitted
-    :type X: array
-    :param y: target values with which to fit the regressor
-        if it is not already fitted
-    :type y: array
+    fits it with the provided data.
+
+    Parameters
+    ----------
+    regressor : object
+        sklearn-style regressor
+    X : array-like
+        Feature matrix with which to fit the regressor if it is not already fitted
+    y : array-like
+        Target values with which to fit the regressor if it is not already fitted
+
+    Returns
+    -------
+    fitted_regressor : object
+        The fitted regressor. If input regressor was already fitted and compatible with
+        the data, returns a deep copy. Otherwise returns a newly fitted regressor.
+
+    Raises
+    ------
+    ValueError
+        If the fitted regressor's coefficients dimensions are incompatible with the
+        target space.
     """
     try:
         check_is_fitted(regressor)
@@ -32,18 +44,18 @@ def check_lr_fit(regressor, X, y):
         # Check compatibility with y
         if fitted_regressor.coef_.ndim != y.ndim:
             raise ValueError(
-                "The regressor coefficients have a dimension incompatible "
-                "with the supplied target space. "
-                "The coefficients have dimension %d and the targets "
-                "have dimension %d" % (fitted_regressor.coef_.ndim, y.ndim)
+                "The regressor coefficients have a dimension incompatible with the "
+                "supplied target space. The coefficients have dimension "
+                f"{fitted_regressor.coef_.ndim} and the targets have dimension "
+                f"{y.ndim}"
             )
         elif y.ndim == 2:
             if fitted_regressor.coef_.shape[0] != y.shape[1]:
                 raise ValueError(
-                    "The regressor coefficients have a shape incompatible "
-                    "with the supplied target space. "
-                    "The coefficients have shape %r and the targets "
-                    "have shape %r" % (fitted_regressor.coef_.shape, y.shape)
+                    "The regressor coefficients have a shape incompatible with the "
+                    "supplied target space. The coefficients have shape "
+                    f"{fitted_regressor.coef_.shape} and the targets have shape "
+                    f"{y.shape}"
                 )
 
     except NotFittedError:
@@ -54,20 +66,37 @@ def check_lr_fit(regressor, X, y):
 
 
 def check_krr_fit(regressor, K, X, y):
-    r"""
+    """
     Checks that a (kernel ridge) regressor is fitted, and if not,
-    fits it with the provided data
-
-    :param regressor: sklearn-style regressor
-    :type regressor: object
-    :param K: kernel matrix with which to fit the regressor
-        if it is not already fitted
-    :type K: array
-    :param X: feature matrix with which to check the regressor
-    :type X: array
-    :param y: target values with which to fit the regressor
-        if it is not already fitted
-    :type y: array
+    fits it with the provided data.
+
+    Parameters
+    ----------
+    regressor : object
+        sklearn-style regressor
+    K : array-like
+        Kernel matrix with which to fit the regressor if it is not already fitted
+    X : array-like
+        Feature matrix with which to check the regressor
+    y : array-like
+        Target values with which to fit the regressor if it is not already fitted
+
+    Returns
+    -------
+    fitted_regressor : object
+        The fitted regressor. If input regressor was already fitted and compatible with
+        the data, returns a deep copy. Otherwise returns a newly fitted regressor.
+
+    Raises
+    ------
+    ValueError
+        If the fitted regressor's coefficients dimensions are incompatible with the
+        target space.
+
+    Notes
+    -----
+    For unfitted regressors, sets the kernel to "precomputed" before fitting with the
+    provided kernel matrix K to avoid recomputation.
     """
     try:
         check_is_fitted(regressor)
@@ -79,18 +108,18 @@ def check_krr_fit(regressor, K, X, y):
         # Check compatibility with y
         if fitted_regressor.dual_coef_.ndim != y.ndim:
             raise ValueError(
-                "The regressor coefficients have a dimension incompatible "
-                "with the supplied target space. "
-                "The coefficients have dimension %d and the targets "
-                "have dimension %d" % (fitted_regressor.dual_coef_.ndim, y.ndim)
+                "The regressor coefficients have a dimension incompatible with the "
+                "supplied target space. The coefficients have dimension "
+                f"{fitted_regressor.dual_coef_.ndim} and the targets have dimension "
+                f"{y.ndim}"
             )
         elif y.ndim == 2:
             if fitted_regressor.dual_coef_.shape[1] != y.shape[1]:
                 raise ValueError(
-                    "The regressor coefficients have a shape incompatible "
-                    "with the supplied target space. "
-                    "The coefficients have shape %r and the targets "
-                    "have shape %r" % (fitted_regressor.dual_coef_.shape, y.shape)
+                    "The regressor coefficients have a shape incompatible with the "
+                    "supplied target space. The coefficients have shape "
+                    f"{fitted_regressor.dual_coef_.shape} and the targets have shape "
+                    f"{y.shape}"
                 )
 
     except NotFittedError:
diff --git a/tests/.coverage.tsf-492-wpa-0-247.epfl.ch.11311.XpjwIfdx b/tests/.coverage.tsf-492-wpa-0-247.epfl.ch.11311.XpjwIfdx
deleted file mode 100644
index 5d637632b..000000000
Binary files a/tests/.coverage.tsf-492-wpa-0-247.epfl.ch.11311.XpjwIfdx and /dev/null differ
diff --git a/tests/test_feature_pcov_fps.py b/tests/test_feature_pcov_fps.py
index 321cc78ee..e6910f9a1 100644
--- a/tests/test_feature_pcov_fps.py
+++ b/tests/test_feature_pcov_fps.py
@@ -24,11 +24,12 @@ def test_restart(self):
 
     def test_no_mixing_1(self):
         """Check that the model throws an error when mixing = 1.0."""
+        selector = PCovFPS(n_to_select=1, mixing=1.0)
         with self.assertRaises(ValueError) as cm:
-            _ = PCovFPS(n_to_select=1, mixing=1.0)
+            selector.fit(self.X, y=self.y)
         self.assertEqual(
             str(cm.exception),
-            "Mixing = 1.0 corresponds to traditional FPS." "Please use the FPS class.",
+            "Mixing = 1.0 corresponds to traditional FPS. Please use the FPS class.",
         )
 
 
diff --git a/tests/test_greedy_selector.py b/tests/test_greedy_selector.py
index 0bfe6de99..f85dbeb34 100644
--- a/tests/test_greedy_selector.py
+++ b/tests/test_greedy_selector.py
@@ -73,7 +73,7 @@ def test_bad_transform(self):
             _ = selector.transform(self.X[:, :3])
         self.assertEqual(
             str(cm.exception),
-            "X has a different shape than during fitting. Reshape your data.",
+            "X has 3 features, but GreedyTester is expecting 10 features as input.",
         )
 
     def test_no_nfeatures(self):
@@ -124,8 +124,8 @@ def test_size_input(self):
             selector_feature.fit(X)
         self.assertEqual(
             str(cm.exception),
-            f"Found array with 1 feature(s) (shape={X.shape})"
-            " while a minimum of 2 is required.",
+            f"Found array with 1 feature(s) (shape={X.shape}) while a minimum of 2 is "
+            "required by GreedyTester.",
         )
 
         X = X.reshape(1, -1)
@@ -135,7 +135,7 @@ def test_size_input(self):
         self.assertEqual(
             str(cm.exception),
             f"Found array with 1 sample(s) (shape={X.shape}) while a minimum of 2 is "
-            "required.",
+            "required by GreedyTester.",
         )
 
 
diff --git a/tests/test_kernel_pcovr.py b/tests/test_kernel_pcovr.py
index e4bbda52e..80adf584e 100644
--- a/tests/test_kernel_pcovr.py
+++ b/tests/test_kernel_pcovr.py
@@ -182,6 +182,8 @@ def test_centerer(self):
         self.assertTrue(hasattr(kpcovr, "centerer_"))
         _ = kpcovr.predict(self.X)
         _ = kpcovr.transform(self.X)
+
+        print(self.Y.shape)
         _ = kpcovr.score(self.X, self.Y)
 
     def test_prefit_regressor(self):
@@ -255,7 +257,7 @@ def test_incompatible_coef_shape(self):
 
         # Dimension mismatch
         with self.assertRaises(ValueError) as cm:
-            kpcovr.fit(self.X, self.Y[:, 0])
+            kpcovr.fit(self.X, np.zeros(self.Y.shape + (2,)))
         self.assertTrue(
             str(cm.exception),
             "The regressor coefficients have a dimension incompatible "
diff --git a/tests/test_pcovr.py b/tests/test_pcovr.py
index e589978d2..2059eed44 100644
--- a/tests/test_pcovr.py
+++ b/tests/test_pcovr.py
@@ -491,13 +491,12 @@ def test_incompatible_coef_shape(self):
 
         # Dimension mismatch
         with self.assertRaises(ValueError) as cm:
-            pcovr.fit(self.X, self.Y.squeeze())
+            pcovr.fit(self.X, np.zeros((self.Y.shape[0], 2)))
         self.assertEqual(
             str(cm.exception),
-            "The regressor coefficients have a dimension incompatible "
-            "with the supplied target space. "
-            "The coefficients have dimension %d and the targets "
-            "have dimension %d" % (regressor.coef_.ndim, self.Y.squeeze().ndim),
+            "The regressor coefficients have a dimension incompatible with the "
+            "supplied target space. The coefficients have dimension 1 and the targets "
+            "have dimension 2",
         )
 
         # Shape mismatch (number of targets)
diff --git a/tests/test_sample_pcov_fps.py b/tests/test_sample_pcov_fps.py
index 7679abb0a..b6ed08662 100644
--- a/tests/test_sample_pcov_fps.py
+++ b/tests/test_sample_pcov_fps.py
@@ -24,11 +24,12 @@ def test_restart(self):
 
     def test_no_mixing_1(self):
         """Check that the model throws an error when mixing = 1.0."""
+        selector = PCovFPS(n_to_select=1, mixing=1.0)
         with self.assertRaises(ValueError) as cm:
-            _ = PCovFPS(n_to_select=1, mixing=1.0)
+            selector.fit(self.X, y=self.y)
         self.assertEqual(
             str(cm.exception),
-            "Mixing = 1.0 corresponds to traditional FPS." "Please use the FPS class.",
+            "Mixing = 1.0 corresponds to traditional FPS. Please use the FPS class.",
         )