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": 14,
   "id": "e7d69a97",
   "metadata": {},
   "outputs": [],
   "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,y = data\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",
    "    -beta*np.ones(n),\n",
    "    -beta*np.ones(n)\n",
    "  )\n",
    "\n",
    "  data_smoothed = tridiagonal.solve(*iba, y)\n",
    "\n",
    "  return data_smoothed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "2b9f2b79",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHgCAYAAABZ+0ykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACSs0lEQVR4nOzdd3hT5dvA8W+S7pXuAZS27E2RDbIERBkKyBBRpr6KICKCjJ8oOFiKIooiqOBCEBkiIEOGypCNIHu0tEALlJZ0r+S8f4QGQndpm6a9P165IOc855z7nFZy55kqRVEUhBBCCCHKCbWlAxBCCCGEKE6S3AghhBCiXJHkRgghhBDliiQ3QgghhChXJLkRQgghRLkiyY0QQgghyhVJboQQQghRrkhyI4QQQohyRZIbIYQQQpQrktwIIUyGDRtGcHCwpcMAIDw8HJVKxbJlywp1XMeOHenYsWOJxGRtgoODGTZsmKXDEKLUSXIjxANatmwZKpXK9HJwcKBSpUp069aNBQsWkJCQUORz7927l+nTp3P79u3iC1hUWKdOnWL69OmEh4dbOhQhSpSNpQMQorx45513CAkJISMjg+joaHbt2sW4ceP46KOPWL9+PY0aNSr0Offu3cuMGTMYNmwY7u7uxR90GRYUFERKSgq2traFOm7r1q0lFJH1O3XqFDNmzKBjx45lpoZOiJIgyY0QxeTxxx+nWbNmpvdTpkxhx44d9OzZkyeeeILTp0/j6OhowQitS1YtWGHZ2dmVQDRFYzAYSE9PL9J9CCGKTpqlhChBjzzyCNOmTePy5cv88MMPpu3Hjx9n2LBhVKtWDQcHB/z9/RkxYgS3bt0ylZk+fToTJ04EICQkxNTsldWksHTpUh555BF8fX2xt7enXr16fPHFFwWObd26dTRo0AAHBwcaNGjA2rVrcyxnMBiYP38+9evXx8HBAT8/P1588UXi4uLMygUHB9OzZ092795NixYtcHBwoFq1anz33XfZznnp0iX69++Pp6cnTk5OtGrVio0bN5qVyanPTXR0NMOHD6dKlSrY29sTEBDAk08+adbMcn+fm127dqFSqfj55595//33qVKlCg4ODnTu3JkLFy5ki23hwoVUq1YNR0dHWrRowd9//13gfjwqlYoxY8bw448/Ur9+fezt7dm8eTMAV69eZcSIEfj5+WFvb0/9+vX55ptvsp3j008/pX79+jg5OeHh4UGzZs1Yvny5aX9u/aKmT5+OSqXKNbZly5bRv39/ADp16mT6fdq1axcAhw4dolu3bnh7e+Po6EhISAgjRozI956FKIuk5kaIEvbcc88xdepUtm7dygsvvADAtm3buHTpEsOHD8ff35+TJ0+yePFiTp48yT///INKpaJv376cO3eOn376iY8//hhvb28AfHx8APjiiy+oX78+TzzxBDY2Nvz222+8/PLLGAwGRo8enWdMW7du5amnnqJevXrMmjWLW7dumZKG+7344ossW7aM4cOHM3bsWMLCwvjss884evQoe/bsMWs2unDhAv369WPkyJEMHTqUb775hmHDhtG0aVPq168PwPXr12nTpg3JycmMHTsWLy8vvv32W5544gl++eUX+vTpk2vcTz31FCdPnuSVV14hODiYGzdusG3bNiIiIvJtZpk9ezZqtZoJEyag0+mYO3cugwcPZv/+/aYyX3zxBWPGjKFdu3a89tprhIeH07t3bzw8PHJ8NjnZsWMHP//8M2PGjMHb25vg4GCuX79Oq1atTMmPj48Pv//+OyNHjiQ+Pp5x48YBsGTJEsaOHUu/fv149dVXSU1N5fjx4+zfv59nnnmmQNfPTfv27Rk7diwLFixg6tSp1K1bF4C6dety48YNHn30UXx8fJg8eTLu7u6Eh4ezZs2aB7qmEBajCCEeyNKlSxVAOXjwYK5ltFqt0qRJE9P75OTkbGV++uknBVD++usv07YPPvhAAZSwsLBs5XM6R7du3ZRq1arlG3NoaKgSEBCg3L5927Rt69atCqAEBQWZtv39998KoPz4449mx2/evDnb9qCgoGzx37hxQ7G3t1def/1107Zx48YpgPL333+btiUkJCghISFKcHCwotfrFUVRlLCwMAVQli5dqiiKosTFxSmA8sEHH+R5bx06dFA6dOhger9z504FUOrWraukpaWZtn/yyScKoJw4cUJRFEVJS0tTvLy8lObNmysZGRmmcsuWLVMAs3PmBlDUarVy8uRJs+0jR45UAgIClJiYGLPtTz/9tKLVak0/yyeffFKpX79+ntcYOnSo2c8oy9tvv63c/096UFCQMnToUNP7VatWKYCyc+dOs3Jr167N93dYCGsizVJClAIXFxezUVP39r1JTU0lJiaGVq1aAXDkyJECnfPec+h0OmJiYujQoQOXLl1Cp9PlelxUVBTHjh1j6NChaLVa0/auXbtSr149s7KrVq1Cq9XStWtXYmJiTK+mTZvi4uLCzp07zcrXq1ePdu3amd77+PhQu3ZtLl26ZNq2adMmWrRowcMPP2za5uLiwv/93/8RHh7OqVOncr1fOzs7du3ala1JrCCGDx9u1h8nK86s2A4dOsStW7d44YUXsLG5W6k9ePBgPDw8CnydDh06mD1HRVFYvXo1vXr1QlEUs+fYrVs3dDqd6Wfu7u7OlStXOHjwYKHv70FkdVbfsGEDGRkZpXptIUqCJDdClILExERcXV1N72NjY3n11Vfx8/PD0dERHx8fQkJCAPJMTO61Z88eunTpgrOzM+7u7vj4+DB16tR8z3H58mUAatasmW1f7dq1zd6fP38enU6Hr68vPj4+Zq/ExERu3LhhVr5q1arZzunh4WGWjFy+fDnbdQBTM0lWfPezt7dnzpw5/P777/j5+dG+fXvmzp1LdHR0rveaV2xZCUtWbFnXrVGjhlk5GxubQo0syvo5Zrl58ya3b99m8eLF2Z7h8OHDAUzPcdKkSbi4uNCiRQtq1qzJ6NGj2bNnT4GvXVQdOnTgqaeeYsaMGXh7e/Pkk0+ydOlS0tLSSvzaQpQE6XMjRAm7cuUKOp3O7ENzwIAB7N27l4kTJxIaGoqLiwsGg4HHHnsMg8GQ7zkvXrxI586dqVOnDh999BGBgYHY2dmxadMmPv744wKdoyAMBgO+vr78+OOPOe7P6v+TRaPR5FhOUZRiiWfcuHH06tWLdevWsWXLFqZNm8asWbPYsWMHTZo0yfPYko4ty/0j4rJ+Fs8++yxDhw7N8ZisaQLq1q3L2bNn2bBhA5s3b2b16tV8/vnnvPXWW8yYMQMg107Der2+yDGrVCp++eUX/vnnH3777Te2bNnCiBEjmDdvHv/88w8uLi5FPrcQliDJjRAl7PvvvwegW7dugLGmYPv27cyYMYO33nrLVO78+fPZjs3tg+y3334jLS2N9evXm9VI3N9MlJOgoKBcr3f27Fmz99WrV+ePP/6gbdu2xTaMPSgoKNt1AM6cOWMWX26qV6/O66+/zuuvv8758+cJDQ1l3rx5ZqPRihoXGDtFd+rUybQ9MzOT8PDwIs1TBMYE0NXVFb1eT5cuXfIt7+zszMCBAxk4cCDp6en07duX999/nylTpuDg4ICHh0eOkzrmVuN1r7xGUwG0atWKVq1a8f7777N8+XIGDx7MihUreP755/M9txBliTRLCVGCduzYwbvvvktISAiDBw8G7tYg3F9jMH/+/GzHOzs7A2T7MMvpHDqdjqVLl+YbU0BAAKGhoXz77bdmzVfbtm3L1t9lwIAB6PV63n333WznyczMLNLMyd27d+fAgQPs27fPtC0pKYnFixcTHBycrd9PluTkZFJTU822Va9eHVdX12JpPmnWrBleXl4sWbKEzMxM0/Yff/yxSH18smg0Gp566ilWr17Nf//9l23/zZs3TX+/dyoAMM7ZU69ePRRFMfWFqV69OjqdjuPHj5vKRUVF5TqU/165/T7FxcVl+30MDQ0FkKYpYZWk5kaIYvL7779z5swZMjMzuX79Ojt27GDbtm0EBQWxfv1600Rubm5upv4iGRkZVK5cma1btxIWFpbtnE2bNgXgf//7H08//TS2trb
      "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": 20,
   "id": "d7f7dbdf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gU1dfA8e/upnfSAyQkVOk1Ih0EDCggYEFAurw/KdKkSxcJgiCKKIo0QSQqRZTeEUF6ldATQgkJIZBed+f9I2ZlSYckm3I+z7OP7MydO2fHZQ/3zp17VYqiKAghhBBCCCGEEEIIIUo9tbEDEEIIIYQQQgghhBBCFA3SWSiEEEIIIYQQQgghhACks1AIIYQQQgghhBBCCPEv6SwUQgghhBBCCCGEEEIA0lkohBBCCCGEEEIIIYT4l3QWCiGEEEIIIYQQQgghAOksFEIIIYQQQgghhBBC/Es6C4UQQgghhBBCCCGEEIB0FgohhBBCCCGEEEIIIf4lnYVCCCGEEEIIIYQQQghAOguFEEIIIYQQQgghhBD/ks5CIYQQQpRYM2bMQKVSGTsMvVWrVqFSqQgODi7QY0qyknQ9FEVhwYIFVKlSBTMzM8qVK8e0adNQFMXYoQkhhBCiFJPOQiGEEEJkkN4hk/6ysLCgbNmy+Pn58eWXXxITE/NM9R45coQZM2bw+PHj/A1YlFrF+Ts1btw4xo0bR4sWLVi8eDFNmjTh448/5pdffjF2aEIIIYQoxaSzUAghhBBZmjVrFmvWrOGbb77hgw8+AGDUqFHUrl2b8+fP57m+I0eOMHPmzGLZsZMf+vTpQ0JCAhUqVCjQY0qT4vqdCgoKYtGiRUyZMoUVK1bwv//9j59++gkbGxv+/vtvY4cnhBBCiFLMxNgBCCGEEKLo6tixI40aNdK/nzRpEvv27aNTp0506dKFwMBALC0tjRhh8aLRaNBoNAV+TEGKi4vD2tra2GEUe5s2bUJRFH0nPICJSdo/zeXvlBBCCCGMSUYWCiGEECJPXn75ZaZOncqtW7dYu3YtALdu3WLo0KFUq1YNS0tLnJyceOuttwzmlZsxYwbjxo0DwMfHR/+Ic3qZ3NSRncOHD+Pr64uFhQWVKlXi22+/zbLs3bt3GThwIG5ubpibm1OzZk1WrFhhUCZ9vsPr16/Tv39/HBwcsLe3Z8CAAcTHx2eo88yZM3Ts2BE7OztsbGxo27ZthhFimc23FxMTw6hRo/D29sbc3BxXV1fat2/P6dOnszwmL7EdOHCARo0aGVyX3M7lmF7u0qVL9OrVizJlytC8efM8XcecPh9A//798fb2zvL82cWX3XcqN+fOTG6/i89aP8CxY8d44YUXcHFx0W87d+4csbGx1KlTJ8fjhRBCCCEKiowsFEIIIUSe9enTh8mTJ7Nr1y4GDx7MiRMnOHLkCO+88w7ly5cnODiYb775htatW3Pp0iWsrKzo3r07V69e5aeffuLzzz/H2dkZQN9Zkps6snLhwgVeeeUVXFxcmDFjBqmpqUyfPh03N7cMZcPCwnjppZdQqVQMHz4cFxcXtm/fzqBBg4iOjmbUqFEG5d9++218fHzw9/fn9OnTfP/997i6uvLpp5/qy/zzzz+0aNECOzs7xo8fj6mpKd9++y2tW7fm4MGDNG7cOMvY33//fX799VeGDx9OjRo1ePjwIYcPHyYwMJAGDRpk+/8hp9jOnDlDhw4d8PDwYObMmWi1WmbNmmXQQZUbb731FlWqVGHOnDn6xTdyex2f5/PlJKfv1LOeO7ffxef5bBcuXKBevXoG2+bPn4+FhQXt27fP1eefOnUqXl5eDB48OMsy0dHR2NnZGWxLSUnB1NQ0V+cQQgghRCmkCCGEEEI8ZeXKlQqgnDhxIssy9vb2Sv369RVFUZT4+PgM+48ePaoAyg8//KDfNn/+fAVQgoKCMpTPbR2Z6dq1q2JhYaHcunVLv+3SpUuKRqNRnv7nzqBBgxQPDw8lIiLCYPs777yj2Nvb6+OYPn26AigDBw40KNetWzfFyckpw/nNzMyUGzdu6Lfdu3dPsbW1VVq2bKnfln5dn/z89vb2yrBhw7L8bJkdk9vYOnfurFhZWSl3797Vb7t27ZpiYmKS4bpkJv08PXv2zLAvt9cxp8+nKIrSr18/pUKFClme/0lPX4/svlO5OXdmcvtdfNb6k5KSFBMTE8Xf31+Ji4tT/v77b6Vv374KoEybNi3X9XzwwQeKSqVSVq5cmWHf33//rVSoUEEBlNatWysPHjxQLl++rNSqVUtRqVSKr6+vcvPmzTzHLoQQQoiSTx5DFkIIIcQzsbGx0a+K/OQcaykpKTx8+JDKlSvj4OCQq0cyn6cOrVbLzp076dq1K15eXvrt1atXx8/Pz6Csoihs2LCBzp07oygKERER+pefnx9RUVEZzvX+++8bvG/RogUPHz4kOjpaf/5du3bRtWtXKlasqC/n4eFBr169OHz4sL5sZhwcHDh27Bj37t3L5upkLrvYtFote/bsoWvXrpQtW1ZfpnLlynTs2PG5zpOX6/g8n+95Peu5c/tdfNb6AwMDSU1NpU6dOixYsICXXnqJH374gWrVqjFy5Eh9ucTExGxf8+bNo1+/fgwaNIh169bpj1MUhffee49vv/2WiIgInJ2d6devH2+//TZdu3YlLCyM999/n//97395ilsIIYQQpYN0FgohhBDimcTGxmJrawtAQkIC06ZNw9PTE3Nzc5ydnXFxceHx48dERUXlqr5nrePBgwckJCRQpUqVDPuqVauWoezjx4/57rvvcHFxMXgNGDAAgPDwcINjnuyABChTpgwAjx490tcZHx+f4VyQ1mGp0+m4fft2lvHPmzePixcv4unpyYsvvsiMGTO4efNmluVzG1t4eDgJCQlUrlw5w3GZbcuOj4+Pwfu8XMfn+XzP61nPndvv4rPWf+HCBQDq1KlDt27dCAgIYOzYsdy9e5fGjRuTkJBAbGwslpaWOb5WrVqFTqejb9++3L9/H4A7d+7g5eWFn58fTk5OrFmzhsDAQHx8fPj4449xcXFh4MCBPHjwgNTU1Ge8ukIIIYQoqWTOQiGEEELk2Z07d4iKitJ3On3wwQesXLmSUaNG0aRJE+zt7VGpVLzzzjvodLpc1ZkfdeQkvZ53332Xfv36ZVrm6cUlslqJWPl37r7n9fbbb9OiRQs2bdrErl27mD9/Pp9++ikbN27McQRgQceW7unVefNyHXPz+bJaxESr1T5X3M96bXP7XXzW+i9evIijoyPly5enfPny1KpVi7fffps6derQt29fjh8/TrNmzVi5cmWOn3Hnzp2sX7+e7t27G8xF+eQ1tbCwwNXVlYcPH6IoisG+/P6uCCGEEKL4k85CIYQQQuTZmjVrAPSP+f7666/069ePBQsW6MskJiby+PFjg+OyW9k2t3U8zcXFBUtLS65du5Zh35UrVzKUtbW1RavV0q5du2zrzS0XFxesrKwynAvg8uXLqNVqPD09s63Dw8ODoUOHMnToUMLDw2nQoAGffPJJnh8XfpKrqysWFhZcv349w77MtuVFXq9jTp+vTJkymf5/vnXrVo5157Sq87Nc27x8F5+l/gsXLlC7du0M29M7Rz08PDAxMaF///7Zfrbdu3ezadMmunbtyrp16/Sdx+XLl+fWrVvs2LGDF198kdmzZ+Pg4MD9+/cZPXo0s2bN4vfff8fR0VEWOhFCCCFEBvIYshBCCCHyZN++fXz88cf4+PjQu3dvIG2E29MjlBYvXpxhZJi1tTVApp0uua0js+P8/PzYvHkzISEh+u2BgYHs3LkzQ9k33niDDRs2cPHixQx1PXjwINtzZXX+V155hd9++43g4GD99rCwMNatW0fz5s0zrEabTqvVZnjE2tXVlbJly5KUlJTnWJ6Oq127dmzevNlgTr3r16+zffv25647N9cxt5+vUqVKREVFcf78ef220NBQNm3alGMsWX2nnufa5ua7+Dz1X7hwgfv37xuUS05OZsmSJdSrV4+qVatme3y62bNn065dOwICAjAx+W8MgEqlYsWKFXzwwQc4OTnx119/sWr
      "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(4**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
}