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": 17,
   "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": 18,
   "id": "2b9f2b79",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHgCAYAAABZ+0ykAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACM70lEQVR4nOzdd3QU5dfA8e/spmzqhiSkQUhC71JCly5dpCkg/iiKvqIoSpGiIkUFUVCwYQUsCCjNQu8d6QiETiCUQEgghfTszvvHmpUlPSTZlPs5Zw/szDMzd3YDc/NURVVVFSGEEEKIUkJj7QCEEEIIIQqSJDdCCCGEKFUkuRFCCCFEqSLJjRBCCCFKFUluhBBCCFGqSHIjhBBCiFJFkhshhBBClCqS3AghhBCiVJHkRgghhBCliiQ3QjykRYsWoSiK+aXT6fDx8aF9+/bMnDmTiIiIfJ87JCSEqVOncvny5YILOBvbt29HURS2b99eJNfLSbt27WjXrl2ejilu92BNmX0Wa9euZerUqVaLSYiiIMmNEAVk4cKF7Nu3j02bNvHFF1/QoEEDZs2aRa1atdi8eXO+zhkSEsK0adOKLLkpbr788ku+/PLLPB3TqFEj9u3bR6NGjQopqpJt7dq1TJs2zdphCFGobKwdgBClRd26dQkODja/79evH6NHj+bRRx+lb9++nD9/Hm9vbytGWPLUrl07z8e4urrSvHnzQogmfxITE9HpdCiKYu1QhCgzpOZGiEJUqVIl5syZQ1xcHF9//bV5+6FDhxg4cCCBgYE4ODgQGBjI008/zZUrV8xlFi1axFNPPQVA+/btzc1eixYtAmDTpk306tWLihUrotPpqFq1Ki+++CKRkZG5iu3MmTN07doVR0dHPD09GTFiBHFxcZmW3bx5Mx07dsTV1RVHR0datWrFli1bLMpMnToVRVE4deoUTz/9NHq9Hm9vb5577jliYmIsyiYlJTFp0iSCgoKws7OjQoUKjBw5kujoaItymTVLzZ8/n0ceeQRnZ2dcXFyoWbMmb775pnl/Zk0xw4YNw9nZmQsXLtC9e3ecnZ3x9/dn7NixJCcnW5z/2rVrPPnkk7i4uODm5sYzzzzDwYMHLT77rKQ3UW7cuJHnnnuO8uXL4+joaL7GsmXLaNGiBU5OTjg7O9OlSxeOHj1qcY5Lly4xcOBA/Pz8sLe3x9vbm44dO3Ls2DFzGUVRMm1aCgwMZNiwYVnGN2zYML744gvzOdJf6TWDv/32G82aNUOv1+Po6EjlypV57rnnsr1nIYojqbkRopB1794drVbLzp07zdsuX75MjRo1GDhwIO7u7oSHhzN//nyaNGlCSEgInp6e9OjRgxkzZvDmm2/yxRdfmJtZqlSpAsDFixdp0aIFzz//PHq9nsuXL/Pxxx/z6KOPcuLECWxtbbOM6datW7Rt2xZbW1u+/PJLvL29Wbx4Ma+88kqGsj///DNDhgyhV69e/PDDD9ja2vL111/TpUsXNmzYQMeOHS3K9+vXjwEDBjB8+HBOnDjBpEmTAFiwYAEAqqrSu3dvtmzZwqRJk2jdujX//PMPU6ZMYd++fezbtw97e/tM4166dCkvv/wyr776KrNnz0aj0XDhwgVCQkJy/B5SU1N54oknGD58OGPHjmXnzp28++676PV63nnnHQDi4+Np3749d+7cYdasWVStWpX169czYMCAHM9/v+eee44ePXrw008/ER8fj62tLTNmzODtt9/m2Wef5e233yYlJYWPPvqI1q1bc+DAAXMtVffu3TEYDHz44YdUqlSJyMhI9u7dmyHxy4/JkycTHx/P8uXL2bdvn3m7r68v+/btY8CAAQwYMICpU6ei0+m4cuUKW7dufejrClHkVCHEQ1m4cKEKqAcPHsyyjLe3t1qrVq0s96elpan37t1TnZyc1Hnz5pm3//bbbyqgbtu2LdsYjEajmpqaql65ckUF1N9//z3b8hMmTFAVRVGPHTtmsb1Tp04W14uPj1fd3d3Vnj17WpQzGAzqI488ojZt2tS8bcqUKSqgfvjhhxZlX375ZVWn06lGo1FVVVVdv359puWWLVumAuo333xj3ta2bVu1bdu25vevvPKK6ubmlu29bdu2LcNnNnToUBVQf/31V4uy3bt3V2vUqGF+/8UXX6iAum7dOotyL774ogqoCxcuzPba6T8LQ4YMsdgeFham2tjYqK+++qrF9ri4ONXHx0ft37+/qqqqGhkZqQLq3Llzs70OoE6ZMiXD9oCAAHXo0KHm95l9FiNHjlQz+69/9uzZKqBGR0dne20hSgJplhKiCKiqavH+3r17TJgwgapVq2JjY4ONjQ3Ozs7Ex8dz+vTpXJ0zIiKCESNG4O/vj42NDba2tgQEBADkeI5t27ZRp04dHnnkEYvtgwYNsni/d+9e7ty5w9ChQ0lLSzO/jEYjXbt25eDBg8THx1sc88QTT1i8r1+/PklJSeZRY+k1AQ82nzz11FM4OTllaO66X9OmTYmOjubpp5/m999/z3UTHJiaYXr27JkhtvubAnfs2IGLiwtdu3a1KPf000/n+jpgqr2634YNG0hLS2PIkCEWn6NOp6Nt27bmJjR3d3eqVKnCRx99xMcff8zRo0cxGo15unZ+NWnSBID+/fvz66+/cv369SK5rhCFQZIbIQpZfHw8UVFR+Pn5mbcNGjSIzz//nOeff54NGzZw4MABDh48SPny5UlMTMzxnEajkc6dO7Ny5UrGjx/Pli1bOHDgAPv37wfI8RxRUVH4+Phk2P7gtlu3bgHw5JNPYmtra/GaNWsWqqpy584di2M8PDws3qc3MaXHFBUVhY2NDeXLl7copygKPj4+REVFZRn34MGDWbBgAVeuXKFfv354eXnRrFkzNm3alO39Ajg6OqLT6TLElpSUZH4fFRWVaafvvHYE9/X1tXif/jk2adIkw+e4bNkyc5KmKApbtmyhS5cufPjhhzRq1Ijy5cszatSoLPtDFZQ2bdqwevVqcxJWsWJF6taty5IlSwr1ukIUBulzI0QhW7NmDQaDwdwxNiYmhr/++ospU6YwceJEc7nk5OQMiUJWTp48yfHjx1m0aBFDhw41b79w4UKujvfw8ODmzZsZtj+4zdPTE4DPPvssyxFIeX3we3h4kJaWxu3bty0SHFVVuXnzprkGISvPPvsszz77LPHx8ezcuZMpU6bw+OOPc+7cOXPNVX55eHhw4MCBDNsz+6yy8+DIqPTPcfny5TnGGBAQwPfffw/AuXPn+PXXX5k6dSopKSl89dVXgCkpe7AjNJBtYpgbvXr1olevXiQnJ7N//35mzpzJoEGDCAwMpEWLFg91biGKktTcCFGIwsLCGDduHHq9nhdffBEwPfhUVc3Qafa7777DYDBYbHuw1iNd+sPzwXPcPyIrO+3bt+fUqVMcP37cYvsvv/xi8b5Vq1a4ubkREhJCcHBwpi87O7tcXTNdegfkn3/+2WL7ihUriI+Pz9BBOStOTk5069aNt956i5SUFE6dOpWnODLTtm1b4uLiWLduncX2pUuXPtR5u3Tpgo2NDRcvXszyc8xM9erVefvtt6lXrx5Hjhwxbw8MDOSff/6xKLt161bu3buXYyxZ/Uw9WKZt27bMmjULIMOILiGKO6m5EaKAnDx50tyXIiIigl27drFw4UK0Wi2rVq0y11K4urrSpk0bPvroIzw9PQkMDGTHjh18//33uLm5WZyzbt26AHzzzTe4uLig0+kICgqiZs2aVKlShYkTJ6KqKu7u7vz555+5ap4BeP3111mwYAE9evTgvffeM4+WOnPmjEU5Z2dnPvvsM4YOHcqdO3d48skn8fLy4vbt2xw/fpzbt28zf/78PH1OnTp1okuXLkyYMIHY2FhatWplHi3VsGFDBg8enOWxL7zwAg4ODrRq1QpfX19u3rzJzJkz0ev1Odb45MbQoUP55JNP+N///sd7771H1apVWbduHRs2bABAo8nf74OBgYFMnz6dt956i0uXLtG1a1fKlSvHrVu3OHDgAE5OTkybNo1//vmHV155haeeeopq1aphZ2fH1q1b+eeffyxq+QYPHszkyZN
      "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(beta, 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), beta)\n",
    "\n",
    "  ax.set_title(f\"{beta=} 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",
    "beta = 10\n",
    "\n",
    "fig, ax = plot.subplots()\n",
    "fig.suptitle(f\"Data denoising results\")\n",
    "\n",
    "# Plot the data\n",
    "plot_denoise_results(beta, 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": 19,
   "id": "d7f7dbdf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQsAAAPLCAYAAAD4+99NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gU1dfA8e/spndIAoQSEnow9AAC0kQ6SFNCeaUI/kQUlA4qVaWDoAIqChFFQCEgCtKkg0BoShcwEEoghJLednfeP0JWlvSQZFPO53n2MTtzZ+bMsO7ZuXOLoqqqihBCCCGEEEIIIYQQotjTmDsAIYQQQgghhBBCCCFEwSCVhUIIIYQQQgghhBBCCEAqC4UQQgghhBBCCCGEEI9JZaEQQgghhBBCCCGEEAKQykIhhBBCCCGEEEIIIcRjUlkohBBCCCGEEEIIIYQApLJQCCGEEEIIIYQQQgjxmFQWCiGEEEIIIYQQQgghAKksFEIIIYQQQgghhBBCPCaVhUIIIYQQQgghhBBCCEAqC4UQQgghhBBCCCGEEI9JZaEQQgghUgkICEBRFOPLxsaGMmXK0Lp1a2bNmkVYWFiO93348GGmTZvGo0ePci/gDEybNg1FUfLlWJlJua7Xrl3L022KsrSuR35/pnJLbGws06dPp0aNGtjY2ODq6krfvn15+PChuUMTQgghRDEmlYVCCCGESNfKlSv5888/2blzJ0uWLKFu3brMmTMHHx8fdu3alaN9Hj58mOnTpxe6ip3c0LlzZ/788088PDzydJvipjB+plRVpW/fvixatIg333yT33//nSlTprBu3ToWLFhg7vCEEEIIUYxZmDsAIYQQQhRcvr6++Pn5Gd/36tWLUaNG8cILL9CzZ08uX75M6dKlzRhh4eLu7o67u3ueb5PXYmNjsbOzM3cYhdq+ffvYvHkzP//8M6+88gqAseVuTEyMmaMTQgghRHEmLQuFEEIIkS2enp4sWLCAqKgovvrqKwCuXLnC4MGDqVq1KnZ2dpQrV46uXbty5swZk22nTZvGuHHjAPD29jZ2c967d2+W95GRLVu2ULduXaytrfH29mb+/Plplrt8+TL9+vWjVKlSWFtb4+Pjw5IlS1LFqigK586do2/fvjg7O1O6dGlef/11IiIiUu3z4MGDtGnTBkdHR+zs7GjatClbtmwxKZNWF9p79+7xv//9jwoVKmBtbY27uzvNmjUzttxMa5vsxvbLL79Qu3ZtrK2tqVSpEosXL85y9+yUcidPnuSVV16hRIkSVK5cOcvXMSvnOGjQILy8vNI9dmbxpfeZyuy46cnqZzGn+wf4+eefKVGiBD169DAu279/P3fv3uXFF1/MdPsUn376KZs2bcqwjE6nS7VMr9dn+RhCCCGEKF6kZaEQQgghsq1Tp05otVr2798PwO3bt3F1dWX27Nm4u7vz4MEDvvvuOxo3bsypU6eoXr06AEOHDuXBgwd8/vnnBAYGGrvW1qxZk9OnT2dpH+n5448/6NatG02aNGHt2rXo9Xrmzp3L3bt3TcqdP3+epk2bGis9y5Qpw/bt2xk5ciTh4eFMnTrVpHyvXr3w9/dnyJAhnDlzhkmTJgGwYsUKY5l9+/bRtm1bateuzbfffou1tTVLly6la9eurFmzBn9//3Tjfu211zh58iSffPIJ1apV49GjR5w8eZL79+9n+u+Qldi2bdtGz549adGiBevWrUOn0zF//vxU1yUzPXv2pE+fPgwbNoyYmJhsXcdnOcfMZPSZ6t27d46Om9XP87Oc1+HDh2ncuLHxeNu3b2f8+PG0adOGLl26ZPn8jx8/zoQJE/jpp5/o3r27ybqTJ0/Sr18/Ll++zEsvvcTq1auJiorilVde4fTp09SvX58ff/yRqlWrZvl4QgghhCgGVCGEEEKIp6xcuVIF1KCgoHTLlC5dWvXx8UlznU6nUxMTE9WqVauqo0aNMlk3b948FVCDg4MzjCGjfaSlcePGatmyZdW4uDjjssjISLVkyZLqkz952rdvr5YvX16NiIgw2f6dd95RbWxs1AcPHqiqqqpTp05VAXXu3Lkm5YYPH67a2NioBoPBuOz5559XS5UqpUZFRZnE7+vrq5YvX95YNuW6PnnuDg4O6nvvvZfueaW1TXZia9iwoVqhQgU1ISHBuCwqKkp1dXVVs/JTMOVYU6ZMMVme1euYlXMcOHCgWrFixXSP/aS0rkd6n6nMjptV6X0Wc7r/uLg41cLCQp06dao6Y8YMFVAB1dPTU71x40a2Y+vXr59qaWmpbty40WRd7dq11W3btqn3799X/f391fbt26sNGzZUJ0yYoN6/f19dt26d2qxZs2zHL4QQQoiiTbohCyGEECJHVFU1/q3T6Zg5cyY1a9bEysoKCwsLrKysuHz5MhcuXMjS/p5lHzExMQQFBdGzZ09sbGyMyx0dHenatavxfXx8PH/88Qc9evTAzs4OnU5nfHXq1In4+HiOHDlisu+XX37Z5H3t2rWJj483zggdExPD0aNHeeWVV3BwcDCW02q1vPbaa9y8eZNLly6lG3ujRo0ICAjg448/5siRIyQlJWV+sbIR2/Hjx+nevTtWVlbGcg4ODibXJSt69epl/Du71/FZzvFZ5PS4Wf0s5nT/J0+eRKfT0ahRI/r378/27duZPn06UVFRtGjRgujoaADCw8NNZiVP62VhYcGPP/5IUlISvXv3NrYYvXnzJuXLl6d9+/aULFmS7777jitXrhhbTJYsWZLevXsDGI8nhBBCCAHSDVkIIYQQORATE8P9+/epVasWAKNHj2bJkiVMmDCBli1bUqJECTQaDUOHDiUuLi5L+3yWfTx8+BCDwUCZMmVSrXty2f3799HpdHz++ed8/vnnae4rPDzc5L2rq6vJe2trawBjTA8fPkRV1TRnKy5btqzxuOlZt24dH3/8Md988w2TJ0/GwcGBHj16MHfu3DTPJyexpTUJTXYnpnny/LJ7HZ/lHJ9FTo+b1c9iTvd/7NgxILmy0c3NjUqVKtGuXTuqVatG3759OXLkCC+99BKOjo4sX7480/Pctm0bGzZsoFu3bsbPhKqqaDT/tQuwtrbGw8OD+/fvYzAYTNY9WfEvhBBCCCGVhUIIIYTIti1btqDX62nVqhUAP/zwAwMGDGDmzJkm5cLDw3FxccnSPp9lHyVKlEBRFO7cuZNq3ZPLSpQoYWzx9/bbb6e5L29v7yzF++Q+NRoNoaGhqdbdvn0bADc3t3S3d3NzY9GiRSxatIiQkBA2b97MxIkTCQsLY9u2bdmKJa3YFEVJc3zCtK5VRp6caCS71zGzc7SxsSEhISHVPp6uuM2unF7brH4Wc7r/Y8eOUalSpXQ/FymVzNbW1gwdOjTDc9yyZQu//fYbr7zyCmvWrMHCIvnnffny5bl+/To7duygSZMmzJw5E0tLS+Li4hgxYgRz585l165d6PV6HB0dMzyGEEIIIYoXqSwUQgghRLaEhIQwduxYnJ2defPNN4HkiqSUVm0ptmzZwq1bt6hSpYrJ8qdbv6XIzj6eZm9vT6NGjQgMDGTevHnGrshRUVH8+uuvxnJ2dna0bt2aU6dOUbt2bZOuuTllb29P48aNCQwMZP78+dja2gJgMBj44YcfKF++PNWqVcvSvjw9PXnnnXf4448/OHToUK7E5ufnx6ZNm5g/f77xfKOjo/ntt99yvN9nuY5pnaOXlxdhYWHcvXvX2OIxMTGR7du3Z2mf6X2mMjtuenLyWczO/o8dO5aqJaqqqnzzzTf4+vpSs2bNDLd/0rx584wT6aRUFKacQ0BAAP379+fixYs8//zzBAYGkpiYSK9evXBwcKB69eqsX78+y8cSQgghRPEglYVCCCGESNfZs2eNY9GFhYVx4MABVq5ciVarZePGjbi7uwPQpUsXAgICqFGjBrVr1+bEiRPMmzeP8uXLp9pnStflxYsXM3DgQCwtLalevXq29pGWjz76iA4dOtC2bVvGjBmDXq9nzpw52Nvb8+DBA2O5xYsX88ILL9C8eXPeeustvLy8iIqK4sqVK/z666/s3r0729dp1qxZtG3bltatWzN
      "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}