notebooks/anomaly_detection/FFT_Outlier_Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Transform Based Outlier Detection in Time Series data using CuPy and CLX libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Authors:\n",
    "   Gorkem Batmaz (Nvidia)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "* Introduction\n",
    "* Implementation using CuPy and CLX\n",
    "* References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outlier detection is a  technique for detecting outliers in the observed data or in time series. In real life, signals of interest can have many sinusoidal components and a lot of noise. Therefore it is not always straightforward to create rules to catch anomalies. We will demonstrate how we can detect outliers in time series data using Fourier transform.\n",
    "We quote from -Brighan, 1974, The FFT- \"..., the essence of the Fourier Transform of a waveform is to decompose or separate the waveform into a sum of sinusoids of different frequencies. If these sinusoids sum to the original waveform then we have determined the Fourier Transform of the waveform.\" Mathematically, this transform is stated as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$X(\\omega)=\\int_{-\\infty}^{\\infty} x(t) e^{-t \\omega t} d t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use discrete Fourier transform for performing Fourier analysis of discrete-time sequences. It will help us denoise the signal and allow us to employ statistical methods to detect anomalies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Discrete Fourier Transform $ X[k] $ of $ x[n] $ is the discrete-time signal is  \n",
    "\n",
    "$ X[k]=\\sum_{k=0}^{N−1} x[n].e^{−J.2pi.kn/N} $\n",
    "\n",
    "where $ x[n] $ is a discrete-time signal with Period N.\n",
    "\n",
    "The Inverse Discrete Fourier transform is\n",
    "\n",
    "$ x[n]=(1/N)\\sum_{k=0}^{N−1} X[k].e^{J.2pi.kn/N} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do this, we will first move the signal into the frequency domain. Then we will apply kind of a filter to remove components that have too high frequencies. Once we have a more precise signal, it will be possible to use techniques such as rolling z score to detect anomalies.\n",
    "In this notebook, we will use an alternative dataset not to disclose internal information. We will use a dataset of daily minimum temperature measurements taken in Melbourne, Australia between 1981 and 1990. This dataset is available to download at https://www.kaggle.com/paulbrabban/daily-minimum-temperatures-in-melbourne"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Importing the libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cupy as cp\n",
    "import cudf\n",
    "from clx.analytics.stats import rzscore\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pyplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read the data into GPU memory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf=cudf.read_csv(\"daily-min-temperatures.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0       20.7\n",
       "1       17.9\n",
       "2       18.8\n",
       "3       14.6\n",
       "4       15.8\n",
       "        ... \n",
       "3645    14.0\n",
       "3646    13.6\n",
       "3647    13.5\n",
       "3648    15.7\n",
       "3649    13.0\n",
       "Name: Temp, Length: 3650, dtype: float64"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf['Temp']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Convert cudf series into CuPy to utilise the Fast Fourier Transform functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/rapids/lib/python3.7/site-packages/cudf/io/dlpack.py:74: UserWarning: WARNING: cuDF to_dlpack() produces column-major (Fortran order) output. If the output tensor needs to be row major, transpose the output of this function.\n",
      "  return libdlpack.to_dlpack(gdf_cols)\n"
     ]
    }
   ],
   "source": [
    "a = cp.fromDlpack(gdf['Temp'].to_dlpack())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([20.7, 17.9, 18.8, ..., 13.5, 15.7, 13. ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualise the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f23610f1b90>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7LUlEQVR4nO2deXgURfrHv5UDAiTcSUAghBs5FDCAiIAgIgIuoq7Xqniiu+jqqvtb1lvXI6ur7nqsioLiBbuuByiIHIrc9w0B5AhnINyEI9dM/f6Y6aQzmZ7prq6u6p6pz/PwkEy6u97prnr7rbfeel9CKYVCoVAovEeCbAEUCoVCwYZS4AqFQuFRlAJXKBQKj6IUuEKhUHgUpcAVCoXCoySJbKxx48Y0OztbZJMKhULheVatWnWEUpoe+rlQBZ6dnY2VK1eKbFKhUCg8DyFkd7jPlQtFoVAoPIpS4AqFQuFRlAJXKBQKj6IUuEKhUHgUpcAVCoXCoygFrlAoFB5FKXCFQqHwKEqBW6CouAxT1+6XLYZCoVAAELyRx+uM+3oDpq8vQJv0VHRpVk+2OAqFIs5RFrgFDp4sBgAUl/kkS6JQKBRKgVtCVS9SKBRuQilwBgiRLYFCoVAoBa5QKBSeRSlwhUKh8ChKgSsUCoVHiarACSEtCCE/E0LyCCGbCCEPBT9/lhCynxCyNvhvmPPiyqG4zIfDRSVQS5gKhcJNmIkDLwfwKKV0NSEkDcAqQsjs4N/eoJT+wznx3MEdHy3H0p3H0D2rfvATtYqpUCjkE9UCp5QWUEpXB38uApAHoJnTgoWy68gZPDttEwpPFeOvX69HablfWNtLdx4DAKgoQoVC4SYs7cQkhGQD6A5gGYC+AB4ghNwOYCUCVvrxMOeMATAGALKyspgFvf/TVdh6qAiLdxzBtkOncUmbxrj6wvOYr6dQKBRex/QiJiEkFcBXAB6mlJ4C8C6ANgC6ASgA8Fq48yil4ymlOZTSnPT0ajU5TbP1UBEAYNuh04HrMl+JHRX/rVAo3IQpBU4ISUZAeX9OKf0aACilhyilPkqpH8AHAHo5J6Y7UC4UhULhJsxEoRAAEwDkUUpf133eVHfYKAAb+YvnTpQlrlAo3IAZH3hfALcB2EAIWRv87HEANxNCuiHgzcgHcJ8D8rkKZYDzZ+vBIpT5/Cq7o0LBQFQFTildiPBxczP4i+MNlAHOjyv/OR8AkJ87XLIkCoX3UDsxFYo4oKTch+xx0/HVqn2yRVFwxLMKXKV2VSjMc/R0KQDgH7O2SpPB56f4bt0BNXY54lkFrlAozOMGlTlx4S48OHkNvl6tyhLyQilwKyjLQeFxZK7fHDwVqGh17EypRCliC6XAGSAqjlDhMdzgttBEUMOHH0qBKxRxRKjxQSnFqz9uwY7Dpx1vm7rCkRNbKAVuAdX9FLFGYVEJ3vl5B26fsNzxtiotcGWC80IpcIUijtGUarlfXHZPpb75YSkbYbzjAjeiQmGZbYeKMH7+zrB/U24Nb6MUuAW0zq4sCIWXuOH9JThxtiziMURAr9YWUpUHhR/KhWKBjftPAVC+cIW30Bc/cYPydIEIMYNS4AyoDqiIFUS6BZXhwx+lwBkQbcVMWpyPi1+aK7ZRF3GquAxlPnGLbPGCz09x8lzAtSKyT6soFH4oBc6ACH+hnmembarYxRaPXPDsLIz9fLVsMWKOZ6ZtxFX/WiCsPRUEwB+lwKNQXOar9hkhwPj5O7Bx/0kJEsUnszYfki1CTKA3fqetPVD5uYC2tSCABGWAc8OzCvzZaZuEtDNl+Z6wn780YwtGvLVQiAwaHy7YiW7PzxLapsL7GFm+iYI1aYUcYVwoJ8+Vwed33kR/aMoaPDh5jePtiMKzCvx4lLAoXpQL6FRmeWF6XtRwMIXCLAmCfdFGI6nM58eFz83Ck99ucFyGqWsP4Lt1B6If6BE8q8AV8UdJuQ9/+34ziorVS4wV/fqNXn+LXFgMbancF1DtKs2sdZQCV3iG/67chwkLd+GN2b/KFiUmEB0NYuTKcVKMs6XlKCmvvo4VKygFHoVwnW79PrmLl/EaUucLfm+fwLwdsYB+u7xeWSYKD+cLvxNTG2NOOCs7Pf0jrnh9vgNXdgdKgTPw+DfO++oi8dbc+LRARa1GrMw/hi+WhV+8jiVkRYOEhuFqL5jScr8jecv3HDvL/ZpuwdO5UIrLfEhJTpQthnBiuUOGEm5AOz31v/69JQCAW3pnOdqObBIT3eFC0X9+4GQxmtWvJUagGMDTFniXZ350vA03Zmvbe/ycbBGEoTZ/2Ed/D/UqOzmhcviL8KYYVeRRj5gdTytwN4X4iWTV7uOyRRCGE0/4cFEJXpm5BX7B/Wfhr0fw7ZrIkRazNx/Cj5sOCpFHdBy4RmirfvWWZsbTLhQRiN42r6iKEz7RcV+tx9wthejbtjH6tm3M/fpG3DphGQDgmu7NDI+595OVAID83OGOy5OUKNZ+o1EWMQGVKM4qnrDAl+48Kq1tN7pQ4gn93eely0uC6VXj3fLT62+hLpRQNR3fj8EWnlDg3693786pMyXlskWIaeJdyfJGZibA0Cc5f9thfLfuQJVnLEu87YVFGD9/h5zGbeAJF4roLb96oukPbRdZrFBUXIYEQlCnpju6htLf9tHfQgJg7Oer0aReSpVjxFTk0QkB4PaJgULKqbq+Jsplec+kFfhwdM+K30f9ezGKistxZ99WSObkWlqy4ygenLwa8/48sMp35ElUSQkhLQghPxNC8gghmwghDwU/b0gImU0I+TX4fwNHJITyi4li0fYj6PrsLHR+5kccO1MqWxyFQ0zfUIAJC3dJW98JbfW0bhYrymU5J6+wyu9nSwO7NXnekTdmb8OR06WOZi0186opB/AopfR8ABcDGEsI6QRgHIC5lNJ2AOYGf3cElQBeDMt2Hav4+XBRiURJKnHCAnfDusaHC3Zie2FR1OMmLtyFbYeiH+cFvlmzL+oxsmZcTiyW10wOqNeScud2DkdV4JTSAkrp6uDPRQDyADQDMBLApOBhkwBc45CMUuv4yR/q4nDja1Kmsi0u8zmS4pRSihem52Hk24uiHvv895sx4k3n0xaLGGParYxkkIle8zhbWnUNi2frmiumTKYC10MIyQbQHcAyAJmU0gIgoOQBZHCXTmvXlaol9nDjREc/nkWr8o5PzcT9n63ifl3tO50pNZdkqZRn7psqGQjDfuw4kdoSqb+X7DiKTk//iAW/Hnakb2lrdz4Hv5RpBU4ISQXwFYCHKaWnLJw3hhCykhCy8vDhwywySq3gEU+LaPoXpRvcDEB4pS3yRTPbgUpAEeoaOIOJRynSTfnVamNXisjxtiI/4DK8bcLyyoRaHNvX9JYT7pmKNswcRAhJRkB5f04p/Tr48SFCSNPg35sCKAx3LqV0PKU0h1Kak56eziSkGy3DWER/n2W9uJ77bhMe+3KdTg53vEhyf9iCB77gU5dT+04yurUbhtLiHcb7OtxiOPBAs8Cd3PBrJgqFAJgAII9S+rruT9MAjA7+PBrAVP7iBZAZRhit6aKSspjc2h6qN7cXFmH/CedzsHy0KB//W1Vpoek7f6mDvsRovPfLDny/vkBa+07gBmUeiuzsGDxfIFqqAif9+mYs8L4AbgMwiBCyNvhvGIBcAFcQQn4FcEXwd2dwsQvlD5+vxnXvLq4SCuVV9K6qL1ftrfK3wa/PR9/cnwRLhCrT/7/P3AIAno/KMDucRc4+RA8xI9eUyEXMcN+ZZ/Oa8edkrc+o0eWU0oUwfr6X8xUnPG5exNSKO5SV+4GakoWxid4Put8lGQ/DWURerwtamZXPuX7t91OU+vxISU4MKejgjrGk5XwJhVLvp4ku9/lBUek5cPKd5Imt9FIXMU3aSy4ZFzFHuM7vErc4M2b7lJ3v+cL0PHR8ambk6k0u7LR5BafQ8amZQtJnOPX1++T+hE5Pz3S8HcAjCtyF/UxhE7+f4l9zjCsLfbZ0N7YdKgqr6jyuv3VJnZxj8vJARaHQdQPDNl0yxrRdiz9tCRsT4QkOF5WgzCdmOdYdCS+iIHMRM57RisE6sZFk8Y6jeGPONsO/P/ntRgDAqicHV/ubWyJTnMbOtzQTUeTGUSU7eZnXupY3LHCJbYt+oH4/xeLtR6QoKf2LctmuY+jw5ExMXXsAvxae5t5WmcnCxOHuguxBLoKjp0uQV2B6uwUTv+oWg92izJ1a71u1+1i1z8KtB3gtjNETCjyefCifLduNWz5chpkbxVRl0aO/zSfPBRYK1+w5IVwOPbGorI1Ki+kZ+q8FGPEW+8xHu3To3dO3aXYnqEicet7XvbukWhistpHHy3hCgctcxDQLr0iZ/COBgsV7jp1FOc8t1MxID8yt/pFgkV6akcf1epsLgpFLEVIR200mRioiIMzfLJ+fSu9zhkUfOHC6uGqo77yt1XeG57wwB3uOshcNF12mzxMK3M1hhLzRLKSXf9iCtk/8ILbtMJ/JNoDdsIg5fv5Ortd766ftXK/HA0II+r/yMzoJKBQeCSddh2bcI2dLfej/6s/Mbfxxyhrmc1nwhAL3ggUeC++YcFN6J6a0eQWn8OXKvdEPhFEYIR+ZDpw4hw8X8FXOboPC/Atv/4lzFVErExbuwt5j7JYoK7J3YtpF9G5dT0ShRHOBvzQjDz4/xVMjOokRyAH8fhpWWZ4T6KcMN9OJNM1n5ap/LTB9bDg/8JWdm3CR4y9fbQAAXH3hecismxLlaH6ImNUYDRmj2az+08NFJfjb95vx+bLd+OnRy3iLFhEn1zxkzyadwBMWeLTdY+Pn78SEhbsESeMMt09cHtZlslXgtvFwt/mbNfuFtR+OI6er+4IbpfLd8lou2OwT2Vqo0jITD6Ap0aJi8ekhnHwUSoFLQnQQypTle0xVS+HJwu1HAFS3nBJdFoHz1txfHc3tIIJQKzQm48qNwlBMoN0OOa5LuT5wr+ENBS7YwTzu6w0Y+s/ANN/04LbRN9bvO2H4twRBT8jnp6YslNdmb8PszeJDHJ3EC/rb6ktGGzFmXRJVNv5AS3dLAq49gS/s0yXOuQy98Jyt4gkFLsMSsDqttvN2f2Xm1oqfQw3uRIMvv2YP3xS2bR6fgRdNhsuVOuAXl4noWHMWi591LYLlLL14v/twGVo/PoOpbY2N+0/ijMlsnedKvZ/VUyQxsYjpFGU+P96UHPJl5EJZvusYumc1ECxNANlRQbyb17+rv1oVvfCuDFhfMpRSyy8MfcWgJTuNiy+Y4VypDyPeWoh2Gam2ruNFRLjmPGGBy4oD/3LlPtNFBOw8q0gvKKMFXJk2sJ3n4UZ/s16mR3XVgJxrr/Jnv9+6gjVDxUYeg8+rHa97pprLhMeo0+p5mk3HMCcvkMTKCaPt5R/4bsgyS9ynk42Ekwrh8W82OHZtI0IHmMvWMAHYkyk3WJRBJpsOnKzyu8w12daPz8Bz322OetzyXceY+jrL8OC6SO2i9/Wi7fZmE6w8/J+1juW18YQCj+RfdotB5xIxhGDnnfLJ4t3c5GDhXKkPx6sVhJD79D5enB/1mNsnLsd/TW5+AipfstTk6oz+paxVUXdL8YdYYO3eE45c1xMKPBLD3jS/KSSUP05eg8+XyVUoQOSBMu6r9QIlMYcT45rFV8tCuCyIoi1w1gXvrQfNZ4WsiCJkaKo8xhapYxlPLGJGYstB9njtaesOYNq6A/hd75a25bCjfCLpwxX54aNNyiQW+HUib0Crv87ALb2zzLVuo/lwp4qPQqn+2cmzZThVHLlUXHG59RC7aht5TJzzyH/XBo7l8JhluCHdiFNdzBMWuFvcJCIwO2Zem21cDMFpnIpC+WLZHmcuHAU39K85eYew/XBkC9uKnJWLmObi+/VsOhDw1+5jqIt64mwpxs/fUWHQTN8gNjeImxDRrTyhwEXCakm/OJ19hZvV0skrOIXhby7AaZMxtryINd+oG3KOP/rlOtz50YqIx0xevgf/WWHtJWd2K72dmayex7/ZgJdmbMGyXd7Pte0FPKHA5Q+v6HzNK2eIBd34yswt2HTgFJbvsre6bnWnHav69vspzpW5r4jAyXNl2CYw54wdtARcZuE1diLVL9Vz6lzAmDAbfhsLHDhxDtOjZCF0ahu/JxS4SGQYYzLs2TmbD1Xs5txsMcSJ1QCXmXBsRf4xzNtaGHYY/e7DZRjyxnxH2t199Ey16BGR2Qh5LQx/ulT+Yr8oBnXMsHT8Jbk/YewXq7EomM9IJJ5YxBSpVGVY+3qXxPu/mM9PbUfWez5ZCQDIzx1uuXQXqwLnlVkxKdG63fHb95YAANY9M6Ta37T+NeIt9ogmI0a+swgnzpbhhpwWle0J6GUVYYQM2QjDIXr3rUwnXXIiW+uhJdtEoCzwENy4UzAaoneqsrbHS8okG9pk8Ou/GP5t437+my1OVIs5F0tod2bt3mYVv5uWR2YwLqCu3XsC2eOmY7kH/PieUOA+kxXMebwBj50ttX0Nq7ioz5vCTYPUKnZrTbIi3jCojEKpKgfb1RJMPnReX7PU58d7v+xAmY0and+tO8B03qFTgT7yy7ZCS+dt3F91h6+I6bwnFPixM+asmJM2rZ2Fvx5Brxfn2rpGPMAahcJL8XtwklQFoS7BUAuc8Tqi39lT1x5A7g9bbIWW2o0usjrT/GSJ8TrBE99sxAkHjENPKHBRYV73fbpSSDu8qPbGF4RsA9yLifn1XVhoRZ5qcrC1bvalrR1W5vMje9x0prb0nLGRXtZtL/ofNvLPo+8JBW42uY6Xp/YsHDlt/Y3u91PbhXxZ7zMvX/2Js2X4aNEuT61XCHegVCxiVm2Z1Riy+szPcqrlaucR273nXtAnURU4IWQiIaSQELJR99mzhJD9hJC1wX/DnBTS56GB6nbm5B3CCzY2HQHAJMkJqV79cSue+24z1jiUIMgJ9IpUxOJYZUWeqp+z5n1JZoj8kY1dteEB/W3KAv8YwNAwn79BKe0W/GevZEcUfB5JrnPPpBWWp42UUqwX6ArhsZFmTt4hpvN4WzRy88F4hapjZ7vJvNyhHA1TXFrP4aISZI+bjgW/io+FNsJtMzQnXghRFTildD4AqfE0Zi1w2VMeLRl9NHYdOYOpawM7N6es2Cs0MiKWtsF76buIrURPURjsU7x02KkoFep3HTnDp6EQWJXw6ZJyzN1iLYqkGjb7l4g8MHY28jxACLkdwEoAj1JKw6bNI4SMATAGALKyzGWbC0VkUVURXPnGfJT6/Fiz54RrS3gp+CLSGFyyozK1gtdHDut9s5ObyEuwOrbeBdAGQDcABQBeMzqQUjqeUppDKc1JT09nasxqgWG3o5WZ+nhxPopEJ6IK+X3gP+YxX2vboSKUMKQ45YWHDPBqVYCcRJ92VtSLI9RSNltCLep1Gc87dU7uBqpwONFfmRQ4pfQQpdRHKfUD+ABAL75iVcW0C8WCl2lFvvt3WYmAdep75HQJhrwxH098szH6wUG8pHB5M+rfi6W0Kyvk8s255pJfRWNHlBS7hnDoa2YucUDC9nk9TAqcENJU9+soAOZHMQP3XNqK+zW13BixwLR1B/Dkt+ay1PFSokVBn+hKSy9Cvho8jt8Hpnlo8loh7Ti1HjF1LdtuSlHc/MFSqe2bCSOcDGAJgA6EkH2EkLsBvEII2UAIWQ9gIIA/OSlk96wGTl4emw84U3BUFN+s2Y/PlhrvWCs8VYzComKubbphhT+eLXqz8EogNuKCptEPchE8uoaZ/lV4ynwAghM5i6IuYlJKbw7z8QTuknCAdUAPe3MBpoy5mK8wYaCU4jMJVWd6vRRID5CfO5x7J7JiefFXuOI1+PEzpWhQp4bwdmVTw2Nx4DxmBGbGimwjwltPxUFYykdZZcP+k3jqW0e9TXGFjMGj1Yt0M05YevLnW+5E9iQwphT45OV7kD1uOlNSK575jv/w+aqwn7uhSolMi0F2Z+fBsTPis1W6ATe4zKzAYzzLtq7NEFMK/KNF+QDY0sqaTZdphhkbwiet8UKHMIsbhnMM3U7X44bnbQUefcNMDiYrrpoyk2mxrRBTClyDRVEmJBCc4ZSAxwg37BzkJcHoicstX88FX98+MfElrBPNAHfSQt977Kxj147EF8v5rlc5MQP3pAK/s282/nh5O8O/s2RcE1Eyyg1D//efr+ZyHRFrBtFwwwsxXpi27oA0N8prs7ZaPqeQQ3oK3ikulAIP8szVnZGeahwJ8J6FupIaPF0oRthVOF7XV7wX11bvPo5v1ohNRUAATF27HwUn5b/ARDNlxd7oBzmA2Y3YZ0vL8enS3Sjz+bFYl07ADtEqAlnp0Y1Sa9oTJgyeKGocjkjKcB1DmlG7ObLNYNfK73xeXUfqNnqV57/fDAC4rH2GsNC+Mp8fD01Zi+xGtTHvzwOFtOkk9Wsnm67bKTIdgJ5oM+pynx8PfLEGB08VY+3eE5i1iV/hhA8X7MLvL2tjfICFMZ1Zl78C94wF/vGdPav8HskaZdk+vHrPCcvnWOU3by+ydX40C3bk29aqy3PBBbOCaLlyztqo6hKK1tTBU+Y2Rrl51pRW05r9FmmzmJN8v74gYprmbYdOY+amg1gbNNx4prQ9Vcwvp4oTHijPKPDLOmRU+T2SMvNYxBM31u2LbCE5Me13sX6qQFtw5YHmB3Yi1tosu46cwezNbDnZ9dzSmy07KAAs23mUaabrBF4JjXVCLXlGgYcSyR0Rqwo8gQCpFq0mPbLzv8gaaCvyw2Y6ZmLLwcDWdJlKY9Br83DvJxzqtxL2sXLj+KUY+Y69GSUAdM+qj7EDI7goTODmWY4eKQUd3Eqkh3a2tBzZ46bjs6VyS3/xpmWjOkhJTmQ+/+BJvvlQ3EIkl9mg1+Y50qaIRW8jeBkoPKMiWBfoKQX+fGVHzHlkADdZRGLle/dvz5ZOOxIeVuDGN+54cFFGxMKkKP42sjNevrYr3JaSwlIuFAflMGLnYWcqxZj9Lk5+Z7thfSyhdjM3hq8yY1eWthmpts6PV1ymDsxjZmCEW9wqOHkOY3hMPwVzRacmqFMzCTWT2C1wJzCzW81xZIhgogNSSpmLCJvB8N6bfGskEGJZ8d7/GZ99BBqJIjZgRODbsX1tnS/bfeNZBW5mChuub77641bMsrgA9MaNF+LZqzvhi3t6WzqPJwnBJxUajROJMp8fu49WWqBOdLZdR87gf6v2mcpB49TGG7fp754vzsGYT1Y6qrx5kJRAbN+74jIfPlq0CzeOZ8uL/c8bu9mUwN6Ccr1ayTbblouH48CjH8OyIzOU9pmpGNW9ue3r2EV7YbVONz/VfHbaJnyuS19b5nNGo0xYuAsAsPf4WdSrXc+RNtxGpJfR4aISzNp8yPGdi3avzsOPX1Rcjjdmb2M6t3aNRLRoWNu2DPGMZy1wUQrcLZgdbPpp9aLt/OJhI5FXENhcJGs6KeMxnzRRc5FnPLITsK6n/DdkR2Yx42IoL/eJnX5n9yUrO52DpxT4S6O6YsLoHADWXSgDXv0Zd3+8wntp1YIkmuwo+45XJv6R3blCcUocO3Ufn726E0dJqnLa4YLVdlPbsirQ//tqfZXfWddBvry/D9N5vGjRsJbU9nngKQV+S+8sXH5+punj9d1q99GzmLulkL9QgiAeeFKyNrfY8TXf0Zd/vVUNp22F3sFKSxrT1h3AsTOlFp4C4SIk61Pv2KRuld8/vD1HWPtDOmXi69/bW8B0Ax5QC+FhXcS021+b1kuxeQU2WPyVotVpqIgvTt+M5bsCRY9/2FBQka+dN24tNiBSroMni/HHyWtw36fmI6wiGeBPDj/f9HWipTIwS5dmbOsnLDO7fu3TkZ5mPzeJ7DmuZxW4mYeWnGj/9vIagyxVgvSYdaFMWbEX2eOmc83/wUJJuQ8fLNiFG94P7P7klcY2HJQCd3y0HEP/Od+xNsIR7T6LfK9om3IKLGzWSiAErdLrVPz+xb2VUVb1a4uv+ynS48ejqRe+34yjkis0eVaBm7FIC04WY9cRZzZyWOXC52fZOt9s53533g4A1qpl80Iv4yP/WSesXT+lmLf1cMU291W7j2PjfuuZ81hCygpPlWDmxoMoDJPc6vhZd5df81GKj+/sFfZvMixL9jbZzzTzjp216WDYPEIfBqOvZOJZBZ7TsoGp4wb+Y56tdtwSySJz6zYL+mRL7Z/4wdG2Qmfw1727GCPesp6ZcVDHjOgHhVDq8+P+z1bhpg+qx0E/991my9djhaV7fLFsDxrWqYHerRoGrsGgCLl2S4Fd3MoC7phPV+EaDnlfnMCzCjyjrhxftCw1yhQxwEnYmknmusneY5VWSrmu/l9plKT4oTw8uHq1pTsuycZLo7qGPZ6Xr5nldmkJ/91QocgqDWobzzhkmC3RjJTscdPx7Zr9AIBnpm5E56dnAmB7iZh1SWocOlWCdXtPIHvcdPz16w3WG3QIzypwUdSxkf2PJ1b1N88B+PoN3Uwdd+8nKyummrx3IQ6/oCl+mxN+Q5VRW/qQSlMwKALtPVVa7sf/Vu2TvqC67/g500nLFo0bBMA9kbVmbv/78wP5jSYt2W2rhi2L0p+yIrApbjLnWpl2iDsFbnWAhbMGZWA1pptSymRRtmhYq1rnTk0x/xI7etq+39cosZHR9zF6plf9c4GldllcCHoX22NfrsMSTqW8rKK/BeNMWoi1a1R9rvrnLsUHbiqyrPqzZpGVZUZb7tBOZjvEhQK3U8hgUEfzcedugqWr3XFJNibd2Qu7Xh6O1o3rRD8hDHZ9ovm5wy0vJhp91yKLG2lYZPeFKBRWq/Chy9sJTez06BXtI/7dKVXVrL7x5hmRL41w9/qmni2qfaafzbgxt01cKPA+L//E7VqidjdaSVplhFVZn/1N54pcK/pTZbsFNIy+j92F5ruCm3kIArlvAKCJyTUWntkYRd7nBy/XzSzd8XhNvUDD3SKWMdm6cfVZXu51F1T77Jugzx0AzugMgl+2HbbcphN4WoFb2WygsXTnMQck4Y++hJxZZcITbUHp+ouau2V8G1pofpt1Ce7smw0AuKFnC3w7ti/mPNIfcx81V2AgVIGzvt4JAR4fZr0/84QY/CysfTMulJDe+M7P26PKGnrZhX8ZiK7NrW8amqkrlsyzTJ8dPK3A7+nX2vI5ZovRGnFefeeV6fCuTav8vvTxyy1fg1J7g/D92y7C7X1a4pXrLnCFhRbpu+gHtZ/BIm7RsDbyc4ejZ3ZD1K6RhLYZaaYXr1naCwelbP1ZBNmNxGQMZLHAX/1xa9RzfnlsYJXfmzcQlwGxUR1nN0RF7aWEkIkARgAopJR2CX7WEMB/AGQDyAdwA6WUX+FBF/PytV0x+HVnd/xd2q6xo9c3Q+v0VDw/skuVz2okJUQtw2UnH8rUKMn1jQa4flC/8/N25vZZCN1G7rFwfQDhk4HpP+G5ByGSm8hMK6FrDkD0e57VqDbeuaUHatVIQIsQ5V3fZj7waDi9rmHGAv8YwNCQz8YBmEspbQdgbvD3uMDpsML7+rcOu5giE22Am6mh+MmSfOZ2LmxRHwBwftNAkqMPQpIbGU2x9XG5i3aITeEaqlC8qMA1nF7fWfPUFbbbZy2RN/yCphjUMRPtMtOqfN4o1X4+FJlEVeCU0vkAQh3HIwFMCv48CcA1fMVyL05n3GudXofTQKLc0plaWVubEswV3bwBe6rOxqk1kZ87HFd0ysRFJnbcbmDYNm+FSDIsDsm5ftfH/Mr1nScpcRpQ1RqO9Pj/73/rI/y1Kg2iuBM8/O6TBqs5mUkpLQAASmkBIcRwDzIhZAyAMQCQlZXF2Fz8YOYF8dbN3fHg5DURj/l5y2FLiY0iwRIc4ZLAFceZtHi3Y9fOqJuCA5yeYST0z2ruowNQeKoE+0+c0/3d+GH+xDFFs9fSRZjlq99fgto1nKll6/giJqV0PKU0h1Kak56e7nRzjlM/wvZjI7YGkyzx4rIO0e/jizPyuLZpFb0CMEOPrPphP7/n0kCIXxsLpeR4E0mtWE0TYAUr0Ud2ilpoEBK4z33aNAq5thjM6u8OTzqXW2f5E9YDBiJBSGAGp7kFecNqgR8ihDQNWt9NAXi3UoIBr99wIYrLqg/OlOREZNatiUOnSjD4/EzMyYteIPlKK2lOTXRi0ZV2RAxgo7WFq7o2RX7ucNPXccLF5WbL0OenOHqmBANenSdbFFPw6EslIWsxPGd7GWny3FYssFrg0wCMDv48GsBUPuK4h2t7NMctvcO7fLTAAzPj+pyNfA1uweoGk+Iy739nUXRpFrDMtDt8ZefKnb852dH9/8fPluK5afayHkZ7ujwVZKTZo4vfk67FTBjhZACXAWhMCNkH4BkAuQD+Swi5G8AeAL91Uki3oSk0MxFCZiI39Jjpw6L7udXx+7Jk941IEhOIrd2YocrxzZu74+jpUpT5/MhqWBsvTI98Lx+asgaLtkfPv/K/+/vg+veWRDzGMNcMxznY8yO74M6+rcK6IpMT2OxJNy+3OB30EFWBU0pvNvgTX2eRh9DGq5mpdSKHqkCy0ZRMxyZpFUUTImF3s5TriPAIu5xXF+v2sUfBVMzmgr/XTErEeRHyhYRiRnkDkcNfNV97SnL4hTaeFnhyYgLah4TyaSQkEOTnDseN7y/Bsl3e2DEtG3fkSvUYmsUVTYHP2FBguvCEhpmxYnaqObxrU0zfUGCp/fCY+74aLJsXemU3tHyOG+CVCsVpKzLSIlrudV1x+fkZzDUpeWN1zWHn4dO22vvx4f6ur55khFLgDDw8uB2e+25zxNCgJ7/dgM+W7qni0zSDfjq++qkrLOcB11OLc+iS2XE1Y8PB6AeFMHZgW8vnuAG7yazckCgsLSUZ1/YIn2tdBlZ94XdPshd736FJ+BkBD5z263s6F4pTTBlzccS/39m3FfJzhyMp0fj2fbY0kPTdajhdmS4srWGdGmGLy5r1q/1v1T5LbRuhTb95VPE2IoHTlmOe/lozuKXknpOI/oqyFzPHDmzD7VpNHN6MpRR4GC5u3Sj6QQ7RijEPt5Nc2rYx/n5dVzw5vJNsUaJiVtl8eX8ffHJX+IK+oUTSJ2bWBCIRB/rfMk4v/EVjTD9+Cjw0HQRvlAJ3GCudcc4jA9CvXfRNOqItFEIIbuyZVcUl086gag4LN+bwy/1idvGrTXoq+rePfK//OMh5t47IGcPKJwcznSfazZMke+HfZvMZuplqY4dzrSgFboM/XdEOdaOUG7OibI1KiYUiy2rTD2Se/vW/X189kb4beGRIB0ubiMygr0gzrGsTfHp3b1zZORN3XpLNtZ1wOK1MePG3kCyYXuM3F54nrC2lwG2QkZaC8Q5PkdyE/sXh9am/LBtPH6Hz799dhMy6KXj/tpyoiZ5kIvpRt2hYm1sO8kEdDdM0GcJjhut0HnANFYViEzdvs3aS8+qnOJ4F0ElkPTaRdS/DUSs5EcNCCoZEom/bRthRyJbC1Q7fPXgpTpeU2y6H+O/f9bB8Do8n9POfLxOyC1tZ4ArTaFZ38wa1kJbibCJ8p7GyNsFT2d8VTM4li7y/DcVrN1xo+vhMSblB0lKS0bQee0piDaPNSU5TNyUZmQJKIXpegU+8IwejujeT1n60MLL1NnbpGSE6VC603TiddHChYZiwULcjq78BgQgoViaMZnNvik4WZwfPK/BBHTOlht7J8AXXiBB/7iRUl0LA6z5wWU5wyR4U02gL1kTys7ajSy8/39omuoo22ZsUjucVOICokSBOIsM6ibSByEm02YaXOrgRsowsrxh3lS9rdyeLindiYhHz1otbIiGB4Ompm2xf68v7+1g7IY56t/ZVCSFcXlwzH+6HwlMltq/DghU9ynNjCaVAm/Q62MFY25EFo7TIkah4WXP66lraXC/glZcsECMWeFJiAm7vk83lWj0tJlVqmymvUoxd/jK0I2b/qb/p47VFmbs5LcR1bFI36mYap7Di52zJKaRN49uxfbFo3CCu14xEY4aQtoqXNfi4UP4zxqJhJBHZO0GtEBMKXIP3pgszeK2Ch572malo3sC8ckqtmYT83OG49eKWnp95JFlwRmdxVOApNRKRlpJcZUNPJP7xW/MRIzypcKEkADwedqR0tgp2YkqBK6yRQEjFdNFL00YeJFtYR7jn0tZ47jedubR7mcUZx/UX2c8SyKJ+K6Or7FvgX9zb294FGJj76ADmc700FpQCj2d0HdXqIPW4AW5pQ02NpASM5rTVnSVEbdnj9mqnsChg7RQeUTOXtGEPBWShVeM6Uotgi0QpcA6w7PYK5U+D23OQxBrnSn3M1oYb8liL5tO7eyGzrvh8Ipl1U5CfOxzv33YR0/ksC85Ut4jpZBrhaLC88DxkQNtGKXAOWK26E46HBrdjPreDQYmqaJSUV2719dK0URb92qVj2eODcf8AfulG3Yo+5n/SXb3wmiRfvAy8NBaUAtfBnEVM8gPv2YrtBUJp5Yp7vLlQvIjIkmdaTvyhXZogs24KruPgi2dBxtBSUSge5c2buzOdJ/uBs3oz/JTd2nhsSAe2E12A3ST7sraWm41cCYWlf3Rokob83OHC/ddckKx/ReYJimsFXq8Wnxudkiz3NrKqE0opc19v0dBeaF2/dmyKgcczu6IT2xZrDVaXlcI6MtwZdtpMSiC4/7LW/ISJQlwr8K//cAn6trVfPk1mZr6berJXs/HiOuQvf74MC/8yUKoMMpOnscDrMc977DLL5yz9q70IGhbs6nw75991aSvUTBKXATGuFXi9Wsn4/J7IBYzNIiqBeyi5111gw4VCPZV5DQDq165hafORRq3kRIzu05LLZi+v3TNeZDMkjXO6qG84ZD6fVMEbluJagSdyfNByjVm21mskJTjuLuQxw+HBqB7N8JzHS3Wx4sWZlh1sW+A29MJ9A8S5T4A4V+AJLsntyZJsSA/rAM1p2dBxH+NdfeUWMNCQocRGXGC+8o0iPHKiUNjo07qRUPcJECPZCFnhWd7KzsaWl0Z1tXzOGzdeiJX5x4NtW2/zozt6IqtRbVBKcWXnTNzcy95LRKFwC/Hk4YprBZ6cGHjS7916EbIb24uq6N2qEWZuOmj5PCtJlfSM6t4co7oHYnOZdtpVVNcheP827xVmTkogKPe72zfgFl+5jLDHr/9wCc6UlNu+jox7yLw7WcJ9tuVCIYTkE0I2EELWEkJW8hJKFNp0Z2iXJujYxF6+4n/e1I3pPJFFkW/v01JYWxpGX8/uwEy1WMRjaJcmttozok268cLe1Rc0xbLHL8d/7/NOKlUeJBCgR1YD9GsnJ1WwXVj7pgw3HQ8f+EBKaTdKqSvMuKdHdJLSLmvxVB7622zH0af0FNXZRPsEjRjgUN7x85uGf/Hn5w7HkM6BXYy9WlnLMc+LihQPDj7rhrroqzH9Awt4rTkmkqrPEPcve2OdSGJuEfNq1u3wkuBhgWvjs0mUKtgy1mwvaeOOKBTZTL73Yky+l0/Iqhle++2FGGxzw1I0OjZJww05lfsQtPh4ntFdz420nsa3e1Z9bu1bgdWIs4NdBU4BzCKErCKEjAl3ACFkDCFkJSFk5eHDh202Fx2XuB1Nw0OpatZ0NB+c3jIRZYG7xQ8smz5tGqGPwJcZNfiZJwmE4LEhlVk0y32BlpIS+T3ztJRk0znRe2YHZhwsSt8qd/bNdrwNM9hdxOxLKT1ACMkAMJsQsoVSOl9/AKV0PIDxAJCTk+O42vCauuBhgfv8fgDRLR+nFlma1ktBwcliS+fwSmOgCE9KcgKuv6g5pq8vwB2ccpnr6dqsHl6+tmuVAttlwX7oVNHtBBLI32PEl/dfwq2twednYk7eoWqfN62XgiXB3aUfLcrn1h4rtu40pfRA8P9CAN8A6MVDKDt4zeLj4R8tLgsMHCN/rIbPX/kzT1U+dWxfU8fVTAp0tz9f2QEvXCNuU82cR9irs0TDjf3tr1d1xLAuTdE4tSa+e/BSnMeYBCsSn97dq1p2RM1f3cMhF4bIe/3h6PBLeiKDDszArMAJIXUIIWnazwCGANjISzBWjG5v7RruWEzTM6B9Ot6+xX4xiBEXBjaMPBul7JdTRRgyDHzvmuV3bY+quUPu6ttKqAXeNsO56ixuLGxx34A2jm9SC7dQ2Do9Fd89cCkeH3Y+57YCuGHfXST9fU138etvdizwTAALCSHrACwHMJ1SOpOPWOyEe0N+cU9vrHhicJXPmjfgb5VY5fymdVGLw4tlxAXnIT93eNQMgVX8ogIUj/ZCef2GbsjPHc51feIdDi8+hQ0MnmXX5vUs1Rs1Q2V5N7EaPNxO2oYRch5p+zJEwnynKaU7KaUXBv91ppS+yFMwZsI844QE4srFTdEytWJIRuQEPHzxfds2xkd39uQgjYIFGeOJ585pM7x9Sw/k5w7HxDsq3Sl2c8nzJubCCMN1rMQE4jrfFaCv/C2Ga3s0k5Y1EajMi8LLQhvYIQNpgrO/aSx/4nIs+L/KtLbM1ZwssuKJwYYV1xvUTkbn8+xtSDOLjNEkawwP6lgZjpkZJVRXNLGnwMN8luhSC5x1Gz0rBATds+zX72Tl/4Z2RH7ucK5T7A3PXVkx1b2ys7Nxz3oy0lLQomHtioU1US/j9LSaaGngKlv91BX47oFLhcghY/HWDT5wtxF7CjxMx0oKY4F3a1FfkERVee/WSt/t2IFthbatvwXuW3qzjwwLzU2LmIQQYRk2Q1upkeS8KjH6bu0yUjHz4X6Ot+9GYi6ZVbhH3Kx+rWqD+x+SqmwP7VK5MFK7htjbT+C9jU5mCKdCG9apgWNnSoXLIgK3hS7+9OgA1BUQVWS0zyH3ugts5zLyKrGnwMM840apNeHX7QDIaljbkW2vgzpmwOen+GVb+B2nTsXHmkU/8F1kOHJD/+zr104WqsBF3k43qG+9QcQz90kkwr241j0zREhI6oOD2qK03G/4d1lRUbGnwA26twijZeIdPfHlyr2GClzrgP/+XQ/mCuN2IADOb5KG2ZsPIaNuTeHtO43+2X96d2889e1G3NtPUIUUgRrcDQZ4gkDnq2ZshC6dDOmUKWw/waNDOkT8+3BJxTtiT4E7lL7ULGbG8bCuch42IcBDg9tjQIcM9JC4mOkU+vDEZvVrYeId4sIMReaCdoMLhWfCKje36XbiRoELI8I4li0aIQSJBLioJR/lPbpPS8NdmDJIqykvv0osuqQiITomG3BPCUQ3EXtRKBHU5O8vayNQkupIf7lw5rmRXYRH0oQlqDz7tmssT4Q4U+BeqpQTy8SeAo/wkP8QVOA39mxhfJBNZJRVisZtF4uvxCMDZaDFNqGRZG5R6LI2kwGxqMAj/C0tJRnbX7yqQpHHC8+P7IztL14lWwzHcMNLU/Su2lh+nkaE+sDdUnln7TNDpLUdewo8yms5KTHB0emfG6fShBDHcjRrfPV7frmYWZExoGU9bqefpxGz/tQff7+uq5S2Q33gbrHAZawHaMTeImbI76J3aLlQfwuB18IoC258aYpg4V8GglKg3ys/C2uzfWYa2memCWsPqJxhVbPAXaLAZRJ7Clz3ULu1qI8OgjtbJNwy5VNwhFb5TyjNG0ROHxxrVLPA1XiKbRfKt2P7uiJmtgKBoqh82YpYQy1SVyfmFLhsIk3nrxa4W0vWzjCZyHhXu2EBNV5wky3mFmJSgXdrUR//uqmbbDGq8PGdPXFrjIfz/SopMkKmD1xre0D7dCnpEeKC4D2ulm1SskIf0781br04S6oMMecDBwKuE1kYWWQ1kpyNfnEDvEtpWUXm3W1YpwYWjRuE7HHTJUoRm5QHE9GF9i/Zo4l37U8WYtICV8QXMt0Y8RoBIxJf8CaHhuvFukFkhpi0wN2IWjGPbbSnO+KCpq6KfIoFfL6AAq8RtMDTaiahqKQctR1ICe01lAJXeB4t1bsMiyzU+n9bUvRP7Rqxq8w0F0pSYuD5Du3SBE3rpeBuUamCXYxS4JwJnVJ3bJKGLQeLhK2gLxo3COdKfWIacwnUYIotpu3A/zJn8z881A+NUuUVq3aaIZ0zMSfvEDo0ScO8rYeRQAgeiZKfO15QPnCHEW0VNqtfC20zxFRICcfoPuIjbfwVUQrCm9Yhr/Hzm9ZFRpp70vry5oacFtjyt6Fo1agOABW6qUcpcM4Yda148YA/c3Vn4W1qiaRk5ItWqkQMKcmJFWGEauG4EqXAHcZNVctjFV/QBJdRlV5DBUQIIHiP/WpIVaAUuMMUFZcDgCNFlN2IlN2QEl0o2nNNUvu8HafCAlfzngrUIqZD9GvXGH3bNsZ/VuwFAKSlxMetlhEJUuFCkdD287/pjOxGtXFZhwzhbccb2tNVk9pKlAXOm2DvympYG/cPaINyvx8AkCSyjHeccWXnJgCA7MZ1hLfdoE4NPDqkg9Sc0PGCNoSUW7KS+DALJaAZg9omhMRENcCd4vY+LXH9Rc1RR2JpK4Xz9GuXjga1k3GPiv+uwJZZSAgZSgjZSgjZTggZx0soLxNqG5RpmxCUheYYhBClvOOAxqk1sebpIejSrJ5sUVwDc68nhCQCeAfAFQD2AVhBCJlGKd3MSzgvo22dH9W9GcbP34m6KcmSJYoP5jzSX3pSLYVCFHbMll4AtlNKdwIAIWQKgJEA4lqBa75uTYmMG9oRD13eDrVieKtzKG0zUjFG0jS3bYbKQ6KIH+wo8GYA9up+3wegd+hBhJAxAMYAQFaW3Ny5Irj+oubYffQMHhjUFkBgc4ms6f1nd/fG0TMlwtud88gA4W0qFPGIHc0SzqlbbXmYUjoewHgAyMnJifnl4xpJCfirC/IEA8Cl7RrLFkGhUDiIHWfhPgAtdL83B3DAnjgKhUKhMIsdBb4CQDtCSCtCSA0ANwGYxkcshUKhUESD2YVCKS0nhDwA4EcAiQAmUko3cZNMoVAoFBGxtbpGKZ0BYAYnWRQKhUJhARUwq1AoFB5FKXCFQqHwKEqBKxQKhUdRClyhUCg8ChGZmpEQchjAbsbTGwM4wlEcp1By8sMLMgJKTp54QUZAvJwtKaXpoR8KVeB2IISspJTmyJYjGkpOfnhBRkDJyRMvyAi4R07lQlEoFAqPohS4QqFQeBQvKfDxsgUwiZKTH16QEVBy8sQLMgIukdMzPnCFQqFQVMVLFrhCoVAodCgFrlAoFB7FEwrcTcWTCSH5hJANhJC1hJCVwc8aEkJmE0J+Df7fQHf8X4NybyWEXOmgXBMJIYWEkI26zyzLRQi5KPj9thNC3iSEcK3GbCDns4SQ/cF7upYQMkymnISQFoSQnwkheYSQTYSQh4Kfu+p+RpDTNfeTEJJCCFlOCFkXlPG54Oduu5dGcrrmXoaFUurqfwikqt0BoDWAGgDWAegkUZ58AI1DPnsFwLjgz+MA/D34c6egvDUBtAp+j0SH5OoPoAeAjXbkArAcQB8EKi79AOAqAXI+C+CxMMdKkRNAUwA9gj+nAdgWlMVV9zOCnK65n8HrpQZ/TgawDMDFLryXRnK65l6G++cFC7yieDKltBSAVjzZTYwEMCn48yQA1+g+n0IpLaGU7gKwHYHvwx1K6XwAx+zIRQhpCqAupXQJDfTET3TnOCmnEVLkpJQWUEpXB38uApCHQA1YV93PCHIaIVxOGuB08Nfk4D8K991LIzmNkDaG9HhBgYcrnhypkzoNBTCLELKKBAo2A0AmpbQACAwqABnBz2XLblWuZsGfQz8XwQOEkPVBF4s2nZYuJyEkG0B3BCwy197PEDkBF91PQkgiIWQtgEIAsymlrryXBnICLrqXoXhBgZsqniyQvpTSHgCuAjCWENI/wrFuk13DSC5Z8r4LoA2AbgAKALwW/FyqnISQVABfAXiYUnoq0qEG8siS01X3k1Lqo5R2Q6Bubi9CSJcIh0u7lwZyuupehuIFBe6q4smU0gPB/wsBfIOAS+RQcOqE4P+FwcNly25Vrn3Bn0M/dxRK6aHg4PED+ACVbiZpchJCkhFQip9TSr8Ofuy6+xlOTjfez6BcJwDMAzAULryX4eR0673U8IICd03xZEJIHUJImvYzgCEANgblGR08bDSAqcGfpwG4iRBSkxDSCkA7BBY4RGFJruBUtogQcnFw5fx23TmOoQ3kIKMQuKfS5AxecwKAPErp67o/uep+GsnppvtJCEknhNQP/lwLwGAAW+C+exlWTjfdy7A4tTrK8x+AYQissO8A8IREOVojsPK8DsAmTRYAjQDMBfBr8P+GunOeCMq9FQ6uRgOYjMAUrwwBK+BuFrkA5CDQSXcAeBvB3boOy/kpgA0A1iMwMJrKlBPApQhMe9cDWBv8N8xt9zOCnK65nwAuALAmKMtGAE+zjhmH76WRnK65l+H+qa30CoVC4VG84EJRKBQKRRiUAlcoFAqPohS4QqFQeBSlwBUKhcKjKAWuUCgUHkUpcIVCofAoSoErFAqFR/l/bC51JybMv88AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.plot(cp.asnumpy(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The raw signal seems to have multiple components, yet it's noisy, and it's not straightforward to detect any outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transfer the data to the Frequency domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a very efficient algorithm to calculate the DFT called FFT. We quote from -MULGREW, 1997, DSP- \"The FFT is based on the observation that there are many symmetries in the DFT matrix and that many multiplications were being needlessly repeated.\"  The equation changes into: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\begin{aligned} X(k) &=\\sum_{n=0}^{N-1} x(n) W_{M}^{k n}, \\quad k=0,1, \\ldots, N-1 \\\\ &=\\sum_{n \\text { even }} x(n) W_{N}^{k n}+\\sum_{n \\text { odd }} x(n) W_{M}^{k x} \\\\ &=\\sum_{m=0}^{(M / 2)-1} x(2 m) W_{M}^{2 m k}+\\sum_{m=0}^{(M / 2)-1} x(2 m+1) W_{M}^{k(2 m+1)} \\end{aligned}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we are computing the one dimensional discrete Fourier Transform of a real-valued signal, we can use a function called `rfft` in CuPy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "rft= cp.fft.rfft(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting the first 15 frequencies to zero was enough to clean our signal when we applied the inverse the Fast Fourier Transform. This number needs to be adjusted with other datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_num=15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rft[freq_num:] = 0\n",
    "y_smooth = cp.fft.irfft(rft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([16.13184904, 16.15809865, 16.18304895, ..., 16.04538145,\n",
       "       16.07548236, 16.10430765])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_smooth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f23602dd050>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABWHUlEQVR4nO29eZRk2V3f+b2x70uulVlrV6/qVXS3hAQIhIQGWdjIBsYj2R5kj+bo4DEwtscLHHzAPjMHM9izMMfY4zbSiO2IYUAYjYdFQgiETEuiutWreq/uqspaco19j3h3/njvvojMjOW9e383IiPqfs7p01WRWXlfvIz3e7/3vb/f98c45zAYDAbD/BGY9QEYDAaDQQ4TwA0Gg2FOMQHcYDAY5hQTwA0Gg2FOMQHcYDAY5pTQNBdbWVnhFy5cmOaSBoPBMPc89dRTe5zz1aOvTzWAX7hwAZcuXZrmkgaDwTD3MMauDHvdSCgGg8Ewp5gAbjAYDHOKCeAGg8Ewp5gAbjAYDHOKCeAGg8Ewp5gAbjAYDHOKCeAGg8Ewp8xNAOec47NPb+GpK4VZH4rBYDCcCOYmgP/209fxD3/zWXz0ia9ip9yc9eEYDAbDzJmbAP6bl64hwIB2z8Lnnr0x68MxGAyGmTMXAbzS7OAv3jrAj373XbhnPYU/e21v1odkMBgMM2cuAvgL18vgHHj0fB6PnV/CN64WYFlmFJzBMK80Oz184lcu4eOf/gs0O71ZH87cMjGAM8Y+xRjbYYy9cOT1H2OMvcIYe5Ex9vP6DhF4bqsIAHjodBZvP5tFudnF1YO6ziUNBoNG/uM3ruPz39zGF1/ewW89tTXrw9FKt2eh2upq+dleMvBPA/jg4AuMse8G8GEAD3POHwDwr+kPrc+tchNnl+JYTkVx11oaAPDGblXnkgaDQSO/+8wN3L2Wwl1rKfze8zdnfThaeelmBQ/+zB/ij765Tf6zJwZwzvmXARwcefnvAvg5znnL+Z4d8iMb4Gf+ygP4o3/4XQCAu9ZSAIDXd0wANxjmkXbXwtNXC3jP3av4rntWcelKAe2uNbX1r+zX8Df+w1fxm39xbSrrCbXgdD5O/rNlNfB7ALyHMfY1xtifMsbeQXlQw4iGggCAbDyM1XTUBHCDYU554UYJra6Fd96RxyNnc2h3Lby6XZna+j/3+y/jz9/Yxz/73RdQanS0r3floAYAOLuUIP/ZsgE8BCAP4F0A/jGA32SMsWHfyBj7BGPsEmPs0u7uruRyh7mwnMC1gtHADYZ55JVbdrB+YDOLR85kAQDPXy9NZe1mp4cvvbKDR85k0e5a+OOX6WWNo1zdr2M5GUEqSj8/RzaAbwH4LLf5OgALwMqwb+ScP8E5f5xz/vjq6rGJQFJs5uK4XmyQ/CyDwTBd3tipIhYO4HQujjP5BCKhAN7cq01l7aevFtDsWPjx99+NdCyEv3hLf2f31YM6zi3TZ9+AfAD/jwDeBwCMsXsARABMrTh7MxfHrVITvSmVEhbrbfzil17HZbNxajAo8/puFRdXUggEGIIBhgvLiakF8Jdu2tn/I2dzeORMDs9cLWpf88p+Hec0yCeAtzLCzwB4EsC9jLEtxtjHAXwKwEWntPA3AHyMcz61wuzNXBydHsdetTWV9X7mcy/iX/3hK/iRX3sKU3ybBsNC8vpO1S1GAIALy8mpBfCXb5axkopgJRXFg6ezeG2ngm5P3wZqz+K4WWrgbH5GAZxz/lHO+QbnPMw5P8M5/yTnvM05/1uc8wc5549yzv9Yy9GN4EzO3s3dKuiXUSrNDn7v+ZtIR0N4dbs6Na3OYFhEmp0erhcbuLiadF+7YzWJq/v1qTxRv7JdwX2nMgCAi6tJdHpcaxzZr7ZgcWA9E9Xy8+eiE/Mom04AvzEFHfypKwV0ehw/+wMPAYBp4zcYFLhVaoJz4MxARnp+KYl2z8K2ZpM6zjneGMj+73RuIpf39Emj22VbJVjLxLT8/LkM4Bs5+2TcLOkP4JfeKiAYYHjffWu4uJJ0u0KnRc/iptXYoJU/eOEmPv/iramsdcO5Zjez/YC24fz5luYAXmp0UGv3cMapx764Ygfyy7v65Judiv2e1k0A75OOhhANBbBXbWtf64UbJdy9lkIyGsKDp7N4fmt6EoplcfxX//5JPPY/fsEtvTIYKPnjl7fxI7/2ND7xq0/hP7+u/+nyZtEOaOIpGgBOiQBe0hvAhVQiAng+GUE6FtIqobgZeNpIKC6MMaymo9it6N/EvLxbw53OI9e9p9K4UWqipsnX4Ch/+touLl0poNbu4Zf+7PJU1jTcXvzKk1ewkooilwjj03/+lvb1hOx5akgGflNzABelx6dzfflmM6u3JFlk4KsmgB9mJRXVXoXS7PSwVajjzlU7gN+xYmtm09ox/9NXdhELB/B9D23gj17aNhUwBlLq7S6+8toefvDR0/j+Rzbxldf2tLe03yg1sZyMIBYOuq9l42HEwgHc0iyJXncy7cGW9o1cTKsUu11uYSUVQTioJ9TObQCfRgZ+Zb8Oi/c3O0QAf2t/OgH8q5f38Y4LS3jP3Sso1Dt4a990nxroeG6rhK7F8a0Xl/Btd66g0elp3+O5WWq4e1gCxhg2snHcmEIGHg8HkU+E3dc2snFX1tHBbqWJ1bQe/RswAXwsonFHZOAXlp0MXOOmh6DTs3B5t4b7NzN4+7kcAOCZa2YeqIEOMV/2W87m8ej5HAD9Le03i01sZI+bOp3KxLRr4NcLDZzOxzHo+rGZjWG/1tZWKLBdbmkrIQTmOYCnojiot9HRWIQvNjeECU08EsRKKurupOvkyn4N7Z6Fe9fTuHM1hVCA4bVt0wlqoOOF6yVcWE4gn4xgLR3DSiqKb94oa13zRqlxqAJFsJbRn5DdLDddvV2w4Wym6rp57FZaWE2ZAH6MlXQUnAMHNX2VKDdLTSQjQWRifROajWxM+2YLALxyyw7W955KIxwM4NxyQmu5k+H2443dquuvDwD3b2bwosYA3uz0UGl2h27oraSi2Ne8p7VXaWHlSDDVuYHKOcdBvY1lE8CPI+5qOu/at8oNnMrGDj1yncrqf9QDjss3F1dSZoiFgYyexfHWft3d3wGAu9dSuLxX1bZZvu8kW0eDKAAspyKotXtotPVIGZzb1hsrqcih18XNREdBRL3dQ7trYSkZnvzNksxvAE/bvwidlSg3S8f1uo1sbCodoDdKDayk+rv1d64m7U1VMwvUQMD1QgPtruUmCIBt09zsWNjRlBSJDHtYRrqS1BdIAaDa6qLVtY5l/+JmomNdoQ7kE5EJ3ynP3AbwnHNSinV9huw3i81D9aqAnYGXm13tteDXi02cHmh2OJ2Po92zpmbgdTuyX23ht5/aQr09nTr/WSKe5gY9Sc47m/RvaSqT3a+KDPx4QFt2XtvXJInuVYdn/7l4GMEA03JdFer2mktJE8CPseQEcF0aeLdnYadyfNNj08nIdbf9Xi/UD3WriXV1l1qdBHYrLXz51d2pPm1wzvHxX76E/+H/eRb/7HdemPwP5pzLTpC+uHrYFRAArmgaGL7rBMnhEor9mi4dfG/E2oEAw1Iy4t5cKBGxKWcy8ONk4mEwZnt162DXcREbloEDett+Oee4UWweCuCidnYa8s0saXV7+Ov//kn88Ke+jn/7J69Pbd1nt0p45loRkVAAn3v2hvYNtVmzVagjGTlcE72ZiyEUYLiiqc9BBMnlIRm4yMp1BFLA3sC01xm+gWoy8CkTDDBk42EUNEkoYlf6aAauUzMTFOsdNDq9QwH89BQdGGfJF1/awZt7NYQCDJ/8yptay0QH+aNvbiMYYPj0334HuhZfeNfJm8UmNnKHa6JDwQBO5+O4eqDnM7ZXbSERCSIROT5abNnRwHd1Z+Dp4TePXS0ZuB2blkwGPpylRAQHmjLwnfJwFzGRKeg00up7NvTXzsbDiIeDuKGxa+wk8PkXbyGfCOMXPvItKNQ7eOZacSrrPnWlgPs3MvjWi8tIx0L42pv7U1l3VtwsNY4lJ4DdULOt6elyv3q8jE8QjwSRjAS1ZeC71TYYGx5MdZUwFmptBAMM6Rj9LEzBXAfwXCKsTUIZtemRiYURCjAc1DSWL7rZfz8DZ4xhU7Nvw0ngycv7eM/dq/iOu1cQYJiKQ17P4njmWhGPnsshGGB4+9kcnpui6+QsuFFquvsqg2xkY9oa1faq7aHyiWA5FcW+putqr9rCUiKC0BBPkpVUBHvVFnn55EG9jXwijEBg6Lx3EuY6gC8lIyjU9EgoIhM4WgKkc9PDXbsmHvcO3zxOZWPaTe9nyU6lie1yCw+fySIbD+PO1RRemMIEpK1CHY1ODw+ctiekv20jg9d2qlpHbc2SVreH3UrrmCcJAJzKxrFdbmrZQN4bk4EDtjau67rar7ZG3jyWU1E0OxZqxDXohVpbawkhMOcBPJeIuBsF1BzUWsjGw4iEjp+i5VRUq4Qifvbykc2PFc3rzhoRrB8+kwNgdwbqbu0G+iV1oib63vU02l1raqZl02a7ZCcIg3ssgo1sDJ0e11LOt1dtDy0hFOQTERQbej7fxXpnZDWIu69FXP9+UGsjr3EDE5jzAJ5PhLUF8L1a+1gAFSwnI9oe9QA7+09FQ4csN+119VvoAnbL82e+fhWv70y381NMDL9/055Z+MBmBjdKTRQ02iUAwBs7dqAWXYn3nrLby19dUO+Z/lScIaZSmqqsLIvjoDY+A8/Fw9r6OkqNzqGKm0HEdU4dSwr1ttYNTGDeA3gygmbH0tJ+ezBGr9P5qAfYj5pDS63SEdTbPe2NJv/y917CT372eXzkia9Otanl8m4NpzIxpKL2po+YXXhZs//6G7tVLCcjboZ2btk2L7uqqR561oh9lGESSt8bhFYHLzY6sPjxp8pBsokwSpoCeLHeQS4+fO2sE9ipbx4HtY7JwMch9CUdWfh+rTWyfnM5GdVqorU/IlPpP+rpW7vW6uL/vnQNF1eS2Ku28Llnbmhb6yhX9ms4v9yflnJuyc6Irx7oD+CDHYmZWBi5RBjXFjaADy+RBQYycOK9FretfExAy8UjqLS6WkpHC/U2ciMycBFHKOUbzrmdgWv0QQFMAB/JfnW0i9hyKoJqq6vNQ3i/Oly+cQ28NMoof/baLpodCz/7Aw/hdC6OP355R9taR3lrv+4OzQCAs0txMAa8tac3kF4vNHB2YEo6AJxbSixsBr5baSE5oh57JRlFKMDI3flKDTu7zcZHBzQRYMX3UtHs9NDqWm6mfWxd55goCyIqrS56Fh+Z9VMx5wGc/sQDdllZod7GyhgNHNDr2zDU8GcKTURPXSkgEgrg0XN5fOc9K3jyjf2ptLRXmh3sVVuuHwcARENBbGbj2joDAft3vT2kIuPsUmJhM/C9avtYhZMgEGBYTkXIN/TKPgI4tZQhft6oihC3q5vwxuHl/VIw3wFc0+ZDsd6GxUe3wIrXdRT/9zd7hmvg9rr6JJRvXC3iwc0MIqEA3n42h0qrO5VqjCvOuLgLy4cz4fPLCW3eHIBdutiz+LGKjHNLCWwVGugtoPvjMF/sQXS0lnvLwCPO99JvJgL9TPsowQBDJhZGiTCOVJr23pHOJh5gzgO4+DCUm7R3bJFZj5JQxI2D+lEPGL/Zs6zZcrNncTx/vYS3n80DAB506qJ1j9kC4Ga7YvqR4JzmTFh0th6tyDidi6NrcXeq+CIxzBd7kJVUlPzpUlwr44ydRIDVlYGPklDs46K15RAZeMZk4KMRAZw6kO6PqMMW6Pqg2WsPb+IBgEgogGw8rC2AbxXqaHUt3OeU0d2znkY4yNzyPp0IzfVoJnwqG8Neta1tWrrwljkqoZxyLBS2y4tnajWpoWYlFSWXUMS1khmTkeqSUERGP06PziUipBKKycA9EA0FEAkGUG7QlrqJGu9RGbhbdqQhAxcblCMrYDSWMIpyvTucioxwMIBzSwm8uae/HvpWuYlIKHCsVldUSujKhN2SuuzxGweg13VyFnR6Fgr1zoQAHsFetU3aWl5qdJCKhoa2sgtEgKW+rlwNfExFiF2DTnddCVUgHZtxBs4Y+xRjbIcxdswkmTH2jxhjnDG2oufwJh4bMvEQuYRSmPALdzN/DdUv4sM2KoDnNXafipmbFwcqQe5YSeFNzXXYgB0oN46MrwPs1m7xdR3cKDaRioaOZYZrziTxRZNQRDnfqE1MwM7A2z0L5SZdYlRqdCZu6KVjITBGf12J63lcBp5P0DYRiQx83BMHBV4y8E8D+ODRFxljZwF8AMBV4mPyRSYeJpdQXP1qxN0zGgoiHg5q0cAnbfbkibW6Qd7cqyITCx26eVxcTeKtKYxyu1VqHnN+BPpShq5B0jdLx+eeAnY5XTDAFi4DFzNkV8dp4Gn6TfpSozNRDw44FtHkGXijjUgogFh4TPZPnBhVTkoGzjn/MoCDIV/63wD8EwAz3abPxMJuwKWi3OggGgoca2UfJEd8xxYUJ2QL2XhEmwPjm3s13LGaOhTM7lhJot21tDnUCW6Vj08/AvpShi4Tr91KC2tDstFAgGEtHV04DXzUZJpB+uWqdJ+zUqM9sgpkEB3t9KV6B7l4+NhN+tC6iTAqzS6ZgVm52UUsHBjqpUSJ1E9njH0/gOuc82c9fO8nGGOXGGOXdnd3ZZYbSyYeJn3UA2z9alK2oCNTAOxMJRIcnS1QP+oNcmW/fqyM74I7J1FfJQjnHLfKTTfbHiQTCyERCWrLwPdr7ZHBbD2zeO6Po2ySB9HRb+BFQgGALPFmImCXEU5yBcy5FW00saTS7GjPvgGJAM4YSwD4KQA/7eX7OedPcM4f55w/vrq66ne5iWTjOjLw7kTtKhvX49tQarSRTYzOFvLJCBqdHnkXqGVxbJebx6pA3ElAGjPwQr2Ddtc6Nr4OsPc5TmVj2qSMcTXRpxYygI+uchKI80EtoXgJ4HkNHv/FemdsCSFA31NSbnS1V6AAchn4nQDuAPAsY+wtAGcAPM0YO0V5YF7JxEL0AdxDBp5L0GvvwOQPuq5Sq71aC50ePyZjrGfti1nnKDdRCTIsAwdsC4Fd4rI2AGi0e6i1e0PHbAHAeiaqfXg1YFeG/OKXXsdvPbWlfa29SguxcADJyGh5MJ+wOxMpx4yVGpODKKBHEvVy88i6pcFEAbzZGbmHRonvWwTn/HkAa+LvThB/nHM+kyGCWWcTk3M+VuPyQ7kx2jtYkItHUGwUSdYbpOjodaMY9H8ZlrHKcrN4fAoQYG/Yrqaj7td1ILLc9RHvZyUdxUsafMHdbDQ5PBtdy8RQadqeN+P2Q1T5xS+9jv/9j16zjyUVwXvvXZvwL+QRNeDjrpVQMIClRIRMQml2emh2LE8ZeDoWcis4qKg0uxODqbjeqUqSK80TkoEzxj4D4EkA9zLGthhjH9d+VD7IxMPoWhwNQkmh3Ox6ysB1bWKOzcCF8Q7xY+Y4h7pNjWO2gL674rDNRABYSdIFk0Prjhl0C/RLOXU6T7a7Fn71ySt4772r2MzG8MmvvKltLcDW/Ef1NwxiT7uikhO8+4LYe1od0hr0crMzMZiKr1OVJHt5iqfASxXKRznnG5zzMOf8DOf8k0e+fmFW2TfQL/WjbOYpNzoTNfBMPIxW1yLXoic9aopMgfrm0W9oGeYRHdcqoeyJxqkRmfBKKopys4tWl/ZcT9rQc03LNHrPPHOtiP1aGx95x1n8lUc28dXL+24Jmg6K9Q6WPEgZ+WSE7MblxQdFkI6F0OlxNDs01SCcc1Rbk7PhfhyhOfd21n8CMvCTDnU7Pefck2amy/py0tqiuYg6gN8q2Z2QwxqINnIx3Cw1yYe+CvarbSQjQcRH6LIiY6TOhCeV1Il1dU5f+s+v7yHAgHdfXMF7711Dp8fxtcvDqnZpsH2xJ1ucUlY7+QngIpBS3cRq7R44n9zS3s/AaRJBOwk8ARn4SScTp330aXR66Fp8soQSp8+EOz0L1VZ3QseYHgfGGyM6IQG7EqXe7pFbFgjsgbPjW7sB+kEWwu9jVNerWFdnBv61N/fxwGYW2UQYj5zNIhhgeG6rqG290gSJTpBPRHBA9BnzFcCJDeq8NtTEwkFEQwGSDLzVtf3HT4QGftKhfvQRQWrypgftrrW9tvigj/7Fx8JBxMIB8lKrW6XGyCoQsVmqSwe3ddnRNy0R3PeIM+G9agvp2PHZowLXNlhTBs45xzdvlPHQGdv1MREJ4Z71NJ7Z0uP+2OlZqLS6nial55N2wxjFU5fbnOZBuqHOhKvOzxFj+sYh9HdV+kZWJgOfCLWEIn6BmTFBFOgHeEoJRTQweKmAoW6n3y63Rla1rKXt13W5IO5V2yP1b6A/iYjaIW+v1nZ/9jBS0RAioYC2DPxmqYlys4u3bWTc1x45k8XzW0UtclXf0tVbPXanZ+vHVOv6kVDIEjIfroB2SbL6+3V9UCbEEArmPoC7j1xkGfh4HxSB+EBQljx5/aDnNDQ7jLMYXXWqQ3TUYgO2hDLOn1pk59Qe1ZMGGzDGnAoYPQH85Vt2aeTbHPtewLbwLdQ7WqY9+cmE84Sb5eJz7SUjzRBfV+IG5CmAk2XgzvuNmgx8IhnqzQc3A/cawOkyYdHZOanhgbqJqN7uot7uzSSA2xOIxksoiYgtG5Fn4NXW2HUBW7450CShCJ/1ewcC+J1rKQDAGzv0Fr6uL7anTUy6EspSwy7jCwYm92nMSgMX30MRR1wZ9iSUEZ50QkG7q4xMQml4s4FMu7vl08/A7W41unXdARYjglkyYrsv6gjg5WYHXYuPlVAYY1qmxBTrnbFT0gFbB9c1+/TKfg1r6eih4HKXE8Bf36UP4H2jNG9lhADNZrnXNnqA/snWlwYeC6FCEEf6Nw0joXgiQ+iH4jUDj4QCiIYCqBBohAIhi3hp+6X0QBdDJEbpwYwxrKbp5yQCfV3dSyZMub7XclGdAzSu7Ndx7sgIuY1MDIlIEK9ryMALviQUuoaxUqPjaU0AiIeDCAUYaT02MF0JpWwCuD9EOz0FQsbwUsOZjoVJJZS+CfyEDJzYA31SBg7YMsqulgA+2R0PsLsxKQNprW2Xi07KRsWAXx2bitcO6jh3xP0xEGC4YyWpZYiGSBC8SCiiAqdQo9HAvWbgjDGkY3RDWiqtLhgDkhEvGbj9ZKv6u+5vYhoJxROZGF1AKzc7iIeDnnx87WlAdBl4tdW1x8RNWDsTC6Pe7qFD5F3sxSNal6GUl5sH4MwsJNy49SpXLScjaHUt1Nq0XaCtbg83y81jGTgAnM0ncL1AX7JZanQQYEDak5wQRoBNX0IB7MBHJaFUmh2kIiEEPOnvIbR7FlqK81fLTfumkfJw01BlMQJ4nM4Ap9zoei7/sTNwugBe8dDyC/TLk6jW3p8whxOw/UK0BPAJbfQC6klEJY9ygts4RayDbxUa4Bw4v3w8gJ/Jx52v02b9hXob2XjYUzALBBjZlBq/ATxN6DBabXaR8ihlUJUwlp35n17OsyqLEcBjYVRahCY0HgvwM7EQqYRSbXY9bbZkiUsn96rtsQ0tALCaiqFQ75Bl/YNrM4Zjw4yPQu2DXnQqMryYlgH01gVX9+0BGcMy8DP5OBqdnp5NWw/yiSCXCCtLKJxzlOr+jJ0yhImRH1dAqgoYL+6HVCxEAE8TFeAD/lzEqK0vqy1v2YKbKRDdPHarrbENLUC/lJB6Q69Qs7PCcdPKAfpAKm5+42wL7HXFpHTa9331wA7gZ4cGcPu1LWIZxasnt2CJIANvdiy0e9bE8zwIpQZuG1l5T8gAoKQYS7y4H1KxEAHc1sxoLCi9TOMRpKO0m5heM3Bxg6HS/fc91EPrqgX3Mu4KoPeAKXqsue9XY9Bm4NeLDURDgaE3zjNLtif7VoF2jF2h7m0upSCXUHck9NOFKaDNwDuerimgXxqsnoFPx8gKWJAAno6FYHGQbDTNMgOvtLpIeejecg28iJ469qqj50IKRICnLiX0qo/mCMvaxLrA5Jrovn0vvf/6qTHmYQB9Bu5XQllKqjsSygTwNOFUHj8SSta9rlQ18OkMcwAWJIBT+ieUfNhApp1qEKpJ1tWWt0evLJFWJ/CSgS9rGm5QrHurEaZs7QZs35lQgCExZrQYMDhqi9h7ptTE+gjzsHQsjFwiTJ6BlzzMhhxEOBKqPNl67W0YJBMP2WWeBNeV18IAYFCaVEuMKq3pDHMAFiWAEwU0zrnt4+u5CsX+PgrDH8CHhEJ4w+r0LBTqnYkZOPXQV0GxMRsJRWT+k8bwRUIBpKIh8vd9q9wc6f4IiIHKdE87wonQjxadS0TQ7lpK0678GGgJhJRBcV35mQ5P5as0rXFqwIIEcKr221q7B4t7zxYo237F5BAvm5iJSBDBACPJwEV53PKElvJ0NIRwkGnJwP1IKJRDBrxmo9Tj8zjndgAfM9N0PRNzZ4VSIBNIKfR/OQ2cRiLs9Cw0O5ZnDTwaCiASDChdz5xzU4XiF8r6zcGfNwmqTQ8AaHUtdHrc04eNMYZMLESyienVwpYxhjxRXbCg27NQaXY9BZVY2PZjoarH9jrYALCzf0oNvFjvoN21RkooALCeiZIGcD9OhIIcQQ28+Iz6kRSorquqjzZ6wLmu4moVMPV2Dz2LmwzcD+LDoZoJe/VBcdclzMD92F4Cwv9FfV1xYXuRMZaI29mF1ui1MoKymafU6HheN0fcRHTLCczjJJT1TAy7lRbZ/oqfNnpBnuCpp9zogHns/hRQTdkS15TXDBwQ7fQK79eH+yEFCxHAqSZKe53G01+XzpHQj2saQGdo1b+wvW0kUmbgfoMKZTt9sdH2nIFTt/G7ATw7et9hPRODxek80P04EQrEvodKDXzRKQrw05UofLSpEjI/wTQdV7OUneYwB2DRAjiVhOJzE5OiFrziM4CrZgqCog+NcolwUvmhtT0+1ueTYbpNzHrH840jnwi7x0rBdskO4OMlFPtrt0o0Moo4fl+dmHEaDdyPbAPQ7S35lVAAKEuTfvzHKViIAB4N2Yb/ZBKK5wycTkIRVgCefRuIjLR8ZeBJTX4kvjJh9fV7Fke52fUsleUc98eeReNNctMJymJU3TDWM3Z2TqWDu+V8Ehp4UVED97OBCdAlRn6sZAUZRYdRr/MEqFiIAA6IaRo0TQfeG3nEo556UBHZgncfFhoHxmLdrof2kvkvJaMo1ttkgazgU0KxNXD1DFz8vvxIKJzTec9sl5tYSUXGuk4KfZwugHt3IhREQvawFNUM3H8Ad8oIVTNwd1/Jn5GWSkJmNHBJKAaSin/v9Y7tDnUg3MT0pYETSSi5xOR6aABYSoRhcboWfr+6bD4RQanRgaV4A/Hahemum3Q284je926lhdUx2TdgD7AIMJDVghcbbeQSEd8OeblEREkD92tkBdANSxE3al+bmIrXVdlo4HJQZODlZgfJSBDhCcZKx9clDOA+qlBaXUvZnc9XOR1xN2bR5xNPLhGBxdU3q10fFM8SCm0T0V6tPXaIMwAEA/YUJMoM3M8GpkC1Bl4mAwdohqX4mUjvrhsNodW10OrKXVcVnzKsKgsTwDOKu8cAnC5MfyeeylLW/yYmjf5eqLc9SxhLxN2YpXobGY/DbgE6Yym/TS39JiKa931Qa431XhecysTcihVVij7b6AUqlUdibJ3MjSMTU9/jqba6CAcZoh6Gs7jrKpYkV5r+11RhcQI4wUBSP17gAipDK78fNir7ANvgyNt7XtKQgfupSxaZnKqE46fyBhgc6kAjoRxU254C+GqabgqSV8uCo6hk4HVnbJ1cBq5+XYk2ei/y4OC69r+VW7vslE36WVOFidGCMfYpxtgOY+yFgdf+FWPsZcbYc4yx32GM5bQepQcopAw/03gG16XwbBA+KF5/8VS+DfYjrr8MnCyA+7h5AHQBvOSzfNHNwAk08Ganh1q7N9F7BrBH3FHWgctkwipdqDJt9AIKCcWrt9Agql3d0/RBAbxl4J8G8MEjr30BwIOc84cBvArgJ4mPyzeqLbCAXAaeitJIKF59UATiOJWz0Xrbs5QgMjjKDDw7gwy85NMhT8yHpJBQRED2koEvp+y6e9VNW0BFQglLbxz73WsYhCYD9x9MlTNwH3bUFEwM4JzzLwM4OPLa5znn4h1+FcAZDcfmi0wsjLbipl5JQgNPx0LK5U6A/YHx4gUucL2LFdZuO8N6vWZmsXAQiQilH4m/AQOUGXg8HEQ0NN5KVhAIMGTjNCWMB1XvAXwlFUXP4srrdnoWqq2ulISSVdg49vukM0iaYG/J9tf3mYErSpMnMQOfxH8D4PdHfZEx9gnG2CXG2KXd3V2C5YZDsaln61f+Tn6KTAPv+POLIDDwEuVhfjrllpK2RzQFooTRK7Sykb+gkktElEdtAf0hzpOqUOzvccbYEU3F8dsRCahtHKtLKBQZ+HSbiMqNjmsFMA2UAjhj7KcAdAH8+qjv4Zw/wTl/nHP++OrqqspyY1G9c1oWR6XlvTtPkI6GUG13lR9zfUsoBJuY/cnsfgyOaNrpLct/hUIsbHfcqstG/gN4Nh4mkVAOXAnFmwYOAHuKG5kyQxUEKj7sZaUAHlIeluJ1QMog/SRBvgplWjXggEIAZ4x9DMBfBvA3OcUwSkUyisZS1XYX3IcXuCAdC4NzoNZW923w87gnvItVmpeKEplZPhkhaacvNzv2+fb5WJ+Nh90bjyx+h/sCjh8Kwfve9yWh2N+zqzjGrihxoxaIz4bMORdPeLIZOKA21EFGzkhFQmBMPgP3M0CCAqkAzhj7IIB/CuD7Oee0c58kUTW08usFLhBZs2olit8MnMK7uN8J6c9ilCITlXHHA5wAPiMJhWIy/X6tjXCQeZLqXAlF0cK3bxks08gjn4GXGh0EPdo0HEV1M5FzjqpEAA84xyuzt9Tt2XtK02riAbyVEX4GwJMA7mWMbTHGPg7g3wBIA/gCY+wZxtj/qfk4J6IqKbgmNL7LCJ0ATqHXSZQ8qejBBR9GVoJ8IkKyiSmT/QOOlKEYSGWaS2wJRT0DF008XspFs/EwQgGmPEjaPdc+btQCVQ3cy9i6YYhrQfZ6bnYsdC3uqzBAkJHs6vbr6U/BxJU45x8d8vInNRyLEqoSil8nQkEqSlMN0up6H/0kUPUu7mvg/gJ4udlFt2ch5MNy4Ch+XBAHycbDuF5U606U0cBzCXtTTfV9H9TanvRvwM4GKYZoyDgRClRKKEuNrpR8AqgbWvVtXeWyfxlpsp8EnqAMfF4gk1BkndMUJJSaTx8UgW3gpVaF4vcRl8rYqV+h4C8rVDUbanV7aHR6vm8cImNXlW/2qu2J80cHWUlF1TPwui1lyFicihJKmacPmbJcgaqEUlHIhmUtZcsKNw1ZFiaAqw76LUlq4BTexTKjnwDh/6Kmged8PuKKqgRVHVxmRiOgroHLlrb1p9OoBXA7A/cRwNMEAdyZPiTb3i3rh1Kqe596dBT3umrJ12MP/hw/yHrtq6wpy8IEcMaYUveWrA2kCLoqGrjsL97WwNWqUPxXY4huTBpHQBktutrqSpeXyT5piUCkqoMf1NpY9lADLlhJRrCnKKEUJNvoBVnJChxZJ0JAfVxhf0Th9Nr4ZWVYFRYmgANqm3ri38kW/qtIKDLG84C6fUDRZyck0M+YVbsDi4020tGQbz05525Wy53vflOLP+lGfH9JYQO12emh2ur6k1CcDFylUtceHycfVKQzcEknQoBAQlGQM2SlSXcepgngcqiMGSs37U5Ir9amgqRTN6qymVgV49QkqlBU7ANsMyl/gcy1lFXtDpT05hD/RlZGkfXncOdDKjx5+GniEaykImh1LaUEQQxzkEXGkVA0aslm4LFw0O5zkK3HlpQlAfvprNry35znd6YuBQsVwNNReQcz24nQ/4ctEGBIRdT8UFwvcMmuMdksRcbgqN+ZpyYlFHyYaA2i6ofidxqPwNX+FTTwAx9GVgKKWvBCTU1CkXEkrLa7sCQa4wZRkURVsuF0LARLojnPr6c/BQsVwDNx+bFq5ab/tltBStF4x5VQfGfgarWy9iOuv8wsHgkiGgqob2JKrA0MatFy68tm4OmY/aRVUnjfwtPEiw+KYFm00ytsZJZ8+q4fJZ8Io9buod31vu9QUnAiFKgEcJFQJaPeDMsGkS1JLjc7SESCSmWmflmsAK4wVk1mGo8gFQ2paeCKGbiMXtd2HstlsuClpLofirSEQpSB+/1du+V0Chn4vhOE/WXg9vfKBnDhRKiigWclKo9kz/MgKp7gFYVgKvai/MaSioQdtSoLFcBVHMxKDfmTr+pdXGl2EWBAPOwvW3AdCSXWFheYbHu1qoRSlB21pehIWGrYT1p+9zoAW3ZRqUIRN71lHxr4qpuBqz1xqG1i+u/GLDfU11XKwFvytq5Cw/a79rStZIEFC+CZuJ0J9yScAVVcxFKxsNIE7WrL3zQegesJLhHMRDWFXzMpQN0PxbK4r0ESg1Bk4LJBJaswHxKwJZRQgPn6nIn6c1kNvORaBqtIKP4zcFmrhEFUPMErEtN4+uvKJQnTHuYALFoAV2i/LStm4FUFDVzGtxgYzMAlnOIk67ABO6ioeIKLDS4ZDTwaUrOUVamMENNpZBGzMP3cqMPBAHKJsOsj7heV37MgJ5GByxilHUXlibqs4AooO1vAZOCKpCU39XqSXuDuulE1CaXa6kiXOwFy3sUqj9aq1qruBpdkdpaLRxTKCOW7A1UllH2P0+iPsqzgh1IgkFByUhm4vJWsYFYSiqwGrpIEyrJQATwj+XgtMnbZD5vyJqZPK1mB6wkukYG7ToQSGdKSU1YmO8RCNStUaaeXqbwR5BQG/AK2hOKnC1OwrOCHIo5XZpyaQEYDLzU6iIQCiIXlQ4wYGC4ricoHcJOBzwTZE99vgZW/Y6tMD5GZng0MeIJLaeBOEE3KbWLKzkkEBke5yQUV1QAu+6SVddwfZQIK4M+JcJCVVER6rJrKXEpBPBxEJBTwZeNbkvDZOYooq5UZliJ7TQFOE1Eo4Ou64pxLS6EqLFQAl9WEVUueRPZca8l1RFYkM3BAlE7KSSjBAPNdew70HQllK1FUBgwA9u9JRsrgnCttYuYUu0APfDoRCpaT8hl4oe44TkbkM0PGmC2b+ehCLSq27wNq7fSqk3H8XletroV2z5pqFyawYAE8K1lipmpCo+qcVpUY5uCuLWmvquJQ1ze0ki1tk/enBpxMWOI919s9dHpcYRNT3omx1e2h4tMHRbCSiqJY76Aj8YQnHCcDEmWTg/j1Q1HZLBb0Da3872nV2j2ljshMzJ/PUN9K1mTg0khLKIoeBiL4qmy4yH7Y/H7QBCoOdaqWsrLdkAJZCUW2jd5dNyHvhe620Utp4PL+MzJ2CcPw6wlebHR8e70fRfZ6ppiMk477q4BxhzkYDVyelOQYpv7JV5weIrGR2bM46u2evIQimY2qONSp+qEUG3aXXDTkv80ZsINJrd3znZHKeoEL3KEOEu9bVJH4aeIR9Lsx/QfwQr2ttIEp8J2BS9b5DyLrtS+uQ5WKEL+OhOIYTR24AqFgwB5I6rOsrqx48lOSHzRAfpiDQFoDV3CoczVwaQlFzVxJBAa/Ny7VzF9lwK/bhSlZhQJAqhbcdpxUDyr5pD8bgSKphOK3GsRx91TIhv1O5XHnCRgJRQ2Z7q1yowPG/JtJCVIKEopqtiBbhaISRFPREEIBJt2VWGq0pTpABbLdmKoVGeJ8yWygiuArWwcOyHVj2nXv6hm4KKH04kve7lqot3tKN2lg0KzNp4RCMBknHfNnTV12n+6MhKKEjKFV2dlElN3oke3cAuSNrPprh9GS8ASXNZMCnKqEpHxbuWoGLh/A1ZpLMvEwGJPTwPsSinwGLlOJUiDKwHPxMDo9e3NwEv2hGbPZxKSwdc3E/WbgZhOTBBlLWRUnQqAffGU0cNlhDgIZT/BOz0Kl1VVu7pAdbqD6eC3es99AKjuNR2APBg5LWcoeCB8UKcuEECLBgG8NvNmxBzjnJW4aR3H3PTzIZio+O4PEwgGEAsy/hCI54WqQdDSEZsfybKGruo8my8IF8HQs7Lucr6xoAxkP2wOVZTxYZIc5CGQ8wSkypJyCsVOx3nZ1dBlky0VF7XsyIrd5CjjTaSSrUPLJiNRTHmMMy6mIa0frFdVN20H81MBTrdufc+vf1hVQk1D6iZG3tcvNDsJBptR5KsPCBXB799hfILW789QaHVJROec02WEOAhl7VdXNPMBup5cJ4Jxzp8ljNhq4andgLh6Wqr7Zk2ziESxLdGMWCNroBSKL9/I7pzDQEsgYWlFMh0/71N+FD4rKZ0uGxQvgcQkNvNFVfvRJRUNSlrLqGrj/jR4Si9GkXCCrtXvoWpxGA/e5PkVzSTYRkZRQ5IysBMvJqO8MXLXjdRA/figUHuQCGUOrqqS//iAZn/p7uSlvhqfCwgVw8Qv3M8WbwsdX1jmNoowQ8JeBC+1arZTPe1XCIKL5R+XijoQCiIeDUhm4alNLXkFCEZuRMiynIr41cNWO10FEJYuX5q0ioXQjK6GkFbNhNwP3+DRvZ+DTrUABFjCAZ2JhtznGKxQ2kLYnuLwGnpT0qnAlFB8fcnGBqTxaLyUi6PS4743bvnyj9lifk/DmpsjAZS1l92tqEsqK40jo54ZZqKv/ngWuJ7iHjeuSKMsl2NCTlVBUXQH9auCVGQxzADwEcMbYpxhjO4yxFwZeW2KMfYEx9prz/7zew/SO3+L/Ts9Crd1TNqER1pd+EW308iWM/j3B3cyMoJnGbzCjeqyXaacv1mkklHKz48uRsNXtodLsKkkoK6kIWl3LUxmfoEgYwMPBANLRkCdHwlK9jUwsLDW27igyT7ZlAldAv7MFyk11GVYGLxn4pwF88MhrPwHgi5zzuwF80fn7iUAEYq8nXnw4VC9s6U1MBdtLwC61CgeZvwy83kGAqW3yiGDk19BK1UpWkJEI4FQZOOdykpWqBg7Alw5erLeVPbkHySW9PX0UFRwfj5KWuK5sCYUqA/choUy5iQfwEMA5518GcHDk5Q8D+GXnz78M4K/SHpY8fjVhqpKndExuqIPsMAcBY3Ztsa8qFKeNXsWhTratnGqDy28G3rM4yk216hdg4MnDx9qiC3NFoo1esCzhh2L7oNBVRnj1Q6F40hGIJ1s/0lFFwd1TkIqEwJj3OKIywk0F2VvzOuf8JgA4/18b9Y2MsU8wxi4xxi7t7u5KLucdvw5mVAE85bP1VlBRcCIUZOL+/FBUnAgFeUkJhep8+w3glWYHnKuXtsk4MbpOhBJGVoKVlEwG3iGRTwRZjyWU9uAKmnXTsRAsDl/SUaWlnoEHAnZpsJfrqtXtodmxFnMTk3P+BOf8cc7546urq7qX872ppzrMQZCOhtDuWmh1/bW0Vyke93w6p6k4EQqkJZR6G/FwEDGFEi/AfwAny/xlMvCqCODTzcApM2HAvnl5KaGkDeD+2+mrRJNxvNpyiGTxRG5ijmCbMbYBAM7/d+gOSQ2/EkqZTEJxLGUldsxpMnAfmmxd3onQXTMWRoD59wQvENw8APsJoN7ueb5hioCruraMpaxowFGpQllyDa28Z+BUVrKCfMJHBk60rt8n6v5oM/Vs2OsGqjtPYI4klM8B+Jjz548B+F2aw1HHbwcVmYQSlfNDURnmIPCtgRME0UCAeX6kPro2RVYofDa8BtJ+5Y2qBu5f+z+otRB0zpcs0VAQmVjIVzcm5WYiYL/3SRU4jbbtvyIzuGIYfj3Bmx0LXYvTZeAerivXSvYkbmIyxj4D4EkA9zLGthhjHwfwcwA+wBh7DcAHnL+fCNyBpD4lFIpNTEBiekhTbRMTcAy8fKxbrLelJ7MPkk9GcOAzAy811I3+Af+T0qkc8sTnxI/2f1CzM2HVsWYrPqbT25YF6k9ag+QSdgXOOOlKbNiqPG0MIgKx1883hRe4IBP3loG7wxxmkIFPfJec84+O+NL7iY+FDPvO6bH8p9lBJBhANKS2HZCSCOCWxVFtq+t1fjLwtlNLTNNeHfEtoRTrHdy1liJZ2/553tan8uewHQlDvvR3VR8UgW1o5e391pz5nxS/Z8HgOR+lcYuSSSrpxq9Vc3+wAoWEEka5WZm8ZmP+NPATTcZH+62wklUttUpHHQ3ch4RS7/TAubyRlSAT9+4J3q/DpsmCD3xaylJp4CIT9pqBUxh4CXI+b1x2Gz1BAE9GPU/lobAsOErOw1OPm4GTSSj+OyLtf0eQgXvVwGeYgS9kAE/7KKuzmztoNjwAf7vlVI97frKUkluNQeMR7SeQcc7taTxE8g3gIwNvtJGOhhAKqn/kcx438wT71ZaSD4rAjx8KReniUbw89VCv61ea7DsREmjgzlCHSTXo4umX4qbhl4UM4H7K6ii68wC5oQ6uEyFBBg54K50sEJXTAfA9lafuPNZTauBey/lUJhAdJRv3Z2i1TyahRFGot9H1MMx5v6ZeungUL8Os3QBOJKEkIrbXvtfEiGIivcBrDXq5afvMJxR85mVZ0ADufRxSudElc00D/GngFaIPm5/SSZE9UU0qb3YsNDw2WfRNtNTPdzwcRCQY8HwDoazI8FoPDdhTcSqtrlIXpmA1FQHn3mQjoZVTrCtwa+AnZOChACOryBBe+17Lc/sSCk0VCjD5urLtqENT9wIHFjWA+6jKKCmOUxNEQ3ZA8RPAKYavAoP+L5PXptSC+5UgXjcSaUr5APvCziXCKHrU4KkqbwB/U3n60+gpJBTv0+kPFIYojyITCyE4YZi1mDxEGcz8GFpRzMPsr+vND4XCjlqWhQzgaR8ZOJWEAtgyStXHOLf+h029CgXwmIE7m5gUcxL91kRTGv0DjgbvwR0PoPECF+ScLlDLgyOhyjDjo/iZTr9ftY2sKAKZgDGG5WQEe5XR6+8TNvEI7GoQf1UoFO/bqzEehTWFLAsZwDMxeyDppC49yzE4ojPe8Wd96Q40JnJO86qBhxTnQgpEdud1uDF1AM/62EwU49Ro1rWlDC/nu1+VQZeBe6kFF/7j1I/1a5kodirNkV+nbKMX+BnqUGl2kIqGiKxsvVXAlIjr7f2wmAHcoxVktd0F5zRyAgBfWh1AM7sP8OcJLpo7KC5s3xJKg05/F+t7qULpz+Gky8ABb808lFr0ig8/FLvyhT6orKVj2KmMvoHsVlpYSdNVvgDCUta7hEJVDeKOK5xwXVGVxsqwkAHc64aiKKmjqt9Mx/zNxVSdxiPw4wlOGsgkJRSqG6Zdxjj5PffncBLdOJLeK2AoM/BMLIxQgHnyQ7EzYdpACgBr6Si2y8PX55xjp9LEGnUAj4VQ8ShNUniB99f19mRbJPac8cNCBnCvmjCVE6EgFfU3/qna6iLplEmp4McTXHhEU+BnzBZgf9Bj4YCyE6Egm7AHDEyq06WcDQn4mw+5X20jGgqQSFaBAPPcjUnV/XmUtbTdTDSslLHS6qLZsTQEcO/XVYXIiRDwNnWq27NQbnZNBk6JV01YfJ0qI8z43MSksr0EvHuC22ZSNBd2OBhAOhbynIEf1DqkG1z5RATtnjVx/ilVG71AXKxe2un3qm2spKJkWrTXbswDxRmco1jLxMA5hppq7TiZ+VpGQwbucVA5hTmcQIyRG2eZLK45s4lJiFcJhcpKVpDyvYmpbmQl8OrPYZv8033Ylnw08xzUWmQudYD3Zp6+kRVRGaFo4/fgDLhfo9WivXRj1ttdUkfAQUR2vTNERhGbm2vpGOmaaWdQecODVQSlBg7Y1VrjnrTEZ5+iqkuGhQzg/iUUml+42MT0Ov6p7OyYU5BNRLwFcCI3QEEuEfFcCUKty4oniUmBlLz6Je5DA6/SVmWspCZn4O7GqQ4NPGMH5+3y8UqUXWdzU4cGDnhrkqsQjzazHTdH/577njMmgJPhVUIRgYdqAyIdC6NrcTQ7k1udATsDp8oWlhLhiYGs0bZHP1F+2PIe1hXYNcKU7njeqkFcAy+iJ62QIx15q0JpuQOJKZhUhw3oaaMXuBn4kEqUvoRCnYF79xmyp8MTZuATPt/U8pxfFjKAJyNBBNjkO3bBmdpN5WHgWsp61MGrhI97+WRkYiClGK57lCWPg24B+gzcNbSa0MwjzgvVJibgdGNOeN+cc+zX2qTnezkVRaPTQ709+rPtdmFqkFDEbM5hteA7lSaioQD5bMiMR0/wVreHdtcilVCWEpGxGjh1EuiXhQzgjDG7e2vCI26hRju1W3xwvdaCU264LCUiqLS6aHdHZ//9rkC6IJrzWMrX7PRQb/dI9WBXi56w/n6tjVQ0hGiIzmxoKRmdOB2n1u6h1bVI3/PqGA1a4GbCxFIGAERCASwnI0NLCXcqLaxl6DZsBV4lFMo2esEkDZy6wskvCxnAAfuXPumOXSCe2i0+OH5KnlTb6AVe7FVdpzjijcTqhBvHobUJH+tdc6UJgVRHd+BqKupqvqMQ9dqUN8yNrC1P3BqiQQtEcF3VEMABYCMXw41i49jrt0pNrBNvYALe582607WIZbpauzfSa79Y77hDPmbBwgZwL46EhRptAb6fuZiWxWk1cDElfkwA36vSjrsCvPty6wjg0VAQiUhw4mailgCejkxsaRcBnrIz8ZQI4KXRAXyn0sRSMkL6xDHI2XwC1wr1Y69vFRo4u5QgX8+rBu5WGxGVyQKDn+/haxcbbeQIBsLIsrgBPB7y0AJLe2H7mR5Sa9P5FgN9DW6cXkfpjHd0XS8yBkB78xDrT9T+NTS1rKSiOKi1xw74FZnwOmFd9Clng/BG6XgGPLiuDvlEcHYpga1C45CZV7dn4Va5idO5OPl6njuriRvzgL6v+ajrqkDoMy/DwgZw28HMg4sY4cn3V+5Eq9d5MZbar9mbthRdgQLRVj7uxmF/veV8P3UgjWBvFhJKOgqLj3/fOuqik9EQMrHQxAx8nbgSZJCz+TjaXQu7A08gt8pN9CyOM3n6AJ6MhMCYF1Mp2r4OYPIT5izb6IEFDuCZCe23lsXHDmeVQQRwLxKK+B6qRh43kI6RMvarbawQO9R5HS5Maas6yHIqir0xWjTn3A7gxBUZohpjnA6+XW4hHGSkjVMAsJGN4+a4AK45Az/jyCRbAzLKVsF+IjiTp5dQAgF7qMOkPa1+wxb9IOdR19VeRU/Hq1cWNoBn4+PLvMrNDixOW4Cf9LGJSTm7DxiQMsZKKLSdkIPrjrtx2Gu3nc0e2mC2koqMbWyptLpo9yzyi0xsEI7TwW1jpxi5PnoqGxuZgfcsjt1qS3MGbgfpawd9GacfwOkzcGByQgYMbGKSZuDju273qvTui35Y2AC+nIqM3T0Wmu1Sku6XHQ4GEAsH/GXglL4NE3xJbI9o2g+bF8MfQJhoRRAg8GkeZCUVxX61PXK4wkGVfrivWBeYEMDLLXJfEMCuRBmVge/XWuhZnFR3P4oI0tcO+hn4lf0aggGGjZyeG0cqOtkTvFjvIBkJIkwwuFowbo+n27NwUG+7n4VZsLgBPDl+80G8Tt0C63UaUH92H2HTwYRmHh2bebGwXQkyi41EwJZQuhYfaSOga/NUNOeMk1B0WKsCdga+V20NLd3U1Q05SCwcxOlcHK/vVt3XXt2u4PxyQlvli5dhKZTTtQTjDK0K9Q44t2eVzorFDeBifuAI4x8hr5CPf/JoPk81kX6QfGK8b8N+raWlvdped/obiUA/kI6SUXSULwL27y0WDkyQUPRIGaIWfJgfiXhNpwYOAPedSuPlmxX376/tVHH3Wkrbel48wanm2x5lJR09tGErcMtyTQZOj7hg9yZc2NQ7yF7HqlWJJtIPMs63od62vZp1fNjyyfBECWVP04SYVXczcdSTFv1wX8Du9l1JRUc6AzY7PRTrHU0ZuC1hDJNRrjsNNqc1adGC+zbSeGO3ila3h1a3hyv7ddyznta2nhdP8HJDz2Sc1XQUu0M6T0UANxKKBtzMbMQFJrTiHKEGDojBxpMDuNhRV53GM0g+Odq3QVcVCODUYk/IwHcqLXKbUWDypHZXQtHkCzJKQnGd+TRk4MM0aMFWoYFoKODe2HTxto0MuhbHa9tVvHKrgp7Fce8pnQF8cmJUbLTJJRTAfpoZ5v3SD+Czk1Bm0/85BcSFfTDqwnamdqcJJQwASEfD2KvUJn5fpdlBOhoi3dRbTdmPepzzY5UP4hFwJa0ngF/ZPx5MBNVWF/V2T8uGnjsnckQg3a+2EQ8HkSC8UQpW01FcHfG++zXg9O/5TD4OxoArQwL4tYO683W9nYEPn84BAJ66UkDX2UB+7Hxe23pib2nYZ1ugQwMHxBzQnWOvC1fIua1CYYz9A8bYi4yxFxhjn2GM6ds58UkyEkQ0FBiZge9WW1glnJQiSHmcoK1Dr1vLxNDuWkM39HYcbVSHJrueiWK73Bzpg76jUZfNJyIIsNGDfrfLTW0VGZvZ2MiOyBtF+z1vZOmljGgoiM1sHFf2jycKW4WGllrso5xbTuD8cgJ/+uouvvLaLs4uxbW8V0E6FkKnx9Ea47lTanS0+HKvZ6Kot3vHnqz3qi1EgvRJoB+kAzhj7DSAHwfwOOf8QQBBAB+hOjBVJmmUu5WWFrOfVNTbYONyo0sfwMd4NYu6YT0BPIZW1xppXbDjGv3Trx0IMMcZcHgGbjvk6ckrNnNxVJrdoR2/oi5alxZ9fjkx9Klnq1DXVot9lA+8bR1//PIOvvTKLj704IbWtURmParaqNmxve61ZOAZ4QB5WEbZrbawkqJtjPOLqgYeAhBnjIUAJADcUD8kOpaSo5s8dAXwjKOBj6pLFpQbHWSJJgEJxo272q7YXYHUVTfAgMHSCIc8N4BryoRXUpGRm5g7ZX1t5ZuO78fN4rDNxDqy8TBpldEg55eTuHpEQqm2uijUO1oMpYbx377nIlZSESwnI/g7336H1rVcq4gRey1lDT4oApF4HLXQ3S43sZ6dreggHcA559cB/GsAVwHcBFDinH/+6Pcxxj7BGLvEGLu0u7srf6QSjJvgvVfVlIHHQuAcqE+Y31dudsi7EtfHjLvaLttdgdSNNIPrjgzgmkvbVkeUeXHOsV1uYV3TuiKAD5NRrhcaWoydBOeXEziotQ9l/0JSOTsFCQWwb9x/9k/eh6/80/e5N3FduB2/I67nvhOhnk1M4PgQixvFpvsZmBUqEkoewIcB3AFgE0CSMfa3jn4f5/wJzvnjnPPHV1dX5Y9UguVkdGhVRrdnOZNS6C9sP97F1I977qPeiHFXujLgU2NuHID9tBMJBbQ83gL2DWR7SEldpWUP99X1vjedrsNh3tjXiw2tpXwXlpMAgDd3+zr46zt2Y83d6/rqsY8SjwQRJzRHG4WoIhrVb1DU0EYvEBLcYMUR5xw3ig1szmsGDuB7ALzJOd/lnHcAfBbAt9EcFg0rKduz+ejm2kGtbXdQadLAgcnOaWUNm5iJSAipaGhoydN2WY/ZPtC/cQwLooB9Q9GxYSzYzMawU2mi2zu8wbXj2rlqet/pGIIBdiyAc861Z+D3OSV7L98qu6+9tl1FMMDc4L5ITPL6KWjq6wBsWTQeDrob04DdhdnqWlo3br2gEsCvAngXYyzB7Cvz/QBeojksGlbT0aGbayJD1VEr25+LOToD7/Qs1No9bRsuQzVwjdUY0VAQ+UR4jAbe1JYFA8BGLg6L2zr/oXVd6UZPAA8GGE5lYocubMB+uqq1e1o3E88tJZCMBPHSQDfkq9sVXFhOIBJavPYO4eg4aoydO8hZQ002Ywznlg4PsRA37bmVUDjnXwPwWwCeBvC887OeIDouEoQud7N8OEMSeqmuTUxgvCOh+JqOMUzDmg4a7R7Kza7WDZf1TGzonEQAjg6tb23RWn7zSCa8XRGVN/puHqdzcbf7USBc+nQG8ECA4b6NDL55o5+Bv3yrorUbcpaEgrYENyoDP9DkeSM4u5Q41DjVD+DzK6GAc/4znPP7OOcPcs7/a875+BlTU2ZjRMux2yWnRUKZrIHrmN0nWEsfD6QiwGxoDOCnsrGhGriQE3RmKv3NxMPr3yrplVAA4MzS8Xrsy3u2Fn1xVa8W/cBmBt+8WUa3Z2Gv2sLVgzrefjandc1ZspQc7fWzX20jEQkiFtajx59bSuDqQd2VY0WZ6Nxm4PNAPzM7fGHfLDbBmJ6yNi/z+9ySJ+IqFMD+QN0qNQ+N+hKPfuc0lpetp2NDJZRCvYNGp6d1Q29UBn6tUEc+EXZ92nVw52oK2+XWoSaPN3aqCDC7UkQn77q4jGqri2e3inj6SgEA8KjGbshZs5SMjOysPqjp8doRnFuKo97uuVLN5b0qMrHQTIc5AAsewNfSUQQYcOtImdf1Yh2rqagW68uUh6k8OoznBeeWEmj3rEPZsHj001kffCYfx26ldcx/XUxt0SknpGNhpKOhY09a1w7qWm9agB3AAeDygLXqG7s1nFvSZ60qePfFZTAGfPnVPXzplR0kI0E8dDqrdc1Zkk9EcDBiZOB+rU3u+T7IOedmLJqnLu/WcHE1NdMmHmDBA3goGMBaOnbs0fp6Ud8jfSoyWQMXtbs6mg5EwBps8rh2UEcsrNfgSHzAjxosXRcdiZofNTeHaNFbhYY7/ksXd63ZFR9vHArgVTew6ySfjODdF5fxya+8ic8+fR3vf9u6NgnhJLCUHK+B68yG71q19xZeuWVvGk/rdzyJhQ7gALCRi+HmkQz8RrGp7ZFezO8bF8B1Z+DA4QB+9aCOs/mE1mxBrHu0vVsEVd3t3eeWE3hrr69F9yxbe9fd1HJuKYlggLk12O2uhcu7Ndw1pVrsH33fXai3u7A4x498151TWXNWLDl9HcM8d3T5zQvOLsWRiYXwwo0S9qstbJdbuGeK9fajWFg3QsFGNoaXb/VLrSyL43qxgQ/cv65tzVQ0hOoY83lR1qhDA9/I2bXJ1w4F8Ib29urzTu3xUYe8rUIDyUhQWxOP4OJqEn/6yi56FkcwwLBdbqLds3B2Se+NIxIK4O61FJ7bKgGw67LbPQuPnMlpXVfwbXeu4D/92HsQCTHctbaYFSiCpWQYbacEd9CigHPujAvUF8AZY3jwdBYvXi/h2a0iAOCRE7BhvPAZ+GY2juuFhutNslezR1HpfKSf5F1canQQceZnUhMOBrCZi+EtJxPuWRyXd6u4c1Vvc0c+Yft+HJVQrh3UcXZJb/YPAHeupNDuWa7mLrLx80v6m1oePZ/HM9eKsCyOZ51A/vCZ6WnR929mFj54A6Pb6WvtHtpdS2sGDgBvP5vDizfK+MI3txEKsBOx37DwAfzOtRRaXct9lBeP+Do3tyYNdbCtZEPagtq962m8fNOuD35rv4ZW18K9pzJa1hKIZoejJXXT0govrh7Wol/Ztp+6dA4ZEDx2Lo9Ks4tXdyr489f3cCoT0675346sjvAk6Q+u1hvAv/eBU+haHJ/5+jV8+10rWqubvLL4AdwJHuLCFlrlXVrn94XdUsFhlJt6ZvcJ7t/M4o3dKpqdnrvpct8UAtn55YSb+QNAq9vD1YO69uwf6P+eX9t2AvitCpaSkalMS/mOu1fAGPDZp6/jz17bw3fftzbz6oRF5JQ7C/RwKeFu1Q7ougcrPHwmi/feu4pIMID/7r0nY79h9rcQzdzpZmY1vPdeu0Y3GgpoLcDPJ8JDzfYFpbqeySGC+zcysLjdmffstSIiwYDWG5bg3lNp/MGLt1Bvd5GIhHBlvw6L209BusknIziTj7v65Es3y7hnfTplXuuZGL7jrhU88eXLAIAfeuy09jVvR06NcL0UDVunNDZsAfZT5ic/9g60uj0tE55kWPgMfDkVRT4RdjPv13aquLiaQlCDraogn4iMtL0EnB1zDaY7AtGN97XL+/jq5X08cjY7lfKy+zcy4M6NA+iXXE2r3Oqx83k8daWAWquLF26UtY74Osq/+P4H8Nj5PP7ed9+Jx84vTW3d24lsPIxIKHBssILoedAdwAHb/+akBG/gNgjgAHDfqQxeuF4C5xzPbhXx0Gm9evByMoJKq4v2iPFPhbrekqdT2RjetpHBL33lTTx3vYT33D0dG9/7N+3zKvw5ntsqIhIKTM2f4/ELS9gut/CZr19Fz+J45x3LU1kXsNvmf/vvfhv+8ffeN7U1bzcYs83Djmbg2+UmIqGAlon0J53bIoC/444lvHijhOe2SijWO3j0nN7MLD9megjnXHvNKgB89J1nsVtpgQH4wcfOaF1LcDoXRzYexvNOJcaz10p4YDMzNXe8771/HcEAw//0/72EXCKMd100mfCisZ6JuuMBBbccp83bcd/htgjg77pjCRYH/sX/+yIA4N136s3MRHAeNkyi0emh1bXcIK+Lv/mt5/Gzf+0h/NrHv3VqFRGMMbzr4hK+8voeaq0untkqar9ZDrKWieGH330eAPDj77tbeyu7YfqsZ2LHBpbcKjWnIp+cRE6OmKORd96xhM1sDE9fLeLRczm36UQX7vy+IQFcBHWdGjhga3V/41vPaV1jGN91zxr+8MVt/OKXXke7a+H9b1ub6vo//Zfvx9//nnu0Nw4ZZsN6JoYvvrQDzrmbcW+Xm3jgBNRkz4LbIgMPBQP4hY9+C77v4Q38yx94WPt6IoAPM58vOGY8ujPwWfF9D20gFQ3h3/7JGzi7FMc7L0xXxmCMmeC9wJzKxNDo9NxuZsviuFlqzny02ay4LTJwAHjHhSW8Y0rBxB3/NEQDFzP9lpKLGWSyiTB+/ocexq8+eQX/4AP3IBS8LXIEw5RwXQEPang4kcN2pYlW19L+VH1SuW0C+DRxxz8NKSXUObvvpPChhzbwoYc2Zn0YhgXkjhVnmPNeDQ+fybmd1bq9108qJj3SgDv+aVgGXptO26/BsIicW0qAMTuAA3Ab5hZxkLMXTADXhD09ZHgADzA9ToQGw6ITCwexmY0PBPA6wkGmdVzgScYEcE2MDOD1NnKJCAIaO0ENhkXm4moSl3ftAP7KrQouLCdv272W2/NdTwF7/NPxAL5XaU3FYMlgWFTu38zg5VtlNDs9PH+9dCJsXWeFCeCaWElFsDdkE3O70tI6Jd1gWHQePZdHp8fxhW9uY6fSwoMmgBuoWcvEsF9rodM77IeyW266vsYGg8E/j53PgzHgZz5nd1a/5+6VGR/R7DABXBOnMjFwDuwOtP1yzrFbbWEtbTJwg0GWlVQU771nFQe1Nu47lZ6KVfJJxdSBa+JU1s6yb5Wbrvd4od5Bp8exnjEZuMGgws/+wEP41FfexH/5+Nnb0sRKYAK4JkSWvT3gnCZGQZkM3GBQYyMbx0993/2zPoyZYyQUTfTHP/UDuBgFtWYycIPBQIBSAGeM5Rhjv8UYe5kx9hJj7N1UBzbvLCUiCAcZbg3M7xPBfN1k4AaDgQBVCeUXAPwB5/yHGGMRALenIcEQAgGGtXTsUAa+dVBHgAEbORPADQaDOtIBnDGWAfCdAP42AHDO2wBGD4K8DdnMxXC90HD/fq3QwEY2jvBt2jVmMBhoUYkkFwHsAvi/GGPfYIz9EmPsmKMMY+wTjLFLjLFLu7u7CsvNHxeWk3hzYDr91YM6zi2ZhxSDwUCDSgAPAXgUwL/jnH8LgBqAnzj6TZzzJzjnj3POH19dnc5w3ZPChZUkdistVFu2+fzVgzrOLk1nvJnBYFh8VAL4FoAtzvnXnL//FuyAbnC46HgXv7VXQ73dxW6lZTJwg8FAhrQGzjm/xRi7xhi7l3P+CoD3A/gm3aHNPxecAH55r4a201J/93p6lodkMBgWCNUqlB8D8OtOBcplAH9H/ZAWhztXU4iEAnh+q4hK056Fef9GZsZHZTAYFgWlAM45fwbA4zSHsnhEQgE8dDqLb1wtotToIBML4UzeaOAGg4EGU8+mmcfO5/HMtSL+03M38e13rdzWvg0Gg4EWE8A18+G3b6JrcdTbPXzfw2bQr8FgoMOYWWnmgc0sfv4HH8Z2uYkPPWgCuMFgoMME8Cnw199xdtaHYDAYFhAjoRgMBsOcYgK4wWAwzCkmgBsMBsOcYgK4wWAwzCkmgBsMBsOcYgK4wWAwzCkmgBsMBsOcYgK4wWAwzCmMcz69xRjbBXBF8p+vANgjPBxdmOOkYx6OETDHSck8HCMw/eM8zzk/NhFnqgFcBcbYJc75iXc+NMdJxzwcI2COk5J5OEbg5BynkVAMBoNhTjEB3GAwGOaUeQrgT8z6ADxijpOOeThGwBwnJfNwjMAJOc650cANBoPBcJh5ysANBoPBMIAJ4AaDwTCnzEUAZ4x9kDH2CmPsdcbYT8z4WN5ijD3PGHuGMXbJeW2JMfYFxthrzv/zA9//k85xv8IY+16Nx/UpxtgOY+yFgdd8Hxdj7DHn/b3OGPs/GPEQzxHH+c8ZY9edc/oMY+xDszxOxthZxtiXGGMvMcZeZIz9987rJ+p8jjnOE3M+GWMxxtjXGWPPOsf4L5zXT9q5HHWcJ+ZcDoVzfqL/AxAE8AaAiwAiAJ4FcP8Mj+ctACtHXvt5AD/h/PknAPzPzp/vd443CuAO530ENR3XdwJ4FMALKscF4OsA3g2AAfh9AH9pCsf5zwH8oyHfO5PjBLAB4FHnz2kArzrHcqLO55jjPDHn0/l5KefPYQBfA/CuE3guRx3niTmXw/6bhwz8nQBe55xf5py3AfwGgA/P+JiO8mEAv+z8+ZcB/NWB13+Dc97inL8J4HXY74cczvmXARyoHBdjbANAhnP+JLc/ib8y8G90HucoZnKcnPObnPOnnT9XALwE4DRO2Pkcc5yjmPpxcpuq89ew8x/HyTuXo45zFDO7hgaZhwB+GsC1gb9vYfyHVDccwOcZY08xxj7hvLbOOb8J2BcVgDXn9Vkfu9/jOu38+ejr0+BHGWPPORKLeJye+XEyxi4A+BbYGdmJPZ9HjhM4QeeTMRZkjD0DYAfAFzjnJ/JcjjhO4ASdy6PMQwAfph/Nsvbx2znnjwL4SwD+HmPsO8d870k7dsGo45rV8f47AHcCeDuAmwD+F+f1mR4nYywF4LcB/H3OeXnct444nlkd54k6n5zzHuf87QDOwM5SHxzz7TM7lyOO80Sdy6PMQwDfAjA41v0MgBszOhZwzm84/98B8DuwJZFt59EJzv93nG+f9bH7Pa4t589HX9cK53zbuXgsAP8BfZlpZsfJGAvDDoq/zjn/rPPyiTufw47zJJ5P57iKAP4EwAdxAs/lsOM8qedSMA8B/C8A3M0Yu4MxFgHwEQCfm8WBMMaSjLG0+DOA/wLAC87xfMz5to8B+F3nz58D8BHGWJQxdgeAu2FvcEwLX8flPMpWGGPvcnbOf3jg32hDXMgOfw32OZ3ZcTo/85MAXuKc/68DXzpR53PUcZ6k88kYW2WM5Zw/xwF8D4CXcfLO5dDjPEnncii6dkcp/wPwIdg77G8A+KkZHsdF2DvPzwJ4URwLgGUAXwTwmvP/pYF/81POcb8CjbvRAD4D+xGvAzsL+LjMcQF4HPaH9A0A/wZOt67m4/xVAM8DeA72hbExy+ME8B2wH3ufA/CM89+HTtr5HHOcJ+Z8AngYwDecY3kBwE/LXjOaz+Wo4zwx53LYf6aV3mAwGOaUeZBQDAaDwTAEE8ANBoNhTjEB3GAwGOYUE8ANBoNhTjEB3GAwGOYUE8ANBoNhTjEB3GAwGOaU/x+ORNIJ5szl5wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.plot(cp.asnumpy(y_smooth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([20.7, 17.9, 18.8, ..., 13.5, 15.7, 13. ])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the difference between the new smooth signal and the original."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "err=(abs(y_smooth-a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the data into a cudf series so that we can use the clx analytics library to calculate the rolling Z-score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We select rolling window as 5. This should be adjusted for other datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "rwindow=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "rz=rzscore(cudf.Series(cp.ascontiguousarray((err))),rwindow)[rwindow-1:] #Remove the NA rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4      -1.280745\n",
       "5      -1.101274\n",
       "6      -0.737006\n",
       "7       0.631617\n",
       "8       1.982907\n",
       "          ...   \n",
       "3645   -0.526520\n",
       "3646   -0.325860\n",
       "3647    0.507178\n",
       "3648   -1.702751\n",
       "3649    1.101080\n",
       "Length: 3646, dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f23196ae650>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0LElEQVR4nO3dd3wUdfoH8M+TEHqHUKQFBKSJgAFERUGKgAXFcngqnj89DtvpeZ6Hh2c5G+edDVERxcbZC4oSpUuTFpBeA0SJCRB6qCHJ8/tjZzezm6k7M7ubzPN+vfLKltn5fnd25pnvfNsQM0MIIUTFlxTvDAghhIgNCfhCCOETEvCFEMInJOALIYRPSMAXQgifqBTvDBhp2LAhp6WlxTsbQghRbqxatWo/M6dqvZfQAT8tLQ2ZmZnxzoYQQpQbRPSL3ntSpSOEED4hAV8IIXxCAr4QQviEBHwhhPAJCfhCCOETjgM+EbUgovlEtJmINhLR/RrLEBFNIKIsIlpHRD2cpiuEEMIeN7plFgH4KzOvJqJaAFYR0Wxm3qRaZiiAdspfbwBvKP+FEELEiOMSPjPnMfNq5XEBgM0AmkUsNhzABxywDEBdImrqNG0RX8t2HkDWvmPxzkZMbd9bgJISmVJclE+u1uETURqA7gCWR7zVDMBu1fMclD0pBNcxmogyiSgzPz8/qnws3JaPKYt3gZnx4qytyNpXENV6KoplOw/g659/03xv294CLN6+P6r1jpy8DANfXOAka46UlDB+2qGf99+9uRSfrvy1zOtTl/2CdTmHbae3MfcIBr20EG8s2KG7THEJ48lvN2LPkVO21+8lZsbCbflY9ctBDHl5IU4WFsc7S644WViMwycKDZd58tuNSH96dtRpLDc4fvS8t2QXLv3P/KjT9IprAZ+IagL4EsADzHw08m2Nj2gWk5h5MjOnM3N6aqrm6GBTo95Zgae+24TcI6cwYV4WRk5eFtV6KoKCU2cwcvIyPPDpGs33B7+0ELdMWY4XZm1FPG6Gc9s7K9D5sR803/ssczfOfXwmiksYh44XYsTrS5B7+GTo/XeW7MLv31qOuZv3an5++a6D+PuX68u8/s+vN+DqiUts5/W3Q4G0f/71sO4yS3ccwLtLsvHwl+tsr9+K2Zv2Ys6mvWVeO3rqjOHnvluXh1HvrMB1byzFlj0F2JQXeYjGxm+HT+LO9zNxorDIcLkzxSXIWJ9nuk9eNXExuv3LOJi/uyQb+48V4sVZW7Ep1/73/p3O8bP610P4bl2u5mee+HYTfjlwIvT8RGER+v/3Rzz93SakjZ2BvUfjUyBwJeATUQoCwf5DZv5KY5EcAC1Uz5sD0N5SLrpywiIAQFGUl+CfZ+7Gou3RXWWYOV1UjOOnAzv9V6tzkDZ2hq1SITPjno9W44cNe1BYVBL2Xt6Rk0gbOwPT1+bi3CdmWVrfq/OycPSU8UHohQXb8nFcp7T5+DcbUXC6CKeLivHl6hys/vUw3l60K/T+rv3HAQC5DkvTBafOIG3sDHyzJlCKKywqweSFO3CmOLBdt+w5iv3HTpuuZ/LCHcg5FDjIvar2+eMHmbjzg9LpRnIOncAfP8jE/R//HHptw29HylzBqE+UkXbkH8OynQds5WPv0VO2CwifZe7GRePnYc7mvZitOmltzjuK4ojtNXFeFu7+cHXYcloiqxR3HzyBvCOB77r/2GkcPF5a+p8wLwvDJiwyDLbfr8/DjHV5lr7PiNd/wr0f/Wy+IIBNuUexa/9xvL04sP++uyQ7VPOQe/gkhr2yCI9+vR6vzc+ytL5oudFLhwBMAbCZmV/UWWw6gFFKb50LABxhZmtb1YFDJwKlnsMnziD96Tm2P/+3L9bh1ikr3M4WAODa135C58dnAgC+WJUDIHDgAcDczXuRNnYG8gv0g0xRCWPGujyM+d8q3DolvAbt+jeWAgD+/LG1nTGICGEHSNDugycwZfEujU9E78jJM8hWAna01CHi3SW7cO9Hq6Naz+6DgQDxxo+Bqpr3ftqFZzO24L0l2QCAIS8vMq262nv0FJ7N2IKxX5W9oojWsdNF+N+yXwwD66kzgZPlLwdLS5NXvroYV09cgjW7D2PWxj2anysu4dBnB7ywwNZV8IjXl6D3s3PxeWZO6DVmxus/ZoUF08MnCjF9bWm57uEvyl71pI2dgaGvLMKr87aHvR4M2odMqmsi9X1+Pvo8Nw8AkP70HPR4qmzpP3i8RVq7+zDu+nA17olyP7Jj0oIdGPjiQhw8XogLx8/Dpryj+N+yX/GfmVs9TdeNEv5FAG4FcBkRrVH+hhHRGCIaoyyTAWAngCwAbwG424V0bbFSQrPj/Z+ykTZ2Bg5pBEgrtC6pg8f1ez9l6y6jZfmug2HPgwdLpIXbjK9Wvl2bix5PzcbPvx4Kez1YRbZtb4FrjbRXvboY/f77o6Vl1fFurU7d+5PfbsJ3FktmZo6dDgTC46pqh8MnjKtMIkuobnhy+kY8+vUG/LTDXuk76JrXlmD01FUAAidztXs+Wo0O/9SuStOz9+gpnCwsxmqlSmvG+jzsOXIKaWNn4JOVu/H8D1vDTrr3ffwz/vzxz/hVVbWhZ33OEVt5cUtxCYcKOfO37nO0rlNnim23jWidkLzkRi+dxcxMzNyVmbspfxnMPImZJynLMDPfw8xnM/O5zFzup8D8eEWgMTAvojrhpdnbkDZ2BoqKS7Q+pil4MB4+WYi0sTOwKMpG1CC90DPqHeOrlWBgiQycBUr98OCXFrrWSPvrQfMgoA5SpDxZ9cuhUDWbVsOQlmDVmVk6W/YUYPLC0gbZU2dKsFszn+bB/djpIhwxOUkAwJETZ/D2op2hUnxxCYceBwPRCY0gsjO/7Il36Y4DeOATa1d1elePR06cCZX8I/V+di5+/3bplcCCbfmYuiwbQOAKCwCOny79bLAaqbA49g3ERlfHn6z8NVRQe2HWVvR4ajYOOCwQnikuQdcnZ6GjTntUovDVSFu9HTka2oEAeGVu4NK00CDgR16ikxK6tu89Zrhc+Gfct0YpuU1ZvAvjv98SalizWlW73eQKgJnxypzt+M2gPlnzcxHPcw7Z+3znx2fqNq4B4SeWZzO2hB5PWrADfZ8v7WnxrXIinLPZvCS4ZvdhnPcv8/aTcV+vx9MzNmPZzoMoKi7B2f/IwHPfbymTr0hb94T3PDt68gxumbIcX69x1jR23r9mocM/f8DUpdma70c2WL82X7/HUlAc+gLgD+/qF252HzyJ7krJeqZS5aVVlWnHta8vCbWlfbVau8ooEfgq4B89qV3iOna6CLdOWa4bxLUEGxoz1ufZPpGszD5kvpCipIRDjYdeUwfiSQt2hE5eVg0yuQLYkX8ML83ZhjFKNYNbPl5etutlJKPqLLJ4+vx2rfv9DI4o+2RhcUmokDB1afh05lYaR/cfK3RcrZQ2dkbo8Yuzt0W1jk15Rw1L11acKCzCzI3GjbVmNlrsjRO8ctTbcj+qqnnyC04jbewMLNDYlzb8Vpreg5+t1UjHUnY856uAr2fWxj1YtH1/VDv5xPlZeC5js63PFJWEB3CjneEf09aj3bjvw17bmX8M50TUv6r70lstUW3MNa43DdZHurWzBuORW1dawa/ppIvh9r0FGPG6/S6a6r7fh08Uhqr4tDAzZm3cE0VADg9G6s8zAj29Js6z1qvD6kkt6NCJM6b7h55Vv4QXaKzsP3O37MM+pcF38EsLQydCryUpedM6Zjb8dgR/eHdl6Pn63w4DCPSxtyp4so7HVY4WXwX8yNb3/ILTGPXOitDOpVeSOnS8EF+tztF9/7fDZbt5HTtV5LhfOwP4ZOXuMq9/uTqnTPBYmX2wzHJmrpiw2Fo+HO6sWfsKwq5SrK4uGCdW/XIIT31XOlMHs/5vddPkZbb6Wv/7h61lu4Va+MLqvt9//WwtXtUJvOOmrUfrRzIweuoqTF64E/9b9gs265ygIpMNBso/TV2Fqct+wdiIvv3vLcl2XIVj5IoJi6Par8b8bxU+XF56laK1OfOOnMKKiM4G/5gW6OFkt8rOieCJsEQjk8EutpEYwMR51q5+E21Qtq8CfmRVyluLdmLhtnw8+e0mnU8EDHppAR78bC1aP5JhOa1ez85F3+fnlxkFeOpMMZZk2W+UDQb4dxbv0hz4E9xh521xdims5iTQl5QwduQfQ96Rkxj44sKwgG3X0xGf/ce09eiqM75g6c4DeHz6hqjTisZ+g/rfD1XVTbmHT+LRrzdg6CuLyiyn7ievVSL+ePmv+HyVuhskcNLGlVK0V2k3TFpq6QQaeQUxbprxbzD++y248c2lYa9pXQEx648h+NFhrxrAeLtE1qSqv+N/Z1mrDYjHYEYjvgr4ZvR+mv3HzBp0tD+Zc+gkHonolz1swiLdhq7FESeC21WXk2/8GChB/uu7TZrd9IIHi7ou0SsPf1G2jlLt6KkzeGH2Vgx4YQGWZAXyqi7NWY09wfrV7RoNwQUmPW+s0jzg41DhGrmf5B4+aTroKFas9IXfuteb6UveWrQTF46fh//M3FLmvYc+L90P9UZbW6UVl4t1grWdGB48h9nZpdLGzgg1JrvN1wE/8jf4Zk1uVGfkwmLWvUwvUI1e3ZR7FDvz9QcbRdZ9qmUfOKE7YAQoDfhJLsYq1jmRfZYZno/Ifv9dn5gVOqn9pnF5bnULx610ZDHdglNnbHW/nbpM997SZZj1vSeyHnge/GyN5XTd5MbI7R3K8aIuJM3fsq/M+Jc73nfW01trXy+z/0VxbOkdQ2Yy1nszLtWN6ZErlKU7D2DSgp3o2KSW7jI7IvpAL9yWj4Xb8rHo4f5lllXXDeqVlMjC6f+LVTmWAr6VddkxYe52HDDpsnbzW8sx76F+mu+pv7/VnBWcOoPTRc56Jlmpe16fc0SzumKCxYZQq9NWeMHOufCr1b+h0xW1o07r5reXa+7bZuz21rGy7x45eQa3v7cS57eqB6edk9PGzkDL+tUBaG/PyHr9AzpX+kbTISRYjY7/SvhDXl5o+P6J08VYuC0fby7cqbvMgBe0ux5qjcb8accB/LAhcLa++e3ISUQDjOY5sSp46ZjkcsC30nNp5/7juvOT2O3aCQD9/vMj0p+eozvHTtBHOt0xd+Qfxw2Tlmq+p3bVxMW2xwTEy+Y98ZnsLGj1r/pXn2aemmGvF5uR4BXVrv3HXal5MzopRV68BauQIrtlGk2HcLqopEzjtBVenSh8F/C37ClwdQCWmt7lm1mjsBvTFezIP4YPl//iapWOHb2fnWt5WbNBLmZXFE7Xn+hOFBaXmRAvMgBEW1UQD2ZTeqjFut0ieNKI3L5rdh8OayOI1nMZm3Hjm0uxKS8xpmj3XcAHgA7//AHHThd5M1w1ThZsy8e4aRtcLeF7Vcoo7wHZK+rNHeyi6Ba3q/q8sGb3YfzxA+O6ePU2cvMbRZ5A31fms3JqszIi+p9f2+s55tXP5cuAD+iPup1kcHMLM3qDW/KOnNK8EYcXnrEwCGxl9kFLc+LEowzpdAZNPUdOnrE1kjrR2R1IVVzirE0kcoI+L5jdyAQAnlOmvtD69j2fsT8jbnA9kYWbSm5dKidYJb5vG20vHD8PdaunlHk906CnjBmjy2ytG3HEy1PfbUrYWxNanUHTrpkb92Lmxr3IHn+FJ+uPNbtVOuo5gqKh117iJrLQ9ehLg3lqopnSQW9qhc8NOkjYkVjh3sclfMB8ylu7CqLohmalVCP8IdEG6bjJyhTJdr+/G9UepSV8b7Z9ov2kvg74btPrhWPE7PZs8WRUqrN7Z54texKj0aoiSbBYYsjKFMl2piE4cLzQdrWWpmCjrfM1lQsS8H3m1JliVxq7vL4zj1+oG1PtNKwmWsnRTKLmV68O3y3R9qbyaoZct+5p+w4R7SMizaZoIupHREdUd8R6zI104yHaaWMTxWiXpyYWzqirEipylY7W/EFOuVKlo6zkd2+aj9uIJa0JGd3gVqPtewAmAvjAYJlFzHylS+nFzYQoBhIlkoXb8nFe8zpxSbsc9AyMmcEvLUCfNg3inY2YKUq0aSMVp4sCVU1FJYyX57h/bCfaOdyVEj4zLwTgfb8t4Yp47YNnivVT/iyz7DTQFdm2vcfwfsSNTjbbHJxT3gsfieDUGW9vLhT1ec6jM0Us6/D7ENFaIvqeiDrrLUREo4kok4gy8/Otj9AT1q2L0w2j7zSY4OrhL9aVuW2f30ROF2wkwQqOcaE1h31F4dU3i1XAXw2gFTOfB+BVAF/rLcjMk5k5nZnTU1NTY5Q9EQtm89ZEc+ep8i7aEbAVub7fqr1Hnd1KMRb0ZtE1U67n0mHmo8x8THmcASCFiBrGIm1RfphNllYRSeC2btF2/1zxezVXUkwCPhE1IaUoQ0S9lHSNJ/wWwgcWbbd/9zO/uu/jn+OdhZjxqhzgSi8dIvoYQD8ADYkoB8DjAFIAgJknAbgewF1EVATgJICRLEUbIaJ2/LQ7V0PP/+BsyoVYkojhnCsBn5lvMnl/IgLdNoWIqw+WZsc7C654e5H+/RrseP3H6CcLFN4p13X4QiSKx77ZGO8suKIi91DR46dKgfLeS0cIIRxx4x655YVXJzcJ+EKUQ9kWZp8UIpIEfCGE8AkJ+EIIkWCk0VYIIXyiXA+8EkIIEX8S8IUQIsFIlY4QQviE9MMXQgifkH74QgjhE1LCF0IIv5A6fCGEEE5IwBdCiASzc/9xT9YrAV8IIXxCAr4QQviEKwGfiN4hon1EtEHnfSKiCUSURUTriKiHG+kKIYSwzq0S/nsAhhi8PxRAO+VvNIA3XEpXCCGERa4EfGZeCOCgwSLDAXzAAcsA1CWipm6kLYQQwppY1eE3A7Bb9TxHea0MIhpNRJlElJmfnx+TzAkhhB/EKuCTxmuaQwuYeTIzpzNzempqqsfZEkII/4hVwM8B0EL1vDmA3BilLYQQArEL+NMBjFJ661wA4Agz58UobSGEEAAqubESIvoYQD8ADYkoB8DjAFIAgJknAcgAMAxAFoATAG53I10hhBDWuRLwmfkmk/cZwD1upCWEECI6MtJWCCF8QgK+EEL4hAR8IYTwCQn4QgjhExLwhRDCJyTgCyGET0jAF0IIn5CAL4QQPiEBXwghfEICvhBC+IQEfCGE8AkJ+EII4RMS8IUQwick4AshhE9IwBdCCJ9wJeAT0RAi2kpEWUQ0VuP9fkR0hIjWKH+PuZGuEEII6xzfAIWIkgG8BmAQAveuXUlE05l5U8Sii5j5SqfpCSGEiI4bJfxeALKYeSczFwL4BMBwF9YrhBDCRW4E/GYAdque5yivRepDRGuJ6Hsi6uxCukIIIWxw4562pPEaRzxfDaAVMx8jomEAvgbQTnNlRKMBjAaAli1bupA9IYQQgDsl/BwALVTPmwPIVS/AzEeZ+ZjyOANAChE11FoZM09m5nRmTk9NTXUhe0IIIQB3Av5KAO2IqDURVQYwEsB09QJE1ISISHncS0n3gAtpCyGEsMhxlQ4zFxHRvQBmAkgG8A4zbySiMcr7kwBcD+AuIioCcBLASGaOrPYRQgjhITfq8IPVNBkRr01SPZ4IYKIbaQkhhIiOjLQVQgifkIAvhBA+IQFfCCF8QgK+EEL4hAR8IYTwCQn4QgjhExLwhRDCJyTgCyGET0jAF0IIn5CAL4QQPiEBXwghfEICvhBC+IQEfCGE8AkJ+EII4RMS8IUQwick4AshhE+4EvCJaAgRbSWiLCIaq/E+EdEE5f11RNTDjXSFEEJY5zjgE1EygNcADAXQCcBNRNQpYrGhANopf6MBvOE0XSGEEPa4UcLvBSCLmXcycyGATwAMj1hmOIAPOGAZgLpE1NSFtIUQQljkRsBvBmC36nmO8prdZQAARDSaiDKJKDM/P9+F7AkhhADcCfik8RpHsUzgRebJzJzOzOmpqamOMyeEECLAjYCfA6CF6nlzALlRLCOEEMJDbgT8lQDaEVFrIqoMYCSA6RHLTAcwSumtcwGAI8yc50LaQgghLKrkdAXMXERE9wKYCSAZwDvMvJGIxijvTwKQAWAYgCwAJwDc7jRdIYQQ9jgO+ADAzBkIBHX1a5NUjxnAPW6kJYQQIjoy0lYIIXxCAr4QQviEBHwhhPAJCfhCCOETEvCFEMInJOALIYRPSMAXQgifkIAvhBA+IQFfCCF8QgK+EEL4hAR8IYTwCQn4QgjhExLwhRDCJyTgCyGET0jAF0IIn5CAL4QQPuHoBihEVB/ApwDSAGQDuJGZD2kslw2gAEAxgCJmTneSrhBCCPuclvDHApjLzO0AzFWe6+nPzN0k2AshRHw4DfjDAbyvPH4fwDUO1yeEEMIjTgN+Y2bOAwDlfyOd5RjALCJaRUSjjVZIRKOJKJOIMvPz8x1mTwghRJBpHT4RzQHQROOtcTbSuYiZc4moEYDZRLSFmRdqLcjMkwFMBoD09HS2kYYQQggDpgGfmQfqvUdEe4moKTPnEVFTAPt01pGr/N9HRNMA9AKgGfCFEEJ4w2mVznQAtymPbwPwTeQCRFSDiGoFHwMYDGCDw3SFEELY5DTgjwcwiIi2AxikPAcRnUVEGcoyjQEsJqK1AFYAmMHMPzhMVwghhE2O+uEz8wEAAzRezwUwTHm8E8B5TtIRQgjhnIy0FUIIn5CAL4QQPiEBXwghfEICvhBC+IQEfCGE8AkJ+EII4RMS8IUQwick4AshhE9IwBdCCJ+QgC+EED4hAV8IIXxCAr4QQviEBHwhhPAJCfhCCOETEvCFEMInJOALIYRPOAr4RHQDEW0kohIiSjdYbggRbSWiLCIa6yRNIYQQ0XFawt8AYAQMbkhORMkAXgMwFEAnADcRUSeH6QohhLDJ6S0ONwMAERkt1gtAlnKrQxDRJwCGA9jkJG0hhBD2xKIOvxmA3arnOcprmohoNBFlElFmfn6+55kTwis/PtQv3lkQIoxpwCeiOUS0QeNvuMU0tIr/rLcwM09m5nRmTk9NTbWYhL8N6NAo3lkQGurVqBzvLAgRxrRKh5kHOkwjB0AL1fPmAHIdrlOotG1cE3O37It3NkQk3WKNEPERiyqdlQDaEVFrIqoMYCSA6TFIVwghhIrTbpnXElEOgD4AZhDRTOX1s4goAwCYuQjAvQBmAtgM4DNm3ugs28bqy6W0EEKU4SjgM/M0Zm7OzFWYuTEzX668nsvMw1TLZTBze2Y+m5mfcZppM0nGvYYqHNJoJln0cP845EQIkcgq5EhbZn9Vnmqd3yrCOa9G5WTbn6lcqULu0iJGerSsG+8seKpCHh3+CvfaTMZGlAvtm9Sy/ZmP7uztQU6ARrWqeLJeIWKpQgb8kgQt4bdrVNOT9WqF9vIf7qP7Dl798vde1tbW8lefdxZqVXU0rlEAqFc9JabpJWbkcE/FDPglsf3ZOp9V29Jyoy9p43FOSlWAAn6Yc5vVsbScV+f6i9s21Hy9tk5Qn3BTdyQlufcjNK9XzbV1WVGrSmKcrC7r0DjeWTBVno61ChnwY32WvvWCVpaWG9ypCepUi02JRashtzz6dPQF+MvA9vj2vostLd+yfnWPcxSuZ1r9mKQT644Ifx3cPqbp6enSzFphKp4StEJBU8UM+KofYMylZ8cvIxFqVEl2rTQw58FLDd8vT6UOPQygd5sGuH9gO8ufaVKnqml9u9194uM/XqDbJnJ5lyahx1avQoR1I3u2DHtePYqGfFGqggb80ojPMSjvWw2ulZKT0L1FXVfSrKSqLtDspeNKKuVTSrLxbt2+sb22lD5nN9B8/b83nIcb00sHkcdiX/ObahEB3uv9ujyV1qNRIQN+jKvwcXnnJuYLKSb+vocraaqDfHorjWoFl46MccM62v7MdT2a49lrz3WcdiKdtLS6+ga7gH59z0Xep1+BTyZ/u/wcy8s2qJl4vaWiuZqOVxtJhQz4w85tWvokBsdJZCnESA0Pfui+7Rpi/IjwAOtWHX40DY+Na1fx5NL7kaEdXFmP29VdyTGoP2tWN7aNtkYm3WJeaPm/i1pbXl8VG2Mn2qTWsLxsNGJ1Wh3UOT6N0RUy4L9w43lxS3vdE4Pjkm5kY7CLHURs8yr+/SmB2mNizW5jf7UU7+q6h3RparrM+a3qYcJN3S2tr3vLeobvqxvizapcrul2lqU0dUUkQOReQUPt39d1dX2dVlTIgK/mxhn7/Fb18MrIbrrvq0vTtauWHpjPX9/VVkknWloNivEceJXoU1tEc/WjtT0jq3m8rP+1u26nV5JO9x+9j/dpU9oeEqzWOL9VacB32smii8sN59ufHupqQeO923ti57PDTNuZvFLxA74LR2GT2lXL/EDbnh4aKnXp7dw3prfAY1eV77s52jnsg71U4h3u7fzmeoPh7PS4Cf7+idTg5/Sc69X0JI9e2TFUhbP0HwOw49lhYe+3qG9cdWWWK60T1dy/XqrZXfdsC9VDVk58wU3V1sLASiJydXyGXRU+4HuFCJj5wCV4/eYeMTtbEwFT7+iFq887C83rhe/AkQdCPHapUMOixWhjNmDNq6sU9Wr1RtB+Y6MhNnhF49UI7/LaFVF7BDiF9tVKSYTkGAW/tIbe1v0DgYJhoqvwAd/LUleTOlXDG4g9tuu5K9C3XSom3NTd9EAxq1bp264h3ru9p5vZC21rq4dwp6bxHVRTNSUJgzppN55ZKYUFv29SUvhzt31514W2q2jshtEFf+tn8xPRqVY52bCY7sU2JACv3tQdU25LN102gS7SPFHhA36sdGtRFy//rlvUn7+gTfQjNgkaB4rJEX9td93bCoevxkbkYHsFfPPLc4vpxqKdxGhenGCbgJ0SvvqeDZWTk1DZ4CqxXaOaeOKqznhwkPXRr3Yvjlo1CC8BO67DB9BbY59u3bBG6EpQKwmnAVevKqpOtRQM6GjcM+btUeYnhPKuwgf84M//p0udzWNjdix/fc9FuEYJoj3TjBt5tXRo4qy0G8tJpvS+W3ATWW0UNdumWgHhpd+V7YHlpJ3Eal4bGvT/Dl4M2AlW36mmimjZoLpmcAwiItSpnoI/D7A+4thr3S1MI9yoVngVx7f3Br5z6ZWgVsQvuxXTlUbdP1yYpplOtFNiN1ZVwaTWqoKBnRqXSd5KjVOC91EI4/SOVzcQ0UYiKiEi3dMjEWUT0XoiWkNEmU7SjFZKkvZXbdUgurlXjH7jz8dciOHd9EvQk24533I6F+qM8iyzXMTkXmY7odXCqNZq9L5bsHSllfZ9GnXlHZvan/742u7NNV+/Mb05+p8TftP7L8b0wduj0vHDA33x+s3uDHjTQlHU4Z/lQb/6P6u2sdNxGOqSslZAnXb3RfjP9fa6FlZKDuTpzVvPxwVt6iMl2Voeq6QE0m/XuKZmCT7aZoCJv++Bvw8x7nJp5UrHqDCQaJyOAtoAYASANy0s25+Z9ztMz7Kpd/RCtZRkZKzfA0A/AFo6Rgno2ty97l5DupQdmXt3/7Px3k/ZZV5/4+bz8dXPOYbdzbS+m5WukW7XVx44XhjIj8Z7l3VohFfnZaFG5WR8+qc+aNuoJpbuOGC4PjtB6/nrS0v+we/VtG41pCuTm3VoUhsDOzbCnM37wg5iO6UzovD9JVg1ESrhe1QBHE08c7PU2alpbazZfdi19Q3o2Fi3esVoE+rtD0lhv6f1Lsr1a1TGiB7N8O8fthikauz567qidrVKGPO/1WXWfVA5HhKJ01scbmbmrW5lxk1926WGDnbAea+VFh7Pwhh5+RtUp3oKbr+ote1ZGdXfN3v8FXgjooTrRWzKLzgdSFtjYwcbmYuZ0aVZHVS1MjBI50czu2ex3cZjqxY81F+zodvrXjrxqDJQB0m97+VFLyqjTag3vYSTcR9O95Ube7bQbFC3sz69qiovxKoOnwHMIqJVRDTaaEEiGk1EmUSUmZ+f70LCzg/C5gk0rD1aWr1OrOyUkQd15OWrVpUTEeHS9uHVK6VBsfS1aH8breqhaNiNzy0bVEe/cxqV+bzXAd/MXwaWbczV+22j6Y1TrDM5lVG3WnfPBcYrU6fFzLi7X/hAKa1Pd1R6iEU2ILs1Z5FeFaJWXp64urMraVphWqVDRHMAaM0ONo6Zv7GYzkXMnEtEjQDMJqItzLxQa0FmngxgMgCkp6e7dwRFsQcmUeCHKw83YYhk9nUb1Khse9f+YkwftIxo89Aa8k8E1KtRGYM7NcasTXsBqEphFhK9t39bTJyfZTN35rwq+Qe/WklJlJ83Xb/FRnALnzG7OtKiF/A7GnSrjfbcp66j1+p9pbXeJrWrouDUsdDzh4d0QHEJ482FO3XT+f7+vmHrC1YXWc33ynEDDdsgerex1vYWa6YlfGYeyMxdNP6sBnswc67yfx+AaQB6RZ9le5we5EO6NPXkxtgNa9o/8PRoHdyRl7mRS/SLaOC0Ij2tfpmqp+GquUuC26lpncAyD6saxIJt5upSsN7BdWHbBpp5tku33Ub1vpOSaHAwTzCI3trH2o1wnNCaRM1eO4T9LxyLK5e3RqVj5gOXhJ20gr2vzBplp97R23HXXLPNom7cXfv4YKTWqoK61QO/eyKNsDbjeZUOEdUgolrBxwAGI9DYGxM1lfo1vRkt4zXt7JKxl+GmXi3NF/SI1QPfbLFU1c1G/ti3Nabclo5rlF48bRvVRBslKGpVe+geKCb9+dWfszO1rtt6KJN+1ahSCdnjr7A8D0zkmAs7e+Dcvxrf+ObT0RfYWJs1RUoJ/+Eh54T1BDJi97wyqFNjnNOkluY+8dDgczCyZwtc10O7h1aTOlXLHMdWt6nV5e5SVRMZTWT38u+64elrulhca+w57ZZ5LRHlAOgDYAYRzVReP4uIMpTFGgNYTERrAawAMIOZf3CSrh33XtYWfx/SATecr72zxEuVSsm2poX1gttVG41rV8WAjo3DTibBUm8TpdRv1vg8onszW/35o7llZHCtzO5NI23Hvf3D+9MTrJ+AtRq7g+0qDWtW8aS7Z/Ae0UM6N8GDg6M/wUZ7NVWvRmWMv66rpYZ+u1cwoa7EUeUsmGbp42u6N8MtFm95Gg9Oe+lMY+bmzFyFmRsz8+XK67nMPEx5vJOZz1P+OjPzM25k3KqqKcm4q9/ZqKTUKbjVmBTP2Sj1qKdmtpK9i9o2xM29ja8yzFdTuoRW6ez2i1oje/wVqF01Bd/f3xdT/mA8ncOLqtHK0W5is6s2q+v97r6LPZlyIDL9JCJHk5WN7NkCr4zsZjnQjOzZwnwhlWIlb0bTeUT2AouW2VYw+221tqPR7106OtxeHb5dF0eMk4mXCj/SNiSK4OF1UHd7RsLaVVNCO3eZOnyN75KSnIRnXLgzVZDZ9+nYtHaoig2Ivmuo1Z8lsvSulZ56XZFTF3RpVqfMlANuCo62dbKbdWtRF0lJhOHdmlmaiIwAjL+uK6beYb0ZLdgYbbT+oTGcU0pLPK7UrJr4++7xzkKIfwK+iYw/98Xax+Jz8xIvRO7+Tg6Hd2/viU88qBsOGtixUdjzWDWCqQNtchLZmrrAbkP+wr/1L/Oa0/sGrHlskO7v4mZZpaHSTqP1nWM9yjSafcPoZBA5/5MXu96VXR3elMVF8bmxYjxo/JKz/nIJ7nh/JYBA426wYbdG5WQcLyw2XF0ilCciR34C+gdE1LeGI0L/cxqZLwf7B4veFYHR5FqBz0Wfpt567HhqeGfdG5vrUXdnDbUhhL4nRXU1GewloubFyfKtUedj0bb9moMDFz3cP1TlEy66I8SrefjNJGANrSf8E/Aj1KxSCe0bh8/lUrlSEsYN64hL2qfi8pc1hwkACIxcdYPTXfu7+y7G3M37wl4b0KER5m7ZVyaAtEmtiXVPDEbXJ2bZSuNKG5fqbh2rhpNrRUjRqGYwnZTN4en61j5pjj4fZNZl+LM/9UHG+jxX0orG52P64PjpIjSqVRXX6XR60Ov9FpyjqmPT2ticd9S1POn9tg2Ubs52G/HL9O4x2Hk+urM3zuiMSSgv/BPwLR7jf7ykDc4URzmCxqHb+rTC7E17kXvklKXlO59VB53PCp9j57WbeyC/4LRm32X17RetqmcyUCdslKPttYfWYrpuPSN0uuppfT5epUc9ZtNJ92pdH71aRz9tthbL7R9k3qNKy4gezfDgoPahG/R8dGdvdH9qdmCdFg7CaH+i0Ze0QWqtKqGum6P6tMJkZeCVpUZbC3mLnKAwKNghRO/kp6aXl7rVU3BL71ahAYfvunyviiDf1uEHD36jHczrq7zItJ8c3gU/PTIg9Hza3RfaXmfVlOSo5v0ZqjGhm112A6pJN3xTf7gwLapBcVpVRlXj2EWWyLvRv5GCd2drVd+4ii+awJs9/gq8eGO3sLux1atRGe0bm9/6L6h/B+PqQ71eOinJSbgxvUVoCpHm9aqjeb2yXVTrVU8J61OvHoQXrd6t6+MvA9tj/IjoO0CseWwwHlKNKbFajWqXf0r4CShYCr/qPO1Gne7KwB6vrXtiMKqlJKPduO9tf9ZJoLJ7wxTL6zV5Xz2MPiU5CX8f0gEDOnpzgBkJnXhUt/1za52RLm2fikm3nB8K+C0bVA9V8UUz6torVu4L68TPER0z3OiHn5REuH+gO/cq6H9OKuZvdT6HmB4J+BqCvSfcvpyO9OCgc3CmhPHYlfG90Xk0VT1BbVKjP0CDc5FUixhQY3alELz7VF2Tm77oHcSRJ5i7+lkbIeu2aE94r9/cw3z+HWWJutVT0CutPsZd0bFMlUPtqilY9ehA1NS4m1dFaMS0c5XidRfs4L5qdGczAJhyW0+dRnB3+C7gEwIlwMgfOLKL3g8P9EWLet5OiVynegqedbEfvJreDlwtJRl/GeTenZPqVEvBHRe3xpTFu2xXA/Q7pxH+fFlb/N/FrTGqT6syk3TpfYfrejTH6aIS3JiuPYDonMa1kF9w2pM5kKJRSaf/eumIYnus3kd5ym3p6HRWbTStoz/6tkE5unmHFyJ32eA+PLBjY6zNOexaOs9cey7Ob1XPtBCZlERI8rAy2XcBv0blSrilTytcrVSjvHpTd7w2Pys04VeQ01sOJqrNTw2xtNzE33fHidPGXVODot09k5MoNFRffe8Cs/NGUhIZjip9/ZYe2JBzpEy3xfS0+pi/Nb/Mb+2lV0Z2w7mRN68Jzq6pKuJ7cYib3cNV6PeUeujy9q7GgDrVAve1iDffBHx1YVE98133lvXw9m3etIiXV7WqVIpqsIhrE9E5nOG0dtUUzR4Vd116NoZ2aWJY4nWb5u0gle+nFWzevb0nbn93ZdTpuVUbkOpiyd/NGgr3azvCf4SOTWth694C1KhcMUNjYlzzijDVLXTv8pJ6Th4rQqMUXYv3xgOvopWURGiTWjOB6qe9+Z5O1/nWqHTN23A6ZTVPvaLoDhptupEn3edGdMXnY/p4foe7eKmYp7EE1Lqh9ZGuK8YN1L3pRCxYacBa9siA0Bz3I3u1xEfLf8UVXd2ZTyXYk0Tr1nEVSbD/du2qKaFqppoOv7OTvWbs0A5o16hm3KuCPrijFwpOFWm+Z+f7tWpQHb8dPmnYllPaLVPpQFA5OarxB1r+e8N5oQFoWmY+cAlqV4vtPl6xjyiV6kpPkLHDjO9S74Vpd19oaxIupwd9LDRR1YOfnVoTG/9lrW3AiovOboiHBrdP6GlmHVHOp12b18GjV3TEtd2boXrlSuhzdgPXgk00BXyr8/nbdW6zOti+75jl/bpqSrK1ex6beOPm87Ei+6Du/aIBoGX96mhYszLGDevoOL1I15tMyX5Ok1qG73sh8SOLSyolJ7k2JYJdsepP79ScBy/BPuVG5PGUlES49zL3ehJFCs76OLiz+9UWz1zbBe0aWTuQiQh39m0Teh6s7//q7gtRX2OeHCsaKROdjb4kPl1NtTw74lzc2qeVK3P133dZW4x6ZwW+GHMhrnvjJ8OuuXWqp2BQJ+Orlaopych8dJDjfJUXjgI+Ef0HwFUACgHsAHA7Mx/WWG4IgFcAJAN4m5nHO0k3Xr4Y0yfeWbDNzh2h2jaqhbYWg1V5lpKchOX/GIB6GkF1VJ9W6OPgfqQ393Z+VdLDQQEhePetRFI1Jdm1Qk/fdqnY9Vzg+yXa9ywPyMn8IkQ0GMA8Zi4ion8DADP/PWKZZADbAAwCkANgJYCbmHmT2frT09M5MzMz6vwJkUjSxs4AAGx48vJyUW0nyiciWsXM6VrvOb3j1SxmDrauLAOgVWnVC0CWcuerQgCfABjuJF0hyjMJ9iJe3Nzz/g/ApxqvNwOwW/U8B0BvF9MVolx44qpO6OnxdB1CGDEN+EQ0B4BW69Y4Zv5GWWYcgCIAH2qtQuM13XokIhoNYDQAtGxpfL9VIcqTPyTASEvhb6YBn5kHGr1PRLcBuBLAANZuEMgBoJ70pDmAXIP0JgOYDATq8M3yJ4QQwhpHdfhK75u/A7iamU/oLLYSQDsiak1ElQGMBDDdSbpCCCHsczq1wkQAtQDMJqI1RDQJAIjoLCLKAAClUfdeADMBbAbwGTNvdJiuEEIImxw12jJzW53XcwEMUz3PAJDhJC0hhBDOyORpQgjhExLwhRDCJyTgCyGET0jAF0IIn3A0l47XiCgfwC9RfrwhgP0uZscr5SGf5SGPgOTTbeUhn+Uhj0Bs89mKmVO13kjogO8EEWXqTSCUSMpDPstDHgHJp9vKQz7LQx6BxMmnVOkIIYRPSMAXQgifqMgBf3K8M2BRechnecgjIPl0W3nIZ3nII5Ag+aywdfhCCCHCVeQSvhBCCBUJ+EII4RMVLuAT0RAi2kpEWUQ0NgHyk01E65XZRDOV1+oT0Wwi2q78r6da/hEl71uJ6HIP8/UOEe0jog2q12zni4jOV75fFhFNICKtG964mccniOg3ZXuuIaJhqvdinkdl/S2IaD4RbSaijUR0v/J6om1PvXwmzDYloqpEtIKI1ip5fFJ5PdG2pV4+E2ZbamLmCvMHIBnADgBtAFQGsBZApzjnKRtAw4jXngcwVnk8FsC/lcedlDxXAdBa+S7JHuXrEgA9AGxwki8AKwD0QeDOZt8DGOpxHp8A8JDGsnHJo7L+pgB6KI9rAdim5CfRtqdePhNmmyrrq6k8TgGwHMAFCbgt9fKZMNtS66+ilfDLyw3ThwN4X3n8PoBrVK9/wsynmXkXgCwEvpPrmHkhgINO8kVETQHUZualHNhzP1B9xqs86olLHpV85jHzauVxAQL3fWiGxNueevnUE/N8csAx5WmK8sdIvG2pl089cds/1SpawNe6YbrRDh0LDGAWEa2iwP16AaAxM+cBgYMQQCPl9Xjn326+mimPI1/32r1EtE6p8gle2idEHokoDUB3BEp8Cbs9I/IJJNA2JaJkIloDYB+A2cyckNtSJ59AAm3LSBUt4Nu6YXqMXMTMPQAMBXAPEV1isGwi5h/Qz1c88vsGgLMBdAOQB+AF5fW455GIagL4EsADzHzUaFGdPMUkrxr5TKhtyszFzNwNgftf9yKiLgaLx21b6uQzobZlpIoW8G3dMD0WOHD3LzDzPgDTEKii2atcykH5v09ZPN75t5uvHOVx5OueYea9yoFWAuAtlFZ5xTWPRJSCQBD9kJm/Ul5OuO2plc9E3abMfBjAjwCGIAG3pVY+E3VbBlW0gJ9QN0wnohpEVCv4GMBgABuUPN2mLHYbgG+Ux9MBjCSiKkTUGkA7BBp0YsVWvpRL6wIiukDpWTBK9RlPBA96xbUIbM+45lFZ7xQAm5n5RdVbCbU99fKZSNuUiFKJqK7yuBqAgQC2IPG2pWY+E2lbavKqNThefwjcS3cbAq3g4+KclzYItMyvBbAxmB8ADQDMBbBd+V9f9ZlxSt63wsPWegAfI3DJeQaBUsYd0eQLQDoCO/UOBG5qTx7ncSqA9QDWIXAQNY1nHpX1X4zAZfg6AGuUv2EJuD318pkw2xRAVwA/K3nZAOCxaI8Zj7elXj4TZltq/cnUCkII4RMVrUpHCCGEDgn4QgjhExLwhRDCJyTgCyGET0jAF0IIn5CAL4QQPiEBXwghfOL/AT9t2DLjKs43AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.plot(cp.asnumpy(cp.asarray(rz)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this plot, we can see that the rolling Z-scores mostly vary between -2 and 2. This view may help us select a threshold."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the data back into a CuPy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/envs/rapids/lib/python3.7/site-packages/cudf/io/dlpack.py:74: UserWarning: WARNING: cuDF to_dlpack() produces column-major (Fortran order) output. If the output tensor needs to be row major, transpose the output of this function.\n",
      "  return libdlpack.to_dlpack(gdf_cols)\n"
     ]
    }
   ],
   "source": [
    "rzn=cp.fromDlpack(rz.to_dlpack())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.28074544, -1.10127388, -0.73700588, ...,  0.50717802,\n",
       "       -1.70275054,  1.10107965])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rzn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Select a rolling z score threshold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This threshold should be adjusted according to the nature of each dataset. For this dataset, we will select the rolling z-score threshold as 1.99, and will mark any value that exceeds this as an anomaly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "rz_threshold=1.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(rzn[abs(rzn)>rz_threshold])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indices of the outliers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1079],\n",
       "       [1806]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cp.argwhere(abs(rzn)>1.99)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have identified two outliers when rolling z-score threshold is selected 1.99 and returned the indices of the outliers. This work can be expanded as a monitoring tool for time-series data in the cyber domain and to highlight the anomalies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Rasheed, Peng, Alhajj, Rokne, Jon: Fourier Transform Based Spatial Outlier Mining, 2009\n",
    "##### BRIGHAN, 1974, The Fast Fourier Transform\n",
    "##### MULGREW, GRANT, THOMPSON, 1997, Digital Signal Processing\n",
    "##### OPPENHEIM, 1997, Signals & Systems"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}