diff --git a/CHANGELOG.md b/CHANGELOG.md index 914b444b80..86dd9c9570 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ ## Optimizations +- Variables are now post-processed using CasADi ([#1316](https://github.com/pybamm-team/PyBaMM/pull/1316)) - Operations such as `1*x` and `0+x` now directly return `x` ([#1252](https://github.com/pybamm-team/PyBaMM/pull/1252)) ## Bug fixes diff --git a/examples/notebooks/Getting Started/Tutorial 6 - Managing simulation outputs.ipynb b/examples/notebooks/Getting Started/Tutorial 6 - Managing simulation outputs.ipynb index a343548e53..9bb1039b41 100644 --- a/examples/notebooks/Getting Started/Tutorial 6 - Managing simulation outputs.ipynb +++ b/examples/notebooks/Getting Started/Tutorial 6 - Managing simulation outputs.ipynb @@ -22,23 +22,23 @@ "metadata": {}, "outputs": [ { - "output_type": "stream", "name": "stdout", + "output_type": "stream", "text": [ - "\u001b[33mWARNING: You are using pip version 20.2.1; however, version 20.2.4 is available.\n", + "\u001b[33mWARNING: You are using pip version 20.2.4; however, version 20.3.3 is available.\n", "You should consider upgrading via the '/Users/vsulzer/Documents/Energy_storage/PyBaMM/.tox/dev/bin/python -m pip install --upgrade pip' command.\u001b[0m\n", "Note: you may need to restart the kernel to use updated packages.\n" ] }, { - "output_type": "execute_result", "data": { "text/plain": [ - "" + "" ] }, + "execution_count": 1, "metadata": {}, - "execution_count": 1 + "output_type": "execute_result" } ], "source": [ @@ -102,33 +102,33 @@ "metadata": {}, "outputs": [ { - "output_type": "execute_result", "data": { "text/plain": [ - "array([3.77047806, 3.75305163, 3.74567013, 3.74038819, 3.73581198,\n", - " 3.73153388, 3.72742394, 3.72343938, 3.71956644, 3.71580196,\n", - " 3.71214617, 3.70860034, 3.70516557, 3.70184247, 3.69863116,\n", - " 3.69553115, 3.69254136, 3.6896602 , 3.68688564, 3.68421527,\n", + "array([3.77047806, 3.75305182, 3.74567027, 3.74038822, 3.73581196,\n", + " 3.73153391, 3.72742393, 3.72343929, 3.71956623, 3.71580184,\n", + " 3.71214621, 3.7086004 , 3.70516561, 3.70184253, 3.69863121,\n", + " 3.69553118, 3.69254137, 3.68966018, 3.68688562, 3.68421526,\n", " 3.68164637, 3.67917591, 3.6768006 , 3.67451688, 3.67232094,\n", - " 3.6702087 , 3.66817572, 3.66621717, 3.66432763, 3.66250091,\n", - " 3.66072975, 3.65900537, 3.65731692, 3.65565067, 3.65398896,\n", - " 3.65230898, 3.65058136, 3.6487688 , 3.64682545, 3.64469796,\n", - " 3.64232964, 3.63966968, 3.63668791, 3.63339298, 3.62984705,\n", - " 3.62616685, 3.62250444, 3.61901236, 3.61580864, 3.61295718,\n", - " 3.61046845, 3.60831404, 3.60644483, 3.60480596, 3.60334607,\n", + " 3.67020869, 3.66817572, 3.66621717, 3.66432762, 3.6625009 ,\n", + " 3.66072974, 3.65900536, 3.65731692, 3.65565066, 3.65398895,\n", + " 3.65230898, 3.65058135, 3.6487688 , 3.64682546, 3.64469798,\n", + " 3.64232968, 3.63966973, 3.63668796, 3.63339303, 3.62984711,\n", + " 3.62616692, 3.6225045 , 3.61901241, 3.61580868, 3.6129572 ,\n", + " 3.61046847, 3.60831405, 3.60644483, 3.60480596, 3.60334607,\n", " 3.60202167, 3.60079822, 3.5996495 , 3.59855637, 3.59750531,\n", " 3.59648723, 3.59549638, 3.59452954, 3.59358541, 3.59266405,\n", " 3.59176646, 3.59089417, 3.59004885, 3.58923192, 3.58844407,\n", " 3.58768477, 3.58695179, 3.58624057, 3.58554372, 3.58485045,\n", " 3.58414611, 3.58341187, 3.58262441, 3.58175587, 3.58077378,\n", " 3.57964098, 3.57831538, 3.5767492 , 3.57488745, 3.57266504,\n", - " 3.5700019 , 3.56679523, 3.56290767, 3.5581495 , 3.55225276,\n", - " 3.54483362, 3.53533853, 3.52296795, 3.50656968, 3.48449277,\n", - " 3.45439366, 3.41299183, 3.35578872, 3.27680073, 3.16842637])" + " 3.5700019 , 3.56679523, 3.56290766, 3.5581495 , 3.55225276,\n", + " 3.54483361, 3.53533853, 3.52296795, 3.50656968, 3.48449277,\n", + " 3.45439366, 3.41299182, 3.35578871, 3.27680072, 3.16842636])" ] }, + "execution_count": 4, "metadata": {}, - "execution_count": 4 + "output_type": "execute_result" } ], "source": [ @@ -148,7 +148,6 @@ "metadata": {}, "outputs": [ { - "output_type": "execute_result", "data": { "text/plain": [ "array([ 0. , 36.36363636, 72.72727273, 109.09090909,\n", @@ -178,8 +177,9 @@ " 3490.90909091, 3527.27272727, 3563.63636364, 3600. ])" ] }, + "execution_count": 5, "metadata": {}, - "execution_count": 5 + "output_type": "execute_result" } ], "source": [ @@ -199,14 +199,14 @@ "metadata": {}, "outputs": [ { - "output_type": "execute_result", "data": { "text/plain": [ - "array([3.72947891, 3.70860034, 3.67810702, 3.65400558])" + "array([3.72947892, 3.7086004 , 3.67810702, 3.65400557])" ] }, + "execution_count": 6, "metadata": {}, - "execution_count": 6 + "output_type": "execute_result" } ], "source": [ @@ -265,16 +265,18 @@ "metadata": {}, "outputs": [ { - "output_type": "display_data", "data": { - "text/plain": "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…", "application/vnd.jupyter.widget-view+json": { + "model_id": "8ceaea70446149079f0194450d7828dc", "version_major": 2, - "version_minor": 0, - "model_id": "0b4ebac3fdd947218f9444b2b381cf04" - } + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…" + ] }, - "metadata": {} + "metadata": {}, + "output_type": "display_data" } ], "source": [ @@ -311,26 +313,28 @@ "metadata": {}, "outputs": [ { - "output_type": "display_data", "data": { - "text/plain": "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…", "application/vnd.jupyter.widget-view+json": { + "model_id": "9c9e516a7aef46688f03aaea77505636", "version_major": 2, - "version_minor": 0, - "model_id": "f4a1b65b2bf945099197135c5598084b" - } + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…" + ] }, - "metadata": {} + "metadata": {}, + "output_type": "display_data" }, { - "output_type": "execute_result", "data": { "text/plain": [ - "" + "" ] }, + "execution_count": 11, "metadata": {}, - "execution_count": 11 + "output_type": "execute_result" } ], "source": [ @@ -425,9 +429,22 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5-final" + "version": "3.8.6" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 -} \ No newline at end of file +} diff --git a/examples/notebooks/models/pouch-cell-model.ipynb b/examples/notebooks/models/pouch-cell-model.ipynb index 5e60f086eb..afb6b7c8cd 100644 --- a/examples/notebooks/models/pouch-cell-model.ipynb +++ b/examples/notebooks/models/pouch-cell-model.ipynb @@ -356,7 +356,8 @@ "comsol_solution = pybamm.Solution(solutions[\"1+1D DFN\"].t, solutions[\"1+1D DFN\"].y)\n", "comsol_model.timescale = simulations[\"1+1D DFN\"].model.timescale\n", "comsol_model.length_scales = simulations[\"1+1D DFN\"].model.length_scales\n", - "comsol_solution.model = comsol_model" + "comsol_solution.model = comsol_model\n", + "comsol_solution.inputs = {}" ] }, { diff --git a/examples/notebooks/parameter-values.ipynb b/examples/notebooks/parameter-values.ipynb index d89916f930..99644cc0e5 100644 --- a/examples/notebooks/parameter-values.ipynb +++ b/examples/notebooks/parameter-values.ipynb @@ -18,9 +18,19 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[33mWARNING: You are using pip version 20.2.4; however, version 20.3.3 is available.\n", + "You should consider upgrading via the '/Users/vsulzer/Documents/Energy_storage/PyBaMM/.tox/dev/bin/python -m pip install --upgrade pip' command.\u001b[0m\n", + "Note: you may need to restart the kernel to use updated packages.\n" + ] + } + ], "source": [ "%pip install pybamm -q # install PyBaMM if it is not installed\n", "import pybamm\n", @@ -41,7 +51,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -69,7 +79,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -146,7 +156,7 @@ "parameter values are {'a': 4,\n", " 'b': 5,\n", " 'c': 6,\n", - " 'cube function': }\n" + " 'cube function': }\n" ] } ], @@ -176,8 +186,8 @@ "parameter values are {'a': 4,\n", " 'b': 5,\n", " 'c': 6,\n", - " 'cube function': ,\n", - " 'square function': }\n" + " 'cube function': ,\n", + " 'square function': }\n" ] } ], @@ -191,7 +201,8 @@ ")\n", "f.close()\n", "parameter_values.update({\"square function\": pybamm.load_function(\"squared.py\")}, check_already_exists=False)\n", - "print(\"parameter values are {}\".format(parameter_values))" + "print(\"parameter values are {}\".format(parameter_values))\n", + "os.remove(\"squared.py\")" ] }, { @@ -341,7 +352,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3xUZfbH8c9JgYQWupSAICI9gEZQbCAWitgX2+5acO2Luru6VkSsPxsua1tU1LUjioirqIu94BoEIiIKIiXUCBJaEkJyfn/MoEmYQIBJ7iT5vl+vvGbmee7MPclcuHPmPs95zN0RERERERERqWhxQQcgIiIiIiIiNYMSUBEREREREakUSkBFRERERESkUigBFRERERERkUqhBFREREREREQqhRJQERERERERqRRKQEVqIDM7z8w+DToOERGRWGdmbmb7Bx2HSHWhBFSkijKzq81skZltMLMVZjbWzBKCjms7M2tuZi+GY8sxs8/MrG/QcYmISM1iZrXMbJKZLQ4nk/2Djqk0M7vGzOaa2UYz+8nMrgk6JpGKogRUpOp6AzjQ3RsA3YGewMhgQyqhHvAVcBDQGHgG+I+Z1Qs0KhERqYk+BX4PrAo6kDIY8EegETAIuMLMzgw2JJGKoQRUpAKY2XVm9mP4m8x5ZnZKtPfh7j+6+/rtuwSKgN0ZImRm9lD46uR8MxsY5fgWufsD7r7S3QvdfTxQC+gUzf2IiEjVURnnx9Lcfau7P+junwKFe/gyQ8Kjjn42s3vNLKqfod39Hnf/2t23ufv3wBTgsGjuQyRWKAEVqRg/AkcAKcCtwHNm1jLShmZ2tpmt38lP27J2En7uBuBnQldA/7UbMfYNx9kUuAV4zcwal7GfN3cS35vl2ZmZ9SKUgC7cjRhFRKR6qZTzYwU4BUgHDgROAi6oqJjNzAj9jb6N5i8gEivM3YOOQaTaM7PZwC3uPqWCXr8joaE7D7v7LocXmdl5wJ1Aaw//J2Bm/wP+6e7PVkB8DYDPgBfc/a5ov76IiFRNFX1+jLC/LOD37v7hbjzHgcHuPi38+DLgNHeP6sihYvu7FTgZ6OPu+RWxD5Eg6QqoSAUwsz+a2ezt33gSmqPZtKL25+4LCH1T+shuPG25l/wGagnQKqqBAWaWDEwFZij5FBGp2Sr6/Ghmbc1s0/afaL0usKzY/Qo5XwKY2RWEvlAequRTqisloCJRZmb7Ao8DVwBN3L0hMJfQPM1I259T/GQZ4ae8Q4wSgA67EWrr8DCf7doCK8qI8e2dxPd2WTsws9rA60AWcPFuxCYiItVMZZwf3X2pu9fb/hPF8NsUu7+z8+Uen9PN7ALgOmCgu2dFMXaRmKIEVCT66gIOZAOY2fmEvuGNyN2fL36yjPCzNNLzzOxCM2sevt8VuB6YXqz/QzMbvZM4mwMjzSzRzH4HdAHeKiPGwTuJb3AZ8SUCk4Bc4Fx3L9pJLCIiUv1VyvkxEjOrbWZJ4Ye1zCxp+5ewFlobe/EuXuIaM2tkZm2AK4GXoxmzmZ1DaGrMse6+qLy/l0hVpARUJMrcfR5wP/AFsBroQWj+Y7QdBnxjZpsJJY5vATcU62+zi/1+CXQkVMDoDuB0d18bxfj6AScAxwHri337e0QU9yEiIlVEJZ4fI/me0BeirYF3wvf3Dfft6nwJoaq0M4HZwH+AJ6Mc3+1AE+CrYufLx6K8D5GYoCJEItWQmaUCE929X9CxiIiIxDIzexe40t2/CzoWkZpACaiIiIiIiIhUCg3BFRERERERkUqhBFREREREREQqhRJQERERERERqRQJQQcQSdOmTb1du3ZBhyEiIlKmmTNn/uzuzSp7vzpHiohIVVDWeTImE9B27dqRkZERdBgiIiJlMrMlQexX50gREakKyjpPagiuiIiIiIiIVAoloCIiIiIiIlIplICKiIiIiIhIpYjJOaAiIjVJQUEBWVlZ5OXlBR2KRJCUlERqaiqJiYlBh1ImHUPRUxXebxGRqkwJqIhIwLKysqhfvz7t2rXDzIIOR4pxd9auXUtWVhbt27cPOpwy6RiKjqryfouIVGW7HIJrZm3M7AMzm2dm35rZlRG2MTMbZ2YLzSzTzA4s1neumS0I/5wb7V+gTJkTYWx3GN0wdJs5sdJ2LSKyO/Ly8mjSpIkShxhkZjRp0iTmryzqGIqOqvJ+i4hEVSXnTeW5AroN+Ku7f21m9YGZZvaeu88rts1goGP4py/wKNDXzBoDtwDpgIef+4a7/xLV36K0zIkwdSQU5IYe5ywLPQZIG16huxYR2RNKHGJXVXlvqkqcsU5/RxGpUQLIm3Z5BdTdV7r71+H7G4HvgNalNjsJ+LeHzAAamllL4HjgPXdfF0463wMGRfU3iGT6mN/+iNsV5IbaRUREREREJJC8abeq4JpZO6A38GWprtbAsmKPs8JtZbVHeu2LzCzDzDKys7N3J6wd5WTtXruISA22bNkyBgwYQNeuXenWrRv/+Mc/Im7n7nz44Yd8+OGHuHuFxrRy5UpOOOGEiH0333wzaWlp9OrVi+OOO44VK1ZEjBVg9OjRJR5vt3jxYl544YVfH3/zzTecd9550Qm+Bor1Y2j27Nm89dZbv/a9+eabjBo1qkL3LyJSJQSQN5U7ATWzesCrwFXuviHagbj7eHdPd/f0Zs2a7d2LpaTuXruISA2WkJDA/fffz7x585gxYwYPP/ww8+bNK7FNbm4u5513Ht9++y1z587lvPPOIzc3t4xX3HsPPPAAf/rTnyL2XXPNNWRmZjJ79mxOOOEExozZ8VvaG2+8kSlTprB27VpGjhzJnDlzSvSXTkB79OhBVlYWS5cuje4vUkPE+jFUOgEdOnQoU6dOZcuWLRW2fxGRqsBTIl4brNC8qVwJqJklEko+n3f31yJsshxoU+xxaritrPaKNXAUJCaXbEtMDrWLiEgJLVu25MADQ7Xj6tevT5cuXVi+vOR/1cnJyTz66KNMmDCBp556ikcffZTk5JL/z27evJkLLriAPn360Lt3b6ZMmQLA2LFjueCCC4DQlcbu3buzZcsWRo8ezR/+8AcOPfRQOnbsyOOPP/7ra7366qsMGhR5xkaDBg1K7DPSnL0777yTadOm8dxzz3H55ZfTq1evEv3XXXcdn3zyCb169WLs2LEADBs2jJdeeqlcfzMpKZaPoa1btzJq1ChefvllevXqxcsvv4yZ0b9/f958882K/LOIiMQ0d+fFeuezxWuV7KjgvGmXRYgsdGZ/EvjO3R8oY7M3gCvM7CVCRYhy3H2lmb0D3GlmjcLbHQdcH4W4d277hNnpYyjKySLbmrLPsDtVgEhEYt6tU79l3oroDjLp2qoBtwzrVq5tFy9ezKxZs+jbt2+J9tzcXC6//HLOP/98AC6//HIeeeSREgnEHXfcwdFHH82ECRNYv349ffr04ZhjjuHKK6+kf//+TJ48mTvuuIN//etf1KlTB4DMzExmzJjB5s2b6d27N0OHDiU/P59GjRpRu3btMuO88cYb+fe//01KSgoffPDBDv033XQTgwYNIiEhgYcffpgRI0bQs2fPX/vvvvtu7rvvvhIJSHp6OnfffTfXXnttuf5WsUrH0I7H0JgxY8jIyOChhx76dV/p6el88sknDB+uzwYiUjM99tEi/u/HLqT0uImhax4PDbtNSQ0lnxWYN5WnCu5hwB+Ab8xsdrjtBqAtgLs/BrwFDAEWAluA88N968zsNuCr8PPGuPu66IW/E2nDIW04z36+mFve+Jb39jmSjpWyYxGRqmnTpk2cdtppPPjggyWuMkLo6tWECRP46KOPgFDyUPrK47vvvssbb7zBfffdB4SWBlm6dCldunTh6aefJi0tjYsvvpjDDjvs1+ecdNJJJCcnk5yczIABA/jf//5H8+bN2dVUjDvuuIM77riDu+66i4ceeohbb721RP9tt92GmTFr1ixGjx5drvmGzZs3jzifVMqvKh1Der9FpCZ765uV/N+0+Qzr2YohZw6BHVfarDC7TEDd/VNgpzXJPXRmv7yMvgnAhD2KLgoGd2/B6Knf8p9vVnLVPvWDCkNEpFzKe5Up2goKCjjttNM455xzOPXUUyNus33YYlncnVdffZVOnTrt0LdgwQLq1au3wwf+0gmImZGcnFxiHcbzzz+fWbNm0apVqxLz+ADOOecchgwZskMCuv11txchKs/SGnl5eTsMCa2KdAzteAxFUl3ebxGR3TVzyTquenk2B7ZtyL2np1X68lO7VQW3KmreIImD2zXmrW9WBh2KiEhMcndGjBhBly5d+Mtf/rLHr3P88cfzz3/+89erjbNmzQIgJyeHkSNH8vHHH7N27VomTZr063OmTJlCXl4ea9eu5cMPP+Tggw/mgAMOYPHixb9u89RTT5UoIrNgwYISz+/cufNux1q/fn02btxYou2HH36ge/fuu/1aEvvHkN5vEZGQRdmbuPCZDFo3TOaJcw8mKTG+0mOo9gkowNAeLflh9SYWrN64641FRGqYzz77jGeffZb333+fXr160atXrx2uNJbHzTffTEFBAWlpaXTr1o2bb74ZgKuvvprLL7+cAw44gCeffJLrrruONWvWAJCWlsaAAQM45JBDuPnmm2nVqhV169alQ4cOLFy4MOJ+rrvuOrp3705aWhrvvvtumUt+7ExaWhrx8fH07Nnz1yJEH3zwAUOHDt3t15LYP4YGDBjAvHnzfi1CBHq/RaTm+XlTPuc99RVxZjx9/sE0rltr10+qAFbR63DtifT0dM/IyIja663ZkEffu6Zz5cCOXHXMAVF7XRGRaPjuu+/o0qVL0GFUutGjR1OvXj3+9re/7dA3efJkZs6cye23314pseTn53PUUUfx6aefkpCw4+yUSO+Rmc109/RKCbCYSOdIHUO7dwytXr2as88+m+nTp0d83Zr69xSR6it3ayFnPT6D+as28OKfDqF320a7ftJeKus8WSOugG4fhvufTA3DFRGpCk455RTatWtXaftbunQpd999d8TkM1rMbIKZrTGzuRH6/mpmbmZNKyyAGmZnx9DSpUu5//77KzcgEZGAFBY5V740izlZ6/nHmb0rJfncmYo708aYoT1acssb37Jg9UY6qhiRiEjgthcIKsuFF15YOYEAHTt2pGPHCq+V/jTwEPDv4o1m1obQMmVLKzqA6mZPj6GDDz64AqIREYk97s5tb87j3XmrGT2sK8d3axF0SDXjCiiEquGawZu6CioiIgFw94+BSEuRjQWuBWJvToyIiFRpT376E09/vpgLD2/PeYe1DzocoAYloM0bJNGnXWOmzllRrvXgREREKpqZnQQsd/c5QcciIiLVy1vfrOSOt75jcPcW3DAkdua115gEFOCkXq1Z9PNm5i7fEHQoIiJSw5lZHeAGYFQ5tr3IzDLMLCM7O7vigxMRkaopcyKM7Y6PbkjPSYczstksxp7Ri7i4yl3rc2dqVAI6pEcLEuON12cvDzoUERGRDkB7YI6ZLQZSga/NbIcJOu4+3t3T3T29WbNmlRymiIhUCZkTYepIyFmG4bS2n7lqy0Mkffdq0JGVUKMS0IZ1atG/U3OmzllBYZGG4YqI7Morr7xCt27diIuLI5rLYwm4+zfu3tzd27l7OyALONDdVwUcWlRdc801dO7cmbS0NE455RTWr18fdEgiItXT9DFQkFuiybblhtpjSI1KQAFO7tWaNRvzmbFobdChiIjsmfDwGkY3DN1mTqywXXXv3p3XXnuNI488ssL2UVOY2YvAF0AnM8sysxGBBVOJx9Cxxx7L3LlzyczM5IADDuCuu+6qsH2JiNRknpMVuaOs9oDUuAR0YJfm1KudwOuzNAxXRKqgYsNrwEO3U0fudQKxePFiunfv/uvj++67j9GjR9OlSxc6deq0l0ELgLuf5e4t3T3R3VPd/clS/e3c/ecKD6SSj6Hjjjvu1/VVDznkELKyYuuDkIhIdbAxr4A1cWUsJZ2SWrnB7EKNS0CTEuM5vlsLps1dRV5BYdDhiIjsngjDayiIveE1EsMCPIYmTJjA4MGDK3w/IiI1SV5BIRf9eyZ3bx1OYXxSyc7EZBi4y1p3larGJaAAJ/duxcb8bXwwf03QoYiI7J4qMrxGYlhAx9Add9xBQkIC55xzToXuR0SkJikscq5+eTZfLFrLkaddRvxJ/4SUNoCFboeNg7ThQYdZQkLQAQShX4emNK1Xm9dnL2dwj5ZBhyMiUn4pqeGhkxHa90JCQgJFRUW/Ps7Ly9ur15MYFsAx9PTTT/Pmm28yffp0zGJnKQARkarM3Rk1ZS5vz13FTUO7cErvVGB4zCWcpdXIK6Dxccawni35YH42OVsKgg5HRKT8Bo4KDacpLgrDa/bZZx/WrFnD2rVryc/P580339yr15MYVsnH0LRp07jnnnt44403qFOnzl7tQ0REfjP2vR94/sulXHJUBy48Yr+gwym3XSagZjbBzNaY2dwy+q8xs9nhn7lmVmhmjcN9i83sm3BfTNXvP7lXa7YWFvH23JVBhyIiUn5pw0PDaaI8vCYxMZFRo0bRp08fjj32WDp37gzA5MmTSU1N5YsvvmDo0KEcf/zxUfglJFCVfAxdccUVbNy4kWOPPZZevXpxySWXROGXEBGp2cZ//CPj3l/I8PRU/j6oahULNPedr4dpZkcCm4B/u3v3XWw7DLja3Y8OP14MpO9uVb/09HSv6PXm3J2BD3xE07q1mXjJoRW6LxGRnfnuu+/o0qVL0GHITkR6j8xsprunV3Yskc6ROoaiS39PEYllz3+5hBsnz2Voj5aMO6s38XGxObWhrPPkLq+AuvvHwLpy7ucs4MXdjC0QZsbpB6Xyv8XrWPzz5qDDERERERER2anXZy3nptfncnTn5ow9o1fMJp87E7U5oGZWBxgEvFqs2YF3zWymmV20i+dfZGYZZpaRnZ0drbB26tTeqcQZvPq1qkeKiIiIiEjsevfbVfz1lTn0bd+YR845kFoJVbOcTzSjHgZ85u7Fr5Ye7u4HAoOBy8PDeSNy9/Hunu7u6c2aNYtiWGVrkZLE4R2b8erMLIqKdj4UWUSkIu1qOoQEp6q8N1Ulzlinv6OIxKJPFmRzxQuz6NE6hSfOPZikxPigQ9pj0UxAz6TU8Ft3Xx6+XQNMBvpEcX9R8buDUlmRk8cXi9YGHYqI1FBJSUmsXbtWH3xjkLuzdu1akpKSdr1xgHQMRUdVeb9FpAbInAhju8PohuTf25Upzz7Ifs3q8vT5B1OvdtVeSTMq0ZtZCnAU8PtibXWBOHffGL5/HDAmGvuLpmO77kP9pAReyVjGYfs3DTocEamBUlNTycrKorKmH8juSUpKIjV179bIrGg6hqKnKrzfIlLNZU6EqSOhIBeA2puXc3vc42zt140GdWoFHNze22UCamYvAv2BpmaWBdwCJAK4+2PhzU4B3nX34tV89gEmhxecTgBecPdp0Qs9OpIS4zmxZyte/TqLDXkFNEhKDDokEalhEhMTad++fdBhSBWmY0hEpBqZPubX5HO7JPJJ+uxO6HN2QEFFzy4TUHc/qxzbPA08XaptEdBzTwOrTKcflMrzXy7lrcyVnNmnbdDhiIiIiIhITZVTRoHUstqrmKpZOinKerVpyP7N6zFpZvV4U0VEREREpGoqqNcqckdK9ZgeoASU39YEzVjyC4uyNwUdjoiIiIiI1ECLf97MmNzTyaXUXM/EZBg4KpigokwJaNgpvVsTH2dMzNBVUBERERERqVxL1m7mrMdn8B+O4JeB90NKG8BCt8PGQdrwoEOMiqpdwzeK9mmQxNGdmzNp5jL+cuwBVXZhVxERERERqVqWrdvCWeNnkFtQyIt/OoRWLRvAEX8MOqwKoSyrmKuaz+b1rZeQeHvj0Lo7mRODDklERERERKqxrF+2cOb4GWzeWsjzF/alS8sGQYdUoXQFdLvMiXSdeRMWFy55nLMstP4OVJvL3SIiIiIiEjtWrM/lrMdnsDGvgBf+dAjdWqUEHVKF0xXQ7aaPwUqtt0NBbmgdHhERERERkSjannyu31zAsyP60r119U8+QVdAf1PN19sREREREZHYsGzdFs5+IpR8PjOiDz3bNAw6pEqjK6DblbWuTjVZb0dERERERIK3ZO1mzvjXF+RsKeC5C/tyYNtGQYdUqXQFdLuBo0JzPosNw90Wn0RCNVlvR0REREREApA5MTStLyeLgnqteCL3dHL9cF68qGbM+SxNV0C3SxseWl8npQ2OsZJmjE+5UgWIREQkKsxsgpmtMbO5xdruNbP5ZpZpZpPNrOaMwRIRqQkyJ4YucuUsA5zETcu5ofBR/jNgVY1MPkEJaElpw+Hqudjo9Uw8/G3uWdGTpWu3BB2ViIhUD08Dg0q1vQd0d/c04Afg+soOSkREKtD0MSVGWAIks5VWGfcEFFDwlICW4YyD2xAfZzz/5ZKgQxERkWrA3T8G1pVqe9fdt4UfzgBUeEBEpDpRodMdKAEtQ4uUJI7vtg8vfbWM3K2FQYcjIiLV3wXA25E6zOwiM8sws4zs7OxKDktERPbU1rqtInfU4EKnSkB34txD25GTW8Abc5YHHYqIiFRjZnYjsA14PlK/u49393R3T2/WrFnlBiciInvkix/XctOmU8mjdsmOxORQAdQaSgnoTvRp35jOLerz9OdLcPegwxERkWrIzM4DTgDOcZ1sRESqhenfrebcp/7HrJRjyRs0FlLaABa6HTauRhc63eUyLGY2gdCJcY27d4/Q3x+YAvwUbnrN3ceE+wYB/wDigSfc/e4oxV0pzIxz+7Xj+te+IWPJLxzcrnHQIYmISDUSPk9eCxzl7qp6JyJSDUyZvZy/TpxD11YNePr8PjSsWwsOOSfosGJGea6APs2OVftK+8Tde4V/tief8cDDwGCgK3CWmXXdm2CDcFKvVjRISuDpzxcHHYqIiFRhZvYi8AXQycyyzGwE8BBQH3jPzGab2WOBBikiInvl2RlLuOrl2Ry0byOev7AvjevWCjqkmLPLK6Du/rGZtduD1+4DLHT3RQBm9hJwEjBvD14rMHVqJXDGwW2Y8NliVuXk0SIlKeiQRESkCnL3syI0P1npgYiISIV45MOF3DPtewZ2bs7D5xxIUmJ80CHFpGjNAT3UzOaY2dtm1i3c1hpYVmybrHBbRLFc4e8Ph7SjyJ0XtCSLiIiIiIgU4+7c/fZ87pn2PSf1asVjfzhIyedORCMB/RrY1917Av8EXt+TF4nlCn9tm9Th6E7NeeF/S8nfpiVZREREREQECoucG1+fy2Mf/cjvD2nL2OG9SIxXnded2eu/jrtvcPdN4ftvAYlm1hRYDrQptmlquK1KOrdfO37etJWpc1YGHYqIiIiIiAQsf1shV740ixe+XMpl/Ttw20ndiYuzoMOKebucA7orZtYCWO3ubmZ9CCW1a4H1QEcza08o8TwTOHtv9xeUIzo2pdM+9Xnik0WcdmBrzHRwiYiIiIjUKJkTYfoYPCeLDfHNiMs9nesHX8jFR3UIOrIqozzLsLwI9AeamlkWcAuQCODujwGnA5ea2TYgFzgzvI7ZNjO7AniH0DIsE9z92wr5LSqBmXHhEe25ZlImnyz4mSMPiK1hwiIiIiIiUoEyJ8LUkVCQiwHNCtfwQPIEEhr1BpSAlld5quBGqtpXvP8hQmXkI/W9Bby1Z6HFnhN7teKed77n8U8WKQEVEREREalJpo+BgtwSTQmFeaH2tOEBBVX1aIbsbqidEM95/drxyYKf+W7lhqDDERERERGRSuI5WZE7ymqXiJSA7qZz+rYlOTGeJz75KehQRERERESkErw/fzUrvEnkzpTUyg2milMCupsa1qnF8PRU3piznNUb8oIOR0REREREKtDEjGX86d8zeaHeuXhCcsnOxGQYOCqYwKooJaB74ILD21NY5Dz9+eKgQxERERERkQrg7jz8wUKunZRJvw5NuHTkDdiJ4yClDWCh22HjNP9zN+31Miw10b5N6nJ8txY8P2MJlw/Yn3q19WcUEREREakuCgqLuGnyXF7OWMbJvVpxz+k9qZUQF0o2lXDuFV0B3UMXH9WBDXnbeOHLJUGHIiIiIiIiUbIhr4Dzn/qKlzOW8eej92fsGb1CyadEhf6Se6hXm4Yctn8THv/kJ/IKCoMOR0RERERE9lLWL1s4/dHPmbFoLfecnsZfj+uEmQUdVrWiBHQvXD5gf7I35vPKTJVeFhERERGpyuYsW8/JD3/Oypw8/n1BH4antwk6pGpJCeheOHS/JvRu25B/ffQjBYVFQYcjIiIiIiJ74J1vV3HG+C9ISoxj8mX96Ld/06BDqrZUPWcvmBlXDNifKc8+yNb7LiMxd1VoHaCBozQ5WUREREQkVmVOhOlj8JwsNiW14K2Np9Cp1Qk88cd0mtWvHXR01ZoS0L10dMFHHFbrSZJy80MNOctg6sjQfSWhIiIiIiKxJXNi6PN6QS4G1M9byT21nsQO7Umt+ocFHV21pyG4e8mmjyGJ/JKNBbkwfUwwAYmIiIiISNmmjwl9Xi+mtudT66PbAwqoZlECurdyyihAVFa7iIiIiIgExvX5PVBKQPdWSurutYuIiIiISCA++iGblTSJ3KnP75VCCejeGjgKEpNLNHlicqhdREREREQC5+48+elPnP/U/3i2zrkUJZT8/I4+v1caJaB7K204DBsHKW1wjKyipmT2HqMCRCIiIiIiMSB/WyF/fzWT296cx3FdW3DFlTcQd2Lo8ztY6HbYOH1+ryS7rIJrZhOAE4A17t49Qv85wN8BAzYCl7r7nHDf4nBbIbDN3dOjF3oMSRsOacMpLCzijw9+TOL3cbw9yImLs6AjExGRGBHpfGpmjYGXgXbAYmC4u/8SVIwiItXNz5vyufS5mXy1+BdGDuzIVQM7hj6jhz+/S+UrzxXQp4FBO+n/CTjK3XsAtwHjS/UPcPde1Tb5LCYhPo4rB3bk+9UbeWvuyqDDERGR2PI0O55PrwOmu3tHYHr4sYiIRMHc5Tmc9NBnfLM8h4fO7s1fjj1AF4hiwC4TUHf/GFi3k/7Pi31bOwOo0bN3T0hrRcfm9fjHfxdQWORBhyMiIjGijPPpScAz4fvPACdXalAiItXU5FlZnPbo56mtF8gAACAASURBVBS588rF/TghrVXQIUlYtOeAjgDeLvbYgXfNbKaZXbSzJ5rZRWaWYWYZ2dnZUQ6r8sTHGVcdcwAL1mzizcwVQYcjIiKxbR933z5kZhWwT6SNqss5UkSkohUUFjH6jW+5+uU59G7bkKl/PpweqSlBhyXF7HIOaHmZ2QBCCejhxZoPd/flZtYceM/M5oe/Ad6Bu48nPHw3PT29Sl86HNy9BZ1b1OfB/y5gSI+WJMar1pOIiOycu7uZRTz/VadzpIhI1GVOhOlj8Jws1sc1Y13e6Vx4+B+4bnBnEvQ5POZE5R0xszTgCeAkd1+7vd3dl4dv1wCTgT7R2F+si4sz/nZcJ376eTMTM5YFHY6IiMSu1WbWEiB8uybgeEREqpbMiTB1JOQsw3CaFa3hgeQJ3NR2rpLPGLXX74qZtQVeA/7g7j8Ua69rZvW33weOA+bu7f6qioFdmnNwu0Y8+N8FbNm6LehwREQkNr0BnBu+fy4wJcBYRESqnuljoCC3RFNCYV6oXWLSLhNQM3sR+ALoZGZZZjbCzC4xs0vCm4wCmgCPmNlsM8sIt+8DfGpmc4D/Af9x92kV8DvEJDPjusGdyd6Yz4RPfwo6HBERCVik8ylwN3CsmS0Ajgk/FhGRcsgrKMRzsiJ3ltUugdvlHFB3P2sX/RcCF0ZoXwT03PPQqr6D9m3McV334bGPFnF2331pXLdW0CGJiEhAdnI+HVipgYiIVANL127hshdm8lhRE1Ljft5xg5QavTBHTNPA6Ap27aBObNm6jYfeXxh0KCIiIiIiVd60uasY+s9PWLYul18OvQ4Sk0tukJgMA0cFE5zskhLQCrZ/8/oMT2/DszMWs2zdlqDDERERERGpkrZuK+K2N+dxyXMz2a9pXd788+H0GPwnGDYOUtoAFrodNg7ShgcdrpQhasuwSNmuOuYAXp+9nPve/Z5/nNk76HBERERERKqU5etzueKFr5m1dD3n9WvH9UM6UzshPtSZNlwJZxWiK6CVoEVKEhcevh9TZq/g66W/BB2OiIiIiEiV8cH8NQwd9wkLVm/i4bMPZPSJ3X5LPqXKUQJaSS7t34Hm9Wtz69R5FBVpDXERERERkZ3ZVljEPdPmc/7TX9EyJZmpfz6coWktgw5L9pKG4FaSurUT+Pugzvz1lTm8Pns5px6oylwiIiIiIr/KnBhavzMni231W/Nw3Nk8sroXZx7chtEndiMpUVc9qwNdAa1Ep/RuTc/UFP5v2nw2528LOhwRERERkdiQORGmjoScZYCTsDGLi9Y/yMuHLuPu09KUfFYjSkArUVycMWpYV1ZvyOexj34MOhwRERERkdgwfQwU5JZoSrat9F30UEABSUVRAlrJDtq3MSf2bMX4jxeR9YuWZRERERER8ZysyB1ltUuVpQQ0ANcN7swJcZ+S/HAvGN0QxnYPDTsQEREREalB3J3nZixhhTeJvEGK6qZUN0pAA9Bq6VTuTniCJttWAx4a6z51pJJQEREREakx1m/ZyqXPfc1Nr89lSpMReEJyyQ0Sk2HgqGCCkwqjBDQI08eQWJRXsq0gNzT2XURERESkmpuxaC1D/vEJ0+ev5oYhnbnkiuuxE8dBShvAQrfDxkHa8KBDlSjTMixB0Bh3EREREamB8rcVcv+7P/D4J4vYt3EdXr20H2mpDUOdacOVcNYASkCDkJIaLjEdoV1EREREpBr6buUGrn55NvNXbeTsvm25cUgX6tZWOlLT6B0PwsBRoTmfxUpN51ttag0chQUYloiIiIhIVGRODE0vy8nCU1KZ3upiLsvcnwbJiUw4L52jO+8TdIQSEM0BDULa8NCY9vAY941JLbkmfwTT7IigIxMRERER2TuZE0MXW3KWAY7lLKPfvDFc23oO71x1hJLPGq5cCaiZTTCzNWY2t4x+M7NxZrbQzDLN7MBifeea2YLwz7nRCrzKSxsOV8+F0etJvmYeC/cZwq1T57Epf1vQkYmIiIiI7LnpY0qM9AOoY1sZkf8cTerVDigoiRXlvQL6NDBoJ/2DgY7hn4uARwHMrDFwC9AX6APcYmaN9jTY6iohPo47TunO6o153DNtftDhiIiIiIjsMS+jsKap4KZQzgTU3T8G1u1kk5OAf3vIDKChmbUEjgfec/d17v4L8B47T2RrrN5tG3F+v/b8+4slfLlobdDhiIiIiIjsFnfnjTkrWEmTyBuo4KYQvTmgrYHiZV2zwm1lte/AzC4yswwzy8jOzo5SWFXL344/gLaN6/D3VzPJ3VoYdDgiIiIiIuWyZmMelzw3k5EvzuLF+udTlJBccoPE5FAhTqnxYqYIkbuPd/d0d09v1qxZ0OEEok6tBO4+rQeL127hgfe+DzocEREREZGdcnemzF7OcWM/5oPvs7l+cGeuvOoG4k78reAmKW1CBTi1xqcQvWVYlgNtij1ODbctB/qXav8wSvuslvp1aMo5fdvy5Kc/MaRHS3q31ZRZEREREYk9azbkcePrc3lv3mp6t23Ivaf3ZP/m9UKdacOVcEpE0boC+gbwx3A13EOAHHdfCbwDHGdmjcLFh44Lt8lOXDe4My0aJHHtpEzyt2korohIdWdmV5vZt2Y218xeNLOkoGMSESmLuzN5VhbHjv2Yj3/I5sYhXZh0Sb/fkk+RnSjXFVAze5HQlcymZpZFqLJtIoC7Pwa8BQwBFgJbgPPDfevM7Dbgq/BLjXH3nRUzEqB+UiJ3nNqD85/6in/8dwHXDuocdEgiIlJBzKw1MBLo6u65ZjYROJNQBXoRkeBlTgwtrZKTxbb6rXmq9h+4I6sHB+3biHtOT6NDMyWeUn7lSkDd/axd9DtweRl9E4AJux9azTagU3OGp6fy2Ec/cnTn5qS3axx0SCIiUnESgGQzKwDqACsCjkdEJCRzIkwd+eu6ngkbs/j9hvvpcOAtHHX6EOLjLOAApaqJmSJEsqNRw7qR2qgOV0+czca8gqDDERGRCuDuy4H7gKXASkLTWN4tvo0qxYtIYKaP+TX53C7ZtnL08seUfMoeUQIaw+rVTmDsGT1Z/ksut06dF3Q4IiJSAcI1Ek4C2gOtgLpm9vvi26hSvIgEIa+gEM/JitxZVrvILigBjXEH7duYy/rvz6SZWUybuzLocEREJPqOAX5y92x3LwBeA/oFHJOI1HCf//gzgx78mOVFTSJvkJJauQFJtaEEtAq48piO9GidwvWvfcOaDXlBhyMiItG1FDjEzOqYmQEDge8CjklEaqhfNm/lmlfmcPbjX+LA5iNugMTkkhslJsPAUYHEJ1VftNYBlQqUGB/H2DN6ccI/P2HS0w9waeELWE5W6JungaO0xpKISBXm7l+a2STga2AbMAsYH2xUIlLTFBU5r36dxd1vzycnt4DL+ndg5MCOJCUOgH3q/1oFV58/ZW8pAa0i9m9ejyd6/8SBc8ZitjXUmLMsVJUM9J+AiEgV5u63EFriTESk0n23cgM3vz6XjCW/cNC+jbj95O50adngtw3ShuuzpkSNEtAq5LAlj/yWfG5XkBv6Rkr/KYiIiIjIzhRbz5OUVHKPvJF7V/TkmS8Wk5KcyD2np3H6ganEqbqtVCAloFWIqQqZiIiIiOyJUut5hkbSXcnaghGccfDZXHt8JxrWqRVsjFIjKAGtSlJSQ/9ZRGoXERERESlLpPU8yefehlOodcpdAQUlNZGq4FYlA0ftUIWsIC5JVchEREREZKfKWs+z1uYVlRyJ1HRKQKuStOEwbByktMEx1ibswzVbRzCj3sCgIxMRERGRGFRU5LySsYxVaD1PiQ1KQKuatOFw9Vxs9Hpq/e1bvml0HJc//zUr1ufu+rkiIiIiUmPMXPILJz/yGddMyuSF+udTGK/1PCV4SkCrsPpJifzrD+nkbyvikudmkldQGHRIIiIiIhKwlTm5XPXSLE579HNWb8jjgeE9ufrqG4k/KTSSDix0O2ycVlKQSqciRFXc/s3rMfaMXvzp3xncOHku9/0uDTOVzhYRERGpafIKCnn840U88uGPFLpzxYD9ubR/B+rWDn/k13qeEgOUgFYDx3bdh6uO6ciD/11AWmoK5/ZrF3RIIiIiIlKRiq3p6SmpzO44kivm7s/y9bkM7t6CG4Z0oU3jOkFHKbIDJaDVxMijOzJ3+QZue3MeB+xTn0M7lDHRXERERESqtlJrelrOMjp9dSPD6vyZI/90Gf06NA04QJGylWsOqJkNMrPvzWyhmV0XoX+smc0O//xgZuuL9RUW63sjmsHLb+LijLFn9KRd07pc8txMfszeFHRIIiIiIlIRIqzpWce28vdaLyv5lJi3ywTUzOKBh4HBQFfgLDPrWnwbd7/a3Xu5ey/gn8Brxbpzt/e5+4lRjF1KqZ+UyFPnHUxCnHHB01+xbvPWoEMSERERkShat3lrmWt6Ws7ySo5GZPeV5wpoH2Chuy9y963AS8BJO9n+LODFaAQnu69N4zqM/2M6K3PyuPjZDPK3qTKuiIiISFWXV1DIox/+yFH3fMBy15qeUnWVJwFtDSwr9jgr3LYDM9sXaA+8X6w5ycwyzGyGmZ1c1k7M7KLwdhnZ2dnlCEvKctC+jbj/dz35avEvXPfqN7h70CGJiIiIyB4oKnJe+zqLo+/7kP+bNp+++zUm4dhbQmt4Fqc1PaWKiHYRojOBSe5e/LLbvu6+3Mz2A943s2/c/cfST3T38cB4gPT0dGVMe2lYz1YsWbuZ+979gWO2fcTQNY9DTlbom7GBo1SCW0RERCSGuTvvz1/Dve98z/xVG0lLTeH+4b3ChSYPhgZJv1bB1ec7qUrKk4AuB9oUe5wabovkTODy4g3uvjx8u8jMPgR6AzskoBJ9lw/Yn8Y/TmHAD/8HFp4PmrMsVDUN9J+UiIiISNCKLaeyPZH8ou5A7n1nPl8vXU+7JnUYd1ZvTujRkri4Ymu9a01PqaLKk4B+BXQ0s/aEEs8zgbNLb2RmnYFGwBfF2hoBW9w938yaAocB90QjcNk1M+OsTU9hVqoYUUFu6D86/aclIiIiEpxSy6mQs4z8yVfwYv4IVtQ7hrtO7cHpB6WSGF+uhStEqoRdJqDuvs3MrgDeAeKBCe7+rZmNATLcffvSKmcCL3nJCYddgH+ZWRGh+aZ3u/u86P4KsjNlVkMro3qaiIiIiFSSCMup1PZ87mgwmcS/3k5SYnxAgYlUnHLNAXX3t4C3SrWNKvV4dITnfQ702Iv4ZG+lpIaG3UZqFxEREZHAeE4WFqG9ft4qUPIp1ZSu51d3A0ftUCUtl1qsTL82oIBEREREarasX7Zw/WvfsELLqUgNpAS0uksbDsPGQUobwCio35o74y5l2Eet+DF7U9DRiYiIiNQY2xPPAfd9yKszs/hs38soStByKlKzRHsZFolFxaqkJQLnrtnI2+NncPbjM3j5okNp17RusPGJiIiIVGPL1m3hkQ8X8kpGFnFmnHlwWy7t34FWDQdDZhstpyI1ihLQGmj/5vV5/sJDOHP8F6Ek9OJDadO4TtBhiYjUWGbWEHgC6A44cIG7f7HzZ4lIzCm1pMraQ67j3hVpTJoZSjzP7htKPFumFLvqqeVUpIbRENwaqlOL+jx3YV82by3krMdnsHx97q6fJCIiFeUfwDR37wz0BL4LOB4R2V3bl1TJWQY45CwjedrVbJ31Muf0bctH1/ZnzEndSyafIjWQEtAarFurFJ4b0Zec3AKGP/YFS9ZuDjokEZEax8xSgCOBJwHcfau7rw82KhHZbRGWVKljW7mn0evcqsRT5FdKQGu4HqkpvHDhIWzZuo3fPfYFC1ZvDDokEZGapj2QDTxlZrPM7AkzKzE538wuMrMMM8vIzs4OJkoRKdOspb/gZayxnrCxjDXZRWooJaBCj9QUXr74UBw4Y/wM5i7PCTokEZGaJAE4EHjU3XsDm4Hrim/g7uPdPd3d05s1axZEjCJSirvz6YKfOfvxGZzyyOesREuqiJSHElAB4IB96vPKxYeSnBjPWY/PYOaSdUGHJCJSU2QBWe7+ZfjxJEIJqYjEoKIiZ9rcVZz88Gf8/skvWbhmEzcO6ULjE+/YYe11LakisiNVwZVftWtal4mXHMrvn/iSl568n271XiNpy0qVBBcRqUDuvsrMlplZJ3f/HhgIzAs6LpEar1RF220DbmZK4WE8+tGPLFyziX2b1OHOU3pw2kGtqZ0QD+wHCXFaUkVkF5SASgmtGybz+pHLqf3W4yRtyQ815iwLVXUD/ScqIlIx/gw8b2a1gEXA+QHHI1Kzba9ou72oUM4yCl7/Mx9tHUFC88GMO6s3Q7q3ICG+1GBCLakisktKQGUHKZ/dBeSXbCzIDX2jp/9URUSizt1nA+lBxyEiYREq2iaTzz0NX6f2lXdiZgEFJlL1aQ6o7KiMKm5lVXcTERERqS6+3klF26QtK5V8iuwlXQGVHaWkhhdRLmltfDOS8rdRr7YOGxEREak+Couc9+at4vFPfmLmkl/4PKkJrfh5xw1V0VZkr+kKqOxo4Kgdqrhti0vi9vzfcdojn7Ns3ZaAAhMRERGJng15BTz56U8cff+HXPLc16zekMctw7qqoq1IBdKlLNnR9nmexaq4JQwcxenJR3PZ8zM58aFPeeScgzi0QxnrXYmIiIjEilLVbBk4ioUtBvPM50t49esstmwt5MC2Dbn2+M4c322fcGGh9qpoK1JBzN13vZHZIOAfQDzwhLvfXar/POBeYHm46SF3fyLcdy5wU7j9dnd/Zlf7S09P94yMjPL+DlKJfvp5Mxc+8xVL1m7h1pO6cU7ffYMOSUQkEGY2090rvXCQzpEiu6F0NVsg32pzTf4IptmRDOvZivP6taNHakqAQYpUT2WdJ3d5BdTM4oGHgWMJLZb9lZm94e6l1yh72d2vKPXcxsAthCr7OTAz/Nxf9vD3kIC1b1qXyZcfxsgXZ3Hj5Ll8u2IDo07oSlJifNChiYiIiJQUoZptbc/n9vqvMerPo2lar3ZAgYnUXOWZA9oHWOjui9x9K/AScFI5X/944D13XxdOOt8DBu1ZqBIrGiQl8uS5B3PxUfvxwpdL+d1jX2heqIiIiMSU71dtLLOabYP81Uo+RQJSngS0NVC8JGpWuK2008ws08wmmVmb3XwuZnaRmWWYWUZ2dnY5wpIgxccZ1w/uwvg/HMTitZsZOu4T/jtvdWioy9juMLph6DZzYtChioiISA2RV1DIqzOzOO3Rzzn+wY9Z4WXUq1A1W5HARKsI0VTgRXfPN7OLgWeAo3fnBdx9PDAeQvNbohSXVLDjurXgPy0acOnzM3njuQc5KmkCiUV5oc6cZaF5F6BJ+yIiIlJhFq7ZxAtfLuXVr7PIyS1gv6Z1uWloF1KSboN3/1JyGK6q2YoEqjwJ6HKgTbHHqfxWbAgAd19b7OETwD3Fntu/1HM/3N0gJba1bVKHVy/tR+49I0gsyCvZWZAbmn+hBFRERET2VrGKtp7Smpn7j+TeFWl8+dM6EuON47u14Oy+bTl0vyaYGbAf1E5QNVuRGFKeBPQroKOZtSeUUJ4JnF18AzNr6e4rww9PBL4L338HuNPMGoUfHwdcv9dRS8xJSownqWBN5M4y5l+IiIiIlFupiraWk0XXjJvoWusy+g/6I79LT408rzNtuBJOkRiyywTU3beZ2RWEksl4YIK7f2tmY4AMd38DGGlmJwLbgHXAeeHnrjOz2wglsQBj3H1dBfweEgtSUkPDbkspatC6XJONRURERCLJ3VqIvz2KOqUq2taxrYyq8yrW/9aAIhOR3VWuOaDu/hbwVqm2UcXuX08ZVzbdfQIwYS9ilKpi4Kgd1tra4rX4v82nc/yPP9OvQ9MAgxMREZGqxN2ZvWw9EzOyeHPOCuawEmzH7UwjrUSqlGgVIRL5bXhLsXkW2b3+xscz2/HM419yTt+2XDe4M/WTEoONU0RERGJW9sZ8Xp+1nIkZy1iwZhNJiXEM6dGSrUtakbR5xY5PUEVbkSpFCahEV6l5FvsCbx1WyAPvfc+Tn/7EB/PXcOepPejfqXlwMYqIiEiwihUTIiWVbQNu5oNa/ZmYsYwP5q9hW5FzYNuG3HVqD05Iaxn68jrz1h1GWqmirUjVowRUKlxyrXhuHNqVIT1acu2kTM576itOPbA1o07oSsM6tYIOT0RERCpTqWJC5Cyj4PU/M3XrCGbVGciIw9vzu/RU9m9ev+TzIoy0UkVbkarH3GNvyc309HTPyMgIOgypAPnbCnno/YU8+uGPpCQnct3gzpyW+AVx7+tkIiJVi5nNdPf0yt6vzpFS1W27vxsJG3ect5lbpxUJf/2WxHiVLhSpDso6T+pfuFSq2gnx/PW4TrxxxeHs26QOn7z2CFsnXxGunuuh26kjQ9+OioiISLWQk1vAy18t5Yx/fUHchshFg5K3rFTyKVIDaAiuBKJrqwZMuqQfufeOICk3v2RnQW5oeI2ugoqIiMS+UvM5t49k2rqtiI9+yGbyrCz++90atm4rYr+mddmU1IIG+at2fB0VExKpEZSASmDi4oy6uRFOQIDnZEWqtC4iIiKxJMJ8zsIpI3n1q2XctbwHv2wpoEndWpzdpy2n9G5NWmoK9s1tKiYkUoMpAZVgpaSGh9+WtMaa8t33a1QtV0REJJZNH1MykQTiC3M5bOkjHN55Eqf2bs3hHZuWHFqrYkIiNZoSUAnWwFE7fAu6LT6Jf8X9nglPfcURHZtyw5AudGnZIMAgRUREpDh3Z05WDj3LGLHUytbyz7N6l/0CpZZtE5GaQwmoBCvCt6AJA0dxXdfTaT1jCeOmL2DIuE84pXdrrhzYkX2b1A02XhGRCmJm8UAGsNzdTwg6HpHS3J1vV2zgzcyV/OebFSxbl8tntZvQ2n7eYVvTfE4RKYMSUAlehG9BawEjDm/P6Qem8vCHC3nm88VMmb2C3x2UyhVH709qozplFj0QEamirgS+AzTkQ4JT6tzqA0fxTePjmDZ3FW/PXcVPP28mIc44bP+mjDy6I424Hd65WvM5RaTclIBKTEupk8gNQ7pw4eHteeTDH3nhy6W8+nUWd3T4jtNX3Evctt+KHjB1ZOi+klARqWLMLBUYCtwB/CXgcKSmilBQKO+1K3hi6wj+wxEcul8TLjpyPwZ1a0GjurXCTzoLasXrC2ERKTdz96Bj2IEW2ZayrMzJ5eEPFnLprJMjDvkhpQ1cPbfyAxORGqesBbb38LUmAXcB9YG/7WwIrs6RUhHytxVS9EB3kres2KFvc3JLCv6cScM6tSI8U0QksrLOk7oCKlVKy5Rkbj+5Bz57beQNciIvbi0iEqvM7ARgjbvPNLP+ZWxzEXARQNu2bSsxOqnONudv46Mfspk2dxUfzF/DHFYQqaJQ3dxVoORTRKJECahUSVbG8i3Z8c347odsjujYFDOtJCoiVcJhwIlmNgRIAhqY2XPu/vvtG7j7eGA8hK6ABhOmVCll1ElYlZPH+/PXMP271Xy68GfytxXRuG4thvRoSf6iVhGvgKKCQiISRUpApWqKsHxLQVwS4/wsnp3wPzo0q8u5/dpx6oGp1Kutw1xEYpe7Xw9cDxC+Avq34smnyG6LMJez4PU/8+A783l47UEApDZK5qw+bTm+WwsObteIhPg4yLx1h3OrCgqJSLSV65O5mQ0C/gHEA0+4+92l+v8CXAhsA7KBC9x9SbivEPgmvOlSdz8xSrFLTRZh+ZbEgaO4qetp9M5cyTOfL2bUlG+5d9r3nHZQKuf2a0f7pnVVOVdERKq9ov/eSlzxJBJILMrj3NxnqTvobAZ23ocD9qm340ihCOdWnSdFJNp2WYQovC7ZD8CxQBbwFXCWu88rts0A4Et332JmlwL93f2McN8md6+3O0GpwIJEw6ylv/DM54v5zzcrKSh0/t4qkz/lPEhCYd5vGyUmw7BxOrmKyG6LZhGi3aFzpESyesNvQ2vHLzqWOIv0+c5g9PpKj01Eaqa9KULUB1jo7ovCL/QScBLwawLq7h8U234GoKFDErjebRvRu20jbhjahRe/XMbJn44kgbySGxXkhr7pVQIqIiKxKsLona1dTydjyTo++iGbj77PZv6qjUBoaO2GWs1pWLB6x9fRXE4RiQHlSUBbA8WrvWQBfXey/Qjg7WKPk8wsg9Dw3Lvd/fVIT1KFP6kozesnceUxHfFPI1fO9ZwsNuUVUD8psZIjExER2YUI8znzJ1/BTZPm8MrWfiTGG+n7Nua6wZ0Z0Kl5aGjtN7drLqeIxKyoVmcxs98D6cBRxZr3dfflZrYf8L6ZfePuP5Z+rir8SUUrq3Lu8qImHH37fxnYuTkn9WpF/07NSUqMDyBCERGR3+RuLYRpt5Bcaj5nbc/nhtqvcNyZIzm0Q5Mdi+1pLqeIxLDyJKDLgTbFHqeG20ows2OAG4Gj3D1/e7u7Lw/fLjKzD4HewA4JqEiFi1A51xOT2dbvJs7e2JY3M1fw9txV1K+dwKDuLTixVysO3a8JCd9O0klcRESiYyfF8IqKnPmrNvL5/7d359Fx1XUfx9/fmclkabO16ZakC4WylrIKCJyHKiCr4kYpLg8uj4gLKp7zPIoVKOXBx6MeWY6eR3iAoyiKiigVgbKKIhZaNlso0NIWSBpamjZJ2zTbzO/5Y27aSXInmaTJzb3p53XOnEzuvb+Z7+/emfnOd+69v/vGVv62divPrG9kTdz/2pyVnVs48/ApuZ9n3gLlKhEJpXwK0BXAHDM7gEzhuRD4RPYCZnYMcAtwtnNuS9b0SqDVOdduZlVkrnX2g+EKXmRQfH4RttOvZta8BSwGvnveYfxzfSP3vbiJh1a/w++fq+OTxcu5hltIdv+m0vx2pojNfjwREZF8+BxOm176NZa/0civd5/I0280sm1XBwAHThrHp06aSccr1RTp2pwiMoYMOAougHdx7BvJXIblDufc9Wa2BFjpnFtqZo8CRwINXpO3nHMfMrOTyRSmaSAGFHYqXAAAFmhJREFU3Oicu32g59MIfzLa2jpTPPHqFk647zQmdvUdyCFdVkvsmy+PQmQiEhYaBVcG7Ya5vqeC1KWr+FjRLZxyYBWnHFTFyQdNZFp5cWZm76IVNIK7iETCvoyCi3PuAeCBXtOuzrp/Ro52T5MpTEUipaggzjlHToM/bPFfoLmef7/jWeYfPIn5h0zigKpxfa+nJiIi+72d7V2s2LCNp9ZtZVFzHTGfZWpijSy/8nT/PKLzOUVkjBnWQYhExpwcAxe1FE6hblsrS+5/hSX3w/QJxZx28CROO3gyJx84kXGFiX7P8xERkYjL8Rnf0tbJyo3beGb9NpZv2Mbq+mZSaUcyEeOywklMSvX9YdPKa6G/HzF1PqeIjCEqQEX64zNwEQXFVJx/HY/Pm89bja08uTZzDbZ7n6/nV8vfoiBufH3Si3yx5SYK0t51R3XuqIjI2OFzLmfHHy/nxmWv8rNtx5F2UBA3jqqt4LLTZvPe2VUcP6uSojXX6/IoIrLfUwEq0p8BDn2aMbGET0+cyadPmkl7V4rnNm7nydff5cIVX6PAtfV8rM7dtC9bTOyIj1MQ9zsIS0REwm7rznZKHryGkl6XRkm6Nj7b9ksS71/ISQdM4JgZlRQne13SS4fTiojkNwhR0DTAgkTe4gqg73sr7Yy57m6Om1nJsTMqOXZmJUdPr6C8uCCzgA7bFYkMDUI0huT47HXO8WZjKyvf3M5zb25nxcZtrNuyk/WFnyDme8SsweKmoKMXEQmlfRqESEQGKce5o+0l07jw0Fqe2bCNmx9fi3OZ034OmjSez5Wt4MJNPyShw3ZFRILjczht558u5xd/X8/Pth/H1p2Zy6KUFiU4bmYlHzu2ls4VNRTu6nNJdF0aRUQkDypARUZCjnNHi8+5lmvnzQVgR1sn/6pr5vk3t/P8W9s57c3/JUHfw3ZbH7yaTVPP5YCqccT9fnLXXlMRkSHZurOdcQ9dQ3Gvw2kL0m2cv/U2XjnsbI6fOYHjZlYyZ/J4Yt2fwRMW61xOEZEhUgEqMhLyOM+ntKiAUw7KXPMNwC1u9H2ootZ3OOPHT1KSjHP4tDLm1pRzRHUZh1eXMWfzgyQf+EaPX+6111RE9lv9/CDXmUrz2js7ePHtJl54q4nn3tzGxsZW1hduAp/f9qa6rfx4wdH+z6NzOUVEhkwFqMhIGeSw+ZbjsN1UaTU/+tBRrK5v5uVNzfxu5du0dqQAeKpwEbXW85d7OnfjHrsWG+i5tedURMYSn0Npu+67nKUv1HNX64msrm+mvSsNwMRxSY6dWcnFJ8yg89khHk6rS6OIiAyJClCRsMhx2G7BBxbz8Xm1fPy4zJehdNqxoXEXrzbsoOZe/72mrqmes254kjlTSpldNY4Dsm4VJUnfL2racyoiUdXU2kHRQ9dQ1OtQ2kSqjRPX/4RfTz2RT500k6OnV3D09ApqK4ux7utuVi7W4bQiIgFSASoSFnke0hWLGQdOGs+Bk8bDY/57TXcUTmHGhBJW1zfz4KoG0lkD8laWFLCMRUxO++05XdL/nlPtNRWRoOT4vNnZ3sWahhZW1TXzUl0TL73d1O+htNXWyD1fOjn38+hwWhGRQKkAFQmTwR7SlWOvafn513HbvPcA0NGV5u3trWx4dxcbtu5iQ+Muql561/fhXFMd59z4N2orS6itLPZumfuzGx6gZNkVQ9trqsJVRAbD5yiNjj9+lR/8ZQ23txxP9xXkppYVcdT0ci56zww6VlRTtGtTn4eyfEam1eG0IiKBUQEqEmV5/HKfTMT27jHtttF/z2lLcjK1lcXUbW/ln29sZZd3rinAU8nvUhLru9d05wNXszxxGlPLi5hSVsTEccm9I0XC0A/3VdEqMnbk8X7evquDlze1sKq+mQVPfZeJXT0/b5KunctSv6L09E9wZG0Zc6vLmVxWtHeBCdfqUFoRkQgw59zASwVMF9kWGWG9i0LIfFH74M17vhQ652je3Und9t3UbW/lrHsOw+j7eZF2xuz2u/Y+TNyYXFrkFaSFfG/jxVR0bu7brqwWu2L13vOwBhlfzn6paJWA5LrA9kiLXI70eT+nE8W8cvx1/DU5n9X1Laze1Ezd9r3z1xd9kpjP5w0YLG7q/7n0GSAiEgq58qT2gIrsj/LYc2pmVJQkqShJMremHB7x32uaLqvhTwtO4Z3mNja3tPFOSxubmzN/X31nB2WdW/xjaK7n0Kseomp8IVXjk0wcX8jEcUkmjE/y1RevprSz797W9KPXEsv1ZXJfBlbSl1aREdHWmSL28GKSvd7Psa7dVDz9fX7UUc2siSUcNb2CT500k7nV5cytKSN2i//njUamFRGJPhWgIvurYTrfNHHmNRw9vQKm52h3g/8XyV1FU7nkqFls3dnO1p0dbG5p45VNLTTuaudbiQbfwURorueQ7z5IZUmSipICyosLqCgpoLIkyZWvXUW5T9Ha8fBi6qvPY3xhgtKiBIWJWM+9rkEfIqxiV6Isx+vXOUdDcxuvvtPCmoYdrGloYU1DCxu27mJdst73/VwTa2TV4g9QWlTQd2aOzxsdTisiEn0qQEUkP0MdKTLHF8nS85bwnXmH9VncOYe7oRZa6vrM21k0hUuOmkVTawdNrZ007e5k49ZWXmht4nsdm32/5CZ2bOJ9P/rr3qeOG6VFBXsK0p83LWKSz4jALX+5int3vIeSZILiZJxxhXGKCxKUJONM3riUKU/+F7GuIRStQRW7QRbIKqr3DzkGBrrl8bXc3vIemlo79yxaW1nMYdPKOO/IabS/UE3xbv/BgXyLT9DItCIiY1he54Ca2dnATUAcuM059/1e8wuBO4HjgEbgIufcRm/elcDngRTwNefcsoGeL3Lnt4hI/wZboAzlHNAb5vruad1dUs1DZz7CjrauPbed7Z2Zv21d/N/GM33PNet9bmu2p5Jfoza2tc/0Bqq4qOQ2ChMxCgtiFCXiFBbEKEzEKSqI8d8bLmZCV9/zYXcUTeO++cu8dnGKvL+FiRhT31zKjH9cubfYBVyimF1n/RjmLaAgbhTEYv0P/JTP+htqu6E+V3fbIIrdESqQh+scUDObTiaHTgEccKtz7qZcyw9bjhxgvTTv7mTt5h28vnknr2/ewVde+jCTUn0Pqd8Sm8wNc+/lsGmlHDatjEOmllKWXVjuy2tEREQiK1eeHLAANbM48DpwJlAHrAAuds69krXMl4F5zrnLzGwh8BHn3EVmdjjwG+AEoBp4FDjYOZfq/TzZVICKSCBFK+QsXNNltTR98QV2tXexuzNFa0eK1o4uWttTnP67Q3wHZHIYVxz6BO1dadq70rR1prz7Kdo60zzccsGwFbt16SpO7bh5z/+JmJGIGwXxGMv4MtX0bbMlNpnLp9xJMhGjIB6jIG4k4jESMSMeM65au4BKnwGjmpNTuf34pcRjMRLxzLLdbS78+zmMb2vo02Z3cTV/P/8JEnEjZkYiFsu089pPfOM+pj/1bWKpngPTvPu+H7D7kI8RjxlmEI9l2sfMKFrzB8Y//E2sVzGeOv8mbN4CYkbfQa1GsPgZxgJ0GjDNOfe8mZUCzwEfzs6z2YYlR/qsl65YEffN+BZ/7DqFtVt2sLmlfc+8kmSc1bGFQxsYqPv5tDdTRGS/si+DEJ0ArHPOrfce6G7gAiA7MV4ALPbu3wP8xDLfAi4A7nbOtQMbzGyd93j/HGpHRGQ/MdhzVIf5EOHYGdcwYVySCeOSfduU+5/XauW13LjwmNzPleN8WMprePYLp9PeubdYbe9K0d6ZpuZXjb4PVRNrZNG5h9GRStOVcnSm0nSm0nSk0kx7zr/NpPS7mEFrR4rOVCcdXZk2aQdd6TTlOQaMKm3fzM2Pr/Odd0mh//m6ha0NXPrL5/zXA/BU8jpisb4D03Quu5b5f56Yo81VlPZqY127eefe73Dq3ZnLDJlB3CtYYzF4LHYlNdb3EGseWxKaAsg51wA0ePd3mNkaoIaeeXZ4Pbak52seSKTbOGnDT/l51QmcclAVB08p5eAp45kzuZSaimJiNw1xYCDQ4EAiIrJHPgVoDZCdceqAE3Mt45zrMrNmYKI3fXmvtjV+T2JmlwKXAsyYMSOf2EVEehrKl9yhFK5DHSCln2J3cmmRf5t+it0v/Nts/zbrcre5+9L35o4vR4Ecq6hl4xXnkU47utKOtMv8TaUc7mf+5+umSqu5/9JTSWW3STlSaUdnOk3Nb3IX1j9ecBRpB+m0I+UybdNpR82y3G2uOOPgzHLeLZXOnE9c/ax/G5r7xhwGZjYLOAZ4ptf04c2ROfpfbY38+fJT/dtoYCARERkGoRmEyDl3K3ArZA4vGuVwRGR/EtTe1qCK3WEukLvbxWJGMtZrd+cZ1/i2KfjA4szle3Lpp7D+6LE59qgtz93m62fM8W/z2j7stQuYmY0H/gB8wznXkj1v2HNkP+s/Jw0MJCIiwyCfArSenhdYqPWm+S1TZ2YJoJzMYET5tBURiZ6hHlIYRLEbZIE8zIc+D3thHZG9dmZWQKb4vMs5d++IP+FQ14sOpRURkX2UzyBECTKDEJ1OpnhcAXzCOfdy1jJfAY7MGoToo865BWZ2BPBr9g5C9BgwR4MQiYiIRsHd8zgG/ALY5pz7xkDLBzUKroiIyL4Y8ii4XuNzgRvJXIblDufc9Wa2BFjpnFtqZkXAL8mct7INWJg1aNEi4HNAF5nDih4c6PlUgIqISNgNYwF6KvB3YBWQ9iZ/xzn3gN/yypEiIhIF+zIKLl4SfKDXtKuz7rcBF+Zoez1w/aCiFRER2U84557CdyxhERGRsSc22gGIiIiIiIjI/kEFqIiIiIiIiARCBaiIiIiIiIgEQgWoiIiIiIiIBCKvUXCDZmbvAm8O08NVAVuH6bFGS9T7EPX4Ifp9iHr8EP0+RD1+iH4fhjv+mc65ScP4eHkZ5hwJ2q5hEPU+RD1+UB/CIOrxQ/T7EEieDGUBOpzMbOVwDJM/mqLeh6jHD9HvQ9Tjh+j3IerxQ/T7EPX4R0rU10vU44fo9yHq8YP6EAZRjx+i34eg4tchuCIiIiIiIhIIFaAiIiIiIiISiP2hAL11tAMYBlHvQ9Tjh+j3IerxQ/T7EPX4Ifp9iHr8IyXq6yXq8UP0+xD1+EF9CIOoxw/R70Mg8Y/5c0BFREREREQkHPaHPaAiIiIiIiISAipARUREREREJBCRLkDN7Gwze83M1pnZt33mF5rZb735z5jZrKx5V3rTXzOzs4KMOyuGgeL/ppm9Ymb/MrPHzGxm1ryUmb3o3ZYGG3mPGAfqw2fM7N2sWP8ja94lZrbWu10SbOR7Yhgo/huyYn/dzJqy5o36NjCzO8xsi5mtzjHfzOxmr3//MrNjs+aN+vr34hioD5/0Yl9lZk+b2VFZ8zZ60180s5XBRd0jvoHin29mzVmvlauz5vX7+gtKHn34z6z4V3uv/QnevDBsg+lm9oT3efmymX3dZ5nQvxeGW9RzpBdHpPNk1HOkF4fy5Oh+T4l0jvTiUJ5UnuzJORfJGxAH3gBmA0ngJeDwXst8GfiZd38h8Fvv/uHe8oXAAd7jxEMY//uAEu/+l7rj9/7fGZFt8BngJz5tJwDrvb+V3v3KsMXfa/nLgTtCtg3+DTgWWJ1j/rnAg4ABJwHPhGX9D6IPJ3fHBpzT3Qfv/41AVci3wXzg/n19/Y1mH3ot+0Hg8ZBtg2nAsd79UuB1n8+i0L8XhnmdRDpHDqIPoc2Tecb/GUKaI/PtQ6/llSeDjz/UOTLPPsxHeXKk4w9VnozyHtATgHXOufXOuQ7gbuCCXstcAPzCu38PcLqZmTf9budcu3NuA7DOe7wgDRi/c+4J51yr9+9yoDbgGAeSzzbI5SzgEefcNufcduAR4OwRijOXwcZ/MfCbQCLLk3Pub8C2fha5ALjTZSwHKsxsGuFY/8DAfXDOPe3FCCF8H+SxDXLZl/fPsBpkH8L4Pmhwzj3v3d8BrAFqei0W+vfCMIt6joTo58mo50hQnhz1bRD1HAnKk2EQtjwZ5QK0Bng76/86+q7IPcs457qAZmBinm1H2mBj+DyZXyW6FZnZSjNbbmYfHokA85BvHz7m7cq/x8ymD7LtSMo7Bu+wrgOAx7Mmh2EbDCRXH8Ow/oei9/vAAQ+b2XNmdukoxZSP95rZS2b2oJkd4U2L3DYwsxIySecPWZNDtQ0scxjpMcAzvWaNtffCQKKeIxlCHGHLk1HPkYOKQ3kyFKKaI0F5MjBhyJOJfWkswTCzTwHHA6dlTZ7pnKs3s9nA42a2yjn3xuhE2K8/A79xzrWb2RfJ/Nr+/lGOaSgWAvc451JZ06KyDcYEM3sfmeR6atbkU71tMBl4xMxe9X6lDJPnybxWdprZucCfgDmjHNNQfRD4h3Mu+1fg0GwDMxtPJul/wznXMhoxyOiIcJ4cKzkSlCdHVYRzJChPBiYseTLKe0DrgelZ/9d603yXMbMEUA405tl2pOUVg5mdASwCPuSca++e7pyr9/6uB/5K5peMoA3YB+dcY1bctwHH5ds2AIOJYSG9DqcIyTYYSK4+hmH9583M5pF5/VzgnGvsnp61DbYAf2R0DhPsl3OuxTm307v/AFBgZlVEbBt4+nsfjOo2MLMCMkn1LufcvT6LjIn3wiBEPUeSbxwhzpNRz5GDjUN5cpREOUeC8mRQQpUn3SieELsvNzJ7b9eTOdyj+8TkI3ot8xV6DrDwO+/+EfQcYGE9wQ9ClE/8x5A5+XpOr+mVQKF3vwpYyyiclJ1nH6Zl3f8IsNy7PwHY4PWl0rs/IWzxe8sdSuYEcgvbNvCefxa5T+w/j54nlD8blvU/iD7MIHMO2sm9po8DSrPuPw2cHcL4p3a/dsgknbe87ZHX6y8MffDml5M5/2Vc2LaBtz7vBG7sZ5lIvBeGcZ1EOkcOog+hzZN5xh/aHJlvH7zllCdHL/7Q58g8+qA8OfKxhypPjsoGHMaVeS6ZUZzeABZ505aQ+RUUoAj4vffGfBaYndV2kdfuNeCckMb/KLAZeNG7LfWmnwys8t6Iq4DPh3gb/A/wshfrE8ChWW0/522bdcBnwxi/9/9i4Pu92oViG5D5la0B6CRzTP7ngcuAy7z5BvzU698q4Pgwrf88+3AbsD3rfbDSmz7bW/8vea+xRSGN/6tZ74HlZH1J8Hv9hbEP3jKfITMwTXa7sGyDU8mcY/OvrNfJuVF7L4zAeol0jsyzD6HOk3nEH+ocmU8fvP8Xozw5WvGHOkfm2QflyZGPP1R5svvXBhEREREREZERFeVzQEVERERERCRCVICKiIiIiIhIIFSAioiIiIiISCBUgIqIiIiIiEggVICKiIiIiIhIIFSAioxhZlZhZl8e7ThERETCSHlSJHgqQEXGtgpAiVVERMSf8qRIwFSAioxt3wcONLMXzeyHox2MiIhIyChPigTMnHOjHYOIjBAzmwXc75ybO8qhiIiIhI7ypEjwtAdUREREREREAqECVERERERERAKhAlRkbNsBlI52ECIiIiGlPCkSMBWgImOYc64R+IeZrdbgCiIiIj0pT4oET4MQiYiIiIiISCC0B1REREREREQCoQJUREREREREAqECVERERERERAKhAlREREREREQCoQJUREREREREAqECVERERERERAKhAlREREREREQC8f/tjayVNCsiOgAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6AAAAEYCAYAAABCw5uAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABWV0lEQVR4nO3dd3xUVfrH8c8z6YEQeg1NRXo1gmJDsYFiWV1su2vBdbGhbrUtIq6rP8uirm1RUdeOHVys2AtKEAiIqICUhN5CS8/5/TEDTpJJASZzZ5Lv+/XKKzO3zH3uzCRnnrnnPMecc4iIiIiIiIjUNZ/XAYiIiIiIiEjDoARUREREREREIkIJqIiIiIiIiESEElARERERERGJCCWgIiIiIiIiEhFKQEVERERERCQilICKNEBmdpGZfe51HCIiItHOzJyZHeR1HCL1hRJQkRhlZteZ2TIz22Zmq81skpnFex3XbmbW2sxeCMSWZ2ZfmNkQr+MSEZGGxcwSzewVM1seSCaHeR1TRWb2FzNbaGbbzexnM/uL1zGJ1BUloCKxaxowyDnXBOgD9AfGeRtSOY2B2cAhQHPgaeB/ZtbY06hERKQh+hz4DbDW60CqYMDvgGbAycBVZnautyGJ1A0loCJ1wMyuN7OlgW8yF5nZmeE+hnNuqXNu6+5DAmXA3nQRMjN7MHB1crGZDQ9zfMucc/9yzq1xzpU65yYDiUD3cB5HRERiRyTax4qcc0XOufucc58Dpfv4MCMDvY42mtndZhbWz9DOubucc98650qccz8AbwJHhPMYItFCCahI3VgKHAWkA7cCz5pZu1Abmtn5Zra1mp9OVR0ksO82YCP+K6D/2YsYhwTibAncArxmZs2rOM5b1cT3Vm0OZmYD8CegS/YiRhERqV8i0j7WgTOBTGAQcDpwSV3FbGaG/zn6LpwnIBItzDnndQwi9Z6ZzQNucc69WUeP3w1/152HnHM1di8ys4uAfwIdXOCfgJl9A/zbOfdMHcTXBPgCeN45d0e4H19ERGJTXbePIY6XA/zGOffxXuzjgBHOuXcC968AznLOhbXnUNDxbgXOAAY75wrr4hgiXtIVUJE6YGa/M7N5u7/xxD9Gs2VdHc859xP+b0of3ovdcl35b6BWAO3DGhhgZinAdGCWkk8RkYatrttHM+tkZjt2/4TrcYFVQbfrpL0EMLOr8H+hfIqST6mvlICKhJmZdQYeA64CWjjnmgIL8Y/TDLX9BcGNZYif2nYxigcO3ItQOwS6+ezWCVhdRYxvVxPf21UdwMySgDeAHOAPexGbiIjUM5FoH51zK51zjXf/hDH8jkG3q2sv97lNN7NLgOuB4c65nDDGLhJVlICKhF8jwAEbAMzsYvzf8IbknHsuuLEM8bMy1H5mdqmZtQ7c7gXcAMwMWv+xmU2oJs7WwDgzSzCzXwM9gRlVxDiimvhGVBFfAvAKkA9c6JwrqyYWERGp/yLSPoZiZklmlhy4m2hmybu/hDX/3NjLa3iIv5hZMzPrCFwDvBTOmM3sAvxDY05wzi2r7XmJxCIloCJh5pxbBNwLfAWsA/riH/8YbkcAC8xsJ/7EcQZwY9D6jjUc92ugG/4CRrcDZzvnNoUxvqHAqcCJwNagb3+PCuMxREQkRkSwfQzlB/xfiHYA3g3c7hxYV1N7Cf6qtHOAecD/gCfCHN8/gBbA7KD28tEwH0MkKqgIkUg9ZGYZwFTn3FCvYxEREYlmZvYecI1z7nuvYxFpCJSAioiIiIiISESoC66IiIiIiIhEhBJQERERERERiQgloCIiIiIiIhIR8V4HEErLli1dly5dvA5DRESkSnPmzNnonGsV6eOqjRQRkVhQVTsZlQloly5dyMrK8joMERGRKpnZCi+OqzZSRERiQVXtpLrgioiIiIiISEQoARUREREREZGIUAIqIiIiIiIiERGVY0BFRBqS4uJicnJyKCgo8DoUCSE5OZmMjAwSEhK8DqVKeg+FTyy83iIisUwJqIiIx3JyckhLS6NLly6YmdfhSBDnHJs2bSInJ4euXbt6HU6V9B4Kj1h5vUVEYlmNXXDNrKOZfWRmi8zsOzO7JsQ2ZmYPmNkSM8s2s0FB6y40s58CPxeG+wSqlD0VJvWBCU39v7OnRuzQIiJ7o6CggBYtWihxiEJmRosWLaL+yqLeQ+ERK6+3iEhYRThvqs0V0BLgT865b80sDZhjZu875xYFbTMC6Bb4GQI8Agwxs+bALUAm4AL7TnPObQnrWVSUPRWmj4PifP/9vFX++wD9RtfpoUVE9oUSh+gVK69NrMQZ7fQ8ikiD4kHeVOMVUOfcGufct4Hb24HvgQ4VNjsd+K/zmwU0NbN2wEnA+865zYGk833g5LCeQSgzJ/7yJO5WnO9fLiIiIiIiIp7kTXtVBdfMugADga8rrOoArAq6nxNYVtXyUI99mZllmVnWhg0b9iasyvJy9m65iEgDtmrVKo499lh69epF7969uf/++0Nu55zj448/5uOPP8Y5V6cxzZ07lzFjxoRcN2bMGPr370+/fv04++yz2bFjR8hYly9fzlNPPRXyMebNm8eMGTP23H/rrbcYP358WGJviKL9PfTxxx/z5Zdf7ln34IMPMmXKlDo9vohITPAgb6p1AmpmjYFXgWudc9vCHYhzbrJzLtM5l9mqVav9e7D0jL1bLiLSgMXHx3PvvfeyaNEiZs2axUMPPcSiRYvKbZOfn89FF13Ed999x8KFC7nooovIz8+v4hH33z//+U/GjRsXct2kSZOYP38+2dnZdOrUiQcffLDSNmPHjuXzzz9n5cqVjBkzhtzc3HLrKyagp5xyCtOnT2fXrl3hPZEGItrfQxUT0EsuuYR///vfdXZsEZFY4dJDXhus07ypVgmomSXgTz6fc869FmKTXKBj0P2MwLKqltet4eMhIaX8soQU/3IRESmnXbt2DBrkrx2XlpZGz549KyVsKSkpPPLII0yZMoUnn3ySRx55hJSU8v9nd+7cySWXXMLgwYMZOHAgb775JgDXXHMNEyf6u/K8++67HH300ZSVlXHRRRcxduxYMjMzOfjgg3nrrbcA2L59O9nZ2fTv3z9kvE2aNAH8V9Py8/NDjtl7+OGHeeGFF5gyZQp33HEHHTr80sAWFRUxfvx4XnrpJQYMGMBLL72EmTFs2LA9Mcjeieb30PLly3n00UeZNGkSAwYM4LPPPiM1NZUuXbrwzTff1PVTIyIStZxzvND4Yna5xPIr6jhvqrEIkflb9ieA751z/6pis2nAVWb2Iv4iRHnOuTVm9i7wTzNrFtjuROCGMMRdvd0DZmdOpCwvhw3Wkjaj/qkCRCIS9W6d/h2LVoe3k0mv9k24ZVTvWm27fPly5s6dy5AhQ8otz8/P58orr+Tiiy8G4Morr+Thhx8ul0DcfvvtHHfccUyZMoWtW7cyePBgjj/+eO644w4OPfRQjjrqKMaNG8eMGTPw+Xx7jvfNN9+wdOlSjj32WJYsWUJWVhZ9+vSpNs6LL76YGTNm0KtXL+69995K66+66irOO+88li1bxk033cStt95K+/btAUhMTGTixIlkZWWVu3qamZnJZ599xujRsd1W6D1U/j3UpUsXxo4dS+PGjfnzn/+851i7X+/Bgwfv13MjIhKrHv1kGf+3tCfpfW/mlPWP+bvdpmf4k886zJtqUwX3COC3wAIzmxdYdiPQCcA59ygwAxgJLAF2ARcH1m02s9uA2YH9JjrnNoct+ur0Gw39RvPMl8u5Zdp3vN/maLpF5MAiIrFpx44dnHXWWdx33317rjLulpKSwpQpU/jkk08Af/JQ8crje++9x7Rp07jnnnsA/9QgK1eupGfPnjz22GMcffTRTJo0iQMPPHDPPqNHj8bn89GtWzcOOOAAFi9ezJo1a6hpKMaTTz5JaWkpV199NS+99NKepGa3hx9+mBUrVlBSUlLrsZ2tW7dm9erVtdpWQoul91Dr1q1ZvHhxOE5bRCTmzFiwhv97ZzGj+rdn5LkjofJMm3WmxgTUOfc5UG1NcuevJHBlFeumAJ6N9B/Rpy0Tpn/H/xas4do2aV6FISJSK7W9yhRuxcXFnHXWWVxwwQX86le/CrnN7m6qVXHO8eqrr9K9e/dK6xYsWECLFi0qJXgVExAzIyUlpdw8jCeddBLr1q0jMzOTxx9/fM/yuLg4zj33XO66665KCaiZ0aVLFy666KIq462ooKCgUpfQWKT3UOX3UCj15fUWEdlbc1Zs5tqX5jGoU1PuPrtfxKef2qsquLGodZNkDu3SnBkL1ngdiohIVHLOMWbMGHr27Mkf//jHfX6ck046iX//+997qpvOnTsXgBUrVnDvvfcyd+5c3n77bb7++pdC6i+//DJlZWUsXbqUZcuW0b17d3r27MmSJUv2bPPuu+8yb948Hn/8cZxze9Y555g2bRo9evTY61jT0tLYvn17uWU//vhjjV1/JbRofw/p9RYR8Vu2YQeXPp1Fh6YpPH7hoSQnxEU8hnqfgAKc0rcdP67bwU/rtte8sYhIA/PFF1/wzDPP8OGHHzJgwAAGDBhQrkJsbf3973+nuLiYfv360bt3b/7+97/vSUzuuece2rdvzxNPPMGll1665+pUp06dGDx4MCNGjODRRx8lOTmZHj16kJeXVylhAH+ic+GFF9K3b1/69u3LmjVr9mn6lGOPPZZFixbtKUIE8NFHH3HKKafs9WNJ9L+HRo0axeuvv76nCNHumE844YTwPQkiIlFu445CLnpyNj4znrr4UJo3Sqx5pzpQmzGgMU/dcEVEqnbkkUeGZU7GlJQU/vOf/1Ra/sEHH+y5fcghh7BgwYI9948//ngeffTRSvtccsklvPTSS1x66aXllvt8Pr744ov9jrV58+bMnj17z/1169aRn59P37599/uxG6Jofw8dfPDBZGdn71k3d+5cevfuTYsWLfY7ZhGRWJBfVMqlT2exfnsBL/z+MDq3aORZLA3iCujubrj/y1Y3XBGRWHD55ZeTlJQUseOtXLkyZDXdcDKzKWa23swWhlj3JzNzZtayToNoQKp7D23cuJHbbrstwhGJiHijtMxxzYtzmZ+zlfvPHcjATs1q3qkONYgroODvhnvLtO/4ad12uukqqIiI55566qkq1yUnJ/Pb3/42YrEceuihkTjMU8CDwH+DF5pZR/zTlK2MRBD1yb6+h9T1VkQaCucct721iPcWrWPCqF6c1Lut1yE1jCug4O+GawZv6SqoiIh4wDn3KRBqKrJJwF+B/e/DKiIiEuSJz3/mqS+Xc+mRXbnoiK5ehwM0oAS0dZNkBndpzvT5q8MyTkVERGR/mdnpQK5zbr7XsYiISP0yY8Eabp/xPSP6tOXGkT29DmePBpOAApw+oAPLNu5kYe42r0MREZEGzsxSgRuBGsv4mtllZpZlZlkbNmyo++BERCQ2ZU+FSX1wE5rS/5UjGddqLpPOGYDPF9m5PqvToBLQkX3bkhBnvDEv1+tQREREDgS6AvPNbDmQAXxrZpUG6DjnJjvnMp1zma1atYpwmCIiEhOyp8L0cZC3CsPRwTZy7a4HSf7+Va8jK6dBJaBNUxMZ1r010+evprRM3XBFRGry8ssv07t3b3w+H1lZWV6HU6845xY451o757o457oAOcAg59xaj0MLq7/85S/06NGDfv36ceaZZ7J161avQxIRqZ9mToTi/HKLrCTfvzyKNKgEFOCMAR1Yv72QWcs2eR2KiMi+CXSvYUJT/+/sqXV2qD59+vDaa69x9NFH19kxGgozewH4CuhuZjlmNsazYCL4HjrhhBNYuHAh2dnZHHzwwdxxxx11diwRkYbM5eWEXlHVco80uAR0eM/WNE6K54256oYrIjEoqHsNOP/v6eP2O4FYvnw5ffr02XP/nnvuYcKECfTs2ZPu3bvvZ9AC4Jw7zznXzjmX4JzLcM49UWF9F+fcxjoPJMLvoRNPPJH4eP+sb4cddhg5OdH1QUhEpD7YXlDMel8VU0mnZ0Q2mBo0uAQ0OSGOk3q35Z2FaykoLvU6HBGRvROiew3F0de9RqKYh++hKVOmMGLEiDo/johIQ1JQXMpl/53DnUWjKY1LLr8yIQWG11jrLqIaXAIKcMbA9mwvLOGjxeu9DkVEZO/ESPcaiWIevYduv/124uPjueCCC+r0OCIiDUlpmeO6l+bx1bJNHH3WFcSd/m9I7wiY//eoB6DfaK/DLCfe6wC8MPTAlrRsnMQb83IZ0bed1+GIiNReekag62SI5fshPj6esrKyPfcLCgr26/EkinnwHnrqqad46623mDlzJmbRMxWAiEgsc84x/s2FvL1wLTef0pMzB2YAo6Mu4ayoQV4BjfMZo/q346PFG8jbVex1OCIitTd8vL87TbAwdK9p06YN69evZ9OmTRQWFvLWW2/t1+NJFIvwe+idd97hrrvuYtq0aaSmpu7XMURE5BeT3v+R575eydhjDuTSow7wOpxaqzEBNbMpZrbezBZWsf4vZjYv8LPQzErNrHlg3XIzWxBYF1X1+88Y0IGi0jLeXrjG61BERGqv32h/d5owd69JSEhg/PjxDB48mBNOOIEePXoA8Prrr5ORkcFXX33FKaecwkknnRSGkxBPRfg9dNVVV7F9+3ZOOOEEBgwYwNixY8NwEiIiDdvkT5fywIdLGJ2Zwd9Ojq1igeZc9fNhmtnRwA7gv865PjVsOwq4zjl3XOD+ciBzb6v6ZWZmurqeb845x/B/fULLRklMHXt4nR5LRKQ633//PT179vQ6DKlGqNfIzOY45zIjHUuoNlLvofDS8yki0ey5r1dw0+sLOaVvOx44byBxvugc2lBVO1njFVDn3KfA5loe5zzghb2MzRNmxtmHZPDN8s0s37jT63BERERERESq9cbcXG5+YyHH9WjNpHMGRG3yWZ2wjQE1s1TgZODVoMUOeM/M5pjZZTXsf5mZZZlZ1oYNG8IVVrV+NTADn8Gr36p6pIiIiIiIRK/3vlvLn16ez5CuzXn4gkEkxsdmOZ9wRj0K+MI5F3y19Ejn3CBgBHBloDtvSM65yc65TOdcZqtWrcIYVtXapidzZLdWvDonh7Ky6rsii4jUpZqGQ4h3YuW1iZU4o52eRxGJRp/9tIGrnp9L3w7pPH7hoSQnxHkd0j4LZwJ6LhW63zrncgO/1wOvA4PDeLyw+PUhGazOK+CrZZu8DkVEGqjk5GQ2bdqkD75RyDnHpk2bSE5OrnljD+k9FB6x8nqLSAOQPRUm9YEJTSm8uxdvPnMfB7RqxFMXH0rjpNieSTMs0ZtZOnAM8JugZY0An3Nue+D2icDEcBwvnE7o1Ya05HhezlrFEQe19DocEWmAMjIyyMnJIVLDD2TvJCcnk5Gxf3Nk1jW9h8InFl5vEannsqfC9HFQnA9A0s5c/uF7jKKhvWmSmuhxcPuvxgTUzF4AhgEtzSwHuAVIAHDOPRrY7EzgPedccDWfNsDrgQmn44HnnXPvhC/08EhOiOO0/u159dscthUU0yQ5weuQRKSBSUhIoGvXrl6HITFM7yERkXpk5sQ9yeduyRSS/MU/YfD5HgUVPjUmoM6582qxzVPAUxWWLQP672tgkXT2IRk89/VKZmSv4dzBnbwOR0REREREGqq8KgqkVrU8xsRm6aQwG9CxKQe1bswrc+rHiyoiIiIiIrGpuHH70CvS68fwACWg/DInaNaKLSzbsMPrcEREREREpAFavnEnE/PPJp8KYz0TUmD4eG+CCjMloAFnDuxAnM+YmqWroCIiIiIiElkrNu3kvMdm8T+OYsvweyG9I2D+36MegH6jvQ4xLGK7hm8YtWmSzHE9WvPKnFX88YSDY3ZiVxERERERiS2rNu/ivMmzyC8u5YXfH0b7dk3gqN95HVadUJYV5NrW83ijaCwJ/2jun3cne6rXIYmIiIiISD2Ws2UX506exc6iUp67dAg92zXxOqQ6pSugu2VPpdecmzFfoORx3ir//DtQby53i4iIiIhI9Fi9NZ/zHpvF9oJinv/9YfRun+51SHVOV0B3mzkRqzDfDsX5/nl4REREREREwmh38rl1ZzHPjBlCnw71P/kEXQH9RT2fb0dERERERKLDqs27OP9xf/L59JjB9O/Y1OuQIkZXQHeral6dejLfjoiIiIiIeG/Fpp2c85+vyNtVzLOXDmFQp2ZehxRRugK62/Dx/jGfQd1wS+KSia8n8+2IiIiIiIgHsqf6h/Xl5VDcuD2P559NvjuSFy5rGGM+K9IV0N36jfbPr5PeEYexhlZMTr9GBYhERCQszGyKma03s4VBy+42s8Vmlm1mr5tZUw9DFBGRcMue6r/IlbcKcCTsyOXG0kf437FrG2TyCUpAy+s3Gq5biE3YytQj3+au1f1ZuWmX11GJiEj98BRwcoVl7wN9nHP9gB+BGyIdlIiI1KGZE8v1sARIoYj2WXd5FJD3lIBW4ZxDOxLnM577eoXXoYiISD3gnPsU2Fxh2XvOuZLA3VmACg+IiNQnKnRaiRLQKrRNT+ak3m14cfYq8otKvQ5HRETqv0uAt0OtMLPLzCzLzLI2bNgQ4bBERGRfFTVqH3pFAy50qgS0Ghce3oW8/GKmzc/1OhQREanHzOwmoAR4LtR659xk51ymcy6zVatWkQ1ORET2yVdLN3Hzjl9RQFL5FQkp/gKoDZQS0GoM7tqcHm3TeOrLFTjnvA5HRETqITO7CDgVuMCpsRERqRdmfr+OC5/8hrnpJ1Bw8iRI7wiY//eoBxp0odMap2Exsyn4G8b1zrk+IdYPA94Efg4ses05NzGw7mTgfiAOeNw5d2d4wo4MM+PCoV244bUFZK3YwqFdmnsdkoiI1COBdvKvwDHOOVW9ExGpB96cl8ufps6nV/smPHXxYJo2SoTDLvA6rKhRmyugT1G5al9FnznnBgR+diefccBDwAigF3CemfXan2C9cPqA9jRJjuepL5d7HYqIiMQwM3sB+ArobmY5ZjYGeBBIA943s3lm9qinQYqIyH55ZtYKrn1pHod0bsZzlw6heaNEr0OKOjVeAXXOfWpmXfbhsQcDS5xzywDM7EXgdGDRPjyWZ1IT4znn0I5M+WI5a/MKaJue7HVIIiISg5xz54VY/ETEAxERkTrx8MdLuOudHxjeozUPXTCI5IQ4r0OKSuEaA3q4mc03s7fNrHdgWQdgVdA2OYFlIUVzhb/fHtaFMud4XlOyiIiIiIhIEOccd769mLve+YHTB7Tn0d8eouSzGuFIQL8FOjvn+gP/Bt7YlweJ5gp/nVqkclz31jz/zUoKSzQli4iIiIiIQGmZ46Y3FvLoJ0v5zWGdmDR6AAlxqvNanf1+dpxz25xzOwK3ZwAJZtYSyAU6Bm2aEVgWky4c2oWNO4qYPn+N16GIiIiIiIjHCktKuebFuTz/9UquGHYgt53eB5/PvA4r6tU4BrQmZtYWWOecc2Y2GH9SuwnYCnQzs674E89zgfP393heOapbS7q3SePxz5Zx1qAOmOnNJSIiIiLSoGRPhZkTcXk5bItrhS//bG4YcSl/OOZAryOLGbWZhuUFYBjQ0sxygFuABADn3KPA2cDlZlYC5APnBuYxKzGzq4B38U/DMsU5912dnEUEmBmXHtWVv7ySzWc/beTog6Orm7CIiIiIiNSh7KkwfRwU52NAq9L1/CtlCvHNBgJKQGurNlVwQ1XtC17/IP4y8qHWzQBm7Fto0ee0Ae25690feOyzZUpARUREREQakpkToTi/3KL40gL/8n6jPQoq9miE7F5Iio/joqFd+OynjXy/ZpvX4YiIiIiISIS4vJzQK6paLiEpAd1LFwzpREpCHI9/9rPXoYiIiIiISAR8uHgdq12L0CvTMyIbTIxTArqXmqYmMjozg2nzc1m3rcDrcEREREREpA5NzVrF7/87h+cbX4iLTym/MiEFho/3JrAYpQR0H1xyZFdKyxxPfbnc61BERERERKQOOOd46KMl/PWVbIYe2ILLx92InfYApHcEzP971AMa/7mX9nsaloaoc4tGnNS7Lc/NWsGVxx5E4yQ9jSIiIiIi9UVxaRk3v76Ql7JWccaA9tx1dn8S433+ZFMJ537RFdB99IdjDmRbQQnPf73C61BERERERCRMthUUc/GTs3kpaxVXH3cQk84Z4E8+JSz0TO6jAR2bcsRBLXjss58pKC71OhwREREREdlPOVt2cfYjXzJr2SbuOrsffzqxO2bmdVj1ihLQ/XDlsQexYXshL89R6WURERERkVg2f9VWznjoS9bkFfDfSwYzOrOj1yHVS0pA98PhB7RgYKem/OeTpRSXlnkdjoiIiIiI7IN3v1vLOZO/IjnBx+tXDGXoQS29DqneUvWc/WBmXHXsQbz5zH0U3XMFCflr/fMADR+vwckiIiIiItEqeyrMnIjLy2FHcltmbD+T7u1P5fHfZdIqLcnr6Oo1JaD76bjiTzgi8QmS8wv9C/JWwfRx/ttKQkVEREREokv2VP/n9eJ8DEgrWMNdiU9gh/cnMe0Ir6Or99QFdz/ZzIkkU1h+YXE+zJzoTUAiIiIiIlK1mRP9n9eDJLlCEj/5h0cBNSxKQPdXXhUFiKpaLiIiIiIinnH6/O4pJaD7Kz1j75aLiIiIiIgnPvlxA2toEXqlPr9HhBLQ/TV8PCSklFvkElL8y0VERERExHPOOZ74/GcufvIbnkm9kLL48p/f0ef3iFECur/6jYZRD0B6RxxGTllLsgdOVAEiEREREZEoUFhSyt9ezea2txZxYq+2XHXNjfhO839+B/P/HvWAPr9HSI1VcM1sCnAqsN451yfE+guAvwEGbAcud87ND6xbHlhWCpQ45zLDF3oU6Tca+o2mtLSM3933KQk/+Hj7ZIfPZ15HJiIiUSJUe2pmzYGXgC7AcmC0c26LVzGKiNQ3G3cUcvmzc5i9fAvjhnfj2uHd/J/RA5/fJfJqcwX0KeDkatb/DBzjnOsL3AZMrrD+WOfcgHqbfAaJj/NxzfBu/LBuOzMWrvE6HBERiS5PUbk9vR6Y6ZzrBswM3BcRkTBYmJvH6Q9+wYLcPB48fyB/POFgXSCKAjUmoM65T4HN1az/Mujb2llAgx69e2q/9nRr3Zj7P/iJ0jLndTgiIhIlqmhPTweeDtx+GjgjkjGJiNRXr8/N4axHvqTMOV7+w1BO7dfe65AkINxjQMcAbwfdd8B7ZjbHzC6rbkczu8zMsswsa8OGDWEOK3LifMa1xx/MT+t38Fb2aq/DERGR6NbGObe7y8xaoE2ojepLGykiUteKS8uYMO07rntpPgM7NWX61UfSNyPd67AkSI1jQGvLzI7Fn4AeGbT4SOdcrpm1Bt43s8WBb4Arcc5NJtB9NzMzM6YvHY7o05YebdO474OfGNm3HQlxqvUkIiLVc845MwvZ/tWnNlJEJOyyp8LMibi8HLb6WrG54GwuPfK3XD+iB/H6HB51wvKKmFk/4HHgdOfcpt3LnXO5gd/rgdeBweE4XrTz+Yw/n9idnzfuZGrWKq/DERGR6LXOzNoBBH6v9zgeEZHYkj0Vpo+DvFUYjlZl6/lXyhRu7rRQyWeU2u9Xxcw6Aa8Bv3XO/Ri0vJGZpe2+DZwILNzf48WK4T1bc2iXZtz3wU/sKirxOhwREYlO04ALA7cvBN70MBYRkdgzcyIU55dbFF9a4F8uUanGBNTMXgC+ArqbWY6ZjTGzsWY2NrDJeKAF8LCZzTOzrMDyNsDnZjYf+Ab4n3PunTo4h6hkZlw/ogcbthcy5fOfvQ5HREQ8Fqo9Be4ETjCzn4DjA/dFRKQWCopLcXk5oVdWtVw8V+MYUOfceTWsvxS4NMTyZUD/fQ8t9h3SuTkn9mrDo58s4/whnWneKNHrkERExCPVtKfDIxqIiEg9sHLTLq54fg6PlrUgw7ex8gbpDXpijqimjtF17K8nd2dXUQkPfrjE61BERERERGLeOwvXcsq/P2PV5ny2HH49JKSU3yAhBYaP9yY4qZES0Dp2UOs0Rmd25JlZy1m1eZfX4YiIiIiIxKSikjJue2sRY5+dwwEtG/HW1UfSd8TvYdQDkN4RMP/vUQ9Av9FehytVCNs0LFK1a48/mDfm5XLPez9w/7kDvQ5HRERERCSm5G7N56rnv2Xuyq1cNLQLN4zsQVJ8nH9lv9FKOGOIroBGQNv0ZC498gDenLeab1du8TocEREREZGY8dHi9ZzywGf8tG4HD50/iAmn9f4l+ZSYowQ0Qi4fdiCt05K4dfoiyso0h7iIiIiISHVKSsu4653FXPzUbNqlpzD96iM5pV87r8OS/aQuuBHSKCmev53cgz+9PJ835uXyq0GqzCUiIiIiskf2VP/8nXk5lKR14CHf+Ty8bgDnHtqRCaf1JjlBVz3rA10BjaAzB3agf0Y6//fOYnYWlngdjoiIiIhIdMieCtPHQd4qwBG/PYfLtt7HS4ev4s6z+in5rEeUgEaQz2eMH9WLddsKefSTpV6HIyIiIiISHWZOhOL8cotSrIghyx70KCCpK0pAI+yQzs05rX97Jn+6jJwtmpZFRERERMTl5YReUdVyiVlKQD1w/YgenOr7nJSHBsCEpjCpj7/bgYiIiIhIA+Kc49lZK1jtWoTeIF11U+obJaAeaL9yOnfGP06LknWA8/d1nz5OSaiIiIiINBhbdxVx+bPfcvMbC3mzxRhcfEr5DRJSYPh4b4KTOqME1AszJ5JQVlB+WXG+v++7iIiIiEg9N2vZJkbe/xkzF6/jxpE9GHvVDdhpD0B6R8D8v0c9AP1Gex2qhJmmYfGC+riLiIiISANUWFLKve/9yGOfLaNz81RevXwo/TKa+lf2G62EswFQAuqF9IxAiekQy0VERERE6qHv12zjupfmsXjtds4f0ombRvakUZLSkYZGr7gXho/3j/kMKjVdaEkkDh+PeRiWiIiIiEhYZE/1Dy/Ly8GlZzCz/R+4IvsgmqQkMOWiTI7r0cbrCMUjGgPqhX6j/X3aA33ctye34y+FY3jHjvI6MhERERGR/ZM91X+xJW8V4LC8VQxdNJG/dpjPu9cepeSzgatVAmpmU8xsvZktrGK9mdkDZrbEzLLNbFDQugvN7KfAz4XhCjzm9RsN1y2ECVtJ+csilrQZya3TF7GjsMTryERERERE9t3MieV6+gGkWhFjCp+lReMkj4KSaFHbK6BPASdXs34E0C3wcxnwCICZNQduAYYAg4FbzKzZvgZbX8XH+bj9zD6s217AXe8s9jocEREREZF95qoorGkquCnUMgF1zn0KbK5mk9OB/zq/WUBTM2sHnAS875zb7JzbArxP9YlsgzWwUzMuHtqV/361gq+XbfI6HBERERGRveKcY9r81ayhRegNVHBTCN8Y0A5AcFnXnMCyqpZXYmaXmVmWmWVt2LAhTGHFlj+fdDCdmqfyt1ezyS8q9TocEREREZFaWb+9gLHPzmHcC3N5Ie1iyuJTym+QkOIvxCkNXtQUIXLOTXbOZTrnMlu1auV1OJ5ITYznzrP6snzTLv71/g9ehyMiIiIiUi3nHG/Oy+XESZ/y0Q8buGFED6659kZ8p/1ScJP0jv4CnJrjUwjfNCy5QMeg+xmBZbnAsArLPw7TMeuloQe25IIhnXji858Z2bcdAztpyKyIiIiIRJ/12wq46Y2FvL9oHQM7NeXus/tzUOvG/pX9RivhlJDCdQV0GvC7QDXcw4A859wa4F3gRDNrFig+dGJgmVTj+hE9aNskmb++kk1hibriiojUd2Z2nZl9Z2YLzewFM0v2OiYRkao453h9bg4nTPqUT3/cwE0je/LK2KG/JJ8i1ajVFVAzewH/lcyWZpaDv7JtAoBz7lFgBjASWALsAi4OrNtsZrcBswMPNdE5V10xIwHSkhO4/Vd9ufjJ2dz/wU/89eQeXockIiJ1xMw6AOOAXs65fDObCpyLvwK9iIj3sqf6p1bJy6EkrQNPJv2W23P6ckjnZtx1dj8ObKXEU2qvVgmoc+68GtY74Moq1k0Bpux9aA3bsd1bMzozg0c/WcpxPVqT2aW51yGJiEjdiQdSzKwYSAVWexyPiIhf9lSYPm7PvJ7x23P4zbZ7OXDQLRxz9kjifOZxgBJroqYIkVQ2flRvMpqlct3UeWwvKPY6HBERqQPOuVzgHmAlsAb/MJb3grdRpXgR8czMiXuSz91SrIjjch9V8in7RAloFGucFM+kc/qTuyWfW6cv8jocERGpA4EaCacDXYH2QCMz+03wNqoULyJeKCguxeXlhF5Z1XKRGigBjXKHdG7OFcMO4pU5ObyzcI3X4YiISPgdD/zsnNvgnCsGXgOGehyTiDRwXy7dyMn3fUpuWYvQG6RnRDYgqTeUgMaAa47vRt8O6dzw2gLWbyvwOhwREQmvlcBhZpZqZgYMB773OCYRaaC27CziLy/P5/zHvsYBO4+6ERJSym+UkALDx3sSn8S+cM0DKnUoIc7HpHMGcOq/P+OVp/7F5aXPY3k5/m+eho/XHEsiIjHMOfe1mb0CfAuUAHOByd5GJSINTVmZ49Vvc7jz7cXk5RdzxbADGTe8G8kJx0KbtD1VcPX5U/aXEtAYcVDrxjw+8GcGzZ+EWZF/Yd4qf1Uy0D8BEZEY5py7Bf8UZyIiEff9mm38/Y2FZK3YwiGdm/GPM/rQs12TXzboN1qfNSVslIDGkCNWPPxL8rlbcb7/Gyn9UxARERGR6gTN50l6BvlH38Tdq/vz9FfLSU9J4K6z+3H2oAx8qm4rdUgJaAwxVSETERERkX1RYT5Pf0+6a9hUPIZzDj2fv57Unaapid7GKA2CEtBYkp7h/2cRarmIiIiISFVCzedJIXc3fZPEM+/wKChpiFQFN5YMH1+pClmxL1lVyERERESkWlXN55m4c3WEI5GGTgloLOk3GkY9AOkdcRib4tvwl6IxzGo83OvIRERERCQKlZU5Xs5axVo0n6dEByWgsabfaLhuITZhK4l//o4FzU7kyue+ZfXW/Jr3FREREZEGY86KLZzx8Bf85ZVsnk+7mNI4zecp3lMCGsPSkhP4z28zKSwpY+yzcygoLvU6JBERERHx2Jq8fK59cS5nPfIl67YV8K/R/bnuupuIO93fkw7M/3vUA5pJQSJORYhi3EGtGzPpnAH8/r9Z3PT6Qu75dT/MVDpbREREpKEpKC7lsU+X8fDHSyl1jquOPYjLhx1Io6TAR37N5ylRQAloPXBCrzZce3w37vvgJ/plpHPh0C5ehyQiIiIidSloTk+XnsG8buO4auFB5G7NZ0Sfttw4sicdm6d6HaVIJUpA64lxx3VjYe42bntrEQe3SePwA6sYaC4iIiIisa3CnJ6Wt4rus29iVOrVHP37Kxh6YEuPAxSpWq3GgJrZyWb2g5ktMbPrQ6yfZGbzAj8/mtnWoHWlQeumhTF2CeLzGZPO6U+Xlo0Y++wclm7Y4XVIIiIiIlIXQszpmWpF/C3xJSWfEvVqTEDNLA54CBgB9ALOM7Newds4565zzg1wzg0A/g28FrQ6f/c659xp4QtdKkpLTuDJiw4l3mdc8tRsNu8s8jokEREREQmjzTuLqpzT0/JyIxyNyN6rzRXQwcAS59wy51wR8CJwejXbnwe8EI7gZO91bJ7K5N9lsiavgD88k0VhiSrjioiIiMS6guJSHvl4Kcfc9RG5TnN6SuyqTQLaAVgVdD8nsKwSM+sMdAU+DFqcbGZZZjbLzM6o6iBmdllgu6wNGzbUIiypyiGdm3Hvr/sze/kWrn91Ac45r0MSERERkX1QVuZ47dscjrvnY/7vncUMOaA58Sfc4p/DM5jm9JQYEe4iROcCrzjngi+7dXbO5ZrZAcCHZrbAObe04o7OucnAZIDMzExlTPtpVP/2rNi0k3ve+5HjSz7hlPWPQV6O/5ux4eNVgltEREQkijnn+HDxeu5+9wcWr91Ov4x07h09IFBo8lBokrynCq4+30ksqU0Cmgt0DLqfEVgWyrnAlcELnHO5gd/LzOxjYCBQKQGV8Lvy2INovvRNjv3x/8AC40HzVvmrpoH+SYmIiIh4LWg6ld2J5FeNhnP3u4v5duVWurRI5YHzBnJq33b4fEFzvWtOT4lRtUlAZwPdzKwr/sTzXOD8ihuZWQ+gGfBV0LJmwC7nXKGZtQSOAO4KR+BSMzPjvB1PYlahGFFxvv8fnf5piYiIiHinwnQq5K2i8PWreKFwDKsbH88dv+rL2YdkkBBXq4krRGJCjQmoc67EzK4C3gXigCnOue/MbCKQ5ZzbPbXKucCLrvyAw57Af8ysDP940zudc4vCewpSnSqroVVRPU1EREREIiTEdCpJrpDbm7xOwp/+QXJCnEeBidSdWo0Bdc7NAGZUWDa+wv0JIfb7Eui7H/HJ/krP8He7DbVcRERERDzj8nKwEMvTCtaCkk+pp3Q9v74bPr5SlbR8ElmT+VePAhIRERFp2HK27OKG1xawWtOpSAOkBLS+6zcaRj0A6R0BozitA//0Xc6oT9qzdMMOr6MTERERaTB2J57H3vMxr87J4YvOV1AWr+lUpGEJ9zQsEo2CqqQlABeu387bk2dx/mOzeOmyw+nSspG38YmIiIjUY6s27+Lhj5fwclYOPjPOPbQTlw87kPZNR0B2R02nIg2KEtAG6KDWaTx36WGcO/krfxL6h8Pp2DzV67BERBosM2sKPA70ARxwiXPuq2p3EpHoU2FKlU2HXc/dq/vxyhx/4nn+EH/i2S496KqnplORBkZdcBuo7m3TePbSIewsKuW8x2aRuzW/5p1ERKSu3A+845zrAfQHvvc4HhHZW7unVMlbBTjIW0XKO9dRNPclLhjSiU/+OoyJp/cpn3yKNEBKQBuw3u3TeXbMEPLyixn96Fes2LTT65BERBocM0sHjgaeAHDOFTnntnoalIjsvRBTqqRaEXc1e4NblXiK7KEEtIHrm5HO85cexq6iEn796Ff8tG671yGJiDQ0XYENwJNmNtfMHjezcoPzzewyM8sys6wNGzZ4E6WIVGnuyi24KuZYj99exZzsIg2UElChb0Y6L/3hcBxwzuRZLMzN8zokEZGGJB4YBDzinBsI7ASuD97AOTfZOZfpnMts1aqVFzGKSAXOOT7/aSPnPzaLMx/+kjVoShWR2lACKgAc3CaNl/9wOCkJcZz32CzmrNjsdUgiIg1FDpDjnPs6cP8V/AmpiEShsjLHOwvXcsZDX/CbJ75myfod3DSyJ81Pu73S3OuaUkWkMlXBlT26tGzE1LGH85vHv+bFJ+6ld+PXSN61RiXBRUTqkHNurZmtMrPuzrkfgOHAIq/jEmnwKlS0LTn277xZegSPfLKUJet30LlFKv88sy9nHdKBpPg44ACI92lKFZEaKAGVcjo0TeGNo3NJmvEYybsK/QvzVvmruoH+iYqI1I2rgefMLBFYBlzscTwiDdvuira7iwrlraL4jav5pGgM8a1H8MB5AxnZpy3xcRU6E2pKFZEaKQGVStK/uAMoLL+wON//jZ7+qYqIhJ1zbh6Q6XUcIhIQoqJtCoXc1fQNkq75J2bmUWAisU9jQKWyKqq4VVXdTURERKS++LaairbJu9Yo+RTZT7oCKpWlZwQmUS5vU1wrkgtLaJykt42IiIjUH6VljvcXreWxz35mzootfJncgvZsrLyhKtqK7DddAZXKho+vVMWtxJfMPwp/zVkPf8mqzbs8CkxEREQkfLYVFPPE5z9z3L0fM/bZb1m3rYBbRvVSRVuROqRLWVLZ7nGeQVXc4oeP5+yU47jiuTmc9uDnPHzBIRx+YBXzXYmIiIhEiwrVbBk+niVtR/D0lyt49dscdhWVMqhTU/56Ug9O6t0mUFioqyraitQRc87VvJHZycD9QBzwuHPuzgrrLwLuBnIDix50zj0eWHchcHNg+T+cc0/XdLzMzEyXlZVV23OQCPp5404ufXo2Kzbt4tbTe3PBkM5ehyQi4gkzm+Oci3jhILWRInuhYjVboNCS+EvhGN6xoxnVvz0XDe1C34x0D4MUqZ+qaidrvAJqZnHAQ8AJ+CfLnm1m05xzFecoe8k5d1WFfZsDt+Cv7OeAOYF9t+zjeYjHurZsxOtXHsG4F+Zy0+sL+W71Nsaf2ovkhDivQxMREREpL0Q12yRXyD/SXmP81RNo2TjJo8BEGq7ajAEdDCxxzi1zzhUBLwKn1/LxTwLed85tDiSd7wMn71uoEi2aJCfwxIWH8odjDuD5r1fy60e/0rhQERERiSo/rN1eZTXbJoXrlHyKeKQ2CWgHILgkak5gWUVnmVm2mb1iZh33cl/M7DIzyzKzrA0bNtQiLPFSnM+4YURPJv/2EJZv2skpD3zGB4vW+bu6TOoDE5r6f2dP9TpUERERaSAKikt5dU4OZz3yJSfd9ymrXRX1KlTNVsQz4SpCNB14wTlXaGZ/AJ4GjtubB3DOTQYmg398S5jikjp2Yu+2/K9tEy5/bg7Tnr2PY5KnkFBW4F+Zt8o/7gI0aF9ERETqzJL1O3j+65W8+m0OefnFHNCyETef0pP05NvgvT+W74ararYinqpNApoLdAy6n8EvxYYAcM5tCrr7OHBX0L7DKuz78d4GKdGtU4tUXr18KPl3jSGhuKD8yuJ8//gLJaAiIiKyv4Iq2rr0Dsw5aBx3r+7H1z9vJiHOOKl3W84f0onDD2iBmQEHQFK8qtmKRJHaJKCzgW5m1hV/QnkucH7wBmbWzjm3JnD3NOD7wO13gX+aWbPA/ROBG/Y7aok6yQlxJBevD72yivEXIiIiIrVWoaKt5eXQK+tmeiVewbCTf8evMzNCj+vsN1oJp0gUqTEBdc6VmNlV+JPJOGCKc+47M5sIZDnnpgHjzOw0oATYDFwU2Hezmd2GP4kFmOic21wH5yHRID3D3+22grImHWo12FhEREQklPyiUtzb40mtUNE21YoYn/oqNuxWbwITkb1WqzGgzrkZwIwKy8YH3b6BKq5sOuemAFP2I0aJFcPHV5pra5dL5P92ns1JSzcy9MCWHgYnIiIiscQ5x7xVW5malcNb81cznzVglbcz9bQSiSnhKkIk8kv3lqBxFhsG/JlP53Th6ce+5oIhnbh+RA/SkhO8jVNERESi1obthbwxN5epWav4af0OkhN8jOzbjqIV7UneubryDqpoKxJTlIBKeFUYZ9EZmHFEKf96/wee+PxnPlq8nn/+qi/Durf2LkYRERHxVlAxIdIzKDn273yUOIypWav4aPF6Ssocgzo15Y5f9eXUfu38X15n31qpp5Uq2orEHiWgUudSEuO46ZRejOzbjr++ks1FT87mV4M6MP7UXjRNTfQ6PBEREYmkCsWEyFtF8RtXM71oDHNThzPmyK78OjODg1qnld8vRE8rVbQViT3mXPRNuZmZmemysrK8DkPqQGFJKQ9+uIRHPl5KekoC14/owVkJX+H7UI2JiMQWM5vjnMuM9HHVRkqsK7m3N/HbK4/bzE9tT/yfviMhTqULReqDqtpJ/YVLRCXFx/GnE7sz7aoj6dwilc9ee5ii168KVM91/t/Tx/m/HRUREZF6IS+/mJdmr+Sc/3yFb1vookEpu9Yo+RRpANQFVzzRq30TXhk7lPy7x5CcX1h+ZXG+v3uNroKKiIhEvwrjOXf3ZCoqKeOTHzfw+twcPvh+PUUlZRzQshE7ktvSpHBt5cdRMSGRBkEJqHjG5zMa5YdogACXlxOq0rqIiIhEkxDjOUvfHMers1dxR25ftuwqpkWjRM4f3IkzB3agX0Y6tuA2FRMSacCUgIq30jMC3W/LW28t+f6H9aqWKyIiEs1mTiyfSAJxpfkcsfJhjuzxCr8a2IEju7Us37VWxYREGjQloOKt4eMrfQtaEpfMf3y/YcqTszmqW0tuHNmTnu2aeBikiIiIBHPOMT8nj/5V9Fhqb5v493kDq36ACtO2iUjDoQRUvBXiW9D44eO5vtfZdJi1ggdm/sTIBz7jzIEduGZ4Nzq3aORtvCIidcTM4oAsINc5d6rX8YhU5Jzju9XbeCt7Df9bsJpVm/P5IqkFHWxjpW1N4zlFpApKQMV7Ib4FTQTGHNmVswdl8NDHS3j6y+W8OW81vz4kg6uOO4iMZqlVFj0QEYlR1wDfA+ryId6p0La64eNZ0PxE3lm4lrcXruXnjTuJ9xlHHNSSccd1oxn/gHev03hOEak1JaAS1dJTE7hxZE8uPbIrD3+8lOe/Xsmr3+Zw+4Hfc/bqu/GV/FL0gOnj/LeVhIpIjDGzDOAU4Hbgjx6HIw1ViIJCBa9dxeNFY/gfR3H4AS247OgDOLl3W5o1SgzsdB4kxukLYRGpNXPOeR1DJZpkW6qyJi+fhz5awuVzzwjZ5Yf0jnDdwsgHJiINTlUTbO/jY70C3AGkAX+urguu2kipC4UlpZT9qw8pu1ZXWrczpR3FV2fTNDUxxJ4iIqFV1U7qCqjElHbpKfzjjL64eZtCb5AXenJrEZFoZWanAuudc3PMbFgV21wGXAbQqVOnyAUn9drOwhI++XED7yxcy0eL1zOf1YSqKNQofy0o+RSRMFECKjHJqpi+ZUNcK77/cQNHdWuJmWYSFZGYcARwmpmNBJKBJmb2rHPuN7s3cM5NBiaD/wqoN2FKTKmiTsLavAI+XLyemd+v4/MlGyksKaN5o0RG9m1H4bL2Ia+AooJCIhJGSkAlNoWYvqXYl8wD7jyemfINB7ZqxIVDu/CrQRk0TtLbXESil3PuBuAGgMAV0D8HJ58iey3EWM7iN67mvncX89CmQwDIaJbCeYM7cVLvthzapRnxcT7IvrVS26qCQiISbrX6ZG5mJwP3A3HA4865Oyus/yNwKVACbAAucc6tCKwrBRYENl3pnDstTLFLQxZi+paE4eO5uddZDMxew9NfLmf8m99x9zs/cNYhGVw4tAtdWzZS5VwREan3yj64FV9wEgkklBVwYf4zNDr5fIb3aMPBbRpX7ikUom1VOyki4VZjEaLAvGQ/AicAOcBs4Dzn3KKgbY4FvnbO7TKzy4FhzrlzAut2OOca701QKrAg4TB35Rae/nI5/1uwhuJSx9/aZ/P7vPuILy34ZaOEFBj1gBpXEdlr4SxCtDfURkoo67b90rV28rIT8Fmoz3cGE7ZGOjQRaaD2pwjRYGCJc25Z4IFeBE4H9iSgzrmPgrafBajrkHhuYKdmDOzUjBtP6ckLX6/ijM/HEU9B+Y2K8/3f9CoBFRGRaBWi905Rr7PJWrGZT37cwCc/bGDx2u2Av2vttsTWNC1eV/lxNJZTRKJAbRLQDkBwtZccYEg1248B3g66n2xmWfi7597pnHsj1E6q8Cd1pXVaMtcc3w33eejKuS4vhx0FxaQlJ0Q4MhERkRqEGM9Z+PpV3PzKfF4uGkpCnJHZuTnXj+jBsd1b+7vWLviHxnKKSNQKa3UWM/sNkAkcE7S4s3Mu18wOAD40swXOuaUV91WFP6lrVVXOzS1rwXH/+IDhPVpz+oD2DOvemuSEOA8iFBER+UV+USm8cwspFcZzJrlCbkx6mRPPHcfhB7aoXGxPYzlFJIrVJgHNBToG3c8ILCvHzI4HbgKOcc4V7l7unMsN/F5mZh8DA4FKCahInQtROdclpFAy9GbO396Jt7JX8/bCtaQlxXNyn7acNqA9hx/QgvjvXlEjLiIi4VFNMbyyMsfitdv5culGPv1pI18v28T3caHn5mxWvJ4TerWp+jj9RqutEpGoVJsEdDbQzcy64k88zwXOD97AzAYC/wFOds6tD1reDNjlnCs0s5b45zq7K1zBi+yVEN8I2/DxdOk3mgnAzaf05Ktlm3hz3mreWbiWl+fkcEHKLG7hPyTu/k4lb5U/iQ1+PBERkdoI0Z22bNo4Zi3dxPP5Q/hy6SY27ywC4MBWjfjNYZ0pWtSeZM3NKSL1SI1VcAECk2Pfh38alinOudvNbCKQ5ZybZmYfAH2BNYFdVjrnTjOzofgT0zLAB9znnHuipuOpwp94raC4lI8Wr2fwm8fQoqRyIYeyJhn4/vidB5GJSLRQFVzZa5P6hBwKklPWkrOS/8MRB7bkiINaMvSgFrRLT/GvrJi0giq4i0hM2J8quDjnZgAzKiwbH3T7+Cr2+xJ/YioSU5IT4hjRtx28uj70Bnm5/G7KNww7uBXDureia8tGledTExGRBm9HYQmzf97M50s2clNeDr4Q23TwbWLWDcNDtyMazyki9UxYixCJ1DtVFC7altSGnM27mPjWIia+BR2bp3DMwa045uDWDD2wBY2S4qsd5yMiIjGuiv/x2wqKyVq+ma+XbWbWz5tZmJtHaZkjMd7H2KRWtCqt/MWmpWdAdV9iajyniNQjSkBFqhOicBEJKTQ99TY+7DeMlZt28clP/jnYXvs2l2dnrSQhzrim1Tz+sO1+EsoC845q7KiISP0RYixn0etXc9+7i3l08yGUOUiIM/pnNGXsMQdw+AEtyezSjOTvb9f0KCLS4CkBFalODV2fOrVI5bctOvPbwzpTWFLKnOVb+OTHDfx69jgSXEH5xyrOp/DdCfh6n01CXKhOWCIiEu027igk9e1bSK0wNUqiK+DigmeIP+5cDuvanIGdmpGSWGFKL3WnFRGpXRGiSFOBBYl5E5oClf+2ypzRx73IIZ2bMahTMwZ1bsaAjk1JT0nwb6BuuyIxQ0WI6pEq/vc651ixaRdZK7YwZ8UWZi/fzJL1O1iWdD6+kD1mDSZsjXDwIiLRab+KEInIXqpi7Ghhajt+3SODr3/ezAMf/oRz/mE/B7VqzCVNZvPr1XcTr267IiKRE6I7bfEbV/P0Z8t4dMshbNzhnxYlLTmeQzo346xBGRTP7kDSzkpTomtqFBGRWlACKlIXqhg7mjLiVm7t1weA7QXFZOfk8e2KLXy7cgvHrHiEeCp329319nhWtx1J15aNiAv1lbuumoqI7JONOwpp9M4tpFToTptQVsCpGx9nUc+TyezcnEM6N6Nb68b4dv8Pbj5BYzlFRPaRElCRulCLcT5pyQkccZB/zjcAN2FTyIdK3rWW4//1CamJcfRq14Q+HdLp3b4Jvdo3odu6t0mccW25b+511VREGqxqvpArLi3jh7XbmbdqK3NXbmXOis0s37SLZUmrIcR3e23dRv41ekDo42gsp4jIPlMCKlJX9rJsvlXRbbc0rT33nNafhbl5fLc6j6lZq9hVVArA50k3kWHlv7mnOB8381aspmPryqmI1CchutKWvHk10+bm8tyuISzMzaOwpAyAFo0SGdS5GecN7kTxN/vYnVZTo4iI7BMloCLRoopuuwknTuDsfhmcfYj/w1BZmePnTTtZvGY7HV4LfdXUbc3lpEmf0K1NGge0bETXoJ+mqYkhP6jpyqmIxKqtu4pIfucWkit0pY0vLWDIsgd5vu0QfnNYZwZ0bMqAjk3JaJaC7Z53s9kEdacVEYkgJaAi0aKWXbp8PuPAVo05sFVjmBn6qun2pDZ0ap7Kwtw83l6whrKggrzNUhN4l5toXRbqyunE6q+c6qqpiERKFf9vdhSW8P2abSzIyWN+zlbmr9pabVfa9raJVy4fWvVx1J1WRCSilICKRJO97dJVxVXT9FNv4/F+hwJQVFLGqi27+HnDTn7euJOfN+2k5fwNIR/Obc1hxH2fktEslYxmKYEf/+0D1swg9d3r9u2qqRJXEdkbIXppFL1+FXf973ue2JbJ7hnk2jZJpn/HdM45tBNFs9uTvHN1pYey2lSmVXdaEZGIUQIqEstq8c19Yrzvlyumuy0PfeV0W2JrMpqlkLNlF18t3cjOwFhTgM8TbybVV/mq6Y4Z45kVfwxt05Np0ySZFo0Sf6kUCfve3VdJq0j9UYu/5y07i/hu9TYW5OYx+vObaVFS/v9NoitkbOmzpA0/n74ZTejTPp3WTZJ/2aD5repKKyISA8w5V/NWEaZJtkXqWMWkEPwf1EY9sOdDoXOOvPxicrbkk7NlFye90hOj8v+LMmccUPjcLw8TZ7ROSw4kpEn8c/l5NC1eV3m/JhnYdQt/GYe1l/FVeV5KWiVCqppgu67FXBsZ4u+5LD6FRZm38XHiMBbmbmPh6jxytvyyflnyBfhC/L8Bgwlbqz+W/geIiESFqtpJXQEVaYhqceXUzGiamkjT1ET6dEiH90NfNS1r0oE3Rh/B2rwC1m0rYO22Atbl+X8vXrudJsXrQ8eQl0uPv79Dy8ZJtGycSIvGSbRolEjzxolcNW88acWVr7aWfXArvqo+TO5PYSV9aBWpEwXFpfjem0Bihb9nX0k+Tb+8k3uK2tOlRSr9OzblN4d1pk/7dPp0aILvP6H/36gyrYhI7FMCKtJQhWm8afwJtzCgY1PoWMV+k0J/kNyZ3JYL+3dh445CNu4oYt22Ahat3samnYX8LX5NyGIi5OXS/ea3aZaaSNPUBNJTEmiamkCz1ERu+OHvpIdIWovem0Bu+1NonBRPWnI8SfG+8lddI91FWMmuxLIq3r/OOdbkFbB47Ta+X7Od79ds4/s12/h5406WJOaG/Hvu4NvEggknkpacUHllFf9v1J1WRCT2KQEVkdrZ10qRVXyQTDtlIjf261lpc+ccblIGbMuptG5Hchsu7N+FrbuK2LqrmK35xSzfuIu5u7byz6J1IT/kxm9fzbH3fPzLoeOMtOSEPQnpU1tvolWIisDb/vd3Xtt+KKmJ8aQkxtEoKY6UhHhSE+NovXwabT75K76SfUhaI5XsRjJBVlLdMFRRGOg/H/7EE9sOZeuu4j2bZjRLoWe7JpzStx2Fc9uTkh+6OFDI5BNUmVZEpB6r1RhQMzsZuB+IAx53zt1ZYX0S8F/gEGATcI5zbnlg3Q3AGKAUGOece7em48Xc+BYRqd7eJij7MgZ0Up+QV1rzU9vzzgnvs72gZM/PjsJi/++CEh5bfkLIsWYVx7YG+zxxHBm+jZWWr6El56Q+TlK8j6QEH8nxcSQl+EiKjyM5wcc/fj6P5iWVx8NuT27Hm8PeDewXR3Lgd1K8j7YrptHpixt+SXYBF5/CzpP+Bf1GkxBnJPh81Rd+qs3zt6/77euxdu8biWS3jhLkcI0BNbOO+NvQNoADJjvn7q9q+7C1kTU8L3n5xfy0bjs/rtvBj+u2c+X8M2hVWrlL/Xpfayb1eY2e7dLo2a4J3dum0SQ4sdyf94iIiMSsqtrJGhNQM4sDfgROAHKA2cB5zrlFQdtcAfRzzo01s3OBM51z55hZL+AFYDDQHvgAONg5V1rxOMGUgIpIRJJWqDJxLWuSwdY/zGVnYQn5xaXsKiplV1EJuwpLGT61e8iCTA7juh4fUVhSRmFJGQXFpYHbpRQUl/HettPDluzmlLXkyKIH9tyP9xnxcUZCnI93uYL2VN5nva81V7f5L4nxPhLifCTEGfFxPuJ9RpzP+PtPo2kWomBUXmJbnsicRpzPR3ycf9vd+/z6sxE0LlhTaZ/8lPZ8dupHxMcZPjPifT7/foH9Wyx9k46fX4+vtHxhmg3H3kV+97OI8xlmEOfz7+8zI/n7V2n83h+xCsl46an3Y/1G4zMqF7Wqw+QnjAloO6Cdc+5bM0sD5gBnBLezwcLSRoZ4Xkp8ybzZ6W+8XnIEP63fzrpthXvWpSbGsdB37r4VBtp9PF3NFBFpUPanCNFgYIlzblnggV4ETgeCG8bTgQmB268AD5r/U8DpwIvOuULgZzNbEni8r/b1RESkgdjbMaph7iLsO/4WmjdKpHmjxMr7pIce12rpGdx37sCqj1XFeFjSO/DN74dTWPxLslpYUkphcRkdnt0U8qE6+DZx08ieFJWWUVLqKC4to7i0jKLSMtrNCb1Pq7INmMGuolKKS4spKvHvU+agpKyM9CoKRqUVruOBD5eEXHdhUujxukm71nDZM3NCPw/A54m34fNVLkxT/O6tDJveoop9/k5ahX2sJJ+1r93IkS/6pxkyg7hAwurzwUzfDXSwyl2smTkxahIg59waYE3g9nYz+x7oQPl2NrxmTiz/ngfiywo47OeHeKrlYI44qCUHt0nj4DaN6dY6jQ5NU/Ddv4+FgUDFgUREZI/aJKAdgOAWJwcYUtU2zrkSM8sDWgSWz6qwb4dQBzGzy4DLADp16lSb2EVEytuXD7n7krjua4GUapLd1mnJofepJtn9/dEHhN5nSdX7vHjZ4VXHV0WC7GuawfLrTqGszFFS5ihz/t+lpQ73aOjxuqVp7XnrsiMpDd6n1FFa5iguK6PDC1Un1v8a3Z8yB2VljlLn37eszNHh3ar3ue74g/3bBX5Ky/zjidt/E3of8irHHA3MrAswEPi6wvLwtpFVnH9728T0q48MvY8KA4mISBhETREi59xkYDL4uxd5HI6INCSRutoaqWQ3zAny7v18PiPRV+Fy5/G3hNwn4cQJ/ul7qlJNYv2rQVVcUZtV9T7XHN8t9D4/7MdVuwgzs8bAq8C1zrltwevC3kZW8/xXSYWBREQkDGqTgOZSfoKFjMCyUNvkmFk8kI6/GFFt9hURiT372qUwEsluJBPkMHd9DntiHSNX7cwsAX/y+Zxz7rU6P+C+Pi/qSisiIvupNkWI4vEXIRqOP3mcDZzvnPsuaJsrgb5BRYh+5ZwbbWa9gef5pQjRTKCbihCJiIiq4O55HAOeBjY7566taftIVcEVERHZH/tcBTew80jgPvzTsExxzt1uZhOBLOfcNDNLBp7BP25lM3BuUNGim4BLgBL83Yrerul4SkBFRCTahTEBPRL4DFgAlAUW3+icmxFqe7WRIiISC/anCi6BRnBGhWXjg24XAL+uYt/bgdv3KloREZEGwjn3OSFrCYuIiNQ/Pq8DEBERERERkYZBCaiIiIiIiIhEhBJQERERERERiQgloCIiIiIiIhIRtaqCG2lmtgFYEaaHawlsDNNjeSXWzyHW44fYP4dYjx9i/xxiPX6I/XMId/ydnXOtwvh4tRLmNhL0ukaDWD+HWI8fdA7RINbjh9g/h4i0k1GZgIaTmWWFo0y+l2L9HGI9foj9c4j1+CH2zyHW44fYP4dYj7+uxPrzEuvxQ+yfQ6zHDzqHaBDr8UPsn0Ok4lcXXBEREREREYkIJaAiIiIiIiISEQ0hAZ3sdQBhEOvnEOvxQ+yfQ6zHD7F/DrEeP8T+OcR6/HUl1p+XWI8fYv8cYj1+0DlEg1iPH2L/HCISf70fAyoiIiIiIiLRoSFcARUREREREZEooARUREREREREIiKmE1AzO9nMfjCzJWZ2fYj1SWb2UmD912bWJWjdDYHlP5jZSREN/JcYaor/j2a2yMyyzWymmXUOWldqZvMCP9MiG3m5GGs6h4vMbENQrJcGrbvQzH4K/FwY2cj3xFBT/JOCYv/RzLYGrfP8NTCzKWa23swWVrHezOyBwPllm9mgoHWeP/+BOGo6hwsCsS8wsy/NrH/QuuWB5fPMLCtyUZeLr6b4h5lZXtB7ZXzQumrff5FSi3P4S1D8CwPv/eaBddHwGnQ0s48C/y+/M7NrQmwT9X8L4RbrbWQgjphuJ2O9jQzEoXbS288pMd1GBuJQO6l2sjznXEz+AHHAUuAAIBGYD/SqsM0VwKOB2+cCLwVu9wpsnwR0DTxOXBTGfyyQGrh9+e74A/d3xMhrcBHwYIh9mwPLAr+bBW43i7b4K2x/NTAlyl6Do4FBwMIq1o8E3gYMOAz4Olqe/704h6G7YwNG7D6HwP3lQMsofw2GAW/t7/vPy3OosO0o4MMoew3aAYMCt9OAH0P8L4r6v4UwPycx3UbuxTlEbTtZy/gvIkrbyNqeQ4Xt1U5GPv6obiNreQ7DUDtZ1/FHVTsZy1dABwNLnHPLnHNFwIvA6RW2OR14OnD7FWC4mVlg+YvOuULn3M/AksDjRVKN8TvnPnLO7QrcnQVkRDjGmtTmNajKScD7zrnNzrktwPvAyXUUZ1X2Nv7zgBciElktOec+BTZXs8npwH+d3yygqZm1Izqef6Dmc3DOfRmIEaLw76AWr0FV9ufvJ6z28hyi8e9gjXPu28Dt7cD3QIcKm0X930KYxXobCbHfTsZ6GwlqJz1/DWK9jQS1k9Eg2trJWE5AOwCrgu7nUPmJ3LONc64EyANa1HLfura3MYzB/63EbslmlmVms8zsjDqIrzZqew5nBS7lv2JmHfdy37pU6xgC3bq6Ah8GLY6G16AmVZ1jNDz/+6Li34ED3jOzOWZ2mUcx1cbhZjbfzN42s96BZTH3GphZKv5G59WgxVH1Gpi/G+lA4OsKq+rb30JNYr2NZB/iiLZ2MtbbyL2KQ+1kVIjVNhLUTkZMNLST8fuzs0SGmf0GyASOCVrc2TmXa2YHAB+a2QLn3FJvIqzWdOAF51yhmf0B/7ftx3kc0744F3jFOVcatCxWXoN6wcyOxd+4Hhm0+MjAa9AaeN/MFge+pYwm3+J/r+wws5HAG0A3b0PaZ6OAL5xzwd8CR81rYGaN8Tf61zrntnkRg3gjhtvJ+tJGgtpJT8VwGwlqJyMmWtrJWL4Cmgt0DLqfEVgWchsziwfSgU213Leu1SoGMzseuAk4zTlXuHu5cy438HsZ8DH+bzIircZzcM5tCor7ceCQ2u4bAXsTw7lU6E4RJa9BTao6x2h4/mvNzPrhf/+c7pzbtHt50GuwHngdb7oJVss5t805tyNwewaQYGYtibHXIKC6vwNPXwMzS8DfqD7nnHstxCb14m9hL8R6G0lt44jidjLW28i9jUPtpEdiuY0EtZORElXtpPNwQOz+/OC/ersMf3eP3QOTe1fY5krKF1iYGrjdm/IFFpYR+SJEtYl/IP7B190qLG8GJAVutwR+woNB2bU8h3ZBt88EZgVuNwd+DpxLs8Dt5tEWf2C7HvgHkFu0vQaB43eh6oH9p1B+QPk30fL878U5dMI/Bm1oheWNgLSg218CJ0dh/G13v3fwNzorA69Hrd5/0XAOgfXp+Me/NIq21yDwfP4XuK+abWLibyGMz0lMt5F7cQ5R207WMv6obSNrew6B7dROehd/1LeRtTgHtZN1H3tUtZOevIBhfDJH4q/itBS4KbBsIv5vQQGSgZcDf5jfAAcE7XtTYL8fgBFRGv8HwDpgXuBnWmD5UGBB4A9xATAmil+DO4DvArF+BPQI2veSwGuzBLg4GuMP3J8A3Flhv6h4DfB/y7YGKMbfJ38MMBYYG1hvwEOB81sAZEbT81/Lc3gc2BL0d5AVWH5A4PmfH3iP3RSl8V8V9Dcwi6APCaHef9F4DoFtLsJfmCZ4v2h5DY7EP8YmO+h9MjLW/hbq4HmJ6TaylucQ1e1kLeKP6jayNucQuD8BtZNexR/VbWQtz0HtZN3HH1Xt5O5vG0RERERERETqVCyPARUREREREZEYogRUREREREREIkIJqIiIiIiIiESEElARERERERGJCCWgIiIiIiIiEhFKQEXqMTNramZXeB2HiIhINFI7KRJ5SkBF6remgBpWERGR0JqidlIkopSAitRvdwIHmtk8M7vb62BERESijNpJkQgz55zXMYhIHTGzLsBbzrk+XsciIiISbdROikSeroCKiIiIiIhIRCgBFRERERERkYhQAipSv20H0rwOQkREJEqpnRSJMCWgIvWYc24T8IWZLVRxBRERkfLUTopEnooQiYiIiIiISEToCqiIiIiIiIhEhBJQERERERERiQgloCIiIiIiIhIRSkBFREREREQkIpSAioiIiIiISEQoARUREREREZGIUAIqIiIiIiIiEfH/t+PNofZo2RwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] @@ -413,7 +424,7 @@ { "data": { "text/plain": [ - "{Variable(-0x5d1a09af22d01535, u, children=[], domain=[], auxiliary_domains={}): Multiplication(0x196c827dd508e63d, *, children=['-a', 'y[0:1]'], domain=[], auxiliary_domains={})}" + "{Variable(0x2d386efb3bed2d38, u, children=[], domain=[], auxiliary_domains={}): Multiplication(0x7851652896e183c2, *, children=['-a', 'y[0:1]'], domain=[], auxiliary_domains={})}" ] }, "execution_count": 13, @@ -493,7 +504,20 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.8" + "version": "3.8.6" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": true } }, "nbformat": 4, diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py index 9c50dbc7ef..8acb88e1dc 100644 --- a/examples/scripts/DFN.py +++ b/examples/scripts/DFN.py @@ -30,7 +30,7 @@ # solve model t_eval = np.linspace(0, 3600, 100) -solver = pybamm.CasadiSolver(mode="fast", atol=1e-6, rtol=1e-3) +solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-3) solution = solver.solve(model, t_eval) # plot diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 0bb685fefe..9a272781eb 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -56,6 +56,7 @@ def version(formatted=False): ABSOLUTE_PATH = root_dir() PARAMETER_PATH = [ + root_dir(), os.getcwd(), os.path.join(root_dir(), "pybamm", "input", "parameters"), ] diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py index f65b74bc43..984047ee53 100644 --- a/pybamm/models/base_model.py +++ b/pybamm/models/base_model.py @@ -109,6 +109,7 @@ def __init__(self, name="Unnamed model"): self.external_variables = [] self._parameters = None self._input_parameters = None + self._variables_casadi = {} # Default behaviour is to use the jacobian and simplify self.use_jacobian = True @@ -117,6 +118,7 @@ def __init__(self, name="Unnamed model"): # Model is not initially discretised self.is_discretised = False + self.y_slices = None # Default timescale is 1 second self.timescale = pybamm.Scalar(1) @@ -327,6 +329,14 @@ def new_empty_copy(self): new_model.convert_to_format = self.convert_to_format new_model.timescale = self.timescale new_model.length_scales = self.length_scales + + # Variables from discretisation + new_model.is_discretised = self.is_discretised + new_model.y_slices = self.y_slices + new_model.concatenated_rhs = self.concatenated_rhs + new_model.concatenated_algebraic = self.concatenated_algebraic + new_model.concatenated_initial_conditions = self.concatenated_initial_conditions + return new_model def new_copy(self): @@ -412,6 +422,31 @@ def set_initial_conditions_from(self, solution, inplace=True): "Variable must have type 'Variable' or 'Concatenation'" ) + # Also update the concatenated initial conditions if the model is already + # discretised + if model.is_discretised: + # Unpack slices for sorting + y_slices = {var.id: slce for var, slce in model.y_slices.items()} + slices = [] + for symbol in model.initial_conditions.keys(): + if isinstance(symbol, pybamm.Concatenation): + # must append the slice for the whole concatenation, so that + # equations get sorted correctly + slices.append( + slice( + y_slices[symbol.children[0].id][0].start, + y_slices[symbol.children[-1].id][0].stop, + ) + ) + else: + slices.append(y_slices[symbol.id][0]) + equations = list(model.initial_conditions.values()) + # sort equations according to slices + sorted_equations = [eq for _, eq in sorted(zip(slices, equations))] + model.concatenated_initial_conditions = pybamm.NumpyConcatenation( + *sorted_equations + ) + return model def check_and_combine_dict(self, dict1, dict2): @@ -888,10 +923,7 @@ def default_spatial_methods(self): @property def default_solver(self): "Return default solver based on whether model is ODE model or DAE model" - if len(self.algebraic) == 0: - return pybamm.ScipySolver() - else: - return pybamm.CasadiSolver(mode="safe") + return pybamm.CasadiSolver(mode="safe") # helper functions for finding symbols diff --git a/pybamm/simulation.py b/pybamm/simulation.py index e33fa0b76a..b3abee7022 100644 --- a/pybamm/simulation.py +++ b/pybamm/simulation.py @@ -517,35 +517,6 @@ def step(self, dt, solver=None, npts=2, save=True, **kwargs): return self.solution - def get_variable_array(self, *variables): - """ - A helper function to easily obtain a dictionary of arrays of values - for a list of variables at the latest timestep. - - Parameters - ---------- - variable: str - The name of the variable/variables you wish to obtain the arrays for. - - Returns - ------- - variable_arrays: dict - A dictionary of the variable names and their corresponding - arrays. - """ - - variable_arrays = [ - self.built_model.variables[var].evaluate( - self.solution.t[-1], self.solution.y[:, -1] - ) - for var in variables - ] - - if len(variable_arrays) == 1: - return variable_arrays[0] - else: - return tuple(variable_arrays) - def plot(self, output_variables=None, quick_plot_vars=None, **kwargs): """ A method to quickly plot the outputs of the simulation. Creates a @@ -695,6 +666,8 @@ def save(self, filename): and self._solver.integrator_specs != {} ): self._solver.integrator_specs = {} + if self.solution is not None: + self.solution.clear_casadi_attributes() with open(filename, "wb") as f: pickle.dump(self, f, pickle.HIGHEST_PROTOCOL) diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 1f43263dad..270a179f4e 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -920,7 +920,11 @@ def get_termination_reason(self, solution, events): # causes an error later in ProcessedVariable) if solution.t_event - solution._t[-1] > self.atol: solution._t = np.concatenate((solution._t, solution.t_event)) - solution._y = np.concatenate((solution._y, solution.y_event), axis=1) + if isinstance(solution.y, casadi.DM): + solution._y = casadi.horzcat(solution.y, solution.y_event) + else: + solution._y = np.hstack((solution._y, solution.y_event)) + for name, inp in solution.inputs.items(): solution._inputs[name] = np.c_[inp, inp[:, -1]] @@ -965,7 +969,6 @@ def __init__(self, function, name, model): self.timescale = self.model.timescale_eval def __call__(self, t, y, inputs): - y = y.reshape(-1, 1) if self.name in ["RHS", "algebraic", "residuals"]: pybamm.logger.debug( "Evaluating {} for {} at t={}".format( diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index d67ecf2cea..54ef2be791 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -135,8 +135,6 @@ def _integrate(self, model, t_eval, inputs=None): return solution elif self.mode in ["safe", "safe without grid"]: y0 = model.y0 - if isinstance(y0, casadi.DM): - y0 = y0.full().flatten() # Step-and-check t = t_eval[0] t_f = t_eval[-1] @@ -180,7 +178,7 @@ def _integrate(self, model, t_eval, inputs=None): # to avoid having to create several times self.create_integrator(model, inputs) # Initialize solution - solution = pybamm.Solution(np.array([t]), y0[:, np.newaxis]) + solution = pybamm.Solution(np.array([t]), y0) solution.solve_time = 0 solution.integration_time = 0 else: @@ -460,7 +458,7 @@ def _run_integrator(self, model, y0, inputs, t_eval): x0=y0_diff, z0=y0_alg, p=inputs, **self.extra_options_call ) integration_time = timer.time() - y_sol = np.concatenate([sol["xf"].full(), sol["zf"].full()]) + y_sol = casadi.vertcat(sol["xf"], sol["zf"]) sol = pybamm.Solution(t_eval, y_sol) sol.integration_time = integration_time return sol diff --git a/pybamm/solvers/processed_symbolic_variable.py b/pybamm/solvers/processed_symbolic_variable.py index f3a18e9ee7..443e39ec0a 100644 --- a/pybamm/solvers/processed_symbolic_variable.py +++ b/pybamm/solvers/processed_symbolic_variable.py @@ -38,7 +38,6 @@ def __init__(self, base_variable, solution): symbolic_inputs_dict[key] = value all_inputs_as_MX = casadi.vertcat(*[p for p in all_inputs_as_MX_dict.values()]) - all_inputs = casadi.vertcat(*[p for p in solution.inputs.values()]) # The symbolic_inputs dictionary will be used for sensitivity symbolic_inputs = casadi.vertcat(*[p for p in symbolic_inputs_dict.values()]) var = base_variable.to_casadi(t_MX, y_MX, inputs=all_inputs_as_MX_dict) @@ -52,10 +51,13 @@ def __init__(self, base_variable, solution): self.mesh = base_variable.mesh self.symbolic_inputs_dict = symbolic_inputs_dict self.symbolic_inputs_total_shape = symbolic_inputs.shape[0] - self.inputs = all_inputs + self.inputs = solution.inputs self.domain = base_variable.domain - self.base_eval = self.base_variable(solution.t[0], solution.y[:, 0], all_inputs) + self.inputs = casadi.vertcat(*[p for p in solution.inputs.values()]) + self.base_eval = self.base_variable( + solution.t[0], solution.y[:, 0], self.inputs + ) if ( isinstance(self.base_eval, numbers.Number) diff --git a/pybamm/solvers/processed_variable.py b/pybamm/solvers/processed_variable.py index 966c5026a6..dfc68377dd 100644 --- a/pybamm/solvers/processed_variable.py +++ b/pybamm/solvers/processed_variable.py @@ -1,6 +1,7 @@ # # Processed Variable class # +import casadi import numbers import numpy as np import pybamm @@ -39,24 +40,25 @@ class ProcessedVariable(object): variable. Note that this can be any kind of node in the expression tree, not just a :class:`pybamm.Variable`. When evaluated, returns an array of size (m,n) + base_variable_casadi : :class:`casadi.Function` + A casadi function. When evaluated, returns the same thing as + `base_Variable.evaluate` (but more efficiently). solution : :class:`pybamm.Solution` The solution object to be used to create the processed variables - known_evals : dict - Dictionary of known evaluations, to be used to speed up finding the solution warn : bool, optional Whether to raise warnings when trying to evaluate time and length scales. Default is True. """ - def __init__(self, base_variable, solution, known_evals=None, warn=True): + def __init__(self, base_variable, base_variable_casadi, solution, warn=True): self.base_variable = base_variable + self.base_variable_casadi = base_variable_casadi self.t_sol = solution.t self.u_sol = solution.y self.mesh = base_variable.mesh self.inputs = solution.inputs self.domain = base_variable.domain self.auxiliary_domains = base_variable.auxiliary_domains - self.known_evals = known_evals self.warn = warn # Set timescale @@ -68,19 +70,10 @@ def __init__(self, base_variable, solution, known_evals=None, warn=True): self.length_scales = solution.length_scales_eval # Evaluate base variable at initial time - if self.known_evals: - self.base_eval, self.known_evals[solution.t[0]] = base_variable.evaluate( - solution.t[0], - solution.y[:, 0], - inputs={name: inp[:, 0] for name, inp in solution.inputs.items()}, - known_evals=self.known_evals[solution.t[0]], - ) - else: - self.base_eval = base_variable.evaluate( - solution.t[0], - solution.y[:, 0], - inputs={name: inp[:, 0] for name, inp in solution.inputs.items()}, - ) + inputs = casadi.vertcat(*[inp[:, 0] for inp in self.inputs.values()]) + self.base_eval = self.base_variable_casadi( + solution.t[0], solution.y[:, 0], inputs + ).full() # handle 2D (in space) finite element variables differently if ( @@ -128,13 +121,8 @@ def initialise_0D(self): for idx in range(len(self.t_sol)): t = self.t_sol[idx] u = self.u_sol[:, idx] - inputs = {name: inp[:, idx] for name, inp in self.inputs.items()} - if self.known_evals: - entries[idx], self.known_evals[t] = self.base_variable.evaluate( - t, u, inputs=inputs, known_evals=self.known_evals[t] - ) - else: - entries[idx] = self.base_variable.evaluate(t, u, inputs=inputs) + inputs = casadi.vertcat(*[inp[:, idx] for inp in self.inputs.values()]) + entries[idx] = self.base_variable_casadi(t, u, inputs).full()[0, 0] # set up interpolation if len(self.t_sol) == 1: @@ -164,15 +152,8 @@ def initialise_1D(self, fixed_t=False): for idx in range(len(self.t_sol)): t = self.t_sol[idx] u = self.u_sol[:, idx] - inputs = {name: inp[:, idx] for name, inp in self.inputs.items()} - if self.known_evals: - eval_and_known_evals = self.base_variable.evaluate( - t, u, inputs=inputs, known_evals=self.known_evals[t] - ) - entries[:, idx] = eval_and_known_evals[0][:, 0] - self.known_evals[t] = eval_and_known_evals[1] - else: - entries[:, idx] = self.base_variable.evaluate(t, u, inputs=inputs)[:, 0] + inputs = casadi.vertcat(*[inp[:, idx] for inp in self.inputs.values()]) + entries[:, idx] = self.base_variable_casadi(t, u, inputs).full()[:, 0] # Get node and edge values nodes = self.mesh.nodes @@ -274,23 +255,12 @@ def initialise_2D(self): for idx in range(len(self.t_sol)): t = self.t_sol[idx] u = self.u_sol[:, idx] - inputs = {name: inp[:, idx] for name, inp in self.inputs.items()} - if self.known_evals: - eval_and_known_evals = self.base_variable.evaluate( - t, u, inputs=inputs, known_evals=self.known_evals[t] - ) - entries[:, :, idx] = np.reshape( - eval_and_known_evals[0], - [first_dim_size, second_dim_size], - order="F", - ) - self.known_evals[t] = eval_and_known_evals[1] - else: - entries[:, :, idx] = np.reshape( - self.base_variable.evaluate(t, u, inputs=inputs), - [first_dim_size, second_dim_size], - order="F", - ) + inputs = casadi.vertcat(*[inp[:, idx] for inp in self.inputs.values()]) + entries[:, :, idx] = np.reshape( + self.base_variable_casadi(t, u, inputs).full(), + [first_dim_size, second_dim_size], + order="F", + ) # add points outside first dimension domain for extrapolation to # boundaries @@ -423,22 +393,13 @@ def initialise_2D_scikit_fem(self): for idx in range(len(self.t_sol)): t = self.t_sol[idx] u = self.u_sol[:, idx] - inputs = {name: inp[:, idx] for name, inp in self.inputs.items()} + inputs = casadi.vertcat(*[inp[:, idx] for inp in self.inputs.values()]) - if self.known_evals: - eval_and_known_evals = self.base_variable.evaluate( - t, u, inputs=inputs, known_evals=self.known_evals[t] - ) - entries[:, :, idx] = np.reshape( - eval_and_known_evals[0], [len_y, len_z], order="F" - ) - self.known_evals[t] = eval_and_known_evals[1] - else: - entries[:, :, idx] = np.reshape( - self.base_variable.evaluate(t, u, inputs=inputs), - [len_y, len_z], - order="F", - ) + entries[:, :, idx] = np.reshape( + self.base_variable_casadi(t, u, inputs).full(), + [len_y, len_z], + order="F", + ) # assign attributes for reference self.entries = entries diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py index 3ddbf5514e..f92c014dbd 100644 --- a/pybamm/solvers/solution.py +++ b/pybamm/solvers/solution.py @@ -8,7 +8,6 @@ import pickle import pybamm import pandas as pd -from collections import defaultdict from scipy.io import savemat @@ -40,10 +39,8 @@ class _BaseSolution(object): def __init__( self, t, y, t_event=None, y_event=None, termination="final time", copy_this=None ): - self._t = t - if isinstance(y, casadi.DM): - y = y.full() - self._y = y + self.t = t + self.y = y self._t_event = t_event self._y_event = y_event self._termination = termination @@ -67,21 +64,24 @@ def __init__( self._variables = pybamm.FuzzyDict() self.data = pybamm.FuzzyDict() - # initialize empty known evals - self._known_evals = defaultdict(dict) - for time in t: - self._known_evals[time] = {} - @property def t(self): "Times at which the solution is evaluated" return self._t + @t.setter + def t(self, t): + self._t = t + @property def y(self): "Values of the solution" return self._y + @y.setter + def y(self, y): + self._y = y + @property def model(self): "Model used for solution" @@ -118,7 +118,13 @@ def inputs(self, inputs): # If there are symbolic inputs, just store them as given if any(isinstance(v, casadi.MX) for v in inputs.values()): self.has_symbolic_inputs = True - self._inputs = inputs + self._inputs = {} + for name, inp in inputs.items(): + if isinstance(inp, numbers.Number): + self._inputs[name] = casadi.DM([inp]) + else: + self._inputs[name] = inp + # Otherwise, make them the same size as the time vector else: self.has_symbolic_inputs = False @@ -181,13 +187,35 @@ def update(self, variables): # Otherwise a standard ProcessedVariable is ok else: - var = pybamm.ProcessedVariable( - self.model.variables[key], self, self._known_evals - ) + var_pybamm = self.model.variables[key] + + if key in self.model._variables_casadi: + var_casadi = self.model._variables_casadi[key] + else: + self._t_MX = casadi.MX.sym("t") + self._y_MX = casadi.MX.sym("y", self.y.shape[0]) + self._symbolic_inputs_dict = { + key: casadi.MX.sym("input", value.shape[0]) + for key, value in self._inputs.items() + } + self._symbolic_inputs = casadi.vertcat( + *[p for p in self._symbolic_inputs_dict.values()] + ) - # Update known_evals in order to process any other variables faster - for t in var.known_evals: - self._known_evals[t].update(var.known_evals[t]) + # Convert variable to casadi + # Make all inputs symbolic first for converting to casadi + var_sym = var_pybamm.to_casadi( + self._t_MX, self._y_MX, inputs=self._symbolic_inputs_dict + ) + + var_casadi = casadi.Function( + "variable", + [self._t_MX, self._y_MX, self._symbolic_inputs], + [var_sym], + ) + self.model._variables_casadi[key] = var_casadi + + var = pybamm.ProcessedVariable(var_pybamm, var_casadi, self) # Save variable and data self._variables[key] = var @@ -234,10 +262,20 @@ def plot(self, output_variables=None, **kwargs): """ return pybamm.dynamic_plot(self, output_variables=output_variables, **kwargs) + def clear_casadi_attributes(self): + "Remove casadi objects for pickling, will be computed again automatically" + self._t_MX = None + self._y_MX = None + self._symbolic_inputs = None + self._symbolic_inputs_dict = None + def save(self, filename): """Save the whole solution using pickle""" # No warning here if len(self.data)==0 as solution can be loaded # and used to process new variables + + self.clear_casadi_attributes() + # Pickle with open(filename, "wb") as f: pickle.dump(self, f, pickle.HIGHEST_PROTOCOL) @@ -392,7 +430,10 @@ def append(self, solution, start_index=1, create_sub_solutions=False): # Update t, y and inputs self._t = np.concatenate((self._t, solution.t[start_index:])) - self._y = np.concatenate((self._y, solution.y[:, start_index:]), axis=1) + if isinstance(self.y, casadi.DM) and isinstance(solution.y, casadi.DM): + self._y = casadi.horzcat(self.y, solution.y[:, start_index:]) + else: + self._y = np.hstack((self._y, solution.y[:, start_index:])) for name, inp in self.inputs.items(): solution_inp = solution.inputs[name] self.inputs[name] = np.c_[inp, solution_inp[:, start_index:]] @@ -404,9 +445,6 @@ def append(self, solution, start_index=1, create_sub_solutions=False): self._t_event = solution._t_event self._y_event = solution._y_event - # Update known_evals - for t, evals in solution._known_evals.items(): - self._known_evals[t].update(evals) # Recompute existing variables for var in self._variables.keys(): self.update(var) diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py index 08798ca215..577b2bbd15 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_external/test_external_current_collector.py @@ -7,7 +7,6 @@ class TestExternalCC(unittest.TestCase): - @unittest.skipIf(not pybamm.have_idaklu(), "idaklu solver is not installed") def test_2p1d(self): model_options = { "current collector": "potential pair", @@ -25,8 +24,7 @@ def test_2p1d(self): pybamm.standard_spatial_vars.y: yz_pts, pybamm.standard_spatial_vars.z: yz_pts, } - solver = pybamm.IDAKLUSolver() - sim = pybamm.Simulation(model, var_pts=var_pts, solver=solver) + sim = pybamm.Simulation(model, var_pts=var_pts) # Simulate 100 seconds t_eval = np.linspace(0, 100, 3) @@ -45,7 +43,7 @@ def test_2p1d(self): sim.step(dt, external_variables=external_variables) # obtain phi_s_n from the pybamm solution at the current time - phi_s_p = sim.get_variable_array("Positive current collector potential") + phi_s_p = sim.solution["Positive current collector potential"].data[:, -1] self.assertTrue(phi_s_p.shape, (yz_pts ** 2, 1)) diff --git a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py index 3bc73e97ac..29d388a3c4 100644 --- a/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py +++ b/tests/integration/test_models/test_submodels/test_external_circuit/test_function_control.py @@ -23,7 +23,6 @@ def constant_current(variables): # First model: 1A charge params[0]["Current function [A]"] = -1 - params[1]["Current function [A]"] = -1 # set parameters and discretise models for i, model in enumerate(models): diff --git a/tests/integration/test_solvers/test_external_variables.py b/tests/integration/test_solvers/test_external_variables.py index 29ef229f6b..8ad4759c50 100644 --- a/tests/integration/test_solvers/test_external_variables.py +++ b/tests/integration/test_solvers/test_external_variables.py @@ -44,12 +44,8 @@ def test_external_variables_SPMe(self): external_variables = {"Volume-averaged cell temperature": T_av} T_av += 1 sim.step(dt, external_variables=external_variables) - var = "Terminal voltage [V]" - t = sim.solution.t[-1] - y = sim.solution.y[:, -1] - inputs = external_variables - sim.built_model.variables[var].evaluate(t, y, inputs=inputs) - sim.solution[var](t) + V = sim.solution["Terminal voltage [V]"].data + np.testing.assert_array_less(np.diff(V), 0) # test generate with external variable sim.built_model.generate("test.c", ["Volume-averaged cell temperature"]) os.remove("test.c") diff --git a/tests/integration/test_solvers/test_solution.py b/tests/integration/test_solvers/test_solution.py index aa3c3caf62..fc92bab4f3 100644 --- a/tests/integration/test_solvers/test_solution.py +++ b/tests/integration/test_solvers/test_solution.py @@ -43,15 +43,6 @@ def test_append(self): step_solution.update("Terminal voltage") old_t = t - # Step solution should have been updated as we go along so be quicker to - # calculate - timer = pybamm.Timer() - step_solution.update("Terminal voltage") - step_sol_time = timer.time() - timer.reset() - solution.update("Terminal voltage") - sol_time = timer.time() - self.assertLess(step_sol_time.value, sol_time.value) # Check both give the same answer np.testing.assert_array_almost_equal( solution["Terminal voltage"](solution.t[:-1] * model.timescale_eval), diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index fc3d53be3a..4c3a523a53 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/tests/unit/test_expression_tree/test_unary_operators.py @@ -16,11 +16,10 @@ def test_unary_operator(self): self.assertEqual(un.domain, a.domain) # with number - absval = pybamm.AbsoluteValue(-10) - self.assertEqual(absval.evaluate(), 10) - - log = pybamm.log(10) - self.assertEqual(log.evaluate(), np.log(10)) + a = pybamm.InputParameter("a") + absval = pybamm.AbsoluteValue(-a) + self.assertEqual(absval.evaluate(inputs={"a": 10}), 10) + self.assertEqual(absval.evaluate(inputs={"a": 10}, known_evals={})[0], 10) def test_negation(self): a = pybamm.Symbol("a") diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py index b023f1e439..b9c780be60 100644 --- a/tests/unit/test_expression_tree/test_variable.py +++ b/tests/unit/test_expression_tree/test_variable.py @@ -91,6 +91,9 @@ def test_external_variable_vector(self): a_test = 2 * np.ones((10, 1)) np.testing.assert_array_equal(a.evaluate(inputs={"a": a_test}), a_test) + np.testing.assert_array_equal( + a.evaluate(inputs={"a": a_test.flatten()}), a_test + ) np.testing.assert_array_equal(a.evaluate(inputs={"a": 2}), a_test) diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py index 9763b1dc87..3635e8f8ad 100644 --- a/tests/unit/test_models/test_base_model.py +++ b/tests/unit/test_models/test_base_model.py @@ -660,7 +660,7 @@ def test_set_initial_conditions(self): self.assertEqual(model.initial_conditions[var_2D].value, 1) self.assertEqual(model.initial_conditions[var_concat].value, 1) - # Update initial conditions + # Discretise var = pybamm.standard_spatial_vars geometry = { "negative electrode": {var.x_n: {"min": 0, "max": 1}}, @@ -685,15 +685,20 @@ def test_set_initial_conditions(self): y = np.tile(3 * t, (1 + 30 + 50, 1)) sol = pybamm.Solution(t, y) sol.model = model_disc + sol.inputs = {} # Update out-of-place first, since otherwise we'll have already modified the # model new_model = model.set_initial_conditions_from(sol, inplace=False) # Make sure original model is unchanged - self.assertEqual(model.initial_conditions[var_scalar].value, 1) - self.assertEqual(model.initial_conditions[var_1D].value, 1) - self.assertEqual(model.initial_conditions[var_2D].value, 1) - self.assertEqual(model.initial_conditions[var_concat].value, 1) + np.testing.assert_array_equal( + model.initial_conditions[var_scalar].evaluate(), 1 + ) + np.testing.assert_array_equal(model.initial_conditions[var_1D].evaluate(), 1) + np.testing.assert_array_equal(model.initial_conditions[var_2D].evaluate(), 1) + np.testing.assert_array_equal( + model.initial_conditions[var_concat].evaluate(), 1 + ) # Now update inplace model.set_initial_conditions_from(sol) @@ -719,6 +724,43 @@ def test_set_initial_conditions(self): self.assertEqual(mdl.initial_conditions[var_concat].shape, (20, 1)) np.testing.assert_array_equal(mdl.initial_conditions[var_concat].entries, 3) + # Test updating a discretised model (out-of-place) + new_model_disc = model_disc.set_initial_conditions_from(sol, inplace=False) + + # Test new initial conditions + var_scalar = list(new_model_disc.initial_conditions.keys())[0] + self.assertIsInstance( + new_model_disc.initial_conditions[var_scalar], pybamm.Vector + ) + self.assertEqual(new_model_disc.initial_conditions[var_scalar].entries, 3) + + var_1D = list(new_model_disc.initial_conditions.keys())[1] + self.assertIsInstance(new_model_disc.initial_conditions[var_1D], pybamm.Vector) + self.assertEqual(new_model_disc.initial_conditions[var_1D].shape, (10, 1)) + np.testing.assert_array_equal( + new_model_disc.initial_conditions[var_1D].entries, 3 + ) + + var_2D = list(new_model_disc.initial_conditions.keys())[2] + self.assertIsInstance(new_model_disc.initial_conditions[var_2D], pybamm.Vector) + self.assertEqual(new_model_disc.initial_conditions[var_2D].shape, (50, 1)) + np.testing.assert_array_equal( + new_model_disc.initial_conditions[var_2D].entries, 3 + ) + + var_concat = list(new_model_disc.initial_conditions.keys())[3] + self.assertIsInstance( + new_model_disc.initial_conditions[var_concat], pybamm.Vector + ) + self.assertEqual(new_model_disc.initial_conditions[var_concat].shape, (20, 1)) + np.testing.assert_array_equal( + new_model_disc.initial_conditions[var_concat].entries, 3 + ) + + np.testing.assert_array_equal( + new_model_disc.concatenated_initial_conditions.evaluate(), 3 + ) + # Test updating a new model with a different model new_model = pybamm.BaseModel() new_var_scalar = pybamm.Variable("var_scalar") @@ -818,6 +860,44 @@ def test_set_initial_conditions(self): new_model.initial_conditions[var_concat].entries, 5 ) + # Test updating a discretised model (out-of-place) + model_disc = disc.process_model(model, inplace=False) + new_model_disc = model_disc.set_initial_conditions_from(sol_dict, inplace=False) + + # Test new initial conditions + var_scalar = list(new_model_disc.initial_conditions.keys())[0] + self.assertIsInstance( + new_model_disc.initial_conditions[var_scalar], pybamm.Vector + ) + self.assertEqual(new_model_disc.initial_conditions[var_scalar].entries, 5) + + var_1D = list(new_model_disc.initial_conditions.keys())[1] + self.assertIsInstance(new_model_disc.initial_conditions[var_1D], pybamm.Vector) + self.assertEqual(new_model_disc.initial_conditions[var_1D].shape, (10, 1)) + np.testing.assert_array_equal( + new_model_disc.initial_conditions[var_1D].entries, 5 + ) + + var_2D = list(new_model_disc.initial_conditions.keys())[2] + self.assertIsInstance(new_model_disc.initial_conditions[var_2D], pybamm.Vector) + self.assertEqual(new_model_disc.initial_conditions[var_2D].shape, (50, 1)) + np.testing.assert_array_equal( + new_model_disc.initial_conditions[var_2D].entries, 5 + ) + + var_concat = list(new_model_disc.initial_conditions.keys())[3] + self.assertIsInstance( + new_model_disc.initial_conditions[var_concat], pybamm.Vector + ) + self.assertEqual(new_model_disc.initial_conditions[var_concat].shape, (20, 1)) + np.testing.assert_array_equal( + new_model_disc.initial_conditions[var_concat].entries, 5 + ) + + np.testing.assert_array_equal( + new_model_disc.concatenated_initial_conditions.evaluate(), 5 + ) + def test_set_initial_condition_errors(self): model = pybamm.BaseModel() var = pybamm.Scalar(1) @@ -855,16 +935,12 @@ def test_set_initial_condition_errors(self): class TestStandardBatteryBaseModel(unittest.TestCase): def test_default_solver(self): model = pybamm.BaseBatteryModel() - self.assertIsInstance( - model.default_solver, (pybamm.ScipySolver, pybamm.ScikitsOdeSolver) - ) + self.assertIsInstance(model.default_solver, pybamm.CasadiSolver) # check that default_solver gives you a new solver, not an internal object solver = model.default_solver solver = pybamm.BaseModel() - self.assertIsInstance( - model.default_solver, (pybamm.ScipySolver, pybamm.ScikitsOdeSolver) - ) + self.assertIsInstance(model.default_solver, pybamm.CasadiSolver) self.assertIsInstance(solver, pybamm.BaseModel) # check that adding algebraic variables gives DAE solver diff --git a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_full.py b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_full.py index dba32e8355..6dae15249a 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_full.py +++ b/tests/unit/test_models/test_full_battery_models/test_lead_acid/test_full.py @@ -47,7 +47,7 @@ def test_well_posed_surface_form_differential(self): options = {"side reactions": ["oxygen"], "surface form": "differential"} model = pybamm.lead_acid.Full(options) model.check_well_posedness() - self.assertIsInstance(model.default_solver, pybamm.ScipySolver) + self.assertIsInstance(model.default_solver, pybamm.CasadiSolver) def test_well_posed_surface_form_algebraic(self): options = {"side reactions": ["oxygen"], "surface form": "algebraic"} diff --git a/tests/unit/test_simulation.py b/tests/unit/test_simulation.py index 35979a8f08..0d5152d9df 100644 --- a/tests/unit/test_simulation.py +++ b/tests/unit/test_simulation.py @@ -15,7 +15,7 @@ def test_simple_model(self): param = pybamm.ParameterValues({"a": 1}) sim = pybamm.Simulation(model, parameter_values=param) sol = sim.solve([0, 1]) - np.testing.assert_array_almost_equal(sol.y[0], np.exp(-sol.t), decimal=5) + np.testing.assert_array_almost_equal(sol.y.full()[0], np.exp(-sol.t), decimal=5) def test_basic_ops(self): @@ -140,22 +140,6 @@ def test_set_crate(self): self.assertEqual(sim.parameter_values["Current function [A]"], 2 * current_1C) self.assertEqual(sim.C_rate, 2) - def test_get_variable_array(self): - - sim = pybamm.Simulation(pybamm.lithium_ion.SPM()) - sim.solve([0, 600]) - - phi_s_n = sim.get_variable_array("Negative electrode potential") - - self.assertIsInstance(phi_s_n, np.ndarray) - - c_s_n_surf, c_e = sim.get_variable_array( - "Negative particle surface concentration", "Electrolyte concentration" - ) - - self.assertIsInstance(c_s_n_surf, np.ndarray) - self.assertIsInstance(c_e, np.ndarray) - def test_set_external_variable(self): model_options = { "thermal": "lumped", @@ -193,18 +177,18 @@ def test_step(self): sim.step(dt) # 1 step stores first two points tau = sim.model.timescale.evaluate() self.assertEqual(sim.solution.t.size, 2) - self.assertEqual(sim.solution.y[0, :].size, 2) + self.assertEqual(sim.solution.y.full()[0, :].size, 2) self.assertEqual(sim.solution.t[0], 0) self.assertEqual(sim.solution.t[1], dt / tau) sim.step(dt) # automatically append the next step self.assertEqual(sim.solution.t.size, 3) - self.assertEqual(sim.solution.y[0, :].size, 3) + self.assertEqual(sim.solution.y.full()[0, :].size, 3) self.assertEqual(sim.solution.t[0], 0) self.assertEqual(sim.solution.t[1], dt / tau) self.assertEqual(sim.solution.t[2], 2 * dt / tau) sim.step(dt, save=False) # now only store the two end step points self.assertEqual(sim.solution.t.size, 2) - self.assertEqual(sim.solution.y[0, :].size, 2) + self.assertEqual(sim.solution.y.full()[0, :].size, 2) self.assertEqual(sim.solution.t[0], 2 * dt / tau) self.assertEqual(sim.solution.t[1], 3 * dt / tau) @@ -227,7 +211,7 @@ def test_step_with_inputs(self): ) # 1 step stores first two points tau = sim.model.timescale.evaluate() self.assertEqual(sim.solution.t.size, 2) - self.assertEqual(sim.solution.y[0, :].size, 2) + self.assertEqual(sim.solution.y.full()[0, :].size, 2) self.assertEqual(sim.solution.t[0], 0) self.assertEqual(sim.solution.t[1], dt / tau) np.testing.assert_array_equal(sim.solution.inputs["Current function [A]"], 1) @@ -235,7 +219,7 @@ def test_step_with_inputs(self): dt, inputs={"Current function [A]": 2} ) # automatically append the next step self.assertEqual(sim.solution.t.size, 3) - self.assertEqual(sim.solution.y[0, :].size, 3) + self.assertEqual(sim.solution.y.full()[0, :].size, 3) self.assertEqual(sim.solution.t[0], 0) self.assertEqual(sim.solution.t[1], dt / tau) self.assertEqual(sim.solution.t[2], 2 * dt / tau) diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py index 1fb2c81c13..d7055f14d3 100644 --- a/tests/unit/test_solvers/test_base_solver.py +++ b/tests/unit/test_solvers/test_base_solver.py @@ -171,7 +171,7 @@ def algebraic_eval(self, t, y, inputs): np.testing.assert_array_almost_equal(init_cond, vec) # with casadi init_cond = solver_with_casadi.calculate_consistent_state(model) - np.testing.assert_array_almost_equal(init_cond, vec) + np.testing.assert_array_almost_equal(init_cond.full().flatten(), vec) # With jacobian def jac_dense(t, y, inputs): diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py index 862ea298ca..9d60a8d164 100644 --- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py +++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py @@ -96,14 +96,8 @@ def test_model_solver_with_time(self): sol = np.vstack((3 * t_eval, 6 * t_eval)) np.testing.assert_array_almost_equal(solution.y, sol) - np.testing.assert_array_almost_equal( - model.variables["var1"].evaluate(t=t_eval, y=solution.y).flatten(), - sol[0, :], - ) - np.testing.assert_array_almost_equal( - model.variables["var2"].evaluate(t=t_eval, y=solution.y).flatten(), - sol[1, :], - ) + np.testing.assert_array_almost_equal(solution["var1"].data.flatten(), sol[0, :]) + np.testing.assert_array_almost_equal(solution["var2"].data.flatten(), sol[1, :]) def test_model_solver_with_time_not_changing(self): # Create model @@ -137,9 +131,7 @@ def test_model_solver_with_bounds(self): # Solve solver = pybamm.CasadiAlgebraicSolver(tol=1e-12) solution = solver.solve(model) - np.testing.assert_array_almost_equal( - model.variables["var1"].evaluate(t=None, y=solution.y), 3 * np.pi / 2 - ) + np.testing.assert_array_almost_equal(solution["var1"].data, 3 * np.pi / 2) def test_solve_with_input(self): # Simple system: a single algebraic equation diff --git a/tests/unit/test_solvers/test_casadi_solver.py b/tests/unit/test_solvers/test_casadi_solver.py index f620efbf13..f7b3cf21c2 100644 --- a/tests/unit/test_solvers/test_casadi_solver.py +++ b/tests/unit/test_solvers/test_casadi_solver.py @@ -31,7 +31,7 @@ def test_model_solver(self): solution = solver.solve(model_disc, t_eval) np.testing.assert_array_equal(solution.t, t_eval) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) # Safe mode (enforce events that won't be triggered) @@ -41,7 +41,7 @@ def test_model_solver(self): solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) # Safe mode, without grid (enforce events that won't be triggered) @@ -49,7 +49,7 @@ def test_model_solver(self): solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) def test_model_solver_python(self): @@ -71,7 +71,7 @@ def test_model_solver_python(self): solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) pybamm.set_logging_level("WARNING") @@ -121,13 +121,13 @@ def test_model_solver_events(self): solver = pybamm.CasadiSolver(mode="safe", rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval) - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-10) + np.testing.assert_array_less(solution.y.full()[0], 1.5) + np.testing.assert_array_less(solution.y.full()[-1], 2.5 + 1e-10) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) np.testing.assert_array_almost_equal( - solution.y[-1], 2 * np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[-1], 2 * np.exp(0.1 * solution.t), decimal=5 ) # Solve using "safe" mode with debug off @@ -135,15 +135,15 @@ def test_model_solver_events(self): solver = pybamm.CasadiSolver(mode="safe", rtol=1e-8, atol=1e-8, dt_max=1) t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval) - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[-1], 2.5 + 1e-10) + np.testing.assert_array_less(solution.y.full()[0], 1.5) + np.testing.assert_array_less(solution.y.full()[-1], 2.5 + 1e-10) # test the last entry is exactly 2.5 np.testing.assert_array_almost_equal(solution.y[-1, -1], 2.5, decimal=2) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) np.testing.assert_array_almost_equal( - solution.y[-1], 2 * np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[-1], 2 * np.exp(0.1 * solution.t), decimal=5 ) pybamm.settings.debug_mode = True @@ -151,13 +151,13 @@ def test_model_solver_events(self): solver = pybamm.CasadiSolver(dt_max=0, rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 5, 100) solution = solver.solve(model, t_eval) - np.testing.assert_array_less(solution.y[0], 1.5) - np.testing.assert_array_less(solution.y[-1], 2.5) + np.testing.assert_array_less(solution.y.full()[0], 1.5) + np.testing.assert_array_less(solution.y.full()[-1], 2.5) np.testing.assert_array_almost_equal( - solution.y[0], np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[0], np.exp(0.1 * solution.t), decimal=5 ) np.testing.assert_array_almost_equal( - solution.y[-1], 2 * np.exp(0.1 * solution.t), decimal=5 + solution.y.full()[-1], 2 * np.exp(0.1 * solution.t), decimal=5 ) # Test when an event returns nan @@ -173,7 +173,7 @@ def test_model_solver_events(self): disc.process_model(model) solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8) solution = solver.solve(model, t_eval) - np.testing.assert_array_less(solution.y[0], 1.02 + 1e-10) + np.testing.assert_array_less(solution.y.full()[0], 1.02 + 1e-10) np.testing.assert_array_almost_equal(solution.y[0, -1], 1.02, decimal=2) def test_model_step(self): @@ -197,7 +197,9 @@ def test_model_step(self): dt = 1 step_sol = solver.step(None, model, dt) np.testing.assert_array_equal(step_sol.t, [0, dt]) - np.testing.assert_array_almost_equal(step_sol.y[0], np.exp(0.1 * step_sol.t)) + np.testing.assert_array_almost_equal( + step_sol.y.full()[0], np.exp(0.1 * step_sol.t) + ) # Step again (return 5 points) step_sol_2 = solver.step(step_sol, model, dt, npts=5) @@ -205,13 +207,13 @@ def test_model_step(self): step_sol_2.t, np.concatenate([np.array([0]), np.linspace(dt, 2 * dt, 5)]) ) np.testing.assert_array_almost_equal( - step_sol_2.y[0], np.exp(0.1 * step_sol_2.t) + step_sol_2.y.full()[0], np.exp(0.1 * step_sol_2.t) ) # Check steps give same solution as solve t_eval = step_sol.t solution = solver.solve(model, t_eval) - np.testing.assert_array_almost_equal(solution.y[0], step_sol.y[0]) + np.testing.assert_array_almost_equal(solution.y.full()[0], step_sol.y.full()[0]) def test_model_step_with_input(self): # Create model @@ -233,7 +235,7 @@ def test_model_step_with_input(self): dt = 0.1 step_sol = solver.step(None, model, dt, npts=5, inputs={"a": 0.1}) np.testing.assert_array_equal(step_sol.t, np.linspace(0, dt, 5)) - np.testing.assert_allclose(step_sol.y[0], np.exp(0.1 * step_sol.t)) + np.testing.assert_allclose(step_sol.y.full()[0], np.exp(0.1 * step_sol.t)) # Step again with different inputs step_sol_2 = solver.step(step_sol, model, dt, npts=5, inputs={"a": -1}) @@ -242,7 +244,7 @@ def test_model_step_with_input(self): step_sol_2["a"].entries, np.array([0.1, 0.1, 0.1, 0.1, 0.1, -1, -1, -1, -1]) ) np.testing.assert_allclose( - step_sol_2.y[0], + step_sol_2.y.full()[0], np.concatenate( [ np.exp(0.1 * step_sol.t[:5]), @@ -275,13 +277,13 @@ def test_model_step_events(self): while time < end_time: step_solution = step_solver.step(step_solution, model, dt=dt, npts=10) time += dt - np.testing.assert_array_less(step_solution.y[0], 1.5) - np.testing.assert_array_less(step_solution.y[-1], 2.5001) + np.testing.assert_array_less(step_solution.y.full()[0], 1.5) + np.testing.assert_array_less(step_solution.y.full()[-1], 2.5001) np.testing.assert_array_almost_equal( - step_solution.y[0], np.exp(0.1 * step_solution.t), decimal=5 + step_solution.y.full()[0], np.exp(0.1 * step_solution.t), decimal=5 ) np.testing.assert_array_almost_equal( - step_solution.y[-1], 2 * np.exp(0.1 * step_solution.t), decimal=4 + step_solution.y.full()[-1], 2 * np.exp(0.1 * step_solution.t), decimal=4 ) def test_model_solver_with_inputs(self): @@ -305,17 +307,23 @@ def test_model_solver_with_inputs(self): t_eval = np.linspace(0, 10, 100) solution = solver.solve(model, t_eval, inputs={"rate": 0.1}) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t), rtol=1e-04) + np.testing.assert_allclose( + solution.y.full()[0], np.exp(-0.1 * solution.t), rtol=1e-04 + ) # Without grid solver = pybamm.CasadiSolver(mode="safe without grid", rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 10, 100) solution = solver.solve(model, t_eval, inputs={"rate": 0.1}) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_allclose(solution.y[0], np.exp(-0.1 * solution.t), rtol=1e-04) + np.testing.assert_allclose( + solution.y.full()[0], np.exp(-0.1 * solution.t), rtol=1e-04 + ) solution = solver.solve(model, t_eval, inputs={"rate": 1.1}) self.assertLess(len(solution.t), len(t_eval)) - np.testing.assert_allclose(solution.y[0], np.exp(-1.1 * solution.t), rtol=1e-04) + np.testing.assert_allclose( + solution.y.full()[0], np.exp(-1.1 * solution.t), rtol=1e-04 + ) def test_model_solver_dae_inputs_in_initial_conditions(self): # Create model @@ -336,10 +344,10 @@ def test_model_solver_dae_inputs_in_initial_conditions(self): model, t_eval, inputs={"rate": -1, "ic 1": 0.1, "ic 2": 2} ) np.testing.assert_array_almost_equal( - solution.y[0], 0.1 * np.exp(-solution.t), decimal=5 + solution.y.full()[0], 0.1 * np.exp(-solution.t), decimal=5 ) np.testing.assert_array_almost_equal( - solution.y[-1], 0.1 * np.exp(-solution.t), decimal=5 + solution.y.full()[-1], 0.1 * np.exp(-solution.t), decimal=5 ) # Solve again with different initial conditions @@ -347,10 +355,10 @@ def test_model_solver_dae_inputs_in_initial_conditions(self): model, t_eval, inputs={"rate": -0.1, "ic 1": 1, "ic 2": 3} ) np.testing.assert_array_almost_equal( - solution.y[0], 1 * np.exp(-0.1 * solution.t), decimal=5 + solution.y.full()[0], 1 * np.exp(-0.1 * solution.t), decimal=5 ) np.testing.assert_array_almost_equal( - solution.y[-1], 1 * np.exp(-0.1 * solution.t), decimal=5 + solution.y.full()[-1], 1 * np.exp(-0.1 * solution.t), decimal=5 ) def test_model_solver_with_external(self): @@ -375,7 +383,9 @@ def test_model_solver_with_external(self): solver = pybamm.CasadiSolver(rtol=1e-8, atol=1e-8) t_eval = np.linspace(0, 10, 100) solution = solver.solve(model, t_eval, external_variables={"var2": 0.5}) - np.testing.assert_allclose(solution.y[0], 1 - 0.5 * solution.t, rtol=1e-06) + np.testing.assert_allclose( + solution.y.full()[0], 1 - 0.5 * solution.t, rtol=1e-06 + ) def test_model_solver_with_non_identity_mass(self): model = pybamm.BaseModel() @@ -403,8 +413,8 @@ def test_model_solver_with_non_identity_mass(self): t_eval = np.linspace(0, 1, 100) solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) - np.testing.assert_allclose(solution.y[0], np.exp(0.1 * solution.t)) - np.testing.assert_allclose(solution.y[-1], 2 * np.exp(0.1 * solution.t)) + np.testing.assert_allclose(solution.y.full()[0], np.exp(0.1 * solution.t)) + np.testing.assert_allclose(solution.y.full()[-1], 2 * np.exp(0.1 * solution.t)) def test_dae_solver_algebraic_model(self): model = pybamm.BaseModel() diff --git a/tests/unit/test_solvers/test_processed_variable.py b/tests/unit/test_solvers/test_processed_variable.py index 951168b8e6..6fb5d750da 100644 --- a/tests/unit/test_solvers/test_processed_variable.py +++ b/tests/unit/test_solvers/test_processed_variable.py @@ -1,6 +1,7 @@ # # Tests for the Processed Variable class # +import casadi import pybamm import tests @@ -8,6 +9,23 @@ import unittest +def to_casadi(var_pybamm, y, inputs=None): + t_MX = casadi.MX.sym("t") + y_MX = casadi.MX.sym("y", y.shape[0]) + + symbolic_inputs_dict = {} + inputs = inputs or {} + for key, value in inputs.items(): + symbolic_inputs_dict[key] = casadi.MX.sym("input", value.shape[0]) + + symbolic_inputs = casadi.vertcat(*[p for p in symbolic_inputs_dict.values()]) + + var_sym = var_pybamm.to_casadi(t_MX, y_MX, inputs=symbolic_inputs_dict) + + var_casadi = casadi.Function("variable", [t_MX, y_MX, symbolic_inputs], [var_sym]) + return var_casadi + + class TestProcessedVariable(unittest.TestCase): def test_processed_variable_0D(self): # without space @@ -17,8 +35,9 @@ def test_processed_variable_0D(self): var.mesh = None t_sol = np.linspace(0, 1) y_sol = np.array([np.linspace(0, 5)]) + var_casadi = to_casadi(var, y_sol) processed_var = pybamm.ProcessedVariable( - var, pybamm.Solution(t_sol, y_sol), warn=False + var, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal(processed_var.entries, t_sol * y_sol[0]) @@ -27,8 +46,9 @@ def test_processed_variable_0D(self): var.mesh = None t_sol = np.array([0]) y_sol = np.array([1])[:, np.newaxis] + var_casadi = to_casadi(var, y_sol) processed_var = pybamm.ProcessedVariable( - var, pybamm.Solution(t_sol, y_sol), warn=False + var, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal(processed_var.entries, y_sol[0]) @@ -47,13 +67,15 @@ def test_processed_variable_1D(self): t_sol = np.linspace(0, 1) y_sol = np.ones_like(x_sol)[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal(processed_var.entries, y_sol) np.testing.assert_array_equal(processed_var(t_sol, x_sol), y_sol) + eqn_casadi = to_casadi(eqn_sol, y_sol) processed_eqn = pybamm.ProcessedVariable( - eqn_sol, pybamm.Solution(t_sol, y_sol), warn=False + eqn_sol, eqn_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( processed_eqn(t_sol, x_sol), t_sol * y_sol + x_sol[:, np.newaxis] @@ -68,8 +90,9 @@ def test_processed_variable_1D(self): # On edges x_s_edge = pybamm.Matrix(disc.mesh["separator"].edges, domain="separator") x_s_edge.mesh = disc.mesh["separator"] + x_s_casadi = to_casadi(x_s_edge, y_sol) processed_x_s_edge = pybamm.ProcessedVariable( - x_s_edge, pybamm.Solution(t_sol, y_sol), warn=False + x_s_edge, x_s_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( x_s_edge.entries[:, 0], processed_x_s_edge.entries[:, 0] @@ -80,8 +103,9 @@ def test_processed_variable_1D(self): eqn_sol = disc.process_symbol(eqn) t_sol = np.array([0]) y_sol = np.ones_like(x_sol)[:, np.newaxis] + eqn_casadi = to_casadi(eqn_sol, y_sol) processed_eqn2 = pybamm.ProcessedVariable( - eqn_sol, pybamm.Solution(t_sol, y_sol), warn=False + eqn_sol, eqn_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( processed_eqn2.entries, y_sol + x_sol[:, np.newaxis] @@ -99,9 +123,10 @@ def test_processed_variable_1D_unknown_domain(self): nt = 100 + y_sol = np.zeros((var_pts[x], nt)) solution = pybamm.Solution( np.linspace(0, 1, nt), - np.zeros((var_pts[x], nt)), + y_sol, np.linspace(0, 1, 1), np.zeros((var_pts[x])), "test", @@ -109,7 +134,8 @@ def test_processed_variable_1D_unknown_domain(self): c = pybamm.StateVector(slice(0, var_pts[x]), domain=["SEI layer"]) c.mesh = mesh["SEI layer"] - pybamm.ProcessedVariable(c, solution, warn=False) + c_casadi = to_casadi(c, y_sol) + pybamm.ProcessedVariable(c, c_casadi, solution, warn=False) def test_processed_variable_2D_x_r(self): var = pybamm.Variable( @@ -134,8 +160,9 @@ def test_processed_variable_2D_x_r(self): t_sol = np.linspace(0, 1) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( processed_var.entries, @@ -165,8 +192,9 @@ def test_processed_variable_2D_x_z(self): t_sol = np.linspace(0, 1) y_sol = np.ones(len(x_sol) * len(z_sol))[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( processed_var.entries, @@ -181,8 +209,9 @@ def test_processed_variable_2D_x_z(self): ) x_s_edge.mesh = disc.mesh["separator"] x_s_edge.secondary_mesh = disc.mesh["current collector"] + x_s_casadi = to_casadi(x_s_edge, y_sol) processed_x_s_edge = pybamm.ProcessedVariable( - x_s_edge, pybamm.Solution(t_sol, y_sol), warn=False + x_s_edge, x_s_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( x_s_edge.entries.flatten(), processed_x_s_edge.entries[:, :, 0].T.flatten() @@ -211,8 +240,9 @@ def test_processed_variable_2D_space_only(self): t_sol = np.array([0]) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal( processed_var.entries, @@ -231,8 +261,9 @@ def test_processed_variable_2D_scikit(self): t_sol = np.linspace(0, 1) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, u_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, u_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, u_sol), warn=False ) np.testing.assert_array_equal( processed_var.entries, np.reshape(u_sol, [len(y), len(z), len(t_sol)]) @@ -250,8 +281,9 @@ def test_processed_variable_2D_fixed_t_scikit(self): t_sol = np.array([0]) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] + var_casadi = to_casadi(var_sol, u_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, u_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, u_sol), warn=False ) np.testing.assert_array_equal( processed_var.entries, np.reshape(u_sol, [len(y), len(z), len(t_sol)]) @@ -268,8 +300,9 @@ def test_processed_var_0D_interpolation(self): t_sol = np.linspace(0, 1, 1000) y_sol = np.array([np.linspace(0, 5, 1000)]) + var_casadi = to_casadi(var, y_sol) processed_var = pybamm.ProcessedVariable( - var, pybamm.Solution(t_sol, y_sol), warn=False + var, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # vector np.testing.assert_array_equal(processed_var(t_sol), y_sol[0]) @@ -277,8 +310,9 @@ def test_processed_var_0D_interpolation(self): np.testing.assert_array_equal(processed_var(0.5), 2.5) np.testing.assert_array_equal(processed_var(0.7), 3.5) + eqn_casadi = to_casadi(eqn, y_sol) processed_eqn = pybamm.ProcessedVariable( - eqn, pybamm.Solution(t_sol, y_sol), warn=False + eqn, eqn_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal(processed_eqn(t_sol), t_sol * y_sol[0]) np.testing.assert_array_almost_equal(processed_eqn(0.5), 0.5 * 2.5) @@ -297,8 +331,9 @@ def test_processed_var_0D_fixed_t_interpolation(self): t_sol = np.array([10]) y_sol = np.array([[100]]) + eqn_casadi = to_casadi(eqn, y_sol) processed_var = pybamm.ProcessedVariable( - eqn, pybamm.Solution(t_sol, y_sol), warn=False + eqn, eqn_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal(processed_var(), 200) @@ -317,8 +352,9 @@ def test_processed_var_1D_interpolation(self): t_sol = np.linspace(0, 1) y_sol = x_sol[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 2 vectors np.testing.assert_array_almost_equal(processed_var(t_sol, x_sol), y_sol) @@ -333,8 +369,9 @@ def test_processed_var_1D_interpolation(self): np.testing.assert_array_almost_equal( processed_var(0.5, x_sol[-1]), 2.5 * x_sol[-1] ) + eqn_casadi = to_casadi(eqn_sol, y_sol) processed_eqn = pybamm.ProcessedVariable( - eqn_sol, pybamm.Solution(t_sol, y_sol), warn=False + eqn_sol, eqn_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 2 vectors np.testing.assert_array_almost_equal( @@ -347,8 +384,11 @@ def test_processed_var_1D_interpolation(self): self.assertEqual(processed_eqn(0.5, x_sol[-1]).shape, (1,)) # test x + x_disc = disc.process_symbol(x) + x_casadi = to_casadi(x_disc, y_sol) + processed_x = pybamm.ProcessedVariable( - disc.process_symbol(x), pybamm.Solution(t_sol, y_sol), warn=False + x_disc, x_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_almost_equal(processed_x(x=x_sol), x_sol[:, np.newaxis]) @@ -357,13 +397,14 @@ def test_processed_var_1D_interpolation(self): disc.mesh["negative particle"].nodes, domain="negative particle" ) r_n.mesh = disc.mesh["negative particle"] + r_n_casadi = to_casadi(r_n, y_sol) processed_r_n = pybamm.ProcessedVariable( - r_n, pybamm.Solution(t_sol, y_sol), warn=False + r_n, r_n_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) np.testing.assert_array_equal(r_n.entries[:, 0], processed_r_n.entries[:, 0]) - # np.testing.assert_array_almost_equal( - # processed_r_n(0, r=np.linspace(0, 1))[:, 0], np.linspace(0, 1) - # ) + np.testing.assert_array_almost_equal( + processed_r_n(0, r=np.linspace(0, 1))[:, 0], np.linspace(0, 1) + ) def test_processed_var_1D_fixed_t_interpolation(self): var = pybamm.Variable("var", domain=["negative electrode", "separator"]) @@ -377,8 +418,9 @@ def test_processed_var_1D_fixed_t_interpolation(self): t_sol = np.array([1]) y_sol = x_sol[:, np.newaxis] + eqn_casadi = to_casadi(eqn_sol, y_sol) processed_var = pybamm.ProcessedVariable( - eqn_sol, pybamm.Solution(t_sol, y_sol), warn=False + eqn_sol, eqn_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # vector @@ -411,8 +453,9 @@ def test_processed_var_2D_interpolation(self): t_sol = np.linspace(0, 1) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 3 vectors np.testing.assert_array_equal( @@ -455,8 +498,9 @@ def test_processed_var_2D_interpolation(self): t_sol = np.linspace(0, 1) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 3 vectors np.testing.assert_array_equal( @@ -486,8 +530,9 @@ def test_processed_var_2D_fixed_t_interpolation(self): t_sol = np.array([0]) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 2 vectors np.testing.assert_array_equal(processed_var(x=x_sol, r=r_sol).shape, (10, 40)) @@ -511,8 +556,9 @@ def test_processed_var_2D_secondary_broadcast(self): t_sol = np.linspace(0, 1) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 3 vectors np.testing.assert_array_equal( @@ -546,8 +592,9 @@ def test_processed_var_2D_secondary_broadcast(self): t_sol = np.linspace(0, 1) y_sol = np.ones(len(x_sol) * len(r_sol))[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) # 3 vectors np.testing.assert_array_equal( @@ -566,8 +613,9 @@ def test_processed_var_2D_scikit_interpolation(self): t_sol = np.linspace(0, 1) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, u_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, u_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, u_sol), warn=False ) # 3 vectors np.testing.assert_array_equal( @@ -606,8 +654,9 @@ def test_processed_var_2D_fixed_t_scikit_interpolation(self): t_sol = np.array([0]) u_sol = np.ones(var_sol.shape[0])[:, np.newaxis] + var_casadi = to_casadi(var_sol, u_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, u_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, u_sol), warn=False ) # 2 vectors np.testing.assert_array_equal(processed_var(y=y_sol, z=z_sol).shape, (15, 15)) @@ -674,8 +723,9 @@ def test_call_failure(self): t_sol = np.linspace(0, 1) y_sol = x_sol[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) with self.assertRaisesRegex(ValueError, "x cannot be None"): processed_var(0) @@ -693,8 +743,9 @@ def test_call_failure(self): var_sol = disc.process_symbol(var) y_sol = r_sol[:, np.newaxis] * np.linspace(0, 5) + var_casadi = to_casadi(var_sol, y_sol) processed_var = pybamm.ProcessedVariable( - var_sol, pybamm.Solution(t_sol, y_sol), warn=False + var_sol, var_casadi, pybamm.Solution(t_sol, y_sol), warn=False ) with self.assertRaisesRegex(ValueError, "r cannot be None"): processed_var(0) @@ -717,9 +768,12 @@ def test_3D_raises_error(self): var_sol = disc.process_symbol(var) t_sol = np.array([0, 1, 2]) u_sol = np.ones(var_sol.shape[0] * 3)[:, np.newaxis] + var_casadi = to_casadi(var_sol, u_sol) with self.assertRaisesRegex(NotImplementedError, "Shape not recognized"): - pybamm.ProcessedVariable(var_sol, pybamm.Solution(t_sol, u_sol), warn=False) + pybamm.ProcessedVariable( + var_sol, var_casadi, pybamm.Solution(t_sol, u_sol), warn=False + ) if __name__ == "__main__": diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py index e32ce3f69d..48415af6fa 100644 --- a/tests/unit/test_solvers/test_scikits_solvers.py +++ b/tests/unit/test_solvers/test_scikits_solvers.py @@ -79,7 +79,7 @@ def test_dae_integrate_bad_ics(self): solver.set_up(model) solver._set_initial_conditions(model, {}, True) # check y0 - np.testing.assert_array_equal(model.y0, [0, 0]) + np.testing.assert_array_equal(model.y0.full().flatten(), [0, 0]) # check dae solutions solution = solver.solve(model, t_eval) np.testing.assert_array_equal(solution.t, t_eval) diff --git a/tests/unit/test_solvers/test_solution.py b/tests/unit/test_solvers/test_solution.py index 9285d253c6..ecd5f8a2ee 100644 --- a/tests/unit/test_solvers/test_solution.py +++ b/tests/unit/test_solvers/test_solution.py @@ -87,7 +87,7 @@ def test_cycles(self): self.assertEqual(cycle_sub_solution, sub_solution) def test_total_time(self): - sol = pybamm.Solution([], None) + sol = pybamm.Solution(np.array([0]), np.array([[1, 2]])) sol.set_up_time = 0.5 sol.solve_time = 1.2 self.assertEqual(sol.total_time, 1.7)