phys512/lcg.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "879902dd",
   "metadata": {},
   "source": [
    "# Linear Congruential Generators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bf6fefd",
   "metadata": {},
   "source": [
    "These widely-used generators are of the form \n",
    "\n",
    "$$X_{n+1} = (aX_n + c)\\ \\mathrm{mod}\\ m.$$\n",
    "\n",
    "Notes:\n",
    "\n",
    "- The starting value of $X$ is the *seed*. If the sequence ever produces the seed again at a later iteration, it will start over with the same sequence of numbers. How long this takes is the *period* of the generator.\n",
    "- The modulus $m$ sets the maximum period of the random number sequence (although the period can be much smaller for bad choices of $a$ and $c$).\n",
    "- To return a floating point number between 0 and 1, we can divide by $m$.\n",
    "- A list of different choices of parameters is given on the [Wikipedia page for linear congruential generators](https://en.wikipedia.org/wiki/Linear_congruential_generator) -- we explore some of these below and compare with the build in numpy generator. (Note that you should always use a well-tested generator like the numpy one rather than coding your own, the tests below should show you why!)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3598ae23",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "73fab248",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the implementation of an LCG from Wikipedia\n",
    "# https://en.wikipedia.org/wiki/Linear_congruential_generator\n",
    "\n",
    "from collections.abc import Generator\n",
    "\n",
    "def lcg(modulus: int, a: int, c: int, seed: int) -> Generator[int, None, None]:\n",
    "    \"\"\"Linear congruential generator.\"\"\"\n",
    "    while True:\n",
    "        seed = (a * seed + c) % modulus\n",
    "        yield seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "d2b53d68",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9d5hkZ3nmj38q59A5554saaKyCEISIpu1ARFMWBZ8gW28WJb5GRu8trFXa7D5AgYTjL1ggsGsTZYlhOJEzUynmek00znHyqkr/v5o3ofqnp4kzYxG0vu5Ll2arq4659Sp0O997ue5H0OhUCig0Wg0Go1Go9FoNJrLhvH5PgCNRqPRaDQajUajebGhhZZGo9FoNBqNRqPRXGa00NJoNBqNRqPRaDSay4wWWhqNRqPRaDQajUZzmdFCS6PRaDQajUaj0WguM1poaTQajUaj0Wg0Gs1lRgstjUaj0Wg0Go1Go7nMaKGl0Wg0Go1Go9FoNJcZLbQ0Go1Go9FoNBqN5jKjhZZGo9FoNBqNRqPRXGa00NJoNBrNFcNgMFzwv7/4i7/g0UcfxWAw8Jd/+ZdnbWN0dBSn08lb3vKWNbf/8pe/5M4776S8vBy/389NN93Et771rbMe/+Uvf5m3vvWtNDY2YjAYeN/73nelnq5Go9FoNIL5+T4AjUaj0bx42Uj4KP7iL/6C4eFhbr75Zu655x7e+c538uCDD/KOd7yDzZs3y/1+93d/F4vFwhe+8AW57Sc/+QlvfvObufXWW/mLv/gLDAYD//7v/8573vMelpaW+MM//EO579/+7d8SjUa56aabmJ2dvTJPVKPRaDSadRgKhULh+T4IjUaj0by0+PrXv84HP/hBPvKRj4iAWlhYYOvWrezatYvHH38cgO9973u84x3v4Atf+AIf+chH5PGvfvWr6e3tZWRkBJvNBkA2m2Xr1q24XC56enrkvuPj4+Jmud1u3vKWt/CNb3zj6j1ZjUaj0bwk0aWDGo1Go7mq9Pb28gd/8Afs3r2bz3zmM3J7ZWUlf/u3f8sTTzzBN7/5TUKhEH/4h3/IjTfeyO/93u+t2UYkEqGkpEREFoDZbKa8vByHw7Hmvk1NTRgMhiv7pDQajUajWYcWWhqNRqO5aiQSCd72trdhMpn43ve+t0YoAXzgAx/g9ttv54EHHuB3f/d3WVxc5Ktf/SpG49o/V6985Svp7e3lk5/8JENDQwwPD/OpT32K48eP87GPfexqPiWNRqPRaDZE92hpNBqN5qrxkY98hL6+Pr75zW+u6cNSGAwGvvrVr7J7927+7d/+jY9+9KPs3r37rPt98pOfZHR0lL/5m7/hr//6rwFwOp38x3/8B7/xG79xxZ+HRqPRaDQXQjtaGo1Go7kqfPe73+Vf/uVfePe738173vOec97P6/VitVqB1V6sjbDZbGzevJm3vOUt/Nu//Rvf/va32bdvH7/927/NkSNHrsjxazQajUZzKegwDI1Go9Fccc6cOcPevXupqamho6MDt9t9zvv+xm/8Bk888QSlpaU4nU56enqwWCxr7vOhD32II0eO0NnZKWWFmUyGHTt2UFJSwjPPPLPhtnUYhkaj0WiuFtrR0mg0Gs0VZWVlhfvuu490Os33vve984qs//zP/+QnP/kJn/rUp/jHf/xH+vv71wRmAKTTaf75n/+Z17/+9Wt6tywWC6997Ws5fvw46XT6ij0fjUaj0WguBi20NBqNRnNFeeCBB+jq6uLTn/70hv1Wimg0yh/8wR+wZ88efv/3f5/Xve51/NZv/RZ//dd/zejoqNxveXmZbDZLLpc7axuZTIZ8Pr/h7zQajUajuZpooaXRaDSaK8YPf/hDvvjFL/KmN72JP/iDPzjvfT/xiU8wOzvLV7/6VUwmEwCf//znMZlM/P7v/77cr7KyEr/fzw9/+MM1zlUsFuOnP/0pW7duPSviXaPRaDSaq41OHdRoNBrNFWF2dpb/8T/+ByaTibvuuotvf/vbG96vra0Nq9XKl770JX7v936Pffv2ye/q6ur4q7/6K+6//37+4z/+g9/6rd/CZDLxwAMP8IlPfIJbbrmF97znPeRyOf75n/+Zqamps/bz05/+VAYYZzIZTpw4IUmFb3rTm7jhhhuu0BnQaDQazUsZHYah0Wg0mivCk08+yZ133nnB+/32b/82fX19zM3N0d/fj9frXfP7XC7Hvn37WFpaor+/X3q8vvvd7/L5z3+e06dPs7Kywg033MAf//Ef81u/9VtrHv++972Pb37zmxvu+//+3//L+973vmf3BDUajUajOQ9aaGk0Go1Go9FoNBrNZUb3aGk0Go1Go9FoNBrNZUYLLY1Go9FoNBqNRqO5zGihpdFoNBqNRqPRaDSXGS20NBqNRqPRaDQajeYyo4WWRqPRaDQajUaj0VxmXhBztPL5PDMzM3g8HgwGw/N9OBrNS5pCoUA0GqW2thaj8YVxrUZ/h2g01w4vxO8QjUajeTa8IITWzMwMDQ0Nz/dhaDSaIiYnJ6mvr3++D+Oi0N8hGs21xwvpO0Sj0WieDS8IoeXxeIDVL+X1gyw1Gs3VJRKJ0NDQIJ/LFwL6O0SjuXZ4IX6HaDQazbPhBSG0VKmP1+vViySN5hrhhVSCp79DNJprjxfSd4hGo9E8G3RxtEaj0Wg0Go1Go9FcZrTQ0mg0Go1Go9FoNJrLjBZaGo1Go9FoNBqNRnOZ0UJLo9FoNBqNRqPRaC4zWmhpNJqrytNPP80b3/hGamtrMRgM/OhHP7rgY5588kn27NmDzWajvb2db3zjG1f8ODUajUaj0WieC1poaTSaq0o8Hmfnzp186Utfuqj7j46O8vrXv54777yT7u5uPvrRj/KBD3yARx555AofqUaj0Wg0Gs2z5wUR767RaF48vPa1r+W1r33tRd//K1/5Ci0tLfz93/89ANu2bePAgQP8f//f/8e99957pQ5To9FoNM8DhUKBRCJBPp/HarVisVgwGrUvoHlhooWW5gXFdChJMJ6mxGWlzu94vg9HcxU4fPgwd99995rb7r33Xj760Y+e8zErKyusrKzIz5FI5EodnkZzWcjn8zzyyCP813/9F0tLS8zMzFBTU0MsFqOhoYHJyUncbjcjIyPs2bOHSCRCJBLB7XYzOzuL2+3G7XYzPz+Py+WisbERo9FIKBTC4/HQ2dnJ3r17MZlMZLNZurq62LNnDxUVFYyMjOB2u+ns7GTPnj2Ew2Gmp6dxuVwYDAa8Xi+zs7M4HA6mpqa4/fbbicfjtLa2Mjw8TC6XY25ujtraWt71rnfxmte8Ri+MNc+KQqFAOp0mnU6Ty+XIZDIYjUbMZjNWqxWz2azfW5oXFFpoaV4wTIeS3P33T5HM5HBYTPzyj16hxdZLgLm5OaqqqtbcVlVVRSQSIZlM4nCc/R548MEH+cu//MurdYia81AoFOjt7eUTn/gELpcLr9fL1NQUhUIBr9eL0WiURb7D4WD79u1EIhEMBgP5fJ5oNEpDQwPRaJRcLofBYKBQKGA0GmltbaWrqwuAxsZGKioqCAQCLC0tEY1GaWxsBGBsbIxEIoHT6cTn8xEKhZiamuK2224jGo0CMDU1hdfrxePxrBmkm8/n5Xg8Hg9zc3M4nU459nA4jMfjIRaL4fF48Pl8dHV1sWvXLgYGBnC5XESjUXw+Hx6PB4BwOEw8Hqe2tlb2MzQ0xPLyMj6fj8nJSQwGA52dnXi9Xg4dOoTD4SAYDOL1ehkcHMThcFAoFEgmk5hMJvL5PIVCAZPJhNFopLu7W5yAQCCA2+1meHiYsrIylpaWcDqdjI6O4na7sdvtzMzM4PV6+bd/+zesViu5XI50Oo3NZiOdTmM2m+U5fuc736GhoYGnnnoKi8VCMBjEZDIxOTlJoVDAYrFwzz33XK23mOZFghJZuVwOQARWoVAgk8mQTqcxGo1YLBYsFosWXZoXBFpoXeNoB+fXBONpkpkcH3lVO//w+BDBePolf040G/Pxj3+c+++/X36ORCI0NDRc1n0UCgX6+/v56le/ytLSEqFQCIDFxUXKy8tpbW2lrKyMzs5O8vk8sVgMg8GwxqW
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d3icZ5nv/5muGY2mSCNp1Ht1kSX3XuLEMUlIgSQEEsJCloVsY9mFJbuUPcs5nLD8WGBZ6gJpwJICMXHsuNtxL7Ilq1i9ayTNjKSRpvf5/eEzz9qO7diJEyfh/VyXLlszb3neMqPn+973/b1liUQigYSEhISEhISEhISEhMQNQ36zByAhISEhISEhISEhIfFBQxJaEhISEhISEhISEhISNxhJaElISEhISEhISEhISNxgJKElISEhISEhISEhISFxg5GEloSEhISEhISEhISExA1GEloSEhISEhISEhISEhI3GEloSUhISEhISEhISEhI3GAkoSUhISEhISEhISEhIXGDkYSWhISEhISEhISEhITEDUYSWhISEhIS7woymYx/+Zd/eUe2PTg4iEwm4+mnn35Htn+9PP3008hkMgYHB2/2UCQkJCQkbhKS0JKQkJB4D5CcmKekpGCz2d7w/rp165g7d+5NGJmEhISEhITEW0ESWhISEhLvIUKhEE8++eTNHsY7QiAQ4Ktf/erNHsa7wiOPPEIgEKCoqOhmD0VCQkJC4iYhCS0JCQmJ9xALFizgv/7rvxgbG7vZQ7khxONxgsEgACkpKSiVyps8ojfi8/lu+DYVCgUpKSnIZLIbvm0JCQkJifcHktCSkJCQeA/xT//0T8RisTeNal2tJunSWqh/+Zd/QSaT0d3dzcMPP4zRaCQzM5Ovfe1rJBIJRkZGuPvuuzEYDFitVr773e++YZuhUIhvfOMblJeXo9FoKCgo4Mtf/jKhUOgN+/6rv/orfvOb3zBnzhw0Gg07duy47LgAbDYbn/nMZ8jNzUWj0VBSUsLnP/95wuEwANPT0/zDP/wD8+bNQ6/XYzAY2Lx5M2fPnr2Gs/lGkimar7/+Oo8//jhZWVnk5+eL91977TVWr15NamoqaWlp3HHHHbS3t79hO52dnTzwwANkZmai1Wqpqqrin//5n9+wnwtrtIqLi7nzzjvZtWsXCxYsICUlhdraWv7whz+8YfszMzN84QtfoKCgAI1GQ3l5Od/+9reJx+MXLfe73/2OhQsXkpaWhsFgYN68efzgBz94S+dGQkJCQuLG8t57tCghISHxJ0xJSQmf/OQn+a//+i++8pWvkJube8O2/eCDD1JTU8OTTz7Jtm3b+N//+3+Tnp7Oz372MzZs2MC3v/1tfvOb3/AP//APLF68mDVr1gDno1If/vCHOXz4MJ/97GepqamhtbWV733ve3R3d7Nly5aL9rNv3z5eeOEF/uqv/gqLxUJxcfFlxzM2NsaSJUuYmZnhs5/9LNXV1dhsNl566SX8fj9qtZr+/n62bNnC/fffT0lJCXa7nZ/97GesXbuWc+fOveXz8/jjj5OZmcnXv/51EdF67rnnePTRR9m0aRPf/va38fv9/OQnP2HVqlU0NTWJ42hpaWH16tWoVCo++9nPUlxcTF9fH1u3buX//J//c9X99vT08OCDD/K5z32ORx99lKeeeor777+fHTt2cOuttwLg9/tZu3YtNpuNv/iLv6CwsJCjR4/yxBNPMD4+zve//30Adu/ezUMPPcQtt9zCt7/9bQA6Ojo4cuQIf/u3f/uWzouEhISExA0kISEhISFx03nqqacSQOLUqVOJvr6+hFKpTPzN3/yNeH/t2rWJOXPmiN8HBgYSQOKpp556w7aAxDe+8Q3x+ze+8Y0EkPjsZz8rXotGo4n8/PyETCZLPPnkk+J1l8uV0Gq1iUcffVS89txzzyXkcnni0KFDF+3npz/9aQJIHDly5KJ9y+XyRHt7+5uO65Of/GRCLpcnTp069YZl4/F4IpFIJILBYCIWi1303sDAQEKj0ST+9V//9ZrOx4Ukz/OqVasS0WhUvO7xeBImkynx53/+5xctPzExkTAajRe9vmbNmkRaWlpiaGjosmO+cD8DAwPitaKiogSQ+P3vfy9em52dTeTk5CTq6+vFa9/85jcTqampie7u7ou2/5WvfCWhUCgSw8PDiUQikfjbv/3bhMFguOg4JCQkJCTeO0ipgxISEhLvMUpLS3nkkUf4+c9/zvj4+A3b7mOPPSb+r1AoWLRoEYlEgs985jPidZPJRFVVFf39/eK1F198kZqaGqqrq5mcnBQ/GzZsAGD//v0X7Wft2rXU1tZedSzxeJwtW7Zw1113sWjRoje8n6xt0mg0yOXn/1TFYjGmpqbQ6/VUVVVx5syZ6zwD/8Of//mfo1AoxO+7d+9mZmaGhx566KJjVCgULF26VByj0+nk4MGDfPrTn6awsPCyY74aubm53HvvveJ3g8HAJz/5SZqampiYmADOn+/Vq1djNpsvGsvGjRuJxWIcPHgQOH+tfD4fu3fvfsvnQUJCQkLinUNKHZSQkJB4D/LVr36V5557jieffPKG1dxcKgyMRiMpKSlYLJY3vD41NSV+7+npoaOjg8zMzMtu1+FwXPR7SUnJm47F6XTidrvf1LI+Ho/zgx/8gB//+McMDAwQi8XEexkZGW+6nytx6Rh7enoAhHi8FIPBACAE6Fu12i8vL3+DIKusrATO191ZrVZ6enpoaWl50/P9+OOP88ILL7B582by8vK47bbbeOCBB7j99tvf0tgkJCQkJG4sktCSkJCQeA9SWlrKww8/zM9//nO+8pWvvOH9K0VPLhQil3JhBOdqrwEkEgnx/3g8zrx58/j3f//3yy5bUFBw0e9arfaKY7hevvWtb/G1r32NT3/603zzm98kPT0duVzOF77whTcYQ1wPl44xua3nnnsOq9X6huXfTbfEeDzOrbfeype//OXLvp8UZllZWTQ3N7Nz505ee+01XnvtNZ566ik++clP8swzz7xr45WQkJCQuDyS0JKQkJB4j/LVr36VX//618Lo4ELMZjNw3p3uQoaGhm74OMrKyjh79iy33HLLDbMrz8zMxGAw0NbWdtXlXnrpJdavX88vf/nLi16fmZl5QyTu7VBWVgacFy8bN2684nKlpaUAbzruK9Hb20sikbjoPHZ3dwMIs42ysjK8Xu9Vx5FErVZz1113cddddxGPx3n88cf52c9+xte+9jXKy8vf0hglJCQkJG4MUo2WhISExHuUsrIyHn74YX72s5+J+p0kBoMBi8Ui6nWS/PjHP77h43jggQew2Wz813/91xveCwQCb6kPlVwu55577mHr1q00Nja+4f1kRE2hUFwUXYPzNUw2m+2693k1Nm3ahMFg4Fvf+haRSOQN7zudTuC8QFyzZg2/+tWvGB4evuyYr8bY2Bgvv/yy+N3tdvPss8+yYMECEUl74IEHOHbsGDt37nzD+jMzM0SjUYCL0jvh/DmdP38+wBts9yUkJCQk3n2kiJaEhITEe5h//ud/5rnnnqOrq4s5c+Zc9N5jjz3Gk08+yWOPPcaiRYs4ePCgiI7cSB555BFeeOEFPve5z7F//35WrlxJLBajs7OTF154gZ07d17W0OLN+Na3vsWuXbtYu3atsI0fHx/nxRdf5PDhw5hMJu68807+9V//lT/7sz9jxYoVtLa28pvf/EZElm4UBoOBn/zkJzzyyCM0NDTwsY99jMzMTIaHh9m2bRsrV67kP//zPwH4j//4D1atWkVDQwOf/exnKSkpYXBwkG3bttHc3HzV/VRWVvKZz3yGU6dOkZ2dza9+9SvsdjtPPfWUWOZLX/oSr7zyCnfeeSef+tSnWLhwIT6fj9bWVl566SUGBwexWCw89thjTE9Ps2HDBvLz8xkaGuKHP/whCxYsoKam5oaeHwkJCQmJ60cSWhISEhLvYcrLy3n44YcvW3Pz9a9/HafTyUsvvSRMEV577TWysrJu6Bjkcjlbtmzhe9/7Hs8++ywvv/wyOp2O0tJS/vZv/1bUDF0veXl5nDhxgq997Wv85je/we12k5eXx+b
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d3yb53nuf2FzA9ykKJIYBJdIyeJeWpZ34tZN06bJOWkdy2mWkzjJLzkn5zSzadJm1dlxbCduUjtJkzq2kziWZFmL4h6SuEAQxCDFBS6AxCDm+/tD57nzgkOWZGvZz/fz8cc2AAIP3hcgn+u97vu6JYIgCOBwOBwOh8PhcDgczhuG9EYvgMPhcDgcDofD4XDebHChxeFwOBwOh8PhcDhvMFxocTgcDofD4XA4HM4bDBdaHA6Hw+FwOBwOh/MGw4UWh8PhcDgcDofD4bzBcKHF4XA4HA6Hw+FwOG8wXGhxOBwOh8PhcDgczhsMF1ocDofD4XA4HA6H8wbDhRaHw+FwOBwOh8PhvMFwocXhcDgczhvM/v37sX///hu9DA6Hw+HcQLjQ4nA4nBvMD3/4Q0gkEtTX19/opdw0vPTSS/jiF794o5dxS/HrX/8a//N//k8YjUZIJBIu9DgcDucGw4UWh8Ph3GCeeeYZaLVadHV1wWKx3Ojl3BS89NJL+NKXvnSjl3FL8aMf/QgvvPAC8vPzkZqaeqOXw+FwOG95uNDicDicG4jNZkNbWxu+/e1vIzMzE88888yNXhLnFuUXv/gF3G43Xn31VWzbtu1GL4fD4XDe8nChxeFwODeQZ555BqmpqXjb296Gd77znZsKLbvdDolEgm9+85v4yU9+AoPBAJVKhdraWnR3d8c89vz583jwwQeh1+sRFxeHnJwcPPTQQ1hcXNzwvCdOnEBNTQ3i4uJgMBjw+OOP44tf/CIkEknM4yQSCR555BE8//zzqKiogEqlwo4dO/Dyyy9veM7+/n7ce++9SElJQVJSEg4ePIiOjo6Yx4RCIXzpS1+C0WhEXFwc0tPT0dLSgqNHjwIAHnzwQfzgBz+g12b/MLxeLz71qU8hPz8fKpUKJSUl+OY3vwlBEDZd9zPPPIOSkhLExcWhuroap06duqp1b3ZsAODpp5+GRCKB3W7fcJ8Yp9OJQ4cOITs7G3Fxcdi1axf+4z/+Y8PjfvWrX6G6uhrJyclISUlBZWUlvvOd71zyuQEgPz8fUin/s87hcDg3C/IbvQAOh8N5K/PMM8/gHe94B5RKJd797nfjRz/6Ebq7u1FbW7vhsc8++yxWV1fxgQ98ABKJBF//+tfxjne8A1arFQqFAgBw9OhRWK1WvO9970NOTg6Ghobwk5/8BENDQ+jo6CCh0N/fj3vuuQe5ubn40pe+hEgkgi9/+cvIzMzcdJ2tra147rnn8OEPfxjJycn47ne/i7/+67/GxMQE0tPTAQBDQ0PYs2cPUlJS8JnPfAYKhQKPP/449u/fj5MnT1IP2he/+EV87Wtfw8MPP4y6ujqsrKygp6cHfX19uPPOO/GBD3wA09PTOHr0KH7xi1/ErEMQBPzFX/wFjh8/jkOHDuG2227D4cOH8elPfxpTU1P493//95jHnzx5Er/+9a/xsY99DCqVCj/84Q9xzz33oKurCxUVFVe07teD3+/H/v37YbFY8Mgjj0Cn0+E3v/kNHnzwQbhcLnz84x+n8/fud78bBw8exL/9278BAEZGRnDmzBl6DIfD4XBuEQQOh8Ph3BB6enoEAMLRo0cFQRCEaDQqbN++Xfj4xz8e8zibzSYAENLT04WlpSW6/YUXXhAACL///e/pNp/Pt+F1fvnLXwoAhFOnTtFt999/v5CQkCBMTU3RbWNjY4JcLhfW/2kAICiVSsFisdBt586dEwAI3/ve9+i2Bx54QFAqlcL4+DjdNj09LSQnJwt79+6l23bt2iW87W1vu+Sx+chHPrJhHYIgCM8//7wAQPjKV74Sc/s73/lOQSKRxKwRgABA6OnpodscDocQFxcn/NVf/dUVr/sLX/jCpmv62c9+JgAQbDYb3bZv3z5h37599P+PPfaYAED4z//8T7otGAwKjY2NQlJSkrCysiIIgiB8/OMfF1JSUoRwOHypw/Oa7NixI+b1ORwOh3P94TUGHA6Hc4N45plnkJ2djQMHDgC4WOr2rne9C7/61a8QiUQ2PP5d73pXTMjBnj17AABWq5Vui4+Pp/9eW1vDwsICGhoaAAB9fX0AgEgkgldeeQUPPPBATC9PUVER7r333k3Xescdd8BgMND/79y5EykpKfTakUgER44cwQMPPAC9Xk+Py83NxXve8x60trZiZWUFAKDRaDA0NISxsbHLOUwxvPTSS5DJZPjYxz4Wc/unPvUpCIKAP/3pTzG3NzY2orq6mv6/oKAAf/mXf4nDhw8jEolc0bpfDy+99BJycnLw7ne/m25TKBT42Mc+Bo/Hg5MnTwK4eGy8Xi+VUXI4HA7n1oULLQ6Hw7kBRCIR/OpXv8KBAwdgs9lgsVhgsVhQX1+Pubk5HDt2bMPPFBQUxPw/E13Ly8t029LSEj7+8Y8jOzsb8fHxyMzMhE6nAwC43W4AF3uF/H4/ioqKNrzGZrdt9trs9dlrz8/Pw+fzoaSkZMPjysrKEI1GMTk5CQD48pe/DJfLheLiYlRWVuLTn/40zp8/v+nrrsfhcGDbtm1ITk7e8BrsfjFGo3HDcxQXF8Pn82F+fv6K1v16cDgcMBqNG3qo1q/7wx/+MIqLi3Hvvfdi+/bteOihhzbtheNwOBzOzQ8XWhwOh3MDePXVVzEzM4Nf/epXMBqN9M/f/u3fAsCmoRgymWzT5xJEIRB/+7d/iyeeeAIf/OAH8dxzz+HIkSO0UY9Go1e93st57ctl7969GB8fx09/+lNUVFTgySefRFVVFZ588smrXt/1YLMgDACbuo9XS1ZWFs6ePYsXX3yRetHuvfde/MM//MMb9hocDofDuT7wMAwOh8O5ATzzzDPIysqidD0xzz33HH73u9/hxz/+cUwp4GuxvLyMY8eO4Utf+hI+//nP0+3rS/SysrIQFxe36cyuq53jlZmZiYSEBIyOjm64z2QyQSqVIj8/n25LS0vD+973Przvfe+Dx+PB3r178cUvfhEPP/wwgK1FTWFhIV555RWsrq7GuFomk4nuF7NZeaLZbEZCQgIFf1zuupmD6HK5oNFo6HHrXbSt1n3+/HlEo9EYV2uzdSuVStx///24//77EY1G8eEPfxiPP/44Pve5z23pOHI4HA7n5oM7WhwOh3Od8fv9eO655/D2t78d73znOzf888gjj2B1dRUvvvjiFT0vc53Wu0yPPfbYhsfdcccdeP755zE9PU23WyyWDT1OV/Lad911F1544YWYmPO5uTk8++yzaGlpQUpKCgBsiJpPSkpCUVERAoEA3ZaYmAjgoqgRc9999yESieD73/9+zO3//u//DolEsqHHrL29nXrTAGBychIvvPAC7rrrLshksitaN+tRE8fDe73eTSPa13PfffdhdnYWv/71r+m2cDiM733ve0hKSsK+ffs2PTZSqRQ7d+4EgJjjw+FwOJybH+5ocTgcznXmxRdfxOrqKv7iL/5i0/sbGhpoePG73vWuy37elJQU7N27F1//+tcRCoWQl5eHI0eOwGazbXjsF7/4RRw5cgTNzc340Ic+ROKloqICZ8+evar39ZWvfAVHjx5FS0sLPvzhD0Mul+Pxxx9HIBDA17/+dXpceXk59u/fj+rqaqSlpaGnpwe//e1v8cgjj9BjWIDFxz72Mdx9992QyWT4u7/7O9x///04cOAA/u///b+w2+3YtWsXjhw5ghdeeAGPPvpoTGAHAFRUVODuu++OiXcHgC996UtXvO677roLBQUFOHToED796U9DJpPhpz/9KTIzMzExMXHJY/OP//iPePzxx/Hggw+it7cXWq0Wv/3tb3HmzBk89thj5M49/PDDWFpawu23347t27fD4XDge9/7Hm677Tbq59qKU6dOkQicn5+H1+vFV77yFQAXyzX37t17yZ/ncDgczhvMjQ095HA4nLc
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d3gc53mvfW+vWJTFoi2w6JUgCRYQhNgpiaRESbQV2XK37Pgcd8u2jp04zueSOMmJbVnOUdwSx5FtyVazrEZK7L2CJED03usCi7KL7WW+P3j2PSRFUqREiZI893XhkjCY8s7MzvL9zfM8v0chSZKEjIyMjIyMjIyMjIyMzA1DebMHICMjIyMjIyMjIyMj815DFloyMjIyMjIyMjIyMjI3GFloycjIyMjIyMjIyMjI3GBkoSUjIyMjIyMjIyMjI3ODkYWWjIyMjIyMjIyMjIzMDUYWWjIyMjIyMjIyMjIyMjcYWWjJyMjIyMjIyMjIyMjcYGShJSMjIyMjIyMjIyMjc4ORhZaMjIyMjIyMjIyMjMwNRhZaMjIyMjIyN5j169ezfv36mz0MGRkZGZmbiCy0ZGRkZG4yP//5z1EoFNTU1Nzsobxj2LFjB9/73vdu9jDeNbhcLn70ox+xdu1abDYbSUlJrFy5kqeeeupmD01GRkbmLxZZaMnIyMjcZJ544gny8vI4deoU3d3dN3s47wh27NjB97///Zs9jHcNx48f59vf/jYpKSn8/d//Pf/0T/+E0WjkQx/6EN/97ndv9vBkZGRk/iKRhZaMjIzMTaSvr49jx47xk5/8BJvNxhNPPHGzhyTzLmTBggV0dXXx/PPP8+CDD/LFL36RvXv3snHjRv71X/8Vr9d7s4coIyMj8xeHLLRkZGRkbiJPPPEEycnJbN26lfvuu++yQqu/vx+FQsGPf/xj/uM//oPCwkJ0Oh3V1dXU1dVdtG5jYyMPPPAABQUF6PV6MjIy+PSnP43L5XrNfg8cOMDy5cvR6/UUFhbyq1/9iu9973soFIqL1lMoFHzpS1/i+eefp7KyEp1Ox4IFC3j11Vdfs8/6+nruuOMOLBYLZrOZW2+9lRMnTly0Tjgc5vvf/z7FxcXo9XqsViurV69m9+7dADzwwAP87Gc/E8eO/8Txer089NBD5OTkoNPpKC0t5cc//jGSJF123E888QSlpaXo9XqWLVvGoUOH3tC4L3dtAB577DEUCgX9/f2v+duFOJ1O/vqv/5r09HT0ej2LFy/mt7/97WvWe/LJJ1m2bBkJCQlYLBYWLlzIv/3bv1113/n5+eTm5r7m/N/3vvcRDAbp7e296vYyMjIyMjce9c0egIyMjMxfMk888QT33nsvWq2WD3/4w/ziF7+grq6O6urq16z7hz/8AY/Hw2c/+1kUCgU//OEPuffee+nt7UWj0QCwe/duent7+dSnPkVGRgYtLS38x3/8By0tLZw4cUIIhfr6erZs2UJmZibf//73iUaj/MM//AM2m+2y4zxy5AjPPfccX/jCF0hISOD//J//w1/91V8xODiI1WoFoKWlhTVr1mCxWPjmN7+JRqPhV7/6FevXr+fgwYOiBu173/se//Iv/8JnPvMZVqxYgdvt5vTp05w9e5bbb7+dz372s4yOjrJ7925+//vfXzQOSZK455572L9/P3/9139NVVUVO3fu5Bvf+AYjIyM88sgjF61/8OBBnnrqKb7yla+g0+n4+c9/zpYtWzh16hSVlZXXNe43g9/vZ/369XR3d/OlL32J/Px8nnnmGR544AFmZ2d58MEHxf378Ic/zK233sq//uu/AtDW1sbRo0fFOtfD+Pg4AKmpqW/6HGRkZGRkrhNJRkZGRuamcPr0aQmQdu/eLUmSJMViMSk7O1t68MEHL1qvr69PAiSr1SpNT0+L5S+88IIESC+99JJY5vP5XnOcP/7xjxIgHTp0SCy7++67JaPRKI2MjIhlXV1dklqtli79pwGQtFqt1N3dLZadO3dOAqRHH31ULHvf+94nabVaqaenRywbHR2VEhISpLVr14plixcvlrZu3XrVa/PFL37xNeOQJEl6/vnnJUD6wQ9+cNHy++67T1IoFBeNEZAA6fTp02LZwMCApNfrpfe///3XPe7vfve7lx3Tf//3f0uA1NfXJ5atW7dOWrdunfj9pz/9qQRIjz/+uFgWCoWk2tpayWw2S263W5IkSXrwwQcli8UiRSKRq12ea8LlcklpaWnSmjVr3vS+ZGRkZGSuHzl1UEZGRuYm8cQTT5Cens6GDRuA86le999/P08++STRaPQ1699///0kJyeL39esWQNwUVqYwWAQ/x8IBJiammLlypUAnD17FoBoNMqePXt43/veR1ZWlli/qKiIO+6447Jjve222ygsLBS/L1q0CIvFIo4djUbZtWsX73vf+ygoKBDrZWZm8pGPfIQjR47gdrsBSEpKoqWlha6urmu5TBexY8cOVCoVX/nKVy5a/tBDDyFJEq+88spFy2tra1m2bJn43eFwsG3bNnbu3Ek0Gr2ucb8ZduzYQUZGBh/+8IfFMo1Gw1e+8hXm5+c5ePAgcP7aeL1ekUb5RonFYnz0ox9ldnaWRx999E3tS0ZGRkbmjSELLRkZGZmbQDQa5cknn2TDhg309fXR3d1Nd3c3NTU1TExMsHfv3tds43A4Lvo9LrpmZmbEsunpaR588EHS09MxGAzYbDby8/MBmJubA87XCvn9foqKil5zjMstu9yx48ePH3tychKfz0dpaelr1isvLycWizE0NATAP/zDPzA7O0tJSQkLFy7kG9/4Bo2NjZc97qUMDAyQlZVFQkLCa44R//uFFBcXv2YfJSUl+Hw+Jicnr2vcb4aBgQGKi4tRKi/+Z/fScX/hC1+gpKSEO+64g+zsbD796U9fthbu9fjyl7/Mq6++yq9//WsWL178pscvIyMjI3P9yEJLRkZG5iawb98+xsbGePLJJykuLhY/H/zgBwEua4qhUqkuuy/pAhOID37wg/znf/4nn/vc53juuefYtWuXmKjHYrE3PN5rOfa1snbtWnp6evjNb35DZWUlv/71r1m6dCm//vWv3/D43g4uZ4QBXDb6+EZJS0ujoaGBF198UdSi3XHHHXzyk5+85n18//vf5+c//zn/+3//bz7+8Y/fsLHJyMjIyFwfshmGjIyMzE3giSeeIC0tTbjrXchzzz3Hn//8Z375y19elAr4eszMzLB3716+//3v853vfEcsvzRFLy0tDb1ef9meXW+0j5fNZsNoNNLR0fGav7W3t6NUKsnJyRHLUlJS+NSnPsWnPvUp5ufnWbt2Ld/73vf4zGc+A1xZ1OTm5rJnzx48Hs9FUa329nbx9wu5XHpiZ2cnRqNRGH9c67jjEcTZ2VmSkpLEepdG0a407sbGRmKx2EVRrcuNW6vVcvfdd3P33XcTi8X4whe+wK9+9Sv+v//v/7tixDHOz372M773ve/x1a9+lb/5m7953XHJyMjIyLx1yBEtGRkZmbcZv9/Pc889x1133cV99933mp8vfelLeDweXnzxxevabzzqdGmU6ac//elr1rvtttt4/vnnGR0dFcu7u7tfU+N0PcfetGkTL7zwwkU25xMTE/zhD39g9erVWCwWgNdYzZvNZoqKiggGg2KZyWQCzouaC7nzzjuJRqP8+7//+0XLH3nkERQKxWtqzI4fPy5q0wCGhoZ44YUX2LRpEyqV6rrGHa9Ru9Ae3uv1Xtai/VLuvPNOxsfHeeqpp8SySCTCo48+itlsZt26dZe9NkqlkkWLFgFcdH0uR9xd8aMf/Sg/+clPXndMMjIyMjJvLXJES0ZGRuZt5sUXX8Tj8XDPPfdc9u8rV64UzYvvv//+a96vxWJh7dq1/PCHPyQcDmO329m1axd9fX2vWfd73/seu3btYtWqVXz+858X4qWyspKGhoY3dF4/+MEP2L17N6tXr+YLX/gCarWaX/3qVwSDQX74wx+K9SoqKli/fj3Lli0jJSWF06dP8+yzz/KlL31JrBM3sPjKV77C5s2bUalUfOhDH+Luu+9mw4YNfPvb36a/v5/Fixeza9c
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdZXhU1/r38e9IJu7unkASYoQEkuBWrC1tKZCW4lJaatSo0FO3Q93wAAUq0ELx4JBABGJEIO6uE09Gnhd9sv9AgVI5p6fn7M91cbUkM3vW7GSG/Zu11n1LtFqtFpFIJBKJRCKRSCQS/Wmkf/UARCKRSCQSiUQikei/jRi0RCKRSCQSiUQikehPJgYtkUgkEolEIpFIJPqTiUFLJBKJRCKRSCQSif5kYtASiUQikUgkEolEoj+ZGLREIpFIJBKJRCKR6E8mBi2RSCQSiUQikUgk+pOJQUskEolEIpFIJBKJ/mRi0BKJRCKRSCQSiUSiP5kYtEQikUj0XyU2NhaJREJJSYnwtVGjRjFq1Ki/bEwikUgk+t8jBi2RSCQS/UJxcTGPPvooPj4+GBgYYGBggJ+fH4888giZmZl/9fD+Ul9++SUzZszAxcUFiUTCvHnz/uohiUQikeg/kPyvHoBIJBKJ/rPs37+fmTNnIpfLeeCBBwgKCkIqlXL58mV++OEHvvzyS4qLi3F1df2rh3rb4uLi/rRjvfvuu7S1tREeHk51dfWfdlyRSCQS/XcRg5ZIJBKJBIWFhcyaNQtXV1eOHz+Ovb39Nd9/9913+eKLL5BK//MWRHR0dGBoaHjD7ykUij/tcU6fPi3MZhkZGf1pxxWJRCLRf5f/vH8pRSKRSPSXee+99+jo6GDz5s2/CFkAcrmcxx57DGdn52u+fvnyZe677z4sLCzQ09MjLCyMn3766ZrbSCSSm/65ej/V7Ryrfx/W6dOnWb58OTY2Njg5Od30ed1oj9ann36Kv78/BgYGmJubExYWxo4dO371HLm6uiKRSH71diKRSCT63ybOaIlEIpFIsH//fry8vIiIiLjt+2RnZxMVFYWjoyPPP/88hoaGfPfdd9x9993s3r2b6dOnA7Bt27Zf3Pell16irq5OmBm63WP1W758OdbW1qxevZqOjo7bHvP69et57LHHuO+++3j88cfp7u4mMzOTpKQkYmJibvs4IpFIJBLdjBi0RCKRSASAUqmkqqqKu++++xffa2lpQaVSCX83NDREX18fgMcffxwXFxdSUlLQ1dUFfg5A0dHRPPfcc0I4evDBB6855vvvv09paSlbt27FysrqNx2rn4WFBcePH0cmk/2m53rgwAH8/f35/vvvf9P9RCKRSCS6XeLSQZFIJBIBPwct4Ib7jkaNGoW1tbXw5/PPPwegqamJEydOcP/999PW1kZDQwMNDQ00NjYyceJE8vPzqays/MXxTp48yapVq1ixYgVz5sz53cdavHjxbw5ZAGZmZlRUVJCSkvKb7ysSiUQi0e0QZ7REIpFIBICxsTEA7e3tv/je2rVraWtro7a29pqZqYKCArRaLS+//DIvv/zyDY9bV1eHo6Oj8PeKigpmzpxJVFQUH3zwwR86lru7+297kv/fc889x7FjxwgPD8fLy4sJEyYQExNDVFTU7zqeSCQSiUTXE4OWSCQSiQAwNTXF3t6erKysX3yvf8/W1UUrADQaDQBPP/00EydOvOFxvby8hP/v7e3lvvvuQ1dXl++++w65/P/+GfqtxwKE5Yu/1cCBA7ly5Qr79+/n8OHD7N69my+++ILVq1fz6quv/q5jikQikUh0NTFoiUQikUgwZcoUNmzYQHJyMuHh4b96ew8PDwB0dHQYN27cr97+scceIz09nTNnzmBra/uHjvVHGRoaMnPmTGbOnElvby/33HMPb775JqtWrUJPT+9f/vgikUgk+u8m7tESiUQikeDZZ5/FwMCABQsWUFtb+4vva7Xaa/5uY2PDqFGjWLt27Q2b99bX1wv/v3nzZtauXcvnn39+wxD3W471RzU2Nl7zd4VCgZ+fH1qtlr6+vj/tcUQikUj0v0uc0RKJRCKRwNvbmx07djB79mx8fX154IEHCAoKQqvVUlxczI4dO5BKpdf0rPr888+Jjo5m0KBBLF68GA8PD2prazl//jwVFRVkZGTQ0NDA8uXL8fPzQ1dXl6+//vqax50+fTqGhoa3daw/w4QJE7CzsyMqKgpbW1tyc3P57LPPmDJlirBX7Wb27dsnjKOvr4/MzEzeeOMNAO68804CAwP/lDGKRCKR6O9NDFoikUgkusZdd93FpUuXWLNmDXFxcWzatAmJRIKrqytTpkxh2bJlBAUFCbf38/PjwoULvPrqq8TGxtLY2IiNjQ0hISGsXr0a+LnARnd3Nzk5OUKVwasVFxdjaGh4W8f6MyxdupTt27fzwQcf0N7ejpOTE4899hgvvfTSr9539+7dbNmyRfh7WloaaWlpADg5OYlBSyQSiUQASLTXrwMRiUQikUgkEolEItEfIu7REolEIpFIJBKJRKI/mRi0RCKRSCQSiUQikehPJgYtkUgkEolEIpFIJPqTiUFLJBKJRCKRSCQSif5kYtASiUQikUgkEolEoj/Z36K8u0ajoaqqCmNjYyQSyV89HJHof5pWq6WtrQ0HBwek0r/HZzXie4hI9J/j7/geIhKJRL/H3yJoVVVV4ezs/FcPQyQSXaW8vPyaprX/ycT3EJHoP8/f6T1EJBKJfo+/RdAyNjYGfn5TNjEx+YtHIxL9b1MqlTg7Owuvy78D8T1EJPrP8Xd8DxGJRKLf428RtPqX+piYmIgXSSLRf4i/0xI88T1EJPrP83d6DxGJRKLfQ1wcLRKJRCKRSCQSiUR/MjFoiUQikUgkEolEItGfTAxaIpFIJBKJRCKRSPQnE4OWSCQSiUQikUgkEv3JxKAlEon+rc6cOcO0adNwcHBAIpGwZ8+eX73PqVOnCA0NRVdXFy8vL2JjY//l4xSJRCKRSCT6I8SgJRKJ/q06OjoICgri888/v63bFxcXM2XKFEaPHk16ejpPPPEEixYt4siRI//ikYpEIpFIJBL9fn+L8u4ikei/x6RJk5g0adJt3/6rr77C3d2dNWvWADBw4EDi4+P58MMPmThx4r9qmCKRSCT6C2i1Wjo7O9FoNCgUCnR0dJBKxXkB0d/T/8xvbmVLF1mVrVS2dP3VQxGJRL/B+fPnGTdu3DVfmzhxIufPn7/pfXp6elAqldf8uV29vb3MmzePf/zjH2zYsAG1Wo1arSY2Nlb4/40bN7Jx40Z6e3uFrwOo1WrWrVvH/Pnz6e3t/cWxrz5Ob28vCxcupKuri9jYWLq6upg/fz7r1q275nhXHx+gr6+Pp59+mu7ubmJjY+nt7WXjxo2sX7+e/fv38/bbb6NSqW743LRaLfn5+Wi12l+Mqbe3l02bNnH48GE0Gs0vnkv/8163bh0bN268ZkxXf2/t2rXMnTv3F8/j6vtef976z/vChQuF83aj534jNzqn15/7q5/jjX4+/eP78ssvhbFfP0aNRsOhQ4d44403rvm9uPp4Go2Gw4cPs2nTphs+x5v9DK4fY/+5uvpcq9Xqa34+NzpHXV1dREdH097efs2xrj7Gzc5pb2/vTX//rj+OVqvlypUrxMXF0dfXx8aNG1mzZg0GBgbMnj37hr/7ItHt0Gq19Pb20tvbS09PDx0dHSiVStrb24XXmEj0dyLRXv9u/x9IqVRiampKa2vr72o2WtnSxbg1p+nqU6OvI+PYypE4mun/C0YqEv33+6Ovx6tJJBJ+/PFH7r777pvexsfHh/nz57Nq1SrhawcPHmTKlCl0dnair//L1/I//vEPXn311V98/XbGvHDhQtRqNWfPnsXd3R1nZ2eio6OJjo7mvffeIyIigri4OIKCgjhy5AgbN27k/PnzzJkzhyVLllBdXc3p06cZOHAgS5YsAUClUrF9+3bmzJnDyJEjWbRoER0dHTzyyCO8++677Nu3j6lTp+Lp6Ymenh5NTU088MADJCYm4u3tTWFhIRs2bABg8uTJBAQEsGPHDo4ePcqyZcuYNGk
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9eXic5Xn2/ZtNs49Go5FG0miXLFmSbVm2vBuwwcaAE5aEbAWy0KRJ22xN+6bN2zTfmx7v2/RLS5t8adOkWSCkJC/QsAUMeMMYL7It2ZKsfV9GGmk0M5oZzb5+f7hz1xgbbDCY5fkdhw6wNM8z92zSfT7XdZ2nLJPJZJCQkJCQkJCQkJCQkJC4asiv9QIkJCQkJCQkJCQkJCTeb0hCS0JCQkJCQkJCQkJC4iojCS0JCQkJCQkJCQkJCYmrjCS0JCQkJCQkJCQkJCQkrjKS0JKQkJCQkJCQkJCQkLjKSEJLQkJCQkJCQkJCQkLiKiMJLQkJCQkJCQkJCQkJiauMJLQkJCQkJCQkJCQkJCSuMpLQkpCQkJCQkJCQkJCQuMpIQktCQkJC4n3FQw89hEwmY2JiQnxv27ZtbNu27ZqtSUJCQkLig4cktCQkJCQkXsP4+Dhf/vKXqaurQ6fTodPpaGxs5E//9E/p7u6+1su7ZkxPT/Pd736X9evXk5eXh9VqZdu2bezfv/9aL01CQkJC4l2G8lovQEJCQkLi3cWzzz7LJz7xCZRKJffccw/Nzc3I5XIGBgZ44okn+Ld/+zfGx8epqKi41ku9bPbu3XtVzvP000/z//6//y933nknn/nMZ0gmkzz88MPs3LmTX/7yl3zuc5+7KvcjISEhIfHeR5bJZDLXehESEhISEu8ORkdHaW5upry8nAMHDlBcXPyqnyeTSX784x9z1113UVZWdo1WeXFCoRB6vZ6HHnqIz33uc4yPj1NZWXlV76O3txebzYbVahXfi8VirF69mmAwyPT09FW9PwkJCQmJ9y5S66CEhISEhOD73/8+oVCIBx988DUiC0CpVPLVr371NSJrYGCAu+++G4vFgkajobW1lWeeeeZVt5HJZJf8On+e6nLOlZ3Devnll/mTP/kTCgsLKS0tveTjutiM1o9+9COamprQ6XTk5eXR2trKb37zm9d9fpqaml4lsgDUajW33XYbDoeDpaWl1z1eQkJCQuKDg9Q6KCEhISEhePbZZ6mtrWXDhg2XfUxvby9btmzBbrfzV3/1V+j1eh577DHuvPNOfve733HXXXcB8Otf//o1x37729/G5XJhMBiu6FxZ/uRP/oSCggK+853vEAqFLnvNP/vZz/jqV7/K3Xffzde+9jWi0Sjd3d2cOHGCP/iDP7js82SZm5sTs2wSEhISEhIgCS0JCQkJif8iEAgwOzvLnXfe+Zqf+Xw+ksmk+Lder0er1QLwta99jfLyck6dOoVarQbOCaCtW7fyl3/5l0Ic3Xvvva865z/8wz8wOTnJww8/LKpEl3uuLBaLhQMHDqBQKK7osT733HM0NTXx+OOPX9FxF2NkZIQnnniCj33sY1e8DgkJCQmJ9y9S66CEhISEBHBOaAGiunQ+27Zto6CgQHz967/+KwBer5eDBw/y8Y9/nKWlJdxuN263G4/Hw65duxgeHmZmZuY153vppZf41re+xVe+8hXuu+++N32uL3zhC29K3JjNZhwOB6dOnbriY88nHA7zsY99DK1Wy9///d+/pXNJSEhISLy/kCpaEhISEhIAGI1GAILB4Gt+9tOf/pSlpSXm5+dfVZkaGRkhk8nwN3/zN/zN3/zNRc/rcrmw2+3i3w6Hg0984hNs2bKFf/qnf3pL56qqqrqyB/lf/OVf/iX79+9n/fr11NbWcvPNN/MHf/AHbNmy5bLPkUql+OQnP0lfXx/PP/88JSUlb2otEhISEhLvTyShJSEhISEBQG5uLsXFxfT09LzmZ9mZrfNNKwDS6TQAf/EXf8GuXbsuet7a2lrx//F4nLvvvhu1Ws1jjz2GUvnff4au9FyAaF+8UhoaGhgcHOTZZ5/lhRde4He/+x0//vGP+c53vsN3v/vdyzrHF77wBZ599lkeeeQRbrzxxje1DgkJCQmJ9y+S0JKQkJCQEOzevZuf//znnDx5kvXr17/h7aurqwFQqVTs2LHjDW//1a9+lc7OTg4fPozNZntL53qr6PV6PvGJT/CJT3yCeDzORz7yEf7P//k/fOtb30Kj0bzusf/jf/wPHnzwQX7wgx/wqU996m1fq4SEhITEew9pRktCQkJCQvDNb34TnU7H/fffz/z8/Gt+fmH0YmFhIdu2beOnP/0pTqfzNbdfWFgQ///ggw/y05/+lH/913+9qIi7knO9VTwez6v+nZOTQ2NjI5lMhkQi8brH/sM//AP/+I//yP/8n/+Tr33ta1dtTRISEhIS7y+kipaEhISEhGDZsmX85je/4VOf+hT19fXcc889NDc3k8lkGB8f5ze/+Q1yufxVmVX/+q//ytatW1m5ciVf+MIXqK6uZn5+nuPHj+NwOOjq6sLtdvMnf/InNDY2olar+Y//+I9X3e9dd92FXq+/rHNdDW6++WaKiorYsmULNpuN/v5+/uVf/oXdu3eLWbWL8eSTT/LNb36TZcuW0dDQ8JrHsXPnztdU6iQkJCQkPphIQktCQkJC4lXccccdnD17lgceeIC9e/fyy1/+EplMRkVFBbt37+ZLX/oSzc3N4vaNjY20t7fz3e9+l4ceegiPx0NhYSEtLS185zvfAc4ZbESjUfr6+oTL4PmMj4+j1+sv61xXgy9+8Ys88sgj/NM//RPBYJDS0lK++tWv8u1vf/t1j8sKveHh4Ys+jpdeekkSWhISEhISAMgyF/aBSEhISEhISEhISEhISLwlpBktCQkJCQkJCQkJCQmJq4wktCQkJCQkJCQkJCQkJK4yktCSkJCQkJCQkJCQkJC4ykhCS0JCQkJCQkJCQkJC4iojCS0JCQkJCQkJCQkJCYmrzHvC3j2dTjM7O4vRaEQmk13r5UhIfKDJZDIsLS1RUlKCXP7euFYj/Q6RkHj38F78HSIhISHxZnhPCK3Z2VnKysqu9TIkJCTOY3p6+lWhte9mpN8hEhLvPt5Lv0MkJCQk3gzvCaFlNBqBc7+UTSbTNV6NhMQHm0AgQFlZmfhcvheQfodISLx7eC/+DpGQkJB4M7wnhFa21cdkMkmbJAmJdwnvpRY86XeIhMS7j/fS7xAJCQmJN4PUHC0hISEhISEhISEhIXGVkYSWhISEhISEhISEhITEVUYSWhISEhISEhISEhISElcZSWhJSEhISEhISEhISEhcZSShJSEh8Y5y+PBhPvzhD1NSUoJMJuOpp556w2MOHTrEmjVrUKvV1NbW8tBDD73t65SQkJCQkJCQeCtIQktCQuIdJRQK0dzczL/+679e1u3Hx8fZvXs327dvp7Ozk69//et8/vOf58UXX3ybVyohISEhISEh8eZ5T9i7S0hIvH+49dZbufXWWy/79j/5yU+oqqrigQceAKChoYEjR47wz//8z+zatevtWqaEhISExDUgk8kQDodJp9Pk5OSgUqmQy6W6gMR7E0loSbwjzPgiLIbi5OlzsJu113o5Eu8hjh8/zo4dO171vV27dvH1r3/9ksfEYjFisZj4dyAQuOz7y2QyDA0NMTExQVlZGW1tbWzatAm5XE46nebIkSPMzc1RXFzMvffey7e//W2+973voVKpyGQyDA8Pk8lkyGQyTExMiP8/c+YMLS0t1NTUsGzZMkZGRshkMgCkUimefvppPvzhD3Pq1CnuueceHnnkETZu3Egmk+GJJ54gk8nQ2trKzp07OXDgAJlMhsrKSuRyObW1ta86H5zLKKqurubhhx+mpKSEG2+8kX/8x3/EarVSWloqMoxkMhkVFRVMTEwwMzOD3W6nvLycX/7yl/z93/89MpmM73//+1gsFjweD3feeSczMzPceOONYh0VFRVMTU1RUVEB8KrH7XQ62bRpE9PT01RUVIj7TafTHD16lEwmQ1lZGTt37mRoaIgnn3yS1atXU1lZyfHjxykpKaG
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d3xc9Znv/57eRxr1Lo2aJdmSJXfZxgYXDJhOYIGEALvZJDc32SWbm703v2xJyE2ySTbJ7mZTdtPYEEiAUA0YNzC427J67200KqORZkbT2+8P3/mubWxjG5t63q+XXuAzZ77ne86Zkb6f8zzP55HF4/E4EhISEhISEhISEhISElcM+fs9AQkJCQkJCQkJCQkJiY8aktCSkJCQkJCQkJCQkJC4wkhCS0JCQkJCQkJCQkJC4gojCS0JCQkJCQkJCQkJCYkrjCS0JCQkJCQkJCQkJCQkrjCS0JKQkJCQkJCQkJCQkLjCSEJLQkJCQkJCQkJCQkLiCiMJLQkJCQkJCQkJCQkJiSuMJLQkJCQkJCQkJCQkJCSuMJLQkpCQkJB434jFYixZsoRvf/vb7/mxh4eHkclk/PM///N7etyHHnqIoqKiM7bJZDK+8Y1viH8/9thjyGQyhoeHr+ix7733Xu65554rOqaEhISExLmRhJaEhITEVSSxYJbJZBw8ePBtr8fjcfLz85HJZNx8881nvCaTyfjiF78o/p0QBjKZjGefffZtY33jG99AJpPhcDjO2L5jxw42btxIRkYGer2e4uJi7rnnHl577TUArr32WjHuhX4SQqCoqEhsk8vlJCcnU11dzWc/+1mOHTt2SdfnD3/4A2NjY2ec5+nXTCaTodVqKS8v54tf/CJTU1OXNP7VJDG/z3zmM+d8/etf/7rY5+x78n7xv//3/+bZZ5+lpaXl/Z6KhISExEce5fs9AQkJCYmPA1qtlieffJL169efsf3NN99kfHwcjUZzSeM9+uij3Hnnnchksgvu98///M989atfZePGjXzta19Dr9fT39/P3r17+eMf/8gNN9zA17/+9TPEwokTJ/i3f/s3/r//7/+jsrJSbK+pqRH/X1tby1e+8hUAPB4PXV1dPPPMM/zyl7/ky1/+Mj/60Y8u6jx+8IMfcO+995KUlHTOc7RarQQCAQ4ePMjPf/5zXn31Vdrb29Hr9Rc1/tVGq9Xy7LPP8rOf/Qy1Wn3Ga3/4wx/QarUEAoEztv/yl78kFou9l9MU1NXVsWLFCn74wx/yu9/97n2Zg4SEhMTHBUloSUhISLwH3HTTTTzzzDP827/9G0rlf//qffLJJ1m+fPklRTxqa2tpbm7m+eef58477zzvfpFIhG9961ts3bqV3bt3v+316elpALZu3XrGdq1Wy7/927+xdetWrr322nOOnZuby6c+9akztn3ve9/j/vvv58c//jFlZWX8j//xPy54Hk1NTbS0tPDDH/7wnK/feOONrFixAoDPfOYzpKam8qMf/YgXX3yR++6774JjX4hYLEYoFLrs95/ODTfcwEsvvcTOnTu57bbbxPbDhw8zNDTEXXfd9bboo0qluiLHvhS8Xi8GgwGAe+65h3/8x3/kZz/7GUaj8T2fi4SEhMTHBSl1UEJCQuI94L777mN2dpY9e/aIbaFQiD/96U/cf//9lzTWvffeS3l5OY8++ijxePy8+zkcDtxuN+vWrTvn6xkZGZd03HdCp9Px+OOPk5KSwre//e0Lzg3ghRdeQK1Ws2HDhosaf9OmTQAMDQ0Bp6J1a9euJTU1FZ1Ox/Lly/nTn/70tvclUjCfeOIJFi9ejEajEWmTZxOPx/nsZz+LWq3mueeee8c55ebmsmHDBp588skztj/xxBNUV1ezZMmSt73nXDVaF8vOnTu55pprMBgMmEwmtm/fTkdHx9vGNxqNDAwMcNNNN2EymfjkJz8pXt+6dSter/eMz6KEhISExJVHEloSEhIS7wFFRUXU19fzhz/8QWzbuXMnLpeLe++995LGUigU/N3f/R0tLS08//zz590vIyMDnU7Hjh07cDqdlz33S8FoNHLHHXdgs9no7Oy84L6HDx9myZIlFx3hGRgYACA1NRWAf/3Xf6Wuro5HH32U73znOyiVSu6++25eeeWVt7339ddf58tf/jJ/9md/xr/+67+eU+hEo1Eeeughfve7371jtPB07r//fnbs2MHCwgJwKpL4zDPPXLKAficef/xxtm/fjtFo5Hvf+x5///d/T2dnJ+vXr3+baUYkEmHbtm1kZGTwz//8z9x1113itaqqKnQ6HYcOHbqi85OQkJCQOBMpdVBCQkLiPeL+++/na1/7Gn6/H51OxxNPPMHGjRvJycm5rLG+9a1v8eijj3LHHXecs1ZLLpfz1a9+lUcffZSCggI2bNjA+vXrueGGG1i2bNmVOKVzkojiDAwMsHjx4vPu193dzerVq8/7usvlwuFwEAgEOHToEI8++ig6nU6YhvT29qLT6cT+X/ziF1m2bBk/+tGP2L59+xlj9fT00NbWRlVVldh2ujiJRCJ86lOf4qWXXuKll17i+uuvv+jz/cQnPsEXv/hFXnjhBT71qU+xe/duHA4H9913H7/97W8vepwLsbCwwF/91V/xmc98hv/8z/8U2x988EEWLVrEd77znTO2B4NB7r77br773e++bSylUkl+fv47CmEJCQkJiXeHFNGSkJCQeI+455578Pv9vPzyy3g8Hl5++eXLjnqcHtV64YUXzrvfN7/5TZ588knq6urYtWsXX//611m+fDnLli2jq6vrMs/kwiTqfjwezwX3m52dxWKxnPf1LVu2kJ6eTn5+Pvfeey9Go5Hnn3+e3NxcgDNE1tzcHC6Xi2uuuYbGxsa3jbVx48YzRNbphEIh7r77bl5++WVeffXVSxJZABaLhRtuuEFEK5988knWrl1LYWHhJY1zIfbs2cP8/Dz33XcfDodD/CgUClavXs0bb7zxtvdcqEbOYrF8YJwQJSQkJD6qSBEtCQkJifeI9PR0tmzZwpNPPonP5yMajfKJT3zissf75Cc/KaJat99++3n3u++++7jvvvtwu90cO3aMxx57jCeffJJbbrmF9vZ2tFrtZc/hXCRS6Ewm0zvue6E6rp/+9KeUl5ejVCrJzMxk0aJFyOX//Xzw5Zdf5v/+3/9Lc3MzwWBQbD9XdM9qtZ73ON/97ndZWFhg586d5zX/eCfuv/9+HnjgAUZHR3nhhRf4/ve/f1njnI++vj7gv+vUzsZsNp/xb6VSSV5e3nnHi8fj7+hYKSEhISHx7pCEloSEhMR7yP33389f/uVfMjk5yY033khycvJlj5WIaj300EO8+OKL77i/2Wxm69atbN26FZVKxX/9139x7NgxNm7ceNlzOBft7e0AlJaWXnC/1NRU5ubmzvv6qlWrhOvg2Rw4cIBbb72VDRs28LOf/Yzs7GxUKhW//e1v32ZMAWdGv85m27ZtvPbaa3z/+9/n2muvvSzheeutt6LRaHjwwQcJBoNXvClwwg7+8ccfJysr622vn+5kCaDRaM4QpWczNzdHWVnZFZ2jhISEhMSZSKmDEhISEu8hd9xxB3K5nKNHj14Rs4RPfepTlJaW8s1vfvMdXf5OJyFg7Hb7u57D6SwsLPD888+Tn59/Rg+uc1FRUSEcBC+VZ599Fq1Wy65du/jzP/9zbrzxRrZs2XJZY61Zs4YXXniBw4cPc/fddxOJRC55DJ1Ox+23387+/fvZunUraWlplzWX81FSUgKcMjjZsmXL234uJRIXiUQYGxt7x/sjISEhIfHukISWhISExHuI0Wjk5z//Od/4xje45ZZb3vV4iahWc3MzL7300hmv+Xw+jhw5cs737dy5E4BFixa96zkk8Pv9PPDAAzidTr7+9a+/Y2pafX097e3tZ6T9XSwKhQKZTEY0GhXbhoeHL1ivdiG2bNnCH//4R1577TUeeOCBy2oo/L/+1//iH//xH/n7v//7y5rDhdi2bRtms5nvfOc7hMPht70+MzNz0WN1dnYSCARYu3btlZyihISEhMRZSKmDEhISEu8xDz744BUdL1Gr1dzcfMZ
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9d3Rc93nnj7+mVwxmUGYw6J0AUQiwgKRISmwSVW1Jliz3nqzjk3zjknU2yUliJ9l1dh0nysZyyU9JZFuKZVvF6hYpURTFBoJoJHqvAwyAwTRMb78/uPMJSQEUKVEWJd/XOTgkLm753DtzZz7v+zzP+5GlUqkUEhISEhISEhISEhISEtcM+Xs9AAkJCQkJCQkJCQkJiQ8aktCSkJCQkJCQkJCQkJC4xkhCS0JCQkJCQkJCQkJC4hojCS0JCQkJCQkJCQkJCYlrjCS0JCQkJCQkJCQkJCQkrjGS0JKQkJCQkJCQkJCQkLjGSEJLQkJCQkJCQkJCQkLiGiMJLQkJCQkJCQkJCQkJiWuMJLQkJCQkJCQkJCQkJCSuMZLQkpCQkJCQkJCQkJCQuMZIQktCQkJC4m3zyCOPIJPJxI9SqaSgoIDPfe5zzM7OrrrND37wA2QyGVu3bl1zv+n9fe9731vzmGfOnBHLvvWtb100Dr1eT3FxMXfddRf/8R//QSQSedN+du/eTX19/arHX1paQiaT8a1vfestroCEhISEhMTqKN/rAUhISEhIvP/5m7/5G8rKygiHw5w6dYpHHnmEY8eO0dPTg1arvWjdxx57jNLSUk6fPs3IyAiVlZVr7ve73/0uf/AHf4Ber7+icfzwhz/EaDQSiUSYnZ3l5Zdf5gtf+AIPPvggzz//PEVFRe/oPCUkJCQkJK4UKaIlISEhIfGOue222/jUpz7Fl770JR5++GH+5E/+hNHRUZ599tmL1hsfH+fEiRP84z/+I7m5uTz22GNr7rOpqQmn08mPfvSjKx7Hfffdx6c+9Sm++MUv8ld/9VccP36cRx99lJ6eHu6///63fX4SEhISEhJXiyS0JCQkJCSuObt27QJgdHT0ouWPPfYYFouFO+64g/vuu++yQmvHjh3s3buX//N//g+hUOhtj+WTn/wkX/rSl2htbeXQoUNvez8SEhISEhJXgyS0JCQkJCSuORMTEwBYLJaLlj/22GPce++9qNVqPv7xjzM8PExbW9ua+/nWt76F0+nkhz/84Tsaz6c//WkADh48+I72IyEhISEhcaVIQktCQkJC4h3j9XpZWlpiZmaGJ598km9/+9toNBruvPNOsU57ezsDAwN87GMfA2Dnzp0UFhZeNqq1a9cu9uzZw3e/+913FNVKm15cGmGTkJCQkJB4t5CEloSEhITEO2b//v3k5uZSVFTEfffdh8Fg4Nlnn6WwsFCs89hjj2Gz2dizZw9w3lnwgQce4PHHHyeRSKy5729961vMz89fVa3WpRiNRgD8fv/b3oeEhISEhMTVIAktCQkJCYl3zEMPPcShQ4d44oknuP3221laWkKj0Yi/JxIJHn/8cfbs2cP4+DgjIyOMjIywdetWnE4nr7766pr7vvHGG9mzZ887qtVaWVkBICMj46q2k8lkb+t4EhISEhISkr27hISEhMQ7pqWlhc2bNwNw9913s3PnTj7xiU8wODiI0Wjk8OHDzM3N8fjjj/P444+/afvHHnuMW265Zc39//Vf/zW7d+/mxz/+MWaz+arH19PTA3CRlbxWq11TuAWDQbGOhISEhITE20GKaElISEhIXFMUCgXf+c53cDgcfP/73wfOCymr1cqvfvWrN/18/OMf5+mnn75stOqmm25i9+7d/O///b/fVlTrZz/7GQAHDhwQy0pKSpienl51f4ODg2IdCQkJCQmJt4MktCQkJCQkrjm7d++mpaWFBx98kFAoxFNPPcWdd97Jfffd96afP/zDP8Tv97+p59alpGu1/vVf//WqxvKf//mfPPzww2zfvp19+/aJ5bfffjuxWIwf//jHF62fTCb54Q9/iFqtvmh9CQkJCQmJq0FKHZSQkJCQeFf47//9v3P//ffzf/7P/8Hv9/OhD31o1fW2bdsmmhc/8MADa+7vpptu4qabbuL1119fc50nnngCo9FINBpldnaWl19+mePHj7NhwwZ+9atfXbTuXXfdxS233MLXvvY1Tp8+zQ033EAwGOTZZ5/l+PHj/N3f/R25ublv7+QlJCQkJH7nkYSWhISEhMS7wr333ktFRQXf/va30Wq13HzzzauuJ5fLueOOO3jsscdwuVxkZ2evuc9vfetbwrVwNf7gD/4AOF9blZOTQ1NTE//+7//OJz7xiYvMOdLHffbZZ/n7v/97Hn/8cZ566imUSiUNDQ08+uijfPKTn3wbZy0hISEhIXEeWSqVSr3Xg5CQkJCQkJCQkJCQkPggIdVoSUhISEhISEhISEhIXGMkoSUhISEhISEhISEhIXGNkYSWhISEhISEhISEhITENUYSWhISEhISEhISEhISEtcYSWhJSEhISEhISEhISEhcY94X9u7JZBKHw0FGRgYymey9Ho6ExO80qVQKv99Pfn4+cvn741mN9BkiIXH98H78DJGQkJB4O7wvhJbD4aCoqOi9HoaEhMQFTE9PU1hY+F4P44qQPkMkJK4/3k+fIRISEhJvh/eF0MrIyADOfyibTKb3eDQSEr/b+Hw+ioqKxH35fkD6DJGQuH54P36GSEhISLwd3hdCK53qYzKZpEmShMR1wvspBU/6DJGQuP54P32GSEhISLwdpORoCQkJCQkJCQkJCQmJa4wktCQkJCQkJCQkJCQkJK4xktCSkJCQkJCQkJCQkJC4xkhCS0JCQkJCQkJCQkJC4hojCS0JCYnfKkePHuWuu+4iPz8fmUzGr3/967fc5siRI2zcuBGNRkNlZSWPPPLIuz5OCQkJCQkJCYl3giS0JCQkfqsEAgE2bNjAQw89dEXrj4+Pc8cdd7Bnzx66urr46le/ype+9CVefvnld3mkEhISEhISEhJvn/eFvbuEhMQHh9tuu43bbrvtitf/0Y9+RFlZGd/73vcAqK2t5dixY/zTP/0TBw4ceLeGKSEhISHxHpBKpQgGgySTSdRqNSqVCrlcigtIvD+RhJaExNtk1hPCHYgCYDGoKTDr3uMRfTA5efIk+/fvv2jZgQMH+OpXv7rmNpFIhEgkIn73+XxXdcxUKsXIyAgVFRWMjIyQTCYZHx/H4XCwY8cOqqqqeOWVV0ilUhQXF3Py5EkKCgrYv38/r776KqlUitLSUtEnSCaTUVVVRSqV4tChQ+LvAJOTk5SWllJVVcXw8DATExMUFhby7LPPctddd9Ha2kp+fj5lZWUkk0l+/etf09jYyNzcHHK5nM9+9rOMjY2RTCaZmJhAJpOxd+9efvazn5Gfn8++ffvE//fv38+hQ4fo6uriG9/4BkeOHKGkpITKykoOHz7M3r17GRkZYWxsjJmZGQCKioq4+eabGRsbo6ysjJ/+9KfY7XZKS0uZmJjA4XCwfft2JiYm6OjoAODuu+/m9OnTbNu2jampKWQyGTfffDMAhw4dIpFIMDs7i0wmY+vWrfzkJz/hf/7P/8nrr79OcXExAGNjY3R2drJx40Zx7s888wwf+tCHmJqaYnZ2loKCAjGO9O/l5eVUVVUxMjJCIpHg5MmTbN++HYBTp07xqU99isOHDwOwd+9eHn30UbZv305VVRUHDx6ks7OTD33oQzz77LPcc8891NTUkEgk+Id/+Af+5E/+BJlMxk9+8hPsdjtlZWUATE1NsXfvXg4dOkRnZyd33XUXjzzyCJ///OeZnp4WY0skErz66qvs27cPpVIp3gPp99b27duZnp6muLiYyclJ8Vo++uijbNu2DblcTjKZ5NSpU3z6059mZGSEkydP8pnPfAaZTMYrr7xCSUkJ1dXVAAwPD5NMJpmcnKSkpIRUKiW2HRsbIx6P8/TTT2O1Wrnhhhs4ffo0n/70pxkfH6e8vJxXX32V4uLiN72PAYaGhpicnGTfvn2MjY1RWVkp9cWSeMekUimi0SjRaJREIkEsFkMul6N
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1oAAAEqCAYAAAAMHyF2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9d3hc53mnAd/TMDOYjjYY1EEjCkkQLGAXSYmkJVrNRba8iS1HsRI72WQ/K944sXfX6W0Tf2sntiOXdWwVb2Q5liyrkyIpsZMAiN47MABmBoPpvX1/6Js3JEVKlESKKue+Ll4SBueceU/F+zvP8/weWTabzSIhISEhISEhISEhISFxzZDf6AFISEhISEhISEhISEh80JCEloSEhISEhISEhISExDVGEloSEhISEhISEhISEhLXGEloSUhISEhISEhISEhIXGMkoSUhISEhISEhISEhIXGNkYSWhISEhISEhISEhITENUYSWhISEhISEhISEhISEtcYSWhJSEhISEhISEhISEhcYyShJSEhISEhISEhISEhcY2RhJaEhISEhISEhISEhMQ1RhJaEhISEhJvm5/85CfIZDI0Gg0Oh+N1v9+zZw9r1qy5ASOTkJCQkJC4sUhCS0JCQkLiHROPx/n7v//7Gz0MCQkJCQmJ9wyS0JKQkJCQeMe0tbXxwx/+kIWFhRs9FAkJCQkJifcEktCSkJCQkHjHfP3rXyedTr9hVGt6ehqZTMZPfvKT1/1OJpPx53/+5+LnP//zP0cmkzE6OspnP/tZTCYTxcXF/K//9b/IZrPMzc1x9913YzQaKS0t5Zvf/OZF2zt69CgymYzHH3+cr3/965SWlqLT6bjrrruYm5sTy/3Zn/0ZKpUKt9v9ujH97u/+LmazmVgs9tYPiISEhITEhx5JaElISEhIvGNqamq47777rnlU69577yWTyfD3f//3bNmyhb/+67/mW9/6Fvv376e8vJx/+Id/oL6+nv/+3/87r7766uvW/5u/+RueffZZ/uRP/oT/9t/+GwcPHmTfvn1Eo1EAPve5z5FKpXj88ccvWi+RSPCLX/yCT37yk2g0mmu2PxISEhISHx4koSUhISEhcU34H//jf5BKpfiHf/iHa7bNzZs387Of/Yzf+73f41e/+hUVFRV85Stf4f777+d73/sev/d7v8czzzyDVqvlxz/+8evWX1lZ4fjx4zz44IP83d/9HQ8//DCjo6P88Ic/BKC+vp5t27bx6KOPXrTes88+i9fr5XOf+9w12xcJCQkJiQ8XktCSkJCQkLgm1NbW8rnPfY4f/OAHLC4uXpNtPvDAA+L/FQoFmzZtIpvN8oUvfEF8bjabaWxsZHJy8nXr33fffRgMBvHzPffcg81m47nnnrtomTNnzjAxMSE+e+yxx6isrGT37t3XZD8kJCQkJD58SEJLQkJCQuKa8T//5/8klUpdMwfCqqqqi342mUxoNBqKiope97nX633d+g0NDRf9LJPJqK+vZ3p6Wnx27733olareeyxxwDw+/0888wz/OZv/iYymeya7IeEhISExIcPSWhJSEhISFwzamtr+exnP3vZqNaVREs6nb7i9hQKxVV9BpDNZt/CSP8Ti8XCHXfcIYTWL37xC+LxOJ/97Gff1vYkJCQkJCRAEloSEhISEteYXFTr0loti8UCgM/nu+jzmZmZ6zaWsbGxi37OZrOMj49jt9sv+vy+++5jdHSUc+fO8dhjj7F+/XpWr1593cYlISEhIfHBRxJaEhISEhLXlLq6Oj772c/y/e9/n6WlJfG50WikqKjode6A3/ve967bWB5++GGCwaD4+Re/+AWLi4scOHDgouUOHDhAUVER//AP/8Arr7wiRbMkJCQkJN4xyhs9AAkJCQmJDx7/43/8Dx555BFGRkYuigw98MAD/P3f/z0PPPAAmzZt4tVXX2V0dPS6jaOgoICdO3dy//3343Q6+da3vkV9fT2/8zu/c9FyKpWKz3zmM3znO99BoVDwX/7Lf7luY5KQkJCQ+HAgRbQkJCQkJK459fX1l40KfeMb3+ALX/gCv/jFL/jqV79KOp3m+eefv27j+PrXv87tt9/O3/3d3/Htb3+bvXv38vLLL5Ofn/+6Ze+77z4A9u7di81mu25jkpCQkJD4cCDLvt3qYQkJCQkJifcoR48e5eabb+aJJ57gnnvuuap1enp6aGtr4+GHH5b6Z0lISEhIvGOkiJaEhISEhATwwx/+EL1ezyc+8YkbPRQJCQkJiQ8AUo2WhISEhMSHml//+tcMDg7ygx/8gD/4gz9Ap9Pd6CFJSEhISHwAkISWhISEhMSHmj/8wz/E6XTy0Y9+lL/4i7+40cORkJCQkPiAINVoSUhISEhISEhISEhIXGOkGi0JCQkJCQkJCQkJCYlrzPsidTCTybCwsIDBYEAmk93o4UhIfKjJZrMEg0HKysqQy98f72qkZ4iExHuH9+MzREJCQuLt8L4QWgsLC1RWVt7oYUhISFzA3NwcFRUVN3oYV4X0DJGQeO/xfnqGSEhISLwd3hdCy2AwAK89lI1G4w0ejYTEh5tAIEBlZaW4L98PSM8QCYn3Du/HZ4iEhITE2+F9IbRyqT5Go1GaJElIvEd4P6XgSc8QCYn3Hu+nZ4iEhITE20FKjpaQkJCQkJCQkJCQkLjGSEJLQkJCQkJCQkJCQkLiGiMJLQkJCQkJCQkJCQkJiWuMJLQkJCQkJCQkJCQkJCSuMZLQkpCQeFd59dVXufPOOykrK0Mmk/HUU0+96TpHjx5lw4YNqNVq6uvr+clPfnLdxykhISEhISEh8U6QhJaEhMS7SjgcZt26dXz3u9+9quWnpqa4/fbbufnmm+nu7ubLX/4yDzzwAC+++OJ1HqmEhISEhISExNvnfWHvLiEh8cHhwIEDHDhw4KqXf+ihh6ipqeGb3/wmAM3NzRw/fpz/83/+D7feeuv1GqaEhISExA0gm80SiUTIZDLk5eWhUqmQy6W4gMT7E0loSUi8h3H4onjDCSy6PMrN2hs9nBvCqVOn2Ldv30Wf3XrrrXz5y1++4jrxeJx4PC5+DgQCb/v7s9ks4+Pj1NbW8vLLL1NZWcnMzAwymYz9+/eTzWZ5+OGH2bp1KwqFgoaGBmQyGel0mp/+9KfYbDZqa2upq6vj4YcfxmazIZPJsNvtoo+QTCajvr6esbExpqensdvtNDQ0MDo6yokTJygrK6OmpgaA6elpZDIZe/fu5ciRI+zevZt/+qd/oqioiMrKSuRyOfv27WNiYoJMJsPs7Cx79+5lfHycmZkZ9u3bh0wmY2xsjFQqxVNPPcW6detQKBTY7faLvmPnzp3ceuutHDx4kLy8PP7t3/6NpaUlNmzYQHV1NU899RQlJSXs2LGD+fl59u7dy+joKE899RT//b//d6ampshms6TTaZ5++mm+8pWvMD09LY5ldXU1q1atIpPJ8NOf/pSysjL279/P+Pg4U1NTAOI4yWQy6urqOHToEAD79+8X+5HNZi86jrnjZrPZUCgU4ljdcsstTE5OUlNTwyOPPMK2bdtoaGjg8OHD7N27F4CXX36ZvXv3IpfLyWazjIyMcOzYMQYGBli9ejW/9Vu/xeHDhwG45ZZbePTRR9myZQvT09P09PSI/c4d29bWVrq7uyktLeW+++7jscceY+vWrQCcPHmSsrIyqqurOXXq1EX/b7PZkMvlVFVVcfr0abZt2ybGNDMzQ2Vl5UXL1dTUUFtbyz/+4z9SUFDA8vIyGzZsQC6XX3QMc9dn7rqur68XP194LHPHs7a2lkceeYQtW7YwNzdHdXU1crlcXK8Xnie5XC62LyHxdshmsyQSCRKJBOl0mmQyiVwuR6lUkpeXh1KplESXxPsKWfbCp+p7lEAggMlkwu/3S81GJT40OHxR9n3zFaLJNFqVgkNf2f2eEFvX8n6UyWQ8+eSTfOxjH7viMqtWreL+++/na1/7mvjsueee4/bbbycSiaDVvv6Y/Pmf/zl/8Rd/8brPr2bMmUyGgwcPks1msdvtTE5O8txzz6FQKKiqqsLv95NMJslkMiwtLQHQ0NDA4OAgDQ0
      "text/plain": [
       "<Figure size 1000x300 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make some plots to compare different choices of LCG parameters\n",
    "\n",
    "params = []\n",
    "# the parameters are:   m,  a,  c,  seed,  name\n",
    "params.append([2**16+1, 75, 74, 33, \"ZX81\"])\n",
    "params.append([2**32, 1664525, 1013904223, 33, \"Numerical recipes\"])\n",
    "params.append([2**17, 1277, 0.0, 13337, \"Anagnostopoulos 1\"])\n",
    "params.append([2**31-1, 7**5, 0.0, 323412, \"Anagnostopoulos 2\"])\n",
    "params.append([2**32, 2**18 + 1, 7, 33, \"Gezerlis 1\"])\n",
    "params.append([2**32, 1812433253, 7, 33, \"Gezerlis 2\"])\n",
    "params.append([2**31-1, 7**5, 0, 33, \"MINSTD (Park Miller)\"])\n",
    "params.append([2**31, 65539, 0, 33, 'RANDU'])\n",
    "params.append([-1, -1, -1, 56, \"Numpy\"])\n",
    "\n",
    "for m, a, c, seed, name in params:\n",
    "    if m == -1:\n",
    "        samples = np.random.default_rng(seed).uniform(size = 10**5)\n",
    "    else:   \n",
    "        samples = np.fromiter(lcg(m, a, c, seed), float, count = 10**5)\n",
    "        samples = samples/m\n",
    "\n",
    "    plt.clf()    \n",
    "    fig = plt.figure(figsize=(10,3))\n",
    "\n",
    "    # 1D distribution\n",
    "    ax1 = fig.add_subplot(131)\n",
    "    ax1.hist(samples, density=True, bins=100, histtype = 'step')\n",
    "\n",
    "    # 2D x_i against x_{i+1}\n",
    "    ax2 = fig.add_subplot(132)\n",
    "    ax2.plot(samples[1:], samples[:-1], 'k.', ms=0.2)\n",
    "    plt.title(name)\n",
    "\n",
    "    # 3D x_i, x_{i+1} and x_{i+2}\n",
    "    ax3 = fig.add_subplot(133, projection='3d')\n",
    "    \n",
    "    # zoom in on part of the data\n",
    "    #cutoff = 0.05\n",
    "    #samples = samples[np.where(samples<cutoff)]\n",
    "    #ax3.set_xlim(0,cutoff)\n",
    "    #ax3.set_ylim(0,cutoff)\n",
    "    #ax3.set_zlim(0,cutoff)\n",
    "    \n",
    "    ax3.scatter(samples[:-2], samples[1:-1], samples[2:], marker='.')\n",
    "    ax3.view_init(-140, 60)    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "161c9b75",
   "metadata": {},
   "source": [
    "Notes:\n",
    "\n",
    "- You can often spot the bad generators from the bin-to-bin variations of the 1D distribution. (E.g. ZX81, Anagnostopoulos 1 and Gezerlis 1 all have smaller than expected variance.) One other thing we could look at is how the bin-to-bin variations scale with the number of samples (expect the standard deviation from a uniform distribution to go down as $1/\\sqrt{N}$.\n",
    "- In other cases, the 1D distribution looks acceptable, but you can see correlations in the 2d plot, ie. successive samples are correlated. In others, e.g. RANDU, these emerge in the 3D plots.\n",
    "- These generators all have the property that when used to make a $n$-d plot, the points do not fill up the whole space, but instead lie on $(n-1)$d hyperplanes (so 2d planes in the 3d plot or lines in the 2d plot). There are at most about $(n! m)^{1/n}$ planes, which is not that large, e.g. $m=2^{31}$ then for $n=3$ we get $(6m)^{1/3}\\sim 2300$. \n",
    "- A related issue with these generators is that often the least significant bits are much less random then the most significant bits. E.g. change the 2D plots to log-log and look at the behavior of small numbers, for example in MINSTD."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b6736b6",
   "metadata": {},
   "source": [
    "**Further reading**\n",
    "\n",
    "- [Linear congruential generator](https://en.wikipedia.org/wiki/Linear_congruential_generator)\n",
    "- More on [RANDU](https://en.wikipedia.org/wiki/RANDU), which was used a lot in the 1960s and 1970s, but as we saw above is not a good generator, with significant correlations.\n",
    "- [Remarks on choosing and implementing random number generators](https://www.firstpr.com.au/dsp/rand31/p105-crawford.pdf)\n",
    "- Documentation for NumPy's random sampling routines can be found at [`numpy.random`](https://numpy.org/doc/stable/reference/random/index.html) including the [PCG64](https://numpy.org/doc/stable/reference/random/bit_generators/pcg64.html#numpy.random.PCG64) generator which is the default generator in numpy. Note that is has a period of $2^{128}$ ($\\sim 10^{38}$), much larger than the examples above!\n",
    "- Having a high quality random number generator is particularly important for cryptography. A starting point is the wikipedia article on [Cryptographically secure pseudorandom number generator](https://en.wikipedia.org/wiki/Cryptographically_secure_pseudorandom_number_generator).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4b2e1f5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}