{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "## Equivalent Circuit Bayesian Inference\n", "\n", "In this notebook, we introduce the Monte Carlo sampling methods for an example Thevenin model. This notebook includes importing experimental data for a Tesla 4680 NCA/Gr cell, and identifying the resistance elements of a two-RC Thevenin model. First, we import PyBOP and the required packages," ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Note: you may need to restart the kernel to use updated packages.\n", "Note: you may need to restart the kernel to use updated packages.\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "%pip install --upgrade pip ipywidgets -q\n", "%pip install pybop -q\n", "%pip install pandas -q\n", "\n", "import sys\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import pybamm\n", "\n", "import pybop\n", "\n", "pybop.plot.PlotlyManager().pio.renderers.default = \"notebook_connected\"\n", "parallel = True if sys.platform != \"win32\" else False" ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "Next, we will import and process the experimental data. The data for this notebook was obtained from the supplemental data provided in: \n", "\n", "[1] M. Ank et al., ‘Lithium-Ion Cells in Automotive Applications: Tesla 4680 Cylindrical Cell Teardown and Characterization’, doi: 10.1149/1945-7111/ad14d0.\n", "\n", "A portion of this data is retained in the PyBOP repository for use in parameterisation and optimisation examples. However, users' are pointed to the official location https://mediatum.ub.tum.de/1725661 for more information. First, we import the three-electrode pOCV dataset,\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [], "source": [ "ocv_df = pd.read_csv(\n", " \"../data/Tesla_4680/T-cell_pOCV_data.txt\",\n", " sep=\"\\t\",\n", " decimal=\",\",\n", ")" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "Next, let's plot the terminal voltage to visualise the protocol. " ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAGwCAYAAACdGa6FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVdElEQVR4nO3dd1zU9R8H8NcdB8e8YynbPXCAe4CZ5jYzNa30R4JmlpZlpWW0lyuzsiw0d6ZRmmaZIzXRTERAUBQ3IKgMReGOdcDd9/cHckUuxsH3xuv5eNyj7nvfu+8bv+q9/EyJIAgCiIiIiIyAVOwCiIiIiCoxmBAREZHRYDAhIiIio8FgQkREREaDwYSIiIiMBoMJERERGQ0GEyIiIjIaMrELqA6dToerV6/CyckJEolE7HKIiIioGgRBgFqthre3N6TS6rWFmEQwuXr1Kvz8/MQug4iIiGohIyMDvr6+1TrXJIKJk5MTgIofTKFQiFwNERERVYdKpYKfn5/+e7w6TCKYVHbfKBQKBhMiIiITU5NhGBz8SkREREaDwYSIiIiMBoMJERERGQ2TGGNCRESmS6vVoqysTOwyqB5YW1vDysrKoJ/JYEJERPVCEARkZWUhLy9P7FKoHjk7O8PT09Ng64wxmBARUb2oDCWNGzeGvb09F8g0M4IgoKioCDk5OQAALy8vg3wugwkRERmcVqvVhxI3Nzexy6F6YmdnBwDIyclB48aNDdKtw8GvRERkcJVjSuzt7UWuhOpb5T021DgiBhMiIqo37L4xf4a+xwwmREREZDTqFEwWLFgAiUSCl19++a7nrFixAn379oWLiwtcXFwwaNAgHD16tC6XJSIiIjNV62ASGxuL5cuXIzAw8J7nRUVFYcKECdi/fz+io6Ph5+eHIUOG4MqVK7W9NBEREZmpWgWTgoIChISEYMWKFXBxcbnnuRs2bMDzzz+Pzp07w9/fHytXroROp8O+fftqVTARmRdBEKAp14pdBpHepEmTIJFIbnsMGzasQa6flpZ2x+tLJBIcOXLEoNc6cOAA/Pz8MHLkyLv+fH/99RckEglOnDhh0GvfTa2mC7/wwgsYMWIEBg0ahI8//rhG7y0qKkJZWRlcXV3veo5Go4FGo9E/V6lUtSmTiEzAB78lY2NMOra/9ADaeFR/a3Si+jRs2DCsWbOmyjG5XN6gNezduxcdOnSocszQU6+3bduGkSNHYsiQIRg7diwuX74MX1/fKuesWbMG3bt3v28PiaHUuMUkMjISx44dw/z582t1wTlz5sDb2xuDBg266znz58+HUqnUP/z8/Gp1LSIybuez1Vh7OA2lWh0Onb8udjlUzwRBQFFpuSgPQRBqVKtcLoenp2eVh4uLC2bPno1HHnlEf94XX3wBiUSCXbt26Y+1atUKK1eu1D9fuXIl2rVrB1tbW/j7++Obb76pVg1ubm631WBtba1//bfffkOPHj1ga2sLd3d3jBkzRv+aRqPB7Nmz4ePjAwcHB/Tq1QtRUVG3XePXX3/Fo48+ikceeQSNGjXC2rVrq7xeUFCATZs2YcqUKdWq2RBq1GKSkZGBmTNnYs+ePbC1ta3xxRYsWIDIyEhERUXd8/3h4eF49dVX9c9VKhXDCZEZmrvjtP7/nWy53qO5Ky7Tov27u0W5dvKHQ2FvU/ffY/369cPKlSuh1WphZWWFAwcOwN3dHVFRURg2bBiuXLmCixcvon///gAqhjO8++67WLp0Kbp06YKEhARMnToVDg4OCAsLq3Udv//+O8aMGYO33noL3333HUpLS7Fjxw796zNmzEBycjIiIyPh7e2NrVu3YtiwYUhKSkLr1q0BAKdOnUJOTg4GDBgAmUyG0NBQrF27Fm+99ZZ+CvCmTZug1WoxYcKE2v+i1VCN7lJ8fDxycnLQtWtX/TGtVouDBw9i6dKl0Gg0d1317dNPP8WCBQuwd+/e+zYHyeXyBm8yI6KGdeDcNUSdvaZ/XsN/0BLVq+3bt8PR0bHKsTfffBPPP/881Go1EhIS0K1bNxw8eBCvvfYafvnlFwAVEz58fHzQqlUrAMB7772HxYsX47HHHgMANG/eHMnJyVi+fPl9g0lwcDCk0qodGwUFBQCAuXPnYvz48fjggw/0r3Xq1AkAkJ6ejjVr1iA9PR3e3t4AgNmzZ2PXrl1Ys2YN5s2bB6CiG2fo0KGwsbEBADz99NNYtGgRDhw4oA9Wa9aswdixY6FUKmv061cXNQomAwcORFJSUpVjkydPhr+/P+bMmXPXUPLJJ59g7ty52L17N7p37177aonILJSUaRG2uuqyATomE7NnZ22F5A+HinbtmnjooYcQERFR5ZirqyucnZ3RqVMnREVFwcbGBjY2Nnj22Wfx3nvvoaCgAAcOHEC/fv0AAIWFhbh48SKmTJmCqVOn6j+nvLxc/0U/fPhw/PXXXwCApk2b4tSpU/rzfvzxR7Rr1+6O9SUmJlb5zH9LSkqCVqtFmzZtqhzXaDRVxqhs27YNM2bM0D/39/dHcHAwVq9ejf79++PChQv466+/8OGHH97318uQahRMnJyc0LFjxyrHHBwc4Obmpj8eGhoKHx8f/RiUhQsX4t1338XGjRvRrFkzZGVlAQAcHR1vS6NEZBm+iboIAJBJJejezAVHUm5Ay2Bi9iQSiUG6UxqCg4ODvtXjv/r374+oqCjI5XL069cPrq6uaNeuHQ4dOoQDBw5g1qxZAP5p3VixYgV69epV5TMq/yG/cuVKFBcXA0CV8SMA4Ofnd9caKveouZOCggJYWVkhPj7+tgaDyu/dzMxMJCQkYMSIEVVenzJlCl588UV8/fXXWLNmDVq2bKkPWg3F4L9D0tPTqzQ9RUREoLS0FOPGjaty3nvvvYf333/f0JcnIiN3OlOFpX+eBwCM7eoLVUnF/ho65hIyEf369cPq1ashk8n0U2z79++PH374AefOndN3g3h4eMDb2xspKSkICQm542f5+PjUqobAwEDs27cPkydPvu21Ll26QKvVIicnB3379r3j+3/77TcEBwffNkP2iSeewMyZM7Fx40Z89913mD59eoNvK1DnYPLfUb7/fZ6WllbXSxCRmSjT6jDrp+PQCcCDbRphwdgAzNiYAADQMZmQEdFoNPoW/koymQzu7u548MEHoVarsX37dixYsABARTAZN24cvLy8qnShfPDBB3jppZegVCoxbNgwaDQaxMXF4ebNm1UmedxJbm7ubTU4OzvD1tYW7733HgYOHIiWLVti/PjxKC8vx44dOzBnzhy0adMGISEhCA0NxeLFi9GlSxdcu3YN+/btQ2BgIEaMGKGfjfNfjo6OePLJJxEeHg6VSoVJkybV8lew9rhXDhE1mE93n0Vypgou9tZY/HgnSCQSSKUV/xrjGBMyJrt27YKXl1eVxwMPPAAAcHFxQUBAABo1agR/f38AwIMPPgidTndbt8czzzyDlStXYs2aNQgICEC/fv2wdu1aNG/e/L41DBo06LYaKgfZ9u/fH5s2bcKvv/6Kzp07Y8CAAVW2e1mzZg1CQ0Mxa9YstG3bFqNHj0ZsbCyaNGmCwsJC7Nu3747BBKjozrl58yaGDh2qHzzbkCRCTSd3i0ClUkGpVCI/Px8KhULscoioFrYlXsHMyEQAwDchXfFwgBcAYGZkArYlXsXbI9rhmb4tRKyQDKmkpASpqalo3rx5rZaXoPqzZcsWvP3220hOTjbI593rXtfm+5stJkRU75Iu5+P1zRXLWU/v31IfSgBAeqv/WsuuHKIG4ejoiIULF4pdxl2ZxvBoIjJZOeoSPLc+DppyHfq3bYTZQ9pWed1TWfEvrEs3isQoj8jiDBkyROwS7oktJkRUb/KKSjFpdSyu5peghbsDlozvAitp1RH+7bwqmneTLueLUSIRGRm2mBBRvbheoEH3j/cCANwdbbB6Ug8o7axvO69384rpiiev5uN6gQbujlz12ZyYwDBGqiND32O2mBCRwWXcKMITy6MBVKy4ueGZ3mjm7nDHcxsrbNHRRwFBAHYmZTZkmVSPKhcLKypiF525q7zH/10grrbYYkJEBhWTkotn18cjv7gMNjIpNk7thbaeTvd8z2NdfHHySjI2xKTjqd5NG3xBJzI8KysrODs7IycnBwBgb2/P+2pmBEFAUVERcnJy4OzsfNdtaWqKwYSIDEKrE/DtwRR8+sdZaHUC/D2dsGZyD3gp7750dqWxXX2xcNcZnMlSI/7STXRv5nrf95Dx8/T0BAB9OCHz5OzsrL/XhsBgQkR19vuJTLy++TgKS7UAgEc7eWPB2IBq74uitLfG6M4++DEuA1/sPY/vn+l1/zeR0ZNIJPDy8kLjxo1RVlYmdjlUD6ytrQ3WUlKJwYSIaq2otBzf7L+IpfsvAAAcbKzwziPt8WQPvxo320/r3xI/xWfg0IXrOJp6Az2bs9XEXFhZWRn8y4vMF4MJEdVYSZkWz284hj/P/NNEr7Szxo6ZfeHjfP+umztp7u6AJ7v7ITI2A+9uO4ntLz4AmRXH5xNZGv6pJ6Jqyy8uw4qDKXjwk/36UOLrYocl4zsj8d3BtQ4llV4f5g+FrQxnstT4Yu95Q5RMRCaGLSZEdF/puUVYfyQNkbEZUJeU64+P7OSNReMCYWttmGZ6VwcbzBrSFu/9egqr/07F2G6+aH6XacZEZJ4YTIjojkrLddh5MhNrD6chIT1Pf7y5uwPCgpriyR5NYGdj+HEDT/VuijV/pyIttwgv/nAMm6cFGyz4EJHxYzAhoirOZ6vx1taTOJp2o8rxns1dMbVvCwzwb3zbsvKGZCWVIPLZIAxbchAnr6jw1taT+PTxQK6BQWQhGEyICGVaHf44lY31R9JwJOWfQOLuKMfQDh4IC26GNh73XiTNkDyVtvhyfBdMWnMUPx+7jMYKOeYM82+w6xOReBhMiCzYzcJSfBd9Cd/HXMI1tQYAIJUAPZq5Iiy4GYa09xBtZsyDbRph1pC2WLT7LCKiLsJRLsMLD7USpRYiajgMJkQWKLdAg/AtSfgjOVt/zN1Rjv/19MP4nk3gXcfZNYbywkOtoNMJWLznHBbtPgtVSRleH+pfr11JRCQuBhMiC6LVCZgZmYDtJ/7ZLK9lIwc8378VRnbyho3M+FYQeHFga0ilEizafRbLD6Qg6XI+vniyMxorbMUujYjqgUQwgT2pVSoVlEol8vPzoVAoxC6HyCRdzSvGSz8kIO7STf2x90e2R1hwM5MYWPrr8at44+cTKCrVwt3RBu8/2gEjArxMonYiS1Wb728GEyILsCMpE+FbkpBfXLFfybR+LTF7SBuTW1n1Qk4Bpn8fj/M5BQCAgf6N8e7I9mjqxrVOiIwRgwkRVaHVCfhk1xksP5gCAOjkq8SS8V3QzIQXLSst1+Hr/Rfw9f4LKNf989fXrzP6INDXWbzCiOg2DCZEpJdfXIaXfkjAgXPXAABT+zbH68P8YW1irSR3cyFHjXd+OYXolFz9sRGBXnhjmD/8XO1FrIyIKjGYEBEA4EpeMSauikHKtULYWkuxaFwnjOzkLXZZBicIAr7cdwGf7z1X5XiPZi54JNAbE3s3hZQzeIhEw2BCREi+qsLDX/4FAPBS2mJFaHd09FGKXFX925Z4BcsOpOBMlgqVf6u52FvDz9UeHzzaAV2auIhbIJEFYjAhsnDH0m8ibNVRqDXl8FLa4ufpwUazJklDSc8twsTVMbiUW1TleNcmzhjdxQdjuvjAydZapOqILAuDCZEFi0nJxdNrY1FYqkVHHwVWh/Ww6LU+Ssq0+C46DfN2nLnttREBXgjp3QRBLdw43ZioHjGYEFmolX+l4OPfTwMAglu6YWVYd9jbcP3ESjmqEmxLvIq1h9NwJa9Yf7yDtwJP9W6Kx7r6QC7jDsZEhsZgQmSB/h1KujRxxg9Te8PWml+yd6LTCUjIyMOmuAxsTbgCTblO/1qgrxI/Tw82m1lLRMaAwYTIwvx5JhtPr43TPz/38XCjXFbeGN0sLMXPxy5jxV8pyFZVbGDYurEjPni0A4JbuYtcHZF5YDAhsiCHzl/H0+tiUXrrX/1nPhrGlpJa0JRrMf7bI0hIz9MfG97RE/PGBMDFwUa8wojMQG2+v/lPKyITdCQlF898VxFK+rZ2x9mPGUpqSy6zwtbn+yDx3cEY38MPALDzZBbGLTuMk1fyRa6OyPIwmBCZmPhLN/D02liUlOnQv20jrAzrzoGbBuBsb4MFYwPx3dM9AQAXrxVizDd/46fYDJErI7IsDCZEJuTklXxMWhOLolIt+rRyw7KnujGUGNiDbRrh2DuDMaS9B8q0Al7/+QTm7zgNnc7oe72JzAKDCZGJSLqcj0e+OgR1STl6NHPBitDu7L6pJ64ONlj2VDe8PKg1AGD5wRS8+lOifjwPEdUfBhMiE3AuW43Q1TEAKqa1rprUg+uU1DOpVIKXB7XB4sc7QSaV4JfEqwhdHYNCTbnYpRGZNQYTIiOXer0QIStjcLOoDK0aO2L1pB5QcEn1BjO2my9WT+oBBxsrHEm5ga4f7UFeUanYZRGZLQYTIiN2+WYRQlYcwTW1Bv6eTtg8LQjujnKxy7I4D7ZphHW3BsVqynUY8eUhhhOiesJgQmSkclQlCFkZg6v5JWjRyAHrp/SCsz3X1RBL92au+OOVB2FjJcWVvGKErT4KVUmZ2GURmR0GEyIjlKOuCCWXcovg52qHDc/0QiMntpSIrY2HE36aFgQnWxmOX87H1HVxKCnTil0WkVmpUzBZsGABJBIJXn755Xuet2nTJvj7+8PW1hYBAQHYsWNHXS5LZNauqTUYv/wIzucUwFNhi43P9IaX0k7ssuiWzn4V+xE5ymWISb2Bl35IQLmWs3WIDKXWwSQ2NhbLly9HYGDgPc87fPgwJkyYgClTpiAhIQGjR4/G6NGjcfLkydpemshsXVNr8OS30Ui5XghvpS0in+0NP1d7scui/+joo8SK0O6wkUnxR3I23tp6EiawuweRSahVMCkoKEBISAhWrFgBFxeXe567ZMkSDBs2DK+99hratWuHjz76CF27dsXSpUtrVTCRubpZWIqw1UeRcq0Q7o42WP9MLzRzdxC7LLqLoJZu+GpCF0glwI9xGVi466zYJRGZhVoFkxdeeAEjRozAoEGD7ntudHT0becNHToU0dHRd32PRqOBSqWq8iAyZ/lFZfjfyhgkZ6rg6mCDzdOC0bKRo9hl0X0M7eCJ+Y8FAACWHbiIFQdTRK6IyPTVOJhERkbi2LFjmD9/frXOz8rKgoeHR5VjHh4eyMrKuut75s+fD6VSqX/4+fnVtEwik1GgKcektUdxOlMFd0c5fny2N1tKTMiTPZpgzjB/AMDcHaexOf6yyBURmbYaBZOMjAzMnDkTGzZsgK2tbX3VhPDwcOTn5+sfGRncRIvMk04n4OXIBCSk50FpZ421k3ugtYeT2GVRDU3r1wJT+zYHAMz5+QT2JmeLXBGR6apRMImPj0dOTg66du0KmUwGmUyGAwcO4Msvv4RMJoNWe/u0OU9PT2RnV/1Dmp2dDU9Pz7teRy6XQ6FQVHkQmaPFe85i7+kc2MikWPd0T3T0UYpdEtWCRCLBmw+3w9iuvtDqBLyw8Rhi026IXRaRSapRMBk4cCCSkpKQmJiof3Tv3h0hISFITEyEldXtG4oFBQVh3759VY7t2bMHQUFBdaucyMT9dvwqvt5/EQCwcGwAOvs5i1sQ1YlEIsHCsQEY6N8YmnIdnl4Ti9TrhWKXRWRyahRMnJyc0LFjxyoPBwcHuLm5oWPHjgCA0NBQhIeH698zc+ZM7Nq1C4sXL8aZM2fw/vvvIy4uDjNmzDDsT0JkQk5dzcdrm48DAJ59sAXGdPEVuSIyBJmVFF9O6IL2XgqoNeWYsi4W+cVcHZaoJgy+8mt6ejoyMzP1z4ODg7Fx40Z8++236NSpEzZv3oxffvlFH2SILI2mXItXfkxESZkO/do00g+cJPPgIJdh7dM94KW0Rcq1QrzIBdiIakQimMCqQCqVCkqlEvn5+RxvQibvsz3n8OW+83BzsMGeV/vB1YH735ijk1fy8fiyaBSXaTG5TzO8N7KD2CURNbjafH9zrxyiBpR8VYVv9l8AAHw4qiNDiRnr6KPEZ090AgCs+TsNG2IuiVwRkWlgMCFqIGVaHV7bfBzlOgHDOnji4YC7z0wj8zA8wAuvDm4DAHhr60kcTeVMHaL7YTAhaiDrDqfh1FUVlHbW+HB0B0gkErFLogbw4oBWaHJrv6OXfkjA9QKNyBURGTcGE6IGcCm3EIv/OAcAeGO4Pxo71d8ChWRcJBIJdszsiyau9shSleDlyETodEY/tI9INAwmRPVMqxMw66fjKC7ToncLVzzZnVssWBpHuQzfhnaDXCbFoQvXsepQqtglERktBhOierbswEXEXboJR7kMi8Z1glTKLhxL5O+pwLsj2wMAPtl9Biev5ItcEZFxYjAhqkcxKbn4fE9FF877j3aA362xBmSZ/tezCYa090CZVsDsTcehKb99Gw8iS8dgQlRPLt8swgsbj6FcJ2BkJ2+M7eojdkkkMolEgvmPBcDNwQZnstT4ct95sUsiMjoMJkT1IK+oFFPWxuF6QSn8PZ3wydhAzsIhAICboxwfjqpY+XrZgRSczlSJXBGRcWEwITKwQk05pqyLw9lsNRo5ybF6Ug/Y2dy+wSVZrhGBXhjS3gNanYC3tiZxlg7RvzCYEBlQUWk5Jq+NRfylm1DYyrB+Sk94O9uJXRYZofcf7QAHGyscS8/DD7HpYpdDZDQYTIgMpKi0HE+vjcXR1Btwksuw7ume8Pfk3k50Z97Odpg1pC0A4LM/zkFdwl2IiQAGEyKDyC8qw1MrY3Ak5QYc5TKsm9ITXZq4iF0WGbmnejdFC3cH5BaWIiLqotjlEBkFBhOiOsq4UYSxyw7jWHoelHbW+G5KT3RlKKFqsJFJMWe4PwBg5aFUXMotFLkiIvExmBDVwYnLeRjzzWFcyCmAp8IWG6f2YiihGhnS3gPBLd1QWq7DJ7vPil0OkegYTIhqad/pbDy5/AiuF2jg7+mErS8Eo4O3UuyyyMRIJBK8PaI9JBLg9xOZiL90U+ySiETFYEJUC+uj0zD1uzgUl2nRt7U7Nk0LgpeSs2+odtp7KzCuqy8A4J1fTnL6MFk0BhOiGijX6vDR9mS8s+0UdALwZHc/rJ7UA0621mKXRibujeH+cJLLkJypwvojl8Quh0g0DCZE1aQuKcPU7+L0O8POGtwGC8YGwNqKf4yo7twc5XhlcBsAwKe7zyIrv0TkiojEwb9RiarhfLYaI786hP1nr8HWWopvQrrixYGtucw8GVRYcDN08nOGWlOOWZsSUa7ViV0SUYNjMCG6jw0xlzD484NIyy2Cj7MdIp8NwsMBXmKXRWbISirBp+MCYW9jhb8v5GLBzjNil0TU4BhMiO5CECr2MXlr60kAQK/mrvh1Rh909nMWtzAya609nLD48U4AKtY22XLsssgVETUsBhOiO9CUa/H65hPYEFOxh0mXJs74bkpPuDnKRa6MLMHwAC+8OKAVAOCNLUk4cTlP3IKIGhCDCdF/XFNrELIiBpviL0MiqZh5s2V6MOQy7hBMDeeVQW0w0L8xSst1eG59PDJuFIldElGDYDAh+pfz2Wo8uvQQ4i7dhL2NFdZO7omF4wI5yJUanFQqwefjO6NlIwdk5pfgfyuP4PJNhhMyfwwmRLesOJiCwZ8fRGZ+CXyc7bBpWhD6tWkkdllkwRS21tjwTG80cbVHxo1iPL4sGhdy1GKXRVSvGEyIAKw6lIq5O04DADr6KLBtRh8uL09GwVNpi5+eC0Krxo7IzC/B2IhoHE29IXZZRPWGwYQsmlYn4MPfkvHR9mQAQI9mLvh5ejDcOciVjIin0hY/Ptsb7b0UyC8uwxPLo7EhhqvDknliMCGLVVyqxXPr47H674qVXEcEeOGn54I4yJWMkpujHD9PD8YA/8YAgLe2nsTne85BELivDpkXBhOySOm5RRj02QHsPZ0NG5kUX03ogq9DunKQKxk1OxsrrAjtjtGdvQEAS/adx+hvDqNQUy5yZUSGw2BCFmff6Ww8uGg/ruQVw8XeGhue6YWRnbzFLouoWqykEnwxvgsWPBYAADiekYexEYdxJa9Y5MqIDIPBhCzKH6eyMGVdHADAwcYKW57vgx7NXEWuiqjmxvdsgiXjOwMAzmSp0WfBnzhw7pq4RREZAIMJWYxfEq5g2vfx+ud/zRmA5u4OIlZEVDejOvvg8BsD0MzNHgAQtvooxkYc5rgTMmkMJmQR1h+5hFd/SoROAFq4O+D4e0Pg6mAjdllEdebtbIedMx/EsA6eAID4Szfx+LJoFHDcCZkoBhMya4Ig4NPdZ/HOLyehE4CQXk2w99V+UNpZi10akcHY2Vgh4qmueL5/SwBA3KWbGBdxGGnXC0WujKjmGEzIbJVrdQjfkoSl+y8AqNh75OPRHSGVcuYNmR+JRILXh/ljy/MV6/CcyVJj5NJD+PNMttilEdUIgwmZpTKtDi//mIjI2AxIJcC8MQGYOag1pwOT2evaxAW/v/QAujZxhrqkHFPWxeHzPeeg03HcCZkGBhMyOyVlWjy/4Ri2n8iEtZUE34R0xf96NRG7LKIG46GwReSzQZjYuykEoWK9k2fXx0FVUiZ2aUT3xWBCZqW4VIup38VhT3LFwmnLJ3bDsI5eYpdF1OBsZFJ8NLojFj/eCTYyKfaezsGopX8j+apK7NKI7onBhMyGuqQMYauP4q/z12FvY4W1k3pggL+H2GURiWpsN19snhYEL6UtUq8X4uEv/+K4EzJqDCZkFm4WliJkZQyOpt2Ak60M66f0QnArd7HLIjIKgb7O+HXGA2hxa92ep9fGYd6O01zvhIxSjYJJREQEAgMDoVAooFAoEBQUhJ07d97zPV988QXatm0LOzs7+Pn54ZVXXkFJSUmdiib6txx1CcZ/ewQnLufD1cEGP0ztjW5NXcQui8ioNHKSY9fLD2JCz4rxVt8eTEHz8B24WVgqcmVEVdUomPj6+mLBggWIj49HXFwcBgwYgFGjRuHUqVN3PH/jxo1444038N577+H06dNYtWoVfvzxR7z55psGKZ7oal4xnlx+BGez1WjsJMePz/ZGRx+l2GURGSUbmRTzxnTESwNa6Y89tSoGWfn8xyIZD4lQx7Y8V1dXLFq0CFOmTLnttRkzZuD06dPYt2+f/tisWbMQExODQ4cO3fUzNRoNNBqN/rlKpYKfnx/y8/OhUCjqUi6ZkUu5hfjfihhcySuGj7MdNk7thaZuXGKeqDp+isvA65tPAAA8FHKsn9ILbTycRK6KzI1KpYJSqazR93etx5hotVpERkaisLAQQUFBdzwnODgY8fHxOHr0KAAgJSUFO3bswMMPP3zPz54/fz6USqX+4efnV9syyUydz1bj8WXRuJJXjObuDtg0LYihhKgGnujuh4OvPYTWjR2RrdJg9Nd/c1AsGYUat5gkJSUhKCgIJSUlcHR0xMaNG+8ZNL788kvMnj0bgiCgvLwc06ZNQ0RExD2vwRYTupeTV/IRuvoobhSWoq2HE9Y/0xONnWzFLovIJN0sLMW07+MRk3oDVlIJ5o3piCd7cN0fMowGaTFp27YtEhMTERMTg+nTpyMsLAzJycl3PDcqKgrz5s3DN998g2PHjmHLli34/fff8dFHH93zGnK5XD/AtvJBBFRsUDZhxRHcKCxFoK8Skc/2ZighqgMXBxt8/0wvPNbFB1qdgDk/J+GplTGcsUOiqfMYk0GDBqFly5ZYvnz5ba/17dsXvXv3xqJFi/THvv/+ezz77LMoKCiAVFq9XFSbxEXm5/CF63jmuzgUlWrRs5krVk3qDidbbsZHZAiCIOCzPefw1Z8Ve0tN6OmHuaMDuLcU1UmDjjGppNPpqnS7/FtRUdFt4cPKygoAmMapRv48k41Ja2NRVKpF39buWPd0T4YSIgOSSCSYNaQtHg7wBAD8cDQDk9fGoqRMK3JlZGlqFEzCw8Nx8OBBpKWlISkpCeHh4YiKikJISAgAIDQ0FOHh4frzR44ciYiICERGRiI1NRV79uzBO++8g5EjR+oDCtH9/H4iE89+F4/Sch0Gt/fAyrDusLPh7x+i+vBNSDd89kQnAMCBc9cQtvoo99ihBiWryck5OTkIDQ1FZmYmlEolAgMDsXv3bgwePBgAkJ6eXqWF5O2334ZEIsHbb7+NK1euoFGjRhg5ciTmzp1r2J+CzNbm+Mt4ffNx6ARgVGdvfPp4J1hbccFiovr0WFdfeCpt8ex3FYNiH4+Ixrqne8JTyfFcVP/qPMakIXCMiWVaH52Gd7ZVLN43vocf5o4JgBX7u4kazKmr+Zi0JhbX1Br4ONthwzO90Myd0/Kp+kQZY0JUH5YfuKgPJZP7NMP8xxhKiBpaB28ltkwPRhNXe1zJK8bjy6Nx8kq+2GWRmWMwIaNSOTNg/s4zAIAXB7TCu4+0h0TCUEIkBj9Xe/z4XG+0buyIa2oNQlbGIP7SDbHLIjPGYEJGQxAEzP39NL7cdx4A8Pqwtpg1pC1DCZHIvJR22DwtGJ38nJFfXIbQVUfx94XrYpdFZorBhIxCuVaHN7eexMpDqQCADx7tgOf7t7rPu4iooSjtrfHD1F54oJU7Cku1eHptLJewp3rBYEKiKynT4pnv4vDD0XRIJcAn4wIRFtxM7LKI6D/sbWRYGdYdA/0bQ1Ouw9Nr4/DJrjNil0VmhsGERJVfXIaw1UcRdfYarK0k+Pp/XfFEd27aSGSsbK2tsGxiN/1CbN9EXcQvCVdErorMCYMJiSbjRhHGRRxGTOoNOMplWDe5J4YHeIldFhHdh7WVFEvGd9E/f+WnRHx/5JKIFZE5YTAhURzPyMOYbw7jfE4BPBW2+PG53ghu5S52WURUTdZWUqTMexihQU0hCMDbv5zE/J2nxS6LzECNVn4lMoTdp7IwMzIBJWU6tPNSYPWk7vBS2oldFhHVkFQqwQePdkBxqRab4i9j+YEUuNrb4Ll+LcUujUwYgwk1GEEQsOpQKubuOA1BAPq3bYSl/+sKRzl/GxKZKolEgk/GBcJBLsPaw2mYv/MMCjXleGVwG071p1rhNwI1iHKtDh9uT8Z30RX90CG9muCDRztAxn1viEyeRCLB+492gIu9DT7few5f/nkBJeU6hA/3ZzihGmMwoXqXV1SK8d8ewZksNSQS4M3h7fBM3+b8C4vIzMwc1BoKOxk++C0Z3x5MgY2VFLOHthW7LDIxDCZUr5Iu52P6hnhcvlkMuUyKL57szJk3RGZscp/mkAB4/7dkLN1/AQ5yGab355gTqj62o1O9EAQBG2PSMTbiMC7fLIaTXIblE7sxlBBZgEl9muPNh/0BAAt3ncGczSeg0xn9RvZkJNhiQgZXVFqOt7eexJZbiy4NaueBxU90gtLOWuTKiKihTO3bAhk3irH+yCX8GJeBMp0Oi8Z14i7hdF8MJmRQF68V4Pnvj+FsthpSCfD6MH8892ALjichsjASiQQfje4IhZ0MX++/iC3HruDQ+ev4a85DkMusxC6PjBi7csggBEHAjI3HMHDxAZzNVsPdUY6NU3tjWr+WDCVEFuy1of74JqQrACBHrcGILw+hpEwrclVkzBhMqM7Sc4sw7Iu/sP1EJgCgV3NX7HjpAfRu4SZyZURkDB4O8ML7I9sDAC7kFOCZdXEo1JSLXBUZKwYTqjVBEPDKj4l4cNF+nM1WAwBCg5rih6m90VhhK3J1RGRMJvVpjk/GBgIADl24jomrYpBfXCZyVWSMGEyoVrJVJQhbE4uttwa4eips8duMB/DhqI6QcnAbEd3BEz38sOX5YCjtrHEsPQ9hq48ynNBtGEyoRnS6imXlB3wahYPnrgEAwoKa4q85DyHAVylydURk7Lo2ccGGZ3pBaWeNxIw8PLEsGtmqErHLIiMiEQTB6CeXq1QqKJVK5OfnQ6FQiF2OxTqXrcbU7+JwKbcIANDJzxmLHw9Eq8ZOIldGRKbmdKYKYauPIketQRNXe6yf0hNN3RzELosMrDbf32wxofsqKdNi0e4zeHjJX7iUWwSpBHhvZHv8PC2IoYSIaqWdlwI/Tw9GUzd7pN8owphvDuNMlkrsssgIMJjQPf11/hqGfH4QX++/iHKdgIfaNsK2Fx7A5D7NuQEfEdWJn6s9Nj0XhA7eCtwoLMWopX9jZ1Km2GWRyNiVQ3d0Ta3Bx78nY1viVQAVg1s/GNUBQzt4ilwZEZmb/KIyPL0uFvGXbgIAvprQBSM7eYtcFRkCu3KoznQ6AT8cTcfAxVHYlngVUgkwuU8z7J3Vj6GEiOqF0t4aG57pBQ+FHADw4g8J+GTXGZGrIrFwSXrSi027gceXReufd/RRYN6YAAT6OotXFBFZBFtrK/w9ZwCmfX8Me09n45uoiyjT6vDG8HbcX8fCsMWEkFdUijmbT+hDib2NFd55pD1+eb4PQwkRNRiZlRQrQrvpu3FW/JWK59bHoYCrxFoUjjGxYDqdgF+PX8VH25ORW1gKAAj0VWLZU93g7WwncnVEZMl+PX4Vr206Dk25Dv6eTlgZ1h2+LvZil0U1VJvvbwYTC3XxWgHCf07C0bQbAIDWjR0x77EA9GjmKnJlREQVEtJvYup38bheoIG7ow2WT+yObk1dxC6LaoCDX6lafj1+FSO/OoSjaTdgb2OF2UPa4LcXH2AoISKj0qWJC36d0QftvBS4XlCKCSuO4Jdb22CQ+WKLiQURBAHTvo/H7lPZAICgFm749IlO8GG3DREZsUJNOV7+MRF7kiv+7npxQCu8MqgN9+UyAWwxoXtasu+8PpQ83ac5vn+mF0MJERk9B7kMy5/qhuf6tQAAfPXnBcz44RiKS7UiV0b1gcHEQmyMSccXe88DAKY80BzvjmzPKXhEZDKkUgnCh7fDonGBsLaSYEdSFp78NhpX84rFLo0MjMHEAvxxKgtv/5IEAHhpQCu880h7kSsiIqqdx7v7YcMzveFib40Tl/MRvOBPJF/lHjvmhMHEzB1NvYGZkYnQCcDj3XzxyuA2YpdERFQnPZu7YtsLD8DexgoAMDbiMLafuCpyVWQoDCZmLP7STUxecxTFZVo81LYR5j8WAImE3TdEZPqauNnjr9cfQlALNxSXaTFjYwLe23YSpeU6sUujOmIwMVPRF3MRtvooCku16N3CFRFPdeNuwERkVtwc5Vg/pSeee7BiUOy66Eto8/ZOZNwoErkyqgt+U5mhHUmZCFtzFAWacgS1cMOqsB6wtbYSuywiIoOTWUkR/nA7LP1fF/2xMd/8jeiLuSJWRXXBYGJmJq05iuc3HENpuQ792zbC6kk94CDnXo1EZN4eCfTG9hcfgKNcpl+M7ct956HTGf1SXfQfNQomERERCAwMhEKhgEKhQFBQEHbu3HnP9+Tl5eGFF16Al5cX5HI52rRpgx07dtSpaLqdIAj4aHsyos5eA1Ax0HVVWA/Y2bClhIgsQ0cfJWLeHIjB7T0AAJ/tOYfHIg4jK79E5MqoJmoUTHx9fbFgwQLEx8cjLi4OAwYMwKhRo3Dq1Kk7nl9aWorBgwcjLS0NmzdvxtmzZ7FixQr4+PgYpHiqUBFKTmPVoVQAFaHkk3GBXKeEiCyOg1yGFaHdMW9MAGytpUjMyEPv+ftw8kq+2KVRNdV5SXpXV1csWrQIU6ZMue21ZcuWYdGiRThz5gysra1rfQ0uSX9v83acxrcHUwAAc8d0REivpiJXREQkvpNX8vHE8mgUlWohl0kx46FWeHFga7HLsigNuiS9VqtFZGQkCgsLERQUdMdzfv31VwQFBeGFF16Ah4cHOnbsiHnz5kGrvfcywhqNBiqVqsqD7mzd4TR9KPloNEMJEVGljj5KRL8xEA+1bQRNuQ6L95zDW1uTUKbllGJjVuNgkpSUBEdHR8jlckybNg1bt25F+/Z3Xkk0JSUFmzdvhlarxY4dO/DOO+9g8eLF+Pjjj+95jfnz50OpVOoffn5+NS3TIvxxKgsf/FbRjRYW1BQTezOUEBH9m9LeGqvCeuj32dkQk46nVsbgmlojcmV0NzXuyiktLUV6ejry8/OxefNmrFy5EgcOHLhjOGnTpg1KSkqQmpoKK6uKQZifffYZFi1ahMzMzLteQ6PRQKP55zeNSqWCn58fu3L+5WjqDUxcFQNNuQ5PdvfDgrFcPI2I6F52n8rCrJ+Oo0BTDgD44NEOCAtuJm5RZq42XTk1nkdqY2ODVq1aAQC6deuG2NhYLFmyBMuXL7/tXC8vL1hbW+tDCQC0a9cOWVlZKC0thY2NzR2vIZfLIZfLa1qaxYi/dANT1sVCU67DoHYemDumI0MJEdF9DO3giVYzHPHc+nhcyCnAe7+egrO9NUZ15oQMY1LndUx0Ol2V1o1/69OnDy5cuACd7p/+vHPnzsHLy+uuoYTu7UyWCpPXxEJdUo7uTV2w9H9duKIrEVE1tWzkiF9e6KN/PjMyEZ/9cZbrnRiRGn2jhYeH4+DBg0hLS0NSUhLCw8MRFRWFkJAQAEBoaCjCw8P150+fPh03btzAzJkzce7cOfz++++YN28eXnjhBcP+FBbi4rUCPLUyBqqScnT2c8bap3tyRVciohpylMuQMu9h/biTL/+8gBd/SEBx6b0nZlDDqFFXTk5ODkJDQ5GZmQmlUonAwEDs3r0bgwcPBgCkp6dDKv0n6/j5+WH37t145ZVXEBgYCB8fH8ycORNz5swx7E9hAVKuFWD8t0dwvaAU7b0UWDe5Jxy5oisRUa1IpRKED2+HVo0c8ebWJPyelImMm0VYEdodHgpbscuzaHVex6QhWPo6JlfzijEu4jCu5pegdWNHRD7bG26OHINDRGQIMSm5mPZ9PG4WlcFTYYuVYd3R0UcpdllmoUHXMaGGkVugwVOrYnA1vwQtGjngB4YSIiKD6tXCDb+80AetGjsiS1WCx5dFY9fJLLHLslgMJkasQFOOKevikHKtEJ4KW6yf0gvuDCVERAbX1M0BW54PRt/W7igu02La9/FYspebAIqBwcRIlZbr8Nz6OCRm5EFpZ43vpvSEj7Od2GUREZktha011kzqgbCgisUqP997Di3e3IH84jKRK7MsDCZGqFyrw8s/JuDvC7lwsLHCuqd7oo2Hk9hlERGZPZmVFB+M6ohPxgbqj43/9ghy1NyhuKEwmBiZcq0Oszcdx46kLFhbSbBsYjd09nMWuywiIovyRA8/fDKuIpyczlThsW8OI+VagchVWQYGEyNSrtXhlZ+O45fEq7CSSvDVhK7o27qR2GUREVmkJ7r7IWp2fzR1s8flm8UYG3EYCek3xS7L7DGYGIkyrQ4vRSbgt+NXIZNK8PX/umJYR0+xyyIismjN3B3w8/RgBPoqcbOoDI8vi0ZMSq7YZZk1BhMjUFquw4yNx/TdNxFPdWMoISIyEu6OcvwwtTd6NXdFuU7AM+vicCZLJXZZZovBRGSaci2e3xCP3aeyYSOT4tuJ3TG4vYfYZRER0b84yGVYM7kHOvs5Q60px9NrYpGj4oDY+sBgIqLcAg36LNiPvadzIJdJsSK0Ox7ybyx2WUREdAf2NjKsndwDLRo54Gp+Caauj4emnPvrGBqDiUhSrxei28d7cb2gYmfm1ZN6oF8bDnQlIjJmzvY2WDOpBxS2MhzPyMP8HWfELsnsMJiIYP/ZHDz2zd/65xue6YU+rdxFrIiIiKqrqZsDvhjfGQCw9nAadp/i8vWGxGDSgMq1Ony6+ywmr4nFzaIyBPoqcfStgQwlREQmZoC/B6b2bQ4AmL3pOC7lFopckflgMGkgV/KK8cTyaCzdfwEAMKFnE2yaFoTGTtxem4jIFL0+zB/dmrpAXVKOad8fQ3Epx5sYAoNJPRMEAdsSr2Dg4igcS8+DrbUUnz/ZCfMfC4BcZiV2eUREVEvWVlJ8/b+ucHOwwelMFV6KTICWm/7VGYNJPcpWlWDqd/GYGZmIkjIdOvs5Y9fMBzGmi6/YpRERkQF4Km2xbGI32Mik2JOcjYW7OBi2rhhM6oEgCIg8mo5Bnx3A3tPZsLaS4KWBrbFpWhCauTuIXR4RERlQj2auWHRrX51vD6bg8z3nRK7ItMnELsDcpOcW4Y0tJ3D4YsWSxZ38nPHJ2EC09eTuwERE5mpUZx9cvFaIL/edx5J95zGsoyfaeSnELsskscXEQLQ6AasOpWLoFwdx+GIubK2leOvhdtgyPZihhIjIArw8sDVa3GoVf+XHRJSUcTBsbTCYGMCFnAKMW3YYH21PRnGZFr1buGLXzAcx9cEWsJJKxC6PiIgagFQqwU/TguDmYIMzWWp8xi6dWmEwqQOtTsATy6Mx6LMDSEjPg6NchrljOmLjM705loSIyAK5O8qxcOw/401OXM4TtyATxGBSS7kFGkxcFYOjqTcAAH1aueGPVx5ESK+mkLKVhIjIYg1q74Hht3aIf2fbKeg4hbhGGExqQVVShvHfHtEPcJ3Yuym+n9IL3s52IldGRETG4INRHeAor9hP57voNLHLMSkMJrWwcOcZnM8pgIdCjj2vPIiPRneERMJWEiIiqtDYyRavDW0LAJi/8wwu5KhFrsh0MJjUUGzaDWyISQcALBnfBa09OOOGiIhuFxrUFMEt3aAp1+HZ7+KRoyoRuySTwGBSQ8uiLgIAnujui94t3ESuhoiIjJVEIsGS8V3grbRFyvVCjPjqEGJScsUuy+gxmNTANbUG+8/mAACe69dS5GqIiMjYNXKSY8PU3mjd2BHX1BqMX3EEH21PhrqkTOzSjBaDSQ3sPpUFnQB08lWiZSNHscshIiIT0NzdAb+80AfDO3pCEIBVh1IR8P4f+OFoOgSBM3b+i8GkBnaezAQADA/wErkSIiIyJQ5yGSKe6oY1k3voj4VvSULfT/ZzrZP/YDCppvziMhxJqVizpHJ+OhERUU081LYxTn4wFKM7ewMALt8sxuiv/8YrPyZycOwtDCbVdPjCdWh1Alo2ckBTN67qSkREteMol+GL8V2wc2ZfdPBWQCcAWxOu4OEv/8LWhMsW373DYFJNB85dAwD0a9NY5EqIiMgctPNS4PeX+iLy2d4AgOsFpXjlx+P434oYXMgpELk68TCYVIMgCP8Ek7aNRK6GiIjMSe8Wbjj38XC8NrQt5DIpolNyMeizA/h8zzmL3KGYwaQazucUIDO/BHKZFL2au4pdDhERmRkbmRQvPNQKe1/th+5NXQAAS/adxyNfHcLxjDxxi2tgDCbVcOBsRWtJ7xZusLW2ErkaIiIyV36u9tg0LQhvj2gHd0cbXMgpwJhv/sb7v56C1kI2A2QwqYaY1IqV+vq04kqvRERUvyQSCZ7p2wJ7XumHkZ28oROAtYfT0PmDPyxi5g6DyX3odAKOplZME+7VnMGEiIgahouDDb6a0AULHgsAAKg15Xh06d9INPOuHQaT+zibrYaqpBz2Nlbo4K0QuxwiIrIw43s2wbqne6KRkxxZqhKM/vpvfBedJnZZ9YbB5D4qW0u6NXWBzIq/XERE1PD6tWmEfbP6IbhlRcv9h78lY9/pbJGrqh/8pr2Pf7pxOBuHiIjEo7C1xtrJPQEA5ToBL/6QYJbL2TOY3IMgCIi5FUx6cnwJERGJzEYmxYW5w9G3tTuKSrV4dOnfZrdTcY2CSUREBAIDA6FQKKBQKBAUFISdO3dW672RkZGQSCQYPXp0beoURer1Qlwv0MBGJkWgr1LscoiIiCCzkuLrkK7651/9eUHEagyvRsHE19cXCxYsQHx8POLi4jBgwACMGjUKp06duuf70tLSMHv2bPTt27dOxTa0ym6czn7OXL+EiIiMhsLWGssndgMArP07zaymEdcomIwcORIPP/wwWrdujTZt2mDu3LlwdHTEkSNH7voerVaLkJAQfPDBB2jRokWdC25Ix9JvAgB6NHMRuRIiIqKqhrT3QNcmzijV6rD2cJrY5RhMrceYaLVaREZGorCwEEFBQXc978MPP0Tjxo0xZcqUan+2RqOBSqWq8hDDicv5AIBOvs6iXJ+IiOhuJBIJnn2wJQDgx9gMlGl1IldkGLKaviEpKQlBQUEoKSmBo6Mjtm7divbt29/x3EOHDmHVqlVITEys0TXmz5+PDz74oKalGVRRaTnOZasBAJ38nEWthYiI6E4GtmsMNwcb5BaW4tCF63iobWOxS6qzGreYtG3bFomJiYiJicH06dMRFhaG5OTk285Tq9WYOHEiVqxYAXd39xpdIzw8HPn5+fpHRkZGTcuss+SrKugEoLGTHB4K2wa/PhER0f1YW0nxcIAXAOCPU1kiV2MYNW4xsbGxQatWrQAA3bp1Q2xsLJYsWYLly5dXOe/ixYtIS0vDyJEj9cd0uopmJplMhrNnz6Jly5Z3vIZcLodcLq9paQZ1/FY3TiC7cYiIyIj1a9MI649cQvTFXLFLMYgaB5P/0ul00Gg0tx339/dHUlJSlWNvv/021Go1lixZAj8/v7peul5VLlrTidOEiYjIiPVs4QqpBEjLLcLVvGJ4O9uJXVKd1CiYhIeHY/jw4WjSpAnUajU2btyIqKgo7N69GwAQGhoKHx8fzJ8/H7a2tujYsWOV9zs7OwPAbceNUeXA1wAGEyIiMmIKW2sE+DrjeEYeoi/mYmw3X7FLqpMaBZOcnByEhoYiMzMTSqUSgYGB2L17NwYPHgwASE9Ph1Rq+ovJ5heXIfV6IQB25RARkfELbumG4xl5OGxpwWTVqlX3fD0qKuqer69du7YmlxPNySsVrSV+rnZwdbARuRoiIqJ7C27phoioi4i+eB2CIEAikYhdUq2ZfvNGPTh+a3xJoI+zqHUQERFVR/emrrC2kuBqfgnSbxSJXU6dMJjcQZJ+Rg7HlxARkfGzs7FCF7+KVcoPm/jsHAaTOzh59dbAVx8GEyIiMg1BLd0AMJiYHXVJGTJuFAMA2nkpRK6GiIioeoJvBZPoi7kQBEHkamqPweQ/zmRVLEPvqbCFCwe+EhGRiejcxBlymRTXCzRIyzXdcSYMJv9xOrNiw8B2Xk4iV0JERFR9cpkV2ntXtPRXLhJqihhM/uOfYMJuHCIiMi2dbq29lZiRJ2oddcFg8h9nb3Xl+DOYEBGRienkVzFpo3L1clPEYPIvgiDgXHYBAMDfk105RERkWipXKz91NR/lWp24xdQSg8m/XM0vQYGmHNZWEjRzcxC7HCIiohpp7uYAJ7kMJWU6nM1Wi11OrTCY/Mu5W904zd0dYCPjLw0REZkWqVSCjrfW4Dp1VSVyNbXDb99/OXcrXbb2YDcOERGZpg63ZuYkM5iYvpRrFTsKt2rkKHIlREREtVM5ZfjUVdMcAMtg8i+p1yuCSYtGHF9CRESmqXK5i7NZapNcAZbB5F9SbgWT5u4MJkREZJpaNHKAlVQCVUk5ctQascupMQaTW1QlZbheUHEDGUyIiMhUyWVWaOpmD+CftblMCYPJLam3xpc0cpLDydZa5GqIiIhqr+2tSRznTHDKMIPJLansxiEiIjNRObv0/K1FQ00Jg8ktleNLWjCYEBGRiatsMTHFRdYYTG5hiwkREZmLNh4Vy16czza9mTkMJrekXq9o7mIwISIiU9fM3QHWVhIUlmpxJa9Y7HJqhMEEFZv3VQ5+bcHF1YiIyMRZW0nRwr2y1cS0xpkwmADIUWtQWKqFVAI0cbUXuxwiIqI6a32rO8fUxpkwmAC4lFsEAPB2tuPmfUREZBZMdcowv4UBXMmrCCa+LnYiV0JERGQYrRlMTNeVmxUDg3xd2I1DRETmoa1nRTC5kFMArc50ZuYwmAD6Ecs+zmwxISIi89DE1R5ymRQlZTpk3CgSu5xqYzABcPlWi4kPu3KIiMhMWEkl+iUwUnMLRa6m+hhM8K+uHLaYEBGRGWnmVhFMLl1nMDEZgiD805XDFhMiIjIjTd0rxk6m5bIrx2RcLyiFplwHiQTwUjKYEBGR+ahsMUljV47pqGwt8XCy5RomRERkVpq6VbSYXGKLiem4fLPiZrEbh4iIzE1li0nGjSKUa3UiV1M9Fh9MMvNKAFSs+kpERGROPBUVvQHlOgFXb33fGTuLDybXCjQAgMZOcpErISIiMiypVIKmrpUDYE1jnInFB5Pr6opg0ojBhIiIzFDTyinDDCamobLFpJEjgwkREZmfZm6mNWWYwYQtJkREZMaaurPFxKQwmBARkTlrrl/LhC0mRq9cq8ONolIADCZERGSeKtcySc8tMoldhi06mNwoLIUgAFIJ4GJvI3Y5REREBuftbAdrKwlKtTpkqYx/ynCNgklERAQCAwOhUCigUCgQFBSEnTt33vX8FStWoG/fvnBxcYGLiwsGDRqEo0eP1rloQ8m51Y3j5iiHlVQicjVERESGZyWVwO/WlGFT2MyvRsHE19cXCxYsQHx8POLi4jBgwACMGjUKp06duuP5UVFRmDBhAvbv34/o6Gj4+flhyJAhuHLlikGKr6vcwopuHDcHtpYQEZH5amZC40xkNTl55MiRVZ7PnTsXEREROHLkCDp06HDb+Rs2bKjyfOXKlfj555+xb98+hIaG3vU6Go0GGo1G/1ylUtWkzGpTFZcBABR21vXy+URERMagqZvpLLJW6zEmWq0WkZGRKCwsRFBQULXeU1RUhLKyMri6ut7zvPnz50OpVOoffn5+tS3zntQl5QAAhW2N8hkREZFJ8XOpCCZXbhaLXMn91TiYJCUlwdHREXK5HNOmTcPWrVvRvn37ar13zpw58Pb2xqBBg+55Xnh4OPLz8/WPjIyMmpZZLeqSihYTJ1u2mBARkfnyVNoCgEkMfq1xU0Hbtm2RmJiI/Px8bN68GWFhYThw4MB9w8mCBQsQGRmJqKgo2Nra3vNcuVwOubz+p+9Wtpg4scWEiIjMmIfiVjDJN8NgYmNjg1atWgEAunXrhtjYWCxZsgTLly+/63s+/fRTLFiwAHv37kVgYGDtqzWwyhYTBVtMiIjIjHndajHJVpVApxMgNeKZqHVex0Sn01UZqPpfn3zyCT766CPs2rUL3bt3r+vlDEqtqWgxcWSLCRERmbFGTnJIJEC5TtDPSDVWNfpGDg8Px/Dhw9GkSROo1Wps3LgRUVFR2L17NwAgNDQUPj4+mD9/PgBg4cKFePfdd7Fx40Y0a9YMWVlZAABHR0c4Ojoa+EepucoV8GRGnByJiIjqytpKCndHOa6pNchWlRj1auc1ajHJyclBaGgo2rZti4EDByI2Nha7d+/G4MGDAQDp6enIzMzUnx8REYHS0lKMGzcOXl5e+senn35q2J+ilipX5pVKGEyIiMi8VXbnZBr5OJMatZisWrXqnq9HRUVVeZ6WllbTehqUTqhIJmwwISIic1cxADbf6GfmWPReOcKtYCJhiwkREZk5z1szc7KNvMXEwoNJxX/ZYkJERObO00S6ciw6mOjYYkJERBZC32LCrhzjxcGvRERkKUxl9VeLDiYCB78SEZGF0AcTduUYL4EtJkREZCEqu3IKNOUouLXAqDGy6GBSOcYEzCVERGTmHOQyOMkrVgkx5lYTCw8mFf9liwkREVkCU+jOsfBgwjEmRERkOUxhAKxFBxOOMSEiIkviYQJThi07mKByHRORCyEiImoA/+yXUyxyJXdn0cFEp6v4L1tMiIjIElS2mGTla0Su5O4sO5gIbDEhIiLLYQqrv1p0MOEYEyIisiSmsF+OZQcTcFYOERFZjspgkluoQZlWJ3I1d2bRwaRyHRNu4kdERJbA1d4G1lYSCAKQozbOcSYWHkwqW0wYTIiIyPxJpZJ/DYA1zpk5MrELENO0fi1xvUADf08nsUshIiJqEJ4KW1y+WWy0M3MsOpgM7eApdglEREQNysPIV3+16K4cIiIiS+Nl5F05DCZEREQW5J/9coyzK4fBhIiIyILo98sx0rVMGEyIiIgsiLHvMMxgQkREZEHcHGwAADcKS0Wu5M4YTIiIiCyIm6McAFCgKUdJmVbkam7HYEJERGRBFLYyWFtVLCxqjK0mDCZEREQWRCKRwPVWd05uAYMJERERiczNoaI7J7fQ+KYMM5gQERFZGDdHtpgQERGRkaicmcMWEyIiIhKdq74rhy0mREREJDJ25RAREZHRMOZF1hhMiIiILEzlImu5BRxjQkRERCLTd+WwxYSIiIjE5sYF1oiIiMhYVHblFJdpUVRaLnI1VTGYEBERWRgHGyv9fjl5RWUiV1MVgwkREZGFkUgkUNpVdOcwmBAREZHonO2tAQB5xcY1zqRGwSQiIgKBgYFQKBRQKBQICgrCzp077/meTZs2wd/fH7a2tggICMCOHTvqVDARERHVnbNdRTDJN+UWE19fXyxYsADx8fGIi4vDgAEDMGrUKJw6deqO5x8+fBgTJkzAlClTkJCQgNGjR2P06NE4efKkQYonIiKi2vmnxcS4golEEAShLh/g6uqKRYsWYcqUKbe99uSTT6KwsBDbt2/XH+vduzc6d+6MZcuWVfsaKpUKSqUS+fn5UCgUdSmXiIiIAMz66Th+PnYZc4b5Y3r/lvVyjdp8f9d6jIlWq0VkZCQKCwsRFBR0x3Oio6MxaNCgKseGDh2K6Ojoe362RqOBSqWq8iAiIiLDMYsxJgCQlJQER0dHyOVyTJs2DVu3bkX79u3veG5WVhY8PDyqHPPw8EBWVtY9rzF//nwolUr9w8/Pr6ZlEhER0T2YxRgTAGjbti0SExMRExOD6dOnIywsDMnJyQYtKjw8HPn5+fpHRkaGQT+fiIjI0jk7GOd0YVlN32BjY4NWrVoBALp164bY2FgsWbIEy5cvv+1cT09PZGdnVzmWnZ0NT0/Pe15DLpdDLpfXtDQiIiKqpsoWE5PvyvkvnU4HjebOuxMGBQVh3759VY7t2bPnrmNSiIiIqGHox5iYcotJeHg4hg8fjiZNmkCtVmPjxo2IiorC7t27AQChoaHw8fHB/PnzAQAzZ85Ev379sHjxYowYMQKRkZGIi4vDt99+a/ifhIiIiKrN+dbKr/lGNl24RsEkJycHoaGhyMzMhFKpRGBgIHbv3o3BgwcDANLT0yGV/tMIExwcjI0bN+Ltt9/Gm2++idatW+OXX35Bx44dDftTEBERUY0Ya4tJndcxaQhcx4SIiMiwVCVlCHz/DwDAmY+GwdbayvDXaMh1TIiIiMh0OcllsJJW7DCsMqLuHAYTIiIiC1Sxw3BFd85NI+rOYTAhIiKyUArbiqGm6hIGEyIiIhKZ4laLiYrBhIiIiMSmsL0VTIrLRa7kHwwmREREFsrpVlcOW0yIiIhIdJUtJuoStpgQERGRyBR2t1pMOF2YiIiIxKYfY8KuHCIiIhKbflYOB78SERGR2PRdOWwxISIiIrE5yStbTBhMiIiISGSVXTmclUNERESiq+zKyWeLCREREYnNqXIdEw1bTIiIiEhkjvKKFpPSch005VqRq6nAYEJERGShKoMJABQYyTgTBhMiIiILZSWV6P//WoFGxEr+wWBCREREOJGRL3YJABhMiIiILNr0/i3RzkuBRzp5iV0KAEAiCIIgdhH3o1KpoFQqkZ+fD4VCIXY5REREVA21+f5miwkREREZDQYTIiIiMhoMJkRERGQ0GEyIiIjIaDCYEBERkdFgMCEiIiKjwWBCRERERoPBhIiIiIwGgwkREREZDQYTIiIiMhoMJkRERGQ0GEyIiIjIaDCYEBERkdFgMCEiIiKjIRO7gOoQBAFAxfbJREREZBoqv7crv8erwySCiVqtBgD4+fmJXAkRERHVlFqthlKprNa5EqEmMUYkOp0OV69ehZOTEyQSidjlGAWVSgU/Pz9kZGRAoVCIXQ79B++PceP9MW68P8avuvdIEASo1Wp4e3tDKq3e6BGTaDGRSqXw9fUVuwyjpFAo+AfXiPH+GDfeH+PG+2P8qnOPqttSUomDX4mIiMhoMJgQERGR0WAwMVFyuRzvvfce5HK52KXQHfD+GDfeH+PG+2P86vMemcTgVyIiIrIMbDEhIiIio8FgQkREREaDwYSIiIiMBoMJERERGQ0GEyN35coVPPXUU3Bzc4OdnR0CAgIQFxenf10QBLz77rvw8vKCnZ0dBg0ahPPnz4tYseXQarV455130Lx5c9jZ2aFly5b46KOPquwJwfvTsA4ePIiRI0fC29sbEokEv/zyS5XXq3M/bty4gZCQECgUCjg7O2PKlCkoKChowJ/CfN3r/pSVlWHOnDkICAiAg4MDvL29ERoaiqtXr1b5DN6f+nO/Pz//Nm3aNEgkEnzxxRdVjhvi/jCYGLGbN2+iT58+sLa2xs6dO5GcnIzFixfDxcVFf84nn3yCL7/8EsuWLUNMTAwcHBwwdOhQlJSUiFi5ZVi4cCEiIiKwdOlSnD59GgsXLsQnn3yCr776Sn8O70/DKiwsRKdOnfD111/f8fXq3I+QkBCcOnUKe/bswfbt23Hw4EE8++yzDfUjmLV73Z+ioiIcO3YM77zzDo4dO4YtW7bg7NmzePTRR6ucx/tTf+7356fS1q1bceTIEXh7e9/2mkHuj0BGa86cOcIDDzxw19d1Op3g6ekpLFq0SH8sLy9PkMvlwg8//NAQJVq0ESNGCE8//XSVY4899pgQEhIiCALvj9gACFu3btU/r879SE5OFgAIsbGx+nN27twpSCQS4cqVKw1WuyX47/25k6NHjwoAhEuXLgmCwPvTkO52fy5fviz4+PgIJ0+eFJo2bSp8/vnn+tcMdX/YYmLEfv31V3Tv3h2PP/44GjdujC5dumDFihX611NTU5GVlYVBgwbpjymVSvTq1QvR0dFilGxRgoODsW/fPpw7dw4AcPz4cRw6dAjDhw8HwPtjbKpzP6Kjo+Hs7Izu3bvrzxk0aBCkUiliYmIavGZLl5+fD4lEAmdnZwC8P2LT6XSYOHEiXnvtNXTo0OG21w11f0xiEz9LlZKSgoiICLz66qt48803ERsbi5deegk2NjYICwtDVlYWAMDDw6PK+zw8PPSvUf154403oFKp4O/vDysrK2i1WsydOxchISEAwPtjZKpzP7KystC4ceMqr8tkMri6uvKeNbCSkhLMmTMHEyZM0G8Sx/sjroULF0Imk+Gll1664+uGuj8MJkZMp9Ohe/fumDdvHgCgS5cuOHnyJJYtW4awsDCRq6OffvoJGzZswMaNG9GhQwckJibi5Zdfhre3N+8PUR2UlZXhiSeegCAIiIiIELscAhAfH48lS5bg2LFjkEgk9XotduUYMS8vL7Rv377KsXbt2iE9PR0A4OnpCQDIzs6uck52drb+Nao/r732Gt544w2MHz8eAQEBmDhxIl555RXMnz8fAO+PsanO/fD09EROTk6V18vLy3Hjxg3eswZSGUouXbqEPXv26FtLAN4fMf3111/IyclBkyZNIJPJIJPJcOnSJcyaNQvNmjUDYLj7w2BixPr06YOzZ89WOXbu3Dk0bdoUANC8eXN4enpi3759+tdVKhViYmIQFBTUoLVaoqKiIkilVf8IWVlZQafTAeD9MTbVuR9BQUHIy8tDfHy8/pw///wTOp0OvXr1avCaLU1lKDl//jz27t0LNze3Kq/z/ohn4sSJOHHiBBITE/UPb29vvPbaa9i9ezcAA96fOgzapXp29OhRQSaTCXPnzhXOnz8vbNiwQbC3txe+//57/TkLFiwQnJ2dhW3btgknTpwQRo0aJTRv3lwoLi4WsXLLEBYWJvj4+Ajbt28XUlNThS1btgju7u7C66+/rj+H96dhqdVqISEhQUhISBAACJ999pmQkJCgn9VRnfsxbNgwoUuXLkJMTIxw6NAhoXXr1sKECRPE+pHMyr3uT2lpqfDoo48Kvr6+QmJiopCZmal/aDQa/Wfw/tSf+/35+a//zsoRBMPcHwYTI/fbb78JHTt2FORyueDv7y98++23VV7X6XTCO++8I3h4eAhyuVwYOHCgcPbsWZGqtSwqlUqYOXOm0KRJE8HW1lZo0aKF8NZbb1X5S5T3p2Ht379fAHDbIywsTBCE6t2P3NxcYcKECYKjo6OgUCiEyZMnC2q1WoSfxvzc6/6kpqbe8TUAwv79+/WfwftTf+735+e/7hRMDHF/JILwr2UqiYiIiETEMSZERERkNBhMiIiIyGgwmBAREZHRYDAhIiIio8FgQkREREaDwYSIiIiMBoMJERERGQ0GEyIiIjIaDCZEVCdRUVGQSCTIy8tr8Gv3798fL7/8coNfl4jqD4MJEdXIf8NAcHAwMjMzoVQqG+T6Bw4cgJ+fX4Nci4gankzsAojItNnY2DTolvPbtm3DyJEjG+x6RNSw2GJCRNU2adIkHDhwAEuWLIFEIoFEIsHatWurdOWsXbsWzs7O2L59O9q2bQt7e3uMGzcORUVFWLduHZo1awYXFxe89NJL0Gq1+s/WaDSYPXs2fHx84ODggF69eiEqKuq2Gn799Vc8+uij+uc6nQ6vv/46XF1d4enpiffff7+efxWIqD6xxYSIqm3JkiU4d+4cOnbsiA8//BAAcOrUqdvOKyoqwpdffonIyEio1Wo89thjGDNmDJydnbFjxw6kpKRg7Nix6NOnD5588kkAwIwZM5CcnIzIyEh4e3tj69atGDZsGJKSktC6dWv9tXJycjBgwAD9tdatW4dXX30VMTExiI6OxqRJk9CnTx8MHjy4AX5FiMjQGEyIqNqUSiVsbGxgb2+v7745c+bMbeeVlZUhIiICLVu2BACMGzcO69evR3Z2NhwdHdG+fXs89NBD2L9/P5588kmkp6djzZo1SE9Ph7e3NwBg9uzZ2LVrF9asWYN58+YBqOjGGTp0KGxsbPTXCgwMxHvvvQcAaN26NZYuXYp9+/YxmBCZKAYTIjI4e3t7fSgBAA8PDzRr1gyOjo5VjuXk5AAAkpKSoNVq0aZNmyqfo9Fo4Obmpn++bds2zJgxo8o5gYGBVZ57eXnpP5eITA+DCREZnLW1dZXnEonkjsd0Oh0AoKCgAFZWVoiPj4eVlVWV8yrDTGZmJhISEjBixIj7Xqvyc4nI9DCYEFGN2NjYVBm0aghdunSBVqtFTk4O+vbte8dzfvvtNwQHB8PV1dWg1yYi48JZOURUI82aNUNMTAzS0tJw/fp1g7ROtGnTBiEhIQgNDcWWLVuQmpqKo0ePYv78+fj9998B3D4bh4jME4MJEdXI7NmzYWVlhfbt26NRo0ZIT083yOeuWbMGoaGhmDVrFtq2bYvRo0cjNjYWTZo0QWFhIfbt28dgQmQBJIIgCGIXQUR0L1u2bMHbb7+N5ORksUshonrGFhMiMnqOjo5YuHCh2GUQUQNgiwkREREZDbaYEBERkdFgMCEiIiKjwWBCRERERoPBhIiIiIwGgwkREREZDQYTIiIiMhoMJkRERGQ0GEyIiIjIaDCYEBERkdH4Pyvfw6/K7wrbAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ocv_df.plot(y=\"Ewe-Ece/V\", x=\"time/h\", kind=\"line\")" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "As we are aiming to identify the resistance elements in an equivalent circuit model, we will construct an OCV function from the discharge portion of the experimental data. This is done by filtering the dataframe and appending an additional point to ensure the OCV function has data across the operating region (i.e. up to 4.2V)." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "sto_measured = (\n", " 1\n", " - ocv_df.loc[(ocv_df[\"I/mA\"] <= 0.0), \"Capacity/mA.h\"]\n", " / ocv_df.loc[(ocv_df[\"I/mA\"] <= 0.0), \"Capacity/mA.h\"].iloc[-1]\n", ")\n", "V_measured = ocv_df.loc[(ocv_df[\"I/mA\"] <= 0.0), \"Ewe-Ece/V\"]\n", "V_measured.iloc[0] = 4.2 # Extend for improved interpolation (a bit of a fudge)" ] }, { "cell_type": "code", "execution_count": null, "id": "8", "metadata": {}, "outputs": [], "source": [ "def ocv(sto):\n", " name = \"OCV\"\n", " x = np.flip(sto_measured.to_numpy())\n", " y = np.flip(V_measured.to_numpy())\n", " return pybamm.Interpolant(x, y, sto, name)" ] }, { "cell_type": "markdown", "id": "9", "metadata": {}, "source": [ "Next, let's construct the parameter set and update it with the known information and placeholder values for this cell." ] }, { "cell_type": "code", "execution_count": null, "id": "10", "metadata": {}, "outputs": [], "source": [ "parameter_set = pybop.ParameterSet.pybamm(\"ECM_Example\")\n", "parameter_set.update(\n", " {\n", " \"Cell capacity [A.h]\": 22.651, # 083/828 - C/20\n", " \"Nominal cell capacity [A.h]\": 22.651,\n", " \"Current function [A]\": 22.651,\n", " \"Initial SoC\": 1.0,\n", " \"Upper voltage cut-off [V]\": 4.25, # Extended to avoid hitting event\n", " \"Lower voltage cut-off [V]\": 2.5,\n", " \"Open-circuit voltage [V]\": ocv,\n", " \"R2 [Ohm]\": 1e-4, # placeholder\n", " \"C2 [F]\": 4e5, # placeholder\n", " \"Element-2 initial overpotential [V]\": 0,\n", " },\n", " check_already_exists=False,\n", ")" ] }, { "cell_type": "markdown", "id": "11", "metadata": {}, "source": [ "## Import Cycling Data\n", "Now we import the corresponding constant current tests for this cell and construct a dataset from a subset of this time-series. This is done to reduce the computational time for inference, but also to improve the robustness of the inference task." ] }, { "cell_type": "code", "execution_count": null, "id": "12", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAG0CAYAAAARqnxaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKt0lEQVR4nO3deVxU5f4H8M+ZGWZYh03ZVwVFRRQUF9wqySXMLOual65my23Bn1rdUrstt3stqMxrWpnaTb1l2uqSpcV1Qc0NQRRERUUFERgQ2WWbOb8/0LFJUZaBMzN83q/XvF5y5uHMdx59zXx8zvM8RxBFUQQRERGRCZNJXQARERHRnTCwEBERkcljYCEiIiKTx8BCREREJo+BhYiIiEweAwsRERGZPAYWIiIiMnkMLERERGTyGFiIiIjI5DGwEBERkclrU2BJSEiAIAiYM2dOk22OHz+OyZMnIyAgAIIgYPHixTe10Wq1eP311xEYGAgbGxt0794d//rXv8C7BhAREREAKFr7i8nJyVi+fDnCwsJu2666uhrdunXDI488ghdeeOGWbd59910sW7YMa9asQZ8+fXD48GHMmDEDjo6OmDVrVrPq0el0uHTpEhwcHCAIQovfDxEREXU8URRRUVEBLy8vyGRNj6O0KrBUVlYiNjYWK1euxIIFC27bNjIyEpGRkQCAefPm3bLNvn378MADDyAmJgYAEBAQgHXr1uHQoUPNrunSpUvw9fVtdnsiIiIyHbm5ufDx8Wny+VYFlri4OMTExCA6OvqOgaU5oqKisGLFCmRlZaFHjx44evQo9u7di0WLFjX5O7W1taitrdX/fP3yUW5uLtRqdZtrIiIiovZXXl4OX19fODg43LZdiwPL+vXrkZqaiuTk5FYX90fz5s1DeXk5QkJCIJfLodVq8fbbbyM2NrbJ34mPj8dbb71103G1Ws3AQkREZGbuNJ2jRZNuc3NzMXv2bKxduxbW1tZtKuz3vvnmG6xduxZfffUVUlNTsWbNGixcuBBr1qxp8nfmz5+PsrIy/SM3N9do9RAREZFpadEIS0pKCjQaDSIiIvTHtFotdu/ejY8++gi1tbWQy+UtLuLll1/GvHnz8OijjwIA+vbtiwsXLiA+Ph7Tp0+/5e+oVCqoVKoWvxYRERGZnxYFltGjRyM9Pd3g2IwZMxASEoK5c+e2KqwAjSuJ/jgzWC6XQ6fTtep8REREZFlaFFgcHBwQGhpqcMzOzg6urq7649OmTYO3tzfi4+MBAHV1dcjMzNT/OS8vD2lpabC3t0dQUBAA4P7778fbb78NPz8/9OnTB0eOHMGiRYvwxBNPtPkNEhERXafValFfXy91GZ2KlZVVqwc0fq/V+7A0JScnx2C05NKlSwgPD9f/vHDhQixcuBCjRo3Crl27AABLly7F66+/jueffx4ajQZeXl545pln8MYbbxi7PCIi6oREUURBQQFKS0ulLqVTcnJygoeHR5v2SRNEC9lOtry8HI6OjigrK+MqISIiMpCfn4/S0lK4ubnB1taWG4x2EFEUUV1dDY1GAycnJ3h6et7Uprnf30YfYSEiIjIlWq1WH1ZcXV2lLqfTsbGxAQBoNBq4ubm1+vIQb35IREQW7fqcFVtbW4kr6byu931b5g8xsBARUafAy0DSMUbfM7AQERGRyWNgISIiIpPHwEJERNRJBQQEYPHixfqfBUHAxo0bJavndhhY7qDsaj1OF1agXstdd4mIqGM9/vjjEAThpse4ceM6rIarV6/Czs4O7777LpydnVFTU3NTm+rqaqjVaixZsqTd6mBguYOo+O2499+7cfHKValLISKiTmjcuHHIz883eKxbt67DXj8xMRH+/v6YMWMGqqqq8MMPP9zU5rvvvkNdXR0ee+yxdquDgeUO3NSNd6UuLL85URIRkfkRRRHVdQ2SPFqzV6tKpYKHh4fBw9nZGQBQWlqKZ555Bu7u7rC2tkZoaCi2bNmi/929e/dixIgRsLGxga+vL2bNmoWqqqoWvf6mTZswceJEuLm54f7778fnn39+U5vPP/8ckyZNgouLS4vfX3Nx47g7cHNQ4VxxFQMLEZGFuFqvRe83fpHktTP/ORa2SuN89ep0OowfPx4VFRX48ssv0b17d2RmZuo3Zjt79izGjRuHBQsW4PPPP0dRURFmzpyJmTNnYtWqVc1+jS1btujntTz55JOYMGECLly4AH9/fwBAdnY2du/ejV9+ad8+ZWC5A/drIyya8lqJKyEios5oy5YtsLe3Nzj26quvYuDAgTh06BBOnDiBHj16AAC6deumbxMfH4/Y2FjMmTMHABAcHIwlS5Zg1KhRWLZsGaytre/42gcOHAAADB48GAAwduxYeHl5YdWqVfjHP/4BAFi9ejV8fX0xevTotr7V22JguQN3tQoALwkREVkKGys5Mv85VrLXbqm7774by5YtMzjm4uKCzz77DD4+Pvqw8kdHjx7FsWPHsHbtWv0xURSh0+lw7tw59OrV646vvWnTJkyYMEF/U2O5XI7p06dj9erVePPNNyGKItasWYMZM2YY3Pi4PTCw3IF+hKWCIyxERJZAEASjXZbpCHZ2dggKCrrp+PV79DSlsrISzzzzDGbNmnXTc35+fs167c2bNyMhIcHg2BNPPIH4+Hjs2LEDOp0Oubm5mDFjRrPO1xbm8zcmEU66JSIiUxQWFoaLFy8iKyvrlqMsERERyMzMvGXYaY7Tp0/jwoULuPfeew2Od+/eHaNGjcLnn38OURQRHR2tn8/SnhhY7sDdofGSEEdYiIhICrW1tSgoKDA4plAoMGrUKIwcORKTJ0/GokWLEBQUhJMnT+r3aZk7dy6GDBmCmTNn4qmnnoKdnR0yMzORmJiIjz766I6vu2nTJkRHR9/yppFPPvkknn76aQCNc1g6Apc134H770ZYWrMcjYiIqC22bdsGT09Pg8fw4cMBAN9//z0iIyMxdepU9O7dG6+88gq0Wi2AxhGYpKQkZGVlYcSIEQgPD8cbb7wBLy+vZr3u9eXMtzJ58mSoVCrY2tpi0qRJRnmfdyKIFvItXF5eDkdHR5SVlUGtVhvtvNV1Dfrlb+n/GAMHayujnZuIiNpfTU0Nzp07h8DAwGatjCGguLgYnp6euHjxItzd3dt8vtv9HTT3+5sjLHdgq1TAQdV45ayQS5uJiKgTKCkpwaJFi4wSVoyFc1iawU2tQkVRAzTlNQhys7/zLxAREZmxHj16NLlcWiocYWkGLm0mIiKSFgNLM7hzaTMRkdmzkCmbZskYfc/A0gxu+t1uOcJCRGRurKwaF0tUV1dLXEnndb3vr/9dtAbnsDSDu8O1EZYKjrAQEZkbuVwOJycnaDQaAICtrS0EQZC4qs5BFEVUV1dDo9HAyclJf2PG1mBgaYbrIywaXhIiIjJLHh4eAKAPLdSxnJyc9H8HrcXA0gw35rDwkhARkTkSBAGenp5wc3NDfX291OV0KlZWVm0aWbmOgaUZ9JeEru12y6FEIiLzJJfLjfLlSR2Pk26b4folodoGHcprGiSuhoiIqPNhYGkGays5HG0aZzZzHgsREVHHY2BpJncubSYiIpIMA0szcfM4IiIi6TCwNJMb92IhIiKSDANLM93Yi4WXhIiIiDoaA0szuTtcCywcYSEiIupwDCzNxM3jiIiIpMPA0kxunHRLREQkGQaWZnL/3RwW3qKciIioYzGwNFPXa3NY6rQ6lFbzPhREREQdiYGlmVQKOVzslAC4tJmIiKijMbC0gJsDd7slIiKSAgNLC1yfeMv7CREREXUsBpYWuLEXC0dYiIiIOhIDSwvwfkJERETSYGBpgRt3bGZgISIi6kgMLC3gxt1uiYiIJMHA0gLunHRLREQkCQaWFtDvdltRC52Ou90SERF1FAaWFuhir4IgAA06EVeq66Quh4iIqNNgYGkBK7kMrtd3u+U8FiIiog7DwNJCbg7XJt5ye34iIqIOw8DSQjfu2szAQkRE1FEYWFrInUubiYiIOhwDSwu5cbdbIiKiDsfA0kI3drvlCAsREVFHaVNgSUhIgCAImDNnTpNtjh8/jsmTJyMgIACCIGDx4sW3bJeXl4fHHnsMrq6usLGxQd++fXH48OG2lNcurk+6LeKkWyIiog7T6sCSnJyM5cuXIyws7Lbtqqur0a1bNyQkJMDDw+OWba5cuYJhw4bBysoKW7duRWZmJj744AM4Ozu3trx2wxEWIiKijqdozS9VVlYiNjYWK1euxIIFC27bNjIyEpGRkQCAefPm3bLNu+++C19fX6xatUp/LDAwsDWltbvrk26LKmuh1YmQywSJKyIiIrJ8rRphiYuLQ0xMDKKjo41SxObNmzFw4EA88sgjcHNzQ3h4OFauXHnb36mtrUV5ebnBoyO42ikhEwCtTsTlKo6yEBERdYQWB5b169cjNTUV8fHxRisiOzsby5YtQ3BwMH755Rc899xzmDVrFtasWdPk78THx8PR0VH/8PX1NVo9t6OQy9DF/vpeLAwsREREHaFFgSU3NxezZ8/G2rVrYW1tbbQidDodIiIi8M477yA8PBx//etf8fTTT+PTTz9t8nfmz5+PsrIy/SM3N9do9dyJO5c2ExERdagWBZaUlBRoNBpERERAoVBAoVAgKSkJS5YsgUKhgFarbVURnp6e6N27t8GxXr16IScnp8nfUalUUKvVBo+Owom3REREHatFk25Hjx6N9PR0g2MzZsxASEgI5s6dC7lc3qoihg0bhlOnThkcy8rKgr+/f6vO196ubx6n4dJmIiKiDtGiwOLg4IDQ0FCDY3Z2dnB1ddUfnzZtGry9vfVzXOrq6pCZman/c15eHtLS0mBvb4+goCAAwAsvvICoqCi88847+NOf/oRDhw5hxYoVWLFiRZvfYHtwc+AICxERUUcy+k63OTk5yM/P1/986dIlhIeHIzw8HPn5+Vi4cCHCw8Px1FNP6dtERkZiw4YNWLduHUJDQ/Gvf/0LixcvRmxsrLHLM4rrc1h4A0QiIqKO0ap9WH5v165dt/05ICAAoije8TwTJkzAhAkT2lpOh9DPYeElISIiog7Bewm1wvXt+XlJiIiIqGMwsLTC9UtCxZW1aNDqJK6GiIjI8jGwtIKrnRJymQBRBIor66Quh4iIyOIxsLSCTCboVwpxaTMREVH7Y2BpJTc157EQERF1FAaWVrqxFwtHWIiIiNobA0srXV/azL1YiIiI2h8DSyu5c2kzERFRh2FgaSX9HZs56ZaIiKjdMbC0kpv+khBHWIiIiNobA0srufOOzURERB2GgaWVbux2W4d67nZLRETUrhhYWsnZ1gpWcgEAUFTBy0JERETtiYGllQRB+N1NEHlZiIiIqD0xsLTB9Ym3XNpMRETUvhhY2uD6XiyceEtERNS+GFjawNvZBgBwMLtE4kqIiIgsGwNLGzw8wAcA8FN6Pk4WlEtcDRERkeViYGmDXp5qxPT1BAAsTjwtcTVERESWi4GljeZEB0MQgG3HC5CRVyZ1OURERBaJgaWNgt0dMLGfFwBg8f84ykJERNQeGFiMYNboYMgE4H8nCnE0t1TqcoiIiCwOA4sRdO9qjwfDGyfg/vt/WRJXQ0REZHkYWIxk1uggyGUCdp0qQsqFK1KXQ0REZFEYWIzE39UOD0c0jrIs5igLERGRUTGwGNHMe4JgJRew53QxDp3jZnJERETGwsBiRL4utvjTQF8AwKLEUxJXQ0REZDkYWIws7u4gKOUyHMguwb6zxVKXQ0REZBEYWIzMy8kGUwc1jrL8OzELoihKXBEREZH5Y2BpB8/fHQSVQobk81ew9wxHWYiIiNqKgaUduKut8dgQfwDAwl85ykJERNRWDCzt5NlR3WGnlONobik2pV2SuhwiIiKzxsDSTro6qPD83UEAgIStJ1Fd1yBxRUREROaLgaUdPTk8ED7ONigor8GnSdlSl0NERGS2GFjakbWVHK/e1wsAsDzpLPJKr0pcERERkXliYGln40M9MCjQBbUNOry79aTU5RAREZklBpZ2JggC3pjQG4IAbD56CSkXuGU/ERFRSzGwdIBQb0dMubZl/1s/ZkKn4zJnIiKilmBg6SAvjekJe5UCxy6W4YcjeVKXQ0REZFYYWDpIVwcVZt7TuMz5vW0nUVXLZc5ERETNxcDSgWYMC4C/qy00FbVYtuus1OUQERGZDQaWDqRS3FjmvGJPNnJLqiWuiIiIyDwwsHSwMb3dEdXdFXUNOiRwmTMREVGzMLB0MEEQ8PqE3pAJwE/p+djHuzkTERHdEQOLBHp5qhE7uPFuzvM3pKOmXitxRURERKaNgUUir4zrCQ+1NS5crsa//5cldTlEREQmjYFFIg7WVlgwKRQA8Nmec8jIK5O4IiIiItPFwCKh6N7umBDmCa1OxCvfHUO9Vid1SURERCaJgUVib97fB442VsjML8dne85JXQ4REZFJYmCRWFcHFV6f0BsAsPh/WThXXCVxRURERKaHgcUETI7wxojgLqht0GHe98d4c0QiIqI/YGAxAYIg4J0H+8LGSo6D50rw9eFcqUsiIiIyKQwsJsLXxRYvjekBAHjn5xMoLK+RuCIiIiLT0abAkpCQAEEQMGfOnCbbHD9+HJMnT0ZAQAAEQcDixYvbfE5LNWNYIPr5OKKipgFvbMqQuhwiIiKT0erAkpycjOXLlyMsLOy27aqrq9GtWzckJCTAw8PDKOe0VHKZgITJYVDIBPxyvBBb0/OlLomIiMgktCqwVFZWIjY2FitXroSzs/Nt20ZGRuL999/Ho48+CpVKZZRzWrJenmo8d1d3AMDfN2ZAU8FLQ0RERK0KLHFxcYiJiUF0dLTRCmnpOWtra1FeXm7wsBQz7wlCb081Sqrq8PK3xyCKXDVERESdW4sDy/r165Gamor4+HijFdGac8bHx8PR0VH/8PX1NVo9UlMp5Pjw0f5QKWRIyirCFwcuSF0SERGRpFoUWHJzczF79mysXbsW1tbWRimgteecP38+ysrK9I/cXMtaChzs7oD540MAAG//dAKnCyskroiIiEg6LQosKSkp0Gg0iIiIgEKhgEKhQFJSEpYsWQKFQgGtVtviAlp7TpVKBbVabfCwNNOjAjCqR1fUNugwe30a6hp4ryEiIuqcFC1pPHr0aKSnpxscmzFjBkJCQjB37lzI5fIWF9Ae57QUgiDg/YfDMO7DPcjML8cHiacwf3wvqcsiIiLqcC0KLA4ODggNDTU4ZmdnB1dXV/3xadOmwdvbWz8fpa6uDpmZmfo/5+XlIS0tDfb29ggKCmrWOTszN7U14h/qi2e+SMGK3dm4q4cbhnZ3lbosIiKiDmX0nW5zcnKQn39j/5BLly4hPDwc4eHhyM/Px8KFCxEeHo6nnnrK2C9tscb28cCjkb4QReClb9JQVl0vdUlEREQdShAtZM1seXk5HB0dUVZWZpHzWapqGxCzZA/OX67GxH5eWDI1XOqSiIiI2qy539+8l5CZsFMp8O8p/SGXCdh89BI2peVJXRIREVGHYWAxI+F+zpg9OhgA8NqGDJwrrpK4IiIioo7BwGJmnr+rOyIDnFFR24BnvjiMqtoGqUsiIiJqdwwsZkYhl+HjP0fAzUGFrMJKvPI9t+4nIiLLx8BihtzU1lj2WASs5AJ+OpaPlXuypS6JiIioXTGwmKkB/i544/4+AICErSex70yxxBURERG1HwYWM/bYYD88PMAHOhGYue4I8kqvSl0SERFRu2BgMWOCIGDBpFCEeqtRUlWHZ79IQU19y+/nREREZOoYWMyctZUcnz42AM62VkjPK8PrGzM4CZeIiCwOA4sF8HG2xdKpEZAJwLcpF7H2YI7UJRERERkVA4uFGB7cBa+MCwEAvPXjcaRcKJG4IiIiIuNhYLEgz4zshpi+nqjXinjmi1Tkl3ESLhERWQYGFgsiCALeezgMIR4OKK6sxTOchEtERBaCgcXC2KkUWDltIJxsrXDsYhnm/5DOSbhERGT2GFgskK+LLT75cwTkMgEbjuThP3vPSV0SERFRmzCwWKiooC54PaYXAOCdn09gd1aRxBURERG1HgOLBZseFYA/Dby2E+5XqThfXCV1SURERK3CwGLBBEHAvyaFItzPCeU1DXjqv4dRUVMvdVlEREQtxsBi4VQKOZY/NgDuahXOaCrxwtdp0Ok4CZeIiMwLA0sn4Ka2xvK/DIRSIcP/TmjwQeIpqUsiIiJqEQaWTqK/rxMSHuoLAPh451msPXhB4oqIiIiaj4GlE3kowgez7gkCALy+MQO/HC+QuCIiIqLmYWDpZF64twcejfSFTgT+b90RJJ/nPYeIiMj0MbB0MoIgYMGkUET3ckddgw5Prk5GVmGF1GURERHdFgNLJ6SQy7B0ajgG+DujvKYB0z8/hEulvFEiERGZLgaWTspGKcd/pg9EkJs98stqMP3zQyir5h4tRERkmhhYOjEnWyXWPDEIHmprnNZU4qn/JvPuzkREZJIYWDo5bycbrHliENTWCiSfv4L/W3cEDVqd1GUREREZYGAh9PRwwGfTI6FUyJCYWYhXN6RDFLkbLhERmQ4GFgIADAp0wdKp4ZAJwDeHL+Kdn08wtBARkclgYCG9sX088O7kMADAyj3n8MmusxJXRERE1IiBhQw8MtAXr0/oDQB4/5dT+PIAt/AnIiLpMbDQTZ4cHoj/u76F/6YM/Hj0ksQVERFRZ8fAQrf04r098Jch/hBF4IWv07DrlEbqkoiIqBNjYKFbEgQBb03sgwf6e6FBJ+LZL1NwmPcdIiIiiTCwUJNkMgELH+mHe0LcUFOvw4zVyci8VC51WURE1AkxsNBtWcll+PjPEYgMcEZFTQMe+89BnCxgaCEioo7FwEJ3ZKOU4z+PR6KfjyNKqurw55UHcaqAd3gmIqKOw8BCzaK2tsJ/nxyMMH1oOYCsQoYWIiLqGAws1GyONlb44onBCPVW4/K10HKaoYWIiDoAAwu1iKOtFb58cjB6e6pRXFmHqSsP4oymUuqyiIjIwjGwUIs52Sqx9qnB6OWpRnFlLaauPICzRQwtRETUfhhYqFWc7RpDS4iHA4oqajF1xQFkM7QQEVE7YWChVnOxU+Krp4cgxMMBmorGkRaGFiIiag8MLNQmLtdGWnq6O6CwvBaPrjjAOS1ERGR0DCzUZq72Knz19GCDkZYzGq4eIiIi42FgIaNoDC1D9HNaHl1xkEueiYjIaBhYyGhc7JRY9/SQa0ueG0dauCMuEREZAwMLGdX11UN9vBr3afnzygO89xAREbUZAwsZ3fXQcmNH3IM4kc/QQkRErcfAQu3CyVaJtU8OQV/vG/ceysgrk7osIiIyUwws1G4cba3w5VOD0c/HEVeq6zF15QEcPl8idVlERGSGGFioXTnaWOGLpwZjUIALKmoa8Jf/HMLurCKpyyIiIjPDwELtTm1thTVPDMKoHl1xtV6LJ9ckY2t6vtRlERGRGWlTYElISIAgCJgzZ06TbY4fP47JkycjICAAgiBg8eLFN7WJj49HZGQkHBwc4ObmhkmTJuHUqVNtKY1MjI1SjpXTBiKmryfqtSLivkrFdykXpS6LiIjMRKsDS3JyMpYvX46wsLDbtquurka3bt2QkJAADw+PW7ZJSkpCXFwcDhw4gMTERNTX12PMmDGoqqpqbXlkgpQKGZZMDcefBvpAJwJ/+/YoVv92TuqyiIjIDCha80uVlZWIjY3FypUrsWDBgtu2jYyMRGRkJABg3rx5t2yzbds2g59Xr14NNzc3pKSkYOTIkbf8ndraWtTW1up/Li/nsllzIJcJeHdyGBysrfCfvefwjx8zUVHTgJn3BEEQBKnLIyIiE9WqEZa4uDjExMQgOjra2PUAAMrKGpe/uri4NNkmPj4ejo6O+oevr2+71ELGJwgCXovphTnRwQCADxKzsOCnE9DpRIkrIyIiU9XiwLJ+/XqkpqYiPj6+PeqBTqfDnDlzMGzYMISGhjbZbv78+SgrK9M/cnNz26Ueah+CIGBOdA+8PqE3AOA/e89h5rpU1NRrJa6MiIhMUYsuCeXm5mL27NlITEyEtbV1uxQUFxeHjIwM7N2797btVCoVVCpVu9RAHefJ4YHoYq/E3749ip/TC1BYfhArpw2Ei51S6tKIiMiEtGiEJSUlBRqNBhEREVAoFFAoFEhKSsKSJUugUCig1bbtf8czZ87Eli1bsHPnTvj4+LTpXGQ+Hujvjf8+MRhqawVSLlzBQ5/8hvPFnHBNREQ3tCiwjB49Gunp6UhLS9M/Bg4ciNjYWKSlpUEul7eqCFEUMXPmTGzYsAE7duxAYGBgq85D5mtod1f88HwUvJ1scP5yNR5atg8pF65IXRYREZmIFl0ScnBwuGleiZ2dHVxdXfXHp02bBm9vb/0cl7q6OmRmZur/nJeXh7S0NNjb2yMoKAhA42Wgr776Cps2bYKDgwMKCgoAAI6OjrCxsWnbOySzEeTmgA1xUXhy9WGk55XhzysPYPGU/hjf11Pq0oiISGJG3+k2JycH+fk3djG9dOkSwsPDER4ejvz8fCxcuBDh4eF46qmn9G2WLVuGsrIy3HXXXfD09NQ/vv76a2OXRybOzcEaXz8zBNG93FDboMPzX6Xisz3ZEEWuICIi6swE0UK+CcrLy+Ho6IiysjKo1Wqpy6E20upE/GPzcXxx4AIA4E8DffDPB0JhbdW6y45ERGSamvv9zXsJkUmSywT884E+eC2mF2QC8M3hi5iy4gAKymqkLo2IiCTAwEImSxAEPDWiG9Y8MQiONlY4mluKCUv34vD5EqlLIyKiDsbAQiZvRHBX/DhzOEI8HFBcWYupKw9g7cELUpdFREQdiIGFzIKfqy1+eD5Kf7fnv2/IwPwf0lHbwJ1xiYg6AwYWMhu2SgU++nM45o4LgSAA6w7lYOqKA9CUc14LEZGlY2AhsyIIAp67qztWPR4JtbUCqTmN81pSc7jJHBGRJWNgIbN0V083bJ45HD3c7aGpqMWjyw/g6+QcqcsiIqJ2wsBCZiugix1+eH4YxvXxQJ1Wh7nfp+P1jRmoa9BJXRoRERkZAwuZNXuVAp/ERuCle3tAEIAvDlxA7GcHUFRRK3VpRERkRAwsZPZkMgH/NzoYn00bCAeVAsnnr+D+pXtxNLdU6tKIiMhIGFjIYozu5Y6NM4ehW1c7FJTX4JHl+/Ht4VypyyIiIiNgYCGL0r2rPTbGDUN0LzfUNejw8nfHMP+HY6ip534tRETmjIGFLI7a2gor/jIQc6KDr+3XkouHP92H3JJqqUsjIqJWYmAhiySTCZgT3QOrZwyCs60VMvLKEbNkD7afKJS6NCIiagUGFrJoo3p0xZZZI9Df1wnlNQ14cs1hvLftJBq0XPpMRGROGFjI4nk72eCbZ4Zi+lB/AMAnu87iL/85xKXPRERmhIGFOgWlQoa3HgjFkqnhsFXKsT/7MmKW7EHy+RKpSyMiomZgYKFOZWI/L2yeOQxBbte29F9xAJ/tyYYoilKXRkREt8HAQp1OkJsDNsUNw8R+XtDqRCz46QSeX5uKipp6qUsjIqImMLBQp2SnUuDDR/vjnw/0gZVcwNaMAkz86DecLCiXujQiIroFBhbqtARBwLShAfjmmaHwcrTGueIqTPr4N3yfclHq0oiI6A8YWKjTC/dzxk+zRmBUj66oqdfhpW+PYv4P6dwdl4jIhDCwEAFwtlNi1eOReCG6x7XdcXMw6ePfkJFXJnVpREQEBhYiPZlMwOzoYKyZMQiudkqcLKjAAx//hkW/nkJdAzeaIyKSEgML0R+M7NEVv74wEjFhntDqRCzZcQYTP9rL0RYiIgkxsBDdgqu9Ch//OQKfxEZwtIWIyAQwsBDdxn19PTnaQkRkAhhYiO7gVqMtkz7+DUu3n+ZNFImIOggDC1EzXR9tua+vBxp0Ij5IzMKflu/HhctVUpdGRGTxGFiIWuD6aMu/p/SDg0qB1JxS3PfhHnyTnMv7ERERtSMGFqIWEgQBD4b7YOucERgU6IKqOi1e+f4YnvkiBSVVdVKXR0RkkRhYiFrJx9kW654egrnjQmAlF/BrZiHGLt6Nnac0UpdGRGRxGFiI2kAuE/DcXd2x4flhCHKzR1FFLWasSsarG9JRWdsgdXlERBaDgYXICEK9HbHl/4bj8agAAMBXB3Mw9t+78duZYmkLIyKyEAwsREZibSXHPyb2wVdPDYaPsw3ySq8i9rOD+DtHW4iI2oyBhcjIooK64Jc5I/GXIf4AgLUcbSEiajMGFqJ2YKdS4F+TQjnaQkRkJAwsRO2oqdGWpKwiiSsjIjIvDCxE7exWoy3TPz+El745itJq7ttCRNQcDCxEHeT6aMuMYQEQBOD71IuIXrQbW9PzpS6NiMjkMbAQdSA7lQJv3t8H3z0bhSA3exRX1uK5tal49osUaMprpC6PiMhkMbAQSWCAvzN+mjUcs+4JgkImYNvxAkQvSsI3h3lPIiKiW2FgIZKISiHHi2N64sf/G46+3o4or2nAK98dw6SPf8P+s5elLo+IyKQIooX8d668vByOjo4oKyuDWq2WuhyiFmnQ6vD5b+fw4f9Oo6pOCwAYHeKGueND0MPdQeLqiIjaT3O/vxlYiExIcWUtlmw/ja8O5qBBJ0ImAH8a6IsX7u0Bd7W11OURERkdAwuRGcsuqsT7v5zC1owCAIC1lQxPj+iGZ0Z1h71KIXF1RETGw8BCZAFSLpTgnZ9PIuXCFQCAh9oabz3QB2P7eEhcGRGRcTCwEFkIURTxa2Yh3vn5BC5crgYAjO3jjrcmhsLDkZeJiMi8Nff7m6uEiEycIAgY28cDv8wZibi7u0MhE/DL8UJEL0rCF/vPQ6eziP9zEBHdFgMLkZmwtpLj5bEh2DJrOPr7OqGytgGvbzqOhz/dh1MFFVKXR0TUrhhYiMxMiIca3z8XhX8+0Af2KgVSc0oRs2QP3t12kneCJiKLxcBCZIbkMgHThgYg8cWRuLe3Oxp0IpbtOou73t+JLw5cQL1WJ3WJRERGxUm3RBbgl+MFiP/5BM5fm5TbrYsdXhkXgrF93CEIgsTVERE1rUMm3SYkJEAQBMyZM6fJNsePH8fkyZMREBAAQRCwePHiW7b7+OOPERAQAGtrawwePBiHDh1qS2lEncrYPh5IfHEU/vlAH7jaKZFdXIVnv0zBw5/uR8qFEqnLIyJqs1YHluTkZCxfvhxhYWG3bVddXY1u3bohISEBHh633jvi66+/xosvvog333wTqamp6NevH8aOHQuNRtPa8og6HSu5DNOGBmDXy3fh/+4JgrWVDCkXrmDysv149osUnC2qlLpEIqJWa1VgqaysRGxsLFauXAlnZ+fbto2MjMT777+PRx99FCqV6pZtFi1ahKeffhozZsxA79698emnn8LW1haff/55k+etra1FeXm5wYOIAAdrK7w0pieSXr4bj0b6QiYA244XYMy/d2Pe98eQX3ZV6hKJiFqsVYElLi4OMTExiI6ObnMBdXV1SElJMTiXTCZDdHQ09u/f3+TvxcfHw9HRUf/w9fVtcy1ElsRdbY2EyWHYNmckonu5Q6sTsT45F6Pe34W3f8rElao6qUskImq2FgeW9evXIzU1FfHx8UYpoLi4GFqtFu7u7gbH3d3dUVBQ0OTvzZ8/H2VlZfpHbm6uUeohsjQ93B3w2fSB+P65KAwKdEFdgw4r95zDyPd2Ysn206jiUmgiMgMtCiy5ubmYPXs21q5dC2trabcEV6lUUKvVBg8iatoAf2d8/dchWD0jEn281KiobcCixCyMen8n1uw7j7oGLoUmItPVosCSkpICjUaDiIgIKBQKKBQKJCUlYcmSJVAoFNBqtS0uoEuXLpDL5SgsLDQ4XlhY2OQkXSJqHUEQcFdPN/w4cziWTg1HYBc7FFfW4c3NxzHm30n46Vg+LGSnAyKyMC0KLKNHj0Z6ejrS0tL0j4EDByI2NhZpaWmQy+UtLkCpVGLAgAHYvn27/phOp8P27dsxdOjQFp+PiO5MJhNwfz8v/PrCSLz9YCi62Ktw/nI14r5KxYOf7MOhc1wKTUSmRdGSxg4ODggNDTU4ZmdnB1dXV/3xadOmwdvbWz/Hpa6uDpmZmfo/5+XlIS0tDfb29ggKCgIAvPjii5g+fToGDhyIQYMGYfHixaiqqsKMGTPa/AaJqGlWchliB/tjUn9vfLbnHJbvPou03FL8afl+RPdyx7zxPRHk5iB1mURELQsszZGTkwOZ7MbAzaVLlxAeHq7/eeHChVi4cCFGjRqFXbt2AQCmTJmCoqIivPHGGygoKED//v2xbdu2mybiElH7sFMpMDs6GFMH+2LJ9tNYdygX/ztRiB0nCzEl0g9zooPhrpZ23hoRdW7cmp+IbnJGU4n3tp3Er5mNc8usrWR4PCoQz43qDkdbK4mrIyJL0tzvbwYWImpS8vkSvLv1JA5fuAIAcLBW4NlR3TFjWABslUYfoCWiToiBhYiMQhRF7DylwXvbTuFkQQUAoIu9CrNHB2FKpB+UCt70nYhaj4GFiIxKpxOx+eglLErMQk5J412h/Vxs8X/3BGFSuDes5AwuRNRyDCxE1C7qGnT4OjkHS3acQVFFLQDA28kGz47qhkcG+sLaquXbGxBR58XAQkTtqrquAV/sv4CVe86huLIxuHSxV+HpEYGIHeIPexXnuBDRnTGwEFGHqKnX4pvDuVielI280sY7QTvaWOHxqAA8HhUAZzulxBUSkSljYCGiDlWv1WHjkTws23UW2cVVAABbpRyPDfHHU8MD4cZ9XIjoFhhYiEgSWp2IbRkF+HjnGWTmlwMAlAoZ/jTQB8+M7A5fF1uJKyQiU8LAQkSSEkURu04V4aOdZ5BybR8XuUzAA/298Pxd3bnlPxEBYGCRuhwiukYURRw8V4KPd57BntPFAABBAMb18UDc3UEI9XaUuEIikhIDCxGZnKO5pfh45xn9lv8AEN3LHXOigxlciDopBhYiMllZhRX4aMcZ/HjsEq5/Ao0OccPs6GCE+ThJWhsRdSwGFiIyeWc0lfhox2lsPnoJumufRHf37IrZ0T3Q39dJ0tqIqGMwsBCR2cguqsRHO89g45E8fXAZHtQFT40IxKgeXSEIgrQFElG7YWAhIrNzvrgKH+08gw1H8qC9llyC3ezx5PBATAr35rb/RBaIgYWIzFZuSTVW7zuPr5NzUVnbAABwtVPisSH+eGyIP7o6qCSukIiMhYGFiMxeeU09vknOxarfzuu3/VcqZJjU3wt/GRKAvj5cWURk7hhYiMhiNGh1+OV4IVbuyUZabqn+eD8fR8QO8cf9YV6wUfJyEZE5YmAhIouUcuEK/rv/PLamF6BOqwMAqK0VeHiAL2KH+KF7V3uJKySilmBgISKLVlxZi28PX8RXhy4gt+Sq/nhUd1dMifTFvb3dYatUSFghETUHAwsRdQo6nYik00VYe+ACdpzU6JdF2yrlGNfHA5PCvRHV3RUKuUzaQonolhhYiKjTySu9iq+Tc7HxSB5ySqr1x7s6qDCxnxceDPdGHy8193UhMiEMLETUaYmiiNScUmw8koctxy7hSnW9/rkQDwc8O6o7JoR5ctSFyAQwsBARAahr0GHP6SL8cCQP/8ssRG1D40Rdf1dbPDuqOx6K8IZKwRVGRFJhYCEi+oOyq/X4Yv95fP7beZRU1QEAPNTWeHpkN0wd5MtJukQSYGAhImpCdV0D1h3Kxcrd2SgorwEAuNgp8cSwAPx5sD9c7JQSV0jUeTCwEBHdQW2DFj+k5mHZrrP6SbpKhQwxfT3x2BA/RPg5c4IuUTtjYCEiaqYGrQ4/pedj5Z5sZOSV64+HeDjgsSH+mBTuDXsVLxcRtQcGFiKiVjiaW4ovD1zA5qOX9BN07VUKTAr3wqORflwWTWRkDCxERG1QVl2P71IvYu2BC8gurtIfD3azx6RwbzzQ3ws+zrYSVkhkGRhYiIiMQBRF7D97GWsP5iDxRCHqro26AMCgQBc8GO6N+0I94WhrJWGVROaLgYWIyMjKa+qxLb0AG47k4cC5y7j+6amUy3BPiBsmD/DBXT27woob0hE1GwMLEVE7ulR6FZuPXsLGI3k4WVChP97FXolJ/b3xyEBf9PRwkLBCIvPAwEJE1EFO5Jfjh9SL2HAkD8WVdfrjfb0d8fAAHzzQ3wtOttzbhehWGFiIiDpYvVaH3VlF+PbwRWw/WYh6bePH6/VLRpPCvXF3SFfeCoDodxhYiIgkVFJVh01pefgu5SKOX7qxt4vaWoGYME9M6u+NyAAXyGRcIk2dGwMLEZGJOJFfjo1pedh05JL+VgAA4O1kgwf6e+GhCG8EuXG+C3VODCxERCZGqxNx8NxlbDySh63pBaiobdA/N8DfGY9G+iImzJM3YaROhYGFiMiE1dRrsf2EBhuO5GHnKQ20usaPYgeVAhP7N+6q29fHUeIqidofAwsRkZnQlNfg25SL+Do5V38TRgDo46XGo5G+uK+vJ1ztVRJWSNR+GFiIiMyMTifiQPZlrEvOxS8ZBajTNu6qKxMad9UdH+qJcaEecFdbS1wpkfEwsBARmbGSqjpsOJKHDUcuGtxBGmic7zI+1ANj+3jA14X3MyLzxsBCRGQhckuqsS2jAFsz8pGaU2rwXLifEx4e4IMJYV5wtOH9jMj8MLAQEVmg/LKr+CWjAFszCpB8vgTX5upCpZBhbB8PPDzAB8OCukDO/V3ITDCwEBFZOE1FDTYduYRvU3KRVVipP+7paI2HIrzxUIQPune1l7BCojtjYCEi6iREUUR6Xhm+S7mITWmXUHa1Xv9ciIcDxoV6YHyoJ3q420MQOPJCpoWBhYioE6ptaNzf5dvDudhzuhgNuhsf8d262OnDS6i3muGFTAIDCxFRJ1daXYf/ndBga3o+9pwu1i+TBgAfZxuMD/XAuFBPhPs68Z5GJBkGFiIi0quoqceOkxpsyyjAzlMa1NTfCC8eamuM7eOOcaGeGBTowgm71KEYWIiI6Jaq6xqQdKoIWzMKsOOkBpW/u6dRF3sl7u3tgfGhHhja3RVWcpmElVJnwMBCRER3VFOvxW9nirE1owCJmYUGE3Ydbaxwb293jA/1wPDgLlAp5BJWSpaKgYWIiFqkXqvDgezL2JpRgF+PF6C4sk7/nL1KgXtC3DA+1AN39XSDjZLhhYyDgYWIiFpNqxORfL4E2zIKsC2jAAXlNfrnrK1kGBzoihHBXTCyR1cEu3G5NLVec7+/23RxMiEhAYIgYM6cObdt9+233yIkJATW1tbo27cvfv75Z4PnKysrMXPmTPj4+MDGxga9e/fGp59+2pbSiIioDeQyAUO6ueIfE/tg37x78MPzUfjryG7wcbZBTb0OSVlFWPDTCYz5924Mid+Ov317FJvS8nC5slbq0slCKVr7i8nJyVi+fDnCwsJu227fvn2YOnUq4uPjMWHCBHz11VeYNGkSUlNTERoaCgB48cUXsWPHDnz55ZcICAjAr7/+iueffx5eXl6YOHFia0skIiIjkMkERPg5I8LPGfPHh+BUYQX2ni7G7tPFOJh9GYXltfgu5SK+S7kIAIjwc8L9/bwQ09cTbryzNBlJqy4JVVZWIiIiAp988gkWLFiA/v37Y/HixbdsO2XKFFRVVWHLli36Y0OGDEH//v31oyihoaGYMmUKXn/9dX2bAQMGYPz48ViwYEGzauIlISKijldTr8Xh81ew53QRdp8uxon8G3eWFgRgcKALJoR5YXyoB1ztVRJWSqaqXS8JxcXFISYmBtHR0Xdsu3///pvajR07Fvv379f/HBUVhc2bNyMvLw+iKGLnzp3IysrCmDFjmjxvbW0tysvLDR5ERNSxrK3kGB7cBfPv64Wts0fg4Kuj8eb9vRHh5wRRBA5kl+C1jRkY9M52/OU/B/F1cg40FTV3PjHRH7T4ktD69euRmpqK5OTkZrUvKCiAu7u7wTF3d3cUFBTof166dCn++te/wsfHBwqFAjKZDCtXrsTIkSObPG98fDzeeuutlpZPRETtyF1tjRnDAjFjWCAuXqnGT8fyseVYPtLzyrDndDH2nC4GAPT1dsTdPbvirhA39PNx4mZ1dEctCiy5ubmYPXs2EhMTYW1tvOuSS5cuxYEDB7B582b4+/tj9+7diIuLg5eXV5OjOPPnz8eLL76o/7m8vBy+vr5Gq4mIiNrGx9kWz4zqjmdGdcf54ipsOXYJiZmFOHqxDOl5jY8lO87AxU6JUT264u4QN4wI6gJnO6XUpZMJatEclo0bN+LBBx+EXH5j/b1Wq4UgCJDJZKitrTV4DgD8/Pzw4osvGqwkevPNN7Fx40YcPXoUV69ehaOjIzZs2ICYmBh9m6eeegoXL17Etm3bmlUb57AQEZmHoopaJGUVYecpDXZnFaGi5sZOu4LQOPoyIrgLRgR3RYSfM5QK7rZryZr7/d2iEZbRo0cjPT3d4NiMGTMQEhKCuXPn3hRWAGDo0KHYvn27QWBJTEzE0KFDAQD19fWor6+HTGb4D1Iul0On04GIiCxLVwcVHh7gg4cH+KBeq0PqhSvYcUqDXSeLcKqwAsculuHYxTJ8vPMsbJVyDO3WuOfLXT3dENDFTurySSItCiwODg76pcjX2dnZwdXVVX982rRp8Pb2Rnx8PABg9uzZGDVqFD744APExMRg/fr1OHz4MFasWAEAUKvVGDVqFF5++WXY2NjA398fSUlJ+O9//4tFixYZ4z0SEZGJspLLMLibKwZ3c8X88b1QWF6DvaeLsed0EfaeKUZxZR22n9Rg+0kN8GMmQjwcMD7UE+P7enDDuk6m1fuwNCUnJ8dgtCQqKgpfffUVXnvtNbz66qsIDg7Gxo0bDYLP+vXrMX/+fMTGxqKkpAT+/v54++238eyzzxq7PCIiMmHuamtMHuCDyQN8oNOJOFFQjj2ni7E7qwgHz5XgZEEFThZU4N//y0K3rna4L9QT40I90MdLzfBi4bg1PxERmYXS6jokZhZiW0YB9pwuRp32xrQBH2cbjAjugqHdu2BoN1d0deCeL+aC9xIiIiKLVVFTjx0nNdiaXoBdWRrU1BvOeQx2s0dUd1cM7e6KId1c4WTLlUemioGFiIg6heq6BhzIvoz9Zy9j39nLyMwvx++/2QQB6OnugEGBLogMcMGgQBe485YBJoOBhYiIOqXS6jocyC7B/rPF2Hf2Mk5rKm9q4+9q2xheAlwwuJsL/F25+kgqDCxERERo3Pcl+XwJDp0rQfL5kptGYADAz8VWv/dLVJAr1NZW0hTbCTGwEBER3UJ5TT1SLlxB8rnGEJOWW4oG3Y2vQrlMQH9fJ32ACfNxhJWcm9e1FwYWIiKiZqisbcDB7MuNy6dPFyG7qMrgeRsrOSL8nfSXkML9nGGjvHmjVGodBhYiIqJWuHil+trmdcX47WwxSqvrDZ5XyAT09XHUz3+JDHCBAy8htRoDCxERURvpdCLOFFXi4LkS/SWkgvIagzZymYC+3o6I6u6KqO5dMMCfIzAtwcBCRERkZKIo4uKVqzh0LbwcOHcZFy5XG7RRymXo7+eEqO6uuKunG8K8HSGTcRfepjCwEBERdYC80qvX9oApxv6zl5FfZjgC42qnxKieXXFPiBtGBHeFow0vH/0eAwsREVEHE0URFy5XY3/2Zew5XYQ9WcWoqG3QPy+XCRjg54y7Q9wwpJsL+ng5Qqno3CuQGFiIiIgkVq/V4fD5K9h5SoOdJzU3bWKnVMjQ19sREX5OiPBzRoS/c6fbhZeBhYiIyMTkllRj1ykNkrKKkHLhCq78YQUSAHg72WCAv7N+HxgPR8sOMAwsREREJkwURZy/XI3UC1eQmnMFR3JKcbKgHLo/fCv3cLfHiOCuGB7cBYMDXWCrVEhTcDthYCEiIjIzVbUNOHqxFAfOXsbu08U4drHUIMAo5TIMDHDGiOCuGBHcBb091Wa/AomBhYiIyMyVVtdh39nL2J1VhD2ni5FXetXgeVc7JYZfu3Q0IriLWc5/YWAhIiKyIKIoIru4CruzirD3dDH2Z19GdZ3WoE1PdwdEBbliUIALIgNd0MVeJVG1zcfAQkREZMHqGnRIzbnSuHz6dDHS88puugt1t652GBTggkGBjbcQ8HG2gSCY1iUkBhYiIqJOpKSqDr+dKcahcyVIPl+CkwUVN7XxcrTGkO6uGNrNFVFBXeDtZCNBpYYYWIiIiDqx0uo6HD5/BcnnS3DofAnSL5ah4Q9LkPxcbBHV3RVDr4UYNwnmwDCwEBERkV51XQNSL5Rif3Yx9p29jGMXy6D9Q4AJcrNvHH3p7ooh3VzhbKds97oYWIiIiKhJFTX1OHz+CvadbQwwmfnlBnNgBAHo5aHG0O6NAWZQoAscrI1/HyQGFiIiImq20uo6HMguwf5rAeaPtxGQywRsihuGUG9Ho75uc7+/LWu7PCIiImoVJ1slxoV6YFyoBwBAU1GjDzD7z15GYXkterg7SFYfAwsRERHdxM3BGhP7eWFiPy8AjauQpLyzdOe+pzURERE1i0sHTMC9HQYWIiIiMnkMLERERGTyGFiIiIjI5DGwEBERkcljYCEiIiKTx8BCREREJo+BhYiIiEweAwsRERGZPAYWIiIiMnkMLERERGTyGFiIiIjI5DGwEBERkcljYCEiIiKTp5C6AGMRRREAUF5eLnElRERE1FzXv7evf483xWICS0VFBQDA19dX4kqIiIiopSoqKuDo6Njk84J4p0hjJnQ6HS5dugQHBwcIgmC085aXl8PX1xe5ublQq9VGO6+lYT81D/upedhPzcN+ah72U/NI1U+iKKKiogJeXl6QyZqeqWIxIywymQw+Pj7tdn61Ws1/6M3Afmoe9lPzsJ+ah/3UPOyn5pGin243snIdJ90SERGRyWNgISIiIpPHwHIHKpUKb775JlQqldSlmDT2U/Own5qH/dQ87KfmYT81j6n3k8VMuiUiIiLLxREWIiIiMnkMLERERGTyGFiIiIjI5DGwEBERkcljYLmDjz/+GAEBAbC2tsbgwYNx6NAhqUtqN7t378b9998PLy8vCIKAjRs3GjwviiLeeOMNeHp6wsbGBtHR0Th9+rRBm5KSEsTGxkKtVsPJyQlPPvkkKisrDdocO3YMI0aMgLW1NXx9ffHee++191szqvj4eERGRsLBwQFubm6YNGkSTp06ZdCmpqYGcXFxcHV1hb29PSZPnozCwkKDNjk5OYiJiYGtrS3c3Nzw8ssvo6GhwaDNrl27EBERAZVKhaCgIKxevbq9355RLFu2DGFhYfoNqIYOHYqtW7fqn+/s/dOUhIQECIKAOXPm6I+xr4B//OMfEATB4BESEqJ/nn10Q15eHh577DG4urrCxsYGffv2xeHDh/XPm/XnuEhNWr9+vahUKsXPP/9cPH78uPj000+LTk5OYmFhodSltYuff/5Z/Pvf/y7+8MMPIgBxw4YNBs8nJCSIjo6O4saNG8WjR4+KEydOFAMDA8WrV6/q24wbN07s16+feODAAXHPnj1iUFCQOHXqVP3zZWVloru7uxgbGytmZGSI69atE21sbMTly5d31Ntss7Fjx4qrVq0SMzIyxLS0NPG+++4T/fz8xMrKSn2bZ599VvT19RW3b98uHj58WBwyZIgYFRWlf76hoUEMDQ0Vo6OjxSNHjog///yz2KVLF3H+/Pn6NtnZ2aKtra344osvipmZmeLSpUtFuVwubtu2rUPfb2ts3rxZ/Omnn8SsrCzx1KlT4quvvipaWVmJGRkZoiiyf27l0KFDYkBAgBgWFibOnj1bf5x9JYpvvvmm2KdPHzE/P1//KCoq0j/PPmpUUlIi+vv7i48//rh48OBBMTs7W/zll1/EM2fO6NuY8+c4A8ttDBo0SIyLi9P/rNVqRS8vLzE+Pl7CqjrGHwOLTqcTPTw8xPfff19/rLS0VFSpVOK6detEURTFzMxMEYCYnJysb7N161ZREAQxLy9PFEVR/OSTT0RnZ2extrZW32bu3Lliz5492/kdtR+NRiMCEJOSkkRRbOwXKysr8dtvv9W3OXHihAhA3L9/vyiKjeFQJpOJBQUF+jbLli0T1Wq1vm9eeeUVsU+fPgavNWXKFHHs2LHt/ZbahbOzs/jZZ5+xf26hoqJCDA4OFhMTE8VRo0bpAwv7qtGbb74p9uvX75bPsY9umDt3rjh8+PAmnzf3z3FeEmpCXV0dUlJSEB0drT8mk8kQHR2N/fv3S1iZNM6dO4eCggKD/nB0dMTgwYP1/bF//344OTlh4MCB+jbR0dGQyWQ4ePCgvs3IkSOhVCr1bcaOHYtTp07hypUrHfRujKusrAwA4OLiAgBISUlBfX29QV+FhITAz8/PoK/69u0Ld3d3fZuxY8eivLwcx48f17f5/TmutzG3f39arRbr169HVVUVhg4dyv65hbi4OMTExNz0fthXN5w+fRpeXl7o1q0bYmNjkZOTA4B99HubN2/GwIED8cgjj8DNzQ3h4eFYuXKl/nlz/xxnYGlCcXExtFqtwT9wAHB3d0dBQYFEVUnn+nu+XX8UFBTAzc3N4HmFQgEXFxeDNrc6x+9fw5zodDrMmTMHw4YNQ2hoKIDG96FUKuHk5GTQ9o99dad+aKpNeXk5rl692h5vx6jS09Nhb28PlUqFZ599Fhs2bEDv3r3ZP3+wfv16pKamIj4+/qbn2FeNBg8ejNWrV2Pbtm1YtmwZzp07hxEjRqCiooJ99DvZ2dlYtmwZgoOD8csvv+C5557DrFmzsGbNGgDm/zluMXdrJpJCXFwcMjIysHfvXqlLMTk9e/ZEWloaysrK8N1332H69OlISkqSuiyTkpubi9mzZyMxMRHW1tZSl2Oyxo8fr/9zWFgYBg8eDH9/f3zzzTewsbGRsDLTotPpMHDgQLzzzjsAgPDwcGRkZODTTz/F9OnTJa6u7TjC0oQuXbpALpffNNO8sLAQHh4eElUlnevv+Xb94eHhAY1GY/B8Q0MDSkpKDNrc6hy/fw1zMXPmTGzZsgU7d+6Ej4+P/riHhwfq6upQWlpq0P6PfXWnfmiqjVqtNosPaaVSiaCgIAwYMADx8fHo168fPvzwQ/bP76SkpECj0SAiIgIKhQIKhQJJSUlYsmQJFAoF3N3d2Ve34OTkhB49euDMmTP89/Q7np6e6N27t8GxXr166S+fmfvnOANLE5RKJQYMGIDt27frj+l0Omzfvh1Dhw6VsDJpBAYGwsPDw6A/ysvLcfDgQX1/DB06FKWlpUhJSdG32bFjB3Q6HQYPHqxvs3v3btTX1+vbJCYmomfPnnB2du6gd9M2oihi5syZ2LBhA3bs2IHAwECD5wcMGAArKyuDvjp16hRycnIM+io9Pd3ggyExMRFqtVr/gTN06FCDc1xvY67//nQ6HWpra9k/vzN69Gikp6cjLS1N/xg4cCBiY2P1f2Zf3ayyshJnz56Fp6cn/z39zrBhw27aYiErKwv+/v4ALOBzvF2n9Jq59evXiyqVSly9erWYmZkp/vWvfxWdnJwMZppbkoqKCvHIkSPikSNHRADiokWLxCNHjogXLlwQRbFxOZyTk5O4adMm8dixY+IDDzxwy+Vw4eHh4sGDB8W9e/eKwcHBBsvhSktLRXd3d/Evf/mLmJGRIa5fv160tbU1q2XNzz33nOjo6Cju2rXLYJlldXW1vs2zzz4r+vn5iTt27BAPHz4sDh06VBw6dKj++evLLMeMGSOmpaWJ27ZtE7t27XrLZZYvv/yyeOLECfHjjz82m2WW8+bNE5OSksRz586Jx44dE+fNmycKgiD++uuvoiiyf27n96uERJF9JYqi+NJLL4m7du0Sz507J/72229idHS02KVLF1Gj0YiiyD667tChQ6JCoRDffvtt8fTp0+LatWtFW1tb8csvv9S3MefPcQaWO1i6dKno5+cnKpVKcdCgQeKBAwekLqnd7Ny5UwRw02P69OmiKDYuiXv99ddFd3d3UaVSiaNHjxZPnTplcI7Lly+LU6dOFe3t7UW1Wi3OmDFDrKioMGhz9OhRcfjw4aJKpRK9vb3FhISEjnqLRnGrPgIgrlq1St/m6tWr4vPPPy86OzuLtra24oMPPijm5+cbnOf8+fPi+PHjRRsbG7FLly7iSy+9JNbX1xu02blzp9i/f39RqVSK3bp1M3gNU/bEE0+I/v7+olKpFLt27SqOHj1aH1ZEkf1zO38MLOyrxuXFnp6eolKpFL29vcUpU6YY7C3CPrrhxx9/FENDQ0WVSiWGhISIK1asMHjenD/HBVEUxfYbvyEiIiJqO85hISIiIpPHwEJEREQmj4GFiIiITB4DCxEREZk8BhYiIiIyeQwsREREZPIYWIiIiMjkMbAQERGRyWNgIaJ2t2vXLgiCcNMN6oiImouBhYiM7q677sKcOXP0P0dFRSE/Px+Ojo4d8vpJSUnw9fXtkNcioo6hkLoAIrJ8SqWyXW87/0ebNm3C/fff32GvR0TtjyMsRGRUjz/+OJKSkvDhhx9CEAQIgoDVq1cbXBJavXo1nJycsGXLFvTs2RO2trZ4+OGHUV1djTVr1iAgIADOzs6YNWsWtFqt/ty1tbX429/+Bm9vb9jZ2WHw4MHYtWvXTTVs3rwZEydOBAB899136Nu3L2xsbODq6oro6GhUVVV1RFcQkRFxhIWIjOrDDz9EVlYWQkND8c9//hMAcPz48ZvaVVdXY8mSJVi/fj0qKirw0EMP4cEHH4STkxN+/vlnZGdnY/LkyRg2bBimTJkCAJg5cyYyMzOxfv16eHl5YcOGDRg3bhzS09MRHBysfy2NRoN77rkH+fn5mDp1Kt577z08+OCDqKiowJ49e8B7vhKZHwYWIjIqR0dHKJVK2Nra6i8DnTx58qZ29fX1WLZsGbp37w4AePjhh/HFF1+gsLAQ9vb26N27N+6++27s3LkTU6ZMQU5ODlatWoWcnBx4eXkBAP72t79h27ZtWLVqFd555x0AjZeDxo4dC6VSifz8fDQ0NOChhx6Cv78/AKBv374d0Q1EZGQMLEQkCVtbW31YAQB3d3cEBATA3t7e4JhGowEApKenQ6vVokePHgbnqa2thaurq/7nTZs2YebMmQCAfv36YfTo0ejbty/Gjh2LMWPG4OGHH4azs3N7vjUiagcMLEQkCSsrK4OfBUG45TGdTgcAqKyshFwuR0pKCuRyuUG76yEnPz8fR44cQUxMDABALpcjMTER+/btw6+//oqlS5fi73//Ow4ePIjAwMD2emtE1A446ZaIjE6pVBpMljWG8PBwaLVaaDQaBAUFGTyuX3r68ccfERUVBRcXF/3vCYKAYcOG4a233sKRI0egVCqxYcMGo9ZGRO2PIyxEZHQBAQE4ePAgzp8/D3t7e/0oSVv06NEDsbGxmDZtGj744AOEh4ejqKgI27dvR1hYGGJiYgxWBwHAwYMHsX37dowZMwZubm44ePAgioqK0KtXrzbXQ0QdiyMsRGR0f/vb3yCXy9G7d2907doVOTk5RjnvqlWrMG3aNLz00kvo2bMnJk2ahOTkZPj5+aGqqgrbt283CCxqtRq7d+/Gfffdhx49euC1117DBx98gPHjxxulHiLqOILI9X1EZAF++OEHvPbaa8jMzJS6FCJqBxxhISKLYG9vj3fffVfqMoionXCEhYiIiEweR1iIiIjI5DGwEBERkcljYCEiIiKTx8BCREREJo+BhYiIiEweAwsRERGZPAYWIiIiMnkMLERERGTyGFiIiIjI5P0/mY/wTQj4GigAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cycling_df = pd.read_csv(\n", " \"../data/Tesla_4680/601-828_Capacity_03_MB_CB1_subset.txt\",\n", " sep=\"\\t\",\n", ")\n", "filter_cycling = cycling_df.loc[54811:61000].copy() # Full cycle is [54811:127689]\n", "filter_cycling[\"time/s\"] = filter_cycling[\"time/s\"] - filter_cycling[\"time/s\"].iloc[0]\n", "\n", "# Take every 100th point\n", "filtered_cycling_df = filter_cycling.iloc[\n", " [i for i in range(len(filter_cycling)) if not i % 100 != 0]\n", "]\n", "filtered_cycling_df.plot(x=\"time/s\", y=\"Ecell/V\", kind=\"line\")" ] }, { "cell_type": "markdown", "id": "13", "metadata": {}, "source": [ "### Defining the parameters for inference" ] }, { "cell_type": "code", "execution_count": null, "id": "14", "metadata": {}, "outputs": [], "source": [ "parameters = pybop.Parameters(\n", " pybop.Parameter(\n", " \"R0 [Ohm]\",\n", " prior=pybop.Gaussian(3e-3, 1e-3),\n", " bounds=[1e-5, 1e-2],\n", " ),\n", " pybop.Parameter(\n", " \"R1 [Ohm]\",\n", " prior=pybop.Gaussian(5e-3, 1e-3),\n", " bounds=[1e-5, 1e-2],\n", " ),\n", " pybop.Parameter(\n", " \"R2 [Ohm]\",\n", " prior=pybop.Gaussian(1e-4, 5e-5),\n", " bounds=[1e-5, 1e-2],\n", " ),\n", ")" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "### Form the dataset using the filtered data\n", "In this parameter inference task, we use the filtered time-series data to reduce the computation time. This choice can be validated by comparing the inference results to the full time-series inference, which isn't covered in this example." ] }, { "cell_type": "code", "execution_count": null, "id": "16", "metadata": {}, "outputs": [], "source": [ "dataset_hundred = pybop.Dataset(\n", " {\n", " \"Time [s]\": filtered_cycling_df[\"time/s\"].to_numpy(),\n", " \"Current function [A]\": -filtered_cycling_df[\"I/mA\"].to_numpy()\n", " / 1000, # Convert mA to A\n", " \"Voltage [V]\": filtered_cycling_df[\"Ecell/V\"].to_numpy(),\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "### Construct the model, problem, and non-scaled posterior\n", "\n", "We now construct a two-pair RC model and build the model with an initial SOC based on the first voltage point in the experimental data. \n", "The `FittingProblem` and likelihood classes are constructed with posterior built from the GaussianLogLikelihood." ] }, { "cell_type": "code", "execution_count": null, "id": "18", "metadata": {}, "outputs": [], "source": [ "model = pybop.empirical.Thevenin(\n", " parameter_set=parameter_set, options={\"number of rc elements\": 2}\n", ")\n", "model.build(\n", " initial_state={\n", " \"Initial open-circuit voltage [V]\": dataset_hundred[\"Voltage [V]\"][0]\n", " }\n", ")\n", "\n", "# Generate problem, likelihood, and sampler\n", "problem = pybop.FittingProblem(model, parameters, dataset_hundred)\n", "likelihood = pybop.GaussianLogLikelihood(problem)\n", "posterior = pybop.LogPosterior(likelihood)" ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "### Let's find the Maximum a Posteriori values\n", "Below we identify the parameters using the Maximum a Posterior estimate and the Covariance Matrix Adaptation Evolution Strategy (CMAES). We will use this estimate later in the Bayesian inference task as a starting position for the chains." ] }, { "cell_type": "code", "execution_count": null, "id": "20", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/engs2510/.pyenv/versions/3.12.4/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning:\n", "\n", "os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "OptimisationResult:\n", " Optimised parameters: [4.15690229e-03 1.64236988e-03 8.78939716e-05 2.78474178e-04]\n", " Final cost: 431.8611377082031\n", " Number of iterations: 78\n", " SciPy result available: No\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/engs2510/.pyenv/versions/pybop/lib/python3.12/site-packages/pybamm/solvers/base_solver.py:762: SolverWarning:\n", "\n", "Explicit interpolation times not implemented for CasADi solver with 'safe' mode\n", "\n" ] }, { "data": { "text/html": [ " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "optim = pybop.CMAES(\n", " posterior,\n", " sigma0=[1e-4, 1e-4, 1e-4, 1e-4],\n", " max_iterations=200,\n", " max_unchanged_iterations=40,\n", " parallel=parallel,\n", ")\n", "results_optim = optim.run()\n", "print(results_optim)\n", "pybop.plot.quick(problem, results_optim.x)\n", "pybop.plot.convergence(optim)\n", "pybop.plot.parameters(optim);" ] }, { "cell_type": "markdown", "id": "21", "metadata": {}, "source": [ "### Monte Carlo Sampling\n", "Below we construct the Monte Carlo sampler, specifically the Hamiltionian sampler, which samples from the posterior using Hamiltonion dynamics with gradient information. To minimise the time to execute this notebook, we greatly limit the number of samples; however, this should be increased for actual inference tasks. We initialise the chains from the Maximum a Posteriori point-estimate obtained from the CMAES optimisation above, as this should be close to the mode of the posterior." ] }, { "cell_type": "code", "execution_count": null, "id": "22", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Using Haario-Bardenet adaptive covariance MCMC\n", "Generating 3 chains.\n", "Running in parallel with 12 worker processes.\n", "/Users/engs2510/.pyenv/versions/3.12.4/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning:\n", "\n", "os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n", "\n", "| Iteration: 1 | Iter/s: 0.00 |\n", "| Iteration: 2 | Iter/s: 28.66 |\n", "| Iteration: 3 | Iter/s: 25.26 |\n", "| Iteration: 4 | Iter/s: 74.38 |\n", "| Iteration: 5 | Iter/s: 28.77 |\n", "| Iteration: 6 | Iter/s: 25.46 |\n", "| Iteration: 7 | Iter/s: 52.50 |\n", "| Iteration: 8 | Iter/s: 82.83 |\n", "| Iteration: 9 | Iter/s: 63.31 |\n", "| Iteration: 10 | Iter/s: 1.61 |\n", "| Iteration: 50 | Iter/s: 15.56 |\n", "| Iteration: 100 | Iter/s: 16.85 |\n", "| Iteration: 150 | Iter/s: 44.17 |\n", "| Iteration: 200 | Iter/s: 125.50 |\n", "| Iteration: 250 | Iter/s: 145.31 |\n", "Initial phase completed.\n", "| Iteration: 300 | Iter/s: 18.41 |\n", "| Iteration: 350 | Iter/s: 124.50 |\n", "| Iteration: 400 | Iter/s: 49.71 |\n", "| Iteration: 450 | Iter/s: 15.42 |\n", "| Iteration: 500 | Iter/s: 12.31 |\n", "| Iteration: 550 | Iter/s: 12.83 |\n", "| Iteration: 600 | Iter/s: 16.69 |\n", "| Iteration: 650 | Iter/s: 18.65 |\n", "| Iteration: 700 | Iter/s: 23.45 |\n", "| Iteration: 750 | Iter/s: 9.57 |\n", "| Iteration: 800 | Iter/s: 24.36 |\n", "| Iteration: 850 | Iter/s: 6.11 |\n", "| Iteration: 900 | Iter/s: 7.26 |\n", "| Iteration: 950 | Iter/s: 12.78 |\n", "| Iteration: 1000 | Iter/s: 9.71 |\n", "| Iteration: 1050 | Iter/s: 10.11 |\n", "| Iteration: 1100 | Iter/s: 10.57 |\n", "| Iteration: 1150 | Iter/s: 9.25 |\n", "| Iteration: 1200 | Iter/s: 14.29 |\n", "| Iteration: 1250 | Iter/s: 12.10 |\n", "| Iteration: 1300 | Iter/s: 37.08 |\n", "| Iteration: 1350 | Iter/s: 11.71 |\n", "| Iteration: 1400 | Iter/s: 30.06 |\n", "| Iteration: 1450 | Iter/s: 48.85 |\n", "| Iteration: 1500 | Iter/s: 9.15 |\n", "Halting: Maximum number of iterations (1500) reached.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[1.00493998 1.00643067 1.00368448 1.01600809]\n" ] } ], "source": [ "sampler = pybop.HaarioBardenetACMC(\n", " posterior,\n", " chains=3,\n", " x0=results_optim.x, # Initialise at the MAP estimate\n", " max_iterations=1500, # Increase for accurate posteriors\n", " warm_up=350,\n", " verbose=True,\n", " parallel=parallel, # Only supported for macOS/WSL/Linux\n", ")\n", "\n", "chains = sampler.run()\n", "\n", "# Create summary statistics\n", "posterior_summary = pybop.PosteriorSummary(chains)\n", "\n", "# Create rhat of chains\n", "print(posterior_summary.rhat())" ] }, { "cell_type": "markdown", "id": "23", "metadata": {}, "source": [ "### Plotting\n", "Next, we plot the parameter traces as well as the combined posterior (across all chains) and each chain individually. \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "24", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "posterior_summary.plot_trace()\n", "posterior_summary.plot_posterior()\n", "posterior_summary.plot_chains()" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "As expected, these chains haven't fully converged yet, so the values obtained from the posterior will be biased based on the initial conditions for sampling. Increasing the number of iterations for the sampler while also calibrating the `warm_up` period will provide better results. Next, for completeness, we will plot the identified time-series model against the experimental data." ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/engs2510/.pyenv/versions/pybop/lib/python3.12/site-packages/pybamm/solvers/base_solver.py:762: SolverWarning:\n", "\n", "Explicit interpolation times not implemented for CasADi solver with 'safe' mode\n", "\n" ] }, { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pybop.plot.quick(problem, posterior_summary.mean);" ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, "source": [ "Now, we can compare the identified parameters for each method. As the Bayesian sampler provides us with samples from the posterior, we use the statistical moments when comparing to the point-based result from the optimiser." ] }, { "cell_type": "code", "execution_count": null, "id": "28", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The XNES result is: [4.15690229e-03 1.64236988e-03 8.78939716e-05 2.78474178e-04]\n", "The posterior means are:[4.14717807e-03 1.64491239e-03 9.20455193e-05 2.86617819e-04]\n", "The difference between the two methods is: [9.72421663e-06 2.54251113e-06 4.15154776e-06 8.14364117e-06]\n" ] } ], "source": [ "print(f\"The XNES result is: {results_optim.x}\")\n", "print(f\"The posterior means are:{posterior_summary.mean}\")\n", "print(\n", " f\"The difference between the two methods is: {np.abs(results_optim.x - posterior_summary.mean)}\"\n", ")" ] }, { "cell_type": "markdown", "id": "29", "metadata": {}, "source": [ "The sampled posterior also gives us information of the uncertainty in the inference process. This is present in the standard deviations as well as the lower and upper confidence intervals, shown below." ] }, { "cell_type": "code", "execution_count": null, "id": "30", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "posterior_summary.summary_table()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 5 }