diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index acf304f9..b6186c06 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -35,13 +35,12 @@ repos: - numpy - packaging - sympy + - optree - repo: https://github.com/codespell-project/codespell rev: v2.4.1 hooks: - id: codespell - args: ["-LHEP"] - exclude: ^docs/usage/intro.ipynb$ - repo: https://github.com/rbubley/mirrors-prettier rev: "v3.6.2" diff --git a/README.md b/README.md index b8cb42ea..77b8c42b 100644 --- a/README.md +++ b/README.md @@ -57,6 +57,12 @@ Vectors may be included in any of these data types: Each of these "backends" provides the same suite of properties and methods, through a common "compute" library. +### Integrations + +Optionally, the vector package provides integration with other libraries. Currently, this includes: + +- [PyTree integrations](https://vector.readthedocs.io/en/latest/src/pytree.html) using the [optree](https://github.com/metaopt/optree) package. + ### Geometric versus momentum Finally, vectors come in two flavors: @@ -134,6 +140,10 @@ journal = {Journal of Open Source Software} - [Interface for 3D momentum](https://vector.readthedocs.io/en/latest/src/momentum3d.html) - [Interface for 4D momentum](https://vector.readthedocs.io/en/latest/src/momentum4d.html) +### Integrations + +- [PyTree integration API](https://vector.readthedocs.io/en/latest/src/pytree_api.html) + ### More ways to learn - [Papers and talks](https://vector.readthedocs.io/en/latest/src/talks.html) diff --git a/docs/index.md b/docs/index.md index 07614dc8..ca07e1b9 100644 --- a/docs/index.md +++ b/docs/index.md @@ -57,6 +57,12 @@ Vectors may be included in any of these data types: Each of these "backends" provides the same suite of properties and methods, through a common "compute" library. +### Integrations + +Optionally, the vector package provides integration with other libraries. Currently, this includes: + +- [PyTree integrations](https://vector.readthedocs.io/en/latest/src/pytree.html) using the [optree](https://github.com/metaopt/optree) package. + ### Geometric versus momentum Finally, vectors come in two flavors: @@ -133,6 +139,7 @@ src/numpy.ipynb src/awkward.ipynb src/numba.ipynb src/sympy.ipynb +src/pytree.ipynb ``` ```{toctree} @@ -156,6 +163,12 @@ src/momentum3d.md src/momentum4d.md ``` +```{toctree} +:maxdepth: 1 +:caption: Integrations +src/pytree_api.md +``` + ```{toctree} :maxdepth: 1 :caption: More ways to learn diff --git a/docs/src/pytree.ipynb b/docs/src/pytree.ipynb new file mode 100644 index 00000000..4eaf6965 --- /dev/null +++ b/docs/src/pytree.ipynb @@ -0,0 +1,292 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "3e852d83", + "metadata": {}, + "source": [ + "# PyTree integrations\n", + "\n", + "PyTrees are a [powerful mechanism](https://blog.scientific-python.org/pytrees/) for working with\n", + "nested data structures, while allowing algorithms like finite-differences, minimization, and integration routines\n", + "to run on flattened 1D arrays of the the same data. To use PyTrees with vector objects, you need to install\n", + "the [optree](https://github.com/metaopt/optree) package, and then register vector with PyTrees, for example:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "6d3ecebc", + "metadata": {}, + "outputs": [], + "source": [ + "import vector\n", + "\n", + "pytree = vector.register_pytree()\n", + "\n", + "state = {\n", + " \"position\": vector.obj(x=1, y=2, z=3, t=0),\n", + " \"momentum\": vector.obj(x=0, y=10, z=0, t=14),\n", + "}\n", + "flat_state, treedef = pytree.flatten(state)" + ] + }, + { + "cell_type": "markdown", + "id": "92744f6a", + "metadata": {}, + "source": [ + "`flat_state` is now a 1D array of length 8, and `treedef` contains the information needed to reconstruct the original structure:\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "aa4b7219", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0, 10, 0, 14, 1, 2, 3, 0]\n", + "PyTreeSpec({'momentum': CustomTreeNode(VectorObject4D[(, , , )], [(*, *), (*,), (*,)]), 'position': CustomTreeNode(VectorObject4D[(, , , )], [(*, *), (*,), (*,)])}, namespace='vector')\n" + ] + } + ], + "source": [ + "print(flat_state)\n", + "print(treedef)" + ] + }, + { + "cell_type": "markdown", + "id": "7577e134", + "metadata": {}, + "source": [ + "The original structure can be reconstructed with:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "419cd827", + "metadata": {}, + "outputs": [], + "source": [ + "reconstructed_state = pytree.unflatten(treedef, flat_state)" + ] + }, + { + "cell_type": "markdown", + "id": "86fa7add", + "metadata": {}, + "source": [ + "## Example: projectile motion with air resistance\n", + "\n", + "In the following example, we use `scipy.integrate.solve_ivp` to solve the equations of motion for a projectile under the influence of gravity and air resistance.\n", + "\n", + "To start, we wrap the solver, using PyTrees to flatten the state dictionary into a 1D array for the integrator, and then unflatten the result back into the original structure" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "c53ba724", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.integrate import solve_ivp\n", + "\n", + "\n", + "def wrapped_solve(fun, t_span, y0, t_eval):\n", + " flat_y0, treedef = pytree.flatten(y0)\n", + "\n", + " def flat_fun(t, flat_y):\n", + " state = pytree.unflatten(treedef, flat_y)\n", + " dstate_dt = fun(t, state)\n", + " flat_dstate_dt, _ = pytree.flatten(dstate_dt)\n", + " return flat_dstate_dt\n", + "\n", + " flat_solution = solve_ivp(\n", + " fun=flat_fun,\n", + " t_span=t_span,\n", + " y0=flat_y0,\n", + " t_eval=t_eval,\n", + " )\n", + " return pytree.unflatten(treedef, flat_solution.y)" + ] + }, + { + "cell_type": "markdown", + "id": "43585fe2", + "metadata": {}, + "source": [ + "Now, we set up the physical constants and initial conditions for the projectile motion problem, using the [hepunits](https://github.com/scikit-hep/hepunits) library to help us consistently track units throughout the calculation." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "bd5adaf6", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import hepunits as u\n", + "\n", + "gravitational_acceleration = vector.obj(x=0.0, y=-9.81) * (u.m / u.s**2)\n", + "air_density = 1.25 * u.kg / u.meter3\n", + "drag_coefficient = 0.47 # for a sphere\n", + "object_radius = 10 * u.cm\n", + "object_area = np.pi * object_radius**2\n", + "object_mass = 1.0 * u.kg\n", + "\n", + "initial_position = vector.obj(\n", + " x=0 * u.m,\n", + " y=0 * u.m,\n", + ")\n", + "initial_momentum = (\n", + " vector.obj(\n", + " x=1 * u.m / u.s,\n", + " y=20 * u.m / u.s,\n", + " )\n", + " * object_mass\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "2f4ef621", + "metadata": {}, + "source": [ + "The `force` function computes the total force on the object, including both gravity and air resistance. The `tangent` function then uses this to compute the time derivatives of position and momentum." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "69d7a592", + "metadata": {}, + "outputs": [], + "source": [ + "def force(position: vector.Vector2D, momentum: vector.Vector2D):\n", + " gravitational_force = object_mass * gravitational_acceleration\n", + " speed = momentum.rho / object_mass\n", + " drag_magnitude = 0.5 * air_density * speed**2 * drag_coefficient * object_area\n", + " drag_force = -drag_magnitude * momentum.unit()\n", + " return gravitational_force + drag_force\n", + "\n", + "\n", + "def tangent(t: float, state):\n", + " position, momentum = state\n", + " dposition_dt = momentum / object_mass\n", + " dmomentum_dt = force(position, momentum)\n", + " return (dposition_dt, dmomentum_dt)" + ] + }, + { + "cell_type": "markdown", + "id": "6cf2d93d", + "metadata": {}, + "source": [ + "Finally, we solve the equations of motion using our wrapped `solve_ivp`, and plot the resulting trajectory alongside an analytic solution for projectile motion without air resistance." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "cc43f146", + "metadata": {}, + "outputs": [], + "source": [ + "time = np.linspace(0, 4, 100) * u.s\n", + "\n", + "position, momentum = wrapped_solve(\n", + " fun=tangent,\n", + " t_span=(time[0], time[-1]),\n", + " y0=(initial_position, initial_momentum),\n", + " t_eval=time,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "dc7f8595", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAHHCAYAAABHp6kXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcptJREFUeJzt3QdYU1cbB/A/GwQEERmKe+HAvbdV695aR2u1U63WWru07Ve7bW1rh7V2uTqdVVv33nsP3CKigrhYskm+5z0RChYQFbjJzf/3PDE3yU14b0K8L+e85xwbo9FoBBEREZGFs9U6ACIiIqL8wKSGiIiIdIFJDREREekCkxoiIiLSBSY1REREpAtMaoiIiEgXmNQQERGRLjCpISIiIl1gUkNERES6wKSGiO6pXLlyGDZsWKH/3NmzZ8PGxgYXLlzIuK9Nmzbqoid6PCYiLTCpIbJQ6Sf89IuzszOqVKmC0aNH4+rVq7AkH3/8MZYsWQI9Cw4OxrvvvpslQSOi/GXDtZ+ILDepeeqpp/D++++jfPnySExMxLZt2/Drr7+ibNmyOHbsGIoUKZIvPyspKQm2trZwcHBAQXBzc0O/fv3UMWWWlpaGlJQUODk5qcRNpLdobNq0CZZk4cKF6N+/PzZu3PifVpnk5GR17ejoqFF0RPpgr3UARPRwOnfujAYNGqjtZ599FsWLF8eUKVOwdOlSDBo0KNvn3L59G66urnn+GZJUaMHOzk5d9I7JDFH+YPcTkc488sgj6jokJERdSy2MtIScO3cOXbp0gbu7Ox5//PGM5OaVV15B6dKlVeJStWpVfP7557i7ATe7mpqoqCiMHTs247mVKlXCp59+CoPBkGU/uf31118jKChIdZGVKFECnTp1wr59+9Tj0gIjccyZMyejKy39Z2VXU5NTS9LEiRNVDBKLxPT666+r++9FWk1q1qyJI0eOoHXr1qp1S15HWlbE5s2b0bhxY7i4uKj3Z926df95jYMHD6rksmjRouq9bteuHXbt2pXxuByHtNKItm3bZhxnemtTdjU1kZGReOaZZ+Dr66vet9q1a6v3KDN5X+R15DP78ccfUbFiRXX8DRs2xN69e+957ER6w5YaIp2R5EVIi0261NRUdOzYES1atFAnQDlxS+LSo0cP1R0iJ886depg9erVeO2113D58mV8+eWXOf6M+Ph4lQDIfsOHD0eZMmWwY8cOTJgwAeHh4fjqq68y9pXXlpO6nPSlJUli2bp1qzrpSwuTdJfJ/Y0aNcLzzz+vniMn57ySpEmOQ7re5PnVqlXD0aNHVfynT5/OU63OrVu30K1bNwwcOFAlH9OnT1fbv//+u0rcRowYgcGDB+Ozzz5T3WRhYWEqORTHjx9Hy5YtVUIjiZR00f3www8qSUlPiFq1aoUxY8bgm2++wZtvvqliFOnXd0tISFDPP3v2rKqRku7FBQsWqGRPksmXXnopy/5//PEHYmNj1WchSc7kyZPRp08fnD9/vsC6DInMktTUEJHlmTVrljSnGNetW2e8du2aMSwszDh37lxj8eLFjS4uLsZLly6p/YYOHar2Gz9+fJbnL1myRN3/4YcfZrm/X79+RhsbG+PZs2cz7itbtqx6nXQffPCB0dXV1Xj69Oksz5WfYWdnZ7x48aK6vWHDBvUzxowZ85/4DQZDxra8VubXv/sYQ0JCMu5r3bq1uqT79ddfjba2tsatW7dmee7333+vnrt9+/Zc30d5Ldnvjz/+yLjv5MmT6j553V27dmXcv3r1anW/xJWuV69eRkdHR+O5c+cy7rty5YrR3d3d2KpVq4z7FixYoJ67cePGbGPIfExfffWV2ve3337LuC85OdnYtGlTo5ubmzEmJkbdJ++L7Cef+c2bNzP2Xbp0qbr/n3/+yfXYifSG3U9EFq59+/aqS0e6XKR1Qbo/Fi9ejFKlSmXZb+TIkVlur1ixQtWrSAtCZtIdJa04K1euzPFnSquBtE4UK1YM169fz7hILFLcu2XLFrXfokWLVMuBdA3dLb3w92FJLNLiERgYmCWW9G44aYm6F3nP5L1LJ91Mnp6e6nWlpSVd+ra0gAg51jVr1qBXr16oUKFCxn7+/v6qZUdaj2JiYu77mOSz8fPzy1ITJS0u8lnFxcWpFqDMBgwYoD6LdPLZZI6TyFqw+4nIwk2bNk0N5ba3t1f1F3JClpFKmcljAQEBWe4LDQ1FyZIlM7pR0qV3icjjOTlz5oyqQZFkKjtSD5LeFSY/w8vLCwVFYjlx4sQ9Y8mNvDd3J1keHh4qUbz7vvTuKnHt2jXVFSfv+d3kfZSuMemqqlGjxn0dk7z3lStX/s/nmNNnI91/maUnOOlxElkLJjVEFk5qUdJHP+VEikfvPkE+DDlZd+jQQdWQZEeSrMIisUgRsoz4ys7diUl2chphldP95jYThqXESVTQmNQQWSmZy0ZG8kiBaebWmpMnT2Y8nhMp5JVuEOluyo3sJ8XHN2/ezLW15mG6ouRnHD58WI04yq8urbyS1iEpuj516tR/HpP3URLJ9KTqfmKT915awiRhy5yM5uWzIbJmrKkhslIyvFtqQr799tss98uoITkBy2ilnDz22GPYuXOnSljuJqNzZIST6Nu3r2oteO+993JtRZA5c+R5D0JikVFYP/30U7ajiGS4eEG2kDz66KNqTqDMw85lRmcZkSSjzWRUlEifFygvxymfTUREBObNm5dxn7ynU6dOVfU/MvKMiP6LLTVEVqp79+5qzpS33npLnZBlHhQpepUTtAxjzm1YtQz7/vvvv9UwaBlmXL9+fZU8yFBqmd9FXs/b21u9/pAhQ9RQZql9kflppPVBhnTLYzJcWcjzpdVIupCkBkeGMGcu0M2NvP78+fPVsGspCm7evLlK1qRVQ+6XxOte3XMP48MPP8TatWtVAvPCCy+o+iUZ0i1z5MjQ6nQyZF6SIJnLJzo6WnUJSjGzj4/Pf15ThqbLa8h7u3//fjVPkLyv27dvV8Pl766DIiITJjVEVkq6NSQxeeedd1SLwKxZs9TJU+ZikRFQuZEuFxmBI2s2yeijX375RbVISC2NtMqkF9QKed1atWphxowZKhmSxyTJaNasWcY+kszIifztt99WrStDhw7Nc1IjxyFz0UgLk8QhI78kPhmNJPO5FHR9jxQBS5Imc/RMmjRJJW0S+2+//ZblGGQ00/fff6/2kbl7JPGSJCy7pEYm+pOJ+caPH68m3JMRVFKMLO+lFguLElkKrv1ERPckdSEyed/PP/+sdShERDliTQ0R5UoWlLxx44bqTiIiMmfsfiKiHEk9yty5c1WXkIwuIiIyZ+x+IqIcSTGvrD8ksxHLmkVEROaMSQ0RERHpAmtqiIiISBeY1BAREZEuWFWhsMwfceXKFTVxVWFPp05EREQPRiplZEkXmZwzt3XsrCqpkYQmL4vbERERkfmRVe8DAgJyfNyqkpr0qcXlTUlfj4WIiIjMm8yqLY0S91oixKqSmvQuJ0lomNQQERFZlnuVjrBQmIiIiHSBSQ0RERHpApMaIiIi0gUmNURERKQLTGqIiIhIF5jUEBERkS4wqSEiIiJdYFJDREREusCkhoiIiHSBSQ0RERHpApMaIiIi0gUmNURERKQLVrWgJRHpmNEoq92ZtlOTgLhI07bcZ2MLOBQBHN0AO/63R6RX/HYTkWW4dQE48Q8QEw7EXb1ziQTirwPJt4EOHwCNnzftG34YmNEh+9exdwZavQa0etV0+/Z1YOsUwN0XKFYO8KoIeFUAHIsU3rERUb5gUkNE5sFgAG6eByIOAxFHgfAjQP2hQPWe/yY1a97O+fnJcf9uS8uMJC/SeiOMaYAh1bSdmgjY2v27761QYNe0/75e0VKAd2Wg3lCgZp/8OUYi0n9SM2nSJPz11184efIkXFxc0KxZM3z66aeoWrVqxj6JiYl45ZVXMHfuXCQlJaFjx4747rvv4Ovrq2nsRPQQYiOAA78AoTuAS3uzJiaiZN1/kxrvqkDNvoBHAODqA7j5Am4+gKs34OQOuHj9+7yABsDbV7O+Vmqy6fWTYk37pytSDGj2oqkF6FYIcOMckBgFxFw2Xap0+nffyJPAyteAMs2Ass2AgIZs0SEyIzZGY/qfMtrp1KkTBg4ciIYNGyI1NRVvvvkmjh07huDgYLi6uqp9Ro4cieXLl2P27Nnw8PDA6NGjYWtri+3bt+f558TExKjnRkdHo2jRogV4RET0H/JfTcQR+W8H8K9luu9mCPBNnX/3sXcBfGuYHvcLAsq2AEpUKfxY42+akpvIYKBsc8C7kun+Q38CS0b8u5+tPVCyHlCpPVC5A+BfB7Dl+Aui/JbX87dZJDV3u3btGnx8fLB582a0atVKHUSJEiXwxx9/oF+/fmofadWpVq0adu7ciSZNmuTpdZnUEBUyKdg9t8FUC3N2nakOJrAbMPB30+Py38+ysYBvTaBMU8CnWtauIXMTfRk4s9rUsiQXacnJ7LFfgeo9tIqOSLfyev42i+6nu0nQwsvL1Jy8f/9+pKSkoH379hn7BAYGokyZMrkmNdJNJZfMbwoRFYJTq4DjfwGnVgJJmb53Dq6Ao6n1NWNkUvevYTE8SgENnjZdJCGLCgXObwbOrAEubAMqtP53313fm7rUpB6nYjvAwVnLyImsgtklNQaDAWPHjkXz5s1Rs2ZNdV9ERAQcHR3h6emZZV+pp5HHcqvVee+99wo8ZiK6y5bPgMv7TNvu/qa6mKqdTa0x9k7QBUnIZLRUfbkMBQxpWVuZDv1u6m47thBwKgrU6A3UexIoVf/foedEpO+kZtSoUaqeZtu2bQ/9WhMmTMC4ceOytNSULl36oV+XiDLVnhz+Ezi2CBiyGHD2MN3f4ClTEW2NXkBAI+uoM7m726zrFCB4CXB8samb6sAc06VEoKmlp/FwrSIl0i2zSmqk+HfZsmXYsmULAgICMu738/NDcnIyoqKisrTWXL16VT2WEycnJ3UhonwmQ653/wAcXWAaIi0ksZGTtaj7hKbhmYXSDU0XmT8ndDtw8DcgeClw7aSpHodJDZE+kxqpVX7xxRexePFibNq0CeXLl8/yeP369eHg4ID169ejb9++6r5Tp07h4sWLaNq0qUZRE1nhPDKnlptqRUIztaTKKKX6w0zdK/Rf0kpVvqXp0mUycOwvwKf6v4/L3DxLRpkmDgzszhmPiR6Cvbl0OcnIpqVLl8Ld3T2jTkYqnWXeGrl+5plnVFeSFA9L5bMkQZLQ5HXkExE9JCmKnf8kYDQANnZAte5A4xFAmSasEckr6Z6TrrnMdv8IXNxhuhQrD7R4Gag9CLB31CpKIotlFkO6bXL4D3HWrFkYNmxYlsn3/vzzzyyT7+XW/XQ3Dukmug8pCcCF7UDlf0cdYtnLgLMn0PBZ00ggengy6d++GcDeGUDCzX9nM242xlRYzMn9iGDR89QUFCY1RHmQlgIc/BXYPNk0r8zofUDxilpHpX9JccD+2cCOqUDcnVGdMrpq9H52SZHVi7HkeWqISAMyJFmKfTd+bFouQBQNMI3cYVJT8JzcgGajgUbPmYaDb/sKqNkva0IjdU3WMJKM6AExqSEi06y/q98GIo+bbruWAFq+aqr/0Mu8MpZC3m8ZRVZ3iKnVLF3IVmD5K0CH94EqHVnHRJQNJjVE1i4xGpg/DEiKNhWyNn8JaDTc1HJA2rFzMF3Sbf8auH4K+HMAUPERoOMkwCdQywiJzA5raoistX4jc9KyazpwKxRo/TpQJNNq12Q+EqKAbVNMn1VasmkEWsNngDYT+JmR7sWwUPi/mNSQ1ZOvu0yYt/pNoNd008rSZFlkXps1/wNOLjPddikGdPkcCDIt9ktkzedvVpwRWYvrZ4FfegJ/PQfcvgbsm6V1RPQgvCqYVjl/cqlpEr+EW6bV0ImINTVEuicnvG1fAlu/MHVb2DsDrV4zzYNClqtCG2D4FtPaUkH9/73/6nFT4uPgomV0RJpgUkOkZ1cOAkteACKDTbcrtTd1VXhlXYqELJQUEtd67N/biTHAb30BR1dT92LpRlpGR1To2P1EpGe3LpgSmiLeQL9ZwOMLmdDomcwvJHVTN84CMzuaam9S7iw4SmQFWChMpMflDTJ3PcgClNI94Vpcy6iosEiNzaoJwOE/Tbe9qwC9vgcC6msdGdEDY6EwkTXOCLx1CjC1ARB37d/7m4xgQmNNZDRU7++BQXMBN1/g+mlgRntg0yemGYmJdIxJDZEeRF82jWxa/x4Qcwk4/IfWEZHWqnYGXthlaqWTldWvneQsxKR7LBQmsnTBS4G/xwCJUYCDK9D5U6DuE1pHReZAJuXr+zNQtQtQse2/SU2qjIJz1Do6onzHpIbIkodqS+3Evhmm2yXrAn1ncPFJ+q+aff7dljLKRU8DTh5Al88AxyJaRkaUr5jUEFmqzZP/TWhavAy0fSvrWkFE2Qk/DJxcbuqSCj8EPPYLE2HSDdbUEFkqWXiydGPg8UVA+3eZ0FDelKxjmo1YVmK/egz4obWpC5NIB5jUEFnS6KbDc03dB8K5KPD0aqBye60jI0tTvhUwfCtQphmQHAvMfxJY9SaQlqJ1ZEQPhUkNkaXMPfJ7P2DxcNMqzek4moUeVFF/YOg/phY/sWsasOhZraMieihMaojMXeQJ4Me2wLkNgEMRwN1P64hIL+zsgQ7vAwP/AFy8gMbDtY6I6KGwUJjInElB51/PA8lxgEcZYNAfgF+Q1lGR3gR2NXVJObn/e1/sVcDdV8uoiO4bW2qIzJHUzWz5HJg72JTQlGsJPL+JCQ0VnMwJzdVgYGp9YMNHnIWYLAqTGiJzJItQbvzYtN3oeWDIYi51QIXn7FpTAfGWycBfz3JRTLIY7H4iMke+NYCuXwCGVKDRc1pHQ9ZGiodlZfd/xgDHFpmW4ZC6GybWZOa4SjeRubgZYpoQjROhkbk4vxmYNwRIigaKlQceXwh4V9I6KrJCMVylm8iCXD4A/NzeNGw7/qbW0RCZVGgNPLMG8CwD3AoxrfYto/GIzBSTGiKtnV0PzO4GxF8HHN2AtGStIyL6l08g8Ox6oFR9wL8O4MWWRDJfrKkh0tKRBcCSEabamQptgAG/ZR2FQmQO3HyAoctM3aPpq3tL5QInfyQzw5YaIq3s/M40skQSmpp9gcELmNCQ+ZLVvJ3c/k1oZIX4rVP+XbaDyAywpYZIC3t+AlZPMG03HgF0nATYms/fGIkpaYhPTkNCShoSklORmGKAna0N7G1t1LWjvS08XBzg5mQPG/61bn0ubAV231mu4/Y14NGPzOr3l6wXkxoiLVTvCez+AagzCGgxrtCb8SVpOX/tNs5fj8O5yNsIuR6HiJhERMYm4VpMEmKTUvP0OpLkeBZxgJerIwKKFUEZryIIKOaC8t6uqOZfFP4ezkx69EhmH5ZEZs1bwK7vTBNEdvsKsLXTOjKychzSTVRY7q5BSEkAHFwK4ccacSYyDgdCb+HwpSgcCovG6auxSDPc+6vvZG8LF0c7dS27y3NS0wxITDUgOfXeM81Ka041f3fUDvBEo/JeaFDWCx5FHPLpyEhzh/4Elr5gqrUJ6g/0mg7Y8fMl7c7fTGqICkNaCrB4BFDxEaDu4wX+4yJjE7Hl9HVsP3sd285ex7XYpGwTjoolXFGxhBsqlHBDSU9n+Lg7w6eoE7zdnFTXknQ15dbaExWfglvxyer1w27FI+xmAsJuxuNsZBzOXYtD6l2Jk+R0gX5F0aJScTwS6IsG5YrBwY7dFhbt+GLT6t5SGxbYDeg3E7B30joq0hkmNdlgUkOaSE0CFgwDTq0A7J2Blw4XyErbkkysPh6BlccicODirSz1m84OtqhbuhjqlPFE7QAP1C7tCb+iBds1lJSahjNX4xAcHoODF29h9/mbOH/9dpZ93J3t0bpKCXQJ8scjgT5wdmD3hUU6tQqY/6RpOoInl5rmtyHKR0xqssGkhjRJaGRG1jOrTQmNDNmu3CHfXj46IQX/HL6Chfsv4VBYVJbHagV4oGVlbzSv5I16ZYqZRcIgLUiS3Gw6dQ0bT0Xi5u1/5+Rxd7JHx5p+6FWnFJpWLJ5rKxGZoXMbgZgrhdISSdYnhknNfzGpIe0SGhdg8FzTXDQPSb6yu87fxJ97LqqWmaQ7tS2SAzQuXxydavrh0Rq+8Pco+HqdhyH1OZKIrQmOwLLD4bgclZDxWGkvFwxuVBaPNQhAcTd2ZVikuGuAo6tpKDjRQ2JSkw0mNWTJCY3UsPx96Apm7biAE+ExGfdX9XVH/wYB6FmnFEq4W2YCYDAYsS/0FpYeuqxanmISTaOvHO1s0SXID8+2rICapTy0DpPyKvYqMKcbULQkMGhuoRTEk74xqckGkxoqNAd/N40KyYeERrqY5uy4oC437nTXuDjYoXe9UhjYsDSCSnnoati0JG+S2Py2KxSHL0Vn3C9daSNbV1RdU3o6Xl0K2wv80hNIuQ1UaAsM+pOJDT0UJjXZYFJDhUa+Vhs+BMq3fOCEJio+GTO3hWDW9gsZ88aU8nTBk03LYmDDMlYxNPropWjM2HYe/xwJzxiCLkXOr3SoopIcJjdmLHQn8FtfU2Ijo/4GSmLjrHVUZKGY1GSDSQ0VKEOa6ZK+Ns4DiktKxU9bzuPnredxOzlN3VfF1w2j2lZC1yB/2FvhEGgZ2SXvx9y9YRk1RI3Le+H1TlVRv6yX1uFRTkJ3AL/1MyU2VToDA37lPDb0QJjUZINJDRUYgwH4+0Ug7qrpP+4HaGpPSTOok/bX607jepypm0lm5R3zSCV0rOEHW44GwvW4JEzfdA6/7grNmPyvfTUfvNmlmpprh8xQyFbg935AaqJpjbM+P3HmYbpvTGqywaSGCoR8hVa+Duz5EbCxBYYsue95OjacvIoPl53ImMelXPEieK1joCqSZRfLf12JSsA3689gwf5LqlvKwc4GTzUvjxcfqQR3Z7YEmJ3Ta4C5gwGPUsDTawB3X60jIgvDpCYbTGoo38nXZ927wPav5OsE9P4BqD3gvrpV3vsnGOtOXFW3i7s64qX2lTGoURnOtJsHMmvxh8uCsfHUNXVbZkIe3zkQfeuVYjJobk6vBvxqAUX9tY6ELBCTmmwwqaF8t+1LU1Ijun0JNHg6z7Pt/rD5PKZtPKtqRGRhyKdbsKXhQW08GYkPlgVntHQ1r1QcH/cOQtnirlqHRjm5cQ4oXlHrKMhCMKnJBpMaylcHfgX+Hm3afvRDoNmLeXqaTDj3+sLDOH01Tt1uUsEL7/esiSq+7gUZre5Jjc3M7SH4at1pJKYY1CKcL3eogmdblLfK4mqzdngusHQU0HES0Ph5raMhHZ2/+U0nehAJt4DVb5m2m4/NU0Ij869MWnECfb7brhIabzdHfD2wDv58rgkTmnzgaG+LEa0rYvXYVqqlRlrAPll5En2m71ALbJIZibpoWgBTatGO/aV1NKQjbKkhelCXDwDHFplaae5RvyGtM+PmHcroHulVpyTe6V4DXq4PN/ybsif/rcl6WB8uP6EmL5RWGxkhNaRJWY4iMwdy2lnxGrBXRkI5AI8vACq21ToqMmPsfsoGkxrKl6Hbtnlv4JSROdM3ncWX686obd+iTvioVxDaV+foj8IQEZ2I1xYextYz19VtmbDv8/614VuUk8BpTuZ0Wvg0ELwEcHQDhi0DStbVOioyU+x+IspvUWHADy2BS/vztPulW/EY9OMufL7mtEpoutXyx5qxrZnQFCI/D2f88nQjvN+zBpwdbFVy0+Xrrdh6xjRaijQkc9X0+REo3wpIjjNN0ifFw0QPgUkNUV4kRAG/9weuHjPVAdyjgXPN8Qh0/nor9ly4CVdHO3zRvzamDqprFUsbmBsZ2v1k03JY9mJLNZmhrJ/15Mw9mLLmVMbSC6QReydgwO+Af20g/jpw6A+tIyILx+4nontJTQZ+7wuEbAHc/YFn1wEeAdnvmmbAZ6tP4Yct59XtOqU9VTEwhxabBynWfn9ZMP7YfTFj5NnUQfUsdnVz3YiNAI4uAJqOvmd9GlmnGNbU/BeTGrpv8vVYPAI4MtfU7//USsC/Vra7RsYkYvQfB1XrjHimRXk1ERwn0TM/Sw9dxpt/HVVra/l7OOOHIfVRK8BT67Aoc72NTGZ5H/VrpG+sqSHKD5smmRIaGzug/5wcE5r9oTfR5ZttKqFxc7LH9Mfr4X/dqjOhMVM965TC3y+2QMUSrgiPTkT/73di8cFLWodFIjkemP8ksG6i1pGQBeL/uEQ5OfEPsPnTf2cLrtw+293m7wvDwB93qcUWA/3c8c+LLdA5iFPBm7uKJdyweFRztAv0UXPavDzvMD5aHsw6G61JN+/JZcCOb4C9P2sdDVkYJjVEOan4CBDYDWj5ClB/aLb1MzI1/+sLjyAlzYjONf2waGQzlPdm/YylKOrsgJ+ebIDRbSup2z9tDcGI3/YjPjlV69CsV9VOQNu3TdsrXgfObdA6IrIgrKkhute8NOKuvv3YxBSM+uMgtpw2DQ0e274yxjxSmRO7WbBlR65g3PzDarmFWgEe+HloA/i4cz4bTchpaclI4PCfgJOHqTi/RBWtoyINsaaG6EEkRgN7fvp3yLYkM3clNOHRCaoGQxIaFwc7fPd4PYxtX4UJjYXrVqsk/nyuMYoVccCRS9HoPW0HTl+N1Tos6yQjoLp/DZRuAiRFA388BsSbCvCJcsOkhihdWiqw4ClgxavAmjvN33c5GRGjTnYnI2LVMOD5w5uiC+tndKN+WS8sfqG56kK8HJWAftN3YN+d0WykxRw2vwEeZYBbIcBfXPiSLCip2bJlC7p3746SJUuqybKWLFmS5fFhw4ap+zNfOnXqpFm8pENr3gLOrQccigBB/f/z8Paz19F/+k5ExCSiko8b/hrZDEEBHpqESgWnnLer+mzrly2GmMRUPDFjNzaeitQ6LOvkVgIYPBfwqgC0fkPraMgCmE1Sc/v2bdSuXRvTpk3LcR9JYsLDwzMuf/75Z6HGSDp24Fdg9/em7d4/ACXrZHl4+ZFwDJu1B7FJqWhc3guLRjRDaa8i2sRKBa6YqyN+e6Yx2lQtgcQUA56bs0/NbUMa8K0BjNoLlG6odSRkAexhJjp37qwuuXFycoKfn1+hxURWImwvsHycabvtW0D1HlkenrvnIt5cfBQy0rdrLX9Meaw2nOzttImVCo2Lo50aGfXqgsNYeugKxs47hFhpuWlSVuvQrI9dplNVxFEgKQ4o21TLiMhMmU1LTV5s2rQJPj4+qFq1KkaOHIkbN27kun9SUpKqmM58IcoiJhyY9wSQlgxU6w60fDXLwz9sPofxf5kSmsGNy+CbgXWZ0FgRmTzxy8fqYGjTsqp2/O0lxzBzW4jWYVmvi7uBGY8C8x43LTBLZKlJjXQ9/fLLL1i/fj0+/fRTbN68WbXspKXJdNrZmzRpkhoCln4pXbp0ocZMFuDSXtNCeiWqAb2mZ4x0kpkOJq86iUkrT6rbI9tUxEe9asKOI5ysjoxqe7dHDbzQpqK6LWtH/biFq0lrwi8IKF4JiL9hSmxSErSOiMyMWc5TI0XAixcvRq9evXLc5/z586hYsSLWrVuHdu3a5dhSI5d00lIjiQ3nqaEsQncAbr5AcdNJS74SHy0/gZ/v/EX+RqdAldSQdZPfiy/XncE368+o2691rIpRdybto0IUdRH4sY0psak1wFQDx0UwdS9G7/PUVKhQAd7e3jh79myuNThy8JkvRFkm1RNlm2VJaOQv8fSE5oOeNZjQUMYfW+M6VFEXIauxT9uY8/8/VEA8y5jWYZP12I7MA3bmPLiErI/FJjWXLl1SNTX+/pwjhB6gMHh6M+Dq8Sx3S0Lz7t/HMWv7BXX7495BGNK0nEZBkrka064yXu9UNSOxmcEam8JXviXQ8WPT9tp3gAvbtI6IzITZJDVxcXE4dOiQuoiQkBC1ffHiRfXYa6+9hl27duHChQuqrqZnz56oVKkSOnbsqHXoZEluXwcWDAWunQC2f50lofnf0mOYszNUtWRP7ltLFQYTZeeFNpXwcntTi42s//XbrlCtQ7I+jYebup+MacD+2VpHQ2bCbIZ079u3D23bts24PW6caYjt0KFDMX36dBw5cgRz5sxBVFSUmqDv0UcfxQcffKC6mIjyxJAGLHwaiLkMFK8MdPk8I6H5YNkJ/LbrokpoPutXG/3qB2gdLZm5Me0qISElDd9vPqdGRTk72PH3pjDJl7XbV4B/HVOCQ2SuhcIFhQtaWrl17wHbpphmDH5uA+BTTd39+epT+PZObcTkfrXwWAOOkqO8kf8+3/snGLN3XIAMjJs2uB46c9kMonyn+0JhovtycrkpoRE9pmYkNN9uOJOR0EhRMBMaut/i4Yndq2NQo9JqLqOX5h7CjnPXtQ7LOqUmAcvGAccXax0JaYhJDenfzRBg8UjTduMRQFA/tSkFnp+vOa223+pSjUXB9MCJzYe9gtCphh+S0wx4/pf9OHY5WuuwrM++WcC+GcDS0cA10/earA+TGtI/6W6StZwCGgEdPlB3Ldx/SRV4Chmi+1yrChoHSZZMJmX8amAdNKnghbikVLVO2IXrt7UOy7o0fBYo2wJIjgMWDOPEfFaKSQ3pn7svMGQxMHgeYO+I9Seu4o1FR9RDz7Usjxcf4QRq9PCkUFjWiqruXxTX45IxZOZuXI/7d/JPKoT1ofrNBFxLAJHHgZVc1dsaMakh/YqN+Hfb1g4o4oW9F27ihd8PIM1gRN96AZjQuZrqPiDKD+7ODpjzdCOU8SqCsJsJeHbOPiSm5LyUCxXAHzB9fpJOQeDAHODIAq0jokLGpIb06dYF4NtGwLKXTQWEUiscEYNnZu9FUqoBjwT64JO+QWpdH6L8VMLdCbOfagjPIg44FBaFsXMPwSBVxFQ4KrYFWr1m2l42FrjOWZ+tCZMa0p/UZGDBU0BSNBBxFLCxxZWoBAyduQcxiamoX7aYGnorKzATFYQKJdzw45AGcLSzxarjEZi08oTWIVmXNuNN9TXSQhvN1bytCf9XJ/1Z9y5w5QDg7An0m4WYFOCpWXtxNSYJlX3cMGNoA7g42mkdJelco/Je+Kx/LbX909YQ/LrTtPwGFQJJZvrNAIZvNbXckNVgUkP6cmYtsOvOAne9v0eKeymM+v0ATl2NVd0Cs1S3gKPWUZKV6FmnlFrNW7z7TzC2n+UcNoXG3Q8oVjbrjOKke0xqSD9irwKLR5i2Gw2HsUonvLX4KLaeuQ4XBzvMHNoQAcWKaB0lWZkX2lREn7qlVHG6FKmHcKh34Tu9GpjWGIi+pHUkVMCY1JA+yGofS0YC8dcB35pAh/cxbeNZzN93SU1f/+3guggK8NA6SrJCMrru4z5BqFvGE9EJKXh2zl7EJKZoHZb1MBiAzZ8CN84Afz3PFhudY1JD+iDDshs+AxQNAPrOwPITtzJmC36vRw20q+ardYRk5XPY/DCkPvw9nHHu2m28+MdB1XJDhcDW1jTM29ENCN0ObL2zXArpEpMa0o/ArsCYgzia7I9XFhxSdz3dvDyXPyCz4OPurCbnc3awxebT1zB59UmtQ7IexSsCXT43bW+aBITt0ToiKiBMasiyJcUBMVcybkbGG/DcLzLhmQGtq5TAm10CNQ2PKLOapTzwWb/aavuHzeex8mi41iFZj9oDgaD+gDENWPQMkMj1ufSISQ1ZttVvAtObqVFPMnOrJDQRMYmoWMIVUwfXhT3noiEz0712STx/Z62xVxccxtnIWK1Dsp4u6q5fAJ5lgKiLwPJXtY6ICgD/xyfLdXKFaSr0hCgY7RzVek6HL0WrmVxnDG2Ios4OWkdIlK3XO1ZF0wrFcTs5Dc//uh+xLBwuHM4equYONnaAnSOQlqp1RJTPmNSQZYqLBP5+0bTddBR+vlQaSw9dgb2tDb57vB7KebtqHSFRjqQFUVoSpXD4/LXbqsXGKCP4qOCVbgS8sAvoNc20CCbpCpMasjzyn78kNDJ826cGdpZ7IWMa+v91q45mFb21jpDonrzdnDD9ifpqKYXVx6/i560hWodkPUpUyfr/CRNK3WBSQ5Zn/2zg9CrVfBzRYSpemHccMjpWVt1+smmmGUSJzFyd0p74X/fqavvTVSexP/SW1iFZl9gI4Pf+wO4ftI6E8gmTGrIsN86ZioMBpLR5G8+uiset+BQElfLAR71rqonOiCzJE43LoFstf6QajHjxjwO4dTtZ65Csx6kVwNm1wLqJQCSH2OsBkxqyvPVcag+EsXwrjL/SEscux8DL1RHfD6mvJjgjsjSSiE/qE4Ty3q64Ep2IcfMPwcCJ+QpH/aeASu2B1ETgr+eAVCaUlo5JDVkWR1eg25eYV/lLLDoYDjtbG7UEQilPF60jI3pg7s4OmDa4HpzsbbHx1DX8sOW81iFZB2nZ7TkNcPECIo4AW+9M0EcWi0kNWYbbN0xruAA4djka76w4o7ZlBWQWBpMeVC9ZFO/2qKG2P19zCgcvsr6m0Fp/Zf4aseVz4MpBrSOih8CkhsxfWgrwW2/g156IuRqKkb/vR3KqAe2r+eL5lqZJzIj0YGDD0qq+RtaFemnuIc5fU1hq9gGq9zLNNrx4JJCapHVE9ICY1JD5kwXowg/DGH4EHyw/ibCbCQgo5oIv+teGrSzBTaSj+pqPegep7tSLN+PxztLjWodkPbpOAVxLAGlJQMxlraOhB8Skhsxb+GFgy2S1ubHC61hwOlXN6zH98frwKMIZg0l/PFwc8PXAOpB8ffHBy1h88JLWIVkH1+LAE4uAEdsAL7YAWyomNWS+pAlYmoINqbhVrgueP2RabVvm9QgK8NA6OqIC06CcF15qZ5og7n9LjiP0xm2tQ7IO/rVNgxHIYjGpIfMlRXuRx2Eo4o0hEY8h1QBVbyDzehDp3ehHKqFROS/EJaXi5XmHkJpmKpSnQiCDEnZ+B6z/QOtI6D4xqSHzFH4E2DZFbc4oOhrHohxRxqsIPu4TxAn2yCrIdAVTBtSGu5M9DlyM4jDvwhS2C1g9Adj6BRC2V+to6D4wqSEzZQSKV8JFvw746EIVtVDlN4PqcuVtsioBxYpg4p1h3l+uPa2mM6BCULYZUGug6f+hpaM4GsqCMKkh8+RfG2d6L0efSwMz5qORdXKIrE3feqXQsYavWkZBuqESU9K0Dsk6dJpkGg11/RSw5TOto6E8YlJD5uXOBHvyH/foecG4nuqCVlVK4DnOR0NWSrpbP+4dpFb1PhMZh89Xn9I6JOtQxOvfSfnUtBJHtI6I8oBJDZkPQxowp7vqx5684hhOXY1V/5FzPhqydsXdnDC5X5DanrE9BLvO39A6JOtQvSdQrYdpUj7phpKJQMmsMakh87FrOhC6DalbpmD5TtNfRZ/3r4US7k5aR0akuUcCfdWMw0Yj8PrCI4hPTtU6JOvQ5XPA2ROIDAbCdmsdDd0DkxoyD7cuABs+VJufGp7AVXhhWLNyaFPVR+vIiMzGW12roaSHs5ptePIqdkMVCndfoM+PwPObgHIttI6G7oFJDWlP/vRc9jKQmoBTzrXx0+2WqOzjhvGdA7WOjMjsVvOe1LeW2p694wJ2sxuqcFTpCPiZuv/IvDGpIe0dXQCc24A0W0eMiH4SDna2+GpgHTg72GkdGZHZaV2lBAY0KK22X190BAnJHA1VqCKOAieXax0F5YBJDWnr9g1g1Xi1+W1qb4QY/fHqo1VRoySXQSDKyVvdqsHfwxmhN+IxefVJrcOxHqE7gR9am5ZviY3QOhrKBpMa0talPTAmxSLUvhy+Te6KJhW8OHyb6B5kEspJfYIyuqEOXLyldUjWIaChqRsqKRpY+YbW0VA2mNSQtqp2xvwGczHi9nA4OTnjs34cvk2UF1JE36deKVWSNn7RESTL4mhUsOzsgR7fADZ2QPAS4NRKrSOiuzCpIU2dvhqL/21LxgljWfyvWzWU9iqidUhEFuN/XaujuKsjTl+Nw/ebz2kdjvWs5N10lGl7+atAUqzWEVEmTGpIG3tnIPXiHrwy/zCS0wxoW7UEHrtT/EhEeVPM1RHvdK+utr/dcBZnI+O0Dsk6tJkAeJYFYi5lTEVB5oFJDRW+yBPAytdhN/NRpFw5Cg8XB3zStxZX3yZ6AD1ql0SbqiXUHwcT/joCg8GodUj651gE6PalaXvPj0DEMa0jojuY1JAGc9KMAwypWGeoj5PGMni/Zw34FnXWOjIiiyR/DHzYqyaKONph74Vb+GPPRa1Dsg6V2gH1nwI6TgJKcE4tc8GkhgrXoT+AizuQCCdMTH4SnWr4qb80iejBBRQroqZCEJ+uOolrsUlah2Qdun8FNBlhKiAms8CkhgpP/E1gzdtqc0pKH8QX8ccHvWqy24koHwxtVg5BpTwQm5iKj5YHax2O9UlJABI4tF5rTGqo8Kx9B0i4iVPG0piZ1hnvdq/BxSqJ8omdrQ0+6i1/JABLDl3B9rPXtQ7JelzcBXzXFFjxmtaRWD0mNVQ4rhwCDv6qNt9MfhqtA0uiZx12OxHlp1oBnhjSpKza/t+SY0hK5RIKhcLO0bQoryz5cn6T1tFYNSY1VDj8amFLtXfwQ2pXnHasgY96B7HbiagAvNqxqmoBPX/9Nn7YfF7rcKxDqXpAo+dM2zIQIiVR64isFpMaKhQXbibguaPVMSn1cbzVtRr8PDjaiaigllB4u2s1tf3txrMIvXFb65CswyNvA26+wM1zwI6pWkdjtZjUUMGKv6nWdprw11EkpRrQvFJxDGjISfaICpKMKGxRyVstnfD+PywaLhTOHsCjH5m2t34BRHFovRaY1FDBWvkG4r9sAJsLm+HsYItJvTnJHlFBk+/Yuz1qwN7WButPRmL9iatah2QdgvoBZVsAqQnAqglaR2OVmNRQwbmwHTg6Hy6JVxFrLIKX21dBmeJc24moMFTyccMzLcqr7feXBSMxhUXDBU7+YOvymWnBy7Rk1tZogEkNFYy0FGDFq2rzj9RHkOZXJ+M/WCIqHC+2qwwfdyeE3ojHz1tZNFwofKsDI3cAg+cDDqwdLGxMaqhg7PkJiAzGTaMbvkgbgE/6BsHejr9uRIXJzcleFeanFw1fjkrQOiTr4BNoarWhQsezDOW/2AgYN5oK5j5NHYQ+zYPU/BlEpE3RcKNyXkhMMXCm4cJ2+zqw7GUgKkzrSKwGkxrKf+vehU1yHA4ZKmK7WyeM61BF64iIrLpo+L2eNWBrA6w4GoFd529oHZL1WDoa2DcTWPOW1pFYDbNJarZs2YLu3bujZMmS6ku4ZMmSLI8bjUa888478Pf3h4uLC9q3b48zZ85oFi/lIC0VMbGxMBht8E7KMLzfOwiuTlzsjUhL1fyLYnDjMmr7vX+CkWYwah2S9cxdY2MLBC8FQrZoHY1VMJuk5vbt26hduzamTZuW7eOTJ0/GN998g++//x67d++Gq6srOnbsiMREVpebE4ONHZ6+PQodkiejZPXmeCTQV+uQiAjAuA5V4e5sjxPhMZi/j90hhcKvJtDgGdP2yvHqjz6ykqSmc+fO+PDDD9G7d+//PCatNF999RXefvtt9OzZE7Vq1cIvv/yCK1eu/KdFh7S1cP8l7Au9hXCHMnine3WtwyGiO7xcHfFSu8pq+/PVpxCTmKJ1SNah7ZuAsycQeRzYP0vraHTPbJKa3ISEhCAiIkJ1OaXz8PBA48aNsXPnzhyfl5SUhJiYmCwXKiCJMUhcMhY/r9imbsqcNCU9XbSOiogyebJpOVTwdsWN28n4dsNZrcOxDkW8TN1QQgZQxN/UOiJds4ikRhIa4eubtStDbqc/lp1Jkyap5Cf9Uro0p+cvMFsmw/nQLHyR9gkCfd0wrHk5rSMiors42tvi7W6mId6ztocg5DrXhSoU9Z8CfKoDCbdMSyiQdSc1D2rChAmIjo7OuISFsR+5QFw/C8Ou79XmF6n98WHvIDhwThois9S2qg9aVSmBlDQjPl15UutwrIOdPdDpE6DR80DLV7SORtcs4szj5+enrq9ezbp+idxOfyw7Tk5OKFq0aJYL5T/DqvGwNaRgQ1od+NbrgQblvLQOiYhyIKNLZRVvGeK96ngE9l5gd0ihqNDatISCdEeRdSc15cuXV8nL+vXrM+6T+hgZBdW0aVNNY7N6Z9bB9uxaJBvt8JX9U3ijc6DWERHRPVTxdceAhqbu+A+Xn1CDMagQyfsdd03rKHTJbJKauLg4HDp0SF3Si4Nl++LFi+ovi7Fjx6rRUX///TeOHj2KJ598Us1p06tXL61Dt15pqUhd9abanJPWEY91bKtGWBCR+ZNi/iKOdjgcFoVlR8K1Dsd6RF8CfukBzOgApCZpHY3umE1Ss2/fPtStW1ddxLhx49S2TLgnXn/9dbz44ot4/vnn0bBhQ5UErVq1Cs7OXDBMMwd/hf2NU2p9p7XeT2JQI9PkXkRk/nyKOmN4q4pq+9NVJ5GUylW8C4UM7752GrgVAuz5UetodMfGaEXtjtJlJaOgpGiY9TUP7/DZi9gx+02EGn3R//m3UL8s+4qJLEl8cirafLYJkbFJeKtLNTzXqoLWIVmHg78BS0cBTh7AmAOAq7fWEenm/G02LTVkWWSa9bdXXlQLVqbUfpIJDZEFKuJoj1cfraq2p244g6j4ZK1Dsg61BwF+QUBSNLBpktbR6AqTGrp/SbGYv/cijl6OhruTPd7obPpPkYgsT9/6Aajq646YxFRM33RO63Csg60d0PFOMrNvFhDJofX5hUkN3beUhcNRfWU/VLW5iLEdqsDHnXVNRJbKztYm4w+TWTsu4EpUgtYhWYfyLYHAboAxjat45yMmNXR/LmyDw5nlqGE8i9JebniyaVmtIyKifJiQr1E5LySnGvDVutNah2M9OrwP2DoAt0K5fEI+sb+fnaOiorB48WJs3boVoaGhiI+PR4kSJdQoJVkxu1mzZvkVF5kjgwGJy8dD2mX+SGuHJ3t24szBRDog02bIHFN9p+9Qi9I+17ICKvu6ax2W/hWvCAz9GwhoCNg5aB2NLuTpjCSrYT/77LPw9/dXc8UkJCSgTp06aNeuHQICArBx40Z06NAB1atXx7x58wo+atLG0QVwvnYUMUYXHKgwXE21TkT6UL9sMTxa3RcGI/DZ6lNah2M9yjZjQlPYLTXSEjN06FDs379fJS7ZkURnyZIl+Oqrr9QaS6+++mp+xklaS0lE4up3VSvNj4YeeKkHW+WI9Ob1TlWx7sRVrAm+iv2ht1SiQ4UkNRk4MMc0MsrJTeto9J3UBAcHo3jx4rnu4+LigkGDBqnLjRs38is+MhOpO6fDOf4Krhi9YGg8EuW9XbUOiYjyWSUfd/SvXxrz9oVh8qqTmPt8E9U1RYXgzwHAuQ1A/A2gzXito9F399O9EpqH3Z/MnNGIaweWqs3vbQdjRPuaWkdERAVkbIfKcLS3xe6Qm9h29rrW4ViPek+arrd/A8RmXbyZCqhQOHONzbZt2xAZGQmDwZDlsTFjxjzIS5IZu3E7GZ2i3kDz5B1o1fMZFHVm/y+RXvl7uODxxmUwa/sFfL76FFpU8mZrTWGo3stUMHxpL7DpY6D711pHZB1JzezZszF8+HA4OjqqFpnMv+yyzaRGf75efwbRiQZc8H8UUxtyCDeR3r3QphLm7gnD4UvRWBt8FY/W8NM6JP2Tc2mHD4BZnYADvwCNRwI+gVpHZXHuezzu//73P7XIpKy/cOHCBbWadvrl/PnzBRMlaSZ85zws3H1Wbb/dtZqaqIuI9K2EuxOeal5ObU9ZexoGGRJFBa9s0zsT8hmAdRO1jsY6khqZm2bgwIGwteX8JLoXthf+q5/HavtX0aWqO5pV4qJrRNZCVvB2d7bHyYhY/HPkitbhWI/27wI2dsDpVUDIFq2jsTj3nZk888wzWLBgQcFEQ+bDaET0PxPU5i5jdYzrWk/riIioEHkUccDzLU2rdn+17gxS0rLWT1IB8a4MNHgaqNQBcPXROhqLY2M0Gu+rXTEtLQ3dunVT89IEBQXBwSFr0eiUKVNg6UuXE5B2chXs5g5AktEB39ach1f6t9M6JCIqZHFJqWg1eSNu3k7G5H618FiD0lqHZB3SUgG7BxrHo1t5PX/f97s2adIkrF69GlWrmhZAu7tQmHTAkIa45W/DA8DvNp3xVJeWWkdERBpwc7LHiNYV8PGKk5i64Qx61y3FpVEKAxOaB3bf79wXX3yBmTNnYtiwYQ/+U8msJR38Ex6xZxBtLALbli/Dy9VR65CISCNPNCmLH7ecR9jNBPx14BIGNCyjdUjWQ+ar2TIZ8K/97zw2lKv7TrmdnJzQvHnz+30aWYqURCSt+UBt/u7QDwNb19Y6IiLSUBFHaa2pqLanbjirVvKmQhK8BNj7M7DhIyA5Xuto9JnUvPTSS5g6dWrBREOau3HzOnYnlka40QtlurwMZwc7rUMiIo093rgsvN2ccOlWAhYduKR1ONaj/jDAsywQFwHsnq51NPrsftqzZw82bNiAZcuWoUaNGv8pFP7rr7/yMz4qZN/sisKcpJfRxN8Wf9Qpr3U4RGQGXBztMLJNRXywLBjfbjiLvvUC1FIKVMDsnYBH3gb+eg7Y9hVQ/ymgiJfWUZm1+/6t9PT0RJ8+fdC6dWt4e3urauTMF7JcF67fxu+7L6rtMV0bwpYT7RHRHbJ0go+7Ey5HJWD+vjCtw7EeNfsBvkFAUgyw9Quto9FfS82sWbMKJhLSVlwkwn59BSUMHVC1aiAn2iOiLKQrWlpr3vsnGNM3nVPDu9laUwhkoluZkO/3vsCen4AmIwGPAK2jMlv8jSTl2oqP0DL6b3ztOA1vdOJ6I0T0X4MalVFLKEhrjYyEokJSqR1QriWQlgRsNd+54CwmqenUqRN27dp1z/1iY2Px6aefYtq0afkRGxUS461QFAv+TW0fKDcc1fw5MSERZd9aM7yVaZbh7zadQypnGS4cMgdcu3eAJqOAtm9qHY3ldz/1798fffv2VTUz3bt3R4MGDVCyZEk4Ozvj1q1bCA4OxrZt27BixQp07doVn332WcFHTvkm/O/3UBKp2Gmoge59BmkdDhGZscGNy6iE5uLNeCw9dAV967MrpFCUbmS60MMnNbLe0xNPPKHWfJo3bx5+/PFHNVVx+izC1atXR8eOHbF3715Uq1YtLy9JZsIQeQq+IYvV9umaL6Opp4vWIRGRmc9b82zL8pi86hSmbTyLXnVLwY6DCgqfzFvjWETrKCx/7ad0ktTI+k/Fixf/z7Buc8W1n/7r8o+PodSV1diABqj3+kp4FuHswUR07zWhmn+yAdEJKfhmUF30qF1S65Csx7VTwMrXAXtnYPA8WIuYPJ6/H7hQWF7cz8/PYhIa+q+UsIMqoTEYbXC1/qtMaIgoz2tCPd3cNI/VtxvOwGB4oL+N6UHY2AEhW4HTq4CwPVpHY3Y4+smK/XXeFt+l9sAS23bo2bGD1uEQkQUZ1rwc3J3scfpqHNYEX9U6HOvhXQmoM9i0vf594ME6W3SLSY2Vik9OxefbrmFy6kDEPfqF6icnIsorDxcHPNmsrNqevuksHrCSgR5E6zcAO0fgwlbg/CatozErTGqskdGIWdtCcC02CaW9XDCQq+4S0QN4qnl5ONnb4vClaOw4d0PrcKyHZ2mgwTOm7Q0fsLUmEyY1Viju5AY03jIEjWxOYFyHKpwVlIgeiCxyKRPyCRkJRYWo5TjAoQhweT9werXW0ZiN+z6bDR06FFu2bCmYaKjgGY2IWv4uGuAEHnc/iB61S2kdERFZsOdaVYC9rY1qqTkUFqV1ONbDzQdo9Jxp++CvWkdjuUmNDKdq3749KleujI8//hiXL18umMioQEQdXYmAuCNINDqgWMfxnF+CiB5KKU8XNVeN+I6tNYWr2UtA96+B/rO1jsRyk5olS5aoRGbkyJFqIr5y5cqhc+fOWLhwIVJSUgomSsofRiNur3pPba4u0g0t69XUOiIi0oERrSuqmfxlFNTpq7Fah2M9XIsD9YcBdpxaJd0DFVOUKFEC48aNw+HDh7F7925UqlQJQ4YMUUsnvPzyyzhz5syDvCwVsBsHlqJU/EncNjrBr8t4NRs0EdHDquTjho7V/dT295vOaR2OdUpLAW6GwNo9VIVoeHg41q5dqy52dnbo0qULjh49qpZN+PLLL/MvSnp4BgOS1n6gNte69USjmlW1joiIdOSFthXV9d+Hr6hVvKkQhR8BptYH/hgAGNJgze47qZEupkWLFqFbt24oW7asWg9q7NixuHLlCubMmYN169Zh/vz5eP/99wsmYnog1w78jZKJZxFrdEGZbmylIaL8VSvAE80qFkeqwYgZW9liUKiKlQUSo4Hrp4CjC2HN7nvGNX9/fxgMBgwaNAh79uxBnTp1/rNP27Zt4enpmV8xUj6YfK40DMkjUMvXAUOrmf6iIiLKT8NbV1SjoObuvYgx7Spx6ZXC4uwBNB9jmmF48ydAzb6AnXVOqHrfLTXSrSStMtOmTcs2oRGS0ISEMFM3F2cj47DoUAQWGVqhTp9XtQ6HiHSqVWVvBPq5Iz45Db/tCtU6HOvSaDhQpDhw8zxwxHoWunzopEYKgp2dnQsmGsp/BgOmrj0OWW+uQ3Vf1C7NFjQiKhjSrS0jocTsHReQmGLd9R2FyskNaP6SaXvLZ0BaKqwRp5LVucs75+K104PRx3aLmj2YiKggda3lr+auuR6XjEUHLmkdjnVp+KypteZWCHB0PqwRkxo9Mxhgs/lTBNhcRzu/RFTzL6p1RESkcw52tnimRXm1/dOW80iTZmIqHI6uQLMxpu3QHbBGTGp07NKOuSiZfAExxiKo2us1rcMhIisxoGFptYr3hRvxWHM8QutwrK+15qlVQM9vYY2Y1Oi8lUZsK94flcoEaB0REVkJVyd7PNHEtNDlz9s4aKTQa2vKNoW1YlKjU2E7/kSplPRWmte1DoeIrMzQpuXgaGeL/aG3cODiLa3DsU5x14DLB2BNmNTokcEA282T1eY27/6oyFYaIipkPkWd0aNOSbXNyfg0ELIF+LoW8NdzVjUSikmNDoUcWP1vK01PttIQkTbSC4ZXHgtH2M14rcOxLiXrAvbOwI2zwLFFsBZManTo4+AS6J30Hv4q+QpbaYhIMzLismVlbzVP1qztF7QOx7o4uQPNRpu2t0y2mjWhmNTozLHL0VgbfBWHUBkteo/QOhwisnLPtqygruftvYjohBStw7EujZ4HnD1NrTXBS2ANmNToidGIWWv3qs0etUuiko+b1hERkZWTpROq+LrhdnIa5u65qHU41tda0+QF0/aWz1W9pd4xqdGRi3v+xkchA/G6/Vy8+EglrcMhIlJLJzzbwtRaM2fHBaSm6f/EalYaDwecigKRwcCpFdA7JjV6YTQiZeOncLZJQaC3Iyr5uGsdERGRIqOgirs64kp0IlYfv6p1ONbFxdPUDWXnaFo+QeeY1OjE5UOrUTHxOJKMDijd7Q2twyEiyuDsYIfHG5sm45u5Xf8nVrPTbDTw0mGg2YvQOyY1OnF77SfqeodHZ1SuxIUrici8PNGkLBzsbNRkfIfDorQOx7q4FAOKmuYM0jsmNTpw5ehGVIk/iGSjHfy7TNA6HCKibCfj617LdGKdxdYa7Vw5CIQfhl4xqdGB2NUfq+ud7o8iMLC61uEQEWXrqeamyfiWHQnH1ZhErcOxPnt/Bn5sA6x6E3plMUnNu+++q6roM18CAwNh7a5cPIdysQeRarSFd6fxWodDRJSjoAAPNCxXDKkGI37bFap1ONanSifA1gEI3QaE7oQeWUxSI2rUqIHw8PCMy7Zt22Dtpu6LR4ukr/BD8TdQo2YdrcMhIspTa83vuy8iMcU6Zrk1Gx4BQJ3Bpu2tX0CPLCqpsbe3h5+fX8bF29sb1iwiOhGL9l/CNRRDw+7Pax0OEdE9PVrdF6U8XXDzdjL+PnxF63CsT4uxgI0tcHatLmtrLCqpOXPmDEqWLIkKFSrg8ccfx8WLuc9OmZSUhJiYmCwXPZm3bgeS0wxoVM4Ljcp7aR0OEdE92dvZqpFQ6ZPxGY1GrUOyLl4VgJp9TdvbvoTeWExS07hxY8yePRurVq3C9OnTERISgpYtWyI2NjbH50yaNAkeHh4Zl9KlS0Mvoq6cxQtH+uEXh0l4sRUXrSQiyzGwYWk42dvi+JUYHLh4S+twrE+Ll03Xx5cA189CTywmqencuTP69++PWrVqoWPHjlixYgWioqIwf/78HJ8zYcIEREdHZ1zCwsKgFyF/fwIHmzS4O9uhRTUmNURkOYq5OqJnHdPw7tk7WDBc6HxrAFU6A+5+QJS+3n97WChPT09UqVIFZ8/mnGU6OTmpi97EXb+E6hGmFVdTmo5VI8GIiCzJk03LYf6+S1h5NByRXaupeWyoEHX/2rSEgr2+zpEW01Jzt7i4OJw7dw7+/v6wNqeXToYTUhBsVxUNWvfUOhwiovtWs5QHGpQ1De+WkVBUyNx9dZfQWFRS8+qrr2Lz5s24cOECduzYgd69e8POzg6DBg2CNUmIvoEqYaYut6j6L8LWzmI+QiKiLJ5sVk5d/7HnIpJTuXq3JgxpwNGFwO3r0AOLOSNeunRJJTBVq1bFY489huLFi2PXrl0oUaIErMnJf76AGxJw1qYsGj5qXQkdEelLpxp+8HF3wrXYJKw8Fq51ONZp4dPAomeAXdOhBxaT1MydOxdXrlxRw7QlwZHbFStWhDVJTkmD67kVajsiaAQc7C22JIqICI72thh8Z/XuX3bqq2DVYtS8M7x7z09AouVPe2IxSQ0Bfx8JR7eEd/GO3Rg06PqM1uEQET20wY3KwN7WtHp38BXLP6lanMBugHcVICka2D8blo5JjYUwGIz4YfM5JMMBJVsNg7MOR3URkfWRUU8da/ip7d92s7Wm0NnaAs3GmLZ3fQekJsGSMamxENsOHMa5yBi4O9lnNNcSEelB+gzDSw5eRmxiitbhWJ9ajwHu/kBsOHAk57nfLAGTGgtgNBhQetVT2OD4Cl4JikdRZwetQyIiyjdNKnihko8b4pPTsPjgZa3DsT72TkCTF0zb27+WrgFYKiY1FuD0jqUon3oeJWyi0LVVE63DISLKVzKB6BN3WqB/3RnK9aC0UH8Y4OQBFPECbl+DpWJSYwEMW02Ljh0s0QslfKxvskEi0r8+9QPg4mCHM5Fx2B1yU+twrI9zUeCFncAza0wT81koJjVmLuTwZlRLOowUox3Kdn1V63CIiAqEdKv3qltKbf+2iwXDmvAwvf+WjEmNmYtZ+7m63u/RHqXLV9E6HCKiAvNEE1MX1KpjEYiMTdQ6HOsVfxM4sgCWiEmNGQs/dxRBsVvVdvFHX9M6HCKiAlWjpAfqlvFU60Et2HdJ63CsU0IU8HVt4K9ngYhjsDRMaszY8fW/wdbGiAPOjVG5ZkOtwyEiKpTJ+MTcvRfV/FxUyFw8gYqPmLZ3TIWlYVJjpm7dTsaLYW3QN2kibB75n9bhEBEVim61SsLd2R5hNxOw9aw+Flm0OM1fMl0fWwhEW9YQeyY1ZkoK5RJSDEjwa4g6DVtoHQ4RUaFwcbRD33oBavsPzjCsjVL1gLItAEMqsNuyFrpkUmOGEuPj8NcOU1/m860qqDkciIisRfqs6etOROJqDAuGNdHsRdP1vtlAYjQsBZMaM3Rs2TQsSx2O11xXomstzktDRNaliq87GpQthjSDEfP3hmkdjnWq/CjgXRVIjgX2z4GlYFJjZgypqSh5YiZcbZJQp2IpONjxIyIi622tmbs3TCU3pMVCl6MBW3sg3nJqm3jGNDNH1/+GksYI3II7avcYpXU4RESa6BLkDw8XB1yOSsCW05Y7bb9FqzUAeOkI0OF9WAomNebEaESRvd+pzeBS/eHmVlTriIiINOHsYId+9U0Fw3/uuah1ONa70KWHZc0yzKTGjJzaswaVU08hyeiAKt3HaR0OEZGmBjYsra7Xn4zkDMNau3ociDgKc8ekxowkbP5KXR/06oQSfqYvMxGRtars6476dwqGF+7nDMOa2TsDmN4MWPM2zB2TGjMRFhaKarf3qG3fjq9oHQ4RkVkYcKe1Zt7eMBiNLBjWRKX2gI0tcH6T2S+dwKTGTPx8IBZtkr7EjOKvoHxgXa3DISIyC91q+cPNyR6hN+Kx8/wNrcOxTsXKAtV7mrZ3meo+zRWTGjMQHZ+C+fsuIRzFEdj5Ba3DISIyG0Uc7dGjTsmM1hrSSNPRpusj84HYCJgrJjVmYMGOYCSkpCHQzx3NKhbXOhwiIrMsGF55LAJR8clah2OdAhoAAY0AQwqw92eYKyY1GktOSkSPbb0w02EyRjd055IIRER3CSrlger+RZGcasDig5a1wKKuNB31b+FwcjzMEZMajR1ZNQM+uImadhfRoUGg1uEQEZkd+WNvYCMWDGsusBvgWQawtQOun4I5YlKjIaPBAK8jP6ntc+UHw8nJReuQiIjMUs86peBkb4uTEbE4etlyFljUFTt7YPB8YOwxoKR5DmhhUqOh4B3LUSEtBPFGJ1TrOkbrcIiIzJYsmdCppp/anr+PBcOa8akGODjDXDGp0VDa9m/U9ZESXeHp7at1OEREZu2xBqYuqKWHriAxJU3rcKybwQCE7oC5YVKjkbDTh1ArYQ8MRhuU7Piy1uEQEZm9phWKo5SnC2ITU7H6uPkOK9a91GTguybArM5A+BGYEyY1GglbO01dH3ZtgjKVa2kdDhGR2bO1tUH/BqZFLtkFpSF7R8C3hllOxsekRqPJ9l6I6I43Up6DXatXtQ6HiMhiyMrdMvPF9rM3EHbTPIcVW9Xw7mOLgNirMBdMajQwd+9FRKXY4XCJHghq3E7rcIiILEZAsSJoXtFbbXORS60n42sIpCUD+2bCXDCpKWSpqan4Zft5tf10i/KcbI+I6D6ld0FJUmMwcM4azTQZabreNwNITYI5YFJTyI6s+QVzEsfgCZdd6FHbtJ4JERHlXccafijqbI/LUQnYcY6LXGqmWg+gaCng9jVTN5QZYFJTyFwP/ohKtlfQNSAezg52WodDRGRx5P/O9EUuFx1gF5Rm7ByAhs+ats9vgjlgUlOITh/YhKopJ5BstEPlri9pHQ4RkcXqW8/UBbXyWDhiE1O0Dsd61R8GDFsB9P4B5oBJTSGK2TRVXR/2bAdvvzJah0NEZLHqlPZExRKuSEwxYOVRzlmjmSJeQLnmskAXzAGTmkJy7XIIakdvVNuebbkkAhHRw5BBFn3r/1swTGYgMQaIv6lpCExqCsm5FV/BwSYNwQ41UblOS63DISKyeH3qBsDWBthz4SZCb9zWOhzrtn82MKU6sPULTcNgUlMIEhNuo+rlv0zb9Z/TOhwiIl3w83BG80qmOWsWHbisdTjWzb0kkBwLnFlrWhdKI0xqCsE/x27gheQX8Y9de9RqN1jrcIiIdDXDsPjrAOes0VSl9sBjvwIjtsl6FpqFwaSmgBmNRszeGYqdhhq43Goy7B0ctQ6JiEhXc9a4O9nj0q0E7A7Rtp7DqtnaAtV7mNaF0jIMTX+6FdgfegvHr8TAyd4WAxqU1jocIiLdzVnTrba/2mbBMDGpKWAJi8dggv3vGFLDCcVc2UpDRJTf+tyZs2bVsXAkJKdpHQ5piElNAbp66RyaRi3HcPvlGFSziNbhEBHpUoOyxVDaywW3k9OwJphz1lgzJjUF6PzKqbC3MSDYMQgVg5poHQ4RkW7nrOldN71gmKOgrBmTmgIcxh142bTAV2K9O2tjEBFRgehdt5S63nrmGiJjErUOhzTCpKaAHF09C8UQgwh4cxg3EVEBK+/tinplPCGjuv8+fEXrcEgjTGoKgNFgQLGjM9V2SPmBHMZNRFQIet8pGGYXlPViUlMATu3bgEpp55BodEBgl9Fah0NEZBW6BfnDwc4GweExOBkRo3U4pAEmNQVg7qlUzEjtjD3Fe6FYCdP8CUREVLBk2oxHAn3U9mK21lglJjX57FpsEn4/acAHqUPg1Vfbhb2IiKxN+iioJYcuI43LJlgdJjX5bO6ei0hJM6JuGU/ULOWhdThERFalbWAJeLg44GpMEnafv6F1OFTImNTko9SUZJTePgGNbU7gySZltA6HiMjqONnboUuQX0ZrDVkXJjX56Mj6P9DLsBbTnb5Gl+rFtQ6HiMgq9axjmrNm5dEIJKZw2QRrwqQmHzkdMA3jPlWqD5ycuSwCEZEWGpXzQkkPZ8QmpWLDyUitw6FCZHFJzbRp01CuXDk4OzujcePG2LNnD8xB6In9qJF8GGlGG5TvyGHcRERasbW1QY87rTVLDrILyppYVFIzb948jBs3DhMnTsSBAwdQu3ZtdOzYEZGR2mfiEeu/VddH3JrDr0xlrcMhIrJqveqWVNebTl1DdHyK1uFQIbGopGbKlCl47rnn8NRTT6F69er4/vvvUaRIEcycaer20Ups9E3UvLZCbTs0eV7TWIiICAj0K4pAP3ckpxmw4li41uFQIbGYpCY5ORn79+9H+/btM+6ztbVVt3fu3KlpbMGrfoSrTSJCbQNQo3l3TWMhIqKsBcPsgrIeFpPUXL9+HWlpafD19c1yv9yOiIjI9jlJSUmIiYnJcslvRqMR6y4k45zBHxFVnoCNrcW8pUREutajjqkLanfITVyJStA6HCoEuj4DT5o0CR4eHhmX0qVL5/vPSEwxILV6Xwxxnorq3V/K99cnIqIHU8rTBY3Ke6ntZUe4crc1sJikxtvbG3Z2drh69WqW++W2n59poqW7TZgwAdHR0RmXsLCwfI/LxdEOE7vXwNbx7eHuymHcRETmpEdtU2vN34eZ1FgDi0lqHB0dUb9+faxfvz7jPoPBoG43bdo02+c4OTmhaNGiWS4Fxc7WpsBem4iIHkyXIH/Y29rg2OUYnLsWp3U4VMAsJqkRMpz7p59+wpw5c3DixAmMHDkSt2/fVqOhiIiI7ubl6ogWlb3V9t+H2Fqjd/awIAMGDMC1a9fwzjvvqOLgOnXqYNWqVf8pHiYiIsrcBSXz1fxz+ArGtq8MGxu2rOuVjVGG71gJGf0kBcNSX1OQXVFERGQ+4pJSUf+DtUhKNWDZiy1Qs5SH1iFRAZ2/Lar7iYiI6H65OdmjXTUftc2CYX1jUkNERFYzCkq6oAwGq+mgsDpMaoiISPfaVPWBu5M9wqMTsS/0ltbhUAFhUkNERLrn7GCHjjX9MlprSJ+Y1BARkVXoVstfXa88Fo7UNIPW4VABYFJDRERWoXklbxQr4oDrcclqPSjSHyY1RERkFRzsbNHpThcU14LSJyY1RERkNbrVMo2CWnksAinsgtIdJjVERGQ1Gpf3grebI6LiU7D97HWtw6F8xqSGiIishr2dLTrXNBUMLzsSrnU4lM+Y1BARkVWOglp9PAJJqWlah0P5iEkNERFZlYblvOBb1AmxianYeppdUHrCpIaIiKyKra0NugSld0FxFJSeMKkhIiKrHQW17kQkElPYBaUXTGqIiMjq1C3tCX8PZ8QlpWLL6Wtah0P5hEkNERFZHemCSh8FteIoR0Hphb3WAZgbg8GA5ORkrcMgui8ODg6ws7PTOgwii9K1lh9mbg/J6IKSRS/JsjGpyUSSmZCQEJXYEFkaT09P+Pn5wcbGRutQiCxC3dLF4FfUGRExidh65jo6VPfVOiR6SExq7jAajQgPD1d/7ZYuXRq2tuyZI8v53Y2Pj0dkZKS67e9valInojx0QQX5Ydb2C6oLikmN5WNSc0dqaqo6MZQsWRJFihTROhyi++Li4qKuJbHx8fFhVxRRHnUN8ldJzbrgq2oiPid7fncsGZsj7khLMw3pc3R01DoUogeSnoynpKRoHQqRxahXxtQFFZvEifj0gEnNXViPQJaKv7tED9YF1ammn9rmKCjLx6TGisyePVsVk+bl5LhkyZJ8/dmbNm1SrxsVFQWtXbhwQcVy6NAhrUMhIjPQ9c5aUGvvdEGR5WJSY0UGDBiA06dPZ9x+9913UadOnUL52c2aNVOF2B4eHtCaFIJLLDVr1rznvkyAiPSvfpliprWgklKx4+wNrcOhh8CkxsqKSaWIVAtSq5TbcGOpacrLUHoZ6SNF3Q9DimglFnt71skTkakLqmMNdkHpAZMaC7Zs2TLVnZRe5CytCZI0jB8/PmOfZ599Fk888cR/up9k+7333sPhw4fVc+Qi96W7fv06evfurYpPK1eujL///jvXWH799Vc0aNAA7u7uKmEYPHhwxhDj7Lqf0mOR161evTqcnJxw8eLF/7xu+vNWrlyJ+vXrq/22bdumEqBJkyahfPnyKlmrXbs2Fi5cmPG8W7du4fHHH0eJEiXU43IMs2bNyrb1Jbd95fVF3bp11XPatGmjbu/duxcdOnSAt7e3an1q3bo1Dhw4kCV22f/nn3/O9X08fvw4unXrhqJFi6r3rmXLljh37lzG4/L8atWqwdnZGYGBgfjuu+9y/RyI6MGk19WsPXEVKWmcq8xS8U/VXFoEEjRa5MzFwS5PRZ9yAoyNjcXBgwdVQrF582Z1kpVEIJ3c98Ybb2TbFXXs2DGsWrUK69atU/dl7hqShGfy5Mn47LPPMHXqVHXSDw0NhZeXV7axyIibDz74AFWrVlXJzLhx4zBs2DCsWLEix/hlCP2nn36qTtzFixfPtRVJErXPP/8cFSpUQLFixVRC89tvv+H7779XycKWLVtU8iaJiSQY//vf/xAcHKySIXlPzp49i4SEhGxfO7d99+zZg0aNGqn3qEaNGhmj4+R9Hzp0qHpv5Hfliy++QJcuXXDmzBmVnOTlfbx8+TJatWqlEqUNGzaoxGb79u0ZLVG///473nnnHXz77bcqqZLP+bnnnoOrq6v62USUfxqV84KXqyNu3k7G7vM30aKyt9Yh0QNgUpMDSWiqv7Nak58d/H5HFHG890cjSYjUxEgSI0mNXL/88svqRBoXF4fo6Gh1gpaT/N2kRcLNzU11wUjLyt0kIRk0aJDa/vjjj/HNN9+oE3ynTp2yjeXpp5/O2JbEQ/Zv2LChikN+Tk6JkLQ8SCvLvbz//vuqZUQkJSWpmCTRaNq0acbPlBacH374QR2vtPpIIiDviyhXrlyOr53bvpIkCUm6Mr9PjzzySJbX+PHHH1XLkySR0vKSl/dx2rRp6jOcO3euWuZAVKlSJeO5EydOVMlSnz59MlqNJPmSY2RSQ5S/7O1s8Wh1X8zdG4aVx8KZ1Fgodj9ZODmBSzIjrQVbt25VJ0DprpATvJxgZTJBacm4X7Vq1crYlpYBaUXI3J10t/3796N79+4oU6aMaqlIT6Sy61JKJ60emX9ObtITDiGJmrTySJIjCVP65Zdffsnouhk5cqRKFiTpe/3117Fjx44cX/t+9k139epV1Woi760kJvL+SAJ39/Hm9j5K95e0tqUnNJndvn1bHcszzzyT5Rg//PDDLN1TRJR/OgeZRkGtPn4VaQaj1uHQA2BLTS5dQNJiotXPzivpupg5c6aqjZGTo9RdyH2S6EitSHatNHlx94lWusNyKuSVE3DHjh3VRbpMpHVDTu5yO7fFQaW1KK9zq0hCkE6SB7F8+XKUKlUqy35ScyM6d+6sunmk+2vt2rVo164dRo0apbqw7nY/+6aTlpIbN27g66+/RtmyZdXPlVaju483t/cxfRbg7KQf408//YTGjRtneYyzBRMVjKYViqOosz2uxyVh34WbaFyhuNYh0X1iUpMDOfnkpQtIa+l1NV9++WVGAiNJzSeffKKSmldeeSXXlpL0IuOHcfLkSXWCl58pw6XFvn37UFAyFxbnlrRJciXJh1zkfXrttddyTFRy2je9hubu90lqX6TrTOpoRFhYmCquvh/SijNnzhzVDXd38uPr66ta2c6fP6/qcIio4Dna26J9dV/8deAyVh6LYFJjgdj9ZOGkaFZOjtJCkj4yR4pPZSSOzEmT20lfakdkVXLpBpETstSqPAjpcpKTvxTCyklYRvhI0XBBke6tV199VdUPSVIg3TFyvPLz5baQAtulS5eqrioZYSQjxaRbLju57SvFy9KiIgXV0uUkdUpCup1kxNeJEyewe/dulXjk1vKSndGjRyMmJgYDBw5USaAUGctrnjp1Sj0utVFSEC11OPJZHj16VI3KmjJlykO+g0SUk84107ugImBgF5TFYVKjA5K4SEtCelIjI2ukNUMKW2U0Uk769u2rClbbtm2rWir+/PPPB/r58lwZor1gwQL1c6XFJreum/wgSZOMWpKTviQgchzSHZU+BFuSrAkTJqiET5I86bKRupns5LavFFJLUiHFudJy0rNnT3X/jBkzVEtYvXr1MGTIEIwZM+a+5wCS4mMZ9SRdTfIZypB16W5Kb7WR4fgyMkwSmaCgILWPvM/px0hE+a9lZW+4OtohPDoRhy9pPwM63R8bo1SYWgn5q1iKOuWvbSnYzCwxMVG1WsgJQ+YEIbI0/B0myh+j/ziAZUfCMaJ1RYzvHKh1OITcz9+ZsaWGiIgok/TZhaULyor+7tcFJjVERESZtA30gaOdLUKu38aZSNNIRLIMTGqIiIgycXOyz5h8b/WxCK3DofvApIaIiOgune50Qa06zqTGkjCpISIiuku7aj6wtQGOX4lB2M14rcOhPGJSQ0REdJfibk5oVN4ro2CYLAOTGiIiolxGQa05flXrUCiPmNQQERHlktTsDb2Ja7EPNuM6FS4mNURERNko6emCWgEekKlq1p1ga40lYFJDD+Tdd99FnTp1YA5k6QBPT0+twyAiHbfWrA1mUmMJmNTQA5EFJdevXw9zMGDAALXgY14wASKi+9Ghuq+63nb2OuKSUrUOh+6BSQ09EDc3N7UgY06Sk5Pz9Dp53S83sjr2/S4mSUSUF5V93FCueBEkpxqw5fQ1rcOhe2BSY+FkZW5ZIfr1119Xq3PLytzSNZTZxYsX1erSkojIQmCPPfYYrl7NvSn1jTfeQJUqVVCkSBFUqFBBrYidkpKSY/fTsGHD0KtXL3z00UdqNeucVgdPf56sPp154cWoqCi1KrWs+C0xPvLIIzh8+HDG82RbVhN3d3dXj8uK1vv27cu29SWnfTdt2oSnnnpKLYhmY2OjLunv1a+//ooGDRqo58h7OHjwYERGRma8pjxX9pfWKdlP3pdmzZrh1KlTWY7vn3/+QcOGDdVxeXt7o3fv3hmPJSUlqRauUqVKwdXVFY0bN1avS0TmS773j7ILymIwqbmX5Ns5X1IS72PfhLzt+wDmzJmjTpK7d+/G5MmT8f7772Pt2rXqMYPBoBKamzdvYvPmzer+8+fPqy6b3MjJXZKF4OBgfP311/jpp5/w5Zdf5vocOeHLSV5+xrJly3Lc7+zZs1i0aBH++usvHDp0SN3Xv39/lUSsXLkS+/fvR7169dCuXTsVt3j88ccREBCAvXv3qsfHjx8PBweHbF8/p30lCfnqq69UohMeHq4ukmQISdg++OADlRAtWbIEFy5cUIna3d566y188cUXKkmyt7fH008/nfHY8uXLVRLTpUsXHDx4UL0fjRo1ynh89OjR2LlzJ+bOnYsjR46oY+7UqRPOnDmT6/tKRObRBbX+xFWkpBm0DodyY7Qi0dHRstyqur5bQkKCMTg4WF1nMbFozpff+mXd90O/nPed2SXrvp+Wz36/+9S6dWtjixYtstzXsGFD4xtvvKG216xZY7SzszNevHgx4/Hjx4+r92HPnj15/jmfffaZsX79+hm3J06caKxdu3bG7aFDhxp9fX2NSUlJub6OPM/BwcEYGRmZcd/WrVuNRYsWNSYmJmbZt2LFisYffvhBbbu7uxtnz56d7WvOmjXL6OHhkXH7fvbNyd69e9V7FBsbq25v3LhR3V63bl3GPsuXL1f3pf/ONG3a1Pj4449n+3qhoaHqc7h8+XKW+9u1a2ecMGGCMT/k+DtMRA8lNc1grPf+GmPZN5YZt525pnU4Vik6l/N3Zmyp0YFatWplue3v75/RdXLixAmULl1aXdJVr15dddfIYzmZN28emjdvrrpipNvq7bffVt1YuQkKCoKjo+M94y1btqzqZkonrSNxcXGqRkd+VvolJCQE586dU/uMGzdOdU+1b98en3zyScb92bmffdNJi0737t1RpkwZ1UrVunVrdf/dx5z5vZb3WaS/19LqJK1L2Tl69CjS0tJUl17mY5TWs7zER0TasbO1QftqptYadkGZN3utAzB7b17J+TEbu6y3Xzuby7535Y9jjyK/3N0NI33A0u30oKSLRLpw3nvvPXTs2BEeHh6qy0S6XXIjXWB5cfd+ktBIgpBdfUl6rYzUvkidi3TxSBfVxIkTVUyZa1bS3c++4vbt2+o45fL777+rhEuSGbl9dyFz5vda3meR/l5LwXJO5Bjt7OxU8iTXmUlyQ0Tm3wU1b18Y1hyPwMTu1TO+/2RemNTci6Or9vs+hGrVqiEsLExd0ltrpE5GCnOlxSY7O3bsUK0pUj+SLjQ0tMBilPqZiIgIVaNSrly5HPeTVg65vPzyyxg0aBBmzZqVY6KS077SkiQtJpmdPHkSN27cUK066e9RehHy/ZBWHKmjkWLku9WtW1f9XGnVadmy5X2/NhFpq0Vlb7g42OFKdKJa5LJmKQ+tQ6JssPtJ56QLRrqFpOXlwIED2LNnD5588knVvSKjeLJTuXJl1VIhrRvSNfLNN99g8eLFBRpj06ZN1eipNWvWqCJdSawkqZLkIiEhQRXZSkuOJFfbt29XRcCSsN3tXvtK0iStJpJ8XL9+HfHx8arLSZKdqVOnqiLqv//+WxUN3y9pEfrzzz/VtXTtSZfTp59+qh6TBEs+A3nvpUBautbks5g0aZJqUSIi8+bsYIfWVUzd5tJaQ+aJSY3OSRPp0qVLUaxYMbRq1UolEDJEW2pmctKjRw/VwiHJgQy/lgRDhnQXZIwrVqxQ8UkrhyQAAwcOVEmJr6+v6q6RlhRJCOQxGZLeuXNn1T12t3vtKyOgRowYoUZ/STeTjBaTaxnptWDBAtV6JS02n3/++QMNr5fXkKRI3jcZli6JSzppLZK4XnnlFTXkXZI4SbgkqSIiy+iCcrSzRUwiJ+EzVzZSLQwrERMTo+pDZJ4SGdabWWJiovrrOfPcKUSWhL/DRAUrITkNqQYD3J2zn06CtDl/Z8aaGiIiojxwcZQi/7sGiJBZYfcTERER6YLFJDVS4Jk+tX36RWofiIiIiCyu+0mm/3/uuecybsskaUREREQWl9SkLzZIREREZLHdT0K6m2QqfZnI7LPPPkNqav4Pq7OiwWCkM/zdJSJrZzEtNWPGjFEzz3p5eal5UyZMmKBWWZ4yZUqOz0lKSlKXzEPCcpI+db1Mi5/bdPdE5komEhQ5rV5ORKR3ms5TM378+IwZV3MiM7MGBgb+5/6ZM2di+PDhanZYJyenbJ8rawBlN0FbduPc5W2QWXRTUlJQsmRJ2NpaVCMWWTH53ZWERpZgkLWy0hfaJCKytnlqNE1qrl27pmZ/zY3Mfpvdys/Hjx9HzZo11bo9MjtrXltqZG2fnN4UaaWRycseZjFIIq1IQiM1Z1xoj4j0xiIm35Pp6eXyIA4dOqRaU3x8fHLcR1pwcmrFyY4kT7Lu0d0rMxOZO+lyunv1byIia2MRNTU7d+7E7t270bZtWzUCSm7L2kRPPPGEWtMoP0mixCnmiYiILI9FJDXS2iIrRkuNjHQnydo2ktSMGzdO69CIiIjITFhEUiOjnnbt2qV1GERERGTGOMSHiIiIdMEiWmryS/pAr9zmqyEiIiLzkn7evteAbatKamJjY9W1DOsmIiIiyzuPy9Bus5ynprDJ/DNXrlxRI6jycy6P9PlvwsLCch0/b8n0fow8Psun92Pk8Vk+vR9jTAEen6QqktDca3Jcq2qpkTciICCgwF5fPkQ9/qJa0zHy+Cyf3o+Rx2f59H6MRQvo+HJroUnHQmEiIiLSBSY1REREpAtMavJpcsCJEyfe15IMlkbvx8jjs3x6P0Yen+XT+zE6mcHxWVWhMBEREekXW2qIiIhIF5jUEBERkS4wqSEiIiJdYFJDREREusCkJo+mTZuGcuXKwdnZGY0bN8aePXty3X/BggUIDAxU+wcFBWHFihXQy/HNnj1bzcic+SLPM1dbtmxB9+7d1UyUEuuSJUvu+ZxNmzap1eGlir9SpUrqmM3Z/R6jHN/dn6FcIiIiYI4mTZqEhg0bqtnAfXx80KtXL5w6deqez7OU7+GDHJ+lfQ+nT5+OWrVqZUzM1rRpU6xcuVIXn9+DHJ+lfX53++STT1TMY8eOhTl9hkxq8mDevHkYN26cGqp24MAB1K5dGx07dkRkZGS2++/YsQODBg3CM888g4MHD6r/oORy7Ngx6OH4hHxpw8PDMy6hoaEwV7dv31bHJIlbXoSEhKBr165o27YtDh06pL60zz77LFavXg29HGM6OXFm/hzlhGqONm/ejFGjRmHXrl1Yu3YtUlJS8Oijj6rjzoklfQ8f5Pgs7Xsos7nLiXD//v3Yt28fHnnkEfTs2RPHjx+3+M/vQY7P0j6/zPbu3YsffvhBJXG50eQzlCHdlLtGjRoZR40alXE7LS3NWLJkSeOkSZOy3f+xxx4zdu3aNct9jRs3Ng4fPtyoh+ObNWuW0cPDw2iJ5Fd+8eLFue7z+uuvG2vUqJHlvgEDBhg7duxo1Msxbty4Ue1369YtoyWKjIxU8W/evDnHfSzte3i/x2fJ38N0xYoVM/7888+6+/zycnyW+vnFxsYaK1eubFy7dq2xdevWxpdeeinHfbX4DNlScw/Jyckq827fvn2WNaTk9s6dO7N9jtyfeX8hLR857W9pxyfi4uJQtmxZtXjZvf4asTSW9Pk9rDp16sDf3x8dOnTA9u3bYSmio6PVtZeXly4/x7wcnyV/D9PS0jB37lzVEiXdNHr7/PJyfJb6+Y0aNUq1ZN/92ZjLZ8ik5h6uX7+ufkF9fX2z3C+3c6o/kPvvZ39LO76qVati5syZWLp0KX777Te1+nmzZs1w6dIl6EFOn5+sQJuQkAA9kETm+++/x6JFi9RF/lNt06aN6n40d/L7Jl2CzZs3R82aNXPcz5K+hw9yfJb4PTx69Cjc3NxUrdqIESOwePFiVK9eXTef3/0cnyV+fnPnzlX/R0gNWF5o8Rla1SrdlD/kL4/Mf33IF7FatWqqj/WDDz7QNDbKG/kPVS6ZP8Nz587hyy+/xK+//gpz/0tR+uS3bdsGPcrr8Vni91B+56ROTVqiFi5ciKFDh6p6opxO/Jbmfo7P0j6/sLAwvPTSS6rmy5wLmpnU3IO3tzfs7Oxw9erVLPfLbT8/v2yfI/ffz/6Wdnx3c3BwQN26dXH27FnoQU6fnxT1ubi4QK8aNWpk9onC6NGjsWzZMjXaSwozc2NJ38MHOT5L/B46Ojqq0YSifv36quD066+/VidyPXx+93N8lvb57d+/Xw0ekVGh6aSVX35Xv/32WyQlJalzidafIbuf8vBLKr+c69evz7hPmgnldk59pXJ/5v2FZLe59a1a0vHdTX6xpdlVujT0wJI+v/wkf2Ga62co9c9ywpfm/A0bNqB8+fK6+hwf5Pj08D2U/2vkZGjpn9+DHJ+lfX7t2rVT8cn/E+mXBg0a4PHHH1fbdyc0mn2GBVaCrCNz5841Ojk5GWfPnm0MDg42Pv/880ZPT09jRESEenzIkCHG8ePHZ+y/fft2o729vfHzzz83njhxwjhx4kSjg4OD8ejRo0Y9HN97771nXL16tfHcuXPG/fv3GwcOHGh0dnY2Hj9+3Giu1foHDx5UF/mVnzJlitoODQ1Vj8uxyTGmO3/+vLFIkSLG1157TX1+06ZNM9rZ2RlXrVplNFf3e4xffvmlccmSJcYzZ86o30sZwWBra2tct26d0RyNHDlSjRTZtGmTMTw8POMSHx+fsY8lfw8f5Pgs7XsosctorpCQEOORI0fUbRsbG+OaNWss/vN7kOOztM8vO3ePfjKHz5BJTR5NnTrVWKZMGaOjo6MaAr1r164sH+zQoUOz7D9//nxjlSpV1P4yPHj58uVGvRzf2LFjM/b19fU1dunSxXjgwAGjuUofvnz3Jf2Y5FqO8e7n1KlTRx1jhQoV1PBLc3a/x/jpp58aK1asqP4T9fLyMrZp08a4YcMGo7nK7tjkkvlzseTv4YMcn6V9D59++mlj2bJlVbwlSpQwtmvXLuOEb+mf34Mcn6V9fnlJaszhM7SRfwquHYiIiIiocLCmhoiIiHSBSQ0RERHpApMaIiIi0gUmNURERKQLTGqIiIhIF5jUEBERkS4wqSEiIiJdYFJDREREusCkhogs1owZM/Doo48+1Gtcv34dPj4+uHTpUr7FRUTa4IzCRGSREhMTUaFCBSxYsADNmzd/qNd69dVXcevWLZUkEZHlYksNEVmkhQsXomjRog+d0IinnnoKv//+O27evJkvsRGRNpjUEJGmrl27Bj8/P3z88ccZ9+3YsQOOjo5Yv359js+bO3cuunfvnuW+YcOGoVevXuq1fH194enpiffffx+pqal47bXX4OXlhYCAAMyaNSvL82rUqIGSJUti8eLFBXCERFRYmNQQkaZKlCiBmTNn4t1338W+ffsQGxuLIUOGYPTo0WjXrl2Oz9u2bRsaNGjwn/s3bNiAK1euYMuWLZgyZQomTpyIbt26oVixYti9ezdGjBiB4cOH/6eGplGjRti6dWuBHCMRFQ7W1BCRWRg1ahTWrVunEpWjR49i7969cHJyynbfqKgolaRI4tKyZcssLTWbNm3C+fPnYWtr+pstMDBQFQLLviItLQ0eHh74+eefMXDgwIznjhs3DgcPHsTGjRsL/FiJqGDYF9DrEhHdl88//xw1a9ZUhb/79+/PMaERCQkJ6trZ2fk/j0lXUnpCI6QbSl43nZ2dHYoXL47IyMgsz3NxcUF8fHw+HQ0RaYHdT0RkFs6dO6e6jQwGAy5cuJDrvpKU2NjYqBFLd3NwcMhyW/bL7j75OZlJkbB0hRGR5WJSQ0SaS05OxhNPPIEBAwbggw8+wLPPPvuflpTMpIi4evXqCA4OzrcYjh07hrp16+bb6xFR4WNSQ0Sae+uttxAdHY1vvvkGb7zxBqpUqYKnn3461+d07NhRFQvnB+l2ki6vh53Ij4i0xaSGiDQlhb1fffUVfv31VzXvjNTDyLaMRJo+fXqOz3vmmWewYsUKlQw9rKVLl6JMmTJZio6JyPJw9BMRWaz+/fujXr16mDBhwkO9TpMmTTBmzBgMHjw432IjosLHlhoislifffYZ3NzcHnrtpz59+mDQoEH5FhcRaYMtNURERKQLbKkhIiIiXWBSQ0RERLrApIaIiIh0gUkNERER6QKTGiIiItIFJjVERESkC0xqiIiISBeY1BAREZEuMKkhIiIi6MH/AQF0t3JGZc+XAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "fig, ax = plt.subplots()\n", + "ax.plot(position.x / u.m, position.y / u.m, label=\"with air resistance\")\n", + "\n", + "position_no_drag = (\n", + " initial_momentum / object_mass * time + 0.5 * gravitational_acceleration * time**2\n", + ")\n", + "\n", + "ax.plot(\n", + " position_no_drag.x / u.m,\n", + " position_no_drag.y / u.m,\n", + " ls=\"--\",\n", + " label=\"no air resistance\",\n", + ")\n", + "\n", + "ax.set_xlabel(\"x (m)\")\n", + "ax.set_ylabel(\"y (m)\")\n", + "ax.set_title(\"Projectile motion\")\n", + "ax.legend()" + ] + }, + { + "cell_type": "markdown", + "id": "3b8d2d92", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".env", + "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.13.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/src/pytree_api.md b/docs/src/pytree_api.md new file mode 100644 index 00000000..b142c065 --- /dev/null +++ b/docs/src/pytree_api.md @@ -0,0 +1,5 @@ +# PyTree API + +```{eval-rst} +.. autofunction:: vector.register_pytree +``` diff --git a/pyproject.toml b/pyproject.toml index 81852680..f246b8bf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -78,6 +78,7 @@ test-extras = [ "jax", "dask_awkward", "spark-parser", + "optree>=0.16", ] docs = [ "awkward>=2", @@ -89,6 +90,8 @@ docs = [ "sphinx-book-theme>=0.0.42", "sphinx-copybutton", "sphinx-math-dollar", + "hepunits", + "matplotlib", ] test-all = [ { include-group = "test"}, @@ -148,6 +151,9 @@ isort.required-imports = [ "src/vector/_methods.py" = [ "PLC0415", ] +"src/vector/_pytree.py" = [ + "PLC0415", +] "src/vector/backends/_numba_object.py" = [ "PGH003", ] @@ -275,3 +281,12 @@ module = [ ignore_missing_imports = true disallow_untyped_defs = false disallow_untyped_calls = false + +[[tool.mypy.overrides]] +module = "vector._pytree" +# optree register_node typing requires vectors to be Collection, unnecessarily +disable_error_code = ["call-arg", "arg-type"] + +[tool.codespell] +ignore-words-list = "HEP" +ignore-regex = "[A-Za-z0-9+/]{100,}" diff --git a/src/vector/__init__.py b/src/vector/__init__.py index 2e55438f..66d09e02 100644 --- a/src/vector/__init__.py +++ b/src/vector/__init__.py @@ -32,6 +32,7 @@ Vector4D, dim, ) +from vector._pytree import register_pytree from vector._version import version as __version__ from vector.backends.awkward_constructors import Array, zip from vector.backends.awkward_constructors import Array as awk @@ -95,7 +96,6 @@ def _import_awkward() -> None: VectorSympy4D, ) - __all__: tuple[str, ...] = ( "Array", "Azimuthal", @@ -148,6 +148,7 @@ def _import_awkward() -> None: "obj", "register_awkward", "register_numba", + "register_pytree", "zip", ) diff --git a/src/vector/_pytree.py b/src/vector/_pytree.py new file mode 100644 index 00000000..0a7320d7 --- /dev/null +++ b/src/vector/_pytree.py @@ -0,0 +1,222 @@ +from __future__ import annotations + +from functools import partial +from typing import TYPE_CHECKING, Any + +import numpy + +if TYPE_CHECKING: + from optree.pytree import ReexportedPyTreeModule + +from vector._methods import ( + Vector2D, + Vector3D, + Vector4D, +) +from vector.backends.numpy import ( + MomentumNumpy2D, + MomentumNumpy3D, + MomentumNumpy4D, + VectorNumpy, + VectorNumpy2D, + VectorNumpy3D, + VectorNumpy4D, +) +from vector.backends.object import ( + MomentumObject2D, + MomentumObject3D, + MomentumObject4D, + VectorObject2D, + VectorObject3D, + VectorObject4D, +) + +Children = tuple[Any, ...] +MetaData = tuple[type, ...] + + +def _flatten2D(v: Vector2D) -> tuple[Children, MetaData]: + children = v.azimuthal.elements + metadata = type(v), type(v.azimuthal) + return children, metadata + + +def _unflatten2D(metadata: MetaData, children: Children) -> Vector2D: + backend, azimuthal = metadata + return backend(azimuthal=azimuthal(*children)) + + +def _flatten3D(v: Vector3D) -> tuple[Children, MetaData]: + children = v.azimuthal.elements, v.longitudinal.elements + metadata = type(v), type(v.azimuthal), type(v.longitudinal) + return children, metadata + + +def _unflatten3D(metadata: MetaData, children: Children) -> Vector3D: + coords_azimuthal, coords_longitudinal = children + backend, azimuthal, longitudinal = metadata + return backend( + azimuthal=azimuthal(*coords_azimuthal), + longitudinal=longitudinal(*coords_longitudinal), + ) + + +def _flatten4D(v: Vector4D) -> tuple[Children, MetaData]: + children = ( + v.azimuthal.elements, + v.longitudinal.elements, + v.temporal.elements, + ) + metadata = type(v), type(v.azimuthal), type(v.longitudinal), type(v.temporal) + return children, metadata + + +def _unflatten4D(metadata: MetaData, children: Children) -> Vector4D: + coords_azimuthal, coords_longitudinal, coords_temporal = children + backend, azimuthal, longitudinal, temporal = metadata + return backend( + azimuthal=azimuthal(*coords_azimuthal), + longitudinal=longitudinal(*coords_longitudinal), + temporal=temporal(*coords_temporal), + ) + + +def _flattenAoSdata(v: VectorNumpy) -> tuple[Children, tuple[type, numpy.dtype]]: + assert v.dtype.fields is not None + field_dtypes = [dt for dt, *_ in v.dtype.fields.values()] + target_dtype = field_dtypes[0] + if not all(fd == target_dtype for fd in field_dtypes): + raise ValueError("All fields must have the same dtype to flatten") + array = numpy.array(v).view(target_dtype) + children = (array,) + metadata = (type(v), v.dtype) + return children, metadata + + +def _unflattenAoSdata( + metadata: tuple[type, numpy.dtype], children: Children +) -> VectorNumpy: + (array,) = children + (vtype, dtype) = metadata + return array.view(dtype).view(vtype) + + +_MODULE: list[ReexportedPyTreeModule] = [] + + +def register_pytree() -> ReexportedPyTreeModule: + """Register Optree PyTree operations for vector objects. + + This module defines how vector objects are handled with the optree package. + See https://blog.scientific-python.org/pytrees/ for the rationale for these functions. + + After calling this function, + + >>> import vector + >>> vector.register_pytree() + + + the following classes can be flattened and unflattened with the `optree` package: + + - VectorObject*D + - MomentumObject*D + - VectorNumpy*D + - MomentumNumpy*D + + For example: + + >>> import optree + >>> vec = vector.obj(x=1, y=2) + >>> leaves, treedef = optree.tree_flatten(vec, namespace="vector") + >>> vec2 = optree.tree_unflatten(treedef, leaves) + >>> assert vec == vec2 + + As a convenience, we return a re-exported module that can be used without the ``namespace`` + argument. For example: + + >>> pytree = vector.register_pytree() + >>> vec = vector.obj(x=1, y=2) + >>> leaves, treedef = pytree.flatten(vec) + >>> vec2 = pytree.unflatten(treedef, leaves) + >>> assert vec == vec2 + + A ravel function is also added to the returned PyTree module, + which can be used to flatten VectorNumpy arrays into a 1D array and reconstruct them. + + >>> import numpy as np + >>> vec = vector.array({"x": np.ones(10), "y": np.ones(10)}) + >>> flat, unravel = pytree.ravel(vec) + >>> assert flat.shape == (20,) + >>> vec2 = unravel(flat) + >>> assert (vec == vec2).all() + + Note that this function requires the `optree` package to be installed. + """ + if _MODULE: + return _MODULE[0] + try: + import optree.pytree + from optree import GetAttrEntry + from optree.integrations.numpy import tree_ravel + except ImportError as e: + raise ImportError("Please install optree to use vector.pytree") from e + + pytree = optree.pytree.reexport(namespace="vector", module="vector.pytree") + + pytree.register_node( + VectorObject2D, + flatten_func=_flatten2D, + unflatten_func=_unflatten2D, + path_entry_type=GetAttrEntry, + ) + pytree.register_node( + MomentumObject2D, + flatten_func=_flatten2D, + unflatten_func=_unflatten2D, + path_entry_type=GetAttrEntry, + ) + pytree.register_node( + VectorObject3D, + flatten_func=_flatten3D, + unflatten_func=_unflatten3D, + path_entry_type=GetAttrEntry, + ) + pytree.register_node( + MomentumObject3D, + flatten_func=_flatten3D, + unflatten_func=_unflatten3D, + path_entry_type=GetAttrEntry, + ) + pytree.register_node( + VectorObject4D, + flatten_func=_flatten4D, + unflatten_func=_unflatten4D, + path_entry_type=GetAttrEntry, + ) + pytree.register_node( + MomentumObject4D, + flatten_func=_flatten4D, + unflatten_func=_unflatten4D, + path_entry_type=GetAttrEntry, + ) + + for cls in ( + VectorNumpy2D, + MomentumNumpy2D, + VectorNumpy3D, + MomentumNumpy3D, + VectorNumpy4D, + MomentumNumpy4D, + ): + pytree.register_node( + cls, + flatten_func=_flattenAoSdata, + unflatten_func=_unflattenAoSdata, + path_entry_type=GetAttrEntry, + ) + + # A convenience function + pytree.ravel = partial(tree_ravel, namespace="vector") # type: ignore[attr-defined] + + _MODULE.append(pytree) + return pytree diff --git a/tests/test_pytree.py b/tests/test_pytree.py new file mode 100644 index 00000000..cbe6edcb --- /dev/null +++ b/tests/test_pytree.py @@ -0,0 +1,173 @@ +# Copyright (c) 2025, Nick Smith +# +# Distributed under the 3-clause BSD license, see accompanying file LICENSE +# or https://github.com/scikit-hep/vector for details. + +from __future__ import annotations + +import numpy as np +import pytest + +import vector + +try: + pytree = vector.register_pytree() +except ImportError: + pytest.skip("optree is not installed", allow_module_level=True) + + +def test_pytree_roundtrip_VectorObject2D(): + vec = vector.obj(x=1, y=2) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + + vec = vector.obj(rho=1, phi=2) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + + +def test_pytree_roundtrip_MomentumObject2D(): + vec = vector.obj(px=1, py=2) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + + vec = vector.obj(pt=1, phi=2) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + + +def test_pytree_roundtrip_VectorObject3D(): + vec = vector.obj(x=1, y=2, z=3) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + + vec = vector.obj(rho=1, phi=2, z=3) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + + vec = vector.obj(x=1, y=2, theta=3) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + + +def test_pytree_roundtrip_MomentumObject3D(): + vec = vector.obj(px=1, py=2, pz=3) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + + vec = vector.obj(pt=1, phi=2, pz=3) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + + vec = vector.obj(px=1, py=2, theta=3) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + + +def test_pytree_roundtrip_VectorObject4D(): + vec = vector.obj(x=1, y=2, z=3, t=4) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3, 4] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + + vec = vector.obj(rho=1, phi=2, z=3, t=4) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3, 4] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + assert type(vec.temporal) is type(vec2.temporal) + + vec = vector.obj(x=1, y=2, theta=3, t=4) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3, 4] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + assert type(vec.temporal) is type(vec2.temporal) + + vec = vector.obj(pt=1, phi=2, eta=3, mass=4) + leaves, treedef = pytree.flatten(vec) + assert leaves == [1, 2, 3, 4] + vec2 = pytree.unflatten(treedef, leaves) + assert vec == vec2 + assert type(vec.azimuthal) is type(vec2.azimuthal) + assert type(vec.longitudinal) is type(vec2.longitudinal) + assert type(vec.temporal) is type(vec2.temporal) + + +def test_pytree_roundtrip_SoA(): + vec = vector.obj(x=1, y=1) * np.ones(10) + flat, unravel = pytree.ravel(vec) + assert flat.shape == (20,) + assert (unravel(flat) == vec).all() + + vec = vector.obj(pt=1, eta=1, phi=1, mass=1) * np.ones(10) + flat, unravel = pytree.ravel(vec) + assert flat.shape == (40,) + assert (unravel(flat) == vec).all() + + +def test_pytree_roundtrip_AoS(): + vec = vector.array( + { + "x": np.ones(10), + "y": np.ones(10), + } + ) + flat, unravel = pytree.ravel(vec) + assert flat.shape == (20,) + assert flat.dtype == np.float64 + assert (unravel(flat) == vec).all() + + vec = vector.array( + { + "pt": np.ones(10), + "eta": np.ones(10), + "phi": np.ones(10), + "mass": np.ones(10), + } + ) + flat, unravel = pytree.ravel(vec) + assert flat.shape == (40,) + assert flat.dtype == np.float64 + assert (unravel(flat) == vec).all() + + +def test_run_once(): + # Calling register_pytree multiple times returns the same object + pytree2 = vector.register_pytree() + assert pytree is pytree2