{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "262a5ec8-2553-4237-ab62-319b6ca22089",
   "metadata": {},
   "source": [
    "# Example-45: Twiss (Beta functions sensitivity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f0ee0901-68d3-44d9-b416-cda8902591d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this example effect of systematic quadruple errors on Twiss parameters is illustrated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b59da5a2-9129-4ccd-8891-8b189daaec2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "\n",
    "from random import random\n",
    "from pprint import pprint\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "matplotlib.rcParams['text.usetex'] = True\n",
    "\n",
    "from model.library.line import Line\n",
    "\n",
    "from model.command.external import load_sdds\n",
    "from model.command.external import load_lattice\n",
    "from model.command.build import build\n",
    "from model.command.tune import tune\n",
    "from model.command.twiss import twiss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e9a3332a-0f78-4886-a240-e102bdae584f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load ELEGANT twiss\n",
    "\n",
    "path = Path('ic.twiss')\n",
    "parameters, columns = load_sdds(path)\n",
    "\n",
    "nu_qx:Tensor = torch.tensor(parameters['nux'] % 1, dtype=torch.float64)\n",
    "nu_qy:Tensor = torch.tensor(parameters['nuy'] % 1, dtype=torch.float64)\n",
    "\n",
    "psi_qx:Tensor = torch.tensor(parameters['dnux/dp'], dtype=torch.float64)\n",
    "psi_qy:Tensor = torch.tensor(parameters['dnuy/dp'], dtype=torch.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1c3b6e50-1d2b-41da-9d1c-c70711318ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build and setup lattice\n",
    "\n",
    "# Note, sextupoles are turned off and dipoles are linear \n",
    "\n",
    "# Load ELEGANT table\n",
    "\n",
    "path = Path('ic.lte')\n",
    "data = load_lattice(path)\n",
    "\n",
    "# Build ELEGANT table\n",
    "\n",
    "ring:Line = build('RING', 'ELEGANT', data)\n",
    "ring.flatten()\n",
    "\n",
    "# Merge drifts\n",
    "\n",
    "ring.merge()\n",
    "\n",
    "# Split BPMs\n",
    "\n",
    "ring.split((None, ['BPM'], None, None))\n",
    "\n",
    "# Roll lattice start\n",
    "\n",
    "ring.roll(1)\n",
    "\n",
    "# Set linear dipoles\n",
    "\n",
    "for element in ring:\n",
    "    if element.__class__.__name__ == 'Dipole':\n",
    "        element.linear = True\n",
    "\n",
    "# Split lattice into lines by BPMs\n",
    "\n",
    "ring.splice()\n",
    "\n",
    "# Set number of elements of different kinds\n",
    "\n",
    "nb = ring.describe['BPM']\n",
    "nq = ring.describe['Quadrupole']\n",
    "ns = ring.describe['Sextupole']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2002817d-2e07-4195-8067-28da06a5aa8a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# Compute tunes (fractional part)\n",
    "\n",
    "guess = torch.tensor(4*[0.0], dtype=torch.float64)\n",
    "nuqx, nuqy = tune(ring, [], alignment=False, matched=True, guess=guess, limit=8, epsilon=1.0E-9)\n",
    "\n",
    "# Compare with elegant\n",
    "\n",
    "print(torch.allclose(nu_qx, nuqx))\n",
    "print(torch.allclose(nu_qy, nuqy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3942deef-8682-4e3e-9f1b-f58134d9a2ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([17])\n",
      "torch.Size([17])\n"
     ]
    }
   ],
   "source": [
    "# Compute nominal beta functions at BPMs\n",
    "\n",
    "kn = torch.zeros(nq, dtype=torch.float64)\n",
    "\n",
    "_, bqx, _, bqy = twiss(ring, [kn], ('kn', ['Quadrupole'], None, None), matched=True, advance=True, limit=1, epsilon=None).T\n",
    "\n",
    "print(bqx.shape)\n",
    "print(bqy.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "472c7222-a00c-4666-bfe3-5360e27055b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([17, 8192])\n",
      "torch.Size([17, 8192])\n"
     ]
    }
   ],
   "source": [
    "# Compute beta functions at BPMs using MC\n",
    "\n",
    "kns = 0.01*torch.randn((8192, nq), dtype=torch.float64)\n",
    "\n",
    "_, bqxs, _, bqys = torch.vmap(lambda kn: twiss(ring, [kn], ('kn', ['Quadrupole'], None, None), matched=True, advance=True, limit=1, epsilon=None), chunk_size=1024)(kns).swapaxes(0, -1)\n",
    "\n",
    "print(bqxs.shape)\n",
    "print(bqys.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fce256c4-0643-48ca-9ab2-4eab7ae695b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([17, 8192])\n",
      "torch.Size([17, 8192])\n"
     ]
    }
   ],
   "source": [
    "# Set beta beating\n",
    "\n",
    "dqxs = 100*((bqx - bqxs.T)/bqx).T\n",
    "dqys = 100*((bqy - bqys.T)/bqy).T\n",
    "\n",
    "print(dqxs.shape)\n",
    "print(dqys.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5021ad7b-5384-43e0-94d3-2c16789ec07e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAADT4AAAC+CAYAAACm5fGbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzj0lEQVR4nO3dv24cR54H8J8FA4xWnuVli1My9wY86gmOAjaZjDw9gehcwRKKDooE+Q3kfQKbzBSa+wSU+AaaxMCGFE1HxALLCwYzniFnmt0z1dPV3Z8PIEBky1arWN+p6vrX39zd3d0FAAAAAAAAAAAAAAAAQEaeNH0DAAAAAAAAAAAAAAAAAPfZ+AQAAAAAAAAAAAAAAABkx8YnAAAAAAAAAAAAAAAAIDs2PgEAAAAAAAAAAAAAAADZsfEJAAAAAAAAAAAAAAAAyI6NTwAAAAAAAAAAAAAAAEB2bHwCAAAAAAAAAAAAAAAAsmPjEwAAAAAAAAAAAAAAAJCdb5u+gSL//ve/45///Gf86U9/im+++abp24GZu7u7+P333+Mvf/lLPHnSjv2D8kSO2pilCHkiT23MkyyRozZmKUKeyFMb8yRL5KiNWYqQJ/LUxjzJEjlqY5Yi5Ik8yROk08Y8yRI5amOWIuSJPLUxT7JEjtqYpQh5Ik9tzJMskaM2ZilCnshTlTxlvfHpn//8Zzx79qzp24CVfv311/jP//zPpm+jFHkiZ23KUoQ8kbc25UmWyFmbshQhT+Rtkzydn59HRMT19XVcXFzEy5cvY29vLyIixuNxnJ2dxXA4jPF4HMfHxzEYDB69VkSWyJm2CdJpU55kiZy1KUsR8kTe5AnSaVOeZImctSlLEfJE3tqUJ1kiZ23KUoQ8kbc25UmWyFmbshQhT+StTJ6y3vj0pz/9KSIm/5CnT582fDfwh5ubm3j27NmsjraBPJGjNmYpQp7IUxvzJEvkqI1ZipAn8pQiT0dHR/GPf/wjDg4O4urqKo6OjuLLly+za58/f46IyUanV69exenp6aPXisgSOdI2QTptzJMskaM2ZilCnsiTPEE6bcyTLJGjNmYpQp7IUxvzJEvkqI1ZipAn8tTGPMkSOWpjliLkiTxVyVPWG5+mr1F7+vSpgJGlNr3qT57IWZuyFCFP5G2TPG37rRqyRM60TZDOJnk6PT2dtUURsdD2zBsOh7N2rOha2XuVJXKkbYJ02pQnWSJnbcpShDyRN3mCdNqUJ1kiZ23KUoQ8kbc25UmWyFmbshQhT+StTXmSJXLWpixFyBN5K5OnrDc+AQCwPdt+qwYAPObg4GD2+9PT0/j+++8jYrJZd3d3d+HP7u7uxuXlZXz69GnltflNVAAAAADztn04GAAAAABQjo1PAABExPbfqgEAZVxeXsZPP/0UL168iOPj44iYLEBa5urqqvDafbe3t3F7ezv7+ubmZuP7BQAAANrJ4WAAAAAAkKcnTd8AAAB5WOetGkXXACCFvb29ePPmTXz58iXOzs4K/+yqTU+rrr179y6+++672a9nz55teLcAAABAWzkcDAAAAADy5I1PAADMeKsGbG66sOH6+jouLi7i5cuXswUT4/E4zs7OYjgcxng8juPj44UFFKuuQd8NBoM4OjqKFy9exNevX2MwGDxoa66urmIwGBReu+/Nmzfx+vXr2dc3Nzc2PwEAAEBPrXM42KdPn1Zem99EBQAAAACszxufemI0mvwCNidLkIa2KU/eqtF+ctW8o6Oj2N3djcPDw/iv//qvODo6Wrj2t7/9LQ4PD+Pw8DBevXpV6hr10R7l6/z8PP785z/Pvh4OhxEx2SQ4vxBp3v7+fuG1+3Z2duLp06cLv6iXvEE12ilIS6agOpmBtLRF+bu8vIyTk5PaDge7ublZ+EW95A3KkxdIT98vH+fn53F+fh5nZ2dxcnISl5eXs2vj8Th++OGHODs7ix9++GGhf1d0je2SJUhPOwXpyNJ2eOMTAAALvFUDNnN6erpwmuv8G53mDYfD2duhiq5BX+3u7i5sYrq8vIzBYLD0tOTxeBz7+/uztmnVNQAAAIAie3t7MRwO4+TkJM7OzuLw8HDln13ncLC3b98muEsAAKo4OjqKf/zjH3FwcBBXV1dxdHQUX758mV37/PlzREzmlF69ehWnp6ePXgMAYLsqbXyaLry7vr6Oi4uLePny5WzB0Xg8jrOzsxgOhzEej+P4+Hhhgd+qawCwCW1Tu0x3tn/82Ox98ND5+XkcHR3F169fI+LhWzU+fPjw4L/Z39+P4XC48tp9Ozs7sbOzk/jOIT/zGzVOT0/j+++/j4hJznZ3dxf+7O7ublxeXsanT59WXlu2yYP6abOat7e3Fy9fvowff/wxIiJ++eWX2eRSxCRfJycn8fz587i4uFiYaCq6Rr1kBwAAgLZzOBhsxvwtALlxcCUAQPtV2vhk53u3jEYWIgHtp22CNLxVA9K6vLyMn376KV68eBHHx8cRsfoE2Kurq8Jry9ze3sbt7e3s65ubm43uF89HuZo/VXmapanhcBjv379/8Oceu0bzbI7aPguOgD6btjsAAEUcDgbpmL8FIDcOrgQAaL9KG5/sfAcgN9omSMNbNdpnfpOGDRv52dvbi+FwGCcnJ3F2dla48WLVpqeia+/evYu3b99ueJcA9IUFR+1kswYAAGyPw8Haz0Er+TB/C+k40AjScXAlpKFtAqAplTY+2fkOQG60TZCOt2pAWoPBII6OjuLFixfx9evXGAwGDwbCr66uZgskVl1b5s2bN/H69evZ1zc3N/Hs2bPk/waAFGzQbZ4FR+1iwxPQBxZItIuF5HmTJ0jD4WDdYRyieeZv20u/Lz8ONIJ0HFwJaWibAGhKpY1PEfXufLfrvTkG/4A2cyoLALk4Pz+Po6Oj+Pr1a0RMFolHTAbuDg4O4sOHDw/+m/39/RgOhyuvLbOzsxM7OzsJ7xyALrPgCIDcWCAB6cgTpONwMEjH/C2k4UAjSMvBle1mg24etE0AEw6z3L7KG5/q3Plu1/t2CFoenMAH6TiVBeirMv06g3/btbu7u7C4/PLyMgaDwdIF4uPxOPb392cD56uukReZAtrKgiNIw5gepGGBBKQjTwDkyPxte1hDlDcHGrWDuaO8ObgS0tI2AdCUyhufIurb+W7XO33iBD5Iy6ksAORgb28vXr58GT/++GNERPzyyy+zvlvEZODv5OQknj9/HhcXFwv9uKJrNMuEFdAFFhxBGsb0IA0LJCAdeQIgV+ZvIY06DzRymFF6y+aURiNzTE1ycCWkp21qp/n2yOb35jloD6orvfFpGzvf7Xqvn8YqH07gaweLXPPmVBYAcjS/kHw6yDc1HA7j/fv3D/7cY9cAIAULjmBzxvQgHW8jhHTkqd0shAW6xPwtpFfngUYOM6IPHFwJ6WmbYHMO2oPqnpT9g0U736cDFVPzu9uLrkGfrXMCX9E1mjMa2VTYFG0TAABAsfPz8/jzn/88+/r+gqNl9vf3C68ts7OzE0+fPl34BV1kTA/S2dvbizdv3sSXL1/i7Oys8M+u+zbC7777bvbLhly6TJ4AyIX5W6jH9ECjo6OjuL6+Tnag0Zs3b+K3336b/fr111/r/GdAYw4PD+P4+DiOj4/j9PR0od2ZHk55eHgY79+/X8hK0TXqZz1e3rRNsBkH7UF1pd/4ZOd7++j45c8rP7vHG6K2S9sE9JE+HtRHXw7S0mbloWjB0X3zi4ruTzZZcNQsbwHIizG99pOpfHgbYXt5fsqPPLVLmbZIewXrMybRHPO33aG/17xtvEHN29PSKtP+6ONBsbIZ0U41Q9sE6axz0N6nT59WXls29wtdU3rjU0QsvI5wOqE7Nd3dfv/PPXYN+swrP9vPoHnztE0AANAck7T5s+AI0jOmB5vZxgKJCIsk6Ad56iZzT0Dbmb+FNBxo1G5FfTobNoC20ja1j7ncvNV50F6Ew/bonkobn4D06jqBz+l7QJsY2IN6zA+oG8wAgP6y4AjSM6YH67NAAtKRJ6CPzClB/WzCzYMDjQCW0041R9sEadV50F6Ew/boHhufoAFe+Qn0nUEIAHKgPYL13d9cG2HBETTB5vbtMqYHaVgg0Q6el9pBnoC+0T4BfeNAI6Dv9P/yo23Kn7mjdqnroL0Ih+3RPTY+dZDOXv6cwNdOsgX1kC1IQ5YAyEmZdsmgO9A2xvS6xcbdZlkg0U5ykyd5AgAAAGAV64nys42D9iIctkf32PgEDXACHwBQJ4MWUExGAMiJTYDtYUyvPWzOgO2SOUjHmAUAbWE8AwCge/TxtsNBe+0lI82y8Qka4gS+djPxBNXJDQAAAG1nTA8AWJcxcmiezbqwGYv8AABgcw7ag/XY+AQAAJkzGQvNmV+UZIESFJMRAHJQtT2ycA82ow8IAADwOM9OAORAe5QPB+11l3V+9bHxqSNShETQAACArjOQB3kxFgHr0Z4BAAB9ZwM7bI9xCAAAgAnPR82x8QkgEY0ZAAAAAABdZywc0pMrWI/swHqqHEYkZwAAAOTAxqeeM0ABAAAAAO1hPA/KkxdYj+zAdsgaAADAQ1Wflaps6IW+MOYAaclUHp40fQMAAI/RcQQAyMc6fbPRSJ8OAAAAqM54AtRHvgAAusN8LNB13vgEcI9TIAAAuqfKAJ/BQACaUtQGaZ8AaAPj6wAAAMC2GDcHgP7wxieAEjwkwXbIGtTL6S4AAN2jjwdAKtoUyI9cAgAAAMB2rTMeZxyvft74BADUSmcOAAAAAIC+MkYOAMAyo5G35dJv82+N9twEzdEeQTFtVD688aljUoRLQAEA8qFvBgAAALSBMQwAgO5zijkAKWhLAICqvPEJAMjK/IkuANAGTkACIBeep2A99xdayBAAqVjMB2nIEmzOmAEAAGXoNwK58sYnAADoGZPEsDk5AiA32iYAAICHlr2dxvMTfTNf59V/2JwcAQDA9tn41AEepiA/cgmbkyP6TP0HAOiXZQvxAAAASMczFwAAQHmeoei7TedvZSi9b5u+AYCcaGggLzJJn6n/AORG2wRAzlK3U9P/38ePaf+/AAAAAFPGHyAtmQKA7vLGJwCgNhbHQv7kFACgGd7yBABAWfqNkJ5nMticHAEA5EGfDOgDG58AAAAAaJTBeAAAAAAA2sbYNgBNshEd6JNvm74BgFzpEAKQG20TAAAAAAAAAACQivVIQBvY+NRidTY00//3x4/1/R0AAAAAAAAAAADQJhaIAwDAdtn4BAAAAAAAACQzGjlcD1KbX1wrXwAA22ezE9QjVbZkFIDceBFNWjY+AQBJ6KQBAHSXySIAAIBybPyDzRmHAKBv9CEBAJpnPCJvT5q+AQAAAAAAgHWZiAIAYF36kgAAAMuNRp6ZgHzY+EQhDRYAQDcZnIC0ZAoAAGCR5yQAAAAAACAFG58AamJSFwCg2/T3AMiV9gmAnGiX6CNjBgAAAADQHsby8vdt0zcAAAAAsMzl5WWcn59HRMTFxUX8/e9/j8FgEBER4/E4zs7OYjgcxng8juPj41LXAACqGI0iPn5s+i6gPtPJXPUcAADgIc9MAACQBxufAAAAgCydn5/H3/72t4iI+OGHH+J//ud/4vPnzxERcXR0NPv9eDyOV69exenp6aPXAJpkoQQAAAAAAABAv5gn3pyNTwDhFYXQBjp+9fNWDQBycnl5Ge/evZttfDo8PIyTk5MYj8cP/uxwOJy1Yfevz18DyIVxCMibjAIAAABAesbdoL2s3QOaZuNTyzTRcGisAKjKQEU7easG0DWeZdptb28v/v73v8++vr6+joiI3d3d+Pnnn2N3d3fhz+/u7sbl5WV8+vRp5bW9vb3a75vNjUZyCwAAwPY5HAwAAAAA8mTjE8AWWHQL5M5bNSAdCyQgncPDw9nvf/rppzg4OIjBYDDbBHXf1dVV4bX7bm9v4/b2dvb1zc3NRvfbFfPPL9vc0G7zPAAAAE1yOBgA0CRriwAAYDUbnwBqZvEe1MfbANLxVg1IxwIJSO/6+jrOzs5mGSn6c1WuvXv3Lt6+fbvh3QEAAABt53AwACAX1kGk5+BKAID2q7TxSQcQ6BIDBdBOslsfb9WAzVkgAfU4OTmJX375ZTaWMBgMHrQ1V1dXMRgMCq/d9+bNm3j9+vXs65ubm3j27Fny+wcAAADy5nAwSMfaIgBy4+BKoIscyA/0TaWNTzqAAOTGwDmk560a/WVj4eYskGjedHBPXe6OH374IU5OTmI4HM4ydXBwEB8+fHjwZ/f392M4HK68dt/Ozk7s7Owkv2cAAGA5Yw9AzhwO1ixtRHdYWwRAThxc2U/6lnSZDU9AX5Xe+KQDmBcNF8CEgfP8aKPaz1s1+kdu06p7gUSERRL0x9nZWezt7c02Pf38889LN6yPx+PY39+ftU2rrtFONjRuzoERQBdZvABAboyxdYfDwWB91hZBWsb1YHMOroS0tE0ANKX0xicdQAByY+Ac0vNWDUinrgUSERZJlNHUYiMLcNMZj8dxdHS08L3BYBDHx8cREXF6ehonJyfx/PnzuLi4WNjEXnQN+siBEQDkxgKJ/NT9DGVDSH3kCdJzOBisz9oiSMu4HqTh4EpIR9sEdI2x6/YovfEpwqvdISUTUXnQYLWbgXNIy1s1IK26FkhEWCRBPwyHw7i7uyu8/v79+4hYHK947Br0jQMj8tHkptwIG3NTMaYHaVggAenIE6TlcDDYnMXlkIZxvbxYX9QNDq6EzWibAGhSpY1PU17tDpszEQVpGDiHNLxVIw8GzLujzgUSERZJANuhXeoGB0Y0r+kNT6RlTC8v6nk7WSCRDxlqP3nqL5vb6+FwsP6SqXpYXN4vcpSecT1Iz8GV/aJtSk/b1F/GEdNz0F6/yFAaa2188mp32IyJKEjPwHmzdMzaz1s1IB0LJPKjnQL6zoERkIYxPUjDAglIR54gHYeDQXoWl/eHMfj61D2uZ0yPPnFwJaShbYI0HLTXX6ORTbnrqrzxyavdYXMmoiA9A+dA25kU6gYLJADImQMjYDPG9PLg2akbbMqFdOQJ0nA4WHP077rJ4vL+soivHnWN6xnToy8cXAnpaZtgfQ7ag/VU2vikAwjp2PkO6Rg4ByAXFkgAkDMHRsDmjOlBWjblQjry1BybNgAWWVsE6dU1rmdMjz5wcGU+PDt1i7YJ1uegPVhP6Y1POoBQDzvfYTMGzgEAAB7nwAhIy5gepGFTLqQjT0AXWAzbftYWQXp1jusZ06MPHFzZLG8C7CZtE2zOW9yhutIbn3QAoR52vvfPdMDeQ93mDJw3y+AEAAB18dyUlgMjID1jev2iXaqHTbmQjjw1xyYNgEXWFuVFO9V+xvWALsihPTK+l462CdLyFncor/TGJ5hnsXsadr7DZgyc91cOgyKQivoMANTJgRGQnjE92JwFEpCOPAEA8yzs7g7jegDkRtsE6XmLO5Rn4xM0xERUP1ncDgAA5MwzS/c4MALSMqbXLO1UN1ggAenIEwAwb/6ZyfNT+xnXAyA32iZIy1vcoRobn6ABJqKI8OY0usGAOQAAAH1hTA/SsEAC0pEnAADoJm9vA4Buc9BeHqx/bRcbn6ABJqIA0jHgBwAAwDYY08NBRnSFyVwAAAAAgGY4aK95xsjbycYnShNyAAAAAAAAAAAAAACozkF7sJ4nTd8AAAAApOCwBugeuQYA6J7RSD8PUpElSEumAAAAAPJk4xMAAGxZzgt8cr43AACApnlmAgAAAACgCcamoRvMNa3HxqcWmFZuFRwAAAAAAACA+8wlQ1oyBQDQPfp4ANBe3zZ9AwAAAAAAAED9LPABAKBu0z7nx4/N3gcAQFcY0wPwxicAAAAAAAAAAAAAoGNsGAGAbrDxCQAAAAAAANgqC48AAAAAoF1GI+N6QDNsfAIAAAAA6DiTULSZiVQAAAAAAMiHMXtg275t+gYopmEAAAAo5rkJ+mE0ivj4sem7AAAAyItxEdpM/QUgR9P2yXg0AADkwxufAACAB0w4AwB0l74eAAAAAOTPOB4AQFo59a9yupc28MYnAKATvAEAAGB9BtQAAAAAAOAPxs1hffIDQG60Te1n4xPQeRorSMfmIgAAAAAA+s7cE2wm5wzlfG/QNuaWAQAASOVJ0zdAe41GBv1gU3IE0G/aANhc7jnS3wPoN20AAAAAALSH8TwAAMiTNz4BAABAQtNJMSdZAvRXrgskcr0vAPrL8xNtor4CAAAAADTDG58AAAAAAAAAAGCOt7nD5uQINidHAP2mDQCYsPEJAAAAamYwEgAAALrBMz4AOdEuAZAbbRMAlDfd5K79fJyNTwAAAFADgxIAAEAO2vBs0oZ7BAAAustiUwAAyJuNTwAAsAUGywEA+kGfD/pB1gHIkTFIgH7y+Q8AAEDXfdv0DQAA7WLQHAAAAAAAAAAAmLKeCPpnmvuPH5u9D6AfbHzKlE4gAAAAAAAAAAAAAKRhowYAtJONTwAAAAAAANBhDtwDAABot9HIRg2APjGeB7DIxicA4FEepAAAAAAWOR0WAAAAAAAA6vek6RsAAEhlNLJJi/ypowDkRtsEQM60UwAAAAAT1kQAANBX3vgEkAGvowYAqMakDgC5aVvb5E01AORG2wQAUF7bxiEi9PcAAABYnzc+AQAAAAAAAAAAAABAw9q4yR2gbt74BAAAAEArOBUWAKA8CyQAAAAAAOgrY+Td4o1PAMBKo1E7O39tvGcAAADou7aOQwBp+RwgJ9omAACAbvKsBwDt4o1PAJlwcjmkJVOwOTkCmnZ5eRmvXr2Kz58/L3x/PB7H2dlZDIfDGI/HcXx8HIPB4NFrAAAAAAAAAABAu9j4BHSWUxkAqrG4HICcTNuXy8vLB9eOjo5m7dV4PI5Xr17F6enpo9dY5JkJAACguxxqVJ0xcgCgLUYj/TwAgK6ZX8Ohr/eQjU8kYeAcoHsshO0Xi8shHQskWMVzUzWHh4dLvz8ejxe+Hg6HcX5+/ug1AKBdjEsAkBPtUvcZIweA/tLX6z7zt/WTIwCgbpU2PukA1ssiOIDqtE2QhsXlkIYFElC/8/Pz2N3dXfje7u5uXF5exqdPn1Ze29vb2+ZtUjOnWUI68gQA0D0WHVZnjBzSMX9bP5/zAOWZvwXYDmvggbqV3vikAwhAbrRNUD+Ly6EaCyS2w6Ruv11fXy/9/tXVVeG1ZW5vb+P29nb29c3Nzaa3B1my4AiA3GibIB15gnoZI0/HoQv9YP4W0tLXg82Zv6UqmzeKaZvqZS0EwGqlNz7pAG5P2xqutt1vLnQAYXPaJqifxeWQhgUSUL9V7VLRtXfv3sXbt2/ruSHIhAVHkJYxve0w5txt2iZIR56gfsbIoRrzt6xy/znPovLH6etBvczfQnXaJkjHfBNUV3rj0yo6gFCdDiDUS9sE9bO4vF+cxLm5lAskIiySoN8Gg8GDfFxdXcVgMCi8tsybN2/i9evXs69vbm7i2bNnye85J11aUNClf0udLDiCdIzpQRraJkhHnqA5xsjXY4N7f5m/her09aBeqedvoQ+0TZCG+SZYz5NN/wepTzi6ublZ+AVddHh4uHSATgcQ0qhjcbn2ib5Kvbj8t99+m/369ddf67ptaI11FkhETBZJfPfdd7NfXd+oETFZFDH9Rb8dHBws/f7+/n7htWV2dnbi6dOnC7+gL4oWHBVdY5G2qV+M6UG9tE2QjjxBOsbIN9eF56a2338OzN9COvp6UK9152+1Te00Pwetz7c+bRNUY76pXj7Tu2vjjU+rrHvCUd8W78E8HUCol8Xl/aUzuz6LyyGNlAskIvq7SKIrtEnVzffVhsPhwrXxeBz7+/sxGAwKrwGLLDiCdIzpQRraJkhHniAdY+SbMQ7GY8zfQnUOKYc0Us/fapvoM20Tnv3SMN+UlnrZPRtvfHLCEaSjA5hOmxusNt97Liwuh81YXA7ppVwgEdHPRRL0z/n5eZycnETEZLLo7Oxsdu309DROTk7i7OwsPnz4sPD69qJrwOMsOILqjOlx/1RYB7CkpW2CdOSpn7RJ1Rkjh/TM36bRlWeNLvwbcuSQcqgm9fyttgke0jZBNQ4zgmIbb3xywhHUTwcQqrG4nHkGWMqxuJzHGLCszgIJ2MzBwUG8f/8+7u7u4v3793F4eDi7NhwOZ997//79Ql6KrgF/sOCIefp69TCmB9Vom9Jq++e6tmkz8gTVGSPnMdqmzZi/hXQcUg7rq3P+VttEn2mbuM/zU1oOM4KJb9f5j66vr2eNTlEn737DZAEfFEvdAXz9+vXs65ubG40WnaZtSmM0ivj4sem7oAkHBwezBeb3TReQR8TCovPHrvEHD/P9cX5+Hr/88ktETAYQnj9/PsvGdBHE8+fP4+LiYukCiWXXAChHX7aag4OD+PDhw4Pv7+/vx3A4XHltlZ2dndjZ2Ul6j9AWxvSY5/lvfdomSEeeoDpj5JCe+VuoR8q+Xp/7ecaT+8P8LdRP2wRp1HGYkTknuqT0xicdQKifDiBUo20CIDcWSACQMwuO6mOjQ78Z04P1aZt4jMWI5ckTwPZMnwG1UQ+Zv4V66OuxjPaoHPO3UA9tU3rmmnCYERQrvfFJBxDqoQPIMgYnytE2AQAAFLPgCOphTA/Wp22qlwUS/SJP9TBHA7A+87eQjr5e/Tw/AVSjbaIM7Ws15ptYxRjlQ6U3PkFZgvY4HUAAAABgGyw4Ss+ETX8Z04M0tE2QjjwBAHSXvh4AudE2QRrmm2A9Nj5lwoKJftEBBHKmTQKKjEY2uAN0gT4fQHXG9AAAANgmczIA5KyLc01d/DcB+THfBOux8QkAAGpiUAw219UceVMuAAAAAADA9thMCAAA7WXjE0l1dVEiAAAAkCdjEQAAAAAAAAAA3WXjEwAAAAAA0AgbWAEAuqfLfbwu/9sA+qCrn+PTf5c3WgG0T1fbJiAdfb0JG58AMjbfqe17gwVAPjxMsS0G+ICqtFEAAN2jj0dTjEsAAAA85FkJgFx1sY3q4r9pXTY+UZvRyCQU2zE/6ekDHihiMyEAObKID4BcGM8DaDfj4wAA3dP1Pp6xCLahDzkCAICue9L0DQAAQBcZYAYAoI1GI31ZAAAAAAAAgJz0fQ7XG58a1vcKCAAAAGyHMQgAAIB+mn8e9FYNAABgXeaaAICmeOMT0BkerIAqfGZQh+np+H2pX335d0Kd5AiAXGmjYHNyBNA+PrsBAACI8HwI0CZ9Wq/XZ9741CABAwAAAAAAYBPmm2BzfcrR9N/qzU8A5ELbBAD0VdfHI0YjfTwgHRufgFbrescPAABgU56bAMiVNgoAgLazYQMAAACgfk+avgEAoHle9QkAzdMWUwf1CgCgm/o4nte3fy9A20zbpj59Xvfp3wpAO2ibAADoKhufGtCnwb4+/VsBAACAPBiLAAAAALbJWAQA0GV9XQfa13839VKnANZj4xNAS+jwAgAAANAFFgwAAAAA22D8AQAAuuHbpm+gbzxMAQB0jz4eAABdNO3nfvzY7H3QLX17fpIjAACgqr49NwEA9E2f+nvGyIFUbHwCgJ7r04MUAEBf6ONBWiZlYHPaJgCA7tLXAwDoLn09IIXRyDwb9dBO9ceTpm+AfvChQkrqE5DKaOQzBVKRJVLy+Qyk4LOEOqhTAJAHbTIAAACwLcYhSMHcJZDC9LOkj58n3vgEAACsrY8PUdRLnYLNyREAAECx+ecmpw0DkIu+j+t5CwAAQLdN+7v6fMA6bHzaAg/mQCo6fqTU94FzAMiR/h6kZ1wGAOiyvo3x9e3fC3UzDsEmfCY/JFOwOTmCzc3nSKYAgC6yBqCfbHzaEoN+kJZMGZwAAAAAgLYwCUdKxscBAABW8ww+4dkRAOgyfZ3+sfEJaA2NFKQjT4sMfEIaNuUCAAAAAABUZ/4W0pKpReZxqUqGAIDc2PjE1niAAgDoB5sJAciNMQk2YYIXgFxokxbp45GKbAEAAABA/ozj9duTpm+A/vGhAwB0xWikbwMpyNFDyoSq1BkAyJ9nSAAAgHI8OwGQK20UbE6OgHXY+ATQQjp+AEDXWAQKAHnTVrMJdQcAAAAAgD4xLg6Q1rdN3wBAGTqBkIYsQTryBNRt+jnz8WOz9wEA80YjbRNsSo4AyMH98U3jEJCWTAFsl7lboG4+Z1iH5wKAdLzxiUY4IRaAXGmfAIA28XxdzbS8lBmrqB/AtvnMgc3JEWxOjqjCcxMAAAAA2+aNTwAAAAAAc7yNZZFFjVQxf4KlurOasqEsp8ICAAAAQLsZDwbq0Lf5AxufaqKRgs1ZZFTO/OeN8mIZbRJsrm8PCbAt2qhyfAYBNEM7BQBAXxmL4D51Aupj/AHYJmuxIC2ZAoDtsfEJyJpBPqBJNhZyn0ErAJqmf5KOdh3WZ9EjwHbM9/2MlRfTNlFEfqpRXrA5OYI0jIUCbJ9+TDnLysm8E1RjPA/S6UsbZONTYjp+1Wm8AMiBhSRQj748WFGez1ioj3zBemQH6iFb5RkjBwDoNv097vO8tD6boQDS0y5BOtbIwObu58gzNVNPmr4BANLxIApQv9HI520ZygnSkynuUyc2p72CzcgPUz5PIR1ZgrRkCtKSqeqUGUB68+MQxiSqU2ZEqAdQJ/mqTplRhjpSrA858sanBLpeSQBoFydHQBpOiwBy5LMJAMiZsXKgafffAOAZCmBz+ngAAFCOcQiozjMnRdSP6rq8ftjGpzXpoKQx/4HU5aBRjYYqDZkCIBfzfWf9aAAAgMcZ24PNmWvoNz//NLRH/SVDkJ75kc34XAIAAOi3rWx8Go/HcXZ2FsPhMMbjcRwfH8dgMNjGXw2d0ocsGazanDIsp+t5Ug/Ypi7mSYbqcX/Te4QJvnldzBI0pWt5stCsfsp4ua5liXro15UjT5TlebRYV7Pk504Tupon2DZZoor7byZkkTxRlTG91eQJ0uhqlpbN27O+ZWV4/3vaq+7mCbatq1m6/7mpfWKZrWx8Ojo6is+fP0fEJHCvXr2K09PTbfzVtROstCyUKCZLVCVTq3UtTzJUv1VlLF/dy9N98sW2dD1LsE3yxDq8ofChrmZJ/64+FvGt1rU8+VnTlK5laRVtFdvQtTzNL3yWobSWtfsWmv+hi1li+2Rqoit5eixHcpaWcbzlupKnebJTP3l6qItZgqZ0OU/aqHRWje9om/7Q5SyRXtc24ta+8Wk8Hi98PRwO4/z8vO6/lpbTYD0kS2xCphZ1KU8enGhal/JEcwyidzNL9/sf2qy07i+GkKM/dClPTtxrXp8XHnUpSxH9/llui8+p1bqWp/v87NMpOtHQZ1g3syQ/NKWLeWI7fG4t6kKWPCs1R9kv6kKeIBd9ypO+yWaUX7GuZcnPmyZ1LU/Uy+fVarJEWV3N0ZO6/4Lz8/PY3d1d+N7u7m5cXl7W/VcnNxr98Qu2rUtZiljMk0xtX9/Lve156vvPj7y0MU/ykxeL+ifamKUq+vyz3YZVOeprubcpT56J8rRs0Xkff0ZtytIqfpbkogt5mpKjvPTtZ9HmLMlO85T/ojbnKWIxU6t+T/2WlXffyr+tWVr1rNS3n18OHvsM69PPpW15kqP83F+D0uefRdvydJ95jub5GUy0PUuQky7mqc+fj01Q3hNdyNKyfruf7/a1tcxrf+PT9fX10u9fXV09+N7t7W3c3t7Ovv7tt98iIuLm5qaWe1vmf//3j9///PPi1zTjr3+d/CwiJj+Px36/DdM6eXd3t7W/s0qWIprP06rsTH9O//rXVm6DJf761+W/32aGpprIUkQ78nT/M27q55//yM/8z4/m3P853P+5bStbbcjTtrO0qn/wr38t//zT78uDtmlRTm1TEflpxqq+wLQa/Otff/z+fn9iG9qQpyazdD83+nb5u7mRpYjm26bH+tnzfcD7/T62a1XZN9HHm2pDnraRpSqfZavGJ9i+x/p+29SGLEU039ebz5fxvDwsK/8m26UIeSpijLw97v8cGhqeMn/7CM9K7bHsZzNfRepeG6FtemjZWMSqsVfyYs5pUW5j5BHLxyjut1ParGZYW5R32zTvsXbIWr08lMnRNtYZtSFPTbdNU2XWv2qjtmtV/1vblF/bFGF+KVfL+nhl5t3rVCVPtW98WmVZ+N69exdv37598P1nz55t4Y4e+u67Rv5alpj/WZT5/bb8/vvv8V3DFWVVQ5ZbnqbkKl9N/mxyyFJEfnlaViQZFBOPuP8z2vbPLOc8NZGlMkWRQXGxgrYpv7aJdnjsOUnb9AdZogr9vIlc2qaiosigmCiQw88n5zzllKUqf4bmeG7Kp226L4OioYRcfk7ytFwOz7RU1/TPKIc85ZalqaZ/NqyviTGJHLIUkVeemh4bYjOenfIYhyiSQRFRgizl1TaVkUGRcc9jP5O+9/VyzdJUBsXWW6vKXtuUb9uUQdFQQi7z7mXyVPvGp8Fg8GAn4dXVVQwGgwd/9s2bN/H69evZ1//+97/j6uoq/uM//iO++eabB3/+5uYmnj17Fr/++ms8ffo0+b13jfIq77Gyuru7i99//z3+8pe/bO2eqmQpQp7qpKyqKSqvJrIUUW+e1I9qlFd5ObZNEfp6uVBW1WibKKK8ytM2qR9FlFU12iaKKK/ytE3qRxFlVU3f2qYIdaQKZVVN3/KkflSjvMrLsa+nbcqHsqpG20QR5VVejm1ThHGIXCirarRNFFFe5Wmb1I8iyqoabRNFlFd5Kdum2jc+HRwcxIcPHx58f39//8H3dnZ2YmdnZ+F7q8I47+nTpypNBcqrvKKy2vYu3SpZipCnbVBW1awqryZ2vG8jT+pHNcqrvJzapgh9vdwoq2q0TRRRXuVpmyiirKrRNlFEeZWnbaKIsqqmb21ThDpShbKqpm95Uj+qUV7l5dTX0zblR1lVo22iiPIqL6e2KcI4RG6UVTXaJooor/K0TRRRVtVomyiivMpL0TY9SXlDywyHw4Wvx+Nx7O/vl2qIgD/IEqQjT5COPEEasgTpyBOkIUuQjjxBGrIE6cgTpCFLkI48QTryBGnIEqQjT5CGLNF3tb/xKSLi9PQ0Tk5O4vnz53FxcRGnp6fb+Guhc2QJ0pEnSEeeIA1ZgnTkCdKQJUhHniANWYJ05AnSkCVIR54gHXmCNGQJ0pEnSEOW6LOtbHwaDofx/v37iIg4PDxM9v/d2dmJ//u//3vwGjaWU17l5VpWdWUpIt9/c46UVTW5lpe2KQ/Kq7ycy0qemqesqsm1vGQpD8qrvJzLSp6ap6yqybW8ZCkPyqu8nMtKnpqnrKrJtbyMkedBWVWTa3lpm/KgvMrLtay0TXlQVtXkWl7apjwor/JyLit5ap6yqibX8pKlPCiv8nIuK3lqnrKqJtfykqU8KK/yUpbVN3d3d3cJ7gkAAAAAAAAAAAAAAAAgmSdN3wAAAAAAAAAAAAAAAADAfTY+AQAAAAAAAAAAAAAAANn5tukbSOXy8jIiIvb29mI8Hsf19XXs7e01fFd5GI/HcXZ2FsPhMMbjcRwfH8dgMGj6trKlLimDIvJUnnqkDIrIUjXqkjIoIk/lqUfKoIgsVaMuKYMi8lSeejShHJaTpWrUI2VQRJ7KU48mlMNq8lSeeqQMishSNeqSMigiT+WpR8qgiCxVoy4pgyLyVJ56pAyKyFI16pIyKCJP5alHE8phOVmqJnk9uuuI4+Pju4i4i4i7g4ODu69fvzZ9S9nY29ub/f7Lly93h4eHDd5N/tQlZVBEnspTj5RBEVmqRl1SBkXkqTz1SBkUkaVq1CVlUESeylOPJpTDcrJUjXqkDIrIU3nq0YRyWE2eylOPlEERWapGXVIGReSpPPVIGRSRpWrUJWVQRJ7KU4+UQRFZqkZdUgZF5Kk89WhCOSwnS9WkrkedeePTf//3f8fXr18jIuycmzMejxe+Hg6HcX5+3tDdtIO6pAxWkadq1CNlsIosVacuKYNV5Kka9UgZrCJL1alLymAVeapGPZpQDg/JUnXqkTJYRZ6qUY8mlMNy8lSNeqQMVpGl6tQlZbCKPFWjHimDVWSpOnVJGawiT9WoR8pgFVmqTl1SBqvIUzXq0YRyeEiWqktdj55s/H/IyGAwEK57zs/PY3d3d+F7u7u7s1eHsZy6pAyWkafq1CNlsIwsrUddUgbLyFN16pEyWEaW1qMuKYNl5Kk69WhCOSySpfWoR8pgGXmqTj2aUA4PyVN16pEyWEaW1qMuKYNl5Kk69UgZLCNL61GXlMEy8lSdeqQMlpGl9ahLymAZeapOPZpQDotkaT0p61Fn3vh0fX0dZ2dnERFxcXER33//fQyHw4bvqnnX19dLv391dbXdG2kRdUkZrCJP1ahHymAVWapOXVIGq8hTNeqRMlhFlqpTl5TBKvJUjXo0oRwekqXq1CNlsIo8VaMeTSiH5eSpGvVIGawiS9WpS8pgFXmqRj1SBqvIUnXqkjJYRZ6qUY+UwSqyVJ26pAxWkadq1KMJ5fCQLFWXuh51ZuPT8fHxbDfYcDiMFy9exJcvX5q9qYytCh/qUoQyqEqellOPlEFVsrSauqQMqpKn5dQjZVCVLK2mLimDquRpOfVoQjmUJ0urqUfKoCp5Wk49mlAO1cjTcuqRMqhKllZTl5RBVfK0nHqkDKqSpdXUJWVQlTwtpx4pg6pkaTV1SRlUJU/LqUcTyqE8WVotdT3KeuPTjz/+WPiPe/HiRRwcHERExHg8jr29vYiYFMx4PI7xeNz73YWDweDBTsKrqyuvnivQ1bokT5uTp2q6Wo9kaXOyVF1X65I8bU6equlqPZKlzclSdV2tS/K0OXmqpsv1SJ42I0vVdbUeydLm5KmaLtcjedqcPFXT1XokS5uTpeq6WpfkaXPyVE1X65EsbU6WqutqXZKnzclTNV2tR7K0OVmqrqt1SZ42J0/VdLkeydNmZKm65PXorgM+f/58NxgMZl9//fr1LiLuvn792txNZeLLly93e3t7C98bDAbKZgV1SRkUkafy1CNlUESWqlGXlEEReSpPPVIGRWSpGnVJGRSRp/LUownlsJwsVaMeKYMi8lSeejShHFaTp/LUI2VQRJaqUZeUQRF5Kk89UgZFZKkadUkZFJGn8tQjZVBElqpRl5RBEXkqTz2aUA7LyVI1ddSjJ+ttl8rLcDiM9+/fz74+Pz+Pw8NDO+giHuyIG4/Hsb+/r2xWUJeUQRF5Kk89UgZFZKkadUkZFJGn8tQjZVBElqpRl5RBEXkqTz2aUA7LyVI16pEyKCJP5alHE8phNXkqTz1SBkVkqRp1SRkUkafy1CNlUESWqlGXlEEReSpPPVIGRWSpGnVJGRSRp/LUownlsJwsVVNHPfrm7u7uLsG9Ne7y8jLOz89jMBjEly9fFgqq78bjcXz48CGeP38eFxcX8ebNGyEroC4pgyLyVJ56pAyKyFI16pIyKCJP5alHyqCILFWjLimDIvJUnno0oRyWk6Vq1CNlUESeylOPJpTDavJUnnqkDIrIUjXqkjIoIk/lqUfKoIgsVaMuKYMi8lSeeqQMishSNeqSMigiT+WpRxPKYTlZqiZ1PerMxicAAAAAAAAAAAAAAACgO540fQMAAAAAAAAAAAAAAAAA99n4BAAAAAAAAAAAAAAAAGTHxicAAAAAAAAAAAAAAAAgOzY+AQAAAAAAAAAAAAAAANmx8QkAAAAAAAAAAAAAAADIjo1PAAAAAAAAAAAAAAAAQHZsfAIAAAAAAAAAAAAAAACyY+MTAAAAAAAAAAAAAAAAkB0bnwAAAAAAAAAAAAAAAIDs2PgEAAAAAAAAAAAAAAAAZMfGJwAAAAAAAAAAAAAAACA7Nj4BAAAAAAAAAAAAAAAA2fl/FlhhhtfpPc0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 3400x200 with 16 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAD5YAAADvCAYAAACD1NktAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABNlElEQVR4nO3dv24bSZso7neMARitzeHJFsdJhyfkR13B0MAmysTPVyA5VzCEogNHgnwHtq/AIjOH5lyBbN6BOzHOhrTMLyIWWP8C/8QVRZGiJHaz/zwPIIzIamtapXrZ3VX1Vv328+fPnwEAAAAAAAAAAAAAAAAAAEBlPdn1CQAAAAAAAAAAAAAAAAAAAJAtieUAAAAAAAAAAAAAAAAAAAAVJ7EcAAAAAAAAAAAAAAAAAACg4iSWAwAAAAAAAAAAAAAAAAAAVJzEcgAAAAAAAAAAAAAAAAAAgIqTWA4AAAAAAAAAAAAAAAAAAFBxEssBAAAAAAAAAAAAAAAAAAAqTmI5AAAAAAAAAAAAAAAAAABAxUksBwAAAAAAAAAAAAAAAAAAqLjfd30CAAAAAABQB/1+P05OTqLZbEZERJqmMRwOI0mSSNM0jo6ONioDAAAAAAAAAIBdMg8Gyuu3nz9//tz1SQAAAAAAQJWNx+P4xz/+Ed+/f58PjP3jH/+IL1++RMSvAbR+vx+DweDOMgAAAAAAAAAA2BXzYKDcCr1j+X//93/Hf/7nf8a//du/xW+//bbr04FC+/nzZ/zrX/+Kf//3f48nT548+ufltWqMOIfNbDvG8yLGYTNiHKpNjEO1iXGotm3GeJqmkSTJwuvrkiSJ0Wh0Z9kmxDhsxnUcqk2MQ/WVMc7FOGxOjEO1iXGoNjEO1SbGodrEOFRbGefBiHHY3H1ivNCJ5f/5n/8Zz58/3/VpQKl8+/Yt/vf//t+P+hnj8TjevHkTJycn8/d6vd7CyjCHh4fzlWHWld1FnMP9bCPG8yTG4X7EOFSbGIdqE+NQbY+N8eFwGAcHB9Hv9+fvjUajaLVaC8e1Wq0Yj8fx+fPnlWXtdnvp589ms5jNZvPX/+///b/4P//n/zz4fKFuXMeh2sQ4VN+24jyPxdfFONxfma7lYhzuT4xDtYlxqDYxDtUmxqHaijwPxhwYeLxNYrzQieX/9m//FhG/fpGnT5/u+Gyg2KbTaTx//nweN4+R5+5J4hw2s80Yz5MYh82Icag2MQ7VJsah2rYR45eXl7cmmFxeXt56/GQyWVt2m9PT03j9+vXS+2Ic1nMdh2oT41B924zzvBZfF+OwuTJey8U4bE6MQ7WJcag2MQ7VJsah2sowD8YcGHi4+8R4oRPLf/vtt4iIePr0qcCHDV3FzUPlvXvSv/71r4gQ57Cpx8Z43lzL4X7EOFSbGIdqE+NQbY+J8fPz8zg6Otr4+FWDaevKTk5O4vj4eP76apBAjMNmXMeh2sQ4VN824jyvxdfFONxfma7lYhzuT4xDtYlxqDYxDtUmxqHaijwPxhwYeLxNYvxJDucBlEReuyc9e/Zs/vX8+fMHni0AAAAAFNtoNIp//vOft5Y1m82lPrTJZBLNZnNt2W0ajcZ8AM1AGgAA3M/V4uvXrVt8fV3ZTbPZLKbT6cIXAAAAAACUUR7zYMyBgXxILAfmzs/Po9vtbnz8Q3dP+vHjx/zr27dv9zxLAAAAACiP8/PzePfuXbx79y7SNI3T09MYj8cr++E6nc7aMgAAYHuyXnzdwusAAAAAAFSJeTBQDb/v+gSAYshz96RGo7GVcwYAAACAIrs5MPbq1at49epVJEmydGyaptHpdOZ9bqvKAACA7Tk/P4+jo6ONj7/v4usnJydxfHw8fz2dTiWXAwAAAABQSubBQHVILAfmzs/P599frRrz8uXL6Ha78fbt26XjO51OJEmysgwAAAAA+JVg8u7du4iIODs7i1evXkW73Y7BYBD9fj/29vbi4uIiBoPB/N+sKwMAAB4vj8XXLbwOAAAAAEDVmAcD5SexHIgIq8YAAAAAQFaazWb89ddf8ddffy28nyRJnJ2dRUTEwcHBxmUAAMB2WHwdAAAAAADuxzwYKD+J5cACq8YAAAAAAAAAUHUWXwcAAAAAAKCOJJYDC6waAwAAAAAAAEBdWHwdAAAAAACAOpFYzoPs7//678ePuz0P4PE2iWcxD+V1n/gV61Ae4rV+RqNRRPya5HpxcREvX76MdrsdEb92RBoOh5EkSaRpGkdHR/PdkdaVUR5iHupL/EN1XI9nsQ31IuahmsR2tVh8vd7EM1TDqlgW41BNYhuI+PVZ4HMAykXcQvW5V4d68xlQLhLLAQAAoOB6vV78/fff0e12YzKZRK/Xi69fv87Lvnz5EhG/EskPDw/nOyStK6N8dLoBQPlcXb+BanBPDgDl53oO9eYzAAAAIF/9fj9OTk422ijFJioAkB+J5QAAAFBwg8FgvkN5RCx0pl+XJMl8d/N1ZQAAAGTHrisAAAAAkD2LxgAU23g8jjdv3sTJycn8PZuoAEAxSCwHAACAgut2u/PvB4NBvHr1KiIiRqNRtFqthWNbrVaMx+P4/PnzyrLrSepXZrNZzGaz+evpdLrNX4F7MPgNAADld3VfDwAAAADAoquNES4vL+Pi4iJevnw5n8tiN2OojjRNI0mShdfX2UQFAHZHYjlLTGAHACgfne1QfePxOD58+BAvXryIo6OjiPgV87eZTCZry25zenoar1+/3sapAgDUnoRSAAAAAACA2/V6vfj777+j2+3GZDKJXq8XX79+nZfZzRjKbzgcxsHBQfT7/fl729pExQYqAPB4EssBAKACdLazKQtJlVe73Y4kSaLf78873ldZlVS+ruzk5CSOj4/nr6fTaTx//vyhpwvACq7FAFB+rucAAAAA2bLBAlTbYDBYSBS9HsPX2c0Yyuny8vLW6++2NlGxgQoAPJ7EcgAAqACd7VAPzWYzer1evHjxIr5//x7NZnOp83wymUSz2VxbdptGoxGNRiOrUwfgDhLUAAAAAAAAfrHBAlRbt9udfz8YDOLVq1cRYTdjqIrz8/M4Ojra+Pj7bqJiAxUAeLwnuz4Bym1//38mvQL1IfYBiuchne3rym4zm81iOp0ufFEsrs/VNBqN4o8//pi/TpIkIn4Ngl+P/es6nc7aMsrBfTcAAJSDe3cAAACA7bLBAlTfeDyOfr8fL168mCegbnM342fPns2/JJwWx1396frby280GsU///nPW8u2tYlKo9GIp0+fLnxRLGIZoPgklrOSCznUj7gHKLcsO9sjdLjDrrRarYUk8fF4HM1mM9rt9jzJ/EqaptHpdKLZbK4tozrcwwMAAAAAAABVk/UGCzZXgN1rt9txcnISX79+jeFwuPbYh+xm/OPHj/nXt2/fHnm2wH2cn5/Hu3fv4t27d5GmaZyensZ4PLaJCgAUyO+7PgGq6WpS+8ePuz0PYHOSUQDK7yrJtN/vx3A4jIODg5XH3rezPeJXh/vx8fH89XQ6lVwOOWi32/Hy5ct49+5dRER8+vQpvnz5Mi8fDAbR7/djb28vLi4uYjAYbFRGubl/r5arnRIuLy/j4uIiXr58Od+BIU3TGA6HkSRJpGkaR0dHCzsyrCqjPG7G8/6+PjUAACgCY95QL2IeAKB4xuNxfPjwIbPdjF+/fr2tUyUDN+/R3bNXU7PZjF6vFy9evIjv379vdTfjRqOR5akDK9xMEH/16lW8evVqaZOUiMWNUm7Gsk1UACBbEsu5FxPXAQCKLavO9ggd7rBL1xeKuBowv5IkSZydnS0dd1cZ5eXZvHp6vV78/fff0e12YzKZRK/Xi69fv87LrhaTSNM0Dg8P54tErCsDACB797k3N/EVACB/FnTkin51gHLKcoMFmyvA7oxGo+j1evH9+/eIiHmyaZqm0e124+3bt0v/ptPpRJIkK8uA4rm8vJxvpHJ2dhavXr2KdrttE5WaMk4GUDwSy7mTCzhANRg0Z1Ou/eWjsx2AK67j5TQYDOb35hGxcC9+XZIk8/v6dWUAAOyGZBWoNs/cAOVjQUcAKD+7GUP1tFqthR2Nx+NxNJvNhTHzK3YzLrf79KfpX6+eZrMZf/31V/z1118L79tEBQCKQWI5ANSEQXOoLp3t1XdXJ7uOdYByu34dHwwG8erVq4j4tXhMq9VaOLbVasV4PI7Pnz+vLLvtHmA2m8VsNpu/nk6n2/wVyJgEFgAAAID7s6AjVJsNFqC6bLBQT8bD6qPdbsfLly/nOxl/+vRpPkc1Yv2OxXYzLidz26AexDpAuUgsB2DJqg46N/vlZtAcqktnOwCU33g8jg8fPsSLFy/i6OgoIn5NhrvNZDJZW3ab09PTeP369TZOlXsyCQag2kxiBwCA4sp6QUeLOcJu2WChflbNXdMPXz02WIDqu74T8dX4+BW7GQMAQPYklgNATdgFkZsMrFWLznZuEuNQbz4DyqfdbkeSJNHv92M4HK69Lq9KKl9XdnJyEsfHx/PX0+k0nj9//tDTBQD+fyaxAwBAsWW5oKPFHMtrf1//eRXYYAGqywYL9WYDJCi/+8TxuoVj3LMDAGTnXonldl1gFRPWAcrBLogAAFBczWYzer1evHjxIr5//x7NZnPp3nsymcx3XFhVdptGoxGNRiOrUwdgQybEVY9J7AAAUGxZLuhoMUfYraw3WLC5AuyWDRYAAAAgO/dKLLfrAkC9mOhaPXZBBKgf13OAYhuNRtHr9eL79+8R8WuyS8SvPrRutxtv375d+jedTieSJFlZBkD+HnLfbcHWash6EnuEiewAAPBYWS3oaDFH2L0sN1iwuQIAAAAAVfXkPgfbdQEAyu9q0LzX68Xl5eXWd0F8+vTpwhcAALtjcYnia7VaCwlp4/E4ms3mfFGo69I0jU6nE81mc20ZAJCv8Xgc/X4/k0nsEb8msj979mz+ZSFHgGzt73ueBqiC0WgUf/zxx/z1zQUdb9PpdNaWUU2u/eXWbrfj5OQkvn79GsPhcO2x991g4eTkJH78+DH/+vbt2yPPFgAAAACK4V47lme964IdFwAgG3ZB5DqD4lBtYhygXNrtdrx8+TLevXsXERGfPn2KL1++zMsHg0H0+/3Y29uLi4uLGAwGG5UBAPm5WhCm3+/HcDiMg4ODlcfedxJ7xK+J7MfHx/PX0+lUcjlAjvS3AZTTugUdb7q+aOPNhRst6AjFd7XBwosXL+L79+9b22Ch0WhEo9HI8tQBAAAqTf86QHHdK7E84lcn+4cPHzLZdeH09DRev35931MiJy7oAOVl0BwAAIrrevLZVX/blSRJ4uzsbOm4u8oAgHxlNYk9wkT2IjBGBgD1s78f8fHjrs+Cx7CgI1SbDRYA6ktfHQAAwOPdO7E8y10X7LgAANkwaM4VHesAAABQLFfP6pJWysckdgCAcrgtSXzVmJmxtGqxoGN9bRLL4r3cbLDAJiwUA+XhugzcxhgaAEB27p1YHpHdrgt2XKgeHXMAxWHQHAAAAAC2xyR2AAAA2A0bLBAhERWqQBwDQL1YMAKgODZOLLfrAgBAeeh0B4B6ce2HahHTAOVgEjsAQHl41gaoHhss1INrOAAgERUAisMmxNWxcWK5XRcAAACKxSA6VIuYBgAoH5PYAerH8zsAAAAAAPyizxygnDZOLLfrAgAAAAAAAAAAUAd2RAOAcrnt2u16DgDFJBEVAGC3Nk4sj7DrAgBAWemEAwAAAAAAAAAAAAAAgHq7V2I5AAAAAAAAANTBaDSKiIjLy8u4uLiIly9fRrvdjoiINE1jOBxGkiSRpmkcHR1Fs9m8swwAAAAAAAAAdkliOVtlN1QA2C3XYgAAKI5t359f/byPH7f7cwGA/Liel0uv14u///47ut1uTCaT6PV68fXr13nZly9fIuJXIvnh4WEMBoM7ywAAAAAAAACqQg5LOUksB6gZF2wAAKgXiSsAAAAPMxgM5juUR8TCjuTXJUky3918XRkAANtnHgwAlI/rNwAAwG5JLAcAAAAAAACAG7rd7vz7wWAQr169ioiI0WgUrVZr4dhWqxXj8Tg+f/68sux6kvqV2WwWs9ls/no6nW7zVwAAAACovJuJ6hZdBwAAWO/Jrk8AAAAAAACAX/b37dYCUCTj8Tj6/X68ePEijo6OIiLi8vLy1mMnk8nastucnp7Gs2fP5l/Pnz/fxmmTM9dvAAAAAAAAoCzsWA4AwEpXE+Gs4grFYHIqAABU1837fc/kAMXQbrcjSZLo9/sxHA7j4OBg5bGrksrXlZ2cnMTx8fH89XQ6lVwOAADwCMbVAQAAIHuev8tNYjkAAEANSVIBAAAA2Eyz2YxerxcvXryI79+/R7PZXNqBfDKZRLPZXFt2m0ajEY1GI6tTBwAAACg1ySoAUF37++awQhm5R6+GJ7s+AQAAAAAAAIA62d834F50o9Eo/vjjj/nrJEkiIiJN0+h2u7f+m06ns7aM6hHLUA9XsS7modrEOAAAlIN7dwCAx7NjOQBAzelgAwCA/LkPBwAotlartZAkPh6Po9lsRrvdXjo2TdPodDrzHctXlQEAAAAAAAAU1dWcto8fd3seZE9iOQAAAAAAZGQ0GkVExOXlZVxcXMTLly/nCWlpmsZwOIwkSSJN0zg6Oponna0rAwCy12634+XLl/Hu3buIiPj06VN8+fJlXj4YDKLf78fe3l5cXFzEYDDYqAwAgOxYzBFYxcR4AMhe1vfj7vehWMRkfZkHUw+eo6tPYjm58YECANnxcA4AAFBMvV4v/v777+h2uzGZTKLX68XXr1/nZVcJammaxuHh4TzxbF0ZAJCPg4OD+fdHR0cLZUmSxNnZ2dJxd5VRTPv7xrEBAAAAACAP5r2Xm3kwUA0Sy4E5q8YAcJ2HdgAAgMcbDAbzPraIWOhTuy5Jknn/3LoyAAAAAAAAAADYBfNgoBoklgNzVo0BAADIn8VcAKqt2+3Ovx8MBvHq1auI+LXIY6vVWji21WrFeDyOz58/ryy7Pjh3ZTabxWw2m7+eTqfb/BUAAKCyLL4OUF761gEAAADyl/U8GHNgIB8Sy5nT2Y5VYwAAisd9OgBA+Y3H4/jw4UO8ePEijo6OIuJX4sptJpPJ2rLbnJ6exuvXr7dxqgAAtaUfrp4svg4AAFA9V8/4Hz/u9jwAAKoqy3kw5sBAPp7s+gSA4njIqjHrym4zm81iOp0ufAEAAACwuf19SS9l02634+TkJL5+/RrD4XDtsasG09aVnZycxI8fP+Zf3759e8TZAgBAfVh8HaB89I0BAAAA7FaW82DMgYF82LEcWGD3JAAAAADYvmazGb1eL168eBHfv3+PZrO51Ic2mUyi2WyuLbtNo9GIRqOR1akDAEBlPWTx9c+fP68su56kHvFr4fXZbDZ/beF1AACgziwOAwBQHVnNgzEHpvyu7vs/ftztebCeHcuBBXZP4j508sHu7Wo1dqvAAwAA3G00GsUff/wxf50kSUT82uHwegLLdZ1OZ20ZANWkrw1gd8bjcfT7/UwWXz89PY1nz57Nv54/f76t02YDxrMAAAAAALbHPBioDjuWA0vsngQAsDt5r9K2v29FOACArLRarYXBsfF4HM1mc2kXw4hfg2ydTmfe57aqDAAA2K52ux1JkkS/34/hcBgHBwcrj73v4usnJydxfHw8fz2dTiWXA5SEnZUAAGAzFnUDqA/zYKA6JJYDEfFr1Zherxffv3+PiOVVY96+fbv0bzqdTiRJsrKMYvCwDgAAALAb7XY7Xr58Ge/evYuIiE+fPsWXL1/m5YPBIPr9fuzt7cXFxUUMBoONygAAgO3KavF1C68DAAAAAFAV5sHUl9y06pFYDkSEVWMAAAAAsqRzvb6u73Z4dHS0UJYkSZydnS0dd1cZAADweBZfB2ATdi6HfIk5AIDyG41GERFxeXkZFxcX8fLly3leSpqmMRwOI0mSSNM0jo6O5rkn68qA4jAPBqpBYjmZM2m2HKwaAwAAUG0m4gDAdun7BoDycz2vN4uvAwAAAMD29Xq9+Pvvv6Pb7cZkMolerxdfv36dl13lqaRpGoeHh/P8k3VlAMB2SSwH5qwaAwBQT5JNAWA7JKUAQHV4VoZqc+9OhMXXqyyvGHe/AAAAALBsMBgsLN54fUfy65Ikme9uvq4MgPLSj15cEssBAAAAAAAAgMqTUM5NFl8HKA/XcQAAgHLodrvz7weDQbx69SoiIkajUbRarYVjW61WjMfj+Pz588qy60nqAMB2SCwHAAAAAAAAiOVkFckrAAAAAAAA9zMej+PDhw/x4sWL+YKOl5eXtx47mUzWlt00m81iNpvNX0+n00efLwDbZ6y92J7s+gQAAAAA6mZ/X6cZAAAAAJSRvj0AKC/XcQCoDtf1Ymu323FychJfv36N4XC49thVSeWryk5PT+PZs2fzr+fPnz/ybLkPsQdQDRLLAQAAAAAAAAAeyYQ6AAC4H/fQUG1iHKDems1m9Hq96PV6cXl5Gc1mc2kH8slkEs1mc23ZTScnJ/Hjx4/517dv37L8NQCgkiSWAwAAAACUkMk4AAAAAABUkf5vAIByGo1G8ccff8xfJ0kSERFpmka3273133Q6nbVlNzUajXj69OnCFwDZ8HxeXb/v+gQAAAAA6kqHGwAAAAAAFId+ewDYPtdXgPpotVoLSeLj8TiazWa02+2lY9M0jU6nM9+xfFUZALB9EsvxsA4AAAAAAAAAAAAAAAA8WLvdjpcvX8a7d+8iIuLTp0/x5cuXeflgMIh+vx97e3txcXERg8FgozIAYLsklgPwKFeLU3z8uNvzALIl1iF/FoACAACA6vP8DwAAAAAAVMnBwcH8+6Ojo4WyJEni7Oxs6bi7yqguc9QBdkNiOQBASZl0CgAAAPkxoA0AbMp9AwAAlIN7dwAAAJCbUkdPdn0C7IZgBwAAAAAAAAAAAAAAAACA+pBYXmP7+xLMAQAAAAAAAAAAAACoNvPmYXvEEwBAuf2+6xOgfq4eIj5+3O15AAAAAAAAAAAAAFAvkuEAAACAOpNYDlBhOsABAKphPB7HaDSKiIiLi4t4//59NJvNiIhI0zSGw2EkSRJpmsbR0dFGZQBUx/6+RRyhysQ4AAAAAAAAAACwLRLLAQAAoOBGo1H89ddfERHx5s2b+PPPP+PLly8REdHr9ebfp2kah4eHMRgM7iwDAACgWK4Wi7WYBAA8jIXXAQAAAAAA7iaxHAAAYIdMGucu4/E4Tk9P54nlBwcH0e/3I03TpWOTJJnvbH6z/HoZ+SlijBfxnIDHE9sAAAAAAJSVBWIA6mU8Hs/nsFxcXMT79++j2WxGxK/5LsPhMJIkiTRN4+joaKMyAABgcxLLAQCgAnS2Q3W12+14//79/PXl5WVERLRarTg/P49Wq7VwfKvVivF4HJ8/f15Z1m63l/4/s9ksZrPZ/PV0Ot3ibwEAAAAAAFBe+/sWdQTIk3kwUG2j0Wi+wcKbN2/izz//jC9fvkRERK/Xm3+fpmkcHh7GYDC4swwAANicxHIAAKgAne1QbQcHB/PvP3z4EN1uN5rN5jzJ/KbJZLK27Danp6fx+vXrx54qAEDl2T0JALgviWgAAAD3Yx4MVNd4PI7T09N5jB8cHES/3480TZeOTZJkvsjEzfLrZQCU39U4vL50gHw82fUJAABwP/v7u5vEvsv/N6tddbZfOTg4iPF4HGmaru1Q19nOKmK9uC4vL2M4HN458L0qqXxd2cnJSfz48WP+9e3bt0ecKQAAAAAAAADcn3kwUG3tdjvev38/f301j6XVasVoNIpWq7VwfKvVivF4vLbsptlsFtPpdOGL6jLXDQDg/uxYDsBWWSkKIH/rOtvPz89Xdqh//vx5ZVm73V76/8xms5jNZvPXOtwhf/1+Pz59+hTNZjMiIprN5tIO5JPJJJrN5tqy2zQajWg0GlmcNnCHq0HwiIiLi4t4//79PFbTNI3hcBhJkkSapnF0dLRRGfXluRwAAIC62fXk8f19z+EAANuUxzwYc2Bgtw4ODubff/jwIbrdbjSbzZWbJUwmk7VlN52ensbr16+3caqssevncQAAHk5iOQDUhGQVqLasO9sjdLhnTUc7d3nz5k30+/1IkmQev91uN96+fbt0bKfTiSRJVpYBxTIajeKvv/6KiF+x/ueff8aXL18iIqLX682/T9M0Dg8PYzAY3FkGAMBmLMoCAFBdxsgBoLwknbIN+v6K7/LyMobD4Xzce91x9yk7OTmJ4+Pj+evpdBrPnz9/6GkCAEDlSCwHYCskwxWfZBWoh6w62yN0uMMuDYfDaLfb86Ty8/PzWyeypWkanU5nvmP5qjLyZ8CaVcbjcZyens7v1Q8ODqLf70eapkvHJkkynwh7s/x6GQCQH4kqAABQXMbIodo8k0M9SDqFauv3+/Hp06f5tbjZbC4tBjGZTObzYFaV3dRoNKLRaGR12gC1lFe+iHl2APmQWF4zEj8B6kmyCtRHVp3tETrcYVfSNI1er7fwXrPZjKOjo4iIGAwG0e/3Y29vLy4uLhYmt60rA4qh3W7H+/fv56+vJra0Wq04Pz+PVqu1cHyr1YrxeByfP39eWdZut5f+P7PZLGaz2fz1dDrd4m8BAPUlUQUAAIrJGDlUn2fycjOXlU1JOoXqevPmTfT7/fkmCxER3W433r59u3Rsp9OJJElWlgEAUC43+wUsJpG/J/c5eDwex5s3b+LNmzfR6/UWVnBL0zTevHkTw+Ew3rx5s3EZ+djf1xEHUGfrklVGo9HKhJR1ZbeZzWYxnU4XvoD83Oxsv7y8jG63e+uxnU5nbRlQHEmSxM+fPxe+vn//vlB+dnYWBwcHcXZ2tjAovq4MKI6Dg4P59x8+fIhutxvNZnNlH9pkMllbdpvT09N49uzZ/MuOCwDweFeJKlcODg5iPB5HmqZrk1EkqgDkz3g5QP3kMUZufBx2xzM51IN5MGyLfoHiGQ6H0W635/F9fn4ezWYzkiRZOC5N0+h0OneWAQBQbO7Ji+deO5Zb4ZFtuvowsKIEQD7ySlZ5/fr1Y08VeIDbOtuPjo6WOs6vd6ivK6OePLAD7Nbl5WUMh8N5P9q64+5bdnJyEsfHx/PX0+lUcjkAPNK6RJXz8/OVySifP39eWdZutzM/b27nmRiqSWwD1FvWY+TGx2F38ngmn81mMZvN5q8tHlFt+/vmshaNeTBQXWmaRq/XW3iv2WzG0dFRREQMBoPo9/uxt7cXFxcXC3kn68oAqBY5Z5CPvMfSjN0Vx8aJ5VcrPF4llh8cHES/319awTHCCo8Au+ZCyzqSVaB6dLaXUxmu1zrmYPvKEPvsTr/fj0+fPs0ntzSbzaUJq5PJZD4xZlXZbRqNRjQajSxOu9ZcKwHIYzFHE9kBAOBxshojNz4Ou2XxCKgu82Cqwdg4qyRJEj9//lxbfnZ2FhGL1/u7ygDYLnNiAKpt48RyKzwCQDVIVoHq0dkOAOX35s2b6Pf7810XIiK63W68fft26dhOpxNJkqwsAwDyl+VijiayA5SXSfQAxZDVGLnxcSgGi0dA9ZgHAwAAANl6cp+D81jh8dmzZ/MvnXAAsF03k1UuLy+j2+3eemyn01lbBgAAbMdwOIx2uz2/Tz8/P49msxlJkiwcl6ZpdDqdO8sAgPxluZjjyclJ/PjxY/717du3LH4FAACoJGPk5be/X7zFWop4TnWW5eIRT58+XfgCAAAAgCq4V2L5lasVHgeDwZ3H3afMxBgAyI5kFQAAbmMC3G6laRq9Xi9evHgRv/32W/zxxx/R7/fn5YPBIPr9fgyHw3j79u1Cf9y6MvIljgDqLetEFRPZAQDgYYyRl5s+NzZh8YhyEtsAAAAAu/X7Q/5Rlis8NhqNh5wSAAVzNQDw8eNuz4NfrpJVrms2m3F0dBQR/5OQsre3FxcXF7cmq9xWBgAAPE6SJPHz58+15WdnZxERcXBwsHEZAJCP2xJVjo6OlsbBriejrCsDAAC2wxg5VJ9ncgCA/FkgBgCgGu6dWH5zhceIiG63G2/fvl06ttPpRJIkK8sAgHxIVgEAyJeBNACA6pOoAtU3Ho9jNBpFRMTFxUW8f/9+nnCSpmkMh8NIkiTSNF1IYFlXBgDkwxg5VJtncgCAfJkHAwDV4/peb/dKLLfCIwAA1109THz8uNvzAAAez3UdAOB+JKpA9Y1Go/jrr78i4tcC7H/++Wd8+fIlIiJ6vd78+zRN4/DwcJ6Qsq4MAAB4PM/kAAAA+ZB4ClBNGyeWW+ERAAAAAACoGwPlAPU0Ho/j9PR0nlh+cHAQ/X4/0jRdOjZJkvnO5jfLr5cBAAAAVJ0+dQAAgOLbOLHcCo8AALth91AAAACoN30DAPlrt9vx/v37+evLy8uIiGi1WnF+fh6tVmvh+FarFePxOD5//ryyrN1uL/1/ZrNZzGaz+evpdLrF3wIAAKA4JJsCAAAAFMPGieUAAAAAAAAAUBfXF03/8OFDdLvdaDab8yTzmyaTydqy25yensbr168fe6oAAAAAbODmIhcWdIVqsVgzAMBmJJYDAPBoOuMAoHzsCgEAAMWkr614Li8vYzgcxpcvX+487r5lJycncXx8PH89nU7j+fPnDzlNAAAypl8dAAAAgLLTx0WExHIAAAAAAMjMeDyO0WgUEREXFxfx/v37aDabERGRpmkMh8NIkiTSNI2jo6ONygCAfPX7/fj06dP8WtxsNpd2IJ9MJtFsNteW3abRaESj0cjitAEAAAAAAGCrzIOBapBYDsy5uAMUm5WhAACgGNybcx+j0Sj++uuviIh48+ZN/Pnnn/PdTnu93vz7NE3j8PAwBoPBnWUAQH7evHkT/X4/kiSZ7zre7Xbj7du3S8d2Op1IkmRlGfV29Rzx8eNuzwMAAAAAqsg4PkA+zIOBapBYTmHs7xtE3zUXdwAAAADYnvF4HKenp/M+t4ODg+j3+5Gm6dKxSZLMF328WX69jHqTkAaQr+FwGO12e55Ufn5+fusCy2maRqfTme9YvqoMKB6LrwMAAAAAwGbMg4HqkFgORISLOwAAAABsW7vdjvfv389fX+1y2mq14vz8PFqt1sLxrVYrxuNxfP78eWVZu91e+v/MZrOYzWbz19PpdIu/BQDUU5qm0ev1Ft5rNptxdHQUERGDwSD6/X7s7e3FxcXFwqLL68qAYrH4OgAAAADAoqsFz+GmPObBmAMD+Xiy6xMAimHdxX00Gq28gK8rAwAAAKgag2fc18HBwfz7Dx8+RLfbjWazOe9/u2kymawtu83p6Wk8e/Zs/vX8+fPHnjYA1F6SJPHz58+Fr+/fvy+Un52dxcHBQZydnS3sVLyuDCiOq8XXrxwcHMR4PI40TdcusG7xdYDs6YMDAAAAKKas58GYAwP5kFgOzOUxyXU2m8V0Ol34AgCgHPb3TeQBAHioy8vLGA6Hd+5iuKq/bV3ZyclJ/PjxY/717du3R5wpQLV5tgXgSh6LrxsfBwAAAACgirKaB2MODOTj912fAJi8UzxXF/cvX77cedx9y05PT+P169ePODsAAAAANnXV9/bx427Pg4h+vx+fPn2a71jabDaXFmecTCbRbDbXlt2m0WhEo9HI4rQBAErJGDSbymNnFePjAAAAQFnpZwNglazmwZgDU0/muOVPYjmwJMtJricnJ3F8fDx/PZ1O4/nz51s9f4Cy0xEHFI3PJQAoH9fvYnnz5k30+/1IkmSehNLtduPt27dLx3Y6nUiSZGUZANtRhmulwXOA/GS1+LrxcQCA8vAcDgCwHePxOEajUUREXFxcxPv37+f5JWmaxnA4jCRJIk3TODo62qiM+nKfXkzmwUD5SSwHFmR9cbdyTP24kQcAAADqbDgcRrvdnve3nZ+f3zoAnqZpdDqd+WKOq8rITxmSTgGA8jF2Vkx2VmFbxDgAAABQd6PRKP7666+I+JWf8ueff84XdOz1evPv0zSNw8PDGAwGd5YBxWEeDFSDxHJgzsW9/Ex2BQAAACiONE2j1+stvNdsNuPo6CgiIgaDQfT7/djb24uLi4uFQfF1ZQAAwPbYWQUAAAAAtmM8Hsfp6ek8sfzg4CD6/X6kabp0bJIk853Nb5ZfLwOKwzwYqA6J5UBEuLgDAAAAwLYlSRI/f/5cW352dhYRvwbUNy0DAAC2w+Lr5WbhdagecQ0UydVn0sePuz0PAIAyabfb8f79+/nrq4UcW61WnJ+fR6vVWji+1WrFeDyOz58/ryxrt9sL789ms5jNZvPX0+l0y78FsIp5MFAdEssrSocW9+XiDgAAAAAAAEBdWHwdAAAAALbvek7Jhw8fotvtRrPZnCeZ3zSZTNaW3XR6ehqvX7/exqkCQG1JLAeoACs2AwAAAAAAQHlYLH73LL5eXkUcHy/iOQHZcz0HAABY7fLyMobDYXz58uXO4+5TdnJyEsfHx/PX0+k0nj9//tDTBIBaklgOAAAAAAAAAFSC5E4AAOAxPFMAAGxHv9+PT58+RbPZjIiIZrO5tAP5ZDKJZrO5tuymRqMRjUYjq9MGgFp4susTAAAAAAAAAAAAAAAAAHis/X2Lxezamzdvot/vR5IkcXl5GZeXl9Htdm89ttPprC0DALbPjuUUytXN+8ePuz0PAAAAAKgS/W4AAFSdiaIAAACQP2NQANw0HA6j3W7Pk8rPz8/j6OhoaffxNE2j0+nMdyxfVQYAbJ/E8oozeA4A5Gl/3yABAJSd6zkAFJ+JegAAAFAcntMBoLxcxwG2K03T6PV6C+81m804OjqKiIjBYBD9fj/29vbi4uIiBoPB/Lh1ZQDUg/mr+ZFYDgAAAAAAAFBQJrcCAABQZjZIAgCojyRJ4ufPn2vLz87OIiLi4OBg4zIAYLsklgMAAAAAAAAAAAA7IekUAAAAds/zOUB9SCwHAAAAAKiJ/X27nQIAAAAAAPA4Es8AAKBc3MNzncRyAACAjF09iEviAsrC5xYs07EOAAAAAAAAAAAAlN2TXZ8AAAAAQBVIOgUAgGLa33e/DgAAAAAAm9KvDgBQbXYsBwAAAAAAAAAAAAC2RjIaAAAA7J7nc24jsRwAACBHHs4BAAAAAAAAAAAAAIBdeLLrEwAAAKB89vclyQMAAADAtuhvAwAAAAAAAPJgx3IAcnE1Eebjx92eBwDskomhQNm4jwcAoEo8lwMAAABANvS9AQAAlIcdywEAAAAAAAAAAAAAAAAAACrOjuUAAAVh1VagLHxeAQAAUAaeXwEAoJiu7tU/ftzteQBswmcWAAAAUDUSywEAyISBNQDqQrIKAFTH/r7nWAAAAAAAWMX4OFRb1WLcPFaoDvEMsF0SywEA2KqqdSwCwH24DgIAAFkxYQYAAIAyMF4GABSVfnYA6qSsz+eu1/mQWA4AsENuegEAAICH0q8AAAAAAAAA65U1sQ4AICtPdn0CANSLB3MAAAAAAAD4ZX/f+BkA9eU6CAAUhfsSACgH12uA7ZBYDgAAAAAAAAAAAABwB8mnAAAAQNn9vusTgNtcdbp9/Ljb84Ci00ENAAC7434cKCt9b7CeazxQdPv7ruMAAAAAAAAAwMPYsRwAAAAAAAAAACissu0MWrbzBbbLZwAAAAAAUGQSywEAAAAAAAAAAAAAAAAAACru912fAAAAAOV1tdvCx4+7PQ8A4P5cxwEAACgbOwADQHG5TgMAQDm5lweoH4nlFbK/bxIoAABAFY3H4zg8PIwvX74svJ+maQyHw0iSJNI0jaOjo2g2m3eWAQAAAFSFCW8AUB6u2wBQXnW5jluYGQCAsqvLvTuPI7G8Aq4Hu8AHAKgnSadQXVdxOh6Pl8p6vd487tM0jcPDwxgMBneWAXA3/WwAAAAAAAC7YR4MVJsYZ9csIgEA1J3EcgAAKDlJp1BtBwcHt76fpunC6yRJYjQa3VkGAED1mPwCAAAAANmyIG1+zIOBahPjxeQ6BwBQLxLLAUrIwzsPYYVHdsXk9uxJOoV6Go1G0Wq1Ft5rtVoxHo/j8+fPK8va7fatP282m8VsNpu/nk6n2z9p4Fbu1QEAALiiTx1ge/S7AUA5mQcD1SbGAQBg9ySWA0ANWOER6knSKVTb5eXlre9PJpO1Zaucnp7G69evt3BmwH24VweAcpOoArA7VV+IeX9fcjnAY+h3g+rzTA71s+15MECxbDPGzXMDAID17pVYriMOAMrJCo9QT5JOoZ5WxfddZScnJ3F8fDx/PZ1O4/nz51s8M+A27tUB2CY7nOZLogoAAJtyr54//W7FUvUFYa7U5fcsAs/kFIHFoPK3zXkwkk6heLYZ4+a5AQDAehsnluuIA4DqsYprcRhgJk+STvMjtslSs9lcGiCbTCbRbDbXlq3SaDSi0WhkcaqVI7bJw7bv1U2OAYDtk6gCAADlYxdE8iTpNDueyYHrHjIPRtLp/RgjZ5ceEuPmuQEAwHpPNj3w4ODg1g5yHXEAUF7b3s14NpvFdDpd+AJ2J4uk06dPny58cbv9/foNqtXxd961brd76/udTmdtGVAO275XPz09jWfPns2/DJoDQHbWJaqsK1tFnxsAAGzPtndB1OcGxbLNZ3LP41A825wHc3JyEj9+/Jh/ffv2LavTBja0zRg3zw1gM+Z9AtTXxonlq+iIA4DqeehuxgbOoVgknUL1XL8OJ0myUJamaXQ6nWg2m2vLgHJ76L26yTEAkB8LxJAHE33yMR6P4x//+MfS+2maxps3b2I4HMabN28WYntdGQD145pdHg/dBVGfGxSLxSOg2rY5D0bSKRSPuW4A3Je+N4CH+/2xP2DbHXGvX79+7CkBABva9m7GJycncXx8PH89nU4NrEHOLi8v53G6LrH0ZixLOoXiGo1G8enTp4j49dy8t7cXBwcHERExGAyi3+/H3t5eXFxcxGAwmP+7dWVA8W37Xr3RaESj0cjiVKmAq0G2jx93ex4AVfeYBWL0uUH+hsNhJEly66LpvV4vvnz5EhG/+tUODw/nz93rygC4PxND2bZt74Kozw3K4aGLR3geh90zD2a33I+TNTFeDGIdACgic9qy9ejE8lV0xAFA8XW73Xj79u3S+51OJ5IkWVm2ioFzNuEGf/sknUK1dbvd6Ha7cXZ2tlSWJMn8/au436QMKL5t36sDAPmxQAxUw6pn6TRNF14nSRKj0ejOMgCgGPS7QbVZPAKqwTwYqDYxDsBjWRAD4PEenViuI44sSTwD2D4rPEL1SDoFgGpwrw7VNB6P4/DwcL5z6ZU0Tec7oaZpGkdHR/M4XlcGFJ9Eld0zmYQsjUajaLVaC++1Wq0Yj8fx+fPnlWXtdjvP06w0MQ7Afel3g/rwTA7VYB4MVJsYB4ByMw+muIyhcR+PTizXEUceJJjnw8UdqssKjwDZcr8K9aDTjSy4V4dqu+o3G4/HS2W9Xm/eD5emaRweHs5jeV0ZUEwSVdgVfRL5u7y8vPX9yWSytmyV2WwWs9ls/no6nT7m9IB7MkYO1aXfDerDMzkAAADkwzwYduHmmLgx8u14UGK5jjioHhd3qDYrPAIA8Bg64rLjXr14LCLBNq2KzzRNF14nSRKj0ejOMqBYJKoA161KKr+r7PT0NF6/fr39EwLuZIwcqk2/G1SbZ3IAAADIn3kwUB0bJ5briINqc3EvB5PbobwkYwFAebkPB+rG80v2RqNRtFqthfdarVaMx+P4/PnzyrJ2u53nadaK6z0PIVEF6qnZbC7tQD6ZTOYLr68qW+Xk5CSOj4/nr6fTaTx//nyr5wzczhg5AJSXZ3IAAAAoDvNgoHw2TizXEQf1tO2L+2w2i9lsNn89nU63f9IAAOyMRDQAgLut2rV0MpmsLVtFnxs37e+7JwfISrfbjbdv3y693+l0IkmSlWWrNBqNaDQaWz1H4HFMgAMA2JwFGwEAuGk8Hsfh4WF8+fJl4f00TWM4HEaSJJGmaRwdHc0XZl1XBhTfNufBmAPzMHV4Pq/D75injRPLgXra9iTX09PTeP369RbODAAAAACqZVV/211l+twAfjGQTFYuLy/nE9iSJFkoS9M0Op3OfMfyVWU8nhgnLybAAbBNFmUGgPJyHQe4v6vk8PF4vFTW6/XmyeZpmsbh4WEMBoM7y4Dyesg8GHNgIB8Sy0usjgPnHtCL46GTXE9OTuL4+Hj+ejqdxvPnz7d4ZgAAAAD3V8e+Nnan2WwuJZ5MJpN5QtqqslX0uQHA9o1Go/j06VNE/JrAsre3FwcHBxERMRgMot/vx97eXlxcXCxMbltXBpSbCXAAAABAlRgjJwtX/eg3pWm68DpJkhiNRneWkQ3xz7Ztcx6MOTCQD4nlwFrbnuTaaDSi0WhkcaqV5+YdqkVMA1VnUSgAgNW63W68fft26f1OpxNJkqwsW0WfGwBsX7fbjW63G2dnZ0tlSZLM3785SW5dGVAOJsCRJX3nAABQDtfn97mPB3i80WgUrVZr4b1WqxXj8Tg+f/68sqzdbi/9rNlsFrPZbP56Op1mc9LAxrY5D8YcGMjHk12fAFBs3W731vc7nc7aMgAAiLCQBgDAles7GCZJslCWpml0Op1oNptrywAAgOxtc4y80WjE06dPF74AAKiO/X1j4gAAm7g+Xn7dZDJZW3ab09PTePbs2fzLQo6wG+bBQLnZsRxYcnl5Ob9Ar7uA37yIu7gDAHVnwBiqR1wDdeIzb/tGo1F8+vQpIn4Nbu/t7c13Lh0MBtHv92Nvby8uLi5iMBjM/926MgAA6sOOaPkxRg4AAAAA+VuVVL6u7OTkJI6Pj+evp9Op5PKa05eeH/NgoDoklgMR4eIOQPYkqQBA+emEB7ifbrcb3W43zs7OlsqSJJm/f9UPt0kZbMp1GwBgPWPkxWZcCQAAAMja/r6xtLw0m82lHcgnk8l8McdVZbdpNBrRaDSyOlVgDfNgoDoklpeQwTOy4OJOnkxqpU5ctwEAAKC4PLcDANSXMXIAAAAAyEe32423b98uvd/pdCJJkpVlAKxn3gsPJbEcoOBc5AEAAAAAYDX96AAAAMAu2WgFAMrLdTw7l5eX813HkyRZKEvTNDqdznzH8lVlAEA2JJYDAAAAAAAAAAA7Z8EYqB5xDdUktgEAuM1oNIpPnz5FRMTp6Wns7e3FwcFBREQMBoPo9/uxt7cXFxcXMRgM5v9uXRkAsH0SywEAAAAAAABKzI4qAAAAAGyLxSOg2sQ4Wep2u9HtduPs7GypLEmS+ftXyeablLE94h+AK092fQIAAAAAABTP/r5BRQAAAAAAAACqz/g4AFAnEsspJTftUA1imTrS7pepDwAAAAAAAAAAys68IAAAAKAMJJaXgI6m1a7qRv0AAAAAAAAAUEXGxKF6xDTUj+s5AAAAu+SZFIDrft/1CQDA1UPKx4+7PQ8gXzooAAAAAAAAAKrJeDBUk9gGAACKSl4K1IuYfxyJ5QAAAADciw45AMrOBFgAAAAAAFik7xwAAKAeJJYDAAAA8CASzKkCE2QAAAAAAACAuxgfp+yMjQMAcEViOQBARnTCAQAAAADA4+hrB4Byce0GgPKQKA0AAFBPT3Z9AgAAAAAAFJfJwFTJ/n5923Sdf3eqTbsGAADKxPM5AAAAALBrdiwHAAB4IJM+NmOFawAAAAAAAAAAAAAA2D2J5QAAAAA3WDgCAKgK9zVUlbYN1JEFHKFaxDRVoS0DQPnpawMAgHJw7862SCwvEYEPAADF4f4c4H/s75s0CAAARSCpBwAAAACA68xzA+rGeBnA3SSWAxRMnR/e3cADAAAAAABwpc7jZgBQNq7bAP/DPDgAKC/XcQAoFxsjPYzEcgAAdk5HHAAAAACbkKyyGf1tAAAAAAAAGFsD4DZPdn0CAAAAALuyv7/YeX7zNVBt4h0AAAAAAMiaMUioFjEN9SDWofzEMFXgekRW7FgOAAAAwKPZEZEy0E6hfsQ9AJSLazfUl4lxAFBsrtUAq+nPAAAAykZiOQCFo5MN6kv8Q7WJcQAAdsXEVwCgivS3AQAAcF/6y6E+xDsAUBfGzO5PYjlAQXh4BwAAAIpK5zsAAHkybrae+gGA8tPfBgDl5/kcAAAoK4nlVI5OdwB2TYcxAAAUm3t2AAAAAAC4m/70x1F/AFAd8lQAgCqRWA4AAHAPBn4fTyc7AAAA9+V5HGC1/X19bZSTvmIoN/foAAAAQNHpgwS4ncRyAAAAoPZMgAPYjAE3AAAAtq3ufXOetSmbusfstvkMgHqxKBQAAHnx/A7VIqbZNonlBSbg70d9QfXoSKdsXIug2sQ4AABQBSasAwBAMblXBwAA4LE8W26feYOUjTYLwCae7PoEAAAAAAAAANi+/X0TiAAoDtclAAAAyIfnbwAA1rFjOcCOeXBfz8p5lIE4BngY13kAAAAAAIDyMDYO9SHeoXrENQAAFJt7dvIksRwAAACArbnZuWnxCHZJZzvUm88AgEUWeAOAcnHtBgAAgGLyzA4AxeZafTeJ5VSWDwCKzsRWgLu5ngOQFffju7O/79oOVWARCQCypl8IAKgrfZdQT/rOAQCgHDy3U1Ta5nrGHgEWSSwvKBd0gEVu5AHIi2sOAABQB559KCpjZACP4xoPAMBjeTaHavB8CAAA1M3NPg3PRatJLC8IHXHZ8QEAAMBDuEffHTsyAAAAwHbp54DqE+eUifYKsJ75blBtN++FxDpUx/X4FtsAAECRSSwHANiQwdv83LVSlKRTqBafr+TJpNX8iXHyYqIGAFSTe/jtc49OkWiPUA+u5wCb8XkJwGN4xgYAALid56VlEsupDR8AUA1imSIwmAsAcH/uoaDaPK8DQPm4flMknhkBVnPNpghcq6EexPpuuNYDAPAY7uMBeAiJ5QAAtzBoA1AsPpd5CO0G6sug2e747IXyEbcUiWs4VNv1a454BwC4P8/wUG1iHAAAALK1v++5+4rE8pzdHCDXEAGg2AzaQLWI6WK7/rxkcjHbJPaLz9+Ix3DNKCad8OyCzwMAKAfX7Hx41gbgMVxHoJpuxrZ782LyGQzVsCpnQIzzWK7fxSKmAXgM1xF2SWI5AACltOpBygMWDyHxCerB4Fpx+FvwGO73ALLhuQhwn8W2aVPF4u8B5eeenV3QlwvVIJah3nwGQDWIZeAmnwsAPIbE8pysumC7kO+GQXPyps0BbMe6eyf3VdzHbe1FG4LyE8dQL2K++FbtfqN/BIpPvJI31/Xduln/PgOgWnzGskvaH5SHeC0mfxeyom0BAEC5uad/uJt1ZzyMXRDDxVD3cXGJ5QA5qvtFB8rMzXtx+dsA4FpQPZ6duI1YBwCAcnIvD/Uj7vOjH41tELPlc/1vJv4BuMm1vZjWLe7ovh5JjuXnb0heXOe3T/xSBGKbXcglsTxN0xgOh5EkSaRpGkdHR9FsNvP4X++Eh7tiW7c7pb/Zw9Qtxu9j1cXdRX/7xHF26hbj4pO6qXKM7+8v744JdVPlGAfEOOXhXuxhxPjd9AdRZnWL8VXx6hpRfD5rH6ZuMU65iOvHE+OLtCmqRoxD9dUlzl2jqau6xPgV/WvVIrHtbmIcqq1uMX4fPg+oAjFOUaxb9Om211WVS2J5r9eLL1++RMSvD4HDw8MYDAZ5/K93qi6NqMzuSvr1t9tMXWP8NtpOftY9HPk7bJcYh2oT45RNXR/eH0qMQ7VVMcZXfc5DHVUxxoH/UfcYd40vB3+nh6tDjN/WJ6PNlMv1hTm5nzrE+EPoq82Hz9rsiXHKZtXngs/j1aoS56691XfbYvr+3nerSowDtxPjlMX1+3TX782Jcai2use4vl2KJvPE8jRNF14nSRKj0Sjr/22u7uqwEfjlscnfSgfdojrE+EOI+91R99tVtRg30a1erv+9Xb9vV7UYv3I9rsV4PYn5X6oa41fEdzWtWwCu7jF9kxinzEx2vVvVY/whNllkkMe7T1267364OsW4+KwO1+/N1SnGr4j18rMj2ubqGOP35TNht9ynP05dYlyc1oPPg9tVIc5X7aq16fGUi7/f/ZQ1xtfNabv5Oa5NVI+/6ebKGuM3bTL/Qbuolvvem9d1jkxVYhy4XR1jXN9M+dRt/mrmieWj0ShardbCe61WK8bjcbTb7YX3Z7NZzGaz+esfP35ERMR0Os36NDf2z3+uLvuP/8jvPMjHf/xHxPn54nv/9V+//rtJs7xqLzd/Rhau4uTnz5/Z/8+uuU+MR5Qjzu/r+t/5qn1QDFfNKs9YzEqdY/w+f7+bx958fRWjrtn1cv3vfbMp3ry329XnRJ1jfFPX/1Y3Y5x6uvlZfr0Zrrpnf8g15cpjPh/qHOMPqcdVfycxX0/XY/1m3Bfl/r7OMf4YYprbrGuOt10f7noG3M45ifGs3XXtv/lMT7426UN5zH33rtUlxtfF2X3/Tvf5N7dd7/XLVc/1fviixX1dYjwrm96zi+tyW/X3u3r/trje9P7t+vubHntfu4jzqsT4ptZ9vnu2L5ab8Xw9jst0j36dGF9038/f28rEbT3dp2/dXLfdzGfd5HnK/DSurIrp2z7j77qPyzLW6xDjj5l/cOVm//e6eQ+u4/V023P7LmP7Sh1i/DHu+lv8138t9724vtfDzb/7qrZyWxt5zFzq+xLju7Fqnp17gPzd5/pbRmJ8vU374O4Tm8bSyuu2+atZjXlty31iPPPE8svLy1vfn0wmS++dnp7G69evl95//vz5tk8LNvbs2f3ef+yxj/Wvf/0rnuX4P7xPjEdUO87z/DuzmZt/kyr8jeoc44/53K3C357tuKst7Lqt1DnG72PXfyeK6bZ2sY17+cf8m5vEeHGfoyinorURMQ6Pt0kIbXLNzyIUxXj2tnnvRr6q8LerS4zf57npvj+H+rreHoraNuoS45CFdaFzn3uArO8X8ozzOsZ4UT/f2UwZrtV3EeOLHvKZWta/Pdl4aB9cVop8v55njG9aBeKZmx7z+Z9He6pDjD/m19tkfEPcc1MRYvtKHWL8MYp230Vx3Gd8+zFj4Y9tX2J8t3w+FEsV/x5ifL0qzIMgG2VpG5vEeOaJ5avc9oFwcnISx8fH89f//d//HZPJJP7X//pf8dtvv638WdPpNJ4/fx7fvn2Lp0+fZnG6XKO+87dJnf/8+TP+9a9/xb//+7/nfHa3W3XRf0ica3P5Ut/5E+PaXJ7Ud/7EuDaXJ/WdPzGuzeVNnedLjGtveVPn+RLj2lve1Hm+xLj2ljd1ni8xrr3lTZ3na9P6LlKci/FyU+f5EuPaW97Ueb7KGOMR25vPqr3lT53nS4xrb3lT5/mrSr+bvJRyUN/5E+PaXJ7Ud/7EuDaXJ/Wdv23HeOaJ5c1mc2n1iMlkEs1mc+nYRqMRjUZj6d9v6unTpxpijtR3/u6q8zxXi7lynxiPeFyca3P5Ut/5E+PaXJ7Ud/7EuDaXJ/WdPzGuzeVNnedLjGtveVPn+RLj2lve1Hm+xLj2ljd1ni8xrr3lTZ3na5P6zjvOxXi1qfN8iXHtLW/qPF9FjPGI/Oazam/5U+f5EuPaW97Uef7K3u8mL6Vc1Hf+xLg2lyf1nT8xrs3lSX3nb1sx/mRbJ7RKt9u99f1Op5P1/xrIgRiHahPjUG1iHKpNjEO1iXGoNjEO1SbGodrEOFSbGIdqE+NQfeIcqk2MQ7WJcag2MQ7VJsaheDJPLE+SZOF1mqbR6XTutVIEUFxiHKpNjEO1iXGoNjEO1SbGodrEOFSbGIdqE+NQbWIcqk2MQ/WJc6g2MQ7VJsah2sQ4VJsYh+L5PY//yWAwiH6/H3t7e3FxcRGDwWCrP7/RaMT//b//NxqNxlZ/LrdT3/krep2L8WpR3/krep2L8WpR3/krep2L8WpR3/krep2L8epR5/kqen2L8epR5/kqen2L8epR5/kqen2L8epR5/kqen2L8epR5/kqen2L8epR5/kqen2L8epR5/kqQ31nGedl+P2rRp3nqwz1LcarRZ3nr+h17n69WtR3/ope52K8WtR3/ope52K8WtR3/rZd57/9/Pnz51Z+EgAAAAAAAAAAAAAAAAAAAIX0ZNcnAAAAAAAAAAAAAAAAAAAAQLYklgMAAAAAAAAAAAAAAAAAAFScxHIAAAAAAAAAAAAAAAAAAICK+33XJ7Bt4/E4IiLa7XakaRqXl5fRbrd3fFbVkaZpDIfDSJIk0jSNo6OjaDabuz6tytKel6mTbInxfGnPy9RJtsR4vrTnZeokW2I8f9r0IvWRLTGeP216kfrIlhjPnza9SH1kS4znT5tepD6yJcbzp00vUh/ZEuP506YXqY9sifH8adPL1El2xHj+tOdl6iQ7Yjx/2vMydZItcZ4v7XmZOsmWGM+X9rxMnWRLjOdLe16mTrIlxvOVWXv+WTFHR0c/I+JnRPzsdrs/v3//vutTqpR2uz3//uvXrz8PDg52eDbVpz0vUyfZEuP50p6XqZNsifF8ac/L1Em2xHj+tOlF6iNbYjx/2vQi9ZEtMZ4/bXqR+siWGM+fNr1IfWRLjOdPm16kPrIlxvOnTS9SH9kS4/nTppepk+yI8fxpz8vUSXbEeP6052XqJFviPF/a8zJ1ki0xni/teZk6yZYYz5f2vEydZEuM5yur9ly5Hcv/8Y9/xPfv3yMirHSwZWmaLrxOkiRGo9GOzqYetOdl6iQ7Yjx/2vMydZIdMZ4/7XmZOsmOGN8NbXqR+siOGN8NbXqR+siOGN8NbXqR+siOGN8NbXqR+siOGN8NbXqR+siOGN8NbXqR+siOGN8NbXqZOsmGGN8N7XmZOsmGGN8N7XmZOsmOOM+f9rxMnWRHjOdPe16mTrIjxvOnPS9TJ9kR4/nLqj0/2dpPKpBmsynoMzAajaLVai2812q1Yjwe7+iM6kF7XqZOsiHGd0N7XqZOsiHGd0N7XqZOsiHGd0ebXqQ+siHGd0ebXqQ+siHGd0ebXqQ+siHGd0ebXqQ+siHGd0ebXqQ+siHGd0ebXqQ+siHGd0ebXqZOtk+M7472vEydbJ8Y3x3teZk6yYY43w3teZk6yYYY3w3teZk6yYYY3w3teZk6yYYY340s2nPldiy/vLyM4XAYEREXFxfx6tWrSJJkx2dVDZeXl7e+P5lM8j2RGtGel6mT7Ijx/GnPy9RJdsR4/rTnZeokO2J8N7TpReojO2J8N7TpReojO2J8N7TpReojO2J8N7TpReojO2J8N7TpReojO2J8N7TpReojO2J8N7TpZeokG2J8N7TnZeokG2J8N7TnZeokO+I8f9rzMnWSHTGeP+15mTrJjhjPn/a8TJ1kR4znL6v2XLnE8qOjo3n2fZIk8eLFi/j69etuT6riVn0g8Hja8zJ1kj8xnh3teZk6yZ8Yz472vEyd5E+MZ0ubXqQ+8ifGs6VNL1If+RPj2dKmF6mP/InxbGnTi9RH/sR4trTpReojf2I8W9r0IvWRPzGeLW16mTrJlxjPlva8TJ3kS4xnS3tepk7yJ86zoz0vUyf5E+PZ0Z6XqZP8ifHsaM/L1En+xHh2smrPpUgsf/fu3dpf9sWLF9HtdiMiIk3TaLfbEfGrotI0jTRNrSqxBc1mc2n1iMlkMm+YbF9d2rMYLwYxnr+6tGcxXgxiPH91ac9ivBjE+G7UoU2L8WIQ47tRhzYtxotBjO9GHdq0GC8GMb4bdWjTYrwYxPhu1KFNi/FiEOO7UYc2LcaLQYzvRl3atDjfPTG+G3Vpz2J898T4btSlPYvxYhDn+atLexbjxSDG81eX9izGi0GM568u7VmMF4MYz19W7bkUieVHR0cbHTcej+PPP/+M79+/L7zfarWyOK3a6Xa78fbt26X3O53ODs6m+urUnsV4MYjxfNWpPYvxYhDj+apTexbjxSDG81eXNi3Gi0GM568ubVqMF4MYz19d2rQYLwYxnr+6tGkxXgxiPH91adNivBjEeP7q0qbFeDGI8fzVqU2L890T4/mrU3sW47snxvNXp/YsxotBnOerTu1ZjBeDGM9XndqzGC8GMZ6vOrVnMV4MYjxfWbbnJ4/+CQWSJEmcnZ3NX49Gozg4OLDiwZbcXMUgTdPodDrqNyPa8zJ1ki0xni/teZk6yZYYz5f2vEydZEuM50+bXqQ+siXG86dNL1If2RLj+dOmF6mPbInx/GnTi9RHtsR4/rTpReojW2I8f9r0IvWRLTGeP216mTrJjhjPn/a8TJ1kR4znT3tepk6yJc7zpT0vUyfZEuP50p6XqZNsifF8ac/L1Em2xHi+smzPv/38+fPno39KgYzH4xiNRtFsNuPr168LFcfjpWkab9++jb29vbi4uIiTkxOBnyHteZk6yZYYz5f2vEydZEuM50t7XqZOsiXG86dNL1If2RLj+dOmF6mPbInx/GnTi9RHtsR4/rTpReojW2I8f9r0IvWRLTGeP216kfrIlhjPnza9TJ1kR4znT3tepk6yI8bzpz0vUyfZEuf50p6XqZNsifF8ac/L1Em2xHi+tOdl6iRbYjxfWbXnyiWWAwAAAAAAAAAAAAAAAAAAsOjJrk8AAAAAAAAAAAAAAAAAAACAbEksBwAAAAAAAAAAAAAAAAAAqDiJ5QAAAAAAAAAAAAAAAAAAABUnsRwAAAAAAAAAAAAAAAAAAKDiJJYDAAAAAAAAAAAAAAAAAABUnMRyAAAAAAAAAAAAAAAAAACAipNYDgAAAAAAAAAAAAAAAAAAUHESywEAAAAAAAAAAAAAAAAAACpOYjkAAAAAAAAAAAAAAAAAAEDFSSwHAAAAAAAAAAAAAAAAAACoOInlAAAAAAAAAAAAAAAAAAAAFff/ASiaN3b1YVWQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 4000x250 with 16 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot histograms\n",
    "\n",
    "fig, axs = plt.subplots(1, len(ring), figsize=(17*2, 2))\n",
    "for dqx, ax in zip(dqxs, axs):\n",
    "    ax.hist((dqx - dqx.mean()).cpu().numpy(), bins=100, range=(-5.0, +5.0), color='blue', alpha=0.7)\n",
    "plt.tight_layout() \n",
    "plt.show()    \n",
    "\n",
    "fig, axs = plt.subplots(1, len(ring), figsize=(len(ring)*2.5, 2.5))\n",
    "for dqy, ax in zip(dqys, axs):\n",
    "    ax.hist((dqy - dqy.mean()).cpu().numpy(), bins=100, range=(-5.0, +5.0), color='blue', alpha=0.7)    \n",
    "plt.tight_layout() \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9969486d-a58a-41f9-a77f-9bb376257db3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAC+CAYAAAB9JI0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATuElEQVR4nO3dTU4jSRoG4I9WS6yazvbsWlObnBsY+wAjXDcwR4Daz8IWq1GtLHwDuk7Qbd+gUpoD0OQNyA1SL01CrZBG7VmUsDD+wQb3UImfR7JERgZZUSGCdPglI3bG4/E4AAAAAAAAvnHfvXYDAAAAAAAAViHUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEr4/rUbsI4///wz/vjjj/jhhx9iZ2fntZsDAAAAAABswHg8ji9fvsTPP/8c3323+HmMSoUaf/zxR7x79+61mwEAAAAAAPwFrq6u4u9///vC85UKNX744YeI+Pqf2tvbe+XWAAAAAAAAm3B7exvv3r2b5ACLVCrUuF9yam9vT6gBAAAAAABvzFNbT9goHAAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCVUaqNwAODb889/vnYLXsd//vPaLQAAAIDtI9QA/lI+7AQA4K/ivSa8fcY5AI9ZfgoAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAS7KkBUBHbupZshPVkAQD+atv6XtP7TACoHqEGAABvkg/oNkdfAgAA3wqhBjyyrZP2CBN3tsu2jnXj/Nu2rT+XEX422R7GObx9xjm8fcY5vC6hxhvhlykAAAAAUDXb+rmmzzSfT6gBAAAAAG/ctn5wHOHDY3hrvnvtBgAAAAAAAKxCqAEAAAAAAFTCWstP5XkeR0dHcXFxsbReURQxHA4jTdMoiiKOj48jSZInzwEAAAAAACyycqhxH0Tkef5k3cPDw0nwURRFHB0dxWAwePIcAAAAAADAIiuHGu12e6V6RVFMHadpGlmWPXkOAAAAAABgmY3vqZFlWdRqtamyWq0WeZ4vPQcAAAAAALDMWntqrKIsy7nlo9Fo6bl57u7u4u7ubnJ8e3v70uYBAAAAAAAVtfEnNRZZFGgsO9fr9eLHH3+cvN69e/fXNA4AAAAAAPjmbTzUSJJk5smL0WgUSZIsPTfPyclJ3NzcTF5XV1ebbi4AAAAAAFARGw81Wq3W3PJGo7H03Dy7u7uxt7c39QIAAAAAALbTs0KNx8tF5XkeRVFERESaplPniqKIRqMRSZIsPQcAAAAAALDMyhuFZ1kWnz9/joive100m81ot9tTx51OJyIiBoNBdLvdaDabcX5+HoPBYHKdZecAAAAAAAAWWTnUaLVa0Wq14vT0dObc42AiTdNJvfvgY5VzAAAAAAAAi2x8Tw0AAAAAAIC/glADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAS1go1iqKIfr8fw+Ew+v1+lGW5sO5wOIyyLOfWyfM88jyfXPP+awAAAAAAgEXWCjUODw+j0+lEu92OdrsdR0dHS+v+9NNP8dNPP8XOzk7s7OxEv9+PiIizs7PY39+PnZ2d+PDhQ6Rp+rL/BQAAAAAA8OZ9v2rFoiimjtM0jSzL5tYtyzIGg0G02+1JWb/fj06nExER+/v7cX19HRERSZKs22YAAAAAAGALrfykRpZlUavVpspqtdrCpaMeBhrD4XDqOOJrmCHQAAAAAAAAVrXykxqL9s8YjUYzZQ/DirIsYzQaTS0xVZZlDIfDiIg4Pz9fuATV3d1d3N3dTY5vb29XbS4AAAAAAPDGrBxqLLJss/CIiG63G6enp1Nlx8fHk+AjTdN4//59XF5eznxvr9eLjx8/vrSJAAAAAADAG7Dy8lNJksw8lTEajZYuIVWWZWRZNlPn4f4caZpGURQze3ZERJycnMTNzc3kdXV1tWpzAQAAAACAN2blUKPVas0tbzQaC7/n999/nwk08jyPg4ODmbqP9+uIiNjd3Y29vb2pFwAAAAAAsJ1WDjUe73lRFEU0Go1JaJHn+czTFnmez4QVaZpOLUeVZVm0222bhgMAAAAAAEuttafGYDCIbrcbzWYzzs/PYzAYTM71er1oNpvR6XSmvudxGJIkSTQajej3+5EkSVxeXk5dBwAAAAAAYJ61Qo2HT1m02+2pc/OCiccBx716vR71en2dfxoAAAAAANhyKy8/BQAAAAAA8JqEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlfD9OpWLoojhcBhpmkZRFHF8fBxJksytm+d5RETU6/UoiiLKsox6vb72dQAAAAAAACLWDDUODw/j4uIiIr4GE0dHRzEYDObWPTs7i19++SUiIlqt1lS9da4DAAAAAAAQsUaoURTF1HGappFl2cL6+/v7cX19HREx9RTGutcBAAAAAACIWGNPjSzLolarTZXVarXJMlPzJEkys6zUc64DAAAAAACw8pMaZVnOLR+NRgvrD4fDiIg4Pz+PDx8+RJqma13n7u4u7u7uJse3t7erNhcAAAAAAHhj1tpTY55FIcXDzb/TNI3379/H5eXlWtfp9Xrx8ePHlzYRAAAAAAB4A1ZefipJkpmnKUaj0czyUvce7p2RpmkURRFFUax1nZOTk7i5uZm8rq6uVm0uAAAAAADwxqwcarRarbnljUZjpizP8zg4OJgpr9Vqa11nd3c39vb2pl4AAAAAAMB2Wnn5qTRNp46LoohGozF5wiLP80iSJNI0jTRN4/T0dFI3y7Jot9tzNw5/fB0AAAAAAIB51tpTYzAYRLfbjWazGefn5zEYDCbner1eNJvN6HQ6kSRJNBqN6Pf7kSRJXF5eTtVddh0AAAAAAIB51go1Hj6B0W63p849Dibq9XrU6/W1rwMAAAAAADDPyntqAAAAAAAAvCahBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJXy/TuWiKGI4HEaaplEURRwfH0eSJHPr5nkeWZZFRMT5+Xl8+vRpUjfP84iIqNfrURRFlGUZ9Xr9+f8LAAAAAADgzVsr1Dg8PIyLi4uI+BpwHB0dxWAwmFs3y7LodDoREdHv9+Pg4GDyvWdnZ/HLL79ERESr1Vp4DQAAAAAAgHsrLz9VFMXUcZqmkycxHsvzPHq93uS43W5HnueTa+zv78f19XVcX1/H58+fFz7tAQAAAAAAcG/lUCPLsqjValNltVptspTUQ/V6PT59+jQ5LstyUv9ekiTCDAAAAAAAYGUrLz91H0w8NhqN5pa32+3J17/++mu0Wq1JiFGWZQyHw4j4ut/Ghw8fIk3TmWvc3d3F3d3d5Pj29nbV5gIAAAAAAG/MWntqzLMo7Hh4fjgcTvbTiIipDcbTNI3379/H5eXlzPf2er34+PHjS5sIAAAAAAC8ASsvP5UkycxTGaPR6MklpLrd7sy+GQ/350jTNIqimNmzIyLi5OQkbm5uJq+rq6tVmwsAAAAAALwxK4carVZrbnmj0Vj4Pf1+P7rdbqRpGmVZRlmWked5HBwczNR9vF9HRMTu7m7s7e1NvQAAAAAAgO20cqjxeM+Loiii0WhMnsDI83zqaYvhcBj1en0SaPz222+RJEmkaRqnp6eTelmWRbvdtmk4AAAAAACw1Fp7agwGg+h2u9FsNuP8/DwGg8HkXK/Xi2azGZ1OJ4qiiMPDw6nvTZJkspdGo9GIfr8fSZLE5eXl1HUAAAAAAADmWSvUePiURbvdnjr3MJhI0zTG4/HC69Tr9ajX6+v80wAAAAAAwJZbefkpAAAAAACA1yTUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACohO/XqVwURQyHw0jTNIqiiOPj40iSZO2661wHAAAAAAAgYs1Q4/DwMC4uLiLiazBxdHQUg8Fg7brrXAcAAAAAACBijeWniqKYOk7TNLIsW7vuOtcBAAAAAAC4t3KokWVZ1Gq1qbJarRZ5nq9Vd53rAAAAAAAA3Ft5+amyLOeWj0ajtequc527u7u4u7ubHN/c3ERExO3t7fLGbqH//ve1W/B6Nv3joC83a1v7U19ulv7cHH25Ofpys/Tn5ujLzdGXm6U/N0dfbo6+3Cz9uTn6crN8drQ5fjY3x0fcs+4/9x+Px0vrrbWnxjyLQop168471+v14uPHjzPl7969W/nf5O378cfXbsHboS83R19ulv7cHH25Ofpys/Tn5ujLzdGXm6U/N0dfbo6+3Cz9uTn6crP05+boy83Rl4t9+fIlflzSQSuHGkmSzDxNMRqNIkmStequc52Tk5P417/+NTn+888/YzQaxd/+9rfY2dlZten8hW5vb+Pdu3dxdXUVe3t7r90cqBxjCJ7P+IHnM37gZYwheD7jB17GGOItG4/H8eXLl/j555+X1ls51Gi1WnF2djZT3mg01qqbpunK19nd3Y3d3d2psnnhB69vb2/PL1J4AWMIns/4geczfuBljCF4PuMHXsYY4q1a9oTGvZVDjTRNp46LoohGozEJGfI8jyRJIk3TpXUfhxKPrwMAAAAAADDPWntqDAaD6Ha70Ww24/z8PAaDweRcr9eLZrMZnU7nybrLzgEAAAAAAMyzVqiRpmmcnp5GRES73Z469ziYWFZ32TmqZXd3N/7973/PLBMGrMYYguczfuD5jB94GWMIns/4gZcxhiBiZzwej1+7EQAAAAAAAE/57rUbAAAAAAAAsAqhBgAAAAAAUAlr7anB9iqKIobDYaRpGkVRxPHxcSRJ8uK6sA3yPI8syyIi4vz8PD59+rRwTOR5HhER9Xo9iqKIsiyjXq//v5oK35x1xoT7D8waDofRarUiIp4cD+5B8HUcHB0dxcXFxVS5+RA8bdH4MR+Cpy0bPxHmQzBjDCuo1+uTry8vL8ftdnsjdWEbnJ6eTn39cIw8dnx8PI6IcUSMW63W+Pr6+v/QQvh2rTMm3H9g1v34efh6eF96yD2IbTcYDMYXFxfjedNk8yFYbtn4MR+C5ZaNH/MhmM+TGjypKIqp4zRNJ39l8ZK6sA3yPI9erxedTiciItrtdnS73SiKItI0nam/v78f19fXEfH0X9TCNlh1TLj/wKyyLGMwGES73Z6U9fv9yT3pMfcgtt3DsfKQ+RA8bdH4MR+Cpy0aPxHmQ7CIPTV4UpZlUavVpspqtdrkEbjn1oVtUK/X49OnT5PjsiwjImbGyUNJkngDDw+sMibcf2C+h5Pk4XC4dNIc4R4E85gPwfOZD8HLmQ/BLE9q8KT7Nx2PjUajF9WFbfHwA6Rff/01Wq3WwjckZVnGcDiMiK/rzX748GHuXzDBtlh1TLj/wKyH95qyLGM0Gi29p7gHwXzmQ/Ay5kPwfOZDMJ9Qg2db9AvzpXXhrbp/M/J446+HHm7klaZpvH//Pi4vL/9PLYRvz0vHhPsPfNXtduP09HRpHfcgWI/5EKzHfAjWZz4E81l+iiclSTKT7I5Go7l/WbFOXdg23W43Pn/+vPI6mGmaRlEUM2tjwjZZdUy4/8BiZVlGlmVPjgf3IJjPfAg2w3wI1mc+BPMJNXhSq9WaW95oNF5UF7ZJv9+PbrcbaZpGWZZz/1oiz/M4ODiYKV+23iy8ZeuMCfcfWOz3339/ckLrHgSLmQ/By5kPwfrMh2AxoQZPerxWX1EU0Wg0JpPjPM8nKfFTdWEbDYfDqNfrkzfwv/3228Lx83BpkCzLot1uGz9srafGhPsPrCbP87mTX/cgWOzhB67mQ7Cex4GF+RCs7vH9x3wI5tsZj8fj124E376iKOLs7CyazWacn5/HycnJ5Bfj4eFhNJvN6HQ6T9aFbVMURfzjH/+YKkuSJK6vryNidvzkeT5ZIuTy8vLJ9c/hrVs2Jtx/YDX9fj8uLy/j7Oxsqtw9CKZlWRafP3+Ofr8fnU4nms3mZINj8yFYbtH4MR+Cpy27/5gPwXxCDQAAAAAAoBIsPwUAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEv4H3stktjxfnN0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAC+CAYAAAB9JI0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATzklEQVR4nO3dT04jSdoH4JdRS6yayvbsWsMm5wbGPkCrXDcwR4Daz8IWq1GtLHwDuk7Qbd+gUuoD0OQNyA1SL00WtUIatb9FCX8Y/8EGd1OJn0eyRGaGo4JQBengR0bsjMfjcQAAAAAAAHzj/vHSDQAAAAAAAFiFUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAq4buXbsA6/vzzz/jjjz/i+++/j52dnZduDgAAAAAAsAHj8Ti+fPkSP/74Y/zjH4ufx6hUqPHHH3/E/v7+SzcDAAAAAAD4C1xdXcW//vWvhdcrFWp8//33EfH1m9rb23vh1gAAAAAAAJtwc3MT+/v7kxxgkUqFGndLTu3t7Qk1AAAAAADglXls6wkbhQMAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKiESm0UDrDNfvrppVvwcn77bfN1bmt//hV9CQBUn89GAEBVCDUAAAAANmRbA6IIIREAfw+hBvCX2tYP9D7Ms02McwAAAODvItQAAACgkgTr8PoZ5wA8ZKNwAAAAAACgEoQaAAAAAABAJQg1AAAAAACASlhrT408z+Po6CguLi6WliuKIobDYaRpGkVRxPHxcSRJ8ug1AAAAAACARVYONe6CiDzPHy17eHg4CT6Kooijo6MYDAaPXgMAAAAAAFhk5VCj3W6vVK4oiqnjNE0jy7JHrwEAAAAAACyz8T01siyLWq02da5Wq0We50uvAQAAAAAALLPWnhqrKMty7vnRaLT02jy3t7dxe3s7Ob65uXlu8wAA2BI//fTSLXgZv/22+Tr1JQDA/9vWz0YRPh/xbdh4qLHIokBj2bVerxcfPnz4axr0yvhhCgAAAABUzbb+XtPvNJ9u48tPJUky8+TFaDSKJEmWXpvn5OQkPn/+PHldXV1turkAAAAAAEBFbDzUaLVac883Go2l1+bZ3d2Nvb29qRcAAAAAALCdnhRqPFwuKs/zKIoiIiLSNJ26VhRFNBqNSJJk6TUAAAAAAIBlVt5TI8uy+PTpU0R83eui2WxGu92eOu50OhERMRgMotvtRrPZjPPz8xgMBpN6ll0DAAAAAABYZOVQo9VqRavVitPT05lrD4OJNE0n5e6Cj1WuAQAAAAAALLLxPTUAAAAAAAD+CkINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBK+e+kGAAAAbIuffnrpFryc33576RYAAPAaCDXgARNNAAAAAIBvk+WnAAAAAACASvCkBgDAN8LTggAAALCcUAMAAAAAXjl/QAO8FpafAgAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCWsFWoURRH9fj+Gw2H0+/0oy3Jh2eFwGGVZzi2T53nkeT6p8+5rAAAAAACARdYKNQ4PD6PT6US73Y52ux1HR0dLy/7www/xww8/xM7OTuzs7ES/34+IiLOzszg4OIidnZ14//59pGn6vO8CAAAAAAB49b5btWBRFFPHaZpGlmVzy5ZlGYPBINrt9uRcv9+PTqcTEREHBwdxfX0dERFJkqzbZgAAAAAAYAut/KRGlmVRq9WmztVqtYVLR90PNIbD4dRxxNcwQ6ABAAAAAACsauUnNRbtnzEajWbO3Q8ryrKM0Wg0tcRUWZYxHA4jIuL8/NwSVAAAAAAAwKNWDjUWWbZZeEREt9uN09PTqXPHx8eT4CNN03j37l1cXl7OvPf29jZub28nxzc3N89tLgAAAAAAUFErLz+VJMnMUxmj0WjpElJlWUaWZTNl7u/PkaZpFEUxs2dHRESv14s3b95MXvv7+6s2FwAAAAAAeGVWDjVardbc841GY+F7fv/995lAI8/zePv27UzZh/t1REScnJzE58+fJ6+rq6tVmwsAAAAAALwyKy8/9XDPi6IootFoTEKLPM8jSZKpcnmez4QVaZpOLUeVZVm02+25T3zs7u7G7u7uqk0EAAAAAABesbX21BgMBtHtdqPZbMb5+XkMBoPJtV6vF81mMzqdztR7HoYhSZJEo9GIfr8fSZLE5eXlVD0AAAAAAADzrBVq3H/Kot1uT12bF0w8DDju1Ov1qNfr6/zTAAAAAADAllt5Tw0AAAAAAICXJNQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJ361TuCiKGA6HkaZpFEURx8fHkSTJ3LJ5nkdERL1ej6IooizLqNfra9cDAAAAAAAQsWaocXh4GBcXFxHxNZg4OjqKwWAwt+zZ2Vn8/PPPERHRarWmyq1TDwAAAAAAQMQaoUZRFFPHaZpGlmULyx8cHMT19XVExNRTGOvWAwAAAAAAELHGnhpZlkWtVps6V6vVJstMzZMkycyyUk+pBwAAAAAAYOUnNcqynHt+NBotLD8cDiMi4vz8PN6/fx9pmq5dDwAAAAAAQMSae2rMsyikuL/5d5qm8e7du7i8vFyrntvb27i9vZ0c39zcPKepAAAAAABAha28/FSSJDNPU4xGo5nlpe7c3zsjTdMoiiKKolirnl6vF2/evJm89vf3V20uAAAAAADwyqwcarRarbnnG43GzLk8z+Pt27cz52u12lr1nJycxOfPnyevq6urVZsLAAAAAAC8MisvP5Wm6dRxURTRaDQmT1jkeR5JkkSappGmaZyenk7KZlkW7XZ77sbhD+u5b3d3N3Z3d1f/bgAAAAAAgFdrrT01BoNBdLvdaDabcX5+HoPBYHKt1+tFs9mMTqcTSZJEo9GIfr8fSZLE5eXlVNll9QAAAAAAAMyzVqhx/wmMdrs9de1hMFGv16Ner69dDwAAAAAAwDwr76kBAAAAAADwkoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUwnfrFC6KIobDYaRpGkVRxPHxcSRJMrdsnueRZVlERJyfn8fHjx8nZfM8j4iIer0eRVFEWZZRr9ef/l0AAAAAAACv3lqhxuHhYVxcXETE14Dj6OgoBoPB3LJZlkWn04mIiH6/H2/fvp289+zsLH7++eeIiGi1WgvrAAAAAAAAuLPy8lNFUUwdp2k6eRLjoTzPo9frTY7b7XbkeT6p4+DgIK6vr+P6+jo+ffq08GkPAAAAAACAOyuHGlmWRa1WmzpXq9UmS0ndV6/X4+PHj5Pjsiwn5e8kSSLMAAAAAAAAVrby8lN3wcRDo9Fo7vl2uz35+pdffolWqzUJMcqyjOFwGBFf99t4//59pGk6U8ft7W3c3t5Ojm9ublZtLgAAAAAA8MqstafGPIvCjvvXh8PhZD+NiJjaYDxN03j37l1cXl7OvLfX68WHDx+e20QAAAAAAOAVWHn5qSRJZp7KGI1Gjy4h1e12Z/bNuL8/R5qmURTFzJ4dEREnJyfx+fPnyevq6mrV5gIAAAAAAK/MyqFGq9Wae77RaCx8T7/fj263G2maRlmWUZZl5Hkeb9++nSn7cL+OiIjd3d3Y29ubegEAAAAAANtp5VDj4Z4XRVFEo9GYPIGR5/nU0xbD4TDq9fok0Pj1118jSZJI0zROT08n5bIsi3a7bdNwAAAAAABgqbX21BgMBtHtdqPZbMb5+XkMBoPJtV6vF81mMzqdThRFEYeHh1PvTZJkspdGo9GIfr8fSZLE5eXlVD0AAAAAAADzrBVq3H/Kot1uT127H0ykaRrj8XhhPfV6Per1+jr/NAAAAAAAsOVWXn4KAAAAAADgJQk1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUwnfrFC6KIobDYaRpGkVRxPHxcSRJsnbZdeoBAAAAAACIWDPUODw8jIuLi4j4GkwcHR3FYDBYu+w69QAAAAAAAESssfxUURRTx2maRpZla5ddpx4AAAAAAIA7K4caWZZFrVabOler1SLP87XKrlMPAAAAAADAnZWXnyrLcu750Wi0Vtl16rm9vY3b29vJ8efPnyMi4ubmZnljt9D//vfSLXg5m/7voC83a1v7U19ulv7cHH25Ofpys/Tn5ujLzdGXm6U/N0dfbo6+3Cz9uTn6crP87mhz/N/cHL/innX3e//xeLy03Fp7asyzKKRYt+y8a71eLz58+DBzfn9/f+V/k9fvzZuXbsHroS83R19ulv7cHH25Ofpys/Tn5ujLzdGXm6U/N0dfbo6+3Cz9uTn6crP05+boy83Rl4t9+fIl3izpoJVDjSRJZp6mGI1GkSTJWmXXqefk5CT+85//TI7//PPPGI1G8c9//jN2dnZWbTp/oZubm9jf34+rq6vY29t76eZA5RhD8HTGDzyd8QPPYwzB0xk/8DzGEK/ZeDyOL1++xI8//ri03MqhRqvVirOzs5nzjUZjrbJpmq5cz+7ubuzu7k6dmxd+8PL29vb8IIVnMIbg6YwfeDrjB57HGIKnM37geYwhXqtlT2jcWTnUSNN06rgoimg0GpOQIc/zSJIk0jRdWvZhKPGwHgAAAAAAgHnW2lNjMBhEt9uNZrMZ5+fnMRgMJtd6vV40m83odDqPll12DQAAAAAAYJ61Qo00TeP09DQiItrt9tS1h8HEsrLLrlEtu7u78d///ndmmTBgNcYQPJ3xA09n/MDzGEPwdMYPPI8xBBE74/F4/NKNAAAAAAAAeMw/XroBAAAAAAAAqxBqAAAAAAAAlbDWnhpsr6IoYjgcRpqmURRFHB8fR5Ikzy4L2yDP88iyLCIizs/P4+PHjwvHRJ7nERFRr9ejKIooyzLq9frf1VT45qwzJtx/YNZwOIxWqxUR8eh4cA+Cr+Pg6OgoLi4ups6bD8HjFo0f8yF43LLxE2E+BDPGsIJ6vT75+vLyctxutzdSFrbB6enp1Nf3x8hDx8fH44gYR8S41WqNr6+v/4YWwrdrnTHh/gOz7sbP/df9+9J97kFsu8FgML64uBjPmyabD8Fyy8aP+RAst2z8mA/BfJ7U4FFFUUwdp2k6+SuL55SFbZDnefR6veh0OhER0W63o9vtRlEUkabpTPmDg4O4vr6OiMf/oha2wapjwv0HZpVlGYPBINrt9uRcv9+f3JMecg9i290fK/eZD8HjFo0f8yF43KLxE2E+BIvYU4NHZVkWtVpt6lytVps8AvfUsrAN6vV6fPz4cXJclmVExMw4uS9JEh/g4Z5VxoT7D8x3f5I8HA6XTpoj3INgHvMheDrzIXg+8yGY5UkNHnX3oeOh0Wj0rLKwLe7/AumXX36JVqu18ANJWZYxHA4j4ut6s+/fv5/7F0ywLVYdE+4/MOv+vaYsyxiNRkvvKe5BMJ/5EDyP+RA8nfkQzCfU4MkW/cB8bll4re4+jDzc+Ou++xt5pWka7969i8vLy7+phfDtee6YcP+Br7rdbpyeni4t4x4E6zEfgvWYD8H6zIdgPstP8agkSWaS3dFoNPcvK9YpC9um2+3Gp0+fVl4HM03TKIpiZm1M2Carjgn3H1isLMvIsuzR8eAeBPOZD8FmmA/B+syHYD6hBo9qtVpzzzcajWeVhW3S7/ej2+1GmqZRluXcv5bI8zzevn07c37ZerPwmq0zJtx/YLHff//90QmtexAsZj4Ez2c+BOszH4LFhBo86uFafUVRRKPRmEyO8zyfpMSPlYVtNBwOo16vTz7A//rrrwvHz/2lQbIsi3a7bfywtR4bE+4/sJo8z+dOft2DYLH7v3A1H4L1PAwszIdgdQ/vP+ZDMN/OeDwev3Qj+PYVRRFnZ2fRbDbj/Pw8Tk5OJj8YDw8Po9lsRqfTebQsbJuiKOLf//731LkkSeL6+joiZsdPnueTJUIuLy8fXf8cXrtlY8L9B1bT7/fj8vIyzs7Ops67B8G0LMvi06dP0e/3o9PpRLPZnGxwbD4Eyy0aP+ZD8Lhl9x/zIZhPqAEAAAAAAFSC5acAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVML/AWyuYO2BIugMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute and plot spreads\n",
    "\n",
    "sigma_dqxs = dqxs.std(1)\n",
    "sigma_dqys = dqys.std(1)\n",
    "\n",
    "plt.figure(figsize=(16, 2))\n",
    "plt.bar(range(len(sigma_dqxs)), sigma_dqxs.cpu().numpy(), color='blue', alpha=0.75, width=0.75)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(16, 2))\n",
    "plt.bar(range(len(sigma_dqys)), sigma_dqys.cpu().numpy(), color='blue', alpha=0.75, width=0.75)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5f7ddba6-e8b3-40bc-811f-ce47f3a0e96c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([17, 28])\n",
      "torch.Size([17, 28])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAC+CAYAAAB9JI0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATtUlEQVR4nO3dT04jSdoH4JdWS6yazvbsWlObnBsY+wAjXDcwR4Daz8IWq1GtLHwDuk7Qbd+gUvoOQJM3IDdIvTQJtUIatb9FCQvjP9jgbirx80iWyMggiQpVkA7/yIid8Xg8DgAAAAAAgG/cd6/dAAAAAAAAgFUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKASvn/tBqzjzz//jD/++CN++OGH2NnZee3mAAAAAAAAGzAej+PLly/x888/x3ffLX4eo1Khxh9//BHv3r177WYAAAAAAAB/gaurq/jnP/+58HylQo0ffvghIr7+o/b29l65NQAAAAAAwCbc3t7Gu3fvJjnAIpUKNe6XnNrb2xNqAAAAAADAG/PU1hM2CgcAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJldooHAAAAO79+9+v3YLX8X//99otAAB4PUIN4C9loglvn3EOAAAA/F2EGgAAAAB8k/wBDQCP2VMDAAAAAACoBKEGAAAAAABQCZafAqiIbX3sOsKj1wAAf7Vtfa/5V7zP3Na+jPC+HYC/h1ADgK20rZNNE00AANhO2zoHijAPgrdGqAEA8I0w0dysbe1Pfbk5+nKzfKAEwFvhfg6vS6jxRvhlujn6EgAAAAD+Htv6WZzP4Z7PRuEAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVsNaeGnmex9HRUVxcXCytVxRFDIfDSNM0iqKI4+PjSJLkyXMAAAAAAACLrBxq3AcReZ4/Wffw8HASfBRFEUdHRzEYDJ48BwAAAAAAsMjKoUa73V6pXlEUU8dpmkaWZU+eAwAAAAAAWGbje2pkWRa1Wm2qrFarRZ7nS88BAAAAAAAss9aeGqsoy3Ju+Wg0Wnpunru7u7i7u5sc397evrR5AAAAAABARW38SY1FFgUay871er348ccfJ6937979NY0DAAAAAAC+eRsPNZIkmXnyYjQaRZIkS8/Nc3JyEjc3N5PX1dXVppsLAAAAAABUxMZDjVarNbe80WgsPTfP7u5u7O3tTb0AAAAAAIDt9KxQ4/FyUXmeR1EUERGRpunUuaIootFoRJIkS88BAAAAAAAss/JG4VmWxefPnyPi614XzWYz2u321HGn04mIiMFgEN1uN5rNZpyfn8dgMJhcZ9k5AAAAAACARVYONVqtVrRarTg9PZ059ziYSNN0Uu8++FjlHAAAAAAAwCIb31MDAAAAAADgryDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACohLVCjaIoot/vx3A4jH6/H2VZLqw7HA6jLMu5dfI8jzzPJ9e8/xoAAAAAAGCRtUKNw8PD6HQ60W63o91ux9HR0dK6P/30U/z000+xs7MTOzs70e/3IyLi7Ows9vf3Y2dnJz58+BBpmr7sXwEAAAAAALx5369asSiKqeM0TSPLsrl1y7KMwWAQ7XZ7Utbv96PT6URExP7+flxfX0dERJIk67YZAAAAAADYQis/qZFlWdRqtamyWq22cOmoh4HGcDicOo74GmYINAAAAAAAgFWt/KTGov0zRqPRTNnDsKIsyxiNRlNLTJVlGcPhMCIizs/PFy5BdXd3F3d3d5Pj29vbVZsLAAAAAAC8MSuHGoss2yw8IqLb7cbp6elU2fHx8ST4SNM03r9/H5eXlzPf2+v14uPHjy9tIgAAAAAA8AasvPxUkiQzT2WMRqOlS0iVZRlZls3Uebg/R5qmURTFzJ4dEREnJydxc3MzeV1dXa3aXAAAAAAA4I1ZOdRotVpzyxuNxsLv+f3332cCjTzP4+DgYKbu4/06IiJ2d3djb29v6gUAAAAAAGynlUONx3teFEURjUZjElrkeT7ztEWe5zNhRZqmU8tRZVkW7XbbpuEAAAAAAMBSa+2pMRgMotvtRrPZjPPz8xgMBpNzvV4vms1mdDqdqe95HIYkSRKNRiP6/X4kSRKXl5dT1wEAAAAAAJhnrVDj4VMW7XZ76ty8YOJxwHGvXq9HvV5f50cDAAAAAABbbuXlpwAAAAAAAF6TUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBK+X6dyURQxHA4jTdMoiiKOj48jSZK5dfM8j4iIer0eRVFEWZZRr9fXvg4AAAAAAEDEmqHG4eFhXFxcRMTXYOLo6CgGg8HcumdnZ/HLL79ERESr1Zqqt851AAAAAAAAItYINYqimDpO0zSyLFtYf39/P66vryMipp7CWPc6AAAAAAAAEWvsqZFlWdRqtamyWq02WWZqniRJZpaVes51AAAAAAAAVn5SoyzLueWj0Whh/eFwGBER5+fn8eHDh0jTdK3r3N3dxd3d3eT49vZ21eYCAAAAAABvzFp7asyzKKR4uPl3mqbx/v37uLy8XOs6vV4vPn78+NImAgAAAAAAb8DKy08lSTLzNMVoNJpZXurew70z0jSNoiiiKIq1rnNychI3NzeT19XV1arNBQAAAAAA3piVQ41WqzW3vNFozJTleR4HBwcz5bVaba3r7O7uxt7e3tQLAAAAAADYTisvP5Wm6dRxURTRaDQmT1jkeR5JkkSappGmaZyenk7qZlkW7XZ77sbhj68DAAAAAAAwz1p7agwGg+h2u9FsNuP8/DwGg8HkXK/Xi2azGZ1OJ5IkiUajEf1+P5IkicvLy6m6y64DAAAAAAAwz1qhxsMnMNrt9tS5x8FEvV6Per2+9nUAAAAAAADmWXlPDQAAAAAAgNck1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqITv16lcFEUMh8NI0zSKoojj4+NIkmRu3TzPI8uyiIg4Pz+PT58+TermeR4REfV6PYqiiLIso16vP/9fAQAAAAAAvHlrhRqHh4dxcXEREV8DjqOjoxgMBnPrZlkWnU4nIiL6/X4cHBxMvvfs7Cx++eWXiIhotVoLrwEAAAAAAHBv5eWniqKYOk7TdPIkxmN5nkev15sct9vtyPN8co39/f24vr6O6+vr+Pz588KnPQAAAAAAAO6tHGpkWRa1Wm2qrFarTZaSeqher8enT58mx2VZTurfS5JEmAEAAAAAAKxs5eWn7oOJx0aj0dzydrs9+frXX3+NVqs1CTHKsozhcBgRX/fb+PDhQ6RpOnONu7u7uLu7mxzf3t6u2lwAAAAAAOCNWWtPjXkWhR0Pzw+Hw8l+GhExtcF4mqbx/v37uLy8nPneXq8XHz9+fGkTAQAAAACAN2Dl5aeSJJl5KmM0Gj25hFS3253ZN+Ph/hxpmkZRFDN7dkREnJycxM3NzeR1dXW1anMBAAAAAIA3ZuVQo9VqzS1vNBoLv6ff70e32400TaMsyyjLMvI8j4ODg5m6j/friIjY3d2Nvb29qRcAAAAAALCdVg41Hu95URRFNBqNyRMYeZ5PPW0xHA6jXq9PAo3ffvstkiSJNE3j9PR0Ui/Lsmi32zYNBwAAAAAAllprT43BYBDdbjeazWacn5/HYDCYnOv1etFsNqPT6URRFHF4eDj1vUmSTPbSaDQa0e/3I0mSuLy8nLoOAAAAAADAPGuFGg+fsmi321PnHgYTaZrGeDxeeJ16vR71en2dHw0AAAAAAGy5lZefAgAAAAAAeE1CDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACASvh+ncpFUcRwOIw0TaMoijg+Po4kSdauu851AAAAAAAAItYMNQ4PD+Pi4iIivgYTR0dHMRgM1q67znUAAAAAAAAi1lh+qiiKqeM0TSPLsrXrrnMdAAAAAACAeyuHGlmWRa1Wmyqr1WqR5/ladde5DgAAAAAAwL2Vl58qy3Ju+Wg0WqvuOte5u7uLu7u7yfHNzU1ERNze3i5v7Bb63/9euwWvZ9P/HfTlZm1rf+rLzdKfm6MvN0dfbpb+3Bx9uTn6crP05+boy83Rl5ulPzdHX26Wz442x//NzfER96z7z/3H4/HSemvtqTHPopBi3brzzvV6vfj48eNM+bt371b+mbx9P/742i14O/Tl5ujLzdKfm6MvN0dfbpb+3Bx9uTn6crP05+boy83Rl5ulPzdHX26W/twcfbk5+nKxL1++xI9LOmjlUCNJkpmnKUajUSRJslbdda5zcnIS//nPfybHf/75Z4xGo/jHP/4ROzs7qzadv9Dt7W28e/curq6uYm9v77WbA5VjDMHzGT/wfMYPvIwxBM9n/MDLGEO8ZePxOL58+RI///zz0norhxqtVivOzs5myhuNxlp10zRd+Tq7u7uxu7s7VTYv/OD17e3t+UUKL2AMwfMZP/B8xg+8jDEEz2f8wMsYQ7xVy57QuLdyqJGm6dRxURTRaDQmIUOe55EkSaRpurTu41Di8XUAAAAAAADmWWtPjcFgEN1uN5rNZpyfn8dgMJic6/V60Ww2o9PpPFl32TkAAAAAAIB51go10jSN09PTiIhot9tT5x4HE8vqLjtHtezu7sZ///vfmWXCgNUYQ/B8xg88n/EDL2MMwfMZP/AyxhBE7IzH4/FrNwIAAAAAAOAp3712AwAAAAAAAFYh1AAAAAAAACphrT012F5FUcRwOIw0TaMoijg+Po4kSV5cF7ZBnueRZVlERJyfn8enT58Wjok8zyMiol6vR1EUUZZl1Ov1v6up8M1ZZ0y4/8Cs4XAYrVYrIuLJ8eAeBF/HwdHRUVxcXEyVmw/B0xaNH/MheNqy8RNhPgQzxrCCer0++fry8nLcbrc3Uhe2wenp6dTXD8fIY8fHx+OIGEfEuNVqja+vr/+GFsK3a50x4f4Ds+7Hz8PXw/vSQ+5BbLvBYDC+uLgYz5smmw/BcsvGj/kQLLds/JgPwXye1OBJRVFMHadpOvkri5fUhW2Q53n0er3odDoREdFut6Pb7UZRFJGm6Uz9/f39uL6+join/6IWtsGqY8L9B2aVZRmDwSDa7fakrN/vT+5Jj7kHse0ejpWHzIfgaYvGj/kQPG3R+IkwH4JF7KnBk7Isi1qtNlVWq9Umj8A9ty5sg3q9Hp8+fZocl2UZETEzTh5KksQbeHhglTHh/gPzPZwkD4fDpZPmCPcgmMd8CJ7PfAheznwIZnlSgyfdv+l4bDQavagubIuHHyD9+uuv0Wq1Fr4hKcsyhsNhRHxdb/bDhw9z/4IJtsWqY8L9B2Y9vNeUZRmj0WjpPcU9COYzH4KXMR+C5zMfgvmEGjzbol+YL60Lb9X9m5HHG3899HAjrzRN4/3793F5efk3tRC+PS8dE+4/8FW3243T09OlddyDYD3mQ7Ae8yFYn/kQzGf5KZ6UJMlMsjsajeb+ZcU6dWHbdLvd+Pz588rrYKZpGkVRzKyNCdtk1THh/gOLlWUZWZY9OR7cg2A+8yHYDPMhWJ/5EMwn1OBJrVZrbnmj0XhRXdgm/X4/ut1upGkaZVnO/WuJPM/j4OBgpnzZerPwlq0zJtx/YLHff//9yQmtexAsZj4EL2c+BOszH4LFhBo86fFafUVRRKPRmEyO8zyfpMRP1YVtNBwOo16vT97A//bbbwvHz8OlQbIsi3a7bfywtZ4aE+4/sJo8z+dOft2DYLGHH7iaD8F6HgcW5kOwusf3H/MhmG9nPB6PX7sRfPuKooizs7NoNptxfn4eJycnk1+Mh4eH0Ww2o9PpPFkXtk1RFPGvf/1rqixJkri+vo6I2fGT5/lkiZDLy8sn1z+Ht27ZmHD/gdX0+/24vLyMs7OzqXL3IJiWZVl8/vw5+v1+dDqdaDabkw2OzYdguUXjx3wInrbs/mM+BPMJNQAAAAAAgEqw/BQAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACASvh/lFxktrPDS0gAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAC+CAYAAAB9JI0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATyklEQVR4nO3dT04bSf8H4C+jkVgN6fG7G71s+r2BsQ8wwrmBOQJk/y5ssXqVlYVvwOQEM/YN0tIcgNA3oDdIszQdskL6afxbRFgY/8EGZ0jD80iW0tXlpiilaBcfumprPB6PAwAAAAAA4Dv3w3M3AAAAAAAAYBVCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKiEH5+7Aev4+++/46+//oqffvoptra2nrs5AAAAAADABozH4/jy5Uv88ssv8cMPi5/HqFSo8ddff8Xu7u5zNwMAAAAAAPgGLi8v49///vfC85UKNX766aeI+PpN7ezsPHNrAAAAAACATbi+vo7d3d1JDrBIpUKN2yWndnZ2hBoAAAAAAPDCPLT1hI3CAQAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVEKlNgoHgE359dfnbsHz+PPP524BAPA98tkIAKgKoQZARbzWiWaEySYAANXhczsAfFtCDeCbeq0f6H2Y5zUxzgEAAIB/ilADAACAShKsw8tnnANwn43CAQAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqIQf16mc53kcHh7G+fn50npFUcRwOIw0TaMoijg6OookSR48BwAAAAAAsMjKocZtEJHn+YN1Dw4OJsFHURRxeHgYg8HgwXMAAAAAAACLrBxqtNvtleoVRTF1nKZpZFn24DkAAAAAgO/dr78+dwuez59/PncLYM3lp1aRZVnUarWpslqtFnmex6dPnxaeq9frm27Kq+KHKQDAtNf6+ehbfDbSlwAAfCs+a7KujYcaZVnOLR+NRkvPzXNzcxM3NzeT4+vr66c2DwAAAAAAqKgf/qkvtCjQWHau1+vFmzdvJq/d3d1v0zgAAAAAAOC7t/FQI0mSmScvRqNRJEmy9Nw8x8fH8fnz58nr8vJy080FAAAAAAAqYuOhRqvVmlveaDSWnptne3s7dnZ2pl4AAAAAAMDr9KhQ4/5yUXmeR1EUERGRpunUuaIootFoRJIkS88BAAAAAAAss/JG4VmWxcePHyPi614XzWYz2u321HGn04mIiMFgEN1uN5rNZpydncVgMJhcZ9k5AAAAAACARVYONVqtVrRarTg5OZk5dz+YSNN0Uu82+FjlHAAAAAAAwCIb31MDAAAAAADgWxBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEr48bkbAADAV7/++twteD5//vncLQAAAKAKhBoAAAD/EOElAAA8jVADAAAAAF44wTrwUgg14B43eQAAAACA75ONwgEAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlrBVqFEUR/X4/hsNh9Pv9KMtyYd3hcBhlWc6tk+d55Hk+uebtvwEAAAAAABZZK9Q4ODiITqcT7XY72u12HB4eLq37888/x88//xxbW1uxtbUV/X4/IiJOT09jb28vtra24t27d5Gm6dO+CwAAAAAA4MX7cdWKRVFMHadpGlmWza1blmUMBoNot9uTsn6/H51OJyIi9vb24urqKiIikiRZt80AAAAAAMArtPKTGlmWRa1Wmyqr1WoLl466G2gMh8Op44ivYYZAAwAAAAAAWNXKT2os2j9jNBrNlN0NK8qyjNFoNLXEVFmWMRwOIyLi7Oxs4RJUNzc3cXNzMzm+vr5etbkAAAAAAMALs3KosciyzcIjIrrdbpycnEyVHR0dTYKPNE3j7du3cXFxMfPeXq8X79+/f2oTAQAAAACAF2Dl5aeSJJl5KmM0Gi1dQqosy8iybKbO3f050jSNoihm9uyIiDg+Po7Pnz9PXpeXl6s2FwAAAAAAeGFWDjVardbc8kajsfA9nz59mgk08jyP/f39mbr39+uIiNje3o6dnZ2pFwAAAAAA8DqtHGrc3/OiKIpoNBqT0CLP85mnLfI8nwkr0jSdWo4qy7Jot9s2DQcAAAAAAJZaa0+NwWAQ3W43ms1mnJ2dxWAwmJzr9XrRbDaj0+lMved+GJIkSTQajej3+5EkSVxcXExdBwAAAAAAYJ61Qo27T1m02+2pc/OCifsBx616vR71en2dLw0AAAAAALxyKy8/BQAAAAAA8JyEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJUg1AAAAAAAACpBqAEAAAAAAFSCUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlfDjOpWLoojhcBhpmkZRFHF0dBRJksytm+d5RETU6/UoiiLKsox6vb72dQAAAAAAACLWDDUODg7i/Pw8Ir4GE4eHhzEYDObWPT09jd9++y0iIlqt1lS9da4DAAAAAAAQsUaoURTF1HGappFl2cL6e3t7cXV1FREx9RTGutcBAAAAAACIWGNPjSzLolarTZXVarXJMlPzJEkys6zUY64DAAAAAACw8pMaZVnOLR+NRgvrD4fDiIg4OzuLd+/eRZqma18HAAAAAAAgYs09NeZZFFLc3fw7TdN4+/ZtXFxcrHWdm5ubuLm5mRxfX18/pakAAAAAAECFrbz8VJIkM09TjEajmeWlbt3dOyNN0yiKIoqiWOs6vV4v3rx5M3nt7u6u2lwAAAAAAOCFWTnUaLVac8sbjcZMWZ7nsb+/P1Neq9XWus7x8XF8/vx58rq8vFy1uQAAAAAAwAuz8vJTaZpOHRdFEY1GY/KERZ7nkSRJpGkaaZrGycnJpG6WZdFut+duHH7/Ondtb2/H9vb26t8NAAAAAADwYq21p8ZgMIhutxvNZjPOzs5iMBhMzvV6vWg2m9HpdCJJkmg0GtHv9yNJkri4uJiqu+w6AAAAAAAA86wVatx9AqPdbk+dux9M1Ov1qNfra18HAAAAAABgnpX31AAAAAAAAHhOQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCUINQAAAAAAgEoQagAAAAAAAJXw4zqVi6KI4XAYaZpGURRxdHQUSZLMrZvneWRZFhERZ2dn8eHDh0ndPM8jIqJer0dRFFGWZdTr9cd/FwAAAAAAwIu3VqhxcHAQ5+fnEfE14Dg8PIzBYDC3bpZl0el0IiKi3+/H/v7+5L2np6fx22+/RUREq9VaeA0AAAAAAIBbKy8/VRTF1HGappMnMe7L8zx6vd7kuN1uR57nk2vs7e3F1dVVXF1dxcePHxc+7QEAAAAAAHBr5VAjy7Ko1WpTZbVabbKU1F31ej0+fPgwOS7LclL/VpIkwgwAAAAAAGBlKy8/dRtM3DcajeaWt9vtyb9///33aLVakxCjLMsYDocR8XW/jXfv3kWapjPXuLm5iZubm8nx9fX1qs0FAAAAAABemLX21JhnUdhx9/xwOJzspxERUxuMp2kab9++jYuLi5n39nq9eP/+/VObCAAAAAAAvAArLz+VJMnMUxmj0ejBJaS63e7Mvhl39+dI0zSKopjZsyMi4vj4OD5//jx5XV5ertpcAAAAAADghVk51Gi1WnPLG43Gwvf0+/3odruRpmmUZRllWUae57G/vz9T9/5+HRER29vbsbOzM/UCAAAAAABep5VDjft7XhRFEY1GY/IERp7nU09bDIfDqNfrk0Djjz/+iCRJIk3TODk5mdTLsiza7bZNwwEAAAAAgKXW2lNjMBhEt9uNZrMZZ2dnMRgMJud6vV40m83odDpRFEUcHBxMvTdJksleGo1GI/r9fiRJEhcXF1PXAQAAAAAAmGetUOPuUxbtdnvq3N1gIk3TGI/HC69Tr9ejXq+v86UBAAAAAIBXbuXlpwAAAAAAAJ6TUAMAAAAAAKgEoQYAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACAShBqAAAAAAAAlSDUAAAAAAAAKkGoAQAAAAAAVIJQAwAAAAAAqAShBgAAAAAAUAlCDQAAAAAAoBKEGgAAAAAAQCX8uE7loihiOBxGmqZRFEUcHR1FkiRr113nOgAAAAAAABFrhhoHBwdxfn4eEV+DicPDwxgMBmvXXec6AAAAAAAAEWssP1UUxdRxmqaRZdnadde5DgAAAAAAwK2VQ40sy6JWq02V1Wq1yPN8rbrrXAcAAAAAAODWystPlWU5t3w0Gq1Vd53r3NzcxM3NzeT48+fPERFxfX29vLGv0P/933O34Pls+r+Dvtys19qf+nKz9Ofm6MvN0ZebpT83R19ujr7cLP25Ofpyc/TlZunPzdGXm+V3R5vj/+bm+BX3rNvf+4/H46X11tpTY55FIcW6deed6/V68f79+5ny3d3dlb8mL9+bN8/dgpdDX26Ovtws/bk5+nJz9OVm6c/N0Zeboy83S39ujr7cHH25Wfpzc/TlZunPzdGXm6MvF/vy5Uu8WdJBK4caSZLMPE0xGo0iSZK16q5znePj4/jvf/87Of77779jNBrFv/71r9ja2lq16XxD19fXsbu7G5eXl7Gzs/PczYHKMYbg8YwfeDzjB57GGILHM37gaYwhXrLxeBxfvnyJX375ZWm9lUONVqsVp6enM+WNRmOtummarnyd7e3t2N7eniqbF37w/HZ2dvwghScwhuDxjB94POMHnsYYgsczfuBpjCFeqmVPaNxaOdRI03TquCiKaDQak5Ahz/NIkiTSNF1a934ocf86AAAAAAAA86y1p8ZgMIhutxvNZjPOzs5iMBhMzvV6vWg2m9HpdB6su+wcAAAAAADAPGuFGmmaxsnJSUREtNvtqXP3g4lldZedo1q2t7fjf//738wyYcBqjCF4POMHHs/4gacxhuDxjB94GmMIIrbG4/H4uRsBAAAAAADwkB+euwEAAAAAAACrEGoAAAAAAACVsNaeGrxeRVHEcDiMNE2jKIo4OjqKJEmeXBdegzzPI8uyiIg4OzuLDx8+LBwTeZ5HRES9Xo+iKKIsy6jX6/9UU+G7s86YcP+BWcPhMFqtVkTEg+PBPQi+joPDw8M4Pz+fKjcfgoctGj/mQ/CwZeMnwnwIZoxhBfV6ffLvi4uLcbvd3khdeA1OTk6m/n13jNx3dHQ0johxRIxbrdb46urqH2ghfL/WGRPuPzDrdvzcfd29L93lHsRrNxgMxufn5+N502TzIVhu2fgxH4Lllo0f8yGYz5MaPKgoiqnjNE0nf2XxlLrwGuR5Hr1eLzqdTkREtNvt6Ha7URRFpGk6U39vby+urq4i4uG/qIXXYNUx4f4Ds8qyjMFgEO12e1LW7/cn96T73IN47e6OlbvMh+Bhi8aP+RA8bNH4iTAfgkXsqcGDsiyLWq02VVar1SaPwD22LrwG9Xo9Pnz4MDkuyzIiYmac3JUkiQ/wcMcqY8L9B+a7O0keDodLJ80R7kEwj/kQPJ75EDyd+RDM8qQGD7r90HHfaDR6Ul14Le7+Aun333+PVqu18ANJWZYxHA4j4ut6s+/evZv7F0zwWqw6Jtx/YNbde01ZljEajZbeU9yDYD7zIXga8yF4PPMhmE+owaMt+oH51LrwUt1+GLm/8ddddzfyStM03r59GxcXF/9QC+H789Qx4f4DX3W73Tg5OVlaxz0I1mM+BOsxH4L1mQ/BfJaf4kFJkswku6PRaO5fVqxTF16bbrcbHz9+XHkdzDRNoyiKmbUx4TVZdUy4/8BiZVlGlmUPjgf3IJjPfAg2w3wI1mc+BPMJNXhQq9WaW95oNJ5UF16Tfr8f3W430jSNsizn/rVEnuexv78/U75svVl4ydYZE+4/sNinT58enNC6B8Fi5kPwdOZDsD7zIVhMqMGD7q/VVxRFNBqNyeQ4z/NJSvxQXXiNhsNh1Ov1yQf4P/74Y+H4ubs0SJZl0W63jR9erYfGhPsPrCbP87mTX/cgWOzuL1zNh2A99wML8yFY3f37j/kQzLc1Ho/Hz90Ivn9FUcTp6Wk0m804OzuL4+PjyQ/Gg4ODaDab0el0HqwLr01RFPGf//xnqixJkri6uoqI2fGT5/lkiZCLi4sH1z+Hl27ZmHD/gdX0+/24uLiI09PTqXL3IJiWZVl8/Pgx+v1+dDqdaDabkw2OzYdguUXjx3wIHrbs/mM+BPMJNQAAAAAAgEqw/BQAAAAAAFAJQg0AAAAAAKAShBoAAAAAAEAlCDUAAAAAAIBKEGoAAAAAAACVINQAAAAAAAAqQagBAAAAAABUglADAAAAAACoBKEGAAAAAABQCUINAAAAAACgEoQaAAAAAABAJQg1AAAAAACASvh/sTNkSGC47c4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute twiss derivatives and estimate spread from linear surrogate model using MC\n",
    "\n",
    "# Compute derivatives\n",
    "\n",
    "kn = torch.zeros(nq, dtype=torch.float64)\n",
    "_, dbqx_dk, _, dbqy_dk = torch.func.jacrev(lambda kn: twiss(ring, [kn], ('kn', ['Quadrupole'], None, None), matched=True, advance=True, limit=1, epsilon=None), chunk_size=1024)(kn).swapaxes(0, 1)\n",
    "\n",
    "print(dbqx_dk.shape)\n",
    "print(dbqy_dk.shape)\n",
    "\n",
    "# Sample\n",
    "\n",
    "kns = 0.01*torch.randn((8192, nq), dtype=torch.float64)\n",
    "\n",
    "dqxs = 100*((dbqx_dk @ kns.T).T/bqx).T\n",
    "dqys = 100*((dbqy_dk @ kns.T).T/bqy).T\n",
    "\n",
    "# Compute and plot spreads\n",
    "\n",
    "sigma_dqxs = dqxs.std(1)\n",
    "sigma_dqys = dqys.std(1)\n",
    "\n",
    "plt.figure(figsize=(16, 2))\n",
    "plt.bar(range(len(sigma_dqxs)), sigma_dqxs.cpu().numpy(), color='blue', alpha=0.75, width=0.75)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(16, 2))\n",
    "plt.bar(range(len(sigma_dqys)), sigma_dqys.cpu().numpy(), color='blue', alpha=0.75, width=0.75)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "59e9c8e2-a7cc-4ec3-ad7b-905acdb5d006",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAC+CAYAAAB9JI0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATsUlEQVR4nO3dT04jSdoH4JdWS6yayvbsWsMm5wYu+wCjct3AHAFqPwtbrEa1svAN6DpBt32DSuk7AE3egNwg9dJkUSukUftblLAw/oMT3E0nPI9kicwMJ1GhCsLBj4zYmU6n0wAAAAAAAPib++65KwAAAAAAALAJoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUwvfPXYEq/vjjj/j999/jhx9+iJ2dneeuDgAAAAAAsAXT6TS+fv0aP/30U3z33ernMWoVavz++++xv7//3NUAAAAAAAD+BJeXl/HPf/5z5fVahRo//PBDRHz7R+3t7T1zbQAAAAAAgG24vr6O/f39WQ6wSq1Cjdslp/b29oQaAAAAAADwwjy09YSNwgEAAAAAgFoQagAAAAAAALUg1AAAAAAAAGpBqAEAAAAAANRCrTYKB+rn3/9+7ho8j//7v+euAfx19HMAnosxCF4+/RyA+zypAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBZsFA5QE691g7yIP2eTvNfanjYcBACW8dloe15rW0b4rAnAX0Oo8UL40LQ92hKA52IM2q7X2p7acnu05Xb5rAnwvIxB26Mtt+u1tqfPRo9n+SkAAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1UGmj8DzP4/DwMM7Pz9eWK4oixuNxpGkaRVHE0dFRJEny4DUAAAAAAIBVNg41boOIPM8fLHtwcDALPoqiiMPDwxiNRg9eAwAAAAAAWGXjUKPb7W5UriiKueM0TSPLsgevAQAAAAAArLP1PTWyLItGozF3rtFoRJ7na68BAAAAAACsU2lPjU2UZbn0/GQyWXttmZubm7i5uZkdX19fP7V6AAAAAABATW39SY1VVgUa664NBoN48+bN7LW/v//nVA4AAAAAAPjb23qokSTJwpMXk8kkkiRZe22Z4+Pj+PLly+x1eXm57eoCAAAAAAA1sfVQo9PpLD3farXWXltmd3c39vb25l4AAAAAAMDr9KhQ4/5yUXmeR1EUERGRpunctaIootVqRZIka68BAAAAAACss/FG4VmWxefPnyPi214X7XY7ut3u3HGv14uIiNFoFP1+P9rtdpydncVoNJrdZ901AAAAAACAVTYONTqdTnQ6nTg5OVm4dj+YSNN0Vu42+NjkGgAAAAAAwCpb31MDAAAAAADgzyDUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBaEGgAAAAAAQC0INQAAAAAAgFoQagAAAAAAALUg1AAAAAAAAGpBqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACohUqhRlEUMRwOYzwex3A4jLIsV5Ydj8dRluXSMnmeR57ns3vefg0AAAAAALBKpVDj4OAger1edLvd6Ha7cXh4uLbsjz/+GD/++GPs7OzEzs5ODIfDiIg4PT2Nt2/fxs7OTnz48CHSNH3avwIAAAAAAHjxvt+0YFEUc8dpmkaWZUvLlmUZo9Eout3u7NxwOIxerxcREW/fvo2rq6uIiEiSpGqdAQAAAACAV2jjJzWyLItGozF3rtForFw66m6gMR6P544jvoUZAg0AAAAAAGBTGz+psWr/jMlksnDublhRlmVMJpO5JabKsozxeBwREWdnZ5agAgAAAAAAHrRxqLHKus3CIyL6/X6cnJzMnTs6OpoFH2maxvv37+Pi4mLhvTc3N3FzczM7vr6+fmp1AQAAAACAmtp4+akkSRaeyphMJmuXkCrLMrIsWyhzd3+ONE2jKIqFPTsiIgaDQbx582b22t/f37S6AAAAAADAC7NxqNHpdJaeb7VaK9/z22+/LQQaeZ7Hu3fvFsre368jIuL4+Di+fPkye11eXm5aXQAAAAAA4IXZePmp+3teFEURrVZrFlrkeR5JksyVy/N8IaxI03RuOaosy6Lb7S594mN3dzd2d3c3rSIAAAAAAPCCVdpTYzQaRb/fj3a7HWdnZzEajWbXBoNBtNvt6PV6c++5H4YkSRKtViuGw2EkSRIXFxdz9wEAAAAAAFimUqhx9ymLbrc7d21ZMHE/4LjVbDaj2WxW+dYAAAAAAMArt/GeGgAAAAAAAM9JqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACoBaEGAAAAAABQC0INAAAAAACgFoQaAAAAAABALQg1AAAAAACAWhBqAAAAAAAAtSDUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBa+r1K4KIoYj8eRpmkURRFHR0eRJMnSsnmeR0REs9mMoiiiLMtoNpuV7wMAAAAAABBRMdQ4ODiI8/PziPgWTBweHsZoNFpa9vT0NH7++eeIiOh0OnPlqtwHAAAAAAAgokKoURTF3HGappFl2cryb9++jaurq4iIuacwqt4HAAAAAAAgosKeGlmWRaPRmDvXaDRmy0wtkyTJwrJSj7kPAAAAAADAxk9qlGW59PxkMllZfjweR0TE2dlZfPjwIdI0rXSfm5ubuLm5mR1fX19vWl0AAAAAAOCFqbSnxjKrQoq7m3+naRrv37+Pi4uLSvcZDAbx8ePHp1YRAAAAAAB4ATZefipJkoWnKSaTycLyUrfu7p2RpmkURRFFUVS6z/HxcXz58mX2ury83LS6AAAAAADAC7NxqNHpdJaeb7VaC+fyPI93794tnG80GpXus7u7G3t7e3MvAAAAAADgddp4+ak0TeeOi6KIVqs1e8Iiz/NIkiTSNI00TePk5GRWNsuy6Ha7SzcOv38fAAAAAACAZSrtqTEajaLf70e73Y6zs7MYjUaza4PBINrtdvR6vUiSJFqtVgyHw0iSJC4uLubKrrsPAAAAAADAMpVCjbtPYHS73blr94OJZrMZzWaz8n0AAAAAAACW2XhPDQAAAAAAgOck1AAAAAAAAGpBqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACoBaEGAAAAAABQC0INAAAAAACgFoQaAAAAAABALQg1AAAAAACAWhBqAAAAAAAAtSDUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqIXvqxQuiiLG43GkaRpFUcTR0VEkSbK0bJ7nkWVZREScnZ3Fp0+fZmXzPI+IiGazGUVRRFmW0Ww2H/+vAAAAAAAAXrxKocbBwUGcn59HxLeA4/DwMEaj0dKyWZZFr9eLiIjhcBjv3r2bvff09DR+/vnniIjodDor7wEAAAAAAHBr4+WniqKYO07TdPYkxn15nsdgMJgdd7vdyPN8do+3b9/G1dVVXF1dxefPn1c+7QEAAAAAAHBr41Ajy7JoNBpz5xqNxmwpqbuazWZ8+vRpdlyW5az8rSRJhBkAAAAAAMDGNl5+6jaYuG8ymSw93+12Z1//8ssv0el0ZiFGWZYxHo8j4tt+Gx8+fIg0TRfucXNzEzc3N7Pj6+vrTasLAAAAAAC8MJX21FhmVdhx9/p4PJ7tpxERcxuMp2ka79+/j4uLi4X3DgaD+Pjx41OrCAAAAAAAvAAbLz+VJMnCUxmTyeTBJaT6/f7Cvhl39+dI0zSKoljYsyMi4vj4OL58+TJ7XV5eblpdAAAAAADghdk41Oh0OkvPt1qtle8ZDofR7/cjTdMoyzLKsow8z+Pdu3cLZe/v1xERsbu7G3t7e3MvAAAAAADgddo41Li/50VRFNFqtWZPYOR5Pve0xXg8jmazOQs0fv3110iSJNI0jZOTk1m5LMui2+3aNBwAAAAAAFir0p4ao9Eo+v1+tNvtODs7i9FoNLs2GAyi3W5Hr9eLoiji4OBg7r1Jksz20mi1WjEcDiNJkri4uJi7DwAAAAAAwDKVQo27T1l0u925a3eDiTRNYzqdrrxPs9mMZrNZ5VsDAAAAAACv3MbLTwEAAAAAADwnoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACoBaEGAAAAAABQC0INAAAAAACgFoQaAAAAAABALQg1AAAAAACAWhBqAAAAAAAAtSDUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBaEGgAAAAAAQC18X6VwURQxHo8jTdMoiiKOjo4iSZLKZavcBwAAAAAAIKJiqHFwcBDn5+cR8S2YODw8jNFoVLlslfsAAAAAAABEVFh+qiiKueM0TSPLssplq9wHAAAAAADg1sahRpZl0Wg05s41Go3I87xS2Sr3AQAAAAAAuLXx8lNlWS49P5lMKpWtcp+bm5u4ubmZHX/58iUiIq6vr9dX9hX63/+euwbPZ9v/HbTldr3W9tSW26U9t0dbbo+23C7tuT3acnu05XZpz+3RltujLbdLe26PttwuvzvaHv83t8evuBfd/t5/Op2uLVdpT41lVoUUVcsuuzYYDOLjx48L5/f39zf+nrx8b948dw1eDm25Pdpyu7Tn9mjL7dGW26U9t0dbbo+23C7tuT3acnu05XZpz+3RltulPbdHW26Ptlzt69ev8WZNA20caiRJsvA0xWQyiSRJKpWtcp/j4+P4z3/+Mzv+448/YjKZxD/+8Y/Y2dnZtOr8ia6vr2N/fz8uLy9jb2/vuasDtaMPwePpP/B4+g88jT4Ej6f/wNPoQ7xk0+k0vn79Gj/99NPachuHGp1OJ05PTxfOt1qtSmXTNN34Pru7u7G7uzt3bln4wfPb29vzgxSeQB+Cx9N/4PH0H3gafQgeT/+Bp9GHeKnWPaFxa+NQI03TueOiKKLVas1ChjzPI0mSSNN0bdn7ocT9+wAAAAAAACxTaU+N0WgU/X4/2u12nJ2dxWg0ml0bDAbRbrej1+s9WHbdNQAAAAAAgGUqhRppmsbJyUlERHS73blr94OJdWXXXaNednd347///e/CMmHAZvQheDz9Bx5P/4Gn0Yfg8fQfeBp9CCJ2ptPp9LkrAQAAAAAA8JDvnrsCAAAAAAAAmxBqAAAAAAAAtVBpTw1er6IoYjweR5qmURRFHB0dRZIkTy4Lr0Ge55FlWUREnJ2dxadPn1b2iTzPIyKi2WxGURRRlmU0m82/qqrwt1OlTxh/YNF4PI5OpxMR8WB/MAbBt35weHgY5+fnc+fNh+Bhq/qP+RA8bF3/iTAfggVT2ECz2Zx9fXFxMe12u1spC6/BycnJ3Nd3+8h9R0dH04iYRsS00+lMr66u/oIawt9XlT5h/IFFt/3n7uvuuHSXMYjXbjQaTc/Pz6fLpsnmQ7Deuv5jPgTrres/5kOwnCc1eFBRFHPHaZrO/sriKWXhNcjzPAaDQfR6vYiI6Ha70e/3oyiKSNN0ofzbt2/j6uoqIh7+i1p4DTbtE8YfWFSWZYxGo+h2u7Nzw+FwNibdZwzitbvbV+4yH4KHreo/5kPwsFX9J8J8CFaxpwYPyrIsGo3G3LlGozF7BO6xZeE1aDab8enTp9lxWZYREQv95K4kSXyAhzs26RPGH1ju7iR5PB6vnTRHGINgGfMheDzzIXg68yFY5EkNHnT7oeO+yWTypLLwWtz9BdIvv/wSnU5n5QeSsixjPB5HxLf1Zj98+LD0L5jgtdi0Txh/YNHdsaYsy5hMJmvHFGMQLGc+BE9jPgSPZz4Eywk1eLRVPzCfWhZeqtsPI/c3/rrr7kZeaZrG+/fv4+Li4i+qIfz9PLVPGH/gm36/HycnJ2vLGIOgGvMhqMZ8CKozH4LlLD/Fg5IkWUh2J5PJ0r+sqFIWXpt+vx+fP3/eeB3MNE2jKIqFtTHhNdm0Txh/YLWyLCPLsgf7gzEIljMfgu0wH4LqzIdgOaEGD+p0OkvPt1qtJ5WF12Q4HEa/3480TaMsy6V/LZHnebx7927h/Lr1ZuElq9InjD+w2m+//fbghNYYBKuZD8HTmQ9BdeZDsJpQgwfdX6uvKIpotVqzyXGe57OU+KGy8BqNx+NoNpuzD/C//vrryv5zd2mQLMui2+3qP7xaD/UJ4w9sJs/zpZNfYxCsdvcXruZDUM39wMJ8CDZ3f/wxH4LldqbT6fS5K8HfX1EUcXp6Gu12O87OzuL4+Hj2g/Hg4CDa7Xb0er0Hy8JrUxRF/Otf/5o7lyRJXF1dRcRi/8nzfLZEyMXFxYPrn8NLt65PGH9gM8PhMC4uLuL09HTuvDEI5mVZFp8/f47hcBi9Xi/a7fZsg2PzIVhvVf8xH4KHrRt/zIdgOaEGAAAAAABQC5afAgAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBaEGgAAAAAAQC0INQAAAAAAgFoQagAAAAAAALUg1AAAAAAAAGpBqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFAL/w9Rm2F0A6GmlgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAC+CAYAAAB9JI0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAATxElEQVR4nO3dT04bSf8H4C+jkVgN0+N3N3qz6fcGjn2AUZwbmCNA9u/CFqtXWVn4BkxOMGPfIC3NAQh9A3qDNEvTISukn8a/RYSF8R/c4AzT8DySJbq7KCqlFO3yh67amU6n0wAAAAAAAPiH++6pGwAAAAAAALAJoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUwvdP3YAq/vrrr/jzzz/jhx9+iJ2dnaduDgAAAAAAsAXT6TS+fPkSP//8c3z33ernMWoVavz555/x6tWrp24GAAAAAADwDVxcXMS///3vlddrFWr88MMPEfH1H7W3t/fErQEAAAAAALbh6uoqXr16NcsBVqlVqHGz5NTe3p5QAwAAAAAAnpn7tp6wUTgAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaqNVG4QAAAMD2/fLLU7fgafzxx/brfKl9GfFt+hMA7hJqANSEydF2vdT+NNEEAAAA6kyoAQA8ioAIAIBvxXtNAO4SagDflDegAAB8K95rAgC8PDYKBwAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBa+r1I4z/M4ODiIs7OzteWKoojxeBxpmkZRFHF4eBhJktx7DQAAAAAAYJWNQ42bICLP83vL7u/vz4KPoiji4OAgRqPRvdcAAAAAAABW2TjU6Ha7G5UrimLuOE3TyLLs3msAAAAAAP90v/zy1C14On/88dQtgIrLT20iy7JoNBpz5xqNRuR5Hp8+fVp5rdlsbrspL4pfpgAA817q+6Nv8d5IXwIA8K14r0lVWw81yrJcen4ymay9tsz19XVcX1/Pjq+urh7bPAAAAAAAoKa++7t+0KpAY921wWAQP/744+z16tWrb9M4AAAAAADgH2/roUaSJAtPXkwmk0iSZO21ZY6OjuLz58+z18XFxbabCwAAAAAA1MTWQ41Op7P0fKvVWnttmd3d3djb25t7AQAAAAAAL9ODQo27y0XleR5FUURERJqmc9eKoohWqxVJkqy9BgAAAAAAsM7GG4VnWRYfP36MiK97XbTb7eh2u3PHvV4vIiJGo1H0+/1ot9txenoao9FoVs+6awAAAAAAAKtsHGp0Op3odDpxfHy8cO1uMJGm6azcTfCxyTUAAAAAAIBVtr6nBgAAAAAAwLcg1AAAAAAAAGpBqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC18P1TNwAAAOCl+OWXp27B0/njj6duAQAAz4FQAwDgH8KHnQAAALCeUAPu8IESAAAA8Nz4vAN4LuypAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBaEGgAAAAAAQC0INQAAAAAAgFoQagAAAAAAALUg1AAAAAAAAGpBqAEAAAAAANRCpVCjKIoYDocxHo9jOBxGWZYry47H4yjLcmmZPM8jz/NZnTdfAwAAAAAArFIp1Njf349erxfdbje63W4cHBysLfvTTz/FTz/9FDs7O7GzsxPD4TAiIk5OTuL169exs7MT7969izRNH/evAAAAAAAAnr3vNy1YFMXccZqmkWXZ0rJlWcZoNIputzs7NxwOo9frRUTE69ev4/LyMiIikiSp2mYAAAAAAOAF2vhJjSzLotFozJ1rNBorl466HWiMx+O544ivYYZAAwAAAAAA2NTGT2qs2j9jMpksnLsdVpRlGZPJZG6JqbIsYzweR0TE6enpyiWorq+v4/r6enZ8dXW1aXMBAAAAAIBnZuNQY5V1m4VHRPT7/Tg+Pp47d3h4OAs+0jSNt2/fxvn5+cL3DgaDeP/+/WObCAAAAAAAPAMbLz+VJMnCUxmTyWTtElJlWUaWZQtlbu/PkaZpFEWxsGdHRMTR0VF8/vx59rq4uNi0uQAAAAAAwDOzcajR6XSWnm+1Wiu/59OnTwuBRp7n8ebNm4Wyd/friIjY3d2Nvb29uRcAAAAAAPAybRxq3N3zoiiKaLVas9Aiz/OFpy3yPF8IK9I0nVuOKsuy6Ha7Ng0HAAAAAADWqrSnxmg0in6/H+12O05PT2M0Gs2uDQaDaLfb0ev15r7nbhiSJEm0Wq0YDoeRJEmcn5/P1QMAAAAAALBMpVDj9lMW3W537tqyYOJuwHGj2WxGs9ms8qMBAAAAAIAXbuPlpwAAAAAAAJ6SUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACoBaEGAAAAAABQC0INAAAAAACgFoQaAAAAAABALQg1AAAAAACAWhBqAAAAAAAAtSDUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBa+r1K4KIoYj8eRpmkURRGHh4eRJMnSsnmeR0REs9mMoiiiLMtoNpuV6wEAAAAAAIioGGrs7+/H2dlZRHwNJg4ODmI0Gi0te3JyEr/++mtERHQ6nblyVeoBAAAAAACIqBBqFEUxd5ymaWRZtrL869ev4/LyMiJi7imMqvUAAAAAAABEVNhTI8uyaDQac+cajcZsmallkiRZWFbqIfUAAAAAAABs/KRGWZZLz08mk5Xlx+NxREScnp7Gu3fvIk3TSvVcX1/H9fX17Pjq6mrT5gIAAAAAAM9MpT01llkVUtze/DtN03j79m2cn59XqmcwGMT79+8f20QAAAAAAOAZ2Hj5qSRJFp6mmEwmC8tL3bi9d0aaplEURRRFUameo6Oj+Pz58+x1cXGxaXMBAAAAAIBnZuNQo9PpLD3farUWzuV5Hm/evFk432g0KtWzu7sbe3t7cy8AAAAAAOBl2nj5qTRN546LoohWqzV7wiLP80iSJNI0jTRN4/j4eFY2y7LodrtLNw6/Ww8AAAAAAMAylfbUGI1G0e/3o91ux+npaYxGo9m1wWAQ7XY7er1eJEkSrVYrhsNhJEkS5+fnc2XX1QMAAAAAALBMpVDj9hMY3W537trdYKLZbEaz2axcDwAAAAAAwDIb76kBAAAAAADwlIQaAAAAAABALQg1AAAAAACAWhBqAAAAAAAAtSDUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBaEGgAAAAAAQC0INQAAAAAAgFoQagAAAAAAALUg1AAAAAAAAGpBqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABq4fsqhYuiiPF4HGmaRlEUcXh4GEmSLC2b53lkWRYREaenp/Hhw4dZ2TzPIyKi2WxGURRRlmU0m82H/ysAAAAAAIBnr1Kosb+/H2dnZxHxNeA4ODiI0Wi0tGyWZdHr9SIiYjgcxps3b2bfe3JyEr/++mtERHQ6nZV1AAAAAAAA3Nh4+amiKOaO0zSdPYlxV57nMRgMZsfdbjfyPJ/V8fr167i8vIzLy8v4+PHjyqc9AAAAAAAAbmwcamRZFo1GY+5co9GYLSV1W7PZjA8fPsyOy7Kclb+RJIkwAwAAAAAA2NjGy0/dBBN3TSaTpee73e7s699++y06nc4sxCjLMsbjcUR83W/j3bt3kabpQh3X19dxfX09O766utq0uQAAAAAAwDNTaU+NZVaFHbevj8fj2X4aETG3wXiapvH27ds4Pz9f+N7BYBDv379/bBMBAAAAAIBnYOPlp5IkWXgqYzKZ3LuEVL/fX9g34/b+HGmaRlEUC3t2REQcHR3F58+fZ6+Li4tNmwsAAAAAADwzG4canU5n6flWq7Xye4bDYfT7/UjTNMqyjLIsI8/zePPmzULZu/t1RETs7u7G3t7e3AsAAAAAAHiZNg417u55URRFtFqt2RMYeZ7PPW0xHo+j2WzOAo3ff/89kiSJNE3j+Ph4Vi7Lsuh2uzYNBwAAAAAA1qq0p8ZoNIp+vx/tdjtOT09jNBrNrg0Gg2i329Hr9aIoitjf35/73iRJZntptFqtGA6HkSRJnJ+fz9UDAAAAAACwTKVQ4/ZTFt1ud+7a7WAiTdOYTqcr62k2m9FsNqv8aAAAAAAA4IXbePkpAAAAAACApyTUAAAAAAAAakGoAQAAAAAA1IJQAwAAAAAAqAWhBgAAAAAAUAtCDQAAAAAAoBaEGgAAAAAAQC0INQAAAAAAgFoQagAAAAAAALUg1AAAAAAAAGpBqAEAAAAAANSCUAMAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACoBaEGAAAAAABQC99XKVwURYzH40jTNIqiiMPDw0iSpHLZKvUAAAAAAABEVAw19vf34+zsLCK+BhMHBwcxGo0ql61SDwAAAAAAQESF5aeKopg7TtM0siyrXLZKPQAAAAAAADc2DjWyLItGozF3rtFoRJ7nlcpWqQcAAAAAAODGxstPlWW59PxkMqlUtko919fXcX19PTv+/PlzRERcXV2tb+wL9H//99QteDrb/u+gL7frpfanvtwu/bk9+nJ79OV26c/t0Zfboy+3S39uj77cHn25Xfpze/TldvnsaHv839weH3Evuvncfzqdri1XaU+NZVaFFFXLLrs2GAzi/fv3C+dfvXq18c/k+fvxx6duwfOhL7dHX26X/twefbk9+nK79Of26Mvt0ZfbpT+3R19uj77cLv25Pfpyu/Tn9ujL7dGXq3358iV+XNNBG4caSZIsPE0xmUwiSZJKZavUc3R0FP/9739nx3/99VdMJpP417/+FTs7O5s2nW/o6uoqXr16FRcXF7G3t/fUzYHaMYbg4YwfeDjjBx7HGIKHM37gcYwhnrPpdBpfvnyJn3/+eW25jUONTqcTJycnC+dbrValsmmablzP7u5u7O7uzp1bFn7w9Pb29vwihUcwhuDhjB94OOMHHscYgoczfuBxjCGeq3VPaNzYONRI03TuuCiKaLVas5Ahz/NIkiTSNF1b9m4ocbceAAAAAACAZSrtqTEajaLf70e73Y7T09MYjUaza4PBINrtdvR6vXvLrrsGAAAAAACwTKVQI03TOD4+joiIbrc7d+1uMLGu7Lpr1Mvu7m7873//W1gmDNiMMQQPZ/zAwxk/8DjGEDyc8QOPYwxBxM50Op0+dSMAAAAAAADu891TNwAAAAAAAGATQg0AAAAAAKAWKu2pwctVFEWMx+NI0zSKoojDw8NIkuTRZeElyPM8siyLiIjT09P48OHDyjGR53lERDSbzSiKIsqyjGaz+Xc1Ff5xqowJ9x9YNB6Po9PpRETcOx7cg+DrODg4OIizs7O58+ZDcL9V48d8CO63bvxEmA/BgilsoNlszr4+Pz+fdrvdrZSFl+D4+Hju69tj5K7Dw8NpREwjYtrpdKaXl5d/Qwvhn6vKmHD/gUU34+f26/Z96Tb3IF660Wg0PTs7my6bJpsPwXrrxo/5EKy3bvyYD8FyntTgXkVRzB2naTr7K4vHlIWXIM/zGAwG0ev1IiKi2+1Gv9+PoigiTdOF8q9fv47Ly8uIuP8vauEl2HRMuP/AorIsYzQaRbfbnZ0bDoeze9Jd7kG8dLfHym3mQ3C/VePHfAjut2r8RJgPwSr21OBeWZZFo9GYO9doNGaPwD20LLwEzWYzPnz4MDsuyzIiYmGc3JYkiTfwcMsmY8L9B5a7PUkej8drJ80R7kGwjPkQPJz5EDye+RAs8qQG97p503HXZDJ5VFl4KW5/gPTbb79Fp9NZ+YakLMsYj8cR8XW92Xfv3i39CyZ4KTYdE+4/sOj2vaYsy5hMJmvvKe5BsJz5EDyO+RA8nPkQLCfU4MFW/cJ8bFl4rm7ejNzd+Ou22xt5pWkab9++jfPz87+phfDP89gx4f4DX/X7/Tg+Pl5bxj0IqjEfgmrMh6A68yFYzvJT3CtJkoVkdzKZLP3Liipl4aXp9/vx8ePHjdfBTNM0iqJYWBsTXpJNx4T7D6xWlmVkWXbveHAPguXMh2A7zIegOvMhWE6owb06nc7S861W61Fl4SUZDofR7/cjTdMoy3LpX0vkeR5v3rxZOL9uvVl4zqqMCfcfWO3Tp0/3Tmjdg2A18yF4PPMhqM58CFYTanCvu2v1FUURrVZrNjnO83yWEt9XFl6i8XgczWZz9gb+999/Xzl+bi8NkmVZdLtd44cX674x4f4Dm8nzfOnk1z0IVrv9gav5EFRzN7AwH4LN3b3/mA/BcjvT6XT61I3gn68oijg5OYl2ux2np6dxdHQ0+8W4v78f7XY7er3evWXhpSmKIv7zn//MnUuSJC4vLyNicfzkeT5bIuT8/Pze9c/huVs3Jtx/YDPD4TDOz8/j5ORk7rx7EMzLsiw+fvwYw+Ewer1etNvt2QbH5kOw3qrxYz4E91t3/zEfguWEGgAAAAAAQC1YfgoAAAAAAKgFoQYAAAAAAFALQg0AAAAAAKAWhBoAAAAAAEAtCDUAAAAAAIBaEGoAAAAAAAC1INQAAAAAAABqQagBAAAAAADUglADAAAAAACoBaEGAAAAAABQC0INAAAAAACgFoQaAAAAAABALfw/w15lfvE0ouYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1600x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute spread using error propagation\n",
    "\n",
    "sigma_dqxs = 100*(dbqx_dk @ (0.01*torch.eye(nq,  dtype=torch.float64))**2 @ dbqx_dk.T).diag().sqrt()/bqx\n",
    "sigma_dqys = 100*(dbqy_dk @ (0.01*torch.eye(nq,  dtype=torch.float64))**2 @ dbqy_dk.T).diag().sqrt()/bqy\n",
    "\n",
    "plt.figure(figsize=(16, 2))\n",
    "plt.bar(range(len(sigma_dqxs)), sigma_dqxs.cpu().numpy(), color='blue', alpha=0.75, width=0.75)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(16, 2))\n",
    "plt.bar(range(len(sigma_dqys)), sigma_dqys.cpu().numpy(), color='blue', alpha=0.75, width=0.75)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "myt0_gMIOq7b",
    "5d97819c"
   ],
   "name": "03_frequency.ipynb",
   "provenance": []
  },
  "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.12.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}