From 7d72b6d11d00b83ff3597911fd96f39642cda17d Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Sat, 12 Sep 2026 12:48:45 -0600 Subject: [PATCH 01/32] WIP: array api standard for torch cpu - staging ground to have array api use - testing for torch cpu only right now - testing regressor only right now - test to match numpy and torch r2 score - test that predictions are in same namespace - many issues with sklearn conformance in classifier - more details to come about onehotencoder issues --- pyproject.toml | 2 + src/gfdl/model.py | 17 ++++-- src/gfdl/tests/test_regression.py | 99 +++++++++++++++++++++++++++++++ 3 files changed, 114 insertions(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 49be48a..c789b3a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,8 @@ classifiers = [ dependencies = ["numpy>=2.0.0", "scikit-learn>=1.5.0", "scipy>=1.13.0", + "array-api-compat", + "torch>=2.12.0,<=2.14.0", "packaging>=24.0"] [project.urls] diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 2f4f216..3a9c67e 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -2,6 +2,7 @@ Estimators for gradient free deep learning. """ +import array_api_compat import numpy as np import scipy from scipy.special import logsumexp @@ -44,12 +45,14 @@ def __init__( self.rtol = rtol def fit(self, X, Y): + xp = array_api_compat.get_namespace(X, Y) + # Assumption : X, Y have been pre-processed. # X shape: (n_samples, n_features) # Y shape: (n_samples, n_classes-1) if self.reg_alpha is not None and self.reg_alpha < 0.0: raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).") - hidden_layer_sizes = np.asarray(self.hidden_layer_sizes) + hidden_layer_sizes = xp.asarray(self.hidden_layer_sizes) if hidden_layer_sizes.min() < 1: raise ValueError("hidden_layer_sizes must be > 0, " f"got {hidden_layer_sizes}") @@ -87,6 +90,8 @@ def fit(self, X, Y): Hs = [] H_prev = X for w, b in zip(self.W_, self.b_, strict=False): + w = xp.asarray(w, device=H_prev.device,) + b = xp.asarray(b, device=H_prev.device,) Z = H_prev @ w.T + b # (n_samples, n_hidden) H_prev = self._activation_fn(Z) Hs.append(H_prev) @@ -95,7 +100,7 @@ def fit(self, X, Y): # or (n_samples, sum_hidden) if self.direct_links: Hs.append(X) - D = np.hstack(Hs) + D = xp.concat(Hs, axis=1) # beta shape: (sum_hidden+n_features, n_classes-1) # or (sum_hidden, n_classes-1) @@ -103,7 +108,7 @@ def fit(self, X, Y): # If reg_alpha is None, use direct solve using # MoorePenrose Pseudo-Inverse, otherwise use ridge regularized form. if self.reg_alpha is None: - self.coeff_ = np.linalg.pinv(D, rtol=self.rtol) @ Y + self.coeff_ = xp.linalg.pinv(D, rtol=self.rtol) @ Y else: ridge = Ridge(alpha=self.reg_alpha, fit_intercept=False) ridge.fit(D, Y) @@ -229,16 +234,20 @@ def partial_fit(self, X, y): def predict(self, X): check_is_fitted(self) + xp = array_api_compat.get_namespace(X) + Hs = [] H_prev = X for W, b in zip(self.W_, self.b_, strict=False): + W = xp.asarray(W, device=H_prev.device) + b = xp.asarray(b, device=H_prev.device) Z = H_prev @ W.T + b # (n, m) H_prev = self._activation_fn(Z) Hs.append(H_prev) if self.direct_links: Hs.append(X) - D = np.hstack(Hs) + D = xp.concat(Hs, axis=1) out = D @ self.coeff_ return out diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index c0cab55..5a759fc 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -1,5 +1,11 @@ +import os + +import array_api_compat import numpy as np import pytest +import torch +from numpy.testing import assert_allclose +from sklearn import config_context from sklearn.base import clone from sklearn.datasets import fetch_openml, make_regression from sklearn.metrics import r2_score @@ -257,3 +263,96 @@ def test_preserve_class_inputs(): for k, v in actual.items(): assert v == expected[k] assert isinstance(v, type(expected[k])) + + +@pytest.mark.parametrize( + "estimator, reg_alpha, n_features", + [ + (GFDLRegressor, 1e-1, 40), + (GFDLRegressor, 1e-1, 400), + (GFDLRegressor, 2, 40), + (GFDLRegressor, 2, 400), + (GFDLRegressor, None, 40), + (GFDLRegressor, None, 400), + ] +) +def test_torch_cpu_matches_numpy(estimator, + reg_alpha, + n_features, + ): + # numpy and torch accuracies should be close up to atol + os.environ["SCIPY_ARRAY_API"] = "1" + report = r2_score + rng = np.random.default_rng(seed=42) + acc_np_s = [] + acc_torch_s = [] + for _ in range(10): + random_state = rng.integers(low=0, high=100, size=1)[0] + X_np, y_np = make_regression( + n_samples=10_000, + n_features=n_features, + n_informative=int(n_features / 10), + n_targets=1, + random_state=random_state, + ) + X_torch = torch.asarray(X_np, device="cpu",) + y_torch = torch.asarray(y_np, device="cpu",) + with config_context(array_api_dispatch=True): + model = estimator(reg_alpha=reg_alpha) + model.fit(X_np, y_np) + y_pred = model.predict(X_torch) + acc_torch = report( + y_torch, + y_pred, + ) + acc_torch_s.append(acc_torch) + # test against numpy which is current usage + with config_context(array_api_dispatch=True): + model = estimator(reg_alpha=reg_alpha) + model.fit(X_np, y_np) + y_pred = model.predict(X_np) + acc_np = report( + y_np, + y_pred, + ) + acc_np_s.append(acc_np) + assert_allclose(acc_torch_s, acc_np_s, atol=1e-3) + + +@pytest.mark.parametrize( + "estimator, namespace, device, X_dtype, y_dtype", + [ + (GFDLRegressor, np, "cpu", np.float64, np.float64), + (GFDLRegressor, torch, "cpu", torch.float64, torch.float64), + ] +) +def test_predict_same_namespace(estimator, + namespace, + device, + X_dtype, + y_dtype, + ): + # ValueError thrown .predict has + # different namespace from X and y + os.environ["SCIPY_ARRAY_API"] = "1" + + X, y = make_regression( + n_samples=1_000, + n_features=10, + n_informative=4, + n_targets=1, + random_state=42, + ) + X = namespace.asarray(X, dtype=X_dtype, device=device) + y = namespace.asarray(y, dtype=y_dtype, device=device) + xp = array_api_compat.get_namespace(X, y) + with config_context(array_api_dispatch=True): + model = estimator() + model.fit(X, y) + y_pred = model.predict(X) + + xp_y_pred = array_api_compat.get_namespace(y_pred) + if xp_y_pred != xp: + raise ValueError( + ".predict output and X, y are not in the same namespace" + ) From aca08e973ca298f4e4f838f8a67fee5e7e1384c2 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Sun, 13 Sep 2026 22:32:02 -0600 Subject: [PATCH 02/32] WIP: ensure float, check and convert arrays - Two private utility functions for the following: - Check and/or convert all to X's dtype, namespace, & device - Ensure that X is of a float type - Test function to ensure predictions have same context as X - Test function that fitted attributes have same context as X - Test function that ints are handled gracefulyl --- src/gfdl/model.py | 75 +++++++++++-- src/gfdl/tests/test_regression.py | 172 ++++++++++++++++++++++++------ 2 files changed, 204 insertions(+), 43 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 3a9c67e..892d43f 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -45,7 +45,9 @@ def __init__( self.rtol = rtol def fit(self, X, Y): - xp = array_api_compat.get_namespace(X, Y) + xp = array_api_compat.get_namespace(X) + X = _ensure_float_X(xp, X) + Y = _check_and_convert_array(X, Y) # Assumption : X, Y have been pre-processed. # X shape: (n_samples, n_features) @@ -68,30 +70,40 @@ def fit(self, X, Y): rng = self.get_generator(self.seed) self.W_.append( + _check_and_convert_array( + X, self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) - ) + ), + ) ) self.b_.append( + _check_and_convert_array( + X, self._weight_mode(1, hidden_layer_sizes[0], rng=rng) - .reshape(-1) + .reshape(-1), + ) ) for i, layer in enumerate(hidden_layer_sizes[1:]): # (n_hidden, n_features) self.W_.append( - self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,) + _check_and_convert_array( + X, + self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), + ) ) # (n_hidden,) self.b_.append( - self._weight_mode(1, layer, rng=rng,).reshape(-1) + _check_and_convert_array( + X, + self._weight_mode(1, layer, rng=rng,).reshape(-1), + ) ) # hypothesis space shape: (n_layers,) Hs = [] H_prev = X for w, b in zip(self.W_, self.b_, strict=False): - w = xp.asarray(w, device=H_prev.device,) - b = xp.asarray(b, device=H_prev.device,) Z = H_prev @ w.T + b # (n_samples, n_hidden) H_prev = self._activation_fn(Z) Hs.append(H_prev) @@ -239,8 +251,8 @@ def predict(self, X): Hs = [] H_prev = X for W, b in zip(self.W_, self.b_, strict=False): - W = xp.asarray(W, device=H_prev.device) - b = xp.asarray(b, device=H_prev.device) + W = _check_and_convert_array(X, W) + b = _check_and_convert_array(X, b) Z = H_prev @ W.T + b # (n, m) H_prev = self._activation_fn(Z) Hs.append(H_prev) @@ -248,7 +260,10 @@ def predict(self, X): if self.direct_links: Hs.append(X) D = xp.concat(Hs, axis=1) - out = D @ self.coeff_ + out = D @ _check_and_convert_array( + X, + self.coeff_ + ) return out @@ -1275,3 +1290,43 @@ def predict(self, X): check_is_fitted(self) X = validate_data(self, X, reset=False) return super().predict(X) + + +def _check_and_convert_array(X1, X2): + """Convert second array to namespace, device, dtype of first if not already""" + xp = array_api_compat.get_namespace(X1) + kwargs = {} + + if hasattr(X1, "dtype") and getattr(X2, "dtype", None) != X1.dtype: + kwargs["dtype"] = X1.dtype + + x1_device = getattr(X1, "device", None) + x2_device = getattr(X2, "device", None) + + if x1_device is not None and x1_device != x2_device: + kwargs["device"] = x1_device + + if kwargs: + return xp.asarray(X2, **kwargs) + + return X2 + + +def _ensure_float_X(xp, X): + """Make design matrix floating point numbers""" + if xp.isdtype(X.dtype, "real floating"): + return X + + device = getattr(X, "device", None) + device_str = str(device).lower() + + if "cuda" in device_str or "gpu" in device_str or "tpu" in device_str: + dtype = xp.float32 + else: + dtype = xp.float64 + + kwargs = {"dtype": dtype} + if device is not None: + kwargs["device"] = device + + return xp.asarray(X, **kwargs) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 5a759fc..a7172e0 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -266,22 +266,29 @@ def test_preserve_class_inputs(): @pytest.mark.parametrize( - "estimator, reg_alpha, n_features", + "reg_alpha, n_features, hidden_layer_sizes", [ - (GFDLRegressor, 1e-1, 40), - (GFDLRegressor, 1e-1, 400), - (GFDLRegressor, 2, 40), - (GFDLRegressor, 2, 400), - (GFDLRegressor, None, 40), - (GFDLRegressor, None, 400), + (1e-1, 40, (100,)), + (1e-1, 400, (100,)), + (2, 40, (100,)), + (2, 400, (100,)), + (None, 40, (100,)), + (None, 400, (100,)), + (1e-1, 40, (100, 100,)), + (1e-1, 400, (100, 100,)), + (2, 40, (100, 100,)), + (2, 400, (100, 100,)), + (None, 40, (100, 100,)), + (None, 400, (100, 100,)), ] ) -def test_torch_cpu_matches_numpy(estimator, - reg_alpha, - n_features, - ): - # numpy and torch accuracies should be close up to atol +def test_torch_matches_numpy(reg_alpha, + n_features, + hidden_layer_sizes, + ): + """NumPy and Torch predictions should be close up to atol""" os.environ["SCIPY_ARRAY_API"] = "1" + estimator = GFDLRegressor report = r2_score rng = np.random.default_rng(seed=42) acc_np_s = [] @@ -298,7 +305,10 @@ def test_torch_cpu_matches_numpy(estimator, X_torch = torch.asarray(X_np, device="cpu",) y_torch = torch.asarray(y_np, device="cpu",) with config_context(array_api_dispatch=True): - model = estimator(reg_alpha=reg_alpha) + model = estimator(reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + seed=random_state, + ) model.fit(X_np, y_np) y_pred = model.predict(X_torch) acc_torch = report( @@ -308,7 +318,10 @@ def test_torch_cpu_matches_numpy(estimator, acc_torch_s.append(acc_torch) # test against numpy which is current usage with config_context(array_api_dispatch=True): - model = estimator(reg_alpha=reg_alpha) + model = estimator(reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + seed=random_state, + ) model.fit(X_np, y_np) y_pred = model.predict(X_np) acc_np = report( @@ -320,21 +333,24 @@ def test_torch_cpu_matches_numpy(estimator, @pytest.mark.parametrize( - "estimator, namespace, device, X_dtype, y_dtype", + "namespace, fit_dtype, predict_dtype", [ - (GFDLRegressor, np, "cpu", np.float64, np.float64), - (GFDLRegressor, torch, "cpu", torch.float64, torch.float64), + (np, np.float64, np.float32), + (torch, torch.float64, torch.float32), + (np, np.float32, np.float64), + (torch, torch.float32, torch.float64), ] ) -def test_predict_same_namespace(estimator, - namespace, - device, - X_dtype, - y_dtype, - ): +def test_predictor_context(namespace, + fit_dtype, + predict_dtype, + ): # ValueError thrown .predict has - # different namespace from X and y + # different namespace, dtype + # from its arguments os.environ["SCIPY_ARRAY_API"] = "1" + estimator = GFDLRegressor + device = "cpu" X, y = make_regression( n_samples=1_000, @@ -343,16 +359,106 @@ def test_predict_same_namespace(estimator, n_targets=1, random_state=42, ) - X = namespace.asarray(X, dtype=X_dtype, device=device) - y = namespace.asarray(y, dtype=y_dtype, device=device) - xp = array_api_compat.get_namespace(X, y) + X1 = namespace.asarray(X, dtype=fit_dtype, device=device) + X2 = namespace.asarray(X, dtype=predict_dtype, device=device) + y = namespace.asarray(y, dtype=fit_dtype, device=device) with config_context(array_api_dispatch=True): - model = estimator() - model.fit(X, y) - y_pred = model.predict(X) + y_pred = estimator().fit(X1, y).predict(X2) - xp_y_pred = array_api_compat.get_namespace(y_pred) - if xp_y_pred != xp: + xp_ypred = array_api_compat.get_namespace(y_pred) + xp_X2 = array_api_compat.get_namespace(X2) + if xp_ypred != xp_X2: raise ValueError( - ".predict output and X, y are not in the same namespace" + ".predict output and input are not in the same namespace" + ) + if predict_dtype != y_pred.dtype: + raise ValueError( + ".predict output and input are not the same dtype" + ) + + +@pytest.mark.parametrize( + "namespace, X_dtype, y_dtype", + [ + (np, np.float64, np.int64), + (np, np.int64, np.float64), + (torch, torch.float64, torch.int64), + ] +) +def test_int_array_api(namespace, + X_dtype, + y_dtype, + ): + """Integer arrays shall be handled gracefully""" + os.environ["SCIPY_ARRAY_API"] = "1" + estimator = GFDLRegressor + + X, y = make_regression( + n_samples=1_000, + n_features=10, + n_informative=4, + n_targets=1, + random_state=42, + ) + X = namespace.asarray(X, dtype=X_dtype, device="cpu") + y = namespace.asarray(y, dtype=y_dtype, device="cpu") + with config_context(array_api_dispatch=True): + estimator().fit(X, y) + + +@pytest.mark.parametrize( + "namespace, X_dtype, y_dtype", + [ + (np, np.float64, np.float64), + (np, np.float32, np.float32), + (torch, torch.float64, torch.float64), + (torch, torch.float32, torch.float32), + (torch, torch.float64, torch.float32), + ] +) +def test_fit_attr_context(namespace, + X_dtype, + y_dtype, + ): + """Fitted attributes shall be same as design matrix""" + os.environ["SCIPY_ARRAY_API"] = "1" + estimator = GFDLRegressor + device = "cpu" + + X, y = make_regression( + n_samples=1_000, + n_features=10, + n_informative=4, + n_targets=1, + random_state=42, + ) + X = namespace.asarray(X, dtype=X_dtype, device=device) + y = namespace.asarray(y, dtype=y_dtype, device=device) + with config_context(array_api_dispatch=True): + model = estimator().fit(X, y) + + dtypes = [] + devices = [] + spaces = [] + W_ = model.W_ + b_ = model.b_ + for w, b in zip(W_, b_, strict=False): + dtypes.extend((b.dtype, w.dtype)) + devices.extend((b.device, w.device)) + spaces.extend( + (array_api_compat.get_namespace(w), + array_api_compat.get_namespace(b),) ) + + coeff_ = model.coeff_ + dtypes.append(coeff_.dtype) + devices.append(coeff_.device) + spaces.append( + array_api_compat.get_namespace(coeff_) + ) + + assert all(dtype == X.dtype for dtype in dtypes) + assert all(device == X.device for device in devices) + + X_space = array_api_compat.get_namespace(X) + assert all(space == X_space for space in spaces) From 4b6db3eafc1e90914a99175e4204b2d091e03d14 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Sat, 19 Sep 2026 11:14:16 -0700 Subject: [PATCH 03/32] WIP: move_to from private utility funcs - Replace custom private utility functions with move_to - Is existing utility from sklearn - Achieves same functionality and more - array_api_compat no longer used in imports - Torch download removed from TOML - Previous tests pass pre PR (on local Windows) - Torch matches NumPy test passes - Other new tests related to namespace, dtype, device fail - Those will be replaced in next commit with sklearn equiv - Remaining ruff errors not related to this PR scope --- pyproject.toml | 1 - src/gfdl/model.py | 83 ++++++++++++++--------------------------------- 2 files changed, 25 insertions(+), 59 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index c789b3a..294b784 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,7 +29,6 @@ dependencies = ["numpy>=2.0.0", "scikit-learn>=1.5.0", "scipy>=1.13.0", "array-api-compat", - "torch>=2.12.0,<=2.14.0", "packaging>=24.0"] [project.urls] diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 892d43f..aff1dd8 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -2,7 +2,6 @@ Estimators for gradient free deep learning. """ -import array_api_compat import numpy as np import scipy from scipy.special import logsumexp @@ -16,6 +15,10 @@ from sklearn.linear_model import Ridge from sklearn.preprocessing import OneHotEncoder from sklearn.utils import column_or_1d +from sklearn.utils._array_api import ( + get_namespace_and_device, + move_to, +) from sklearn.utils.metaestimators import available_if from sklearn.utils.multiclass import check_classification_targets, unique_labels from sklearn.utils.validation import check_is_fitted, validate_data @@ -45,9 +48,8 @@ def __init__( self.rtol = rtol def fit(self, X, Y): - xp = array_api_compat.get_namespace(X) - X = _ensure_float_X(xp, X) - Y = _check_and_convert_array(X, Y) + xp, _, device = get_namespace_and_device(X) + Y = move_to(Y, xp=xp, device=device) # Assumption : X, Y have been pre-processed. # X shape: (n_samples, n_features) @@ -70,33 +72,37 @@ def fit(self, X, Y): rng = self.get_generator(self.seed) self.W_.append( - _check_and_convert_array( - X, + move_to( self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) ), + xp=xp, + device=device, ) ) self.b_.append( - _check_and_convert_array( - X, + move_to( self._weight_mode(1, hidden_layer_sizes[0], rng=rng) .reshape(-1), + xp=xp, + device=device, ) ) for i, layer in enumerate(hidden_layer_sizes[1:]): # (n_hidden, n_features) self.W_.append( - _check_and_convert_array( - X, + move_to( self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), + xp=xp, + device=device, ) ) # (n_hidden,) self.b_.append( - _check_and_convert_array( - X, + move_to( self._weight_mode(1, layer, rng=rng,).reshape(-1), + xp=xp, + device=device, ) ) @@ -246,13 +252,13 @@ def partial_fit(self, X, y): def predict(self, X): check_is_fitted(self) - xp = array_api_compat.get_namespace(X) + xp, _, device = get_namespace_and_device(X) Hs = [] H_prev = X for W, b in zip(self.W_, self.b_, strict=False): - W = _check_and_convert_array(X, W) - b = _check_and_convert_array(X, b) + W = move_to(W, xp=xp, device=device) + b = move_to(b, xp=xp, device=device) Z = H_prev @ W.T + b # (n, m) H_prev = self._activation_fn(Z) Hs.append(H_prev) @@ -260,9 +266,10 @@ def predict(self, X): if self.direct_links: Hs.append(X) D = xp.concat(Hs, axis=1) - out = D @ _check_and_convert_array( - X, - self.coeff_ + out = D @ move_to( + self.coeff_, + xp=xp, + device=device, ) return out @@ -1290,43 +1297,3 @@ def predict(self, X): check_is_fitted(self) X = validate_data(self, X, reset=False) return super().predict(X) - - -def _check_and_convert_array(X1, X2): - """Convert second array to namespace, device, dtype of first if not already""" - xp = array_api_compat.get_namespace(X1) - kwargs = {} - - if hasattr(X1, "dtype") and getattr(X2, "dtype", None) != X1.dtype: - kwargs["dtype"] = X1.dtype - - x1_device = getattr(X1, "device", None) - x2_device = getattr(X2, "device", None) - - if x1_device is not None and x1_device != x2_device: - kwargs["device"] = x1_device - - if kwargs: - return xp.asarray(X2, **kwargs) - - return X2 - - -def _ensure_float_X(xp, X): - """Make design matrix floating point numbers""" - if xp.isdtype(X.dtype, "real floating"): - return X - - device = getattr(X, "device", None) - device_str = str(device).lower() - - if "cuda" in device_str or "gpu" in device_str or "tpu" in device_str: - dtype = xp.float32 - else: - dtype = xp.float64 - - kwargs = {"dtype": dtype} - if device is not None: - kwargs["device"] = device - - return xp.asarray(X, **kwargs) From 16535037cc763ee5cb5dbc0162136a1f0840b4ef Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 21 Sep 2026 16:37:03 -0600 Subject: [PATCH 04/32] WIP: activation functions with xp - Replace np with xp in some activations - I think scipy uses array API standard - When SCIPY_ARRAY_API is set - Hence no change to scipy.special funcs - Responds to eiviani-lanl comment - No speed change from np.tanh to xp.tanh --- src/gfdl/activations.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/gfdl/activations.py b/src/gfdl/activations.py index e4aae04..0cf50a3 100644 --- a/src/gfdl/activations.py +++ b/src/gfdl/activations.py @@ -4,6 +4,7 @@ import numpy as np import scipy +import array_api_compat def relu(z): @@ -17,10 +18,11 @@ def relu(z): Returns ------- - numpy.ndarray + array The output array with only positive values. """ - return np.maximum(0, z) + xp = array_api_compat.array_namespace(z) + return xp.maximum(0, z) def tanh(z): @@ -34,14 +36,15 @@ def tanh(z): Returns ------- - numpy.ndarray + array The output array with hyperbolic tangent values. See Also -------- numpy.tanh : The hyperbolic tangent function. """ - return np.tanh(z) + xp = array_api_compat.array_namespace(z) + return xp.tanh(z) def sigmoid(z): @@ -55,7 +58,7 @@ def sigmoid(z): Returns ------- - numpy.ndarray + array The output array with the function values. See Also @@ -76,7 +79,7 @@ def identity(z): Returns ------- - numpy.ndarray + array The input array is returned unchanged. """ return z @@ -93,7 +96,7 @@ def softmax(z): Returns ------- - numpy.ndarray + array The output array with the function values. See Also @@ -116,7 +119,7 @@ def softmin(z): Returns ------- - numpy.ndarray + array The output array with the function values. See Also @@ -137,7 +140,7 @@ def log_sigmoid(z): Returns ------- - numpy.ndarray + array The output array with the function values. See Also @@ -158,7 +161,7 @@ def log_softmax(z): Returns ------- - numpy.ndarray + array The output array with the function values. See Also From f5283ca78720bfe3ec7ed0ed7dad30d6879aea0e Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 21 Sep 2026 17:36:06 -0600 Subject: [PATCH 05/32] WIP: zero must be tensor in max func - xp.maximum requires both to be Tensor - 0 as scalar only works for np - numpy not needed in activations.py - use sklearn way to get namespace --- src/gfdl/activations.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/gfdl/activations.py b/src/gfdl/activations.py index 0cf50a3..9192536 100644 --- a/src/gfdl/activations.py +++ b/src/gfdl/activations.py @@ -2,9 +2,8 @@ Activation functions for Gradient Free Deep Learning estimators. """ -import numpy as np import scipy -import array_api_compat +from sklearn.utils._array_api import get_namespace_and_device def relu(z): @@ -21,8 +20,9 @@ def relu(z): array The output array with only positive values. """ - xp = array_api_compat.array_namespace(z) - return xp.maximum(0, z) + xp, _, device = get_namespace_and_device(z) + zero = xp.asarray(0, dtype=z.dtype, device=device) + return xp.maximum(zero, z) def tanh(z): @@ -43,7 +43,7 @@ def tanh(z): -------- numpy.tanh : The hyperbolic tangent function. """ - xp = array_api_compat.array_namespace(z) + xp, _, device = get_namespace_and_device(z) return xp.tanh(z) From 6f0a8f61025c1f8d66fccf4a29b47b39b0fef74a Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Thu, 24 Sep 2026 17:19:16 -0600 Subject: [PATCH 06/32] WIP,ENH: array api dtypes, sklearn tags, and fixtures - Control the dtypes with `xp.astype` - `xp.astype` is used >200 times in sklearn - `move_to` is used over >200 times in sklearn - I.e., sklearn learn uses these for device, dtype - `conftest.py` sets global behavior for array API fixture - sklearn tag tells check estimator that array API is supported - All `test_regression.py` functions pass on local Mac OS - Custom tests are yet to be removed - Seeking alternative existing tests from sklearn to replace --- src/gfdl/model.py | 53 ++++++++++++++++++++++++-- src/gfdl/tests/conftest.py | 14 +++++++ src/gfdl/tests/test_regression.py | 62 ++++++++++++++----------------- 3 files changed, 90 insertions(+), 39 deletions(-) create mode 100644 src/gfdl/tests/conftest.py diff --git a/src/gfdl/model.py b/src/gfdl/model.py index aff1dd8..f33829f 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -18,6 +18,7 @@ from sklearn.utils._array_api import ( get_namespace_and_device, move_to, + _max_precision_float_dtype, ) from sklearn.utils.metaestimators import available_if from sklearn.utils.multiclass import check_classification_targets, unique_labels @@ -47,9 +48,22 @@ def __init__( self.reg_alpha = reg_alpha self.rtol = rtol + def __sklearn_tags__(self): + # Fetch default tags or construct them + tags = super().__sklearn_tags__() + tags.array_api_support = True + return tags + def fit(self, X, Y): xp, _, device = get_namespace_and_device(X) - Y = move_to(Y, xp=xp, device=device) + default_float_dtype = _max_precision_float_dtype(xp, device) + if xp.isdtype(X.dtype, "integral"): + X = xp.astype(X, default_float_dtype) + Y = xp.astype( + move_to(Y, xp=xp, device=device), + X.dtype, + copy=False, + ) # Assumption : X, Y have been pre-processed. # X shape: (n_samples, n_features) @@ -72,37 +86,53 @@ def fit(self, X, Y): rng = self.get_generator(self.seed) self.W_.append( + xp.astype( move_to( self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) ), xp=xp, device=device, + ), + X.dtype, + copy=False, ) ) self.b_.append( + xp.astype( move_to( self._weight_mode(1, hidden_layer_sizes[0], rng=rng) .reshape(-1), xp=xp, device=device, + ), + X.dtype, + copy=False, ) ) for i, layer in enumerate(hidden_layer_sizes[1:]): # (n_hidden, n_features) self.W_.append( + xp.astype( move_to( self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), xp=xp, device=device, + ), + X.dtype, + copy=False, ) ) # (n_hidden,) self.b_.append( + xp.astype( move_to( self._weight_mode(1, layer, rng=rng,).reshape(-1), xp=xp, device=device, + ), + X.dtype, + copy=False, ) ) @@ -253,12 +283,23 @@ def partial_fit(self, X, y): def predict(self, X): check_is_fitted(self) xp, _, device = get_namespace_and_device(X) + default_float_dtype = _max_precision_float_dtype(xp, device) + if xp.isdtype(X.dtype, "integral"): + X = xp.astype(X, default_float_dtype) Hs = [] H_prev = X for W, b in zip(self.W_, self.b_, strict=False): - W = move_to(W, xp=xp, device=device) - b = move_to(b, xp=xp, device=device) + W = xp.astype( + move_to(W, xp=xp, device=device), + X.dtype, + copy=False, + ) + b = xp.astype( + move_to(b, xp=xp, device=device), + X.dtype, + copy=False, + ) Z = H_prev @ W.T + b # (n, m) H_prev = self._activation_fn(Z) Hs.append(H_prev) @@ -266,10 +307,14 @@ def predict(self, X): if self.direct_links: Hs.append(X) D = xp.concat(Hs, axis=1) - out = D @ move_to( + out = D @ xp.astype( + move_to( self.coeff_, xp=xp, device=device, + ), + X.dtype, + copy=False, ) return out diff --git a/src/gfdl/tests/conftest.py b/src/gfdl/tests/conftest.py new file mode 100644 index 0000000..c511f0d --- /dev/null +++ b/src/gfdl/tests/conftest.py @@ -0,0 +1,14 @@ +import os +import pytest +from sklearn import set_config + +def pytest_configure(config): + os.environ["SCIPY_ARRAY_API"] = "1" + +@pytest.fixture(autouse=True) +def enable_array_api(): + # Turn it on before the test runs + set_config(array_api_dispatch=True) + yield + # Reset it back to default after the test finishes + set_config(array_api_dispatch=False) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index a7172e0..1ae0870 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -111,7 +111,8 @@ def test_regression_against_grafo(n_samples, n_targets, hidden_layer_sizes, @parametrize_with_checks([GFDLRegressor()]) def test_sklearn_api_conformance(estimator, check): - check(estimator) + with config_context(array_api_dispatch=False): + check(estimator) @pytest.mark.parametrize("reg_alpha, expected", [ @@ -287,7 +288,6 @@ def test_torch_matches_numpy(reg_alpha, hidden_layer_sizes, ): """NumPy and Torch predictions should be close up to atol""" - os.environ["SCIPY_ARRAY_API"] = "1" estimator = GFDLRegressor report = r2_score rng = np.random.default_rng(seed=42) @@ -304,31 +304,29 @@ def test_torch_matches_numpy(reg_alpha, ) X_torch = torch.asarray(X_np, device="cpu",) y_torch = torch.asarray(y_np, device="cpu",) - with config_context(array_api_dispatch=True): - model = estimator(reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - seed=random_state, - ) - model.fit(X_np, y_np) - y_pred = model.predict(X_torch) - acc_torch = report( - y_torch, - y_pred, - ) - acc_torch_s.append(acc_torch) + model = estimator(reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + seed=random_state, + ) + model.fit(X_np, y_np) + y_pred = model.predict(X_torch) + acc_torch = report( + y_torch, + y_pred, + ) + acc_torch_s.append(acc_torch) # test against numpy which is current usage - with config_context(array_api_dispatch=True): - model = estimator(reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - seed=random_state, - ) - model.fit(X_np, y_np) - y_pred = model.predict(X_np) - acc_np = report( - y_np, - y_pred, - ) - acc_np_s.append(acc_np) + model = estimator(reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + seed=random_state, + ) + model.fit(X_np, y_np) + y_pred = model.predict(X_np) + acc_np = report( + y_np, + y_pred, + ) + acc_np_s.append(acc_np) assert_allclose(acc_torch_s, acc_np_s, atol=1e-3) @@ -348,7 +346,6 @@ def test_predictor_context(namespace, # ValueError thrown .predict has # different namespace, dtype # from its arguments - os.environ["SCIPY_ARRAY_API"] = "1" estimator = GFDLRegressor device = "cpu" @@ -362,8 +359,7 @@ def test_predictor_context(namespace, X1 = namespace.asarray(X, dtype=fit_dtype, device=device) X2 = namespace.asarray(X, dtype=predict_dtype, device=device) y = namespace.asarray(y, dtype=fit_dtype, device=device) - with config_context(array_api_dispatch=True): - y_pred = estimator().fit(X1, y).predict(X2) + y_pred = estimator().fit(X1, y).predict(X2) xp_ypred = array_api_compat.get_namespace(y_pred) xp_X2 = array_api_compat.get_namespace(X2) @@ -390,7 +386,6 @@ def test_int_array_api(namespace, y_dtype, ): """Integer arrays shall be handled gracefully""" - os.environ["SCIPY_ARRAY_API"] = "1" estimator = GFDLRegressor X, y = make_regression( @@ -402,8 +397,7 @@ def test_int_array_api(namespace, ) X = namespace.asarray(X, dtype=X_dtype, device="cpu") y = namespace.asarray(y, dtype=y_dtype, device="cpu") - with config_context(array_api_dispatch=True): - estimator().fit(X, y) + estimator().fit(X, y) @pytest.mark.parametrize( @@ -421,7 +415,6 @@ def test_fit_attr_context(namespace, y_dtype, ): """Fitted attributes shall be same as design matrix""" - os.environ["SCIPY_ARRAY_API"] = "1" estimator = GFDLRegressor device = "cpu" @@ -434,8 +427,7 @@ def test_fit_attr_context(namespace, ) X = namespace.asarray(X, dtype=X_dtype, device=device) y = namespace.asarray(y, dtype=y_dtype, device=device) - with config_context(array_api_dispatch=True): - model = estimator().fit(X, y) + model = estimator().fit(X, y) dtypes = [] devices = [] From 1833bc483be4ec28904f450b063ad09dd05ac6f5 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 28 Sep 2026 14:57:55 -0600 Subject: [PATCH 07/32] WIP,ENH: sklearn-based tests * Added `check_array_api_input_and_values` function. This is a public utility that is not covered by `check(estimator)` and will verify if public attributes match to NumPy runs given flexible choice of array namespace. * I checked many combinations over `GFDLRegressor` inputs except for Torch and NumPy seem to handle backends of bounded weight schemes with slight numerical differences. * `check_array_api_input_and_values` does not have enough flexibility like `rtol` or `atol` to control slight backend differences * Modified `test_sklearn_api_conformance` to check both array API and non-array API version. Made a pytest.xfail for array API when _NotAnArray is put in. This is an sklearn specific object that is incompatible with the array API tests. * Two of my custom tests were already covered by `test_sklearn_api_conformance` with the `array_api_dispatch=True` flag * sklearn.linear_model favors `check_same_namespace` as function to ensure `.predict` input matches the `.fit` input * I liked the idea of converting to what the user puts into `.predict` but I will revert to the standard that sklearn has set * This removes `move_to` instances in `GFDL().predict` * For the time being, I commented out a custom test about integer arrays. sklearn does not address this case. --- src/gfdl/activations.py | 2 +- src/gfdl/model.py | 29 +--- src/gfdl/tests/conftest.py | 3 + src/gfdl/tests/test_regression.py | 269 +++++++++--------------------- 4 files changed, 90 insertions(+), 213 deletions(-) diff --git a/src/gfdl/activations.py b/src/gfdl/activations.py index 9192536..cd143e9 100644 --- a/src/gfdl/activations.py +++ b/src/gfdl/activations.py @@ -43,7 +43,7 @@ def tanh(z): -------- numpy.tanh : The hyperbolic tangent function. """ - xp, _, device = get_namespace_and_device(z) + xp, _, _device = get_namespace_and_device(z) return xp.tanh(z) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index f33829f..d6e5f31 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -16,9 +16,10 @@ from sklearn.preprocessing import OneHotEncoder from sklearn.utils import column_or_1d from sklearn.utils._array_api import ( + _max_precision_float_dtype, + check_same_namespace, get_namespace_and_device, move_to, - _max_precision_float_dtype, ) from sklearn.utils.metaestimators import available_if from sklearn.utils.multiclass import check_classification_targets, unique_labels @@ -158,7 +159,10 @@ def fit(self, X, Y): if self.reg_alpha is None: self.coeff_ = xp.linalg.pinv(D, rtol=self.rtol) @ Y else: - ridge = Ridge(alpha=self.reg_alpha, fit_intercept=False) + ridge = Ridge(alpha=self.reg_alpha, + fit_intercept=False, + solver="svd", + ) ridge.fit(D, Y) self.coeff_ = ridge.coef_.T return self @@ -283,6 +287,7 @@ def partial_fit(self, X, y): def predict(self, X): check_is_fitted(self) xp, _, device = get_namespace_and_device(X) + check_same_namespace(X, self, attribute="coeff_", method="predict") default_float_dtype = _max_precision_float_dtype(xp, device) if xp.isdtype(X.dtype, "integral"): X = xp.astype(X, default_float_dtype) @@ -290,16 +295,6 @@ def predict(self, X): Hs = [] H_prev = X for W, b in zip(self.W_, self.b_, strict=False): - W = xp.astype( - move_to(W, xp=xp, device=device), - X.dtype, - copy=False, - ) - b = xp.astype( - move_to(b, xp=xp, device=device), - X.dtype, - copy=False, - ) Z = H_prev @ W.T + b # (n, m) H_prev = self._activation_fn(Z) Hs.append(H_prev) @@ -307,15 +302,7 @@ def predict(self, X): if self.direct_links: Hs.append(X) D = xp.concat(Hs, axis=1) - out = D @ xp.astype( - move_to( - self.coeff_, - xp=xp, - device=device, - ), - X.dtype, - copy=False, - ) + out = D @ self.coeff_ return out diff --git a/src/gfdl/tests/conftest.py b/src/gfdl/tests/conftest.py index c511f0d..b7e3acc 100644 --- a/src/gfdl/tests/conftest.py +++ b/src/gfdl/tests/conftest.py @@ -1,10 +1,13 @@ import os + import pytest from sklearn import set_config + def pytest_configure(config): os.environ["SCIPY_ARRAY_API"] = "1" + @pytest.fixture(autouse=True) def enable_array_api(): # Turn it on before the test runs diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 1ae0870..cc5250e 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -1,17 +1,16 @@ -import os -import array_api_compat import numpy as np import pytest -import torch -from numpy.testing import assert_allclose from sklearn import config_context from sklearn.base import clone from sklearn.datasets import fetch_openml, make_regression from sklearn.metrics import r2_score from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler -from sklearn.utils.estimator_checks import parametrize_with_checks +from sklearn.utils.estimator_checks import ( + check_array_api_input_and_values, + parametrize_with_checks, +) from gfdl.activations import ACTIVATIONS from gfdl.model import GFDLRegressor @@ -110,9 +109,21 @@ def test_regression_against_grafo(n_samples, n_targets, hidden_layer_sizes, @parametrize_with_checks([GFDLRegressor()]) -def test_sklearn_api_conformance(estimator, check): - with config_context(array_api_dispatch=False): - check(estimator) +@pytest.mark.parametrize("array_api_dispatch", + [False, True], + ids=["dispatch_false", "dispatch_true"], + ) +def test_sklearn_api_conformance(estimator, check, array_api_dispatch): + with config_context(array_api_dispatch=array_api_dispatch): + try: + check(estimator) + except Exception as exc: + if array_api_dispatch and "NotAnArray" in repr(exc): + pytest.xfail( + "Full check_estimator under array_api_dispatch=True is expected " + "to fail on sklearn's NotAnArray sentinel." + ) + raise @pytest.mark.parametrize("reg_alpha, expected", [ @@ -266,191 +277,67 @@ def test_preserve_class_inputs(): assert isinstance(v, type(expected[k])) -@pytest.mark.parametrize( - "reg_alpha, n_features, hidden_layer_sizes", - [ - (1e-1, 40, (100,)), - (1e-1, 400, (100,)), - (2, 40, (100,)), - (2, 400, (100,)), - (None, 40, (100,)), - (None, 400, (100,)), - (1e-1, 40, (100, 100,)), - (1e-1, 400, (100, 100,)), - (2, 40, (100, 100,)), - (2, 400, (100, 100,)), - (None, 40, (100, 100,)), - (None, 400, (100, 100,)), - ] -) -def test_torch_matches_numpy(reg_alpha, - n_features, - hidden_layer_sizes, - ): - """NumPy and Torch predictions should be close up to atol""" - estimator = GFDLRegressor - report = r2_score - rng = np.random.default_rng(seed=42) - acc_np_s = [] - acc_torch_s = [] - for _ in range(10): - random_state = rng.integers(low=0, high=100, size=1)[0] - X_np, y_np = make_regression( - n_samples=10_000, - n_features=n_features, - n_informative=int(n_features / 10), - n_targets=1, - random_state=random_state, - ) - X_torch = torch.asarray(X_np, device="cpu",) - y_torch = torch.asarray(y_np, device="cpu",) - model = estimator(reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - seed=random_state, - ) - model.fit(X_np, y_np) - y_pred = model.predict(X_torch) - acc_torch = report( - y_torch, - y_pred, - ) - acc_torch_s.append(acc_torch) - # test against numpy which is current usage - model = estimator(reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - seed=random_state, - ) - model.fit(X_np, y_np) - y_pred = model.predict(X_np) - acc_np = report( - y_np, - y_pred, - ) - acc_np_s.append(acc_np) - assert_allclose(acc_torch_s, acc_np_s, atol=1e-3) - - -@pytest.mark.parametrize( - "namespace, fit_dtype, predict_dtype", - [ - (np, np.float64, np.float32), - (torch, torch.float64, torch.float32), - (np, np.float32, np.float64), - (torch, torch.float32, torch.float64), - ] -) -def test_predictor_context(namespace, - fit_dtype, - predict_dtype, - ): - # ValueError thrown .predict has - # different namespace, dtype - # from its arguments - estimator = GFDLRegressor - device = "cpu" - - X, y = make_regression( - n_samples=1_000, - n_features=10, - n_informative=4, - n_targets=1, - random_state=42, - ) - X1 = namespace.asarray(X, dtype=fit_dtype, device=device) - X2 = namespace.asarray(X, dtype=predict_dtype, device=device) - y = namespace.asarray(y, dtype=fit_dtype, device=device) - y_pred = estimator().fit(X1, y).predict(X2) - - xp_ypred = array_api_compat.get_namespace(y_pred) - xp_X2 = array_api_compat.get_namespace(X2) - if xp_ypred != xp_X2: - raise ValueError( - ".predict output and input are not in the same namespace" - ) - if predict_dtype != y_pred.dtype: - raise ValueError( - ".predict output and input are not the same dtype" - ) - - -@pytest.mark.parametrize( - "namespace, X_dtype, y_dtype", - [ - (np, np.float64, np.int64), - (np, np.int64, np.float64), - (torch, torch.float64, torch.int64), - ] -) -def test_int_array_api(namespace, - X_dtype, - y_dtype, - ): - """Integer arrays shall be handled gracefully""" - estimator = GFDLRegressor - - X, y = make_regression( - n_samples=1_000, - n_features=10, - n_informative=4, - n_targets=1, - random_state=42, - ) - X = namespace.asarray(X, dtype=X_dtype, device="cpu") - y = namespace.asarray(y, dtype=y_dtype, device="cpu") - estimator().fit(X, y) - - -@pytest.mark.parametrize( - "namespace, X_dtype, y_dtype", - [ - (np, np.float64, np.float64), - (np, np.float32, np.float32), - (torch, torch.float64, torch.float64), - (torch, torch.float32, torch.float32), - (torch, torch.float64, torch.float32), - ] -) -def test_fit_attr_context(namespace, - X_dtype, - y_dtype, +@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("direct_links", [True, False]) +@pytest.mark.parametrize("activation", ["identity", + "relu", + "tanh", + "log_sigmoid", + "softmax", + "softmin", + "sigmoid", + ]) +@pytest.mark.parametrize("weight_scheme", ["normal", + "zeros", + ]) +@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) +@pytest.mark.parametrize("namespace", ["numpy", "torch"]) +def test_array_api_values(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + namespace, ): - """Fitted attributes shall be same as design matrix""" - estimator = GFDLRegressor - device = "cpu" + """NumPy and Array API values and predictions should be close""" - X, y = make_regression( - n_samples=1_000, - n_features=10, - n_informative=4, - n_targets=1, - random_state=42, + estimator = GFDLRegressor( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=42, ) - X = namespace.asarray(X, dtype=X_dtype, device=device) - y = namespace.asarray(y, dtype=y_dtype, device=device) - model = estimator().fit(X, y) - - dtypes = [] - devices = [] - spaces = [] - W_ = model.W_ - b_ = model.b_ - for w, b in zip(W_, b_, strict=False): - dtypes.extend((b.dtype, w.dtype)) - devices.extend((b.device, w.device)) - spaces.extend( - (array_api_compat.get_namespace(w), - array_api_compat.get_namespace(b),) - ) - - coeff_ = model.coeff_ - dtypes.append(coeff_.dtype) - devices.append(coeff_.device) - spaces.append( - array_api_compat.get_namespace(coeff_) + check_array_api_input_and_values( + estimator.__class__.__name__, + estimator, + array_namespace=namespace, + device_name="cpu", ) - assert all(dtype == X.dtype for dtype in dtypes) - assert all(device == X.device for device in devices) - - X_space = array_api_compat.get_namespace(X) - assert all(space == X_space for space in spaces) +# @pytest.mark.parametrize( +# "namespace, X_dtype, y_dtype", +# [ +# (np, np.float64, np.int64), +# (np, np.int64, np.float64), +# (torch, torch.float64, torch.int64), +# ] +# ) +# def test_int_array_api(namespace, +# X_dtype, +# y_dtype, +# ): +# """Integer arrays shall be handled gracefully""" +# estimator = GFDLRegressor + +# X, y = make_regression( +# n_samples=1_000, +# n_features=10, +# n_informative=4, +# n_targets=1, +# random_state=42, +# ) +# X = namespace.asarray(X, dtype=X_dtype, device="cpu") +# y = namespace.asarray(y, dtype=y_dtype, device="cpu") +# estimator().fit(X, y) From 553ac1aa2bf32ec73540abf01f7b1ecddd323892 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 28 Sep 2026 16:48:20 -0600 Subject: [PATCH 08/32] WIP,ENH: array api checks local to some tests * Removed conftest.py file which globally set some array API settings * Reverted back to original check_estimator test * Made separate check_estimator test for array API Includes skips if array API not set up * Use pytest.importorskip for "torch" test * array_api_strict bawks at hidden_layer_sizes.min() * Removed array_api_compat dependency --- pyproject.toml | 1 - src/gfdl/model.py | 2 +- src/gfdl/tests/conftest.py | 17 ------ src/gfdl/tests/test_regression.py | 99 ++++++++++++++----------------- 4 files changed, 45 insertions(+), 74 deletions(-) delete mode 100644 src/gfdl/tests/conftest.py diff --git a/pyproject.toml b/pyproject.toml index 294b784..49be48a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,7 +28,6 @@ classifiers = [ dependencies = ["numpy>=2.0.0", "scikit-learn>=1.5.0", "scipy>=1.13.0", - "array-api-compat", "packaging>=24.0"] [project.urls] diff --git a/src/gfdl/model.py b/src/gfdl/model.py index d6e5f31..03132e5 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -72,7 +72,7 @@ def fit(self, X, Y): if self.reg_alpha is not None and self.reg_alpha < 0.0: raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).") hidden_layer_sizes = xp.asarray(self.hidden_layer_sizes) - if hidden_layer_sizes.min() < 1: + if min(hidden_layer_sizes) < 1: raise ValueError("hidden_layer_sizes must be > 0, " f"got {hidden_layer_sizes}") fn = resolve_activation(self.activation)[1] diff --git a/src/gfdl/tests/conftest.py b/src/gfdl/tests/conftest.py deleted file mode 100644 index b7e3acc..0000000 --- a/src/gfdl/tests/conftest.py +++ /dev/null @@ -1,17 +0,0 @@ -import os - -import pytest -from sklearn import set_config - - -def pytest_configure(config): - os.environ["SCIPY_ARRAY_API"] = "1" - - -@pytest.fixture(autouse=True) -def enable_array_api(): - # Turn it on before the test runs - set_config(array_api_dispatch=True) - yield - # Reset it back to default after the test finishes - set_config(array_api_dispatch=False) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index cc5250e..c790195 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -1,4 +1,6 @@ +import os + import numpy as np import pytest from sklearn import config_context @@ -7,6 +9,9 @@ from sklearn.metrics import r2_score from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler +from sklearn.utils._testing import ( + skip_if_array_api_compat_not_configured, +) from sklearn.utils.estimator_checks import ( check_array_api_input_and_values, parametrize_with_checks, @@ -109,21 +114,8 @@ def test_regression_against_grafo(n_samples, n_targets, hidden_layer_sizes, @parametrize_with_checks([GFDLRegressor()]) -@pytest.mark.parametrize("array_api_dispatch", - [False, True], - ids=["dispatch_false", "dispatch_true"], - ) -def test_sklearn_api_conformance(estimator, check, array_api_dispatch): - with config_context(array_api_dispatch=array_api_dispatch): - try: - check(estimator) - except Exception as exc: - if array_api_dispatch and "NotAnArray" in repr(exc): - pytest.xfail( - "Full check_estimator under array_api_dispatch=True is expected " - "to fail on sklearn's NotAnArray sentinel." - ) - raise +def test_sklearn_api_conformance(estimator, check): + check(estimator) @pytest.mark.parametrize("reg_alpha, expected", [ @@ -277,6 +269,28 @@ def test_preserve_class_inputs(): assert isinstance(v, type(expected[k])) +@skip_if_array_api_compat_not_configured +@pytest.mark.skipif( + os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." +) +@parametrize_with_checks([GFDLRegressor()]) +def test_sklearn_array_api_conformance(estimator, check,): + with config_context(array_api_dispatch=True): + try: + check(estimator) + except Exception as exc: + if "NotAnArray" in repr(exc): + pytest.xfail( + "Full check_estimator under array_api_dispatch=True is expected " + "to fail on sklearn's NotAnArray sentinel." + ) + raise + + +@skip_if_array_api_compat_not_configured +@pytest.mark.skipif( + os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." +) @pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ["identity", @@ -291,16 +305,16 @@ def test_preserve_class_inputs(): "zeros", ]) @pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -@pytest.mark.parametrize("namespace", ["numpy", "torch"]) -def test_array_api_values(hidden_layer_sizes, - direct_links, - activation, - weight_scheme, - reg_alpha, - namespace, - ): +def test_torch_array_api(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + ): """NumPy and Array API values and predictions should be close""" + pytest.importorskip("torch") + estimator = GFDLRegressor( reg_alpha=reg_alpha, hidden_layer_sizes=hidden_layer_sizes, @@ -309,35 +323,10 @@ def test_array_api_values(hidden_layer_sizes, direct_links=direct_links, seed=42, ) - check_array_api_input_and_values( - estimator.__class__.__name__, - estimator, - array_namespace=namespace, - device_name="cpu", - ) - -# @pytest.mark.parametrize( -# "namespace, X_dtype, y_dtype", -# [ -# (np, np.float64, np.int64), -# (np, np.int64, np.float64), -# (torch, torch.float64, torch.int64), -# ] -# ) -# def test_int_array_api(namespace, -# X_dtype, -# y_dtype, -# ): -# """Integer arrays shall be handled gracefully""" -# estimator = GFDLRegressor - -# X, y = make_regression( -# n_samples=1_000, -# n_features=10, -# n_informative=4, -# n_targets=1, -# random_state=42, -# ) -# X = namespace.asarray(X, dtype=X_dtype, device="cpu") -# y = namespace.asarray(y, dtype=y_dtype, device="cpu") -# estimator().fit(X, y) + with config_context(array_api_dispatch=True): + check_array_api_input_and_values( + estimator.__class__.__name__, + estimator, + array_namespace="torch", + device_name="cpu", + ) From c78321f4b8bae4f9c6e0027ff57741ea67ac21f5 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 28 Sep 2026 18:57:10 -0600 Subject: [PATCH 09/32] WIP,ENH: make_regression test for torch * `check_array_api_input_and_values` uses `make_classification` which is unstable for some weight schemes between torch, np * Shrunk number of cases for that test which calls `check_array_api_input_and_values` * Added custom `make_regression`-based test that checks public attributes and predictions --- src/gfdl/tests/test_regression.py | 76 +++++++++++++++++++++++++++++-- 1 file changed, 71 insertions(+), 5 deletions(-) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index c790195..f82efc5 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -3,6 +3,7 @@ import numpy as np import pytest +from numpy.testing import assert_allclose from sklearn import config_context from sklearn.base import clone from sklearn.datasets import fetch_openml, make_regression @@ -296,10 +297,6 @@ def test_sklearn_array_api_conformance(estimator, check,): @pytest.mark.parametrize("activation", ["identity", "relu", "tanh", - "log_sigmoid", - "softmax", - "softmin", - "sigmoid", ]) @pytest.mark.parametrize("weight_scheme", ["normal", "zeros", @@ -311,7 +308,7 @@ def test_torch_array_api(hidden_layer_sizes, weight_scheme, reg_alpha, ): - """NumPy and Array API values and predictions should be close""" + """NumPy and Array API public attributes and predictions should be close""" pytest.importorskip("torch") @@ -324,9 +321,78 @@ def test_torch_array_api(hidden_layer_sizes, seed=42, ) with config_context(array_api_dispatch=True): + # Simulates from X, y = make_classification(n_samples=30, n_features=10) + # Unstable across backends for some weight schemes check_array_api_input_and_values( estimator.__class__.__name__, estimator, array_namespace="torch", device_name="cpu", ) + + +@skip_if_array_api_compat_not_configured +@pytest.mark.skipif( + os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." +) +@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("direct_links", [True, False]) +@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) +@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) +@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) +def test_torch_make_regression(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + ): + """Custom make_regression()-based NumPy and PyTorch matching""" + pytest.importorskip("torch") + import torch + estimator = GFDLRegressor + rng = np.random.default_rng(seed=42) + + pred_torch = [] + pred_numpy = [] + coef_torch = [] + coef_numpy = [] + for _ in range(3): + random_state = rng.integers(low=0, high=100, size=1)[0] + X_np, y_np = make_regression( + n_samples=250, + n_features=20, + n_informative=5, + n_targets=1, + random_state=random_state, + ) + + # NumPy .fit and .predict + model = estimator( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=random_state,) + model.fit(X_np, y_np) + y_pred_np = model.predict(X_np) + pred_numpy.append(y_pred_np) + coef_numpy.append(model.coeff_) + + # PyTorch .fit and .predict + X_torch = torch.asarray(X_np, device="cpu",) + y_torch = torch.asarray(y_np, device="cpu",) + model = estimator( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=random_state,) + model.fit(X_torch, y_torch) + y_pred_torch = model.predict(X_torch) + pred_torch.append(y_pred_torch) + coef_torch.append(model.coeff_) + + assert_allclose(pred_numpy, pred_torch) + assert_allclose(coef_numpy, coef_torch) From 45b7090baf6aac8c459ed15158807189cfdd6892 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 28 Sep 2026 21:32:40 -0600 Subject: [PATCH 10/32] WIP,MAINT: xp.min instead of min or .min() * Applied to `hidden_layer_sizes` --- src/gfdl/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 03132e5..792ba54 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -72,7 +72,7 @@ def fit(self, X, Y): if self.reg_alpha is not None and self.reg_alpha < 0.0: raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).") hidden_layer_sizes = xp.asarray(self.hidden_layer_sizes) - if min(hidden_layer_sizes) < 1: + if xp.min(hidden_layer_sizes) < 1: raise ValueError("hidden_layer_sizes must be > 0, " f"got {hidden_layer_sizes}") fn = resolve_activation(self.activation)[1] From 4112b4323bd9114c396f117f74f2ac11ffc131f8 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Tue, 29 Sep 2026 10:22:10 -0600 Subject: [PATCH 11/32] WIP,ENH,TST: xp flexible tests for check values * Implemented argument array_namespace instead of torch only * Used `clone()` in custom test to cut redundancy. This is what `sklearn` does in the related test `_check_array_api_core`. * Removed the 3x for loop of different random states. * The custom `xp_make_regression` test remains but ideally it could be removed if we could get an MR change through `sklearn` that allows for control of the simulating data. --- src/gfdl/tests/test_regression.py | 113 ++++++++++++++++-------------- 1 file changed, 59 insertions(+), 54 deletions(-) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index f82efc5..d280a6f 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -302,15 +302,22 @@ def test_sklearn_array_api_conformance(estimator, check,): "zeros", ]) @pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -def test_torch_array_api(hidden_layer_sizes, +@pytest.mark.parametrize("array_namespace", + [ + "numpy", + "torch" + ], + ) +def test_xp_check_values(hidden_layer_sizes, direct_links, activation, weight_scheme, reg_alpha, + array_namespace, ): """NumPy and Array API public attributes and predictions should be close""" - pytest.importorskip("torch") + pytest.importorskip(array_namespace) estimator = GFDLRegressor( reg_alpha=reg_alpha, @@ -326,7 +333,7 @@ def test_torch_array_api(hidden_layer_sizes, check_array_api_input_and_values( estimator.__class__.__name__, estimator, - array_namespace="torch", + array_namespace=array_namespace, device_name="cpu", ) @@ -340,59 +347,57 @@ def test_torch_array_api(hidden_layer_sizes, @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -def test_torch_make_regression(hidden_layer_sizes, - direct_links, - activation, - weight_scheme, - reg_alpha, - ): - """Custom make_regression()-based NumPy and PyTorch matching""" - pytest.importorskip("torch") - import torch +@pytest.mark.parametrize("array_namespace", + [ + "numpy", + "torch" + ], + ) +def test_xp_make_regression(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + array_namespace, + ): + """Custom make_regression()-based NumPy and array API matching""" + xp = pytest.importorskip(array_namespace) estimator = GFDLRegressor - rng = np.random.default_rng(seed=42) - pred_torch = [] + pred_xp = [] pred_numpy = [] - coef_torch = [] + coef_xp = [] coef_numpy = [] - for _ in range(3): - random_state = rng.integers(low=0, high=100, size=1)[0] - X_np, y_np = make_regression( - n_samples=250, - n_features=20, - n_informative=5, - n_targets=1, - random_state=random_state, - ) + random_state = 42 + X_np, y_np = make_regression( + n_samples=250, + n_features=20, + n_informative=5, + n_targets=1, + random_state=random_state, + ) - # NumPy .fit and .predict - model = estimator( - reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - activation=activation, - weight_scheme=weight_scheme, - direct_links=direct_links, - seed=random_state,) - model.fit(X_np, y_np) - y_pred_np = model.predict(X_np) - pred_numpy.append(y_pred_np) - coef_numpy.append(model.coeff_) - - # PyTorch .fit and .predict - X_torch = torch.asarray(X_np, device="cpu",) - y_torch = torch.asarray(y_np, device="cpu",) - model = estimator( - reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - activation=activation, - weight_scheme=weight_scheme, - direct_links=direct_links, - seed=random_state,) - model.fit(X_torch, y_torch) - y_pred_torch = model.predict(X_torch) - pred_torch.append(y_pred_torch) - coef_torch.append(model.coeff_) - - assert_allclose(pred_numpy, pred_torch) - assert_allclose(coef_numpy, coef_torch) + # NumPy .fit and .predict + model = estimator( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=random_state,) + model.fit(X_np, y_np) + y_pred_np = model.predict(X_np) + pred_numpy.append(y_pred_np) + coef_numpy.append(model.coeff_) + + # xp (array API) .fit and .predict + X_xp = xp.asarray(X_np, device="cpu",) + y_xp = xp.asarray(y_np, device="cpu",) + model_xp = clone(model) + model_xp.fit(X_xp, y_xp) + y_pred_xp = model_xp.predict(X_xp) + pred_xp.append(y_pred_xp) + coef_xp.append(model_xp.coeff_) + + assert_allclose(pred_numpy, pred_xp) + assert_allclose(coef_numpy, coef_xp) From 8a05709c89e40320abf76dbc4f181524a164254e Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Tue, 29 Sep 2026 12:47:31 -0600 Subject: [PATCH 12/32] WIP,ENH,TST,BUG: need with context, atol and rtol control * In `test_xp_make_regression` we need to use `with config_context()` otherwise `sklearn` utilities revert to NumPy. This was a bug. * No handling of dtype for the time being. All tests are xp.float64. * Needed to control `atol` and `rtol` * `move_to` to `assert_allclose` on same namespace (NumPy) * Small code change to use the array_api_compat vendor of the other namespaces (PyTorch) --- src/gfdl/model.py | 6 ++-- src/gfdl/tests/test_regression.py | 47 ++++++++++++++++++------------- 2 files changed, 29 insertions(+), 24 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 792ba54..c2fd95f 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -58,8 +58,7 @@ def __sklearn_tags__(self): def fit(self, X, Y): xp, _, device = get_namespace_and_device(X) default_float_dtype = _max_precision_float_dtype(xp, device) - if xp.isdtype(X.dtype, "integral"): - X = xp.astype(X, default_float_dtype) + X = xp.astype(X, default_float_dtype) Y = xp.astype( move_to(Y, xp=xp, device=device), X.dtype, @@ -289,8 +288,7 @@ def predict(self, X): xp, _, device = get_namespace_and_device(X) check_same_namespace(X, self, attribute="coeff_", method="predict") default_float_dtype = _max_precision_float_dtype(xp, device) - if xp.isdtype(X.dtype, "integral"): - X = xp.astype(X, default_float_dtype) + X = xp.astype(X, default_float_dtype) Hs = [] H_prev = X diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index d280a6f..faf266c 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -10,6 +10,11 @@ from sklearn.metrics import r2_score from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler +from sklearn.utils._array_api import ( + get_namespace, + get_namespace_and_device, + move_to, +) from sklearn.utils._testing import ( skip_if_array_api_compat_not_configured, ) @@ -342,7 +347,7 @@ def test_xp_check_values(hidden_layer_sizes, @pytest.mark.skipif( os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." ) -@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @@ -362,12 +367,11 @@ def test_xp_make_regression(hidden_layer_sizes, ): """Custom make_regression()-based NumPy and array API matching""" xp = pytest.importorskip(array_namespace) + with config_context(array_api_dispatch=True): + xp, _ = get_namespace(xp.zeros(2)) + estimator = GFDLRegressor - pred_xp = [] - pred_numpy = [] - coef_xp = [] - coef_numpy = [] random_state = 42 X_np, y_np = make_regression( n_samples=250, @@ -378,26 +382,29 @@ def test_xp_make_regression(hidden_layer_sizes, ) # NumPy .fit and .predict - model = estimator( + model_np = estimator( reg_alpha=reg_alpha, hidden_layer_sizes=hidden_layer_sizes, activation=activation, weight_scheme=weight_scheme, direct_links=direct_links, seed=random_state,) - model.fit(X_np, y_np) - y_pred_np = model.predict(X_np) - pred_numpy.append(y_pred_np) - coef_numpy.append(model.coeff_) + with config_context(array_api_dispatch=False): + model_np.fit(X_np, y_np) + y_pred_np = model_np.predict(X_np) + xp_np, _, device_np = get_namespace_and_device( + X_np + ) # xp (array API) .fit and .predict - X_xp = xp.asarray(X_np, device="cpu",) - y_xp = xp.asarray(y_np, device="cpu",) - model_xp = clone(model) - model_xp.fit(X_xp, y_xp) - y_pred_xp = model_xp.predict(X_xp) - pred_xp.append(y_pred_xp) - coef_xp.append(model_xp.coeff_) - - assert_allclose(pred_numpy, pred_xp) - assert_allclose(coef_numpy, coef_xp) + model_xp = clone(model_np) + with config_context(array_api_dispatch=True): + X_xp = xp.asarray(X_np,) + y_xp = xp.asarray(y_np,) + model_xp.fit(X_xp, y_xp) + y_pred_xp = model_xp.predict(X_xp) + + y_pred_xp = move_to(y_pred_xp, xp=xp_np, device=device_np) + model_xp_coeff_ = move_to(model_xp.coeff_, xp=xp_np, device=device_np) + assert_allclose(y_pred_np, y_pred_xp, atol=1e-5, rtol=1e-5) + assert_allclose(model_np.coeff_, model_xp_coeff_, atol=1e-5, rtol=1e-5) From 403928c36432cb3c0376d61a36f7d7d235202102 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 12 Oct 2026 17:56:44 -0600 Subject: [PATCH 13/32] WIP,ENH,TST: use yield for namespace, device, dtype * `check_array_api_input` only * Values are evaluated with custom test that is based on `check_array_api_core` but gives control over simulating data, `atol` and `rtol` * Removed `pytest.importorskip("torch")` * Instead, use the utility from sklearn `yield_namespace_device_dtype_combinations()` * `if` block to scope to torch, cpu, float64 * Don't force to max dtype precision * All `WEIGHTS.keys()` and `ACTIVATIONS.keys()` --- src/gfdl/model.py | 10 +- src/gfdl/tests/test_regression.py | 157 +++++++++++++++++------------- 2 files changed, 97 insertions(+), 70 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index c2fd95f..66b7817 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -57,8 +57,9 @@ def __sklearn_tags__(self): def fit(self, X, Y): xp, _, device = get_namespace_and_device(X) - default_float_dtype = _max_precision_float_dtype(xp, device) - X = xp.astype(X, default_float_dtype) + if not xp.isdtype(X.dtype, "real floating"): + default_float_dtype = _max_precision_float_dtype(xp, device) + X = xp.astype(X, default_float_dtype) Y = xp.astype( move_to(Y, xp=xp, device=device), X.dtype, @@ -287,8 +288,9 @@ def predict(self, X): check_is_fitted(self) xp, _, device = get_namespace_and_device(X) check_same_namespace(X, self, attribute="coeff_", method="predict") - default_float_dtype = _max_precision_float_dtype(xp, device) - X = xp.astype(X, default_float_dtype) + if not xp.isdtype(X.dtype, "real floating"): + default_float_dtype = _max_precision_float_dtype(xp, device) + X = xp.astype(X, default_float_dtype) Hs = [] H_prev = X diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index faf266c..34e0c16 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -1,5 +1,4 @@ -import os import numpy as np import pytest @@ -11,15 +10,15 @@ from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.utils._array_api import ( - get_namespace, - get_namespace_and_device, move_to, + yield_namespace_device_dtype_combinations, ) from sklearn.utils._testing import ( + _array_api_for_tests, skip_if_array_api_compat_not_configured, ) from sklearn.utils.estimator_checks import ( - check_array_api_input_and_values, + check_array_api_input, parametrize_with_checks, ) @@ -276,9 +275,6 @@ def test_preserve_class_inputs(): @skip_if_array_api_compat_not_configured -@pytest.mark.skipif( - os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." -) @parametrize_with_checks([GFDLRegressor()]) def test_sklearn_array_api_conformance(estimator, check,): with config_context(array_api_dispatch=True): @@ -294,35 +290,30 @@ def test_sklearn_array_api_conformance(estimator, check,): @skip_if_array_api_compat_not_configured -@pytest.mark.skipif( - os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." -) -@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) @pytest.mark.parametrize("direct_links", [True, False]) -@pytest.mark.parametrize("activation", ["identity", - "relu", - "tanh", - ]) -@pytest.mark.parametrize("weight_scheme", ["normal", - "zeros", - ]) -@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -@pytest.mark.parametrize("array_namespace", - [ - "numpy", - "torch" - ], - ) -def test_xp_check_values(hidden_layer_sizes, +@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) +@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) +@pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize( + "array_namespace, device_name, dtype_name", + yield_namespace_device_dtype_combinations(), +) +def test_xp_check_inputs(hidden_layer_sizes, direct_links, activation, weight_scheme, reg_alpha, array_namespace, + device_name, + dtype_name, ): - """NumPy and Array API public attributes and predictions should be close""" - - pytest.importorskip(array_namespace) + # NumPy and Array API public attributes and predictions should be close + if (array_namespace not in ["numpy", "torch"] or + device_name not in [None, "cpu"] or + dtype_name not in [None, "float64"] + ): + pytest.skip("unsupported configuration") estimator = GFDLRegressor( reg_alpha=reg_alpha, @@ -335,54 +326,67 @@ def test_xp_check_values(hidden_layer_sizes, with config_context(array_api_dispatch=True): # Simulates from X, y = make_classification(n_samples=30, n_features=10) # Unstable across backends for some weight schemes - check_array_api_input_and_values( + check_array_api_input( estimator.__class__.__name__, estimator, array_namespace=array_namespace, - device_name="cpu", + device_name=device_name, ) @skip_if_array_api_compat_not_configured -@pytest.mark.skipif( - os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." -) @pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -@pytest.mark.parametrize("array_namespace", - [ - "numpy", - "torch" - ], - ) +@pytest.mark.parametrize( + "array_namespace, device_name, dtype_name", + yield_namespace_device_dtype_combinations(), +) def test_xp_make_regression(hidden_layer_sizes, direct_links, activation, weight_scheme, reg_alpha, array_namespace, + device_name, + dtype_name, ): - """Custom make_regression()-based NumPy and array API matching""" - xp = pytest.importorskip(array_namespace) - with config_context(array_api_dispatch=True): - xp, _ = get_namespace(xp.zeros(2)) - + # Custom make_regression()-based NumPy and array API matching + # This test is closely based on _check_array_api_core from sklearn. + # Here, we have better control of `atol`, `rtol`, and simulating data. + # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117 + # Other Array API estimators like LogisticRegression also write custom tests to + # control `atol`, `rtol`, and simulating data. + # See: `test_logistic_regression_array_api_compliance` function + # `test_xp_check_inputs` uses sklearn utility already to check + # flanking test infracture. + if (array_namespace not in ["numpy", "torch"] or + device_name not in [None, "cpu"] or + dtype_name not in [None, "float64"] + ): + pytest.skip("unsupported configuration") + xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) + + atol = 1e-3 + rtol = 1e-6 + n_samples = 500 + random_state = 42 estimator = GFDLRegressor - random_state = 42 - X_np, y_np = make_regression( - n_samples=250, + X, y = make_regression( + n_samples=n_samples, n_features=20, n_informative=5, n_targets=1, random_state=random_state, ) + X = X.astype(dtype_name, copy=False) + X_xp = xp.asarray(X, device=device) - # NumPy .fit and .predict - model_np = estimator( + # NumPy .fit + est = estimator( reg_alpha=reg_alpha, hidden_layer_sizes=hidden_layer_sizes, activation=activation, @@ -390,21 +394,42 @@ def test_xp_make_regression(hidden_layer_sizes, direct_links=direct_links, seed=random_state,) with config_context(array_api_dispatch=False): - model_np.fit(X_np, y_np) - y_pred_np = model_np.predict(X_np) - xp_np, _, device_np = get_namespace_and_device( - X_np - ) + est.fit(X, y) - # xp (array API) .fit and .predict - model_xp = clone(model_np) + # xp (array API) .fit + est_xp = clone(est) with config_context(array_api_dispatch=True): - X_xp = xp.asarray(X_np,) - y_xp = xp.asarray(y_np,) - model_xp.fit(X_xp, y_xp) - y_pred_xp = model_xp.predict(X_xp) - - y_pred_xp = move_to(y_pred_xp, xp=xp_np, device=device_np) - model_xp_coeff_ = move_to(model_xp.coeff_, xp=xp_np, device=device_np) - assert_allclose(y_pred_np, y_pred_xp, atol=1e-5, rtol=1e-5) - assert_allclose(model_np.coeff_, model_xp_coeff_, atol=1e-5, rtol=1e-5) + est_xp.fit(X_xp, y) + + # Check fitted attributes are close + array_attributes = { + key: value for key, value in vars(est).items() if isinstance(value, np.ndarray) + } + for attribute_name, attribute_value in array_attributes.items(): + est_xp_attr = getattr(est_xp, attribute_name) + est_xp_attr_np = move_to(est_xp_attr, xp=np, device="cpu") + assert_allclose( + attribute_value, + est_xp_attr_np, + err_msg=f"{attribute_name} not the same", + atol=atol, + rtol=rtol, + ) + + # Check method outputs are close + methods = [ + "predict", + ] + for method_name in methods: + with config_context(array_api_dispatch=False): + result = getattr(est, method_name)(X) + with config_context(array_api_dispatch=True): + result_xp = getattr(est_xp, method_name)(X_xp) + result_xp_np = move_to(result_xp, xp=np, device="cpu") + assert_allclose( + result, + result_xp_np, + err_msg=f"{method_name}() did not the return the same result", + atol=atol, + rtol=rtol, + ) From fd5b32a1de233f010e7e82053bf0a71180cbf63a Mon Sep 17 00:00:00 2001 From: sdtemple Date: Wed, 14 Oct 2026 15:26:02 -0600 Subject: [PATCH 14/32] WIP: revise for strict Array API * AI disclosure: I consulted AI but wrote the small edits myself * ridge.coef_.T not strict compliant for 1-d arrays. Revised to not call .T if the array is 1-d * Revert from xp.asarray to np.asarray for hidden_layer_sizes. GaussianNB and others in sklearn have precedent for not all arrays being in xp. This just loops over layer sizes. In experimentation, I find that if this is xp.asarray it creates issues for weight sampling in cupy. No diff from main. * Shortened and revised some comments * Shrunk num of pytest marks by 2 * Added dtype_name to sklearn utility invocation --- src/gfdl/model.py | 7 ++++--- src/gfdl/tests/test_regression.py | 21 +++++++-------------- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 66b7817..fb96545 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -71,8 +71,8 @@ def fit(self, X, Y): # Y shape: (n_samples, n_classes-1) if self.reg_alpha is not None and self.reg_alpha < 0.0: raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).") - hidden_layer_sizes = xp.asarray(self.hidden_layer_sizes) - if xp.min(hidden_layer_sizes) < 1: + hidden_layer_sizes = np.asarray(self.hidden_layer_sizes) + if hidden_layer_sizes.min() < 1: raise ValueError("hidden_layer_sizes must be > 0, " f"got {hidden_layer_sizes}") fn = resolve_activation(self.activation)[1] @@ -164,7 +164,8 @@ def fit(self, X, Y): solver="svd", ) ridge.fit(D, Y) - self.coeff_ = ridge.coef_.T + coef = ridge.coef_ + self.coeff_ = coef.T if coef.ndim == 2 else coef return self def partial_fit(self, X, y): diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 34e0c16..d937aa4 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -281,16 +281,16 @@ def test_sklearn_array_api_conformance(estimator, check,): try: check(estimator) except Exception as exc: - if "NotAnArray" in repr(exc): + if "_NotAnArray" in repr(exc): pytest.xfail( "Full check_estimator under array_api_dispatch=True is expected " - "to fail on sklearn's NotAnArray sentinel." + "to fail on sklearn's _NotAnArray sentinel." ) raise @skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @@ -308,7 +308,6 @@ def test_xp_check_inputs(hidden_layer_sizes, device_name, dtype_name, ): - # NumPy and Array API public attributes and predictions should be close if (array_namespace not in ["numpy", "torch"] or device_name not in [None, "cpu"] or dtype_name not in [None, "float64"] @@ -324,13 +323,12 @@ def test_xp_check_inputs(hidden_layer_sizes, seed=42, ) with config_context(array_api_dispatch=True): - # Simulates from X, y = make_classification(n_samples=30, n_features=10) - # Unstable across backends for some weight schemes check_array_api_input( estimator.__class__.__name__, estimator, array_namespace=array_namespace, device_name=device_name, + dtype_name=dtype_name, ) @@ -339,7 +337,7 @@ def test_xp_check_inputs(hidden_layer_sizes, @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) -@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) +@pytest.mark.parametrize("reg_alpha", [None, 0.1]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), @@ -354,14 +352,9 @@ def test_xp_make_regression(hidden_layer_sizes, dtype_name, ): # Custom make_regression()-based NumPy and array API matching - # This test is closely based on _check_array_api_core from sklearn. + # This test is closely based on `_check_array_api_core` from sklearn. # Here, we have better control of `atol`, `rtol`, and simulating data. - # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117 - # Other Array API estimators like LogisticRegression also write custom tests to - # control `atol`, `rtol`, and simulating data. - # See: `test_logistic_regression_array_api_compliance` function - # `test_xp_check_inputs` uses sklearn utility already to check - # flanking test infracture. + # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. if (array_namespace not in ["numpy", "torch"] or device_name not in [None, "cpu"] or dtype_name not in [None, "float64"] From 0f0a168187596a6da41d642c69ccd91238c2fe2b Mon Sep 17 00:00:00 2001 From: sdtemple Date: Wed, 14 Oct 2026 17:46:25 -0600 Subject: [PATCH 15/32] WIP,TST: no check array api for classifier * AI disclosure: I did not use AI for this * Skip on `check_array_api` for classifier tests * We will remove those skips once we address the classifier in PR #130 --- src/gfdl/tests/test_model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 66ce1af..3225d41 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -436,7 +436,9 @@ def test_classification_against_grafo(hidden_layer_sizes, n_classes, activation, @parametrize_with_checks([GFDLClassifier(), EnsembleGFDLClassifier()]) -def test_sklearn_api_conformance(estimator, check): +def test_sklearn_api_conformance(estimator, check, request): + if "check_array_api" in request.node.name: + pytest.skip("Array API checks run separately") check(estimator) From ac565e59972f1f353614b94ba285b7514bb0da5a Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 19 Oct 2026 12:35:22 -0600 Subject: [PATCH 16/32] WIP,MAINT: CI for array API * SI disclosure: I did not use SI for this. This is my first set up of a CI, so hopefully it works smoothly. * Added new YML file with SCIPY_ARRAY_API=1 as the environment file * Only macos and ubuntu and python 3.14 * Only torch and array-api-strict installed --- .github/workflows/array-api.yml | 34 +++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 .github/workflows/array-api.yml diff --git a/.github/workflows/array-api.yml b/.github/workflows/array-api.yml new file mode 100644 index 0000000..700278d --- /dev/null +++ b/.github/workflows/array-api.yml @@ -0,0 +1,34 @@ +name: Array API CI + +on: + push: + branches: + - maintenance/** + pull_request: + branches: + - main + - maintenance/** + +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, macos-latest] + python-version: ["3.14"] + + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install -v . + python -m pip install -v --group test + python -m pip install torch array-api-strict array-api-compat + - name: test + run: | + cd $RUNNER_TEMP + SCIPY_ARRAY_API=1 python -m pytest --pyargs gfdl --cov=gfdl --cov-report=term-missing From b5672ab787e9099b9ccf9d8154b5db329d67ec08 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 19 Oct 2026 17:08:42 -0600 Subject: [PATCH 17/32] WIP,MAINT: use xp for `sigmoid`, `log_sigmoid` * Replace SciPy version of `log_sigmoid` and `sigmoid` with an `xp` version * `log_sigmoid` and `sigmoid` fail some array-api-strict regression tests. It is known from SciPy that `log_expit` is not yet supported. * Update `atol` and `rtol` values for `float32` tests. This is a pain identifying `atol` and `rtol` values suitable for 1e4 tests across platforms, backends, devices, and dtypes. The vast majority of tests pass with stricter than `atol=1e-0` and `rtol=1e-2` for `float32`. For tests that wouldn't pass a stricter threshold, I have verified that the coeffs are qualitatively the same (often a big coefficient 1e2 that may be off by a little over 1e-1). You can also inspect the singular values (manually) to see that they are also qualitatively the same. --- src/gfdl/activations.py | 6 ++++-- src/gfdl/tests/test_regression.py | 14 ++------------ 2 files changed, 6 insertions(+), 14 deletions(-) diff --git a/src/gfdl/activations.py b/src/gfdl/activations.py index cd143e9..3e46d81 100644 --- a/src/gfdl/activations.py +++ b/src/gfdl/activations.py @@ -65,7 +65,8 @@ def sigmoid(z): -------- scipy.special.expit : The logistic sigmoid function. """ - return scipy.special.expit(z) + xp, _, _ = get_namespace_and_device(z) + return xp.exp(-xp.logaddexp(xp.zeros_like(z), -z)) def identity(z): @@ -147,7 +148,8 @@ def log_sigmoid(z): -------- scipy.special.log_expit : The logistic sigmoid function. """ - return scipy.special.log_expit(z) + xp, _, _ = get_namespace_and_device(z) + return -xp.logaddexp(xp.zeros_like(z), -z) def log_softmax(z): diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index d937aa4..743ab99 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -308,11 +308,6 @@ def test_xp_check_inputs(hidden_layer_sizes, device_name, dtype_name, ): - if (array_namespace not in ["numpy", "torch"] or - device_name not in [None, "cpu"] or - dtype_name not in [None, "float64"] - ): - pytest.skip("unsupported configuration") estimator = GFDLRegressor( reg_alpha=reg_alpha, @@ -355,15 +350,10 @@ def test_xp_make_regression(hidden_layer_sizes, # This test is closely based on `_check_array_api_core` from sklearn. # Here, we have better control of `atol`, `rtol`, and simulating data. # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. - if (array_namespace not in ["numpy", "torch"] or - device_name not in [None, "cpu"] or - dtype_name not in [None, "float64"] - ): - pytest.skip("unsupported configuration") xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) - atol = 1e-3 - rtol = 1e-6 + atol = 1e-3 if dtype_name == "float64" else 1e-0 + rtol = 1e-6 if dtype_name == "float64" else 1e-2 n_samples = 500 random_state = 42 estimator = GFDLRegressor From 5c590e60536a23dc46c3a546eb564e3b28e908a9 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 19 Oct 2026 18:19:18 -0600 Subject: [PATCH 18/32] WIP,ENH: `partial_fit` for array API * SI disclosure: I did not use SI * Implemented `xp` for base estimator's `partial_fit` method * Be wary of checking namespace without running `.fit()`, hence why you do not put this in `partial_fit` * First draft at adding `partial_fit` boolean to the custom test function * `xp.linalg.lstsq` does SVD as in `Ridge().fit()` but is not supported currently --- src/gfdl/model.py | 79 ++++++++++++++++++++++++++----- src/gfdl/tests/test_regression.py | 23 ++++++++- 2 files changed, 88 insertions(+), 14 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index fb96545..682d7fc 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -190,6 +190,15 @@ def partial_fit(self, X, y): # # It's evident that it is indeed possible to update normal equations # through summation + xp, _, device = get_namespace_and_device(X) + if not xp.isdtype(X.dtype, "real floating"): + default_float_dtype = _max_precision_float_dtype(xp, device) + X = xp.astype(X, default_float_dtype) + y = xp.astype( + move_to(y, xp=xp, device=device), + X.dtype, + copy=False, + ) if self.reg_alpha is not None and self.reg_alpha < 0.0: raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).") @@ -212,22 +221,54 @@ def partial_fit(self, X, y): rng = self.get_generator(self.seed) self.W_.append( + xp.astype( + move_to( self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) - ) + ), + xp=xp, + device=device, + ), + X.dtype, + copy=False, + ) ) self.b_.append( + xp.astype( + move_to( self._weight_mode(1, hidden_layer_sizes[0], rng=rng) - .reshape(-1) + .reshape(-1), + xp=xp, + device=device, + ), + X.dtype, + copy=False, + ) ) for i, layer in enumerate(hidden_layer_sizes[1:]): # (n_hidden, n_features) self.W_.append( - self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,) + xp.astype( + move_to( + self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), + xp=xp, + device=device, + ), + X.dtype, + copy=False, + ) ) # (n_hidden,) self.b_.append( - self._weight_mode(1, layer, rng=rng,).reshape(-1) + xp.astype( + move_to( + self._weight_mode(1, layer, rng=rng,).reshape(-1), + xp=xp, + device=device, + ), + X.dtype, + copy=False, + ) ) # design matrix shape: (n_samples, sum_hidden+n_features) @@ -243,11 +284,14 @@ def partial_fit(self, X, y): # or (n_samples, sum_hidden) if self.direct_links: Hs.append(X) - D = np.hstack(Hs) + D = xp.concat(Hs, axis=1) if not hasattr(self, "A"): - self.A = np.zeros((D.shape[1], D.shape[1])) - self.B = np.zeros((D.shape[1], y.shape[1])) + self.A = xp.zeros((D.shape[1], D.shape[1]), + dtype=X.dtype, device=device) + self.B = xp.zeros((D.shape[1], y.shape[1]), + dtype=X.dtype, device=device + ) self.A += D.T @ D self.B += D.T @ y @@ -271,7 +315,8 @@ def partial_fit(self, X, y): # = (D.T @ D)^-1 @ D.T @ y if self.reg_alpha is None: - self.coeff_ = np.linalg.pinv(self.A, rtol=self.rtol) @ self.B + coef = xp.linalg.pinv(self.A, rtol=self.rtol) @ self.B + self.coeff_ = xp.squeeze(coef, axis=-1) else: # scipy.linalg.solve(self.A + reg_mat, self.B) # is equivalent to @@ -279,10 +324,20 @@ def partial_fit(self, X, y): # alpha=self.reg_alpha, fit_intercept=False, solver='cholesky' # ) # ridge.fit(D_all[:current_batch], Y_all[:current_batch]) - reg_mat = np.eye(self.A.shape[0]) * self.reg_alpha - self.coeff_ = scipy.linalg.solve(self.A + reg_mat, self.B) - if self.coeff_.ndim == 2 and self.coeff_.shape[1] == 1: - self.coeff_ = self.coeff_.ravel() + reg_mat = xp.astype( + move_to( + xp.eye(self.A.shape[0]) * self.reg_alpha, + xp=xp, + device=device, + ), + X.dtype, + copy=False, + ) + + coef = xp.linalg.solve( + self.A + reg_mat, self.B + ) + self.coeff_ = xp.squeeze(coef, axis=-1) return self def predict(self, X): diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 743ab99..9b51948 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -333,6 +333,7 @@ def test_xp_check_inputs(hidden_layer_sizes, @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize("partial_fit", [False, True]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), @@ -345,6 +346,7 @@ def test_xp_make_regression(hidden_layer_sizes, array_namespace, device_name, dtype_name, + partial_fit, ): # Custom make_regression()-based NumPy and array API matching # This test is closely based on `_check_array_api_core` from sklearn. @@ -354,10 +356,21 @@ def test_xp_make_regression(hidden_layer_sizes, atol = 1e-3 if dtype_name == "float64" else 1e-0 rtol = 1e-6 if dtype_name == "float64" else 1e-2 + if partial_fit: + atol *= 1e1 + rtol *= 1e1 n_samples = 500 random_state = 42 estimator = GFDLRegressor + def run_partial_fit(X, y, model): + batch = 25 + for start in range(0, X.shape[0], batch): + end = min(start + batch, X.shape[0]) + Xb = X[start:end, :] + yb = y[start:end] + model.partial_fit(Xb, yb) + X, y = make_regression( n_samples=n_samples, n_features=20, @@ -377,12 +390,18 @@ def test_xp_make_regression(hidden_layer_sizes, direct_links=direct_links, seed=random_state,) with config_context(array_api_dispatch=False): - est.fit(X, y) + if partial_fit: + run_partial_fit(X, y, est) + else: + est.fit(X, y) # xp (array API) .fit est_xp = clone(est) with config_context(array_api_dispatch=True): - est_xp.fit(X_xp, y) + if partial_fit: + run_partial_fit(X_xp, y, est_xp) + else: + est_xp.fit(X_xp, y) # Check fitted attributes are close array_attributes = { From a93a3799fc1235c15fe76977c7c14637575d6e05 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Tue, 20 Oct 2026 11:35:05 -0600 Subject: [PATCH 19/32] WIP,ENH: softmax; tighter atol, rtol bounds for float64 * SI disclosure: I did not use SI * Replaced `scipy.special` versions of softmax functions with `xp` versions based on Emma's `sklearn` pull request * Updated `atol` and `rtol` bounds for tighter control in `xp_make_regression` test * Added quantitative comments about why we can't tighten further without failing the test (very few cases out of 720 tests over configurations) * We expect some classifier tests to fail, as the classifier will be addressed in #130 which is a descendant of this PR. The vast majority that fail right now concern `partial_fit` in the `GFDLClassifier` --- src/gfdl/activations.py | 12 +++++++++--- src/gfdl/tests/test_regression.py | 17 +++++++++++------ 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/gfdl/activations.py b/src/gfdl/activations.py index 3e46d81..e2b7b35 100644 --- a/src/gfdl/activations.py +++ b/src/gfdl/activations.py @@ -2,7 +2,6 @@ Activation functions for Gradient Free Deep Learning estimators. """ -import scipy from sklearn.utils._array_api import get_namespace_and_device @@ -104,7 +103,10 @@ def softmax(z): -------- scipy.special.softmax : The softmax function. """ - return scipy.special.softmax(z, axis=-1) + xp, _, _ = get_namespace_and_device(z) + tmp = z - xp.max(z, axis=1, keepdims=True) + exp_tmp = xp.exp(tmp) + return exp_tmp / xp.sum(exp_tmp, axis=1, keepdims=True) def softmin(z): @@ -170,7 +172,11 @@ def log_softmax(z): -------- scipy.special.log_softmax : The logarithm of softmax function. """ - return scipy.special.log_softmax(z, axis=-1) + xp, _, _ = get_namespace_and_device(z) + shifted = z - xp.max(z, axis=1, keepdims=True) + return shifted - xp.log( + xp.sum(xp.exp(shifted), axis=1, keepdims=True) + ) ACTIVATIONS = { diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 9b51948..ac48582 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -328,7 +328,7 @@ def test_xp_check_inputs(hidden_layer_sizes, @skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @@ -354,17 +354,22 @@ def test_xp_make_regression(hidden_layer_sizes, # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) - atol = 1e-3 if dtype_name == "float64" else 1e-0 - rtol = 1e-6 if dtype_name == "float64" else 1e-2 + # torch float64 (MacOS) .fit() : 13/720 fail with atol=1e-11 & rtol=1e-14 + # torch float32 (MacOS) .fit() : 5/720 fail with atol=1e-2 & rtol=1e-4 + # torch float64 (MacOS) .partial_fit() : 23/720 fail with atol=1e-7 & rtol=1e-9 + # torch float32 (MacOS) .partial_fit() : 10/720 fail with atol=1e1 & rtol=1e-3 + # Increased by order of magnitude from notes above for Linux + atol = 1e-9 if dtype_name == "float64" else 1e-1 + rtol = 1e-12 if dtype_name == "float64" else 1e-3 if partial_fit: - atol *= 1e1 - rtol *= 1e1 + atol = 1e-6 if dtype_name == "float64" else 1e-0 + rtol = 1e-8 if dtype_name == "float64" else 1e-1 n_samples = 500 random_state = 42 estimator = GFDLRegressor def run_partial_fit(X, y, model): - batch = 25 + batch = 100 for start in range(0, X.shape[0], batch): end = min(start + batch, X.shape[0]) Xb = X[start:end, :] From d53840ffdfa2988a9ebdd677cf6703b52c8fc869 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Tue, 20 Oct 2026 12:03:16 -0600 Subject: [PATCH 20/32] WIP,TST: `sklearn` api compliance already handles array API * Removed the `_NotAnArray` try/except block. The `sklearn` compliance function already handles which tests to run if SCIPY_ARRAY_API=1 * Removed the related test --- src/gfdl/tests/test_regression.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index ac48582..7fd2e70 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -274,21 +274,6 @@ def test_preserve_class_inputs(): assert isinstance(v, type(expected[k])) -@skip_if_array_api_compat_not_configured -@parametrize_with_checks([GFDLRegressor()]) -def test_sklearn_array_api_conformance(estimator, check,): - with config_context(array_api_dispatch=True): - try: - check(estimator) - except Exception as exc: - if "_NotAnArray" in repr(exc): - pytest.xfail( - "Full check_estimator under array_api_dispatch=True is expected " - "to fail on sklearn's _NotAnArray sentinel." - ) - raise - - @skip_if_array_api_compat_not_configured @pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) From ee40304fd788ff879d07c15850578ac974958010 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Fri, 23 Oct 2026 15:14:36 -0600 Subject: [PATCH 21/32] WIP,ENH,TST: revisions after Emma's PR * SI disclosure: I did not use SI, except I did use GPT 6.1 Sol when reviewing Emma's work * There were some Array API styles that Emma used in her sklearn PR that I preferred over mine. I also identified some redundant tests in my sklearn PR review. Here I am porting over those suggested changes. * https://github.com/eiviani-lanl/scikit-learn/pull/2 * Use `xp.asarray` instead of `move_to` and `xp.astype` * Use `dtype=[xp.float64, xp.float32]` in `validate_data` * Only `xp.squeeze` when needed. Fixes a bug with classifier which will go upstream when I rebase #130 * Removed redundant tests which were identified during review --- src/gfdl/model.py | 98 +++++++++++++++---------------- src/gfdl/tests/test_regression.py | 64 +------------------- 2 files changed, 51 insertions(+), 111 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 682d7fc..567c1ab 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -16,7 +16,6 @@ from sklearn.preprocessing import OneHotEncoder from sklearn.utils import column_or_1d from sklearn.utils._array_api import ( - _max_precision_float_dtype, check_same_namespace, get_namespace_and_device, move_to, @@ -57,9 +56,6 @@ def __sklearn_tags__(self): def fit(self, X, Y): xp, _, device = get_namespace_and_device(X) - if not xp.isdtype(X.dtype, "real floating"): - default_float_dtype = _max_precision_float_dtype(xp, device) - X = xp.astype(X, default_float_dtype) Y = xp.astype( move_to(Y, xp=xp, device=device), X.dtype, @@ -87,53 +83,37 @@ def fit(self, X, Y): rng = self.get_generator(self.seed) self.W_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) ), - xp=xp, + dtype=X.dtype, device=device, - ), - X.dtype, - copy=False, - ) ) + ) self.b_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode(1, hidden_layer_sizes[0], rng=rng) .reshape(-1), - xp=xp, + dtype=X.dtype, device=device, - ), - X.dtype, - copy=False, ) ) for i, layer in enumerate(hidden_layer_sizes[1:]): # (n_hidden, n_features) self.W_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), - xp=xp, + dtype=X.dtype, device=device, - ), - X.dtype, - copy=False, ) ) # (n_hidden,) self.b_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode(1, layer, rng=rng,).reshape(-1), - xp=xp, + dtype=X.dtype, device=device, - ), - X.dtype, - copy=False, ) ) @@ -191,9 +171,6 @@ def partial_fit(self, X, y): # It's evident that it is indeed possible to update normal equations # through summation xp, _, device = get_namespace_and_device(X) - if not xp.isdtype(X.dtype, "real floating"): - default_float_dtype = _max_precision_float_dtype(xp, device) - X = xp.astype(X, default_float_dtype) y = xp.astype( move_to(y, xp=xp, device=device), X.dtype, @@ -316,7 +293,9 @@ def partial_fit(self, X, y): if self.reg_alpha is None: coef = xp.linalg.pinv(self.A, rtol=self.rtol) @ self.B - self.coeff_ = xp.squeeze(coef, axis=-1) + if coef.ndim == 2 and coef.shape[1] == 1: + coef = xp.squeeze(coef, axis=-1) + self.coeff_ = coef else: # scipy.linalg.solve(self.A + reg_mat, self.B) # is equivalent to @@ -337,16 +316,15 @@ def partial_fit(self, X, y): coef = xp.linalg.solve( self.A + reg_mat, self.B ) - self.coeff_ = xp.squeeze(coef, axis=-1) + if coef.ndim == 2 and coef.shape[1] == 1: + coef = xp.squeeze(coef, axis=-1) + self.coeff_ = coef return self def predict(self, X): check_is_fitted(self) - xp, _, device = get_namespace_and_device(X) + xp, _, _ = get_namespace_and_device(X) check_same_namespace(X, self, attribute="coeff_", method="predict") - if not xp.isdtype(X.dtype, "real floating"): - default_float_dtype = _max_precision_float_dtype(xp, device) - X = xp.astype(X, default_float_dtype) Hs = [] H_prev = X @@ -530,7 +508,9 @@ def fit(self, X, y): Fitted estimator. """ # shape: (n_samples, n_features) - X, Y = validate_data(self, X, y) + xp, _, _ = get_namespace_and_device(X) + X, Y = validate_data(self, X, y, + dtype=[xp.float64, xp.float32]) self.classes_ = unique_labels(Y) # onehot y @@ -594,7 +574,9 @@ def partial_fit(self, X, y, classes=None): numerical stability and avoid loss of information. """ # shape: (n_samples, n_features) - X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_")) + xp, _, _ = get_namespace_and_device(X) + X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"), + dtype=[xp.float64, xp.float32]) Y = column_or_1d(Y, warn=True) check_classification_targets(Y) @@ -630,7 +612,9 @@ def predict(self, X): The predicted classes, with shape (n_samples,) or (n_samples, n_outputs). """ check_is_fitted(self) - X = validate_data(self, X, reset=False) + xp, _, _ = get_namespace_and_device(X) + X = validate_data(self, X, reset=False, + dtype=[xp.float64, xp.float32]) out = self.predict_proba(X) y_hat = self.classes_[np.argmax(out, axis=1)] return y_hat @@ -652,7 +636,9 @@ def predict_proba(self, X): have shape (n_samples, n_classes). """ check_is_fitted(self) - X = validate_data(self, X, reset=False) + xp, _, _ = get_namespace_and_device(X) + X = validate_data(self, X, reset=False, + dtype=[xp.float64, xp.float32]) out = super().predict(X) out = np.exp(out - logsumexp(out, axis=1, keepdims=True)) return out @@ -995,7 +981,9 @@ def fit(self, X, y): The fitted estimator. """ # shape: (n_samples, n_features) - X, Y = validate_data(self, X, y) + xp, _, _ = get_namespace_and_device(X) + X, Y = validate_data(self, X, y, + dtype=[xp.float64, xp.float32]) self.classes_ = unique_labels(Y) # onehot y @@ -1059,7 +1047,9 @@ def partial_fit(self, X, y, classes=None): numerical stability and avoid loss of information. """ # shape: (n_samples, n_features) - X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_")) + xp, _, _ = get_namespace_and_device(X) + X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"), + dtype=[xp.float64, xp.float32]) Y = column_or_1d(Y, warn=True) check_classification_targets(Y) @@ -1107,7 +1097,9 @@ def predict_proba(self, X): have shape (n_samples, n_classes). """ check_is_fitted(self) - X = validate_data(self, X, reset=False) + xp, _, _ = get_namespace_and_device(X) + X = validate_data(self, X, reset=False, + dtype=[xp.float64, xp.float32]) outs = self._forward(X) probs = [] @@ -1133,7 +1125,9 @@ def predict(self, X): The predicted classes, with shape (n_samples,) or (n_samples, n_outputs). """ check_is_fitted(self) - X = validate_data(self, X, reset=False) + xp, _, _ = get_namespace_and_device(X) + X = validate_data(self, X, reset=False, + dtype=[xp.float64, xp.float32]) if self.voting == "soft": P = self.predict_proba(X) @@ -1312,7 +1306,9 @@ def fit(self, X, y): object The fitted estimator. """ - X, Y = validate_data(self, X, y, multi_output=True) + xp, _, _ = get_namespace_and_device(X) + X, Y = validate_data(self, X, y, multi_output=True, + dtype=[xp.float64, xp.float32]) super().fit(X, Y) return self @@ -1360,7 +1356,9 @@ def partial_fit(self, X, y): the condition number is squared. This may require a lower ``rtol`` to maintain numerical stability and avoid loss of information. """ - X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_")) + xp, _, _ = get_namespace_and_device(X) + X, Y = validate_data(self, X, y, reset=not hasattr(self, "n_features_in_"), + dtype=[xp.float64, xp.float32]) if Y.ndim == 1: Y = Y.reshape(-1, 1) super().partial_fit(X, Y) @@ -1383,5 +1381,7 @@ def predict(self, X): (n_samples, n_outputs). """ check_is_fitted(self) - X = validate_data(self, X, reset=False) + xp, _, _ = get_namespace_and_device(X) + X = validate_data(self, X, reset=False, + dtype=[xp.float64, xp.float32]) return super().predict(X) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 7fd2e70..9d244ca 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -18,7 +18,6 @@ skip_if_array_api_compat_not_configured, ) from sklearn.utils.estimator_checks import ( - check_array_api_input, parametrize_with_checks, ) @@ -274,51 +273,12 @@ def test_preserve_class_inputs(): assert isinstance(v, type(expected[k])) -@skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) -@pytest.mark.parametrize("direct_links", [True, False]) -@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) -@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) -@pytest.mark.parametrize("reg_alpha", [None, 0.1]) -@pytest.mark.parametrize( - "array_namespace, device_name, dtype_name", - yield_namespace_device_dtype_combinations(), -) -def test_xp_check_inputs(hidden_layer_sizes, - direct_links, - activation, - weight_scheme, - reg_alpha, - array_namespace, - device_name, - dtype_name, - ): - - estimator = GFDLRegressor( - reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - activation=activation, - weight_scheme=weight_scheme, - direct_links=direct_links, - seed=42, - ) - with config_context(array_api_dispatch=True): - check_array_api_input( - estimator.__class__.__name__, - estimator, - array_namespace=array_namespace, - device_name=device_name, - dtype_name=dtype_name, - ) - - @skip_if_array_api_compat_not_configured @pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @pytest.mark.parametrize("reg_alpha", [None, 0.1]) -@pytest.mark.parametrize("partial_fit", [False, True]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), @@ -331,7 +291,6 @@ def test_xp_make_regression(hidden_layer_sizes, array_namespace, device_name, dtype_name, - partial_fit, ): # Custom make_regression()-based NumPy and array API matching # This test is closely based on `_check_array_api_core` from sklearn. @@ -341,26 +300,13 @@ def test_xp_make_regression(hidden_layer_sizes, # torch float64 (MacOS) .fit() : 13/720 fail with atol=1e-11 & rtol=1e-14 # torch float32 (MacOS) .fit() : 5/720 fail with atol=1e-2 & rtol=1e-4 - # torch float64 (MacOS) .partial_fit() : 23/720 fail with atol=1e-7 & rtol=1e-9 - # torch float32 (MacOS) .partial_fit() : 10/720 fail with atol=1e1 & rtol=1e-3 # Increased by order of magnitude from notes above for Linux atol = 1e-9 if dtype_name == "float64" else 1e-1 rtol = 1e-12 if dtype_name == "float64" else 1e-3 - if partial_fit: - atol = 1e-6 if dtype_name == "float64" else 1e-0 - rtol = 1e-8 if dtype_name == "float64" else 1e-1 n_samples = 500 random_state = 42 estimator = GFDLRegressor - def run_partial_fit(X, y, model): - batch = 100 - for start in range(0, X.shape[0], batch): - end = min(start + batch, X.shape[0]) - Xb = X[start:end, :] - yb = y[start:end] - model.partial_fit(Xb, yb) - X, y = make_regression( n_samples=n_samples, n_features=20, @@ -380,18 +326,12 @@ def run_partial_fit(X, y, model): direct_links=direct_links, seed=random_state,) with config_context(array_api_dispatch=False): - if partial_fit: - run_partial_fit(X, y, est) - else: - est.fit(X, y) + est.fit(X, y) # xp (array API) .fit est_xp = clone(est) with config_context(array_api_dispatch=True): - if partial_fit: - run_partial_fit(X_xp, y, est_xp) - else: - est_xp.fit(X_xp, y) + est_xp.fit(X_xp, y) # Check fitted attributes are close array_attributes = { From 01c1228675f18aaf1fef20b4bb66a6da04f66af6 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 26 Oct 2026 14:41:50 -0600 Subject: [PATCH 22/32] WIP,TST: `xpx` assert, small `partial_fit` test * SI disclosure: I used SI to write a simple `partial_fit` test function * Replaced `numpy.testing.assert_allclose` with `xpx.testing.assert_close` for array API tests * Decreased number of Array API test cases --- src/gfdl/tests/test_regression.py | 102 ++++++++++++++++++++++-------- 1 file changed, 74 insertions(+), 28 deletions(-) diff --git a/src/gfdl/tests/test_regression.py b/src/gfdl/tests/test_regression.py index 9d244ca..a0b6c89 100644 --- a/src/gfdl/tests/test_regression.py +++ b/src/gfdl/tests/test_regression.py @@ -2,14 +2,15 @@ import numpy as np import pytest -from numpy.testing import assert_allclose from sklearn import config_context from sklearn.base import clone from sklearn.datasets import fetch_openml, make_regression +from sklearn.externals import array_api_extra as xpx from sklearn.metrics import r2_score from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.utils._array_api import ( + _atol_for_type, move_to, yield_namespace_device_dtype_combinations, ) @@ -274,16 +275,25 @@ def test_preserve_class_inputs(): @skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) -@pytest.mark.parametrize("direct_links", [True, False]) -@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) -@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) -@pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize( + "reg_alpha, hidden_layer_sizes, activation, weight_scheme, direct_links", + [ + (2, (10,), "relu", "normal", True), + (2, (10, 10), "sigmoid", "zeros", True), + (2, (10,), "softmax", "lecun_normal", True), + (2, (10,), "tanh", "he_normal", False), + (None, (10,), "softmin", "glorot_normal", True), + (None, (10, 10), "log_softmax", "uniform", False), + (None, (10,), "log_sigmoid", "lecun_uniform", True), + (2, (10,), "identity", "glorot_uniform", True), + (2, (10,), "identity", "he_uniform", True), + ] +) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), ) -def test_xp_make_regression(hidden_layer_sizes, +def test_xp_regressor_fit(hidden_layer_sizes, direct_links, activation, weight_scheme, @@ -298,21 +308,11 @@ def test_xp_make_regression(hidden_layer_sizes, # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) - # torch float64 (MacOS) .fit() : 13/720 fail with atol=1e-11 & rtol=1e-14 - # torch float32 (MacOS) .fit() : 5/720 fail with atol=1e-2 & rtol=1e-4 - # Increased by order of magnitude from notes above for Linux - atol = 1e-9 if dtype_name == "float64" else 1e-1 - rtol = 1e-12 if dtype_name == "float64" else 1e-3 - n_samples = 500 - random_state = 42 + atol = _atol_for_type(dtype_name) * 10 estimator = GFDLRegressor X, y = make_regression( - n_samples=n_samples, - n_features=20, - n_informative=5, - n_targets=1, - random_state=random_state, + n_samples=500, n_features=20, n_informative=5, n_targets=1, random_state=0, ) X = X.astype(dtype_name, copy=False) X_xp = xp.asarray(X, device=device) @@ -324,7 +324,7 @@ def test_xp_make_regression(hidden_layer_sizes, activation=activation, weight_scheme=weight_scheme, direct_links=direct_links, - seed=random_state,) + seed=0,) with config_context(array_api_dispatch=False): est.fit(X, y) @@ -340,12 +340,10 @@ def test_xp_make_regression(hidden_layer_sizes, for attribute_name, attribute_value in array_attributes.items(): est_xp_attr = getattr(est_xp, attribute_name) est_xp_attr_np = move_to(est_xp_attr, xp=np, device="cpu") - assert_allclose( - attribute_value, - est_xp_attr_np, + xpx.testing.assert_close( + attribute_value, est_xp_attr_np, err_msg=f"{attribute_name} not the same", atol=atol, - rtol=rtol, ) # Check method outputs are close @@ -358,10 +356,58 @@ def test_xp_make_regression(hidden_layer_sizes, with config_context(array_api_dispatch=True): result_xp = getattr(est_xp, method_name)(X_xp) result_xp_np = move_to(result_xp, xp=np, device="cpu") - assert_allclose( - result, - result_xp_np, + xpx.testing.assert_close( + result, result_xp_np, err_msg=f"{method_name}() did not the return the same result", atol=atol, - rtol=rtol, ) + + +@skip_if_array_api_compat_not_configured +@pytest.mark.parametrize("reg_alpha", [None, 2]) +@pytest.mark.parametrize("hidden_layer_sizes", [(10,), (10, 10)]) +@pytest.mark.parametrize( + "array_namespace, device_name, dtype_name", + yield_namespace_device_dtype_combinations(), +) +def test_xp_regressor_partial_fit(hidden_layer_sizes, + reg_alpha, + array_namespace, + device_name, + dtype_name, + ): + + atol = _atol_for_type(dtype_name) * 10 + batch = 2 + + xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) + X_np, y_np = make_regression( + n_samples=6, n_features=20, n_informative=5, random_state=0, + ) + X_np = X_np.astype(dtype_name, copy=False) + X_xp = xp.asarray(X_np, device=device) + + params = { + "reg_alpha": reg_alpha, + "hidden_layer_sizes": hidden_layer_sizes, + "seed": 0 + } + estimator = GFDLRegressor + + predictions = [] + for X, dispatch in [(X_np, False), (X_xp, True)]: + with config_context(array_api_dispatch=dispatch): + model = estimator(**params) + for start in range(0, X.shape[0], batch): + model.partial_fit( + X[start:start + batch, ...], y_np[start:start + batch], + ) + predictions.append( + move_to(model.predict(X), xp=np, device="cpu") + ) + + xpx.testing.assert_close( + predictions[0], predictions[1], + err_msg="predict() did not the return the same result", + atol=atol, + ) From 1e84cdf600dc6ac26db32dbe86133b183e60b29b Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 26 Oct 2026 14:44:46 -0600 Subject: [PATCH 23/32] WIP,ENH: `xp.asarray` in `partial_fit` * SI disclosure: I did not use SI in this commit * Copy and pasted the `xp.asarray` set up for the `self.W_` and `self.b_` initialization to `partial_fit` * Simplified some lines to decrease line number --- src/gfdl/model.py | 76 +++++++++++------------------------------------ 1 file changed, 18 insertions(+), 58 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 567c1ab..2e88977 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -58,8 +58,7 @@ def fit(self, X, Y): xp, _, device = get_namespace_and_device(X) Y = xp.astype( move_to(Y, xp=xp, device=device), - X.dtype, - copy=False, + X.dtype, copy=False ) # Assumption : X, Y have been pre-processed. @@ -87,16 +86,13 @@ def fit(self, X, Y): self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) ), - dtype=X.dtype, - device=device, + dtype=X.dtype, device=device ) ) self.b_.append( xp.asarray( self._weight_mode(1, hidden_layer_sizes[0], rng=rng) - .reshape(-1), - dtype=X.dtype, - device=device, + .reshape(-1), dtype=X.dtype, device=device ) ) for i, layer in enumerate(hidden_layer_sizes[1:]): @@ -104,16 +100,14 @@ def fit(self, X, Y): self.W_.append( xp.asarray( self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), - dtype=X.dtype, - device=device, + dtype=X.dtype, device=device ) ) # (n_hidden,) self.b_.append( xp.asarray( self._weight_mode(1, layer, rng=rng,).reshape(-1), - dtype=X.dtype, - device=device, + dtype=X.dtype, device=device ) ) @@ -198,54 +192,28 @@ def partial_fit(self, X, y): rng = self.get_generator(self.seed) self.W_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode( self._N, hidden_layer_sizes[0], rng=self.get_generator(self.seed) - ), - xp=xp, - device=device, - ), - X.dtype, - copy=False, - ) + ), dtype=X.dtype, device=device) ) self.b_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode(1, hidden_layer_sizes[0], rng=rng) - .reshape(-1), - xp=xp, - device=device, - ), - X.dtype, - copy=False, - ) + .reshape(-1), dtype=X.dtype, device=device) ) for i, layer in enumerate(hidden_layer_sizes[1:]): # (n_hidden, n_features) self.W_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode(hidden_layer_sizes[i], layer, rng=rng,), - xp=xp, - device=device, - ), - X.dtype, - copy=False, - ) + dtype=X.dtype, device=device) ) # (n_hidden,) self.b_.append( - xp.astype( - move_to( + xp.asarray( self._weight_mode(1, layer, rng=rng,).reshape(-1), - xp=xp, - device=device, - ), - X.dtype, - copy=False, - ) + dtype=X.dtype, device=device) ) # design matrix shape: (n_samples, sum_hidden+n_features) @@ -267,8 +235,7 @@ def partial_fit(self, X, y): self.A = xp.zeros((D.shape[1], D.shape[1]), dtype=X.dtype, device=device) self.B = xp.zeros((D.shape[1], y.shape[1]), - dtype=X.dtype, device=device - ) + dtype=X.dtype, device=device) self.A += D.T @ D self.B += D.T @ y @@ -303,19 +270,12 @@ def partial_fit(self, X, y): # alpha=self.reg_alpha, fit_intercept=False, solver='cholesky' # ) # ridge.fit(D_all[:current_batch], Y_all[:current_batch]) - reg_mat = xp.astype( - move_to( - xp.eye(self.A.shape[0]) * self.reg_alpha, - xp=xp, - device=device, - ), - X.dtype, - copy=False, + reg_mat = ( + xp.eye(self.A.shape[0], dtype=self.A.dtype, device=device) + * self.reg_alpha ) - coef = xp.linalg.solve( - self.A + reg_mat, self.B - ) + coef = xp.linalg.solve(self.A + reg_mat, self.B) if coef.ndim == 2 and coef.shape[1] == 1: coef = xp.squeeze(coef, axis=-1) self.coeff_ = coef From d1fa627bc8474db88deb26e7254cacf6b9ced35d Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 14 Sep 2026 12:37:17 -0600 Subject: [PATCH 24/32] WIP, ENH: torch CPU for GFDLClassifier - private utility function for OneHotEncoder - one line call to it in .fit() - GFDLRegressor only has to invoked base GFDL - GFDLClassifier has that 1 line change in .fit() - Test function for numpy vs torch matching - Copy, paste, modify from test_regression.py - Other test functions in regressor PR would be redundant --- src/gfdl/model.py | 23 ++++++++++++ src/gfdl/tests/test_model.py | 72 ++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 2e88977..81014de 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -1345,3 +1345,26 @@ def predict(self, X): X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) return super().predict(X) + + +def _to_numpy_cpu_y(y): + """Convert label array y to a 1D NumPy array on CPU.""" + # if y is None: + # raise ValueError("y cannot be None for a supervised estimator.") + + # PyTorch: handles CPU and CUDA tensors. + if hasattr(y, "detach") and hasattr(y, "cpu") and hasattr(y, "numpy"): + y = y.detach().cpu().numpy() + + # CuPy: GPU -> CPU NumPy. + elif hasattr(y, "get"): + y = y.get() + + # Generic fallback: NumPy, JAX CPU/device arrays, pandas, lists, etc. + else: + y = np.asarray(y) + + # Normalize shape for sklearn classifier utilities. + y = column_or_1d(y, warn=True) + + return y diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 3225d41..54be0b5 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -1,6 +1,10 @@ +import os + import numpy as np import pytest +import torch from numpy.testing import assert_allclose +from sklearn import config_context from sklearn.base import clone from sklearn.datasets import load_breast_cancer, load_digits, make_classification from sklearn.metrics import accuracy_score, roc_auc_score @@ -824,3 +828,71 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): clf = classifier(hidden_layer_sizes=hidden_layer_sizes, seed=0) with pytest.raises(ValueError, match="must be > 0"): clf.fit(X, y) + + +@pytest.mark.parametrize( + "reg_alpha, n_features, hidden_layer_sizes, n_classes", + [ + (1e-1, 40, (100,), 2), + (1e-1, 400, (100,), 2), + (2, 40, (100,), 2), + (2, 400, (100,), 2), + (None, 40, (100,), 2), + (None, 400, (100,), 2), + (1e-1, 40, (100, 100,), 3), + (1e-1, 400, (100, 100,), 3), + (2, 40, (100, 100,), 3), + (2, 400, (100, 100,), 3), + (None, 40, (100, 100,), 3), + (None, 400, (100, 100,), 3), + ] +) +def test_torch_matches_numpy(reg_alpha, + n_features, + hidden_layer_sizes, + n_classes, + ): + """NumPy and Torch predictions should be close up to atol""" + os.environ["SCIPY_ARRAY_API"] = "1" + estimator = GFDLClassifier + report = accuracy_score + rng = np.random.default_rng(seed=42) + acc_np_s = [] + acc_torch_s = [] + for _ in range(10): + random_state = rng.integers(low=0, high=100, size=1)[0] + X_np, y_np = make_classification( + n_samples=10_000, + n_features=n_features, + n_informative=int(n_features / 10), + n_classes=n_classes, + random_state=random_state, + ) + X_torch = torch.asarray(X_np, device="cpu",) + y_torch = torch.asarray(y_np, device="cpu",) + with config_context(array_api_dispatch=True): + model = estimator(reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + seed=random_state, + ) + model.fit(X_np, y_np) + y_pred = model.predict(X_torch) + acc_torch = report( + y_torch, + y_pred, + ) + acc_torch_s.append(acc_torch) + # test against numpy which is current usage + with config_context(array_api_dispatch=True): + model = estimator(reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + seed=random_state, + ) + model.fit(X_np, y_np) + y_pred = model.predict(X_np) + acc_np = report( + y_np, + y_pred, + ) + acc_np_s.append(acc_np) + assert_allclose(acc_torch_s, acc_np_s, atol=1e-3) From 29e15e08206466f0ff9672309f1d1d14fc652d64 Mon Sep 17 00:00:00 2001 From: Seth Temple Date: Mon, 28 Sep 2026 21:13:06 -0600 Subject: [PATCH 25/32] WIP,ENH: fixes after rebase from regressor * Order of `unique_labels` mattered for conformance * Test with larger dataset for `make_classification`. Tyler doesn't like this test but for the commit it confirms similar results for torch versus numpy for ALL weight schemes * Change .predict_proba and .predict in classifier to xp array versions * scipy.special.logsumexp is supported for array API * 6 tests fail with array-api-strict versus 0 for compat. Issue appears at `_to_numpy_cpu_y` * We want to remove that private utility eventually --- src/gfdl/model.py | 7 +- src/gfdl/tests/test_model.py | 187 ++++++++++++++++++++++++----------- 2 files changed, 132 insertions(+), 62 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 81014de..81e71db 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -472,6 +472,7 @@ def fit(self, X, y): X, Y = validate_data(self, X, y, dtype=[xp.float64, xp.float32]) self.classes_ = unique_labels(Y) + Y = _to_numpy_cpu_y(Y) # onehot y # (this is necessary for everything beyond binary classification) @@ -576,7 +577,8 @@ def predict(self, X): X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) out = self.predict_proba(X) - y_hat = self.classes_[np.argmax(out, axis=1)] + xp, _, _device = get_namespace_and_device(self.coeff_) + y_hat = self.classes_[xp.argmax(out, axis=1)] return y_hat def predict_proba(self, X): @@ -600,7 +602,8 @@ def predict_proba(self, X): X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) out = super().predict(X) - out = np.exp(out - logsumexp(out, axis=1, keepdims=True)) + xp, _, _device = get_namespace_and_device(self.coeff_) + out = xp.exp(out - scipy.special.logsumexp(out, axis=1, keepdims=True)) return out diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 54be0b5..39a6995 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -2,7 +2,6 @@ import numpy as np import pytest -import torch from numpy.testing import assert_allclose from sklearn import config_context from sklearn.base import clone @@ -10,7 +9,13 @@ from sklearn.metrics import accuracy_score, roc_auc_score from sklearn.model_selection import StratifiedKFold, train_test_split from sklearn.preprocessing import OneHotEncoder, StandardScaler -from sklearn.utils.estimator_checks import parametrize_with_checks +from sklearn.utils._testing import ( + skip_if_array_api_compat_not_configured, +) +from sklearn.utils.estimator_checks import ( + check_array_api_input_and_values, + parametrize_with_checks, +) from ucimlrepo import fetch_ucirepo from gfdl.activations import ACTIVATIONS @@ -334,10 +339,10 @@ def test_soft_and_hard_can_differ(alpha): def test_invalid_activation_weight(Classifier): X = np.zeros((30, 4)) y = np.zeros((30,)) - invalid_act = Classifier(hidden_layer_sizes=100, + invalid_act = Classifier(hidden_layer_sizes=(100,), activation="bogus_activation", weight_scheme="uniform") - invalid_weight = Classifier(hidden_layer_sizes=100, + invalid_weight = Classifier(hidden_layer_sizes=(100,), activation="identity", weight_scheme="bogus_weight") # the sklearn estimator API bans input validation in __init__, @@ -830,69 +835,131 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): clf.fit(X, y) -@pytest.mark.parametrize( - "reg_alpha, n_features, hidden_layer_sizes, n_classes", - [ - (1e-1, 40, (100,), 2), - (1e-1, 400, (100,), 2), - (2, 40, (100,), 2), - (2, 400, (100,), 2), - (None, 40, (100,), 2), - (None, 400, (100,), 2), - (1e-1, 40, (100, 100,), 3), - (1e-1, 400, (100, 100,), 3), - (2, 40, (100, 100,), 3), - (2, 400, (100, 100,), 3), - (None, 40, (100, 100,), 3), - (None, 400, (100, 100,), 3), - ] +@skip_if_array_api_compat_not_configured +@pytest.mark.skipif( + os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." +) +@parametrize_with_checks([GFDLClassifier()]) +def test_sklearn_array_api_conformance(estimator, check,): + with config_context(array_api_dispatch=True): + try: + check(estimator) + except Exception as exc: + if "NotAnArray" in repr(exc): + pytest.xfail( + "Full check_estimator under array_api_dispatch=True is expected " + "to fail on sklearn's NotAnArray sentinel." + ) + raise + + +@skip_if_array_api_compat_not_configured +@pytest.mark.skipif( + os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." +) +@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("direct_links", [True, False]) +@pytest.mark.parametrize("activation", ["identity", + "relu", + "tanh", + ]) +@pytest.mark.parametrize("weight_scheme", ["normal", + "zeros", + ]) +@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) +def test_torch_array_api(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + ): + """NumPy and Array API public attributes and predictions should be close""" + + pytest.importorskip("torch") + + estimator = GFDLClassifier( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=42, + ) + with config_context(array_api_dispatch=True): + # Simulates from X, y = make_classification(n_samples=30, n_features=10) + # Unstable across backends for some weight schemes + check_array_api_input_and_values( + estimator.__class__.__name__, + estimator, + array_namespace="torch", + device_name="cpu", + dtype_name="float64", + ) + +@skip_if_array_api_compat_not_configured +@pytest.mark.skipif( + os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." ) -def test_torch_matches_numpy(reg_alpha, - n_features, - hidden_layer_sizes, - n_classes, - ): - """NumPy and Torch predictions should be close up to atol""" - os.environ["SCIPY_ARRAY_API"] = "1" +@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("direct_links", [True, False]) +@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) +@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) +@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) +@pytest.mark.parametrize("n_classes", [2, 3]) +def test_torch_make_classificatino(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + n_classes, + ): + """Custom make_classification()-based NumPy and PyTorch matching""" + pytest.importorskip("torch") + import torch estimator = GFDLClassifier - report = accuracy_score rng = np.random.default_rng(seed=42) - acc_np_s = [] - acc_torch_s = [] - for _ in range(10): + + pred_torch = [] + pred_numpy = [] + coef_torch = [] + coef_numpy = [] + for _ in range(3): random_state = rng.integers(low=0, high=100, size=1)[0] X_np, y_np = make_classification( - n_samples=10_000, - n_features=n_features, - n_informative=int(n_features / 10), + n_samples=250, + n_features=20, + n_informative=5, n_classes=n_classes, random_state=random_state, ) + + # NumPy .fit and .predict + model = estimator( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=random_state,) + model.fit(X_np, y_np) + y_pred_np = model.predict(X_np) + pred_numpy.append(y_pred_np) + coef_numpy.append(model.coeff_) + + # PyTorch .fit and .predict X_torch = torch.asarray(X_np, device="cpu",) y_torch = torch.asarray(y_np, device="cpu",) - with config_context(array_api_dispatch=True): - model = estimator(reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - seed=random_state, - ) - model.fit(X_np, y_np) - y_pred = model.predict(X_torch) - acc_torch = report( - y_torch, - y_pred, - ) - acc_torch_s.append(acc_torch) - # test against numpy which is current usage - with config_context(array_api_dispatch=True): - model = estimator(reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - seed=random_state, - ) - model.fit(X_np, y_np) - y_pred = model.predict(X_np) - acc_np = report( - y_np, - y_pred, - ) - acc_np_s.append(acc_np) - assert_allclose(acc_torch_s, acc_np_s, atol=1e-3) + model = estimator( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=random_state,) + model.fit(X_torch, y_torch) + y_pred_torch = model.predict(X_torch) + pred_torch.append(y_pred_torch) + coef_torch.append(model.coeff_) + + assert_allclose(pred_numpy, pred_torch) + assert_allclose(coef_numpy, coef_torch) \ No newline at end of file From 4e24e9cb6c022a2eb0e31870120defa7455358f4 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Wed, 14 Oct 2026 19:27:17 -0600 Subject: [PATCH 26/32] WIP,ENH,TST: revised tests, move_to classes * AI disclosure: I consulted AI but wrote the edits myself * Copy, paste, and modify the tests from `test_regression.py` in #129 * Use `move_to` on the `self.classes_` * `check_same_namespace` called individually in each estimator's methods --- src/gfdl/model.py | 35 ++---- src/gfdl/tests/test_model.py | 207 ++++++++++++++++++++--------------- 2 files changed, 127 insertions(+), 115 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 81e71db..eb68e45 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -471,8 +471,13 @@ def fit(self, X, y): xp, _, _ = get_namespace_and_device(X) X, Y = validate_data(self, X, y, dtype=[xp.float64, xp.float32]) - self.classes_ = unique_labels(Y) - Y = _to_numpy_cpu_y(Y) + xp, _, device = get_namespace_and_device(X) + self.classes_ = move_to( + unique_labels(Y), + xp=xp, + device=device, + ) + Y = move_to(Y, xp=np, device="cpu") # onehot y # (this is necessary for everything beyond binary classification) @@ -573,6 +578,7 @@ def predict(self, X): The predicted classes, with shape (n_samples,) or (n_samples, n_outputs). """ check_is_fitted(self) + check_same_namespace(X, self, attribute="coeff_", method="predict") xp, _, _ = get_namespace_and_device(X) X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) @@ -598,6 +604,7 @@ def predict_proba(self, X): have shape (n_samples, n_classes). """ check_is_fitted(self) + check_same_namespace(X, self, attribute="coeff_", method="predict_proba") xp, _, _ = get_namespace_and_device(X) X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) @@ -1344,30 +1351,8 @@ def predict(self, X): (n_samples, n_outputs). """ check_is_fitted(self) + check_same_namespace(X, self, attribute="coeff_", method="predict") xp, _, _ = get_namespace_and_device(X) X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) return super().predict(X) - - -def _to_numpy_cpu_y(y): - """Convert label array y to a 1D NumPy array on CPU.""" - # if y is None: - # raise ValueError("y cannot be None for a supervised estimator.") - - # PyTorch: handles CPU and CUDA tensors. - if hasattr(y, "detach") and hasattr(y, "cpu") and hasattr(y, "numpy"): - y = y.detach().cpu().numpy() - - # CuPy: GPU -> CPU NumPy. - elif hasattr(y, "get"): - y = y.get() - - # Generic fallback: NumPy, JAX CPU/device arrays, pandas, lists, etc. - else: - y = np.asarray(y) - - # Normalize shape for sklearn classifier utilities. - y = column_or_1d(y, warn=True) - - return y diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 39a6995..9277ff9 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -1,4 +1,3 @@ -import os import numpy as np import pytest @@ -9,11 +8,16 @@ from sklearn.metrics import accuracy_score, roc_auc_score from sklearn.model_selection import StratifiedKFold, train_test_split from sklearn.preprocessing import OneHotEncoder, StandardScaler +from sklearn.utils._array_api import ( + move_to, + yield_namespace_device_dtype_combinations, +) from sklearn.utils._testing import ( + _array_api_for_tests, skip_if_array_api_compat_not_configured, ) from sklearn.utils.estimator_checks import ( - check_array_api_input_and_values, + check_array_api_input, parametrize_with_checks, ) from ucimlrepo import fetch_ucirepo @@ -445,9 +449,7 @@ def test_classification_against_grafo(hidden_layer_sizes, n_classes, activation, @parametrize_with_checks([GFDLClassifier(), EnsembleGFDLClassifier()]) -def test_sklearn_api_conformance(estimator, check, request): - if "check_array_api" in request.node.name: - pytest.skip("Array API checks run separately") +def test_sklearn_api_conformance(estimator, check): check(estimator) @@ -836,9 +838,6 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): @skip_if_array_api_compat_not_configured -@pytest.mark.skipif( - os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." -) @parametrize_with_checks([GFDLClassifier()]) def test_sklearn_array_api_conformance(estimator, check,): with config_context(array_api_dispatch=True): @@ -854,28 +853,29 @@ def test_sklearn_array_api_conformance(estimator, check,): @skip_if_array_api_compat_not_configured -@pytest.mark.skipif( - os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." -) -@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) -@pytest.mark.parametrize("activation", ["identity", - "relu", - "tanh", - ]) -@pytest.mark.parametrize("weight_scheme", ["normal", - "zeros", - ]) -@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -def test_torch_array_api(hidden_layer_sizes, +@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) +@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) +@pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize( + "array_namespace, device_name, dtype_name", + yield_namespace_device_dtype_combinations(), +) +def test_xp_check_inputs(hidden_layer_sizes, direct_links, activation, weight_scheme, reg_alpha, + array_namespace, + device_name, + dtype_name, ): - """NumPy and Array API public attributes and predictions should be close""" - - pytest.importorskip("torch") + if (array_namespace not in ["numpy", "torch"] or + device_name not in [None, "cpu"] or + dtype_name not in [None, "float64"] + ): + pytest.skip("unsupported configuration") estimator = GFDLClassifier( reg_alpha=reg_alpha, @@ -886,80 +886,107 @@ def test_torch_array_api(hidden_layer_sizes, seed=42, ) with config_context(array_api_dispatch=True): - # Simulates from X, y = make_classification(n_samples=30, n_features=10) - # Unstable across backends for some weight schemes - check_array_api_input_and_values( + check_array_api_input( estimator.__class__.__name__, estimator, - array_namespace="torch", - device_name="cpu", - dtype_name="float64", + array_namespace=array_namespace, + device_name=device_name, + dtype_name=dtype_name, ) + @skip_if_array_api_compat_not_configured -@pytest.mark.skipif( - os.environ.get("SCIPY_ARRAY_API") != "1", reason="SCIPY_ARRAY_API not set to 1." -) -@pytest.mark.parametrize("hidden_layer_sizes", [(100,), (100, 100)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) -@pytest.mark.parametrize("reg_alpha", [None, 0.1, 0.5, 1]) -@pytest.mark.parametrize("n_classes", [2, 3]) -def test_torch_make_classificatino(hidden_layer_sizes, - direct_links, - activation, - weight_scheme, - reg_alpha, - n_classes, - ): - """Custom make_classification()-based NumPy and PyTorch matching""" - pytest.importorskip("torch") - import torch +@pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize( + "array_namespace, device_name, dtype_name", + yield_namespace_device_dtype_combinations(), +) +def test_xp_make_classification(hidden_layer_sizes, + direct_links, + activation, + weight_scheme, + reg_alpha, + array_namespace, + device_name, + dtype_name, + ): + # Custom make_regression()-based NumPy and array API matching + # This test is closely based on `_check_array_api_core` from sklearn. + # Here, we have better control of `atol`, `rtol`, and simulating data. + # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. + if (array_namespace not in ["numpy", "torch"] or + device_name not in [None, "cpu"] or + dtype_name not in [None, "float64"] + ): + pytest.skip("unsupported configuration") + xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) + + atol = 1e-3 + rtol = 1e-6 + n_samples = 500 + random_state = 42 estimator = GFDLClassifier - rng = np.random.default_rng(seed=42) - - pred_torch = [] - pred_numpy = [] - coef_torch = [] - coef_numpy = [] - for _ in range(3): - random_state = rng.integers(low=0, high=100, size=1)[0] - X_np, y_np = make_classification( - n_samples=250, - n_features=20, - n_informative=5, - n_classes=n_classes, - random_state=random_state, + + X, y = make_classification( + n_samples=n_samples, + n_features=20, + n_informative=5, + n_classes=2, + random_state=random_state, + ) + X = X.astype(dtype_name, copy=False) + X_xp = xp.asarray(X, device=device) + + # NumPy .fit + est = estimator( + reg_alpha=reg_alpha, + hidden_layer_sizes=hidden_layer_sizes, + activation=activation, + weight_scheme=weight_scheme, + direct_links=direct_links, + seed=random_state,) + with config_context(array_api_dispatch=False): + est.fit(X, y) + + # xp (array API) .fit + est_xp = clone(est) + with config_context(array_api_dispatch=True): + est_xp.fit(X_xp, y) + + # Check fitted attributes are close + array_attributes = { + key: value for key, value in vars(est).items() if isinstance(value, np.ndarray) + } + for attribute_name, attribute_value in array_attributes.items(): + est_xp_attr = getattr(est_xp, attribute_name) + est_xp_attr_np = move_to(est_xp_attr, xp=np, device="cpu") + assert_allclose( + attribute_value, + est_xp_attr_np, + err_msg=f"{attribute_name} not the same", + atol=atol, + rtol=rtol, ) - # NumPy .fit and .predict - model = estimator( - reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - activation=activation, - weight_scheme=weight_scheme, - direct_links=direct_links, - seed=random_state,) - model.fit(X_np, y_np) - y_pred_np = model.predict(X_np) - pred_numpy.append(y_pred_np) - coef_numpy.append(model.coeff_) - - # PyTorch .fit and .predict - X_torch = torch.asarray(X_np, device="cpu",) - y_torch = torch.asarray(y_np, device="cpu",) - model = estimator( - reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - activation=activation, - weight_scheme=weight_scheme, - direct_links=direct_links, - seed=random_state,) - model.fit(X_torch, y_torch) - y_pred_torch = model.predict(X_torch) - pred_torch.append(y_pred_torch) - coef_torch.append(model.coeff_) - - assert_allclose(pred_numpy, pred_torch) - assert_allclose(coef_numpy, coef_torch) \ No newline at end of file + # Check method outputs are close + methods = [ + "predict", + "predict_proba", + ] + for method_name in methods: + with config_context(array_api_dispatch=False): + result = getattr(est, method_name)(X) + with config_context(array_api_dispatch=True): + result_xp = getattr(est_xp, method_name)(X_xp) + result_xp_np = move_to(result_xp, xp=np, device="cpu") + assert_allclose( + result, + result_xp_np, + err_msg=f"{method_name}() did not the return the same result", + atol=atol, + rtol=rtol, + ) From 44b6c8eae741fd4c2e0d5df3240a90959ebb9188 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Tue, 20 Oct 2026 14:55:33 -0600 Subject: [PATCH 27/32] WIP,ENH: `partial_fit` for classifier * SI disclosure: I did not use SI for this * Main edit to to put `classes` in `predict` on the same namespace as `X` so that you can index on `xp.argmax` * Revise `xp_make_classification` test to allow for `partial_fit` results --- src/gfdl/model.py | 6 ++-- src/gfdl/tests/test_model.py | 57 ++++++++++++++++-------------------- 2 files changed, 29 insertions(+), 34 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index eb68e45..221ee76 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -549,7 +549,7 @@ def partial_fit(self, X, y, classes=None): if not hasattr(self, "classes_"): if classes is None: raise TypeError("Classes must not be None for first partial_fit call") - self.classes_ = classes + self.classes_ = unique_labels(classes) self.enc_ = OneHotEncoder(handle_unknown="ignore", sparse_output=False) self.enc_.fit(self.classes_.reshape(-1, 1)) @@ -583,8 +583,8 @@ def predict(self, X): X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) out = self.predict_proba(X) - xp, _, _device = get_namespace_and_device(self.coeff_) - y_hat = self.classes_[xp.argmax(out, axis=1)] + classes = move_to(self.classes_, xp=xp, device=device) + y_hat = classes[xp.argmax(out, axis=1)] return y_hat def predict_proba(self, X): diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 9277ff9..87d002a 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -613,7 +613,7 @@ def test_partial_fit( pf_model.partial_fit(Xb, yb) assert_allclose( - getattr(pf_model, attr), getattr(ff_model, attr), rtol=1e-5, atol=1e-3 + getattr(pf_model, attr), getattr(ff_model, attr), rtol=1e-5, atol=2e-3 ) @@ -837,21 +837,6 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): clf.fit(X, y) -@skip_if_array_api_compat_not_configured -@parametrize_with_checks([GFDLClassifier()]) -def test_sklearn_array_api_conformance(estimator, check,): - with config_context(array_api_dispatch=True): - try: - check(estimator) - except Exception as exc: - if "NotAnArray" in repr(exc): - pytest.xfail( - "Full check_estimator under array_api_dispatch=True is expected " - "to fail on sklearn's NotAnArray sentinel." - ) - raise - - @skip_if_array_api_compat_not_configured @pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) @pytest.mark.parametrize("direct_links", [True, False]) @@ -871,11 +856,6 @@ def test_xp_check_inputs(hidden_layer_sizes, device_name, dtype_name, ): - if (array_namespace not in ["numpy", "torch"] or - device_name not in [None, "cpu"] or - dtype_name not in [None, "float64"] - ): - pytest.skip("unsupported configuration") estimator = GFDLClassifier( reg_alpha=reg_alpha, @@ -896,11 +876,12 @@ def test_xp_check_inputs(hidden_layer_sizes, @skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) +@pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize("partial_fit", [False, True]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), @@ -913,24 +894,32 @@ def test_xp_make_classification(hidden_layer_sizes, array_namespace, device_name, dtype_name, + partial_fit, ): # Custom make_regression()-based NumPy and array API matching # This test is closely based on `_check_array_api_core` from sklearn. # Here, we have better control of `atol`, `rtol`, and simulating data. # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. - if (array_namespace not in ["numpy", "torch"] or - device_name not in [None, "cpu"] or - dtype_name not in [None, "float64"] - ): - pytest.skip("unsupported configuration") xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) - atol = 1e-3 - rtol = 1e-6 + atol = 1e-10 if dtype_name == "float64" else 2e-2 + rtol = 1e-13 if dtype_name == "float64" else 1e-4 + if partial_fit: + atol = 1e-7 if dtype_name == "float64" else 1e-0 + rtol = 1e-9 if dtype_name == "float64" else 1e-1 n_samples = 500 random_state = 42 estimator = GFDLClassifier + def run_partial_fit(X, y, model): + batch = 100 + classes = np.unique(y) + for start in range(0, X.shape[0], batch): + end = min(start + batch, X.shape[0]) + Xb = X[start:end, :] + yb = y[start:end] + model.partial_fit(Xb, yb, classes=classes) + X, y = make_classification( n_samples=n_samples, n_features=20, @@ -950,12 +939,18 @@ def test_xp_make_classification(hidden_layer_sizes, direct_links=direct_links, seed=random_state,) with config_context(array_api_dispatch=False): - est.fit(X, y) + if partial_fit: + run_partial_fit(X, y, est) + else: + est.fit(X, y) # xp (array API) .fit est_xp = clone(est) with config_context(array_api_dispatch=True): - est_xp.fit(X_xp, y) + if partial_fit: + run_partial_fit(X_xp, y, est_xp) + else: + est_xp.fit(X_xp, y) # Check fitted attributes are close array_attributes = { From 3cf5506ec4a10668d9f6479d822945818fe7b18f Mon Sep 17 00:00:00 2001 From: sdtemple Date: Tue, 20 Oct 2026 16:27:18 -0600 Subject: [PATCH 28/32] WIP,ENH: ensure classes for `partial_fit` are on `np` `cpu` * SI disclosure: I did not use SI * Consistent with the classifier `.fit`, we want `classes` to be on `np` and `cpu` for `OneHotEncoder` if this is the first pass when the encoder is fit. --- src/gfdl/model.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 221ee76..00ee603 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -549,7 +549,10 @@ def partial_fit(self, X, y, classes=None): if not hasattr(self, "classes_"): if classes is None: raise TypeError("Classes must not be None for first partial_fit call") - self.classes_ = unique_labels(classes) + classes = unique_labels( + move_to(classes, xp=np, device="cpu") + ) + self.classes_ = classes self.enc_ = OneHotEncoder(handle_unknown="ignore", sparse_output=False) self.enc_.fit(self.classes_.reshape(-1, 1)) From b2a3f2bcfcd8b2a27335a534e43dda6f553673d8 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Wed, 21 Oct 2026 10:11:49 -0600 Subject: [PATCH 29/32] WIP,ENH: `softmax` instead of `scipy.linalg.solve` * `xp.linalg.solve` documation for `partial_fit` in non-ensemble GFDL estimators * `softmax` from @activations.py achieves the same thing as in `partial_fit` custom written with `scipy.special.logsumexp` * These would eventually need to be addressed in an ensemble branch in the same way * Reverted to `rtol=1e-3` for `test_partial_fit` which had `pinv` instability in `EnsembleGFDLClassifier`. It looks like Tyler will `xfail` this case in #138 --- src/gfdl/model.py | 10 ++++------ src/gfdl/tests/test_model.py | 2 +- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 00ee603..128c122 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -24,7 +24,7 @@ from sklearn.utils.multiclass import check_classification_targets, unique_labels from sklearn.utils.validation import check_is_fitted, validate_data -from gfdl.activations import resolve_activation +from gfdl.activations import resolve_activation, softmax from gfdl.weights import resolve_weight @@ -264,7 +264,7 @@ def partial_fit(self, X, y): coef = xp.squeeze(coef, axis=-1) self.coeff_ = coef else: - # scipy.linalg.solve(self.A + reg_mat, self.B) + # xp.linalg.solve(self.A + reg_mat, self.B) # is equivalent to # ridge = Ridge( # alpha=self.reg_alpha, fit_intercept=False, solver='cholesky' @@ -531,7 +531,7 @@ def partial_fit(self, X, y, classes=None): ``Ridge()`` as we do in full fit. In ``partial_fit()`` we never have access to the full design matrix, we only persist the aggregate quantities ``D.T @ D``, ``D.T @ y``. Instead, for the regularized case, we directly solve - the system using ``scipy.linalg.solve(A + reg_mat, B)``, which is mathematically + the system using ``xp.linalg.solve(A + reg_mat, B)``, which is mathematically equivalent to scikit-learn's ``Ridge(solver='cholesky')``. One other difference between full fit and partial fit arises in the direct @@ -612,9 +612,7 @@ def predict_proba(self, X): X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) out = super().predict(X) - xp, _, _device = get_namespace_and_device(self.coeff_) - out = xp.exp(out - scipy.special.logsumexp(out, axis=1, keepdims=True)) - return out + return softmax(out) class EnsembleGFDL(BaseEstimator): diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 87d002a..0d418c6 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -613,7 +613,7 @@ def test_partial_fit( pf_model.partial_fit(Xb, yb) assert_allclose( - getattr(pf_model, attr), getattr(ff_model, attr), rtol=1e-5, atol=2e-3 + getattr(pf_model, attr), getattr(ff_model, attr), rtol=1e-5, atol=1e-3 ) From 2b87cae83aa00f04a1c57de490a96b6edd54ecbe Mon Sep 17 00:00:00 2001 From: sdtemple Date: Fri, 23 Oct 2026 17:27:19 -0600 Subject: [PATCH 30/32] WIP,ENH: more revisions after Emma's PR * SI disclosure: I did not use SI, but the basis of the changes from review of Emma's PR did use GPT 6.1 Sol * Use `LabelBinarizer` instead of `OneHotEncoder`. The `OneHotEncoder` is designed for features, not labels. See the `sklearn` documentation. The `LabelBinarizer` has Array API support, which means we don't have to do device transfers all of the time. * Add `ValueError` if only 1 class is specified in a classification problem. Had to modify the `np.zeros` to `np.arange() % 2` in an invalid activation/weight test so that the `ValueError` doesn't trigger first. * `rtol` spec in a `np.assert_allclose` test --- src/gfdl/model.py | 63 +++++++++++++++++++------------- src/gfdl/tests/test_model.py | 71 ++++-------------------------------- 2 files changed, 45 insertions(+), 89 deletions(-) diff --git a/src/gfdl/model.py b/src/gfdl/model.py index 128c122..b85fb4d 100644 --- a/src/gfdl/model.py +++ b/src/gfdl/model.py @@ -13,10 +13,11 @@ RegressorMixin, ) from sklearn.linear_model import Ridge -from sklearn.preprocessing import OneHotEncoder +from sklearn.preprocessing import LabelBinarizer, OneHotEncoder from sklearn.utils import column_or_1d from sklearn.utils._array_api import ( check_same_namespace, + get_namespace, get_namespace_and_device, move_to, ) @@ -471,19 +472,21 @@ def fit(self, X, y): xp, _, _ = get_namespace_and_device(X) X, Y = validate_data(self, X, y, dtype=[xp.float64, xp.float32]) - xp, _, device = get_namespace_and_device(X) - self.classes_ = move_to( - unique_labels(Y), - xp=xp, - device=device, - ) - Y = move_to(Y, xp=np, device="cpu") - # onehot y - # (this is necessary for everything beyond binary classification) - self.enc_ = OneHotEncoder(handle_unknown="ignore", sparse_output=False) - # shape: (n_samples, n_classes-1) - Y = self.enc_.fit_transform(Y.reshape(-1, 1)) + Y = column_or_1d(Y, warn=True) + check_classification_targets(Y) + + self.enc_ = LabelBinarizer(sparse_output=False) + Y = self.enc_.fit_transform(Y) + self.classes_ = self.enc_.classes_ + + # Preserve two columns for binary classification. + if self.classes_.shape[0] == 2: + xp_y, _ = get_namespace(Y) + Y = xp_y.concat((1 - Y, Y), axis=1) + + if self.classes_.shape[0] < 2: + raise ValueError("y must contain more than 1 class.") # call base fit method super().fit(X, Y) @@ -549,20 +552,30 @@ def partial_fit(self, X, y, classes=None): if not hasattr(self, "classes_"): if classes is None: raise TypeError("Classes must not be None for first partial_fit call") - classes = unique_labels( - move_to(classes, xp=np, device="cpu") - ) - self.classes_ = classes - self.enc_ = OneHotEncoder(handle_unknown="ignore", sparse_output=False) - self.enc_.fit(self.classes_.reshape(-1, 1)) - - if not set(unique_labels(Y)) <= set(self.classes_): + xp_y, _, device_y_ = get_namespace_and_device(Y) + classes = move_to(classes, xp=xp_y, device=device_y_) + self.enc_ = LabelBinarizer(sparse_output=False).fit(classes) + self.classes_ = self.enc_.classes_ + if self.classes_.shape[0] < 2: + raise ValueError("classes must contain more than 1 class") + + # Check membership without Python sets or converting labels to NumPy. + xp_y, _, device_y = get_namespace_and_device(self.classes_) + Y = move_to(Y, xp=xp_y, device=device_y) + matches = xp_y.expand_dims(Y, axis=1) == xp_y.expand_dims( + self.classes_, axis=0 + ) + if not xp_y.all(xp_y.any(matches, axis=1)): raise ValueError( - f"Expected only labels in classes_ = {list(self.classes_)!r}, " - f"but got {unique_labels(Y)!r}." + f"Expected only labels in classes_ = {self.classes_!r}," + f"but got {unique_labels(y)!r}." ) Y = self.enc_.transform(Y.reshape(-1, 1)) + # Preserve two output columns for binary classification. + if self.classes_.shape[0] == 2: + Y = xp_y.concat((1 - Y, Y), axis=1) + super().partial_fit(X, Y) return self @@ -582,7 +595,7 @@ def predict(self, X): """ check_is_fitted(self) check_same_namespace(X, self, attribute="coeff_", method="predict") - xp, _, _ = get_namespace_and_device(X) + xp, _, device = get_namespace_and_device(X) X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) out = self.predict_proba(X) @@ -608,7 +621,7 @@ def predict_proba(self, X): """ check_is_fitted(self) check_same_namespace(X, self, attribute="coeff_", method="predict_proba") - xp, _, _ = get_namespace_and_device(X) + xp, _ = get_namespace(X) X = validate_data(self, X, reset=False, dtype=[xp.float64, xp.float32]) out = super().predict(X) diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 0d418c6..2618b6a 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -17,7 +17,6 @@ skip_if_array_api_compat_not_configured, ) from sklearn.utils.estimator_checks import ( - check_array_api_input, parametrize_with_checks, ) from ucimlrepo import fetch_ucirepo @@ -342,7 +341,7 @@ def test_soft_and_hard_can_differ(alpha): @pytest.mark.parametrize("Classifier", [GFDLClassifier, EnsembleGFDLClassifier]) def test_invalid_activation_weight(Classifier): X = np.zeros((30, 4)) - y = np.zeros((30,)) + y = np.arange(30) % 2 invalid_act = Classifier(hidden_layer_sizes=(100,), activation="bogus_activation", weight_scheme="uniform") @@ -364,7 +363,7 @@ def test_invalid_alpha(Classifier): # so we need to call fit() for error handling to kick in: # https://scikit-learn.org/stable/developers/develop.html#developing-scikit-learn-estimators X = np.zeros((30, 4)) - y = np.zeros((30,)) + y = np.arange(30) % 2 bad_est = Classifier(hidden_layer_sizes=100, activation="identity", weight_scheme="uniform", @@ -509,7 +508,9 @@ def test_rtol_classifier(reg_alpha, rtol, expected_acc, expected_roc): roc_cur = roc_auc_score(y_test, y_hat_cur_proba, multi_class="ovo") np.testing.assert_allclose(acc_cur, expected_acc) - np.testing.assert_allclose(roc_cur, expected_roc) + np.testing.assert_allclose(roc_cur, expected_roc, + rtol=1e-4 + ) @pytest.mark.parametrize("reg_alpha, rtol, expected_acc, expected_roc", [ @@ -837,51 +838,12 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): clf.fit(X, y) -@skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(4,), (4, 4)]) -@pytest.mark.parametrize("direct_links", [True, False]) -@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) -@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) -@pytest.mark.parametrize("reg_alpha", [None, 0.1]) -@pytest.mark.parametrize( - "array_namespace, device_name, dtype_name", - yield_namespace_device_dtype_combinations(), -) -def test_xp_check_inputs(hidden_layer_sizes, - direct_links, - activation, - weight_scheme, - reg_alpha, - array_namespace, - device_name, - dtype_name, - ): - - estimator = GFDLClassifier( - reg_alpha=reg_alpha, - hidden_layer_sizes=hidden_layer_sizes, - activation=activation, - weight_scheme=weight_scheme, - direct_links=direct_links, - seed=42, - ) - with config_context(array_api_dispatch=True): - check_array_api_input( - estimator.__class__.__name__, - estimator, - array_namespace=array_namespace, - device_name=device_name, - dtype_name=dtype_name, - ) - - @skip_if_array_api_compat_not_configured @pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) @pytest.mark.parametrize("direct_links", [True, False]) @pytest.mark.parametrize("activation", ACTIVATIONS.keys()) @pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) @pytest.mark.parametrize("reg_alpha", [None, 0.1]) -@pytest.mark.parametrize("partial_fit", [False, True]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), @@ -894,7 +856,6 @@ def test_xp_make_classification(hidden_layer_sizes, array_namespace, device_name, dtype_name, - partial_fit, ): # Custom make_regression()-based NumPy and array API matching # This test is closely based on `_check_array_api_core` from sklearn. @@ -904,22 +865,10 @@ def test_xp_make_classification(hidden_layer_sizes, atol = 1e-10 if dtype_name == "float64" else 2e-2 rtol = 1e-13 if dtype_name == "float64" else 1e-4 - if partial_fit: - atol = 1e-7 if dtype_name == "float64" else 1e-0 - rtol = 1e-9 if dtype_name == "float64" else 1e-1 n_samples = 500 random_state = 42 estimator = GFDLClassifier - def run_partial_fit(X, y, model): - batch = 100 - classes = np.unique(y) - for start in range(0, X.shape[0], batch): - end = min(start + batch, X.shape[0]) - Xb = X[start:end, :] - yb = y[start:end] - model.partial_fit(Xb, yb, classes=classes) - X, y = make_classification( n_samples=n_samples, n_features=20, @@ -939,18 +888,12 @@ def run_partial_fit(X, y, model): direct_links=direct_links, seed=random_state,) with config_context(array_api_dispatch=False): - if partial_fit: - run_partial_fit(X, y, est) - else: - est.fit(X, y) + est.fit(X, y) # xp (array API) .fit est_xp = clone(est) with config_context(array_api_dispatch=True): - if partial_fit: - run_partial_fit(X_xp, y, est_xp) - else: - est_xp.fit(X_xp, y) + est_xp.fit(X_xp, y) # Check fitted attributes are close array_attributes = { From c10ca08f247406a6471c753a1a560364d7f688f8 Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 26 Oct 2026 15:28:17 -0600 Subject: [PATCH 31/32] WIP,TST: modified `test_model.py` for Array API * SI disclosure: I did not use SI * Simplify/shorten some code lines in classifier fit test * Use `xpx.testing.assert_close` instead * Include a partial fit Array API test for classifier which only checks `predict` and `predict_proba` --- src/gfdl/tests/test_model.py | 107 +++++++++++++++++++++++++++-------- 1 file changed, 82 insertions(+), 25 deletions(-) diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index 2618b6a..bf565c1 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -5,10 +5,12 @@ from sklearn import config_context from sklearn.base import clone from sklearn.datasets import load_breast_cancer, load_digits, make_classification +from sklearn.externals import array_api_extra as xpx from sklearn.metrics import accuracy_score, roc_auc_score from sklearn.model_selection import StratifiedKFold, train_test_split from sklearn.preprocessing import OneHotEncoder, StandardScaler from sklearn.utils._array_api import ( + _atol_for_type, move_to, yield_namespace_device_dtype_combinations, ) @@ -839,16 +841,25 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): @skip_if_array_api_compat_not_configured -@pytest.mark.parametrize("hidden_layer_sizes", [(2,), (2, 2)]) -@pytest.mark.parametrize("direct_links", [True, False]) -@pytest.mark.parametrize("activation", ACTIVATIONS.keys()) -@pytest.mark.parametrize("weight_scheme", WEIGHTS.keys()) -@pytest.mark.parametrize("reg_alpha", [None, 0.1]) +@pytest.mark.parametrize( + "reg_alpha, hidden_layer_sizes, activation, weight_scheme, direct_links", + [ + (2, (10,), "relu", "normal", True), + (2, (10, 10), "sigmoid", "zeros", True), + (2, (10,), "softmax", "lecun_normal", True), + (2, (10,), "tanh", "he_normal", False), + (None, (10,), "softmin", "glorot_normal", True), + (None, (10, 10), "log_softmax", "uniform", False), + (None, (10,), "log_sigmoid", "lecun_uniform", True), + (2, (10,), "identity", "glorot_uniform", True), + (2, (10,), "identity", "he_uniform", True), + ] +) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), ) -def test_xp_make_classification(hidden_layer_sizes, +def test_xp_classifier_fit(hidden_layer_sizes, direct_links, activation, weight_scheme, @@ -863,18 +874,11 @@ def test_xp_make_classification(hidden_layer_sizes, # https://github.com/scikit-learn/scikit-learn/pull/33215#issuecomment-5953082117. xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) - atol = 1e-10 if dtype_name == "float64" else 2e-2 - rtol = 1e-13 if dtype_name == "float64" else 1e-4 - n_samples = 500 - random_state = 42 + atol = _atol_for_type(dtype_name) * 10 estimator = GFDLClassifier - X, y = make_classification( - n_samples=n_samples, - n_features=20, - n_informative=5, - n_classes=2, - random_state=random_state, + X, y = make_classification(n_samples=500, n_features=20, + n_informative=5, n_classes=2, random_state=0, ) X = X.astype(dtype_name, copy=False) X_xp = xp.asarray(X, device=device) @@ -886,7 +890,7 @@ def test_xp_make_classification(hidden_layer_sizes, activation=activation, weight_scheme=weight_scheme, direct_links=direct_links, - seed=random_state,) + seed=0,) with config_context(array_api_dispatch=False): est.fit(X, y) @@ -902,12 +906,10 @@ def test_xp_make_classification(hidden_layer_sizes, for attribute_name, attribute_value in array_attributes.items(): est_xp_attr = getattr(est_xp, attribute_name) est_xp_attr_np = move_to(est_xp_attr, xp=np, device="cpu") - assert_allclose( - attribute_value, - est_xp_attr_np, + xpx.testing.assert_close( + attribute_value, est_xp_attr_np, err_msg=f"{attribute_name} not the same", atol=atol, - rtol=rtol, ) # Check method outputs are close @@ -921,10 +923,65 @@ def test_xp_make_classification(hidden_layer_sizes, with config_context(array_api_dispatch=True): result_xp = getattr(est_xp, method_name)(X_xp) result_xp_np = move_to(result_xp, xp=np, device="cpu") - assert_allclose( - result, - result_xp_np, + xpx.testing.assert_close( + result, result_xp_np, err_msg=f"{method_name}() did not the return the same result", atol=atol, - rtol=rtol, ) + + +@skip_if_array_api_compat_not_configured +@pytest.mark.parametrize("reg_alpha", [None, 2]) +@pytest.mark.parametrize("hidden_layer_sizes", [(10,), (10, 10)]) +@pytest.mark.parametrize( + "array_namespace, device_name, dtype_name", + yield_namespace_device_dtype_combinations(), +) +def test_xp_classifier_partial_fit(hidden_layer_sizes, + reg_alpha, + array_namespace, + device_name, + dtype_name, + ): + + atol = _atol_for_type(dtype_name) * 10 + batch = 2 + + xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) + X_np, y_np = make_classification( + n_samples=6, n_features=20, n_informative=5, random_state=0, n_classes=2 + ) + X_np = X_np.astype(dtype_name, copy=False) + X_xp = xp.asarray(X_np, device=device) + + params = { + "reg_alpha": reg_alpha, + "hidden_layer_sizes": hidden_layer_sizes, + "seed": 0 + } + estimator = GFDLClassifier + + predictions = [] + for X, dispatch in [(X_np, False), (X_xp, True)]: + with config_context(array_api_dispatch=dispatch): + model = estimator(**params) + for start in range(0, X.shape[0], batch): + model.partial_fit( + X[start:start + batch, ...], y_np[start:start + batch], + classes=np.unique(y_np) + ) + predictions.extend( + (move_to(model.predict(X), xp=np, device="cpu"), + move_to(model.predict_proba(X), xp=np, device="cpu")) + ) + + xpx.testing.assert_close( + predictions[0], predictions[2], + err_msg="predict() did not the return the same result", + atol=atol, + ) + xpx.testing.assert_close( + predictions[1], predictions[3], + err_msg="predict() did not the return the same result", + atol=atol, + ) From d1d67087469e23323b7f9bbf67682582a21c6f9d Mon Sep 17 00:00:00 2001 From: sdtemple Date: Mon, 26 Oct 2026 17:41:35 -0600 Subject: [PATCH 32/32] WIP,TST: include multiclass in `xp` classifier test * SI disclosure: I did not use SI * Include `n_classes` as a parameter for Array API classifier test --- src/gfdl/tests/test_model.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/gfdl/tests/test_model.py b/src/gfdl/tests/test_model.py index bf565c1..50e4dd4 100644 --- a/src/gfdl/tests/test_model.py +++ b/src/gfdl/tests/test_model.py @@ -855,6 +855,7 @@ def test_gh_85_classifiers(hidden_layer_sizes, classifier): (2, (10,), "identity", "he_uniform", True), ] ) +@pytest.mark.parametrize("n_classes", [2, 5]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), @@ -864,6 +865,7 @@ def test_xp_classifier_fit(hidden_layer_sizes, activation, weight_scheme, reg_alpha, + n_classes, array_namespace, device_name, dtype_name, @@ -878,7 +880,7 @@ def test_xp_classifier_fit(hidden_layer_sizes, estimator = GFDLClassifier X, y = make_classification(n_samples=500, n_features=20, - n_informative=5, n_classes=2, random_state=0, + n_informative=5, n_classes=n_classes, random_state=0, ) X = X.astype(dtype_name, copy=False) X_xp = xp.asarray(X, device=device) @@ -933,12 +935,14 @@ def test_xp_classifier_fit(hidden_layer_sizes, @skip_if_array_api_compat_not_configured @pytest.mark.parametrize("reg_alpha", [None, 2]) @pytest.mark.parametrize("hidden_layer_sizes", [(10,), (10, 10)]) +@pytest.mark.parametrize("n_classes", [2, 5]) @pytest.mark.parametrize( "array_namespace, device_name, dtype_name", yield_namespace_device_dtype_combinations(), ) def test_xp_classifier_partial_fit(hidden_layer_sizes, reg_alpha, + n_classes, array_namespace, device_name, dtype_name, @@ -949,7 +953,8 @@ def test_xp_classifier_partial_fit(hidden_layer_sizes, xp, device = _array_api_for_tests(array_namespace, device_name, dtype_name) X_np, y_np = make_classification( - n_samples=6, n_features=20, n_informative=5, random_state=0, n_classes=2 + n_samples=6, n_features=20, n_informative=5, random_state=0, + n_classes=n_classes ) X_np = X_np.astype(dtype_name, copy=False) X_xp = xp.asarray(X_np, device=device)