solar-conflux/python/denoising/Main.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9785abf6",
   "metadata": {},
   "source": [
    "We have:\n",
    "$$\n",
    "\\nabla f = (I + \\beta A) \\vec s - \\vec y\n",
    "    = \\vec s - \\vec y + \\beta A \\vec s\n",
    "$$ \n",
    "We can substitute in the partial derivatives:\n",
    "$$ \n",
    "\\left(\\begin{matrix}\n",
    "s_0 - y_0 + \\beta (s_0 - s_1) \\\\\n",
    "s_1 - y_1 + \\beta (2s_1 - s_0 - s_2) \\\\\n",
    "\\vdots \\\\\n",
    "s_N - y_N + \\beta (s_N - s_{N - 1})\n",
    "\\end{matrix}\\right) \n",
    "= \\vec s - \\vec y + \\beta A \\vec s\n",
    "$$\n",
    "\n",
    "Adding $\\vec y - \\vec s$ to both sides and then dividing by $\\beta$ yields:\n",
    "$$ \n",
    "\\left(\\begin{matrix}\n",
    "s_0 - s_1 \\\\\n",
    "2s_1 - s_0 - s_2 \\\\\n",
    "\\vdots \\\\\n",
    "s_N - s_{N - 1}\n",
    "\\end{matrix}\\right) \n",
    "= A \\vec s\n",
    "$$\n",
    "\n",
    "Solving for $A$ is trivial and yields:\n",
    "$$ \n",
    "A = \\left(\\begin{matrix}\n",
    "1  & -1  &   &  &  &  &  \\\\\n",
    "-1 & 2   & -1 &    & &     &  \\\\\n",
    "  & -1  & 2  & -1 & &     &  \\\\\n",
    " & & \\ddots & \\ddots & \\ddots & & \\\\\n",
    " &     & & -1 & 2  & -1 &  \\\\\n",
    " &     & &    & -1 & 2  & -1\\\\\n",
    "  &  &  &  &  & -1 & 1\n",
    "\\end{matrix}\\right) \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e7d69a97",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m=[[1. 2. 0. 0.]\n",
      " [3. 1. 2. 0.]\n",
      " [0. 3. 1. 2.]\n",
      " [0. 0. 3. 1.]]\n",
      "rhs=array([1., 1., 1., 1.])\n",
      "result=array([-0.26315789,  0.63157895,  0.57894737, -0.73684211])\n",
      "[1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import modules.tridiagonal as tridiagonal\n",
    "\n",
    "def sample_with_noise(points, vectorized_g):\n",
    "  \"\"\"\n",
    "  Sample a function at a given set of points,\n",
    "  introducing noise in the process.\n",
    "\n",
    "  Returns a tuple of the real values of the function\n",
    "  and the values with introduced noise.\n",
    "\n",
    "  For performance reasons, this function takes in\n",
    "  a vectorized version of `g`. If you want to pass\n",
    "  a non vectorized function use `np.vectorize`. Eg:\n",
    "  ```py\n",
    "  sample_with_noise(points, np.vectorize(g))\n",
    "  ```\n",
    "  \"\"\"\n",
    "  gs = vectorized_g(points)\n",
    "  ys = gs + 0.05 * np.random.randn(len(gs))\n",
    "\n",
    "  return (gs, ys)\n",
    "\n",
    "def p5(x):\n",
    "  \"\"\"\n",
    "  Vectorized version of the function\n",
    "  described in the pdf.\n",
    "  \"\"\"\n",
    "  n = 5\n",
    "  result = 0\n",
    "\n",
    "  for k in range(n + 1):\n",
    "    result += x**k\n",
    "\n",
    "  return result/(n + 1)\n",
    "\n",
    "def smooth(data, beta):\n",
    "  \"\"\"\n",
    "  Smooth noisy data (y) by means of solving \n",
    "  a minimization problem.\n",
    "\n",
    "  Args:\n",
    "    data (numpy.array): noisy data to be smoothed (y)\n",
    "    beta (float): parameter >= 0 that balances fit and smoothing\n",
    "\n",
    "  Returns:\n",
    "    numpy array of smoothed data (s)\n",
    "  \"\"\"\n",
    "  x = data[0]\n",
    "  y = data[1]\n",
    "\n",
    "  # We need to solve for s in:\n",
    "  #      (I + βA)s - y = 0\n",
    "  #  <=> (I + βA)s = y\n",
    "  #\n",
    "  # Because A and I are tridiagonal, I + βA\n",
    "  # is tridiagonal as well (linear combinations\n",
    "  # of tridiagonal matrices are themselves\n",
    "  # tridiagonal)\n",
    "  #\n",
    "  # This means we can use our existing implementation\n",
    "  # of tridiagonal linear equation system solving.\n",
    "  # \n",
    "  # We start by constructing \n",
    "  #   iba := I + βA\n",
    "  n = len(x) - 1\n",
    "\n",
    "  # The a component of the a-c-e tridiagonal representation of iba\n",
    "  iba_a = np.ones(n + 1) * (2 * beta + 1)\n",
    "  iba_a[0] = 1 + beta\n",
    "  iba_a[-1] = 1 + beta\n",
    "\n",
    "  # I + βA\n",
    "  iba = (\n",
    "    iba_a,\n",
    "    -np.ones(n),\n",
    "    -np.ones(n)\n",
    "  )\n",
    "\n",
    "  data_smoothed = tridiagonal.solve(*iba, y)\n",
    "\n",
    "  return data_smoothed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2b9f2b79",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHgCAYAAABZ+0ykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACPSUlEQVR4nOzdd3xT9frA8U+S7pXS3dJJKRssskGGgIICMhyoXBmi96ooKoqAioJeBVQUN8pV8KcgqCCibEFU9pa9SksZbemAdK/k/P4IDaRN27S0Tcfzfr3yanLyPec8Jymcp9+pUhRFQQghhBCinlDbOgAhhBBCiKokyY0QQggh6hVJboQQQghRr0hyI4QQQoh6RZIbIYQQQtQrktwIIYQQol6R5EYIIYQQ9YokN0IIIYSoVyS5EUIIIUS9IsmNEMJk7NixhIeH2zoMAOLi4lCpVCxatKhC+/Xp04c+ffpUS0x1TXh4OGPHjrV1GELUOEluhLhJixYtQqVSmR5OTk4EBQUxYMAAPvroIzIyMip97O3btzNjxgyuXr1adQGLBuvYsWPMmDGDuLg4W4ciRLWys3UAQtQXb7zxBhERERQUFJCYmMiWLVt47rnneP/991m1ahXt2rWr8DG3b9/OzJkzGTt2LJ6enlUfdC0WFhZGTk4O9vb2Fdpvw4YN1RRR3Xfs2DFmzpxJnz59ak0NnRDVQZIbIarIXXfdRceOHU2vp02bxubNmxk8eDD33HMPx48fx9nZ2YYR1i1FtWAV5eDgUA3RVI7BYCA/P79S1yGEqDxplhKiGvXt25fp06dz7tw5vvvuO9P2Q4cOMXbsWJo0aYKTkxMBAQE8+uijpKammsrMmDGDyZMnAxAREWFq9ipqUli4cCF9+/bFz88PR0dHWrVqxeeff251bCtXrqRNmzY4OTnRpk0bfv75Z4vlDAYD8+bNo3Xr1jg5OeHv789//vMfrly5YlYuPDycwYMHs3XrVjp37oyTkxNNmjTh//7v/0oc8+zZs9x///14eXnh4uJC165dWb16tVkZS31uEhMTGTduHMHBwTg6OhIYGMjQoUPNmlmK97nZsmULKpWKH374gbfeeovg4GCcnJzo168fZ86cKRHbp59+SpMmTXB2dqZz5878/fffVvfjUalUPP300yxevJjWrVvj6OjIunXrALh48SKPPvoo/v7+ODo60rp1a77++usSx/j4449p3bo1Li4uNGrUiI4dO7JkyRLT+6X1i5oxYwYqlarU2BYtWsT9998PwO233276fdqyZQsAe/fuZcCAAfj4+ODs7ExERASPPvpoudcsRG0kNTdCVLNHHnmEl19+mQ0bNvD4448DsHHjRs6ePcu4ceMICAjg6NGjfPnllxw9epSdO3eiUqkYMWIEp06d4vvvv+eDDz7Ax8cHAF9fXwA+//xzWrduzT333IOdnR2//vorTz31FAaDgQkTJpQZ04YNG7j33ntp1aoVs2bNIjU11ZQ0FPef//yHRYsWMW7cOCZOnEhsbCyffPIJBw4cYNu2bWbNRmfOnOG+++5j/PjxjBkzhq+//pqxY8fSoUMHWrduDUBSUhLdu3cnOzubiRMn4u3tzTfffMM999zDTz/9xPDhw0uN+9577+Xo0aM888wzhIeHc/nyZTZu3Eh8fHy5zSyzZ89GrVbz4osvotPpeOeddxg1ahS7du0ylfn88895+umn6dmzJ88//zxxcXEMGzaMRo0aWfxsLNm8eTM//PADTz/9ND4+PoSHh5OUlETXrl1NyY+vry9r165l/PjxpKen89xzzwGwYMECJk6cyH333cezzz5Lbm4uhw4dYteuXTz88MNWnb80vXr1YuLEiXz00Ue8/PLLtGzZEoCWLVty+fJl7rzzTnx9fZk6dSqenp7ExcWxYsWKmzqnEDajCCFuysKFCxVA2bNnT6lltFqt0r59e9Pr7OzsEmW+//57BVD++usv07Z3331XAZTY2NgS5S0dY8CAAUqTJk3KjTk6OloJDAxUrl69atq2YcMGBVDCwsJM2/7++28FUBYvXmy2/7p160psDwsLKxH/5cuXFUdHR+WFF14wbXvuuecUQPn7779N2zIyMpSIiAglPDxc0ev1iqIoSmxsrAIoCxcuVBRFUa5cuaIAyrvvvlvmtfXu3Vvp3bu36fUff/yhAErLli2VvLw80/YPP/xQAZTDhw8riqIoeXl5ire3t9KpUyeloKDAVG7RokUKYHbM0gCKWq1Wjh49arZ9/PjxSmBgoJKSkmK2/cEHH1S0Wq3puxw6dKjSunXrMs8xZswYs++oyOuvv64U/y89LCxMGTNmjOn1jz/+qADKH3/8YVbu559/Lvd3WIi6RJqlhKgBbm5uZqOmbux7k5ubS0pKCl27dgVg//79Vh3zxmPodDpSUlLo3bs3Z8+eRafTlbpfQkICBw8eZMyYMWi1WtP2O+64g1atWpmV/fHHH9Fqtdxxxx2kpKSYHh06dMDNzY0//vjDrHyrVq3o2bOn6bWvry/Nmzfn7Nmzpm1r1qyhc+fO3HbbbaZtbm5u/Pvf/yYuLo5jx46Ver0ODg5s2bKlRJOYNcaNG2fWH6cozqLY9u7dS2pqKo8//jh2dtcrtUeNGkWjRo2sPk/v3r3NPkdFUVi+fDlDhgxBURSzz3HAgAHodDrTd+7p6cmFCxfYs2dPha/vZhR1Vv/tt98oKCio0XMLUR0kuRGiBmRmZuLu7m56nZaWxrPPPou/vz/Ozs74+voSEREBUGZicqNt27bRv39/XF1d8fT0xNfXl5dffrncY5w7dw6AqKioEu81b97c7PXp06fR6XT4+fnh6+tr9sjMzOTy5ctm5UNDQ0scs1GjRmbJyLlz50qcBzA1kxTFV5yjoyNz5sxh7dq1+Pv706tXL9555x0SExNLvdayYitKWIpiKzpv06ZNzcrZ2dlVaGRR0fdYJDk5matXr/Lll1+W+AzHjRsHYPocp0yZgpubG507dyYqKooJEyawbds2q89dWb179+bee+9l5syZ+Pj4MHToUBYuXEheXl61n1uI6iB9boSoZhcuXECn05ndNB944AG2b9/O5MmTiY6Oxs3NDYPBwMCBAzEYDOUeMyYmhn79+tGiRQvef/99QkJCcHBwYM2aNXzwwQdWHcMaBoMBPz8/Fi9ebPH9ov4/RTQajcVyiqJUSTzPPfccQ4YMYeXKlaxfv57p06cza9YsNm/eTPv27cvct7pjK1J8RFzRd/Gvf/2LMWPGWNynaJqAli1bcvLkSX777TfWrVvH8uXL+eyzz3jttdeYOXMmQKmdhvV6faVjVqlU/PTTT+zcuZNff/2V9evX8+ijjzJ37lx27tyJm5tbpY8thC1IciNENfv2228BGDBgAGCsKdi0aRMzZ87ktddeM5U7ffp0iX1Lu5H9+uuv5OXlsWrVKrMaieLNRJaEhYWVer6TJ0+avY6MjOT333+nR48eVTaMPSwsrMR5AE6cOGEWX2kiIyN54YUXeOGFFzh9+jTR0dHMnTvXbDRaZeMCY6fo22+/3bS9sLCQuLi4Ss1TBMYE0N3dHb1eT//+/cst7+rqysiRIxk5ciT5+fmMGDGCt956i2nTpuHk5ESjRo0sTupYWo3XjcoaTQXQtWtXunbtyltvvcWSJUsYNWoUS5cu5bHHHiv32ELUJtIsJUQ12rx5M2+++SYRERGMGjUKuF6DULzGYN68eSX2d3V1BShxM7N0DJ1Ox8KFC8uNKTAwkOjoaL755huz5quNGzeW6O/ywAMPoNfrefPNN0scp7CwsFIzJ999993s3r2bHTt2mLZlZWXx5ZdfEh4eXqLfT5Hs7Gxyc3PNtkVGRuLu7l4lzScdO3bE29ubBQsWUFhYaNq+ePHiSvXxKaLRaLj33ntZvnw5R44cKfF+cnKy6fmNUwGAcc6eVq1aoSiKqS9MZGQkOp2OQ4cOmcolJCSUOpT/RqX9Pl25cqXE72N0dDSANE2JOklqboSoImvXruXEiRMUFhaSlJTE5s2b2bhxI2FhYaxatco0kZuHh4epv0hBQQGNGzdmw4YNxMbGljhmhw4dAHjllVd
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plot\n",
    "\n",
    "# Reusable function to plot the data\n",
    "def plot_denoise_results(β, N, ax):\n",
    "  # Create the data\n",
    "  x = np.linspace(-1, 1, num=N+1)\n",
    "  g, y = sample_with_noise(x, p5)\n",
    "  s = smooth((x, y), β)\n",
    "\n",
    "  ax.set_title(f\"{β=} and {N=}\")\n",
    "  ax.set_xlabel(\"x\")\n",
    "  ax.set_ylabel(\"y\")\n",
    "  ax.plot(x, y, \".\", label=\"raw data points\")\n",
    "  ax.plot(x, g, \"-\", label=\"original function\")\n",
    "  ax.plot(x, s, \"-\", label=\"denoised values\")\n",
    "  ax.legend()\n",
    "\n",
    "\n",
    "N = 100\n",
    "β = 10\n",
    "\n",
    "fig, ax = plot.subplots()\n",
    "fig.suptitle(f\"Data denoising results\")\n",
    "\n",
    "# Plot the data\n",
    "plot_denoise_results(β, N, ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2e5c860",
   "metadata": {},
   "source": [
    "We notice that finding a point where $\\nabla f(s) = 0$ is the same as finding a point where $\\frac 1 \\beta \\nabla f(s) = \\nabla \\frac 1 \\beta f(s) = 0$. We can use this to our advantage, by defining $h(s) = \\frac 1 \\beta f(s)$ and finding the points where $\\nabla h(s) = 0$ instead. When $\\beta \\to \\infty$, we have:\n",
    "$$\n",
    "\\lim_{\\beta \\to \\inf} h _\\beta (s) =\n",
    "\\lim_{\\beta \\to \\inf} \n",
    "  \\frac 1 {2\\beta} \\sum_{k=0}^n (y_k - s_k)^2\n",
    " +\\frac 1 2 \\sum_{k=1}^n (s_k - s_{k - 1})^2\n",
    " = \\frac 1 2 \\sum_{k=1}^n (s_k - s_{k - 1})^2\n",
    "$$\n",
    "\n",
    "It now becomes trivial to verify that $\\nabla f(s) = As$, which means the solution to this will simply be the zero vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d7f7dbdf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gU1dvG8e9mUwkpQBIINaFK7yDSFY0KKGBBQQRR3580FQQFlWYhKKIooihKUUFQQSwgoCAoRUCKgBQBE4IQSigJKaTszvvHmpUlPSTZkNyf68rF7syZmWfHNU/OmVNMhmEYiIiIiIiIiIiISKnn4uwAREREREREREREpHhQY6GIiIiIiIiIiIgAaiwUERERERERERGRf6mxUERERERERERERAA1FoqIiIiIiIiIiMi/1FgoIiIiIiIiIiIigBoLRURERERERERE5F9qLBQRERERERERERFAjYUiIiIiIiIiIiLyLzUWioiIiIiIiIiICKDGQhEREREREREREfmXGgtFRESkxJo0aRImk8nZYdjNnz8fk8lEZGRkoR5TkpWk+2EYBtOnT6dOnTq4u7tTpUoVJkyYgGEYzg5NRERESjE1FoqIiEgG6Q0y6T+enp5UrlyZsLAw3nnnHS5dupSv827evJlJkyZx8eLFgg1YSq3r+Ts1ZswYxowZQ8eOHZk5cybt2rXj5Zdf5ssvv3R2aCIiIlKKqbFQREREsvTSSy/x6aef8v777zNixAgAnn76aRo3bsyePXvyfL7NmzczefLk67JhpyAMGDCApKQkatSoUajHlCbX63cqIiKCGTNm8OKLLzJ37lz+97//8fnnn1O2bFl+++03Z4cnIiIipZirswMQERGR4uuOO+6gVatW9vfjxo1j3bp19OjRg7vuuosDBw7g5eXlxAivL2azGbPZXOjHFKaEhAS8vb2dHcZ17+uvv8YwDHsjPICrq+1Pc/0/JSIiIs6knoUiIiKSJzfffDPjx4/n2LFjfPbZZwAcO3aMoUOHUq9ePby8vKhQoQL33Xefw7xykyZNYsyYMQCEhobahzinl8nNObKzceNGWrdujaenJ7Vq1eKDDz7IsuyJEycYPHgwFStWxMPDg4YNGzJ37lyHMunzHR45coRBgwbh7++Pn58fjzzyCImJiRnOuWvXLu644w58fX0pW7Yst9xyS4YeYpnNt3fp0iWefvppQkJC8PDwICgoiFtvvZWdO3dmeUxeYlu/fj2tWrVyuC+5ncsxvdz+/fvp168f5cqVo0OHDnm6jzl9PoBBgwYREhKS5fWziy+771Rurp2Z3H4X83t+gK1bt3LDDTcQGBho3/bHH38QHx9PkyZNcjxeREREpLCoZ6GIiIjk2YABA3j++edZs2YNjz/+ONu3b2fz5s088MADVK1alcjISN5//326dOnC/v37KVOmDH369OGvv/7i888/56233iIgIADA3liSm3NkZe/evdx2220EBgYyadIk0tLSmDhxIhUrVsxQ9vTp09x4442YTCaGDx9OYGAgP/zwA48++ihxcXE8/fTTDuXvv/9+QkNDCQ8PZ+fOnXz00UcEBQXx2muv2cv8+eefdOzYEV9fX5599lnc3Nz44IMP6NKlCxs2bKBt27ZZxv7EE0/w1VdfMXz4cBo0aMC5c+fYuHEjBw4coEWLFtn+d8gptl27dnH77bcTHBzM5MmTsVgsvPTSSw4NVLlx3333UadOHaZMmWJffCO39/FaPl9OcvpO5ffauf0uXstn27t3L82aNXPYNm3aNDw9Pbn11ltz9fnHjx9P9erVefzxx7MsExcXh6+vr8O21NRU3NzccnUNERERKYUMERERkavMmzfPAIzt27dnWcbPz89o3ry5YRiGkZiYmGH/li1bDMD45JNP7NumTZtmAEZERESG8rk9R2Z69epleHp6GseOHbNv279/v2E2m42r/9x59NFHjeDgYCMmJsZh+wMPPGD4+fnZ45g4caIBGIMHD3Yo17t3b6NChQoZru/u7m4cPXrUvu3kyZOGj4+P0alTJ/u29Pt65ef38/Mzhg0bluVny+yY3MbWs2dPo0yZMsaJEyfs2w4fPmy4urpmuC+ZSb/Ogw8+mGFfbu9jTp/PMAxj4MCBRo0aNbK8/pWuvh/Zfadyc+3M5Pa7mN/zJycnG66urkZ4eLiRkJBg/Pbbb8bDDz9sAMaECRNyfZ4RI0YYJpPJmDdvXoZ9v/32m1GjRg0DMLp06WKcPXvWOHjwoNGoUSPDZDIZrVu3Nv7+++88xy4iIiIln4Yhi4iISL6ULVvWvirylXOspaamcu7cOWrXro2/v3+uhmReyzksFgurV6+mV69eVK9e3b69fv36hIWFOZQ1DIOlS5fSs2dPDMMgJibG/hMWFkZsbGyGaz3xxBMO7zt27Mi5c+eIi4uzX3/NmjX06tWLmjVr2ssFBwfTr18/Nm7caC+bGX9/f7Zu3crJkyezuTuZyy42i8XCTz/9RK9evahcubK9TO3atbnjjjuu6Tp5uY/X8vmuVX6vndvvYn7Pf+DAAdLS0mjSpAnTp0/nxhtv5JNPPqFevXo89dRT9nKXL1/O9uf1119n4MCBPProoyxatMh+nGEYPPbYY3zwwQfExMQQEBDAwIEDuf/+++nVqxenT5/miSee4H//+1+e4hYREZHSQY2FIiIiki/x8fH4+PgAkJSUxIQJE6hWrRoeHh4EBAQQGBjIxYsXiY2NzdX58nuOs2fPkpSURJ06dTLsq1evXoayFy9e5MMPPyQwMNDh55FHHgHgzJkzDsdc2QAJUK5cOQAuXLhgP2diYmKGa4GtwdJqtXL8+PEs43/99dfZt28f1apVo02bNkyaNIm///47y/K5je3MmTMkJSVRu3btDMdlti07oaGhDu/zch+v5fNdq/xeO7ffxfyef+/evQA0adKE3r17s2TJEkaPHs2JEydo27YtSUlJxMfH4+XllePP/PnzsVqtPPzww5w6dQqAf/75h+rVqxMWFkaFChX49NNPOXDgAKGhobz88ssEBgYyePBgzp49S1paWj7vroiIiJRUmrNQRERE8uyff/4hNjbW3ug0YsQI5s2bx9NPP027du3w8/PDZDLxwAMPYLVac3XOgjhHTtLP89BDDzFw4MBMy1y9uERWKxEb/87dd63uv/9+OnbsyNdff82aNWuYNm0ar732GsuWLcuxB2Bhx5bu6tV583Ifc/P5slrExGKxXFPc+b23uf0u5vf8+/bto3z58lStWpWqVavSqFEj7r//fpo0acLDDz/Mtm3baN++PfPmzcvxM65evZrFixfTp08fh7kor7ynnp6eBAUFce7cOQzDcNhX0N8VERERuf6psVBERETy7NNPPwWwD/P96quvGDhwINOnT7eXuXz5MhcvXnQ4LruVbXN7jqsFBgbi5eXF4cOHM+w7dOhQhrI+Pj5YLBa6deuW7XlzKzAwkDJlymS4FsDBgwdxcXGhWrVq2Z4jODiYoUOHMnToUM6cOUOLFi149dVX8zxc+EpBQUF4enpy5MiRDPsy25YXeb2POX2+cuXKZfrf+dixYzmeO6dVnfNzb/PyXczP+ffu3Uvjxo0zbE9vHA0ODsbV1ZVBgwZl+9l+/PFHvv76a3r16sWiRYvsjcdVq1bl2LFjrFq1ijZt2vDKK6/g7+/PqVOnGDlyJC+99BLfffcd5cuX10InIiIikoGGIYuIiEierFu3jpdffpnQ0FD69+8P2Hq4Xd1DaebMmRl6hnl7ewNk2uiS23NkdlxYWBjLly8nKirKvv3AgQOsXr06Q9l77rmHpUuXsm/fvgznOnv2bLbXyur6t912G9988w2RkZH27adPn2bRokV06NAhw2q06SwWS4Yh1kFBQVSuXJnk5OQ8x3J1XN26dWP58uUOc+odOXKEH3744ZrPnZv7mNvPV6tWLWJjY9mzZ499W3R0NF9//XWOsWT1nbqWe5ub7+K1nH/v3r2cOnXKoVxKSgqzZs2iWbNm1K1bN9vj073yyit069aNJUuW4Or6Xx8Ak8nE3LlzGTFiBBUqVGDTpk3Mnz+fpUuX8tNPP+H
      "text/plain": [
       "<Figure size 1280x960 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plot.subplots(\n",
    "    3,\n",
    "    3,\n",
    "    layout=\"constrained\",\n",
    "    figsize=(12.8, 9.6)\n",
    ")\n",
    "\n",
    "fig.suptitle(r\"Data denoising results as $\\beta \\to \\infty$\")\n",
    "\n",
    "N = 100\n",
    "\n",
    "for i in range(len(axs.flat)):\n",
    "    ax = axs.flat[i]\n",
    "    plot_denoise_results(2**i, N, ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfc485d5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}