docs: def _plot_residuals_to_ax(
      data_all, model_ML, ax, e_unit=u.eV, sed=True, errorbar_opts={}
  ):
      """Function to compute and plot residuals in units of the uncertainty"""
      if "group" not in data_all.keys():
          data_all["group"] = np.zeros(len(data_all))
  
      groups = np.unique(data_all["group"])
  
      MLf_unit, MLsedf = sed_conversion(model_ML[0], model_ML[1].unit, sed)
      MLene = model_ML[0].to(e_unit)
      MLflux = (model_ML[1] * MLsedf).to(MLf_unit)
  
      ax.axhline(0, color="k", lw=1, ls="--")
  
      interp = False
      if data_all["energy"].size != MLene.size or not np.allclose(
          data_all["energy"].value, MLene.value
      ):
          interp = True
          from scipy.interpolate import interp1d
  
          modelfunc = interp1d(MLene.value, MLflux.value, bounds_error=False)
  
      for g in groups:
          groupidx = np.where(data_all["group"] == g)
  
          data = data_all[groupidx]
  
          notul = ~data["ul"]
          df_unit, dsedf = sed_conversion(data["energy"], data["flux"].unit, sed)
          ene = data["energy"].to(e_unit)
          xerr = u.Quantity((data["energy_error_lo"], data["energy_error_hi"]))
          flux = (data["flux"] * dsedf).to(df_unit)
          dflux = (data["flux_error_lo"] + data["flux_error_hi"]) / 2.0
          dflux = (dflux * dsedf).to(df_unit)[notul]
  
          if interp:
              difference = flux[notul] - modelfunc(ene[notul]) * flux.unit
          else:
              difference = flux[notul] - MLflux[groupidx][notul]
  
          # wrap around color and marker cycles
          color = color_cycle[int(g) % len(color_cycle)]
          marker = marker_cycle[int(g) % len(marker_cycle)]
  
          opts = dict(
              zorder=100,
              marker=marker,
              ls="",
              elinewidth=2,
              capsize=0,
              mec=color,
              mew=0.1,
              ms=6,
              color=color,
          )
          opts.update(errorbar_opts)
  
          ax.errorbar(
              ene[notul].value,
              (difference / dflux).decompose().value,
              yerr=(dflux / dflux).decompose().value,
              xerr=xerr[:, notul].to(e_unit).value,
              **opts
          )
  
      from matplotlib.ticker import MaxNLocator
  
      ax.yaxis.set_major_locator(
          MaxNLocator(5, integer="True", prune="upper", symmetric=True)
      )
  
      ax.set_ylabel(r"$\Delta\sigma$")
      ax.set_xscale("log")