{
 "cells": [
  {
   "cell_type": "raw",
   "id": "0",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Invert a Synthetic Scene with MART\n",
    "==================================\n",
    "In this tutorial, we'll use the Multiplicative Algebraic Reconstruction Technique (MART) :cite:p:`Gordon1970` to invert the :func:`~ctis.scenes.gaussians` sample scene developed by Amy R. Winebarger.\n",
    "We will use a simple instrument with four projections to image this scene and then use the MART algorithm to reconstruct the original scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "1",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "import IPython.display\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.units as u\n",
    "import astropy.visualization\n",
    "import named_arrays as na\n",
    "import ctis"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "2",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Start by defining a grid of Doppler velocities on which to reconstruct the scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "velocity = na.linspace(-500, 500, axis=\"wavelength\", num=21) * u.km / u.s"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "4",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Define the rest wavelength for converting between velocity and wavelength."
   ]
  },
  {
   "cell_type": "code",
   "id": "5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "wavelength_rest = 171 * u.AA"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "6",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Now define a grid of positions on which to reconstruct the scene,"
   ]
  },
  {
   "cell_type": "code",
   "id": "7",
   "metadata": {},
   "source": [
    "position_scene = na.Cartesian2dVectorLinearSpace(\n",
    "    start=-10 * u.arcsec,\n",
    "    stop=10 * u.arcsec,\n",
    "    axis=na.Cartesian2dVectorArray(\"scene_x\", \"scene_y\"),\n",
    "    num=na.Cartesian2dVectorArray(64 + 1, 64 + 1),\n",
    ")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "8",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "and a grid of positions on the sensor representing the vertices of each pixel."
   ]
  },
  {
   "cell_type": "code",
   "id": "9",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "position_sensor = na.Cartesian2dVectorArray(\n",
    "    x=na.arange(0, 128 + 1, axis=\"sensor_x\") * u.pix,\n",
    "    y=na.arange(0, 64 + 1, axis=\"sensor_y\") * u.pix,\n",
    ")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "10",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Combine the 1D velocity grid and the 2D position grid into a single 3D grid for both the scene and sensor coordinates."
   ]
  },
  {
   "cell_type": "code",
   "id": "11",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "coordinates_scene = na.DopplerPositionalVectorArray.from_velocity(\n",
    "    velocity=velocity,\n",
    "    wavelength_rest=wavelength_rest,\n",
    "    position=position_scene,\n",
    ")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "12",
   "metadata": {},
   "source": [
    "coordinates_sensor = na.DopplerPositionalVectorArray.from_velocity(\n",
    "    velocity=velocity,\n",
    "    wavelength_rest=wavelength_rest,\n",
    "    position=position_sensor,\n",
    ")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "13",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Create a synthetic scene composed of spatial/spectral 3D Gaussians with various Doppler shifts."
   ]
  },
  {
   "cell_type": "code",
   "id": "14",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "scene = ctis.scenes.gaussians(coordinates_scene)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "15",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Add a small background equal to 1 percent of the maximum value of the scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "16",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "scene = scene + scene.max() / 100"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "17",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Display the scene as a false-color image."
   ]
  },
  {
   "cell_type": "code",
   "id": "18",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "with astropy.visualization.quantity_support():\n",
    "    fig, axs = plt.subplots(\n",
    "        ncols=2,\n",
    "        gridspec_kw=dict(width_ratios=[.9,.1]),\n",
    "        constrained_layout=True,\n",
    "    )\n",
    "    ax, cax = axs\n",
    "    colorbar = na.plt.rgbmesh(\n",
    "        C=scene,\n",
    "        axis_wavelength=\"wavelength\",\n",
    "        ax=ax,\n",
    "        vmin=0,\n",
    "        vmax=scene.outputs.max(),\n",
    "    )\n",
    "    na.plt.pcolormesh(\n",
    "        C=colorbar,\n",
    "        axis_rgb=\"wavelength\",\n",
    "        ax=cax,\n",
    "    )\n",
    "    ax.set_aspect(\"equal\")\n",
    "    ax.set_xlabel(f\"scene $x$ ({ax.get_xlabel()})\")\n",
    "    ax.set_ylabel(f\"scene $y$ ({ax.get_ylabel()})\")\n",
    "    cax.xaxis.set_ticks_position(\"top\")\n",
    "    cax.xaxis.set_label_position(\"top\")\n",
    "    cax.yaxis.tick_right()\n",
    "    cax.yaxis.set_label_position(\"right\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "19",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Compute the average spectrum of the scene"
   ]
  },
  {
   "cell_type": "code",
   "id": "20",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "spectrum = scene.outputs.mean((\"scene_x\", \"scene_y\"))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "21",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Plot the average spectrum of the scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "22",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "with astropy.visualization.quantity_support():\n",
    "    fig, ax = plt.subplots(constrained_layout=True)\n",
    "    ax2 = ax.twiny()\n",
    "    na.plt.stairs(\n",
    "        velocity,\n",
    "        spectrum,\n",
    "        ax=ax,\n",
    "    )\n",
    "    na.plt.stairs(\n",
    "        scene.inputs.wavelength,\n",
    "        spectrum,\n",
    "        ax=ax2\n",
    "    )\n",
    "    ax.set_xlabel(f\"Doppler velocity ({ax.get_xlabel()})\")\n",
    "    ax2.set_xlabel(f\"wavelength ({ax2.get_xlabel()})\")\n",
    "    ax.set_ylabel(f\"average radiance ({ax.get_ylabel()})\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "23",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Define the dispersion angles for our instrument.\n",
    "In this case we'll define four channels, each each separated by :math:`90^\\circ` degrees."
   ]
  },
  {
   "cell_type": "code",
   "id": "24",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "angle = na.linspace(0, 360, num=4, axis=\"channel\", endpoint=False) * u.deg + 5.64 * u.deg"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "25",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Define the magnitude of dispersion for our instrument in terms of Doppler velocity and then convert to wavelength units."
   ]
  },
  {
   "cell_type": "code",
   "id": "26",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "dispersion = 10 * u.km / u.s\n",
    "dispersion = dispersion.to(u.AA, equivalencies=u.doppler_optical(wavelength_rest))\n",
    "dispersion = (dispersion - wavelength_rest) / u.pix\n",
    "dispersion.to(u.mAA / u.pix)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "27",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Create an ideal CTIS using the dispersion magnitude and angles."
   ]
  },
  {
   "cell_type": "code",
   "id": "28",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "instrument = ctis.instruments.IdealInstrument(\n",
    "    area_effective=1 * u.cm ** 2,\n",
    "    timedelta_exposure=20 * u.s,\n",
    "    plate_scale=.4 * u.arcsec / u.pix,\n",
    "    dispersion=dispersion,\n",
    "    angle=angle,\n",
    "    wavelength_ref=wavelength_rest,\n",
    "    position_ref=na.Cartesian2dVectorArray(64, 32) * u.pix,\n",
    "    coordinates_scene=coordinates_scene,\n",
    "    coordinates_sensor=coordinates_sensor,\n",
    "    channel=\"dispersion angle = \" + angle.to_string_array(\"%03d\"),\n",
    "    axis_channel=\"channel\",\n",
    "    axis_wavelength=\"wavelength\",\n",
    "    axis_scene_xy=(\"scene_x\", \"scene_y\"),\n",
    "    axis_sensor_xy=(\"sensor_x\", \"sensor_y\"),\n",
    ")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "29",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Apply the forward model of this instrument to the scene to calculate the observed images."
   ]
  },
  {
   "cell_type": "code",
   "id": "30",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": "images = instrument.image(scene)",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "31",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Display the images as an animation, where each frame represents a different channel / dispersion direction."
   ]
  },
  {
   "cell_type": "code",
   "id": "32",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "with astropy.visualization.quantity_support():\n",
    "    fig, ax = plt.subplots(\n",
    "        constrained_layout=True,\n",
    "        figsize=(9.2, 4),\n",
    "    )\n",
    "    norm = plt.Normalize(\n",
    "        vmin=0,\n",
    "        vmax=images.outputs.value.ndarray.max(),\n",
    "    )\n",
    "    colorizer = plt.Colorizer(\n",
    "        cmap=\"gray\",\n",
    "        norm=norm,\n",
    "    )\n",
    "    ani = na.plt.pcolormovie(\n",
    "        instrument.channel,\n",
    "        images.inputs.position.x,\n",
    "        images.inputs.position.y,\n",
    "        C=images.outputs.value,\n",
    "        axis_time=\"channel\",\n",
    "        ax=ax,\n",
    "        kwargs_pcolormesh=dict(\n",
    "            colorizer=colorizer,\n",
    "        ),\n",
    "    )\n",
    "    plt.colorbar(\n",
    "        mappable=plt.cm.ScalarMappable(colorizer=colorizer),\n",
    "        ax=ax,\n",
    "        label=f\"signal ({images.outputs.unit:latex_inline})\",\n",
    "    )\n",
    "    ax.set_aspect(\"equal\")\n",
    "    ax.set_xlabel(f\"sensor $x$ ({images.inputs.position.x.unit})\")\n",
    "    ax.set_ylabel(f\"sensor $y$ ({images.inputs.position.y.unit})\")\n",
    "\n",
    "result = ani.to_jshtml(fps=2)\n",
    "result = IPython.display.HTML(result)\n",
    "\n",
    "plt.close(ani._fig)\n",
    "\n",
    "result"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "33",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Initialize the MART inversion algorithm with the instrument model. We'll also enable saving intermediate results so that we can visualize the behavior of the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "id": "34",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "mart = ctis.inverters.MartInverter(\n",
    "    instrument=instrument,\n",
    "    intermediate=True,\n",
    ")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "35",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Invert the images using our instance of MART and the initial guess."
   ]
  },
  {
   "cell_type": "code",
   "id": "36",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": "inversion = mart(images)",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "37",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Display the results as a false-color movie, where each frame represents subsequent iterations of the MART algorithm."
   ]
  },
  {
   "cell_type": "code",
   "id": "38",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "with astropy.visualization.quantity_support():\n",
    "    fig, axs = plt.subplots(\n",
    "        ncols=3,\n",
    "        gridspec_kw=dict(width_ratios=[.5, .5, .1]),\n",
    "        constrained_layout=True,\n",
    "        figsize=(10, 4.5),\n",
    "    )\n",
    "    ax1, ax2, cax = axs\n",
    "    ax2.set_yticklabels([])\n",
    "    na.plt.rgbmesh(\n",
    "        C=scene,\n",
    "        axis_wavelength=\"wavelength\",\n",
    "        ax=ax1,\n",
    "        vmin=0,\n",
    "        vmax=scene.outputs.max(),\n",
    "    )\n",
    "    label = \"iteration = \" + inversion.iteration.to_string_array(\"%d\") + \"\\n\"\n",
    "    name = r\"$\\langle \\chi^2 \\rangle$\"\n",
    "    label = label + f\"{name} = \" + inversion.mean_chi_squared.mean(instrument.axis_channel).to_string_array()\n",
    "    ani, colorbar = na.plt.rgbmovie(\n",
    "        label,\n",
    "        scene.inputs.wavelength,\n",
    "        scene.inputs.position.x,\n",
    "        scene.inputs.position.y,\n",
    "        C=inversion.solutions.outputs,\n",
    "        axis_time=inversion.inverter.axis_iteration,\n",
    "        axis_wavelength=\"wavelength\",\n",
    "        ax=ax2,\n",
    "        vmin=0,\n",
    "        vmax=scene.outputs.max(),\n",
    "    )\n",
    "    na.plt.pcolormesh(\n",
    "        C=colorbar,\n",
    "        axis_rgb=\"wavelength\",\n",
    "        ax=cax,\n",
    "    )\n",
    "    ax1.set_title(\"original\")\n",
    "    ax2.set_title(\"reconstructed\")\n",
    "    unit_x = scene.inputs.position.x.unit\n",
    "    unit_y = scene.inputs.position.y.unit\n",
    "    ax1.set_xlabel(f\"scene $x$ ({unit_x:latex_inline})\")\n",
    "    ax2.set_xlabel(f\"scene $x$ ({unit_x:latex_inline})\")\n",
    "    ax1.set_ylabel(f\"scene $y$ ({unit_y:latex_inline})\")\n",
    "    cax.xaxis.set_ticks_position(\"top\")\n",
    "    cax.xaxis.set_label_position(\"top\")\n",
    "    cax.yaxis.tick_right()\n",
    "    cax.yaxis.set_label_position(\"right\")\n",
    "\n",
    "result = ani.to_jshtml(fps=20)\n",
    "result = IPython.display.HTML(result)\n",
    "\n",
    "plt.close(ani._fig)\n",
    "\n",
    "result"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "39",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Isolate the solution array from the inversion result object."
   ]
  },
  {
   "cell_type": "code",
   "id": "40",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "solution = inversion.solution"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "41",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Compute the average spectrum of the reconstructed scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "42",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "spectrum_inverted = solution.outputs.mean((\"scene_x\", \"scene_y\"))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "43",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Plot the average spectrum of the original scene vs. the average spectrum of the reconstructed scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "44",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "with astropy.visualization.quantity_support():\n",
    "    fig, ax = plt.subplots(constrained_layout=True)\n",
    "    na.plt.stairs(\n",
    "        scene.inputs.wavelength,\n",
    "        spectrum,\n",
    "        ax=ax,\n",
    "        label=\"original\",\n",
    "    )\n",
    "    na.plt.stairs(\n",
    "        scene.inputs.wavelength,\n",
    "        spectrum_inverted,\n",
    "        ax=ax,\n",
    "        label=\"reconstructed\",\n",
    "    )\n",
    "    ax.set_xlabel(f\"wavelength ({ax.get_xlabel()})\")\n",
    "    ax2.set_xlabel(f\"wavelength ({ax2.get_xlabel()})\")\n",
    "    ax.set_ylabel(f\"average radiance ({ax.get_ylabel()})\")\n",
    "    ax.legend()"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "45",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Plot 2D histograms of the true vs. reconstructed value of the total radiance, median (Doppler shift), and interquartile range (Doppler width) for every pixel in the scene."
   ]
  },
  {
   "cell_type": "code",
   "id": "46",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "inversion.plot_moments(scene);"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "raw",
   "id": "47",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/x-rst",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Plot :math:`\\langle \\chi^2 \\rangle` and the signal-correlated residual as a function of iteration."
   ]
  },
  {
   "cell_type": "code",
   "id": "48",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "fig, ax = plt.subplots(\n",
    "    nrows=2,\n",
    "    sharex=True,\n",
    "    constrained_layout=True,\n",
    ")\n",
    "na.plt.plot(\n",
    "    inversion.iteration,\n",
    "    inversion.mean_chi_squared,\n",
    "    ax=ax[0],\n",
    "    axis=inversion.inverter.axis_iteration,\n",
    "    label=instrument.channel,\n",
    ")\n",
    "na.plt.plot(\n",
    "    inversion.iteration,\n",
    "    inversion.correlation_residual,\n",
    "    ax=ax[1],\n",
    "    axis=inversion.inverter.axis_iteration,\n",
    "    label=instrument.channel,\n",
    ")\n",
    "ax[0].set_ylabel(r\"$\\langle \\chi^2 \\rangle$\")\n",
    "ax[1].set_xlabel(\"iteration\")\n",
    "ax[1].set_ylabel(\"signal-correlated residual\")\n",
    "ax[0].set_yscale(\"log\")\n",
    "ax[0].legend();"
   ],
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}