docs: def mvarica(x, var, cl=None, reducedim=0.99, optimize_var=False, backend=None,
              varfit='ensemble', random_state=None):
      """ Performs joint VAR model fitting and ICA source separation.
      
      This function implements the MVARICA procedure [1]_.
      
      Parameters
      ----------
      x : array-like, shape = [n_trials, n_channels, n_samples] or [n_channels, n_samples]
          data set
      var : :class:`~scot.var.VARBase`-like object
          Vector autoregressive model (VAR) object that is used for model fitting.
      cl : list of valid dict keys, optional
          Class labels associated with each trial.
      reducedim : {int, float, 'no_pca', None}, optional
          A number of less than 1 is interpreted as the fraction of variance that should remain in the data. All
          components that describe in total less than `1-reducedim` of the variance are removed by the PCA step.
          An integer number of 1 or greater is interpreted as the number of components to keep after applying PCA.
          If set to None, all PCA components are retained. If set to 'no_pca', the PCA step is skipped.
      optimize_var : bool, optional
          Whether to call automatic optimization of the VAR fitting routine.
      backend : dict-like, optional
          Specify backend to use. When set to None the backend configured in config.backend is used.
      varfit : string
          Determines how to calculate the residuals for source decomposition.
          'ensemble' (default) fits one model to the whole data set,
          'class' fits a new model for each class, and
          'trial' fits a new model for each individual trial.
          
      Returns
      -------
      result : class
          A class with the following attributes is returned:
              
          +---------------+----------------------------------------------------------+
          | mixing        | Source mixing matrix                                     |
          +---------------+----------------------------------------------------------+
          | unmixing      | Source unmixing matrix                                   |
          +---------------+----------------------------------------------------------+
          | residuals     | Residuals of the VAR model(s) in source space            |
          +---------------+----------------------------------------------------------+
          | var_residuals | Residuals of the VAR model(s) in EEG space (before ICA)  |
          +---------------+----------------------------------------------------------+
          | c             | Noise covariance of the VAR model(s) in source space     |
          +---------------+----------------------------------------------------------+
          | b             | VAR model coefficients (source space)                    |
          +---------------+----------------------------------------------------------+
          | a             | VAR model coefficients (EEG space)                       |
          +---------------+----------------------------------------------------------+
          
      Notes
      -----
      MVARICA is performed with the following steps:        
      1. Optional dimensionality reduction with PCA
      2. Fitting a VAR model tho the data
      3. Decomposing the VAR model residuals with ICA
      4. Correcting the VAR coefficients
          
      References
      ----------
      .. [1] G. Gomez-Herrero et al. "Measuring directional coupling between EEG sources", NeuroImage, 2008
      """
  
      x = atleast_3d(x)
      t, m, l = np.shape(x)
  
      if backend is None:
          backend = scotbackend
  
      # pre-transform the data with PCA
      if reducedim == 'no_pca':
          c = np.eye(m)
          d = np.eye(m)
          xpca = x
      else:
          c, d, xpca = backend['pca'](x, reducedim)
  
      if optimize_var:
          var.optimize(xpca)
  
      if varfit == 'trial':
          r = np.zeros(xpca.shape)
          for i in range(t):
              # fit MVAR model
              a = var.fit(xpca[i, :, :])
              # residuals
              r[i, :, :] = xpca[i, :, :] - var.predict(xpca[i, :, :])[0, :, :]
      elif varfit == 'class':
          r = np.zeros(xpca.shape)
          for i in np.unique(cl):
              mask = cl == i
              a = var.fit(xpca[mask, :, :])
              r[mask, :, :] = xpca[mask, :, :] - var.predict(xpca[mask, :, :])
      elif varfit == 'ensemble':
          # fit MVAR model
          a = var.fit(xpca)
          # residuals
          r = xpca - var.predict(xpca)
      else:
          raise ValueError('unknown VAR fitting mode: {}'.format(varfit))
  
      # run on residuals ICA to estimate volume conduction    
      mx, ux = backend['ica'](cat_trials(r), random_state=random_state)
  
      # driving process
      e = dot_special(ux.T, r)
  
      # correct AR coefficients
      b = a.copy()
      for k in range(0, a.p):
          b.coef[:, k::a.p] = mx.dot(a.coef[:, k::a.p].transpose()).dot(ux).transpose()
  
      # correct (un)mixing matrix estimatees
      mx = mx.dot(d)
      ux = c.dot(ux)
  
      class Result:
          unmixing = ux
          mixing = mx
          residuals = e
          var_residuals = r
          c = np.cov(cat_trials(e).T, rowvar=False)
  
      Result.b = b
      Result.a = a
      Result.xpca = xpca
          
      return Result