docs: def optimize_order(self, data, min_p=1, max_p=None):
          """Determine optimal model order by minimizing the mean squared
          generalization error.
  
          Parameters
          ----------
          data : array, shape (n_trials, n_channels, n_samples)
              Epoched data set on which to optimize the model order. At least two
              trials are required.
          min_p : int
              Minimal model order to check.
          max_p : int
              Maximum model order to check
          """
          data = np.asarray(data)
          if data.shape[0] < 2:
              raise ValueError("At least two trials are required.")
  
          msge, prange = [], []
  
          par, func = parallel_loop(_get_msge_with_gradient, n_jobs=self.n_jobs,
                                    verbose=self.verbose)
          if self.n_jobs is None:
              npar = 1
          elif self.n_jobs < 0:
                  npar = 4  # is this a sane default?
          else:
              npar = self.n_jobs
  
          p = min_p
          while True:
              result = par(func(data, self.delta, self.xvschema, 1, p_)
                           for p_ in range(p, p + npar))
              j, k = zip(*result)
              prange.extend(range(p, p + npar))
              msge.extend(j)
              p += npar
              if max_p is None:
                  if len(msge) >= 2 and msge[-1] > msge[-2]:
                      break
              else:
                  if prange[-1] >= max_p:
                      i = prange.index(max_p) + 1
                      prange = prange[:i]
                      msge = msge[:i]
                      break
          self.p = prange[np.argmin(msge)]
          return zip(prange, msge)