diff --git a/HISTORY.rst b/HISTORY.rst index 553a0fd43..9ef3c2297 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -4,7 +4,8 @@ History ##### (##########) ------------------ - +* Refactor MapieRegressor and ConformityScore to add the possibility to use X in ConformityScore. +* Separate the handling of the estimator from MapieRegressor into a new class called EnsembleEstimator. * Fix an unfixed random state in one of the classification tests 0.6.5 (2023-06-06) diff --git a/examples/regression/1-quickstart/plot_compare_conformity_scores.py b/examples/regression/1-quickstart/plot_compare_conformity_scores.py index 7d345b634..e4b79c701 100644 --- a/examples/regression/1-quickstart/plot_compare_conformity_scores.py +++ b/examples/regression/1-quickstart/plot_compare_conformity_scores.py @@ -39,7 +39,7 @@ from mapie.metrics import regression_coverage_score from mapie.regression import MapieRegressor -np.random.seed(0) +random_state = 42 # Parameters features = [ @@ -50,7 +50,7 @@ "GarageArea", ] alpha = 0.05 -rf_kwargs = {"n_estimators": 10, "random_state": 0} +rf_kwargs = {"n_estimators": 10, "random_state": random_state} model = RandomForestRegressor(**rf_kwargs) ############################################################################## @@ -66,7 +66,7 @@ X, y = fetch_openml(name="house_prices", return_X_y=True) X_train, X_test, y_train, y_test = train_test_split( - X[features], y, test_size=0.2 + X[features], y, test_size=0.2, random_state=random_state ) ############################################################################## @@ -87,9 +87,11 @@ ############################################################################## # First, train model with # :class:`~mapie.conformity_scores.AbsoluteConformityScore`. -mapie = MapieRegressor(model) +mapie = MapieRegressor(model, random_state=random_state) mapie.fit(X_train, y_train) -y_pred_absconfscore, y_pis_absconfscore = mapie.predict(X_test, alpha=alpha) +y_pred_absconfscore, y_pis_absconfscore = mapie.predict( + X_test, alpha=alpha, ensemble=True +) coverage_absconfscore = regression_coverage_score( y_test, y_pis_absconfscore[:, 0, 0], y_pis_absconfscore[:, 1, 0] @@ -118,10 +120,12 @@ def get_yerr(y_pred, y_pis): ############################################################################## # Then, train the model with # :class:`~mapie.conformity_scores.GammaConformityScore`. -mapie = MapieRegressor(model, conformity_score=GammaConformityScore()) +mapie = MapieRegressor( + model, conformity_score=GammaConformityScore(), random_state=random_state +) mapie.fit(X_train, y_train) y_pred_gammaconfscore, y_pis_gammaconfscore = mapie.predict( - X_test, alpha=[alpha] + X_test, alpha=[alpha], ensemble=True ) coverage_gammaconfscore = regression_coverage_score( diff --git a/mapie/conformity_scores/conformity_scores.py b/mapie/conformity_scores/conformity_scores.py index c5f689746..ef4a79ade 100644 --- a/mapie/conformity_scores/conformity_scores.py +++ b/mapie/conformity_scores/conformity_scores.py @@ -1,8 +1,11 @@ from abc import ABCMeta, abstractmethod import numpy as np +from typing import Tuple +from mapie._compatibility import np_nanquantile from mapie._typing import ArrayLike, NDArray +from mapie.estimator.interface import EnsembleEstimator class ConformityScore(metaclass=ABCMeta): @@ -31,7 +34,9 @@ class ConformityScore(metaclass=ABCMeta): - ``get_signed_conformity_scores`` The following equality must be verified: ``self.get_estimation_distribution( - y_pred, self.get_conformity_scores(y, y_pred) + X, + y_pred, + self.get_conformity_scores(X, y, y_pred) ) == y`` It should be specified if ``consistency_check==True``. @@ -51,6 +56,7 @@ def __init__( @abstractmethod def get_signed_conformity_scores( self, + X: ArrayLike, y: ArrayLike, y_pred: ArrayLike, ) -> NDArray: @@ -63,47 +69,62 @@ def get_signed_conformity_scores( Parameters ---------- - y: NDArray - Observed values. + X: ArrayLike of shape (n_samples, n_features) + Observed feature values. + + y: ArrayLike of shape (n_samples,) + Observed target values. - y_pred: NDArray - Predicted values. + y_pred: ArrayLike of shape (n_samples,) + Predicted target values. Returns ------- - NDArray - Unsigned conformity scores. + NDArray of shape (n_samples,) + Signed conformity scores. """ @abstractmethod def get_estimation_distribution( self, + X: ArrayLike, y_pred: ArrayLike, - conformity_scores: ArrayLike, + conformity_scores: ArrayLike ) -> NDArray: """ Placeholder for ``get_estimation_distribution``. Subclasses should implement this method! Compute samples of the estimation distribution from the predicted - values and the conformity scores. + targets and ``conformity_scores`` that can be either the conformity + scores or the quantile of the conformity scores. Parameters ---------- - y_pred: NDArray - Predicted values. + X: ArrayLike of shape (n_samples, n_features) + Observed feature values. - conformity_scores: NDArray - Conformity scores. + y_pred: ArrayLike + The shape is either (n_samples, n_references): when the + method is called in ``get_bounds`` it needs a prediction per train + sample for each test sample to compute the bounds. + Or (n_samples,): when it is called in ``check_consistency`` + + conformity_scores: ArrayLike + The shape is either (n_samples, 1) when it is the + conformity scores themselves or (1, n_alpha) when it is only the + quantile of the conformity scores. Returns ------- - NDArray + NDArray of shape (n_samples, n_alpha) or + (n_samples, n_references) according to the shape of ``y_pred`` Observed values. """ def check_consistency( self, + X: ArrayLike, y: ArrayLike, y_pred: ArrayLike, conformity_scores: ArrayLike, @@ -114,16 +135,24 @@ def check_consistency( The following equality should be verified: ``self.get_estimation_distribution( - y_pred, self.get_conformity_scores(y, y_pred) + X, + y_pred, + self.get_conformity_scores(X, y, y_pred) ) == y`` Parameters ---------- - y: NDArray - Observed values. + X: ArrayLike of shape (n_samples, n_features) + Observed feature values. + + y: ArrayLike of shape (n_samples,) + Observed target values. - y_pred: NDArray - Predicted values. + y_pred: ArrayLike of shape (n_samples,) + Predicted target values. + + conformity_scores: ArrayLike of shape (n_samples,) + Conformity scores. Raises ------ @@ -131,7 +160,7 @@ def check_consistency( If the two methods are not consistent. """ score_distribution = self.get_estimation_distribution( - y_pred, conformity_scores + X, y_pred, conformity_scores ) abs_conformity_scores = np.abs(np.subtract(score_distribution, y)) max_conf_score = np.max(abs_conformity_scores) @@ -141,8 +170,8 @@ def check_consistency( "get_estimation_distribution of the ConformityScore class " "are not consistent. " "The following equation must be verified: " - "self.get_estimation_distribution(y_pred, " - "self.get_conformity_scores(y, y_pred)) == y. " # noqa: E501 + "self.get_estimation_distribution(X, y_pred, " + "self.get_conformity_scores(X, y, y_pred)) == y" # noqa: E501 f"The maximum conformity score is {max_conf_score}." "The eps attribute may need to be increased if you are " "sure that the two methods are consistent." @@ -150,6 +179,7 @@ def check_consistency( def get_conformity_scores( self, + X: ArrayLike, y: ArrayLike, y_pred: ArrayLike, ) -> NDArray: @@ -158,20 +188,154 @@ def get_conformity_scores( Parameters ---------- - y: NDArray - Observed values. + X: NDArray of shape (n_samples, n_features) + Observed feature values. + + y: NDArray of shape (n_samples,) + Observed target values. - y_pred: NDArray - Predicted values. + y_pred: NDArray of shape (n_samples,) + Predicted target values. Returns ------- - NDArray + NDArray of shape (n_samples,) Conformity scores. """ - conformity_scores = self.get_signed_conformity_scores(y, y_pred) + conformity_scores = self.get_signed_conformity_scores(X, y, y_pred) if self.consistency_check: - self.check_consistency(y, y_pred, conformity_scores) + self.check_consistency(X, y, y_pred, conformity_scores) if self.sym: conformity_scores = np.abs(conformity_scores) return conformity_scores + + @staticmethod + def get_quantile( + conformity_scores: NDArray, + alpha_np: NDArray, + axis: int, + method: str + ) -> NDArray: + """ + Compute the alpha quantile of the conformity scores or the conformity + scores aggregated with the predictions. + + Parameters + ---------- + conformity_scores: NDArray of shape (n_samples,) or + (n_samples, n_references) + Values from which the quantile is computed, it can be the + conformity scores or the conformity scores aggregated with + the predictions. + + alpha_np: NDArray of shape (n_alpha,) + NDArray of floats between ``0`` and ``1``, represents the + uncertainty of the confidence interval. + + axis: int + The axis from which to compute the quantile. + + method: str + ``"higher"`` or ``"lower"`` the method to compute the quantile. + + Returns + ------- + NDArray of shape (1, n_alpha) or (n_samples, n_alpha) + The quantile of the conformity scores. + """ + quantile = np.column_stack([ + np_nanquantile( + conformity_scores.astype(float), + _alpha, + axis=axis, + method=method + ) + for _alpha in alpha_np + ]) + return quantile + + def get_bounds( + self, + X: ArrayLike, + estimator: EnsembleEstimator, + conformity_scores: NDArray, + alpha_np: NDArray, + ensemble: bool, + method: str + ) -> Tuple[NDArray, NDArray, NDArray]: + """ + Compute bounds of the prediction intervals from the observed values, + the estimator of type ``EnsembleEstimator`` and the conformity scores. + + Parameters + ---------- + X: ArrayLike of shape (n_samples, n_features) + Observed feature values. + + estimator: EnsembleEstimator + Estimator that is fitted to predict y from X. + + conformity_scores: ArrayLike of shape (n_samples,) + Conformity scores. + + alpha_np: NDArray of shape (n_alpha,) + NDArray of floats between ``0`` and ``1``, represents the + uncertainty of the confidence interval. + + ensemble: bool + Boolean determining whether the predictions are ensembled or not. + + method: str + Method to choose for prediction interval estimates. + The ``"plus"`` method implies that the quantile is calculated + after estimating the bounds, whereas the other methods + (among the ``"naive"``, ``"base"`` or ``"minmax"`` methods, + for example) do the opposite. + + Returns + ------- + Tuple[NDArray, NDArray, NDArray] + - The predictions itself. (y_pred) of shape (n_samples,). + - The lower bounds of the prediction intervals of shape + (n_samples, n_alpha). + - The upper bounds of the prediction intervals of shape + (n_samples, n_alpha). + """ + y_pred, y_pred_low, y_pred_up = estimator.predict(X, ensemble) + signed = -1 if self.sym else 1 + + if method == "plus": + alpha_low = alpha_np if self.sym else alpha_np / 2 + alpha_up = 1 - alpha_np if self.sym else 1 - alpha_np / 2 + + conformity_scores_low = self.get_estimation_distribution( + X, y_pred_low, signed * conformity_scores + ) + conformity_scores_up = self.get_estimation_distribution( + X, y_pred_up, conformity_scores + ) + bound_low = self.get_quantile( + conformity_scores_low, alpha_low, axis=1, method="lower" + ) + bound_up = self.get_quantile( + conformity_scores_up, alpha_up, axis=1, method="higher" + ) + else: + quantile_search = "higher" if self.sym else "lower" + alpha_low = 1 - alpha_np if self.sym else alpha_np / 2 + alpha_up = 1 - alpha_np if self.sym else 1 - alpha_np / 2 + + quantile_low = self.get_quantile( + conformity_scores, alpha_low, axis=0, method=quantile_search + ) + quantile_up = self.get_quantile( + conformity_scores, alpha_up, axis=0, method="higher" + ) + bound_low = self.get_estimation_distribution( + X, y_pred_low, signed * quantile_low + ) + bound_up = self.get_estimation_distribution( + X, y_pred_up, quantile_up + ) + + return y_pred, bound_low, bound_up diff --git a/mapie/conformity_scores/residual_conformity_scores.py b/mapie/conformity_scores/residual_conformity_scores.py index 033baa594..8dbff9fd4 100644 --- a/mapie/conformity_scores/residual_conformity_scores.py +++ b/mapie/conformity_scores/residual_conformity_scores.py @@ -24,6 +24,7 @@ def __init__( def get_signed_conformity_scores( self, + X: ArrayLike, y: ArrayLike, y_pred: ArrayLike, ) -> NDArray: @@ -36,14 +37,18 @@ def get_signed_conformity_scores( def get_estimation_distribution( self, + X: ArrayLike, y_pred: ArrayLike, - conformity_scores: ArrayLike, + conformity_scores: ArrayLike ) -> NDArray: """ Compute samples of the estimation distribution from the predicted values and the conformity scores, from the following formula: signed conformity score = y - y_pred <=> y = y_pred + signed conformity score + + ``conformity_scores`` can be either the conformity scores or + the quantile of the conformity scores. """ return np.add(y_pred, conformity_scores) @@ -89,22 +94,21 @@ def _check_predicted_data( "in conformity with the Gamma distribution support." ) + @staticmethod def _all_strictly_positive( - self, y: ArrayLike, ) -> bool: - if np.any(np.less_equal(y, 0)): - return False - return True + return not np.any(np.less_equal(y, 0)) def get_signed_conformity_scores( self, + X: ArrayLike, y: ArrayLike, y_pred: ArrayLike, ) -> NDArray: """ - Compute samples of the estimation distribution from the predicted - values and the conformity scores, from the following formula: + Compute the signed conformity scores from the observed values + and the predicted ones, from the following formula: signed conformity score = (y - y_pred) / y_pred """ self._check_observed_data(y) @@ -113,14 +117,18 @@ def get_signed_conformity_scores( def get_estimation_distribution( self, + X: ArrayLike, y_pred: ArrayLike, - conformity_scores: ArrayLike, + conformity_scores: ArrayLike ) -> NDArray: """ Compute samples of the estimation distribution from the predicted values and the conformity scores, from the following formula: signed conformity score = (y - y_pred) / y_pred <=> y = y_pred * (1 + signed conformity score) + + ``conformity_scores`` can be either the conformity scores or + the quantile of the conformity scores. """ self._check_predicted_data(y_pred) return np.multiply(y_pred, np.add(1, conformity_scores)) diff --git a/mapie/estimator/__init__.py b/mapie/estimator/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/mapie/estimator/estimator.py b/mapie/estimator/estimator.py new file mode 100644 index 000000000..33bda7b3d --- /dev/null +++ b/mapie/estimator/estimator.py @@ -0,0 +1,509 @@ +from __future__ import annotations + +from typing import List, Optional, Tuple, Union, cast + +import numpy as np +from joblib import Parallel, delayed +from sklearn.base import RegressorMixin, clone +from sklearn.model_selection import BaseCrossValidator, ShuffleSplit +from sklearn.utils import _safe_indexing +from sklearn.utils.validation import (_num_samples, check_is_fitted) + +from mapie._typing import ArrayLike, NDArray +from mapie.aggregation_functions import aggregate_all, phi2D +from mapie.utils import (check_nan_in_aposteriori_prediction, + fit_estimator) +from mapie.estimator.interface import EnsembleEstimator + + +class EnsembleRegressor(EnsembleEstimator): + """ + This class implements methods to handle the training and usage of the + estimator. This estimator can be unique or composed by cross validated + estimators. + + Parameters + ---------- + estimator: Optional[RegressorMixin] + Any regressor with scikit-learn API + (i.e. with ``fit`` and ``predict`` methods). + If ``None``, estimator defaults to a ``LinearRegression`` instance. + + By default ``None``. + + method: str + Method to choose for prediction interval estimates. + Choose among: + + - ``"naive"``, based on training set conformity scores, + - ``"base"``, based on validation sets conformity scores, + - ``"plus"``, based on validation conformity scores and + testing predictions, + - ``"minmax"``, based on validation conformity scores and + testing predictions (min/max among cross-validation clones). + + By default ``"plus"``. + + cv: Optional[Union[int, str, BaseCrossValidator]] + The cross-validation strategy for computing conformity scores. + It directly drives the distinction between jackknife and cv variants. + Choose among: + + - ``None``, to use the default 5-fold cross-validation + - integer, to specify the number of folds. + If equal to ``-1``, equivalent to + ``sklearn.model_selection.LeaveOneOut()``. + - CV splitter: any ``sklearn.model_selection.BaseCrossValidator`` + Main variants are: + - ``sklearn.model_selection.LeaveOneOut`` (jackknife), + - ``sklearn.model_selection.KFold`` (cross-validation), + - ``subsample.Subsample`` object (bootstrap). + - ``"split"``, does not involve cross-validation but a division + of the data into training and calibration subsets. The splitter + used is the following: ``sklearn.model_selection.ShuffleSplit``. + - ``"prefit"``, assumes that ``estimator`` has been fitted already, + and the ``method`` parameter is ignored. + All data provided in the ``fit`` method is then used + for computing conformity scores only. + At prediction time, quantiles of these conformity scores are used + to provide a prediction interval with fixed width. + The user has to take care manually that data for model fitting and + conformity scores estimate are disjoint. + + By default ``None``. + + test_size: Optional[Union[int, float]] + If ``float``, should be between ``0.0`` and ``1.0`` and represent the + proportion of the dataset to include in the test split. If ``int``, + represents the absolute number of test samples. If ``None``, + it will be set to ``0.1``. + + If cv is not ``"split"``, ``test_size`` is ignored. + + By default ``None``. + + n_jobs: Optional[int] + Number of jobs for parallel processing using joblib + via the "locky" backend. + If ``-1`` all CPUs are used. + If ``1`` is given, no parallel computing code is used at all, + which is useful for debugging. + For ``n_jobs`` below ``-1``, ``(n_cpus + 1 - n_jobs)`` are used. + ``None`` is a marker for `unset` that will be interpreted as + ``n_jobs=1`` (sequential execution). + + By default ``None``. + + agg_function: Optional[str] + Determines how to aggregate predictions from perturbed models, both at + training and prediction time. + + If ``None``, it is ignored except if ``cv`` class is ``Subsample``, + in which case an error is raised. + If ``"mean"`` or ``"median"``, returns the mean or median of the + predictions computed from the out-of-folds models. + Note: if you plan to set the ``ensemble`` argument to ``True`` in the + ``predict`` method, you have to specify an aggregation function. + Otherwise an error would be raised. + + The Jackknife+ interval can be interpreted as an interval around the + median prediction, and is guaranteed to lie inside the interval, + unlike the single estimator predictions. + + When the cross-validation strategy is ``Subsample`` (i.e. for the + Jackknife+-after-Bootstrap method), this function is also used to + aggregate the training set in-sample predictions. + + If ``cv`` is ``"prefit"`` or ``"split"``, ``agg_function`` is ignored. + + By default ``"mean"``. + + verbose: int + The verbosity level, used with joblib for multiprocessing. + The frequency of the messages increases with the verbosity level. + If it more than ``10``, all iterations are reported. + Above ``50``, the output is sent to stdout. + + By default ``0``. + + random_state: Optional[Union[int, RandomState]] + Pseudo random number generator state used for random sampling. + Pass an int for reproducible output across multiple function calls. + + By default ``None``. + + Attributes + ---------- + single_estimator_: sklearn.RegressorMixin + Estimator fitted on the whole training set. + + estimators_: list + List of out-of-folds estimators. + + k_: ArrayLike + - Array of nans, of shape (len(y), 1) if ``cv`` is ``"prefit"`` + (defined but not used) + - Dummy array of folds containing each training sample, otherwise. + Of shape (n_samples_train, cv.get_n_splits(X_train, y_train)). + """ + no_agg_cv_ = ["prefit", "split"] + no_agg_methods_ = ["naive", "base"] + fit_attributes = [ + "single_estimator_", + "estimators_", + "k_", + ] + + def __init__( + self, + estimator: Optional[RegressorMixin], + method: str, + cv: Optional[Union[int, str, BaseCrossValidator]], + agg_function: Optional[str], + n_jobs: Optional[int], + random_state: Optional[Union[int, np.random.RandomState]], + test_size: Optional[Union[int, float]], + verbose: int + ): + self.estimator = estimator + self.method = method + self.cv = cv + self.agg_function = agg_function + self.n_jobs = n_jobs + self.random_state = random_state + self.test_size = test_size + self.verbose = verbose + + @staticmethod + def _fit_oof_estimator( + estimator: RegressorMixin, + X: ArrayLike, + y: ArrayLike, + train_index: ArrayLike, + sample_weight: Optional[ArrayLike] = None, + ) -> RegressorMixin: + """ + Fit a single out-of-fold model on a given training set. + + Parameters + ---------- + estimator: RegressorMixin + Estimator to train. + + X: ArrayLike of shape (n_samples, n_features) + Input data. + + y: ArrayLike of shape (n_samples,) + Input labels. + + train_index: ArrayLike of shape (n_samples_train) + Training data indices. + + sample_weight: Optional[ArrayLike] of shape (n_samples,) + Sample weights. If None, then samples are equally weighted. + By default ``None``. + + Returns + ------- + RegressorMixin + Fitted estimator. + """ + X_train = _safe_indexing(X, train_index) + y_train = _safe_indexing(y, train_index) + if not (sample_weight is None): + sample_weight = _safe_indexing(sample_weight, train_index) + sample_weight = cast(NDArray, sample_weight) + + estimator = fit_estimator( + estimator, X_train, y_train, sample_weight=sample_weight + ) + return estimator + + @staticmethod + def _predict_oof_estimator( + estimator: RegressorMixin, + X: ArrayLike, + val_index: ArrayLike, + ) -> Tuple[NDArray, ArrayLike]: + """ + Perform predictions on a single out-of-fold model on a validation set. + + Parameters + ---------- + estimator: RegressorMixin + Estimator to train. + + X: ArrayLike of shape (n_samples, n_features) + Input data. + + val_index: ArrayLike of shape (n_samples_val) + Validation data indices. + + Returns + ------- + Tuple[NDArray, ArrayLike] + Predictions of estimator from val_index of X. + """ + X_val = _safe_indexing(X, val_index) + if _num_samples(X_val) > 0: + y_pred = estimator.predict(X_val) + else: + y_pred = np.array([]) + return y_pred, val_index + + def _aggregate_with_mask( + self, + x: NDArray, + k: NDArray + ) -> NDArray: + """ + Take the array of predictions, made by the refitted estimators, + on the testing set, and the 1-or-nan array indicating for each training + sample which one to integrate, and aggregate to produce phi-{t}(x_t) + for each training sample x_t. + + Parameters + ---------- + x: ArrayLike of shape (n_samples_test, n_estimators) + Array of predictions, made by the refitted estimators, + for each sample of the testing set. + + k: ArrayLike of shape (n_samples_training, n_estimators) + 1-or-nan array: indicates whether to integrate the prediction + of a given estimator into the aggregation, for each training + sample. + + Returns + ------- + ArrayLike of shape (n_samples_test,) + Array of aggregated predictions for each testing sample. + """ + if self.method in self.no_agg_methods_ or self.cv in self.no_agg_cv_: + raise ValueError( + "There should not be aggregation of predictions " + f"if cv is in '{self.no_agg_cv_}' " + f"or if method is in '{self.no_agg_methods_}'." + ) + elif self.agg_function == "median": + return phi2D(A=x, B=k, fun=lambda x: np.nanmedian(x, axis=1)) + # To aggregate with mean() the aggregation coud be done + # with phi2D(A=x, B=k, fun=lambda x: np.nanmean(x, axis=1). + # However, phi2D contains a np.apply_along_axis loop which + # is much slower than the matrices multiplication that can + # be used to compute the means. + elif self.agg_function in ["mean", None]: + K = np.nan_to_num(k, nan=0.0) + return np.matmul(x, (K / (K.sum(axis=1, keepdims=True))).T) + else: + raise ValueError("The value of self.agg_function is not correct") + + def _pred_multi(self, X: ArrayLike) -> NDArray: + """ + Return a prediction per train sample for each test sample, by + aggregation with matrix ``k_``. + + Parameters + ---------- + X: ArrayLike of shape (n_samples_test, n_features) + Input data + + Returns + ------- + NDArray of shape (n_samples_test, n_samples_train) + """ + y_pred_multi = np.column_stack( + [e.predict(X) for e in self.estimators_] + ) + # At this point, y_pred_multi is of shape + # (n_samples_test, n_estimators_). The method + # ``_aggregate_with_mask`` fits it to the right size + # thanks to the shape of k_. + y_pred_multi = self._aggregate_with_mask(y_pred_multi, self.k_) + return y_pred_multi + + def predict_calib(self, X: ArrayLike) -> NDArray: + """ + Perform predictions on X : the calibration set. + + Parameters + ---------- + X: ArrayLike of shape (n_samples_test, n_features) + Input data + + Returns + ------- + NDArray of shape (n_samples_test, 1) + The predictions. + """ + check_is_fitted(self, self.fit_attributes) + + if self.cv == "prefit": + y_pred = self.single_estimator_.predict(X) + else: + if self.method == "naive": + y_pred = self.single_estimator_.predict(X) + else: + cv = cast(BaseCrossValidator, self.cv) + outputs = Parallel(n_jobs=self.n_jobs, verbose=self.verbose)( + delayed(self._predict_oof_estimator)( + estimator, X, calib_index, + ) + for (_, calib_index), estimator in zip(cv.split(X), + self.estimators_) + ) + predictions, indices = map( + list, zip(*outputs) + ) + n_samples = _num_samples(X) + pred_matrix = np.full( + shape=(n_samples, cv.get_n_splits(X)), + fill_value=np.nan, + dtype=float, + ) + for i, ind in enumerate(indices): + pred_matrix[ind, i] = np.array( + predictions[i], dtype=float + ) + self.k_[ind, i] = 1 + check_nan_in_aposteriori_prediction(pred_matrix) + + y_pred = aggregate_all(self.agg_function, pred_matrix) + + return y_pred + + def fit( + self, + X: ArrayLike, + y: ArrayLike, + sample_weight: Optional[ArrayLike] = None, + ) -> EnsembleRegressor: + """ + Fit the base estimator under the ``single_estimator_`` attribute. + Fit all cross-validated estimator clones + and rearrange them into a list, the ``estimators_`` attribute. + Out-of-fold conformity scores are stored under + the ``conformity_scores_`` attribute. + + Parameters + ---------- + X: ArrayLike of shape (n_samples, n_features) + Input data. + + y: ArrayLike of shape (n_samples,) + Input labels. + + sample_weight: Optional[ArrayLike] of shape (n_samples,) + Sample weights. If None, then samples are equally weighted. + By default ``None``. + + Returns + ------- + EnsembleRegressor + The estimator fitted. + """ + # Initialization + single_estimator_: RegressorMixin + estimators_: List[RegressorMixin] = [] + full_indexes = np.arange(_num_samples(X)) + cv = self.cv + estimator = self.estimator + n_samples = _num_samples(y) + + # Computation + if cv == "prefit": + single_estimator_ = estimator + self.k_ = np.full( + shape=(n_samples, 1), fill_value=np.nan, dtype=float + ) + else: + single_estimator_ = self._fit_oof_estimator( + clone(estimator), X, y, full_indexes, sample_weight + ) + cv = cast(BaseCrossValidator, cv) + self.k_ = np.full( + shape=(n_samples, cv.get_n_splits(X, y)), + fill_value=np.nan, + dtype=float, + ) + if self.method == "naive": + estimators_ = [single_estimator_] + else: + estimators_ = Parallel(self.n_jobs, verbose=self.verbose)( + delayed(self._fit_oof_estimator)( + clone(estimator), X, y, train_index, sample_weight + ) + for train_index, _ in cv.split(X) + ) + if isinstance(cv, ShuffleSplit): + single_estimator_ = estimators_[0] + + self.single_estimator_ = single_estimator_ + self.estimators_ = estimators_ + + return self + + def predict( + self, + X: ArrayLike, + ensemble: bool = False, + return_multi_pred: bool = True + ) -> Union[NDArray, Tuple[NDArray, NDArray, NDArray]]: + """ + Predict target from X. It also computes the prediction per train sample + for each test sample according to ``self.method``. + + Parameters + ---------- + X: ArrayLike of shape (n_samples, n_features) + Test data. + + ensemble: bool + Boolean determining whether the predictions are ensembled or not. + If ``False``, predictions are those of the model trained on the + whole training set. + If ``True``, predictions from perturbed models are aggregated by + the aggregation function specified in the ``agg_function`` + attribute. + + If ``cv`` is ``"prefit"`` or ``"split"``, ``ensemble`` is ignored. + + By default ``False``. + + return_multi_pred: bool + If ``True`` the method returns the predictions and the multiple + predictions (3 arrays). If ``False`` the method return the + simple predictions only. + + Returns + ------- + Tuple[NDArray, NDArray, NDArray] + - Predictions + - The multiple predictions for the lower bound of the intervals. + - The multiple predictions for the upper bound of the intervals. + """ + check_is_fitted(self, self.fit_attributes) + + y_pred = self.single_estimator_.predict(X) + if not return_multi_pred and not ensemble: + return y_pred + + if self.method in self.no_agg_methods_ or self.cv in self.no_agg_cv_: + y_pred_multi_low = y_pred[:, np.newaxis] + y_pred_multi_up = y_pred[:, np.newaxis] + else: + y_pred_multi = self._pred_multi(X) + + if self.method == "minmax": + y_pred_multi_low = np.min(y_pred_multi, axis=1, keepdims=True) + y_pred_multi_up = np.max(y_pred_multi, axis=1, keepdims=True) + else: + y_pred_multi_low = y_pred_multi + y_pred_multi_up = y_pred_multi + + if ensemble: + y_pred = aggregate_all(self.agg_function, y_pred_multi) + + if return_multi_pred: + return y_pred, y_pred_multi_low, y_pred_multi_up + else: + return y_pred diff --git a/mapie/estimator/interface.py b/mapie/estimator/interface.py new file mode 100644 index 000000000..468f33185 --- /dev/null +++ b/mapie/estimator/interface.py @@ -0,0 +1,86 @@ +from __future__ import annotations + +from typing import Optional, Tuple, Union + +from sklearn.base import RegressorMixin + +from mapie._typing import ArrayLike, NDArray + + +class EnsembleEstimator(RegressorMixin): + """ + This class implements methods to handle the training and usage of the + estimator. This estimator can be unique or composed by cross validated + estimators. + """ + + def fit( + self, + X: ArrayLike, + y: ArrayLike, + sample_weight: Optional[ArrayLike] = None, + ) -> EnsembleEstimator: + """ + Fit the base estimator under the ``single_estimator_`` attribute. + Fit all cross-validated estimator clones + and rearrange them into a list, the ``estimators_`` attribute. + Out-of-fold conformity scores are stored under + the ``conformity_scores_`` attribute. + + Parameters + ---------- + X: ArrayLike of shape (n_samples, n_features) + Input data. + + y: ArrayLike of shape (n_samples,) + Input labels. + + sample_weight: Optional[ArrayLike] of shape (n_samples,) + Sample weights. If None, then samples are equally weighted. + By default ``None``. + + Returns + ------- + EnsembleRegressor + The estimator fitted. + """ + + def predict( + self, + X: ArrayLike, + ensemble: bool = False, + return_multi_pred: bool = True + ) -> Union[NDArray, Tuple[NDArray, NDArray, NDArray]]: + """ + Predict target from X. It also computes the prediction per train sample + for each test sample according to ``self.method``. + + Parameters + ---------- + X: ArrayLike of shape (n_samples, n_features) + Test data. + + ensemble: bool + Boolean determining whether the predictions are ensembled or not. + If ``False``, predictions are those of the model trained on the + whole training set. + If ``True``, predictions from perturbed models are aggregated by + the aggregation function specified in the ``agg_function`` + attribute. + + If ``cv`` is ``"prefit"`` or ``"split"``, ``ensemble`` is ignored. + + By default ``False``. + + return_multi_pred: bool + If ``True`` the method returns the predictions and the multiple + predictions (3 arrays). If ``False`` the method return the + simple predictions only. + + Returns + ------- + Tuple[NDArray, NDArray, NDArray] + - Predictions + - The multiple predictions for the lower bound of the intervals. + - The multiple predictions for the upper bound of the intervals. + """ diff --git a/mapie/regression/regression.py b/mapie/regression/regression.py index a7f17316f..e51fbecce 100644 --- a/mapie/regression/regression.py +++ b/mapie/regression/regression.py @@ -1,26 +1,23 @@ from __future__ import annotations -from typing import Iterable, List, Optional, Tuple, Union, cast +from typing import Iterable, Optional, Tuple, Union, cast import numpy as np -from joblib import Parallel, delayed -from sklearn.base import BaseEstimator, RegressorMixin, clone +from sklearn.base import BaseEstimator, RegressorMixin from sklearn.linear_model import LinearRegression -from sklearn.model_selection import BaseCrossValidator, ShuffleSplit +from sklearn.model_selection import BaseCrossValidator from sklearn.pipeline import Pipeline -from sklearn.utils import _safe_indexing, check_random_state -from sklearn.utils.validation import (_check_y, _num_samples, check_is_fitted, +from sklearn.utils import check_random_state +from sklearn.utils.validation import (_check_y, check_is_fitted, indexable) -from mapie._compatibility import np_nanquantile from mapie._typing import ArrayLike, NDArray -from mapie.aggregation_functions import aggregate_all, phi2D from mapie.conformity_scores import ConformityScore +from mapie.estimator.estimator import EnsembleRegressor from mapie.utils import (check_alpha, check_alpha_and_n_samples, check_conformity_score, check_cv, check_estimator_fit_predict, check_n_features_in, - check_n_jobs, check_nan_in_aposteriori_prediction, - check_null_weight, check_verbose, fit_estimator) + check_n_jobs, check_null_weight, check_verbose) class MapieRegressor(BaseEstimator, RegressorMixin): @@ -162,21 +159,12 @@ class MapieRegressor(BaseEstimator, RegressorMixin): valid_methods_: List[str] List of all valid methods. - single_estimator_: sklearn.RegressorMixin - Estimator fitted on the whole training set. - - estimators_: list - List of out-of-folds estimators. + estimator_: EnsembleRegressor + Sklearn estimator that handle all that is related to the estimator. conformity_scores_: ArrayLike of shape (n_samples_train,) Conformity scores between ``y_train`` and ``y_pred``. - k_: ArrayLike - - Array of nans, of shape (len(y), 1) if ``cv`` is ``"prefit"`` - (defined but not used) - - Dummy array of folds containing each training sample, otherwise. - Of shape (n_samples_train, cv.get_n_splits(X_train, y_train)). - n_features_in_: int Number of features passed to the ``fit`` method. @@ -220,9 +208,7 @@ class MapieRegressor(BaseEstimator, RegressorMixin): valid_agg_functions_ = [None, "median", "mean"] ensemble_agg_functions_ = ["median", "mean"] fit_attributes = [ - "single_estimator_", - "estimators_", - "k_", + "estimator_", "conformity_scores_", "conformity_score_function_", "n_features_in_", @@ -396,139 +382,37 @@ def _check_ensemble( f"in '{self.ensemble_agg_functions_}'." ) - def _fit_and_predict_oof_model( + def _check_fit_parameters( self, - estimator: RegressorMixin, X: ArrayLike, y: ArrayLike, - train_index: ArrayLike, - val_index: ArrayLike, sample_weight: Optional[ArrayLike] = None, - ) -> Tuple[RegressorMixin, NDArray, ArrayLike]: - """ - Fit a single out-of-fold model on a given training set and - perform predictions on a test set. - - Parameters - ---------- - estimator: RegressorMixin - Estimator to train. - - X: ArrayLike of shape (n_samples, n_features) - Input data. - - y: ArrayLike of shape (n_samples,) - Input labels. - - train_index: ArrayLike of shape (n_samples_train) - Training data indices. - - val_index: ArrayLike of shape (n_samples_val) - Validation data indices. - - sample_weight: Optional[ArrayLike] of shape (n_samples,) - Sample weights. If ``None``, then samples are equally weighted. - By default ``None``. - - Returns - ------- - Tuple[RegressorMixin, NDArray, ArrayLike] - - - [0]: RegressorMixin, fitted estimator - - [1]: NDArray of shape (n_samples_val,), - estimator predictions on the validation fold. - - [2]: ArrayLike of shape (n_samples_val,), - validation data indices. - """ - X_train = _safe_indexing(X, train_index) - y_train = _safe_indexing(y, train_index) - X_val = _safe_indexing(X, val_index) - if sample_weight is None: - estimator = fit_estimator(estimator, X_train, y_train) - else: - sample_weight_train = _safe_indexing(sample_weight, train_index) - estimator = fit_estimator( - estimator, X_train, y_train, sample_weight_train - ) - if _num_samples(X_val) > 0: - y_pred = estimator.predict(X_val) - else: - y_pred = np.array([]) - return estimator, y_pred, val_index - - def _aggregate_with_mask( - self, - x: NDArray, - k: NDArray - ) -> NDArray: - """ - Take the array of predictions, made by the refitted estimators, - on the testing set, and the 1-or-nan array indicating for each training - sample which one to integrate, and aggregate to produce phi-{t}(x_t) - for each training sample x_t. - - Parameters: - ----------- - x: ArrayLike of shape (n_samples_test, n_estimators) - Array of predictions, made by the refitted estimators, - for each sample of the testing set. - - k: ArrayLike of shape (n_samples_training, n_estimators) - 1-or-nan array: indicates whether to integrate the prediction - of a given estimator into the aggregation, for each training - sample. - - Returns: - -------- - ArrayLike of shape (n_samples_test,) - Array of aggregated predictions for each testing sample. - """ - if self.method in self.no_agg_methods_ \ - or self.cv in self.no_agg_cv_: - raise ValueError( - "There should not be aggregation of predictions " - f"if cv is in '{self.no_agg_cv_}' " - f"or if method is in '{self.no_agg_methods_}'." - ) - elif self.agg_function == "median": - return phi2D(A=x, B=k, fun=lambda x: np.nanmedian(x, axis=1)) - # To aggregate with mean() the aggregation coud be done - # with phi2D(A=x, B=k, fun=lambda x: np.nanmean(x, axis=1). - # However, phi2D contains a np.apply_along_axis loop which - # is much slower than the matrices multiplication that can - # be used to compute the means. - elif self.agg_function in ["mean", None]: - K = np.nan_to_num(k, nan=0.0) - return np.matmul(x, (K / (K.sum(axis=1, keepdims=True))).T) - else: - raise ValueError("The value of self.agg_function is not correct") - - def _pred_multi( - self, - X: ArrayLike - ) -> NDArray: - """ - Return a prediction per train sample for each test sample, by - aggregation with matrix ``k_``. - - Parameters - ---------- - X: NDArray of shape (n_samples_test, n_features) - Input data - - Returns - ------- - NDArray of shape (n_samples_test, n_samples_train) - """ - y_pred_multi = np.column_stack( - [e.predict(X) for e in self.estimators_] + ): + # Checking + self._check_parameters() + cv = check_cv( + self.cv, test_size=self.test_size, random_state=self.random_state + ) + estimator = self._check_estimator(self.estimator) + agg_function = self._check_agg_function(self.agg_function) + cs_estimator = check_conformity_score( + self.conformity_score ) - # At this point, y_pred_multi is of shape - # (n_samples_test, n_estimators_). The method - # ``_aggregate_with_mask`` fits it to the right size - # thanks to the shape of k_. - y_pred_multi = self._aggregate_with_mask(y_pred_multi, self.k_) - return y_pred_multi + X, y = indexable(X, y) + y = _check_y(y) + sample_weight, X, y = check_null_weight(sample_weight, X, y) + self.n_features_in_ = check_n_features_in(X) + + # Casting + cv = cast(BaseCrossValidator, cv) + estimator = cast(RegressorMixin, estimator) + cs_estimator = cast(ConformityScore, cs_estimator) + agg_function = cast(Optional[str], agg_function) + X = cast(NDArray, X) + y = cast(NDArray, y) + sample_weight = cast(Optional[NDArray], sample_weight) + + return estimator, cs_estimator, agg_function, cv, X, y, sample_weight def fit( self, @@ -539,11 +423,9 @@ def fit( """ Fit estimator and compute conformity scores used for prediction intervals. - Fit the base estimator under the ``single_estimator_`` attribute. - Fit all cross-validated estimator clones - and rearrange them into a list, the ``estimators_`` attribute. - Out-of-fold conformity scores are stored under - the ``conformity_scores_`` attribute. + + All the types of estimator (single or cross validated ones) are + encapsulated under EnsembleRegressor. Parameters ---------- @@ -570,84 +452,34 @@ def fit( The model itself. """ # Checks - self._check_parameters() - cv = check_cv( - self.cv, test_size=self.test_size, random_state=self.random_state - ) - estimator = self._check_estimator(self.estimator) - agg_function = self._check_agg_function(self.agg_function) - X, y = indexable(X, y) - y = _check_y(y) - sample_weight = cast(Optional[NDArray], sample_weight) - self.n_features_in_ = check_n_features_in(X, cv, estimator) - sample_weight, X, y = check_null_weight(sample_weight, X, y) - self.conformity_score_function_ = check_conformity_score( - self.conformity_score + (estimator, + self.conformity_score_function_, + agg_function, + cv, + X, + y, + sample_weight) = self._check_fit_parameters(X, y, sample_weight) + + self.estimator_ = EnsembleRegressor( + estimator, + self.method, + cv, + agg_function, + self.n_jobs, + self.random_state, + self.test_size, + self.verbose ) - y = cast(NDArray, y) - n_samples = _num_samples(y) - - # Initialization - self.estimators_: List[RegressorMixin] = [] - - # Work - if cv == "prefit": - self.single_estimator_ = estimator - y_pred = self.single_estimator_.predict(X) - self.k_ = np.full( - shape=(n_samples, 1), fill_value=np.nan, dtype=float - ) - else: - cv = cast(BaseCrossValidator, cv) - self.k_ = np.full( - shape=(n_samples, cv.get_n_splits(X, y)), - fill_value=np.nan, - dtype=float, - ) - - self.single_estimator_ = fit_estimator( - clone(estimator), X, y, sample_weight + # Fit the prediction function + self.estimator_ = self.estimator_.fit(X, y, sample_weight) + y_pred = self.estimator_.predict_calib(X) + + # Compute the conformity scores (manage jk-ab case) + self.conformity_scores_ = \ + self.conformity_score_function_.get_conformity_scores( + X, y, y_pred ) - if self.method == "naive": - y_pred = self.single_estimator_.predict(X) - else: - outputs = Parallel(n_jobs=self.n_jobs, verbose=self.verbose)( - delayed(self._fit_and_predict_oof_model)( - clone(estimator), - X, - y, - train_index, - val_index, - sample_weight, - ) - for train_index, val_index in cv.split(X) - ) - self.estimators_, predictions, val_indices = map( - list, zip(*outputs) - ) - - pred_matrix = np.full( - shape=(n_samples, cv.get_n_splits(X, y)), - fill_value=np.nan, - dtype=float, - ) - for i, val_ind in enumerate(val_indices): - pred_matrix[val_ind, i] = np.array( - predictions[i], dtype=float - ) - self.k_[val_ind, i] = 1 - check_nan_in_aposteriori_prediction(pred_matrix) - - y_pred = aggregate_all(agg_function, pred_matrix) - - self.conformity_scores_ = ( - self.conformity_score_function_.get_conformity_scores(y, y_pred) - ) - - if isinstance(cv, ShuffleSplit): - self.single_estimator_ = self.estimators_[0] - return self def predict( @@ -708,74 +540,24 @@ def predict( self._check_ensemble(ensemble) alpha = cast(Optional[NDArray], check_alpha(alpha)) - y_pred = self.single_estimator_.predict(X) - n = len(self.conformity_scores_) - if alpha is None: + y_pred = self.estimator_.predict( + X, ensemble, return_multi_pred=False + ) return np.array(y_pred) - alpha_np = cast(NDArray, alpha) - check_alpha_and_n_samples(alpha_np, n) - - if self.method in self.no_agg_methods_ \ - or self.cv in self.no_agg_cv_: - y_pred_multi_low = y_pred[:, np.newaxis] - y_pred_multi_up = y_pred[:, np.newaxis] else: - y_pred_multi = self._pred_multi(X) - - if self.method == "minmax": - y_pred_multi_low = np.min(y_pred_multi, axis=1, keepdims=True) - y_pred_multi_up = np.max(y_pred_multi, axis=1, keepdims=True) - else: - y_pred_multi_low = y_pred_multi - y_pred_multi_up = y_pred_multi - - if ensemble: - y_pred = aggregate_all(self.agg_function, y_pred_multi) - - # compute distributions of lower and upper bounds - if self.conformity_score_function_.sym: - conformity_scores_low = -self.conformity_scores_ - conformity_scores_up = self.conformity_scores_ - else: - conformity_scores_low = self.conformity_scores_ - conformity_scores_up = self.conformity_scores_ - alpha_np = alpha_np / 2 - - lower_bounds = ( - self.conformity_score_function_.get_estimation_distribution( - y_pred_multi_low, conformity_scores_low - ) - ) - upper_bounds = ( - self.conformity_score_function_.get_estimation_distribution( - y_pred_multi_up, conformity_scores_up - ) - ) - - # get desired confidence intervals according to alpha - y_pred_low = np.column_stack( - [ - np_nanquantile( - lower_bounds.astype(float), - _alpha, - axis=1, - method="lower", - ) - for _alpha in alpha_np - ] - ) - y_pred_up = np.column_stack( - [ - np_nanquantile( - upper_bounds.astype(float), - 1 - _alpha, - axis=1, - method="higher", + n = len(self.conformity_scores_) + alpha_np = cast(NDArray, alpha) + check_alpha_and_n_samples(alpha_np, n) + + y_pred, y_pred_low, y_pred_up = \ + self.conformity_score_function_.get_bounds( + X, + self.estimator_, + self.conformity_scores_, + alpha_np, + ensemble, + self.method ) - for _alpha in alpha_np - ] - ) - - return y_pred, np.stack([y_pred_low, y_pred_up], axis=1) + return np.array(y_pred), np.stack([y_pred_low, y_pred_up], axis=1) diff --git a/mapie/regression/time_series_regression.py b/mapie/regression/time_series_regression.py index 080eb3e44..9e7354b98 100644 --- a/mapie/regression/time_series_regression.py +++ b/mapie/regression/time_series_regression.py @@ -283,7 +283,7 @@ def predict( check_is_fitted(self, self.fit_attributes) self._check_ensemble(ensemble) alpha = cast(Optional[NDArray], check_alpha(alpha)) - y_pred = self.single_estimator_.predict(X) + y_pred = self.estimator_.single_estimator_.predict(X) n = len(self.conformity_scores_) if alpha is None: @@ -321,7 +321,7 @@ def predict( y_pred_low = y_pred[:, np.newaxis] + lower_quantiles y_pred_up = y_pred[:, np.newaxis] + higher_quantiles else: - y_pred_multi = self._pred_multi(X) + y_pred_multi = self.estimator_._pred_multi(X) pred = aggregate_all(self.agg_function, y_pred_multi) lower_bounds, upper_bounds = pred, pred diff --git a/mapie/tests/test_common.py b/mapie/tests/test_common.py index 0f1999824..45379bc24 100644 --- a/mapie/tests/test_common.py +++ b/mapie/tests/test_common.py @@ -108,7 +108,12 @@ def test_none_estimator(pack: Tuple[BaseEstimator, BaseEstimator]) -> None: MapieEstimator, DefaultEstimator = pack mapie_estimator = MapieEstimator(estimator=None) mapie_estimator.fit(X_toy, y_toy) - assert isinstance(mapie_estimator.single_estimator_, DefaultEstimator) + if isinstance(mapie_estimator, MapieClassifier): + assert isinstance(mapie_estimator.single_estimator_, DefaultEstimator) + if isinstance(mapie_estimator, MapieRegressor): + assert isinstance( + mapie_estimator.estimator_.single_estimator_, DefaultEstimator + ) @pytest.mark.parametrize("estimator", [0, "a", KFold(), ["a", "b"]]) diff --git a/mapie/tests/test_conformity_scores.py b/mapie/tests/test_conformity_scores.py index 7441302fb..a47cfbf40 100644 --- a/mapie/tests/test_conformity_scores.py +++ b/mapie/tests/test_conformity_scores.py @@ -10,6 +10,7 @@ y_pred_list = [4, 7, 10, 12, 13, 12] conf_scores_list = [1, 0, -1, -1, 0, 3] conf_scores_gamma_list = [1 / 4, 0, -1 / 10, -1 / 12, 0, 3 / 12] +random_state = 42 class DummyConformityScore(ConformityScore): @@ -17,12 +18,12 @@ def __init__(self) -> None: super().__init__(sym=True, consistency_check=True) def get_signed_conformity_scores( - self, y: ArrayLike, y_pred: ArrayLike, + self, X: ArrayLike, y: ArrayLike, y_pred: ArrayLike, ) -> NDArray: return np.subtract(y, y_pred) def get_estimation_distribution( - self, y_pred: ArrayLike, conformity_scores: ArrayLike + self, X: ArrayLike, y_pred: ArrayLike, conformity_scores: ArrayLike ) -> NDArray: """ A positive constant is added to the sum between predictions and @@ -45,10 +46,12 @@ def test_absolute_conformity_score_get_conformity_scores( """Test conformity score computation for AbsoluteConformityScore.""" abs_conf_score = AbsoluteConformityScore() signed_conf_scores = abs_conf_score.get_signed_conformity_scores( - y_toy, y_pred + X_toy, y_toy, y_pred ) - conf_scores = abs_conf_score.get_conformity_scores(y_toy, y_pred) - expected_signed_conf_scores = np.array([1, 0, -1, -1, 0, 3]) + conf_scores = abs_conf_score.get_conformity_scores( + X_toy, y_toy, y_pred + ) + expected_signed_conf_scores = np.array(conf_scores_list) expected_conf_scores = np.abs(expected_signed_conf_scores) np.testing.assert_allclose(signed_conf_scores, expected_signed_conf_scores) np.testing.assert_allclose(conf_scores, expected_conf_scores) @@ -63,7 +66,9 @@ def test_absolute_conformity_score_get_estimation_distribution( ) -> None: """Test conformity observed value computation for AbsoluteConformityScore.""" # noqa: E501 abs_conf_score = AbsoluteConformityScore() - y_obs = abs_conf_score.get_estimation_distribution(y_pred, conf_scores) + y_obs = abs_conf_score.get_estimation_distribution( + X_toy, y_pred, conf_scores + ) np.testing.assert_allclose(y_obs, y_toy) @@ -72,10 +77,10 @@ def test_absolute_conformity_score_consistency(y_pred: NDArray) -> None: """Test methods consistency for AbsoluteConformityScore.""" abs_conf_score = AbsoluteConformityScore() signed_conf_scores = abs_conf_score.get_signed_conformity_scores( - y_toy, y_pred + X_toy, y_toy, y_pred ) y_obs = abs_conf_score.get_estimation_distribution( - y_pred, signed_conf_scores + X_toy, y_pred, signed_conf_scores ) np.testing.assert_allclose(y_obs, y_toy) @@ -86,7 +91,9 @@ def test_gamma_conformity_score_get_conformity_scores( ) -> None: """Test conformity score computation for GammaConformityScore.""" gamma_conf_score = GammaConformityScore() - conf_scores = gamma_conf_score.get_conformity_scores(y_toy, y_pred) + conf_scores = gamma_conf_score.get_conformity_scores( + X_toy, y_toy, y_pred + ) expected_signed_conf_scores = np.array(conf_scores_gamma_list) np.testing.assert_allclose(conf_scores, expected_signed_conf_scores) @@ -104,7 +111,9 @@ def test_gamma_conformity_score_get_estimation_distribution( ) -> None: """Test conformity observed value computation for GammaConformityScore.""" # noqa: E501 gamma_conf_score = GammaConformityScore() - y_obs = gamma_conf_score.get_estimation_distribution(y_pred, conf_scores) + y_obs = gamma_conf_score.get_estimation_distribution( + X_toy, y_pred, conf_scores + ) np.testing.assert_allclose(y_obs, y_toy) @@ -113,10 +122,10 @@ def test_gamma_conformity_score_consistency(y_pred: NDArray) -> None: """Test methods consistency for GammaConformityScore.""" gamma_conf_score = GammaConformityScore() signed_conf_scores = gamma_conf_score.get_signed_conformity_scores( - y_toy, y_pred + X_toy, y_toy, y_pred ) y_obs = gamma_conf_score.get_estimation_distribution( - y_pred, signed_conf_scores + X_toy, y_pred, signed_conf_scores ) np.testing.assert_allclose(y_obs, y_toy) @@ -136,7 +145,9 @@ def test_gamma_conformity_score_check_oberved_value( """Test methods consistency for GammaConformityScore.""" gamma_conf_score = GammaConformityScore() with pytest.raises(ValueError): - gamma_conf_score.get_signed_conformity_scores(y_toy, y_pred) + gamma_conf_score.get_signed_conformity_scores( + [], y_toy, y_pred + ) @pytest.mark.parametrize( @@ -147,6 +158,14 @@ def test_gamma_conformity_score_check_oberved_value( [1, -7, 10, 12, 13, 12], ], ) +@pytest.mark.parametrize( + "X_toy", + [ + np.array([0, -7, 10, 12, 0, 12]).reshape(-1, 1), + np.array([0, 7, -10, 12, 1, -12]).reshape(-1, 1), + np.array([12, -7, 0, 12, 13, 2]).reshape(-1, 1), + ], +) @pytest.mark.parametrize( "conf_scores", [ @@ -155,7 +174,7 @@ def test_gamma_conformity_score_check_oberved_value( ], ) def test_gamma_conformity_score_check_predicted_value( - y_pred: NDArray, conf_scores: NDArray + y_pred: NDArray, conf_scores: NDArray, X_toy: NDArray ) -> None: """Test methods consistency for GammaConformityScore.""" gamma_conf_score = GammaConformityScore() @@ -163,27 +182,31 @@ def test_gamma_conformity_score_check_predicted_value( ValueError, match=r".*At least one of the predicted target is negative.*" ): - gamma_conf_score.get_signed_conformity_scores(y_toy, y_pred) + gamma_conf_score.get_signed_conformity_scores( + X_toy, y_toy, y_pred + ) with pytest.raises( ValueError, match=r".*At least one of the predicted target is negative.*" ): - gamma_conf_score.get_estimation_distribution(y_pred, conf_scores) + gamma_conf_score.get_estimation_distribution( + X_toy, y_pred, conf_scores + ) def test_check_consistency() -> None: """ - Test that a dummy ConformityScore class that gives inconsistent conformity - scores and distributions raises an error. + Test that a dummy ConformityScore class that gives inconsistent scores + and distributions raises an error. """ dummy_conf_score = DummyConformityScore() conformity_scores = dummy_conf_score.get_signed_conformity_scores( - y_toy, y_pred_list + X_toy, y_toy, y_pred_list ) with pytest.raises( ValueError, match=r".*The two functions get_conformity_scores.*" ): dummy_conf_score.check_consistency( - y_toy, y_pred_list, conformity_scores + X_toy, y_toy, y_pred_list, conformity_scores ) diff --git a/mapie/tests/test_regression.py b/mapie/tests/test_regression.py index c089d0e24..4a9d9975d 100644 --- a/mapie/tests/test_regression.py +++ b/mapie/tests/test_regression.py @@ -18,12 +18,13 @@ from sklearn.utils.validation import check_is_fitted from typing_extensions import TypedDict -from mapie._typing import ArrayLike, NDArray +from mapie._typing import NDArray from mapie.aggregation_functions import aggregate_all from mapie.conformity_scores import (AbsoluteConformityScore, ConformityScore, GammaConformityScore) from mapie.metrics import regression_coverage_score from mapie.regression import MapieRegressor +from mapie.estimator.estimator import EnsembleRegressor from mapie.subsample import Subsample X_toy = np.array([0, 1, 2, 3, 4, 5]).reshape(-1, 1) @@ -173,8 +174,8 @@ def test_valid_estimator(strategy: str) -> None: estimator=DummyRegressor(), **STRATEGIES[strategy] ) mapie_reg.fit(X_toy, y_toy) - assert isinstance(mapie_reg.single_estimator_, DummyRegressor) - for estimator in mapie_reg.estimators_: + assert isinstance(mapie_reg.estimator_.single_estimator_, DummyRegressor) + for estimator in mapie_reg.estimator_.estimators_: assert isinstance(estimator, DummyRegressor) @@ -502,30 +503,42 @@ def test_aggregate_with_mask_with_prefit() -> None: """ Test ``_aggregate_with_mask`` in case ``cv`` is ``"prefit"``. """ - mapie_reg = MapieRegressor(cv="prefit") + mapie_reg = MapieRegressor(LinearRegression().fit(X, y), cv="prefit") + mapie_reg = mapie_reg.fit(X, y) with pytest.raises( ValueError, match=r".*There should not be aggregation of predictions if cv is*", ): - mapie_reg._aggregate_with_mask(k, k) - - mapie_reg = MapieRegressor(agg_function="nonsense") + mapie_reg.estimator_._aggregate_with_mask(k, k) + + +def test_aggregate_with_mask_with_invalid_agg_function() -> None: + """Test ``_aggregate_with_mask`` in case ``agg_function`` is invalid.""" + ens_reg = EnsembleRegressor( + LinearRegression(), + "plus", + KFold(n_splits=5, random_state=None, shuffle=True), + "nonsense", + None, + random_state, + 0.20, + False + ) with pytest.raises( ValueError, match=r".*The value of self.agg_function is not correct*", ): - mapie_reg._aggregate_with_mask(k, k) + ens_reg._aggregate_with_mask(k, k) def test_pred_loof_isnan() -> None: """Test that if validation set is empty then prediction is empty.""" mapie_reg = MapieRegressor() - y_pred: ArrayLike - _, y_pred, _ = mapie_reg._fit_and_predict_oof_model( + y_pred: NDArray + mapie_reg = mapie_reg.fit(X, y) + y_pred, _ = mapie_reg.estimator_._predict_oof_estimator( estimator=LinearRegression(), X=X_toy, - y=y_toy, - train_index=[0, 1, 2, 3, 4], val_index=[], ) assert len(y_pred) == 0 @@ -575,3 +588,27 @@ def test_conformity_score( ) mapie_reg.fit(X, y + 1e3) mapie_reg.predict(X, alpha=0.05) + + +@pytest.mark.parametrize("ensemble", [True, False]) +def test_return_only_ypred(ensemble: bool) -> None: + """Test that if return_multi_pred is False it only returns y_pred.""" + mapie_reg = MapieRegressor() + mapie_reg.fit(X_toy, y_toy) + output = mapie_reg.estimator_.predict( + X_toy, ensemble=ensemble, return_multi_pred=False + ) + assert len(output) == len(X_toy) + + +@pytest.mark.parametrize("ensemble", [True, False]) +def test_return_multi_pred(ensemble: bool) -> None: + """ + Test that if return_multi_pred is True it returns y_pred and multi_pred. + """ + mapie_reg = MapieRegressor() + mapie_reg.fit(X_toy, y_toy) + output = mapie_reg.estimator_.predict( + X_toy, ensemble=ensemble, return_multi_pred=True + ) + assert len(output) == 3 diff --git a/mapie/tests/test_time_series_regression.py b/mapie/tests/test_time_series_regression.py index 1bdd5beea..110b423a2 100644 --- a/mapie/tests/test_time_series_regression.py +++ b/mapie/tests/test_time_series_regression.py @@ -317,11 +317,10 @@ def test_invalid_aggregate_all() -> None: def test_pred_loof_isnan() -> None: """Test that if validation set is empty then prediction is empty.""" mapie_ts_reg = MapieTimeSeriesRegressor() - _, y_pred, _ = mapie_ts_reg._fit_and_predict_oof_model( - estimator=mapie_ts_reg, + mapie_ts_reg.fit(X_toy, y_toy) + y_pred, _ = mapie_ts_reg.estimator_._predict_oof_estimator( + estimator=mapie_ts_reg.estimator_.estimators_[0], X=X_toy, - y=y_toy, - train_index=[0, 1, 2, 3, 4], val_index=[], ) assert len(y_pred) == 0