text: """
  ===============================
  Stability selection transformer
  ===============================
  
  This module contains a scikit-learn compatible implementation of
  stability selection [1]_ .
  
  References
  ----------
  .. [1] Meinshausen, N. and Buhlmann, P., 2010. Stability selection.
      Journal of the Royal Statistical Society: Series B
      (Statistical Methodology), 72(4), pp.417-473.
  
  .. [2] Shah, R.D. and Samworth, R.J., 2013. Variable selection with
     error control: another look at stability selection. Journal
     of the Royal Statistical Society: Series B (Statistical Methodology),
      75(1), pp.55-80.
  """
  
  from warnings import warn
  
  import matplotlib.pyplot as plt
  import numpy as np
  from sklearn.base import BaseEstimator, TransformerMixin, clone
  from sklearn.externals.joblib import Parallel, delayed
  from sklearn.feature_selection import SelectFromModel
  from sklearn.linear_model import LogisticRegression
  from sklearn.pipeline import Pipeline
  from sklearn.utils import check_array, check_random_state, check_X_y, safe_mask
  from sklearn.utils.validation import check_is_fitted
  
  from .bootstrap import (bootstrap_without_replacement,
                          complementary_pairs_bootstrap, stratified_bootstrap)
  
  __all__ = ['StabilitySelection', 'plot_stability_path']
  
  BOOTSTRAP_FUNC_MAPPING = {
      'subsample': bootstrap_without_replacement,
      'complementary_pairs': complementary_pairs_bootstrap,
      'stratified': stratified_bootstrap
  }
  
  
  def _return_estimator_from_pipeline(pipeline):
      """Returns the final estimator in a Pipeline, or the estimator
      if it is not"""
      if isinstance(pipeline, Pipeline):
          return pipeline._final_estimator
      else:
          return pipeline
  
  
  def _return_support_from_pipeline(pipeline, variable_selector):
      """Returns the support of a Pipeline after variable selection"""
      support = variable_selector.get_support()
      if isinstance(pipeline, Pipeline):
          reverse_iter = reversed(list(pipeline._iter()))
          for idx, name, transform in reverse_iter:
              if hasattr(transform, 'get_support'):
                  temp = transform.get_support()
                  temp[transform.get_support(indices=True)] = support
                  support = temp
      return support
  
  
  def _bootstrap_generator(n_bootstrap_iterations, bootstrap_func, y,
                           n_subsamples, random_state=None):
      for _ in range(n_bootstrap_iterations):
          subsample = bootstrap_func(y, n_subsamples, random_state)
          if isinstance(subsample, tuple):
              for item in subsample:
                  yield item
          else:
              yield subsample
  
  
  def _fit_bootstrap_sample(base_estimator, X, y, lambda_name, lambda_value,
                            threshold=None):
      """
      Fits base_estimator on a bootstrap sample of the original data,
      and returns a mas of the variables that are selected by the fitted model.
  
      Parameters
      ----------
      base_estimator : Estimator
          Estimator to be fitted on the data
  
      X : {array-like, sparse matrix}, shape = [n_samples, n_features]
          The training input samples.
  
      y : array-like, shape = [n_samples]
          The target values.
  
      lambda_name : str
          Name of the penalization parameter of base_estimator
  
      lambda_value : float
          Value of the penalization parameter
  
      threshold : string, float, optional default None
          The threshold value to use for feature selection. Features whose
          importance is greater or equal are kept while the others are
          discarded. If "median" (resp. "mean"), then the ``threshold`` value is
          the median (resp. the mean) of the feature importances. A scaling
          factor (e.g., "1.25*mean") may also be used. If None and if the
          estimator has a parameter penalty set to l1, either explicitly
          or implicitly (e.g, Lasso), the threshold used is 1e-5.
          Otherwise, "mean" is used by default.
  
      Returns
      -------
      selected_variables : array-like, shape = [n_features]
          Boolean mask of selected variables.
      """
  
      base_estimator.set_params(**{lambda_name: lambda_value})
      base_estimator.fit(X, y)
  
      # TODO: Reconsider if we really want to use SelectFromModel here or not
      selector_model = _return_estimator_from_pipeline(base_estimator)
      variable_selector = SelectFromModel(estimator=selector_model,
                                          threshold=threshold,
                                          prefit=True)
      return _return_support_from_pipeline(base_estimator, variable_selector)
  
  
  def plot_stability_path(stability_selection, threshold_highlight=None,
                          **kwargs):
      """Plots stability path.
  
      Parameters
      ----------
      stability_selection : StabilitySelection
          Fitted instance of StabilitySelection.
  
      threshold_highlight : float
          Threshold defining the cutoff for the stability scores for the
          variables that need to be highlighted.
  
      kwargs : dict
          Arguments passed to matplotlib plot function.
      """
      check_is_fitted(stability_selection, 'stability_scores_')
  
      threshold = stability_selection.threshold if threshold_highlight is None else threshold_highlight
      paths_to_highlight = stability_selection.get_support(threshold=threshold)
  
      x_grid = stability_selection.lambda_grid / np.max(stability_selection.lambda_grid)
  
      fig, ax = plt.subplots(1, 1, **kwargs)
      if not paths_to_highlight.all():
          ax.plot(x_grid, stability_selection.stability_scores_[~paths_to_highlight].T,
                  'k:', linewidth=0.5)
  
      if paths_to_highlight.any():
          ax.plot(x_grid, stability_selection.stability_scores_[paths_to_highlight].T,
                  'r-', linewidth=0.5)
  
      if threshold is not None:
          ax.plot(x_grid, threshold * np.ones_like(stability_selection.lambda_grid),
                  'b--', linewidth=0.5)
  
      ax.set_ylabel('Stability score')
      ax.set_xlabel('Lambda / max(Lambda)')
  
      fig.tight_layout()
  
      return fig, ax
  
  
  class StabilitySelection(BaseEstimator, TransformerMixin):
      """Stability selection [1]_ fits the estimator `base_estimator` on
      bootstrap samples of the original data set, for different values of
      the regularization parameter for `base_estimator`. Variables that
      reliably get selected by the model in these bootstrap samples are
      considered to be stable variables.
  
      Parameters
      ----------
      base_estimator : object.
          The base estimator used for stability selection. The estimator
          must have either a ``feature_importances_`` or ``coef_``
          attribute after fitting.
  
      lambda_name : str.
          The name of the penalization parameter for the estimator
          `base_estimator`.
  
      lambda_grid : array-like.
          Grid of values of the penalization parameter to iterate over.
  
      n_bootstrap_iterations : integer.
          Number of bootstrap samples to create.
  
      sample_fraction : float, optional
          The fraction of samples to be used in each bootstrap sample.
          Should be between 0 and 1. If 1, all samples are used.
  
      threshold : float.
          Threshold defining the minimum cutoff value for the stability scores.
  
      bootstrap_func : str or callable fun (default=bootstrap_without_replacement)
          The function used to subsample the data. This parameter can be:
              - A string, which must be one of
                  - 'subsample': For subsampling without replacement.
                  - 'complementary_pairs': For complementary pairs subsampling [2]_ .
                  - 'stratified': For stratified bootstrapping in imbalanced
                     classification.
              - A function that takes y, and a random state
                as inputs and returns a list of sample indices in the range
                (0, len(y)-1). By default, indices are uniformly subsampled.
  
      bootstrap_threshold : string, float, optional default None
          The threshold value to use for feature selection. Features whose
          importance is greater or equal are kept while the others are
          discarded. If "median" (resp. "mean"), then the ``threshold`` value is
          the median (resp. the mean) of the feature importances. A scaling
          factor (e.g., "1.25*mean") may also be used. If None and if the
          estimator has a parameter penalty set to l1, either explicitly
          or implicitly (e.g, Lasso), the threshold used is 1e-5.
          Otherwise, "mean" is used by default.
  
      verbose : integer.
          Controls the verbosity: the higher, the more messages.
  
      n_jobs : int, default=1
          Number of jobs to run in parallel.
  
      pre_dispatch : int, or string, optional
          Controls the number of jobs that get dispatched during parallel
          execution. Reducing this number can be useful to avoid an
          explosion of memory consumption when more jobs get dispatched
          than CPUs can process. This parameter can be:
              - None, in which case all the jobs are immediately
                created and spawned. Use this for lightweight and
                fast-running jobs, to avoid delays due to on-demand
                spawning of the jobs
              - An int, giving the exact number of total jobs that are
                spawned
              - A string, giving an expression as a function of n_jobs,
                as in '2*n_jobs'
  
      random_state : int, RandomState instance or None, optional, default=None
          Pseudo random number generator state used for random uniform sampling
          from lists of possible values instead of scipy.stats distributions.
          If int, random_state is the seed used by the random number generator;
          If RandomState instance, random_state is the random number generator;
          If None, the random number generator is the RandomState instance used
          by `np.random`.
  
      Attributes
      ----------
      stability_scores_ : array, shape = [n_features, n_alphas]
          Array of stability scores for each feature for each value of the
          penalization parameter.
  
      References
      ----------
  
      .. [1] Meinshausen, N. and Buhlmann, P., 2010. Stability selection.
             Journal of the Royal Statistical Society: Series B
             (Statistical Methodology), 72(4), pp.417-473.
      .. [2] Shah, R.D. and Samworth, R.J., 2013. Variable selection with
             error control: another look at stability selection. Journal
             of the Royal Statistical Society: Series B (Statistical Methodology),
              75(1), pp.55-80.
      """
      def __init__(self, base_estimator=LogisticRegression(penalty='l1'), lambda_name='C',
                   lambda_grid=np.logspace(-5, -2, 25), n_bootstrap_iterations=100,
                   sample_fraction=0.5, threshold=0.6, bootstrap_func=bootstrap_without_replacement,
                   bootstrap_threshold=None, verbose=0, n_jobs=1, pre_dispatch='2*n_jobs',
                   random_state=None):
          self.base_estimator = base_estimator
          self.lambda_name = lambda_name
          self.lambda_grid = lambda_grid
          self.n_bootstrap_iterations = n_bootstrap_iterations
          self.sample_fraction = sample_fraction
          self.threshold = threshold
          self.bootstrap_func = bootstrap_func
          self.bootstrap_threshold = bootstrap_threshold
          self.verbose = verbose
          self.n_jobs = n_jobs
          self.pre_dispatch = pre_dispatch
          self.random_state = random_state
  
      def _validate_input(self):
          if not isinstance(self.n_bootstrap_iterations, int) or self.n_bootstrap_iterations <= 0:
              raise ValueError('n_bootstrap_iterations should be a positive integer, got %s' %
                               self.n_bootstrap_iterations)
  
          if not isinstance(self.sample_fraction, float) or not (0.0 < self.sample_fraction <= 1.0):
              raise ValueError('sample_fraction should be a float in (0, 1], got %s' % self.sample_fraction)
  
          if not isinstance(self.threshold, float) or not (0.0 < self.threshold <= 1.0):
              raise ValueError('threshold should be a float in (0, 1], got %s' % self.threshold)
  
          if self.lambda_name not in self.base_estimator.get_params().keys():
              raise ValueError('lambda_name is set to %s, but base_estimator %s '
                               'does not have a parameter '
                               'with that name' % (self.lambda_name,
                                                   self.base_estimator.__class__.__name__))
  
          if isinstance(self.bootstrap_func, str):
              if self.bootstrap_func not in BOOTSTRAP_FUNC_MAPPING.keys():
                  raise ValueError('bootstrap_func is set to %s, but must be one of '
                                   '%s or a callable' %
                                   (self.bootstrap_func, BOOTSTRAP_FUNC_MAPPING.keys()))
  
              self.bootstrap_func = BOOTSTRAP_FUNC_MAPPING[self.bootstrap_func]
          elif not callable(self.bootstrap_func):
              raise ValueError('bootstrap_func must be one of %s or a callable' %
                               BOOTSTRAP_FUNC_MAPPING.keys())
  
      def fit(self, X, y):
          """Fit the stability selection model on the given data.
  
          Parameters
          ----------
          X : {array-like, sparse matrix}, shape = [n_samples, n_features]
              The training input samples.
  
          y : array-like, shape = [n_samples]
              The target values.
          """
  
          self._validate_input()
  
          X, y = check_X_y(X, y, accept_sparse='csr')
  
          n_samples, n_variables = X.shape
          n_subsamples = np.floor(self.sample_fraction * n_samples).astype(int)
          n_lambdas = self.lambda_grid.shape[0]
  
          base_estimator = clone(self.base_estimator)
          random_state = check_random_state(self.random_state)
          stability_scores = np.zeros((n_variables, n_lambdas))
  
          for idx, lambda_value in enumerate(self.lambda_grid):
              if self.verbose > 0:
                  print("Fitting estimator for lambda = %.5f (%d / %d) on %d bootstrap samples" %
                        (lambda_value, idx + 1, n_lambdas, self.n_bootstrap_iterations))
  
              bootstrap_samples = _bootstrap_generator(self.n_bootstrap_iterations,
                                                       self.bootstrap_func, y,
                                                       n_subsamples, random_state=random_state)
  
              selected_variables = Parallel(
                  n_jobs=self.n_jobs, verbose=self.verbose,
                  pre_dispatch=self.pre_dispatch
              )(delayed(_fit_bootstrap_sample)(clone(base_estimator),
                                               X=X[safe_mask(X, subsample), :],
                                               y=y[subsample],
                                               lambda_name=self.lambda_name,
                                               lambda_value=lambda_value,
                                               threshold=self.bootstrap_threshold)
                for subsample in bootstrap_samples)
  
              stability_scores[:, idx] = np.vstack(selected_variables).mean(axis=0)
  
          self.stability_scores_ = stability_scores
          return self
  
      def get_support(self, indices=False, threshold=None):
          """Get a mask, or integer index, of the features selected
  
          Parameters
          ----------
          indices : boolean (default False)
              If True, the return value will be an array of integers,
              rather than a boolean mask.
  
          threshold: float.
              Threshold defining the minimum cutoff value for the
              stability scores.
  
          Returns
          -------
          support : array
              An index that selects the retained features from a feature vector.
              If `indices` is False, this is a boolean array of shape
              [# input features], in which an element is True iff its
              corresponding feature is selected for retention. If `indices` is
              True, this is an integer array of shape [# output features] whose
              values are indices into the input feature vector.
          """
  
          if threshold is not None and (not isinstance(threshold, float)
                                        or not (0.0 < threshold <= 1.0)):
              raise ValueError('threshold should be a float in (0, 1], '
                               'got %s' % self.threshold)
  
          cutoff = self.threshold if threshold is None else threshold
          mask = (self.stability_scores_.max(axis=1) > cutoff)
  
          return mask if not indices else np.where(mask)[0]
  
      def transform(self, X, threshold=None):
          """Reduce X to the selected features.
  
          Parameters
          ----------
          X : array of shape [n_samples, n_features]
              The input samples.
  
          threshold: float.
              Threshold defining the minimum cutoff value for the
              stability scores.
  
          Returns
          -------
          X_r : array of shape [n_samples, n_selected_features]
              The input samples with only the selected features.
          """
          X = check_array(X, accept_sparse='csr')
          mask = self.get_support(threshold=threshold)
  
          check_is_fitted(self, 'stability_scores_')
  
          if len(mask) != X.shape[1]:
              raise ValueError("X has a different shape than during fitting.")
  
          if not mask.any():
              warn("No features were selected: either the data is"
                   " too noisy or the selection test too strict.",
                   UserWarning)
              return np.empty(0).reshape((X.shape[0], 0))
  
          return X[:, safe_mask(X, mask)]
  