From 4d0b2bbd7a1317c39ec5f4d88768506881f24621 Mon Sep 17 00:00:00 2001 From: Valentin Sulzer Date: Sun, 1 Nov 2020 21:12:45 -0500 Subject: [PATCH] #1219 start script on speeding up solvers --- CHANGELOG.md | 3 +- examples/notebooks/change-input-current.ipynb | 15 +- examples/notebooks/speed-up-solver.ipynb | 341 ++++++++++++++++++ examples/scripts/cycling_ageing_yang.py | 56 +-- pybamm/expression_tree/binary_operators.py | 18 +- .../interface/sei/ec_reaction_limited.py | 3 +- pybamm/solvers/algebraic_solver.py | 12 +- pybamm/solvers/base_solver.py | 8 +- pybamm/solvers/casadi_algebraic_solver.py | 8 +- pybamm/solvers/casadi_solver.py | 15 +- pybamm/solvers/dummy_solver.py | 4 +- pybamm/solvers/idaklu_solver.py | 6 +- pybamm/solvers/jax_solver.py | 6 +- pybamm/solvers/scikits_dae_solver.py | 6 +- pybamm/solvers/scikits_ode_solver.py | 6 +- pybamm/solvers/scipy_solver.py | 6 +- pybamm/solvers/solution.py | 3 + .../test_binary_operators.py | 22 +- 18 files changed, 474 insertions(+), 64 deletions(-) create mode 100644 examples/notebooks/speed-up-solver.ipynb diff --git a/CHANGELOG.md b/CHANGELOG.md index 8898f5c788..87fc5d6acc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Added `Solution.integration_time`, which is the time taken just by the integration subroutine, without extra setups ([#1223](https://github.com/pybamm-team/PyBaMM/pull/1223)) - Added parameter set for an A123 LFP cell ([#1209](https://github.com/pybamm-team/PyBaMM/pull/1209)) - Added variables related to equivalent circuit models ([#1204](https://github.com/pybamm-team/PyBaMM/pull/1204)) - Added an example script to check conservation of lithium ([#1186](https://github.com/pybamm-team/PyBaMM/pull/1186)) @@ -17,7 +18,7 @@ ## Bug fixes -- Raise error if saving to matlab with variable names that matlab can't read, and give option of providing alternative variable names ([#1206](https://github.com/pybamm-team/PyBaMM/pull/1206)) +- Raise error if saving to MATLAB with variable names that MATLAB can't read, and give option of providing alternative variable names ([#1206](https://github.com/pybamm-team/PyBaMM/pull/1206)) - Raise error if the boundary condition at the origin in a spherical domain is other than no-flux ([#1175](https://github.com/pybamm-team/PyBaMM/pull/1175)) - Fix boundary conditions at r = 0 for Creating Models notebooks ([#1173](https://github.com/pybamm-team/PyBaMM/pull/1173)) diff --git a/examples/notebooks/change-input-current.ipynb b/examples/notebooks/change-input-current.ipynb index d0cf216944..477cde0eb4 100644 --- a/examples/notebooks/change-input-current.ipynb +++ b/examples/notebooks/change-input-current.ipynb @@ -327,7 +327,20 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.8.5" + }, + "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/notebooks/speed-up-solver.ipynb b/examples/notebooks/speed-up-solver.ipynb new file mode 100644 index 0000000000..70e0921b5e --- /dev/null +++ b/examples/notebooks/speed-up-solver.ipynb @@ -0,0 +1,341 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Speeding up the solvers" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook contains a collection of tips on how to speed up the solvers" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[33mWARNING: You are using pip version 20.2.1; however, version 20.2.4 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", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Choosing a solver" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since it is very easy to switch which solver is used for the model, we recommend you try different solvers for your particular use case. In general, the `CasadiSolver` is the fastest.\n", + "\n", + "Once you have found a good solver, you can further improve performance by trying out different values for the `method`, `rtol`, and `atol` arguments. Further options are sometimes available, but are solver specific. See [solver API docs](https://pybamm.readthedocs.io/en/latest/source/solvers/index.html) for details." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Choosing and optimizing CasadiSolver settings" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Handling instabilities" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If the solver is taking a lot of steps, possibly failing with a `max_steps` error, and the error persists with different solvers, this suggests a problem with the model itself. This can be due to a few things:\n", + "\n", + "- A singularity in the model (such as division by zero). Solve up to the time where the model fails, and plot some variables to see if they are going to infinity. You can then narrow down the source of the problem.\n", + "- High model stiffness. Again, plot different variables to identify which variables or parameters may be causing problems. To reduce stiffness, all dimensionless parameter values should be as close to 1 as possible.\n", + "- Non-differentiable functions (see [below](#Smooth-approximations-to-non-differentiable-functions))\n", + "\n", + "If none of these fixes work, we are interested in finding out why - please get in touch!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Smooth approximations to non-differentiable functions" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Some functions, such as `minimum`, `maximum`, `heaviside`, and `abs`, are discontinuous and/or non-differentiable (their derivative is discontinuous). Adaptive solvers can deal with this discontinuity, but will take many more steps close to the discontinuity in order to resolve it. Therefore, using smooth approximations instead can reduce the number of steps taken by the solver, and hence the integration time. See [this post](https://discourse.julialang.org/t/handling-instability-when-solving-ode-problems/9019/5) for more details.\n", + "\n", + "Here is an example using the `maximum` function. The function `maximum(x,1)` is continuous but non-differentiable at `x=1`, where its derivative jumps from 0 to 1. However, we can approximate it using the [`softplus` function](https://en.wikipedia.org/wiki/Rectifier_(neural_networks)#Softplus), which is smooth everywhere and is sometimes used in neural networks as a smooth approximation to the RELU activation function. The `softplus` function is given by\n", + "$$\n", + "s(x,y;k) = \\frac{\\log(\\exp(kx)+\\exp(ky))}{k},\n", + "$$\n", + "where `k` is a strictly positive smoothing (or sharpness) parameter. The larger the value of `k`, the better the approximation but the stiffer the term (exp blows up quickly!). Usually, a value of `k=10` is a good middle ground.\n", + "\n", + "In PyBaMM, you can either call the `softplus` function directly, or change `pybamm.settings.max_smoothing` to automatically replace all your calls to `pybamm.maximum` with `softplus`." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exact maximum: maximum(x, y)\n", + "Softplus (k=10): log(exp(10.0 * x) + exp(10.0 * y)) / 10.0\n", + "Softplus (k=20): log(exp(20.0 * x) + exp(20.0 * y)) / 20.0\n", + "Softplus (k=30): log(exp(30.0 * x) + exp(30.0 * y)) / 30.0\n", + "Exact maximum: maximum(x, y)\n" + ] + } + ], + "source": [ + "x = pybamm.Variable(\"x\")\n", + "y = pybamm.Variable(\"y\")\n", + "\n", + "# Normal maximum\n", + "print(\"Exact maximum:\", pybamm.maximum(x,y))\n", + "\n", + "# Softplus\n", + "print(\"Softplus (k=10):\", pybamm.softplus(x,y,10))\n", + "\n", + "# Changing the setting to call softplus automatically\n", + "pybamm.settings.max_smoothing = 20\n", + "print(\"Softplus (k=20):\", pybamm.maximum(x,y))\n", + "\n", + "# All smoothing parameters can be changed at once\n", + "pybamm.settings.set_smoothing_parameters(30)\n", + "print(\"Softplus (k=30):\", pybamm.maximum(x,y))\n", + "\n", + "# Change back\n", + "pybamm.settings.set_smoothing_parameters(\"exact\")\n", + "print(\"Exact maximum:\", pybamm.maximum(x,y))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here is the plot of softplus with different values of `k`" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEvCAYAAABhSUTPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABXnElEQVR4nO3dd1xWdf/H8deXDQIqgnvvAU40RzlSS3On5t62s8wybyvNxl027l/d1Z1WZmqZ2dYsLXPkyG1uU1FRQRyIIkPmdX5/XAgaDoxxMd7Px4OHXN+zPl8ukDff8z3nGMuyEBEREZF/xsnRBYiIiIgUZApTIiIiItmgMCUiIiKSDQpTIiIiItmgMCUiIiKSDQpTIiIiItng4qgD+/v7W1WrVnXU4UVERESybNu2bZGWZQVca5nDwlTVqlXZunWrow4vIiIikmXGmGPXW6bTfCIiIiLZoDAlIiIikg0KUyIiIiLZ4LA5U9eSnJxMWFgYCQkJji5FcpCHhwcVK1bE1dXV0aWIiIjkuHwVpsLCwvDx8aFq1aoYYxxdjuQAy7I4d+4cYWFhVKtWzdHliIiI5Lh8dZovISGBUqVKKUgVIsYYSpUqpdFGEREptPJVmAIUpAohvaciIlKY3TRMGWMqGWNWGWP2GWP2GmOeuMY6xhjzrjEmxBizyxjTNHfKLdh++OEH9u3b5+gyREREJAdlZWQqBXjKsqz6QEvgUWNM/b+t0xWolfbxADAjR6ssJBSmRERECp+bhinLsiIsy9qe9nkMsB+o8LfVegHzLLuNQAljTLkcrzaPfP7557Ro0YLGjRvz4IMPsmnTJho2bEhCQgJxcXE0aNCAPXv2EBsbS8eOHWnatClBQUEsWrQofR/z5s2jYcOGNGrUiGHDhvHHH3+wePFiJk6cSOPGjTl8+LADeygiIlIIJCewaONfHDod49AybulqPmNMVaAJsOlviyoAJ654HZbWFpGd4hxh//79LFy4kPXr1+Pq6sojjzzCgQMH6NmzJ88//zyXLl1i6NChBAYGkpKSwvfff4+vry+RkZG0bNmSnj17sm/fPl555RX++OMP/P39iYqKws/Pj549e9K9e3f69evn6G6KiIgUbEnxhD7emUpulxhT+kV+fKoLxT0dcwueLIcpY4w38C0w3rKsi//kYMaYB7CfBqRy5co3XLfqv376J4e4qdDp3W64fMWKFWzbto3mzZsDcOnSJUqXLs3UqVNp3rw5Hh4evPvuu4D9sv9nn32WNWvW4OTkRHh4OKdPn2blypX0798ff39/APz8/HKlLyIiIkXVT+/NofrqKLzcbYyfZhwWpCCLYcoY44o9SM23LOu7a6wSDlS64nXFtLarWJb1EfARQHBwsHXL1eYBy7IYMWIEr7322lXtERERxMbGkpycTEJCAsWKFWP+/PmcPXuWbdu24erqStWqVXULABERkVz20OK3WRlVjFfKVKN0+9u4t093h9Zz0zBl7Ne1fwLstyzr/66z2mLgMWPMl8BtQLRlWdk6xXezEaTc0rFjR3r16sWTTz5J6dKliYqKIiYmhnHjxvHyyy9z9OhRJk2axPvvv090dDSlS5fG1dWVVatWceyY/YHSd955J3369GHChAmUKlUq/TSfj48PMTGOPa8rIiJSUNnizvHvH6azPmUZrtV8SOz+Ce1b1XN0WVkamWoDDAN2G2N2pLU9C1QGsCxrJvAzcA8QAsQDo3K80jxSv359XnnlFe666y5sNhuurq706tULV1dXBg8eTGpqKq1bt2blypUMGTKEHj16EBQURHBwMHXr1gWgQYMGPPfcc7Rr1w5nZ2eaNGnCnDlzGDhwIPfffz/vvvsu33zzDTVq1HBwb0VERAoGW2Isx0d1YERsNL92q8mddUcwLB8EKQBjWY452xYcHGxt3br1qrb9+/dTr17++MJIztJ7KyIi/5TNZuO/81bS+c0ncLbZODXlZToMztuLuYwx2yzLCr7Wsnz1bD4RERGRK6WkptL3q2f463wMf7Qaw+PNS+V5kLoZhSkRERHJl1LPHmblhw9yuEIkriWg18NjuKPZbY4uKxOFKREREcl3UpKSOTyqL1WPX2JMl4p49nuKofkwSEE+fNCxiIiIFG0JyUk89tVytsdWIjnVmduDH+ehFvc4uqzr0siUiIiI5BtxF8/S7ftnOJuynw3t7+eT2xvQvFNLR5d1QxqZEhERkXwh8eBazj/cAis2FCeTyvN96xKcz4MUKExly9q1a2nQoAGNGzdmw4YN/PzzzzfdJjQ0lMDAwBw5/jvvvMO8efMAaN++PX+/1URWzJkzh4CAABo3bkzjxo2ZNWsWAGfPnqVLly45UqeIiMjNJCSnsn38VGK2efDar3G80eoj+ge1cXRZWaIwlQ3z589n8uTJ7NixgwMHDmQpTOWUlJQUZs+ezeDBg7O9rwEDBrBjxw527NjB2LFjAQgICKBcuXKsX78+2/sXERG5kXPxMfT6bDqvletFrJcXJSa8R9d6TRxdVpYpTF0hLi6Obt260ahRIwIDA1m4cCFgf/hxkyZNCAoKYvTo0SQmJjJr1iy++uorpkyZwqBBg5g6dSoLFy6kcePGLFy4kGnTpjFs2DBatWpFrVq1+PjjjzMdb86cOTz22GPpr7t3787q1atJTU1l5MiRBAYGEhQUxNtvv51p25UrV9K0aVNcXK6e9maz2Rg5ciTPP/98tr8evXv3Zv78+dnej4iIyPXEH1hJzy9HcdL5S443OoHvol+p3y5/XrV3PZqAfoVly5ZRvnx5fvrpJwCio6NJSEhg5MiRrFixgtq1azN8+HBmzJjB+PHjWbduHd27d6dfv37MmTOHrVu38v777wMwbdo0du3axcaNG4mLi6NJkyZ065a15w3u2LGD8PBw9uzZA8CFCxcyrbN+/XqaNWt2VVtKSgpDhgwhMDCQ5557DrCPOh04cCDT9hMmTGD48OEAfPvtt6xZs4batWvz9ttvU6mS/ZnVwcHBORLKREREriV287ecmzSR8c18eKWeH+/2GE2dSqUcXdYty98jU9OK2z+u9MUAe9uBpRltWz+1ty1+PKPtYoS97a06WT5cUFAQy5cvZ9KkSaxdu5bixYtz4MABqlWrRu3atQEYMWIEa9asydL+evXqhaenJ/7+/nTo0IHNmzdnabvq1atz5MgRxo0bx7Jly/D19c20TkREBAEBAVe1Pfjgg1cFKYCFCxemn8K78uNykOrRowehoaHs2rWLzp07M2LEiPRtS5cuzcmTJ7NUs4iIyK24EJ/EV+8sJz7CnerrXfim+/e0rdbA0WX9I/k7TOWx2rVrs337doKCgnj++ed56aWXsrU/Y8wNX7u4uGCz2dJfJyQkAFCyZEl27txJ+/btmTlzZvo8pit5enqmr39Z69atWbVq1VXtAwYMSJ9cfuXH5YnrpUqVwt3dHYCxY8eybdu2q+rx9PT8J10XERG5riNRp+n0xRCmV2/KH3VbEjBrIbXK+Dm6rH8sf5/mmxaduW3wwsxtwaPsH1fyLXft7W/g5MmT+Pn5MXToUEqUKMGsWbN45plnCA0NJSQkhJo1a/LZZ5/Rrl27TNv6+PgQExNzVduiRYuYPHkycXFxrF69munTp5OUlJS+vGrVqnzwwQfYbDbCw8PTR64iIyNxc3Ojb9++1KlTh6FDh2Y6Xr169QgJCbmqbcyYMaxZs4b77ruP7777DhcXl/R5X9cTERFBuXLlAFi8ePFVDyM+ePBgjl15KCIiAhC99P+YfHgNiT6HKVY1ia5PL6B8CS9Hl5Ut+TtM5bHdu3czceJEnJyccHV1ZcaMGXh4ePDpp5/Sv39/UlJSaN68OQ899FCmbTt06MD06dNp3LgxkydPBqBhw4Z06NCByMhIpkyZQvny5QkNDU3fpk2bNlSrVo369etTr149mjZtCkB4eDijRo1KH7V67bXXMh2va9euDBs2LFP7hAkTiI6OZtiwYcyfPx8npxsPPr777rssXrwYFxcX/Pz8mDNnTvqyVatWZXmel4iIyM2c2bCUmBc+4M1SyYzq0Y73+71a4IMUgLEsyyEHDg4Otv5+X6T9+/dfNTJSkE2bNg1vb2+efvrpXDtGnz59eOONN6hVq1au7L9t27YsWrSIkiVLZntfhem9FRGRW7fv9EneePMnnlr6Dpd8van63RJKlwu4+Yb5hDFmm2VZwddapjlTBdj06dOJiIjIlX2fPXuWCRMm5EiQEhGRIsyy2Hp0NwOXDGJDxc3M6Pk0Nb5eXKCC1M3oNF8umTZtWq4fo06dOtSpk/WrFW9FQEAAvXv3zpV9i4hIEWFZRH38EKcObcCq74SP72mmju5LqWtcpV6QKUyJiIhIrji8dRvJH/xOrVSL3t5deHTCNMr6FK4gBTrNJyIiIrlg+aEdDPxtJ9tK1uFU6Uo88/jLlPUpnFNHNDIlIiIiOScliXUb5jPh0ExspWHJoEnMGnAXxYp7O7qyXKMwJSIiIjnDlsrZ13tTdcMeSt1dkxSfsswc0oViHgX/9gc3otN82bB27VoaNGhA48aN2bBhAz///PNNtwkNDc2xG2G+88476Xcyb9++PX+/1URWrFmzJv2Byd98881Vy+bOnUutWrWoVasWc+fOTW/v1KkT58+fz17xIiJS6Gw7fJawb8KJCSnG2L+CWDroE4oX8iAFClPZMn/+fCZPnsyOHTs4cOBAlsJUTklJSWH27NkMHjw4W/upXLkyc+bMybSfqKgoXnzxRTZt2sTmzZt58cUX0wPUsGHD+OCDD7J1XBERKVzmbV/B8J+nMaXZaEIa3sZ9/3kDH/ei8UgyhakrxMXF0a1bNxo1akRgYGD6o1hWrFhBkyZNCAoKYvTo0SQmJjJr1iy++uorpkyZwqBBg5g6dSoLFy6kcePGLFy4kGnTpjFs2DBatWpFrVq1+PjjjzMdb86cOTz22GPpr7t3787q1atJTU1l5MiRBAYGEhQUxNtvv51p25UrV6aPKF3JZrMxcuRInn/++Sz1uWrVqjRs2DDTndJ/+eUXOnfujJ+fHyVLlqRz584sW7YMgJ49e7JgwYIs7V9ERAq5S+c5+slw3tzxLE4l1lPsTme6LpiNm4e7oyvLM5ozdYVly5ZRvnx5fvrpJwCio6NJSEhg5MiRrFixgtq1azN8+HBmzJjB+PHjWbduHd27d6dfv37MmTOHrVu38v777wP2+0zt2rWLjRs3EhcXR5MmTbL8aJYdO3YQHh7Onj17ALhw4UKmddavX0+zZs2uaktJSWHIkCEEBgby3HPPAfYHHR84cCDT9hMmTGD48OHXrSE8PJxKlSqlv65YsSLh4eGA/UHMiYmJnDt3jlKlSmWpTyIiUjiF/XsoyUsO8eSdpfimWUc+u+8hXJyL1lhNvu5t0NwgguYGXdX22IrHCJobxOoTq9Pbvj74NUFzg5j2x7T0tjPxZwiaG8SdX92Z9eMFBbF8+XImTZrE2rVrKV68OAcOHKBatWrUrl0bgBEjRrBmzZos7a9Xr154enri7+9Phw4d0h9kfDPVq1fnyJEjjBs3jmXLluF7jZubRUREEBBw9d1jH3zwwauCFMDChQvZsWNHpo8bBamsKF26NCdPnszWPkREpGD7aXcoyzd7kZrgTEB0M5YMehs3l6I3TpOvw1Req127Ntu3bycoKIjnn3+el156KVv7M8bc8LWLi0v6w4wBEhISAPvIz86dO2nfvj0zZ85k7Nixmfbt6emZvv5lrVu3ZtWqVVe1DxgwgMaNG2f6uDxx/XoqVKjAiRMn0l+HhYVRoUKFq2r19Cwa58JFRORvki/x+pqFTNo8iFdv68jm/g/T7eN3cXIyN9+2EMrX8XH3iN2Z2t7v+H6mtv61+9O/dv+r2kp7lb7m9jdy8uRJ/Pz8GDp0KCVKlGDWrFk888wzhIaGEhISQs2aNfnss89o165dpm19fHyIiYm5qm3RokVMnjyZuLg4Vq9ezfTp00lKSkpfXrVqVT744ANsNhvh4eHpI1eRkZG4ubnRt29f6tSpw9ChQzMdr169eoSEhFzVNmbMGNasWcN9993Hd999h4uLS/q8r1t199138+yzz6ZPOv/111957bXXALAsi1OnTlG1atV/tG8RESnAoo5w7pXurKhWDuMbS4uGkQzv91ymAYOiJF+Hqby2e/duJk6ciJOTE66ursyYMQMPDw8+/fRT+vfvT0pKCs2bN+ehhx7KtG2HDh2YPn06jRs3ZvLkyQA0bNiQDh06EBkZyZQpUyhfvjyhoaHp27Rp04Zq1apRv3596tWrR9OmTQH7fKVRo0alj1pdDjFX6tq1K8OGDcvUPmHCBKKjoxk2bBjz58/PNLH877Zs2UKfPn04f/48P/74Iy+88AJ79+7Fz8+PKVOm0Lx5cwCmTp2Kn58fANu2baNly5aZJr+LiEjht+PtN3FfauONymf4+MFHea/3A0U6SAEYy7IccuDg4GDr7/dF2r9/P/Xq1XNIPTlt2rRpeHt78/TTT+faMfr06cMbb7xBrVq1cu0Y1/LEE0/Qs2dPOnbsmOVtCtN7KyJSVL226kfWLghl+voPCbv7Xu75v+xNhylIjDHbLMsKvtYyDS0UYNOnTyciIiLPw1RgYOAtBSkRESngTu7gmW2/sTRqLkn1m7O16yxG92np6KryDYWpXDJt2rRcP0adOnWoU6dOrh/n7+6///48P6aIiDjIiS2cndyfOwK8WdrQjc41GzO6u4LUlRSmRERE5Lp+WLiFOls9qeZq4/G73uL+rlm/5VBRoTAlIiIimdhsNkb/8AbrLhbnuQr1qdjlLgWp61CYEhERkavYdn/H61uXsc1pAx5VvPHsMYfOzfN+WklBoTAlIiIi6WxnD3Lm+fEMueTE4u5B9Gg0hIEKUjekO6Bnw9q1a2nQoAGNGzdmw4YN/PzzzzfdJjQ0lMDAwBw5/jvvvJN+J/P27dvz91tNZMWaNWvSH5j8zTffXLVs7ty51KpVi1q1ajF37tz09m3bthEUFETNmjV5/PHHuXx7jaeffpqVK1dmo0ciIuJIKampvPHjcc4c9iXmhBevVHySZ9sNdHRZ+Z7CVDbMnz+fyZMns2PHDg4cOJClMJVTUlJSmD17NoMHD87WfipXrsycOXMy7ScqKooXX3yRTZs2sXnzZl588cX0u6E//PDDfPzxxxw6dIhDhw6xbNkyAMaNG8f06dOzVY+IiDhG8qUYen45njkxn/F8m4c4/+yrdLxXc6SyQmHqCnFxcXTr1o1GjRoRGBiY/iiWFStW0KRJE4KCghg9ejSJiYnMmjWLr776iilTpjBo0CCmTp3KwoULady4MQsXLmTatGkMGzaMVq1aUatWLT7++ONMx5szZw6PPfZY+uvu3buzevVqUlNTGTlyJIGBgQQFBfH2229n2nblypXpI0pXstlsjBw5kueffz5Lfa5atSoNGzbMdKf0X375hc6dO+Pn50fJkiXp3Lkzy5YtIyIigosXL9KyZUuMMQwfPpwffvgBgCpVqnDu3DlOnTqVpWOLiEg+YFmkrpzOn6+04UTSelx8dzP44RbcPqy3oysrMDRn6grLli2jfPny/PTTTwBER0eTkJDAyJEjWbFiBbVr12b48OHMmDGD8ePHs27dOrp3706/fv2YM2cOW7du5f337c8OnDZtGrt27WLjxo3ExcXRpEkTunXrlqU6duzYQXh4OHv27AHgwoULmdZZv349zZo1u6otJSWFIUOGEBgYyHPPPQfYH3R84MCBTNtPmDCB4cOHX7eG8PBwKlWqlP66YsWKhIeHEx4eTsWKFTO1X9a0aVPWr19P3759s9RXERFxrKTY8xx/Yx7Fw1Lp0rkx9YcOZ2DjFo4uq0C56ciUMWa2MeaMMWbPdZYXN8b8aIzZaYzZa4wZlVPF7a9bj/11r34EyYmHHmZ/3XrErFyV3nZ+4Vfsr1uPiClT09uST59hf916HLqjbZaPFxQUxPLly5k0aRJr166lePHiHDhwgGrVqlG7dm0ARowYwZo1a7K0v169euHp6Ym/vz8dOnRIf5DxzVSvXp0jR44wbtw4li1bhq+vb6Z1IiIiCAgIuKrtwQcfvCpIASxcuJAdO3Zk+rhRkMqO0qVLc/LkyVzZt4iI5Ky4xETu/3YNa5PrkYAbAzuOYlSzTo4uq8DJymm+OUCXGyx/FNhnWVYjoD3wH2OMW/ZLy3u1a9dm+/btBAUF8fzzz/PSS9l75tDfH/z499cuLi7pDzMGSEhIAKBkyZLs3LmT9u3bM3PmTMaOHZtp356enunrX9a6dWtWrVp1VfuAAQNo3Lhxpo/LE9evp0KFCpw4cSL9dVhYGBUqVKBChQqEhYVlar+yD56enjfct4iIOJgtldi/fqHLl6PZZvs373TogjVzHs26d3B0ZQXSTU/zWZa1xhhT9UarAD7GnhS8gSggJSeKq/fX/kxtlWbOyNRWcsB9lBxw31VtrmVKX3P7Gzl58iR+fn4MHTqUEiVKMGvWLJ555hlCQ0MJCQmhZs2afPbZZ7Rr1y7Ttj4+PsTExFzVtmjRIiZPnkxcXByrV69m+vTpJCUlpS+vWrUqH3zwATabjfDw8PSRq8jISNzc3Ojbty916tRh6NChmY5Xr149QkJCrmobM2YMa9as4b777uO7777DxcUlfd7Xrbr77rt59tln0yed//rrr7z22mv4+fnh6+vLxo0bue2225g3bx7jxo1L3+7gwYP079//Hx1TRETygGWR/OVo4n/4Dc/b63LB24mX721E4/pNHF1ZgZUTc6beBxYDJwEfYIBlWbYbb5I/7d69m4kTJ+Lk5ISrqyszZszAw8ODTz/9lP79+5OSkkLz5s156KGHMm3boUMHpk+fTuPGjZk8eTIADRs2pEOHDkRGRjJlyhTKly9PaGho+jZt2rShWrVq1K9fn3r16tG0aVPAPl9p1KhR6aNWr732Wqbjde3alWHDhmVqnzBhAtHR0QwbNoz58+dnmlj+d1u2bKFPnz6cP3+eH3/8kRdeeIG9e/fi5+fHlClTaN68OQBTp07Fz88PgA8++ICRI0dy6dIlunbtSteuXQFITk4mJCSE4OBrPlRbRETygfjkVP787AR+R30ZcakE5d5+jztr1nd0WQWauXyPoBuuZB+ZWmJZVqYbJBlj+gFtgAlADWA50MiyrIvXWPcB4AGAypUrNzt27NhVy/fv30+9evX+vlmBNG3aNLy9vXn66adz7Rh9+vThjTfeoFatWrl2jFvx/fffs337dl5++eVMywrTeysiUlCdjr3AoC//g9vmMkz+80vKvv8/6tzWyNFlFQjGmG2WZV1ztCAnbo0wCvjOsgsBjgJ1r7WiZVkfWZYVbFlW8N8nT8utmz59OhEREY4uI11KSgpPPfWUo8sQEZG/S7hIwjcP0WfhKM66/sDphoeosHiJglQOyYnTfMeBjsBaY0wZoA5wJAf2W6BNmzYt149Rp04d6tTJP7f411wpEZH86dJXT3Bm1u+Ma+LP9KCyzOz1GDXLlXB0WYXGTcOUMWYB9qv0/I0xYcALgCuAZVkzgZeBOcaY3YABJlmWFZlrFYuIiEiWnYtN5OtfvGl3xp0yW4qxaMo3VPUv7uiyCpWsXM036CbLTwJ35VRBlmVluoWAFGxZmZcnIiI5LCmOA9EXGLToYS7W6EJishP9/j2BCgpSOS5f3QHdw8ODc+fOUapUKQWqQsKyLM6dO4eHh4ejSxERKTqijpD4QQ9e9q1OcvFQfCv/ysCnFlCmuO4DmBvyVZiqWLEiYWFhnD171tGlSA7y8PC46hE0IiKSuyLX/ciFH5J4psRRHrn3TuYOeE5BKhflqzDl6upKtWrVHF2GiIhIgfVneChvrvXmmQQPziWWZdG9UylVyt/RZRVq+SpMiYiIyD8Uto0/YuJ4aMOzJFWqzsy+k3jt0a6ULF3K0ZUVegpTIiIiBd3xjVx6ry+XkktgNXLH1zuGl0b3pKS3r6MrKxIUpkRERAq4QyedSPnVlzKpFnf4D+GlJx/Av5iPo8sqMnLiDugiIiLiIIv3b+K+lX+yNqAhxyvU581xDypI5TGNTImIiBREf37OxrCDPHd2KVZZi9WDn2XWfV3w8vZydGVFjsKUiIhIQXP2IHGfPEnFEC98ujTCtWQpPh7YDS93d0dXViQpTImIiBQwG6N8cP2jIl7xCQw7ehsjH52Ip6ubo8sqshSmRERECgLLgoQLfLRnI//d+BUVmt/Pg3H7GfvGRFwVpBxKYUpERCS/s9lg6UTO7PqN90p74OQbj3/7pvQe8C7OTnr8mqMpTImIiOR3SbGc+XUV51clcGe7phxrXZbP7ntUQSqf0K0RRERE8rlv9kfwzYF62JKcaJtSju8Hvoqzs36F5xcamRIREcmPkuLh4DKmnYnhm9D3SWgxEu872jN08hiM0YhUfqIwJSIikt/YbPDFfcRu3MyWuo0wJRO4o3Esw/o8oyCVD2mMUEREJL9xcmL7wYqcWO3H0EXudC45kU/vnaQglU8pTImIiOQzU5cvZGJUI2JdPfFp1Zb/6znc0SXJDeg0n4iISH5wPhSWPMl4n6asiP6K5PqN2N/jE4Z3a+royuQmNDIlIiKSHyx7lnNLNtHst1+wUj3oUbudglQBoZEpERERB7Msi2/C2hG480/qu8Tyr54zGHpXK0eXJVmkMCUiIuIo0WHYfMoz5NtpbIkrzWNVWhDYrYOCVAGjMCUiIuIIISuwLRjCG2XuYo/7n3hU8qJSr/l0alLT0ZXJLdKcKREREQdIvXiaM5vcaL9gH0Q1YUzdZ+mrIFUgaWRKREQkjyWlJPPadn+6H/PDMzWRD2qM5o427R1dlvxDClMiIiJ5Zfs8EqrcTo+lr3HiUiwb2z7E1LaVuaNne0dXJtmgMCUiIpIX/pyP9cM49iVWJKKOFy7eFmMeeopWQcGOrkyySWFKREQkDyRWu4tjGyvhGmbR4FIX+oy9m34KUoWCJqCLiIjkFpsNLIvohFhG/rCSn01TEpzdmdy9CwMb3u7o6iSHaGRKREQkN6QkwXf3c9G/Ll2ObyfGOsL+DvfTqeMkAps3cHR1koMUpkRERHLDiY2k7lrMxf1rSGrTHFPcgzf7tiawloJUYaMwJSIikgtiy7dm27bmlD56gjGJxWjyv3m0qlLD0WVJLlCYEhERySmxZyE5juOmGEO+no5L1Xv514VvuOvFp6mpIFVoKUyJiIjkhOhwmNeT5OREBhSvRqxbKG6BSdR96QeqBHg7ujrJRbqaT0REJCd4+JKQXIzD39mouaUWzklVmX3v0wpSRYDClIiISA6ISHTiw+3NsCIthu3ezc99v6RRuSqOLkvygE7ziYiI/FOh6+H4H/xZ815GLX2Y2JrdSbG6M+SVJyjrV8zR1UkeUZgSERH5J+Kj4IsBJF+I49/lt5BaMoISFdcw4qnP8ffxcHR1kod0mk9EROSf8PLjTP2J7FtRiZ7fJ+Ef241F/T9SkCqCFKZERERuRUI0AOuP/cXDfxTnUrIrZUjmu/7PUKG4n4OLE0dQmBIREckKy4Lf34APWrNi51IeWjGSvZWXMnfAM7T89nNKllaQKqpuGqaMMbONMWeMMXtusE57Y8wOY8xeY8zvOVuiiIhIPpCaDCEruHTiDJs/no9FKn5ezvz7ye4U9y/p6OrEgbIyMjUH6HK9hcaYEsAHQE/LshoA/XOkMhERkfzExY19Df7N/hUVuGflTjqeHM6yQbPx89JVe0XdTcOUZVlrgKgbrDIY+M6yrONp65/JodpEREQcK/kS/Pk5WBYLdq1iwKoNrCzflBPVGvDa+DH4emqyueTMrRFqA67GmNWAD/Bfy7LmXWtFY8wDwAMAlStXzoFDi4iI5BLLggUD4chq1h7fy7+jlkFZw85BUxh3X1fcFaQkTU5MQHcBmgHdgLuBKcaY2tda0bKsjyzLCrYsKzggICAHDi0iIpJLjIGmI7gQXY6Id/dgooIo79yaWUN6KEjJVXIiTIUBv1iWFWdZViSwBmiUA/sVERHJe5aV/ula05Kdv/sRFB7C+ON1WDL4HdxddL9ruVpOhKlFwO3GGBdjjBdwG7A/B/YrIiKSt87sh4/aQ9RR3lr/JQ+ufooXbhvB3tbdGPHWJNwUpOQabvpdYYxZALQH/I0xYcALgCuAZVkzLcvab4xZBuwCbMAsy7KuexsFERGRfGv1axCxg4MLJjPH4xjO3rFU7dCee/u9iZOTcXR1kk/dNExZljUoC+u8CbyZIxWJiIg4Sq//cXBXKrEf7KHe7d1w7ejDp/0fwBgFKbk+3QFdRESKtrBt6fOk5u44zLd/FsMzNYkhbi7M7/8vBSm5KZ38FRGRomvLLPjpaWj7NE8mebM8Yh6XWoygUse2DB4/WEFKskRhSkREii7fCuDkzI7fj/KbTyomIJmuTd0Y2mOooyuTAkSn+UREpOiq05UV1lO4L/yDZxafp1/pV3i/x8OOrkoKGIUpEREpOpITYPE4OHsQy7IY/9OHvHjGn/Pu3vh3upMXuvZydIVSAOk0n4iIFB1r/wPb50H4dkaU6cyfMd+SXDeQE33nMuDO+o6uTgoohSkRESk6bn8S25l9rFznT+L2RGy3+TA0sC8DOihIyT+nMCUiIoVbzCkoVhqcnLBcPZkb1pWWv/2bfzm7cnTQHO7t0NTRFUoBpzlTIiJSeEXshJl3wPIpJKUm0ePLx3gl/gyLarYl+olnufdOBSnJPo1MiYhI4RV/Di6dxxaxi4e+eJdjtjV4VthO4KsLaN+wqqOrk0JCI1MiIlJ41biT1KE/sGR5RTp99DtOZ1vzdKPX6a4gJTlIYUpERAqXrbPh1B4AYhJjmfJbLKX2bKN+1DFm1hvMyOC2Di5QChud5hMRkcJjz7ew5EnwKceFsavovugxziUbtrZ/kFc6VKFVt3aOrlAKIYUpEREpPOp0g+rtSa7dh9f/9y3ny4Tj7OHChIdb0aKubn8guUNhSkRECrbkS+DsBk7O4OrBpb4LWHHfGAYd2U3YnX0Y9sQg7qqtICW5R3OmRESk4IqPgrk94eeJYFmEXzzLoC8XsgsfElzdmdSvu4KU5DqNTImISMF1LsR+L6mYCE6dCaHXz4+Q4BTJsbb306frROo1rOXoCqUIUJgSEZGCq1ILGPQFF10rsuTJ90hsUgmnEoaPB3ejXqUqjq5OigiFKRERKViOrAZ3H6jQDIDoCu34pd8o2hzeBtEtaDXrcxqUK+3YGqVIUZgSEZGC4/gm+LwfeBSHh9ayMz6GB398E+dq7ZgcfY7O05+jioKU5DGFKRERKTgqNIVqbSGgLmdSvRn502hS3MLxqF+Cxi9/S8WSXo6uUIoghSkREcnfbKlg2cDZ1f4x6EtOnTjD7u79qFevGQdblOHrAc8rSInDKEyJiEj+lXwJvnsA3Lyh9wdgDEcuxLDghZnce+4EI/dDo7e+p0yJYo6uVIowhSkREcm/oo5AyG/g5ALnJ7IsKoJn1j5FXM3+ODv3Y8xLjxCgICUOpjAlIiL5V5kGMOBz8ClLaITFS78uwPKLo3S5gzz09HuU8HJzdIUiugO6iIjkMyd3wLENGa9rduTwSSfChw3l8W9DqRgzjJ8G/1dBSvINhSkREck/Tu+FT++BBQPsp/iA7/avYtzCLVg2GyXcnFg4+EH8ink4uFCRDDrNJyIi+UdAXajeHty9wbcis//8gbd3vkBKpRrMHzyJNx/ohHcJX0dXKXIVhSkREXEsy4LUZHBxAydn6P8pOLuxZ/Vm/li4FVtzTyp41eE/Y3vi5a5fW5L/6LtSREQcJyUJFo+D1EToOxucnMDFnR3bDpD8+MM8kZKIf8XxTJs0Gg9XZ0dXK3JNClMiIuI4F47DgZ/tN+aMPACl6/HK2hl8uT6eAVVa0Ng5jhefGoa7gpTkYwpTIiLiOP417bc+8CgOpevxydZlLDzyAZRx4fjAN5nUvy1u7rpqT/I3hSkREclb4dshIRpqdLC/rt4OgE1fLcXzfzNx6taSehUbMnNgR5ydjAMLFcka3RpBRETyzpn9MKcbLBwGkYcAsFk2ftl1lNg3XqXJ6YNMiKjO/IGPKUhJgaGRKRERyTv+daBOV3DxgBJVSLWl8sDSf7Hh+CFKtBrBkylHGfrmMwpSUqAoTImISO5KTYHUJHDzsl+t1+dD+7P2jOG7X39n06n1OHkk0OLO0vTrPQZjFKSkYFGYEhGR3JMYA9+MBuMEA7+w30fK2RWAle/OpfbMt2jcshule7fgPz16KEhJgaQwJSIiuSfuLIRtAQxEHQX/miSkJPDe+pXErNjKIFsK95UsRo8efRSkpMBSmBIRkdzjVx0GfwVepaBUDZJSk+j7/RiOxe3lUrNh1O3yCoMf7uvoKkWyRWFKRERy1v4lYAzU7WZ/XalF+qKl//2CsCQvrAAv7m/djMGd2jumRpEcdNNbIxhjZhtjzhhj9txkvebGmBRjTL+cK09ERAqUsK2wcCh8O9Z+Wu8Ki19+n9ofvc5zP57k0Rrv8S8FKSkksnKfqTlAlxutYIxxBl4Hfs2BmkREpKCq0AwaD4G2T0PJqgCcjT9Lv6/H8+pZT854lqB4jx480r6ZY+sUyUE3Pc1nWdYaY0zVm6w2DvgWaJ4TRYmISAESH2U/redZ0v5vr/ft/wKWZTFk0eNEJO0huXocZ++bS7/bazu4YJGcle07oBtjKgB9gBnZL0dERAqUc4dhVif4agSkJtvb0oKUzWZj0bjnKftrZVJj6zD19vEKUlIo5cQE9HeASZZl2W52Wasx5gHgAYDKlSvnwKFFRMShnF0h8aJ9dOrSBfAOACAuKY7P//sDbX/7jgnOrpwetYC7gxs4tlaRXJITYSoY+DItSPkD9xhjUizL+uHvK1qW9RHwEUBwcLCVA8cWERFHKlEZhi+2/+vuDcDes3sZsfQhLiR052TdznTs1Za72ypISeGV7TBlWVa1y58bY+YAS64VpEREpBBITYFfn4MKwdCwv72tTP30xSnJKbz80wISzQXc/HbQbvo7tK1bxkHFiuSNm4YpY8wCoD3gb4wJA14AXAEsy5qZq9WJiEj+cnApbJoJ7r5Qq5N90nma5MQklg5/jP4njnOgWx/e6/cw7esoSEnhl5Wr+QZldWeWZY3MVjUiIpK/1e0Otz8JtbtcFaQ2hG/iqyXn6HVwN77J8XzU7C6a1ynnwEJF8o7ugC4iIjd2eBWUDYJi/vYr9TpNu2rxr0dX8PSap0hOrczuOx/izY5Vad7ldsfUKuIA2b41goiIFGK7v4HP74WFwyAlKdPi+Nh4Vs/dTWqyFy7JVXjriZ40uaedAwoVcRyNTImIyPVVvR18ykG1tuB09a+M+PhEVvUbweDQvZxtO4jHp44jqGIJx9Qp4kAKUyIicrWEaPAobv/cpyw8shE8fNMXW5bFjB0fsWRbEk3dy1HK4xgTRnajnoKUFFE6zSciIhlObIb3m8P2eRltVwQpgNXH1zNj1/scd57NL7fdQamFX1OvddM8LlQk/1CYEhGRDGcPQOxp2Ps9WJnvrXz+zDmOPfUZzidvx+PCYL564B5q1anigEJF8g+d5hMRkQxNh4G7D9Ttlv6MPYBkWzJnYmJYN+pxbju8nfFxzeg4/1Eq+Xk5sFiR/EFhSkSkKIs9Az8/DV1eB9+0+0I16H3VKompiTyx4im2nDiBc41eTI6Po+PbLylIiaRRmBIRKcqWTYZ9i+yn9AZ8ds1VQiLC2BC2i1STgF91D1q9/CVli3vkcaEi+ZfClIhIUdZlOlg2+7/XcPLgMc4MG0vTui0JbdaYhWPupbSPgpTIlTQBXUSkKElJhJ1fZkwu9w6A/p+Cz9XP0Iu8FMlX+5Yy89VPKRt9ioFHdvH1iN4KUiLXoJEpEZGiwrLg874QuhaS4qD5mGuuFpccx5AlwzgZF058jZG4ewzl4RcewK+45kiJXIvClIhIUWEMNB0O549BhevfF+rUnuPEhdYi1ddGg1L1GTexI74ernlYqEjBojAlIlKYWRZcOA4l0+4F1fA+qNsd3DKPMlmWxaHNuzn30P085VGcbwY+yYcPdcLbXb8qRG5Ec6ZERAqrpDhYOBQ+vhOiwzParxGktp7aypAlo3nkh50kGGfcfbz5cEQbBSmRLNBPiYhIYeXiYQ9UqckQdQSKV7jmaqm2VKaue5kTcUdI9PPjm+FTeHNMW7x8iuVxwSIFk8KUiEhhk5oMzq7g5Ax9Z0HiRfCrft3Vdy/fQJWl9QmpX5nbA+7jP0OC8XB1zsOCRQo2hSkRkcIiKR5+ngipSXDvR/YJ58X87R/XEHI+hAuHUzFPPcajKUmUr/Ysk4Y2x81FM0BEboXClIhIYXHxZNoDim1wLgT8a1131Rk7ZzBjxwxSTw2mZ832tHSJYdLE+xSkRP4BhSkRkcLCvyb0/RhKVLlhkAIIO5eAzYIkWwKJg0fS7d4gXFx0ak/kn9CfICIiBVXyJVj8OBz6LaOtbjcoG3jDzdZ/voimU37FHHiAvrXv5a1+jRSkRLJBI1MiIgXVjvmwfS6ErIDH/wQXt+uuGp8cz7t/vktt594U++/b1I85zb+SYhjUJwgnJ5OHRYsUPgpTIiIFVbPRcHovBI++YZAC+Pemf7P48GJSYv7Eu9VYnnE6wqA3JylIieQAneYTESkoEmPgl+fs/wI4OUH3t6Fs0E03DTx/B6nxVUk804X+9wTT7+0pODnpV4BITtDIlIhIQfHDw7D/R0iIhl7v33T1mKQYfNx8WP7WxwR98g4tm9xHw9FteKJjLYzRiJRITlGYEhEpKO6cCrFnoPXjN11199ndPLbyMVqVGEHqH38x2rJxb2U3uneqnQeFihQtClMiIvnVxQgI+Q2aDrO/DqgNo3+x34zzJnZH7iYqIYofDi0jofZwmvW4k4GjuuVywSJFk8KUiEh+lBgLH7WD2NNQohJUb29vz+LpOc8fYnC+1I+YpCa81CuQga2q5lqpIkWdwpSISH7k7g3Nx8KxPyCgbpY2+e7Qd7Sv1J7Vr31GvYUzmVayEpf+bziDFKREcpXClIhIfhG+HYwTlG9sf33HU3DH0/ar9m5i/v75TN88HX/Xmly8cC+vF/PHu/8Aereqlrs1i4hujSAiki8cXgmfdIZvx9gfWAzg5JylIAVwV5W78HEux4nQZlwoVor4mfPoNGFMLhYsIpcpTImI5AeVW4N/bajZOcvzopJTkwGw2Wz88cTrBP7UAmKb8v6gJvRsrhEpkbyi03wiIo5gWbDnW6jX0373clcPuH8luHpmafMz8Wd45LdH6F/7Ps78mEyntT9RxcmFgY8O5s6gcrlcvIhcSSNTIiKO8OMT9lN6v7+e0ZbFIAWw5dQWDpw/wDubZ/Ofc97MCepO/JRXufP2+rlQrIjciEamREQcoeEA+OsnKBv4jza/q0JnFqbsYe2Rani4utHztUm0rumfw0WKSFZoZEpEJC9cPAn7FmW8rtoGxu+CBn2yvIvlx5ZzOu40SQmJLB38IPd99CsBKc7MGdWCNgpSIg6jkSkRkdwWFwkftILkePsk89L17O1uxbK8iyVHljB57WRql6xDhfD76REagl9iDB92KkuT6qVyqXARyQqFKRGR3FbMH+r3hJjT4FnyH+2iTfk2VPGtyqWoxiw+lsyujo/x7l2VadipVQ4XKyK3SmFKRCSn2WywZRbUuBP8a9rb7nkLnN2yfNsDgMTURNyc3DDG4JboSvDmLsyxqlLSy5X3x9xBYIXiudQBEbkVmjMlIpLT1v0Hlk6ExePst0AAcHG/pSB1Ku4Ug38azLx984iNS2BN32EM+P59ep7dyYIHWipIieQjNw1TxpjZxpgzxpg911k+xBizyxiz2xjzhzGmUc6XKSJSgASPgbJB0OqRWwpQV9oduZuD5w/y9YFvGPHZJpYXr8l5T1/GPdSDumV9c7hgEckOY13+q+l6KxjTFogF5lmWlekaXmNMa2C/ZVnnjTFdgWmWZd12swMHBwdbW7du/Ydli4jkI6HrYNdC6PFuRniyrH8cpC77+q9FzFvhye4TyZQr7sH8gQ2oXq1sDhQsIrfKGLPNsqzgay276ciUZVlrgKgbLP/DsqzzaS83AhX/UZUiIgVRUjx8PRK2z7Pf0fyyWwxSlmXxxf4vOBFzAoBz4WdIfOZ7joVcoGJJT756sJWClEg+ldMT0McAS3N4nyIi+YvNZv/XyQncvKDLdDgXAvV6/ONdLvhrAa9tfo2vD37NjA6fs+r+CTQ7so1nkpLo8MInVCiR9buji0jeyrEwZYzpgD1M3X6DdR4AHgCoXLlyTh1aRCTvnNoNSyZAsxHQZKi9Lahftnfbo0YPFh9ezL3VhzB01lbO1e7Gv1KS6Pi/6ZRTkBLJ13Lkaj5jTENgFtDLsqxz11vPsqyPLMsKtiwrOCAgICcOLSKSt87sh7DNsOF/GSNU/9CpuFNcnrfq4+bDW8H/Y+bPPoSciaV0lQrc+e1nlKuhPzxF8rtshyljTGXgO2CYZVkHs1+SiEg+YrPB2Sv+awvqD11eh9HL7Kf5/qF14eu4d9G9fLTrIwCO7wvhSI/e1Nm6gvrlfFnwQEsCfNyzW72I5IGbnuYzxiwA2gP+xpgw4AXAFcCyrJnAVKAU8IGxT7hMud5sdxGRAuXSefjsXog6Ao//CV5+9onlLR/K9q5tlo3Y5FgOnD/A0cgYZr71JSNjIukdsY0WoyZTsphbDnRARPLCTcOUZVmDbrJ8LDA2xyoSEckvPEqAuw+4etoDlZdfju26bcW2zOkyBx9Tk0EfbeZU6SZ43eXJE1NGUcLXK8eOIyK5T3dAFxG5LCke1rwJMafsr42B3jPgsa1QMXsD7hGxEYz5ZQxHo4+mtxULdeHBd1dw6mICLar58dSbT1AiIOcCm4jkDYUpEZHLlj4DK1+xf1xWvAK4e2d71x/u+pDNpzbz1ta3ANi3/k8uPDCap359j44V3Jkzqjne7npcqkhBpJ9cESnaUhLtz80DaDPefrVewwE5fphnmj+Dm7MbjzZ+lN1h0Tz642Ged/XCKuXPeyNb4eWm/45FCqqbPk4mt+hxMiLiUFFH4OeJ9jlR/efk+O7jkuNYeGAhIxuMxMlknATYfvw8I2ZvJiYhhZ6VPXhjZGs8vHQfKZH87kaPk9GfQiJSNDm72Z+p5+wOcZFQzD/Hdm1ZFo/89gjbz2wnxZbCAw0fAGDbklXMW7CKmEq30TWwLG8NbIKbi2ZbiBR0+ikWkaIh5hRs/jjjdfGK0H+u/ZYHORikAIwxPNjoQWqXrM1dVe4CYMOWvzCTxvPQtq952CeK9wYpSIkUFhqZEpHCLzUZPmwHsacgoA5Ua2tvr9Mlxw6RmJrIX1F/0SigEQCty7fmtu634ezkzJqDZ7l/0VG61utCB5cLPPXMIFycFaRECguFKREpnOKjwLOk/fYGzq7QfCyc2gneZXP8ULFJsYz+ZTShF0P5stuXVC9RHQBnJ2dW7j3JQ1/sJCnVhufgIXTvFYizgpRIoaKfaBEpfDZ8AO80hL+WZLS1fRoGfA4BtXP8cMVci1G9RHX8Pf1JtiWnt//+ydfEjxmOZ/xFhreqwr97BylIiRRCGpkSkcLH2RWSYuDoWqjXw95mf9xVjolPjichNQE/Dz+MMUxtOZVUKxUfNx8AfvzzBOajmdSIPslk11AG9LwPk8M1iEj+oDAlIgXbheOw/r9Qvgk0GWpvazLM/jqbdy2/nsMXDjN+1XjKFivLh50/xMk44eWa8QiY77aH8fTXuyjeaizPeZzgvlefVpASKcQUpkSkYAvfBltmQYkq0GgwODmBq0euBSmA4u7FiU6MxsXJhaiEKPw9M64GXPTDOp7aFI1lwcjuwfTtOFBBSqSQU5gSkYLDsuDo7xB7Fhr2t7fV6wmtHrOPSjnl3nykCwkXKO5eHGMM/p7+fHzXx1TxrYKHi0f6Oste/R815/2Pexr2psEjo3mkfc1cq0dE8g/NhBSRguPkdpjXy/4MvaQ4e5uTM9z9byhdL9cOu+zoMu757h5+PPJjelsdvzpXBanZ646yZEc4Tlh0q1NKQUqkCNHIlIjkXzGnIXwr1O1mf12+KdTuYj+Fl4ePwkpMTSQmOYZ14evoWaNnpuUzfz/M9KV/QfU2tOvbif4D7syz2kTE8RSmRCR/io+C/za0h6YJ++x3KTcGBi/M9UNblsXp+NOULWa/J1XPGj0J8AygVflWmdb97oX/8mF0aYynL6/2CaJ/i8q5Xp+I5C8KUyKSPyTGQuhaqNPV/trLD6q1S7vNQWyOP/LleuKT43l23bNsP72d73t9TynPUhhjaF2h9VXrWZbFN1PfIfDrj3jFtxyX3v+UfgpSIkWSwpSIOF5qMrzXFGJPw7jtUKqGvX3gF+Cct/9Nebp4EpscS5ItiYPnD9LKM/NolGVZTF/6F1/GluO14uXwHD6Kbi2q5GmdIpJ/KEyJSN67dAEOLIVGAzMe91KjI0QdgcSLGevlUZA6EHWA0l6lKelREmMML7V+CSfjlH6a70o2m42Xluxnzh+huHj6Yj6cx52NK+ZJnSKSPylMiUjesiz4sC1cOGYfgarUwt7e47/g4pbn5Sw+vJip66fSrXo3/n37vwEo713+muumpqSyZPR4oi4Vw612O/43pCmd65fJy3JFJB/SrRFEJHed2gMrXoKUJPtrYyDwXvt8KK64maUDghRA44DGuDi54O3qjc2yXXe9VJvFO/+3kNqbf2P0vp+Z1a2ygpSIABqZEpHc9t39cGYfVGiWcYuDji/k+LPysmpP5B7Wh6/nwUYPAlDZtzJL711KgFfAdbdJSbXx9Nc7+SGqOMeb9mPwvW1o27pBXpUsIvmcwpSI5JyNM2HvdzDoS/vVeADNx8LpPeBXPWM9BwWp6MRoRv8ymksplwguG0yzMs0AbhikkhISmTxvAz8cicPLzZlBr4zntuql8qpkESkAFKZE5J+7cAJKVMp4fXApnNgEfy2BpsPtbc3HOKa2NLFJsRRzLYYxhuLuxbk/6H5ikmOo61f3ptsmxF9i+cCx3HUyjE0dx/HfMXfSrIpfHlQtIgWJ5kyJyK2zLJjdBd4JtF+Bd9ntT0L/ORDYz2GlXWnBXwu465u72HByQ3rb/Q3vZ0KzCRRzLXbDbROSU3lyzgacT0XglxTLB10qKUiJyDUpTInIjdlscGAZ/PKc/XOwn6bzrQDuvnD2YMa61dtDgz7g5uWQUv/uUsolYpJj+D3s91vbLimVsXO3svT4JV7v9Che//uQhh1uy6UqRaSgM1YePt/qSsHBwdbWrVsdcmwRuQHLgphT4Fsu4/XbDeBiONy/Cio0tbfHngGP4uDi7rharxCXHMeXf31JZd/KdK7SGYCElAR2nN1By3Its7yfmPMXeefl2XziVgt/bzfmj21JnbI+uVW2iBQQxphtlmUFX2uZ5kyJSIbEWJjRCuLOwaRQ++0KjIHbHoTkS+BdOmPdKz/PB1YcX8E729+hqm9V7qx0J85Ozni4eNxSkLoYl8D6vkPpd/IQF1sO4sE3JlCztHcuVi0ihYHClEhRFXsW1r8DCdHQ6317m7s3uHlDUpx9LlTptEnabZ5wWJnXE3I+hJNxJ2lbsS0AXat1ZV3YOnrX7I2TufUZDNHxyQyfs43ypRsxJDqSh8ffR3UFKRHJAp3mEykKEqLhSNq8ofo9M9perwrGCSYdswcpgIsR4F0GnPLvlMq95/YycMlA/D39+bXvr7g6u2Zrf+fjkhj6ySb2nrxIxZKefDG0EZUr6PYHIpLhRqf58u//liLyz9hS7XcdjzyU0RaxE74aBmveyGjzKA5d34Ch310978m3XL4LUhGxEaw4viL9dX2/+jQMaEjHyh2JT4nP1r7PHDvJr31HcPJIGFVLefHVg60UpETklug0n0hBZlkQfQLcfcCzpL1t4wz49TloNtL+vDuACsH2K+2qtLFvc/mmmS3ud0TVt+Rs/Fm6fNcFd2d3VvZfibebN8YYPuv62T86nXelMxcT+PXhSTQL3cFEY+j40mzK+HrkUOUiUlQoTIkUFDYbnD8KJSrD5dNaPz4O2+fZQ1Ozkfa2isFQvHJGuAL7rQqGL8rzkm9Vii2FzRGb2Re1j7FBYwH73clblWuFj5sPscmxeLvZT0dmN0hFRF9i8MebOFenB5OModOM1xWkROQfUZgSyY8unYfzoVC+SUbbxx0gYgc8uAbKNbK3BdQFTz/7lXaXVboNntydl9Vmi2VZmLSRsksplxi3chzJtmR6VO9BmWL2Bwl/0OmDbIenKx0/epIhX//FiahL1K9ajq6vzMKvmGMetCwiBZ/ClIgjJcba5zPZUqB6O3tbUhy8Xg2cXOC5iIxRqFI1IPY0xEVmbN/8fmj5yNXPunPQc+9u1fmE87yy8RXCYsP4stuXGGPwcfNhQN0BeLt64+KU8d9TTgapoxu2c+bhBwmu2QG/9j2ZN/o2intlbwK7iBRtClMieSV8Gxz6DSo2g5qd7G0RO2HOPVC+KTywyt7mVgxK17NPCo8/Bz5l7e29Z2S+QaZLwRhNSUpN4s8zfxKVEEXXal0B8HHzYdOpTUQnRhMWE0YlX/sz/p5p/kyu1XH4bCzvz1rGAwmxtI85yrRRwQpSIpJtClMi2ZUUDxdPgn/NjLYVL0PIb9D9/6BCM3vbsQ2w+lX7aNLlMBVQ1z45/PJdxS97+I/MI0z55E7jWRGbFEt0UjQVvCsAEB4bzthfx+Ln4cfdVe/GyTjh4uTC63e8ThXfKlT0qZjrNR08HcPgjzcRGdCIYr1LMem5YfgU0xwpEck+hSmRG7HZMm4TEB0Gu78GjxIQPCpteSpMr2T/9/nTGYEn6oh9flPkoYwwVaU1tBkPVe/I2H+xUnB/xiX/6QrIqTqA5NRkElMT0yeGrwlbw2MrHqNNhTbM6DQDgKq+VWldvjU1S9QkISUBL1f7s/vaVGiT6/VZlsX+GbN58lgxIo0Pt9f057nhXfB0c871Y4tI0ZC/biYjkhdSU65+/dfPsP6/9mfNXbb+XXitMqx/O6Mt5hT8Ng22zs5oc3KGElXAr5r9lNxlt4+HMcuhzj0ZbRWaQucXoVannOxNnrEsi9Nxp4lPzriv0+f7PqfFFy34dO+n6W21StTC2ckZi4wbAhtj+LDzh0xsPjE9SOWVvR/Mxrz7FuNXzuTOmiWZNSJYQUpEcpRGpqTgi4u0Bx2fcvaRHrDftHL3V+BfG5oMtbddOg//bZR2x+/QjO3/eBeOb7CPIF1+3pyTMyRGQ8zpjPVKVoNWj4F/rauPP25b5pGky1fbFVC7z+7meMxxulbrmj75+8nVT9qff9f+HTpW6QiAv5c/KbYUzl3KCJJli5Vlw6ANeLg4/hTa9uPneSisFM+VqMjB9r2ZMfI23F0UpEQkZylMSd67/AijywEk5hRcOG5/hEnJKva2ixH2+yd5loTbHsjYdl4vuHDCfmrs8n2Ufp0CO7+AXv/LCE7nj9pHm+rck9Hm7guJMYCB1OSMq+Tq9bRPAC92xYN7mwyFRoOuvldTsVJw978z96eAnJKzWTaiE6MBKOlh79epuFN8uOtDPJw9mNRiUvq641aO41zCOZqVaUbZYvYJ8BW8K1DcvTixybHp67Wv2J6NgzdSzLVYepsxxqFByrLZiF25kv3l6zLqyz3EWa788tir/HdQM1ydNRgvIjnvpv+zGGNmG2POGGP2XGe5Mca8a4wJMcbsMsY0vdZ6ks9ZVubTX+cOw+l99nlDl4Vtgz3f2QPNZRG77BOu93yb0ZYQDZ/3hQWDr97nnO7wsj+Eb89o2z4PPukMf36W0RZ/zj5Ze9unV28fdQSiDttHmS4rWcU+kfvKX+BlAqHjVGg6PKPNyRkmHoYpkRlBCqDVI9DlVQiondHmURy8/PJ9ULJZNnae3cn68PVXtX+x/wueWfMMf0X9ld72+b7PabuwLR/u+vCq7b85+A2/hv561fZ3VLyDzlU6k2LL+J4Y32w8awespVfNXultHi4eVwWp/ODoxH8R9tg4vnnhXeKSUunVuDzvKkiJSC7KysjUHOB9YN51lncFaqV93AbMSPu34LPZwEq1nxZySjs1kJoCiRftrz2KZ6wbHWYf7ShRJWPC8oUTcCkKiley/2IG+7ycM/uhmD+UaZC2z2T46yf7L+76Gb+o+OsniA6Het3Bt7y97fgmCFluvzFjrc72tpjTsPJl+zE6v5Sx/Q+PwoVj0OdDKG6/qoq1/7GHlzueyggaR36HeT3tE6NHLsnYfkYbSLkEz0bY76ANsOE92Ps99P0EStgvZefsX7D2LQjsB4F97W3GyX41W9qk5Ku/rimQFJPxukQV+xVt3mUy2nzKQduJGf2+bOACexDyveLqr/b/sn9cya+avY9/d/l9yEGJqYkkpSbh4eKBq5M9pEUnRhN5KRJfN18CvAIAiE+OZ8upLbg6u9K6fOv07b8++DWRlyIZXHcwxd3t31PfH/qeVSdW0admHzpU7gDApohNTFg9gaalm/Jex/cAMBhGLB1BqpXK9qHb0x/4u/X0VpYfW07Hyh2p61cXAH9Pf3zdfDFkBMQArwCeve1Zyhe7+uv8cpuXM/Xzct/yE1tcHEnHj2M8PDhXsiwf/n6YgxfK8LBncaJcvRjUojKv9A7E2Sl/h2IRKdhuGqYsy1pjjKl6g1V6AfMsy7KAjcaYEsaYcpZlReRUkf/ExiPnOPjZBErs20R8lRocavsmAD7ntlPj55dJ9XRjd7+F6es3/PY+nOMTOXrXJC6Usf+iq7NuAl5HQ4ip04SQFi8CUObY95Rf8ykpvj7s7jU/ffvGX/TGpNrY0/cjkr3sp0XqL78f91OniWrSmWOB4wCovG8GpbYtJal0AHvv/gQA58TzNPzqaSxj2DE0Yz5O0JLncTkfQ0Sbk5yqfh8ANba+g+/+LcRX/p0D7ewByTdyG9WX/oLNw5VdiQMz+rRoNU7xSYSmLk3vU+11v+F1NJ7YyMWEnLQ//Drg2DrK7y9Nyslw9vrvS9++wV/+kGJx8IcNJHuVs9e/8SyeEdWIWraFs8fq249/6AD+IXWJi73AaQ/79rbkOKocaYTNyZnwJRn79N7tR7Fzt3FqRRhmn73d6a8j+O9x4vypgySfSds+5gRlfllNkpcP0adaZmy/bCoesRc402kCTiWqA2B2zMbv0FaiGrTHqm//OqWe3k7Amo+I8ytDYscX07f3/f4RnFKSier1H5zcfO2NG17HN+wvopr1xaV6dwCSQn/Gb/MCYsrXwNw+zb7PlFh8vh2DzdmZhH5fpO8zZcXjeEeeJPqOx3Av3x6A+H0f47f7F6JrNMc92H76LOniQXyXPkuSly8remRMYrd++g+esRf5z91lcC1hDz6Xti+m0qFN/NYA1gXaQ2bC6d30Xn2eeL89vHwp42s66vcSuKRYvOq9BSf3EgBU3ezJM+FB7Io9zeYq9nV9w2w8s7010X9583La19mkphD003bOO+3k5e4B6fus/scyvCMjONLqbmID7EHL//Beyu/dTGT1+pwMtP+95BZ3kborviXJy4e/OvVL377Oyu9xj73Agfa9SfS1n1Ist28rASG7OVWvGWdqNQSg2LlT1Fi/lDi/Mhy+PWOyfuBPn+GcnMSee4aS6ma/QrLG+qX4H91PSJuunKtWz9627mea/PAJ+5rdyb+qdiM51QL/OpR84r+M71yXwApX/MEjIpJLcmLOVAXginM+hKW1ZQpTxpgHgAcAKleunAOHvr69Jy8Se+QozXcattgi+MTpKACNk0K4e6cTkcVT+GTd0fT1Pz2QStkoFz4vd4iNh+zBYerRcwTtdGWl20k+SbKve09sKHfsdOVo+cSrtv92rzNeic68XCmE0y72R3u8ejSO+gdd2eIbxtwL9nUHno9gxE5Xdta8lL59MdtF+ux0JckFxl2xz3dCUqkR5soi/5MsP2lvHx0RTf+drqxMjecTZ3tbjcSTvL/ThdPFzVU1zTpgqBDlwntlT3HQ097+8JFkeu5yYb5rEp8n29saxsTw+k4X9pXnqu2/2Z1KsUSYUvEkcS4JAEw6FE2Lg4n8n+dZlp+zr9v+3Akmbb3I+hon+cTF3uZmi2PR5rMkucBLZTP2+frOvTQMi+c5r51sP2z/Jdv35Hbu3RzGj7EWn8S1B6Bc4iFmbwznVHEnRvlmbP/R9kNUikrhoeL7OOZpH224/+huOuw8zdzUnXwZ1TytT/t4fdNZ9pWP4Sn3K/q09TTFEqFfwEHiXOyT1SftO0z7gxf5P9f9LD/ZIK1Phxi8OZ71NU7wCpf7FM+iLSkkuaTQ68o+7b5AwzCL53xC2X7E3n5v+Gm6b7FYHH+WGQn2tnKXTjJ7i8Up3wRGlczY/sM/U6l8weLBEuEcL2YPDmP+SqX7XxazE5KZn/a90yAylre2WOwp5cFEz4ztv94UjXdKAv1KnSTOzX768+mtoQSH/cmbl+qy8oS9n3eE7+SuLUtYU74hnyRWtfcpNZlFa5aQ6OTC+BIZt2x4bd0qakeG8GlqBXYEJALQO2Qr7fcsYV9YNJ9csM8vKxcXyew1S4jwKsUnHs0y+rT2VyrHnOE/LnU57mv/42L03g20ObSadZE2vj7tk9anI3Rdt4Q9parxCfUy+rR2Gd4pCUz2bUWcm2dan7bTJGw7C50qsSLcfkq36Rkn/HzKsCPWiVSbRc9G5Xm0Q03qlPVBRCSvGOvyZOAbrWQfmVpiWVbgNZYtAaZblrUu7fUKYJJlWVtvtM/g4GBr69YbrpIte8KjWb/0C1z2rSW5aiAuTewjNqkXI3BeOQOblzdOd2XcaTn1t7dxirtIarsxuJSwn0JK3P0jrkd2kFy3Ne517FcvJZ89hPMfX5LqVwbXOzImRqf+/DomOQmry5M4u9tPbSVt/RKX8EMkN74L9yr2v+STTmzDZfvPJJerhnsL+8RoKzkRlr6J5eSEU/dn0/eZtO4TXKMiSG7ZD7fS9hGLxJDfcdu3luQqQbg1sp8STI05jdOqj7F5euPceXxGTSvexynuAqntRuGSdpovac8SXI/sJLluS9xqZ/TJdcNXJJcsg+sdY9O3t35+E5OSRGqX8Ti72efFJG9diGvEYZIbdsK1Sgt7W9ifuP25jKQy1XBtMTC9T07L3sZycoZuE9P3aVs/F5eoUyTd1huX0nUASAlZh/v+P0iq0gDnht3s68WexW31XGye3tg6PpS+vdOqWTjHXSSh7RCcfe0jNtbeFXgc3U1CnRaYWvYROCvyGF6bfiS5ZAAprQekb+/16yc4paRwsdMwnNJOXbrt+A3PiGPEBrYhtZL96+xy8jC+u/4gsUwl4pq0t2+cnEyZld+BsyunO92bvs9SW9fiHhVJZLM2JJWyh4xix0IofmA3cZWrE13XfmWfc3wsZdYtJ8XTizN33J2+fen1y3GJi+VM646keNtHy3wP7sE79BAXazUgtpp9LpdbVCT+W9eSVLIUkc3bpm9fdtUSnJKTiejQDcvVfkf0Eru34hVxgvOBzbhU3v6Hi+epMEru2sKlMhU438j+3pmUFMqt/BHL2ZmIjj3zvk/nI/HfspakEn5EtmiX0afVP2FSUjjVrmt6n3xC9uF+7gwx1euSGFCWKzkZQ/s6AVQPuMZpZRGRHGCM2WZZVvA1l+VAmPoQWG1Z1oK01weA9jc7zZfbYUpEREQkp9woTOXE5S2LgeFpV/W1BKIdPV9KREREJK/cdM6UMWYB0B7wN8aEAS8ArgCWZc0EfgbuAUKAeGBUbhUrIiIikt9k5Wq+QTdZbgGP5lhFIiIiIgWI7mInIiIikg0KUyIiIiLZoDAlIiIikg0KUyIiIiLZoDAlIiIikg0KUyIiIiLZoDAlIiIikg1ZepxMrhzYmLPAsTw4lD8QmQfHyY/U96KrKPe/KPcdinb/1feiKy/6X8WyrIBrLXBYmMorxpit13uWTmGnvhfNvkPR7n9R7jsU7f6r70Wz7+D4/us0n4iIiEg2KEyJiIiIZENRCFMfOboAB1Lfi66i3P+i3Hco2v1X34suh/a/0M+ZEhEREclNRWFkSkRERCTXFNgwZYzpYow5YIwJMcb86xrL3Y0xC9OWbzLGVL1i2eS09gPGmLvztPAckIW+TzDG7DPG7DLGrDDGVLliWaoxZkfax+K8rTxnZKH/I40xZ6/o59grlo0wxhxK+xiRt5VnXxb6/vYV/T5ojLlwxbIC/d4bY2YbY84YY/ZcZ7kxxryb9rXZZYxpesWyAv2+Q5b6PySt37uNMX8YYxpdsSw0rX2HMWZr3lWdM7LQ9/bGmOgrvr+nXrHshj8z+V0W+j7xin7vSfs590tbVqDfdwBjTCVjzKq032l7jTFPXGMdx//sW5ZV4D4AZ+AwUB1wA3YC9f+2ziPAzLTPBwIL0z6vn7a+O1AtbT/Oju5TDve9A+CV9vnDl/ue9jrW0X3Ig/6PBN6/xrZ+wJG0f0umfV7S0X3Kyb7/bf1xwOxC9N63BZoCe66z/B5gKWCAlsCmwvC+30L/W1/uF9D1cv/TXocC/o7uQy72vT2w5Brtt/Qzkx8/btb3v63bA1hZWN73tD6UA5qmfe4DHLzG//kO/9kvqCNTLYAQy7KOWJaVBHwJ9PrbOr2AuWmffwN0NMaYtPYvLctKtCzrKBCStr+C4qZ9tyxrlWVZ8WkvNwIV87jG3JSV9/567gaWW5YVZVnWeWA50CWX6swNt9r3QcCCPKksD1iWtQaIusEqvYB5lt1GoIQxphwF/30Hbt5/y7L+SOsfFLKf+yy899eTnf8v8oVb7Huh+pkHsCwrwrKs7WmfxwD7gQp/W83hP/sFNUxVAE5c8TqMzF/c9HUsy0oBooFSWdw2P7vV+sdgT+yXeRhjthpjNhpjeudCfbktq/3vmzbc+40xptItbptfZbn+tFO71YCVVzQX9Pf+Zq739Sno7/s/8fefewv41RizzRjzgINqym2tjDE7jTFLjTEN0tqKzHtvjPHCHhS+vaK5UL3vxj5dpwmw6W+LHP6z75IbO5X8wRgzFAgG2l3RXMWyrHBjTHVgpTFmt2VZhx1TYa75EVhgWVaiMeZB7COUdzq4prw2EPjGsqzUK9qKwntf5BljOmAPU7df0Xx72ntfGlhujPkrbcSjsNiO/fs71hhzD/ADUMuxJeW5HsB6y7KuHMUqNO+7McYbe1Acb1nWRUfX83cFdWQqHKh0xeuKaW3XXMcY4wIUB85lcdv8LEv1G2M6Ac8BPS3LSrzcbllWeNq/R4DV2FN+QXLT/luWde6KPs8CmmV123zuVuofyN+G+wvBe38z1/v6FPT3PcuMMQ2xf8/3sizr3OX2K977M8D3FKypDTdlWdZFy7Ji0z7/GXA1xvhThN57bvwzX6Dfd2OMK/YgNd+yrO+usYrjf/bzejJZTnxgH1E7gv00xuVJhQ3+ts6jXD0B/au0zxtw9QT0IxSsCehZ6XsT7JMua/2tvSTgnva5P3CIgjcZMyv9L3fF532AjWmf+wFH074OJdM+93N0n3Ky72nr1cU+8dQUpvc+rfaqXH8ScjeunoS6uTC877fQ/8rY54C2/lt7McDnis//ALo4ui853Peyl7/fsQeG42nfB1n6mcnvHzfqe9ry4tjnVRUrhO+7AeYB79xgHYf/7BfI03yWZaUYYx4DfsF+tcZsy7L2GmNeArZalrUY+AT4zBgTgv2bbGDatnuNMV8B+4AU4FHr6lMh+VoW+/4m4A18bZ9zz3HLsnoC9YAPjTE27KOS0y3L2ueQjvxDWez/48aYntjf3yjsV/dhWVaUMeZlYEva7l6yrh4Sz9ey2Hewf69/aaX9b5KmwL/3xpgF2K/a8jfGhAEvAK4AlmXNBH7GflVPCBAPjEpbVqDf98uy0P+p2OeFfpD2c59i2R/8Wgb4Pq3NBfjCsqxled6BbMhC3/sBDxtjUoBLwMC07/9r/sw4oAv/WBb6DvY/Gn+1LCvuik0L/Puepg0wDNhtjNmR1vYs9j8e8s3Pvu6ALiIiIpINBXXOlIiIiEi+oDAlIiIikg0KUyIiIiLZoDAlIiIikg0KUyIiIiLZoDAlIiIikg0KUyIiIiLZoDAlIiIikg3/D3/OOSguuRC5AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "pts = pybamm.linspace(0, 2, 100)\n", + "\n", + "fig, ax = plt.subplots(figsize=(10,5))\n", + "ax.plot(pts.evaluate(), pybamm.maximum(pts,1).evaluate(), lw=2, label=\"exact\")\n", + "ax.plot(pts.evaluate(), pybamm.softplus(pts,1,5).evaluate(), \":\", lw=2, label=\"softplus (k=5)\")\n", + "ax.plot(pts.evaluate(), pybamm.softplus(pts,1,10).evaluate(), \":\", lw=2, label=\"softplus (k=10)\")\n", + "ax.plot(pts.evaluate(), pybamm.softplus(pts,1,100).evaluate(), \":\", lw=2, label=\"softplus (k=100)\")\n", + "ax.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Solving a model with the exact maximum, and smooth approximations, demonstrates a clear speed-up even for a very simple model" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Exact: 262.878 us\n", + "Smooth, k=5: 259.492 us\n", + "Smooth, k=10: 221.944 us\n", + "Smooth, k=100: 262.987 us\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "ed38eaf261704906bb260e3c3a4f353c", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "interactive(children=(FloatSlider(value=0.0, description='t', max=2.0, step=0.02), Output()), _dom_classes=('w…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "model_exact = pybamm.BaseModel()\n", + "model_exact.rhs = {x: pybamm.maximum(x, 1)}\n", + "model_exact.initial_conditions = {x: 0.5}\n", + "model_exact.variables = {\"x\": x, \"max(x,1)\": pybamm.maximum(x, 1)}\n", + "\n", + "model_smooth = pybamm.BaseModel()\n", + "k = pybamm.InputParameter(\"k\")\n", + "model_smooth.rhs = {x: pybamm.softplus(x, 1, k)}\n", + "model_smooth.initial_conditions = {x: 0.5}\n", + "model_smooth.variables = {\"x\": x, \"max(x,1)\": pybamm.softplus(x, 1, k)}\n", + "\n", + "solver = pybamm.CasadiSolver(mode=\"fast\")\n", + "\n", + "# Exact solution\n", + "timer = pybamm.Timer()\n", + "time = 0\n", + "for _ in range(100):\n", + " exact_sol = solver.solve(model_exact, [0, 2])\n", + " # Report integration time, which is the time spent actually doing the integration\n", + " time += exact_sol.integration_time\n", + "print(\"Exact:\", timer.format(time/100))\n", + "sols = [exact_sol]\n", + "\n", + "ks = [5, 10, 100]\n", + "for k in ks:\n", + " time = 0\n", + " for _ in range(100):\n", + " sol = solver.solve(model_smooth, [0, 2], inputs={\"k\": k})\n", + " time += sol.integration_time\n", + " print(f\"Smooth, k={k}:\", timer.format(time/100))\n", + " sols.append(sol)\n", + "\n", + "pybamm.dynamic_plot(sols, [\"x\", \"max(x,1)\"], labels=[\"exact\"] + [f\"smooth (k={k})\" for k in ks]);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Other smooth approximations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here are the other smooth approximations for the other non-smooth functions:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Smooth minimum (softminus):\t log(exp(-10.0 * x) + exp(-10.0 * y)) / -10.0\n", + "Smooth heaviside (sigmoid):\t (1.0 + tanh(10.0 * (y - x))) / 2.0\n", + "Smooth absolute value: \t\t x * (exp(10.0 * x) - exp(-10.0 * x)) / (exp(10.0 * x) + exp(-10.0 * x))\n" + ] + } + ], + "source": [ + "pybamm.settings.set_smoothing_parameters(10)\n", + "print(\"Smooth minimum (softminus):\\t {!s}\".format(pybamm.minimum(x,y)))\n", + "print(\"Smooth heaviside (sigmoid):\\t {!s}\".format(x < y))\n", + "print(\"Smooth absolute value: \\t\\t {!s}\".format(abs(x)))\n", + "pybamm.settings.set_smoothing_parameters(\"exact\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + }, + "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": 4 +} diff --git a/examples/scripts/cycling_ageing_yang.py b/examples/scripts/cycling_ageing_yang.py index e82d3f4e38..13c8ef22aa 100644 --- a/examples/scripts/cycling_ageing_yang.py +++ b/examples/scripts/cycling_ageing_yang.py @@ -1,51 +1,29 @@ import pybamm as pb pb.set_logging_level("INFO") -options = {"sei": "ec reaction limited", "sei porosity change": True} +options = { + "sei": "ec reaction limited", + "sei porosity change": True, + "thermal": "x-lumped", +} param = pb.ParameterValues(chemistry=pb.parameter_sets.Ramadass2004) +param.update( + { + "Separator density [kg.m-3]": 397, + "Separator specific heat capacity [J.kg-1.K-1]": 700, + "Separator thermal conductivity [W.m-1.K-1]": 0.16, + }, + check_already_exists=False, +) model = pb.lithium_ion.DFN(options) experiment = pb.Experiment( [ - "Charge at 1 C until 4.2 V", - "Hold at 4.2 V until C/10", - "Rest for 5 minutes", - "Discharge at 2 C until 2.8 V", + "Charge at 0.3 C until 4.2 V", "Rest for 5 minutes", - ] - * 2 - + [ - "Charge at 1 C until 4.2 V", - "Hold at 4.2 V until C/20", - "Rest for 30 minutes", - "Discharge at C/3 until 2.8 V", - "Charge at 1 C until 4.2 V", - "Hold at 4.2 V until C/20", - "Rest for 30 minutes", "Discharge at 1 C until 2.8 V", - "Charge at 1 C until 4.2 V", - "Hold at 4.2 V until C/20", - "Rest for 30 minutes", - "Discharge at 2 C until 2.8 V", - "Charge at 1 C until 4.2 V", - "Hold at 4.2 V until C/20", - "Rest for 30 minutes", - "Discharge at 3 C until 2.8 V", + "Rest for 5 minutes", ] + * 5 ) sim = pb.Simulation(model, experiment=experiment, parameter_values=param) -sim.solve(solver=pb.CasadiSolver(mode="safe", dt_max=120)) -sim.plot( - [ - "Current [A]", - "Total current density [A.m-2]", - "Terminal voltage [V]", - "Discharge capacity [A.h]", - "Electrolyte potential [V]", - "Electrolyte concentration [mol.m-3]", - "Total negative electrode sei thickness", - "Negative electrode porosity", - "X-averaged negative electrode porosity", - "Negative electrode sei interfacial current density [A.m-2]", - "X-averaged total negative electrode sei thickness [m]", - ] -) +sim.solve(solver=pb.CasadiSolver(mode="safe", dt_max=120)) \ No newline at end of file diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py index b726cbebdc..b3db032b7a 100644 --- a/pybamm/expression_tree/binary_operators.py +++ b/pybamm/expression_tree/binary_operators.py @@ -138,7 +138,23 @@ def format(self, left, right): def __str__(self): """ See :meth:`pybamm.Symbol.__str__()`. """ - return "{!s} {} {!s}".format(self.left, self.name, self.right) + # Possibly add brackets for clarity + if isinstance(self.left, pybamm.BinaryOperator) and not ( + (self.left.name == self.name) + or (self.left.name == "*" and self.name == "/") + ): + left_str = "({!s})".format(self.left) + else: + left_str = "{!s}".format(self.left) + if isinstance(self.right, pybamm.BinaryOperator) and not ( + (self.name == "*" and self.right.name == "*") + or (self.name == "+" and self.right.name == "+") + or (self.name == "*" and self.right.name == "/") + ): + right_str = "({!s})".format(self.right) + else: + right_str = "{!s}".format(self.right) + return "{} {} {}".format(left_str, self.name, right_str) def get_children_domains(self, ldomain, rdomain): "Combine domains from children in appropriate way" diff --git a/pybamm/models/submodels/interface/sei/ec_reaction_limited.py b/pybamm/models/submodels/interface/sei/ec_reaction_limited.py index 72495bb27b..ae3bdbbd4b 100644 --- a/pybamm/models/submodels/interface/sei/ec_reaction_limited.py +++ b/pybamm/models/submodels/interface/sei/ec_reaction_limited.py @@ -93,7 +93,8 @@ def set_algebraic(self, variables): + self.domain.lower() + " electrode interfacial current density" ] - except KeyError: + except KeyError as e: + print(e) j = variables[ "X-averaged " + self.domain.lower() diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py index dd4e0d0b20..2bcebeb17f 100644 --- a/pybamm/solvers/algebraic_solver.py +++ b/pybamm/solvers/algebraic_solver.py @@ -82,6 +82,8 @@ def _integrate(self, model, t_eval, inputs=None): y_alg = np.empty((len(y0_alg), len(t_eval))) + timer = pybamm.Timer() + integration_time = 0 for idx, t in enumerate(t_eval): def root_fun(y_alg): @@ -135,6 +137,7 @@ def jac_fn(y_alg): method = self.method[5:] if jac_fn is None: jac_fn = "2-point" + timer.reset() sol = optimize.least_squares( root_fun, y0_alg, @@ -144,6 +147,7 @@ def jac_fn(y_alg): bounds=model.bounds, **self.extra_options, ) + integration_time += timer.time() # Methods which use minimize are specified as either "minimize", which # uses the default method, or with "minimize__methodname" elif self.method.startswith("minimize"): @@ -170,6 +174,7 @@ def jac_norm(y): (lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1]) ] extra_options["bounds"] = bounds + timer.reset() sol = optimize.minimize( root_norm, y0_alg, @@ -178,7 +183,9 @@ def jac_norm(y): jac=jac_norm, **extra_options, ) + integration_time += timer.time() else: + timer.reset() sol = optimize.root( root_fun, y0_alg, @@ -187,6 +194,7 @@ def jac_norm(y): jac=jac_fn, options=self.extra_options, ) + integration_time += timer.time() if sol.success and np.all(abs(sol.fun) < self.tol): # update initial guess for the next iteration @@ -210,4 +218,6 @@ def jac_norm(y): y_diff = np.r_[[y0_diff] * len(t_eval)].T y_sol = np.r_[y_diff, y_alg] # Return solution object (no events, so pass None to t_event, y_event) - return pybamm.Solution(t_eval, y_sol, termination="success") + sol = pybamm.Solution(t_eval, y_sol, termination="success") + sol.integration_time = integration_time + return sol diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 981b53d2f2..99ce106c08 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -578,6 +578,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None): "initial conditions" ] = model.concatenated_initial_conditions set_up_time = timer.time() + timer.reset() # (Re-)calculate consistent initial conditions self._set_initial_conditions(model, ext_and_inputs, update_rhs=True) @@ -647,7 +648,6 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None): t_eval_dimensionless[end_index - 1] * model.timescale_eval, ) ) - timer.reset() new_solution = self._integrate( model, t_eval_dimensionless[start_index:end_index], ext_and_inputs ) @@ -671,13 +671,13 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None): model, t_eval_dimensionless[end_index], ext_and_inputs ) - # restore old y0 - model.y0 = old_y0 - # Assign times solution.set_up_time = set_up_time solution.solve_time = timer.time() + # restore old y0 + model.y0 = old_y0 + # Add model and inputs to solution solution.model = model solution.inputs = ext_and_inputs diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py index 0b484edfbe..af86ccbcd6 100644 --- a/pybamm/solvers/casadi_algebraic_solver.py +++ b/pybamm/solvers/casadi_algebraic_solver.py @@ -118,6 +118,8 @@ def _integrate(self, model, t_eval, inputs=None): "constraints": list(constraints[len_rhs:]), }, ) + timer = pybamm.Timer() + integration_time = 0 for idx, t in enumerate(t_eval): # Evaluate algebraic with new t and previous y0, if it's already close # enough then keep it @@ -137,7 +139,9 @@ def _integrate(self, model, t_eval, inputs=None): t_eval_inputs_sym = casadi.vertcat(t, symbolic_inputs) # Solve try: + timer.reset() y_alg_sol = roots(y0_alg, t_eval_inputs_sym) + integration_time += timer.time() success = True message = None # Check final output @@ -179,4 +183,6 @@ def _integrate(self, model, t_eval, inputs=None): y_diff = casadi.horzcat(*[y0_diff] * len(t_eval)) y_sol = casadi.vertcat(y_diff, y_alg) # Return solution object (no events, so pass None to t_event, y_event) - return pybamm.Solution(t_eval, y_sol, termination="success") + sol = pybamm.Solution(t_eval, y_sol, termination="success") + sol.integration_time = integration_time + return sol diff --git a/pybamm/solvers/casadi_solver.py b/pybamm/solvers/casadi_solver.py index 48b74bd8a7..2d657c0cc0 100644 --- a/pybamm/solvers/casadi_solver.py +++ b/pybamm/solvers/casadi_solver.py @@ -396,11 +396,15 @@ def _run_integrator(self, model, y0, inputs, t_eval): # Try solving if use_grid is True: # Call the integrator once, with the grid + timer = pybamm.Timer() sol = integrator( 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()]) - return pybamm.Solution(t_eval, y_sol) + sol = pybamm.Solution(t_eval, y_sol) + sol.integration_time = integration_time + return sol else: # Repeated calls to the integrator x = y0_diff @@ -411,19 +415,24 @@ def _run_integrator(self, model, y0, inputs, t_eval): t_min = t_eval[i] t_max = t_eval[i + 1] inputs_with_tlims = casadi.vertcat(inputs, t_min, t_max) + timer = pybamm.Timer() sol = integrator( x0=x, z0=z, p=inputs_with_tlims, **self.extra_options_call ) + integration_time = timer.time() x = sol["xf"] z = sol["zf"] y_diff = casadi.horzcat(y_diff, x) if not z.is_empty(): y_alg = casadi.horzcat(y_alg, z) if z.is_empty(): - return pybamm.Solution(t_eval, y_diff) + sol = pybamm.Solution(t_eval, y_diff) else: y_sol = casadi.vertcat(y_diff, y_alg) - return pybamm.Solution(t_eval, y_sol) + sol = pybamm.Solution(t_eval, y_sol) + + sol.integration_time = integration_time + return sol except RuntimeError as e: # If it doesn't work raise error raise pybamm.SolverError(e.args[0]) diff --git a/pybamm/solvers/dummy_solver.py b/pybamm/solvers/dummy_solver.py index 483bba8a77..e4c9f5ae02 100644 --- a/pybamm/solvers/dummy_solver.py +++ b/pybamm/solvers/dummy_solver.py @@ -33,4 +33,6 @@ def _integrate(self, model, t_eval, inputs=None): """ y_sol = np.zeros((1, t_eval.size)) - return pybamm.Solution(t_eval, y_sol, termination="final time") + sol = pybamm.Solution(t_eval, y_sol, termination="final time") + sol.integration_time = 0 + return sol diff --git a/pybamm/solvers/idaklu_solver.py b/pybamm/solvers/idaklu_solver.py index 6e7cc7fcc3..93a8a3731b 100644 --- a/pybamm/solvers/idaklu_solver.py +++ b/pybamm/solvers/idaklu_solver.py @@ -234,6 +234,7 @@ def rootfn(t, y): ids = np.concatenate((rhs_ids, alg_ids)) # solve + timer = pybamm.Timer() sol = idaklu.solve( t_eval, y0, @@ -251,6 +252,7 @@ def rootfn(t, y): atol, rtol, ) + integration_time = timer.time() t = sol.t number_of_timesteps = t.size @@ -266,12 +268,14 @@ def rootfn(t, y): elif sol.flag == 2: termination = "event" - return pybamm.Solution( + sol = pybamm.Solution( sol.t, np.transpose(y_out), t[-1], np.transpose(y_out[-1])[:, np.newaxis], termination, ) + sol.integration_time = integration_time + return sol else: raise pybamm.SolverError(sol.message) diff --git a/pybamm/solvers/jax_solver.py b/pybamm/solvers/jax_solver.py index 77595b0a7a..715a82f05f 100644 --- a/pybamm/solvers/jax_solver.py +++ b/pybamm/solvers/jax_solver.py @@ -188,10 +188,12 @@ def _integrate(self, model, t_eval, inputs=None): various diagnostic messages. """ + timer = pybamm.Timer() if model not in self._cached_solves: self._cached_solves[model] = self.create_solve(model, t_eval) y = self._cached_solves[model](inputs) + integration_time = timer.time() # note - the actual solve is not done until this line! y = onp.array(y) @@ -199,4 +201,6 @@ def _integrate(self, model, t_eval, inputs=None): termination = "final time" t_event = None y_event = onp.array(None) - return pybamm.Solution(t_eval, y, t_event, y_event, termination) + sol = pybamm.Solution(t_eval, y, t_event, y_event, termination) + sol.integration_time = integration_time + return sol diff --git a/pybamm/solvers/scikits_dae_solver.py b/pybamm/solvers/scikits_dae_solver.py index 6f3b56ec8f..df2272e14c 100644 --- a/pybamm/solvers/scikits_dae_solver.py +++ b/pybamm/solvers/scikits_dae_solver.py @@ -130,7 +130,9 @@ def jacfn(t, y, ydot, residuals, cj, J): # set up and solve dae_solver = scikits_odes.dae(self.method, eqsres, **extra_options) + timer = pybamm.Timer() sol = dae_solver.solve(t_eval, y0, ydot0) + integration_time = timer.time() # return solution, we need to tranpose y to match scipy's interface if sol.flag in [0, 2]: @@ -144,12 +146,14 @@ def jacfn(t, y, ydot, residuals, cj, J): t_root = None else: t_root = sol.roots.t - return pybamm.Solution( + sol = pybamm.Solution( sol.values.t, np.transpose(sol.values.y), t_root, np.transpose(sol.roots.y), termination, ) + sol.integration_time = integration_time + return sol else: raise pybamm.SolverError(sol.message) diff --git a/pybamm/solvers/scikits_ode_solver.py b/pybamm/solvers/scikits_ode_solver.py index 457e5b520c..0c4d6b9cd4 100644 --- a/pybamm/solvers/scikits_ode_solver.py +++ b/pybamm/solvers/scikits_ode_solver.py @@ -147,7 +147,9 @@ def jac_times_setupfn(t, y, fy, userdata): extra_options.update({"rootfn": rootfn, "nr_rootfns": len(events)}) ode_solver = scikits_odes.ode(self.method, eqsydot, **extra_options) + timer = pybamm.Timer() sol = ode_solver.solve(t_eval, y0) + integration_time = timer.time() # return solution, we need to tranpose y to match scipy's ivp interface if sol.flag in [0, 2]: @@ -161,12 +163,14 @@ def jac_times_setupfn(t, y, fy, userdata): t_root = None else: t_root = sol.roots.t - return pybamm.Solution( + sol = pybamm.Solution( sol.values.t, np.transpose(sol.values.y), t_root, np.transpose(sol.roots.y), termination, ) + sol.integration_time = integration_time + return sol else: raise pybamm.SolverError(sol.message) diff --git a/pybamm/solvers/scipy_solver.py b/pybamm/solvers/scipy_solver.py index 613ae72a51..41eb69838a 100644 --- a/pybamm/solvers/scipy_solver.py +++ b/pybamm/solvers/scipy_solver.py @@ -83,6 +83,7 @@ def event_fn(t, y): events = [event_wrapper(event) for event in model.terminate_events_eval] extra_options.update({"events": events}) + timer = pybamm.Timer() sol = it.solve_ivp( lambda t, y: model.rhs_eval(t, y, inputs), (t_eval[0], t_eval[-1]), @@ -92,6 +93,7 @@ def event_fn(t, y): dense_output=True, **extra_options ) + integration_time = timer.time() if sol.success: # Set the reason for termination @@ -107,6 +109,8 @@ def event_fn(t, y): termination = "final time" t_event = None y_event = np.array(None) - return pybamm.Solution(sol.t, sol.y, t_event, y_event, termination) + sol = pybamm.Solution(sol.t, sol.y, t_event, y_event, termination) + sol.integration_time = integration_time + return sol else: raise pybamm.SolverError(sol.message) diff --git a/pybamm/solvers/solution.py b/pybamm/solvers/solution.py index 5844606445..2c131397a4 100644 --- a/pybamm/solvers/solution.py +++ b/pybamm/solvers/solution.py @@ -53,12 +53,14 @@ def __init__( self._model = pybamm.BaseModel() self.set_up_time = None self.solve_time = None + self.integration_time = None self.has_symbolic_inputs = False else: self._inputs = copy.copy(copy_this.inputs) self._model = copy_this.model self.set_up_time = copy_this.set_up_time self.solve_time = copy_this.solve_time + self.integration_time = copy_this.integration_time self.has_symbolic_inputs = copy_this.has_symbolic_inputs # initiaize empty variables and data @@ -396,6 +398,7 @@ def append(self, solution, start_index=1, create_sub_solutions=False): self.inputs[name] = np.c_[inp, solution_inp[:, start_index:]] # Update solution time self.solve_time += solution.solve_time + self.integration_time += solution.integration_time # Update termination self._termination = solution.termination self._t_event = solution._t_event diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py index 7b409c66ce..b758816c67 100644 --- a/tests/unit/test_expression_tree/test_binary_operators.py +++ b/tests/unit/test_expression_tree/test_binary_operators.py @@ -143,12 +143,22 @@ def test_diff(self): self.assertEqual((a / a).diff(a).evaluate(y=y), 0) self.assertEqual((a / a).diff(b).evaluate(y=y), 0) - def test_addition_printing(self): + def test_printing(self): a = pybamm.Symbol("a") b = pybamm.Symbol("b") - summ = pybamm.Addition(a, b) - self.assertEqual(summ.name, "+") - self.assertEqual(str(summ), "a + b") + c = pybamm.Symbol("c") + d = pybamm.Symbol("d") + self.assertEqual(str(a + b), "a + b") + self.assertEqual(str(a + b + c + d), "a + b + c + d") + self.assertEqual(str((a + b) + (c + d)), "a + b + c + d") + self.assertEqual(str((a + b) * (c + d)), "(a + b) * (c + d)") + self.assertEqual(str(a * b * (c + d)), "a * b * (c + d)") + self.assertEqual(str((a * b) * (c + d)), "a * b * (c + d)") + self.assertEqual(str(a * (b * (c + d))), "a * b * (c + d)") + self.assertEqual(str((a + b) / (c + d)), "(a + b) / (c + d)") + self.assertEqual(str(a * b / (c + d)), "a * b / (c + d)") + self.assertEqual(str((a * b) / (c + d)), "a * b / (c + d)") + self.assertEqual(str(a * (b / (c + d))), "a * b / (c + d)") def test_id(self): a = pybamm.Scalar(4) @@ -310,13 +320,13 @@ def test_sigmoid(self): self.assertAlmostEqual(sigm.evaluate(y=np.array([2]))[0, 0], 1) self.assertEqual(sigm.evaluate(y=np.array([1])), 0.5) self.assertAlmostEqual(sigm.evaluate(y=np.array([0]))[0, 0], 0) - self.assertEqual(str(sigm), "1.0 + tanh(10.0 * y[0:1] - 1.0) / 2.0") + self.assertEqual(str(sigm), "(1.0 + tanh(10.0 * (y[0:1] - 1.0))) / 2.0") sigm = pybamm.sigmoid(b, a, 10) self.assertAlmostEqual(sigm.evaluate(y=np.array([2]))[0, 0], 0) self.assertEqual(sigm.evaluate(y=np.array([1])), 0.5) self.assertAlmostEqual(sigm.evaluate(y=np.array([0]))[0, 0], 1) - self.assertEqual(str(sigm), "1.0 + tanh(10.0 * 1.0 - y[0:1]) / 2.0") + self.assertEqual(str(sigm), "(1.0 + tanh(10.0 * (1.0 - y[0:1]))) / 2.0") def test_modulo(self): a = pybamm.StateVector(slice(0, 1))