{ "cells": [ { "cell_type": "markdown", "id": "8026ec99", "metadata": {}, "source": [ "\n", "*This notebook contains material from [CBE60499](https://ndcbe.github.io/CBE60499);\n", "content is available [on Github](git@github.com:ndcbe/CBE60499.git).*\n" ] }, { "cell_type": "markdown", "id": "62b667e1", "metadata": {}, "source": [ "\n", "< [2.7 Stochastic Programming](https://ndcbe.github.io/CBE60499/02.07-SP.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [2.9 Supplementary material: data for parmest tutorial](https://ndcbe.github.io/CBE60499/02.09-Parmest-generate-data.html) >
"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0b64e5d5",
"metadata": {},
"outputs": [],
"source": [
"# IMPORT DATA FILES USED BY THIS NOTEBOOK\n",
"import os, requests\n",
"\n",
"file_links = [(\"data/parmest_20210609_data_exp1.csv\", \"https://ndcbe.github.io/CBE60499/data/parmest_20210609_data_exp1.csv\"),\n",
" (\"data/parmest_log_file.csv\", \"https://ndcbe.github.io/CBE60499/data/parmest_log_file.csv\")]\n",
"\n",
"# This cell has been added by nbpages. Run this cell to download data files required for this notebook.\n",
"\n",
"for filepath, fileurl in file_links:\n",
" stem, filename = os.path.split(filepath)\n",
" if stem:\n",
" if not os.path.exists(stem):\n",
" os.mkdir(stem)\n",
" if not os.path.isfile(filepath):\n",
" with open(filepath, 'wb') as f:\n",
" response = requests.get(fileurl)\n",
" f.write(response.content)\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[2.8 Parameter estimation with `parmest`](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8-Parameter-estimation-with-`parmest`)",
"section": "2.8 Parameter estimation with `parmest`"
}
},
"source": [
"# 2.8 Parameter estimation with `parmest`\n",
"\n",
"Created by [Kanishka Ghosh](https://github.com/kanishka-ghosh), [Jialu Wang](https://github.com/jialuw96), and [Prof. Alex Dowling](https://github.com/adowling2/) at the University of Notre Dame."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpages": {
"level": 1,
"link": "[2.8 Parameter estimation with `parmest`](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8-Parameter-estimation-with-`parmest`)",
"section": "2.8 Parameter estimation with `parmest`"
}
},
"outputs": [],
"source": [
"# This code cell installs packages on Colab\n",
"\n",
"import sys\n",
"if \"google.colab\" in sys.modules:\n",
" !wget \"https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py\"\n",
" import helper\n",
" helper.install_idaes()\n",
" helper.install_ipopt()\n",
" helper.download_data(['parmest_20210609_data_exp{:d}.csv'.format(i) for i in range(1,17)])\n",
" helper.download_data(['parmest_log_file.csv'])\n",
" !pyomo build-extensions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbpages": {
"level": 1,
"link": "[2.8 Parameter estimation with `parmest`](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8-Parameter-estimation-with-`parmest`)",
"section": "2.8 Parameter estimation with `parmest`"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"from pyomo.environ import *\n",
"from pyomo.dae import *\n",
"\n",
"# Define the directory to save/read the data files\n",
"data_dir = './data/'"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[2.8.1 What is parameter estimation?](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.1-What-is-parameter-estimation?)",
"section": "2.8.1 What is parameter estimation?"
}
},
"source": [
"## 2.8.1 What is parameter estimation?\n",
"\n",
"Given a function $f(x,\\theta)$ where $x$ is the input or array of inputs, $\\theta$ is the vector of unknown model parameters, and $y$ is the array of observed output, parameter estimation is performed to determine the values of $\\theta$ to minimize the error between $f(x,\\theta)$ and $y$. Commonly, parameter estimation is set up as a least squares objective problem:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\\begin{split}\n",
" \\min_{\\hat{\\theta}} \\quad & \\sum_{i}^{} (y_i - f(x_i,\\hat{\\theta}))^2\\\\\n",
" \\textrm{s.t.} \\quad & \\mathrm{bounds \\ on} \\ \\theta\\\\\n",
" & \\mathrm{other \\ physical \\ constraints}\\\\\n",
"\\end{split}\n",
"\\end{align}\n",
"$$\n",
"\n",
"where $i$ is used to index the datapoints in a dataset and $\\hat{\\theta}$ is the optimal set of parameter values that minimizes the prediction error."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[2.8.2 What is `parmest`?](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.2-What-is-`parmest`?)",
"section": "2.8.2 What is `parmest`?"
}
},
"source": [
"## 2.8.2 What is `parmest`?\n",
"\n",
"`parmest` is a Python package built on the [Pyomo optimization modeling language](http://www.pyomo.org/) to support parameter estimation using experimental data along with confidence regions and subsequent creation of scenarios for [PySP](https://github.com/Pyomo/pysp). `parmest` supports scenario generation for multiple 'experiments' and can be used to characterize estimate uncertainties through, for example, confidence region generations. `parmest` requires the following positional arguments in order solve the optimization problem:\n",
"1. Function that accepts an 'experimental' dataset or a list of 'experimental' datasets, each defined as a dictionary, as it's argument and returns the Pyomo model.\n",
" \n",
" Later in this tutorial, that function is defined above as `create_model()`\n",
"\n",
"\n",
"2. List of datasets where each dataset is a dictionary\n",
"\n",
" Later in this tutorial, the list of datasets is generated using the function `create_data_dict()` (defined below) and is stored in `data_dict_overall`\n",
"\n",
"\n",
"3. List of parameter names (as they appear in the Pyomo model definition) that are being estimated\n",
"\n",
" later in this tutorial, the list of parameter names to be estimated is defined below by `theta_names`\n",
"\n",
"\n",
"4. Optional keyword argument to define the verbosity of solver output. Default: False\n",
"\n",
"More information about the `parmest` package can be found [here](https://pyomo.readthedocs.io/en/stable/contributed_packages/parmest/index.html).\n",
"\n",
"Detailed explanation of the various methods in `parmest` can be found [here](https://www.osti.gov/servlets/purl/1761797).\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[2.8.3 Example: Reaction Kinetics](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3-Example:-Reaction-Kinetics)",
"section": "2.8.3 Example: Reaction Kinetics"
}
},
"source": [
"## 2.8.3 Example: Reaction Kinetics"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[2.8.3 Example: Reaction Kinetics](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3-Example:-Reaction-Kinetics)",
"section": "2.8.3 Example: Reaction Kinetics"
}
},
"source": [
"Consider two chemical reactions that convert molecule $A$ to desired product $B$ and a less valuable side-product $C$.\n",
"\n",
"$$A \\overset{k_1}{\\rightarrow} B \\overset{k_2}{\\rightarrow} C$$\n",
"\n",
"Our ultimate goals is to design a large-scale continous reactor that maximizes the production of $B$. This general sequential reactions problem is widely applicable to CO$_2$ capture and industry more broadly (petrochemicals, pharmasuticals, etc.).\n",
"\n",
"The rate laws for these two chemical reactions are:\n",
"\n",
"$$r_A = -k_1 C_A$$\n",
"\n",
"$$r_B = k_1 C_A - k_2 C_B$$\n",
"\n",
"$$r_C = k_2 C_B$$\n",
"\n",
"Here, $C_A$, $C_B$, and $C_C$ are the concentrations of each species. The rate constants $k_1$ and $k_2$ depend on temperature as follows:\n",
"\n",
"$$k_1 = A_1 \\exp{\\frac{-E_1}{R T}}$$\n",
"\n",
"$$k_2 = A_2 \\exp{\\frac{-E_2}{R T}}$$\n",
"\n",
"$A_1, A_2, E_1$, and $E_2$ are fitted model parameters. $R$ is the ideal-gas constant and $T$ is absolute temperature."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.8.3.1 Batch Reactor](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3.1-Batch-Reactor)",
"section": "2.8.3.1 Batch Reactor"
}
},
"source": [
"### 2.8.3.1 Batch Reactor"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.8.3.1 Batch Reactor](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3.1-Batch-Reactor)",
"section": "2.8.3.1 Batch Reactor"
}
},
"source": [
"The concentrations in a batch reactor evolve with time per the following differential equations:\n",
"\n",
"$$ \\frac{d C_A}{dt} = r_A = -k_1 C_A $$\n",
"\n",
"$$ \\frac{d C_B}{dt} = r_B = k_1 C_A - k_2 C_B $$\n",
"\n",
"$$ \\frac{d C_C}{dt} = r_C = k_2 C_B $$\n",
"\n",
"This is a linear system of differential equations. Assuming the feed is only species $A$, i.e., \n",
"\n",
"$$C_A(t=0) = C_{A0} \\quad C_B(t=0) = 0 \\quad C_C(t=0) = 0$$\n",
"\n",
"leads to the following analytic solution:\n",
"\n",
"$$C_A(t) = C_{A,0} \\exp(-k_1 t)$$\n",
"\n",
"$$C_B(t) = \\frac{k_1}{k_2 - k_1} C_{A,0} \\left[\\exp(-k_1 t) - \\exp(-k_2 t) \\right]$$\n",
"\n",
"$$C_C(t) = C_{A,0} - \\frac{k_2}{k_2 - k_1} C_{A,0} \\exp(-k_1 t) + \\frac{k_1}{k_2 - k_1} \\exp(-k_2 t) C_{A,0} = C_{A,0} - C_{A}(t) - C_{B}(t)$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.8.3.1 Batch Reactor](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3.1-Batch-Reactor)",
"section": "2.8.3.1 Batch Reactor"
}
},
"source": [
"The following Python code simulates and plots this model."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.8.3.1 Batch Reactor](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3.1-Batch-Reactor)",
"section": "2.8.3.1 Batch Reactor"
}
},
"outputs": [],
"source": [
"def kinetics(A, E, T):\n",
" ''' Computes kinetics from Arrhenius equation\n",
" \n",
" Arguments:\n",
" A: pre-exponential factor, [1 / hr]\n",
" E: activation energy, [kJ / mol]\n",
" T: temperature, [K]\n",
" \n",
" Returns:\n",
" k: reaction rate coefficient, [1 / hr]\n",
" \n",
" '''\n",
" R = 8.31446261815324 # J / K / mole\n",
" \n",
" return A * np.exp(-E*1000/(R*T))\n",
"\n",
"def concentrations(t,k,CA0):\n",
" '''\n",
" Returns concentrations at time t\n",
" \n",
" Arguments:\n",
" t: time, [hr]\n",
" k: reaction rate coefficient, [1 / hr]\n",
" CA0: initial concentration of A, [mol / L]\n",
" \n",
" Returns:\n",
" CA, CB, CC: concentrations of A, B, and C at time t, [mol / L]\n",
" '''\n",
" CA = CA0 * np.exp(-k[0]*t);\n",
" CB = k[0]*CA0/(k[1]-k[0]) * (np.exp(-k[0]*t) - np.exp(-k[1]*t));\n",
" CC = CA0 - CA - CB;\n",
" \n",
" return CA, CB, CC"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"nbpages": {
"level": 3,
"link": "[2.8.3.1 Batch Reactor](https://ndcbe.github.io/CBE60499/02.08-Parmest-tutorial.html#2.8.3.1-Batch-Reactor)",
"section": "2.8.3.1 Batch Reactor"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3gU5fbA8e9JCITQQu8QmjRFpEsRvEgVxQJYUVRE7F302sBybVfFn6CoiFzEigVBKaIUlSIgvUhvAYRQA4RAyvv742xICCHZhN1skj2f55kn2ZnZnTOUPfN2cc5hjDEmeIUEOgBjjDGBZYnAGGOCnCUCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgpwlApPvichWEbks0HEUBCIyVkRe8vJc+3MvICwRGL/wfEkcF5GjInJQRH4SkepevjdKRJyIFPJDXANEJMkTV6yILBeRXr6+TrprdhKRaB9/5gDPn9Fb6fZf5dk/1pfXMwWbJQLjT1c454oDlYE9wLsBjifFfE9ckcB7wJciEhngmM4qk4S4Cbgu3fFbgPX+j8oUJJYIjN855+KBb4BGKftE5HIRWep5Kt8hIkPTvOU3z89Dnif3iz3vuVNE1orIERFZIyLN0rynqYisEJHDIvKViIR7EVcy8ClQDKjnuUYREfmviGwXkT0iMkpEinqOlRaRH0UkxlPK+VFEqqW5pzIi8omI7PIcnygixYCpQBXPvRwVkSqe6wz3nLvL83sRz+d0EpFoERkiIv8An5zlFv4BVgLdUq4PtAUmpT1JRK4UkdUickhEZotIwzTHLhKRJZ4/06+A8HTv7SUiyzzvnSciTbL6czX5jyUC43ciEgFcByxIs/sY+vQaCVwO3C0iV3mOXeL5GemcK+6cmy8ifYGhnveUBK4E9qf5vH5Ad6AW0AQY4EVcocBtQAKwzbP7NeA8oClQF6gKPOc5FoJ+KdcEagDHgRFpPvJTIAJoDFQA3nbOHQN6ALs891LcObcLeBpo47nOhUAr4Jk0n1UJKOO51qBMbmMc+mcCcD3wA3AizT2eB3wBPASUB6YAk0WksIgUBiZ64i4DTACuTfPeZsAY4C6gLPABMCklYZkCxDlnm20+34CtwFHgEJAI7AIuyOT84egXJ0AU4IBCaY5PBx7M5Fo3p3n9OjDqLOcO8MRzCE0Ax4F+nmOCJqg6ac6/GNhyls9qChz0/F4ZSAZKZ3BeJyA63b5NQM80r7sBW9OcfxIIz+TPawDwB1AUrXYrhSbadsBLwFjPec8CX6d5Xwiw03ONSzx/L5Lm+DzgJc/v7wMvprvuOqBjmj/3ywL9b822c9+sRGD86SrnXCRQBLgPmCMilQBEpLWIzPJUsxwGBgPlMvms6uiX59n8k+b3OKB4Jucu8MRVGq1G6eDZXx59ov/LUxVyCJjm2Y+IRIjIByKyTURi0SqsSE/JojpwwDl3MJPrplWF1FIInt+rpHkd47RKLVPOuePAT2hpopxzbm5m13FaHbYDLelUAXY6z7d6mjhS1AQeTfmz8Px5VE8XpykALBEYv3POJTnnvgOSgPae3Z+jX8LVnXOlgFHoEzloaSC9HUAdH8d1FLgH6C8iFwH70BJCY+dcpGcr5bRhGeBRoD7Q2jlXktQqLPHEV+Ysjc4Z3c8u9Is2RQ3PvszeczbjPLF9mtV1RETQL/OdwG6gqmdf2jhS7ABeTvNnEemci3DOfZGN2Ew+YInA+J2o3ugT+FrP7hLoE3S8iLQCbkzzlhi0mqV2mn2jgcdEpLnn8+qKSNov0hxxzu33fPZznqflj4C3RaSCJ/aqItItTczH0UbsMsDzaT5nN9oo/J6nUTlMRFISxR6grIiUSnPpL4BnRKS8iJRD2yHG5/A25gBdyLhX1tfA5SLSWUTC0IRxAq0Cmo9Wkz0gIoVE5Bq0rSLFR8BgT+lNRKSYp5G/RA7jNHmUJQLjT5NF5CgQC7wM3OqcW+05dg/wgogcQb8Ev055k3MuznP+XE+VRBvn3ATPvs+BI2gjZxkfxTkc6OnpETME2Ags8FT//IKWAlLOK4qWHBag1UZp9UfbHf4G9qINtDjn/ka/+Dd77qcKWo+/GFiB9vxZ4tmXbU796pw7kMGxdcDNaJLYB1yBdus96Zw7CVyDtjccRBv0v0vz3sXAnWiD+EHPn8uAnMRo8jY5vXrQGGNMsLESgTHGBDlLBMYYE+QsERhjTJCzRGCMMUHO57M7+lu5cuVcVFRUoMMwxph85a+//trnnCuf0bF8lwiioqJYvHhxoMMwxph8RUS2ne2YVQ0ZY0yQs0RgjDFBzhKBMcYEuXzXRmCMMf6SkJBAdHQ08fFZTvyaZ4WHh1OtWjXCwsK8fo8lAmOM8YiOjqZEiRJERUVx+qSs+YNzjv379xMdHU2tWrW8fp/fqoZEZIyI7BWRVWc5LiLyfyKyUXSJwWYZnWeMMbklPj6esmXL5sskACAilC1bNtslGn+2EYxFlw48mx7oOrH10KX43vdjLMYY45X8mgRS5CR+vyUC59xvwBnT4qbRGxjnmUJ3AbrSU2V/xbNpEzz0ECQk+OsKxhiTPwWy11BVdAWkFNGefWcQkUEislhEFsfExOToYmvXwjvvwJgxOXq7McYUWIFMBBmVXzJcHME596FzroVzrkX58hmOkM7S5ZdDu3bwwgtw/HiOPsIYYwqkQCaCaHTt1BTVOH3NVp8SgVdegV27YMQIf13FGGN8Y8KECbRu3ZomTZpQt25dhg0b5rdrBTIRTAJu8fQeagMc9qz76jcdOkCPHpoQDh3y55WMMSbn/ve///Haa6/x7bffsmLFCpYtW0ZERITfrufP7qNfoItj1xeRaBG5Q0QGi8hgzylTgM3oOqgfoWvY+t1//gMHD8J//5sbVzPGmOyJjY3lkUce4euvv6ZatWoAFC9enMcff9xv1/TbgDLn3A1ZHHfAvf66/tk0bQrXXw9vvw333QeVKuV2BMaY/OChh2DZMt9+ZtOmMHx45ud8//33tG7dmtq1a/v24pkIyrmGXngBTpyAl18OdCTGGHO61atX07Rp07Mev+6663jzzTd9es2gnGKiXj0YOBA++AAeeQSyMRLbGBMksnpy95dixYpx/CxdG3/44Qd69erFL7/84tNrBmWJAODZZyE0FJ5/PtCRGGNMqp49ezJhwgT27NkDwIkTJ/joo4+Ij49nwoQJ9O/fn8OHD/v0mkGbCKpWhfvvh/HjYVWGsyEZY0zua9myJUOHDqVbt25ccMEFNG3alL179/LGG29w9OhRBg8ezOrVq89aasiJoE0EAE8+CSVLwtNPBzoSY4xJ1b9/f5YtW8bKlStZu3Yt/fv3Z+vWrUycOJFRo0Zx8803s2LFCp9dL6gTQZky8PjjMGkSzJ0b6GiMMSZjNWrU4OOPPz71+vnnn6d169Y++/ygTgSgXcSqVIFHHwWX4QQXxhhTsAV9IihWTLuR/vknfPVVoKMxxpjcF/SJAOCWW3Sgx5AhNiGdMSb4WCIAQkLgzTdh+3adqtoYY4KJJQKPf/0LrrxS5yLauzfQ0RhjTO6xRJDG669r1ZANMjPGBBNLBGnUrw933w0ffgirVwc6GmOMyR2WCNJ5/nkdZPbYY4GOxBhjcoclgnTKltV5iKZNg+nTAx2NMcb4nyWCDNx7L9Spo4PMEhMDHY0xJhgFy1KVeVaRIvDaa9pO8NFHgY7GGBNscnupSnH5bF6FFi1auMWLF/v9Os5B5866QtH69VCunN8vaYwJsLVr19KwYcOAxhAbG0utWrVYtGhRjlcpy+g+ROQv51yLjM4PyoVpvCEC776rI47//W/tSWSMCS6dxnbK8pxe5/XisbaPnTp/QNMBDGg6gH1x++jzdZ/Tzp09YHaWn2dLVeYxjRvDAw/A6NGwaFGgozHGBIPMlqps0aIF9957Lx07dmS1D/u4W4kgC88/D59/rg3ICxbodBTGmODgzRP82c4vF1Eu2++Hsy9VuWPHDlq1asXIkSN56623iI6OpnHjxtn+/IzY11oWSpaEN97QEsGYMYGOxhhT0J1tqcq//vqL9evXc/vttzNz5ky6devms2taIvDCTTdB+/a6otmBA4GOxhhTkJ1tqcq//vqLN998kzFjxlCkSBGOHTvms2ta1ZAXRGDECGjWTAebjRwZ6IiMMQVZ//796d+//2n7evbsyb59+wgJCaFhw4YUK1bMZ9ezROClCy+Ee+6B996DgQPhoosCHZExJphMmTLFb59tVUPZ8OKLOgXFvfdCcnKgozHGGN+wRJANkZE64nj+fPjf/wIdjTHG+IYlgmy69VZo21ZnJ42JCXQ0xhhz7iwRZFNIiM4/dOQIPPJIoKMxxphzZ4kgBxo1gqeegvHj4eefAx2NMcacG0sEOfTUU7qi2eDB4MPuvMYYk+v8mghEpLuIrBORjSLyZAbHS4nIZBFZLiKrReQ2f8bjS+Hh8MEHsGUL+HGacGOM8Tu/JQIRCQVGAj2ARsANItIo3Wn3AmuccxcCnYA3RaSwv2LytY4ddUzBW2/B0qWBjsYYY3LGnyWCVsBG59xm59xJ4Eugd7pzHFBCRAQoDhwA8tWaYK+/rmsV3HknJCUFOhpjjMk+fyaCqsCONK+jPfvSGgE0BHYBK4EHnXNnDNUSkUEislhEFsfksT6bpUvDO+/AX3/p+gXGGOMLBWWpSslgX/rl0LoBy4AqQFNghIiUPONNzn3onGvhnGtRvnx530d6jvr1g5494ZlnYNu2QEdjjMnvcnupSn/ONRQNVE/zuhr65J/WbcCrTtfL3CgiW4AGwEI/xuVzIjoHUePGMGgQTJum+4wx+VynTmfu69dPJx6Li9MnwPQGDNBt3z7oc/oKZcyeneUlY2NjeeSRR1i0aBHVqlUDoHjx4jz++OPZjd5r/iwRLALqiUgtTwPw9cCkdOdsBzoDiEhFoD6w2Y8x+U3Nmtpe8PPPuqKZMcbkRCCWqvRbicA5lygi9wHTgVBgjHNutYgM9hwfBbwIjBWRlWhV0hDn3D5/xeRvgwfDd9/piOOuXTU5GGPyscye4CMiMj9erpxXJYD0Mluq8quvvmLu3LkkJydTrFgxXnvttWx/fkb8Og21c24KMCXdvlFpft8FdPVnDLkpJAQ+/hjOPx9uvx1mzLClLY0x2XO2pSrnz5/Pn3/+yf/93/8BcPLkSZ9dM9OvKREp48UW6bNoCoCaNXVcwcyZOuDMGGOy42xLVY4dO5aHHnro1HmFC/tuyFVWJYJdni2zps9QoIbPIioABg6Eb76Bxx+Hbt0gF6v6jDH5XNqlKpOSkkhMTOTmm28mPj6eQoVSv7KTkpIIDQ31yTWzSgRrnXOZrsUlIjamNh0RbTBOqSKaOdOqiIwx3stoqcrVq1fz6KOPUr58eY4cOcLbb79NZKRvKmSySgQXe/EZ3pwTdKpXh+HDNRGMGAEPPBDoiIwx+Vnjxo354osv/PLZmT6nOufiz3ZMRLZndU6wGzBAuxk/+SRs2BDoaIwxJmPnUmFhQ6ayIKKL2BQpAv37Q0JCoCMyxpgznUsiSD9dhMlAlSowahT8+Se8+GKgozHGmDNl2kYgImdbjDFltlDjheuug6lT4eWXoUsX6NAh0BEZY87GOYfk4zlidMae7MmqRFDiLFtx4J1sXy2IvfsuREXBzTfDoUOBjsYYk5Hw8HD279+foy/TvMA5x/79+wkPD8/W+7LqNbQBmO6c25/jyAwAJUrA559Du3Y6X9Vnn9nEdMbkNdWqVSM6Opq8Nt19doSHh5+arM5bWSWCGsAEEQkDfgWmAgtdfk2XAda6NQwdCs8+Cz16aAOyMSbvCAsLo1atWoEOI9eJN9/pIlICuAzojq48thaYhpYW9vg1wnRatGjhFi9enJuX9KmkJJ3ZdvlyWLbMRh0bY3KHiPzlnGuR0TGveg0554445753zt3lGWn8ElAeGOfDOINCaCiMH68jjW+6CRLz1cKcxpiCKKtJ55pltAHhwCznXLfcCbNgqVlTu5QuWAAvvBDoaIwxwS6rNoI3MznmgH/5MJagcv31upLZSy9pd9IuXQIdkTEmWGWaCJxzl+ZWIMFo5EhYvFiriJYuhapVAx2RMSYYedVGICJhIvKAiHzj2e7z9CQy56BYMZgwQZc+vf56ay8wxgSGt1NMvA80B97zbM09+8w5atgQPvwQ/vgDnnkm0NEYY4KRt0tVtnTOXZjm9UwRWe6PgILRjTfCb7/Ba69B+/bQq1egIzLGBBNvSwRJIlIn5YWI1AaS/BNScBo+HJo2hVtugW3bAh2NMSaYeJsIHgdmichsEZkDzAQe9V9YwSc8XNsLkpKgXz/w4brUxhiTKW8HlP0K1AMe8Gz1nXOz/BlYMKpbF8aMgYUL4bHHAh2NMSZYeNVGICKhQDcgyvOeziKCc+4tP8YWlK69Fh56SKuKWrTQqiJjjPEnbxuLJwPxwEog2X/hGIDXX9e5iAYNgkaNNCEYY4y/eJsIqjnnmvg1EnNKWBh89RW0bAlXX62DzipWDHRUJrc454hLiOPwicMcij902nbkxBHiEuKoXbo2vRv0BuD+KfdzcfWLufGCG4lPjKfvhL4kJCVwMukkSS7p1Nz6Dnfq9z6N+vBQm4c4kXiC3l/25o6L7qBv477sPbaXh6c/THhoOOGFTt+KFy5OscLFKF64OM0rN6d+ufqcTDrJpgObqF6qOsUL21pV+ZW3iWCqiHR1zv3s12jMKeXLw/ff6/oFffrAr79C4cKBjsqci/jEeHbG7mTnkZ0kJSdxaS0duP/EjCcoU7QMT7Z/EoAK/63Avrh9mX7WNQ2vOZUIZm2dRdmIsgAIws7YnYSFhlE4tDChEkpISAjiWWJcRBCEwqH6jynZJXMo/hDxifEAxCXE8Wf0n8QnxnMi6QTxifHEJ8aTmHz6aMc3urxB/XL12XZoG43ea8SnV3/KzU1u5o/tf9D7y95EhkdSpmiZ1C1cf5aNKEu5iHJ0rNmR6qWqk5CUQJJLIrxQ9hZSMb7lbSJYAHwvIiFAArpUpXPOlfRbZIaLLoKPP9ZxBg8/rFNSmLxv/o75rNy7kk0HNrHp4CY2H9zMjtgdp325X1DhAlbcvQKAzQc3c+zksVPHHrv4MUSEUkVKERkeeWorFV6KkkVKEhEWQbGwYqfOX3XPqlO/FylUhCV3LfE61qJhRVkwcMGp11GRUWx8YOMZ5yUmJ3Ls5DGOnjzKsYRjlClaBoAKxSrwxbVf0KZaGwDKRZTj+sbXc/jEYQ4cP8CB4wfYcnAL+4/v5+DxgzjPUuff9vuW6qWqM2vrLLqN78bc2+fStnpbZmyawYdLPqRSsUpUKq5b5RKVqVqiKlVLVqVs0bL5ehnJvMrb9Qg2A1cBKwO9KE1+X48gJ554At54Az76CAYODHQ05kTiCTYd3ESj8o0AeHPem8zeNpvJN0wGoOunXZmxeQZhIWHUKl2L2qVrU7NUTaqVrHZqq1GqBueVPS+Qt5HrUkof++L2Ual4JUoWKcmG/Rv4evXXDGw2kIrFK/LVqq8YOmco/xz9h0PxZ67pWji0MFVKVGHaTdOoX64+C3cuZN6OedzV/C6KhhUlLiGO8ELhhIi3PeODR2brEXhbItgArAp0EghWr7yijcf33KONx23bBjqi4BF7Ipbl/yxnye4lLPlnCUt3L2VNzBqSXBKHhhyiVHgpwkLDCAsJO7Xo+YieIwgvFE7VElUJDQkN9C3kGSEScqqqKEW9svV4+pKnT72+7vzruO786wCtSttzdA+7juxi15Fd7Dyy89TPlKqw6Run8/zs57mn5T0ADJkxhA/++oDqpapTs1RNakbW1J+lahIVGUWt0rWoVrIahUK8/eoLDt6WCMYCtdGlKk+k7A9E99FgLBEAHDigjcfHjuk4gxo1Ah1RwTVvxzxGLR7FgugFbDiw4dT+isUq0rxKc5pWbMr5Fc7nyvpXUqxwsUw+yfibc46D8QdPJZepG6YyZ9scth3exrZD29h2eBu7j+w+VSUFEBEWwZGnjhAiIXy6/FMOxh/kgdYPAHDs5DEiwiIKZPWTL0oEWzxbYc/m7YW7A+8AocBo59yrGZzTCRgOhAH7nHMdvf38YFKmDEyeDBdfrHMR/fEHlLQWGp9YE7OGh6Y9xKuXvUqzys3YdWQXP2/6mTbV2nDLhbfQrHIzLqp0EZVLVA50qCYdETmthNGjXg961Otx2jknk06y4/AOth7aypZDWzgUf+hU1dHk9ZPZfnj7qUTQcWxHNhzYQJ3Sdahbpu6prV6ZetQtU5dKxSsVzCThr9oezyC09UAXIBpYBNzgnFuT5pxIYB7Q3Tm3XUQqOOf2Zva5wVoiSDFjhi5836WLJoZCVsL1WlJyEn/t/otfN//KzK0z6dOwD3e1uItdR3bR87OevNn1TTrX7kxSchIhElIg/8ObMyUmJ56qKvpg8Qes3LuSzQc3s/HARrYc2nJaj6liYcXo27gvn/T+BIBJ6yZRt0zdU+1FeVmOSwQiMtQ5NzSH57QCNjrnNnvO+xLoDaxJc86NwHfOue0AWSUBowng/fd1sNmDD8KIEWDfV2cXcyyGKRum8OOGH5mxaQaHTxwG4PwK55/6z1+lRBWWDV526j1Wrx9c0rYX3NXirtOOJSYnsv3wdjbs38DGAxvZcGADUZFRgFZL3fDtDQxqNoi3u7/NicQTdB7XmXpl61G/bH0alGtA/bL1qVOmzqnuunlVVs+TA0UkNpPjAlwPDM3gWFVgR5rX0UDrdOecB4SJyGygBPCOc27cGRcRGQQMAqhhlePceSds2KA9ierV0ykpzOnenv82E9ZMYEH0AhyOysUr06dRHy6rfRmXRl1KxeI2Qs9krVBIIWqXrk3t0rXpxplLtC+9aylFQosAcDD+IGGhYUzfOJ2xy8aeOidUQqldujYNyjWgYbmGXNvoWlpVbZVbt+CVrBLBR+gXdFbnZCSj59T09VCF0EVuOgNFgfkissA5t/60Nzn3IfAhaNVQFvEEhVdfhY0b4ZFHoE4duOKKQEcUWGti1jBryyzubXUvADO3ziQhOYGhnYbS67xeXFTpIqvqMT4lIqd1Aa5UvBKzbtW5OGNPxLJ+/3rW7VvH3/v+5u/9f7M2Zi3TNk6jbpm6tKraipV7VtJtfDc+u+YzLq11KbuO7GLdvnU0Kt+ICsUq5Oq/16zWLB52Dp8dDVRP87oasCuDc/Y5544Bx0TkN+BCtG3BZCIkBMaPh44d4YYb4PffdQBasHDOsWrvKuqUqUNEWAQT/57I87Ofp1/jfpQvVp7v+n1HWKitpmoCo2SRkrSo0oIWVU6vkk9MTiQpWZdyCQsNo1vdblQrWQ2AKRumcOfkOwEoU7QMjcs3plH5Rqk/KzSmUvFKfonXn43FhdAv9M7ATrSx+Ebn3Oo05zQERqAzmxYGFgLXO+dWnfmJKtgbi9PbvRtat9b1jufPh5o1Ax2Rf+09tpfxK8bzybJPWLV3FV/3+Zq+jfty4PgBkl0y5SLKBTpEY3Jkf9x+luxewpqYNayJWcPqmNWsjll9amBd2tHoOeGL7qPZ5pxLFJH7gOlo99ExzrnVIjLYc3yUc26tiEwDVqCzmo7OLAmYM1WuDFOmQIcO0K2bdistV8C+CxOSEpi6cSqfLPuEH9f/SGJyIq2rtmZkz5F0iuoEcFoXQmPyo7IRZelSpwtd6nQ5tc85xz9H/2FNzBpOJvlvtSq/lQj8xUoEGfv9d+jaFZo00QnqiheAiSC3HdrG+4vfZ+yysew5toeKxSrSv0l/brvotnzRXc+YvOScSwQiUh64k9SFaQBwzt3uiwDNuevQAb78Eq65RmcrnTxZp7POz/pM6MPS3Uu5ov4V3N70drrX7W71/sb4gbdVQz8AvwO/YIvW51m9e8MHH2j30ttug3HjtFE5v5i1ZRZP/vok02+eTmR4JO9f/j4Vi1WkeqnqWb/ZGJNj3iaCCOfcEL9GYnxi4EDYuxeefhoqVIA338zbA872xe3jROIJqpasSqnwUiQmJ7LryC4iwyPP6HFhjPEPb58XfxSRnn6NxPjMU0/B/ffD22/roLO8KDo2moenPUzN4TV58lddkKVZ5WYsvnOx1f8bk8u8LRE8CPxbRE6iC9OALUyTZ4nA8OFaMhgyRCenGzw40FGpDfs38Nrc1xi3fBzJLpmbmtzEkHaphU0b9GVM7vMqETjnshpdbPKYkBBtIzh6FO6+G4oWhVtvDVw82w9v54U5LzB22VjCQsMY1HwQj7V97NS8LcaYwPF6HIGIXAlc4nk52zn3o39CMr5SuDB8841OP3H77RAeDtddl7sx7Ivbx0u/vcT7i98H4P5W9/Nk+ydtrh9j8hBvu4++CrQEPvPselBE2jvnnvRbZMYnwsNh4kSduvqmm6BIEbjqqty7/rGTx/hoyUf0b9Kf5zo+R41SNmmgMXmNtyuUrQCaOueSPa9DgaXOuSZ+ju8MNqAsZ44c0Smsly6FH36A7t39d62xy8by86af+fzazwEdOp+ytKAxJjAyG1CWnV7mkWl+L3VuIZncVqIETJ2qax5ffTXMnOm/ax04foDdR3cTlxAHYEnAmDzO20TwCrBURMaKyP+Av4D/+C8s4w+lS+sKZynTVs+Z45vP3XZoG30n9OXT5Z8C8GDrB5l5y0wiwiJ8cwFjjF95lQicc18AbYDvPNvFzrkv/RmY8Y9y5eCXXyAqStsNfv01558VnxjP0NlDaTCyAT+t/+nU6l+hIaHWDdSYfCTTRCAiDTw/mwGV0fUDdgBVPPtMPlSpEsyaBXXrQq9eMH169j9jztY5XDjqQobNGcbVDa5m3X3ruK/Vfb4P1hjjd1n1GnoEXSLyzQyOOeBfPo/I5IoKFbSdoEsXuPJK+PZbTQpZORR/iCdmPMFHSz6idunazOg/g8tqX+b/gI0xfpPVCmWDPL/2cM7Fpz0mIuF+i8rkinLltGqoWzedtfTrrzPvWvrd2u+4b8p97Dm2h8fbPs7QTkOtHcCYAsDbxuJ5Xu4z+UyZMtpm0Lw59O0LEyZkfF5iciLPz36eSsUrsejORbze5XVLAsYUEJmWCESkElAVKCoiF5G6IH1JwL4FCohSpbSd4PLLdf3juBcPSXEAACAASURBVLjU6SimbphKuxrtKFmkJFNvmkql4pUoFOK3he2MMQGQ1f/obsAAdOH5t9LsPwL8208xmQAoWVLHGVxzDQwYAPv3w1W3babXF714psMzDLt02KlFto0xBUtWbQT/A/4nItc6577NpZhMgBQvriubXXtbNI8+Wo2YmNr8dPsULq3VKdChGWP8yNvZR78VkcuBxkB4mv0v+Cswk/sSkxN5fcEr/NzwRXrdP41XX/0X+/Z1o8uoQEdmjPEnbyedG4W2CVwKjAb6AAv9GJfJZXuO7uGGb29g1tZZ3HjBjYzo0Yy3SsFLL8GBA/DZZzqBnTGm4PG211Bb59wtwEHn3DDgYsAWki0g/tj+B80+bMb86PmMuXIMn13zGaWLRvLii7rK2XffaUNybGygIzXG+IO3iSBlDEGciFRBVymr5Z+QTG5xzvH2/LfpNLYTRQsVZf4d87ntottOO+ehh3SBm99+g/btITo6QMEaY/zG20QwWUQigTeAJcBW4At/BWX8L/ZELP2+6ccjPz/CFfWvYPGgxTSt1DTDc/v3h59+gq1boXVrWLYsd2M1xvhXlolAREKAX51zhzw9h2oCDZxzz/k9OuM3vb/szfdrv+f1y17nu37fERkemen5XbvC3LkQGgodOmhXU2NMweDtwjTznXMX50I8WbKFaXxj7va5JLkkLql5SdYnp7Frl85JtGIFjBwJd93lpwCNMT7li4VpfhaRa8XmFs7XRi4cyfOzngegXY122U4CAFWqaHtBt24weDAMGQLJyb6O1BiTm7ydK+ARoBiQKCLx6FQTzjlX0m+RGZ9yzrHsn2XsjdtLUnISoSGhOf6s4sV1ucsHHoDXX4f167VBuUQJHwZsjMk13g4os//i+dTh+MPExMVQt0xd3rv8PUIk5JySQIpChbRqqH59eOQRaNtWk0Pt2j4I2hiTq7yqGhKRM9axymifyVu2HtpK2zFt6fV5LxKTEwkLDfNJEkghAg8+qBPW7dwJLVv6dy1kY4x/ZLVCWbiIlAHKiUhpESnj2aKAKrkRoMmZlXtW0vbjtuw+spv3L3/frzOGXnYZLFyoK5917QojRoAXfRCMMXlEViWCu9CF6ht4fqZsPwAjs/pwEekuIutEZKOIPJnJeS1FJElE+ngfujmbudvncsnYSwiREH6/7XcurXWp369Zty7Mnw89e8L998OgQXDypN8va4zxgUwTgXPuHedcLeAx51xt51wtz3ahc25EZu8VkVA0WfQAGgE3iEijs5z3GpCDlXNNej+t/4kun3ahQrEKzL19Lo0rNM61a5csCRMnwtNPw+jRcMklsGNHrl3eGJNDXrUROOfeFZG2InKjiNySsmXxtlbARufcZufcSeBLoHcG590PfAvszVbk5gzjV4yn95e9aVS+EX/c9gc1I2vmegwhITpR3TffwJo1cNFF8PPPuR6GMSYbvG0s/hT4L9AeaOnZMhyYkEZVIO3zYLRnX9rPrQpcDWQ60bGIDBKRxSKyOCYmxpuQg87wBcPp/31/OkZ1ZOatMylfrHxA47n2Wli8WMcddO8Ow4bZeANj8ipvWxBbAI2cN8OQU2U0+Cz9+4cDQ5xzSZmNVXPOfQh8CDqyOBsxBI3th7dzTcNr+OyazwgvlDfmiz7vPFiwAO6+G4YO1TaE8eOhXLlAR2aMScvbRLAKqATszsZnR3P6VNXVgF3pzmkBfOlJAuWAniKS6JybmI3rBLV9cfsoF1GO/3b9L845n3YP9YWICBg7Vmcuvf9+aNYMvv4a2rQJdGTGmBTeTjFRDlgjItNFZFLKlsV7FgH1RKSWiBQGrgdOe4+n4TnKORcFfAPcY0nAe2/MfYMm7zchOjbaZwPF/EEE7rwT5s3TgWjt28N//gNJSYGOzBgD3pcIhmb3g51ziSJyH9obKBQY45xbLSKDPcdtAcRz1LNeT/Yc20OVEvljSEezZrB0qc5R9PTT8Msv8OmnULVq1u81xviPV7OPAohITaCec+4XEYkAQp1zR/waXQZs9lH4ZfMvdK7Vmfw6B6Bz8L//wX33QZEiMGYM9M6oP5kxRv/DrF8Pe/Zon+wcOufZR0XkTrTq5gPPrqqAVeEEwNDZQ+nyaRe+XfttoEPJMREYMACWLIGoKLjqKrj3Xjh+PNCRGZPHPPYY1KgBDRr4dc53b9sI7gXaAbEAzrkNQAV/BWUy9uKcFxk2Zxi3Nb2NaxpeE+hwztl552m7waOPwnvvQfPm2uXUmKBz7Jiu9vTIIzo8P8Xx43DxxfDBB/Djj367vLdtBCeccydTqiJEpBBndgU1fjRi4Qiem/0ct154K6OvHE2IeJvD87YiReC//9U5im6/XXsTPfOMtiGEhQU6OmP87Kef4I039IkoIUH/Q7RvD0eP6nzvI7OcyccnvP02mSMi/waKikgXYAIw2X9hmbS+WPkFD0x9gN71exeoJJBW166wciXccIMOPmvTBlavDnRUxvjQ5s36ZN+nD6xdq/uOHYPYWHjoIR2Cf/Cg9qIoXjxXQ/N2qcoQ4A6gKzpQbDowOpsDzHwi2BqLp22cxhVfXEG76u2YetNUioYVDXRIfvfdd9qzKDZWp6t4+GFdK9mYfGf3bnjhBZgxAzZt0n3Vq+tkXF275moomTUWe5sIigHxzrkkz+tQoIhzLs6nkXohmBLB/B3zuezTyziv7HnMvnU2pcJLBTqkXLN3r7aNTZyoi96MHg0NGwY6KmMykZgIixbpk32dOnDzzXD4sPaI6NABunTRrX597TGRy3yxZvGvQNpH0aLAL+camDm7ZJfMwMkDqVKiCtNumhZUSQCgQgUtGYwbB3//DU2b6oOVTW1t8pxPPtHJtcqV06eWYcO0zh+gVCnYtw8mTdKh9Q0aBCQJZMXbEsEy51zTrPblhmAqEWw+uBlBqFW6VqBDCai9e7UK9YsvoHFjLR3YFBUmII4cgdmz9enk8cd132WXaT//rl2hWzfo3BnKlAlomBnxRYngmIg0S/OBzQHr9e0HB48f5M15b5LskqldunbQJwHQ0sHnn8PkyVrSbttWl8g8ejTQkZmgsH49vPIKdOqkX/BXXgkvvghxnprx776Dbdv0CaVv3zyZBLLibSJ4CJggIr+LyO/AV8B9/gsreI1fMZ6nfn2K1Xuty0x6vXrpGgf33gvvvqttBt9+a8tiGh/bs0fnPtm3T19Pmwb//rc+hTz6KPz6K8TE6IyKoCsy5cHqnuzIzhQTYUB9tNfQ3865BH8GdjYFvWrIOcfafWtpVP6MxdxMGvPm6fTWK1Zoafzdd6FevUBHZfKlhASYOxemT9dt6VLd//nn2p/54EE4cUIX5c7HzrnXkOdD2gJRpBmE5pwb54sAs6OgJoLRS0ZzcbWLc3VpyfwuMVHH2zz7rP4/HTIEnnoKihb8HrbmXG3ZAvHxWqzcuhVq1dKpcdu21SeL7t21h0JIwRmz44u5hnKyQpnx0sS/JzJo8iDenP9moEPJVwoV0raCdet0jM6LL2pj8uTJVl1k0omL0ykcHnxQu2/Wrg3PP6/HoqL02P79MGeOVgM1a1agkkBWvO01tJbsr1DmFwWtRLBk9xI6fNKBxuUbM3vAbCLCIgIdUr41e7a2H6xZo92133oLzj8/0FGZgHBOB3NV8UzR3rKlTmQVHg6XXqpP/T166IRXQSKzEoE/VygzWdgZu5MrvriCskXLMumGSZYEzlGnTrBsGbz/vi6NeeGFOiht2DAoH9glnE1uOHIEZs3Sp/tp07Sxd/9+KFwYnntOf15yidUdZsDbRJCyQtlC4ETKTufclX6JKgjEJ8Zz1VdXEXsilrm3z6VS8fzdEJVXhIXBAw/ATTdpAnjvPW3ze/ZZHc9TuHCgIzQ+45xuISHadfOee7Tht3hx7dvfvXvqMnhXXBHYWPM4b6uGOma03zk3x+cRZaGgVA0NnDSQj5d+zMTrJtK7ga3K4i9r1+rMvtOm6aj/l1/Wrt5BVP1bsMTG6qRsKU/9H3+sA7mWLIGvvtIv/3btLONn4Jwbiz1f+H8DJTzb2kAkgYJi9JLRfLz0Y57u8LQlAT9r2FC/M6ZM0RqB66+HVq30u8TkI3v2aN1+2bI6ncPXX+tfZMmSerxZM3jtNT3HkkC2edtrqB+wEOgL9AP+FJE+/gysoNp2aBv3TrmXLrW7MKzTsECHEzR69ND2g7FjdSxQly6pD5ImjzlyRGcbHDRIJ5gCncfHOR3QNXu21v9/+63NNeIj3lYNLQe6OOf2el6XB35xzl3o5/jOUBCqhj5f+Tld63SlXES5QIcSlOLjtUH55Ze1LfG667RxuUGDQEcW5EaP1gmlfv9d6/pLlIDbboN33gl0ZAWCL+YaCklJAh77s/FeAyQlJ7Fu3zoAbrzgRksCARQermscbNqkK6H9+KOOP7j5Zp1WxuSCY8d0wMeTT6YO+vjjD51h8OGHtffPvn2WBHKJtyWCN4AmwBeeXdcBK51zT/gxtgzl1xLBK7+/wgu/vcCqu1dRp0ydQIdj0oiJ0dUCR47U0sLNN2svo7p1Ax1ZARMdDd98ow02c+bonOLFiulMntWqaSnA1if1G180Fj8OfIAmgwuBDwORBPKz2y+6nVc6v2JJIA8qXx5ef11XEnzoIW2HbNBAayU2bAh0dPnY8ePas2f7dn3955/6tB8drX15f/lF6+aqVdPjlgQCJtMSgYjUBSo65+am238JsNM5t8nP8Z0hv5UI9sftJzI8ktAQW2sxv/jnH+2AMmqUPrT26aNzGDXN9dU38qFt2/SJ/6efYOZMTQavvKJVQHFx2vunlk2tHgjnUiIYDhzJYH+c55jJRGJyIld+eSVXf3V1oEMx2VCpErz9ts5F9sQT2v30oovg8st1kkqTRkIC7Nihv8fF6ZQN99yjAzgGDkyd3wd02mZLAnlSVokgyjm3Iv1O59xidCZSk4kX57zIvB3zuOH8GwIdismBihX1YXb7du1htHAhtG+vsxRMngzJyYGOMED++Uf74fbtq906r79e90dEwPjxWue/cSP83//pAC+b0iHPyyoRhGdyzP52M/H7tt956feXuPXCW7nhAksE+VlkpE5IuW2bdmLZulUXqWrUCD74IHWhqgIrbfXxAw9A5cragDJvHvTrp/N/p+jbN2CLs5ucyyoRLBKRO9PvFJE7gL/8E1L+d/D4QW767iZql67Nuz3eDXQ4xkciIvR7cNMm7e5eogQMHgw1auicZnv2BDpCHzp0SKdsuPVWncEzJkb3d+wIL72ki7dER8NHH2lWNPlaVo3FFYHvgZOkfvG3AAoDVzvn/vF7hOnk9cZi5xx9J/Rl0rpJzLtjHi2q2LINBZVzOvbpzTe1qigsTGtJ7r8fWuTXv/alS7VOf948nbCtTBkdlv3ii1a/n8/luLHYObfHOdcWGAZs9WzDnHMXByIJ5Aejl4zm27Xf8vK/XrYkUMCJaHvBDz9otfidd+o65i1b6swH48fryml5Vsqgrrvv1sABSpfWKR6efFJbxvfu1RuxJFCgeb1UZY4+XKQ78A4QCox2zr2a7vhNQEoF41Hgbufc8sw+My+XCNbGrKX5h81pV6Md02+eTojY4OtgExsL48bBiBG6clr58pogBg2CmjUDHR1ajBk5Urt3zpqlmapYMV2t6/HHAx2d8SNfTDGRk4uGAiOBHkAj4AYRSb8i+xago3OuCfAi8KG/4skN7y16j2KFizHuqnGWBIJUyZJw333ae3LGDF0C99VX9YG6Z0+dSy0hIRcDOnlSv/DHjtXXItrCvWmTlgR+/lkHdVkSCGp+KxGIyMXAUOdcN8/rpwCcc6+c5fzSwCrnXNXMPjcvlwiSkpPYfHAz9crWC3QoJg/Ztg3GjNGp83fuTO10M3Cgn2pc9uxJHdT1889a1VO2rO4PDdWG4MhIP1zY5GUBKREAVYEdaV5He/adzR3A1IwOiMggEVksIotjUnov5CFrY9ay+8huQkNCLQmYM9Ssqaulbd0KkyZB8+ZaSqhTRxfSGj/+HLugJifrIIeUBol33oHbb4cFC7T1euJEvXioZ3S7JQGTjj9LBH2Bbs65gZ7X/YFWzrn7Mzj3UuA9oL1zbn9mn5vXSgTOOdp83Ia4hDhWDF6BWP9p44XoaC0ljB0LW7ZoV9R+/WDAAF1gK8t/RocOwfTp+uQ/dap275w+XRdZ2LYNDh7URZvt36Px8MXi9TkRDVRP87oasCv9SSLSBBgN9MgqCeRFIsK4q8ax99heSwLGa9Wq6diDZ57R2ZfHjoUvv9Tqozp1oH9/uPFGqJdSwHROp0YtWhRWrdKJj1K6d3bvrg0QrVrpuTVr5pGWaZNf+LNEUAhYD3QGdgKLgBudc6vTnFMDmAnc4pyb583n5qUSwcHjB4kMj7QEYHzi6FHtxTl2rC7CVdQd4656MxlQfgqNtk2h0NVXwrvvagJ46SVdZq1169QqH2MykVmJwN/dR3uik9OFAmOccy+LyGAA59woERkNXAts87wl8WyBpsgrieBE4gmaf9icS6Mu5d2eNnrY+Fbctf0p/MPXFEo6yVGK8QtdWH7+TVR9sA9XX61tv8ZkR6CqhnDOTQGmpNs3Ks3vA4GB/ozBX4bNGcbqmNW80eWNQIdi8rMTJ+C337SHz8qVOke/CBENa0LNe+Hyy9lZvj1Lvy3C55/Dxjt1Wot//Uun9bn6ap33zZhz4dcSgT/khRLBop2LaPNxG2698FbG9B4T0FhMPjVnDrz1Fvz6q47wLVIELr1U5/cpWTLDtzgHy5bpwjkTJuhQgNBQfVufPjrlT+XKuXwfJt8IVPfRAulE4gkG/DCAysUr81a3twIdjskPEhK00v+JJ3QuCtD1eJcv10ndfvwRDhzQ3j9nSQKgHYAuukinxt6wAZYs0Y/culVLCVWqwMUX66I669blyp2ZAsJKBNn00m8v8eysZ/npxp/oWa9nwOIwedyxY9oNaMoUre6JjdVZ6T75BG66Sfv+i/ike6dzsHq1DheYOBH+8kwP2aAB9O4NvXrp3EeF/FoRbPK6gDUW+0MgE8HGAxs5/73zubL+lXzd9+uAxGDyqIQEmD8fEhO1Aj9lNG/Fijp7Z8+e0LmzDhjws+3bdeDa999r80Nios4l1727JoXu3bXXqQkulgh8wDlH98+6syB6AWvvXUuVElVyPQaTx0RH6+Ls06bpxEKxsdChg377go4Ui4oK6KCuw4d1lomfftLCSUwMhIRoFVK3bro1b249UIOBJQIf+H7t91zz9TW82+Nd7mt1X65f3+QBJ09qvX7Llvq6e3cdzVutmv7eo4fOGZFJPX8gJSfDokWaFKZO1Sok57R00KWLJoWuXaFqprN9mfzKEoEPHE84zsdLP+buFncTGmKPT0Fjyxb9sp82LbWHz549Or/0X39pb5/GjfPlVA4xMdp8MX26bv94Vhhp0EBrtzp3hk6drBqpoLBEcI4SkxMpFGItbUEhZfa3iAj49FO45RZ9HRWlT/zdu+tjc3hmy3nnP87pMIaff9Z899tv+kchAs2aaWLo1Anat8+zBR6TBUsE52DhzoXc+O2NfH/d91xQ8YJcu67JJc7BmjWpT/2//QajRunsbzt26JwP3bvDeefly6f+nDp5Uic0/fVX3RYs0PbwkBDtwtqxo24dOmhDtMn7AjayuCBwzhEVGUXNSJvEq8BIStLW0dhYrdaJjtb9jRrBvffqhG4A1avr+r1BqHBhffpv314XL4uL005Rc+boNnKkjocT0T/Cdu1St1q1gipnFghWIjAFX0KCPtJOn651HzVqwDff6LGHH9YE0L27fvEbr8THw59/agFq7lxNErGxeqxyZU0IF1+s4xeaNStwNWn5kpUIcmBn7E6GLxjOcx2fo0QR//f9Nn7y2GO6NOPRo1oKaN1av6FSvP124GLLx8LDU6uHQAtZq1ZpUkjZUnJtWJgWstq00T/+Vq10qu0Qm9cgz7ASwVncOvFWvlz1JWvuWUOdMnX8fj1zjmJitDJ7xgz9Flq2TL+thg/X+Ra6dtVJeWx1rlzzzz9aaliwQLdFi7TTFUCpUtCihW4tW+rPGjWsSsmfrESQTQt3LmTc8nEMaTfEkkBeN2OGTrizbJm+jozULi4HDujkOw89FNj4glilSjrFRe/e+joxUdvlFy3SbfFibWdISNDjZctqQ3SzZqk/69a1kkNusBJBOs452o1px+aDm9lw/warFsorEhL02yOlG8uTT2q9/p9/6u9duuhgLhsmm6/Ex8OKFfpXu3SpTqS3alVqciheHJo00VU3U7bzz9f9JnusRJANX6z6gvnR8/n4yo8tCeQFBw9qX/45c3T+HhGtcD55Uo+3bg2zZgU2RpNj4eHaZpCyyiboX+3q1amJYfly+OwzeP99PS6ibQxNmmhSOP987blUr562R5jssxJBGsdOHqPByAZUKFaBRXcuIkSsTJprnNO6/Jkz9Yu9dm2dTzk5WbugNG2qQ10vvdSW5wpCzsG2bZoUVqzQn6tW6XTcycl6TliYjopu3BgaNtStUSNNEIULBzb+vMBKBF56Y94bRMdG8/k1n1sSyE2PPQaffw67d+vr6tVTV20PCdG+iSaoiejg7qio1DYH0Kqlv//WpLBqlY6OXrBAZwBPERqqJYiGDaF+/dM3W91NWSLw2HF4B6/PfZ1+jfvRoWaHQIdT8KQ88f/2my7SsnatlvtF9H9qp066de6spQHrPmK8EB6uhcWUMYAp4uL0n9uaNfpPbc0aTRhTpqS2P4DOo3TeefrcUa+eNk6n/F6qVO7eSyBZIvCYtG4SDsfrl70e6FAKhpTyekgIjB0LQ4bA3r26r3Jl/dI/elTn53/ttUBFaQqoiAjteXTRRafvT0zUKqZ161K39eu1NvLTT08/t1w5fSapU0e3lN9r19YOaQWpN5O1EaSxM3YnVUvaHLw5cuKEdv344w/d5s7VuY7btNEpLseNg0su0a1ePXviN3nO8eO6DvSGDbpt3AibN+u+7dtTn21A2xxq1NDpNKKiUn/WrKlb5cp5L1FYG0Emkl0yGw9s5Lyy51kSyI59+/R/RoUKWsXTtq0mA9DK2D59Uvv4XXaZbsbkYUWLpvZCSi8hQUsSmzfrtmWLrhW9ZYsuDxoTc/r5YWG6TEVKYqheXbdq1VJ/L1Uq7zwPBX0i+HLVl9zy/S3Mu2Meraq2yvoNwcg5rWidNy91W7cO/v1vePll7apx3306FWXbtjpXvzEFSFiYth/UrZvx8WPHNDFs364JI2Xbvl2HvezadXqJAvQ5qWrV07dq1fRnlSq6VayYO11ig75qKOZYDGOWjuHxdo9bT6EU+/bpQK3ERO2ikZysrWqHD2vFadu2unXvriN8jDGZSkzUTnE7duhktzt26LZzp27R0Xo8MfH094noc1WVKlrddMMN0L9/zmKwqqFMlC9WniHthwQ6jMD79FOt01+4UCtFQb/ke/fWys6vvtKKUKvfNybbChVKrRI6m+Rk7U+xc6cmhV27dEv5ffdunTnFL/H552PzvqMnj3L9N9fzfMfnaVm1ZaDDyR3Hj+tonCVLdJnFzZu13CqiSeD333UGsEGDdMRu8+ap7+3WLXBxGxMEQkJ0fqZKlXL/2kGbCIYvGM5PG37imUueCXQo/nHokH7pt26t6+q+9ho8/bTOFww6OrdlS63cLF4cPvlEzzPGBJ2gTAT74vbxxrw3uKrBVbSp1ibQ4fjG+vVavbN8uW7bt+v+xYv1yb51a23cbdZMt+rVT6/isSRgTNAKykTwyu+vcPTkUV669KVAh+K9pCTthrBmjc7IlTKm/j//0UXVd+zQ3+vX14bce+7ROv7zztP3p4zcNcaYdIIuEWw/vJ2Ri0Zyy4W30LhC40CHc6aDB3Uky/r1Oia+XTvtnbNxo3bTTFGtmnZ4TplN65JLdKRu0aKBidsYk28FXSIYNnsYDsfQjkMDE0BCgvYV27pVn/ArV9aG2BMntAPx/v2p54aEwLPPaiKoXRtGj06dXjH9SlthYTYHrzEmR/yaCESkO/AOEAqMds69mu64eI73BOKAAc65Jf6KZ23MWsYuH8uDrR+kZmRN31/gxAldny+lc/CuXVC6tM6nD1o3v3z56SNLrr5aE0GRInD77TqCpE6d1BmwUuruw8Lgjjt8H7MxJuj5LRGISCgwEugCRAOLRGSSc25NmtN6APU8W2vgfc9Pv3h65tMUCyvGU+2fyvgE5/TLPC5Ot2PH9HWTJnr8xx91hO3Bgzroau9eHWg1Zoweb99eG2fTatcuNRFccQX06qVjzlMmJqlRI/Xc123CO2NM7vNniaAVsNE5txlARL4EegNpE0FvYJzT4c0LRCRSRCo753b7Opj9cftp8/HPfLQ2nLKj22jja2Ki1rFv3qwn3XyzzoufVsWK+pQPWjXzww86bXL58rpFRKSe+8gjmkCqVEkdM16mTOrxYcN8fVvGGHPO/JkIqgI70ryO5syn/YzOqQqclghEZBAwCKBG2ifobCgbUZb7r32dsJmzIKyIDvULDT2922S/fnDBBfrlnrKl/SL/5BNNHBERGY+uveGGHMVmjDGB5M9EkNE8BOknNvLmHJxzHwIfgs41lNOAit51D9x1z9lP6N379OWP0itdOqeXNsaYPMufs6xFA2ln1qgG7MrBOcYYY/zIn4lgEVBPRGqJSGHgemBSunMmAbeIagMc9kf7gDHGmLPzW9WQcy5RRO4DpqPdR8c451aLyGDP8VHAFLTr6Ea0++ht/orHGGNMxvw6jsA5NwX9sk+7b1Sa3x1wrz9jMMYYkzlbicUYY4KcJQJjjAlylgiMMSbIWSIwxpggl+8WrxeRGGBbDt9eDtjnw3DyA7vn4GD3HBzO5Z5rOufKZ3Qg3yWCcyEii51zLQIdR26yew4Ods/BwV/3bFVDxhgT5CwRGGNMkAu2RPBhoAMIALvn4GD3HBz8cs9BbTZM1AAABfhJREFU1UZgjDHmTMFWIjDGGJOOJQJjjAlyBTIRiEh3EVknIhtF5MkMjouI/J/n+AoRaRaIOH3Ji3u+yXOvK0RknohcGIg4fSmre05zXksRSRKRPrkZnz94c88i0klElonIahGZk9sx+poX/7ZLichkEVnuued8PYuxiIwRkb0isuosx33//eWcK1AbOuX1JqA2UBhYDjRKd05PYCq6Qlob4M9Ax50L99wWKO35vUcw3HOa82ais+D2CXTcufD3HImuC17D87pCoOPOhXv+N/Ca5/fywAGgcKBjP4d7vgRoBqw6y3Gff38VxBJBK2Cjc26zc+4k8CWQfv3J3sA4pxYAkSJSObcD9aEs79k5N885d9DzcgG6Glx+5s3fM8D9wLfA3twMzk+8uecbge+cc9sBnHP5/b69uWcHlBARAYqjiSAxd8P0Hefcb+g9nI3Pv78KYiKoCuxI8zrasy+75+Qn2b2fO9Anivwsy3sWkarA1cAoCgZv/p7PA0qLyGwR+UtEbsm16PzDm3seATREl7ldCTzonEvOnfACwuffX35dmCZAJIN96fvIenNOfuL1/YjIpWgiaO/XiPzPm3seDgxxziXpw2K+5809FwKaA52BosB8EVngnFvv7+D8xJt77gYsA/4F1AFmiMjvzrlYfwcXID7//iqIiSAaqJ7mdTX0SSG75+QnXt2PiDQBRgM9nHP7cyk2f/HmnlsAX3qSQDmgp4gkOucm5k6IPuftv+19zrljwDER+Q24EMivicCbe74NeNVpBfpGEdkCNAAW5k6Iuc7n318FsWpoEVBPRGqJSGHgemBSunMmAbd4Wt/bAIedc7tzO1AfyvKeRaQG8B3QPx8/HaaV5T0752o556Kcc1HAN8A9+TgJgHf/tn8AOohIIRGJAFoDa3M5Tl/y5p63oyUgRKQiUB/YnKtR5i6ff38VuBKBcy5RRO4DpqM9DsY451aLyGDP8VFoD5KewEYgDn2iyLe8vOfngLLAe54n5ESXj2du9PKeCxRv7tk5t1ZEpgErgGRgtHMuw26I+YGXf88vAmNFZCVabTLEOZdvp6cWkS+ATkA5EYkGngfCwH/fXzbFhDHGBLmCWDVkjDEmGywRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+QsERhjTJCzRGAKFBEp65mCeZmI/CMiOz2/HxWR9/xwvbEisiWlX7vnda5Ndy0i13mmI/4xt65pCp4CN6DMBDfP1BlNAURkKHDUOfdfP1/2cefcN/68gIiEOueS0u93zn0lInuAx/x5fVOwWYnABAXPYi0/en4fKiL/E5GfRWSriFwjIq+LyEoRmSYiYZ7zmovIHM8sntOzMdXvJaKL/2xOKR14pgN4Q0RWea5zXfq4PK9HiMgAz+9bReQ5EfkD6CsiD4jIGs9iJF/68I/HBDkrEZhgVQe4FGgEzAeudc49ISLfA5eLyE/Au0Bv51yM54v7ZeB2Lz67Mjq7awN0XphvgGvQksqF6AR4izwTwmUl3jnXHkBEdgG1nHMnRCQyG/dqTKYsEZhgNdU5l+CZnyYUmObZvxKIQicuOx+d0hjPOd5O7DXRMx/+Gs8kaKCJ4QtP9c4e0SUkWwJZTZX8VZrfVwCfichEID9PnmfyGEsEJlidAHDOJYtIgkuddCsZ/X8hwGrn3MU5/WwPSfczvUROr6INT3f8WJrfL0eXMbwSeFZEGjvn8u1KXCbvsDYCYzK2DigvIhcDiEiYiDQ+h8/7DbhOREJFpDz6hb4Q2AY0EpEiIlIKz3TK6YlICFDdOTcLeAJdm7j4OcRjzClWIjAmA+7/27t/HIKCIADj34RL6RzAccQFHECiUCtUGo17qFQuoEOrGcVbhUIkT3hhv1+/u1PtZPZv5rVs9M7LBN2n+fFs37LLDTCg+Xw9gUlmHgEiYk2z7HMAdk/a94BViSWAWWaeW8YiPfAZaukNEbEEtp8+PvoihiEwzsxRVzHot7k0JL3nAkzvF8q+rZxmWgCnLsbXf7AikKTKWRFIUuVMBJJUOROBJFXORCBJlbsBoPZyIFSl/z4AAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}