code: """ Principal Component Analysis
  """
  
  # Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
  #         Olivier Grisel <olivier.grisel@ensta.org>
  #         Mathieu Blondel <mathieu@mblondel.org>
  #         Denis A. Engemann <d.engemann@fz-juelich.de>
  #         Michael Eickenberg <michael.eickenberg@inria.fr>
  #
  # License: BSD 3 clause
  
  from math import log, sqrt
  
  import numpy as np
  from scipy import linalg
  from scipy.special import gammaln
  
  from ..base import BaseEstimator, TransformerMixin
  from ..utils import check_random_state, as_float_array
  from ..utils import check_array
  from ..utils.extmath import fast_dot, fast_logdet, randomized_svd
  from ..utils.validation import check_is_fitted
  
  
  def _assess_dimension_(spectrum, rank, n_samples, n_features):
      """Compute the likelihood of a rank ``rank`` dataset
  
      The dataset is assumed to be embedded in gaussian noise of shape(n,
      dimf) having spectrum ``spectrum``.
  
      Parameters
      ----------
      spectrum: array of shape (n)
          Data spectrum.
      rank: int
          Tested rank value.
      n_samples: int
          Number of samples.
      n_features: int
          Number of features.
  
      Returns
      -------
      ll: float,
          The log-likelihood
  
      Notes
      -----
      This implements the method of `Thomas P. Minka:
      Automatic Choice of Dimensionality for PCA. NIPS 2000: 598-604`
      """
      if rank > len(spectrum):
          raise ValueError("The tested rank cannot exceed the rank of the"
                           " dataset")
  
      pu = -rank * log(2.)
      for i in range(rank):
          pu += (gammaln((n_features - i) / 2.)
                 - log(np.pi) * (n_features - i) / 2.)
  
      pl = np.sum(np.log(spectrum[:rank]))
      pl = -pl * n_samples / 2.
  
      if rank == n_features:
          pv = 0
          v = 1
      else:
          v = np.sum(spectrum[rank:]) / (n_features - rank)
          pv = -np.log(v) * n_samples * (n_features - rank) / 2.
  
      m = n_features * rank - rank * (rank + 1.) / 2.
      pp = log(2. * np.pi) * (m + rank + 1.) / 2.
  
      pa = 0.
      spectrum_ = spectrum.copy()
      spectrum_[rank:n_features] = v
      for i in range(rank):
          for j in range(i + 1, len(spectrum)):
              pa += log((spectrum[i] - spectrum[j]) *
                        (1. / spectrum_[j] - 1. / spectrum_[i])) + log(n_samples)
  
      ll = pu + pl + pv + pp - pa / 2. - rank * log(n_samples) / 2.
  
      return ll
  
  
  def _infer_dimension_(spectrum, n_samples, n_features):
      """Infers the dimension of a dataset of shape (n_samples, n_features)
  
      The dataset is described by its spectrum `spectrum`.
      """
      n_spectrum = len(spectrum)
      ll = np.empty(n_spectrum)
      for rank in range(n_spectrum):
          ll[rank] = _assess_dimension_(spectrum, rank, n_samples, n_features)
      return ll.argmax()
  
  
  class PCA(BaseEstimator, TransformerMixin):
      """Principal component analysis (PCA)
  
      Linear dimensionality reduction using Singular Value Decomposition of the
      data and keeping only the most significant singular vectors to project the
      data to a lower dimensional space.
  
      This implementation uses the scipy.linalg implementation of the singular
      value decomposition. It only works for dense arrays and is not scalable to
      large dimensional data.
  
      The time complexity of this implementation is ``O(n ** 3)`` assuming
      n ~ n_samples ~ n_features.
  
      Read more in the :ref:`User Guide <PCA>`.
  
      Parameters
      ----------
      n_components : int, None or string
          Number of components to keep.
          if n_components is not set all components are kept::
  
              n_components == min(n_samples, n_features)
  
          if n_components == 'mle', Minka\'s MLE is used to guess the dimension
          if ``0 < n_components < 1``, select the number of components such that
          the amount of variance that needs to be explained is greater than the
          percentage specified by n_components
  
      copy : bool
          If False, data passed to fit are overwritten and running
          fit(X).transform(X) will not yield the expected results,
          use fit_transform(X) instead.
  
      whiten : bool, optional
          When True (False by default) the `components_` vectors are divided
          by n_samples times singular values to ensure uncorrelated outputs
          with unit component-wise variances.
  
          Whitening will remove some information from the transformed signal
          (the relative variance scales of the components) but can sometime
          improve the predictive accuracy of the downstream estimators by
          making there data respect some hard-wired assumptions.
  
      Attributes
      ----------
      components_ : array, [n_components, n_features]
          Components with maximum variance.
  
      explained_variance_ratio_ : array, [n_components]
          Percentage of variance explained by each of the selected components. \
          k is not set then all components are stored and the sum of explained \
          variances is equal to 1.0
  
      mean_ : array, [n_features]
          Per-feature empirical mean, estimated from the training set.
  
      n_components_ : int
          The estimated number of components. Relevant when n_components is set
          to 'mle' or a number between 0 and 1 to select using explained
          variance.
  
      noise_variance_ : float
          The estimated noise covariance following the Probabilistic PCA model
          from Tipping and Bishop 1999. See "Pattern Recognition and
          Machine Learning" by C. Bishop, 12.2.1 p. 574 or
          http://www.miketipping.com/papers/met-mppca.pdf. It is required to
          computed the estimated data covariance and score samples.
  
      Notes
      -----
      For n_components='mle', this class uses the method of `Thomas P. Minka:
      Automatic Choice of Dimensionality for PCA. NIPS 2000: 598-604`
  
      Implements the probabilistic PCA model from:
      M. Tipping and C. Bishop, Probabilistic Principal Component Analysis,
      Journal of the Royal Statistical Society, Series B, 61, Part 3, pp. 611-622
      via the score and score_samples methods.
      See http://www.miketipping.com/papers/met-mppca.pdf
  
      Due to implementation subtleties of the Singular Value Decomposition (SVD),
      which is used in this implementation, running fit twice on the same matrix
      can lead to principal components with signs flipped (change in direction).
      For this reason, it is important to always use the same estimator object to
      transform data in a consistent fashion.
  
      Examples
      --------
  
      >>> import numpy as np
      >>> from sklearn.decomposition import PCA
      >>> X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
      >>> pca = PCA(n_components=2)
      >>> pca.fit(X)
      PCA(copy=True, n_components=2, whiten=False)
      >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS
      [ 0.99244...  0.00755...]
  
      See also
      --------
      RandomizedPCA
      KernelPCA
      SparsePCA
      TruncatedSVD
      """
      def __init__(self, n_components=None, copy=True, whiten=False):
          self.n_components = n_components
          self.copy = copy
          self.whiten = whiten
  
      def fit(self, X, y=None):
          """Fit the model with X.
  
          Parameters
          ----------
          X: array-like, shape (n_samples, n_features)
              Training data, where n_samples in the number of samples
              and n_features is the number of features.
  
          Returns
          -------
          self : object
              Returns the instance itself.
          """
          self._fit(X)
          return self
  
      def fit_transform(self, X, y=None):
          """Fit the model with X and apply the dimensionality reduction on X.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              Training data, where n_samples is the number of samples
              and n_features is the number of features.
  
          Returns
          -------
          X_new : array-like, shape (n_samples, n_components)
  
          """
          U, S, V = self._fit(X)
          U = U[:, :self.n_components_]
  
          if self.whiten:
              # X_new = X * V / S * sqrt(n_samples) = U * sqrt(n_samples)
              U *= sqrt(X.shape[0])
          else:
              # X_new = X * V = U * S * V^T * V = U * S
              U *= S[:self.n_components_]
  
          return U
  
      def _fit(self, X):
          """Fit the model on X
  
          Parameters
          ----------
          X: array-like, shape (n_samples, n_features)
              Training vector, where n_samples in the number of samples and
              n_features is the number of features.
  
          Returns
          -------
          U, s, V : ndarrays
              The SVD of the input data, copied and centered when
              requested.
          """
          X = check_array(X)
          n_samples, n_features = X.shape
          X = as_float_array(X, copy=self.copy)
          # Center data
          self.mean_ = np.mean(X, axis=0)
          X -= self.mean_
          U, S, V = linalg.svd(X, full_matrices=False)
          explained_variance_ = (S ** 2) / n_samples
          explained_variance_ratio_ = (explained_variance_ /
                                       explained_variance_.sum())
  
          components_ = V
  
          n_components = self.n_components
          if n_components is None:
              n_components = n_features
          elif n_components == 'mle':
              if n_samples < n_features:
                  raise ValueError("n_components='mle' is only supported "
                                   "if n_samples >= n_features")
  
              n_components = _infer_dimension_(explained_variance_,
                                               n_samples, n_features)
          elif not 0 <= n_components <= n_features:
              raise ValueError("n_components=%r invalid for n_features=%d"
                               % (n_components, n_features))
  
          if 0 < n_components < 1.0:
              # number of components for which the cumulated explained variance
              # percentage is superior to the desired threshold
              ratio_cumsum = explained_variance_ratio_.cumsum()
              n_components = np.sum(ratio_cumsum < n_components) + 1
  
          # Compute noise covariance using Probabilistic PCA model
          # The sigma2 maximum likelihood (cf. eq. 12.46)
          if n_components < n_features:
              self.noise_variance_ = explained_variance_[n_components:].mean()
          else:
              self.noise_variance_ = 0.
  
          # store n_samples to revert whitening when getting covariance
          self.n_samples_ = n_samples
  
          self.components_ = components_[:n_components]
          self.explained_variance_ = explained_variance_[:n_components]
          explained_variance_ratio_ = explained_variance_ratio_[:n_components]
          self.explained_variance_ratio_ = explained_variance_ratio_
          self.n_components_ = n_components
  
          return (U, S, V)
  
      def get_covariance(self):
          """Compute data covariance with the generative model.
  
          ``cov = components_.T * S**2 * components_ + sigma2 * eye(n_features)``
          where  S**2 contains the explained variances.
  
          Returns
          -------
          cov : array, shape=(n_features, n_features)
              Estimated covariance of data.
          """
          components_ = self.components_
          exp_var = self.explained_variance_
          if self.whiten:
              components_ = components_ * np.sqrt(exp_var[:, np.newaxis])
          exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.)
          cov = np.dot(components_.T * exp_var_diff, components_)
          cov.flat[::len(cov) + 1] += self.noise_variance_  # modify diag inplace
          return cov
  
      def get_precision(self):
          """Compute data precision matrix with the generative model.
  
          Equals the inverse of the covariance but computed with
          the matrix inversion lemma for efficiency.
  
          Returns
          -------
          precision : array, shape=(n_features, n_features)
              Estimated precision of data.
          """
          n_features = self.components_.shape[1]
  
          # handle corner cases first
          if self.n_components_ == 0:
              return np.eye(n_features) / self.noise_variance_
          if self.n_components_ == n_features:
              return linalg.inv(self.get_covariance())
  
          # Get precision using matrix inversion lemma
          components_ = self.components_
          exp_var = self.explained_variance_
          exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.)
          precision = np.dot(components_, components_.T) / self.noise_variance_
          precision.flat[::len(precision) + 1] += 1. / exp_var_diff
          precision = np.dot(components_.T,
                             np.dot(linalg.inv(precision), components_))
          precision /= -(self.noise_variance_ ** 2)
          precision.flat[::len(precision) + 1] += 1. / self.noise_variance_
          return precision
  
      def transform(self, X):
          """Apply the dimensionality reduction on X.
  
          X is projected on the first principal components previous extracted
          from a training set.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              New data, where n_samples is the number of samples
              and n_features is the number of features.
  
          Returns
          -------
          X_new : array-like, shape (n_samples, n_components)
  
          """
          check_is_fitted(self, 'mean_')
  
          X = check_array(X)
          if self.mean_ is not None:
              X = X - self.mean_
          X_transformed = fast_dot(X, self.components_.T)
          if self.whiten:
              X_transformed /= np.sqrt(self.explained_variance_)
          return X_transformed
  
      def inverse_transform(self, X):
          """Transform data back to its original space, i.e.,
          return an input X_original whose transform would be X
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_components)
              New data, where n_samples is the number of samples
              and n_components is the number of components.
  
          Returns
          -------
          X_original array-like, shape (n_samples, n_features)
          """
          check_is_fitted(self, 'mean_')
  
          if self.whiten:
              return fast_dot(
                  X,
                  np.sqrt(self.explained_variance_[:, np.newaxis]) *
                  self.components_) + self.mean_
          else:
              return fast_dot(X, self.components_) + self.mean_
  
      def score_samples(self, X):
          """Return the log-likelihood of each sample
  
          See. "Pattern Recognition and Machine Learning"
          by C. Bishop, 12.2.1 p. 574
          or http://www.miketipping.com/papers/met-mppca.pdf
  
          Parameters
          ----------
          X: array, shape(n_samples, n_features)
              The data.
  
          Returns
          -------
          ll: array, shape (n_samples,)
              Log-likelihood of each sample under the current model
          """
          check_is_fitted(self, 'mean_')
  
          X = check_array(X)
          Xr = X - self.mean_
          n_features = X.shape[1]
          log_like = np.zeros(X.shape[0])
          precision = self.get_precision()
          log_like = -.5 * (Xr * (np.dot(Xr, precision))).sum(axis=1)
          log_like -= .5 * (n_features * log(2. * np.pi)
                            - fast_logdet(precision))
          return log_like
  
      def score(self, X, y=None):
          """Return the average log-likelihood of all samples
  
          See. "Pattern Recognition and Machine Learning"
          by C. Bishop, 12.2.1 p. 574
          or http://www.miketipping.com/papers/met-mppca.pdf
  
          Parameters
          ----------
          X: array, shape(n_samples, n_features)
              The data.
  
          Returns
          -------
          ll: float
              Average log-likelihood of the samples under the current model
          """
          return np.mean(self.score_samples(X))
  
  
  class RandomizedPCA(BaseEstimator, TransformerMixin):
      """Principal component analysis (PCA) using randomized SVD
  
      Linear dimensionality reduction using approximated Singular Value
      Decomposition of the data and keeping only the most significant
      singular vectors to project the data to a lower dimensional space.
  
      Read more in the :ref:`User Guide <RandomizedPCA>`.
  
      Parameters
      ----------
      n_components : int, optional
          Maximum number of components to keep. When not given or None, this
          is set to n_features (the second dimension of the training data).
  
      copy : bool
          If False, data passed to fit are overwritten and running
          fit(X).transform(X) will not yield the expected results,
          use fit_transform(X) instead.
  
      iterated_power : int, optional
          Number of iterations for the power method. 3 by default.
  
      whiten : bool, optional
          When True (False by default) the `components_` vectors are divided
          by the singular values to ensure uncorrelated outputs with unit
          component-wise variances.
  
          Whitening will remove some information from the transformed signal
          (the relative variance scales of the components) but can sometime
          improve the predictive accuracy of the downstream estimators by
          making their data respect some hard-wired assumptions.
  
      random_state : int or RandomState instance or None (default)
          Pseudo Random Number generator seed control. If None, use the
          numpy.random singleton.
  
      Attributes
      ----------
      components_ : array, [n_components, n_features]
          Components with maximum variance.
  
      explained_variance_ratio_ : array, [n_components]
          Percentage of variance explained by each of the selected components. \
          k is not set then all components are stored and the sum of explained \
          variances is equal to 1.0
  
      mean_ : array, [n_features]
          Per-feature empirical mean, estimated from the training set.
  
      Examples
      --------
      >>> import numpy as np
      >>> from sklearn.decomposition import RandomizedPCA
      >>> X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
      >>> pca = RandomizedPCA(n_components=2)
      >>> pca.fit(X)                 # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
      RandomizedPCA(copy=True, iterated_power=3, n_components=2,
             random_state=None, whiten=False)
      >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS
      [ 0.99244...  0.00755...]
  
      See also
      --------
      PCA
      TruncatedSVD
  
      References
      ----------
  
      .. [Halko2009] `Finding structure with randomness: Stochastic algorithms
        for constructing approximate matrix decompositions Halko, et al., 2009
        (arXiv:909)`
  
      .. [MRT] `A randomized algorithm for the decomposition of matrices
        Per-Gunnar Martinsson, Vladimir Rokhlin and Mark Tygert`
  
      """
  
      def __init__(self, n_components=None, copy=True, iterated_power=3,
                   whiten=False, random_state=None):
          self.n_components = n_components
          self.copy = copy
          self.iterated_power = iterated_power
          self.whiten = whiten
          self.random_state = random_state
  
      def fit(self, X, y=None):
          """Fit the model with X by extracting the first principal components.
  
          Parameters
          ----------
          X: array-like, shape (n_samples, n_features)
              Training data, where n_samples in the number of samples
              and n_features is the number of features.
  
          Returns
          -------
          self : object
              Returns the instance itself.
          """
          self._fit(check_array(X))
          return self
  
      def _fit(self, X):
          """Fit the model to the data X.
  
          Parameters
          ----------
          X: array-like, shape (n_samples, n_features)
              Training vector, where n_samples in the number of samples and
              n_features is the number of features.
  
          Returns
          -------
          X : ndarray, shape (n_samples, n_features)
              The input data, copied, centered and whitened when requested.
          """
          random_state = check_random_state(self.random_state)
          X = np.atleast_2d(as_float_array(X, copy=self.copy))
  
          n_samples = X.shape[0]
  
          # Center data
          self.mean_ = np.mean(X, axis=0)
          X -= self.mean_
          if self.n_components is None:
              n_components = X.shape[1]
          else:
              n_components = self.n_components
  
          U, S, V = randomized_svd(X, n_components,
                                   n_iter=self.iterated_power,
                                   random_state=random_state)
  
          self.explained_variance_ = exp_var = (S ** 2) / n_samples
          full_var = np.var(X, axis=0).sum()
          self.explained_variance_ratio_ = exp_var / full_var
  
          if self.whiten:
              self.components_ = V / S[:, np.newaxis] * sqrt(n_samples)
          else:
              self.components_ = V
  
          return X
  
      def transform(self, X, y=None):
          """Apply dimensionality reduction on X.
  
          X is projected on the first principal components previous extracted
          from a training set.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              New data, where n_samples in the number of samples
              and n_features is the number of features.
  
          Returns
          -------
          X_new : array-like, shape (n_samples, n_components)
  
          """
          check_is_fitted(self, 'mean_')
  
          X = check_array(X)
          if self.mean_ is not None:
              X = X - self.mean_
  
          X = fast_dot(X, self.components_.T)
          return X
  
      def fit_transform(self, X, y=None):
          """Fit the model with X and apply the dimensionality reduction on X.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_features)
              New data, where n_samples in the number of samples
              and n_features is the number of features.
  
          Returns
          -------
          X_new : array-like, shape (n_samples, n_components)
  
          """
          X = check_array(X)
          X = self._fit(X)
          return fast_dot(X, self.components_.T)
  
      def inverse_transform(self, X, y=None):
          """Transform data back to its original space.
  
          Returns an array X_original whose transform would be X.
  
          Parameters
          ----------
          X : array-like, shape (n_samples, n_components)
              New data, where n_samples in the number of samples
              and n_components is the number of components.
  
          Returns
          -------
          X_original array-like, shape (n_samples, n_features)
  
          Notes
          -----
          If whitening is enabled, inverse_transform does not compute the
          exact inverse operation of transform.
          """
          check_is_fitted(self, 'mean_')
  
          X_original = fast_dot(X, self.components_)
          if self.mean_ is not None:
              X_original = X_original + self.mean_
          return X_original
  
apis: numpy.sqrt numpy.sum numpy.eye numpy.empty numpy.zeros scipy.special.gammaln numpy.var scipy.linalg.inv numpy.log numpy.maximum numpy.dot numpy.mean scipy.linalg.svd