Skip to content

Instantly share code, notes, and snippets.

@aphearin

aphearin/jax_odeint_demo.ipynb

Created November 30, 2021 16:39
Show Gist options
  • Save aphearin/795323e7c2e1339551c7bb6cf58981d2 to your computer and use it in GitHub Desktop.
Save aphearin/795323e7c2e1339551c7bb6cf58981d2 to your computer and use it in GitHub Desktop.
Download ZIP
Demo of using JAX to solve an ODE and differentiate through the solution
Display the source blob
Display the rendered blob
Raw
jax_odeint_demo.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "demonstrated-singles",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "executive-subscriber",
"metadata": {},
"outputs": [],
"source": [
"from matplotlib import pyplot as plt\n",
"from matplotlib import lines as mlines"
]
},
{
"cell_type": "markdown",
"id": "basic-tucson",
"metadata": {},
"source": [
"## Define the ODE\n",
"\n",
"Let's use jax and scipy to solve the following ordinary differential equation:\n",
"$$dy/dt = -3t^2e^y$$\n",
"\n",
"### Analytical solution\n",
"Using symbolic differentiation, you can convince yourself that the following analytical solution is correct:\n",
"$$y(t) = -\\log(t^3 - A),$$ where $A$ is some constant whose value is determined by the boundary condition.\n",
"\n",
"### Boundary value\n",
"Let's use $t_0=0$ to define the boundary condition of our ODE, so that our boundary value problem becomes:\n",
"$$y_0 = -\\log(-A),$$\n",
"or equivalently:\n",
"$$A = -e^{-y_0}$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "characteristic-obligation",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"def dy_dt_numpy(y, t):\n",
" return -3*t*t*np.exp(y)\n",
"\n",
"def get_a_from_y0(y0):\n",
" return -np.exp(-y0)\n",
"\n",
"def analytical_ode_result(y0, t):\n",
" A = get_a_from_y0(y0)\n",
" return -np.log(t**3 - A)"
]
},
{
"cell_type": "markdown",
"id": "underlying-praise",
"metadata": {},
"source": [
"# Scipy solver\n",
"\n",
"Now let's use scipy to numerically integrate this ODE and compare the result to the analytical solution."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "waiting-trash",
"metadata": {},
"outputs": [],
"source": [
"from scipy.integrate import odeint as spodeint\n",
"\n",
"t0 = 0.0\n",
"t_domain = np.linspace(t0, 15, 100)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "joint-parks",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEhCAYAAABhpec9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8uklEQVR4nO3deVxU1fvA8c9hc00HXDJ3wdTUyljazCUdcinXWCqXbBG0srRMs0WtTAOzNK0E26yskFFb7JvFYGa2CpSaZaljbqipOO4mwvn9wTA/REBQmDsDz/v14lVz5y7PnXHuc89yz1Faa4QQQojieBkdgBBCCPcmiUIIIUSJJFEIIYQokSQKIYQQJZJEIYQQokQ+RgcghLtQSk0A7ECWY1GA1jqxmPUmOdZLcCyuB5gAu9Z6YjHbxDr2n1TgrSAgFAjWWqvyOI9Cxw0G4gC01uHlvX9RNSjpHivciVIqDjBprWNdeMxA8i6mI7XW9vMtd7yXDNgKJwWlVAQQW9RFubhtCrw3UWttu/gzOmffwUBcWROFUspU+Lwdy13+HQljSYlCuJuk869S7lKA8MIXRa21TSmVACQDRV1kDxZeoLW2KKUClFIJZbyQziCvRFIR7Be4XRRwTokKY74jYSBpoxBuRWudobXOcNXxHFVCGcXdyWutrUCAo6RQKo7qqihHieR8xzc5tskAzru+ixVZAnH1dySMJ4lCVHWxwNrzrJPmWK8s0oDSJBdzgf+3lvEYFSa/esnoOIR7kKonUWGUUjFA/p26iQKNw473zmo0dtyBJ0Bew2uBhlg7/99oHAgEFaznd+xromO9yPzSgVIqBQggr42huDvgQOB8d8fp5FXDlIUNCCtphUKfAUW1BxSxfqk/z2L2EQwsALLy2ywKfH5xju8hwrH/QEeJCyBRa20v/B0Vcy4F4yr4Hc4g7/swAWFFtdUIN6W1lj/5K/c/IAYILPDaRN6FCPIuHDEF3gvMfw0EAykF3jMDW8lrPC24LKWI4yUUWhZRcLsiYgwENHk9jko6l4i8n8pZy5KBCSVsEwekF7FNOjAh/30gohw+z+Qi3iv4GQYWen3WZ+xYllDoOzlnneLeK8XxzY5zDSy0jdnof6fyV7o/qXoSFclZXaPz7paTHHXyMfrsO94I8rqJwrkNr1nk9RRyLtd57QaBSilzgWWJFHHXr0u4S9f/3y5xvraBwCLiKo2i2j2sWut4nXc33bOM+yvq8wwm7wJsK/SezXGXX1r2C1m3lMfPIi9hF/w8bLhfm4wohiQKUSEcF+5ApZRWSqUopWJ0XvWPmUIX0AIXzrLIIO/OtiBrfqOzIyHZS7Gf0lywgshrcyiLQIpOFE75F9T81wUTXxHrFvd5hhZznK1ASBljvhClPb7dBbGICiKJQlQIRx/8SMCfvGqWSEdX04o0g/+/6zY7Sh7nk0AxvXsKMPP/bSSlZaYU3Uj12W0npuLWK+HzLHYb8toDykUJPbhKe/ysYtcSbk8ShagokyDvrllrbdV5DZ/5DcfnXHTyu4mWQTCFegk5LrqBpemWWmCbeMc2hUsn+XFFkPe0taW0+3RsY9Vl70JaUtzFfZ7WYrYLIu/5kKLYOTeJmM4TW5GfzwUeX3gYSRSiwhRRR25z1FMvLuK9YqtdgNCCicRxIS6uL3+C468sF+lwIK5wgnEkj1iKb0uoV3hBgW1GluH4BYcPKWmdoj7PDPLaA4ILrGcCQnUxPZ9wtBkUWhZYaFnhKrkiYyvD8Ysq3RSOQbgp6R4rKspB8i4g+QnAxP+PORSrlJpQoLunSec90Zw/ZEaoow4+/0KT5liWv59ARzVMURLJe8q61ENhONYNd8RUcLtAXfRQHBNwtLUU6D5acH8lbZOllCr4RLdzrCdKrgIr6fOMdMSef2EPxJHcCn2mExztQXalVELBz5+8u/9YpZRNa20ptE5WfomqqO/oPMcPJq80FJh/fEeijwDsSqmMUlYRCgPJWE/CrV3IOEVKqYiyVBUJIUomVU+iUlFKmSVJCFG+JFEId3fenjuOKpL8OnJTxYYjRNUjiUK4rQINw6FFtQUUkIxj4D4pTQhR/qSNQgghRImkRCGEEKJEla57bP369XXLli2NDkMIITxKenr6Aa11g6Leq3SJomXLlqSllXVYHiGEqNqUUtuLe0+qnoQQQpRIEoUQQogSSaIQQghRIkkUQgghSiSJQgghRIkkUQghhCiRJAohhBAlqnTPUVyoHYd3MHL2SP5Z8w+dYzpTo2YNTvxzAvtmO/2H9OeaZtfQsWFHfLzkIxNCVC1y1XM4+t9R1v6xlsNrD3Oo2yFy/XLJWpGF/kbzSZ1PwBt8fvah5paajJ83ntE3jaZ+zfpGhy2EEBWu0g0KGBoaqsvryewzOWfYuW8np3xP8dve33jr7bdYl7aOA+YD+Hn7EfxvMA/f9jB33nJnuRxPCCGMopRK11qHFvmeJIqy2/jvRhLWJvDa8NfIbZjLuNnjeO7m56jtV7tCjyuEEBWlpEQhjdkXoEPDDrx666vYNtoY/uRwXvnpFdq82Ibp70w3OjQhhCh3kiguQotGLVg4bCE/3PsDp9ec5qmYp5j+pSQLIUTlIomiHNzQ7Ab+SPqDGybcwFO/PMXYFWPJyc0xOiwhhCgXkijKSUP/hnz3/HeMvW4scz6aQ6uurTh1+pTRYQlRKdlsNiIjI7FarUaHUiVIoihH3l7evNL7FQbWG8jOLTu5J+kecnWu0WEJUalYrVZsNhs2m83oUKoMeY6iAix7eRlP936aF358gUtXXMrLvV7Gy0tyshDlwWw2AxAQEGBwJFWHJIoK8nz489hP25nzzBw2Lt1IyoIUo0MSLjJ2xVh+2/ubIcfu1KgTs3vPLtM2FosFm82G2WwmODgYgPDwcFJS5N+syCOJooIopZjdZzZfzPwC6yYry/5cxqArBhkdlhBnsVqtmM1mZsyYQVpaGsHBwdjt9jJPJzxx4sTzrhMdHe1MRMKzSKKoQD7ePmxM2Ui3hd24+5O7uaLBFbSr387osEQFK+sdvZECAwMxmUxYrVZiY2MBSEtLc1bvWCwWAgMDSUtLIyYmptj9xMXFuSReYQxJFBWspl9NlkYt5eppVxPWI4w/U/+kaYOmRoclBJCXKGw2G3a7ncDAQABSUlIIDw/HarWSlZVFREQEAPHx8UyYMKHcY0hMTCQ9Pb3EdSZOnOiMT7ieJAoXaFa3GVNCpjA2YSzD3h7GygkrUUoZHZYQAGRkZDhLEICzdJGQkEB4eDgAJpOJlJSUYhPFxVQ9lVRSEe5BEoWLPHLHI/zX6D8mfjuRN9Le4IGwB4wOSQgAsrKyCAoKAsBut2Oz2QgMDMRutzvXCQgIICsrq9h9SNVT5SZ9Nl3o8W6P0zuoN49Me4Rl3y4zOhwhAIiKimLr1q1YLBZmzJhBaGjeuHAmk8mZLLKystymO2pGRgbx8fGkpaURFxdHYmKi0SFVelKicCGlFLO7zab96Pbcf/h+eqf0poZvDaPDEoKEhAQA1q5d62zUDg8Pdz7UZrPZnNVQRgsODiY4OLhC2ktE0aRE4WJtm7XlzU/fJKtbFo+nPG50OKKKs9ls9OzZE8irdsrIyHA2XpvNZux2O1arlYyMDLkwV2EyH4VBxq0Yx+xVs3mj8xuM6jfK6HBEFWaxWIC8pBETE4PJZDI2IGEImbjIDZ3MPkn9K+vz3/7/2Ld9H/Vq1zM6JCFEFSYTF7mhGr41eGP2G+RG5PJ4qlRBCSHclzRmG2h47+Fs8t3EjDUz6NuyLxFXRxgdkhBCnENKFAab0m0Kjf9qTHT3aLZlbjM6HCGEOIdbJwqllEkpNcHxl6yUqnSPcFbzqUb83fHoVpqnv3na6HCEEOIc7l71NElr7RwbQCm1VSmF1rpSPWEzpNcQ/vD5g+lrpjN8y3B6te5ldEhCCOHktiUKpZQJKDwKWAJw/kFlPNAz3Z4hyCuI2++8ncwDmUaHI4QQTm6bKBzMSqmCycLOucmjUqjuU50JV07g+IbjPPCmjAMlhHAfbpsotNZ2rbW/1rrgxLjhQKWdTT1mYAyjF43m0/8+ZfX21UaHI4QQgBsnisIcVVFmzlP1lJmZiVLK+Td16lRXhFduZvabSStTK+6aeRdZR4ofrVMIIVzF3RuzC1oARGqtM0paqXHjxmRmem4dfy2/WjzZ9klGjh1JRHYEKxNWGh2SEG4lPj6egwcPEh0dTVZWFsnJyc5BDUXFcGmiUEpFANHnWS1Lax1baLsJQILWutJWOxV0f+/7+fiJj/nG9xt+3fMr11x2jdEhCeFWEhMTSUxMxGw2s2DBAqPDqfRcmii01hbAUpZtHMklIz9JKKXMVSFhJE9Opv3r7RmxdAQ/3f8TNarJcORCQN48GYcOHTI6jCrFrauelFJmIACwOtooAoBgKnGDdj7/Gv7Ed4lneL/hDFo/iBUJK4wOSZRW+lg49Jsxx/bvBCGzy7SJxWLBZrNhNpudU5WGh4eTkpJS/vGVo4yMDEwmk8yl7QJumygciSH/X2rBCsgylUg82bBrh/Fs8LOknkjlrwN/0bZ+W6NDEpWM1WrFbDYzY8YM0tLSCA4Oxm63U9YRmC9mzuwLYbFYMJvNWK1WEhISZCrWCibDjLu5vcf20v619rRv0J7V96zGS3lMRzXhAfLnxw4JCSE5OZnAwEDnxTc5ORmLxUJgYCBpaWnExLjnCDpBQUEkJCRgNpuNDsWjlTTMuNuWKESeRrUb8UqvVxjxzAiGbh3Khy98aHRIohIJDAzEZrNht9udVTgpKSmEh4djtVrJyspyzngXHx9fIbPcJSYmkp6eXuI6EydOdMaXkZFxVskkODiYlJQUSRQVSBKFBxh21TDG7xjP4q2LeWH8C7Tyb2V0SKISycjIOOsia7VaiY2NJSEhwTlPtslkIiUlpdhEcTFVT2UpqWRkZNCzZ8+zGrPtdjtBQUGl3ocoO0kUHsDLy4tVy1dx3fvXEbM8hq+Hfo1SyuiwRCWRlZXlvNDa7XZndZTdbneuExAQQFZW8Q+AuqqNIDg4+Jxj2Ww2oqKiXHL8qkoShYfo0KwDs26ZxSjLKB5PfJyXYl8yOiRRSURFRTFx4kQsFgtr164lNDSvmtpkMjmTRVZWFgEBAQZG+f9CQ0OJj4/HZDKxdetWkpOTZZ7vCiaJwoPEhMQw+aHJzHp1Fvfddh9XNLnC6JBEJZH/ZPPatWuJjc173jU8PBybLW+oNZvN5qyGMlpwcHC59Z4SpSNdaDyIUork+clUv6c6474ZR2XrsSZcz2az0bNnTyCv2ikjI8PZeG02m7Hb7VitVjIyMiqkIVt4Buke64Hm/TKPMV+OYe7Nc3mo60NGhyM8nMWS92iSzWYjJiZGqnGqqJK6x0qi8EC5OpcOMR3YtGgT3//8PTdeeaPRIQkhPFxJiUKqnjyQl/Ii4cEEfDr6MGnNJHJ1rtEhCSEqMUkUHqprp67Me30eq/9dzfy0+UaHI4SoxCRReLCYkBi6+nfl4eEPk5qWanQ4QohKShKFB1NK8VL4S+TuyeXBhQ9yJveM0SEJISohSRQeLuyKMN5JfYe/6v9F3BoZQVMIUf4kUVQCw4OHE90hmsmJk1n01SKjwxFCVDKSKCoBpRSzesyC5TD6ydGczD5pdEhCiEpEEkUl0SSgCQuSF3C071EeT3nc6HBEFZeRkUFISEipRpV1hYqKx93Os6JIoqhE7g2/l0dvepTXfnyNuZ/NNTocUYUFBwc7x4wqi8TExHOWRUZGFrncFfFU1H4r6jwriiSKSmZ6z+n4p/jzyF2PsDlzs9HhiCqsrKPN5k+gVFhsbGy5TEpUUaPfutt5VgRJFJVMNZ9qfPDKB/gO8uWB1AfkqW3hMYqb08JsNjtnt6sMPPE8ZZjxSqjvDX2Z5zePmOUxTEuZxuRbJhsdUpUyduxYfvvtN0OO3alTJ2bPnl2mbTIyMoC8O92UlBTi4uIwmUxYrVYmTpxIdHT0WdOkFp6WtKhtC7NYLM7t8uePiIyMxGazsWDBArKysrDZbNhsNkwmE4GBgZjNZjIyMhg5ciRms9l5gbXZbCQkJBAWFobJZCIgIMA57Hhp4ylOYmKic9KmlJQU5/Drdrvd+V5WVpYzvsLyP7PQ0FASEhLIyMhwnndCQgJWq7XU51ncMUvzvZQ3SRSV1P3B9/PB5x8wZcAUGn3UiJiBpZ9uUlQtkZGRJCcnO4cXnzhxIgkJCZjNZmJjY0lOTiYlJQXIm8AoISHBeTErbtvC8t9PSkpyXrhjY2MJDQ11vs6/yBecGjW/DWDr1q3OZeHh4aSnp2MymZwX2Pw5t0sbT1Hi4+Mxm83OpJM/FwdAz549z5rXOzIy8qwElS//M8tfNz/+pKQk5/ulPc/ijlma76W8SaKopJRSvDfqPTp+0ZFn1z1LRK8IAmq4xwxllV1Z7+iNln/RBQgMDDzrAgmcdUceEBDgvDiVZtuCIiIiGDlyJHa73Tl7Xmnu9gMCApwXUIvFgslkcm4XHBxMaur/D19TlngKCwwMZOTIkcTGxhIVFeW8kFsslnPu1KOjo5kxYwbJycml3v/5FD7P8x2zpO+lvEkbRSXWolELvvn0G/Z77WfEJyPIzZX2ClG0/KlQ09LSzpkb+3yNtSVtW1hUVBSJiYnY7fYLqiax2WznxFM42ZQlnoIiIiKYNGkSycnJ+Pv7O7u8rl27tshjliUJlVVpjunKqWklUVRyoY1Dmd51Op/P+JzB4wYbHY5wM3a7nZCQECZNmkRERIRzvuz898p72/yqIKvVWuJ0pvmTKRUWHBxc7MX/Ys4F8toXIiIiSElJQWtNWloaNpuNoKCgc45ZlkRXUrIq7jwv9pjlTRJFFfBol0dpXKMxn//5OT/s/MHocIQbSUtLO6sqp+Ac2ee7Yz7ftllZWedcoAMDAzGZTEVePAMDAzl48OA5ywvuJ78BuWBs+c8eXEg8BaWkpJy13/xjxcTEONsV8iUlJTFp0qRz4ss/j4LnV7hKqDTneb5jupokiirAy8uLDd9soNVtrYhMjuTf4/8aHZJwE2azmdDQUBITE513+aGhoc473eTkZKxWq3Pe7PyePBaLpUzbFpTfBlBYREQENpuNxMTEsxq4C+8nNTWVhIQELBaLM46ynktRgoKCnO9bLBbCwsLO6VVksViIj48nNjaW4ODgIuMzm80EBAQ44wsPD8dqtToTWmnPs7THLPy9VASZCrUKWbd3Hdc+dy0BGwL4J/UfqvlVMzokUQVZLBZnryThPmQqVAHA1Y2uZnij4ezdsJdxyeOMDkdUIbGxsc5nCEpqmxDuSRJFFbNg8gLuSbyHN7a8wZI/lhgdjqgiIiMjsdvtZGRkuO3Tx6J4UvVUBf135j+6v9uddEs6SZOSGNRtkNEhCSEMJlVP4izVfKqxwLyAnB9zuPe5ezl08pDRIQkh3JgkiiqqY8uOLLUu5Vi3Y9y55E6Zb1sIUSxJFFXYgLABvHHbG3y17ivCHwo3OhwhhJuSRFHF3R98P9ftvY5VC1bx4vIXjQ5HCOGGJFEIVr21is4vdOaZX5/h23++NTocIYSbkUQhqO5XnS8e/oLWAa25beptpPxScaNQCiE8jyQKAUDd6nVJ6pfE8SXHGfzAYA6eOHcsGiFE1SSJQjhd1fwq3lv2Hqf7nmZg0kBOnTlldEhCCDcgiUKcZah5KO9Fvcca2xq6PNCFMznSbVaIqk4ShThHdMdooogibUEaQ18ZanQ4wg1FRkY6R0MVlZ9MhSqK9NGzH5Hjn0PS0SRu+OkGHrn+EaNDEm4kNjZWxmyqQqREIYrk5eVF0tgkBl8xmLHvjeXxVx83OiThRsxmsySKKkQShSiWt5c3Hwz6AP8f/Hlp6kt8uelLo0PyGN27d+fdd98FIDs7m+7du/PBBx8AcOLECbp3705SUhIAhw8fpnv37ixduhSAAwcO0L17dz7//HMA9u7dS/fu3VmxYgUAO3fupHv37s4Jbmw2G927d+fbby/sGZj8iX4sFguxsbHO5TabzTlxTv4EOZA3yU5ISIhzTmmr1UpISAixsbHOyXpiY2Ox2WxYLBaCgoIIDw93zt4WGRlJSEjIOTO4CfflUVVPSqkErXXs+dcU5aWGbw3S/5dOn7f7ELk0kpV3r+TaJtcaHZYoJ/Hx8ZjNZuccEQWnAg0PDyc9PR2TyURGRgYjR44kPT2d4OBgYmNj2bp1K5BXuoiOjmbt2rXOCYnMZjMhISHOdZKSkpyzucXGxhIaGup8LTyA1toj/oA4IOV864WEhGhR/jKPZOqWr7TUNW+rqT/97lOjwxHlJDk5WQcHB+uEhAR96NAhfejQobOWF5T/Xv77EyZMcL6Oi4s767XWWpvNZp2QkKC11tpkMp21b+F+gDRdzHXVI6qelFIyJZbBLrvkMpJvS+bUt6e466m72HZom9EhiXIQERHBpEmTSE5Oxt/f31mdZLPZCAgIOGvdspYAAgMDnSWKqKgoEhMTsdvt0rbhgTwiUQChgIwrYbDQ1qH875v/4dvblx7v9WDXkV1GhyQuktVqJSIigpSUFLTWpKWlOacrzcrKuqh922w2goKCAJg4cSIJCQlYrVaZCtUDuX2iUEpFAIuNjkPk6RXSi6+Hf82BAwdo37M9v9t+NzokcRFSUlLOapcwm81n/bfgewWfm8jKynI2TucruK7dbsdmsxETEwPklS5MJtNFJx9hDLduzFZKmQC71tqulCrVNpmZmRRcd8qUKUydOrVC4quqwpqEMStsFrGzYuk/pz9rZ6ylXs16RoclLkBQUBBWq5XAwEDsdjthYWHOqqHU1FRmzJhBWFgY8P/JIyMjg+TkZGw2G1ar1bkcwGKxALB27VpSUs6uBIiNjSUqKsoVpyXKmVvPma2UitFaJzr+3wxM1FqXOMOOzJntOp+u+5Toz6O5osEVpAxNoX6t+kaHJAwSHx/PwYMHiYuLK3Ydi8Xi7BUl3E9Jc2a7tEThqEaKPs9qWVrrWEcDttUFYYkLNODqAXxa61P6Pd+PlvEtWbd6HUFNgowOS7iR2NhYIiMjCQwMlLYJD+bSRKG1tgCWUq4eAJgLVCOFAYFKqQmARWttK3ZL4TK9WvdicrfJTP5mMoOTBvPN6G8IqBFw/g1FpWG1WklKSsJutxMeHn5WVVRkZCR2u52MjAwpTXgwt656KkgpFQPEaq1DSlpPqp6M8cXfXzB48WCuqHcFSwYsIaixlCyE8CQlVT25fa8ncCaJSBwlCkcjt3Ajt7a5lc/u+IzfF/5Oh7AO2PZJgU+IysIjEoXWOlFrHa619tdax2ut7UbHJM7Vq3Uv4sbEkds+l1stt5J5NNPokIQQ5cAjEoXwHI8NeYzUN1PZdWQX18Vfx7e/XdhAdUII9yGJQpS7Li26YB1mZc87ezD3NfP7PnkoTwhPJolCVIjrml7H4g8WUyeyDt0XdictUzoYCOGpJFGICjO4+2B+nvIztf1qc9NDNzHz/ZlGhySEuACSKESFah3QmlXDVqHWKia+OBHLH6V9jEYI4S7KlCiUUvcrpQYrpepUVECi8mlZryV//PwH1z10HVHJUbz6w6tGhySEKIOylij8gScBu1Jqs1LqDUkcojRaNW5F6v2p9A3syyNDH6Hz0M7k6lyjwxJClEKZEoXWeqbjyT1/YBRwmLzEcciROMZL0hDFqelbE0u0hfbt2vPDiR8Yvmw4p3NOGx2WEOI8LqiNQmt9WGudqrV+wpE4woBUoDWwTSk1qDyDFJVHdb/q/P7F77ww7gUWbVjEjc/dyLZMmS1PCHdW5kRRVIlBa51B3nyro4Ag4E6lVI9yiE9UQkopnuzyJG/2fpP0Welc1fcqbIdkyA8h3FVZG7PnA/8opQ462icKJoMgAK21XWsdBciYwqJE9113H6++/Srefby5/s3r+XHnj0aHJIQoQllLFCla6wAgHFCARSmVo5TKAQ4CKKU6OdaV+gRxXmOixvDL+F+oW70uN915E6OnjzY6JCFEIWVNFHal1GCtdYbWepQjaQRorb211i851rEopd4gr8FbiPNqU68Nq4asos7BOsz/Yj5Ppj4pPaKEcCNlmrhIa52qlKqrlOqhtV7pWHa40DqtlVJ1Cy8XoiRNApqw69ddPLLiEWasmUHaH2m8c+c7NKnfxOjQhKjyytyY7ejxtPJ861x4SKKqqlW9FgsGLOCV8FdImZZCmxvbsPngZqPDEqLKkyE8hFtRSjH2xrHMnDUT7x7eXPvmtXy15SujwxKiSpNEIdzS+KHjWTdjHc3rNqf3uN7cEnMLZ3LOGB2WEFWSJArhtlr5t+KHe3+g1fFWpHyTQsTiCI78d8TosISociRRCLdWy68WW6xbiFsQx/LNywmeE8zSVUuNDkuIKkVprY2OoVyFhobqtLQLmCTn9GFsacn8sXkHvboF4+vrw9Z/9rBp6y56dw/B29uHLf/sYcs/e+h1cxjKyxvb9r3szDxAt84hoHzYvS+LA4eOcfWVHcHLF/vRU/yXrbi0cTPwrgFefqBU+Z90FbF6+2p6R/Xm5PqTvP7164zuIs9cCFFelFLpjiGZzlGm7rGV2vFtfPbGSMZ9AIcSwVQLLJ/DEx/DiXeghh98uBSmLIGc90F5wduL4cXP4cz7ebuY/SG8lpK3PsCU92Dhd2BfkPd6zEL4LEOxPaEBeNfiecsJMmzZLHvhevCty6KV+9lzSDM+5laoVp/ftx0j18fEVaFdoPql4OVrzGfjJrq26MpPH/7EXXPu4oGVD5BxKIPZvWZTq1oto0MTolKTEkW+MyfYt/k7du7eR6cr2+Dj48OevQfYlbmPkKvb4uWl2JW5j5279nJ9aAcUmm3/7GL7rj10v/FKyD3Dxk1b2bY9k9t6XA252fyQtpG/bbsYMeAayDnJsq9/Y+Pm3Tw94ko4c4yX3v+V9VsO8t5jTeG0neGzdrLun9Osm5EXUv9ZsOMA/OZ4PSyhGqd1dZKeuwlqNsXy/RFq+Telz60DoXYQVG9YJUosZ3LPMPmbycz4YAbVrdVZ/ulyegb3NDosITxaSSUKSRRuRufmoM4cgVMH2PDrzxyz7+aGDv5wcg8vzv8fp0/YmRxVC07spNOjB2hWDz4fn7dtt2ledAwM4LUnekCdtnz162latb+BNsG9wLu6sSdWAWa8N4PJkyfjO9yX+bfPZ/jVw40OSQiPJYmikjp++F+O/buJS2seg2NbeSrufZqbThLb/RT62FbqxWiir4c37vOCSy4nJvE0/Xp1pt/AOyEgGGo0MvoULtruI7sZsnQI3277lk5bOrF8znJ5mluICyCJogrSZ07xx9oVVDu9g9b+BziW+StX37OCh8xnGNcHjp2CDhO9eXF0J+6MHkyu/7Ucr96BSwIuMzr0MjuTe4aRc0fy7rh3aXh3Qz6f9jnXNrnW6LCE8CiSKIST/s+OOryePZu+4fEX3uX+bmfo3nIX63fANU/C0qdaMmDAbRyrGcqRalfSOMhzRotftHIRk9ZNYs+xPTzU6iHiouPw8/UzOiwhPIIkClGy04fYsW45b731JiNvzqWp+hXL98eJfBXSXm5OSOfe7PcOgYZdaNDsCqOjLdGhk4cY8dEIPnvoMy674TJ+XPIjLUwtjA5LCLcniUKUTW42/6z7kqUfv8mY8Bx8D61h2uIjTFkCB5M6Ymrdm399r6Nu4C1Uq+l+U6Tn5uYy6oVRLDq0CJ+GPswJn8Pwa4bj5SXPlwpRHEkU4uLk5rDxh2S+t35MTJcjsH8Nw+Zl891fim2L+6Ka3sZJ/57UqH+50ZGexXbIxvBlw/k+4XuaqCakp6Rz6SWXGh2WEG5JEoUoX9nHSF06h51/rGRE8DY4vo1rn4F2LU289/JD0HQABIS4xTMdObk59B/TnxXrV9BgQAMS+yXSv21/o8MSwu1IohAVR2v04T95efoELvPewl1XbebMmVzCpvgy7u5uDB/1JDToAl7GDgKwft96hi8bzrrf1xH0ZxAp76fQqnErQ2MSwp2UlCik0lZcHKVQpvY8Fr+cu2ZsgkH7ONTuVVo2qY/p2GpI7cG/7zbk/v5t+fu7tyA325Awr7r0Kn4Z+Qv9avZj6y9buWHBDSz/e7khsQjhaSRRiPJVvT4Nwsaw7LtM+s/Igi5LWH/kGpJS/ubU9/fD0kZssUSRmjyTnOzTLg3Nz9uPz+I/Y3XGaho0bEC/D/sRdl8Ym3fKLHpClEQShag4PrWg2WDMY1PZf+AQVw5ZBo37kPjeJ/S9awJHP24KaQ9j35qCzs11WVhdLu9Cekw6o1uMJu3dNILHBLN442IqWzWsEOVF2iiEy508mkXG12/QueGvsHs5A2b+x6FT1Vj94XhoNQzqtHVZLMtWL+P5P57n132/0tW3Ky8NeomwK8Jcdnwh3IU0Zgv3lX2ERfMe5+Su77g/9C90bi5D3wzg9sGDGDzyRahev8JDOJN7hld+fIWJt09E+SheXfoqo0JH4e3lXeHHFsJdSKIQnuFEJlnr3qbH0BcY1e0Uo8J9+K9Bb/63rT23Dn0Gvxq1K/TwK9NXMmn5JH7hF8IuDWNS+0kM6jqoQo8phLuQRCE8Tu7B3/Da/gHLkt5icLydr56uwy0D7yGnxXC861fc+FNaaxZtWMSop0dx/Ivj3JNwD3OHz6WWn0yOJCo3SRTCY505fYrU5JmYm67He89nzPzsNB/9UoPvkp+nVvt7oFpAhRx3887NDH9uOD81+YlmdZsxudNk7u12rwwDIioteY5CeCwfv+r0GvIM3t2SYVAmzcNGENK6JrX+HA/LGvPuU9exauks0OXba+ryZpfz44IfWXPvGmpn12Zk35G0i2jH1qyt5XocITyBlCiEZ8r6ldwtb3F5/ze4PiiXRRNaQdB9HPIfiH+TDuV6qBOnTnDH43eQ6p1KTv0cxl4zlgldJhBQp2JKM0IYQUoUovIJuAava+exccsBZs1+DWq3Ytc3T9OoZUcWPtEJdn0OuWfK5VA1q9fks7mfsfn5zdze/nbino3j0qBL+ejXj+TZC1ElSKIQHq16bX8ahT0APVPx7fsjjw2/gS7Nd8Pq/mS83Jgpo7pxaNeGcjlW40sas2jwIl559BXqd6nPXZ/dxS0f3ML3m74vl/0L4a7cPlEopUxKqTilVIzjz3OmXBMudWnQ9Ux/6wcCR2VCl2Ws3laPVxauxnvFVbDqVvakv8uZ06cu+jhj7xjLTstO5vaZy8+//cxNV95E7yd6c/DEwXI4CyHcj1u3USilTECy1jrc8ToGCNdaRxa3jbRRiIIOZ26g7v7FYHuL8Ml7OJHtx/eLn4Cg+6FWs4ve/987/yZyTCQb2m7AFGBi/NXjGdt9LDWr1yyH6IVwHY/tHquUSgaStNYWx2sTEKC1thW3jSQKUaTcM3z6zjOc2v4l0e3Wo4EHkpoz7N6HuLHfOLjIp7A37NvAuK/GkTotleqnq5O0Iol+bfuh3GBODiFKw5MThQb8gQDApLXOON82kijEeR3bxvbVL3HdkPnE35HL8F7NOdVsBIf8B3FZYKcL3q3WmimJU1iwZgF7W++lR6sejAkaw8DOA8stdCEqikf2eirQFhFaYFmyo1QhxIWr3YoWfV9jR+Zh7njsQ7jkcpIWPEfzNtfw+8JbYG/qBT2XoZTiudjn2PHuDub2mUvamjQG3TSIm5+6mR2Hd1TAiQjhGm6bKIDA/P/RWtscpYkkYEFJG2VmZqKUcv5NnTq1gsMUnsqvRm38Wt8JPa10GW3l2VE30aFmGqw0s2DMZbw8sT/61IEy79fX25eHrn2IDc9voOuIrvzo+yNt5rbhnsR72LpbHtgTnsdtq56UUmYgBfDXWtsdy4KBdK11sRW/UvUkLkrOKdhh4c57x3HgwAFSnqoGzaPY5teflsGDURcwhMeOwzuY/M1kFo5eiLe3N9M+msbD1z9MTV9p8Bbuw23aKJRSEUD0eVbL0lrHKqUCga0Fk0L+Mgokj8IkUYjycnLPWmrsWsjhjQtpMuoYYwc2ZNrzU6HlEPCtU+b9LVm1hJmpM/nZ52ca1WhEn1N9mDdhnvSQEm7BbRJFWSmlDgGtCpUoUrXW/sVtI4lClLeTR/bzwdzHuM5/LVeZNrEtqyZxK1vw1LMv0+zK3mXe33fbvyPmpRg2zdtEo5GNmDVmFnd0vAMv5c41waKy88jGbIcZQFSB19GOZUK4TI06DRj51HtcNfoP6PUL6Ueu48MVf+L1bR/46nr2/jyHE0dK35bRpUUXNs7ZyPSF02lwTQOGLB1Cy5iWTJo3iVwXTgkrRGm5dYkCQCk1oeBrrXV8SetLiUK4wsnDmdTYkwxb5jNkxiZ+2OzFluUP4t12NNS9otT7ydW5JG9M5p7B93Dy9EmCJwXzfI/n6dO6jzyDIVzKY6ueLoQkCuFSWvPdZ/OwrV3E3R0zIDebsZYm3NLvLvoOfx68q5VqN6dOn2L+t/N59Y9X2bZ3G3WX1+W5Z59jzKAxkjCES0iiEMIVTv3L4XVvEDxwGqN7nGH84PrkthzBNp/eBF3ds1S7yM7J5tmPn+XFMS+SE5nD9ddez1M3PEXfdn1l0iRRoTy5jUIIz1G9IXWvm8LmXScZ88Ln0LArXye/TOtOZla+FAw7LJCbXeIufL19mTZkGkf2HGH+yPnsPrKbfiP7UbdNXT774zMZ1lwYQkoUQlSgvbbfWDh3AuNu3IRf9k4+WluHDUc6MHXmO/gFtD3v9qdzThM7LZalq5dypOsRQi4LYUiDIYzpPwYfbx8XnIGoKqREIYRBGgV2YuIrX+N3+zbo9gW/7m3AipU/4vtlO1jZiz9S55D934lit/fz9uOdKe9w4OsDvNX/Lfbv2c+jEY/SZGATPtrwETm5OS48G1FVSYlCCBc7bbfht/N9sv9aQPORuzFfVZ3354yD1vdD7cAStz11+hSPznqUlDMpbMndQvP/mtPNqxvzJs2jTq2yPwQoRD5pzBbCDeVkn+bLD6fR4MQqrjN9z8EjuQx9qx7TJo8npM+j4O1X7La5OpdPNn3Cg489yF7rXho904jxPccTExLDJdUuceFZiMpCqp6EcEPevn7cdvdzXDd6NQzYztZLYvjzHzu+v0+CT5qw68tYNq39X5HbeikvBl8xmN3Ld/POl+9wRfMrGJ8ynnrX1yN8dDj/Hv/XxWcjKjMpUQjhRnLPZOP170rYmsijM5bxWopm3+IbMF09CppHgE/x40J9v+17ou6IIvOSTKrfXJ3hVw0nukk0PYJ7uPAMhKeSEoUQHsLLxxca94IuS5j4+nqSXh6OyecA/HQ394abeOa+UMgqev6uzq06s/vn3fz50Z8Mu2oY7yx9h54hPbnpqZv4fsf30rVWXDApUQjh7rRG/7ua+++/j8a+//B8RA74d2LRxqvoO2wq/pe2KnKzDbYNPPT8Q6wPWo89x067w+3o07QPLz78In4+xbd/iKpJGrOFqCxOH4J/PmTjynl0HL2JuSN8eOj+SHJa3o1q1BOvIp6tOH76OO/+9i4TRk7gROYJmk9qztgbxnJPp3sw1TC5/hyEW5KqJyEqCz9/aPMgHUb9SfrKDxg6bARkfslnL/Um8LLqbP3fI3B8+1mb1PKrxYPXPsihtEMkfJxAC/8WPPrlo9RrVY8eD/Vga5bMuidKJiUKITxdzim+TZ7OGwmJfHDvPny8FUv/7siJup0Z8tBLKN9a52ySujGVkQ+OZHvT7ejLNb2b9eZm35t5bMhjMqZUFSVVT0JUFce3g+1dbrs3jv2HTvLzjDrQ4g52+Pam2dUDzpnKNfNoJm+sfYPZc2dz7JNjBD4RyOO3P87Qq4ZS26+2QSchjCCJQogqJjfnDPv//JRLj3zKf7ZkGo8+xbAe/syePgFaDYOaTc5a337MzpQFU/iu5nf8uvdXqv1QjY6XdOTD+R/Spn4bg85CuJIkCiGqsP+O7WfRvAm0r5nB9fXXs+cQ3LuwHjOmPkqnW8ae9WyG1pofd/3IiJgRbNm1BR2h6RXUi141e/Fgvwfx85XeUpWVJAohRJ6jW/h+6QsMefx9vpqQQ9vml7BZ3cIev27cdNvos3pNZR7J5M1f3+T1Va+z77l91O1ZlwnPTOC+a+7j0tqXGngSoiJIohBCnEXn5qD2r4Zt7/Hwsx+wIPUM+xY2p06HuznVKIrqDTs61z1x6gRT5k9hTfYafjrxEz77fbjs18uIj4snunO0zMBXSUiiEEIU6/jhf8lYMZcuDX6BfVb6vJhLHVMASQueg+bRUL2+c92/DvzFY7Mf43+v/g89WtO2WVtuv/R2Rlw/gsubXW7gWYiLJYlCCFE6J3bz8rOjqHEsg9GdM9F4M/rjptwxZATdI54A7+oAHDt1jCWblpCQnsCP8T/Cv3Dnm3cSExJDtxbdpJThgSRRCCHK7tA6dv/8OtcOeZPnbs/lvvA6nGgwkJ8OXEW3AWPwdjRsL1m1hHe+eYc1NdZw+NRhan1Si34D+jH78dnSluFBJFEIIS5YTvZpcvasxG/3x3z8cRJ3vnqKb5+vR9e+I8huEolPgzCUlxcnsk/w7k/vMun+SRxpcwSfa33o06oP15y4hkkjJlHdr7rRpyJKIIlCCFEuThw5yJcfTmdQ27/x2vcV05Zk89HPfqQtmUiNdsPhktYA/Ln/T97+9W0SP0jkyPtHCIgJIDYqlhFXj5DnMtyUJAohRPn77yBLFjzF6tTPmROZCcD0r5tQP/BGYsbPhpqNOXHqBM+/9Tzr665nhW0FuWtyqbO9DnHvxjEkeIjMxudGJFEIISrW8Z3o7R/T487naGk6xjuxChp2Y+mmtnQdOI76TduSeTSTsS+O5auUrzhy6xFq+tak095ORNwYwSPRj+ClZIwpI0miEEK4zOkDG/Dbs5TdGe/RdISNaVGKp0b1IqdpFEfr9qBug+b8vPtn3k5/mzfveRPdVNPsnmYMu2oYvRr2ouuVXY0+hSpJEoUQwuV0bi7r1iTT6PQaGh1fzqq1/9ArDqwv30iX2x6Apv3JOpnN4ozFfLr7U7769Sv0y5qWd7Vk/MPjuaPjHdSrWc/o06gyJFEIIYylNVvWLmH+3Bd5tm8mtfQe3v/eh2UbGrDw9Re4pG0kf+zZwYT4CWxuuJm/9d947/amQUYDnp7+NPfefC81fGsYfRaVmkxcJIQwllK0vjaCl95Po9YduyB8DcfrdOHgwYPUXncvLG3Arwvu5p4Ol/PXpHR+i/2NWy+7lf3b9/PQqodoNKsRt8bdysz3Z5J9Jtvos6lypEQhhDCOzoX938OOxXQemkCdatl8+WR1uKwP32Z24Moe95F+dDOLNixi0ZOLOLPnDJc9dRl3XnknvRr1wnylWSZaKidS9SSEcHs52afJ2vwVDU6kcPTvxTQYsY+YHl68+kRfaB7BdjqRvD6FNafW8MXfX3BmzhkuaX0J414cx51X3km7+u2MPgWPJolCCOFRcnPO8OOXC2hw+mfaeKWyacsurnwCFk++mkF3PcgBU2cemxPHbyd+Y0OdDejTmlof1yLqgSieufcZWvm3MvoUPI4kCiGE59KaneuXM++VaTzcbR9Nqm3ni98U8Svq8MGcsfi1v5W5qz9lzhNzONb1GLSEq/2upsWeFkx7dBpXtrzS6DPwCJIohBCVg9Zg38AnC1/gpQXLWTnxBH4+8NH6VthOBDH8sfF8tPM35s1/jZ2LdsI46NyxMz1MPejfvj+h7Yq8DgokUQghKqsjm2HXMkZPmMUvG/8l/QXgkjYs/as9JwLasD2oNkl/JLMhcQNsguvnXE/0VdHcFngbrRu2Njp6tyKJQghR6Z3K2kr1g1+jdyylRbSVG1pD0oRG0HQAb/+cw8osO+sb/sWGfzfA29CgaQMmzpzI7e1vp6WppdHhG04ShRCiSrH/ux3735/Qku85sOl/XBpznKcH+/LsIwPYfUkIo+d+wW/Z29kZtBNyoc6SOvS/qz/PxD5Dm3pVc3RbSRRCiCrr9MmjpC55hTY1/yRIfcva3/dww1RYPu1Kwnr1491th5jy+HscDzkOHaFtjbY0/L0hkx6eRO+Q3lVmtr6SEoWPq4MRQghX8qtxCX2GTs57oXNp2PJzntgxk7CWB6m3aTqXroFm2X58Mnwov19Sj+c/tfLdB9/xXY3vaPl9S3qYetCpRidiB8Xi55jVr6qRRCGEqDqUFy2uGcC0BQPyXh/bRv0T8XTcuJTLj35M26Nn+CerBu3NTbk3ug8L9m7n3YR3yf0ul6lTpjKg0wC6+HdhQKcBBNQJMPZcXEiqnoQQAuD0Ydj7Nc8+N411Gzax9JHToLyZ+mVjDvvWxrt3IAu2rubIe0dR/yr6ze3HwHYDuaXlLTTxb2J09BdN2iiEEKIscnPg4M/oXZ/Tbdhc2jY4zoKRoGu1YsyiHI6aqvP9FcfZemQ3zIfLOl7GuKnj6N+2P23rtzU6+gsio8cKIURZeHlDgxtR18xg9e/HeH3ZFgh7HbvX5bz96Q46HP2bzY2z2Bd8I7d3akCTppoJ1gm0m9OO6s2r0+/pfqz6ZxXZOZVjpFu3L1EopSYAdsdLk9Y6vqT1pUQhhKhIJ49mkZ35DXWOruK7r5fQ9ak9fDIO+nZvTcp+f6Yv+J2M1v9x8vJc6pysg2mNiQfHP8h9t9zn1hMxeWzVk1JqQsHEoJQKBswlJQtJFEIIV9G5ufyVvoLmPhuoeSiV+YtWMvrtHP6YVQNTmw4sWGvnzTe3sHMweDXwouPpjjTa34gp46dwQ+sb3KrrrSd3j40GnElBa52hlJpkYDxCCOGkvLxoF9YX6AtMZOg1e2nZfT7t2v6L2rMC+9pt7N8PBzs35a/qjZm55C++XLKezvW/pln9Zlx98mqu8r+KSfdOona12kafTrHcvUSRAti11pGO1zFAltbaUtw2UqIQQrgFrbGtT2Xddx8xqONe2PcNfWecZGcWrHj1CqwnvZg59w82/qXxG+dHtxbdaG1vTb9r+9H7Wtc/6OfJVU+BQLrj5QzAVlKSAEkUQgg3lXOKb5bNxf7PGgZdsRUObyRoHFzbxpfHxl7O4qzDvPPsbg7UgpYPtKRP6z60ONKCu3vdTaOARhUenscmCnA2ZocDZiBeaz2xpPUbN26s9+zZ43w9ZcoUpk6dWqExCiFEWeUe28nC16bQ2GcrvVpu5Pjhg9QbBRMia3PNrZfy1o6dWKeeJjvUi24ju9ErqBdB2UEM7ja4QqZ/9dg2CqVUApCgtY5XSpmBZKVUYH5VVFEaN25MZmam64IUQogL4FW7GfdMfDvvhc7lzD+rmTx6OuFtDhLmtYE2XtmE5sK8602cOvMnCz/5hnWzoE50HQbeMZCbm9xMWP0wOrTqUOGxujRRKKUiyGugLkmW1jrW0cPJrrXOANBaW5VSrYBtFR2nEEK4lPKibqvuPDmne97r7GPUWJvE6MjXuK3jcZr5/U2LHBhZG17tVI203RaeXfYe/7wH7Z5sx8CbB9L50s70aNODmtVqlnt4Lk0UjvaFEtsYCggADhba3q6UspZ7YEII4U58axN44328fON9ea9P7qFJjQQGbf2QqKDjDMk+QcNTkNBA8XzLLJZtiOfN13K5ZkA4X4/4utzDcdsns7XWVvLaJpyUUibAZkhAQghhlBqXcd2gqbz52d/4RuyCW/+kQ/cYBnRtRYTpPxZdmkvEUUiod6ZCDu/WbRRArFIqDtiav+B8jdlCCFGpKQV129F/dAL9R5M3LtWhDPofeZ1WHbpVzCHdvddTWUn3WCGEKDsZFFAIIcQFk0QhhBCiRJIohBBClEgShRBCiBJJoihEhvtwP/KduB/5TtxTRX0v0uupEKUUle0z8XTynbgf+U7c08V8L9LrSQghxAWTRCGEEKJEla7qSSm1H9h+EbtoDMjws+5FvhP3I9+Je7qY76WF1rpBUW9UukQhhBCifEnVkxBCiBJJohBCCFEiSRRCCCFK5O7DjLuMY25uG3kTJqG1TjQ2oqpNKRUDhADJjkWRQJzWWuYjcSHHHDAxQL2ihviX343rlfSdVNTvRhIF4JjzYq1jBj6UUnFKqYj818IwUeT9IDKAkZIkXMsxT70JCCrmffnduNj5vhOHcv/dSKLIE1MoMycBcZR+2lZRAbTW/kbHUJU5ZplEKRVG3sWpMPnduFgpvpMK+d1U+TYKpVRwEYvtgNnFoQjhMeR3U7VIiSKvbjWr0LLCr4UBHPWtWUj9tzuS342bqojfjSSKYopvkNdopLW2uy4UUUAaYM+vX1VKJSulsqT+222YintDfjeGqpDfTZWveiKvuBxQaFnh18LFtNYZhRrh1gKTjIpHnMOO/G7cTkX9biRR5BXRTIWWmQDkrsg4jt4dBdmAourFhTHkd+OGKup3U+UThdY6g7y7o4ICAKvroxEASqlAIMXRX7wg6R7rJuR3434q8ndT5ROFw2KlVESB1+FAglHBVHWOovPEQnem0eR1vRTuQ343bqQifzcyeqyD4wnTDCAQpIeN0Rx3R/kXoXrAVvlOXMvRBdYMxDoWJQBWR2kifx353bjQ+b6TivrdSKIQQghRIql6EkIIUSJJFEIIIUokiUIIIUSJJFEIIYQokSQKIYQQJZJEIYQQokSSKIQQQpRIEoUQLuR4IEoIjyKJQgjXSi407IUQbk+ezBbChZRSGvCXEVaFJ5EShRAu4hgC2iZJQngaKVEIUcEcCSKcvMHa7OQNxb1WZusTnkIShRAuopRKB2ZIghCeRhKFEC4i7RPCU0kbhRAuIO0TwpNJohDCNcKRaUKFh5JEIYRrmIF0AKWUyVHCEMIjSKIQwjVM/H+JIkZrLaUL4TGkMVsIF3A8jR0I2Mib49hubERClJ4kCiGEECWSqichhBAlkkQhhBCiRJIohBBClEgShRBCiBJJohBCCFEiSRRCCCFKJIlCCCFEiSRRCCGEKJEkCiGEECX6P6hB094wij92AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"y0 = -1.0\n",
"scipy_result = spodeint(dy_dt_numpy, y0, t_domain)\n",
"correct_result = analytical_ode_result(y0, t_domain)\n",
"__=ax.plot(t_domain, correct_result, color='green')\n",
"__=ax.plot(t_domain, scipy_result, ':', color='k')\n",
"\n",
"y0 = -5.0\n",
"scipy_result = spodeint(dy_dt_numpy, y0, t_domain)\n",
"correct_result = analytical_ode_result(y0, t_domain)\n",
"__=ax.plot(t_domain, correct_result, color='orange')\n",
"__=ax.plot(t_domain, scipy_result, ':', color='k')\n",
"\n",
"solid_line=mlines.Line2D([],[],ls='-',c='k',label=r'${\\rm analytical\\ solution}$')\n",
"scipy_line=mlines.Line2D([],[],ls=':',c='k',label=r'${\\rm scipy}$')\n",
"green_line=mlines.Line2D([],[],ls='-',c='green',label=r'$y_0 = -1$')\n",
"orange_line=mlines.Line2D([],[],ls='-',c='orange',label=r'$y_0 = -5$')\n",
"leg=ax.legend(handles=[green_line, orange_line, solid_line, scipy_line])\n",
"\n",
"xlabel = ax.set_xlabel(r'$t$')\n",
"ylabel = ax.set_ylabel(r'$y$')\n",
"title = ax.set_title(r'${\\rm scipy\\ ODE\\ solution}$')"
]
},
{
"cell_type": "markdown",
"id": "appropriate-customs",
"metadata": {},
"source": [
"# jax solver\n",
"\n",
"Now let's use jax instead of scipy."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "biblical-princess",
"metadata": {},
"outputs": [
{
"ename": "TracerArrayConversionError",
"evalue": "The numpy.ndarray conversion method __array__() was called on the JAX Tracer object Traced<ShapedArray(float32[], weak_type=True)>with<DynamicJaxprTrace(level=1/0)> (https://jax.readthedocs.io/en/latest/errors.html#jax._src.errors.TracerArrayConversionError)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTracerArrayConversionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-7-4a3937a0a2ef>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mjax\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexperimental\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mode\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0modeint\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mjodeint\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mjax_result\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjodeint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdy_dt_numpy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_domain\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/experimental/ode.py\u001b[0m in \u001b[0;36modeint\u001b[0;34m(func, y0, t, rtol, atol, mxstep, *args)\u001b[0m\n\u001b[1;32m 170\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 171\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 172\u001b[0;31m \u001b[0mconverted\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconsts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcustom_derivatives\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclosure_convert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 173\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0m_odeint_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconverted\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0matol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmxstep\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mconsts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/custom_derivatives.py\u001b[0m in \u001b[0;36mclosure_convert\u001b[0;34m(fun, *example_args)\u001b[0m\n\u001b[1;32m 917\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0m_closure_convert_for_avals\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__wrapped__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_tree\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_avals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 918\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 919\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0m_closure_convert_for_avals\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_tree\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_avals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 920\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 921\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mcache\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/_src/util.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 196\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 197\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 198\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mcached\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbool\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx64_enabled\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 199\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 200\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcache_clear\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcached\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcache_clear\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/_src/util.py\u001b[0m in \u001b[0;36mcached\u001b[0;34m(_, *args, **kwargs)\u001b[0m\n\u001b[1;32m 189\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mfunctools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlru_cache\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax_size\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 190\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mcached\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 191\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 192\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 193\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mfunctools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwraps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/custom_derivatives.py\u001b[0m in \u001b[0;36m_closure_convert_for_avals\u001b[0;34m(fun, in_tree, in_avals)\u001b[0m\n\u001b[1;32m 923\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0momnistaging_enabled\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 924\u001b[0m \u001b[0mwrapped_fun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_tree\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mflatten_fun_nokwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwrap_init\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_tree\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 925\u001b[0;31m \u001b[0mjaxpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_pvals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconsts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpe\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrace_to_jaxpr_dynamic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapped_fun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_avals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 926\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 927\u001b[0m \u001b[0min_pvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mpe\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mPartialVal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munknown\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maval\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0maval\u001b[0m \u001b[0;32min\u001b[0m \u001b[0min_avals\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/interpreters/partial_eval.py\u001b[0m in \u001b[0;36mtrace_to_jaxpr_dynamic\u001b[0;34m(fun, in_avals)\u001b[0m\n\u001b[1;32m 1188\u001b[0m \u001b[0mmain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msource_info\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun_sourceinfo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# type: ignore\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1189\u001b[0m \u001b[0mmain\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjaxpr_stack\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# type: ignore\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1190\u001b[0;31m \u001b[0mjaxpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_avals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconsts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtrace_to_subjaxpr_dynamic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_avals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1191\u001b[0m \u001b[0;32mdel\u001b[0m \u001b[0mmain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1192\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mjaxpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_avals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconsts\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/interpreters/partial_eval.py\u001b[0m in \u001b[0;36mtrace_to_subjaxpr_dynamic\u001b[0;34m(fun, main, in_avals)\u001b[0m\n\u001b[1;32m 1198\u001b[0m \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mDynamicJaxprTrace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcur_sublevel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1199\u001b[0m \u001b[0min_tracers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnew_arg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0min_avals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1200\u001b[0;31m \u001b[0mans\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall_wrapped\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0min_tracers\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1201\u001b[0m \u001b[0mout_tracers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfull_raise\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mans\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1202\u001b[0m \u001b[0mjaxpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_avals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconsts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mframe\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_jaxpr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0min_tracers\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout_tracers\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/linear_util.py\u001b[0m in \u001b[0;36mcall_wrapped\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 164\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 165\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 166\u001b[0;31m \u001b[0mans\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 167\u001b[0m \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 168\u001b[0m \u001b[0;31m# Some transformations yield from inside context managers, so we have to\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-3-85fc1b6d88ac>\u001b[0m in \u001b[0;36mdy_dt_numpy\u001b[0;34m(y, t)\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mdy_dt_numpy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mget_a_from_y0\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/core.py\u001b[0m in \u001b[0;36m__array__\u001b[0;34m(self, *args, **kw)\u001b[0m\n\u001b[1;32m 486\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 487\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__array__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 488\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mTracerArrayConversionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 489\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 490\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__index__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mTracerArrayConversionError\u001b[0m: The numpy.ndarray conversion method __array__() was called on the JAX Tracer object Traced<ShapedArray(float32[], weak_type=True)>with<DynamicJaxprTrace(level=1/0)> (https://jax.readthedocs.io/en/latest/errors.html#jax._src.errors.TracerArrayConversionError)"
]
}
],
"source": [
"from jax.experimental.ode import odeint as jodeint\n",
"jax_result = jodeint(dy_dt_numpy, y0, t_domain)"
]
},
{
"cell_type": "markdown",
"id": "potential-action",
"metadata": {},
"source": [
"# What happened?\n",
"\n",
"In order to use the implementation of `odeint` in jax, the ODE must be implemented in jax. So let's try again after implementing $dy/dt$ in jax instead of numpy."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "senior-ethics",
"metadata": {},
"outputs": [],
"source": [
"from jax import numpy as jnp\n",
"from jax import jit as jjit\n",
"\n",
"@jjit\n",
"def dy_dt_jax(y, t):\n",
" return -3*t*t*jnp.exp(y)\n",
"\n",
"@jjit\n",
"def get_a_from_y0_jax(y0):\n",
" return -jnp.exp(-y0)\n",
"\n",
"@jjit\n",
"def analytical_ode_result_jax(y0, t):\n",
" A = get_a_from_y0_jax(y0)\n",
" return -jnp.log(t**3 - A)\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "weighted-sugar",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEhCAYAAABhpec9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7+0lEQVR4nO3deVxU1f/H8ddhc00HXHMXzF0zwKzMJR1yKddAK5dsEbSytEyyRS1NE7M0rQRatLRCRm2xbxaDmmmbOG5pljLuuOO4myz39wfD/BBhBAPuDHyej4ePRzNz79zPZZr7nnPOvfcoTdMQQggh8uOhdwFCCCFcmwSFEEIIpyQohBBCOCVBIYQQwikJCiGEEE5JUAghhHDKS+8CRNmglDICIUA4YAVGappmyWdZg6ZptnzeIyzHe4RommbN8Xq0/TUTMCO/989nmxMAG5Bqf8pP07SYfJabaF8u2v50NcAA2DRNi8xnnQj7+8fleCkACAYCNU1TBa21IJRSgcBMAE3TQoryvUXZo+Q6ClGSlFIJQHxeB2H766FAhLODmz0QgjVNC8r1fKD9+TzfO5/38ifrgDoyZzjl97z9tXjAmjsUnNWe3zo5XovMGXpFITssbiQo8gprpdRMwKBpWkQRlSjchHQ9CVfjDxiVUob8FrAfqPyVUuG5XjIWJiTsEsg6SNtybcNKVoshPp/1TuVRlwmItwdZYcwgq0VS1Gz/Yd1BeTwXx/+3okQZIkEhXIY9HCz2f3kdqHIKA6KzA8X+a75QIWHvErLk90te0zQz4Gd/7wKxB9Uge4vkets32NexkBWQruSaVoimaZbCdOeJ0kOCQrgSo/3gPIOsPv182ZczAbHZB+W8xjWuIwLYeJ1lkq5XSz7rFCRcjDn+21zIbRSb7C4mvesQrkMGs4XL0TTNpJSKV0r5X6fffiSwl6yuoxvpN/cnq/XizCau37rJzQq0d7aAvdsse+DcacjZl83+OxjIMdCe633yHIDP8T6BQCyQmj1uYV8/kqyxjBh768lAVtfeBPuqMYAf9m6nnGMezrafY0DdRlb4+9nfu31eYzXCdUmLQrgEe6sg50E7huu3Kmz25QY5G9Nwsj3IcbDORyqF/3VtI++uJKNSaoL9F3uBgs1+IDZrmma2t6LMZJ0tlT0IbtY0zWQfH1lqP1kgT/Zuo5G5noshR2vG/j7RZA28R9n/2eyBnXvw3un27dubaf9b2Oz7YCIrhHK2poSLk6AQriJ36yH7VNd82cNhI1ldPbGF2ViObV1vbMCfGxsUzqslZLYfeCOB7oV4L0eo2MMxzv5r/aq/mf01ax6D/NdjK+xyhdh+KllnSuX8e1hxvTEZ4YR0PQlXEWI/+FxFKZU9bpGXcE3TopRSFiD5OsvmpSAHrACygqgw/Mk7KBw0TbMppXJeA5Jn7fbuoHillEbWL/94+3M5u6NySgaC8ni+qAUXYvu2Yq9GFCsJCqE7e0DMyOO8/Wpk/Zq+5gBq77owQVbrQCkVSdaprL6F2HQ0WWf3RDlZxkiuLpcCMFKAFkOuM4gMeS1jv54hzN56CgYilVJBZB2Q8+NX8FKvL5+xIkMhtn+97j3h4qTrSRQrpVRyrv7ovA5i/vkM5saR/9lDubs9osjq9ijwef72dfzzasmA45Rbm71fvUDs65hv4DTS/Fo2E+21Zvfxh9iXNeezTgBZ14bkx8a1n4HhOrXl9fe50e0LNyRBIYqbBXvXjf1XsT+wNPtF+3N5/gK2H2ytua9jUEpF53N2TxgQnt+BPx8hwMzc1z3Y3yOC/FsG1XI/kWOdkdcunr8ctw/J7/XcYw7WHH+bwBzLGbj+lel5Dc7753oud5fcNbUVcvt5fb65axAuTLqeRHGLJOtsHxtZB+Xu2a0He0tjJmBQSiXl/hVuP0D68f8HchtZB+JApVR8Hn362S2XRKXUDCDmetdW2FslIfazkXJ2r/jncyuOCfbtWHOcPprz/Zytk6qUynlFt+NeT+RxgZvdKfu2svfNwP/fwynMXnf2Qd0fe7DluAVJsFJqgr31lD02Ep1jjMNAVgsgQilltZ/BlHOZVPvpyjnfL1zTtBhn27fXEEhWi8g/uwZ76IcCNqWUpZBjSkIncq8nIYQQTknXkxBCCKckKIQQQjglQSGEEMIpCQohhBBOlbqznqpXr641atRI7zKEEMKtbNq06aSmaTXyeq3UBUWjRo1ISirsHReEEKJsU0rtz+816XoSQgjhlASFEEIIpyQohBBCOCVBIYQQwikJCiGEEE5JUAghhHBKgkIIIYRTpe46iht14MwBRs4Zyb71++gY3pEKFStwcd9FbLtt9B3Sl9vq30brmq3x8pA/mRCibJGjnt25f8+xcedGzmw8w+kup8n0ySR1VSraGo2vqnwFnuD1uxcV91Rk/PzxjL57NNUrVte7bCGEKHalbj6K4OBgraiuzE7PSOfgsYNc9r7MlqNb+Ojjj9iatJWTxpP4ePoQeDyQZ+5/hofufahItieEEHpRSm3SNC04z9ckKApvx/EdRG+M5r3h75FZM5Nxc8bx+j2vU9mncrFuVwghiouzoJDB7BvQqmYr3r3vXaw7rAx/aTjv/PYOTd9syvRPputdmhBCFDkJiv+gYe2GLBq2iF8e+4Ur66/wcvjLTP9ewkIIUbpIUBSBO+vfyc64ndw54U5e/uNlxq4aS0Zmht5lCSFEkZCgKCI1fWvy89SfGdthLHO/mEvjzo25fOWy3mUJUSpZrVbCwsIwm816l1ImSFAUIU8PT97p+Q79q/Xn4J6DPBr3KJlapt5lCVGqmM1mrFYrVqtV71LKDLmOohiseHsFr/R8hTd+fYNaq2rxdo+38fCQTBaiKBiNRgD8/Px0rqTskKAoJlNDpmK7YmPuq3PZsXwHCbEJepckSsjYVWPZcnSLLttuV7sdc3rOKdQ6JpMJq9WK0WgkMDAQgJCQEBIS5P9ZkUWCopgopZjTaw7fzfoO8y4zK/5awYAWA/QuS4irmM1mjEYjM2bMICkpicDAQGw2W6GnE46MjLzuMoMHD3YEkXAvEhTFyMvTix0JO+iyqAuPfPUILWq0oHn15nqXJYpZYX/R68nf3x+DwYDZbCYiIgKApKQkR/eOyWTC39+fpKQkwsPD832fmTNnlki9Qh8SFMWsok9Flg9azq3TbqV9t/b8lfgX9WrU07ssIYCsoLBardhsNvz9/QFISEggJCQEs9lMamoqoaGhAERFRTFhwoQiryEmJoZNmzY5XSYyMtJRnyh5EhQloH7V+kwOmszY6LEM+3gYqyesRimld1lCAGCxWBwtCMDRuoiOjiYkJAQAg8FAQkJCvkHxX7qenLVUhGuQoCghzz74LP/W/pfInyL5IOkDnmz/pN4lCQFAamoqAQEBANhsNqxWK/7+/thsNscyfn5+pKam5vse0vVUusk5myXohS4v0DOgJ89Oe5YVP63QuxwhABg0aBDJycmYTCZmzJhBcHDWfeEMBoMjLFJTU13mdFSLxUJUVBRJSUnMnDmTmJgYvUsq9aRFUYKUUszpMoeWo1vyxJkn6JnQkwreFfQuSwiio6MB2Lhxo2NQOyQkxHFRm9VqdXRD6S0wMJDAwMBiGS8ReZMWRQlrVr8ZH379IaldUnkh4QW9yxFlnNVqpXv37kBWt5PFYnEMXhuNRmw2G2azGYvFIgfmMkzmo9DJuFXjmLN2Dh90/IBRfUbpXY4ow0wmE5AVGuHh4RgMBn0LErqQiYtc0KW0S1RvU51/T/zLsf3HqFa5mt4lCSHKMJm4yAVV8K7AB3M+IDM0kxcSpQtKCOG6ZDBbR8N7DmeX9y5mrJ9B70a9Cb01VO+ShBDiGtKi0NnkLpOp83cdBncdzN6UvXqXI4QQ13DpoFBKGZRSE+z/4pVSpe4SznJe5Yh6JAqtscYra17RuxwhhLiGq3c9TdQ0zXFvAKVUslIKTdNK1RU2Q3oMYafXTqavn87wPcPp0aSH3iUJIYSDy7YolFIGIPddwKKB699Uxg292uVVAjwCeOChB0g5maJ3OUII4eCyQWFnVErlDAsb14ZHqVDeqzwT2kzgwvYLPPmh3AdKCOE6XDYoNE2zaZrmq2lazolxQ4BSO5t6eP9wRi8Zzdf/fs26/ev0LkcIIQAXDorc7F1RRq7T9ZSSkoJSyvFvypQpJVFekZnVZxaNDY15eNbDpJ7N/26dQghRUlx9MDunWCBM0zSLs4Xq1KlDSor79vFX8qnES81eYuTYkYSmhbI6erXeJQnhUqKiojh16hSDBw8mNTWV+Ph4x00NRfEo0aBQSoUCg6+zWKqmaRG51psARGuaVmq7nXJ6oucTfPnil6zxXsPmI5u57ebb9C5JCJcSExNDTEwMRqOR2NhYvcsp9Uo0KDRNMwGmwqxjDxdLdkgopYxlITDiJ8XT8v2WjFg+gt+e+I0K5eR25EJA1jwZp0+f1ruMMsWlu56UUkbADzDbxyj8gEBK8YB2Nt8KvkR1imJ4n+EM2DaAVdGr9C5JFNSmsXB6iz7b9m0HQXMKtYrJZMJqtWI0Gh1TlYaEhJCQkFD09RUhi8WCwWCQubRLgMsGhT0Ysv9PzdkBWagWiTsbdvswXgt8jcSLifx98m+aVW+md0milDGbzRiNRmbMmEFSUhKBgYHYbDYKewfm/zJn9o0wmUwYjUbMZjPR0dEyFWsxk9uMu7ij54/S8r2WtKzRknWPrsNDuc2JasINZM+PHRQURHx8PP7+/o6Db3x8PCaTCX9/f5KSkggPd8076AQEBBAdHY3RaNS7FLfm7DbjLtuiEFlqV67NOz3eYcSrIxiaPJTP3/hc75JEKeLv74/VasVmszm6cBISEggJCcFsNpOamuqY8S4qKqpYZrmLiYlh06ZNTpeJjIx01GexWK5qmQQGBpKQkCBBUYwkKNzAsLbDGH9gPEuTl/LG+Ddo7NtY75JEKWKxWK46yJrNZiIiIoiOjnbMk20wGEhISMg3KP5L11NhWioWi4Xu3btfNZhts9kICAgo8HuIwpOgcAMeHh6sXbmWDp91IHxlOD8O/RGllN5liVIiNTXVcaC12WyO7iibzeZYxs/Pj9TU/C8ALakxgsDAwGu2ZbVaGTRoUIlsv6ySoHATreq3Yva9sxllGsULMS/wVsRbepckSolBgwYRGRmJyWRi48aNBAdndVMbDAZHWKSmpuLn56djlf8vODiYqKgoDAYDycnJxMfHyzzfxUyCwo2EB4Uz6elJzH53No/f/zgt6rbQuyRRSmRf2bxx40YiIrKudw0JCcFqzbrVmtVqdXRD6S0wMLDIzp4SBSOn0LgRpRTxC+Ip/2h5xq0ZR2k7Y02UPKvVSvfu3YGsbieLxeIYvDYajdhsNsxmMxaLpVgGsoV7kNNj3dD8P+Yz5vsxzLtnHk93flrvcoSbM5myLk2yWq2Eh4dLN04Z5ez0WAkKN5SpZdIqvBW7luxiw+8buKvNXXqXJIRwc86CQrqe3JCH8iD6qWi8Wnsxcf1EMrVMvUsSQpRiEhRuqnO7zsx/fz7rjq9jQdICvcsRQpRiEhRuLDwonM6+nXlm+DMkJiXqXY4QopSSoHBjSineCnmLzCOZPLXoKdIz0/UuSQhRCklQuLn2LdrzSeIn/F39b2aulztoCiGKngRFKTA8cDiDWw1mUswklvywRO9yhBCljARFKaCUYna32bASRr80mktpl/QuSQhRikhQlBJ1/eoSGx/Lud7neCHhBb3LEWWcxWIhKCioQHeVLQnFVY+r7WdxkaAoRR4LeYzn7n6O9359j3nfzNO7HFGGBQYGOu4ZVRgxMTHXPBcWFpbn8yVRT3G9b3HtZ3GRoChlpnefjm+CL88+/Cy7U3brXY4owwp7t9nsCZRyi4iIKJJJiYrr7reutp/FQYKilCnnVY7F7yzGe4A3TyY+KVdtC7eR35wWRqPRMbtdaeCO+ym3GS+Fet/Zm/k+8wlfGc60hGlMuneS3iWVKWPHjmXLli26bLtdu3bMmTOnUOtYLBYg65duQkICM2fOxGAwYDabiYyMZPDgwVdNk5p7WtK81s3NZDI51suePyIsLAyr1UpsbCypqalYrVasVisGgwF/f3+MRiMWi4WRI0diNBodB1ir1Up0dDTt27fHYDDg5+fnuO14QevJT0xMjGPSpoSEBMft1202m+O11NRUR325Zf/NgoODiY6OxmKxOPY7Ojoas9lc4P3Mb5sF+VyKmgRFKfVE4BMs/nYxk/tNpvYXtQnvX/DpJkXZEhYWRnx8vOP24pGRkURHR2M0GomIiCA+Pp6EhAQgawKj6Ohox8Esv3Vzy349Li7OceCOiIggODjY8Tj7IJ9zatTsMYDk5GTHcyEhIWzatAmDweA4wGbPuV3QevISFRWF0Wh0hE72XBwA3bt3v2pe77CwsKsCKlv23yx72ez64+LiHK8XdD/z22ZBPpeiJkFRSiml+HTUp7T+rjWvbX2N0B6h+FVwjRnKSrvC/qLXW/ZBF8Df3/+qAyRw1S9yPz8/x8GpIOvmFBoaysiRI7HZbI7Z8wrya9/Pz89xADWZTBgMBsd6gYGBJCb+/+1rClNPbv7+/owcOZKIiAgGDRrkOJCbTKZrfqkPHjyYGTNmEB8fX+D3v57c+3m9bTr7XIqajFGUYg1rN2TN12s44XGCEV+NIDNTxitE3rKnQk1KSrpmbuzrDdY6Wze3QYMGERMTg81mu6FuEqvVek09ucOmMPXkFBoaysSJE4mPj8fX19dxyuvGjRvz3GZhQqiwCrLNkpyaVoKilAuuE8z0ztP5dsa3DBw3UO9yhIux2WwEBQUxceJEQkNDHfNlZ79W1OtmdwWZzWan05lmT6aUW2BgYL4H//+yL5A1vhAaGkpCQgKappGUlITVaiUgIOCabRYm6JyFVX77+V+3WdQkKMqA5zo9R50Kdfj2r2/55eAvepcjXEhSUtJVXTk558i+3i/m662bmpp6zQHa398fg8GQ58HT39+fU6dOXfN8zvfJHkDOWVv2tQc3Uk9OCQkJV71v9rbCw8Md4wrZ4uLimDhx4jX1Ze9Hzv3L3SVUkP283jZLmgRFGeDh4cH2NdtpfH9jwuLDOH7huN4lCRdhNBoJDg4mJibG8Ss/ODjY8Us3Pj4es9nsmDc7+0wek8lUqHVzyh4DyC00NBSr1UpMTMxVA9y53ycxMZHo6GhMJpOjjsLuS14CAgIcr5tMJtq3b3/NWUUmk4moqCgiIiIIDAzMsz6j0Yifn5+jvpCQEMxmsyPQCrqfBd1m7s+lOMhUqGXI1qNbuf312/Hb7se+xH2U8ymnd0miDDKZTI6zkoTrkKlQBQC31r6V4bWHc3T7UcbFj9O7HFGGREREOK4hcDY2IVyTBEUZEzsplkdjHuWDPR+wbOcyvcsRZURYWBg2mw2LxeKyVx+L/EnXUxn0b/q/dF3YlU2mTcRNjGNAlwF6lySE0Jl0PYmrlPMqR6wxloxfM3js9cc4fem03iUJIVyYBEUZ1bpRa5abl3O+y3keWvaQzLcthMiXBEUZ1q99Pz64/wN+2PoDIU+H6F2OEMJFSVCUcU8EPkGHox1YG7uWN1e+qXc5QggXJEEhWPvRWjq+0ZFXN7/KT/t+0rscIYSLkaAQlPcpz3fPfEcTvybcP+V+Ev4ovrtQCiHcjwSFAKBq+arE9YnjwrILDHxyIKcuXnsvGiFE2SRBIRzaNmjLpys+5UrvK/SP68/l9Mt6lySEcAESFOIqQ41D+XTQp6y3rqfTk51Iz5DTZoUo6yQoxDUGtx7MIAaRFJvE0HeG6l2OKAFms5mAgAC9yxAuSqZCFXn64rUvyPDNIO5cHHf+difP3vGs3iWJYhQcHFzguaVF2SMtCpEnDw8P4sbGMbDFQMZ+OpYX3n1B75JEMTIYDI55HYTITYJC5MvTw5PFAxbj+4svb015i+93fa93SW6ja9euLFy4EIC0tDS6du3K4sWLAbh48SJdu3YlLi4OgDNnztC1a1eWL18OwMmTJ+natSvffvstAEePHqVr166sWrUKgIMHD9K1a1fHBDdWq5WuXbvy0083dg2MzWYjIiICpZTjOYvF4pgIJyIiwjHzmslkIigoiICAAGw2G1arFaUUERERxTqHtNCXWwWFUkraxiWsgncFNv1vE82ea0bY8jD+OPyH3iWJImYwGK7pdgoLCwOyZmMLCQkhMjLS8TgxMdGxnJ+fH9HR0URHR8vtw0sxt7nNuFJqJhCoaZrTmxLJbcaLx5FzR7jro7s4nnicLyK/oO/dffUuSRQxX19fTp/OupOwzWa7aprOyMjIq+Z+NplMREdHExYWRnh4uB7liiLm9rcZV0rJlFg6u/mmm4m/P57LP13m4ZcfZu/pvXqXJIpZ9nzNSUlJpKamXvWaTGVatrhFUADBgNxXQmfBTYL535r/4d3Tm26fduPQ2UN6lySKgc1mIygoiIkTJxIaGkpwcPBVr8H/tzJmzpwpYxNlgMsHhVIqFFiqdx0iS4+gHvw4/EdOnjxJy+4t+dP6p94liSKWlJSEwWBwdD1lB4HVasVqtWKz2UhKSsJoNDq6n0Tp5tJBoZQyADZN02wFXSclJQWllOPflClTiqu8Mqt93fbMbj+bczvO0XduX7kvVCljNBoJDg4mJiYGs9lMYGAgwcHBmEwmzGYzQUFBJCcnA1mD2RaLhbCwMCwWi86Vi+Li0oPZSqlwTdNi7P9tBCJlMNt1fL31awZ/O5gWNVqQMDSB6pWq612S+A9yDmaLssfZYHaJXplt70YafJ3FUjVNi7APYJtLoCxxg/rd2o+vK31Nn6l9aBTViK3rthJQV24D4U5MJhMJCQlyeqtwymVbFPYWRM6zndrbH0cDJk3T8hxBkxZFyZv28TQmTZlEm7FtWDN6DX4V/PQuSRSQzWYjJiYGyOpyCgyUEwzLKmctCpcNityUUuFAhKZpQc6Wk6DQx3f/fMfApQNpUa0Fy/otI6COtCyEcCel4TqKcCAM8FdKTbAPcgsXcl/T+/jmwW/4c9GftGrfCusxOWVSiNLCLYJC07QYTdNCNE3z1TQtqjBnQYmS06NJD2aOmUlmy0zuM91HyrkUvUsSQhQBtwgK4T6eH/I8iR8mcujsITpEdeCnLTd2ozohhOuQoBBFrlPDTpiHmTnyyRGMvY38eUwuyhPCnUlQiGLRoV4Hli5eSpWwKnRd1JWkFDnBQAh3JUEhis3ArgP5ffLvVPapzN1P382sz2bpXZIQ4gZIUIhi1cSvCWuHrUVtVES+GYlpp0nvkoQQhVSooFBKPaGUGqiUqlJcBYnSp1G1Ruz8fScdnu7AoPhBvPvLu3qXJIQohMK2KHyBlwCbUmq3UuoDCQ5REI3rNCbxiUR6+/fm2aHP0nFoRzK1TL3LEkIUQKGCQtO0WfYr93yBUcAZsoLjtD04xktoiPxU9K6IabCJls1b8svFXxi+YjhXMq7oXZYQ4jpuaIxC07QzmqYlapr2oj042gOJQBNgr1JqQFEWKUqP8j7l+fO7P3lj3Bss2b6Eu16/i70pMlueEK6s0EGRV4tB0zQLkKRp2iggAHhIKdWtCOoTpZBSipc6vcSHPT9k0+xNtO3dFutpueWHEK6qsIPZC4B9SqlT9vGJnGEQAKBpmk3TtEFcfedXIa7xeIfHeffjd/Hs5ckdH97Brwd/1bskIUQeCtuiSNA0zQ8IARRgUkplKKUygFMASql29mWlP0Fc15hBY/hj/B9ULV+Vux+6m9HTR+tdkhAil8IGhU0pNVDTNIumaaPsoeGnaZqnpmlv2ZcxKaU+IGvAW4jralqtKWuHrKXKqSos+G4BLyW+JGdECeFCCjXDnaZpiUqpqkqpbpqmrbY/dybXMk2UUlVzPy+EM3X96nJo8yGeXfUsM9bPIGlnEp889Al1q9fVuzQhyrxCD2bbz3hafb1lbrwkUVZVKl+J2H6xvBPyDgnTEmh6V1N2n9qtd1lClHlyCw/hUpRSjL1rLLNmz8Kzmye3f3g7P+z5Qe+yhCjTJCiESxo/dDxbZ2ylQdUG9BzXk3vD7yU9I13vsoQokyQohMtq7NuYXx77hcYXGpOwJoHQpaGc/fes3mUJUeZIUAiXVsmnEnvMe5gZO5OVu1cSODeQ5WuX612WEGWK0jRN7xqKVHBwsJaUdAOT5Fw5gzUpnp27D9CjSyDe3l4k7zvCruRD9OwahKenF3v2HWHPviP0uKc9ysMT6/6jHEw5SZeOQaC8OHwslZOnz3Nrm9bg4Y3t3GX+TVPUqlMfPCuAhw8oVfQ7XUas27+OnoN6cmnbJd7/8X1Gd5JrLoQoKkqpTfZbMl2jUKfHlmoX9vLNByMZtxhOx4ChEpi+hRe/hIufQAUf+Hw5TF4GGZ+B8oCPl8Kb30L6Z1lvMedzeC8ha3mAyZ/Cop/BFpv1eMwi+Mai2B9dAzwrMdV0EYs1jRVv3AHeVVmy+gRHTmuMD78PylXnz73nyfQy0Da4E5SvBR7e+vxtXETnhp357fPfeHjuwzy5+kkspy3M6TGHSuUq6V2aEKWatCiypV/k2O6fOXj4GO3aNMXLy4sjR09yKOUYQbc2w8NDcSjlGAcPHeWO4FYoNPbuO8T+Q0foelcbyExnx65k9u5P4f5ut0JmGr8k7eAf6yFG9LsNMi6x4sct7Nh9mFdGtIH087z12Wa27TnFp8/Xgys2hs8+yNZ9V9g6I6ukvrPhwEnYYn88LLocV7TyxL1+N1Ssh2nDWSr51qPXff2hcgCUr1kmWizpmelMWjOJGYtnUN5cnpVfr6R7YHe9yxLCrTlrUUhQuBgtMwOVfhYun2T75t85bzvMna184dIR3lzwP65ctDFpUCW4eJB2z52kfjX4dnzWul2medDa34/3XuwGVZrxw+YrNG55J00De4BneX13rBjM+HQGkyZNwnu4NwseWMDwW4frXZIQbkuCopS6cOY454/volbF83A+mZdnfkYDwyUiul5GO59MtXCNwXfAB497wE23EB5zhT49OtKn/0PgFwgVauu9C//Z4bOHGbJ8CD/t/Yl2e9qxcu5KuZpbiBsgQVEGaemX2blxFeWuHKCJ70nOp2zm1kdX8bQxnXG94PxlaBXpyZuj2/HQ4IFk+t7OhfKtuMnvZr1LL7T0zHRGzhvJwnELqflITb6d9i23171d77KEcCsSFMJB+9eGOrONI7vW8MIbC3miSzpdGx1i2wG47SVY/nIj+vW7n/MVgzlbrg11AtznbvFLVi9h4taJHDl/hKcbP83MwTPx8fbRuywh3IIEhXDuymkObF3JRx99yMh7MqmnNmPacIGwdyHp7QYEdezJCc8gqNmJGvVb6F2tU6cvnWbEFyP45ulvuPnOm/l12a80NDTUuywhXJ4EhSiczDT2bf2e5V9+yJiQDLxPr2fa0rNMXgan4lpjaNKT494dqOp/L+Uqut4U6ZmZmYx6YxRLTi/Bq6YXc0PmMvy24Xh4yPWlQuRHgkL8N5kZ7Pglng3mLwnvdBZOrGfY/DR+/luxd2lvVL37ueTbnQrVb9G70qtYT1sZvmI4G6I3UFfVZVPCJmrdVEvvsoRwSRIUomilnSdx+VwO7lzNiMC9cGEvHSZBs4YGPn37aajXD/yCXOKajozMDPqO6cuqbauo0a8GMX1i6Nusr95lCeFyJChE8dE0tDN/8fb0CdzsuYeH2+4mPT2T9pO9GfdIF4aPeglqdAIPfW8CsO3YNoavGM7WP7cS8FcACZ8l0LhOY11rEsKVOAsK6bQV/41SKENLno9aycMzdsGAY5xu/i6N6lbHcH4dJHbj+MKaPNG3Gf/8/BFkpulSZttabflj5B/0qdiH5D+SuTP2Tlb+s1KXWoRwNxIUomiVr06N9mNY8XMKfWekQqdlbDt7G3EJ/3B5wxOwvDZ7TINIjJ9FRtqVEi3Nx9OHb6K+YZ1lHTVq1qDP531o/3h7dh+UWfSEcEaCQhQfr0pQfyDGsYmcOHmaNkNWQJ1exHz6Fb0fnsC5L+tB0jPYkhPQMjNLrKxOt3RiU/gmRjccTdLCJALHBLJ0x1JKWzesEEVFxihEibt0LhXLjx/QseZmOLySfrP+5fTlcqz7fDw0HgZVmpVYLSvWrWDqzqlsPraZzt6deWvAW7Rv0b7Eti+Eq5DBbOG60s6yZP4LXDr0M08E/42WmcnQD/14YOAABo58E8pXL/YS0jPTeefXd4h8IBLlpXh3+buMCh6Fp4dnsW9bCFchQSHcw8UUUrd+TLehbzCqy2VGhXjxb42e/G9vS+4b+io+FSoX6+ZXb1rNxJUT+YM/aF+rPRNbTmRA5wHFuk0hXIUEhXA7mae24LF/MSviPmJglI0fXqnCvf0fJaPhcDyrF9/9pzRNY8n2JYx6ZRQXvrvAo9GPMm/4PCr5yORIonSToBBuK/3KZRLjZ2Gstw3PI98w65srfPFHBX6On0qllo9COb9i2e7ug7sZ/vpwfqv7G/Wr1mdSu0k81uUxuQ2IKLXkOgrhtrx8ytNjyKt4domHASk0aD+CoCYVqfTXeFhRh4Uvd2Dt8tmgFe1ZU7fUv4VfY39l/WPrqZxWmZG9R9I8tDnJqclFuh0h3IG0KIR7St1M5p6PuKXvB9wRkMmSCY0h4HFO+/bHt26rIt3UxcsXefCFB0n0TCSjegZjbxvLhE4T8KtSPK0ZIfQgLQpR+vjdhsft89mx5ySz57wHlRtzaM0r1G7UmkUvtoND30JmepFsqmL5inwz7xt2T93NAy0fYOZrM6kVUIsvNn8h116IMkGCQri18pV9qd3+SeieiHfvX3l++J10anAY1vXF8nYdJo/qwulD24tkW3VuqsOSgUt457l3qN6pOg9/8zD3Lr6XDbs2FMn7C+GqXD4olFIGpdRMpVS4/Z/7TLkmSlStgDuY/tEv+I9KgU4rWLe3Gu8sWofnqraw9j6ObFpI+pXL/3k7Yx8cy0HTQeb1msfvW37n7jZ30/PFnpy6eKoI9kII1+PSYxRKKQMQr2laiP1xOBCiaVpYfuvIGIXI6UzKdqqeWArWjwiZdISLaT5sWPoiBDwBler/5/f/5+A/hI0JY3uz7Rj8DIy/dTxju46lYvmKRVC9ECXHbU+PVUrFA3Gappnsjw2An6Zp1vzWkaAQecpM5+tPXuXy/u8Z3HwbGvBkXAOGPfY0d/UZB//xKuztx7Yz7odxJE5LpPyV8sStiqNPsz4oF5iTQ4iCcOeg0ABfwA8waJpmud46EhTius7vZf+6t+gwZAFRD2YyvEcDLtcfwWnfAdzs3+6G31bTNCbHTCZ2fSxHmxylW+NujAkYQ/+O/YusdCGKi1ue9ZRjLCI4x3Px9laFEDeucmMa9n6PAylnePD5z+GmW4iLfZ0GTW/jz0X3wtHEG7ouQynF6xGvc2DhAeb1mkfS+iQG3D2Ae16+hwNnDhTDjghRMlw2KAD/7P/QNM1qb03EAbHOVkpJSUEp5fg3ZcqUYi5TuCufCpXxafIQdDfTabSZ10bdTauKSbDaSOyYm3k7si/a5ZOFfl9vT2+evv1ptk/dTucRnfnV+1eazmvKozGPknxYLtgT7sdlu56UUkYgAfDVNM1mfy4Q2KRpWr4dv9L1JP6TjMtwwMRDj43j5MmTJLxcDhoMYq9PXxoFDkTdwC08Dpw5wKQ1k1g0ehGenp5M+2Iaz9zxDBW9ZcBbuA6XGaNQSoUCg6+zWKqmaRFKKX8gOWcoZD9HjvDITYJCFJVLRzZS4dAizuxYRN1R5xnbvybTpk6BRkPAu0qh32/Z2mXMSpzF716/U7tCbXpd7sX8CfPlDCnhElwmKApLKXUaaJyrRZGoaZpvfutIUIiidunsCRbPe54Ovhtpa9jF3tSKzFzdkJdfe5v6bXoW+v1+3v8z4W+Fs2v+LmqPrM3sMbN5sPWDeChX7gkWpZ1bDmbbzQAG5Xg82P6cECWmQpUajHz5U9qO3gk9/mDT2Q58vuovPH7qBT/cwdHf53LxbMHHMjo17MSOuTuYvmg6NW6rwZDlQ2gU3oiJ8yeSWYJTwgpRUC7dogBQSk3I+VjTtChny0uLQpSES2dSqHAkHvYsYMiMXfyy24M9K5/Cs9loqNqiwO+TqWUSvyOeRwc+yqUrlwicGMjUblPp1aSXXIMhSpTbdj3dCAkKUaI0jZ+/mY914xIeaW2BzDTGmupyb5+H6T18KniWK9DbXL5ymQU/LeDdne+y9+heqq6syuuvvc6YAWMkMESJkKAQoiRcPs6ZrR8Q2H8ao7ulM35gdTIbjWCvV08Cbu1eoLdIy0jjtS9f480xb5IRlsEdt9/By3e+TO/mvWXSJFGs3HmMQgj3Ub4mVTtMZvehS4x541uo2Zkf49+mSTsjq98KhAMmyExz+hbent5MGzKNs0fOsmDkAg6fPUyfkX2o2qwq3+z8Rm5rLnQhLQohitFR6xYWzZvAuLt24ZN2kC82VmH72VZMmfUJPn7Nrrv+lYwrREyLYPlPyznb5SxBNwcxpMYQxvQdg5enVwnsgSgrpEUhhE5q+7cj8p0f8XlgL3T5js1Ha7Bq9a94f98cVvdgZ+Jc0v69mO/6Pp4+fDL5E04mnOSjvh9x4sgJngt9jrr96/LF9i/IyMwowb0RZZW0KIQoYVdsVnwOfkba37E0GHkYY9vyfDZ3HDR5Air7O1338pXLPDf7ORLSE9iTuYcG/zagi0cX5k+cT5VKhb8IUIhsMpgthAvKSLvC959Po8bFtXQwbODU2UyGflSNaZPGE9TrOfD0yXfdTC2Tr3Z9xVPPP8VR81Fqv1qb8d3HEx4Uzk3lbirBvRClhXQ9CeGCPL19uP+R1+kweh3020/yTeH8tc+G958T4au6HPo+gl0b/5fnuh7Kg4EtBnJ45WE++f4TWjRowfiE8VS7oxoho0M4fuF4Ce+NKM2kRSGEC8lMT8Pj+GpIjuG5GSt4L0Hj2NI7Mdw6ChqEglf+94XasHcDgx4cRMpNKZS/pzzD2w5ncN3BdAvsVoJ7INyVtCiEcBMeXt5Qpwd0Wkbk+9uIe3s4Bq+T8NsjPBZi4NXHgyE17/m7OjbuyOHfD/PXF38xrO0wPln+Cd2DunP3y3ez4cAGObVW3DBpUQjh6jQN7fg6nnjicep472NqaAb4tmPJjrb0HjYF31qN81xtu3U7T099mm0B27Bl2Gh+pjm96vXizWfexMcr//EPUTbJYLYQpcWV07Dvc3asnk/r0buYN8KLp58II6PRI6ja3fHI49qKC1cusHDLQiaMnMDFlIs0mNiAsXeO5dF2j2KoYCj5fRAuSbqehCgtfHyh6VO0GvUXm1YvZuiwEZDyPd+81RP/m8uT/L9n4cL+q1ap5FOJp25/itNJp4n+MpqGvg157vvnqNa4Gt2e7kZyqsy6J5yTFoUQ7i7jMj/FT+eD6BgWP3YML0/F8n9ac7FqR4Y8/RbKu9I1qyTuSGTkUyPZX28/2i0aPev35B7ve3h+yPNyT6kySrqehCgrLuwH60Luf2wmJ05f4vcZVaDhgxzw7kn9W/tdM5VryrkUPtj4AXPmzeH8V+fxf9GfFx54gaFth1LZp7JOOyH0IEEhRBmTmZHOib++ptbZr/nXGk+d0ZcZ1s2XOdMnQONhULHuVcvbztuYHDuZnyv+zOajmyn3Szla39Sazxd8TtPqTXXaC1GSJCiEKMP+PX+CJfMn0LKihTuqb+PIaXhsUTVmTHmOdveOveraDE3T+PXQr4wIH8GeQ3vQQjV6BPSgR8UePNXnKXy85Wyp0kqCQgiR5dweNix/gyEvfMYPEzJo1uAmdqt7OeLThbvvH33VWVMpZ1P4cPOHvL/2fY69foyq3asy4dUJPH7b49SqXEvHnRDFQYJCCHEVLTMDdWId7P2UZ15bTGxiOscWNaBKq0e4XHsQ5Wu2dix78fJFJi+YzPq09fx28Te8Tnhx8+abiZoZxeCOg2UGvlJCgkIIka8LZ45jWTWPTjX+gGNmer2ZSRWDH3Gxr0ODwVC+umPZv0/+zfNznud/7/4PbbRGs/rNeKDWA4y4YwS31L9Fx70Q/5UEhRCiYC4e5u3XRlHhvIXRHVPQ8GT0l/V4cMgIuoa+CJ7lATh/+TzLdi0jelM0v0b9CsfhoQ8fIjwonC4Nu0grww1JUAghCu/0Vg7//j63D/mQ1x/I5PGQKlys0Z/fTralS78xeNoHtpetXcYnaz5hfYX1nLl8hkpfVaJPvz7MeWGOjGW4EQkKIcQNy0i7QsaR1fgc/pIvv4zjoXcv89PUanTuPYK0umF41WiP8vDgYtpFFv62kIlPTORs07N43e5Fr8a9uO3ibUwcMZHyPuX13hXhhASFEKJIXDx7iu8/n86AZv/gcewHpi1L44vffUhaFkmF5sPhpiYA/HXiLz7e/DExi2M4+9lZ/ML9iBgUwYhbR8h1GS5KgkIIUfT+PcWy2JdZl/gtc8NSAJj+Y12q+99F+Pg5ULEOFy9fZOpHU9lWdRurrKvIXJ9Jlf1VmLlwJkMCh8hsfC5EgkIIUbwuHETb/yXdHnqdRobzfBKhoGYXlu9qRuf+46herxkp51IY++ZYfkj4gbP3naWid0XaHW1H6F2hPDv4WTyU3GNKTxIUQogSc+XkdnyOLOew5VPqjbAybZDi5VE9yKg3iHNVu1G1RgN+P/w7H2/6mA8f/RCtnkb9R+szrO0wetTsQec2nfXehTJJgkIIUeK0zEy2ro+n9pX11L6wkrUb99FjJpjfvotO9z8J9fqSeimNpZalfH34a37Y/APa2xqNHm7E+GfG82DrB6lWsZreu1FmSFAIIfSlaezZuIwF897ktd4pVNKO8NkGL1Zsr8Gi99/gpmZh7DxygAlRE9hdczf/aP/gediTGpYavDL9FR675zEqeFfQey9KNZm4SAihL6Vocnsob32WRKUHD0HIei5U6cSpU6eovPUxWF6DzbGP8GirW/h74ia2RGzhvpvv48T+Ezy99mlqz67NfTPvY9Zns0hLT9N7b8ocaVEIIfSjZcKJDXBgKR2HRlOlXBrfv1Qebu7FTymtaNPtcTad282S7UtY8tIS0o+kc/PLN/NQm4foUbsHxjZGmWipiEjXkxDC5WWkXSF19w/UuJjAuX+WUmPEMcK7efDui72hQSj7aUf8tgTWX17Pd/98R/rcdG5qchPj3hzHQ20eonn15nrvgluToBBCuJXMjHR+/T6WGld+p6lHIrv2HKLNi7B00q0MePgpTho68vzcmWy5uIXtVbajXdGo9GUlBj05iFcfe5XGvo313gW3I0EhhHBfmsbBbSuZ/840nulyjLrl9vPdFkXUqiosnjsWn5b3MW/d18x9cS7nO5+HRnCrz600PNKQac9No02jNnrvgVuQoBBClA6aBrbtfLXoDd6KXcnqyIv4eMEX2xpjvRjA8OfH88XBLcxf8B4HlxyEcdCxdUe6GbrRp0Uf2rdor/ceuCwJCiFE6XR2NxxawegJs/ljx3E2vQHc1JTlf7fkol9T9gdUJm5nPNtjtsMuuGPuHQxuO5j7/e+nSc0melfvUiQohBCl3uXUZMqf+hHtwHIaDjZzZxOIm1Ab6vXj498zWJ1qY1vNv9l+fDt8DDXq1SByViQPtHyARoZGepevOwkKIUSZYju+H9s/X9GIDZzc9T9qhV/glYHevPZsPw7fFMToed+xJW0/BwMOQiZUWVaFvg/35dWIV2larWze3VaCQghRZl25dI7EZe/QtOJfBKif2PjnEe6cAiuntaF9jz4s3HuayS98yoWgC9AamlVoRs0/azLxmYn0DOpZZmbrcxYUXiVdjBBClCSfCjfRa+ikrAdaJjUbfcuLB2bRvtEpqu2aTq31UD/Nh6+GD+XPm6ox9WszPy/+mZ8r/EyjDY3oZuhGuwrtiBgQgY99Vr+yRoJCCFF2KA8a3taPabH9sh6f30v1i1G03rGcW859SbNz6exLrUBLYz0eG9yL2KP7WRi9kMyfM5kyeQr92vWjk28n+rXrh18VP333pQRJ15MQQgBcOQNHf+S116exdfsulj97BZQnU76vwxnvynj29Cc2eR1nPz2HOq7oM68P/Zv3595G91LXt67e1f9nMkYhhBCFkZkBp35HO/QtXYbNo1mNC8SOBK1SY8YsyeCcoTwbWlwg+exhWAA3t76ZcVPG0bdZX5pVb6Z39TdE7h4rhBCF4eEJNe5C3TaDdX+e5/0Ve6D9+9g8buHjrw/Q6tw/7K6TyrHAu3igXQ3q1tOYYJ5A87nNKd+gPH1e6cPafWtJyygdd7p1+RaFUmoCYLM/NGiaFuVseWlRCCGK06VzqaSlrKHKubX8/OMyOr98hK/GQe+uTUg44cv02D+xNPmXS7dkUuVSFQzrDTw1/ikev/dxl56IyW27npRSE3IGg1IqEDA6CwsJCiFESdEyM/l70yoaeG2n4ulEFixZzeiPM9g5uwKGpq2I3Wjjww/3cHAgeNTwoPWV1tQ+UZvJ4ydzZ5M7XerUW3c+PXYw4AgFTdMsSqmJOtYjhBAOysOD5u17A72BSIbedpRGXRfQvNlx1JFV2Dbu5cQJONWxHn+Xr8OsZX/z/bJtdKz+I/Wr1+fWS7fS1rctEx+bSOVylfXenXy5eosiAbBpmhZmfxwOpGqaZspvHWlRCCFcgqaRvNXMtvVfMqD1UTi2ht4zLnEwFVa92wLzJQ9mzdvJjr81fMb50KVhF5rYmtDn9j70vL3kL/Rz564nf2CT/eEMwOosJECCQgjhojIus2bFPGz71jOgRTKc2UHAOLi9qTfPj72Fpaln+OS1w5ysBI2ebESvJr1oeLYhj/R4hNp+tYu9PLcNCnAMZocARiBK07RIZ8vXqVNHO3LkiOPx5MmTmTJlSrHWKIQQhZV5/iCL3ptMHa9kejTawYUzp6g2CiaEVea2+2rx0YGDmKdcIS3Ygy4ju9AjoAcBaQEM7DKwWKZ/ddsxCqVUNBCtaVqUUsoIxCul/LO7ovJSp04dUlJSSq5IIYS4AR6V6/No5MdZD7RM0vetY9Lo6YQ0PUV7j+009UgjOBPm32HgcvpfLPpqDVtnQ5XBVej/YH/uqXsP7au3p1XjVsVea4kGhVIqlKwBamdSNU2LsJ/hZNM0zQKgaZpZKdUY2FvcdQohRIlSHlRt3JWX5nbNepx2ngob4xgd9h73t75AfZ9/aJgBIyvDu+3KkXTYxGsrPmXfp9D8peb0v6c/d9e6m3ua3kPFchWLvLwSDQr7+ILTMYYc/IBTuda3KaXMRV6YEEK4Eu/K+N/1OG/f9XjW40tHqFshmgHJnzMo4AJD0i5S8zJE11BMbZTKiu1RjHjvTQL738sPj/xQ5OW47JXZmqaZyRqbcFBKGQCrLgUJIYReKtxMhwFT+PCbf/AOPQT3/UWrruH069yYUMO/LKmVSeg5WOBXPFeCu/QYBRChlJoJJGc/cb3BbCGEKNWUgqrN6Ts6mr6jybov1WkLfc++T+NWXYpnk65+1lNhyemxQghReHJTQCGEEDdMgkIIIYRTEhRCCCGckqAQQgjhlARFLnK7D9cjn4nrkc/ENRXX5yJnPeWilKK0/U3cnXwmrkc+E9f0Xz4XOetJCCHEDZOgEEII4VSp63pSSp0A9v+Ht6gDyO1nXYt8Jq5HPhPX9F8+l4aaptXI64VSFxRCCCGKlnQ9CSGEcEqCQgghhFMSFEIIIZxy9duMlxj73NxWsiZMQtO0GH0rKtuUUuFAEBBvfyoMmKlpmsxHUoLsc8CEA9XyusW/fG9KnrPPpLi+NxIUgH3Oi432GfhQSs1USoVmPxa6GUTWF8ICjJSQKFn2eeoNQEA+r8v3poRd7zOxK/LvjQRFlvBcyRwHzKTg07aKYqBpmq/eNZRl9lkmUUq1J+vglJt8b0pYAT6TYvnelPkxCqVUYB5P2wBjCZcihNuQ703ZIi2KrL7V1FzP5X4sdGDvb01F+r9dkXxvXFRxfG8kKPJpvkHWoJGmabaSK0XkkATYsvtXlVLxSqlU6f92GYb8XpDvja6K5XtT5rueyGou++V6LvdjUcI0TbPkGoTbCEzUqx5xDRvyvXE5xfW9kaDIaqIZcj1nAJBfRfqxn92RkxXIq19c6EO+Ny6ouL43ZT4oNE2zkPXrKCc/wFzy1QgApZQ/kGA/XzwnOT3WRcj3xvUU5/emzAeF3VKlVGiOxyFAtF7FlHX2pnNkrl+mg8k69VK4DvneuJDi/N7I3WPt7FeYWgB/kDNs9Gb/dZR9EKoGJMtnUrLsp8AagQj7U9GA2d6ayF5Gvjcl6HqfSXF9byQohBBCOCVdT0IIIZySoBBCCOGUBIUQQginJCiEEEI4JUEhhBDCKQkKIYQQTklQCCGEcEqCQogSZL8gSgi3IkEhRMmKz3XbCyFcnlyZLUQJUkppgK/cYVW4E2lRCFFC7LeAtkpICHcjLQohipk9IELIulmbjaxbcW+U2fqEu5CgEKKEKKU2ATMkIIS7kaAQooTI+IRwVzJGIUQJkPEJ4c4kKIQoGSHINKHCTUlQCFEyjMAmAKWUwd7CEMItSFAIUTIM/H+LIlzTNGldCLchg9lClAD71dj+gJWsOY5t+lYkRMFJUAghhHBKup6EEEI4JUEhhBDCKQkKIYQQTklQCCGEcEqCQgghhFMSFEIIIZySoBBCCOGUBIUQQginJCiEEEI49X+djC4t8NfzfwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"y0 = -1.0\n",
"jax_result = jodeint(dy_dt_jax, y0, t_domain)\n",
"correct_result = analytical_ode_result_jax(y0, t_domain)\n",
"__=ax.plot(t_domain, correct_result, color='green')\n",
"__=ax.plot(t_domain, jax_result, ':', color='k')\n",
"\n",
"y0 = -5.0\n",
"jax_result = jodeint(dy_dt_jax, y0, t_domain)\n",
"correct_result = analytical_ode_result_jax(y0, t_domain)\n",
"__=ax.plot(t_domain, correct_result, color='orange')\n",
"__=ax.plot(t_domain, jax_result, ':', color='k')\n",
"\n",
"solid_line=mlines.Line2D([],[],ls='-',c='k',label=r'${\\rm analytical\\ solution}$')\n",
"jax_line=mlines.Line2D([],[],ls=':',c='k',label=r'${\\rm jax}$')\n",
"green_line=mlines.Line2D([],[],ls='-',c='green',label=r'$y_0 = -1$')\n",
"orange_line=mlines.Line2D([],[],ls='-',c='orange',label=r'$y_0 = -5$')\n",
"leg=ax.legend(handles=[green_line, orange_line, solid_line, jax_line])\n",
"\n",
"xlabel = ax.set_xlabel(r'$t$')\n",
"ylabel = ax.set_ylabel(r'$y$')\n",
"title = ax.set_title(r'${\\rm JAX\\ ODE\\ solution}$')"
]
},
{
"cell_type": "markdown",
"id": "lasting-source",
"metadata": {},
"source": [
"# Differentiating through the solver\n",
"\n",
"Besides raw performance, the reason to use JAX instead of scipy is to put autodiff to use. So let's use JAX to calculate $\\partial y(t) / \\partial y_0,$ the derivative of the solution with respect to the initial condition, $y_0.$ First we'll work out the correct result analytically, and then we'll use JAX in two different ways to compute the result.\n",
"\n",
"$$y(t) = -\\log(t^3 - A)$$\n",
"\n",
"$$A = -e^{-y_0}$$\n",
"\n",
"$$\\frac{\\partial y(t)}{\\partial y_0} = \\frac{\\partial y(t)}{\\partial A}\\cdot\\frac{\\partial A}{\\partial y_0}=\\dots=\\frac{A}{A-t^3}$$\n",
"\n",
"Let's first implement this solution analytically, and then move on to using JAX to calculate the same result based on automatic differentiation."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "brave-processing",
"metadata": {},
"outputs": [],
"source": [
"@jjit\n",
"def analytical_dy_dy0(y0, t):\n",
" a = get_a_from_y0_jax(y0)\n",
" return a / (a - t**3)"
]
},
{
"cell_type": "markdown",
"id": "offshore-diary",
"metadata": {},
"source": [
"## Calculating $\\partial y(t) / \\partial y_0$ by differentiating through the solver\n",
"\n",
"When using the `grad` function to transform a function into its derivative, JAX requires that the function being differentiated returns a scalar value. So if you want your derivative function to accept/return an array, then you need to first define the scalar version of the function, operate on it with `grad`, and then operate on that with `vmap`. Here I'll just implement this directly without further comment, but this is a standard rigamarole that you can read more about in [the JAX docs](https://jax.readthedocs.io/en/latest/) and/or in [this gist](https://gist.github.com/aphearin/208f7a53cb156a49f72040b2b1d54320). "
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "familiar-saint",
"metadata": {},
"outputs": [],
"source": [
"from jax import grad\n",
"from jax import vmap\n",
"\n",
"def ode_result_jax(y_init, t):\n",
" return jodeint(dy_dt_jax, y_init, t)[-1]\n",
"\n",
"ode_result_jax_grad = grad(ode_result_jax, argnums=0)\n",
"ode_result_jax_grad_vmap = vmap(ode_result_jax_grad, in_axes=[None, 0])"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "different-macro",
"metadata": {},
"outputs": [],
"source": [
"t_domain_vmap = np.zeros((t_domain.size, 2))\n",
"t_domain_vmap[:, 1] = t_domain"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "painful-whale",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEiCAYAAADnMZWTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABGp0lEQVR4nO3deVyU5f7/8dfFjogMuC+5DC5paglotmo6tNp2Amz9Vp6E6rQvktVJf51TJp1TtpyT4KlTnTaFds0KrKysDBg1M5dkNDXccUQFZLt+f8xAOLLLzD0Dn+fjwaNm7vue+z3g3J+5r+u+7ktprRFCCCEa4md0ACGEEN5NCoUQQohGSaEQQgjRKCkUQgghGiWFQgghRKOkUAghWkUpZTqR5fWsbz6RPMJ9pFB0IEopi1JqrlLqgFIqXykV08i6pkZeI10ppZVSBa4f7jrLMht7/ROhlIpRSmUrpbLd8fqeopQyK6VmOH/c8ruqsy+T82+X0BYHZKVUstba3sRqSc0tFs71Uhta1pbZRcspGUfR8TgPsJla64wGlicAKVrr+EZeIx2I01rHujwf43y+3tduK879zG0sYwPbmeo7wCml5gImrXVKG0VsTpZ8YDKQDKC1TnPTfizO/81z/jcJiG3te1VKJQM5Wmuby/Pmep6bq7WutwC4rDcDsGqtc9yZXbSOnFGI+pgBS2PfBp0fVLPzoFGXxd1Fwsneyu2SGnh+IZDeytdsMefv1qS1tmut09xYJExANmBz7ste8/dxfiFozetFuxYEp8x6XjO7mfuJr6dItGl20XpSKMQxnB9Oq/OnoYNqjUQgvaagOD+8nigSJ6LeMxCttVVrbfVgjiigvoNtm3KePSXWc2DPB1rzrTwJR1GtTwxwzMHeefCf2tgLOpuTjvvduyG7aCUpFMKVxfnhnkMTH0bnelnAgpq242a0WxumpnnJ6ByeprXOamBRa9r74+srqM4mIlsDf/+iJvoWUmjgbK6Ns4tWCjA6gPBOWussZ4f0ce3OLqYDW4DU1rQbO5uual7fBETVaV5IBoqcy6IaatJy9lcsAIpq+iyc26bi6MfIcJ7tmHA0l81wbpqhtbY7D2LpzvcdX+d1691/Tf8IjuavOTjODkzA2Ga2x1twnNnUZLHVHBCb2OcCHG31mTgOlPFa68Rm7i+GP37PNuoUTOdyk/P1Uuo8PwOo29xTu0097yUBsDuLca7LAT4bsNDw2WZMQ//GmsrekvziBGit5aeD/eD44CbX87wZMNd5nI7jQNvU680FDuBoc29JjmSX/Zlq9ofzYOiyLNsla93HMXUf18mf3Ng6DS1rxv4tOJpAzC7bWJr53s315G32Pp3LZjRjP+mu6+E4qGc792cCEpzPH/M3dO4rweV3VO+/B9d16/ndNrSdpaH30VT2Or+jZuWXn9b/SNOTqMv17CEd5xU5DXF+y8zF8U13QSv2WfsNUDuaLRY6vz0fk8W5zFZP53lj7K1Zt5n7L8JxUKr7+7LRyiaRlu5TOzvBm3jNZBxXoLmuZ8VxgM52Ls+quTzXuc+av6trn0MUsL+B3R3XP1GHnYZ/L4nUc6bRzOw0N7/zkuqEmp8GsogGSNOTqCu+vuv5lVI1/Rb1SdZapymlrEBBE+seQzuahDKVUhrHBzrT+Vzd5qi6CoDYep5va3HN3L/dgH22pAM8nfr7mWoO2ov0H30KKcAi1zz62D4HM/W85yb6J8BZ4BpYFtXAdk1mh9p+Mmg6/0ztbKZzjvXJaSSvcCGFQgC132jnuH54lFJdcXwIjzv4Ow8QWQBaa5tSKhVHc0ZkM/dp0lonOr/9xQGpSqlYHAfHhkQ157Wbuf+G+l9Mzdx/UYNrtVxz92lvzovVKfj1Fe0UnP0zdZ5LwvHtvkZ8PdsWUf/vv75166r3Ci/nN/vjrqBqRXZoJL9zP3X3X4DjrKShjnLhQpqeOgDnCGpLnafq+7CbG/iGtRBHu3B9XJtK0nA0lTR3PMJM53Z2rXWOdnQkm3F8wOtrqojmjyYHV3aOf1+mJvbf0Gjo1uz/RLlln66F0Nlxb6HOKOg6YzrqHpjrNu/UsDvzuKrpO6kdRe2y3ET9BW6qbviqpmZlb2Z+M8c2mdmBsQ3tVxxPCkXHYMU5stX5oTJT5zTd+Vy939S141JIm2u7rlIqXdd/NUkikNzcW1LU0+dgq7PPmDrrmWh8xHd9zRtml+dc+xDs9b1QC/Zf3+/MNUOztGCfzXr9Rl4vE5hc90tBPWeRMdTf55BH/cXMVGfd5HqaHs04+rHq7qPB99GS7C3M75pZNJM0PXUMqThGWttxnJLXftic3/7mAialVJ52uUbeeSCPAuY6v9HZcZz+xyilMus5KNR8m1ymlJpD/c0ENfbjOCDUbGNyZsHZJDWjzvX3Zhy3u6j5ZjkXiFNKzdCOkc12Z9tzTf+GCcc3yhSllE1rneWyTpH+45LUuq+XrLXOaGL/MTjOhsw1+3cW0ppLRI+7FYXL7/SY/DW/o2bu85htGtqHUyyOv1tNU1409Q9gA0ez31z+6BM5rs/B+fsz1bctkKCUslH/JbDxHH8fpyQaHwnfkuzNyd/VZX07otnkXk9CiGM4mw4L6ruqylmkcly/UDTxepnaZbyHUipbt/A+XS3Y3zH5nUW8doyL69gV0TRpehKig3P2YcU4/9+E49t+Q018GTRxSw6X17bg0mHtPGNqs9uXNJXfWRDqNplF03izlHAhhUKIDsx5YK17hrAAmN5Qs5bz+eNuL9+I+Hq+uSfQRjdgbEH+dOcYCguOQYyuy0UjpOlJiA5O/XHbEDP13Oq7gW1mNHPA3yLXg3JbNzu1Jr9oGSkUQohWUQ3M7dHU8kbGrwgvJYVCCCFEo9rd5bHdunXTAwcONDqGEEL4lPz8/H1a6+71LWt3hWLgwIHk5eU1vaIQQohaSqnfGlomVz0JIYRolBQKIYQQjZJCIYQQolEe76NwDpBJBrrq5k0bOQPHKM4ocMxh4NaAQgghjuHRMwrnqEgLjiH0pmasPxfnPVmcBSJaZqcSQgjP8mihcM45kEXz79yY7DL8fyH1z3olhBDCTby2j6KB+Qzs/HEbayGEEB7gzeMoojh+qsm2nHryGPn5T1Ja8BohAcGEBISwOOcg+4uCePqus0AF8PbnNnbsK+PBaRMhMJzl1p0cPurPJRdfCKG9OKoiCQrrhvLz2torhBCt4s2FwtTQgsbuMVNYWIhSqvbxrFmzmD17dpM7O3LgF6KPbqaqtBqAzevAaoOjl+0g2M+PTz8vxrqlkgfHfgnAs8/A1r1wSfBTAFz1NBSV+PPdvFMhPJr3cv2I6H0Kky+9ETqdBHUyCSGEL/HmQmHn+Kkm652us64+ffpQWFjY4p2da3kDeAOtNfYyO5H2J9n+w6uErN/HlKFTeOvztwgPDofqSqg8TPr4AkqKd0OvYCjdxdVXL+aIfSeEhECRlceeKcDcAyZXPwahfXh8SVdOjTuHy2+YAWEDWpxPCCGMYshNAZ1XM5m01g12TDv7KPK11qqx51zFxcXptrqFx9HKo8z7YR4PPfMQp/qdyuqs1c3e9vCBXdi3raRf6A6qdn1NdEIWSeOqSbsWtGkM877qxeXX3Y/51MltklUIIU6EUipfax1X3zKvbVB3TkRid3k6Cg/OTBUcEEzq2anE+cexJm8Nn6z/pNnbdo7sRb9TL4ehf8H/3IVs2VXB4y+vhjFPs2mn5r6nl/LFixbImUjFplcpPeS27hchhDghXlUolFJml3ESi1wex9NGM2O1RM6rOQy9byi3fXYbxUeLW/Uays+PkJ6nwvAHGPbnVez4NZdrkv8flOzg3Rdupn+/bmxafC+U29s2vBBCnCBPD7iLcY60TgAsSqkZLpfBJlBnnISzacqslLI4Z7EqMGJC9IiwCF698lW279nOFQ9f0Sav2XdwHGFjH4NLNzFkyktcNaEfg+3z4IP+WN+exsG929pkP0IIcaLa3cRFbdlH4ersG89mxesreOPzN7gu/rq238GBNVSueYLBiZkM7xfI0oXPwuBk8Ats+30J4cNsNhupqamkpKRgscjQqrbgk30U3ujdee8Sek8oy8uWu2cHkacSMHER7y16g79PHw15d1D24WhWfrrAPfsTwgfl5ORgs9mw2WQ2VU+RQtECPSN7kjghkUXrFlFaUeq2/cScdx2xt+XCuR8w772dnHFxMhuzboCKw27bpxC+wmKxYLFYiIpq8mp50Ua8eRyFV7pq0FW8/uTrzOo8i7Q70ty3I6Wg3+X85blf6DfkJoaVvwmffMPBkQuIiI53337FCbvn03tYvWu1Ifs+rddpzLtwXou2ycrKwmazYbFYiIlxdBnGx8eTnZ3thoTCF0mhaKELR1xIwNYAPv7+Y/cWCqfwqD5c/9jnsPc7Nr87lXHXnc+Lj03l2gffAiUnhOLE5OTkYLFYmDNnDnl5ecTExGC321s8nXBqapMzBjB16tTaQiR8ixSKFgoKDOK+/93HP1f+k12Hd9Grcy/P7Lj7mXS/agVXLZnIOaEL4csDcNZbENzVM/sXzdbSb/RGMpvNmEwmcnJySElxXHCYl5dX20GclZWF2WwmLy+P5OTkBl9n7ty5HskrjCGFohVujLmRtB/SeGPVGzxwzgMe229E9/4s+KgAChagc+/g/900hKR7XmPE6Zd6LINoX8xmMzabDbvdjtlsBiA7O5v4+HhycnIoKioiIcExlCktLY0ZM2a0eYaMjAzy8/MbXSc1NbU2n/A8KRStMKL7CHqu6MljrzzGA795rlAAjr6LwcnsKuvN/NuvgMoEZj/3AfS5yLM5RLthtVqPucS05uwiPT2d+HhHf5jJZCI7O7vBQnEiTU+NnakI7yCFopUuOPsCXi97nfwd+cT2i/X4/nuPvJTV+Svpsf4WWD6FI6c8R9joOzyeQ/i+oqIioqOjAbDb7dhsNsxmM3a7vXadqKgoiooavs2MND21b9Ib2kpp96XBZMjeatyVIb3McfhdsIJ9oecx5oI7mffQFGhnAyiF+yUlJVFQUEBWVhZz5swhLs4x5spkMtUWi6KiIq+5HNVqtZKWlkZeXh5z584lIyPD6EjtnhSKVurZuScjuo5gad5SY4MEhBFx0UdMGBvNuJAlsOoBKRaixdLT02v7Imo6tePj42vPImw2W20zlNFiYmKYMWMGBw4cIDs7W5quPEAKxQkI+jKIr1O/pqSsxNAcgcGdWPDhJs6cchdseIbc1xPR1dWGZhK+wWazMXmy41b3drsdq9VaWzAsFgt2u52cnBysVqtbOrKFb5B7PZ2ApxY+xcx3ZvLV018xYfAEj+yzUVqz6p2bibvuNebddw53Pr1cZtYTTcrKctxn02azkZycjMlkMjaQMERj93qSzuwTMG3KNGZumMn3u773jkKhFKdNfYWX1hVynTkb8u+G2OekWIhG1ZxBCNEQaXo6AT3CejAsZBjvffqe0VFqKT8/kv/2GWGn3Uf5Ly/wztwrjI4khPBxUihOUNC3QeQ+nWt4P8UxlIIx/yBj7VlcM/Mj8jLvNTqREMKHSaE4QdOSp8ENsGr3KqOjHEspbnvyC3KePYe4inmw5X9GJxJC+CgpFCfo2knXwknwzfZvjI5yHP/AICbfmQ09J/Fz5k3kvDPH6EhCCB8kheIE9QjrgbnMzFtvvWV0lPr5B8O5H3DXm534y32PULnPy858hBBeTwpFGwj7OYy1r6yl9Kj7JjM6IYHh/O/dL8ie3Z2Ab6+A0l1GJxJC+BApFG3gL/f9Be6CNXvXGB2lQX0Hj6V/wlJ02V7+99hZlBTvMzqSEMJHSKFoA1eOuxLC4Ntt3xodpXFRMayJfIIb/2lj/iPxcqsPIUSzSKFoAz3CetC1oCuLXl9kdJQmnXbBvXz9WjL3nL4aNvzT6DhCCB8gI7PbSNCGIH7a/ZPRMZrl7Ovnw7dFFH07g81bghl38Z1GRxKi2dLS0ti/fz9Tp06lqKiIzMxM0tPTjY7VrskZRRv58+N/5ujNRyk+Wmx0lKYpBeP/y80vd+aK6++mbN96oxMJ0SIZGRlMnjyZ9PR0mQvDA+SMoo2cGX0mrIRVO1cxYaAX3PepKYGdeSb9A/YuvpyQ3Ovg/O/AP8ToVEI0yWQyceDAAaNjdChSKNrI6O6j4TN4JeAVJjziA4UCiD51EtFd34KvL+PX929mSMLbRkdqH/LvgQOrjdl35GkQO69Fm2RlZWGz2bBYLLVTlcbHx5OdbdykXM1htVoxmUwyl7YHSKFoI31NffH/1R/raqvRUVqm36V8deRqLDe8w6IXe/Kn2+YZnUh4UE5ODhaLhTlz5pCXl0dMTAx2u52W3qr/RObMbo2srCwsFgs5OTnS/OQJWut29RMbG6uNMuXNKfrkF082bP+tVXG0RD/+f/118X9DtbavMzqO8KCCggKttdYxMTG1/5+dna0TEhK01lpnZmbq/Px8nZ6ebljGppjNZp2dnW10DJ8H5OkGjquGnFEopWYANiDKWawanfTWub7d+dCktU5za8BWiusTx5Jfl3C4/DCdgzobHafZAoJC+etL38HSMVQtT6Jy8rcEh5mMjiU8wGw2Y7PZsNvttU042dnZxMfHk5OTQ1FRUe18FWlpaW6Z5S4jI4P8/PxG10lNTa3NZ7VajzkziYmJITs7G4vF0ubZhIPHC4VSai6Qq7XOqnmslEqoeVzP+jPqFgalVIzrc96id3lv9Fuat2PeZvql042O0zKd+lIR9zIXXHwZw0ecw78y1xqdSHiI1Wo95iCbk5NDSkoK6enptfNkm0wmsrOzGywUJ9L01JI5r61WK5MnTz6mM9tutxMdHd3s1xAtZ8QZRbLWuu6/qoXAXKDeQgFMBWqLgtbaqpSa6cZ8rRY3MA6KIL8gn+n4WKEAAgdcytlnno7ZbyVs/wBOusLoSMIDioqKag+0drsdm82G2WzGbrfXrhMVFUVRUVGDr+GpPoKYmJjj9mWz2UhKSvLI/jsqjxYKpVR9PVl2oLFzxiKlVKbWOtH5Gsk4iovXiRkaQ6+HelE60EtvDtgMj89fDp+fCSunQVQshJ1kdCThZklJSaSmppKVlUVubi5xcY5pk00mU22xKCoqIioqysCUf4iLiyMtLQ2TyURBQQGZmZkyz7ebefqMIgpw/VrS8NcUhxQgXyl1AJgD2BpqpvIGsb1jyS9svL3Vq/kHw1nv8MmcUfz3+RjeWfY7/oFBRqcSblYzsjk3N5eUlBTAcYmszWYDHN/aa5qhjBYTE9NmV0+J5vH0yGxTQwuUUvUu01rbcBSIPBxNVGMb20FhYSFKqdqf2bNntzpsa4RsCmHd7HXste/16H7bVJch7DZdy+Zt+9j/wxNGpxFuZLPZmDx5MuBodrJarbWd1xaLBbvdTk5ODlar1S0d2cI3KO3BO4gqpSxAptY6ss5zZqAAiNRa2+vZJh1Id/ZNWIBMIKemKcpVXFycbuk14G3p8VceZ1baLD7834dcNvYyw3KcKF1dTeXyqwjcvQQuyIXIU42OJNwkK8txgm6z2UhOTpZmnA5KKZWvtY6rb5mnm56KOP6swgTQQJGIAexaa6tznRyl1CBgi1tTnoCbr7qZWdtnsU1vMzrKCVF+fgSesYCyD0by3H3nc/eLGwmRS2bbpZozCCEa4tGmJ+cB3+7ydBSQ08AmUcB+l9ewN7K+4fp16UfX0K6s3rna6CgnLqQbPwTcx8zX9rDkX9cbnUYIYRAj7h67SClV9ytMPFB7j2CllLlmudY6x7mcOstNOAbreSWlFEFLg8hMzTQ6SpuYmDCDdQuv5ap+S2D3l0bHEUIYwOOFQmudApiVUhbnpa4FLlcxJeC40qlGinNQXrJz/SSXcRheZ9gpwyjtU0q1rjY6SpsYfuUCCB/C5vev51BRodFxhBAe5tHObE8wujMbYEH+ApIXJ1NwVwHmyPZxZ8v9G5cyaMzF3DRlBM8vWmd0HCFEG2usM1smLnKDUT1HgQbrdh+7k2wjug67iBceupgZ5/wCv39idBwhhAdJoXCD6PBomAvzX5hvdJQ2deMj79HPPBJ+vIWKI3uMjiOE8BApFG7QPaI7EedEUNmr0ugobcs/GD3+v/zfM7uYnnSm0WmEEB4ihcJNzr35XPb29uHR2Q1QXeMYMvpcojsVoHcsNjqO8FJWq5XY2Nhm3VXWE9yVx9vep7tIoXCTUT1GseG3DRwpO2J0lDb31+c/4683jUTlpkC5zF0sjhcTE1N7z6iWyMg4fmqaxMTEep/3RB53va673qe7SKFwk9K1pVSnVbNkxRKjo7Q9/2A441V+WLuLWbedZ3Qa4aVaerfZmgmUXKWkpLTJpETuuvutt71Pd5BC4SYXT7gYzofderfRUdwjKpbF28fz6odrOPDLIqPTiHagoTktLBZL7ex27YEvvk9DpkLtCCaMnkDg2YHs0DuMjuI2f31uKTPix9Jl0/0w5EII7GJ0JK9wzz33sHr1akP2fdpppzFv3rwWbWO1Oi7jttlsZGdnM3fuXEwmEzk5OaSmpjJ16tRjpkl1nZa0vm1dZWVl1W5XM39EYmIiNpuNBQsWUFRUhM1mw2azYTKZMJvNWCwWrFYr06dPx2Kx1B5gbTYb6enpjB07FpPJRFRUVO1tx5ubpyEZGRm1kzZlZ2fX3n7dbrfXLisqKqrN56rmdxYXF0d6ejpWq7X2faenp5OTk9Ps99nQPpvzd2lrUijcJNA/kKHhQ/l+9fcuNyFpP4I7dSF40utUf3YGXy24nkm3f2R0JNEKiYmJZGZm1t4cMDU1lfT0dCwWCykpKWRmZpKdnQ04JjBKT0+vPZg1tK2rmuULFy6sPXCnpKQQFxdX+7jmIF93atSaPoCCgoLa5+Lj48nPz8dkMtUeYGvm3G5unvqkpaVhsVhqi07NXBwAkydPPmZe78TExGMKVI2a31nNujX5Fy5cWLu8ue+zoX025+/S1qRQuNGR94+w4acN8KDRSdyo2+n8a81E7nr6Y/IGpxN7ftt3GPqaln6jN1rNQRfAbDYfc4AEjvlGHhUVVXtwas62dSUkJDB9+nTsdnvt7HnN+bYfFRVVewDNysrCZDLVbhcTE8OyZctalceV2Wxm+vTppKSkkJSUVHsgz8rKOu6b+tSpU5kzZw6ZmW13TzfX99nUPhv7u7Q16aNwo4uuuYiqi6ooKmlqEj/f9udHF/H2/T2JKZ0LlSVGxxGtUDMVal5e3nFzYzfVWdvYtq6SkpLIyMjAbre3qpnEZrMdl8e12LQkT10JCQnMnDmTzMxMIiMjay95zc3NrXefLSlCLdWcfXpyalopFG40ZfIUGAbr9rbveyN16tKNq2csQh3ZwtG8mUbHES1gt9uJjY1l5syZJCQk1M6XXbOsrbetaQrKyclpdDrTmsmUXMXExDR48D+R9wKO/oWEhASys7PRWpOXl4fNZiM6Ovq4fbak0DVWrBp6nye6z7YmhcKNRnYfCbsgJ9drp89oOz3OJa/iKgZd9DwrP3vZ6DSimfLy8o5pyqk7R3ZT35ib2raoqOi4A7TZbMZkMtV78DSbzezfv/+45+u+Tk0Hct1sNWMPWpOnruzs7GNet2ZfycnJtf0KNRYuXMjMmTOPy1fzPuq+P9cmoea8z6b26WlSKNyoX5d+qFcVWa/U/62hvRk6ZR7jhoYQtnkOVB01Oo5oBovFQlxcHBkZGbXf8uPi4mq/6WZmZpKTk1M7b3bNlTxZWVkt2raumj4AVwkJCdhsNjIyMo7p4HZ9nWXLlpGenk5WVlZtjpa+l/pER0fXLs/KymLs2LHHXVWUlZVFWloaKSkpxMTE1JvPYrEQFRVVmy8+Pp6cnJzagtbc99ncfbr+XdxBbjPuZqMeGEVI9xByU3ONjuIZv38Cyy+BkbNg9Gyj0wgvlJWVJdOveiG5zbiBzpp4FgUU0N4KcoP6XkxJr6tJffRx1q541+g0wkukpKTUjiForG9CeCcpFG42KHgQB6wH2LBjg9FRPKZk6Gz+uxw+/e+9UF1ldBzhBRITE7Hb7VitVq8dfSwaJuMo3Cx4fzAsgg8u/oDhNw43Oo5HdOs3jPVfpdN1fTJseh5OvtfoSMJg3noPI9E8ckbhZldMugKSwX+Av9FRPKrrabdAnyls/mwmv/3yjdFxhBAnQAqFmw3sMZCeQ3qy8eBGo6N4llKUjX6Ws2eVc2/KVdBR+miEaIekUHhAv0P9+PLDL42O4XEhUYN59Z938MLUvWD7r9FxhBCtJIXCAypXV7LlzS1UVrWzqVGb4cJp8+h78rlgvY+jB7YYHUcI0QpSKDzghjtugHthe/F2o6N4nvKDcQu4+V+Huebys6QJSggfJIXCA848+UwIbf/3fGpQl6GcOv4iYnrupPq3trvbphDCM6RQeMCI7iNgJbz7cccdgHbPk+/z6LQY/Kx3wdH2fTddIdobKRQeEBESgf/3/nyztANfJuoXAONf4fu1+5h9+ySj0wghWkAKhYecN/c8wq4KMzqGsSJPZcmO8fz3wzUc+EWaoITwFYYUCqXUDKVUglIqWSmV3Iz1TUqpuTXrK6V87mYxYwaNYcP+DVRWd7wrn+r663Of8PO/hxC56QGoOGR0HOGUk5NDdHS00TGEl/J4oVBKzQVsWussrXUGEK2UavBWkkopE5CptU51rg/gc7PjdCvtRvnicpavWW50FEMFd+pC+HmvUX14G1+9/H9GxxFOcXFxzZ5bWnQ8RpxRJGut6940fSHQ2ETLC4C6/4IXAanuCOZOfYP7ghW+XvW10VGM1/0Mnl91Lufd9gFrvnzF6DQCxzSbcj8m0RCPFooGmozsQGP/QhOAHKWUWSkVo7W2a63dN1mtm1wx8QqYCX5m6RYCmP5YJu880J3RxU+2y3m2J06cyKuvvgpARUUFEydO5I033gCgpKSEiRMnsnDhQgAOHjzIxIkTee+99wDYt28fEydO5OOPPwZg165dTJw4kU8//RSA7du3M3HixNoJbmw2GxMnTmT58tadrdrtdlJSUlBK1T5ntVprJ8JJSUmpnXktKyuL2NhYoqOjsdvt2Gw2lFKkpKS4dQ5pYSxPH7WiANdrIxu8VrJOYYmr81ymsznKp4QFhzG422DW7llrdBSvEBbRg6kPLEQdKaAs1+daEtsVk8l0XLNTYmIi4JiNLT4+ntTU1NrHy5Ytq10vKiqK9PR00tPT5fbh7ZnW2mM/OM4OClyeMwEaMDWwvgYsLs9lNrSP3r17a+c2GtCzZs3S3iLujjjd5YwuRsfwKrmv/kn3MqG//yTD6Cgdnslkqv3/AwcO1P5/fn6+tlgsx6ybmZmpLRaLTk9P91Q84WZAnm7guOrpMwo7jrOKulwfu64PUHduUxuOYlGvPn36HPMGZ8+e3YqY7tGltAvFW4opLi02OorXGHbpC4wfFkq4bQ5UlRkdR9RRM19zXl4eRUXHnvjLVKYdi6cLRRGOM4i6TABaa3s969vqWWaH2quhfMod998Bt8LGog52y/FGhEf14f33P+CUyC2wdrbRcQSOPovY2FhmzpxJQkICcXFxxywDRx9Gamoqc+fOlb6JDsCjhUJrbeWPs4QaUUBOA+vbALtLUTAB9gYKi1cb1XMUAD/t/sngJF6m9/mU9r2JBx6dy8rPXjY6TYeXl5eHyWTCZDIB1BYCm82GzWbDbreTl5eHxWIhPT29tj9DtF9GXIKzyGXcRDx1Ln91Xt1Ud/kcIKnO46nO53yOOdKM/4f+vPy8HAxdVYz8G5k/+vPFWzOlCcpgFouFuLg4MjIyyMnJISYmhri4OLKyssjJySE2NpaCggLA0ZlttVpJTEzEarUanFy4i9IG3PZZKTUDsAJmAP3HQLqaZfFa63iX52pprdMaeu24uDidl5fX0GLDdY3rSpd+XdjygczN4OrgxveJyP8TDJ8BY+YaHafDiYyM5MCBA0bHEAZRSuVrrePqWxbg6TDQ+IHeuSytnufahatmXcV7699Da33MdesCIoZdCfbp/Prl0xyxn8Jp58nIbXfLysoiOztbLm8VjZLRXx42uudo9pfuZ+fhnUZH8UrVp6Vx+bP+3HprCrriiNFx2j2LxUJ0dDRpaWksWLDA6DjCSxlyRtGRRRyKgJfgtYGvMfNGGWjmyi/YxP/+8wK9fr0NteZhiHvO6EjtmslkYsaMGU2vKDo0KRQedtbwsyActhZvNTqK14q94Fbo+gtsep69oRPpfsqVRkcSokOTpicPM/cx0+/2fhzpK80qjTrtKR56L5K48xIp3r/D6DRCdGhSKAwwuudo1uxaY3QM7xbQiStuSSP5vGo6rX/U6DRCdGjS9GSAyh8q+fk/P3P45sN0Du1sdByvNf7CWxjfeyusewK97TJU/z8ZHUmIDknOKAwQNyoORsGa7XJW0aSRj7Gq6GTOip/Kri0yol0II0ihMMA1l18DF8PWsq1GR/F+/kGEjHuaosNVFGbfDgYMEBWio5NCYYBhXYcRoALI/y3f6Cg+Yfi4KazLfpaYzitg83yj4wjR4bSoUCilblFK/Ukp1cVdgTqCQP9AQt4M4fVHXjc6is/wH34n1T3ieeZvd7H+x8VGxxGiQ2npGUUk8DCOO7r+qpR6SQpH64yOH0358HKMuNeWT1J+7B/6DHM+rOK1p6dB1VGjEwnRYbSoUGitn3beNCoSuBU4iKNwHHAWjgekaDTPNf93DYdGHOL3Q78bHcVndO8/krzP/8OcK/bCmoeNjiNEh9GqPgqt9UGt9TKt9UPOwjEWWAYMBrYopWQobRNiesdAGXy57kujo/iUAadPQw29nT0rnyF36QtGxxGiQ2hxoajvjME5IVGe1vpWIBq4Rik1qQ3ytVtDuwyFNPjP/P8YHcX3jPkH18wPI+mme6gollHbQrhbSzuz5wNblVL7nf0TdYtBNDimLdVaJwExbZiz3ekW0Y1eSb2oHlptdBTfExDKC+lv8tEDAQTm/Rm0/A6FcKeWnlFka62jcMxKp4AspVSVUqoK2A+glDrNua7MzNOE+KR4CkIKjI7hk0aMv5xRV74Auz5n+xePGB1HiHatpYXCrpT6k9baqrW+1Vk0orTW/lrrfzjXyVJKvYSjw1s0YnS30ezctJNNv28yOopvip7Owl/HE33BU+R+LtPLCuEuLb3qaRmwrG6Tk9b6oMs6g4GHtNbS+N6EkF0hkAFvLnnT6Ci+SSkuvP0d7r8iglMOzIaj+41OJES71OLObOcVT180tU7rI3UcCZMTIBHKepYZHcVnRXQfwJz0bDrpPZR/fT3VVZVGRxKi3ZFbeBioV1Qvhp07jA2lG4yO4tu6jqV46BNMuP1T0h6YYnQaIdodKRQGGx48nBWfrTA6hs8LP+0+Rg03M1h/DruXGx1HiHZFCoXBKn+qZP+r+1m3ZZ3RUXya8vMj493VJEweAiuuRpcUGh1JiHZDCoXBbrz+RkiBrUe3Gh3F9wWGwznv8v6KIi6ecArlpYeNTiREuyCFwmDxp8ZDb1i9Z7XRUdoH00jKB92K3W6neMW9RqcRol2QQmGwiJAI+u7vy0fvf2R0lHZj6h3P8e2bd9Ntz39gyxtGxxHC50mh8AL+Vn/yX5NJjNqSf9zTlJnOJvmWm8j/QsapCHEipFB4gVtSb6Hq1iq2H9xudJT2wy+QQ6PS+XwtfPf2HVC62+hEQvgsKRRe4OLYiyEYvtv+ndFR2pXuJ43g57zl3Gk5Ct9eJZMdCdFKhhQKpdQMpVSCUipZKZXcwm3T3ZXLKKN7jiZoVRAv/1fuV9TWOvc/C8b/l5/yV/BY8hnoarnTrBAt5fFCoZSaC9i01lla6wwgWimV0IJtzW4NaIBA/0A6bezEjzk/Gh2lfRowlfe2ns3L769iz3d/NzqNED7HiDOKZK11Vp3HC4GUpjZSSrXr+S1uefoWDl91mCPlR4yO0i499q8vWf3aFHpumw07PjY6jhA+xaOFooGDvR2wNGPzOCC7TQN5kYlDJ1Klq8grzDM6Srvk5x9A94sWoiPH8MzDCaz+eqHRkYTwGZ4+o4gCilyec318HGfT1CK3JPIS4/qMgyXw4oIXjY7SfgV0wj76fzz7SRWvPn0LlPxudCIhfIKnC4WpoQVKqXqXOZ+3a63tzdlBYWEhSqnan9mzZ7c8pQG6d+5OyO4Q1vyyxugo7Vpk3xF8/9USnrlOw1cXQ0Wx0ZGE8HqeLhR2HGcVdbk+dpWktc5p7g769OmD1rr2x1cKBcDVz11N0RlFaK2NjtKu9Rt5AX4T3mPf9nU8eP1ouSeUEE3wdKEo4vizChNAfWcMzj6NZhcJX3dW/7PYX7qfX4t+NTpK+9f7fL4oT+bF939j1ZtJIMVZiAZ5tFBora04zirqiqLhYhAFJDjHXczAcXWU2fm43V0mO6LzCHgdnv3Ps0ZH6RCS7vw3vy65n9M7LYXVqUbHEcJrBRiwz0VKqYQ6l8jGA7WD6JwFIMY5ziKHOkXEOTjPrLVO82hiDxkXPQ5/7c/GvRuNjtJh9LM8DXllfPL202xauIF7npKbMwrhyuPjKLTWNWcFFueBv8BlXEUC9YyrcK6byB9nFCaPBPagAP8Azv/b+ew2y32JPEYpiHuet9cM4I1FH1O+foHRiYTwOqq9dZzGxcXpvDzfHYsw99u5PJTzEDvu3UHfiL5Gx+kwyksPc3TZZYQfWg5nLYL+VxkdSQiPUkrla63j6lsmNwX0MiODR8Kz8MS/nzA6SocSFNqZ8As+otJ0OrfclMgnr/8/oyMJ4TWkUHiZ+NPiCRochK3KZnSUjiewM6XjFrF6eyg/Lf0b7PrC6ERCeAUpFF4mKDCIyx66jF/CfpHxFAYI79qPb/M28dD1w2H5pVTt/MroSEIYTgqFF7IMsrB9z3ZW/bbK6CgdUkhEX5iUTUFxL0bGTeKbj+S2KqJjk0LhhUaGjIS5kPZSu7wK2DeE9qLT5HfpZgohYuMM2PON0YmEMIwUCi90xogzMF1iYm/XvUZH6dB6m0/j69zNjB4+AL66iMI1WU1vJEQ7JIXCC/n5+ZEwPYH86nyqqquMjtOhqU59YPKXvJNnYvC4RHI/ec7oSEJ4nBQKLzVp4CQObj3IJ/mfGB1FhPbCctcybp/Sg1OLHoDt7xmdSAiPkkLhpU7rchqkwwsLXjA6igC69RvGP97eSFCPcRzJSeSjjHuMjiSEx0ih8FLDBw5n0G2DOHKKTI3qNYJMcN5npH05kD/d9hwFnz5kdCIhPEIKhRdLuDKBvEN5FB+VyXW8RmBnHp6/ms/+OZHoorlgfQB0tdGphHArKRRe7BLzJZSvKefZRXLbcW8S3CmcyXcvg6F38P3H/+SWy4ZytESKuWi/pFB4sTMHnInfp368+t9XjY4iXCk/iH2eH0sv4eu8Ag4tvRCONjn9uxA+SQqFFwsMCCTpn0nsnbyXssoyo+MIV0px91OLWf3Ff+lWkU/1Z2dQuOk7o1MJ0eakUHi5GyfdyJHKIyyzLTM6imhAp+E3waQc/v7mNkbHnc2OVZlGRxKiTUmh8HKTBk0iODeYWX+bZXQU0Zge53DNQx9x92Vd6bv+GtgsEyCJ9kMKhZcL8g+iV3Ev1lrXyihtLzdkTDx/fWUzqnc825cm8/DNsZSXHjY6lhAnTAqFD5j7r7mUJ5Xz7bZvjY4imhIUARMW8+GOSbz4jpVtCydA6U6jUwlxQqRQ+IBLhl1CsH8w762XW0f4BD9/7pi7jI1fv8Tg0A2wNIbf8uWGgsJ3SaHwAZ2DOjN402Bemv4S1dUyuMtX9B57K1ywko/yFYNPT2TZK8kyOE/4JCkUPuK8UedR0b2C7PXZRkcRLWEayYR7fiD1uuGcE7AAll8GR/cbnUqIFpFC4SOeuvspwpPCeWfzO0ZHES0U0b0/f391HUHjX6R02+dcdEZflr8/z+hYQjSbFAofERYUxtUjr+adb99hV9Euo+OIllIKhv6FXaPeZeteTekP98Lqh6Cq3OhkQjRJCoUPmRA6gbJ/lDHzuZlGRxGtNGjMpaz9dQ8X/mk6/DKXdx87mXXff2h0LCEaJYXCh1wz6Rp6XNWDnzv/bHQUcQICQiLg9AzKT1/IA//5jZl3XAm/PA0yTkZ4KSkUPsTPz48H7n2AvJI8Nu7baHQccYKCopNYmfsT6Y9eCKtnsP+98az/cbHRsYQ4jhQKH3PDqTfgt82P1OdTjY4i2kCP/qfQ+4olcMYbPJqxlnETLuXAikehusLoaELUUlprz+9UqRmADYgC0FpnNLKuCUh2PhwLZDe2flxcnM7Ly2u7sF6oV2wv9m3fR+nOUgL9A42OI9rI7q0/8/X/kkmM/h5Mo/mt1ywGxPzJ6Fiig1BK5Wut4+pb5vEzCqXUXMCmtc5yHvCjlVIJjWwyU2ud5vxJBFKVUsmNrN/uPfmPJ6maXkXWLzLatz3pOXAkiX/9Ds79gPz1u4geexXv/H2yzHMhDGdE01Oy1rruEW4hkFLfis6zCbPL0+lAh253uWniTZzc52TSvkuTkdrtUb/LGTZtNQ9PO52Len8Fi4ex78dnqKqQS2mFMTxaKJRSMfU8bQcsjWxmUUrVLRZ2ji8eHYqf8mPaoGmsfmI1c1+fa3Qc4QadI3vz+IIfiPiTFR1+Mok33c/F47vB3hVGRxMdkKfPKKIA1/PoBs+rtdZ2rXWk1tpW5+l4IMcd4XxJyoQUAv0DeWe1jNRu1yJPBctyUm67kxvPDYDss9HfJLFjoxQM4TmeLhSmhhY4m5ka5VzHQiNNT4WFhSilan9mz57d8pQ+oEtYF5544wl+ivyJ/MJ8o+MIN1J+flx95/Nc+/R2GDmLjz76kOiRZ/P9q9fIfaOER3i6UNhxXulUh+vjxiwAErXW1oZW6NOnD1rr2p/2WigAkmOTCQ8M56HXHzI6ivCEgDAYPZuxt37H/deOZqz/O/CRmXUf3MnhA3JbF+E+ni4URRx/VmECRzNTYxs6L6lN11p3+GanGhEhEcTYYsh5OIfsXLmrbEfRJzqWJ19bQ8Cla6nqNpErb3uRKyf2hw3PQmWJ0fFEO+TRQuE8E7C7PB1FE30OzstnrTVFQinVWOd3h/KvR/5FSEIIL/76otFRhKeZRuJ/3oe89koGj908Cqz3cfRdMy/NSqSkeJ/R6UQ7YsTlsYtcxk3E47jkFQCllLnucmdRiALylFIm5xVQ9V091SGdMugU/nrnX/no149YvnW50XGEAc64aDrn3JMPluUs/qUbtz+exQ/PDIZ1T0H5QaPjiXbAyJHZVpyXudYdae1cFq+1jnd2Xh+o5yWynIPvjtMRRma7Kq0o5aSUk6i0VrIvfx8B/gFGRxIGyvt8AbEB76J2fcbz2SHs9o/hb88twi+sr9HRhBfzqpHZAM5R1jla6wzX23E4l8U7/9+utVb1/NRbJDqq0MBQkoYkcXDvQeYvn290HGGwuPOnoyZ9Chfmsf5AH9au+g6/j83wwzR2b/rK6HjCBxlyRuFOHfGMAqCyqpLx/xnPriO72HTnJjoFdjI6kvASlQc2ErD5BYp+epmT/lLG/7txMA/89R/QZwr4+RsdT3gJrzujEG0vwD+AeRfN4/f9vzP1b1ONjiO8SEDkMBj7IoGXr+ex2y/kwhGH4Osr2JzRn6fuvRD7rl+Njii8nBSKduTs/mcTtyWOxX9fzNtfv210HOFlwrsNJPWZpYy8cwecnclnv4Tx2AufUfbBCPj2ag5seJ/qqkqjYwovJE1P7czuA7sZM2sMIdEh/HTbT3QO6mx0JOHFdm78mt6H3oWt/yPpHwfYWhTEyqxHUNE3QtgAo+MJD5Kmpw6kZ2RPFt6/kK32rdyZdafRcYSX6z3sXIh7Dq4sJOH6u5h28UDUz7Pgw4HcfllfPnn5XqgoNjqmMJhcR9kOnTPgHBI6JfDqja8ysGQgs26ZZXQk4e38Q0i64znH/x/eyoE1GXye+zRDI+Zxcef5VPS8hE8KhnDB1amEhJkMjSo8T5qe2qmi4iKGXTaMsnPKyL07l5O7nWx0JOFjdHU1FTu/JWhnJp+8/zqXPFnMxzNCmXL5VRzpdin+fc+XotGOSNNTBxTVJYr8D/MJ7RLKZW9dRmFRodGRhI9Rfn4E9T0X4l4g/rGdfPrWE5x/6VQoXMLLf59K926R7P7ocvhtIVpGgLdrckbRzn299WsmXjaRrqorv+f9TlBgkNGRhK+rKmfl0pdY8v7rPH7JDijbw71v+PHboW68u+BRVL8p0HmQ0SlFC8kZRQd27sBzufbSa9nXax+3f3I71VqmThUnyD+I06fczeMv58MVhWD5mj5Dz6F/ZDnKehd8ZCblwkiee3AS7MqBqqNGJxYnSDqzO4A3nniDgV8M5IlvnqDyUCWvXPsKfn7yHUG0AT9/6HEOD877yvG4eBPVOxaz4/kn6bZ5OXzxJdovhP97OYprEy/hoqQ7IWIkKGVobNEycrToIP523t+47eTbeO221zjj+jNob02Owkt0GYrfiPtY8sM+nlh0ECZ8zN6oa1n5y15+z10An4zmwP96co1lAPnvPwzFm0D+LXo9KRQdhFKKFxNfZPwV4/kx8kduW3IbldUyCle4UWBn6DuFHue/zKbfy5n2zBY4/RVsFbF8vWoHJavmwOJhrPpnd6ZO7k9B9iw4sBqqq4xOLlxIoehA/Pz8+O7173joyodIz08n5i8x7Ni7w+hYooPwCx8I0TcT++el7Nhbwdn3b4Cx8/m9eiTf/1RI582Pw9IxLLq/M+fHdWX/1zOg8FP00SKjo3d4Uig6GKUUcyxzeHLMk6xdsJaYm2LYat9qdCzRwSg/P1TEMBiSwpQZX7FtTwU9b9gCZ7xBVdezKSktI3L7P+Cri3h0aldGDQym6pvrYOML7Px5MSXF+41+Cx2KXB7bgb3w3gs8+sujqADFc5Oe48ZxNxodSYg/VByC/bm8/fp8cvOsPHP1ESjbxVXz4OcdsPE/oyAqls/XhWLqeyrjJl8DgV2MTu2zGrs8VgpFB2c7YOPqd64m9/FcTjn9FL5f+D3hweFGxxLieFpD6e/kfPAy+3esYerYMjhgZfQ9u+kXBZ/MAMIG8dCiAEaMGMH/3XADmEahw8womfWxSVIoRKNKj5Zy/i3n823Vt5w09iSeveBZrjz5SrmEVviEnbbVHPo9l6FRe9FFqxk77SPOGXyUZ29w1JYBd8MtF/XisdsmQ8QIluQe4bTTz6fvsDPBL9Do+F5DCoVolu+2f8eti29l7dK1dN/Wnc8++IwxA8cYHUuIFtMVR1CHNlC+ZxUPPzGfc4ZVc/nIfRzctx3TdHjqaki9LICSQDPXPlfMnddNYHL8+VSGmtlf3o0e/UegOtgXJSkUotkqqyu59uFree/D99BXa/4c82ceGPcAQ3sONTqaECesoqSI1Ss+pldoESd13s3WjVYueXA5f0+s4srYKtbtgJGp8PZdIVx94clsP9KHfy3Zzy3XX8rgU86gMuQk6HQSAUEhRr+VNieFQrTYzkM7mfPtHOZ/N5+K5yuYcOMEMmZlMLSrFAzRDlVXwpHf2G37kYWLsrjs9HAGhu/hixU/ceHs31n+VzhjCHy6Bqb8A757qjfjYk7m510RvLvCzm03XUmPAaMo8++FX9hJBIX63oRhUihEq622reaa269hc//NVPatxNLDwnkh5/HA9Q/IDQZFh1BVUQ6lO/Av2cr6NSt4M+sz7rmiJ90Cd/H64l+48QU7256Hk7pCxheQ8jLsyOhO35MGkv1LEIt/PMwT9yfQubuZncUhHKmOwDxiPH5BYUa/tWNIoRAnbPfh3byU9xJPP/U0JZ+V0PvR3kw7bxrXDL+GU3qfYnQ8IQxztKSYoMo9qJJt5K9czuLPvuaR6wYRcPR3nl+4hllv7mbvSxDgD3/NhCc/hKOvQUCIiReWhbAkv5yl8y5CderDivVlFB7wIzHhCgjpSRkRBIf38kh/iRQK0Wbsh+38c+E/yQvJ4/OCz6l+v5qwI2HMyJjBVSOuYkT3ESi54ZsQxyo/CKW/s27VCtausXL15H5QupN/vfUtn6z4jSWPREDpTqbNLyd7LWx/wbHZjfPh241QML8vhPTgHx+WsOeQP2n3XwQh3Vm++gA6MIKJ550Hwd0gtA8EdGpVRCkUwi12FO/g/qfuZ+XGlWwbvQ2NJvS9UEbFjGLGgzOYMHAC3Tp1MzqmEL5Baw7u/Y39Ozdh7hUIpTt5/+NsduzYwZ1X9oOyvdzxTC47dh/ig/s0VJVheRJKyuG72c7XiH0eht3Zqt1LoRBut/PQTj7a8BF/u+tv7Ivcx9HTj0I1hP4vlLMSzuKG627g9L6nEx0ZTYAMfhLixGgNlUfYYfuJ0oM7GdI3FI7ug65jIWJ4q15SCoXwqIqqCnILc/lk9SfMf2Q+R045QtnJZVAM/BtG/HkE8VPiGR4xnJADIVx81sV0j+hudGwhOjQpFMJQVdVVbNi3gSW5S3jtX6/hP86fgqACSjaVwGvADTAgZgB9Svtw+LvDXDntSsaeMpbewb05qctJ9IjsYfRbEKLd87pCoZSaAdiAKACtdUZbrS+FwjdUVVeRuzmX9z97nwBzAFvKtpD7RS6bX94MKTj+0quBD8D0oIlB0YMI/T2UYmsxl992OYP7Dsb/kD+BZYGcPfZsenXpRYCfNGkJ0VpeVSiUUnOBXK11Vn2PT3R9KRS+rbq6mr0le9lq38ryH5ez7NNl9L+wP4VlhaxevJqdH+xE36shCPga+AJ4BAiE0NxQqvKqiPt7HN06d+PwmsMcth3m8jsuJzIkkuJtxVQeqmTS5ElEhERAGXQJ6UKfbn3wUx3rdg1CuPK2QnFAax1Z53EMMFdrHd8W60uhaP9KK0opPFRI/rp88lfnM2D8APYc2cP3n37Pxu83Yp5mZn/pfmxZNg7nHob7nBt+DGwAHnQ+/hDYDNwPYYFh6KUavUdz8n0nExYUxp5P91BVXMWZyWfSKbATW77cAhVwbsK5hAaGUrCygEC/QM6YdAYhASFsW7+NkKAQRp06imD/YPbv3E9YSBj9+/UnOCCYirIKwkLC6BzaWQqT8DqNFQqPnqs7D/Ku7IClLdYXHUNoYCjRUdFEnxNN0jlJfyyY6LLirY7/lFWWcaD0ABuu2MC2ndvoObgnB8sO8kP3HyjcUcjJ557MofJDrNy+kqLCIvp26cuR8iMcLDpI6f5Svtn2DUfKj1C0tIiqw1V83vVzxwu/ClTDcweeczx+GccnqmZajwygE3C98/FLgAm4BvyVP9Xp1QT0DCDimggC/QIpeqmIkJNC6J3Qm0C/QH779290GdyFgZcNJMAvgPX/Xk/Xk7sSfWE0AX4BrJq/ip4jezJk0hAC/AJYmb6Sfqf1Y+jZQ/FX/qx4eQUDxwxk8OmD8dN+fPvGt0THRBM9Jhqq4LvM7xgcM5hBpwxCV2p++OgHhowZwoBhA6g8WknuZ7kMPW0o/cz9KC8tZ/Xy1QwZPYTe/XtztOQoa39Yy5CRQ+jZpydlR8pYb13PkFOG0K1HN8qOlLFp7SYGDx9MZNdISg6VYNtowzzUTIQpgtLDpfxW8BuDhgwiPDyc0iOl/P7b7wwwD6BTWCdKj5Sya8cuThp4Ep06deLI4SPs272Pfv37ERIcQsmREvbv20/vPr0JDg6m9EgpxfZievbuSWBAIGWlZRwqPkT3Ht0J8A+grKyMksMlRHWNwt/fn6NlRykrLSMyMhI/Pz/Kj5ZTXl5OeHg4fsqPiooKKioqCOsUhp+fH5UVlVRXVxMaEgo4znx1tSYwMBClFLpao7UmICAAxR/jiNrLmCJPN+pGAa7zGjY2z2FL1xfiOCEBIfQO703vU3vDqX88P3Xk1GNXvMBlQ9d5nB4ErTXlVeWUVpZSeFMhh0sPE941nLLKMtaesZbyynL6DulLWWUZP/b8Ee2viY6JpryqnG8qviGwUyCDzxxMRXUF32z7hmBTMINHDKaiqoIVI1bQuXdnBvYYSEV1BUURRYR1DiPQL5DK6kqOlhyl+EgxO4p3UFldyS7bLg6FH2LPtj1UVldSuLKQbXob+aZ8KqsqKV5czJqDa/DTflSWV1L9ejXf7PjG8VXrKI7CFQ+cBZQA/wAuBMYDh4B/ApcAY3FsMw+4HBgD7ANeBP4EjAZ2O18vETgFKMRRKK8GTga24yik1wFDgC04LmS4ERiE48zuDWAa0B/YCLwNTAf6Ar8Ai3AU/17AWuBd4C9Ad2AVjjPEu4FIIA9YjONssguwEljq+BsSBnwHfA48BIQA3wDLgEdxHBW/cv7MAhSQ49zmMee/hc+c+3jE+fgT4Cfn68HxZ68fON6zulehlKL63WoohIC7HIfgqkVV6H2awNsdhafy7Ur0IU1wSjAA5W+UQxkE3+J4fPR15+XnN4eilKL0lVJUuOLfC/7N9NjptDVPFwpTQwuUUiattf0E16ewsPCYKj5r1ixmz57dwphC1E8pRXBAMMEBwZgGmI5ZNqb3sbdkv3L4lcc8/su4vxz7YpNcXvwyl8cudYybXR7f5vL4XpfHM10eP+64iKBaV1NZXcmhBw/hH+CPf6A/lVWV7L1tL8GdggkJDaGisoLCmwrp3KUznTp34mj5UbYmbaVr96507tKZktIStlyxhZ59ehIWEUZpSSmbLtpEnwF96GLqwuFDh9lw3gYGDh1IeGQ4xQeL2XDOBqJHRBMeGc7BooOsP3M9Q0YPITwynAP7DvDLuF8YETeCzhGd2bd7HxtiNzBy/Eg6denE3jP2snHMRkadNcrxeOxeNozewJjzxhDaOZSdp+5k46iNxMXHERwWzO/Df+fXkb8y7qJxBIUGsWPwDjafspnxU8YTGBzItoHbKBhRwJkXn4l/gD+/9fuNLSO2cFb8WfgF+LGl5xa2jtjKuZZzASiIKmD7yO2cO+lcNJpfw39l5+idnH3e2Wit+TX0V3bF7OKsCWcBsDFwI/u27+PMCWeitWaD2kDRziLGnzMegHWV6zi09xDjzhyHRvNz6c+U2EsYe8ZYNJq1xWspO1xG3DhHS9Ca/WuoLK8kZlwMWmtW71qNRnNqrOObj3WbleCwYEZ0H4E7eLSPQillATJd+hzMQAEQ6Xrgb+n6IH0UQgjRGo31UXi6R62I488STAD1HfRbsb4QQog25tFCobW24mjtrCsKRwvgCa8vhBCi7Rlxjd4ipVRCncfxQHrNA6WU2WV5o+sLIYRwL48XCq11CmBWSlmUUslAgcvguQQcY3Obu74QQgg3MmTUj9Y6TWudo7XOcL0dh3NZfHPXb0tydZR3kr+L95G/iXdy199FbgpYh1KK9vb7aA/k7+J95G/inU7k7+JNVz0JIYTwMVIohBBCNKrdNT0ppfYCv7Vy8z44bj4gvIv8XbyP/E2804n8XQZoreudQazdFQohhBBtS5qehBBCNEoKhRBCiEZJoRBCCNEomWSYls/hLdzPOQo/Fsh0PpWIY2ZDm3GpOhallAlIBrpqrVPrWS6fGwM09ndx1+emwxeK+ubkVkolyG1CvEISjg+EFZguRcJznLf4NwHRDSyXz40Bmvq7OLX556bDFwog2aUqLwTmAvIP3mB15yERnqW1zgFQSo2l/gnE5HNjgGb8XdzyuenQfRQyJ7cQLSefm46no59RyJzcXszZ3lqEtIF7G/nceDF3fG46eqEwNbSgoTm5hcfkAfaa9lWlVKZSqkjawL2CqaEF8rkxnFs+Nx266QnH6XKUy3Ouj4UBtNZWl064XGCmUXnEMezI58Yruetz09ELhczJ7aWcV3fUZQPqaxsXniefGy/lrs9Nhy4UMie3d1JKmYFs5/XidcnlsV5APjfeyZ2fmw5dKJxkTm4v4zx1TnX5djoVx+WXwjvI58bLuPNzI3ePpXaEqRUwg1xd4w2c345qDkRdccyVLn8XD3FeAmvhj/nr04Ec59lEzTryufGwpv4u7vrcSKEQQgjRKGl6EkII0SgpFEIIIRolhUIIIUSjpFAIIYRolBQKIYQQjZJCIYQQolFSKIQQQjRKCoUQHuIcDCWEz5FCIYTnZLrc9kIInyAjs4XwEKWUBiLlDqvC18gZhRAe4Lz9s02KhPBFckYhhBs5C0Q8jhu12XHcijtXZuoTvkQKhRAeoJTKB+ZIgRC+SAqFEB4g/RPCl0kfhRBuJv0TwtdJoRDC/eKRaUKFD5NCIYT7WYB8AKWUyXmGIYTPkEIhhPuZ+OOMIllrLWcXwqdIZ7YQbuYcjW0GbDjmN7Ybm0iIlpFCIYQQolHS9CSEEKJRUiiEEEI0SgqFEEKIRkmhEEII0SgpFEIIIRolhUIIIUSjpFAIIYRolBQKIYQQjZJCIYQQolH/H+den4171nLgAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"y0 = -1.0\n",
"__=ax.plot(t_domain, analytical_dy_dy0(y0, t_domain), color='green')\n",
"__=ax.plot(t_domain, ode_result_jax_grad_vmap(y0, t_domain_vmap), ':', color='k')\n",
"\n",
"y0 = -5.0\n",
"__=ax.plot(t_domain, analytical_dy_dy0(y0, t_domain), color='orange')\n",
"__=ax.plot(t_domain, ode_result_jax_grad_vmap(y0, t_domain_vmap), ':', color='k')\n",
"\n",
"leg=ax.legend(handles=[green_line, orange_line, solid_line, jax_line])\n",
"\n",
"xlabel = ax.set_xlabel(r'$t$')\n",
"ylabel = ax.set_ylabel(r'$y$')\n",
"title = ax.set_title(r'${\\rm JAX\\ solution\\ for}\\ \\partial y(t)/ \\partial y_0$')\n"
]
},
{
"cell_type": "markdown",
"id": "provincial-roulette",
"metadata": {},
"source": [
"## Calculating $\\partial y(t) / \\partial y_0$ by differentiating through analytical solution\n",
"\n",
"Since $\\partial y(t) / \\partial y_0$ admits an analytical solution, we can also use `grad` to operate on the analytical implementation. This last calculation is really just for fun: the magic of autodiff is that we will get the same result within working precision no matter how we do the calculation."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "boxed-instrument",
"metadata": {},
"outputs": [],
"source": [
"analytical_ode_result_jax_grad = grad(analytical_ode_result_jax, argnums=0)\n",
"analytical_ode_result_jax_grad_vmap = vmap(analytical_ode_result_jax_grad, in_axes=[None, 0])"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "breeding-single",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEiCAYAAADnMZWTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABJDElEQVR4nO3deXhTVf7H8fdJN8qaFhAoIJC6IajYFkFRAUlVUBC1BdcRF1r3wY0O6m9Axxls3bcZUtTRcYXWFfeWEdwGpA244kLDvkMattJSmvP7I2ktJelGk5u039fz5NHk3uR+0pD7zT3n3nOU1hohhBDCH5PRAYQQQoQ2KRRCCCHqJYVCCCFEvaRQCCGEqJcUCiGEEPWSQiGEaBallPlIlvtY33IkeUTgSKFoQ5RSVqVUtlKqVClVrJRKqmddcz2vYVNKaaVUSd0vd61lefW9/pFQSiUppQqUUgWBeP1gUUpZlFLTvbeA/K1qbcvs/ezSWmKHrJTK0Fq7GlhtUmOLhXe9LH/LWjK7aDol11G0Pd4dbJ7WOtfP8jQgU2udWs9r2IAUrXVynceTvI/7fO2W4t1Odn0Z/TzP7GsHp5TKBsxa68wWitiYLMXAGCADQGudE6DtWL3/W+T97yQgubnvVSmVARRqrR11Hrf4eCxba+2zANRZbzpg11oXBjK7aB45ohC+WABrfb8GvV9Ui3enUZs10EXCy9XM503y8/g8wNbM12wy79/WrLV2aa1zAlgkzEAB4PBuy1X9+Xh/EDTn9RLrFgSvPB+vWdDI7aT6KBItml00nxQKcQjvl9PuvfnbqVZLB2zVBcX75Q1GkTgSPo9AtNZ2rbU9iDniAV872xblPXpK97FjLwaa86t8Ep6i6ksScMjO3rvzn1zfC3qbkw772wcgu2gmKRSiLqv3yz2bBr6M3vXygbnVbceNaLc2THXzktE5gk1rne9nUXPa+1N9FVRvE5HDz+fvbKBvIRM/R3MtnF00U6TRAURo0lrnezukD2t3rmMqsBrIak67sbfpqvr1zUB8reaFDMDpXRbvr0nL218xF3BW91l4n5uFpx8j13u0Y8bTXDbd+9RcrbXLuxOzed93aq3X9bn96v4RPM1fs/EcHZiBoY1sj7fiObKpzuKo3iE2sM25eNrq8/DsKFO11umN3F4Sf/ydHdQqmN7lZu/rZdZ6fDpQu7mn5jk+3ksa4PIW42V1dvAFgBX/R5tJ/v6NNZS9KfnFEdBay62N3fB8cTN8PG4BLLXu2/DsaBt6vWygFE+be1NyZNTZnrl6e3h3hnWWFdTJWvt+Uu37tfJn1LeOv2WN2L4VTxOIpc5zrI187xYfeRu9Te+y6Y3Yjq3uenh26gXe7ZmBNO/jh3yG3m2l1fkb+fz3UHddH39bf8+z+nsfDWWv9TdqVH65Nf8mTU+itrpHDza8Z+T44/2VuQzPL925zdhmzS9A7Wm2mOf99XxIFu8yh4/O8/q4mrNuI7fvxLNTqv33ctDMJpGmblN7O8EbeM0MPGeg1V3PjmcHXeBdnl99eq53m9Wfa90+h3hgp5/NHdY/UYsL/3+XdHwcaTQyO43N7z2lOq365ieL8EOankRtqb7O51dKVfdb+JKhtc5RStmBkgbWPYT2NAnlKaU0ni90nvex2s1RtZUAyT4eb2kpjdy+y4BtNqUD3IbvfqbqnfZ8/UefQiYwv24efWifgwUf77mB/gnwFjg/y+L9PK/B7FDTTwYN55+hvc103mt9CuvJK+qQQiGAml+0s+t+eZRSXfF8CQ/b+Xt3EPkAWmuHUioLT3NGXCO3adZap3t//aUAWUqpZDw7R3/iG/Pajdy+v/4XcyO37/S7VtM1dpuuxrxYrYLvq2hn4u2fqfXYJDy/7qul+niuE99/f1/r1ubzDC/vL/vDzqBqRnaoJ793O7W3X4LnqMRfR7moQ5qe2gDvFdTWWg/5+rJb/PzCmoenXdiXuk0lOXiaShp7PcIM7/NcWutC7elItuD5gvtqqkjkjyaHulwc/r7MDWzf39XQzdn+kQrINusWQm/HvZVaV0HXuqaj9o65dvNONZc3T13VfSc1V1HXWW7Gd4GbrP2f1dSo7I3Mb+HQJjMXMNTfdsXhpFC0DXa8V7Z6v1QWah2mex/z+Utde06FdNRt11VK2bTvs0nSgYzGDknho8/BUWubSbXWM1P/Fd++mjcsdR6r24fg8vVCTdi+r79Z3QyN0oRtNur163m9PGBM7R8FPo4ik/Dd51CE72JmrrVuho+mRwuefqza2/D7PpqSvYn562YWjSRNT21DFp4rrV14DslrvmzeX3/ZgFkpVaTrnCPv3ZHHA9neX3QuPIf/SUqpPB87hepfkwuVUrPx3UxQbSeeHUL1c8zeLHibpKbXOv/egme4i+pfltlAilJquvZc2ezytj1X92+Y8fyizFRKObTW+XXWceo/Tkmt/XoZWuvcBrafhOdoyFK9fW8hrT5F9LChKOr8TQ/JX/03auQ2D3mOv214JeP53Kqb8hLxfQEbeJr9svmjT+SwPgfv38/s67lAmlLKge9TYFM5fBynSdR/JXxTsjcmf9c667sQjSZjPQkhDuFtOizxdVaVt0gV1v1B0cDr5ek613sopQp0E8fpasL2DsnvLeI117jUvXZFNEyanoRo47x9WEne/zfj+bXvr4kvlwaG5Kjz2lbqdFh7j5habPiShvJ7C0LtJrNE6m+WEnVIoRCiDfPuWGsfIcwFpvpr1vI+ftjw8vVI9fHLPY0WGoCxCflt3msorHguYqy7XNRDmp6EaOPUH8OGWPAx1Lef50xv5AV/8+vulFu62ak5+UXTSKEQQjSL8jO3R0PL67l+RYQoKRRCCCHq1epOj+3WrZvu37+/0TGEECKsFBcX79Bad/e1rNUViv79+1NUVNTwikIIIWoopdb6WyZnPQkhhKiXFAohhBD1kkIhhBCiXkHvo/BeIJMBdNWNmzZyOp6rOOPBM4dBQAMKIYQ4RFCPKLxXRVrxXEJvbsT62XjHZPEWiESZnUoIIYIrqIXCO+dAPo0fuTGjzuX/8/A965UQQogACdk+Cj/zGbj4YxhrIYQQQRDK11HEc/hUky059eQhiov/ge2h2axef5CIiAi2ODVuHcPEkQN48MbTeeMzB+8sXs3gY3vRqVNnNjvLiYiKZeKF5zJsxDlUqDiiO3RDmUK29gohRLOEcqEw+1tQ3xgzmzZtQilVc3/mzJnMmjWrwY3tK/0Z19a9/OydJsVVBgeqymhv2sWDY9fwyWe7+WjpQfIKfjnkeYWfvMuyv8Glj8BXv0FiQiyJfbtiijYzaPDJXJtxO32OPQ1qZRJCiHASyoXCxeFTTfqcrrO2hIQENm3a1OSNnW19lbOtr+J2u1m7ZS3T7p/GB/M/YKnZzfjKM3j9s9d5OaYTBw+Us9e1lZLffmTrpjUkxEdCYkcuu+wDqt77Gn2wDPvKLZRs2QCf/kjxotd5974EHvywK3voya13zaTfwBFNzieEEEYxZFBA79lMZq21345pbx9FsdZa1fdYXSkpKbqlhvDYvW83T379JDPnzOQkfRLfv/N9o5+7ZX0JSz/PJ6H9TpKOWkv/S/LYsNPzt05KjKVHrz7cdc99jJlwTYtkFUKII6GUKtZap/haFrIN6t6JSFx1Ho4niDNTde7Qmb+e+1dO3HsiPyz4gb+/8PdGP7dn30Qu+lMWQ9NyiDh7Hmu3VrLiy3d45O4LcKP4+KvfsV40hTfvO5HK315i/56Adb8IIcQRCalCoZSy1LlOYn6d+6m00MxYTbHw9YVEdY1i5j0z2bB9Q7NewxQRwSlnTuTuRz5g+ap9LPvyQx64xYp14H7eeuZaevToyh1Xj8C1ze+4XEIIYYhgX3CX5L3SOg2wKqWm1zkNNo1a10l4m6YsSimrdxarEiMmRO/ZrSfP5T5HlauKYWOGtchrppw5jr8+W0C3K0s49sJ/0b9XR5589RuO7t+f6y46hXUlP7bIdoQQ4ohprVvVLTk5WQdK70G9NaBnPTUrIK+/4ot5etLovhrQURHoVx+5RuuqAwHZlhDhrKSkRKelpemCggKjo7QaQJH2s18NqaanUPfFx19AFLz/v/cD8vqnnDWJef9dx6tz/sZJlo5ox8uUv3cySz+ZG5DtCRGOCgsLcTgcOBwym2qwSKFoAktfC2n/SqPklBL2V+4P2HauzLyfol92ceWMd3jy7c0MH5vB1eNOYG/ploBtU4hwYbVasVqtxMc3eLa8aCGhfB1FSLo65Wryp+bz1w5/5ZHbHgnYdpTJBH0ncstTP/P1mhG89smvfHl8H+Y+9xip6X8O2HbFkZv2yTRWbFlhyLaH9BzCk+c/2aTn5Ofn43A4sFqtJCV5ugxTU1MpKCgIQEIRjqRQNNH5J56P+lFhe9QW0EJRrVN8Agu+XM03H9qYfM0tnDtpGpPPf5HXPyjGFCEfnzgyhYWFWK1WZs+eTVFREUlJSbhcriZPJ5yV1eCMAUyePLmmEIkw46/zIlxvgezMrjZs0jCNCf19yfcB31Zta1f9oAf06qCVQn/73DCty3cEdfui9SkpKdFaa52UlFTz/wUFBTotLU1rrXVeXp4uLi7WNpvNsIz+WK1W6cxuQdTTmS0/SZvhvtvvY8L8Cdz3xH28/0xgOrZ9OTpxMCUbdlP07v2kVDzKrCnHMm7qvzjtnMlByyBaF4vFgsPhwOVyYbFYACgoKCA1NZXCwkKcTidpaZ5LmXJycpg+fXqLZ8jNzaW4uLjedbKysmryieCTQtEM40eMJ7JzJB/YPoBngrttZTIx9JJ/sPnH03lq6kVkv3UZb73wK+Ou/mtwg4hWw263Y7X+MXp/YWEhmZmZ2Gw2UlNTATCbzRQUFPgtFEfS9JSRkdHM5CJYpFA0U9KZSXz70be8WfAml6VeFvTt9xo8noJP3mfqNemMv2YmT61aya0PvBH0HCL8OZ1OEhMTAXC5XDgcDiwWCy6Xq2ad+Ph4nE7/w8xkZ2cHOqYwkJwe20wvPPsCnAX2UrthGVLOvJCv7GsZfWpXbnvwTS4cdQLa7TYsjwhPkyZNoqSkhPz8fGbPnk1KimdcOLPZXFMsnE5nyJyOarfbycnJoaioiOzsbHJzc42O1OrJEUUzDR4wmBMnncjSTUsNzdGhy1G8v7iEk0/ow4eLf+Wzf6Zz3i35Mv+FaBKbzTOE2rJly8jM9Iyik5qaWnNRm8PhqGmGMlpSUhJJSUkB6S8RvskRxRGILIzki3u+YNeeXYbmaN+xC7+tLSXv4bGcG/c2y/6TLkcWolEcDgdjxowBPM1Odru9pvPaarXicrkoLCzEbrfLjrkNM2Q+ikBqyfkoGnLlHVfy+pOv89z857g5/eagbLNeWrP8zWtJvuJlUof15ZNv1sjUrKJB+fmecTYdDgcZGRmYzWZjAwlD1DcfhRSKI/DT6p8YbBnMeZnn8cmcT4KyzYZot5uJo4/j/S9KuG3yyTz1+nIpFkKIBoXlxEXhYNCAQUT3iOZ/Bf8zOkoNZTLx7ue/ceeVyTwz73smnnO80ZGEEGFOCsUR6tiuI7sduw3vp6hNmUw8+p9vGXFKD95fvIo7rhlpdCQhRBiTQnGExl82HoA3C980OMmhlMnEf5es5uxTurFi2Rcc/P0loyMJIcKUFIojdP/t98MNsCNuh9FRDhPdLpbPvinh48dG8ss711P45myjIwkhwpAUiiN0TMIxWLpZmDd/ntFRfIpp35l21gXc9HIsE6fcy/Iv5hsdSQgRZqRQtADT1yZ+yP2B3Xt3Gx3Ft6hOPGGbR8dYExPTrmDL6u+NTiSECCNSKFrA6HNGQxXM/2/o/lpPOfMCPnr7P+zYXcWIM05jx5Z1RkcSQoQJKRQt4M/XeGace/fTd40N0oCk0Vfy0F+uw7GlgtSzT5Grt4UQjSKFogUMGjAIUwcT//s0dK6n8OeOWS+QkT6U0lIXW7+aZXQcIUQYkEEBW0hUVBSla0uNjtEoc95cwq5PLsb960N8u7c7p427zehIQjRaTk4OO3fuZPLkyTidTvLy8moGNRSBIUcULWT4hOHoGM22XduMjtIgZTJhTn2Na+Z24JxLbmdl0WdGRxKiSXJzcxkzZgw2m03mwggCKRQt5K6su+BuWFm60ugojRPVkbtn2TApuPKyiZTvcxmdSIhGMZvNlJaWUlpaSl5engxiGATS9NRCkhKS4DN4MepFRt4XHkNmjBx3Ba/9cyUTrnuI6y4eyuuf/W50pNaheBqUrjBm23FDIPnJJj0lPz8fh8OB1Wqtmao0NTWVgoKCls/Xgux2O2azWebSDgIpFC2kt7k3arniI+dHcJ/RaRpv/LV/47J5+bzx6S/0uO58nngxNEbBFcFRWFiI1Wpl9uzZFBUVkZSUhMvloqkjMB/JnNnNkZ+fj9VqpbCwUJqfgkFr3apuycnJ2ig9knro6B7Rhm2/ufbv2637HRWju3dG79tUZHQcEUQlJSVaa62TkpJq/r+goECnpaVprbXOy8vTxcXF2mazGZaxIRaLRRcUFBgdI+wBRdrPftWQIwql1HTAAcR7i1W9k95613d575q11jkBDdhMA08eyCL7Ijbv3Eyvrr2MjtNo7dp34ptvvqZ0QSoxy66hYsxXxHQwGx1LBIHFYsHhcOByuWqacAoKCkhNTaWwsBCn01kz411OTk5AZrnLzc2luLi43nWysrJq8tnt9kOOTJKSkigoKMBqtbZ4NuER9EKhlMoGlmmt86vvK6XSqu/7WH967cKglEqq+1ioGDRgEItYxBPPP0FOVsjFq1dCYjLdL34Z69gJdO2ZTP5/S4yOJILEbrcfspMtLCwkMzMTm81WM0+22WymoKDAb6E4kqanjIyMJmUdM2YMpaV/nIrucrlITExs9GuIpjPiiCJDa137X9U8IBvwWSiAyUDNXldrbVdKzQhgvmYbP3o8z818jq+XfG10lGaJ6jcec49E3vq8hHefn8HEG2S02bbA6XTW7GhdLhcOhwOLxYLL5apZJz4+HqfT6fc1gtVHkJSUdNi2HA4HkyZNCsr226qgFgqllK+eLBdQ3zGjUymVp7VO975GBp7iEnLOO+s8oo+Jxt0rfIfGmPfRcs44uSfXTcsm+ayJ9D1+mNGRRIBNmjSJrKws8vPzWbZsGSkpntkwzWZzTbFwOp3Ex8cbmPIPKSkp5OTkYDabKSkpkVNkgyDYRxTxQN2fJf5/pnhkAsVKqVJgNuDw10wVClIfSGWNa43RMZotOrYTb8x/l5OHnUvSsDPZuGU30e1ijY4lAqz6yuZly5aRmZkJeE6RdTgcgOdXe3UzlNGSkpJa7Owp0TjBvuDO7G+BUsrnMq21A0+BKMLTRDW0vg1s2rQJpVTNbdasWc0O2xztfmvHTzN/Yqtza1C325KOPTWVyy86kx27DvL3OycYHUcEkMPhYMyYMYCn2clut9d0XlutVlwuF4WFhdjt9oB0ZIvwEOwjChfeM51qqfd4VillA2xa6xyllBXIU0pZqpui6kpISGDTpk0tErY5ulR0ASe8/M7LTL8+fL9Yz7+5CDUhkSGxn0PpdxB3itGRRABYLBZmzJhRc9FdXl7eIcuri4OcUdS2BbtQODn8qMIMoLV21V3Z26fh0lrbvesUKqUGAKsDmvIIZFyVwYsPv8iyX5YZHeWImCIieOGtIsrfHUz2nefy52d/pZ2cMtsqVR9BCOFPUJuevDt8V52H44FCP0+JB3bWeQ1XPesbbujAoaj2iu+/awWzyLXrxpLIO/nLS9u47IJko9MIIQxixKCA85VStX/CpAI1YwQrpSzVy7XWhd7l1FpuxnOxXkgymUxERUdR8mXruA5hVNp0rjj/ON5b7ODztx43Oo4QwgDKc+V2kDfqudLaDljg0CuzvctStdap3vsWPGc+1ex567uSOyUlRTd1nJqW1ufkPmz8eSOVFZVERoT/cFplu3cwZGACZeVuvrX/TEK/44yOJIRoYUqpYq11is9lRhSKQAqFQjHt2Wk8Ne8pfpr3EycmnGholpby8RuPMe6Kuxk0wMyPjvCYoEkI0Xj1FQqZjyIALpt4GYyBn3f8bHSUFjP28rsYd6aFlWtc/PblXKPjCCGCSApFACR2SoSH4ansp4yO0qLe/nQ53zxq4bgdM6ncF/oz+QkhWoYUigDo3qU7JpOJX1f8anSUFhXTvjOnXTOPPz2+hUnjfB6hCiFaISkUAdJ9YHd2O3cbHaPFqa4puDsP4r0v1vPByzONjiNClN1uJzk5uVGjygZDoPKE2vsMFCkUAZJ4QiIVWytw7mpoKKvw80Le1wzqF0PmnQ9RujVkr30UBkpKSqoZM6opcnMPP6ExPT3d5+PByBOo1w3U+wwUKRQBYo4wg4bn5z1vdJQWF9O+My+9+DxbXW5GnykX4gnfmjrabPUESnVlZma2yBAigRr9NtTeZyBIoQiQiydcDMCPjh8NThIYyedcxbCTevPdqlLm20JyehARZvzNaWG1Wmtmt2sNwvF9hv/VYCHqiguvIHNcJu1OaGd0lID5+PPljErpy94fbFA5A6I6Gx0pJEybNo0VK1YYsu0hQ4bw5JNPNuk5drsd8PzSLSgoIDs7G7PZTGFhIVlZWUyePPmQaVLrTkvq67l15efn1zyvev6I9PR0HA4Hc+fOxel04nA4cDgcmM1mLBYLVqsVu93O1KlTsVqtNTtYh8OBzWZj6NChmM1m4uPja4Ydb2wef3Jzc2smbSooKKgZft3lctUsczqdNfnqqv6bpaSkYLPZsNvtNe/bZrNRWFjY6Pfpb5uN+VxamhSKAGnfrj0DUwfyy/ZfjI4SMJ3julP0v8+hYAT/nXsV59z8vtGRRDOkp6eTl5dXMzhgVlYWNpsNq9VKZmYmeXl5FBQUAJ4JjGw2W83OzN9z66pePm/evJodd2ZmJikpKTX3q3fytadGre4DKCn5Y0ic1NRUiouLMZvNNTvY6jm3G5vHl5ycHKxWa03RqZ6LA2DMmDGHzOudnp5+SIGqVv03q163Ov+8efNqljf2ffrbZmM+l5YmhSKAdr25i5U/roR7jE4SOKajTufpFSP586MLsJXdRcbdjxkdyXBN/UVvtOqdLniGHa+9gwQO+UUeHx9fs3NqzHNrS0tLY+rUqbhcrprZ8xrzaz8+Pr5mB5qfn4/ZbK55XlJSEgsXLmxWnrosFgtTp04lMzOTSZMm1ezI8/PzD/ulPnnyZGbPnn3YsOxHou77bGib9X0uLU36KAKoT68+uPe5cWwM2TEMW8RV97xCT3ME2Y8+TdnuHUbHEc1QPRVqUVHRYXNjN9RZW99z65o0aRK5ubm4XK5mNZM4HI7D8tQtNk3JU1taWhozZswgLy+PuLi4mlNely1b5nObTSlCTdWYbQZzalopFAF0/rjzAfjom48MThJY8Uf1Yd5LT+DYepB7bzzX6DiiCVwuF8nJycyYMYO0tLSa+bKrl7X0c6ubggoLC+udzjQ/3/dsx0lJSX53/kfyXsDTv5CWlkZBQQFaa4qKinA4HCQmJh62zaYUuvqKlb/3eaTbbGlSKALovNPPA+Cjz1p3oQA4+6LbuPQcC0+9sZwXnrrP6DiikYqKig5pyqk9R3ZDv5gbeq7T6TxsB22xWDCbzT53nhaLhZ07dx72eO3Xqe5Arp2t+tqD5uSpraCg4JDXrd5WRkZGTb9CtXnz5jFjxozD8lW/j9rvr26TUGPeZ0PbDDqtdau6JScn61Bx8OBBDei4vnFGRwmK9at/1u2i0cf0itIH9u8xOo5opIyMDG2z2XRBQYEuKSnRGRkZOjs7WxcXF2ur1aotFosuKCg45H5eXl6TnlubzWbTpaWlPrOkpaXVvJ7W2ufrlJaW6unTp+u8vDydl5enS0pKmvxefLHZbDXPrX7taiUlJTXbzM7OrjdfdY7q17DZbNpsNmubzdak99nYbfr6XJoDKNJ+9qsyzHiAdejTAVOUiT2r9xgdJSgK3vg7u5fez6XXzYSTZxkdR4Sg/Px8mX41BMkw4wYaO30spnQTra0g+5N6+X2MnXAZWfc/yHdfttwZISK8ZWZm1lxDUF/fhAhNUigCbOjAoexeu5tfNrTe6ynqKjtuFv8s0Fw86WqqKg8YHUeEgPT0dFwuF3a7PWSvPhb+SaEIMPc6N8wH2yuNu+inNejW53gen3kDq7dU8PT/TTI6jggBVquVtLQ0aXIKU1IoAuzc0z2ni67asMrgJMF1w3QbF55+FPc+/h5ffjrP6DhCiCMghSLAkk9MxtTRxG+//WZ0lKBSJhNP2PKpqIRLJ1+NdruNjiSEaCYpFEHQoVsH1n23zugYQXfMSWdx69Wj2bG7km/f+avRcYQQzSSFIghiVAwVOyo40AY7dp/892esmJPEMPezVJTKJEdChCMpFEEw6qJRAHz7y7fGBjGAKSKSky97g2uf28sF1tOgjZwmLERrIoUiCG7MvBFmgCvGZXQUY3Q+jnY9Ulho30HenLuNTiOEaCIpFEGQ0j8FVsBbC94yOophnnllEUnHxHLbvU/i3FzS4PpCiNAhhSIIurTrgvpcscC2wOgohomMbseLL77Ejt1uzjl7qNFxhBBNIIUiSOKPiWff/n1GxzDUKWdNYvjJffhuVSl5NhlhVohwYUihUEpNV0qlKaUylFIZjVjfrJTKrl5fKRV2g8VYBloo31xO+YFyo6MY6qP/2jnVEsO+n+ZAZdsYKDEcFBYWkpiYaHQMEaKCXiiUUtmAQ2udr7XOBRKVUn6v61dKmYE8rXWWd30AgwZlb74BCQOgCt788E2joxiqc1x3ir5ZyJ9Oc7LohT8ZHUd4paSkNHpuadH2GHFEkaG1rj2t0zwgs5715wK1/wXPB7ICESyQTjrmJAA+Kmj9kxg1xNRjBE/Zz2L0Te/y0tNhV/NbJbPZXDNRjxB1BbVQ+GkycgH1/QtNAwqVUhalVJLW2qW1DrtJqDOvzIRYqOhQYXSUkHD19Fc4qouJvz/8WKucZ3vUqFG89NJLAFRWVjJq1CheffVVAMrKyhg1ahTz5nnGwNq1axejRo3i7bffBmDHjh2MGjWKBQs8Jz9s2bKFUaNG8cknnwCwfv16Ro0aRWFhIeCZwW3UqFEsXry4WVldLheZmZkopWoes9vt2O128vPzyczMrJl5LT8/n+TkZBITE3G5XDgcDpRSZGZmBnQOaWGsYB9RxAN150D0O6FsrcKSUuuxPG9zVFjpHtedxIcTiRoaZXSUkNCtZz/efOERVm2uZEamzLNtJLPZfFizU3p6OgBpaWmkpqaSlZVVc3/hwoU168XHx2Oz2bDZbDJ8eGvmb+q7QNzwHB2U1HnMDGjA7Gd9DVjrPJbnbxu9evXS3udoQM+cObPZUwO2tJRbU3Tn0zsbHSOkXDp6gAb03Mf/YnSUNs9sNtf8f+2pSqun2qwtLy9PW63WQ6b3FOGNeqZCDfYRhQvPUUVtde/XXR+g9tymDjzFwqeEhIRD3uCsWbOaETMwKlZWsHvZbrbs3GJ0lJDx1IsfExut+Ef2Y1SWy1lQoSQrK4v8/HyKiopwOg898Jd5JdqWYBcKJ54jiNrMAFprl4/1HT6WuaDmbKiwMvaisXAQPl3yqdFRQkbv/sez4JXZPHJZJVG//N3oOAJPn0VycjIzZswgLS2NlJSUQ5aBpw8jKyuL7Oxs6ZtoA4JaKLTWdv44SqgWDxT6Wd8BuOoUBTPg8lNYQlr1JEafL/3c4CShZcykLMZdPIW778/m649yG36CCKiioiLMZjNmsxmgphA4HA4cDgcul4uioiKsVis2m62mP0O0XkacHju/znUTqdQ6/dV7dlPt5bOB2vNpTvY+FnZGnjoSTPBpvhxR1FU5+G/8e7Ei/apbKN/nMjpOm2a1WklJSSE3N5fCwkKSkpJISUkhPz+fwsJCkpOTKSnxjNcVHx+P3W4nPT0du91ucHIRKEobMOyzUmo6YAcsAPqPC+mql6VqrVPrPFZDa53j77VTUlJ0UVGRv8WGi+wYSVRsFPu37zc6Ssh554WZXHLDg9zzp9PIeXmp0XHanLi4OEpLS42OIQyilCrWWqf4XGZEoQikUC8UyVck85PjJ/b/b/8h560Lj4yJA3nh/V/49zP/x59uedDoOK1efn4+BQUF2Gw2kpOTKS4uNjqSMEh9hUIGBQyya6ddS8XYCjbv3Wx0lJCUM/dTTCbIvOMh9u3abnScVs9qtZKYmEhOTg5z5841Oo4IUVIogqzLni7wL/h3/r+NjhKSzN2P5tl/TKOySrPw+euMjtPqmc1mpk+fzvTp00lKCruxNkWQSKEIsqT+SbADFrzXduemaEjm9CdY9c41TOj1Adt/esfoOEK0eVIogmxQ4iAiOkWwtXSr0VFCWv9x/+Qvb8cxeEQaG9asNDqOEG2aFAoDxPePZ0uJXJ1dr8j2DJ+QxY7dbu6/ZYLRaYRo06RQGKBjREfKN5bj3OV3PEQBTPxTFjOuH8HLH63i7dywG1leiFZDCoUBTj71ZAA++OIDg5OEvr8+/Qkn9I1h8s05fP/tf42OI0SbJIXCAJlTM2E47DC1vnkYWlp0bEceeeQJqqrg5uvT0G630ZGEaHOkUBjAmmwlcmwkG9wbjI4SFi6cfBPPzUxnclIprJpjdBwh2pwmXZmtlLoBzwiwhVrr3QFLdQRC/crsah2P7UiUKYrSX2XIhEbRbtwLz+fJ/3zO+be8zYnDxhudSIhWpSWvzI4D7sUzouvvSql/KaUuUUp1PuKUbUyHqA7sWrsLtzSlNI4ysfO4x5mVf5Cx4y+loiwkf6eErcLCQhITEwP2+na7nZycHHJyckhPT6+ZxlWEhyYVCq31I96KEwfcCOzCUzhKvYXjbikajXP2BWejKzTFv8rYOo3V/ejBPPHgbazbXsm9GecYHadVSUlJOWw61JZUWFhYcwX43LlzZbTZMNOsPgqt9S6t9UKt9V+8hWMosBA4BlitlLq4JUO2RtYzrAC8vuB1g5OEl+vvfpqb0wbz+GvFPPP3m42O02qYzWasVmtAXttutzN79h8zA5jNZlJSUuSoIow0efRYpVRnX/0TSqkbtNbPeycZygXmaK2Dfj5juPRRrN20lv69+3P0yUez9ru1RscJK/v3OOnevRv7KzVrf7XT55hTjY50mFGjRh322KRJk7j55pspKytj3Lhxhy2fMmUKU6ZMYceOHT6nGr3pppuYPHky69ev5+qrrz5s+V133cX48U3vu3G5XGRlZZGbm1s9L33Nr32Hw0FBQQHZ2dmYzWby8/OZPXs2LpeL4uJinE4niYmJZGRkkJWVhcVi8bmNwsLCQwpRYmIiWVlZZGRkNDmvCIwW66NQSs0B1iildnr7J2of/yeCZ9pSrfUkQEYYq0e/hH5Edo5kT4XME91UsZ3iee2FJ4kywdtPXg5a+nmOhNlsPqzZqXrWurS0NFJTU8nKyqq5v3Dhwpr14uPjsdls2Gw2v0UCOKRIOBwOnE4nkyZN8ru+CDFa60bfgEu9/00C5uA5A6rKe7vbu2xI7XWDfUtOTtbhYvS9o3XXW7oaHSNslSx8SOvX0OsK/2J0lFbBbDbX/H9paWnN/xcXF2ur1XrIunl5edpqtWqbzdbk7VitVl1cXNzsnCIwgCLtZ7/a1D4Kl1LqEq21XWt9o9Y6HojXWkdorR/1rpOvlPoXng5vUY9x549jZ+VOftv4m9FRwpJl9L3M+304lnMf5uVn7jc6TquTlZVFfn4+RUVFOJ2HDjfjq2msMXJycsjKypIhzcONvwri7wZ0Ac5paJ2mvm5L3cLpiOLx1x7XgL7lwVuMjhK2nFscumsnk+7bLULv2Pib0XHCWvURRWlpqbZYLDVHFcXFxTopKalmWfVjBQUF2mKx6JKSkka9fl5e3iFHEo19nggOWvCIovqMp3o7qbXWu5pRs9qc8SM8HY8rVq4wNkgYi+sxgI/fepGtriquungEBysPGB0p7BUVFWE2mzGbzYCnT6H6vw6HA5fLRVFREVarFZvNVtOfUZ/CwkLMZnPNkYTL5ZLTY8OIDOFhoGP6HUP0UdGU/F5idJSwNjT1Gv5x10Q++XY7F44ebHScsGe1WklJSSE3N5fCwkKSkpJISUkhPz+fwsJCkpOTKSnx/JuNj4/HbrfXe12Ew+EgNTWV1NRUlFIopYiLi6u381uEliafHhvqwuX02Gp9h/Vl04+bqNpXZXSUsOauquLYo7uwevM+lnz4L04be6PRkcJOXFwcpaUypExb1ZJDeIgW1imqE+4yN18s+8LoKGHNFBHBih9+4+9XdePUPbPQZZuMjhQW8vPzyczMBJBf+MIvKRQGm3zFZAAWrVhkbJBWoFN8AjOe/pwPlpZiPWMge2ViqAZZrVYSExPJyclh7ty5RscRIUoKhcEyL8+Em/BeriiOmHkwpT2u5aufdnPHlLOMThPyzGZzzRhMcsqq8EcKhcF6xvWkd2RvPlzwodFRWo3rps/hpvRTef7dn3n1iZuMjiNE2JNCEQIqP69k2b+WGR2jVXnk319x5uDOTLl7Dq/NecjoOEKENSkUISDljBR0hWbpT0uNjtJqRMW0x/Yfz5zkt979V7au+dHgREKELykUIeDyiy8H4I2P3zA4Sety4qlnsXjBi5x5vCKy+HqoqjA6khBhyZBCoZSarpRKU0plKKWaNM6wUipws6sY5JJRl0AEvPP6O0ZHaXVGjLuWBW+/zsbfv+W+G4bjrpLrVYRoqqAXCqVUNuDQWudrrXOBRKVUo0YY8z631Z3s3b5deyKiI9j4y0ajo7RO/Sbzxm/DefSVFTw07QKj0wgRdow4osjQWufXuj8PyGzoSUqpVn3uXsqFKVS1r2JXmQyTFQh/sy3m3GE9mfXcpyz49/8ZHUeIsBLUQuFnZ+8CGjMHYwpQ0KKBQsi9s+6F22DFthVGR2mVIqOimffxDyQltict4yHmvZBtdCQhwkawjyji8Ux2VFuDl896m6bmByRRiDj96NPhQ3h27rNGR2m12nfuxn/e/JAqN9xw6ww2rgqfMcGEMFKwC4XZ3wLvXNv+HndprV2N2cCmTZtqRqhUSjFr1qympzRA947difglgs9e+8zoKK3aicmj+HDecyT2UOz7/CqoPGz6dyFEHcEuFC48RxW11b1f1yStdWFjN5CQkHDIhBvhUigABgwdwJ7te3C7ZQ7oQDov7WbsX39IvHsVd195EmV7XEZHEiKkBbtQODn8qMIM4OuIwdun0egiEe7OOOMMdJmmYFmr7YoJGabe57OwYipP5q8j/bzBaCnOQvgV1EKhtbbjOaqoLR7/xSAeSPNedzEdz9lRFu/9Vnea7LmnnwvAw48+bHCStmHy7f9i2tVn8NH/NvKX6043Oo4QISvSgG3OV0ql1TpFNhWouYjOWwCSvNdZFFKriHgvzrNorXOCmjhI0lPTucp0Fb/++qvRUdqMR/79JeX7h5Dz8rf8snYQ733+k9GRhAg5Qb+OQmtdfVRg9e74S+pcV5GGj+sqvOum88cRhTkogYMoOiqafqn9qEqUq4eDRZlMPP2GnX492/P+op95/h/XGh1JiJBjyBAeWuscrXWh1jrXe3V23WWpPp6Tq7VO1VrHeddxBS1wEN10/01sO2UbG3fJVdrBYoqI5Mdf1jP+9G702f0yrHvL6EhChBQZFDDEDI4ZDI/D/z0qVw8HU8cu8by/eDXWkcO5YUo6rz59t9GRhAgZUihCjPUUK1TApx9+anSUtieqI/tPm8/nK6O47s7H+G/+Y0YnEiIkSKEIMTHRMfRO6s3W1VvlegoDdOrahy++Keb4PjGMv+puFr37tNGRhDCcFIoQNHL0SKpcVby36D2jo7RJvQcMpuDzb+kVF8U5l/yZZ2dPMzqSEIaSQhGCLjrrIgAeelim8DRKzwEn89a7H9EuSjF91lOsX57f8JOEaKWkUISgS8+9FFOsifWb1hsdpU07ZZiVkpV2Hr+hB31/n8Km76RYiLZJCkUIioiI4Jxp51B2VhlVbrmmwki9LEO48dEVvFlkxpKSzlMPNjh1ihCtjhSKEHXdZdexT+/jo+KPjI4iYntivX0hPeOjuOfBXN7OnW50IiGCSgpFiBrSeQjY4IFZDxgdRQDd+hzPih9+Z+jxnZh00yP8ecq5RkcSImikUISogf0HEt0tmt9//93oKMLLfFQ/Pv3qN/r1iOXplwv4yw2jjY4kRFBIoQhhp4w+hd2O3WzYvsHoKMKrY1xPVvy8nlGndmfl8kVULbsTtFzvIlo3KRQh7E9pf4IqmDZrmtFRRC2dzF1ZuGwT856+kW8/eoIrz7Owu3S70bGECBgpFCFs6sSpoODj/I+NjiLqMEVE0u6Mf/K/fWN56/O1nH9mIs7NJUbHEiIgpFCEsJjoGBLPTGR/5X7KDpQZHUfUpRR35nzEnNmZFP+2hzOGDuSbwvlGpxKixUmhCHEPP/kw+lbN52s+NzqK8GPK3XMofOsZ1m2v5MxzJ5M3d5bRkYRoUVIoQtyEUyYQY49h5t9mGh1F1OOsCbfy3vx/07VTBLP/8QDu32wNP0mIMCGFIsRFR0TTcVVH7Pl2Kg5UGB1H1CP1oimsWuXgvYdHsvHTG7n7qiHs3eU0OpYQR0wKRRi4eurV6HLNnLfnGB1FNKBL96Ppm76Qd9eP5qk3vmPU0KPZ7FhhdCwhjogUijCQdW0WRMJLb75kdBTRGKYIbsv5L8/87XpWrttHUnIy+S/8zehUQjSbFIow0DO+J516dmLFByuoqpJBAsPFjfc+z9L/vo1SivQb/kpG+jDcVQeNjiVEk0mhCBPDRgyDKnjytSeNjiKaYPAZF7Nk2QoG9uvE2598y/b3xkLFTqNjCdEkUijCxOtzXify1Ei+cH5hdBTRREcnDuYnh4vi9x+k897FnD88AdsjdxodS4hGk0IRJrqbu3PtA9fy2cbP2OLcYnQc0UTKZKLf6P9jy0lvsWJNFTdNf4K/XDucA/v3Gh1NiAZJoQgjI2NHUv5oOTfcc4PRUUQzDTh1PL+uWsfUiSeQ/dJSBlq68t5rzxgdS4h6SaEII5eNvoyIThEsWrjI6CjiCHTpmoDtnZXM/+ddrNl6gEuuvp0XH5oEMpuhCFFSKMJIREQE5119HvvW7uPjJTJQYLhLv+lRvi9azLjhPTmuMo+dbw9n5bcfGB1LiMNIoQgzs++cDQpuv+t2o6OIFjAo6WwWfL2JM6e8yv25P3DKGeO5P2MUlRUyCKQIHUprHfyNKjUdcADxAFrr3HrWNQMZ3rtDgYL61k9JSdFFRUUtFzYExZhjOLDnAGVlZcTGxBodR7SQzY7vSZ94Dl//sJOTB7Rj5l//yiVTZhgdS7QRSqlirXWKr2VBP6JQSmUDDq11vneHn6iUSqvnKTO01jneWzqQpZTKqGf9Vi/znkwwwePvPm50FNGCellO5qvvd/Du8zPY6jzApdfey5ihvWWeC2E4I5qeMrTW+bXuzwMyfa3oPZqw1HnYBmQFJlp4ePwvj3PsQ8eS78rH7ZZpOFubi67/B8u//4Xhg3uy9IdNON8ayo5vH6eq8oDR0UQbFdRCoZRK8vGwC7DW8zSrUqp2sXBxePFoUyIjIpl60lRWPLSCB2wPGB1HBECvo4/lfz9sZt3Pi0k8bhDpU+4isU9Hvv7gX0ZHE21QsI8o4oG64y77HYdZa+3SWsdprR21Hk4FCgMRLpxkjsxEORWPz5bmp9Ys3nI2WBcz+crr2VN2kDPH38zkc/rytcykJ4Io2IXC7G+Bt5mpXt51rNTT9LRp0yaUUjW3WbNmNT1lGOjcoTNjp4xl7/q9vPrJq0bHEQGkTCZuvO951m3YwswbR/Lelxs4M3Uyl59/Aru3rzE6nmgDgnrWk1LKCuRpreNqPWYBSoA4rbWrgefnATattd8jirZw1lO1dVvX0a9fP3oc34Mt38mwHm1F8dcfc92UK9m6vZTfn+7Eus7X0G/kfXSM62l0NBHGQumsJyeHH1WYwdPMVN8TvafU1lsk2pqjexxN3+P6svX7rbyY/6LRcUSQJI8Yy3e/O/nt+//Rvu9oJt74LEf3TeCJv0ykbPcOo+OJViiohUJrbcfTGV1bPA30OXhPn7VXFwnvkYkAXn3hVYiEB56TTu22pvPRw4kY/R5PPv4oAxI6cmf2eyT260H6eUPYsWWd0fFEK2LE6bHz61w3kYrnlFfA0xRVe7m3KMQDRUops7epytfZU23S2UPPJvPpTNaNWMfiNYuNjiMMcMEVd1H8224Wv/s03cyx5H/2Hccm9ufXD+6CA7uMjidaASOvzLbjPc219pXW3mWpWutUb+d1qY+XyPdefHeYttRHUW1/5X76ZvalsqiS7cXbiY6KNjqSMNDcx6azZOHr5F6xkecWtmPJVguPPPMKvRPl95XwL5T6KADwXmVdqLXOrTsch3dZqvf/XVpr5ePms0i0VbFRsYztPpbdP+/m+v+73ug4wmBT78rhhY82EDGuiJ929uLtgp+xnJDMdRcex6IF0pclms6QI4pAaotHFAAHKg8Qd2wc5a5ytq7dSrcu3YyOJELEr/YCnnn4Ll587wf2H4BjEmJ5dc4DDLvgTjBFGB1PhIiQO6IQLS86KppHHnsE9y43p513mtFxRAg5PimVZ+d/z88/2Dn/zOPYW1bBgWXTWZV7NPdmjOb3H5cYHVGEOCkUrcjNl95M92O7s3rpamY+M9PoOCLE9D/uVD7+8lfWb93HmdfP59OfO5D9/CIGDTmdy8YczTsvzuKgjCclfJBC0cosWbiEiM4RPP7a4+w9IPMxi8NFRrdD9Uvnlqd/4/MP/8PN6UP4bOkGLrn+ATq2j+Hxe86HfWuNjilCiBSKVsbS18L8z+az97y93JZ/m9FxRIg7e+zVPPnGcjZtdXL7NefRr2csvxR/Cu/156bxCcy46QJ279xgdExhMOnMbqUmPTKJvBl5XHPHNbz0yEtGxxFhpGpXCbt/fIGTL8phw84q2kXBBWf05tjBw8h64DnMXWWokNaovs5sKRStlHO3k579e1JZVsmHn3/IuNPHGR1JhBl3VRVffWgj77U5vPLBj+wq08RGw6LnzmfQ2dcS0ftc2nUwGx1TtBA566kNiu8cz+IvFqMiFRMnTmTl6pVGRxJhxhQRwdkTbuaZed+zaWspf//LFK44/zgGxy7hhYcm08Ucx8SzejPv2T+za8d6o+OKAJIjilbumfnPcPvk24nqFMWubbuIbSdzbIsjVHWApR//i3tn5fDjqs1s26UxKejWJYqH7ryUqXf8AzoOMDqlaCI5omjDbpt0G6edfxqVeyoZeeNI3FqmThVHKCKaYRf+mYVFG9m0vZwv3nuGs5L6cKCyivfefhPet5B5fhwTRx9D4fwcKsp2G51YHCEpFG3A0o+XcsFdF7AsYRnXvXadzLMtWkxEVDRnTbiVRUXrKd1bxbyPluMe8hirdyg++KKE1MlZxMV1oXe3dmRMOh3Hio+hlbVitAVSKNqIBY8s4KZTbuLlG1/GMswixUIERIdeQzCdeCefFTlx7djMgpf+j8vPH8h2VwVz85ZQ+Ow4Sl/pwcSzenPfTWP5zf4pWv4thjzpo2hD3G43icMTWbNsDQPPH4j9PTvtotsZHUu0EWt+/hpzxQpK7B9g/fOnuPZ59j3dOiliY9txbfpIpt01g7h+I2QMKgPI6bGihtvtZsSVI1jy5hI69u3It4u+ZaBloNGxRBvjrqri9xUFLProFd54u5Cvvt9GlRs2Pwdf/N6Ov78fyXlnn8w51rGcNuYK4ntajI7c6kmhEIe55NZLeOe5d4joGMHibxYz4qQRRkcSbZh2u9nw2xL6xqzmjdde4u4nFrGl9CBu7+4pOhKuPNfCC4/9mS0HLXTqM4yO5u7Ghm5lpFAIn26deSv/fPafdMzoyDMXP8M1p11jdCQhauxxbmLZf9/A9vxLFP/g4IReVXxwZwWXPgkffwfHJLQjeVAfYjv14JTk4Vx+7TQ6d+1jdOywJYVC+PXrtl+5+q2rWfbAMnp064F9kZ2EbglGxxLicFrD/o0UvvsCL74yH5dzO/bfdrLV5ekMv2AIfDBrAH+ZH8m60mgunTCWk06zYhk8isjoGGOzhwEpFKJe+yv2M2zCMH747AcizBHcMesOHr71YSIipENRhL4VSwpY9sUCju1exshjdpM05T2+W32g5ixcBfSIi+TvN6Zw3ZXjeX/pHrr0GMgZ1jSiYtobmj2USKEQjWJ7x8a0W6dRvqmc6E7RzH9vPheNvsjoWEI02b5d2/ml+FOWL13InFc+Rh8sY9q4KCYMdmKe6lknMgISe0aza79iZMoA7r39Mk48dSTby+Pp0W8Qpjb2Q0kKhWi08gPlDB87nO8WfQfd4Ia5N3DPsHs4rsdxRkcT4ohVljn54pM3+d7+FVs3Olj+k4NFy7fjdsO82+H4XjA4C2Kj4ISjY+nVvQs79yrGp57G9VOuotvRp0D7vkS2wtPKpVCIJvtu1XfM/mg2b29/m8qnKumd2JsX5rzAecPOMzqaEC3u4IFy9N41ONcv54lnX2DL5o1s276T70ucbNxZBcDSB8G5F8Y9Aj3NEZzQrxNdOnfGVWbi4rFnMOVPlxNtTkS170NM+04Gv6Omk0Ihmm2FYwXj08ezYcUGcEO3wd248NwL+eff/ykDDIo24UD5fjasKqJ3p304flnGQ0/No6JsN5u3u/h57d6aCwc3PQsLlkPmC3BUZxP9esYSE9OOPeWKtPOTueeWy3FWdGLzrgiOP3kEHbp0M/idHUoKhThiPzp+5NYHb2XxG4vhAJg6m7jjlTu4dui1DOo1yOh4Qhhmd+l2tq39HstRmuXLvuQx29u0i6hg4xYn35e42FJaRVQElL8EM9+Ch971PM/cQdEuOoLKg3DxqP7MfeBivl5ZzuLl2zjj9OH06HMMXY4aQM++JwSlv0QKhWgxW3Zs4Zb7b8G+zs664etwv+Mmak0UZ158Jn+e8mfGnzkek0mGEBOitj07N9Apcjc/Lf+a1+e/T/vIcjZv2cpX9nWs3bKXE/tG8vVMzXVzDvD6N1BReejzR5wQzVePD+LR98pYsLSUoYP70r17d0r3KbrEdWfc+VZOST4dYhMgsnlnckmhEAGxYfcG7px9JwteXED5tnIAVJSi24Bu3HbfbdyUdhPd2ofW4bUQIUtrdm1fy/Jli9HlO9m8oYTCxctYv2kbI0+J5/7LE7j18WXMX7yDPfvdlNcqJmefAIv/D0h+Go6/rVmbl0IhAm7F7yt49KVHeWvOW5TvLocRwEho93I7usV1Y+ylY7nYejFjkscQHRVtdFwhwpp2u9m3axvf279i8/oSBvSMJek4M3QdCl2aN3abFAoRVGXlZXzt+JrFJYt59s5n2bVmFxz8Y3n7hPZcePeFnDP8HNRmReqwVAb0lhnRhDCSFAphqAOVB/j4fx/z6luvsvCdhZTrcvQETfm+cnjZs05EXAQdu3YkojKCoWOGck3mNRzX4zj6dOpDj/gexr4BIdqAkCsUSqnpgAOIB9Ba57bU+lIowkOVu4plq5Yx96W5rFq/irVr1rLl1y1UbK/wrDANWAO8CypW0a5bO9rFtkNVKE46+ySuuPkKYitj2bd5H0NPHsogyyCZW0OIIxBShUIplQ0s01rn+7p/pOtLoQhvVVVVrFy3kj0Re/iy6EvmzZ1HmbuMndt2Urq6lIO7DnoG77kf+Br4b60nRwFu6HxMZ0bNGsXe7/aydvFaLIMtdI3rijqgiFJRnH3O2Zw+7HQoh2ii6Z/Qn8iISEPerxChItQKRanWOq7W/SQgW2ud2hLrS6Fo/Ur3lOKsdFL8UzEL3llAZVQlW7ZuwfG9g53rd9K+X3sSbkjAke9g7+K9h/SPANAFuAN4D/geqAKiAQ0KRZfBXRh822C2fbKNbUXbiO8bT0xsDGXbyjCZTJxqPZURE0ZQsrSE9T+tZ8DxA+gQ2wHXVhft27cnZVgKx59wPDs372SPcw99+/SlU/tOKLfC3NlMfJd4KUwi5NRXKIL6r9W7k6/LBVhbYn3RNsR1iiOOOBLPSmTSWZP8r3ij5z+uvS7WbF7Dih9XULKmBHOCmT6D+rCk+xKWfLqEaHM0e/fsZcPKDezfvZ+4nnFEmaLY5dzFXude9rr2UlVRhd6vwQ2rt6/m7XZvw0vABg4vRMcCVwK5wFY8hai2QRAxKQK3zY3erlHRChWhcJe5MUWZiBsZR8KEBNY8t4aKTRW069YOU4SJsk1lRHeKJnFcIpYxFoqfK2bfln10TuhMREQE23/bTsfuHRkyYQiJwxL5Ys4X7HPuI753PBGmCNb/tJ6uCV05bcJp9D2hL4teWUT53nK69eqGQrH2x7Uc1e8oho8dzlG9j+LLd76ksrySbj27UXWwinW/rKNXv14Msw6jU5dOLC1YinZr4rrGcfDAQTY4NtC7f2+SRyTTrl07ln+zHJPJRBdzFyr2V7B5w2Z6H92bk5JOIsIUwU/LfyI6OpoOHTtwoOIA2zdvJ6FvAscPOh7t1qxauYqYmBjad2jP/rL9uJwuEnon0N/Sn4OVB1m7ei2xsbHExMRQvr+cfXv20b1Hd/r06UPlgUq2bNpC+w7tiYqMoqKigvKycszxZo466igOVByg1FlK+9j2REREUHmgksqDlXTq1IkunbtQWVnJ3r17iYmOISIigqqDVWitaR/bnpiYGNxuN5UHKomMjPRcN6RBa01UVBQRtaZxVUod0b/1UBHsnzXxgLPOY3XvH8n6QhzG3NHMkGOHMOTYIYc8PnnwZLi5nif6mMfJ7Xaz/8B+Kqlk05RNrCpZhSnGxJ59e/j5p5/ZX76ffgP70ef4Pnzb41t+tv9Ml6O6UF5Rzu8rfgcTHH/m8Rwz4hi+XPcla79bS8ejOnKw8iBbVm4hIiaC/kf3p298X3Z22UmVswpTpAl3lZuqyioqyirY4dqB2+lm6+qtVJZWssu5C12lqdpVxe5tu3H2cvJlzJfs/mS3p4gpwA1oWP/9elboFbAe+Lfvt/3pjk/hOOAx38tfXvUy9AKe9fN3uwLoiKdQ+jLFk6X6RIbDZAB7gDf8LL8N2AS85Wf5XcCvwAd+lt8LFAGf+Vk+C1jkvdUVA8wAPgWW4HkftXXybn8BsILDfyR0A9NtJtxvuWFlreXac1N9FFGZURx84yDuVW7P58Yf65mOMRFzTQwV/6nAvcb9x+tWgjIrbAU2piZP9fPGmi/YhcLsb4FSyqy1dh3h+mzatOmQKj5z5kxmzZrVxJhC+GYymejQrgMA5n5mTux34h8Lzz103YsHXlz/i53TwMYmN7D8pgaWzzj8ocqDlbi1Gzdudt21i/KKcpRJcaDyANu2bSMyJpJu3bsRGR3JmivWUFZWRkxsDBUVFaxbu44u5i70P7Y/EVER/Hzuz1QeqCS6fTTl+8tZ61hL/FHxnHjqiZhMJoqSi9BoYjvEsnfPXtaVrKNn354MOWMIbrebpQOXokyKdh3bsXfXXtatWkefxD4MHT2U/fv3s6TfEqJjo4mOiWZX6S42r95Mn2P7MHzscPbu2ss3Pb6hfZf2REVH4drhYsuaLfQ7sR8jLhrB9o3bWdJ1CZ3jOxMRFUHp1lK2b9hO/8H9OXvc2Ww8biNLOy8lrmccyqRwbnGyc+NOBpwygFGpo1jdYzXLOi6jW59uaDQ7N+7EtdXFsUOPZcQ5I/i90++s6LSCrn27orVm5/qd7CndwzFDj+GMkWfwa9Sv/NDlh5rlO9buYP/e/Rw77FhOO/M0fjr4E7+YfyGuTxxaa7av2U5leSXHjTiOpOFJ/LD7B37/5nfienta3beVbMPtdnP8WcdzytBTWLFlBauLV2PuZQZg66qtmBPMnNj9xMM/9BYQ1D4KpZQVyKvT52ABSoC4ujv+pq4P0kchhBDNUV8fRbAH5XFy+FGCGcDXTr8Z6wshhGhhQS0UWms7ns7o2uKBwpZYXwghRMszYpjP+UqptFr3UwFb9R2llKXO8nrXF0IIEVhBLxRa60zAopSyKqUygJI6F8+lAZlNWF8IIUQAGTJxgNY6R2tdqLXOrTsch3dZamPXb0lydlRoks8l9MhnEpoC9bnIoIC1KKVobX+P1kA+l9Ajn0loOpLPJZTOehJCCBFmpFAIIYSoV6trelJKbQfWNvPpCXgGBxChRT6X0COfSWg6ks+ln9a6u68Fra5QCCGEaFnS9CSEEKJeUiiEEELUSwqFEEKIesk0WzR9Dm8ReN6r8JOBPO9D6XhmNnQYl6ptUUqZ8cwO0VVrneVjuXxvDFDf5xKo702bLxS+5uRWSqXJMCEhYRKeL4QdmCpFIni8Q/ybgUQ/y+V7Y4CGPhevFv/etPlCAWTUqcrzgGxA/sEbrPY8JCK4tNaFAEqpofieQEy+NwZoxOcSkO9Nm+6jkDm5hWg6+d60PW39iELm5A5h3vZWJ9IGHmrkexPCAvG9aeuFwuxvgb85uUXQFAGu6vZVpVSeUsopbeAhwexvgXxvDBeQ702bbnrCc7gcX+exuveFAbTW9jqdcMuAGUblEYdwId+bkBSo701bLxQyJ3eI8p7dUZsD8NU2LoJPvjchKlDfmzZdKGRO7tCklLIABd7zxWuT02NDgHxvQlMgvzdtulB4yZzcIcZ76JxV59fpZDynX4rQIN+bEBPI742MHkvNFaZ2wAJydk0o8P46qt4RdcUzV7p8LkHiPQXWyh/z19uAQu/RRPU68r0JsoY+l0B9b6RQCCGEqJc0PQkhhKiXFAohhBD1kkIhhBCiXlIohBBC1EsKhRBCiHpJoRBCCFEvKRRCCCHqJYVCiCDxXgwlRNiRQiFE8OTVGfZCiLAgV2YLESRKKQ3EyQirItzIEYUQQeAd/tkhRUKEIzmiECKAvAUiFc9AbS48Q3Evk5n6RDiRQiFEECilioHZUiBEOJJCIUQQSP+ECGfSRyFEgEn/hAh3UiiECLxUZJpQEcakUAgReFagGEApZfYeYQgRNqRQCBF4Zv44osjQWsvRhQgr0pktRIB5r8a2AA488xu7jE0kRNNIoRBCCFEvaXoSQghRLykUQggh6iWFQgghRL2kUAghhKiXFAohhBD1kkIhhBCiXlIohBBC1EsKhRBCiHpJoRBCCFGv/wfKjmu8MxjVHgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1)\n",
"\n",
"y0 = -1.0\n",
"__=ax.plot(t_domain, analytical_dy_dy0(y0, t_domain), color='green')\n",
"__=ax.plot(t_domain, ode_result_jax_grad_vmap(y0, t_domain_vmap), ':', color='k')\n",
"__=ax.plot(t_domain, analytical_ode_result_jax_grad_vmap(y0, t_domain), '--', color='k')\n",
"\n",
"y0 = -5.0\n",
"__=ax.plot(t_domain, analytical_dy_dy0(y0, t_domain), color='orange')\n",
"__=ax.plot(t_domain, ode_result_jax_grad_vmap(y0, t_domain_vmap), ':', color='k')\n",
"__=ax.plot(t_domain, analytical_ode_result_jax_grad_vmap(y0, t_domain), '--', color='k')\n",
"\n",
"dashed_line=mlines.Line2D([],[],ls='--',c='k',label=r'${\\rm jax\\ 2}$')\n",
"leg=ax.legend(handles=[green_line, orange_line, solid_line, jax_line, dashed_line])\n",
"\n",
"xlabel = ax.set_xlabel(r'$t$')\n",
"ylabel = ax.set_ylabel(r'$y$')\n",
"title = ax.set_title(r'${\\rm JAX\\ solution\\ for}\\ \\partial y(t)/ \\partial y_0$')\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "specified-amendment",
"metadata": {},
"outputs": [],
"source": []
}
],
"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.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment