code: """
  Generalized Linear Models.
  """
  
  # Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
  # Fabian Pedregosa <fabian.pedregosa@inria.fr>
  # Olivier Grisel <olivier.grisel@ensta.org>
  #         Vincent Michel <vincent.michel@inria.fr>
  #         Peter Prettenhofer <peter.prettenhofer@gmail.com>
  #         Mathieu Blondel <mathieu@mblondel.org>
  #         Lars Buitinck
  #         Maryan Morel <maryan.morel@polytechnique.edu>
  #         Giorgio Patrini <giorgio.patrini@anu.edu.au>
  # License: BSD 3 clause
  
  from abc import ABCMeta, abstractmethod
  import numbers
  import warnings
  
  import numpy as np
  import scipy.sparse as sp
  from scipy import linalg
  from scipy import optimize
  from scipy import sparse
  from scipy.special import expit
  from joblib import Parallel
  
  from ..base import (BaseEstimator, ClassifierMixin, RegressorMixin,
                      MultiOutputMixin)
  from ..utils import check_array
  from ..utils.validation import FLOAT_DTYPES
  from ..utils.validation import _deprecate_positional_args
  from ..utils import check_random_state
  from ..utils.extmath import safe_sparse_dot
  from ..utils.sparsefuncs import mean_variance_axis, inplace_column_scale
  from ..utils.fixes import sparse_lsqr
  from ..utils._seq_dataset import ArrayDataset32, CSRDataset32
  from ..utils._seq_dataset import ArrayDataset64, CSRDataset64
  from ..utils.validation import check_is_fitted, _check_sample_weight
  from ..utils.fixes import delayed
  from ..preprocessing import normalize as f_normalize
  
  # TODO: bayesian_ridge_regression and bayesian_regression_ard
  # should be squashed into its respective objects.
  
  SPARSE_INTERCEPT_DECAY = 0.01
  # For sparse data intercept updates are scaled by this decay factor to avoid
  # intercept oscillation.
  
  
  def make_dataset(X, y, sample_weight, random_state=None):
      """Create ``Dataset`` abstraction for sparse and dense inputs.
  
      This also returns the ``intercept_decay`` which is different
      for sparse datasets.
  
      Parameters
      ----------
      X : array-like, shape (n_samples, n_features)
          Training data
  
      y : array-like, shape (n_samples, )
          Target values.
  
      sample_weight : numpy array of shape (n_samples,)
          The weight of each sample
  
      random_state : int, RandomState instance or None (default)
          Determines random number generation for dataset shuffling and noise.
          Pass an int for reproducible output across multiple function calls.
          See :term:`Glossary <random_state>`.
  
      Returns
      -------
      dataset
          The ``Dataset`` abstraction
      intercept_decay
          The intercept decay
      """
  
      rng = check_random_state(random_state)
      # seed should never be 0 in SequentialDataset64
      seed = rng.randint(1, np.iinfo(np.int32).max)
  
      if X.dtype == np.float32:
          CSRData = CSRDataset32
          ArrayData = ArrayDataset32
      else:
          CSRData = CSRDataset64
          ArrayData = ArrayDataset64
  
      if sp.issparse(X):
          dataset = CSRData(X.data, X.indptr, X.indices, y, sample_weight,
                            seed=seed)
          intercept_decay = SPARSE_INTERCEPT_DECAY
      else:
          X = np.ascontiguousarray(X)
          dataset = ArrayData(X, y, sample_weight, seed=seed)
          intercept_decay = 1.0
  
      return dataset, intercept_decay
  
  
  def _preprocess_data(X, y, fit_intercept, normalize=False, copy=True,
                       sample_weight=None, return_mean=False, check_input=True):
      """Center and scale data.
  
      Centers data to have mean zero along axis 0. If fit_intercept=False or if
      the X is a sparse matrix, no centering is done, but normalization can still
      be applied. The function returns the statistics necessary to reconstruct
      the input data, which are X_offset, y_offset, X_scale, such that the output
  
          X = (X - X_offset) / X_scale
  
      X_scale is the L2 norm of X - X_offset. If sample_weight is not None,
      then the weighted mean of X and y is zero, and not the mean itself. If
      return_mean=True, the mean, eventually weighted, is returned, independently
      of whether X was centered (option used for optimization with sparse data in
      coordinate_descend).
  
      This is here because nearly all linear models will want their data to be
      centered. This function also systematically makes y consistent with X.dtype
      """
      if isinstance(sample_weight, numbers.Number):
          sample_weight = None
      if sample_weight is not None:
          sample_weight = np.asarray(sample_weight)
  
      if check_input:
          X = check_array(X, copy=copy, accept_sparse=['csr', 'csc'],
                          dtype=FLOAT_DTYPES)
      elif copy:
          if sp.issparse(X):
              X = X.copy()
          else:
              X = X.copy(order='K')
  
      y = np.asarray(y, dtype=X.dtype)
  
      if fit_intercept:
          if sp.issparse(X):
              X_offset, X_var = mean_variance_axis(X, axis=0)
              if not return_mean:
                  X_offset[:] = X.dtype.type(0)
  
              if normalize:
  
                  # TODO: f_normalize could be used here as well but the function
                  # inplace_csr_row_normalize_l2 must be changed such that it
                  # can return also the norms computed internally
  
                  # transform variance to norm in-place
                  X_var *= X.shape[0]
                  X_scale = np.sqrt(X_var, X_var)
                  del X_var
                  X_scale[X_scale == 0] = 1
                  inplace_column_scale(X, 1. / X_scale)
              else:
                  X_scale = np.ones(X.shape[1], dtype=X.dtype)
  
          else:
              X_offset = np.average(X, axis=0, weights=sample_weight)
              X -= X_offset
              if normalize:
                  X, X_scale = f_normalize(X, axis=0, copy=False,
                                           return_norm=True)
              else:
                  X_scale = np.ones(X.shape[1], dtype=X.dtype)
          y_offset = np.average(y, axis=0, weights=sample_weight)
          y = y - y_offset
      else:
          X_offset = np.zeros(X.shape[1], dtype=X.dtype)
          X_scale = np.ones(X.shape[1], dtype=X.dtype)
          if y.ndim == 1:
              y_offset = X.dtype.type(0)
          else:
              y_offset = np.zeros(y.shape[1], dtype=X.dtype)
  
      return X, y, X_offset, y_offset, X_scale
  
  
  # TODO: _rescale_data should be factored into _preprocess_data.
  # Currently, the fact that sag implements its own way to deal with
  # sample_weight makes the refactoring tricky.
  
  def _rescale_data(X, y, sample_weight):
      """Rescale data sample-wise by square root of sample_weight.
  
      For many linear models, this enables easy support for sample_weight.
  
      Returns
      -------
      X_rescaled : {array-like, sparse matrix}
  
      y_rescaled : {array-like, sparse matrix}
      """
      n_samples = X.shape[0]
      sample_weight = np.asarray(sample_weight)
      if sample_weight.ndim == 0:
          sample_weight = np.full(n_samples, sample_weight,
                                  dtype=sample_weight.dtype)
      sample_weight = np.sqrt(sample_weight)
      sw_matrix = sparse.dia_matrix((sample_weight, 0),
                                    shape=(n_samples, n_samples))
      X = safe_sparse_dot(sw_matrix, X)
      y = safe_sparse_dot(sw_matrix, y)
      return X, y
  
  
  class LinearModel(BaseEstimator, metaclass=ABCMeta):
      """Base class for Linear Models"""
  
      @abstractmethod
      def fit(self, X, y):
          """Fit model."""
  
      def _decision_function(self, X):
          check_is_fitted(self)
  
          X = check_array(X, accept_sparse=['csr', 'csc', 'coo'])
          return safe_sparse_dot(X, self.coef_.T,
                                 dense_output=True) + self.intercept_
  
      def predict(self, X):
          """
          Predict using the linear model.
  
          Parameters
          ----------
          X : array-like or sparse matrix, shape (n_samples, n_features)
              Samples.
  
          Returns
          -------
          C : array, shape (n_samples,)
              Returns predicted values.
          """
          return self._decision_function(X)
  
      _preprocess_data = staticmethod(_preprocess_data)
  
      def _set_intercept(self, X_offset, y_offset, X_scale):
          """Set the intercept_
          """
          if self.fit_intercept:
              self.coef_ = self.coef_ / X_scale
              self.intercept_ = y_offset - np.dot(X_offset, self.coef_.T)
          else:
              self.intercept_ = 0.
  
      def _more_tags(self):
          return {'requires_y': True}
  
  
  # XXX Should this derive from LinearModel? It should be a mixin, not an ABC.
  # Maybe the n_features checking can be moved to LinearModel.
  class LinearClassifierMixin(ClassifierMixin):
      """Mixin for linear classifiers.
  
      Handles prediction for sparse and dense X.
      """
  
      def decision_function(self, X):
          """
          Predict confidence scores for samples.
  
          The confidence score for a sample is proportional to the signed
          distance of that sample to the hyperplane.
  
          Parameters
          ----------
          X : array-like or sparse matrix, shape (n_samples, n_features)
              Samples.
  
          Returns
          -------
          array, shape=(n_samples,) if n_classes == 2 else (n_samples, n_classes)
              Confidence scores per (sample, class) combination. In the binary
              case, confidence score for self.classes_[1] where >0 means this
              class would be predicted.
          """
          check_is_fitted(self)
  
          X = check_array(X, accept_sparse='csr')
  
          n_features = self.coef_.shape[1]
          if X.shape[1] != n_features:
              raise ValueError("X has %d features per sample; expecting %d"
                               % (X.shape[1], n_features))
  
          scores = safe_sparse_dot(X, self.coef_.T,
                                   dense_output=True) + self.intercept_
          return scores.ravel() if scores.shape[1] == 1 else scores
  
      def predict(self, X):
          """
          Predict class labels for samples in X.
  
          Parameters
          ----------
          X : array-like or sparse matrix, shape (n_samples, n_features)
              Samples.
  
          Returns
          -------
          C : array, shape [n_samples]
              Predicted class label per sample.
          """
          scores = self.decision_function(X)
          if len(scores.shape) == 1:
              indices = (scores > 0).astype(int)
          else:
              indices = scores.argmax(axis=1)
          return self.classes_[indices]
  
      def _predict_proba_lr(self, X):
          """Probability estimation for OvR logistic regression.
  
          Positive class probabilities are computed as
          1. / (1. + np.exp(-self.decision_function(X)));
          multiclass is handled by normalizing that over all classes.
          """
          prob = self.decision_function(X)
          expit(prob, out=prob)
          if prob.ndim == 1:
              return np.vstack([1 - prob, prob]).T
          else:
              # OvR normalization, like LibLinear's predict_probability
              prob /= prob.sum(axis=1).reshape((prob.shape[0], -1))
              return prob
  
  
  class SparseCoefMixin:
      """Mixin for converting coef_ to and from CSR format.
  
      L1-regularizing estimators should inherit this.
      """
  
      def densify(self):
          """
          Convert coefficient matrix to dense array format.
  
          Converts the ``coef_`` member (back) to a numpy.ndarray. This is the
          default format of ``coef_`` and is required for fitting, so calling
          this method is only required on models that have previously been
          sparsified; otherwise, it is a no-op.
  
          Returns
          -------
          self
              Fitted estimator.
          """
          msg = "Estimator, %(name)s, must be fitted before densifying."
          check_is_fitted(self, msg=msg)
          if sp.issparse(self.coef_):
              self.coef_ = self.coef_.toarray()
          return self
  
      def sparsify(self):
          """
          Convert coefficient matrix to sparse format.
  
          Converts the ``coef_`` member to a scipy.sparse matrix, which for
          L1-regularized models can be much more memory- and storage-efficient
          than the usual numpy.ndarray representation.
  
          The ``intercept_`` member is not converted.
  
          Returns
          -------
          self
              Fitted estimator.
  
          Notes
          -----
          For non-sparse models, i.e. when there are not many zeros in ``coef_``,
          this may actually *increase* memory usage, so use this method with
          care. A rule of thumb is that the number of zero elements, which can
          be computed with ``(coef_ == 0).sum()``, must be more than 50% for this
          to provide significant benefits.
  
          After calling this method, further fitting with the partial_fit
          method (if any) will not work until you call densify.
          """
          msg = "Estimator, %(name)s, must be fitted before sparsifying."
          check_is_fitted(self, msg=msg)
          self.coef_ = sp.csr_matrix(self.coef_)
          return self
  
  
  class LinearRegression(MultiOutputMixin, RegressorMixin, LinearModel):
      """
      Ordinary least squares Linear Regression.
  
      LinearRegression fits a linear model with coefficients w = (w1, ..., wp)
      to minimize the residual sum of squares between the observed targets in
      the dataset, and the targets predicted by the linear approximation.
  
      Parameters
      ----------
      fit_intercept : bool, default=True
          Whether to calculate the intercept for this model. If set
          to False, no intercept will be used in calculations
          (i.e. data is expected to be centered).
  
      normalize : bool, default=False
          This parameter is ignored when ``fit_intercept`` is set to False.
          If True, the regressors X will be normalized before regression by
          subtracting the mean and dividing by the l2-norm.
          If you wish to standardize, please use
          :class:`~sklearn.preprocessing.StandardScaler` before calling ``fit``
          on an estimator with ``normalize=False``.
  
      copy_X : bool, default=True
          If True, X will be copied; else, it may be overwritten.
  
      n_jobs : int, default=None
          The number of jobs to use for the computation. This will only provide
          speedup for n_targets > 1 and sufficient large problems.
          ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
          ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
          for more details.
  
      positive : bool, default=False
          When set to ``True``, forces the coefficients to be positive. This
          option is only supported for dense arrays.
  
          .. versionadded:: 0.24
  
      Attributes
      ----------
      coef_ : array of shape (n_features, ) or (n_targets, n_features)
          Estimated coefficients for the linear regression problem.
          If multiple targets are passed during the fit (y 2D), this
          is a 2D array of shape (n_targets, n_features), while if only
          one target is passed, this is a 1D array of length n_features.
  
      rank_ : int
          Rank of matrix `X`. Only available when `X` is dense.
  
      singular_ : array of shape (min(X, y),)
          Singular values of `X`. Only available when `X` is dense.
  
      intercept_ : float or array of shape (n_targets,)
          Independent term in the linear model. Set to 0.0 if
          `fit_intercept = False`.
  
      See Also
      --------
      Ridge : Ridge regression addresses some of the
          problems of Ordinary Least Squares by imposing a penalty on the
          size of the coefficients with l2 regularization.
      Lasso : The Lasso is a linear model that estimates
          sparse coefficients with l1 regularization.
      ElasticNet : Elastic-Net is a linear regression
          model trained with both l1 and l2 -norm regularization of the
          coefficients.
  
      Notes
      -----
      From the implementation point of view, this is just plain Ordinary
      Least Squares (scipy.linalg.lstsq) or Non Negative Least Squares
      (scipy.optimize.nnls) wrapped as a predictor object.
  
      Examples
      --------
      >>> import numpy as np
      >>> from sklearn.linear_model import LinearRegression
      >>> X = np.array([[1, 1], [1, 2], [2, 2], [2, 3]])
      >>> # y = 1 * x_0 + 2 * x_1 + 3
      >>> y = np.dot(X, np.array([1, 2])) + 3
      >>> reg = LinearRegression().fit(X, y)
      >>> reg.score(X, y)
      1.0
      >>> reg.coef_
      array([1., 2.])
      >>> reg.intercept_
      3.0...
      >>> reg.predict(np.array([[3, 5]]))
      array([16.])
      """
      @_deprecate_positional_args
      def __init__(self, *, fit_intercept=True, normalize=False, copy_X=True,
                   n_jobs=None, positive=False):
          self.fit_intercept = fit_intercept
          self.normalize = normalize
          self.copy_X = copy_X
          self.n_jobs = n_jobs
          self.positive = positive
  
      def fit(self, X, y, sample_weight=None):
          """
          Fit linear model.
  
          Parameters
          ----------
          X : {array-like, sparse matrix} of shape (n_samples, n_features)
              Training data
  
          y : array-like of shape (n_samples,) or (n_samples, n_targets)
              Target values. Will be cast to X's dtype if necessary
  
          sample_weight : array-like of shape (n_samples,), default=None
              Individual weights for each sample
  
              .. versionadded:: 0.17
                 parameter *sample_weight* support to LinearRegression.
  
          Returns
          -------
          self : returns an instance of self.
          """
  
          n_jobs_ = self.n_jobs
  
          accept_sparse = False if self.positive else ['csr', 'csc', 'coo']
  
          X, y = self._validate_data(X, y, accept_sparse=accept_sparse,
                                     y_numeric=True, multi_output=True)
  
          if sample_weight is not None:
              sample_weight = _check_sample_weight(sample_weight, X,
                                                   dtype=X.dtype)
  
          X, y, X_offset, y_offset, X_scale = self._preprocess_data(
              X, y, fit_intercept=self.fit_intercept, normalize=self.normalize,
              copy=self.copy_X, sample_weight=sample_weight,
              return_mean=True)
  
          if sample_weight is not None:
              # Sample weight can be implemented via a simple rescaling.
              X, y = _rescale_data(X, y, sample_weight)
  
          if self.positive:
              if y.ndim < 2:
                  self.coef_, self._residues = optimize.nnls(X, y)
              else:
                  # scipy.optimize.nnls cannot handle y with shape (M, K)
                  outs = Parallel(n_jobs=n_jobs_)(
                          delayed(optimize.nnls)(X, y[:, j])
                          for j in range(y.shape[1]))
                  self.coef_, self._residues = map(np.vstack, zip(*outs))
          elif sp.issparse(X):
              X_offset_scale = X_offset / X_scale
  
              def matvec(b):
                  return X.dot(b) - b.dot(X_offset_scale)
  
              def rmatvec(b):
                  return X.T.dot(b) - X_offset_scale * np.sum(b)
  
              X_centered = sparse.linalg.LinearOperator(shape=X.shape,
                                                        matvec=matvec,
                                                        rmatvec=rmatvec)
  
              if y.ndim < 2:
                  out = sparse_lsqr(X_centered, y)
                  self.coef_ = out[0]
                  self._residues = out[3]
              else:
                  # sparse_lstsq cannot handle y with shape (M, K)
                  outs = Parallel(n_jobs=n_jobs_)(
                      delayed(sparse_lsqr)(X_centered, y[:, j].ravel())
                      for j in range(y.shape[1]))
                  self.coef_ = np.vstack([out[0] for out in outs])
                  self._residues = np.vstack([out[3] for out in outs])
          else:
              self.coef_, self._residues, self.rank_, self.singular_ = \
                  linalg.lstsq(X, y)
              self.coef_ = self.coef_.T
  
          if y.ndim == 1:
              self.coef_ = np.ravel(self.coef_)
          self._set_intercept(X_offset, y_offset, X_scale)
          return self
  
  
  def _pre_fit(X, y, Xy, precompute, normalize, fit_intercept, copy,
               check_input=True, sample_weight=None):
      """Aux function used at beginning of fit in linear models
  
      Parameters
      ----------
      order : 'F', 'C' or None, default=None
          Whether X and y will be forced to be fortran or c-style. Only relevant
          if sample_weight is not None.
      """
      n_samples, n_features = X.shape
  
      if sparse.isspmatrix(X):
          # copy is not needed here as X is not modified inplace when X is sparse
          precompute = False
          X, y, X_offset, y_offset, X_scale = _preprocess_data(
              X, y, fit_intercept=fit_intercept, normalize=normalize,
              copy=False, return_mean=True, check_input=check_input)
      else:
          # copy was done in fit if necessary
          X, y, X_offset, y_offset, X_scale = _preprocess_data(
              X, y, fit_intercept=fit_intercept, normalize=normalize, copy=copy,
              check_input=check_input, sample_weight=sample_weight)
      if sample_weight is not None:
          X, y = _rescale_data(X, y, sample_weight=sample_weight)
      if hasattr(precompute, '__array__') and (
          fit_intercept and not np.allclose(X_offset, np.zeros(n_features)) or
              normalize and not np.allclose(X_scale, np.ones(n_features))):
          warnings.warn("Gram matrix was provided but X was centered"
                        " to fit intercept, "
                        "or X was normalized : recomputing Gram matrix.",
                        UserWarning)
          # recompute Gram
          precompute = 'auto'
          Xy = None
  
      # precompute if n_samples > n_features
      if isinstance(precompute, str) and precompute == 'auto':
          precompute = (n_samples > n_features)
  
      if precompute is True:
          # make sure that the 'precompute' array is contiguous.
          precompute = np.empty(shape=(n_features, n_features), dtype=X.dtype,
                                order='C')
          np.dot(X.T, X, out=precompute)
  
      if not hasattr(precompute, '__array__'):
          Xy = None  # cannot use Xy if precompute is not Gram
  
      if hasattr(precompute, '__array__') and Xy is None:
          common_dtype = np.find_common_type([X.dtype, y.dtype], [])
          if y.ndim == 1:
              # Xy is 1d, make sure it is contiguous.
              Xy = np.empty(shape=n_features, dtype=common_dtype, order='C')
              np.dot(X.T, y, out=Xy)
          else:
              # Make sure that Xy is always F contiguous even if X or y are not
              # contiguous: the goal is to make it fast to extract the data for a
              # specific target.
              n_targets = y.shape[1]
              Xy = np.empty(shape=(n_features, n_targets), dtype=common_dtype,
                            order='F')
              np.dot(y.T, X, out=Xy.T)
  
      return X, y, X_offset, y_offset, X_scale, precompute, Xy
   """
  Testing for the partial dependence module.
  """
  
  import numpy as np
  import pytest
  
  import sklearn
  from sklearn.inspection import partial_dependence
  from sklearn.inspection._partial_dependence import (
      _grid_from_X,
      _partial_dependence_brute,
      _partial_dependence_recursion
  )
  from sklearn.ensemble import GradientBoostingClassifier
  from sklearn.ensemble import GradientBoostingRegressor
  from sklearn.ensemble import RandomForestRegressor
  from sklearn.experimental import enable_hist_gradient_boosting  # noqa
  from sklearn.ensemble import HistGradientBoostingClassifier
  from sklearn.ensemble import HistGradientBoostingRegressor
  from sklearn.linear_model import LinearRegression
  from sklearn.linear_model import LogisticRegression
  from sklearn.linear_model import MultiTaskLasso
  from sklearn.tree import DecisionTreeRegressor
  from sklearn.datasets import load_iris
  from sklearn.datasets import make_classification, make_regression
  from sklearn.cluster import KMeans
  from sklearn.compose import make_column_transformer
  from sklearn.metrics import r2_score
  from sklearn.preprocessing import PolynomialFeatures
  from sklearn.preprocessing import StandardScaler
  from sklearn.preprocessing import RobustScaler
  from sklearn.preprocessing import scale
  from sklearn.pipeline import make_pipeline
  from sklearn.dummy import DummyClassifier
  from sklearn.base import BaseEstimator, ClassifierMixin, clone
  from sklearn.exceptions import NotFittedError
  from sklearn.utils._testing import assert_allclose
  from sklearn.utils._testing import assert_array_equal
  from sklearn.utils import _IS_32BIT
  from sklearn.utils.validation import check_random_state
  from sklearn.tree.tests.test_tree import assert_is_subtree
  
  
  # toy sample
  X = [[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]]
  y = [-1, -1, -1, 1, 1, 1]
  
  
  # (X, y), n_targets  <-- as expected in the output of partial_dep()
  binary_classification_data = (make_classification(n_samples=50,
                                                    random_state=0), 1)
  multiclass_classification_data = (make_classification(n_samples=50,
                                                        n_classes=3,
                                                        n_clusters_per_class=1,
                                                        random_state=0), 3)
  regression_data = (make_regression(n_samples=50, random_state=0), 1)
  multioutput_regression_data = (make_regression(n_samples=50, n_targets=2,
                                                 random_state=0), 2)
  
  # iris
  iris = load_iris()
  
  
  @pytest.mark.filterwarnings("ignore:A Bunch will be returned")
  @pytest.mark.parametrize('Estimator, method, data', [
      (GradientBoostingClassifier, 'auto', binary_classification_data),
      (GradientBoostingClassifier, 'auto', multiclass_classification_data),
      (GradientBoostingClassifier, 'brute', binary_classification_data),
      (GradientBoostingClassifier, 'brute', multiclass_classification_data),
      (GradientBoostingRegressor, 'auto', regression_data),
      (GradientBoostingRegressor, 'brute', regression_data),
      (DecisionTreeRegressor, 'brute', regression_data),
      (LinearRegression, 'brute', regression_data),
      (LinearRegression, 'brute', multioutput_regression_data),
      (LogisticRegression, 'brute', binary_classification_data),
      (LogisticRegression, 'brute', multiclass_classification_data),
      (MultiTaskLasso, 'brute', multioutput_regression_data),
      ])
  @pytest.mark.parametrize('grid_resolution', (5, 10))
  @pytest.mark.parametrize('features', ([1], [1, 2]))
  @pytest.mark.parametrize('kind', ('legacy', 'average', 'individual', 'both'))
  def test_output_shape(Estimator, method, data, grid_resolution,
                        features, kind):
      # Check that partial_dependence has consistent output shape for different
      # kinds of estimators:
      # - classifiers with binary and multiclass settings
      # - regressors
      # - multi-task regressors
  
      est = Estimator()
  
      # n_target corresponds to the number of classes (1 for binary classif) or
      # the number of tasks / outputs in multi task settings. It's equal to 1 for
      # classical regression_data.
      (X, y), n_targets = data
      n_instances = X.shape[0]
  
      est.fit(X, y)
      result = partial_dependence(
          est, X=X, features=features, method=method, kind=kind,
          grid_resolution=grid_resolution
      )
      # FIXME: Remove 'legacy' support in 1.1
      pdp, axes = result if kind == 'legacy' else (result, result["values"])
  
      expected_pdp_shape = (n_targets,
                            *[grid_resolution for _ in range(len(features))])
      expected_ice_shape = (n_targets, n_instances,
                            *[grid_resolution for _ in range(len(features))])
      if kind == 'legacy':
          assert pdp.shape == expected_pdp_shape
      elif kind == 'average':
          assert pdp.average.shape == expected_pdp_shape
      elif kind == 'individual':
          assert pdp.individual.shape == expected_ice_shape
      else:  # 'both'
          assert pdp.average.shape == expected_pdp_shape
          assert pdp.individual.shape == expected_ice_shape
  
      expected_axes_shape = (len(features), grid_resolution)
      assert axes is not None
      assert np.asarray(axes).shape == expected_axes_shape
  
  
  def test_grid_from_X():
      # tests for _grid_from_X: sanity check for output, and for shapes.
  
      # Make sure that the grid is a cartesian product of the input (it will use
      # the unique values instead of the percentiles)
      percentiles = (.05, .95)
      grid_resolution = 100
      X = np.asarray([[1, 2],
                      [3, 4]])
      grid, axes = _grid_from_X(X, percentiles, grid_resolution)
      assert_array_equal(grid, [[1, 2],
                                [1, 4],
                                [3, 2],
                                [3, 4]])
      assert_array_equal(axes, X.T)
  
      # test shapes of returned objects depending on the number of unique values
      # for a feature.
      rng = np.random.RandomState(0)
      grid_resolution = 15
  
      # n_unique_values > grid_resolution
      X = rng.normal(size=(20, 2))
      grid, axes = _grid_from_X(X, percentiles, grid_resolution=grid_resolution)
      assert grid.shape == (grid_resolution * grid_resolution, X.shape[1])
      assert np.asarray(axes).shape == (2, grid_resolution)
  
      # n_unique_values < grid_resolution, will use actual values
      n_unique_values = 12
      X[n_unique_values - 1:, 0] = 12345
      rng.shuffle(X)  # just to make sure the order is irrelevant
      grid, axes = _grid_from_X(X, percentiles, grid_resolution=grid_resolution)
      assert grid.shape == (n_unique_values * grid_resolution, X.shape[1])
      # axes is a list of arrays of different shapes
      assert axes[0].shape == (n_unique_values,)
      assert axes[1].shape == (grid_resolution,)
  
  
  @pytest.mark.parametrize(
      "grid_resolution, percentiles, err_msg",
      [(2, (0, 0.0001), "percentiles are too close"),
       (100, (1, 2, 3, 4), "'percentiles' must be a sequence of 2 elements"),
       (100, 12345, "'percentiles' must be a sequence of 2 elements"),
       (100, (-1, .95), r"'percentiles' values must be in \[0, 1\]"),
       (100, (.05, 2), r"'percentiles' values must be in \[0, 1\]"),
       (100, (.9, .1), r"percentiles\[0\] must be strictly less than"),
       (1, (0.05, 0.95), "'grid_resolution' must be strictly greater than 1")]
  )
  def test_grid_from_X_error(grid_resolution, percentiles, err_msg):
      X = np.asarray([[1, 2], [3, 4]])
      with pytest.raises(ValueError, match=err_msg):
          _grid_from_X(
              X, grid_resolution=grid_resolution, percentiles=percentiles
          )
  
  
  @pytest.mark.parametrize('target_feature', range(5))
  @pytest.mark.parametrize('est, method', [
      (LinearRegression(), 'brute'),
      (GradientBoostingRegressor(random_state=0), 'brute'),
      (GradientBoostingRegressor(random_state=0), 'recursion'),
      (HistGradientBoostingRegressor(random_state=0), 'brute'),
      (HistGradientBoostingRegressor(random_state=0), 'recursion')]
  )
  def test_partial_dependence_helpers(est, method, target_feature):
      # Check that what is returned by _partial_dependence_brute or
      # _partial_dependence_recursion is equivalent to manually setting a target
      # feature to a given value, and computing the average prediction over all
      # samples.
      # This also checks that the brute and recursion methods give the same
      # output.
      # Note that even on the trainset, the brute and the recursion methods
      # aren't always strictly equivalent, in particular when the slow method
      # generates unrealistic samples that have low mass in the joint
      # distribution of the input features, and when some of the features are
      # dependent. Hence the high tolerance on the checks.
  
      X, y = make_regression(random_state=0, n_features=5, n_informative=5)
      # The 'init' estimator for GBDT (here the average prediction) isn't taken
      # into account with the recursion method, for technical reasons. We set
      # the mean to 0 to that this 'bug' doesn't have any effect.
      y = y - y.mean()
      est.fit(X, y)
  
      # target feature will be set to .5 and then to 123
      features = np.array([target_feature], dtype=np.int32)
      grid = np.array([[.5],
                       [123]])
  
      if method == 'brute':
          pdp, predictions = _partial_dependence_brute(est, grid, features, X,
                                                       response_method='auto')
      else:
          pdp = _partial_dependence_recursion(est, grid, features)
  
      mean_predictions = []
      for val in (.5, 123):
          X_ = X.copy()
          X_[:, target_feature] = val
          mean_predictions.append(est.predict(X_).mean())
  
      pdp = pdp[0]  # (shape is (1, 2) so make it (2,))
  
      # allow for greater margin for error with recursion method
      rtol = 1e-1 if method == 'recursion' else 1e-3
      assert np.allclose(pdp, mean_predictions, rtol=rtol)
  
  
  @pytest.mark.parametrize('seed', range(1))
  def test_recursion_decision_tree_vs_forest_and_gbdt(seed):
      # Make sure that the recursion method gives the same results on a
      # DecisionTreeRegressor and a GradientBoostingRegressor or a
      # RandomForestRegressor with 1 tree and equivalent parameters.
  
      rng = np.random.RandomState(seed)
  
      # Purely random dataset to avoid correlated features
      n_samples = 1000
      n_features = 5
      X = rng.randn(n_samples, n_features)
      y = rng.randn(n_samples) * 10
  
      # The 'init' estimator for GBDT (here the average prediction) isn't taken
      # into account with the recursion method, for technical reasons. We set
      # the mean to 0 to that this 'bug' doesn't have any effect.
      y = y - y.mean()
  
      # set max_depth not too high to avoid splits with same gain but different
      # features
      max_depth = 5
  
      tree_seed = 0
      forest = RandomForestRegressor(n_estimators=1, max_features=None,
                                     bootstrap=False, max_depth=max_depth,
                                     random_state=tree_seed)
      # The forest will use ensemble.base._set_random_states to set the
      # random_state of the tree sub-estimator. We simulate this here to have
      # equivalent estimators.
      equiv_random_state = check_random_state(tree_seed).randint(
          np.iinfo(np.int32).max)
      gbdt = GradientBoostingRegressor(n_estimators=1, learning_rate=1,
                                       criterion='mse', max_depth=max_depth,
                                       random_state=equiv_random_state)
      tree = DecisionTreeRegressor(max_depth=max_depth,
                                   random_state=equiv_random_state)
  
      forest.fit(X, y)
      gbdt.fit(X, y)
      tree.fit(X, y)
  
      # sanity check: if the trees aren't the same, the PD values won't be equal
      try:
          assert_is_subtree(tree.tree_, gbdt[0, 0].tree_)
          assert_is_subtree(tree.tree_, forest[0].tree_)
      except AssertionError:
          # For some reason the trees aren't exactly equal on 32bits, so the PDs
          # cannot be equal either. See
          # https://github.com/scikit-learn/scikit-learn/issues/8853
          assert _IS_32BIT, "this should only fail on 32 bit platforms"
          return
  
      grid = rng.randn(50).reshape(-1, 1)
      for f in range(n_features):
          features = np.array([f], dtype=np.int32)
  
          pdp_forest = _partial_dependence_recursion(forest, grid, features)
          pdp_gbdt = _partial_dependence_recursion(gbdt, grid, features)
          pdp_tree = _partial_dependence_recursion(tree, grid, features)
  
          np.testing.assert_allclose(pdp_gbdt, pdp_tree)
          np.testing.assert_allclose(pdp_forest, pdp_tree)
  
  
  @pytest.mark.parametrize('est', (
      GradientBoostingClassifier(random_state=0),
      HistGradientBoostingClassifier(random_state=0),
  ))
  @pytest.mark.parametrize('target_feature', (0, 1, 2, 3, 4, 5))
  def test_recursion_decision_function(est, target_feature):
      # Make sure the recursion method (implicitly uses decision_function) has
      # the same result as using brute method with
      # response_method=decision_function
  
      X, y = make_classification(n_classes=2, n_clusters_per_class=1,
                                 random_state=1)
      assert np.mean(y) == .5  # make sure the init estimator predicts 0 anyway
  
      est.fit(X, y)
  
      preds_1 = partial_dependence(
          est, X, [target_feature], response_method='decision_function',
          method='recursion', kind='average'
      )
      preds_2 = partial_dependence(
          est, X, [target_feature], response_method='decision_function',
          method='brute', kind='average'
      )
  
      assert_allclose(preds_1['average'], preds_2['average'], atol=1e-7)
  
  
  @pytest.mark.parametrize('est', (
      LinearRegression(),
      GradientBoostingRegressor(random_state=0),
      HistGradientBoostingRegressor(random_state=0, min_samples_leaf=1,
                                    max_leaf_nodes=None, max_iter=1),
      DecisionTreeRegressor(random_state=0),
  ))
  @pytest.mark.parametrize('power', (1, 2))
  def test_partial_dependence_easy_target(est, power):
      # If the target y only depends on one feature in an obvious way (linear or
      # quadratic) then the partial dependence for that feature should reflect
      # it.
      # We here fit a linear regression_data model (with polynomial features if
      # needed) and compute r_squared to check that the partial dependence
      # correctly reflects the target.
  
      rng = np.random.RandomState(0)
      n_samples = 200
      target_variable = 2
      X = rng.normal(size=(n_samples, 5))
      y = X[:, target_variable]**power
  
      est.fit(X, y)
  
      pdp = partial_dependence(
          est, features=[target_variable], X=X, grid_resolution=1000,
          kind='average'
      )
  
      new_X = pdp["values"][0].reshape(-1, 1)
      new_y = pdp['average'][0]
      # add polynomial features if needed
      new_X = PolynomialFeatures(degree=power).fit_transform(new_X)
  
      lr = LinearRegression().fit(new_X, new_y)
      r2 = r2_score(new_y, lr.predict(new_X))
  
      assert r2 > .99
  
  
  @pytest.mark.parametrize('Estimator',
                           (sklearn.tree.DecisionTreeClassifier,
                            sklearn.tree.ExtraTreeClassifier,
                            sklearn.ensemble.ExtraTreesClassifier,
                            sklearn.neighbors.KNeighborsClassifier,
                            sklearn.neighbors.RadiusNeighborsClassifier,
                            sklearn.ensemble.RandomForestClassifier))
  def test_multiclass_multioutput(Estimator):
      # Make sure error is raised for multiclass-multioutput classifiers
  
      # make multiclass-multioutput dataset
      X, y = make_classification(n_classes=3, n_clusters_per_class=1,
                                 random_state=0)
      y = np.array([y, y]).T
  
      est = Estimator()
      est.fit(X, y)
  
      with pytest.raises(
              ValueError,
              match="Multiclass-multioutput estimators are not supported"):
          partial_dependence(est, X, [0])
  
  
  class NoPredictProbaNoDecisionFunction(ClassifierMixin, BaseEstimator):
      def fit(self, X, y):
          # simulate that we have some classes
          self.classes_ = [0, 1]
          return self
  
  
  @pytest.mark.filterwarnings("ignore:A Bunch will be returned")
  @pytest.mark.parametrize(
      "estimator, params, err_msg",
      [(KMeans(),
        {'features': [0]},
        "'estimator' must be a fitted regressor or classifier"),
       (LinearRegression(),
        {'features': [0], 'response_method': 'predict_proba'},
        'The response_method parameter is ignored for regressors'),
       (GradientBoostingClassifier(random_state=0),
        {'features': [0], 'response_method': 'predict_proba',
         'method': 'recursion'},
        "'recursion' method, the response_method must be 'decision_function'"),
       (GradientBoostingClassifier(random_state=0),
        {'features': [0], 'response_method': 'predict_proba', 'method': 'auto'},
        "'recursion' method, the response_method must be 'decision_function'"),
       (GradientBoostingClassifier(random_state=0),
        {'features': [0], 'response_method': 'blahblah'},
        'response_method blahblah is invalid. Accepted response_method'),
       (NoPredictProbaNoDecisionFunction(),
        {'features': [0], 'response_method': 'auto'},
        'The estimator has no predict_proba and no decision_function method'),
       (NoPredictProbaNoDecisionFunction(),
        {'features': [0], 'response_method': 'predict_proba'},
        'The estimator has no predict_proba method.'),
       (NoPredictProbaNoDecisionFunction(),
        {'features': [0], 'response_method': 'decision_function'},
        'The estimator has no decision_function method.'),
       (LinearRegression(),
        {'features': [0], 'method': 'blahblah'},
        'blahblah is invalid. Accepted method names are brute, recursion, auto'),
       (LinearRegression(),
        {'features': [0], 'method': 'recursion', 'kind': 'individual'},
        "The 'recursion' method only applies when 'kind' is set to 'average'"),
       (LinearRegression(),
        {'features': [0], 'method': 'recursion', 'kind': 'both'},
        "The 'recursion' method only applies when 'kind' is set to 'average'"),
       (LinearRegression(),
        {'features': [0], 'method': 'recursion'},
        "Only the following estimators support the 'recursion' method:")]
  )
  def test_partial_dependence_error(estimator, params, err_msg):
      X, y = make_classification(random_state=0)
      estimator.fit(X, y)
  
      with pytest.raises(ValueError, match=err_msg):
          partial_dependence(estimator, X, **params)
  
  
  @pytest.mark.parametrize(
      "with_dataframe, err_msg",
      [(True, "Only array-like or scalar are supported"),
       (False, "Only array-like or scalar are supported")]
  )
  def test_partial_dependence_slice_error(with_dataframe, err_msg):
      X, y = make_classification(random_state=0)
      if with_dataframe:
          pd = pytest.importorskip('pandas')
          X = pd.DataFrame(X)
      estimator = LogisticRegression().fit(X, y)
  
      with pytest.raises(TypeError, match=err_msg):
          partial_dependence(estimator, X, features=slice(0, 2, 1))
  
  
  @pytest.mark.parametrize(
      'estimator',
      [LinearRegression(), GradientBoostingClassifier(random_state=0)]
  )
  @pytest.mark.parametrize('features', [-1, 10000])
  def test_partial_dependence_unknown_feature_indices(estimator, features):
      X, y = make_classification(random_state=0)
      estimator.fit(X, y)
  
      err_msg = 'all features must be in'
      with pytest.raises(ValueError, match=err_msg):
          partial_dependence(estimator, X, [features])
  
  
  @pytest.mark.parametrize(
      'estimator',
      [LinearRegression(), GradientBoostingClassifier(random_state=0)]
  )
  def test_partial_dependence_unknown_feature_string(estimator):
      pd = pytest.importorskip("pandas")
      X, y = make_classification(random_state=0)
      df = pd.DataFrame(X)
      estimator.fit(df, y)
  
      features = ['random']
      err_msg = 'A given column is not a column of the dataframe'
      with pytest.raises(ValueError, match=err_msg):
          partial_dependence(estimator, df, features)
  
  
  @pytest.mark.parametrize(
      'estimator',
      [LinearRegression(), GradientBoostingClassifier(random_state=0)]
  )
  def test_partial_dependence_X_list(estimator):
      # check that array-like objects are accepted
      X, y = make_classification(random_state=0)
      estimator.fit(X, y)
      partial_dependence(estimator, list(X), [0], kind='average')
  
  
  def test_warning_recursion_non_constant_init():
      # make sure that passing a non-constant init parameter to a GBDT and using
      # recursion method yields a warning.
  
      gbc = GradientBoostingClassifier(init=DummyClassifier(), random_state=0)
      gbc.fit(X, y)
  
      with pytest.warns(
              UserWarning,
              match='Using recursion method with a non-constant init predictor'):
          partial_dependence(gbc, X, [0], method='recursion', kind='average')
  
      with pytest.warns(
              UserWarning,
              match='Using recursion method with a non-constant init predictor'):
          partial_dependence(gbc, X, [0], method='recursion', kind='average')
  
  
  def test_partial_dependence_sample_weight():
      # Test near perfect correlation between partial dependence and diagonal
      # when sample weights emphasize y = x predictions
      # non-regression test for #13193
      # TODO: extend to HistGradientBoosting once sample_weight is supported
      N = 1000
      rng = np.random.RandomState(123456)
      mask = rng.randint(2, size=N, dtype=bool)
  
      x = rng.rand(N)
      # set y = x on mask and y = -x outside
      y = x.copy()
      y[~mask] = -y[~mask]
      X = np.c_[mask, x]
      # sample weights to emphasize data points where y = x
      sample_weight = np.ones(N)
      sample_weight[mask] = 1000.
  
      clf = GradientBoostingRegressor(n_estimators=10, random_state=1)
      clf.fit(X, y, sample_weight=sample_weight)
  
      pdp = partial_dependence(clf, X, features=[1], kind='average')
  
      assert np.corrcoef(pdp['average'], pdp["values"])[0, 1] > 0.99
  
  
  def test_hist_gbdt_sw_not_supported():
      # TODO: remove/fix when PDP supports HGBT with sample weights
      clf = HistGradientBoostingRegressor(random_state=1)
      clf.fit(X, y, sample_weight=np.ones(len(X)))
  
      with pytest.raises(NotImplementedError,
                         match="does not support partial dependence"):
          partial_dependence(clf, X, features=[1])
  
  
  def test_partial_dependence_pipeline():
      # check that the partial dependence support pipeline
      iris = load_iris()
  
      scaler = StandardScaler()
      clf = DummyClassifier(random_state=42)
      pipe = make_pipeline(scaler, clf)
  
      clf.fit(scaler.fit_transform(iris.data), iris.target)
      pipe.fit(iris.data, iris.target)
  
      features = 0
      pdp_pipe = partial_dependence(
          pipe, iris.data, features=[features], grid_resolution=10,
          kind='average'
      )
      pdp_clf = partial_dependence(
          clf, scaler.transform(iris.data), features=[features],
          grid_resolution=10, kind='average'
      )
      assert_allclose(pdp_pipe['average'], pdp_clf['average'])
      assert_allclose(
          pdp_pipe["values"][0],
          pdp_clf["values"][0] * scaler.scale_[features] + scaler.mean_[features]
      )
  
  
  @pytest.mark.parametrize(
      "estimator",
      [LogisticRegression(max_iter=1000, random_state=0),
       GradientBoostingClassifier(random_state=0, n_estimators=5)],
      ids=['estimator-brute', 'estimator-recursion']
  )
  @pytest.mark.parametrize(
      "preprocessor",
      [None,
       make_column_transformer(
           (StandardScaler(), [iris.feature_names[i] for i in (0, 2)]),
           (RobustScaler(), [iris.feature_names[i] for i in (1, 3)])),
       make_column_transformer(
           (StandardScaler(), [iris.feature_names[i] for i in (0, 2)]),
           remainder='passthrough')],
      ids=['None', 'column-transformer', 'column-transformer-passthrough']
  )
  @pytest.mark.parametrize(
      "features",
      [[0, 2], [iris.feature_names[i] for i in (0, 2)]],
      ids=['features-integer', 'features-string']
  )
  def test_partial_dependence_dataframe(estimator, preprocessor, features):
      # check that the partial dependence support dataframe and pipeline
      # including a column transformer
      pd = pytest.importorskip("pandas")
      df = pd.DataFrame(scale(iris.data), columns=iris.feature_names)
  
      pipe = make_pipeline(preprocessor, estimator)
      pipe.fit(df, iris.target)
      pdp_pipe = partial_dependence(
          pipe, df, features=features, grid_resolution=10, kind='average'
      )
  
      # the column transformer will reorder the column when transforming
      # we mixed the index to be sure that we are computing the partial
      # dependence of the right columns
      if preprocessor is not None:
          X_proc = clone(preprocessor).fit_transform(df)
          features_clf = [0, 1]
      else:
          X_proc = df
          features_clf = [0, 2]
  
      clf = clone(estimator).fit(X_proc, iris.target)
      pdp_clf = partial_dependence(
          clf, X_proc, features=features_clf, method='brute', grid_resolution=10,
          kind='average'
      )
  
      assert_allclose(pdp_pipe['average'], pdp_clf['average'])
      if preprocessor is not None:
          scaler = preprocessor.named_transformers_['standardscaler']
          assert_allclose(
              pdp_pipe["values"][1],
              pdp_clf["values"][1] * scaler.scale_[1] + scaler.mean_[1]
          )
      else:
          assert_allclose(pdp_pipe["values"][1], pdp_clf["values"][1])
  
  
  @pytest.mark.parametrize(
      "features, expected_pd_shape",
      [(0, (3, 10)),
       (iris.feature_names[0], (3, 10)),
       ([0, 2], (3, 10, 10)),
       ([iris.feature_names[i] for i in (0, 2)], (3, 10, 10)),
       ([True, False, True, False], (3, 10, 10))],
      ids=['scalar-int', 'scalar-str', 'list-int', 'list-str', 'mask']
  )
  def test_partial_dependence_feature_type(features, expected_pd_shape):
      # check all possible features type supported in PDP
      pd = pytest.importorskip("pandas")
      df = pd.DataFrame(iris.data, columns=iris.feature_names)
  
      preprocessor = make_column_transformer(
          (StandardScaler(), [iris.feature_names[i] for i in (0, 2)]),
          (RobustScaler(), [iris.feature_names[i] for i in (1, 3)])
      )
      pipe = make_pipeline(
          preprocessor, LogisticRegression(max_iter=1000, random_state=0)
      )
      pipe.fit(df, iris.target)
      pdp_pipe = partial_dependence(
          pipe, df, features=features, grid_resolution=10, kind='average'
      )
      assert pdp_pipe['average'].shape == expected_pd_shape
      assert len(pdp_pipe["values"]) == len(pdp_pipe['average'].shape) - 1
  
  
  @pytest.mark.parametrize(
      "estimator", [LinearRegression(), LogisticRegression(),
                    GradientBoostingRegressor(), GradientBoostingClassifier()]
  )
  def test_partial_dependence_unfitted(estimator):
      X = iris.data
      preprocessor = make_column_transformer(
          (StandardScaler(), [0, 2]), (RobustScaler(), [1, 3])
      )
      pipe = make_pipeline(preprocessor, estimator)
      with pytest.raises(NotFittedError, match="is not fitted yet"):
          partial_dependence(pipe, X, features=[0, 2], grid_resolution=10)
      with pytest.raises(NotFittedError, match="is not fitted yet"):
          partial_dependence(estimator, X, features=[0, 2], grid_resolution=10)
  
  
  @pytest.mark.parametrize('Estimator, data', [
      (LinearRegression, multioutput_regression_data),
      (LogisticRegression, binary_classification_data)])
  def test_kind_average_and_average_of_individual(Estimator, data):
      est = Estimator()
      (X, y), n_targets = data
      est.fit(X, y)
  
      pdp_avg = partial_dependence(
              est, X=X, features=[1, 2], kind='average'
      )
      pdp_ind = partial_dependence(
          est, X=X, features=[1, 2], kind='individual'
      )
      avg_ind = np.mean(pdp_ind['individual'], axis=1)
      assert_allclose(avg_ind, pdp_avg['average'])
  
  
  def test_warning_for_kind_legacy():
      est = LogisticRegression()
      (X, y), n_targets = binary_classification_data
      est.fit(X, y)
  
      err_msg = ("A Bunch will be returned in place of 'predictions' from "
                 "version 1.1")
      with pytest.warns(FutureWarning, match=err_msg):
          partial_dependence(est, X=X, features=[1, 2])
  
      with pytest.warns(FutureWarning, match=err_msg):
          partial_dependence(est, X=X, features=[1, 2], kind='legacy')
  
apis: numpy.ones numpy.sum numpy.find_common_type scipy.linalg.lstsq numpy.asarray scipy.sparse.linalg.LinearOperator numpy.ascontiguousarray numpy.vstack scipy.sparse.dia_matrix scipy.optimize.nnls numpy.average numpy.zeros numpy.empty scipy.sparse.issparse scipy.sparse.csr_matrix numpy.ravel numpy.iinfo scipy.special.expit scipy.sparse.isspmatrix numpy.sqrt numpy.dot numpy.full
apis: numpy.ones sklearn.inspection._partial_dependence._partial_dependence_brute sklearn.inspection.partial_dependence sklearn.tree.DecisionTreeRegressor numpy.asarray sklearn.datasets.make_regression numpy.random.RandomState sklearn.ensemble.RandomForestRegressor sklearn.cluster.KMeans sklearn.linear_model.LogisticRegression sklearn.inspection._partial_dependence._partial_dependence_recursion sklearn.ensemble.GradientBoostingRegressor sklearn.datasets.load_iris numpy.allclose sklearn.base.clone sklearn.utils._testing.assert_array_equal sklearn.utils.validation.check_random_state sklearn.linear_model.LinearRegression sklearn.preprocessing.scale sklearn.tree.tests.test_tree.assert_is_subtree sklearn.preprocessing.RobustScaler sklearn.ensemble.HistGradientBoostingClassifier sklearn.dummy.DummyClassifier numpy.mean numpy.corrcoef sklearn.utils._testing.assert_allclose sklearn.pipeline.make_pipeline sklearn.preprocessing.StandardScaler sklearn.ensemble.HistGradientBoostingRegressor sklearn.datasets.make_classification sklearn.inspection._partial_dependence._grid_from_X numpy.iinfo numpy.testing.assert_allclose numpy.array sklearn.ensemble.GradientBoostingClassifier sklearn.preprocessing.PolynomialFeatures