Tutorial_Estimating_Energy_Expenditure_During_Exercise.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Tutorial: Estimating Energy Expenditure During Exercise",
      "provenance": [],
      "collapsed_sections": [
        "PBzVrpDXG2s5"
      ],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/stanfordnmbl/mobilize-tutorials/blob/main/Tutorial_Estimating_Energy_Expenditure_During_Exercise.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X38L6tanrnrB"
      },
      "source": [
        "### Mobilize Center & Restore Center @ Stanford Tutorial\n",
        "# Estimating Energy Expenditure During Exercise \n",
        "\n",
        "Physical inactivity is the fourth leading cause of global mortality. Health organizations have requested a tool to objectively measure physical activity. Respirometry and doubly labeled water accurately estimate energy expenditure, but are infeasible for everyday use. Smartwatches are portable, but have significant errors. Existing wearable methods poorly estimate time-varying activity, which comprises 40% of daily steps. We developed a Wearable System that estimates metabolic energy expenditure in real-time during common steady-state and time-varying activities with substantially lower error than state-of-the-art methods. \n",
        "\n",
        "The Wearable System uses inertial measurement units worn on the shank and thigh as they distinguish lower-limb activity better than wrist or trunk kinematics and converge more quickly than physiological signals. When evaluated with a diverse group of new subjects, the Wearable System has a cumulative error of 13% across common activities, significantly less than 42% for a smartwatch and 44% for an activity-specific smartwatch. This approach enables accurate physical activity monitoring which could enable new energy balance systems for weight management or large-scale activity monitoring. We will walk through how we process data from the IMUs and perform estimation of energy expenditure.\n",
        "\n",
        "##Tutorial Overview\n",
        "**In this notebook we illustrate how to analyze human motion measured by two inertial measurement units worn on one leg to estimate the calories burned during exercise.** The notebook provides a pipeline that you can adapt to get started with IMU-based energy estimation projects.\n",
        "\n",
        "This work is designed to capture energy expenditure during cyclic motions, for example, walking, running, stair climbing, and biking. These cyclic motions allow for understanding periodic and repetitive movements associated with energy expenditure. For tasks with one-off movements like weight lifting, this approach is not well-suited. Here we focus on lower-limb movement, but this approach could capture upper-limb movement with additional sensors worn on the appropriate limb. We also are only able to estimate aerobic activities because we rely on respirometry data as a ground truth signal, which is inaccurate for measuring for anaerobic activites such as sprinting or weight lifting.\n",
        "\n",
        "The notebook is a part of the [Mobilize Center](https://mobilize.stanford.edu) webinar series, and is jointly offered with the [Restore Center](https://restore.stanford.edu). The Mobilize Center is an NIH-funded Biomedical Technology Resource Center which provides tools and training to help researchers produce insights from wearables, video, medical images, and other data sources. The Restore Center is an NIH Medical Rehabilitation Research Resource Network Center which is creating a worldwide collaboration to advance the use of real-world data in rehabilitation outcomes for those with movement impairments.\n",
        "\n",
        "##Background and Citation\n",
        "Most of the [code](https://github.com/pslade2/EnergyExpenditure) in this notebook comes from our study on estimating energy expenditure during exercise. To cite our work, please use:\n",
        "\n",
        "*Slade, P., Kochenderfer, M.J., Delp, S.L., Collins, S.. Sensing leg movement enhances wearable monitoring of energy expenditure. Nat Commun 12, 4312 (2021).*\n",
        "\n",
        "Learn more about this work:\n",
        "* [Read our publication](https://www.nature.com/articles/s41467-021-24173-x)\n",
        "* [See our hardware documentation](https://spectrum.ieee.org/a-diy-calorie-counter-more-accurate-than-a-smartphone)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PBzVrpDXG2s5"
      },
      "source": [
        "## Learning Goals\n",
        "\n",
        "In this tutorial notebook, you will learn how to:\n",
        "\n",
        "* Perform the data collection protocol for this approach\n",
        "  - Select sensors\n",
        "  - Place sensors\n",
        "* Interpret inertial measurement unit data for lower-extremity movement\n",
        "* Preprocess movement data for the energy expenditure algorithm\n",
        "  - Segment cyclic movement by gait cycle\n",
        "  - Filter lower-limb movement data\n",
        "  - Format data for each gait cycle\n",
        "* Estimate energy expenditure\n",
        "  - Use the model to estimate energy expenditure\n",
        "  - Compare estimates against ground truth\n",
        "  - Develop intuition into the model \n",
        "  - Explore time-varying estimation\n",
        " "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xM_Vf8l4f5iE"
      },
      "source": [
        "# Setup\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GsH_4sdoQPgB"
      },
      "source": [
        "\n",
        "## Introduction to Google Colab\n",
        "\n",
        "Google Colab is a cloud-based environment for running Python code interactively (via Jupyter notebooks, for those who are familiar with those). If you are new to Colab, you can learn about the key features in [this tutorial](https://colab.research.google.com/notebooks/basic_features_overview.ipynb). For the purposes of our tutorial, you only need to understand how to interact with the \"code cells.\" "
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "## Installation\n",
        "Run the code below to install and load the necessary packages and the data. The data for this tutorial are stored in a [Github repository](https://github.com/stanfordnmbl/mobilize-tutorials-data). It is done when you see `Unpacking objects: 100% (19/19), done.`"
      ],
      "metadata": {
        "id": "DQ8KtFKspj10"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Installing and loading necessary packages and data files\n",
        "from scipy import signal\n",
        "import numpy as np\n",
        "import time\n",
        "import os\n",
        "import sys\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "from google.colab import drive\n",
        "!git clone https://github.com/stanfordnmbl/mobilize-tutorials-data.git\n",
        "drive_path = '/content/mobilize-tutorials-data/' # path to folder containing code and data"
      ],
      "metadata": {
        "id": "QqLcHEddgJlE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Protocol for Collecting Data"
      ],
      "metadata": {
        "id": "R84zIq7dgjmk"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Sensors and Hardware\n",
        "\n",
        "The wearable system we used to estimate energy expenditure consisted of two inertial measurement unit (IMU) sensors, a microcontroller (Raspberry Pi), and rechargeable battery. The details of the components we selected are available in [this article](https://spectrum.ieee.org/a-diy-calorie-counter-more-accurate-than-a-smartphone). \n",
        "\n",
        "The components we have selected can be replaced with others, but may require modifications to the code we provide. In order to provide sufficient sensing to capture dynamic movement like running, you should select inertial measurement units that can measure +/-16G of acceleration, +/- 2000 degrees per second of rotation, and that can sample data at a rate of at least 50 Hz.\n",
        "\n",
        "## Placement of Sensors\n",
        "\n",
        "The IMUs were placed on the lateral side of the left leg, one in the middle of the thigh, and the other on the middle of the shank. The inertial measurement units consist of a three-axis accelerometer (accel) and a three-axis gyroscope (gyro). The IMUs were oriented so that the X, Y, and Z directions represent the fore-aft, mediolateral, and vertical axes, respectively. These orientations are assigned during quiet standing, but the sensors are attached to different body parts that move relative to the global reference frame. \n",
        "\n",
        "It is important to note that the accuracy of the calorie estimates are sensitive to the placement of the IMUs. For example, misalignment of the three axes changed estimates by a large amount, so the training data for the data-driven model was augmented with random small rotations to provide a model that was more robust to small errors in placement. Translation of the sensor along the vertical axis resulted in only very small changes in error. \n",
        "\n",
        "For use with other patient populations, we recommend following the axes alignment as closely as possible. Making distal or proximal translations to sensor placement should have little impact on the device performance. "
      ],
      "metadata": {
        "id": "FuEnMjXMkSne"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Interpreting Inertial Measurement Unit Data for Lower-Extremity Movement\n",
        "\n"
      ],
      "metadata": {
        "id": "w4EW9wzigqlE"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Loading and Formatting Time Series IMU Data\n",
        "\n",
        "We will now load pre-recorded IMU data into the variable `raw_IMU_data`. This data is the raw time series data recorded at 100 Hz during walking. We recommend a minimum sampling rate of 50 Hz to capture dynamic activities like running. The data should be stored in the following format:"
      ],
      "metadata": {
        "id": "pHFF-akq8-wv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "![cap_format.PNG](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABCkAAAN0CAYAAACZQXLAAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAEnQAABJ0Ad5mH3gAAIrzSURBVHhe7f0/q3RZGz949VvoN2BgpxP2gIiBgdCTzYhiJ4KpjfgGOvqBE3QgA8OM9KiImDQa22g2zjRMpraJgtKJghp0MozCGP1avs/Dfn7rXs+3qvbeVXVWnXM+wYf7PnWu/W/tVXtd69q76nz1L//lv/wLAAAAYDVFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEA/M1PP/3011dfffUP+bnFAcAz/Pbbb1+MQ999912N42NTpAD45H7//fe/vvnmGxPTBX755Ze/vv7667+1e/799ddfa9xbUaQAuO3PP//82+R5u1Z+//33f3utxXKMIgWhSAHwyY0Fik2ShBbL4/zxxx//1O4pVKxMdBUpAG774YcfvrhWul4+jiIFoUgBcEImmN9+++3DBtGs78cff/zizkxk0prXfv7557/FtGXvNW5vI9l6vjkR26wsEClSANw2j9VhMv0YihSEIgXAAbnLPU/k7hlEs752R+aSZ0waW7KVj4C0WP7Z2H5Hzk/O/fZRj01+brFvRZECWOHsdXSV+VoZuZnQYjlGkYJQpAC4IZPJfHdAPnM6DpyjM4NoBuJ5krpHihptfWelILEliPnoR461xdHdk1znOyi2j9vkyZzVxSFFCmCF91akGG8wZBzPk5AtjuMUKQhFCoAbxuTpkqODaAoBbR2ZtCb52eIyWLcnLSREr2MsNL33Sb0iBbDCR7qOch9FCkKRAuCGRxcp5sf88/9bj4m2osazvqOCY8ZzokgBcJzrDhtFCkKRAuCGRxcpIo/1pzgRex/xz9MT4zYf/bEPzhnPiSIFwHGuO2wUKQhFCoATHjGIZh1HvoNg/pOV+S6DFsfbyVMx4zlRpAA45qNdR7mPIgWhSAFwwqpBdPyzp9FiVkvb5CmP7QshI//PF4/md22ZV5OkOYnyeAw5x/NHbOZ+oEhxTvti2rT3rS9xTftnufHjU1ku+z1+twvvz3Zut/Oa9+Krvb8+c7/d/mz2OCblePYc/+yjXUdfUfpVPlaa8zO29Sv2u7k/ZB9bHB+bIgXACasG0TnBeMRfgzh6LEloWmySnPkjKc09H1PJ8c6JceTnrDfH0pZrLh1HtjEWJ0YtMd9r3PZm3Ic4kpxnX5J0ZrIz789WFDo6Wbhnf6651NbpM3OfnuXY5gQ6P8+Tw9mRj1LNtrad9+1su46yT2mPrHvuZ2f68ShfvJvl5/5wZL339IG5vW5t71LspevInmvHW1wjPkO/TTuN6xh/l22MBZamHf/okdfRox7xPolx2Wgxl2z9dDxXadOjfXQ07su4jvz/0pi2uaffpT23Y5n7RV7L79I327LN3PeyjhbXPOrcsp4iBcAJ9wyi9xgTmnjEgHv0WMbEPQN/XksCMicF1xyZ+ETWPx/7JYnbkxC1CUiWu5R85/XE7N2P2bjtzbgPsaddcr5uTXRGOS97E8Qz+7PHuN6tzyQhvtTWs7FPZuKzt69l/dcmSrPEJqFu65ol6T/y/ssE8ch7JAn13n1PW+5d9639vqcPzO+NW+3TYq+1fybH8zo2z75GfLZ+m5hxuWwj0i/H168Zj3+291zN2rr2euT7JOZlWswsbXjk+n1JW/f4+0zW81reM+Pr1xzpd4nL++NW8WN07f07mvvetX60efS5ZT1FCoATzgyijzAPwi3mqKPHMk9ixqQ9/+b34wQgyVJLHvZOnI9MCjaJz3JtfZvxOBKf164lOVu7PDK5ntsyP7e4yPGc3XaOa0/yeWR/jpjX2/rMtn/pF3N8ZJKfmPEcjXfoLi2XSdW4L5fM697r1t3p9P8jifxoz76feX/EpXM7t+GluGbun0cmeZnApK3G12aX1vcW14j4TP12Hhfm63j2b2y//L9Nvi9t55HX0T3O9JG41v/n2BYzOrsPTVv/+Pvs91hQSntvhYu4p99lubPHcavfxdGc5BnnlvUUKQBOODqIPkKS0XGbmfi0uKOOHsuc2GzJQf4dk9ZR9n2eqO25q9KSjyRRSba2yUHkGMaELLK9MWY2H8f4c5bdJheJTULXjq2tY4655sjy2YcxNra22CY8kf2c2yIyMRrX19x7PJfM693OaSY9475fWyb9cpskZflLyW5eH5dLbIsb5TzPE71sL+ua2zb7NffJS/0+2v7kXKTPjv0z5zHbHGNjnFg0837n53mZ+f2Rfdjb7kf6wLz/2W6L24yxaZOxXbO/Y7teaofEzOdje188+hqxbecz9Nu01xg7LnvpGGIuVOTnFje6p8/t9ej3SYzrixazaedqHmfyb36e45p5/TH+fjxfOaYWH/OTFnv6XXsfZT3j+307lvxujL3VpjH3vfTpFrd5xrllPUUKgBOODqKPkAF/3Oa1xOOIo8cyJ5RxK9mPef9vbSfrGxOt/H9Mgpp5G9eS3fk4tm0lqb51LJt7k+ujy2/xOfe3Eqw5+Yxbx3Xv8VwyrzeSSN7an/H8j3KeW/xmToxv9Zv0xTH+VgFtnhjf6sv5feLTDreOeZ7kpZ1aXGQ/9sZG+kz25dKEP+7pA3M73mr3MXZrz/x7bf9Gb32NiM/Sb7Ptcd2ba4WNSB+bl2lxo3v63B7PeJ/EuM5oMZv5GK/1o7w+Xgf2jvXj+je3+kQc7Xfb+y7L3erTiZ3Xf6tgPve9a/30WeeW9RQpAE44Mog+yjzQ30oW9zp6LHOylWRlz74kWRmXixa3SSIzxu493nH/sm8tJubjiD0TkNG8jvzc4i45unz27ch5z/GM67+VUN57PJfM68152dPO84Q99uzT3HeuLZNkdYzdk9TH/L65dl7yu739qr1PLhWk5na9dX73mNe5p703uXaMy96a7IyxmyPH8NbXiM/Ub+e42Htu5uvOrX4wt/OetjpiXv8j3icxrjNazGYevy+9pzfpZ+lvW/yeSfW4/thb3DjS7za3zulo7qtpixa3mfvetZzkWeeW9RQpAE44Mog+QhLQZ23v6LHck1COy0WLiXmiduSpkXnZS8ndfBxxpAAQ97RF3Lv8LfP6k4y2uM2z9ufseufl9vb7+f1ybXvjhOpW8jwbl73VtkfkOLf1xt4+vGcic8vZcxXzft+ayIyxkcl9i2tWXCP2tsW83Hvst/eMcUf7wdl23mte/yPeJzGuM1pMzE+X7G3LsdiVc9ZiRuM2osU0z27/GNcfeQ+2uDjS9551bllPkQLghHsSuKOS4Ix3VOJW0nfE0WO5J6GZk9cWE3Oyfuuu02xMwi/t33wcRyY5m3vaIu5d/pa3PLfXnF3v2eXm4740+U1Raow7ehduvAP5yGvA3uPO/o5x6ffXkv89zrZ5zO/vW9epMTaOvM9XXCP2tsXZ5V6p3x69dozm47/VD862117PeJ/EuM5oMTG35d6xZm6XFjMaY/fEb+b9e3T7x5Frw5G+96xzy3qKFAAnHBlE7zUP7mcm09ccPZZ7Esr5WFpMjHFH7xLGuPylRP+e49jcu45H7MM1b3lurzm73rPL7T3uef1HJ7rj8tc+NnDU3uNOMj4/Rp4kPcff4vfYu+1mfn/f2o8xNlrMJa98jTi73Cv126PXjtG8f7f6wdn22usZ75MY1xctJua2VKR4XJHiWeeW9RQpAE44MojeY75buPcz0UccPZZ7Eso5UWkxMT85co9Lx3PPcWzuXccj9mGT85j+krukOeYkauO645nn9pqz6z273N4+PT5O/QhtG9fkvZx9zXFl4pL9nBPuuHbcuave3i9ZT/rCPRPY2Nvmkf0fl82xtbjNGBst5pJXvkacXe6V+u3efWnm47/VD8621xGPfp/EvK4Wsxnjcm1uMbP5PLeY0Ri7J34zn+t72j/tmI9cZB3Z//Sb1u7X+sTRvveMc8t6ihQAJ9yTwO01byOuDexnHT2WexLKrHtctsXEGHOvS8dzz3Fs7l3HPcsn8UpRYm7Ta555bq85u96zy+3t00fabo+2jVkS6iTOrYh0ya3jzjpbcWOTbeWx6D0FzrNtHnN73rpejbHRYi6Zl73Hpf5xti3OLvdK/fbouDCaj/9WPzjbXkc98n0S8/ItZjMXHG61Sa7vY3z2u8WNxvhoMc18ro+2f4oSKbJea9vZteM/0/cefW5ZT5EC4IQzg+gR7c7AsxK3o8dyT0I5J9ctJsaYeyU5bNu45zg2967jzPJJXs/eSX3mub3m7HrPLre3T8/98R55v7ZtbLJPRwoTo73HnST8WqKefby1rrNtHnN75phb3GaMjRZzybzsPR59jTi73Cv126Pjwmg+/lv94Gx7nfWI90nMy7WYzdyeuRZcmjDn9fla8ej9Gc37trf9UyQ/UpgYXesT9/S9R51b1lOkADjhnkH0lpagPHL9s6PHkgF+jD8y4Gfd47ItJsaYJBwt5l73HMfm3nUcXT4J2Bg/SjvlblbWkXOaQtdbnttrzq737HJ7j3vuj896LDjnZdzOKO/1PFmRhD/7nX04e9ybnPus81Kyfm2CdM+25/bM8bS4zRgbLeaScblXu0acXe6V+u3Ra8doPv5b/eBse93rnvdJzPEtZjQfZybMuaZv5y/bapPs/HxtPzbjMtFimvlc32r/lqOM0leyjhxL1p34I9eGe/re5t5zy3qKFAAnPGIQbdrg/+zB9Oix3JNQzolKi4kkb3vi7nHPcWzuXceR5VuBIglYJreXJilveW6vObves8vtPe6jj2Cf0QoU2Z+cz0vv67PH3eRR7Pl9F894wmjezq32HGOjxVzyyteIs8u9Ur/duy/NfPy39u9sez3S0fdJzLEtZnatYNmkn2fC3dY1m5dtMc18rq+1f8tRIsd17c9/Hrk23NP3mjPnlvUUKQBOePQgGm3wT4LyrLu7m6PHck9COScKLSbmuGvJz1mPSIzvXcfe5dM35klZksIWO3rLc3vN2fWeXW7vcc/rz523FnfWvB+R4kSLHZ097mtakatNfu7Z9vy+zfG3uM0YGy3mkle+Rpxd7lX6bRy9dozm/bvVD8621zPsfZ/EHNdiRpcm+JekzY+M//PyLaaZz/W19p/P1d4iypFrwz1975oj55b1FCkATnjGIDrfYTlyB+UeR4/lnoRyTlRaTOTpgDFuz4T8qHuOY3PvOvYunwnYGLe3v73lub3m7HrPLrf3uPP+GuPyZEqLOyuTx3H9zz7uW/asd+5rR+42zu/vnIcWtxljo8Vc8srXiLPLvUq/jaPXjtF8/Lf6wdn2epa9+zPGRIsZjQWK9NdMmvP+GgvQicl141abNds6Ni2mmc/1tfafiyx79/PItWHenyN975ZX62tcpkgBcMKjB9G5QBFvVeE/eiz3DPJzotJiInePxrh4dHs8IlmZJ3RHJ0p792GO23M3Pt7y3F5zdr1nlzty3PNnljP5bXFnzP19713Rs8d9y9wubb1zTCYlc0wzT5wj62qxmzm+xVzyyteIs8u9Sr+NI/sym4//Vj+49zr6aPOxXzp/Y0y0mM1YVMu5u/RRr3uM+xItptl7vDHGHSmOzdfCa33inr53y5FjZS1FCoATHjmItgLF3knoIxw9ljkBPTLIz4lKi9nM7ZLJ0iMTu3uOYzO33d4J3WbvPsxxe/pHJnHzXa9nnttrzq737HJH+vT8CPAjn2Ca+/ueIkX2ff5oz6Xjzvshv9v7vpiPta23Tf5vrT/tNe9z5Fha/GaObzHXvOo14uxyr9Jv48i+zObjv9UP5m0dvY7e8oz3SYwx0WI247XgGR/PiXFfosU0c/tf669j3N4ixfzUU1zrE0f63rPOLespUgCccGQQvWZ1gSKOHksG9TH+yCA/JmrRYjZJOuaJT5LXvRO9W/t1z3Fsso/jOuLIRGHvPhz9uEfi26Txmef2mrPrPbvc0T4998u0XdbRYkfpi9fu+h79uMd8vLeW264fmSzc2t/01blodamvznfpr02qLhUo4tY+zfEt5ppXvUacXe5V+m0c3ZfRfPy39une6+gtz3qfjDHRYjbjuTrSlkeM+xItppnP9bX+OrfNtfZMW85f8rpnuXl/rrXXs84t6ylSAJxwZBC9ZBtcR29doIijx3I2AY8xUYsWM5r3LZKIp+3mhCTJRu7YjEnItcnKPccxmo8pyVKKBNvvs1+X1r13H9pkLNsdt5OY/Dzvz7jcrXM73/F61Lefn23rs8sd7dOtfSPHnzbN77fY9Km8T8fk+9L7tvXf9N0xMd7WNxYG5n1px51lxpjIceb1sd9n3/PanJxfa5N253Pe7/x/vIZl/fNy83t0NsZGi7nlFa8R773fxtF9Gc3Hf6sfRNY/LnPkOnrNM98nY1y0mM1csMx5GN9PjzCuP1pMM5/ra+08n9v0v7w2tuV2rsa+OffTa30iy4+x6QtjX94889yyniIFwAXzwH1WW/eY3G8yiGfQnAfSW+4daI8mo3OSci2hmWXd47ItZpaEpSXie6Sd2zrjnuMY7e0nbTJ0ZB/axPGa7c7S2Oa3zu2cHG7Ots3mbFufXe5on44kskffe5u0dVtnjJPCPbKv6Svja+24s79H173Jcbakf3O0Lbb1ze2en9v6N2NstJg9Xu0a8RH67Zl92czHf6sfxLy9S9p19Jpnvk/m+BazyXrO9NGco7R92vTWsc/Ltphmbvtr/TXHkX0a429JgSYFp/G1W32itdXcB595bllPkQLggr1J0y1t3S3urCPJYzMf5631nU3AI+sel20xTZKJVti5JElUJvXXkpB7jmPW7ujMEjMvd3Qf9hQqktxlPduxj9vY01fmcxRJBFvsXkePc3N2uaN9epRt7J1MJC4J+LXJw96+mz479pHxPFw77hxrO2eX5FzuSc4z8d8zGRk/CjIXuW5NRMbYaDF7vdI14uxyr9Rv792Xcdlb/WBz9jq6xzPeJ/MyLWa09z11zfh+m82xLaaZz/Wt/pp+s6cwlvYez/34u1t9Yu5DkX7bYp9xbllPkQLggnngPqutu8WddSR5bObjvLW+OXm4ldCM5kSixVyT5CgTi6xnTsjzWiYp4yPC19xzHE32Ldsfk7fsY/Yr624TgjP7kPUkUZ2TxGzn1qRrj22yN7bvZypSRNogk6Ec9zypSLvn9fz+SFtnn9Ku4/rSxtu62jJHZAKU9rn03kifSUxb9pIc3/Z+G9eXNsixtD59xLjOaDFHvcI14iP023v2ZT7+rKvFNWeuo0c88n0yLhstZjS2S85Ptjcuv9elQsUc12Ka+Vzv7a9bXxvbMceV83fknF+S/Rj7cf5/re8+4xrIOooUAAAAT5KJ+zZhzv9bzCwT8hTUUggYJ9xxb7EGXp0iBQAAwBOM38dw62mAS+anLh7xpAK8MkUKAACAJxg/srD3oxSz8aMiihR8BooUAAAATzAWF85+/8z4cZFQpOCjU6QAAAB4grG4sPf7KEb5ssdxHXkyo8XBR6JIAQAA8ATjxz0if32mxTX5Pov5L1U84q8BwatTpAAAAHiCFCXGIkOkcJHX28c28lp+Nxc34syTGPAeKVIAAAA8yfydEmec/dJNeI8UKQAAAJ4oH9NoT0fckgLHH3/8UdcJH5UiBQAAwBvI90z8+OOPf3333Xf/9H0TkddTmEhR488//6zrgI9OkQIAAAB4CYoUAAAAwEtQpAAAAABegiIFAAAA8BIUKQAAAICXoEgBAAAAvARFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAAAMBLUKQAAAAAXoIiBQAAAPASFCkAAACAl6BIAQAAALwERQoAAADgJShSAAAAAC9BkQIAAAB4CYoUAAAAwEtQpAAAAABegiIFAAAA8BIUKQAAAICXoEgBAAAAvARFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAAAMBLUKQAAAAAXoIiBQAAAPASFCkAAACAl6BIAQAAALwERQoArvrzzz//+uabb/766quv/uHXX3+tsbPvvvvui+V+/vnnGvco8/Z++eWXGnfNDz/88MU6/vjjjxr3SPN+H/Xbb7/V9X406T/jcX/77bc1bpb2GZdLf06/brGP0LbX4q5JvxvXkX7Z4h5p3u+j0o/bet87/a7HPop+1+l3PfZR9LvXpkgBwE1nkp55mbcY0FM8Gbf5/fff17hrvv76638svzcpvJcixX5zW+0pfM3LvEV7pe+M2/z9999r3CXzBGVvYfAekvbL9Lvn0e8u0++eR797bYoUAOySCf84QP/44481bnNv0nJWCijjdo88CTEXOfYkhI8wJ5VHPPsu2atJPxqPP0Wla+f4EYWrM/IUz7jdo3cGx/dPjrHFPNq9SftbvV9W0O+eR7+7TL97Hv3utSlSALBLJsJJHsZB+lLhYU5Yfvrppxr3DNnWuO0jiUQSq3HZIwWOt5D9mc/BW9wlezXzOb6WiI9Fq7TdWxV05vfLkcQ753lbLt7i0eejUqQc9/Ez3FXU79bT7/S7FT5jv1tNkQKA3ea7NG2gTrIyJkpvfad/Tniy/RbXjElW7uy0mJXmJy5uPc3ykc1P6rRiTQpUY8xb3/mai157vyNl3u+871rcKvMdyLxv8r5rsR+NfreOfvevjl2/ezufud+tpEgBwCHzRHlOROa7PisSjtxpGvdhz0dN5gLM3gTrrcyJXJLWtyz+vJo5cZyLUfOdvRVFp/lR7b2PXo8TkiNFtrcwFyHj1d4rz6TfraHf6XcrfPZ+t5IiBQCH5A7CmAzl/xnI87s5UdqbpDzaXHDY88TBfBdoO6ZXMLd57Cm8fHTzI7gpkG2/m4tlq9prTMDjVr/KuR7jX+1pmbnNV73HV9Lv3p5+p9+toN+to0gBwGHzXf0tsRgn+plUJwGZl30r492P7EuLGb1CceWSOfEbk9PPLAnwfJ7T514p8c1dt3Ffbj2CPb+3XqkYNRf/0t6vVMx7K/rd29Lv/k6/e1v63VqKFACcMk+c5wH9VnLybPOdpexfi4t531/pcc75ONLuLe6zms9dCkxjsSxJ/crEMttOcrvtz63zN76vsu8tZoX5OOLae+qj0+/ehn73Jf3ubeh36ylSAHBK7niMA/joVmLyFua7S9eejpifAFmZ5I1aG7/SnaZXkXM7t9PmFRLLsX/FpXM499lXemJmbuNr76fPQr97Pv3un+l3z6ffradIAcBp813+Tfvm8RXmZOlS8WG8Y5JlWswK412meKUk7pW0u17xKn8mbu/j2POjz1muxb21+e7tKxXyVtLvnku/6/S759LvXoMiBQCntSLFKzxFsZm/Eb19BGVOSF7hTlTMbfsqCegrupS0v8q5jPGv4mRfW8xYlHqV91Fr21cpQq6m3z2PfneZfvc8+t3rUKQA4JT2UYTNK32nQz7nuu1XS4Ryl2f7/aVk6q3NbZv9epW7TK/o0uPPr/QZ5/kL5eYJRZLj8ferv9NlM7ftyi/lezX63fPod5fpd8+j370ORQoATpnvgox3H/L/JCFtubd265HSsYjxCh/1SLuN+xSvksC9oiS/Y1uNd/DilT4iM75H5r42J/WvUJSa3zvpl6/yvl5Nv3se/e4y/e559LvXokgBwGHzYJ7HIefkKXck2rJvLUnGuF/jnZH5iYVXeKxzfLIjkoS2OP65oLO11VhAi0tf3PbW5nM7JsDjHbxXOOeZNIyTjPDY89/pd8+j312m3z2Pfvd6FCkAOGQezMc7JEk2xkF+fsRzlezjtk9J8rbXx+99GF9fJUnRtj+Rdn6FO0yvak6Ct7aa2zFJ/LzsCtm/cb9yN3H73aXXV5nfy2nrFvcZ6XfPo99dpt89j373ehQpADhk/szmlijFnJRkkj3ePVllTuK2O03jHajVScl8lyxeIXl7VfNTMPPjxGNhKl7lMegxGd6eNpqfQlr9npmflMr75BXex69Av3se/e4y/e559LvXpEgBwG5zctESofHphNgSk0uSDGwJVooazyoWjAWAbCPb3X6OvY/IPmt/5yTzVrt9dmOBqRXD8nNeH9v01jnO77ekOv3lGUWirHPcp+xn+tD285Hz/oz9zTqPtttn8ox+l+vqdn3K+p/R3o/sd8/YX/3uOtc717vPRpECgF3mJKglSlvclsBu8iTDHLeZJ+fxjLtA492SLcEZf27LNM/Y37n4c6lt+bu5EHap/ec7ZElwW9xm7rdxre+eNb6P0g/H7R5JvJ+xv+NkKJ7xXnyvntHv5qfP4lnv/0f0u2ftr3532TP6Xbje6XevTJECgF3Gux9xLblI4jDGJrm4lMSOcWN8i71Htj9vY/t/kru2TDOuY3PP/ma/xiQuUrRosfx9kjS216223+68bS6d67nPblKUavH3GN9LYz/McbX45hn7O0+GksC3uM/oWf1ubvPNkcnbXo/od8/YX/3uMte7v3O9+3wUKQC4aU4Qbt2hifmJgyQqLW6MGSU5a/H3mPfpzLba8nF2f/O4a1vfEXvOx0cxJ+Hpmy1uk/MyxicxbufqUhJ8a1JwxrxPmyMJ96P3N484t/Uddet8vFfP6neXJv3PmCw+ot89en/1u+tc7/7u0fur370+RQoAbhrvfsSez2y2JwTagH5pkv6Mwb8lJkfvnjx6f9u6jvosRYrcrR2Pe+/nmeeJ1aX2mvvrpsXeq/Wjo0/QPHJ/L00+j/qISfsz+92lydKz3tP39rtH769+d5nr3Zdc7z4XRQoArpoH80tPRDTzdy2k2DF/7CM/547KGBfPGvznz6Ae+ahHPHp/5/Wc8awJzSuZi175f7tDeMme854JWNpyjIs57hHm90aOp8Vd88j9lbR3b9Hv0hfmQvCz3tOP6HeP3F/9rnO9+2eud5+LIgUAL2NMQJKQtJh7zXd0kgy2uD3eYn9ZYzuvZ5LpPdJftm3EkeJf8+z95fky4dnO49675kc9st+9xf7yNrbz6HrHq1CkAOBljHfn2u8fIUnNto17E+u32F/e3vgZ6rN3iG/Jnc1tG3FPkest9pfnGx/vz53eFnOvR/a7t9hfns/1jlekSAHASxjvyp39ErZb5kdOk2S3uD3eYn9ZY3wU+J4+cs34OHaKXS1mr7fYX55vfDLryKP9Rzyy373F/vJ8rne8IkUKAJbLRy7GJOZZn/Mck+p7HhN9q/3l7eUO3/a0Tf695+NAl4wFrmifF9/rLfaX58tka+sPz7o7/Mh+9xb7y/O53vGqFCkAWGqe8N/7WdVL5qcozj6e/Fb7y9sbE+BIn2lx9xr7zz2J9lvtL881TvhzPt/iKYp7+t1b7S/P5XrHK1OkAGCZ3FEZk44kMy3uXnNyczZReqv95W2lL6TYtJ3XeFbxaf7LMGeKZW+5vzxPJvfj013xrInXI/rdW+4vz+N6x3ugSAHAMuNf2kgyc6ZwMMs6UpTY/p+7fmNhIc4m1s/YX17D2EfueRx5lH649ZH8f+w/cU+R6xn7y9tKn9jOYc7noz429qx+96z95e253vHqFCkAWCbFgnyJ1iPvxo2PIjf33IF5xv7yGtIvcpc4yXX7/RnpK60PRpLuex6Tf8b+8vYycXt0wfOZ/e4Z+8vbc73j1SlSAPChzI+XjtyB4a0kIW99MDLRuydhh0v0O1bQ73g0RQoAPpQUKcZHQ5Mg5S6MJIm3lEfh0/e2fpg+mcef86RPi4dH0O9YQb/j0RQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAAAMBLUKQAAAAAXoIiBQAAAPASFCkAAACAl6BIAQAAALwERQoAAADgJShSAAAAAC9BkQIAAAB4CYoUAAAAwEtQpADg6b766qsvtBj+ld9+++2L9vrpp59qHPcb2zlazEf0WY/7Fczv7++++67GfTSf9bhXybgxtrdx5Lb0ybHNWgxvQ5ECgKcbB30D/22KFG9nbOdoMff6888/v0h+v//++7+91mLfynjM0WJ4jrearP/+++9/ffPNN//YzurriCLF21KkOE6R4nUoUgDwdOOgb+C/TZHi7YztHC3mXj/88MM/bWf1OZ33p8XwHG81WR8LFJtsu8W+BUWKt6VIcZwixetQpADg6cZB38B/myLF2xnbOVrMvebEN1ZP0Ob9aTE8x1tN1sdtbFZeSxQp3pYixXGKFK9DkQKApxsHfQP/bYoU+81t1WKuGZc9s/we82Qhfv755xr7Vub9aTFcdk+/e6vJ+jzhinwEpMW+BUWKt6VIcZwixetQpADg6cZB38B/myLFfu+hSJHvn9g+8vH111//9eOPP9a4t/QWx/2RvYciRQoS26QrH/345ZdfatxbUaR4W4oUxylSvA5FCgCebhz0Dfy3KVLslycSxrZqMdeMy55Z/r36rMf9KPf0u886WVekeFuKFMcpUrwORQoAnm4c9A38tylS7Dcn4i3mmnHZM8u/V5/1uB/lnn6nSPG5jnsVRYrjFClehyIFAE83DvoG/tsUKfZTpDjnsx73oyhSHKdI8bYUKY5TpHgdihQAPN046Ocz+S2Gf0WRYr98v8PYVi3mmnHZz9Q3P+txP8o9/W5+f3///fc17qP5rMe9ylykWP1lve+BIsXrUKQA4OnGQf8j3z3LF9Ml8c6kbzveb7/99m9fmvjHH3/UZRpFiv3uTSrHZT/Tnd3PetyPck+/+6zvb9e1tzUXKdL+Le69yxcTpwAzvyfzc9ogv2/LNfe8r3ksRQoALvr111//NsHORHscuLeJ996kZ1x2nhBl8p5EYk4Oso2jCcYo603ikqLBvP/5pvu8fvTb7sd9HBPsbGve/2bv9s4m8/M+5By1uPds+4sFYyFoj0ttOMZkvePvHtk3x+WjxVzS3ofpw7mbf6T4NRr35ZnH/VE8st/den9nW+26m+0fuWbN25nP8y2P7ndvddzvXdp2m3SP/S3/z2v53Z73YNp3bMd5vM75zTiYczpuI6/NsUdkvekj8/5HXjvTf8Z1jPuW/4/732Qf0rfG9V2S/RuXbTGzuV8f2R6XKVIA8E8ywM6J4iVJEG4lNGN8koDt9SRbcxIzOzrgZ1+SZLV1NTnOvQnTmMAk0cpr2d6tYxjdaquYk5420ZkluR+XGdv5I5nbZq9LbTjGPLNvzsu2mNmR9+Ell/rBpZhnvCc/gkf2u0vv70w8c10Zf9ekT+yZpM7b2XtNeFa/e6vjfs/SJu3Ym7TZtfaY15X2z+sZ73J+xt81GVPmdV6S/cj2bhUMRrnWtHU143IpguS1LD++fk2uWXv6ztwuLWaU98p8vdzamfsoUgDwhTbo7rElnM0YlyQgycKeJGmT/blVSMh+H1nnaG/yO64//8/dve3nJGf5eVtP/s3Pc1smbl7v7FIyf8m4H/GRk/m5bfa61IZjzLP6ZszLtZjRfE7PyrG09c8xzzruj+KR/a69v3P9OjLBy3t8Xu9s3s6lvjB6Zr97q+N+r+ZC8x55H1661qd9x9i0fxuTrtmK8ddkO0fWOcr+tHXOxmWyvbGt0te2wkXkujQfe+wpuszXwBazSbunP47xe4+H2xQpAPjCPOjm5zEBiCQ7Y5Jwa8Iyri8J6ZiUzo+Wnk0wslxbJvs+7luS4nHfN3uSsTGBGZOytNGlRDHb2+I2ea3Fbloy3+Jijr11Lj6a8dijxVwzLvusvhnzMi1mk/4xJ/0tEd/z1EObLMYY88zj/qjmtmgxl8zv2Vw/xvOYth6vEfl/ezrs1oRo3s6lvrB5dr97q+N+j9LG83GmPcZref6fYx/fq9fGrfk9m3Oy/T/tnnWN688+zON/jDHNPJ5m/9JHxmtIxsd53yP7cWv9MS+z/f/aNWh+0iLLtbjR2EbRYjZz7J4cgv0UKQD4hySF46CbhKXFbZJczElsM65zk4Th2nJJaOZlLhUCNltSlsTlVuLTHhW9tf45KYkk0S12NCdx2c8Wt5mT+UvxbVIxJvmfwXjs0WKumZePZ/TNOb7FbOZ+di0RTz8fE//sW4ubjevfPOO4P6q5HVrMJfP7e5Pr7bX373wduXV9PlqkeHa/e6vjfo/mYsytMTXj163zORcpNtnWpfdtXk/7jvG3Jt9ZJteO9Idb/SCxc6Fiz+R+jN+kDVrsaN7WWDhp5vdAi4m5T+7JAzhGkQKAf5iTmr0TnlvGdUYSmj2T6TlZupW4JQE6Mkmf13/reOcEJstfSvZG2e9xuVvJ5Z4iRUsmH3W+3pPx+KPFXDMv/6y+OcZGi4lM/sa4JNktbjT2l+z/rQJdjNvYlnvGcX9UYxtEi7mkTdb3XEvy+3m5a8vM27l23XmLfvdWx/0ejce2p+33aEWKa4WnzTxe5Ry1uFHObXu9mde/53jH+NhzHJECyLjcpYL/Zk+RIuPsGLOnD3OcIgUA/zAnNY+agIzrjL2FhKMJxlHz8d66ozMnMHsTsyOThZjj23HPd94e3TbvxdgG0WKumZd/Vt8cY6PFxJwA7z2vmSRuy5y5M/kq78n3YmyDaDGXzO/vTNL2TnLmItG1a9CR685b9Lu3Ou73aDy2PUWBPXIOx/XuvdvfikIt7h7z+m/1gzm+xTRzG9zq17eKFLlOjr9P/1egeA5FCgD+4Vl3CMZ1Rotp5o9kPHpCdCSJj7NFihiXu7dIMU8U995V+ojGdogWc83Z5Y/2zTE2WkzMSfXep2PGvnmrf8W4jWgxzbPfk+/F2AbRYi45et0ZzcXJa9egI9t5i373Vsf9HqVgMx7f3va/Zj6nR96rY/EpWsw9xn4Tt87nGBstppn73K02mPdr/F2eFBrbJf/fW9zlOEUKAP4hBYk5WUqh4t6EcFxftJjmaIJx1Lz+W0nz0cRqNC53azvXjvtZhaT3amyLaDHXnF3+aN8cY6PFRNYzxr1akeLZ78n3YmyDaDGXHL3ujOb+ce0adGQ7b9Hv3uq436O5+Bc55nuu7XOb5ecW14znNVrMPeb13zqfY2y0mGbuc7fa4NJx5zzMT/M8opDEZYoUAHwhdwbmuyiR4kXu4O/5vPtsXleLaY4mGNdkXUkEcwxJROaEI24lzUcTq9G43K3tXDru+VHT+OyJ0tweLeaas8sf7ZtjbLSYmNd76xH6zfievdW/YtxGtJjm6HF/VGMbRIu5ZG7DPedrk/Yel826Wlwc2c5b9Lu3Ou73av4yxk2eIjlznZ/b7Mh7NedmXLbF7JF8IR8bzbZzHFlvyy9unc85vsU0c5+71QaXjnt+kif50Lgcj6dIAcA/yWR4fqJilAl+kqa9d3nm5VtMczTBGCU5SlFiTjquuZU0z+s6kiiPy93aTjvu+VHTTV7zJMW/0mKuObv80b45xkaL2Yznec/5zXtxXPeej/+M8dFimnvekx/J2AbRYi6Z2/DW9WCU9h6XvXYNOrqdZ/e7tzru92z+KN8o5yRtvPcjBnObHXmvzmNdi7kkRYns57UcYnbrfM7xLaaZ+9ytNmjHfemcHGlPjlOkAOCiJKHXEo0kTXsG6nm5FtMcTTAik/n5rsdet5LmOYE5kiiPy93aznzcSZLakx+bHG9bz2cwt0WLuebs8kf75hgbLWaTdY2x1yZ/7cmnPf1yjI8W05x5T35EYxtEi7lkbsNb14PR3Deuneuj23l2v3ur437vMoal7cfjnaXtbhUr5jY78l7N+sdlW8wsNwWOFCZGt87nHN9imrnP3WqD+bjnQtxsb8GI4xQpALgpA3EmypcSkFvfizDHt5jmaIJxLaHIvifxyzqy3hzTvP5bSfOcwNxKrEbjcre2M+/XLMc5F2Ie9ZdY3puxDaLFXHN2+aN9c4yNFjOa+1reYznH2/ssE5lsc54o7i1YjctEi2mOHvdHNbZBtJhLjl53Rmnvcdlr16Az23lmv3ur4/4o0ubtWj/K79uyMbfZkffq3A9azCb7ea2InnVl29nXnLfEz+u/dT7H2GgxzdznbrXBvF+jHOO8vrzW1sP9FCkAOCQJaxvIryWpc2yLaY4kGK1AkcJE7u4ksW7LzOvPcbW4zdHEajQud2s7836NtjZIojdOFHKs20TiMxnbJlrMNWeXP9I3Y4yNFjO6lfg3id/bB+ZlW0xz9Lg/qrENosVccvS6M0p7j8teuwad2c4z+91bHfdHdKk4FJfu5s9tduS9Oo91LSYu9ZfcELhWOJ/Xf+t8jrHRYpq5z91qg3m/Nmn3rY/P7ZocY14P91OkAOCUVhS4lCzNcS2m2ZtgJHmYk7drjypvjibNRxOr0bjcre3M+7WZjynJ0fj7vV9295GMxx8t5pqzyx9NfsfYaDGjvJfahOSSnPs9E8XNvHyLaY4e90c1tkG0mEuOXndG8wTp2jXozHae2e/e6rg/slYUuHSDYG6zI+/VeaxrMTFvI33nUh4wmtd/63yOsdFimrnP3WqDeb9iPqacg/Gp0vz+0o0QzlOkAOC0vUnQGBMtptmbYOSOzRi3N/k9mjQfTaxG43K3tjPvVyQxbbFzwnpknz6C8dijxVxzdvm9fXMzxkaL2SThHSeKKUbF3P/yc7Z7JkEe1xMtpjl63B/V2AbRYi6Z2/DW9WCU9h6XvfZ+P7qdZ/e7tzruj64V5Vvc3GZH3qvzOW8xcXb8mdd/a7kxNlpMM/e5W20w71e0J0LO5hzsp0gBwGl7E4AxJlpMs3f9czJ27XO6o3n9txKNo4nVaFzu1nbm/UoieOlu5Rz72T72MSfJRydO47LRYpq9fXMzxkaL2YyfQd/zWf8zxn2JFtMcPe6P6p5+N7fhkQnOfK3LulpcHN3Os/vdWx33ZzCPRS1mbrMj79U9648xJmNPi2nm9d86n2NstJhm7nO32mDer2u5xBybgl6L4xxFCgD+IZPbDOJ7J7nzRz4uJQBjTLSYZm+CkdfHuD1Fikwq5olGko4WuzmaWI3G5W5t52hilUeuj8R/JPM5ufY56GZcNlpMc/QcjbHRYjZj3NHj2WvcRrSY5uhxf1T39Lu5DW9dD0bzte7aNejodsbYZ/S7tzru9yhj1pFjGj9uEC1mbrMj79W5f7eYGGP2FikymR+Xi1vHPse3mGbuc7faYO9xx/zkUf7/mW4QPJsiBQD/sP3ZsyQbt5KGDMbzJH/Vd1IcffQy8fPjsnuWmxOYI0nluNyt7RxNrHIu5uPZ87ngj2Au0By9AzwuGy2mOXqOxthoMZsx7sjE4ohxG9FimqPH/VHd0+/mNrx1PRilvcdlr12Djm5njH3GeX2r435vxu8B2fMdH/PNgUt9b26zI+d0HutaTMw5wLXzkuMan9bZu1zM8S2mmfvcrTbYe9ybuY2f8QTSZ6VIAcDfzIlPZMDO6+OjzEk08tqRpxDGuGgxzd4Eo03Ssz/j3cDE5Oc5CRmXu3YMMS97JFEel7u1naOJVcznL+enxX00c1tFEv2tz27nfewLo3nZFtMcPUdjbLSYzfzeyrqPfJxgj3H90WKaM33zI7qn383L3roejNLe47JZV4uLo9t5dr97q+N+T9JP5nbPmJS+NB9jihlzcSwutcXcZkfeq/NY12Ji3kb2fe432e+8No6183h963yOsdFimrnP3WqDvcc9mp9quTTWcIwiBQB/k2Tp0l2OW5JkZfm23pjjW0xzJMFoj5Besz0tMiYlt5LmOYG5lViNxuVubedoYrWZ92/vcu/dfNzNpb/2Mse1mOboORpjo8Vs5nXvlfdh3sN5L1x7P8a8bItpzvbNj+hsv5vb8Nb1YJT2HpfNulpcHN3Os/vdWx33e5NJ7Txp3yvF6bbOmNvsyHt17tstJnK+50n6LSm05JjH126dzzE2Wkwz97lbbbD3uEfzNnIub11/uU2RAoAvZMCdB+prkpy+yoRoT6Fiu9Oz7fOYyN1Kmud2uZVYjcblbm3n6HFvcsdqXC7H+ug78K8o53K+Gzm79HnpOa7FNEfP0RgbLWY0J/FH5dxfm8DM8S2mOds3P6Kz/W5uw1vXg9F4vYpr16Az23lmv3ur436P0pdyjGm/8TgvSdytO/Zzmx15r85jXYvZZIy59T6IrHM8b+Pvbp3PMTZaTDP3uVttcOS4R/MNnhRiWhz7KVIAUG2PaGbQnhOnvJZBODFt2dm4bLSY5miCEUmYsm9z0pR93nOH+Zo5gTmSKI/LZT0tZnPmuDeJHZe9ta2PJJOjHO/YX9MPcjf70sRpbKtoMc3RczTGRovZpI+OSW+OYc8koLk0kZnjWkxzT9/8qI72u7kNj7xH5/f3tWvQ0e08u9+91XG/ZzkH6TPpO3Pbp+CV85Pf7xnH5jY78l7NuRmXbTGz7Ff2b3wfZJ9zLPeer3FfosU0c5+71QZnjjtyPuY86aP20beiSAEA8CKS7I6TkzbJbVKcS2yb2LR4GOl3wCtRpAAAeBHjF+PlrmSLuabd0WtxMNLvgFeiSAEA8AJyV3qc5J19XPjsI8t8Tvod8GoUKQAAXsD8+elMHlvcLXnUflxPi4GNfge8GkUKAIAXME8W934vwGj+CzftT2DCSL8DXo0iBQDAC5gfu89n/I88ej9/m3+cvSvO56HfAa9GkQIA4EWMfwJyk8/65+72/Cd/82WFmUzmSw/nLy2MM3fE+Zz0O+CVKFIAALyI+U9BnpGJ46+//lrXD41+B7wSRQoAgBeSCWMeoW93qa9JfJbL8m29cI1+B7wKRQoAgBeUSV8enc+XEM5/3nGbHOb1PHbvDjaPot8BqylSAAAAAC9BkQIAAAB4CYoUAAAAwEtQpAAAAABegiIFAAAA8BIUKQAAAICXoEgBAAAAvARFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAACAl/TTTz/99dVXX/1Dfm5x8Ej63VqKFAAAwEsyWWQF/W4tRQoAAOAlmSyygn63liIFAADwkkwWWUG/W0uRAgAAeEkmi6yg362lSAEAALwkk0VW0O/WUqQAAABekskiK+h3aylSAAAAL8lkkRX0u7UUKQAAgJdkssgK+t1aihQAAMBLMllkBf1uLUUKANjhjz/++OuHH3746+uvv/5bwpJ/v//++xoLz/TLL7/89e233/4jec7/f/311xoL753JIivod2spUgDwVJk8ZXI/Tqq++eabv3788ce/TfzbMrM5Wcj6Wtw1mdiN68jPc8x33333j9+PCUlit+LEKMcxLt/kGHOs47q3ZVPkyLr//PPPuiyPk/Pw888//+08jOdyOwdtmWY+j0eW3aT/juuY3we//fbbF7/fXk8/mbe/2bMfWe/8Xoz8nNcVOs5Ju23v8fk6kdeOXOsu2fpv+muuHfM28rsz23jWen///fe/Hfelvpa+2JZrTBa7t+h3uebk2nLpupFt5Fy3Za951nq3/pzjH9c5jrdtuUa/W0uRAoCnSIIxJx9HJdHIupLQjK8nIZu3d0sSlHEdrTAwJjZbQpJEcFxulHXO69gkWZoTpUtyPEms2nq435xsnrFNqub+cK0PXDJOKPIemX9/qUhx7f10LaHP+uYJ6CWJOzKB/Kxy/Ui/2tuuceY9nu1kwtbW12Sf2rVt9qz1HrnuJW7PJHp+/+bnFvcZvFW/iyw3Fz8uyXVwb0HkGes90p/3XuP0u7UUKQB4uNytGAf3s7YiReSOy/i7I3dEksCMy16aWI7J9ZYcXUumLiUtSYD2JmGjbHPPRIB90pb3Fso2Y1I7TxD2JucxFznaBCLbGmO2u9Lja7N5HZusv8XfIiG/LG1z5v0dR65bOe9HJqOb8brZPHO9R9sl8dcKbGGy+Hdv1e9y3RzHwr2yb9fGr2eu98x1/lab6HdrKVIA8FAtUU1iksnZFpNJ3Z67KWNSPE/cLhUamrlocik5GROo/H/+eTyGHGebnOb1bZlNJgTZ5hif45kLLzEeM/eZ2zf9LYnmeB5yTtOXxrhmLFLMyWsrNFwy71PrQ3NfH7c3H0MS9HHfRq1AsfXjLenPv/l57OsbSXk3n8O8v9PW43lIu+Y9PxcDcv7aOZ9dmvBn2+N2sq55O9cm/W+53qxz7GuRbbT2G2Nm8/vts/bLt+h3Wb5N+PdcN65dB99yvVl2Hm/TP9Nv5j56rU/rd2spUgDwUGNyEUmsWlwkiRiTqSQWLW4zJ157kq6YJ6FbQjQb931MZvYmJ1lvS9QvbS+SYM7LSIbul3Yd2zRtvDchTdJ77Zyl343rTr9scc14rrOdFtP2fYu/dgyjeR1xLdmP/H5eZpwA8Xfb+zzn/dY1K7HzdStPxbTYTVvm1rnPMrnW3JrQPWu9Y7/O/2/1m7TbuB/Xrnnje/NW7Ee2tfOz+l3MY2W2l6JBi92kb9y6afCs9c75xrV+GunrY1/N8i0u9Lu1FCkAeJgzk7dxMpXk4VrhIQnIuP5bCUkkWRuXuVY0mROeW/GzxI7LXpqEztqEcm8Bhm4+F7eS+hjP/62Efk669xQPkpSPy1zap9Yf8t7YW6CIMxOUmBPzPe/hz+hI8WY+77fadP5oT859rmMt9oi3Wu/efjr2texLiwmTxX/lmf1ujo8j27vkrda7Jx+I+fp6qb/qd2spUgDwMEfujo2SoG7LXJtMzQWHPUWAeZ+S2LS4mIsUSer2JvFzgSaOFBrmSfWR4gj/bJ6kt5jZXAS7du7nBHlPEWA+x5fWPyfRsTcBj7nPHy00zG23p8DDdWN7xqVzP1/j4to1a6+3Wu+R69a87KX9MVk8b2y3uNTvIuPpGLvnmrbHW6z36DVuXPbS/uh3aylSAPAw86C+d3IzFgeuPX4Z80Tv1l278Y73tbt1MRcpjkzO5gnurcdUZzmOcflb+8p1Y1vuTWDn4sCtu33jZH7P+RqLcdf6x7wfRxPw+SmPIwWOuLcv88/ma8ulvjW3/Z5C7B5vtd6jT4CNk8VLk0CTxfP29rt5/IlrBY293mq9RwupKUxsy17KOfS7tRQpAHiYeVB/RpFinsBduytz9C7f3oSumZc9c/d5voN9ZPt8aWzHZxUp5v5+7c70/OTFtf4x78et98RsXDaOThzbU0Etjv32XlvmuKMFpkveYr1Hi2kxLn+pGGayeN7efje38aMKk2+13qPXuHH5SwXmeRv5ucXxHIoUADzMPLna+1jneIc5SVWLGe29g52J4BYXtyadexO6ZjyGo8tu7r0Dzr8yn8s9d+/mu8K3zuE8mb+WgI9PAKWvXNuf+X205z2xaXcuW9wtc3/OelscX0qfSEEqE5r0h5y7uS3jUt96xHWkeav13uNSPzdZvO3efjePPY9q47da773aNvS7tRQpAHioMTHK/29NDudCwp7PNM/LXLqDPSYye+7yzRPbI4n8uFy0mFskRY8zFxz2tOX8JMueu3Pzx48u9ffxfXGrj6ffjes8UqS4Z9nRPe+FzybXn5zTuf9cc6k957gWc8Zbrfcel/qq62L3yH73rPf7W633Xm0b+t1aihQAPNQ8sF+bkOXu7Hy3Z08Sk4nguFyKES1mXO+epzruSajG5aLF3CIpepy5j+T/154GmIsNlyZMs/SRcbn29EsmE2PMpaLaZl7n3n2Je5Yd3fNe+Cxyro9MEEeX2nOOazFnvNV679Gu4+G6+KVn9Ltnvd/far33yNjQtqHfraVIAcDDzQlEvhwtE7PtLnPuUGfAnwsUl5LUZp5Uznew54nhnsfV70moxuWixdwiKXqsnL+xPbc2HftC+klLeI98vGGcMLQvJBy/pO1SQjya9zv71+Kae5Yd3fNe+OhyrRm/8HGWtks/yxNfabfE723PMSZazBlvsd68D1rMvVwX/+6Z/e5Z7/e3Wu+ep96O0u/WUqQA4OFuJVNN4rNcW1+TSeS4fBKz8fdjEWNv8nxPQjUXXI5Mcjfz52wlRfebPxq0x9yXbpk/WjInzGMR49qTRZv0u3F96ZctrpnfF3uKIs3cnx81uXjvLl3bcl6vPSGz99ryrHZ/q/W2mHuZLD6/3z1r7Hmr9T7j+qTfraVIAcDDZaI0J6/X5E7zkQLFZkza8v/xd+P29yYX9xQp5mWPTnRjfoT3GYnXZzM/cXNN2v9Mm6fvjusZP1o0Fw32rD8x4zJHihQxLhtH7zImfl5Hi/uM5olLrjN7CpJ7ry2PuI40b7XeaxPms0wWn9/v5vXvKabu8Vbr3fsl3UfM28jPLY7nUKQA4KEywRkLBLnLHHOylJ8z6N/zmOZ8l3xb1zzJ27uNvQldM99NP/LRlXjUHXD+lfFjFiliZQKVJHksBuX/ee3eSdtYDMk6t9fHRHd8/Zp7ixTzXcb0zRZ3yb19+SOb72bvvUacnSwePfeXPGu9c1951CR0NO97fm5xH9mz+918zcn4c+bGwexZ653Hy73X1iP0u7UUKQB4qHGC9OzJTZKduSCS1+fJ6bzcJXsTuqbdfT5SgBknufGMZP8zmZPYPXcd7zEn49v2xsnF3rt987qOTijn4t3RBD7x4/KPuuv+EYztcqRd915b2nXkyHXokrdc76PfayaLz+93Mb/vH9XOb7Xeo8XYW/S7tRQpAHiocVB/xqO/s3FyvxUkxonhkcTlniJFzIWGvZPLbGdcLo4UOPhnY4J5pFB1jzFpTkFi/hjI3snb3B+OFini7MRgTsyznhb3WZ1pm/lpg7h2bZmvQ+m/e+8+X4t71nrn696R9e5hsvg2/W5u59wA2HvNuna+n7XeuRh7ZL176HdrKVIA8FBvPagnKRm3OU/wjkz27y1SJKEan+yIJPDXEq1sY17mGZ+v/WzGBDPt+8hJ0yXjpCATiTGJPjLZn/vwmSLFvI64VbBrk5q3KDS+J2MBNK5dI9Ln5o/e7Fku16z5mpDt3rqWbefv0rqftd523duz3sg6b40TJotv0+9i3k7O661l8vvEXTsvz1rvPGbvWW+kb956WlG/W0uRAoCHmpORDOx7ktV7jNvMZHD7/9GPm9xbpIgsM64jtgnr2A6Jm+9ARo5lXB/nZHI9t+uzJ9yZHIzbHPvirQLBaO5DZ4oUMSfZ27rSDlvRJv/m57nvh2LZP5vbdJtEje/tFE7zWn43xo3L3bq2zHeJI+vIORnvFme7iR37Wv6/nd/Zs9bbrntZb65x87FmO3k/jNfta2PEXDx79scIX9Fb9busY14mch7H6+el68bYh0bPWm/i23rTR7Kesb9ufXos4OTncX0j/W4tRQoAHqolq3skYU0SkMTgUiJ8SUu841oC0syJ0a2E7pIkRy1xuiXHf/TYuaxNvG/JectymbSdOf+t8BTXJmGz+T2U/Wlxe8yTm72yXFvfZ5f35zhx3yN9KdeE8bU9fevSde2W7N+lSV08a72XJqJ75H3T1hlZb1vmM/XRt+x3Z89jlrk25j5rvWmbseB1RNq0rTP0u7UUKQB4uDkxOupWUtK05OfohP9RRYrIpHTvJDn7fuROO/vk/J8pVIyy/JF+1BLbo0/HPLJIEVnf3glO4u7p959B3tt7JkU5b2Nbjr/b28bpT0f6cCame/rrs9abmEuFuib9bU9huu3rZ7uz/Zb97uh5zDb3FGKftd6YnyK5JnHp07fWrd+to0gBwEMlCckgvg3oSarO3uVIsaNto0nCMS57JpGYE5JHTNaSBGXf5nUnOc8+phhzK0HnnPGJlvybc7A3iR0dLTLM/f1oASr9blw++93ijsp6M0GY9y8/5/Uj7zf+/kRC3sNjn8r7Om35iGvHKEWFa9eRM0+gxbPWm+telm3vubx2tL9tk9txXWeu8R/BW/a7S+cx/89rKQzsLSKMnrXe9JOtfdIm23oj17m8fmTM1e/WUaQA4GEyoI8ToCQDLW6WZCSx8+QpSUaLb5LwjMuacH1u6U9bX0i/2puUJsmfC16xty9Hkthx2b3bBgAUKQB4oHFyd+ZuQyZz4x2LaHHNODHMOloMn0OKXmM/yt3iFndN7uRty8eRzyG76wYA5ylSAPAQmRiOk7qzj53mUc9xPS1mluLGuEwez2xxfA45/1tfSH9qMbek/459am+RIk/wjMsdeQIDAFCkAOBB5kndmc+Txvw50hYzm+96n7lzzscxFrrOFqzGj4vE3iLFuG1P9ADAcYoUADzEXKQ4cwd5/l6JPRPM+dH+s3fO+TjGQsGR7zXZ5MmcuVi2p+g2P0Vx5CMiAMDfKVIA8BDzxz1SODjykY/5aYi4NTHMZHL+ss2zHzPh4xg/7hH5ee+XV+YpnLlP7SmWZbmxWJb/790mAPCvKFIA8DDzXzWI3NXOUxXzRzAygUtBIV+2OU7uNu1JjLEAkbvW82Qy6xrj+ZzS18Z+Eelj6R/pQ3PxIPHpb+MTGJv2l0Hy89af8/8sO/dhf10GAM5RpADgYTJhmwsHR2Wy1yZ485Mas2x3XobPK4WD1k+OSNFtLlDsWbdiGQCcp0gBwENlUpePbrSnI65JfJZrk8K4NjE88jg/n0eemmhPR9ySZa59bGj+OMko36vSlgEA9lGkAOApUjRIYSETujZRTFEir+eu855H4zP5G7/MMP/Puq9NJiHy0YwUwNLf5i/EjDyFk6cm0sf2fEFm+t1YhMvy6cd7lgUArlOkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAAAMBLUKQAAAAAXoIiBQAAAPASFCkAAACAl6BIAQAAALwERQoAAADgJShSAAAAAC9BkQIAAAB4CYoUAAAAwEtQpAAAAABegiIFAAAA8BIUKQAAAICXoEgBAAAAvARFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAAAMBLUKQAAAAAXoIiBQAAAPASFCkAAACAl6BIAQAAALwERQoAAADgJShSAAAAAC9BkQIAAAB4CYoUAAAAwEtQpAAAAABegiIFAAAA8BIUKQAAAICXoEgBAAAAvARFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAv4D/69//tv/6Df+tf/+t/+F/5L3xI7ZgBgNfXxvWP4H/y3/wv//Uf/rv/oh4zaylSACz2H/+P/u06eH4k7bgBgNfXxvWPRKHi9ShSACz287/5bR00P5J23ADA62vj+kfy7/0b/1o9btZRpABYrA2YH007bgDg9bVx/aNpx806ihQAi7XB8qNpxw0AvL42rn807bhZR5ECYLE2WH407bgBgNfXxvWPph036yhSACzWBstosQAAnNPyrWixrKNIAbBYGyyjxQIAcE7Lt6LFso4iBcBibbCMFgsAwDkt34oWyzqKFACLtcEyWiwAAOe0fCtaLOsoUgAs1gbLaLEAAJzT8q1osayjSAGwWBsso8UCAHBOy7eixbKOIgXAYm2wjBYLAMA5Ld+KFss6ihQAi7XBMlosAADntHwrWizrKFIALNYGy2ixAACc0/KtaLGso0gBsFgbLKPFAgBwTsu3osWyjiIFwGJtsIwWCwDAOS3fihbLOooUAIu1wTJaLAAA57R8K1os6yhSACzWBstosQAAnNPyrWixrKNIAbBYGyyjxQIAcE7Lt6LFso4iBcBibbCMFgsAwDkt34oWyzqKFACLtcEyWiwAAOe0fCtaLOsoUgAs1gbLaLEAAJzT8q1osayjSAGwWBsso8UCAHBOy7eixbKOIgXAYm2wjBYLAMA5Ld+KFss6ihQAi7XBMlosAADntHwrWizrKFIALNYGy2ixAACc0/KtaLGso0gBsFgbLKPFAgBwTsu3osWyjiIFwGJtsIwWCwDAOS3fihbLOooUAIu1wTJaLAAA57R8K1os6yhSACzWBstosQAAnNPyrWixrKNIAbBYGyyjxQIAcE7Lt6LFso4iBcBibbCMFgsAwDkt34oWyzqKFACLtcEyWiwAAOe0fCtaLOsoUgAs1gbLaLEAAJzT8q1osayjSAGwWBsso8UCAHBOy7eixbKOIgXAYm2wjBYLAMA5Ld+KFss6ihQAi7XBMlosAADntHwrWizrKFIALNYGy2ixAACc0/KtaLGso0gBsFgbLKPFAgBwTsu3osWyjiIFwGJtsIwWCwDAOS3fihbLOooUAIu1wTJaLAAA57R8K1os6yhSACzWBstosQAAnNPyrWixrKNIAbBYGyyjxQIAcE7Lt6LFso4iBcBibbCMFgsAwDkt34oWyzqKFACLtcEyWiwAAOe0fCtaLOsoUgAs1gbLaLEAAJzT8q1osayjSAGwWBsso8UCAHBOy7eixbKOIgXAYm2wjBYLAMA5Ld+KFss6ihQAi7XBMlosAADntHwrWizrKFIALNYGy2ixAACc0/KtaLGso0gBsFgbLKPFAgBwTsu3osWyjiIFwGJtsIwWCwDAOS3fihbLOooUAIu1wTJaLAAA57R8K1os6yhSACzWBstosQAAnNPyrWixrKNIAbBYGyyjxQIAcE7Lt6LFso4iBcBibbCMFgsAwDkt34oWyzqKFACLtcEyWizn/fzzz3999dVX//Dtt9/WuNlvv/32xXLffPPNX3/++WeN5TLtv5b25zPT/9m0fCtaLOsoUgAs1gbLaLHc57vvvvsi4Uzi2uJG8zJJWlsct2n/tbQ/n5n+T7R8K1os6yhSACzWBstosdzn999//yLh/Prrr//6448/amz8+uuvX8R///33NY59tP9a2p/PTP8nWr4VLZZ1FCkAFmuDZbRY7vfTTz/tTjzzaO8Wl4TWY7730/5raX8+M/2flm9Fi2UdRQqAxdpgGS2Wx8jnkcdEtT3CO3+Gec+jweyj/dfS/nxm+v/n1vKtaLGso0gBsFgbLKPF8hhJSscENHfMxt/njlnunG2/3/sla+yj/dfS/nxm+v/n1vKtaLGso0gBsFgbLKPF8jg//vjjF4lqHgPefjc/EpzPMo/Lcj/tv5b25zPT/z+vlm9Fi2UdRQqAxdpgGS2Wx8ndsvkzx/kStRgT1CSzbXnuo/3X0v58Zvr/59XyrWixrKNIAbBYGyyjxfJY7dvbf/jhh3/8nCQ2yWxblvtp/7W0P5+Z/v85tXwrWizrKFIALNYGy2ixPF4S0zFRHSWJbcvwONp/Le3PZ6b/fz4t34oWyzqKFACLtcEyWiyPlztl45ekbb777rsaz2Np/7W0P5+Z/v/5tHwrWizrKFIALNYGy2ixPN6lJNVdtLeh/dfS/nxm+v/n0/KtaLGso0gBsFgbLKPF8niXHvfN55FbPI+l/dfS/nxm+v/n0/KtaLGso0gBsFgbLKPF8li5WzYmpnnEd/x5/LN0PJ72X0v785np/59Ty7eixbKOIgXAYm2wjBbL4+Qx39wtGxPUvP7tt99+kaj6G/nPof3X0v58Zvr/59XyrWixrKNIAbBYGyyjxfI4+fv3YzKav4+f13/77bcvXk/SOi/L/bT/Wtqfz0z//7xavhUtlnUUKQAWa4NltFgeI3fHxkQ0fxt//P34t/LDY7+Ppf3X0v58Zvr/59byrWixrKNIAbBYGyyjxfIY4yO9+Wb3PPo7/j4/z9/47rHfx9H+a2n/9fJ9CNvHDfJRg6Ptm/ht+ZzPM8tv/SDr+UznV///3Fq+FS2WdRQpABZrg2W0WO6Xu2Jj8nnpLtnPP//8Rdz2mWXuo/3X0v7r5aMFY9tGJs4ttmmT6DbZvuTe5d8z/Z+Wb0WLZR1FCoDF2mAZLZb7ZHIwJue5g9jiNklMx0Q1iWuLYx/tv5b2fw3zRHmz9279/L0Jm7ze4meXls/THS3+o9D/iZZvRYtlHUUKgMXaYBktlvvMSeetpH6+45kEN6+1WG7T/mtp/9dwb5Fi/tOZm3uLFL/88kuN/yj0f6LlW9FiWUeRAmCxNlhGi+W8JOBjwvn999/XuNk8oUii2+K4Tvuvpf1fxzz5DR/3eC79n03Lt6LFso4iBcBibbCMFss5c1Ke/x+5IzZ+0Vp47PcY7b+W9n89vjjz7ej/jFq+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktlvN+/vnnv7766qt/+Pbbb2vc7LfffvtiuW+++eavP//8s8YCAK+r5VvRYllHkQJgsTZYRovlPt99990XBYcULlrcaF4mRYsWBwC8tpZvRYtlHUUKgMXaYBktlvv8/vvvXxQcvv7667/++OOPGhu//vrrF/Hff/99jQMAXl/Lt6LFso4iBcBibbCMFsv9fvrpp92Fh3y0Y4tLQcPHPADg/Wr5VrRY1lGkAFisDZbRYnmMfB/FWKhoH+GYv8Niz0dDAIDX1fKtaLGso0gBsFgbLKPF8hjtyzDH3+eJiTw5sf1+75dsAgCvq+Vb0WJZR5ECYLE2WEaL5XF+/PHHLwoV+RjI9rv5IyH5LotxWQDg/Wn5VrRY1lGkAFisDZbRYnmcPC0xf+dEvkQzxgJFihlteQDgfWn5VrRY1lGkAFisDZbRYnms9tc7fvjhh3/8nCKGL8sEgI+h5VvRYllHkQJgsTZYRovl8VKYGAsVoxQx2jIAwPvT8q1osayjSAGwWBsso8XyePOXZG6+++67Gg8AvE8t34oWyzqKFACLtcEyWiyPd6lI4SkKAPhYWr4VLZZ1FCkAFmuDZbRYHu/Sxz3mP0sKALxvLd+KFss6ihQAi7XBMlosjzV/cWY+4jH+PP5ZUgDgfWv5VrRY1lGkAFisDZbRYnmc+U+Qbt9B8e23335RqPj999//aVkA4P1p+Va0WNZRpABYrA2W0WJ5nB9//PGLYsQff/zxt9d/++23L15P0WJeFgB4f1q+FS2WdRQpABZrg2W0WB4jT0eMhYgffvjhi9/n5/H3PvYBAO9fy7eixbKOIgXAYm2wjBbLY4wf6chf9shHP8bft7/44WMfAPC+tXwrWizrKFIALNYGy2ix3C9PRYzFh0tPSfz8889fxG3fWQEAvE8t34oWyzqKFACLtcEyWiz3yfdOjE9I3Pozo/Nf+0jhosUBAK+v5VvRYllHkQJgsTZYRovlPnPRIV+S2eI2KWqM8SlwbF+wCQC8Ly3fihbLOooUAIu1wTJaLOf98ssvXxQcvv/++xo3mz8e4mMfAPA+tXwrWizrKFIALNYGy2ixnDN/EebRJyLGL9oMH/sAgPen5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFjO+/nnn//66quv/uHbb7+tcbPffvvti+W++eabv/78888aCwC8rpZvRYtlHUUKgMXaYBktlvt89913XxQcUrhocaN5mRQtWhwA8NpavhUtlnUUKQAWa4NltFju8/vvv39RcPj666//+uOPP2ps/Prrr1/Ef//99zUOAHh9Ld+KFss6ihQAi7XBMlos9/vpp592Fx7y0Y4tLgUNH/MAgPer5VvRYllHkQJgsTZYRovlMfJ9FGOhon2EY/4Oiz0fDQEAXlfLt6LFso4iBcBibbCMFstjtC/DHH+fJyby5MT2+71fsgkAvK6Wb0WLZR1FCoDF2mAZLZbH+fHHH78oVORjINvv5o+E5LssxmUBgPen5VvRYllHkQJgsTZYRovlcfK0xPydE/kSzRgLFClmtOUBgPel5VvRYllHkQJgsTZYRovlsdpf7/jhhx/+8XOKGL4sEwA+hpZvRYtlHUUKgMXaYBktlsdLYWIsVIxSxGjLAADvT8u3osWyjiIFwGJtsIwWy+PNX5K5+e6772o8APA+tXwrWizrKFIALNYGy2ixPN6lIoWnKADgY2n5VrRY1lGkAFisDZbRYnm8Sx/3mP8sKQDwvrV8K1os6yhSACzWBstosTzW/MWZ+YjH+PP4Z0kBgPet5VvRYllHkQJgsTZYRovlceY/Qbp9B8W33377RaHi999//6dlAYD3p+Vb0WJZR5ECYLE2WEaL5XF+/PHHL4oRf/zxx99e/+233754PUWLeVkA4P1p+Va0WNZRpABYrA2W0WJ5jDwdMRYifvjhhy9+n5/H3/vYBwC8fy3fihbLOooUAIu1wTJaLI8xfqQjf9kjH/0Yf9/+4oePfQDA+9byrWixrKNIAbBYGyyjxXK/PBUxFh8uPSXx888/fxG3fWcFAPA+tXwrWizrKFIALNYGy2ix3CffOzE+IXHrz4zOf+0jhYsWBwC8vpZvRYtlHUUKgMXaYBktlvvMRYd8SWaL26SoMcanwLF9wSYA8L60fCtaLOsoUgAs1gbLaLGc98svv3xRcPj+++9r3Gz+eIiPfQDA+9TyrWixrKNIAbBYGyyjxXLO/EWYR5+IGL9oM3zsAwDen5ZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5GCl/fdd9/99dVXX/3Db7/9VuNYK+flhx9++Oubb775x7n6+uuv/3b+fv7557oMf9cGy2ixnJd+OF5Lvv322xo3S98el0sf//PPP2ssAPC6Wr4VLZZ1FCl4eYoUr+/HH3/84hzN9k4GP6s2WEaL5T7z9WRPAc01CAA+hpZvRYtlHUUKXp4Jwmv79ddfvzg/TYoYbVn+rg2W0WK5z++///5F38zTPn/88UeNjbl/f//99zUOAHh9Ld+KFss6ihQ3zI/5xtlJ8k8//fSQ9Xw2ihSvbfx4R2QSt0368kh8JnnXJoEYMN/afC2+VniYP77kYx4A8H61fCtaLOsoUtzQihRnH11XpDhHkeJ1ze+PTOhaHNe1wTJaLI+R6/jYd9t1Zf4OC9+tAgDvW8u3osWyjiLFDa1IcTZZVaQ4R5Hidc19Oj+3OK5rg2W0WB7jVoEtT0zkyYnt975XBQDev5ZvRYtlHUWKGy4VKc489qtIcY4ixeuavzAzH+1ocVzXBstosTzO3H/HItt8vc53WYzLAgDvT8u3osWyjiLFDZeKFHH0ywAVKc5RpHhdzs1jtMEyWiyPk0Lz/J0T+f6UGPu1L34FgI+h5VvRYllHkeKGuUgx33k7cndNkeIcE+HX5dw8Rhsso8XyWO2vd/zwww//+DlFDF+WCQAfQ8u3osWyjiLFDXORIgnt+DnlTNLaco0ixTnzRNhj169DkeIx2mAZLZbHS2Fi7McjH2ECgI+j5VvRYllHkeKGuUiRn+dvfN+bxL7XIsUvv/zyT0l8Jqd5vcVvcnxZbi7qpB2O3JmcJ8It5r1Le6RdxkfPc9yv/qc7FSkeow2W0WJ5vLz/xuvUJv27xQMA71PLt6LFso4ixQ2tSJHXx8nk3seBjxYpxthoMc28z7cS7XG/xtgc0zwJneUb7+djz8/X7kxGJgR7n4iY92H8XSbx+QjO/OcEs0xev3eSn31s68/PeST8yKT8UjtnG2N/Gj2qSLG109yW2W7OVQpOe/pwZD2X9rc50kafVRsso8XyeOn7rUjhKQoA+FhavhUtlnUUKW6YJ/yZaOb1TOra69e8hyLF9mf2MnFuSXszrj/J/jyhvyTr3zMxnifWeS3LjZ8bv2bPuZllUj9v95LE7SkmjO28tVmWu9TOeX1ex1FHjiPb2/Ondfeub6NIcVsbLKPF8niXiqopxrV4AOB9avlWtFjWUaS44VKRIsa7yZng3ZqovociRYwFivyb32/FhBzjHB/bnfixQDE+yXBpuRQaxn1p5knx/L0ge+zZzuZIgWaT+CzX1rcZjz/xee1aQefWebsl/eDocUQmbNeKR4oUj9cGy2ixPFauJ2N/nft33rdtOQDg/Wn5VrRY1lGkuOFakWL+XSZ347KzeZJ+a/I2xkaLaeb9Olqk2Ca2mUBfKrzMy2Qb24Q7y1/6vor5CZRtsn7NtUlxfpdJxjapzr/5uS2Tosm87lkrUKTAMW4j0sbzkxwpWl2b3M9tNv6cZcePXKTdbxU9rsmy47bGbYzntB1H3Oozm7mdb/VpujZYRovlcfJ+G4vNW7+fi4f3vBcBgNfR8q1osayjSHHDPOHPxHL8/ZFJ2jxJvTWhG2OjxTTzPt+acM77FUnSr02449Jd+ksFis04KYhb7TC3cWQdt5Zrk+9LRZfI8Y7HlP/f2sZcdJn7x2hu521bt55cOGo+jkhbXNtGjnNe5tqxbBQpHqMNltFieZwULsf+u10f0o/H13M9nJcFAN6flm9Fi2UdRYob5mR1nrglqR1/fy2ZnSeptyZ0Y2y0mGbe56NFikxW90ya2+e490xs54nBrWXmiXB+3jupn5fNZL3Fxbxfe++eju2XtmsxMbdz7CkGHTUXZ/ZOsOZ+E9eKOqFI8RhtsIwWy2PMTxvN14b5fbTn2gYAvLaWb0WLZR1FihtuFSlinnxe+vLBOe5VixR7k/F5uVvb2cx/wvXW9u6ZCM9tcamIkELBGHetmDGbl81HQ1rc3F7x6MfI56JZ3Co0jOaJ2a12UKR4jDZYRovlMVK82/ptK8zm5/npIh/7AID3reVb0WJZR5HihnmS2ybUczLbEt6YJ6m3JnRjbLSYZt7ntypS7F1u3r9b3+Vx70R4nmi05efCyZGJfYwTnkvtMLfXkULIXvNx3Grb2Xx3OW3X4jaKFI/RBstosdxv77Vrfj/dupYCAK+t5VvRYllHkeKGeUJ9KZmdv5ugfUnjnBjfmtCNsdFimnmfX71IcWv/7p0Iz8u3J13GmHzfxfz7W8blLxUGzrbXEfOx3vp+kObId4bce274uzZYRovlPilAjoXLW+/3PdcPAOB9aPlWtFjWUaS4YZ5QX5tYzsns/GjwPEm9NaEbY6PFNPM+Z79a3Obs5Pnsckf3796J8J79nJ+2uMel4znbXkfseWrklvm7Rq5Nyu49N/xdGyyjxXKfo312/ghV3mNHn7QCAF5Dy7eixbKOIsUN84T62sRyjp0nq/Mk9VZyPMZGi2lu7cfs7OT57HJH9+/eifCe/Rx/f69Lx3O2vY4Y1x8t5pYj+3nvueHv2mAZLZbz5ife9n4can5P3LpmAQCvqeVb0WJZR5HihnlCfWtiOX/x4Pi4/Zzo3prQjbHRYpqjRYAjk9LR2eWO7t+9E+E9+zn+/l4rP+4xrj9azC1H9lOR4jHaYBktlnPadwcdeSJi/N6Z8LEPAHh/Wr4VLZZ1FClumCfUtyaWcyI8ft55nvzdmtCNsdFimqNFgCOT0tHZ5Y7u31sXKc58J8UeZ9vriHH90WJuObKfihSP0QbLaLEAAJzT8q1osayjSHHDPKHeM7GcJ3nb0xT505Tj67cmdGNstJjmaBHgyKR0dHa5o/t370R4z3cszN/lMP/+Ec621xHzcZz5k4lze13bT0WKx2iDZbRYAADOaflWtFjWUaS4YZ5Q751Yjn8hIf/PExbzum5N6MbYaDHNvJ3PXqQYz8Wl5edtpKA0x9zrbHsdMR+Hv+7xPrTBMlosAADntHwrWizrKFLcME+oz07Es9z82q0J3fwZ6L13xeft3CoCnJ08n13u6P7dMxFOm43L5kmDFpenK8a4fLdIi7vH2fY6Yj6OvV8MuNnbXhtFisdog2W0WAAAzmn5VrRY1lGkuGGeUB+ZWI6PzWeyd/TjHvMEcO/d/flx/c9cpJi/yPRS8WH+M4Nx5qMS15xtryPacRz5csC97bVRpHiMNlhGiwUA4JyWb0WLZR1FihvuKVLME8Z5knprQjdPGH/88ccaN5qXiY9WpNj7rfrzduLahH1uuzzJko/ptNgzzrbXUfNx3GrfzdH2CkWKx2iDZbRYAADOaflWtFjWUaS4YZ64HZ1YzhPT0a0J3fw3/fM0xpFJ9uajFSkix3qtgJBtzF8ieavIM/9llkihYs+TCNnereM/215HteN4RnuFIsVjtMEyWiwAAOe0fCtaLOsoUtwwT6iPTiwzMZy/iHCzZ0I3TxyzrhQvtgln/s3P4zbmnz9ikSLSNtnm+LGMrLsVa/Y+FTHvW2Q7Wed8vrLdPNUxfnfItYLG2fY6ox3H1nfGfbzWXuP6LlGkeIw2WEaLBQDgnJZvRYtlHUWKG+bJ3pmJ5fxExGbPhO7SspckPsuNk8ePVKTI923ME+Nbjn5sI8WHuTi0Vyb8bZ1xtr3OyneYnDmOtPHe9prPhSLFOW2wjBYLAMA5Ld+KFss6ihQ3PKJIEW1ivXdCl8fu52VnuUs+PlEwfnnmrSLA2cnz2eXuKVJs28gTDHsm4Gm7IwWKTZZpTxhckvbPPl3b1tn2ukeemmh9r0l77v2+j828bkWKc9pgGS0WAIBzWr4VLZZ1FClueFSRIgWEcT1xZEKX2Eyax49xZFKZSeL29MRZZyfPZ5c7WqS4JAWBHHsKMmPBIutLcWLPd0ncknVk4p51zkWRvJZzsvevrpxtr0fIcaRNss/jPqQ/pf3SjmeKOfP6FCnOaYNltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpADelf/3/+X/9Nd/8j/7d/76n/+3/2sfRhssox0/AADntHwr/sN/9198GMmT/9P/5/+jHv97oUgBvBv/t//4f1sHlo+qtQEAAOe0fOuj+r//H/6T2gbvgSIF8G78L/4739WL8EfV2gAAgHNavvVR/Xv/xr9W2+A9UKQA3o1/57/6X6wX4Y+qtQEAAOe0fOsja23wHihSAO/GL//d/3q9AH9UrQ0AADin5Vsf1f/qv//fqm3wHihSAO/G/+v//H/863/83/gv1QvxR9TaAACAc1q+9RHlox6+kwLgjfz//r//n7/+r/+7/81f/+t/8d/7MNrgEu34AQA4p+Vb8b//X/5PP4x80fx//p/9p/X43wtFCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCnghJ9++umvr7766h/yc4uDPdpgGS0WAIBzWr4VLZZ1FCngBEUKHqkNltFiAQA4p+Vb0WJZR5ECTlCk4JHaYBktFgCAc1q+FS2WdRQp4ARFCh6pDZbRYgEAOKflW9FiWUeRAk5QpOCR2mAZLRYAgHNavhUtlnUUKeAERQoeqQ2W0WIBADin5VvRYllHkQJOUKTgkdpgGS0WAIBzWr4VLZZ1FCngBEUKHqkNltFiAQA4p+Vb0WJZR5ECTlCk4JHaYBktFgCAc1q+FS2WdRQp4ARFCh6pDZbRYgEAOKflW9FiWUeRgof5888///r555//+vbbb/8xef/666//+v777//67bff6jKbX3755a/vvvvui+Xyc15v8aspUvBIbbCMFgsAwDkt34oWyzqKFC8gk/tMxjOZ/+abb76Y/Gai/uOPP/71xx9/1GWv+f333/82ec465vWmkPDDDz/cLB7MxkLCODHPtuZtzLK9cV3bcmNRo8nv00bzskdtbZwCyLbu7HNeO1oMeUSRIseeczsf/9lzE+knKRTlPI3HuRV98rtHtCWP1QbLaLEAAJzT8q1osayjSLHYr7/++sVk8ppMMvcUKzLhvjXxH2VCvHfiOhYpMsHOa9neuL5r5sLGkWPflrskk/pxme317N+tAkokJvs0rvOSe4oUOYdjO16z95zHvE/X5NwpVryONlhGiwUA4JyWb0WLZR1FioWOTO5Hl+6wp+CxZzLetKccmnFynf+Px5Bt5+ctNhP+PKWw/X6TSfdYoMi/4x3+S8uN627mIkXaY28xYJTl2vpHZ4sURwozm8TfKp7k/LVlr8l6FSpeQxsso8UCAHBOy7eixbKOIsUimajPk9VMzMfJaCaQ80Q7T0iM6xnNRY+sP3fMM3kfJ6OXJu97JufjcuP+X3saY95WYrdlrz29ME+80z4tbjMXKUbZXgoJ47by/7w2n4f8fOvphTNFimxv3laOMe0+tl2OYz72tNOl9s3yY2xkf8ZjyP/TP8Yi1vYkDOu1wTJaLAAA57R8K1os6yhSLJInB8ZJ5a0JYyauKVDcmjynILBNyC9Najfz0wrXCiCbVtzIhLrFbjI5n5eJ7OelAkXkWOdlWtzmUpHi1kcbLhUPWuzmaJEi2x+3kf9nf1vsZi46XdrGfB5vFZu276xov2ONNlhGiwUA4JyWb0WLZR1FikXmyf6t4sNemXDfKk5sEjfuw579mPd772R3LgLErUl6zN+tcW2ZuUiRbd76iMimFTiutcXRIkUKJWP8teLMaNxOjqfFjOvNkxIthtfWBstosTzG9p7M+8pTRQDwObR8K1os6yhSLDJP9vcWFh5t3o9bd+Hn+D2FhpiXuzWp38xPCRwpUuwtoGzmfcwTBy0ujhQp5mLQrac0RvOy7fyMv9/zNAyvpw2W0WK53/z+jb3XJADg/Wr5VrRY1lGkWGT+zoEjE9dHOjLZjkcVKfYuN+/ftcLBvUWK+SM415Y/0m7zeo8+NTM+TdK2M37PROx9eoTX0QbLaLHcb74exdHrBQDw/rR8K1os6yhSLNK+pyFPDTzqYx97HZlsx9liw9nljuzfvUWKeflLH6+II/s1HvuZj2OMy7cvD52LIJH9WfV0Dse1wTJaLPebP0YWR68XAMD70/KtaLGso0ix0PzFiJsk0Jl4PnKSmXVlEp7Ja57aSEI+34GPa5PtGCfM8ZGKFDEuHy0mjuxX+z6Osy4d0/xkziZFDU9WvL42WEaL5X7z+zeuvYcBgI+h5VvRYllHkWKxTCCvTWIzybz1PRGX5GmNfCFcu2t4ya1EXZHi747s1xh3r2vHNH855yh9LIWMvV/Yydtqg2W0WB7DF2cCwOfT8q1osayjSPEC8pRDJrnXihV56mFvsSKT9SOFidG1yXYoUvzdkf0a4+7VPu4xyseFLj1VsUm7KFa8ljZYRosFAOCclm9Fi2UdRYoXk0JEJpmXChb5XVtuc22CmsJF7hjmoySZ0GdCe2SyHYoUf3dkv8a4t/oToSl85SmdFDXG7Y98DOR1tMEyWiwAAOe0fCtaLOsoUryobZLZvjciRYa2TCtQZKKe9Vz6fosjk+34yEWK9mWmLS6O7NdccGoxz7QVo1rhyxMVr6ENltFiAQA4p+Vb0WJZR5HiHZiLD5lszjHzBD323Ck/MtmOj1ykSHvtXf7Ifs3HfvY7Ru6VQtX8MaBbHx/hbbTBMlosAADntHwrWizrKFK8E7cm+fOXJt4qNmyOTLbjbLHh7HJH9u/eIsX80Yhr2zqyX/OfCL31kZ1nSqFi9ZMd/LM2WEaLBQDgnJZvRYtlHUWKd2KeFM+T/LkIkEf8x99fcmSyHR+1SJH2GpeNa214ZL/auld+zGI+Fy2Gt9UGy2ixAACc0/KtaLGso0ixSCapl75bopk/8vGIIkXWMd9V/0hFiny0ocU18/7d+hjEkf2K+fxl3y59T8hR+ZjK3vaM+XtOWgxvqw2W0WIBADin5VvRYllHkWKB8bsBMjm+VVBIQWMsJrTvpDj6cY95kr13ufdUpIi087WnFnIu2nd+3DonR/Yr2scssm97i0mX1j/2jfSBW4WP+Xs3fCfFa2iDZbRYAADOaflWtFjWUaRYYJ4UR17LFyqOk8xMYPO0xTy5bRPWNkHPOscJetaXSep4J33PukfvrUixmdu3tcUmr8/rns3fM7Fnst/2L+2ffZvbY3vSZitmxVzQyLGMv9/Wl2JFW99cyIq954HnaoNltFgAAM5p+Va0WNZRpFggE8Z5crlXJrRtnTF/8eMt21Mc42sfqUhxpo33FCgi57Atf6v9stxcGNqrnfsUXs6ub++x8nxtsIwWCwDAOS3fihbLOooUC126k3/JrQlw7qxnItuWHWWb4wR1LCB8pCJFYlOEmbfdpE327tOmrXfPExV7z9Mm+5anKsanbOb15Vj3FisSt+rPoNK1wTJaLAAA57R8K1os6yhSvIBMGPMo/jzpzWQyr2WCOj/qf00m25kEjwWQrCsT6Hvvnp8tNpxd7t4ixfa7PMEwt0n+f0+bbMWGsTiwp0ixyTnNuU3bzAWGvJZ1HykmZH9yLFlufopkPNZLxQ7WaYNltFgAAM5p+Va0WNZRpODDuFakgFfWBstosQAAnNPyrWixrKNIwYcxFynydEqLg1fTBstosQAAnNPyrWixrKNIwYcxFynykYkWB6+mDZbRYgEAOKflW9FiWUeRgg9DkYL3qg2W0WIBADin5VvRYllHkYIPQ5GC96oNltFiAQA4p+Vb0WJZR5GCD0ORgveqDZbRYgEAOKflW9FiWUeRgg9DkYL3qg2W0WIBADin5VvRYllHkYIPQ5GC96oNltFiAQA4p+Vb0WJZR5GCD0ORgveqDZbRYgEAOKflW9FiWUeRgg9DkYL3qg2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WB7jxx9//Ourr7766+uvv/7b/1sMAPCxtHwrWizrKFIALNYGy2ix3O+nn376W4FilNdaLADwcbR8K1os6yhSACzWBstosdzvu++++6ciRV5rsQDAx9HyrWixrKNIAbBYGyyjxXK/b7/9VpECAD6hlm9Fi2UdRQqAxdpgGS2W+/m4BwB8Ti3fihbLOooUAIu1wTJaLI/hizMB4PNp+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFieYwff/zxr6+++uqvr7/++m//bzEAwMfS8q1osayjSAGwWBsso8Vyv59++ulvBYpRXmuxAMDH0fKtaLGso0gBsFgbLKPFcr/vvvvun4oUea3FAgAfR8u3osWyjiIFwGJtsIwWy/2+/fZbRQoA+IRavhUtlnUUKQAWa4NltFju5+MeAPA5tXwrWizrKFIALNYGy2ixPIYvzgSAz6flW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WM77+eef//rqq6/+4dtvv61xs99+++2L5b755pu//vzzzxoLALyulm9Fi2UdRQqAxdpgGS2W+3z33XdfFBxSuGhxo3mZFC1aHADw2lq+FS2WdRQpABZrg2W0WO7z+++/f1Fw+Prrr//6448/amz8+uuvX8R///33NQ4AeH0t34oWyzqKFACLtcEyWiz3++mnn3YXHvLRji0uBQ0f8wCA96vlW9FiWUeRAmCxNlhGi+Ux8n0UY6GifYRj/g6LPR8NAQBeV8u3osWyjiIFwGJtsIwWy2O0L8Mcf58nJvLkxPb7vV+yCQC8rpZvRYtlHUUKgMXaYBktlsf58ccfvyhU5GMg2+/mj4TkuyzGZQGA96flW9FiWUeRAmCxNlhGi+Vx8rTE/J0T+RLNGAsUKWa05QGA96XlW9FiWUeRAmCxNlhGi+Wx2l/v+OGHH/7xc4oYviwTAD6Glm9Fi2UdRQqAxdpgGS2Wx0thYixUjFLEaMsAAO9Py7eixbKOIgXAYm2wjBbL481fkrn57rvvajwA8D61fCtaLOsoUgAs1gbLaLE83qUihacoAOBjaflWtFjWUaQAWKwNltFiebxLH/eY/ywpAPC+tXwrWizrKFIALNYGy2ixPNb8xZn5iMf48/hnSQGA963lW9FiWUeRAmCxNlhGi+Vx5j9Bun0HxbfffvtFoeL333//p2UBgPen5VvRYllHkQJgsTZYRovlcX788ccvihF//PHH317/7bffvng9RYt5WQDg/Wn5VrRY1lGkAFisDZbRYnmMPB0xFiJ++OGHL36fn8ff+9gHALx/Ld+KFss6ihQAi7XBMlosjzF+pCN/2SMf/Rh/3/7ih499AMD71vKtaLGso0gBsFgbLKPFcr88FTEWHy49JfHzzz9/Ebd9ZwUA8D61fCtaLOsoUgAs1gbLaLHcJ987MT4hcevPjM5/7SOFixYHALy+lm9Fi2UdRQqAxdpgGS2W+8xFh3xJZovbpKgxxqfAsX3BJgDwvrR8K1os6yhSACzWBstosZz3yy+/fFFw+P7772vcbP54iI99AMD71PKtaLGso0gBsFgbLKPFcs78RZhHn4gYv2gzfOwDAN6flm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYBktFgCAc1q+FS2WdRQpABZrg2W0WAAAzmn5VrRY1lGkAFisDZbRYgEAOKflW9FiWUeRAmCxNlhGiwUA4JyWb0WLZR1FCoDF2mAZLRYAgHNavhUtlnUUKQAWa4NltFgAAM5p+Va0WNZRpABYrA2W0WIBADin5VvRYllHkQJgsTZYRosFAOCclm9Fi2UdRQqAxdpgGS0WAIBzWr4VLZZ1FCkAFmuDZbRYAADOaflWtFjWUaQAWKwNltFiAQA4p+Vb0WJZR5ECYLE2WEaLBQDgnJZvRYtlHUUKgMXaYPnRtOMGAF5fG9c/mnbcrKNIAbBYGyw/mnbcAMDra+P6R9OOm3UUKQAWa4PlR9OOGwB4fW1c/2jacbOOIgXAYj//m9/WAfMjaccNALy+Nq5/JP/ev/Gv1eNmHUUKgMX+o3//f1AHzY+kHTcA8PrauP6R/If/7r+ox806ihQAL+A/+vf/7b/+g3/rX6+D50fQjhkAeH1tXP8I8gRFChT/+X/2n9bjZh1FCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwEhQpAAAAgJegSAEAAAC8BEUKAAAA4CUoUgAAAAAvQZECAAAAeAmKFAAAAMBLUKQAAAAAXoIiBQAAAPASFCkAAACAl6BIAQAAALwERQoAAADgJShSAAAAAC9BkQIAAAB4CYoUAAAAwEtQpAAAAABegiIFAAAA8BIUKQAAAICXoEgBAAAAvARFCgAAAOAlKFIAAAAAL0GRAgAAAHgJihQAAADAS1CkAAAAAF6CIgUAAADwAv7lX/9/soCIHgSrOFAAAAAASUVORK5CYII=)"
      ],
      "metadata": {
        "id": "qqWn6Ayl7hUz"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# You can edit this first variable to choose which activity data you want to run (walk, run, bike, stairs)\n",
        "# Options for dataset_name in the line below: [walk_sample.csv, run_sample.csv, stairs_sample.csv, bike_sample.csv]\n",
        "dataset_name = 'walk_sample.csv'\n",
        "data_path = drive_path+'energyExpenditureData/' + dataset_name\n",
        "raw_IMU_data = np.loadtxt(data_path, delimiter=',')\n",
        "print(raw_IMU_data.shape)"
      ],
      "metadata": {
        "id": "cs5H1YH3NHKn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now we have loaded the CSV file containing the walking data sample. \n",
        "\n",
        "We have also printed out the size of the matrix containing the data. The size indicates that there are 13 time-series signals (columns) with 35,500 samples (rows) for each signal. These signals, in order, are:\n",
        "- The three signals from the gyroscope on the shank (x,y,z axes signals)\n",
        "- The three signals from the gyroscope on the thigh (x,y,z axes signals)\n",
        "- The three signals from the accelerometer on the thigh (x,y,z axes signals) \n",
        "- The three signals from the accelerometer on the shank (x,y,z axes signals)\n",
        "- One signal to indicate corrupted messages \n",
        "\n",
        "Since we sampled data at 100 Hz (or 100 samples/sec), that means this file is approximately 355 seconds or about 6 minutes. \n",
        "\n",
        "If you would like to use sensor data from a different sensor, format the samples and signals into the same rows and columns format defined here to be able to use the following processing code."
      ],
      "metadata": {
        "id": "iRwsQrItNHaX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Understanding IMU Signals During Lower-limb Movement\n",
        "\n",
        "Now we will plot the different signals to visualize the raw data. "
      ],
      "metadata": {
        "id": "20HT8mkq9GoW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plot_start = 5000 # Sample number to start plotting from\n",
        "plot_length = 500 # Number of samples to plot (seconds * 100)\n",
        "\n",
        "plt.figure() \n",
        "plt.title('Gyroscope signals on the shank')\n",
        "plt.plot(raw_IMU_data[plot_start:plot_start+plot_length,0:3])\n",
        "plt.xlabel('Samples (100 Hz)')\n",
        "plt.ylabel('Rate of rotation (deg/s)')\n",
        "plt.legend(['x','y','z'])\n",
        "plt.show()\n",
        "\n",
        "plt.figure() \n",
        "plt.title('Gyroscope signals on the thigh')\n",
        "plt.plot(raw_IMU_data[plot_start:plot_start+plot_length,3:6])\n",
        "plt.xlabel('Samples (100 Hz)')\n",
        "plt.ylabel('Rate of rotation (deg/s)')\n",
        "plt.legend(['x','y','z'])\n",
        "plt.show()\n",
        "\n",
        "plt.figure() \n",
        "plt.title('Accelerometer signals on the thigh')\n",
        "plt.plot(raw_IMU_data[plot_start:plot_start+plot_length,6:9])\n",
        "plt.xlabel('Samples (100 Hz)')\n",
        "plt.ylabel('Accelerations (m/s^2)')\n",
        "plt.legend(['x','y','z'])\n",
        "plt.show()\n",
        "\n",
        "plt.figure() \n",
        "plt.title('Accelerometer signals on the shank')\n",
        "plt.plot(raw_IMU_data[plot_start:plot_start+plot_length,9:12])\n",
        "plt.xlabel('Samples (100 Hz)')\n",
        "plt.ylabel('Accelerations (m/s^2)')\n",
        "plt.legend(['x','y','z'])\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "DL3F5d_WD7df"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Notice that the signals are cyclic, repeating approximately every 120 samples. This is because as the person walks, their leg moves in the same way with each step. Motion in the sagittal plane is associated with the z-axis measurements. We can see that these signals have the largest rates of rotation. The vertical axis measurements of acceleration (y-axis) show large spikes when the leg comes in contact with the ground."
      ],
      "metadata": {
        "id": "ShHpfOoOtFVc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Preprocessing Movement Data\n",
        "\n",
        "Next we will look at steps to preprocess our IMU data to extract the essential information about the movements. This will include filtering the signals to remove noise and segmenting the periodic motions into strides, also called gait cycles."
      ],
      "metadata": {
        "id": "vtDDm_QBglnM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Filter Lower-Limb Movement Data\n",
        "\n",
        "You'll notice that the raw IMU signals do not appear very smooth. Most have small perturbations, which is often called sensor noise. This noise is due to many factors including the sensitivity of the sensors (even resting sensors have noise) and high frequency components of motion, such as slight vibration or shifting of the sensor when in contact with the person due to shifting clothing or skin.\n",
        "\n",
        "We will use a low-pass filter to eliminate these higher frequencies of noise and isolate only a portion of lower-frequency signal associated with movement. Here, we use a common low-pass Butterworth filter to isolate a signal of 6 Hz or less, which is standard for filtering movement data for walking or running. This cutoff works well for our tasks of walking, running, biking, and stair climbing. You may want to pick a higher value for higher speed tasks like sprinting."
      ],
      "metadata": {
        "id": "kifafd1e9q_m"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Design the filter\n",
        "wn = 6 # Crossover frequency for low-pass filter (Hz)\n",
        "filt_order = 4 # Fourth order filter\n",
        "sample_freq = 100 # in Hz *Change this value if you use a different sample frequency*\n",
        "b,a = signal.butter(filt_order, wn, fs=sample_freq) # params for low-pass filter\n",
        "\n",
        "# Filter data\n",
        "filtered_data = signal.filtfilt(b,a,raw_IMU_data[:,:-1], axis=0)\n",
        "\n",
        "# Plot filtered and raw gyroscope shank z-axis signal\n",
        "shank_gyro_z_index = 2 # column number for z-axis of shank gyroscope readings used to segment gait cycle\n",
        "\n",
        "plt.figure() \n",
        "plt.title('Gyroscope shank z-axis signal')\n",
        "plt.plot(raw_IMU_data[plot_start:plot_start+plot_length,shank_gyro_z_index])\n",
        "plt.plot(filtered_data[plot_start:plot_start+plot_length,shank_gyro_z_index])\n",
        "plt.xlabel('Samples (100 Hz)')\n",
        "plt.ylabel('Rate of rotation (deg/s)')\n",
        "plt.legend(['raw','filtered'])\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "PM0iXYRO9tt_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "Notice in the plot above, the filtered signal (in orange) appears much smoother but retains the same general shape of the unfiltered signal (in blue)."
      ],
      "metadata": {
        "id": "WkvYwVnZwjt9"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Segmenting Cyclic Movement by Gait Cycle\n",
        "\n",
        "Next, we'll segment periodic movements by gait cycle. This is often done with a measure of ground reaction forces when the heel of the foot strikes the ground. \n",
        "\n",
        "With wearable sensors like IMUs, we typically do not have that data. We will instead use change in rotation of the leg right before heelstrike, specifically, the sagittal plane rotation of the shank IMU (the gyroscope z-axis signal shown in the plot above). Immediately before heel strike the leg stops swinging forward and starts to move backwards. This change in rotational direction creates a peak in the angular velocity of the shank.\n",
        "\n",
        "To detect the peak, we use a simple `find_peaks` function. The function depends on two parameters to automatically segment this movement:\n",
        "\n",
        "**Threshold selection:** We selected a threshold of 70 deg/s which reliably segmented gait cycles from slow walking at 0.7 m/s to running at 3.25 m/s and worked well for cycling and stair climbing. If you plan to use this approach for participants moving very slowly, you may need to adjust this threshold. \n",
        "\n",
        "**Minimum distance between peaks:** The minimum distance between peaks corresponds to the fastest expected gait cycle. Here we selected a minimum gait cycle of 0.6 seconds, which was a smaller value than the shortest gait cycle we recorded. If you plan to evaluate motion that is much faster or slower, you may need to adjust this value. \n",
        "\n",
        "Visual inspection was performed on motion data for all the activities we intended to estimate energy expenditure. If you would like to include additional activities, we recommend using trial and error to tune the threshold and minimum distance between peaks to determine what parameters work best for segmenting gait cycles."
      ],
      "metadata": {
        "id": "D-0fXeGr9lwP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "peak_height_thresh = 70 # minimum value of the shank IMU gyro-Z reading in deg/s\n",
        "peak_min_dist = int(0.6*sample_freq) # min number of samples between peaks\n",
        "\n",
        "# Find peaks\n",
        "peak_index_list = signal.find_peaks(filtered_data[:,shank_gyro_z_index], height=peak_height_thresh, distance=peak_min_dist)[0]\n",
        "peak_index_list_plot = signal.find_peaks(filtered_data[plot_start:plot_start+plot_length,shank_gyro_z_index], height=peak_height_thresh, distance=peak_min_dist)[0]\n",
        "\n",
        "# Plot\n",
        "plt.figure() \n",
        "plt.title('Gyroscope shank z-axis signal')\n",
        "plt.plot(filtered_data[plot_start:plot_start+plot_length,shank_gyro_z_index])\n",
        "axes = plt.gca()\n",
        "ymin, ymax = axes.get_ylim()\n",
        "xmin, xmax = axes.get_xlim()\n",
        "plt.plot([xmin,xmax],[peak_height_thresh, peak_height_thresh], color='k',ls='--')\n",
        "for i in peak_index_list_plot:\n",
        "    plt.plot([i,i],[ymin,ymax],color='k')\n",
        "plt.xlabel('Samples (100 Hz)')\n",
        "plt.ylabel('Rate of rotation (deg/s)')\n",
        "plt.legend(['filtered signal','threshold','detected heelstrikes'],loc='lower right')\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "T3QP__Pk9qsn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Format Data for Each Gait Cycle\n",
        "\n",
        "Now that we have segmented movement by gait cycle, we will format the data to pass into our estimation model. The duration of gait cycles changes depending on walking speed and activity, but our estimation model uses a fixed input size. We perform a time-invariant formatting to adjust varying lengths of gait cycles into one fixed size, by averaging the measurements across a gait cycle into a set number of discrete bins. \n",
        "\n",
        "In this case we use 30 bins to capture sufficient detail about the movement. Our experiments found that increasing the number of bins did not improve model performance for our application. Fewer bins resulted in worse performance, likely due to a loss of detailed information about movements.\n",
        "\n",
        "We also use a special condition where a gait cycle is ignored if it takes longer than the 4 second time window used to detect gait cycles. This was to prevent fringe cases where motion was not cyclic, such as someone pausing mid-step for a long time, or outside the activities we intended to capture. In this case, the slowest speed of stair walking had a gait cycle of about 2.4 seconds, so we were able to detect all gait cycles within 4 seconds.\n",
        "\n",
        "Run the following cell to process the time-series data into a time-invariant format. "
      ],
      "metadata": {
        "id": "4T-1yTgJ9uGG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Variables to change depending on your data\n",
        "weight = 68 # participant bodyweight in kg\n",
        "height = 1.74 # participant height in m\n",
        "num_bins = 30 # number of bins to discretize one gait cycle of data from each signal into a fixed size\n",
        "stride_detect_window = 4*sample_freq # maximum time length for strides (seconds * samples/second)\n",
        "\n",
        "# Initializing variables and useful parameters\n",
        "gait_cnt = 0 # keep track of the number of gait cycles\n",
        "gait_data = np.zeros((3+12*num_bins, len(peak_index_list)-1)) # storage for gait cycles of formatted data \n",
        "deg2rad = 0.0174533 # rad / deg\n",
        "\n",
        "# process time series data, format units, and discretize each gait cycle into a fixed number of discrete bins\n",
        "def processRawGait(data_array, start_ind, end_ind, b, a, weight, height, num_bins=30):\n",
        "    gait_data = data_array[start_ind:end_ind, :] # crop to the gait cycle\n",
        "    gait_data = gait_data*np.array([deg2rad,-deg2rad,-deg2rad,deg2rad,-deg2rad,-deg2rad,1,-1,-1,1,-1,-1]) # flip y & z, convert to rad/s\n",
        "    bin_gait = signal.resample(gait_data, num_bins, axis=0) # resample gait cycle of data into fixed number of bins along axis = 0 by default\n",
        "    shift_flip_bin_gait = bin_gait.transpose() # get in shape of [feats x bins] for correct flattening\n",
        "    model_input = shift_flip_bin_gait.flatten()\n",
        "    model_input = np.insert(model_input, 0, [1.0, weight, height]) # adding a 1 for the bias term at start\n",
        "    return model_input\n",
        "\n",
        "for i in range(len(peak_index_list)-1): # looping through each gait cycle of data\n",
        "    gait_start_index = peak_index_list[i] # index at start of gait\n",
        "    gait_stop_index = peak_index_list[i+1] # index at end of gait\n",
        "    if (gait_stop_index - gait_start_index) <= stride_detect_window: # if gait cycle within maximum time allowed\n",
        "        gait_data[:,gait_cnt] = processRawGait(filtered_data, gait_start_index, gait_stop_index, b, a, weight, height, num_bins)\n",
        "        gait_cnt += 1 # increment number of gait cycles of data stored\n",
        "gait_data = gait_data[:,:gait_cnt] # get rid of empty storage\n",
        "\n",
        "# Plotting discretized and continuous sampled data for one gait cycle\n",
        "gait_cycle_to_plot = 1\n",
        "plt.figure() \n",
        "plt.title('Discretized and sampled data for gyroscope shank z-axis signal')\n",
        "start_gait_index = peak_index_list[gait_cycle_to_plot]\n",
        "stop_gait_index = peak_index_list[gait_cycle_to_plot+1]\n",
        "plt.plot(filtered_data[start_gait_index:stop_gait_index,shank_gyro_z_index])\n",
        "\n",
        "# plotting discretized data\n",
        "gait_length = stop_gait_index - start_gait_index\n",
        "x = np.linspace(start=0,stop=gait_length,num=num_bins+1)\n",
        "bin = 2\n",
        "plt.plot(x[:-1], gait_data[3+num_bins*bin:3+num_bins*(bin+1),gait_cycle_to_plot]*(-1/deg2rad))\n",
        "plt.xticks(ticks=[0,gait_length],labels=[0,100])\n",
        "plt.xlabel('Gait cycle (%)')\n",
        "plt.ylabel('Rate of rotation (deg/s)')\n",
        "plt.legend(['sampled data', 'discretized data'])\n",
        "plt.show()\n",
        "\n",
        "print(\"Format for processed gait data:\",gait_data.shape)"
      ],
      "metadata": {
        "id": "jQGbRGFwuZ5C"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Notice the shape of the processed gait data is (363, X). The length of processed data for each gait cycle is 363 values. The first three are the model bias parameter, the subject weight, and the subject height. Then the next 360 values are comprised of the 30 bins recording one gait cycle of data for the 12 signals of the thigh and shank IMUs. The value for X is the number of gait cycles recorded during the condition."
      ],
      "metadata": {
        "id": "1tIcMsqbPjbY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Estimating Energy Expenditure\n",
        "\n",
        "We will now load in the model weights and pass the formatted data into the data-driven model to estimate the energy expended with each gait cycle."
      ],
      "metadata": {
        "id": "WiVPg349gtnE"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Load in model weights\n",
        "model_dir = drive_path+'energyExpenditureData/full_fold_aug/' # path to model weights\n",
        "model_weights = np.loadtxt(model_dir + 'weights.csv', delimiter=',') # model weight vector\n",
        "weights_unnorm = model_weights[3:] # get rid of bias and height/weight offset\n",
        "unnorm_weights = np.reshape(weights_unnorm,(-1,num_bins))\n",
        "\n",
        "# Apply model to data to estimate energy expended\n",
        "estimates = np.zeros(gait_cnt)\n",
        "for i in range(gait_cnt):\n",
        "    estimates[i] = np.dot(gait_data[:,i], model_weights)\n",
        "\n",
        "# Plot energy estimates\n",
        "plt.figure()\n",
        "plt.plot(estimates)\n",
        "plt.title(\"Energy expenditure per gait cycle\")\n",
        "plt.xlabel(\"Gait cycles\")\n",
        "plt.ylabel(\"Energy expenditure (W)\")\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "5mchOYxn-PxO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We will now compare the energy expenditure estimates to ground truth respirometry measurements collected for walking in the lab. Note that if you re-run with another activity of data this plot will still compare it to the respirometry recorded during walking."
      ],
      "metadata": {
        "id": "gKW3rGGuTzOg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "metabolics_breaths = np.loadtxt(drive_path+'energyExpenditureData/walking_metabolics.csv', delimiter=',', skiprows=1) \n",
        "avg_metabolics = np.mean(metabolics_breaths[len(metabolics_breaths)//2:])\n",
        "\n",
        "plt.figure()\n",
        "plt.plot(estimates)\n",
        "plt.title(\"Energy expenditure per gait cycle\")\n",
        "plt.xlabel(\"Gait cycles\")\n",
        "plt.ylabel(\"Energy expenditure (W)\")\n",
        "plt.plot(np.linspace(start=0,stop=len(estimates),num=len(metabolics_breaths)),metabolics_breaths, c='k', alpha=0.5)\n",
        "plt.plot([0,len(estimates)],[avg_metabolics, avg_metabolics], c='k')\n",
        "plt.legend(['estimates per gait cycle', 'respirometry per breath','steady-state respirometry'])\n",
        "avg_energy_expenditure_estimate = np.mean(estimates)\n",
        "plt.plot([0,len(estimates)],[avg_energy_expenditure_estimate, avg_energy_expenditure_estimate], c='b')\n",
        "plt.legend(['estimates per gait cycle', 'respirometry per breath','steady-state respirometry', 'steady-state estimates'])\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "Pj6y4X4tzPMk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "We see that the steady-state respirometry (in black) is quite close to the estimated energy expenditure from our model (in blue). You will also notice a slight delay in the respirometry. "
      ],
      "metadata": {
        "id": "pqHzrtqx48Eb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Intuition Behind the Data-Driven Model\n",
        "\n",
        "The Wearable System used a linear regression model to estimate energy expenditure with inputs from the IMU signals worn on the shank and thigh. A larger magnitude weight indicated a more informative input. The input signals are shown in descending order of importance based on contribution to total model weight. Gyroscope inputs were more informative than accelerations. The X, Y, and Z directions represent the fore-aft, mediolateral, and vertical axes, respectively. Run the following cell to visualize the model weights."
      ],
      "metadata": {
        "id": "aZRdjfLdgYT0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Plotting model weights\n",
        "pelv_x_dir = \" Z\"# dir\"\n",
        "pelv_y_dir = \" Y\"# dir\"\n",
        "pelv_z_dir = \" -X\"# dir\"\n",
        "thigh_x_dir = \" X\"# dir\"\n",
        "thigh_y_dir = \" -Z\"#-dir\"\n",
        "thigh_z_dir = \" Y\"#-dir\"\n",
        "shank_x_dir = \" X\"#-dir\"\n",
        "shank_y_dir = \" -Z\"#-dir\"\n",
        "shank_z_dir = \" Y\"#-dir\"\n",
        "thigh = \"Thigh\"\n",
        "shank = \"Shank\"\n",
        "accel = \" accel\"\n",
        "gryo = \" gyro\"\n",
        "\n",
        "full_features = ['bias','weight','height',\n",
        "                 shank+gryo+shank_x_dir,shank+gryo+shank_y_dir,shank+gryo+shank_z_dir,\n",
        "                 thigh+gryo+thigh_x_dir,thigh+gryo+thigh_y_dir,thigh+gryo+thigh_z_dir,\n",
        "                 thigh+accel+thigh_x_dir,thigh+accel+thigh_y_dir,thigh+accel+thigh_z_dir,\n",
        "                 shank+accel+shank_x_dir,shank+accel+shank_y_dir,shank+accel+shank_z_dir]\n",
        "\n",
        "plt.figure()\n",
        "plt.imshow(unnorm_weights, cmap='RdBu')\n",
        "\n",
        "xticks = [0,6.5,14,21.5,29]\n",
        "yticks=[0,1,2,3,4,5,6,7,8,9,10,11]\n",
        "gait_labels = ['0','25','50','75','100']\n",
        "plt.yticks(ticks=yticks, labels=full_features[3:])\n",
        "plt.xticks(ticks=xticks, labels=gait_labels)\n",
        "plt.ylabel(\"Input signals\")\n",
        "plt.xlabel(\"Gait cycle (%)\")\n",
        "plt.title(\"Estimation model weights\")\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "C2u1KMgF-GH_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Here you can see how different signals are weighted during the gait cycle. For example, the thigh gyro (z-direction) has a large negative weight (red) during stance (around 25% of the gait cycle), but large positive weight (blue) during swing (around 75% of the gait cycle)."
      ],
      "metadata": {
        "id": "BEWCGaMnNmXu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Exploring Time-Varying Activities\n",
        "\n",
        "This model can be used as a person transitions between different activities. In this example, we'll specifically look at a time-varying condition where participants transitioned between walking and running on a treadmill, following a step function speed profile with a period of 30 seconds.  Run the code block below to plot the continous measurement of energy expenditure.  "
      ],
      "metadata": {
        "id": "BvJ1dkOd-Kjl"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# @title Energy expenditure when transitioning between walking and running\n",
        "walking_data_path = drive_path+'energyExpenditureData/walk_to_run_sample.csv'\n",
        "raw_IMU_data = np.loadtxt(walking_data_path, delimiter=',')\n",
        "\n",
        "b,a = signal.butter(filt_order, wn, fs=sample_freq) # params for low-pass filter\n",
        "filtered_data = signal.filtfilt(b,a,raw_IMU_data[:,:-1], axis=0)\n",
        "shank_gyro_z_index = 2 # column number for z-axis of shank gyroscope readings used to segment gait cycle\n",
        "\n",
        "peak_index_list = signal.find_peaks(filtered_data[:,shank_gyro_z_index], height=peak_height_thresh, distance=peak_min_dist)[0]\n",
        "peak_index_list_plot = signal.find_peaks(filtered_data[plot_start:plot_start+plot_length,shank_gyro_z_index], height=peak_height_thresh, distance=peak_min_dist)[0]\n",
        "\n",
        "gait_data = np.zeros((3+12*num_bins, 1000)) # storage for up to 1000 gait cycles of formatted data \n",
        "gait_cnt = 0\n",
        "for i in range(len(peak_index_list)-1): # looping through each gait cycle of data\n",
        "    gait_start_index = peak_index_list[i] # index at start of gait\n",
        "    gait_stop_index = peak_index_list[i+1] # index at end of gait\n",
        "    if (gait_stop_index - gait_start_index) <= stride_detect_window: # if gait cycle within maximum time allowed\n",
        "        gait_data[:,gait_cnt] = processRawGait(filtered_data, gait_start_index, gait_stop_index, b, a, weight, height, num_bins)\n",
        "        gait_cnt += 1 # increment number of gait cycles of data stored\n",
        "gait_data = gait_data[:,:gait_cnt] # get rid of empty storage\n",
        "\n",
        "estimates = np.zeros(gait_cnt)\n",
        "for i in range(gait_cnt):\n",
        "    estimates[i] = np.dot(gait_data[:,i], model_weights)\n",
        "\n",
        "plt.figure()\n",
        "plt.plot(estimates)\n",
        "plt.plot([0,len(estimates)],[300,300])\n",
        "plt.plot([0,len(estimates)],[750,750])\n",
        "\n",
        "plt.title(\"Energy expenditure per gait cycle\")\n",
        "plt.xlabel(\"Gait cycles\")\n",
        "plt.ylabel(\"Energy expenditure (W)\")\n",
        "plt.legend(['Per step estimates','Steady-state walking energy','Steady-state running energy'])\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "QuoM8cee-P-n"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Notice how the energy expenditure changes smoothly between between walking and running states."
      ],
      "metadata": {
        "id": "DEl1VQNYQOM2"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exploring the Code on Your Own\n",
        "1. What happens if you filter the data at 20 Hz, instead of 6 Hz?\n",
        "2. What is the minimum threshold you can use to separate the gait cycles during walking?\n",
        "3. Now you can test out data from other activities. Go up to the \"Loading and Formatting Time Series IMU Data\" cell and edit the variable named 'dataset_name' to select running, biking, or stair climbing data. You can change the threshold to separate gait cycles and the filter value to explore how those impact the movement data of other activities."
      ],
      "metadata": {
        "id": "M7ZPIxcY3Rq6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Open-Ended Questions\n",
        "\n",
        "1.   How does the model depend on subject specific parameters like height and weight? \n",
        "2.   How could you use the model to determine the importance of a single signal?\n",
        "3.   What would happen if two activities had the exact same motion (e.g., biking with a different resistance level but at the same pedaling frequency)?\n",
        "4.   What additional information could you incorporate to differentiate between activities with similar motions? \n",
        "5.   How might you include that additional information in the model?\n",
        "\n"
      ],
      "metadata": {
        "id": "egrL0leI-Qie"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JE3IfCEDo8qo"
      },
      "source": [
        "#Feedback\n",
        "\n",
        "This notebook is a work-in-progress and we welcome your feedback on how to increase its usefulness. Email comments to us at [mobilize-center@stanford.edu](mailto:mobilize-center@stanford.edu) or submit an issue on GitHub."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ifTa38Spblu9"
      },
      "source": [
        "\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j9gidtZddPYg"
      },
      "source": [
        "*Version* 1.05\n",
        "\n",
        "Creator: Patrick Slade | Contributors: Patrick Slade, Joy Ku, Matt Petrucci  \n",
        "Last Updated on May 4, 2022\n",
        "\n",
        "This notebook is made available under the [Apache 2.0 license](https://www.apache.org/licenses/LICENSE-2.0)."
      ]
    }
  ]
}