docs: def plot_fit(
      sampler,
      modelidx=0,
      label=None,
      sed=True,
      last_step=False,
      n_samples=100,
      confs=None,
      ML_info=False,
      figure=None,
      plotdata=None,
      plotresiduals=None,
      e_unit=None,
      e_range=None,
      e_npoints=100,
      threads=None,
      xlabel=None,
      ylabel=None,
      ulim_opts={},
      errorbar_opts={},
  ):
      """
      Plot data with fit confidence regions.
  
      Parameters
      ----------
      sampler : `emcee.EnsembleSampler`
          Sampler with a stored chain.
      modelidx : int, optional
          Model index to plot.
      label : str, optional
          Label for the title of the plot.
      sed : bool, optional
          Whether to plot SED or differential spectrum.
      last_step : bool, optional
          Whether to use only the samples of the last step in the run when
          showing either the model samples or the confidence intervals.
      n_samples : int, optional
          If not ``None``, number of sample models to plot. If ``None``,
          confidence bands will be plotted instead of samples. Default is 100.
      confs : list, optional
          List of confidence levels (in sigma) to use for generating the
          confidence intervals. Default is to plot sample models instead of
          confidence bands.
      ML_info : bool, optional
          Whether to plot information about the maximum likelihood parameters and
          the standard deviation of their distributions. Default is True.
      figure : `matplotlib.figure.Figure`, optional
          `matplotlib` figure to plot on. If omitted a new one will be generated.
      plotdata : bool, optional
          Wheter to plot data on top of model confidence intervals. Default is
          True if the physical types of the data and the model match.
      plotresiduals : bool, optional
          Wheter to plot the residuals with respect to the maximum likelihood
          model. Default is True if ``plotdata`` is True and either ``confs`` or
          ``n_samples`` are set.
      e_unit : `~astropy.units.Unit`, optional
          Units for the energy axis of the plot. The default is to use the units
          of the energy array of the observed data.
      e_range : list of `~astropy.units.Quantity`, length 2, optional
          Limits in energy for the computation of the model samples and ML model.
          Note that setting this parameter will mean that the samples for the
          model are recomputed and depending on the model speed might be quite
          slow.
      e_npoints : int, optional
          How many points to compute for the model samples and ML model if
          `e_range` is set.
      threads : int, optional
          How many parallel processing threads to use when computing the samples.
          Defaults to the number of available cores.
      xlabel : str, optional
          Label for the ``x`` axis of the plot.
      ylabel : str, optional
          Label for the ``y`` axis of the plot.
      ulim_opts : dict
          Option for upper-limit plotting. Available options are capsize (arrow
          width) and height_fraction (arrow length in fraction of flux value).
      errorbar_opts : dict
          Addtional options to pass to `matplotlib.plt.errorbar` for plotting the
          spectral flux points.
  
      """
      import matplotlib.pyplot as plt
  
      ML, MLp, MLerr, model_ML = find_ML(sampler, modelidx)
      infostr = "Maximum log probability: {0:.3g}\n".format(ML)
      infostr += "Maximum Likelihood values:\n"
      maxlen = np.max([len(ilabel) for ilabel in sampler.labels])
      vartemplate = "{{2:>{0}}}: {{0:>8.3g}} +/- {{1:<8.3g}}\n".format(maxlen)
      for p, v, ilabel in zip(MLp, MLerr, sampler.labels):
          infostr += vartemplate.format(p, v, ilabel)
  
      # log.info(infostr)
  
      data = sampler.data
  
      if e_range is None and not hasattr(sampler, "blobs"):
          e_range = data["energy"][[0, -1]] * np.array((1.0 / 3.0, 3.0))
  
      if plotdata is None and len(model_ML[0]) == len(data["energy"]):
          model_unit, _ = sed_conversion(model_ML[0], model_ML[1].unit, sed)
          data_unit, _ = sed_conversion(data["energy"], data["flux"].unit, sed)
          plotdata = model_unit.is_equivalent(data_unit)
      elif plotdata is None:
          plotdata = False
  
      if plotresiduals is None and plotdata and (confs is not None or n_samples):
          plotresiduals = True
  
      if confs is None and not n_samples and plotdata and not plotresiduals:
          # We actually only want to plot the data, so let's go there
          return plot_data(
              sampler.data,
              xlabel=xlabel,
              ylabel=ylabel,
              sed=sed,
              figure=figure,
              e_unit=e_unit,
              ulim_opts=ulim_opts,
              errorbar_opts=errorbar_opts,
          )
  
      if figure is None:
          f = plt.figure()
      else:
          f = figure
  
      if plotdata and plotresiduals:
          ax1 = plt.subplot2grid((4, 1), (0, 0), rowspan=3)
          ax2 = plt.subplot2grid((4, 1), (3, 0), sharex=ax1)
          for subp in [ax1, ax2]:
              f.add_subplot(subp)
      else:
          ax1 = f.add_subplot(111)
  
      if e_unit is None:
          e_unit = data["energy"].unit
  
      if confs is not None:
          plot_CI(
              ax1,
              sampler,
              modelidx,
              sed=sed,
              confs=confs,
              e_unit=e_unit,
              label=label,
              e_range=e_range,
              e_npoints=e_npoints,
              last_step=last_step,
              threads=threads,
          )
      elif n_samples:
          plot_samples(
              ax1,
              sampler,
              modelidx,
              sed=sed,
              n_samples=n_samples,
              e_unit=e_unit,
              label=label,
              e_range=e_range,
              e_npoints=e_npoints,
              last_step=last_step,
              threads=threads,
          )
      else:
          # plot only ML model
          _plot_MLmodel(ax1, sampler, modelidx, e_range, e_npoints, e_unit, sed)
  
      xlaxis = ax1
      if plotdata:
          _plot_data_to_ax(
              data,
              ax1,
              e_unit=e_unit,
              sed=sed,
              ylabel=ylabel,
              ulim_opts=ulim_opts,
              errorbar_opts=errorbar_opts,
          )
          if plotresiduals:
              _plot_residuals_to_ax(
                  data,
                  model_ML,
                  ax2,
                  e_unit=e_unit,
                  sed=sed,
                  errorbar_opts=errorbar_opts,
              )
              xlaxis = ax2
              for tl in ax1.get_xticklabels():
                  tl.set_visible(False)
          xmin = 10 ** np.floor(
              np.log10(
                  np.min(data["energy"] - data["energy_error_lo"])
                  .to(e_unit)
                  .value
              )
          )
          xmax = 10 ** np.ceil(
              np.log10(
                  np.max(data["energy"] + data["energy_error_hi"])
                  .to(e_unit)
                  .value
              )
          )
          ax1.set_xlim(xmin, xmax)
      else:
          ax1.set_xscale("log")
          ax1.set_yscale("log")
          if sed:
              ndecades = 10
          else:
              ndecades = 20
          # restrict y axis to ndecades to avoid autoscaling deep exponentials
          xmin, xmax, ymin, ymax = ax1.axis()
          ymin = max(ymin, ymax / 10 ** ndecades)
          ax1.set_ylim(bottom=ymin)
          # scale x axis to largest model_ML x point within ndecades decades of
          # maximum
          f_unit, sedf = sed_conversion(model_ML[0], model_ML[1].unit, sed)
          hi = np.where((model_ML[1] * sedf).to(f_unit).value > ymin)
          xmax = np.max(model_ML[0][hi])
          ax1.set_xlim(right=10 ** np.ceil(np.log10(xmax.to(e_unit).value)))
  
      if e_range:
          # ensure that xmin/xmax contains e_range
          xmin, xmax, ymin, ymax = ax1.axis()
          xmin = min(xmin, e_range[0].to(e_unit).value)
          xmax = max(xmax, e_range[1].to(e_unit).value)
          ax1.set_xlim(xmin, xmax)
  
      if ML_info and (confs is not None or n_samples):
          ax1.text(
              0.05,
              0.05,
              infostr,
              ha="left",
              va="bottom",
              transform=ax1.transAxes,
              family="monospace",
          )
  
      if label is not None:
          ax1.set_title(label)
  
      if xlabel is None:
          xlaxis.set_xlabel(
              "Energy [{0}]".format(e_unit.to_string("latex_inline"))
          )
      else:
          xlaxis.set_xlabel(xlabel)
  
      f.subplots_adjust(hspace=0)
  
      return f