diff --git a/docs/tutorials/vacf_testing_examples.ipynb b/docs/tutorials/vacf_testing_examples.ipynb new file mode 100644 index 0000000..af6b630 --- /dev/null +++ b/docs/tutorials/vacf_testing_examples.ipynb @@ -0,0 +1,286 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5bc2ea6e-8449-4160-9d7b-6181da09c487", + "metadata": {}, + "source": [ + "## Calculation Testing Examples" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "54857ea4-5d70-4a91-843f-dace861c3e0a", + "metadata": {}, + "outputs": [], + "source": [ + "# imports\n", + "import MDAnalysis as mda\n", + "from transport_analysis.velocityautocorr import VelocityAutocorr as VACF\n", + "\n", + "# needed for the test examples in this notebook\n", + "import numpy as np\n", + "\n", + "# test data for this example\n", + "from MDAnalysis.tests.datafiles import PRM_NCBOX, TRJ_NCBOX" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "8b59ebc1-20f7-497d-864f-75e8c3a17382", + "metadata": {}, + "outputs": [], + "source": [ + "# test data with velocities\n", + "u = mda.Universe(PRM_NCBOX, TRJ_NCBOX)\n", + "ag_vels = u.select_atoms(\"name O and resname WAT and resid 1-10\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "a1753ee3-7db9-405f-82eb-1486d23e7916", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "v:\n", + "[ 4.22902895e+01 -2.69143315e+00 6.98534787e-01 2.97003597e+00\n", + " 6.34654795e-01 1.23400236e+00 -2.78565552e+00 7.25828978e-01\n", + " -1.14432591e-02 -6.46827198e+00]\n", + "float64\n" + ] + } + ], + "source": [ + "# confirm that analysis runs for test data with velocities (PRM_NCBOX, TRJ_NCBOX)\n", + "v = VACF(ag_vels, fft=True)\n", + "v.run()\n", + "\n", + "print(\"v:\")\n", + "print(v.results.timeseries)\n", + "\n", + "# confirm analysis uses float64\n", + "print(v.results.timeseries.dtype)" + ] + }, + { + "cell_type": "markdown", + "id": "871875ee-0424-4869-86b9-3ec230063251", + "metadata": {}, + "source": [ + "### Unit Velocity Trajectory Tests" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "d1126928-b847-4b7b-bde3-0fd8d7f85a85", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "v_simple:\n", + "[ 8.5500000e+01 8.0000000e+01 7.3500000e+01 6.6000000e+01\n", + " 5.7500000e+01 4.8000000e+01 3.7500000e+01 2.6000000e+01\n", + " 1.3500000e+01 -1.0658141e-14]\n", + "float64\n", + "poly:\n", + "[85.5 80. 73.5 66. 57.5 48. 37.5 26. 13.5 0. ]\n", + "float64\n", + "v_window:\n", + "[85.5 80. 73.5 66. 57.5 48. 37.5 26. 13.5 0. ]\n", + "float64\n" + ] + } + ], + "source": [ + "# set unit velocity trajectory to have 10 frames\n", + "NSTEP = 10\n", + "\n", + "# Step trajectory of unit velocities for a single particle i.e. v = 0 at t = 0,\n", + "# v = 1 at t = 1, v = 2 at t = 2, etc. for all components x, y, z\n", + "v_test = np.arange(NSTEP)\n", + "velocities = np.vstack([v_test, v_test, v_test]).T\n", + "# NSTEP frames x 1 atom x 3 dimensions, where NSTEP = 10\n", + "velocities_reshape = velocities.reshape([NSTEP, 1, 3])\n", + "step_vtraj = mda.Universe.empty(1, n_frames=NSTEP, velocities=True)\n", + "for i, ts in enumerate(step_vtraj.trajectory):\n", + " step_vtraj.atoms.velocities = velocities_reshape[i]\n", + "\n", + "# Expected VACF results for unit velocity trajectory\n", + "# At time t, VACF is:\n", + "# sum_{x=0}^{N - 1 - t} x*(x + t) * n_dim / n_frames\n", + "# n_dim = 3 (typically) and n_frames = total_frames - t\n", + "def characteristic_poly(last, n_dim, first=0, step=1):\n", + " diff = last - first\n", + " frames_used = diff // step + 1 if diff % step != 0 else diff / step\n", + " frames_used = int(frames_used)\n", + " result = np.zeros(frames_used)\n", + " for t in range(first, last, step):\n", + " sum = 0\n", + " sum = np.dtype(\"float64\").type(sum)\n", + " lagtime = t - first\n", + " for x in range(first, (last - lagtime), step):\n", + " sum += x * (x + lagtime)\n", + " current_index = int(lagtime / step)\n", + " vacf = sum * n_dim / (frames_used - current_index)\n", + " result[current_index] = vacf\n", + " return result\n", + "\n", + "# the following test examples set n_dim to 3\n", + "# all of the following calculations should give the same result\n", + "# and use float64\n", + "\n", + "poly = characteristic_poly(NSTEP, 3)\n", + "\n", + "# result for analysis run on unit velocity trajectory WITH FFT\n", + "v_simple = VACF(step_vtraj.atoms, fft=True)\n", + "v_simple.run()\n", + "print(\"v_simple:\")\n", + "print(v_simple.results.timeseries)\n", + "print(v_simple.results.timeseries.dtype)\n", + "\n", + "# result from characteristic_poly (expected VACF results)\n", + "print(\"poly:\")\n", + "print(poly)\n", + "print(poly.dtype)\n", + "\n", + "# result for analysis run on unit velocity trajectory WITHOUT FFT\n", + "v_window = VACF(step_vtraj.atoms, fft=False)\n", + "v_window.run()\n", + "print(\"v_window:\")\n", + "print(v_window.results.timeseries)\n", + "print(v_window.results.timeseries.dtype)" + ] + }, + { + "cell_type": "markdown", + "id": "448dc507-3cfd-4a9b-8f6a-6271aa047dc2", + "metadata": {}, + "source": [ + "## Plotting Examples" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "091388e6-cf1b-4de5-a693-00ccddbe2eda", + "metadata": {}, + "outputs": [], + "source": [ + "# if using only the terminal, the plots may not appear automatically like they do here\n", + "# you will then need the following:\n", + "# import matplotlib.pyplot as plt\n", + "\n", + "# then after running the desired plotting function, run:\n", + "# plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "08ff9ecb-3969-4545-a6d5-fae778097705", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAGwCAYAAAC+Qv9QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYnElEQVR4nO3dd1xV9eMG8OdcxmUjIvcCgjIdiHtbuTE1SzPL1BypDUdqVpotrW9q2i8zR5aWq1zlKFvmxpWJGqm4WCKiLJG97z2/P5CrhIN777kcLvd5v168vnIuHB7jqzx+zmcIoiiKICIiIjJjCrkDEBERERmLhYaIiIjMHgsNERERmT0WGiIiIjJ7LDRERERk9lhoiIiIyOyx0BAREZHZs5Y7QHXRarW4fv06nJ2dIQiC3HGIiIioCkRRRE5ODry9vaFQ3H8cxmIKzfXr1+Hr6yt3DCIiIjJAYmIifHx87vu6xRQaZ2dnAGX/QVxcXGROQ0RERFWRnZ0NX19f3c/x+7GYQlP+mMnFxYWFhoiIyMw8bLoIJwUTERGR2WOhISIiIrPHQkNERERmj4WGiIiIzB4LDREREZk9FhoiIiIyeyw0REREZPZYaIiIiMjssdAQERGR2WOhISIiIrPHQkNERERmj4WGiIiIzB4LjZE0WhFXb+YjJbtQ7ihEREQWi4XGSB/8fA5dPz2A748nyB2FiIjIYrHQGCnAwwkAEJOaK3MSIiIiy8VCY6QgVVmhiWahISIikg0LjZGCbxeaK+l5KNFoZU5DRERkmVhojOTlagdHWyuUakUk3MyXOw4REZFFYqExkiAICFSVz6PJkTkNERGRZWKhkUCQihODiYiI5MRCIwEWGiIiInmx0EggyIMrnYiIiORkbcwnl5SUIDk5Gfn5+fDw8EDdunWlymVWgtXOAIDYtFxotSIUCkHmRERERJZF7xGa3NxcfP311+jevTtcXV3h5+eHkJAQeHh4oGHDhnjppZcQERFhiqw1lq+bPWytFCgs0SIps0DuOERERBZHr0Lz+eefw8/PD6tWrULPnj2xfft2REZG4tKlS/jrr78we/ZslJaWIiwsDH379kV0dLSpctco1lYKBHg4AuA8GiIiIjno9cjp2LFjOHDgAJo3b37P1zt06ICxY8fiq6++wrfffovw8HAEBwdLErSmC1Q54WJyDmJSc9GjiUruOERERBZFrxGaH3/88b5l5m5KpRITJ07E+PHj9Q40f/58CIKAadOm6a6Joog5c+bA29sb9vb26N69O6KiovS+tyndmRjMvWiIiIiqW41a5RQREYGVK1eiRYsWFa4vXLgQixYtwrJlyxAREQFPT0+EhYUhJ6fmlIdgNZduExERycWgQnPjxg18//33+P3331FcXFzhtby8PHz00Ud63zM3NxcjRozAqlWr4ObmprsuiiIWL16Md999F4MHD0ZoaCjWrVuH/Px8bNy40ZD4JnH3IZWiKMqchoiIyLLoXWgiIiIQEhKCSZMmYciQIQgNDa3w+Cc3Nxcffvih3kEmTZqEJ554Ar17965wPT4+HsnJyejTp4/umlKpRLdu3XDs2LH73q+oqAjZ2dkV3kzJv54jFAKQU1iKtJwik34tIiIiqkjvQvPOO+9g8ODBuHXrFlJSUhAWFoZu3brhn3/+MTjE5s2bcfr0acyfP7/Sa8nJyQAAtVpd4bparda9di/z58+Hq6ur7s3X19fgfFWhtLZCQ3eudCIiIpKD3oXm1KlTmDlzJhQKBZydnbF8+XLMmDEDvXr1Mmj/mcTEREydOhXff/897Ozs7vtxglBxszpRFCtdu9usWbOQlZWle0tMTNQ7m74CuWMwERGRLAzaKbiwsLDC+zNmzIBCoUCfPn2wevVqve516tQppKamom3btrprGo0Ghw4dwrJly3Dp0iUAZSM1Xl5euo9JTU2tNGpzN6VSCaVSqVcWYwWrnbD3QgpHaIiIiKqZ3oUmNDQUx44dq7QS6c0334Qoihg2bJhe9+vVqxfOnj1b4dqLL76IJk2aYObMmQgICICnpyf27NmD1q1bAwCKi4sRHh6OBQsW6BvfpMqXbrPQEBERVS+9C82oUaMQHh6OV199tdJrb731FkRRxIoVK6p8P2dnZ4SGhla45ujoCHd3d931adOmYd68eQgODkZwcDDmzZsHBwcHDB8+XN/4JnX3SiciIiKqPnoXmvHjxz9ww7wZM2ZgxowZRoW61z0LCgowceJE3Lp1Cx07dsTu3bvh7Ows6dcxVuDtQpOeW4TM/GLUcbCVOREREZFlEEQL2TQlOzsbrq6uyMrKgouLi8m+Tpf5+3A9qxBbX+2Mdn6Wefo4ERGRVKr689vgnYKPHDli6KfWauWjNJxHQ0REVH0MKjS//vorBg8eLHWWWiFYVfYYjIWGiIio+uhdaL777jsMHz68Rh07UJNwYjAREVH106vQLF68GC+//DI2bNhQ6YgCKsNDKomIiKqfXqucpk+fjqVLl+LJJ580VR6zV74XTVJmAfKLS+Fga9DehURERKQHvUZounTpguXLl+PmzZumymP23Bxt4e5Ytlw7NjVP5jRERESWQa9Cs3fvXvj7+yMsLMzkp1ebs/J5NDFpOTInISIisgx6FRo7Ozvs3LkTISEhePzxx02VyezpJgancB4NERFRddB7lZOVlRW+//57dOjQwRR5aoUg7kVDRERUrQzeWO+LL76QMketotuLJo2FhoiIqDoYXGjo/spHaBJu5qO4VCtzGiIiotrP4EKza9euCscfLF++HK1atcLw4cNx69YtScKZK7WLEs5Ka2i0Iq7c5EonIiIiUzO40Lz11lu6lU5nz57FG2+8gf79+yMuLg7Tp0+XLKA5EgSBZzoRERFVI4N3fYuPj0dISAgAYNu2bRgwYADmzZuH06dPo3///pIFNFdBKidEJmaWrXRqLncaIiKi2s3gERpbW1vk5+cDKNufpk+fPgCAunXrco8aAMG6vWg4QkNERGRqBo/QPProo5g+fToeeeQRnDhxAlu2bAEAXL58GT4+PpIFNFd39qLh5npERESmZvAIzbJly2BtbY2tW7dixYoVqF+/PgDgjz/+QN++fSULaK7KC01ceh40WlHmNERERLWbwSM0DRo0wK+//lrp+ueff25UoNrCx80BSmsFikq1uHYrHw3dHeWOREREVGsZdRS0RqPBjh07cOHCBQiCgCZNmmDQoEGwtuYJ01YKAQEeTrhwIxvRKbksNERERCZkcPM4d+4cnnrqKaSkpKBx48YAyubPeHh4YOfOnWjenEt7glVlhSYmLRe9oZY7DhERUa1l8Bya8ePHIzQ0FNeuXcPp06dx+vRpJCYmokWLFnj55ZelzGi2eEglERFR9TB4hObff//FyZMn4ebmprvm5uaGuXPnon379pKEM3dBXLpNRERULQweoWncuDFSUlIqXU9NTUVQUJBRoWqL8r1oYlNzIYpc6URERGQqBheaefPmYcqUKdi6dSuuXbuGa9euYevWrZg2bRoWLFiA7Oxs3ZulaujuCCuFgNyiUiRnF8odh4iIqNYy+JHTgAEDAADPPfccBEEAAN0oxJNPPql7XxAEaDQaY3OaJVtrBRq6OyAuLQ8xqbnwcrWXOxIREVGtZHChOXDggJQ5aq1glZOu0DwW7CF3HCIiolrJ4ELTrVs3KXPUWkEqJ/wZlYJonrpNRERkMgbPoaGqCVY5AwBiWGiIiIhMhoXGxHRLt1loiIiITIaFxsQCPMqOPMjIK0ZGXrHMaYiIiGonvQtNbi5HGvThYGsNH7ey1U0cpSEiIjINvQtNvXr10K9fP6xYsQJJSUmmyFTr6I5ASM2ROQkREVHtpHehuXTpEvr3749t27YhICAA7du3x//+9z+cOXPGFPlqhSAPzqMhIiIyJb0LTcOGDfHaa69h7969SE1NxfTp0xEVFYWuXbvC398fU6dOxf79+y12M717CVaz0BAREZmSUZOCXV1dMWzYMGzevBnp6en4+uuvodVq8eKLL8LDwwMbNmyQKqdZ40onIiIi0zJ4Y71KN7K2Rp8+fdCnTx8sXboU//zzD0pLS6W6vVkL8ijbi+ZGViFyCkvgbGcjcyIiIqLaRbJC81+tW7c21a3NjquDDTyclUjLKUJsWh5a+daROxIREVGtwn1oqgknBhMREZkOC0014cRgIiIi02GhqSZ3JgZzLxoiIiKpsdBUEz5yIiIiMh2DCs2///6Ljz/+GF9++SXS09MrvJadnY2xY8dKEq42Cbr9yOlqRj4KS7hHDxERkZT0LjS7d+9Ghw4dsHnzZixYsABNmzbFgQMHdK8XFBRg3bp1koasDTyclHCxs4ZWBOLT8+SOQ0REVKvoXWjmzJmDN998E+fOncOVK1cwY8YMPPXUU9i1a5cp8tUagiAgWF22Hw0fOxEREUlL70ITFRWle6QkCALeeustrFy5EkOGDMEvv/wiecDapHweTTQLDRERkaT03lhPqVQiMzOzwrVhw4ZBoVDg+eefx2effSZVtlqnfKVTLAsNERGRpPQuNK1atcKBAwfQtm3bCteHDh0KrVaL0aNHSxautgniXjREREQmoXehmTBhAg4dOnTP14YNGwYAWLlypXGpaqnyR05x6bko1WhhbcVV80RERFIQRFEU5Q5RHbKzs+Hq6oqsrCy4uLjIkkGrFdFs9p8oKNFg/xvdEHC74BAREdG9VfXnN4cIqpFCISBQ5QiAE4OJiIikZHCh+eGHH6TMYTG4YzAREZH0DCo0X331FaZNmyZxFMtQvhcNVzoRERFJR+9JwR9//DEWL16Mffv2mSJPrRfIvWiIiIgkp1ehmTZtGtasWYPdu3ejZcuWpspUq+n2oknLhVYrQqEQZE5ERERk/vR65LRkyRJ89tln6Nixo6ny1HoN3R1gYyUgv1iDG9mFcschIiKqFfQqNM888wxmz56NuLg4U+Wp9WysFPBzv73SKSVH5jRERES1g16F5ocffsCAAQPQq1cvJCUlmSpTrRfMHYOJiIgkpVehEQQBX3/9NZ5//nn07NnTVJlqPS7dJiIikpZBy7bnz5+PCRMmSBZixYoVaNGiBVxcXODi4oLOnTvjjz/+0L0uiiLmzJkDb29v2Nvbo3v37oiKipLs61e3QBULDRERkZQM3lhPyn1ofHx88Mknn+DkyZM4efIkevbsiYEDB+pKy8KFC7Fo0SIsW7YMERER8PT0RFhYGHJyzHMOSrCqbC+a6NRcWMjJE0RERCYl+dEHJSUl+PHHH/X6nCeffBL9+/dHo0aN0KhRI8ydOxdOTk44fvw4RFHE4sWL8e6772Lw4MEIDQ3FunXrkJ+fj40bN0odv1oEeDhCEICsghKk5xbLHYeIiMjsSVZo/vnnH0yZMgWenp4YMWKEwffRaDTYvHkz8vLy0LlzZ8THxyM5ORl9+vTRfYxSqUS3bt1w7Nix+96nqKgI2dnZFd5qCjsbK/i6OQDgYyciIiIpGFVobt68iS+++AKtWrVCu3btcP78eSxYsADXr1/X+15nz56Fk5MTlEolXn31VezYsQMhISFITk4GAKjV6gofr1arda/dy/z58+Hq6qp78/X11TuTKQWXz6NJY6EhIiIyVpUKTXR0NJ5//nlMmjQJ6enp+PXXXzF48GB4e3tj3bp1GDlyJARBwJIlSzB+/HjUq1dP7yCNGzdGZGQkjh8/jgkTJmD06NE4f/687nVBqLijriiKla7dbdasWcjKytK9JSYm6p3JlMp3DI7hXjRERERGq9LRBy+88AKGDBkCPz8/eHl5oX79+hg2bBg+/vhjhISEAADefvtto4LY2toiKCgIANCuXTtERETgiy++wMyZMwEAycnJ8PLy0n18ampqpVGbuymVSiiVSqMymVIQR2iIiIgkU6URmrS0NLRr1w4dOnSAKIro2rUrevbsiaZNm5osmCiKKCoqgr+/Pzw9PbFnzx7da8XFxQgPD0eXLl1M9vVNrbzQRKew0BARERmrSiM0n3/+OcaMGQN3d3ds3rwZx44dwwsvvABra2sMGzYMI0aMeODjn4d555130K9fP/j6+iInJwebN2/GwYMHsWvXLgiCgGnTpmHevHkIDg5GcHAw5s2bBwcHBwwfPtzgrym38r1oUnOKkF1YAhc7G5kTERERma8qFZqBAwdi4MCBuveHDBmChQsX4pdffsHatWvRsWNHlJaWYt26dRg3bhwaNWqkV4iUlBSMHDkSN27cgKurK1q0aIFdu3YhLCwMADBjxgwUFBRg4sSJuHXrFjp27Ijdu3fD2dlZr69Tk7jY2cDTxQ7J2YWISc1FmwZuckciIiIyW4Iowc5uqampWL9+PdauXYvz588jNDQUZ86ckSKfZLKzs+Hq6oqsrCy4uLjIHQcA8MI3f+NITDoWPtMCz7WvWauwiIiIaoKq/vyWZB8alUqFN998E+fOncPx48fx6KOPSnHbWo8Tg4mIiKRRpUdO+ujQoQM6dOgg9W1rpTsTg7l0m4iIyBh6jdD07dv3gbvzlsvJycGCBQuwfPlyg4NZAo7QEBERSUOvEZpnn30Wzz33HJydnfHUU0+hXbt28Pb2hp2dHW7duoXz58/jyJEj+P333zFgwAB8+umnpspdK5TvFnztVgEKijWwt7WSOREREZF50qvQjBs3DiNHjsTWrVuxZcsWrFq1CpmZmQDKdvINCQnB448/jlOnTqFx48amyFuruDsp4eZgg1v5JYhNy0VofVe5IxEREZklvefQ2NraYvjw4bo9YLKyslBQUAB3d3fY2HAvFX0FqZwQceUWCw0REZERjF7l5OrqCk9PT5YZAwWpyvbS4anbREREhpNk2TYZjkcgEBERGY+FRmZc6URERGQ8FhqZla90upKehxKNVuY0RERE5omFRmZernZwtLVCqVZEws08ueMQERGZJaN3Ci4uLkZqaiq02oqjCw0aNDD21hZBEAQEqZzw77UsxKTm6iYJExERUdUZXGiio6MxduzYSjsHi6IIQRCg0WiMDmcpAm8XmuiUXPQNlTsNERGR+TG40IwZMwbW1tb49ddf4eXlBUEQpMxlUTgxmIiIyDgGF5rIyEicOnUKTZo0kTKPRQrmXjRERERGMXhScEhICNLT06XMYrHKR2hi03Kh1YoypyEiIjI/BheaBQsWYMaMGTh48CBu3ryJ7OzsCm9Udb5u9rC1UqCwRIukzAK54xAREZkdgx859e7dGwDQq1evCtc5KVh/1lYKBHg44mJyDqJTc+Bb10HuSERERGbF4EJz4MABKXNYvECVEy4m5yAmNRc9m6jljkNERGRWDC403bp1kzKHxSvfMZgTg4mIiPRn1MZ6mZmZ+Pbbb3HhwgUIgoCQkBCMHTsWrq6uUuWzGLpDKlloiIiI9GbwpOCTJ08iMDAQn3/+OTIyMpCeno5FixYhMDAQp0+fljKjRQi6a4RGFLnSiYiISB8Gj9C8/vrreOqpp7Bq1SpYW5fdprS0FOPHj8e0adNw6NAhyUJaAv96jlAIQE5hKdJyiqBysZM7EhERkdkwaoRm5syZujIDANbW1pgxYwZOnjwpSThLorS2QkN3RwB87ERERKQvgwuNi4sLrl69Wul6YmIinJ15wKIhAj04MZiIiMgQBheaoUOHYty4cdiyZQsSExNx7do1bN68GePHj8ewYcOkzGgxgtXlE4NzZE5CRERkXgyeQ/N///d/EAQBo0aNQmlpKQDAxsYGEyZMwCeffCJZQEsSxBEaIiIigxhcaGxtbfHFF19g/vz5iI2NhSiKCAoKgoMDd7k1VPkITUxqnsxJiIiIzItR+9AAgIODA5o3by5FFotXPocmPbcImfnFqONgK3MiIiIi86BXoRk8eDDWrl0LFxcXDB48+IEfu337dqOCWSJHpTW8Xe1wPasQMam5aOdXV+5IREREZkGvQuPq6gpBEACUrXIq/zVJJ0jtjOtZhYhmoSEiIqoyvQrNmjVrdL9eu3at1FkIZRODD11O48RgIiIiPRi8bLtnz57IzMysdD07Oxs9e/Y0JpNFC+IhlURERHozuNAcPHgQxcXFla4XFhbi8OHDRoWyZHdWOrHQEBERVZXeq5zOnDmj+/X58+eRnJyse1+j0WDXrl2oX7++NOksUPleNEmZBcgrKoWj0uiFaERERLWe3j8tW7VqBUEQIAjCPR8t2dvbY+nSpZKEs0Rujrao52SL9NxixKXlobmPq9yRiIiIajy9C018fDxEUURAQABOnDgBDw8P3Wu2trZQqVSwsrKSNKSlCfRwQnpuBqJTc1hoiIiIqkDvQtOwYUMAgFarlTwMlQlSOeHv+AzOoyEiIqoigycFz58/H6tXr650ffXq1ViwYIFRoSxdsKr8kEoWGiIioqowuNB8/fXXaNKkSaXrzZo1w1dffWVUKEsXpHIGAMSy0BAREVWJwYUmOTkZXl5ela57eHjgxo0bRoWydOV70SRk5KOoVCNzGiIioprP4ELj6+uLo0ePVrp+9OhReHt7GxXK0qldlHBWWkOjFXElPV/uOERERDWewZucjB8/HtOmTUNJSYlu+fa+ffswY8YMvPHGG5IFtESCICBQ5YTIxEzEpOaisaez3JGIiIhqNIMLzYwZM5CRkYGJEyfqdgy2s7PDzJkzMWvWLMkCWqrguwoNERERPZjBhUYQBCxYsADvv/8+Lly4AHt7ewQHB0OpVEqZz2IF6VY65cichIiIqOYzel99JycntG/fXoosdBceUklERFR1BheavLw8fPLJJ9i3bx9SU1MrbbQXFxdndDhLFnx76XZceh40WhFWCkHmRERERDWXUZOCw8PDMXLkSHh5eUEQ+ANXSvXd7KG0VqCoVIvEjHz41XOUOxIREVGNZXCh+eOPP/Dbb7/hkUcekTIP3WalEBDg4YQLN7IRk5rLQkNERPQABu9D4+bmhrp160qZhf6DRyAQERFVjcGF5n//+x8++OAD5Odz4zdT4cRgIiKiqjH4kdNnn32G2NhYqNVq+Pn5wcbGpsLrp0+fNjqcpSsfoYlJY6EhIiJ6EIMLzaBBgySMQfdSPkITm5oLURQ58ZqIiOg+DC40s2fPljIH3UNDd0dYKQTkFpUiObsQXq72ckciIiKqkQyeQ0OmZ2utgJ+7AwAgOoWPnYiIiO7H4EKjUChgZWV13zeSBicGExERPZzBj5x27NhR4f2SkhL8888/WLduHT788EOjg1GZIJUT/oxK4cRgIiKiBzC40AwcOLDStSFDhqBZs2bYsmULxo0bZ1QwKlN+BEIMHzkRERHdl+RzaDp27Ii9e/fq9Tnz589H+/bt4ezsDJVKhUGDBuHSpUsVPkYURcyZMwfe3t6wt7dH9+7dERUVJWX0GimIS7eJiIgeStJCU1BQgKVLl8LHx0evzwsPD8ekSZNw/Phx7NmzB6WlpejTpw/y8vJ0H7Nw4UIsWrQIy5YtQ0REBDw9PREWFoacnBwpfws1TqCHEwQByMgrxs3cIrnjEBER1UgGP3Jyc3OrsC+KKIrIycmBg4MDvv/+e73utWvXrgrvr1mzBiqVCqdOnULXrl0hiiIWL16Md999F4MHDwYArFu3Dmq1Ghs3bsQrr7xi6G+jxrO3tUL9Ova4dqsAMam5cHdSyh2JiIioxjG40CxevLjC+wqFAh4eHujYsSPc3NyMCpWVlQUAurOi4uPjkZycjD59+ug+RqlUolu3bjh27Ng9C01RURGKiu6MaGRnZxuVSU5BKqeyQpOWi44B7nLHISIiqnH0LjSrV6/GiBEjMHr0aFPkgSiKmD59Oh599FGEhoYCAJKTkwEAarW6wseq1WokJCTc8z7z58+vNautglVOOHgpjXvREBER3Yfec2heeukl3QgKAHh7e+PKlSuSBZo8eTLOnDmDTZs2VXrtv1v/P+g4gFmzZiErK0v3lpiYKFnG6qY7AoETg4mIiO5J7xEaURQrvJ+TkwOtVitJmNdeew07d+7EoUOHKkws9vT0BFA2UuPl5aW7npqaWmnUppxSqYRSWTvmm3BzPSIiogerEUcfiKKIyZMnY/v27di/fz/8/f0rvO7v7w9PT0/s2bNHd624uBjh4eHo0qVLdcetdkEeZXvR3MgqRE5hicxpiIiIah69R2gEQajwmOe/7xti0qRJ2LhxI37++Wc4Ozvr5sy4urrC3t4egiBg2rRpmDdvHoKDgxEcHIx58+bBwcEBw4cPN+prmwNXBxt4OCuRllOE2LQ8tPKtI3ckIiKiGsWgR06NGjXSlZjc3Fy0bt0aCkXFwZ6MjIwq33PFihUAgO7du1e4vmbNGowZMwYAMGPGDBQUFGDixIm4desWOnbsiN27d8PZ2Vnf34JZClY5IS2nCNEpOSw0RERE/6F3oVmzZo3kIf47L+deBEHAnDlzMGfOHMm/vjkIUjnhWOxN7hhMRER0D3oXGlMt16YH06104sRgIiKiSmrEpGB6uPJCE81CQ0REVAkLjZkoLzSJGfkoLNHInIaIiKhmYaExEx5OSrjYWUMrAvHpeQ//BCIiIgvCQmMmBEFAsLpsRRcfOxEREVXEQmNGgjy4YzAREdG9GHzatkajwdq1a7Fv3z6kpqZWOv5g//79RoejioLV5YUmR+YkRERENYvBhWbq1KlYu3YtnnjiCYSGhhq9WzA9XCDPdCIiIrongwvN5s2b8cMPP6B///5S5qEHKH/kFJ+eh1KNFtZWfGJIREQEGDGHxtbWFkFBQVJmoYeoX8ce9jZWKNGISMjIlzsOERFRjWFwoXnjjTfwxRdfVOnYApKGQiEgUOUIgI+diIiI7mbwI6cjR47gwIED+OOPP9CsWTPY2NhUeH379u1Gh6PKgjyccC4pGzGpuXi8mdxpiIiIagaDC02dOnXw9NNPS5mFqqB8LxqO0BAREd1hcKExxanb9HCB3IuGiIioEoMLTbm0tDRcunQJgiCgUaNG8PDwkCIX3cedvWhyodWKUCi4XJ6IiMjgScF5eXkYO3YsvLy80LVrVzz22GPw9vbGuHHjkJ/PFTim0rCuA2ysBBSUaHA9q0DuOERERDWCwYVm+vTpCA8Pxy+//ILMzExkZmbi559/Rnh4ON544w0pM9JdrK0U8HPnSiciIqK7GVxotm3bhm+//Rb9+vWDi4sLXFxc0L9/f6xatQpbt26VMiP9x92PnYiIiMiIQpOfnw+1Wl3pukql4iMnE+MhlURERBUZXGg6d+6M2bNno7CwUHetoKAAH374ITp37ixJOLq38jOdolloiIiIABixyumLL75A37594ePjg5YtW0IQBERGRsLOzg5//vmnlBnpP4JVd/aiEUWRB4MSEZHFM7jQhIaGIjo6Gt9//z0uXrwIURTx/PPPY8SIEbC3t5cyI/1HgIcjBAHIKihBem4xPJyVckciIiKSlVH70Njb2+Oll16SKgtVkZ2NFRrUdUDCzXxEp+aw0BARkcXTq9Ds3LkT/fr1g42NDXbu3PnAj33qqaeMCkYPFuThhISb+YhNzUWXwHpyxyEiIpKVXoVm0KBBSE5OhkqlwqBBg+77cYIgQKPRGJuNHiBI5YR9F1O50omIiAh6FhqtVnvPX1P1C+JKJyIiIh2Dl22vX78eRUVFla4XFxdj/fr1RoWihysvNByhISIiMqLQvPjii8jKyqp0PScnBy+++KJRoejhyveiSc0pQlZBicxpiIiI5GVwobnf/ifXrl2Dq6urUaHo4VzsbODpYgeAozRERER6L9tu3bo1BEGAIAjo1asXrK3v3EKj0SA+Ph59+/aVNCTdW5DKCcnZhYhNzUXbhm5yxyEiIpKN3oWmfHVTZGQkHn/8cTg5Oeles7W1hZ+fH5555hnJAtL9BamccCQmHdGpOXJHISIikpXehWb27NkAAD8/PwwdOhR2dnaSh6Kq4cRgIiKiMgbvFDx69Ggpc5ABdIUmjYWGiIgsm8GFRqPR4PPPP8cPP/yAq1evori4uMLrGRkZRoejBwu+XWiu3SpAQbEG9rZWMiciIiKSh8GrnD788EMsWrQIzz33HLKysjB9+nQMHjwYCoUCc+bMkTAi3Y+7kxJuDjYQRSCWozRERGTBDC40GzZswKpVq/Dmm2/C2toaw4YNwzfffIMPPvgAx48flzIjPQDn0RARERlRaJKTk9G8eXMAgJOTk26TvQEDBuC3336TJh09VJDKGQALDRERWTaDC42Pjw9u3LgBAAgKCsLu3bsBABEREVAqldKko4fiCA0REZERhebpp5/Gvn37AABTp07F+++/j+DgYIwaNQpjx46VLCA9WLDukEruRUNERJbL4FVOn3zyie7XQ4YMgY+PD44dO4agoCA89dRTkoSjhysfoUm4mY/iUi1srQ3uqERERGbL4ELzX506dUKnTp2kuh1VkZerHRxtrZBXrEHCzTwEq53ljkRERFTt9Co0O3furPLHcpSmegiCgCCVE/69loWY1FwWGiIiskh6FZryc5weRhAEaDQaQ/KQAQLvKjRERESWSK9Co9VqTZWDjBCkmxjMQkNERJZJkhmkhYWFUtyGDBTMvWiIiMjCGVxoNBoN/ve//6F+/fpwcnJCXFwcAOD999/Ht99+K1lAerjyEZrYtFxotKLMaYiIiKqfwYVm7ty5WLt2LRYuXAhbW1vd9ebNm+Obb76RJBxVja+bPWytFSgq1SLpVoHccYiIiKqdwYVm/fr1WLlyJUaMGAErqzunPLdo0QIXL16UJBxVjbWVAgH1HAEAMWncYI+IiCyPwYUmKSkJQUFBla5rtVqUlJQYFYr0F1g+MTiF82iIiMjyGFxomjVrhsOHD1e6/uOPP6J169ZGhSL9BfNMJyIismAG7xQ8e/ZsjBw5EklJSdBqtdi+fTsuXbqE9evX49dff5UyI1WB7pDKNBYaIiKyPAaP0Dz55JPYsmULfv/9dwiCgA8++AAXLlzAL7/8grCwMCkzUhXoCk1KLkSRK52IiMiyGDRCU1pairlz52Ls2LEIDw+XOhMZwL+eIxQCkFNUitScIqhd7OSOREREVG0MGqGxtrbGp59+yuMNahCltRUaut9e6cR5NEREZGEMfuTUu3dvHDx4UMIoZCzdEQgpXLpNRESWxeBJwf369cOsWbNw7tw5tG3bFo6OjhVe52nb1S9I5YQ951M4MZiIiCyOwYVmwoQJAIBFixZVeo2nbcsjyIN70RARkWUyuNDw5O2aJ1h950wnIiIiS2LQHJrS0lJYW1vj3LlzkoQ4dOgQnnzySXh7e0MQBPz0008VXhdFEXPmzIG3tzfs7e3RvXt3REVFSfK1a5PA2yM06bnFuJVXLHMaIiKi6mPwKqeGDRtK9lgpLy8PLVu2xLJly+75+sKFC7Fo0SIsW7YMERER8PT0RFhYGHJyOPn1bo5Ka3i7li3X5jwaIiKyJAavcnrvvfcwa9YsZGRkGB2iX79++PjjjzF48OBKr4miiMWLF+Pdd9/F4MGDERoainXr1iE/Px8bN2687z2LioqQnZ1d4c0SBKmdAXDpNhERWRaDC82SJUtw+PBheHt7o3HjxmjTpk2FN6nEx8cjOTkZffr00V1TKpXo1q0bjh07dt/Pmz9/PlxdXXVvvr6+kmWqyconBrPQEBGRJTF4UvCgQYMkjHF/ycnJAAC1Wl3hulqtRkJCwn0/b9asWZg+fbru/ezsbIsoNeUTg6NZaIiIyIIYdThldRIEocL7oihWunY3pVIJpVJp6lg1TvnmerEsNEREZEEMLjTlTp06hQsXLkAQBISEhKB169ZS5NLx9PQEUDZS4+XlpbuemppaadSG7jxySsosQF5RKRyVRn+LiYiIajyD59CkpqaiZ8+eaN++PaZMmYLJkyejbdu26NWrF9LS0iQL6O/vD09PT+zZs0d3rbi4GOHh4ejSpYtkX6e2cHO0RT0nWwDcj4aIiCyHwYXmtddeQ3Z2NqKiopCRkYFbt27h3LlzyM7OxpQpU/S6V25uLiIjIxEZGQmgbCJwZGQkrl69CkEQMG3aNMybNw87duzAuXPnMGbMGDg4OGD48OGGxq/VAjkxmIiILIzBzyN27dqFvXv3omnTprprISEhWL58eYUVSVVx8uRJ9OjRQ/d++WTe0aNHY+3atZgxYwYKCgowceJE3Lp1Cx07dsTu3bvh7OxsaPxaLUjlhL/jMzgxmIiILIZRRx/Y2NhUum5jY6P3sQjdu3eHKIr3fV0QBMyZMwdz5szRN6ZFClZxhIaIiCyLwY+cevbsialTp+L69eu6a0lJSXj99dfRq1cvScKRYYJU3FyPiIgsi8GFZtmyZcjJyYGfnx8CAwMRFBQEf39/5OTkYOnSpVJmJD2V70WTcDMPRaU89ZyIiGo/gx85+fr64vTp09izZw8uXrwIURQREhKC3r17S5mPDKByVsJZaY2colJcSc9HY0/ONSIiotrN6E1KwsLCEBYWJkUWkoggCAhUOSEyMRPRqTksNEREVOsZ/MhpypQpWLJkSaXry5Ytw7Rp04zJRBLgxGAiIrIkBheabdu24ZFHHql0vUuXLti6datRoch4QSw0RERkQQwuNDdv3oSrq2ul6y4uLkhPTzcqFBmPhYaIiCyJwYUmKCgIu3btqnT9jz/+QEBAgFGhyHjBt5dux6XnoVSj375ARERE5sbgScHTp0/H5MmTkZaWhp49ewIA9u3bh88++wyLFy+WKh8ZqL6bPZTWChSVapF4qwD+9RzljkRERGQyBheasWPHoqioCHPnzsX//vc/AICfnx9WrFiBUaNGSRaQDGOlEBDo4YTzN7IRk5rLQkNERLWawY+cAGDChAm4du0aUlJSkJ2djbi4OJaZGoTzaIiIyFIYvQ9NWloaLl26BEEQ0LhxY9SrV0+KXCSB8kITnZojcxIiIiLTMniEJi8vD2PHjoWXlxe6du2Kxx57DF5eXhg3bhzy8/OlzEgGKt+LJpYjNEREVMsZXGimT5+O8PBw/PLLL8jMzERmZiZ+/vlnhIeH44033pAyIxno7kdODzrNnIiIyNwZ/Mhp27Zt2Lp1K7p376671r9/f9jb2+O5557DihUrpMhHRmjo7ggrhYC8Yg1uZBXCu4693JGIiIhMwuARmvz8fKjV6krXVSoVHznVELbWCvi5OwDgxGAiIqrdDC40nTt3xuzZs1FYWKi7VlBQgA8//BCdO3eWJBwZ787EYBYaIiKqvQx+5LR48WL069cPPj4+aNmyJQRBQGRkJOzs7PDnn39KmZGMEKxyxp9RKRyhISKiWs3gQtO8eXNER0fj+++/x8WLFyGKIp5//nmMGDEC9vacq1FTBHGlExERWQCDC82hQ4fQpUsXvPTSSxWul5aW4tChQ+jatavR4ch43IuGiIgsgcFzaHr06IGMjIxK17OystCjRw+jQpF0Aj2cIAjArfwS3MwtkjsOERGRSRhcaERRhCAIla7fvHkTjo48N6imsLe1Qv3by7U5MZiIiGorvR85DR48GAAgCALGjBkDpVKpe02j0eDMmTPo0qWLdAnJaEEqJ1y7VYCY1Fx0CnCXOw4REZHk9C40rq6uAMpGaJydnStMALa1tUWnTp0qzasheQWrnHDwUhpXOhERUa2ld6FZs2YNAMDPzw9vvvkmHy+ZAZ66TUREtZ3Bq5xmz54tZQ4yoSCVMwAWGiIiqr0MLjT+/v73nBRcLi4uztBbk8TKR2iSswuRU1gCZzsbmRMRERFJy+BCM23atArvl5SU4J9//sGuXbvw1ltvGZuLJORqbwMPZyXScooQk5qL1g3c5I5EREQkKYMLzdSpU+95ffny5Th58qTBgcg0glVOLDRERFRrGbwPzf3069cP27Ztk/q2ZCRODCYiotpM8kKzdetW1K1bV+rbkpGCWWiIiKgWM/iRU+vWrStMChZFEcnJyUhLS8OXX34pSTiSTmB5oUljoSEiotrH4EIzaNCgCu8rFAp4eHige/fuaNKkibG5SGLlj5yuZuSjsEQDOxsrmRMRERFJxyT70KSlpcHDw8PQW5MJeDgp4Wpvg6yCEsSl5SHE20XuSERERJKRbA6NKIr4/fffMXjwYPj4+Eh1W5KIIAh3JgbzsRMREdUyRheauLg4vPfee2jQoAFGjBgBBwcHbN68WYpsJLEgj9uFJiVH5iRERETSMuiRU2FhIbZu3YpvvvkGx48fR1hYGG7cuIHIyEiEhoZKnZEkEqzmCA0REdVOeo/QTJw4Ed7e3li+fDmeffZZJCUl4ZdffoEgCFAoJF8FThIqX+kUncJCQ0REtYveIzQrV67EzJkz8fbbb8PZ2dkUmchEyveiuXIzD6UaLaytWECJiKh20Psn2vr163HixAl4eXlh6NCh+PXXX1FaWmqKbCQxb1d72NtYoUQjIiEjX+449BAZecUoLNHIHYOIyCzoPUIzfPhwDB8+HFeuXMGaNWswadIk5OfnQ6vV4vz58wgJCTFFTpKAQiEgUOWIc0nZiE7JReDtScJUM4iiiMspudgdlYzd51NwNikLdjYKdG+kQt9QT/RoooKrPU9KJyK6F0EURdGYG4iiiD///BOrV6/Gzp07Ua9ePQwePBhLliyRKqMksrOz4erqiqysLLi4WO4eLK9vicSOf5Lw1uONMalHkNxxLJ5GK+L01Vu6EpNw8/4jZzZWAroE1kO/UE+Ehajh7qSsxqRERPKo6s9vgzfWKycIAvr27Yu+ffsiIyMD69evx5o1a4y9LZkID6mUX2GJBkdj0rE7KgV7L6TgZl6x7jVbawW6BtdDnxBP9GyqQkp2IXadS8auc8mITs1F+OU0hF9Owzs7zqK9X130DfXE48084V3HXsbfERGR/IweoTEXHKEps+tcMl79/hRC67vg19cekzuOxcjKL8H+SynYHZWC8MtpyC++MzfG1d4GvZqo0KeZGo8Fe8BRee9/Z8Sk5uLPqLJyczYpq8JrLX3roG8zT/QN9YR/PUeT/l6IiKpTVX9+s9BYmNi0XPT6LBz2NlaI+vBxKBTCwz+JDHI9swB7zqdg9/lkHI/LgEZ754+at6sd+jTzRJ8QNdr714WNnivOrt3Kx59RKfjzXDIiEjJw95/iJp7OePx2uWni6VzhEFkiInPDQvMfLDRlSjVaNP1gF0o0Ig7P6AHfug5yR6o1RFHEpZQc7IlK0U3qvVsTT2f0CVGjTzNPNPN2kaxopOUUYff5spGbv2JvovSu4uTn7oDHQz3Rt5knWvrUYYElIrPDQvMfLDR39Pk8HJdTcrHmxfbo0VgldxyzptGKOJVwZ1Lv1buWwwsC0L5hXfRppkZYiBoN3U3/KCgrvwR7L6RgV1QyDl1OQ1GpVveap4sdHm+mxuOhnujgV5f7EBGRWai2ScFkfoJUTrickovY1FwWGgMUlmhwJDodu88nY9+F1AdO6q1XzSuRXB1s8ExbHzzT1gd5RaUIv5yGP84lY/+FFCRnF2LdXwlY91cC6jraIqypGn1DPdElyB1Ka6tqzUlEJDWDC42fnx/Gjh2LMWPGoEGDBlJmIhMrP6SSRyBUXWZ+MfZfTNVN6i0o0X9Sb3VzVFqjf3Mv9G/uhcISDY7FpmPXuWTsOZ+CjLxibDmZiC0nE+GstEbPpir0beaJbo094GBbM/ITEenD4L+53njjDaxduxYfffQRevTogXHjxuHpp5+GUsm9MWq6IHXZkRU8pPLBkjILsOf2o6S/46Wd1Fvd7Gys0LOJGj2bqFGq0eJEfAZ2RSXjz6hkpGQX4efI6/g58jqU1gp0a+SBfs090bOJmhv5EZHZMHoOzb///ovVq1dj06ZNKC0txfDhwzF27Fi0adNGqoyS4ByaO85fz0b/JYfhYmeNf2f34SqY28on9e6OKluZdC4pu8LrpprUKyetVkTktUzdXjd3zwGyVgjoEnRnI7/qfnxGRATIMCm4pKQEX375JWbOnImSkhKEhoZi6tSpePHFF2vEX/wsNHcUlmjQ9INdEEXgxLu9oHK2kzuSbGrapF45iaKICzdysOvcDeyKSsblux5JKgSgnV9d9G3micdDPVHfgjby02pFZBWU4GZeMdwdbeHmaCt3JCKLUm2FpqSkBDt27MCaNWuwZ88edOrUCePGjcP169exbNky9OjRAxs3bjTmS0iChaaibp8eQMLNfGx8qSO6BNaTO061qsmTemuS2LQ7G/mdufafjfx8XHXLwQPM7EwwjVbErfxiZOQV42ZuMW7mFel+nZF3+/rtaxl5xbiVX6J73GhrrcCLj/hhYvcgPo4jqiYmLzSnT5/GmjVrsGnTJlhZWWHkyJEYP348mjRpovuYiIgIdO3aFQUFBYZ8CUmx0FQ0bm0E9l1MxUcDm2FUZz+545icOU7qrUmSMgvw57lk7IpKRsSVihv5NVI7oW+oF/o280RTr+rfyK9Eo8WtvGLc1JWRYmTklhWS9LxiZORWLCmZBSUw5G89R1sr5N3e4dnNwQZTegVjRMeGsLWu2fOniMydyQuNlZUVwsLCMG7cOAwaNAg2NpX/tZKXl4fJkyfXiLOdWGgqmv/HBXwdHodRnRvio4Ghcscxido0qbcmScspwp7zZXvdHItJr7CRX4O6DrrzpVr7GraRX1Gp5h4jJsXIuGskRVdecouQXVhq0O/DzcEGdR1t4e6oRF1HW9R1soW7o23Zr++6Xs+p7DGTtULAgUupmPf7Rd1ZaH7uDpjZtwn6hnrWiEfrRLWRyQtNQkICGjZsaHDA6sZCU9GPJxPx1tYz6BLojo0vdZI7jmTSc4uw+cRV7IqyjEm9cis/o+qPs8kI/89GfmoXZdkRDM080cDd4a7Rk3uUlLw7j4Byi/QvKAoBuiJSoaQ42sLdqfJ1NwcbgzcWLNVoseVkIj7fE4303CIAQNuGbnj3iaZo08DNoHuSfnKLSrHp76v491om6jkpoXJRQuVsB7WLEmoXO6id7eBib80/47WEyQtNQEAAIiIi4O7uXuF6ZmYm2rRpg7i4OENuazIsNBX9c/UWnv7yGDyclYh4t7fccYyWmJGPlYfi8MPJRN0PVUub1Cu3/OJShF9Kw66osrlJhhSTctYKAW6OZSMmZYVE+Z/Rk7vLihJ17G2q/ViH3KJSrAyPxcrDcSgsKfv/3BPNvTCjb2P+f81EsvJLsPbYFaw5Fo/M/JIHfqzSWlFWblyUUN0uOeWFR1VefFzs4MRHzDWeyQuNQqFAcnIyVKqKO82mpKSgQYMGKCoqMuS2JsNCU1FOYQmaz9kNAPj3gz5wdTDPCY7nr2fjq/BY/Hb2hu6RUksfVwzv2AC9mnKpsVyKSjU4FnOzbCO/CynILSytMGJSVkiU/xk9uTOKYk7/uk7OKsSiPZfw46lrEEXAxkrAyE5+eK1nEFdESeRmbhG+PRKP9X8l6Iqyfz1HDLm9I3ZqThFSsguRml2ElJzCh5aduznaWlUqOSrn8hJ055q9LXfTlovJCs3OnTsBAIMGDcK6devg6uqqe02j0WDfvn3Ys2cPLl26ZGD0+/vyyy/x6aef4saNG2jWrBkWL16Mxx57rEqfy0JTWad5+5CcXYhtEzqjbcO6csepMlEUcTwuA1+FxyL8cpruetdGHni1WwA6B7ibzQ9DS1D+V0xt/55cuJGN+X9cxKHb/590trPGaz2DMKqzH+xs+MPQEMlZhVh5KA4bTyToRsGaeDpjUo8g9G/uBav7jMoVlmiQdrvkpGTf/t+c24Unu1BXfnL0GEV0trPWjfione3KCo+uBJU98lK5KHmMiAmYrNAoFGXPnQVBwH8/1cbGBn5+fvjss88wYMAAA2Lf35YtWzBy5Eh8+eWXeOSRR/D111/jm2++wfnz56t09AILTWUvfPM3jsSkY8EzzTG0fc0/vkKrFbH7fAq+Co9FZGImgLK5E0+08MYrXQMQWt/1wTcgqgaHLqdh3u8XcDE5BwBQv449ZvRtjCdbePO08ypKzMjHivBYbD15DcWasiLT0scVk3sGo1cTlWT/He8e3dGN8GQXIkU34lNWiO5eFfkwbg42t0d87ozw/HeOj4ezkgsR9GDyR07+/v6IiIhAvXrVs4dJx44d0aZNG6xYsUJ3rWnTphg0aBDmz59f6eOLiooqPPbKzs6Gr68vC81d5uyMwtpjV/DSY/5494kQuePcV3GpFj/9k4SvDsUiLi0PQNnz8Wfb+eClxwI4X4FqHI1WxPbT1/B/uy8hJbvs76EWPq54p39TdApwf8hnW66Y1Fx8eTAGP0de1z1C7uBfF5N7BOGx4HqyjPKJoojcolKkZBeVFZycO6M+qXeN/qRkF6H4rknxDyIIgLujra7kBHg4YWL3QLjzEfk9mfy07fj4eEM/VW/FxcU4deoU3n777QrX+/Tpg2PHjt3zc+bPn48PP/ywOuKZrSDV7UMqU2vmmU7lKxm+PRKP5OxCAICLnTVGdm6IMV384eHMP/xUM1kpBDzbzhcDWnjjm8Nx+Co8FmeuZeH5lcfRu6kab/drovvzR2Vz4ZYfiMHv527o9gh6LLgeJvcIQkeZC6AgCHC2s4Gznc0Dv2eiWLajdMrdj7XuGv0pL0SpOUUo1YpIzy1Gem4xzt8ADlxKw74LKVg3tgP/gWYEvQrNkiVL8PLLL8POzg5Llix54MdOmTLFqGB3S09Ph0ajgVqtrnBdrVYjOTn5np8za9YsTJ8+Xfd++QgN3RF8+w9nTA0rNOm5RVh79ArW/3VFt8eI2kWJ8Y8GYFjHBlyVQGbD3tYKr/UKxvMdGmDx3svYHJGIvRdScOBSKoZ3aICpvYMteuL6P1dvYfmBGOy9kKq7FhaixuQeQWjpW0e+YAYQBAF1HGxRx8EWjT2d7/txWq2IjPxi3QTm5KxCLNsfgys38/HMimNYPaY9WvjUqb7gtYhej5z8/f1x8uRJuLu7w9/f//43FQRJl21fv34d9evXx7Fjx9C5c2fd9blz5+K7777DxYsXH3oPzqGp7GZuEdp+vBeCAER9+DgcbOUtCldv5mPl4Vj8ePKabul1gIcjXu0aiIGtvTnZjsxeTGoOPvnjou4HuJPSGhO6B2LsI/4Ws4pGFEX8HZ+BZftjcCQmHUDZI5gBLbwxqUcgmnha3t/PqdmFGLMmAudvZMPB1grLR7RBj8aqh3+ihTDJI6e7HzNV5yOnevXqwcrKqtJoTGpqaqVRG6o6dycl3BxscCu/BHFpebJNqo26noWvwuPw25nrKN90tqVvHUzoFog+IWpOpKRaI0jljG9Gt8dfsTcx7/cLOJuUhU//vITv/krAm483xtOt69935Y65E0UR4ZfTsGx/DE4m3AJQtt/QoNb1MaF7IALN7EwwKalc7LDllU6YuOE0DkenY/y6k5j/dHM8155PFfRhFmP3tra2aNu2Lfbs2YOnn35ad33Pnj0YOHCgjMnMX7DKGSeuZCAmNbdaC40oivgr7ia+Co/TLXMFgG6NPPBqt0B0Cqhb65f5kuXqHOiOnyc9gl/OXMfCXZeQlFmAN3/8F98eice7/Zvi0eDac2Bs+erE5QdicDap7JBTWysFnmvvg1e6BsK3roPMCWsGZzsbfDu6Pd7efgbbTydhxrYzuJFViCm9gvh3YRUZXGiGDBmCdu3aVZqo++mnn+LEiRP48ccfjQ53t+nTp2PkyJFo164dOnfujJUrV+Lq1at49dVXJf06liZQ5YQTVzIQnZpTLV+v7C+3ZKwIj8O/dy29HtDCG690C0Azby69JsugUAgY2Ko+Hm/miXXHrmDZgRhcuJGNF779G90aeeCd/k0fOBejpivVaPHb2RtYfiAGl1PK5unZ21hhRMcGeKlrANQudjInrHlsrRX47NmW8HK1w/IDsfh872UkZxfgfwNDDT6qw5IYXGjCw8Mxe/bsStf79u2L//u//zMq1L0MHToUN2/exEcffYQbN24gNDQUv//+u1mdJ1UTBVXTxOCiUg1++icJX4fHIS79ztLr59r54qXHAtDAnf9KI8tkZ2OFV7oF4tl2vli6Pxrf/ZWA8MtpOBydhmfb+mJ6n0Zm9cO/uFSLHf9cw4qDsbhyMx8A4Ky0xqguDTH2EX8uTX4IQRDw1uNN4Olqj9k/n8OmE4lIyS7CsuGtZZ/nWNMZvA+Nvb09IiMj0bhx4wrXL168iNatW6OgoECSgFLhpOB7O3Q5DaNWn0CghyP2vdFd8vvnFJZg04mypdfl+3G42FljdBc/jO7iZ9ErPIju5Up6Hhb+eRG/ny2bM2hvY4WXuwbg5a4BcKzBK/wKSzT44WQivjoYi+tZZdssuDnYYOwj/hjVxQ+u9uZ5vIqc/oxKxpRN/6CoVIuWPq74dkx7i/w70+Qb67Vv3x5PPvkkPvjggwrX58yZg19++QWnTp0y5LYmw0Jzb9czC9Dlk/2wUgi48FFf2FpLM6yZllOENUfj8d3xBOTcXnrt6WKH8Y/54/kOXHpN9DCnEjIw97cLOH01EwDg4azE9LBGeLatT416/JBXVIoNfydg5aF43enjHs5KvPxYAIZ3bFCjS5g5OJWQgXHrTiIzvwQN3R2w7sUO8KtnWXvVmLzQ7Ny5E8888wyGDx+Onj17AgD27duHTZs24ccff8SgQYMMCm4qLDT3JooiQmf/ibxiDfa83hXBauOe2SfczMPKQ3H48dQ13a6ZgR6OeKVbIAa1qi9ZYSKyBKIo4o9zyViw6yISbj++CVY54Z3+TdG9sYesk0WzCkqw7tgVrD565+Tr+nXs8Wq3ADzbzpfnV0koNi0Xo1efwLVbBXB3tMXqMe3Nbp8eY5i80ADAb7/9hnnz5iEyMhL29vZo0aIFZs+ejW7duhl6S5Nhobm/gcuO4N9rWVgxog36Nfcy6B7nkrLwVXgsfj97Q7f0unWDOni1WyDCmnLpNZExiku1+P54Apbsj9aVhy6B7ninf9Nq326h/OTr7/5K0B3u6OfugIndgzCoNf/RYiqpOYV4cU0Eoq5nw97GCstHtEbPJpaxbUm1FBpzwkJzf9N/iMT200mYHtYIU3oFV/nzRFHEX7E3sSI8Foej03XXuzf2wIRugejgz6XXRFLKKijBlwdisOboFRRrtBAE4OlW9fHm443hXcfepF87Jfv2ydd/X9Ud1thY7YyJPQLxRHOvGvUYrLbKLSrFhO9P4XB0OqwUAuYOCsXzHWr+wcLGMvlZTuVOnTqFCxcuQBAEhISEoHXr1sbekqpZsKrsMVNVVzpptCJ2RyVjxe3zaYCys2uebOGFV7oFoqkXCyORKbja22BW/6Z4oVND/N/uS/g58jq2/5OE387ewNhH/TGheyBc7KSdfJuYkY+vwst28C4/+bp5fVdM7hnE0ddq5qS0xuox7fH2trPYdvoa3t5+FjeyCjGtdzD/8QgjRmhSU1Px/PPP4+DBg6hTp07ZwVxZWejRowc2b94MDw8PqbMahSM097fnfApeWn8STb1c8MfUx+77cUWlGmw/nYSVh+IQf3vptZ2NAkPb+WL8YwHcIIuomp25lom5v13A3/EZAIC6jraY1jsYwzo0gI2RIyaxabn48kAsfopM0p183d7PDZN7BqOrTCdfUxlRFLFoz2Us3R8DABjazhdzn669e9WY/JHT0KFDERsbi++++w5NmzYFAJw/fx6jR49GUFAQNm3aZFhyE2Ghub/49Dz0+L+DUForcP6jvpW2Xs8pLMGGv69i9ZF4pOaUrWJwtbfB6M4NMbqLH/eVIJKRKIrYeyEV8/+4gLi0sn9oBNRzxIy+TfB4M7XexePCjbKTr387W/NOvqaKNvydgPd/OgetCPRo7IHlI9rUyr1qTF5oXF1dsXfvXrRv377C9RMnTqBPnz7IzMw05LYmw0Jzf6UaLUJm/4niUi0OvdVDt8ldak4h1hy9gu/vWnrt5WqHcY/6Y1gHLsckqklKNFpsjkjE4j2XcTOvGEDZiMo7/ZuidQO3h35+ZGImlu2Pwd4LKbprvZuqMblnEFpZ0Ioac7PnfApe23QahSVatPBxxepauFeNyefQaLVa2NhUflZrY2MDrVZr6G1JBtZWCgTUc8TF5BxEp+ZAK4pYeTgOW+9aeh2kcsKr3QLxVEtvrmIgqoFsrBQY2akhBrXyxtfhcfjmSBwirtzC018ew4AWXpjxeJN77sj9d9xNLDsQo5vYLwhA/+ZemNQ9CCHe/MdfTRcWosbGlzph3NoInLmWhWdWHLPIvWoAI0ZoBg4ciMzMTGzatAne3t4AgKSkJIwYMQJubm7YsWOHpEGNxRGaB5u88TR+PXMDfu4OuJqRr1t63aZBHUzoHoReTVSc/EdkRm5kFeCz3Zex7fQ1iCJgYyVgdGc/TO4ZBFd7GxyKTsey/dGIuFJ28rWVQsCgVmUnX5cfiULmIy4tF6PXnEBiRgHq3t6rpraMrJn8kVNiYiIGDhyIc+fOwdfXF4Ig4OrVq2jevDl+/vln+Pj4GBzeFFhoHmzx3stYvDda937PJiq82i0Q7f3cOPmPyIydv56Neb9fwJGYshEYV3sb1K9jj/M3sgGUnXz9bDsfvNqNJ1+bu7ScIoxdG4GzSVmwt7HCsuGt0aup+e9VU2370OzZswcXL16EKIoICQlB7969jbmdybDQPNiV9DxM2HAaTTyd8Uq3ADTx5H8jotpCFEWEX07D/N8v4lJKDoCyFYrDOzTEy10D4OlqPodf0oPlFZVi4obTCL+cBoUAzH26OYaZ+V413FjvP1hoiMjSabQifvonCcnZhRja3rfWTR6lMiUaLd7ZfhY/nroGAJjSMwivhzUy29F2k0wKXrJkSZU/dsqUKfrcmoiITMxKIeCZtjVrOgBJz8ZKgYVDWsDL1Q5L9sdgyf4Y3MgqxLzBzY3en6gm02uExt/fv2o3FQTExcUZHMoUOEJDRESWZuPfV/HeT2ehFcuOpVk+vI3ZbbnBR07/wUJDRESWaO/5FEy+vVdN8/ple9V4OJvP48aq/vw2euypuLgYly5dQmlpqbG3IiIiIon1DlFj00udUNfRFmeTsjB4xVHEpVXt7D5zYnChyc/Px7hx4+Dg4IBmzZrh6tWrAMrmznzyySeSBSQiIiLjtG7ghm0TuqBBXQckZhTgmRXHcPrqLbljScrgQjNr1iz8+++/OHjwIOzs7iz56927N7Zs2SJJOCIiIpKGfz1HbJvQBS18XHErvwTDVx3H3vMpD/9EM2Fwofnpp5+wbNkyPProoxWWgoWEhCA2NlaScERERCQdD2clNr3UCd0be6CwRIuXvzuJDX8nyB1LEgYXmrS0NKhUqkrX8/LyzHatOxERUW3nqLTGqlHt8Fw7H2hF4N0d5/DZ7ksw9zVCBhea9u3b47ffftO9X15iVq1ahc6dOxufjIiIiEzCxkqBBc+0wNRewQCApftj8NbWMyjRmO/h0gYvRp8/fz769u2L8+fPo7S0FF988QWioqLw119/ITw8XMqMREREJDFBEPB6WCN4udrh3Z/OYeupa0jNKcKKEea3Vw1gwAhNZGQkAKBLly44evQo8vPzERgYiN27d0OtVuOvv/5C27Ztpc5JREREJvB8hwZYNaot7G2scOhyGoau/AupOYVyx9Kb3hvrKRQKtG7dGuPHj8fw4cPh6upqqmyS4sZ6RERE9xeZmIlxayNwM68YPm72WDe2AwI9nOSOZbqN9Y4ePYo2bdrg7bffhpeXF0aOHIkDBw4YFZaIiIjk1cq3DrZN6IKG7g64dqsAQ1Ycw6kE89mrRu9C07lzZ6xatQrJyclYsWIFEhMT0bt3bwQGBmLu3Lm4du2aKXISERGRifnd3qum5V171eyOSpY7VpVIcpZTbGws1qxZg/Xr1+PGjRsICwvD77//LkU+yfCRExERUdXkF5di0obTOHApDQoB+HBgKEZ2aihLlmo/nDI3NxcbNmzAO++8g8zMTGg0GiluKxkWGiIioqor1Wjx3k/nsDkiEQAwqUcg3uzTuNr3mqu2wynDw8MxevRoeHp6YsaMGRg8eDCOHj1q7G2JiIhIRtZWCswf3Byv924EAFh+IBZv/Phvjd2rxqCF5omJiVi7di3Wrl2L+Ph4dOnSBUuXLsVzzz0HR0dHqTMSERGRDARBwNTewfB0VeKdHeew/XQS0nKKsOKFtnCqYXvV6J0mLCwMBw4cgIeHB0aNGoWxY8eicePGpshGRERENcDQ9g2gcrbDxA2ncTg6HUO//gtrxrSHysXu4Z9cTfR+5GRvb49t27bh2rVrWLBgAcsMERGRBejRRIXNL3eCu6Mtoq5n4+kvjyEmNVfuWDqSTQqu6TgpmIiIyHgJN/MwevUJXLmZjzoONvh2dDu0bVjXZF+v2iYFExERkeVo6H57rxrfOsjML8HwVX/jzxqwVw0LDREREenF3UmJTS91RK8mKhSVajHh+1P47q8rsmZioSEiIiK9Odha4+uRbTGsQwNoReD9n6OwaM9l2fKw0BAREZFBrK0UmPd0KN4IawSFADT1dJYvi2xfmYiIiMyeIAh4rVcw+jX3RJBKvkLDERoiIiIympxlBmChISIiolqAhYaIiIjMHgsNERERmT0WGiIiIjJ7LDRERERk9lhoiIiIyOyx0BAREZHZY6EhIiIis8dCQ0RERGaPhYaIiIjMHgsNERERmT0WGiIiIjJ7LDRERERk9qzlDlBdRFEEAGRnZ8uchIiIiKqq/Od2+c/x+7GYQpOTkwMA8PX1lTkJERER6SsnJweurq73fV0QH1Z5agmtVovr16/D2dkZgiDIHadGys7Ohq+vLxITE+Hi4iJ3HIvH70fNwu9HzcLvR81iyu+HKIrIycmBt7c3FIr7z5SxmBEahUIBHx8fuWOYBRcXF/4FUYPw+1Gz8PtRs/D7UbOY6vvxoJGZcpwUTERERGaPhYaIiIjMHgsN6SiVSsyePRtKpVLuKAR+P2oafj9qFn4/apaa8P2wmEnBREREVHtxhIaIiIjMHgsNERERmT0WGiIiIjJ7LDRERERk9lhoLNz8+fPRvn17ODs7Q6VSYdCgQbh06ZLcsei2+fPnQxAETJs2Te4oFi0pKQkvvPAC3N3d4eDggFatWuHUqVNyx7JIpaWleO+99+Dv7w97e3sEBATgo48+glarlTuaRTh06BCefPJJeHt7QxAE/PTTTxVeF0URc+bMgbe3N+zt7dG9e3dERUVVSzYWGgsXHh6OSZMm4fjx49izZw9KS0vRp08f5OXlyR3N4kVERGDlypVo0aKF3FEs2q1bt/DII4/AxsYGf/zxB86fP4/PPvsMderUkTuaRVqwYAG++uorLFu2DBcuXMDChQvx6aefYunSpXJHswh5eXlo2bIlli1bds/XFy5ciEWLFmHZsmWIiIiAp6cnwsLCdOcpmhKXbVMFaWlpUKlUCA8PR9euXeWOY7Fyc3PRpk0bfPnll/j444/RqlUrLF68WO5YFuntt9/G0aNHcfjwYbmjEIABAwZArVbj22+/1V175pln4ODggO+++07GZJZHEATs2LEDgwYNAlA2OuPt7Y1p06Zh5syZAICioiKo1WosWLAAr7zyiknzcISGKsjKygIA1K1bV+Yklm3SpEl44okn0Lt3b7mjWLydO3eiXbt2ePbZZ6FSqdC6dWusWrVK7lgW69FHH8W+fftw+fJlAMC///6LI0eOoH///jIno/j4eCQnJ6NPnz66a0qlEt26dcOxY8dM/vUt5nBKejhRFDF9+nQ8+uijCA0NlTuOxdq8eTNOnz6NiIgIuaMQgLi4OKxYsQLTp0/HO++8gxMnTmDKlClQKpUYNWqU3PEszsyZM5GVlYUmTZrAysoKGo0Gc+fOxbBhw+SOZvGSk5MBAGq1usJ1tVqNhIQEk399FhrSmTx5Ms6cOYMjR47IHcViJSYmYurUqdi9ezfs7OzkjkMAtFot2rVrh3nz5gEAWrdujaioKKxYsYKFRgZbtmzB999/j40bN6JZs2aIjIzEtGnT4O3tjdGjR8sdj1D2KOpuoihWumYKLDQEAHjttdewc+dOHDp0CD4+PnLHsVinTp1Camoq2rZtq7um0Whw6NAhLFu2DEVFRbCyspIxoeXx8vJCSEhIhWtNmzbFtm3bZEpk2d566y28/fbbeP755wEAzZs3R0JCAubPn89CIzNPT08AZSM1Xl5euuupqamVRm1MgXNoLJwoipg8eTK2b9+O/fv3w9/fX+5IFq1Xr144e/YsIiMjdW/t2rXDiBEjEBkZyTIjg0ceeaTSVgaXL19Gw4YNZUpk2fLz86FQVPzRZWVlxWXbNYC/vz88PT2xZ88e3bXi4mKEh4ejS5cuJv/6HKGxcJMmTcLGjRvx888/w9nZWfcM1NXVFfb29jKnszzOzs6V5i85OjrC3d2d85pk8vrrr6NLly6YN28ennvuOZw4cQIrV67EypUr5Y5mkZ588knMnTsXDRo0QLNmzfDPP/9g0aJFGDt2rNzRLEJubi5iYmJ078fHxyMyMhJ169ZFgwYNMG3aNMybNw/BwcEIDg7GvHnz4ODggOHDh5s+nEgWDcA939asWSN3NLqtW7du4tSpU+WOYdF++eUXMTQ0VFQqlWKTJk3ElStXyh3JYmVnZ4tTp04VGzRoINrZ2YkBAQHiu+++KxYVFckdzSIcOHDgnj8zRo8eLYqiKGq1WnH27Nmip6enqFQqxa5du4pnz56tlmzch4aIiIjMHufQEBERkdljoSEiIiKzx0JDREREZo+FhoiIiMweCw0RERGZPRYaIiIiMnssNERERGT2WGiIiIjI7LHQEFG1mDNnDlq1aiXb13///ffx8ssvS3a/IUOGYNGiRZLdj4iMw52CichogiA88PXRo0frTgt3d3evplR3pKSkIDg4GGfOnIGfn58k9zxz5gx69OiB+Ph4uLi4SHJPIjIcD6ckIqPduHFD9+stW7bggw8+qHBCtb29PZycnODk5CRHPHz77bfo3LmzZGUGAFq0aAE/Pz9s2LABEyZMkOy+RGQYPnIiIqN5enrq3lxdXSEIQqVr/33kNGbMGAwaNAjz5s2DWq1GnTp18OGHH6K0tBRvvfUW6tatCx8fH6xevbrC10pKSsLQoUPh5uYGd3d3DBw4EFeuXHlgvs2bN+Opp56qcK179+6YPHkyJk+ejDp16sDd3R3vvfce7h60/vLLLxEcHAw7Ozuo1WoMGTKkwj2eeuopbNq0ybD/aEQkKRYaIpLN/v37cf36dRw6dAiLFi3CnDlzMGDAALi5ueHvv//Gq6++ildffRWJiYkAgPz8fPTo0QNOTk44dOgQjhw5AicnJ/Tt2xfFxcX3/Bq3bt3CuXPn0K5du0qvrVu3DtbW1vj777+xZMkSfP755/jmm28AACdPnsSUKVPw0Ucf4dKlS9i1axe6du1a4fM7dOiAEydOoKioSOL/MkSkLxYaIpJN3bp1sWTJEjRu3Bhjx45F48aNkZ+fj3feeQfBwcGYNWsWbG1tcfToUQBlIy0KhQLffPMNmjdvjqZNm2LNmjW4evUqDh48eM+vkZCQAFEU4e3tXek1X19ffP7552jcuDFGjBiB1157DZ9//jkA4OrVq3B0dMSAAQPQsGFDtG7dGlOmTKnw+fXr10dRURGSk5Ol/Q9DRHpjoSEi2TRr1gwKxZ2/htRqNZo3b65738rKCu7u7khNTQUAnDp1CjExMXB2dtbNyalbty4KCwsRGxt7z69RUFAAALCzs6v0WqdOnSpMaO7cuTOio6Oh0WgQFhaGhg0bIiAgACNHjsSGDRuQn59f4fPt7e0BoNJ1Iqp+nBRMRLKxsbGp8L4gCPe8ptVqAQBarRZt27bFhg0bKt3Lw8Pjnl+jXr16AMoePd3vY+7F2dkZp0+fxsGDB7F792588MEHmDNnDiIiIlCnTh0AQEZGxgO/NhFVH47QEJHZaNOmDaKjo6FSqRAUFFThzdXV9Z6fExgYCBcXF5w/f77Sa8ePH6/0fnBwMKysrAAA1tbW6N27NxYuXIgzZ87gypUr2L9/v+7jz507Bx8fH11pIiL5sNAQkdkYMWIE6tWrh4EDB+Lw4cOIj49HeHg4pk6dimvXrt3zcxQKBXr37o0jR45Uei0xMRHTp0/HpUuXsGnTJixduhRTp04FAPz6669YsmQJIiMjkZCQgPXr10Or1aJx48a6zz98+DD69Oljmt8sEemFhYaIzIaDgwMOHTqEBg0aYPDgwWjatCnGjh2LgoKCB25u9/LLL2Pz5s26R1flRo0ahYKCAnTo0AGTJk3Ca6+9pttNuE6dOti+fTt69uyJpk2b4quvvsKmTZvQrFkzAEBhYSF27NiBl156yXS/YSKqMu4UTES1niiK6NSpE6ZNm4Zhw4YBKNuHplWrVli8eLFB91y+fDl+/vln7N69W8KkRGQojtAQUa0nCAJWrlyJ0tJSye5pY2ODpUuXSnY/IjIOVzkRkUVo2bIlWrZsKdn9pDzokoiMx0dOREREZPb4yImIiIjMHgsNERERmT0WGiIiIjJ7LDRERERk9lhoiIiIyOyx0BAREZHZY6EhIiIis8dCQ0RERGbv/wGpnzKrOhwxhgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# plot vacf for test data with velocities (PRM_NCBOX, TRJ_NCBOX)\n", + "v_vacf_plot = v.plot_vacf()" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "3e9151ef-54bb-44ef-9f53-76c90baf6e34", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGwCAYAAACKOz5MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQwklEQVR4nO3deVhUZfsH8O8wwLAPO4gssuMuiplL7tpiLvmWZVYW+ZZpbuRb2eKaUb25ZJa5ZZap9assq7dcU0NzQTDTZBMURNlhhnWAmfP7AxklUJlh4Awz3891zRVzzpkzt5LMzfPcz/1IBEEQQERERGTELMQOgIiIiOhOmLAQERGR0WPCQkREREaPCQsREREZPSYsREREZPSYsBAREZHRY8JCRERERs9S7AAMQaPR4OrVq3B0dIREIhE7HCIiImoGQRBQWloKHx8fWFjcfgzFJBKWq1evws/PT+wwiIiISA9ZWVnw9fW97TUmkbA4OjoCqPsDOzk5iRwNERERNYdSqYSfn5/2c/x2TCJhqZ8GcnJyYsJCRETUzjSnnINFt0RERGT0jCJhqa2txRtvvIHAwEDY2toiKCgIS5cuhUajETs0IiIiMgJGMSX07rvv4pNPPsHWrVvRtWtXxMfH45lnnoFcLsecOXPEDo+IiIhEZhQJyx9//IHx48djzJgxAIBOnTphx44diI+PFzkyIiIiMgZGMSU0aNAgHDhwACkpKQCAP//8E3FxcXjggQeavF6lUkGpVDZ4EBERkekyihGWV155BQqFAhEREZBKpVCr1Vi+fDkmT57c5PWxsbFYsmRJG0dJREREYjGKEZavvvoK27Ztw/bt25GQkICtW7fi/fffx9atW5u8fsGCBVAoFNpHVlZWG0dMREREbUkiCIIgdhB+fn549dVXMXPmTO2xt956C9u2bUNSUtIdX69UKiGXy6FQKNiHhYiIqJ3Q5fPbKEZYKioqGu0hIJVKuayZiIiIABhJDcvYsWOxfPly+Pv7o2vXrkhMTMTKlSsRHR0tdmhERERkBIxiSqi0tBRvvvkmdu3ahby8PPj4+GDy5MlYuHAhrK2t7/h6TgkRERG1P7p8fhtFwtJSTFiIiIjan3ZXw0JEJKbqWg1Kq2rEDoOIbsMoaliIiNqSIAhIyinF0bQC/J5agJMZRaisUSPI3R6R/i7oHeCM3v4uCPNyhNTizrvIElHrY8JCRGYhR1GF31PzcTStAHFphSgoUzW6Jr2gHOkF5fg24QoAwEFmiZ5+cvT2d0FvfxdE+jvD2e7OdXVEZHhMWIjIJJVW1eBEehHi0goQl1aAtLyyBudtraS4K9AV94S6Y1CoOzwdbfBnVgkSMouRkFmMM5klKFPV4mhaIY6mFWpfF+Rhr01gegc4I9STozBEbYFFt0RkEmrUGvyZVVKXoKQW4ExWCWo1N368WUiA7r7OGBTihkEhHugd4AyZpfSW91NrBKTkltYlMJdLkJhZjPSC8kbXcRSGSH9cJUREJk8QBFzML0dcaj7i0gpxPL0QZaraBtcEuNlhUIg77gl1R/8gd8jtrFr0nsXl1UjMqktgEjKL8WdWCcqr1Y2u4ygMUfMwYSEik5RfqsKxi3WFskfTCnBNUdXgvIudFQaEuGPQ9Yefq12rxqPWCEjOKdVOIyVmliCjiVEYR5klevo5o7e/MyIDXNDbz6XFyRORKWDCQkQmobJajRMZhdrVPEk5pQ3OW1taoG8nFwwK8cCgEHd09XGChcgjGUXl1Ui8nsAkXC7Bn1dKUNHEKExw/ShMQN1ITKing+ixE7U1JixE1C6pNQLOZSsQl1aA31PzkXC5BNXqhnuKdfVxqhtBCXVH306usLG6dR2KMahVa5CcW4qEzBIkXq5LZC4VVjS6zlFmiV7+znXLqq//V27LURgybUxYiKjduFxYri2UPXaxEIrKhg3cOjrbYlCIOwaGumNgsBvcHGQiRWo4hWUqnKlfkXSbUZgQTwf09nfWjsSEeHAUhkwLExYiMlrF5dU4drEQcWn5iEsrQFZRZYPzjjJL9A92wz2h7hgY4o5Ad3tIJKb9Id3sURgbS/Tyu5HA9PJz5igMtWtMWIjIaFTVqJFwuRi/Xx9FOXdVgZt/6lhaSNA7wEU7zdOjoxyWUu4aUlimQmLmjb4wf2YpUFnTcBRGIgFCPBy0q5F6+7sgmKMw1I4wYSEi0Wg0Ai7kKBGXWtew7dSlIlTVNKxDCfNyqCuUDXVDv0A32MvYw/JOatUaJOWUXi/orUtkLjcxCuNsZ4UhYR4YHuGJIWEe7AlDRo0JCxG1qasllYhLLcDvaQU4llaAwvLqBuc9HWXaEZRBIe7wdLIRKVLTUnDzKMzlYpy90nAUxkICRAW4YliEJ0Z09kSop4PJT69R+8KEhYhaXXF5NTbFpeOXv3IadYC1s5bi7iA3DLzetI0flG2jVq3BmawSHEzKw8GkvEbLwH1dbDE8whPDIzxxd5Cb0a+wItPHhIWIWo2yqgabf8/A5rgMbWdZCwnQ088Z94TUFcpG+rvA2pJ1KGK7UlyB35LzcfBCLo5eLER17Y2pOVsrKQaGuGNEZ08MC/eEt5yjXtT2mLAQkcFVVNdi67HLWH/kIkoq6pYed+7ghBlDgzE4zIOrVYxcRXUtjqUV4mByHg5eyEOOsmGX4K4+ThgR4Ynhnb3Qo6OchbvUJpiwEJHBVNWoseNkJj767SIKylQA6rq0xowKx/3dvPnB1g4JgoC/rynxW1IeDiTl4UxWSYOVW+4O1hga7okREZ4YFOoORxsmo9Q6mLAQUYvVqDX4v/gr+PBgqnbPHn9XO8wdGYrxvTpyMz8TUlimwqHkfBxMysORlHyU3rSJpJVUgrsCXTEs3BMjOnsh0N1exEjJ1DBhISK9qTUCfjiTjdX7U5FZVLdstoPcBrOGh+KRKF9YsUeKSatRa3DqUhEOXsjDweQ8pOc3LKgOcrevW3UU4YmoTq6sVaIWYcJCRDrTaAT8ej4HK/elIC2vDEDd1MCMoSF4vJ8/V5SYqYyCchxMysNvSXk4kVGIGvWNjwxHmSXuCXPH8AgvDA33gLsJbJtAbYsJCxE1myAI+C05Dyv2puD8VSUAQG5rhelDgjF1QADsrNnUjeqUVtUgLrWgLoFJzkNB2Y1+OxIJ0NPXGSMiPDEswhNdfZzMdil7rVqDgrJq5CqroBYE9PR15hTqLTBhIaJmOZZWgPf3JiMhswQA4CCzxLODAvHsPYFwYqEl3YZGI+BstuJ6z5dcnMtWNjjv7WSDYdd7vgwMcTOJxFcQBBSVVyNXqUJuaRXylFXIUdz4OlepQo6yCgVlqgZFzCM7e2Lt4705StkEJixEdFunLxfh/T0p+CO9EABgY2WBqQM6YfrgYLjYs5U76S5XWaVddRSXWtCg4661pQUGBLtpR198XexEjLQxQRBQpqpFrlJVl4RcTz5ylVXIK61CjqLueX6pCtVqzZ1vCEBqIYGnowyF5dWortWgX6ArNk2N4oqrf2DCQkRNOpetwIq9yfgtOR8AYC21wOP9/DFjWDA8Hdk4jAyjqkaNExlFOHghFweS8nCluOGO3OFejhjeuW70JdLPuVU3u6yqUSPv+ohI7k2JSP0j7/qoSEW1+s43u87N3hpeTjbwcpLBy8kGnk428G7wXAY3exmkFhKcSC/EtK3xKFXVoltHJ2x95i64sdZHiwkLETWQkluKVftS8Mu5HAB1v/1NivLFi8ND0dHZVuToyJQJgoC0vDIcSKprWHc6sxhqzY2PHX03a7y5TiRHeWNK5sbzuiSlvslhczjaWMLrevLheT358HK8/l+5DbycbODhINN5ZdS5bAWmfnoSheXVCPKwx7Zn+8GH/+4AMGEROxwio3GpoByr96fghz+vQhDqCiMn9OqIOSNC0Yn9NEgEJRXVOJxS1/PlUHI+FJU3EgqphQR9/F0wvLMn+nZygbKqttl1Ircjs7SAt9wGXo51iUjdaMiNpKQ+QWnNOpv0/DI8sekEriqq4CO3wRfT+iHYw6HV3q+9YMJCZOaySyrx4YFU/N/pK9rfZu/v5o15o8IQ5uUocnREdWrVGiRmleDAhbpl08m5pXd+0U0sr9eJeN40HXPjceO5k42lUaxYulpSiSc2n0B6fjnc7K2xNfoudOsoFzssUTFhITJTecoqfPRbGnaczNIWBw4L98BLo8PN/gcjGb+sogr8lnx9p+lrpXBzuFEr4ulYl3x4y2987WZv3e62higsU2HqlpM4l62Eo8wSm6ZGoV+Qm9hhiYYJC5GZKS6vxieHL2LrH5dQVVOXqPQPcsP8e8PQJ8BV5OiI6GalVTWYtjUeJzKKILO0wLonemN4hJfYYYmCCQuRmVBW1WDT7xn4NC4DZdf3f+nt74z5o8MxIMRd5OiI6FaqatR4cXsC9l/Ig6WFBCsm9cT4Xh3FDqvN6fL53f47+RCZoYrqWnx27BLWH07XFi129XHC/NHhGBruYRTz9UR0azZWUqx7og9e/uYsdiVmY+5XZ6CorMFT/TuJHZrRMoqEpVOnTrh8+XKj4zNmzMBHH30kQkRExqmqRo3tJzLx8aE0bVv0UE8HxIwKw71dvdvdfD6RObOSWmDFIz0ht7XCZ8cuYeEP56GoqMGLw0P4S0cTjCJhOXXqFNTqG017zp07h1GjRuGRRx4RMSoi41Gj1uDr+Cx8eCANOcoqAECAmx3mjgzFuJ4duU8JUTtlYSHBorFdILe1wgcHUrFiXwpKKmvw+gOd+QvIPxhFwuLh4dHg+TvvvIPg4GAMGTJEpIiIjINaI+D7xGx8cCAVmUUVAIAOchvMHhGKh/v4wqoVO4QSUduQSCSYNyoMclsrLP3pb2yOy4CisgbvTOzeql2A2xu9E5ZLly7h999/x6VLl1BRUQEPDw9ERkaif//+sLHRv8V3dXU1tm3bhpiYmFsOialUKqhUKu1zpVLZ5HVE7ZVGI+CXczlYuS8ZF/PLAQDuDjK8OCwYj93lz03UiExQ9KBAyG2t8PK3Z/HN6StQVtZgzeRI/nu/TueEZfv27VizZg1OnjwJT09PdOzYEba2tigqKsLFixdhY2ODKVOm4JVXXkFAQIDOAX3//fcoKSnB008/fctrYmNjsWTJEp3vTWTsBEHAwaQ8rNibgr+v1SXiznZWmD4kGE/1DzCJHW+J6Nb+1ccXTrZWmLk9AXv/zkX0Z6ew4akoOMj4b1+nZc29e/eGhYUFnn76aYwbNw7+/v4NzqtUKvzxxx/YuXMnvv32W3z88cc616Hce++9sLa2xo8//njLa5oaYfHz8+OyZmq3BEHAsYuFeH9vMhIzSwAADjJLTLsnENGDAuHEHV6JzMqxiwX499Z4lFer0dNXjs+eucskd1JvtT4sP//8M8aMGdOsawsKCpCRkYG+ffs29/a4fPkygoKC8N1332H8+PHNfh37sFB7Fn+pCO/vTcbx9CIAgI2VBZ4eEIjnBweZ5A8oImqes1dKMPXTkyiuqEGopwO+eLYfvOWmtat6u20ct3jxYqxfvx5ZWVmwtGz+8BcTFmqP/rqiwIp9yTiUnA8AsJZa4PF+/pgxLBiejqb1Q4mI9JOWV4onNp1EjrIKHZ1tsW1aPwSa0Malunx+611+nJCQgL/++kv7/IcffsCECRPw2muvobq6Wuf7aTQabNmyBVOnTtUpWSFqb1JzSzH9i9MYuzYOh5LzYWkhweS7/HHoP0OxeFxXJitEpBXi6YhvXuiPQHd7ZJdU4pFP/sDfV81zoYneCcvzzz+PlJQUAEB6ejoee+wx2NnZ4f/+7//w8ssv63y//fv3IzMzE9HR0fqGRGTUsooqEPP1Gdy7+gh+PZ8DiQSYGNkRB14agtiJ3eHjbCt2iERkhHxd7PD18/3RpYMTCspUeHTDH4i/VCR2WG1O7ykhuVyOhIQEBAcH491338XBgwexZ88eHD16FI899hiysrIMHestcUqIjFleaRU+OpiG7SczUaOu++d2X1dvvDQ6DKFejiJHR0TthbKqBs9+dgqnLhXDxsoCnzzRB0PDPcUOq0XaZC8hQRCg0dTtCrt//348+OCDAAA/Pz8UFBToe1sik6GorMH6wxex5eglVNbUdXK+J9Qd80eHo6efs7jBEVG742Rjhc+j++GFL0/jUHI+pm2Nx6pHe2FsTx+xQ2sTeicsUVFReOuttzBy5EgcPnwY69atAwBkZGTAy8s8t8kmAm5sTPjJoYtQVtXtoNzLzxkv3xeOAcHcQZmI9GdrLcWGJ6Mw///+xO4/r2L2zkQoq2owpZ/ufc/aG70TltWrV2PKlCn4/vvv8frrryMkJAQA8M0332DAgAEGC5Covaiu1WDnqUx8eDAN+aV1fYLCvBwwf3Q4RnXx4mZmRGQQ1pYWWP1oLzjZWmLb8Uy8vuscFJU1eGFIsEn/nDH4suaqqipIpVJYWbVdoyvWsJCY1BoBP5zJxqr9KcgqqgQA+LnaImZUGDcmJKJWIwgCVu5LwYcH0wAAzw8Owqv3R7SrpKVNaljqxcfH48KFC5BIJIiIiEBUVFRLb0nULgiCgL1/52LF3mSk5JYBADwcZZg9IhSPRvnB2pKblhFR65FIJHhpdDjktlZ46+cLWH8kHSUVNXh7YneT/EVJ74TlypUrmDx5Mo4ePQpnZ2cAQElJCQYMGIAdO3bAz8/PUDESGZ1jaQV4b08yzmSVAADktnX7/Tw9oBNsrblRGRG1nWn3BMHJxgqvfncWX8VnQVlVg9WP9YLM0rR+Fun9K2B0dDRqampw4cIFFBUVoaioCBcuXIAgCHj22WcNGSOR0TiTVYInNp3A45tO4ExWCWytpJg5LBhHXh6GF4YGM1khIlFM6uuHj6f0hrXUAr+cy8G0rfEoV9WKHZZB6V3DYmtri2PHjiEyMrLB8YSEBAwcOBCVlZUGCbA5WMNCrS01txTv703GnvO5AAArqQRT+gVg5rAQeDjKRI6OiKjO0bQC/PvzeFRUqxHp74wtT/eFs53x7knWJjUs/v7+qKmpaXS8trYWHTt21Pe2REYlq6gCq/an4PvEbGgEwEICPBTpi7kjQ+Hnaid2eEREDQwMcceX0/rh6S2nkJhZgkfXH8cXz94FT6f2v+WH3lNC7733HmbNmoX4+HjUD9LEx8djzpw5eP/99w0WILWcqlaNBd/9hcW7z+NQch6qrjcxo1vLK63Coh/OYfiKQ/guoS5Zua+rN/bMHYwVk3oyWSEioxXp74L/m94fXk4yJOeW4uFP/kBmYYXYYbWY3lNCLi4uqKioQG1trXazwvqv7e0b7iRZVNS6ex5wSuj2fvzzKmbtSNQ+t7GyQP8gNwwN98TQcA8EuJnOzp8tpaiowfoj7E5LRO1fVlEFnth8ApcLK+DpKMPnz96FCG/j+oxskymh1atX6/tSamNJOXU7e3Zys4OqVoNriir8lpyP35LzAQCB7vYYEuaBoeEeuDvIDTZW5lc4yu60RGRq/Fzt8H/T++OpzSeRlFOKR9cfx5Zn+qK3v4vYoenF4I3jxMARltubtvUU9l/Iw+KxXTB1QCck55biUHI+DiXnIf5SMWo1N/4XMLfRF3anJSJTp6ioQfTWUzh9uRi2VlJseKoP7gn1EDssALp9fjNhMQOD3j2IK8WV2Pnc3bg7yK3BudKqGhxNK8Sh5DwcSs5HjrKqwflAd3sMDffA0HBP9At0NZnRF7VGwPeJdd1prxSzOy0RmbaK6lpM35aAIyn5sJJK8MFjkXigewexw2LCQjeUqWrRbdEeAMCZhaNuu7xNEATt6MtvSXk4fdn0Rl/YnZaIzFV1rQbzvj6Dn89eg4UEiJ3YHY/29Rc1pjZtzU/GLTmnFADg5SS741p8iUSCCG8nRHg7YfqQ4OujLwXXp4/qRl9urn0JcrfHkHY0+nL0enfaP9mdlojMkLWlBdY8FgknGyvsOJmJV779CyUVNXh+SLDYoTULExYTV5+whOtRGe5oY4X7unXAfd06QBAEJOXcqH05fbkY6QXlSC8ox5ajl7SjL8MiPDE0zBP+bsaz7PdMVgn+uycJR9MKAQC2VlJED+qE5wYHQ27bdpt0EhGJTWohwdsPdYPc1gqfHL6I2F+SUFJZg5fvDTf6mj2dExYfHx+MHz8e48aNw4gRI2Btbbwd9AhIvr5CKMLbsUX3kUgk6NzBCZ07OOGFocFQVtXg2PXRl9+S85CrVN00+nLeKEZfUnJL8f6eZOz9m91piYjqSSQSvHp/BJztrPDOL0lYd+giFJU1WDa+m1HX7+mcsGzfvh0//vgjZs+ejdzcXNx7770YN24cxowZA1dX19aIkVogqX6ExatlCcs/Od1m9CVe5NGX+u60uxKzIbA7LRFRk6YPqRtlfm3XX9h+IhPKyhqsnNTLaGv5WlR0e/78eezevRs//PADEhMT0b9/f+3oS3Bw282Jsei2aYIgIHLZPpRU1OCnWYPQraO8Td5XWVWDo6nXa19S6kZfbtZaoy95pVVYezANO05mokZd97/1fV298dLoMIQaOGEjIjIVP5+9hrlfJaJGLWBImAc+eaJPm9X1ibJKKDc3F7t378bu3btx4MABBAUF4d1338WYMWMMcfvbYsLStBxFFe6OPQCphQTnl9wryrTMzaMvv12vfVH/Y+XRgGD3uqXTeo6+sDstEVHLHEnJx/NfnEZljRpRAS7Y/HTfNqnxE31Zc0VFBfbs2QNHR0eMHDnS0LdvhAlL0w4l5+HpLacQ7GGPAy8NFTscAM0bfalfNn3XHUZfKqprseXoJaw/zO60REQtdfpyMZ7ZchLKqlp07uCErdF94enYupsmip6wtDUmLE1bf70CfEz3DvhoSm+xw2lEEARcuFaKQyl1Tev+OfpiayVF/2C3RqMv9d1p1xxIQ0EZu9MSERnKhWtKPPXpSeSXqtDJzQ5fPNuvVWv/mLAQACDmqzP4LjEbMaPCMHtEqNjh3JGisr7vS10Ck1f6j9EXD3sMCHbDoeR8dqclImollwvL8cTmE8gqqoSXkwzbnu3XanWATFgIAPDAB7/j72tKrH+yD+7t6i12ODppMPqSlI/TmQ1HX9idloio9eQqq/Dk5hNIyS2Ds50VPnvmLvRqhZpAJiyEWrUGXRbtQXWtBof/M7RdttG/Wf3oy/H0Qvi72mFKvwB2pyUiakUlFdV4esspnMkqgb21FBufisKAEMPWB+ry+c1fTU3UpcJyVNdqYGcthZ9L++89Ire1wgPdO2Dp+G6Ydk8QkxUiolbmbGeNL6f1w6AQd5RXq/HhwTSIOcahV8Jy7do1bNu2Df/73/9QXV3d4Fx5eTmWLl1qkOBIf/UN40K9HGHB2g4iItKDvcwSm5+OwvNDgvDJE31EXdSgc8Jy6tQpdOnSBTNnzsTDDz+Mbt264fz589rzZWVlWLJkiUGDJN3V7yEUwYZpRETUAjJLKRbc3xlyO3H3XtM5YXnttdcwceJEFBcXIzc3F6NGjcKQIUOQmJjYGvGRnupHWCI6MGEhIqL2T+e9hE6fPo2PPvoIFhYWcHR0xEcffYSAgACMGDECe/bsgb+/f2vESTq6sUszExYiImr/dE5YAKCqqqrB85dffhkWFhYYPXo0Pv30U4MERvorV9Uis6gCABDhzVVTRETU/umcsHTr1g3Hjh1Djx49GhyfP38+BEHA5MmTDRYc6Sclt250xcNRBld7a5GjISIiajmda1ieeuopHD16tMlz//nPf7B06VK9poWys7PxxBNPwM3NDXZ2dujVqxdOnz6t833opoJbTgcREZGJ0DlhmTZtGr744otbnn/55ZeRkZGh0z2Li4sxcOBAWFlZ4ZdffsHff/+NFStWwNnZWdfwCDcKbsO5QoiIiEyEXjUshvbuu+/Cz88PW7Zs0R7r1KnTLa9XqVRQqW7sM6NUKlszvHaHBbdERGRq9O50GxcXZ7Agdu/ejaioKDzyyCPw9PREZGQkNm7ceMvrY2NjIZfLtQ8/Pz+DxdLeCYKApJy6BI4Ft0REZCr0Slh++uknTJw40WBBpKenY926dQgNDcWePXswffp0zJ49G59//nmT1y9YsAAKhUL7yMrKMlgs7V1+qQrFFTWwkAChXg5ih0NERGQQOk8JffHFF5g5cya+++47gwWh0WgQFRWFt99+GwAQGRmJ8+fPY926dXjqqacaXS+TySCTyQz2/qakvn6lk5s9bKy43w4REZkGnUZYVq9ejeeeew5ffvklRo4cabAgOnTogC5dujQ41rlzZ2RmZhrsPcwF61eIiMgU6TTCEhMTgw8//BBjx441aBADBw5EcnJyg2MpKSkICAgw6PuYgyQmLEREZIJ0GmEZMGAAPvroIxQWFho0iHnz5uH48eN4++23kZaWhu3bt2PDhg2YOXOmQd/HHCTn1hfcMmEhIiLToVPCsn//fgQGBmLUqFEGXUrct29f7Nq1Czt27EC3bt2wbNkyrF69GlOmTDHYe5gDtUZAam4ZACCcK4SIiMiESARBEHR5gVqtxtSpU3Hx4kX88ccfrRWXTpRKJeRyORQKBZyczPeD+mJ+GUasOAwbKwucX3IfpBYSsUMiIiK6JV0+v3Ve1iyVSrFt2zbcddddegdIrSP5pg63TFaIiMiU6N047oMPPjBkHGQALLglIiJTpXfCcjunTp1qjdvSHSRf73DL+hUiIjI1BktY8vLysHLlSnTr1g133323oW5LOuAuzUREZKpalLCo1Wr88MMPmDBhAvz8/LBx40ZMmDAB8fHxhoqPmqmiuhaXiyoAcEqIiIhMT7Max+Xn52PlypVwdXXF3LlzkZycjC1btmDbtm0AgEmTJkGj0eDbb79t1LGW2kZqbhkEAXB3sIa7A7ctICIi09KsEZbHH38cFRV1v7137NgRd999N65evYpPP/0UV69exYcfftiqQdKdsSU/ERGZsmYlLElJSZgyZQqio6NRVFSE5557DkuXLsWYMWMglXKDPWOgXSHkxYJbIiIyPc1KWN544w089NBDGDJkCN555x1cunQJ3bp1Q79+/bB27Vrk5+e3dpx0B0k5bMlPRESmq1k1LM8//zymTJkCmUwGKysrAHV1Ldu2bcPGjRsxb948aDQa7Nu3D35+fnB05IdmW+OUEBERmbJmrxJycHDQJisA4OHhgXnz5uHPP//E8ePH8cILL2DZsmXw9PTEuHHjWiVYalp+qQqF5dWQSIAwLyYsRERkegzSh6VPnz5Yu3Ytrl69iq1bt6K2ttYQt6Vmqh9dCXC1g601a4qIiMj0GLTTrbW1NSZNmoT//e9/hrwt3UGStsMtR1eIiMg06ZSwTJ8+HVlZWc269quvvsKXX36pV1Ckmxv1K1whREREpqlZRbf1PDw80K1bNwwYMADjxo1DVFQUfHx8YGNjg+LiYvz999+Ii4vDjh074Ovriw0bNrRW3HST5Fy25CciItOmU8KybNkyzJo1C5s3b8Ynn3yCc+fONTjv6OiIkSNHYvPmzRg9erRBA6WmqTUCUnK5QoiIiEybRBAEQd8Xl5SU4PLly6isrIS7uzuCg4MhkUgMGV+zKJVKyOVyKBQKODmZ17RIRkE5hr1/CDJLC/y99D5ILdr+75+IiEgfunx+6zTC8k/Ozs5wdnZuyS2ohZKvF9yGeTkyWSEiIpNl0FVC1PaS2DCOiIjMABOWdq5+hRALbomIyJQxYWnn2JKfiIjMAROWdqyqRo1LheUAmLAQEZFp0zlhOXjwIFvvG4nU3DJoBMDV3hoeDjKxwyEiImo1Oicso0aNQlFRkfb53XffjezsbIMGRc2jbcnv5SjKcnIiIqK2onPC8s+2LefPn4dKpTJYQNR8XCFERETmgjUs7RhXCBERkbnQOWGRSCQNph/++ZzaDkdYiIjIXOjc6VYQBIwYMQKWlnUvraiowNixY2Ftbd3guoSEBMNESE0qLFOhoKxuKi7MiwkLERGZNp0TlkWLFjV4Pn78eIMFQ81XPx3k72oHe1mLdlggIiIyei1OWEgcnA4iIiJzonMNS1VVFXbv3o3S0tJG55RKJXbv3s1VQ22ABbdERGROdE5Y1q9fjw8++ACOjo0/KJ2cnLBmzRps3LjRIMHRrSXl1icst9+Om4iIyBTonLB8+eWXmDt37i3Pz507F59//rlO91y8eLF2tVH9w9vbW9fQzIZGIyA1l1NCRERkPnSuYUlNTUXPnj1veb5Hjx5ITU3VOZCuXbti//792udSqVTne5iLrOIKVFSrYW1pgU5udmKHQ0RE1Op0Tlhqa2uRn58Pf3//Js/n5+frtdeQpaVls0dVVCpVgzoZpVKp8/u1Z/UFt6GeDrCUsvcfERGZPp0/7f45EvJP+/btQ9euXXUOJDU1FT4+PggMDMRjjz2G9PT0W14bGxsLuVyuffj5+en8fu1ZMlcIERGRmdE5YYmOjsayZcvw008/NTr3448/4q233kJ0dLRO9+zXrx8+//xz7NmzBxs3bkROTg4GDBiAwsLCJq9fsGABFAqF9pGVlaXrH6Nd4wohIiIyNzpPCT333HM4cuQIxo0bh4iICISHh0MikeDChQtISUnBpEmT8Nxzz+l0z/vvv1/7dffu3dG/f38EBwdj69atiImJaXS9TCaDTCbTNXSTod2lmSuEiIjITOhVALFt2zbs3LkTYWFhSElJQVJSEsLDw7Fjxw7s2LGjxUHZ29uje/fuehXvmrqqGjUyCsoBcISFiIjMh9493SdNmoRJkyY1Oq5Wq/Hjjz9iwoQJegelUqlw4cIF3HPPPXrfw1Sl5ZVBIwDOdlbwdDTfUSYiIjIvBltikpSUhJdffhk+Pj5NJjK3M3/+fBw+fBgZGRk4ceIEHn74YSiVSkydOtVQ4ZkMbUt+L0fukk1ERGajRQlLeXk5Pv30UwwcOBBdu3ZFQkICli9fjqtXr+p0nytXrmDy5MkIDw/HxIkTYW1tjePHjyMgIKAl4Zmk5Ov1K5wOIiIic6LXlNAff/yBTZs24euvv0ZoaCimTJmCEydOYM2aNejSpYvO99u5c6c+YZilG5sesuCWiIjMh84JS5cuXVBRUYHHH38cJ06c0CYor776qsGDo8bYg4WIiMyRzlNCaWlpGDx4MIYNG4bOnTu3Rkx0C8Xl1cgrrevwy4SFiIjMic4JS0ZGBsLDw/HCCy/A19cX8+fPR2JiIgtA20D9dJCviy0cZHov8CIiImp3dE5YOnbsiNdffx1paWn44osvkJOTg4EDB6K2thafffYZUlJSWiNOAgtuiYjIfLVoldDw4cOxbds2XLt2DWvXrsXBgwcRERGBHj16GCo+uklybn1LfhbcEhGRedE5YZkwYQJ++uknaDQa7TG5XI4ZM2YgPj4eCQkJGDp0qCFjpOuSWHBLRERmSueEpbKyEhMmTICvry9ee+21Ru3ze/XqhTVr1hgsQKqj0QhI4aaHRERkpnROWPbs2YNLly7hhRdewNdff42IiAgMHjwYn3/+OSorK1sjRgKQXVKJ8mo1rKUW6ORuL3Y4REREbUqvGhZfX1+8+eabSEtLw/79+xEQEIAZM2bA29sbzz//PE6cOGHoOM1e/XRQsKcDrKQG21GBiIioXWjxJ9+wYcPwxRdf4Nq1a3jvvffwzTffYODAgYaIjW7CFUJERGTODNLMIz09HZ999hk+++wzKBQKjBw50hC3pZuw4JaIiMyZ3iMslZWV+PzzzzFs2DCEhobiiy++wLRp05CRkYFff/3VkDES2JKfiIjMm84jLMeOHcOWLVvw9ddfo7q6GhMmTMCePXs4qtKKVLVqpBeUA+CUEBERmSedE5ZBgwahZ8+eWL58OaZMmQIXF5fWiItukpZXBrVGgJONJbydbMQOh4iIqM3pnLDEx8ejd+/erREL3UJyzo0Ot9yziYiIzJHONSxMVtoe61eIiMjcsaFHO8AVQkREZO6YsLQDyWzJT0REZo4Ji5FTVNQgR1kFAAhjwkJERGaqRQlLbW0t9u/fj/Xr16O0tG4U4OrVqygrKzNIcAQkXe9w29HZFk42ViJHQ0REJA69O91evnwZ9913HzIzM6FSqTBq1Cg4OjrivffeQ1VVFT755BNDxmm2knM5HURERKT3CMucOXMQFRWF4uJi2Nraao8/9NBDOHDggEGCIxbcEhERAS0YYYmLi8PRo0dhbW3d4HhAQACys7NbHBjV4ZJmIiKiFoywaDQaqNXqRsevXLkCR0d+uBqCIAhIualpHBERkbnSO2EZNWoUVq9erX0ukUhQVlaGRYsW4YEHHjBEbGYvu6QSpapaWEklCPKwFzscIiIi0eg9JbRq1SoMGzYMXbp0QVVVFR5//HGkpqbC3d0dO3bsMGSMZqt+OijYwwFWUq5AJyIi86V3wuLj44MzZ85gx44dSEhIgEajwbPPPospU6Y0KMIl/bHgloiIqI7eCQsA2NraIjo6GtHR0YaKh27CglsiIqI6LUpYUlJScOjQIeTl5UGj0TQ4t3DhwhYFRjeaxrEHCxERmTu9E5aNGzfihRdegLu7O7y9vSGRSLTnJBIJE5YWqq7VID2/HAAQzhVCRERk5vROWN566y0sX74cr7zyiiHjoesu5pehViPA0cYSPnIbscMhIiISld5LT4qLi/HII48YMha6ibZ+xcuxwegVERGROdI7YXnkkUewd+9eQ8YCAIiNjYVEIsHcuXMNfu/2hCuEiIiIbtBpSmjNmjXar0NCQvDmm2/i+PHj6N69O6ysGu4kPHv2bJ2DOXXqFDZs2IAePXro/FpTk8yCWyIiIi2dEpZVq1Y1eO7g4IDDhw/j8OHDDY5LJBKdE5aysjJMmTIFGzduxFtvvXXba1UqFVQqlfa5UqnU6b3agxtLmllwS0REpFPCkpGR0VpxYObMmRgzZgxGjhx5x4QlNjYWS5YsabVYxKaorMFVRRUATgkREREBLahhWbp0KSoqKhodr6ysxNKlS3W6186dO5GQkIDY2NhmXb9gwQIoFArtIysrS6f3M3YpuXWjKz5yG8htre5wNRERkenTO2FZsmQJysrKGh2vqKjQafQjKysLc+bMwbZt22Bj07zluzKZDE5OTg0epoQFt0RERA3p3YdFEIQml9v++eefcHV1bfZ9Tp8+jby8PPTp00d7TK1W48iRI1i7di1UKhWkUqm+YbZL9QW3rF8hIiKqo3PC4uLiAolEAolEgrCwsAZJi1qtRllZGaZPn97s+40YMQJ//fVXg2PPPPMMIiIi8Morr5hdsgLcKLjlCiEiIqI6Oicsq1evhiAIiI6OxpIlSyCXy7XnrK2t0alTJ/Tv37/Z93N0dES3bt0aHLO3t4ebm1uj4+ZAEAROCREREf2DzgnL1KlTAQCBgYEYOHAgLC1btH8i/cM1RRVKq2phaSFBsIeD2OEQEREZBb2zjSFDhhgyjgYOHTrUavc2dvXTQUEe9rC21LsmmoiIyKTwE9HIXGDBLRERUSNMWIwMC26JiIga0ylhOXv2LDQaTWvFQmi4SzMRERHV0SlhiYyMREFBAQAgKCgIhYWFrRKUuapRa3Axv64ZH1cIERER3aBTwuLs7KzdT+jSpUscbTGw9Pxy1KgFOMgs4etiK3Y4RERERkOnVUL/+te/MGTIEHTo0AESiQRRUVG3bOyWnp5ukADNSdL1gtswL4cmuwgTERGZK50Slg0bNmDixIlIS0vD7Nmz8e9//xuOjpy6MBRt/QpXCBERETWgcx+W++67D0DdHkBz5sxhwmJA9QlL5w78OyUiIrqZ3o3jtmzZov36ypUrkEgk6Nixo0GCMldJXCFERETUJL37sGg0GixduhRyuRwBAQHw9/eHs7Mzli1bxmJcPZRW1SC7pBIAEMEpISIiogb0HmF5/fXXsXnzZrzzzjsYOHAgBEHA0aNHsXjxYlRVVWH58uWGjNPkpeTWja54O9lAbmclcjRERETGRe+EZevWrdi0aRPGjRunPdazZ0907NgRM2bMYMKiI+7QTEREdGt6TwkVFRUhIiKi0fGIiAgUFRW1KChzxJb8REREt6Z3wtKzZ0+sXbu20fG1a9eiZ8+eLQrKHHGEhYiI6Nb0nhJ67733MGbMGOzfvx/9+/eHRCLBsWPHkJWVhf/973+GjNHkCYJwUw8WJixERET/pPcIy5AhQ5CSkoKHHnoIJSUlKCoqwsSJE5GcnIx77rnHkDGavFylCorKGkgtJAjxdBA7HCIiIqOj9wgLAPj4+LC41gAuXG/JH+huD5ll01sdEBERmTO9R1jIcDgdREREdHtMWIyAdoUQO9wSERE1iQmLEeAKISIiottjwiKyGrUGF/PKALAlPxER0a0wYRHZpYJyVKs1sLOWwtfFVuxwiIiIjJJOq4QiIyMhkUiadW1CQoJeAZmb+umgMC9HWFg07++WiIjI3OiUsEyYMKGVwjBf9QW3nTuwfoWIiOhWdEpYFi1a1FpxmC1twS1XCBEREd0Sa1hElpxb1zQunAW3REREt6R3p1u1Wo1Vq1bh66+/RmZmJqqrqxuc547Nd1amqkVWUSUA7tJMRER0O3qPsCxZsgQrV67EpEmToFAoEBMTg4kTJ8LCwgKLFy82YIimKyW3bjrI01EGF3trkaMhIiIyXnonLF9++SU2btyI+fPnw9LSEpMnT8amTZuwcOFCHD9+3JAxmiy25CciImoevROWnJwcdO/eHQDg4OAAhUIBAHjwwQfx888/GyY6E6dtyc+EhYiI6Lb0Tlh8fX1x7do1AEBISAj27t0LADh16hRkMplhojNxSTksuCUiImoOvROWhx56CAcOHAAAzJkzB2+++SZCQ0Px1FNPITo62mABmipBEDjCQkRE1Ex6rxJ65513tF8//PDD8PPzw9GjRxESEoJx48bpdK9169Zh3bp1uHTpEgCga9euWLhwIe6//359wzN6eaUqFFfUwEIChHg6iB0OERGRUdNrhKWmpgbPPPMM0tPTtcf69euHmJgYnZMVoG566Z133kF8fDzi4+MxfPhwjB8/HufPn9cnvHahvmFcJ3d72FhJRY6GiIjIuOmVsFhZWWHXrl0GC2Ls2LF44IEHEBYWhrCwMCxfvhwODg4mvdoo+Xr9CqeDiIiI7qxFNSzff/+9AUOpo1arsXPnTpSXl6N///5NXqNSqaBUKhs82psbLflZcEtERHQnetewhISEYNmyZTh27Bj69OkDe3v7Budnz56t0/3++usv9O/fH1VVVXBwcMCuXbvQpUuXJq+NjY3FkiVL9A3dKLAHCxERUfNJBEEQ9HlhYGDgrW8qkTSob2mO6upqZGZmoqSkBN9++y02bdqEw4cPN5m0qFQqqFQq7XOlUgk/Pz8oFAo4ORn/iEWtWoMui/agulaDQ/OHopO7/Z1fREREZGKUSiXkcnmzPr/1HmHJyMjQ96VNsra2RkhICAAgKioKp06dwgcffID169c3ulYmk7XrXi+XCitQXauBrZUU/q52YodDRERk9Ix2t2ZBEBqMopiS+umgMG9HWFhIRI6GiIjI+Ok9whITE9PkcYlEAhsbG4SEhGD8+PFwdXW9471ee+013H///fDz80NpaSl27tyJQ4cO4ddff9U3PKOmXSHkxfoVIiKi5tA7YUlMTERCQgLUajXCw8MhCAJSU1MhlUoRERGBjz/+GC+99BLi4uJuWTxbLzc3F08++SSuXbsGuVyOHj164Ndff8WoUaP0Dc+oJbHgloiISCd6Jyz1oydbtmzRFsoolUo8++yzGDRoEP7973/j8ccfx7x587Bnz57b3mvz5s36htEuJeeyJT8REZEu9F4l1LFjR+zbt6/R6Mn58+cxevRoZGdnIyEhAaNHj0ZBQYFBgr0VXaqMxVZRXYuui/ZAEIDTb4yEm0P7LR4mIiJqCV0+v/UuulUoFMjLy2t0PD8/X9vIzdnZGdXV1fq+hUlKyS2DIADuDjImK0RERM2kd8Iyfvx4REdHY9euXbhy5Qqys7Oxa9cuPPvss5gwYQIA4OTJkwgLCzNUrCaBLfmJiIh0p3cNy/r16zFv3jw89thjqK2trbuZpSWmTp2KVatWAQAiIiKwadMmw0RqIlhwS0REpDu9ExYHBwds3LgRq1atQnp6OgRBQHBwMBwcHLTX9OrVyxAxmhS25CciItJdixvH5eTk4Nq1awgLC4ODgwP0rOE1C4IgaEdYOCVERETUfHonLIWFhRgxYgTCwsLwwAMP4Nq1awCAadOm4aWXXjJYgKYkv0yFovJqSCRAqCcTFiIioubSO2GZN28erKyskJmZCTu7G/vhPProoybbobal6qeDOrnZw9ZaKnI0RERE7YfeNSx79+7Fnj174Ovr2+B4aGgoLl++3OLATJG2foUt+YmIiHSi9whLeXl5g5GVegUFBe16J+XWxBVCRERE+tE7YRk8eDA+//xz7XOJRAKNRoP//ve/GDZsmEGCMzXJLLglIiLSi95TQv/9738xdOhQxMfHo7q6Gi+//DLOnz+PoqIiHD161JAxmgS1RkBK/R5CHYx7+wAiIiJjo/cIS5cuXXD27FncddddGDVqFMrLyzFx4kQkJiYiODjYkDGahMuF5VDVamBjZQF/18ZTaURERHRreo+wZGZmws/PD0uWLGnynL+/f4sCMzX100FhXo6QWkhEjoaIiKh90XuEJTAwEPn5+Y2OFxYWIjAwsEVBmaIkrhAiIiLSm94JiyAIkEgajxSUlZXBxsamRUGZIrbkJyIi0p/OU0IxMTEA6lYFvfnmmw2WNqvVapw4cYJ7CDUhub7g1psFt0RERLrSOWFJTEwEUDfC8tdff8Ha2lp7ztraGj179sT8+fMNF6EJqKxW41JhOQCOsBAREelD54Tlt99+AwA888wz+OCDD+DkxBGDO0nNK4UgAG721vBwZFM9IiIiXem9SmjLli2GjMOkscMtERFRy+idsJSXl+Odd97BgQMHkJeXB41G0+B8enp6i4MzFUnXmLAQERG1hN4Jy7Rp03D48GE8+eST6NChQ5MrhqhOcq4SAFvyExER6UvvhOWXX37Bzz//jIEDBxoyHpN0Y0kz632IiIj0oXcfFhcXF7i6uhoyFpNUUKZCQVk1JBIgzMtB7HCIiIjaJb0TlmXLlmHhwoWoqKgwZDwmp350xd/VDnbWeg9oERERmTW9P0FXrFiBixcvwsvLC506dYKVlVWD8wkJCS0OzhSwJT8REVHL6Z2wTJgwwYBhmK7knOsFtx1Yv0JERKQvvROWRYsWGTIOk1U/JcQVQkRERPrTu4aF7kyjEZCSWwaAPViIiIhaQucRFhcXl2b1XCkqKtIrIFOSWVSByho1ZJYW6ORmL3Y4RERE7ZbOCcvq1atbIQzTVF9wG+rlAKkFG+sRERHpS+eEZerUqa0Rh0nSNozzYsEtERFRSxhFDUtsbCz69u0LR0dHeHp6YsKECUhOThY7rBZjS34iIiLDMIqE5fDhw5g5cyaOHz+Offv2oba2FqNHj0Z5ebnYobUId2kmIiIyDKNovfrrr782eL5lyxZ4enri9OnTGDx4sEhRtUxVjRqXCuoSLo6wEBERtYxRJCz/pFAoAOCWexWpVCqoVCrtc6VS2SZx6SItrwwaAXCxs4KHo0zscIiIiNo1o5gSupkgCIiJicGgQYPQrVu3Jq+JjY2FXC7XPvz8/No4yju7cK0uiQr3dmzWMnAiIiK6Nb1HWGJiYpo8LpFIYGNjg5CQEIwfP17nHZ1ffPFFnD17FnFxcbe8ZsGCBQ3eX6lUGl3ScqPDLVcIERERtZTeCUtiYiISEhKgVqsRHh4OQRCQmpoKqVSKiIgIfPzxx3jppZcQFxeHLl26NOues2bNwu7du3HkyBH4+vre8jqZTAaZzLinWZJzWXBLRERkKHpPCY0fPx4jR47E1atXcfr0aSQkJCA7OxujRo3C5MmTkZ2djcGDB2PevHl3vJcgCHjxxRfx3Xff4eDBgwgMDNQ3LKPBFUJERESGIxEEQdDnhR07dsS+ffsajZ6cP38eo0ePRnZ2NhISEjB69GgUFBTc9l4zZszA9u3b8cMPPyA8PFx7XC6Xw9bW9o6xKJVKyOVyKBQKODmJPwVTVF6N3sv2AQDOLbkXDjKjrG0mIiISlS6f33qPsCgUCuTl5TU6np+fr1214+zsjOrq6jvea926dVAoFBg6dCg6dOigfXz11Vf6hieqpJy6P7+fqy2TFSIiIgPQ+9N0/PjxiI6OxooVK9C3b19IJBKcPHkS8+fPx4QJEwAAJ0+eRFhY2B3vpecgj9FiwS0REZFh6Z2wrF+/HvPmzcNjjz2G2trauptZWmLq1KlYtWoVACAiIgKbNm0yTKTtyI2EhfUrREREhqB3wuLg4ICNGzdi1apVSE9PhyAICA4OhoODg/aaXr16GSLGdocFt0RERIbV4gILBwcH9OjRwxCxmASNRkBKLkdYiIiIDEnvhKW8vBzvvPMODhw4gLy8PGg0mgbn09PTWxxce3SluBIV1WpYW1qgk5u92OEQERGZBL0TlmnTpuHw4cN48skn0aFDB7afv65+hVCIhwMspUa38wEREVG7pHfC8ssvv+Dnn3/GwIEDDRlPu8eCWyIiIsPTewjAxcVF532CzEESW/ITEREZnN4Jy7Jly7Bw4UJUVFQYMp52L5krhIiIiAxO7ymhFStW4OLFi/Dy8kKnTp1gZWXV4HxCQkKLg2tvqmrUyCgoB8CmcURERIakd8JS382WbkjLK4NaI0BuawUvJ+PeTZqIiKg90TthWbRokSHjMAk3Twdx1RQREZHhcN2tASWzYRwREVGr0GmExdXVFSkpKXB3d4eLi8ttRxGKiopaHFx7w5b8RERErUOnhGXVqlVwdKz7MF69enVrxNOuJV9vGscRFiIiIsPSKWGZOnVqk18TUFJRjVylCgAQ5sWEhYiIyJBatPmhRqNBWlpak3sJDR48uEWBtTf100G+LrZwtLG6w9VERESkC70TluPHj+Pxxx/H5cuXIQhCg3MSiQRqtbrFwbUnbMlPRETUevROWKZPn46oqCj8/PPP3PwQLLglIiJqTXonLKmpqfjmm28QEhJiyHjarfqC23B2uCUiIjI4vfuw9OvXD2lpaYaMpd0SBAEpuWUAOCVERETUGvQeYZk1axZeeukl5OTkoHv37o32EurRo0eLg2svrhRXokxVCyupBIHu9mKHQ0REZHL0Tlj+9a9/AQCio6O1xyQSCQRBMLui2/qC22APB1hJ2TyYiIjI0PROWDIyMgwZR7vGlvxEREStS++EJSAgwJBxtGs3Vgix4JaIiKg1tKhxXEpKCg4dOtRk47iFCxe2KLD2hC35iYiIWpfeCcvGjRvxwgsvwN3dHd7e3g36sEgkErNJWFS1alzMLwfAHixEREStRe+E5a233sLy5cvxyiuvGDKedudiXjnUGgGONpboILcROxwiIiKTpPeSluLiYjzyyCOGjKVdSs69MR1k7t1+iYiIWoveCcsjjzyCvXv3GjKWdokt+YmIiFqf3lNCISEhePPNN3H8+PEmG8fNnj27xcG1B8lcIURERNTqJMI/t1pupsDAwFvfVCJBenq63kHpSqlUQi6XQ6FQwMmpbROH/rEHcE1RhW+m90dUJ9c2fW8iIqL2TJfPbzaOawFFRQ2uKaoAAGGcEiIiImo1RtFH/siRIxg7dix8fHwgkUjw/fffix1Ss9R3uO3obAsnG6s7XE1ERET60nuE5eY9hJry6aefNvte5eXl6NmzJ5555hntHkXtQX3DOBbcEhERtS69E5bi4uIGz2tqanDu3DmUlJRg+PDhOt3r/vvvx/33369vKKLhCiEiIqK2oXfCsmvXrkbHNBoNZsyYgaCgoBYFdScqlQoqlUr7XKlUtur73Ur9CiG25CciImpdBq1hsbCwwLx587Bq1SpD3raR2NhYyOVy7cPPz69V368pgiBoa1g4wkJERNS6DF50e/HiRdTW1hr6tg0sWLAACoVC+8jKymrV92vKVUUVSqtqYWkhQZC7Q5u/PxERkTnRe0ooJiamwXNBEHDt2jX8/PPPmDp1aosDux2ZTAaZTNaq73En9QW3wR4OsLY0isVWREREJkvvhCUxMbHBcwsLC3h4eGDFihV3XEFkClhwS0RE1Hb0Tlh+++23W57Lzs5Gx44dm32vsrIypKWlaZ9nZGTgzJkzcHV1hb+/v74htqqka0xYiIiI2opB5zJycnIwa9YshISE6PS6+Ph4REZGIjIyEkDddFNkZCQWLlxoyPAMiiuEiIiI2o7OCUtJSQmmTJkCDw8P+Pj4YM2aNdBoNFi4cCGCgoJw/PhxnZrGAcDQoUMhCEKjx2effaZreG2iulaDi/llADjCQkRE1BZ0nhJ67bXXcOTIEUydOhW//vor5s2bh19//RVVVVX45ZdfMGTIkNaI06ikF5ShViPAUWaJjs62YodDRERk8nROWH7++Wds2bIFI0eOxIwZMxASEoKwsDCsXr26FcIzTvXTQWHejpBIJCJHQ0REZPp0nhK6evUqunTpAgAICgqCjY0Npk2bZvDAjFkS61eIiIjalM4Ji0ajgZXVjZ2JpVIp7O3tDRqUsWPBLRERUdvSeUpIEAQ8/fTT2sZtVVVVmD59eqOk5bvvvjNMhEYoWduDxUnkSIiIiMyDzgnLP7vYPvHEEwYLpj1QVtUgu6QSABDuxREWIiKitqBzwrJly5bWiKPdSLk+utJBbgO5ndUdriYiIiJD4CY4OmJLfiIiorbHhEVHyUxYiIiI2hwTFh1xhRAREVHbY8KiA0EQkJSjBACEe3GFEBERUVthwqKDHGUVlFW1kFpIEOxpXr1niIiIxMSERQf1BbdB7vaQWUpFjoaIiMh8MGHRQdI1FtwSERGJgQmLDpKv16+w4JaIiKhtMWHRQRJb8hMREYmCCUsz1ag1uJhfBoAjLERERG2NCUszZRSUo0YtwN5aio7OtmKHQ0REZFaYsDTTzS35LSwkIkdDRERkXpiwNFN9wS3rV4iIiNoeE5ZmYkt+IiIi8TBhaSbu0kxERCQeJizNUKaqxZXiSgAcYSEiIhIDE5ZmqJ8O8nKSwdnOWuRoiIiIzA8TlmZIZsM4IiIiUTFhaQa25CciIhIXE5Zm0BbcejFhISIiEgMTljsQBAHJuVwhREREJCYmLHeQV6pCSUUNpBYShHg6iB0OERGRWWLCcgcXrtXVr3Rys4ONlVTkaIiIiMwTE5Y7uNHhliuEiIiIxMKE5Q6S2eGWiIhIdEaTsHz88ccIDAyEjY0N+vTpg99//13skACwJT8REZExMIqE5auvvsLcuXPx+uuvIzExEffccw/uv/9+ZGZmihpXrVqDtPwyAOzBQkREJCajSFhWrlyJZ599FtOmTUPnzp2xevVq+Pn5Yd26daLGdamwHNW1GthZS+HnYidqLEREROZM9ISluroap0+fxujRoxscHz16NI4dO9bka1QqFZRKZYNHa6ifDgrzcoSFhaRV3oOIiIjuTPSEpaCgAGq1Gl5eXg2Oe3l5IScnp8nXxMbGQi6Xax9+fn6tEpuisgb21lJOBxEREYnMUuwA6kkkDUcwBEFodKzeggULEBMTo32uVCpbJWmZ0i8Ak/v6o6pWbfB7ExERUfOJnrC4u7tDKpU2Gk3Jy8trNOpSTyaTQSaTtUV4sLCQwM5a9L8mIiIisyb6lJC1tTX69OmDffv2NTi+b98+DBgwQKSoiIiIyJgYxdBBTEwMnnzySURFRaF///7YsGEDMjMzMX36dLFDIyIiIiNgFAnLo48+isLCQixduhTXrl1Dt27d8L///Q8BAQFih0ZERERGQCIIgiB2EC2lVCohl8uhUCjg5MQ9f4iIiNoDXT6/Ra9hISIiIroTJixERERk9JiwEBERkdFjwkJERERGjwkLERERGT0mLERERGT0mLAQERGR0WPCQkREREaPCQsREREZPaNozd9S9c16lUqlyJEQERFRc9V/bjen6b5JJCylpaUAAD8/P5EjISIiIl2VlpZCLpff9hqT2EtIo9Hg6tWrcHR0hEQiETsco6RUKuHn54esrCzut2QE+P0wLvx+GBd+P4xPa31PBEFAaWkpfHx8YGFx+yoVkxhhsbCwgK+vr9hhtAtOTk78AWBE+P0wLvx+GBd+P4xPa3xP7jSyUo9Ft0RERGT0mLAQERGR0WPCYiZkMhkWLVoEmUwmdigEfj+MDb8fxoXfD+NjDN8Tkyi6JSIiItPGERYiIiIyekxYiIiIyOgxYSEiIiKjx4SFiIiIjB4TFhMWGxuLvn37wtHREZ6enpgwYQKSk5PFDouui42NhUQiwdy5c8UOxaxlZ2fjiSeegJubG+zs7NCrVy+cPn1a7LDMUm1tLd544w0EBgbC1tYWQUFBWLp0KTQajdihmYUjR45g7Nix8PHxgUQiwffff9/gvCAIWLx4MXx8fGBra4uhQ4fi/PnzbRYfExYTdvjwYcycORPHjx/Hvn37UFtbi9GjR6O8vFzs0MzeqVOnsGHDBvTo0UPsUMxacXExBg4cCCsrK/zyyy/4+++/sWLFCjg7O4sdmll699138cknn2Dt2rW4cOEC3nvvPfz3v//Fhx9+KHZoZqG8vBw9e/bE2rVrmzz/3nvvYeXKlVi7di1OnToFb29vjBo1SrufX2vjsmYzkp+fD09PTxw+fBiDBw8WOxyzVVZWht69e+Pjjz/GW2+9hV69emH16tVih2WWXn31VRw9ehS///672KEQgAcffBBeXl7YvHmz9ti//vUv2NnZ4YsvvhAxMvMjkUiwa9cuTJgwAUDd6IqPjw/mzp2LV155BQCgUqng5eWFd999F88//3yrx8QRFjOiUCgAAK6uriJHYt5mzpyJMWPGYOTIkWKHYvZ2796NqKgoPPLII/D09ERkZCQ2btwodlhma9CgQThw4ABSUlIAAH/++Sfi4uLwwAMPiBwZZWRkICcnB6NHj9Yek8lkGDJkCI4dO9YmMZjE5od0Z4IgICYmBoMGDUK3bt3EDsds7dy5EwkJCTh16pTYoRCA9PR0rFu3DjExMXjttddw8uRJzJ49GzKZDE899ZTY4ZmdV155BQqFAhEREZBKpVCr1Vi+fDkmT54sdmhmLycnBwDg5eXV4LiXlxcuX77cJjEwYTETL774Is6ePYu4uDixQzFbWVlZmDNnDvbu3QsbGxuxwyEAGo0GUVFRePvttwEAkZGROH/+PNatW8eERQRfffUVtm3bhu3bt6Nr1644c+YM5s6dCx8fH0ydOlXs8Ah1U0U3EwSh0bHWwoTFDMyaNQu7d+/GkSNH4OvrK3Y4Zuv06dPIy8tDnz59tMfUajWOHDmCtWvXQqVSQSqVihih+enQoQO6dOnS4Fjnzp3x7bffihSRefvPf/6DV199FY899hgAoHv37rh8+TJiY2OZsIjM29sbQN1IS4cOHbTH8/LyGo26tBbWsJgwQRDw4osv4rvvvsPBgwcRGBgodkhmbcSIEfjrr79w5swZ7SMqKgpTpkzBmTNnmKyIYODAgY2W+qekpCAgIECkiMxbRUUFLCwafixJpVIuazYCgYGB8Pb2xr59+7THqqurcfjwYQwYMKBNYuAIiwmbOXMmtm/fjh9++AGOjo7aOUi5XA5bW1uRozM/jo6OjeqH7O3t4ebmxroikcybNw8DBgzA22+/jUmTJuHkyZPYsGEDNmzYIHZoZmns2LFYvnw5/P390bVrVyQmJmLlypWIjo4WOzSzUFZWhrS0NO3zjIwMnDlzBq6urvD398fcuXPx9ttvIzQ0FKGhoXj77bdhZ2eHxx9/vG0CFMhkAWjysWXLFrFDo+uGDBkizJkzR+wwzNqPP/4odOvWTZDJZEJERISwYcMGsUMyW0qlUpgzZ47g7+8v2NjYCEFBQcLrr78uqFQqsUMzC7/99luTnxlTp04VBEEQNBqNsGjRIsHb21uQyWTC4MGDhb/++qvN4mMfFiIiIjJ6rGEhIiIio8eEhYiIiIweExYiIiIyekxYiIiIyOgxYSEiIiKjx4SFiIiIjB4TFiIiIjJ6TFiIiIjI6DFhIaIWW7x4MXr16iXa+7/55pt47rnnDHa/hx9+GCtXrjTY/Yio5djplohu605bx0+dOlW727Sbm1sbRXVDbm4uQkNDcfbsWXTq1Mkg9zx79iyGDRuGjIwMODk5GeSeRNQy3PyQiG7r2rVr2q+/+uorLFy4sMEOx7a2tnBwcICDg4MY4WHz5s3o37+/wZIVAOjRowc6deqEL7/8Ei+88ILB7ktE+uOUEBHdlre3t/Yhl8shkUgaHfvnlNDTTz+NCRMm4O2334aXlxecnZ2xZMkS1NbW4j//+Q9cXV3h6+uLTz/9tMF7ZWdn49FHH4WLiwvc3Nwwfvx4XLp06bbx7dy5E+PGjWtwbOjQoXjxxRfx4osvwtnZGW5ubnjjjTdw84Dyxx9/jNDQUNjY2MDLywsPP/xwg3uMGzcOO3bs0O8vjYgMjgkLEbWKgwcP4urVqzhy5AhWrlyJxYsX48EHH4SLiwtOnDiB6dOnY/r06cjKygIAVFRUYNiwYXBwcMCRI0cQFxcHBwcH3Hfffaiurm7yPYqLi3Hu3DlERUU1Ord161ZYWlrixIkTWLNmDVatWoVNmzYBAOLj4zF79mwsXboUycnJ+PXXXzF48OAGr7/rrrtw8uRJqFQqA//NEJE+mLAQUatwdXXFmjVrEB4ejujoaISHh6OiogKvvfYaQkNDsWDBAlhbW+Po0aMA6kZKLCwssGnTJnTv3h2dO3fGli1bkJmZiUOHDjX5HpcvX4YgCPDx8Wl0zs/PD6tWrUJ4eDimTJmCWbNmYdWqVQCAzMxM2Nvb48EHH0RAQAAiIyMxe/bsBq/v2LEjVCoVcnJyDPsXQ0R6YcJCRK2ia9eusLC48SPGy8sL3bt31z6XSqVwc3NDXl4eAOD06dNIS0uDo6OjtibG1dUVVVVVuHjxYpPvUVlZCQCwsbFpdO7uu+9uUDDcv39/pKamQq1WY9SoUQgICEBQUBCefPJJfPnll6ioqGjweltbWwBodJyIxMGiWyJqFVZWVg2eSySSJo9pNBoAgEajQZ8+ffDll182upeHh0eT7+Hu7g6gbmroVtc0xdHREQkJCTh06BD27t2LhQsXYvHixTh16hScnZ0BAEVFRbd9byJqWxxhISKj0Lt3b6SmpsLT0xMhISENHnK5vMnXBAcHw8nJCX///Xejc8ePH2/0PDQ0FFKpFABgaWmJkSNH4r333sPZs2dx6dIlHDx4UHv9uXPn4Ovrq02KiEhcTFiIyChMmTIF7u7uGD9+PH7//XdkZGTg8OHDmDNnDq5cudLkaywsLDBy5EjExcU1OpeVlYWYmBgkJydjx44d+PDDDzFnzhwAwE8//YQ1a9bgzJkzuHz5Mj7//HNoNBqEh4drX//7779j9OjRrfOHJSKdMWEhIqNgZ2eHI0eOwN/fHxMnTkTnzp0RHR2NysrK2zZve+6557Bz507t1FK9p556CpWVlbjrrrswc+ZMzJo1S9sN19nZGd999x2GDx+Ozp0745NPPsGOHTvQtWtXAEBVVRV27dqFf//73633ByYinbDTLRG1a4Ig4O6778bcuXMxefJkAHV9WHr16oXVq1frdc+PPvoIP/zwA/bu3WvASImoJTjCQkTtmkQiwYYNG1BbW2uwe1pZWeHDDz802P2IqOW4SoiI2r2ePXuiZ8+eBrufITdSJCLD4JQQERERGT1OCREREZHRY8JCRERERo8JCxERERk9JixERERk9JiwEBERkdFjwkJERERGjwkLERERGT0mLERERGT0/h8Y4teZCBdSggAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# plot running integral for test data with velocities (PRM_NCBOX, TRJ_NCBOX)\n", + "v_running_integral_plot = v.plot_running_integral()" + ] + }, + { + "cell_type": "markdown", + "id": "0adb88ec-8fc2-44d2-a49a-5d69ed68361e", + "metadata": {}, + "source": [ + "## Exceptions" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "1c7a8eba-420c-46bb-814f-f84b7966779f", + "metadata": {}, + "outputs": [ + { + "ename": "NoDataError", + "evalue": "VACF computation requires velocities in the trajectory", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNoDataError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[18], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;66;03m# Running without velocities raises NoDataError\u001b[39;00m\n\u001b[1;32m 2\u001b[0m u_no_vels \u001b[38;5;241m=\u001b[39m mda\u001b[38;5;241m.\u001b[39mUniverse\u001b[38;5;241m.\u001b[39mempty(\u001b[38;5;241m10\u001b[39m, n_frames\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m5\u001b[39m, velocities\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[0;32m----> 3\u001b[0m v \u001b[38;5;241m=\u001b[39m VACF(u_no_vels\u001b[38;5;241m.\u001b[39matoms)\u001b[38;5;241m.\u001b[39mrun()\n", + "File \u001b[0;32m~/Projects/mdanalysis/package/MDAnalysis/analysis/base.py:448\u001b[0m, in \u001b[0;36mAnalysisBase.run\u001b[0;34m(self, start, stop, step, frames, verbose, progressbar_kwargs)\u001b[0m\n\u001b[1;32m 446\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mframes[i] \u001b[38;5;241m=\u001b[39m ts\u001b[38;5;241m.\u001b[39mframe\n\u001b[1;32m 447\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mtimes[i] \u001b[38;5;241m=\u001b[39m ts\u001b[38;5;241m.\u001b[39mtime\n\u001b[0;32m--> 448\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_single_frame()\n\u001b[1;32m 449\u001b[0m logger\u001b[38;5;241m.\u001b[39minfo(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFinishing up\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 450\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_conclude()\n", + "File \u001b[0;32m~/Projects/transport-analysis/transport_analysis/velocityautocorr.py:187\u001b[0m, in \u001b[0;36mVelocityAutocorr._single_frame\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 180\u001b[0m \u001b[38;5;66;03m# This runs once for each frame of the trajectory\u001b[39;00m\n\u001b[1;32m 181\u001b[0m \n\u001b[1;32m 182\u001b[0m \u001b[38;5;66;03m# The trajectory positions update automatically\u001b[39;00m\n\u001b[1;32m 183\u001b[0m \u001b[38;5;66;03m# You can access the frame number using self._frame_index\u001b[39;00m\n\u001b[1;32m 184\u001b[0m \n\u001b[1;32m 185\u001b[0m \u001b[38;5;66;03m# trajectory must have velocity information\u001b[39;00m\n\u001b[1;32m 186\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ts\u001b[38;5;241m.\u001b[39mhas_velocities:\n\u001b[0;32m--> 187\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m NoDataError(\n\u001b[1;32m 188\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mVACF computation requires velocities in the trajectory\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 189\u001b[0m )\n\u001b[1;32m 191\u001b[0m \u001b[38;5;66;03m# set shape of velocity array\u001b[39;00m\n\u001b[1;32m 192\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_velocities[\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_frame_index] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39matomgroup\u001b[38;5;241m.\u001b[39mvelocities[\n\u001b[1;32m 193\u001b[0m :, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_dim\n\u001b[1;32m 194\u001b[0m ]\n", + "\u001b[0;31mNoDataError\u001b[0m: VACF computation requires velocities in the trajectory" + ] + } + ], + "source": [ + "# Running without velocities raises NoDataError\n", + "u_no_vels = mda.Universe.empty(10, n_frames=5, velocities=False)\n", + "v = VACF(u_no_vels.atoms).run()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}