{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c957bd5f-edbc-48d3-a89b-2a078738c246",
   "metadata": {},
   "source": [
    "# Center of Mass estimator: tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc74374a-7475-4b70-b1c8-433cb80f16fb",
   "metadata": {},
   "source": [
    "If the center_of_mass package is not installed, run the following lines:\n",
    "```\n",
    "import sys\n",
    "sys.path.insert(0, '../src')\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7d65325f-4993-477f-8767-b2f80e9060f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import center_of_mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4e9bf0d9-a7f9-49d1-8deb-9cf7ab0f3282",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy, pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4f09dbb-038f-468f-bc23-100cf632add6",
   "metadata": {},
   "source": [
    "The example data file is loaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a38d719f-73f0-40be-aae7-83879e85a2da",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file     = '..\\examples\\example_data.pkl'\n",
    "pickle_file    = open(input_file,'rb')\n",
    "data           = pickle.load(pickle_file)\n",
    "pickle_file.close()\n",
    "\n",
    "sex                 = data['sex']                # 'female' or 'male'\n",
    "Labels              = data['Labels']             # list of marker labels\n",
    "Kinematic_frequency = data['Position_frequency'] # in Hertz\n",
    "Position            = data['Position']           # dictionary with, for each marker, the position in meters (numpy.ndarray of shape (3,duration*Position_frequency)\n",
    "Force_frequency     = data['Force_frequency']    # in Hertz\n",
    "GroundReactionForce = data['GroundReactionForce']# in Newtons (numpy.ndarray of shape (3,duration*Force_frequency)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd605c35-07f2-433a-a2a1-fe0ffed44340",
   "metadata": {},
   "source": [
    "The center of mass is calculated (in this example we use the method of Dumas et al. 2007)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7d519aa-0b2b-4e29-8e60-01ff4a4ea27f",
   "metadata": {},
   "outputs": [],
   "source": [
    "kinematics    = center_of_mass.Kinematics(Position, Labels, sex)\n",
    "Kinematic_com = kinematics.calculate_CoM()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b706e15-c7d6-4a00-b71b-6e0f8c50382e",
   "metadata": {},
   "source": [
    "The person's mass is determined as the median vertical ground reaction force during the initial 1.5 seconds of quiet standing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f8d445eb-a837-4087-b1d7-d9a96231b5cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mass = numpy.median(GroundReactionForce[2,:int(1.5*Force_frequency)])/9.81"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e88759-bde1-4f21-8c19-d7695c8b49d2",
   "metadata": {},
   "source": [
    "The position and velocity estimates are calculated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "71ac05df-5ea7-4bd0-a68f-5c7bcf12cef1",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = (GroundReactionForce, Force_frequency, Kinematic_com, Kinematic_frequency, mass)\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args)"
   ]
  },
  {
   "attachments": {
    "afed88e6-49cc-4d93-9895-86132a5b8436.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABMYAAAHhCAYAAACFnLLlAAAgAElEQVR4Xuzdf5QV1Z33+683mQQVI79iJJrY0jAkF1wcFMGJrdA+OrAgxsbMGuiszNxuNAFxXKtbsiQLdYBRSHCp3c8YESdK97PGOw25K9LGwMDElQbBNfJLDle5CYHGNoPBcfghT9Awz2OWt7/H6bbrUOecqjp1qnbteu9/OpH6sffrW3RxPmfXrvM+6m1CQwABBBBAAAEEEEAAAQQQQAABBBBAIGUC5xGMpaziDBcBBBBAAAEEEEAAAQQQQAABBBBAICdAMMaFgAACCCCAAAIIIIAAAggggAACCCCQSgGCsVSWnUEjgAACCCCAAAIIIIAAAggggAACCBCMcQ0ggAACCCCAAAIIIIAAAggggAACCKRSgGAslWVn0AgggAACCCCAAAIIIIAAAggggAACBGNcAwgggAACCCCAAAIIIIAAAggggAACqRQgGEtl2Rk0AggggAACCCCAAAIIIIAAAggggADBGNcAAggggAACCCCAAAIIIIAAAggggEAqBQjGUll2Bo0AAggggAACCCCAAAIIIIAAAgggQDDGNYAAAggggAACCCCAAAIIIIAAAgggkEoBgrFUlp1BI4AAAggggAACCCCAAAIIIIAAAggQjHENIIAAAggggAACCCCAAAIIIIAAAgikUoBgLJVlZ9AIIIAAAggggAACCCCAAAIIIIAAAgRjXAMIIIAAAggggAACCCCAAAIIIIAAAqkUIBhLZdkZNAIIIIAAAggggAACCCCAAAIIIIAAwRjXAAIIIIAAAggggAACCCCAAAIIIIBAKgUIxlJZdgaNAAIIIIAAAggggAACCCCAAAIIIEAwxjWAAAIIIIAAAggggAACCCCAAAIIIJBKAYKxVJadQSOAAAIIIIAAAggggAACCCCAAAIIEIxxDSCAAAIIIIAAAggggAACCCCAAAIIpFKAYCyVZWfQCCCAAAIIIIAAAggggAACCCCAAAIEY1wDCCCAAAIIIIAAAggggAACCCCAAAKpFCAYS2XZGTQCCCCAAAIIIIAAAggggAACCCCAAMEY10CkAj09PVJVVRXpOTkZAggggAACfQI9H7wrVRdcAggCCCCAAAIIIIAAAjmB2IOxp59+WhYsWJDrzNChQ6Wrq0smTJhgRHl27Ngh3/jGN+TUqVO5/qxZs0bmz59vRN+S2In33ntPJk6cKG+++WYSu0+fEUAg5QInTpyQKVOmSHd3d07irrvuktWrVxuhkt+3OXPmyLp164zom0mdeO9/vy8Tt94rb97ytEndoi8IIIBA6AImf8bK79vPfvYzqampCd2AAyKAAAJeBWINxvQf8mPGjOkPnrTT+oFj1KhRXvtf0e2OHDki1dXV/efQ4O7QoUMyfPjwip7X1oMvW7ZMli9fLkuXLhX93zQEEEAgSQKPPPKILF68uL/LJgVj2qnRo0f3h3b6/7dv384HjbwLbNmv18nyg+tl6dg5suwrc5N0+dFXBBBAwLOA6Z+x9Iub+vr6/vFMnz5dNm/e7Hl8bIgAAgiELRBrMGb6hwzFXrhwoTz11FP97qtWrZL77rsv7DpYfzydLXbllVeK/hwyZIgjDLV+8AwQAQQSL2D6hwwF5oNG8ctMZ4td+Yv5oj+H/MmFcmrmc4m/LhkAAggg4CaQhM9YfJnDtYsAAiYJxBaMuX3IyGazxjxG2VckfZzyhhtu6K8Zs8aCXb59s8X69mbWWDBH9kIAgXgE8j9kBHlUUWch79q1Szo7O3Mzu/bs2eO4t0yePFluuukm+Yu/+ItAM6fd7qvMGvvkeumbLdZ/H2LWWDx/mTgrAghUVCCsz1ibNm2SV155Rfbu3Zu7d/UtLaOdnzRpklx77bVy4403yi233BLoaZr8+yqzxip6WXBwBBAoIRBbMJb/zbb+gt29e3coBdNftD/84Q/7f4GX+4s2/xuNjo4OmTuXRzC8FmvgbLG+fZg15lWP7RBAwASB/PvAxo0bZebMmZ66pl+w/OhHP5L169d72l430sc0H3roId8fNvJnOZv2uKdngJA3HDhbrP8+xKyxkJU5HAIImCBQzmcsDdWeffZZx+eoUmPSSQM/+MEPfK/DnL9kjZ7HpCV1So2bP0cAAbsEYgvG9FuGgd+Wh7GwvX74aG5udhy3r1wfffRR4MoNXCBSDxJmiBe4UwnaMX+2WF/XmTWWoCLSVQRSLKDfms+aNatfQD8EnDx5sqSIfsC4++67fQViAw8a5IU0+/fvl0wm4+jb8ePHfQdsJQeXsA3yZ4v134eYNZawStJdBBAoJRD0M5be67797W8HXu4kyBcxOtFg4JdGS5YskRUrVpQaIn+OAAIIhC4QSzAW9j/c9cPH448/LitXriwIVE4w5vaNhomPfYZ+dYRwQLfZYn2HZdZYCMAcAgEEKi4Q9B/uM2bMkC1btpTVPw3H9EskPy+lyZ/dFsYXT2UNIuad3WaL9d+HmDUWc3U4PQIIhCkQ9DOW2359/dL7kD7qr00/cw2c2JDfd7/BVv7sNq9fPIVpxrEQQAABFYglGMt/1CPIWi195dNZYg0NDY43cbmVtpxgTI+X/+2L31/8ab3cCs0W6/Ng1lharwzGjUAyBMr5YuS8885zDFL/wa8h27e+9S356le/2j+LSz+Q/OQnPyn45Y7f5QDuv/9+x7HSPsu50Gyx/vsQs8aS8ZeRXiKAQEmBoJ+x8tf70hPpDLD58+efs/5z3+OWA9/SPLBjfh6H1GONGDHCMS4/SxWUBGEDBBBAwKNALMHYsGHDHNN0g3yb7WWW2ECDcoOx/Mcpq6ur5fDhwx6Z07lZsdlifSLMGkvntcGoEUiKQP7vfj/fZvcFY7rP97//fbnjjjuKPtKoAVltba3rYyx+Zim7ffPv54NKUmrjpZ/FZov134eYNeaFkm0QQCABAkE/Yw0MxnTCgj6FU2qmcv4yA308ficP5M/KDvJIZgJKQxcRQMBwgciDMbdfon7/wV5slpj+Ms9/25fWoNxgzO2Dhp8PKoZfBxXpXqnZYn0nZdZYRfg5KAIIhCCQP1vYzz/Y9VFKbatXry75AaOvq2F90Aj64SgEMqMOUWq2WP99iFljRtWNziCAgH+Bcj5j6Werb3zjG/Lcc895frGM9jB/hpr+N7+TB8r5Asq/EnsggAAC7gKRB2P5j3j4/eXp9lhL3y/h9vZ2qampEbd1XcoNxvQc+eu2+P1GJE0XoZfZYn0ezBpL05XBWBFIjoDbIx5RvJU4P4xTMb+PU+Z/WEnj45ReZov134eYNZacv5j0FAEEXAXK/YwVhLXQ2mR+Pne5fbbjccog1WAfBBAoRyDyYCw/XPLz7bsOVL/RuOGGG/rH3PeIyn333df/3yoVjPFBw/ul5nW2WN8RmTXm3ZYtEUAgGoH8RYH1rH5nOAfpqdtaL36/RHLre9reTul1tlj/fYhZY0EuV/ZBAAFDBMr9jBV0GPnraepx/ARjuj2TD4Lqsx8CCIQlEGkw5vaNgN9v3wcGYxqqPfTQQ+es2VKpYCx/qm9UH5LCKnZUx/EzW6yvT8wai6o6nAcBBLwK5H8Z4jec8nqe/O3cgjG/HzTcvsX3e78N2n8T9vMzW6z/PsSsMRNKRx8QQCCAQBifsQKcNrdLGMFY/jpjUd1vg46Z/RBAwD6BSIMxt2+wg6zTpR8avva1r+Uem3RrlQrG0v5Bw+vl73e2WN9xmTXmVZjtEEAgCoH8b7DLeYOyn/66BWNBHoXMX2fM7wxtP302bVu/s8X670PMGjOtlPQHAQQ8CIT1GcvDqc7ZJIxgjMkHQeTZBwEEwhSINBjLf/ZdB+J3qq2XwVcqGNNz5//yZ50xZ0WCzBbrOwKzxrxc3WyDAAJRCLitL7Zq1SoZ+Nh+pfrhFoz5XWNM+5Z/LwwSrlVqjJU8bpDZYv33IWaNVbI0HBsBBCokENVnrPzuu00aCHKvyV8qR8/DOmMVulg4LAIIuApEGozl/yM9yD/0vdSxksFY/qLIlRqDl3GauE3Q2WJ9Y2HWmIlVpU8IpE/A7R/p27dvLzhTOUwht7d8BZntFdcHpTAtghwr6Gyx/vsQs8aCsLMPAgjEKBDVZ6z8Ibrdr4JMGojzy6gYy8apEUDAIIFIg7H82VaVeiylksFYJY9t0HURqCvlzBbrOyGzxgLRsxMCCIQs4DZrK8ij/0G6lf8Ipx4jyPpgbmOIKtwLMu4w9ilntlj/fYhZY2GUgmMggECEAlF9xho4pEJvpAz6kpr8MTD5IMILiFMhgIBEFoy5LQpZqcdSKhleuX3QSNubvgr9vSl3tljfcZk1xm8mBBCIW8DtW/BKPPqfP85CHzSC3GfS+GhKubPF+u9DzBqL+68g50cAAY8CUX7G6uuS3qtqa2vl1KlTjl6W89kurY//eywzmyGAQIUFIgvG3P6BXs4vz2IuUQdjtn8D7+UaDGO2WN95mDXmRZxtEECgkgKVvI8U63dYj1HqOaK871ayFl6PHcZssf77ELPGvLKzHQIIxCwQ5e96DcR0ofynnnrqnFEHeeR/4EHiuu/GXD5OjwAChghEFoy5vS2lUoFSJX+xbtq0SWbNmuUoH4tDioQ1W6wPllljhvyGoBsIpFQg/3HGKB7pKDRbLOi90m3NliBrvyTlEghrtlj/fYhZY0kpPf1EINUCYX/G0uO1t7efY7plyxZX56FDh8oPfvADmT9/fll1cFsXM8hs6bI6wc4IIJBagciCsSjXOqlkMBbltzJJuipbW1tFZ425teXLlxccigZgbk1njTU1NSWJgL4igIBFAnGsdeJ27yr3G/g4xhHXZdDa/aLorDHX+9DB9YXvQ70BmOt9qHfWWFP1rXENh/MigAACngTC/ozldi9y60h1dbV897vflTvuuEOGDx/uqa/FNgp7HGV3iAMggECqBAjGfJabYMwnWO/m+R/MBh4hijV7/PeYPRBAIO0CUQdK+mjKggULHOz6LfyhQ4fK+sAR9ThMvW7Oe2F2wa59dNsGU7tNvxBAAIGSAmEHSl6Dsb6OhRWQhT2OknBsgAACCAwQIBjzeTkQjPkEIxjzD8YeCCAQu0CUgVKhRYzDeEw/ynHEXrQiHSAYM7k69A0BBMoRCDtQcjuel/7plznPPfeczJw508vm52wT9jgCdYKdEEAgtQIEYz5LTzDmE4xgzD8YeyCAQOwCUQVKug6Yfju/Z88ex5jDWgssqnHEXrASHSAYM71C9A8BBIIKRBUo6f3qV7/6lfzzP/+zrF+/Xrq7u1273NHRIXPnzvU9nKjG4btj7IAAAqkQIBjzWWaCMZ9gBGP+wdgDAQRiF4gqUHJ7ZGXSpEmye/fuUAyiGkcona3gQQjGKojLoRFAIFaBuAIlXaRf36R86tSpc8Yf5KUxcY0j1uJxcgQQMEaAYMxnKQjGfIIRjPkHYw8EEIhdIIpAST9Q5L/yXtdq2blzZ1nrig3Ei2IcsRfLQwcIxjwgsQkCCCRSIM5AqdBSAHovO3z4sC/POMfhq6NsjAACVgoQjPksK8GYTzCCMf9g7IEAArELVDpQcgvFdH2Wrq4umTBhQmjjr/Q4QutohQ9EMFZhYA6PAAKxCcQdKOnMsfr6+nPG73edzLjHEVsBOTECCBghQDDmswwEYz7BCMb8g7EHAgjELlDJQKnQGyjDDsUUsZLjiL1IPjpAMOYDi00RQCBRAiYESqNHjz5nzTG/a2WaMI5EFZ7OIoBAqAKRBWNuHwT8fpPgdeRua7Z89NFHXncvut2mTZtk1qxZjm0qNY5QOmzAQfI/mA3sUlh1MWCYdAEBBCwSGDZsmGPdlLDW/XL7Zr0SM8W0FLpQ8ogRIxxV8ftBxZaSEozZUknGgQAC+QJRfsYqpO82C3r69OmyefNmzwW7//77ZeXKlY7tjx8/HtrSAp47woYIIJBKgciCsShnWlUyGOPbDP9/TwjG/JuxBwIIxCtQiftIocdNgr7Bq5RQlPfdUn2J+88JxuKuAOdHAIFKCZjwu97t85HfYKwS991KmXNcBBCwTyCyYEwXZ8xkMg7BVatWyX333Re6aiV/sbr94tfXFY8aNSr0cdhyQIIxWyrJOBBIj4C+al5fRz+wlTPDNepQTPvt9mGpUiGc6VcGwZjpFaJ/CCAQVCDKz1iF+liJYCzIAv5BDdkPAQQQiCwYU+qo1jqpZDBWyWPbejkSjNlaWcaFgL0CYc4OdvvQonKVDqnCHEPSK00wlvQK0n8EECgmENVnrEJ9cPsyye8EiLjHwBWGAALpFog0GMtfmDGsNVvyS1jJ8Cp/DH6nCafxciMYS2PVGTMCyRZwm+EVZD3JQq+yr3Qopvpua76UM+styRUlGEty9eg7AgiUEojqM5ZbP9zWs9Tt/Nwzjxw5IjpDbGDzG6yVMuLPEUAAgaJfMPT+IzmcVek9OIf9aEqhU1YyGMsPee666y5ZvXq1h9GndxOCsfTWnpEjkFSBMB5NKRSKrVmzRubPn19xmvx7YZofSyEYq/jlxgkQQCBGgXI+Y/Xdq5577jmZOXOm71G4LZqvL5U5efKk52Px6L9nKjZEAIEKCUQ6Y8ztrSnbt2+XmpqaUIdXqWCMX9rBykQwFsyNvRBAIF6B/DdT+pkhrN+gjxkzxvFmSx1NlF+m8EXOJ9cPwVi8f5c4OwIIVFagnM9YAx+71/vcAw884Pmzmdt5daR+Z3uxhnNlrw+OjgACpQUiDcbcvoEP8s25fjOxd+/egqPbtWvXOR9G9Be9WxsyZIg8+eSTnl4F7PbLP5vNyoQJE0pLp3gLgrEUF5+hI5Bggfxv4P18A17ow0K5HF6/TArrfltuf03Zn2DMlErQDwQQqIRAOb/zCy2c39DQILfccovrZ6RNmzbJ3//938uWLVvOGY4ulbN582ZPn636di7nflsJT46JAALpE4g0GFPe/G/g58yZI7qWi59WLGjxc5y+bb1+0Mj/pZ3mx1L8OCc5GCs068Pr+Cu1jp7X87MdAggEFyjnyxC3DxrBe/LJnl6/hXfre5rfoJzUYOzrrz4sG//9ky8Cr754lOyd9ligS+marYvktdNH+vdtqr5VWsbPC3QsdkIAAfMEgn7GcnsiZuDo9POOrmHW19zCsL4/0y+Qurq6fE8ayO97lLOrzaskPUIAgTgEIg/GwvhGIK5gjF/awS7RJAdjpf6x4EWEWYVelNgGAfMEwv4GPowReg3G+CLHqZ3UYGzqjgfk5RMHHIN5dFyDLBp9m6/L6fljr8o3d61y7HPnFbfIjzMLfR2HjRFAwFyBcj5jub2sxe9I9cvg9evXy6hRo3zt6navjeIFNb46ycYIIGC9QOTBmNubvvwGB9dee63s2bMnlOLoNxt6rFK/xN1+aft520oonU3oQZIcjLnV3W8Z/F7ffo/P9gggUDmB/Dd9eZ3lrI+ZfPvb3z7nsf5ye+p1hnP+FzlLliyRFStWlHv6xO6f1GDsscMvyPcOtDvcZ33hGvn5dQ/4qkX9nsdk3ds7HPv8dPJiuX3kdb6Ow8YIIGCuQLmfsfTL4B/96Ee5cMtP0xllDz/8sGgwF6S5zXA+fvy4r0cxg5yXfRBAAIGBApEHY26v9PX6DXicpct/LMbPWjNx9tuEcyc5GFM/vdG/9957nijd1reL8MWvnvrIRggg4F3A7W1bpv+ddpvpmvaAPqnBmF6pF22slzMfnu2/aAd/epD8flaH94u4d8svbpknx86e6t9n9IUj5dDNvFHbFyIbI2C4QFifsfQ4v/jFL+T111/Premc/+hk36OVN910k3zta1/zvFB/Ib78l6Z5/QLK8HLQPQQQSJhA5MFYX9Aw8NuIJKzDlD9LjWffvV/pSQ/GvI7UbXYZ14lXPbZDwEyBJM4Wzg/zWA9TJMnBmNtsr78d+5ey/Cv1nv7SuD1GyfpinujYCIHECeQ/Tmn6Zyy3MI/HKBN32dFhBKwQiCUY00dMZs2a5QA0eVHgI0eOiH6wGNh4jNL79Z+WYCw/PNVZhYcOHWIquPdLhS0RMFIgaV+M5D/+mfbHKPWiSnIw5hZs3Th8nGyredjT3xe3YE331WPQEEDALoGkfcZye/yTxyjtuiYZDQJJEYglGFOc/H+4m/w4Zf5jlHz77u/yTkMw5vYGOsJTf9cJWyNgqkD+P9xNDr3dHqM0+YunqGqe5GBMjca8tFAOv3/MwXX45qek+sJLSxLmP0bpJ1QreXA2QAAB4wSS9Bkr/zFKnrQw7nKiQwikRiC2YCxJYVOSbjAmXrm2B2NuMwqnT58umzdvNrEc9AkBBHwK6KMeY8aMcSykb+qjHvlvFuN30cfFTnow1vzGWmntftFx5Xp5nPLZt16SO7NP+t7P518RNkcAAYMEkvIZy+3fz15fMGMQN11BAAFLBGILxtyeKTfxl2H+t+8mzxQw9Zq0PRjL/7bL65tOTa0X/UIAgXMF8tftMjFwcgvwTLyvxnF9JT0Y637/HRn90l0OOi8L6H/91Ydl47/v7d9PF+7PTmvxNNMsjjpxTgQQKF8gKZ+x8gM8E++r5VeDIyCAQFIEYgvGFCj/m20Tp8/mL2LJWi3+L22bgzG3V0yb/Fiw/+qxBwIIqIDbN9umPaKY//vI9EWXo7yykh6MqdXUHQ/IyycOONhKrRWW/0bLWV+4Rn5+3QNR0nMuBBCIQSAJn7GGDRvmmInNEiQxXCicEgEE+gViDcZM/6CR3z9miwX7m2NrMOY2O4MPosGuEfZCIAkCpn/QyH/sn9lin1xVNgRjjx1+Qb53oN3xV6XY2yXdtn90XIMsGn1bEv660UcEEChDwPTPWPlrdzJbrIxisysCCIQiEGswpiMY+A23Bk8/+9nPpKamJpTBlXuQ/Mco16xZI/Pnzy/3sKnb39ZgLH82oRY2m83KhAkTUldjBoxAGgQ0DJ8yZYroTDFtc+bMEf3HvQktP6g3qW8m+NgQjKlj/gywkYOGyu+mr3Ulzn+Msti2JtSIPiCAQLgCJn/GGvgYpWmf/8KtAkdDAIGkCMQejCUFin4GF7AxGHN7HTaP2Qa/RtgTAQQQqKSALcFY/Z7HZN3bOxxUP528WG4fed05fPkh2tzLaqRj0qJKMnNsBBBAAAEEEEAgkQIEY4ksW7I6bVswlj9rRKtRXV0tO3fulOHDhyerOPQWAQQQSIGALcHY88delW/uWuWomFvg5fYYZaEALQXlZ4gIIIAAAggggEBRAYIxLpCKC9gWjOW/nU4BWcun4pcRJ0AAAQQCC9gSjCnAmJcWyuH3j/VbuD0imf8YpZc3WAbGZUcEEEAAAQQQQCDhAgRjCS9gErpvUzCWv+6c+rOWTxKuQvqIAAJpFrApGGt+Y620dr/oKGf+ovr5j1EWW6Q/zdcFY0cAAQQQQAABBFSAYIzroOICNgVj+W99402lFb98OAECCCBQtoBNwVj3++/I6JfucpjM+sI18vPrHsj9t6W/7pC/O/gTx59vq3lYbhw+rmxHDoAAAggggAACCNgoQDBmY1UNG5MtwdjAN+j0EXd0dIi+nZKGAAIIIGCugE3BmCpP3fGAvHziQD/44E8Pkt/P6sj9//w/00BMgzEaAggggAACCCCAgLsAwRhXRsUFbAjG9u/fL5lMxmE1ffp02bx5c8X9OAECCCCAQHkCtgVjbovr/+3Yv5S//lLtObPJ9L8v/0p9eYDsjQACCCCAAAIIWCxAMGZxcU0Zmg3B2IwZM2TLli39pPoI5Z49e2TUqFGmMNMPBBBAAIECArYFYzrM/HXEdGbYtBHjHI9R6kyy7LQWqb7wUq4NBBBAAAEEEEAAgUL/Vvyot6GDQCUFkh6MPf3007JgwQIH0apVq+S+++6rJBvHRgABBBAIScDGYKx+z2Oy7u0dDqGrLx4lr50+0v/fBq49FhIlh0EAAQQQQAABBKwTYMaYdSU1b0BJDsZOnDghY8aMkVOnTvXDTpo0SXbv3m0eND1CAAEEEHAVsDEYe/7Yq/LNXauKVjz/bZVcHggggAACCCCAAALnChCMcVVUXCDJwZjbbLHq6mrRt1N6bQ0NDSzQ7xWL7RBAAIEKCNgYjCnTmJcWyuH3j7mKjRw0VH43fW0FNDkkAggggAACCCBglwDBmF31NHI0SQ7G3N5E6ReZRfr9irE9AgggEK6ArcFY8xtrpbX7RVesuZfVSMekReFCcjQEEEAAAQQQQMBCAYIxC4tq2pAIxnh7pWnXJP1BAIF0CdgajHW//845b6Hsq+xPJy+W20del65CM1oEEEAAAQQQQCCAAMFYADR28SeQ5GBsx44d8o1vfMOxxpi/0YssWbJEVqxY4Xc3tkcAAQQQCEnA1mBMedwW4R994Ug5dPPqkPQ4DAIIIIAAAgggYLcAwZjd9TVidEkOxowApBMIIIAAAmUJ2ByMPfvWS3Jn9kmHT1P1rdIyfl5ZZuyMAAIIIIAAAgikRYBgLC2VjnGcBGMx4nNqBBBAAAGxORj7+qsPy8Z/3+uo8raah+XG4eOoPAIIIIAAAggggIAHAYIxD0hsUp4AwVh5fuyNAAIIIFCegK3BmNsaYxqIaTBGQwABBBBAAAEEEPAmQDDmzYmtyhAgGCsDj10RQAABBMoWsDUYc3sr5d+O/UtZ/pX6ss04AAIIIIAAAgggkBYBgrG0VDrGcRKMxYjPqRFAAAEErH2U8otb5smxs6ccFT5881NSfeGlVB0BBBBAAAEEEEDAowDBmEcoNgsuQDAW3I49EUAAAQTKF7BxxpjbovuzvnCN/Py6B8oH4wgIIIAAAggggECKBAjGUlTsuIZKMBaXPOdFAAEEEFABG4Mxt0X3Hx3XIItG30bREUAAAQQQQAABBHwIEIz5wGLTYAIEY8Hc2AsBBBBAIBwB24IxXXQ/s7VZznx4th9o8KcHye9ndYQDxlEQQAABBBBAAIEUCRCMpajYcQ2VYCwuec6LAAIIIKACtgVjS3/dIX938CeO4s69rEY6Ji2i4AgggAACCNZAc2IAACAASURBVCCAAAI+BQjGfIKxuX8BgjH/ZuyBAAIIIBCegO3B2I3Dx8lDX60X/UlDAAEEEEAAAQQQ8CdAMObPi60DCBCMBUBjFwQQQACB0ARsC8ZCg+FACCCAAAIIIIAAAkIwxkVQcQGCsYoTcwIEEEAAgSICBGNcHggggAACCCCAAAKFBAjGuDYqLkAwVnFiToAAAgggQDDGNYAAAggggAACCCAQQIBgLAAau/gTIBjz58XWCCCAAALhCjBjLFxPjoYAAggggAACCNgkQDBmUzUNHQvBmKGFoVsIIIBASgQIxlJSaIaJAAIIIIAAAggEECAYC4DGLv4ECMb8ebE1AggggEC4AgRj4XpyNAQQQAABBBBAwCYBgjGbqmnoWAjGDC0M3UIAAQRSIkAwlpJCM0wEEEAAAQQQQCCAAMFYADR28SdAMObPi60RQAABBMIVIBgL15OjIYAAAggggAACNgkQjNlUTUPHQjBmaGHoFgIIIJASAYKxlBSaYSKAAAIIIIAAAgEECMYCoLGLPwGCMX9ebI0AAgggEK4AwVi4nhwNAQQQQAABBBCwSYBgzKZqGjoWgjFDC0O3EEAAgZQIEIylpNAMEwEEEEAAAQQQCCBAMBYAjV38CRCM+fNiawQQQACBcAUIxsL15GgIIIAAAggggIBNAgRjNlXT0LEQjBlaGLqFAAIIpESAYCwlhWaYCCCAAAIIIIBAAAGCsQBo7OJPgGDMnxdbI4AAAgiEK0AwFq4nR0MAAQQQQAABBGwSIBizqZqGjoVgzNDC0C0EEEAgJQIEYykpNMNEAAEEEEAAAQQCCBCMBUBjF38CBGP+vNgaAQQQQCBcAYKxcD05GgIIIIAAAgggYJMAwZhN1TR0LARjhhaGbiGAAAIpESAYS0mhGSYCCCCAAAIIIBBAgGAsABq7+BMgGPPnxdYIIIAAAuEKEIyF68nREEAAAQQQQAABmwQIxmyqpqFjIRgztDB0CwEEEEiJAMFYSgrNMBFAAAEEEEAAgQACBGMB0NjFnwDBmD8vtkYAAQQQCFeAYCxcT46GAAIIIIAAAgjYJEAwZlM1DR0LwZihhaFbCCCAQEoECMZSUmiGiQACCCCAAAIIBBAgGAuAxi7+BAjG/HmxNQIIIIBAuAIEY+F6cjQEEEAAAQQQMENg+5HTMud//H9y7H/+r/4OXfOli2TPvVd77uB3f/Ib+fG/Huvffsznz5ffLJnseX8bNiQYs6GKho+BYMzwAtE9BBBAwHIBgjHLC8zwEEAAAQQQSLHAn67cJYf+4w8OgXunXS6P3VZdUmXtznfkjnUHHdvVX32J/NNffbXkvjZtQDBmUzUNHQvBmKGFoVsIIIBASgQIxlJSaIaJAAIIIIBACgXcwq3Bn/2UbPruVXLDqIuLinxx6b86ZpuN/NxnZPs9GakecX6qJAnGUlXueAZLMBaPO2dFAAEEEPhYgGCMKwEBBBBAAAEEbBb41j/+Sjpee9cxxKnVQ2Tr30woOOz8Ryh1w8frqqV56uU2U7mOjWAsdSWPfsAEY9Gbc0YEEEAAgU8ECMa4GhBAAAEEEEDAZoHu43+QG57IOmZ/FQu6dG2yG3u3H9hKBWk2+xGM2VxdQ8ZGMGZIIegGAgggkFIBgrGUFp5hI4AAAgggkCKBlm1H5d7ObseICz0aOenx12Tvv/2+f1t99DL7vWtS9whlHwDBWIr+osQ1VIKxuOQ5LwIIIICAChCMcR0ggAACCCCAQBoEpv1ov2zrfs8x1O/82Uj5h7/80/7/tuiFbnl861HHNl4X67fVkGDM1soaNC6CMYOKQVcQQACBFAoQjKWw6AwZAQQQQACBFAroI5WZR/fKmf/8o2P0z88bJ7OvGiFuf37Nly6SPfdenUKtT4ZMMJbq8kczeIKxaJw5CwIIIICAuwDBGFcGAggggAACCKRFwG1GWF/45Taj7OXet1CWentlUu22Hn5Ppo0eUrL7BGMlidigXAGCsXIF2R8BBBBAoBwBgrFy9NgXAQQQQAABBJImkL+GmPb/6+OGy88PnHAMJf8xy6SN062/7/3hQ2ntXW/thTc+Huu+3rXTSjWCsVJC/HnZAgRjZRNyAAQQQACBMgQIxsrAY1cEEEAAAQQQSJzAhtePy+1rDxTt95jPny+/WTI5cWMr1GGdHaZhWPuud0TDMW1NUy+XlrrqkmMkGCtJxAblChCMlSvI/ggggAAC5QgQjJWjx74IIIAAAgggkESB7/7kN/Ljfz1WsOvPzh0r86ZcmsSh9fd54Oyw7NtnzhnLht611ep611Yr1QjGSgnx52ULEIyVTcgBEEAAAQTKECAYKwOPXRFAAAEEEEAgsQJfXPqvcux//q9z+q+PVb545/hEjkvDMA3B/vvLb4vOEuubHeY2mDcfnCJVwwaVHCfBWEkiNihXgGCsXEH2RwABBBAoR4BgrBw99kUAAQQQQACBJAoUe5xy8Gc/JdnetbeqR5yfmKGVmh2WP5DMZYM9rS+m+xGMJeYySG5HCcaSWzt6jgACCNggQDBmQxUZAwIIIIAAAgj4EfjTlbvk0H/8oeAuU6uHyNa/meDnkJFv2zc7bPmWt3KzxIrNDsvv3NLpV8iyGVWe+kww5omJjcoRIBgrR499EUAAAQTKFSAYK1eQ/RFAAAEEEEAgSQJu64uN/Nxnznms0tR1xnpOns0tov8/dv+76P8O0rruniDTRg/xtCvBmCcmNipHgGCsHD32RQABBBAoV4BgrFxB9kcAAQQQQACBpAhsP3JaZv7D63LmP//Y32WdHfbX135B7lh30DEMDcu235Mx4pHKvtlhzZ3duTDMz+yw/NoMOf/Tcmrl9Z5LRjDmmYoNgwoQjAWVYz8EEEAAgTAECMbCUOQYCCCAAAIIIJAEgUmPvyZ7/+33/V0duJ7Yrc+8IT8/cMIxjPqrL5F/+quvxj40XUi/9sn9ofRD30Spb6T02gjGvEqxXWABgrHAdOyIAAIIIBCCAMFYCIgcAgEEEEAAAQSMF1j0Qrc8vvWoo5/3TrtcHrutOvffuo//QTKP7nXMJtP//nxviDS7N0yKu4UVjrXVj5WGyZd6Hg7BmGcqNgwqQDAWVI79EEAAAQTCECAYC0ORYyCAAAIIIICAyQJuodeYz58vv1ky2dFtt/DMbbu4xqqL7OvMsXIepXzzwSlSNWyQ5yEQjHmmYsOgAgRjQeXYDwEEEEAgDAGCsTAUOQYCCCCAAAIImCww7Uf7ZVv3e44uFlpcP/9xS91p4MyyuMepoZiGYxqS+W2ZywbLvu9d42s3gjFfXGwcRIBgLIga+yCAAAIIhCVAMBaWJMdBAAEEEEAAARMFWrYdlXt7F60f2L4+bri8eOd41+5ueP243L72gOPPdC2yTd+9Sm4YdbERQwwajjVNvVxa6j5+dNRrIxjzKsV2gQUIxgLTsSMCCCCAQAgCBGMhIHIIBBBAAAEEEDBSwO0RyoEL7hfq9Lf+8VfS8dq7jj++5ksXyZ57rzZinEEfqdRF93XxfT+NYMyPFtsGEiAYC8TGTggggAACIQkQjIUEyWEQQAABBBBAwDgBt4Br2YwqWTr9iqJ9LbQQ/+O9s62ae2ddxdl6Tp7NPUqpP/20Ied/WnR9Mf3ppxGM+dFi20ACBGOB2NgJAQQQQCAkAYKxkCA5DAIIIIAAAggYJbB25ztyx7qDjj75mfXl9giml9lmlUQI+gil9mna6CHSdfcE390jGPNNxg5+BQjG/IqxPQIIIIBAmAIEY2FqciwEEEAAAQQQMEVg+Za3ZNnmnv7ujPzcZ+TJvxgjs308Spg/4yzIMcLy0FCsseOgdPaugRak6Sw5nS3ntxGM+RVje98CBGO+ydgBAQQQQCBEAYKxEDE5FAIIIIAAAgggUAEBr6GYzgrbetj59s2+7uhsMf1zv41gzK8Y2/sWIBjzTcYOCCCAAAIhChCMhYjJoRBAAAEEEEAAgQoI6Eyx9l3vFD1y5rLBuUclde2xiY/udWyr64qdWnl9oJ4RjAViYyc/AgRjfrTYFgEEEEAgbAGCsbBFOR4CCCCAAAIIIBCeQOu2o9Lc2e0pFOtbWD//rZX6Jkp9I2WQRjAWRI19fAkQjPniYmMEEEAAgZAFCMZCBuVwCCCAAAIIIIBASAJeQrGqYYNyM8X058A2cKF+DcU0HAvSCMaCqLGPLwGCMV9cbIwAAgggELIAwVjIoBwOAQQQQAABBBAIQUAfndSZYhpwFWo6Q0xDMX2M0q3pvrPXHpC2+rHnBGdeu0gw5lWK7QILEIwFpmNHBBBAAIEQBAjGQkDkEAgggAACCCCAQIgC+Y9Cuh1aQzENvILOBPPaXYIxr1JsF1iAYCwwHTsigAACCIQgQDAWAiKHQAABBBBAAAEEQhIwKRTTIRGMhVRYDlNYgGCMqwMBBBBAIE4BgrE49Tk3AggggAACCCDwiYC+UbL2yf25N0sWay111dI09fJI6AjGImFO90kIxtJdf0aPAAIIxC1AMBZ3BTg/AggggAACCCAgubXENBTTGWPFmgZiGoxF1QjGopJO8XkIxlJcfIaOAAIIGCBAMGZAEegCAggggAACCKRaQEOxxo6D0vn6caNCMe0MwViqL81oBk8wFo0zZ0EAAQQQcBcgGOPKQAABBBBAAAEE4hPwGoo1TL40N1NMF92PshGMRamd0nMRjKW08AwbAQQQMESAYMyQQtANBBBAAAEEEEilgM4Ua9/1TtGxZy4bLF13T4g8FNNOEYyl8rKMdtAEY9F6czYEEEAAAacAwRhXBAIIIIAAAgggEI9A67aj0tzZbWwoRjAWz3WRurMSjKWu5AwYAQQQMEqAYMyoctAZBBBAAAEEEEiJgJdQrGrYoNxMMf0ZV2PGWFzyKTovwViKis1QEUAAAQMFCMYMLApdQgABBBBAAAGrBfTRSZ0ppuuLFWq6lpiGYvoYZZyNYCxO/ZScm2AsJYVmmAgggIChAgRjhhaGbiGAAAIIIICAlQLZt89I7ZP7S4ZibfVjpe6qEbEbEIzFXgL7O0AwZn+NGSECCCBgsgDBmMnVoW8IIIAAAgggYJOAl1BMx6uhmL6F0oRGMGZCFSzvA8GY5QVmeAgggIDhAgRjhheI7iGAAAIIIICAFQI9J8/mZorpz2Ktpa5amqZebsyYCcaMKYW9HSEYs7e2jAwBBBBIggDBWBKqRB8RQAABBBBAIMkCupaYhmI6Y6xY00BMgzGTGsGYSdWwtC8EY5YWlmEhgAACCREgGEtIoegmAggggAACCCRSQEOxxo6D0vn68cSFYtphgrFEXnbJ6jTBWLLqRW8RQAAB2wQIxmyrKONBAAEEEEAAAVMEvIZiup6YzhTTN1Ga1gjGTKuIhf0hGLOwqAwJAQQQSJAAwViCikVXEUAAAQQQQCBRAs2d3dK67WjRPmcuGyxdd08wMhTTjhOMJeqSS2ZnCcaSWTd6jQACCNgiQDBmSyUZBwIIIIAAAgiYJKCBmAZjxZrpoRjBmElXlMV9IRizuLgMDQEEEEiAAMFYAopEFxFAAAEEEEAgUQJeQrGqYYNyM8X0p8mNGWMmV8eSvhGMWVJIhoEAAggkVIBgLKGFo9sIIIAAAgggYKRA+653cjPFdH2xQk3XEtswb5xMGz3EyDEM7BTBmPElSn4HCcaSX0NGgAACCCRZgGAsydWj7wgggAACCCBgkkD27TNS++T+kqFYW/1YqbtqhEldL9gXgrFElCnZnSQYS3b96D0CCCCQdAGCsaRXkP4jgAACCCCAgAkCXkIx7aeGYvoWyqQ0grGkVCrB/SQYS3Dx6DoCCCBggQDBmAVFZAgIIIAAAgggEKtAz8mzuZli+rNYa6mrlqapl8faV78nJxjzK8b2vgUIxnyTsQMCCCCAQIgCBGMhYnIoBBBAAAEEEEidgK4lpqGYzhgr1jQQ02AsaY1gLGkVS2B/CcYSWDS6jAACCFgkQDBmUTEZCgIIIIAAAghEKqChWGPHQel8/biVoZgOimAs0ksqnScjGEtn3Rk1AgggYIoAwZgplaAfCCCAAAIIIJA0AQ3F9C2UxZquJ6YzxfRNlElsBGNJrFrC+kwwlrCC0V0EEEDAMgGCMcsKynAQQAABBBBAIBKB5s5uad12tOi5MpcNlq67JyQ2FNPBEYxFcjml+yQEY+muP6NHAAEE4hYgGIu7ApwfAQQQQAABBJImoIGYBmPFmg2hGMFY0q7MhPaXYCyhhaPbCCCAgCUCBGOWFJJhIIAAAggggEAkAl5CMX1sct/3rpGqYYMi6VMlT8KMsUrqcuycAMEYFwICCCCAQJwCBGNx6nNuBBBAAAEEEEiSgK4npjPFdNH9Qk1DsQ3zxsm00UOSNLSCfSUYs6KMZg+CYMzs+tA7BBBAwHYBgjHbK8z4EEAAAQQQQCAMgezbZ6T2yf0lQ7G2+rFSd9WIME5pxDEIxowog92dIBizu76MDgEEEDBdgGDM9ArRPwQQQAABBBCIW8BLKKZ91FBM30JpUyMYs6maho6FYMzQwtAtBBBAICUCBGMpKTTDRAABBBBAAIFAAvrY5MRH90rPybNF92+pq5amqZcHOofJOxGMmVwdS/pGMGZJIRkGAgggkFABgrGEFo5uI4AAAggggEDFBTQUm732gGw9/F7Rc2kgpsGYjY1gzMaqGjYmgjHDCkJ3EEAAgZQJEIylrOAMFwEEEEAAAQQ8CWgo1thxUDpfP57aUEwHTjDm6XJho3IECMbK0WNfBBBAAIFyBQjGyhVkfwQQQAABBBCwUUBDMX0LZbGmi+zrumL6JkpbG8GYrZU1aFwEYwYVg64ggAACKRQgGEth0RkyAggggAACCBQVaO7sltZtR4tuk7lssHTdPcHqUEwBCMb4y1JxAYKxihNzAgQQQACBIgIEY1weCCCAAAIIIIDAJwIaiGkwVqxVDRsk+753jfWhGMEYfzMiESAYi4SZkyCAAAIIFBAgGOPSQAABBBBAAAEEPhbwEorpY5M6U0xnjKWhMWMsDVWOeYwEYzEXgNMjgAACKRcgGEv5BcDwEUAAAQQQQCAnoIvs67piuuh+oaahmK4ppmuLpaURjKWl0jGOk2AsRnxOjQACCCAgBGNcBAgggAACCCCQdoHs22ek9sn9hGIuFwLBWNr/dkQwfoKxCJA5BQIIIIBAQQGCMS4OBBBAAAEEEEizgJdQTH10pljD5EtTR0UwlrqSRz9ggrHozTkjAggggMAnAgRjXA0IIIAAAgggkFYBfWxSZ4ppOFasLZ1+hSybUZVKJoKxVJY92kETjEXrzdkQQAABBJwCBGNcEQgggAACCCCQRgENxWavPSBbD79XdPhNUy+XlrrqNBLlxkwwltrSRzdwgrHorDkTAggggMC5AgRjXBUIIIAAAgggkDYBDcV0oX1dcL9Y00cnNRTTRffT2gjG0lr5CMdNMBYhNqcKRWDdunVSX1/ff6zq6mr5l3/5Fxk1apTv4z/yyCOyePHi/v3mzJkjenwaAghEJ0AwFp01ZwpP4DvZ1fLMW7/oP+DgTw+SZWPnyqLRt3k6ycsnDsjcPY/JsbOn+rd/dFyD5/09nYSNEEAAAQSMFdBQrH3XO0X7l7lssHTdPSHVoZgCEYwZexnb0zGCMXtqmZaRPP3007JgwQLHcJcsWSIrVqzwTTBs2DA5deqTDyXTp0+XzZs3+z4OOyCAQHABgrHgduwZn8Bjh1+Q7x1od3Rg5KChsr1mpVRfWHph5Gu2LpLXTh9xBGsbr3tAbhw+Lr5BcWYEEEAAgUgElm3ukeVb3iIU86hNMOYRis2CCxCMBbdjz3gETpw4ISNGjHCcXGeNHT582FeHNm3aJLNmzXLss2bNGpk/f76v47AxAgiUJ0AwVp4fe8cnkB9uaU/mXlYjHZMWFe3U0l93yN8d/Iljm6bqW6Vl/Lz4BsOZEUAAAQQiEWjddlSaO7uLnqtq2KDcTDH9SWPGGNdABAIEYxEgc4rQBebOnSvr1693HDebzcqECRM8n2vhwoXy1FNPObY/fvy4DB8+3PMx2BABBMoXIBgr35AjxCOgj0POevVhOfPhWUcHfjp5sdw+8jrXTnW//45ktjY79hl94Ug5dPPqeAbBWRFAAAEEIhPwEorpWmIaiuljlLSPBZgxxpVQcQGCsYoTc4IKCLjN9rrrrrtk9WrvHyzyH6NkfbEKFIpDIuBBgGDMAxKbGCvQ/MZaae1+0dG/YkHX13uDtI3/vtdzkGbswOkYAggggIAvAV1kX9cV00X3CzUNxdrqx0rdVc6nY3ydyMKNCcYsLKppQyIYM60i9MerQH6wNXToUDl58qSn3d2CtY6ODtGZaDQEEIhWgGAsWm/OFr7AmJcWyuH3jzkO7PZo5LNvvSR3Zp90bHfnFbfIjzMLw+8UR0QAAQQQMEYg+/YZqX1yP6FYwIoQjAWEYzfvAgRj3q3Y0iyB+++/X1auXOno1MaNG2XmzJklO5r/GKWfUK3kwdkAAQR8CRCM+eJiYwMFnj/2qnxz1ypHz/QtldlpLf0L8esjlDfsWOJ4C6WfxfoNHDZdQgABBBDwINBz8mwuFNOfxZrOFGuYXPrlLR5Oad0mBGPWldS8ARGMmVcTeuRNYP/+/ZLJZBwbe3mcUhfvHzNmjONtlF7289YrtkIAAb8CBGN+xdjeRIHvZFfLM2/9wtE1fcPktpqHc//N7c+fydwtd1xxs4nDoU8IIIAAAiEI6GOTGorpjLFirWnq5dJSVx3CGe08BMGYnXU1alQEY0aVg874FLj22mtlz549jr1KLaC/bt06qa+vd+yzfft2qamp8Xl2NkcAgTAECMbCUOQYcQu4zQjTPmn4NfQzg8+ZUTbrC9fIz697IO5uc34EEEAAgQoJaCima4rp2mKEYuUhE4yV58feHgQIxjwgsYmxAk8//bQsWLDA0b9Sa4Xlv9GyurpaDh8+bOwY6RgCtgsQjNle4fSMz20NMX1ccuRnh8prp4/0Q+Q/ZpkeIUaKAAIIpEPAayimj07qTDFddJ9WWIBgjKuj4gIEYxUn5gQVFNDHIkeMcL61pdjbJd22X7JkiaxYsaKCveTQCCBQTIBgjOvDJoGpOx6Ql08cKDqkR8c1yKLRt9k0bMaCAAIIIDBAQGeKte96p6hJ5rLB0nX3BEIxD1cOwZgHJDYpT4BgrDw/9o5fIH8GmPaou7tbRo0adU7n3B6jLLRt/COjBwikQ4BgLB11Tsso9ZHKzNZmOfOh+yLLV188SvZOeywtHIwTAQQQSJ3Ass09snzLW4RiIVaeYCxETA7lLkAwxpWRdIFNmzbJrFmzHMNYs2aNzJ8//5yh5YdokyZNkt27dyedgP4jkGgBgrFEl4/Ouwgs/XWH/N3Bn7ja6GL8uig/DQEEEEDAPoHWbUelubO76MCqhg3KzRTTnzRvAgRj3pzYqgwBgrEy8NjVGIFhw4Y53jLpFni5PUZZKEAzZmB0BIEUCBCMpaDIKRtisccp77ziFvlxZmHKRBguAgggYL+APjqpoZiuL1ao6VpiGorpY5Q07wIEY96t2DKgAMFYQDh2M0rg/vvvl5UrVzr6lP+IpNtjlKXeYGnUIOkMApYKEIxZWtiUDuuxwy/I9w60Fx39TycvlttHXpdSIYaNAAII2CeQffuM1D65v2Qo1lY/Vuqucq6PbJ9G+CMiGAvflCPmCRCMcUnYILB//37JZDKOoeQvqp//GGWxRfptMGEMCCRFgGAsKZWin6UE3NYX07dSHjt7yrEr64yVkuTPEUAAgeQIEIpVvlYEY5U3Tv0ZCMZSfwlYA3DttdfKnj17+sdTXV0thw8fzv3/I0eOiP7/ga2jo0M0LKMhgEC8AgRj8fpz9vAE6vc8Juve3uE4oM4OW3Hw/5HXTh9x/Pem6lulZfy88E7OkRBAAAEEIhfoOXk2N1NMfxZrOlOsYfKlkffPlhMSjNlSSYPHQTBmcHHomi+Bp59+WhYsWODYJ5vNyoQJEyT/z4YOHSonT570dXw2RgCByggQjFXGlaNGK/DsWy/JndknHSftW0/s5RMHRNcdG9gGf3qQbLzuARbij7ZMnA0BBBAITUDXEtNQTGeMFWtNUy+XljrnF/ShdSIlByIYS0mh4xwmwVic+pw7TAG3xfXvuusuWb16teTPJuv772Gen2MhgEAwAYKxYG7sZY6APkJ5w44ljkcm9RHK7TUrpfrCj2cIfCe7Wp556xeOTuvbKfUtlTQEEEAAgWQJaCjW2HFQOl8/TigWQekIxiJATvspCMbSfgXYNf78dcR0Zpg+Xpn/GOX27dulpqbGrsEzGgQSKkAwltDC0e1+AbfQ69FxDbJo9G3927iFZ/qH+dvBigACCCBgtoDXUEwfndSZYvomSlp5AgRj5fmxtwcBgjEPSGySGIFNmzbJrFmzHP2dPn26bNmypf+/DVx7LDEDo6MIWCxAMGZxcVMwtOePvSrf3LXKMdJCi+u7PW6ZP7MsBWQMEQEEEEi0gM4Ua9/1TtExZC4bLF13TyAUC6nSBGMhQXKYwgIEY1wdtgkMGzZMTp1yvgFs4Bjz31Zp2/gZDwJJEyAYS1rF6O9AgTEvLZTD7x9zoOjjkfqYpFvTtcZ0zbGBbe5lNdIxaRGwCCCAAAKGC7RuOyrNnd2EYhHXiWAsYvA0no5gLI1Vt3vM999/v6xcubLgILu7u2XUqFF2IzA6BBIkQDCWoGLRVYeA2yOUfQvuF6LSUGzWqw/LmQ+dbzB7JnO33HHFzQgjgAACCBgq4CUUqxo2KDdTTH/SwhMgGAvPkiMVECAY49KwTWD//v2SyWRchzVp0iTZvXu3bUNmPAgkWoBgLNHlS23n3QIufSzyd9PXljRxC9RGXzhSDt28uuS+bIAAAgggEL2APjqpM8V0fbFCTdcS01BMH6OkhStAMBauJ0dzESAY47KwUSD/LZR9Y1yzZo3Mnz/fxiEzJgQSK0AwltjSpbrjnHMTBQAAIABJREFU12xdJK+dPuIw8DPry+0RzFKzzVINzuARQACBmASyb5+R2if3lwzF2urHSt1VI2Lqpd2nJRizu75GjI5gzIgy0ImQBdwW4ddTHD9+XIYPHx7y2TgcAgiUI0AwVo4e+8Yl8MUt8+TY2U/Ws/S7Tpgu2v83/+8/lHWMuMbOeRFAAIG0CBCKmVFpgjEz6mB1LwjGrC5vagd34sQJGTNmjGMR/jlz5si6detSa8LAETBVgGDM1MrQLwQQQAABBNIr0HPybG6mmP4s1lrqqqVp6uXphYpg5ARjESCn/RQEY2m/AuwcvwZg9fX1jsF1dHTI3Llz7Rwwo0IgwQIEYwkuHl1HAAEEEEDAQgFdS0xDMZ0xVqxpIKbBGK2yAgRjlfXl6L0CBGNcBjYK5K8xNnToUDl58qSNQ2VMCCRegGAs8SVkAAgggAACCFgjoKFYY8dB6Xz9OKGYIVUlGDOkEDZ3g2DM5uqmc2xub6W86667ZPVq3vaVziuCUZsuQDBmeoXoHwIIIIAAAukQ8BqKNUy+NDdTTN9ESau8AMFY5Y1TfwaCsdRfAtYBLFy4UJ566inHuDZu3CgzZ860bqwMCAEbBAjGbKgiY0AAAQQQQCD5As2d3dK67WjRgWQuGyxdd08gFIuw3ARjEWKn9VQEY2mtvJ3jdlt0v7q6Wg4fPmzngBkVAhYIEIxZUESGgAACCCCAQMIFNBDTYKxYIxSLp8gEY/G4p+qsBGOpKrf1g3VbdH/JkiWyYsUK68fOABFIqgDBWFIrR78RQAABBBCwQ8BLKFY1bFBuppj+pEUrQDAWrXcqz0YwlsqyWzvoGTNmyJYtWxzjy2azMmHCBGvHzMAQSLoAwVjSK0j/EUAAAQQQSK5A+653cjPFdH2xQk3XEtNQTGeM0aIXIBiL3jx1ZyQYS13JrR3wkSNHRB+bHNgmTZoku3fvtnbMDAwBGwQIxmyoImNAAAEEEEAgeQLZt89I7ZP7S4ZibfVjpe6qEckboCU9JhizpJAmD4NgzOTq0Dc/Avo2ytraWjl16lRuNw3FWlpapKamxs9h2BYBBCIWIBiLGJzTIYAAAggggIB4CcWUSUMxfQslLT4BgrH47FNzZoKx1JSagSKAAAJGChCMGVkWOoUAAggggIC1Aj0nz+ZmiunPYq2lrlqapl5urUNSBkYwlpRKJbifBGMJLh5dRwABBCwQIBizoIgMAQEEEEAAgYQI6FpiGorpjLFiTQMxDcZo8QsQjMVfA+t7QDBmfYkZIAIIIGC0AMGY0eWhcwgggAACCFgjoKFYY8dB6Xz9OKFYgqpKMJagYiW1qwRjSa0c/UYAAQTsECAYs6OOjAIBBBBAAAGTBbyGYrqemM4U0zdR0swQIBgzow5W94JgzOryMjgEEEDAeAGCMeNLRAcRQAABBBBIvEBzZ7e0bjtadByZywZL190TCMUMqzbBmGEFsbE7BGM2VpUxIYAAAskRIBhLTq3oKQIIIIAAAkkU0EBMg7FijVDM3MoSjJlbG2t6RjBmTSkZCAIIIJBIAYKxRJaNTiOAAAIIIJAIAS+hmD42ue9710jVsEGJGFPaOkkwlraKxzBegrEY0DllIgRmz54tLS0tUlVVlYj+0kkEkipAMJbUytHvSgvM3vVDaRk/T6ouuKTSp+L4CRXQNZP0Q79+qNc36NEQQMAp0L7rndxMMf27Uqjp358N88bJtNFD4DNUgGDM0MLY1C2CMZuqyVjCFOj7u1FXV0dAFiYsx0IgT4BgjEsCAXeBvr8bdSOnEJBxkTgE+gKx//7y27kP/EunXyHLZlShhAACAwSyb5+R2if3lwzF2urHSt1VI7AzWIBgzODi2NI1gjFbKsk4whbI/7tBQBa2MMdD4GMBgjGuBASKB2N9f0pAxpWSH4j1iRCMcW0g4BTwEorpHhqK6VsoaWYLEIyZXR8rekcwZkUZGUQFBAr93SAgqwA2h0y1AMFYqsvP4IsIFPq7QUCWvsumUCBGMJa+a4ERlxboOXk2N1NMfxZrLXXVPIJcmtOILQjGjCiD3Z0gGLO7vowuuECxvxt6VAKy4LbsicBAAYIxrgcE3AWK/d3I3Yd4xNL6S6dUIEYwZv0lwAB9CujfmdlrD8jWw+8V3VPX5NNgjJYMAYKxZNQp0b0kGEt0+eh8BQVKBWN9pyYgq2AROHQqBAjGUlFmBhlAoFQw1n8fIiALoGv2Ll4DMYIxs+tI76IV0L83jR0HpfP144Ri0dJX/GwEYxUn5gQEY1wDCLgLeA3GCMi4ghAoT4BgrDw/9rZXwGswRkBmzzXgNxAjGLOn9oykfAENxfQtlMWaLrKv64rpmyhpyREgGEtOrRLbU78f/hM7UDqOQEQCzCCLCJrTWCPg98O/NQNnIAhUSIBHLCsEW8HDBg3EKtglDo2AdQKZywZL190TCMUSWFmCsQQWLUldfu+992To0KFJ6jJ9RSAxAgRkiSkVHY1ZgGAs5gJwemsFCMjMLy2BmPk1ood2CBCKJbuOBGPJrp/xvW9sbJT29nbj+0kHEUiyAAFZkqtH36MQIBiLQplzpFmAgMy86hOImVcTemSvgD42qTPFNByjJVOAYCyZdUtEr5ctWybLly/v7+ugQYNk3Lhxsnfv3kT0n04ikDSBlpYWaWpqSlq36S8CFRVY9ut1svzg+v5zfPq8T8mHH/2xoufk4AikVaBl/Dxpqr41rcM3Ztwail350E7RnzQEEKiswKA/+T9kTubzUjVskCybUVXZk3H0igkQjFWMNt0Hbm1tlebm5nQjMHoEIhKYNm2aLF26VPQnDQEEPhFo7X5Rmt9YCwkCCFRYYNqI8bJ07BzRnzRzBFq3HZXmzm5zOkRPELBc4KOWqZaP0N7hEYzZW9vYRkYoFhs9J06ZAIFYygrOcH0JEIr54mJjBAIJEIgFYot8JwKyyMk5YUoFCMaSW3iCseTWzsieZ7NZqa2tFV10n4YAApURIBCrjCtHtUcge/pNqX3lQXnvf79vz6AYCQIGCRCIGVQMH10hIPOBxaYIBBAgGAuAZsguBGOGFMKGbngJxYYMGSJtbW2ii4XTEEi7wHnnneeLgEDMFxcbp1TASyg25E8ulLaJ94guGE5DIM0Cfl9MQSBmx9USJCBbOv0K1k+yo/yMokyB85q3FTwCwViZuDHuTjAWI75Np+7p6cnNFNOfxZqGYg0NDTYNnbEgEFjAazBGIBaYmB1TJtDzwbu5mWL6s+i9qDcUa/jyTSnTYbgInCvgNRgjELPz6vETkBGM2XkNMCr/AgRj/s2SsAfBWBKqZHgf9bFJDcV0xlixpm/L07fm0RBA4GOBUsEYgRhXCgLeBfSxSQ3FdMZY0XtR7xvz9M15NAQQ6L0PvTC7KAOBWDquEi8BGcFYOq4FRllagGCstFEStyAYS2LVDOqzhmKNjY3S2dlJKGZQXehKMgQKBWMEYsmoH700R0BDscZ9T0jnsZ2EYuaUhZ4kQKBQMEYgloDiVaCLxQIygrEKgHPIRAoQjCWybCU7TTBWkogNCgl4DcX00UmdKabri9EQQOATgfxgjECMqwMB/wJeQzF9dFJniun6YjQEEPhYID8YIxDjylABt4CMYIxrA4H/+r3JGmNWXgoEY1aWNZpB6Uyx9vb2oifLZDLS1dVFKBZNSThLwgT6gjECsYQVju4aJaAzxdp/+8vi96KLr5Su6x8iFDOqcnTGBIG+YIxAzIRqmNeHgQEZwZh59aFH8QgwYywe90qflWCs0sKWHr+1tVWam5sJxSytL8OKRkDX5lu6dKloMEZDAAH/Aq3dL0rzG2sJxfzTsQcCOQFdl2/p2DmiwRgNgUICGpC994cPeSsllwgCvQIEY3ZeBgRjdta1oqPyEopVVVXlZorpTxoCCCCAAAJhC3gJxaouuCQ3U0x/0hBAAAEEEEAAgXIFCMbKFTRzf4IxM+tibK/00UmdKabrixVqupaYhmL6GCUNAQQQQACBsAX00UmdKabrixW8F/WuJaahWKb3MUoaAggggAACCCAQhgDBWBiK5h2DYMy8mhjbo2w2K/roV6lQrK2tTerq6owdBx1DAAEEEEiuQPb0m7nHv0qFYm0T75G6kVOSO1B6jgACCCCAAALGCRCMGVeSUDpEMBYKo/0HIRSzv8aMEAEEEDBdgFDM9ArRPwQQQAABBOwWIBizs74EY3bWNdRR9fT05GaK6c9iTWeKNTQ0hHpuDoYAAggggIAK9Hzwbm6mmP4sei/qnSnW8OWbQEMAAQQQQAABBEIXIBgLndSIAxKMGVEGczuhj01qKKYzxoq1pqYmaWlpMXcg9AwBBBBAILEC+tikhmI6Y6zovaj6VmkZPy+x46TjCCCAAAIIIGC2AMGY2fUJ2juCsaByKdhPQ7HGxkbp7OwkFEtBvRkiAgggYKKAhmKN+56QzmM7CcVMLBB9QgABBBBAIEUCBGN2FptgzM66lj0qr6GYPjqpM8X0TZQ0BBBAAAEEwhTwGorpo5M6U2xI75soaQgggAACCCCAQKUECMYqJRvvcQnG4vU39uw6U6y9vb1o/zKZjHR1dRGKGVtFOoYAAggkW0BnirX/9pfF70UXXyld1z9EKJbsUtN7BBBAAAEEEiFAMJaIMvnuJMGYbzL7d2htbZXm5mZCMftLzQgRQAABYwVau1+U5jfWEooZWyE6hgACCCCAQPoECMbsrDnBmJ11DTwqL6FYVVVVbqaY/qQhgAACCCAQtoCXUKzqgktyM8X0Jw0BBBBAAAEEEIhCgGAsCuXoz0EwFr25sWfURyd1ppiuL1ao6VpiGorpY5Q0BBBAAAEEwhbQRyd1ppiuL1bwXtS7lpiGYpnexyhpCCCAAAIIIIBAVAIEY1FJR3segrFovY09Wzabldra2pKhWFtbm9TV1Rk7DjqGAAIIIJBcgezpN6X2lQdLhmJtE++RupFTkjtQeo4AAggggAACiRQgGEtk2Up2mmCsJJH9GxCK2V9jRogAAgiYLkAoZnqF6B8CCCCAAAIIEIzZeQ0QjNlZV8+j6unpyc0U05/FWktLizQ1NXk+LhsigAACCCDgVaDng3dzM8X0Z9F70fh50lR9q9fDsh0CCCCAAAIIIBCqAMFYqJzGHIxgzJhSRN8RXUtMQzGdMVasaSCmwRgNAQQQQACBsAV0LTENxXTGWNF7UW8g1tIbjNEQQAABBBBAAIG4BAjG4pKv7HkJxirra+zRNRRrbGyUzs5OQjFjq0THEEAAAbsFNBRr3PeEdB7bSShmd6kZHQIIIIAAAlYIEIxZUcZzBkEwZmddi47KayjW0NCQmymmb6KkIYAAAgggEKaA11Cs4cs35WaKDel9EyUNAQQQQAABBBCIU4BgLE79yp2bYKxytsYeubm5WVpbW4v2L5PJSFdXF6GYsVWkYwgggECyBZrfWCut3S8WvxddfKV0Xf8QoViyS03vEUAAAQQQsEaAYMyaUjoGQjBmZ10LjkoDMQ3GijVCsZRdFAwXAQQQiFhAAzENxoreiwjFIq4Kp0MAAQQQQACBUgIEY6WEkvnnBGPJrFugXnsJxaqqqnIzxfQnDQEEEEAAgbAFvIRiVRdckpsppj9pCCCAAAIIIICAKQIEY6ZUItx+EIyF62ns0drb23MzxXR9sUJN1xLTUExnjNEQQAABBBAIW6D9t7/MzRTT9cUK3ot61xLTUCzTO2OMhgACCCCAAAIImCRAMGZSNcLrC8FYeJbGHimbzUptbW3JUKytrU3q6uqMHQcdQwABBBBIrkD29JtS+8qDJUOxton3SN3IKckdKD1HAAEEEEAAAWsFCMbsLC3BmJ117R+Vl1BMN9ZQTN9CSUMAAQQQQCBsAS+hWO5e1BuK6VsoaQgggAACCCCAgIkCBGMmVqX8PhGMlW9o7BF6enpyM8X0Z7HW0tIiTU1Nxo6DjiGAAAIIJFeg54N3czPF9GfRe9H4edJUfWtyB0rPEUAAAQQQQMB6AYIxO0tMMGZnXXOPTWoopjPGijUNxDQYoyGAAAIIIBC2gK4lpqGYzhgrei/qDcRaeoMxGgIIIIAAAgggYLIAwZjJ1QneN+uDsUceeUQWL17cLzR06FB57rnnZObMmZ7Ujhw5In/+538u3d3d/dtv3LjR8/6eThLyRhqKNTY2SmdnJ6FYyLYcDgEEEAgi0P3+O5LZ2ixnPjwbZHcZfeFIOXTz6kD7xrWThmKN+56QzmM7CcXiKgLnRQABBAwRWL7lLVm2uae/NyM/9xlZ/3/9n3LDqIt99/C7P/mN/Phfj/XvN7V6iGz9mwm+j8MOCAQRIBgLomb+PtYHY5s2bZJZs2Y5KlFdXS07d+6U4cOHl6zQjBkzZMuWLf3babC2Z88eGTVqVMl949jAayim64npTDF9EyUNAQQQQKCyAs++9ZLcmX2yrJP8dPJiuX3kdWUdI6qdvYZiup6YzhQb0vsmShoCCCCAgL0Ci17olse3HnUMsP7qS+Sf/uqrvgd90fd3yJn//GP/ftd86SLZc+/Vvo/DDggEESAYC6Jm/j7WB2NagvxwS//bkiVLZMWKFUUrtG7dOqmvr3dss2bNGpk/f76xlW1ubpbW1tai/ctkMtLV1UUoZmwV6RgCCNgm8PyxV+Wbu1aVNawkBWPNb6yV1u4Xi9+LLr5Suq5/iFCsrKuCnRFAAIFkCHQf/4OMXrHL0VmdNfa75X/mawAt247KvZ2fPMmjO9877XJ57LZqX8dhYwSCChCMBZUze79UBGP6OOSkSZPk1KlTjmro+lsTJrhPuz1x4oSMGTPGsY8eY/fu3cZWVAMxDcaKNUIxY8tHxxBAwHKBqTse6H2U8g+eRvmb9393zmOXH922wdO+cW+kgZgGY0XvRYRicZeJ8yOAAAKRC0z70X7Z1v2e47zPzxsns68a4bkvtz7zhvz8wAnH9ofvnyzVI873fAw2RKAcAYKxcvTM3TcVwZjyP/3007JgwQJHJYoFXQsXLpSnnnrKc5AWd4m9hGL62OS+ffukqqoq7u5yfgQQQACBAgJus8tmfeEa+fl1Dxhv5iUU08cm9017XKouuMT48dBBBBBAAIHwBNxme3193HB58c7xnk+S/xgl64t5pmPDkAQIxkKCNOwwqQnG1P3aa6/NrQ82sLk9Grljxw654YYbHNutWrVK7rvvPsPK93F32tvbczPFdH2xQk1DsQ0bNsi0adOMHAOdQgABBBD4WGDMSwvl8PufLCo8+NODJDutRaovvNRoovbf/jI3U0zXFyt4L+oNxTZM/r5MG+H9Q5DRg6ZzCCCAAAK+BPKDrcGf/ZT8/oc1no7hFqwtm1ElS6df4Wl/NkIgDAGCsTAUzTtGqoKx/fv3iz5KOLDpYvqHDh3qX4hfH6GcMmWK4y2Ufhbrj7rE+jhobW1tyVCsra1N6urqou4e50MAAQQQ8CHwnexqeeatXzj2eHRcgywafZuPo0S/afb0m1L7yoMlQ7G2ifdI3cgp0XeQMyKAAAIIGCHwrX/8lXS89q6jL4/XVUvz1MtL9i//MUo/oVrJg7MBAh4FCMY8QiVss1QFY1qbRx55RBYvXuwo05w5c0QX2i/059u3b5eaGm/fZERZfy+hmPZHQzF9CyUNAQQQQMBcgZdPHBBdh2xgu/riUbJ32mPmdrq3Z15Csdy9qDcU07dQ0hBAAAEE0iuw4fXjcvvaAw4AL49T6uL9mUf3Ot5G6WW/9Eoz8koJEIxVSjbe46YuGHObEaYl0PDroosuOmdG2V133SWrV6+Ot0ouZ+/p6cnNFNOfxVpLS4s0NTUZ1386hAACCCDgFLhm6yJ57fSR/v+oj1Bu7F1X7Mbh44yl6vng3dxMMf1Z9F40fp40Vd9q7DjoGAIIIIBAdAJ/unKXHPoP58toSi2gv3zLW7Jss/Nzz7Nzx8q8KWYvMxCdKmeKSoBgLCrpaM+TumBMed3WENPHJUePHi1btmzpr0D+Y5bRlqbw2XQtsdmzZ8vWrVuLdkkDMQ3GaAgggAACZgvo2ly6cP3AducVt8iPMwuN7biuJTZ71w9l6/E3it+LegOxlt5gjIYAAggggIAKLHqhWx7fetSBUWqtsPw3Wo783Gfkd8v/DFAEIhcgGIucPJITpjIYU9m5c+fK+vXriyJv3LhRZs6cGUkhvJ5EQ7HGxkbp7OwkFPOKxnYIIICAwQLd778jma3NcubDs/29HH3hSDl0s3mzlfs6qKFY474npPPYTkIxg68tuoYAAgiYKKCPRY5escvRtWJvl3Tbvv7qS+Sf/uqrJg6PPlkuQDBmZ4FTG4zpI5VjxoyRU6dOuVZ2+vTpsnnzZuOqrqGYvoWyWNP1xHSmmL6JkoYAAgggYLaArium64sNbD+dvFhuH3mdsR3XUEzfQln0XtS7npjOFBvS+yZKGgIIIIAAAgMF8meA6Z+9fE9Gbhh18TlQbo9RFtoWZQQqLUAwVmnheI6f2mBMuXXB/fr6elf57u5uGTVqVDxVKXDW5uZmaW1tLdonfetmV1cXoZhRlaMzCCCAgLvAY4dfkO8daHf84dzLaqRj0iJjydwe+8zvbObiK6Xr+ocIxYytIh1DAAEE4hVo2XZU7u3sdnTi3mmXy2O3VZ/TsfwQbcznz5ffLJkc7wA4e2oFCMbsLH2qg7Fij1OuWrVK7rvvPmOqroGYBmPFGqGYMeWiIwgggEBJAX2E8oYdS+TY2U9mLo8cNFS216yU6gvNXExY10HTYKzovYhQrGTt2QABBBBAQOSi7+9wvGXSLfBye4yyUICGKQJRCBCMRaEc/TlSG4xt2rRJZs2aVVQ8m83KhAkToq9K3hm9hGL62OS+ffukqqoq9v7SAQQQQACB0gL1ex6TdW/vcGz4TOZuueOKm0vvHMMWXkIxfWxy37THpeqCS2LoIadEAAEEEEiSwLf+8VfS8Zrzrcb5j0i6PUZZ6g2WSTKgr8kTIBhLXs289DiVwZjb+mL6Vkp9fHJgM2GdMV1PTGeK6aL7hZqGYhs2bJBp06Z5qTnbIIAAAgjELPDsWy/JndknHb24cfg42VbzcMw9cz+9riemM8V00f2C96LeUGzD5O/LtBHjjRwDnUIAAQQQMEtgw+vH5fa1zjU28xfVz3+Mstgi/WaNjt7YKkAwZmdlUxmM3X///bJy5UpHRXV22OLFi2XLli2O/75mzRqZP39+LNXXPtXW1pYMxdra2qSuri6WPnJSBBBAAAH/Al/cMs/xCOXgTw+S7LQWIx+hzJ5+U2pfebBkKNY28R6pGznFPwZ7IIAAAgikVuBPV+6SQ//xh/7xj/zcZ+R3y/8s9/+3HzktNz6Rddgsm1ElS6dfkVovBh6/AMFY/DWoRA9SF4zt2LFDbrjhBodl33piR44cEZ05NrANHTpU9uzZE/lC/F5CMe2nhmL6FkoaAggggEAyBL6TXS3PvPULR2f/duxfyvKvuL8MJs5ReQnFcvei3lCsofctlDQEEEAAAQT8CCx6oVse33rUscvz88bJ7KtGSP6fDf7sp+T3P6zxc3i2RSB0AYKx0EmNOGCqgjF9hHLKlCmORyY1CNu5c6cMHz48V5BHHnkkN3NsYJszZ07uDZZRNX1scuLEidLT01P0lC0tLdLU1BRVtzgPAggggECZAs8fe1W+uWuV4yhXXzxK9k57rMwjh7+7PjY5ceu90vOBc/2X/DO1jJ8nTdW3ht8BjogAAgggYL2A2+L6Xx83XF68c7zkzybr++/WozBAowUIxowuT+DOpSoYcwu9Nm7cKDNnzuwHdAvP9A/ztwssXmJHDcVmz54tW7duLbqlBmIajNEQQAABBJIjcM3WRfLa6SP9HdZHKDde94Do+mImNQ3FZu/6oWw9/kbxe1FvIKbBGA0BBBBAAIGgAvnriOnMsE3fveqcxyifnTtW5k0x863NQcfOfskTIBhLXs289Dg1wdj+/fslk8k4TAotru/2uGX+zDIvuH630VCssbFROjs7CcX84rE9AgggYLiALl6vb3Yc2O684hb5cWahUT3XUKxx3xPSeWwnoZhRlaEzCCCAgJ0CLduOyr2dzpegXfOli2Tvv/2+f8AD1x6zU4FRJUWAYCwplfLXz9QEY9dee21urbCBTd9COWrUKFexuXPnyvr16x1/tmTJElmxYoU/YR9bayimb6Es1nSRfV1XTN9ESUMAAQQQSIZA9/vvSGZrs5z58Gx/h0dfOFIO3bzauAFoKKZvoSx6L+pdZF/XFRvS+yZKGgIIIIAAAuUKXPT9HXLmP/9Y8DD5b6ss93zsj0BQAYKxoHJm75eKYMztEcq+BfcLlUcX4p80aZKcOnXKscn27dulpib8RR+bm5ultbW16NWiM966uroIxcz+O0XvEEAAgXME3GaLjRw0VEZ+dqhnra9fOqniC/S79TO/g5mLr5Su6x8iFPNcOTZEAAEEECgl8K1//JV0vFZ4TcuX78nIDaMuLnUY/hyBigsQjFWcOJYTWB+MuQVc+ljk4cOHS4K7BWoalu3evbvkvn420EBMg7FiraqqSvbt20co5geWbRFAAAFDBNzeROm3a5VepF8f89RgrOi96IJLZN+0xwnF/BaP7RFAAAEEigpseP243L72gOs2Yz5/vvxmyWQEETBCgGDMiDKE3gnrg7EZM2bIli1bHHB+Zn25PYJZaraZnyp5CcX0sUmdKZa/Rpqf87AtAggggEB8AqYHY15CMX1sUmeK6YwxGgIIIIAAAmEL5L+Fsu/49067XB67rTrs03E8BAIJEIwFYjN+J+uDsdGjR4uuJdbX/K4Tpov2f/Ob3yzrGIWuAl1kX9cV00X3CzUNxXRNMV1bjIYAAgggkEyBZ996SZreeNaxxpjfkcy9rEY6Ji3yu1vJ7XWRfV1XTBfdL3gv6g3FdE2xut61xWgIIIAAAghUQsBtEX49z+H7J0v1iPMrcUqOiYBvAYIx32SJ2MH6YMzUKmSzWamtrSUUM7XaqQ1IAAAfjElEQVRA9AsBBBBIgUD29JtS+8qDhGIpqDVDRAABBEwX6D7+B8k8utexCP/U6iGy9W8mmN51+pciAYIxO4tNMBZDXXt6emTixIlFQzHtls4Ua2hoiKGHnBIBBBBAwHaBng/elYlb7y0aiuXuRb0zxRq+fJPtHIwPAQQQQCBmgeVb3pJlm3scvVg2o0qWTr8i5p5xegQ+ESAYs/NqIBiLuK762KTOFNMZY8Xa0qVLZdmyZRH3jtMhgAACCKRBQB+b1JliOmOs6L1o7BxZ9pW5aSBhjAgggAACMQvkrzE2+LOfkt//sCbmXnF6BJwCBGN2XhEEYxHWVUMxXVNM1xYr1pqamqSlpSXCnnEqBBBAAIG0CGgopmuK6dpiRe9F1bdKy/h5aWFhnAgggAACMQq4vZXy6+OGy4t3jo+xV5wagXMFCMbsvCoIxiKqq59QTGeL6aL7NAQQQAABBMIU8BOKLe2dLaZvoqQhgAACCCBQaYFbn3lDfn7ghOM0j9dVS/PUyyt9ao6PgC8BgjFfXInZmGAsolLpTLH29vaiZ8tkMtLV1UUoFlFNOA0CCCCQNgGdKdb+218WvxddfKV0Xf8QoVjaLg7GiwACCMQk4Lbo/sjPfUZ+t/zPYuoRp0WgsADBmJ1XB8FYBHXVtcKWL19OKBaBNadAAAEEEHAXWPbrdbL84HpCMS4QBBBAAAGjBNwW3a+/+hL5p7/6qlH9pDMIqADBmJ3XAcFYheva2toqzc3NRc9SVVWVmymmP2kIIIAAAgiELdDa/aI0v7G2+L3ogktyM8Wqen/SEEAAAQQQiEpg0uOvyd5/+73jdM/PGyezrxoRVRc4DwKeBQjGPFMlakOCsQqWS0MxnSmm64sVarqWmIZi+hglDQEEEEAAgbAFNBTTmWK6vljBe1HvWmIaimV6H6OkIYAAAgggEJXA9iOn5cYnso7Tjfn8+fKbJZOj6gLnQcCXAMGYL67EbEwwVqFSZbNZqa2tLRmKtbW1SV1dXYV6wWERQAABBNIskD39ptS+8mDJUKxt4j1SN3JKmqkYOwIIIIBADAL6Nsq//r9/LWf+84+5s2so9v3/9mWZN+XSGHrDKREoLUAwVtooiVsQjFWgaoRiFUDlkAgggAACvgQIxXxxsTECCCCAAAIIIFBSgGCsJFEiNyAYC7lsPT09uZli+rNY05liDQ0NIZ+dwyGAAAIIICDS88G7uZli+rPovah3pljDl2+CDAEEEEAAAQQQQMCDAMGYB6QEbkIwFmLRdC0xDcV0xlix1tTUJC0tLSGemUMhgAACCCDwsYCuJaahmM4YK3ovqr5VWsbPgw0BBBBAAAEEEEDAowDBmEeohG1GMBZSwTQUa2xslM7OTkKxkEw5DAIIIICAPwENxRr3PSGdx3YSivmjY2sEEEAAAQQQQKCkAMFYSaJEbkAwFkLZvIZi+uikzhTTN1HSEEAAAQQQCFPAayimj07qTLEhvW+ipCGAAAIIIIAAAgh4FyAY826VpC0JxkKols4Ua29vL3qkTCYjXV1dhGIheHMIBBBAAIFzBXSmWPtvf1n8XnTxldJ1/UOEYlxACCCAAAIIIIBAAAGCsQBoCdiFYKzMIrW2tkpzczOhWJmO7I4AAgggEFygtftFaX5jLaFYcEL2RAABBBBAAAEESgoQjJUkSuQGBGNllM1LKFZVVZWbKaY/aQgggAACCIQt4CUUq7rgktxMMf1JQwABBBBAAAEEEAgmQDAWzM30vQjGAlZIH53UmWK6vlihpmuJaSimj1HSEEAAAQQQCFtAH53UmWK6vljBe1HvWmIaimV6H6OkIYAAAggggAACCAQXIBgLbmfyngRjAaqTzWaltra2ZCjW1tYmdXV1Ac7ALggggAACCBQXyJ5+U2pfebBkKNY28R6pGzkFTgQQQAABBBBAAIEyBQjGygQ0dHeCMZ+FIRTzCcbmCCCAAAKhCxCKhU7KARFAAAEEEEAAgZICBGMliRK5AcGYj7L19PTkZorpz2JNZ4o1NDT4ODKbIoAAAggg4E2g54N3czPF9GfRe1HvTLGGL9/k7aBshQACCCCAAAIIIFBSgGCsJFEiNyAY81g2XUtMQzGdMVasNTU1SUtLi8ejshkCCCCAAALeBXQtMQ3FdMZY0XtR9a3SMn6e9wOzJQIIIIAAAggggEBJAYKxkkSJ3IBgzEPZNBRrbGyUzs5OQjEPXmyCAAIIIBC+gIZijfuekM5jOwnFwufliAgggAACCCCAQEkBgrGSRIncgGCsRNm8hmL66KTOFNM3UdIQQAABBBAIU8BrKKaPTupMsSG9b6KkIYAAAggggAACCIQrQDAWrqcpRyMYK1EJnSnW3t5edKtMJiNdXV2EYqZc1fQDAQQQsExAZ4q1//aXxe9FF18pXdc/RChmWe0ZDgIIIIAAAgiYI0AwZk4twuwJwVgRzdbWVmlubiYUC/OK41gIIIAAAr4EWrtflOY31hKK+VJjYwQQQAABBBBAIHwBgrHwTU04IsFYgSp4CcWqqqpyM8X0Jw0BBBBAAIGwBbyEYlUXXJKbKaY/aQgggAACCCCAAAKVEyAYq5xtnEcmGHPR10cndaaYri9WqOlaYhqK6WOUNAQQQAABBMIW0EcndaaYri9W8F7Uu5aYhmKZ3scoaQgggAACCCCAAAKVFSAYq6xvXEcnGMuTz2azUltbWzIUa2trk7q6urjqxnkRQAABBCwWyJ5+U2pfebBkKNY28R6pGznFYgmGhgACCCCAAAIImCNAMGZOLcLsCcHYAE1CsTAvLY6FAAIIIBBEgFAsiBr7IIAAAggggAAClRcgGKu8cRxnIBj7L/Wenp7cTDH9Way1tLRIU1NTHLXinAgggAAClgv0fPBubqaY/ix6Lxo/T5qqb7Vcg+EhgAACCCCAAAJmCRCMmVWPsHpDMNYrqWuJaSimM8aKNQ3ENBijIYAAAgggELaAriWmoZjOGCt6L+oNxFp6gzEaAggggAACCCCAQLQCBGPRekd1ttQHYxqKNTY2SmdnJ6FYVFcd50EAAQQQcAhoKNa47wnpPLaTUIxrAwEEEEAAAQQQMFSAYMzQwpTZrVQHY15DsYaGhtxMMX0TJQ0BBBBAAIEwBbyGYg1fvik3U2xI75soaQgggAACCCCAAALRCxCMRW8exRmtDMY08PISYjU3N0tra2tR50wmI11dXZ6OF0XBOAcCCCCAQDIENPDyEmI1v7FWWrtfLH4vuvhK6br+IU/HS4YOvUQAAQQQQAABBJInQDCWvJp56bGVwZiGXRqOLVu2rKCBbqPBWLFGKOblEmIbBBBAAAE3AQ27NBxb9pW5he9FvdtoMFb0XkQoxgWGAAIIIIAAAggYIUAwZkQZQu+ElcHY7Nmzc2uG1dXVyYYNG85B8xKKVVVV5WaK6U8aAggggAACfgVm7/phbs2wupFTZMPk7597L/IQilVdcEluppj+pCGAAAIIIIAAAgjEK0AwFq9/pc5uXTCmM8WuvPLK3IwxbdOmTcsFXH2tvb09N1Os78/dYPUxTN1HZ4zREEAAAQQQ8CugM8Wu/MX83Iyx3L1oxPhcwNV/L/rtL3Mzxfr+3PVe1LuWmO6T6Z0xRkMAAQQQQAABBBCIX4BgLP4aVKIH1gVjOlNMZ4wNbH2PRPb09EhtbW3JUKytrS0324yGAAIIIIBAEAGdKaYzxhz3ov96JLLng3el9pUHS4ZibRPvyc02oyGAAAIIIIAAAgiYIUAwZkYdwu6FdcGYriu2fPnyc5z0kUidJVZqphihWNiXGMdDAAEE0iew7NfrZPnB9efei3ofidRZYqVmihGKpe+aYcQIIIAAAgggYL4AwZj5NQrSQ+uCMX2MUmeGBWktLS3S1NQUZFf2QQABBBBAoF9AH6PUmWFBWsv4edJUfWuQXdkHAQQQQAABBBBAoIICBGMVxI3x0FYFY9lsViZOnBiIUwMxDcZoCCCAAAIIlCOQPf2mTNx6b6BDaCCmwRgNAQQQQAABBBBAwDwBgjHzahJGj6wKxnRh/cbGRt8uhGK+ydgBAQQQQKCAQHvvwvqN+57w7UMo5puMHRBAAAEEEEAAgUgFCMYi5Y7sZFYFY7rovi6+76c1NDTkZorpmyhpCCCAAAIIlCugi+7r4vt+WsOXb8rNFBvS+yZKGgIIIIAAAggggICZAgRjZtal3F5ZFYz9/+3dQYhd1RkH8AMWbYoaW6MY0ToaVCQJnSSSIMVm0o0uLIkubC1IQqxZtaQTXQVkEoSuJBmoK6WzKVjrQiuCtKtERCTFyEiVNmmnpigoNpSxRlO10uY+M0NGE+e9d8+7890zv7eZRd499zu//wsDf86708/zxZwWq/sRcj0BAgQInCnQz/PFnBbzGSJAgAABAgQIxBdQjMXPqJ8JiynG6jxfbGRkJB04cKAfP9cQIECAAIFZgTrPFxtZtiod+O7DNAkQIECAAAECBIIKKMaCBlNzrGKKsfHx8TQ6Oto3h3KsbzoXEiBAgMBpgfGp59Lo6xN9eyjH+qZzIQECBAgQIEBg4AKKsYETL8gNiinG+nm+2BfFh4eHOyfHPG9sQT6LbkqAAIHWC/TzfLEv/S5aem3n5JjnjbX+42ADBAgQIECAQGECirHCAj29nSKKsenp6VQ9X6z6Wfc1NDTUKceqn14ECBAgQKBbgelPP0zV88Wqn3VfQ9+4vFOOVT+9CBAgQIAAAQIEYggoxmLkkHuKIoqx6i9RVifG6ryqU2I7d+5MW7ZsSdXJMS8CBAgQINCLQPWXKKsTY3Ve1SmxndfdkbYs35CGT50c8yJAgAABAgQIEIgjoBiLk0XOSYooxvbs2ZP27t3bs0tVhlUl2NatWzuFmK9Q9kzoAgIECBA4LbDnL0+mvUd+27NHVYZVJdjWqzd1CjFfoeyZ0AUECBAgQIAAgUYEFGONMDd+kyKKseprlMeOHesaz+mwrqm8kQABAgS6FKi+Rnnso/e6fHfqFGBOh3XN5Y0ECBAgQIAAgQUXUIwteAQDGaD1xdjk5GRas2bNvDgzp8PGxsY6p8ScDpuXzBsIECBAoEuByfffTGsO7pr33TOnw8Zu/GHnlJjTYfOSeQMBAgQIECBAIIyAYixMFFkHaX0xNj4+nkZHR8+JUj1Ev/qq5LZt2zxQP+tHx2IECBAgMCMwPvVcGn194ty/i049RL/6quS2b3/fA/V9bAgQIECAAAECLRVQjLU0uHnGbn0xVj10v3r4/pmvmdNh+/fv75RhToeV+eG1KwIECEQRqB66Xz18f87votPPDtu/anunDHM6LEpa5iBAgAABAgQI9CegGOvPLfpVrS/Gzny+2MjISNq8eXPndJgyLPpHz3wECBAoR+DM54uNLFuVNl+xvnM6TBlWTsZ2QoAAAQIECBBQjJX5GWh1MVY9X2zTpk2pKsRmnh1WZkx2RYAAAQJRBarni2166aFUFWIzzw6LOqu5CBAgQIAAAQIE+hdQjPVvF/nKVhdj09PTHVunwyJ/xMxGgACBsgWmP/3w899Fp7466UWAAAECBAgQIFCugGKszGxbXYyVGYldESBAgAABAgQIECBAgAABAtEEFGPREskzj2Isj6NVCBAgQIAAAQIECBAgQIAAgYIFFGNlhqsYKzNXuyJAgAABAgQIECBAgAABAgQyCijGMmIGWkoxFigMoxAgQIAAAQIECBAgQIAAAQIxBRRjMXOpO5VirK6g6wkQIECAAAECBAgQIECAAIHiBRRjZUasGCszV7siQIAAAQIECBAgQIAAAQIEMgooxjJiBlpKMRYoDKMQIECAAAECBAgQIECAAAECMQUUYzFzqTuVYqyuoOsJECBAgAABAgQIECBAgACB4gUUY2VGrBgrM1e7IkCAAAECBAgQIECAAAECBDIKKMYyYgZaSjEWKAyjECBAgAABAgQIECBAgAABAjEFFGMxc6k7lWKsrqDrCRAgQIAAAQIECBAgQIAAgeIFFGNlRqwYKzNXuyJAgAABAgQIECBAgAABAgQyCijGMmIGWkoxFigMoxAgQIAAAQIECBAgQIAAAQIxBRRjMXOpO5VirK6g6wkQIECAAAECBAgQIECAAIHiBRRjZUasGCszV7siQIAAAQIECBAgQIAAAQIEMgooxjJiBlpKMRYoDKMQIECAAAECBAgQIECAAAECMQUUYzFzqTuVYqyuoOsJECBAgAABAgQIECBAgACB4gUUY2VGrBgrM1e7IkCAAAECBAgQIECAAAECBDIKKMYyYgZaSjEWKAyjECBAgAABAgQIECBAgAABAjEFFGMxc6k7lWKsrqDrCRAgQIAAAQIECBAgQIAAgeIFFGNlRqwYKzNXuyJAgAABAgQIECBAgAABAgQyCijGMmIGWkoxFigMoxAgQIAAAQIECBAgQIAAAQIxBRRjMXOpO5VirK6g6wkQIECAAAECBAgQIECAAIHiBRRjZUasGCszV7siQIAAAQIECBAgQIAAAQIEMgooxjJiBlpKMRYoDKMQIECAAAECBAgQIECAAAECMQUUYzFzqTuVYqyuoOsJECBAgAABAgQIECBAgACB4gUUY2VGrBgrM1e7IkCAAAECBAgQIECAAAECBDIKKMYyYgZaSjEWKAyjECBAgAABAgQIECBAgAABAjEFFGMxc6k7lWKsrqDrCRAgQIAAAQIECBAgQIAAgeIFFGNlRqwYKzNXuyJAgAABAgQIECBAgAABAgQyCijGMmIGWkoxFigMoxAgQIAAAQIECBAgQIAAAQIxBRRjMXOpO5VirK6g6wkQIECAAAECBAgQIECAAIHiBRRjZUasGCszV7siQIAAAQIECBAgQIAAAQIEMgooxjJiBlpKMRYoDKMQIECAAAECBAgQIECAAAECMQUUYzFzqTuVYqyuoOsJECBAgAABAgQIECBAgACB4gUUY2VGrBgrM1e7IkCAAAECBAgQIECAAAECBDIKKMYyYgZaSjEWKAyjECBAgAABAgQIECBAgAABAjEFFGMxc6k7lWKsrqDrCRAgQIAAAQIECBAgQIAAgeIFFGNlRqwYKzNXuyJAgAABAgQIECBAgAABAgQyCijGMmIGWkoxFigMoxAgQIAAAQIECBAgQIAAAQIxBRRjMXOpO5VirK6g6wkQIECAAAECBAgQIECAAIHiBRRjZUasGCszV7siQIAAAQIECBAgQIAAAQIEMgooxjJiBlpKMRYoDKMQIECAAAECBAgQIECAAAECMQUUYzFzqTuVYqyuoOsJECBAgAABAgQIECBAgACB4gUUY2VGrBgrM1e7IkCAAAECBAgQIECAAAECBDIKKMYyYgZaSjEWKAyjECBAgAABAgQIECBAgAABAjEFFGMxc6k7lWKsrqDrCRAgQIAAAQIECBAgQIAAgeIFFGNlRqwYKzNXuyJAgAABAgQIECBAgAABAgQyCijGMmIGWkoxFigMoxAgQIAAAQIECBAgQIAAAQIxBdpcjE0dP5mGHzmcTnz8WV+411+2JB3dvb6va6NfpBiLnpD5CBAgQIAAAQIECBAgQIAAgQUXaHMxNnHo3XTfk0dqGT69fWW6c/WyWmtEvFgxFjEVMxEgQIAAAQIECBAgQIAAAQKhBNpcjD3zp+Pprok3ankqxmrxuZgAAQIECBAgQIAAAQIECBAg0F6BNhdjlfrIo6+lE59091XKI+999KWvXf5v/8b2hvcVkzsxVmSsNkWAAAECBAgQIECAAAECBAjkFGh7MdatxdlOl92x8tL03E9WdbtEq96nGGtVXIYlQIAAAQIECBAgQIAAAQIEFkJgsRRjN/zij+mv/zw5S3zhBeelyQfXpRXLliwE+8DvqRgbOLEbECBAgAABAgQIECBAgAABAm0XWAzF2I6njqbHX35nTlT7tqxIoxuvant855xfMVZstDZGgAABAgQIECBAgAABAgQI5BIovRh78e/vp+/9cnIO17qrL0qv7FqbizDkOoqxkLEYigABAgQIECBAgAABAgQIEIgkUHoxdvO+V9Phtz6YJa++Qvn8jtXp1uuWRooh+yyKseykFiRAgAABAgQIECBAgAABAgRKEyi5GHvg2am07+DbcyK7/5bl6bG7bygtxi/tRzFWfMQ2SIAAAQIECBAgQIAAAQIECNQVKLUYmzp+Mg0/cjid+PizWaLrL1uSju5eX5esFdcrxloRkyEJECBAgAABAgQIECBAgACBhRQotRgbefS19MLU9Bzap7evTHeuXraQ3I3dWzHWGLUbESBAgAABAgQIECBAgAABAm0VKLEY2//C22nX76bmRHLP2svTE/fe1NaYep5bMdYzmQsIECBAgAABAgQIECBAgACBxSZQWjFWfYXy1lN/hfKdf38yG+Xyi89PL/5sOK1YtmTRxKsYWzRR2ygBAgQIECBAgAABAgQIECDQr0BpxdiPf/3n9JtX35vD8asf3Zi2b7iiX6JWXqcYa2VshiZAgAABAgQIECBAgAABAgSaFCipGJs49G6678kjc/g2rrgkHfzpd5okDXEvxViIGAxBgAABAgQIECBAgAABAgQIRBYoqRi7cuzlOV+hvPCC89Lkg+sW1VcoZz5rirHI/+vMRoAAAQIECBAgQIAAAQIECIQQKKUY2/HU0fT4y+/MMd1z+1Aau+2aEM5ND6EYa1rc/QgQIECAAAECBAgQIECAAIHWCZRQjD3zp+Pprok35tivu/qi9Mquta3LI9fAirFcktYhQIAAAQIECBAgQIAAAQIEihUooRi7ed+r6fBbH8xmVH2F8vkdq9Ot1y0tNrf5NqYYm0/IvxMgQIAAAQIECBAgQIAAAQKLXqDtxdgDz06lfQffnpPj/bcsT4/dfcOizlYxtqjjt3kCBAgQIECAAAECBAgQIECgG4E2F2NTx0+m4UcOpxMffza71esvW5KO7l7fzdaLfo9irOh4bY4AAQIECBAgQIAAAQIECBDIIdDmYuxsp8WWX3x+unLpBV3T/GDlpUU+oF8x1vVHwBsJECBAgAABAgQIECBAgACBxSrQ5mLsbH+JstccS31Iv2Ks10+C9xMgQIAAAQIECBAgQIAAAQKLTkAxVuZfr1SMLbr/yjZMgAABAgQIECBAgAABAgQI9CrQ5mJs4tC7aeczf5vzjLFe93/P2svTE/fe1Otl4d+vGAsfkQEJECBAgAABAgQIECBAgACBhRZoczG20HaR768Yi5yO2QgQIECAAAECBAgQIECAAIEQAoqxEDFkH0Ixlp3UggQIECBAgAABAgQIECBAgEBpAoqx0hL9fD+KsTJztSsCBAgQIECAAAECBAgQIEAgo4BiLCNmoKUUY4HCMAoBAgQIECBAgAABAgQIECAQU0AxFjOXulMpxuoKup4AAQIECBAgQIAAAQIECBAoXkAxVmbEirEyc7UrAgQIECBAgAABAgQIECBAoEeBPb8/ds4r9v7hH+f8t7Hbrjnrv12y5Gvp5xuv6nEKb29SQDHWpLZ7ESBAgAABAgQIECBAgAABAmEFqmLsqwqwXgevCrM9tw/1epn3NyigGGsQ260IECBAgAABAgQIECBAgACBuALTJ/+bvrn7pSwDVqfF3nxoQ6p+esUVUIzFzcZkBAgQIECAAAECBAgQIECAQMMCuU6NOS3WcHB93k4x1iecywgQIECAAAECBAgQIECAAIHyBHKcGnNarD2fC8VYe7IyKQECBAgQIECAAAECBAgQINCAQN1TY06LNRBSplsoxjJBWoYAAQIECBAgQIAAAQIECBAoQ6DOqTGnxdr1GVCMtSsv0xIgQIAAAQIECBAgQIAAAQINCPR7asxpsQbCyXgLxVhGTEsRIECAAAECBAgQIECAAAECZQj0c2rMabH2Za8Ya19mJiZAgAABAgQIECBAgAABAgQaEOj11JjTYg2EkvkWirHMoJYjQIAAAQIECBAgQIAAAQIEyhDo5dSY02LtzFwx1s7cTE2AAAECBAgQIECAAAECBAg0INDtqTGnxRoIYwC3UIwNANWSBAgQIECAAAECBAgQIECAQBkC3Zwac1qsvVkrxtqbnckJECBAgAABAgQIECBAgACBBgTmOzXmtFgDIQzoFoqxAcFalgABAgQIECBAgAABAgQIEChD4KtOjTkt1u6MFWPtzs/0BAgQIECAAAECBAgQIECAQAMC5zo15rRYA/gDvIVibIC4liZAgAABAgQIECBAgAABAgTKEDjbqTGnxdqfrWKs/RnaAQECBAgQIECAAAECBAgQINCAwBdPjTkt1gD6gG+hGBswsOUJECBAgAABAgQIECBAgACBMgTOPDXmtFgZmSrGysjRLggQIECAAAECBAgQIECAAIEGBGZOjTkt1gB2A7dQjDWA7BYECBAgQIAAAQIECBAgQIBAGQLVqbFrHz6U3nxoQ6pOjXm1W0Ax1u78TE+AAAECBAgQIECAAAECBAg0LHDsX/9JQ9/6esN3dbtBCCjGBqFqTQIECBAgQIAAAQIECBAgQIAAgfACirHwERmQAAECBAgQIECAAAECBAgQIEBgEAKKsUGoWpMAAQIECBAgQIAAAQIECBAgQCC8gGIsfEQGJECAAAECBAgQIECAAAECBAgQGISAYmwQqtYkQIAAAQIECBAgQIAAAQIECBAIL6AYCx+RAQkQIECAAAECBAgQIECAAAECBAYhoBgbhKo1CRAgQIAAAQIECBAgQIAAAQIEwgsoxsJHZEACBAgQIECAAAECBAgQIECAAIFBCCjGBqFqTQIECBAgQIAAAQIECBAgQIAAgfACirHwERmQAAECBAgQIECAAAECBAgQIEBgEAL/Bw2IfC6tDqq+AAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "21f9d344-82a5-4301-b7c0-0e831e49e08d",
   "metadata": {},
   "source": [
    "## Option 1: Frame of reference\n",
    "The frame of reference for the Ground reaction force and the kinematics should be the same. By default, the code assumes that the first two dimensions (x,y) are horizontal, and the third dimension (z) is vertical, oriented upwards, as in (1). \n",
    "![image.png](attachment:afed88e6-49cc-4d93-9895-86132a5b8436.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "504549d0-e7bd-4d70-8f8d-d9016ba7a524",
   "metadata": {},
   "source": [
    "If this is not the case, then the direction of gravity must be specified as gravity_direction: for example [0,-1,0] in (2), or [0,0,1] in (3). The code then calculates the CoM acceleration by subtracting the person's weight from the Ground reaction force."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bf58d44d-5a56-4f7b-a7cb-34db6b5ff725",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example in (2): Rotation of the Ground reaction force and kinematic CoM\n",
    "rotated_GroundReactionForce = numpy.zeros((numpy.shape(GroundReactionForce)))\n",
    "rotated_Kinematic_com       = numpy.zeros((numpy.shape(Kinematic_com)))\n",
    "rotation_indices            = [1,2,0]\n",
    "\n",
    "for i, index in enumerate(rotation_indices):\n",
    "    rotated_GroundReactionForce[i] = GroundReactionForce[index]\n",
    "    rotated_Kinematic_com[i]       = Kinematic_com [index]\n",
    "\n",
    "# The row corresponding to the vertical direction is specified when calculating the estimate of CoM position and velocity\n",
    "gravity_direction = [0,-1,0]\n",
    "rotated_args      = (rotated_GroundReactionForce, Force_frequency, rotated_Kinematic_com, Kinematic_frequency, mass)\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*rotated_args,gravity_direction = gravity_direction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b835db98-9459-4d1f-b938-804d671d37a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example in (3): Inversion of the Ground reaction force and kinematic CoM\n",
    "inverted_GroundReactionForce = - GroundReactionForce\n",
    "inverted_Kinematic_com       = - Kinematic_com\n",
    "# The direction of the vertical dimension is specified to -1 (downwards)\n",
    "gravity_direction = [0,0,1]\n",
    "inverted_args = (inverted_GroundReactionForce, Force_frequency, inverted_Kinematic_com, Kinematic_frequency, mass)\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*inverted_args,gravity_direction = gravity_direction)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc637a94-748a-433a-9104-8508a67e6bdc",
   "metadata": {},
   "source": [
    "## Option 2: Sub-sampling frequency\n",
    "In the example data, the kinematics are sampled at 200 Hz and the force at 1000 Hz. The code automatically selects a sub-sampling frequency of 200 Hz for the two signals (corresponding to the greatest common divisor). However, the user can specify any common divisor as the sub-sampling frequency through the optional argument \"sub_frequency\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ee1b8fff-545f-4e25-9131-279b0032721f",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_frequency = 100 # in Hertz\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,sub_frequency = sub_frequency)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "619f5b8c-053b-4b25-8743-5c97177997a1",
   "metadata": {},
   "source": [
    "## Option 3: Initial conditions\n",
    "By default, the initial position and speed for the estimator are obtained by a least-squares fit of the first and last ten samples of the kinematic Center of Mass. This number of samples can be modified by specifying the optional argument \"initial_samples\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1b4f96ba-5d93-4ae7-805c-eea06e6fd9a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_samples = 20\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,initial_samples = initial_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bb9eea5-b8dc-49f1-9228-8f826fbba08e",
   "metadata": {},
   "source": [
    "Alternately, the user may wish to use a custom method to determine the initial conditions. In this case, they can specify the desired initial conditions for both the forwards estimator using the optional argument \"Initial_conditions\", and for the backwards estimator using the optional argument \"Final_conditions\".\n",
    "\n",
    "In the example data, the subject is initially static. We can specify an initial velocity of zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ab2ab83b-7377-4dcb-be25-305b89cbcad5",
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_position = Kinematic_com[:,0]\n",
    "initial_velocity = numpy.zeros((3))\n",
    "Initial_conditions = numpy.array([initial_position,initial_velocity]).T\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,Initial_conditions = Initial_conditions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5fa6daac-0a67-40fd-95d6-38c6d00c66d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,Initial_conditions = Initial_conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1730a0c-2ba5-4263-b91f-9b7bd3022ccb",
   "metadata": {},
   "source": [
    "## Option 4: Measurement noise\n",
    "\n",
    "Different values of the position and measurement noise can be specified. These are used for the calculation of the estimator gains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ce9d2993-716a-4c2c-afb1-5bfa8431e0f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "Position_std = 0.0035 # standard deviation of the error in CoM position obtained from the kinematics alone (in m)\n",
    "Force_std    = 2      # standard deviation of the error in Ground reaction force (in N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "79228bb9-12ce-4bac-9855-636e6ec82561",
   "metadata": {},
   "outputs": [],
   "source": [
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,Position_std = Position_std, Force_std = Force_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b2b5e8c-8a23-4ce8-b6ef-b7ef97b47457",
   "metadata": {},
   "source": [
    "The measurement noise can be specified in each direction for either position, force or both:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "86497d60-61d1-4cbd-8b3c-8d07a54946a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "Position_std = 0.0035 # standard deviation of the error in CoM position obtained from the kinematics alone (in m)\n",
    "Force_std    = numpy.array([2,2,4]) # standard deviation of the error in Ground reaction force in each direction (in N)\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,Position_std = Position_std, Force_std = Force_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c0fd7a52-8d43-4652-9301-6dbf323e2d7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Position_std = numpy.array([0.002,0.002,0.004]) # standard deviation of the error in CoM position in each direction (in m)\n",
    "Force_std    = 2                                # standard deviation of the error in Ground reaction force(in N)\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,Position_std = Position_std, Force_std = Force_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "a7178a27-869c-4354-bb0e-252b80c16f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "Position_std = numpy.array([0.002,0.002,0.004]) # standard deviation of the error in CoM position in each direction (in m)\n",
    "Force_std    = numpy.array([2,2,4])             # standard deviation of the error in Ground reaction force in each direction (in N)\n",
    "Pos_estimate, Vel_estimate, Frequency = center_of_mass.estimator.optimal_combination(*args,Position_std = Position_std, Force_std = Force_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0cec802-c7a0-434e-ac10-a9727999bfb7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}