diffusion_model_book_3_VE_SDE.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyOJ2slfa1Nttr8j+spvFY1x",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/shu65/diffusion-model-book/blob/main/diffusion_model_book_3_VE_SDE.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QEb0eFhy-8Fu"
      },
      "outputs": [],
      "source": [
        "from tqdm import tqdm_notebook as tqdm\n",
        "import torch\n",
        "\n",
        "device=\"cpu\""
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "n_samples = int(1e6)\n",
        "sigma = 0.01\n",
        "\n",
        "dist0 = torch.distributions.MultivariateNormal(torch.tensor([-2, -2], dtype=torch.float).to(device), sigma*torch.eye(2, dtype=torch.float).to(device))\n",
        "samples0 = dist0.sample(torch.Size([n_samples//2]))\n",
        "    \n",
        "dist1 = torch.distributions.MultivariateNormal(torch.tensor([2, 2], dtype=torch.float).to(device), sigma*torch.eye(2, dtype=torch.float).to(device))\n",
        "samples1 = dist1.sample(torch.Size([n_samples//2]))\n",
        "samples = torch.vstack((samples0, samples1))\n",
        "\n",
        "mean = torch.mean(samples, dim=0)\n",
        "std = torch.std(samples, dim=0)\n",
        "\n",
        "normalized_samples = (samples - mean[None, :])/std[None, :]\n",
        "\n",
        "dataset = torch.utils.data.TensorDataset((normalized_samples))"
      ],
      "metadata": {
        "id": "OUg5wTin-_yP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "%matplotlib inline\n",
        "import numpy as np\n",
        "import torch\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plot_samples = normalized_samples.cpu().numpy()\n",
        "\n",
        "plt.hist2d(plot_samples[:,0], plot_samples[:,1], range=((-2, 2), (-2, 2)), cmap='viridis', rasterized=False, bins=100, density=True)\n",
        "plt.gca().set_aspect('equal', adjustable='box')\n",
        "plt.xlim([-2, 2])\n",
        "plt.ylim([-2, 2])\n",
        "plt.title('Sample Density')\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 281
        },
        "id": "8SRAyXtK_Nj4",
        "outputId": "6a4687cf-fb91-4318-dbbb-6841227bf603"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAEICAYAAACnA7rCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAS4klEQVR4nO3de5RdZX3G8e8zQ0hEEkLMhZhEQI20eCuIBBa4pAUqpAreaKHeQF1ZWrFarZdq6+2PSletWpVVSsWKQAEXig0lComokWoogYZroKYUzYQkExIIYLhkZn79Y7+THA7n5J2wT/bZZ/J81pp19t7nnfO+GzLPvGfvd85PEYGZ2a70dXsAZlZ/Dgozy3JQmFmWg8LMshwUZpbloDCzLAeFtSXpc5Iu7fY4doek10i6t9vjGG8cFDUk6XhJv5C0VdIWSf8p6dXdHtfukHS/pMclPSrp4XQ+75O0R//NRcTPI+KwpnGctCf73Bs4KGpG0hTgP4CvA9OAOcDngSe7Oa5n6Q0RMRk4GDgP+ARwUXeHZM+Gg6J+XgIQEZdHxHBEPB4R10fE7QCSXiTpBkmbJT0o6TJJU0e/Of0G/Zik2yX9VtJFkmZJ+mH67b5M0oGp7SGSQtIiSQ9IWi/pL9sNTNIxaWbwsKTbJJ0wlhOKiK0RsRj4E+Bdkl6WXm+ipC9J+o2kjZIukPSc9NwJkgYkfVTSYBrbOQ1jWSjp7nRO60bHPfp9afsS4AXANZIek/RxSddK+mDTed0u6U1jOZe9VkT4q0ZfwBRgM3AxcCpwYNPzLwZOBiYCM4DlwFcbnr8fWAHMopiNDAK3AkcAk4AbgM+mtocAAVwOPBd4ObAJOCk9/zng0rQ9J41rIcUvmJPT/ow253H/6Os0Hf8N8P60/RVgMcXMaTJwDfDF9NwJwBDwBWBC6nfb6H8PYD3wmrR9IHBkw/cNtBsH8MfATQ37r0znsW+3/9/X+cszipqJiEeA4yl+gP8F2CRpsaRZ6fk1EbE0Ip6MiE3Al4HXNr3M1yNiY0SsA35O8YPx3xHxBHA1RWg0+nxE/DYi7gD+FTirxdDeDiyJiCURMRIRS4GVFD/Au+MBYJokAYuAv4iILRHxKPC3wJkNbbcDX4iI7RGxBHgMOKzhucMlTYmIhyLi1jH2vxh4iaT5af8dwJUR8dRunsdexUFRQxGxOiLOjoi5wMuA5wNfBUhvI65I0+1HgEuB6U0vsbFh+/EW+/s3tV/bsP3r1F+zg4Ez0tuOhyU9TBFos3fr5IqZyRaK2dB+wC0Nr/ejdHzU5ogYatjf1jD2t1CE1K8l/UzSsWPpPIXllcDb04XVs4BLdvMc9joOipqLiHuAb1MEBhS/dQN4eURMofhNr5LdzGvYfgHFb/1ma4FLImJqw9dzI+K8sXaS7tzMAW4EHqQIrZc2vN4BEdEcYi1FxM0RcTowE/gB8N12TVscuxh4G3AisC0ifjnWc9hbOShqRtLvpAt4c9P+PIrfeitSk8kUU/CtkuYAH+tAt38jaT9JLwXOofiN2+xS4A2SXiepX9KkdOFw7hjOaYqk1wNXUFzzuCMiRijeWn1F0szUbo6k143h9faV9DZJB0TEduARYKRN843ACxsPpGAYAf4BzybGxEFRP48CC4CbJP2WIiDuBD6anv88cCSwFbgW+H4H+vwZsAb4MfCliLi+uUFErAVOBz5FccFzLUVI7erf0DWSHk1tP01xPeWchuc/kfpdkd5GLWPnNYicdwD3p+97H8UMoZUvAn+d3t403tH5DsXF255aUNYtSld+bS8k6RDg/4AJTdcCxj1J7wQWRcTx3R5LL/CMwvY6kvYD/gy4sNtj6RWlg0LSPEk/SYtf7pL0oRZtJOlrktakxS1Hlu3X7NlI10A2UVy7+LcuD6dnlH7rIWk2MDsibpU0GbgFeGNE3N3QZiHwQYrbWQuAf4yIBaU6NrPKlJ5RRMT60cUuadHMaopbYI1OB74ThRXA1BQwZtYD9unki6WLY0cANzU9NYenL+oZSMfWt3iNRRQr9uin/1X7MaWTQzSzBo/y0IMRMSPXrmNBIWl/4HvAh9My5GclIi4kXWSaommxQCd2aIRm1mxZXPXrsbTryF0PSRMoQuKyiGh1X38dT1/9NzcdM7Me0Im7HqL4jIHVEfHlNs0WA+9Mdz+OAbZGxDPedphZPXXircdxFKvk7pC0Kh37FMXfDBARFwBLKO54rKH4w55znvkyZlZXpYMiIm4k80dJUdyD/UDZvsysO7wy08yyHBRmluWgMLMsB4WZZTkozCzLQWFmWQ4KM8tyUJhZloPCzLIcFGaW5aAwsywHhZllOSjMLMtBYWZZDgozy3JQmFmWg8LMshwUZpbVqU/h/pakQUl3tnn+BElbJa1KX5/pRL9mVo1O1fX4NvANilLy7fw8Il7fof7MrEIdmVFExHJgSydey8zqp8prFMdKuk3SDyW9tMJ+zaykjtYe3YVbgYMj4rFU2fwHwPxWDRtrj05iv4qGZ2a7UsmMIiIeiYjH0vYSYIKk6W3aXhgRR0XEUROYWMXwzCyjkqCQdFAqPYiko1O/m6vo28zK68hbD0mXAycA0yUNAJ8FJsCOkoJvBd4vaQh4HDgzVQ8zsx7QkaCIiLMyz3+D4vapmfUgr8w0sywHhZllOSjMLMtBYWZZDgozy3JQmFmWg8LMshwUZpbloDCzLAeFmWU5KMwsy0FhZlkOCjPLquoTrszqT8XvTfX3P+2xUQwPF49D26sbVw14RmFmWZ5RmKWZRN+k4qMX+2Y8D4DhmQfubDIyAkD/4MMAjGwqPqBt5MknqhplV3lGYWZZnlHYXm/0WkTfrJkAbDnu+QBsPHZkZ5thAXDQLycDMPXG4nhsGCwex/k1C88ozCyrqtqjkvQ1SWsk3S7pyE70a2bVqKr26KkUBX/mAwuAf0qPZl2nCcWPwfDMKQBseG1xC/TuP9r5edDbYgiAV/d9BIApvyoudCpd1PRbjzEYQ+3R04HvRGEFMFXS7E70bWZ7XlUXM+cAaxv2B9Kx9c0NXVLQuiX60oKrfYqLmM/Rzkp1w6NlaNRUjqZPlYyt22p3MdMlBc3qp6oZxTpgXsP+3HTMrPvSsux9NjwMwIzlxe3RQ0feu7NNmkgctCItztpUtB3ePlTJELutqhnFYuCd6e7HMcDWiHjG2w4zq6eqao8uARYCa4BtwDmd6NesE0ZGZwWDmwCY/tNid+o9U3e0UbpG0b9hQ/E9e8ndjlFV1R4N4AOd6MvMqucl3GZR3OUYebz4A68YeACAvjTDKA4WM4rRaxKjf26+t6jdXQ8zqx/PKMxGpZlFDI0+7h3XH8bCMwozy3JQmFmWg8LMshwUZpbloDCzLAeFmWU5KMwsy0FhZlkOCjPLclCYWZaDwsyyHBRmluWgMLMsB4WZZTkozCyrUyUFT5F0byoZ+MkWz58taZOkVenrva1ex8zqqfQH10jqB84HTqYo7HOzpMURcXdT0ysj4tyy/ZlZ9ToxozgaWBMR90XEU8AVFCUEzWyc6ERQtCsX2OwtqZL5VZLmtXgeKEoKSlopaeV2nuzA8MysrKouZl4DHBIRrwCWAhe3a+iSgmb104mgyJYLjIjNETE6Pfgm8KoO9GtmFelEUNwMzJd0qKR9gTMpSgjuIGl2w+5pwOoO9GtmFSl91yMihiSdC1wH9APfioi7JH0BWBkRi4E/l3QaMARsAc4u26+ZVUeRKiDV0RRNiwU6sdvDMBu3lsVVt0TEUbl2XplpZlkOCjPLclCYWZaDwsyyHBRmluWgMLMsB4WZZTkozCzLQWFmWQ4KM8tyUJhZloPCzLIcFGaW5aAwsywHhZllOSjMLMtBYWZZDgozy6qqpOBESVem52+SdEgn+jWzapQOioaSgqcChwNnSTq8qdl7gIci4sXAV4C/K9uvmVWnqpKCp7Oz6M9VwImS1IG+zawCVZUU3NEmIoaArcDzWr2YSwqa1U/tLma6pKBZ/VRSUrCxjaR9gAOAzR3o28wqUElJwbT/rrT9VuCGqHPlITN7mqpKCl4EXCJpDUVJwTPL9mtm1SkdFAARsQRY0nTsMw3bTwBndKIvM6te7S5mmln9OCjMLMtBYWZZDgozy3JQmFmWg8LMshwUZpbloDCzLAeFmWU5KMwsy0FhZlkOCjPLclCYWZaDwsyyHBRmluWgMLMsB4WZZTkozCyrVFBImiZpqaRfpccD27QblrQqfTV/8K6Z1VzZGcUngR9HxHzgx2m/lccj4vfS12kl+zSzipUNisZSgRcDbyz5emZWQ2WDYlZErE/bG4BZbdpNSmUCV0h6465e0CUFzeon+3H9kpYBB7V46tONOxERktoV9Tk4ItZJeiFwg6Q7IuJ/WzWMiAuBCwGmaJqLBJnVQDYoIuKkds9J2ihpdkSslzQbGGzzGuvS432SfgocAbQMCjOrn7JvPRpLBb4L+PfmBpIOlDQxbU8HjgPuLtmvmVWobFCcB5ws6VfASWkfSUdJ+mZq87vASkm3AT8BzosIB4VZDylVUjAiNgMntji+Enhv2v4F8PIy/ZhZd3llppllOSjMLMtBYWZZDgozy3JQmFmWg8LMshwUZpbloDCzLAeFmWU5KMwsy0FhZlkOCjPLclCYWVapvx61XdDODFZ/f/E4If3nlgCIp7YXj8PDxfEYqW58ZrvBMwozy/KMYg/pmzRx5/bM6QAMH1SUPYk0o9hncCsAI+s3Fo9PNHyYsGcXViOeUZhZlmcUHaZ9JgDQN+N5O45tOX4uABteO5zaFLOFGctnAzD9p8WHjcfAAzu+J4Y8o7D6KFtS8AxJd0kakXTULtqdIuleSWsktasmZmY1VXZGcSfwZuCf2zWQ1A+cD5wMDAA3S1o8Xj9gd/QOx8iMqTuObTy2mB3ctvBrAOzfV1y/eNHIIgCm3lu07RvctON7Ymj7nh6q2ZiV/XDd1QBKF+faOBpYExH3pbZXUJQiHJdBYTYeVXGNYg6wtmF/AFjQrrGkRcAigEnst2dHVhGNFEG6nWJmsW3kqeKJGH0+bYQLo1k9lSopGBHPKPhTlksKmtVPqZKCY7QOmNewPzcdM7MeUcVbj5uB+ZIOpQiIM4E/raDfrhhdjt2/8aEdxw76xf4AvFofKQ70RTpe3HTq37ABgOHtQ1UN02y3lL09+iZJA8CxwLWSrkvHny9pCUBEDAHnAtcBq4HvRsRd5YZtZlVS1PgC2hRNiwV6RsXCntA3cdLO7bT4amTW05dw9w8Ws46RTZuLRy/htooti6tuiYi2a6BGeQm3mWV5CfceMvLkEzu2Y8MgAEozB/UVM4rRaxJeXGV15xmFmWV5RlGB0RmDZw7WqzyjMLMsB4WZZTkozCzLQWFmWQ4KM8tyUJhZloPCzLIcFGaW5aAwsywHhZllOSjMLMtBYWZZDgozy3JQmFlWVSUF75d0h6RVklaW6dPMqrfHSwo2+P2IeLBkf2bWBVWUFDSzHlfVNYoArpd0SyoZaGY9pKqSgsdHxDpJM4Glku6JiOVt+ht3tUfNel0VJQWJiHXpcVDS1RQVzlsGhWuPmtXPHn/rIem5kiaPbgN/SHER1Mx6xB4vKQjMAm6UdBvwX8C1EfGjMv2aWbXK3vW4Gri6xfEHgIVp+z7glWX6MbPu8spMM8tyUJhZloPCzLIcFGaW5aAwsywHhZllOSjMLMtBYWZZDgozy3JQmFmWg8LMshwUZpbloDCzLAeFmWU5KMwsy0FhZlkOCjPLclCYWZaDwsyyyn647t9LukfS7ZKuljS1TbtTJN0raY2kT5bp08yqV3ZGsRR4WUS8Avgf4K+aG0jqB84HTgUOB86SdHjJfs2sQqWCIiKuj4ihtLsCmNui2dHAmoi4LyKeAq4ATi/Tr5lVq2w180bvBq5scXwOsLZhfwBY0O5FGksKAk8ui6vGY7Gg6cB4rOw+Xs8Lxu+5HTaWRh2pPSrp08AQcNnujLCVxpKCklZGxFFlX7NufF69Z7yem6SVY2lXuvaopLOB1wMnRkSrWqHrgHkN+3PTMTPrEWXvepwCfBw4LSK2tWl2MzBf0qGS9gXOBBaX6dfMqlX2rsc3gMnAUkmrJF0AT689mi52ngtcB6wGvhsRd43x9S8sOb668nn1nvF6bmM6L7V+t2BmtpNXZppZloPCzLJqHRRjXSLeiySdIekuSSOSev6223hdpi/pW5IGJY2r9TyS5kn6iaS707/DD+2qfa2DgjEsEe9hdwJvBpZ3eyBljfNl+t8GTun2IPaAIeCjEXE4cAzwgV39P6t1UIxxiXhPiojVEXFvt8fRIeN2mX5ELAe2dHscnRYR6yPi1rT9KMUdyTnt2tc6KJq8G/hhtwdhLbVapt/2H53Vi6RDgCOAm9q16eTfejwrVS8Rr9JYzs2smyTtD3wP+HBEPNKuXdeDogNLxGsrd27jiJfp9yBJEyhC4rKI+P6u2tb6rccYl4hb93mZfo+RJOAiYHVEfDnXvtZBQZsl4uOBpDdJGgCOBa6VdF23x/RslVymX2uSLgd+CRwmaUDSe7o9pg45DngH8AfpZ2uVpIXtGnsJt5ll1X1GYWY14KAwsywHhZllOSjMLMtBYWZZDgozy3JQmFnW/wMQi1L6azDj+gAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "sigma_min = torch.tensor(0.01, device=device)\n",
        "sigma_max = torch.tensor(10.0, device=device) \n",
        "\n",
        "def sigma(t, sigma_min=sigma_min, sigma_max=sigma_max):\n",
        "  return sigma_min * (sigma_max / sigma_min) ** t "
      ],
      "metadata": {
        "id": "m0NnI2oi_ZK5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plot_t = torch.linspace(0, 1, 1000, device=device)\n",
        "plot_log_sigma = sigma(plot_t)\n",
        "\n",
        "plt.plot(plot_t.cpu().numpy(), plot_log_sigma.cpu().numpy())\n",
        "plt.show()\n",
        "plt.plot(plot_t.cpu().numpy(), np.log(plot_log_sigma.cpu().numpy()))\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 513
        },
        "id": "8Up4sG3lAYEf",
        "outputId": "67edbeae-42f9-493c-86d3-b6baee99e2d5"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcFklEQVR4nO3de3zU9Z3v8dcn9wshISQQCISAchFQBCOC1lYL9dZW3banR1rd2nrW7f2yPdvTy+l6Tj17zp5Hu93a1t2Wdq11W7XV2hVrVy2oaKkiQUBuEm4BQkLuCblPZua7f8xAYyQQMpffzOT9fDzymN9t5vf5MvDml+985/sz5xwiIpJ80rwuQERExkYBLiKSpBTgIiJJSgEuIpKkFOAiIkkqI54nKykpcZWVlfE8pYhI0tu6dWuLc650+Pa4BnhlZSXV1dXxPKWISNIzsyNn2q4uFBGRJKUAFxFJUgpwEZEkpQAXEUlSCnARkSR1zgA3swfMrMnMdg3ZVmxmfzCz/eHHSbEtU0REhhvNFfiDwA3Dtn0V2OCcmwtsCK+LiEgcnTPAnXMvAW3DNt8C/Dy8/HPg1uiWJSKSGg42d/Pd5/bRdLI/6q891j7wqc65hvDyCWDqSAea2d1mVm1m1c3NzWM8nYhIcnq5ppnvP3+AwWD0770Q8YeYLnRHiBErc86tdc5VOeeqSkvf9k1QEZGUVn2knWmFOZQX5Ub9tcca4I1mNg0g/NgUvZJERFLH1iPtXDYrNuM8xhrg64CPhZc/BjwZnXJERFJHfUcfDZ39VHkV4Gb2CPAKMN/M6szsLuAfgPeY2X5gdXhdRESGqD7SDsBls4pj8vrnnI3QObdmhF2rolyLiEhK2VrbRl5WOhdNK4jJ6+ubmCIiMVJ9pJ1LZxaRkR6bqFWAi4jEQPeAn70NJ2PW/w0KcBGRmNh+tIOgg8sqY9P/DQpwEZGY2HqkHTNYWlEUs3MowEVEYqD6SBvzpxYwMSczZudQgIuIRFkg6Nh2tCNmX+A5RQEuIhJlextO0j3gp6pSAS4iklQ2Hw5N4HrF7MkxPY8CXEQkyl473MrM4lymx2ACq6EU4CIiURQMOl473Bbzq29QgIuIRNX+pm7aewdZPjt2479PUYCLiETR5sOtAKzQFbiISHLZfLiNaYU5zCyObf83KMBFRKLGOcfmQ21cMbsYM4v5+RTgIiJRcqilh5buAa6YE/vuE1CAi4hEzeZDofHf8fgAExTgIiJR89rhVkomZDOnJD8u51OAi4hEgXOOzYfbuGJOfPq/QQEuIhIVR9t6aejsZ0Wcuk9AAS4iEhWbDoTGf6+8oCRu51SAi4hEwaaDLUydmM0FpfHp/wYFuIhIxIJBxysHW7nqgpK49X+DAlxEJGJvnuiircfHVRfGr/sEFOAiIhHbdKAFQAEuIpJsNh1sYU5pPmWFOXE9rwJcRCQCPn+Q1w63cVUcR5+cogAXEYnAjroOen0BrrowPvOfDKUAFxGJwB/3t2AGK+foClxEJKn86WALF5cXUpiXGfdzK8BFRMaoZ8DPtqMdXOlB/zcowEVExuzVQ634g453xHn44CkRBbiZfcnMdpvZLjN7xMziO4ZGRMRDG2uayc1Mp6pykifnH3OAm1k58Hmgyjm3GEgHbotWYSIiiW5jTTMrL5hMTma6J+ePtAslA8g1swwgD6iPvCQRkcRX29LDkdZe3jWv1LMaxhzgzrnjwHeAo0AD0Omce274cWZ2t5lVm1l1c3Pz2CsVEUkgG2tCeXbN/CQMcDObBNwCzAamA/lmdvvw45xza51zVc65qtJS7xoqIhJNG2uaqZycx6zJ8Zs+drhIulBWA4edc83OuUHgCeDK6JQlIpK4+gcDvHKw1dPuE4gswI8CK8wsz0IT4K4C9kanLBGRxFVd207fYIB3edh9ApH1gW8GHgdeB3aGX2ttlOoSEUlYG2uayEpPY8Wc+M9/MlRGJE92zt0D3BOlWkREksLGmmaWzy4mLyuiCI2YvokpInIe6jv6qGns9rz/GxTgIiLnZcObTYC3wwdPUYCLiJyHDXsbmTU5jwunTPC6FAW4iMho9fr8/OlgK6sWTI3r3edHogAXERmll/e34PMHWX3RFK9LARTgIiKjtn5PIwU5GVw+u9jrUgAFuIjIqASDjhf2NXHN/ClkpidGdCZGFSIiCW57XQct3b6E6T4BBbiIyKhs2NtIeppxzTwFuIhIUlm/p4nLKyd5cvPikSjARUTO4VhbL/sau1h90VSvS3kLBbiIyDk8u/sEgAJcRCTZPLPrBAvKCqgs8e7mDWeiABcROYumk/1sPdrOjYuneV3K2yjARUTO4tk9jTgHNywu87qUt1GAi4icxTO7GphTks+8qd5PXjWcAlxEZATtPT5ePdTGDYvLEmLyquEU4CIiI/jDnkYCQZeQ/d+gABcRGdEzu09QXpTL4vKJXpdyRgpwEZEz6Oof5I/7WxK2+wQU4CIiZ7RhbxO+QJAbE3D0ySkKcBGRM1i3o55phTksq5jkdSkjUoCLiAzT3uPjpZpm3r9kOmlpidl9AgpwEZG3eWb3CfxBx81LpntdylkpwEVEhlm3vZ7ZJfksmp6Yo09OUYCLiAzReLKfVw+38v4l0xN29MkpCnARkSGefqMB50j47hNQgIuIvMW6HfUsnDaRC6ck3twnwynARUTCjrb2sv1YBzdfmvhX36AAFxE57cntxwF43yWJOffJcBEFuJkVmdnjZvamme01s5XRKkxEJJ6cc/zm9TpWzClmxqQ8r8sZlUivwO8DnnHOLQCWAHsjL0lEJP62HmmntrWXD1020+tSRi1jrE80s0LgncCdAM45H+CLTlkiIvH1+NY68rLSE3ruk+EiuQKfDTQDPzOzbWb2UzN72x0/zexuM6s2s+rm5uYITiciEht9vgBPv9HADYvLyM8e83Vt3EUS4BnAMuBfnHNLgR7gq8MPcs6tdc5VOeeqSktLIzidiEhsPLfnBF0Dfj502QyvSzkvkQR4HVDnnNscXn+cUKCLiCSVx7fWUV6Uy4rZk70u5byMOcCdcyeAY2Y2P7xpFbAnKlWJiMTJic5+Nh1o4QPLyhN65sEzibSz53PAL80sCzgEfDzykkRE4ueJbXUEHXxwWXJ1n0CEAe6c2w5URacUEZH4cs7xWHUdVbMmUVnytjEYCU/fxBSRceuVg60cbulhzfIKr0sZEwW4iIxbv3ztKIW5mbw3Sb46P5wCXETGpZbuAZ7bfYIPLCsnJzPd63LGRAEuIuPSY9V1DAYcH70iObtPQAEuIuNQMOh45LWjLK8s5sIpBV6XM2YKcBEZdzYdbOFoWy8fSeKrb1CAi8g49PDmo0zKy+SGJJq46kwU4CIyrjSe7OcPexr54LIZSfvh5SkKcBEZV37x6hECznH7illelxIxBbiIjBv9gwEe3nyUd8+fkpTfvBxOAS4i48ZTO+pp7fHxiXfM9rqUqFCAi8i44JzjZ5tqmT+1gCsvSK5pY0eiABeRceG1w23saTjJnVdVYpZc08aORAEuIuPCzzbVUpSXya2XlntdStQowEUk5R1r6+W5PSdYs7yC3KzkHjo4lAJcRFLezzbVYmbckQJDB4dSgItISmvv8fHolqPcsmQ604tyvS4nqhTgIpLSHnrlCL2+AH/9rgu8LiXqFOAikrJ6fX4e/NNhVl80hfllyTvr4EgU4CKSsn615RjtvYN86prUu/oGBbiIpKjBQJCfvHSI5ZXFXDar2OtyYkIBLiIpad32euo7+1P26hsU4CKSggJBxz+/eIAFZQVcM7/U63JiRgEuIinnqR31HGzu4fOr5qbM1+bPRAEuIinFHwjy/Q37WVBWwA2LkvuOO+eiABeRlLJuRz2HWnr44uq5pKWl7tU3KMBFJIWcuvq+aNpErluY2lffoAAXkRTy223HqW3t5Uvj4OobFOAikiIGA0F+8PwBFpdP5D0Lp3pdTlwowEUkJTy65RhH23r50up5KT3yZCgFuIgkve4BP/etr2H57GLevWCK1+XETcQBbmbpZrbNzH4XjYJERM7XT146REu3j6/duGDcXH1DdK7AvwDsjcLriIict6aufn7y8iHee/E0llZM8rqcuIoowM1sBvBe4KfRKUdE5Pzct34/Pn+Qv71+vtelxF2kV+DfA74CBEc6wMzuNrNqM6tubm6O8HQiIn92sLmbR7cc46NXVFBZku91OXE35gA3s/cBTc65rWc7zjm31jlX5ZyrKi1N3UllRCT+/v7pveRmpvO5VXO9LsUTkVyBXwXcbGa1wKPAu83sF1GpSkTkHJ5/s5Hn32ziC6vmUjIh2+tyPDHmAHfOfc05N8M5VwncBjzvnLs9apWJiIxgwB/gW0/t4YLSfD52ZaXX5Xgmw+sCRETO17/+8TC1rb089InlZGWM36+zRCXAnXMvAi9G47VERM7mRGc/P3z+ANctnMo7543vz9XG739dIpKU/u/v9+IPOv7nexd6XYrnFOAikjQ21jSzbkc9n3rXBVRMzvO6HM8pwEUkKfT6/HzjtzuZU5rPp69N3RsVnw99iCkiSeF76/dT197Hr/96JdkZ6V6XkxB0BS4iCW/X8U5++vIh1iyvYPnsYq/LSRgKcBFJaP5AkK8+8QaTJ2Tz1RsXeF1OQlEXiogktB9tPMiu4yf5548uozA30+tyEoquwEUkYe2u7+S+Dft53yXTuOniaV6Xk3AU4CKSkAb8Af7mVzuYlJfFvbcs9rqchKQuFBFJSP/0h/3sa+ziZ3dezqT8LK/LSUi6AheRhFNd28balw6yZnkF146je1yeLwW4iCSUzt5BvvDodmZMyuMb773I63ISmrpQRCRhOOf4ym920NTVz+OfvJIJ2Yqos9EVuIgkjIdeOcKzuxv5HzcsYMnMIq/LSXgKcBFJCLuOd/L3T+/l3QumcNc7ZntdTlJQgIuI5zr7Bvnsw69TnJ/Fd/7LEszM65KSgjqYRMRTgaDji49uo669j0fvXkGxhgyOmq7ARcRT31tfwwv7mrnn5kVUVWqiqvOhABcRzzyzq4EfPH+A/1o1k9uvqPC6nKSjABcRT9Q0dvHlX+/g0plFfOvWRer3HgMFuIjEXXPXAJ94cAu5WRn86PbLdIOGMVKAi0hc9fr83PXzLbR2+3jgzirKCnO8LilpKcBFJG4CQcfnH9nOruOd/GDNUi6ZUeR1SUlNAS4icXPv7/awfm8j97x/EasXTvW6nKSnABeRuLj/hQM8+Kda/ts7ZvOxKyu9LiclKMBFJOZ+/qdavv3sPv5iaTlfv0kzDEaLAlxEYurxrXXcs2431y2cyrc/dAlpaRouGC0KcBGJmd/vbOArj+/g6rkl/OAjS8lIV+REk/40RSQmnn6jgc8/so2lFZP48R0a6x0LCnARibontx/nc4+8ztKKIh78+OXkZWnevFgYc4Cb2Uwze8HM9pjZbjP7QjQLE5Hk9Fj1Mb74q+1cMXsyD358OQU5mV6XlLIi+W/RD3zZOfe6mRUAW83sD865PVGqTUSSzL+9Uss3n9zN1XNLWHtHFblZ6jaJpTEHuHOuAWgIL3eZ2V6gHFCAi4wzzjm+89w+7n/hIKsvmsIPP7KMnEyFd6xFpWPKzCqBpcDmaLyeiCSPwUCQrz+xk8e21rFm+UzuvWWxRpvEScQBbmYTgN8AX3TOnTzD/ruBuwEqKjTfr0gq6Rnw85mHX+fFfc18cfVcvrBqrqaFjaOIAtzMMgmF9y+dc0+c6Rjn3FpgLUBVVZWL5HwikjiOtfXyVw9VU9PYxf/7wMWsWa4LtHgbc4Bb6L/ZfwX2Oue+G72SRCTRvXqolU//8nUGA0EeuPNyrpk/xeuSxqVIOqquAu4A3m1m28M/N0WpLhFJUL949Qi3/3QzRXmZPPmZqxTeHopkFMofAXV2iYwTfb4Af/fkLh7bWsc180v5/pqlTNQYb0/p61Eick77G7v49C9f50BzN5+99kK+9J55pGtSKs8pwEXkrB6rPsY3n9xFflYGD31iOVfPLfW6JAlTgIvIGbX3+Pjmk7v43RsNrJhTzPdvW8qUibp/ZSJRgIvI26zf08jXfruTjl4f//26eXzqmgvVZZKAFOAictrJ/kHufWoPj22tY0FZAQ9+/HIWTS/0uiwZgQJcRHDO8R+7TvC/n9pNc9cAn7n2Aj6/aq7m8E5wCnCRce5oay/ffHIXG2uaWThtIj++o4pLZxZ5XZaMggJcZJzqHwzwk5cO8cMXDpCRZvzd+xbylytnaSKqJKIAFxlngkHHkzuO8+1n9lHf2c+Ni8u45/2LKCvUCJNkowAXGUdeO9zG/3l6D2/UdbK4fCL/+OFLWXnBZK/LkjFSgIuMA2/UdfC99ft5/s0mphXm8N0PL+HWS8tJ09DApKYAF0lhO+s6uW9DDev3NlGYm8nfXj+fT1w1W7c6SxEKcJEU45xj65F2frTxEOv3NjIxJ4Mvv2ced15VqRsMpxgFuEiKCAQdz+0+wdqXD7HtaAeFuZl8afU8Pv6OSs0amKIU4CJJ7mT/IE9sreOBTbUcbeulojiPb92yiA9dNoO8LP0TT2V6d0WSkHOOHXWdPLz5CE/taKBvMMDSiiK+duMCrltUpnlLxgkFuEgS6ewd5Hc763l481F2158kLyudW5dO5yPLZ3HxDM1ZMt4owEUSXP9ggBfebOLftx/nhTeb8QWCXDRtIvfeuphbL52uDybHMQW4SALy+YNsPtzK73Y08PtdDXT1+yktyOaOlbO49dJyFpdPJHRfcRnPFOAiCaLX52fjvmae3X2CDW820dXvJz8rnesXl/EXS8tZOWey5imRt1CAi3joWFsvG2uaeXFfEy/vb2HAH2RSXibXLyrjuoVTuXpuqb50IyNSgIvEUfeAn1cOtvJSTTMv72+mtrUXgPKiXNYsr+D6RWVcXjlJV9oyKgpwkRjq6PVRXdvOlto2XqttY2ddJ/6gIy8rnZVzJnPnlZVcPa+UOSX56tOW86YAF4kS5xxHWnvZUdfBlto2thxuZ19jFwBZ6WlcMqOQu985h6vnlrJsVpHudiMRU4CLjIFzjrr2PnYe7+SNuk52Hu9gZ10nJ/v9AEzIzmDZrEm8f8k0Lq8sZsnMInIyFdgSXQpwkXM42T/I/sYu9p3opqaxizdPnGTfiS7aewcByEw3FpRN5H1LpnNJeSEXzyhk/tQC9WNLzCnARQjdpebEyX5qW3o43NpDbUsPB5q62Xeii/rO/tPH5WelM6+sgOsXlbG4vJBLZhQyv6xA3SHiCQW4jBu9Pj/1Hf3Ud/RR197HkdYeDrf0UNvaw5HWXgb8wdPHZmekMbskn+Wzi5lXVsD8qQXMLyugvChXHzZKwlCAS0ro9flpOjlAc/cAjSdDIV3f0c/xjr7wct/pLo9TsjLSmFWcR2VJPtfMn8KsyXnMnpxPZUk+ZRNzdLcaSXgKcElIgaCjs2+Q9l4fHb2DdPT6aO320dw9QNPJfpq7B2ju+vNPjy/wtteYkJ1BeVEu04tyuHRmEeWTcsProZ9pCmlJcgpwiQnnHH2DAbr7/XQN+Onu99M94Kcr/NjdP0hnn5+OvlBAt/f6aA8HdUfvICf7B3HuzK9dkJPBlIJsSguyuXhGEaUTQsuntk2ZmM30olzdxEBSXkQBbmY3APcB6cBPnXP/EJWqJOaccwwGHL5AkIHBAH2DAfoHA/T5gvSF1/t84W3h5T8fE6B3MEB/eFuvL0DPkHDu6h+ke8BPcIQAHqogO4Oi/EyKcrMoystkVnEeRXmZFOVlMSkvk0l5WRSGHyfnZ1FakK3heCJhYw5wM0sH7gfeA9QBW8xsnXNuT7SKiwfnHM6BA4LOEQyvh5Y5ve6GrA89xrnQr/tB5/AHHYGgwx8Yuh7EHwhvD771MbQcfPu+QJCAI/TcoCMQeOt+nz/IYCCIzx/EFwj/DN0WXh4I7x88vd+d3u8LBM/5Z3MmZpCXmU5uVjo5menkhpcLcjKYPCGPCdmZFORkMCE7gwnhx9Pr4W0F2Zmhx5wMMjXUTmTMIrkCXw4ccM4dAjCzR4FbgKgH+Nd/u5PNh1r/HJqEgzQ4PFgBhqwH3xrGLrzPDQvnZJGeZqSbkZWRRlZGGpnpFn5MIys9jexTyxlp5GdnnD4uK7w/a8j+rCHPfUsYnyGcT23PyQq9jkZhiCSGSAK8HDg2ZL0OuGL4QWZ2N3A3QEVFxdhOVJTLgmkTSTPDgDQjtGx2ejktDWDIuhHef2r5DM8JH2MjPCctvO9Mz0k7/bzQekaahQI27dRy2ultpx/TQ8/LSEs7vf6W/eHt6cOfF35UcIrIUDH/ENM5txZYC1BVVTWm693PXHthVGsSEUkFkXRAHgdmDlmfEd4mIiJxEEmAbwHmmtlsM8sCbgPWRacsERE5lzF3oTjn/Gb2WeBZQsMIH3DO7Y5aZSIiclYR9YE7534P/D5KtYiIyHnQIFwRkSSlABcRSVIKcBGRJKUAFxFJUubi+F1yM2sGjozx6SVASxTLSQZq8/igNqe+SNs7yzlXOnxjXAM8EmZW7Zyr8rqOeFKbxwe1OfXFqr3qQhERSVIKcBGRJJVMAb7W6wI8oDaPD2pz6otJe5OmD1xERN4qma7ARURkCAW4iEiSSrgAN7MbzGyfmR0ws6+eYX+2mf0qvH+zmVV6UGZUjaLNf2Nme8zsDTPbYGazvKgzms7V5iHHfdDMnJkl9ZCz0bTXzD4cfp93m9nD8a4x2kbx97rCzF4ws23hv9s3eVFnNJnZA2bWZGa7RthvZvb98J/JG2a2LKIThm7qmxg/hKalPQjMAbKAHcDCYcd8GvhRePk24Fde1x2HNl8L5IWXPzUe2hw+rgB4CXgVqPK67hi/x3OBbcCk8PoUr+uOQ5vXAp8KLy8Ear2uOwrtfiewDNg1wv6bgP8ADFgBbI7kfIl2BX76RsnOOR9w6kbJQ90C/Dy8/DiwypL7ZpHnbLNz7gXnXG949VVCdz9KZqN5nwHuBf4/0B/P4mJgNO39K+B+51w7gHOuKc41Rtto2uyAieHlQqA+jvXFhHPuJaDtLIfcAjzkQl4Fisxs2ljPl2gBfqYbJZePdIxzzg90ApPjUl1sjKbNQ91F6H/wZHbONod/tZzpnHs6noXFyGje43nAPDPbZGavmtkNcasuNkbT5v8F3G5mdYTuK/C5+JTmqfP9935WMb+psUSPmd0OVAHv8rqWWDKzNOC7wJ0elxJPGYS6Ua4h9BvWS2Z2sXOuw8uiYmwN8KBz7h/NbCXwb2a22DkX9LqwZJFoV+CjuVHy6WPMLIPQr16tcakuNkZ1c2gzWw18A7jZOTcQp9pi5VxtLgAWAy+aWS2hvsJ1SfxB5mje4zpgnXNu0Dl3GKghFOjJajRtvgv4NYBz7hUgh9CkT6ksqjeDT7QAH82NktcBHwsvfwh43oU/HUhS52yzmS0FfkwovJO9bxTO0WbnXKdzrsQ5V+mcqyTU73+zc67am3IjNpq/1/9O6OobMysh1KVyKI41Rtto2nwUWAVgZhcRCvDmuFYZf+uAvwyPRlkBdDrnGsb8al5/ajvCp7Q1hD7B/kZ427cI/QOG0Jv8GHAAeA2Y43XNcWjzeqAR2B7+Wed1zbFu87BjXySJR6GM8j02Qt1Ge4CdwG1e1xyHNi8ENhEaobIduM7rmqPQ5keABmCQ0G9VdwGfBD455H2+P/xnsjPSv9f6Kr2ISJJKtC4UEREZJQW4iEiSUoCLiCQpBbiISJJSgIuIJCkFuIhIklKAi4gkqf8E16yIKUDCap0AAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfuUlEQVR4nO3dd3RVZb7G8e9LJxBCCR1C6AESaqg6isooAoKA3rF3Ued6504TgqCiWLCN41hGsTOjoyOhCiii2FAQUDgpJBBC7zUJKaSc9/6RrLkMQwmcnbNPeT5rsdZp7P3snORhZ5+9fxhrLSIiEryquR1ARER8oyIXEQlyKnIRkSCnIhcRCXIqchGRIFfDjZVGR0fb2NhYN1YtIhK01q5de9Ba2/Tkx10p8tjYWNasWePGqkVEgpYxZtupHtehFRGRIKciFxEJcipyEZEgpyIXEQlyKnIRkSCnIhcRCXIqchGRIKciFxHxgyP5xTy6MI3cohLHl+3KBUEiIuHCWsvilL08siCVowUlXNAxmmHdmzu6DhW5iEgV2Z9bxNR5qSxN30dC6yj+dudAurVs4Ph6VOQiIg6z1vLxmp1MX5ROcamXyVfGceeF7alRvWqOZqvIRUQctP1QAZPneliRdYgB7Rvz9PietI+uV6XrVJGLiDigzGt59/utPPdZJtWrGR6/Op4bBsRQrZqp8nWryEVEfLRpXx4Tkz38vP0ol3RtyhNjE2jVsK7f1u9zkRtj2gKzgOaABWZaa1/0dbkiIoGuuNTLa19v5uUvs6hXuzp//lVvxvRuhTFVvxd+Iif2yEuBP1hrfzLGRAJrjTGfW2vTHVi2iEhA8uw8ysTZHjL25nFVr1Y8clV3ouvXdiWLz0Vurd0D7Km4nWeM2QC0BlTkIhJyCovL+POyjbzxbTZNI2vzxi2J/NLh88LPlaPHyI0xsUAfYNUpnpsATACIiYlxcrUiIn6xMvsQScketh4q4PoBbZk8ohsN6tR0O5ZzRW6MqQ8kA7+11uae/Ly1diYwEyAxMdE6tV4RkaqWV1TCjCUZvL9qOzGNI/jgroEM6RTtdqx/caTIjTE1KS/x9621c5xYpohIIPgyYx9T5qayL7eIuy5szx8u70rdWtXdjvVvnDhrxQBvARustX/yPZKIiPsO5xfz2MI05q3bTZfm9Xn1xiH0iWnkdqxTcmKP/ALgZiDFGLOu4rEHrbWLHVi2iIhfWWtZ6NnDtAVp5BWV8Nthnfn10E7UqhG4w2KdOGvlO8C/J02KiFSBvTlFTJ2XwrIN++nVtiHPjO9J1xaRbsc6K13ZKSJhz1rLh6t38OSiDZR4vUwd2Y3bL2hPdT9cXu8EFbmIhLVth/JJSk7hh+xDDO7QhBnjE2jXpGqHXDlNRS4iYanMa3lnxRaeW5pJzWrVeGpcAtf1b+v3y+udoCIXkbCTubd8yNX6HUcZ1q0Zj1+dQIuoOm7HOm8qchEJG8WlXl5ZnsWrX2XRoE5NXrq+D6N6tgzKvfATqchFJCys23GUibPXs3HfMa7u3YqHr+pB43q13I7lCBW5iIS0wuIynl+aydsrttC8QR3evi2RS+PcHXLlNBW5iISs7zcfJCk5he2HC7hxYAxJV8YRGQBDrpymIheRkJNbVMJTizfwjx93ENskgg8nDGJQhyZux6oyKnIRCSmfp+9j6rwUDuQd556LOvDbYV0CbsiV01TkIhISDh47zrQFaXzi2UNci0jeuCWRnm0auh3LL1TkIhLUrLXMX7ebRxemkX+8jD/8sgv3XNwxoIdcOU1FLiJBa/fRQqbOS+XLjP30iSkfctW5eeAPuXKailxEgo7Xa/ngx+3MWJJBmdfy8Kju3DokNmiGXDlNRS4iQWXLwXySkj2s2nKYCztF89S4BNo2jnA7lqtU5CISFErLvLz13Rb+9PlGatWoxjPje3JtYpugv7zeCSpyEQl46btzmZTsIWVXDpd3b870q+Np3iB4h1w5TUUuIgHreGkZL3+ZxV+/2kzDiJq8ckNfRiS00F74SVTkIhKQ1m47wqRkD1n7jzGub2seGtmdRiEy5MppjhS5MeZtYBSw31ob78QyRSQ8FRSX8uxnmbz7/VZaNqjDO7f355KuzdyOFdCc2iN/F3gZmOXQ8kQkDH236SBJczzsPFLILYPbMXF4HPVr68DB2TjyFbLWfmOMiXViWSISfnIKSnhicTr/XLOTDtH1+Oc9gxnQvrHbsYKG/qkTEVd9mrqXh+ancji/mPuGduR/L+tMnZqhPeTKaX4rcmPMBGACQExMjL9WKyIB6kBe+ZCrRSl76N6yAe/c1p/41lFuxwpKfitya+1MYCZAYmKi9dd6RSSwWGuZ89MuHvskncLiMh64oisTLupAzerhM+TKaTq0IiJ+s+toIQ/OSeHrjQfo164RT4/vSadm9d2OFfScOv3wH8BQINoYsxN4xFr7lhPLFpHg5/Va/r5qG08vycACj47uwc2D2lEtTIdcOc2ps1aud2I5IhJ6Nh84RlKyh9Vbj/CLztE8OVZDrpymQysiUiVKyry88W02f162ibo1q/Pctb0Y37e1Lq+vAipyEXFc6q4cJiV7SNudy5XxLXh0TA+aRWrIVVVRkYuIY4pKynjpy0289nU2jSJq8dcb+3JlQku3Y4U8FbmIOGLN1sNMTPaQfSCfa/q1YerIbjSM0JArf1CRi4hPjh0v5dlPM5i1chutouoy644BXNSlqduxwoqKXETO29cbD/DgnBR25xRy6+BYHriiK/U05Mrv9BUXkXN2tKCY6Z9sIPmnnXRsWo+P7xlMYqyGXLlFRS4i52RJyh4emp/GkYJi7r+kE/df2klDrlymIheRStmfW8TD89P4NG0vPVo14L07+tOjlYZcBQIVuYickbWW2Wt3Mv2TdIpKvUwaHsfdv2hPDQ25ChgqchE5rR2HC3hwbgrfbjrIgNjGzBifQIemGnIVaFTkIvIfyryWWT9s5dnPMjHA9DE9uHGghlwFKhW5iPybrP15TEpOYe22I1zcpSlPjkugdcO6bseSM1CRiwhQPuTq9a8385cvsoioXZ0//VcvxvbRkKtgoCIXEVJ35fDAbA8b9uQysmdLpl3Vg6aRtd2OJZWkIhcJY0UlZfx52Sbe+DabJvVq8frN/biiRwu3Y8k5UpGLhKlV2YdImpPCloP5/CqxLQ+O7EZU3Zpux5LzoCIXCTN5RSU882kmf1u5jbaN6/L+XQO5oFO027HEBypykTCyPHM/U+aksCe3iDsuaM8fr+hCRC3VQLDTOygSBo7kFzP9k3Tm/LyLzs3qk3zfEPrGNHI7ljjEkSI3xgwHXgSqA29aa2c4sVwR8Y21lkUpe3hkfho5hSX85tJO/PelnahdQ0OuQonPRW6MqQ68AvwS2AmsNsYssNam+7psETl/+3KLmDovlc/T99GzTRR/v2sg3Vo2cDuWVAEn9sgHAFnW2mwAY8yHwBhARS7iAmst/1yzg8cXbaC41MuDI+K44wINuQplThR5a2DHCfd3AgNPfpExZgIwASAmJsaB1YrIybYfKiBpjofvNx9iYPvGPD2+J7HR9dyOJVXMbx92WmtnAjMBEhMTrb/WKxIOyryWd7/fynOfZVK9muGJsfFc3z9GQ67ChBNFvgtoe8L9NhWPiYgfbNyXx8TZHtbtOMqlcc14Ymw8LaM05CqcOFHkq4HOxpj2lBf4dcANDixXRM6guNTLX7/azMvLN1G/dg1evK43o3u10pCrMORzkVtrS40x9wOfUX764dvW2jSfk4nIaa3fcZRJyR4y9uYxulcrHrmqO03qa8hVuHLkGLm1djGw2IllicjpFRaX8cKyjbz5bTbNIuvw5i2JDOve3O1Y4jJd2SkSJH7YfIjJczxsPVTA9QNimDwijgZ1NORKVOQiAS+3qIQZSzL4YNV22jWJ4IO7BzKko4Zcyf9TkYsEsC827GPK3FT25xVx9y/a8/tfdqVuLV1eL/9ORS4SgA4dO86jC9NZsH43XZtH8trN/ejdtqHbsSRAqchFAoi1lgXrd/PownTyikr43bAu3De0I7Vq6PJ6OT0VuUiA2JNTyNS5qXyRsZ9ebRvyzPiedG0R6XYsCQIqchGXeb2WD1fv4KnFGyjxepk6shu3X9Ce6rq8XipJRS7ioq0H80ma42Fl9mEGd2jCjPEJtGuiIVdyblTkIi4oLfPyzoqtPP95JjWrVWPGuAR+1b+tLq+X86IiF/GzjL25TJrtYf3OHIZ1a87jV8fTIqqO27EkiKnIRfzkeGkZryzfzKvLs4iqW5OXru/DqJ4ttRcuPlORi/jBz9uPMCnZw8Z9xxjbpzUPjepO43q13I4lIUJFLlKFCopLeX7pRt5esYUWDerw9m2JXBqnIVfiLBW5SBX5PusgSXNS2H64gJsGxTBpeByRGnIlVUBFLuKwnMISnlq8gQ9X7yC2SQQfThjEoA5N3I4lIUxFLuKgpWl7mTovlYPHjnPPxR343bAu1KmpIVdStVTkIg44eOw40xak8YlnD3EtInnz1kR6tmnodiwJEypyER9Ya5m3bhePLkyn4HgZf/hlF+4d2pGa1TXkSvxHRS5ynnYfLWTK3BSWZx6gT0z5kKvOzTXkSvxPRS5yjrxey/s/bufpJRmUeS0Pj+rOrUNiNeRKXONTkRtjrgWmAd2AAdbaNU6EEglU2QeOkZScwo9bD3Nhp2ieGpdA28YRbseSMOfrHnkqMA543YEsIgGrtMzLm99t4YXPN1K7RjWeuaYn1/Zro8vrJSD4VOTW2g2AvpklpKXvzmVi8npSd+VyRY/mTB8TT7MGGnIlgcNvx8iNMROACQAxMTH+Wq3IeTteWsbLX2bx16820zCiJq/e2Jcr41tox0UCzlmL3BizDGhxiqemWGvnV3ZF1tqZwEyAxMREW+mEIi5Yu618yFXW/mOM69uah0Z2p5GGXEmAOmuRW2uH+SOISCDIP17Kc0szeff7rbSKqsu7t/dnaNdmbscSOSOdfihS4dtNB5g8J4WdRwq5dXA7HhgeR/3a+hGRwOfr6YdjgZeApsAiY8w6a+0VjiQT8ZOcghIeX5TOx2t30qFpPT6+dzD9Yxu7HUuk0nw9a2UuMNehLCJ+92nqXh6an8rh/GJ+PbQjv7mss4ZcSdDR740SlvbnFTFtQRqLU/bSvWUD3rmtP/Gto9yOJXJeVOQSVqy1JP+0i+mfpFNYUsYDV3RlwkUdNORKgpqKXMLGziMFPDg3lW82HqBfu0Y8Pb4nnZrVdzuWiM9U5BLyvF7L31Zu4+lPMwB4dHQPbh7UjmoaciUhQkUuIW3zgWNMmu1hzbYjXNSlKU+OjadNIw25ktCiIpeQVFLmZeY32bz4xSbq1qzOc9f2Ynzf1rq8XkKSilxCTuquHCYle0jbncuIhBZMG92DZpEaciWhS0UuIaOopIy/fLGJ17/JplFELV67qS/D41u6HUukyqnIJSSs3nqYSbM9ZB/M59p+bZg6sjtRETXdjiXiFypyCWrHjpfyzKcZzPphG20a1WXWHQO4qEtTt2OJ+JWKXILW1xsP8OCcFHbnFHLbkFgeuKIr9TTkSsKQvusl6BwtKOaxT9KZ89MuOjatx+x7B9OvnYZcSfhSkUvQsNayJHUvD89P5WhBCfdf0on7L+2kIVcS9lTkEhT25xbx0PxUPkvbR3zrBrx3xwB6tNKQKxFQkUuAs9by8dqdPP5JOsdLvSRdGcddF7anhoZcifyLilwC1o7DBUyek8J3WQcZENuYGeMT6NBUQ65ETqYil4BT5rXM+mErz3yaSTUD06+O58YBMRpyJXIaKnIJKFn785g428NP248ytGtTnhibQOuGdd2OJRLQVOQSEErKvLz21WZe+jKLiNrVeeFXvbi6t4ZciVSGr//58rPAVUAxsBm43Vp71IFcEkZSdubwwOz1ZOzNY2TPljw6ugfR9Wu7HUskaPj60f/nQLy1tiewEZjseyQJF0UlZTy1ZANjXvmOw/nFvH5zP165oa9KXOQc+bRHbq1desLdlcA1vsWRcLEq+xBJc1LYcjCf6/q3ZfKIbkTV1ZArkfPh5DHyO4CPTvekMWYCMAEgJibGwdVKMMkrKuHpTzP4+8rttG1cl/fvGsgFnaLdjiUS1M5a5MaYZUCLUzw1xVo7v+I1U4BS4P3TLcdaOxOYCZCYmGjPK60EteUZ+3lwbgp7c4u488L2/OHyLkTU0uftIr4660+RtXbYmZ43xtwGjAIus9aqoOU/HM4v5rGFacxbt5vOzeqTfN8Q+sY0cjuWSMjw9ayV4cBE4GJrbYEzkSRUWGv5xLOHaQvSyCks4TeXdea/L+lI7RoaciXiJF9/r30ZqA18XnG+70pr7b0+p5Kgty+3iClzU1m2YR8920Tx97sG0q1lA7djiYQkX89a6eRUEAkN1lo+Wr2DJxZvoLjUy5QR3bj9glgNuRKpQvqkSRyz/VABSXM8fL/5EAPbN+bp8T2Jja7ndiyRkKciF5+VeS3vrNjCc0szqVGtGk+OTeC6/m015ErET1Tk4pPMvXlMTPawfsdRLo1rxhNj42kZpSFXIv6kIpfzUlzq5dWvsnhleRaRdWry4nW9Gd2rlYZcibhARS7nbP2Oo0yc7SFzXx5jerfi4VHdaaL5KCKuUZFLpRUWl/GnzzN567stNIusw5u3JDKse3O3Y4mEPRW5VMoPmw+RNMfDtkMF3DAwhqQr42hQR0OuRAKBilzOKLeohKcWZ/CPH7fTrkkEH9w9kCEdNeRKJJCoyOW0lqXvY8q8FA7kHWfCRR343bAu1K2ly+tFAo2KXP7DoWPHeXRhOgvW7yauRSQzb06kV9uGbscSkdNQkcu/WGtZsH430xakcex4Kb8b1oX7hnakVg1dXi8SyFTkAsCenEKmzk3li4z99G7bkGeu6UmX5pFuxxKRSlCRhzmv1/KP1dt5anEGpV4vU0d24/YL2lNdl9eLBA0VeRjbcjCfpGQPq7YcZkjHJswY15OYJhFuxxKRc6QiD0OlZV7eXrGF55dupFb1aswYl8Cv+rfV5fUiQUpFHmY27MllUrIHz84chnVrzuNXx9Miqo7bsUTEByryMHG8tIxXlm/m1eVZRNWtycs39GFkQkvthYuEABV5GPhp+xEmzfawaf8xxvZpzcOjutOoXi23Y4mIQ1TkIayguJTnl27k7RVbaNGgDu/c1p9L4pq5HUtEHKYiD1Ersg6SNMfDjsOF3DQohknD44jUkCuRkORTkRtjpgNjAC+wH7jNWrvbiWByfnIKS3hy0QY+WrOD9tH1+GjCIAZ2aOJ2LBGpQr7ukT9rrX0IwBjzG+Bh4F6fU8l5WZq2l6nzUjmUX8y9F3fkt8M6U6emhlyJhDqfitxam3vC3XqA9S2OnI8DeceZtjCNRZ49dGvZgLdu7U9Cmyi3Y4mIn/h8jNwY8wRwC5ADXHKG100AJgDExMT4ulqhfMjV3J938dgn6RQcL+OPl3fhnos7UrO6hlyJhBNj7Zl3oo0xy4AWp3hqirV2/gmvmwzUsdY+craVJiYm2jVr1pxrVjnBrqOFTJmbwleZB+gbUz7kqlMzDbkSCWXGmLXW2sSTHz/rHrm1dlgl1/E+sBg4a5HL+fN6Le+v2saMJRl4LTxyVXduGRyrIVciYczXs1Y6W2s3VdwdA2T4HklOJ/vAMZKSU/hx62Eu7BTNU+MSaNtYQ65Ewp2vx8hnGGO6Un764TZ0xkqVKC3z8sa3W3hh2Ubq1KjGM9f05Np+bXR5vYgAvp+1Mt6pIHJq6btzmZi8ntRduVzRoznTx8TTrIGGXInI/9OVnQGqqKSMl7/M4rWvN9MwohZ/vbEvVya0dDuWiAQgFXkAWrvtMBNne9h8IJ/xfdvw0KhuNIzQkCsROTUVeQDJP17Ks59l8t4PW2kVVZf37hjAxV2auh1LRAKcijxAfLPxAJPnpLA7p5BbBrXjgeFx1K+tt0dEzk5N4bKcghKmL0pn9tqddGhaj3/eM5j+sY3djiUiQURF7qJPU/fw0Pw0DucX8+uhHfnNZRpyJSLnTkXugv15RTwyP40lqXvp3rIB79zWn/jWGnIlIudHRe5H1lpmr93J44s2UFhSxgNXdGXCRR005EpEfKIi95Mdhwt4cG4K3246SGK7RswY35NOzeq7HUtEQoCKvIp5vZZZP2zlmc8yMcBjY3pw08B2VNOQKxFxiIq8CmXtP0ZSsoc1245wUZemPDk2njaNNORKRJylIq8CJWVeZn6TzYvLNlG3VnWev7YX4/q21pArEakSKnKHpe7KYeJsD+l7chmR0IJHR8fTNLK227FEJISpyB1SVFLGi19sYuY32TSuV4vXburL8HgNuRKRqqcid8DqrYeZNNtD9sF8ru3XhqkjuxMVUdPtWCISJlTkPjh2vJRnPs1g1g/baNOoLn+7cwC/6KwhVyLiXyry8/RV5n6mzE1ld04ht18Qyx8v70o9DbkSEReoec7Rkfxipi9KZ85Pu+jUrD6z7x1Cv3aN3I4lImFMRV5J1loWp+zlkQWpHC0o4X8u7cT9l3aidg0NuRIRdzlS5MaYPwDPAU2ttQedWGYg2Z9bxNR5qSxN30dC6yhm3TGQ7q0auB1LRARwoMiNMW2By4HtvscJLNZaPl6zk+mL0iku9ZJ0ZRx3XdieGhpyJSIBxIk98heAicB8B5YVMHYcLmDynBS+yzrIgPaNmTEugQ5NNeRKRAKPT0VujBkD7LLWrg+Vy8/LvJb3vt/Ks59lUr2a4fGr47lhQIyGXIlIwDprkRtjlgEtTvHUFOBByg+rnJUxZgIwASAmJuYcIvrPpn15TEz28PP2owzt2pQnxybQqmFdt2OJiJyRsdae3180JgH4AiioeKgNsBsYYK3de6a/m5iYaNesWXNe660KxaVeXvt6My9/mUW92tV55KoejOndSkOuRCSgGGPWWmsTT378vA+tWGtTgGYnrGArkBhsZ614dh5l4mwPGXvzGNWzJdNG9yC6voZciUjwCNvzyItKynjh84288W02TSNrM/Pmflze41RHkEREAptjRW6tjXVqWVVtZfYhkpI9bD1UwPUD2pJ0ZTei6mrIlYgEp7DaI88rKmHGkgzeX7WdmMYRfHDXQIZ0inY7loiIT8KmyL/M2MeUuansyy3irgvb8/vLuxBRK2w2X0RCWMg32eH8Yh5bmMa8dbvp3Kw+r943hD4xGnIlIqEjZIvcWstCzx6mLUgjt7CE/72sM7++pKOGXIlIyAnJIt+bUz7katmGffRqE8XTdw8kroWGXIlIaAqpIrfW8uHqHTy5aAMlXi9TRnTjjgvbU12X14tICAuZIt92KJ+k5BR+yD7EoA6NmTGuJ7HR9dyOJSJS5YK+yMu8lndWbOG5pZnUrFaNJ8cmcF3/thpyJSJhI6iLPHNv+ZCr9TuOcllcMx4fG0/LKA25EpHwEpRFXlzq5dWvsnhleRaRdWry4nW9Gd1LQ65EJDwFXZGv23GUSbM9ZO7LY0zvVjw8qjtNNORKRMJYUBX5S19s4oVlG2kWWYe3bk3ksm7N3Y4kIuK6oCrymCYRXDcghqQr42hQR0OuREQgyIp8TO/WjOnd2u0YIiIBRf8dvIhIkFORi4gEORW5iEiQU5GLiAQ5FbmISJBTkYuIBDkVuYhIkFORi4gEOWOt9f9KjTkAbDvPvx4NHHQwTjDQNocHbXN48GWb21lrm578oCtF7gtjzBprbaLbOfxJ2xwetM3hoSq2WYdWRESCnIpcRCTIBWORz3Q7gAu0zeFB2xweHN/moDtGLiIi/y4Y98hFROQEKnIRkSAXsEVujBlujMk0xmQZY5JO8XxtY8xHFc+vMsbEuhDTUZXY5t8bY9KNMR5jzBfGmHZu5HTS2bb5hNeNN8ZYY0xQn6pWme01xvxXxfucZoz5wN8ZnVaJ7+sYY8xyY8zPFd/bI9zI6SRjzNvGmP3GmNTTPG+MMX+p+Jp4jDF9fVqhtTbg/gDVgc1AB6AWsB7oftJrfg28VnH7OuAjt3P7YZsvASIqbt8XDttc8bpI4BtgJZDodu4qfo87Az8DjSruN3M7tx+2eSZwX8Xt7sBWt3M7sN0XAX2B1NM8PwJYAhhgELDKl/UF6h75ACDLWpttrS0GPgTGnPSaMcB7FbdnA5cZY4wfMzrtrNtsrV1urS2ouLsSaOPnjE6rzPsMMB14GijyZ7gqUJntvRt4xVp7BMBau9/PGZ1WmW22QIOK21HAbj/mqxLW2m+Aw2d4yRhgli23EmhojGl5vusL1CJvDew44f7OisdO+RprbSmQAzTxS7qqUZltPtGdlP+LHszOus0Vv3K2tdYu8mewKlKZ97gL0MUYs8IYs9IYM9xv6apGZbZ5GnCTMWYnsBj4H/9Ec9W5/ryfUVD958tSzhhzE5AIXOx2lqpkjKkG/Am4zeUo/lSD8sMrQyn/jesbY0yCtfaom6Gq2PXAu9ba540xg4G/GWPirbVet4MFi0DdI98FtD3hfpuKx075GmNMDcp/JTvkl3RVozLbjDFmGDAFGG2tPe6nbFXlbNscCcQDXxljtlJ+LHFBEH/gWZn3eCewwFpbYq3dAmykvNiDVWW2+U7gnwDW2h+AOpQPlgpllfp5r6xALfLVQGdjTHtjTC3KP8xccNJrFgC3Vty+BvjSVnyKEKTOus3GmD7A65SXeLAfO4WzbLO1NsdaG22tjbXWxlL+ucBoa+0ad+L6rDLf1/Mo3xvHGBNN+aGWbD9mdFpltnk7cBmAMaYb5UV+wK8p/W8BcEvF2SuDgBxr7Z7zXprbn+6e4VPfEZTvjWwGplQ89hjlP8hQ/mZ/DGQBPwId3M7sh21eBuwD1lX8WeB25qre5pNe+xVBfNZKJd9jQ/nhpHQgBbjO7cx+2ObuwArKz2hZB1zudmYHtvkfwB6ghPLfsu4E7gXuPeF9fqXia5Li6/e1LtEXEQlygXpoRUREKklFLiIS5FTkIiJBTkUuIhLkVOQiIkFORS4iEuRU5CIiQe7/ABvYWgecHO6/AAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def ve_sde_marginal_prob_statistics(x, t, sigma_min, sigma_max):\n",
        "  mean = x\n",
        "  std = sigma(t=t, sigma_min=sigma_min, sigma_max=sigma_max)\n",
        "  return mean, std"
      ],
      "metadata": {
        "id": "KxspcTDE6UUW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def ve_sde_drift(t, sigma_min, sigma_max):\n",
        "  drift = torch.zeros_like(t)\n",
        "  return drift\n",
        "\n",
        "def ve_sde_diffusion(t, sigma_min, sigma_max):\n",
        "  std = sigma(t=t, sigma_min=sigma_min, sigma_max=sigma_max)\n",
        "  diffusion = std * torch.sqrt(2 * (torch.log(sigma_max) - torch.log(sigma_min))) # (30)\n",
        "  return diffusion"
      ],
      "metadata": {
        "id": "MVd11yVV9aRA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def dsm_loss(score_model, samples, sigma_min, sigma_max):\n",
        "  eps = 1.0e-8\n",
        "  t = torch.distributions.uniform.Uniform(torch.tensor([eps], device=samples.device), torch.tensor([1], device=samples.device)).sample([samples.shape[0]]) \n",
        "  z = torch.randn_like(samples)\n",
        "  mean, std = ve_sde_marginal_prob_statistics(x=samples, t=t, sigma_min=sigma_min, sigma_max=sigma_max)\n",
        "  noise = z * std\n",
        "  perturbed_samples = mean + z * std\n",
        "  scores = score_model(perturbed_samples, t)\n",
        "  target = - 1 / (std ** 2) * noise\n",
        "  \n",
        "  target = target.view(target.shape[0], -1)\n",
        "  scores = scores.view(scores.shape[0], -1)\n",
        "  g = ve_sde_diffusion(t=t, sigma_min=sigma_min, sigma_max=sigma_max)\n",
        "  lmd = g ** 2\n",
        "  loss = torch.sqrt(((scores - target) ** 2).sum(dim=-1)) * lmd\n",
        "  return loss.mean()"
      ],
      "metadata": {
        "id": "S5IND5euA58Z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "\n",
        "class ScoreModel(nn.Module):\n",
        "  def __init__(self, sigma_min, sigma_max, n_channels=2):\n",
        "    super(ScoreModel, self).__init__()\n",
        "    self.sigma_min = sigma_min\n",
        "    self.sigma_max = sigma_max\n",
        "    self.model = nn.Sequential(\n",
        "        nn.Linear(n_channels, 2*n_channels),\n",
        "        nn.ELU(),\n",
        "        nn.Linear(2*n_channels, 16*n_channels),\n",
        "        nn.ELU(),\n",
        "        nn.Linear(16*n_channels, 2*n_channels),\n",
        "        nn.ELU(),\n",
        "        nn.Linear(2*n_channels, n_channels),\n",
        "    )\n",
        "\n",
        "  def forward(self, x, t):\n",
        "    y = self.model(x)\n",
        "    sigma_t = sigma(t=t, sigma_min=self.sigma_min, sigma_max=self.sigma_max)\n",
        "    return y/sigma_t"
      ],
      "metadata": {
        "id": "AgbS5rmk-_70"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import torch\n",
        "\n",
        "batch_size = 512\n",
        "n_steps = 100000\n",
        "\n",
        "dataloader = torch.utils.data.DataLoader(dataset, batch_size=512, shuffle=True, num_workers=0)\n",
        "dataloader_iter = iter(dataloader)\n",
        "\n",
        "score_model = ScoreModel(sigma_min=sigma_min, sigma_max=sigma_max).to(device)\n",
        "\n",
        "optimizer = torch.optim.Adam(score_model.parameters())\n",
        "lr_scheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=0.001, total_steps=n_steps)\n",
        "\n",
        "for i in range(n_steps):\n",
        "  try:\n",
        "    x = next(dataloader_iter)[0]\n",
        "  except StopIteration:\n",
        "    dataloader_iter = iter(dataloader)\n",
        "    x = next(dataloader_iter)[0]\n",
        "  x = x.to(device)\n",
        "\n",
        "  optimizer.zero_grad()\n",
        "  loss = dsm_loss(score_model, x, sigma_min=sigma_min, sigma_max=sigma_max)\n",
        "  loss.backward()\n",
        "  optimizer.step()\n",
        "  lr_scheduler.step()\n",
        "  if (i % 1000) == 0:\n",
        "    print(f\"{i} steps loss:{loss}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "zVbQT-Iv_el2",
        "outputId": "cff8cda6-5cbd-4aa3-bd70-cccdb955dcc7"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0 steps loss:1859.2049560546875\n",
            "1000 steps loss:1944.5516357421875\n",
            "2000 steps loss:2025.0797119140625\n",
            "3000 steps loss:1731.193359375\n",
            "4000 steps loss:1460.9530029296875\n",
            "5000 steps loss:1305.2705078125\n",
            "6000 steps loss:1971.225341796875\n",
            "7000 steps loss:1768.1434326171875\n",
            "8000 steps loss:1572.624755859375\n",
            "9000 steps loss:1639.954345703125\n",
            "10000 steps loss:1911.223876953125\n",
            "11000 steps loss:1662.96435546875\n",
            "12000 steps loss:1700.100341796875\n",
            "13000 steps loss:1971.902587890625\n",
            "14000 steps loss:1711.3748779296875\n",
            "15000 steps loss:1776.391845703125\n",
            "16000 steps loss:1484.6676025390625\n",
            "17000 steps loss:1919.545654296875\n",
            "18000 steps loss:1619.9324951171875\n",
            "19000 steps loss:1688.3319091796875\n",
            "20000 steps loss:1337.708740234375\n",
            "21000 steps loss:1704.1593017578125\n",
            "22000 steps loss:1349.9967041015625\n",
            "23000 steps loss:1664.606201171875\n",
            "24000 steps loss:1819.019287109375\n",
            "25000 steps loss:1611.3419189453125\n",
            "26000 steps loss:1434.19580078125\n",
            "27000 steps loss:1729.96923828125\n",
            "28000 steps loss:1770.2747802734375\n",
            "29000 steps loss:1499.20166015625\n",
            "30000 steps loss:1685.4002685546875\n",
            "31000 steps loss:1660.073974609375\n",
            "32000 steps loss:1777.1656494140625\n",
            "33000 steps loss:1501.560791015625\n",
            "34000 steps loss:1613.48681640625\n",
            "35000 steps loss:1724.016357421875\n",
            "36000 steps loss:1239.4561767578125\n",
            "37000 steps loss:1382.387451171875\n",
            "38000 steps loss:1752.8114013671875\n",
            "39000 steps loss:1778.257080078125\n",
            "40000 steps loss:1354.485595703125\n",
            "41000 steps loss:1895.1209716796875\n",
            "42000 steps loss:1586.7813720703125\n",
            "43000 steps loss:1223.6361083984375\n",
            "44000 steps loss:1156.4461669921875\n",
            "45000 steps loss:1445.104736328125\n",
            "46000 steps loss:1626.667724609375\n",
            "47000 steps loss:1659.5614013671875\n",
            "48000 steps loss:1508.3685302734375\n",
            "49000 steps loss:1189.2852783203125\n",
            "50000 steps loss:1416.574462890625\n",
            "51000 steps loss:1378.7445068359375\n",
            "52000 steps loss:1555.0146484375\n",
            "53000 steps loss:1449.140380859375\n",
            "54000 steps loss:1639.216064453125\n",
            "55000 steps loss:1417.510498046875\n",
            "56000 steps loss:1543.335205078125\n",
            "57000 steps loss:1600.5015869140625\n",
            "58000 steps loss:1736.684326171875\n",
            "59000 steps loss:1602.336669921875\n",
            "60000 steps loss:1650.90234375\n",
            "61000 steps loss:1914.06884765625\n",
            "62000 steps loss:1724.775390625\n",
            "63000 steps loss:1662.1180419921875\n",
            "64000 steps loss:1596.939453125\n",
            "65000 steps loss:1654.9027099609375\n",
            "66000 steps loss:1425.41015625\n",
            "67000 steps loss:1750.7001953125\n",
            "68000 steps loss:1770.9078369140625\n",
            "69000 steps loss:1326.52978515625\n",
            "70000 steps loss:1599.9464111328125\n",
            "71000 steps loss:1204.662109375\n",
            "72000 steps loss:1262.2047119140625\n",
            "73000 steps loss:1526.552978515625\n",
            "74000 steps loss:1537.0894775390625\n",
            "75000 steps loss:1506.67578125\n",
            "76000 steps loss:1488.7884521484375\n",
            "77000 steps loss:1979.8375244140625\n",
            "78000 steps loss:1668.7127685546875\n",
            "79000 steps loss:1291.64990234375\n",
            "80000 steps loss:1474.0928955078125\n",
            "81000 steps loss:1243.7567138671875\n",
            "82000 steps loss:1886.950439453125\n",
            "83000 steps loss:1724.3499755859375\n",
            "84000 steps loss:1750.46484375\n",
            "85000 steps loss:1466.8629150390625\n",
            "86000 steps loss:1885.470458984375\n",
            "87000 steps loss:1604.2086181640625\n",
            "88000 steps loss:1576.0745849609375\n",
            "89000 steps loss:1479.397216796875\n",
            "90000 steps loss:2024.5301513671875\n",
            "91000 steps loss:1208.362060546875\n",
            "92000 steps loss:1726.1832275390625\n",
            "93000 steps loss:1445.5631103515625\n",
            "94000 steps loss:1562.2728271484375\n",
            "95000 steps loss:1294.272216796875\n",
            "96000 steps loss:1663.67822265625\n",
            "97000 steps loss:1811.1575927734375\n",
            "98000 steps loss:1532.77783203125\n",
            "99000 steps loss:1943.4501953125\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def euler_maruyama_sample(n_samples, score_model, device=device, n=1000):\n",
        "  with torch.no_grad():\n",
        "    x = torch.randn(n_samples, 2, device=device)\n",
        "    dt = torch.tensor(1.0 / n, device=x.device)\n",
        "    for t in range(n, 0, -1):\n",
        "      if (t % 10) == 0:\n",
        "        print(f\"t:{t}\")\n",
        "      assert not torch.any(torch.isnan(x))\n",
        "      t_tensor = torch.full((n_samples, 1), t/n, device=device)\n",
        "      z = torch.randn(n_samples, 2)\n",
        "      f = ve_sde_drift(t_tensor, score_model.sigma_min, score_model.sigma_max)\n",
        "      g = ve_sde_diffusion(t_tensor, score_model.sigma_min, score_model.sigma_max)\n",
        "      g2 = g ** 2\n",
        "      score = score_model(x, t_tensor)\n",
        "      x = x - (f*x - g2 * score) * dt + g * torch.sqrt(dt) * z\n",
        "    return x"
      ],
      "metadata": {
        "id": "cZdfLqp8AoWC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "samples_pred = euler_maruyama_sample(n_samples=100000, score_model=score_model)\n",
        "samples_pred = samples_pred.cpu().numpy()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HbUOnUjbArk3",
        "outputId": "e08f551a-c2f4-48a6-fbfb-a400b6966bf9"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "t:1000\n",
            "t:990\n",
            "t:980\n",
            "t:970\n",
            "t:960\n",
            "t:950\n",
            "t:940\n",
            "t:930\n",
            "t:920\n",
            "t:910\n",
            "t:900\n",
            "t:890\n",
            "t:880\n",
            "t:870\n",
            "t:860\n",
            "t:850\n",
            "t:840\n",
            "t:830\n",
            "t:820\n",
            "t:810\n",
            "t:800\n",
            "t:790\n",
            "t:780\n",
            "t:770\n",
            "t:760\n",
            "t:750\n",
            "t:740\n",
            "t:730\n",
            "t:720\n",
            "t:710\n",
            "t:700\n",
            "t:690\n",
            "t:680\n",
            "t:670\n",
            "t:660\n",
            "t:650\n",
            "t:640\n",
            "t:630\n",
            "t:620\n",
            "t:610\n",
            "t:600\n",
            "t:590\n",
            "t:580\n",
            "t:570\n",
            "t:560\n",
            "t:550\n",
            "t:540\n",
            "t:530\n",
            "t:520\n",
            "t:510\n",
            "t:500\n",
            "t:490\n",
            "t:480\n",
            "t:470\n",
            "t:460\n",
            "t:450\n",
            "t:440\n",
            "t:430\n",
            "t:420\n",
            "t:410\n",
            "t:400\n",
            "t:390\n",
            "t:380\n",
            "t:370\n",
            "t:360\n",
            "t:350\n",
            "t:340\n",
            "t:330\n",
            "t:320\n",
            "t:310\n",
            "t:300\n",
            "t:290\n",
            "t:280\n",
            "t:270\n",
            "t:260\n",
            "t:250\n",
            "t:240\n",
            "t:230\n",
            "t:220\n",
            "t:210\n",
            "t:200\n",
            "t:190\n",
            "t:180\n",
            "t:170\n",
            "t:160\n",
            "t:150\n",
            "t:140\n",
            "t:130\n",
            "t:120\n",
            "t:110\n",
            "t:100\n",
            "t:90\n",
            "t:80\n",
            "t:70\n",
            "t:60\n",
            "t:50\n",
            "t:40\n",
            "t:30\n",
            "t:20\n",
            "t:10\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "%matplotlib inline\n",
        "import numpy as np\n",
        "import torch\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.hist2d(samples_pred[:,0], samples_pred[:,1], range=((-2, 2), (-2, 2)), cmap='viridis', rasterized=False, bins=100, density=True)\n",
        "plt.gca().set_aspect('equal', adjustable='box')\n",
        "plt.xlim([-2, 2])\n",
        "plt.ylim([-2, 2])\n",
        "plt.title('Predicted Sample Density')\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 281
        },
        "id": "qup4YWWkIcIr",
        "outputId": "dac3201c-1d72-4775-a885-5cde6be0cb5f"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAEICAYAAACnA7rCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXXElEQVR4nO3de7ScVXnH8e/vnBwSciPGIGASQYVS8YpSAkvbWtEWqAVrtcXWC1XrqpWl7XK1tbWrrra6qst6qQsq4qVqaxULoqgggkrRKkqg8QIhJVIxCeEWEkhISHLOefrH3nPynsnM2SeZOe/MnPw+a80672XPu/fMmXnmefe7Z7YiAjOzqQz1ugFm1v8cKMysyIHCzIocKMysyIHCzIocKMysyIGiRyR9UtI78/IvS1pXU70h6fg66jpYks6X9J1et+NASHqCpB2ShnvdlpngQDEFST+TtCu/AO7Nb+6F3a4nIr4dESdOoz0z+gaS9FRJX5f0oKRtkm6WdPZM1TcTJF0v6VFJ2yU9nB/D2yTNncl6I+LnEbEwIsYq7Xj9TNZZJweKst+KiIXAs4FTgL9pLiBpTu2tmhlfBq4FjgYeB7wZeLinLTo4F0TEIuAY4K3AecBVktTbZg0uB4ppiohNwNXA02AihX+TpDuAO/K2F0takz+NvyvpGY37SzpZ0i35k+5SYF5l3/Mlbaysr5T0BUn3S9oi6UJJTwEuBk7PGc62XHaupH+S9POc9Vws6fDKsf5c0mZJd0t6bbvHJ2kZ8ETgoxGxJ9/+OyK+k/c/RtJXcpu25uUVlftfL+md+XHvkPRlSY+V9Jn8yX6TpOMq5UPSmyXdKekBSe+V1PL1KOkXJV2bM511kn53mv+zRyLieuAc4HTgN/PxhnKW8dP8/H5e0tK877jcttfk5/QBSW+vtOVUSavzY7pX0vub7jdH0ruAXwYuzM/FhZIukvS+psd1paQ/m85j6bmI8K3NDfgZ8MK8vBK4FfiHvB6kT9+lwOHAycB9wCpgGHhNvv9c4DDgLuDPgBHgZcBe4J35WM8HNublYeCHwAeABaSA8ry873zgO01t/ABwZW7HIlJW8I9535nAvaTgtgD4j9zu41s8VpEC3leAlwBHNe1/LPA7wPxcz38CX6zsvx5YDzwZOAK4Dfhf4IXAHODTwL9WygfwrdzuJ+Syr29+nLndG4A/zMc5GXgAOKnN/+z6xnGatt8AvCcvvwW4EViR/z8fAT6b9x2X2/bR/H99JrAbeEre/z3gVXl5IXBa0/3mtGoHcCpwNzCU15cBO5uf53699bwB/XwjvdF3ANtIb/R/AQ7P+wJ4QaXsh8lBpLJtHfCrwK/kF4kq+75L60BxOnB/4wXXdLyJN1BeF/AI8OTKttOB/8vLnwDeXdn3C7QJFHn/CuBC4KfAeH5zndCm7LOArZX164G3V9bfB1xdWf8tYE1lPYAzK+t/Anyj+XECvwd8u6nujwDvaNOuSW/QyvbPkbIlgLXAGZV9x5AC95zKG35FZf8PgPPy8g3A3wHLmo7fuF/LQFGp90V5+QLgql6/xqd786lH2UsiYklEHBsRfxIRuyr7NlSWjwXemk87tuVTg5XA4/NtU+RXSHZXm/pWAndFxOg02nYk6RP+5kqdX8vbyfVW29iuTgAiYmNEXBART86P5xFSJoCk+ZI+IukuSQ+T3jBLmnr5760s72qx3twR3Ny2x7do1rHAqqbn9Q9I/SgHYjnwYOWYV1SOtxYYA46qlL+nsryz0vbXkQLu7fl06sUH0IZPAa/My68E/u2AHkEPOVB0pvrG3wC8KweVxm1+RHwW2Awsb+pMe0KbY24AntCmg7T5q74PkN6AT63UeUSkzldyvSunUef+FUVsAC4i98mQOgVPBFZFxGJSlgQpqzlYzW27u0WZDcB/NT2vCyPijdOtRNJK4DnAtyvHPKvpmPMi9UNNKSLuiIhXkDp73wNcJmlBq6Ittv07cK6kZwJPAb443cfQaw4U3fNR4I8lrVKyQNJvSlpEOq8dBd4saUTSS0nnrK38gPQGf3c+xjxJz8377gVWSDoMICLGc70fkPQ4AEnLJf1GLv954HxJJ0maD7yjXeNzZ+XfSTo+d/YtA15LOpeH1C+xC9iWO/7aHusA/HmudyWp3+DSFmW+AvyCpFfl525E0i/lzt0p5SzoV4EvkZ7Xq/Kui4F3STo2lztS0rnTabCkV0o6Mj/32/Lm8RZF7wWeVN0QERuBm0iZxOVN2Wlfc6DokohYDfwR6Rx/K6lj7/y8bw/w0rz+IOm8+wttjjNGOp8/Hvg5sDGXB/gmqUP1HkkP5G1/meu6MZ8SXEf65CcirgY+mO+3Pv9tZw/pPPs60iXRn5A68c7P+z9I6tx7gBQ8vjbFsabrS8DNwBrgq8DHmwtExHbg10mXOO8mnRK8h9QJ2c6FkraT3qwfBC4n9Yc03tD/TOoA/noudyOpE3o6zgRulbQjH+e8Nm/4fwZelq8Qfaiy/VPA0xmg0w7InWtmdZMUpI7S9b1uS50k/QrpFOTYGKA3nzMKs5pIGiGdYn1skIIEdCFQKA0O+pak2yTdKuktLcpI0ockrZf0I0nP7rRes0GS+1S2kS7FfrCnjTkI3Rh6PAq8NSJuyR13N0u6NiJuq5Q5Czgh31aRxhxM95zQZqGIOKSGU0fEWtLgsYHUcUYREZsj4pa8vJ10TXp5U7FzgU9HciPp+vsxndZtZvXo6peZlMbynwx8v2nXciYPrtmYt21ucYw3AG8AGGb4OfNZ3M0mmlnFdrY+EBFHlsp1LVAoff36cuBPI+Kgv3EYEZcAlwAs1tJYpTO61EIza3ZdXDblaN2Grlz1yL25lwOfiYhW4wM2MXkU3oq8zcwGQDeueog0UGZtRLy/TbErgVfnqx+nAQ9FxH6nHWbWn7px6vFc4FXAjyWtydv+mvy9goi4mDR09mzS6MCdpK8Mm9mA6DhQRPphkykvdeXBJW/qtC4z6w2PzDSzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMyvq1q9wf0LSfZJ+0mb/8yU9JGlNvv1tN+o1s3p0a16PTwIXAp+eosy3I+LFXarPzGrUlYwiIm4AHuzGscys/9TZR3G6pB9KulrSU2us18w61NW5R6dwC3BsROyQdDbwRdLM5vupzj06j/k1Nc/MplJLRhERD0fEjrx8FTAiaVmbspdExCkRccoIc+tonpkV1BIoJB2dpx5E0qm53i111G1mnevKqYekzwLPB5ZJ2gi8AxiBiSkFXwa8UdIosAs4L88eZmYDoCuBIiJeUdh/IenyqZkNII/MNLMiBwozK3KgMJtpGkq3ATbYrTezWtQ14Mrs0BXjvW5Bx5xRmFmRMwo7dLXrNziQDKBxjOnc50DK9hlnFGZW5EBhZkU+9bBDQ4vTDA0JgBiPpvWhSeutyzR9A2Gqy5+NU40BPOVocEZhZkXOKOzQlbMADTdWU7agOe2zj4n1sTGgkllUsoVW2wadMwozK3JGYYeGyqe7hocn7Wqsa97cSX+Ze1ilUM4oRlMmwZ49AIw/sjMdfs/eyhFzmQG+HNrMGYWZFTmjsENDq6sejUzi8Hnp75IjABg9ZgkAu46aN1F2dH66//zNKZOYu/nhdJ/cV8Ho6ETZGGssDH4m0eCMwsyKnFHYoWHSVYl8tSOvK/dFjB69BIB7Vi0AYPuzdk/cZ3huShPm35R+Gf7o3anMnK0PpWNWqpoYazHGrOGMwsyK6pp7VJI+JGm9pB9JenY36jWzetQ19+hZpAl/TgBWAR/Of81qt+/UIJ8bzElvg91Hpsuij57+CABfPu3iifs8GqnMy3e/EYCl69LpypyRkVSgesl1bBadc2R1zT16LvDpSG4Elkg6pht1m9nMq6szczmwobK+MW/b3FzQUwraTGtkEhNDtfMUM6Pz0vpjF6eM4qkj+15/d+zdkYruzV2gjWHa47mTVJOHeKdtHnA1YzyloFn/qSuj2ASsrKyvyNvMem93GkS18K6USWz478cB8PRHf3+iyOhY+kxddHvqkxjZnoZuN7IRxgY/a5hKXRnFlcCr89WP04CHImK/0w4z6091zT16FXA2sB7YCfxhN+o1OyhNw7nHd+0CYPi+NCx7xTfTFYwd6xdNlIncBbFg86MAjNyTytL4Mli1H2Li+LPn6kddc48G8KZu1GVm9fMQbjt0NbKA/IWu2JKu8A8/nLKFJT/b96UwhnKW0BgjkTOJ8Z072x93FlztaOi7qx5m1n+cUdihocXXzCd+sk7pb+xK/Q/ak/ooYscj7Y/T+Gm8xk/itRqN2VznAGcYzijMrMiBwsyKfOphh7zm04YYzZc8p5iro+38HrOUMwozK3JGYdbQ3Nk4Refjfr9eNZ2ZwgaYMwozK3JGYbNbq0/6mfiEbzmEe/aYfY/IzLrOGYXNLu0+zevsJ5huXdW29nk/hjMKMytyRmGzS59/Mk8yQG11RmFmRQ4UNtg0NLuuMvTp4+m/FplZ33GgMLOibk0peKakdXnKwLe12H++pPslrcm313ejXjNifP/BTn2Yuk9b8+PpEx1f9ZA0DFwEvIg0sc9Nkq6MiNuail4aERd0Wp+Z1a8bofdUYH1E3BkRe4DPkaYQNKtfn34iD7puBIp20wU2+508k/llkla22A+kKQUlrZa0ei+7u9A8M+tUXSdzXwaOi4hnANcCn2pX0FMKmvWfbgSK4nSBEbElIhrpwceA53ShXjOrSTcCxU3ACZKeKOkw4DzSFIITJB1TWT0HWNuFes2sJh1f9YiIUUkXANcAw8AnIuJWSX8PrI6IK4E3SzoHGAUeBM7vtF4zq48i+vfHQRdraazSGb1uhvVaY1yEr2Z03XVx2c0RcUqp3ACPTDGzujhQmFmRf4/C+p9POXrOGYWZFTlQmFmRA4WZFTlQmFmRA4WZFTlQmFmRA4WZFTlQmFmRA4WZFTlQmFmRA4WZFTlQmFmRA4WZFTlQmFmRA4WZFdU1peBcSZfm/d+XdFw36jWzenQcKCpTCp4FnAS8QtJJTcVeB2yNiOOBDwDv6bReM6tPXVMKnsu+SX8uA86QpC7UbWY1qGtKwYkyETEKPAQ8ttXBPKWgWf/pu85MTylo1n9qmVKwWkbSHOAIYEsX6jazGtQypWBef01efhnwzejnmYfMbJK6phT8OPBvktaTphQ8r9N6zaw+XZnXIyKuAq5q2va3leVHgZd3oy4zq1/fdWaaWf9xoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzIgcKMytyoDCzoo4ChaSlkq6VdEf++5g25cYkrcm35h/eNbM+12lG8TbgGxFxAvCNvN7Kroh4Vr6d02GdZlazTgNFdarATwEv6fB4ZtaHOg0UR0XE5rx8D3BUm3Lz8jSBN0p6yVQH9JSCZv2n+HP9kq4Djm6x6+3VlYgISe0m9Tk2IjZJehLwTUk/joiftioYEZcAlwAs1lJPEmTWB4qBIiJe2G6fpHslHRMRmyUdA9zX5hib8t87JV0PnAy0DBRm1n86PfWoThX4GuBLzQUkPUbS3Ly8DHgucFuH9ZpZjToNFO8GXiTpDuCFeR1Jp0j6WC7zFGC1pB8C3wLeHREOFGYDpKMpBSNiC3BGi+2rgdfn5e8CT++kHjPrLY/MnA00lG5mM8SvLjMrcqA4WL36FHf2YD3gV5yZFXXUmXlIi/F662vOIqrrdbfFDjnOKMysyBnFwWp8ojd/mrfqP+jkE9/9EdYH/Co0syJnFAerkSU0feJreLhFWeW7xOT7NmuRPWhomvc1m0HOKMysyIHCzIp86tEljVOExt/J0umI8lnJxGlEQ6vTiXwaEmNj3Wqi2UFzRmFmRc4oOjXRqTmS/6bYq8NG9pUZmfw0DzVnC6OjaX3Pnn2HdSZhfcQZhZkVOaPoVPPl0Zw9aOGCiW2xZDEAo8vStvE56T4j23YBMHTv1lTw4e37DpSzixgf7X6bzQ6QMwozK3JG0W2NAVdHLJrY9NAzlgJw/8m5/yJ3Pyz+v7kALLslXSkZ2rt34j6xN2USEwOu3GdhPdTplIIvl3SrpHFJp0xR7kxJ6yStl9RuNjEz61OdZhQ/AV4KfKRdAUnDwEXAi4CNwE2Srpx1P7A7cfUjZQBjRxw+sWvL01I8fubz7gBgjlLZ/xk7EYDHrDsMgKGh/b86vt+YC7Me6PTHddcCSK0GGU04FVgfEXfmsp8jTUU4uwKF2SxWRx/FcmBDZX0jsKpdYUlvAN4AMI/5M9uybmqMn8hXPWK4mh2kP9t2zwPg/h0LAZizM+9uMZrTmYT1k46mFIyI/Sb86ZSnFDTrPx1NKThNm4CVlfUVeZuZDYg6Tj1uAk6Q9ERSgDgP+P0a6q3VxJfBcn/N0KP7LnUuvT0lRnfvXjnpPovuTtvnbM+ztocTKOtPnV4e/W1JG4HTga9KuiZvf7ykqwAiYhS4ALgGWAt8PiJu7azZZlYnRR9/ii3W0lil/WYs7C+5E3PosHSJUwvSZdFJQ7gXpE7ZsSWHT7rr0J40iGpo645UbsvWiX3jj6SezomBVv5lK5sB18VlN0dE2zFQDR7CbWZFHsJ9sBqXQ5subcajub9hz74+CnY8AsCcLfnpHpucHUSLr5mb9RNnFGZW5IyiWxp9CKN56HV1X/6CV2Nbo99hv2ykOsiqcTz3TVgfcEZhZkXOKA5Um5m7xveONhXbly3s9+MzjS98eRYwGxB+pZpZkTOKA9U0Q1i7GbxiOr8zM1X/g7MN6yN+NZpZkQOFmRX51ONgzfRlS18WtT7ijMLMihwozKzIgcLMihwozKzIgcLMihwozKzIgcLMiuqaUvBnkn4saY2k1Z3UaWb1m/EpBSt+LSIe6LA+M+uBOqYUNLMBV1cfRQBfl3RznjLQzAZIXVMKPi8iNkl6HHCtpNsj4oY29Q3m3KNms1gdUwoSEZvy3/skXUGa4bxloPDco2b9Z8ZPPSQtkLSosQz8OqkT1MwGxIxPKQgcBXxH0g+BHwBfjYivdVKvmdWr06seVwBXtNh+N3B2Xr4TeGYn9ZhZb3lkppkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVOVCYWZEDhZkVdfrjuu+VdLukH0m6QtKSNuXOlLRO0npJb+ukTjOrX6cZxbXA0yLiGcD/An/VXEDSMHARcBZwEvAKSSd1WK+Z1aijQBERX4+I0bx6I7CiRbFTgfURcWdE7AE+B5zbSb1mVq9OZzOvei1waYvty4ENlfWNwKp2B6lOKQjsvi4um42TBS0DZuPM7rP1ccHsfWwnTqdQV+YelfR2YBT4zIG0sJXqlIKSVkfEKZ0es9/4cQ2e2frYJK2eTrmO5x6VdD7wYuCMiGg1V+gmYGVlfUXeZmYDotOrHmcCfwGcExE72xS7CThB0hMlHQacB1zZSb1mVq9Or3pcCCwCrpW0RtLFMHnu0dzZeQFwDbAW+HxE3DrN41/SYfv6lR/X4Jmtj21aj0utzxbMzPbxyEwzK3KgMLOivg4U0x0iPogkvVzSrZLGJQ38ZbfZOkxf0ick3SdpVo3nkbRS0rck3ZZfh2+ZqnxfBwqmMUR8gP0EeClwQ68b0qlZPkz/k8CZvW7EDBgF3hoRJwGnAW+a6n/W14FimkPEB1JErI2Idb1uR5fM2mH6EXED8GCv29FtEbE5Im7Jy9tJVySXtyvf14GiyWuBq3vdCGup1TD9ti866y+SjgNOBr7frkw3v+txUOoeIl6n6Tw2s16StBC4HPjTiHi4XbmeB4ouDBHvW6XHNot4mP4AkjRCChKfiYgvTFW2r089pjlE3HrPw/QHjCQBHwfWRsT7S+X7OlDQZoj4bCDptyVtBE4Hvirpml636WB1OEy/r0n6LPA94ERJGyW9rtdt6pLnAq8CXpDfW2sknd2usIdwm1lRv2cUZtYHHCjMrMiBwsyKHCjMrMiBwsyKHCjMrMiBwsyK/h+MhDv3Nz6aFAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}