1_DFT_and_uncertainty_propagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1 Discrete Fourier Transform and uncertainty propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Before usage of this notebook, please download folder `Files_for_notebooks.zip` from this link https://github.com/harislulic/ZeMA-machine-learning-tutorials/releases/tag/v0.1.2\n",
    "and store the files in the same folder which is location for this notebook.\n",
    "It is also necessary to install pip in your environment and using this, package PyDynamic:*\n",
    "\n",
    "pip install PyDynamic\n",
    "\n",
    "*In order to see interactive diagrams, write:* \n",
    "\n",
    "pip install ipywidgets\n",
    "\n",
    "jupyter nbextension enable --py widgetsnbextension   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Introduction "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Signals are represented mathematically as functions of one or more independent variables. The independent variable in the mathematical representation of a signal may be either continuous or discrete. Continuous time signals are defined along a continuum of times and thus are represented by a continuous independent variable. Continuous time signals are often referred to as analog signals. Discrete time signals are defined at discrete times and, the independent variable has discrete values; i.e. discrete time signals are represented as sequences of numbers. In practical setting, such sequences can often arise from periodic sampling of an analog signal. In this case, the numeric value of the *n*th number in the sequence is equal to the value of the analog signal $x_{a}(t)$, at time *nT*; i.e.:\n",
    "$$x[n] =x_{a}(nT), {-\\infty}<n<{\\infty} $$\n",
    "\n",
    "The quantity *T* is called the *sampling period* and its reciprocal is the *sampling frequency*.[8]\n",
    "\n",
    "A discrete - time system is defined mathematically as transformation or operator that maps an input sequence with values *x[n]* into an output sequence with values *y[n]*. This can be denoted as:\n",
    "\n",
    "$$y[n] =T[x[n]]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with continuous time signals, discrete time signals may be represented in a number of different ways. For example, sinusoidal and complex exponential sequences play a particulary important role in representing discrete-time signals. This is because complex exponential sequences are eigenfunctions of linear time-invariant systems and the response to a sinusoidal input is sinusoidal with the same frequency as the input and with amplitude and phase determined by system.\n",
    "\n",
    "Linear time invariant systems (more about in [8]) can be described by their impulse response, and the response to any time-distributed signal is described by convolving the signal with the impulse response. Convolution is a commutative and associative operation, but it is also quite computationally and conceptually expensive. For discrete signals, with the development of efficient Fast Fourier transform (FFT) algorithms, almost always, it is faster to implement a convolution operation in the frequency domain than in the time domain.\n",
    "\n",
    "Properties of linear time invariant systems and their combinations are much better described and manipulated after Fourier transformation, for example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fourier transform of a signal tells you what frequencies are present in your signal and in what proportions.\n",
    "\n",
    "*Have you ever noticed that each of your phone's number buttons sounds different when you press during a call and that it sounds the same for every phone model? That's because they're each composed of two different sinusoids which can be used to uniquely identify the button. When you use your phone to enter combinations to navigate a menu, the way that the other party knows what keys you pressed is by doing a Fourier transform of the input and looking at the frequencies present.*[9]\n",
    "\n",
    "The magnitude square of the Fourier transform, $|X(f)|^2$ instantly tells us how much power the signal x(t) has at a particular frequency f."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1.1 Discrete Fourier Transform (DFT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Discrete Fourier Transform (DFT) is the equivalent of the continuous Fourier Transform for signals known only at instants separated by sample times (i.e.a ﬁnite sequence of data). The DFT corresponds to samples, equally spaced in frequency of the Fourier transform of the signal. In addition to its theoretical importance as a Fourier representation of sequnces, the DFT plays cetral role in the implementation of a variety of digital signal processing algorithms. This is because efficient algorithms exist for the computation of the DFT. These efficient algorithms are called *Fast Fourier transform (FFT)   algorithms* [8]. In particular for dynamic metrology, where the involved quantities have time varying values, the DFT is considered an essential tool. It is employed, for instance, to assess the frequency content of a measured time series, design suitable digital filters or carry out input estimation.[7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time taken to evaluate a DFT on a digital computer depends principally on the number of multiplications involved, since these are the slowest operations. With the DFT, this number is directly related to $N^2$ (matrix multiplication of a vector), where *N* is the length of the transform. For most problems, *N* is chosen to be at least 256 in order to get a reasonable approximation for the spectrum of the sequence under consideration – hence computational speed becomes a major consideration. Fast Fourier Transform (FFT) algorithms rely on the fact that the standard DFT involves a lot of redundant calculations and reduces their number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2. The structure of tutorials"
   ]
  },
  {
   "attachments": {
    "sensors.PNG": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGzCAYAAAAmH71NAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAFiUAABYlAUlSJPAAAGdmSURBVHhe7d0FeBRn/gfw+99dvb3rtVfvXf3qQmmLO8UDBCvuVrxQoIIE91Lc3Ys7hQSH0uIUd+IhxN32+5/fm1kIYXczkwDZTL6f5/k9yc7Mzs5slp0v7/vOzF9AREREZFEMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6REREZFkMOkRERGRZDDpERERkWQw6RERE2WSz2RAfH4/U1FR9CrkbBh0iIqJsSktLQ58+fXD8+HF9CrkbBh0iIqIcSE5OVoGH3BODDhEREVkWgw4RERFZFoMOERERWRaDDhEREVkWgw4REd1BTpteuXI5mjZugK/btkSHdubra62aN2uEdm1aOJyfdbVC2zbNtXU01tbVysF8I9UKLVs0QeuWTXOwjpZor70HzZs2VD8dzc+6WqF1q2Zo0byx+t3xMq5Ltr+F9l60bdVce5y9ddj3I7t/Uyn5e7Zs3gS7d+3QPy3ujUGHiIjuIGcR/fjDd2jTqDw2LR2EDYu8TNXGxQOwel5fFPzodUwd1Ul7PNDhcq5q89LBGNG3JYoXeg9rF/RT63S0nKva/Mtg1PUohnbNKmPzMvP7ISWvO3f8N/jo3VewcvaP2dqOjUsGonv7mqj65ef4dfkQh8u4KnnN9Qu9UKrohxj4XeNsr2POuG7p+zEne/uxYdEA9beoX7M4pk2ZrH9a3BuDDhER3UFadAYO6IcZP3UFYncDEdvNVeQOpIVsQZ3qJXDp8Hwgaqfj5VyULW43/tg6AS0aVATCtqp1OlrOaUVuR1rCXgz6rhlmjeuevf2Q0l436NQyVK9UGEmBm7R9MbkdWtmid+GXWT+iZ6d6QOJetW2OlnNa8n6G+6Bt0yrYsXYMEK+tw9Fyrkr24+RStR+pQZvVY4fLuSrZjrBtGNqnJWbMmK5/Wtwbgw4REd3BHnSkNQYRPrBd/9VUIXQrEvw3oFa1Yjh3YBZwY5vD5VyW9rr7No5F03pfIjlICxjaOh0u56LStAOzV8/GmD5GC2zZ2A8peV3fY4tQrcIXiL22LlvbYQvzxuJp36F7hzpAzA7Hy7goeU0Jjq0aVYL3yhEqbDlazlXJOq4dXajtRyHE+a7P1n7Ic1KDN2Pg980xY/o0/dPi3hh0iIjoDgw6t0pel0EnveQ5DDpERJTnMejcKnldBp30kucw6BARUZ7HoHOr5HUZdNJLnsOgQ0REeR6Dzq2S12XQSS95DoMOERHleQw6t0pel0EnveQ5DDpERJTnMejcKnldBp30kucw6BARUZ7HoHOr5HWdBZ20EG2+tt4s949BJ9cw6BAR0R0YdG6VvK7TFp3I7Vg13wt+xxcD4d63Pe+2YtDJNQw6RER0BwadWyWv6zDoaMEjTVvn2mVD4XdyqQozmZ97sxh0cg2DDhER3YFB51bJ6zpt0dHCB+J2pbfmhNz+vNuKQSfXMOgQEdEdGHRulbwuByOnlzyHQYeIiPI8dws6KcHZCDqhWtDRAoEKOj91y1HQ8Tuew6ATLkHn+1tBR9s2h8s5KXnNtOtb0LqxFnRWZT/o+B5LDzrxDDpERJSf2YPO9NFdgOgdQJi3uQr3QUrQJtSuVhwX/pirQgbCtpkqm/a6v20Zh+b1v9QO8ptV95AEJsOlrSMtdhcG9GqSHnS0cOBwuaxK2x//E4tvBQTZP0fLuSgZtLxk+g/oIUFHurpkHx0s57RkX7RwJEHHZ/VILSxp63C0nKvStjs9sBVCot+GbO2HPEcC1+AfWmAm715ORER5lQSdAV790LNDLfyxdSJ+2/yzudICyq61Y1Ci8EeYOrIzNizywvoF/U3VpiWDMLJ/KxT86HVMHNFBCwrfYdGUnoZr8bTemDfpW5Qv+TGqffk5pozoiGmjO6tWKjM1fUwXfNe5Ll7777MY8n0zLTR1cbics5LXnDyyI2p7FEXhgm9jkrYdsk5HyzqrGdryowe0wUfvvYpWDb/UtqGr6X2R1+yt7cfrrzyHoT+a3w+paWM6Y8LQ9vjKszx+WbZU/7S4NwYdIiK6g2rR6d8Hn374OprWK48mJkue07BOWbz8wjMoXLgIypcrh3Jly5qq8uW0n2XKoHjxYihSpEi2qmjRovji88/x6acFUfCzz1Gw4GfZrk8//VRbh+N5WZZ67YL6OrK5Hdrz1PO19eRkX3K6HwW097JAgU+xatUK/dPi3hh0iIjoDvauqykjO6ePbZExJSYKN+xjdErgxPGjSEpMQEJ8rMmKQ1xsDCLCbyA8mxUREYagwED4+fkhKCgIgdrv2a2rV686nG6k5LX9/f3h6+ub7e2Q5127dk2tJyf7ktP9kPdyzJgxWtBZpX9a3BuDDhER3cEedKR7BJHb1SBUU3VjG5IC0oPOxQtn9bWSFaSlpWHgwIGYNo2DkYmIKI+yBx0Zl5Gzs65K4NzZ0/payQpSU1MxaNAgzJgxQ5/i3hh0iIjoDgw65AyDDhER5XkMOuQMgw4REeV5DDrkDIMOERHleQw65AyDDhER5XkMOuQMgw4REeV5DDrkDIMOERHleQw65AyDDhER5XkMOuQMgw4REeV5DDrkDIMOERHleQw65AyDDhGRm7KlpcCWFKs/IlduBZ3s3+sqMWADansw6FgNgw4RkZtJi72BmAOTcWNZA6RGB+pTyZXbgk6UFnRuSHgxUWHbkBywSQs6pXDu3Bl9rWQFDDpERG4gLS4UiRe8EbG+MwJHvwn/gU8iZEZpOYLrS5ArEnQGePVFlXKfYWifFhj0fXNTNfiH5vDq1QTvv/M6Llxg0LESBh0iolxlQ/yZdVqoKQX/QU/Bz+txFXL8BjyBSJ8B+jKUFXvQqVe9BKb91A1TRnUxVVO1mjCsIwp89CYunGfQsRIGHSKiXGVD4uVdCJ70uQo3Enb8B/1LVcJFH30Zyoq962r2uG+ApH1AzC5zFbcbaWHeqFezNMfoWAyDDhFRbrOlIWJzT/j1fQj+A/6hBZ5/qOCTGh2gL0BZuTVGJydnXa3nWVcWxKBDRJSbtAN0zMGZCBj6nAo5gWPfge8Pf0XYytawJcfrC1FW7k7Q4enlVsSgQ0SUa2xIOLsBASP+C3+vJxC6qC4Sr+7DjaUNEHtotr4MGcGgQ84w6BAR5Qot5JzbiKCxb8OvzwMImVUeyWGX0uckxcKWGKV+J2MYdMgZBh0iolyQ5Ps7gsZ9CF8t5ARPLoTk4JP6HMoOBh1yhkGHiOg+Swk5jZCpxeDX7xEETyqIRN/9+hzKLgYdcoZBh4joPpKQc31GGfj2+TsCx32AhPPaQZZyjEGHnGHQISK6T1Ij/RG6oDp8+z6IgGEvIfboQn0O5RSDDjnDoENEdB+kxUcgfGVL+Hs9ioDhWsg5Ml+fQ3cDgw45w6BDRHSP2VISEL6+E/z6PwL/gf9E9G+T9Dl0tzDokDMMOkRE95ItFVE+A+Hf/3H4e/0DEVt/kOSjz6S7hUGHnGHQISK6V9KSEbNvnLrqsZ+EnM09kZYUrc+ku4lBh5xh0CEiuidsiPljOgKHvahOIw9b2xFpvAjgPWMPOtNGdwaidgA3tpmrMG8kBW5EbQ8GHath0CEiuutsiD+5AgHDX9ZCzsMInV8VqbHX9Xl0L9iDzpiBrRF1ZQ3CLqw0VeEXVyHo5DJ4VCqGc+fO6GslK2DQISK6q+TWDpsROPpNdTfy6/M9kBJxTZ9H94oEnQFeffD+2/9BHY8SqmXGTMlzPKsUxcsvPotLF9miYyUMOkREd1Hi5Z0IGveRCjkhM8sgOfSsPofuJdWi49UPg79vgmvHF+Hy4Xmm6sqR+Ti7fxYqli3EriuLYdAhIrpLUsMvI3jSF/Dr/xiCxr6HpIDD+hy61+xdVzPGdAVid6oByaYqcjtSgzejtkdJBh2LYdAhIroLUsOvIGRGafj1eRBBP7+PxKv79Dl0P9iDDs+6oswYdIiIcig1JhihC2rAt99DCBz5KuLPrNfn0P3CoEPOMOgQEeWALSEKYSuap1/1eOhziDmYN75MrYZBh5xh0CEiS4iKikJQUJD6UrtfbEmxiFjfGX5ej6mLAkb/PkWmps+k+4pBh5xh0CGiPE0OcFu2bMEHH3yA+vXrIzw8XJ9zb9kSYxC5oZsWch5XVz2O2j5Ym5imz6X7jUGHnGHQIaI87eTJk3jrrbfQqVMnhISEqAPePacFmpg9o+E/8F/w93oC4Ru786rHuYxBh5xh0CGiPEsObgMGDECJEiWQkpJ+o0yZdubMGRw4cOAete7YEPP7NPgPeU6dRh6+rhNsyTH6PMotDDrkDIMOEeVJsbGxiI+PR79+/fDKK69g//792Lt3Lxo0aIAnn3wSTz31lApA/v7++jPuhjTEHpqtbu3g2/ch3PilCdLiQvV5lJsYdMgZBh0iynMCAwPRrFkzdOnSBQcPHkSBAgXwl7/8BY899hgqVaqEefPmYenSpSoASRCSg2DO2RB/dgMCR7+hzrC6PrcKUsIu6fMotzHokDMMOkSUp8gBberUqXj00UexZMkSJCcn4/Lly1i/fj127dqF6Ojom8t5eHigbdu26nFOJV7ZhaCx76qWnOCpRZF846I+h9wBgw45w6BDRHmGhBpvb2+0bNkSNWvW1KfeSQ56y5cvx/PPP49NmzbpU7MvOfAYgsZ9AD+5IODY97THR/U55C4YdMgZBh0iyjNkTM7HH3+MBx54QI2/kS+wzI4cOYIOHTrgueeew4gRIxwuY4bcv+r69FJayHlYhZ3489pBkdwOgw45w6BDRHmKnE5eunRpFXak1SYzOduqadOmaoxOTkNOSthldWsHv36PImDUq4g7s06fQ+7mVtDprIIOJLyYKS3oJPpvZNCxIAYdIsoTpNvKTs6kKlWqFJ544gnMnz9fn3qL/VTznLAlROLG4rqqJcd/yLOIO7ZIn0PuyB50po3uAkTvAMK2matwbyQHbuTdyy2IQYeI3Jp0RfXs2RN169bFyJEjceHCBTU9ODhYjdP561//iqFDh6oD3d2SlhiD8LUd0kPOsOcR/dt4fQ65KxV0vPqhdrVimDyyEyYM72iqJg7vhLGD26PAh2/hwvkz+lrJChh0iMgtyYFr7dq1ePnll9Xp44UKFVKnkL/33nvYuHGjWkbubyXdVBJ2Bg8erKblVFpCJCI2favuXyUXBYzaNVymps8kt+Tn54f58+eiz4/fwbNyEYwe0AYj+rc2VSP7t8HQH1viw/dew4VzDDpWwqBDRG7pxIkTeOONN9CrVy/VeiO3d5CA89JLL6mLAcpYHCGnk7du3VpdVyenbGnJiPQZAP+BT8Kv/6MIX9cZtpQEfS65q5iYGPz555/w6t8Hc8Z1BxL2pHdfmamYnUgL3Ya6NUqx68piGHSIyO1Ia47c2uF///ufarXJ6Pfff8eLL76oLgwoX2BCls/5uBwbYg/OQMCgp+Df/3GErWwNW2qiPo/cneq6yvFZV+s5GNmCGHSIyO3IQat9+/aoWLEiEhNvDxsyb/To0aq7Si4UeHdoIefYQgQMfhZ+/R5B6KI6SI27oc+jvODuBB2eXm5FDDpE5HbkoDV27Fh1VpUMRs5I5q1cuRIPP/wwTp++GwckG+LPbUbgqNfUVY+vz6mIlPAr+jzKKxh0yBkGHSJyCzLWRsbgnDt3Tj0+e/as6rp69dVX1Xgdu7S0NHWPq3feeUddQDCnkq7uRfCEAvDr+yBCZpRG8nUORM2LJOgMGzIQs37upo/R2em4YhxMU9N3IU0LO3Wrl8K1q7yHmZUw6BBRrgsNDVX3pXrrrbewbNkyNU0OXPK7DDyWM6+mTJmibv/QtWtXddXjDRs2qOVyIsn/DwRPLADffg8jeEphJLnZrR1SU1KwccM6zJ41A3PnzMpWzZs7CwvnOZ5nvGZjwbzZmD93toN5xmqe9tz07cj+OuZrz5e6Y562j4sXzUeVypXRwLMUFk7ujbkTvr2t5qifPTF/Ui/Mm9jzjvkyb8ZP3+C9t1/F0CGD1Pbe8TpayfQFOdwPKVmH/G0czTNSc7TPxJrVK5CUyMHyWWHQIaJcJYFGwoucQi6tOJnt2LFDXRxQTi1/6KGH1EBkR1dENisl5BSuzygJv75y/6q3kXDRR5/jPqKjI1C2dFG0aVIBg79riv49G5mqAb0bo1XDitrB+7/w6tkEXr0aO1zOVQ3QnvN916/weYH/oVGdMhj8QzOHy7mqQdq2t2lcCZ988Dp6dqqLgb2bOFzOVQ3s3RRVyn2OMsU/1NYn+3L7fK+ejfFDtwb4sXtDbZ42X/Y3Qw3Qpn3ftb62H2+jcd3y2n40v21+f71kHf2+vX3d9pL9aN24otqPXtncD/kbyHv66cdvaqGstPY3Mr8O2Y5OLauieOFPERIcqH9ayBkGHSLKVTdu3MCbb76JyZMnq8cSfORU4TFjxmDmzJlqMLKcefXHH39g69atCAgIUMvlRGqEL0Lne2gh5+8IHP0m4s/mvHXoXoiICEeDetXw5+4pQOR2h2NLXJVc8fe3zT+jUe2y2hu9VY1DcbScq5J1RF9dg65tamLL8iGqi8fRcq5Ktv3gtglaSKiM62e1kBpufgyNLcIHU0Z1UQEFEdsd7otMc7aPsh9RV9aic+ua2LpimOrCcrici/dI9uOPrePRtkllhJ6T/fB2uJyrkvWnaX+Lji09sGHxoGy9FzIG6cqhOWhYrwr8/Xz1Tws5w6BDRLkqNjYWBQsWRJEiRbBgwQJ1Q85//vOf6sKATz/9NGrXrn3b7R9yKjU6CDcW1oLciTxg+IuIPb5Un+N+JOg0+soDp/ZOA6J23DyQGy1bmDcO/DoejeuW04LOtvRysJzL0g7mMdfWoVs7z/SAELfb8XIuyqZt+yHvSWjbtApunF9xM6iYKZsWMqaN6YahfVqkvxdm90X24+o6Fdi8Vw7P9n4c3DYR7Ztp+3FB9sPH4XIuS9vutPBt6NSqOjYt0YKjtl/2EGq4tOdcOTxP+2xUY9AxgEGHiHKVtOBs2rQJX3zxBf7973+jbNmymDNnDq5du4Z58+apsLNz50596ZyxJccjbFUb+A14XN2/Kvq3CdpEfaYbsgedk7unZqtFx6YdVH/bPC496OgHSYfLuSjVonNlLbq29cSvy4dmq0XHpgUbCQhtm1TJUYvO1NFdMeRHLehIODC5L6pF5/JadNGCzjYXLTouS9uPP7ZOQDstsOW4RUcLOhsXD1ZhydFyrkqec/nQXO2zUZVBxwAGHSK67y5evKhuxildVHZBQUHqjKvw8HB9CrBv3z41Jmf//v36lOyzpSQicktv+Hk9Dr9BTyF6zxgtZLn3rR0YdDIUg87NYtAxh0GHiO4bOTV8/fr1eOGFF1RLzd/+9jd069YNCQm3nzkiVzmWa+QULVoUtWrV0qfmgBZoonYOg7/XP+DX/3FEeveFLQ/cv4pBJ0Mx6NwsBh1zGHSI6L7x9/fHBx98gBYtWmDz5s1o166dCjsNGjRQ97MS0pU1fPhwdf2cypUr4+rVq2p69tkQ+/t0BAx5Dn5eT6iuK1tKnD7PvTHoZCgGnZvFoGMOgw4R3Rdy48WDBw/io48+UmdaCfkCkjOrJOxUqVLlZtiRs6tmzZqlHegj1OPssyHu2CIEDHsJ/l6PIWxFC6TFherz3B+DToZi0LlZDDrmMOgQ0T0np4jXrVsXpUuXRrVq1W67Aae04Mh1ceQaOSVLlrxLt3UQNsSfXIHAUW/A98e/I3SBJ1Jj8tY1Rxh0MpSLoKP2LdL1Ohl08i8GHSK652RszoQJE/Dss8+q08ZPnTqlz7ll7dq1auxO06ZN1fI5lXh1N4J+fhd+fR5AyMwySL5+58UI3R2DToZyFnS0n3GBG3HlzyVICd6iHt/xXK0YdPKvpKQk9OjRAz/99JM+xb0x6BDlYTIuR86ikqsgX7p05/2E5KKAR4/m/DYMSX5ya4eC8O37MEJmlERyyJ3BKi9g0MlQToIOwrxx7fQyTBzWEfEBm9T23vFctRyDTn4lLcj9+vVTt5HJCxh0iPIIaZU5dOgQtmzZolpwpItKyDVznn/+eXVTzgMHDqhpd1NS0AmETC8J3z5/R9Dkz5Hk+5s+J+9h0MlQToKO7fpWJAZvQZjveqTK49AMz8lQDDr5l3RdDRkyRI0HzAsYdIjygMjISLRp0wbPPPMMXnrpJTX+pnv37jcHIcs4nHfffVfN27VLO3DqISinUsKv4PqcCuomnem3dtioz8mbGHQylNOgo22jlJOWHHsx6ORfHKNDRHeVNBN///33+Pjjj1VrjrTa9OrVS92Uc+zYsfpSwJkzZ9QZWMWLF1f3ssqptOhA3FhUC/79H4X/iJcQf3KlPifvYtDJUC6CjpFi0Mm/GHSI6K6SKxzLNXBWrVqlHktrjbTmFCtWTF0ROSO5W7m06OT0XlZpidG48Utj+PZ7CP6Dn0bMH1owsAAGnQzFoHOzGHTMYdAhorvq8uXLeOqpp3DkyBEVcuQLRs60sp82LtfK+fXXXxEfH68e55QtOQERG7rBr/+jCBj8b0TtHAJb2q3T1/MyCTqN63vgzP7pQOxO7cC6zVSlH5gnokm98urgqEo7OJuqqO2I89uAb76uDZ9VI4DEvdp6tOkmyha7G0e3T0X75tUQfkkLwFrIuLk9RkqCkbb/M8d+g+F9W2nvxS4t7DhYzlVF70jfj3a1sH3NSG0/9jlezkXZ4nbjyPYp+Frbj4jLq83vh5QW0tK091QC15Zf0oOjw+VclfZe+B5boH02eFNPIxh0iOiuCgkJUWNz5HTOESNG4K233sLu3bvVPAk+y5Ytw4cffnjz4oA5YUuKQ5R3f/hrAUeuehzpMwC2lLsToNxBpBZ0vqpTBbvXjUbYpdWqNcRMhZ5biS3LhsKzSlFcOTwfV4/Mx+XD80yVHFAPbBmPWtWKY+zg9rh4cA5O7p1uqs4dmIM543ugcrnPsH7hQBzdMRkHt04wXtsm4JD3RHRtWwMtGlRIf7xtouNlndSR7ZOwdfkwVCr7uRaWWuLPXdMcLueqju2Yiskju6By+c+wcdEgLbyZ3A+tZLv3b/oZHhULY8gPzXDYZ5LD5VzVUW1fVs7piyYNayM4KED/tJAzDDpElGPSiuPj46PuWSVfKl5eXvj73/+OBx98EL/9duusJ7mdg4zd+fbbb9VyOWND1O4RCBj0NPz6PYrwdZ2Rlhyrz7OG2OhwlChWGOVKfIyGtcugvmcpU9WwVml8WaoAnvn3v1Ci8CcoUcR8lSxaAAU/fhcvvvAs3nzjv/iswLso+Im5kue8/darWgB+Gv978xV8+P6b+OA94/Wh/Hz3Tbz04nN4/vl/4/1330ifZqLkNd/536t4+ql/4b8vP4+P3n/L4XKu6kPtOa+9+rL2fj6Vrf2Qku1+9+3X8az2Xvzn5RccLpNVyTpee+UlFC1aBOE3ruufFnKGQYeIsk1aaNatW6dOF69fvz78/PzU9OvXr6Nq1apqALLcz+rEiRPq9g+ffPIJypUrdxdac2yIOzIf/oOfhV//R9T4nLTESH2edURFhqNurSrYvGQgfI8vwsWDc03VlSMLtP/590GVyhVx+PAhHDl8GIcPHTJVR7TnHdDC6t69u7Fnzy7s2LEDO3earx07tqsbum7YuEnVxo2bDdcG+blJ+7lhI9Zp61C/Z1omq0p/zU03t0HW4Wg5l6Ved5P6zNt/d7ici5Lt3rRpi7YOfTsyzTdU2mvPX7AQP/74I0JD884tTXILgw4RZZucOfX222/jhx9+uOMUcRmDI3cmf+WVV1TgeeSRR1CzZk11Y8+c0ULO8SUIGPGf9JCztD5SY0P0edZiH4x8Zp+M0dkFuTiemVJjdH4dh6aN6+trJKuQ/1Q0bNjwLvx7sj4GHSLKFgk2ckl1ufCf/d5Vcqn1FStW4Oeff1ZXORbHjx/HypUr1anmcmPPnLEh/sx6BI5+C379HkLItBLq2jlWdVfOuto0Fo0b1tPWlvPbapD7kO7iRo0aMegYwKBDRNkiQWf69Ol44IEHVLCZM2eOOoVcBiL/85//VPetklPH76bkoOMIGvueGpMTPOkzJIfcrRuAuicGHXKGQcc4Bh0iyraIiAg0adJEXeH4zTffVE3p+/fvV9fSkTuRV6pUSV8y55KDTyB44qfw7fsggn7+AIm+d//2Ee6GQYecYdAxjkGHiHJEuq1koPG5c+f0KdohNS0NjRs3hqenpz4lZ1IjruL67Arw6/cIAke/gfizG/Q51sagQ84w6BjHoENEd5WEnHnz5qkuLBmXk1NpsaEIW9ZQhZyAES8j7sQv+hzrY9AhZxh0jGPQIaI7yOnh48ePR0CAuYuR/fnnn6hcubIanzNmzBh9avalxYcj7Jcm6qrH/kOeQczBvHH34buFQYecYdAxjkGHiG4jg4wl5Mgp4RUqVLh5bRwj5KrIcrHAPXv26FOyz5YQifANXeDn9Rj8B/4LUbtHaduWvw7WDDrkDIOOcQw6RHQb+VL4+uuvUaBAAXXtGxlULKeI31epKYjyGQB/ryfgP+CfiPTpD1tynD4z/2DQIWcYdIxj0CGi24SFhalwI6eOr1mzRt3G4Y033sDRo0f1Je4xWyqidg6D34An4d//cYRv/MZS968yg0GHnGHQMY5Bh4huI7dnkC4rudS/WLVqFV5++WUVdu72dXHuZEPMH9MRMOQ5+PV/DDd+aYq0+Ah9Xv7DoEPOMOgYx6BDRLeR08UPHTqkBiTbybVx5FYO//73v9W9guxOnTqFK1fu1pWJbYg/tRKBI19Nv7XDkvpIjQnS5+VPDDrkDIOOcQw6RPmc/fYNWfn999/x7rvvqqseb968GQcOHFAXCZw2bZq+RM4kXNqurpHj2+cBXJ9TEanR+TvkCAYdcoZBxzgGHaJ8bPfu3fjqq6/UzTmNkC/X8uXL469//SuefvppNG3aFNHR0frc7EsOOKKudizXygme/AVSwi7qc/K3iIgwNKlfHRcPztbepP1AzE5TZYvfjaM7JqFJo6+0tTHoWAmDjnEMOkT5lJxGXrduXXUaecGCBXH6tLH7RslFAB966CE1jufq1av61OyT+1eFTC2aHnKmFEai72/6HIqKCke1yuUwekBrrFzghRVz+pqsfhj0fRPUq1tHW9vtd5envI1BxzgGHaJ8KjY2FmXLlkWzZs3w3HPP4bXXXlMtPK4EBgaiWrVqKFeunBq0nFMp18/i+pwK8Ov7AAJ/+h8SLt/rwc55S4wWdMqXKYbGtcvguy5f4duOdUxVr051Ua9GCdTyrK6tjS06VsKgYxyDDlE+dfLkSTXG5sSJE/Dx8cGzzz6LJ5980mnYkRYgGY/z9ttv49KlS/rU7JNbO4QuqA6/fg+n39rh1Gp9DtnJGJ2G9arh6M5JSAnbhqSgTeYqZAv2bBijj9Fhi46VMOgYx6BDlE+dP39efVHab/Mgp45/8MEHKvCsXLlSTcssNDT0rnRXpcWEIGxFCy3kPISAUa8i9ugCfQ5lZB+MfGrvNCBqh/YH2GqqbGHeOLD5Zw5GtiAGHeMYdIjyKfnHL7d3SEpK0qdA3YH8k08+UWdWLVy4UE2T5eTigWfPnlWPc0quixO+vkt6S87wFxHz+xR9DmXGs67IGQYd4xh0iOg28gX65Zdf4uGHH1ZfDAMHDlRdWnfj/lVIS0GUd3/4D/wH/Af8A5Hb+gH57P5VZjDokDMMOsYx6BDRHeRigXLauZyRJfe7Wrx4sT4nB9KSEb17hLp3lZ/XPxCxuac2keNGXGHQIWcYdIxj0CGiO8jA49GjR+OBBx7AhAkT1BdFztgQ/dtELeTI/aseQ9jKlkhLyvn1d6yOQYecYdAxjkGHiG4jIWfTpk145plnMG7cOPU4Z2yIO74EAcNeVONyQhd6IjWaX85GMOiQMww6xjHoENFtJNgMGDAA7du3R1paTg+ONiScWY+gn96Bb5+/I3R+daREXtPnUVYYdMgZBh3jGHSI6A4JCQl3oSUHSLzog6Cx76qWnJBpJZEcek6fQ0Yw6JAzDDrGMegQ0T2RHHgMwZM+h1+fBxA0sQCS/A7qc8goBh1yhkHHOAYdIrrrkkPOIGRGaXVrh6DxHyHx6l59DpnBoEPOMOgYx6BDRHdVyo1LuD67AuQmnQEj/oP4c5v1OWQWgw45w6BjHIMOEd01aYnRuLGkHvz6P4aA4S8h9vBcfQ5lB4MOOcOgYxyDDhHdFbbEGISvbq0GHvsPfhrR+8frcyi77EHn1B7e64pux6BjHIMOEeWYLSVRXenYz+txddXjyG19ZaI+l7JLBZ36Hji7X/uCjtutTfAxVbbIHTi4dTwaN/pKWxuDjpUw6BjHoENEOZSGmH1j4T/oKXXl44iNPZCWGKXPo5yIiQ5HudJFUb9mSfToUAfd2nmaq/a14FmtKGrVrKGtLadXtyZ3wqBjHIMOEeWADTG/TUTA0OdVl1XYqja8tcNdFBUVjmpVymHs4LZYv2gQVs/zMlVrFw7AsB+bo17dWtra2KJjJQw6xjHoEFE22RB3dH76rR36PoDQJfWRFntdn0d3g73r6sz+6UDsLiDc21xFbccf0nXFMTqWw6BjHIMOEWWDDQlnNyBw1Bvp969a4IGU8Iv6PLpbeNYVOcOgYxyDDhGZpm7tIPev6vcQQmaWQQpv7XBPMOiQMww6xjHoEJEpib4HEDzxU/j2eQCB4z5Uj+neYNAhZxh0jGPQISLDUm5cwPVZ5eDX7yEEjddCDm/tcE8x6JAzDDrGMegQkSGp0UEInVsp/YKAw19EwkUffQ7dKww65AyDjnEMOkSUpbTEGNxY2kC15PgPex6xh+foc+heYtAhZxh0jGPQISKX0uIjEL66Hfy8HkXA8BcRc2AKYLPpc+leYtAhZxh0jGPQISKn5ArHkZt7wn/gk/Af9C9E7xsHWypv7XC/MOiQMww6xjHoEJETNkRtH6gFnCfh1/9xhG/oAltasj6P7gcGHXKGQcc4Bh0iupMWaGJ+n6pacfwG/APh6zrClhKvz6T7hUGHnGHQMY5Bh4huZ0tRg40Dhr4Iv/6PIWxlS6Ty1g65gkGHnGHQMY5Bh4gysN/a4XV1Gvn12RWREh2oz6P7jUGHnGHQMY5Bh4huSry6G4Gj5f5Vj6Tf2iGKX6K5SYJO4/rVcfY37Qs6brc2wcdU2aJ24NC2CWjc6CttbQw6VsKgYxyDDhEpSdf2IWjsO/Dt+xCCpxRCcshJfQ7llsjIcNStVRmblw3G1RNLcPnwfHN1ZCFWzumDhvUl6PCSAFbCoGMcgw4RaaHmtBZuiqibdAaOfRcJl3z0OZSb4uNjULpkMVQuVxCtGlZCs6++NFXN61dA+ZKfwLNGNW1tqekrJUtg0DGOQYcon0uN8EXovKrw7fswAn96CwkXtupzKDfZbDbs3rUT5csUgffKYfD9cwmuHFlgqq4dW4TVc6VFp66sMX3FZAkMOsYx6BDlY6lR/gid76HG5AQMexFxf67Q51Buk6AjX8wlinyMiwfnALG7gHBvU2WL2o6DW8dzMLIFMegYx6BDlE+lxYWpU8f9+mshZ+hziDkwSZvIg6E7iY+PR7NGnvhz1xSedUW3YdAxjkGHKB+yJccjYmN3+Hk9Bv/BT6ff2iGNt3ZwNzy9nJxh0DGOQYcon1EhZ3MvdTFAf6/HEekzkCHHTTHokDMMOsYx6BDlKzZE7x6pWnH8BzyJ8DXtkZYUrc8jdxMZGYFG9avj1J5pQPQOIHSrqbKFeePA5p/RRF1Hh6yEQcc4Bh2ifMOG2EOzEDD0eTX4+MYvjZGWEKnPu3+SkhIRHOQPf79rCPD3zXYFB15DYIDjeUZKnhscKL/73THPWPkhSFtHkFpH9svZfgQF+uPs6T9RvVpZeK8ajpCzy+H/5xKHFXhyKYJOLXUwfRnWLfBCjRrVEOgvr+NsX/3UeyH743i+sZJ9cTTdTMk2xsfH6Z8WcoZBxzgGHaJ8wYb4UysRMOwl+PZ9ENdnV0BqTJA+7/6RM4k2b1yPD997CxXKfI6qXxbVqojpqly+MF5/7WWULl7A4fysqlqFYihR9BO8+fp/UKlcYVStYH47PCoVU/tR4KO3s/V8qUrlCqn9KFP8U4fzK2r7+elHb6F8yQLaNhdC1fJf3FlffoFSRT9CoU/fRrUvMy2jzStb4hNtHf+7bb2Zq4pWb7/1Cop8/iE8KhZzuExWVb7053j1lZdQsay2DQ7mGyl5Pz549w3MmTVT/8SQMww6xjHoEOUDSVf3IHD0m2pcTsjMcki+flqfc39J0Fm5Yjm6tfVExIUViPNdj5ira01V7LV1iLy0CjUqFcEfW8cjwc/8OhIDN2H3ujGoW6MkwtV2rHO4nKtKCt6Cft82xOgBbZHgvxGxDpZxVXHafoRfWInqFQvjsPcEp/sRr71HsqyjeVJxfhuxbMaP6Kq9p8nBm+/YDnm/ZB0Zp2UsmR93bS1aNKyI1fP6ISlos8PlXJWs4+Ifc1Gh9KcIPr3M5fY6K1lHzNXVGNCrCSZNnKB/YsgZBh3jGHSILC7xyi6ETPocfn0fRPDUokgKPK7Puf8k6KxauRy9OtXTtmMjEObtcGyJy7qxDSnaAb2ORwmc3DM9e+uQe0B5T0SjOmWRFKBth7ZOh8u5KFvMTgz9sTkmjugIRGx3uIzL0l4zKXAzansUx+n9M7O3H1rZwn2wYnZfdGlTE4je6XCZrCpNq3bNqmHLsiHZW4e2L0GnlqhWJAktCDP/fso6bKG/YkS/VlrQmah/YsgZBh3jGHSILCwtKQZhK5rD9/u/IOjn95B4dZ8+J3dkDDqJARvSD26OzhZyUXJQTA7apILOn7u1oJOddURux8FtE1TQSfTXtkNbp6PlXFb0Dgz5oTkmDNeCjhY2HC7jouQ1EwM2pQedfdoXcDb2Q5UWkG4GHS3AOVzGRUGrtOvpQWfz0sHZW4e2L4EnF6ugE31ljbYv5t/P9L/BFgzvy6BjBIOOcQw6RBaXHHgMN5Y1RMLZ9fqU3MOgc6sYdG4vBh1zGHSMY9Ahyg+0gOEOGHRuFYPO7cWgYw6DjnEMOkR03zDo3CoGnduLQcccBh3jGHSI6L5h0LlVDDq3F4OOOQw6xjHoENF9w6BzqwwFnZBMjx0Vg06+dPXqVTRp0gQBAQH6FHKGQYeI7hsGnVuVVdBJC9mivcZ2daq2+j3DvNuKQSffkX9Hp0+fRrVq1XDx4kV9KjnDoENE9w2Dzq1yGXQk2MTswthhHXBijzYv0sX6GXTypXPnzsHT0xNXrlzRp5AzDDpEdN8w6Nwql0FHSgswZ3+fg/BLq13vI4NOviQBp3Hjxuy6MoBBh4juGwadW5VV0JHuKtlO6brKOP2OYtDJlzgY2TgGHSK6bxh0blVWQUcVByOTEww6xjHoENF9c7eDzkkZv6Id6GWaqdIO5nKvq8Z1y6l1qXVo22KmbHG7MKxPC0we1Vlb3870+zuZqXBvpARvQZ3qJXDuwGwgwsfxclmUTQttq+f1R7d2nkDcbofLZFVpYT74uoUHti4fCsTucriMy9L2JeTMMnX3dLlRKyK099PRcq5KW4ctbKu619XECbypZ1YYdIxj0CGi+8YedL7rWh9poVpokRtISveMmZIWB+3AWK9GyfSWEBUQ5MBqomJ24sj2yWhYqwyirq5VN7VMDfnVcKVpgSst0gcDejfFmIHtkBiyRd0RPSFwo+FKDNqMiEurUaNyURz2mYTE4M0Ol3NV6jW15y2Y8p0KKonafiRqwc3Rss5KtiPGdwNaNaqE1fP7I0ULcYkOlnNV8poXD85FpTKfaYHnl+zti7YdcQHrMaxva0yfNk3/xJAzDDrGMegQ0X0jQWfN6hWoUr4Qpo7uhNnjv8Wscd1NVg9M/6kbPv3oLXRtVx0j+jTHsB/N1ah+LdGhRRW8/spz6NSqGvr2qI8fv/nKcPXpXh8/dK2nbcMbKPjx6+japjq6t/fEN+1qGi5Zvk2TCnjhuSdRq0oR7bHj5VyVrKNbmxooVfh9vPX6c+jSxgM9TG5Hj6898XWzytp78SzKl/xYbUd29qWhZ0k8++8n0KJ+eXTX1uloOVclr9upZVWUKV4A8+fN1T8x5AyDjnEMOkR030jQWasFnWKF3kevznW1gNFICw0NTdcPXRvg7Tf+A89atdGmTRu0btUqG9UazZs1RYOGDVG/QQPTJc+rV+8r1NK2oVatOmpbzFcd7bm1UKt2dp+vlfZceb6sR9bncBkXVdNTfsrz0386WsZIyXuQ032pXqMmSpcuq1r9yDX7WVfXr1/Xp5AzaWlpGDp0KIMOEd176V1XK1TXlS27XVfR0nW1FfVqlMbJP4/JWrX1ppkqyM+0VCQmxCMhPi5blZgYj6ioSISFhWk/oxAZGZmtunHjhsPpRkteOzw8XK0nu9shz5Pny3pya1/s+zFhwgQsWrQo/QNDTkmLTp06dbBr1y6cP39eXVeH5bjk4oqdOnXCtDzSJcqgQ5SH3ZXByDe2IkUGI1cvhVMnT+hrJiuQz8ewYcMwfvx4fQo5I7eAKFy4MKpWrYqmTZuq1h3WnSW3yZAuvrJly2L16tX6u+feGHSI8rC7EnRC9bOutKBz8s/j+prJCuTzMXz4cEyaNEmfQs5Ii85XX32FY8eOqe6r4OBglosKDQ1FUlKS/u65NwYdojyMQYdcYdAxzj4YOSgoSJ9CVsGgQ5SHMeiQKww6xvGsK+ti0CHKwxh0yBUGHeMYdKyLQYcoD2PQIVcYdIxj0LEuBh2iPIxBh1xh0DGOQce6LBt0ZDS4lFzYyBn5EkhMTERycrI+xTEj65J1yLrkipH3k4x8//PPP11u2/2UkpKi3gcpec/kfZGf9mlSOWH0b5aZvD9ZPU/Wbf9by+/O2PfJ0Xsuf//M+5vV+nJC1sugQ87I54NBxxgGHeuyVNCRf9SHDx9Gw4YN8corr+DRRx9F8eLF8csvv9wWQORgvHnzZpQoUQL/+Mc/8PTTT6srkMpFkDI6cuSIOt3wiSeewD//+U8UK1YM27Zt0+emCwgIQPfu3fHss8+qZd599111tUg5uN1rsr/t27fHX/7yF2zfvl2faoxcSOzkyZOIj4/Xp+ScHPjlvZD39F//+pd6Px5//HFVTz75pKrXXnvNdEixW79+vfqbyXplXbVr18aZM2f0uY7J32HOnDn4+OOP1XOee+45tGrVSp0emZGPjw/Kly+vlpHtl2tpHDp0SJ+b7uzZs+qLUP7W8pn4/PPPsWrVqpshRt5LDw+P2/ZXSj4TBw8eVMvcbQw65AqDjnEMOtZlqaATExODatWqoVKlSvj5558xePBgvPHGG+rAs3btWn0pYOXKlXjsscdQuXJljBs3Dv369cOLL76IIkWK3Dy1UEKAHJQ//fRTjBgxAqNGjVK/P/PMMzh69KhaJiIiQh145bl9+vRRVyCV137ooYfu25VIN2zYgB49esDX11efYoxc0fLf//43/vjjD31KzknQkcAlFyeT9cvFpSSE/fDDD5gyZQomT56sQmd2Wp8kmD744IMoV66c+pt5eXnh9ddfx0cffaSupOuIvI687v/93/+pUCRf9t9++60KYfI5sYe8vXv33gzFY8eOVZ+bd955B//73/9ufuldu3YNb7/9Nt566y0MGTIEY8aMQalSpdS67CFTrkQrnwfZxqlTp6rXltecP3/+HcHqbmHQIVcYdIxj0LEuSwUd+d/7iRMnEBsbq0+BunKjHGyl5UPIvAoVKqBAgQLq8up2s2fPxl//+lfMnDlTPZaWieeffx6///67eix+++039T90uReQWLhwoQo18ly7kJAQfPbZZ2r92Tmg3y9yEH7hhRfUPt0rckVWCZTSepQT8mUtIfSTTz657RoX8+bNU+HH/jfLTILof//7XxVq7C1ssi4JtrJdW7duVY8l8Eogli86+zLyuZFlRo8eraYNGDBA/a29vb3VYyGXQpeWQ/lylL91dHS0CjoDBw7Ul7j3ZFsZdMgZ+Xww6BjDoGNdlh+MLFe5lKDToUMH9Vj+Zy4Hx8wHIwlAf/vb31SIkS8H+Z+7HCAzdrPIdAkw0hUhvzdr1kyFhYw3gZPp0sIirQhyk7jMZH6XLl3Qtm1bNbZGXkPWUahQIWzatElfKp1sk7QMSLfLSy+9pF5X7hgr43LsZs2ahQ8//FDtl5D9lXVt3LhRHZyltUkO9iNHjrw5tkTuUSLdL3//+9/xn//8R7VcSKuDkO4Zufy57L+0akjXXXYCW1xcnOoOlEuqS0uHnWyDBExHJS0zjsY4BQYGqr/ZN998o09JJ92Gsm81a9bUp9xOup7kbz99+nR9Srrjx4+r9UkrnIyheeqpp1TrU0byHsv7Kq02QlrqXn31VfW7nWyrfAakVUf2V1oUJej06tXr5vgcIy5duoR27drhzTffVK1Gsj9GuxTl88SgQ87I54NBxxgGHeuyfNCRLiQJHdLVIOQgJwc/+4HdTgKABIm6deuqA7P8771Fixb63HTypSEHRAkBclCWe33IgTaziRMnqgNpxtYgO1mH3DhODo4SmuRAKQdcCTPSgmB/jhzo6tevr7rdmjdvrgKObI88r3r16jcPonIHWVnm1KlT6rF0w8gYE+nWkTEkffv2VS0Wss8SiuT15TVat26t1iXdcuvWrcPFixfVfku3S9GiRbF8+XIsXrwY33//vRrTZJa0bEmAkvEwGcOLjHGScJW5JMh98MEHuHDhgr7kLTK+RcbESJdVRrIvEgykBc0R6aJ85JFH1NiejCSUSKiVcOHn56e68CQUZiThTrqyJHzI30JalKpUqaLPTSevL++vfL6kNUdKQqa0BMqyNWrUQP/+/VUrozMSZj09PVW36LJly1R99913t7VKuiLbcDPoBG4EwrxVcDFVWjhKCd7MoGNB8vlg0DGGQce6LB10pOviiy++UAcr+xgWGasjB33p9shIgk7BggVVuJCWFmntkICRkXxpSIuPBCJp+ZB1Z/5fvpCuFBnz4Wj8i6xDwpQMil2yZIk+Nf1gLi07EgzEli1b1DoytzzJAVm2XwbPCvkSk+40+0Bq6YqS+XLwlLAhpMVJgpQceO1kwLS09mQccCthScYlyTiVnJLgJO+hjFvKSL5MJAjK9MwlLS+OurlkwK8EFkdB57333lOBzhHZj4cfflgFuYwkxDzwwAP4+uuvVeiTYJg56Mi6S5cufXOcjvzNZYByRrKMjOeRoCN3i5bWKunGlNeVMWK9e/dWf1MJnXJHZEfk/ZB1y5ij7JBtkLuX5zjoBG1GXQYdy5HPB4OOMQw61mXZoCP/wLt166YOaHKgtJPBu85adKRlQLpqJDTIQdpRi450OcmBScZnyP/eHQUdaTmRkOKsRUdCiBycM/6vXQ6S0s0jY0WEDJqVlhq5o25GEp7k4PnTTz+px5mDzr59+9TBPWOQk9eU/ZLttnM0RsfeiiQtHNL9dv78eX2OeStWrFBjnn799Vd9SvZJOJWgk/kOzLJfEnSctejI8vJeZG7RsQcdadGRgClBJ3OglHWXKVNGhWTpInv//fdVN2NGsow96EiodkQGrktrVIMGDZCQkKBPvUU+dxKe5cy/jh07qrBppqtQtmHdmhV4/51XUatqUdT3LJ2tqluzFF568RmcOc2gYyXy+WDQMYZBx7osG3Tkf9RyGnLmgarSEiNBR1o0MpLuGemSkgOSHHyk6ynzuA350pDAIMtJKJCuK+lyyTyuRA6wciCVbrLMZB0SdKQ1SFoB7GS6tCBIK4uQ7ZCurIxjhISsU8KVdDkJR0FHQlbGICfrlvVJcLJzNhhZuq/krCJpAZIzzCRwGe1GsZPXky4v2U45ey0jaV2S8CMtVhlLpkkrlaPXkrAgQce+z3byd5JxRNLF5Ii05EgXpHRhZSSvIX/fzp07q8HNEjKkuygjed9lvRJwpJtQLi1QsmRJfW462U/pdpS/dca/ZUbyuZIgJqFYxho5ItsjwVUGW0vIlMCTcVyTK7INa1avRNOvKmDfxp9wyHsi/tg6wVQd3DYRv20eh7IlPsGpP513s1HeI58PBh1jGHSsy5JBR7pB5GAv41oyjy+RQbsy+FRaezKSAaESgGRMhXw5yAFOBqJmPODIGAwJOfI/fdGyZUt1YMrY8mEfoCrhw9FBW9YtQUcOfhnXLS0CcqCTVh0hp2RL0LGPvbGTMS4ZA5yzoJO5RSdz0JFTvSXoHDhwQJ9yO+lCknFN0iIiA6LNkNeTcUQSFGSQcUZymriMDZJtzliyTxKMpKUsMzk1WwKLvK8ZSeiT8TDSReSIBCT5m8p7ZCfbJsFKWlnkPZQBxBLoMrfWyHZIl5OcrSd/U7k+jnxuMg4Slr+vjMWR1jlnA4/l8yetQnKNnqzOPpPPg/xd5O8kg8eNkP2RMTq9O3+FlOBN2gfJBwjbZq7CvWG7vgV1a7Drymrk88GgYwyDjnVZKujIP2o5UMiBvmfPnmrQqX26lJDuA2mVkYPr/v371XQ52EmXjZyJZO9ukkG+cnCVlh9ZRroT5AtD/ve+YMECtYx0g8mgVmkZkNYFWU6myUE08xlCdrKMvL6EGDnNWQ6Wsm7Zblm3tLSInTt3qmXk4C7LyPNkzI0ckGX99n+M2Q06c+fOVeFCgoe8vr0ytk7JwV5amOS9ENI6I6+fuSsoM3k9CQXypSG/ZyTvtYS3zCXrlov/OQoMsg4JThKQpNVHHkvgkPFM8v7LmWZC1iGXDrC3Zsn7JgO+JWjI+yPPk+AlQVVCnnRJybQmTZqo91q6yOSx/C3l8gLSfSmDu4W00kloktYuWUbeK/kcSBCUeTJNznST6wfJPHksnzUJ3fI86RqTaZnJshnDuHRVSsDK3MLkjKyTZ12RM/L5YNAxhkHHuiwVdGRgrXRDyIFFuobkTBnpMpDWEznzSLqthAQcGWQqQUbOMJIzp+R/9XKQspNWBOlKkmXsp5RLoOjatas6OAn5KQdbCRbS0iOvI2FC/vfv7AJx8sVTr1499XrSgiODf6X1SNYhB2AJAvblJGBIYJCDsqxb9k1+zziIWcaISECyBx05MMv+Z7y2j6xLTlmWMTN28o9atkHCg1x0T86w2rNnj9pn+ccu41fkLCgpGVgsZNyNrNt+Bpszsi2ynD0g3Q3SeiN/J3mf5Awled9kQLcERDs5U0xeN+PYKgkf9qtWS8uLdDVKeMs4QFn2T/5+0j0mfxMJRvKcjK0qEqwqVqyoPg9yppf8zeS9k1Y9aemT91iCj7TwyTgr6eqS7ZW/mVyp29kYHhlzJX93CaLSemR/fRkQb4S8LoMOOSOfDwYdYxh0rMtSQUe6pWTMhIxxkIOGHKztJS08Ga8eLKcxS9eWzJOuDwk/mUlXg/yPXJaRVps1a9aoL46M5LGEBLk2jSwn3SGuxlfI8tJ1JQFLulbkS0i2Vc44cnSFX7lAnQwMlnXL9sotLjKSlh/ZfnsXkeyjjKvJvJwEHxk3k5E8V7rw5L2R1htpMZLgIK8lg65l2+Qfv50c+OVAntXVlKUbUF7L6MHaKNkWaRmR7ZNtzjy+SK5bJPuTeQC0vM+yPfI86ZrMPG5IyJebXOBQlpHWHDlLKvPfWloIpRtPrskkfzMJfhm7sqQFR/5ecsq5rEdeU7bFHowdkb+btBrKGWDynkuLkaPuO2dkGxl0yBn5fDDoGMOgY12WHYzsruxBR1oXpCUgr5BBu9KSImGL3AeDDrnCoGMcg451MejcZ/LFI9djkavuOuvOcEfSaiXjmDK2YFDuY9AhVxh0jGPQsS4GnVwgF8z78ccfDd8igMgZBh1yhUHHOAYd62LQIcrDGHTIFQYd4xh0rItBhygPY9AhVxh0jGPQsS4GHaI8jEGHXGHQMY5Bx7oYdIjyMAYdcoVBxzgGHeti0CHKwxh0yBUGHeMYdKyLQYcoD2PQIVcYdIxj0LEuBh2iPIxBh1xh0DGOQce6GHSI8ohTV69j2NKd2Hvyqj7FHnRW8O7l5BCDjnEMOtbFoEOUB2jHK0xY9xvebz0W30zZqE9NP5CtWb1SCymlsW7BAGz5ZSg2Lx1irpYNxYZFA1Gs8Mc4deqEvmayAgYd4xh0rItBhygPkKAzesUefNJ+AjpNuHXndTmQrV29AgU+fAstG1RA+2ZV0a5pFXOlPad1kyp49ZUXcOYUW3SshEHHOAYd62LQIcoDJOiM+GU3Cnw9AR3Gr9Wnph/IpOvqu671YQv9FYjeAURuN1dRO2C7sRX1apZm15XFMOgYx6BjXQw6RHlAWpoNgxdtx6cdJqLN2NVISE5V09ODTvpg5KTAjeljbkK3mqsb29T4Hg5Gth4GHeMYdKyLQYcoD4iNT0TvmZvxWcdJaDJyOYIjYtT0jEGHZ11RZgw6xjHoWBeDDlEeEBoZi3Y/r0ahLpNRZ9BinPULVdPlQLZ21XL07lQHqcEbgQhvIGyruQqXcLQZdTxK4AwHI1sKg45xDDrWxaBDbiMhKRkbfz+LLQfPsTLUtiMXMGPjQVT4fjZKdJ+Osj1nYuQvu+B95KI27yKGTF2IFl27YMO2X/DrzlXYsmOlqfp150ps8l4Oj8bNMXuND349dN7hduSF2nDgDILCo/VPFDHoGMegY10MOuQ2RizbjdLfTtdqBkpqB/SSPVj2Kv7NNBTtNlX9LKb9LKb9LKXPK9VjBkp9O1P7XXvfclCl1TrufO08UdrnRd4P+fx0m7IBKWlp+qcqf2PQMY5Bx7oYdMgtxCYmofaARSjadSo8vRai5ZgVaD5qOetmrVBjc4Yu2YkWo+W9SX9/Wo5ZiVr9ZqNyj7FoNnQOmg+fm82ag3JdRqHBoIVortbvaBvct+Q9qT90iQo95XrNQlAYW3UEg45xDDrWxaBDbuHgOX9U+mEOKv84F5sPnkVAaBR8r0eyMtTVkAikpKbdNs0vNBKzFi5Fl87NcenkGvid3wTfsxtMld/5jbh8ei0861TD9n2/q3VmfI28UP7aNh8854f6Q5agxDfTsen3s/onK39j0DGOQce6GHTILaw/cFZ1wTSTM4o4xsIwOZBxMHK6+MRk9Jq5BYW7TMa0jX/oU/M3Bh3jGHSsi0GH3MLcbUdQuOsUdJiwBjEJifpUyoocyHh6ebqUlDQMXboTn3WahMGLt+tT8zcGHeMYdKyLQYfcwvg1v6lrxPSYtlF9OZMxDDq3yMdm/Or96qKKPaZv0qfmbww6xjHoWBeDDrmF0cv3qgNU71lb9ClkBIPOLRJ0pqw/gI/bj7/tNhn5GYOOcQw61sWgQ27Ba7433mw+Bt2n3bozN2WNQecWaQdctP0YyvWaiZ5s0VEYdIxj0LEuBh1yCzuOXcL4Nfux9fAFfQoZcVvQUfe68lbBxVRp4SgleDPq1iiZ528BERQRg8Pn/XFBv3J0fsegYxyDjnUx6BDlYXIgW7N6lRZ0vkJCwEbYbnjDdn2ruQrVgk7QZtSpwXtdWQ2DjnEMOtbFoEOUh6WmpmLJooUo9Ok76NjSA9+0r4Vu7TxN1TftPdGljSfeeO1lnDnFoGMlDDrGMehYF4MOUR4WHR2FenVro+qXX+CX2X2wcl5/rJzTz1zN7YdfZvVB4c/ew2kGHUth0DGOQce6GHSI8rDk5GSMHz8e33X5Cgj3BmJ3aelnh7mK2QlbhDfqe5Zh15XFMOgYx6BjXQw6RHmYHMjWrVuD3jzrihxg0DGOQce6GHSI8jA5kPH0cnKGQcc4Bh3rYtAhysPkQLZ65Qp11lX66eXbbj913EjJ6eVa0KmrBZ1TJ/Puva7oTgw6xjHoWBeDDuVYapoNySmpiE1MxuXgcBy5EACfoxexYvefmLXlECau/Q0jf9kNrwU++G7mFvSYtgnfTNmAThPWovvUjTh5JVhfE5llDzo9O9ZDUpAWdGScjhZcHJaEIPV7eri5Nd0bqcGb1E09T51ki46VMOgYx6BjXQw6ZEqaFmqCw2Nw1jcUmw+eUxf5662FlyYjlqPCd7NRtudMdRfyEt2noWi3qSjcZQoKdZ6MzztNUveykts8SBX4egI+bjcexbVl9py8oq/9ltTUFMydPQsjhg3C2DHD8dNo8zV61DCMHDHY4TyjJc8fOXwwxo4e4XC+kRo+dCDGjBqa7f2QGj1yiLYtQ+6YPnHcaDRr0ggli3yEUV6tMXZwe/w0sN0dNUarEf1aa8u0vXP+oPbque+//Qp69uiqbeed+zpWqzEjh6m/x0+jbp9npkY52Q+jJdsxSnv+pInjEBcbo39ayBkGHeMuM+hYFoMOGZKQlAKfIxcxavlutBy9EqW/nYFCXSanBxnt5xdamJG7j3+phZ1aAxaiycjl+Hr8GtV688PsrRiwwAfDlu7E2JV7MW71PkxYu1+19MzafAjXQiL0V7klLjYapUsURscWlTB+aDuMHdDKVP08qA26tauOd996GWO8WmLsQC0EOFjOacny2joqly2ASqU/xoTh7R0v56J+HtQaQ79vgnfeeAHftK+h7Udbh8tlVeOHtEXtqkVQttiHGDe4zW378pP8HKStd/DX2k9tGwdp71WmGjekPYb3bYkvPn0HLRpUwPjhHe9Y5ietftbW8dPANunrzFSy7d3aVsfH7/4XA3o2xPjB5vdF9qNWFW0/imv7MeT2/TBaso4futbGZwXehe+1OwMy3Y5BxzgGHeti0KEs+d2IUneDLtNzJj7rOBFFu05Fye7TUHvgItWaM2ndAazcexL7Tl7FicvBuOB/A77XI3E9MhZRcYmIS0xWXVs2dTciY2Kio9CskSfO7J8BpOxPP23aTCXswel9M9C4bjl1+jTidjtezllpy9vid2PSiE6YMKyjtkW/O17OVWnPTw3ejOb1v8TZ32Znbz+0siXtw8IpvTGsbytA+930vmjvhQw67t25HrxXjdCOfgccL+eqtG0/+9ssdGjpgagra9K3w9FyLkr2Y/6kXhjeL5v7IaU978b55WjaoBquXrmc/mEhpxh0jGPQsS4GHXJJuqoGLdqOQl2mqK6pLpPWY/H2YzjvH4qY+CQkSYAxnl8Mk6DTuEENHPGZAERuv+NMoaxKxqMc8pmE+p6lkRayRQ26dbScs1LL39iqWjpGD2irDrKOlnNVMhZGQkGjOmVxdMeUbO2HqggfzPr5Gwzo3RTQfje9L9p7EXdtHbp/XQsblwwGonc6XM5VybYf3T4ZbZtUQfCpZUC4j8PlXJa27TPGfoOB3zXL1n5IyfP8ji9Eo3pVtKDDFp2sMOgYx6BjXQw65FJ8YgrK9JyluqWW7DimT733GHQyFIPOzWLQMYdBxzgGHeti0CGXJOiU7z1bDSg+fCFAn3rvMehkKAadm8WgYw6DjnEMOtbFoEMuZQw6f5zz06feeww6GYpB52Yx6JjDoGMcg451MeiQSww6DDpSDDp5E4OOcQw61sWgQy4x6DDoSDHo5E0MOsYx6FgXgw65FJ+UrAWdWSroHDx3/74AYqIjtaBTUws6E7MddA56T0TDWmWQGpy9oCNXDR4zsK0qxGQn6GxD5OXVKugc25n9oCPPm/FTNwz6XgsIcuXjbAadbu08sXnZkGwHncPek9CuaVUEn85e0JF1TBvTFUN+aJ6t/ZCSoON7bCEa1q0M32tX9U8LOSNBZ8iQIZg6dao+hZyRoNOwYUMEBQXpU8gqGHTIJQYd9w460mokrU2u9o9BJ/9KTk5G9+7dVasOucagY10MOuQSg44bBx1t/cGXVmPNokGID9ykutoyP1eKQSf/YouOcQw61sWgQy4x6Lhv0JH139DWv3HJICQEbFSPMz9XLcegk29J0Bk6dCimTZumTyFnGHSsi0GHXGLQcfOuK9lO7eCfcfnMxaCTfzHoGMegY10MOuQSgw4HI0sx6ORNDDrGMehYF4MOucSgw6AjxaCTNzHoGMegY10MOuRSbgadpg09cXKv9gWduBeI2mGuYnfhxK5paCJ3L9cOsOoO5tHadKMly2vrGD+sA8YN7QCk/OZ4OVelPT8+YCOa1iuPk/tnZm8/pBL2Yu6EbzGsT8v0O35rQcXhcs5K247E4M3o0aE2tq4cASTvd7ycq9K2/dju6ejQohpCL63WtmOP4+VcVeIezBrXHSPk7uXx2dgPqfg9CDyzDE3qV8W1a7yOTlYYdIxj0LEuBh1yKbeCTmxMFOp4VsH8Sd/i0I4p+GPrBFN10HuS9txeKFP8I2xZOghbfxmMX5cZL1nee/kQtGr4JWpWLow18/pi1dw+WDnnR8Mlz5k2qhM+fu8VfN+1LhZP640F2v6YqYWTpXqibvXiqFC6gBZ4umOB9tjRss5q8ZSemDi8Awp+9AbaNK6IpdO/c7icq1oyvTd6daqDzz5+AyP7t8Tiqeb2xb4ftasVRaUyBTBvovn9kFo0pRdGaa9fs3olBAbcvwtY5lUMOsYx6FgXgw65lFtBJyE+Cl+WK4malQqhY8vqaN+8qqn6ukU11KlWHC+9+Dw8qlZAda08qhiv6lrV0J5TpnQJFC70ufazOMqW0Up+Gq0yJVCqZDF8+OEH+OyzT1G06BcoWqRQtqpAgU/wyccfa7+bX0cx7XVlH2Q7ChYsgGIOlsmqZB2yDx988AEKffGZeuxouawqJ/shJa9b4JOP8OWX5REawgNSVuxBZ9asWfoUckauiCxXRmbQsR4GHXIpN1t0GtTzgM+qYbh+aTVCzv5iqq6fXwHvFUNRw6MKLl+8iKtXr+LKlcumSm4xcPnSRZw5exqnTp3KVp0+fRqHDx/S6kiO6tChw6oOHznqcH5WdeToUe35sh3aOhzMN1RHjuCgrCOb2yCV0/04duw4fv11q7oIXmBgoP5pIVe8vLxQsWJFddHAwYMHsxyUhMEuXbqgUKFCCAkJ0d85sgoGHXIpN8foNGlYEyd2TQHidqUPXjVTUXJvpklo3LCevkayCgk4DRo0wLVr1/Qp5MqAAQPg6empuq/knlcsxzV+/HjMmzcPiYmJ+jtHVsGgQy7lZtDJ8VlX28ajYf06SEtN1tdKVuDn56fGUvj6+upTyBl719XMmTP1KUT5D4MOucSgQ+6GQcc4DkYmYtChLDDokLth0DGOQYeIQYeywKBD7oZBxzgGHSIGHcoCgw65GwYd4xh0iBh0KAsMOuRuGHSMY9AhYtChLDDo5B+LFy9GVFSU/sh9MegYx6BDxKBDWWDQyT9WrFiB6Oho/ZH7YtAxjkGHiEGHssCgQ+6GQcc4Bh0iBh3KAoMOuRsGHeMYdIgYdCgLDDrkbhh0jGPQIWLQoSww6JC7YdAxjkGHiEGHspCbQadJQ0/8uXuqthF7tLDjY66id+DYjslo1KAug47FMOgYx6BDxKBDWcitoBMXG4WqlcuhZ8famDqmKyaO6GSqJo/sjO5f14JnTQ/YGHQshUHHOAYdIgYdykJuBp0qlcqiY8tqGDOwHUb0a22qRnq1wdfNq6BmjaoMOhbDoGMcgw4Rgw5lIde7rnZN0VLPbiDCx1xJ19X2SWjEMTqWw6BjHIMOEYMOZSG3ByMflcHIUduB0F/NVfg2HOJgZEti0DGOQYeIQYeykNtBh2ddUWYMOsYx6BAx6FAWGHTI3TDoGMegQ8SgQ1lg0CF3w6BjHIMOEYMOZYFBh9wNg45xDDpEDDqUBQYdcjcMOsYx6BAx6FAWGHTI3TDoGMegQ8SgQ1lg0CF3w6BjHIMOEYMOZYFBh9wNg45xDDpEDDqUBQYdcjcMOsYx6BAx6FAWGHTI3TDoGMegQ8SgQ1lg0CF3w6BjHIMOEYMOZSE3g04TLegc3zVZe7BLCy7e5koLR0e2T2TQsSAGHeMYdIgYdCgLuRd0otDoq+rYvW4UIv3WI+LiKlMVeWUNdq4ZgQZfMehYDYOOcQw6RAw6lIXcCjoJcVGoUK4Ualctgq5tPdGpVXVT1bl1DdSsXBgeHpVhS03S10pWwKBjHIMOEYMOZSG3gk5sTBRq16yM6WM6Y9fGsdixdrSp2rVuDKaM7Ii6tauzRcdiGHSMY9AhYtChLOTmGJ2mDT1xap/2BZ20F4jeYa7iduHP3VPQqAG7rqyGQcc4Bh0iBh3KQm4GHZ51RY4w6BjHoEPEoENZYNAhd8OgYxyDDhGDDmWBQYfcDYOOcQw6RAw6lAUGHXI3DDrGMegQMehQFhh0yN0w6BjHoEPEoENZYNAhd8OgYxyDDhGDDmWBQYfcDYOOcQw6RAw6lAUGHXI3DDrGMegQMehQFhh0yN0w6BjHoEPEoGNIWloaAgICEBsbq0/JPxh0yN0w6BjHoENkMuhERERg27ZtWLRoEX755RecOnUKqamp+lzgxo0b8PHxUfO2bNmiwkFmEhqOHTuGxYsXq+VOnDihz7lTVFQU1q1bp5bPTcHBwfjLX/6CcePG6VPyDwYdcjcMOsYx6BCZCDq//fYbypUrh7fffhulSpVCgQIFULp06ZtfNjt27EDJkiXx/PPP45///KcKBsWKFVPPs4uJicEPP/yAt956C59//jnee+89vPrqqxg8eDBSUlL0pdKdPHkSnp6eaj1NmzbVp+YOe9AZM2aMPiX/YNAhd8OgYxyDDpHBoCNfKP/5z39UwDl9+jQiIyPVwV/CSHx8vFpm/PjxaNWqlWrluXTpkmr9kHDw5ZdfIi4uTi2zYMEC/O1vf8OECRNU609gYCC6dOmCRx99VAUluyVLlqjXq1ixIv7v//4PzZs31+fkDgYdBh1yHww6xjHoEBkMOhJGnnnmGZw9e1afcifpkspMwsr//vc/XLt2TT2uUKHCbcFHXL16FU899RTatm2rHicnJ6NRo0YYNGgQEhIS8Mgjj6BZs2ZqXlakG02CmISm/fv3Izw8XE3/448/cO7cOfV7RjLt6NGj6stASGg7fPgwtm/fjiNHjtwck+Ms6Mj8Q4cOYevWrarlKuN+2Z0/fx7e3t6qy+/ChQt3tFy5u9wMOnL38jP7ZwAp+7U3e5e5StiDU3umolGDugw6FsOgYxyDDpGBoJOYmIhChQqhWrVq+hTjXnrpJXzyyScqcMg/uI8//hh16tRBUlKSvkS61157TXV7SciR5UJDQ/U5wMMPP2wo6Mj4oW+//RbPPfccHn/8cfW8smXL4sqVK6hSpQoKFy6sWpHs5LU++ugj1R0nIU0C11dffaW63f7xj3+ogNWhQwe1XtmezEHnzJkzqFWrFt544w0ULFgQL7/8MurXr69as4Tsh4xB+uCDD/Dhhx+qn9IVJ9uTl+RW0ImLjYJHlfL4odtXmD2+B6b/1M1UzRz7DXp3rgPPmh5a0Ln980Z5G4OOcQw6RAaCjrSQvPnmm6pVZ8qUKfDw8FDhQB7/+eef+lJ3OnjwoOqSqlevngoS8g9OgsF///tfFRLspBVFApCEBXsLTEZGg87w4cNVGPnxxx9VK82ePXvQpEkT1b0mg6cffPBB/Prrr/rSUK02Tz75JObPn69aWWTbZHvlsQx+Xr16NXr06KECkASkjEEnOjpaBTYZY3TgwAFcv35d7a/sm4Qd2V9pjZKWqk6dOqkuOn9/fxw/flyNU8pLcivoxMdFo3LFsmjbpDKG92mFQd83M1VDfmiBNo0qoEb1KrCxRcdSGHSMY9AhMhB0pAtHBgy/+OKLKgzIP5o+ffqoLilpickYWuykC0cCkRzoJSjZSRfPs88+q1pBOnfujK+//lq1tMhy0mokrSeZGQk60mUl65XXzEimS+iQrjP7a9pJUJPXlfFG0v0kXXMDBw7U56azPz9z15W8JzJ2aObMmaplSEKN/JSB1q+//rrqopLHDzzwANq0aaPOHsurcq/rKgpNZIzO9omwRW5H2vVfTZXtxjYc3MoxOlbEoGMcgw6RwaAjLRVeXl76lHS7du3CY489hlGjRulTbpEgJF1Acgp5ZtJa8v3336sBxhI2ZPxKkSJF1EDnjKeq2xkJOtJdJF1N0prjiKxXXku6l6TVSEqCVYsWLdR8afF54okn1CnxjmQOOkuXLlWPpVtOxhxJF5m0ckkYfOGFF26eMj99+nT1WJaTFif7WKW8JDfH6EjQObpjEhC1Awjdaq7CvXHIewKDjgUx6BjHoENkIOhIiJBTyjt27KhPSSddOu+//z5at26tT0k3ceJE/Otf/zJ8zRnpypHWlHbt2ulTbmck6Mggael2ytwik9HGjRvV2J1NmzapbisJaXLNHzFnzhw1b+fOnepxZpmDzrJly9Tj7t27q1Yde82dOxcbNmxQXVt20qUlwU9CloQrV9197ig3gw7PuiJHGHSMY9AhMhB0ZOCwtLi8++67qhvHbu/evXjooYduHvxl3uTJk1VLjhzYjRo5cqTLkCFBp2XLlvojx2SMjbTI1KhRQ5+STqbbW4mk++iLL75QA5Nl0LK0sti7ymQ8j4SzAQMGqMd29udnDjpyFpc8XrFihXpshIwVkucMGTJEn5I3MOiQu2HQMY5Bh8hA0BGjR49WB+lvvvlGnbYtg3plbI101Vy+fFmFnNmzZ6vuI+nG+f3331VLhgQCWT4oKEitR6ZLy4ecii3T+/Xrp1pWMl8QUFqRZJCvPF8GEVetWlWtT7rRnJFuq7///e+qi02eK9so65UBwHYSZOQ6PhKspPVHvgSEnFlWpkwZdbaVdLfJdsq1fCQQyVlSYWFhtwUdGZxcqVIlNU5JXkcGGstp5PIcadGyLyMtRnIKu8xfuHAh/vrXv6prCeUlDDrkbhh0jGPQITIYdKQrRs4eklYPaWGRlpzy5cvfvMifDD6WMCJhQEKEBAZpYZGSZe3dWOvXr1dXTpZlpKtJrrPTu3dv1WJiJ6FJxu9IAJLnyzqlZJ0y0NcZWYdcsFCeJ+uWliW54KB8KdpJGJEzrWQQcsYrNgtpcZFT6CUsyfbJMjKuR1qC7KeX//TTT/rS6cvLGWXSJSWtQ/ayf6FIwJIgJFeBlrPKZF8lOEloyksYdMjdMOgYx6BDZDDoCOnGkbEw0loiLSsZz5CScCKtMNLqIq0hEiLsJcvL6dVCusFkOVlGpss4H0ckkMh8eb69ZUh+l7OjXJH1y1le8lwJGo7OdpIzxaT7SlpxMpPlZbC0PF9aYuxXfZZ9l9eW08gzkoAnV4KWbjzZxoz39pKzrmQdst0yX967zNcPygsYdMjdMOgYx6BDZCLoWIFcnVhahuRsKDKGQYfcDYOOcQw6RPks6NSsWRNPP/20GjNDxjDokLth0DGOQYconwWdffv2qTOsyDgGHXI3DDrGMegQ5bOgQ+alaV+UEnS+6DQJf5xl0KHcx6BjHIMOEYMOuRAVl4gdxy+jTM+ZKNl9Oo5eTB9Ufj8w6JAzDDrGMegQMehQBmlpNlwJCsf2o5cwad1v6DxpPcr1moUvOk9G54nrEB6dfhba/cCgQ84w6BjHoEPEoGMZ0sUUm5B0W0mLTEx8kqro+ET1OCImAf43onDiShB2Hr+MZTtP4OdV+9BzxiY0Hr4M1fst0MLNTHzeaZKqcr1nofvUTTh1NUR/pfvDHnTUva6idwA3tpqrcG8c5r2uLIlBxzgGHSIGHUtISU3D7F8PofGIX9B81HKtVqDF6JXwHLAQjYYvRcNhS1F30GLU6L8AX/aejZI9pt+q7tNRQqvi30xTj6v+OE9bfhk6TliHsSv3Yv/p3LkRqdy9XILOoW0TtNDig9TgzaYq7fpW/P7rz1rQqcugYzESdBo1asSgYwCDDhGDjiVExSZoIWY+inSbimIZqkT3aSrA2EseS5gp1WNGen2bXqXlpza9wnez0W3yBizafgx+oXdebPF+igi7jtIli6Jm5SLo3LoGOrTwMFUdW3rAo+IXqFypAhLjY/S1khVcuHBBXfRTLshJrjHoEDHoWEJqmg2/7DqBtj+vRocJa2+WtMy0GbsKrX9ahaYjl6P+0KXw6DsfX2qBRs6kki6qMj3TQ449DMk8aQlqMmI5vOZ7Y+vhC6rF6H6LigpHTY8KmDSiI3asGY1tK4abKp9VIzBhaHvU9qyG1OQEfa1kBXJF9Tp16qj77JFrDDpEDDqWkpqWdltJQEnVSn7aKzklFRGxCbgcFI7DFwKw6Y9zmL3lEAYv2o6OWjhqMGwpqvw4F4W7TMFnnSapFqC2Y1fj0Pn7e5FFGaPTpKEn/tw9BYjfDUT4mKvoHTi+YxIaNWDXldWw68o4Bh0iBh3KJCgsGr+f9cMin2P4cc5WVP5hrrpYYLufVyM0MlZf6t7jWVfkDAcjG8egQ8SgQy4kJKWoVp+yPWeqMT9HtN/vFwYdcoZBxzgGHSIGHcqCdIHdvDIybwFBboBBxzgGHSIGHcoC73VF7oZBxzgGHSIGHcoCgw65GwYd4xh0iBh0KAsMOuRuJOg0btwYISH392rdedXIkSMZdChfY9Ahlxh0yN34+/vDw8MD8+fPx/bt27Ft2zaWk9q6dSsaNGiACRMm6O8eUf7DoEMuMeiQu5GgU6JECdWq061bN3Tp0oXlpDp37ozmzZvDx8dHf/eI8h8GHXKJQYfcjXRd1a9fH+fPn0d8fDzi4uJYLiohIQFpaff/6uZE7oJBh1xi0CF3Yx+MHBBw/67rRER5F4MOucSgQ+6GZ10RkRkMOuQSgw65GwYdIjKDQYdcYtAhd8OgQ0RmMOiQS7kZdJo29MSpfdOApL3qbuSmKm6XuvM5715uPQw6RGQGgw65lFtBJy42Cp4eFTF6QGusWzIIa+Z7maq1CwdgZL9WqF2rJmypSfpayQoYdIjIDAYdcim3gk5CXDQqlCuFOh7F8E37Wujcpqap6trWE55ViqB6tcqwpTHoWAmDDhGZwaBDLuVe11UUmjTwxJHtk2CLkjE6W81VmDcOqTE6tdl1ZTEMOkRkBoMOuZTbg5GPakEHUTuA0K3mKlwLOt4TOBjZghh0iMgMBh1yKbeDDs+6oswYdIjIDAYdcolBh9wNgw4RmcGgQy4x6JC7YdAhIjMYdMglBh1yNww6RGQGgw65xKBD7oZBh4jMYNAhlxh0yN0w6BCRGQw65BKDDrkbBh0iMoNBh1xi0CF3w6BDRGYw6JBLDDrkbhh0iMgMBh1yiUGH3A2DDhGZwaBDLjHokLth0CEiMxh0yKXcDjrZvteVuqkn73VlRQw6RGQGgw65lJtBp2lDT5zcM1XbiD2qVcdUxezE8R2T0ahBPQYdi2HQISIzGHTIpdwKOvFx0fiyXCnUqFQYHVtWR/vm1UxVhxbV4FHhC1SvVgW21CR9rWQFDDpEZAaDDrmUW0EnLjYKntUr4ufB7fDr8mHYsHiQqdq8dAjGDGiN2rWqa0GHLTpWwqBDRGYw6JBLuT1G5/jOyaobSgYXm6pIHxzZPpFjdCyIQYeIzGDQIZdyO+jwrCvKjEGHiMxg0CGXGHTI3TDoEJEZDDrkEoMOuRsGHSIyg0GHXGLQIXfDoENEZjDokEsMOuRuGHSIyAwGHXKJQYfcDYMOEZnBoEMuMeiQu2HQISIzGHTIJQYdcjcMOkRkBoMOucSgQ+6GQYeIzGDQIZcYdMjdMOgQkRkMOuQSgw65GwYdIjKDQYdcYtAhd8OgQ0RmMOiQS7kZdJo0rInju6YAsbuAcG9zpYWjo9snMehYEIMOEZnBoEMu5VrQiYlE/boe2LJsMILO/gK/E4tNVeCppdi0eCDq16vNoGMxDDpEZAaDDrmUW0EnIT4aZUsVR8nC7+OrGqVQ26O4qapbvQSKF3oPVatUhC01SV8rWQGDDhGZwaBDLuVe11UUGtSrjm0rhiLswioEn15mqq6fXa5ag+p/VYstOhbDoENEZjDokEu5OkanQU0c2zkZiN4J3NhmriJ8cNhnIsfoWBCDDhGZwaBDLuVm0OFZV+QIgw4RmcGgQy4x6JC7YdAhIjMYdMglBh1yNww6RGQGgw65xKBD7oZBh4jMYNAhlxh0yN0w6BCRGQw65BKDDrkbBh0iMoNBh1xi0CF3w6BDRGYw6JBLDDrkbhh0iMgMBh1yiUGH3A2DDhGZwaBDLjHokLth0CEiMxh0yCUGHXI3DDpEZAaDDrnEoEPuhkGHiMxg0CGXcjPoNG3oiTP7pwMp+7UJO81V/B6c2jMVjRrUZdCxGAYdIjKDQYdcyq2gExcbDY+qX6Lftw0xf0pvzJ3wramaN6kX+nxTH7U8q8OWlqSvlayAQYeIzGDQIZdyK+jEx0WhUoUyaF6/PLx6NUWfHo1MVb9vG6NpvbKo4VEZNrboWAqDDhGZwaBDLuVe11UUmsgYne0TYYsyP0bHFrYNh9QYndrsurIYBh0iMoNBh1zKzTE6EnSObp8ERO0AQreaq3BvHPKewMHIFsSgQ0RmMOiQS7kZdHjWFTnCoENEZjDokEsMOuRuGHSIyAwGHXKJQYfcDYMOEZnBoEMuMeiQu2HQISIzGHTIJQYdcjcMOkRkBoMOucSgQ+6GQYeIzGDQIZcYdMjdMOgQkRkMOuQSgw65GwYdIjKDQYdcYtAhd8OgQ0RmMOiQSxmDzpELgfrUe49Bh5xh0CEiMxh0yCUJOmV7zUSRrlMwd+th2PTp9xqDDjnDoENEZjDokEs2mw2jV+xBoS6TUfrbGWgxegWmbfwdhy8EIDg8BtFxiUhOSdOXvnsYdMgZBh0iMoNBh7IkgWbAAh9U+G42Pu80SbXuFO06FR5956HLpPUYvnQXZv96CBt+P4M9J6/iiBaCzvhex7WQSFyPiEVkbALiEpORkmo8ENmDzt24qactLUVfK1mBv78/gw4RGcagQ4YkayHlj7N+mLbxD3SYsFaN25HAU7jLFHzReTI+0wKQhCBp9anaZx7qDFqMZiOXo93Pq9F50jr0mLYR383agn7ztsFLC00DFm7HoEXb0X++NzYfPKe/yi0SdBrW88D+TT8hOXgz4nzXmqqEgA3Yu2E06taugdiYaCQkJCA+Pp6VxyspKQkXLlxAvXr1GHSIyBAGHTItKi4RgWHR2H/qmhq3M0gLLe3Hr0FNr4Uq5FT8fg7K9ZqJUj2mo/g301Qgkq4vCUIFO0zEp1oV+HqCqreaj1GBJ7PEhDhUrlgeRT5/Hx6ViqNqhWLmqmIxlCz2Cd5443VUrVoF1apVY1mgPDw8ULp0aXzxxRcICQnRPy1ERM4x6NBdk5qWhuCIGJz1C8XvZ3yx5eB5LN99EvO3HVEtQWNX7VPdXBJspCus3zxvfDt9E1btO6Wv4RYZGxQaGgy/a1fg53s1W+XvdxUB/r5qTAfLWhUcHIw07fNGRJQVBh0iIiKyLAYdIiIisiwGHSIiIrIsBh0iIiKyLAYdIiIisiwGHSIiIrIsBh0iIiKyKOD/Ab7dSjo0GnhDAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous file, ZeMA testbed dataset was imported and visualized. It contains data from these 11 sensors available for estimation of remaining lifetime.\n",
    "![sensors.PNG](attachment:sensors.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Machine learning methods thet are applied at ZeMA are separated into steps:\n",
    "\n",
    "1) Feature Extraction\n",
    "\n",
    "2) Feature Selection\n",
    "\n",
    "3) Classification\n",
    "\n",
    "4) Evaluation."
   ]
  },
  {
   "attachments": {
    "workflow.PNG": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHACAYAAABKwtdzAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAFiUAABYlAUlSJPAAAHlESURBVHhe7b1vqL1fetY3r6WpFvHFvJgxvmgKDaEx5pUlMgwpRZ030YEkQygdkqEGYtrRdsowNClmSiR0dCDI2E76YqyFoakEGyIRAtNUpamBaG2b2o6F2iaWqqnamGj8029/nzO9fnPN/bvX8+ecvdfaZ+/rhg977/XnXvdaz/mu+/qu5zlnv+tNLBaLxWKx2INZBFAsFovFYrGHswigWCwWi8ViD2cRQLFYLBaLxR7OIoBisVgsFos9nEUAxWKxWCwWeziLAIrFYrFYLPZwFgEUi8VisVjs4SwCKBaLxWKx2MNZBFAsFovFYrGHswigWCwWi8ViD2cRQLFYLBaLxR7OIoBisVgsFos9nEUAxWKxWCwWeziLAIrFYrFYLPZwFgEUi8VisVjs4SwCKBaLxWKx2MNZBFAsFovFYrGHswigWCwWi8ViD2cRQLFYLBaLxR7OIoBisVgsFos9nEUAxWKxWCwWeziLAIrFYrFYLPZwdpMC6L/4d78qhHAjrLIf+z3vCiHcCPdoEUAhhE1WWbcJhxDWcI8WARRC2GSVdZtwCGEN92gRQCGETVZZtwmHENZwjxYBFELYZJV1m3AIYQ33aBFAIYRNVlm3CYcQ1nCPFgEUQthklXWbcAhhDfdoEUAhhE1WWbcJhxDWcI8WARRC2GSVdZtwCGEN92gRQCGETVZZtwmHENZwjxYBFELYZJV1m3AIYQ33aBFAIYRNVlm3CYcQ1nCPFgEUQthklXWbcAhhDfdoEUAhhE1WWbcJhxDWcI8WARRC2GSVdZtwCGEN92gRQCGETVZZtwmHENZwjxYBFELYZJV1m3AIYQ33aBFAIYRNVlm3CYcQ1nCPFgEUQthklXWbcAhhDfdoEUAhhE1WWbcJhxDWcI8WARRC2GSVdZtwCGEN92gRQCGETVZZtwmHENZwjxYBFELYZJV1m3AIYQ33aBFAIYRNVlm3CYcQ1nCPFgEUQthklXWbcAhhDfdoEUAhhE1WWbcJhxDWcI8WARRC2GSVdZtwCGEN92gRQCGETVZZtwmHENZwjxYBFELYZJV1m3AIYQ33aBFAIYRNVlm3CYcQ1nCPFgEUQthklXWbcAhhDfdoEUAhhE1WWbcJhxDWcI8WARRC2GSVdZtwCGEN92gRQCGETVZZtwmHENZwjxYBFELYZJV1m3AIYQ33aBFAIYRNVlm3CYcQ1nCPFgEUQthklXWbcAhhDfdoEUAhXIif+c/+9Se6utfMKus24fDY/L3/7b97+tn4i5/61rZ+NsRxK7Fcm3u0hxZAJKst+7t/479v+4VQ8Z+lexNBq6zbhGeiZDuyR0l8t8Qv/43/5Wnt/8p/8r1t/Uy4/rJH+Fm4R3toAfSX/8t/52m8f/KPfvXNL/+tv/YO/sb/+BNtv0uj5BnB9WW4JtDVreKv/9znn64Tr7Xuv/rj/+rbMfO+1r9mVlm3Cc9EyfZX//b/8fS+8l9/4pvafpfmf//p//QpDl67+keCdcdmCaCttef6/5N/+CtPzPpZWMk9WgTQW4bY6epncStx3BKyrm4VCGJsljC+FVZZtwnPZHayHfF//sU//RQHr139IzH7mmTtv8w9WgTQWxYBdHvIurpVRADNtW4TnkkE0O0RAbSOe7QIoLfsqPDgVhW3qbjNgfH6S3/9L775M3/4X3xHW8pIlL/6d3/hqS3G+7/2F/74V7QbmcfEe6x7tkRJmbl4ObHRjzjU/9d+5Ze+4vbMmfls0c3z53/qD39FG41R49TtP2Ljs2LtTH10KwpfrKd8uzCh3H3hn/rR3Gg/ulb6OenMxyQO0Gen+ice1rq7XcY1wdSPthi+/68vfuH09Xkpq6zbhGdyNtn+1f/8D73dB/u1/+dvDxMnt0x+6X/+b57ayHjmyMfi/cjcL8ZtGH12FI+X6bYO/v/XP/PDT32xGiufuf0n4z1z9DZb8FyMP0fFOHyu7f7sv/Gb27XonqvZuiZn1h+Ye50fZdQdXXvdAtNnp/onHubZ3S7DB7GzFvjXNaHPrdz6vEeLAHrLSJRdvUNCx0hCJC4SnxIVZZ6UeK+kRRvakrhoh3nS5D3+MCVpUPIFJfIqHoC2XR1GwlXSxQfjq92Z+Wwh/7zig3lq7nxWO+aDUdf1V1y005ww3gv1Ub2vG350HTUPxQSKqbvW8oNRr/XGiAuRgg+V8yq/vu4y9w2aIyb/+lngtQpbjePzUxnG/Lz9tVll3SY8k61kW/mbf+WnntqSsHjv4gE/3hZ/SpwkRNryKtN4JErqFAevfAaPSabPjvp6Gf0xjUnMxOpxKnZeac+caIfx2f11+PMxWg+JnK4dhuhRO5XVtdd8avmZ9Qdfb+rlF8P3S9deY2PyrznxWsUdRh9fd/dxCyLoHi0C6C0jCSkxOTrFQAzQpktWOo0gWamMhEnSq23raYdQHIzp5YJyzJOtIAl3dbIu5rPzGaGxa1sXgH7CoXnQj89bY8lqOWhcrOvL2tcTqFFMLgR9LWiPmPMyjav4KzIvUx9i8nHxT+xYve5aJ8wfuNbPD+a+rs0q6zbhmfzy/5/8SES8d/wkg4SodvwP3n3oBEQnC8B7ymtbncyQnL2cpIvx6uVC1tVpDl4mf1gdy+trHfFKBHWnGI58+Lw7FF9tp9+wqoJJ7V2EnF1/ToqwKkToi3jysuesvfoQj68T/iW8mIf3kRGTz00/E3UdVnCPFgG0YUo+Or0gIVYfgJFAu7qKzMuuKYC6PpeajwQFybzWKS5P4CRt/GLExXvo+stqOcj3aL1GdOuo0yI/cRuhcXnt6mVepjXqrgPzlnVC0ddOKN7O37VYZd0mPBMl285IVGqnpNbdHlJydsE0Qm1rcnxOEhaag5fJXx1HSORUMQHqu3cisRczIA4whEJXr9hdkKjMRcLZ9e9E0Yi9eci8TOvnMQrWVObiSFZPhkAnR7V8NvdoEUBv2V4iRSioHcmvIqv9SHCcMOy13YujS9xCPmudzMvEc+fjePLufChm3ns/PrvVkxoh6+rko/qusCa0EZ0YUVk9CevYG1emz4gabEtMdtf2yPU+ItguxSrrNuGZdMm2gwSOkYhJlM7of/xAAiQJ77WlDuPVy4Wsq9McvGzLnydo6ivyx/va12FuStz06cSJTmIQDHUc6ISExveys+svv53YqOAD47Wrl+mzRJ0L5Eo3B5m3E2rvgmkF92gRQG8ZCaerF0pIW1aTHIlKpx2dedu9OI4kxFon8zLxnPlUFPOe1VMMF05bY8i6uj0hgqiSsOnM10rm/UfsjSvT573rCt21PXK9RzFcg1XWbcIz6RJVxxHzEwi/FdKZJ2s4m4QdzcHLtvwx1yN25JkUErZOWzCEhwshxbFneydAR8zXX6bPW5xde61fvYbO1hy8nejar+Ae7SZn1W3C1+BIggKdmBxNOnqupHsWReZle3EcSYi1TuZl4ux8OlzIdPUjNLYEyigGWVenOXd9dZsNOCUhTtV16yiReuSZmq1xQabPemZnS+h1MR253qMYrsEq6zbhmRxNPDqBOJqg9MAuSdmT+yh5nk3CjubgZVv+/ASo1j0XhBBj6URIIkgnQHW+W3TX5Oz6K44jJypn157riV3jBOjo/K7FPVoE0Fs2Eh6CUwyM5N3VV/ScRnd7R+Zle3EoIXb+Rs+EyLxMnJ3PCIkHFxlbaJ56IFj9u9tPsloOWyJga25bYuNazwBtCSwXkXtCTezFcA1WWbcJz+Ro4tEpR3ebp0MJuD5jwzgY43r50SRc/fFZt3u8fM/fKL6Xwu0+TKcxEgujZ4A6umtydv3l41rPAG0JLBeYvr4ybyuO/hxem3u0CKC3bCQ8hD/L0SVr6j1Z6bd7apLSw8eYl+ukgFMRLxc6NSFOJUpeNQ5Wk6XMy8TZ+YzQ+N1vYoGLCuKtz9tIrCCI1E5sCYctEaA1rteUMeXT56b2xOZjES/+fX0U72i+Mi/TGjFHFzm8l9Cp/iKAvmTdJjyTXz6YeJTYSeSdaKC/J8Pu2RYYPQM0+u0woRMQTpZUxngqx7z9XlJXHKPxjgiHDq2Tx6m16HyyllXUdNfk7PqrPWN7OX1ZEz+V21t7mZdp/Wo8vFf81Z/My8TRn8Nrc48WAfSW7QkgcPFCwiIBAYkN84Qkv5jaKqkpCbtvkDjglbYeE4lZ/VSvzxq/JkuZlzln5jOCJK72xCUf+KvzlIirJzPqj7jwcnxg+NF8VccYWBcjMWls+ikeTOV1rRQDprWXeVsJR4w+4DHI9Bl8jbDqnzoXRqD6Gidszf1arLJuE57JmcSjpMf//vV3aHjtxI4EiLelnU4OqgAiSctIquDiRbeSMOoUN/40vvvT+O7DIVHjB6M/7YA5Ksaun0MMoPn5+rjA4L3PW2PpVKeuheZWr8mZ9QfND1OsMm+7t/YyfQZfP6z6p86FEci8TIzmPJt7tIcWQDp5IeF09RUSEreclEgx+tbkDdyuqomPMt2yqu2JRe3xXxMc9Z448UOZTiVqssQHeFnlzHxGkLyJ1eeKGEJ06PRE68w4Ndl7nZfTTmuF8V51mvMoTnx6X2IjxtFaaQ4SoRjr4CdYgjK149VvS47WvPOvmOp6gGLX+jn0wc5co5eyyrpNeCZKwkcTDycLv2yJjmSMj+7WDElUyVntJAbwUdvjW+15rT757P4QBCRZzcHb6lRj60Fm+hKjJ3L849cFzIja1+dY2yIy8Kv4MfoiYKpQ0Hw6P2fWX/PzMemLj9p2a+0ZA7w9dP6ZE2V1TjDyA5rTkXW/JvdoDy2AQgj7rLJuEw4hrOEeLQIohLDJKus24RDCGu7RIoBCCJussm4TDiGs4R4tAiiEsMkq6zbhEMIa7tEigEIIm6yybhMOIazhHi0CKISwySrrNuEQwhru0SKAQgibrLJuEw4hrOEeLQIohLDJKus24RDCGu7RIoBCCJussm4TDiGs4R4tAiiEsMkq6zbhEMIa7tEigEIIm6yybhMOIazhHi0CKISwySrrNuEQwhru0SKAQgibrLJuEw4hrOEeLQIohLDJKus24RDCGu7RIoBCCJussm4TDiGs4R4tAiiEsMkq6zbhEMIa7tEigEIIm6yybhMOIazhHi0CKISwySrrNuEQwhru0SKAQgibrLJuEw4hrOEeLQIohLDJKus24RDCGu7RIoBCCJussm4TDiGs4R4tAiiEsMkq6zbhEMIa7tEigEIIm6yybhMOIazhHi0CKISwySrrNuEQwhru0SKAQgibrLJuEw4hrOEeLQIohLDJKus24RDCGu7RIoBCCJussm4TDiGs4R4tAiiEsMkq6zbhEMIa7tEigEIIm6yybhMOIazhHi0CKISwySrrNuEQwhru0SKAQgibrLJuEw4hrOEeLQIohLDJKus24RDCGu7RIoBCCJussm4TDiGs4R4tAiiEsMkq6zbhEMIa7tEigEIIm6yybhMOIazhHi0CKISwySrrNuEQwhru0SKAQgibrLJuEw4hrOEeLQIohLDJKus24RDCGu7RbnJWv/iXPhdCuBFW2f/95/7tEMKNcI92swLob/5PPxZCWMxqAfT3f+7fDyEsJgJoov3SX/uzb/7WX/3TT68hhDXo3+Aq+9X/4VNvfuUv/aGn1xDCGvRv8B4tAiiE0BIBFEKIAJpsvvmGENYQARRCiACabL75hhDWEAEUQogAmmy++YYQ1hABFEKIAJpsvvmGENYQARRCiACabL75hhDWEAEUQogAmmy++YYQ1hABFEKIAJpsvvmGENYQARRCiACabL75hhDWEAEUQogAmmy++YYQ1hABFEKIAJpsvvmGENYQARRCiACabL75hhDWEAEUQogAmmy++T4CP/3jn3miq3vtvOtd73rz/m/6xrbuHmGuzLmre21EAK3jZ/7UH3yiqzvKP/fP/ronurrV/Cv/8r/w9O+kqwu3RQTQZPPNdxb8YzzKT3z+U62P54Dwkd97FEHMKwLodfKIAuin/sT3PF2/3/Ke39jWww9/3wef2nzwd/7Wth7U5pMf/d1t/RYIH/rCS0SQfHR1q4kAej1EAE0233xnQeJyvvq97376B8prrbukUPm5L3zuzW/49V/1BO+7Nq8Z1pA16+oq3/KB9z21nykEtfZd3XNgrsyhqzvKp37ge5988NrVz+JRT4BYe/j5n/xEW6/kDb/w5/v4EEfUf+6HvqOt34JxdXoziuEIirGr65Dw2hJ2lyIC6PUQATTZfPNdxfd/7Duf/oHy2tWHY7CGRwWQxMMlT9j2YDzo6p7DJQTQrfzsPaoAknjhFKfWIXj0MwM/9pnvekcbQLxQPxJIM1CMXV2HTr8QJ139JYkAej1EAE0233xXEQF0GVjDCKBzRAB9yXzznQmnNqx/dxKC4KHuG772PU+vH/nW3/6ONjOFxBbEAF1dRwRQ6IgAmmy++Z7h8z/yyTdf/3Vf8/QPi9saH/3ub3/zxZ/90bbtHltJSM/tUMd73S7zRE8S55YOcSieD3/oA+1tru42jG6D4AeUWAG/Z26XnY2FsVg35qc+zHHrlgxttQ704TM++Ozr0kFb2nVUMURbXWPgPdfd23hbxVTb+npWqp8O/1kDraf81vbU0cbj4Zr4/HivugpzOePrErxWAcStHJ3iAO/P3EqiLf26B4g/9pFvfqrbelaI536oqydInAYhmOhDPRBb95yPboHV8jq3CnGprcp4T0wal9cam9pWqhjS+Drh4pU5jU66EIwSi0Bb1jcC6PUQATTZfPM9CglJ/8gcklTXfg8lZU88QolKokKigTLqJV5IUIgwfKietlWUUa6+QuOT6NQPH7zyGd9HxN1zYmHNlNz9PXQiSDECvr2vymofh/UkLiV0/PEZXKS5X+qYj/rw2X1Sp7HlS++p/+ynP/70njbqL9xPB33Vj1h8PRWPt2d+1IHm5msm4cJcqZM/j11tjvq6BK9RAOn5Ga2HoOyMCFLSdkGhcgkTtakCRsndyz0uBASCBDGgsjoOZeBlEl30QYi5qCEWPvsc5UNx8qr34CJI8VCOTz6DP8Ok0y/Gpy31EmOUVRGkkzRgXK0LbRW3tw+3SQTQZPPN9yhKjh2jE4ItSCz05bXWkWTkm3GriFBi8jJAMNGHBOrl8uVlGh/wpzF4VZKtfjqeGwsJ1pOpCylvr7WgvT+8rERNHYnc+4xQ4u+SuNaD+Xi5r4eLJT7XWDtoB11dB+NpXr52jO0/g96Hdqx3/TnRmtY5aa7ddTvr6yW8RgGkE5oO6ro+HfJDkleZnv8h6XsbFxJqQ4JXGSj514eiEUlde8rAyyQ2XCwxHoICvC3IB3Xeh3gpr2NKYNVTH/BxquCTP0RRbU+5zxmB5iJM5eF2iQCabL75HkX/oDq6RLLHVhLypF8T0RYjn4rTy9SWZOfloGTXxXaUvVi638RS4vcyEi5l3cmQbhVeQgBJ5HTrrbl4DHw+cn1oB11dBwKE9t11YSzFWes69HNU12d0bbYY+XoJr1EASWh0dIl9hIQJyVplOgGR4OkEg9q42CLpV1+OYnZhwWfwdnyuogW6/iAftRwkTrxsSwDpNGckIqnDpz6rvcSigzjKCdDrIQJosvnmexQlno5rnQDtJRsEAEkZH6AEX30qTi87Mv6Z/+2/JBahPn7SorJOtAB1e+skRr4QFpSD4ne6uegWGCKIco/Zkd+urgNfdSxHsXR1xICAoi9IPNb12RsDjvp6Ca9RAOk2TseZEyCQSCBh81m+XVBUIaFTIYSQyiSKSPq6teRIDPgpDZ/dr8qeI4C8TKiP3zLbEkCaF3XdHOpYKuPV/QiN39WF2yICaLL55nuUFc8AjZINYkOJsKP6VLmXvWR85xKxCPlxgaLbPtcUQJrvHvUUis8ujBEJ9URIdV62hYRVd11Ac/AyxpRA6ajrs3Xtz/p6Ca9RAPmzNg5lnuyPIMGjWziIDz/lAN2WkuDR7R1v4wJhi70ToE7oHLkFVstBvlx0bQkgtd/CY5BgigB6/UQATTbffM/gv5nD//xJVjXhHWUrCe0JEMVAXz95GPmkDLzsJeM7l4hFKLm7QOnKHOqOJuWRLz8B8vKj8HMh31UQn/W7dV1A43iZRBO3zRCkKh9dx60xzvp6Ca9RAAECQcIEeH9W/IBu4yCEdBvLn3PxNiR8tam3fXQC1AmLEbQHL9NYCA3GQ1zo9Egizel8iLMCaE/QVCT6IoBePxFAk80331VsJaGtZENSoq4mWhj5pAy87LnjO5eKRSi5u0AhEVNWT19A4x9Nyp1/oeePnitoQadBLhz4DN5uC+ZJ+zPPAI1iH13HrWt/1tdLeK0C6FJwusKaIjIkPqrQkOjxNv5QNCDIKB89A9RBe9BnnfTgw8Udn/12m1N9OGcFkB50RgjVug61r2IQ8gzQ6yICaLL55ruKrSS0lWw4ZaGOROVJivd+GuN9KAMve+74zqViEZ1A0UPBjOEnTJy6KFkfTcoSU/isdbrtw2utg65PReLE10IxeuxbaE3Bny1DVGlNwftoXF830Jzq+khkdXM96+slPLoAAt3S0mt3kqQ6iYquzdZJDYKgihjagj5LRHWCYkT14XQCSGMQq7cFCT1EmN+C83r3pfbgc6Ov1gtUHm6XCKDJ5pvvKl4iQCQUSFb057YFiVbJ9ojoeMn4ziViEfJVk68nftooSav9kThBYsr9aCwXbZoLkPg1F/mhD20QVGqnmKqokHDAB2Py6vUdug0FxORxaRxvz/gag76KR3HX9XGRJf/0eY6vlxAB9OWTDOiEAej2EJDcuzYkfsQDbRAf3BoCnebUUxf58zIfx6EvdQgpb696LxOdAAIJNV5p43HphAu4Fag5SNDU2111XdQO3xrH24fbJAJosvnmuwr9L1yJx9Gtne42CJCsPTGToPgs4VJ90g68TON3t5bkZzS+c4lYBMmVPn4LCbox5JeyM0lZAg0/JH4/mWEc/EpwaCzG9pjoIzGjdi4iHHzq5AmOrCngi7HpwzjE4L722tOOmEfrgxhUe179tOmsr+cSAfTlUxEY3f7xNlUEOJyKIByU/AFR0IkXxBLoM/USELSX+OC9ynnd8uFIANXTHD+hoW+dD4IJ0UYdbQBf9bafoL/mSx/mz1wk/Lo+4baIAJpsvvmGENYQAXQ76CRqJLAkWuqJTggvJQJosvnmG0JYQwTQ7YDwQeB0DzxzoqJTlu75oxBeQgTQZPPNN4Swhgig24GTHQQOcPsIQQTcUtLtqNEtuhBeQgTQZPPNN4Swhgig26J7/ob3lI1+FT6ElxIBNNl88w0hrCECKIQQATTZfPMNIawhAiiEEAE02XzzDSGsIQIohBABNNl88w0hrCECKIQQATTZfPMNIawhAiiEEAE02XzzDSGsIQIovGb4+0j6LbmuPhwjAmiy+eYbQlhDBNDrR1/TcU0RgNDgV/R57epXob+e3X0B7TVhLba+TsTh7zjR7tbWzokAmmy++YYQ1hAB9PrRH1D0LzW9NIira4us58Bfx0aAdHUCkUSbvXYIlPo9bqxp/ftL/CVufFHH2kvgdCKMvtRtCaRbIAJosvnmG0JYQwTQ6+dSAmhLJPAXqBnjlv4StU6+ECBdPYJGwk107UCihjYSNy6EXNzoK0v8K0nUR5+B8bvyWyQCaLL55htCWEME0OvnUgJIyb6ru0UkyrqTF9ZEgkbttubG2lFfxRS+KceXbmFJVHk7+tf15zP9XsN3t0UATTbffEMIa4gAev08qgDSMzidwGAtOMFhbfi8NTf6U+cix9E4ElpHBNCqZ5OeSwTQZPPN9975R7/jG9+6Cu968/c+94Nt/aX55T/27735u3/6j7V1r4lf+cTve/NP3/2bntbu//1nft2bv/+Df6Btdyl+9d/6157G4rWr72CdWe+u7jXw6AKIBKUEB7wfnSj47RSSJZ+VYB1PkDp9oL2eJdGtG26l8F63WqqAoT1l1AHtumdJtgQQc+l8eKL3+oraKKHzqjKBgODkhDmqHz67dfS5q5/6sPbdeo6gD2N2dfjxOWoMbyOIkzquW1evW146HdJnF176eeA9c/TPr4EIoMnmm+9MECFv/fS++ce/7Wvb+mvAWIw5I1EqiSMYuvrXAmulefyD7/y9b/7ht/2uJ0HUtb0UjMOYvHb1HcRHnzOi6ZZ4ZAHkyZeE7ULAE7HfQiFJ0270fAjIj/zTx8WLBAsJkkQJ6iMfGpN+vKcvfiijrdq5v1ouIUY/+oPi9rbEr6QOagtqo3ovAz3kC/jGr+IEiQZR5857+vh6IiC8T8doziPku6sbzU1ojhpLJ0Z8Jg5dKwlc5s/cutOkWyUCaLL55jsTJdZ7FUCckjDWP/nnf3Nb/1pA8DCPa4se5zkCiHWmTz2dmn3q91weVQApqZGoPOGStEhoKlOipR3v1Q50KlKTHYmRcuhOQeQTSJY1Uaq+q5Oocb9qX8UA/ZWUBf4kNvwEAygDLxOdSMCXREw9GSIm1fm61bl7DJpbFU0dozmP0JhdnYTqSAB1Y7GuWkdedT20Tj7n10AE0GTzzXcm9y6A7oUVa/YcATTitVzzRxVAnZDo2EuOJPHqRwJoL6FW4SQ0ZhUvoL7EX8uOigHFV5M0ZeBlQond58Sct8bVyYgLGsWKaKhz1+2xI/NQPIzR1VdoC13d0eu1F5faKaYqkrrbh7dCBNBk8813JhFAr4MIoDk8qgBSYtq73bLXrhMGI4Eh9hKqRBXCAb+OxJH33fNHvfvQnGp8lIGXCfpRx6vKFMsosXdx7cVKHfF1dU4Xzxa0ha5OYnjkay9m0Mka1473uk2Gb/pvidpbIAJosvnmewZuNei2A89fcKvk7/y5P9m27RgJID07Qz0ogQG3M0YPFdNWtzuAh3ZJoB5TlwwVB31V5lBXYwR868FgYC389ouecer8ElPXv8YrWF9iUD8970L/s8+8sH6/9oH3ve0D8F1vcTGO6h36eTtHfbqYdG26W2mq03WRH15ZR7+uvO9+BtRGt7pG8YNff7XVzzLUazmDRxVAJCPo6py9dl0iJlFSVgWGOCIC9jhyAuQnEB01PpV7mbjUPI/OvatzFM8lToC6uTl7MYMEjoSyTr/8Fh+f/brdEhFAk80336Po+ZYKiaNr3yHhUcWFEhdJmlclf0/6VSSQVKkDYvD2+FM7yinrBFCNQ3R1ei6GcvyD3qvNyC/JW7FpbsBnIP46P5VDfQ9HRZBfN9aRcd0Pa662xM98JNKo4/PWWPJf5wzdGIJ1AH1mHNoqNq2T1q37GdAa6tqO4gcXUL6m1HFt1YfPandtHlUA6fmU+hxMRe1WnACN+lc6f8yL2IFbVX67aRQfZeBlopunTk4ufQIEXZ2DuNvyU9nyq7mNxInqR2JLsXQ/A96O63o03tlEAE0233yP4kmzcvR/ziOBQNKRL5KWEh2vSkx+iuCCop4ukKy9rCZJGMUhujrKiMXLKiO/isHnBsxD61oTL2XAPD125kf5XizAWFqnKmLwqTr3D92abSE/PjetBVDv7SWaXBid/RmAUZxb8WscHxt8nO606Ro8qgBSgtp7BkhJftROYsVvbcj3SMDsiQCNefR2SedPD2h3SXsUH2XgZaITQBqD0w9vK3QK4jFcSgDt+als+dWzR4jFrr67xgJhSb8ah9a4lh2NdzYRQJPNN9+jvPUTNeTo/5qVFKtAUFLiloaXgxK+j6GymsQ6umQ4ikN0dZSRyD3JVzq/JNOtvhIDVdBQBt1vMklw1PKKTslG89SpVl3Hbs220K0oF8J+Ysarz0N1LmjO/gzAKM6t+CVyumuhGKpYvBaPKoD0AC+3iPwUiIRGkteJj9qR5OopkBJ8fWZFya8KDLGXvDUmiddPbgT9PebOn3zUMZgDc+niU3l3KtYJINpRBlUg6lSk+tubu/p0dQ7rcrQt7LXVrcIq5jTveo0F82Dd6prpZ8OvX+f/VogAmmy++R5FiaPjUidANcGB+niS5j1lRxJVlwxHcYiuTkkb8UGc3SlB51cipEvsgnrYKxOa095Jxd46jdZB/n3NttAc/RpxsgUSeKyf6joRcvZnAEZxjsoZj3KNU1E/3nu/a/GoAggQGCQlIJGBPrs40PMdQB8lPT53wkh+qsAQeyIANCb+SaYkYl6VqN13508nEyqnv/usPsDrNUfVdQIIJLSA2OinGKEKo725q19XV9H16wQbp1OMIeTXy7y94oI6D9ahXmPQCVidI8gfp3m819rWNb8VIoAmm2++R1Eiq5DkuvYdo4RLwqG8Szxdn1GC6+jajuIQozqEhAtBErIn8c7v1tyE/O2ViaPz32s3Woej/oWEhW516TOiR+/1c6ITsSoIz/4MwCjOUbn87JEToOuj0x5P2CS+LqFRpoQL9CGpdclXt7C6pAmUUz965kQwJvFoTBIxfeqtmJE/yhULED/zVeKuyZj18PbuT3149T6AH++nOLtkvzd39e/qKjpl6a6XBMcWtU9dL+IYXWMJTOprnSAu/WzxevSW5goigCabb75nQATpmRWSnRJc17ZjlMjOJj/dcqnPhHR0yXAUh9iqA9ZBfl0Adn5JppRd4wTI59Shdbr2CRBoLG516URIJ4OqQ/xoPeq1myGAJMbAy1fxyAIovG50m20kpsJxIoAmm2++MxklsrPJT+3rLZGOLhnKZ3d6RQKv443QaZCeb+lilT/aqszRyVqNhTLwMtHNqUNiY7ROuq3nt6fgqH9HogdfjMd71bnoke8qnGcIINDzU2eE+7WIAAqvGZ2w+LM24TwRQJPNN9+ZjBLZ2eSn2yhQkxxJ1k8XumToJwE6pQDaKEHWGDvqsyyj+aldFSKII9XVExHKwMtEN6cOX6fqX8IL6rNER/07Ggshx5x8DVQnv3V94JICaOuEUOJsJAq7PtciAii8ZvRsUncbLBwnAmiy+eY7k1Eie07y0+kFkHSpl3hxP6Mk6f1pgw/5Upna0pekTmLFN3SiZhQrQkex8erjVR9CdbUcRnPqIKHLl4SJYocu4Z/x7/ic6m03H7O7Jfecn4FRnD5nzVdtEKuKk3JdT66BrpH7uiYRQOE1o2dx6kPN4RwRQJPNN9+ZjG4v6RZJlxiV/BAftY5E50mX9yQzv73hz6Z4X6Ctkp6SIeUSKWrHCQaf1dbHUhvQ/LpY8eFJFvA5OnGgHXR19KN/N6cO1lDrAPjlcxUOYmvNttB17PqyVhq7njjBc34GtuJE4GqtuVY+Jj8fxOM/O1x/rs/ZOb+ECKAQQgTQZPPNN4SwhgigEEIE0GTzzTeEsIYIoBBCBNBk8803hLCGCKAQQgTQZPPNN4SwhgigEEIE0GTzzTeEsIYIoBBCBNBk8803hLCGCKAQQgTQZPPNN4Swhgigc9Tvy8rfnwn3QATQZPPNN4Swhgigc+jLSRFB/BXifA9VuAcigCabb74hhDVEAB2H0x/ED98/1dXfA0e/qT7cFxFAk8033xDCGiKAjvNTf+J77v621yPMMbyTCKDJ5pvvPcJXJ3RftxDCLREBdJwIoHCvRABNNt987w19p9Toe7RCuBUeXQCR8P2hZr5Yk8+UezvVV374+z74Fe0q3FLq/HsbfFBH29oeIVJjEWqLT/n+yLf+9qdbdV37Oldu5fEck9qrvOJiSP15/7GPfPPTe8b9sc9811OZvp29WxeJK9p4ufukH3HxWfGpHe+9TmNuoXjqmLXev03ex4Fv+Nr3HBrrNRMBNNl88703/v4P/oEnAcQXXXb1IdwKjyyAlMCV5Ej0nvhqUkRcUK7EDIgQ9+n8/E9+4kkcAGPJB/29HeX4lX/a11jqOCRktaUfPiQkKKsiiLnIl+ZKOz7TlzaKjzKfo68D/ahXO3yprXy4T2ckgKpP/OGX98Da6TOvXrcnTFg32tU1F1oDrhWf9TNBTMQJel/73hMRQJPNN98QwhoeVQApGZMA6wmLTmSqkFAfEqK3H0HSpL0LiA61A0SAjylR4GNST2xQhZFip5/KJMS6WGjvZXtzpJz6zhdoLp1g2BNA4CdHmgsQv4sdCRVEn8pGSEiORKT74PNILN0zEUCTzTffEMIaHlUASVh0iRp0ynBGHFS2xICjdi5aBGKHOpAwIiY+IwJqe6AOwaDPnSgacVQAjea0Nec9AeTiR1AOVbwg6ij3eY5QTHW99DPg11g+XYQ+AhFAk80335no+RweUoZ//Nu+9ukz/KPf8Y3DB5cp/7UPvO/puR61p++vfOL3vaPt3/vcD77tz8v/zp/7k62Pbsx/8J2/9+kWmtrxnltrtV0IL+FRBdDoVEB0ifysAPKTF/qMbtd0YzmKlfH53N2mcagD9Vei7wRG5agAUiwVjc9rrdsTQJ1P1dVyoHxU54zEEp+r2NHaUk6cujV270QATTbffGeCsHjrJ/xJiPCKGEGESJT803f/pieh4n30TI/qae/iBF/eHmFFOe28XH1oTxz/8Nt+11MZ7bt2vKod41LGZ28bwkt4VAG0lzy7RL4nDjpIoHo2BxAzVQh1YzlVIOjzFp7stwRGZW+Oe7625iLftW7Lp+pqOVA+qqvIj9Zet7+6UzGEokSn2tz7iVAE0GTzzXcmEkCAEJHY4VUiw091KJc44vRI5YBwUZ2LmE4AqayKpYoLNC/3+PLr9eFSPKoAQiCQ3K55AuQghPClcV0EdWM5oxOgUfuKBFj3zE5lb45bYgW25iLftW7Lp+pqOVA+qqvotqEEj07FRvMArpHG55Zo1+ZeiACabL75zkQCo96eAt0e81MWxBBl9TRHcDpDvQuWLQE08iMkcuopFCj2KsRCeC6PKoD2RAEJj3oXKnvi4AhKxP7grURD90yPbt+AyvRMT9e+Q/67047K3hy3xApszUVxrxBAepZKt7x4RVh2bSsSoCOxfA9EAE0233xnIhHhIkd0pzS8p2wkOs6IHd3aQuQgrKrI4TP1QHwV/KnO+4XwXB5VAEmIkAhrYtMJS02QlxRALhAkGmo5iVpCzcWLP9PSJWXqXUy4iKoig3hcBOKPdiNxsCeAdGuJ2BhX5RI/sEIAgU595POogJQAuufbYBFAk80335kgHt76aW5FRCdcJDqo87bijABC4HBipNtmvHoc6rdHToDCpXhUAQRKiMCJD4mRxM3nTlycFUC0AxItSV/jVd8SQDp1or7G4mICJKQAv/gA+agiQ6IO6lxrWyV8XjUH1fGeupEAgjoP+RvFtuVTdbUcKB/VdUicibqmjE+siE6tp2J3AXqPRABNNt98Z3JWAHGrjLJLnAAJhBAnQLrdhShSOZ+h9gnhGjyyAAKEhBIzkPBIdjU5gk5H/PbVFiRQ940goG8VVrSjnleStAuIUSxAwsafhAwgGEa/7VXnynvGrCcbxOcx0EZ1OpGqc3DwJ7GncZiX1s/9wZbPLQFEbNDVjZCgwW+tY50p9/XUGtW290YE0GTzzXcmZwWQngsaPbysZ4AkYmBPAAkeZqYdJ0Eq0+lQvT0WwjV4dAF0C5BgSbaPkGjDbRIBNNl8853JWQEkkQL1b/74r8f7b2adFUA8G6QyPXM0Elw1hhBeQgTQeiKAwmoigCabb74zOSuAQL8JBty2ol63r6CKks4P4yF0EDa892eBvD8nP/6wNG2BfmqvtiG8lAig9UQAhdVEAE0233xnolta3TM9Ei7dr8hTp+eBADHCZ8pr2+4vQXNahCBSf+Bz99edEUESTGqLGEIE4bu2D+G5RACtR78hNXp2J4RrEwE02XzzDSGsIQIohBABNNl88w0hrCECKIQQATTZfPMNIawhAiiEEAE02XzzDSGsIQIohBABNNl88w0hrCECKIQQATTZfPMNIawhAig4+kvIXV24XyKAJptvviGENUQABaGvqjjzvVf6G0Zb3w0WxvjXjqz8W1ARQJPNN98QwhoigILQ93dtfc9XJQLo+fCdafreMX1h7tFvqL80EUCTzTffEMIaHl0A+TeE39MfItSpwpnbWbTly0K7uhERQM9HP3tnTtyuRQTQZPPNN4SwhkcXQPyPWwKIb0vv2rw2OFnQnI4Kmucm4wig56O147Wrn0kE0GTzzTeEsIZHF0A89CuxcObWzy2DGNGcPvg7f2vbpiIheFYERgA9nwigORYBFEJoeVQBpOddOhBFavfzP/mJp7Z6VoMTFRKW3zrz5I/goIxXBJUElguRzie33/QQMlQxxqkO49JWbRQLdfKruo4tkaI4RyKQcs0NaI8/JfHON2vkAlPx1nbAHHxNKt6PNvilj/zjm/mrjeKVP17xr7WqnG0/AgHpDzUTF348NtZK9ZXR+l+bCKDJ5ptvCGENjyqAPDFXdApCMlJCrHi5J0mVVf98PuNT/vb6gHxvJVYSsfusqF1X5+Pjh/H0Wcm+CiCdKNFeD/iqreIVrJ/8IRZoq/WjnM/un3J8uT/aqY3EKWXyJ/FGWRU1Z9uPcOFDTL5OvErcIIZ8jrzyGc4KrksRATTZfPMNIazhUQUQIHSUsLrfviF5q14PSJNkldTAhQUJTuUV+T/ik0QqnyRL1VGuJE+5+1J78PKjCbXzI5SoEQcqwy+f1c8FCu8pI946voSF32pjHShDAHhbjSvhICgDFxWC8Sjv6jROnceZ9iOInbbM2U97fJ2YT9enznsFEUCTzTffEMIaHlkAeQKvz764OKoJ0G9/kdC7PupHAiQh8nrUp9d5jEqsvLovkrfaU6dyF2dbuGCpdfLnYzgSNC6AFDNz8ragsXzdJHQ0PyGBIKEoKAMfU2hdOkEL1PlczrYfIdHZxcS1pw58jhFAcywCKEzl577wuTc/8flPtXXhtnhkAUTCV2Kq//vfqlPiAk9eJFGVe4IXWz512gBK+C5mtiCJy89InG0hUVJPKPbqQGvhiV/zRAhR70gcub/nCiAvE7oG+Kxjy5/3Pdu+Y08kgubo6yT/vHrbFUQATTbffGfy0z/+mbd/qDu+5QPva/uthJi6WAVz6vqt4jf8+q96iutTP/C9bX24HR5ZAPm/IS/3/7F3SU3JDDyheXkVOHs+dZIC8umnPA79lbBHogGOJlaJnO4ESIKK8WodaDxfB42/hYsz+ahCR+tZ11I+vEz4NRjh63+2fYfWb7RGoHF8nTTvo9fpmkQATTbffGfCyQQ/dCTp93/TN76D7//Yd7b9LgXjQlc3griI+eu/7mveES9w4tL1WwVxEu/nf+STbX24HR5VAClpQU1cW3X1VMbrVNYlzOf4VIKE7nZShxIteLLdwsVZrVPcNWbRCSCdAB0d3+evUyPNg8+1vdrWctCJzlFRcbZ9BwINH1tCSfPxNdHavWTsSxEBNNl8852JBBDCoau/NowNXd0IYqVPbiuFS/OoAshvOZEEvU5JH2rid4HhJyaexDuxsOXTT3/cpxIkHBVAJGH1Qdh0bTrUp5bvJXfF7oldZWdEG/4RO4qfZ2pGwoB66Op0Xes1HXG2/QjFXU/kwAWmX5MIoDkWAWREAIXwZR5VALm40O0Y3WpxMUNi4zO4UAE/nfBnb7pk+hyfLtIQRrSnnCTKeAiHmjzdF/OhrfptsXVqowd8fV749di9n27d4dMTvqCtx8R6MEbXtkNjdnX4pQ6f9daZ6j3Ws+1HcN26OfNeormeZkUAzbG7EkDcVtEtFm4lffS7v/3NF3/2R9u2HWcEELeWPvyhD7z56ve++6kP8DxOJ0SIgbZ6/kVj6PaURExH9VU5K4Bo588NEdMobvnuniPidmA3Lv7ox5zVnzXSXDV255Pr52tBv+62o57Voo73ugb0rW3D83lUAdQ9X+O/NaWkVVGiA39mxQWVP5TsjHx6ufskeZKYvW2lJm2JFWcUj4O4GbV1cYd/xUtsI+GkdaINvlkfXhWft3f/Dr7xU+eoei9z/NrSn7FBsVbBcbZ9B9dK7YE18utKXRV4+D3q/9pEAE0233yPQvLUD5SDIOrad5DM6bOXSGlHogeEDYmYV41ZRYFEmdoizChTu89++uNP5erPe+F+OiQY6pgdjKsxGJ++LuCIw9tv+Va8tU6+NWd8sE5qN/Kp2IiH9/h3H96WvpQjpnQd5NfbhZfxqAIISMgSGCRmT/4kKz/h4D1J209lPDGr7egkAUY+vawKCU4gSMoubBiD5NqJFcb2RIzvmng7JEIYa1Rf/TKW1qPGDcTnIoC46Ycvb8eY8inxAS60WAe15zO4j4rWlXYaH38uMF/SvoN1Jm6/VqwZZd010NqdGeNaRABNNt98j6Jk2XH0gVslViXUivwgFEi+9XSJ32yiP0JHZfLpZSMUb1c3grjoI0HjEKPa+dyq+FDc1Pmc5Lu2hy0BJF/dKU/nU7Exh7qmOjFycab2oz7hMjyyALoVlDD3kvo1IUGviAERxbiIja5eJ1OIiK4+XIYIoMnmm+9R+Icw4shJCnhi7dj71W31J8lvlY3QOF3dCAmKDkSI2umEarQWEpAuNF4igLo+0PlUbJ1Q1fqNxFzEz/WIAJoDJxg6ZdJpgE53+DkHkn3tNxPFUk9orgknL1tz51SG+u60K1yOCKDJ5pvvUfxWTuXsCdARsQI816LbV6BEXvtLXBAj7UdJW/F2dSO2RIqj9elOZUCCxgXSSwSQlzmdT60P66e1FN2anr1O4XlEAM1BiX5E94zIbBTj6DbYNWDOuu2k20XgzwtR3vUNlyMCaLL55nuUmc8AIWCUmDtqf9rzXAsnFtTzSnL3NqD+tXyLTlB07PkmHuo9ri3fal/r9sbpfKrPFt0J0N51Ci8jAmgOJHqEhT8fAtz6uaXTDT1309Vdi+45J63NLTwf8whEAE0233zPgAjSaQIiA9Fx5hbJ0cSqB3ZJyn6istefWDgB0mkMfryeMvCyPbZEiiPxdcsnQHtzELSj/d51Ci8jAiiEEAE02XzzncnRxCoxUcXV0f7cOqMdfrycMvCyPbZEitM9TOxIhPjtQvnubiHKXx2XMvAyp4tXvq51qzI8jwigEEIE0GTzzXcmRxOrTnBq8h89A1SRAKq35ySs9DdzjtAJig6ED+0Yo54C6USLeXXljCGxx6vf/qvjqtzLnC5excZ6dCd2tPU14TPt99Y5vIwIoBBCBNBk8813JkcTq27/ICYQCXxGPEjAeH/qSOx6wJf2aldPYyQsqMcHr17f0QmKES5ciEljaMwqjBAdqmd+3l4nRnVcysDLnFG8PnetKa+d2Dx6ncLLiAAKIUQATTbffGeCADiaWCV6lLT1PBDvvT+3dZT0BZ+72z2cfuh2EPC+tqmMbkWNQHRJvABzQHyMTp2Yk8eveepvB9VxmT94maN4q9gCYvOxtK51rXSdjqxPeD4RQCGECKDJ5ptvCGENEUAhrIU/AzD6Q5D6a9H8iYKu/lJEAE0233xDCGuIAAqifjWIf/3EWfQHDEdfCxK+jP4AZSdyuAbUjQTSpYgAmmy++YYQ1hABFIA/Pkii1R8jJCm/RLzo7wld++SiohOTEVVI6Ks24OgfgPQxtvrwF7W77xdjTBeX+jqQkS+JyZcI0j0igCabb74hhDVEAAVA+JxJskraJOeufpUAkpDjBIsYKvUrNyijvdgTff6Xq4H+XTud6oDHMupHG/x6meAPZdLnmn8UMgJosvnmG0JYQwRQACXmrq5DX5sxEgBK9qsEEK9dfUVxSrDsnQLp9Eftu/mrDkHTzR9BU4WYTqK672FDdFF3za8EiQCabL75hhDWEAEUgAQLXV3HvQkghAeCBba+k00nNbTv5q91gTO3EOWvCiOhE7qt2F5CBNBk8803hLCGRxRASlLd7Rv/H/6oriZXkpcSKZAkRwmYsf25EF4Zq7v1RB1+SXryj+8jt6log1+NA/io3zvmcTujW1vQtQd8qY38Ml/WRwmceEjyo0SOaOjW52jif64AIk6d3Iy+m00ihXb6GfI5g3zw6uV76Jbi6JSHNaGeGLr6lxIBNNl88w0hrOERBRDJlGQCtU4JEaFR65SE/H/2unVBe96TeJXsa3KUgPK2Go9EX5M85fhyf7Qj+Xq7ihK1xqKffIAnZ5K9RAPwHrYSLfVK9PhXHxcOmpfGVRwapxNYiluiB59a8259OuhDe167+opiYk0laoi5a+vXfySAiJPy5wgV+kFXd3ZeZ4kAmmy++YYQ1vCIAghG/6NWEgIXOhJNJHKVecKsyVn+XRQoqXs76NoCZUBS9Vi2IA4l4frQLPGqjvdep7G8bIuRABAudjwOnXSAn2Qp7m6uW6dyFQkFCS6Hda7tFafWg358rjEQK+W05/No/pTBkVO6iq5NJ/S0BqNbZC8lAmiy+eYbQljDowqgLqHoBEJJ0RO36jwJ6xSk+9++EmSXdCuj/91TBlWsbIGIok9NzEInVlVMaCwv2+KoAPI1FBJ8Pi/FPUrw1CEQujpHa9nR9VecimUktrRuEqmj+VMGXnaUGouzt94vJQJosvnmG14GXy/hX9lx9BvXQ3hUAaT/0fuJjpKc6vxWSCd2dHuHOhKvo/ZdwuJ0gUSrtkp8vPd2lIGX7aFxO+EBl0rcewlZc+qSOfOkzk+8tPb007o41IH76VBbXrv6So1TJ32IJZ3E+OmU+u2tY3eKs8fWmu2t90uJAJpsvvnOQj+cRzj6vVur0ZeGInz05aL1C1hDGPGoAggkYHS7g88SPapTIkMo1RME7RVb+AkQ4yjJddSkrXIv22MricJe4vayLfYS8lYcnUjZWhfRneBUXiqAoIpInU65z7117Oa9x9aajca7FBFAk80331nwJZyOTk14rXXdF3neIgge5hDRE57DIwsg3e7gVac+ugWjOpIfwoX3JEbvL5F0NNmpPYnUnxEZJW3KwMv20O2l13AC5PPVCVBdg7OM1nJEF6fmptNBPRfk12w0f/nTz9EZEHj07U6P9PP4HL9HiACabL75roITE36oeO3qXwOINebwWk6swm3xyALIb3Xpf/lKhC56lFTrsz4SG90zQBX50wmTM0ralIGX7aFEWcWakNCoifTsWCMBIM4KoEsl+NFajhjFKdEjf/VZrtH89XOEmHHBdAT6QVd3dl5niQCabL75riICKDw6jyyAQIlOYqbWkch0cuN1oGRHffe/dpKkkqDEFv68Le/9ZEjlQFk37hYaB4jP6/QgN9TkrHIv20KCTqcklbMCyNcH395e9Z2vylmhMIpTgkxUkTsSQKDryWsXM76q0NsSyHBGbD+HCKDJ5pvvKrYEELfAVMd73S5DcKgNouNbPvC+p+dvqOP1wx/6wJuf+8LnvsIXUEffL/7sjz75VB/8fuoHvvcd7YGHmb/+677mqZ3aKlbFXsGv+1CMXs/nTjCpHe91a432eqja14Q5ul/e67Yhr17HmjBvjRNuh0cXQDoRgfq/fD0L0tXVNiRufJF4eZWw8gSoZEud2tEPKK9JmzLwsiNImGksxlU8UIURqK6WbyGfGsPFgObaCYCRSPG4dfIGI4HYId8ek1OvI2W0r3EiTBULvrwOtgSQi1r11/gjn5r76OSO9tRX4XopIoAmm2++q5CIkKhwEAjUkcgRAaDTFuoRLbxHlCAW8KF62taETzliRoLG30MVQYgOyiV6gFiAeuKjTMIMkcFn9yMRA4xFfGoP9bkhxY8v9XHRVdeE9/TRe161Llov1eHLxwq3waMLIP3PG+pzM0pysPU/b5KXJzcEDYm29iExkuAkeEhqfNY4NcFLHHnZUfCpUwOPqSZ68ZyxWDslevp6/Bq7O82RSKnrDYobf4qdte3adtSTm0qd41acErfd2Pq5oX+tA641/fznAviM8K1CRmN14lSnY50QuxQRQJPNN99VkNj5wVKCd5TsgeRdBY0Eh5eBTj6quJAvBIGfvriQ8vYSI91pkqN27hMUfx0PXKT4vOQLuoeq65ooNny4sGIN5Jc2EkE6IQq3w6MLoBBuAZ3wIJxqnU6HEE617lJEAE0233xXcUQAVZGwx8gnZdCJAAkELxsJm8qonU5xurkBAoZ6FzryNeqjNUHs1DXRvHVC5SiWesoV1hMBFMJadJo0uv2lU6pr3f6CCKDJ5pvvKkZiBZTsEQW1zkHQkNjxASMRQRl4mVAfP+3RLTBAQIxOT0YCSCcyo37d3Ee+xNaaqK7OG7qxwm0QARTCWnT7q7s9qdtf3DqrdZckAmiy+ea7iq3EvCeAEBYSDB3Vp8q9TIyEB591UgO0q4Jm1Fd9vMzp5j7yJbbWRHXdWnZjhdsgAiiEtfBc0kjg6Jmm0bNblyICaLL55ruKrcS8lexBwoS+fnIz8kkZeJnYEx6IHt1G4naZjzfqu/fcTRdnBNDjEQEUQogAmmy++a5iKzFvJXtEBXWIoFo38kkZeJnYEx6ie5Zm1Hf0MLaQgPPvDYsAejwigEIIEUCTzTffVWwl5q1kzwkMdfUBad77yZD3oQy8TJwVQEdEC8KHcmKsp0D69fjRb56N4thaE9V1a7m1zmEtEUAhhAigyeab7yq2EvNWsgeJBUQE/REViA3deqo+KQMvE53wwA/l+AG1QWC56Or6Cgkm9aOt4uuE0ZYv2FoT1XVrSdmoLqwlAiiEEAE02XzzXYX+Hk6XmHWbq/u1bkCEIDAkKBBCfB4JAdqBlwkJDxckjKvf5HL/9dfPdaurihnBSZBOpdyPP0ck9nxtrcmWANI659fgb48IoBBCBNBk8803hLCGCKAQQgTQZPPNN4SwhgigEEIE0GTzzTeEsIYIoBBCBNBk8803hLCGCKAQQgTQZPPNN4SwhgigEEIE0GTzzTeEsIYIoBBCBNBk8803hLCGCKBwi/D9WNDVhcsTATTZfPMNIazhkQXQN3zte57+PhV87oe+o21T+eDv/K2n+7xmfuwz3/U0Z8SI5s0Xd37sI9/89E3lXZ9LoLG6unB5IoAmm2++IYQ1PLIAUpKFI6cNfCO39/nkR3932+6a/Myf+oNPYyNKuvpL8pFv/e1vz/W3vOc3PgkfUNnoG8wvgcbo6sLliQCabL75hhDWEAH0rrcT/d6Jjk5/1H6FAJIIu6b4AM0RYciYtZ614hSoll8KxoauLlyeCKDJ5ptvCGENEUDvOiQquN2jNggf3t+rAPKTLk6cujbXRuN3deHyRABNNt98QwhreK0CiMTsz+Pw/uwzKZ5kucXD+5EPTjuo5+RjSwD9wp//Q0/l8gc8a0QZdd5Wt7OoY1y/5USfevKiukoVQ1obPbfDK77r+CMUB69d/QjWhlgUF2uwN2/ea618HvLh/QTj+PNb9CfWaz6TdO9EAE0233xDCGt4jQKIROcP5QrKziRB9eM9yZj3o9s6+Abeqy2v3sbj4pWE7oKApO1iQCctLlZo7+IJgaD2jCdxInEBfuuOh5apxx9tqZdQpOyICFIs+OrqOzSGxB50wgbqvEHrpDa898/ChY/WV/Hy6usVjhMBNNl88w0hrOE1CiCdxnSceS5FfXiPeOE9SbS2Q2BQR1LnM6/+WSiJIzxcaOBbidv7+K0m6l28SVDUUxj1qaICGFOCogqBH/6+D7b+OhTTGTFJ/FUwEY9EkPuq8+5Emeq9TOte14r+EobduoR9IoAmm2++IYQ1vEYBJKHRcSYBqo8+S3T4iQpIvCjpdgKIOsoQH11C18kMgkBlEgKU1T66TVTnsyWAJNRGIpC6TuBVaAdd3Vl0rYhbZZrDaK2gi0Fiyn0J/KjPGeEWvkQE0GTzzTes4+e+8Lk3P/H5T7V1r4mf/vHPvPn6r/uatzfB7//Yd7btZsGaKp7f8Ou/6s1WPB/+0Aee2n3xZ3+0rb8mr1EA6X/7Hc89AQKJFBcXStaII5V1Akjiw9tVqAd93hIzQJ0LJtjqo5Mx6oitQh3UfhW1G4mTLYjPx+xEy968QTHoswtMb+fgr44VjhEBNNl8850JiVL/uEZ0/e4VkjNz/tQPfG9b/xpAOGgeH/3ub3+D2OC1a1uhL23p+y0feF/bRvCzQxuNxSufKfd2iErq3v9N3/gkhCRwPvvpj39FO/j8j3zyqW6VYHuNAohkSCJk3RzKzvzvX/28TAlbt5AktjypSkzwulVWoR70eU8I1Paw1UcCYIszJ0BnhATiUWvX4b725g3qp89H+mj+Z+IOXyICaLL55jsTEhL/SJSgOrp+RyAZ4rsmxJVI8I2Su04pSMRd/WtAIgKh0dWPYG381Gjr2msM+Or3vvuprQshv+aIGcoRQirTz5s+g4TbS37mXsprFECAQNEtK+D92Vsf6utlEjIIH5061FOYTuzoGZtLnwB5e9jqoxOgLRF2BAmJo6dpEqTASZifHHWi5IiYoR70WbcEtwRcN1Y4RgTQZPPNdyYSQNdIOvjE9y3dUrrmfG8FCY4zpyiceNEHAaITmtEa+QlTPcXp+koIezvqq38+49eF0mxeqwC6BFwj8DKJHhKtBAXixtt0AkgJuoolodtrPE+ksksLIImwM7cBO3Q7jzU4Iiq3xr2UAALioayLyZ8Bes6tu0cnAmiy+eY7kwig++M5AkhromegttYI0TOqRxxRBxIyRwSQBFh3W2wmEUBfmWTBT5agJtROAIFuAdXftEIcqc4fsH6OANoSWi7eaNfVuxDZQg9+89r1QdBJ8Egw1XkQg0TLJQSQbkcSk18T3kto1bUPx4gAmmy++Z6BWxH+cCnPb5CEurYdZwSBEllNrIzHbRDqSHpKwB0SQ7oVRVveq7/HQVvG1GmDTidGJwRqr7HwiX+th8orPqb6+y0cUf0TD581J0cJnTpgDPWjz9lTDoSB355ibnUtGEf1lW4+jt/yk5/Rz4ROeUbPSWmuEjP6efBYtXa8Jzb/vJIIoHcKIJ3WQJdQRwLIEz6vJGUJic7XcwQQSEzxSl/vLzECjEeMoDhqzCMQFR67j+Vlaqt5U88YEiuXFEA1phpPFUbhOBFAk80336P4cxgOibJr37GX7BwSmMbwZKaHZiWM8Ml7iRoSJp9B/TQuSY/kB0qc1EtA4EMP8qqetlXk6VQCmD9tacdnxcWrkrfEEfipg8YgPpWB5uj+NT+oJxf4pVzjaX6Kib5HhaoLH3y4H14lblhbxtUceNUcLymK5b+ukag/D/q5oR99VC/RxfyYx5kYr8UjCyCSM3R1EhndSYpu+dRbY8ApC8lfiR9I0n7yI3SaM3puiLouPvpJCFBfRQ0CA581hi7eLRAT9HGRIV+c/vitKGLykzPiIy6tlQugvXkDsXdzJyb86vr4WBE/zycCaLL55nsUT4wV/x/9Fkp2StAVFzqgxE4dn0m+fCahezugDXVdotS4wDxq8pNo8jJAMNHHBQcxShBUIYKQ8jKNq/grXcy+RnUuEmo1gWudgLmojlcJpxprh/ywRn4t8CNxVeeiPt36HWFvjRR/XQvRjc/Po/rxqrmr7cjXbB5ZAIUQvkQE0GTzzfcoJI4RR5Ofkt2I7taJhBdJTO+7BHZEAJ39n3+XXCVCEATetkPjnhFAEhqjNfX1UJni7G7rKN4j12hLbLBu1IGLo26NzrC3RhqziwmOjq9xOBHicxVJrFPtc20igEIIEUCTzTffoyhZdJw9ARoluw71EUpglSMCaG9cBBiJkGQK8unJVQLlSMLcG7eLWevciUHoEn5XJhTDnmDTrSNEYlcPXbxbYx9hb4102uZjOkfG10mYTv80VwQjfnVNj/4cX4oIoBBCBNBk8833KCQHkkSFpNK179hLdiN06gEjYdAlZ7E3Lj7Vv8OT69Y4lb1xO18a09s5XcLfEgF7MYgj7bp4t8Y+wnPWyDkyvgSOfnb0XJCfZPF59oPREUAhhAigyeab7xkQQRIj/M+cRHLmltLRZOxIeOlk5DmJcm9czYkkund7p3suaMRzkrtOPC59AjSKQTAe7W7tBEhjjk7cVD86vdHPj8enPt6On4FRDNciAiiEEAE02XzzncnRZCwQVyRkQJhIqHTJUEnNk7PYGleJH9+1rkvuKtu7pQRHk7vHvCewtAae8Ls4xZk1l/hyESj8GSAXvVtjH2EvPq419d3pzCgmr2dO1bfWvZYdWaNLEgEUQogAmmy++c7kTDIGiQEJHj+lqEl6SzhsjYsf+fQkyns/GartwYULMLaPr3g5vfJ2QonY/dCfMuKpp0C6dVP9XUoA6VYR865roVir8Lu2APL19rUlJl3zkRjFZ/ezonX0OW75uRYRQCGECKDJ5pvvTJTs9L/yDiUr3bqoJzNKXvVEQMIB8INIkLDYS7KUU08fEjljECNQXpO7YgDio/+orW7d8Uo7j0HjKk4hIdL557UKoy0Rsjd3x0Wf+ihGqMIItsYewfzkW+MxL5Wxvt5eYygG2mg9uphAJ0cumoTWhJ8h3mu963W4NhFAIYQIoMnmm+9MdCKyhZK7ElxNSiQ7iYruhMQTo8SUxq2iSeCTJKi++OezEmWX3EmsSt4aj3Y1GTO22uHffekEo84Dqn/FVE8zQMm+uzXoyb7WdRA/MWqNYTQ32Bp7hPvuQODUPgjiuh6jmChjrVmvWidYX8XB6+gZomsSARRCiACabL75hhDWEAEUQogAmmy++YYQ1hABFC4FX5HRfX1FuH0igCabb74hhDVEAN0eR74r69ro+7+6ug7F3H15bLh9IoAmm2++IYQ1RADdHke+Lf3anBVA+vZ3hJCX68tQK/kC09siAmiy+eYbQlhDBNAa9M3pVTDAaxRA3PriG9prOSIHP9ThE/RN9sD7iKD1RABNNt98QwhriABagwQGYqfWvTYB9GOf+a6ntt3tLwkgXr0c4Ydo6urCfCKAJptvviGENUQAreGeBNDHPvLNT20/90Pf8Y66kQAC2q+eZ/gSEUCTzTffEMIaHlUA6aFdEvPP/+Qn3r4lBbzXrSlevY5TjtEtG05CJByA2z418UsQdEgMuQBiLMbUaQm3jDrRBLTFP+PK596zNsTst6QYi/U4I4DUtrudp/nWdYCR0NOcfR5+TSoIqW7d65z9mvNe/n186ur6sUbuB1gjvy7y04lAH1f93P/oes4kAmiy+eYbQljDowogJV8Sq5IYCUzvedUDvLz3OpJW9adTEJIn70l2EhaeYBnXkyzJkM9KjmojXx6bJ+YqBujrsdMe1J5YqiDQCYzq1Z7+Gsvbj5CPro55UcfrqI5roDKfB+W0caFRxYLEKfHTFhS7rzvUa651oox6XUPK5Evv3Y9u+QFj0UbXGojX249+1rau52wigCabb74hhDU8ugACkpfEByLBExNJS8LBk7MnLPnqRIYSdD0ZUOKtCR1GsYH81SQrf5R7DPRVcvZEThvNxWPz9qDyLbbaMiZ1VUQwR43vJyyaR10vnaJwbbycWOsJjV9DX7u6rvVadf4rvm4IZK/zOfl1fc71nE0E0GTzzfc14F8nAf51EuG+4as8uq//2GLrK0ZuiUcXQCS8mgiVtElOXg4kKuo8+amsu1WicaqvIwKoi01CgP4qI6FSRvKt7UEnFp7cdfrTzREfEhC1rqJYSexdvdYSf8QM8g2e+DWPkS+t2ZHTkm59FetonbbqhNbN19/RKZLP6+z1XEEE0GTzzffW0fc68YPKd30hfuqXZYb7RN8zxvXv6kfwXWL0m/3lpmd5dAHUJR7V1VML6E40SNiUkfQodySO6jhdghZbsQF1R8WMoB70mdj4zKu3E4qvq3P2YtU4DiKDWKtgdKGm9XMknEZrttd2L1aJF+LDh5/UCF3PevojujH2xqXOr+cKIoAmm2++t46+FX7riy1fO8/5QtFHQNee0z8v3/ty2wigY+ab70y2kpLqSIK1jrJax+c9nnMCtJUwQZ+7mCq1j5L9qI/i6+ocxbp3ArQVm1DbPfwECNEksdPh67u3roCwcX8IHj+12bpu0I2xN67G6upmEQE02XzzvXU48eEH9J5vez3CHC+JvuG++9Z4iAA6Zr75zmQrKamuS9pdQtcJ0Cgpdmwl0rMJUw9rrzgBQhxU387eOI5OgEbzrnBCw2kNcAq2J1T21tUhFvlwcccaU3bpEyDo6mYRATTZfPO9dSKAQiUC6DLmm+9MtpKS6rqk3SV0JcV6S2cLJVfGqnVnE+boAWEhYeGJfEs0ISR0ClLrOmg3antGAGkeo9OkiubAaVat69Z3b107tA46ddKY/oyPo5M1j+ns9VxBBNBk8813BZ/99Me/4qHmr37vu59ucfnDrkpyHXsPt3LrpPqv4oLbJ3q2pLavsTi0VYId+XbqXHlPGXVbc3SfxMmYPA+lsRlXMaqsWxeJqyoG3CfPVOk5K+JjjrShjrVQHe2PPFisMbtbenpAWWvQ1SnW7lYXnzuITW20Hvjxa8s8mCvzUtuVRAC9XADpGRwSd33IFfBXnyeRaKKvl8NzEqYSdU3MJG7V+Vh64BhcuNFeJ1qg8i3Unrhr3RkBBF2sgrX1WLXudZ2YA6dCNaa9de1QPLquvm41RglN8Ov9nOs5mwigyeab72xcDJCoQAmWVyVYkjuJVMmMVz7DVgIj4dFewgRIoJ5EQX5J8GrvsfC5jkPyVJ0eyNZ86OttQb5VrzGB5Lw1Rxcs1DGOj0WcaqP+3kfgq6urPnllXrwHxIvWAv+qo2wkDoWuAf1qnfyzNrUO36DPxFz9MB+/ZlovF1RaD82NdiqD+rOwigiglwsg0MOxJF7+9089r0qgnohByRuIgXZq85yE6UmfV/q6kOlOLHRaAbRVe2JR3LVPh/x0omW0XiN8HsyBfiDB6GuCKKlt/Trw6uu+ta7UMWfG0Zhag7p2fu1ogz+1hboOz7mes4kAmmy++c5EyZjE5ElUJw3U1aSpPrx6+Qglur0k7QnRTyqIRYnTx1Qypq4KI51ceBKWCCCh+6kJfRFPXrY3R+o6X0JzOSuA5NPrNBdgrhqPuDXOkYe18UtbXyutIbjQAa2XCyO1rz8To3IxurbMReV7Px8zeFQBRKLlGnS3gJSwSIK1TrdAumdASHwkOV1fkjD+R7fGEA5K1IgPnRpsxQbyXcvpT7KWTyCeTpgIT/T0oz/CQoKj61PRyUcnsrbWa4Tm4aKC9WG96gkba6VY1Y45aVwXQFvrypisla+dfNW2gF8fl3589vHEc6/nTCKAJptvvjPRKUJNxkCi1A+0J6c9cVDZEgOO2nW3YpSMEQAqk0CjztuCErKfLHSiaMRRATSa09ac9wRQFVSaC9eqCj3mQt2RExTN39dLJ2iK18dWna+XYqG9yrbKhfxv3YLr1mo2jyqAwuVAlPDzvDqJh+cTATTZfPOdBaKGf6j1f/5Ol8jPCiCJF0C0dCcm0I0lFCuoTKdC+CQWpzu9ktgbje/gg7a8dvXUQVcHW3OR71q35ZNyn4vAx6iuIrHkJzqsIegaIXpUp/XqTozqeHtxHFmPI8L02kQAhUvAiQ0/06PTrnDbRABNNt98Z3EkeXaJa08cdNBfgkVjViGylSRBfevnLfxkRGX6vMXeHPd8HUn4tW7LJ+XddcLHqK6iEz0JXn1G9Oi9TtgkOOvJ0mi8vTiOrMeZn6drEQEULoFuG3a3wcLtEwE02XzznQUChH+k1z4BchhTpzOM67fWtpLk1glQ176D8WjvY47Ym2ONpXIk4de6LZ+Ud+ICH6O6Dt1u4jroREi3xFTH+nCrivf1VGY03l4cR9bjOT9PlyYCKFwKPf/U1YXbJgJosvnmO5MtUaBTAfDbIJdIWBJB/kyIkmT3TI+StSdYJeyufYf8H7nVsjdH6qCrg625KO4VAkjryKmProHqXPQofr/uEAF0XfPNN4Swhgigyeab70yUBDlN8WTHeyUtf2YELpGwNK4LBI2HKHORgjjT8yhermReYxckWxd2ao8vL6cvc/FbchIDde6COujqwB8uVmy8at5QxcCWT/mq5XvCo8K8ac+asQ7eT3W6DmfG02kiPr1cyGcE0Lb55htCWEME0GTzzXcmJGXdSlJiU7KCTlycTVgIGnzSHuS/+vZyXpWgea9y9wsSFIyB6MA/rxJMNeFuzdXbSgwAfcDnqzp9rtBfp2uahz4rhhrblk/K8VHLR4JkC1+D+ltZWreuDrbGU1/N19vwnro6Z2BdqfP1XUUEUAghAmiy+eY7G0QIyceTnxJ+FT+g05EuQXZwy8d98x7hUn17kvR4eB3FAn7LBhAajNndfurmSl8/WRKUeQzujzHA21c4FfG4iIkyrV8VA1s+FWct18kLvmvdCI0PfuoFrA3lxOGnZGJrPOokruiPL9XRvhsPNObRn6drEgEUQogAmmy++T4qLoC6+hCuTQRQCCECaLL55vuoRACF1UQAhRAigCabb76PSgRQWE0E0Otn76sWLgF/7Zm/9VO/iiLcBxFAk80330dl6zmREGYQAfT60R8hHH3Z5iXQ915dU2SFdUQATTbffEMIa4gAev1cSgDxXV6j7/PSN77z2tWH100E0GTzzTeEsIYIoNfPpQQQPqCrC/dNBNBk8803hLCGCKDXTwRQeCkRQJPNN98QwhoeXQB97oe+4803fO173k7+vKestkNk6DkY4FYRnymvbdWO97p1RHt9U7oeWv7kR3/30/vf8p7f+PS5ChjaU0Yd0I4+3ga2BBBz6Xz4w8xeX1GbH/6+Dz595lVl4ud/8hNPX4LKHNUPn906+tzVT31Y+249R6i/1g9G14T5djHio7Z9RCKAJptvviGENTyyAPLkSzJ0IeBJVCIGSNK086RbE738yD99XLxIsJCsScigPvKhMenHe/pKqNFW7dxfLZcQox/9QXF7W+KnjnJQW1Ab1XsZINLUD9/4VZzAGnj7Onfe08fXE5HkfTrwo7VjDOLy6+nXDxST2rKmlNV2j0oE0GTzzTeEsIZHFUBK3CRQT7icFJAcVaaETbuaLHUqQt3oRGV0mqR6knD91XLVd3USNe5X7asAor9OnQT+JDbq6Qdl4GWiE0D4Yu6U15MhYlKdr1udu8eguVXR1MH8aV/XR9fEfWjMI34flQigyeabbwhhDY8qgDoh0aFThXryIUji1Y8E0KiPEnIVTkJjVvEC6kv8tawKoBGKj35eThl4megEEHPeGlenWJ0YQYTVuev22NF5dHRrcXZ9HpEIoMnmm+9Mzn6vF/B3el7yt3roy5hnvr8qhBk8qgDSKcje7Za9dp0wGAkMsZeQJaoQDvh1JI7OJHjq3YfmVOOjDLxM0I86XlWmWLrngqCLay9W6oivq+vgBAkhprl16wNaU3zTvhOej0wE0GTzzXcm+iJK/+LKLSRe4LkiiL/0TP/uyz1DWMmjCiD9m+7qnL12nTAg+VJWBYY4IgL2OHICxAmSxE5HjU/lXiYuNc+jc+/qHASMxE5H9U97TqR0W45Xn8ujEwE02XzznclZAcQ3hOtby7tvCz/CLQkgzaWrC4/HowogJcK93wJSuxUnQKP+lc4f8yJ2qCceo/goAy8T3Tx1G/HSJ0DQ1Tm6vUYMfm32/LMOrIeEIX66do9GBNBk8813JmcF0CW4JQFEHNDVhcfjUQWQRADJsKsXe88KSaz48zojgSH2krTG7J4B6uj86WHgLsGP4qMMvEx0Aqh74NiRSPEY9uZOHXR1joRpvZW1518gEGmHn67+0YgAmmy++c4kAigCKHyZRxVACBr+HXAS4KdAJFSSvE4V1I5EWU+BlODrMysjgSH2krTGRFzVBA/095g7f/JRx2AOEg81PpV3p2KdAJKIgCoQEW+q24vVUZ+uztEJTp3D6BmgimJnjbv6RyMCaLL55jsTF0A806MvJAXed7e5RreNav+KvuXdBdAXf/ZH33z4Qx948kfZ13/d17zdroM4aSOfvP/8j3xy2Par3/vuti1jq7xS/VSIuYvjs5/++Ntt9kSe6pm7lxOfx0b8jOVtQM9i6bppnj4eY3A9tLa8Mt7o1mVdL6fOQ9fN2zPWc58LuxUeVQCBTm+AhAn67InVnzWhD+0kFjphJD81OYs9EQAaE/8ILYQHr13i7/whnBQj5fR3n9UHeL3mqLpOAIGEFhAb/RQjVGG0N3f16+ocxePrw7g+Z2/LdWN+vKe92tX4HpUIoMnmm+9MSHr84CuZkyRJdkqaJDiSnfehHLxMCZ1+H/3ub3/yq+SIbz4r8aot9RqHMT2ZdolUMcof46gPn70tdfJLHeg99YgV3ms8tVH9CNZCY5Lwac9Ymof3V1knOCQUXTQpZvxrDTXnKkC0hhI4QBvKqNdv97kv1dO2XlPVaY0QN3wG+nuczEdz0xq4iCU29/2aeGQBpNMeT9gkzi4pUuaCiT4k1O60RLewqjASlFNPu65eMCbxaEySNn3qrbGRP8oVCxA/89WtqyqAWA9v7/7Up3veBz/eT3FW/7A3d/Xv6ip+7TQm/nnvAoj18nUEPh+9xfgIRABNNt98Z0Ly0j8CkpgSoyd6T36g9l6mhO7JDx9Kzt5WyRtI8C4Q5IdYvI/irOUep/vhM+XetoN20NV1MAYxV4GmExkfUyKi/okBYqbc10Vrgm9dA9GJpbqGtQ9js2ZeBp0vxV5Flta8xi+xVH8uujV4bTyyAAohfIkIoMnmm+9MlORIjLVOpwg1kVIGtaxLfEqWLhiUvGlfE/coGdOW8toeukTNZwRG196hHXR1Z6m+XNR4O4QD5S7mJJa623ny49dIZUfm6Git/JqqrAqablwEIGV1TqK73q+JCKAQQgTQZPPNdyZdQhRKgPXUhTKoZWcFUBU5ovoiwVMGxFnRGLxXH91OQiBQ7qdDjvx2dVsQE2LF4+h8Sbj5/HUKwzqoDEFBmU5uHIkjX6+9NRSMizCUr26teE/ZEQHEnCljXvLpaL4+t9dEBFAIIQJosvnmOxOSVk2IYpRkKQMvU2L1RI9IQID4rR7YS97Vv9rvUW/V8FkJGRAS9bREdV62B2vFnNS3UttShiDjs8ScCzzw/iO6E6DRGnIddE06/HrLl/sHxe7rqrI9/OfgNREBFEKIAJpsvvnORAnNE6I4I4B0WwdhoIduJT5GJwtHBZCfAHm7o/hvVtVbN2f96gQEP/V2VedLt4wkeLROdb11AsTaePmIvTWUP8bx06/R9da1wh91OnXCj4tGzX807msnAiiEEAE02XzznckoIcJRAaSTHpKlbu9AJxLgrAACnbjUE5wzKMn76UQ31haaXzevkS+JET+Vqbfktvx2bK0h41DHuLWuu94qQ/RojVjv7sRsy/c9EAEUQogAmmy++c7kEgJISbHeQhnxHAGkEwlevVzUU6YOJXdP6hJWo2eEKoqjrpdOdsDLvU7ipxMPakNdJ/JYM49xaw116sTc3BfvJcY8fokvF4ZbaB27NWeMoyLuFokACiFEAE0233xncgkBBHrouEJf6jwRP0cAefImARMvIEgkYtQW/7QhsaudknYVUBI0+CAeXr2+otjlC9/0kw9eax9iVx8YiTWPRbcReVXsjK22e2uomLRW+MGvYqRMbRE+GsPRGlZBQ3v5YRx8gYTUKKbXQARQCCECaLL55jsTHnAladUHiEFJtp7sKJHqs4sTJW4lXZXzqvZ7J0bUuX/BOPiVTyBJIxz89IITEJKwkjTQh77uD/CpxA2jmBwEgcfAWJTJT9fHxQ1jdm0AcSTxovb47UQI9aN4GYMxtQZaJ11TXws/fdK1A3yrfx2fNcafCyf6c8235nfrRACFECKAJptvvq8NiahOYIDEgp9ghNsBkTMSZggfrt1rPtU5QwRQCJeBv0TNX6Huvr/t1okAmmy++b42ED4kyXpKACRVnRIcfc4mzIVrg0jt6iRuOe3p6u+NCKDbQ9+XxVc9dPVHIRkf/VoJ4Gsk/Os+Xjr+Kli/7itKrg1fXfJa1y0CaLL55vva0G0V4LaJbqH4LRhujXR9w3p0yw2h6rcwdXLHNXwU8RoBdHtcSgDhA7q6in95qn/5atf2ltF3lp0RfpeEcfl+sq7ulokAmmy++b5GEEH+zAjwvnt+JdwWnNK54BH38EzPWSKAbo8VAogvBqUtX/Da1d8S+uJVTqxqnebBSVatm4G+Ub+L7ZaJAJpsvvmGENYQAXR7rBBAjHWJMWegb3Znnbr6lUiAvbbTswigyeabbwhhDRFAt0cE0Da3LIC4lUhsq06gnksE0GTzzTeEsIZHFUDcolDC54FZ3VYB3usWBq9exy2O0W/5kJC9Lc+D8HmUqPGDPz17w7MjPEi7JYDUh7Yax+N1VF/LHY3VIZ8SHN0YEk51jsyJfsRLG58jz+l4W6G2dW7yrbE61IYY+Uw/9w2df4QKZdTV9s+ZA6hd5/NWiQCabL75hhDW8KgCSImfRKmERbLTe179gVqv6/53zy0P6lRPe0+0CBtvT3KkHXXyr/YqJ+l6H4SaYiBu6l1AVRFCGXhZBZ/4YXza8spnUAJXXfUPtBuNzTw0F38PVUB06+HtacMYjKd1Yu6KVb/1petKf/fvayf/mhcwVhUsKlcc/h5GImhrvW6VCKDJ5ptvCGENjy6AlNiUQEmCLlwQGkqMnkT9NES+qKtJz0WUJ1gJB8b2crUH2qgclFirmNKpB3F7ufx42QjFU8eErYSufrVOY9c10fxqrBKQdT1Y83qasxUPZdTRxsvVp57g4V+ips6dMjg6B0G81PM8UFd/i0QATTbffEMIa3h0AUQS84QISuo18QIJlDr/37/KOvEASrAuXCSyuttKSrDuj0RNGb68rVCCd398Bm83QnPu5vASAdTNTyKyK5MQ3eKsANLaVREq9ODySEAenYPYWstbJQJosvnmew9sfcXGPcOvjPMnAc7+6jh/MgC6uksyYxzm/1r/btCjCyBPlLWuS2BdctsSM9D14TN4O9GN70ma8opicFGwNUYFH3VM8RIB5GVC/iR2XKDUth1nBZD+QGEnaAX1sFcm6hycrbW8VSKAJptvvvcAf1eGH3peu/p7hb97xLx57epH0Ae6ukty7XEkfGeIuWsQAfRyAcRn8HZO7cMpxFafbnz52OPWToC8TFR/mm93LTq24ul8bc1NUA97ZeLImmyNd2tEAE0233zvgUcVQPzhQOZd//L13okYddDVnWHWOCP03WGjr9a4dSKAXi6AdDtkxgnQUZEAW2NUuhjFkWRf67bGrv5YNz7XW1AjtuLR2vk66ZbiNU6AuhjyDNBtWQTQBB5VAI3YWw/qoKs7w6xx7pUIoJcLICW8+nCy0DNAnhC3RFP3DJBEwugZoA7aQ1dX6eYllOy7hK65VyGwNXYnHtT+Gs8A7Qksicu6tpSBl4mtGLbqbpUIoMnmm+89sJeIH41ZwmTWOPdKBNDLBZCeMUHUVEGj326qyVcPTjO+P5gr/3UM0HM+ndDCRxUo8uNlI7p5Cc3BY+VVc4Ca7FXuZaITCPKFoPL1QBBRp8+wJThH11VrV31xvUbrShl4mejmICRufR63TgTQZPPNdxb6EkweWt2q1wOt3XdG8f6zn/74O/qOEjFlo+dDNF5X1407+o4x2uob6PfaCo390z/+mXfUES913YO9qtPnegvKvyi24mujMt57/LyObmc51x6H9dMaqb37Fawf9fUZKD0bxXvdJmTtbu174iKAXi6AwMUAJwn4VSLshBGJ3eu9vZJrHQMf3oZ6kCCoc6EMvGzEaF7gsSIWPFadblUhQBl4mdD8vE+3HvJd/UhwAu2ISb5G19XXrvNfhRGorpZDNwdgHpQTk5ffOhFAk80331koWfOt7bVOiYxEx2fEj5IlyYzkRyKTAKjJkM9dOWXgZWIkgCR8eNW4iqX6V3LFF3Wg996uon61nQuLKvS6ZE9/94NoUgyUe1wuPDWGz1XvYU8EXXMcrQ1rznt8qg/jeFutVy1XXPys8Ur/kYhayaMKIBIi1wXxUOuURDsxoNtT/mvwgsTsSZUkSGIlKda2QAwSL0BfTnEU20iM4FOnFurHKU09cSDRg5eN2JoXEJOSPhA3ZepXhcDW2PJDfy/X3OhHPa+MU30D81U75q811tp117X6B2LpTpKAdtDVjeYgcdYJqlsmAmiy+eY7C5ImP5zdiYySnhIibUla9YSkCiVBYqO8JjjKwMuEkqSXyU8VaS7I/GSGzzWWI2gezNHLNT7UUw2tkQuj0bxH5UJjcC1csEikHp3TpceRoGFd6q/261TH578ngGr7W+NRBVAI1wDhw7/57nmpWyYCaLL55jsTJaZ6K6ITFyOU2LxslIi7tqITQIqj+7s6GsNPLfhMcu/a76HTLO+rkwqoQpE62vsajeY9KhfUQRWYoLhqecelx9GJTXerSmLHheGeABrFdStEAIVwOTgxem23vyACaLL55jsT/jdek5hOQ2oSA8QByZBEJmgL3k7lNeF1bUUVQIyl9j6e6JKqTmVI5pQfEXBCyV4nFBofn/IrIYBfPo9OjDymrXJBHXR1mueRuVx6HIk81gafjtbLf072BBD1Xn5rRACFcBl0+6u7dXnrRABNNt98Z0KSrycf3a0dIOmpbUdtSxmvXt61FUqS+qxkukd9boXPOjkCErXmtgXihvYSg/rMq96zNtRJONaxR/MelQvqoKs7Ix4uPY7ab3HmBCgCaNt88w3hNcNzR5wAvabf/hIRQJPNN9/Z1JMP3e5x0SABwIkA71UOlIOXjRJx11ZUAaQTmFH7PYhTPutJzQjmDbzXumgdKGdteK/nX+qpzGjeo3JBHXR1Z8TDpcfRCdBR4RIB9DLzzTeEsIYIoMnmm+9sJG5I6kpgJH9vo4RfxQ9QDl42SsRqW09k+KxTGy/XiVNtfwb57Z57qUj0ME/6eSJ30cNrJ6pG8x6VC+qgq1spgLaue0cE0MvMN98QwhoigCabb74rkEhQoqoJT8KgJlbdCgIvHyVinSjoVhIgKFRe/WjcKshEvU3XobkdEVGajxK/x686rVEnMkbz1m9ZjeZBHXR1Z8TDpcfRnLk+3frR1k/BIoBeZr75hhDWEAE02XzzXYGe+wEEQ61XYgOSKwleSU2nNN5+JAR02gQkVfchoeLtSboSR9TjD4ihjkuMtEG8qJ18jgRBhfFoL/zUSCc/XZ1gTOp49XLvy3zA26jO+4gz4uEa47B2lLPe/Jzgj1etrbfXzwm+3MeZOawkAiiEEAE02XzzXQHJXIKiJm+BeJEYUZKjTKcl3lYnEfUhYaCPkidjkmARHp0foI6YfGz6068KFGLSPKAKgCMoDvzUOiVy/NY62Jo3pymaN69+ysZY3XigMY/cwoNrjINP1QM+WCf3DfSlnjov15oencMqIoBCCBFAk8033xDCGiKAQggRQJPNN98QwhoigEIIEUCTzTffEMIaIoBCCBFAk8033xDCGiKAQggRQJPNN98QwhoigEIIEUCTzTffEMIaIoBCCBFAk8033xDCGiKAQggRQJPNN98QwhoigEIIEUCTzTffEMIaIoBCCBFAk8033xDCGiKAQggRQJPNN99wHH331Nmvu7gUq8e/R1hPvnajq7s2EUAhhAigyeabbzhOBND9EQEUARTCSiKAJptvvrfI1pdoriQC6P6IAIoACmElEUCTzTffW4SkBF3dSiKA7o8IoAigEFYSATTZfPO9RSKAeiKALk8EUARQCCuJAJpsvvnOhuT91e9999si5+u/7mvefP5HPvlURyJSecV9fPbTH/+KtvjD7xd/9ke/ot1P//hnnuqp+7kvfO7Nhz/0gbf7MC6CwtsL/NCW23Dyz5gjAUJ7yvDp/unj7eBTP/C9T/X4ol5juM+z44+g3WitHco8dsYmhm/5wPuePntbfFLGPLwctuJ7zjXjveKvIoWYO3/exqFOvlhXjd35nkUEUAghAmiy+eY7k49+97e/nXBIQKD31JMkea+kpjaqByVlErbqRklSCZk+EhO0UXsgyXofkqLEAH28vco9HtqrnnGoY56dsAE+Uy4xRl/8Kvaz44/YW2vBmtNO/mnDe8bWuN5e8XcxjATQS66Z1oAytdHc8MF7/Gltqj9w4Uu92m71mUEEUAghAmiy+eY7E5INSaurc5SsujqSVj3FcBHCSY/KlUyBfl6npExyVBmQTNUevyrXyQ14gscnbauQ0klGna/8Qx3b64+OP4J2e2uNfwk1hJDKNSeN530UXxfDSADh67nXzNfA67s6XVOfi9ozT79GlGvuEUD9xhxCuD4RQJPNN9+ZKBHVxFWhHXR1I3RKQGJTmZIfibaOKYFSk5+SchU0IBFyRIAAbcHLJCBGSfdS49Nub611+oNwqHX0Uyxervi7GLTeR9dn65qNYtdpTncrT319Pmrf3bIb/QzMIgIohBABNNl8852Jbl2Q3EiS/j9/hzbQ1QmSHT6EknWXTEcJjjr61bLR2PLXJXiSNUnZY+p8qbzzAV0fsTV+5cha78UigeJlW3324qOeOvGca6aTKYSN+wKJHe+rOfgYztZY1yYCKIQQATTZfPOdDf8TV+JTIqv/01edlwlEhvevnEmm6qPPxFHLHPkj2Xo5nxEa6lupbSmrPuC544/YW2uJpJE/iQcv24p/FN81rtkWfgIkweRjONRFAPUbcwjh+kQATTbffFdBUlSCJUl5nRKZlwGnGAgN4PaNJ3P5ek4y3SsT8ucJnnlQxhx47+07X1sCAro+ohv/CKO13otFfbxsq08X36Wv2Z6gqXRjOFtjXZsIoBBCBNBk8813NToZ8Gde+AzeDvQMDCcXte45ybQbRyc5R5/B0YO3VfxA539LQMDZ8c9Q11r+nvMMUHcNuvgufc221rtD7Ymj1uUZoAigEFYTATTZfPNdjZKsnwxIBNTnVvTQbk1YJDL1eakA8udIPCYlfvAEr/ZVlChW8HL5qe3F2fHPUNea9ZVPFxSsp05aQOVAO8pYb78+Ejrg8V36mskf8fn6CPp7XGpf42UeGn801rWJAAohRABNNt98Z0FiIgHzP3ISJCghk/S9rUQACYrkxCvlJDxPWvjwtryeSabUgZeRJOXLx5cfXj3BawwgFurUTv3cP/XVh3N2/I4za63ngABRIeHj7709qFzxybfKPb5LXzPw/sSPT14Vh/sDxSW/aqf13BrrmkQAhRAigCabb76zILGTaJT0gMTkyVKQNHXrAnivOk4OvE4+dALhyU+3OLy/Qx3x1PJuDE4M5K/GTF1NspTJh7dVnN0tGXF2/MqZtQbKJQokKrgGEgi1PXUSIfK9FV83H9o855oJTnYUHxA3fYijtlW8Wg/mqhgpw0/tM4MIoBBCBNBk8803hBEjARQuQwRQCCECaLL55hvCiAig6xIBFEKIAJpsvvmGMCIC6LpEAIUQIoAmm2++IYyIALouEUAhhAigyeabbwhhDRFAIYQIoMnmm28IYQ0RQCGECKDJ5ptvCGENEUAhhAigyeabbwhhDRFAIYQIoMnmm28IYQ0RQCGECKDJ5ptvCGENEUAhhAigyeabbwhhDRFAIYQIoMnmm2+4Tfj+Kr4ji9eu/lLwvWH1y0PDHO5dAP3if9uXhxC+TATQZPPNN9wm+vLQvS8FfSn6gtCtL2cN1+HeBdC3/f4/8pYI+iNtXQjhS0QATTbffMPlOPpN5kfgG9nxxWtXf5S9mPQt9t23qIfrcu8C6F/6PX/0zbd+zx+NCAphgwigyeabb7gc3EpCTPAVEl39Cm4xpvAlHkEARQSFsE0E0GTzzTdcjgigcIZHEUARQSGMiQCabL75zuT7P/adT8m4e95EiZo2Xs4zKiRvHgamTs+sfPV737353MpnP/3xt2/vAO8pq+249YN/tcNvjQF0K4k63tOOzxIW6l9x4aE5HImLuVFf53hmPTRGxWPSs0bMyfvC0bVRW5/XVtvwJR5FAP22b/+PI4JCGBABNNl8850JCZHk2CXGkQCijMSq5OrvoRNBH/7QB96uJ4F7EmcctdNzNiRr3jO2fLtIAMWHYEB0SIhQRj19Na6SP0jcIFgolw/qGFMChs8+Hp+7csqOrgd9t2ICzcHXBc6sDeKnjsEcR88dhS/xKALo/f/mFyKCQhgQATTZfPOdCYmRRMlrrdsSQIBQ8CStExKSrrdXMqa9n2ogQEjkKtN4JHXq1A50KuJCQe1HfbxNFQjAr5vTr5608Jk+dR6jtVIMR9djKyboBJD6HF0b+WCO3jZs8ygC6Jv/wJ//ChH0bb8/IigEEQE02XzznclLBFB3i0anJ17WJegOnYwgmGqdYvETDJUxZid+vM1IbIygD3jZngA6uh57MXUC6OzadD7CPo8kgN4pgvIr8iFABNBk8813Ji8RQF4mupMHTkAo6wSCo9s5JHvGdCQAXDTsCQk40gbxhLDw8egD3k7lvHp511Z067EXk/q4eDm7Njp1U5+9tQ9f4tEEUERQCO8kAmiy+eY7ExIoSZLXWvcSAeTJe6u9o3ZbdCdAIyEBe22Ym05pOmpbyl66Hnsxba3hFr42QH8JJ8DvowohJf7wJVwARQSF8JVEAE02Nul7FUASGHvPoyhZe98tFN9ISMBWG52SMG69tUQ5eNlKAXR2bRxEj06JuBaP+FxQFQCPThVAEUEhfJkIoMnGJr1SAPEwcq3TQ7wvTfgq23sGSM8Kdc+5dOwJCdhqszUe5eBlKwXQ2bXpkAjqfkvv3tlK/OHLvD8iKIQIoNnGJr1CAOkUpJ4MSPzASxM+wocyngXyMXj2Bt+6LaN2nHaMfqPrzLM0gG+NXeskCOr8FAd4+aUE0FZMsLWGR9emQ/N9iYh6rUQAHef9EUHhwYkAmmxs0isEEOj2CiKI5KuHllX+0oQP8gW0Ubva1m/V6G/d8KqYvC3v5U9lHerLq8amXP2BcRlLcem2nfu5lACCUUxbfc6sja4lbXxeIwF170QAneP9EUHhgYkAmmxs0qsEEAlRyVVJklMCnVTUhE9yBS8TSrQ61RE67VGyBtp2t8Uokx9gLG4B1ZMLxVcf/q3QzkWezwefVZxRpltO7kenYvUW0nPWYysmjV37wNG1oczXmvdc40cUPxABdJ73RwSFByUCaLKxSa8SQCHcOxFAzyMiKDwiEUCTjU06AiiE6xAB9HwigsKjEQE02dikI4BCuA4RQC/jSQR9239kIqhPHCHcAxFAk41NOgIohOsQAfQyIoDCIxEBNNnYpCOAQrgOEUDPp94C+4Wf6ZNGCPdCBNBkY5OOAArhOkQAPY+vFD9/NOInPAQRQJONTToCKITrEAF0noif8KhEAE02NukIoBCuQwTQOSJ+wiMTATTZ2KQjgEK4DhFAx3Hx863fE/ETHo8IoMnGJh0BFMJ1iAA6RsRPCBFA041NOgIohOsgARS+RMRPCGMigCYbm3QEUAjXoQqARyfiJ4QxEUCTjU06AiiEtejf4CrzzfcadAIo4ieEryQCaLL55htCWMOjCaCInxDeSQTQZPPNN4SwhkcSQBE/IfREAE0233xDCGt4FAEU8RPCmAigyeabbwhhDY8igCJ+QhgTATTZfPMNIazhUQRQxE8IYyKAJptvviGENTyKAIr4CWFMBNBk8803hLCGRxBAET8hbBMBNNm0+YYQ1rJaAP3qX/6Bq/Ft3/MfvvnFv/AftHUhBCMCaJ7pf6AhhPWsMv0P9Fr84n/bl4cQ3sk92k0KoFgsFovFYrFrWgRQLBaLxWKxh7MIoFgsFovFYg9nEUCxWCwWi8UezN68+f8AuM5PqkAoN4cAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of these steps were previously implemented in ZeMA-machine-learning-tutorials found on this link: https://github.com/harislulic/ZeMA-machine-learning-tutorials/tree/computation_without_uncertaintites .\n",
    "\n",
    "This and the following notebooks are supplements of the previous tutorials, with focus on the uncertainty propagation in steps 1) and 2). Some parts from the previous tutorials are also presented here in order to show the comparison of results with and without uncertainties involved.\n",
    "\n",
    "This notebook will explain the basics of (FFT) and principles of using PyDynamic software for the uncertainty propagation.\n",
    "File [2_Features_extraction_and_selection.ipynb](2_Features_extraction_and_selection.ipynb) contains steps 1) and 2) supplemented by uncertainties. That file is very intensive and because of that, steps 1) and 2) can be found separately in  [2a_Features_extraction.ipynb](2a_Features_extraction.ipynb) and [2b_Features_selection.ipynb](2b_Features_selection.ipynb).\n",
    "\n",
    "\n",
    "The workflow for these notebooks is:\n",
    "![workflow.PNG](attachment:workflow.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.3 Fast Fourier transform (FFT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A first idea for machine learning on the dataset from the previous notebook is to try to use fast Fourier transform to make useful features from time series. That's why we will first examine FFT on these time series."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fourier showed that any periodic signal s(t) can be written as a sum of sine waves with various amplitudies, frequencies and phases [2]:\n",
    "\n",
    "<img src=\"pictures/fft.png\">\n",
    "\n",
    "Results of this are amplitudes and phase shifts for frequency range from 0 to half the sampling frequency and it is covered by (N/2)+1 points where N is the size of the window i.e. number of samples. From Nyquist's theorem we know that the largest frequency component in the original signal must be half the sampling frequency. So, from a 2000 point FFT of signals sampled at 2000Hz we get 1001 unique spectral points covering the range 0 to 1000Hz.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py                                                 #importing the h5 package\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from matplotlib import colors                                #ploting with colors\n",
    "import scipy.fftpack                                         #fourier tranformation\n",
    "from collections import OrderedDict\n",
    "import matplotlib as mpl\n",
    "font = {'family' : 'DejaVu Sans',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 20}\n",
    "mpl.rcParams['figure.figsize'] = (20,10)\n",
    "mpl.rc('font', **font)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.1 Data import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    \n",
      "Input matrices have dimensions: (2000, 6291), where 2000 represents number of measurements in time\n",
      "and 6291 represents number of cycles.\n"
     ]
    }
   ],
   "source": [
    "url = 'https://zenodo.org/record/1326278/files/Sensor_data_2kHz.h5'\n",
    "\n",
    "def get_filename(url):\n",
    "    return url.split('/')[-1]\n",
    "\n",
    "filename = get_filename(url)                    # Data filename.\n",
    "f = h5py.File(filename, 'r')                    # Importing the h5 file. \n",
    "    \n",
    "#print(\"Keys: %s\" % f.keys())\n",
    "a_group_key = list(f.keys())[0]\n",
    "\n",
    "data = list(f[a_group_key])                       # Transforming data into list\n",
    "\n",
    "sensorADC=[]                                      # Initialising a list \"sensor\" and\n",
    "for i in range(len(data)):                        # Filling it with data from all sensors \n",
    "    sensorADC.append(pd.DataFrame(data[i]))\n",
    "\n",
    "for i in range(len(data)):                             \n",
    "    sensorADC[i]=sensorADC[i].iloc[:,:-1]         # Cuting the last cycle because it contains all zero elements.\n",
    "\n",
    "print(\"\"\"    \n",
    "Input matrices have dimensions: %s, where %s represents number of measurements in time\n",
    "and %s represents number of cycles.\"\"\" % (np.shape(sensorADC[0]),np.shape(sensorADC[0])[0],np.shape(sensorADC[0])[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.2 Converting into SI units "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset=[0, 0, 0, 0, 0.00488591, 0.00488591, 0.00488591,  0.00488591, 1.36e-2, 1.5e-2, 1.09e-2]\n",
    "gain=[5.36e-9, 5.36e-9, 5.36e-9, 5.36e-9, 3.29e-4, 3.29e-4, 3.29e-4, 3.29e-4, 8.76e-5, 8.68e-5, 8.65e-5]\n",
    "b=[1, 1, 1, 1, 1, 1, 1, 1, 5.299641744, 5.299641744, 5.299641744]\n",
    "k=[250, 9.81, 98.1, 98.1, 1.25, 100000, 30, 0.5, 2, 2, 2]\n",
    "units=['[Pa]', '[mm/s^2]', '[mm/s^2]', '[mm/s^2]', '[kN]', '[Pa]', '[mm/s]', '[A]', '[A]', '[A]', '[A]']\n",
    "\n",
    "pd.set_option('mode.chained_assignment', None)\n",
    "\n",
    "sensor = list(sensorADC)\n",
    "for i, df in enumerate(sensor):\n",
    "    for row_label, row in df.iterrows():\n",
    "        sensor[i].iloc[row_label,:]=((sensorADC[i].iloc[row_label,:]*gain[i])+offset[i])*b[i]*k[i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### If you have problems with previous step, you can skip conversion into SI units by running next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor=sensorADC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*One cycle of microphone sensor will be analysed with `FFT`.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 0 represents sensor number (microphone)\n",
    "# 3000 represents cycle number from microphone\n",
    "df=sensor[0][3000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time axis is formed based on informations about data (`0_Data_Import_and_Visualisation.ipynb`). Sampling period is 1 second, and number of points is 2000 i.e. length of this time series.\n",
    "\n",
    "After fast Fourier transform, this signal can be represented in the frequency domain. (**Amplitude-Frequency** plot) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = 1                                                                   # sampling period 1 s\n",
    "n_of_sampling_pts=2000                                                     # sampling points\n",
    "\n",
    "x = np.true_divide(df.index, n_of_sampling_pts)                            # time values\n",
    "freq = np.fft.rfftfreq(n_of_sampling_pts, float(time)/n_of_sampling_pts)   # frequency axis\n",
    "amp = np.fft.rfft(df)   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rcParams['figure.figsize'] = (20,10)\n",
    "mpl.rc('font', **font)\n",
    "\n",
    "plt.figure(figsize=(20,12))\n",
    "\n",
    "plt.subplot(211)\n",
    "plt.plot(x.values, df, label=\"Real values\")\n",
    "plt.ylabel(\"Microphone (Pa)\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.plot(freq,np.abs(amp))\n",
    "plt.xlabel(\"Frequency (Hz)\")  \n",
    "plt.ylabel(\"Amplitude (Pa)\")  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of idea to use 10% of these amplitudes for machine learning workflow, which should be enough to reconstruct the signal and make prediction, reconstruction using 10% highest amplitudes will be shown here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Amplitudes are sorted by their absolute values from highest to lowest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coeffs = dict(zip(freq, amp))\n",
    "\n",
    "# take top N % of coefficients  !!!! This can be changed !!!!\n",
    "N=10\n",
    "\n",
    "# taking top N % aplitudes and sorting them by value, t[1] means sorting by second column i.e. by values of dictionary\n",
    "out = OrderedDict(sorted(coeffs.items(), key=lambda t: abs(t[1]), reverse=True)[:round((N/100.0)*len(freq))])  \n",
    "\n",
    "sorted_amp=list(out.values())       #transforming sorted dictionary to the list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keeping only 10% choosen frequencies. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "treeshold_apsolute_amplitude=abs(sorted_amp[-1])      #choosing the treeshold amplitude\n",
    "\n",
    "for key, value in coeffs.items():\n",
    "        if abs(value) <= treeshold_apsolute_amplitude:\n",
    "             coeffs[key]=0+0j\n",
    "                \n",
    "top_10=list(coeffs.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another plot in frequency domain, will show which frequencies out of full spectrum are selected and why."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(freq,np.abs(amp),zorder=1, linewidth=4,label='Full spectrum') \n",
    "plt.ylabel(\"Amplitude (Pa)\") \n",
    "\n",
    "plt.plot(np.abs(top_10),zorder=2, label='10% frequencies')\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next plot shows sorted amplitudes, but without frequency labels.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.stem(np.abs(sorted_amp), markerfmt=' ', use_line_collection=True)\n",
    "plt.ylabel(\"Amplitude (Pa)\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reconstructing using inverse fast Fourier transform will give similar sensor signal, and that is why these frequencies and amplitudes are very good features for comparing different cycles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inverse_fft = np.fft.irfft(top_10)\n",
    "\n",
    "plt.plot(x.values, df, label=\"Real values\")\n",
    "plt.ylabel(\"Microphone (Pa)\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.plot(x.values,inverse_fft, label=\"Reconstructed values\")\n",
    "plt.xlim(0,1)\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This procedure can be done for every sensor and every cycle of each sensor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 PyDynamic for the uncertainty propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to perform uncertainty evaluation, the PyDynamic software is used. The software is based on Discrete Fourier Transform (DFT) with some FFT algorithm.\n",
    "The time domain signal is *x(t).* A white noise will be added to this signal:\n",
    "\n",
    "$$x_{n}(t) = x(t)+\\epsilon$$\n",
    "\n",
    "White noise has normal distribution ${\\mathcal {N}}(0 ,\\sigma ^{2})$ and standard deviation that can be specified by user. \n",
    "\n",
    "The time domain signal that will be analyzed is $x_{n}$ with associated uncertainty  $u_{x}$ as squared standard deviation representing white-noise variance ($\\sigma^2$).\n",
    "White noise is a signal whose spectrum has equal power within any equal interval of frequencies.\n",
    "\n",
    "The result of the DFT is the vector X of real and imaginary parts of the DFT applied to  $x_{n}$ and the associated uncertainty UX ( covariance matrix associated with  or as a variance in the case of white noise.) When only white noise is considered, all off-diagonal elements of UX are equal to zero. [5] \n",
    "\n",
    "Sampling period is 1 second, and number of points is 2000 i.e. length of this time series.\n",
    "\n",
    "*One cycle of microphone sensor will be analysed with `DFT`.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PyDynamic import __version__ as version\n",
    "version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma=float(input(\"Assume standard deviation:\"))                           #signal noise standard deviation\n",
    "sensor_num = int(input(\"Enter the number of sensor (0-11):\"))\n",
    "cycle=int(input(\"Enter the cycle(0-6291):\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if you did not convert the signal to SI units, convert here\n",
    "offset=[0, 0, 0, 0, 0.00488591, 0.00488591, 0.00488591,  0.00488591, 1.36e-2, 1.5e-2, 1.09e-2]\n",
    "gain=[5.36e-9, 5.36e-9, 5.36e-9, 5.36e-9, 3.29e-4, 3.29e-4, 3.29e-4, 3.29e-4, 8.76e-5, 8.68e-5, 8.65e-5]\n",
    "b=[1, 1, 1, 1, 1, 1, 1, 1, 5.299641744, 5.299641744, 5.299641744]\n",
    "k=[250, 1, 10, 10, 1.25, 1, 30, 0.5, 2, 2, 2]\n",
    "units=['[Pa]', '[g]', '[g]', '[g]', '[kN]', '[bar]', '[mm/s]', '[A]', '[A]', '[A]', '[A]']\n",
    "\n",
    "sensor[sensor_num][:][cycle]=((sensor[sensor_num][:][cycle]*gain[sensor_num])+offset[sensor_num])*b[sensor_num]*k[sensor_num]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adding white noise to the signal\n",
    "sensor[sensor_num][:][cycle]=sensor[sensor_num][:][cycle]+np.random.randn(sensor[sensor_num].shape[0])*sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import PyDynamic\n",
    "ux=sigma**2                                                                \n",
    "time = 1                                                                   # sampling period 1 s\n",
    "n_of_sampling_pts=2000                                                     # sampling points\n",
    "time_steps=np.arange(0, time, float(time)/n_of_sampling_pts) \n",
    "\n",
    "# timesteps for the period (0 - 1) s\n",
    "X,UX=PyDynamic.uncertainty.propagate_DFT.GUM_DFT(sensor[sensor_num][:][cycle],ux)\n",
    "print(\"Vector X has dimensions:\",X.shape)                                         # vector of real and imaginary parts\n",
    "print(\"Vector UX has dimensions:\",UX.shape)                                        # vector/covariance matrix of uncertainites\n",
    "freq=PyDynamic.uncertainty.propagate_DFT.GUM_DFTfreq(2000, float(time)/n_of_sampling_pts)\n",
    "# frequency axis (array of length N//2 + 1 containing the sample frequencies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First plot represents time domain signal values. The second and third plot represent real and imaginary parts of the vector X, with their associated uncertainty. \n",
    "*Note: au - arbitrary unit*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mpl.rcParams['figure.figsize'] = (20,10)\n",
    "mpl.rc('font', **font)\n",
    "\n",
    "plt.figure(figsize=(20,12))\n",
    "plt.plot(time_steps, sensor[sensor_num][:][cycle])\n",
    "plt.xlabel(\"time (s)\")\n",
    "plt.ylabel(\"Microphone (Pa)\")\n",
    "plt.figure(figsize=(20,12))\n",
    "plt.subplot(211)\n",
    "plt.errorbar(freq, X[:len(freq)], xerr=np.sqrt(UX[:len(freq)]))\n",
    "plt.ylabel(\"real parts (au) \")\n",
    "plt.xticks([])\n",
    "plt.subplot(212)\n",
    "plt.errorbar(freq, X[len(freq):], xerr=np.sqrt(UX[len(freq):]))\n",
    "plt.ylabel(\"imaginary parts (au)\")\n",
    "plt.xlabel(\"frequency / Hz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.1 Transformation from real and imaginary parts to amplitudes and phases\n",
    "\n",
    "DFT gives complex values obtained from real and imaginary parts as results. Nevertheless, it is possible to calculate amplitudes and phases instead. \n",
    "\n",
    "Amplitudes and phases are calculated on the basis of vector X  and covariance matrix of uncertainties, UAP on the basis of  uncertainties UX. The shape of matrix of uncertainties, UAP is (2M,2M), where *M* corresponds to the number of (*N*/2)+1 points (*N* - number of sampling points). The shape of matrices of amplitudes and phases, respectively A, P is (M,M).\n",
    "\n",
    "*Note: When amplitudes are small relative to the uncertainty associated with real and imaginary parts , the GUM uncertainty propagation becomes unreliable and a Monte Carlo method is recommended instead. Consequently, GUM2DFT does raise a warning to the user and recommends using a Monte Carlo method instead whenever an element of  is below a pre-defined threshold. The default threshold in GUM2DFT is 1.0, but may be adjusted for specific applications.[7]*\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A,P,UAP=PyDynamic.uncertainty.propagate_DFT.DFT2AmpPhase(X, UX, keep_sparse=False, tol=1.0, return_type='separate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned before, because of the white noise, UAP matrix contains lots of zeros (all of non-diagonal elements). In this case, UAP matrix is dense matrix. Operations using standard dense-matrix structures and algorithms are slow and inefficient when applied to large sparse matrices as processing and memory are wasted on the zeroes. If matrix has large dimensions, sparse data is by nature more easily compressed and thus requires significantly less storage [6]. Amplitudes and phases with their uncertainties are shown on the following plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,12))\n",
    "plt.subplot(211)\n",
    "\n",
    "plt.errorbar(freq, A,xerr=np.sqrt((np.diag(UAP)[:1001])), ecolor='black') #other diagonals are zero\n",
    "plt.ylabel(\"Amplitude\") \n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.subplot(212)\n",
    "plt.errorbar(freq, P,xerr=(np.sqrt((np.diag(UAP))[1001:])),ecolor='black')\n",
    "plt.xlabel(\"Frequency (Hz)\")\n",
    "plt.ylabel(\"Phase\")   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4.2 Reconstruction of time domain signals\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reconstruction using 10% highest amplitudes can be made with iDFT function or AmpPhase2Time from PyDynamic.\n",
    "The first one was implemented into function `reconstructed_signal_with_iDFT` with arguments (*parts*,*unc*,*perc*), where *parts* denotes vector of real and imaginary parts, *unc* their uncertainties and *perc* the percentage of the spectrum that will be taken into consideration.\n",
    "\n",
    "From real and imaginary parts contained in vector *parts*, amplitude in every frequency can be calculated as\n",
    "$A=\\sqrt{(Real part)²+(Imaginary part)²}$. \n",
    "\n",
    "Indices of 10% of highest amplitudes are found by function `np.argsort`.\n",
    "\n",
    "Values of real and imaginary parts of 10% of spectrum are assigned to vector *F*, on the indices that were found by `np.argsort`. Other values are considered as zeros.\n",
    "\n",
    "Values of uncertainties of real and imaginary parts of 10% of spectrum are assigned to vector *UF*, on the indices that were found by `np.argsort`. Other values are considered as zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstructed_signal_with_iDFT(parts,unc,perc):\n",
    "    Amplitudes=[]\n",
    "    range_F=unc.shape[0]-2\n",
    "    for i in range(range_F//2):\n",
    "        Amplitudes.append(np.sqrt(parts[i]**2+parts[range_F//2+i]**2)) #calculating amplitudes\n",
    "    Sorted_amplitudes=np.argsort(Amplitudes)[::-1] # indices of amplitudes in descending order\n",
    "    num=round((perc/100.0)*len(Amplitudes))\n",
    "    N_highest_amplitudes=Sorted_amplitudes[:num]\n",
    "    F=np.zeros(len(parts)) #Vector (2*len(parts),)\n",
    "    UF=np.zeros(len(parts)) #Vector (2*len(parts),)\n",
    "    for i in N_highest_amplitudes:\n",
    "        F[i]=parts[i] # assigned real parts on indices of amplitudes \n",
    "        F[i+range_F//2+1]=parts[i+range_F//2+1] # assigned imaginary parts on indices of amplitudes \n",
    "        UF[i]=unc[i] # assigned uncertainties of real parts on indices of amplitudes \n",
    "        UF[i+range_F//2+1]=unc[i+range_F//2+1] # assigned uncertainties of imaginary on indices of amplitudes \n",
    "    x,ux=PyDynamic.uncertainty.propagate_DFT.GUM_iDFT(F, UF, Nx=None, Cc=None, Cs=None, returnC=False)\n",
    "    return x,ux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function is executed on the previously obtained *X*,*UX* and *N*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=int(input(\"Enter the percentage of choosen amplitudes (0-100):\"))\n",
    "x1,ux1=reconstructed_signal_with_iDFT(X,UX,N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_steps, df, label=\"Real values\")\n",
    "plt.ylabel(\"Microphone (Pa)\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.errorbar(time_steps,x1,xerr=np.sqrt(np.diag(ux1)),label=\"Reconstructed values with DFT\")\n",
    "\n",
    "plt.xlim(0,1)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reconstruction using 10% highest amplitudes can be also made with AmpPhase2Time from PyDynamic.\n",
    "This was implemented into function `reconstructed_signal_with_AmpPhase2Time` with arguments (*amplitudes*, *phases*, *uncertainties*, *perc*), where *amplitudes* denotes vector of amplitudes, *phases* denotes vector of phases and *uncertainties* is the uncertainty matrix. *perc* denotes the percentage of the highest amplitudes that will be taken into consideration.\n",
    "\n",
    "Values of amplitudes for 10% of spectrum are assigned to the vector *A1*, on the indices that were found by `np.argsort`. Other values are considered as zeros.\n",
    "\n",
    "Values of phases that correspond to the amplitudes are assigned to the vector *P1*, on the indices that were found by `np.argsort`. Other values are considered as zeros.\n",
    "\n",
    "Values of uncertainties of 10% of spectrum are assigned to vector *UAP1*, on the indices that were found by `np.argsort`. Other values are considered as zeros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reconstructed_signal_with_AmpPhase2Time(amplitudes,phases,uncertainties,perc):\n",
    "    Amplitudes=amplitudes\n",
    "    Phases=phases\n",
    "    Uncert=uncertainties\n",
    "    Index_amplitudes=np.argsort(Amplitudes)[::-1]\n",
    "    num=round((perc/100.0)*len(Amplitudes))# indices of amplitudes in descending order\n",
    "    N_highest_amplitudes=Index_amplitudes[:num]\n",
    "    A1=np.zeros(len(Amplitudes)) \n",
    "    P1=np.zeros(len(Amplitudes))\n",
    "    UAP1=np.zeros(Uncert.shape)\n",
    "    range_F=N = Uncert.shape[0]-2\n",
    "    for i in N_highest_amplitudes:\n",
    "        A1[i]=Amplitudes[i] # assigned values of amplitudes on indices of N% amplitudes \n",
    "        P1[i]=Phases[i]  # assigned values of phases on indices of N% amplitudes \n",
    "        UAP1[i,i]=Uncert[i,i]                             # Covariance matrix contains uncertainties of amplitudes, phases and \n",
    "        UAP1[i,i+range_F//2+1]=Uncert[i,i+range_F//2+1]   # covariance uncertainties of amplitudes and phases.\n",
    "        UAP1[i+range_F//2+1,i+range_F//2+1]=Uncert[i+range_F//2+1,i+range_F//2+1]         \n",
    "    x,ux=PyDynamic.uncertainty.propagate_DFT.AmpPhase2Time(A1,P1,UAP1)\n",
    "    return x,ux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D=int(input(\"Enter the percentage of choosen amplitudes (0-100):\"))\n",
    "x2,ux2=reconstructed_signal_with_AmpPhase2Time(A,P,UAP,D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time_steps, df, label=\"Real values\")\n",
    "plt.ylabel(\"Microphone (Pa)\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.errorbar(time_steps,x2,xerr=np.sqrt(np.diag(ux2)),label=\"Reconstructed values with DFT\")\n",
    "\n",
    "plt.xlim(0,1)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References:\n",
    "\n",
    "[1]  PTB, ZeMA, - Deep dive into the ZeMA machine learning (ppt), January 2019\n",
    "\n",
    "[2]  https://www.nti-audio.com/en/support/know-how/fast-fourier-transform-fft\n",
    "\n",
    "[3]  http://www.sthda.com/english/wiki/correlation-test-between-two-variables-in-r\n",
    "\n",
    "[4]  https://en.wikipedia.org/wiki/Pearson_correlation_coefficient\n",
    "\n",
    "[4]  Edouard Duchesnay, Tommy Löfstedt, - Statistics and Machine Learning in Python, March 2018\n",
    "\n",
    "[5]  https://pydynamic.readthedocs.io/en/master/Getting%20started.html#quick-examples\n",
    "\n",
    "[6]  https://en.wikipedia.org/wiki/Sparse_matrix\n",
    "\n",
    "[7]  S Eichstädt and V Wilkens, - GUM2DFT—a software tool for uncertainty evaluation of transient signals in the frequency domain, March 2016\n",
    "\n",
    "[8]  Allan V. Oppenheim, Ronald W. Schafer, John R.Buck - Discrete Time Signal Processing , 1999\n",
    "\n",
    "[9]  https://dsp.stackexchange.com/questions/69/why-is-the-fourier-transform-so-important\n",
    "\n"
   ]
  }
 ],
 "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": 2
}