profiling obspy read_stream

read_stream_profile.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Profiling obspy read of mseed files\n",
    "\n",
    "I was looking for ways to speed up some geophysical code that uses obspy. I found there are some large overheads in reading time series data that didn't make a lot of sense in using obspy.read(), and obspy.read(format='mseed') compared to calling the _read_mseed function directly from the obspy.io.mseed.core module, so I decided to investigate. \n",
    "\n",
    "My setup:\n",
    "- ubuntu 16.04\n",
    "- python 3.6.2\n",
    "- obspy 1.0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import itertools\n",
    "import os\n",
    "from functools import partial\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from obspy import read\n",
    "from obspy.io.mseed.core import _read_mseed as mread\n",
    "\n",
    "PATHS = {'small': 'small.mseed', 'normal': 'temp.mseed', 'large': 'large.mseed'}\n",
    "READ_FUNCS = {'read': read, 'read_format': partial(read, format='mseed'), \n",
    "             'mread': mread}\n",
    "\n",
    "# make streams of various sizes\n",
    "\n",
    "def make_small():\n",
    "    st = read()\n",
    "    for tr in st:\n",
    "        tr.data = tr.data[:100]\n",
    "    return st\n",
    "\n",
    "\n",
    "def make_large():\n",
    "    st = read()\n",
    "    for tr in st:\n",
    "        tr.data = np.tile(tr.data, 25)\n",
    "    return st\n",
    "\n",
    "\n",
    "MAKE_FUNCS = dict(normal=read, small=make_small, large=make_large)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# make streams and save mseed files\n",
    "for key in PATHS:\n",
    "    st = MAKE_FUNCS[key]()\n",
    "    st.write(PATHS[key], 'mseed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 3: 4.04 ms per loop\n",
      "100 loops, best of 3: 2.53 ms per loop\n",
      "1000 loops, best of 3: 486 µs per loop\n",
      "100 loops, best of 3: 4.19 ms per loop\n",
      "100 loops, best of 3: 2.62 ms per loop\n",
      "1000 loops, best of 3: 571 µs per loop\n",
      "100 loops, best of 3: 6.38 ms per loop\n",
      "100 loops, best of 3: 4.65 ms per loop\n",
      "100 loops, best of 3: 2.62 ms per loop\n"
     ]
    }
   ],
   "source": [
    "# profile each stream type with each function\n",
    "df = pd.DataFrame(columns=list(PATHS), index=list(READ_FUNCS))\n",
    "for stream_type, path in PATHS.items():\n",
    "    for func_name, func in READ_FUNCS.items():\n",
    "        time = %timeit -o func(path)\n",
    "        df.loc[func_name, stream_type] = time.best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>small</th>\n",
       "      <th>normal</th>\n",
       "      <th>large</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>read</th>\n",
       "      <td>0.00404061</td>\n",
       "      <td>0.00419304</td>\n",
       "      <td>0.0063838</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>read_format</th>\n",
       "      <td>0.00253432</td>\n",
       "      <td>0.00262287</td>\n",
       "      <td>0.00464931</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mread</th>\n",
       "      <td>0.000485681</td>\n",
       "      <td>0.000570669</td>\n",
       "      <td>0.00261966</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   small       normal       large\n",
       "read          0.00404061   0.00419304   0.0063838\n",
       "read_format   0.00253432   0.00262287  0.00464931\n",
       "mread        0.000485681  0.000570669  0.00261966"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a dataframe of results\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7efc3f2f64a8>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAD9CAYAAAB9YErCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVfed//HXl10WkV1lR9ncFxSyGBdcALfYCGnUNmZp\naqfNtDPNTKZpp01/07TJTNqJbTpN0zRNE80C2mwKGpcYExPciDuCC7vKooAgO3x/f5wjEoNIiHC5\n8Hk+Hj7iPfece783D+R9v+d8z+ejtNYIIYQQPWFj6QEIIYSwXhIiQgghekxCRAghRI9JiAghhOgx\nCREhhBA9JiEihBCixyREhBBC9JiEiBBCiB6TEBFCCNFjdpYewFfh7e2tQ0JCLD0MIYSwKgcPHqzQ\nWvv0xmtbVYiEhIRw4MABSw9DCCGsilKqoLdeW05nCSGE6DEJESGEED0mISKEEKLHrOqaiBBCADQ3\nN1NcXExDQ4Olh9KvODk5ERAQgL29fZ+9p4SIEMLqFBcX4+bmRkhICEopSw+nX9Bac/HiRYqLiwkN\nDe2z95XTWUIIq9PQ0ICXl5cESAdKKby8vPp8diYhIoSwShIgX2aJ/ydWFSJXmq/Q2tZq6WEIIYQw\nWVWI5F/OJz4tnl/v/TWHyg4h/eGFEAOFq6srAPn5+YwbN87Co+k+q7qwHugWyBS/KWzM3cgbJ99g\npMtIEkITSApNIsIjQqa3QgjRx6xqJjLUYSi/m/U7Prr3I359568JGxbG34//neXvL+fud+/mhcMv\nUHi50NLDFEIMcFeuXGHhwoVMnDiRcePG8dZbbxESEsJPfvITJk2aRExMDFlZWSxYsIBRo0bxwgsv\nAFBbW0t8fDxTpkxh/PjxvPvuuxb+JF+fVc1ErnJ1cGXxqMUsHrWYyoZKthVsIz0vnf879H/88dAf\nGes1lsTQRBaELGC4y3BLD1cI0Yt++f5xTpy7fEtfc8zIofxi8dgbPr9lyxZGjhzJ5s2bAaiurubx\nxx8nKCiIQ4cO8S//8i+sXr2aPXv20NDQwLhx41izZg1OTk68/fbbDB06lIqKCuLi4liyZIlVn0Wx\nqplIZzycPEiJTOGVhFf4YPkHPBbzGBrNsweeZf6G+azesprUnFQqGyotPVQhxAAxfvx4tm3bxuOP\nP87HH3+Mu7s7AEuWLGl/PjY2Fjc3N3x8fHB0dKSqqgqtNU888QQTJkxg7ty5lJSUUFpaasmP8rVZ\n5UzkRoa7DOf+sfdz/9j7KbhcQEZeBhl5GfxX5n/xm72/IW5kHEmhScwJmoOLvYulhyuEuAW6mjH0\nloiICLKyskhPT+dnP/sZ8fHxADg6OgJgY2PT/verj1taWli/fj3l5eUcPHgQe3t7QkJCrP6u+wEV\nIh0FDw1mzcQ1fHfCd8mtzCU9L50teVt44pMncLR15K6Au0gKTWJGwAwcbR1v/oJCCGE6d+4cnp6e\nrFq1imHDhvHSSy9167jq6mp8fX2xt7fnww8/pKCg1yq095kBGyJXKaWI9Iwk0jOSH035EYfLD5Oe\nl87W/K1sK9iGi70L8UHxJIYmEjsiFnubvqs5I4SwTkePHuXf/u3fsLGxwd7enj/96U8sX778pset\nXLmSxYsXM378eGJiYoiKiuqD0fYuZU33WsTExOhb1ZSqpa2FfRf2sSVvC9sLtlPTXIOHowfzQ+aT\nGJrIZN/J2Cirv2QkxICUnZ1NdHS0pYfRL3X2/0YpdVBrHdMb79et35JKqQSlVI5S6rRS6j86eV4p\npX5vPn9EKTWlO8cqpR5VSp1USh1XSv331/843WdnY8ftI2/n/93x/9h17y7Wzl5L7IhY3j39Lqu3\nrGbBxgX89sBvOXHxhNzUKIQQN3DT01lKKVvgj8A8oBjYr5R6T2t9osNuiUC4+ScW+BMQ29WxSqnZ\nwFJgota6USnleys/2FfhYOvAnKA5zAmaQ11zHR8WfUhGXgbrstfxyvFXCBkaQkJoAomhiYS5h1lq\nmEII0e9055rIdOC01vosgFLqTYxf/h1DZCnwqja+smcqpYYppUYAIV0c+z3gaa11I4DWuuzWfKSv\nx9nemYVhC1kYtpDqxmq2FWwjIy+DPx/+My8cfoEozygSQxNJDElkhOsISw9XCCEsqjuns/yBog6P\ni81t3dmnq2MjgBlKqb1KqY+UUtO+ysD7grujO8sjlvPXBX9le/J2Hp/2OA42Dvzvwf9l/sb5fDvj\n27xx8g0u1l+09FCFEMIiLLk6yw7wBOKAaUCqUipMX3cBQin1CPAIQFBQUJ8P8ipfZ19WjVnFqjGr\nKKopYkveFtLz0vn13l/zzL5niB0RS2JoIvFB8bg5uFlsnEII0Ze6MxMpAQI7PA4wt3Vnn66OLQb+\noQ37gDbA+/o311q/qLWO0VrH+Pj4dGO4vS/QLZDvTPgOby99m38s+QcPjnuQwsuF/Oee/2TmWzP5\n4c4fsiV/C/Ut9ZYeqhBC9KruzET2A+FKqVCMAPgmsOK6fd4DfmBe84gFqrXW55VS5V0c+w4wG/hQ\nKRUBOAAVX/cD9bVwj3DCPcJ5dPKjHKs41n4Pys6inTjbOTM7aDZJoUncNuI27G3lHhQhxNc3a9Ys\nnn32WWJiemXV7ldy0xDRWrcopX4AbAVsgZe11seVUmvM518A0oEk4DRQBzzQ1bHmS78MvKyUOgY0\nAfdffyrLmiilGO8znvE+43ks5jEOlh4kPS+dbQXb2Hx2M+6O7swLnkdiSCJT/aZia2Nr6SELISyg\npaUFO7uBc593tz6J1jodIyg6bnuhw9818P3uHmtubwJWfZXBWgtbG1umj5jO9BHT+WnsT/n03Kek\n56Wz+exmNuRuwHeIL/ND5pMUmsQ473FWXcFTiMEoPz+fxMRE7rzzTj799FP8/f159913ycnJYc2a\nNdTV1TFq1ChefvllPDw8mDVrFpMmTeKTTz7hvvvu4+jRowwZMoTPP/+csrIyXn75ZV599VU+++wz\nYmNjeeWVVwD43ve+x/79+6mvr2f58uX88pe/tOwH78TAicN+yt7WnpmBM5kZOJO65jp2F+8mPS+d\nt3LeYl32OgLdAkkIMRprjfYYbenhCmF9Mv4DLhy9ta85fDwkPt3lLqdOneKNN97gL3/5CykpKWzc\nuJH//u//5g9/+AMzZ87k5z//Ob/85S957rnnAGhqauJqxY3Vq1dTWVnJZ599xnvvvceSJUvYs2cP\nL730EtOmTePQoUNMmjSJp556Ck9PT1pbW4mPj+fIkSNMmDDh1n7Wr0lCpA852zuTEJpAQmgCl5su\ns6NgBxl5Gfz12F/5y9G/EO4RTlJoEgkhCQS4BVh6uEKILoSGhjJp0iQApk6dypkzZ6iqqmLmzJkA\n3H///SQnJ7fvf++9937h+MWLFxunwcePx8/Pj/HjxwMwduxY8vPzmTRpEqmpqbz44ou0tLRw/vx5\nTpw4ISEiDEMdhrIsfBnLwpdRUV/BB/kfkJGXwdqstazNWssE7wntjbV8nPvHqjQh+qWbzBh6S8dS\n77a2tlRVVXW5v4vLF9tP3KxsfF5eHs8++yz79+/Hw8OD1atX98uy8VJhsB/wHuLNiugVvJb0Glvv\n2cqPpvyIxtZGntn/DHM3zOXhrQ+zMXcj1Y3Vlh6qEOIG3N3d8fDw4OOPPwbgtddea5+V9MTly5dx\ncXHB3d2d0tJSMjIybtVQbymZifQzI11H8tD4h3ho/EOcrTpLel46GXkZPPnZk/xq76+4c+SdJIYm\nMitwFs72zpYerhCig7///e/tF9bDwsL429/+1uPXmjhxIpMnTyYqKorAwEDuuOOOWzjSW2fQloK3\nJlprTlw6QcbZDDLyMyirK2OI3RBmBcwiITSBO/3vxMHWwdLDFKLPSCn4G+vrUvAyE7ECSinGeo1l\nrNdY/jXmX8kqzSIjL4MPCj4gIz8DNwc35gbNJTE0kenDp8s9KEKIPiMhYmVslA0xw2OIGR7Df8T+\nB5nnMtmSv4UPCj7g7dNv4+XkxYKQBSSGJjLRZ6LcgyKE6FUSIlbM3saeGQEzmBEwg/9s+U8+LvmY\njLwMNuRu4PWTr+Pv6k9CiNEHJcIjQgJFCHHLSYgMEE52TswLnse84HnUNtWys2gn6XnpvHL8Ff56\n7K+EuYeRGJpIUmgSQUMtVw1ZCDGwSIgMQK4OriwZtYQlo5ZwqeES2/K3kZ6Xzh8P/ZE/HvojY73G\nkhiaSEJIAn4ufpYerhDCisnqrEHkwpUL7X1Qsi9lo1BM9ZtKYmgi84Ln4eHkYekhCtEtsjrrxvp6\ndZaEyCCVX51PRn4GGXkZ5FXnYafsiBsZR1JoEnOC5uBi73LzFxHCQvpDiLi6ulJbW2vRMXRGlviK\nPhHiHsL3Jn6PNRPWkFOZQ3peOlvytvDEJ0/gaOvIXQF3kRSaxIyAGTjaOt78BYUQN6S1RmuNjc3A\nKxIiITLIKaWI8owiyjOKH035EYfLD5N+Np0PCj5gW8E2XO1dmRM0h6TQJGJHxGJnIz8yQnRUW1vL\n0qVLqayspLm5mV/96lcsXbqU/Px8FixYQGxsLAcPHiQ9PZ3t27fzzDPPMGzYMCZOnIijoyPPP/88\n5eXlrFmzhsLCQgCee+65fnuH+vXkdJboVEtbC/su7CMjL4MdBTuoaa7Bw9GD+SHzSQxNZLLvZGzU\nwPtWJaxDx1M2z+x7hpOXTt7S14/yjOLx6Y93uc/V01ktLS3U1dUxdOhQKioqiIuL49SpUxQUFBAW\nFsann35KXFwc586d4/bbbycrKws3NzfmzJnDxIkTef7551mxYgX/9E//xJ133klhYSELFiwgOzu7\nR2OX01miX7CzseP2kbdz+8jb+Vncz/ik5BMy8jJ49/S7vJXzFsNdhrffgxLtGS33oIhBS2vNE088\nwe7du7GxsaGkpITS0lIAgoODiYuLA2Dfvn3MnDkTT09PAJKTk8nNzQVg+/btnDhxov01L1++TG1t\nLa6urn38ab46CRFxU462jsQHxRMfFM+V5it8WPQhGXkZrDuxjleOv0LI0BASQxNJDE0k1D3U0sMV\ng8zNZgy9bf369ZSXl3Pw4EHs7e0JCQlpL9l+ffn3G2lrayMzMxMnJ6feHGqvkPMR4itxsXdhUdgi\n/hj/R3bdu4tf3PYLfJ19eeHwCyx5Zwkp76fwt2N/43zteUsPVYg+UV1dja+vL/b29nz44YcUFBR0\nut+0adP46KOPqKyspKWlhY0bN7Y/N3/+fP7whz+0Pz506FCvj/tWkZmI6DF3R3eWRyxnecRyyurK\n2Jq/lYy8DH538Hf87uDvmOw7mcTQROYHz8driJelhytEr1i5ciWLFy9m/PjxxMTEEBUV1el+/v7+\nPPHEE0yfPh1PT0+ioqJwd3cH4Pe//z3f//73mTBhAi0tLdx111288MILffkxekwurItbruhyUfs9\nKKerTmOrbIkdEUtiaCLxQfG4ObhZeojCyvWH+0R64up1jpaWFpYtW8aDDz7IsmXLbul7yM2GXZAQ\nsT65lblk5BmBUlJbgoONAzMCZpAYmsjMgJk42VnfOWBhedYaIo899hjbt2+noaGB+fPns3bt2lu+\nKEVWZ4kBJcIjggiPCP558j9zpOIIW/K2sCV/CzsKd+Bs58zsoNkkhSZx28jbsLext/RwhehVzz77\nrKWHcMtJiIg+oZRios9EJvpM5LGYxzhQeoCMvAy2FWxj89nNuDu6My94HkmhSUzxnSKNtcRNaa1l\nafl1LHFmSU5nCYtqbm1mz7k9pOels6toF/Ut9fgO8WVB6AKSQpMY6zVWflGIL8nLy8PNzQ0vLy/5\n+TBprbl48SI1NTWEhn5xqb3Fr4kopRKAtYAt8JLW+unrnlfm80lAHbBaa53V1bFKqSeB7wDl5ss8\nobVO72ocEiIDW11zHbuLd5Oel84nJZ/Q3NZMoFsgCSEJJIUmMdpjtKWHKPqJ5uZmiouL2+/HEAYn\nJycCAgKwt//iqWGLhohSyhbIBeYBxcB+4D6t9YkO+yQBj2KESCywVmsd29WxZojUaq27fZJQQmTw\nqG6sZmeh0Vhr34V9tOk2wj3CSQpNIiEkgQC3AEsPUQirYekL69OB01rrs+Zg3gSWAic67LMUeFUb\niZSplBqmlBoBhHTjWCG+xN3RnWXhy1gWvoyK+or2e1DWZq1lbdZaJvhMICk0iQUhC/Ae4m3p4Qox\naHXnjnV/oKjD42JzW3f2udmxjyqljiilXlZKSUck0SnvId6sjF7JuqR1bLlnCz+c8kMaWxp5et/T\nxKfF8/DWh9mYu5HqxmpLD1WIQceSZU/+BIQBk4DzwG8720kp9YhS6oBS6kB5eXlnu4hBxN/Vn4fH\nP8yGJRt4Z+k7fGf8dzh/5TxPfvYks1Jn8eiOR0k/m05dc52lhyrEoNCd01klQGCHxwHmtu7sY3+j\nY7XWpVc3KqX+Amzq7M211i8CL4JxTaQb4xWDxKhho/jB5B/w/Unf58TFE0Zjrfwt7CrexRC7IcwK\nmEViaCJ3+N+Bg62DpYcrxIDUnRDZD4QrpUIxAuCbwIrr9nkP+IF5zSMWqNZan1dKld/oWKXUCK31\n1Sp9y4BjX/vTiEFJKcVY77GM9R7Lj2N+zMHSg+33oGTkZ+Dm4Ma84HkkhiYyzW+a3IMixC3U3SW+\nScBzGMt0X9ZaP6WUWgOgtX7BXOL7PJCAscT3Aa31gRsda25/DeNUlgbyge92CJVOyeos8VU0tzXz\n2bnPyMjLYGfhTupa6vBy8mJByAISQxOZ6DNR7jEQg4LF7xPpLyRERE81tDSwu3g3GXkZ7C7eTVNb\nE/6u/u2NtSI8IiRQxIAlIWKSEBG3Qk1TDTsLd5KRl0Hm+UxadSuj3Ee1N9YKGhpk6SEKcUtJiJgk\nRMStdqnhEh/kf0BGXgZZZVkAjPMaR0JoAgkhCfi5+Fl4hEJ8fRIiJgkR0ZsuXLnAlrwtpOelk30p\nG4Viqt/U9sZaw5yGWXqIQvSIhIhJQkT0lbzqvPZAyb+cj52y47aRt5EYmsicoDm42Hevd7YQ/YGE\niElCRPQ1rTUnL500GmvlZ3DhygUcbR2ZGTCTxNBEZgTMwNHW0dLDFKJLEiImCRFhSW26jUNlh8jI\ny+CDgg+41HAJV3tX5gTNYWHoQuJGxmGjLFkEQojOSYiYJEREf9HS1sK+8/tIz0tnR+EOaptr8Xf1\nZ3nEcu4efbcUhRT9ioSISUJE9EeNrY3sLNxJWm4a+y/sx87GjrlBc0mJTCHGL0buPxEWJyFikhAR\n/d3Z6rOk5aTx7pl3qWmqIWRoCCmRKSwZtQR3R3dLD08MUhIiJgkRYS0aWhrYmr+V1NxUjpQfwdHW\nkQUhC0iJTGGC9wSZnYg+JSFikhAR1ujkpZOk5aSx6ewm6lrqiPSIJCUyhYVhC2WpsOgTEiImCRFh\nza40X2Hz2c2k5qSSU5mDs50zC8MWkhKZQpRnlKWHJwYwCRGThIgYCLTWHK04SmpOKlvyt9DY2sgE\n7wkkRyazIGQBQ+yGWHqIYoCREDFJiIiBprqxmvfPvE9qbip51Xm4ObixdNRSkiOSCRsWZunhiQFC\nQsQkISIGKq01B0oPkJaTxrbCbbS0tTBt+DRSIlKID4rH3tbe0kMUVkxCxCQhIgaDi/UXeef0O6Tl\nplFSW4KnkyfLRi/jnoh7CHQLvPkLCHEdCRGThIgYTNp0G5+d+4zUnFR2Fe9Ca83t/reTEpHCXQF3\nYWfTne7WQkiItJMQEYPVhSsX+Mepf7AxdyNl9WX4OvuyPHw53wj/hvQ8ETclIWKSEBGDXUtbCx8V\nf0RaThp7zu3BVtkyM2AmKZEp3DbyNikAKTolIWKSEBHimqKaIjbkbuCd0+9wqeES/q7+JEckc/fo\nu/Ea4mXp4Yl+RELEJCEixJc1tTaxo3AHqTmpHCg9gJ2NHfOC5pEcmSwFIAUgIdIuJshZH3jpXyFq\nMQRMAxuZugvR0dmqs6TlXisAGeoeSkpECotHLZYCkIOYhIgpJtRDH3jAFtqawcUXopIgahGE3gV2\n0l1OiKvqW+rZmr+VtJw0jlQYBSATQhJIiUxhvPd4mZ0MMhIippiYGH3gkx1wahtkvw+nt0NTLTi4\nQfg8iFoI4fPBaailhypEv5F9MZu0XKMAZH1LPVGeUSRHJEsByEHE4iGilEoA1gK2wEta66eve16Z\nzycBdcBqrXVWN4/9MfAs4KO1ruhqHF+6JtLcAHkfwclNcDId6irAxh7CZhozlMgkcJPlj0IA1DbV\nkp6Xzls5b5FbmYuznTOLwhaREplCpGekpYcnepFFQ0QpZQvkAvOAYmA/cJ/W+kSHfZKARzFCJBZY\nq7WOvdmxSqlA4CUgCpj6lUOko7ZWKNpnBsomqMwHFARON2YoUYvAa1TX/zeEGAS01hypOEJqTipb\n87caBSB9JpASkcKCkAU42TlZeojiFrN0iNwGPKm1XmA+/gmA1vo3Hfb5M7BLa/2G+TgHmAWEdHWs\nUmoD8F/Au0DMzULEMzhaP/PaZuKjfQnzcb3xjlpD2QnINgPlwhFju080RC8yQmXEJJDzwmKQq26s\n5r0z75Gak0r+5fxrBSAjkwlzlwKQA4WlQ2Q5kKC1fth8/C0gVmv9gw77bAKe1lp/Yj7eATyOESKd\nHquUWgrM0Vr/UCmVTzdCxD0wSnus/C0AYT4uzI32Iz7Kl6nBHtjZdrFSq7IActKNUCn8FHQbDA0w\nwiR6EQTdDrZSQkIMXlcLQKbmpLK9cLsUgBxgejNELPKbUynlDDwBzO/Gvo8AjwAEBQXx8b/PZufJ\nMrZnl/K3PXm8uPssw5ztmRXhQ3y0HzMjfRjqdN0PvEcwxH3P+HPlIuRmwMnNkPV32PdnGOIBEQnG\nKa9Rc8DBuRc+tRD9l1KKacOnMW34NCrqK3jn9DtsyN3Av+3+t/YCkMsjlhPgFmDpoYp+xiKns4BN\nwA6Mi/AAAcA5YLrW+sKNxnL9NZGahmY+PlXB9uxSPjxZRmVdM3Y2itgwT+Kj/Jgb7UeQVxeB0HQF\nTu8wTnnlboGGarAbAqPjjUCJWADOnl3+/xFioGrTbXx67lNSc1L5qPgjKQBpxSx9OssO4+J4PFCC\ncXF8hdb6eId9FgI/4NqF9d9rrad351jz+Hy6cTqrqwvrrW2azwsr2Z5tzFJOl9UCEO7rytwxfsyN\n9mVSoAe2Nje4DtLaDPmfGDOUk5uh5hwoWwi549pKr2FShlsMTlIA0rr1hyW+ScBzGMt0X9ZaP6WU\nWgOgtX7BXOL7PJCAMbt4QGt94EbHdvL6+XzNELlewcUrbM8uY0d2KfvyLtHSpvF0cWB2pC9zo32Z\nEeGDq+MNvklpDeeyjDDJ3gQVOcb2EZOMQIleBD5RcmFeDDpSANI6WTxE+oue1s6qrm/mo9xydmSX\nsiunnOr6ZhxsbYgb5cXcaF/io/3wH9ZFX+uKU+bS4c1QvN/Y5hlmBErUIinBIgalostFbDh1rQBk\ngGsAyyOWSwHIfkhCxHQrCjC2tLZxoKCSHdmlbM8uI6/iCgBRw92YN8aP+Gg/Jvi7Y3Oj016Xzxsr\nvU5ugrzd0NbSoQTLYgidISVYxKAiBSD7PwkRU29U8T1TXtseKAfyL9GmwcfNkTmRvsRH+3JnuDfO\nDjc47VVfZZRgObnJ+G/zFaMES8R8Y/nw6HlSgkUMKlIAsn+SEDH1din4yitN7MotY3t2Gbtzyqlp\nbMHRzoY7RnsTH+1LfJQfw91vcDfv1RIs2e9DToZRgsXWAUJnmnfMLwRX314buxD9iRSA7F8kREx9\n2U+kqaWN/fmX2J5dyvbsUoou1QMw3t+d+Ghf5kb7MXbk0M7/MbS1QtFe88L8+1BVwLUSLOYd81KC\nRQwSUgDS8iRETJZqSqW15lRZLduzS9mRXUZWYSVaw/ChTsyJ9mVetB+3jfLCyd62s4Oh9Li5dPh9\nuHDU2O475lpNrxETZaWXGPCuLwDpYu/CorBFJEckSwHIXiYhYuovnQ0rahv58GQZO7LL2H2qnLqm\nVobY23JnuDdzo32ZHeWLr9sNTntV5hsVh09ugsLPjBIs7oHXTnlJCRYxwHVWAHKiz0RSIlOYHzxf\nCkD2AgkRU38JkY4aW1rJPHuJ7SdK2ZFdyrnqBgAmBg5jnrl8OGq4W+enva5UGHfKZ2+CMzuhtdEs\nwZJo3IsSNltKsIgB7foCkEMdhrJ09FKSI5IJdQ+19PAGDAkRU38MkY601mSfrzFWe50s43BRFQD+\nw4YYF+aj/YgL88TRrpPTXo21cGaHcdpLSrCIQUZrzf4L+0nNTWVHwQ5adAvTh08nOTKZ+EApAPl1\nSYiY+nuIXK/scoNZLLKMT06X09DchouDLXeZxSJnR/rg5drJPSXtJVjMGxxrzncowbLYuCfFXQrh\niYGpYwHIktoSPJ08+Ub4N7gn/B4pANlDEiImawuRjhqaW/n0TEV7KZbSy40oBVOCPJgbbdT2Gu3r\n+uXTXm1tcO7za822KnKN7SMmmb1RpASLGJha21qNApC5qewu3o3Wmjv87yAlIoUZATOkAORXICFi\nsuYQ6UhrzbGSy8Zqr5OlHCu5DECQp3P78uHpoZ7Yd9YjpTz32gylxPx/4TnK7I2yGPxjpASLGHAu\nXLnAxlMb2Zi7kfL6cvyc/bgn4h6+MVoKQHaHhIhpoITI9c5X17PDnKHsOXORppY23JzsmBnhw9xo\nP2ZF+jDM2eHLB14+Z5Zg2XytBIurn1FxOGoRhN4Fdp0cJ4SVam5rZnfRblJzU/n03KfYKltmBc4i\nJSKFuJFxUgDyBiRETAM1RDqqa2rhk1MVRqicLKOithFbG8XUYA/mRfvduDVwewmW9+HUdqMEi+NQ\nCJ9nBEr4PHB06/sPJEQvKbpcRNqpNN459Q6VjZUEuAaQHJnM3aPvxtNJFqF0JCFiGgwh0lFbm+Zw\ncRU7zB4pJy/UABDm7dK+2iums9bAzQ1wdpdx2qtjCZawWcZpr8gkKcEiBoym1ia2F2wnNTeVg6UH\nsbexZ27wXFIiUpjqN1VKrCAh0m6whcj1iivr2gMl8+xFmls17kPsmR3ZRWvgqyVYsjcZs5SqQowS\nLLHXesyQw+qwAAAgAElEQVR7hlnk8whxq52pOkNabhrvnX6PmuYawtzDSIk0CkAOdRi8xVAlREyD\nPUQ6qm1s4ePccrZnl7HzZGl7a+DpoZ7maq9OWgNrDaXHrjXbKu1YgsWs6SUlWMQAUN9Sz5a8LaTl\npnG04ihOtk4khCaQEpHCOO9xg252IiFikhDpXMfWwDuySznVoTVwfLQf88bcoDVwZf61dsBfKsGy\nCIJukxIswuqduHiCtNw0Np/dTH1LPdGe0SRHJrMwdCHO9oOjIoSEiElCpHt61Br4SoVx/eTk5g4l\nWDwhMtEIlVFzwL6L7o9C9HO1TbVsPruZt3Lf4lTlqUFVAFJCxCQh8tVV1zezO7ec7de1Bo4N82zv\n5Pil1sBXS7Bkb4LcrdBYDfbORpBEL4bw+VKCRVgtrTWHyw+TlpvGlrwtNLU1DfgCkBIiJgmRr6dj\na+Ad2WWc7dAaeK65fHhiwLAvtgZuaYKCT66d9movwXLnteso7v4W+kRCfD1VDVW8d+Y90nLTBnQB\nSAkRk4TIrdVZa2BvV0fio27QGri9BMv7RqBcLcEycrIZKIvAJ1IuzAur01kByNjhsSRHJjMncI7V\nF4CUEDFJiPSeqromduUYp70+6tAa+PZRXswd49d5a+D2EiyboOSgsc1r9LUL81KCRVihqwUg03LS\nOHflHF5OXkYByIh78He1zlm3hIhJQqRvdGwNvCO7jMJLdQCM8x9KfJSxfHic/3WtgS+fu3bKK//j\nL5ZgiV4EIVKCRViXzgpA3ul/JymRKczwn4GtTSctHfopCRGThEjf01pzuqyWbTdoDTw32pfbR3l/\nsTVwfRWc+sCYoXyhBMt8Y5YiJViElTlfe56Npzbyj1P/oLy+nOEuw7kn/B6+Ef4NfJ37f/UHCRGT\nhIjlXaxt5MOccnZkl7I7t5wrZmvgO0YbrYHnRF/XGri5Hs5+ZFxHycmAuosdSrAsMpYQSwkWYSWa\n25r5qOgjUnNS+ez8Z9gqW2YHziY5Mpm4Ef23AKTFQ0QplQCsBWyBl7TWT1/3vDKfTwLqgNVa66yu\njlVK/RewFGgDysxjznU1DgmR/uVqa+Crq71KquoBozXw3Cijtlf0iA6tgdtaoTDTPO3VoQRLUNy1\n6yieA2dFjBjYCi8XsiF3A2+ffpuqxioC3QJJjkhm6eil/a4ApEVDRCllC+QC84BiYD9wn9b6RId9\nkoBHMUIkFlirtY7t6lil1FCt9WXz+H8Gxmit13Q1FgmR/ktrzckLNWw/0c3WwFdLsGSbvVHaS7CM\nNZttLYThE2Sll+j3mlqb2FawjdScVLLKsrC3sWde8DxSIlOY4julX5RYsXSI3AY8qbVeYD7+CYDW\n+jcd9vkzsEtr/Yb5OAeYBYTc7NgO24O01t/raiwSItajrKaBD0+Wse3EF1sDzwj3Ye6YTloDX8oz\neqNkbzJKsKDBPcicoSyUEizCKpyuPE1abhrvn3mfmuYaRrmPIjky2eIFIC0dIsuBBK31w+bjbwGx\nWusfdNhnE/C01voT8/EO4HGMELnhsUqpp4BvA9XAbK11eVdjkRCxTl21Br7ayTG8Y2vg2nLIvVqC\n5cPrSrAsglGzpQSL6NfqmuvYmr+V1JxUjl08hpOtE4mhiaREpjDWa2yfz04GbIh0OP4ngJPW+hed\nvP8jwCMAQUFBUwsKCnr4UUV/8JVbAzfWwOkdxkqv3A+ulWAZHW8ESsQCGOJhwU8kRNeOXzxOWk4a\n6Xnp7QUgUyJTSApN6rMCkJYOkb44nRUEpGutx3U1FpmJDDwXqhvYcbKU7Sc6tAZ2tOOuSB/mXd8a\nuKXJuAfl6v0otRfAxg6C7zBqekUmSQkW0W/VNNUYBSBz3uJ01en2ApApkSlEeET06ntbOkTsMC6O\nxwMlGBfHV2itj3fYZyHwA65dWP+91np6V8cqpcK11qfM4x8FZmqtl3c1FgmRga2r1sBzzYvzo662\nBm5rg3NZxgwlexNcPGVsHznFbLa1GLwj5MK86HeuFoBMzUlla/5WmtqamOQzySgAGTIfR1vHm7/I\nV9QflvgmAc9hLNN9WWv9lFJqDYDW+gVzie/zQALGEt8HtNYHbnSsuX0jEImxxLcAWKO1LulqHBIi\ng0dbm+ZISbWx2qs7rYHLc6/V9PpCCRazppf/VCnBIvqdqoYq3j3zLmm5aRRcLsDd0Z2lo4wCkCHu\nIbfsfSweIv2FhMjgVVxZx86TZWzPLiPzzEWaWttwH2LPrEgf5l7fGri6xFjpdXIT5H9ilmAZDlFJ\nxixFSrCIfkZrzb4L+0jNSWVn4c5bXgBSQsQkISLgi62BP8wp49KVpvbWwPHRfsyN9iXYy8XYub7S\nuCB/chOc3g7NdddKsEQvgtFzpQSL6Fcq6it4+9TbbMjdcMsKQEqImCRExPVa2zSHiirZdqLz1sBz\no32ZHGS2Bm6uh7O7jGsoOelQfwlsHc0SLAuNC/OuPpb8OEK0a21rZc+5PaTlpLG75FoByHsj7+VO\n/zu/UgFICRGThIi4mcKLdWzPNq6jdGwNfPW0111XWwO3tkCRWYIlexNUdyzBYt4xLyVYRD/RWQHI\n5eHL+Ub4N/BxvvkXHwkRk4SI+CqutgbekV3Kh9e1Br7ayTHAw9kowXLhqLl0eJNRjgXAb9y1ml7D\nx8tKL2FxnRWAnBM0h+SIZGJHxN6wAKSEiElCRPTUV2oNfCnv2r0oV0uwDAu6NkMJjJMSLMLiri8A\nGeQW1F4A0sPpizfgSoiYJETErXK2vJYd2WVsyy7lYEElrW0ab1dH5kT5EB/tx4yrrYGvlmDJ3mRc\nT2ltBGcviEg0LsyHzZISLMKiGlsb2VawjbSctPYCkPND5pMSkcJk38kopSRErpIQEb2hs9bADnY2\n3DHKi3hzljLCfYhZgmW7MUPJ3QqNl6+VYIleaoSKBIqwoFOVp9oLQNY21zJ62GiSI5JZOWalhAhI\niIje1+3WwK3NZgmWTXAy3SjBMsQDJq2EmAfBa5SFP4kYzOqa69iSv4XUnFSOXzzOsdXHJERAQkT0\nrautgbdnl7E9u7S9NbDfUEfmRPkxb4zZGthWGYFy4K/GLKWtBcJmw7SHjNNecv1EWNDxiuOM8xkn\nIQISIsKyOmsN7GRvw52jfZg/xo+FE0bg0lQBWa/CwVfgcgm4jYSp98OU+2HoCEt/BDFIyTURk4SI\n6C8aW1rZe/baaa+SqnpcHe1YNtmfVXHBRPoMgdwtxuzkzE5QtsbKrmkPQehMWS4s+pSEiElCRPRH\nWmuyCitZn1nIpqPnaWppY1qIB6vigkkYNxzH6nw4+Df4fJ1RhsVrtHHdZNIK6YUi+oSEiElCRPR3\nl640kXagiNf3FVJwsQ4vFweSYwJZGRtEoJuC4+8Ys5Pi/WA3BMbdA9MeNKoMC9FLJERMEiLCWrS1\naT4+XcG6zAJ2ZJeigZkRPqyKDWZ2lC+2pUeNMDmSBs1XYMQk41TXuOXg0Dfd7sTgISFikhAR1uhc\nVT1v7i/izX2FlNU04j9sCPdNDyRlWiC+9o1w+C0jUMpPgqO7cZor5kHw6d1ud2LwkBAxSYgIa9bc\n2sb2E6Ws21vAntMXsbNRLBg3nFWxwcSFeqAKPzPC5MR70NYMITOM2UnUIvia/STE4CYhYpIQEQPF\nmfJaXt9byIaDxVTXNzPa15WVsUF8Y0oA7q2V15YJVxeBq5+xRHjq/eAeYOmhCyskIWKSEBEDTUNz\nK+8fPse6vYUcLqpiiL0tSyaOZFVcMONHusKpbcbs5NQ2Y1lwRKJxIT5sjrT7Fd0mIWKSEBED2dHi\natbvLeCdQyU0NLcxMcCdlXHBLJ4wkiFXioyZSdZrUFcBHqHGdZPJq8DZ09JDF/2chIhJQkQMBtX1\nzbydVcy6vYWcLqtlqJMdy6cGsjIuiFEe9pD9Pux/yShTb+sIY5cZ104CpslNjKJTEiImCRExmGit\n2Zt3iXWZBWw9foHmVs3to7xYFRfMvDF+2Fdkw4GXjdVdTTXgN94Ik/HJ4Ohq6eGLfkRCxCQhIgar\n8ppGUg8U8freQkqq6vF1c+Sb0wK5LzaIEU4tcDQN9v/V6Mro4AYTv2kEim+0pYcu+gEJEZOEiBjs\nWts0u3LKWJdZwK7cchQQH+3HqrhgZozywqZkv3Eh/vjb0NoEQbcbYRK9GOwcLT18YSESIiYJESGu\nKbpUx+v7CkndX8TFK00EezmzYnoQyTGBeFIDh9YZp7sq88HZG6Z8G6auBo9gSw9d9DEJEZOEiBBf\n1tjSypZjF1ifWci+/Es42NmwcPwIVsUFMSXQHXXmQ2N2krsFtIbw+cbsZPRcsLG19PBFH7B4iCil\nEoC1gC3wktb66eueV+bzSUAdsFprndXVsUqp/wEWA03AGeABrXVVV+OQEBGiazkXali/t4B/ZJVQ\n29hC1HA3VsUFc/dkf1zrz0PW3+Hg3+FKGQwLgqkPwORvgauPpYcuepFFQ0QpZQvkAvOAYmA/cJ/W\n+kSHfZKARzFCJBZYq7WO7epYpdR8YKfWukUp9QyA1vrxrsYiISJE91xpbOHdQ+dYl1nAifOXcXW0\n4+7Jxk2MUd5ORlvfAy8bHRlt7GHMUmN2EnSbLBMegHozRLrTt3M6cFprfdYczJvAUuBEh32WAq9q\nI5EylVLDlFIjgJAbHau1/qDD8ZnA8q/7YYQQBhdHO1bEBnHf9EA+L6piXWYBqQeKWZdZSEywB6vi\nYklctRTHytNGmBx6HY5tAN8xxk2ME+4Fp6GW/hjCCnSnboI/UNThcbG5rTv7dOdYgAeBjM7eXCn1\niFLqgFLqQHl5eTeGK4S4SinFlCAPfpcyib0/ieenSdFU1Dbyo7cOcdtvdvKbA20UTv8F/PgkLP69\nUegx/TH4bRS8/yO4cNTSH0H0c92ZifQqpdRPgRZgfWfPa61fBF4E43RWHw5NiAHFw8WB79wVxkN3\nhrLnjNHr5KWP83hx91nuCvdhVVwCcx7+FrbnPzcuxB9+w+jIGDDdONU15m6wd7L0xxD9THdCpAQI\n7PA4wNzWnX3suzpWKbUaWATEa2taJiaEFbOxUcwI92FGuA/nq+t5c18Rb+4v5DuvHmCkuxP3TQ/i\n3vjf4Tv/V0aQ7P8rvP1d2PITo1ZXzAPgGWbpjyH6ie5cWLfDuDgejxEA+4EVWuvjHfZZCPyAaxfW\nf6+1nt7Vseaqrd8BM7XW3TpPJRfWhegdza1t7MguZV1mIZ+crjB6nYwdzsq4IG4L9UDl7TZmJyfT\nQbfCqHhjdhK+AGwtfkJD3ER/WOKbBDyHsUz3Za31U0qpNQBa6xfMJb7PAwkYS3wf0FofuNGx5vbT\ngCNw0XybTK31mq7GISEiRO87a/Y6STN7nYzycWFlbDD3TA3Avbn8Wq+TmvMwNMC4gXHKt8HNz9JD\nFzdg8RDpLyREhOg7Dc2tbDpynnWZBRwqqsLJ3qa918mEES6Qk2HMTs7uAhs7owPjtIeMjoyyTLhf\nkRAxSYgIYRnHSsxeJ5+fo765lQkB7qyKDWbxxJEMuZxnXID/fB00VIF3pLFMeOI3YcgwSw9dICHS\nTkJECMu63NDM21klrMss4JTZ6+SeqQGsjA1mtIctHPuHMTspOQj2zjDuHmN2MnKypYc+qEmImCRE\nhOgftNbsy7vEur2FbDl2nuZWzW1hRq+T+WP9sC89YoTJ0Q3QXAcjpxhhMvYb4OBs6eEPOhIiJgkR\nIfqf63ud+FztdTI9iJGOjXD4TSNQKnLBaRhMWmmc7vIebemhDxoSIiYJESH6r9Y2zUe5ZazLLOTD\nnDIUMCfKj1VxQdw12hubwj1GmGS/D20tEDrTmJ1EJhl3yoteIyFikhARwjoUXarjjX2FpB4ooqK2\niSBPZ1bEBpESE4hnWyV8/ioceAUuF4PbCJhyP0y9H4aOtPTQByQJEZOEiBDWpamljS3HL7Aus4B9\neZdwsLUhafxwVsUFMzVwKOrUB8bs5PQOUDYQmWjMTkJngU13SvuJ7pAQMUmICGG9cktrWJ9p9Dqp\nMXudrIwLZtlkf1yvFF1bJlx3ETxHGddNJq0AZ09LD93qSYiYJESEsH5XGlt477DR6+T4ucu4ONhy\n92R/VsUFE+3tACfeNWYnRXvBzslY0TXtIfCfKjcx9pCEiElCRIiBQ2vNoaIq1mUWsunIORpb2pga\n7MGquCASx43A6WK2ESZHUqGpFoZPgGkPw/jl4OBi6eFbFQkRk4SIEANTVV0TGw4Ws35vIXkVV/Bw\nticlJpAVsUEEu7TCkbeM5lllJ8DR3bgbftpD4BNp6aFbBQkRk4SIEANbW5vm0zMXWZdZwLbsUlrb\nNHdF+LAqNog5kT7YlewzZicn3oXWJgi+0wiTqEVg52Dp4fdbEiImCREhBo8L1Q28ub+QN/YVUnq5\nkRFmr5NvTgvE16YGDq0zZidVheDia1QSnroahgXe9LUHGwkRk4SIEINPS2sb27PLWL+3gI9PGb1O\n5o/1Y1VsMLeFDkOd2WnMTnK3GhfewxcY105GzZFlwiYJEZOEiBCDW17FFV7fW0DawWKq6poJM3ud\nLJ8SgHvTeaPPSdarcKUcPEJg6gMw+Vvg4mXpoVuUhIhJQkQIAUavk81HzrNubwGfFxq9ThZPMHqd\nTBzhDCffN9r6FuwBWwejP/y0hyFw+qBcJiwhYpIQEUJc7/i5atZlFvLuoRLqmloZ7+/Oqrgglkz0\nZ0hVrnHd5PCb0HgZ/MYZNzFOSAFHN0sPvc9IiJgkRIQQN3K5oZl3Pjd6neSW1uLmZMc9UwJYFRfE\naHcFxzbA/pfgwlFwcIUJ9xoru/zGWnrovU5CxCQhIoS4Ga01+/MrWZdZQIbZ6yQuzNPodRLth8OF\nLONC/LF/QGsjBMYZp7rGLAE7R0sPv1dIiJgkRIQQX0VF7bVeJ8WV9Xi7mr1OYoPwd6iHQ+uNayeV\neeDsDZNXQcwDxkX5AURCxCQhIoToidY2ze7cctZlFrCzvdeJLyvjgpk52gubvF3GtZOcdNAaRs81\nTnWFzwcbW0sP/2uTEDFJiAghvq7iSqPXyVv7jV4ngZ5DWDE9mJSYALxaKyDr73Dw71B7AdyDjD4n\nU74Nrr6WHnqPSYiYJESEELdKU0sbW81eJ3vNXieJZq+TmABXVE66ce0kbzfY2EP0YmN2EnyH1S0T\nlhAxSYgIIXrDqdIa1u8tZOPBYmoaW4j0c2NVXBB3T/bHrTbfONV1aD00VINPFMQ8BBPvBSd3Sw+9\nWyweIkqpBGAtYAu8pLV++rrnlfl8ElAHrNZaZ3V1rFIqGXgSiAama61vmg4SIkKI3lTX1MJ7h86x\nbm8Bx0qMXidLJ/uzKjaYMd52cGyjMTs59znYuxhl6ac9BCMmWnroXbJoiCilbIFcYB5QDOwH7tNa\nn+iwTxLwKEaIxAJrtdaxXR2rlIoG2oA/A49JiAgh+gutNYeLq1mXWcD7h41eJ1OChrEqLpik8SNw\nKjtshMnRjdBSDwHTjNnJ2GVg72Tp4X+JpUPkNuBJrfUC8/FPALTWv+mwz5+BXVrrN8zHOcAsIKQb\nx+5CQkQI0U9d7XXy+t5Czpq9TpJjAlkxPYgQlybjbvj9f4WLp2CIB0xaadwV7zXK0kNv15sh0p0S\nl/5AUYfHxea27uzTnWOFEKLfGubswMMzwtjx45msfziWuDAv/vpJHrOe3cW3Xs9lq9syWr63F779\nHoTeBXtfgD9MgdeWQfYmaG2x9EfoVXaWHsDNKKUeAR4BCAoKsvBohBCDlVKKO0Z7c8dob0ovN/Dm\nviLe2FfId187yAh3J745LYhvJryIX2KVUUn44Cvw1koY6g9TzGXCQ0dY+mPcct2ZiZQAHbu8BJjb\nurNPd47tktb6Ra11jNY6xsfH56scKoQQvcJvqBM/nBvOJ4/P5s/fmspoX1f+d3sutz+9k++9W8Ie\n/wfRPzwM9643Wvju+jU8Nw5Svw1nPzJuaBwgujMT2Q+EK6VCMQLgm8CK6/Z5D/iBUupNjAvr1Vrr\n80qp8m4cK4QQVsnO1oYFY4ezYOxw8iuu8Pq+QtIOFJFx7AJh3i6siB3D8nveYlh9ERz8G3y+zmjt\n6xVuXDeZdJ9xHcWKdXeJbxLwHMYy3Ze11k8ppdYAaK1fMJf4Pg8kYCzxfeDqhfLOjjW3LwP+APgA\nVcChqxfgb0QurAsh+ruG5lbSj55nXWYBWYVVONrZsHii2evEzwF14l1jZVfxfrAbAuPugWkPgv/U\nXhuTxe8T6S8kRIQQ1uTEucus21vAO58bvU7G+Q9lVWwwSyaNxPniCSNMjqRB8xUYOdlYJjzuHnBw\nvqXjkBAxSYgIIaxRTXuvk0JySmvae52sjA0i3L0NDr9lBEr5SeMu+IkrjNNdPhG35P0lREwSIkII\na6a15kCB2evk6AWaWtuIDTV6nSwY44dDSaYRJifeg7ZmCJlh9DqJWgi29j1+XwkRk4SIEGKguFjb\nSOqBYl7fV0DRJaPXyb3TArhvehAB9rXw+Wtw4BWoLgTX4cYS4amrwf2r32onIWKSEBFCDDRtbZqP\nTpWzPrOAnSfLAJgd6cuquGDuGu2J7Zntxuzk1DajenBEolGvK2w22HTnLg0JkXYSIkKIgaykqp43\n9hby5v4iKmobCfAYworYIFJiAvFuPm/cwJj1GtRVgEeocd1k8ipw9uzydSVETBIiQojBoKmljQ9O\nGL1OMs9ewt5WkThuBKvigpkW4Iw6uQn2vwSFn4Gto1H4cdrDEBDTaa8TCRGThIgQYrA5XVbDusxC\nNmYVU9PQQoSfK6viglk22R+36lyj18nht6CpBoaPN5YJj08GR9f215AQMUmICCEGq7qmFt4/fI51\nmYUcLanG2cGWpZP8WRUXxFgvGziaZlQTLj0GjkNhwr3GtRPfaAmRqyREhBACDhdVGb1OjpyjobmN\nyUHDWBUbzMLxw3G6cNC4EH/8bWhtguA7UA9mSIiAhIgQQnRUXdfMhqxi1u8t4Gz5FYY525M8NYAV\nscGEDmmAQ+vgwMuoHx2REAEJESGE6IzWms/OXGTd3gI+OF5KS5tmRrg3K2ODmRvljb29fa+FSL/v\nJyKEEKJrSiluH+3N7aO9KbvcwJv7jV4na9YdZPjQ3m3X2707VYQQQlgF36FO/HN8OB//+2xe/NZU\nIoa79er7SYgIIcQAZGdrw/yxw3n1wem9+j4SIkIIIXpMQkQIIUSPSYgIIYToMQkRIYQQPSYhIoQQ\nosckRIQQQvSYhIgQQogekxARQgjRY1ZVO0spVQPkWHoc3eANVFh6EN0g47x1rGGMIOO81axlnJFa\n6165dd3aamfl9FYRsVtJKXVAxnnrWMM4rWGMIOO81axpnL312nI6SwghRI9JiAghhOgxawuRFy09\ngG6Scd5a1jBOaxgjyDhvtUE/Tqu6sC6EEKJ/sbaZiBBCiH5kQIWIUipEKXXM0uMQQghrp5R6Uin1\n2M3261chogz9akxicOvOFxOl1BtKqSNKqX/pw3H9SCnl3FfvJ6yLUqrPbt+w+C9s8x9pjlLqVeAY\n8C2l1GdKqSylVJpSytXc7+dKqf1KqWNKqReVUsrcPlUpdVgpdRj4vgU/iuiHevuLiVJqODBNaz1B\na/2/3TzmVvwD/xEgITJImL8nTyqlXlFK5Sql1iul5iql9iilTimlppszh9eUUnuA15RStkqp/zF/\nbx5RSn3XfC1XpdQO83fsUaXU0g7v81Pz9T8BIrs1OK21Rf8AIUAbEIdx9+duwMV87nHg5+bfPTsc\n8xqw2Pz7EeAu8+//Axyz9GeSP/3iZyoHeBU4DtwPfAZkAWmAq7nfz4H9GF9eXuTaQpOpwGHzT5c/\nU+bPXz1wCJgBTAIyze1vAx7mfruA54ADwI+BV4A/mfueBWYBLwPZwCsdXv9P5jHHgV+a2/4ZaAKO\nAh9a+v+3/Omzn+kWYDzGl/+D5s+LApYC7wBPmtuHmMc8AvzM/Luj+XMUinGT+VBzuzdw2nydqebP\nlDMw1Nz+2M3GZvGZiKlAa52JESRjgD1KqUMY//iDzX1mK6X2KqWOAnOAsUqpYcAwrfVuc5/X+nrg\not8KB/4PmAk8BMzVWk/B+If0r+Y+z2utp2mtxwFDgEXm9r8Bj2qtJ3bjfZYAZ7TWk7TWH2ME1+Na\n6wkY/yB/0WFfB611jNb6t+ZjD+A24F+A94D/BcYC45VSk8x9fqqNO6InADOVUhO01r8HzgGztdaz\nv8r/FGHV8rTWR7XWbRhfKnZoIwmOYoQMwHta63rz7/OBb5u/S/cCXhj/LhTwa6XUEWA74A/4YXwJ\neltrXae1vozxM3lT/aXsyRXzvwrYprW+r+OTSiknjF8IMVrrIqXUk4BT3w5RWJkCrXWmUmoR176Y\nADhgzErA+GLy7xjfvDyB40qpj/nyF5PE7ryhUsrdPPYjc9PfMWY+V7113SHva621+cWoVGt91Hyd\n4xi/FA4BKUqpRzD+rY4wP8uR7oxHDDiNHf7e1uFxG9d+l1/psI/C+DK0teOLKKVWAz7AVK11s1Iq\nn6/x+7S/zESuygTuUEqNBlBKuSilIrj2ASvMayTLAbTWVUCVUupO8/mVfT1g0W9d/8VkkvlnjNb6\noQ5fTJZrrccDf6H3v5hcue5xx18C1/+CsFNKhQKPAfHmzGZzH4xRDBxbge8ppewBlFIRSikXwB0o\nMwNkNtfO9uwG7lZKDVFKuQGLu/Mm/SpEtNblwGrgDXOq9RkQZYbFXzDOXW/FOI991QPAH80pm+rb\nEQsr0GdfTLTW1UClUmqGuelbwEddHHIzQzGCp1op5ccXZ0Q1QK9UZRUDxkvACSDLXGH4Z4wZy3og\nxpwBfxs4CaC1zsKYLR8GMvji79kbsvjpLK11PjCuw+OdwLRO9vsZ8LNOth8EOp67/vdbP0phrbTW\n5eb0/Q2llKO5+Wda61yl1NUvJhf48heTl5VSGvjgK77l/cAL5vLbs+Zr9XTsh5VSn2P8Iy8C9nR4\n+iNt5i0AAAB6SURBVEVgi1LqnFwXGfg6+T25+kbPddjeBjxh/rnebTd4n6eAp77K2KTsiRBCiB7r\nV6ezhBBCWBeLn84SwhoopRYAz1y3OU9rvcwS4xGiv5DTWUIIIXpMTmcJIYToMQkRIYQQPSYhIoQQ\nosckRIQQQvSYhIgQQoge+/+F6Ae89boe0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc3f2f6198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot results\n",
    "df.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"2\" halign=\"left\">read</th>\n",
       "      <th colspan=\"2\" halign=\"left\">read_format</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>ratio</th>\n",
       "      <th>difference</th>\n",
       "      <th>ratio</th>\n",
       "      <th>difference</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>small</th>\n",
       "      <td>8.31948</td>\n",
       "      <td>0.00355493</td>\n",
       "      <td>5.21809</td>\n",
       "      <td>0.00204864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>normal</th>\n",
       "      <td>7.34758</td>\n",
       "      <td>0.00362237</td>\n",
       "      <td>4.59613</td>\n",
       "      <td>0.0020522</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>large</th>\n",
       "      <td>2.43689</td>\n",
       "      <td>0.00376415</td>\n",
       "      <td>1.77478</td>\n",
       "      <td>0.00202965</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           read             read_format            \n",
       "          ratio  difference       ratio  difference\n",
       "small   8.31948  0.00355493     5.21809  0.00204864\n",
       "normal  7.34758  0.00362237     4.59613   0.0020522\n",
       "large   2.43689  0.00376415     1.77478  0.00202965"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# df comparing ratios and differences of read and read_format to mread\n",
    "# init dataframe\n",
    "tups = itertools.product(('read', 'read_format'), ('ratio', 'difference'))\n",
    "cols = pd.MultiIndex.from_tuples(tuple(tups))\n",
    "df2 = pd.DataFrame(index=list(PATHS), columns=cols)\n",
    "\n",
    "# fill in values\n",
    "funcs = {'ratio': lambda x, y: x / y, 'difference': lambda x, y: x - y}\n",
    "for read_name in ('read', 'read_format'):\n",
    "    for func_name, func in funcs.items():\n",
    "        df2.loc[:, (read_name, func_name)] = func(df.loc[read_name], df.loc['mread'])\n",
    "df2.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "It is somewhat understandable that calling read without any format specified may have some significant overhead as the function has to try to determine what kind of file it is reading. However, when I specify the file type there is still a 2 ms + overhead compared to calling _read_mseed directly. Getting rid of this overhead can mean significant speed ups when batching through mseed files, but calling a the private method _read_mseed outside of obspy seemed dirty (respect privacy)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Profiling\n",
    "Next, I profiled a simple script with the python profile tools and visualized the results with [snakeviz]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " "
     ]
    },
    {
     "data": {
      "text/plain": [
       "3 Trace(s) in Stream:\n",
       "BW.RJOB..EHZ | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples\n",
       "BW.RJOB..EHN | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples\n",
       "BW.RJOB..EHE | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%prun\n",
    "read('temp.mseed', format='mseed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext snakeviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file '/tmp/tmpxl5fn4nj'. \n"
     ]
    }
   ],
   "source": [
    "%snakeviz _ = read('temp.mseed', format='mseed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the sunburst visualization (run the notebook if you want to see it) I found the _read_mseed was only taking about 17 % of the time whereas load_entry_point from the pkg_resources module was taking over 60%. Obspy uses the setuptools plugin system to find various IO methods, but in this common case it has significant overhead.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caching (memoization)\n",
    "A reasonable speed up might be cache all results of load_entry_points so function need only be loaded once rather than every call. This could potentially have huge benefits for the many obspy functions that depend on load_entry_point. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# basic decorator for load_entry_points\n",
    "import functools\n",
    "import inspect\n",
    "import pkg_resources\n",
    "\n",
    "\n",
    "def decorate_load_entry_point(func):\n",
    "    cache = {}\n",
    "    sig = inspect.signature(func)\n",
    "    \n",
    "    @functools.wraps(func)\n",
    "    def wrapper(*args, **kwargs):\n",
    "        b = sig.bind(*args, **kwargs).args\n",
    "        if not b in cache:\n",
    "            cache[b] = func(*args, **kwargs)\n",
    "        return cache[b]\n",
    "    \n",
    "    return wrapper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now we can test if this will actually give us a speed up by monkey patching load_entry_points function in the obspy.core.util.base module. This is generally bad practice but it is a good way to see if it will work, and thus merit refactoring some of obspy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# monkey patch\n",
    "import obspy.core.util.base\n",
    "func = decorate_load_entry_point(pkg_resources.load_entry_point)\n",
    "obspy.core.util.base.load_entry_point = func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 10.46 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "1000 loops, best of 3: 1.03 ms per loop\n",
      "1000 loops, best of 3: 971 µs per loop\n",
      "1000 loops, best of 3: 490 µs per loop\n",
      "1000 loops, best of 3: 1.14 ms per loop\n",
      "1000 loops, best of 3: 1.11 ms per loop\n",
      "1000 loops, best of 3: 591 µs per loop\n",
      "100 loops, best of 3: 3.27 ms per loop\n",
      "100 loops, best of 3: 3.27 ms per loop\n",
      "100 loops, best of 3: 2.69 ms per loop\n"
     ]
    }
   ],
   "source": [
    "# profile each stream type with each function\n",
    "df = pd.DataFrame(columns=list(PATHS), index=list(READ_FUNCS))\n",
    "for stream_type, path in PATHS.items():\n",
    "    for func_name, func in READ_FUNCS.items():\n",
    "        time = %timeit -o func(path)\n",
    "        df.loc[func_name, stream_type] = time.best\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>small</th>\n",
       "      <th>normal</th>\n",
       "      <th>large</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>read</th>\n",
       "      <td>0.00102862</td>\n",
       "      <td>0.00114128</td>\n",
       "      <td>0.00327326</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>read_format</th>\n",
       "      <td>0.000970601</td>\n",
       "      <td>0.00110734</td>\n",
       "      <td>0.00327072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mread</th>\n",
       "      <td>0.000490304</td>\n",
       "      <td>0.00059149</td>\n",
       "      <td>0.00268998</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   small      normal       large\n",
       "read          0.00102862  0.00114128  0.00327326\n",
       "read_format  0.000970601  0.00110734  0.00327072\n",
       "mread        0.000490304  0.00059149  0.00268998"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like this helped quite a bit, speeding up the read_format function from 2.6 ms to 1.1 ms, but it still seems that read_format shouldn't be so much slower than mread. Time to profile again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      "*** Profile stats marshalled to file '/tmp/tmpwz24753k'. \n"
     ]
    }
   ],
   "source": [
    "%snakeviz _ = read('temp.mseed', format='mseed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cache works well; only 2% of the time is spent in the loading the plugin (after the first load) and we could probably make the cache more efficient but that might be over optimizing. \n",
    "\n",
    "Now we can see about 55 % of the time is spent on the _read function, and a large chunk of time (~30%) is spent in the tarfile.is_tarfile function, and about %5 is spent in the zipfile.is_zipfile function. A reasonable compromise might be to add something like \"check_compression\" as a parameter in read that can be switched off to avoid these function calls when one is sure the files are not compressed with tar or zip formats. After doing so the read_format function should be fairly close to the mread function speeds. Let's simulate this by monkey patching the tarfile module and the zipfile module to return false whenever is_tarfile or is_zipfile is called. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import tarfile\n",
    "import zipfile\n",
    "\n",
    "tarfile.is_tarfile = lambda x: False\n",
    "zipfile.is_zipfile = lambda x: False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 5.00 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "1000 loops, best of 3: 655 µs per loop\n",
      "1000 loops, best of 3: 615 µs per loop\n",
      "1000 loops, best of 3: 480 µs per loop\n",
      "1000 loops, best of 3: 737 µs per loop\n",
      "1000 loops, best of 3: 701 µs per loop\n",
      "1000 loops, best of 3: 568 µs per loop\n",
      "100 loops, best of 3: 2.71 ms per loop\n",
      "100 loops, best of 3: 2.67 ms per loop\n",
      "100 loops, best of 3: 2.54 ms per loop\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>small</th>\n",
       "      <th>normal</th>\n",
       "      <th>large</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>read</th>\n",
       "      <td>0.000654831</td>\n",
       "      <td>0.000737354</td>\n",
       "      <td>0.00270822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>read_format</th>\n",
       "      <td>0.000614946</td>\n",
       "      <td>0.000700654</td>\n",
       "      <td>0.00266961</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mread</th>\n",
       "      <td>0.000479516</td>\n",
       "      <td>0.000567598</td>\n",
       "      <td>0.00254124</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   small       normal       large\n",
       "read         0.000654831  0.000737354  0.00270822\n",
       "read_format  0.000614946  0.000700654  0.00266961\n",
       "mread        0.000479516  0.000567598  0.00254124"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# profile each stream type with each function\n",
    "df = pd.DataFrame(columns=list(PATHS), index=list(READ_FUNCS))\n",
    "for stream_type, path in PATHS.items():\n",
    "    for func_name, func in READ_FUNCS.items():\n",
    "        time = %timeit -o func(path)\n",
    "        df.loc[func_name, stream_type] = time.best\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could do some further optimizations of both the read_format and mread functions but it seems fairly reasonable now. We have achieved a speed up from 2.6 ms to 0.7 ms on the normal stream with minimal changes to the obspy code base."
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}