test_GPS_OLS_NUUK.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime as dt\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import glob\n",
    "from tidegauge_functions import read_GPS_SONEL\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datum = 'dNUUK_43013M001_NGL14.neu'\n",
    "dir_in = 'data/GPS'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/GPS\\dNUUK_43013M001_NGL14.neu\n"
     ]
    }
   ],
   "source": [
    "pattern = os.path.join(dir_in, datum)\n",
    "print(pattern)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames = sorted(glob.glob(pattern))  # , key=alphanum_key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "data/GPS\\dNUUK_43013M001_NGL14.neu\n",
      "            North  East  Vertical  NorthSTD  EastSTD  VerticalSTD\n",
      "Year                                                             \n",
      "2008-09-30  -47.0  75.8     -31.1       3.3      2.3          5.4\n",
      "2008-10-02  -47.3  76.9     -28.8       2.1      1.4          3.5\n",
      "2008-10-03  -46.3  76.4     -27.5       2.1      1.4          3.5\n",
      "2009-01-01  -43.1  74.3     -35.0       2.1      1.4          3.4\n",
      "2009-01-02  -42.6  74.4     -36.5       2.1      1.4          3.5\n"
     ]
    }
   ],
   "source": [
    "for f, filepath in enumerate(filenames):\n",
    "    df = read_GPS_SONEL(filepath)\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "def decompose(df, var, site, window):  \n",
    "    \n",
    "    decomposition = seasonal_decompose(df[var], period = window)\n",
    "\n",
    "    trend = decomposition.trend\n",
    "    seasonal = decomposition.seasonal\n",
    "    residual = decomposition.resid\n",
    "    Amplitude = max(seasonal) - min(seasonal)\n",
    "    \n",
    "    decomposition.plot();\n",
    "    df['denoised'] = df[var] - residual\n",
    "    df['trend'] = trend\n",
    "    #print(f'\\n\\n{filepath}')\n",
    "    print(df.head())\n",
    "    #print(f'{site} Seasonal Amplitude =', Amplitude)\n",
    "    #seasonal_amp = pd.DataFrame(list(zip(list({site}), list(Amplitude))),\n",
    "                           # colums = ['Site_Name', 'Seasonal_Amplitude'])\n",
    "    #seasonal_amp.to_csv('GPS_Seasonal_Amplitude_SONEL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Change Point Detection\n",
    "def change_pt_detection(df, var, points):\n",
    "    model = \"l1\"  \n",
    "    algo = rpt.Dynp(model=model, min_size=3, jump=5).fit(points)\n",
    "    my_bkps = algo.predict(n_bkps=10)\n",
    "\n",
    "    \n",
    "def change_pt_plot()\n",
    "    rpt.show.display(points, my_bkps, figsize=(10, 6))\n",
    "    plt.suptitle(f\"{site}Change Point Detection\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_calc_OLS(df, var, start, end):\n",
    "    df = df[start:end]\n",
    "    x, y = np.arange(len(df[var].dropna())), df[var].dropna()\n",
    "    x = sm.add_constant(x)\n",
    "    model = sm.OLS(y, x)\n",
    "    res = model.fit()\n",
    "    return res\n",
    "\n",
    "def test_plot_OLS(df, res, res2, site, var, trend_array, trend_array2, start, end, start2, end2):\n",
    "    df1 = df[start:end]\n",
    "    df2 = df[start2:end2]\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(12,6));\n",
    "    ax.plot(df1[var].dropna().index, df1[var].dropna().values, \n",
    "            label='trend', marker='.', linestyle='', color = 'darkgrey')\n",
    "    ax.plot(df1[var].dropna().index, [res.params.x1*i + res.params.const for i in np.arange(len(df1[var].dropna()))], \n",
    "            marker = '.', linestyle = '', color = 'purple')\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, res.resid.values);\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, np.abs(res.resid.values));\n",
    "    # ax[1].hlines(0, 0, len(res.resid), color='k');\n",
    "    # ax[1].set_title(\"Residuals\");\n",
    "    ax.plot(df2[var].dropna().index, df2[var].dropna().values, \n",
    "            label='trend', marker='.', linestyle='', color = 'darkgrey')\n",
    "    ax.plot(df2[var].dropna().index, [res2.params.x1*i + res2.params.const for i in np.arange(len(df2[var].dropna()))],\n",
    "           marker = '.', linestyle = '', color = 'purple')\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, res.resid.values);\n",
    "    # ax[1].plot(df['Vertical'].dropna().index, np.abs(res.resid.values));\n",
    "    # ax[1].hlines(0, 0, len(res.resid), color='k');\n",
    "    plt.axvline(dt.datetime(2007, 9, 7), color= 'red', linestyle='--')\n",
    "    ax.set_title(f\"Trend {start}:{end} = {trend_array * 1:.2f} mm/yr, Trend {start2}:{end2} = {trend_array2 * 1:.2f} mm/yr\");\n",
    "    # ax[1].set_title(\"Residuals\");\n",
    "    plt.suptitle(f\"{site}\")\n",
    "    plt.savefig(f'figs/test_GPS_OLS_{site}.png')\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "data/GPS\\dNUUK_43013M001_NGL14.neu\n",
      "f: 0\n",
      "            North  East  Vertical  NorthSTD  EastSTD  VerticalSTD  denoised  \\\n",
      "Year                                                                          \n",
      "2008-09-30  -47.0  75.8     -31.1       3.3      2.3          5.4       NaN   \n",
      "2008-10-02  -47.3  76.9     -28.8       2.1      1.4          3.5       NaN   \n",
      "2008-10-03  -46.3  76.4     -27.5       2.1      1.4          3.5       NaN   \n",
      "2009-01-01  -43.1  74.3     -35.0       2.1      1.4          3.4       NaN   \n",
      "2009-01-02  -42.6  74.4     -36.5       2.1      1.4          3.5       NaN   \n",
      "\n",
      "            trend  \n",
      "Year               \n",
      "2008-09-30    NaN  \n",
      "2008-10-02    NaN  \n",
      "2008-10-03    NaN  \n",
      "2009-01-01    NaN  \n",
      "2009-01-02    NaN  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd5hURdaH39M9eYCBgSGHAUGQnJMKKibEFXXNa9ZVP8Oadl3QNe6quK6uYXXNcRXTmlYMiKIkJUrOYchhAIEhTOo+3x/3dk9PT3dPz3RkqPd5+unbVffW/fW93XVuVZ06JaqKwWAwGAzJhiPRAgwGg8FgCIQxUAaDwWBISoyBMhgMBkNSYgyUwWAwGJISY6AMBoPBkJQYA2UwGAyGpMQYKIPhMEVE7haRVyIsI19EVERSoqXLYIgWxkAZDDFARL4RkYcCpI8WkW01NQgicoKIbPJNU9VHVPXaSLUaDMmKMVAGQ2x4A7hMRMQv/TLgHVUtD7cg07oxHKkYA2UwxIZPgVzgeE+CiDQCzgTeEpF0EfmHiGwQke0i8oKIZNr7nSAim0TkzyKyDRgPfAW0FJH99quliDwgIv/xKf84EZkhIntEZKOIXGmnjxKRX0Rkn53+QNyugsEQAcZAGQwxQFUPAR8Al/skXwAsV9UFwGPA0UBvoCPQCrjPZ9/mWAaunV3GSGCLqtazX1t8zycibbGM2LNAnl3ufDv7gF1GQ2AU8H8icnb0vq3BEBuMgTIYYsebwPmelhGWkXjT7vb7PXC7qu5W1SLgEeAin2PdwP2qWmIbu+r4HTBJVcerapmq7lLV+QCq+oOqLlJVt6ouxGqRDY/SdzQYYobp2zYYYoSqThORQmC0iMwCBgDnYrVwsoC5PkNUAjh9Di9U1eIanK4NsCZQhogMAsYB3YE0IB34sAZlGwwJwbSgDIbY8hZWy+kyYKKqbgd2AoeAbqra0H7lqGo9n+P8lxmobtmBjcBRQfLeBT4H2qhqDvAClkE0GJIaY6AMhtjyFnAyVpfemwCq6gZeBv4pIk0BRKSViJwWopztQGMRyQmS/w5wsohcICIpItJYRHrbefWB3apaLCIDgUsi/1oGQ+wxBspgiCGqWgDMALKxWjEe/gysBn4WkX3AJKBziHKWY40drbW99Fr65W8AzgDuBHZjOUj0srNvBB4SkSIsR4wPIv9mBkPsEbNgocFgMBiSEdOCMhgMBkNSYgyUwWAwGJISY6AMBoPBkJQYA2UwGAyGpKROT9Rt0qSJ5ufnJ1qGwWAwGEIwd+7cnaqa559epw1Ufn4+c+bMSbQMgyFsJi/fweyC3dx1epdESzEY4oaIrA+UntAuPhF5TUR2iMhin7RcEflWRFbZ743sdBGRZ0RktYgsFJG+iVNuMMSGq96YzfM/BIxYZDAccSR6DOoN4HS/tDHAd6raCfjO/gxWNOdO9us64N9x0mgwGAyGBJBQA6WqU7BmvfsyGjskjP1+tk/6W2rxM9BQRFrER6nBYDAY4k2iW1CBaKaqWwHs96Z2eiusgJgeNtlplRCR60RkjojMKSwsjLlYg8FgMMSGZDRQwQgUfblKnCZVfUlV+6tq/7y8Kk4hBoPBYDhMSEYDtd3TdWe/77DTN2GteeOhNbAFg8FgMNRJktFAfQ5cYW9fAXzmk3657c03GNjr6Qo0GAwGQ90jofOgRGQ8cALQREQ2Afdjrfz5gYhcA2wAzrd3/xJrOYHVwEHgqrgLNhjihKris9quwXBEklADpaoXB8kaEWBfBW6KrSKDITlwuZUUpzFQhiObZOziMxiOeFxmnTaDwRgogyEZcbsTrcBgSDzGQBkMwJszCli380CiZXgxLSiD4TA0UCJyuoissGPyjan+CIMhNKXlbu7/fAkn/uOHREvxMnPtrkRLCMrUVYUs37avSrqq8tKUNRQWlURUfnGZi/NfmMH6XTV7YNi1v4Si4rKIzm1ILg4rAyUiTuA5rLh8XYGLRaRrYlUZkoHS8tr3iZWUuwJuA7jdyuvT14Ws+Fbv2M/qHUW1Pn8grnlzDsVlFVrcbqWouIzP5m/m8wVbeG7yan5cWREp5cM5G5mwMHqzLlQVl7uiFTfuq+Xkj5nA8z+s5rJXZ3H6U1NRVQp2HqDcZV37aat38siXyxn78SKASvrDpbTczaBHvmN2wa+MfHpq0P0mL9/B5BU7KHe5cbuVFduK6Pe3SfR4YGKNz+mL263sPlDKx/M28fSkVRGVZYgc0cOoK0FEhgAPqOpp9uexAKr6aKD9+/fvr2a5jeTG7Vbembmeto2zaZGTwT+/XUm3lg2YsWYX1w3rQMem9SgsKmF2wW76tWtE37aNqrhfL92yjzOemcorl/fn5K7NqpyjpNzF3oNlPPi/pXy7dDsndsnjN71a8v7sjdx1WhfGfrKQxZsrWgR3nHI0HfKyufndX6qU1blZfb74w3E4RSgqLgeBXg9aleKQDo2589Sj6deuqkaAwqISBjw8CYAXL+vHMc0b0CwnnfQUJy9PWcvDXy6rcswpXZvx6Lk9+L//zGV2wa8Br+EbVw3gytdnA7DsodPJTHNW2cfzPxcRistcdLn360r5LXMy2LK3OGD5NaVL8/q8ckV/jntsMqd1a0Z6ipO/ndOdBhmp3n0Wb97Lqh1FnNSlGc9NXs3OohLK3UpJuYtvlmz37nfZ4HbkZqdxsLScl6euA+CKIe1486eAqzMA8OlNx9K7TUPv9/5p7S6GdGhc5Z6Ulrt58cc1FJWU89KUtUHLy0pzkuIQUp0OUpxCisNBqlNwBkhLcVifU50OUhzizUtxCqlV8vzK8UmrKMfetvMql1P5fE5HRV6g8zkdkrRTF0Rkrqr2r5J+mBmo84DTVfVa+/NlwCBVvdlnn+uwop3Ttm3bfuvXB/8hGxLPOzPXc88ni6vf0Yfc7DTGnduDbq1yyEp18v3yHdz54QJG927J0xf1AaDc5abjPV/FQnLS07R+OqN6tmBHUQm79pfw81r/eMyJ4apj8znlmGZc8srMmJ/rxhOO4vvlO1i+LbKWbfdWDRjcvjHlbqXM5abcpZS5rXeXJ80nr9ztpsx+tz4r5a7KaZ5jPGXFswr2GkSvIXOQGsD4eQ1iICNsH+P0M6ipTqmyX4pDghheK0+BHUUlXD4kv04YqPOB0/wM1EBVvSXQ/uG0oM5+bjpb9hziqQt787+FWzi3b2sG5OdW2mfhpj2sKdzPOX1a10jv9n3FlLncNKmXTkaqE1XlUJkLVXCr8tn8LezcX8KtIzrRfuyXAEy8fRhHN6tfo/P4UlruZveBUprUSyPFGX4PbpnLzT2fLOKWkzrRJjer1uevKS/8uIZxXy2P2/mCkeIQ2uRmVXGUaNUwk87N6/P98h1Bjowvfz27O33aNKRg1wHW7DjAPyetTKieAfmNvC278/u15sO5m+J27vQUByV+XbtHN6vHyu37o1J+29wsJv/xBJyO2LY6XG4fg+ZjACsMmWXgXO6KvHKXm7IAxs+7v/3uW54n3RXIoPodU+6unGcZYd9tt11OVcPr2zUcLusfOzOggTrcVtStUTy+rXuLGfPfhTTPyaB5gwya5WTQwt7OyUxl14FS5m/cA+B9qhs/qyJgutMhvHPtIC566WcAbn9/Ae/+fhDLthZx1dB8RPA2mYvLXLw+vYDOzetx9RtVjWKLnAyOadEgYEX3ns85T/3nFEb1aMElg9ryO1vT5D+ewKZfD/LSlLVcOKANI7u3wOkQtuw5RGaqk5Xbi7j+P3PZc7DyOEnBuFEhLmVlfl67iw/mbGLLnmL+c+2gsI+LlBT7zz/z7hE89vVyzunTig559fh+2XZystIoLnPx1Lcr+evZ3XGI8M2Sbbw3e2M1pVYlr346z17ch9nrdjOqZwv2HCrj6Gb16X7/NwDMGHsSTetnoKoUlZSzv7iclg0zq5Tz7swN3P2JNcZy6eC2ZKWlcNngdkxYtJUhHRqz6ddDpKU4+P1bc2jeIINt+yq6zX7btzV3nd4Zl1v5eN4m/jHRMi5/GNGJri3qc3r3qqvHuN2KYkVKdvhUlN1b5QBwXv/WFBaVcHSzenS97xtO7JzH5BXhRfFf+beRuFXZd6iM7PQUSsvdNMpOC7p/SbmLjbsPMXVVIQPyc70anpu8mpzMVC4d3I7rh3fgqjdms3H3IY7v1IRTuzXn3k+Dt5BvOvEojmnRgDkFv/LGjAIApv35RFo1zORAqYtte4vp0CSbnQdKaJydTkm5i2e/X81VQ/Np2iADgIOl5aQ5HWzdW0yb3Czyx0wIeK5Tuzbj++U7KPepQFs3yuTDG4aQ5nSQk5lao4e6aOF0CE6Hk/TDrTYOgttdYeA8LcVgBvXMZ6eFLOtwa0GlACuxIk1sBmYDl6jqkkD7p7fopG2vfpoyV9XvmJHqoLisbk82qYmBmrqqkMtencVxHZvE1UC9+OMaHv1qOUsfOo2stPD/of7jKP3bNWLO+ooxmjSng5/vHkFGqoP0FGfQp+Ad+4pJS3HQMCt4xRwNSsvdpKXEp/KbuXYX3yzZzmvT1/HB9UMYkF95TGz5tn0cKCmnX7vcEKVEn4Ol5bw/eyMndm5KfpPsKvkbdh1kz6FSerZuGNF59hWXsXVPMZv3HOSkLlXHJF1ujXmryBAeZS43m389RPu8eod/C0pVy0XkZuAbwAm8Fsw4AbTLzWLF30ZS5lJ2FBWzfV8xW/cWs22vtT1vwx5WbS9i2piT+GD2Rv42oeogdSxY9MCp/Gvyavq3y+X3bx3ZThyep9maVhgZqc6ABrimMew8T+GxJl7GCWBQh8YM6tCY+34T2MG1S/MGcdPiS1ZaClcd2z5oftvGWbQl8u7lBhmpNGieSufmgbvKjXFKHlKdjoAPKx4OKwMFoKpfYgWOrZYGmamICGkpQutGWbRuFPzHf+3xHbj2+A4B8zbuPkirhpk4HMK1b85m0jKrm+7OU47my8XbWLa1wgPsuzuHM+KJHwFY9+gZlLuVn9bs4vLXZgFQLz2F+hmpjB15DGC1cjbsOsiwxycDsPjB07zdTsF4/LyeNM/JYNLS7bz503pGdGnKq1cO8Obf8cF8ZibJwHh1ePqrUxzRqcCT1UvJYDDUnMNqHlSiaJOb5e3/z29sWft7zjiGW0Z04qtbj6+0rycfrMoy1elg2NF5vO5jQPxplF3hfpuVWtVF2J/z+rXm+E55DOrQGLCeQnxxHEaVtKcFZR5qDQaDP8ZARZlgFW2guSkVx1QcFI5t8bQS3Pb4oX/jwyEVecmOy+0mJYnnZxgMhsSREAMlIueLyBIRcYtIf7+8sXYYoxUicppP+mER4ihYRZsSoolQ2UCFX1F7bJD/MQ6Rw8hAmTEBg8EQmES1oBYD5wJTfBPtsEUXAd2A04HnRcRZF0IceSrhQF6TtW08eFtQfgWICLWYipAQPC0og8Fg8CchThKqugwCthZGA++pagmwTkRWAwPtvNWqutY+7j1736XxURw5oZwAajtm5LF1/vW7QwIbwmSk3K2V5vcYDAaDh2Qbg2oF+M7C3GSnBUs/bAjVjVXb+jlYC8oRRgvq0a+W8ePKQorLXHS//xu+WFA50KjLrcxdHzj2W6Q8MXEFH9kRB8pdalpQBoMhIDFrQYnIJKB5gKx7VPWzYIcFSFMCG9KAVbBfLL4wlMaHUMt317YF5faOQfmXV9VJ4lCpi09+2cxv+7Vi8eZ9vPjjWl78cS3f3zmc/SXlvD/Hsv9b9x7iiYkrePb71d5jrxyaT5nLjVuV8bM28tDoblwysG21s+7LXW72l5Tz7dLt/OmjhVXy//jhglp8a4PBcKQQMwOlqifX4rBQoYzCCnGkqi8BL4EVi68WGmKCdwwqQF6kY1BCgDEo23rN37iHi176yRs1wxOmx8NJ9pwtD2sKD1QyToA3BI2H+z5bwn2fLeG5S/ry92+Ws37XQRY+cCqrd+zn6UmrKi0DYTAYDLUl2Sbqfg68KyJPAi2BTsAsrJZVJxFpjxXi6CLgkoSprAWhurFq62KtGngO0crtRewrLg8akywScjJT2XvIivl307vzvOk9a7AOzzEtGrBs6z4yUh1kp6XwyLk9oq7TYDAc/iTEQInIOcCzQB4wQUTmq+ppqrpERD7Acn4oB25SVZd9TNghjmKJJ5hmg8zgl65n6xya+4XQ8bSg2gaJFD6ofS6XD8mvlPbsxX1okJnKFXYUitzsNHYfKPXmVzhJVLZQvx6sflXRgnGjuPbNOUxaVrH2zmc3HcsT367kifN70aRemncS7aEyFy9PWcuawv08/7t+qKo3+noomjfI4I+ndUZV6Z+fi1OEto3jFyndYDAc3hxWwWJrSiwWLCxzuflwziYuHNDGa3ROefJH1u08wGc3H0u3ljkBj3O7lce+Xs7p3ZvTp22jkOfocf83XDqkHX8+vQsAx477niuH5nPBgDYUl7loZhu/ouIybnxnHuN+25NWPpG3i4rLOPWfU0h1Orjr9M4MOzqPBhmp7NxfwrqdB6osJ7J9XzGFRSXe6NThUu5yk+J0sHjzXlKcQhs7lFR2XQnLbDAY4kKdWLCwpsRrRd39JeWkOoX0lOrDFBkMBoOhMkekgRKRQqC6JXWbADvjIKc6jI7k0gBGhz9GR3JpgLqjo52q5vkn1mkDFQ4iMieQ5TY6jmwNRofRkewajgQdyTZR12AwGAwGwBgog8FgMCQpxkDZk3qTAKOjgmTQAEaHP0ZHBcmgAeq4jiN+DMpgMBgMyYlpQRkMBoMhKTEGymAwGAxJSZ0zUCLSRkQmi8gye9XeW+30XBH5VkRW2e+N7HQRkWfslXoXikhfn7IeE5HF9uvCGOvoIiI/iUiJiPzRr6xaryYcZR2vicgOEVmcCA3BykmAjgwRmSUiC+xyHkyEDp/ynCLyi4h8kSgdIlIgIotEZL6I1Gh2fJR1NBSRj0RkuV3ekHjrEJHO9nXwvPaJyG0JuBa322UsFpHxIpIR6Jxx0HGrrWFJuNfBi6rWqRfQAuhrb9cHVmKtwvt3YIydPgZ4zN4+A/gKKyDtYGCmnT4K+BYrXmE2MAdoEEMdTYEBwMPAH33KcQJrgA5AGrAA6BpvHXbeMKAvsDjG9yTYtQhYTgJ0CFDP3k4FZgKDE3FP7Pw7gHeBLxJxX+y8AqBJnP6zoXS8CVxrb6cBDRN1X3z+v9uwJqLG8zfaClgHZNqfPwCuTMB/pTvWCupZWHXpJKBTuDrqXAtKVbeq6jx7uwhYhnWzRmP9eLHfz7a3RwNvqcXPQEMRaYF1M35U1XJVPYBlGE6PlQ5V3aGqswH/SK8DsVcTVtVSwLOacLx1oKpTgN3hnjvaGkKUE28dqqr77Y+p9itsb6No3hMRaY31MPVKuOePhY5IiJYOEWmA9RD1qr1fqaruibcOP0YAa1S1uog2sdCQAmSKSAqWgQi4RFGMdRwD/KyqB1W1HPgROCdcHXXOQPkiIvlAH6wn3GaquhWsi49l8SH4ar0LgJEikiUiTYATqbwmVbR1BCNqqwlHqCMqREuDXzlx12F3q80HdgDfqmpCdABPAXcB7tqcP4o6FJgoInPFWjQ0ETo6AIXA63aX5ysikp0AHb5cBIyPtwZV3Qz8A9gAbAX2qmr4a+JESQdW62mYiDQWkSysHquw69E6a6BEpB7wX+A2Vd0XatcAaWrfzC+BGVg/sJ+wlgCJlY4a6UuAjoiJloZIy4mGDlV1qWpvrMUzB4pI93jrEJEzgR2qOremx0ZTh82xqtoXGAncJCLDEqAjBasL+t+q2gc4gNUNFW8dnnLSgLOAD+OtwR4bGg20x1pbL1tELo23DlVdBjyGNVzyNdaDf9j1aJ00UCKSinVR31HVj+3k7XbXHfb7Djs96Cq+qvqwqvZW1VOwDMWqGOoIRqhVhuOpIyKipSFIOXHX4cHuQvqBGnT/RlHHscBZIlKA1fV7koj8JwE6UFXPf2YH8AlW13S8dWwCNvm0Zj/CMljx1uFhJDBPVbdXu2f0NZwMrFPVQlUtAz4GhiZAB6r6qqr2VdVhWEMEYdejdc5AiYhg9UEvU9UnfbI+B66wt68APvNJv1wsBmM1hbfaXTiN7TJ7Aj2BsJvItdARjNnYqwnbT2QX2WXEW0etiZaGEOXEW0eeiDS0tzOxKoPl8dahqmNVtbWq5mP9Lr5X1bCfkqN4PbJFpL5nGzgVq2snrjpUdRuwUUQ620kjsBY/jasOHy6mht17UdSwARhsD1EI1rVYlgAdiEhT+70tcC41uSZaC6+bZH4Bx2F1gS0E5tuvM4DGwHdY1vs7INfeX4DnsDzlFgH97fQMrB/3UuBnoHeMdTTHegLcB+yxtxvYeWdgedGsAe5JoI7xWP3ZZXb6NfHUEKyceF8LrIeVX+xyFgP3Jeqe+JR5AjX34ovW9eiA1XWzAFiS4N9obyyP24XAp0CjBOnIAnYBOQm8Fg9iPTgtBt4G0hOkYypWPboAGFGT62FCHRkMBoMhKalzXXwGg8FgqBsYA2UwGAyGpMQYKIPBYDAkJcZAGQwGgyEpMQbKYDAYDEmJMVAGg8FgSEqMgTIYDAZDUmIMlMFgMBiSEmOgDAaDwZCUGANlMBgMhqTEGCiDwWAwJCXGQBkMBoMhKTEGymAwGAxJiTFQBoPBYEhKjIEyGAwGQ1KSkmgBsaRJkyaan5+faBmGwxRVcKui3s8KCCIV+wggAm4Ft1u9eb7LrKlPeYq1EJvLbZXtVsXtttIdYpUtCA4BsT970j2FqU/Jan9WtRcf9Z5bcToclRSojxi1tSPWu+e7un3ysT+77S8jIt7vW+mz50J43wTvOnP29xGpuFZin1QClKlY3zfV6SlZK13PSqvXacU3q7gOFdsAGakOHL43zJCUzJ07d6eq5vmn12kDlZ+fz5w5cxItwxBjVJX9JeUcKnWhWBVquUspLnNxsNTFgZJyDpS6OFTmQoAUh+BwCMVlLg6UuNi5v4Sd+0soLCqxt0spLCphf0l5XL+HO8rlRUO9AM4geYqfwUhChh7TlFeuGJBoGYZqEJH1gdLrtIEy1H32HirjD+N/4ceVhRGV0yAjhbz66TSpl07Xlg3Iq5dO4+w0stJTLIMmkOJ0VLSqVHG5lXK39Z6R6iQzzYnbbbWEnA7xaWl4WgUOUhxCRqqT7HQn2ekpZKelkJ2eglOEEpeL0nI3JeVuSsrclJS7rO1yN8VlLgCcPmWLCCkOu1ynkGa/pzodlLncHCx1Wa0vHw1Q0SLzfJdyt+IUoV5GCukpDm/5KQ4HmalO0lOslpjLNvzlbjcut1Jmb3tajxWtMCUtxYGIUO5yU+ZSe3835W6l3PuulLndVpkut1fj/pJyCotKcCs4fFqQDrHug8MhiGfb5z0txUGq0/MSSsvd5GSmRvS7MCQWY6AMhy0bdx/kqjdms37XAW45qSPNGmR4KzOPIchMtQxBvfQUMtOsitZTOXoMRW52GukpwdoJ8SS5K1MHQqoTgrepDIboEtRAiUiDUAeq6r7oyzEYwmPu+l+57q05lLuVt64exJCjGidaksFgiDKhWlBLqBhLbQkU2dv1gM1A25irMxgC8MXCLdzxwQJa5GTw2pUDOCqvXqIlGQyGGBDUQKlqGwAReR74WlU/tz//BhgWH3kGQwWqyktT1vLoV8vp364RL13en9zstETLMhgMMSKceVADPcYJQFX/B5wYO0mHDy63smt/CQdLy3G7k92f6fDG5Vbu+2wJj361nFE9W/CfawcZ42Qw1HHCcZLYLSJjgP9gdfldCvwaU1WHCVv3HuK4xyZ7P2fanlyZqU6y0qzteukp5GanWZ5ReDySfDyQHJY31u79pRTsOkBpuZsm9dL569nd6dy8fuK+XBJxqNTFLeN/YdKy7Vw3rANjTu+Cw2HmthgMdZ1wDNQlwIPAV/bnKcDFMVN0GFE/I5UHz+rGQXuOzaHScp9taw7OrwdL2ba32DsJUrEnZvpMgnQr1M9IIb9xFtnpKfy8djejn5vG2JHHcNngdkd0ZbyjqJjfvzmHhZv38uBZ3bhiaH6iJRkMhjgh3hnfSYKIPAD8HvBMbLlbVb+088YC1wAu4A+q+k2osvr376+H40TdHUXF3PnBAqau2snxnZrwxPm9aNogI9Gy4s6KbUVc/cZsdh8o5emLenNqt+aJlmQwGGKAiMxV1f5V0qszUCLSEbgDyMenxaWqp0ZZo+d8DwD7VfUffuldgfHAQCyvwknA0arqClbW4WqgwGphvTNzA3+bsJSstBQeP68nI45plmhZcWPKykJuemcemWlOXr1iAD1a5yRaksFgiBHBDFQ4XXwfAa9ijUEFNQZxYDTwnqqWAOtEZDWWsfopgZpihohw6eB2DO6Qyy3j53PNm3O4fEg77j7jGDJS6/ZEyXdnbuDezxbTqWk9XrtyAC0bZiZaksFgSADhGCi3qj4bcyWVuVlELgfmAHeq6q9AK+Bnn3022WmVEJHrgOsA2rY9/KdqdWxan09vGsrfv17Bq9PWMX31Tm47+WjO7NkCqWNBMN1uZdzXy3lpylpO6JzHvy7pS710E+zEYDhSCcfN/DMRuU5E8kSkgecVyUlFZJKILA7wGg38GzgK6A1sBZ7wHBagqCr9k6r6kqr2V9X+eXlVguMelqSnOLn3zK68efVAr0fbnR8u4GBpfIOZxpJDpS5ufGceL01Zy2WD2/HK5f2NcTIYjnDCqQGutd/v9UlTIogkoaonh7OfiLwMfGF/3AS08cluDWyprYbDkeFH5zH1zyfxr+9X89R3K1m4aS/P/64vRzc7vN3RfT317j2zK1cfm1/nWocGg6HmVNuCUtU2AV4x6zsTkRY+H88BFtvbnwMXiUi6iLQHOgGzYqUjWXE6hFtP7sQ71wxiz8EyzvrXND6Ys5Fk88YMlxXbijjnuRms3L6fFy/txzXHtTfGyWAwAGEYKBHJFJExIvJv+3NHERkZQ01/F5FFIrIQK2LF7QCqugT4AFgKfA3cFMqDr64ztGMTvrz1OPq2bcRdHy3k5vG/UFhUkmhZNWLqqkLO+/cMylxuPrh+iHEjNxgMlQjHzXw8sAi4RFW7i0gWMF1V+8RDYCQczm7m4eJyK//+YTVPf7eKnMw0nrygF7bcKsoAACAASURBVMOOTv6xN+OpZzAYPARzMw/HSaKTqj4ClAGo6kECOywYEoDTIdx8Uie+uOV4GmWlcvlrs3h4wlJKypOzcel2K49+uYy7P1nEcR2b8OENQ4xxMhgMAQnHQJWKSAa2x5w9/lMaU1WGGtO5eX0+v/k4fjeoLS9PXcfof01n7vrkCpno8dR7ccpaLh3cllev6E/9jORepM9gMCSOcAzUQ1hjPq1F5E1gMjA2pqoMtSIzzcnD5/TgtSv7s+tAKee9MIMH/7ckKdzRdxQVc9FLP/HN0m38ZdQx/HV0d1Kc4fz8DAbDkUrIMSix3KmaA+XAUKyuvRmquiM+8iLjSBiDCsb+knIe/3o5b/60nnaNs3j8vF4MbJ+bEC0rtxdx1esmpp7BYAhMrcag1LJeX6hqoap+pqqfHi7G6UinXnoKD47uznvXDUYVLnjxJ+54fz57D5XFVcf01Tv57b9nUOpy8/71g41xMhgMYRNOH8ssEekbzZOKyPkiskRE3CLS3y9vrIisFpEVInKaT/rpdtpqe30qQxgM7tCYr287nhuGH8VnC7Yw8qkpzFi9My7n/uSXTVzx2ixa5GTwyY1D6dm6YVzOazAY6gZBDZSIeKJMHIdlpFaIyDwR+UVE5kV43sXAuVhrS/mesytwEdANOB14XkScIuIEngNGAl2Bi+19DWGQlZbCmJFd+O//DSU91cklr8xk7McL2XMwNr4uqsoLP67h9vcXMCA/lw9vGErrRlkxOZfBYKi7hAp1NAvoC5wd7ZOq6jIgUMSAYBHLAVar6lr7uPfsfZdGW1tdpnebhnz5h+N58tsVvDa9gB9WFPLPC3szuEPjqJ2jzOXm4QnLeGNGAWf2bMETF/QiPaVuR183GAyxIVQXnwCo6ppArxjpaQVs9PnsiVgeLL2qaCuw7RwRmVNYWBholyOazDQn94zqyqc3HktGqpOLX/6Zf3yzgjKXO+Ky95eUc/Ubs3ljRgFXH9ueZy7qY4yTwWCoNaFaUHkickewTFV9MlTBIjIJywPQn3tU9bNghwU6FYENaUD3Q1V9CXgJLC++UBqPZHq0zuGLW47jwf8t4V+TVzNt9U6euagPbRvXritu5/4Srn5jNku27OOx3/bgwgGH/1InBoMhsYQyUE6gHrWMGhFuxHI/QkUsP6IjmceC7PQU/n6eFRpp7MeLOOOZqfz17G6c06d1jcrZsOsgV7w+i617D/HSZf2OqJV/DQZD7AhloLaq6kNxU2LxOfCuiDyJtay7J2K5AJ3sKBabsRwpLomztjrLmT1b0rtNQ25/fz63v7+AH1YU8tDo7uRkVh/lYfHmvVz5+mzK3W7euXYQ/dolZq6VwWCoe4QyUDGLtyci5wDPAnnABBGZr6qnqeoSEfFELC/HJ2K5iNwMfIPVsnvNjm5uiBKtG2Ux/veDef6HNTz93Spmr9vNuN/2pGXDTH49WMruAxWvXz3bB0uZvW43DbPSeO/qQXRsenivS2UwGJKLoJEkRCRXVXfHWU9UOZIjSUTC/I17uP39+azbeSBgflaak0ZZaeRmp9E2N4t7z+xK85yMOKs0GAx1hWCRJIK2oA5342SoPb3bNGTCH45jwsKtpKU4yM1O8xqk3Ow0MlKNZ57BYIg91a4HdTgjIoXA+mp2awLEJ7RCaIyO5NIARoc/RkdyaYC6o6OdqlZZyK5OG6hwEJE5gZqWRseRrcHoMDqSXcORoMOsd2AwGAyGpMQYKIPBYDAkJcZA2VEnkgCjo4Jk0ABGhz9GRwXJoAHquI4jfgzKYDAYDMmJaUEZDAaDISkxBspgMBgMSUmdM1Ai0kZEJovIMnvV3lvt9FwR+VZEVtnvjex0EZFn7JV6F/quHiwij4nIYvt1YYx1dBGRn0SkRET+6FdWrVcTjrKO10Rkh4gsToSGYOUkQEeGiMwSkQV2OQ8mQodPeU6xFhL9IlE6RKRARBaJyHwRqVH4lijraCgiH4nIcru8IfHWISKd7evgee0TkdsScC1ut8tYLCLjRSTscC9R1nGrrWFJuNfBi6rWqRfQAuhrb9cHVmKtwvt3YIydPgZ4zN4+A/gKK/bgYGCmnT4K+BYr2kY2MAdoEEMdTYEBwMPAH33KcQJrgA5AGrAA6BpvHXbeMKxFLBfH+J4EuxYBy0mADgHq2dupwExgcCLuiZ1/B/Au8EUi7oudVwA0idN/NpSON4Fr7e00oGGi7ovP/3cb1kTUeP5GWwHrgEz78wfAlQn4r3THWkE9C6sunQR0CldHnWtBqepWVZ1nbxcBy7Bu1misHy/2u2el4NHAW2rxM9BQRFpg3YwfVbVcVQ9gGYbTY6VDVXeo6mygzK+ogdirCatqKeBZTTjeOlDVKUCNQ2BFS0OIcuKtQ1V1v/0x1X6F7W0UzXsiIq2xHqZeCff8sdARCdHSISINsB6iXrX3K1XVPfHW4ccIYI2qVhfRJhYaUoBMEUnBMhBhL1EURR3HAD+r6kFVLQd+BM4JV0edM1C+iEg+0AfrCbeZqm4F6+JjWXwIvlrvAmCkiGSJSBPgRCqvSRVtHcEIezXhGOuICtHS4FdO3HXY3WrzgR3At6qaEB3AU8BdQERLIkdBhwITRWSuiFyXIB0dgELgdbvL8xURyU6ADl8uAsbHW4Oqbgb+AWwAtgJ7VXVivHVgtZ6GiUhjEcnC6rEKux6tswZKROoB/wVuU9V9oXYNkKb2zfwSmIH1A/sJawmQWOmokb4E6IiYaGmItJxo6FBVl6r2xlo8c6CIdI+3DhE5E9ihqnNremw0ddgcq6p9gZHATSIyLAE6UrC6oP+tqn2AA1jdUPHW4SknDTgL+DDeGuyxodFAe6y19bJF5NJ461DVZcBjWMMlX2M9+Iddj9ZJAyUiqVgX9R1V/dhO3m533WG/77DTg67iq6oPq2pvVT0Fy1CsiqGOYIRaZTieOiIiWhqClBN3HR7sLqQfqEH3bxR1HAucJSIFWF2/J4nIfxKgA1X1/Gd2AJ9gdU3HW8cmYJNPa/YjLIMVbx0eRgLzVHV7AjScDKxT1UJVLQM+BoYmQAeq+qqq9lXVYVhDBGHXo3XOQImIYPVBL1PVJ32yPgeusLevAD7zSb9cLAZjNYW32l04je0yewI9gbCbyLXQEYzZ2KsJ209kF9llxFtHrYmWhhDlxFtHnog0tLczsSqD5fHWoapjVbW1quZj/S6+V9Wwn5KjeD2yRaS+Zxs4FatrJ646VHUbsFFEOttJI7AWP42rDh8upobde1HUsAEYbA9RCNa1WJYAHYhIU/u9LXAuNbkmWguvm2R+AcdhdYEtBObbrzOAxsB3WNb7OyDX3l+A57A85RYB/e30DKwf91LgZ6B3jHU0x3oC3Afssbcb2HlnYHnRrAHuSaCO8Vj92WV2+jXx1BCsnHhfC6yHlV/schYD9yXqnviUeQI19+KL1vXogNV1swBYkuDfaG8sj9uFwKdAowTpyAJ2ATkJvBYPYj04LQbeBtITpGMqVj26ABhRk+thQh0ZDAaDISmpc118BoPBYKgbGANlMBgMhqTEGCiDwWAwJCXGQBkMBoMhKTEGymAwGAxJiTFQBoPBYEhKjIEyGAwGQ1JiDJTBYDAYkhJjoAwGg8GQlBgDZTAYDIakxBgog8FgMCQlxkAZDAaDISlJOgMlIm1EZLKILBORJSJyq53+gIhsFpH59uuMRGs1GAwGQ+yoVTRzEQm5CJjaa9nXSpC1CFYLVZ1nrzEzF2vd+wuA/ar6j9qWbTAYDIbDh5RaHvdEiDwFTqpluai1zr1nzfsiEVkGtKpteQaDwWA4PEnq9aBEJB+YAnQH7gCuxFoQaw5wp6r+GuCY64DrALKzs/t16dIlTmprT5nLTaoz8b2tyaJDgfIk0KIKLlVSHJJQHcmkxeVWRMAhidWRLL9VQ3SYO3fuTlXN80+P2ECJSHegK9YKtACo6lsRFWqVWw/4EXhYVT8WkWbATqz6669Y3YBXhyqjf//+OmfOnEilxJRnvlvFk9+u5Mc/nUC7xtkJ07G2cD8nPfEjfzqtMzed2DFhOgAe+HwJb8woYN69p5CbnZYwHRe++BMz1+2mYNyohGnwcO7z05m3YU/CteSPmYAIrHs0cTq27j3EkEe/5/9OOIo/n578D6CG6hGRuara3z89okcQEbkfeNZ+nQj8HTgrkjLtclOB/wLvqOrHAKq6XVVdquoGXgYGRnqeZGDa6p0AbN1bnFAdm/ccAuCnNbsSqgPg26XbAThQUp5QHTPX7U7o+X2Zt2FPoiV4SXSny86iUgCmripMrBBDzIm0jXweMALYpqpXAb2A9EgKFBEBXgWWqeqTPuktfHY7B1gcyXkMgVES3+XrclsaUpyJ71ozGAyJo7ZOEh4OqapbRMpFpAGwA+gQYZnHApcBi0Rkvp12N3CxiPTG6uIrAK6P8DyGJKXc7QbAmQRjPwaDIXFE2oKaIyINsbrc5gLzgFmRFKiq01RVVLWnqva2X1+q6mWq2sNOP8v29juseHjCUh74fEnAvHh2mzw1aSV/+nBB/E4Ygrs/WcTj3yyvlFbmsltQjvgMgqsqo/81ja8XJ/4n9c7M9Vz5ekR/oaiw91AZwx+fzOLNexOq45cNv3LiP34I2N2b6K5GQ+yJqAZQ1RtVdY+qvgCcAlxhd/UZAvDy1HW8MaOAcpc7oTqemrSKD+duYu+hMm+akJjWyrszN/Dc5DW43RW1jef6xEtRqcvNgk17ueE/tZ6+FzXu+WQxP6woZMe+xI5Jzly7i/W7DnLl67MTquOxr5ezbucBvli4JaE6DIkh4kdUEWklIkOBtkBDERkWuay6zcNfLku0BADOf2FGoiV4eWHKGu92uTtxj8aFRSVV0twJ0DPwke/ifk5fPG7kO/dXvR6J0PHn/y5KqA5DYojUi+8xYDrwF+BP9uuPUdAVNy55+Wcuf20Wswt28/CEpcRjXtjr0wuqpF388s+M/XghP6zYwb9/WFP1oBiwcvv+KmnTV+/i1H/+yJIte7n7k0WUlLviouXvX6/wbnsMVJ+/fsvzP6zm8wVb+HDOxpid2/eW3/1J1Yqww91f8uq0dXw8bxOf/rI5ZjrCIX/MBN6cUcBHczfx+YLYtSp8e1cnL98RUMesdbt5YuIKFmyMnYeh73wr///mki37GP74ZJZt3cfYjxdSXBaf36ohfkTqJHE20FlVE/uYVUuKisuYYbtVT1lpuaxedWx7WjbMTIie8bM2Mn6WVRHfMLwDkqDJkCu372fUM9MA6N2mIRf0bxO3c6uq14sPKhuu82Okw7fe87i4+/PXL5Z6t8/uE5/AJgU7D5DfpOrcuPt9xjHP6tUyJuf27fJ9eepaTuzStMo+F7z4EwDPfr86ZvOzfL1Kl20t4pgW9Svlr991kJFPTwWgXeNsbhh+VEx0GBJDpAZqLZAKJKWBWrR5LwMenkTHvHq0aJhBi5wMWuRk0rKh9b5x98Eqx1z1+mw+vnEoxWUutu8roWvLBjHRlj9mQsj8Z75bzU0nHsX2ohJKylx0yKsXEx2DH/mObSHGO+76aCEjujSlcb109h4sY+OvB+neKicmWqq7Jj+s2MHwo/NYU3iAeukpNM/JCLl/uLj9nsyr0zF99U6GHtWY1Tv2k5OZStMG0dHhzwn/+KHafWas3snQjk1Ysa2IRtmpNK0fJS0+z0Yz1uyq9pos3ryX7q1ymF2wm56tc0hPcUZFhu/DyhnPTA2577ivlvPbvq3JzU5j1rrdDDmqcVQ0GBJHpAbqIDBfRL7Dx0ip6h8iLDdqFBaV0C43i5lrd7NtX3GlH3wgVmwv4sxnp3Go1MW2fcUJm7n/z0krWV24n//Z3Tix0hHKOHk489lpvHfdYK5/ey7LtxWx9pEzcCTABfzK12fz19HduPczqwURrWvib6Cq43evzOSpC3tz2/vzcTqENY8kLrD+Ja/M5LlL+nLTu/PISHWw/K8jo1JuTe/umc9O46MbhnD+Cz9xQf/W/P28XlHRUd3/1Z8BD0/itpM78dSkVbxz7SCO7dgkKjoMiSFSJ4nPscIOzcByM/e8koKGmalMumM4H/3fUKaPOYmVfxvJz2NH8PGNQ3nukuAB2dftPOCtuMPt1161vYj8MROYUxC96AP/q8UYw5Y9h8gfM4HJK6qOG9SWrXuLGf74DyzfVgRQyfsvFBt3HyR/zISozvj3GKea8Nn8zeSPmcD2IMa4Nj4Qt71vTdGrSQW6bW8x+WMm8P3ywN2IteWmdy3vw+Ky8L1Dx321nPwxE4L+vmsTa++8F6wuv+mrw49GctITP3DWv6YFza/NvXlq0irA+h+HQ0m5i/wxE3jmu1U1P5khpkTqZv4mMJ4Kw/SunZYUtMnNomPTiq4xp0NonpNB37aNGNWzBd/fOdyb95sgffm7DpSGda5Rz1p/sgtf+rnafc/t24oePt1k487tUVFOjxaBDgnbk+zW934BrK7K6ji1azM6Na1HYzve3ePn9fTmXdC/ddDjdh0Ir0f3+retZ5XLXq1+Xs/5/VrTq3XFNbl4oDXeJAID2+eGdb5g3PqeZUzenx3Y0cJ38P30bs0Z2b2597NnOyvNWeme1YY7P7R0XP1G9fEhz+3TitO7NefigW29aSd0tmJpNshIoUvz+sEODYsXfrQccVZuLwqY72ufjuvYhCEdqnaXnd078H/GEzYrHNYWHmDhpuBzrXwfADo2rUd2mtV1+KjPf8b3GvkSbqisb5ZYDwxPfrsyrP0N8SOiLj4ROQF4EyuygwBtROQKVZ0SubSg5zwdeBpwAq+o6rjaluU7rvPsxX0CtliOHfc97Ztk8+mNx5KTlRq0rNJy6+k1nCfqh8/uQWaak1XbixARWjfKZMzHi7w6JiyqOmG0w91fkpXm5KcxI0LqmF1QJcB7UMaM7EKHvHrs2l/CmsIDDGyfy58+WgjA38/rxQdzNgEw7Og8rxMJwMlPTiEnM5UpfzoxpJalW/eFreXx860uoV82/EqLnEw27znE+Fkb6dW6IY+e24MRT/xY5Zj8MRN48KxunNu3FfUzAuvw3BeAzb8Grjh9e/heuKwfYI3tdGuVw7z1v/LV4m0MyM/l1pM7ce7zVV3z88dM4Inze3Fy12bkZAa/HjVpWfxucFv6tbMM80ldmnJcxyb8uHIHP6woZHCHxlx1bHsufrnqw1D+mAk8dWFvTu7ajHrpgf/evq70wR7AfJ0k/nPtIFSV75btYHjnPH738kxmFezm4oFtKdh1kPkBvPjyx0xg6l0nklc/nYzUwONRM9dWXI9yl5uUANHJPd2vJ3bO4/WrBrLnYCnLthYx5KjGjLX/M4+e24PxszZUOfbRr5bz6FfLWfTAqWSmOgOWD/BFDL0hDZERaRffE8CpqjpcVYcBpwH/jFxWYETECTwHjMSKoH6xiHSN1fk8rNt5gF4PTeTtnwoqpRcWlXDjO3MDuiaHwhNjrlOz+nRsWq9Sd0qosZ2DpS56PTSRv/l4lAHsKy7j+rfn8GwNuyjSUqzb37heeshWys0BopvvPVRGr4cm8sTEFZXTD5Zxw9tzq6SHS5+2jWiek+Ft1VQ31HX/50vo8cDEKl1Vv2z4lb98uoij//KVN+39IK7qgcaghnZsQk5mqjevurBLd364gF4PTmSbX9DfIvvePDd5degv4oevB+cpXZuRmeb0dndVp+W29+fT/f5vKPObEP714m08+tUyBjw8yZt2w9uBe+T9TyEinNy1GalOBy7PvXGIt/UdiOP/Ppku935dJf1f36/ilalrK/U2fDR3U8AyPA98nu/cMCutxs4PPR6YSMd7vqqUpqrc/v58Pv1lMxN9PDcPlRpX9WQiUieJVFX11kSqutKORB4rBgKrVXUtgIi8B4wGloY8Kkrc+9kS7v1sCesePYPpq3dx6asza1WO069/v6Yx516Zto5Xpq1j3aNnsHDTXkY/Nx2o6KoIl7Qw19MJZbye/X41z36/mnWPnsG8DXv47b+tFsbXNR8qqkS4lbGHLvd+zbl9WvHEBb14bvJq/jGxandNsNZNqEavp4J0SEVLKyvNycEgFdngR60JtusePYMlW/Zxpt31W9N7E+hbV2gRr/t1TmZq0DHBTvd8xdXHtufeM4/hpnfn8eWibVX2GXFMVfdxIOQUB7dW6PDs16NVDouChEXKHzOB8b8fzKD2uXS858uA17tREEPn2Tca60/lj5nAkgdPo9yl9HpoIgCf+M1rc5n4SUlFpAZqjoi8Crxtf76U2DpJtAJ8H4M3AYN8d/BdsLBt28B905HSfuyXQfPCGRvwbyXV1iEulI5QT7YePC2oaBBKS8tauIN7KmPfijK/cRYFu6pODfDw8S+b+TjERNoTO1dZDw2oOgHUF9/K2EPb3Cyvw0gwQl2PBhm1+9u5fVouHlo3ygzptPLa9HW8Nn1d0PxOTQP/XkPZA7e7auu2fjXfKVB3pC/pQX6Lbnd4Ldhw6Xb/NyHza+o1aIgtkdZQ/wcsAf4A3Gpv3xCpqBAE+pVW+kWp6kuq2l9V++flBa6QYkltHsBiMSE3HDfweK1IWhuX9EBdfJFep2B1T6g6KdATfKRP8+Fcj0DncAe4JpFW3LWJnBLomkTa8Ah2vCeyfbxW8E3mFcaPRCL14itR1SdV9VzgGuC7GEeV2AT4hhNoDUQ0wvmXUcfwyuVVFnL0klrNmkRN6lVe/sq/i+CjuZv4enHVrhV/bhh+FB/fODRofkZq6FvlPyDuX299uWgrn/xSuZ8/UAvqifN78Xcfb77a4N+V5v+f/9+CLXw2P3TIoAHtc7mwfxse95lP438nbh3RKWQZLfxabv735p/frmTJlr0h18AacUxTzuvXmgfO6uZN8w+y/qfTOofUkZVW2UnAv4t3xpqdvDatcisnUH08snsLftu3NfeMOsZnv8o71nSFWV/jrKrc88kiduwrDjk37JmL+3DJoLYhJ2wHGrcMraPifPuKy/jThwsod7krjGEAQ/z0Rb15+JzuIcutqV0zDaj48u3S7XwQxLsWIo/F94OINBCRXGA+8LqIPFndcREwG+gkIu1FJA24CGsuVq259vgOnNy1WVj7BpoYOucvJ1f67PtH27r3EH/8cAE3/Gcuew+Gnjs0ZmQX+rZtFDS/3FVRbsG4UVW68BY/eJqfjopty5ljHre/v4AdRRWD+CkB/vS/7de6RqGNCsaNqmIcF9x/qp+WCjE7ioq5Zfwv3Pre/Epa/El1OnjsvJ60yc0KsU+F/rUBJst+5zONACo/Hf+0ZhdPf7eKUc9MC7macXqKk3+c36tS1Ar/p3nPWF6KQygYN4rmfpEllj50etDyd+4v4ZKXZ/LQF0vZ4NN92bpR1e+dkerkiQt6VYoW4X8LfbvKfrn3lEp5vVrnMOEPx1VK8703475azjszNzDwke9Yb2vp2qJqJJX2TbJ55JweIVtvDW3vziuH5vPm1ZUXv376ot6c16/yNAbf32vPByby4dxNnPnsNO89/k3PqtMvRvduxe8GtQuqAaBljhW2bNIdw7l+WOWl6tY9WvU347kepeVuPp63KSGBgo8kfv/WHO7678Kg+ZGOQeWo6j4RuRZ4XVXvF5HgZ4sQVS0XkZuBb7DczF9T1QiH40NjudtWNg5PTlzBM9+vDji2ogrP/7CaJVv2VQqi6RmU7du2Ya10+Ef4nnvvKbw3a4PXPb2qDuXtn9fz/bLt7PEZoxj4sDWI36tNw6h1LS5+8DTenFFQKUacL25V3pi+jmmrd1aKXOHRcmaAyiccfJ+qHbZxGP74ZNbvOsjfzu5epQJ1uZWb351Hl+b1KzlReFzHHw+z5eh/3Tzn8dyhn+8ewUdzN/HHIGtuuVV5Z+Z6Ji7ZXmnMY9jjk73buWGMIUJVY+m7CnGj7DQKxo3yhil64oJeVX5HB0tdXPDCT1w6pB0vTlnrTfe4cD80uhuRIALDj86rpOOUrs2Yvnpnpf0Kdh7gzGenctuIo71pvuN8I44J7yEyGCkOYewZx3D98KPo+9dvbW1Vf//TV+/klanrGHZ0E56bvIYnJq5kdO+WOB2CQwSnQ3y2CZAmOEVwOPzybYcSZ4B0R8DjqZxv5zkcgkOolO453nucvW+iYnlGk0gNVIq9FPsFwD1R0FMtqvolEHwEOkKa1EsPe4mBCwdUdcJYt/NApQCn/hwKc7b/8Z2aMHXVzup3BC4M0OLZub+Uez9dHFxHafWTGPPqpwdcfiIUlw+p+kS7fV8JD/wvuKNlMI84D7lZVmU9skdznptcEendv6vMl+M6NiHVry/O40n3xcLACxNWN87RsqH1QHJq12aVHj5CDQOd27dqYNlfD5ZxzyfB7004tLFbWKd0bcbc9RVz3zyVUqCv0rFp/SoTc9+YUQDArCARUKqr407onMekZdvpkJfNT2urzvMKdE2z0lKqpHuWoLn2rcCTmKsbagvmQKN+0wQ85QTqPYCKSd0ej8TNew7x0pS1uFQPuwUSxc+QOT3GzVFhwCoMX6B9K4yi0z4ukEF2CH77VjbilfetON5zB3Kz01gf5DtEaqAewmrNTFPV2SLSAThs44Usfeg0HCIB525Ei3CWr1j0wKlkpDrp5Dd3I7o6qjeUU+86scZx6mpDdXNPGmWnseC+U6mfkVLZQFVTa9XUOaO6BXxb5GQy/75TyMlM5fFvKh5CvC2oONZgbXKzmH/fKTTISGXcVxUrEnuMdrBvXlOfiuqewn83qC2jerSo4q7t+d0EO7qmT/fV7f/1bcNwubWKl57/dIWannfhA6fSwJ4E7om071LF7bbGNF1uxe1NU580n3yf94Dpvsep4nJT5TiXu3K+W4OnV07z0eabb6dVOc67r+/xeNNcqpSWu/2+LxVlqe/1oHKZVa6XNT776hX96Xtf4OsfkYFS1Q+BD30+rwV+G0mZiSQrzbocd53emTU7DvDfeYEnD9aEs3q1rLRuTzh92p6ondbAFQAAIABJREFUCJcPaUdudpo3tlgkHNuxcaVIBuHUpb4RAG45qSOHSl284jeYXxtO6JzHDysqIlOEM/fEE7HipC5NmbfhV/YcLIu4C+Oa49rzqs/3CWdV4YZ2a+64jk2YXbCbknK31xD6fouaGKvBHXL5eW3NYzh6tAzukMucgl8pdysex8xg16a6VuIVQ9rx5k8Vz7PV7S8iNMpO4ze9WvLOzPWsKbTi33m+frDDPTqPysv2HuPLKV2bBV36JBCe3+qzF/fhgzkbvb0PHgcYj45gBrpJvTR27q8aVaO+z/iqiJDilIif6g3hE6mTRIaI3CQiz4vIa55XtMQlihtP6Mgj54b2DgpEv3aVnRw+vGEIz1zch7+dXVFWTVaLfWh0d247+ejqd/TjFD+nj7euHsg71w7m7jMqPLxq2jK689TO/OXMmgft8J979M61g3jjqoGVPN9qMvfktSsHcHZvq9ssEvM0fcxJ/GXUMXTzWU6lJvbuP9cOYpQ9dhaqEvc3eqd3a17p85tXD+S964bw798FD15cHe9dN4TT7JiBHsMUrCL2aG3XuKoTxtKHTuNuHw/BUOX4k1c/ne/uPMH72XNHg10bT/qwo6tOBSkYN4q/jq75/w+smJpvX1MxNdLzM3d4r0vgLtDrh1VeR+q6YR0oGDeqTozjHM5EOhHmbaA5VoijH7HcvkPPYDxM8HSVBAsiG4h3fz+okgfbgHwrAsOlg9vx3/+zXMhrOxHwpAALxgXjuUv6suC+Ch2eSuC6YUfx6hWWS30kvVH+BjAUL17Wv9I18Sx/cNOJHXnL9u6q6RLvnsgW3WqwVteSB09j9j0VHpetGmYiIkz4w/E1use+eAKoeiZnh3NNn7m4T6V7M9y+NyN7tGDSHcODHVYtx9nXtZMdHDlYa9BTQbvcyoL7T+UTe2pD29wsstJSSE9xUjBulHfibU3nH3l+a94HoGoMZZtGWcz5y8lcYY9dPvAb6yGoeU5GxMuppDrF28XnkRHs+1x1bH4lj9wrhuZHdG5DdIi0tdpRVc8XkdGq+qaIvIs1JnXYk+J08Mu9p1AvIyXsrr70FGfQhdo87rI1rYwBFtx3KplpTrrf/w2lrurHjtJSHEGjRES6jtP8+04hKy2lUoy72mrxeJy53OEvEwFwRo8WzL7nZPLqW3PQqjNUR+Vlk52eEnTMyuP4VtOHh/P6teaEzk1pUs/qahsUIiSU0yG43BryelQ35y4UFw1ow8nHNPO6dw/qEFhLg0zrLz+4Q2NyMlO93dqB9ELN5xG9cnl/DpW6mL7G6mI7KshCm5571j4vmyb10kPMRKsdSx48DRG44/0FfL1kGxn2XDT/7+PpbnY6pNKcxtQErHdmqEqkBsrjv7xHRLoD24D8CMtMGjzxwWbdPaKSU4H/n2nOX06u9unZE7WhNi0oz/jLjLEnsb+4wvvOv6Rf7j2lWgPmaRnWdkDfM+4x+56TOejjCehf3rx7T6G8Gi2e+UO+c7zCxWOcvr19GM183P39v9a0P5/onTgczHPLaXtH1PThQUS8Or6+7fhKc5eq/EbuObnaexNJZA9/LW0CzKMC6/59d+dwWjey5gf5e7d5SHFU7hILF48BHtm9Of/9v6FBp1Wc1681PVrn0KW5Zai8XXFRMgzZ9tjRPy/szS0793sdHfy/zr9/148dRcVVuvLiFWXFEJpIDdRLItIIuBdrwmw9IIg/xuGL/5LengFZT3QH/2gSgYjEQHloUi+90rk8FXym/XQYLOCmL/5zdmqLVRlWaEm3r0mm/R7OPB7PNYnEU7BTs8qx5DxdU57v6Ws0grWgPJVxJJMyPRWtB8+E2awa3JtoVYr+WvzxbdV4LkmwuV21HYIRkSpjsv75vjq9zgy1O11QMtOcdGtZEfHC0/XpeWjJTHPSrnF2leNSImjNGqJHpF58r9ibPwIdQu1bl7jmuPaUlLu58tj8sI/xVAThrkYbDqN7t2TD7oP8flj4l97zRBwqekJtOK9fa7buLeaG4eFr8bRYVm7fHzUdL1/eny8WbgkYgSLYgPeqHdawaU0W2quOUT1asLbwQI3uTSRdfKH44pbjWBZkba5gJnn7Pmv+W21at8F46bJ+QdeoipeHflqKgwfP6uYd+wuGaUElB5EuWNgMeARoqaoj7bWZhqjqq7Us73HgN0ApsAa4SlX3iEg+sAzwTD75WVVjGZQ2JBmpTu44pWbedU3qV9/KqikpTge311CHZ7JptEl1Omp8TVo1zIy6jpYNM7nOzyOrOjo3b8C8DXui+vRem3uTHaTyjpTurXKCxs0L5suQ4hDK3dGdnHqqnwdjJR2ejTh4zYXjABGsO9gQXyL9R7wBvE5FFImVwPtArQwU8C0w1g5p9BgwFviznbdGVXtHoDWhNMhI5atbj6dZg9gYiHBp1zibb28fFnIl3HjRPCeDL245Luh4Sbx49NwenNWrZcRLy0dKqtPBpDuG13o5jtrgsQf+LYYVfxvJ7ILd9Ggd2TL34eIxCP7RQX6595Sot/bDIVrLexgiI9J/QhNV/UBExoI3Vl6tl6RU1Yk+H38Gzvv/9s49SqrqzvefH03TNi2keRrFKHIlGBMZQIJkcHGTEOMro+g81BsnZt1knFnjzEicMAtH10TvjQudZLKc5M7KxGgyZjSG+AiaxKj4iMma4ANsnooDKCIt2sTYwQARaL/3j70Lqquruutx6pyi+X3WqlWn9jm161v7PH778du/XaO+TPjXi6cVvcA/UCTwZj255c9PZefv+4Y0Khy3SYN/v3QGe/b1vTT6i4hdD/7xnJOKjo1UukprrXz7MzPp3t13YuiJ44t7vdWLSWPb+Mu5k7hkVu+wXU1DjNmT0iuTq854P1Lf0FCj2oaVNX6XFD/529N5YkOXz39qEKyW8Cxm9nNC5IhlkmaY2WzgJknVT+g4mPePgSWS7ohdfOsJLbSdwLWSflnie/kLFp76yiulojw5juM4jYCZrZTUZ92jWg3UDOAbwIeAdcA44E8klYxobmaPEib3FnKNpPvjMdcAM4ELJcnMWoAjJb1pZqcCS4EPSio+8nvwt3ZAyTiEOcYC5UVlrS+uo7E0gOsoxHU0lgYYPDqOl9THc6UqA2VmHwZelfS6mQ0F/pLQknoe+CdJlQcWO5j3ZYRVeedJKrq+d2y5fVFS8dDHlf3eimKWO21cR2NpcB2uo9E1HA46qvWl/BbB0w7gDwlOEv8GvAXcUq0YMzuL4BRxXr5xMrNxZtYUtycBk4GXiufiOI7jDAaqdZJoymslXQTcIule4F4zW1WDnv9HmP25LA5S5tzJ5wL/x8z2Az3AX9XSSnMcx3Ean6oNlJkNlbQfmEd0SqgxTySdWCL9XuDeavMdgKpbfAnjOg7SCBrAdRTiOg7SCBpgkOuodgzqGuAcwqDYccCM6MxwInC7pDnJynQcx3EON6r24osu5UcDj0jaFdPeT/C2ey45iY7jOM5hiaRB9QLeBzxBCI20Hrgypo8mRKrYGN9HxXQDvg5sAtYQWoO5vG4iuM+vAy6qs46TgOXAOwQPxfy8ziKEedoELMpQx3eALmBdFhpK5ZOBjiOAZ4DVMZ/rszoncX8T0AH8JMNrYwuwFlgFrMhQRztwD7Ah5veRDK6PKbEccq+dwIIMyuILMY91wF3AERmdkyujhvXllsOB71Zy8KHwIrTqZsTtEYTJvScD/0x8uAOLCBOKIXRV/oxgqGYDT8f0c+MJGAq0ASuAkXXUMR74MHBDwcXeRIhLOAkYRngonpy2jrhvLjCDyg1UUmVRNJ8MdBihpwCgGXgamJ3FOYn7rwK+T+UGKslrYwshskwa92x/Om4HPh+3hwHtWZ2XvPv3dcI8nzSv0QnAy0Br/PxD4LMZ3Cu5ObLDCc/SR4HJ5eoYdCF7JW1X7GKU9DahBjABOJ9w8RLf58ft84HvKfAU0G5mRxNOxpOS9it0Ya4mtGTqokNSl6RnObjGVo5ZwCZJL0naC/wg5pG2DiT9AqjYezIpDf3kk7YOScqFYG+Or7L7ypM8J2Z2LKEydWvhvjR11EJSOsxsJKESdVs8bq+k7rR1FDCPEEe0rJA2CWsYCrTGuarDgdfK0ZCwjg8QvLF3KzjVPQlcUK6OQWeg8okhkqYTarhHSdoOofAJFh9Cob+a97VtMW01cLaZDTezscDHCM3eeukoRSl9aetIhKQ0FOSTug4za4pTKroIob4y0QHcDPwDUNmyxMnrEPCIma2M4cay0DEJ2AF818w6zOxWM+u72FP9deRzMaF7LVUNkjqBrwJbge3Ab9U71mkqOgitp7lmNsbMhhN6rMp+jg5aA2VmRxJc0xeo/5BIxaJCKp7MB4FfES6w5UDfyKvJ6ahIXwY6aiYpDbXmk4QOST0K0fWPBWbFFaVT1WFmnwK6JK2s9LtJ6ojMkTQDOBu4wszmZqBjKKEL+puSpgO7CN1QaevI5TMMOA+4O20NcSHZ84ETgGOANjO7NG0dkl4gjOUvAx4iVPzLfo4OSgNlZs2EQr1T0n0x+Y3YdUd874rp2+ht0Y8lNoUl3SBpmqQzCIZiYx11lKKkvpR11ERSGkrkk7qOHLEL6edU0P2boI45wHlmtoXQ9ftxM7sjAx1Iyt0zXcCPCF3TaevYBmzLa83eQzBYaevIcTbwnKQ3MtDwCeBlSTsk7QPuI0T9SVsHkm6TNEPSXMIQQdnP0UFnoCyEoLgNeEHS1/J2PQBcFrcvA+7PS/+MBWYTmsLbYxfOmJjnVGAqUHYTuQodpXgWmGxmJ8Qa2cUxj7R1VE1SGvrJJ20d48ysPW63Eh4GG9LWIelqScdKmki4Lh6XVHYtOcHyaDOzEblt4JOErp1UdUh6HXjVzKbEpHmE+KCp6sjjEirs3ktQw1ZgdhyiMEJZvJCBDsxsfHw/DriQSspEVXjdNPILOJ3QBbaGg26e5wBjgMcI1vsxYHQ83ghxBDcT3GRnxvQjCBf384S1qabVWcd7CTXAnUB33B4Z951D8KLZTIj6npWOuwj92fti+ufS1FAqn7TLglBZ6Yj5rCMESM7knOTl+VEq9+JLqjwmEbpucm73WV6j0wget2sIqx6MykjHcOBN4D0ZlsX1hIrTOuA/gZaMdPyS8BxdTQgCXnZ51LTchuM4juPUi0HXxec4juMMDtxAOY7jOA2JGyjHcRynIXED5TiO4zQkbqAcx3GchsQNlOM4jtOQuIFyHMdxGhI3UI7jOE5D4gbKcRzHaUjcQDmO4zgNiRsox3EcpyFxA+U4juM0JG6gHMdxnIbEDZTjOI7TkLiBchzHcRqSoVkLqCdjx47VxIkTs5bhOAPSvXsfnd17eDdvfbYhZkxob6V9eHOGyhyn/qxcufLXksYVpmdqoMzsO8CngC5JH4ppo4ElwERgC/Bnkt6KSxD/K2FVx93AZyU911/+EydOZMWKFfX7A46TEHNufJz93Xv6pB/V3sp/Lfp4BoocJz3M7JVi6Vl38f0HcFZB2iLgMUmTCUsKL4rpZwOT4+ty4JspaXScuvNaEePUX7rjHA5kaqAk/QL4TUHy+cDtcft2YH5e+vcUeApoN7Oj01HqOPXlmPbWitId53Ag6xZUMY6StB0gvo+P6ROAV/OO2xbTemFml5vZCjNbsWPHjrqLdZwkWHjmFFqbm3qltTY3sfDMKRkpcpzsaUQDVQorkqY+CdItkmZKmjluXJ8xN8dpSOZPn8DiC09hQnsrBkxob2Xxhacwf3qfOpjjHDY0ohffG2Z2tKTtsQuvK6ZvA96Xd9yxwGupqxukLO3o5CsPv8hr3Xs4pr2VhWdO8YdjysyfPsHL3HHy6NdAmdk3KNJKySHp7xJXBA8AlwE3xvf789L/xsx+AJwG/DbXFegEqjUySzs6ufq+tezZ1wNAZ/cerr5vLYA/MB3HyYyBuvhWACuBI4AZwMb4mgb01PrjZnYXsByYYmbbzOxzBMN0hpltBM6InwEeBF4CNgHfBv661t8fTHz628tZsGQVnd17EAeNzNKOzgG/+5WHXzxgnHLs2dfDVx5+sU5qHcdxBqbfFpSk2wHM7LPAxyTti5//HXik1h+XdEmJXfOKHCvgilp/czBy7dK1/NfmQmfIYGQWLFk1YCsoKxfnpR2dLFiyqk/6pbOP48vzT6nrbzuO0/iUOwZ1DDCCgy7hR8Y0J2OWdnRyx1Nb+z1m6pceYs31hdPNDnJMeyudRYxRPV2cSxkngDue2trrP03wMTHHOSwp10DdCHSY2RPx8/8ErquLIqcsPv3t5UVbTcXY+U7/vbELz5zSawwK6u/ifM2P1pZ9bGf3Hr6wZBULlqyqi7G6dula7nr6VXokmsy45LT3eQvOcRqAsgyUpO+a2c8IzgkAiyS9Xj9ZTiH5DhAlvVaqJPewT9OLb9feyoYwc/+5s3sPC5as4vofr+dLf/TBijQWcyJZ/ODzvPH23gPH9EgHWm9pG6nCSsec/zGaO//iI3X9zWLGGUjVYBeel/09Pb3OCXgr+nDFpNKPOzM7SdIGM5tRbP9AsfCyZubMmcoiFl/+TZ+jlhusktZSKbbceG5N34dkXdEnLvppzXoM+HSZ41XXLl07YFdoYd4vJ1BmpbQUGoCnX3qTjV27+hxbTyNVSZnUa1xwaUcnC+9ezb53y6t23XzRtD7XnE+ROPQxs5WSZhamD9SCuooQ9+5fiuwT4FEsIwMZkVzNPzfuMtANn7vpio0NZUHSruhGP/MXykTAnfEB+8SGHXR276HJjB6JCe2tTBzTyq82/6aq38l9p5yHX7FjoHeL9GMnjWPJM1vZ927v38lvsRWj1DWVf33k/+dyHs7BKKzqo6U/vv/0Vr48/5SS3aHF8mwZOoS9+9/tU2619AYsvLu3049PkRjc9NuCOtRJqwWVRAunnhSrdQ5E/kPkPa3N7Pz9PopVcpvM+Jc/+4OK86+0RZMFN180rejY3OILQ8WivwpE0xCjp8xWwUDkWr/5Rqk/Az/EQOKAYXxiw44DhnLimNbMrtW2YU0Vd+0W4+aLpnHdA+vp3rNvwOPcSB0alGpBlWWgzOxPgYckvW1m1xLmRP1fSR3JSx1Qy1mEZTeagFsl3Vjq2HobqGJdeY1G4U3aX4ugllbb8OYhtDQ30b17X9ndLIeCkWoEjhoxrM+YjFMeo4Y3VzxW6aRPrQZqjaSpZnY6sBj4KvCPkk4b4KuJYmZNwH8TJvBuA54FLpH0fLHj62mgGv3hWqwLcWlHJ1f9cFXRllDS5Foa86dPKMtLrtFboc6hzeTxbSy76qNZy3BKUKuB6pA03cwWA2slfT+XVg+x/ej4CHCdpDPj56sBJC0udvyIESN06qmn1kXLltO+CFYsfm0GxHPY0v0yR794b9FDfjfmJH594qfS1SwROqKs9+9KHPn6c4x95fFehzdUmTqDC6nf+8PJlieffLIqJ4kcnWb2LeATwE1m1kI2kdCLLbnRqxVnZpcTHDtoaWlJT1lWSNg7Ozl+1S0lD3ll2uWoZWT6D38zigahN+N3753BntEnMmrrLzjyzQ0ANO3dSU/Le9LV6BwemPFO+wlZq3AqpNwW1HDCyrdrJW2MUcZPkVRzuKNKiGNhZ0r6fPz858AsSX9b7Ph6dvEl4SZdK0OHGF/90+IOCo3eBVmMyePb2PbW7/vEBXScpEhiuoWTPNW6mQMgabeZdQGnE4LF7o/vaXPYL7lRzryYQ3U8Z2PXLiaPb2P33nf7OHEcigbXcZzaKMtAmdmXgJnAFOC7QDNwBzCnftKK8iww2cxOADqBi4H/lbIGAIYa7B+g8TmypYkRrcMSm8tUjtvsaTcsO6Q9vjZ27Spay/3y/FOYefzokvH7KuWIJuP3PY3rfZkVE0rEZXScLCh3DOoCYDrwHICk18xsRN1UlUDSfjP7G+Bhgpv5dyStT1sHwKbF53LCop+WnIsysqWpV4DWWlo1hXmV4oyv/fyQNk4DkVvQr9bW1FEjhvH0NWcc+Nxf4Np6Mnl8G6dNGtMQUxXyKz9Tv/TQgPEb68mo4c2cO/XoxFvMl84+LtH8nPpT7hjUM5JmmdlzkmaYWRuwXNLU+kusnjQm6i7t6Ow1abDceReFM+CL0TasiRsuKL3sd+GcpnrVfAvnONV7smc54wSVRkMYKKZcEpE7irVwCysmLUOHcNMfTy16Tpd2dHLVklUM9JdyLtNJdOXmTwfIJ0kjVe48rsJ7J8mKgy/h0tjU6mb+RWAyYf7RYuB/A3dJ+nrSQpMkq1h85VJLDLFyDFwSlBOSKenWRyMMZJc6N0nHWUyagSpMtV5z1VTG0sQj0x+a1GSgYgZnAJ8k+A0/LGlZshKTp9ENVCUUPlh2793PW7v7D/VSC5U+dIPBXMOeSgK8laARDJTjOOlRkxcfQDRIy2JmTWb2aUl3JqhxUJBUDa4wFt6uvfvZFwf16zmIXW307Nz4EGQ3puM4zuCiXwNlZiMJy6xPAB4gGKgrgIXAKsANVB6Fg/f5kaoLl1PoL/RKYffdQEExyyWN9YWgt7EqZCAHh+Yspn87jtOQDLQe1P3AW8ByYB4wChgGXCmp4avIaTlJVLt0QCkjNefGx6tqJeUGvO9esTX1he8qpZQHZD3XYXIcpzGptotvkqRTYga3Ar8GjpP0dh00HnIs7ehk4T2rD3S9VUquRVXYnVdui6m9tZm2lqF9BrwbadC6FKW8Do9pb81AjeM4jchABurAk1JSj5m9nIRxiiGLrgM+QAhVtCJv39XA54Ae4O8kPRzTy15mIy2u//H6qo1Tjmq781qbm7juvMbyoKqEhWdOKbrWUm6xP8dxnIEM1B+Y2c64bUBr/GyAJI2s8nfXARcC38pPNLOTCdEhPggcAzxqZu+Pu/+NvGU2zOyBUstspEUSXnRfefjFil3Fm8yKzl05lMhp96W6HccpRb8GSlJTPX5U0gsA1je69vnADyS9A7xsZpuAWXHfJkkvxe/9IB6bqYEaiJwXX6GDRI7J49vYVCS9P0pNrDwUOVS6Ix3HyYZG85kqtpzGhH7S+2Bml5vZCjNbsWPHjroJhTAGVIoJ7a1sXnwOX55/Csuu+iiTx7f12p9zkChnzKXJDIt5Dhbj5DiOMxBlz4OqFDN7FHhvkV3XSLq/1NeKpInihrTo4I+kW4BbooYdZvbKAFLHEpw/KmZI68jRQ0eOm9inKSjp1Z07ttjVO0vGoXkFsL8/kMfxmBWvLEjv7t+545V39+z8zRbggqurUVoRVZfHINMArqMQ19FYGmDw6Di+WGLdDJSkT1Txtf6W06h4mQ1J4wY6xsxWFHNvTBvX0VgaXIfraHQNh4OORuviewC42Mxa4pIak4FnyFtmw8yGERwpHshQp+M4jlNn6taC6g8zuwD4BjAO+KmZrZJ0pqT1ZvZDgvPDfuAKST3xOw2xzIbjOI6TDpkYKEk/An5UYt8NwA1F0h8EHqyDnFvqkGc1uI6DNIIGcB2FuI6DNIIGGOQ6yo5m7jiO4zhp0mhjUI7jOI4DuIFyHMdxGpRBZ6DM7H1m9oSZvWBm683sypg+2syWmdnG+D4qppuZfd3MNpnZGjObkZfXTWa2Lr4uqrOOk8xsuZm9E1cwzs/rLDN7MWpclKGO75hZl5mty0JDqXwy0HGEmT1jZqtjPtdnoSMvvyYz6zCzn2Slw8y2mNlaM1tlZhUtIZCwjnYzu8fMNsT8yg7jn+D1MSWWQ+6108wWZFAWX4h5rDOzu8zsiLTLIu67MmpYX245HEDSoHoBRwMz4vYI4L+Bk4F/BhbF9EXATXH7HOBnhEnCs4GnY/q5hPWvhgJtwApgZB11jAc+THAQ+WJePk3AZmASYamT1cDJaeuI++YCM4B1dT4npcqiaD4Z6DDgyLjdDDwNzM7inMT9VwHfB36SxXmJ+7YAY1O6Z/vTcTvw+bg9DGjP6rzk3b+vA8enfI1OAF4GWuPnHwKfzeBe+RAh9upwwrP0UWByuToGXQtK0nZJz8Xtt4EXCCfrfMLFS3yfH7fPB76nwFNAu5kdTTgZT0raL2kXwTCcVS8dkrokPUteBPnILGIcQkl7gVwcwrR1IOkXQMnoGPXW0E8+aeuQpN/Fj83xVba3UZLnxMyOJVSmbi339+uhoxaS0mFhgdW5wG3xuL2SutPWUcA8YLOkgSLa1EPDUEKA76EEAzFgcIM66PgA8JSk3ZL2A08CF5SrY9AZqHzMbCIwnVDDPUrSdgiFT7D4UDrO32rgbDMbbmZjgY/RO5pF0jpKUXYcwjrrSISkNBTkk7qO2K22CugClknKRAdwM/APwLvV/H6COgQ8YmYrzezyjHRMAnYA341dnreaWdsA36mHjnwuBu5KW4OkTuCrwFZgO/BbSY+krYPQepprZmPMbDihx6rs5+igNVBmdiRwL7BA0s7+Di2SpngyHwR+RbjAlhMmD9dLR0X6MtBRM0lpqDWfJHRI6pE0jRB2a5aZfShtHWb2KaBL0spKv5ukjsgcSTOAs4ErzGxuBjqGErqgvylpOrCL0A2Vto5cPsOA84C709YQx4bOB04gLF3UZmaXpq1DYeWKmwjDJQ8RKv5lP0cHpYEys2ZCod4p6b6Y/EbsuiO+d8X0kvH/JN0gaZqkMwiGYmMddZSiv/iEaeqoiaQ0lMgndR05YhfSz6mg+zdBHXOA88xsC6Hr9+NmdkcGOpCUu2e6CJPwZ/X/jbro2AZsy2vN3kMwWGnryHE28JykNzLQ8AngZUk7JO0D7gP+MAMdSLpN0gxJcwlDBGU/RwedgTIzI/RBvyDpa3m7HgAui9uXAffnpX/GArMJTeHtsQtnTMxzKjAVKLuJXIWOUtQUhzBBHVWTlIZ+8klbxzgza4/brYSHwYa0dUi6WtKxkiYSrovHJZVdS06wPNrMbERuG/gkoWsnVR2SXgdeNbPcsszzqGDNuDrcK5dQYfdeghq2ArPjEIURyuKFDHRgZuOJ4SWMAAAA7klEQVTj+3GEhWrLLxNV4XXTyC/gdEIX2BpgVXydA4wBHiNY78eA0fF4I6zWuxlYC8yM6UcQLu7ngaeAaXXW8V5CDXAn0B23R8Z95xC8aDYTlivJSsddhP7sfTH9c2lqKJVP2mVBqKx0xHzWAf+U1TnJy/OjVO7Fl1R5TCJ03awG1md8jU4jeNyuAZYCozLSMRx4E3hPhmVxPaHitA74T6AlIx2/JDxHVwPzKikPD3XkOI7jNCSDrovPcRzHGRy4gXIcx3EaEjdQjuM4TkPiBspxHMdpSNxAOY7jOA2JGyjHcRynIXED5TiO4zQk/x9PFB172t4ruAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for f, filepath in enumerate(filenames):\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(f'f: {f}')\n",
    "    df = read_GPS_SONEL(filepath)\n",
    "    df = decompose(df, 'Vertical', filepath[-24:-20], 365)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_trend_toyearly(df, res):\n",
    "    \n",
    "    period = df.index.year.value_counts().max()\n",
    "    yearlytrend = res.params.x1 * period\n",
    "    \n",
    "    return yearlytrend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x24f96a176c8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3hURduH79mWXkghgUAIvQZCRwQFUYoKiICABSmKioivHSyIgIqioqiv7UNFgRc7oqB0pKM06U0JEEhICKT33fP9cZJNNrubbJLd1Lmva6/dMzPnnDnZ7G9mnnnmGaEoChKJRCKpW2iqugISiUQiqXyk+EskEkkdRIq/RCKR1EGk+EskEkkdRIq/RCKR1EGk+EskEkkdRFfVFXCUoKAgJSIioqqrIZFIJDWGffv2XVEUJdhWXo0R/4iICPbu3VvV1ZBIJJIagxDinL08afaRSCSSOogUf4lEIqmDSPGXSCSSOogUf4lEIqmDSPGXSCSSOogUf4nEARIzE8k2Zld1NSQSpyHFXyJxgH7f9mPaxmlVXQ2JxGlI8ZdIHGR37O6qroJE4jSk+EskEkkdRIq/RCKR1EGk+EskEkkdRIq/RCKR1EGk+EskEkkdRIq/RCKR1EGk+EskEkkdRIq/RCKR1EGk+EskEkkdRIq/RCKR1EGk+EskEkkdRIq/RCKR1EGcIv5CiM+FEPFCiCNF0gKEEOuFEKfz3+vlpwshxCIhxBkhxCEhRBdn1EEikUgkjuOsnv+XwOBiaTOAjYqitAQ25h8DDAFa5r+mAB85qQ4SiUQicRCniL+iKFuBq8WShwNL8j8vAe4okv6VorIb8BdCNHBGPSQSiUTiGK60+YcoihILkP9ePz89DLhQpFxMfpoVQogpQoi9Qoi9CQkJLqyqRCKR1C2qYsJX2EhTbBVUFOVTRVG6KYrSLTg42MXVkkgkkrqDK8X/coE5J/89Pj89BmhcpFwj4JIL6yGRSCSSYrhS/FcB9+d/vh/4uUj6+Hyvn15AcoF5SCKRSCSVg84ZFxFC/A/oBwQJIWKAl4H5wLdCiMnAeWB0fvE1wK3AGSADmOiMOkgkEonEcZwi/oqijLOTNcBGWQV41Bn3lUgkEkn5kCt8JRKJpA4ixV8ikUjqIFL8JRKJpA4ixV8ikUjqIFL8JRKJpA4ixV8ikUiqIVm5Rnq/vpHNJ+NLL1wOnOLqKZFIJJKKk2s0sfvfRNYejWPp7vMATPziL6Ln3+b0e0nxl0gkkiokIyePracSWHv0MhuPXyYlKw8Pvdbl95XiL5FIJJXMtfQcNhy/zNqjl9l2OoHsPBP+nnoGtg9lYLsQbmgVjF6rofnza+jZNABFURDCVkzM8lP7xX/JUEg4BTc8DaEdIaQ9uHlXda0kEkkd41JSJuuOxrH26GX+jL6K0aTQ0M+dcT3CGdg+hB4RAei0ltOwrjD3FFC7xd+YC2e3qp/XPJ2fKCCwOYRGqo1Bq0FqgyCRSCRORFEUzsSnsTZf8A9fTAagZX1vHrmxOQPbhxAZ5uf0Hr2j1G7xF1p4ZCfoPUGrh7jDEHsI4g7BxX1w9CfY+Ao8HwsGz6qurUQiqeGYTAp/xySx9uhl1h2N498r6QBENfbnucFtGNQ+hGbB1cPyULvFX6Ox7NX7NYLWQwqP/9sb4o/Ca0V2kZydXPp1M5PAzQc0rp+UkUgk1ZuiHjrrj13mcko2Oo3guuaBTLw+glvahRLq517V1bSidot/aQx5A5bcbpk22w+GvQ+rHoPJGyDjCkT0LZwnyM2CN5pAw84wZUtl11gikVQD7Hno3NgqmEEdQripdQh+nvqqrmaJ1G3xb9oXBr4K616wTF/1mPq++ObCtJAO0OJmaH6TenzpAJhM6uhCIpHUeux56NzSLpRB7UPo2zIYD0PNsQbUbfEH6D0N9n8F/o3hzAb75S4fUV873i1Mm1PPMTORRCKpkdjy0GlQiodOTUGKP8C0P0FRYO9i1S30wNdw11ewbFTp5/75GbQfAV5BhWmKAlU0gy+RSMqPPQ+dFvW9efjGZgxqH1qlHjrORKgba1V/unXrpuzdu7dyb5qTDsvHQPS20sve9KLqTXTsZwjrCg9ucn39JJVG5JJIAA7ff7iKayJxNiV56AxqH8rA9iE0ryYeOmVFCLFPUZRutvJkz78kDF4w4Vf1s8kEOalwYjUoJvi52E6Um+YVfr64D3Z+AAHNoM2tlVdfiUuoKR0kieMUeOisO3qZdcfizB46vZpVbw8dZyLF31E0GnD3g6i71WN3f9Do4McHITvFunzRSWSvYBjwMrS5DTwDKqe+EonEgpI8dAa2D2FAm+rvoeNMpPiXl7b5LqIzL8DVf2FRZ2g5SHUJPfKDZdn0BFg1TX0BvJwk5wQkkkqgwENn3bHLbD1V8z10nIkUf2cQ0MzS6+eWOfBhL9VMZItX/NX3Z8/KkUANQEGafWoS9jx0xnZvzKD2ofRoWnM9dJyJFH9X4NdIHRGsfwnaDIUrpwp7/UV5s6m6bmDwG2q8IbliWCIpMwUeOuuOXWbt0TgOxdReDx1nIsXfVQgBA/MngcN7Qpf71M+z/SzL/bMJPuxumRZ5F3QcAy1vRiKRWGPPQ6dTY3+eHdyaQe1Da6yHTmXhcvEXQkQDqYARyFMUpZsQIgD4BogAooG7FEW55uq6VAse2QWf3gjGHPtlDn+rvloNgbtXVF7dJJJqTK7RxJ5/r7L2aFyd9dBxJpXV8++vKMqVIsczgI2KoswXQszIP36ukupStYS0g5cS1M9/LIDN8yColWoaKs6p32DHIrh+euXWUWKBdPWsOmx56LjrNWoMnfahdc5Dx5lUldlnONAv//MSYAt1RfyLcuMz6isrWV0n0HUibF+o9voLWP8SdBoH3sFVV0+JxElk5hhRUPA02Jeea+k5bDwRz9qjcdJDx4VUhvgrwDohhAJ8oijKp0CIoiixAIqixAoh6ldCPaov7n5w6wL188jPQGuAg0sL899qUfh5/CpodmPl1k8isUNsciaBXm4YdLa9Z9Ky88jMMRLkbWD9sctM+XofYL1DlfTQqXwqQ/yvVxTlUr7ArxdCnHD0RCHEFGAKQHh4uKvqV/2440O46QV4p6113lfDoNPd8Pdy9VgGlpNUAX+evUrbBj5c97oaxuT6FoEse6CXOT8+JQtfDz39FmzhSlq21fnxKVmkZOWy9qj00KkqXC7+iqJcyn+PF0L8BPQALgshGuT3+hsA8XbO/RT4FNTYPq6ua7XCt6Eq7CYjfNwH4o8V5hUIP0DCSQhuXfn1q0NIP/9C0rLzOBSTxN2f7bFI33EmkW//usDwzg1p/eLvpV6nx2sbzZ+lh07V4NLAbkIIL0CjKEpq/uf1wBxgAJBYZMI3QFGUZ0u6VpUEdqtunNkIS++0Tpe9f5eSZ8qj89edgbob2K1AJ5rOXOOU613fIlANmiY9dFxKVQZ2CwF+yh+66YDliqL8LoT4C/hWCDEZOA+MdnE9agctBqhCX3ytgETiYprOXMOANhWfmlsxpRcmk0LvFkGlF5a4FJeKv6Io/wKdbKQnovb+JeWhoKdf0AjM9oMnjqoriyVOp7abfZ757m9ualOfIZENbObHp2QBsPGETeusXXzdddzWsSGP3dSCQzFJXN8iCB936ZZZXZArfGsLC/M3qo+8C258VnUfbWRztCeRmFEUhe/2xfDdvhgLD5yYaxl88se/TO3f3Dypa4/f/9OXep4G3PVa9FpBdq6JpbvPcf/1Efjmi31Dfw+XPoek7Ejxr8lMWANfFtsvoGB1MECrweATCkPfq/y6Sao9qVm55OSZrNIVRaHPG5sB+Hr3uRKv8fCNzWkT6muR5mmAxwa0dF5FJS5Bin9NJuJ61QQUsw/+7ybr/FP5XheXDsLk9aDVy1DSEjORs9dZHM//7QTHY1P441RCqeeeff1Wvt8Xw/CoMFdVT+JipPjXBhp1hSlbIDUO/jfWOj/2IMwLhiFvQse7IDtV3YjGt2Fl17RmUktM/muPxvHQ1/uY0DuC9ccuW+V//Mc/pV5jVNdGLBjVESEEo7s1dkU1JZWEFP/aQkPVFZHZyXBxP+x8H47+aFnmt2fVVwGtb4MxX8tQ0rWAuOQsgn3c0GoEGTl5pGTmkZadS1auCV93PV/vjuazbWcB+HJndJmuPa5HY+YM74BerrCtVUjxr42EdYHRX6ivzCR4o4ntcidXw5wAmLobDixVN6GRDUGNIy07j16vb+TeXuH0bRnMQ/khFCpCkLcbe1+8mZw8k93QDZKajfxWazse/jDxN8u0Qa9ZHv+3F+z6oHCOQFJteXX1MdrP+p1f/r7EgfPXSM3KZcWf5wFYuvt8mYQ/LN8Dp3O4P2+NLvTI/vWxPuyaqc4hSeGvvcief12gSW+45wdIvgAo0G0SXPeo9WKxFXdD/XbQbya0G1YlVa2OVBc//+d/OszyParQP/a/AxW6VoFbZ3JmLm46Ddn5Xj+LxnWmQ5hcRFgXkOJfV7C1K9gLl+HPT0AxwYbZalr8Mfg2f9exaXshSLrsVRcKhN9Rmgd70Tm8HqO7NqJDmB9ebjr2nbtKerbRXMbPQ/XDd9drrSJtSmo3UvzrMnp3uP5x9bNfY/hhsmX+B/mLxG5fqI4WJC5nz7+JfLD5DJ9P6G6eYO23YDPRiRkOX2PeHR0Y3a0Rbjrr+ZuuTQKcVldJzcalgd2ciQzsVkkc/B+sfNh+/vOXwOBlmWYy1uqJ4hxjDl2XdgWcH9gtK9fIkp3RnL+awZjujRn2wQ5znrteQ1au9SKs4rw9uhMt6ntzLDaFlvW96RYhBV6iUpWB3SQ1jahx6uYy53fBzkXW+a81hPtWqq6lBm849A38PLVWxxZypc2/3azfMeVfflkxs4494fc0aPnjmf5k5ORR38fdvKtVp8b+LqunpPYhxV9iTZtb1deAWfDlbXDBMnY7X99hfU7CyUoT/yNXjhDkEUSoV2il3K+smEwKl1OzaODnYZUen5qNQafh1OVUBJiFvzQMWg0rHupFRKAXOq3Ij5nj5vS6S+oOUvwl9tHqYXJ+CIAjP8D3Jdj9l94JvabC4NcL01Ivw+Jb4J7vIbiV06o1bvU4NELD3+P/dto1ncnb60/y4WZ1tezZ12/lh/0Xefo7x+vqZdCSnlM4Kbt75gAZ817idKT4Sxyjw0iI6At/vAlX/4V/NlqX2f1f9QXQ4mZofSsknYNd78Ow951aHZNSui3cWZQ0L2YyKWg0hfGSUrNyzcIPMODtP/j3Snqp93hvbBSPrzgIwNE5g8nOM5JnVPBykz9RiWuQ/1kSx/GuD7e9ZZlmb2OZMxvUF8D+r6DpjRDcBkI7uLaOLqa42Ld4YQ2T+zSlV7NAJi+xdkhwRPhfvzOS4VFhhPq6k5SZC4CbTovUfYkrkf9ekoox66oaIqI0CtxIJ2+A+m0hJ00NN13DGPHRTn56pDdzVx/jfGIGJgU+23bWHDfHEZoHe7HuiRvRaiwjrPZsFujs6kokdpHiL6kYGi28EAe/PQc9HlTNQjn5vV1bpqHFRRabjVsB4depISiqMaYis7J/X0ii26sbuJqeU+I543qEs/VUAheTMgnxdWNw+1BeGV6zRz2S2oUUf0nF0XvAsHy30DFfW+adWgfL7WzRXBB++rlz1a4B+ONUAvd//icAi+5uZ5FXmvB/cHdnbu/omnDZqTmpKCj4Ggo3UFEUhS+OfsGgiEGEeZctvr7RZOSh9Q8xOXIy1zW8ztnVlVRjZNQmiWtpNRAe2gqzrqmriG3xRhNY/3Ll1qsEzl5JNws/wPQicXSELsX82VAsxHGAl4Ftz/Z3mfAD9P5fb67/3/UWafEZ8Szct5CpG6aW+XrJOcnsidvDlPVTnFVFSQ1Bir/E9TToBBoNPHEE+jxhu8yOd+GfzZBwyuHLlmV1uqIoVuVPXU61MOn8k5BGxIzV9H9ri+XJorCMR9hSQJ2kPfXqEKb2aw7Azhk3sXPGTTQO8HS4Ts6iwPMpPbf0yeXiaEXZVmZvjdnKzos7y3wfSfVDmn0klcvNs9UXwG8zYM9HhXkFi8dmJzt0qTxTHnqtvtRyiqLQdOYaBrcP5fGbW9K2gS+HYpIY9sEOmgR6cq7UuDmFbqX+PtkM7hLG2O7qKObpga2Z3Kcpgd6uXXCVkJGAQWtw+nXLGt7l0Y2PAs4PcyGpfKT4S6qOIfPVcBKf3GCZPtsPXowHXcmCmmXMKlH8Vx64yLnEDIz5Avf70Th+PxrHbZENWH04FqBE4d/34s3c/M4fJGUXCqSvh453RkaZjzUa4XLhB7jpOxt7NOdjVIx280qjIudKajZ1wuyjmCpvQZCkbGR4NiHp4RPQ/UHLjHn11UYgJ13doN5kgsPfq0Hk8rmQZHuEkJGTx3sbTvOfbw6ycMMpFm08bc5zM2bx26GYUuvVOMCDQG83DswayN4XB5TpmZKykvjyyJdl7lWXl1xTbrnP3RO7p/RCTuTvhL/59uS3lXpPiW1qdc/fmJfH4scfpEO/W+g9+u5Sy8ccO4JvcH18g+vbLZNyJYHY0ydpfV0fZ1aVrLQ0zh0+SHJ8HIGNGtO8a89Sz4n75zR6d3cCw0rfSLtAiLLT0zm08Xe6D70ToXF9238t7hJ+9UPIzcpG7+aGRmtpY/7quemkX7vKU9/8SsK2ZZgQhLgXsV2/pk6e5hi1ZBp1nAzdQYNcd2KDsrjt/U2cnn0vQgiy8nJYuH8BD0Q+QK+56m5WfrnJJOvVRWh6Uy4NsmIZfnk1R73bsCm4Pw2zLtE4M4a//Ltiyrd9H31lkMWq2pQr8WRrsrnv93C0JsGW2y09fdKuJrL315/wub49SboMBjQZwNxdc1h3fj2RwZF0Dela4t8n7p/ThDRtbve7uHz2H4IaW2/DqSgKQqjrBPJMefik6xCewqpccXKzslAUEwYPT5LjL/P2z7Og2DKN1MQruPv4oNXqWHN8FQEBIfRu2Nvm9f498BcNWrTGw8fXZn5xJqxS94q4q/VdDpWXuI4qC+kshBgMvAdogf9TFGV+SeXLE9LZlJ3N4ofHo+h0NOrYWf2xCIEQGoRGIPI/IyA9KYl/9u4GoNPA20BRuHT6BF7+9UiJv0xoi1bk5eRwavd2ADRaLRqNlla9rkfvYT3Jd+3SBRAa6oU2AGH5w85IvkZudjaXTh4nomNnTCYjZ/7abVGmcfuO5GRm4OVfT/1h5f/QY0+dICc7i9BmLcznNGrXgeyMDIIbN0Gj02PKyyX2zEkMHl6ENm+Jopg4tOF3hEZjHgV5+vmTk5VJXnY2rXr1QavXo5hMnNjxh7kOencPmnftYT7Oyczg3/1/5T+/Dr/6IWRnpNO8aw/OHf4bv+D6XLlwjszUFEJbtCIrNZWky7FWfxshNCglhGfQCSNBbumk5LqTYbRv5z7p05hsfBEohGefxS8nk3QDeBXzxLxiqEdQzjW71yl41tysTEKbtyQoPIIjm9fbLevlX48mkVFodLoSyxUQFB6BT0Ag0X8fwK9+COEdOpGefI1/9pat120UClqldIFv2Lodl04eA6BBy9Ykx1/G3duHqxcvlHpu8269zL+D4hj8fclJSuGKXzZByZamriYdO+MTGISiKBzdoq7s9vD1o0X3XhzeuBaA0FatiTt1EoCOAwYjtLU3DLgzMXh4cMPdE8p1bkkhnatE/IUQWuAUcAsQA/wFjFMU5Zi9c8oj/jkZOXz+wMtk555GrzOpgXmVwgC9BZ8UIFfJM5+nz+8F5ipGBAIFBQ0Ck53QvnobHhO5+bZUndBR/Oda9F6O4IYeBbWnl63YH+Jr0aBFg0CQTa7VswDqwwrQKRryhKUA23pGD+Fm/iuZFIUcrO+vQ0seVWM71psU8tCiaJxj2nMT6hxCSX/nAgzoyMGx71IntOTlfwd6oS38PmoJ6n+d2rDnFXk2LRqMBRPm+f97YPs3I7GNQWvg4WXflevc6hjPvwdwRlGUfwGEECuA4YBd8S8PBk8D9cJu4ErqLRaSJoq9gzr8KE5F/j1LOrc81y2oq6OxHYuXs3VPR778on83TQn3r0r7oStkpLx/55Io+jeqzdJX/H+hdH8sSUnolEzXXNclVy2dMKDoGDQGsDJyCyGmAFMAwsPDy3WjoS8NIPFSmkNlM1NzMLjp0Bo0VukePqr5QTEpZCTnoNVrMHjo0GhtD8Nzs4x4+hrQaAWKkm9zV0BRUEcf6jAkPw8Uo2Kefs9Ky8Uv2CP/foX2esWkWFwrKz0XoRG4eepUE5ZGIDSoNvC0XLVuovCaBncdGp3aQ3Pz0pOZmkNOptF8L1eQlZGLVitIT1JtMd4B7mg0gpxcI5PWTkCr6Hg6Yj6N63uTl2fiyz3nGNKtEfN/PkqWUBsfLaBTQC9ymGP4jFeDAtCZ9KS5JfFJ7GUW+/vyp4cH7rleZOiT+Dw2kfeMQ3lk3CSMeQrpuhRe2jILP3c/3h26gAKLU0ZKDlcvpRMc7o1WryXpahqZ2jSCPYPJyzGRk5mHojHxzM6nMCkKHrneKMLEomFvIwSMXzseBQXPHD8yDMnoTHomd5rM35cPsj/2IBPDHyIypAPeAW7m78DNU4coEtPnSvJV/nfhK64Pux4PvQdv71xIniaHntp+HMz+C89cH7J1qkdSuiHZ5sYyGkWDW54nGYZU1PGswn9v+S+BHoG46QzmuQGhAZNRITk+k8y0HAweOh7aMRmjyMNgdOeOFiMY3moYk1ZPpkFqc+J8zuKd7U+6IQnPXF9S3a6Sq80GBEsHLQUFMlNz0btpEcV+B5mpObh56lGEiUfXPYZR5GLU5OGe54lXjj9JHvGF9UdDgHs9gjyCCPIIItAjiCDPIILdgwn0CMTPzc/8DHURrdY1c3NVJf62vkmr/2pFUT4FPgXV7FOeG3n6GvD0ldvaVSZpOWlk5GVQ39P2xPmhmCQmL99HVqN/AXhg5zG+nHgd0Vey+C7hGt/9dg0MoPPbi0fD7wFIPT4f0DFGPxKvoHfM16qXdIGfQi0nvIf4wbOJH9Ds59fh3h+4ENiU6IAjAFz1iqVVvcK9BZp0KAym9vq5l9h8YTMH7zuINn9bylk7ZhHtf8Ti+mGt6gFweWe01bPNjZ5Jy3otifE/jXenXNo2a2D373Q08Sjj96ghLn479RPP93yeS35nAPiFr6AC68XG7hxBv8b9eP8m61DaA9b3ZVybcTzT+RkuHyp8hh/Sv2JEg4Fc9D/FRf+SF9vVi3DHQ1dyp2HlmZUcvXKU8/WOWuW90vsVwrzDuJh2kZjUGC6mXeRi2nn+TtvFlStXLMq6ad1o4NWAMJ8wGnk3Isw7TH3lH/safOt041Beqkr8Y4Civ9hGwKUqqoukgvyT9A8eOg8aejckLj2OW76/BYCVt+7i1q+foWlYMh6Jj/LNlOu4lpFj3qfWW9EghAmEwoQv/rK6boHwA2jcYjn+0n10W2Zpvnwy8iZIOV38VN4MrMd9KamwdCS5eh00Ur2GRq4aaXOB0o+nf2Tzhc0AbLqwiZvDb0YIwU9nfrIqe93y69g2dpvdv8f5FHU7xrxS5nbe+PMNi+PX9rxWYvnSaOHfgjNJZ8zHf1wonLzfd3kfE36fwIrbV5BryuWrY1/hrfe2OP9yxmXzd1caf8b+yY2NbyyxzEs7XrKbF+YdRs8Gtj3asvKyuJR2iZi0/EYh9WJ+43CRQwmHSM1JtSjvrfe2aBDCvAsbiYbeDfHUV/6q65pAVYn/X0BLIURT4CIwFijdF1NS7biUlMkdP6src5v5NePf5H/NebeteABD4FEuZkF6wt+0ejHJ8mQl3ywlTFbjPqGz9M7xavYe3Za9Z3X/PTaEvzgPhlqOQNLT4/Hyskx7eWdhbKEntzzJrOtmMarlKJvXS8tNY/05+14+2cZsAJKzS16pfCD+QIn5ZaWRdyML8S9qItpyYQsAOy4WbhD/37//W+57Tds0jcP3HyY1J5Xe/+vNghsXMDhisDn/WGLJ03c6jX3pcde508y/Gc38m9nMT8lJsWgQCl7nUs6x89JOsoxZFuUD3AMKG4diDUQDrwYOrRKvjVSJ+CuKkieEmAasRTXpfq4oivXYUFLteeb7v80zekWFH0DvW/iVekV8RNrp51HyCv3BtRoNJoxQ3FNIk4F3S8tecXk41+Ueno3bQLzO8t981P9u4LeJ+b1/Oz/8j//+GHet/encZ7c+W+r93/zrTe5rd59FWkZuBqeunSIyKLLU88tKbgkeSgXxf94/4Nwd1c6nqqOcxYcXW4h/afcpSfxLw9fgi2+gL20D21rlKYpCYlai1YghJi2GI1eOsOHcBosRmUZoqO9Z39wwNPJuZG4cwrzDCPYINpsAaxtV5qShKMoaYE1V3V9SMeJTs3h9zQl2nEnEx/o3aBOP8P8j498nuS2yAR/e04Uey14hMy+XxRO6Mvnz4wD8OPU6HtsxjJSSoyY7xLi0g6S6WYdeiNHrefyTNmzy8qSLW32WjNlgVSY+I57ntz9f4TrkmfLIMeag1+jRa/U8t/U5tsRsqfB17d3LHteySl7nUFa6hXRDURTm7JoDwImrJ3hl1yu81OsljHlGehh6MLLdyHwHUGtEvOD4leNOrVNRDBhoSlOa6ptCPdQXauNgUkzkKXkYTUaMipE8Ux5Gxage5xrhKiRdTSKJJAQCrUaLVmjN7zqNzvyuEdUjSIK7uzuNGjVCr3d8FFOrV/hKXEePV21s1FIKWrd4wESnxuqq2wJh6NjIh+j5twEQucR5PeLU3FS7eZu8VDvw/ux4Ir/qWKH7bB+7nSVHl/DZ4c+s8jp/3RmAtgFt+WLwFy4TfrAfqmHLhS388u8vTr1Xak4qKTkpFuad7099z/envueFli/QqVEnDD4GuxOx7QLbVctJWpNiIteYS64plxxTDjnGHPVz/rvRVGyUKkCv1WPQGArfNXoMWvW9MkYNiqKQmJhITEwMTZs2dfg8Kf6SUknLzsPbTWd2OX3+pyOlnGGfoTecYuL1qtAXhBMuCC6Wa3QsRk1Tv6acTXZ820RX8uGAD/Fz82N6l+mk5KTwzclvbJY7fvU4vZb3cuiavRr0Ynes7VW2BbzU6yXm7p5b6rWOJx5nb1zZFkcW5+0b31Vklz4AACAASURBVOapP56ySDt57aTd0USQPshK+D31nmTkFgbRq47CD6oZyE3nhhu2g/UZTUaLxiDHlEOuUX1Pz003m9cK0Gq0Fo2BRSOh1Ttl5CCEIDAwkISEhDKdJ8VfYuZcYjo6rQaNgAZ+qhvfhasZ9H1zs9PusSXhS/RaVUg0+fFsCswVH/39kd3zinI2+Szf3P4NY34dY077854/8dB5cPraae5cdafT6lvAsObD8DX4MrzFcEb/Urgz2Q2NCiOS3t/+frviXxY+G/gZP53+iVk7Z9ktc1fru8ziP+/6eby440Xb5X6teAydgRED4Q/r9KErh9osLxBW4q7XFJojAtxrruu1VqOaf9x11vNBiqJgVIzmxqCggcg15ZKVl0WqKdUq2J9OoytsGIq96zV6YlJjyDJm0bJeyxLrVZ7GVIp/HWHXP4l0auyHp8H6K998Ip6JX1q7Wj55SyuW7TlX4nXbhPoQo2gRwvFwBe/ue5c+YX3Q5K9qK+gtRadEO3yNdoHtGBIxhN+ifwMw+5yX9iMpDwPC+vJqn1fNx2Nbj2XFyRW80vsVi3JFBa6ijGg5ghEtR5iP98btZeLaiRZlbmh0A1tjtuLn5ue0+7oKX4MvOo0OP4OfTeGsDQgh0AkdOo0OD6zXQCiKos4BFRktFLxn5GZYeYcVhJZxFVL8azEmk0KeSeFqeg7jPtvNrZGhTO3Xgu1nrtA9IoCRH+1ErxXkGm3/g72z3v5Cn+UP9qR38yAAor4S9Ajtza5Yx3Z4WnxkMYuPLEYn1H+/PCWPv+L+KtF90hZv3vgmgR6B3N3G0ku4gVcDYtOtA8qVRphJw0UbMYIWbl8GSVnQ+V4Iv44Xer3AC71esCoX6hVa5nsCrB6xmi+OfkHPBj2p72F7YVy7wHZWacUbz/Lyyc2f4Ofux6u7X+XxLo9zKe0Szfybce+ae81l9t27DwWFRfsX0SagjcOT4b5uvnjqPPEx+ODr5ljkT2cghODJJ5/k7bffBuCtt94iLS2N2bNnO3yNLVu2YDAY6N1bjWg6YcIEbr/9dkaNsu0C7Eid9Fp14t9WzAuTYjI7COSYckjMTCTH6ATPBztUj6nqKiIpK4nIJZH8ePrHSr1vReKvl4aiKHyw6TTnEtN54Ku9tHrxN3q9rk7Orjkcx+3vb2f+bycY+ZEq1PaE3xZtG6g/3qZBXmbhz78rHYLam4/shf8tToHLndFkZNLaSSWWfabbMzbTn+vxHI19LVf4Tupg/1qzrrNvSlk+bjOHxh+yShcAR36Ar0fAq6GwZBjs+hB2f6wWiDsC19QR0ndDv2PtyLUlPktxwn3Defm6lxkcMZguIV1slvHUe1otTisY6psUE7c3u93h+xXd/N1D50HvsN60D2zP8tuW07NBT0a0HEHHIHUS/M6WqgnNoDXgpnXjme7P2F25bYv6HvUJ9AisdBu/m5sbP/74o9VqYUfJy8tjy5Yt7NxZeVtWaoQGg9aAt8GbAPeAUldQV5Q63fMv6B0uP76cxj6NqedWjxb1WliUyTPlVcgnuYBrWdfIysvixNUTTN88nZd6vVRqTHOjSeGbvy5wZ5cw3PWFXgObT8bj665jzi/HeGNUR9qEqj9mk0nhtve3czw2hbfWqb12z6bvkZvUldxrpe8/sPLR69lyMp53N1gunPr1sT50CCvdtPBu/3fx0nvRq0Ev9sbtxcfgQ1x6HDc2vpHtF7fzyIZHbJ63//J+m+lTO01l+6XtvNbnNZr4NmHB3gWl1gFgTOsxvLqn0EwTFRzFA5EPYFJM9A/vT3pOOm/ve9vqvOK26I9ufBf3FeOsb3D2D/UF8PtzhenPnqVNQBug0CTjChp4qSEjCrZ1zDHmkJbrWPyqzwZ+Rq8GvTifcp7bfrqNQPdAm+WEEOy5ew9uWuuJT1supT8M+4FXd7/K/njb32Vlo9PpmDJlCgsXLuTVV1+1yDt37hyTJk0iISGB4OBgvvjiC8LDw5kwYQIBAQEcOHCAgIAAduzYgVarZenSpbz/vrpuYevWrbzzzjvExcXx5ptvlnsU4Aiu9hSq0+JfIOrXsq6Ze55Fe1f7L+/n/t/v54tBX9At1GZUVIcZ9MMgMvMKo/PN3T23VPFfezSO5386zIbjl7mlXQjjeoSz48wVJhYJhTD43W1Ez7+N47EpvLzqKMdjUyyuoXWPRRv6q03xD/V1Jy4li9dGRDKuR2OEEEQ19ufxAS35du8FhkdZNjr2KLBLDggv3PGq4O/VOqA1AH3C7Dc+8/bMs5n+SNQjPBJlu8EoCSEE28dup8+KPoxqNYqXr3vZIv+etvegoPDOvnfsXEGlT8QAeHgfRG9To+/9PLXkG7+Z72b38A4+7DOf0ynnuHONjcajAvx1z19mD5G729zN2ui1dAnpQuuA1uy/vJ/O9Tuz4MYF3PXLXTbnUHo1UD2O7PnfF8VeWITi575141u0qteKJUOWWLjqGrQG5q85w7Fi/5MVpV1DX14e2r7Uco8++igdO3bk2WctF+RNmzaN8ePHc//99/P5558zffp0Vq5cCcCpU6fYsGEDWq2W2bNn4+3tzdNPPw3A4sWLiY2NZfv27Zw4cYJhw4a5VPyDPIK4mnnVZdev0+Jf0HOKzyyMMJhjzDGnF/hN747dXSHxN5qMFsJvi53/XCHU151v9l7gkz8sV8puOhHPphPxzPzR9qbZf0VfZfTHuxyuz4m5g8nMMeLjruPPs1fp3SLIIl8IwZjuZYui6siwfumtSy3syOWlfWDpP3w/Nz92jNuBl87LKk+v1TOxw0Ti0uNYfmI5ABG+Eeb8oc2GFvrF+zeGqPw5heY3wXcT4ELJbph8fD0ALQGaWv4d7293P9nGbOp71mfRgUVW8XVKo+hkaZeQLhadlR3jCkM3/DjsR7ostTQhPdn1SfNnZ04kDooYZJVm9uMXVbdw39fXl/Hjx7No0SI8PApNKLt27eLHH1VT73333WfROIwePRptCZvM3HHHHWg0Gtq1a8fly5ddV3mc60Bgizot/rYmyubsmsO8Pvk90TKYKY0mI49vfpw/Yv5gZMuRzO49mz8u/EG4bzjv7beOSVOAoijM+vkoX+8u2aumAEPQehSTgdyrhUG1yiL8AO56rblHX1z4y4OjGwJ1Cu7Ewn4LeWLLE6WW/X7o91Zp3UO781fcXzzU8SGH7lfUtm2LmT1nMqPHDD499Cm3NrvVnP5a39d4ra+NIGu+DWDyWnVf4X82QWocrHnaoboU8HR3tXxydjKLDixymS1cr9Wz/Nbl3L2mcDJ8YoeJJZzhOBF+EXbzFvVfREBKgPm5HOmhu5L//Oc/dOnShYkT7T970e/Ay8u6s1AUtyIrxqtqF0RnUafF39bEq63FQzmm0mfcEzIT+CNGtQP/cPoHnuj6BNM2TSvxnAVrT9Ak0INvoxci9P1Qckv2f/Zu9TJCqwYN+99dMxhlQ/S3PdufxgGeRMxYbZF+5tUhHItNoWMj/1KfpTw4YkYAuLnJzaWWsRV1E+Djmz8mMy/Tqa6NQgge6uRYY2LG4AVt833ctQYI7wXLx6jpl+0vgDvcaDT8sxma9zcLh6N/t/IQGRxpjvT5yS2fWD5C/ug22DO4zNdt6N2QvffupdvSbjTxtdxfuH94f44fd13YhrISEBDAXXfdxeLFi5k0STXt9u7dmxUrVnDfffexbNky+vSxbZL08fEhJcW5Jquy0tinsVPmHG1Rp8X/+1PWvctDVw6x4dwGbm5ys/mH+cWRLyyGzLZY8JflZGR8RrydkoUsPjWb3Kt98Iz4E41bPFmXRuMXupNlI16nfcN6/JOQRj1PAwFeBhRFoeNXM8zndmjkyeMDWvLexsLJ2WUP9KRxgGqnXTimE/4eBqbnr/jXaTUuE35nmhB+HfGr3TyD1mAWrWpD1/vV98cPqu+zS2iYtr2tviavxyesK63rtebRqEddWr1vb/+W9Nx0/N0tv/tQr1DmXT+vxLmYknDTurH01qWE+5Rvk6XK5KmnnuKDDz4wHy9atIhJkyaxYMEC84SvLYYOHcqoUaP4+eefzRO+lY0r3WNrvfj/cOoHTiedZlrUNLwNlvbVtgG2I5I9seUJ9tzt+ObacelxrDu3ziLNkVWmet+jmLJUN0V3vZGg1j+QZPwX3C4B9WgeXFjflBzLHsjqf1fzxC0j6RZRj5/2X+T6FkFcX8SEM6JzI/VD2fYILzfOMF8cGn+o2i77d5jpByDhJCT+A10n8N7WV3j80u+WZRbfgrb7A3zf+TloVDFHgtLQa/X4a203+sNbDK/QtTsFd6rQ+a4kLa3Q+ykkJISMjMLQEhEREWzatMnqnC+//NLiuFWrVhw6VOj627dvX7v3qInUej//H0//yLLjyxj9y2gOJ1iaE7z09u17PZf3tIiNbm8hzbHEYw5vgGELt/qqMOTpLxDiq07wLD682KLM2eSz9Flh2UMriEHet2Uw74yJYmTXRuWugzMoj/nitT6WdvUaL/wAAc2g9RDoPQ3cvLnpFjvuqX/9Hyy+Gb4dDzsWVW4dJRLqgPgnZCYQGRSJUTEy/rfx/N/h/zNH5ivNA2dtdOFinU5fdSJySSSJadl8tSuawzHJfPT3RxbxZSpKga/2+nPrORB/gMTMRE5dO8WwlcOsyho01cP8UZ5Jr88Hfc6MHjMY2nwoM3qopqwxrZ33d6xuvHnDmwxuMhB8bTTQx1fB+pdUc9G5nXD8FzBVbMWuROIItdrsoygKVzKvMLjpYB6IfIA5u+bw3v732HVpF6/1ec1iIZCj3LnkE85dDAOTOz5ty78Tki0upl00f87IzWDEzyO4lm07cqKrJoHKS1l6/t1Du9M9tDug+twPaz4MT13t3WpvSNMhDGk6BPrlLyz75T+wz4ad+YshhZ9fToLaMBKSVFtqfc9/7ci1TGg/AV+DLwtuWMCc3nM4fOUwI38Zae75j2091uHrXfX+DJ/Wsx0uP6uX/XACJfHwhoftCj/YX4BT2ThjstfH4FNrd0uyydB3YWaM6ilkj7XPw6Iu6tyBROICarX4CyEI9gw2L9sXQjCi5Qi+vf1bwrzDzOVuaVJ2m73GULo3D6hxWw7ff5iRLUeW+R4lUZoPe6UjO6llw80HXkqAWVeh22Tr/N3/hav/wPtdVJPQ1rfUieQDy2Dfksqvr6TWUavF3x4RfhEsHbLUfNw5pLP585zec6zKu9nYy9WrecmhAQooMIeUNr/gKNXN3FPTF7pUORot3P4OTFgDz56F4XZMiZvmwoc91BATv0yHrQvg1FpIuVS59ZXUGqqXklQieq2ed/u/S0xqDHqN3mJh0YcHP+RyRuHS7cQT/8G75fxy3afAg6Wsy/iL88FNH9CrYS92XtzJ9M3Tq43oFph9XLlYqU4QoYaEoPM9qsfQF4NLLr+pSDykFxPg50dhwCw1HIUErVZLZGRhnKGxY8cyY8aMEs6wpnhI59pGnez5FzAgfAD3t1cX6Szfc56IGas5GZfK81EfkHlpdJGSjgltv0b97OY92e1J3G2MIBzlhkY34KZ1q7a2cSn+TqTJdfDUSbjrazD4lF5+XjAc/hbe7QA/PKCmJV2AvGzX1rMa4+HhwcGDB82vsgo/UOkhnSubOi3+scmZRMxYTcSM1Tz/k9rzH/TuViYt/oe85K5kX+lPTlI3FKNjvfYhTYewcbTtjc299F5sG7vNfLxmxJoy1bW6+sC7cqehOo1PKLQbBs/HwPSDhemlhbY4/J06R/BuB/huohqNVGJmzpw5dO/enQ4dOjBlyhTzCHrRokW0a9eOjh07MnbsWKKjo/n4449ZuHAhUVFRbNu2rZQr1zzqrNlHURTGfVpydMachMJohUHuQVzJsr0xRMGG4kII6nvWx0fvQ2puqlW5ohEZSxPzw/cfJjEzkX7f9rNd/2omutW1caoVBDSF2ZZb/JF5DS4dhK/vsH/eydXwij88/jfUi3BpFe3y2wyIsx2rqdyERsKQks2wmZmZREVFmY9nzpzJmDFjmDZtGrNmqR549913H7/++itDhw5l/vz5nD17Fjc3N5KSkvD39+fhhx+2COlc26hTPf9co4kdZ64we9VR+r65mejEjBLLPzNIjUX/SL/mJYpb63pquYJeRFFhrog5xJY7Z7Uzr1SvNqju4FEPmvdXzUP3/6I2Dvfa2ZHuvU7w52eVW78qprjZZ8wYdRHh5s2b6dmzJ5GRkWzatImjR9WQ0x07duSee+5h6dKl6HR1o09c658yIyePjcfj2XD8MptPxJOSlYebTkPflkE8dlMLdv97lZ8OXKSep55dM9XNSNYcjiUyzI81h+MA0BYR/oc6PkS/xv0Yt9r+Jh2O7Kla23rK1a5Rqiv4hKovgBYD4M7/gx8fsC635unC8NMvJYK2kn76pfTQK5OsrCymTp3K3r17ady4MbNnzyYrKwuA1atXs3XrVlatWsXcuXPNjUJtxmX/AUKI2cCDQEJ+0vOKoqzJz5sJTAaMwHRFUcq26amDZOTk0W6WeukALwOD2odyc7sQ+rYMwtOgPvqY7uE82r8FwT5u5hj3d3ZRl+EPNCks3HCKIZGhrNmuXvPOlnfS0LuhzfvZEnR7Il9bxLK6mZ/qPB1HQ3aKuufAuZ1w6jfrMnPzt258MR501ts01lYKhD4oKIi0tDS+//57Ro0ahclk4sKFC/Tv358+ffqwfPly0tLSqkVIZ1fi6uZ/oaIobxVNEEK0A8YC7YGGwAYhRCtFUYzOvrlGCNo28KVZkBeLxnVGq7EtuC3q257QbdvAl+j5t1mklSTaBWafLwd/yV2/lrxFY0XFv7qIrtnVs5aNZGo03fMXjV03DfKy4Mj3sOox63Lz6oNvGKRchLH/gza3WpepoRS3+Q8ePJj58+fz4IMPEhkZSUREBN27qyFGjEYj9957L8nJySiKwhNPPIG/v79VSOfiUT1rOlVh9hkOrFAUJRs4K4Q4A/QAyrYdlQO467Wsmd6n0oWpbWBbPHWeZORlVEjkbZ0rRVbiMBoNGDyhy3hoegNsnKs2BEVJyY8ntevDWiX+RqPtvuS8efOYN896z+jt27dbpRUP6VzbcPWE7zQhxCEhxOdCiHr5aWHAhSJlYvLTrBBCTBFC7BVC7E1ISLBVpFScJZaOiHjRe5V239oi4tVlsZmkFOpFwKjFMG2f7fxz21UX0V8eV+MJFf9eE/+BrNprAqmLVEj8hRAbhBBHbLyGAx8BzYEoIBZ4u+A0G5eyqSCKonyqKEo3RVG6BQeXfbs5Z1JWM4t5mz4XiXx1E93aModR6wlqAc+dg+eiYcSnoPOwzN/3pRpP6JV66h7F/2yGmL1q2ic3VEWNJS6iQmYfRVFK35AVEEJ8BhTszxcDFF2D3gioMQFKqmokUZbzKpPqMvcgKQMe+Tt7dRqjvrYvhA2zixVS4O3WlknXzsLiQTBwLqx+CoZ/AA2q725ekpJxmdlHCNGgyOEIoGBn61XAWCGEmxCiKdAS+NNV9ahqqptYu4raYsaqk1z/H7j1rdLLAVzYDYtvgbhD8NlN6lzB1X9dWz+JS3DlhO+bQogoVJNONPAQgKIoR4UQ3wLHgDzgUVd4+lR3NKJ2rK+rbuYnSTkQAno8qL4ATEaYE1D6eaY8dd+BvV/AY3tdW0eJ03GZ+CuKcl8Jea8CZd9GqwbhjAnfmjRqqEl1lZSCRquGl973JVyLhsGvw/ndsPRO2+UTT8PCDupis9ZDIKj2eA3VZmpH97Ma42xzSHUV2epaL0k58QyAvk/CsEVg8FJXD89Ohj5P2C6ffAFi/oKNc8CBFe6upF+/fqxda7lu9N1332Xq1KkOX+O1116zOC5vWOfZs2fz1lsOmtQqGSn+DiLFzTZywreOcfNstRGYnQxDFtgukxwDGYlVFlF03LhxrFixwiJtxYoVjBtnPyRLAYqiYDKZrMS/NoZ2luLvICWJXEl59hqN2rLCtwA54VsH6TkFxq2wnZd0HmIPQsZVdW6gEhk1ahS//vor2dnqfgbR0dFcunSJPn36sGDBArp3707Hjh15+eWXzflt27Zl6tSpdOnShcmTJ5tXCN9zzz0AeHsXRgF48803iYyMpFOnTuZ9Aj777DO6d+9Op06dGDlyJBkZJQeNrA7U+sBuzsZZI4By2/yrmcbKCd86Tush8Nh+0LnDdxMgxtJx7409r3Ii9TwYKraTXVHaBLThuR7P2c0PDAykR48e/P777wwfPpwVK1YwZswY1q9fz+nTp/nzzz9RFIVhw4axdetWwsPDOXnyJF988QX//a+6jeZ3333HwYMHra7922+/sXLlSvbs2YOnpydXr14F4M477+TBB9UJ8xdffJHFixfz2GM2QmpUI2TP30FKEv3y5EkzkqTWENgc/MLggfXQ50nbZfKyK3VnsaKmnwKTz7p161i3bh2dO3emS5cunDhxgtOnTwPQpEkTevXqVep1N2zYwMSJE/H0VMOtBwSoXlFHjhyhb9++REZGsmzZshoRFVT2/B2kvGYfe9QW8a9u5idJFXPzy3D8OPgFQXIMz7W+17qMVzD4NXJpNe644w6efPJJ9u/fT2ZmJl26dGHZsmXMnDmThx56yKJsdHQ0Xl5eDl1XURSbo/YJEyawcuVKOnXqxJdffsmWLVuc8RguRfb8y0hZbdsVsoWXcGp1MbfIDdwlNvEKhpD2tvPSEyA306W39/b2pl+/fkyaNMk80Tto0CA+//xz0tLSALh48SLx8fE2z9fr9eTm5lqlDxw4kM8//9xs0y8w+6SmptKgQQNyc3NZtmyZKx7J6UjxLyO2RLc8wlfeRqG6iqyc8JVYoTVAw862866ezZ8Mdp1b6Lhx4/j7778ZO3YsoAr33XffzXXXXUdkZCSjRo0iNdV6u1WAKVOmmHf3KsrgwYMZNmwY3bp1IyoqyuzGOXfuXHr27Mktt9xCmzZtXPZMzkSafRykvKJbXcXaWVSXEYikGhPUSm0IFBOkxUNuuir6SedAxKhrCjwDQe9R+rXKwIgRI6z+Px9//HEef/xxq7JHjhyxOH7jjTd44403zMcFowWAGTNmmL18CnjkkUd45JFHrK47e/bs8lS9UpDi7yDltvnb0f6K9pSrm629tjdykgpgKGJP98+P6agokJMG6YmQfkU1Bem9wCsQ3OupexFIXIoU/2qM3MxFUmsRAtx81JcxFzKvqg1B0nkQF102GpAUIsXfQZzt6qmpZdMtslGSlButHrxDwKu+OhrIsDUa8FdjDkmchhR/B3G6q2ctEUtp85c4jaKjAd88dTSQcaVwNOBRD7yC5GjASUjxdzHO3sO3gOoiutVt7kFSS9DqwLu+6jKak642AhmJ6rveEzyD1E1p5Gig3Ejxd5CymnZKE/3a5j1UXeslqeEIAW7e6ss8GkiE5PPq5vMe9dS5AYNnVde0xlG7DM9VRImB3WqJeccesucvqTQKRgPBbSCwJbj7qQ3BlZOQcFKdJzCp+0IJIXjqqafMp7711ltldrvcsmWLRTTPCRMm8P3331uVM5lMTJ8+nQ4dOhAZGUn37t05e/YsPXv2JCoqivDwcIKDg4mKiiIqKoro6GgiIiKIjIwkMjKSdu3a8eKLL5oD0VUWsufvII6IXFl6vw4FdiuhTHUT3dreyEmqERajgbAio4EL+aOBANzc3Pjxxx+ZOXMmQUFBZb5FXl4eW7Zswdvbu9RY/t988w2XLl3i0KFDaDQaYmJi8PLyYs+ePQB8+eWX7N27lw8++MDivM2bNxMUFERaWhpTpkxhypQpLFmypMx1LS+y5+9EbAmys80h1c28Ul3mHiR1FKvRgD9kJKLTaphyzwgWvvGqeTRQwLlz5xgwYAAdO3ZkwIABnD9/HlB79k8++ST9+/dnzJgxfPzxxyxcuJCoqCi2bdsGwNatW+nduzfNmjUzjwJiY2Np0KABmvy1CY0aNaJevXoOP4K3tzcff/wxK1euNIeLqAxkz99BaruIV5Ta9jwS5xD32mtkHz/h1Gu6tW1D6PPPWyYWHQ34hYEQPDpxHB37DefZB0ZC5jXIbwOmTZvG+PHjuf/++/n888+ZPn06K1euBODUqVNs2LABrVbL7Nmz8fb25umnnwZg8eLFxMbGsn37dk6cOMGwYcMYNWoUd911F3369GHbtm0MGDCAe++9l86d7YS1sIOvry9Nmzbl9OnT9OzZs8J/I0eQPX8HqQqzT02gupmfJBI0OkDg26wb48ePZ9GSn9RAchlXIOEEu3bu5O6xYwC477772L59u/nU0aNHo9Xa9yC644470Gg0tGvXjsuXLwNqT//kyZO8/vrraDQaBgwYwMaNG8tc7coeRcuevxMpixA60lDUpN50TaqrpPKw6qFXJkLwn6eepUuXLkyccD+4aUFBjTF0+Qj41ge9r0VHrLTQzm5ububPRcXazc2NIUOGMGTIEEJCQli5ciUDBgxwuKqpqalER0fTqlUrx5+vgsiev4OU2zXT2Ru4V9MRQ3Wtl6RuExAQwF133cXiz79QYwwFt6Z3796sWLMdMq+x7NP36NO9k7qauNjG8z4+PnajfhZl//79XLp0CVA9fw4dOkSTJk0crmNaWhpTp07ljjvuKNNcQUWR4u9EyuPvX9ORE76S6s5TTz3FlStX1AMhWPTBh3zx7So6DrqPr3/eyHvznlc3nc9KUhuBnHRQFIYOHcpPP/1kMeFri/j4eIYOHUqHDh3o2LEjOp2OadOmlVqv/v3706FDB3r06EF4eDiffPKJsx7ZIaTZx0EcMek43exTkqunFF2JxC5FQzCHhIRYbKgeERHBpk2bCgsrCuRm8OUnH6gTw1dOgc6DVmGBHDq4P38OAfr27WvzHoMHD2bw4MF26zJhwgQmTJhgkRYdHV3OJ3MeFer5CyFGCyGOCiFMQohuxfJmCiHOCCFOCiEGFUkfnJ92Rggxw/qqtYuq7PlXRgMhJ3wlNR4hH8gqZAAAF9RJREFUVJOQfziEdAC//LDTyTEQdxSunTOPBmoTFe35HwHuBCzGK0KIdsBYoD3QENgghCiYyfgQuAWIAf4SQqxSFOVYBevhcmrbBK2zkTZ/Sa1Ao1WDx3kFQU6G6iGUeU1dSKZzV2MKedYzjwZqMhV6AkVRjoPNH/5wYIWiKNnAWSHEGaBHft4ZRVH+zT9vRX7Zai/+5e3h2hXFWqKV0vwkqbUYPMEQnr+K+Jq6ijglplhMIS915FADcVXzFQbsLnIck58GcKFYeuWsaKhllNYY/Zv0L2eSzlRaPeryqEdSy7EaDSSqIwHzaCBQ3Xymho0GSq2tEGIDEGoj6wVFUX62d5qNNAXbcwx2VUwIMQWYAhAeHl5KTV1LdYnCWdr1LqReIMgjiOE/D3fqfYvz4+kfaeDVgKZ+TR2ql0RSKzB4qi/fhkVGAxch5VKNGw2UKv6KotxcjuvGAI2LHDcCLuV/tpdu696fAp8CdOvWrUbZF6pCDBVF4dYfb+W6Bte5/F4v73wZgPWj1rv8XhJJtcPmaOCa5WjAI0CNPVRNcVXNVgHLhRDvoE74tgT+RB0RtBRCNAUuok4K3+2iOlQP7Jr8nb+Be56SB8Cu2F0VunZ5kBO+kuqEVqslMjLSfDx27FhmzCibc+GWLVswGAzmqJ4nT57koYceIikpiezsbPr27cvIkSN57rnnADhz5gxhDUPxMOjp2KY5k8bewfBJT9KsWVMyMrMICQnh2Wef5fbbb3feg1aACom/EGIE8D4QDKwWQhxUFGWQoihHhRDfok7k5gGPKopizD9nGrAW0AKfK4pytEJPUMcoSWRzjbklnjt57WSEEPzfwP9zWn3khK+kOuLh4cHBgwcrdI3iIZ2nT5/OE088wfDhqkn18OHDREZGMmiQ6sner18/3nrrLbp16wa5mWxZ+wt9e0bx65L3QOvGwX/iuOPuSXh4eJQp9IOrqJCfv6IoPymK0khRFDdFUUIURRlUJO9VRVGaK4rSWlGU34qkr1EUpVV+3qsVuX9NoDLNP9nGkjeD+DPuT/bE7nHqPU9eOwlIm7+kZjBnzhy6d+9Ohw4dmDJlirnzsmjRItq1a0fHjh0ZO3Ys0dHRViGdY2NjadSokflaRUcWVug91FDTBh91/YBGR1REPWZNn8AH77wJ2alVvm6g+hqkqhklbuBeji/RFWKZY8wpU3mjyYhGaBw22WTkZpCWm8bhK4fNaY9teqxM95TULbZ9e4orF9JKL1gGghp70/eukgOgZWZmEhUVZT6eOXMmY8aMYdq0acyaNQtQI3r++uuvDB06lPnz53P27Fnc3NxISkrC39+fhx9+2CKk8xNPPMFNN91E7969GThwIBMnTsTf37/kygqR7w0UCLmZdOkez4KPvoLEM6B1A6+CuQF9xf4o5UDG9nExruoR22pwckyOif+2mG2k5aQR9XUUL+540eF73v7T7Qz4bgD/2fwfh8+RSKqCArNPwWvMGDWE8+bNm+nZsyeRkZFs2rSJo0dVq3PHjh255557WLp0KTqd7T7xxIkTOX78OKNHj2bLli306tWrbFsv6j1QvOuDzs08GiDlElw+ClfPVvpoQPb8HaTEDdxL6Dk7PapnCfVwtOc/deNU8+dV/6xieufphHiFlHpeQmaC3bzk7GSH7i2pW5TWQ69MsrKymDp1Knv37qVx48bMnj2brKwsAFavXs3WrVtZtWoVc+fONTcKxWnYsCGTJk1i0qRJdOjQgSNHjtC1a1eH63DgwAHatm1rMRogIxEyrqqB5SpxNCB7/g5SnWPYXMu6RnJ2coniH7nEvn3yyS1PlnoPU7Fwt1b5lJwvkVQ1BUJfsG9uwTaMJpOJCxcu0L9/f958802SkpJIS0uzCun8+++/k5urOlXExcWRmJhIWFiY9Y3scOjQIebOncujjz5amKj3AL9Gakwh/yaqa2jR0cDVs5Dnmo3dZc/fxdjrqVd0RFC0MbrhmxsqdK1DVw6VmJ+QkcDrf75eYpnSGgeJpDIpbvMfPHgw8+fP58EHHyQyMpKIiAi6d+8OgNFo5N577yU5ORlFUXjiiSfw9/dn6NChjBo1ip9//pn333+fdevW8fjjj+Pu7g7AggULCA21tf61kG3bttG5c2cyMjKoX78+ixYtsu3po9Goq4Q9AyA3S40plJ4/0s5KgoZl2xbSEaT4O0h18WYp3miUdZLXHmeunaFFvRY282767qZSz5cun5LqhNFotJk+b9485s2bZ5VedCvHAlq1asWhQ4Udo759+/LOO+/YveeWLVssjvv160dycjnMoXp3dTTgHaLuOOYipNnHxbiy0UjMTPz/9u48uKo6S+D49xASQxZZQlqQNJuDGDALkRBQY0OpIeIKaCDDjDHaoKKWAzMuuCAtPXQVWDY142ARRohC7KAyYrc0ipQRxQxioOlWO9K0SrdA29EgI5uAcOaPe/N4gZfk5eUlb8n5VFGQu+V3yH0nv/u7v4U5780JyrVmvTOLitoKMp7P4MgJZ+7zrw5/RX5lfgtnOnrGd9wKRMZ0CjGxkJgK5yS3y+Wt5h8EoerqOfnXk6n/vr7N1wHY/d1uXvjkBQC+PfYtCbEJlL5RyoFjB3wePzxlOJ/UOy/Frhl0DTf9w01BKYcxxkv3tJaPCZAl//bWXhV/JWiJHyC9Vzr1R53rxUgMX373JXsO7Wny+MrrKoP2vY0xHc+afYKgI+e1CXYz0ms3OhOz1u6vpe5oHQDvfPkOE16d0OQ5Vw+4OqhlMMZ0PEv+fgp0hG9Hvyh+KPehJve9fcvblA4vJTv1dC+IwT0Gn3Xcv3/Q/Kwb8THxgRfQGBMWrNmnvbSQ89v6tHDnxjt9bu8R73u4eX6/fFITUpk9cjbHTx7nklX+D0zxlhybzD0j7mn5QGNMWLOav5/CZYRvS0b8yHd/4ITYBM+/42LiArr20quWUv2P1fRL8n9gizEdbezYsbz55puNti1evJiZM2c2ccbZFixY0Ojrhpk9W2vevHk89dRTAGzZsoW8vDyys7NJT09n3rx5rFixguzsbLKzs4mLiyMjI4Ps7GwefvhhysvLSU1NZcSIEQwZMoTx48dTXV0dUDl8sZp/EHgGOAWY568ffD2/+fw3fh17Un33X24Q26XxkPCLel1Ebp9cfprx08AKBzx56ZP07tabS/sF9gEwpiMVFxdTWVnpmWoZoLKykkWLFrV4rqqiqixYsIBHHnnEsz0YSbekpISXXnqJrKwsTp48yc6dOxk2bBilpaUADBw4kKqqKnr37g1AeXk5U6ZM4ZlnngGceYkmTZpEVVWVM0VEG1nNPwiO/OD0i0/oerp23fCK4MwnhsdHP05mamajbQvyG9cymnP4xOGztqX3On0jpHZLbbTv3LhzeTD3QXrF92rynA2TNzT7PScOmUh+mn/9/Y0JtZtvvpnXX3/dM+na7t272bdvH5dffjmLFi0iNzeXzMxMnnjiCc/+9PR0Zs6cSU5ODnfccYdnhPC0adMASEpK8lx/4cKFZGRkkJWV5VkgZtmyZeTm5pKVlcXkyZM5cuTIWeWqq6ujb9++gLPYzLBhw1oV17hx45gxYwZlZWWt/0/xwWr+QTArZxaHjh9i5Hkjz9p3ZvIvGlpE0dCigL9Xt67dzvr69ozbeWDTAyTHJiMibLx5I9X7qplbPbfxLyQvFRMqPCt/Jcc5g0imDp1K5c7GXTiXXLkk4LIaU1VeRt1fPg/qNX80YDDjbpvR5P6UlBRGjRrFG2+8wY033khlZSVTpkzhrbfeYteuXWzduhVV5YYbbuDdd9+lf//+7Ny5kxUrVrBkiXO/v/zyyz4Xg1m/fj1r167lgw8+ICEhgf379wMwadIkpk+fDsBjjz3Gc889x333NZ7ufNasWQwdOpSxY8dSWFhISUmJZ6oIf+Xk5LB06dJWndMUq/kHwdBeQ1k1YVWjdvVAm4BiJKbZ/bl9cht9fV7C6dk4R58/2tmWeJ6nXf/MXxYNYmNiPfuS4pKoLq5mTt7p0cLXDb6OqqIqq/GbiNTQ9ANOk09xcTEbNmxgw4YNjBgxgpycHD799FN27doFwIABAxg9enSL1924cSOlpaUkJDif9V69nCfqjz/+mPz8fDIyMqioqPA5K+jcuXOpqamhoKCAF198kcLCwlbHFcxpVKzm76fu8d2pO1pHF2n8+3J4ynCfxx887swG2NoXvu9NfY+Tp5pv1/e2rGAZP5xyavDXD77es/3yfpdzYc8LuSv7Lr+u01D7B8jsnckv8pufyM0YfzRXQ29PN910E7Nnz2b79u0cPXqUnJwcKioqmDNnDnfe2bin3O7du0lMTPTruqrq8zN92223sXbtWrKysigvLz9rnp8GF1xwAXfffTfTp08nNTWV+vp6UlJS/I7LMyV0EFjN309LrlzCo3mPkprQuE29pZGure3nnxyX3GR3zTNNHTqVPol9SEtO46OSjxjXf5xnX/dzurPmhjUM7n52P/7mVBdXs6JwRavOMSbcJCUlMXbsWG6//XaKi4sBGD9+PMuXL+fQIWdlsb1791JXV+fz/NjYWM/0zd4KCgpYvny5p02/odnn4MGD9O3blxMnTlBRUeHzmuvWrfPU3Hft2kVMTEzLK4F52bRpE2VlZZ7mpbaymr+f+iT2YepFU/0+fuu0rWz/+3a/FklpyRVpvqdsvic7+P3tvZ8AjIlkxcXFTJo0ydP8U1BQQG1tLWPGjAGcXxCrVq0iJubsptYZM2aQmZnpeWJoUFhYyI4dOxg5ciRxcXFMmDCBBQsWMH/+fPLy8hgwYAAZGRmN1gFosHLlSmbNmkVCQgJdu3aloqLC5/f2tnr1ajZv3syRI0cYNGgQa9asCVrNXyJlKt6RI0dqTU1NqIvh0bA4ykclH7VwZPPXGNN3DGUFZTz47oPUfFXD20VvNzrmwPcHSIxNJNZrVZ8v/u8Lep7T0+8nBGM6Um1tbdASlPGfr/93Edmmqmf3RMFq/gHrIl3avIDJ+knrSenmtPctvGKhz2N8JfhB3Qe16fsaY4wl/wBVFVW16sWsL2nJ7TddqzHGNMeSf4DOHDRljDmtqV4xpn0E0nxvvX2MMUEVHx9PfX29Le3ZQVSV+vr6Vg8Ya1PNX0RuAeYB6cAoVa1xtw8EaoGd7qFbVPUud98lQDnQDfgtcL/aXWJM1EhLS2PPnj18/fXXoS5KpxEfH09aWuuakdva7PMxMAnwNd74M1XN9rH9WWAGsAUn+RcC69tYDmNMmIiNjWXQIOuUEO7a1OyjqrWqurPlIx0i0hc4V1X/163tvwDY4q/GGNPB2rPNf5CI/E5ENolIwwQx/QDvhWH3uNt8EpEZIlIjIjX2CGmMMcHTYrOPiGwE+vjY9aiqvtbEaX8D+qtqvdvGv1ZEhuN7urMm2/tVtQwoA2eQV0tlNcYY458Wk7+qXtXai6rqMeCY++9tIvIZcCFOTd/7rUQasM+fa27btu0bEflLa8sShnoD34S6EO3MYowOFmPkG9DUjnbp5y8iqcB+VT0pIoOBIcDnqrpfRA6KyGjgA+BW4D/9uaaqprZ8VPgTkZqmhltHC4sxOliM0a1Nbf4iMlFE9gBjgHUi0rBw5hXAH0Tk98ArwF2qut/ddzfw38Cfgc+wnj7GGNPh2lTzV9VXgVd9bF8DrGninBrg4rZ8X2OMMW1jI3w7XnAW4AxvFmN0sBijWMRM6WyMMSZ4rOZvjDGdkCV/Y4zphCz5G2Oimtjc0j5Z8m8HIhLj/h21N100x9ZARLq7f0ft50REhotI6+YCjjzdQl2AcBS1N3UoiMhlIvI88JiI9IrGqapFJE9ElgEPuYP5ooqIdBGRc0XkdeA/AFTbuF5nGBKRTBHZDPwcSAl1edqDiIwWkTXAf4lIQUOlzDgs+QeJO5J5CVCFM6R6vohcG9pSBY+IxIjIL3C6xr0P5ABPiMh5oS1ZcLmJ/iAQC/QTkSkQlbX/x4BXVHWiqu6F6HqaE5GxOJ/H/8FZV+SfgJ6hLFO4ibYbOpQuAWpVtRz4V2AHcJ2I/DikpQqeLsBfgVvcGP8FGE10PlJfhDPfy2Jgmogkq+qpaEiO7pPNBcAhVV3sbrtaRHoA0dRcmQF8qKoVwEqcX+aHQluk8GLJP0DuI+WFXps+BNJE5Meq+i1O7fgAMDEkBQyCM2I8BfxKVf8kIueo6j6cifp6h66Ebecdo1fS+zNwHPjC/VMiIv0jtRnPO0b3yaYOyBeRa0VkLfBvOE1cD7jHRFycPj6P7wG3iMhcYDvQF1jirj5osOTfaiLSQ0TWAW8BRSKS5O76HtgMFLlf7wT+CKRE2gs1XzGq6klVPQDOrK0ikgwMws9ZWcONjxgTvZLeSOA7Vf0E+AR4AnhWRGIjqfnHV4wAqnoQWAHMB5ar6nic+bZGu5MuRoymPo+qugNnlcCBwExVHYtTISsUkfQQFTesRMyNHEYSgTeB+9x/X+Fu/xpnacoMERmlqieBvcBlqvp9SEoauDNjzPdxTB7wiaruE5EkERnSkQUMgqZ+juA0byWLyGrgQWAb8CdVPRFhL3+bi/F1nMTY0A5eA/wddyr2CNLkvaqqW4FUYLe76W0gGTjcsUUMT5b8/SAit4rIT0TkXPflWBnwEk5tf5SI9HOT/Rbgd8Av3RrIcOCvIpIQssL7qYUY80TkfPe4hskAewBfikgpTpOXr/Waw4q/MeIkxFTgK2AEzky0QyOhxuhHjP0AVPUPOM0894pIb5wXohcD9SEqut9aca+eA1QD97inXonTsynSKmPtwub2aYLb/tsHeBGnvfsznJrF/ar6jXvMZTjNPDWqutLr3KdxFqoZANzamnWOO1IrY/xQVVd5nbsSmAY8D/zSTSZhJ9Cfo4j09tqfBMR5TUseVtp4r84GGtbcmKWqf+zg4vulDT/H4TjNdn2AE8C9qlrb8RGEH6v5+yAiMW77bzKwV1WvBGYC+/GaBVBV38d5pBwqIt3ddnBwalR3qGpeGCf+1sZ4kdv/veEdxzqgSFVLwzjxB/pzTFTVb9zurV1U9VAYJ/423auq+jRO0h8fxok/kBh7iEg3971NCXCbql5pif80q/l7cZs0nsTp8vZb4FzgZlUtcfcLzgvOqaq6yd2WhDNQ5lKcmv4ItydMWGpjjJcB/YFsVf1bCIrvF/s5WoycvldzGsYxmMas5u8SkZ/gvNjridPVbz7OY+I4ERkFni5wTwLzvE69FqcW8nsgI8w/TG2NcQdOjOGc+O3niMXI6XvVEn8T2mUN3wh1CnjKq61wBE5XxrnAs8Albje/V3FuwIGquhvn5dFVqvpuaIrdKhajxWgxGsBq/t62AS/J6fk/3gf6qzOaNUZE7nO7+aUBJ90bDVV9LYJuNIvRYrQYDWDJ30NVj6jqMXW6bAJcjdN3H6AUSBdnsq9f4YwYjLhh8BajxRgpOkOMoWbNPmdwaxoKnAf82t18EHgEpx/0Fw3tiBqhb8stRosxUnSGGEPFav5nO4UzCdQ3QKZbu3gcOKWqm6PkBZLFaDFGis4QY0hYV08fxJnfpNr9s0JVnwtxkYLOYowOFqMJlCV/H0QkDfhn4GlVjbS5TvxiMUYHi9EEypK/McZ0Qtbmb4wxnZAlf2OM6YQs+RtjTCdkyd8YYzohS/7G+CCOzSJyjde2IhF5I5TlMiZYrLePMU0QkYuBl3FW84rBmSmyUFU/a8M1u6rqD0EqojEBs+RvTDNEZCHOmq+JwEFVnS8iJThLA8bhDDy6V1VPiUgZkAN0A1ar6pPuNfYAS3EWFF+sqi+HIBRjGrG5fYxp3s9wJg47Dox0nwYmApeq6g9uwp+Ks7zgw6q6312EpEpEXtHTq2MdVtXLQhGAMb5Y8jemGap6WERWA4dU9ZiIXAXkAjXuJJLdgC/dw4tF5A6cz9X5wDCgIfmv7tiSG9M8S/7GtOyU+wdAgOWq+rj3ASIyBLgfGKWqB0RkFRDvdcjhDimpMX6y3j7GtM5GoEhEegOISIqI9MdZX/Yg8J2I9AXGh7CMxrTIav7GtIKqfiQiPwM2ussIngDuAmpwmng+Bj7HWXnKmLBlvX2MMaYTsmYfY4zphCz5G2NMJ2TJ3xhjOiFL/sYY0wlZ8jfGmE7Ikr8xxnRClvyNMaYTsuRvjDGd0P8DoTyaO0tENiEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = read_GPS_SONEL(filepath)\n",
    "df.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "data/GPS\\dNUUK_43013M001_NGL14.neu\n",
      "f: 0\n",
      "Trend 1 = 0.81 mm/yr\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:               Vertical   R-squared:                       0.024\n",
      "Model:                            OLS   Adj. R-squared:                  0.023\n",
      "Method:                 Least Squares   F-statistic:                     43.28\n",
      "Date:                Mon, 03 Aug 2020   Prob (F-statistic):           6.22e-11\n",
      "Time:                        10:31:13   Log-Likelihood:                -6062.3\n",
      "No. Observations:                1781   AIC:                         1.213e+04\n",
      "Df Residuals:                    1779   BIC:                         1.214e+04\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        -11.6568      0.345    -33.789      0.000     -12.333     -10.980\n",
      "x1             0.0022      0.000      6.579      0.000       0.002       0.003\n",
      "==============================================================================\n",
      "Omnibus:                     1664.776   Durbin-Watson:                   0.671\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):           325466.306\n",
      "Skew:                           3.782   Prob(JB):                         0.00\n",
      "Kurtosis:                      68.792   Cond. No.                     2.05e+03\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 2.05e+03. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems.\n",
      "Trend 2 = 5.56 mm/yr\n",
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:               Vertical   R-squared:                       0.008\n",
      "Model:                            OLS   Adj. R-squared:                  0.006\n",
      "Method:                 Least Squares   F-statistic:                     3.324\n",
      "Date:                Mon, 03 Aug 2020   Prob (F-statistic):             0.0690\n",
      "Time:                        10:31:13   Log-Likelihood:                -1791.5\n",
      "No. Observations:                 407   AIC:                             3587.\n",
      "Df Residuals:                     405   BIC:                             3595.\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        -39.1406      1.958    -19.988      0.000     -42.990     -35.291\n",
      "x1             0.0152      0.008      1.823      0.069      -0.001       0.032\n",
      "==============================================================================\n",
      "Omnibus:                      389.391   Durbin-Watson:                   1.283\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             7470.254\n",
      "Skew:                          -4.391   Prob(JB):                         0.00\n",
      "Kurtosis:                      22.062   Cond. No.                         468.\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAGQCAYAAAC+gOY2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXwcd33/8fdHsuUjVmQlKzu+YtlOfMSQKI5DCaTYCeEIkHIToNwt4Si/ll9JfwUK5T7aQpu2v3KElh9HuNJwU84QYggJEDt2DuMDS5Zj+YrWURTfiqXv74+Z2cyOZmdnpV3NSno9H499SLszO/vd2esz3/l8P19zzgkAAABAsoasGwAAAACMBwTOAAAAQAoEzgAAAEAKBM4AAABACgTOAAAAQAoEzgAAAEAKBM4AAABACgTOADDGzKzbzA6Z2Rmh2/7czG43s3Yzc2Y2JXKfL5jZh/3/329mN8Vs15nZef7/t5vZn4eWrTezPjN7ee2eGQBMbATOAJCNKZL+aiweyMyeKek7kt7gnPv6WDwmAExEBM4AkI1/knS9mc2u5YOY2fMk3Szplc65b9fysQBgoiNwBoBsbJR0u6Tra/gY10i6SdJLnHM/rOHjAMCkQOAMANn5e0n/y8zaarT9KyTtlPTrGm0fACYVAmcAyIhz7gFJP5D0ztDNp/2/UyOrT5X0WGidouVmFlx/LHTzeyWdkvQdM5tWjTYDwGRG4AwA2XqfpDdKWuBfPyAv+G2PrLdE0h7//wdLLB+UtC902zFJz5HUIumWUHANABgBAmcAyJBzbpekb0j6S//6oKRvSvqImZ1tZlPN7BWSLpD0I/9uP5a0wsxe7S8/S9JHJd3inDsd2f4RSc+WNF/SV82scUyeGABMQATOAJC9D0o6I3T9rZIelnSfpIckvU3Sc51zhyTJOfeQvJ7kN/nLH5DUL+ktcRt3zj0i6RmSlkv6kpnx3Q8AI2DOuazbAAAAANQ9eh0AAACAFAicAQAAgBQInAEAAIAUCJwBAACAFAicAQAAgBQInAEAAIAUCJwBAACAFAicAQAAgBQInAEAAIAUCJwBAACAFAicAQAAgBQInAEAAIAUCJwBAACAFAicAQAAgBQInAEAAIAUCJwBAACAFAicURVmdpWZdWfdDgBjy8x+Y2avyrodADAWCJwnCDM7GroMmdmJ0PU/zbht55jZ183sgJn1m9mvzOzSyDqvMrM9fnu/ZWazQ8vONrPvmtkxM+s2s2tDy8zM/t7MHjSzR83sq2Y2q0x75pjZ1/y29JnZlxLW/RMzu9PMHvHb/9nw9s3sFWZ2l5kdN7NbyzzuKjP7npn1mtnDZvYjMzs/ss7fmNlBv23/aWZNoWUfNbMHzOy0mb0nZvsl92HMumZmHzOz/f5z+4WZrUpqf+i+XzYzZ2btadavhJmtMbN7/P15t5ldmLDuUjP7sf8aHjSzfzWzxtDy/zKznf7nIfPAzsxek/ReG+E23xv6nJ80s8HQ9Xur+VgjbN8bzWyz/9nca2YfMbOG0PI2M/u+/9nebWYvCS0718x+4L+2zszOiWz7T/3P3gkz+3GKtiSub2ZfMLM/+O+Xl6fYXsn1/WXh7+RTZtZbYjszzez/hb7DNpnZMyLrXO2/l4+Z2a1mtrCC5zXFzD7u78cj/vZnjuR5+cvfaWaH/O+oz5rZ1IRtPdPMNviPuz1m+UH/sx7sp++X2tZImffbcyr0GI8krPtm875fw6/dZZF1rrfHv2e31uJ7MC0zu9LMbsvq8ScrAucJwjk3K7hIelDSNaHbvhJd38ymjGHzZkn6jaSLJZ0l6auS/if48vaDo09J+lNJ50h6TNL/Dd3/M5KOSZoj6bWSPmdmK/1lb5D0ckmXSVog6UxJ/1qmPd+VtFfSIn+b/5KwbrOkD0iaJ2m1pCWSPh5afljSP0v6pzKPKUktkr4jaYWkuZK2SPp2sNDMnivpHZKu8B9nhaS/D91/p6TrJcX9OJbbh1GvkPQqSU+VlJN0t6QvlnsCZrZe0uJy642EmU2T99r8P0mtkr4m6TsJP8yfkbRP3vO9WNJVkt4UWr5Z0pslZR5A+p4j6YeV3KHc59Q596HQ5/5tkn4V+txfVOn2amCapL+QdLakp0i6RtJfhpbfKKlP3ufwzyR93h4/mByU9ANJLyux7cOSPinv85dGufXvkXSdpAdSbq/k+s6510W+k78t6b9LbGeapE5Jl0uaLekjkr5pZvMlyczmSfqGpL+R91n9vaSbKnheH5f3+Vgr7zvoDfK+Hyp+Xmb2fHmv39MkLZX0REl/l7Cto/Je43clrPPM0L66JmG90fhQ6DFKdij4bg+/ds65u4IFZvY2ed+dz5T32/ACSSUD8TFQ9e8UpOCc4zLBLpK6JV0Vue3D8r58vybpiKTXyTtwere8L+28pK9LavXXP0+Sk/QaST2SeiW9M7S9mZK+LO9Hb6ukv5XUXUEbj0m6yP//HyV9KbRshaRT/mOcKe9Lfmlo+dckfdj//zuS/ndo2dMkHZc0vcTjPsd/vg0j3Lcvk7Q55vY3S7q1wm3N8fdxi3/9ZkkfDC1/lqSemPt9XdJ7IreV3IclHvvvJH01dP0iScfKtHeqvCD0Ir/d7VV+3z5H0oOh6yYvML6qxPp/kPejG1z/F0n/EbPebyS9qsxj3yTvQOMn8n7sfynv4Obf5f0wbgver/76PfIOYh7Q48HBXP/+j0r6qaTZofUbJT0k74DgJ5LeEnn830t6nqQp/r59q6RdknZVsP/+XN6Pfvi26f723uK/77f7tz9B0m3yPr/bJL0g8v66wW/nEUm/lrQ4tPy5/r5/RF6wVnb/hu77bkn/7f/fKul0ZNv/Len9kfvM8p/DOSW2+TZJP65gPyWuL2mjpJdXsL3E9eUFqyck/VEF29wp6bn+/38p6bbQstmSBqKfv7jnJe875rikhWkfO+l5SfqWpL+PvBe6U2zrecF7L3L7QUmXV9q2Cp/HsO/LhHVLfo/L+/47KOmpKbf1cUlfkfe7e1ReR8kSSe+T93vbLemK0Pq/kfR+Sb/z1/+WvAPOm+V9p/wm+jrK+964QNJ/SfpIZNnPJL05tJ+vl/dbfbyW+3syXOhxnlxeKK+3t0Xeh/mv5X3xPU3SQnnB7L9F7vMUeUH0syR9INQb9EF5PbZL5QU8r03bCDNbKy8o6vJvWq1Qr6BzboekIUnnywsATzrnukKbuNe/j/ztWHjzkmZIWlbi4Z8saYekm8zssJn9zswuT9t2eftqa9qV/VN5pXrMniYvMO73rxftB///BWbWkuKhkvZhnK9JWmFm55mXDvJaST8q8xjXS7pVKZ6//7wfKXGJvsdKPQcn6X49/lpH3SDpFWY2wz91/WzF9MZX4GWS3imvV8/J+6G6S96P13clfSKy/oskXSlppaQXS/ofSf9HXrAS9LQGLpO0wznXJ69nv5A6YmaX+I8ZbvufSLpUXo9eNTxP0iWSLjazM+X9qP6X/7ivkdfTe15o/VfK6yU8S9IBeWdd5KdL3CzvzEibvAPqtRW0I/z5WSnpiHNuT2h5+LM9UVwrL7j8bXCDeWlEsT3E/nu5XV5QJA3/XDwi76ximv3UIS/oer2fXrHdzN44omcR0xb//8VWJj2ujFvM7CHzUtdKPicze1/Cd8rBMo/xv/3v+41m9idl1n2ymeXNbIeZvcseTy1aIu/g+BIz6zGzLjN7j5lZwrZeKO/M2Gx5vzu3yfudPUfeWYJPRda/Vt730LnyPvu/lvQf8j6HexTq3fdTRGY4534v7zvllUFb/LMVT5X3WQ1v+xnyvs8wCgTOk8sdzrnvO+eGnHMn5J3Wfrdzbp9z7qS8o92Xhb4oJK/356Rz7h55P3jB6d+Xyev17fN/+JLSAgr8IPCLkt7nnDvi3zxLUn9k1UflnQqLW9bvL5O8YO86M1tsXk7v//FvL5XDt1DS1fJ6BM+Rl9bxPTM7K0Xbr5YXULyv3LoB59xq59zN0dvN7Fx5Byl/Hbo5+lyD/5tVXtI+jLNP0p3yeg6PS3q+vGAolpktlneK9/0p2hI879klLn9Z4m7lXuuo2+UFBkfkpd78Wt6p/ZH6pnNus/9Z+I6ko865rzrnBuUdaF4cWf9fnXMPOed6JN0h6S7n3L2h+4fXf64eP6X6bUmrzWypf/3Vkr7unDsdWv+j/mfrxCieT9hHnHOP+Nt7oaQHnHNfcc4NOufulvR9ecF/4Gbn3D3OucfkHWx3+Lf/iaS7nXPf85f9o6SH0zTAzN4iL1i+wb+p0td7vHqtImlQzrk/c879dXRFP13pq5I+7Zzb7d88mv20UF6wN09eitUrJf2jmT2tomfwuFLfUSMNnF8i7yBhiaTfSvqJmcU+L+fcBxK+U86Ju4/vn+R1/pwj6UOSvup33sS5Vd7ZmDnyAs3XS/orf1mQV36lvF7eZ8j7TkwaP/Fz59wv/M/2LfLOoH7Sv/51SSvNbEZo/f90znU75x6W9xu1zTm3IXT/Ut8pv5J3sB90Ar1S3tmH8GfzX5xz+6v4nTJpEThPLnsj18+V9P3gqF1e756T96UhSXLOhY/kj+vxL8h5ke2Fe41imdkZ8nrlfumcC+cEH5X3hRJ2pryAKGmZJH1O3hfKL/32/9y/vcfM1tvwgVIn5J3+/oJz7jHn5X8fknRZifWDtj9FXmrKi5xzneWeaxIzmyPvS/FfnXPhvMfocz0zdHs5JfeTmb02ZvDNB+V9CS+Qdzr/Y5JuM7PpJfbDv6n4YKcWyr3WBX6e3k/kBbQz5fV+niMvP3SkDoX+PxFzPRocVLJ+IRfR/+G6RdKfmjeY8eXy3lth0c/qaIW3t1jS08I9dvKC5nmhdUp97ueHt+UfVOwr9+D+WZf3SXq232MqVfB6V8qKB+cNC1Ar3NY0Kx4sNqf8vQr3PU/SH6k4J7nUulPkBVOHVXxAPZr9FARJHwh1gNwi6eoRPq+S31Fm9oHQtm6Iue8wzrk7/HYdc869X17qzpPT3Dct59wm/yD0Mefcd+U9/xeWWHeXH7gOOee2SPqovOBeenxffsw596j/O/Bf8j7bpUS/E3r9M2nh7Z2RsH7a7xQn6Ut6PIh/lWr/nTJpEThPLi5yvUfSMyJH7tMjwXIpB+WlagTOTVrZzKbLO93dJS9/Myzcky0zWy7vvfkHeae3ZpjZktD6F/n3kd9j9h7n3GLn3CJJ2+V9QRx0zt3uhg+Uuk/D94P8bcWtH6SWfEfSa5xztyc9z3LM7Gx5vRq3OOf+IbK4aD/4/+8LBRpJSu5D59wX3fDBNxdJ+prfA3HaOfef8nqmVpbYD0+X9M/+KdEe/7a7LVThJPI8d0R+lMOXUmcnos/B5J2ujEsNyckL+v/dOTfgnMtL+oKSf8QyYWYLJJ3lnAsfjH1R3kDOZ0rq83t9w2Lfo6MQ3t5eST+NfO5nOefenmI7BxT63PtnpxYk3cE/Nf5/JV3tnAtXVtgu6Uz/7Eug8NkeDVc8OC/t4MFS2zrligeLPVTB3V8rLz+5J2klfz9+Sd5B4Mv9A5JA9HPRIu/gJ81+ui94GtEFI3xecd9Re5xzR51z7wttK817KY5TcepdQSQwj17y1XiMMuv+Xt6A1Wp/Nivm/54+WV7qR+BLkl5iXurXInmdVGGZt3uiIHCe3D4j6aPBD5d5ZdrK5X8Fbpb0bjOb7d//baVWNC+H9lvyTuu9PnTEHbhJ0gvM7Cl+r/QH5Q0gOu6ce1RewP0h88o2/bG8U1Q3+dvOmVeWzMzsCfLyUN8f8xiBb0qaa14Jp0Y/8GuTl8sa1/aL5B3Vv9U5N2z0sr+N6fIGdTX4Pbaxo5b9H7yfyvshHVZOTt4X3xvNbKWfOvIeecFgcP+p/mM1SJriP1bwGS65D0vsh7slXeu/5g1m9np5X6xdJdZfKu90fYe8XFnJC1K/F7eyc25F5Ec5fCn1XrlNUqOZ/YV/yvqv5A0M3RCz7iF5AeBbzCu31SovV7cQnJpZk7+/TNJUf3+l/cGspudqeP74HfIGG/2DhvcMDWNmd1hMCcIR+o68XOdr/fdUk5k92T/YKud7ki41s+eZV+3kb+TlX5Zq97PkVUl5vnNuc3iZ8/K9fyDpg/5ne728PPWvhO4/XV6+uCRN898XwbLUn70061f6fim3vv//qxX6DJfYjsnruVwo6YXOuVORVW6Rt8+v8R/vA5LudM51l3tezrmt8gabvcdv7xPl5eZHg6q0z+tLkt5kZsv9ToB3Jz0//7tlurz3uvnbmuovW2pml/nvwRn++3u6vJSNYSKBefSSK/H4U8zshWZ2hr+fnivppfJSk+LWf46Ztfn/P0Fenv93/cfvl/c79rf+9oL0tdGkh43UlZJ+60JpF84bB/R7eZ+3bzjnBjJo1+Tg6mCEIpfqXlS6qsYXIrc1yvvh+4O803675JXtkfyqGpH175D0Ov//WfJ+4B5Rmaoa8nornbxTvkdDl8tC67xa3oCXY/JyQFtDy3LyfrCPyUsJuTa0bJW8EejH/ef9Vyn2z3o9Xg3hbklPSVj3y/J6GcLtvje0/M/95xa+/Gdo+Y6gvfLKbbnIto5Kmh9a/2/kVV94VN6PaVNo2U0xj/WqNPsw5nnNkPRpeT2Ij0rapFCFijL7L6j80F6D9+4l8srInZA3qv/C0LL3Svp+6PoaeUH1I/IGqX1DUlvk/RrdX7Ej+P19+/7Q9aLR9fJyc0+GrvdIWh+6XjRy37//j/3/v6NQ1YrQOu+XN4Bzcbl9q8gI/BLPIamqRnQ0/mp5gxHz/uVWSatLPJdnK1TdQ15JuV1KUVVD3gHp6cj7/duh5W3yAo/g8/vSmLaHL+HX4M0xyz+TsH8S1/efR3T5kxO2l7i+vLKS/fIGcEXv+wVJN/j/r/DveyKyn14cWv9qPT4e4dbw65nieS3273PMf91eX+Z9VO55vVPed1S/vGoyUxO29eyYbQWfiw49/j18WF6nQkeVv0+myhv70C/ve25zZL8u9x9/jn/93/3ndkxeFZr3SmoMrd8qr/PliLzv2nclPPbHVfxbUFRZRI9Xi8mF9nv4+/wTkdfxefLGJkjeGZy3lfgOcAr9tvq317x6yWS6mL9TAQBV5p9t2S8vED4aWfYGeek/68tso13Sl51zf1yjZgIYR8ysU16aZVfk9mdK+pRz7rz4e6IaSNUAgNo5S9LfxQTNM+Xl+t9YbgPOG6xE0Awg+O74TEzQ3CSv5nfZ7xSMDj3OADCG/DzLm+VVBXmpKx4IBgAVMbMOeSkpGyU9xzl3LOMmTWgEzgAAAEAKpGoAAAAAKRA4AwAAACmUrHlZb3K5nGtvbx/bBx0a8v42cHwBAAAwGWzatCnvnGuLWzZuAuf29nZt3Lgx62YAAABgAjOzPaWW0ZWa5FOf8i4AAACY9Aick9x8s3cBAADApEfgDAAAAKRA4AwAAACkQOAMAAAApEDgDAAAAKQwbsrRZeL227NuAQAAAOoEPc4AAABACgTOST7xCe8CAACASY/AOckPfuBdAAAAMOkROAMAUEfy+by2bdumfD6fdVMARDA4EACAOpHP57VhwwYNDg6qsbFR69atUy6Xy7pZAHz0OAMAUCd6e3s1ODgoSRoaGlJvb2/GLQIQRo9zkhkzsm4BAGASaWtrU2Njo4aGhtTQ0KC2trasmwQghMA5yY9+lHULAACTSC6X07p169Tb26u2tjbSNIA6Q+AMAEAdyeVyBMxAnSLHOcmHPuRdAAAAMOkROCf5+c+9CwAAACY9AmcAAAAgBQJnAAAAIAUCZwAAACAFqmokOfvsrFsAAACAOkHgnOSb38y6BQAAAKgTpGoAAAAAKRA4J3nXu7wLAAAAJj1SNZLcdVfWLQAAAECdoMcZAAAASIHAGQAAAEiBwBkAAABIgRznJAsXZt0CAAAA1AkC5yQ33ZR1CwAAAFAnSNUAAAAAUiBwTvL2t3sXAAAATHqkaiTZsiXrFgAAAKBO0OMMAAAApEDgDAAAAKRA4AwAAACkUJXA2cw+b2YPmdkDodvOMrOfmdkf/L+t/u1mZv9mZrvM7D4zW1ONNtTE8uXeBQAAAJNetXqcvyDp2ZHb3inp58658yX93L8uSVdLOt+/XCfp01VqQ/XdeKN3AQAAwKRXlcDZOfdLSQ9Hbn6+pC/6/39R0gtCt3/JeX4jabaZzatGOwAAAIBaqWWO81zn3AFJ8v/O8W9fIGlvaL0e/7ZhzOw6M9toZht7e3tr2NQSrrvOuwAAAGDSy2JwoMXc5uJWdM7d6Jxb65xb29bWVuNmxdi507sAAABg0qtl4HwoSMHw/z7k394jaVFovYWS9tewHQAAAMCo1TJw/p6k1/r/v1bSd0O3v8avrvFkSf1BSgcAAABQr6oy5baZfU3Sekk5M+uR9D5JH5d0s5n9maQHJb3UX/2Hkp4jaZek45JeX402AAAAALVUlcDZOfeKEoueHrOuk/QX1XjcmuvoyLoFAAAAqBNVCZwnrBtuyLoFAAAAqBNMuQ0AAACkQOCc5FWv8i4AAACY9EjVSNLTk3ULAAAAUCfocQYAAABSIHAGAAAAUiBwBgAAAFIgxznJZZdl3QIAAADUCQLnJB/7WNYtAAAAQJ0gVQMAAABIgcA5yYtf7F0AAAAw6ZGqkeTw4axbAAAAgDpBjzMAAACQAoEzAAAAkAKBMwAAAJACOc5Jnv70rFsAAACAOkHgnOS97826BQAAAKgTpGoAAAAAKRA4J7n6au8CAACASY9UjSQnTmTdAgAAANQJepwBAACAFAicAQAAgBQInAEAAIAUyHFO8rznZd0CAAAA1AkC5yTXX591CwAAAFAnSNUAAAAAUiBwTrJ+vXcBAADApEfgDAAAAKRA4AwAAACkQOAMAAAApEDgDAAAAKRAObokL3tZ1i0AAABAnSBwTvLWt2bdAgAAANQJUjWSHD/uXQAAADDp0eOc5DnP8f7efnumzQAAAED26HEGAAAAUiBwBgAAAFIgcAYAAABSIHAGAAAAUmBwYJLXvS7rFgAAAKBOEDgnIXAGAACAj1SNJPm8dwEAAMCkR49zkpe8xPtLHWcAAIBJjx5nAAAAIAUCZwAAACAFAmcAAAAgBQJnAAAAIIWaDw40s25JRyQNSjrtnFtrZmdJ+oakdkndkl7mnOurdVsq9pa3ZN0CAAAA1Imx6nG+wjnX4Zxb619/p6SfO+fOl/Rz/3r9ufZa7wIAAIBJL6tUjedL+qL//xclvSCjdiTbu9e7AAAAYNIbi8DZSfqpmW0ys+v82+Y65w5Ikv93Ttwdzew6M9toZht7e3vHoKkRr361dwEAAMCkNxYToDzVObffzOZI+pmZbU97R+fcjZJulKS1a9e6WjUQAAAAKKfmPc7Ouf3+34ckfVvSkyQdMrN5kuT/fajW7QAAAABGo6aBs5mdYWbNwf+SninpAUnfk/Raf7XXSvpuLdsBAAAAjFatUzXmSvq2mQWP9VXn3I/N7G5JN5vZn0l6UNJLa9wOAAAAYFRqGjg757okXRRz+2FJT6/lY1fFO96RdQsAAABQJ8ZicOD4dc01WbcAAAAAdYIpt5Ps2OFdAAAAMOnR45zkTW/y/t5+e6bNAAAAQPbocQYAAABSIHAGAAAAUiBwBgAAAFIgcAYAAABSYHBgkve8J+sWAAAAoE4QOCe56qqsWwAAAIA6QapGki1bvAsAAAAmPXqck7z97d5f6jgDAABMevQ4AwAAACkQOAMAAAApEDgDAAAAKRA4AwAAACkwODDJRz+adQsAAABQJwickzzlKVm3AAAAAHWCVI0kd97pXQAAADDp0eOc5N3v9v5SxxkAAGDSo8cZAAAASIHAGQAAAEiBwBkAAABIgcAZAAAASIHBgUluuCHrFgDAhJfP59Xb26u2tjblcrmsmwMAJRE4J+noyLoFADCh5fN5bdiwQYODg2psbNS6desIngHULVI1ktx6q3cBANREb2+vBgcHJUlDQ0Pq7e3NuEUAUBo9zkk+/GHv71VXZdsOAJig2tra1NjYqKGhITU0NKitrS3rJgFASQTOAIDM5HI5rVu3jhxnAOMCgTMAIFO5XI6AGcC4QI4zACBT+Xxe27ZtUz6fz7opAJCIHmcAQGaoqgFgPCFwTvLZz2bdAgCY0OKqahA4A6hXBM5JVqzIugUAMKG1tbWpoaFBQ0NDMjOqaogJYcB7oJ4ROCf5/ve9v9dck207AGACc84V/Z3MSF0ZvZEGnfUSrPIeqG8Ezkk++UnvL4EzANREd3d3UeDc3d09qYOEcOrK4OCgtm7dqtWrV0/qfVKJkQad9RSsRt8DpC/VF6pqAADqxsmTJ7NuQqaCCWEChw4d0oYNG6g4ktJIZ6Kspxksm5qaEq8jWwTOAIDMtLa2Fl3fv3//pA4Sgwlh5s6dW7gt6HVEecGBh5lVNBPlSO9XCwMDA4nXq4lSkJUjVQMAkJm4oGCyp2vkcjktXLhQhw4dKtxGr2M6I52Jsp5msByrHud6Sk8ZT+hxBgBkpq2tTWaWdTPqzlj2Oo5ncT2muVxObW1t6u3tTd2TWi8DAyWpr68v8Xq1xOVSozx6nJN8+ctZtwAAJrRcLqc1a9Zo06ZNkiQzU3t7e7aNqgNB6sDQ0FDmqQP1qlSPaWdnp+655x4551L1pIa3Y2Zas2aNli1bNobPpFg0z79Wef/kUo8MgXOSRYuybgEATHgtLS1FtZxRX6kD9SpuQF9/f3/hIExKV5UivB3nnO655x61tLRkts+nT5+eeL1aOKsxMqRqJPnGN7wLAKBment7NTQ0JMkLXDhl7Mnlclq1alVNA7jxPDgsOqCvqalJ99xzT9E6aSbViaYLZfkeDF6HoD1mNmwAbbVEe5jrLXCu1/cmPc5JPv1p7++112bbDgCYwEhLyEY+n9ftt99e2O/r168fVz3b0V753t7eYZPorFmzpuxzCtKFwukdWbwHoykjkhfEb968uSY94NFAeceOHVqwYEFdvAfqeeAigTMAIFMjSUuop8Fc9ayzs1M9PT1auHDhsLzd7u7uQk//0NDQuKxmksvlitrc2NhYSLuYP3++WlpaUm1n2bJlamlpyfQ9FU0ZCYz0tSn3GYk7OKiX90BcGk49tEsicIKP6sUAACAASURBVAYA1IFoAJSkWr1REy34jj6fzs7OQr5vUNouy0Fv1RZ9vsEBWHd3t3bv3q39+/fr4MGDqXrS6+G9ED7zIo1uCvpyn5Hg+eZyubpLhZDq+ywUgTMAYFypRm9UPZ8KDlQSzIWfj+T1th45cqRonZ6enqLAub29Xbt375ZzbtxVMwk/34aGBrW3t6u9vV25XK5oGvc0vbX18l4In3lpamrS5s2bC4Fjpa9N0mckuu8C9fAeCL/n63VwLIEzACAzI+npq0ZvVD2fCpYqD+bCz0fyZmCMmj179rDbzKwQOI8n0devq6tLe/bs0bp160qWcyuVthJXzzir90L4zMtoUkeampoKr230MxLdd4Gs3wPR93xHR0em7Skls8DZzJ4t6V8lNUr6T+fcx7NqS0m33JJ1CwCg5rI6TZ3Ua5ikGqXa6vlUsFRZYJ/P53Xs2LFCoFTKrl27igZ/xVUzySJgTHr/lVoWvH7hg4XBwUFt375dDz/8cNE2pk+fPixt5ejRo7roooskJZdlyzKFo5L0pbB8Pq8tW7YUDoiCAHTbtm1qa2srVBGJvleyeA8E+7epqUk9PT1FBzDB61VvZ4QyCZzNrFHSf0h6hqQeSXeb2fecc7/Poj0l1cmLBAC1kuVp6rhew+7u7pIBdD6fV3d3tySptbVVx44d07FjxySp4jbXe53kaGBYanKKcGWMcqK9qfVw8JD0/ovrgRwYGCi8XuvWrdP27duLetfjetpbW1vV1dVVdFtQQUKSdu7cWbRs586dhWVj/dmoRqAeHWS4Z88e9fX1FT2P2bNnD5uRcKzfA/fee6927NhRdr2szwJEZdXj/CRJu5xzXZJkZl+X9HxJ9RU4f+EL3t/XvS7LVgBAzWSZshDXaxg97R4OokoFiN3d3SMqpTbSHr1aCPe8BcHheeedpx07dsg5py1btsSWJAtXxkgjHBiN9cFDXKpENE3id7/7nVasWKFly5YNWxaUiwtm95PSTUfd1dUVu15Qq7lUz2vwuNLYfDbiDhSCdpc6E9PZ2VkI/JcvX65ly5apv79/2HYDQRC6dOnSooli5s+fr5UrV47Z56Gzs7Nk0Dx9+vSidJs0tbjHUlaB8wJJe0PXeyT9UXQlM7tO0nWSdO65545Ny8IInAFMcLXodUzbaxYEbvfdd9+wkf3RQCWcVhAVrBusl3UPcqW9hnEHBQ0NDUUBXalet0qnY968ebOWLl06ZtU1gn0xMDBQCJTCFT6iB09Hjx7Vpk2btGfPHjU1NRX2Q3CRvMA2HPSVUyq4bmpq0oEDB4bdHg7Uop+NavQIl9pG9EAh/Bx3796tK664omj9cPqJJG3atEk9PT2F/VvKwMCAVq1apaNHj2rfvn1asGBBIW1lrPT09JRcNm3aNJ06daroIKleDnCl7ALnuAz0YYlZzrkbJd0oSWvXrh15XRYAQKxq9joGqRRBpYY0k2rkcjk1NjYOuz0axJeb1ayzs1MPPPBAYQKLtKfV43p6R7sPKj29H9drHHeQED2oyefzOnjwYEXt6+vrKwRbLS0tNU1FiFb6CAsqfATvvzvvvLPoICB8IFVJybQVK1bokUce0cDAQGJvdHNzc8nge/ny5YX9EP5sSKNP3Uh6f5RKx5G8g4VodZC44LNc0CypcACza9cuDQ4ODst9r5XwAcPs2bNLtjXoMQ+C5noroZhV4NwjaVHo+kJJwxOTAAAjGjxViWqkLMT1mpYrBZbP57V9+/ZhZdOmTZum9vb2Qi9yf39/2VzI48ePF/4PBoldfvnliffp7OwsnP4PjLYu9ObNmys6vZ/P51OlG8RVPEjqhS+np6dHAwMDNU1F6O7ujg2aJa+n/Ec/+pGam5u1cuXK2IOnQCV1hg8fPqwLL7xQv/vd7xLXi77nwh577LHC/+HPxrZt2yp+baOfzaTUqDTvA+nxlJe0Zs6cWfT5kFT0eRocHNTWrVu1evXqmgXP0QOGuXPnlr2Pc04HDhwgcPbdLel8M1siaZ+kl0t6ZUZtAYBMpAmIo/Vcwz249VJ/Viqda/voo48WRvNHJ2D4xS9+EVsF4tSpU4Uf9sbGRs2cObPi9uzfv1+dnZ0lf3Tz+fywoFkaXV3o6IGDmampqSn2+Qf3KdUjK3m9okGA55zTr3/9ay1YsKCQ7xqXI57WwoULC22MK1k2Wvl8ftiAvLCgV/HIkSOxA/pG87i/+MUvYkvvjWabwWcxbVpTqc9mODXKzHT48GHdcccdklSUNtLQ0KCFCxfqwQcfLNw2derUYekZcVpbW4uC8Obm5mGBc9ShQ4d06NAhtbS06JJLLqn690g0DSXta75//37l83lSNZxzp83sbZJ+Iq8c3eedc1uzaAsAZCGpFFup0eZDQ0Pavn27zj77bLW1tQ37McpyutxSubb5fF75fF5mpkWLFunUqVNauHChBgYGUs2MNpKgMNDV1RUbOAc9w3GPP9KBSHG9v83NzdqyZUvJA5to7eWoaHrKqVOn1NXVVZTvum7dOm3atGnYgLA0wgcOHR0dVe9tzopzTidOnCi5PDr4LCw6CUi5yh6llOpZzuVy6ujoKAxYLBVAtre364wzzii6befOnWUPCFpaWoYFzmnSNwL9/f36xS9+MSyfWhr52a2gXGJDQ0PiGZKWlpbY93G9TAMeyKyOs3Puh5J+mNXjp/LD+m4egPGrVCm25ubmxCBo//792r9/f+FHPPxjtHv3bkmlR+DXQvBjWo5zrtB7dujQoYoGfCcFQUlmzJgx7LY0vbzB86lkH8YF2+HXMW5wX6l6uoFTp07F3u6c03333aczzzxTJ0+eHFHQHNT5DezZs6cqOd7V0NzcrNOnT2twcLBsbnsp06ZNKxkcJw2oXLJkSeKAvb6+vmEBbZy2traiz+axY8fU2dmpvr6+whiAUhobG9Xe3j7sdXXO6fTp04mPe95552nz5s1l25ckrp7zSM9uRTsI4tJGAiN5H2eBmQOTjOD0IACkEf1hlbwAOu2PRxBUtLe3F06JO+cKAfhIyrNVqtTUvWmET0GXUy5YKGXevHmF/4MA/9ixY4m9vP39/br//vtrkvoyMDCgjRs3Snr84GbNmjUVVYgIBD35IxXdB8H2Kn3ecb2Q+Xy+4mofYUn5x2kNDQ2V7MEsJW7K6WjgHgS9leyn4HOZ1oIFC9Tb26vDhw8PW5a0b5YuXaqBgYER572HRQ8ER1q2MpznPjQ0VDZlJKoepgGPInBO8qlPeX/f+tZs2wFgwsnlcjrrrLNGFfwEs4BFf5TLDcqrhnw+r61btxb9KIZ7T1taWnTy5MmSvaZjIQh6wnmhaacVjvYQhydfievRT5OaEE6/6erq0vz58zVlSu1/hpNSE6IqmWwifOBkZpo9e7bmzJmjP/zhD1UJ3tIKypeFjST4nj179rBe1ujkKMH7u1zwOJqBm5UcVIYFAeZI897D+vv7i977QapFOB++XOpG+DMzUuEKJ/WCwDnJzTd7fwmcAVRZZ2fnqIJmyftxKzUS/9FHHx3VtpOUSncIn34eyWnXM844ozATYDUEP+7hHl3nXOqeyPDp9fDBye7du4vKZKWtjBGVdoDUaAOhSnt/w2XRkoKj6Ax1fX19I9oPo1Wtg7M5c+YUDeTs7e0tmVKRlAufz+dje4traerUqYVZN9etW6fu7u6KermjgnKB0bNKS5YsKQTo5VI3RnPwEEgq0ZcVAmcAyEAl5aRKSTrFn8/nde+996aa2KCSQT/RkmvVMn/+fEkaFjiPNGgMcsXjBlkGOapHjhwp+cNuZiUDj2ACjr6+PrW2thYGANZKQ0NDRdsPTxoSqCRtYfPmzcOeW1xwFJduNJ7t3LmzkIbR0dGReBBXahxBuRz6QDV6hcMee+wxdXV1qaurS+eee65OnTpVqE1eziWXXKKjR48WfVaCqivRFI3wgVE493vz5s1qbW0t7JegZF5SDn8aI81xryUCZ2AS+darvqX7v3a/VOnvnElLn7FUr/7Jq2vSrskoaQKAatmxY4dmzZqVWAe1kkE/SdNej1Zzc/Ow0+JS5UFjYMqUKSUPLMr19M6aNUtHjx4t+xhdXV2jDgzSCNcVLqe1tVUXX3yx+vv7i6aoruQMQDBYNfzcwqkJ4fJsLS0tmfQyj1Q40I8G/cFzDU/vXcrJkye1cePGYQF0uUopgUre05VMACNVluphZhoYGNCsWbNil0dLHpY6qxDcvnv3bi1atKioDaM5SNi5c+eYTM5SCQJnTFqf+6PPaf/vmHcnFSd1/bRLX37Wlwmeq6SSYGg0tm7dqpaWlsRczLSDfkrVaq6Ghx56KDZQGel+Gk0wlyZoDtQ6aK5Ua2troexZS0uLent71dPTM6L9EQTdkhdkNjU1aePGjWWrQtSr+fPnq7m5udCzmvReLvf8goOv6EDc0dTWLmW0KV3lHDt2bFgu8s6dO4tmdty6dWuqA/1w9ZxA0r7I5XJavHhxybM2cTMmZo3AGZMSQfPIPPirkQ1awXC1zEEOO3nypDZs2FCyJzk8IUPSpA7VGOhTiplpxowZ46rnMmwsep3TCOqBB4IAuqmpaUSVO8J1i5uammqekjJSTU1NOuecc7R3795hr0N4gpd58+bpnnvuib3/E5/4xGG57GnEDcRdvHixJK8azEgH+lXL3LlzEwPeUhU/jhw5op/97GeaM2eOHnvssZq87j3/3qNtd27Tr4Z+VXKdhlkNmvK/pmjt2rVVf/yRInBOcvvtWbcANXLwnoNZN2FcOveP09feRWlpS4lNmTJlxKXYwpIqJQQ9SuVynKsx0KeUWbNmaeXKlTp48GDiY+RyOT388MOFCh5LlixRa2urdu7cWZUSZiNVzaA5PFtgmJlp3rx5sWkmwb4olXcbpOr09PRo2rRp6uvrS7W/Ojo6CvcNTzddKzNnztSqVasqDmAHBgaGBc2zZs3SihUrCr3uwYRBca/VE5/4xMJAuD179lT8PIP66WOR716p1atXa/Xq1SMaLDiSwZ6HvnJID//oYalKJ9SGjg7p/o/dr8Xti3XJdZdUZ6OjROCMSemcNefQ41wJcpyrKs2EIZJ00UUXjagXLE7SbHhBz2SSWo5unzZtmnK5nNavX18YmFaqDevXrx8W5Le0tJScvjuNpEkZxlpbW5uOHz9eOAMQnqmuv78/NnAOV/goZdmyZVq2bJk6Ozu1d+/eVG3p6+srnGk4efJkzXvWV61aVXge7e3t6u7uLqSFNDQ06Pzzz48d7CkNTyt50pOeVHh/hN/b0bzmFStWFO27xYsXa9++fRVV6gh6bau1f2bOnKlzzjlHU6dOLfl80+rt7dWqVaskqez3SKmqNsd3Hte+/9in0wdHfxA/Ur//5u8JnMeFT3zC+3v99dm2A1X3xt++cdKmazQvbNZLb36pFl22KOumTFrhPEgz05lnnlk0eGvq1Km68MILC71gow2cp0+frrvvvlvNzc2aN2/eiGaIq6TnKaiQIXmnfMv1bp555pmSvADn4osvThyAGBfk53I5XXHFFSMuwVUvQXN46vW4MwClDrjSVh7I5/NlB72FBVUaaiVa6SP8Hgte5+i+WLBgQSFlKOjhDc5AnHPOOZo+fXrJnvdge8FzMrPCAWHaahhJqnVQcfz4ce3Zs6eQ8hEVl0Mdd/AXLpkXfe/UQzBciQtefEHWTSggcE7ygx94fwmcJ6Q3/vaNWTcBk1Q0PaK7u7sogGhrayv0gvX395c8fR/V0NCg6dOnD/sBPXnypE6ePKkjR44UeizTzHwWzmtOm5O9YsWKohJ44clH4kRnBkvqeZ4+fXrJ7QR1d8er5uZmXXrppbG9pIG4gWeNjY2JZxPCkuoSj4Voj2ya2SajB0rR6y0tLeru7lZ3d3dhKvqkmeba29u1Z8+eYTn9SdUwpk6dOmaDeQOlKn9IXpAc/j6YP3++pk+frq6uLu35yB4dv+/xz//v9fuxaXCNTD9ruq762FV109ssETgDQCaiAUC4UsGBAwfU2dkpqbhWc3Nzs5YvX64DBw7EnrKvZErbchU0Rlp6bsGCBUXb2LJliyQvaFq+fLl27dpVCFByuZwuvPDC2B7kiy++uCj9Is3Uu2kDyDjVTEMYybbSnAEIH3AFNXorOXPQ1tY2rG1jlaYyf/58zZs3r9BD3NDQoKVLl2rz5s2F6yOZWjk4YArep+Xe17lcTh0dHerp6dHChQuLqmGUUuugOa78Yfjsw8//8efa+amdGjoW/1ncpm01bV8tNcxq0NM+8DSt++t1WTclNQJnAMhYLpfTkiVLCqeQnXO65557NHv27KL1jhw5oi1btqi1tXXUjxn0tnV2dqqrq0szZszQypUrC4HESAcDbt++XWeffXZhMFa4F6+pqUnnnXdeIW8zKf0jnH4hlZ5wInqfFStWpMoLjdbGdc6l7tlPMn36dK1evbowSCw4YCjXprSvaZp89KT7rlmzppCu0djYqMbGxhFtK63W1lYtXbq0cAYlPFgvXDKv0tShsGhlmKampqIZAMOCg7nBwUHl8/lCqcbgUuvSb4E0qRIP6IExaUu1WZOp9dmtmvunc4tunz9/fuG7QdKoX/esEDgDdSR8ajxNoICJo729vajX2TkXW6JtaGgo9pRyJb2c8+fP18qVK9Xf31/o0e7r69P+/ft15ZVXKpfLjXhWuKA3PChlFg1oHnjg8WBgcHAwsUbrSILENIMYGxoatHjx4mFB0pQpUwrtNbMRHTisXr1ay5YtGxYQLliwQHfffXfJwHysZkiLti38HqimpLMJSakXI32scE980myHSXXL494TlYqmSkwYJs184kwt/rv4vOu0goGKUnw60nhA4JxkxoysW4AJLDrNcT6fLzo1vXv3bl1xxRXj9ssFlYnrDVy5cqXmzZunrq4uPfLII5I07BR3uCxb+LakGfeCIC1u2u9wINve3q7Dhw8nzjoXVDIIqhqEZ5obGBgoyuWOy7Ht7u6u6kFiXDpClHNOAwMDw3qnly5dWhRUbt++vSglJqiJXGrmwXCFhnBAGHzWly9fHluuLKl+di2E25bL5Uqm/oxG8L4Zq++v4DmFy+bFpW2EDwjDg+ek0mdADn3lkB7+n4el+qkyVxXWYFpy1ZKiakn5fF633XZbTR7vwIEDyufz4/43jcA5yY9+lHULMEHFTXPc3d1d9GPvnNP27dt1+eWXZ9hSjKW4nspcLleorlHuFHdwW7kJL4KDtDlz5pRcnia/Oei5llTU2xceeBXtUYwObnPOJeakVip6ABIVHFQEbZs1a1Yh3zUc9ErDByMG1T9KCed3B8Kf9SBto6mpSU1NTYVALeuzSytXrtShQ4eqWn84SDdKmrWyFtJO6NN3a58e+upD2nps65i1bUyY1PqEVr3ytleOeL/Pnz+/7IHU0qVLKy6XV+3PelYInIEMxJ0uPHny5LD19u/fPyGO0JFeqVPXaU5xh3vdynHO6fDhw8Nub21t1d13350qTSE8yC8avMflLgan1IMqCEF93mr3tgYHIHfddZdOnDhRuL25uVnt7e1FbQvqG8cJagkHQVjQq19KXFAQ/qw757Rz5866O5MUvC7RHvbRGutAae9de/Xt135bfX94vOd4vOYJRzVMadDqa1frRTe9qGaPkaYkX2trqy6++OLCaxo+8Dx69GhiID3WZ1ZqhcA5yYc+5P1973uzbQcmnHAeZlBLtNQPVlIOKBAnbX5y3KyEleS7Hjp0SPl8viiPtFzOai6X0759+zRlyhRNmTJFq1atqsn7O5fL6YILLih6PsuXLy87UUh0G+EJV5IGTEZP+weiqSP12uuWy+V0+eWXq7OzUzt27BhW5SEqTU59JaXyknzrVd/S/V+7X6rNxJWZiUuVyFJSST7Jez3DQbP0+IFnuIJOVC6X05lnnpn5mZVqIXBO8vOfe38JnBERnDYfSUmoaF1b55z27NlTq6aOWjRFAPUvCPjuu+++mlcJSJrOO869995b6JUaGBio6en88FTT4VSMSsSlmgTBczhwXLJkSckzBdHc9XrudQt668ul6gSlBcOTjxw8eLCQ6x6UUiv1ut76t7fqrhvu0tDABIuGJU1bNk1LP7q0aCbDcjXT60HSwNqlS5cmvp7hoDsYd3HGGWdMyN8NAmegQnGns+JGb3d2dhZ+sCVvFq6BgYFhU5o65xJ7d6ZOnVrlZ5BOPp8vOnXb0NCg9evXT7gvwYkql8vpyiuv1B133FH1gV9hpXpaS9m3b1/R9Vr3wCalYlQqrnpDmhrEcbnr9Sw48Aqm2o57/zQ1NRUN/Mzlctrwzxv06w/+Wo/1P6atmli5w01nNumZ//TMwkQcnZ2dibMwmpnWrFmTqmOlXjonSlV2mT9/vtauXZt432hu+UTpXY5D4AxUKO50VnT0drhX+dChQ4nbC0pjlcoNC6ophNX6izZa4UPyniNpI+NPtQd+NTc369ixY4WexjVr1lT0nliwYEHRe73SwDtrafO5k+43HoTb+/krP6+9v9hbtHw8T7oRGE3ecPhg6PDhw8MOLpYsWZLqgC1uoHhW75O4mSkl6eDBg2XH2oQPKrM+AKg1AmegQk1NTcPy+6I//nFlvuLMnDlTT37yk5XL5XTixAk9+OCDw9Y5cuSIOjs7C1/C5b5oqxFURyt8YPyq5sAvM9Oll14qaeSTFwTTcT/44IM644wzYmv9jhfjLRgOI2949ILXf+PGjcVtSDHLZSCprvRYCwe/4YOBtGeFxvPnoRIEzknOPjvrFqDOBAMgokHl0NCQ+vv7C8HEwoULy/Y0Syoa8X/q1KnYdY4fP17ovV62bFniF221ei/iKnxIGtGUuMheLpfT2WefPerAefbs2UWDAEfqoosuKgTQqJ5NN27ST6//qQaOjM1kKmOpeWGzXnrzS7XoskVZN2WYcPWVSs/CpC2fN1aC4Defz+vQoUN10656QuCc5JvfzLoFqBNBL+6xY8dKnvIOD/5Zt26dzj333Nge5DDnnLZu3arVq1eXDbZ7enq0bNmyxAL+1eq9iNavlbw8t8nQmzBRlToNG1auUsLSpUtr0TSUUCivtqtPmoAngKJ5w+NVtPpKJd+T9ZriUK/tqgcEzkAZaWpbSiqaMa2S0+JBSa+Ojo7CkX6cYJDhvn37CqPdg7xjyfuiG23vRXCA0NraOmxZmjw31K9cLqeOjo7EcnOlgubm5uaKS7kh3q1/e6t+82+/0eDJCTYNna9xWqOe/FdPVsffdEyqoGs0aQr1muJQr+3KGoFzkne9y/v7sY9l2w5kqru7O9XAqqD0kJnpwIEDFT3G4OBgoce6oaFBZ511lk6dOqXHHntMp0+f1vz587Vs2bJCjdWwrq4udXd3FypehEf99/b2Sip/Wj2fzxcmpShVgirr/DuMXqlR83Gampq0cOHCCT06vlr23rVX//2y/9aRniNZN6U2RjgbHe8bTEQEzknuuivrFiAj4TrNQY9uOY2NjZozZ45OnDhRmEq3EuEJEubNm1c0bfKDDz6otrY2dXV1xd43WvHi2LFjuv/++wvLr7zyypI/Yml71KXkOp+of3Gvn5lp0aJFw9KKZsyYUbYE1UT3uT/6nPb/rnal/LKWJm84/P2wYcOGcVGPGKglAmcgolx9zlJOnz49qsFX4WL5bW1tw6b23bZtm44fP564jVJB8H333acrr7yyaL3gNGraHnWpsh5L1J/o6xdMnxucmQh79NFHJ1xqzkSedCMQpEpc9Q9XVWV79VT1AagHBM5ASD6frzhoTjO1cRrOObW2thYGYEXrN5cLmk+fPl0yCA5PsFJJD3MYI6vHv2gOfHj63Oj7uF6nho6a6GkSjU2NuuClF4yo1nA11FvVByBrBM5ASG9vb8U9zWeccYaOHKnOj3ZfX5+2bNmixYsXV9yOgwcPluwRnjFjhrZt26a2trbYCVzKMbOiIAvjU6mR8kFVgO3bt+vAgQOZTw1dqCbxh8pTnsaD6WdN11Ufu2pcVJOgugJQjMA5iV/FAJPHSFIR0gTNM2fOLNtjHAiC2nKlwaKSer37+vrU19dXGHhYKefciPK2UX9KjZTP5XK6/PLLazor5Zef9WV1/TQ+T3+8s0bTE17+hMx6hmuJ6grA4wick9x0U9YtQA1FA4R8Pl9y2uuRamxsVEdHh1paWnTbbbeluk8w61Rvb29FPdmnT58uu87Q0FDJcnfldHd3U2FhEqg0SJrQecMmLX3G0jGZhQ7A+EDgjEkpOsNeR0dH6mmyKzF37txC7dtLLrmkUCXDzDR79uzYXtw1a9ZISteTPZbGS84rRm8iz0A3USbdAJANAuckb3+79/eGG7JtB6ounOc7ODiYOCnEaKxcubLw/7Jly9TS0lLo5e7v7y/5uGlL4I2FaLUPjE8TORieyGkSAOoLgXOSLVuybgFqZCzKqpmZ+vv7i3pow6fB40qASdLOnTvrKkBdsmSJzjjjDAYG1amJmjdMMAygHhE4Y1KKlnobqaQBfM453XPPPWppaYkNOEtNJnLkyBEtX7684sGB1TJt2jQNDAwUKiuQ1zz2JmzesEnzL52vN/72jVm3BABGhMAZk9LChQt16NChqmwrHOCee+652rt3b9EsgKXygpN6vfv6+jILnJ/61KdKEuWnamwi1R8mbxjAZEHgjElp2bJl6u3tHTbNcBrhgDYc2DY0NOi8885TW1tbYRKVpFq4wcQC0ZrKDQ0NkpLLyyWZO3fuiA4KmpubdemllxbV9kVt7L1rrz7/lM9n3YxE1Z6BDgAmAgLnJMuXZ90C1Eg+n9fevXtTrdvc3Ky2tja1trZqYGBATU1N2rJlSyGwjfYur1q1qmgQYKkANJhYYOvWrUWBbnt7u9rb27Vnz57CbF1z584tms47l8upr69vWBvMTHPmzNHq1avV3d2trq50ua9mVhQ0o7a6b+/O5HEJhgFgdAick9x4Y9YtQI10d3enToM4duzYQrUN7AAAEa9JREFUsKAyCIzDQXS46kTaWri5XE6rV69WPp8vbCPIKQ7P1iVJhw4dKqxz4YUXSlLJNgSPv3//fp08eTKxDWamNWvWEDSPofb17dXbGHnDADBmCJyBMoaGhoblKYcD4zS9y0mSpkEOb6/UOkltWLx4cclJXcxMS5YsYfBfBhZdtkhvuPMNiTnO5A0DQP2xLAYfjcTatWvdxo0bx/ZBr7vO+0vP84TT2dmZunazmemKK64Yt8Hlvffeq3379mnBggVasGBBoUY0ATMAAMOZ2Sbn3Nq4ZfQ4J9m5M+sWoEbiKlo0NjZqaGioKIVjIqQxXHTRRbrooosK18fzcwEAIEsEzpiU4ipdBEGzmWn58uVqamqiHBsAAChoyLoBQBZyuZwuueTx3NFozeSmpiatWrWKoBkAABTQ44xJa9myZWWrY1Qqn88zcQgAABMUgXOSjo6sW4Aaq2Z1jHw+rw0bNmhwcFCNjY1at24dwTMAABNIzQJnM3u/pDdK6vVverdz7of+sndJ+jNJg5L+0jn3k1q1Y1RuuCHrFmAMpa29XEpvb29hFsC4EnYAAGB8q3WP87845z4RvsHMLpD0ckmrJc2XdKuZLXfODcZtABgvgim0R5vuAQAA6lMWqRrPl/R159wpSbvNbJekJ0m6K4O2JHvVq7y/N92UbTswLpSayAQAAEwMtQ6c32Zmr5G0UdI7nHN9khZI+k1onR7/tmHM7DpJ10nSueeeW+OmxujpGfvHxLg22nQPAABQv0ZVjs7MbjWzB2Iuz5f0aUnLJHVIOiDpk8HdYjYVO32hc+5G59xa59xaTnsDAAAgS6PqcXbOXZVmPTP7nKQf+Fd7JC0KLV4oaf9o2gEAAADUWs0mQDGzeaGrL5T0gP//9yS93MymmdkSSedL+l2t2gEAAABUQy1znP/RzDrkpWF0S3qTJDnntprZzZJ+L+m0pL+o24oal12WdQsAAABQJyw8zXA9W7t2rdu4cWPWzQAAAMAEZmabnHNr45bVLFUDAAAAmEgInJO8+MXeBQAAAJNeFhOgjB+HD2fdAgAAANQJepwBAACAFAicAQAAgBQInAEAAIAUyHFO8vSnZ90CAAAA1AkC5yTvfW/WLQAAAECdIFUDAAAASIHAOcnVV3sXAAAATHqkaiQ5cSLrFgAAAKBO0OMMAAAApEDgDAAAAKRA4AwAAACkQI5zkuc9L+sWAAAAoE4QOCe5/vqsWwAAAIA6QaoGAAAAkAKBc5L1670LAAAAJj0CZwAAACAFAmcAAAAgBQJnAAAAIAUCZwAAACAFytElednLsm4BAAAA6gSBc5K3vjXrFgAAAKBOkKqR5Phx7wIAAIBJjx7nJM95jvf39tszbQYAAACyR48zAAAAkAKBMwAAAJACgTMAAACQAoEzAAAAkAKDA5O87nVZtwAAAAB1gsA5CYEzAAAAfKRqJMnnvQsAAAAmPXqck7zkJd5f6jgDAABMevQ4AwAAACkQOAMAAAApEDgDAAAAKRA4AwAAACkwODDJW96SdQsAAABQJwick1x7bdYtAAAAQJ0gVSPJ3r3eBQAAAJMePc5JXv1q7y91nAEAACY9epwBAACAFAicAQAAgBQInAEAAIAUCJwBAACAFBgcmOQd78i6BQAAAKgTo+pxNrOXmtlWMxsys7WRZe8ys11mtsPMnhW6/dn+bbvM7J2jefyau+Ya7wIAAIBJb7SpGg9IepGkX4ZvNLMLJL1c0mpJz5b0KTNrNLNGSf8h6WpJF0h6hb9ufdqxw7sAAABg0htVqoZzbpskmVl00fMlfd05d0rSbjPbJelJ/rJdzrku/35f99f9/WjaUTNvepP3lzrOAAAAk16tBgcukBSecq/Hv63U7bHM7Doz22hmG3t7e2vSUAAAACCNsj3OZnarpHNiFv2dc+67pe4Wc5tTfKDuSj22c+5GSTdK0tq1a0uuBwAAANRa2cDZOXfVCLbbI2lR6PpCSfv9/0vdDgAAANStWqVqfE/Sy81smpktkXS+pN9JulvS+Wa2xMya5A0g/F6N2gAAAABUzagGB5rZCyX9u6Q2Sf9jZlucc89yzm01s5vlDfo7LekvnHOD/n3eJuknkholfd45t3VUz6CW3vOerFsAAACAOmHOjY/U4bVr17qNGzdm3QwAAABMYGa2yTm3Nm4ZU24n2bLFuwAAAGDSY8rtJG9/u/eXOs4AAACTHj3OAAAAQAoEzgAAAEAKBM4AAABACgTOAAAAQAoMDkzy0Y9m3QIAAADUCQLnJE95StYtAAAAQJ0gVSPJnXd6FwAAAEx69Dgnefe7vb/UcQYAAJj06HEGAAAAUiBwBgAAAFIgcAYAAABSIHAGAAAAUmBwYJIbbsi6BQAAAKgTBM5JOjqybgEAAADqBKkaSW691bsAAABg0qPHOcmHP+z9veqqbNsBAACAzNHjDAAAAKRA4AwAAACkQOAMAAAApEDgDAAAAKTA4MAkn/1s1i0AAABAnSBwTrJiRdYtQB3J5/Pq7e1VW1ubcrkcbQEAYJIhcE7y/e97f6+5Jtt2IHP5fF4bNmzQ4OCgGhsbtW7duswC1npqCwAAkwk5zkk++Unvgkmvt7dXg4ODkqShoSH19vbW5HHy+by2bdumfD6feVsAAEAxepyBFJqammRmcs6poaFBbW1tVX+MtD3JbW1tamxs1NDQUM3aAgAAhiNwBsrI5/PasmWLnHMyM3V0dNQkNSKuJznucXK5nNatW0eOMwAAY4zAGSgjHNBK0sDAQE0ep5Ke5FwuR8AMAMAYI3AGyhir1Ihq9iRTdQMAgOojcE7y5S9n3QJkIBp0jiSgHWngWo2eZKpuAABQGwTOSRYtyroFGGOlgs5KAtp8Pq/bb7+90EO9fv36MQ1c0+ZKAwCAylCOLsk3vuFdMGlUo9Rbd3e3hoaGCtvo7u6uZhPLClJLzIyqGwAAVBE9zkk+/Wnv77XXZtsOjJlq5DM/+uijiddrjaobAADUBoEzEFKNoDNcgSPu+lhIk1rCAEIAACpD4AxEjHaA3tKlS7Vp06ai66NVLsitNAhmACEAAJUjcAZ81eqBXbZsmSSpp6dHCxcuLFwfTbuSgtyRBMEMIAQAoHIEzoCqXwlj2bJlow6YA+WC3JEEwUzbDQBA5Qick9xyS9YtwBiJq4RRLz2w5YLckQTBDCAEAKByBM5JCCZQB8oFuSMNgpm2GwCAyhA4J/nCF7y/r3tdlq3AGGhvby/0Ojc0NKi9vT3rJhUpF+QSBAMAUHsEzkkInCeNXC6n9evXk7oAAABKInAGfPTaAgCAJEy5DQAAAKRA4AwAAACkQOAMAAAApECOc5If/jDrFgAAAKBOjKrH2cxeamZbzWzIzNaGbm83sxNmtsW/fCa07BIzu9/MdpnZv5mZjaYNNTVzpncBAADApDfaVI0HJL1I0i9jlnU65zr8y5tDt39a0nWSzvcvzx5lG2rnU5/yLgAAAJj0RhU4O+e2Oed2pF3fzOZJOtM5d5dzzkn6kqQXjKYNNXXzzd4FAAAAk14tBwcuMbPNZrbBzP7Yv22BpJ7QOj3+bbHM7Doz22hmG3t7e2vYVAAAACBZ2cGBZnarpHNiFv2dc+67Je52QNK5zrnDZnaJpO+Y2WpJcfnMrtRjO+dulHSjJK39/+3db4xcVR3G8e+TpbVQqAValbBEbFC0MU1pVlLFAAIqrQYkvhATAy80vjFGYoipMdFowgt4QYzGGA3+jYqJ+J9AbIOgibZIC20pXRRKqjQUFjW1wItC8eeLc5q7Njtl/pw7Z+/yfJKbvb0zc3p+z9ydOXvn3LlTUz3vZ2ZmZmbWtlccOEfEFYM2GhFHgCN5fYekfcBbSEeYJ2fddRJ4atD2zczMzMzGrZWpGpJWSprI66tIJwE+EREHgeckrc/fpnEd0OuotZmZmZnZvKF0jt6QD5auAb4OrAQOATsj4v2SPgx8BTgKvAx8KSJ+mx8zBXwfOBm4G/h09NEJSc8Cfx+6s/PLCuCftTsxTziLhrNoOIuGs2g4i8Q5NJxFw1k0Rs3ijRGxcq4bRho423AkbY+IqVe+58LnLBrOouEsGs6i4SwS59BwFg1n0WgzC19y28zMzMysDx44m5mZmZn1wQPnOr5duwPziLNoOIuGs2g4i4azSJxDw1k0nEWjtSw8x9nMzMzMrA8+4mxmZmZm1gcPnM3MzMzM+uCBcwGSzpF0r6RpSY9I+kzefoakLZIeyz9Pz9sl6WuSHpe0W9K6WW3dktuYzveZ6zLl81bhLG6WtCcvH6lV07CGyOKtkrZKOiLpxuPaulLSX3NOm2rUM4rCWXxX0oykPTVqGVWpLHq10yUFs1gi6S+SduV2vlyrpmGV/B3Jt09IekjSneOuZRSFXyv2S3pY0k5J22vUM4rCWSyXdIekR3N776xR07AKvlacn/eHY8thSTcM1JmI8DLiApwFrMvrpwF/A1YDtwCb8vZNwM15fSPp4i8C1gP35+3vAv4ETORlK3Bp7foqZfEBYAvpsvBLge3Astr1tZzF64B3ADcBN85qZwLYB6wCFgO7gNW166uRRb7tYmAdsKd2XZX3iznbqV1fpSwEnJrXFwH3A+tr11cji1ntfRb4CXBn7dpq5QDsB1bUrmmeZPED4BN5fTGwvHZ9tbKY1eYE8DTpYid998VHnAuIiIMR8WBefw6YBs4GribtrOSfH8rrVwM/jGQbsFzSWUAAS0g79WtIbwDPjK2QAgpmsRr4Q0QcjYgXSIPFK8dYysgGzSIiZiLiAeCl45q6EHg8Ip6IiBeBn+Y2OqNgFkTEH4F/j6PfbSiVxQna6YyCWUREPJ//uSgvnTrzveTviKRJ0sGH28bQ9aJK5tB1pbKQtIx0wOE7+X4vRsShsRRRSEv7xeXAvogY6KrUHjgXJulc4ALSEY/XR8RBSE866S8gSE/2k7MedgA4OyK2AvcCB/Pyu4iYHk/PyxslC9JAeYOkUyStAN4DnDOenpfXZxa99Mqok0bMYkEplcVx7XTSqFnkqQk7gRlgS0S8arMAvgp8DvhvS10ciwI5BLBZ0g5Jn2yrn+MwYhargGeB7+XpO7dJWtpid1tV8D3kWuD2Qf9/D5wLknQq8HPghog4fKK7zrEtJJ0HvA2YJA2MLpN0cfmetm/ULCJiM3AX8GfSjr0VOFq8o2MwQBY9m5hjW6eOph1TIIsFo1QWCyHTEjVExMsRsZb0+nmhpLeX7OO4jJqFpA8CMxGxo3jnxqjQfn1RRKwDNgCfehW8n/ZyEml62zcj4gLgBdK0hs4p+Lq5GLgK+Nmgj/XAuRBJi0hP5o8j4hd58zN52gH550zefoD/P3o6CTwFXANsi4jn88eOd5Pm/XZKoSyIiJsiYm1EvJc0eHxsHP0vacAseumZUZcUymJBKJVFj3Y6pfR+kT+Cvo+OTe2CYllcBFwlaT9pWtdlkn7UUpdbUWqfiIhj7yUzwC9J0946peB7yIFZn8LcQRpId0rh14oNwIMRMfB0WA+cC5Ak0tyh6Yi4ddZNvwGuz+vXA7+etf06JeuB/+SPGP4BXCLppLyDXEKax9MZpbLIH7uemdtcA6wBNo+liEKGyKKXB4A3S3pT/iv52txGZxTMovNKZXGCdjqjYBYrJS3P6ycDVwCPlu9xe0plERGfj4jJiDiX9Frx+4j4WAtdbkXBfWKppNOOrQPvAzr1TTwF94mngSclnZ83XQ7sLdzdVrXwHvJRhpimAfhbNUoswLtJH53vBnbmZSNwJnAP6UjpPcAZ+f4CvkH6poSHgalozvD8FmmwvBe4tXZtFbNYkjPYC2wD1taubQxZvIF0ZOAwcCivL8u3bSSdRbwP+ELt2ipncTvpHICX8vaP166vRha92qldX6Us1gAP5Xb2AF+sXVutLI5r81K6960apfaJVaRzZXYBj/h1k7Wkb6faDfwKOL12fRWzOAX4F/DaYfriS26bmZmZmfXBUzXMzMzMzPrggbOZmZmZWR88cDYzMzMz64MHzmZmZmZmffDA2czMzMysDx44m5mZmZn1wQNnMzMzM7M+/A+quiVHS+7iyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sonel_trend_array = np.full(np.shape(filenames)[0], np.nan)\n",
    "sonel_trend_array2 = np.full(np.shape(filenames)[0], np.nan)\n",
    "site_name_array = []\n",
    "\n",
    "\n",
    "for f, filepath in enumerate(filenames):\n",
    "    print(f'\\n\\n{filepath}')\n",
    "    print(f'f: {f}')\n",
    "    ## ToDo - extract and keep other important thing about each site from the file...Lat/Lon, name, etc.\n",
    "    \n",
    "    # Read in data\n",
    "    df = read_GPS_SONEL(filepath)\n",
    "    \n",
    "    # Get trend using linear regression\n",
    "    res2 = test_calc_OLS(df, 'Vertical', '2009-6-21', '2010-8-4')\n",
    "    sonel_trend_array2[f] = convert_trend_toyearly(df, res2)\n",
    "    \n",
    "    res = test_calc_OLS(df, 'Vertical', '2011-11-17','2016-10-15')\n",
    "    sonel_trend_array[f] = convert_trend_toyearly(df, res)\n",
    "    \n",
    "    # Save site name\n",
    "    site_name_array.append(filepath[-24:-20])\n",
    "    \n",
    "    print(f\"Trend 1 = {sonel_trend_array[f] * 1:.2f} mm/yr\")\n",
    "    print(res.summary())\n",
    "    print(f\"Trend 2 = {sonel_trend_array2[f] * 1:.2f} mm/yr\")\n",
    "    print(res2.summary())\n",
    "   # np.savetxt('GPS_trend.csv', {sonel_trend_array[f] * 365:.2f}, delimiter=',')\n",
    "   \n",
    "    # Make Plot\n",
    "    test_plot_OLS(df, res, res2, filepath[-24:-20], 'Vertical', sonel_trend_array[f], sonel_trend_array2[f], '2009-6-21', '2010-8-4', '2011-11-17','2016-10-15')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}