docs: def pca(x, subtract_mean=False, normalize=False, sort_components=True,
          reducedim=None, algorithm=pca_eig):
      """Calculate principal component analysis (PCA).
      
      Parameters
      ----------
      x : ndarray, shape (trials, channels, samples) or (channels, samples)
          Input data.
      subtract_mean : bool, optional
          Subtract sample mean from x.
      normalize : bool, optional
          Normalize variances before applying PCA.
      sort_components : bool, optional
          Sort principal components in order of decreasing eigenvalues.
      reducedim : float or int or None, optional
          A value less than 1 is interpreted as the fraction of variance that
          should be retained in the data. All components that account for less
          than `1 - reducedim` of the variance are removed.
          An integer value of 1 or greater is interpreted as the number of
          (sorted) components to retain.
          If None, do not reduce dimensionality (i.e. keep all components).
      algorithm : func, optional
          Function to use for eigenvalue decomposition
          (:func:`pca_eig` or :func:`pca_svd`).
          
      Returns
      -------
      w : ndarray, shape (channels, components)
          PCA transformation matrix.
      v : ndarray, shape (components, channels)
          Inverse PCA transformation matrix.
      """
  
      x = np.asarray(x)
      if x.ndim == 3:
          x = cat_trials(x)
  
      if reducedim:
          sort_components = True
  
      if subtract_mean:
          x = x - np.mean(x, axis=1, keepdims=True)
  
      k, l = None, None
      if normalize:
          l = np.std(x, axis=1, ddof=1)
          k = np.diag(1.0 / l)
          l = np.diag(l)
          x = np.dot(k, x)
  
      w, latent = algorithm(x)
  
      # PCA is just a rotation, so inverse is equal to transpose
      v = w.T
  
      if normalize:
          w = np.dot(k, w)
          v = np.dot(v, l)
  
      latent /= sum(latent)
  
      if sort_components:
          order = np.argsort(latent)[::-1]
          w = w[:, order]
          v = v[order, :]
          latent = latent[order]
  
      if reducedim is not None:
          if reducedim < 1:
              selected = np.nonzero(np.cumsum(latent) < reducedim)[0]
              try:
                  selected = np.concatenate([selected, [selected[-1] + 1]])
              except IndexError:
                  selected = [0]
              if selected[-1] >= w.shape[1]:
                  selected = selected[0:-1]
              w = w[:, selected]
              v = v[selected, :]
          else:
              w = w[:, :reducedim]
              v = v[:reducedim, :]
  
      return w, v