comparison2scipy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numbalsoda import lsoda_sig, lsoda\n",
    "from numba import njit, cfunc\n",
    "from scipy.integrate import solve_ivp \n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "@cfunc(lsoda_sig)\n",
    "def f(t, u, du, p):\n",
    "    du[0] = u[0]-u[0]*u[1]\n",
    "    du[1] = u[0]*u[1]-u[1]\n",
    "\n",
    "@njit\n",
    "def f_scipy(t, u):\n",
    "    return np.array([u[0]-u[0]*u[1],u[0]*u[1]-u[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "funcptr = f.address\n",
    "u0 = np.array([5.,0.8])\n",
    "data = np.array([1.0])\n",
    "t_eval = np.linspace(0.0,50.0,1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "usol, success = lsoda(funcptr, u0, t_eval, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAFHCAYAAAABXlMnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACGB0lEQVR4nO29eZgkV3Xm/bsRudRe1dVdvXertbT2FQQCxA5GBlmAbDBgwGCw8YYXbMOMbGNj42Hs4RsG2xhjZrAxxoBBgFkMyAIkARYSau1bt9bel9q3zMo17vfHjciMzIqIjMilKiI63uepJ7uzMrPiZtx7zznvec+5QkpJggQJEiRIEAZo630BCRIkSJAggYXEKCVIkCBBgtAgMUoJEiRIkCA0SIxSggQJEiQIDRKjlCBBggQJQoPEKCVIkCBBgtAg1es/sGnTJrlnz55e/5kECRIkSBAR3H333dNSygmn3/XcKO3Zs4d9+/b1+s8kSJAgQYKIQAhxyO13CX2XIEGCBAlCg8QoJUiQIEGC0CAxSgkSJEiQIDRIjFKCBAkSJAgNEqOUIEGCBAlCg8QoJUiQIEGC0KClURJCvF0IIR1+fm0tLjBBggQJEpw+CFKn9FJgxfb/p7p8LQkSJEiQ4DRHEKN0l5RyuWdXkiBBggQJTntEIqd0YmGFlVK18cnDd0JuZn0uKEGCBAkS9ARBjNKTQoiKEOKAEOJXe3ZFTShXDZ77P7/Pr3zG1qqoUoR/fAX83bPX6jLWDeWqwa/+yz6+89CJ+pOVEnz6Z+ChL6/fha0RlosVPv+Tw1SqRv3JShEe+TpUy+t3YWuEqiF5+PhC45Mr83Dvv4KU63JNCRL0En6M0gng/cBbgeuAO4FPCCHe08sLs/DAUbUgf/TEdP3Jqf3qMT8d+4V5y/5Jbnr4FH/1nQP1J5/8Phz8IXxzTW7BuuKzdxzihq88yKdvP1h/8v4vwBffCj/8yLpd11rhvTfez7V/8yOemrIx51/5Ffjab8CJ+9btutYK77vxfn7nC/c2PnnTH8G3/9v6XNAaompIPvKfBzg6l68/KSXs+ydYnly/C+sxWholKeVNUsq/kFL+p5Ty21LKXwS+CPyxEMLx/UKIdwkh9gkh9k1NTXV0gRm9/idq3vLMk/UXrMx19Plhx9PTOQAKZRt9OW/2MqwU1+GK1hZPTqrN2PoeAJh+TD3OPL4OV7S2+Pp9xwF45MRi/cmTD6nHI3etwxWtHQrlKl/cd5Sv3Xe8Pv+lhB9/DO78BORn1/cCe4wHjs7zN99/gnd+2sYSnXoYvvm78E+vXLfr6jXazSndCIwDe5x+KaX8pJTySinllRMTjt3JfeOSnaP8+WsuAmAub9I1K7bJOHewo88POw7PKi9puVipP2l5SdUSGIbDu+IDa/ynFm0GePbpxscYI206ZU9O2oyyZuqT5l0bLccC9gjBmgcs2WjsmN//+47MAzCTs819iyWaeWLtL2iN0KnQYU24s01DWQCml82bY4+OctMO74gPTiwUAFgqVKga5te9fEo9SkNRmDHG1JK655NLhfqT1pjtG1QMsVQos2JGCHP5knpSSsiZTknMx398vn7PZ5bN8edszMvC4TW+orXFsbmV1U/aDXFMc6rtGqWfA6aBNXHVNg5mANvEXJmv/zLm9N1MrlT79+KKOQnthnjx+Bpf0dpi0jRKlnECoGAm/mNO30wv1+/9gnXvi4tQMTfrxXgbpZMLdaM0a62DvE1xu3B0ja9obXHKnPMzuVI9dWFniWKaV2pZpySE+DLwE+ABQAfeYP78tpRyTbijoT51mTUKa2UO0gNQzkNhfi0uYd0wawvdF1bKbBjMQHEJhGZGSvGVxRfK1do9b6AvLaeksgLlFUj3r/3FrQGa7z3Q6JBZEXNMMbVcH3/tu7A7IjF3Sk4tKqMspTJMW0b6mhzyWRjdsT4X10P4iZQOAO8Avgx8CbgQ+EUp5d/28sLsGMoqo5QvmRtTYQHGdqt/xzxSmsuV2bNxAIB5a2MqLcHITvXvwoLLO6OPxYIa74aBNLliBWkpLQsLkB1R/47xxmQxA6P9aeYt+q5kqvD6x2N970FFR5mU2qJq+WT7/Y752p+2GeWlgrX3zddfENO570d994dSyvOklANSyn4p5TOllP+yFhdnYSCjjFLO8pZLOcgOQ99orCdmqWKwXKxwxsZBQOUY1C9yMGoZpfn1ubg1QL6o8ilbRvowJCq/UimqCGn8TPWilXguTKhTVmdNDNroO9Moje5Q9z7GJRFz+RITQ1kyKY2c5ZAWTUM8dkas5z7A4kqFLSMqn56zs0T94+a/4zn3I9HRYTCrA5CzujqU84q+6xuLtVGyIoWtI30A9a4WxeW6UbKH8zGDRdltNse/XKgo6hLUpgSx9RahTtltH+snb937kjn+0V1gVNRaiCnmciU2DKYZzOg1B4XyiqKuhzbHeu1LKVlcKbN9TFHTdaM0X3fIYjr3I2GU+tM6QtgjpTxkBqF/Q6wnphWyW96SpcSitAyDm0DPxJrCse73lmE1/uWibRMe3aUeY+otgnJKdE2wcTBjM0qmNHzEzCXE2CmZy5fZMJBhIJOqR0rlAqT6VbQQ47VfKBuUqkbNKC1Ze1/ZxpLEdO5HwigJIRjMpMgVmyKl2Bsl5SlP2CMlw1BGKTOk6MsYUxjWRrzFipSKFeWQQD3BG1NvEZRTMpRNMZBJ1R2SGn0X/5zicrHCcF+KwawtUqqYwpaYr32LJdk+quZ+zSGvFBVDlB6IrUMSCaMEisKr3Zhy/rSYmLVIyYwU8qVqPVLIDKrJGfNNCWDziMP4a5FCvI3ScF+K/rROqWKoOrVSs1GaX7fr6zXyxQoDGWWU85ZRLhdsa39+Xa+vl7BTt2AzSrXxj8fWIYuMUepP6xQqFoVxutB3amJaOZVaoh/MiTkW64VpLcRxs06tUK6qnAKosacHIR/v+z/Sl2Ygo3KqK2WbUR7eph5j7JTkSlWGslakZHNIU33q/hcXY1tAahmlbaPKKC3XIsUCpLIwsCG2DllkjFI2pdf7X5VzptBhVC3KmCqQFs1IaeNgBk2Y9J1VOKlnzEhpft2ur9ewhC11o2TUjVJ6QCkwi4tub488Fs1Iqc80SvlSRXWIB5Xoh9g6JVJKcsUKAxmdgUyqXqdWKUC6TzmkEFujbBXKW/nkYsWk7qvFek4tiZTWF9m0RrFiKM/IqJib0pD6d7XU+gMiCIu+U95yStFXVTNSSvXVjXJMYUVKGwdtC9OKFNL96v6X4nvupKLv0gyklVEqlAy1KQu9LguO6f0vVQ0qhmQwm6Ivrau1D2ax9EDdKMWUKbEipQ0DGTK6phwyyyFNZaFvJLYOWWSMUl9Kp9h8YzJD6t/FeG5MFn031KcW5krZ5imnMqcFfZdNabWSgAb6Lj2g7n9M7z1Y9F2KfitSKpsOmOWQQGwjZUvUNJjRyegaJcsoVQpq/Nlh9f+YbsxWpDTSnyab0tTct/a+dD9khmM79yNjlLJpTeWUapuyzShZtRsxw1KhwmBGR9cE2ZRGqSLrkZKeVQszxpFCrlSpecqA8pYbIqV4j98SOvSl1TJV9G1ROSR6SuXUYhopWVHyQDZFJqU1RUr9sXdIF1YsliRF1ooUaw55n8kSxHPfi45RSmkqUrKoOj2jbgzEdmIuFcoM96UByKQ0SlWjLnRIZZXYo1qqG+qYIVesMpjVa0apwVtM9cU6UpJS1u5/RlfjL1dNp0RXdGac6UurHGAomzIdMlvxrLUpQ2zHv1ws05/WSekafWmNop0lsM/9GObTo2OULPVdLVLIqE0ZYjsxLU8Z1GGH5UqzUTIpjJiOf7lYYTCTos/sf1YoG3W1lZ6OtbeYK1UxJAz3pWr930rW/U+ZRikzVC+mjRksYcNARq87ZGAKHQbqcz+mTslyscqg2fOzz9r7aspb0yhLm5MWI0THKFmRUgN9F+9NealQqXVIT6eEWph2+i7mRjlv0ncpXSOlCRUpGaZR0tKxjpSsfOJwX7pmlMrVZqM0GNvxW82Xh7KpxpxSeUVtyqfB3B8yc6n1vc+KlPpjbZQjY5RqCpwafZeOPX2XK1VqHdIzumZuSjahQ8zHv8pbLBtQNaXBtUgpnmO3cirDfSnSugCoz/8afTcc20ipllPKqEjRkKgzhSpmm6GY03e5YrXWiHpVpJTK2sYfP6YgMkZJeQvVpkgh3hMzX6wyaE7MtK41Jjsbxh/PjSlfbPQWCxUrUhKg6cpbrBTqhipGsHIqAxmdrEXf1SIlVbdFZjCWmxLU1XdDWTt9Wa1HSmkzUoqpQ6ZYAjX3+9JaU41evIUeETJKVk7JllOI8Y0BFSkNmBMzkzIpjOrpoz7MmW1mQBnlStXMKelK/BFnb9HalAcyKdK6Sd/Vckqqw0ecc0oWfTeQVZJwgFKpqPIoqX7QNHP8MV37trlfaxzQrL6DWI4/MkYpk9IoVyXSXqcU4xsDylu2IqU6fdekvoPYbkzLxTp9mU4JpT4zKiqfBLF2SqxNeTCr1yMFK6eo2yKlGI4d6m11BjP1SKlcsJUDgJlTjJ9DAkroMmhjCcrVJkl4klNaf6Q1xatXyzZJeCqrNqiYGiWrzQo4REp6NtY5JSmlMsrmwkzrpgKrWlY1OhBrpyRXcoiUmoUOMc4p5UsVNKGoq5pRrhklM1KMcU4xbypPQc39clWqZqxQV99BLFmC1HpfgF+kzYlZLRfURVveYjaeCqyqISlWjAb6qjFSygDx3ZSLFdVmZsAeKVYMlVOqRUrx9RatBqQNkZJF39kjpZJZqyLEel1qT2CVAwghajm1ctE0SilbpBRTo6wipfraV/feVqcklLMWx7kfnUjJ9BYrZVtOBWLLK9vpG7BFSs3FoxDL8VvqqyHbwqwYUkVKWnOkFD9v0R4p1XIqFn1nzykhY3n6bN6mvLTGXy2YBihtG38MN2V7M1qATEo0OaRJTikUyJiyWMMKYWveYlyNUn1TAou+ko0dLdL96mjoGHqLdvUZQFo3F6ZRqQsd4pxTsheP6vZIqdSovoNYjr9Z5ANQKduUtxDb4ml7M1qAlGayJA3qO4sliN/4I2OUUpa3VLZ1dIDY0ne5YmOkVGu1YtE3QqifmHqLyw6RUsnqEr8qUorf+HOlKpmURlrX0DRBSjONckObofgWj+dsOZVsSq2BGktid0piOPfzxWaHzMwpVYqAUOtfTykaMzFK6weLvjPsQgc4bSIlS33YUDwJ9bxCzFA3yvbxN0dKMc4plVQzXgv1vEIzfUc8779N5FKLlCpWNw+bUxLLsTfO/Vo3l4rZ98/KH8Z0/BEySiZ9V8upxFsWW9uUbfRVqVbRbjdK8ZyYVk7Frr5bJQmPc07JVtEPNqekuXgWYjv/B20OGUDV6mZiGaWYRko5mxweVE6tXqOXqb8wpuOPkFFyiZRienxBLVKqJXt1qobEsEuCwYyU4pdTao6UajkluyTcUiHFcGHaK/rBjJTKFaU+XEXfxe/+50sOQodm+i47rKIHo7oel9gz5GyFw6DuvSHBsM99SCKl9YZllGSlKdkZ20ihKVJKqUhRlouN3lI2nod9LdciRbvQo0kSLkRsF2au1BgpZVNaPZ/aHCnFMlKsONB3VqRk0bfxbMqab4qUag55pVQfO8T2oL8IGSWLvispxZnlMcSUvqslO5u8RaPsFCnFcfxNOSWrTsveZghiuzDzxeZISTSWA0Csex82ttmxIiWrxZiNvoPY3f9cUzlIbe+rVurUJcRWfRgZo5SxR0oNkcKQ8p6tCComaI6ULG9Rns45pYqVU4o/hdEcKWVSmnJIoLF4FmK3KRuGJF+uNohcwKSvwCZ0iKf6sFaj2JRTk5VSI32X5JTWFyndfmPsm3I8efVV6rsGo9wUKcVwYi4XK6R1UZMDp1PCJVKKZ07NSX1X65BvLxyH2I1/pVxFSptDVqOvbGdpQWwjpeUaS6LGn9JsRtlO32Xj2dEiMkbJCmFVnYZ9UxpQjzG7OflShZQmal5SQ07NHinFtP9Z3kbfgFqYq3JKEFujlCtWa9QtmN5ycz41lVFRU8wonGZJdC1SqjQLHeIpic+vyqea+eRqpWnvS4zSusLalHFSn0Hsbo6SBNc95VoI7xQplHNgGGt9iT3FcrFaK5wFW53SabIwnSIlw37Ao4UYjr+W6G8SOhjNdUpxXfsmS9KfbqLuq2V1jpgFe+/DGCF6RmnVphxPCsM6CtyCNTEbOhpAffwx63+2WhJtHV3RPP74CT0MQ3VIb1bfrcqpQCzzCsu2U2cBUppACDAM8zBHe50SxG/tF5VDopknI9RYkmb6LjOI6n24sg5X2TtEyCjZ6TunSCleC1Mlum2Rkt0or5qYxG78y030XVrXqBrSOVKM2aa0Um6MFMC8/0ZTTgUUhVWO1/itfKoVKQshyOiayieDraNHPOe+6vvXOPcBD4c8XuOPkFGy3xg7fRHPEF5Jgh0iJaOyWoEDsRt/znbAHzTdf62ZvovXoqwVTzYZZWnRV3pTpBizSClna0ZrIaNryiGB+NN3xWoTdWvmlOzdTCAxSusNa1MW1dJqTh3iNzGbIqXaptwsiY6pt2g/4A9s6kOjvNool/OxqupvzqmAmv/SGuMq+jJuc7+xGS9AShdq7kM9WkjH0yjlS40sQZ0lqazOKUHsxh8Zo5Qy+VVRLblESnHblOu9v6ApUnKk7+I1MZeLjeOv07cu449RTq1lpLQqUozXvW8uHAfQNa1ulKzx6ylVSByztZ8rNjpk6drad6CuIXb3PzJGyboxwjVSiNeNyTdJgq1NWThFChC78eea6MvTaWFaOZVmp6RGX63KqcVrU15uakYMav6LapPQAWIbKQ76pa4hdvc/MkbJCmFF86YU0xA+51Q8CYqmOg3oO/tx0NC8MONtlOvNeO30pY2+ivmmnHeIFFO6AGlFis2y6HiNP9fEElgskVr7trFbdVoxyykGNkpCiB1CiGUhhBRCDPXiopxQ25Rk06Yc0xA+33R0gW5OTM2peBRitTDLVYNSxWgwynWnJP45tbxDol/XNDR5ehil5WKVjK7VKWsgrWkqUtLS9fOEIJZCl3yz8jbl4pDHcO1De5HSh4E1nwW6WasgmjdliN3ClFKaIbwtUjJbjQjpUDwKsRp/87EVYDkl8rRYmDkH+i6tCzV2WH3/Y3Z8Q3ONGqhIadW9h1jm1FZR17pbPjl+ax8CGiUhxAuAnwb+v95cjjfSuoaQ1cacCsTOKBUrBoZ0oC8AcRrQd83NWEFtyhpm5XrMF2advrJHSsIWKcXcKDexBGD2f2tmSSB2OTXlkFZXd4jHK1KKz/ghgFESQujA3wJ/Dkz37Io8kNE1tGb6BmIXwtebsdYnpsUra7Jp/KmsOuguVpuSQ6SU0khjSYKbuoQDFOPT/6128miDJFpTDgnEXhacK66OlGqRoqNRis/YixWDqiEdJeFCNuWUUn3qGJ8YjR+CRUq/BvQBf9eja2mJlC7MGxNv+s6peDBl0ldaM30nROwojGUHo5TRNVJYm3K8I4V8qYIm6ucIAaQ1QUqY4487fdukPgM1/1fNfYitQ+okctKa6bva2o/P+MGnURJCbAQ+CPyelJYExvP17xJC7BNC7Juamur0GmtIWxMz5t5SbWJmG4UOOmbTVUejHJ+J2XzyJqhIsWaU4r4pF6sMZlIIW0Jf10U9UnSkL+MUKTaqz8CkL5s3ZYjd2q85pA7lEEK65dTis/bBf6T0P4A7pZTf8vNiKeUnpZRXSimvnJiYaP/qmpCpeUtO9F2MJqZDTiGt2zZlrZHaiNvCrEdKjQWEaVw6GkCsFma+VGmQg4MSujgb5ThGilVn+q6ZvoLYzX2nGjUrp+ScuojX+AFSrV4ghLgIeAfwQiHEmPm0eYgRo0KIqpRyTdrUpnWBVna7MTHalBxyCrpm85RXeUvxmph5hzYzalN2GL+eUfMhRuPPlaqOkYLu5JTEMVIsrY6UUpqZT3aKFCoF80iTlttZ6OHokGn2nFJilAD2Amngxw6/Owp8Cvjlbl6UG1K6hn465JScIiXNJacC8YsUm44uAJVPrOVUVvHq8br/+aJDpKQLW6QYbwWWOuBw9fgVde+w9kF1StdH1+gKe4d8aXU+VdOEckpcc2rxmfvgzyj9CHhJ03M/Dfw34FXAU92+KDeonNLpEMI3njwJamKmhZVTahp/dgiWTqzV5fUc1nHQQ01tltJO9BXEbmHmmhpygiqePV3ou+Y6HTAjJUfq3jb+vugbJUt5aR3wZyGjSQTS2SjHaO2DD6MkpZwGbrU/J4TYY/7zh1LKNXPR0rqHt2AVEDZv2BFErtaQsnEsfZpplJzouxi1GrHUZ33pesoz1RApxpy+LVXZOJhpeK4xUoxv78dK1aBYMVbTd7owHVKHtQ+xGb8TdQ3Qr1trv2nuZ4dgJh5jtxCZ3ndgKnBwoe8gdhOzeWH26V7qu3iMHeodwu3qs5TulVOLlwIpV2w85A3qQhcp9NVtdiA24885KE+hhfIWYjf+Zoc0azmkp0FOqS2jJKX8tJRSrGWUBJARUlX1u07MeNwc9xDeLVKKGX3lQN+k3eqUIHYLM19qPOQN6vSdbJ77qYz6PmIy/rpD1jx+N5YkXms/X3RxSDWHbiYQu7UPEYuUsm4hbAxD+IGMjqaJhudr9J1jTm0ZpFyjK+wtmtusQHOdUrw7euSKq3NKaU3l1GTzpgSxMspOdTpgUve4dHOBWI1fiNUOab/ulk+N19qHiBmljPAIYSE2G5M6dXZ1uq8ewjt5ixLKa6LM7zkcE9265qw+g1htylJKxzodSxIunXKmMfKWczWRS7NTorWIlGK09tOrHdJMLVJy2vvis/YhYkYp6xrCxi+Eb96UALLCg76D2IzfqaI/pQlSIv51WsWKQaWp9xko+jJNFSncIqWYbMoO5QBgCR2M2FP3qnB6tUPa75pTildOESJnlNzUV/HalJcduiQDZFzHHzNvsVh1iJRc6nQgVpGCVdHfrL7STfpyVU4JlAIrJurLnMv4VU7RS+gQj/uvWkw5OKSaRzkExGbtQ1SNkmutQjxuTL5UWUVfgD2nFPOF6XCeTmObHRcFkmGs0RX2Dk7NeKEuCTeEg1GKUaTodGwHqEhZp+q+KcfEKOcdatTAzhLFe+1DxIxSRpweN8Ytp5Rx6hIN8Ru/Q05J00R9/G68eiX6vLpTM16o12nFPafk1CEeVE4xJSurhR6pjGo1FROHVLEEq++xu/I2XmsfomaUNA/6BmJzY9xySu7JzniF8LlidRV9Ax5Cj2x8vOVlr0iJikekFP2xg3PfRzAjJWGoOq1mxCxSdIyUXB2yeK19iJpROk3Ud/lWkVKMhR5VQ7JSrq7alMF2/2PMq7tV9Kc1JXSIO31XM8rp1UY5TQXDKacWo0gxV3J2yDKeylti4ZBZiKZRas4pxKxT9HKx4jIxW+XUoj/+nMumDC2q2iEe47daTDl2CTcwHCOF+GzKbjV6Vk4x/pFixdEhc82nZ+PFEkHkjJJLpBCzTtHWwmyGe04tRpGCy6YMp0dOrd4l2qlLuAd9V46J0MOFJbDUl6dDpNhMXYKftR+P8UPEjFLajb6D2HiLpYpBuSqdJ6ZX8SjEYvxO58lYyJ4GOUW3Oh1LEl5125QByvleX17PkXPJp6Z0DT3mRskqnPaMlFzXfvQdUguRMkqu9B3UvcWIw00SC5B2C+HTA4CIxcJ0y6mAx/2P0cL0rNPxkoRDLO6/dRR8M9I1oxxf+rJUVYXTng5pM0uQ6gOhxWL8FiJllNJuxaMQG29p2aUhI9iFHk70ZTwW5rJLpACnh9Ajb/Y+sx/bAfVIyTWnBPEwyh6Rkjd9F/2x16lrD4e0ee+rrf3oj99CpIxSCrUpO1a1p+NhlNzqVADSbrJQMBfmUi8vbU2Qczjgz0KmZVV79O+/dRS6/dgOsCThVSoxj5TyJeecSloz0IR0py9jMPacy5E14JFTgsQorScsT7mKW61C9G9MvUuyg7fkFsJDbBamW6IfPCLFmlGKg1F2FrnUJOGOcz9GkVLJhb4z576zUYoHS+DLIY3x2rcQKaNkTcyK04G5MbkxtYnpSd/Fd/xuFf2gnBIDDbSmaZvKgtBjMX63OhVdV212nCOlGEWKbkbZnPtVpy0rJupDL4fUtZsJxGbtW4iUUUq1mpgxuDFuvc8A5+OwLcTEW8x5GKW0qFJxihRilFPLFyvOUbJZp+PKEkA8IiVXSbTFksRXfejlkHqzJPGY+xYiZZTS5tEFzpFSPG6MV/Foa/ouDpuSmexNO3iLbpJoiM/4XdrM6OYhf95GKdrzX0rpeMAj2CIltzZDEPnxt++QxmPuW4iUUbKEDhXX/lfLkT+BsbYpOy5MF/UZxCpSdKroBzV+x00ZYjR+56ML1HlSVXeHDCI//mLFoOoqiW7hkELkN+aa0MGrHMJp7cfo6BKInFFqIXSQBlQKa3xV3UXeQ4GTooqBWJ1TgVhFik6LEtT4HTcliI9Rchm/ZtXpxJi+86Kv6tS9l9Aj2vffckgdnRKsAy6TnFKokKoJHeI/Mfsd6Cu1KXtFCtHelMC9QziovIJjlAyxMcp5l+JRMO+/0/hjIvTwpK88135cjLIldHA3yo7lMDGZ+xYiZZQs+qosXYQOEPmJ6UVfpahSli3oq8jTl87qKzhNjHLJWegAKqdYcbr/MRF6eOdTLfou/g6pUz413cooxyB1YSFSRknHhwIn6hOztPoocAuuOQVQ4zcqUC318Op6D7eGlGBtyi7jjwGvLqVU6jOXSElvdf8jbpTrkmgHoUerchCI/PjdOqSDLVJ0mv/WIZfl6B9yCREzStaNKTlGSvHwlvKliiOnDJCSVSputyw243en71Ki4jH+6PPqxYqBIZ0T3WAa5RiPv97Nw119503fRXz8Lh3SQUWKhhSUWW2w4iL0sBAtoyTjzyvniu4T06KvqoZDmB6b8XvQd7JK2TVSiAF95dEhHSkVfRtj9WG9GbFDTkVa9F2MHVKXvn8AOgZldCpVp7WfGKV1g4afnFLEJ2bJa2KqTbniVLkek/HnSs4HHIIaf8ucUoR5dUt95uiUGGpTds0pZocjT18uFz3Ud9bajzF95xUppahQRadSje/atxApo1SbmK68KpG/MW4V7WBOTKm18JaiPn7vSNFT6BFxXr3eId5hjIZHpACxyCktF8oADPU55ZQ8hA4xUh+6UvemQ1b2ZEmiPX4LkTJKOlUMKajI+PKqbg0poR4pVGJK3xmGVJGSw6YEtKCvom+U8x7Fk1TVhu3okEEs6LtcrSGpgyRcWg5pnNWHVUeRB4AuK5TRqSb0XbigS5Xo9t6Uoz0x8x45FTUxU84hfDYGm3K5ipTOiW5QkaJ3pESkF2ateNJp/K3ouxgYpeVihYyukU2tHqPuVQ4CsYgU88WK69zXMaiiU/ai7iNO31qIllEyK/odcyrpAfUY8YXpJQnXqVJFcxE6mEYpwhPTSvQPZR1aqQC6rFJqaZSie/+9Ev21SMlLfRnhsQMsF9zzqZqMv1HOe+SUrEjJkbqPgUNqR6SMkmZGSmWnG6NpkT/oz6pT8Y6UWvHK0TVKSwUP9Rkqr+DpKUOkx++V6PcXKUVb6OGZT7XouxirD3Me5SA6FSpSp+wodEiM0rpBlx6SaIh8CF+qGlRcGlKCNf6UC68c/UihHil5eYspl0hxWD1G+P57HXCIoSIl70gx+kIPr3sPXkYp+pFivuidU4pzPtmOSBklTSpZpKO3AJH3lvIeDRkBNGkpcOJJX/oxSmphxlMWW88pOdF36rtxLByHWHjLXkbJKgcpGfHMKZUqBqWq4b72MdTcd9r7Un0gtEjfezsiZZR0U33lyKtC5L2lnEdDRjDpSzdJuKYrwxThhbnkccAfKKPsXkAYfaOUL1XQBGRTDsvSsIxSnIUe7vSdFnOhh2c+kbpD5pi6iIn60EKkjJImK1RlfOm7nFdOASun5iL0gMgvTCtSGnaRhOtS8epxrdNaNvveCeFQ8lCj71rl1KI9frdIyRd9GfGxg3ONFqi1X3ZjCcA0Sku9urw1RbSMkmEl+mO6KXvlFGgRKUDkx7/cMlLyuP8xiBTyRXflZb1O6TQ1ShZ96UbfRbyjhTX3h92oa8N0SGO69u1oaZSEEK8TQtwuhJgRQhSEEAeEEH8shMisxQU2XIusmq024nlj8l45Beo5NcdkJ0Q+hF9ukVOqRYpO9z/dD4hIj9/r2IoafeeaU7GEHtH1lnNeRtkPfRlh9eFyoYVDZpaDxDWfbofLDGjARuAW4MPAPPBs4APAVuDdvbowJ2heChRQm3IMvCU3SXht/J4TM7rjzxUrpDThnFOhPn7HhRkDXt3r2Ar/OaVojr/WzcPVKFv0nQO1CY3qw8xAby6yh2hF3wlDFY7H1SG1o6VRklL+Q9NTtwghRoDfFEL8lpRr55oIw6KvXDbl7DAUo+spWhNzpM+5eNSi7zxzaoXFXl1ez6GKJ11yKth5dY/xR/z+u+XTLPqu6KU+g8huTLVuHi3G726Uba12ImyU3Og7TZbNfLLH3F8+2avLW1O0m1OaAdacvtNk2Zu+6huB4mKEQ3izIaXLxBRG2Ww1Ek+hx7LHUeigcoqqTs3NKYm2t7hU8Er0q02rYMSz91+uRT7RGn/RbfxZk76MqFNi0XeekVIrliTCLJEdfug7AIQQOpAFngH8NvD33YiSFhcXmZycpFwut37xBe9B7K1yZV+eRx99dPXvx38KrnkWPPqI0u1HDBcOlPm/r97GicNPMpNOs3nzZkZGRmq/16Tqku26KWeinuwtu2/KUtZySo6yWIg8r75cdG9GW+vo0CpSiuqm3CKfSLVMFc25GTPYjFI0mYLW+dQW+fSIO2R2+DZKQA5llAA+A7y30z++uLjIqVOn2LFjB/39/a60TQ3TOsvFKstDZ7J1tM/hCqdh4Qhs2Qv6mgdyHePEwgoDyyUu2D7CysoKx44dA6gZJmFUzGSnV6QYzU0JrES3d6LftU4LIs+rLxUqrvRNjb5z25TT/aClInv/WxVOU5v7blGy6bxFdPy1FlsuOUVhsgTuyuNoz307goQTzwNeAPw+8BrgY24vFEK8SwixTwixb2pqyvUDJycn2bFjBwMDA60NEoAEiUDisilZ0ZHbjQs5qoZEEwIhBAMDA+zYsYPJycna71UI7yELzQ4rTzGi41eRgnM+zdqUK6RalAREM1KUUpo5JZfx1+grlyUrRP3+RxCt1GfKKLm0mILo03emHF7TnPdBYZR9lINEV31oh2+jJKW8R0r5IynlR1D03a8LIc52ee0npZRXSimvnJiYcP3McrlMf39/gMuVDQ+roJlettm8MWowDNBtk7K/v7+B1hRSeYuuBXTZYZQCKZoe07JH6/4afdWqo0dEN6VC2aBqyJb0nWtOBSIt9PFD31WES0cDiL5R8sonosphXBuygk19mO/NBa4h2k283GM+ntnpBfiKkCxIiRTC1SYhzAVrRNQoSYndUWr+bkSrArqIUxg5z4r++smjcVRfLhW9RS51oYPHesmORnb8VuG4V0eHqlc3E2vuR1R96plPBDDK3uUwNaMcTabAjnaN0tXm49PduhB/kIDAVV9Ri5SiSV9VDdkQKa2C0aLViDUxo7owC+69z+r0nYf60FJfRhAWfdVSEl71WLLZ4Ujfe/Cm7wzhUQ7RZzlk0Rz/kpdDRj2n5O6QRXv8dvjp6PAdIcQfCCFeKYR4hRDiz4D/DfyblPLJ3l+iDVKaOSUXWDmlkNF3H//4x7n22mvZuHEjQghuvfVWx9dVpcopOUJKs6OFy8m7EOlISUqreNK795nrybugxl8pQKXUo6vsHVrSVzZJuKtTFmWjbHYz8WozVCXlTl+lskrcFMG5D6ocxMsoUS27N2SFyEeKdviJlO4C3g58CfgicB1wA/DW3l2WG9QNcc3lhZS++8xnPsPs7CzXXHON5+sM6REp1boke9B3EfYWV8pVDOmdUwBUQ9YYGmVLfdXKKLmePAyRFjrkihV0TdCXdtmSjDJVr0gJIk3fevb9QzUOqAoP+rK29hd6cHVrCz8dHd4PvH8NrqU1WkVKmhUphYu+u/3229E0jYceeojPf/7zrq8zDOmqvrE2Ze9IKbq1Gn7oG2ghdLAvzMGN3b7EnmKpRfFk/Th0ZZRTTnqHiG/KgxndPcdcLVMVHjVqEGmjvFxonVMyhFedUrSFHnZErMJUKumrW6gkNEBb00jpxS9+Ma973esanrv11lsRQvDQQw8BoGmtv2YpJVUDdLeUkmHblGJYq7Hc4tgKe06pZU4twuN3azFVF3p4UFjZ6NaptYoUVE7JQxIO8R5/tYwhPNoMnWb0XXggTaGD12s0LXQ5JT+QEiQekZJpaD37X8VgU3ZvSGoTOsSQV2/VYsoafxWP4uHsMFRLUC704hJ7Cq8D/oCa0MHVIENkjVK9Rs3t3lcBqYyS2/gjTN03I1pGCZVV8qwPE3oki0er5qB0D/oCaFHVHl31XauzlKzzdDyFDn3RjxTd6UvllJS9nJK+UfUY0fF70letIgWILH3XMp9qRslS81CeRpglaUaQNkNrgj/7xsM8ctxlYpWWqTBDRRynL+1SRFjOK4ovdSLw375w+wh/et1Fgd/XDRjmZHOPlOr0jSuFoemRLSC1Tt119xZ9SMIjLItdKlbIpjQyLsd2WE6JgfBBXy7CkHvRehjRmr4rY2it6LtoSuJbNWOt3XuRdnfINB3S0T4lwELkIqXWEO4dH0KMlpGSuSlLr6p2ML3F6Clwls3iUT91Si1zahFcmMsFD/oGzER3ChA+kt3RG//iStk9nwZQrSCFRz4NItv7callOYBplDQP5S2Y44/e2m9G6CIl10hFSjhxH7PaOHP6Rs6eGHJ+3exTUCnC5gt6d5E29PX1USo11sXMzs4G/pyWkVLVCuE96CuIrAJrqVXxqCX0kB6y4Ahvyp7HVoDKqWjq93GVxHsbZVPo0NIhW1J7RZBOMeuM1oXT6veIlDtLAOr+R9Aha0YEIyXROqe0hpLwnTt3sn///obnbr755sCfY601d/WdzSjFUIHU0ihV7eozl/Gn+8wCyugtzNY5FRUpAN5tliCSG9NiocxIv0ekZJSRmg9JuFFWBdQRQr1w2lt5GVeHtBmhi5RcYVmiVh6Qpq+pJPz666/nU5/6FO95z3u49tprueWWW7jpppsaXrNv3z4OHjzIkSNHALjtttuYnp5mz549XHnllYCfnJJJ32keBXQQWV59caVMNqWRdSzAwZ8kHiLrLbZqyIlRQZqRUtyakpYqBoWy4X5sB0C1gqENup8lBo2RYjpIo+f1RevC6fra9zTKfSNQiD59F6FIqX4zPM8WFKYkfI1auF977bV86EMf4sYbb+T666/n0KFDfPSjH214zcc+9jFe//rX83u/93sAfOADH+D1r389H/tY/fSP1jklNXGVLDR+Ve0tPWU/ve8gsuNf8jq2AmqRAuBR1W+p76JllJdMOXzLSKnl3I8mfemnQzqA1NItHNJoOmTNiE6kVDNKLeqUrFZD0qj/u8e44YYbuOGGGxqesxvOT3/603z605/2/Ay/OSXRir6LaLJ30UdOATA3plbJ7ugtzOVimWEr0nFCtWIzSvHKqbWkbkHdf83jLC2I7PhrNWotji2RrYQOEXXImhGdSMlG33kGQSFtNdQKVSkRQrg3ZLUiJc1DFgpmTilaixL8qK9Mb1H3kVOLoLfYus1MBWo5JZfxp7KgZyO3MS1akVKL++9LEg6Ru//1GjXvs8TQWqkPRyO59psRHaPUECl50XfhbMraCoYh3ak7qPHKeBXQgVqYpeXIjX+pUGlJ3wAILR27qn6ror9VnU4tUmqV7I7Ypuw/UkpTrkp3+j6iObWlQoWMVz61aq19H22Wyvn66yOK6Bglu9DBM1KK5umzVVkP8pxfYE1MH5ESRO5Y8MVC2R9956tWI1qbcrFiUK56nDoLavy6MtqtuxpEa1NeXPGRU6rWjbLr8CNqlBYLZUY9HTJz7utpb4c0wh1N7IiOUQqaU4pYpNA6UjLHE1MKY3Gl0rJ4Eqxkb7xazdQihRbqM1oJHSCSRtlvpCQ1NT9cI+WItllaWPFnlERLhzSaObVmRMco1WxSq5xSNE+frXodWwF1+k73UasBkVuYS4UyI/3e9BWA0P3k1JbWTH3ZDSz6VJ/RShIOkaQv/Y2/gjDXduvi6WjJopVD5t33DwDd4+gKiHRHEzuiY5QC55QqPb+ibqIqJakWR6EDir5r5SlDpDamQrlKsWL4EjrQSujQN6IckgjRlzX6ylMSbouUYmeUKggBQ24d4kHRdxZ92UroEbFaHT+FwwBCz3irD2PSKTw6Rsm3+i6iOaVW9F3VSvS3koWaFEaEFqZF33h6i76FDtGLFBd85lTqOaV4CR0WV9RR4K2YAmHRd17j7x+L1NwHH/SdvRzED0sSsfvfjOgYJXvxrNfLhDmkiOWUWtN3ZqSUahUpWEZpvmvX1mv4Kp60+n/p6dhRGJZR8s4rVBGtIgVQ9z9C9x5M5aVXlCilyinpymnxzKn2jcHKfFevr9doWQ5h2Oe+l0MSzZxaM6JjlBrUd170nTDPVIqOUTKkxJAS3YdRapns7B9TjxFamIu+Et1lEDop3ePkWYhkVwNr/C1zan6EDv1jauwRmv+tlZdqLForoQOYkdJ89y6ux5BSslio+MqnanqLhqwJfbfWaFTfebYa0lKRWpSW5+dplKr1RL93/6sx9Rihhekrp2LSV2ldtE70Q6QonEVfkVKlFin5u//RGf9SoVWkUM8nQrwipVypStWQLei7ek7JVzlIhO69E6JjlEwjJMy8i/eR6DrIcAgdTpw4wXvf+14uu+wyhoaG2LVrF29729s4fvx47TWGH6Nkq+r2XJSpDKQHIrUw65LgVon+NLomvMffv0E9Rmj8Cytl+tIexZMA1QpC87Ep1yLlue5dYI+xuNIiUrDV6EEroxwt+tK3yAVaU9fpPkj1RWr8ToiOUcJG3+HnSPRwREp33303X/3qV3nTm97EN77xDT784Q9z55138rznPY/lZaUQq0VKnnVKJn2XSnsnekF5ixGamHVJcIuNSU+R0j2Og4f6phyl8bfKKYBK9Os+6KsIRspLxXJrhwQ198GHUY5QpOBXDg8Wfedj7UfIIXNCdBqyyjp9Zz5h+3cTNB3KJeffrTGe//zns3//flKp+lf9jGc8g/POO48vf/nLvO1tb6t3CPcRKWlamkq1xdiitjB9eYtl0BR9Vyj72JQjFCm0VF+Bou9SPjo6RDBSbFmnY9TVZ+DDKBcWwTBatEgJBxbyPqhbi75LpakaRe8PjFhOzQnhv2tNEH4iJU1fM0n4i1/8Yl73utc1PHfrrbcihOChhx5ibGyswSABnHvuuQwMDDA5OQkEzCml0t6eIkTOW1oqVNA1wUDGm75CT5PSNG9ePZWB9GCkxr+w0qJOBaBaDkbfRWRjklKyVGgRKZlzX/MbKSEjU0BbE7n4iBS1Vr3/QDklEZr7ToiOUbI6NPjNKRlrd6ZSUDzwwAPk83kuvPBCwKdRsofwXpsyRM5bstRXolVDWi3VWugAkRx/60ipWtuUfdF3EdmYcqUqhvShPAR1qjCt2iyNqceIjL/e9691Tk1LqfHHSejhhPDRd9/+73DywdXPm8ccD+n9nFUBPaO7n0JbLUG1CJkhXCk+J2y9BF75l21dtl8YhsHv/M7vsHfvXl7xilcAwXJKWquOBmBOzAe6cblrAl85FVN9l9I0700JIrcwF1bK7N3scZYSmMWjKTTRok4pYpGSVaPmHSlZc9+H0KE2/qhESv6UlwBaqn6elqsmpn8DnHqom5e45gifUWoJQYs4yWasPPJO64QbbriBH//4x9x2222k0yYdISWaEN7Fs1UVKaR0PXaR0lKrA/6gllNK6cJ7UwaTwohQTilf9s6pQN0o65p3sjvdr1rtRMQoL674PLYCakIP70ghWsXjVuF0q2NLAHQzUixXDfrSLlapfyxSc98J4TNKbpFKbhoWjpAbOZdD82XO3TLsfmNW5mDuIEycrxZpSPDxj3+cD3/4w3z+85/nqquuqj3fspsD1CTRqVaSaFCRQmm5oTVNmLHYqk4FzJxSipQmWiuQ+sdg9qmuXV8vYRiSpWLFF32HpsZfjRF9aW3KY/0Z9xfVikfjR18urqhztFK6RybFihTTPjp69G+I1Np3QnRySlh1Sj4ueQ2Pr+jr66NUalTDzc7Ornrdl7/8ZX7rt36L//W//hdveMMbGn7Xsu8d1BpypnTNR6Qwph4jQmGoRL/fSElrvSlHiL5bKlaQsoUkGGo5tZQmfNK38926xJ5iPq/WzthAN4UORMYoK4fMx9xHoFs5RT85tYisfSdExyg1ScK9OzqsXVPWnTt3sn///obnbr755ob/33rrrbz5zW/m3e9+N3/wB3+w6jOqRosWQ2Ae8mYm+v3kVCBCG1OZDQMenjLU6Es1fh9GOSIUhq8D7qB2yF+6VZ0WRCpSmjcl0Z5GqZZTCdDRIiJz35fy0nRI0+Ye4SunGJH574Tw0XeusCIl+/9coK1dpHT99dfzqU99ive85z1ce+213HLLLdx000213z/66KO89rWv5fzzz+cNb3gDd9xxR+13ExMTnH322VQN6R2+Qz2npGmqP6UX5Rchb1FKyXy+zKjXpgTqXvqRhIMaf2UFKkV1nEGI4asZK9hyij7p26UT3bnAHmN+xYqUvOg7yyhlAMN7/JlB1WYsIpHCos9yACufCD7oO4iMUXZC9CIlk77z7uhg2to1OFPp2muv5UMf+hA33ngj119/PYcOHeKjH/1o7fd33nknCwsL3H///Vx99dU897nPrf188IMfBJTQoXWkVE/0g88QPgITM1+qUqoarSMlo74pt6QvIzR+f33vqoA0c4parCTxc/kyKU0w6FmjVm9ICi0k4UJEqtXQYqsO6VDLJ6cDrf0kUloDBOx9B2tWQHvDDTdwww03NDxn0YsXX3wxb3/72z3f74++qye6QXlLroKdCEVK8+amvKFVpFQtQ3ZY0VcthQ6WtzgHw1u6cJW9w0KQ3mearoyyL0l8NCKF+XyZsYFM6xo1QE9lgFJroxyhnNriSpkLt414v8ik7nVf9J059yOw9t0QoUjJfKzxd62Or9BC0//OC1JKf0IHW+83aNFqJkLe0lxO0TejXuorqEeKmh9J+Jh6jMDCrNWp+Ej0K/rS5/iLC5GY//P5knc+Cerqs5rQIT45NV8iHxt1Dz6OLoHIGGUnRMcomTVH9gokT4SoKasXDD9976BBEg74O1MpAgtzwXekVKkZ5YrRotVKX3R49Xqk1LpOx6LvfEVKEIm8ihK5+FAeAqmUVacTj0ipVDFYLlYYb0ldN9J3LQ95hEg4pG6IjlGSEoTwFSgBkTlTyVeLIbBJws2J6Xl8RRZS/ZFYmHN5H4luaIiUID6y4IWVMppoVTxZP7bEX/HwmHqMwPjnV8qto2QzUtTTPtrsQGQaEtfk8IM+xt/Akng4JXoaMsORuPduaGmUhBCvF0J8XQhxTAixLIS4WwjxprW4uEYEjJRCdKaSFyzjkmpllMyJmTZD+LjIgufyAXJKetqfUbbnlEIOSw7vmVOp0Xfm0R1+1HcQCafEH33X2Put5dzvG4vEvbfmvu9IydwjfAldInDv3eAnUvo9YBl4D/Bq4Bbgc0KI3+rlha1C8yF/rUIlLbUm6rtOETRS0v1EChAZCmPB9BZbS8KthenDKNcojPkuXGFvMZcvsaGVp2yLlNKa8JdTgUg4JYHoO7+R0sC4GnvImZJZM5/qa/xayp8kHCJVp+cEP+q766SU07b/f18IsR1lrP620wuQUnp7ifVXujdgdUJE6LuKh1FqMLzNkvBWE3NgPBITcy5fZiCje5+6CjYKwwevrumQjYYseDZX8uEp109e1TU/XdLH1WPI73+hXGWlXG1N3VqRUtrqEt5q7m9UpwoUFtQ6CCks+q6lU1LLp/qQhEPkDvlsRstIqckgWbgX2NzpH0+n06ysrPh7sTTpO3PvbjUv1fEVldAeX2HBapnjRN+trKzUmrZaxaNpP7wyqIWZn+nqtfYCvro5QEObIfCxMPtHQ78pA8zlymwY9FE4DLX737J4eGCjesyvbncVJtT63vmk79JmIXRL+q42/nDP/1nLKLU0yiXQMzWWwF9D4vkuXOH6oF2hw/OARzr945s3b+bYsWPk8/nWdJwZKfnOKmlWAW24o6WKIRE0RkpSSvL5PMeOHWPzZtP2V8uqTsVPrQJEyCiVWnczADNSyvhrtQIqWgj5pgxqYxr3k+gGW52Sj00JQn//ayIXP+UAqEhJCD/0XTSMkq8WS6DGr2dsLIEP+nYl/HPfDYGLZ4UQLwNeA7yj0z8+MqKKxo4fP065XPZ+cW4aqmWq04JTCwWK02lOeSmWSjk1KWcfCXW33Pl8iZVSlf1Ljd3M0+k0W7ZsqX1HNUm4n0Q/mEZpNvTHQqucil+j5LPVCkTCKEspmcuV/HWzAP91WnpKUTg5J5IjPJj3LXIxGx6bOcXW1HU0jNJsrsRARnc/7cCCOffrHR18rn0ZMOUREgQySkKIPcDngK9JKT/t8bp3Ae8C2L17t+dnjoyM1DdeL/zbW2D6CWbf/gN+5oM384HrLuTtl5/p/vonvgtf/Xl4x3/C7ktbf/464Tf/9R72n1zke7//DO8X1rpEW5uyDwpDVlURpeU5hxDzK2W2jbU4XkRKk8KoS8Jb0neDm2Dm8S5dZW+wVKxQMaQPo2SThPupUwI1/pBvypZRailyMYtnrWihpdAjIkZpLu/DIQE199MDAdb+JrVfFBfrop8IwbcLLYQYB74NHAbe4vVaKeUnpZRXSimvnJiY6PASTRjVxkR3TEL42ZwP+gZWJzv9eou5cI9/Pl9mzM9ZQsgGCqM1hbMp9PSd1c3CV6Ibave/ZaQEkYgU5/3WqNUiJd2f0CMiObW5XBCWIOOvzRAohwRCHym7wZdREkIMAN8EMsC1UspcT6/KCbVWG/EySnP5UutFCfXiWdNbarkpD4Z//IYhmffjLVr0ldklHPwku8fVYWflQheutDewJMHjLYUO9Y4OabOjRUtYFE6I4bvvoZlTQQjSutZ67mcGVPF4iOc+KOWpv0jJrFH002IMlEMGoR+/G/wUz6aALwF7gVdKKSd7flVOqFX0B6CvIPQ3xpckGIJ1CYdIjH+pWMGQftRXtpyCH0k41L3FfHi9xXpOxW9OSdWptZz7oIxyiMcOyiHL6Br9fnIqmpojuuajIS1EwigHou/sQoeWa9+UwUc0UvKTU/o48Crgd4BxIcRzbL+7V0pZ7MmVNaNaacwptNqU0gOQ6gv1piyl9Fc8CXVJuF9ZaASMkn/6xoqUMv5arUCjtzi6s5PL7BnqkZLPnJKZ7G6Z6AaTvpwJdbJ7wTxHq2Wdou1o77Qf+g5MoxzeuQ8mfdfKIYNapFgvHPfrkIV7/G7wY5ReYT7+tcPvzgQOdu1qvGBUIJVB0wSaH1moEEoWHGJpZK5UpVyVrekbqEvCazmV6EdKNUmsn0POwDTKQXNq4fUW54IUT4JZEuCjTgnU+KslRWFmhzu80t5gernERl/51FLNKOl+DjmE0OfUKlWDxULFn0NqO+DReq8nBsLPEnihpVGSUu5Zg+toDaMM2iCA2f8r+iF8LdEdJIT3uylnBk1ePbwTcyanguyNQz4T3UEk4RHwFmdzJVKaYNirtAFs48/6q1OCRqckpEZpJldk05CPk4GtnBKYknCfa3/+cIdX2DvU82nt0Hc+cmrpgVA7ZF4IbwFLM5pCeH8KpHCH8L7pGzDHn/VPX0HojfL0shp/y42pgb6LT07NErm0pq8so5TxV6cENgVWeMc/s1xq7ZBALVIAYqM+rFPXfliSShN952fth78kwA3RMUrmyaugkp1xCOFn/dI3sLpOJwZGecY0Si03pgb1nUlf+jlTR+ih9hZVOYBP6hZqkaJvhwRCfv+LbBz0ESlV65GS7zqtgY2q/5tFfYYMszmzQ7jvtZ+qpS58O+QhnvteiJBRqntLaT1ACB/iRembvqsVjwao04HQj396uchARmcg45O+sqsvW21MmhZ6ozyX8ysJrkdKVqK/ZVsuS4EV0vGvlKrkSlWfkVI9p5ROaZT8bsoQ2v6Hc3773kFt7UOA1MXgplBT916IjlGy0XeBQviV+dD2v6vRdy0lwbbiUb+SeAi9UZpZLvqnb0Btyn6Lh8GkMMK7MH31vYMGo6Sb97+lTxLySMnKJ27yc/+NSm3tZ3WNciX6kaLFErRWXhqqM0stp+Z370vou97DLB4FK4T3m+yVoe2YO5cvoWuCYa+jsKEh0R9oUx7cFO6cQq7kn76BRqGDXwonxOOf91sO0CD0sO5/i/FnR1RtT0iNco269XX/S7U6pXRK+GdJILQb8/SyT5GPrUYNVKTkyyEN+dr3QnSMUkMI67OAzupqkFufet9WmFlWxXNay1NnbYnuoJtycaG+qYcMU0s+1Vf2TTlITm0wvJFi1ZD+C6ft9J1fBZYQoY6UfSsvoVHkpGuUAhmlcBrlqaUio/1pf+eIgS2n5rdObRzKOSj7PBooRIimUfIbwg5tUY/L4TRK08tFJoaDRQrB6Ktw5xVmciV/9E0DfedTEg6hpu9mcyUMCZtHgtz/+vh9U1gh9ZZ9Ky9htVHyM/ZBs+dm5Nd+3SEB5ZC3FPlAvVYpgmKHCBmlMpiHfKX9KpAGzbOIQjoxp5aCT8yMuSn5W5jhHb9hRgq+PGW7+s5vqxWwHd8Rvpzi1JKKFCb8RopCa+h/5ovCGtwEualOLrNn8E1fQa3FFkAm5VPkNLhJfWchnPtgsQRBHLJ66sK30AFC65R5ITpGqVKsV3X7jpSsTflUDy+sfUwtFf1vSgC6qmlJ68IfhRHiSHF+pUzVkMHoOy0ofbcJlVMMnwJrytyUfTslpqecSaklW/TjlAxvheWTbV9jLzGzXPKnvITG8es+zlMCdfL04ERo176KlPpav9BopO/SvkVe4S8ed0M0jJKUZlW3WsBKFunjxvSNqveEcGJKKZlqg74Dc2H62pQsoxS+8c/UPOVg9JWVU2vZZglC7ZTUIiW/99+2KYPPSGlos3JIWp7qvPbwrbyE2rEtoDZlX2OH+vhDiOllv9R13SED/NepWZHScjgjZS9EwyjZEt2gZJG+NiUhVLQQQgpjYaVMuSrb4pVVrUYQ+i583rIVKWzy280Cggsdhraqx6UQjn/JkkT7jZTq9BXgM1LeCpWCOuwtZPCtvISGSMl3TgnU2g+hQ7JSqrJcrPjPp0G9HMZvQ9phc+6HcO23QsSMkpqYvg76sjC0OZQTM5in3GSU/BYPZwaUNDiE3mK9m0Mw9V0goUPII6XBjM5gq753sGpTBihXAgh9lsI3ft+RAjTklHyr78A0SuGb+9OBqNtm+s6nJDwzCJnhUN77VoiGUaqYm5Jd6BDxiRk40Q0N9J2vnAKo8YcwUrDou0AUhnn6phA+hQ7DIY6U/FK30KA+q0dKPsQbITbKvlsMQcP4s36FDlB3SENGX9byiW2sfd8NeUHR90mk1CM43BhfbXYgvJFSzVsKtimDpUDyO/5wGuXp5RKa8HGWEjScJwRWp2gf47e8xTDe/6VCAKNUakh0A5T8REo1Cidc4w+kvIRVkvBAc98oh07oEjifCE1z369R3hpKh6wVImaUTKFDIPpui9Lqh6wxYz1S8qHAaTZKukap4lPmHFKjPJMrMj6oIp+WWJXs9Xn6KihvMYQL03c5ADQYpWygnFI4I6X5lTIVv8pLaBB6WMeh+3JKQzr+6eUA+cQm9Z3vFmsQ2rnfChEzSvVOwYFCeGTo9PpTy0UyusZIv5+cQpO3lApglIe3hm5RAkwu+uzmAA73PwCFMbQ1lJGi73IAWBUpgM/i2b4x5ciFbGM6uVAAYOuoD4cMGpsxp3y2WQJbSUS45r/lkPpuRgsN6jtfymMw53746MtWiJZRStXpq0DJTgjlxJwYzrY+SwccI6VARrm0DMXlNq+0Nzi5WPC/KVUb6TvfslgIJa9erFRZLFTaipQCqe8s9WnIjPKpJWWUtvjpZgGr6pTAb6QYzjq96eUiGwbqoh1POLQZCsQSlPNQXGrzStcH0TBKFeVZ2CkM/7LQcHY1mFoqsilIohsahB6BhA4Quv5/pxYLbB3xa5RKylM2DbjvNlNg8urhckisFjvBhA5N6rugyf4Q4dSCZZR83H8pG7qEW0bZV6QYVvpuqRSMuoQGSbj/fHo4c4qtEA2j1Fw82pZRCteNCUbfNKnvAimQwicLLlcNppdL/jYlaPCUIWCye3iLakwZIm8xUKIbGuuU9AAdHSCU9O3JRWWUNvvpaGCtfa2RvvQVKWVHINUXuvFPLbdDXduFHgEiJQgdfdsKETFKVqSkbmQgozQYTqM0vVwKtilBk9AhuvTlpLkp+88p1D1lCNAlHmwFtOEZfyCRC5gtthrpu0hHSotFNg5mamPxhLX2bSwB+KzTEiKUXR1OLhTYFmTuQ9MJCUmktP5wkEQX/S7KzABkR0PlLZSrBjO5gHUq0CQJj65RqiW6A9F3NqMUiL4LX1eLU2ak0A59F6ghL6iNKT9Tr/ULAU4tFvxHybUaRfX6miQ+orLoqiE5tVhg21iAuQ82+k4Lpr6DUI3fDyJilKyciplTMiOFlkdCWxjZDovHe3RxwXFqsYCUsN13on91CO97UxrYqPIxIZqY1qbcGX0XgL6CUI3/xMIKuibaou/S7URKEKqc4qkgIpeKmiu1tR90/MNbYelE0EvsGWaWi1QMydbRfn9vaFLfBer9Z6kvQ+SQ+UE0jFKz0CGtDsby7S2N7oCFo724srZwwowUto0FnJjtFM9qGgxvC5VRDiwJtkmiISiFEb5I8cR8gS3DWX81WuCsPvPrlIxsV4+L4dmYVaTk0yDXjJIVKQU0SqM7YeFYaGTRtbXvO1JspC8DzX0hzFql8Mx9P4iGUarRV2ZOKfDC3AGLx3pxZW3h+Lw6DdI3r+zQ/8p3ohtCN/5TiwUyKY0NA+nWLwa1MaXq35UepE6tf4N6b4iM8omFgn+HBJrUdxZ95XNjGt2pHheOBLnEnqFUaUPkAqvUh4HWfjkXmq4OJ4I6ZE1GOVCNJiiHNESRoh9ExChZkVJT/y+/E3N0p+oUbnkd64yatxSEvhOaOiMGyAQJ4cH0FsOxKYFSX20Z8VmjBeq+peqedTqILFaI0BnlEwsr/u89NERKQohgQpeRHeoxJOO32mv5zie6REqBWBIIzfhPLAR0SCvNQo8A+VQI3dz3g4gYpdX0FQSYmCFbmCfmVxjOphju8xkpNOVUAqkPQS3MxePgV7HWY5xcCFCjBKsipUCtVsA0yuGgb6WUnFgosL3NSAmC5hVGVf+/kIz/ZJAaJbAJHZrVh37zyVakGI61f3JBsQTjfo5sAWWUtHTNIQ1UOA51+jIka98PImKUGotHa7Ua5YC8ekgmpqJvgnjKzZtSwBB+dJcybCE5VyqQ+gpWR0q6zyOhLYztCs2mPJsrUawYbURKdQcmkFMiRKiM8mRQkUtTpJQJ0mYJbJFSOMZ/wpSDB2MJ6t9V2uz76VvkNbpLMU0ha7PmhWgYpUojfZdNBw3hTW8pLJHSQoFtftU34LgpVQyJ4ZfCGgnPwpRScmqxGDBSalyYgegrUAtz6WQoZNF16jbo/e/EKQmP0McqnPUvcmmShKeCSsK3KPVpSBzS9liCukNWP3nZ59of26Ue58ND37dCNIxSU5fwtoQOEJqFeWJhhe1BIiVb8SS0w6tbFMb6j39xpcJKuep/U4JVkVJfWqdQ9tklHczxS1haf7FDYJGLUQVZ7ZC+3Rkah+zY3Ap96YAiF2i/zZKmm+rTcIz/xGLAfGKTQ2ZJ4n0LnUImdPGDiBmlRqGD7xuTGVAqrBBMzGKlyvRyKZin3OQtBTq+AGwTc/3Hf2QuD8DODUHH37gwA6kPQ2SUrUjBN31bVkaMdFOkGCRSGjGFPuWC//f0CEfnVti5YSAYfQWr6LvARjkEc98wpIqUAq/9ukPSvlFa/7nvF9ExSlq61pAzsPoO1MIMwcQ8GVR5B+bErE/kQMcXgDLI6YFQTMyjNaM04P9Nzd5iWg9olMJDYRyfL5DWBZv8nrratClDm5EShMIpOzqfD+iQWONvFDoEL4lY/7k/kytRrso21n799X1mjaZvpqBvLFRCFz+IhlGqlBwjhaLfg+7AVKCt/6I8Pq+MUiD1VbnQ6CkHVSDVZNHrPzGPzinPf1cgo7Q6UgxE34WIvj2xsMLW0T40v4WzFTNSsie728kpQSjGf2R2JXiUDLXx96UCbsoQGvVpYDk4rKKurXy6b6NcE7qsv0PmF9EwSk2J/qw5MYNFSuFI9rY3MQurNiVoh8JY//Efmc0znE35O9zQQlOk1Bc0Ukr3qca8IViYJ+YDilwsyi1df0+g88QgNJHSYqHMwko5mEPSVA4SeFMGxZJUS+uuQDs8q1iCXeNBxt8099syyolR6j6aNuXAdUoAY7uhMA+FxS5fXDAcmVVGKVCk5Dr+gN5iGIzS3Ao7xwPkFMCRVy9VDP/qQwiPUZ4LSl9ZkYJdEi+CO2Sw7uM/ZkbJwajbxkgpm9IQAoqBhS6sO31rGaXdQYySg8gHAhqlEJVE+EE0jVI7kcL4mepx7uluXllgHJrNsXWkrza5fKFcaPCU+0yjVPBbpwUwdobq/2YlztcJR4NuykZVHYfdkFNqwykJgVEqlKucXCxwxvig/zfVNmXb/Q8aKaaySho9f8j/e3qAozWjFMQoW5JwtTELIRR9G2T8G/aox3Ve+0dm82wayjCYDcISrBb5QMBIcXSn6hRfyvt/zzoiGkapvKIS9SbaEjpsMI3S7PpOzMMzeXZvDOApgcordOotWeOfW7+NSUrJkdmVgPmkxjYr0CaFMbYb5g+va2POo3N5pIQzgtx/B/VdXyqgJB5g/CyYPRjsPV3Gkdk2lZdaqtbRANooCbCM0nqv/dl8MOoOHHJK7dB3ptAnIhRehIzSavoqkLcQEm/p0GyeM9qamPWF3J9RE3MlyMQcP0s9zj4V7G93EbO5EivlakeJbmgzr7BhjzLu63iExaEZk74JYpQc1Hf9GT3YvQfTKK3fvQcVKfWndf8tdkDlglKN+dfARjkzoGqV1nvtz+SDUXeg5r/uJPKKnkPuF76MkhDiHCHEPwgh7hdCVIUQt/b4uhrRJInOthMp9Y2os4XW8cbkSxWmlorBPGVYZZTrkUIb9OU6bkxHLOVdUE4dOo+UNp6tHmef9P+eLsMySoGcEgf1XV9aC3bvQW1MS8fXlcI5Opdn13h/8Hyi3mjE2h7/Os79ctXg+PxKG0apCzmlEMz9IPAbKV0EvAp4zPxZW7hESoFyCqAm5jp6S7VE58YAOQVYxSv3Z9T4A3nL/RtUc851HP+RmvpoHSKlEESKh2fzDGVTwSIFB/VdX1qnUAoaKVn07cFg7+sijpiFs4HQNPehDfoOzEhx/eb+8fkVDBnQIQOHOqWAfT+hvvZn4mWUviGl3CWlfD3wcC8vyBHllUZJrN7GpgRqYa4jr96WpwwOyc42vCUh1p3Csbo57AjUIbsx0Q3tymJ3qQLsdVyYh2Zy7G5HeQir2ywFqdGDulFeJ6dESsnhmRy7glC3YNYoNhrxbFoPJnQAGN+jTmAt5YK9r0uwHNLga7+5RlPN/UA1mkLA+NnrTt/6hS+jJKVc36qzykoDfWedKRPoxoCKlBaPrltjTitSCETfVStgVFZ5yhBQFgvrTmE8PZVj01DW/5Ed0L1ISdNVXmkdx39oNh+cunVQ3/WndcpVSSUIU7DO9O3kUpFcqcrZm4eCvdEpUgpaPA02oc/BYO/rEuosSXcipcD05cazY0ffrS+aOhqAySu3Q2FIQ6mw1gGHZvKM9KUYGwhA39RyCnVvqS2hAyhvef5I/SiQNcZT0znOnghKXTrklNrh1WFdI0XDkBydXQm+KTmp76yNKYhR7t+gftZp/E9OLgNw1qaARqmpQzqYbabaufewbhTe4Zk8GV1jy3CQI0sqqhmvA0sS2CEfP0uVRITkoFMvRMMoVRol4QADmVTwTXnD+tYqKU+53U25wzolMI1ydd2M8pNTy5w10YanDM61Gm15i0+tiyz85GKBUtUIVqMEjuq7KBrlJ6cVbXZWYKekscUUWJFSG3Mf1m3tPzWdY/fGAf/tpcA29xsPeNREO2v/bOWQr2NO0S96YpSEEO8SQuwTQuybmurCwXLl1SH8QEYnHzhSMr2ldcorHJzOtae8gwZPOaVrpHXRXqQE67IwZ3Ml5vPlDiKl1Zty8JziWVDOr4ss/KC5Ke9pp0ZN6A1ttqzxrwSd/+tI3z41tUx/Wg92lhAotWCTQ9pWTq1/g2pOuk5r/8nJZc4J7JCtnvtCiPaFHhCJvFJPjJKU8pNSyiullFdOTEx09mFGVfV/SjcmSPszevBFObRZTcyp/Z1dUxtYKVU5Mpfn3C3Dwd7oEClAuwWUpjR0+olg7+sCnpxS9M3Z7UZK+ur2/e0vzLXfmB47tQTAOVsCjt/BIasb5YDj33i2onDW4QiLp6ZynLlpMFikAMqJyDQ6MkoSHnDsABPnwfTai4eLlSqHZvPsDXrvHah7aOPoFqjLwiOgwAs/feeyKfen2yggFAImzoepA126OP94cmoZKWFvO4leWL0xZdowSjWj/Giw93UBT7VrlBwk0XVevc2FOf14sPd1AY9PLjPan2ZiyOeRFRYqq/Op/bVIKeD4J85XFM46bMxPTS8Hp+5AqeVWGSU9OH0FyihNPrrm9O3B6TxVQ3JO0LVv1ZRlGt/XVqQ0MG5GimvvkAZF+I1SbVNqDOH726HvADafrzblNZ6Yj08qT3lv0EjJYVOGNgsIhYDNF8Lk2keKT07lyKQ0dgSVBJeUMbNvTHUFUsD7P7ob0oPrEik/fmqZc7cMBZODw6rCcbALHQKOf/MF6nGNx18oVzk6txI8nwgqUnKg7wJHiQATF8DKrDrwcA3xhCnyCGyUyqZ8vWn8bUVKoO7/5No7pEHht6PDgBDidUKI1wE7gAnr/0KIgCR5QFRW51RA5ZQC03egvMWVuTWfmI+dWiatizYkwc4hfH+6zfFvPn9dvMWnppY5c+MgelD6xqorsRmltiMlTVPjP7W2pXZSSh6bXOKczQEdElA5RYd7D23klMbPVn3kJh8Jfh0d4NCM6vkXOJ8IKlpojpRMoYMMOoc3n68e19goPz65hBBtsAQOcx/ajJSgbpTWsf+jH/iNlDYDXzJ/ngNcaPv/5t5cmgkr0d/kLQ5kUuTLleCfN2FOzDX2GB4/tcyZmwZrZyH5hoP6DtpM9oLyFosLsHQi+Hs7wJNTOc7e3CZ9AyrCMWHllALTt7Au3uKMKfIITN2Cuv+rouQ21XepDGzcu+aRshUpBN6UpVSRcnOk0K7QZcKMFNdh/Ls2DAQ7GQBsRqnxewvcJd3C5gvV2l88Hvy9awi/xbMHpZTC5edgT6/QQX0GamEG5tShbpTWOK/0+ORScOoOWoy/zU0Z1nRjzpcqHJzJtRkp5FQ+Ta+3+9c0odSXxTacks0XqcPeltcuUrZEDoFFLqDoq1VCh06M8vlrHik9emIRXRPB6atKAZCOkQK0URIwvBWyo2ueU31icjn42MFmlFaXw7Q39y9Uj2t8/4Mi/DmlinNORdF3bdyY4a2qD9QaTsxCucrh2Xx7nrJDTgWsSKlNXhnWlMI4cHIJKeGi7SPB31zKrfKUAQazKXLt3P+aUV47Cs+KFAKrr0Dd/+zqRDe0sSmDihbmD61pu51HTyxy9sRgG5GClehvnPs1+rIdodPm89c0Uqoakqemc22ufWf6bjCbYrkto2TN/cQodQZX+k4nX64G55UtBd4aTswnJpXyri1PuWgZpcb39qXa6GgBMLgJBjat6cR85IQ67ffCbW0apczqBT2Y0VkuthMpWN7i2jklj51aYqQvxebhgMo7UPffQX0FbQgdwOaUrB1T8MiJRS5o5967JPoHs2r8bW3ME2srdHpyaplSxWg/SoZV938o26bIa2BcHeFxKjFKncGFvurP6EjZBq8MsOViOPUQGGvT0u+R42pTPn9rGxOzpKgfso3vbetMHQtrnFd55Pgiw32pYOcoWSgtr/IUwYyU2tmUhjZD//iaGuWHjy9y/raR4Mo7MCOlpnvfrtAB1txbnsuVOLFQaNMhsTblRqM0ZJ7c2tb933KREjqtUV7loWMLAFyyczT4my2WZJVRbnPug6m+TYxSZ9h+Ofz8Z+qFnyY6WpjbL4fi4pp1Nnjw2ALD2RR7grYYAuUpa6lVCqzBbIp8O/QVwLbLlAJtjXrgPXpikQvb3pRXq6+gg4UphNqYTj4U/L1toFI1eOT4IpfuaGNTAigurfKUrbmfa2fuj5+lRCMn7m/vegLiUTNK7ixSWk1fQZuR0vYr1OOJ+4K/tw08dGyRvrQWXOQBiiUQuuPab2vsoJySqQOqr15IEX6jNLwVLnwN9I81PD1gNiXNtxMtWBPz+L0dXpw/PHBsgYt2jASvZgczUhhSm6kNw9kUS4U2J9b2K1Subg08pqoh2X9yiQvbySeBY/EkKG+5rZwSKKfk1ENr0i3+8cllihWjPU/ZUp815ZQ0TTCUTbHczv3XdDX+Y/cEf28beKQTo+SSUxnqxChtuVht9MfvC/7eNvDQ8QUu3DYSvBQCTIds9dofzKQoVoxgXeItbL9CdchZhwJ6vwi/UXJBf0ZNzLZrlVJ9a2KUylWDR08scunOsfY+oLiavgG1MIsVI9jpuxbW0CgfmsmRL1Xbo2/APaeUTZFrJ6cEsP0Zqvv0GogdHjxq0jftREqVojq2xGH8Q9kUS4U2I93tV8DJB9fEKD9yYpGJ4SwT7eTTWtB3bRnlzICKFtZg7huG5JHji1zcbpRcWl41dqjn1Nqa/zueoR6P7mvvmtYA0TVKndB3elp5TGvgLT12aolSxehgYq6mbwCG+zrwFsfPUi1H1sBbfuh4B54yuC/MjN4+hbHjmerx2N3tvT8AHjy2wFC71K2VU8iu/u6G+zqgcNbQW37g6AIXtxslu9B3Q+bc7yhSPn5vz8UOB2dyLBcrXLy93bXvzhIALLcz/g1nqua0azD320VkjdJgxuLVO1iYJ+7vudjB8pTbzymspm8AhsyD8tryFoVQ418Db/GeQ3P0p/X2RB7gujAHs23WagCM7VYKxDUwyg8cW+DidqnboiVycbr/HRgly1vu8fgX8mWemFzmmWdsaO8DWkRKHdHX+WnVnLaHsByyth1Sh2a0UM+ptZ1T3fHMNaNv20FkjZJ1emlHE7O01PMGhQ8eW2C4LxW8vZCF0mpJMNgWZrEDCmfykZ53jL770ByX7xojFbSThYVy3p2+K1UxjDa83drC7K23aFG3bVF3YKtRc4qU0yy2O/ctb7nHTsk9R+YAeMbuNo2SJYluipSyKY2UJtpXoK0RfX3v4Tn60lp79Wlg1uh5REqdMAVTj9bLTUKGyBqlkX51YxZX2tyUd16pHo/c2aUrcsa+g2pTbkt5Bq6R0khfh97ijmeofMXJB9p7vw/kihUeObHYvqdsGB4URgdCFzAX5oF6NNIDPHhsgVLF4PJdbY7f2jQc7v9wNsVyuzklK1LusVG+99AcmoDLdo219wEFxTLQ10j/CSE6k0VvuRj0bM/X/l0HZ7li14bgrcUseJRDQJuREqi5L401U2AGRWSN0mi/ipQW2jVKm85VFM6h/+riVTViNlfiwKklrjpzvP0PKS2vKpyFOq/eFn0HsOs56rGH47//6DxVQ/LMPZ14yqvbzEAXFubOZ6rPPnpXe+/3gZ88PQvAs9u9/7VIyVno0rZDArD7uaosID/b/me0wN2H57hg20jtXgXGyryKFGwHHFoYyqZYavfep7LKKe3h3F8qlHnk+CLP6mTtF5ccRU51oUOHOdUeG+V2EVmjZNF3i514i2c8r6cT866DasFfddbG9j+kuOScU+g0hB+aUC1nDv6o/WtrgbsPmvRNu5FCYV499o2t+tVgpkOjtOsqVf/19A/be78P3PnUDGdPDLanPAPPnFJHQgeAPc8HJBz+cfuf4YGqIbnv8Hz71B2o+9/v/P7hvg4iJYAzrlaRQmGx/c/wwD2H5zEkPHtPB0ZpZc5x/PW13yZLMLip52u/E0TWKOmaYDibaj9SArUw5w/D/JHuXZgNdz41SzalcWk7NSqg1EEOxZNgz6l1OP7Dd/SsiHbfoTn2bh5idGC1p+sLK8qoeS/MNjem7LCShh/sjVGqGpJ9B+c6c0g8ckpDfSnypWp7tSqgvOVUX882pkeOL5IrVdunbkFFSk31iRY6KiAF2HO1orB6FC385OkZdE1wxe6x9j7AMFyNUscsAcCZL1AOyRqUBQRFZI0SwEh/msWVDr0l6Fm0dOfTM1yxe6x2/k9glFdAVh1DeEsS3jaFAcoolZZ7wi0XK1V+8vQszz27g03ZopYcFqZl6ObzHRjUM1+oVEg9yCs9emKRpWKlM+rW8uId778af9u1Wqks7Hp2z4zyD59QXdifd04H939lzjFKBtovHraw89mgpXtmlO96eo6Ld4y2T12WlpTRdJj7I30dpi5Azf1yPpTS8EgbpeG+DiOlzReqSd+DhbmwUuaRE4s8+8xOFqW5KQ+s/gxLgdTRwrSMcg/Gv+/gHCvlKi/cO9H+h1iR0sDqjX3DQAaAuXwHnt6ZL1BG/1D3Kaw7npoBOsgnAeRnFMXYtzrSHu5UfQmw5wWq3VIP8ko/fGyaC7aNsHm4r/WL3VCYd42URvrTna39zICKFntglJaLFe49MsdzzuqQugNHo5RJaQxnU8zmOpj7Z1wNCHj6B+1/Ro8QaaM02p9uP6cE6iTSs14Mj3+364V0P3hsCinhReduav9DctPq0cEoCSEY7kt1Nv6hCWWYn/he+5/hgh88NkVaF51FSh4Lc3zQNEqdLMxdV4Gegadva/8zXHDLgUn2bh5i22gbTWgt5GfUvXdQblqRckdMwZkvBGTXx58vVdh3aJYX7u1g7oOi71wipY2DGWY6ufeg1v7xeyA309nnNOFHj09Trkpecl4H5596zH2ADYOZzhyygXHYeklP5n6niLRRUvRdh/mQc38alk92ncL6/v5JNgyk25cDg9qUQCUmHTA+mGEu1+n4r1Hc8sp8Z5/ThNsem+KZZ2xon74Az4U52p9GCJjthL5L9yuxy2M3tf8ZDlgqlLnzqVleekGHhzJbRskBGwa7ECnuuFJ1TD/wnfY/wwF3PjVLuSp54bkdRMngGSmND2ZYKlTaa7Nl4dxrFEX2xM3tf4YDbj0wyXBfqsN8Wmuj1FGkBHD2S1VOzfpbIUGkjdJoN4zS3lcAAh7r3sKsVA1uOTDJS87b3F4jRgt5d/oOYONQlunlYvufD3Deq1S90hPf7exzbJhcLLD/5FLnm9LKrErGp1dHG7omGO1PM9/Jpgxw/s/AzOMw9Vhnn2PDDx+fpmJIXnb+ls4+yMMobRpSir6O7r+eUhvz4zd1tWv0rQcm6UtrnW3KlaLKebgYJcsod3T/t10OQ1vhwLfb/4wmSCm55cAkL9w70X59ErQ0SuMD6c4cElBz36jA4901yp0i0kZppC/NfKdGaXCjSvh20Sjde2Se+Xy5C56yO30HsGko07lR2vFMVa/VxYV508MnAXj5BR1uyi7qIwvjA13wFs97pXo88B+dfY4N3330FGMDaZ7RrvLKgqdRUpvy9HIXxr8y1zUVmmFIvv3QSV587ubgJ83aYUXuHvQd0BmFp2nKKD/xva6p0B4+vsipxSIvPq9Th8wHfdcpS7LjmTC0BfZ/s7PP6TIibZQmhrPkS9XOpJGgJubxe7vWC+s/Hz5JShOdRwr5GRCax8LMds6ra7qiMB+/uWsL8xsPnODcLUPtnbZpx8q8opdcMNYNb3F0p/KY93fHKJWrBrfsn+TF506031rJgodRGulLk9IEM506JWe/VOXVDnyrs88xse/QHJNLRV516bbOPmjFXXkJ9ZxiV5yS0lLXxD7fevAEuiZ4yfkdOqQ1ozTm+OvxgQ5zSqCM8nmvUka5x+3GgiDSRsk6XnpyqcOFedHPqscHv9ThFan6lK/ff5wXnzdRk262jfyM2pQ159u0cSjDfL5Mud1aFQsXvhqKC13h1k8tFrjr4CzXXrK9489S4/eIlLrhLYKiMY7e1ZV6tdsOTDGXL3PdZR2O36iqjcnFKGmaYHwww0ynkVJ2WBmmh76s/maH+I8HjpNNabys00156YR6HHY2bl2JlADOeglkR7uy9g1D8rX7jvP8czbV6NW2kZtWnTxSzp+zYTBDvlSl0G6bLQvn/4wqC3my+2KndhFpo7RlRMlNJxc7tPLjZyol1v3/1rEK7/Ynpzm1WOT6K3Z2dk2gJqbLpgQqpwQdKtBAbUqDE3D/5zv7HJSnKCVce+nWjj+LxeMw4u5xb+gGfQdw6evV4/1f6PijvnrvMcYHM11I8i+oJLzH/d80lGUm16FDBnDZG5UR6FCJVakairo7b6IzgQvAommUXO5/LVLqNFJM98FFr4VHvt5xg9J7Ds9xbH6F117RBYds8RiMuH9O1yLFs14Eg5vhvs919jldRKSN0uYRtSmf6jRSArj0DapzbocNSr967zGG+1K8rNN8EiijNOi+uW0a7FJeQU/DJT+vVFgd1KxIKfnivqNcuG2EczZ3SN0ZhtooPRbm5hEl9Ki20yncjg17VM3Off/akVOyWChz86OnuO7SbZ0luQGWVF7OTXkJKlKe6vTeA5z7SlULdV9nTsmtB6aYXCpy/RU7Or+mFpHS2EBGqS+74ZRc9iZ1dlOHuZV/v+8YfWmNn7qwWw6Z+9yfMB3SU5065HoaLnuDyqkvT3X2WV1CpI3SluEuRUoAF12vuPV7PtP2RyyslPnOQye59pJtnSV5ax94ROU8XLCxGwosC5e9EYxyRzTGPYfnefTEIm9+zu7Oryc/o06HHXHf4LaP9VMxJJNLXbj/l78Z5p6GQ7e3/RFfu+84pYrB9c/oQpS8YFKJY+7f5aahLNPdcMjSfYrCfvQbHZUG/Oudh9g8nOVlnQpcQBmlvjFH5SUo9eWW4T6OL3Th3u9+DoydAfd+tu2PWCqU+fd7j/PKi7fVWmB1hMUTnnN/xwb1vRyf78bcf4tS4T3wb51/VhcQaaM00p8im9I6zymBKia7+HUqjG1Tt/+FnxwmX6ry1uee0fn1VMsqhB/b5fqS7WPKKB+bX+n87227VKlx7vxE27mFf73jEEPZFK+9vAue8qIpOvHwFneMqYV5bK4L47/w1Sq3cNf/bevthiH5px89zaU7R7ms3V6HdswfVo+j7vd/22gfpxYL7fe/s+PKX4LKCtzzz229/chsnlsfm+KNz9rVeZQI5qbsTYPt3NDPkdl8539LCHjWO5XY4UR7TMmNdx9luVjh7c/b0/n1VCuqdtKDuraM0rH5Lox/8/mq7dJd/68recVOEWmjJIRg80iWk93wlgCe+xuqNuLu4AuzUjX459sP8pyzxrmo3eOP7Vg8rnIKHp7yttF+0rrg0EwXJibAc98Ns0+1JQ+fXCzwzQdPcP0VOzrPJ4AaP3h7i5ZR6oZRzgzCs94Bj3xNfQcBcetjkzw1neOdzz+z/bOz7Fg4oiL3Ifeo44yNA1QMyYluzP9tlykK845PtKXC/PTtB9GE4A3P7kKUDCpSGvamwXZu6OdoNxwSgGe8TR2TccfHA7/VMCT/fPtBrtg91v7ZUXYsn1Jr38Moj/SlGc6muuOQATzv3YopePQb3fm8DhBpowSwe3yAQ93wlkC13djzArjzHwJLJL/5wAmOLxR45/PP6s61WPSNh6esa4KdGwa64y0CXPBqGN0Nt/9N4NzKJ257iqoh+eUXnNmda/FhlLZ30ygBXPVrqtfc7R8L9DYpJZ+47Sm2jvTxqks6lEJbmDepWxflJcDucXXOVNeckuf9Niwdh4e/EuhtU0tF/vXOQ7z28h01R6FjzB/2pK4Bdo0PcLJbkWL/GFzxFnjwxsClId966AQHZ/L80tXdmvvH1KPH3AcVLXVt7p//MzB+FvzXX3e95VpQRN4onT0xxFOTy8hufZEvep9amAFonFLF4H/ffIALto10LoW1MN86pwBqYR7ullHSU3D1b6tCysf/0/fbJhcL/Oudh/jZK3ZwxsbVB/K1hZknID3gKfQYzKYYG0h3z1sc3qqS3vf+C8wd9P22Wx+b4idPz/LrLz67O9QVmJuyu0MCsHvjAED37v85L4ctl8AtH1IdFXzi//7wKUoVg3e/9JzuXEd+VhWObzrX82U7N/RT7VakCCpaEEKN3ycqVYOP/OdjnLtliGu75ZBMm91FNnp/nzvGuhgparpySo7f09VGAm1dyrr+9S7grE2DLBUrTHUj2Q+qSeU5L4cf/H++k76fu/MQR2ZX+O+vPB+tk7ZCdsw8AUJv6S2eMT7AoZlcd/4mwDPfDuNnw81/4rv1zEdufoyqIbu3KQFMPgoT53tGCqDu/+OTnUl5G/Di/66ipe/+ma+XG4bkf33nALvHB3hTt6grKVXro3HvqHvrSB8ZXePQbJfuv6bBT/0ZzB+Cuz7l6y2HZ/J8+vaDvPbyHZy5qUsOydQB9bjpPM+X7RpXRvnp6S6Nf2w3PPtdKq986mFfb7nx7qM8NZ3j919xXmctxeyYfFQd175hj+fLzpoY5KnpXHciRVCR4sa98J9/3LMz1vwg8kbp7M3qALSnprq4Mb/8A6pO5Ht/3vKlk4sFPnLzYzzv7I2dd0W24+SDalN2KZ6zcNbEIIuFSufSUAt6Wo1/ar9KfLbA3Ydm+cJdR3jn88/sXpQEamFuvqDlyy7YNsKjJxa7FymPbFce48NfgYOtz9n6lzsO8eiJRf7gmvPIpLq0nBaOqPm39RLPl+ma4MxNgzx2sovnQZ3zMlVQettf1mXpHvizbzxMWhP8t1ee371rmDaN0oR3pHTB1hEAHjnRxdNjX/D70DcC33qvKkvwwPRykb/8zn6uPGMDr7iwC4pDC1MHVJSoeSt4L9g2QqlidM8o62l4xV8oh/jOf+jOZ7aByBulc0yjtL+bE3PrJfDc34R9n/I8LltKyfu/9hDFisFfvPbi7iS4LZx8oOWmBNROtX3g6EL3/vYF16lGtd/7M8+kf6Fc5Q+/8hDbRvv47Zft7d7fz81AbtK3UVoqVLrHrYOiMDfsga/9hmdB5ZHZPH/1nf288NwJruu0rY4dJx9Sjz7u/yU7R3nw2EL3jDLAqz6s6Ltv/K5nfuHr9x/ne/sn+d2Xn1srZO8KJvdDql/lNz2wYTDDjrF+Hj7exbU/MK425kP/pda/Bz74zUfIFSv8z5+9pLtrf2q/UsS1wAXbemCUz70G9l4D3/9gV5sUB0HkjdLWkT62jfZx9+H57n7wS/4INpwJX/21+rlGTfjCXUe46eFTvOenzuWsidVHVreN5SmlPvKxKV24bRRdEzx4dL57f18IuO6v1cmcX3mXa37hg998hAOnlvjQz17SHcWdheP3qEcf479gmyrS7erGlBmE1/49zB2C7/w3x425VDH43X+7DwHd35ROPgAIddZVC1y6c5Tp5VL38ioAm/bCS/8YHvs23P1Pji85PJPnD7/yIM88YwNvv3pP9/42wJE7VHlCC+oW4KLtIzx8vIsOGcAVb1VdTm7+E1ca74v7jvC1+47z7pfsZW+nPR7tWDyuIuVtl7V86dkTQ2R0jUe6OfeFgFf/jaoP++q71qUnXuSNkhCCZ56xgbuenu2ut5gZgNf/E+Sm4EtvX8Wx3n1ojj/92sO88NwJfuUFXVLcWThknoa581ktX9qf0Tl3yzB3H+7ymSgj29XkPHoX/MfvrdqYv3jXEf71zsP86gvP6uwwMyc8/QMlh9757JYvvWj7KH1pjR8/2d2D2jjjeYrKufezq6gMKSUf+MbD3H1ojr963aXdU5xZePqHsOViyLZ2dC7bOQao+dhVPOc3VG71W+9bVVC8WCjzrn/Zhybgr994effEHaCOpj9xv/r+feCK3Rt4airXnQJ6C0LAaz4O2RH4wi+s6nLywNF53v/vD3H1ORv5zZec3b2/C3XKeM8LWr40k9K4fNcY//Wks9PcNoa3wqv/VjWp/sbvrLkaL/JGCeCqM8c5uVjgyW7mlQC2X6EihoM/hC//ci3x/8jxRd7x6bvYNtbH37zx8u4lOC088V3V9mXHM329/IXnbuInT892dgqtEy56LbzwvWpj/t6f1ybnzY+c4oavPsgL9m7iD67xTka3hadvUwY5M9DypX1pnavO3MgPH+9Bi5SX/BGcdy3cdIOSCpv4Pzc/xufuPMyvvehsfubSLvQ5s6O4pCKFvS/39fKLd4wyPpjhe4+e6u51aDr83Kdgwxnw+TfB8fsAWClV+dXP3M0Tk8t87Beewc4Nre9RIBz8karR8WmUXmT2GLz1sS7f/5Ft8IbPqsjlsz+ncnzAY6eWeNs//oRNQ1n+5o1XdN4JvhlP36rWvg+WANTaf+jYYufd4ptxwXVq/j/wBRUxmmv/gaPz3XX+HRALo/SKi7YiBHzzgePd//DL36Q45kf+Hb70Nu5+/Ahv/n93MJjR+ew7r2JsINPdv1cpqZNQz36pkmj7wE9dsIVyVXLL/snuXgvAi/9QKfJ+9BH4zg38+90H+fXP3s3F20f4xFue2V0vGWD6CeUp732F/0s8b4Inp3I8dqqLCX9Q9NHPfhJ2Pxe+8isYd/0T//Nbj/I333+CNz5rF+/rhUE+8B3V8uWcn/L1cl0TvPT8zXx//yTFSper8fvH4C1fURHDZ17D8v5beeun7uSOp2f4/15/WedNZ53w4I2qM/wZV/t6+QXbhtk60sdND7UWZQTGrmfBz/+LEh3986t54JFHeeMn7yCta3zuV66qtfnqGsoFVby695qWIgcLLzZZim/3YvwvfC9c+U64/W+Q33ovf/e9R3n1x/6Lbz5wovt/y4ZYGKUtI30858yNfGnf0c6OR3bD834L45q/xNj/LQY/+you7pvi8+96Tk2S2lU8+nVFGV7+Ft9vuWL3BnaN9/MvPz7U/evRNLj2/1B51rvgzr9n47+/mZft0vjsL1/V3TyShXs+raTwl73R91tec/kOMrrGZ+/owfizQ/DmL1E640Vo//G7nHH7DfzSs7fyoesv6Z783467/p/KZe5+ru+3vPbyHSwWKnztvh44ZRvOgLd/k0J2nL4v/CyXH/88f/emy3ltN5quNmN5Up1rddH1kPLn7AkheN0zd3LLgcnuFZHbcd5PI9/wWcqTj7H1336aq9MH+MK7ntNdpamFR/5dRWSX/4Lvt1y0fYQLt43w2TsOdT+CEQKu/d8UnvUbiLv+L8+87Zf4xUv6eMVFXVQaOsCXURJCXCiE+J4QIi+EOC6E+HMhRBc6jnYP73rRWRybX+Fzd3Z/Y9p/cpE3P3gFbyu9j536HJ8pvoczHv44lLuo+AIlKLjlQ6pW4OyX+n6brgneefWZ7Ds015No6QdPzPCqx36G95bfxXNTB/jEwq8y/PDnut8na/4I/OT/wcU/27LFjB3jgxlee8V2vvCTIxzsljzWRKVq8MUHZnnuoV/lH4zX8AupW/jTI+9Ee/K73efaD3xHUXdX/ZqvJL+Fq8/ZyAXbRvjb7z/OSqm792QhX+aD/5Xjqqk/4g7tGfyx/hleddc76grBbuK2v1JNeJ/zm4He9ubn7Cala3z4pgNdv6T9Jxd56w83cG3+T5HpAf628EecdccfKYVoN1HKq7W/5RI480W+3yaE4Jeu3sP+k0tdd0oMQ/L1B07wovtfxu9XfoNnpp7iz46+k+zh7hyI6IaWM18IsQH4LiCB1wB/Dvw+4K+6cI3w4nMneMHeTfzld/Zz96H2j1+wIKXkzqdmePfn7uGVf/1DHjq2wLWvfTOD79mHOO+V8P2/gI9eqtpydKPlu2HAt/4AZp+En/7LQJsSwJuu2s3ezUO878sPcLgLbWcqVYPvPHSSN37yx/ziP/6EQtngVb/4XlK/cTtiy0Xwjd+Gjz0L7v40FLqg/ikuw43vUCftvvSPA7/9919xHtmUxrs/fw8LK53n1paLFT5352Gu+egPeN+ND7B7YoSXvfvj8ItfU9TKv74OPvUKRbd048Te2afg6+9WBaNXviPQW4UQ/Ol1F6oC7q880PlRHqjWTX/57f288MO38I//9TSvuvJcLn3vt5QAYPoAfOJq+Le3qMR8i3oeX3joyypKfNYvw6ZgRdjbRvv59RedzdfvP86/dCFallLyk6dn+a3P38sr//qH3H90nre+5pVsfu9dqj/kPZ+Bj14CN/2Rops7hVFV62n+MFzzPwKv/Z99xk4u2znK+7/2UFeUePlShS/uO8K1f/sjfvvz97JxMMsv/fp/J/3rP0Jsv6Jlp4lOIVqFfEKIG4D3AWdIKRfN594HfADYaj3nhiuvvFLu27evO1fbApOLBd7wyTs4sbDCe685n1949m76M/4Dunypwt2H5rjtwBTf2z/J09M5hrMp3vycM/i1F53VmD86+F+qwPDpHyi66ZyXwdkvU4dmTZyvQl+/mDoAN/2hEji84A/gZe8PMOo6Hj+1xOs+8WOEgD++9kJec/n2QDmfqaUi9x6e45YDk9z8yCTTy0V2jPXzS1fv4a3PPYNsyvwuDUOdPfODDyv5cqofzvtpFd2d+UJ1DIDf8UupvsNvvw+mH1eKxwtf08bo4fv7T/Guz9zNrvEB/vS6C3nh3gnfFJuUkqNzK/z4qRluOzDFLQcmyZeqnL91mN99+V6uuWhrXfZdKaqN6b/+BhYOqxzIBdfBWS+GM54PwwHojUpJHRdy8/vV5vTOm1sWjbrh7255gg/fdIBnnznOn/zMhVy8w39j4FLF4NETi9z59Aw3P3KKfYfmEMA1F23l3S89p7HJcH4W7vh7pUosLqhOCBe8Wt373c9Vxad+sXRK9Vq84+NKbfmLX1NHaQREpWrwK5/Zxy0HpnjTs3fzOy/by9ZR/5+zsFLmgaPzjWu/L8UvXLWbX3/R2Y1rf+qA6vjy0I1KlLHjSlXfc+aLYMczVBGqH0ipyh/+8/2qLuplf6IUn23g6Fyen/v721kqVHjfNefxxmfv9n18jjX37zo4yy0Hprh1/yRLxQrnbB7i3S85h1dftr3rVLUQ4m4p5ZWOv/NhlH4AHJdSvtH23G7gEPBqKaVnW9m1NEqgqqx/74v384PHphjI6Dz3LEVtbB3tY6Q/TUoTlKsGxYrBbK7EyYUCJxZWeOzUMgdnckippJbP3jPOay7fzrWXbmMg45E7mTqg2pI88u/1fmmZIVX4uelcRUMNbVFnw2i6+imvqNqn2afg6D449aA6+vjlf6o8xQ5qXp6ezvE7X7iXB44uMNqf5nlnb2Tv5iG2jvYzkNHJpjQKlSq5YpXFQpmjcyscmc3z1FSuVoA6lE3xovMmuO7S7bz8gs3uCiMp1fXf/3kVMeRM6rB/XBnmTXvV2Ic2K0WR0NTYisvqtdOPw6Efq419ZKeSoJ/zsrbHDnDXwVl+9wv3cWx+hS0jWZ5z1kbOmRhi03C2lgMrlqsUKgZTiwWOLxQ4Opfn0RNLtQhLnQm0mZ+/cheX7xpzr0GqluHJW9Q5NI/dBCVTaDG8TY194151//s31I91lxIK86pbwtR+tRmtzKmN7fpPqPd1gBvvPsqfff1hlooVzpoY5MozNnDGxkE2DGToS2vomiBXrJIvqdZcR2bzHJld4cCppVo+9sJtI7zioi28/spd3nL3Ug4e/aYa/8EfQbUICGWkNl+gHgcn1E8qq36HVDnTxRNK0HLkTpBVFR2+4i9UjVibKFcN/urb+/mn2w9iSMnlu8a4ePsou8cHGO5L0Z/RKVUMChWDpUKZE/MFjs+v8OTUMgdNdiGja1x11jjXXbqd6y7b7u3ULp6AB7+oorwTD6ixaWl1DyfOU/NgcEKdHqzpav5Ximr8cwfV2GeeUGvjmv8JV7y57bGDOvDvPf92H7c/OUN/Wueqs8Y5d8sw20b7GMymSOuCQtmgUK4ymytxbH6FY3Pq3s/n1dzfNJTlJedN8Pord/GsPRu6W39nQ6dGaRL4uJTyA03P54APSCk/7PX+tTZKoCz/XQfn+Oq9x7j70CxPTuVcKY2hbIotI1n2bh7mgm0jXLprlOecuTFQhFXD3CHl9Z98ECYfURNueVItOidkR2HHFSrCuOwXYKg7aibDkNz62CTfevAkP3l6lqNzedwYnQ0DaXaND3DGxkEu2znK5bvGuGTnaD0q8gspVSNJa/xTBxQVmZtGMb8OGNoKu56tvMyLftaXBNwPShWD/3jwODc9dIoHjy24dnsQQhmg7WP9nL91hIt3jHDFrg1csG04+GKsVlTUePBHqkXS9GPq/hfmnV8vNNU1YtdVauznvDwwbeOGhZUyX73nKLccmOLh44uuh0BmdI2dG/rZNT7AuVuGuHzXBp5xxhjbRtuouyoX4OhP4PAdavxTB9SZWAWXwtb0gCoOPvMFStQTkLLzwqGZHP9+73F+8PgUj51aYqng3MNxtD/N9rF+zhgf4JKdo1y8Y5Qrz9jQnoAnP6tKR47fW7//S6fUibZOGJyAbZfDea+Ei39OKR27AIt6/OYDJ7jr4CxPTeccxV+aMBsPjPWzd/MQl+wc5bKdY1y4baQ3Ap4mdGqUysB7pZQfbXr+KPAZKeUfer1/PYxSM0oVg/l8iYWVMlUpSesaGV1jw2CmO6dEesEwYGVWLU6jquS+qayalNnhjqIivyhVDGZyRQplg2KlSjalM5jRGepLeUeB3YBRVafIFhYBqf6fGVTjb4OmaQcrpSrzKyWWzc2pL62TTWtsGMh0X9LejEpJRULW5iyEuu+DE75lv50iV6ywVKiwUq5SNSRD2RSDWZ3BTKr3G1ClqByTagnlnAgVOazR3JdSsliokCuq8Wd0jb60zmBW7/3cBxVN5mcVzWdFUoMTvtWFnaJqSBZXyiwXK1QMSV9aoy+lM9yX6n6NVQB4GSW/d8XJcgmX5xFCvAt4F8Du3V3qnNwBMimNzSN9bO5mfy6/0DQY3KR+1gmZlNae99sNaLqi74a63PUhAPozOv2ZfujC2YuBkcqoHFOQPFOXMZhN9Ua+7wepLIz2QD7uE0IIRvvTjPb7zPN0G5nBjijJTqFrgg2DGTYMro0R7Ab8mMo5YMzh+VFg3ukNUspPSimvlFJeOTHRgwK7BAkSJEgQS/gxSvuBhpa1QohdwKD5uwQJEiRIkKAr8GOUvg1cI4Swt8J9A7AC3NaTq0qQIEGCBKcl/BilTwBF4CtCiJeb+aIPAB9pVaOUIEGCBAkSBEHL7KeUck4I8TLgY8A3UHmk/4MyTAkSJEiQIEHX4EuSI6V8BPDfjC1BggQJEiRoA7HoEp4gQYIECeKBxCglSJAgQYLQIDFKCRIkSJAgNEiMUoIECRIkCA1a9r7r+A8IMYXqKN4JNgHTXbicOCL5btyRfDfeSL4fdyTfjTu68d2cIaV0bPfTc6PUDQgh9rk17zvdkXw37ki+G28k3487ku/GHb3+bhL6LkGCBAkShAaJUUqQIEGCBKFBVIzSJ9f7AkKM5LtxR/LdeCP5ftyRfDfu6Ol3E4mcUoIECRIkOD0QlUgpQYIECRKcBgi1URJCXCiE+J4QIi+EOC6E+HMhxNqcIR0iCCHOEUL8gxDifiFEVQhxq8NrhBDiD4UQR4QQK0KIHwghLl/7q107CCFeL4T4uhDimBBiWQhxtxDiTU2vOe2+FwAhxOuEELcLIWaEEAUhxAEhxB8LITK215yW300zhBA7zPkjhRBDtudPu+9HCPF283to/vk122t6+r2E1igJITYA30Uduf4a4M+B3wf+bD2va51wEfAq4DHzxwn/HXg/8FfAdcAy8F0hxNY1ucL1we+hxvke4NXALcDnhBC/ZXvN6fi9AGxEfR+/DLwS+Efgj4CP2F5zun43zfgwauzNOJ2/n5cCz7X9fMX2u95+L1LKUP4AN6COYh+xPfc+IG9/7nT4ATTbv28Ebm36fR+wAPyJ7blBYAr4i/W+/h5+L5scnvsc8PTp/L14fF//A3X0jEi+m9qYXwDMAn+AcoCHTue5A7zd/j04/L7n30toIyWUd3eTbDxI8AtAP/Ci9bmk9YGU0mjxkucBI8AXbe/Joc6/emUPL21dIaV0qiq/F9hs/vu0/F48MANY9N1p/92YqYC/RbEwzXPptP9+XNDz7yXMRul8YL/9CSnlYVSkdP66XFF4cT5QBR5vev5RTr/v6nnAI+a/T/vvRQihCyEGhBDPB34b+Hup3NvT/rsBfg3l+f+dw+9O9+/nSSFExcxF/qrt+Z5/L74O+VsnbEBRDc2YM3+XoI4NwLKUstr0/BwwIITISClL63BdawrzhOTXAO8wn0q+F8gBWfPfnwHea/77tP5uhBAbgQ8Cb5FSloUQzS85Xb+fE6h80U8AHXgT8AkhxICU8v+wBt9LmI0SKG6zGcLl+dMdbt+V2+9iBSHEHlQ+6WtSyk/bfnVafy+oyHEAeDbwJ8DHgN8wf3c6fzf/A7hTSvktj9ecdt+PlPIm4CbbU98WQmSBPxZC/LX1Moe3du17CbNRmgPGHJ4fxTmCOp0xBwwLIfQmD2YMyEspy+tzWWsDIcQ48G3gMPAW269O6+8FQEp5j/nPHwkhpoF/FkL8b07j70YIcREqmn6hEGLMfHrAfBwVQlQ5jb8fB9wI/DywhzX4XsKcU9pPE0cphNiFUnrsd3zH6Yv9qFD7nKbnV+Xl4gYhxADwTVQC/1oz6WrhtP1eXGAZqDM5vb+bvUAa+DFqk52jnlc6ihI/nM7fjxska/C9hNkofRu4RggxbHvuDcAKcNv6XFJocTuwCLzeesLcrK9DfY+xhBAiBXwJtcm8Uko52fSS0/J78cDV5uPTnN7fzY+AlzT9/JX5u1eh6pZO5++nGT+HUiceYg2+lzDTd59AqYW+IoT4K+As4APAR5pk4rGHedNfZf53BzAihHid+f9vSSnzQoi/BN4vhJhDeSy/h3I6/nbNL3jt8HHU9/I7wLgQ4jm2390rpSycpt8LQojvoIrPH0appa5GFZ//m5TySfM1p+V3Y5YS3Gp/zsxJAvxQSrlsPnfafT9CiC+jRA4PoCKiN5g/v22WpvR+Ta13sVaLQq4Lge+joqMTKLWMvt7XtQ7fwx5U6Oz0s8d8jUBV7B81v68fAles97X3+Hs5mHwvrt/NB4GHUNX28yjq7reAtO01p+V34/J9vZ2motHT8fsBPgQcQJXerAB3A29tek1Pv5ekS3iCBAkSJAgNwpxTSpAgQYIEpxkSo5QgQYIECUKDxCglSJAgQYLQIDFKCRIkSJAgNEiMUoIECRIkCA0So5QgQYIECUKDxCglSNBjCCF+Xgjx9vW+jgQJooCkTilBgh5DCHEj6pTcF6/3tSRIEHYkkVKCBAkSJAgNEqOUIEEPIYT4NKqh5YuEENL8+cD6XlWCBOFFmBuyJkgQB3wQ2I06b8Y6XO/oul1NggQhR2KUEiToIaSUTwohZgFNSnnHel9PggRhR0LfJUiQIEGC0CAxSgkSJEiQIDRIjFKCBAkSJAgNEqOUIEHvUQL61vsiEiSIAhKjlCBB77EfuEQI8VohxJVCiO3rfUEJEoQViVFKkKD3+Djwn8A/AncB71rfy0mQILxI2gwlSJAgQYLQIImUEiRIkCBBaJAYpQQJEiRIEBokRilBggQJEoQGiVFKkCBBggShQWKUEiRIkCBBaJAYpQQJEiRIEBokRilBggQJEoQGiVFKkCBBggShQWKUEiRIkCBBaPD/A+Ek/l6m//qNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams.update({'font.size': 15})\n",
    "fig,ax = plt.subplots(1,1,figsize=[7,5])\n",
    "\n",
    "ax.plot(t_eval,usol[:,0],label='u1')\n",
    "ax.plot(t_eval,usol[:,1],label='u2')\n",
    "ax.legend()\n",
    "ax.set_xlabel('t')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "# Check that numbalsoda and scipy match.\n",
    "\n",
    "# numbalsoda\n",
    "usol, success = lsoda(funcptr, u0, t_eval, data, rtol = 1e-9, atol = 1e-30)\n",
    "\n",
    "# scipy\n",
    "t_span = [min(t_eval),max(t_eval)]\n",
    "sol = solve_ivp(f_scipy, t_span, u0, t_eval = t_eval,\\\n",
    "                rtol = 1e-9, atol = 1e-30, method='LSODA')\n",
    "print(np.all(np.isclose(sol.y.T,usol)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "152 µs ± 256 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)\n",
      "11.1 ms ± 55 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "\n",
      "scipy took 73 times longer than numbalsoda\n"
     ]
    }
   ],
   "source": [
    "t_nb = %timeit -o usol, success = lsoda(funcptr, u0, t_eval, data, rtol = 1e-6, atol = 1e-30)\n",
    "t_sp = %timeit -o sol = solve_ivp(f_scipy, t_span, u0, t_eval = t_eval,\\\n",
    "                                  rtol = 1e-6, atol = 1e-30, method='LSODA')\n",
    "\n",
    "print(\"\\nscipy took \"+'%i'%(t_sp.average/t_nb.average)+\" times longer than numbalsoda\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[5.        , 0.8       ],\n",
       "       [5.02801788, 0.97855571],\n",
       "       [5.00657991, 1.19693864],\n",
       "       ...,\n",
       "       [0.1390893 , 0.1149374 ],\n",
       "       [0.1454087 , 0.1101119 ],\n",
       "       [0.15204821, 0.10552706]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# numbalsoda within jit compiled function works\n",
    "@njit\n",
    "def test():\n",
    "    usol, success = lsoda(funcptr, u0, t_eval, data)\n",
    "    return usol\n",
    "test()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypingError",
     "evalue": "Failed in nopython mode pipeline (step: nopython frontend)\n\u001b[1mUntyped global name 'solve_ivp':\u001b[0m \u001b[1m\u001b[1mCannot determine Numba type of <class 'function'>\u001b[0m\n\u001b[1m\nFile \"../../../../../../../var/folders/sf/43vm953d201c4jw4yg22hnhc0000gn/T/ipykernel_58698/4157945421.py\", line 4:\u001b[0m\n\u001b[1m<source missing, REPL/exec in use?>\u001b[0m\n\u001b[0m",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypingError\u001b[0m                               Traceback (most recent call last)",
      "Input \u001b[0;32mIn [9]\u001b[0m, in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m     sol \u001b[38;5;241m=\u001b[39m solve_ivp(f_scipy, t_span, u0, t_eval \u001b[38;5;241m=\u001b[39m t_eval, method\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mLSODA\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      5\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m sol\n\u001b[0;32m----> 6\u001b[0m \u001b[43mtest_sp\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/miniforge3/lib/python3.9/site-packages/numba/core/dispatcher.py:468\u001b[0m, in \u001b[0;36m_DispatcherBase._compile_for_args\u001b[0;34m(self, *args, **kws)\u001b[0m\n\u001b[1;32m    464\u001b[0m         msg \u001b[38;5;241m=\u001b[39m (\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mstr\u001b[39m(e)\u001b[38;5;241m.\u001b[39mrstrip()\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mThis error may have been caused \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    465\u001b[0m                \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mby the following argument(s):\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00margs_str\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m    466\u001b[0m         e\u001b[38;5;241m.\u001b[39mpatch_message(msg)\n\u001b[0;32m--> 468\u001b[0m     \u001b[43merror_rewrite\u001b[49m\u001b[43m(\u001b[49m\u001b[43me\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mtyping\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m    469\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m errors\u001b[38;5;241m.\u001b[39mUnsupportedError \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[1;32m    470\u001b[0m     \u001b[38;5;66;03m# Something unsupported is present in the user code, add help info\u001b[39;00m\n\u001b[1;32m    471\u001b[0m     error_rewrite(e, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124munsupported_error\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/miniforge3/lib/python3.9/site-packages/numba/core/dispatcher.py:409\u001b[0m, in \u001b[0;36m_DispatcherBase._compile_for_args.<locals>.error_rewrite\u001b[0;34m(e, issue_type)\u001b[0m\n\u001b[1;32m    407\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m e\n\u001b[1;32m    408\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 409\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m e\u001b[38;5;241m.\u001b[39mwith_traceback(\u001b[38;5;28;01mNone\u001b[39;00m)\n",
      "\u001b[0;31mTypingError\u001b[0m: Failed in nopython mode pipeline (step: nopython frontend)\n\u001b[1mUntyped global name 'solve_ivp':\u001b[0m \u001b[1m\u001b[1mCannot determine Numba type of <class 'function'>\u001b[0m\n\u001b[1m\nFile \"../../../../../../../var/folders/sf/43vm953d201c4jw4yg22hnhc0000gn/T/ipykernel_58698/4157945421.py\", line 4:\u001b[0m\n\u001b[1m<source missing, REPL/exec in use?>\u001b[0m\n\u001b[0m"
     ]
    }
   ],
   "source": [
    "# scipy within jit compiled function does not work\n",
    "@njit\n",
    "def test_sp():\n",
    "    sol = solve_ivp(f_scipy, t_span, u0, t_eval = t_eval, method='LSODA')\n",
    "    return sol\n",
    "test_sp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}