docs: def plot_mixture(mixture, i=0, j=1, center_style=dict(s=0.15),
                   cmap='nipy_spectral', cutoff=0.0, ellipse_style=dict(alpha=0.3),
                   solid_edge=True, visualize_weights=False):
      '''Plot the (Gaussian) components of the ``mixture`` density as
      one-sigma ellipses in the ``(i,j)`` plane.
  
      :param center_style:
          If a non-empty ``dict``, plot mean value with the style passed to ``scatter``.
  
      :param cmap:
  
          The color map to which components are mapped in order to
          choose their face color. It is unaffected by the
          ``cutoff``. The meaning depends on ``visualize_weights``.
  
      :param cutoff:
          Ignore components whose weight is below the ``cut off``.
  
      :param ellipse_style:
          Passed on to define the properties of the ``Ellipse``.
  
      :param solid_edge:
          Draw the edge of the ellipse as solid opaque line.
  
      :param visualize_weights:
          Colorize the components according to their weights if ``True``.
          One can do `plt.colorbar()` after this function and the bar allows to read off the weights.
          If ``False``, coloring is based on the component index and the total number of components.
          This option makes it easier to track components by assigning them the same color in subsequent calls to this function.
  
      '''
      # imports inside the function because then "ImportError" is raised on
      # systems without 'matplotlib' only when 'plot_mixture' is called
      import numpy as np
      from matplotlib import pyplot as plt
      from matplotlib.patches import Ellipse
      from matplotlib.cm import get_cmap
  
      assert i >= 0 and j >= 0, 'Invalid submatrix specification (%d, %d)' % (i, j)
      assert i != j, 'Identical dimension given: i=j=%d' % i
      assert mixture.dim >= 2, '1D plot not supported'
  
      cmap = get_cmap(name=cmap)
  
      if visualize_weights:
          # colors according to weight
          renormalized_component_weights  = np.array(mixture.weights)
          colors = [cmap(k) for k in renormalized_component_weights]
      else:
          # colors according to index
          colors = [cmap(k) for k in np.linspace(0, _max_color, len(mixture.components))]
  
      mask = mixture.weights >= cutoff
  
      # plot component means
      means = np.array([c.mu for c in mixture.components])
      x_values = means.T[i]
      y_values = means.T[j]
  
      for k, w in enumerate(mixture.weights):
          # skip components by hand to retain consistent coloring
          if w < cutoff:
              continue
  
          cov = mixture.components[k].sigma
          submatrix = np.array([[cov[i,i], cov[i,j]], \
                                [cov[j,i], cov[j,j]]])
  
          # for idea, check
          # 'Combining error ellipses' by John E. Davis
          correlation = np.array([[1.0, cov[i,j] / np.sqrt(cov[i,i] * cov[j,j])], [0.0, 1.0]])
          correlation[1,0] = correlation[0,1]
  
          assert abs(correlation[0,1]) <= 1, 'Invalid component %d with correlation %g' % (k, correlation[0, 1])
  
          ew, ev = np.linalg.eigh(submatrix)
          assert ew.min() > 0, 'Nonpositive eigenvalue in component %d: %s' % (k, ew)
  
          # rotation angle of major axis with x-axis
          if submatrix[0,0] == submatrix[1,1]:
              theta = np.sign(submatrix[0,1]) * np.pi / 4.
          else:
              theta = 0.5 * np.arctan( 2 * submatrix[0,1] / (submatrix[1,1] - submatrix[0,0]))
  
          # put larger eigen value on y'-axis
          height = np.sqrt(ew.max())
          width = np.sqrt(ew.min())
  
          # but change orientation of coordinates if the other is larger
          if submatrix[0,0] > submatrix[1,1]:
              height = np.sqrt(ew.min())
              width = np.sqrt(ew.max())
  
          # change sign to rotate in right direction
          angle = -theta * 180 / np.pi
  
          # copy keywords but override some
          ellipse_style_clone = dict(ellipse_style)
  
          # overwrite facecolor
          ellipse_style_clone['facecolor'] = colors[k]
  
          ax = plt.gca()
  
          # need full width/height
          e = Ellipse(xy=(x_values[k], y_values[k]),
                                     width=2*width, height=2*height, angle=angle,
                                     **ellipse_style_clone)
          ax.add_patch(e)
  
          if solid_edge:
              ellipse_style_clone['facecolor'] = 'none'
              ellipse_style_clone['edgecolor'] = colors[k]
              ellipse_style_clone['alpha'] = 1
              ax.add_patch(Ellipse(xy=(x_values[k], y_values[k]),
                                         width=2*width, height=2*height, angle=angle,
                                         **ellipse_style_clone))
  
      if center_style:
          plt.scatter(x_values[mask], y_values[mask], **center_style)
  
      if visualize_weights:
          # to enable plt.colorbar()
          mappable = plt.gci()
          mappable.set_array(mixture.weights)
          mappable.set_cmap(cmap)