{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L-curve criterion\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This page shows how to use the L-curve criterion to find the regularization parameter in the case\n",
    "of the example ill-posed problem.\n",
    "\n",
    "The basic theory of the L-curve criterion is described in [this page](../../user/theory/lcurve.ipynb).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition of example inverse problem\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a famouse ill-posed linear equation, Fredholm integral equation is often used:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\int_a^b K(s, t)\\ x(t)\\ \\mathrm{d}t = b(s), \\quad c\\leq s\\leq d.\n",
    "\\label{eq:fredholm}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Here we think of the following situation as the above equation form:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& K(s, t) \\equiv (\\cos(s) + \\cos(t))\\left(\\frac{\\sin(u)}{u}\\right)^2,\\quad u \\equiv \\pi\\left(\\sin(s) + \\sin(t) \\right),\\\\\n",
    "& [a, b] = [c, d] \\equiv \\left[-\\frac{\\pi}{2}, \\frac{\\pi}{2} \\right].\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "And, the true solution $x_\\text{true}(t)$ is assumed as follows:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "x_\\text{true}(t) = 2.0 \\exp\\left[-6(t-0.8)^2 \\right] + \\exp\\left[-2(t+0.5)^2 \\right].\n",
    "\\end{equation}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.sparse import diags\n",
    "\n",
    "from cherab.inversion import Lcurve, compute_svd\n",
    "from cherab.inversion.tools import parse_scientific_notation\n",
    "\n",
    "plt.rcParams[\"figure.dpi\"] = 150"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, let us code the $K(s, t)$ and $x_\\text{true}(t)$ as a function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def kernel(s: np.ndarray, t: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Kernel of Fredholm integral equation of the first kind.\"\"\"\n",
    "    u = np.pi * (np.sin(s) + np.sin(t))\n",
    "    if u == 0:\n",
    "        return np.cos(s) + np.cos(t)\n",
    "    else:\n",
    "        return (np.cos(s) + np.cos(t)) * (np.sin(u) / u) ** 2\n",
    "\n",
    "\n",
    "def x_t_func(t: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"Define the function x_true(t)\"\"\"\n",
    "    return 2.0 * np.exp(-6.0 * (t - 0.8) ** 2) + np.exp(-2.0 * (t + 0.5) ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discretization of the equation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When discretizing the integral equation \\eqref{eq:fredholm} using the trapezoidal integral approximation,\n",
    "the following linear equation is obtained:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{K} \\mathbf{x} = \\mathbf{b},\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{K}\\in\\mathbb{R}^{M\\times N}$ is the discretized kernel matrix,\n",
    "$\\mathbf{x}\\in\\mathbb{R}^N$ is the discretized solution vector and\n",
    "$\\mathbf{b}\\in\\mathbb{R}^M$ is the discretized data vector.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$N$ and $M$ are the number of discretization points of $t$ and $s$, respectively.\n",
    "Here we set $N=M=64$ and generate points evenly spaced in $[-\\pi/2, \\pi/2]$.\n",
    "$x_\\mathrm{true}(t)$ discretized on these points yields the true solution vector $\\mathbf{x}_\\mathrm{true}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# discretize s, t\n",
    "s = np.linspace(-np.pi * 0.5, np.pi * 0.5, num=64, endpoint=True)\n",
    "t = np.linspace(-np.pi * 0.5, np.pi * 0.5, num=64, endpoint=True)\n",
    "\n",
    "# vectorize solution\n",
    "x_t = x_t_func(t)\n",
    "\n",
    "# discretize kernel\n",
    "k_mat = np.zeros((s.size, t.size))\n",
    "k_mat = np.array([[kernel(i, j) for j in t] for i in s])\n",
    "\n",
    "# trapezoidal rule\n",
    "k_mat[:, 0] *= 0.5\n",
    "k_mat[:, -1] *= 0.5\n",
    "k_mat *= t[1] - t[0]\n",
    "\n",
    "print(f\"{k_mat.shape = }\")\n",
    "print(f\"{x_t.shape = }\")\n",
    "print(f\"condition number of K is {np.linalg.cond(k_mat):.4g}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The given data $\\mathbf{b}$ is generated by adding white noise $\\mathbf{e}$ to the true\n",
    "data $\\bar{\\mathbf{b}} = \\mathbf{K}\\mathbf{x}_\\mathrm{true}$, that is,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{b} = \\bar{\\mathbf{b}} + \\mathbf{e}.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "The noise variance is set to $10^{-4}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_bar = k_mat @ x_t\n",
    "rng = np.random.default_rng()\n",
    "noise = rng.normal(0, 1.0e-4, b_bar.size)\n",
    "b = b_bar + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve the inverse problem\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution of the ill-posed linear equation is obtained with the regularization procedure:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{x}_\\lambda = \\left(\\mathbf{K}^\\mathsf{T}\\mathbf{K} + \\lambda\\mathbf{H}\\right)^{-1}\\mathbf{K}^\\mathsf{T}\\mathbf{b},\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{H}$ is the regularization matrix.\n",
    "Here we set $\\mathbf{H} = \\mathbf{D_2}^\\mathsf{T}\\mathbf{D_2}$, where $\\mathbf{D_2}$ is the second-order difference matrix.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmat = diags([1, -2, 1], [-1, 0, 1], shape=(t.size, t.size)).tocsr()\n",
    "print(f\"{dmat.shape = }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create lcurve solver object after calculating the singular value decomposition according\n",
    "to the [series expansion of solution](../../user/theory/inversion.ipynb#Series-expansion-of-the-solution).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s, u, basis = compute_svd(k_mat, dmat.T @ dmat)\n",
    "lcurve = Lcurve(s, u, basis, data=b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us solve the inverse problem.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol, status = lcurve.solve()\n",
    "print(status)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate the L-curve criterion\n",
    "\n",
    "Next we evaluate the solution obtained by the L-curve criterion\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot L-curve\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "fig, ax = lcurve.plot_L_curve(scatter_plot=7)\n",
    "ax.autoscale(axis=\"both\", tight=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The L-curve shown above is limited to the range of $\\lambda$ from $\\sigma_0^2$ to $\\sigma_{r}^2$ and\n",
    "it is enough to find the corner of the L-curve in this range.\n",
    "\n",
    "The below plot shows why it is enough by plotting points of $\\lambda = \\sigma_i^2$ on the L-curve,\n",
    "where $\\sigma_i$ is the $i$-th singular value and $i$ is indicated by the annotation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = lcurve.plot_L_curve(plot_lambda_opt=False)\n",
    "\n",
    "indices = list(range(0, 20)) + [lcurve.s.size - 1]\n",
    "sigmas = lcurve.s[indices]\n",
    "residuals = [lcurve.residual_norm(beta) for beta in sigmas**2]\n",
    "regularizations = [lcurve.regularization_norm(beta) for beta in sigmas**2]\n",
    "ax.scatter(residuals, regularizations, color=\"red\", marker=\".\")\n",
    "ax.legend([\"L-curve\", \"Points at $\\\\lambda = \\\\sigma_i^2$\"])\n",
    "for i, ind in enumerate(indices):\n",
    "    ax.annotate(f\"{ind}\", (residuals[i], regularizations[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"sigma_7^2 : {sigmas[7]**2:.4e}\")\n",
    "print(f\"sigma_8^2 : {sigmas[8]**2:.4e}\")\n",
    "print(f\"sigma_9^2 : {sigmas[9]**2:.4e}\")\n",
    "print(f\"lambda_opt: {lcurve.lambda_opt:.4e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot L-curve's curvature\n",
    "\n",
    "$\\lambda_\\mathrm{opt}$ is the regularization parameter that maximizes the curvature of the L-curve.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = lcurve.plot_curvature()\n",
    "lambda_opt_text = parse_scientific_notation(f\"{lcurve.lambda_opt:.3e}\")\n",
    "ax.set_title(f\"$\\\\lambda_\\\\mathrm{{opt}} = {lambda_opt_text}$\")\n",
    "ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare $\\mathbf{x}_\\lambda$ with $\\mathbf{x}_\\mathrm{true}$\n",
    "\n",
    "Let us compare solutions at different regularization parameters $\\lambda=10^{-9}$,\n",
    "$\\lambda_\\text{opt}$, $10^3$ with the true solution $\\mathbf{x}_\\mathrm{true}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "lambdas = [1.0e-9, lcurve.lambda_opt, 1.0e3]\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(10, 4), sharey=True, layout=\"constrained\")\n",
    "\n",
    "for ax, beta in zip(axes, lambdas, strict=False):\n",
    "    ax.plot(t, x_t, \"--\", label=\"$\\\\mathbf{x}_\\\\mathrm{true}$\")\n",
    "    ax.plot(t, lcurve.solution(beta=beta), label=\"$\\\\mathbf{x}_\\\\lambda$\")\n",
    "\n",
    "    ax.set_xlim(t.min(), t.max())\n",
    "    ax.set_ylim(0, x_t.max() * 1.1)\n",
    "    ax.set_xlabel(\"$t$\")\n",
    "    parsed_lambda = parse_scientific_notation(f\"{beta:.2e}\")\n",
    "    ax.set_title(f\"$\\\\lambda = {parsed_lambda}$\")\n",
    "    ax.tick_params(direction=\"in\", labelsize=10, which=\"both\", top=True, right=True)\n",
    "    ax.legend(loc=\"upper left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the solution at $\\lambda < \\lambda_\\mathrm{opt}$ is perturbed by noise, while\n",
    "the solution at $\\lambda > \\lambda_\\mathrm{opt}$ is smoothed too much.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot norms and curvature as a function of $\\lambda$\n",
    "\n",
    "Let us plot the residual norm $\\sqrt{\\rho}$ and the regularization norm $\\sqrt{\\eta}$ as a function of $\\lambda$.\n",
    "Additionally, we plot the curvature of the L-curve as a function of $\\lambda$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "fig.subplots_adjust(right=0.85)\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax3 = ax1.twinx()\n",
    "\n",
    "ax3.spines.right.set_position((\"axes\", 1.2))\n",
    "\n",
    "# calculation of the values\n",
    "lambdas = np.logspace(-10, 0, num=500)\n",
    "rhos = [lcurve.residual_norm(beta) for beta in lambdas]\n",
    "etas = [lcurve.regularization_norm(beta) for beta in lambdas]\n",
    "kappa = [lcurve.curvature(beta) for beta in lambdas]\n",
    "\n",
    "# plot lines\n",
    "(p1,) = ax1.loglog(lambdas, rhos, color=\"C0\")\n",
    "(p2,) = ax2.loglog(lambdas, etas, color=\"C1\")\n",
    "(p3,) = ax3.semilogx(lambdas, kappa, color=\"C2\")\n",
    "\n",
    "# set axes properties\n",
    "ax1.set(\n",
    "    xlim=(lambdas[0], lambdas[-1]),\n",
    "    xlabel=\"Regularization parameter $\\\\lambda$\",\n",
    "    ylabel=\"Residual norm $\\\\sqrt{\\\\rho}$\",\n",
    ")\n",
    "ax2.set(ylabel=\"Regularization norm $\\\\sqrt{\\\\eta}$\")\n",
    "ax3.set(ylabel=\"curvature of L-curve\")\n",
    "\n",
    "ax1.yaxis.label.set_color(p1.get_color())\n",
    "ax2.yaxis.label.set_color(p2.get_color())\n",
    "ax3.yaxis.label.set_color(p3.get_color())\n",
    "\n",
    "ax1.tick_params(axis=\"x\", which=\"both\", direction=\"in\", top=True)\n",
    "ax1.tick_params(axis=\"y\", which=\"both\", direction=\"in\", colors=p1.get_color())\n",
    "ax2.tick_params(axis=\"y\", which=\"both\", direction=\"in\", colors=p2.get_color())\n",
    "ax3.tick_params(axis=\"y\", which=\"both\", direction=\"in\", colors=p3.get_color())\n",
    "\n",
    "ax3.spines[\"left\"].set_color(p1.get_color())\n",
    "ax2.spines[\"right\"].set_color(p2.get_color())\n",
    "ax3.spines[\"right\"].set_color(p3.get_color())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\sqrt{\\rho}$ is monotonically increasing with $\\lambda$, while $\\sqrt{\\eta}$ is monotonically decreasing with $\\lambda$.\n",
    "This behavior is consistent with the theory of the L-curve criterion.\n",
    "\n",
    "The curvature of the L-curve is maximized at the center region where both are flat.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the relative error\n",
    "\n",
    "The relative error between the solution $\\mathbf{x}_\\lambda$ and the true solution\n",
    "$\\mathbf{x}_\\mathrm{true}$ is defined as follows:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\epsilon_\\mathrm{rel} = \\frac{\\|\\mathbf{x}_\\lambda - \\mathbf{x}_\\mathrm{true}\\|_2}{\\|\\mathbf{x}_\\mathrm{true}\\|_2}.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Let us seek the minimum $\\epsilon_\\mathrm{rel}$ as a function of $\\lambda$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize_scalar\n",
    "\n",
    "X_T_NORM = np.linalg.norm(x_t, axis=0)\n",
    "\n",
    "\n",
    "def relative_error(\n",
    "    log_lambda: float, x_t: np.ndarray = x_t, x_t_norm: float = X_T_NORM, lcurve: Lcurve = lcurve\n",
    ") -> float:\n",
    "    \"\"\"Calculate relative error.\"\"\"\n",
    "    beta = 10**log_lambda\n",
    "    sol = lcurve.solution(beta=beta)\n",
    "    return np.linalg.norm(x_t - sol, axis=0) / x_t_norm\n",
    "\n",
    "\n",
    "# minimize relative error\n",
    "bounds = -10, -1\n",
    "res = minimize_scalar(\n",
    "    relative_error,\n",
    "    bounds=bounds,\n",
    "    method=\"bounded\",\n",
    "    args=(x_t, X_T_NORM, lcurve),\n",
    "    options={\"xatol\": 1.0e-10, \"maxiter\": 1000},\n",
    ")\n",
    "\n",
    "# obtain minimum relative error and lambda\n",
    "error_min = res.fun\n",
    "lambda_min = 10**res.x\n",
    "\n",
    "print(f\"minimum relative error: {error_min:.2%} at lambda = {lambda_min:.4g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot the relative error and curvature as a function of $\\lambda$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# set regularization parameters\n",
    "num = 500\n",
    "lambdas = np.logspace(*bounds, num=num)\n",
    "\n",
    "# calculate errors and curvatures\n",
    "errors = np.asarray([relative_error(log_lambda) for log_lambda in np.linspace(*bounds, num=num)])\n",
    "kappa = np.asarray([lcurve.curvature(beta) for beta in lambdas])\n",
    "\n",
    "# create figure\n",
    "fig, ax1 = plt.subplots()\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "# plot errors and curvatures\n",
    "(p1,) = ax1.loglog(lambdas, errors, color=\"C0\")\n",
    "(p2,) = ax2.semilogx(lambdas, kappa, color=\"C1\")\n",
    "\n",
    "# plot minimum error vertical line and point\n",
    "ax1.axvline(lambda_min, color=\"r\", linestyle=\"--\", linewidth=0.75)\n",
    "ax1.scatter(lambda_min, error_min, color=\"r\", marker=\"o\", s=10, zorder=2)\n",
    "ax1.text(\n",
    "    lambda_min,\n",
    "    1.5e-2,\n",
    "    \"$\\\\lambda_\\\\mathrm{min}$\",\n",
    "    color=\"r\",\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\",\n",
    ")\n",
    "\n",
    "# plot maximum curvature vertical line and point\n",
    "assert lcurve.lambda_opt is not None\n",
    "ax1.axvline(lcurve.lambda_opt, color=\"g\", linestyle=\"--\", linewidth=0.75)\n",
    "ax2.scatter(\n",
    "    lcurve.lambda_opt, lcurve.curvature(lcurve.lambda_opt), color=\"g\", marker=\"o\", s=10, zorder=2\n",
    ")\n",
    "ax1.text(\n",
    "    lcurve.lambda_opt,\n",
    "    1.5e-2,\n",
    "    \"$\\\\lambda_\\\\mathrm{opt}$\",\n",
    "    color=\"g\",\n",
    "    horizontalalignment=\"left\",\n",
    "    verticalalignment=\"center\",\n",
    ")\n",
    "\n",
    "# set axes\n",
    "ax1.set(\n",
    "    xlim=(lambdas[0], lambdas[-1]),\n",
    "    ylim=(0.01, 1),\n",
    "    xlabel=\"$\\\\lambda$\",\n",
    "    ylabel=\"Relative error $\\\\epsilon_\\\\mathrm{rel}$\",\n",
    ")\n",
    "ax2.set(ylabel=\"curvature of L-curve\")\n",
    "\n",
    "ax1.yaxis.label.set_color(p1.get_color())\n",
    "ax2.yaxis.label.set_color(p2.get_color())\n",
    "\n",
    "ax1.tick_params(axis=\"x\", which=\"both\", direction=\"in\", top=True)\n",
    "ax1.tick_params(axis=\"y\", which=\"both\", direction=\"in\", colors=p1.get_color())\n",
    "ax2.tick_params(axis=\"y\", which=\"both\", direction=\"in\", colors=p2.get_color())\n",
    "\n",
    "ax2.spines[\"left\"].set_color(p1.get_color())\n",
    "ax2.spines[\"right\"].set_color(p2.get_color())\n",
    "\n",
    "error_opt = relative_error(np.log10(lcurve.lambda_opt))\n",
    "ax1.set_title(\n",
    "    f\"$\\\\epsilon_\\\\mathrm{{rel}}(\\\\lambda_\\\\mathrm{{opt}})$ = {error_opt:.2%}, \"\n",
    "    + f\"$\\\\epsilon_\\\\mathrm{{rel}}(\\\\lambda_\\\\mathrm{{min}}) = ${error_min:.2%}\"\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the solution at $\\lambda_\\mathrm{opt}$ with $\\lambda_\\mathrm{min}$\n",
    "\n",
    "$\\lambda_\\mathrm{min}$ is the regularization parameter that minimizes the relative error.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(t, x_t, \"k--\", label=\"$\\\\mathbf{x}_\\\\mathrm{true}$\", linewidth=1.0)\n",
    "ax.plot(t, lcurve.solution(lambda_min), label=\"$\\\\lambda_\\\\mathrm{min}$\")\n",
    "ax.plot(t, sol, label=\"$\\\\lambda_\\\\mathrm{opt}$\")\n",
    "ax.set_xlabel(\"$t$\")\n",
    "ax.set_xlim(t.min(), t.max())\n",
    "ax.set_ylim(0, x_t.max() * 1.1)\n",
    "ax.legend()\n",
    "ax.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True, right=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete Picard Plot\n",
    "\n",
    "**Discrete Picard Condition**<cite data-footcite=\"Correia2009-xk\"></cite>\n",
    "\n",
    "> The data vector $\\mathbf{b}$ satisfies the **discrete Picard condition (DPC)** if the data space\n",
    "> coefficients $|\\mathbf{u}_i^\\mathsf{T}\\mathbf{b}|$ decay faster than the singular values $\\sigma_i$.\n",
    "\n",
    "In ill-posed problems, we find that the DPC holds initially and then fails at some point $i_\\mathrm{DPC}$,\n",
    "where the data become dominated by errors (noise).\n",
    "If this is the case, and if the regularization parameter is accurately selected,\n",
    "then the regularized solution should provide a valid solution.\n",
    "Examining $i_\\mathrm{DPC}$ provides a method of characterizing the ill-posedness of the problem.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ub = np.abs(lcurve.U.T @ b)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax2 = ax1.twinx()\n",
    "\n",
    "# DP plot\n",
    "ax1.semilogy(lcurve.s, \".-\", label=\"$\\\\sigma_i$\")\n",
    "ax1.semilogy(ub, \"o\", label=\"$|\\\\mathbf{u}_i^\\\\mathsf{T}\\\\mathbf{b}|$\")\n",
    "ax1.semilogy(ub / lcurve.s, \"s\", label=\"$|\\\\mathbf{u}_i^\\\\mathsf{T}\\\\mathbf{b}|/\\\\sigma_i$\")\n",
    "ax1.axhline(1e-4, color=\"k\", linestyle=\"--\", zorder=0)\n",
    "ax1.axvline(10, color=\"k\", linestyle=\"--\", zorder=0)\n",
    "ax1.set_xlabel(\"$i$\")\n",
    "ax1.legend(loc=\"lower left\")\n",
    "ax1.tick_params(axis=\"both\", which=\"both\", direction=\"in\", top=True)\n",
    "\n",
    "# filter plot\n",
    "assert lcurve.lambda_opt is not None\n",
    "(p2,) = ax2.semilogy(\n",
    "    lcurve.filter(lcurve.lambda_opt), \".-\", color=\"C4\", label=\"$\\\\mathbf{f}_\\\\mathrm{opt}$\"\n",
    ")\n",
    "ax2.set_ylabel(\"$f_{\\\\lambda, i}$\", color=\"C4\")\n",
    "\n",
    "# color axes2\n",
    "ax2.yaxis.label.set_color(p2.get_color())\n",
    "ax2.tick_params(axis=\"y\", which=\"both\", direction=\"in\", colors=p2.get_color())\n",
    "ax2.spines[\"right\"].set_color(p2.get_color())\n",
    "\n",
    "ax1.set_xlim(0, 25);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vertical dashed line in the above figure marks the biginning of $|\\mathbf{u}_i^\\mathsf{T}\\mathbf{b}| < \\sigma_i$,\n",
    "and the horizontal dashed line represents the noise level.\n",
    "DPC is satisfied for $i < i_\\mathrm{DPC} \\simeq 10$.<br>\n",
    "We confirm that the filter factor $f_{\\lambda_\\mathrm{opt}, i}$, the $\\lambda_\\mathrm{opt}$ of which\n",
    "is selected by the L-curve criterion, starts to decrease around $i_\\mathrm{DPC}$. This behavior works\n",
    "as a filter to suppress the noise component and yields a physically meaningful solution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. footbibliography::"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3.10.8 ('cherab-phix-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"
  },
  "vscode": {
   "interpreter": {
    "hash": "2725905a4c02db19e04df9b8fdbbe5ec65a73ea52bebaf9474aa1cc98819834c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}