{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Minimum Fisher Regularization (MFR) tomography\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This page shows how to apply the MFR method to the specific case of tomography.\n", "\n", "The basic theory of the MFR is described in [this page](../../user/theory/mfr.ipynb).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example Tomography Problem\n", "\n", "As the example tomography problem, we consider the bolometer measurement and use the geometry matrix\n", "calculated at `cherab.core` [documentation](https://www.cherab.info/demonstrations/bolometry/geometry_matrix_with_raytransfer.html#bolometer-geometry-raytransfer)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pickle\n", "from pathlib import Path\n", "from pprint import pprint\n", "\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from matplotlib.cm import ScalarMappable\n", "from matplotlib.colors import AsinhNorm, ListedColormap, Normalize\n", "from matplotlib.ticker import (\n", " MultipleLocator,\n", " PercentFormatter,\n", ")\n", "from mpl_toolkits.axes_grid1 import ImageGrid\n", "\n", "from cherab.inversion import Mfr\n", "from cherab.inversion.data import get_sample_data\n", "from cherab.inversion.derivative import derivative_matrix\n", "\n", "plt.rcParams[\"figure.dpi\"] = 150\n", "\n", "# custom Red colormap extracted from \"RdBu_r\"\n", "cmap = plt.get_cmap(\"RdBu_r\")\n", "CMAP_RED = ListedColormap(cmap(np.linspace(0.5, 1.0, 256)))\n", "\n", "# path to the directory to store the MFR results\n", "store_dir = Path().cwd() / \"01-mfr\"\n", "store_dir.mkdir(exist_ok=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sample data can be retrieved by ``get_sample_data`` function with `\"bolo.npz\"` argument." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Load the sample data\n", "grid_data = get_sample_data(\"bolo.npz\")\n", "\n", "# Extract the data having several keys\n", "grid_centres = grid_data[\"grid_centres\"]\n", "voxel_map = grid_data[\"voxel_map\"].squeeze()\n", "mask = grid_data[\"mask\"].squeeze()\n", "gmat = grid_data[\"sensitivity_matrix\"]\n", "\n", "# Extract grid limits\n", "dr = grid_centres[1, 0, 0] - grid_centres[0, 0, 0]\n", "dz = grid_centres[0, 1, 1] - grid_centres[0, 0, 1]\n", "rmin, rmax = grid_centres[0, 0, 0] - 0.5 * dr, grid_centres[-1, 0, 0] + 0.5 * dr\n", "zmin, zmax = grid_centres[0, 0, 1] - 0.5 * dz, grid_centres[0, -1, 1] + 0.5 * dz\n", "\n", "print(f\"{grid_centres.shape = }\")\n", "print(f\"{voxel_map.shape = }\")\n", "print(f\"{mask.shape = }\")\n", "print(f\"{gmat.shape = }\")\n", "print(f\"gmat density: {np.count_nonzero(gmat) / gmat.size:.2%}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The geometry matrix density is 8.35% so this problem is known as the sparse problem as well.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Show the geometry matrix\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sensitivity_2d = np.full(voxel_map.shape, np.nan)\n", "sensitivity_2d[mask] = gmat.sum(axis=0)\n", "\n", "fig, ax = plt.subplots()\n", "image = ax.imshow(\n", " sensitivity_2d.T,\n", " origin=\"lower\",\n", " cmap=CMAP_RED,\n", " extent=(rmin, rmax, zmin, zmax),\n", ")\n", "cbar = plt.colorbar(image, pad=0.0)\n", "cbar.set_label(\"[m$^3$ sr]\")\n", "cbar.ax.ticklabel_format(style=\"sci\", axis=\"y\", useMathText=True)\n", "cbar.ax.yaxis.set_offset_position(\"left\")\n", "ax.set_title(\"Geomtery matrix\")\n", "ax.set_xlabel(\"$R$ [m]\")\n", "ax.set_ylabel(\"$Z$ [m]\")\n", "ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define phantom emission profile\n", "\n", "As a test emission profile, we define the following phantom emission profile which is used in the `cherab.core`\n", "[documentation](https://www.cherab.info/demonstrations/bolometry/inversion_with_raytransfer.html#bolometer-raytransfer-inversion) as well.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "PLASMA_AXIS = np.array([1.5, 1.5])\n", "LCFS_RADIUS = 1\n", "RING_RADIUS = 0.5\n", "\n", "RADIATION_PEAK = 1\n", "CENTRE_PEAK_WIDTH = 0.05\n", "RING_WIDTH = 0.025\n", "\n", "\n", "def emission_function_2d(rz_point: np.ndarray) -> float:\n", " direction = rz_point - PLASMA_AXIS\n", " bearing = np.arctan2(direction[1], direction[0])\n", "\n", " # calculate radius of coordinate from magnetic axis\n", " radius_from_axis = np.hypot(*direction)\n", " closest_ring_point = PLASMA_AXIS + (0.5 * direction / radius_from_axis)\n", " radius_from_ring = np.hypot(*(rz_point - closest_ring_point))\n", "\n", " # evaluate pedestal -> core function\n", " if radius_from_axis <= LCFS_RADIUS:\n", " central_radiatior = RADIATION_PEAK * np.exp(-(radius_from_axis**2) / CENTRE_PEAK_WIDTH)\n", "\n", " ring_radiator = (\n", " RADIATION_PEAK * np.cos(bearing) * np.exp(-(radius_from_ring**2) / RING_WIDTH)\n", " )\n", " ring_radiator = max(0, ring_radiator)\n", "\n", " return central_radiatior + ring_radiator\n", " else:\n", " return 0\n", "\n", "\n", "# Create a phantom vector\n", "phantom = np.zeros(gmat.shape[1])\n", "\n", "nr, nz = voxel_map.shape\n", "for ir, iz in np.ndindex(nr, nz):\n", " index = voxel_map[ir, iz]\n", " if index < 0:\n", " continue\n", " phantom[index] = emission_function_2d(grid_centres[ir, iz])\n", "\n", "# Create a 2D phantom\n", "phantom_2d = np.full(voxel_map.shape, np.nan)\n", "phantom_2d[mask] = phantom" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "image = ax.imshow(\n", " phantom_2d.T,\n", " origin=\"lower\",\n", " cmap=CMAP_RED,\n", " extent=(rmin, rmax, zmin, zmax),\n", ")\n", "cbar = plt.colorbar(image, pad=0.0)\n", "cbar.set_label(\"[W/m$^3$]\")\n", "ax.set_title(\"Phantom emissivity\")\n", "ax.set_xlabel(\"$R$ [m]\")\n", "ax.set_ylabel(\"$Z$ [m]\")\n", "ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compute the measurement data\n", "\n", "The bolometer measurement data should be calculated by the ray-tracing method, however, we generate\n", "the measurement data by the multiplication of the geometry matrix and the phantom profile.\n", "As the noise, we add the Gaussian noise with the standard deviation of 1% of the maximum value of the\n", "measurement data.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gmat /= 4.0 * np.pi # Divide by 4π steradians to use for measuring power in [W].\n", "\n", "data = gmat @ phantom\n", "\n", "rng = np.random.default_rng()\n", "data_w_noise = data + rng.normal(0, data.max() * 1.0e-2, data.size)\n", "data_w_noise = np.clip(data_w_noise, 0, None)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.bar(np.arange(data.size), data, label=\"w/o noise\")\n", "plt.plot(data_w_noise, \".\", c=\"C1\", label=\"w/ noise\")\n", "plt.legend()\n", "plt.xlabel(\"channel of bolometers\")\n", "plt.ylabel(\"Power [W]\")\n", "plt.xlim(-1, data.size)\n", "plt.ticklabel_format(style=\"sci\", axis=\"y\", useMathText=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Solve by the Laplacian regularization\n", "\n", "The Laplacian regularization (or Phillips regularization(Phillips 1962)) is one of the most popular regularization methods.\n", "However, this regularization tends to over-smooth the solution, which is unsuitable for the sparse problem.\n", "Let's look at this behavior by solving the tomography problem with the laplacian regularization.\n", "As the optimization method, we use the L-curve method." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from cherab.inversion import Lcurve, compute_svd\n", "from cherab.inversion.derivative import laplacian_matrix\n", "\n", "lmat = laplacian_matrix(voxel_map.shape, (dr, dz), mask=mask)\n", "s, u, basis = compute_svd(gmat, lmat.T @ lmat)\n", "\n", "lcurve = Lcurve(s, u, basis, data=data_w_noise)\n", "sol, status = lcurve.solve()\n", "lcurve.plot_L_curve()\n", "lcurve.plot_curvature()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comparison of the solution with the phantom reveals over-smoothing, resulting in lost details of the phantom." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vmin = min(sol.min(), phantom.min())\n", "vmax = max(sol.max(), phantom.max())\n", "\n", "fig = plt.figure(figsize=(10, 5))\n", "grids = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.05, cbar_mode=\"single\")\n", "norm = Normalize(vmin=vmin, vmax=vmax)\n", "for ax, value, label in zip(\n", " grids.axes_all, [phantom, sol], [\"Phantom\", \"Reconstruction\"], strict=True\n", "):\n", " image = np.full(voxel_map.shape, np.nan)\n", " image[mask] = value\n", "\n", " ax.imshow(\n", " image.T,\n", " origin=\"lower\",\n", " cmap=CMAP_RED,\n", " extent=(rmin, rmax, zmin, zmax),\n", " norm=norm,\n", " )\n", " ax.set_title(label)\n", " ax.set_xlabel(\"$R$ [m]\")\n", " ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "\n", "grids[0].set_ylabel(\"$Z$ [m]\")\n", "mappable = ScalarMappable(cmap=\"Reds\", norm=norm)\n", "cbar = plt.colorbar(mappable, cax=grids.cbar_axes[0])\n", "cbar.set_label(\"Emissivity [W/m$^3$]\")\n", "ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Quantitative evaluation\n", "\n", "Let's assess the solution quantitatively using the relative error, total power, and negative power values. Each value is calculated as follows:\n", "\n", "$$\n", "\\begin{gather}\n", "\\text{Relative error} = \\frac{\\left\\| \\mathbf{x}_\\mathrm{recon} - \\mathbf{x}_\\mathrm{phan} \\right\\|_2}{\\left\\| \\mathbf{x}_\\mathrm{phan} \\right\\|_2},\\\\\n", "\\text{Total power} = \\sum_i \\mathbf{x}_iv_i,\\\\\n", "\\text{Negative power} = \\sum_i \\min\\left\\{0, \\mathbf{x}_i\\right\\}v_i,\n", "\\end{gather}\n", "$$\n", "\n", "where $\\mathbf{x}_\\mathrm{recon}$ and $\\mathbf{x}_\\mathrm{phan}$ are the reconstructed and phantom profiles, respectively.\n", "$\\mathbf{x}_i$ represents the $i$-th element of the profile, and $v_i$ is the volume of the $i$-th\n", "voxel calculated with the Pappus theorem:\n", "\n", "$$\n", "v_i = \\Delta r \\Delta z \\cdot 2\\pi r_i,\n", "$$\n", "\n", "where $\\Delta r$ and $\\Delta z$ are the radial and vertical voxel sizes, respectively, and $r_i$ is the radial position of the $i$-th voxel.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "volumes = dr * dz * 2.0 * np.pi * grid_centres[:, :, 0]\n", "volumes = volumes[mask]\n", "\n", "print(\n", " \"The relative error :\",\n", " f\"{np.linalg.norm(sol - phantom) / np.linalg.norm(phantom):.2%}\",\n", ")\n", "print(\"--------------------------------------------\")\n", "print(f\"total power of phantom : {phantom @ volumes:.4g} W\")\n", "print(f\"total power of reconstruction: {sol @ volumes:.4g} W\")\n", "print(f\"total negative power of reconstruction: {sol[sol < 0] @ volumes[sol < 0]:.4g} W\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Solve by the MFR regularization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Prepare the derivative matrices\n", "\n", "Before performing the MFR tomography, we need to create derivative matrices. Here $\\mathbf{H}$ is defined as follows:\n", "\n", "$$\n", "\\mathbf{H} \\equiv \\mathbf{D}_r^\\mathsf{T} \\mathbf{W} \\mathbf{D}_r\n", " + \\mathbf{D}_z^\\mathsf{T} \\mathbf{W} \\mathbf{D}_z,\n", "$$\n", "\n", "where $\\mathbf{D}_r$ and $\\mathbf{D}_z$ are derivative matrices along the $r$ and $z$ coordinate direction, respectively.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dmat_r = derivative_matrix(voxel_map.shape, dr, axis=0, scheme=\"forward\", mask=mask)\n", "dmat_z = derivative_matrix(voxel_map.shape, dz, axis=1, scheme=\"forward\", mask=mask)\n", "\n", "dmat_pair = [(dmat_r, dmat_r), (dmat_z, dmat_z)]\n", "pprint(dmat_pair)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reconstruct\n" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "Let us do the MFR tomography. :obj:`cherab.inversion` offers the :obj:`.Mfr` class and simple :obj:`~.Mfr.solve` method.\n", "The :obj:`.store_regularizers` parameter of the :obj:`~.Mfr.solve` method stores the regularizer object :obj:`.Lcurve` instance at each MFR iteration." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "num_mfi = 4 # number of MFR iterations\n", "eps = 1.0e-6 # small positive number to avoid division by zero\n", "tol = 1.0e-3 # tolerance for the convergence criterion\n", "\n", "mfr = Mfr(gmat, dmat_pair, data=data_w_noise)\n", "sol, stats = mfr.solve(\n", " miter=num_mfi, eps=eps, tol=tol, store_regularizers=True, spinner=False, path=store_dir\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluate the tomographic reconstruction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compare the tomographic reconstruction with the phantom profile.\n", "We observe that the MFR tomography captures finer details of the phantom profile than the Laplacian regularization." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "nbsphinx-thumbnail" ] }, "outputs": [], "source": [ "vmin = min(sol.min(), phantom.min())\n", "vmax = max(sol.max(), phantom.max())\n", "\n", "fig = plt.figure(figsize=(10, 5))\n", "grids = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.05, cbar_mode=\"single\")\n", "norm = Normalize(vmin=vmin, vmax=vmax)\n", "for ax, value, label in zip(\n", " grids.axes_all, [phantom, sol], [\"Phantom\", \"Reconstruction\"], strict=True\n", "):\n", " image = np.full(voxel_map.shape, np.nan)\n", " image[mask] = value\n", "\n", " ax.imshow(\n", " image.T,\n", " origin=\"lower\",\n", " cmap=CMAP_RED,\n", " extent=(rmin, rmax, zmin, zmax),\n", " norm=norm,\n", " )\n", " ax.set_title(label)\n", " ax.set_xlabel(\"$R$ [m]\")\n", "\n", " ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "\n", "grids[0].set_ylabel(\"$Z$ [m]\")\n", "mappable = ScalarMappable(cmap=\"Reds\", norm=norm)\n", "cbar = plt.colorbar(mappable, cax=grids.cbar_axes[0])\n", "cbar.set_label(\"Emissivity [W/m$^3$]\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quantitative evaluations are shown as follows." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### L-curve and its curvature plot\n", "\n", "The L-curve of the last solution and its curvature are displayed below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lcurve = stats[\"regularizer\"]\n", "fig, axes = plt.subplots(1, 2, dpi=150, figsize=(9, 4), layout=\"constrained\")\n", "lcurve.plot_L_curve(fig, axes[0])\n", "lcurve.plot_curvature(fig, axes[1]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [solution bases](../../user/theory/inversion.ipynb#Series-expansion-of-the-solution) from 0-th to 18-th are shown below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot solution bases\n", "fig = plt.figure(figsize=(16, 10))\n", "grids = ImageGrid(\n", " fig,\n", " 111,\n", " nrows_ncols=(3, 6),\n", " axes_pad=0.0,\n", " cbar_mode=None,\n", ")\n", "\n", "for i, ax in enumerate(grids.axes_all):\n", " profile2d = np.full(voxel_map.shape, np.nan)\n", " profile2d[mask] = lcurve.basis[:, i]\n", "\n", " absolute = max(abs(lcurve.basis[:, i].min()), abs(lcurve.basis[:, i].max()))\n", " norm = AsinhNorm(vmin=-1 * absolute, vmax=absolute, linear_width=absolute * 1e-1)\n", "\n", " ax.imshow(\n", " profile2d.T,\n", " origin=\"lower\",\n", " cmap=\"RdBu_r\",\n", " extent=(rmin, rmax, zmin, zmax),\n", " norm=norm,\n", " )\n", " ax.set_xlabel(\"$R$ [m]\") if i >= 12 else None\n", " ax.set_ylabel(\"$Z$ [m]\") if i % 6 == 0 else None\n", " ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "\n", " ax.text(\n", " 0.5,\n", " 0.93,\n", " f\"basis {i}\",\n", " horizontalalignment=\"center\",\n", " verticalalignment=\"center\",\n", " transform=ax.transAxes,\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compare the measurement powers" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "back_calculated_measurements = gmat @ sol\n", "\n", "fig, (ax1, ax2) = plt.subplots(2, 1, layout=\"constrained\", figsize=(8, 6), sharex=True)\n", "\n", "# Plot the phantom and the back-calculated measurements\n", "ax1.bar(np.arange(data.size), data, label=\"Matrix-based measurements\")\n", "ax1.plot(back_calculated_measurements, \".\", color=\"C1\", label=\"Back-calculated from inversion\")\n", "ax1.legend()\n", "ax1.set_ylabel(\"Power [W]\")\n", "ax1.ticklabel_format(style=\"sci\", axis=\"y\", useMathText=True)\n", "ax1.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "ax1.set_xlim(-1, data.size)\n", "ax1.xaxis.set_major_locator(MultipleLocator(base=5))\n", "ax2.xaxis.set_minor_locator(MultipleLocator(base=1))\n", "\n", "# Plot the residuals between the measurements b and the back-calculated Tx (b - Tx)\n", "ax2.axhline(0, color=\"k\", linestyle=\"--\")\n", "ax2.bar(np.arange(data.size), data - back_calculated_measurements, color=\"C2\")\n", "ax2.set_xlabel(\"channel of bolometers\")\n", "ax2.set_ylabel(\"Residual [W]\")\n", "ax2.ticklabel_format(style=\"sci\", axis=\"y\", useMathText=True)\n", "ax2.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The measured power calculated by multiplying the geometry matrix by the emission vector and the back-calculated power calculated by multiplying the geometry matrix by the inverted emissivity are all in good agreement" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\n", " \"The relative error of power measurements is\",\n", " f\"{np.linalg.norm(data - back_calculated_measurements) / np.linalg.norm(data):.2%}\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\n", " \"The relative error :\",\n", " f\"{np.linalg.norm(sol - phantom) / np.linalg.norm(phantom):.2%}\",\n", ")\n", "print(\"--------------------------------------------\")\n", "print(f\"total power of phantom : {phantom @ volumes:.4g} W\")\n", "print(f\"total power of reconstruction: {sol @ volumes:.4g} W\")\n", "print(f\"total negative power of reconstruction: {sol[sol < 0] @ volumes[sol < 0]:.4g} W\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Iteration history" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As described in the [MFR definition section](../../user/theory/mfr.ipynb), MFR is the iterative method.\n", "In order to observe the convergence behavior, we will examine the iteration history." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reconstruction profiles history\n", "\n", "Firstly, let us see each iteration solution.\n", "To better understand the negative values, we can visualize the solutions using the arcsinh scale." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Stored regularizer files\n", "reg_files = list(store_dir.glob(\"regularizer_*.pickle\"))\n", "reg_files = sorted(reg_files, key=lambda x: int(x.stem.split(\"_\")[-1]))\n", "\n", "# Load each solution\n", "sols = []\n", "regs = []\n", "for file in reg_files:\n", " with open(file, \"rb\") as f:\n", " reg = pickle.load(f)\n", "\n", " sols.append(reg.solution(reg.lambda_opt))\n", " regs.append(reg)\n", "\n", "profiles = [phantom] + sols\n", "\n", "# set vmin and vmax for all solutions\n", "vmin = min(profile.min() for profile in profiles)\n", "vmax = max(profile.max() for profile in profiles)\n", "\n", "\n", "# Plot the solutions\n", "fig = plt.figure(figsize=(10, 5))\n", "grids = ImageGrid(\n", " fig,\n", " 111,\n", " nrows_ncols=(1, len(profiles)),\n", " axes_pad=0.05,\n", " cbar_mode=\"single\",\n", " cbar_location=\"right\",\n", ")\n", "\n", "absolute = max(abs(vmax), abs(vmin))\n", "linear_width = absolute * 0.1\n", "norm = AsinhNorm(linear_width=linear_width, vmin=-1 * absolute, vmax=absolute)\n", "i = 0\n", "for ax, profile in zip(grids.axes_all, profiles, strict=True):\n", " sol_2d = np.full(voxel_map.shape, np.nan)\n", " sol_2d[mask] = profile\n", "\n", " ax.imshow(\n", " sol_2d.T,\n", " origin=\"lower\",\n", " cmap=\"RdBu_r\",\n", " extent=(rmin, rmax, zmin, zmax),\n", " norm=norm,\n", " )\n", " if i == 0:\n", " ax.set_title(\"Phantom\")\n", " else:\n", " ax.set_title(f\"MFR iteration {i}\")\n", "\n", " ax.set_xlabel(\"$R$ [m]\")\n", " i += 1\n", "\n", " ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "\n", "grids[0].set_ylabel(\"$Z$ [m]\")\n", "mappable = ScalarMappable(norm=norm, cmap=\"RdBu_r\")\n", "cbar = plt.colorbar(mappable, cax=grids.cbar_axes[0])\n", "cbar.set_label(\"Emissivity [W/m$^3$]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Quantitative evaluation history\n", "\n", "The following plots display the quantitative evaluation changes over the course of the iteration." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "relative_errors = []\n", "total_powers = []\n", "negative_powers = []\n", "\n", "for sol in sols:\n", " relative_errors.append(np.linalg.norm(sol - phantom) / np.linalg.norm(phantom))\n", " total_powers.append(sol @ volumes)\n", " negative_powers.append(sol[sol < 0] @ volumes[sol < 0])\n", "\n", "# Append nan value to the last iteration\n", "relative_errors.append(np.nan)\n", "negative_powers.append(np.nan)\n", "\n", "# show each value as a bar plot\n", "fig, (ax1, ax2, ax3) = plt.subplots(3, 1, layout=\"constrained\", sharex=True, figsize=(4, 5))\n", "\n", "x = np.arange(1, len(relative_errors) + 1) # the label locations\n", "rects = ax1.bar(x[:4], total_powers, color=\"C0\", label=\"Reconstruction\")\n", "ax1.bar_label(rects, padding=-15, fmt=\"{:.3f}\", color=\"w\")\n", "rects = ax1.bar(x[-1], np.sum(phantom * volumes), color=\"C1\", label=\"Phantom\")\n", "ax1.bar_label(rects, padding=-15, fmt=\"{:.3f}\", color=\"w\")\n", "ax1.set_ylabel(\"Total power [W]\")\n", "ax1.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "ax1.ticklabel_format(style=\"sci\", axis=\"y\", useMathText=True)\n", "\n", "\n", "rects = ax2.bar(x, negative_powers, color=\"C2\")\n", "ax2.bar_label(rects, padding=3, fmt=\"{:.3f}\")\n", "ax2.set_ylim(ymin=-0.55)\n", "ax2.set_ylabel(\"Negative power [W]\")\n", "ax2.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "ax2.ticklabel_format(style=\"sci\", axis=\"y\", useMathText=True)\n", "\n", "rects = ax3.bar(x, relative_errors, color=\"C3\")\n", "ax3.bar_label(rects, padding=3, fmt=\"{:.1%}\")\n", "ax3.set_ylim(ymax=0.48)\n", "ax3.set_ylabel(\"Relative error\")\n", "ax3.set_xlabel(\"MFR iteration\")\n", "ax3.set_xticks(x)\n", "ax3.set_xticklabels(x.tolist()[:4] + [\"Phantom\"])\n", "ax3.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)\n", "ax3.yaxis.set_major_formatter(PercentFormatter(xmax=1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, the minimum relative error is achieved at the 2nd iteration.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution basis history\n", "\n", "Solution bases are also modified during the iteration and localized to the respective regions as the iteration progresses." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot solution bases\n", "fig = plt.figure(figsize=(16, 10))\n", "\n", "for j, reg in enumerate(regs):\n", " grids = ImageGrid(\n", " fig,\n", " (len(sols), 1, j + 1),\n", " nrows_ncols=(1, 6),\n", " axes_pad=0.0,\n", " cbar_mode=None,\n", " )\n", "\n", " grids[0].set_ylabel(f\"Iteration {j + 1}\")\n", "\n", " for i, ax in enumerate(grids.axes_all):\n", " profile2d = np.full(voxel_map.shape, np.nan)\n", " profile2d[mask] = reg.basis[:, i]\n", "\n", " absolute = max(abs(lcurve.basis[:, i].min()), abs(lcurve.basis[:, i].max()))\n", " norm = AsinhNorm(vmin=-1 * absolute, vmax=absolute, linear_width=absolute * 1e-1)\n", "\n", " ax.imshow(\n", " profile2d.T,\n", " origin=\"lower\",\n", " cmap=\"RdBu_r\",\n", " extent=(rmin, rmax, zmin, zmax),\n", " norm=norm,\n", " )\n", " ax.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)\n", "\n", " ax.text(\n", " 0.5,\n", " 0.93,\n", " f\"basis {i}\",\n", " horizontalalignment=\"center\",\n", " verticalalignment=\"center\",\n", " transform=ax.transAxes,\n", " )\n", "\n", "fig.subplots_adjust(hspace=0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As the iteration progresses, the solution bases expand towards the left and bottom edges. This is likely due to the absence of boundary restrictions from backward difference derivatives.\n", "\n", "To achieve a solution of zero at the boundary, anisotropic smoothing regularization is applied. This involves using numerical derivative matrices based on different coordinate systems, such as the polar coordinate system.\n", "\n", "If you'd like to see an example of anisotropic smoothing regularization using derivative matrices derived from different schemes, check out the [MFR tomography with anisotropic smoothing](02-mfr-aniso-smoothing.ipynb) example." ] } ], "metadata": { "kernelspec": { "display_name": "cherab-inv-dev", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 2 }