notebooks/offset_computation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5184fc16",
   "metadata": {},
   "source": [
    "# Offset computation\n",
    "\n",
    "We present here some code examples for computing accelerometer and gyroscope offsets using a data from a simple multi-point tumble test as presented in Fayat et al., 2021.\n",
    "\n",
    "This notebook must be seen as an illustration of how such analysis can be performed to allow others to use the pipelines rather than code for replicating the exact panels of the paper."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6d5ea90",
   "metadata": {},
   "source": [
    "## Imports and data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba7fa5c9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ax</th>\n",
       "      <th>ay</th>\n",
       "      <th>az</th>\n",
       "      <th>gx</th>\n",
       "      <th>gy</th>\n",
       "      <th>gz</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.017334</td>\n",
       "      <td>-0.082031</td>\n",
       "      <td>0.989258</td>\n",
       "      <td>2.197266</td>\n",
       "      <td>0.946045</td>\n",
       "      <td>-0.091553</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.020264</td>\n",
       "      <td>-0.090332</td>\n",
       "      <td>0.989990</td>\n",
       "      <td>1.983643</td>\n",
       "      <td>0.915527</td>\n",
       "      <td>0.122070</td>\n",
       "      <td>0.003333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.021729</td>\n",
       "      <td>-0.076172</td>\n",
       "      <td>0.984131</td>\n",
       "      <td>2.166748</td>\n",
       "      <td>1.037598</td>\n",
       "      <td>-0.213623</td>\n",
       "      <td>0.006667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.016113</td>\n",
       "      <td>-0.090820</td>\n",
       "      <td>0.980225</td>\n",
       "      <td>2.014160</td>\n",
       "      <td>1.007080</td>\n",
       "      <td>-0.030518</td>\n",
       "      <td>0.010000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.021240</td>\n",
       "      <td>-0.086670</td>\n",
       "      <td>0.980225</td>\n",
       "      <td>2.105713</td>\n",
       "      <td>0.976562</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.013333</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         ax        ay        az        gx        gy        gz      time\n",
       "0  0.017334 -0.082031  0.989258  2.197266  0.946045 -0.091553  0.000000\n",
       "1  0.020264 -0.090332  0.989990  1.983643  0.915527  0.122070  0.003333\n",
       "2  0.021729 -0.076172  0.984131  2.166748  1.037598 -0.213623  0.006667\n",
       "3  0.016113 -0.090820  0.980225  2.014160  1.007080 -0.030518  0.010000\n",
       "4  0.021240 -0.086670  0.980225  2.105713  0.976562  0.000000  0.013333"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from cycler import cycler\n",
    "import numpy as np\n",
    "import sys\n",
    "sys.path.append(\"..\")\n",
    "from imu_helpers import IMU_Calibrator, get_immobility_smooth\n",
    "\n",
    "# Color cycler for the figures with x, y, z curves\n",
    "rgb = [\"#D55E00\", \"#009E73\", \"#0072B2\"]\n",
    "rgb_cycler = cycler(color=rgb)\n",
    "\n",
    "# Column names\n",
    "col_acc = [\"ax\", \"ay\", \"az\"]\n",
    "col_gyr = [\"gx\", \"gy\", \"gz\"]\n",
    "\n",
    "# Load the data as a pandas dataframe\n",
    "df = pd.read_csv(\"data/example_sensor_calibration.csv\", index_col=0)\n",
    "\n",
    "# Default parameter of IMU_Calibrator for selecting the time intervals\n",
    "# used for the offset estimate\n",
    "threshold = 5  # deg/s\n",
    "sr = 300.  # herz\n",
    "\n",
    "# Grab the accelerometer and gyroscope values\n",
    "acc = df[col_acc].values # Accelerometer values, in G\n",
    "gyr = df[col_gyr].values  # Gyroscope values, in deg/s\n",
    "\n",
    "df.head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5452d317",
   "metadata": {},
   "source": [
    "## Example data visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c38d1db0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, '')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgUAAAFBCAYAAADnvdzqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA450lEQVR4nO3dd3gc5bn38e+9q25bsmXZuGJTbIwxIGxsHOAAoSchtECOgQMEOHEIJIeThHAg5A2k8IaQQM5JfWN6EkoIJTiUHFoIJTQ3QrHBDsUVbFmWq/re7x87kteqK2lXs7v6fa5Ll3Zmnpm5NR7v3vvMU8zdEREREYmEHYCIiIhkBiUFIiIiAigpEBERkYCSAhEREQGUFIiIiEhASYGIiIgAOZYUmFnUzBab2SNhxyIiIpJtciopAC4DloYdhIiISDbKmaTAzMYBnwFuCTsWERGRbJQzSQHw38AVQCzkOERERLJSXtgBpIKZnQSsd/eFZnZUF+XmAnMBSkpKZkyaNKl/ApSsYmZEo9Fe79/c3IyGD5dU6eh+XFi1GoAZFePald+wvYGVm2qpGFTAhGHFbNm2neU1TZRQz+TRw3VvDmCvv/56lbuP6KqM5cINYmY/BM4FmoAioBR40N3/rbN9Kisr/fnnn++nCCVbxGLxiqaysrJeH6OmpgYzIxLJpYo4CUNn96PdfjkAfsFP2u3zm5c+4OL732Du7N35zZkH8sSLL3HCg1XMinzAE98+D0D35gBVWlq60N0P7qpMTtwZ7n6Vu49z94nAHOCZrhICEZGsFjOIdf32nf1f9yQMOfH4QERkQFkxExqKO9xkWD8HI7kk55ICd38WeDbkMERE0qeThCBRDjwZlhDkXFIgIjKQmSoKUqKpqYmqqioaGhrCDqXHCgoKqKioIC+v5x/xSgpERHKZqgx6paqqirKyMsrLy7EsyrTcnerqaqqqqhg1alSP98+JhoYiIrIr5QJ909DQkHUJAcS7sJaXl/e6hkNJgYhIDsmuj7DMlm0JQYu+xK2kQEQkB7k6JUovKCkQEckh2frtVjKDkgIREZEMdfrppzNr1iwOOOAAbr75Zj788EOmTJlCVVUVsViMI488kieeeCJl51PvAxGRHKSGhrnhlltuoby8nNraWmbPns3pp5/ON7/5TS655BJmzZrF1KlTOf7441N2PiUFIiI5TdlBX62/++vUrVqS0mMWja9k5Nk3dVvu5z//OQ8//DAAq1atYvny5Vx00UXcf//9zJs3j4ULF6Y0Lj0+EEnCym2bqGtqDDsMERlAnn32WZ5++mleeOEFFi1aRGVlJXV1dezYsYM1a9YAsG3btpSeUzUFIkmY8Mfr+My4fXnkuIvCDkWkS2pmmHrJfKNPhy1btjBs2DBKSkpYtmwZr7zyCgBXXXUVZ511FhMmTOBLX/oS8+fPT9k5VVMgkqRHVy8NOwSRHlPbgux1wgkn0NTUxEEHHcQ111zDIYccwnPPPcdrr73GFVdcwdlnn01BQQF33HFHys6pmgIRkRzUNhdQDUL2KSws5NFHH223/jvf+U7r6/vvvz+l51RNgUg3YvU7wg5BJGkapkD6QkmBSDeat20MOwSRHnM9N5BeUFIgIpJDrN2DAiUHkjwlBSI9YLdfzisbPgw7DJFutUsF9FhBkqCkQKQbbd9cH1mlXgiSudSmQPpCSYGIiIgASgpEkrDrVy814JJsoNtUekNJgYhIDtHjA+kLDV4kkoz6Isivh4irLbdkhdY7VTdsVjv99NNZvXo1dXV1fPWrXyUajfLmm29y003xoZdvueUWli5dyo033piS8+VETYGZjTezv5rZUjN7y8wuCzsmyR0btjfC8kNg3d5hhyIiA8wtt9zCq6++yiuvvMIvfvELTjjhBB555BEaG+MTtN1xxx184QtfSNn5cqWmoAn4hrsvMrMhwEIze9Ld3w47MMl+NfXN8Rfbh4Yah0hPtB+vQHrr66/NZ0n1upQes7J8NDfNPLnbcm2nTl61ahWf/OQnefTRR5kyZQpNTU3sv//+KYsrJ5ICd18HrAtebzWzpcBYQEmB9J0e0koW0oOu7Jc4dXJJSQlHH300dXV1XHjhhVx//fXss88+nH/++Sk9Z04kBYnMbCJwEPBKyKFIjlArbskm7ccz1A3cV8l8o0+HzqZOPuSQQ1i9ejWLFy9m8eLFKT1nTrQpaGFmg4EHgP909y0dbJ9rZgvMbMHGjRrPXpKkmgLJcO9sXs+DH7wBQHUwgdeGuu1tSuk+zjYdTZ3c4owzzuDQQw9l2LBhKT1nztQUmFk+8YTgLnd/sKMy7j4PmAdQWVmp9FmS0vZGcZzlmzcw+cEf8Y9TvsH+5aNDiUukxZQHbwDAL/gJb9V8BMDyzVVhhiQp0NnUyQAvvvgil12W+jb1OVFTYGYG3Aosdfebwo5Hct+DH8a/lf3+nwtDjkSkaxpsK7fU1NSw7777UlxczDHHHJPy4+dKTcFhwLnAG2a2JFj3LXd/LLyQJFe0e0br6vot2cd01+aEoUOHsnRp+uZfyYmkwN1fQA/MJE30VirZSPet9EZOPD4QSa9d8021O5RMpvEJpC+UFIh0o/WZbENJsLxzmylDkEzT5pasqo2FE4dkJSUFIt3p4HP/o9qtALy84cN+DkZkp6ZYc6fbWpLZ857c2l/hSA5QUiDSQ69VreKZdSsA+NtH74UcjWSy6oZaXtu4Om3Hf3H9B2k7tgxMSgpEeqimoZYNddvCDkOywKnP/Y5j/npr2o7fYXdD1+MC6b2c6H0g0p8WV68hpr7fkoQlNfFJdNw9Le1P8iLR+IvGAmgsAqB+1evAnsTqdx3RsMH0HTAbtZ06edSoUVx77bUA1NbW0tDQwIoVK1J2PiUF0qGPttaTFzEqBhWEHUro2n7+KyGQnnLS02c60nLUFQdDcz4AsZbhjWO71hisyCtJQwSSbrfccgvl5eXU1tYye/ZsnnnmGRYujA+aNmfOHI444oiUnk9JgXRo8o0vAbDl2qPCDSQTdPMNb0dTAyV5Sp6kc56mtOCm+W9x6BcntiYEu55zV83qKdNrX5//NkvWtptOp08qx5Ry08lTuy3Xdurk5cuXM3z4cH784x9TXFzMJZdcktK4VJ8k7dy1OLXzhue6+95/PewQpAvuziNrltEc4rP2dNUuPbCs/d+kkQtzR+LUyYsWLaKyspK6ujqefvppHnjgAX71q1+l/Jz9VlNgZuVJFIu5e026YwF4/aPtlF77LADDS/LZuKMRgDFDCli7taG1XNtlgKhBc8L/u+ljhrCjsZllG3bsUm5oUR41dU0ATK4oYe2WerY17NqFaHxZIas21wNQEDVKC/OoCmIBGDkon/XbG+nKcXuXc9ykcu77x8fsaIzx9vqdzxILo/FvB/XNzv67DcLMWLZhO/uNHMSRew7jriUfsaGL45de++wuf2/EYOKwYt6rriViEHM4aMwQFq/d2nq+wrwI+4wYxGurd2bW5cV5NDR7u78foDgvQm1TrPXfYcqIErbVN7N6S/0u5fYdUcLS4BofP6mcjTsaeWfDDrY1NHPI+FJeWbXzfMau35T2GFbE+5vqWv9NygqjbGtobv27KkryGVQQ5cOaul3OuXfFID6oTvh3jRlEdn3TveCFP/CFSTM7vYYSrofXLOW8l//IZZMP5fsHHJeSYy7fupEfL32OebNOS6p8Oj+mNbdB+iXzjT4dOpo6eeXKlVx//fU8+uijFBcXp/yc1l83lJnVAWvpug4t6u6790s8oyY552juJOmFia/D4Jqwo5CB7M0jW1/efc50zr5rEQBb/++nuPTXP+K3qw5gfHk1K68+H/vGnwEYVLyRNd+MJzGRiCqJu7Ny5UqmTJkSagz19fWcfvrprF27lsmTJ7NhwwaOOuoofvnLXzJu3DgARo8ezSOPPNJu32XLlrH77rt+nJaWli5094O7Omd/tilY6u4HdVXAzBb3VzAivbZ9qJICyRjn3rPzbfP3C1dT26ymYrmis6mTv/Od76TtnP1593wiRWVS4sBRg3g+oRFd9Y5Gtjc0U1aUx9P/rOa0/UZ2uf+OhmYcGFQQTW+gIWh5rAIDr6FhLGixXVZW1rru44/WMOrHi8IKSaRLzbGdtb2L1mzmj2vjVd21Dflc8sA/WrfpIYMko9+SAnff5WGtmY0FWj5R17p7U9sy/am8JJ/ykngL3u4SAoCSHEwGEu1ZXszvPr9f2GFkhePHTOZ/T5jbulxTU4OZqYo2RAuq11Db1Mi/jJzI5sY6vrJgPv89/SSGF6amW15dcxOLqtdw6IgJXZb7yoL5/PaDxSz7zNcYU1za4/N0lKQCrY8E2rr55ZWtr6u2DeHXf9cw3NIz/dnQ8Cog392/F6x6CagBCoA7gR/2VyzStYFWO9BXjx53UdghSBsHl49tfV2WX8TvPvH5lB6/KJrXbUIAcM20Yzhqtz17lRCkmjokSjL686vMmcCNCcsb3f0AYD/gM/0Yh0jPtOnfveeQ4bsst44qJ9LGiKJBnDF+WthhSC9la8+OvsTdr/Wb7p447ub/BOuagdT3qxBJEW/zNPZf96jkM+P2DSkakd7ZXjucA3/2aqfb3Z312xo63T7QFBQUUF1dnXWJgbtTXV1NQUHvBlTrz4aGg80s390bAdz9DgAzKwTCr1sT6UTbMes/2rCBDUVqLyDZ58OaOoZ+77mww8gKwwojXHvocPYetopIFj17iTms2NTItX/fyKb6ZT3evz+TgvuB35jZV9x9B4CZDQJ+EWwTyQq3/2MzUAHT4JOj9go7HBFJg031MS7764aww+h3/ZkU/B/gOmClmbU0id0duDXYJpJVvjLqc1wde4G3zzem3LyDSIGegklm++2ZUzl534pOe8bUN8UozFMtWDLSNfOle3zk1yGFqf94Lk1ivL7+7JLYDFxpZt8F9g5Wr3D32v6KQSSVfvFUFV9q+hkRoHlbNZGEFu8imejkfSu63K6EIHnpSAhajpuOhCBZ/XYHmNl0AHevdfc3gp/ajsqIZIubGo6Nvwhxsh0RkVTpz7TwdjMbZmblnf0Qf5TQK2Z2opm9Y2YrzOzKFMYt0qlV0dEAuJICyTBTyjQKp/Rcf9ZRlAEL6XoMjV616jCzKPBL4DhgNfCamc1397d7czyRRF3dsE8VfoJ7Gj/F/1FSIBlmWdFQ2Bx2FJJt+rNNwcQ0Hn4W8fYJ7wGY2b3AKUCnSUHT6tdZ9ZX2PSGLDziJ2n+0n3EqHYqmHkfd20/2y7lyjeUX4Y3dj4rdm2tsRaWsqds5DfMmK4Hh93Za/geDv8ycy/cMZcS44umnY/nF7HjlrtZ10WHjad60KoRoOhYdOobmmrU92idvt0k0fbw8TRG1N/yC29l4+wX9dr6eWtN2RcX8MMKQASBXptMaCyS+C64GDunNgforIQAyOiE49LBL2GNHNXct7vzDMEzJJATQ/TWOAS+WT+Tw6g9aP9Q9ISFI1rSK+QyNbWGIb2dc88fURIYwu+F13sifTC2FHND0LlGaqbVCni74BKW+jX9pWMiSvCnk08S+Te+xITIMxyj0BrZFSrhox4PcWXwKq6KjeD86jv2b3mWv5lU0kE+pb6PY61m9fDdWR3YjVjadUbEqluXtycGNb7K2bCQTmtey0YayPlLOyFg1RdQzqrmK7VbMSwWVOEaJ11FjQzig6V02RUqptUKmNv2TrTaIEbFqfl98MifXPUN1pIwIMXZvXkcD+RTQSKlv57n8GUxoXsuGSDllvpVG8qjwmtbr8mr+/my2wXxm0HO8lF/Jvk3/5PmCGQzyWg5rXMS70YnEiPBW/iQOanybD6Nj2LfpPT636QneKDmM4V7DusgINkTKyfdGRsaqWR3djQ+jY/goMoLxzevYs3k1Q30Lr+bvz3vR8TRafA6Tcc0fMb75IzZGhvJu3kT2a1xOHs1MaXqP6kgZDZbPu9GJzGh8m0f+PJxDSr/PO3l7sFfzKipim/ggOobDGhZT7HW8lzeeIm9gkw3hH/n7UBHbBBj53siE5rX8M288FbEa1kRGcnjjIqY0vcf3B3+ZU+ue4W8FB3No4xJeyj+QMbENrI2M4IiGBWy3Elbk7U4zEbZZCZ+tf5bDGxbxYNGxTGpeyTMFh9BMhCG+nclNH7AmuhurIqMYnnB9u+QGeyyB9ytbV21ramBwXu8GtZGBwbJttKaOmNmZwAnu/u/B8rnALHf/aptyc4H4zDUTdpvBNf/W36FKNmqOwtLDw45CpHtF26BucPz1yPdh5Ep488id26f9LZy4JDNceONCdz+4qyK50v9kNTA+YXkc0K6+0t3nufvB3V0UEZEw7F5W2LsdK1ZCyWbYe2HCyiwahk8yRr8/PjCz04Bn3H1zsDwUOMrd/9SHw74GTDKzPYg/fpsDnN3VDgcOHcXzZ1zTbr27s725kcF5Ba2DU7g7jR6jpqGWQXkFDMorYOX2GmI4EwcNa92vJ/1WOyofcyfmTl4Pptxte5ymWIyoWbtjt5SLuRMJ/qbO4p0w/wa+tNcsrpx6JJFe9MVtqX3q6/XorEzbsr0ZRKSzfWKxGO5OWVlZ6/aXljzLoUu3dn28Gz/b+lpTJ0uLZ9/bxO8Wr+PWz03t1f5JT5086n0Abpp1Ml9/M2hd6O3v71eOvZh9Skfo3hygSi+8sdsyYbQpuMbdH2pZcPcaM7sG+FNvD+juTWb2FeB/gShwm7u/1ZtjmVnrM7eWDwUzo8CijCwa3Fpu90FD2+3X0/O0FTHr8Ydw2+N0llC0lIsk/E2d+fDkK3oUQ3cxpWof6yT2VJ/P2iRV3R3+mq2/BD7bdSEZkI7acxhH7Tms3873tf2O4IroAzQ1F0BF+8am+5SO6LdYJDuFkS52dM4+Jyfu/pi7T3b3vdz9ur4eT6SFdVMN+/n6/+2nSES69/f9FvKNhl+zevRwms6/IexwJMuEkRQsMLObzGwvM9vTzH5KfPwCkYzUWVvcQq/nyIZXGTT1mP4NSCTR3gt2WRxf0MCFWx6n0IxoQs3h3pGq/o5MslAYjw++SnwCpD8Ey08A3w4hDpGkdNY/Z9HGM5l6Z/b33pEsVrwFirbz+fwo9zU2t9nY9jFb/4Ul2avfkwJ33w5oGGLJGtZJVUHp7LP6ORKRNsrXAXDNkzdwDbDjiDPZ8up94cYkWa3fHx+Y2ZNBj4OW5WFmpoeykjX+79afclLRcsZ9+e6wQxHZxY53nqOpemcDw6ZtG0OMRrJRGG0KKtx3Dsnl7puAkSHEIdIrp9T/lV8My9zRKGUAaVOJtf4PO3sONW1azbuXdj1VskhbYSQFMTPbvWXBzCbQ+WNbkdB11H0xXXOpi/RM5/dh9RP/s8uy12+neXt1p+U91kzjhn+mLDLJTmE0NLwaeMHMWsbbPIKWoYdFRCRpY2q3MGJz8pNNfXTVngBEh0+keeMHAFhxGV7b8XSKkSEjiG3dAHkF0NRAdNg4mjetBotQuM+RREt3Y8eru86PYvnFeGMt0bLRNG+Ot3korjyF2iUPtzl4lEjJMGLbqogMGUn+qH2oX/48RVOPpWHlEmLbquITflmEHQvvB6Bw8pE0ffwusR2b8MY68sdOo3HNm8F5d50kreyU77HliZvw2hoACvaYRcP7r7b7G62gBG/YEQ+pZBixHZvIG70vse2biG35KKnr2tEEXkVTjqZu2TMADDn2MrY+9T8QySNvxJ40ffzuzstQMhR37/jfwAwrKsVrNzP0zJ9Q88fLyR+7P41r3ojvO3g4sU4eEeWPP4jGVYuTin+XU4Yx94GZVQCziae5L7l7v/eVqays9Oeff76/TysZrqMR5F5e8jc+8budkyS9VXUyg/Y/kQmXP97hMTSioaRKdyMaXr3t/3F23WNdHmO/YEbFvZs+5OGar3ZZVnLbfr8lY+c+aAbWE5/te6qZHRFSHCLd0oMCyVTeg7vT9JQ2q1lxWfeFUiCMuQ/+HbiM+KRFS4jXGLwEHN3fsYgkw2Pt30zVpkAyQTPRHpUf+7MagA5rsby5kcZ1yygYt38qQpNM9NvSbouEUVNwGTAT+NDdPwkcBGwIIQ6RpBTmKQGQzBTrSU1BXtczMFo0XwmBhJIU1Ll7HYCZFbr7MmCfEOIQScrgfINJrzJi4hMsO0PtBCRznFn/RNghSI4J4x1udTB40Z+AJ83sYSD55rMi/cwBCmuJ5jUwpjjsaER2GuS1nW7Lr5iwy7LaFEgy+j0pcPfT3L3G3a8lPgfCrcCp/R2HSI9Fw+jBK9K1yT/b2W1uz+/t7IKWVzaayT9fH0ZIksVCfZdz9791X0okXK2NCs00q4xknLyy3VpfF02oxPIK8KYGrLCEvNIRCSV170r39IBUpBuqdJVMtdeP3m23bvQFNwOQN1hDHEvPKSkQ6U7iAF8hDPYl0pnCUZParbO8ghAikVyhpECkGx1WuuoxgmQ63aPSC0oKRJKk/yySFVSbJX2g5tQi3RhbWMIFK1/jrEgM+FR8pd54JeNZmyXds9I9ffkR6YaZcfl7z7FX4w4GTzueoUd+kdEXzAs7LJFOdPLhr6cJkgTVFIj0gOXlM+ZCJQSSBdSmQHoh62sKzOzHZrbMzP5hZg8FoyWKpIGqXyULdPJoSymCJCPrkwLgSWCaux8AvAtcFXI8IiIiWSnrkwJ3f8Ldm4LFl4lPySySBvquJVlAjw2kD7I+KWjjQuDxsIMQEQmNesZIH2RFQ0MzewoY1cGmq9394aDM1UATcFcXx5kLzAUYN04VCiKSw1RjIL2QFUmBux/b1XYzOx84CTjGvfM02d3nAfMAKisrlU6LSM7p4i1QpFtZkRR0xcxOBP4LONLdd4Qdj4hIJjC1gZFeyIU2Bb8AhgBPmtkSM/t/YQckIhIedUmU3sv6mgJ33zvsGERE+tMLG88JPuS3dF5IbQqkF7I+KRBJt8Jx0yiefDi7zflJ2KGIADBlzrVEioZ0vFFtCqQPlBSIdCOSX8geVz8fdhgirYaf+PUkSqmmQHouF9oUiIhINzRLoiRDSYGISE7Rh7/0npICEZFcpIaG0gtKCkREcklnsyQqR5AkKCkQEclyQ4/8YvuVygKkF5QUiIhkuTEXzmPqnUENgbokSh8oKRARyUmqKZCeU1IgIjIAKEWQZCgpEBEREUBJgYiIiASUFIiI5BDXLInSB0oKRERykKlLovSCkgIRkVzSpkviyXXPhBSIZCMlBSIiuSioKTir7rGQA5FsoqRARCSXdDbMcT+HIdlJSYGISE5SGiA9p6RARCSneJslJQeSPCUFIiK5qE3vA+ukq6JIIiUFIiI5bNB+x4YdgmQRJQUiIrmkTUPD0kPmxF/oKYIkIWeSAjO73MzczCrCjkVEJHzxLMA1lbL0QE4kBWY2HjgOWBl2LCIi4VKXROm9nEgKgJ8CV9DZ/wYRkYFGwxxLL+SFHUBfmdnJwBp3f72nY303NDSkJyjJWu5OYWFhn44RiURoaGjQ2PPSZ725HwtGTQageOJ0AEaPHAHUcEiFE4lEqK+v1705ACX7GMmy4XmTmT0FjOpg09XAt4Dj3X2zmX0AHOzuVZ0cZy4wN1icBryZhnAlrgLo8N9B+kzXNr10fdNL1ze9urq+E9x9RFc7Z0VS0Bkz2x94GtgRrBoHrAVmuftH3ey7wN0PTnOIA5aub/ro2qaXrm966fqmV1+vb1Y/PnD3N4CRLcvd1RSIiIhI53KloaGIiIj0UVbXFLTl7hN7UHxeuuIQQNc3nXRt00vXN710fdOrT9c3q9sUiIiISOro8YGIiIgASgpEREQkoKRAREREACUFIiIiElBSICIiIoCSAhEREQkoKRARERFASYGIiIgElBSIiIgIoKRAREREAkoKREREBFBSICIiIgElBSIiIgIoKRAREZFAqEmBmd1mZuvN7M2EdeVm9qSZLQ9+D0vYdpWZrTCzd8zshIT1M8zsjWDbz8zM+vtvERERyXZh1xTcAZzYZt2VwNPuPgl4OljGzKYCc4D9gn1+ZWbRYJ9fA3OBScFP22OKiIhIN0JNCtz9OaC6zepTgDuD13cCpyasv9fd6939fWAFMMvMRgOl7v6Suzvw24R9REREJElh1xR0ZDd3XwcQ/B4ZrB8LrEootzpYNzZ43Xa9iIiI9EBe2AH0QEftBLyL9e0PYDaX+GMGSkpKZkyaNCl10UnOMDOi0Wj3BTvR3NxMvNJKpO/6ej8m0r05sL3++utV7j6iqzKZmBR8bGaj3X1d8GhgfbB+NTA+odw4YG2wflwH69tx93nAPIDKykp//vnnUx27ZLlYLAZAWVlZr49RU1ODmRGJZGJFnGSTVNyPiTZv3gyge3OAKi0t/bC7Mpl4Z8wHzg9enw88nLB+jpkVmtkexBsUvho8YthqZrODXgfnJewjIiIiSQq1psDM7gGOAirMbDVwDXA9cJ+ZXQSsBM4EcPe3zOw+4G2gCbjU3ZuDQ32ZeE+GYuDx4EdERER6INSkwN3P6mTTMZ2Uvw64roP1C4BpKQxNRERkwMnENgUiIiKhampqoqqqioaGhrBD6bGCggIqKirIy+v5R7ySAhERkTaqqqooKyujvLycbBok192prq6mqqqKUaNG9Xj/TGxoKCIiEqqGhoasSwgg3oW1vLy81zUcSgpEREQ6kG0JQYu+xK2kQERERAAlBSIiIhJQUiAiIpKhfvCDH7DffvtxwgkncM4553DjjTcye/Zsnn32WQC+9a1v8e1vfztl51PvAxERkQy0YMECHnroIRYsWEBTUxMzZ85k+vTp3Hrrrfzrv/4rP/3pT3niiSf4+9//nrJzKikQERHpwvq7v07dqiUpPWbR+EpGnn1Tl2VefPFFPvvZz1JcXAzAZz7zGQD2228/zjnnHE499VReeOEFCgoKUhaXHh+IiIhkoK5mtHzzzTcZOnQoH3/8cUrPqZoCERGRLnT3jT5dDjvsMC655BKuvPJKmpqaePzxx7nooot46KGH2LhxI3/961855ZRTmD17NkOHDk3JOVVTICIikoFmzpzJSSedxPTp0znjjDOYMWMG7s63vvUtbr75ZiZPnswll1zC1772tZSdU0mBiIhIhvrGN77B22+/zYMPPsi7777L0UcfzdKlSxk/fjwAX/3qV7n99ttTdj49PhAREclQF198MUuXLqWuro5zzz2X6dOnp/V8SgpEREQy1O9///t+PZ8eH4iIiAigpEBEREQCSgpEREQEUFIgIiIiASUFIiIiAigpEBERkYCSApFuxGLOwTc9yz0P3h92KCIywLSdOvmGG25g5syZrduXL1/OrFmzUna+jE0KzOwDM3vDzJaY2YJgXbmZPWlmy4PfwxLKX2VmK8zsHTM7IbzIJddsrW9i4ZqtzH2+iR0rXg47HBEZIBKnTr7//vtZuHAh0WiU0tJSlixZAsAdd9zBeeedl7JzZvrgRZ9096qE5SuBp939ejO7Mlj+LzObCswB9gPGAE+Z2WR3b+7/kCWXNW+vDjsEEelnX39tPkuq16X0mJXlo7lp5sldluls6uSLLrqIO++8k/33358//vGPvPTSSymLK2NrCjpxCnBn8PpO4NSE9fe6e727vw+sAFJXnyIS8KaGsEMQkQGis6mTTz/9dP7yl7/w6KOPMn36dIYPH56yc2ZyTYEDT5iZA79x93nAbu6+DsDd15nZyKDsWCCxXnd1sE4kpbYufIjSGaeGHYaI9KPuvtGnS2dTJxcVFXH88cdz6aWXMm/evJSeM5NrCg5z9+nAp4BLzeyILspaB+vapVhmNtfMFpjZgo0bN6YqThlAiiamdzISkWRsbazjh/94muZYLOxQJI06mjq5rKwMgLPOOgsz4/jjj0/pOTM2KXD3tcHv9cBDxB8HfGxmowGC3+uD4quB8Qm7jwPWdnDMee5+sLsfnMrqFhk4LJLJlWsyUFzx2qN8a+HjPLTyjbBDkTRrO3VyyyyJL774Il/4wheIRqMpPV9GvsOZ2SAg4u5bg9fHA98D5gPnA9cHvx8OdpkP3G1mNxFvaDgJeLXfAxcR6QdbG+sBqGtuCjkSSbeOpk7+3Oc+x3vvvceTTz6Z8vNlZFIA7AY8ZGYQj/Fud/+Lmb0G3GdmFwErgTMB3P0tM7sPeBtoAi5VzwNJtYZIlEs2buDesAMRkQGjo6mTH3jggbSdLyOTAnd/Dziwg/UbgWM62ec64Lo0hyYDWEMkjz/s2KakQDJGJ43TRXqt06TAzMqT2D/m7jWpC0dERLpjHTWtFkmBrmoK1gY/Xd1+UWD3lEYkIiJJ8fadrET6pKukYKm7H9TVzma2OMXxiGScxpiap0hmseC7WmeD24j0VlddEj+RxP7JlBHJah9u2xR2CCK70NMDSZdOkwJ3rwMws73MrDB4fZSZ/YeZDU0sI5LL9F1MRAaKZHofPAAcbGZ7A7cSjAkAfDqdgYlkCn0rk0wT21ET/62UNef94Ac/4J577mHcuHFUVFQwZcoU/vSnP7Vuf+ONN1i+fDkTJkxIyfmSSQpi7t5kZqcB/+3uP1dbAhGR8Gxb9DCMnhZ2GJJmiVMnNzU1MXPmTKZPn87ChQsB+NWvfsVzzz2XsoQAkksKGs3sLOIjCH42WJefsghERKRX5r54P2dOPJDB+YVhh5LTvj7/bZas3ZLSY1aOKeWmk6d2WaazqZNbtt122208++yzKY0rmbkPLiDeoPA6d3/fzPYA2g+xJJKjOur2defy1/ju4idCiEZkp6bllRz666fDDkPSpLPeJevWrWPu3LncfffdDB48OKXn7GrwonnA48BT7v4fCUG+T3zuAZEB6wsv/AGAaw5K7QxlIslobedSN4Q3PlSX2XTr7ht9unQ0dfKFF17InDlz+OEPf8jkyZNTfs6uagpuIz7U8GNm9rSZ/ZeZtRt6WCTXmZoaSoYxNTAcEDqaOnnHjh0sWLCA7373u8yYMYMZM2awdm27SYF7rasuiS+7+7Xu/i/A54lPQPQNM1tiZreZ2edTFoVIBtOocZJplKYOHG2nTj7ppJPYvn07CxcubP0ZM2ZMys6XTJsC3H2ju9/j7ue5eyXwS+LTE4uISAZ5o3oddvvlPLZqaeu6dzdv0MicWeriiy9mxowZzJw5k9NOO43p06en9Xzd9j4ws693sHoz8GjqwxERkW510ADt0VVvM7l0BC9t+BCAP618k0+P35dV22rY58EfcdnUf+G7U47q50ClrzqaOjmdkqkpOBi4GBgb/MwFjgJuNrMr0heaSAZyuPqxpVBXAnWD2F7fFHZEMgDFPAKxXd++T3rqNiY/+KPWRwstacP6um0APPfRe/0XoGStZMYpGA5Md/dtAGZ2DXA/cASwELghfeGJZJitw/m/b60AZgKw+w+eYuP3Tww3JhlwHq85Hzbt7Iq27OOtra/Ndp0sqaVNjKZb7jl3b72e2aQvE2UlU1OwO9CQsNwITHD3WqC+12cWyRKt/71iebB6yi7bqnc09ns8InW+a9/0xWu2QGMBsLMRYssQyC0fEOpF0zMFBQVUV1dn3UyU7k51dTUFBQW92j+ZmoK7gZfN7OFg+bPAPWY2CHi7V2cVySK7vCfEkvkvI9K/zr5rEfAJmPa3dtMqO8DWcmJDlRT0REVFBVVVVWzYsCHsUHqsoKCAioqKXu3b7Tucu3/fzB4DDieehF7s7guCzef06qwiIpIWLbXdLbnsyk118OH+fLBjGxwdWlhZJy8vj1GjRoUdRr9LqksiUAxscff/Bj4MhjoWGRBiWVZ9KANbpE1NQW1DDIC6Hb2rTpaBpdukIGhY+F/AVcGqfDT3gYhIRnqtahUAb9V8DCQMdKTcVpKQTE3BacDJwHYAd18LDElnUCLZZFHV6rBDEGn1h/dfh+YoCzfG78vm7VUAeLMaxUr3kkkKGjxeD+UAQQPDjGNmJ5rZO2a2wsyuDDseyR1Nmz/ucvuMe2/jr+tW9FM0Il2rrxoJSw+HhiLW3noRsYbtAHgsFnJkkg2SSQruM7PfAEPN7IvAU8DN6Q2rZ8wsSnzo5U8BU4GzzCycaa0k5zTVbu26wPsHccHf7uqfYES64rCtqiz+uqGQNS/+jj+veTfYpt4H0r1keh/8xMyOA7YA+wDfcfcn0x5Zz8wCVrj7ewBmdi9wCuoyKamQxOAlH26r7YdARLrRWIg3R+OvIzFmHf4f8E68cre+KSMreSXDWLYNzNARMzsDONHd/z1YPhc4xN2/0tk+0dETveSib8f3J/5spOV3T7Qr383nR/LnaSnZ1Tk9iQ+sxNJdlW0pFy9TF4m/sZQ2N+xSytrcL22P2HZK1123t4/AWs7tbcoa1FmUrXnxFtOlTQ10JJkpZJP5N+3qMtbFCqmvHdHtMYYUf7zL+Xb+bQ7WUaVc5//GItD2Xtppa+1uCSV6cA8Nrmaw2hYMWNuu++JCdz+4qzKd1hSY2Va6eD9199I+xJZqHf2vaP95bTaX+NwNREZNxLywtWDLAbr675XUf71uP4F6/yGwy54py+U6jqcu+JRspqClNUm3p+7p+s63GTjE4gEEcRS12S9FH6Zt33Vbl3cev6kxya5cwf2UeCgj3sDLonkoAcgQlj1fhHa+L3V27/QwKWjOw9za7NNZ6pF7LIv+7cPSaVLg7kMAzOx7wEfA74jfNeeQeb0PVgPjE5bHAWvbFnL3ecA8gMrKSn/+ylP7JTjJHrGgMVZZWVnruhdXrOLwXy/pdt8t150LQE1NDWZGJJLsMCAiHevofgSwb/y55RXk1UNTIUx+mRuidVyxfm/YOG5n4fFvwar9KCiqYf03TwbQvTlAlV53frdlkrkzTnD3X7n7Vnff4u6/Bj7X5+hS6zVgkpntYWYFwBxgfsgxSY6wSH63ZZ780ux+iESkLYNh62Da36Cgnr2P+wLkx6ekKRhUzS2fP5BzYh8FJfUtWbqXTFLQbGbnmFnUzCJmdg6tlbqZwd2bgK8A/wssBe5z97fCjUpyhXX3rWrfFzh2cvdtDkTSIqHWP88iTMmPP+6aUbSViw7ZneNmz4kXC9oJiXQlmaTgbODzwMfBz5nBuozi7o+5+2R338vdrws7HskdHVW1ThqfMC6BnlNKqHbef/mRKEXB7ZofrI8OHg4QtGsR6VoyXRI/IN69T2Rg6qD91XN71jJ6dTO4vn1J5siPRNl5w7a9N3O/IaH0Xac1BUFL/S4lU0Yk200dObjduijO0fl/AZzK4aP7PyiRwJi6mtbXeZEI+UPiNQP5I/cGIBbUGCglkGR0VVNwpZlVdbHdgMsIWvOL5KpIxMBi4Ak5tDsXbvszz0wfTEF0fOc7i6TZWeuWcOPuewIQtQhEghEygjYEk3crhpIaxu69JbQYJXt0lRT8DfhsN/tn2siGIilnZrsMqPDKZYfD4tdbW3NHOxyYSKR/RBJuzkjiiAbB6oKowZ6vM3jI2H6PTbJPV+MUXNCfgYhktp2Vr7N2H8aGxU5zkAxYEsMgi6RLxHdOdGQGeXnxbKC4IH5/tj4+0G0qSdBXHJHeSBjuWe+1EqYIu85+uOe4Zhi1gtNmxseYaxnKPqKsQJKgPioiSYm/oR7U+DYtT9Va5oFQh0QJU+KgRIZx0yGfpSTvL/zbpEogsaGhkgLpXpc1BcFgRZ/vr2BEMlHim+mlO+4BIFJStnOOgxyYVEyyVzShpiAvEmFE0WB+c9gZFAbjEgwKJhWbOHhYKPFJdumypsDdY2b2FeC+fopHJKONiFUDUH70JZTX10PVemoaNG2yhCfx8UFFYfvpkacNG819R53LCWP3wWvr+zM0yULJtCl40swuN7PxZlbe8pP2yEQy0O7N6wCwvHxKZ50JQFX99jBDkgHO3Kk/73qWnPJ1Jg7p+K35zD0OpLSgqMNtIomSaVNwYfD70oR1DuyZ+nBEMls+Ta2vi6LxiZLGlQwNKRqR+OODgmgeB5aPCTsUyQHJDHO8R38EIpINdp2FXl29JCS2c4hta9P7QKQvun18YGYLzOwSMxvaD/GIZLQhM05rfd3SwFCtuqXfJcy5EVH/F0mhZNoUzAHGAgvM7F4zO8E0WosMIDGPQcEOAMb/x4Ot61veitX/W8IULS4NOwTJId0mBe6+wt2vBiYDdwO3ASvN7LtqcCgDQQyHvRbBpFd2Wb/v0JEUR/P57kHHhxSZDFhTn299OeKMH4YYiOSapAYvMrMDgAuATwMPAHcBhwPPAJXpCk4kE+RZFKLN8Z8EQ/KL2HGe3pAlBJEYDP0IakYRKWo/i6dIb3WbFJjZQqAGuBW40t1bOrq+YmaHpTE2kYxQnJcfdgginYpp8CxJoWRqCs509/c62uDup6c4HhERSUo8GYgpJ5AUSiYpOLWDdoWbgYXuviTlEYmISPeCt2UNsy2plEzvg4OBi4n3QBgLzAWOAm42syvSF5qIiHRHNQWSSsnUFAwHprv7NgAzuwa4HzgCWAjckL7wRESkYy2PD5QVSOokU1OwO9CQsNwITHD3WkCza4iIhGFrBQDvbtDcG5I6ySQFdwMvm9k1QS3Bi8A9ZjYIeDvVAZnZtWa2xsyWBD+fTth2lZmtMLN3zOyEhPUzzOyNYNvPNLiSiOS85viohh9v1XczSZ1k5j74vpk9RnxcAgMudvcFweZz0hTXT939J4krzGwq8dEV9wPGAE+Z2WR3bwZ+Tbytw8vAY8CJwONpik1EJHzl62DjOA4coxENJXWSGrzI3RcSbz8QplOAe4NxEt43sxXALDP7ACh195cAzOy3wKkoKRCRXBZtBKA4P9pNQZHkJfP4IAxfMbN/mNltZjYsWDcWWJVQZjU7e0Ss7mC9iEjuGroeBm3i9ANGhR2J5JBQkgIze8rM3uzg5xTijwL2Ij588jrgxpbdOjiUd7G+o/PODWZ9XLBx48a+/yEiImEpqIM9/sHYsuKwI5EcktTjg1Rz92OTKWdmNwOPBIurgfEJm8cBa4P14zpY39F55wHzACorK9WPR3okzzK1Yk0GqsF5hWGHIDkmlKSgK2Y22t3XBYunAW8Gr+cDd5vZTcQbGk4CXnX3ZjPbamazgVeA84Cf93fcktuqz/4eUSUFkmEWnfyfYYcgOSbjkgLgBjOrJP4I4APgSwDu/paZ3Ue8G2QTcGnQ8wDgy8AdQDHxBoZqZCgpNaywJOwQRNqZVDYi7BAkx2RcUuDu53ax7Trgug7WLwCmpTMuERGRXKf6UBEREQGUFIiIiEhASYGIiIgASgpEREQkoKRAREREACUFIiIiElBSICIiIoCSAhEREQkoKRARERFASYGIiIgElBSIiIgIoKRAREREAkoKREREBFBSICIiIgElBSIiIgIoKRAREZGAkgIREREBlBSIiIhIQEmBiIiIAEoKREREJKCkQERERAAlBSIiIhIIJSkwszPN7C0zi5nZwW22XWVmK8zsHTM7IWH9DDN7I9j2MzOzYH2hmf0hWP+KmU3s5z9HREQkJ4RVU/AmcDrwXOJKM5sKzAH2A04EfmVm0WDzr4G5wKTg58Rg/UXAJnffG/gp8KO0Ry8iIpKDQkkK3H2pu7/TwaZTgHvdvd7d3wdWALPMbDRQ6u4vubsDvwVOTdjnzuD1/cAxLbUIIiIikry8sANoYyzwcsLy6mBdY/C67fqWfVYBuHuTmW0GhgNVXZ3I3WloaEhR2JIrYrEYBQUFfT5OQ0MDkYia7EjfxGIxCgsLU3Y8d6exsVH3pnQqbUmBmT0FjOpg09Xu/nBnu3WwzrtY39U+HcU0l/gjCID6ioqKNzuJQ/qugm4SM+k1Xdv0yprraxfeGHYIvZE11zdLdXV9J3S3c9qSAnc/the7rQbGJyyPA9YG68d1sD5xn9VmlgeUAdWdxDQPmAdgZgvc/eCOyknf6fqmj65teun6ppeub3r19fpmWh3SfGBO0KNgD+INCl9193XAVjObHbQXOA94OGGf84PXZwDPBO0OREREpAdCaVNgZqcBPwdGAI+a2RJ3P8Hd3zKz+4C3gSbgUndvDnb7MnAHUAw8HvwA3Ar8zsxWEK8hmNN/f4mIiEjuCCUpcPeHgIc62XYdcF0H6xcA0zpYXwec2Ysw5vViH0merm/66Nqml65veun6plefrq+ppl1EREQg89oUiIiISEgGZFJgZicGwyivMLMrw44nm5nZeDP7q5ktDYauvixYX25mT5rZ8uD3sLBjzWZmFjWzxWb2SLCs65siZjbUzO43s2XBffwJXd/UMLOvBe8Lb5rZPWZWpGvbe2Z2m5mtN7M3E9Z1ej07mzagKwMuKQiGTf4l8ClgKnBWMLyy9E4T8A133xeYDVwaXM8rgafdfRLwdLAsvXcZsDRhWdc3df4H+Iu7TwEOJH6ddX37yMzGAv8BHOzu04Ao8Ybgura9dwc7h/hv0eH17GbagE4NuKQAmAWscPf33L0BuJf4UMnSC+6+zt0XBa+3En9DHcuuw0/fyc5hqaWHzGwc8BngloTVur4pYGalwBHEezHh7g3uXoOub6rkAcXBGDIlxMeX0bXtJXd/jvbj8HR2PTucNqC7cwzEpKB1WORA4pDJ0gfBDJUHAa8AuwXjSxD8HhliaNnuv4ErgFjCOl3f1NgT2ADcHjyeucXMBqHr22fuvgb4CbASWAdsdvcn0LVNtc6uZ68+6wZiUpD0sMiSPDMbDDwA/Ke7bwk7nlxhZicB6919Ydix5Kg8YDrwa3c/CNiOqrNTIni2fQqwBzAGGGRm/xZuVANKrz7rBmJS0NlQytJLZpZPPCG4y90fDFZ/HMxuSfB7fVjxZbnDgJPN7APij7qONrPfo+ubKquB1e7+SrB8P/EkQde3744F3nf3De7eCDwIHIqubap1dj179Vk3EJOC14BJZraHmRUQb4gxP+SYslYw7PStwFJ3vylhU+Lw0+ezc1hq6QF3v8rdx7n7ROL36jPu/m/o+qaEu38ErDKzfYJVxxAfUVXXt+9WArPNrCR4nziGeJsjXdvU6ux6djhtQHcHG5CDF5nZp4k/p40CtwWjKEovmNnhwPPAG+x85v0t4u0K7gN2J/7mcKa7dzhRlSTHzI4CLnf3k8xsOLq+KWFmlcQbcRYA7wEXEP/CpOvbR2b2XeBfifdSWgz8OzAYXdteMbN7gKOIz4T4MXAN8Cc6uZ5mdjVwIfHr/5/u/nj7o7Y5x0BMCkRERKS9gfj4QERERDqgpEBEREQAJQUiIiISUFIgIiIigJICERERCSgpEJE+C2YavCR4PcbM7g87JhHpOXVJFJE+C+a9eCSYDU9EslRe2AGISE64HtjLzJYAy4F93X2amX2B+KxtUWAacCPxQYLOBeqBT7t7tZntRXxK8xHADuCL7r6sv/8IkYFOjw9EJBWuBP7p7pXAN9tsmwacTXza1uuAHcHkQy8B5wVl5gFfdfcZwOXAr/ojaBHZlWoKRCTd/uruW4GtZrYZ+HOw/g3ggGCGzUOBP8aHyAegsP/DFBElBSKSbvUJr2MJyzHi70ERoCaoZRCREOnxgYikwlZgSG92dPctwPtmdibEZ940swNTGZyIJEdJgYj0mbtvBF40szeBH/fiEOcAF5nZ68BbwCmpjE9EkqMuiSIiIgKopkBEREQCSgpEREQEUFIgIiIiASUFIiIiAigpEBERkYCSAhEREQGUFIiIiEhASYGIiIgA8P8BNX8VNynsDPgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gyr_norm = np.linalg.norm(gyr, axis=1)\n",
    "is_immobile = get_immobility_smooth(gyr_norm, threshold=threshold)\n",
    "\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(8, 5))\n",
    "ax1.set_prop_cycle(rgb_cycler)\n",
    "ax2.set_prop_cycle(rgb_cycler)\n",
    "\n",
    "# Plot the acc and gyr data\n",
    "df.plot(y=col_acc, x=\"time\", ax=ax1, ylabel=\"acc [G]\")\n",
    "df.plot(y=col_gyr, x=\"time\", ax=ax2, ylabel=\"gyr [deg/s]\")\n",
    "\n",
    "# Add the intervals below the threshold as a shaded area\n",
    "# These time intervals will be used for the offset estimate\n",
    "for ax, lim in zip([ax1, ax2], [4, 1000]):\n",
    "    ax.fill_between(df.time, -lim, lim,\n",
    "                    color=\"black\", alpha=.05,\n",
    "                    where=is_immobile)\n",
    "    ax.set_ylim([-lim, lim])\n",
    "ax1.set_xlim([df.time.iloc[0], df.time.iloc[-1]])\n",
    "fig.suptitle(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "180f55ef",
   "metadata": {},
   "source": [
    "## Offset estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "2be7b847",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "acc offsets [G]:      [ 0.02793173 -0.07261109 -0.01655679]\n",
      "gyr offsets [deg/s]:  [ 2.13623047  0.9765625  -0.09155273]\n"
     ]
    }
   ],
   "source": [
    "c = IMU_Calibrator(sr=sr, n_static_positions=3)\n",
    "c.compute_offsets(acc, gyr)\n",
    "print(\"acc offsets [G]:     \", c.acc_offsets)\n",
    "print(\"gyr offsets [deg/s]: \", c.gyr_offsets)"
   ]
  }
 ],
 "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}