{ "metadata": { "kernelspec": { "codemirror_mode": { "name": "ipython", "version": 3 }, "display_name": "IPython (Python 3)", "language": "python", "name": "python3" }, "name": "", "signature": "sha256:56ee8b7ea5e38567dafb3c8b98bb3d6b2f93787914bb6ee7f7a978dec06fd473" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "import os\n", "from IPython.core.display import HTML\n", "\n", "with open('creative_commons.txt', 'r') as f:\n", " html = f.read()\n", " \n", "name = '2015-03-16-outlier_detection'\n", "\n", "html = '''\n", "\n", "

This post was written as an IPython notebook.\n", " It is available for download\n", " or as a static html.

\n", "

\n", "%s''' % (name, name, html)\n", "\n", "%matplotlib inline\n", "from matplotlib import style\n", "style.use('ggplot')\n", "\n", "from datetime import datetime\n", "\n", "title = \"3 ways to remove outliers from your data\"\n", "hour = datetime.utcnow().strftime('%H:%M')\n", "comments=\"true\"\n", "\n", "date = '-'.join(name.split('-')[:3])\n", "slug = '-'.join(name.split('-')[3:])\n", "\n", "metadata = dict(title=title,\n", " date=date,\n", " hour=hour,\n", " comments=comments,\n", " slug=slug,\n", " name=name)\n", "\n", "markdown = \"\"\"Title: {title}\n", "date: {date} {hour}\n", "comments: {comments}\n", "slug: {slug}\n", "\n", "{{% notebook {name}.ipynb cells[1:] %}}\n", "\"\"\".format(**metadata)\n", "\n", "content = os.path.abspath(os.path.join(os.getcwd(), os.pardir, os.pardir, '{}.md'.format(name)))\n", "with open('{}'.format(content), 'w') as f:\n", " f.writelines(markdown)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "According to Google Analytics, my post\n", "[\"Dealing with spiky data\"](https://ocefpaf.github.io/python4oceanographers/blog/2013/05/20/spikes),\n", "is by far the most visited on the blog. I think that the reasons are: it is\n", "one of the oldest posts, and it is a real problem that people have to deal everyday.\n", "\n", "Recently I found an amazing series of post writing by Bugra on how to perform\n", "outlier detection using\n", "[FFT, median filtering](http://bugra.github.io/work/notes/2014-03-31/outlier-detection-in-time-series-signals-fft-median-filtering/),\n", "[Gaussian processes](http://bugra.github.io/work/notes/2014-05-11/robust-regression-and-outlier-detection-via-gaussian-processes/),\n", "and [MCMC](http://bugra.github.io/work/notes/2014-04-26/outlier-detection-markov-chain-monte-carlo-via-pymc/)\n", "\n", "I will test out the low hanging fruit (FFT and median filtering) using the same\n", "data from my [original post](https://ocefpaf.github.io/python4oceanographers/blog/2013/05/20/spikes)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from datetime import datetime\n", "from pandas import read_table\n", "\n", "fname = './data/spikey_v.dat'\n", "cols = ['j', 'u', 'v', 'temp', 'sal', 'y', 'mn', 'd', 'h', 'mi']\n", "\n", "df = read_table(fname , delim_whitespace=True, names=cols)\n", "\n", "df.index = [datetime(*x) for x in zip(df['y'], df['mn'], df['d'], df['h'], df['mi'])]\n", "df = df.drop(['y', 'mn', 'd', 'h', 'mi'], axis=1)\n", "\n", "df.head()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
juvtempsal
1994-11-22 04:00:00 1056.1667-0.1 0.0 23.9 34.6
1994-11-22 05:00:00 1056.2083-0.2 0.7 23.9 34.6
1994-11-22 06:00:00 1056.2500-0.1 2.0 23.9 34.6
1994-11-22 07:00:00 1056.2917 0.0 3.1 23.9 34.6
1994-11-22 08:00:00 1056.3333-0.1 2.7 23.9 34.6
\n", "
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 2, "text": [ " j u v temp sal\n", "1994-11-22 04:00:00 1056.1667 -0.1 0.0 23.9 34.6\n", "1994-11-22 05:00:00 1056.2083 -0.2 0.7 23.9 34.6\n", "1994-11-22 06:00:00 1056.2500 -0.1 2.0 23.9 34.6\n", "1994-11-22 07:00:00 1056.2917 0.0 3.1 23.9 34.6\n", "1994-11-22 08:00:00 1056.3333 -0.1 2.7 23.9 34.6" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Starting with Bugra's `get_median_filtered()` we have:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "\n", "\n", "def get_median_filtered(signal, threshold=3):\n", " signal = signal.copy()\n", " difference = np.abs(signal - np.median(signal))\n", " median_difference = np.median(difference)\n", " if median_difference == 0:\n", " s = 0\n", " else:\n", " s = difference / float(median_difference)\n", " mask = s > threshold\n", " signal[mask] = np.median(signal)\n", " return signal" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "import matplotlib.pyplot as plt\n", "\n", "figsize = (7, 2.75)\n", "kw = dict(marker='o', linestyle='none', color='r', alpha=0.3)\n", "\n", "df['u_medf'] = get_median_filtered(df['u'].values, threshold=3)\n", "\n", "outlier_idx = np.where(df['u_medf'].values != df['u'].values)[0]\n", "\n", "fig, ax = plt.subplots(figsize=figsize)\n", "df['u'].plot()\n", "df['u'][outlier_idx].plot(**kw)\n", "_ = ax.set_ylim(-50, 50)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAawAAADJCAYAAAB/hgbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8XOWZ6PHfVNVRtYptWXKTXze54YZtXIIpBkwJkKWF\nQEghZG8gyd1kk81mQ3Y3m8ZdQhoJKSSEQMAQagwGbFNccbdl+5Vly02WZfU+mnbuH2dGlm11S5oZ\n6/l+Pnw0c+bMmcfDnPOct1sMw0AIIYSIdNZwByCEEEL0hCQsIYQQUUESlhBCiKggCUsIIURUkIQl\nhBAiKkjCEkIIERXs4fpgn89v1NQ0h+vjhRBCRKCMDJels9fCVsKy223h+mghhBBRSKoEhRBCRAVJ\nWEIIIaKCJCwhhBBRQRKWEEKIqCAJSwghRFSQhCWEECIqSMISQggRFSRhCSGEiAqSsIQQQkSFi5qa\nSSmVCWwHrgQCwNPBv/uAL2utZTljIYQQ/aLPJSyllAP4DdAEWID/B3xba704+PymfolQCCGE4OKq\nBH8C/BooCz6fpbX+IPh4NbD8YgITQggh2utTwlJK3QdUaK3XBDdZgv+FNALJFxeaEEIIcVZf27Du\nBwyl1HJgBvAnIKPd6y6gtruDZGS4+vjxQgghhpo+JSyt9ZLQY6XUOuBB4CdKqSVa6/eBFcB73R2n\noqKhLx8vhBDiEtVVQaa/FnA0gK8DTymlnMB+YFU/HVsIIYTAYhhh63luSAlLCCFEexG54rAQQgjR\nG5KwhBBCRAVJWEIIIaKCJCwhhBBRQRKWEEKIqCAJSwghRFSQhCWEECIqSMISQggRFSRhCSGEiAr9\nNTXTkGDfthVrdRUAgbR0fLPnhjkiIYQYOqSE1UP2bVuxVlWZsyYaYK2qwvH+Oiz1deEOTQghhgRJ\nWD0UKlkdPF7Dn946SKvHh8Xtxr5zR5gjE0KIoUESVi+dqmyitrGV2kZPuEMRQoghRRJWDwXS0gHw\n+wMA+AIBjNhYfDNnhTMsIYQYMiRh9ZBv9lwzQfnN5Vi89hi8S5ZhJCWHOTIhhBgaJGH1gm/mLFps\nDjz2GKonTA13OEIIMaRIt/ZeMJKS0eNnUWivYUFcYrjDEUKIIUVKWL3k9ZltWB6vP8yRCCHE0CIJ\nq5c8wYTlDXa+EEIIMTgkYfVSqIQV+iuEEGJwSMLqJUlYQggRHpKwesnjM9uuJGEJIcTgkoTVS1LC\nEkKI8OhTt3allAP4A5AHxAD/BRwAngYCwD7gy1pro3/CjBxtvQR90ktQCCEGU19LWHcDFVrrxcC1\nwC+Bx4BvB7dZgJv6J8TIYRhGWy9Bn5SwhBBiUPU1Yb0IfLfdMbzALK31B8Ftq4HlFxlbxAlNywRS\nJSiEEIOtT1WCWusmAKWUCzN5fQf4abtdGoFLbpI9b7tqQI8kLCGEGFR9nppJKTUKeBn4pdb6OaXU\nj9u97AJquztGRoarrx8fFtX17rbHVps16uIXQoho1tdOF1nAGuAhrfW64OadSqklWuv3gRXAe90d\np6KioS8fHzYVtS1tjxubPVEXvxBCRLquCgJ9LWF9G7PK77tKqVBb1sPAE0opJ7AfWNXHY0es9tWA\nXplLUAghBlVf27AexkxQ51t6UdFEuPZtWDKXoBBCDC4ZONwL7XsGSi9BIYQYXJKweqF9laD0EhRC\niMElCasXpIQlhBDhIwmrFyRhCSFE+EjC6oX2qwxLwhJCiMElCasX2vcMlIQlhBCDSxJWL3i9Z5NU\nwDDwByRpCSHEYJGE1QuhJUUsluBzryQsIYQYLJKweiFUDRgfY463lsHDQggxeCRh9UIoYSXEOQBZ\nE0sIIQaTJKxeCA0WToi1n/NcCCHEwJOE1QsXVAlKwhJCiEEjCasXQpPfxsc6gs8lYQkhxGCRhNUL\n3vOqBNvP3i6EEGJg9XnF4aEo1GYlJSwRSQzDoKK2hUMn6zAMmDMpkxiHrcN9G1u8vLC2mBNnGgFI\ndcWQPyqZ/JwU8rJcOOxyDysilySsXriwhCUJS4RHS6uPw6V17CupZuehCipq3W2vvbCumBsuz+Oq\nOaOwhAYNAvuPVvPUG/upa/TgtFuxWCwcK29gV3ElAHabldkqg1sWjyUjJW7Q/01CdEcSVi94fQGs\nFguxTvPuVXoJisHS7PayZX85ew5XUVbdTEVtC4ZhvhbrtHGZymDCqBQamj2s3V7K82uLqWvycNvS\ncVgsFnYequBXf98HwG1Lx3Ht3FysVgvV9W6KS+s4dKKOwqPVbN5fzjZ9hqtmj+LGhWOIcdpobPGy\ndvtJNu0vZ0y2i08uHsuwHiS0Vq+fhiYPXn+A7LT4c5JnJAkEDF79qIQFBdlkpcaHOxzRBUlYveDx\n+XE4rNiD1SZSwhIDrbHFyxsbj7J+Z2nbDZIr3sGEnBTG5ySjRqWgclPPqcr7xKwcfvTXnazecpyq\nejcZKXG8teU4NpuFr94+HZWb2rZvWlIsc5NimTspi4BhsPVAOS+tP8LqLcfZeuAMKYlOjp5uwB8w\nsFktlFc3s02f4eYrxrJiXu4FSajJ7aWwpJoP95Sxv6SaYE4lOy2eK6YNZ/LoNDJT4zhaVs+JM400\ntHixWGDuxCxyMhMH/PvsSHFpHa9vPEqr188dV+aHJQbRM5KwesHrC+C0W9suDjLThRgoHq+f97af\n5I1Nx2hp9ZGeFMPSmSO5fEo2aUmxXb43JTGGb9w5kx//dQdbD5wBwGm38sht5yar81ktFuZPzmZm\nfgavbShhzdYT1DS0kpuVyLzJWSyePoLdhyt5YW0xq9Yf5nBpHXMmZuL1BSg53cChk7WUVjS1HW/M\n8CSy0+Lx+vzsKq7ixfWHgcMdfvYbG48xdkQSi6ePYM7ETOJiBu/SVNfkAcybAxHZJGH1gtcXwGG3\n4rSbVYJer/QSFH3T6vWzv6Sa/cdqqGvy4Pb4GJYcR2qik8o6N/tKqqlpaCUh1s4dnxjPslk5veoQ\nkeqK4fsPzOVkRRMNzV5GpMf3qBoPIMZh4/al41kxLw+HzUqM82wHjvmTs5mcl8aTr+5j56FKdh6q\nbHvN6bAyKS+V/Jxk5k7KYsSwhLbXGlu87C6u5NDJOipqW8jLdjFmeBLJCU7qmzx8tLeMvUeqOHKq\nnlXrD/PwbdMYNzK5x//ei1EvCStqSMLqBY8vQFyMXUpYok/cHh/bdQU7iiooLKnusg3UYbeyYl4u\n112eR0KwV2pvOew2xgxP6mu4JMZ1/LlJCU6+fscMdhRV0uT2YrVYGJWZyKjMROy2jpNqYpyDhQXD\nWVgwvMPXZ0/MpLrezfu7TvHmpmP89PldfOXWAiaNTutz/D3V0GwmrCZJWBFPElYveH0BkhOsOGzS\nhiV6xusLcPB4DdsOnmHrwTO0esxS+fD0eGZNyGD6uGFkpMbhtFuprHNT0+AmPTmOrNS4Ti/+kcBm\ntTJnYma/HjMtKZZbFo8lL9vFk6/u42er9vCDL8zvtgr0YtU3m4lKSliRr18TllLKCvwKmAa0Ap/T\nWndcaR2FvD4/DrsVh8O8kEgvQdGVzYWnefadIprcPgDSk2K4dm4ucydlMjw94YL9Q6WUoW7WhAzu\nuVrx9OqDvLbhKPetmDign9cgVYJRo79LWDcDTq31AqXUPOCx4LaoFwgY+PyG2elCSliiC16fnz/+\n4yCb95cT47Bx1exRzJowjPycFKzWyOzaHWkWFmTz1pbjfLSnjBXzcwe0u3l9sEqw2e0jEDDk/1EE\n6+86h4XAWwBa6y3A7H4+ftiE2qvs7XsJytRM4jyBgMFvX9vP5v3ljBuRxKOfncOdy/NRualyIewF\nm9XKLYvHEjAMXv2wpNP9Glu8VNW5MUKD0s4TMIxOXwsJdbowgOZWX59jFgOvv0tYSUB9u+d+pZRV\na31BUeTNDSUEvD7sNiuNLV4aWrw0NnsJBAwS4x2MHJbA9PHDImaqmFBpymm3ne0lGKUlrJqGVjYV\nnm4bPCr6h2EYPLNGs72ogom5KXz1U9Nx2DueIkl07zKVQW5mIlv2l3Pd/LwLxmmVVTXxo2d3UN/s\nJdUVw5yJmdy0aAxxMXbqGlt5d/tJ1u4oJT8nmX/+ZEGnbYKhNiwwE2BnnU1E+PV3wqoHXO2ed5is\nAJ58eU+3B3PFO7lqbi43Lh5LenJ4p4qx1LaYMSXEkJ1l9ryy2mxkZLi6eltEWv3xCVatP0zBhExm\nqf5tOB+q/AGDX764i/d3nWLsyGQe/eKCtjknRd/df+NUHv3dZt7ccpzvfHZe2/by6mb+94Xd1Dd7\nmZGfwZFTdaz5+ATbdAUjMhLYf6SKgAE2q4U9h6v4yzuH+Prdl11wg+b1+WlpV6pyxDii8pweKvo7\nYW0AVgIvKqXmA51mpX+9dw6l5fX4/AFccQ4S4x244pxYLNDQYo6W37C3jJfXF/Pah4dZVDCcGxeN\nISUxpp9D7pnymmYA/H4/9XXm44amVioqGsISz8U4cdosBJecqGFUmswZd7G8Pj9PvXGAbQfPkJfl\n4iu3FtDU4Kapwd39m0WXctPjGJ+TzJbC02zefZJxI5KpaWjlh89up7LOze3LxrFiXh5eX4DVm4/x\nxqaj1NS7GTcymcunZDFnUhZPvLSHD3aVkhTv4JOLx55z/Or6c/8fnSirIz1BbjTCqasbhv5OWH8H\nrlJKbQg+v7+zHRdOH0FFReeBTRmdxi1XjGVT4Wn+sfkY63edYmPhaW5cOIZr5+ViHeR5ybzeUJVg\n+zasKK0SrG81/za0hjmS6NfQ7OEXL+/l0Mk6JuQk85XbphMfK6NF+ovFYuHWxWP50V938pe3i7hh\nQR5//7CEilo3NywYzYp5eYA5bu3GRWNYNmskhmGOFQt5+LZp/McftvL21uMsmzmSVNfZm96GYHVg\nrNOG2+OXsVgRrl/PLK21AXypv47nsFtZPH0ECwuy+XBPGa98WMKq9Yc5dKKWz6+cPKhVLp52bVg2\nqwWLJXoHDocSlSSsi1Ne3cz/vribMzUtzJ2UyQPXT5I2qwGgclOZrTLYpiv4ZXAC3+Wzc7jlijEX\n7OuKd16wLSHWwY0Lx/D06oO8ueko91yt2l4L9RAcnh5PSVmDdG2PcJHRo6EbNquVpTNG8p8PzGXy\n6FR2H67iJ8/vGtReeqHPsgeXZXDYrW2lrmhiGAbVkrAuWtGJWv7rz9s4U9PC9Zfn8YUbp0iyGkAP\n3jSVb90zi2vn5XLL4rHceWV+r2Z/XzA1m4yUWD7YfYqqurPVgKEegqFxcZKwIltUJKwQV7yTr31q\nBpdPyebY6QaeX1s8aJ99tpegNfjXFpUlrIYWL75g3JKw+iYQMPj5S3twe/zcv2Iity4ZN+hV1EON\n1WohPyeFTy0bz8oFo3u9VIndZuXGhWPw+Q3e2HS0bXv7EhbI9EyRLuoq261WC/deqzh+poF1O0pR\no1KYOylrwD83lLBC7VcOuzUqx2GF2q8AaoZwp4CTZxrZXlTBoZO1NLt9ZKfFM2ZEEvMnZ51TrdTq\n9fPh7lOs3VHKtHHp3HFlPq1eP01uH9PGpXPF9BFh/FeI3rh8SjZvbjoWHIycR2ZKHA1NZoIaISWs\nqBB1CQvM2aQfunkqjz79MX99p4iCsekDvhyB57wSlsNmpTUaE1a7UlWT24fH68fZyXLqlxqvz0/h\n0RrW7yxlz+Gqtu12m5WjpxvYvL+cF9YWMzE3hRHDEqlpcLP3SDWtwVn57Tbzrj5UsnZGyBhB0TNW\nq4WbFo3hN68V8vqGEh64fvLZEtYwSVjRICoTFph1ztfNy+OVj0pYveX4Bd1V+5snmJxC7RQOh5Um\nd/T9uEOlKqfDiscboKaxNepWWfX6/Ow/WsPB4zXUN3nx+PxkJMeRmRZHcryT+Fg7FosFnz9AY4uX\nMzUtHDpZR9GJ2rbkM2FUCssvy2FiXirxMXYq693sOlTJh3tOUXi0hsKjNQBkpsYxf3IWa3eUtpWy\nfeeVtkX0mDMpkzc2HmXjvtNcf/notoSVmhhDrNNGY4vMdBHJojZhAVwzN5f1u0pZs/U4S2eMGNBZ\nnc+/SDls1qjs1h7qcDEmOwl9opbahshPWF5fgPKaZo6XN7DrUOU5pZ7eyEqLZ2b+MOZMzLxg2Y3M\nlDiunjOKq+eMotntpay6mTinneHp5tLumwpPt5Wyz68eFtHDajFLWb96ZR+vbzhKQ5OXGIeNGKeN\nxDhHVN6EDiVRnbBinDZuWTyWP/7jIKveP8wXVk4ZsM+6oErQbsXjC2AYRq8bgMOpOtiGNXaEmbCq\nI7jjhdfn542Nx1i95XhbRxEwSz3m0hzppCfH4rDbOFPTzJmaFhpbvDS7fRiYsxwkxjlIdcUwbkQS\nyT0cdB4f62DciHMXD3TYbbS0eoJxBROWbWhUpV5qZqkMRgxLYOuBcpwOK654c3hMQpyDU5VN3bxb\nhFNUJyyAhVOHs25HKZsLy1kwJZupY9MH5HNCCSu0tEgocfn8Bg579CSsUJVgqIRRG6EJq7rezU+e\n20l5TQuprhgKxqYzIj2eKcG/598kJCc4yc9JGbB4HDZrW9tV6K+UsKKT1WLhmrmj+OM/DtLS6m/r\n0p4Y58DrC9Dq9RMzRNp1o03UJyyr1cJ9Kyby/ae38ee3Nf/5wLxzlvTuq2a3j4/2nKLoZB1lVU1t\nsziHlhZxtJsAN5ouXDUNrSTFO8gILpceqSWs5949RHlNC8tmjeS2JeMGvFNNdxwOa1u1sFQJRr/5\nk7N5+YMj1DV6SAr2Cg1NetvU4pWEFaEuiTMuN8vFivm5VNa5eWaNJhDoejmB7lTUtvDfz2zj+bXF\n7CiqoL7JQ2Ozl/gYO1lpZnuPPQqXGDEMg5qGVlKTYklNMqvHInEs1p7DlWwvqmBCTjL3XDUh7MkK\nzBsVf8DAHwhIwroEOOxWrpo9CqCtSjAxOHOO9BSMXOG/EvSTlQtGs6+kmo37TuPzB/jcDZP7tMT4\niTONPPb8TuqbvSyfncM1c3JJT47FMAwMg7bZnkNVgs++e4gYh5XT1c34/QbTxqUzZ2ImIzMib+XY\nJrcPjy9AmisGV5wDu80ScQnL4/XzlzVFWC0W7rlGRUz7YPv5IyVhXRqWzRzJsdMNzJ+SDUBCnHk5\nlMHDkeuSSVhOh41/uWMmP1u1m60HzlBR6+Yz1ypys3q+VEBNQyuPv2guWXD3VRO48rKcttcsFnP+\nwJDR2S427jvNtoNnALOBH+Do6QZe23CU+ZOzuGXx2Laqt0gQmpk61RWDxWIhJTEm4gYPv7juMJV1\nbq6dm0tOBCX9cxKWtGFdEuJi7Hzp5qltz0NVgo1u6doeqS6ZhAUQH2vna/80gz+9dZDNheV8/+lt\njB2ZRFZqHIGAQUurn4KxaSwsGH7BYFm3x8cTq/ZQ09DK7UvHnZOsOrJ89igWFgynscVcdDI9ORav\nL8DeI1Ws3nyczfvL2VFUwT9dmc/SGSMioqQQKk2FZqtOc8VwqLQOfyCAzRr+i++ew5W8t+MkI4Yl\ncHMHE5uGk7NdwvIEu9Q7+lCCF5GrLWFJCStiXVIJC8xZML6wcgoLpmazat1hDpfWUXyyru31XcWV\nvPpRCbcvG8+CqdlYLBaa3V4ef3EPx8obuGLacK6dl9ujz4qLsZ/TvmK3WZk7KYvZEzPZUljOX98t\n4pm3NfuOVPHFG6cM6owSgYDB4VN1FJZUc6qqmaq6FuqD09CkuczxaimuGAwDdh2qJMUVQ1OLj9rG\nVsqqmqhrNLtwx8bYGTciifxRKWQkx/Zr4tXHayiraqa+2cOZmhZ2F1dit1n4wsrJETf7hpSwLn2J\nwbas594tYu32k2Snx5OZEofHG6Cu2UN5dTOVdS0EAmaNSkZqHFmpcSTFO0lKcJKdFk9etiuialXC\nyeP1c/R0Ax6vn+z0eNKSYi96zs1LLmGFTB2TztQx6Xh9Aarq3W13w+t3lfLu9pP8/s0D7DpUSU5m\nItv1GU5WNDF/chaf7od2E6vFwuVTs5mYl8pTrxey81AlP1u1h6/cNm3Aex+VVTXxwe5TbCosb5uJ\nGsxkarNZzDFJOeYYo9CJFVqyoTPrd5YCkJLopGBsOisXjGbYRZyUXp+fZ98p4oPdZedst9us3HN1\n76pxB0tozJW0YV268nNSWDpzJMdO11NW1UzpeWOynA4rmSlx2Gxmj9FTlU0cO33uAq4WC3zppqnM\nnjh0V/JuaPbw7DtFbNcV+Nt1gIuPsTM+J5lFBcP7/P1YDOPietRdBCNcq/VW1Lbw1Bv7zyl5LZ05\nknuuntDvs257fQGefHUfOw9Vokal8PDt04h1Dsx9wpqtx/nb2mIMzOqNWROGMX38MPKyXG3tVu01\ntnj5+EA5dU0eWr1+EuMcJMU7GZ6eQFqSuX9tYyvFJ+s4dLKWopN11Dd5sNssXDM3l5sWjel1x5b6\nJg+Pv7ibo6cbyM1K5Jq5ubiC3eyHJcdGRNVkR15YW8xbW4/z75+ZzcFjNby4/jCP3D6NaeOGhTs0\nMQAMw6C20UNlXQsxDhuueCfJic5zrg+BgEFtYyuNLV5qG1s5VdnMqxtKMAIG37x71gWzqVzqmtxe\ndh2qZNX6w9Q1eRg5LIEpY9KIi7FTVtXE0bIGztS2AHDDgjxuvmJsh9fbjAxXpxfhIZmwwPyxHThW\ng9ViVo2FBg8OBJ8/wG9fK2SbrmB8TjJfvX16v3bVNgyDVz8q4bUNR0lJdHLHlfnMzM/o9xJAwDDY\nsr+cl98/TFV9K/k5yXzp5qmk9HAGiep6Nz95fhfl1c0snJrNp69REVf115mXPzjCGxuP8s27ZqKP\n1/LKRyX8yx0zmDQ6LdyhiQiyu7iSJ17aQ1K8k0c/O/eclY8vRaWVTWzad5qik7WUnKrHHzCwWS18\ncvFYrpmb29aruv3+P39pD2dqWrjyshzuvmrCBcfsKmHZvve97/X7P6KHvtfc7Ol+rwFisVjITI0j\nIyWuw1VK+5PVamGWyqC8upm9R6rZdvAMPr9Bdnp8v1ywP9xTxt/WFpOREss375rFhFEpbb0W+5PF\nYmFUZiKLp4+gvKaFfUeq2bTvNMPTE8hO63o+wpMVjfzk+Z1U1rpZMS+Xu6+e0KdhB+FyuLSOA8dq\nmDc5izO15mS6S2eMHND5K0X0yU6Lx+mwsqOokprGVmar3lV9BQIGFXVuTp5p5NjpBuqbPLgSnBF3\nrgQMgzUfn+DXr+xDn6ilpqGVvGwXy2bl8OlrFDPyMzpsWkmKdzJvcha7DlVSWFLNgoLsC1aOT0iI\nebSzz71k27Aijc1q5Qsrp5AU72T9rlJeWFfM6i3HePCmqUzKS+3zcavq3Dz/3iHiYmz8y50zGZY8\n8A2+sU47X7ppCu+MTGbV+mJ+tmoPn5g1kn/6RH6HpbrCkmp+9cpeWlr93LpkLNdfPnrAY+xvMg5L\n9NQ1c3LZrivYsr+c+ZOzmD6+62rjitoWtukz7DxUydGyhnPmzQSzg8cV04ZzzzUqIhYKbWzx8tvX\nCtlXUk1SvIM7l09g2rieL/Hkindyw4I8fvfGAdZ8fIK7ll9YyuqMJKxBZLVauOuqCdy4aAzrd5by\n6kcl/PT5ndxztWLZzJG9Pp5hGPxx9QFz5dvrJg5KsgqxWCxcPWcUk/JS+e1rhazdUcrR0w08dPPU\nc0odH+w+xZ/f0litFr544xTmTR74xTYHQkcJyy4JS3QgNF3co3/8mD+9dZC7lk9g8uhU9pVUU1px\ntiOHARw5Vcf+4FI2FgvkZroYMSye9ORYEuOcVNe72VVcyfpdp7DZrNy1PH/Qh8gYhkF9s5fTVU2c\nqmpm9eZjVNa5KRibzgPXT+pTtefcSVm89P4RPtxdxk2LxpBwXimrM5KwwiAxzsENC0YzMTeVX7y8\nh+feLWLy6NReLfMRCBj86a2D7D9aw7Rx6SwqGD6AEXduVGYi3/nM7Laxb//++63cvnQco4e7+HB3\nGet2lpIY5+D/3FowoJPTDrRQL1MpYYmeyMlI5NYl43hxXTG/eqXrXrgTcpJZUDCcGfnD2uY1bG/l\nwtH88NkdvLf9JCmJzgGtoTAMg9LKJg6dqOVIWT2nq5opq2pum0sVwALcvGgMNywc3ecSn91mTo31\nwrpi1u0o5YYFo3v2vj59mugX43OSuedqxa9e2cff3ivmK7dN69H7Wlp9/PEfB9imK8jLdvHA9ZPC\nOjA5xmHj8zdMZsKoFF5cV8yf39Ztr2WlxfPI7dMifs2t7oRm6W8/DktWHBZduXZeLtPHp/PWluOU\nVTczdXQa+TnJ53RESHXFkNnNuZEQ6+Brn5rBD57ZxsvvH2FkRiIzuqlm7AuvL8DTqw+yqfB02zab\n1Wzrn5iXyvD0eLLT4hkzPIkRwy6+k9qSGSN4fWMJ724/yTVzR7VNKN4VSVhhdpnKQI1KYVdxJfuO\nVHW7PMqOogqefaeImoZWJoxK4eHbpkXE5LAWi4WlM0YyfdwwXt9QQqs3wMz8YRSMS78kZr6WcVii\nL4anJ3D/dZMu+jiprhj++ZPT+J+/bOep1wv5zr2z+7Vnc7PbyxMv7aXoRC2js10snTmScSOTyUqN\nG7AOH3ExdpbOGMnqLcfZVFjO4ukjun1Pr690Sqlk4C+AC3ACX9Nab1ZKzQceB3zAGq3193t77KHI\nYrFw5/J8Hn36Y377+n6+cus0xuckX7Bfs9vHs+8UsanwNHabhRsXjub6y0dH3EUz1RXDvddODHcY\n/a6tDcsfwOOTqZnE4MvLdnHfion89vX9/ObVQr7zmdn9lkxeWFdM0YlaZk/M5HPXTxq04SbLZ49i\nzccneGvLcRZNG95tFWNfbs2/CryjtX5CKTUBeA64DHgSuEVrXaKUelMpNUNrvasPxx9ycrNc3HuN\n4pm3i/jxcztZMn0EyYlOEuMdxMfYOVxaz9aD5dQ1ehid7eJzN0zulyK56LlQwvJ4/W3rYkmnCzHY\n5k/JZv/RGj7aW8Y7H59gxfy8iz7myTONfLinjJHDEvjijZMHdfB+qiuG+VOy2LD3NLuLK5mZn9Hl\n/n1JWP8LhNakcAAtSikX4NRalwS3vw0sByRh9dCS4JieX7+yj/d2nLzg9bgYGysXjGblwtERNyZj\nKHC2K2GCQNlfAAANY0lEQVR5fQHsNktEdDEWQ8+nPjGePYcreeWjEmapjItuH35hXTGGYR43HDPN\nXDM3lw17T/PC2mLyc1LoKmV1mbCUUg8Aj5y3+T6t9XalVDbwDPAwkAzUt9unARjbh9iHtIKx6fz0\noQWU17TQ2OKlsdlLY4uX4enxTMxLlUQVRud3a4+0qlgxdCTGObjrqgk8+Wohq9Yd5sufLOjzsfYd\nqWJfSTVTRqcydUx4Zm3JyUhkxfxcVm8+zq/+vpefPLyk0327TFha698Dvz9/u1KqALMq8Ota6w+V\nUkmYbVohSUBtd4FmZETeJKeRIG9UuCMQ53MHx3LaHXYCQIzDLr9fETbXDUvk7Y9PsKu4Enusg1RX\n72dc8QcMXv7TNiwW+OKt08nMDN/chw/eOoPaJi+b9pZ1uV9fOl1MBl4Ebtda7wXQWtcrpTxKqbFA\nCXA18L3ujhXOuQSF6I3GenPSzvoGN+5WHzar/H5FeF0+OYuSU/W8/n4xK+b1vi3rg92nOFpWz6Jp\nw0l0WMP+e7736gk4uqm46Eu9xg8wewc+oZRap5T6e3D7g8CzwBZgh9b64z4cW4iI5HCc263d3oMx\nI0IMpPlTsrHbrHy4u4zeTmLu9vj4+wdHcDqs3HJFZLTexDhs3Lei6yEAvS5haa1v7mT7FuDy3h5P\niGhw/kwXrnhpwxLhlRjn4DKVwZb95RSX1vV4Jhm3x8fPX9pLXZOHGxeObluBPBrIWSdED5ztdOHH\n65dOFyIyXDHNnJJt3Y7SHu1f09DKY8/v4sCxGmbmD4u6iajDP0WCEFHAbrNgATzSS1BEkIl5qYzK\nTGTL/nJuWDC60/GZAcNg3Y5SXnr/MG6PnwVTs7n/uokRu2BqZ6IrWiHCxGKx4LBbaWkNznIhCUtE\nAKvFws2LxmAAr20o6XAfj9fPb14t5Nl3irBaLHzmWsVnr58UdckKpIQlRI857FaaW72ATHwrIseM\n/GHkZbvYeuAM181vIDfr7HCLsqom/vCPAxwurSc/J5mHbikgOYpXQZaEJUQP2e1WmtzmMgtSwhKR\nwmKxcMsVY3j8xT389zPbWVQwnGHJsZyqbGJTYTkBw2D+5Czuv25S1P9uJWEJ0UMOm5U6j6ftsRCR\nomBsOvetmMgbG4+ybufZDhjZafHctnQcM/OHhXUJov4iCUuIHmo/g3W036mKS4vFYmHx9BEsLDAn\nxw0EDFzxTvKyE6OyraozkrCE6KH2pSqZqV1EIpvVSkE3a+pFMznrhOih9qUqKWEJMfjkrBOih85J\nWNKGJcSgk7NOiB6SEpYQ4SVnnRA9dG7CkslvhRhskrCE6CEpYQkRXnLWCdFD7dutpA1LiMEnZ50Q\nPSQlLCHCS846IXrI2a7dSuYSFGLwyVknRA/ZpYQlRFjJWSdED0mVoBDhJWedED0kUzMJEV5y1gnR\nQ06H9BIUIpzkrBOih87p1i4lLCEGnZx1QvSQtGEJEV5y1gnRQzI1kxDh1ef1sJRSE4HNQKbW2qOU\nmg88DviANVrr7/dTjEJEhPZJStqwhBh8fTrrlFJJwGOAu93mXwN3aq0XAfOUUjP6IT4hIoZUCQoR\nXr0+65RSFuA3wLeAluC2JCBGa10S3O1tYHl/BSlEJGifpOw2SxgjEWJo6rJKUCn1APDIeZuPAc9r\nrfcopQAsQBJQ326fBmBsP8YpRNiFqgEddisWiyQsIQZblwlLa/174PfttymlDgEPBJNZNmZpaiXg\nardbElDb3YdnZLi620WIiOEOmH+dDpv8doUIg153utBa54ceK6VKgKuDnS48SqmxQAlwNfC97o5V\nUdHQ248XImwa61sAszpQfrtCDIyubgb73EswyGj3+EHgWcAGvK21/vgijy1ERAm1YUkPQSHC46IS\nltZ6bLvHW4DLLzoiISJUwp7tzDnyMamuGOzb7Phmzw13SEIMKXKrKEQP2LdtxVlbA4aB3WrBWlWF\n4/11WOrrwh2aEEOGJCwhesBaXYXNasFiOVslaHG7se/cEebIhBg6LrYNS4ghw2KxsKhgBCmJMeEO\nRYghSRKWED0QSEvHWlXFrAkZbduM2Fh8M2eFMSohhhapEhSiB3yz52LExrY9N2Jj8S5ZhpGUHMao\nhBhaJGEJ0UO+mbMwYmOlZCVEmFgMw+h+r4FhyOBLIYQQ7WVkuDqd90xKWEIIIaKCJCwhhBBRQRKW\nEEKIqCAJSwghRFQIZ6cLIYQQosekhCWEECIqSMISQggRFSRhCSGEiAqSsIQQQkSFIT/5rVLKAfwB\nyANigP8CDgO/De5yCPic1tofngg711HsWuvXg6/dBfyz1npBGEPsVCff+0ngDaAouNuvtdYvhCfC\nrnUS/xbgKSAFsAD3aq2PhivGznQS+11AdnCXMcBGrfVd4Ymwc53Efgj4HeYK6EWY52vE9SbrJPbj\nwJOAD/Pf8aDW2hO2IDuhlLJh/rYnYH7PDwKtwNNAANgHfHmgv3cpYcHdQIXWejFwLfBLzB/Sv2qt\nFwX3WRmu4Lpxfuy/AFBKzQQ+G87AeqCj730W8JjWelnwv4hMVkEdxf8j4Bmt9RLgu8DUMMbXlQt+\nN1rrO7XWy4BbgBrgq+EMsAsdfe//gXmzdgVmIrg+jPF1paPYnwK+Goy9FHgojPF15QYgELwmfgf4\nAfAY8O3gv8cC3DTQQUjCghcxLy5gfh9e4Fat9UdKKSfmXWdtuILrxgWxK6XSgP8GHsH8EUWqjr73\ny4DrlVLvK6V+p5RKDFt03eso/oXAKKXUO5gXp7Vhiq0758fua/fa94EntNblgx5Vz3T0vbcA6Uop\nC+ACIq6EEtRR7Dla683BbRuBJeEIrDta61eBLwafjsa8qblMa/1BcNtqYPlAxzHkE5bWuklr3aiU\ncmH+oP5Na20opXIxi7npwJ6wBtmJDmL/LmaVw9eAxrAG142OvndgK/B/gyWUI5h3zhGpg/i/g3ki\nV2utr8Ks6vlmGEPsVCffPUqpTOATmNU8EamT2H8B/AzYD2QC74cxxE518ps5opRaHNxlJZAQtgC7\nobX2K6Wexvyun+XcG+JGYMDX2hnyCQtAKTUK8274z1rr5wG01se11hOA3wD/L5zxdaV97Jh14OOB\nXwPPAZOVUlERe/B7/7vWemfw5VeAmWELrgfOi/85oAp4Lfjy68DscMXWnY5+88BtwLOR2P7TXgex\n/wW4Qms9CXgGs6oqInXwm/ks8C2l1LtAOVAZzvi6o7W+D1CYbYax7V5yMQg1UUM+YSmlsoA1wDe0\n1k8Ht72mlBof3KURiLgOF3Bh7Frrj7XWU4NtEXcA+7XWXwtvlB3r6HsH3lJKzQk+vhLYFo7YeqKT\n+D/ibPvJEswSesTpJHYwv/PVYQmqhzqJPR4IrVVUhtnpJeJ0EvsNwN1a6+WYtTlvhym8LimlPq2U\n+lbwaQvmNXGbUipUhbkC+KDDN/ejIT81k1LqZ8DtgG63+d+AH2PWhTdh9jqKuDr9TmJfobV2K6VG\nA3+N4F6CHcX+r5h3x17MC88XtNYRWbXZQfwGcB/mnWcC5t3mXVrrurAE2IVOYr8O8wZhgda6Plyx\ndaeT381PgX8H3Jg91z6vtT4ehvC61Ensj2G2G7ZiVok/EoklXKVUHGZVcTbgAP4HOIjZacSJWR37\n+YGOfcgnLCGEENFhyFcJCiGEiA6SsIQQQkSFIZmwlFJLlVK1Sqmcdtt+qJT6TDjjEkII0bkhmbCC\nWoE/tnsujXlCCBHBhmrCMjDHQlQppb7c/gWl1NeVUluVUhuVUj8MbvtYKZUXfHybUurxwQ9ZCCGG\ntqGasEIjtB8CvqqUGhd87sLsdnp5sDt4vlLqeuD3wL3Bfe7j7MS4QgghBslQTVgAaK2rMefc+xPm\ndxELbG43M/uHwBTgr8BtSqnhQJLWen844hVCiKFsSCcsAK31G5gD+e7DHHg4TyllC06kudjcRdcD\n24HHMefqE0IIMciGasIyOLeTxSOY043UAy8AGzDXNioJzlIM5ojua4C/DWKcQgghgmSmCyGEEFFh\nqJawhBBCRBlJWEIIIaKCPdwBDDal1Dzgh1rrZUqp6cCTmCuuHgIe1Fp7lFJfA+7B7ITx8+C6NaH3\nTwQ2A5la60hd2VQIIS45Q6qEpZT6BmbniZjgpt8BX9VaXwGUAg8ppaZijrmaDywD/i24jg1KqSTM\n5QDcgx27EEIMdUMqYQHFwCc5O3A4R2u9Ofh4I+aie5OA9Vprj9a6FXMRvvnBbu6/Ab6F2aNQCCHE\nIBpSCUtr/TJm9V/IEaXU4uDjlZgrl+4FFiulEpVS6cACzAX5/gN4U2u9J7i/BSGEEINmSCWsDtwP\nfEsp9S5QDlRqrQ8CvwDeAn6OOR6rErgbeEAptQ5z1c2IXMpaCCEuVUOu08V5bgDu1lpXK6WeANYo\npYZhTr+0SCmVjFlVuElrnR96k1KqBLg6PCELIcTQNFRLWKHR0kXAu0qpTcFtf9ZaVwJKKbUVsxT1\nDa11QyfvF0IIMUhkpgshhBBRYaiWsIQQQkQZSVhCCCGigiQsIYQQUUESlhBCiKggCUsIIURUkIQl\nhBAiKkjCEkIIERUkYQkhhIgK/x9Bq3P6z1JYAgAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Not bad. Still, it missed the two lower spikes. Let's go for the\n", "`detect_outlier_position_by_fft()`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def detect_outlier_position_by_fft(signal, threshold_freq=0.1,\n", " frequency_amplitude=.001):\n", " signal = signal.copy()\n", " fft_of_signal = np.fft.fft(signal)\n", " outlier = np.max(signal) if abs(np.max(signal)) > abs(np.min(signal)) else np.min(signal)\n", " if np.any(np.abs(fft_of_signal[threshold_freq:]) > frequency_amplitude):\n", " index_of_outlier = np.where(signal == outlier)\n", " return index_of_outlier[0]\n", " else:\n", " return None" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "outlier_idx = []\n", "\n", "y = df['u'].values\n", "\n", "opt = dict(threshold_freq=0.01, frequency_amplitude=0.001)\n", "\n", "win = 20\n", "for k in range(win*2, y.size, win):\n", " idx = detect_outlier_position_by_fft(y[k-win:k+win], **opt)\n", " if idx is not None:\n", " outlier_idx.append(k + idx[0] - win)\n", "outlier_idx = list(set(outlier_idx))\n", "\n", "fig, ax = plt.subplots(figsize=(7, 2.75))\n", "\n", "df['u'].plot()\n", "df['u'][outlier_idx].plot(**kw)\n", "_ = ax.set_ylim(-50, 50)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAawAAADJCAYAAAB/hgbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VOeZ6PHfmaY66hWERBMvvcg0A6bEuODu2M7GZRNn\nvUmcclPvluxms0l2N5vdxDe9JxsnjmPHLW5ZbGxTjE3vpr0gEAiEEOp9NO3cP86MLEAdSTODnu/n\n449mzjlz5mE8Z57zdsM0TYQQQohoZ4t0AEIIIUR/SMISQggREyRhCSGEiAmSsIQQQsQESVhCCCFi\ngiQsIYQQMcERqTf2+wNmfX1bpN5eCCFEFMrOdhs97YtYCcvhsEfqrYUQQsQgqRIUQggREyRhCSGE\niAmSsIQQQsQESVhCCCFigiQsIYQQMUESlhBCiJggCUsIIURMkIQlhBAiJkjCEkIIEROuaGompVQO\nsBu4HggCj4f+HgQ+o7WW5YyFEEIMiUGXsJRSTuAXQCtgAP8P+Cet9fLQ8zuHJEIhhBCCK6sS/A7w\nM6Ay9LxEa/126PFaYPWVBCaEEEJ0NaiEpZR6GKjWWq8LbTJC/4W1AKlXFpoQQgjxvsG2YX0MMJVS\nq4G5wO+A7C773UBDXyfJznYP8u2FEEKMNoNKWFrrFeHHSqkNwKPAd5RSK7TWm4A1wFt9nae6unkw\nby+EEOIq1VtBZqgWcDSBLwO/Ukq5gMPAc0N0biGEEALDNCPW89yUEpYQQoiuonLFYSGEEGIgJGEJ\nIYSICZKwhBBCxARJWEIIIWKCJCwhhBAxQRKWEEKImCAJSwghREyQhCWEECImSMISQggRE4ZqaqZR\nwbFrB7a6WgCCGZn45y+McERCCDF6SAmrnxy7dmCrrbVmTTTBVluLc9MGjKbGSIcmhBCjgiSsfgqX\nrI6W1/O7147S4fVjeDw49u6JcGRCCDE6SMIaoHM1rTS0dNDQ4o10KEIIMapIwuqnYEYmAIFAEAB/\nMIgZH49/XkkkwxJCiFFDElY/+ecvtBJUwFqOxeeIw7diFWZKaoQjE0KI0UES1gD455XQbnfidcRR\nN2VmpMMRQohRRbq1D4CZkoqeXMIhRz1LEpIjHY4QQowqUsIaIJ/fasPy+gIRjkQIIUYXSVgD5A0l\nLF+o84UQQoiRIQlrgMIlrPBfIYQQI0MS1gBJwhJCiMiQhDVAXr/VdiUJSwghRpYkrAGSEpYQQkTG\noLq1K6WcwP8ARUAc8O/AEeBxIAgcBD6jtTaHJszo0dlL0C+9BIUQYiQNtoT1IFCttV4O3Az8BHgM\n+KfQNgO4c2hCjB6maXb2EvRLCUsIIUbUYBPWs8DXupzDB5Rord8ObVsLrL7C2KJOeFomkCpBIYQY\naYOqEtRatwIopdxYyeurwHe7HNICXHWT7Pm6VAN6JWEJIcSIGvTUTEqpccALwE+01k8ppf67y243\n0NDXObKz3YN9+4ioa/J0PrbZbTEXvxBCxLLBdrrIBdYBn9Zabwht3quUWqG13gSsAd7q6zzV1c2D\nefuIqW5o73zc0uaNufiFECLa9VYQGGwJ65+wqvy+ppQKt2V9HvihUsoFHAaeG+S5o1bXakCfzCUo\nhBAjarBtWJ/HSlCXWnlF0US5rm1YMpegEEKMLBk4PABdewZKL0EhhBhZkrAGoGuVoPQSFEKIkSUJ\nawCkhCWEEJEjCWsAJGEJIUTkSMIagK6rDEvCEkKIkSUJawC69gyUhCWEECNLEtYA+HzvJ6mgaRII\nStISQoiRIglrAMJLihhG6LlPEpYQQowUSVgDEK4GTIyzxlvL4GEhhBg5krAGIJywkhKcgKyJJYQQ\nI0kS1gCEBwsnxTsuei6EEGL4ScIagMuqBCVhCSHEiBn0elijUXjy28R4Z+i5JCwRWxy7dmCrqwUg\nmJGJf/7CCEckRP9JCWsAfJdUCXadvV2IaOfYtQNbbS2YgAm22lqcmzZgNDVGOjQh+kVKWAMQbrOS\nEpaIJqZpUt3QzvGzjZgmLJiWQ5zTftlxtrpa2r1+3jlQSU1oMdLkBCe5R6tIvm0NRblunA65hxXR\nSxLWAFxewpKEJSKjvcPPiYpGDpbVsfd4NdUNns59z2wo5bZri7hhwTiM8KBBoLyqhXU7ymn1+HDa\nbWAYXGho52h1B/vrd+Ow25ivsrl7+USy0xIi8c8SoleSsAbA5w9iMwziXdbdq/QSFCOlzeNj++Eq\nDpyopbKujeqGdkzT2hfvsnONymbKuDSa27ys313B0+tLaWz1cu/KSRiGwd7j1by5u5ZUr5+ls/Ip\nmZKNzTBoCto4nDWJrPogh07Vse1wFbv0BW6YP447lk4gzmWnpd3H+t1n2Xq4igl5bj64fCJZ/Uho\nHb4Aza1efIEgeRmJFyXPaBIMmrz0ThlLZuWRm54Y6XBELyRhDYDXH8DptOEIVZtICUsMt5Z2H69u\nOcXGvRWdN0juRCdTCtKYXJCKGpeGKky/qCrvAyUF/Ncf97J2ezm1TR6y0xJ4bXs59nFT+XphM4Up\n1mVvxscTv2IVJUAJ1nRjO45U8fzGk6zdXs6OIxdIS3Zx6nwzgaCJ3WZQVdfGLn2Bu66byJpFhZcl\noVaPj0NldWw+UMnhsjpCOZW8jESum53P9PEZ5KQncKqyiTMXWmhu92EYsHBqLgU5ySPwiV6utKKR\nV7acosMX4MPXF0ckBtE/krAGwOcP4nLYOn8cZKYLMVy8vgBv7T7Lq1tP097hJzMljpXzxnLtjDwy\nUuJ7fW1achx/f/88/vuPe9hx5AIALoeNL9w7h/x0G+bePQD455Vc9DqbYbB4eh7zirN5+d0y1u04\nQ31zB4W5ySyansvyOWPYf6KGZ9aX8tzGE5yoaGTB1Bx8/iBl55s5fraBiurWzvNNyE8hLyMRnz/A\nvtJant14AjjRbcyvbjnNxDEpLJ8zhgVTc0iIG7mfpsZWL2DdHIjoJglrAHz+IE6HDZfDqhL0+aSX\noBicDl+Aw2V1HD5dT2OrF4/XT1ZqAunJLmoaPRwsq6O+uYOkeAcf/sBkVpUUDKhDRLo7jm8+spCz\n1a00t/kYk5lIVloCJuBbsarX18Y57dy3cjJrFhXhtNuIc73fgWPx9DymF2Xw85cOsvd4DXuP13Tu\nczltTCtKp7gglYXTchmTldS5r6Xdx/7SGo6fbaS6oZ2iPDcT8lNITXLR1Orlnfcqee9kLSfPNfHc\nxhN8/t7ZTBqb2v8P9Ao0ScKKGZKwBsDrD5IQ55ASlhgUj9fPbl3NnmPVHCqr67UN1OmwsWZRIbdc\nW0RSqFfqQDkddibkpww2XJITun/flCQXX/7wXPYcq6HV48NmGIzLSWZcTjIOe/dJNTnBydJZ+Syd\nld/t/vlTc6hr8rBp3zn+svU03316H5+7ZxbTxmcMOv7+am6zElarJKyoJwlrAHz+IKlJNquHFdKG\nJfrm8wc5Wl7PrqMX2HH0Ah1eq1Sen5lIyZRs5kzKIjs9AZfDRk2jh/pmD5mpCeSmJ/T44x8N7DYb\nC6bmDOk5M1LiuXv5RIry3Pz8pYP84LkDfOsTi/usAr1STW1WopISVvQb0oSllLIBPwVmAx3A32qt\nu6+0jkE+fwCnw4bTaf2QSC9B0Ztth87z5BvHaPX4AchMiePmhYUsnJZDfmbSZceHSymjXcmUbB66\nUfH42qO8/O4pHl4zdVjfr1mqBGPGUJew7gJcWuslSqlFwGOhbTEvGDTxB0yr04WUsEQvfP4Av/3f\no2w7XEWc084N88dRMiWL4oI0bLbo7NodbZbOyuO17eW8c6CSNYsLh7W7eVOoSrDN4ycYNOX/URQb\n6jqHpcBrAFrr7cD8IT5/xITbqxxdewnK1EziEsGgyS9fPsy2w1VMGpPCN/5mAfevLkYVpssP4QDY\nbTbuXj6RoGny0uayHo9rafdR2+jBDA9Ku0TQNHvcFxbudGECbR3+Qccsht9Ql7BSgKYuzwNKKZvW\n+rKiyF/eLSPo8+Ow22hp99Hc7qOlzUcwaJKc6GRsVhJzJmdFzVQx4dKUy2F/v5dgjJaw6ps72Hro\nPDcvLJQf0SFkmiZPrNPsPlbN1MI0vvihOTgdl0+RJPrnGpVNYU4y2w9XccviosvGaVXWtvJfT+6h\nqc1HujuOBVNzuHPZBBLiHDS2dPDm7rOs31NBcUEqn/3grB7bBMNtWGAlwJ46m4jIG+qE1QS4uzzv\nNlkB/PyFA32ezJ3o4oaFhdyxfCKZqZGdKsYIzb3mToojL9fqeWWz28nOdvf2sqi0ducZntt4gllT\ncihRQ9twPloFgiY/eXYfm/adY+LYVL7xySWdc06KwfvYHTP5xq+38Zft5Xz1bxZ1bq+qa+N7z+yn\nqc3H3OJsTp5rZN3OM+zS1YzJTuLwyVqCJthtBgdO1PKHN47z5QevuewGzecP0N6lVOWMc8bkNT1a\nDHXCehe4HXhWKbUY6DEr/eNHFlBR1YQ/EMSd4CQ50Yk7wYVhQHO7NVr+3fcqeWFjKS9vPsGyWfnc\nsWwCaclxQxxy/1TVtwEQCARoarQeN7d2UF3dHJF4rsSZ81YhuOxMPeMyZM64K+XzB/jVq0fYdfQC\nRbluPnfPLFqbPbQ2e/p+sehVYWYCkwtS2X7oPNv2n2XSmFTqmzv49pO7qWn0cN+qSaxZVITPH2Tt\nttO8uvUU9U0eJo1N5doZuSyYlssPnz/A2/sqSEl08sHlEy86f13Txf+PzlQ2kpkkNxqR1NsNw1An\nrD8DNyil3g09/1hPBy6dM4bq6p4DmzE+g7uvm8jWQ+f5322n2bjvHFsOneeOpRO4eVEhthGel8zn\nC1cJdm3DitEqwaYO629zR4QjiX3NbV5+/MJ7HD/byJSCVD537xwS42W0yFAxDIN7lk/kv/64lz+8\nfozblhTx581lVDd4uG3JeNYsKgKscWt3LJvAqpKxmKY1Vizs8/fO5l//Zwev7yhn1byxpLvfv+lt\nDlUHxrvseLwBGYsV5Yb0ytJam8Cnhup8ToeN5XPGsHRWHpsPVPLi5jKe23iC42ca+Pjt00e0ysXb\npQ3LbjMwjNgdOBxOVJKwrkxVXRvfe3Y/F+rbWTgth0dunSZtVsNAFaYzX2WzS1fzkz8fBGD1/ALu\nvm7CZce6E12XbUuKd3LH0gk8vvYof9l6ioduVJ37wj0E8zMTKatslq7tUS46ejT0wW6zsXLuWP7t\nkYVMH5/O/hO1fOfpfSPaSy/8Xg6HDcMwcDpsnaWuWGKaJnWSsK7YsTMN/Pvvd3Ghvp1bry3iE3fM\nkGQ1jB69cyZfeaiEmxcVcvfyidx/ffGAZn9fMjOP7LR43t5/jtrG96sBwz0Ew+PiJGFFt5hIWGHu\nRBdf+tBcrp2Rx+nzzTy9vnTE3vv9XoK20F97TJawmtt9+ENxS8IanGDQ5EfPH8DjDfCxNVO5Z8Wk\nEa+iHm1sNoPigjQ+tGoyty8ZP+ClShx2G3csnYA/YPLq1lOd27uWsECmZ4p2MVfZbrMZfORmRfmF\nZjbsqUCNS2PhtNxhf99wwgq3XzkdtpgchxVuvwKoH8WdAs5eaGH3sWqOn22gzeMnLyORCWNSWDw9\n96JqpQ5fgM37z7F+TwWzJ2Xy4euL6fAFaPX4mT0pk+vmjIngv0IMxLUz8vjL1tOhwchF5KQl0Nxq\nJagxUsKKCTGXsMCaTfrTd83kG4/v5I9vHGPWxMxhX47Ae0kJy2m30RGLCatLqarV48frC+DqZjn1\nq5HPH+DQqXo27q3gwInazu0Ou41T55vZdriKZ9aXMrUwjTFZydQ3e3jvZB0doVn5HXbrrj5csnZF\nyRhB0T82m8Gdyybwi5cP8cq7ZTxy6/T3S1hZkrBiQUwmLLDqnG9ZVMSL75Sxdnv5Zd1Vh5o3lJzC\n7RROp41WT+x9ucOlKpfThtcXpL6lI+ZWWfX5Axw+Vc/R8nqaWn14/QGyUxPIyUggNdFFYrwDwzDw\nB4K0tPu4UN/O8bONHDvT0Jl8poxLY/U1BUwtSicxzkFNk4d9x2vYfOAch07Vc+hUPQA56Qksnp7L\n+j0VnaVs/yWlbRE7FkzL4dUtp9hy8Dy3Xju+M2GlJ8cR77LT0i4zXUSzmE1YADctLGTjvgrW7Shn\n5dwxwzqr86U/Uk67LSa7tYc7XEzIS0GfaaChOfoTls8fpKq+jfKqZvYdr7mo1DMQuRmJzCvOYsHU\nnMuW3chJS+DGBeO4ccE42jw+KuvaSHA5yM+0lnbfeuh8Zyn70uphETtshlXK+umLB3nl3VM0t/qI\nc9qJc9lJTnDG5E3oaBLTCSvOZefu5RP57f8e5blNJ/jE7TOG7b0uqxJ02PD6g5imOeAG4EiqC7Vh\nTRxjJay6KO544fMHeHXLadZuL+/sKAJWqcdamiOTzNR4nA47F+rbuFDfTku7jzaPHxNrloPkBCfp\n7jgmjUkhtZ+DzhPjnUwac/HigU6HnfYObyiuUMKyj46q1KtNicpmTFYSO45U4XLacCdaw2OSEpyc\nq2nt49UikmI6YQEsnZnPhj0VbDtUxZIZecycmDks7xNOWOGlRcKJyx8wcTpiJ2GFqwTDJYyGKE1Y\ndU0evvPUXqrq20l3xzFrYiZjMhOZEfp76U1CapKL4oK0YYvHabd1tl2F/0oJKzbZDIObFo7jt/97\nlPaOQGeX9uQEJz5/kA5fgLhR0q4ba2I+YdlsBg+vmco3H9/F71/X/Nsjiy5a0nuw2jx+3jlwjmNn\nG6msbe2cxTm8tIizywS4sfTDVd/cQUqik+w0a0qmaC1hPfXmcarq21lVMpZ7V0wa9k41fXE6bZ3V\nwlIlGPsWT8/jhbdP0tjiJSXUKzQ86W1ru08SVpS6Kq64wlw3axYXUtPo4Yl1mmCw9+UE+lLd0M5/\nPLGLp9eXsudYNU2tXlrafCTGOcjNsNp7HDG4xIhpmtQ3d5CeEk96ilU9Fo1jsQ6cqGH3sWqmFKTy\n0A1TIp6swLpRCQRNAsGgJKyrgNNh44b54wA6qwSTQzPnSE/B6BX5X4IhcvuS8Rwsq2PLwfP4A0H+\n9rbpg1pi/MyFFh57ei9NbT5Wzy/gpgWFZKbGY5ompknnbM/hKsEn3zxOnNPG+bo2AgGT2ZMyWTA1\nh7HZ0bdybKvHj9cfJMMdhzvBicNuRF3C8voC/GHdMWyGwUM3qahpH+w6f6QkrKvDqnljOX2+mcUz\n8gBISrB+DmXwcPS6ahKWy2nn7z48jx88t58dRy5Q3eDhozcrCnP7v1RAfXMH33/WWrLgwRumcP01\nBZ37DMOaPzBsfJ6bLQfPs+voBcBq4Ac4db6Zl989xeLpudy9fGJn1Vs0CM9Mne6OwzAM0pLjom7w\n8LMbTlDT6OHmhYUURFHSvyhhSRvWVSEhzsGn7prZ+TxcJdjika7t0eqqSVgAifEOvvRXc/nda0fZ\ndqiKbz6+i4ljU8hNTyAYNGnvCDBrYgZLZ+VfNljW4/Xzw+cOUN/cwX0rJ12UrLqzev44ls7Kp6Xd\nWnQyMzUenz/IeydrWbutnG2Hq9hzrJq/ur6YlXPHREVJIVyaCs9WneGO43hFI4FgELst8j++B07U\n8Naes4zJSuKubiY2jSRXl4TlDXWpdw6iBC+iV2fCkhJW1LqqEhZYs2B84vYZLJmZx3MbTnCiopHS\ns42d+/eV1vDSO2Xct2oyS2bmYRgGbR4f33/2AKermrludj43Lyrs13slxDkual9x2G0snJbL/Kk5\nbD9UxR/fPMYTr2sOnqzlk3fMGNEZJYJBkxPnGjlUVse52jZqG9tpCk1Dk+G2xqulueMwTdh3vIY0\ndxyt7X4aWjqorG2lscXqwh0f52DSmBSKx6WRnRo/pIlXl9dTWdtGU5uXC/Xt7C+twWE3+MTt06Nu\n9g0pYV39kkNtWU+9eYz1u8+Sl5lITloCXl+QxjYvVXVt1DS2EwxaNSrZ6QnkpieQkugiJclFXkYi\nRXnuqKpViSSvL8Cp8814fQHyMhPJSIm/4jk3r7qEFTZzQiYzJ2Ti8wepbfJ03g1v3FfBm7vP8pu/\nHGHf8RoKcpLZrS9wtrqVxdNz+eshaDexGQbXzsxjalE6v3rlEHuP1/CD5w7wuXtnD3vvo8raVt7e\nf46th6o6Z6IGK5na7YY1JqnAGmMUvrDCSzb0ZOPeCgDSkl3MmpjJ7UvGk3UFF6XPH+DJN47x9v7K\ni7Y77DYeunFg1bgjJTzmStqwrl7FBWmsnDeW0+ebqKxto+KSMVkup42ctATsdqvH6LmaVk6fv3gB\n16mVmr+ak05xQSrBjEz88xeO5D8hKjS3eXnyjWPs1tUEunSAS4xzMLkglWWz8pk/dXArnRumeWU9\n6q6AGanVeqsb2vnVq4cvKnmtnDeWh26cMuSzbvv8QX7+0kH2Hq9BjUvj8/fNJt41PPcJ63aU86f1\npZhY1RslU7KYMzmLolx3Z7tVVy3tPnYeqaKx1UuHL0BygpOURBf5mUlkpFjHN7R0UHq2keNnGzh2\ntpGmVi8Ou8FNCwu5c9mEAXdsaWr18v1n93PqfDOFucnctLAQd6ibfVZqfFRUTXbnmfWlvLajnH/5\n6HyOnq7n2Y0n+MJ9s5k9KSvSoYlhYJomDS1eahrbiXPacSe6SE12XfT7EAyaNLR00NLuo6GlA+87\nW9m37SiYJvesmERuRiJmfDz+eSWYKam9vNvVodXjY9/xGp7beILGVi9js5KYMSGDhDgHlbWtnKps\n5kJDOwC3LSnirusmdvt7m53t7vFHeFQmLLC+bEdO12MzrKqx8ODB4eAPBPnly4fYpauZXJDKF++b\nM6RdtU3T5KV3ynj53VOkJbv48PXFzCvOHvISQNA02X64ihc2naC2qYPiglQ+dddM0vo5g0Rdk4fv\nPL2Pqro2ls7M469vUlFX9deTF94+yatbTvEPD8xDlzfw4jtl/N2H5zJtfEakQxNRwrVuLWXnmnhl\nSxmJcU4eWF1MYrwTMz4e34pVkQ5vWFTUtLL14HmOnW2g7FwTgaCJ3WbwweUTuWlhYWev6q7H/+j5\nA1yob+f6awp48IYpl52zt4Rl//rXvz7k/4h++npbm7fvo4aJYRjkpCeQnZbQ7SqlQ8lmMyhR2VTV\ntfHeyTp2Hb2AP2CSl5k4JD/Ymw9U8qf1pWSnxfMPD5QwZVxaZ6/FoWQYBuNyklk+ZwxV9e0cPFnH\n1oPnyc9MIi+j9/kIz1a38J2n91LT4GHNokIevHHKoIYdRMqJikaOnK5n0fRcLjRYk+munDt2WOev\nFLHFfqKU9OQ4HA4bJyoaaW33MbkgDRwOguP714koGDSpbvRw9kILp88309TqxZ3kirprJWiarNt5\nhp+9eBB9poH65g6K8tysKingr29SzC3O7rZpJSXRxaLpuew7XsOhsjqWzMq7bOX4pKS4b/T0vldt\nG1a0sdtsfOL2GaQkuti4r4JnNpSydvtpHr1zJtOK0gd93tpGD0+/dZyEODt/d/88slKHv8E33uXg\nU3fO4I2xqTy3sZQfPHeAD5SM5a8+UNxtqe5QWR0/ffE92jsC3LNiIrdeO37YYxxqMg5L9CWYkYmt\ntpZ5xdmUnm1En2mguDiPwluX9Pq66oZ2dukL7D1ew6nK5ovmzQSrg8d1s/N56CYVFQuFtrT7+OXL\nhzhYVkdKopP7V09h9qT+L/HkTnRx25Iifv3qEdbtPMMDqy8vZfVEEtYIstkMHrhhCncsm8DGvRW8\n9E4Z3316Lw/dqFg1b+yAz2eaJr9de8Ra+faWqSOSrMIMw+DGBeOYVpTOL18+xPo9FZw638yn75p5\nUanj7f3n+P1rGpvN4JN3zGDR9OFfbHM4dJewHJKwRBf++QtxbtqAzePh+msK+N2m03yvNosHznUw\n3eXjYFkdFdXvd+QwgZPnGjkcWsrGMKAwx82YrEQyU+NJTnBR1+RhX2kNG/edw2638cDq4hEfImOa\nJk1tPs7XtnKuto21205T0+hh1sRMHrl1GilJA6+hWjgtl+c3nWTz/kruXDaBpEtKWT2RhBUByQlO\nblsynqmF6fz4hQM89eYxpo9PH9AyH8Ggye9eO8rhU/XMnpTJsln5wxhxz8blJPPVj87vHPv2L7/Z\nwX0rJzE+383m/ZVs2FtBcoKT/3PPrGGdnHa4hXuZSglL9MY/rwTH3j1kxscz/R7Fvh1V/PTF3nvh\nTilIZcmsfOYWZ3XOa9jV7UvH8+0n9/DW7rOkJbuGtYbCNE0qalo5fqaBk5VNnK9to7K2rXMuVQAD\nuGvZBG5bOn7QJT6H3Zoa65kNpWzYU8FtS8b373WDejcxJCYXpPLQjYqfvniQP71Vyufund3r8Y5d\nO7DV1dLhC/LK0SY2+3IpynPzyK3TIjowOc5p5+O3TWfKuDSe3VDK71/XnftyMxL5wn2zo37Nrb6E\nZ+nvOg5LVhwWlzJTUjs7WFwPTJ8znte2l1NZ18bM8RkUF6Re1BEh3R1HTh/XRlK8ky99aC7femIX\nL2w6ydjsZOZOHvreqT5/kMfXHmXrofOd2+w2q61/alE6+ZmJ5GUkMiE/hTFZV95JbcXcMbyypYw3\nd5/lpoXjOicU740krAi7RmWjxqWxr7SGgydre1wexbFrB7baWk5UNLJxXwUt7T5uzW/m1lvvJX6Y\nO430h2EYrJw7ljmTsnjl3TI6fEHmFWcxa1LmVTHztYzDEoORn5nEx26ZdsXnSXfH8dkPzuY//7Cb\nX71yiK9+ZP6Q9mxu8/j44fPvcexMA+Pz3KycN5ZJY1PJTU8Ytg4fCXEOVs4dy9rt5Ww9VMXyOWP6\nfM2AI1FKpSqlXlFKbVRKbVFKLQ5tX6yU2qaUekcp9bVBxD8qGYbB/auLMQz45SsXjw3ryne+mtd3\nlPPq1lO0d/hZNC2X+xaPxX14/8gG3Id0dxwfuXkqH799OvOn5lwVyQq6tGEFgnj9MjWTGHlFeW4e\nXjOV9o4Av3jp0GWdM67EMxtKOXamgflTc/jHB0tYPmcMY7OShr134ur547DbDF7bXk6wH0OsBhPN\nF4E3tNYrgYeBn4S2/xy4X2u9DFiklJo7iHOPSoW5bj5yk6LN4+e/n9rLk+uO8eqWU2zcV8GOI1U8\n9eZxfr/bc0dXAAAOQUlEQVTuKEfL68lNT+CB1VNYPCMv6rq6Xs3CCcvrC3SuiyWdLsRIWzwjj2Wz\n8im/0MIbO88MyTnPXmhh84FKxmYl8ck7RnZatHR3HItn5HK+ro39pTV9Hj+YKsHvAeE1KZxAu1LK\nDbi01mWh7a8Dq4F9gzj/qLQiNKbnZy8e5K09Zy/bP9uVzHXjXSycltM5G0R4FL0Yfq4uJSyfP4jD\nbkRFF2Mx+nzoA5M5cKKGF98po0RlX3H78DMbSjFN67yRmGnmpoWFvPveeZ5ZX0pxQRrZvRzba8JS\nSj0CfOGSzQ9rrXcrpfKAJ4DPA6lAU5djmoGJg4h9VJs1MZPvfnoJVfXttLT7aGnz0dLuIz8zkalF\nK0l4ZxOGx1oO5GoePR+NLu3WLu1XIlKSE5w8cMMUfv7SIZ7bcILPfHDWoM918GQtB8vqmDE+nZkT\nIjNrS0F2MmsWF7J2Wzk//fN7fOfzK3o8tteEpbX+DfCbS7crpWYBTwFf1lpvVkqlAF1nLE0BGvoK\nNDs7+iY5jQZF43rYccMK2LnTerxgAaTK5zdSPKHmAofTQRCIczrk+ysi5pasZF7feYZ9pTU44p2k\nuwc+40ogaPLC73ZhGPDJe+aQk5MyDJH2z6P3zKWh1cfW9yp7PW7AVYJKqenAs8B9Wuv3ALTWTUop\nr1JqIlAG3Ah8va9zRXIuwdhkgzmLrIdeQD6/EdPSZE3a2dTswdPhx26T76+IrGun51pzF24qZc2i\nogG//u395zhV2cSy2fkkO20R/z5/5MYpOPuouBhMvca3ABfwQ6XUBqXUn0PbHwWeBLYDe7TWOwdx\nbiGiktN5cbd2Rz/GjAgxnMIdrzbvr2Sgk5h7vH7+/PZJXE4bd18XHa03cU47D6/pfQjAgEtYWuu7\neti+Hbh2oOcTIhZcOtOFO1HasERkJSc4uUZls/1wFaUVjf2eScbj9fOj59+jsdXLHUvHd65AHgvk\nqhOiH97vdBHAF5BOFyI6XDfbmpJtw56Kfh1f39zBY0/v48jpeuYVZ8XcRNQy04UQ/eCwGxiAV3oJ\niigytSidcTnJbD9cxW1Lxvc4ZVLQNNmwp4LnN53A4w2wZGYeH7tlatQumNoTSVhC9INhGDgdNto7\nQrNcSMISUcBmGNy1bAI/euE9Xn63jM+ObcNWVwtYy5345y/E6wvwm78cYefRCyTGOfjozYrr5oyJ\nyXGEkrCE6Cenw0Zbhw+QiW9F9JhbnEVRnpvGtzZTuzSH7DRrILGttpbmV9byy8oEDtUFKS5I5dN3\nzyJ1EMuBRAtJWEL0k8Nho9VjLbMgJSwRLQzD4O7rJvDulvU8s6GZ6UXppCRZa2kdLW/AYXexeM2N\nfOyWaTH/vZWEJUQ/Oe02Gr3ezsdCRItZEzOJv6aAnUeqOHCytnN7ujuOm66ZQNHt0yO6BNFQkYQl\nRD91nRQ01u9UxdXFMAymX1PM9KI0yqtaME1r+Y7s/HSCJddgXgXJCiRhCdFvXUtVMlO7iDb++Qtx\nbtrA+Lz3J8e+2uYblatOiH7qWqqSEpaIRv55JZjx8VftSg5SwhKiny5KWNKGJaKQmZJ61ZWqupKr\nToh+khKWEJElV50Q/XRxwpLJb4UYaZKwhOgnKWEJEVly1QnRT13braQNS4iRJ1edEP0kJSwhIkuu\nOiH6ydWl3UrmEhRi5MlVJ0Q/OaSEJUREyVUnRD9JlaAQkSVXnRD9JFMzCRFZctUJ0U8up/QSFCKS\n5KoTop8u6tYuJSwhRpxcdUL0k7RhCRFZctUJ0U8yNZMQkTXo2dqVUlOBbUCO1tqrlFoMfB/wA+u0\n1t8cohiFiApdk5S0YQkx8gZ11SmlUoDHAE+XzT8D7tdaLwMWKaXmDkF8QkQNqRIUIrIGfNUppQzg\nF8BXgPbQthQgTmtdFjrsdWD1UAUpRDTomqQc9qtjyXEhYkmvVYJKqUeAL1yy+TTwtNb6gFIKwABS\ngKYuxzQDE4cwTiEiLlwN6HTYMAxJWEKMtF4Tltb6N8Bvum5TSh0HHgklszys0tTtgLvLYSlAQ19v\nnp3t7usQIaKGJ2j9dTnt8t0VIgIG3OlCa10cfqyUKgNuDHW68CqlJgJlwI3A1/s6V3V180DfXoiI\naWlqB6zqQPnuCjE8ersZHHQvwRCzy+NHgScBO/C61nrnFZ5biKgSbsOSHoJCRMYVJSyt9cQuj7cD\n115xREJEqaQDu1lwcifp7jgcuxz45y+MdEhCjCpyqyhEPzh27cDVUA+micNmYKutxblpA0ZTY6RD\nE2LUkIQlRD/Y6mqx2wwM4/0qQcPjwbF3T4QjE2L0uNI2LCFGDcMwWDZrDGnJcZEORYhRSRKWEP0Q\nzMjEVltLyZTszm1mfDz+eSURjEqI0UWqBIXoB//8hZjx8Z3Pzfh4fCtWYaakRjAqIUYXSVhC9JN/\nXglmfLyUrISIEMM0zb6PGh6mDL4UQgjRVXa2u8d5z6SEJYQQIiZIwhJCCBETJGEJIYSICZKwhBBC\nxIRIdroQQggh+k1KWEIIIWKCJCwhhBAxQRKWEEKImCAJSwghREwY9ZPfKqWcwP8ARUAc8O/ACeCX\noUOOA3+rtQ5EJsKedRe71vqV0L4HgM9qrZdEMMQe9fC5nwVeBY6FDvuZ1vqZyETYux7i3w78CkgD\nDOAjWutTkYqxJz3E/gCQFzpkArBFa/1AZCLsWQ+xHwd+jbUC+jGs6zXqepP1EHs58HPAj/XveFRr\n7Y1YkD1QStmxvttTsD7nR4EO4HEgCBwEPjPcn7uUsOBBoFprvRy4GfgJ1hfpH7XWy0LH3B6p4Ppw\naew/BlBKzQP+JpKB9UN3n3sJ8JjWelXov6hMViHdxf9fwBNa6xXA14CZEYyvN5d9b7TW92utVwF3\nA/XAFyMZYC+6+9z/Fetm7TqsRHBrBOPrTXex/wr4Yij2CuDTEYyvN7cBwdBv4leBbwGPAf8U+vcY\nwJ3DHYQkLHgW68cFrM/DB9yjtX5HKeXCuutsiFRwfbgsdqVUBvAfwBewvkTRqrvP/RrgVqXUJqXU\nr5VSyRGLrm/dxb8UGKeUegPrx2l9hGLry6Wx+7vs+ybwQ6111YhH1T/dfe7tQKZSygDcQNSVUEK6\ni71Aa70ttG0LsCISgfVFa/0S8MnQ0/FYNzXXaK3fDm1bC6we7jhGfcLSWrdqrVuUUm6sL9Q/a61N\npVQhVjE3EzgQ0SB70E3sX8OqcvgS0BLR4PrQ3ecO7AD+b6iEchLrzjkqdRP/V7Eu5Dqt9Q1YVT3/\nEMEQe9TDZ49SKgf4AFY1T1TqIfYfAz8ADgM5wKYIhtijHr4zJ5VSy0OH3A4kRSzAPmitA0qpx7E+\n6ye5+Ia4BRj2tXZGfcICUEqNw7ob/r3W+mkArXW51noK8Avg/0Uyvt50jR2rDnwy8DPgKWC6Uiom\nYg997n/WWu8N7X4RmBex4PrhkvifAmqBl0O7XwHmRyq2vnT3nQfuBZ6MxvafrrqJ/Q/AdVrracAT\nWFVVUamb78zfAF9RSr0JVAE1kYyvL1rrhwGF1WYY32WXmxGoiRr1CUsplQusA/5ea/14aNvLSqnJ\noUNagKjrcAGXx6613qm1nhlqi/gwcFhr/aXIRtm97j534DWl1ILQ4+uBXZGIrT96iP8d3m8/WYFV\nQo86PcQO1me+NiJB9VMPsScC4bWKKrE6vUSdHmK/DXhQa70aqzbn9QiF1yul1F8rpb4SetqO9Zu4\nSykVrsJcA7zd7YuH0Kifmkkp9QPgPkB32fzPwH9j1YW3YvU6iro6/R5iX6O19iilxgN/jOJegt3F\n/o9Yd8c+rB+eT2ito7Jqs5v4TeBhrDvPJKy7zQe01o0RCbAXPcR+C9YNwhKtdVOkYutLD9+b7wL/\nAniweq59XGtdHoHwetVD7I9htRt2YFWJfyEaS7hKqQSsquI8wAn8J3AUq9OIC6s69uPDHfuoT1hC\nCCFiw6ivEhRCCBEbJGEJIYSICaMyYSmlViqlGpRSBV22fVsp9dFIxiWEEKJnozJhhXQAv+3yXBrz\nhBAiio3WhGVijYWoVUp9pusOpdSXlVI7lFJblFLfDm3bqZQqCj2+Vyn1/ZEPWQghRrfRmrDCI7Q/\nDXxRKTUp9NyN1e302lB38GKl1K3Ab4CPhI55mPcnxhVCCDFCRmvCAkBrXYc1597vsD6LeGBbl5nZ\nNwMzgD8C9yql8oEUrfXhSMQrhBCj2ahOWABa61exBvI9jDXwcJFSyh6aSHO5dYhuAnYD38eaq08I\nIcQIG60Jy+TiThZfwJpupAl4BngXa22jstAsxWCN6L4J+NMIximEECJEZroQQggRE0ZrCUsIIUSM\nkYQlhBAiJjgiHcBIU0otAr6ttV6llJoD/BxrxdXjwKNaa69S6kvAQ1idMH4UWrcm/PqpwDYgR2sd\nrSubCiHEVWdUlbCUUn+P1XkiLrTp18AXtdbXARXAp5VSM7HGXC0GVgH/HFrHBqVUCtZyAJ6Rjl0I\nIUa7UZWwgFLgg7w/cLhAa70t9HgL1qJ704CNWmuv1roDaxG+xaFu7r8AvoLVo1AIIcQIGlUJS2v9\nAlb1X9hJpdTy0OPbsVYufQ9YrpRKVkplAkuwFuT7V+AvWusDoeMNhBBCjJhRlbC68THgK0qpN4Eq\noEZrfRT4MfAa8COs8Vg1wIPAI0qpDVirbkblUtZCCHG1GnWdLi5xG/Cg1rpOKfVDYJ1SKgtr+qVl\nSqlUrKrCrVrr4vCLlFJlwI2RCVkIIUan0VrCCo+WPga8qZTaGtr2e611DaCUUjuwSlF/r7Vu7uH1\nQgghRojMdCGEECImjNYSlhBCiBgjCUsIIURMkIQlhBAiJkjCEkIIERMkYQkhhIgJkrCEEELEBElY\nQgghYoIkLCGEEDHh/wNR3d9Ika92IwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Not sure if this method is the best here... Maybe if the signal was\n", "contaminated by high frequency noise this method would perform better.\n", "\n", "Inspired by Bugra's median filter let's try a `rolling_median` filter using pandas." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from pandas import rolling_median\n", "\n", "threshold = 3\n", "df['pandas'] = rolling_median(df['u'], window=3, center=True).fillna(method='bfill').fillna(method='ffill')\n", "\n", "difference = np.abs(df['u'] - df['pandas'])\n", "outlier_idx = difference > threshold\n", "\n", "fig, ax = plt.subplots(figsize=figsize)\n", "df['u'].plot()\n", "df['u'][outlier_idx].plot(**kw)\n", "_ = ax.set_ylim(-50, 50)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAawAAADJCAYAAAB/hgbxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8XGeZ6PHfmaY6qlaxLUtu8uve4hbbcQGnOD0kYSEJ\nkBBKCHtJgLuwsCwL3F2WuguhhlACISSkkYrTsJ3iuMTdluzXki03WZbV+2jauX+cGVm21S1pZqzn\n+/nk45kzZ848msw5z3m7YZomQgghRLSzRToAIYQQoi8kYQkhhIgJkrCEEELEBElYQgghYoIkLCGE\nEDFBEpYQQoiY4IjUB/v9AbOurjVSHy+EECIKZWW5je5ei1gJy+GwR+qjhRBCxCCpEhRCCBETJGEJ\nIYSICZKwhBBCxARJWEIIIWKCJCwhhBAxQRKWEEKImCAJSwghREyQhCWEECImSMISQggREy5qaial\nVDawA/ggEAQeDf27H/i81lqWMxZCCDEoBlzCUko5gYeBFsAA/gf4utZ6Rej5TYMSoRBCCMHFVQn+\nEPgVUBF6Pl9r/Xbo8TpgzcUEJoQQQnQ2oISllLobqNJavx7aZIT+C2sGUi8uNCGEEOKsgbZh3QOY\nSqk1wFzgj0BWp9fdQH1vB8nKcg/w44UQQow0A0pYWuuV4cdKqQ3AfcAPlVIrtdZvAWuBf/R2nKqq\npoF8vBBCiEtUTwWZwVrA0QS+DDyilHIBxcAzg3RsIYQQAsM0I9bz3JQSlhBCiM6icsVhIYQQoj8k\nYQkhhIgJkrCEEELEBElYQgghYoIkLCGEEDFBEpYQQoiYIAlLCCFETJCEJYQQIiZIwhJCCBETBmtq\nphHBsX0bttoaAIIZmfgXLIpwREIIMXJICauPHNu3YaupsWZNNMFWU4PzrQ0YjQ2RDk0IIUYESVh9\nFC5ZHTxexx9fPUi714/h8eDYtTPCkQkhxMggCaufTlW3UN/cTn2zN9KhCCHEiCIJq4+CGZkABAJB\nAPzBIGZ8PP558yMZlhBCjBiSsPrIv2CRlaAC1nIsPkccvpWrMVNSIxyZEEKMDJKw+sE/bz5tdide\nRxy1U2ZGOhwhhBhRpFt7P5gpqejJ8yly1LE0ITnS4QghxIgiJax+8vmtNiyvLxDhSIQQYmSRhNVP\n3lDC8oU6XwghhBgekrD6KVzCCv8rhBBieEjC6idJWEIIERmSsPrJ67fariRhCSHE8JKE1U9SwhJC\niMgYULd2pZQT+D1QAMQB/wkcAB4FgsB+4PNaa3NwwoweHb0E/dJLUAghhtNAS1h3AlVa6xXANcAv\ngB8DXw9tM4CbBifE6GGaZkcvQb+UsIQQYlgNNGE9DXyz0zF8wHyt9duhbeuANRcZW9QJT8sEUiUo\nhBDDbUBVglrrFgCllBsreX0D+FGnXZqBS26SPV+nakCvJCwhhBhWA56aSSk1DngO+IXW+gml1A86\nvewG6ns7RlaWe6AfHxG1jZ6Oxza7LebiF0KIWDbQThc5wOvA/VrrDaHNu5RSK7XWbwFrgX/0dpyq\nqqaBfHzEVNW3dTxubvXGXPxCCBHteioIDLSE9XWsKr9vKqXCbVkPAA8ppVxAMfDMAI8dtTpXA/pk\nLkEhhBhWA23DegArQZ1v1UVFE+U6t2HJXIJCCDG8ZOBwP3TuGSi9BIUQYnhJwuqHzlWC0ktQCCGG\nlySsfpASlhBCRI4krH6QhCWEEJEjCasfOq8yLAlLCCGGlySsfujcM1ASlhBCDC9JWP3g851NUkHT\nJBCUpCWEEMNFElY/hJcUMYzQc58kLCGEGC6SsPohXA2YGGeNt5bBw0IIMXwkYfVDOGElJTgBWRNL\nCCGGkySsfggPFk6Kd5zzXAghxNCThNUPF1QJSsISQohhIwmrH8KT3ybGO0PPJWEJIcRwkYTVD77z\nqgQ7z94uhBBiaA14xeGRKNxmJSUsEU1M06Sqvo2Skw2YJiyclk2c097lvs1tPp5aX8qJM80ApLvj\nKByXSmFeGgU5bpyOkXcP69i+DVttDQDBjEz8CxZFOCLRHUlY/XBhCUsSloiMtnY/h8sb2F9Wy66S\nKqrqPR2vPbWhlOsvL+DKheMwwoMGgeKjtTzycjENzV5cDhuGYXCssondpdUAOOw2Fqgsblkxkay0\nhGH/myLBsX0btpqajue2mhqcb23AP28+ZkpqBCMTXZGE1Q8+fxCbYRDvsu5epZegGC6tHh9biyvZ\ne7iGitpWqurbME3rtXiXnctUFlPGpdHU6mX9jnKeXF9KQ4uX21ZNwjAMdpVU8cu/7QfgtlWTuGZR\nPjabQW2jh9LyBkpONFB0tJYtxZVs12e4csE4blw2gTiXneY2H+t3nGRzcSUTct18aMVERvUhobX7\nAjS1ePEFguRmJJ6TPKOFrbaGoGmytbiSaQXppCXHYXg8OHbtxLdydaTDE+eRhNUPXn8Ap9OGI1Rt\nIiUsMdSa23y8/N5RNu4q77hBcic6mZKXxuS8VNS4NFR++jlVeR+Yn8f3/7KLdVuPU9PoISstgVe3\nHsduN/ji7XNQ+ekd+2akxLMoJZ5F03IImibbDlTy7MYjrNt6nG0HzpCW7OLo6SYCQRO7zaCytpXt\n+gw3XzGRtYvzL0hCLR4fRWW1vLO3guKyWkI5ldyMRK6YPZrp4zPITk/gaEUjJ84009TmwzBg0dQc\n8rKTh/z77EpFTQvbDlTi8wdZMWdMRGIQfSMJqx98/iAuh63j4iAzXYih4vUF+MeOk7y8+Rht7X4y\nU+JYNW8sl8/IJSMlvsf3piXH8ZWPzuMHf9nJtgNnAHA5bDx427nJ6nw2w2DJ9FzmFWbx4qYyXt92\ngrqmdvJzklk8PYcVc8aw53A1T60v5ZmNhzlc3sDCqdn4/EHKTjdRcrKe8qqWjuNNGJ1CbkYiPn+A\n3aU1PL3xMHC4y89++b1jTByTwoo5Y1g4NZuEuOG5NAUzMmk9UQ+Ax+sHwIyPxz9v/rB8vugfSVj9\n4PMHcTpsuBxWlaDPJ70ExcC0+wIUl9VSfKyOhhYvHq+fUakJpCe7qG7wsL+slrqmdpLiHXzkA5NZ\nPT+vXx0i0t1xfOfeRZysaqGp1ceYzMQ+VeMBxDnt3L5qMmsXF+C024hzne3AsWR6LtMLMvj1C/vZ\nVVLNrpLqjtdcThvTCtIpzEtl0bQcxoxK6nituc3HntJqSk42UFXfRkGumwmjU0hNctHY4uXdfRXs\nO1LDkVONPLPxMA/cNptJY4e+Dcm/YBFNu48C4PEGMOPjpSowiknC6gevP0hCnENKWGJAPF4/O3QV\nOw9VUVRW22MbqNNhY+3ifK69vICkUK/U/nI67EwYnTLQcElO6PpzU5JcfPkjc9l5qJoWjw+bYTAu\nO5lx2ck47F0n1eQEJ8tmjWbZrNFdvr5gaja1jR7e2n2KVzYf40dP7uYLt85i2viMAcffV8fyFN4d\n5TQF7VKyinKSsPrB5w+SmmTDaZc2LNE3Pn+Qg8fr2H7wDNsOnqHda5XKR2cmMn9KFnMmjSIrPQGX\nw0Z1g4e6Jg+ZqQnkpCd0e/GPBnabjYVTswf1mBkp8dyyYiIFuW5+/cJ+fvrMXr77mSW9VoFerBoj\nnj0FczidnsDt0jMwqg1qwlJK2YBfArOBduBTWuuuK61jkM8fwOmw4XRaFxLpJSh6sqXoNI+/cYgW\nj9U2kpkSxzWL8lk0LZvRmUkX7B8upYx086dkcddVikfXHeTFTUe5e+3UIf28phYvYFVbiug22CWs\nmwGX1nqpUmox8OPQtpgXDJr4A6bV6SKGS1gySHLo+fwB/vD3g2wpriTOaefKBeOYP2UUhXlp2GzR\n17U7Gi2blcurW4/z7t4K1i7JJyc9ccg+q7HVSlitHj/BoCn/j6LYYNc5LANeBdBabwUWDPLxIybc\nXuXo3EswxqZm6hgkaQLm2UGSRmNDpEO7ZASDJr95sZgtxZVMGpPCtz+5kI+uKUTlp8uFsB/sNhu3\nrJhI0DR54Z2ybvdrbvNR0+DBDA9KO0/QNLt9LawxVMIygdZ2/4BjFkNvsEtYKUBjp+cBpZRNa31B\nUeSVTWUEfX4cdhvNbT6a2nw0t/oIBk2SE52MHZXEnMmjomaqmHBpyuWwn+0lGGMlrHDJqrnNx8Hj\ndcyfkoVNBkkOGtM0eex1zY5DVUzNT+OLH56D09H1FEmid5epLPKzk9laXMm1SwouGKdVUdPC9x/f\nSWOrj3R3HAunZnPT8gkkxDloaG7nzR0nWb+znMK8VP75Q7O6bRNsbD1bFdjc5uu2s4mIvMFOWI2A\nu9PzLpMVwK+f29vrwdyJLq5clM+NKyaSmRrZqWKM+jYrpqQ4cnOsnlc2u52sLHdPb4suqYlgmuwo\nqWLTvgrGj0llUl4aJCRALP0dUSgQNPnF07t5a/cpJo5N5dufXdox56QYuHtunMm3f7uFV7Ye5xuf\nXNyxvbK2lf99ag+NrT7mFmZx5FQDr79/gu26ijFZSRQfqSFogt1msPdwDX9+o4Qv33nZBaVcnz9A\nW6dSlTPOGVvn9Agz2AlrE3AD8LRSagnQbVb6148vpLyyEX8giDvBSXKiE3eCC8OApjZrtPymfRU8\nt7GUF985zPJZo7lx+QTSkuMGOeS+qaxrBSAQCNDYYD1uammnqqopIvEMhMOegK2mhurQ31JZ3Uxm\nVir+6VMxY+jviDY+f4BHXj7A9oNnKMhx84VbZ9HS5KGlydP7m0WP8jMTmJyXytai02zZc5JJY1Kp\na2rne4/voLrBw+2rJ7F2cQE+f5B1W47x8uaj1DV6mDQ2lctn5LBwWg4PPbuXt3eXk5Lo5EMrJp5z\n/NrGc/8fnahoIDNJbjQiqacbhsFOWH8DrlRKbQo9v6e7HZfNGUNVVfeBzRifwS1XTGRz0Wn+vuUY\nG3ef4r2i09y4bALXLM7HNszzkvl84SrBzm1YsVUl6F+wCOdbG2hute4o6wN2qQq8SE2tXn7+3D5K\nTjYwJS+VL9w2h8R4GS0yWAzD4NYVE/n+X3bx59cOcf3SAv72ThlV9R6uXzqetYsLAGvc2o3LJ7B6\n/lhM0xorFvbAbbP5j99v47Vtx1k9byzp7rM3vU2h6sB4lx2PN0CL9BSMaoN6ZmmtTeBzg3U8p8PG\nijljWDYrl3f2VvD8O2U8s/EwJSfq+fQN04e1ysXbqQ3LbjMwjNgcOOyfN5/al4vxOuIoyS1kYaQD\nimGVta3879N7OFPXxqJp2dx73TRpsxoCKj+dBSqL7bqKX4Qm8F2zII9brphwwb7uRNcF25Lindy4\nbAKPrjvIK5uPctdVquO1cA/B0ZmJlFU0Sdf2KBcTt4J2m41Vc8dy2ZQsHn6xiD2Ha/jhk7v5+l3z\nh+0CEe4R6Agty+B02DpKXbEk6E5hU/YM/JlBZgWk6mOgDp2o52fP7qXF4+e6ywu4ZcXEYS/1jyT3\n3TSTw6ca2FVSTUKcg+svL+jX7O9LZ+byyuajvL3nFGsXF5CZag1GDvcQHJ2ZJAkrBkRHF7w+cie6\n+NKH53L5jFyOnW7iyfWlw/bZZ3sJ2kL/2mOyhNXU5sMfiruuqT3C0cSmYNDkZ8/uxeMNcM/aqdy6\ncpIkqyFmsxkU5qXx4dWTuWHp+H4vVeKw27hx2QT8AZOXNx/t2N65hAVIlWCUi4kSVmc2m8HHr1Ec\nP9PEhp3lqHFpLJqWM+SfG05Y4fYrp8MWc+OwAOoazyapuhHcKeDkmWZ2HKqi5GQ9rR4/uRmJTBiT\nwpLpOedUK7X7Aryz5xTrd5Yze1ImH/lgIe2+AC0eP7MnZXKFLEcRMy6fkcsrm4+FBiMXkJ2WQFOL\nlaDGhGYekRJWdIu5hAXWbNL33zyTbz/6Pn954xCzJmYO+XIE3vNKWE67jfZYTFidSlUtHj9eXwBX\nN8upX2p8/gBFR+vYuKucvYfPrjLrsNs4erqJLcWVPLW+lKn5aYwZlUxdk4d9R2ppD83K77Bbd/Xh\nkrUrSsYIir6x2QxuWj6Bh18s4qVNZdx73fSzJaxRkrBiQUwmLLDqnK9dXMDz75axbuvxC7qrDjZv\nKDmF28ycThstntj7cYdLVS6nDa8vSF1z+5BOezMUfP4AxUfrOHi8jsYWH15/gKzUBLIzEkhNdJEY\n78AwDPyBIM1tPs7UtVFysoFDJ+o7ks+UcWmsuSyPqQXpJMY5qG70sLukmnf2nqLoaB1FR+sAyE5P\nYMn0HNbvLO8oZfvPK22L2LFwWjYvv3eU9/af5rrLx3ckrPTkOOJddprbZKaLaBazCQvg6kX5bNxd\nzuvbjrNq7pghndX5/IuU026LuW7tALWhEtaE3BT0iXrqm6I/Yfn8QSrrWjle2cTukupzSj39kZOR\nyLzCUSycmn3BshvZaQlctXAcVy0cR6vHR0VtKwkuB6MzraXdNxed7ihln189LGKHzbBKWb98fj8v\nbTpKU4uPOKedOJed5ARnTN6EjiQxnbDiXHZuWTGRP/z9IM+8dZjP3DBjyD7rgipBhw2vP4hpmv1u\nAI6k2lAb1sQxVsKqjeKOFz5/gJffO8a6rcc7OoqAVeqxlubIJDM1HqfDzpm6Vs7UtdHc5qPV48fE\nmuUgOcFJujuOSWNSSO3joPPEeCeTxpy7zITTYaet3RuKK5Sw7COjKvVSM19lMWZUEtsOVOJy2nAn\nWr1lkxKcnKpu6eXdIpJiOmEBLJs5mg07y9lSVMnSGbnMnJg5JJ8TTljhpUXCicsfMHE6YidhhasE\nwyWM+ihNWLWNHn74xC4q69pId8cxa2ImYzITmRH69/ybhNQkF4V5aUMWj9Nu62i7Cv8rJazYZDMM\nrl40jj/8/SBt7YGOpV6SE5z4/EHafQHiRki7bqyJ+YRlsxncvXYq33l0O396TfP/7l18zpLeA9Xq\n8fPu3lMcOtlARU1LxyzO4aVFnJ0mwI2lC1ddUzspiU6yQsulR2sJ64k3S6isa2P1/LHctnLSkHeq\n6Y3TaeuoFpYqwdi3ZHouz719hIZmLymhXqHhSW9b2nySsKLUJXHG5ee4Wbskn+oGD4+9rgkGe15O\noDdV9W3812PbeXJ9KTsPVdHY4qW51UdinIOcDKu9xxGDS4yYpkldUzvpKfGkp1jVY9E4Fmvv4Wp2\nHKpiSl4qd105JeLJCqwblUDQJBAMSsK6BDgdNq5cMA6go0owOTRzjvQUjF6RvxIMkhuWjmd/WS3v\n7T+NPxDkU9dPH9AS4yfONPPjJ3fR2OpjzYI8rl6YT2ZqPKZpYpp0zPYcrhJ8/M0S4pw2Tte2EgiY\nzJ6UycKp2YzNir6VY1s8frz+IBnuONwJThx2I+oSltcX4M+vH8JmGNx1tYqa9sHO80dKwro0rJ43\nlmOnm1gyIxeApATrciiDh6PXJZOwXE47//KRefz0mT1sO3CGqnoPn7hGkZ/T96UC6pra+cnT1pIF\nd145hQ9eltfxmmFY8weGjc91897+02w/eAawGvgBjp5u4sVNR1kyPYdbVkzsqHqLBuGZqdPdcRiG\nQVpyXNQNHn56w2GqGzxcsyifvChK+uckLGnDuiQkxDn43M0zO56HqwSbPdK1PVpdMgkLIDHewZf+\naS5/fPUgW4oq+c6j25k4NoWc9ASCQZO29gCzJmawbNboCwbLerx+HnpmL3VN7dy+atI5yaoraxaM\nY9ms0TS3WYtOZqbG4/MH2XekhnVbjrOluJKdh6r4pw8WsmrumKgoKYRLU+HZqjPccZSUNxAIBrHb\nIn/x3Xu4mn/sPMmYUUnc3MXEppHk6pSwvKEu9c4BlOBF9OpIWFLCilqXVMICaxaMz9wwg6Uzc3lm\nw2EOlzdQevLsEvC7S6t54d0ybl89maUzczEMg1aPj588vZdjlU1cMXs01yzO79NnJcQ5zmlfcdht\nLJqWw4Kp2WwtquQvbx7isdc0+4/U8NkbZwzrjBLBoMnhUw0UldVyqqaVmoY2GkPT0GS4rfFqae44\nTBN2l1ST5o6jpc1PfXM7FTUtNDRbXbjj4xxMGpNC4bg0slLjBzXx6uN1VNS00tjq5UxdG3tKq3HY\nDT5zw/Som31DSliXvuRQW9YTbx5i/Y6T5GYmkp2WgNcXpKHVS2VtK9UNbQSDVo1KVnoCOekJpCS6\nSElykZuRSEGuO6pqVSLJ6wtw9HQTXl+A3MxEMlLiL3rOzUsuYYXNnJDJzAmZ+PxBaho9HXfDG3eX\n8+aOk/zulQPsLqkmLzuZHfoMJ6taWDI9h48NQruJzTC4fGYuUwvSeeSlInaVVPPTZ/byhdtmD3nv\no4qaFt7ec4rNRZUdM1GDlUztdsMak5RnjTEKn1jhJRu6s3FXOQBpyS5mTczkhqXjGXURJ6XPH+Dx\nNw7x9p6Kc7Y77Dbuuqp/1bjDJTzmStqwLl2FeWmsmjeWY6cbqahppfy8MVkup43stATsdqvH6Knq\nFo6dPnfhU8OAz900kwVTs4cz9KjS1Orl8TcOsUNXEejUAS4xzsHkvFSWzxo94O/HMM2L61F3EcxI\nrdZbVd/GIy8Xn1PyWjVvLHddNWXQZ932+YP8+oX97CqpRo1L44HbZxPvGpr7hNe3Heev60sxsao3\n5k8ZxZzJoyjIcXe0W3XW3Obj/QOVNLR4afcFSE5wkpLoYnRmEhkp1v71ze2Unmyg5GQ9h0420Nji\nxWE3uHpRPjctn9Dvji2NLV5+8vQejp5uIj8nmasX5eMOdbMflRofFVWTXXlqfSmvbjvOv39iAQeP\n1fH0xsM8ePtsZk8aFenQxBAwTZP6Zi/VDW3EOe24E12kJrvOuT4Egyb1ze00t/mob27nVHUrL2wq\nwwyafPXO+RfMpnKpa/H42F1SzTMbD9PQ4mXsqCRmTMggIc5BRU0LRyuaOFPfBsD1Swu4+Yqul+TJ\nynJ3exEekQkLrB/bgWN12Ayraiw8eHAo+ANBfvNiEdt1FZPzUvni7XMGtau2aZq88G4ZL246Slqy\ni498sJB5hVmDXgIImiZbiyt57q3D1DS2U5iXyudunklaH2eQqG308MMnd1NZ28qymbl87GoVdVV/\n3Xnu7SO8/N5RvnrHPPTxep5/t4x/+chcpo3PiHRoIorsKa3moWf3kpLo4tufXHTOyseXovLqFjbv\nP82hk/WUnWokEDSx2ww+tGIiVy/K7+hV3Xn/nz27lzN1bXzwsjzuvHLKBcfsKWHZv/Wtbw36H9FH\n32pt9fa+1xAxDIPs9ASy0hK6XKV0MNlsBvNVFpW1rew7Usv2g2fwB0xyMxMH5YL9zt4K/rq+lKy0\neL56x3ymjEvr6LU4mAzDYFx2MivmjKGyro39R2rZvP80ozOTyM3oeT7Ck1XN/PDJXVTXe1i7OJ87\nr5oyoGEHkXK4vIEDx+pYPD2HM/XWZLqr5o4d0vkrRezJzUjE5bSx81A1dc3tLFD9q/oKBk2qGjyc\nPNPMsdNNNLZ4cSe5ou5cCZomr79/gl89vx99op66pnYKct2snp/Hx65WzC3M6rJpJSXRxeLpOewu\nqaaorJals3IvWDk+KSnu29197iXbhhVt7DYbn7lhBimJLjbuLuepDaWs23qM+26aybSC9AEft6bB\nw5P/KCEhzs6/fHQeo1KHvsE33uXgczfN4I2xqTyzsZSfPrOXD8wfyz99oLDLUl1RWS2/fH4fbe0B\nbl05kesuHz/kMQ42GYcl+urqhfns0FVsLa5kyfQc5kzuudq4qr6N7foMu0qqOVrRdM68mWB18Lhi\n9mjuulpFxUKhzW0+fvNiEfvLaklJdPLRNVOYPanvSzy5E11cv7SA3758gNffP8Eday4sZXVHEtYw\nstkM7rhyCjcun8DGXeW88G4ZP3pyF3ddpVg9b2y/j2eaJn9Yd8Ba+fbaqcOSrMIMw+CqheOYVpDO\nb14sYv3Oco6ebuL+m2eeU+p4e88p/vSqxmYz+OyNM1g8fegX2xwKXSUshyQs0YXwdHHf/sP7/PHV\ng9yxZgrTx6ezv6yW8qqzHTlM4MipBopDS9kYBuRnuxkzKpHM1HiSE1zUNnrYXVrNxt2nsNtt3LGm\ncNiHyJimSWOrj9M1LZyqaWXdlmNUN3iYNTGTe6+bNqBqz0XTcnj2rSO8s6eCm5ZPIOm8UlZ3JGFF\nQHKCk+uXjmdqfjo/f24vT7x5iOnj0/u1zEcwaPLHVw9SfLSO2ZMyWT5r9BBG3L1x2cl84xMLOsa+\n/fvvtnH7qkmMH+3mnT0VbNhVTnKCk/9z66whnZx2qIV7mUoJS/RFXlYyt66cxNMbSvnl8z33wp2S\nl8rSWaOZWziqY17Dzm5YNp7vPb6Tf+w4SVqya0hrKEzTpLy6hZIT9RypaOR0TSsVNa0dc6kCGMDN\nyydw/bLxAy7xOezW1FhPbShlw85yrl86vm/vG9CniUExOS+Vu65S/PL5/fz1H6V84bbZfXpfW7uf\nP/z9ANt1FQW5bu69blpEBybHOe18+vrpTBmXxtMbSvnTa7rjtZyMRB68fXbUr7nVm/As/Z3HYcmK\nw6In1yzOZ87kTF7depyK2lZmjs+gMC/1nI4I6e44sns5N5LinXzpw3P57mPbee6tI4zNSmZuL9WM\nA+HzB3l03UE2F53u2Ga3WW39UwvSGZ2ZSG5GIhNGpzBm1MV3Uls5dwwvvVfGmztOcvWicR0TivdE\nElaEXaayUOPS2F1azf4jNb0uj7LzUBWPv3GIuqZ2poxL44HbZkfF5LCGYbBq7ljmTBrFS5vKaPcF\nmVc4ilmTMi+Jma9lHJYYiNGZSdxz7bSLPk66O45//tBs/vvPO3jkpSK+8fEFg9qzudXj46Fn93Ho\nRD3jc92smjeWSWNTyUlPGLIOHwlxDlbNHcu6rcfZXFTJijljen1Pv690SqlU4M+AG3ABX9Jab1FK\nLQF+AviB17XW3+nvsUciwzD46JpCvv3o+/zmpWK+cOtsJuelXrBfq8fP428cYnPRaRx2gxuXjee6\ny8dH3UUz3R3Hx6+ZGukwBl1HG1YgiNcvUzOJ4VeQ6+butVP5zUvFPPxCEd/4xIJBSyZPbSjl0Il6\nFkzN5lPXTRu24SZrFozj9fdP8OrW4yyfPbrXKsaB3Jp/EXhDa/2QUmoK8ARwGfBr4BatdZlS6hWl\n1Fyt9e7gTGfYAAAOHUlEQVQBHH/Eyc9x8/GrFY+9dogfPLGLlXPGkJrsIjnRSWKcg8PljWw7WElD\ns5fxuW4+df30QSmSi74LJyyvL9CxLpZ0uhDDbcmMXIqP1vHuvgreeP8Ea5cUXPQxT55p5p29FYwd\nlcRnb5w+rIP3091xLJmRw6Z9p9lTWs28wqwe9x9IwvpfILwmhRNoU0q5AZfWuiy0/TVgDSAJq49W\nhsb0/Or5/fxj58kLXk+Is3PD0vHcsGx81I3JGAlcnUpYPn8Qh92Iii7GYuT58Acms/dwNc+/W8Z8\nlXXR7cNPbSjFNK3jRmKmmasX5bNp32meWl9KYV4aPaWsHhOWUupe4MHzNt+ttd6hlMoFHgMeAFKB\nxk77NAETBxD7iDZrYiY/un8plXVtNLf5aG710dzmY3RmIlML0iVRRdD53dqjrSpWjBzJCU7uuHIK\nv36hiGc2HObzH5o14GPtP1LD/rJaZoxPZ+aEyMzakpeVzNol+azbcpxf/m0fP3xgZbf79piwtNa/\nA353/nal1CysqsAva63fUUqlYLVphaUA9b0FmpUVfZOcRoOCcZGOQJzPExrL6XA6CAJxTof8fkXE\nXDsqmdfeP8Hu0moc8U7S3f2fcSUQNHnuj9sxDPjsrXPIzo7c3If33TqX+hYfm/dV9LjfQDpdTAee\nBm7XWu8D0Fo3KqW8SqmJQBlwFfCt3o4VybkEheiP5kZr0s7GJg+edj92m/x+RWRdPj2HslONvPRW\nKWsX978t6+09pzha0cjy2aNJdtoi/nv++FVTcPZScTGQeo3vYvUOfEgptUEp9bfQ9vuAx4GtwE6t\n9fsDOLYQUcnpPLdbu6MPY0aEGEpLZuTisNt4Z08F/Z3E3OP187e3j+By2rjliuhovYlz2rl7bc9D\nAPpdwtJa39zN9q3A5f09nhCx4PyZLtyJ0oYlIis5wcllKoutxZWUljf0eSYZj9fPz57dR0OLlxuX\nje9YgTwWyFknRB+c7XQRwBeQThciOlwx25qSbcPO8j7tX9fUzo+f3M2BY3XMKxwVcxNRR36KBCFi\ngMNuYABe6SUoosjUgnTGZSeztbiS65eO73Z8ZtA02bCznGffOozHG2DpzFzuuXZq1C6Y2p3YilaI\nCDEMA6fDRlt7aJYLSVgiCtgMg5uXT8AEXtxU1uU+Xl+Ah18o4vE3DmEzDD5xjeKT102LuWQFUsIS\nos+cDhut7T5AJr4V0WNu4SgKct1sO3CGa5c0kZ9zdrhFRU0Lv//7AQ6XN1KYl8r9t8wiNYZXQZaE\nJUQfORw2WjzWMgtSwhLRwjAMbrliAj95ei//9dgOls8azajUeE5Vt7C5qJKgabJkeg73XDst5n+3\nkrCE6COn3UaD19vxWIhoMWtiJnevncrL7x1lw66zHTByMxK5bdUk5hWOiugSRINFEpYQfdR5ButY\nv1MVlxbDMFgxZwzLZlmT4waDJu5EFwW5yTHZVtUdSVhC9FHnUpXM1C6ikd1mY1Yva+rFMjnrhOij\nzqUqKWEJMfzkrBOij85JWNKGJcSwk7NOiD6SEpYQkSVnnRB9dG7CkslvhRhukrCE6CMpYQkRWXLW\nCdFHndutpA1LiOEnZ50QfSQlLCEiS846IfrI1andSuYSFGL4yVknRB85pIQlRETJWSdEH0mVoBCR\nJWedEH0kUzMJEVly1gnRRy6n9BIUIpLkrBOij87p1i4lLCGGnZx1QvSRtGEJEVly1gnRRzI1kxCR\nNeD1sJRSU4EtQLbW2quUWgL8BPADr2utvzNIMQoRFTonKWnDEmL4DeisU0qlAD8GPJ02/wr4qNZ6\nObBYKTV3EOITImpIlaAQkdXvs04pZQAPA18D2kLbUoA4rXVZaLfXgDWDFaQQ0aBzknLYjQhGIsTI\n1GOVoFLqXuDB8zYfA57UWu9VSgEYQArQ2GmfJmDiIMYpRMSFqwGdDhuGIQlLiOHWY8LSWv8O+F3n\nbUqpEuDeUDLLxSpN3QC4O+2WAtT39uFZWe7edhEianiC1r8up11+u0JEQL87XWitC8OPlVJlwFWh\nThdepdREoAy4CvhWb8eqqmrq78cLETHNjW2AVR0ov10hhkZPN4MD7iUYYnZ6fB/wOGAHXtNav3+R\nxxYiqoTbsKSHoBCRcVEJS2s9sdPjrcDlFx2REFEqae8OFh55n3R3HI7tDvwLFkU6JCFGFLlVFKIP\nHNu34aqvA9PEYTOw1dTgfGsDRmNDpEMTYsSQhCVEH9hqa7DbDAzjbJWg4fHg2LUzwpEJMXJcbBuW\nECOGYRgsnzWGtOS4SIcixIgkCUuIPghmZGKrqWH+lKyObWZ8PP558yMYlRAji1QJCtEH/gWLMOPj\nO56b8fH4Vq7GTEmNYFRCjCySsIToI/+8+Zjx8VKyEiJCDNM0e99raJgy+FIIIURnWVnubuc9kxKW\nEEKImCAJSwghREyQhCWEECImSMISQggREyLZ6UIIIYToMylhCSGEiAmSsIQQQsQESVhCCCFigiQs\nIYQQMWHET36rlHICvwcKgDjgP4HDwG9Cu5QAn9JaByITYfe6il1r/VLotTuAf9ZaL41giN3q5ns/\nCbwMHArt9iut9VORibBn3cS/FXgESAMM4ONa66ORirE73cR+B5Ab2mUC8J7W+o7IRNi9bmIvAX6L\ntQL6IazzNep6k3UT+3Hg14Af6++4T2vtjViQ3VBK2bF+21Owvuf7gHbgUSAI7Ac+P9Tfu5Sw4E6g\nSmu9ArgG+AXWD+lftdbLQ/vcEKngenF+7D8HUErNAz4ZycD6oKvvfT7wY6316tB/UZmsQrqK//vA\nY1rrlcA3gZkRjK8nF/xutNYf1VqvBm4B6oAvRjLAHnT1vf8H1s3aFViJ4LoIxteTrmJ/BPhiKPZy\n4P4IxteT64Fg6Jr4DeC7wI+Br4f+HgO4aaiDkIQFT2NdXMD6PnzArVrrd5VSLqy7zvpIBdeLC2JX\nSmUA/wU8iPUjilZdfe+XAdcppd5SSv1WKZUcseh611X8y4BxSqk3sC5O6yMUW2/Oj93f6bXvAA9p\nrSuHPaq+6ep7bwMylVIG4AairoQS0lXseVrrLaFt7wErIxFYb7TWLwCfDT0dj3VTc5nW+u3QtnXA\nmqGOY8QnLK11i9a6WSnlxvpB/ZvW2lRK5WMVczOBvRENshtdxP5NrCqHLwHNEQ2uF11978A24P+G\nSihHsO6co1IX8X8D60Su1VpfiVXV89UIhtitbr57lFLZwAewqnmiUjex/xz4KVAMZANvRTDEbnXz\nmzmilFoR2uUGICliAfZCax1QSj2K9V0/zrk3xM3AkK+1M+ITFoBSahzW3fCftNZPAmitj2utpwAP\nA/8Tyfh60jl2rDrwycCvgCeA6UqpmIg99L3/TWu9K/Ty88C8iAXXB+fF/wRQA7wYevklYEGkYutN\nV7954Dbg8Whs/+msi9j/DFyhtZ4GPIZVVRWVuvjNfBL4mlLqTaASqI5kfL3RWt8NKKw2w/hOL7kZ\nhpqoEZ+wlFI5wOvAV7TWj4a2vaiUmhzapRmIug4XcGHsWuv3tdYzQ20RHwGKtdZfimyUXevqewde\nVUotDD3+ILA9ErH1RTfxv8vZ9pOVWCX0qNNN7GB95+siElQfdRN7IhBeq6gCq9NL1Okm9uuBO7XW\na7Bqc16LUHg9Ukp9TCn1tdDTNqxr4nalVLgKcy3wdpdvHkQjfmompdRPgdsB3WnzvwE/wKoLb8Hq\ndRR1dfrdxL5Wa+1RSo0H/hLFvQS7iv1fse6OfVgXns9oraOyarOL+E3gbqw7zySsu807tNYNEQmw\nB93Efi3WDcJSrXVjpGLrTTe/mx8B/w54sHqufVprfTwC4fWom9h/jNVu2I5VJf5gNJZwlVIJWFXF\nuYAT+G/gIFanERdWdeynhzr2EZ+whBBCxIYRXyUohBAiNkjCEkIIERNGZMJSSq1SStUrpfI6bfue\nUuoTkYxLCCFE90ZkwgppB/7Q6bk05gkhRBQbqQnLxBoLUaOU+nznF5RSX1ZKbVNKvaeU+l5o2/tK\nqYLQ49uUUj8Z/pCFEGJkG6kJKzxC+37gi0qpSaHnbqxup5eHuoMXKqWuA34HfDy0z92cnRhXCCHE\nMBmpCQsArXUt1px7f8T6LuKBLZ1mZn8HmAH8BbhNKTUaSNFaF0ciXiGEGMlGdMIC0Fq/jDWQ726s\ngYeLlVL20ESaK6xddCOwA/gJ1lx9QgghhtlITVgm53ayeBBrupFG4ClgE9baRmWhWYrBGtF9NfDX\nYYxTCCFEiMx0IYQQIiaM1BKWEEKIGCMJSwghRExwRDqA4aaUWgx8T2u9Wik1B/g11oqrJcB9Wmuv\nUupLwF1YnTB+Flq3Jvz+qcAWIFtrHa0rmwohxCVnRJWwlFJfweo8ERfa9Fvgi1rrK4By4H6l1Eys\nMVdLgNXAv4XWsUEplYK1HIBnuGMXQoiRbkQlLKAU+BBnBw7naa23hB6/h7Xo3jRgo9baq7Vux1qE\nb0mom/vDwNewehQKIYQYRiMqYWmtn8Oq/gs7opRaEXp8A9bKpfuAFUqpZKVUJrAUa0G+/wBe0Vrv\nDe1vIIQQYtiMqITVhXuAryml3gQqgWqt9UHg58CrwM+wxmNVA3cC9yqlNmCtuhmVS1kLIcSlasR1\nujjP9cCdWutapdRDwOtKqVFY0y8tV0qlYlUVbtZaF4bfpJQqA66KTMhCCDEyjdQSVni09CHgTaXU\n5tC2P2mtqwGllNqGVYr6ita6qZv3CyGEGCYy04UQQoiYMFJLWEIIIWKMJCwhhBAxQRKWEEKImCAJ\nSwghREyQhCWEECImSMISQggREyRhCSGEiAmSsIQQQsSE/w9EMLaWeML7EwAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Update:** A friend, that knows this data, challenged me to use the same\n", "technique on $v$. Here it is:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "df['pandas'] = rolling_median(df['v'], window=3, center=True).fillna(method='bfill').fillna(method='ffill')\n", "\n", "difference = np.abs(df['v'] - df['pandas'])\n", "outlier_2 = difference > 2\n", "outlier_3 = difference > 3\n", "\n", "fig, ax = plt.subplots(figsize=figsize)\n", "df['v'].plot()\n", "df['v'][outlier_2].plot(**kw)\n", "kw.update(color='g', marker='*', alpha=1)\n", "df['v'][outlier_3].plot(**kw)\n", "_ = ax.set_ylim(-50, 60)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAawAAADOCAYAAABigzZJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclNe9+PHPwLCvsoqgKCLHBffduJvFqEljlmZrmjRp\ne7Pc3mxtb9P29qb99famS7rdtkmaNs3SLDVpY6JJ1Bh3iVEQVECPoAKCgOybbMPM749nMKgwwIAO\nk/m+X6+84sw88zxfhuH5Puec73OOyWazIYQQQgx1Xq4OQAghhOgLSVhCCCHcgiQsIYQQbkESlhBC\nCLcgCUsIIYRbkIQlhBDCLZidfaNS6ingBsAH+AOwF3gZsALZwCNaa6mZF0IIMSicamEppZYC87XW\nC4ClQBLwLPB9rfViwAR8aZBiFEIIIZzuErwWOKKUWg9sAN4HZmqtd9lf/wi4ehDiE0IIIQDnuwSj\ngZHAGozW1QaMVlWnRiBsYKEJIYQQn3M2YVUCR7XWFuC4UqoFiO/yeghQO9DghBBCiE7OJqw9wKPA\nr5VSI4BA4BOl1BKt9U7geuATRzuw2Ww2k8nkaBMhhBCep8fEYHJ28lul1M+BZRjjYE8BBcCLgC+Q\nC3yjlypBW0VFg1PHFkII8cUUHR0y+AlrEEjCEkIIcQFHCUtuHBZCCOEWJGEJIYRwC5KwhBBCuAVJ\nWEIIIdyCJCwhhBBuQRKWEEIItyAJSwghhFuQhCWEEMItOL0elhCuYk7fz95iY2GAqxIWY5k1x8UR\nCSGuBElYwq2Y0/fjVVXF/1a9TVt7B1v8J+KzczuW6TOwhcoCAUJ8kUmXoHAre4t3sbrkx6S1HiXd\nepxrC37EnpoMzJkHXR2aEOIyk4Ql3MqiwEn8Kvr+84+f9LqDRQGTXBiREOJKkYQl3Io1IpL1jfv4\nSsdKlpYuZGtHBjZ/fyzTZ7g6NCHEZSZjWMKtWGbNYfyZLSwtjWd/WTnWkdW0L1nm6rCEEFeAtLCE\n21mz9Fu0efvSZvZjkrrP1eEIIa4QaWEJt2MLDePUpNkcai4mxT/I1eEIIa4QaWEJt2TpMBYetVis\nLo5ECHGlSMISbqkzUbV3uGzFbCHEFSYJS7gli9VIWJYOaWEJ4SkkYQm31NnCkoQlhOcYUNGFUioG\nyABWAFbgZfv/s4FHtNbSXyMui/NjWNIlKITHcLqFpZTyAV4AmgAT8Gvg+1rrxfbHXxqUCIXoRru9\nZSVFF0J4joF0Cf4SeA4otT+eobXeZf/3R8DVAwlMCEekS1AIz+NUwlJK3QdUaK232J8y2f/r1AjI\n1NnisuksumiXhCWEx3B2DOtrgE0pdTUwDXgFiO7yeghQ29tOoqNDnDy88Hgm4/rIy9tLvkdCeAin\nEpbWeknnv5VS24EHgV8qpZZorXcC1wOf9LafiooGZw4vBC0tFgDONbfL90iILxBHF6CDNTWTDXgS\neFEp5QvkAu8M0r6FuIQUXQjheQacsLTWXafKXjrQ/QnRF53FFlJ0IYTnkBuHhVvqvP9Kii6E8ByS\nsIRb+rysXW4cFsJTSMISbkm6BIXwPJKwhFtql4QlhMeRhCXcTofVis3eEyhVgkJ4DklYwu10HbeS\n9bCE8BySsITb6doNKF2CQngOSVjC7XTtBrRYrNhs0soSwhNIwhJup+u9VzbAKglLCI8gCUu4nYvv\nvbJYJGEJ4QkkYQm3c/G4lcx2IYRnkIQl3M7FCUsKL4TwDJKwhNu5uAtQ7sUSwjNIwhJu5+IuQItV\nxrCE8ASSsITbuaRLUFpYQngESVjC7XQmLJP9sRRdCOEZJGEJt9NZ1u7v521/LAlLCE8gCUu4nc4u\nQH9f8wWPhRBfbJKwhNvp7AL09/W2P5aiCyE8gSQs4XY6uwAD/IwWVod0CQrhEczOvEkp5QO8BCQC\nfsBPgaPAy4AVyAYe0VrLpa8YdOfHsM63sCRhCeEJnG1h3Q1UaK0XAyuBPwLPAt+3P2cCvjQ4IQpx\nIUvHRWNYkrCE8AjOJqy3gR912Uc7MENrvcv+3EfA1QOMTYhudRZZBPh2VglKQ14IT+BUl6DWuglA\nKRWCkbx+CPyqyyaNQNiAoxOiG+eLLuxjWO1SJSiER3AqYQEopUYC/wL+qLV+Uyn1iy4vhwC1ve0j\nOjrE2cMLD+Zj7wqMHBYIgH+Ar3yXhPAAzhZdxAJbgIe11tvtT2cqpZZorXcC1wOf9LafiooGZw4v\nPFxDYysAVksHALV15+S7JMQXhKOLT2dbWN/H6PL7kVKqcyzrUeD3SilfIBd4x8l9C+HQ+bJ2GcMS\nwqM4O4b1KEaCutjSAUUjRB9cMtOFVAkK4RHkxmHhdj4vurDfhyVFF0J4BElYwu10dgEG+HbOdCFd\ngkJ4AklYwu1YLm5hSZegEB5BEpZwO58XXcgYlhCeRBKWcDudRRe+PrIelhCeRBKWcDvtHTbM3l74\nmI2vrxRdCOEZJGEJt2PpsOJjNmH2NtkfS9GFEJ5AEpZwO5YOK2ZvL8zeXucfCyG++JyeS1AIV+lM\nWN5eJkxIwhKiv06fbeTDfYW0tnUQHODD2sVJDAvxc3VYvZKEJdyOpcOG2duEyWTCbPaShCVEP+w9\nUsprmzVtXcZ+q+pb+PYd0zCZTC6MrHeSsDCu0D/JKOaz3HKSRoSyIDWOpBGhrg5L9KDdYsXf1wcA\ns7cX7ZYrM4ZlTt/P3mJjyberEhZjmTXnihxXiMGy7WAxf99ynAA/bx65YTITEofx5w05HD5RxY6s\nMyybHu/qEB3y6IRls9nIzKtk3fZ8ztY0A1BQ1sC2gyV8fc0EFqTGuThCx44V1rDlwGmsNhthQb5c\nPy+R4RGBrg7rsrN0WPGxj1+ZvU1XpIVlTt+PV1UVz1Qbczp/EDAJn53bsUyfgS1Uln4TQ19uQTVv\nfJxHSKAPT31l5vlzxb0rx/Nff/mMddvymZIUSWSYv4sj7ZnHJqz8kjr+tfMEx4pq8TKZuHpmAqsX\njKagtJ7n38vhza15pI6JJDTI19WhXqK1rYO/fJBLhq644Pmcgmp+cM8st+iLHghLhw2zuTNhXZku\nwb3Fu3im+h32tuQCsLrkx3wv4lYWZkL7kmWX/fhCDER1fQvPrc/GZIJ/v3nyBRe2w0L8+PLyZF7+\n6Bib9hdx9zUpLozUMY+rEmxsbufX67L42WsZHCuqZcrYSH7ywBzuuiaFsCBfpiZHcfPiJJpaLLz5\nSZ6rw72EzWbj5U3HyNAVJMeH8V/3zuKPjy9m7aIxVNe38tu3D9HcanF1mJeNzWYzii68jL52nyuU\nsBYFTuLnw+47//jZ6AdYFDDpsh/3Smputch44BfUzqwzNLVY+PLyZMYlhF/y+oLU4QwL8WP34TM0\nNre7IMK+8bgW1usfHyf7ZDVqZDg3LRqDGjXskm1WzExgX245n+WWs3hKHBNGR7gg0u5ttY+1JceH\n8d27pp8v7V6zYDQ1jW3syCxhw94Cvrw82cWR9o3NZqOovJHwYF/CgntvGXZYjfGq8y0ssxfnrkCC\ntkZE8lbxGywtXQjA2yF7+eGoe7FMn3HZj305WTqs/GvXSbLyKimrPofZ28So2BCunpXAvInDXR2e\nGARWm4207DL8fLxZPGVEt9uYvb24ZtZI1m3PZ0dmCWsWjL6yQfaRRyWsDH2Wz3LLGTsilO/cOR0v\nr+4rYry8TNxzXQo/eTmdd/ecYnzisCFRPVNefY512/IJDfLloZtSzycrAJPJxJ0rxpGZV8H2rBJW\nzU8kOMDHhdE6ZrPZ2H24lM37iyitOgdAYmwItyxJIjUpssf3dc5qYe4yhnUlJr+1zJpD7KF4osuM\nVpX3mNoLugJPnKkjQ1eQd7qWUcNDWD0vkYjQoTsWAEbX8h/XHyH7ZDUBft6MHxVOc1sHhWUN/Pn9\nXGob2lg5d5SrwxQDdLyolqr6Fq5KHY6ffdHT7iyZNoINaafYmlHMdXNGnZ9JZijxmITV1NLOa5s1\nPmYv7l89ocdk1Wn08FCmJUeRlV9JbkENk8a4vpW1Jf00HVYbd109rttxKh+zFyvnjOIf2/LZllHM\njQvHuCDKvnl39yk2phXg7WVi9vgYGpvbOX66lt+9c5j7V09g/qTur+47u6w+L7rwOj+34OU2PO52\nCsw7AIgNv/X881sOnOatLt3HJ87UsyvrDHddkzJkq64sHVZ+8/Yhjp+uZXJSJA+vTcXPPjdjSWUT\nv/5HFuu25wNI0rrCOqxWvEymQbtITssuA2DBZMdFZAF+ZpZMjWfT/iLS9dke/wZdaeil0MskLbuM\n+nPtrFkwmrjIoD6950v2E/57e05hs7l2+p+mlnb2HiklMtSPmSq6x+2WTBtBkL+Zj9NP09I2NMey\nNn1WxMa0AmLCA3jm3+bz0E2pfOfO6Xznzun4+njz4oZc0o+d7fa9ndMwdS266LDasF6B309BkxeH\nEqeSO3Y6hyqMfv71u0/y1id5hAX78h+3TuGPjy/mgdUTCPI38/ctmqz8ysselzM2phVw/HQt08dF\n8a1bJp9PVgDxUUE8dfcMQoN8Wb/7JDUNrQM+XrulgwPHzrI1/TR7DpcO2liZzWaj3dIxKPvqq6aW\ndrJPVtFhHbwLJZvNxiubjvHIb3bxjV/s4JdvZtLWPvCfq7WtgwP6LJGh/qhRl45dXWzJNKPLcM/h\n0gEf+3LwnIR1pAwvk4klU7vvw+1O4vAQpo+LIr+kjpyC6ssYXe92ZZ2hrd3Kipkj8fbq+dfm72vm\nmlkjaWqxsDPrzBWMsG92ZpWwbns+w0L8+PYd0y4ooU0ZGc5/3jUdH7MXb23L6/ZE1Hmi+7zowvj/\nlVjEsaSyidAgXyYkDqO8ppmP9hXy/l4j8X7/KzOZlhxFgJ+ZqybH8ehtU/Hx9uKF93IoKm+47LH1\nR35JHRvSCogM9eOB1RMv6FruFBUewM2Lk2izWPnXzhMDOl67pYPfvn2Y59Zn88bWPF768Chvbx/Y\nPgFqG1t5+m8H+Ldf7eTBX+3gV29lUlh2+T7rhnNtvLLpGE/+YS+/XneIZ14/yNna5kHZ97aDJezM\nOoO/rzcjooI4VlTLa5v1gC+UcwqqaW3rYN6kWLz60GKLjQhkXEIYxwprqBykn20wDZmE1djczp/e\nPcJTf97Hq5s1J87UDdq+iysaKSxvYHJSRL/L1G+8yt7K2u26Vpalw8rWjGJj0HRq7/eGLZ+ZgJ+v\nN5v2Fw2pmcz3Hy3n1U2a4AAfvn3HNKLCAy7ZZlRsCFfPTKC6vpVPMkouef18wurSwur6/OXS0mah\nsq6F+KggxtsLdd7ecQJ/X2+euH0q0Rf9LGPiQvnGDRNpbe/gd+8cdrqVUtfYyt8+PMojv9nFQ7/e\nyb//ZhfffS6NZ508ObdbOvjLxlywwdfXTCTQv+dRgYWT40iIDmZvdhkFZfVOxW/psPLc+hyOFtYw\ndWwkD92UyoioID5OP82+3DKn9glQ09DKz9/I5PTZRpJGhBIzLIDcghp+8vIB/rIxl/Lqc07vuzv1\nTW384s1MdmadITTIl8lJkZwoqefpl/aTXzywc9WZyibWbc8nOMCHH351Fv993yxGDw9hb3YZOw8N\n7KLzWGENAKn9GNJYOCUOG7A32/nfz+UyqAlLKeWllHpeKZWmlNqulBrbl/cVljXw478dIF1XUFXX\nzI7MEp75+0GODlKrprMP96pe+nC709nKOnGmnpxTrmllHTxeQU1DKwsnxxHo33shRXCAD8umx1PX\n2Mbe7KHRtD99tpEXN+Ti7+fNk7dPc9gtu2p+IoF+Zj74tICmlgtLbC8puuhcYuQyJ6zOwpD4qCAm\nJH5eWXrvyvHEDOv+Zu2ZKobblo6lpqGV379zmNa2/nXxZOVX8r0/72P34VIC/bwZPiyQiFA/Oqw2\ncgpq+H+vpPP2jvx+dYdu3n+aszXNrJiZ0G2FbFdeXibuWGFUm76/p6BfsXfamFZAVn4lk0YP4+G1\nk5k9PoZH1qbi7+vNyx8d42xN/xNLu8XKb9Ydorz6HNfPG8UP7pnJTx6Yy5N3TCM+Opi07DK+/+I+\n/rXrpFMxX6yxuZ1fvpVJSUUTK2Yk8My/zeex26bwwOoJtLVb+f0/D1M2gAT5yqZjtFus3LtSMSzE\nDx+zNw+vTSXI38w/d5wY0EXnsaJafMxeJI3o+83ts8fH4Ofjzd4jpU53tbe1d/T7+94Xg93Cugnw\n1VovAL4HPNvbG3ILqnnm9YNU17fwpYVj+NMTS3hkbSomE/zh3SMUVzQOKKAOq5VPc8oI9DMzNTnK\nqX10jmWtd9FY1scHTmMCrp6V0Of3XDt7JGZvLz7aV3hBX7s5fT++Wz7Cd8tHmNP3X4ZoL2XpsPLX\nD3LpsNr45g2TSBwe4nD7IH8fVi9IpKnFwqbPii7al/H5dxZddP7/chdedH4PR0QHMTI2mNSkCFbN\nS2TuxFiH71s5dxSLp8ZRWN7Ab9Zlca6lb+OKFbXNvLghF5vVxlevUzzz4Hz++2uz+ckDc3n2kat4\n8o5pRIb58dG+It7pY/dadX0LGz8tICTQh5sW9a0gZ+LoCMbEhXAov7LfXUQVtc18uK+I8GBfHl47\n+XzVWVxkEPdcp2hrt/LWJ/n92ifAhrRTFFc0smhKHLcuGXu+OGHS6Aievn82D9+USlSYPxvTCtiR\ndWkrvb/Wbc83ktXMBO66ZhxeXkZBxFWT4/jqSmXc2/mPLKfuXyosayCvuI7UpAhmqpjzz0eFBbBo\nygiaWixk5lU42EPPGpvbKa5oJDk+rF8Vf/6+ZmaPj6GyrgVdVNuvY548U8+LG3J49Pd7eOrPnzp1\nQeLIYCesq4BNAFrrz4BZPW1o6bCyI7OE3759iA6rlYfXpvKlhWMwe3sxU8Vw/+oJNLd28Lu3D9Fw\nrs3pgI6cqKausY05E2OdLtMcFRvCzJRoTp6pv2R2icvtREkdJ87UMzU5ith+TLsUHuzHoilxVNS2\nsPeI0cLsnF4IG2ADr6oqfHZux1Q/eN2v3fnw00KKyhtZOCWuzxcNK2YkEBbsy9aM4gt+/593CRon\nKe/za2Jd3oR1prIJgISoYLxMJp748jRuXdp7B4LJZOIr1ypmjY/heHEdv3jDuDhzxNJh5fn3cmhu\ntXD3tSksnR5/ybjlpNER/Ne9s4mLDGTT/iK2pp/uNZZ3dpygrd3KrUvG9qml3mn5jARswPZ+nvzX\nbcvH0mHly8uSCfC7sOtx3sRY1MhwsvIrOXKyqs/7PFVaz4efFhEV5s8dK8ZdUknnZTIxa3wM375j\nOsEBPry+5Ti6qKZfcXd18kw9ew6XkhAdzB0rki853uKpI1izIJHKuhb+6cRY3/bMYsD4vl9s4RSj\nR2i3kwUQnT93X4otejr2nsN965K0dFh5d9dJ/ue1dD7NKSfQ30xtYxu/eitrUIp2Og12wgoFunZ2\ndyiluj3GQz//hFc3a7y9vXjstqkXXF0AzJs4nJsWjaGqvpUX3s/BanWuZfNJhvGHvHyA5cW3LB2L\n2dvEm5/k9an67lxLOzuySvjFGwf5zp/28vRL+3lnx4l+t9A+tp+Irpk9st8xr56fiJ+PN+/sOEFj\nczte1caJoaahlT1HztBu6cDU0oI582C/991XVXUtbEgrYFiIH3csH9fn9/n6eLNqXiKtbR1s2v95\nK+vzoouLWliXueiipMJIWCOi+j9Xo9nbiwdvnMTSaSMoOtvIf/31M7ZnlnRbZWaz2Xh1s+ZUaT3z\nJ8Wy0EE3dnCAD4/fNpXQIF/e3JrHbgfjHaVVTXyWW86o2GCumtK/rvE5E2IIDvBh96HSPlfkHS2s\nIeN4BckJYd22Qk0mE3ddk4LJBG9s7b7A5mIdVit/+/AoVpuNr62acEkS7Co6PIBH1qYC8Mom7dT5\nw2qz8cbW4wDcfc24HoudbrxqDPFRQezKOsOJkr5f/J1raWdfTjlRYf5M7ubewxFRQYyNDyX3VDVV\ndY4vcrpzzN46Gt9L1293xiWEETMsgAxd0WuvQFF5Az97LYMNaQVEhPjz5B3T+NXDC7hp4Rgq61r4\n3TuHnD5/X2yw78OqB7r293hprbu99K2sbWHVgtHcujyF6GGXDr4DfO3GyZRUneNAbjmb0ou5d/XE\nfgVzuryBnIIaUsdGMn3SwCayjY4O4Zbl4/jHx8fZklHCAzem9njM93efZHvG6fN9uFFh/pTVNFO0\nr5CYyCBu6eOJu6yqiXRdwei4UBbNHNnv+zKio0O4e+V4XtqQw8Z9RXwrLJAOSwf/2JZPWXUTwUF+\nLJmeAAEBEO24m85Z/9x9ig6rjfvWTCRxZP/+cG65WrF5fxHbD5Zw9/UTCQv2I6jS6GIICwsgOjqE\nkBCjyjAk1Hh8OdQ3tXG6opHIMH8SRzp/P94TX5lF6rgi/rYhm9c2azbsLWDh1BFMTo4iKT4MP19v\n/rktnz2HS0lOCOPxu2c5PCmD8Tv+6YML+MFze3l50zGCQ/y5dm7iJdu9uS0fG3DXygnExvR/JYKV\n80fzzrY8jhbXs2J27/dlPbvuEAAP3zqVmB6OFx0dwuqrxrBxzyle35rPt78y0+F3fOOekxRXNHHN\nnFEsntV7DNHRIRzMr+Lj/UWcKG9kQQ+zPPRk76EznDxTz6Jp8Syc6fh437p9Ot/74x7e/CSfXz+2\nGO9uKi8v9v6uE7RZrKxemERsbPef0aqrkvi/dVlknqzmzmtVv+LPL6nD18ebOVNG4GPu+Ybhnlw3\nbzSvfXSUYyV1XDdv9CWv22w2/rk9n9c+OorVamPZzAS+uXbK+QkL7r9pMjVN7ezMLCa3uI5lM/t/\n0X2xwU5Ye4EbgLeVUvOAwz1t+Nb/rKK+9hxYLFRU9FztdO+1KRSW1vPOtjyC/bxZ3I+y9Le3agAW\nT45zeIy+WjYljm0Hinhv1wlamo17ugL9zZxraedUaQN7jpSyP7ccGxAZ6s+a+YnMnzSciFB/6pra\n+PHf9vPKh7lEh/j2abqnVz40vgjXzk6gstK5sbx546PZsi+ILZ8VMmFkM6aqSsqqjdbCp4dLSUmO\nxWfidGwVDZRWNVFY3kBUaADDIwMHPFNGbWMrm/cVEhXmz4SEMKd+ByvnjOKNrXk8+/d0Hrl5MpVV\nRuytLe1UVDTQbm/tnq1sIMR38Ite65uMbo26xjaunzdqwN+jGWMjSPr6XD5IK+Szo+Vs3HuKjXtP\nXbDNiKggvnXzZBrrm+nLbz3Yx4snb5/Gr97K4g9vZ+FtszJl7Oddr9X1LWxLP01sRCDJscFO/Qxz\nVBT/3J7H+h35TBnt+MJDF9WQc7KKKWMjCfc3OzzeDfNGcaygml1ZJYQH+XDToqRut2s418ZrHx4l\nwM/M6rl9/z0snRrH1v1FvLXlGMnDg/t10bd+h3Ez+MrZCb0eLybElwWpw0nLLmNL2ilmjY9xuL3V\nauO9XScwe5uYPjaix/2Pjw/Fz8ebLfsKWD4trk+l6QD159ooLGtgQuIwap0cR5o6Zhh/N8FHe08x\nY+yFLUCrzcZbn+SxNb2YiFA/7ls5ntSkSJobW2hu/Lw1uHruSPYcKuG1D3MZHx/a7S0UF3N04TnY\nCetd4Bql1F7746/1tGHXGxUdCfT34bHbpvKz1zJ4dZNmWIhft83nizU2t5N2pIyIUD+mpzhXbHEx\nXx9vvnnjJJ5fn8Om/UVsPlCEt5fpgu6oUTHBrFkwmhkp0RfMphEWZAw8//z1g/x5Qy4/++Y8h1fP\nZ2vOkXakjLjIQOaMdzyw74jZ24uvrZrAr/+Rxe9OBzC1sJpwPzOpSZHsyatlnSmRlDOt7Pooi+yL\nqiCjwvxJTggjdUwEk5MiCQns3y0Bmz4rwtJhZdX8xD59UbuzdHo8mXmVZOZVsn73SeKjgoHP77/q\n3G9bWwcllU2crT5H3bk2Av3M+Pua8TIZM4BEhwcQHuLX5z94MBLuL9/MpLTqHMtmxHPLkj4VvfYq\nPNiPu69N4fYVyeQV13HyTB0lFU10WG0E+Zu54aox/f6sR8WG8PiXp/LM6wd54f0cvn/PLOKjjErM\nLQeMGVKunzuq1xleehIVFsC05Cgy8yo5eabe4XpxGz8tBOjTfHQ+Zm/+/ebJ/M+r6cY9bcMCul3W\n5587T3Ku1cKdK8b169aUuMggZqREk3G8gmOFNX2eF7SksoljRbVMSBzW54kGVs1LJC27jO2ZJb0m\nrMy8Ss7WNLNwShyhDn7XAX5mZk+IYc/hUo4V1jCxj/EfLTDGr7pWtPZXRKg/k8ZEkH2ymhMldYyN\nNyoNbTYbr23W7Mw6Q3xUEE/cPq3HFSKiwgNYMm0E2w6WsOdIKUunDWxoZlATltbaBjw0mPsEGB4R\nyH/cMoVfvpXJ8+9l86N7Z/dagLBuWz6t7R2sXZzk8Ebb/ho7IoyffXMuH6cXk5VXic1mw9/Xm9Fx\noahR4UwaHdHjVVxyfBg3LBjN+j2n2JBWwJeX9TxB7ca0Qqw2GzdeNcbpk0ynMXGh/Pj+Ofz1g6Pk\ntybzxFgLSSNCedcWz6H0YramGwO/KQlhzEiJpraxjeKKRk6V1rMvp5x9OeX4+XrzwKoJvf4hdqpp\naGVHVgnDQvy4agDripm9vXjoplR++ko6G9MKiQz1O/+88X/js3n2H4d6LcENC/Ll9hXJzJ0Q2+uV\ndnV9C798M5PymmaunT2S25dfOuA+UGZvLyYkDhvQSaWrMXGh3L9qAi+8n8Pv3j7Et++cTkVNMx+n\nn2ZYiN+Ap9pZPiOBzLxKPskoJmlE993zuqiGnFPVTEgcRnJ830qpQwN9eey2qfz01Qz+9uEx+6wM\nn38mx0/XsuuQcXJcNqP/J7xV8xPJOF7Bx+nFfU5YOw4aBSbL+3G8EVFBjB8VztHCGs5UNjEiqudE\nt9k+LnvdnN67NhdNiWPP4VJ2Hy7tc8I6fMIYr+7Lxb0jq+clkn2ympc+PMrTX5uNj9mbDWkF7Mw6\nw6jY4PPFLY6sWTCaPYdL+SCtgEVT4gZ0PnabuQSTE8K4b+V4XtyYyx/ePcIP75nV40SOOQXV7DlS\nyqjYYFbMHPy53HzMRkHAqnmXjhX0ZuXcUew5UsrHB06zaEpct1dvJRWNpGWXMSIqiNl9TBC9iQg1\nBkMbmoz/zOBTAAAVOElEQVSxIBuwNr6C9/acYvLYCOZPGk5CdPAF77HZbJRUNJGVX8kHnxbyp/XZ\nrJqXyC1Lkno9ef9rl1GRdueK0QOeRDM4wIdHb5vC3z46xgn7TZqdM7uPHm50mUSHBzA6LoS4iEDC\ngn1pbu04XxzT0tbB2ZpmDuVX8uf3c9mXU84Dqyf02IqprG3mF29mUlnXwur5idy8uPefd6iYOzGW\nitpm/rXrJP/7WgbtFiveXiYeuil1wL+HCaOHMTwikAPHyrl9RfIlLYOqOmPNJS+TibWLu+/a60lc\nZBCPrE3lN+sO8Yd/HeGJ26cxJi6UdksHr2w6hgm47/rxTrXUx8SFMjImmOxTVZxrae+1QrKlzcLe\n7FLCg32ZNq5/vTPLZyRwrKiWHZkl3NXDulL5xXXkl9QxZWzk+VawI8nxYQyPCCRDV9DU0k5QL/Fb\nbTayT1URFuTLqNhgh9v2Ro0axoqZCXySUczLH2n8fLzYkXWGyFB/Hr9tap+GDcKD/Vg4JY5tB0vI\n0BXMmeB8j5H3008/7fSbB+jpc/0sVx8ZE0zDuTYOn6iiqq6FGSnRl5xImlst/O4dY02oR2+dQkTI\n0Jox29vbi4hQfz47Wk5p1TnmTxp+wc9gs9l44f1cKmqbuX/VhD53R/SFyWTC3/fza5QR9ivWSaO7\nnwHEZDIRGuRLyshwpo+LIrewhqy8SgL9zOe7B7pTWNbA37ccJyE6iHtXjh+Uk31IoC+Lpoxg+Yx4\nZqZEo0aFYzKZiI0IZPX8RJbPSGD6uGjGJYQzMiaEpBGhpIwMJ2VkOBNHRzB7fAxzJsZSUtFEzqlq\nDhwtZ1xC+CVdGeU15/jFm5lU1bdy08Ix3LRojNskq04pI8MJ9DeTrito77DytVUTmObkPYhdmUwm\nrDYbh09U0W6xMqXLuEZjczu/XXeIs7Ut3HVNCrNU/y+0osMDiArzZ//Rs+zLLcfX7M2H+wo5UVLP\nipkJLBlAd1Jjczu5BTXERQYxKtZxcc7m/UUcOVnNyjmjmJDYvyKbmGEB7Dp8hoKyepZMi8f3oqGP\nmoZW/vJBLnWNbdx3/XiiwrovOOvKZDLRZukg+2Q1ESH+DrtjwVg1/eP0YmaPj2FGSs/zjvaVGhnO\ngaNn0adrKShrICzIlyfvmHbJ7C6OxA4LZFtGMTUNrb3WIQQF+f24p9fcKmGBcSNjbmE1R05WExzg\nc8Ed3DabjRc3HiW/uI7r545yamaLKyEuMpDCsgayT1XT2t5B6pjP//APHq/gw31FTBkb2eMAtCuE\nBvkyMyWaz46Wc/B4BUkjQontZoYHS4eVF97Poaq+lW/cOKnbbQbCz8ebiFD/C5JIXxNKkL8P81OH\n4+1lIjOvkj1HSunosJEcH4q3lxcZuoLfv3OY+qZ2bl06ljULRrtdsuo0dkQYSSNCmTI2stsxIWfF\nRwWRoSs4crKK4AAf4qODOHDsLP/3zmHO1rawdNqIASX5kTEhJNj3efhEFeU1zSREB/HNGyYNqIUY\nFuTLJxnFdFhtzHPQNVrX1Maf1mcT4Gvmmzf2/5heXiZsNhuH8qvIK6kzup+9TJRUNHHweAV/fj+H\n8ppmFqQO55pZfa/8jQkPYGt6MWdrm1k2Pd7h+3YfOsOxolpWLxjdpxZcbzq7rkODfLnxqjHcsWJc\nn9au6yo4wIfCsgaOFtaQOibC4dI7X6iE5eVlYnJSJPtyysjMq0SNDD9/lbJ5/2k+Tj9NSkIYD6yZ\nOOCxn8vFZDIxZWwUmXkVZOVX4e/rzZi4EArLGvjLxlzaLVb+49YpQ249qwA/MykJ4aRll5FxvAI1\nMpzILl88q9XGixtzOXKympkp0ayeP9p1wfbAZDKhRhnjK0eLajiUX8XW9GJ2Zp1hZ9YZTMA91ymu\nnjXwElxXix0WeEk370CZvb2YMjaS/UfPknHsLB/uKyJDV2Czwc1Lkli7OAmvAY4Zj4gKYkJiBMMj\nArl16Vi+tCgJXyfKsrsKDvAh83gFJ87Uc/XMhB7LvP+xLY8TJfXcvjyZlJH9v+EWYGx8GGXV5zhy\nspqs/ErW7z7F5gOnOWRvmd6+Yhy3Lh3br6Tu72umur6FnFM1RIcHOGwlvr3jBHWNbdy7UjlVzt6d\n0CBfxo8aRnR4gNPn1dBAH9Kyy2ht73A4Fv6FSlhgnDhHDw8lLbuMT3PKaGu3suvQGbYcOE1YkC/f\nvnN6v+7kdwUfsxepYyLOJ95Ps8v4+MBpmls7+PKy5EHpwrkchoX4MSIqiP1Hy/ns6FlGDw8hOjyA\n2sY2Xt18jP1HzzIuIYxHbp7sdGXglRAzzJj6xtJhpbG5ndb2DkbGBPPorVMdLiApjJaqshcXDI8M\nZNZ4Y2aaacnR/arCdCQi1J9xCeGEB/sNWiu3qaWdHAfdguU153jpw6PERQZy36rxTv8sJpOJqclR\n5J2u5VRpA4H+ZmamRLNsejy3L09manKUUz/TqNgQth0sprCsgWUz4rtNHPVNbbz1SR7jEsJY1s3s\nGa4UFeZPxvEK8k7XsWx6fI+V4o4SlsmF6zzZBnpPS/bJKl7edIzqemPqjzFxodx3/XhGxgzuVeXl\nVFVnzO+253ApIYE+fH3NxD5XArlShq7g+fey6bDaCAvypamlHUuHjcThIXznjukOZwEXwhXKa87x\n1Av7SB0TwRO3T7vk9T2HS3npw6Pcc23KoJzs2y1WqupbiB0WMGhJ942tx9maXszd16SwYualMb67\n6yQb0gq48+pxXDMEewm27C/irW35DuOLjg7p8cNy67NKalIk/++BuXycfpqYYQHMmdC3NV+Gksgw\nf+5dOZ61i5Lw8/F2uIT1UDJTRfOdO6ezM6uEnFPVRIT6s2peIgtShw/plpXwXLHDAhkTF0pOQTV1\nja2XjMN0TnA8speijL7yMXsxvB/zf/bFmvmj2X24lH/uPEFqUsQFY8TnWixszSgmOMCHxf2c1eNK\nmZc6nLd3nGDP4VKnEqpbJywwugc716xyZ/1dp2so6KzCE8JdzJ8Uy6nSej47epZrL5qfszNhDUah\nwuUSGuTLV69VvLgxl+fX5/D9e2aeLwzZdrCY5lYLtyxJGrIXvqGBvkwZG0lmXiWFZQ2XrNxQUFbv\ncKYLuRQWQniMORONXphPcy5dnLD4bCNRYf69zt/oavNTh7NwsrFkzYsbcjjXYqGovIEtB04T5G9m\n+RAbu7pY50zwuy6aCX5fThk/eTnd4XuH9m9GCCEGUWigL6lJERw+UXXBbBR1TW3Un2tnWnLfFzp0\npbuvSaG02pgcO6+kjoamdqw2G3euGDfkE+7kpEgiQv3YmXmG+ROHk5wQxtHCGv76wVEC/By3DN2y\nSlAIIZxlMpnI0BWkZZex50gZESF+nGu18Gl2GbPGxwzaVFmXk9nbiwWpw7FhFJ9Fhvrz0E2pDu8x\nGyq8vEwkxoawN7uU3IJqzrVYeP3j49hs8OhtU0kcEfbFrBIUQoj+ard08OKGXEqrznGmymhlXZUa\nx7rt+Tx0U+qgTYd2pVTXtxAS6DNo91xdKe/vOcX6PcZKBaFBvtx7nWJ6SvQXt0pQCCH6y8fszcNr\nJwPw3PpsDhw7y54jxqq+CdFDt+CiJ45mjRjK1iwYTVOLhQA/b66bM6pPXZmSsIQQHmvptBHUb9tN\n1dkjxJhMxJ8Mwho519VheQQvLxN3Xt33VchBqgSFEB4stSKfUd4t7Ijdxe743Zirq/HZuR1Tfd+X\nuhdXjiQsIYTHSivZxUvJr1EQUkSefwGrS37MnpoMzJkHXR2a6IYkLCGEx1oUOInfx33j/ONnox9g\nUcAkF0YkHJEqQSGExzKn7+fnec+ff2wC/nPEPVimz8AW6h73ZH3RSJWgEEJ0wzJrDuPPbGGt30wA\n3m3NoH3JMhdHJXoiLSwhhEcz1dedH7OSlpXrOWph9TthKaXCgL8DIYAv8ITWep9Sah7wW8ACbNFa\n/6SXXUnCEkIIcQFHCcuZoovHgY+11kuB+4A/2p9/HrhTa70QmKuUunTBGSGEEMJJzoxh/QZotf/b\nB2hWSoUAvlrrU/bnNwNXA1kDD1EIIYToJWEppR4AHrvo6fu01hlKqeHAa8CjQBhQ32WbBiBpMAMV\nQgjh2RwmLK31X4G/Xvy8Umoy8CbwpNZ6t1IqFGNMq1MoUNvbwR0t1CWEEEJ01e8uQaXUROBt4Dat\n9REArXW9UqpNKZUEnAKuBZ7ubV9SdCGEEKIrRw0ZZ8awfoZRHfh7pRRArdZ6LfAg8DrgDWzWWh9w\nYt9CCCFEt+Q+LCGEEEPGYJe1CyGEEFecJCwhhBBuQRKWEEIItyAJSwghhFuQhCWEEMItSMISQgjh\nFiRhCSGEcAuSsIQQQrgFSVhCCCHcgiQsIYQQbkESlhBCCLcgCUsIIYRbkIQlhBDCLUjCEkII4RYk\nYQkhhHALkrCEEEK4BUlYQggh3IIkLCGEEG5BEpYQQgi3IAlLCCGEWzA7+0al1HhgHxCjtW5TSs0D\nfgtYgC1a658MUoxCCCGEcy0spVQo8CzQ0uXp54A7tdYLgblKqWmDEJ8QQggBOJGwlFIm4AXgKaDZ\n/lwo4Ke1PmXfbDNw9WAFKYQQQjjsElRKPQA8dtHThcBbWuvDSikAExAK1HfZpgFIGsQ4hRBCeDiH\nCUtr/Vfgr12fU0rlAQ/Yk9lwjNbUDUBIl81CgdreDh4dHdLbJkIIIQQAJpvN5vSblVKnAGUvusgE\nbgFOARuBp7XWBxy83VZR0eD0sYUQQnzxREeHmHp6zekqQbuu2e5B4HXAG9jcS7ISQggh+mVALawB\nkhaWEEKICzhqYcmNw0IIIdyCJCwhhBBuQRKWEEIItyAJSwghhFuQhCWEEMItSMISQgjhFiRhCSGE\ncAuSsIQQQrgFSVhCCCHcgiQsIYQQbsGVUzMJIYQQfSYtLCGEEG5BEpYQQgi3IAlLCCGEWxjoelhu\nTynlA7wEJAJ+wE+BE8Cf7ZvkAV/XWne4JsKedRe71nqD/bW7gH/XWi9wYYg96uFzL8ZY/PO4fbPn\ntNbrXBOhYz3E/xnwIhAOmICvaq0LXBVjT3qI/S6MFcQBxgBpWuu7XBNhz3qIPQ/4C8b6fMcx/l6H\n3OB8D7EXAc8DFoyf40GtdZvLguyBUsob47udgvE5Pwi0Ai8DViAbeORyf+7SwoK7gQqt9WJgJfBH\njC/S97TWC+3b3OCq4Hpxcex/AFBKTQfud2VgfdDd5z4DeFZrvcz+35BMVnbdxf9z4DWt9RLgR0Cq\nC+Nz5JLvjdb6Tq31MmAtUAM87soAHejuc/9vjIu1RRiJYLUL43Oku9hfBB63x14CPOzC+BxZA1jt\n58QfAj8DngW+b/95TMCXLncQkrDgbYyTCxifRztwi9Z6j1LKF+Oqs9ZVwfXiktiVUhHA/wCPYXyJ\nhqruPveZwGql1E6l1F+UUsEui6533cV/FTBSKfUxxslpm4ti683FsVu6vPYT4Pda6/IrHlXfdPe5\nNwORSikTEAIMuRaKXXexJ2it99mfSwOWuCKw3mit3wP+zf5wNMZFzUyt9S77cx8BV1/uODw+YWmt\nm7TWjUqpEIwv1A+01jal1CiMZm4kcNilQfagm9h/hNHl8ATQ6NLgetHd5w7sB75tb6GcxLhyHpK6\nif+HGH/I1VrrazC6ev7ThSH2qIfPHqVUDLAco5tnSOoh9j8AvwNygRhgpwtD7FEP35mTSqnF9k1u\nAIJcFmAvtNYdSqmXMT7r17nwgrgRCLvcMXh8wgJQSo3EuBp+VWv9FoDWukhrnQK8APzalfE50jV2\njD7wZOA54E1golLKLWK3f+7vaq0z7S+vB6a7LLg+uCj+N4Eq4H37yxuAWa6KrTfdfeeBW4HXh+L4\nT1fdxP53YJHWegLwGkZX1ZDUzXfmfuAppdRWoByodGV8vdFa3wcojDFD/y4vhXAFeqI8PmEppWKB\nLcB3tdYv2597XymVbN+kERhyBRdwaexa6wNa61T7WMQdQK7W+gnXRtm97j53YJNSarb93yuAdFfE\n1hc9xL+Hz8dPlmC00IecHmIH4zP/yCVB9VEPsQcCDfZ/l2IUvQw5PcS+Brhba301Rm/OZheF55BS\n6h6l1FP2h80Y58R0pVRnF+b1wK5u3zyIPH6mC6XU74DbAN3l6R8Av8DoC2/CqDoacn36PcR+vda6\nRSk1GnhjCFcJdhf79zCujtsxTjzf1FoPya7NbuK3AfdhXHkGYVxt3qW1rnNJgA70EPsqjAuEBVrr\nelfF1psevje/Av4LaMGoXPuG1rrIBeE51EPsz2KMG7ZidIk/NhRbuEqpAIyu4uGAD/C/wDGMohFf\njO7Yb1zu2D0+YQkhhHAPHt8lKIQQwj1IwhJCCOEWJGEJIYRwCx6ZsJRSS5VStUqphC7PPaOUuteV\ncQkhhOiZRyYsu1bgb10eS/WJEEIMYZ6asGwYN+9VKaUe6fqCUupJpdR+pVSaUuoZ+3MHlFKJ9n/f\nqpT67ZUPWQghPJunJqzOKUUeBh5XSo21Pw7BuE9ivv3+pXFKqdXAX4Gv2re5j89nchdCCHGFeGrC\nAkBrXY0xSewrGJ+FP7Cvy1Iiu4FJwBvArUqpOCBUa53riniFEMKTeXTCAtBab8S48/w+jDvl5yql\nvO0zPy82NtH1QAbwW4zJZYUQQlxhnpqwbFxYZPEYxvxY9cA6YC/GYnyn7NPqgzEFyXXAP65gnEII\nIexkaiYhhBBuwVNbWEIIIdyMJCwhhBBuwezqAK40pdRc4Bmt9TKl1FTgeYwlwvOAB7XWbUqpJ4Cv\nYBRh/J99obXO948H9gExWuuhuhS3EEJ84XhUC0sp9V2M4gk/+1N/AR7XWi8CSoCHlVKpGPdczQOW\nAT+wL7yGUioUY/2alisduxBCeDqPSlhAPnAzn984nKC13mf/dxrGKrETgB1a6zatdSvGqrHz7GXu\nLwBPYVQUCiGEuII8KmFprf+F0f3X6aRSarH93zdgLLV9BFislApWSkUCCzBWkP1v4AOt9WH79iaE\nEEJcMR6VsLrxNeAppdRWoByo1FofA/4AbAL+D+N+rErgbuABpdR2jGWiN7smZCGE8EweV3RxkTXA\n3VrraqXU74EtSqkojOmXFiqlwjC6Cj/VWo/rfJNS6hRwrWtCFkIIz+SpLazOu6WPA1uVUp/an3tV\na10JKKXUfoxW1He11g09vF8IIcQVIjNdCCGEcAue2sISQgjhZiRhCSGEcAuSsIQQQrgFSVhCCCHc\ngiQsIYQQbkESlhBCCLcgCUsIIYRbkIQlhBDCLfx/H5GFWQ4aOa8AAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that I had to reduce the threshold from 3 -> 2 to get them all." ] }, { "cell_type": "code", "collapsed": false, "input": [ "HTML(html)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n", "

This post was written as an IPython notebook.\n", " It is available for download\n", " or as a static html.

\n", "

\n", "
python4oceanographers by Filipe Fernandes is\n", "licensed under a Creative Commons\n", "Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.\n" ], "metadata": {}, "output_type": "pyout", "prompt_number": 9, "text": [ "" ] } ], "prompt_number": 9 } ], "metadata": {} } ] }