{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c9b316b3-8d04-4b55-8395-54be3213bd03",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pykep as pk\n",
    "import numpy as np\n",
    "import json\n",
    "import pickle as pkl\n",
    "import requests\n",
    "import pygmo as pg\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "509d978f-c0b2-422d-9c0b-fbe4cfdc7642",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cascade as csc\n",
    "from copy import deepcopy\n",
    "from tqdm.notebook import tqdm\n",
    "import heyoka as hy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a5f6485",
   "metadata": {},
   "source": [
    "# We obtain atmospheric data from the nrlmsise00 model\n",
    "An account is needed at https://amentum.com.au/atmosphere and an API-Key allowing a limited amounts of queries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fd6d6819",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will contain the density in kg/m3\n",
    "density = []\n",
    "# This is the grid on the altitude\n",
    "altitude = np.linspace(78, 1000, 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27036d7d",
   "metadata": {},
   "source": [
    "... uncomment and run this cell to retreive the data again, else use the pickled file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "32aabb15",
   "metadata": {},
   "outputs": [],
   "source": [
    "#url = \"https://atmosphere.amentum.io/nrlmsise00\"\n",
    "\n",
    "#headers = {\"API-Key\" : \"Pmd5s5e36Ny0zpPGhxvDWRtgrmB2rR6K\"}\n",
    "#for alt in altitude:\n",
    "#    params = {\n",
    "#        'year' : 2021,\n",
    "#        'month' : 2,\n",
    "#        'day' : 1,\n",
    "#        'geodetic_latitude' : 0,\n",
    "#        'geodetic_longitude' : 42,\n",
    "#        'altitude' : alt, # km\n",
    "#        'utc' : 2, # hours\n",
    "#    }\n",
    "#    # handle exceptions\n",
    "#    response = requests.get(url, params=params, headers=headers)\n",
    "#    json_payload = response.json()\n",
    "#    density.append(json_payload['total_mass_density']['value'])\n",
    "\n",
    "#with open('data/atmosphere.pk', 'wb') as file:\n",
    "#    pkl.dump((altitude, density), file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e323550c",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/atmosphere.pk\", \"rb\") as file:\n",
    "    altitude, density = pkl.load(file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f34c0a",
   "metadata": {},
   "source": [
    "# We try and fit the data with some predefined shape\n",
    "It is important that the limit for h that goes to infinity of the density goes to zero. The chosen shape must guarantee this property. We here use:\n",
    "\n",
    "$$\n",
    "\\rho(h) = \\sum_i\\alpha_i\\exp(-\\beta_i(h-\\gamma_i))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b6b0fb6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "class my_udp:\n",
    "    def __init__(self, X, Y, n=4):\n",
    "        self.X = X\n",
    "        self.Y = Y\n",
    "        self.n = n\n",
    "    def fitness(self,x):\n",
    "        p1 = x[:self.n]\n",
    "        p2 = x[self.n:2*self.n]\n",
    "        p3 = x[2*self.n:]\n",
    "        Y_pred = self._fitting_function(self.X,p1,p2, p3)\n",
    "        return [np.mean(np.abs(np.log10(self.Y)-np.log10(Y_pred)))]\n",
    "    def _fitting_function(self, x, p1,p2, p3):\n",
    "        retval = 0.\n",
    "        for alpha,beta, gamma in zip(p1,p2, p3):\n",
    "            retval += alpha*np.exp(-(x-gamma)*beta)\n",
    "        return retval\n",
    "    def get_bounds(self):\n",
    "        lb = [0]*self.n*3\n",
    "        ub = [1.]*self.n*3\n",
    "        ub[self.n*2:] = [1000]*self.n\n",
    "        return (lb, ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3a04c85d",
   "metadata": {},
   "outputs": [],
   "source": [
    "udp = my_udp(altitude, density, n=4)\n",
    "prob = pg.problem(udp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a66d2ddc",
   "metadata": {},
   "outputs": [],
   "source": [
    "uda = pg.sade(50, memory=True)\n",
    "algo = pg.algorithm(uda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "62a34622",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3576922/3502744231.py:15: RuntimeWarning: overflow encountered in exp\n",
      "  retval += alpha*np.exp(-(x-gamma)*beta)\n"
     ]
    }
   ],
   "source": [
    "pop  =pg.population(prob, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "bdfde242",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3576922/3502744231.py:15: RuntimeWarning: overflow encountered in exp\n",
      "  retval += alpha*np.exp(-(x-gamma)*beta)\n",
      "/tmp/ipykernel_3576922/3502744231.py:11: RuntimeWarning: divide by zero encountered in log10\n",
      "  return [np.mean(np.abs(np.log10(self.Y)-np.log10(Y_pred)))]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.000578643324542\n",
      "3.026808000514021\n",
      "2.4888110377559025\n",
      "2.3741149857341592\n",
      "2.312471379267102\n",
      "2.2576780689300335\n",
      "2.0151634575017283\n",
      "2.0084685558333772\n",
      "1.5234155262675437\n",
      "1.3571149607583286\n",
      "1.0855638472547917\n",
      "0.8467652617569728\n",
      "0.7011640993233682\n",
      "0.5934220922537995\n",
      "0.5363040572810933\n",
      "0.5078961266510804\n",
      "0.4920599286034672\n",
      "0.48799479536958273\n",
      "0.479819505196312\n",
      "0.4716551608588495\n",
      "0.46825567427994175\n",
      "0.4484020342723868\n",
      "0.445002529220869\n",
      "0.4377128856788668\n",
      "0.427695253821798\n",
      "0.4229629427606548\n",
      "0.4200896953749397\n",
      "0.4167008624915533\n",
      "0.41361431353847544\n",
      "0.4091325606249492\n",
      "0.40309691405076675\n",
      "0.3980713123493298\n",
      "0.3912029354224048\n",
      "0.3851995593468093\n",
      "0.3788508992534813\n",
      "0.37098361940892416\n",
      "0.36163285492529007\n",
      "0.3561494602323464\n",
      "0.35368312363426146\n",
      "0.33929542194608464\n",
      "0.33364284014919504\n",
      "0.3273158916656956\n",
      "0.3135478041939112\n",
      "0.3027096069851567\n",
      "0.2835032997383042\n",
      "0.27669080196643236\n",
      "0.2713205888237881\n",
      "0.2628590563356096\n",
      "0.2598657268055497\n",
      "0.2550407640760979\n",
      "0.2540268590371742\n",
      "0.24907796851812145\n",
      "0.23959344882534878\n",
      "0.23542166501748768\n",
      "0.23288154786562096\n",
      "0.22394633525698684\n",
      "0.22007610547669385\n",
      "0.21060876115310465\n",
      "0.20443083545469265\n",
      "0.19961557424758344\n",
      "0.1940756870715757\n",
      "0.18651570742293846\n",
      "0.1781016678830613\n",
      "0.17369337142826374\n",
      "0.16698619406646173\n",
      "0.15910149674434648\n",
      "0.15403364035119416\n",
      "0.14231691494683466\n",
      "0.1395002407554455\n",
      "0.13170647530527865\n",
      "0.12324990458467226\n",
      "0.1169945714917992\n",
      "0.10960988588670306\n",
      "0.10644681964187655\n",
      "0.10273929343267044\n",
      "0.10273929343267044\n",
      "0.10273921077719367\n",
      "0.09569923646596111\n",
      "0.09156711673733456\n",
      "0.09054846927156188\n",
      "0.09039879820452842\n",
      "0.09039879820452842\n",
      "0.0878004038254625\n",
      "0.08752095929289397\n",
      "0.08241040954391866\n",
      "0.08218034438515787\n",
      "0.08088739464892686\n",
      "0.07663325219674079\n",
      "0.07470603455588634\n",
      "0.07279300242426047\n",
      "0.07220101532924912\n",
      "0.0713367984490804\n",
      "0.07095171129746897\n",
      "0.0703295143798569\n",
      "0.06906259025990472\n",
      "0.06597661269811318\n",
      "0.06483987753487981\n",
      "0.06470492758421352\n",
      "0.06389135099128224\n",
      "0.06327347705869955\n",
      "0.06301510443269336\n",
      "0.0622108180640035\n",
      "0.06156054380647122\n",
      "0.059915027795085854\n",
      "0.05968036437365327\n",
      "0.05924887238946277\n",
      "0.05907321990401259\n",
      "0.05881136787309061\n",
      "0.05875917339826256\n",
      "0.058732871143027755\n",
      "0.05872052747452336\n",
      "0.058709197730248734\n",
      "0.058709097905236235\n",
      "0.058702553036291086\n",
      "0.05870104977127192\n",
      "0.058691466841453724\n",
      "0.05869059647827932\n",
      "0.05868644283707711\n",
      "0.05868415002462084\n",
      "0.05868300613620628\n",
      "0.05867944608823571\n",
      "0.05867891662290886\n",
      "0.058677703517599766\n",
      "0.058671307114957244\n",
      "0.05866864369385569\n",
      "0.058668374882563164\n",
      "0.058668374778486125\n",
      "0.05866701177078011\n",
      "0.0586665561609856\n",
      "0.058663610768035925\n",
      "0.05866212871312979\n",
      "0.05866212867385527\n",
      "0.05866207750613922\n",
      "0.05866142486875221\n",
      "0.05866142486875221\n",
      "0.05866053690354258\n",
      "0.0586589386661164\n",
      "0.05865696730789123\n",
      "0.05865696730649865\n",
      "0.05865535101479122\n",
      "0.05865393321911934\n",
      "0.05865393321911934\n",
      "0.05865350888918775\n",
      "0.058652849077985285\n",
      "0.058652849077985285\n",
      "0.058652800136009375\n",
      "0.058651379721646905\n",
      "0.05865094543613008\n",
      "0.05864985904593143\n",
      "0.058649859038742536\n",
      "0.058649583459958084\n",
      "0.05864875771847456\n",
      "0.05864875771843256\n",
      "0.05864838156887919\n",
      "0.05864666878944143\n",
      "0.05864666878944143\n",
      "0.05864666878944143\n",
      "0.058644611306578234\n",
      "0.05864461130440735\n",
      "0.0586444820761525\n",
      "0.058643661635015464\n",
      "0.058643453688196706\n",
      "0.05864338348547432\n",
      "0.05864280823726556\n",
      "0.05864218915377721\n",
      "0.05864128413021084\n",
      "0.05862422490000797\n",
      "0.0574803998845354\n",
      "0.05724386433482603\n",
      "0.05680018295286984\n",
      "0.05665319581473068\n",
      "0.055941904113954166\n",
      "0.05589392498939354\n",
      "0.055179323872021\n",
      "0.05497894924505599\n",
      "0.05436619807613646\n",
      "0.05264784067341612\n",
      "0.0524748080396169\n",
      "0.052360250644621055\n",
      "0.05159743519239462\n",
      "0.05058877763562947\n",
      "0.05058800831702628\n",
      "0.05035916889183357\n",
      "0.050231491941899124\n",
      "0.049995200893906494\n",
      "0.049517948744033286\n",
      "0.047351070521415656\n",
      "0.04713214475522915\n",
      "0.04687439174775425\n",
      "0.046404757523906444\n",
      "0.04584436900968516\n",
      "0.04576869111071847\n",
      "0.04447393771780149\n",
      "0.04254979269326746\n",
      "0.04175180061225008\n",
      "0.04066118368592922\n",
      "0.039569957391937985\n",
      "0.03771002904231824\n",
      "0.035035479829954815\n",
      "0.03326565059259057\n",
      "0.03019763774657662\n",
      "0.028231656943684166\n",
      "0.027971907804859367\n",
      "0.027635497613244207\n",
      "0.02609629843384177\n",
      "0.026064155692279743\n",
      "0.025395030573214513\n",
      "0.025232232498192345\n",
      "0.025102533258203302\n",
      "0.024382188721835226\n",
      "0.022860985466878727\n",
      "0.02097931975492627\n",
      "0.020299067881177458\n",
      "0.019011719016095333\n",
      "0.018746802137251297\n",
      "0.017396830053945152\n",
      "0.017009976425025684\n",
      "0.017002557939223263\n",
      "0.016787579636978956\n",
      "0.016706351432038885\n",
      "0.016493265311654737\n",
      "0.016387404819675204\n",
      "0.01573781403197213\n",
      "0.015336488002839443\n",
      "0.015063306051548358\n",
      "0.014918958536443885\n",
      "0.014849704688076955\n",
      "0.014699767901876513\n",
      "0.014640773670458946\n",
      "0.014610665032930105\n",
      "0.014578601029922452\n",
      "0.014521910016482549\n",
      "0.0143297258694309\n",
      "0.014205242717723205\n",
      "0.014138558934870992\n",
      "0.014055728536150678\n",
      "0.013835197883715615\n",
      "0.0137139283754364\n",
      "0.01364279396675939\n",
      "0.013572606381602697\n",
      "0.013462189532158995\n",
      "0.013371623946791047\n",
      "0.013370868318555934\n",
      "0.013331717072994432\n",
      "0.013290775730079973\n",
      "0.01326718850160387\n",
      "0.013179647050852505\n",
      "0.013154076132504508\n",
      "0.013102914616253401\n",
      "0.013074860533430899\n",
      "0.012954243194404764\n",
      "0.012873372920293278\n",
      "0.012850418758729347\n",
      "0.01277105517889769\n",
      "0.012707442643899327\n",
      "0.012655120207127321\n",
      "0.012616762227043523\n",
      "0.012599101378908522\n",
      "0.01259299729166396\n",
      "0.012588890105825987\n",
      "0.01258577805361547\n",
      "0.01257871155758254\n",
      "0.012575610148529757\n",
      "0.012572432689884776\n",
      "0.012566735821092382\n",
      "0.012549295319679245\n",
      "0.012524355958152355\n",
      "0.01249004048503096\n",
      "0.012449934035581754\n",
      "0.012443987538412795\n",
      "0.012437354180827578\n",
      "0.012414698424235598\n",
      "0.012411719018807145\n",
      "0.012401350643362213\n",
      "0.012383839195623256\n",
      "0.012381112754200157\n",
      "0.012359693838732717\n",
      "0.012348994714147974\n",
      "0.01233517340236208\n",
      "0.012326393056010874\n",
      "0.012325767048746163\n",
      "0.012323807642349686\n",
      "0.012320942722633808\n",
      "0.012318723449328734\n",
      "0.01231781931972275\n",
      "0.012311857400181925\n",
      "0.012305015664255797\n",
      "0.012299505573035684\n",
      "0.012295613793504022\n",
      "0.012295084424239891\n",
      "0.012294211214588123\n",
      "0.012294211214588123\n",
      "0.01229368229552369\n",
      "0.012293644254650031\n",
      "0.012293294115909585\n",
      "0.012293186961553513\n",
      "0.012293186961553513\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293186606085608\n",
      "0.012293177702191231\n",
      "0.012293161044504872\n",
      "0.012293161044504872\n",
      "0.012293161044504872\n",
      "0.012293157424806944\n",
      "0.012293155368922823\n",
      "0.012293155368922823\n",
      "0.012293155368922823\n",
      "0.012293147982021564\n",
      "0.012293147982021564\n",
      "0.012293147982021564\n",
      "0.012293147982021564\n",
      "0.012293147982021564\n",
      "0.012293147982021564\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293088416335064\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293078896235788\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293071314443065\n",
      "0.012293067948687\n",
      "0.012293067948687\n",
      "0.012293067948687\n",
      "0.012293067230421872\n",
      "0.012293067230421872\n",
      "0.012293067230421872\n",
      "0.012293067230421872\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065250805242\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293065064350097\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054890810486\n",
      "0.012293054815777306\n",
      "0.012293054815777306\n",
      "0.012293054815777306\n",
      "0.012293053805091185\n",
      "0.012293053805091185\n",
      "0.012293053805091185\n",
      "0.01229305270020193\n",
      "0.01229305270020193\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293046694111407\n",
      "0.012293045729810697\n",
      "0.012293045729810697\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044977439554\n",
      "0.012293044350276689\n",
      "0.012293039029654817\n",
      "0.012293039029654817\n",
      "0.012293039029654817\n",
      "0.012293039029654817\n",
      "0.012293038375444487\n",
      "0.012293038375444487\n",
      "0.012293038375444487\n",
      "0.012293038375444487\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293035745615377\n",
      "0.012293032524686378\n",
      "0.012293032524686378\n",
      "0.012293032524686378\n",
      "0.012293032524686378\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029298193146\n",
      "0.012293029001891291\n",
      "0.012293029001891291\n",
      "0.012293029001891291\n",
      "0.012293025642992178\n",
      "0.012293025642992178\n",
      "0.012293025642992178\n",
      "0.012293025642992178\n",
      "0.012293025642992178\n",
      "0.012293025642992178\n",
      "0.012293025642992178\n",
      "0.012293025305335945\n",
      "0.012293025305335945\n",
      "0.012293025305335945\n",
      "0.012293025305335945\n",
      "0.012293025305335945\n",
      "0.012293025305335945\n",
      "0.012293025305335945\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n",
      "0.012293022270074427\n"
     ]
    }
   ],
   "source": [
    "for i in range(500):\n",
    "    pop = algo.evolve(pop)\n",
    "    print(pop.champion_f[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63fb81cc",
   "metadata": {},
   "source": [
    "# We plot the data and the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e422876e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_density(h_in_km, best_x):\n",
    "    \"\"\"\n",
    "    returns the atmosheric density in kg.m^3\n",
    "    \"\"\"\n",
    "    p1 = np.array(best_x[:4])\n",
    "    p2 = np.array(best_x[4:8])\n",
    "    p3 = np.array(best_x[8:])\n",
    "    retval = 0.\n",
    "    for alpha,beta, gamma in zip(p1,p2, p3):\n",
    "        retval += alpha*np.exp(-(h_in_km-gamma)*beta)\n",
    "    return retval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "846521fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f5718864040>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/jklEQVR4nO3deXzU1b3/8fc3CUkIJANJIJCNzQ2IYglerwgKpYVGTAu1rVsBt1Z+XpeYq62WbtJrY2ur0Hrhihta28ptFao+vEJaQbAoSCAVBBcwEJZgWDNJkCRkzu+PJEMmmYSEzPqd1/PxyEPynW++c05A8uaczznHMsYYAQAAhLmoYDcAAADAFwg1AADAFgg1AADAFgg1AADAFgg1AADAFgg1AADAFgg1AADAFgg1AADAFmKC3YBAcblcOnDggBITE2VZVrCbAwAAusAYo+rqaqWnpysqqvOxmIgJNQcOHFBWVlawmwEAAM7C3r17lZmZ2ek9ERNqEhMTJTV9U5KSkoLcGgAA0BVOp1NZWVnun+OdiZhQ0zLllJSURKgBACDMdKV0hEJhAABgC4QaAABgC4QaAABgCxFTUwMACD5jjE6dOqXGxsZgNwUhpFevXoqOju7xcwg1AICAqK+vV0VFhU6cOBHspiDEWJalzMxM9e3bt0fPIdQAAPzO5XKprKxM0dHRSk9PV2xsLBuhQlLT6N2hQ4e0b98+nXvuuT0asSHUAAD8rr6+Xi6XS1lZWUpISAh2cxBiBgwYoN27d6uhoaFHoYZCYQBAwJxpm3tEJl+N2vGnCwAA2AKhBgCAIJs0aZIKCgp69Izdu3fLsiyVlpb6pE1S0wjKihUrOnzdGKPvf//7Sk5Odr+3L/pytqipAQAgzNx00006fvy4R+DIyspSRUWFUlNTA9aON998U0uXLtWaNWs0fPhwpaam6pVXXlGvXr3c9wwdOlQFBQUBCTqEGl+o2i8d3SUlj5AcGcFuDQAgAkVHR2vQoEEBfc9du3Zp8ODBGj9+vPtacnJyQNvQGtNPPfT5miVyPZ4jPZ8vLciRNr8Q7CYBAHzkzTff1IQJE9SvXz+lpKTo6quv1q5du9yvt0z5vPLKK5o8ebISEhI0ZswYvfvuu+57jhw5ouuvv16ZmZlKSEjQhRdeqD//+c8dvuf8+fN14YUXtruem5urn/70p/r5z3+u559/Xn/7299kWZYsy9KaNWu8Tj99+OGHmj59upKSkpSYmKiJEye62//+++/rq1/9qlJTU+VwOHTllVdq8+bNXf7e3HTTTbrrrrtUXl4uy7I0dOhQSZ5TaZMmTdKePXt07733utvqT4Sanqjar4FrfqgouZo+Ny7ptYKmkRsAQIeMMTpRfyooH8aYLreztrZWhYWFev/99/WPf/xDUVFRmjlzplwul8d98+bN03333afS0lKdd955uv7663Xq1ClJ0smTJ5Wbm6vXX39d27Zt0/e//33NmjVLGzZs8Pqet9xyi7Zv367333/ffe2DDz7Qli1bdNNNN+m+++7Td77zHX3ta19TRUWFKioqPEZKWuzfv19XXHGF4uPj9dZbb6mkpES33HKLu13V1dWaM2eO1q1bp/fee0/nnnuurrrqKlVXV3fpe7Nw4ULNnz9fmZmZqqio8Ghvi1deeUWZmZmaP3++u63+xPRTTxzdJUuef7BlGqWjnzENBQCd+KKhUaN+ujIo7719/jQlxHbtx98111zj8fkzzzyjgQMHavv27crJyXFfv++++zR9+nRJ0kMPPaTRo0dr586duuCCC5SRkaH77rvPfe9dd92lN998U3/5y1906aWXtnvPzMxMTZs2Tc8995wuueQSSdJzzz2nK6+8UsOHD5ck9e7dW3V1dZ1ON/33f/+3HA6HXnrpJXeNy3nnned+/ctf/rLH/U8++aT69++vt99+W1dfffUZvzcOh0OJiYmdTnslJycrOjpaiYmJAZkaY6SmJ5JHyFhtvoVWtJQ8PDjtAQD41K5du3TDDTdo+PDhSkpK0rBhwyRJ5eXlHvdddNFF7l8PHjxYklRZWSlJamxs1MMPP6yLLrpIKSkp6tu3r1atWtXuGa1973vf05///GedPHlSDQ0N+uMf/6hbbrmlW20vLS3VxIkTPYp2W6usrNTcuXN13nnnyeFwyOFwqKamptN2hTpGanrCkSErf6EaX71H0XLJWNGy8hcwSgMAZ9C7V7S2z58WtPfuqvz8fGVlZempp55Senq6XC6XcnJyVF9f73Ff6+DQUjfSMkX129/+Vo8//rgWLFigCy+8UH369FFBQUG7Z7R937i4OC1fvlxxcXGqq6trN2p0Jr179+709ZtuukmHDh3SggULNGTIEMXFxemyyy7rtF2hjlDTU2Nnq3BTij7fvV03fO1KfX3svwW7RQAQ8izL6vIUULAcOXJEO3bs0JNPPqmJEydKkt55551uP2fdunX6xje+oe9+97uSmsLOp59+qpEjR3b4NTExMZozZ46ee+45xcXF6brrrvM4XiI2NvaMJ51fdNFFev7559XQ0OB1tGbdunVatGiRrrrqKknS3r17dfjw4W7370y60lZfYfrJBxIHZus91yh9dCIp2E0BAPhI//79lZKSoiVLlmjnzp166623VFhY2O3nnHPOOSouLtb69eu1Y8cO3X777Tp48OAZv+62227TW2+9pf/7v/9rN/U0dOhQffDBB/r44491+PBhNTQ0tPv6O++8U06nU9ddd502bdqkTz/9VH/4wx/08ccfu9v1hz/8QTt27NCGDRt04403nnF052wMHTpUa9eu1f79+/0Smloj1PjAkOQ+kqQ9R08EuSUAAF+JiorSSy+9pJKSEuXk5Ojee+/Vo48+2u3n/OQnP9HYsWM1bdo0TZo0SYMGDdKMGTPO+HXnnnuuxo8fr/PPP79dQfH3vvc9nX/++Ro3bpwGDBigf/7zn+2+PiUlRW+99ZZqamp05ZVXKjc3V0899ZR71ObZZ5/VsWPH9KUvfUmzZs3S3XffrYEDB3a7f2cyf/587d69WyNGjNCAAQN8/vzWLNOdtW1BFhMT4642HzdunJ5++ukuf63T6ZTD4VBVVZWSknw7orLyw4O6/Q8lujDDodfumuDTZwOAHZw8eVJlZWUaNmyY4uPjg92csGCM0QUXXKDbb7/9rEaIwklnfz668/M7tCc02+jXr59Pz7TwlSEpTfOce47UBrklAAA7qKys1B/+8Aft379fN998c7CbEzbCKtSEquzkplDjPHlKx0/Uq19CbJBbBAAIZ2lpaUpNTdWSJUvUv3//YDcnbPispmbt2rXKz89Xenp6h6d6Llq0yD20lJubq3Xr1nXrPZxOp3JzczVhwgS9/fbbPmp5zyXExmhAYpwkac8R6moAAD1jjNGhQ4d0ww03BLspYcVnIzW1tbUaM2aMbr75Zq9r6ZctW6aCggItWrRIl19+uZ588knl5eVp+/btys7OltR0rkVdXV27r121apXS09O1e/dupaena9u2bZo+fbq2bt3q8/qYszUkOUGHquu05+gJjcnqF+zmAAAQcXwWavLy8pSXl9fh64899phuvfVW3XbbbZKkBQsWaOXKlVq8eLGKiookSSUlJZ2+R3p6uiQpJydHo0aN0ieffKJx48Z5vbeurs4jIDmdzm71p7uyUxK0ac8x7WUFFAAAQRGQJd319fUqKSnR1KlTPa5PnTpV69ev79Izjh075g4p+/bt0/bt291nYHhTVFTk3vbZ4XAoKyvr7DvQBe5l3RQLAwAQFAEJNYcPH1ZjY6PS0tI8rqelpXVpAyJJ2rFjh8aNG6cxY8bo6quv1sKFC5WcnNzh/Q8++KCqqqrcH3v37u1RH84kO6VpwyJqagAACI6Arn5qOQ+jhTGm3bWOjB8/Xlu3bu3ye8XFxSkuLq5b7euJ7OaRmnKmnwAACIqAjNSkpqYqOjq63ahMZWVlu9GbcNWyV42c+1X/6Rqpan9Q2wMAQKQJSKiJjY1Vbm6uiouLPa4XFxdr/PjxgWiC36X0idWs2Lf1Tuzdiv3jN6QFOdLmF4LdLACAH0yaNEkFBQXBbgba8FmoqampUWlpqXvH37KyMpWWlqq8vFySVFhYqKefflrPPvusduzYoXvvvVfl5eWaO3eur5oQVJbzgB6KWqJoq/nUCeOSXitgxAYAItyaNWtkWZaOHz8e7KbYns9qajZt2qTJkye7P285p2LOnDlaunSprr32Wh05ckTz589XRUWFcnJy9MYbb2jIkCG+akJwHd2lKLU5Rss0Skc/kxwZwWkTAAARxGcjNZMmTZIxpt3H0qVL3ffccccd2r17t+rq6lRSUqIrrrjCV28ffMkj5Gr77bSipeSOl50DAM5C1X6pbG3ARsJra2s1e/Zs9e3bV4MHD9Zvf/tbj9dffPFFjRs3TomJiRo0aJBuuOEGVVZWSpJ2797t/gd///79ZVmWbrrpJknSm2++qQkTJqhfv35KSUnR1VdfrV27dgWkT3YVkJqaiODI0Macn+qUaf6WWtFS/gJGaQDAlza/0FSz+Hx+wGoX77//fq1evVrLly/XqlWrtGbNGo/NYuvr6/WLX/xC//rXv7RixQqVlZW5g0tWVpZefvllSdLHH3+siooKLVy4UFJTWCosLNT777+vf/zjH4qKitLMmTPlcrn83ie74kBLHzo1ZpYmbErR+P5VemzuNwk0AOBLVful1+5pqlmUTtcujpjit79va2pq9Mwzz+iFF17QV7/6VUnS888/r8zMTPc9t9xyi/vXw4cP1+9+9zv927/9m2pqatS3b1/3nmoDBw5Uv3793Pe2PVLomWee0cCBA7V9+3bl5OT4pT92x0iNDw1JSdBBpej16nPkSkwPdnMAwF6O7jodaFq01C76ya5du1RfX6/LLrvMfS05OVnnn3+++/MtW7boG9/4hoYMGaLExERNmjRJktwLZTp79g033KDhw4crKSlJw4YN69LXoWOEGh8a7IhXTJSl+lMuHXSeDHZzAMBekkdIVmBrF40xnb5eW1urqVOnqm/fvnrxxRf1/vvva/ny5ZKapqU6k5+fryNHjuipp57Shg0btGHDhi59HTpGqPGhmOgoZfTnuAQA8AtHhpS/sCnISAGpXTznnHPUq1cvvffee+5rx44d0yeffCJJ+uijj3T48GE98sgjmjhxoi644AJ3kXCL2NhYSVJjY6P72pEjR7Rjxw79+Mc/1pQpUzRy5EgdO3bMb/2IFNTU+Fh2coL2HDmh8qO1umxESrCbAwD2MnZ2Uw3N0c+aRmj8XLvYt29f3Xrrrbr//vuVkpKitLQ0zZs3T1FRTWMC2dnZio2N1e9//3vNnTtX27Zt0y9+8QuPZwwZMkSWZen111/XVVddpd69e6t///5KSUnRkiVLNHjwYJWXl+uBBx7wa18iASM1PtZyXAIjNQDgJ44MadjEgC3GePTRR3XFFVfo61//ur7yla9owoQJys3NlSQNGDBAS5cu1V/+8heNGjVKjzzyiH7zm994fH1GRoYeeughPfDAA0pLS9Odd96pqKgovfTSSyopKVFOTo7uvfdePfroowHpj51Z5kwThjbhdDrlcDhUVVWlpKQkv73PU2s/08Nv7ND0iwbrv28Y67f3AYBwcvLkSZWVlWnYsGGKj48PdnMQYjr789Gdn9+M1PhYdvNITTkjNQAABBShxsdOTz/VBrklAABEFkKNj2UnN4Ua58lTqjrREOTWAAAQOQg1PpYQG6MBiXGSpD1HGa0BACBQCDV+MCSZFVAAAAQaocYP3MXCRwk1ANBahCy4RTf56s8FocYPhqb0kSTtPsz0EwBIUq9evSRJJ07wjz2013I0RHR0dI+ew47CfsAGfADgKTo6Wv369XMfIZCQkCDLsoLcKoQCl8ulQ4cOKSEhQTExPYslhBo/GNI8UvPFkXKprKHpELYA7XwJAKFq0KBBktTubCQgKipK2dnZPQ66hBo/GJqSoO9Er1ZR/dPS86bpVNn8hU1nlgBAhLIsS4MHD9bAgQPV0MCWFzgtNjbWfZ5WTxBq/KBfwyEV9Xpa0WoufDIu6bWCpkPYGLEBEOGio6N7XDsBeEOhsD8c3XU60LQwjU2nygIAAL8g1PhD8gi52n5rrWgpeXhw2gMAQAQg1PiDI0Orhj+oU6b522tFS/kLmHoCAMCPqKnxk+pR12vC9kHKzzypebOmE2gAAPAzRmr8ZEhKHx1Uit48cQ6BBgCAACDU+MnQ5g349h/7QvWnXEFuDQAA9keo8ZMBiXHq3StaLiPtP/5FsJsDAIDtEWr8xLIs93EJu49wBhQAAP5GqPEj9xlQHGwJAIDfEWr8qOUMqD1HOdgSAAB/I9T4Ead1AwAQOIQaPxraMlJDTQ0AAH5HqPGj7OSmkZq9R79Qo8uc4W4AANAThBo/Su/XW72iLdU3ulRRxbJuAAD8iVDjR9FRlrKaR2vKqasBAMCvCDV+NiS5Za8aQg0AAP5EqPGz08u6KRYGAMCfCDV+dnoDPkZqAADwJ0KNn7Us6+aoBAAA/ItQ42ctIzXlR0/IGJZ1AwDgL4QaP8vsn6AoSzpR36hDNXXBbg4AALZFqPGz2JgopffrLYll3QAA+BOhJgCGpvTRIB1R7cerpar9wW4OAAC2FBPsBkSCmXpLM+J+peh3jfRelJS/UBo7O9jNAgDAVhip8beq/Zq579eKtpqLhI1Leq2AERsAAHwsbELNxx9/rIsvvtj90bt3b61YsSLYzTqzo7sUJZfnNdMoHf0sOO0BAMCmwmb66fzzz1dpaakkqaamRkOHDtVXv/rV4DaqK5JHyFhRskyrYGNFS8nDg9cmAABsKGxGalp79dVXNWXKFPXp0yfYTTkzR4ZOXfW4Tpmmb7WxoqX8BZIjI7jtAgDAZnwWatauXav8/Hylp6fLsiyvU0OLFi3SsGHDFB8fr9zcXK1bt+6s3ut///d/de211/awxYHT65Kb9J3eT+q6+h+r9JtrKRIGAMAPfBZqamtrNWbMGD3xxBNeX1+2bJkKCgo0b948bdmyRRMnTlReXp7Ky8vd9+Tm5ionJ6fdx4EDB9z3OJ1O/fOf/9RVV13lq6YHhCNtqN5zjdL2E4nBbgoAALbks5qavLw85eXldfj6Y489pltvvVW33XabJGnBggVauXKlFi9erKKiIklSSUnJGd/nb3/7m6ZNm6b4+PhO76urq1Nd3ekdfJ1OZ1e64TcjBvTV6o8PaVclZ0ABAOAPAampqa+vV0lJiaZOnepxferUqVq/fn23ntXVqaeioiI5HA73R1ZWVrfex9fOGdhXkrTzUE1Q2wEAgF0FJNQcPnxYjY2NSktL87ielpamgwcPdvk5VVVV2rhxo6ZNm3bGex988EFVVVW5P/bu3dvtdvvSiOZQs6uSUAMAgD8EdEm3ZVkenxtj2l3rjMPh0Oeff96le+Pi4hQXF9et9vnTOQOaQs3+41/oRP0pJcSGzWp6AADCQkBGalJTUxUdHd1uVKaysrLd6I1d9e8Tq+Q+sZKkzw5RVwMAgK8FJNTExsYqNzdXxcXFHteLi4s1fvz4QDQhJLSM1uyirgYAAJ/z2RxITU2Ndu7c6f68rKxMpaWlSk5OVnZ2tgoLCzVr1iyNGzdOl112mZYsWaLy8nLNnTvXV00IeSMG9tHG3UepqwEAwA98Fmo2bdqkyZMnuz8vLCyUJM2ZM0dLly7VtddeqyNHjmj+/PmqqKhQTk6O3njjDQ0ZMsRXTQh5IwawAgoAAH/xWaiZNGmSjDGd3nPHHXfojjvu8NVbhp3TK6CoqQEAwNfC8uyncNVSU1N2uFaNrs4DIAAA6B5CTQBl9OutuJgo1Te6tPfoiWA3BwAAWyHUBFBUlKXhrIACAMAvCDUB5j4ugRVQAAD4FKEmwEYM6COJkRoAAHyNUBNg5wzsq0E6ovh9/5Sq9ge7OQAA2AahJsDGHnld/4y7W/OPPyizIEfa/EKwmwQAgC0QagKpar8Gr/2hoq2m5dyWcUmvFTBiAwCADxBqAunorqYg05pplI5+Fpz2AABgI4SaQEoeIVltvuVWtJQ8PDjtAQDARgg1geTIkPIXytX8bXcpSspf0HQdAAD0iM/OfkIXjZ2t16ov0J9Xvq2BQ0bqd2OvDnaLAACwBUJNEGQMOUfvuQ4p41jvYDcFAADbYPopCEY0H5Ww//gXqq07FeTWAABgD4SaIOjfJ1YpfWIlNZ3YDQAAeo5QEyQtozWcAQUAgG8QaoLknLSmULOjwhnklgAAYA+EmiAZN6S/JOm9sqNBbgkAAPZAqAmSS4enSJK27a9S9cmGILcGAIDwR6gJkox+vZWdnKABrsPaufH/OP8JAIAeYp+aILq737uaWftrRb9lpNVRUv5CaezsYDcLAICwxEhNsFTt1zUHHnWf2C1O7AYAoEcINcFydJcscWI3AAC+QqgJFk7sBgDApwg1wcKJ3QAA+BSFwsE0drZWnhyl519frYS0c/Xs2BnBbhEAAGGLkZogGzNqtN5zjdKag73YrwYAgB4g1ARZer/eGpKSIJeRNu0+FuzmAAAQtgg1IeDfhzXtLvzhRzuksrUs6wYA4CxQUxMC/n1EssyWF3RH6dNSqWlaFcVGfAAAdAsjNSFg/IA6FcU8rSixER8AAGeLUBMC0hr2n95ZuAUb8QEA0C2EmlCQPMK9X40bG/EBANAthJpQ4MjQvy5+SKdM82+HFc1GfAAAdBOFwiFi8OTva8J7/TQs6nMtvOMaDcxklAYAgO5gpCZEDHLEK3vYuXrXNUrPbq0PdnMAAAg7hJoQ8r2JTaMz/9iwWV98sprVTwAAdAOhJoRMuWCg/sOxXm+aO9T7TzOkBTnS5heC3SwAAMICoSaERFUf0H/WLTq9vJv9agAA6DJCTSg5uktRcnleY78aAAC6hFATSpJHNB2R0IphvxoAALqEUBNKHBlS/sKmICPplInS7sseZr8aAAC6gFATasbOllWwVUuG/U4T6hZq8b6hnNwNAEAXEGpCkSNDX877pq6I/kBFe66Xns9nJRQAAGcQVqHmN7/5jUaPHq2cnBy9+OKLwW6OX50T59QjvZ5hJRQAAF0UNqFm69at+tOf/qSSkhJt2rRJixcv1vHjx4PdLP9hJRQAAN0SNqFmx44dGj9+vOLj4xUfH6+LL75Yb775ZrCb5T+shAIAoFt8FmrWrl2r/Px8paeny7IsrVixot09ixYt0rBhwxQfH6/c3FytW7euy8/PycnR6tWrdfz4cR0/flxvvfWW9u+38VSMl5VQO3LnsxIKAIAO+OyU7traWo0ZM0Y333yzrrnmmnavL1u2TAUFBVq0aJEuv/xyPfnkk8rLy9P27duVnZ0tScrNzVVdXV27r121apVGjRqlu+++W1/+8pflcDh0ySWXKCbG5oeMj50ta8QU/enNNfrdlkZlb0/QHy94W70GnEO4AQCgDcsYY3z+UMvS8uXLNWPGDPe1Sy+9VGPHjtXixYvd10aOHKkZM2aoqKio2+9x2223aebMmZo+fbrX1+vq6jwCktPpVFZWlqqqqpSUlNTt9wum6pMNevxXP9E81/80FQ5bUVL+Qmns7GA3DQAAv3I6nXI4HF36+R2Qmpr6+nqVlJRo6tSpHtenTp2q9evXd/k5lZWVkqSPP/5YGzdu1LRp0zq8t6ioSA6Hw/2RlZV1do0PAYl1lfqJeZKVUAAAdCIg8zeHDx9WY2Oj0tLSPK6npaXp4MGDXX7OjBkzdPz4cfXp00fPPfdcp9NPDz74oAoLC92ft4zUhKWju2R1tBKKaSgAACQFKNS0sCzL43NjTLtrnenOqE5cXJzi4uK6fH9Ia1kJZU4HG2NFy2IlFAAAbgGZfkpNTVV0dHS7UZnKysp2ozfwonkllFqthHqq390ySelBbhgAAKEjIKEmNjZWubm5Ki4u9rheXFys8ePHB6IJ4W/sbKlgqz6f+bImN/5ez1YM14a3VlBXAwBAM59NP9XU1Gjnzp3uz8vKylRaWqrk5GRlZ2ersLBQs2bN0rhx43TZZZdpyZIlKi8v19y5c33VBPtzZChtTIZ+svVRTfn0YUWvMzLvRMliJRQAAL4LNZs2bdLkyZPdn7cU6c6ZM0dLly7VtddeqyNHjmj+/PmqqKhQTk6O3njjDQ0ZMsRXTYgMVfv11V2/lNW8EspqWQk1YgpFwwCAiOaXfWpCUXfWuYe0srVNp3a3Ned1adjEwLcHAAA/Crl9auBDXs6EalSUTPKwIDUIAIDQQKgJN21WQjWaKD3ScK02bHqfomEAQERj+ilcVe2Xjn6mtWtW6fLdv1e0ZWQsioYBAPbC9FMkcGRIycM1sfwJ9/EJFscnAAAiGKEmnB3d1RRkWms5PgEAgAhDqAlnHRQNN/anaBgAEHkINeGsbdGwmoqGi9etZwoKABBxKBS2g+ai4ZL3/qGLP1pA0TAAwDYoFI40zUXDYz9ZSNEwACBiEWrsgqJhAECEI9TYhdeiYUtfHD/IaA0AICIQauyiTdGwS5YsI/X+223Sghxp8wtBbiAAAP5FqLGTsbOlgq3St5bKkqWo5voaUV8DAIgAhBq7cWRIfVJkifoaAEBkIdTYkZf6GpeipOThQWoQAAD+R6ixozb1NadMlH516jrt3bWVKSgAgG2x+Z6dVe2XObpLf3n1NV1z9Ck25QMAhB0230MTR4as5BH69vGn2ZQPAGB7hBq7Y1M+AECEINTYXQeb8jVWVzJaAwCwFUKN3XWwKV/0K7ewKR8AwFYINZGATfkAABGAUBMp2JQPAGBzhJpI4nVTPkuqPcRoDQAg7BFqIkmb+ppGY0lG0l9vpr4GABD2CDWRprm+xnzrOVmWqK8BANgGoSYSOTJk9UlVlNpsJk19DQAgjBFqIpWX+hpDfQ0AIIwRaiKVl/oaI1FfAwAIW4SaSNZcX1M389mm+hpRXwMACF+EmkjnyFBc0gDqawAAYY9Qgw7qa6KorwEAhBVCDdz1NcajvsZQXwMACCuEGjQZO1tWwVZVXf0U9TUAgLBEqMFpjgw5UgZRXwMACEuEGniivgYAEKYINfBEfQ0AIEwRatBec33NkbwlEvU1AIAwQaiBd44MpQwcrGjqawAAYYJQg45RXwMACCOEGnSM+hoAQBgh1KBzzfU1TvavAQCEOEINzsyRoST2rwEAhDhCDbrGa32NRX0NACBkhGSomTlzpvr3769vfetb3XoNfuS1vkbU1wAAQkZIhpq7775bL7zg/YdkZ6/Bz5rra+pnPkt9DQAg5IRkqJk8ebISExO7/RoCwJGh2KQB3utrPlxBsAEABE23Q83atWuVn5+v9PR0WZalFStWtLtn0aJFGjZsmOLj45Wbm6t169b5oq0IFV7raySt+hFTUQCAoOl2qKmtrdWYMWP0xBNPeH192bJlKigo0Lx587RlyxZNnDhReXl5Ki8vd9+Tm5urnJycdh8HDhw4+560UVdXJ6fT6fEBH2mur1FzfY0xktXyGlNRAIAgienuF+Tl5SkvL6/D1x977DHdeuutuu222yRJCxYs0MqVK7V48WIVFRVJkkpKSs6yuV1XVFSkhx56yO/vE7HGzpZGTJE+XCFr1Y88X2tZ6u3ICE7bAAARyac1NfX19SopKdHUqVM9rk+dOlXr16/35Vud0YMPPqiqqir3x969ewP6/hHBkSGNniHDUQoAgBDQ7ZGazhw+fFiNjY1KS0vzuJ6WlqaDBw92+TnTpk3T5s2bVVtbq8zMTC1fvlyXXHLJGV9rLS4uTnFxcT3rEM7MkSErf6HMawWyTKMajSXLMrL+enNT3U3+wqZRHQAA/MynoaaFZVkenxtj2l3rzMqVK8/qNQTJ2NmyRkzR4R3vqP+bt7df6j1iClNRAAC/8+n0U2pqqqKjo9uNylRWVrYbvYHNODKUmjZY0Sz1BgAEiU9DTWxsrHJzc1VcXOxxvbi4WOPHj/flWyEUsdQbABBE3Q41NTU1Ki0tVWlpqSSprKxMpaWl7iXbhYWFevrpp/Xss89qx44duvfee1VeXq65c+f6tOEIQSz1BgAEUbdrajZt2qTJkye7Py8sLJQkzZkzR0uXLtW1116rI0eOaP78+aqoqFBOTo7eeOMNDRkyxHetRuhqXuptPlwua9U8z9dapqJGz6DGBgDgc5Yxxpz5tvDndDrlcDhUVVWlpKSkYDfH/qr2yyzIkWVc7ktGzSM3rIoCAHRRd35+h+TZT7CBlqXeTEUBAAKEUAP/aT7V23nlQ2q3or9l12EAAHyEUAP/cmQoaey32XUYAOB3hBr4X5upqEZjSTLSX29mqTcAwGcINQiM5qmobeMXSpKstrsOM2IDAOghQg0Cx5GhnHOHK9pi12EAgO8RahBY7DoMAPATQg0Ci12HAQB+QqhB4I2dLRVslZn6sPel3ns3BqVZAIDwRqhBcDgyZI2e2W6ptyTp5VuYhgIAdBuhBsHTstS77R9DpqEAAGeBUIPgGjtb1reeaX+dFVEAgG4i1CD4si5lRRQAoMcINQi+tiuixIooAED3EWoQGppXRGnqL9V2QRRTUQCAriDUIHQ4MqTRM5iKAgCcFUINQgtTUQCAs0SoQehhKgoAcBYINQhNTEUBALqJUIPQxVQUAKAbCDUIbUxFAQC6iFCD0MdUFACgCwg1CA9tp6IMU1EAAE+EGoSP1lNRbeeiTKN09LOgNAsAEBoINQgvHUxFNcpSXdXnjNYAQAQj1CD8tJmKchlLlpHiVtwqQ30NAEQsQg3CU8tU1LeWyrIsRVlGkmQZl8yr90jbXmHUBgAiDKEG4cuRIfVJkSWXx2VLLumvN7MqCgAiDKEG4S15RLv6GjdWRQFARCHUILy1qa9phw36ACBiEGoQ/lrV17BBHwBELkIN7MGRIeXM5KwoAIhghBrYC2dFAUDEItTAfjgrCgAiEqEG9tT2rCgxFQUAdkeogX0xFQUAEYVQA3tjKgoAIgahBvbHVBQARARCDSIDU1EAYHuEGkQOpqIAwNYINYgsTEUBgG0RahB5mIoCAFsi1CAydTQVZcRUFACEqZAMNTNnzlT//v31rW99y+N6dXW1LrnkEl188cW68MIL9dRTTwWphbAFb1NRLUM3TEUBQNgJyVBz991364UX2v8rOSEhQW+//bZKS0u1YcMGFRUV6ciRI0FoIWyDqSgAsI2QDDWTJ09WYmJiu+vR0dFKSEiQJJ08eVKNjY0yxgS6ebCbM0xFGaaiACAsdDvUrF27Vvn5+UpPT5dlWVqxYkW7exYtWqRhw4YpPj5eubm5WrdunS/aKkk6fvy4xowZo8zMTP3gBz9Qamqqz56NCNbJVJRlXDKv3iNte4VRGwAIYd0ONbW1tRozZoyeeOIJr68vW7ZMBQUFmjdvnrZs2aKJEycqLy9P5eXl7ntyc3OVk5PT7uPAgQNnfP9+/frpX//6l8rKyvSnP/1Jn3/+eXe7AHjXyVSUJZf015spIAaAEBbT3S/Iy8tTXl5eh68/9thjuvXWW3XbbbdJkhYsWKCVK1dq8eLFKioqkiSVlJScZXNPS0tL00UXXaS1a9fq29/+drvX6+rqVFdX5/7c6XT2+D0RAVqmoop/3FQs3FZLAfGIKU33AgBChk9raurr61VSUqKpU6d6XJ86darWr1/f4+d//vnn7nDidDq1du1anX/++V7vLSoqksPhcH9kZWX1+P0RIdpMRbVDATEAhCSfhprDhw+rsbFRaWlpHtfT0tJ08ODBLj9n2rRp+va3v6033nhDmZmZev/99yVJ+/bt0xVXXKExY8ZowoQJuvPOO3XRRRd5fcaDDz6oqqoq98fevXvPvmOIPC1TUd9ayrEKABAmuj391BWW5VmRYIxpd60zK1eu9Ho9NzdXpaWlXXpGXFyc4uLiuvyeQDuODMkxU6qvbppyMo3ej1VgKgoAQoJPR2pSU1MVHR3dblSmsrKy3egNEDbYywYAwoJPQ01sbKxyc3NVXFzscb24uFjjx4/35VsBgcWxCgAQ8ro9/VRTU6OdO3e6Py8rK1NpaamSk5OVnZ2twsJCzZo1S+PGjdNll12mJUuWqLy8XHPnzvVpw4GAaykgbj0V1epYBfPqPbJi+0pZlzIdBQBBYJlubsm7Zs0aTZ48ud31OXPmaOnSpZKaNt/79a9/rYqKCuXk5Ojxxx/XFVdc4ZMGny2n0ymHw6GqqiolJSUFtS0Ic1X7m6acVv3I++tWVFP4GTs7oM0CADvqzs/vboeacEWogU9V7W+acvK2l43UtBy8YCsjNgDQQ935+R2SZz8BIY+9bAAg5BBqgLPFXjYAEFL8sk8NEDG6spfNq/dIFBADgN8xUgP4Qmd72XAYJgAEBKEG8JUO9rJxaxm12fYKtTYA4AeEGsCXzlRAzKgNAPgNoQbwtU4KiN1azo1ixAYAfIZQA/iDI0PKmcmybwAIIEIN4E+tRm0My74BwK8INYC/NY/aWK1GbYzxsux7X0nQmggAdkCoAQKl9bLvduu+XTLPTGHEBgB6gFADBFIny74tY2RY8g0AZ41QAwRay7JvL//7WSz5BoCzRqgBgmHsbOm2v8vLPFQT42LUBgC6iVADBEtmrpT/uw6XfDNqAwDdQ6gBgqmLG/UxagMAZ8Yp3UCweTnpuy33qI0VJX3l51L6l6TkEZz6DQCtWMYYE+xGBILT6ZTD4VBVVZWSkpKC3RzAu6r90t6N0su3NO1f44VR8x43VlRTwfHY2YFsIQAEVHd+fjNSA4SSLo3aNGvZtC+2r5R1KaM2ACIeNTVAKOpKrY0kTv0GgNMYqQFClZdRG/fUU1uM2gAANTVAWKjaLx39TDqwRfr7z71OS7lRTAzARrrz85tQA4Sb7hYTE3AAhDEKhQE7624xcfFPmy+yWgqAvVEoDISrLhcTN2upu2ETPwA2xUgNEM66MGrjiU38ANgXNTWAXXgpJu5wtZSouwEQHigU9oJQg4jSndVSIuAACF2EGi8INYhYXVgt1RoBB0AoYfUTgNO6WXfTeuWUKf4pAQdA2GCkBogk3ay7aY0RHADBwPSTF4QaoA0CDoAwQKjxglADdMIfAUeSju4i7ADoEUKNF4QaoIt8EHCMrOb7DaM5AHqEUOMFoQY4C91cGt4RpqsAnC1CjReEGqCH/BlwevWRGmoJOgDaIdR4QagBfMjHAYeRHAAdIdR4QagB/MRrwLGa04rpcj1OC3fQsaJkUXgMRDxCjReEGiAAWgJO8vCmz3s4mtMScFzNhccWhcdAxCHUeEGoAYLIR9NVLVqP5jRO/plissYymgPYFKHGC0INECJ8HXCMZFmeozlMXQH2QajxglADhKCWgNMrQWo44bOgI3numSN1MHUlEXaAEEeo8YJQA4SJMxQe9xRhBwgvhBovCDVAGPJx4XFn2A0ZCE1hH2pmzpypNWvWaMqUKfrrX//q8VpMTIxycnIkSePGjdPTTz/dpWcSagAb8fNoTmucbQUEV9iHmtWrV6umpkbPP/98u1CTmpqqw4cPd/uZhBrAps44msPUFRDOuvPzOyZAbeqWyZMna82aNcFuBoBw4MjwDA+ODGnYRCnnGp+HHcv931ZfY1wyxT8l7AAhoNuhZu3atXr00UdVUlKiiooKLV++XDNmzPC4Z9GiRXr00UdVUVGh0aNHa8GCBZo4caJPGux0OpWbm6vevXvr4Ycf1pVXXumT5wKwmRALO17rdDjzCvCpboea2tpajRkzRjfffLOuueaadq8vW7ZMBQUFWrRokS6//HI9+eSTysvL0/bt25WdnS1Jys3NVV1dXbuvXbVqldLT0zt9/927dys9PV3btm3T9OnTtXXrVqaTAHTdWYed7jtT2Dl9I6M6gC/0qKbGsqx2IzWXXnqpxo4dq8WLF7uvjRw5UjNmzFBRUVGXn71mzRo98cQT7WpqWsvLy9MvfvELjRs3rt1rdXV1HsHJ6XQqKyuLmhoA3ROgouQurb6SToed1r8m+MDGglZTU19fr5KSEj3wwAMe16dOnar169f3+PnHjh1TQkKC4uLitG/fPm3fvl3Dhw/3em9RUZEeeuihHr8ngAjXMrLjp6mrFmce1Wm5w3j+mlEewM2noebw4cNqbGxUWlqax/W0tDQdPHiwy8+ZNm2aNm/erNraWmVmZmr58uW65JJLtGPHDt1+++2KioqSZVlauHChkpOTvT7jwQcfVGFhofvzlpEaADhrAazTaXF6mqr113oGHxX/tM3dzWEnf6E0dna33xMIV35Z/WRZHrPFMsa0u9aZlStXer0+fvx4bd26tUvPiIuLU1xcXJffEwDOWgDrdDrXJuy8ViCNmMKIDSKGT0NNamqqoqOj243KVFZWthu9AQDbO1PY8XrmlQ/rdUyjDpXvUGpOerf+YQmEK5+GmtjYWOXm5qq4uFgzZ850Xy8uLtY3vvENX74VAISvtmGnm1NYrffDaV1h09YpE6X8Px7QyYRijU5P0qjBSRqd7tDo9CQNH9BX0VEEHdhLt0NNTU2Ndu7c6f68rKxMpaWlSk5OVnZ2tgoLCzVr1iyNGzdOl112mZYsWaLy8nLNnTvXpw0HAFvpxhSW1frXrYJP67DjUpR+n/AfOtyQqlMnGvTPnUf0z51H3I+P7xWl8wclaXR6y4dDFwxKVHyv6AB1GPC9bi/pXrNmjSZPntzu+pw5c7R06VJJTZvv/frXv1ZFRYVycnL0+OOP64orrvBJg88WxyQAsC1vR0UkD5ccGTrZ0KhPP6/R9ooqfXjAqQ8POLWjwqkT9e1reqIsacSAvu6QM6o58PRLiA1wh4DTwv7sJ38g1ABAk0aX0Z4jte6Q8+GBKm0/4NSR2nqv92f0661R7umrJI3OcCjdEU+dDgKCUOMFoQYAOmaMUWV1nTvgtASe8qMnvN7fL6HX6ZDTXKczLLWPYqKjAtxy2B2hxgtCDQB0n/Nkg3a0CjkfHqjSzsoanXK1/9HRuk6nJfBcMChJvWOp08HZI9R4QagBAN+oO9VUp/PhgaY6ne0HnNp+hjqdlvqcUYObRnX696FOB11DqPGCUAMA/uNyGe0+UqvtFadHdbYfqNLhGu91OoMd8U0hp3nqatTgJGX2702dDtoh1HhBqAGAwKt0ntSHFc7mOp2mkZ09R7zX6Th6n67TGdVcqzNiAHU6kY5Q4wWhBgBCQ/XJBu2oqPYoSv60sloNje1/HMXFROmCQYlNq6+aR3VGUqcTUQg1XhBqACB0tdTpbG81qrP9gFO1HdTpDEvt41511bKnTjJ1OrZEqPGCUAMA4cXlMio/esK96qqlVudwTZ3X+911OoNPj+pQpxP+CDVeEGoAwB4qq0+eXnXVHHh2d1CnkxQfo5GtzrwanZGkEQP6qhd1OmGDUOMFoQYA7KulTmd7yzLzCqc++dx7nU5sTJTOT0t0n3s1Kj1JIwcnKSHWp2c8w0cINV4QagAgstSfcunTymp3MXLLfjo1dafa3Wu1qtM5vVNyklL6xgWh5WiNUOMFoQYA4HIZ7T12wh1yWmp1Kqu91+kMSop3bxzYUpRMnU5gEWq8INQAADpyqLqueePA07sklx2u9XpvYnxM82jO6ZPMzxlInY6/EGq8INQAALqjpu6UPqpodZJ5hVOfHKxRfaOr3b0tdTqjBjcVI7ece9UnjjqdniLUeEGoAQD0VP0pl3ZWtjr3qsKpHQecqu6oTielj0a2mroaNThJAxKp0+kOQo0XhBoAgD8YY7T36Bfu0ZyWkZ3Pnd7rdAYmxnlsGjg6PUnZyQnU6XSAUOMFoQYAEEiHa+rcK69aAk/Z4Vp5+6mbGBejkc0bB7YEnnPTqNORCDVeEWoAAMFWW3dKHx10tgo7Tn18sNp7nU50lM5N6+sxqjNycJL6RlidDqHGC0INACAUNTQ21em0HdWpPum9TmdoSp+mAz5bjerYuU6HUOMFoQYAEC6MMdp37AuPM6+2H3DqoPOk1/sHNNfpjGp1JER2coKiosK/TodQ4wWhBgAQ7lrqdFoXJHdUp9M3LkYjBye6p65GDU7SeWmJio0JrzodQo0XhBoAgB2dqD/lPveqJex8dLBa9afa1+n0irZ07sBE95lXo9MdGjk4UYnxvYLQ8q4h1HhBqAEARIpTjS7tOlTrsUPyhweq5PRSpyNJQ1MS3CGnZZn5wMT4ALfaO0KNF4QaAEAkO12n0zR91XKieUWV9zqd1L5xHieZj053aEgQ6nQINV4QagAAaO9obb3H4Z4fHqjSZ2eo02l99tW5aX0VFxPtt/YRarwg1AAA0DUn6k/po4PV7qmr7Qeq9NHBatV1UKdzTnOdzrgh/XXdv2X7tC3d+fkdWTv4AACAM0qIjdHY7P4am93ffe1Uo0ufHW6u09l/egVW1RcN2lHh1I4Kp/YePeHzUNMdhBoAAHBGMdFROi8tUeelJWrml5quGWO0//gX7o0DM/r1Dm4bg/ruAAAgbFmWpcz+Ccrsn6CpowcFuzkKrx14AAAAOkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAtkCoAQAAthAxp3QbYyRJTqczyC0BAABd1fJzu+XneGciJtRUV1dLkrKysoLcEgAA0F3V1dVyOByd3mOZrkQfG3C5XDpw4IASExNlWVawmxMQTqdTWVlZ2rt3r5KSkoLdnICL5P5Hct+lyO5/JPddiuz+27XvxhhVV1crPT1dUVGdV81EzEhNVFSUMjMzg92MoEhKSrLVH/DuiuT+R3LfpcjufyT3XYrs/tux72caoWlBoTAAALAFQg0AALAFQo2NxcXF6Wc/+5ni4uKC3ZSgiOT+R3LfpcjufyT3XYrs/kdy31tETKEwAACwN0ZqAACALRBqAACALRBqAACALRBqAACALRBqwkxRUZEuueQSJSYmauDAgZoxY4Y+/vhjj3uMMfr5z3+u9PR09e7dW5MmTdKHH37ocU9dXZ3uuusupaamqk+fPvr617+uffv2BbIrPVZUVCTLslRQUOC+Zve+79+/X9/97neVkpKihIQEXXzxxSopKXG/btf+nzp1Sj/+8Y81bNgw9e7dW8OHD9f8+fPlcrnc99ip72vXrlV+fr7S09NlWZZWrFjh8bqv+nrs2DHNmjVLDodDDodDs2bN0vHjx/3cu8511veGhgb98Ic/1IUXXqg+ffooPT1ds2fP1oEDBzyeEa59l878e9/a7bffLsuytGDBAo/r4dz/HjMIK9OmTTPPPfec2bZtmyktLTXTp0832dnZpqamxn3PI488YhITE83LL79stm7daq699lozePBg43Q63ffMnTvXZGRkmOLiYrN582YzefJkM2bMGHPq1KlgdKvbNm7caIYOHWouuugic88997iv27nvR48eNUOGDDE33XST2bBhgykrKzN///vfzc6dO9332LX///Vf/2VSUlLM66+/bsrKysxf/vIX07dvX7NgwQL3PXbq+xtvvGHmzZtnXn75ZSPJLF++3ON1X/X1a1/7msnJyTHr168369evNzk5Oebqq68OVDe96qzvx48fN1/5ylfMsmXLzEcffWTeffddc+mll5rc3FyPZ4Rr34058+99i+XLl5sxY8aY9PR08/jjj3u8Fs797ylCTZirrKw0kszbb79tjDHG5XKZQYMGmUceecR9z8mTJ43D4TD/8z//Y4xp+ouhV69e5qWXXnLfs3//fhMVFWXefPPNwHbgLFRXV5tzzz3XFBcXmyuvvNIdauze9x/+8IdmwoQJHb5u5/5Pnz7d3HLLLR7XvvnNb5rvfve7xhh7973tDzZf9XX79u1Gknnvvffc97z77rtGkvnoo4/83Kuu6eyHeouNGzcaSWbPnj3GGPv03ZiO+79v3z6TkZFhtm3bZoYMGeIRauzU/7PB9FOYq6qqkiQlJydLksrKynTw4EFNnTrVfU9cXJyuvPJKrV+/XpJUUlKihoYGj3vS09OVk5PjvieU/cd//IemT5+ur3zlKx7X7d73V199VePGjdO3v/1tDRw4UF/60pf01FNPuV+3c/8nTJigf/zjH/rkk08kSf/617/0zjvv6KqrrpJk77635au+vvvuu3I4HLr00kvd9/z7v/+7HA5HWH0/qqqqZFmW+vXrJ8n+fXe5XJo1a5buv/9+jR49ut3rdu//mUTMgZZ2ZIxRYWGhJkyYoJycHEnSwYMHJUlpaWke96alpWnPnj3ue2JjY9W/f/9297R8fah66aWXtHnzZr3//vvtXrN73z/77DMtXrxYhYWF+tGPfqSNGzfq7rvvVlxcnGbPnm3r/v/whz9UVVWVLrjgAkVHR6uxsVEPP/ywrr/+ekn2/71vzVd9PXjwoAYOHNju+QMHDgyb78fJkyf1wAMP6IYbbnAf4Gj3vv/qV79STEyM7r77bq+v273/Z0KoCWN33nmnPvjgA73zzjvtXrMsy+NzY0y7a2115Z5g2rt3r+655x6tWrVK8fHxHd5nx75LTf9CGzdunH75y19Kkr70pS/pww8/1OLFizV79mz3fXbs/7Jly/Tiiy/qT3/6k0aPHq3S0lIVFBQoPT1dc+bMcd9nx753xBd99XZ/uHw/GhoadN1118nlcmnRokVnvN8OfS8pKdHChQu1efPmbrfTDv3vCqafwtRdd92lV199VatXr1ZmZqb7+qBBgySpXdqurKx0/8tu0KBBqq+v17Fjxzq8JxSVlJSosrJSubm5iomJUUxMjN5++2397ne/U0xMjLvtduy7JA0ePFijRo3yuDZy5EiVl5dLsvfv/f33368HHnhA1113nS688ELNmjVL9957r4qKiiTZu+9t+aqvgwYN0ueff97u+YcOHQr570dDQ4O+853vqKysTMXFxe5RGsnefV+3bp0qKyuVnZ3t/jtwz549+s///E8NHTpUkr373xWEmjBjjNGdd96pV155RW+99ZaGDRvm8fqwYcM0aNAgFRcXu6/V19fr7bff1vjx4yVJubm56tWrl8c9FRUV2rZtm/ueUDRlyhRt3bpVpaWl7o9x48bpxhtvVGlpqYYPH27bvkvS5Zdf3m75/ieffKIhQ4ZIsvfv/YkTJxQV5fnXVXR0tHtJt5373pav+nrZZZepqqpKGzdudN+zYcMGVVVVhfT3oyXQfPrpp/r73/+ulJQUj9ft3PdZs2bpgw8+8Pg7MD09Xffff79Wrlwpyd7975KAlyajR/7f//t/xuFwmDVr1piKigr3x4kTJ9z3PPLII8bhcJhXXnnFbN261Vx//fVel3tmZmaav//972bz5s3my1/+ckgubT2T1qufjLF33zdu3GhiYmLMww8/bD799FPzxz/+0SQkJJgXX3zRfY9d+z9nzhyTkZHhXtL9yiuvmNTUVPODH/zAfY+d+l5dXW22bNlitmzZYiSZxx57zGzZssW9wsdXff3a175mLrroIvPuu++ad99911x44YVBX9bbWd8bGhrM17/+dZOZmWlKS0s9/g6sq6tzPyNc+27MmX/v22q7+smY8O5/TxFqwowkrx/PPfec+x6Xy2V+9rOfmUGDBpm4uDhzxRVXmK1bt3o854svvjB33nmnSU5ONr179zZXX321KS8vD3Bveq5tqLF731977TWTk5Nj4uLizAUXXGCWLFni8bpd++90Os0999xjsrOzTXx8vBk+fLiZN2+exw8yO/V99erVXv8/nzNnjjHGd309cuSIufHGG01iYqJJTEw0N954ozl27FiAeuldZ30vKyvr8O/A1atXu58Rrn035sy/9215CzXh3P+esowxJhAjQgAAAP5ETQ0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALAFQg0AALCF/w/wpXVgg6Q21QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "altitudes = np.linspace(78, 1500, 100)\n",
    "plt.semilogy(altitudes, compute_density(altitudes, pop.champion_x), label ='analytical fit')\n",
    "plt.semilogy(altitude, density, '.', label='data')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "748e5787",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/best_fit_density.pk\", \"wb\") as file:\n",
    "    pkl.dump(pop.champion_x, file)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cascade",
   "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.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "2347a58ecb642ebeea52b8eeec5565393206a0611cdaf781015dedeac8e4b4a2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}