{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Fast Lomb-Scargle Periodograms in Python\n", "\n", "The Lomb-Scargle Periodogram is a well-known method of finding periodicity in irregularly-sampled time-series data.\n", "The common implementation of the periodogram is relatively slow: for $N$ data points, a frequency grid of $\\sim N$ frequencies is required and the computation scales as $O[N^2]$.\n", "In a 1989 paper, Press and Rybicki presented a faster technique which makes use of fast Fourier transforms to reduce this cost to $O[N\\log N]$ on a regular frequency grid.\n", "The ``gatspy`` package implement this in the ``LombScargleFast`` object, which we'll explore below.\n", "\n", "But first, we'll motivate *why* this algorithm is needed at all.\n", "We'll start this notebook with some standard imports:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# use seaborn's default plotting styles for matplotlib\n", "import seaborn; seaborn.set()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To begin, let's make a function which will create $N$ noisy, irregularly-spaced data points containing a periodic signal, and plot one realization of that data:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAFVCAYAAAApGgzgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmQlNW9//FPzwrMNDiasfz9kgKVawgEMYxaIShuVRqj\nxDIw6IAZCHJTiUtI4kAFr4IblhOR8sYAIWp0cGI5QmmqEqzKvRr1kgylISBEIBovRTCpbCPy055h\nmYbu3x+ke7bunl6e5Zzneb+qqOqF7j5z+unne5bvOU8kmUwmBQAArFDmdwEAAED+CNwAAFiEwA0A\ngEUI3AAAWITADQCARQjcAABYpKTAvWvXLjU3Nw95vK2tTTNnzlRzc7Oam5u1f//+Uj4GAAD8S0Wx\nL3ziiSf085//XDU1NUOe27Nnjx5++GFNmjSppMIBAICBiu5xjxs3TmvWrFGm/Vv27Nmj9evXa968\neXr88cdLKiAAAOhTdOC+6qqrVF5envG5a6+9Vvfff782bNig7du36/XXXy/2YwAAQD+uJKctWLBA\np5xyiiorK3XppZdq7969Of8/u64CAJCfoue4s4nFYrruuuv00ksvaeTIkXrjjTfU2NiY8zWRSERd\nXTGni4J+6uuj1LEHqGf3Ucfuo47dV18fLfq1JQfuSCQiSdq8ebMOHz6sG264QS0tLZo/f76qqqo0\nffp0XXLJJaV+DAAAkBQx5epgtO7cRQvaG9Sz+6hj91HH7iulx80GLAAAWITADQCARQjcAABYhMAN\nAIBFCNwAAFiEwA0AgEUI3AAAWITADQCARQjcAABYhMANAIBFCNwAAFiEwA0AgEUI3AAAWITADQCA\nRQjcAABYhMANAIBFCNwAAFgkkIF76bqtWrpuq9/FAADAcYEM3AAABBWBGzAIo0UAhkPgBgDAIgRu\nAAAsQuAGAMAiBG4AACxC4AYAwCKBC9yPdLylgx8f1cGPj+qRjrf8Lg4AAI4KVOB+pOMt7f3TofT9\nvX86pJa1nTrw95iPpQLyQ6MTQD4CFbj/0C9opxyKHdNjL/zeh9IA+cvW6PzuD3/Dum4AAwQqcAO2\nytbojB3u9aE04cPGN7BJoAL3xDPrhjxWF63W4tlTfCgNAADOC1TgXtI0VXXR6vT9umi1Vt92kcad\nEfWxVMDwsjU6o6OqfCgNAJMFKnBL0uLZU1QWkcoioqcNa2RrdFaUB+4nCqBEgTsrjDsjqrroCNVF\nR9DThlVodCKIyB9wXoXfBQBwUqrRmboNAJkErscNBAXrugFkQuAGDMRmQt6hgQTbWB+4mT8xG99P\ncdhMyBs0kGAj6wM3ABSLBhJsFMjktFW3Tve7CEBJJp5ZN6AnKLGZEICT6HGjKAyBu4vNhLzBbovu\nIn/AHQRuwBCphlBqxIh13e6jgeQe8gfcQ+AGDMVmQt6ggeQO8gfcE8g5bgDIFxvfwDYl9bh37dql\n5ubmIY+/+uqramxsVFNTkzZt2lTKR+TE/AkAmIn8AfcUHbifeOIJ3X333YrH4wMej8fjam1t1dNP\nP6329nY9//zzOnjwYMkFHYz5E38dih3VodjRnP+HhhUQXuQPuKfowD1u3DitWbNGyWRywOP79u3T\n2LFjFY1GVVlZqfPPP1/btm0ruaCDMX9iNhpWAMgfcEfRc9xXXXWV/vKXvwx5vLu7W9FoX4uqpqZG\nsRgn67BIZUZ/+PHQ3niqYbX6tou8LhYAH5A/4A7Hk9Oi0ah6enrS93t6ejRmzJhhX1dfX9iXet45\n9dr5XteAx04bM0J33/z5gt8rLAqpl0Ur/1uS9JO7rxry3PL1W5X410DLYy+8rQe+2bfhTXl55OSN\niKTkkJeqrCwS+O+nmL9v+fqtOvivxk7/Ok3VZ9DrrFBO1wf1PJRTdUHdOs/xwH322WfrwIED+uij\njzRy5Eht27ZNixYtGvZ1XV2F9coXzz5XLWs7dSh2TNLJ+ZNVt0wv6r3CoL4+WlC9nDhxMuoOfs3g\nIfCd73Vp/r2/1OLZUzTujGj6dRPHZd756/avnBvo76fQepZy12nrN74gqe97GLzWO4yKqePh9D/e\nqWNn6zjbuSTsSmnIlLyOOxI52ZravHmzNm7cqMrKSi1btkyLFi1SU1OTGhsbdfrpp5f6MRkVM3/C\njl+lyTe3gMSU/JGv4b9Vt04PdaCGXUrqcX/qU59SR0eHJGnmzJnpxy+//HJdfvnlpZUsD5nmT2gt\nm2Px7Cl6YMO29G24i2MfCAd2TkNBClmbyc5f+WG9K4BCELhRkCVNU1WRSkCTVFEeYQi8REwrIMiY\nhnAegRsFeaTjLR0/0ZcufvxEkvXZDmC9K4B8EbgxRK4dz3IlUrFTWvGYVvAWSaqwGYEbAxS741ns\ncG/G1x0/kXCtrGFDwwhBka3hRIMqP9YHbuZPnDXc0qRsiVT9h8/7vy52uNf5QoYQW8h6g8YRbGB9\n4O5vuB8dP8rSZUukimT5/2Nqq2lYOWC4BhXHduloHPmHnnZhAhO4h/vR8aPMTz5LkzLdZkmTfzi2\nncFGOLBFYAL3cD+6fH+UYW/55bM0adwZ0XQGdOpxljS5K1fDiIADhEtgAjeck8/SpFQGdKGvQ3Fo\nGLmvkFGjsDfw4a/ABO7hfnQM5eav2KVJLGkqzXCJltkaRhzbhcmWD0DjCLYITOAe7kfHjxK2y9Yw\n4tjO33D5AIwauS9bw+njnl4SLPMUmMAtDf+jG/w8w10ICgJOfobLB2DUyF3ZGk4fdR9TvN+eDyRY\n5haowD3cj44fJYKKYxs2yNZwOp7IvA8ECZaZlXRZT4QXa7Nho4ln1g3o8UnkA8A+gepxl4pNLIDg\nyDQVRj6AvzIlUmZDgyo7Ave/sImFM9iCFqYrNR+ABn7xBjecsolERIMqBwL3v7CJxUAEYARVKfkA\nNPBL17/hdPb/GZ3x/3x95iSPS2UX5riRlhpWJGCbi+/GfbnqOFcDf/VtF7lZrMBINZwk6e4FF6hl\nbacOxY5JOhnM66IjNO2zZ/hZROMFLnAPd2JLPZ8a7krdJmkFtiOow0aLZ0/RAxu2SZKio6p8Lo0d\nAhe485FpuKsuWq3oqErFDscl9SWtAEAKDXzn9e+BIz+hnOPONtwliU0sAGRFVjpMEMrAnU1FeRmb\nWADIiV3q4LdQBm4uygAE23BLtkpZNcEudfBbKOe4lzRNHZDJyHw2EBzZlmwtnj2l5EDLtQ2cka3R\nRIJlfkLZ45YY7gKCaOm6rUOSx6Rw78mA4Alt4Ga4ayB2gwIAO4RyqBwDuTm0CHitsrxswCUiJXJY\nTMcQeWEI3IOE8QBiNygEyeiaKiWSSXJYLMdOjtkRuAEETv/duNzoaRNM4KfQznGjD8vjgiPTpSzD\niBwWBBmBG+wGBeSBBM6BaCT6h8ANSSyPg/3cDKxczhMmCfUcN/NUffpv9E9PG7bJFFjLIs5dbYoE\nTpjEmh43wzLBw3cKp2QKrImkFDvc60NpAHdZE7gBwC8kcHqLfILcCNwArJctsC5fcKEj708Cp3fI\nJxgegRsIiDD3UrwIrCRweiNXPgFOInADAUAvxf3AytrwPv0bid9Y9Tq5Kh4jcCOtlGsUw1/0Ugis\nXhncSIyfSOhQ7KhjjUTyCYZH4AZcRva8e6hb72XL4HeqkUg+wfCsCNxhnrsLKr5TZ9FLQZCQT5Cb\n8YGbubvg4Tt1Hr0UeCVTI9HpAMu0R25FBe5EIqEVK1aoqalJzc3Nev/99wc839bWppkzZ6q5uVnN\nzc3av39/0QVk7i54+E7dQS8FXhjcSCyLiADrsaK2PH3llVcUj8fV0dGhXbt2qbW1VevWrUs/v2fP\nHj388MOaNGmSYwUFkBvb1rqP5M2T+l821altZZG/onrcO3bs0IwZMyRJ5513nnbv3j3g+T179mj9\n+vWaN2+eHn/88ZIKaMLcHQkwzjLhO0UwsTLCG6lGYnlZmT7q6SVXxWNFBe7u7m7V1tam75eXlyuR\nSKTvX3vttbr//vu1YcMGbd++Xa+//nrRBWTuLnj4TgH7fdzTq/iJvvM+uSreKWqovLa2Vj09Pen7\niURCZWV9bYAFCxakA/ull16qvXv36rLLLsv5nvX12U/aK/59mlr+83/St3P9XzeUl0ck5S6jDUwq\nv9/fqZsG/y1eHj9BOVaHU18f1fL1W3Xw46OSpMdeeFsPfJOetpOGO4b6B+2UQ7FjWvOzt9W24osl\nf35YjuViFBW4Gxoa9Nprr+lLX/qSdu7cqQkTJqSfi8Viuu666/TSSy9p5MiReuONN9TY2Djse3Z1\nZW+ljakuT8/djakuz/l/3XDiRFJS7jKarr4+alT5/f5O3TK4nh/peEv/PHREkvS9H27Rkqaprn5+\nEI7V4dTXR/W9H24ZsDJh53tdmn/vL7V49hRGbhxQyvkikUg6cvy1fuMLkoJ7LJfSIClqqPzKK69U\nVVWVmpqa1NraqjvvvFObN2/Wxo0bFY1G1dLSovnz5+umm27Spz/9aV1yySVFFxCwFcve3MPKBP9V\nlg8NH+SqeKOoHnckEtF999034LGzzjorfXvmzJmaOXNmaSUDLJcruKy+7SJXPpPELHhldE2VDsWO\nKnFykCedqwL3Gb8BCwAMxsoE/626dbqWL7jQt70Dwrzax5jAbeqXwNacKBbBxT2sTDADO5z5w5jA\nbSLmKFEKgou72CkOYVXUHLcf/Ji782OOMkzCMB/bf4cpgouzTNwpLjVqGIZjG/6xJnADNjIxuACw\nG0PlOTBHCQAwDYE7B+YoAQCmIXAPgwQYAMjOjwu7hH21jxGBO7XnsIlfAssdEBSmLrkECsFqH0MC\n9873utK3w/glAADyw3a3hmaVs+QKQcLSIPeYVLep4dvUbbcvKIPwMqLHDX8wdAo4IwzDt6acL1jt\nY2jgDtuXAMBuDN96h9U+hgTu08aMSN8O45cAAMhf2Ff7GBG4777586H+EhAspgwpwjsM33or7Kt9\njEhO+7dPnWL0tpAmJcAEGfs8w1ZLmqaqZW2nDsWOSeLa1HCXET3uweixAM4K+4YVXgj78C28Y2Tg\nBuCcMGQ8myDsw7fwDoE7pOiBhQcZz0CwELiLYPtQPj0wAIWgoW8WAncI0QMLFzKeUQpTG/p+XNzE\nFMYE7jB/CYCb2LACpaChbx5jAje8Y2IPLChDcaZe6Y6MZyA4jAvcQTmBm8y0HpipQ3GFeqTjLWOv\ndEfGszeCOHJoYkM/7IwK3EE5gdvApB5YUIbigvJ3AP2Z1tCHYYGbE5936IEByJdJDX0YFrgRTkEZ\nigvK3wEMRkPfLEYFbhtOfMzBOy8oQ3FLmqZypTsArjMqcJt+AmcO3j1BGYrjSncA3GZU4JbMPoEz\nB++eoAzFpa50Z+Lf0T/j2fbd/4AwM+Kynv2lTuCp2wAAoI9xgdtkE8+sGzBULpk3B1+IoK03BYAw\nIHAXYEnTVLWs7dSh2DFJfXPwABB0NPTNYdwct+lMnoMHAAQfPe4CMQcPAPATgRtwGEOKGCyVwc+x\nAScYGbg5uMOJ790bqU2EUreXNE31uUQACsEcNxAibCIE2I/ADYQImwgB9iNwAwBgESPnuE3HXCxs\nFbRNhIAwoscNhIjpF/IBMDwCNxAybCLkLS4FDKcRuIGQCcqV2GxAFj/cUFTgTiQSWrFihZqamtTc\n3Kz3339/wPOvvvqqGhsb1dTUpE2bNjlSUACwDVn8cENRgfuVV15RPB5XR0eHlixZotbW1vRz8Xhc\nra2tevrpp9Xe3q7nn39eBw8edKzAAACEWVGBe8eOHZoxY4Yk6bzzztPu3bvTz+3bt09jx45VNBpV\nZWWlzj//fG3bts2Z0iI0lq7bmt4mErDVxDPrhjxGFr+zwniuKGo5WHd3t2pra9P3y8vLlUgkVFZW\npu7ubkWjffNmNTU1isWGn8+pr2euzW021XF5eUSSXWVOsaHMNtevZE+5v/+tS/S1+/9LBz86ucXs\naWNGqG3FF30uVX5sqWPbj+ViFBW4a2tr1dPTk76fCtqSFI1GBzzX09OjMWPGDPueXV0ka7ipvj5q\nVR2fOJGUZN9xYUs921q/kj11nHL7V87VAxu2pW/bUHab6tjWY7mUhkZRQ+UNDQ3asmWLJGnnzp2a\nMGFC+rmzzz5bBw4c0EcffaTe3l5t27ZNn/vc54ouIABI9g6JksUPpxXV477yyivV2dmppqYmSdJD\nDz2kzZs36/Dhw7rhhhu0bNkyLVq0SIlEQo2NjTr99NMdLTSA0rD7H2CvogJ3JBLRfffdN+Cxs846\nK3378ssv1+WXX15ayQAACAgnr8nOBiwAAEn2TkeEDYEbAACLELhhHPZ2BpCPsJ4rCNwwCns7+4uh\nUnesunU6CYEOC/O5gsANo7C3MzLp37Navp6GBcJ9riBwAzDa4J7Vzve6QtOzAjIhcAeUrUOe7O2M\nwcLcs0J2YT5XELhhlCVNU1UXrU7fr4tWa/VtF7HjFOAy26YjbDpXOJ1ER+CGcRbPnqKyiFQWUSha\nz8gtzD0rr9g6HWHDuSJbEt3//uX/Ff2eBG4Yh72d0d/gntVpY0YY27Oyla3TETacK7LV7cqn3iz6\nPQncQIlszScYzOQ1sf17Vnff/Hm/iwP4isANwPg1sf17Vv/2qVP8Lk7gMB3hnmx1W0oDlMANwNqh\nUjiD6Qj3ZEuiK6UBSuAOIJOHPAGYiekI9zidREfgDhg3MhgRfAyVgukI9zidRFfU9bhhrr05MhhX\n3WLPXsmZ9nV28nq2GGhJ01S1rO3UodgxSX3DeXAex7GzwliP9LgBSLJjTSycF5RVEWFCjztgKsvL\nFD+RGPBYqRmMCIfUcF7qtmnC2LMCMqHHHTCja6pUFum770QGI7IjERCA1wjcARQdVcWQpwcyJQJ+\n7f7/MmbtM4BgYqjcQaYknVSUlxk95BkUmdY+H/zoqB574fckdsFKfp+7gszJuqXHDQCARQjcAdJ/\nvvXjnl6/i+MoE+eSM619Pm3MCKYnkJWJxzGc53amPoE7IAbPt8ZPJHQodjQQ862m7qOdaSvDthVf\ntHp6YtWt0xkudYmJxzENCTsRuAMi03xrIqlA7DVt8j7arH1Gvkw7jk1sSCA/BG6gBDZcD9hWbAzi\nLtMaEm4I6jFE4A6IIO81HeS/DeHBcQynELgDItul44LQCwzy34bwMO04piFhLwK3Q0xI8gjyfGuQ\n/zaEh0nHsWkNCeSPwO0AU5I8gjzfGuS/DeFh2nFsUkMiKLzoxBG4HRCGJA8AwWNaQ8J2XnXiCNxA\nEfpnq7L2GYDkXSeOwO0AkjwAZ5mQMwK7BfkYInA7gCQPwDmm5IzAXn4dQ1514gjcDiHJA3AGOSMo\nlV/HkFedOC7r6ZBUkkfqtl+CPNca5L8N4cFxHGyLZ0/RAxu2pW+7gcANwCgTz6wbMMwpkTPipiA2\nJPw8hrzoxDFUHgBB3Y8X4UTOCEoV9GOIwA0UKMjZqqYgZwSlCvIxROAGCkDGszfYGASlCvIxROAG\nCkDGMwC/kZzmoCAmeQAAzEKPGygAu+QB8FvBPe6jR49q6dKl+vDDD1VTU6PW1ladeuqpA/7PypUr\ntWPHDtXU1CgSiWjdunWqra11rNCAX5Y0TVXL2k4dih2T1JetCgApbo++Ftzjfu655zRhwgQ9++yz\nuv766/WjH/1oyP/Zu3evnnrqKbW3t+uZZ54haLuIDGf3DV5uF+RsVQDmK7jHvWPHDn3961+XJM2Y\nMUPr1q0b8HwikdCBAwe0fPlyffDBB2psbNTs2bOdKS0GyJbhvHj2lMBlUQ7W/8pcXjNll7ygC1vO\niJ/HdFAFtS5zBu5NmzbpmWeeGfDYaaedppqaGklSTU2NYrGBy2COHDmi5uZmLVy4UMePH9f8+fM1\nefJkTZgwIWdB6us5ARbqDwcyZziv+dnbalvxxSHPBamOy8sjkrz5mzJ9Vq7PD1I9myqIdezlMZ0P\nU8qBoXIG7jlz5mjOnDkDHvvWt76lnp4eSVJPT49Gjx494PmRI0equblZ1dXVqq6u1rRp0/TOO+8M\nG7i7ulgHW7Bk5ocTieSQ+qyvjwaqjk+cOPnHe/E3ZfqsbJ8ftHo2UVDr2MtjejhBrGPTRjRKaRgV\nPMfd0NCgLVu2SJK2bNmiCy64YMDz+/fv17x585RIJBSPx7V9+3ZNnjy56AIiOzKcAcAMXm49XfAc\n99y5c/W9731P8+bNU1VVlVavXi1Jamtr09ixY3XFFVfo+uuv14033qiKigrNmjVL48ePd7zgtnGj\ntUeGM9xmWi8FQBGBe8SIEfrBD34w5PGvfe1r6dsLFy7UwoULSyoY8uPFJeQAAOZg5zTLkeHsrtRy\nu9TtJU1TJdEDBeAfdk6Ddbxau84FReAVr45pLgEcDARuWMXLYMoFReAFGogoFIEbViGYImg4pt0X\ntB0mCdwIhI+6jzn+niy3A+znxYiG1w0DArcHgtba81OmYFoWkaKjqhz/rCVNU1UXrU7fTy23C0sS\nIMetN2ggusvtEQ0/pjoI3C7z4ktddev00GQ5ZwqmddERqih351AO6wVFmHf1TtgbiLbzY6qDwO0y\n5q+c52UwTS23q4uOCNWJlOPWW2FtIHohiCMaBG5YJ6zBFMHFMe0et0c0/GgYELhdtHTd1oxDuLa3\n9hB8QeylhF2YcxbcHNHwY6qDwO2y0TVVzF/BOsy7BkvYcxbcHtHweqqDwO0B5q/cE+ZehNu8PG7Z\n0ctd5Cy4y+upDvYq9wD7ibvj457e9D7iUl8vYvHsKY7Uc9gDCcctYCZ63LBW/ERiyGP0IoChyFkI\nFnrcsNKqW6drUeurSnr0WYDb3DzOljRNVcvaTh2KndxhMJWzADvR44a16EUA+SPXJjgI3C5xMmmK\nxJ3MyHxG0Lj5W2eteHAwVO6CTEsv3NpPO+wWz56iBzZsS98GgEyCNOVF4HZBpqUXiaRUVhbxoTTB\nRuazu/qf7FI9QadPgKnRqdTtJU1THX1/wAteNgwYKgcyYH24N8K+MQhQDAK3C9xMmmK+230EE+dl\nO27ZGAQoHIHbBSRN2ad/YCGYADAZgdslLL0AhseSPm+tunV6oJK0worA7RKWXtiLYOIdRqf6kFdh\nBxOmKwnchuPHPDynexEEE28xOkVeBQpD4DZYph/zodhRHc+wRzecRTAZyM0GJKNT5FWgMARug2Vb\nDx473OtDacKFYNLHrd6gCUOOgI0I3BZKJMWwOTxTam+Q6Z7huZFXQcOocLbUGYHbIIMPmkw/5hTm\nwJxDYHEPc7f5Ia8ChSBwu6jUpKnBP+bBmAMrTv8GEoFleKX0Bpm7zV8heRW29AzhDgK34VI/ZriD\nwDI8eoPeIK8C+SJwGy71Y64sH/pVsbYYXik2yz7f3jobg8AGpkyrcXUwS4yuqVIimdSh2DFJfb0e\nlGbimXUDhsqlvsBCr6dPrquw5bpq2JKmqWpZ28lxC+tlm1bz41xBj9sirC12HsPA7st03JrScwHy\nZdK0GoHbIsyBuYMGkbsGH7ckBLqPhlH+Uol+NtUZgdtHg7ObbTlogoYGkbdM6rkEEQ2jwn3c0zts\nnZl0DQMCtwH4oQGQ8kvSG66RT8OocPEM20gPrjOTptUI3AYY7odGxm1xMq11ZWTDfyb1XGxDI99f\npkyrEbgRGpz0zGBSz8U2+fSmaRgVLtdy2/4dAFOm1QjcBuCH5g2GEEvj5MiPKT2XIKJhVLjRNVVW\n1RmB2wCF/NAYNodJ8p16GHzcmtJzsU2+jXwaRoWzqc4I3Iaw6aCxFSMbzmLqwXv5NvJpGOWnf8Nz\n0+v/a02dEbgNwQ/NfblOeoxkFI6pB39kauRz0ZHCZWp4Hood1fEMGeamKTpwv/zyy2ppacn43MaN\nGzV79mzdeOONev3114v9iEAju9kfjGx4g0DivFSd0sh3RqaGZyIpxQ73+lCawhS1V/nKlSvV2dmp\nSZMmDXmuq6tL7e3tevHFF3Xs2DHNnTtX06dPV1VVVcmFDYpsQ4yJRFIVGbIbUbhUwyh1e0nTVEm5\n99xGYXLt877mxbd9KhXgLhNG5oqKEg0NDbr33nuVTCaHPPf73/9eDQ0NqqysVG1trcaNG6d33323\n5IIGSbYhRhtaejZg7tUbZC/DZplyXsoiUnRUXyfT1JHRnIF706ZN+vKXvzzg3+7du3XNNddkfU1P\nT4+i0b4fbk1Njbq7u50rMTAM5l69w9QDbJWp4VkXHZEe9TS5A5BzqHzOnDmaM2dOQW9YW1urnp6e\n9P2enh6NHj162NfV14enlX7eOfXa+V7XgMdOGzNCiURS5eUR1+oiNHUckTR0MEhlZSfrtrw8Ism9\n+ghNPevk3/qJU0ZKki449/+mH8+njtvu+WJJnxs2/et0cP1mqu9Sj/Mw1PGKf5+mlv/8n/Tth9p+\nK+nk3/6HA5k7AGt+9rbaVhR/7DrB8etxT5kyRY8++qh6e3t17Ngx7du3T+ecc86wr+vq8r8V45XF\ns88dco3iVbdM19J1W3XiRNKVuqivj4amjieOyzz3evtXzlVXV0wnTpyM6tSzMzLVJ3XsvP51Orh+\nM9V36ze+MOSxfIWljsdUl6dzXsZUlw+sswyNf0lKJJw5R5fSMCo6EyoSiSgSiaTvt7W16dVXX9Un\nPvEJzZ8/X/PmzdOCBQt0xx13kJiWQaYhRpYkOYO5V4SJqfOwtjN534dIMlOGmQ/C0LobLLVcxotg\nHZYWdMqBv8f0wIZtkqTlCy70LGiHrZ6locdx/7nBSWfWpTP6nRLGOs5Wp4PnYaW+4FLKMR+mOs51\nHh48Mrr6tosc+1xfetyAyVjr6g+TE3pslatOScR0l6nJlwRuAI4hkDiPOvWPqR0Ax5PTAIQLeRn+\nybUJDoKLwA3AMQQS5+Wq03FnRF2dhw0DGxueDJUDcAwZ/c4brk5NnYeFewjcPmL5F4KIQOK8XHVq\n6jws3MNQOQKLRpE/uJCL86hT9EfgBgAgCxM7AAyVAwBgEQI3AAAWIXADAGAR5rgBwHImzsPCPfS4\nAQCwCFcHC4kwXe3HT9Sz+6hj91HH7uPqYAAAhASBGwAAixC4AQCwCIEbAACLELgBALAIgRsAAIsQ\nuAEAsAjnwCFXAAAGY0lEQVSBGwAAixC4AQCwCIEbAACLELgBALAIgRsAAIsQuAEAsAiBGwAAixC4\nAQCwCIEbAACLELgBALAIgRsAAIsQuAEAsAiBGwAAixC4AQCwCIEbAACLELgBALAIgRsAAIsQuAEA\nsAiBGwAAixC4AQCwCIEbAACLELgBALBIRbEvfPnll/XLX/5Sq1evHvLcypUrtWPHDtXU1CgSiWjd\nunWqra0tqaAAAKDIwL1y5Up1dnZq0qRJGZ/fu3evnnrqKZ1yyiklFQ4AAAxU1FB5Q0OD7r33XiWT\nySHPJRIJHThwQMuXL9fcuXP1wgsvlFxIAABwUs4e96ZNm/TMM88MeOyhhx7SNddcozfffDPja44c\nOaLm5mYtXLhQx48f1/z58zV58mRNmDDBuVIDABBSOQP3nDlzNGfOnILecOTIkWpublZ1dbWqq6s1\nbdo0vfPOO8MG7vr6aEGfg8JRx96gnt1HHbuPOjaX41nl+/fv17x585RIJBSPx7V9+3ZNnjzZ6Y8B\nACCUis4qj0QiikQi6fttbW0aO3asrrjiCl1//fW68cYbVVFRoVmzZmn8+PGOFBYAgLCLJDNlmAEA\nACOxAQsAABYhcAMAYBECNwAAFiFwAwBgkaKzyp2QSCR077336o9//KMqKyv14IMPauzYsX4WKRDi\n8bj+4z/+Q3/961/V29urW265RePHj9eyZctUVlamc845R/fcc8+AVQEozsGDBzVr1iy1tbWprKyM\nOnbBj3/8Y7322muKx+P66le/qoaGBurZQYlEQnfddZf+9Kc/qaysTA888IDKy8upY4fs2rVLjzzy\niNrb23XgwIGM9bpx40Y9//zzqqio0C233KLLLrss53v62uN+5ZVXFI/H1dHRoSVLlqi1tdXP4gTG\nL37xC5166ql69tln9eSTT+r+++9Xa2ur7rjjDj377LNKJpP61a9+5XcxrRePx7VixQqNHDlSyWRS\nDz30EHXssDfffFNvvfWWOjo61N7erj//+c8cyw77zW9+oyNHjui5557TbbfdpkcffZQ6dsgTTzyh\nu+++W/F4XJIyniO6urrU3t6ujo4O/eQnP9Hq1avV29ub8319Ddw7duzQjBkzJEnnnXeedu/e7Wdx\nAuPqq6/W4sWLJZ1sTVdUVGjv3r268MILJUmXXHKJtm7d6mcRA+Hhhx/W3LlzVV9fL0nUsQs6Ozs1\nYcIE3XrrrfrmN7+pK664Qnv27KGeHTRixAjFYjElk0nFYjFVVlZSxw4ZN26c1qxZk76uR6ZzxNtv\nv62GhgZVVlaqtrZW48aN07vvvpvzfX0N3N3d3QMu91leXq5EIuFjiYJh1KhRqqmpUXd3t7797W/r\nO9/5zoB6HTVqlGKxmI8ltN+LL76oU089VRdffLEkKZlMDrjoDnXsjA8//FC7d+/WY489pvvuu08t\nLS3Us8MaGhrU29urq6++WitWrFBzczN17JCrrrpK5eXl6fv967WmpkaxWEzd3d2KRqMDHu/u7s75\nvr7OcdfW1qqnpyd9P5FIqKyMfDkn/O1vf9Ptt9+um266STNnztSqVavSz/X09Gj06NE+ls5+L774\noiKRiLZu3ap33nlHy5Yt06FDh9LPU8fOqKur0/jx41VRUaGzzjpL1dXV+uc//5l+nnou3ZNPPqmG\nhgZ997vf1d///nfNnz9fx48fTz9PHTunf3zr7u7W6NGjh8TBfOrb1yjZ0NCgLVu2SJJ27tzJFcQc\n8sEHH+jmm2/W0qVLNWvWLEnSxIkT9dvf/laStGXLFl1wwQV+FtF6P/3pT9Xe3q729nZ95jOf0fe/\n/31dfPHF1LHDzj//fP3617+WJP3jH//Q0aNHNW3aNOrZQUeOHFFNTY0kafTo0Tp+/LgmTZpEHbsg\n03l4ypQp+t3vfqfe3l7FYjHt27dP55xzTs738bXHfeWVV6qzs1NNTU2STk7co3Tr169XLBbT2rVr\ntXbtWknSXXfdpQcffFDxeFzjx4/X1Vdf7XMpgyUSiWjZsmVavnw5deygyy67TNu2bVNjY6MSiYTu\nueceffKTn6SeHbRo0SLdeeedmjdvno4fP66WlhZ99rOfpY4dlMrIz3SOiEQimj9/fvriXHfccYeq\nqqpyvx97lQMAYA8mlAEAsAiBGwAAixC4AQCwCIEbAACLELgBALAIgRsAAIsQuAEAsMj/B+1h+ow+\nv0tNAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def create_data(N, period=2.5, err=0.1, rseed=0):\n", " rng = np.random.RandomState(rseed)\n", " t = np.arange(N, dtype=float) + 0.3 * rng.randn(N)\n", " y = np.sin(2 * np.pi * t / period) + err * rng.randn(N)\n", " return t, y, err\n", "\n", "t, y, dy = create_data(100, period=20)\n", "plt.errorbar(t, y, dy, fmt='o');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this, our algorithm should be able to identify any periodicity that is present." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Choosing the Frequency Grid\n", "The Lomb-Scargle Periodogram works by evaluating a power for a set of candidate frequencies $f$. So the first question is, how many candidate frequencies should we choose?\n", "It turns out that this question is *very* important. If you choose the frequency spacing poorly, it may lead you to miss strong periodic signal in the data!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Frequency spacing\n", "First, let's think about the frequency spacing we need in our grid. If you're asking about a candidate frequency $f$, then data with range $T$ contains $T \\cdot f$ complete cycles. If our error in frequency is $\\delta f$, then $T\\cdot\\delta f$ is the error in number of cycles between the endpoints of the data.\n", "If this error is a significant fraction of a cycle, this will cause problems. This givs us the criterion\n", "$$\n", "T\\cdot\\delta f \\ll 1\n", "$$\n", "Commonly, we'll choose some oversampling factor around 5 and use $\\delta f = (5T)^{-1}$ as our frequency grid spacing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Frequency limits\n", "Next, we need to choose the limits of the frequency grid. On the low end, $f=0$ is suitable, but causes some problems – we'll go one step away and use $\\delta f$ as our minimum frequency.\n", "But on the high end, we need to make a choice: what's the highest frequency we'd trust our data to be sensitive to?\n", "At this point, many people are tempted to mis-apply the Nyquist-Shannon sampling theorem, and choose some version of the Nyquist limit for the data.\n", "But this is entirely wrong! The Nyquist frequency applies for regularly-sampled data, but irregularly-sampled data can be sensitive to much, much higher frequencies, and the upper limit should be determined based on what kind of signals you are looking for.\n", "\n", "Still, a common (if dubious) rule-of-thumb is that the high frequency is some multiple of what Press & Rybicki call the \"average\" Nyquist frequency,\n", "$$\n", "\\hat{f}_{Ny} = \\frac{N}{2T}\n", "$$\n", "With this in mind, we'll use the following function to determine a suitable frequency grid:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def freq_grid(t, oversampling=5, nyquist_factor=3):\n", " T = t.max() - t.min()\n", " N = len(t)\n", " \n", " df = 1. / (oversampling * T)\n", " fmax = 0.5 * nyquist_factor * N / T\n", " N = int(fmax // df)\n", " return df + df * np.arange(N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's use the ``gatspy`` tools to plot the periodogram:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "750\n" ] } ], "source": [ "t, y, dy = create_data(100, period=2.5)\n", "freq = freq_grid(t)\n", "print(len(freq))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFVCAYAAADc5IdQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VeW9L/Dv2vOUOSGMYZTBAQQUbSsHS4VKa3u0Ig1q\nqq2393nu1HO91R7bc+RCbz3E+ng8PRWPVXtsiwqtFa1Si4qgIIIMgWCYwphACJmnPWRPa90/1l5r\nT9k7IWRnryy/n+fxgey9SV6WJN/9e9/f+y5BkiQJRERElDWGbA+AiIjoi45hTERElGUMYyIioixj\nGBMREWUZw5iIiCjLGMZERERZNqAwrq6uRkVFRdLj27Ztw/Lly1FeXo7XX399yAdHRET0RWDq7wUv\nvvgi3n77bTidzrjHg8EgKisr8cYbb8Bms2HlypVYvHgxioqKMjZYIiIiPeq3Mp44cSKeffZZJJ4N\ncvr0aZSVlSEnJwdmsxnz58/Hvn37MjZQIiIiveo3jJcuXQqj0Zj0uNvtRk5Ojvqx0+lET0/P0I6O\niIjoC2DQDVw5OTnweDzqxx6PB3l5eWn/DE/eJCIiStbvmnEqU6ZMQV1dHbq6umC327Fv3z489NBD\naf+MIAhoaWH1nGklJTm8zhnGa5x5vMaZx2s8PEpKcvp9zYDDWBAEAMDmzZvh9XqxYsUKPPbYY3jo\noYcgiiKWL1+OUaNGDX60REREX1DCcN+1ie/CMo/vdjOP1zjzeI0zj9d4eAykMuahH0RERFnGMCYi\nIsoyhjEREVGWMYyJdC4YCqPy1SrsPdaU7aEQUQqD3tpERCPDkbMdqD3fidrznVgwqzTbwyGiPrAy\nJtK5sChmewhE1A+GMZHOhUWefEekdQxjIp1jGBNpH8OYSOfCYYYxkdYxjIl0jmvGRNrHMCbSOZHT\n1ESaxzAm0rkQw5hI8xjGRDoXCnOamkjrGMZEOhcMMoyJtI5hTKRzgRDDmEjrGMZEOhcIhrM9BCLq\nB8OYSOeUythoELI8EiJKhWFMpHPBkFwZGxjGRJrFMCbSuUCQlTGR1jGMiXSO09RE2scwJtI5TlMT\naR/DmEjnlMpY4kFcRJrFMCbSOWVrE28YQaRdDGMinVMqY97XmEi7GMZEOhdUwpj3NSbSLIYxkc5F\np6klSFw4JtIkhjGRzsWeTc0sJtImhjGRjomSpE5TA2ziItIqhjGRjgUT7tgU4roxkSYxjIl0LDGM\nRc5TE2kSw5hIxxJvn8iOaiJtYhgT6VggoTLmXmMibWIYE+lYUmXMBi4iTWIYE+lY4poxK2MibWIY\nE+kY14yJRgaGMZGOhRIqYZGVMZEmMYyJdCwxfDlNTaRNDGMiHVP2FQuC/DHDmEibGMZEOqY0T5tN\n8rc6u6mJtIlhTKRjyl2azMZIGLOBi0iTGMZEOqZMS6uVMY/DJNIkhjGRjilrxmoYszIm0iSGMZGO\niWplbATANWMirWIYE+mYUhmbjHI7NfcZE2kTw5hIx5Ql4mg3NcOYSIsYxkQ6pjZwsZuaSNMYxkQ6\nlrxmzDAm0iKGMZGOJXVTs4GLSJMYxkQ6JonxDVysjIm0iWFMpGPhpMqYYUykRQxjIh1T1owtypox\nG7iINIlhTKRjSiFsinRTizwOk0iTGMZEOiYlnk0dZgMXkRYxjIl0LLmbmpUxkRalDWNRFLFq1SqU\nl5ejoqIC9fX1cc9/8MEHuPvuu7F8+XJs2LAhowMlosuXdNcmhjGRJpnSPbl161YEg0Fs3LgR1dXV\nqKysxHPPPac+v3btWrz11luw2+345je/iTvuuAM5OTkZHzQRDYyYeD9jhjGRJqUN46qqKixcuBAA\nMGfOHNTU1MQ9bzab0d3dDYPBAEmSIAhC5kZKRJdNiiwR8xaKRNqWNozdbjdcLpf6sdFohCiKMBjk\nb+zvf//7uPvuu2G327F06dK416ZSUsLKeTjwOmfeSLjGVpsZAFCQ75A/tppGxLgVI2msIxWvsTak\nDWOXywWPx6N+HBvEFy9exKuvvopt27bBbrfj0UcfxZYtW3D77ben/YItLT1DMGxKp6Qkh9c5w0bK\nNfZ4/ACAXl8AAOD2+EfEuIGRc41HMl7j4TGQNzxpG7jmzZuHHTt2AAAOHTqEGTNmqM/5/X4YDAZY\nLBYYDAYUFhaip4f/U4m0JPEErhDXjIk0KW1lvGTJEuzatQvl5eUA5IatzZs3w+v1YsWKFbjrrrtQ\nXl4Oq9WKiRMn4q677hqWQRPRwCTuMxYZxkSalDaMBUHAmjVr4h6bPHmy+vsHH3wQDz74YEYGRkRX\njt3URCMDD/0g0jHljokmnsBFpGkMYyIdY2VMNDIwjIl0TOQJXEQjAsOYSMcSz6ZmAxeRNjGMiXSM\nlTHRyMAwJtKxpPsZM4yJNIlhTKRjSvgaBEBAdNqaiLSFYUykY0r4GgwCDAYBzGIibWIYE+lYtDIW\nIAisjIm0imFMpGNxlbEgQGIYE2kSw5hIx+IrY0E9kYuItIVhTKRjSvO0IAAGA6epibSKYUykY6Io\nqVWxAE5TE2kVw5hIx0RJgiHyXW4wCOA2YyJtYhgT6ZhSGQPyXmNWxkTaxDAm0jFRkiAY5DCWG7gY\nxkRaxDAm0jFRRLQy5qEfRJrFMCbSMVGSYFQrY3ZTE2kVw5hIx+Q1Y/n3PPSDSLsYxkQ6Fr9mDHZT\nE2kUw5hIx+K7qQVOUxNpFMOYSMckKRrGgiBAYmlMpEkMYyIdC4vRBi4e+kGkXQxjIh0TJcStGbOB\ni0ibGMZEOpbYTc3KmEibGMZEOiZJEgzcZ0ykeQxjIh0Tpfhuak5TE2kTw5hIx8JibGUsQBSzPCAi\n6hPDmEjH4s+mZgMXkVYxjIl0TIq5n7HAQz+INIthTKRjiSdwSRKrYyItYhgT6ZQoSZAQM00d2eLE\nKCbSHoYxkU6JkU3FsQ1csY8TkXYwjIl0SpmONsQchyk/nrUhEVEKDGMinVK2MUVvFBF5nGlMpDkM\nYyKdUkI39jhMgA1cRFrEMCbSKTFxmlpdM87akIgoBYYxkU6Fkxq45Mcl9lMTaQ7DmEinJCWMhcTK\nmGFMpDUMYyKdUjJXrYzZTU2kWQxjIp1S9xmrDVyRx5nGRJrDMCbSqWg3dfw0NbOYSHsYxkQ6lXwC\nV/zjRKQdDGMinUq1tYn7jIm0h2FMpFNiQje1ejY1w5hIcxjGRDqldlMra8aG+MeJSDsYxkQ6pVTG\nQuS7XOA0NZFmMYyJdEqZjjYmHYfJMCbSGoYxkU4lrhlzaxORdjGMiXQqsZuat1Ak0i6GMZFOqWvG\nagMXK2MirWIYE+lUtJta/pWVMZF2MYyJdEqpjJMauBjGRJpjSvekKIpYvXo1amtrYTab8cQTT6Cs\nrEx9/vDhw3jyySchSRJKS0vx5JNPwmKxZHzQRNS/5DXjyDS1mLUhEVEKaSvjrVu3IhgMYuPGjXjk\nkUdQWVmpPidJElatWoXKykq89tpr+NKXvoQLFy5kfMBENDBJ3dTqoR+sjIm0Jm1lXFVVhYULFwIA\n5syZg5qaGvW5s2fPIj8/Hy+//DJOnjyJRYsWYcqUKZkdLRENmBK6SkUsgId+EGlV2srY7XbD5XKp\nHxuNRoiiPMfV0dGBgwcP4v7778fLL7+M3bt3Y8+ePZkdLRENWORbNbpmbFDWjLM1IiJKJW1l7HK5\n4PF41I9FUYQhMteVn5+PsrIytRpeuHAhampqcPPNN6f9giUlOVc6ZhoAXufM0/o1dl3sBgDk5tpQ\nUpIDl8sa9/FIMFLGOZLxGmtD2jCeN28etm/fjmXLluHQoUOYMWOG+tyECRPg9XpRX1+PsrIyHDhw\nAMuXL+/3C7a09Fz5qCmtkpIcXucMGwnXuLPTBwDwevxoaemBzxsAAHR0ejU/dmBkXOORjtd4eAzk\nDU/aMF6yZAl27dqF8vJyAMDatWuxefNmeL1erFixAk888QR+/OMfQ5IkzJs3D4sWLRqakRPRFVPX\njBNP4GI3NZHmpA1jQRCwZs2auMcmT56s/v7mm2/G66+/npmREdEVSe6mZgMXkVbx0A8inUq8a5PS\nTc2tTUTawzAm0qlUlTGnqYm0h2FMpFNKASxEvsuVNWNOUxNpD8OYSKfU4zAFnk1NpHUMYyKdSt3A\nlbUhEVEKDGMinVJO2lKPw+QtFIk0i2FMpFNqZRz5LlcqZK4ZE2kPw5hIp6SENWMe+kGkXQxjIp1K\nvJ8xK2Mi7WIYE+lUUgMXu6mJNIthTKRTSgNXpDBW9xszi4m0h2FMpFPRBi5WxkRaxzAm0ikJkbs2\nCYlrxlkbEhGlwDAm0imlazq5m5ppTKQ1DGMinYrez1j+mN3URNrFMCbSqcRuakFdM87akIgoBYYx\nkU4l3ShC7aZmGhNpDcOYSKckZc3YkFgZM4yJtIZhTKRT0coYkV8jYcx5aiLNYRgT6VTicZhKNzUL\nYyLtYRgT6RSPwyQaORjGRDoV3dqkNHCxm5pIqxjGRDqlFMDKN3l0mpppTKQ1DGMinUp1C0VOUxNp\nD8OYSKeUNWMh4dAPZcsTEWkHw5hIp9RbKBriD/1gZUykPQxjIp2SxPh9xgI4TU2kVQxjIp1KWjM2\n8BaKRFrFMCbSqeQbRUQeZxoTaQ7DmEin1DXjhEM/mMVE2sMwJtKp6DS1/LFaGfPUDyLNYRgT6ZRy\nuIcgJK4ZM4yJtIZhTKRTyWvG7KYm0iqGMZFOKbPRgnoLRflXZjGR9jCMiXRKlCQIQsw0NStjIs1i\nGBPplCRKagADMdPUbOAi0hyGMZFOiZKkNm0BnKYm0jKGMZFOiSLiK2MDp6mJtIphTKRTcmUc/ZiH\nfhBpF8OYSKdEKX7N2MBDP4g0i2FMpFOiKKlNW0DM/YxZGhNpDsOYSKckKVoNA7Fbm7I0ICJKiWFM\npFOiJKlNW0D0jGpWxkTawzAm0ikx1T5jhjGR5jCMiXRKSmrgYjc1kVYxjIl0SpQQt7WJt1Ak0i6G\nMZFO9TVNLYBrxkRaxDAm0qnE4zABOZBZGBNpD8OYSKcSK2NAnrZmZUykPQxjIp0SJcQd+gEolTHD\nmEhrGMZEOpV4NjUgd1RzmppIexjGRDolSVIflbF8n2Mi0haGMZFOJd5CEVAqY4YxkdakDWNRFLFq\n1SqUl5ejoqIC9fX1fb7u8ccfx9NPP52RARLR4Eh9TFMLAg/9INKitGG8detWBINBbNy4EY888ggq\nKyuTXrNx40acPHkyaTqMiLKr725qVsZEWpQ2jKuqqrBw4UIAwJw5c1BTU5P0/OHDh/Hd736X2yWI\nNESSJEhINU2dnTERUWppw9jtdsPlcqkfG41GiKIIAGhubsa6deuwatUqBjGRxijVb/KhH2zgItIi\nU7onXS4XPB6P+rEoijBEFqHee+89dHR04Ic//CFaW1vR29uLqVOn4s4770z7BUtKcoZg2NQfXufM\n0/I1DobCAACr1RQ3TpPJCMEgaHrssUbKOEcyXmNtSBvG8+bNw/bt27Fs2TIcOnQIM2bMUJ+rqKhA\nRUUFAODNN9/EmTNn+g1iAGhp6bnCIVN/SkpyeJ0zTOvX2B+UwzgUCseNUxIlhERR02NXaP0a6wGv\n8fAYyBuetGG8ZMkS7Nq1C+Xl5QCAtWvXYvPmzfB6vVixYkXca9nARaQdyp2Z+joOM1I0E5GGpA1j\nQRCwZs2auMcmT56c9Lq77rpraEdFRFdEaeNIDGP5OEwxCyMionR46AeRDikNXIkTVkaDwAYuIg1i\nGBPpUKpuaoMgIMwwJtIchjGRDkkp1ox5P2MibWIYE+mQEriJlbGRJ3ARaRLDmEiHot3U8Y8bDNHn\niEg7GMZEOqSuGfd1HCbDmEhzGMZEOqR2Uyc2cBkYxkRaxDAm0qGUh34IAiSA58kTaQzDmEiHUjVw\nKR+ziYtIWxjGRDokpTj0Qw1jTlUTaQrDmEiH0k1Ty88P+5CIKA2GMZEOpTqb2shpaiJNYhgT6VD0\nOMz4x5Vs5pGYRNrCMCbSoVTT1KyMibSJYUykQylvFMEGLiJNYhgT6ZAStkLKBi6GMZGWMIyJdEjd\nZ5xqaxOnqYk0hWFMpEPp7mcMsDIm0hqGMZEOpbqfsdJdzSwm0haGMZEOKWGbfAKX/C3PrU1E2sIw\nJtKh1NPU8q8Sw5hIUxjGRDqUepqaDVxEWsQwJtIhtTJOsbWJ09RE2sIwJtIh3kKRaGRhGBPpUPQ4\nzPjHjTyBi0iTGMZEOqRUvkJCGgvcZ0ykSQxjIh3q/0YRwz4kIkqDYUykQ6kbuCLPM42JNIVhTKRD\nktrAFf84G7iItIlhTKRD6ppxqn3GrIyJNIVhTKRDKQ/9YAMXkSYxjIl0iPuMiUYWhjGRDqXaZ8wT\nuIi0iWFMpEMpu6lZGRNpEsOYSIdSHfqh7DOWxGEfEhGlwTAm0qFUh34oH3KamkhbGMZEOiSm2mcs\ncJqaSIsYxkQ6lGprE28UQaRNDGMiHWIDF9HIwjAm0iFlTdho5KEfRCMBw5hIh5QwTnnoB8OYSFMY\nxkQ6pIStKaGDK9rANexDIqI0GMZEOtRfZRwWudGYSEsYxkQ6lDqM5V9ZGRNpC8OYSIfESOVrSjyB\nS1BO4GIaE2kJw5hIh/qfpmYYE2kJw5hIh9StTYbE4zC5z5hIixjGRDokpghjIw/9INIkhjGRDnGf\nMdHIwjAm0qFwuO/KOHoC17APiYjSYBgT6ZAyDW1MPPSDlTGRJjGMiXQoHJZL36Rp6siHXDMm0haG\nMZEOpeqm5l2biLSJYUykQ6IoQQD3GRONFKZ0T4qiiNWrV6O2thZmsxlPPPEEysrK1Oc3b96MP/zh\nDzAajZg+fTpWr16t7mMkouwJi1JSEAPRBi6ewEWkLWkr461btyIYDGLjxo145JFHUFlZqT7X29uL\nX/3qV1i/fj02bNgAt9uN7du3Z3zARNS/sCglTVED0WnrMKepiTQlbRhXVVVh4cKFAIA5c+agpqZG\nfc5qteKPf/wjrFYrACAUCsFms2VwqEQ0UKIowWhMDmP1BC5WxkSaknaa2u12w+VyqR8bjUaIogiD\nwQBBEFBYWAgAWL9+PXw+H7785S/3+wVLSnKucMg0ELzOmafpa2wQYDIaksZotlnkXy0mbY8/YiSM\ncaTjNdaGtGHscrng8XjUj5Ugjv34qaeeQl1dHX79618P6Au2tPQMcqg0UCUlObzOGab1axwIhCEg\n+fvN7QsCAHy+oKbHD2j/GusBr/HwGMgbnrTT1PPmzcOOHTsAAIcOHcKMGTPinl+1ahUCgQDWrVun\nTlcTUfbJ09TJ397qPmNOUxNpStrKeMmSJdi1axfKy8sBAGvXrsXmzZvh9Xpx7bXX4o033sANN9yA\n733vewCABx54ALfddlvmR01EaYVFUe2cjsV9xkTalDaMBUHAmjVr4h6bPHmy+vtjx45lZlREdEXC\nogSLqa/KmA1cRFrEQz+IdChVNzUrYyJtYhgT6VDKQz94owgiTWIYE+lQqkM/OE1NpE0MYyIdElOE\nMSCfwsUsJtIWhjGRDqWapgbkxkzeKIJIWxjGRDojSVJkmrrvb2+DgQ1cRFrDMCbSGSVn005TszIm\n0hSGMZHOhEURQPK9jBUGQWBlTKQxDGMinVHWg1NVxgZWxkSawzAm0pl+w1hgGBNpDcOYSGeUME45\nTW3gNDWR1jCMiXRGZGVMNOIwjIl0Jhzub80YPPSDSGMYxkQ6E5b6m6Y2sDIm0hiGMZHORKepUxz6\nIYAncBFpDMM4iyRJQkOrh800NKTCYXmfcbqtTRL/zRFpCsM4iz462IDHX/oM7+6uy/ZQSEf666Y2\n8mxqIs1hGGdR9ek2AMCB2pYsj4T0RJlpSVUZC9zaRKQ5DGMt4M9FGkL9dVPLZ1MP54iIqD8M4yxS\nflRKTGMaQv0e+sF9xkSawzDOIkGI/LDkz0UaQv0f+sFbKBJpDcNYA/hjkYbSQG4UATCQibSEYZxF\namHMn4k0hNQwNqbYZ6yEMaeqiTSDYUykM0rIGoTUa8axryOi7GMYawJ/KNLQCYv9H/ohv47/7oi0\ngmGcYS+8cwSvvl8LAPAHwjjd0KU+pzRw8UciDaXoNHXqrU0A14yJtIRhnGF7jjThw6oLAIBnNx3G\nE+sPoPZ8J4Do1iamMQ2l/rY2mU3yt30gyM3GRFrBMB5GR851AABaOn3yA9zZRBmgbm1KsWZsMRsB\nAIFQeNjGRETpMYyzwBTpclUP/eB0IQ2h/qapraZIGLMyJtIMhvEwee2DWvX3asGSonIhuhL9TVNb\nzJFpalbGRJrBMM6g2Ip364ELWRwJfZH0dz9jrhkTaQ/DOIP62zoSnabO/Fjoi6O/E7isyppxUK6M\nG9s8WPPyPrz4ztHhGSARJWEYZ5By95xESvjyaGrKBGWfceppaqWBS4TPH8K/rD+AuqYe7D5yadjG\nSETxGMYZFO7nPnXCEJyH2RsIsQGM4vR3owiLOk0dRkunD57ekPqcP8B1ZKJsYBhnUGig09SD/Pye\n3iD++7/uwLo3awb5GeL5g2G8trUWzR3eIfl8g3GsrgO/3XwUoTDXMwerv2nq2MrY7QvGPdfa3ZvZ\nwRFRnxjGGZRqmnqoNLbJoVlV2zIkn+8vn5zF1v0X8NJfjw3J5xuMpzYcxK6aS6hvcmdtDCPRqYYu\nfH6mDUD0391AKmMljPOcFgBAWxfDmCgbGMYZ1N809ZWWxv7g0E4pnm/qAQCEQtmpSnsD0elSrz+Y\n5pUUS5Ik/Mv6A3jmT9WQJEk95rLfNeNgGD1e+TpPHJ0DAGhjZUyUFQzjDErVTS1F0ldI+PhyDXVo\nKqNIdVhEprm90QBOnD6l1NQT3SBft3A/W5ui+4yj09STlDBmZUyUFaZsD0DP+uumVuJ4sP1XQ30L\nvGyfld0TE8AeXyjNKynW8fpO9fft3f4BTFNHT+Dq9cuzK5NG5wJgZUyULayMMyhVZayEaEPLla2L\nBq+gyan2fCf+7fVqeHu1U4HGVsOsjAfuWF2H+vu27t7o/YwHcAJXjy8AAJgwygWDILAyJsoShnEG\npVoz3rTjDEJhEfXN6cP43KVu+PypK8TgFUxT/8dbNTh8ug3v7T0ffTDLt3TkNPXgnLkYvS1nW3cv\nQv3czzhaGUfXjHOdFhTkWFgZE2UJwziDUk1Td/T4cbaxW/24r2nqI+fa8fPf7cefPzqd8vOnqoyr\njjejxxtIOzblh3BfnyNb25Zjp6l7Y96E7D5yCb96vTquwWuwfP4Q2nUUOJIkodMd/X/d3t2rvoFz\n2PpehVIr46C8ZmyzGGE2GVCUa0Nnj5/byoiygGGcQemOw4xvvkp+3aXItqXtBxtSfo5gH2cL7zve\njP/74m68GnNjiliNbR6sf++E2nEbWz2F1R/C2Uljty8aKr0xneKvvl+L6tNteHdP/RV9fm9vCA8/\n+wn+73/u7b/TfYTw+cMIhsSYbmg/PJE3Nc6UYRy/z9hlNwMAivJskAC09/gzP3AiisMwzqB0a8Kx\nFelgo6+vqrY+sj1p37HmPv/Mpo/PxAW8PxCGPxDGb94+ghORRqBQhvdHpxI7TR17ElQ48sahoyd1\nRfvhgQt4/Lef4ei59pSvOXiyBYGgCE9vCO3dQxc4oiThYqvnipYNBqvLI/89JoxywWgQ0N7dC09v\nCBaTAebIdHQi5UYR/sg0dY4jGsYAO6qJsoFhnCH+QBjr3++7OgWA3piwGey0cF8//KO3Z+z7z7Qm\n/KD1B8M4Xt+Bz442qW8Khnr/8kDFTVMHY6+PPLIeb9/ryJIkYcPWk2ho8WBbVeqZhKaO6Bag2O1A\nV6rylSr880uf4S+fnB2yzwnIjVlP/GE/Xnn/RMrXdEWmqPNdVhTkWNHW3QtPbzDlFDUAGAQBZpMB\nPd4AQmERLrt84EdRLsOYKFsYxhnSXwPSucaePh//8MAFPPzsJ302bnW5/XHnUPddiaW/+0RiNe0P\nhuFJ6KgOBMPYf7wZx9JUmZng9gYhQL6rkFIZh0VRnY7vSWjwevWDWnx8qAGNbV512v3I2faUa56x\nATxUYdztDeBUQ5f6tYfSS5uP4vTFbmyraki5ja3Lo4SxBYW5NnS5A+j2BuGMTD2nYjEZ1NmB2Glq\ngNubiLKBYZwh/W072rI3uv4ZG7CvflCLLncAJ853xr3+/X3n8fCzu3D0XHQbS19hrPzQTlVsJzZB\n+QNhdHviw9jTG8Jzb9XgqY2H0t6EoscbwGdHm9QgBIBDp1px5mJ3yj8DyH/fvcea0OmOnyru8ckV\nnd0aDeNuT1D9u8Q2pb2/rx4fHriA3285gQMxx4H6g2G8t7cez/ypOukNTXNMZZw4QzBYdZeib6rO\nN7uHpMkMkJv8OmLWbptSnBeuhHGe04J8l1zh+gNhOK3pjxCwmI3qDIg6Tc3KmChrGMYZEriMqd6+\n4i62EpIkCX/dfQ4AcOBEdC04bt05Eoj9TTEnVuz+YBjdCZ3XsSGvVKMdfXTZ/vx3+/Cbt4+gNrLW\n3OUJ4N//fBi/+MP+tCFefboNz//lCNZt+jx+bN4AchwWWC0mdZo6drzKNHZYFPHJ4Ub18b/tqQMA\nXD+tGADwxsdn8PmZNnx2tCnu87d0+mAyyv/khypwlK740YUOiJKEs/28ERmoc5HPqzRhXWjx9Pk6\nZc04z2lFvsuqPt5vZWyOrierlXFuZitjfzCMnYcvMuyJ+sAwzpBX0qwXJ+pyB3AhYc9xbEC+uPmo\nuq6sRFzdpR6cbojuL1U6t9OFcSAYRiCYPE3d40m9DarHF8T+48348bpdeHPnmbjn2iLTnI1tclDU\nRG5UAADNkWngYCh5PEp4n44JLlGS4PaF4HKYYYuZpnbHXAd/IIxgKIyPDl5EpzuAayYXAoiuv8+f\nURL3dWKnjb29Ibh9QVw1Pg8AkqpyAHjl/RP407ZTl3VLSqUyvnXuOADA+RShebnOXpKvzcLZYwEg\n7t9Ha5dPnSXojqwZ57osyItUxkDqbU0K5WYRANQ/ZzEbkeswZyQsj51rx89e2IOX3z2Oyler4qp+\nImIYZ8ypmKAciFX/uTfuB1RXzN7RPUea1A5YpWpd87t9OB/zAzoYErHz8MW4ijFRX+vY/qCI7hSN\nUQDw2geXdNAYAAAW2ElEQVS12Lz7HADENUfFntylNEbF/p1bO3ux5bN6/M9/2xk3TqDvysvnD0GU\nJOTYzbBa5CnUvceasP9E/B2pdh5uxKsf1MJqMeL+pdNhs8gVnsko4IaZo+JeeylmaldZIx5T5ECO\nw4wOd/wbkMY2uflry976pK8Zy9sbQs2ZNvXNz6V2L5w2E66eWAAAuNg6RGEc6Sn4yuwxAIALkc78\nHm8AP/mP3Xj6j4cAAJ3KNLXDgnxnTGVs668yjn7rjy50qL8vyrOhvac3bunhSvn8IbzwzlF0ewK4\nbkoR2rp78cyfqrPSfX45JElCW1cv913TsODZ1BlwOZVVrNjDKBKDUwlqty+I3/71aNKfDYVFvPzu\n8bSfv88wDoTRnaYyjj1qMXZvdEtndKzNfYRxc6cPf9p+CgDwwb7z+O7XpsFhNUEQhLgwDgTDsJiN\n6rYml92sTtE//5cj6utyHGb0eIN4c4dcnf/oO9ehtMCBcSVOnG7oRp7TAqvZiMXzxuGzo00wmQxo\n7eqFPxDGxm0n1b/jqHw7ClxWNCU0cMVuBTtwohk3RoL9yNl2HDrZiuVfnYpgSMRjz++G1x+CNyRh\n3tRCNHf4MGl0DkoLHTAIAi62DS6MD59uw5s7z6B88TRMG5+H0w1dKC10YGzkzYMSxsqJafVNbjR1\neNHlDsBmMcJqMcZVxqn2GCssMdueSgtiwjjXhrONPehyB1CQY+3rj162d/fUocsTwN/fMhnf/sok\nvPy34/jkcCM+OtSAJTdMuOLPL4oSvJH+AFc/0/MD4Q+E8eePT6OqtgUdPX6MLXbif33nOpTGvGmh\nkUOSJNQ19eBsYw/qm+T/mtp9cNhMyHdZkeeUZ5VKCxwYXeTA6EIHinJtKY+TzRSGcQYM9h3/E+sP\n9PuaHm8Ah0+3JT3u6U1uHAqLYtydezx9hHFvIDTgKkiM3J7PIAjqNDQgB6+3N4iLLR6YTQYEQyKq\nT7Wqz3/yeSM++bwRt90wHvfeNj0ujD89cgk9ngCmT8gHALgc5j7HM6bQgR5vFzy9IYwpcmDWJHmK\neurYPJxu6MbUcfL08723Tcc9t07Db94+gkOnWvFx9UV8fOii+nlKCuzIz7GivtkNnz8Ee6TR6WTk\njYTJaMDx+k5IkoSWTp9agea6LAiFRPWH/o6DF1BW7EBYlDC60AGzyYBRBXZcbPFAkiQIQupv5IZW\nD/62pw4TR+dgyQ0TEBZFrH/vONq6/fjlawfxw29djd5AGDeV5UMQBIwvceFYXQc8vUFsP3hB/Tx7\njzWjvbtXDc28mDVjRz+VsTVmzVhp4ALiO6ovJ4wvtXvh9gYxttgZN0Xe0unDe3vPoyDHittvKoMg\nCLjn1qnYf7wZmz89h4Wzx8BmGdyPIUmS8GnNJfz5o9NqI9uCWaPw/WWzYLX0vce6P25fEP/2ejXO\nXOyGy27GzLJ8HK/vxM9/vw//7c5rce3kokF93ssRFkW0dsknqfkDYYwtdiLHYen/D45wHT1+nG7o\ngs1qhNNmhsNmQnGeLeXdxwai9nwn/vzxaZy6EC0UjAYBowrs6A2Ecbaxu8/DmUxGA0YX2jG6MBrQ\npYUOjMq3w2U3p/3+HiyGcQZkcp9uqiq2tSt5q05jmxcn6jsxZWwuJoxyxe3jVXh6Q30GeV8kSZ4+\nz3dZ4s5Dbu7w4VRDFyQAC2ePwbaqhj7fMHx0sAF3L5oaNwX/hy3yHtopY+W7BuU6LJD6eC8zusiB\n2sg31Lhip/r4nQsnw2E14Ws3jAcg3xzBajGiOBIqO6ovxn2eq8bnq28UOnr8sFtNkCQJ5xq7UZJv\nw9SxedhztAkX27zYdyzaAPbWzjMQICDHYUZxng1Hz7ajNtLxPrpIrpjGFjtxqd2Ljh4/CiPNUIkC\nwTCe3ngQne4APq25hNlTilDX1IO2br+8pSsYxgvvyDMfM8rkNyhKGO+uuQSfP4wFs0ahqrYFW/ef\nh9cfwuxpckDkx1bG9vTf2qbIbTILcqxxP1iUKeu6Sz2YNi4PHx64gA/2n0dhjhXfXXyVetJXrHd2\nncWbO+U91rlOC3563zy1inx9+ymEwiLuuXWq+gYgx2HB1xeU4S+fnMWWz+px58Ipacfa0ePHzuqL\naOrwwmg0oDDHCrPJgD1Hm9DQ4oHFZMD104rR1t2LvceacanNi4dXzIl7czIQLZ0+/OrPh3Gx1YMv\nXzsaDy6bCZPRgE9rGvH7LSfw7KbP8bP756OsNPkaDIXWTh92Hm7EzsMX4444FQRgxoR8LLi6FLdc\nN0ZtQtSLc5e68f6+89h3rDkpGItybVi6YAL+bvbYy3qDdaHFjTc+Oo3qyM+h66cVY+70YkwszcHY\nYqd6DeVelSA6e/y41O6V/2vzojHya1+Nk3arCaMK7BiVb5d/LbCjtMCBknw7chzmQf//SfsdK4oi\nVq9ejdraWpjNZjzxxBMoKytTn9+2bRuee+45mEwm3H333bjnnnsGNQi9SRfGJqPhitag2lKcHLV1\n/4Wkx1b9dq/6+3yXBV+NNBn1Jc9pUasLACgb5erzRhZtXb34YN95dbr06kkFOHquA/uPy+us104p\nwoETLern+sd75+LgyVYcOtWK5g4fDp7sez1W2Q41ttjZ5zvVSWNysaO6UX2NwmYx4du3TE56/agC\nOwB5DddpM+FbX5mMaePy4LKbUZIvP9fY5sXYYidau+RTq66eVIiZEwuw52gTjtd1YN/xZlhMBiye\nNx5b9tZDgoSlN8rTqmcbe7CtSr7mY4rk8Uwek4Oq2hacaujCgoQwliR5KvWjgw3odAdQlGtFW7cf\nb+86h0Ckye1nFfPx3Fs1aGqX17pnTJDXocePkj+/cnLa7KlF8PiCOBLZ5jYtMivgsJrUf1/9rRkr\nSxaFCdWv0hT3+Zk2zJtegte3n0JYlNDc4cMzr1fjnyvmozhy/QDgr7vP4c2dZ1GUa8PVkwqw83Aj\nntp4EA+vuB4dPb3Yf6IFU8fm4qarS+O+ztIbJ+DjQw346+46zJte0mfAdbn9+P2W49hZ3djnbInR\nIGDBrFG459ZpKMqzIRQWI3vPL+LXmz7HT1bOjesaT+fouXY8/5cjcPuCWHrjBKxYPA2GyJuUL187\nBlazCeve/Bz//sZhPP7AjchzXn6lKkoS/IEwegNh9AZC8PnD8AVCuNTmxWfHmtTqzW414uZrSpHr\nsMBoFFBb34njkf/e33se9y65algq9P6EwiLOXOzGpXYvGts8aO/2ozjPhnElTowrdmFciTNlMImS\nhOpTrXh/73l1G+fYYie+dE2puuzQ6Q7gYG0LNmw9ibc/OYvF88bjq/PGxe0aSNTa5cNfdp7FpzWX\nIAGYPiEfy2+dqn6PJDIIAnIdFuQ6LEn/BpVz3y+1edDY7kVzh0/+r9OHhhZP3LbGWDaLES67GS67\nGaZIr88zD9/az9XsJ4y3bt2KYDCIjRs3orq6GpWVlXjuuecAAMFgEJWVlXjjjTdgs9mwcuVKLF68\nGEVF2f9Hkm0vvJ28pqtYOGcMtqc5JWqw+qpEY3W6A/gwxdede1UxzCYD9h5rxqh8O+5fOh2lhQ78\n4/O7k177L69Ep9KX3VQGm8WIo+c68MnnclBOGZuLscVONYynjc/DjLICTB6Ti9+8fQS7Ig1mZaUu\n1Dclh/34ElefnbZlo6LfKLFhnMqEUS7191PH5akhKn9t+XPtOXIJHT298EW6saeNy8PMSDX6wf7z\naO7w4YYZJVh2cxnqmnpQkm/HbfMnoC5y5KgyfuUbfUaZHJ4nzndiwaxo+DS0uPHsm9GQddnNWPXg\njah8tQq7j1wCAJQW2DG+xIn/cscsfFpzCeOLneo08fgS+e/SGDmvfNr4fHh8oaQwFgQB+S4LWrt6\n+w1j5RrnJ4RxcZ4d44qdOFbXgT9uO4lASMSDy2YiEAzjta0n8eRrB/F/vjsHowrs2PJZPd74+AwK\nc634x3vnojhfntZ7/aPT+H+/24dgWIQgAOW3XZU0rWe3mvD9b8zCM3+qxkubj+KR8rnIjQRcKCzi\nwwMXsPnTc+qyxNIbJ+CaSYUIixJau3rh9gVxzeTCuDVik9GA7319BgJBEbuPXMJLfz2GH95xtdr8\nKEoSTp7vxIUWD9p7ehEOy8sux+s6caHFDaNBwAO3z8Ci65PftM6fUYK7/m4K3txxBk9vPIT/cde1\nKdeQRUlCQ4sHJ+o7UHepB+db3Gjp9KHXH065/18AMLMsH1+6djQWzCxNqgLbu3vx1z11+OhgA/71\nj9W4ZnIhbr+pDFdPLFCvrT8YxskLnTh2rgNt3b3w+kPo9Ydht5qQ6zQjx2GJHJNqgMko/1qQ70A4\nEEJRng0l+XY4baZ+p2Bbu3z4+NBF7Ky+mLb502YxYmZZAa6ZXIhxxU5IkoSwJOFSmxdbD1xQe02u\nnVyIpQvk/7+JX7vHG8CHBy5gW1UD3vn0HN7dU4f5M0qweN54XDU+D4IgIBQW0eMN4r299dhWdQGh\nsITxJU4sv3UqrptSNOgpZUEQUJAjn2ynLIspRElCZ49fDefmDh9au3xw+4Jwe4Po8QXR0OpBKCxC\nSHUcYoK0YVxVVYWFCxcCAObMmYOamhr1udOnT6OsrAw5OfIPtvnz52Pfvn24/fbbL+svrEeJndT3\nLZmO0gI7nHaz+sM3nVtmj0nbFT1Y3Z4AXHYzfrJyLv7w/gmMKXRg5+FG3HzNaPW0La8/hGunJL+h\nKi10qGECAN/+yiTcuXAKTtRHG7yK82zIdVhQWmDHsboOmIwGdb3nmsmFEASoAXL9tGI1zL5y3Wjs\nOdKE2VOLUJBjjVvPBOTwHV/ijPu4P+NjwliZAlcoYXygtiXuwJDZU4tQkm/HmCKHGnw3zipFjsOC\nR1fOVV83eUwODAIgSvKMgxIik0bnwGI2YHtVA46caUeO04z500dhy956dHsCmDouF8GQiB98YxZy\nHBbcev04bPjwJABg7vQSCIKAqWPzMHVs/Lv4scVOCIK8TJDntKAkz4bZ04qw4cOTsFqMGBdzbfLU\nME4/Ta2EcV/rwtdNLcKWz+qx91gzxpc4cct1Y2AwCPAFwnhzxxn880ufwWE1wdMbQkGOFT9ZOVet\nlpfdPBFFeTa8/O5xlBY4ULF0etLfR/06U4pw69xx+OhgA/7pxT24YeYoBEMiPj/Thh6vfAOLe2+7\nCrfOHRdXYaVrpBIEAQ8um4nWLh/2H29GY6sHi+ePR0uHD/tPNPd52IvJKGD21CLc8eVJKSsoALjj\nSxPR5fZjW1UDVv9uH/7+K5MxZWwuchxmtHX34mKrF7XnO3GiviNu6ccS6SdwWE2wWU2wW02wW4yw\nWUywWY3IdVowZ2px2jX6wlwbKpbOwKI5Y/HHbadw5Gw7jpxtx5giB2wWI7z+cJ+d38q/m4Fy2c2Y\nNDoHk8bkqGvVOXYzenxBnKiX/26nGrogSXKT4NfmjUfZaBfGFDlRmGNFS6cPF1s9ON/sxrH6Thw6\nJc+K9XXNF84egyU3TlDfbPYlx2HBnQunYNlNE/HpkUvYVnUBe481Y++xZtgsRgRDYtxMWlGuDXf9\n3WTcfPXojDZgGQQBhbk2FObaMDOyk+JKpf2OdbvdcLmiF8poNEIURRgMBrjdbjWIAcDpdKKnp++y\n/YvGbjXFnf5UmGtVA84fCKtTyj9ZORe/3HAw6c9/f9lMXGh241zCNMg3bp6IdyMHXKSz4OrRgCSi\nvsmNS+3xJzfdct0YjB/lws/unw9JkrDkxgkYW+xU16Jj12PX/tebAch7mJ02Ex5+dhcAOYjv+PIk\nAMBVE/IxeUwuzjZ247b58rrt128qw9FzHSi/7Sr1c7nsZiyYVYrPjjbJXc/zx2Pzp3UQJQnf/NIk\nfHfxVWrjz+ypRbjp6lIsu6kMDpsJDqs5broxtvs3FafNjKsnFaCjx49brhsT91ye04KJpTlqhQvI\nAav8kP/6gjL87m/HUZxnw+ypyW9MzCYjVtw2Axs/OBG3ncpkNODrN5bhnU/Pye+WO3043SBPv9+3\nZDq+Frk+6v+L2WNw8kIn8pxWfCtyPftiNRux/Nap+Px0G26cVQpBEFBa4MCCWaNQkGONa3AZU+jE\n+Sa3+gYhlbtvnYoNW0+q+5hjzYmEcUm+DT9aPlv9ofatL09CSb4NHx28iJZOH266uhTfuHli0vr4\nglmluG5KEaxmY78/EO9fMh1jihzY9PEZtdEu12nB0hsn4IFvXQu/9/L3I5tNBvzve+bgzx+fxvaq\nBqx/T+5LsJgNuGX2GFwzqRBFuTaYTAJEUd7uZu/nxDJADvr7l87AVePz8fstx9XdAomKcm24flox\nZpQVYOq4XJQWOIYsGMpKc/Doyrk429iN9/bW48CJFhgMAuxWE8YVOzFrUgGumVSoNtJZTAb0BuSD\nfdzeIIIhEcGwiFDkV7vDgsZmN1o7fWjt6kVDqxs1Z9tRk+JoV0GQmyYXXT8WN84clbQMUJhrU2eI\nAHkd/Mi5drR3+2E0CBAMAuwWI26cVXpZU/1WixFfnTsOt14/FrXnO7GtqgEX2zywmY2wmI1yFT6x\nALdeP06dCRlpBCnNPpzKykrMmTMHy5YtAwAsWrQIH3/8MQDgxIkTePrpp/HCCy8AANauXYv58+dj\n6dKlwzBsIiIi/Uj7FmLevHnYsWMHAODQoUOYMWOG+tyUKVNQV1eHrq4uBAIB7Nu3D9dff31mR0tE\nRKRDaStjSZKwevVqnDghT/OsXbsWR44cgdfrxYoVK7B9+3asW7cOoihi+fLluPfee4dt4ERERHqR\nNoyJiIgo80bmSjcREZGOMIyJiIiyjGFMRESUZQxjIiKiLBuWMBZFEatWrUJ5eTkqKipQX18/HF/2\nC6m6uhoVFRXZHoZuBYNBPProo7jvvvtwzz33YNu2bdkeku6Ew2H89Kc/xcqVK3Hvvffi5MmT2R6S\nbrW1tWHRokU4e/ZstoeiS3fddRcqKipQUVGBn/3sZ2lfOyx3bUp3xjUNnRdffBFvv/02nM7+j4uk\nwXnnnXdQWFiIp556Cl1dXbjzzjuxePHibA9LV7Zv3w6DwYANGzZg7969eOaZZ/jzIgOCwSBWrVoF\nu93e/4vpsvn98ulx69evH9Drh6UyTnfGNQ2diRMn4tlnnwV3q2XO7bffjh/96EcA5Bkfo3Fw982l\n1G677Tb8/Oc/BwA0NDQgLy/1edE0eL/85S+xcuVKlJSUZHsounT8+HH4fD489NBDeOCBB1BdXZ32\n9cMSxqnOuKahtXTpUoZDhjkcDjidTrjdbvzDP/wDHn744WwPSZeMRiMee+wx/OIXv8Add9yR7eHo\nzqZNm1BYWIhbbrkFAPgGPgPsdjseeugh/Pa3v8WaNWvwyCOPpM29YQljl8sFjyd6k2blZhNEI1Fj\nYyMeeOAB3HnnnfjmN7+Z7eHoVmVlJd577z08/vjj6O1NvtsSDd6mTZvw6aefoqKiAsePH8djjz2G\n1tbkuyvR4E2aNAnf/va31d/n5+ejpaXv+7kDwxTG6c64JhpJWltb8YMf/ACPPvoovvOd72R7OLr0\n1ltv4Te/+Q0AwGazQRAEvnkfYq+88grWr1+P9evXY+bMmXjyySdRXFyc7WHpyqZNm1BZWQkAaGpq\ngtvtTrskMCwNXEuWLMGuXbtQXl4OQD7jmjJnsDfTpv49//zz6Onpwbp167Bu3ToAwEsvvQSrNfW9\naOny3H777Xjsscdw//33IxQK4Z/+6Z9gsQz8dntEWrB8+XL89Kc/xX333QdAzr10byp5NjUREVGW\nce6HiIgoyxjGREREWcYwJiIiyjKGMRERUZYxjImIiLKMYUxERJRlDGMiIqIs+/8Te8MjrHPZnQAA\nAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from gatspy.periodic import LombScargle\n", "model = LombScargle().fit(t, y, dy)\n", "period = 1. / freq\n", "power = model.periodogram(period)\n", "plt.plot(period, power)\n", "plt.xlim(0, 5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The algorithm finds a strong signal at a period of 2.5.\n", "\n", "To demonstrate explicitly that the Nyquist rate doesn't apply in irregularly-sampled data, let's use a period below the averaged sampling rate and show that we can find it:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFVCAYAAADVDycqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtgVOWBNvBnZjJJJplcSQKKCQYUvKBY8L5Saitat2wr\nC8VYi9V17W7Xb3fbqv3cr58U+pUSa2mtilqtYkEFpHhrVBAERbkGQgjhFkiAJARyTyYzk7me8/0x\nc86cuWYSEvJmeH5/ZXImkzfvTM5z3uvRybIsg4iIiIShH+4CEBERUTCGMxERkWAYzkRERIJhOBMR\nEQmG4UxERCQYhjMREZFg4grn/fv3Y/78+WHf37x5M+bOnYuSkhKsXbt20AtHRER0IUrq6wmvvvoq\nPvzwQ6Snpwd93+12o7S0FOvWrUNqairuu+8+fPOb38SoUaOGrLBEREQXgj5bzuPGjcMLL7yA0L1K\namtrUVRUhIyMDBiNRkybNg3l5eVDVlAiIqILRZ/hfOedd8JgMIR932q1IiMjQ32cnp6Onp6ewS0d\nERHRBWjAE8IyMjJgs9nUxzabDVlZWTF/hjuFEhER9a3PMedoxo8fj1OnTqG7uxsmkwnl5eV4+OGH\nY/6MTqdDaytb10MtPz+D9TzEWMdDj3U89FjH50d+fkbfTwoRdzjrdDoAQFlZGex2O+bNm4cnn3wS\nDz/8MCRJwty5c1FQUNDvAhAREVEw3fm+KxWv0oYer4aHHut46LGOhx7r+PwYSMuZm5AQEREJhuFM\nREQkGIYzDattB85g0RvlcHu8w10UIiJhMJxpWL320WGcOtuDmobu4S4KEZEwGM4kBt1wF4CISBwM\nZxICs5mIKIDhTEJgOBMRBTCciYiIBMNwJiEoO9ARERHDmQTBbCYiCmA4ExERCYbhTEREJBiGMwmB\nY85ERAEMZyIiIsEwnEkIbDgTEQUwnEkIOm5DQkSkYjiTGJjNREQqhjMREZFgGM4kBDaciYgCGM4k\nBqYzEZGK4UxC4IQwIqIAhjMJgUupiIgCGM5ERESCYTiTEGR5uEtARCQOhjMREZFgGM4kBI45ExEF\nMJxJCOzWJiIKYDgTEREJhuFMQmC3NhFRAMOZiIhIMAxnIiIiwTCciYiIBMNwJiIiEgzDmYTApVRE\nRAEMZyIiIsEwnEkIXEpFRBTAcCYhsFubiCiA4UxERCQYhjMJgd3aREQBDGciIiLBMJxJCBxzJiIK\nYDgTEREJhuFMQuCYMxFRAMOZhMBubSKiAIYzERGRYBjOJAR2axMRBTCciYiIBMNwJiIiEgzDmYiI\nSDAMZyIiIsEwnEkIXEpFRBTAcCYiIhIMw5mEwKVUREQBDGcSAru1iYgCYoazJElYsGABSkpKMH/+\nfNTX1wcd37hxI+bMmYO5c+di1apVQ1pQIiKiC0VSrIObNm2C2+3G6tWrsX//fpSWluLFF19Ujy9Z\nsgTvv/8+TCYTvvOd72DWrFnIyMgY8kJT4mG3NhFRQMxwrqiowPTp0wEAU6ZMQXV1ddBxo9EIi8UC\nvV4PWZah4xmWiIjonMUMZ6vVCrPZrD42GAyQJAl6va83/KGHHsKcOXNgMplw5513Bj03mvx8tqzP\nh5FWz9nZaSOuzCOtvCMR63josY7FFDOczWYzbDab+lgbzE1NTXjrrbewefNmmEwmPPHEE1i/fj2+\n/e1vx/yFra09g1BsiiU/P2PE1XNnpx2tppgfR6GMxDoeaVjHQ491fH4M5AIo5oSwqVOnYuvWrQCA\nyspKTJo0ST3mdDqh1+uRnJwMvV6P3Nxc9PTwTaaB4YgIEVFAzKbKzJkzsW3bNpSUlADwTQArKyuD\n3W7HvHnzMHv2bJSUlCAlJQXjxo3D7Nmzz0uhKfFwKRURUUDMcNbpdFi0aFHQ94qLi9WvH3zwQTz4\n4INDUjAiIqILFTchISGwW5uIKIDhTEJgtzYRUQDDmYiISDAMZxICu7WJiAIYzkRERIJhOJMQOOZM\nRBTAcCYiIhIMw5mIiEgwDGcSggz2axMRKRjOREREgmE4ExERCYbhTEREJBiGM4mBQ85ERCqGMxER\nkWAYzkRERIJhOJMQ2KtNRBTAcCYiIhIMw5mIiEgwDGcSA/u1iYhUDGciIiLBMJyJiIgEw3AmIiIS\nDMOZhMC7UhERBTCciYiIBMNwJiIiEgzDmYQgs1ebiEjFcCYiIhIMw5mIiEgwDGciIiLBMJyJiIgE\nw3AmIiISDMOZiIhIMAxnEoLMtVRERCqGMxERkWAYzkRERIJhOBMREQmG4UxERCQYhjMREZFgGM5E\nRESCYTiTELiSiogogOFMREQkGIYzERGRYBjOJAT2ahMRBTCciYiIBMNwJiIiEgzDmYiISDAMZxID\n11IREakYzkRERIJhOBMREQmG4UxCYKc2EVEAw5mIiEgwDGciIiLBMJxJCOzWJiIKSIp1UJIkLFy4\nEDU1NTAajVi8eDGKiorU41VVVXj66achyzJGjx6Np59+GsnJyUNeaCIiokQWs+W8adMmuN1urF69\nGo8//jhKS0vVY7IsY8GCBSgtLcXbb7+NW265BY2NjUNeYCIiokQXs+VcUVGB6dOnAwCmTJmC6upq\n9diJEyeQnZ2N5cuX49ixY5gxYwbGjx8/tKUlIiK6AMRsOVutVpjNZvWxwWCAJEkAgM7OTuzbtw8/\n/OEPsXz5cuzYsQM7d+4c2tJS4uKgMxGRKmbL2Ww2w2azqY8lSYJe78vz7OxsFBUVqa3l6dOno7q6\nGjfffHPMX5ifn3GuZaY4jLR6zsoyjbgyj7TyjkSs46HHOhZTzHCeOnUqtmzZgrvvvhuVlZWYNGmS\neqywsBB2ux319fUoKirC3r17MXfu3D5/YWtrz7mXmmLKz88YcfXc3d07oso8Eut4pGEdDz3W8fkx\nkAugmOE8c+ZMbNu2DSUlJQCAJUuWoKysDHa7HfPmzcPixYvx2GOPQZZlTJ06FTNmzBhYyemCJ7Nf\nm4hIFTOcdTodFi1aFPS94uJi9eubb74Za9euHZqSERERXaC4CQkREZFgGM4kBvZqExGpGM5ERESC\nYTgTEREJhuFMREQkGIYzCYFDzkREAQxnIiIiwTCciYiIBMNwJiHI7NcmIlIxnImIiATDcCYiIhIM\nw5mIiEgwDGcSBAediYgUDGciIiLBMJyJiIgEw3AmIXApFRFRAMOZiIhIMAxnIiIiwTCcSQjs1SYi\nCmA4ExERCYbhTEREJBiGMxERkWAYziQGDjoTEakYzkRERIJhOBMREQmG4UxCkNmvTUSkYjgTEREJ\nhuFMREQkGIYzERGRYBjOJAYOORMRqRjOREREgmE4ExERCYbhTEJgrzYRUQDDmYiISDAMZyIiIsEw\nnEkIMvu1iYhUDGciIiLBMJyJiIgEw3AmIiISDMOZBMFBZyIiBcOZiIhIMAxnIiIiwTCcSQhcSkVE\nFMBwJiIiEgzDmYiISDAMZyIiIsEwnImIiATDcCYiIhIMw5mIiEgwDGcSAldSEREFMJyJiIgEw3Am\nIiISDMOZhCBzizAiIlXMcJYkCQsWLEBJSQnmz5+P+vr6iM976qmnsHTp0iEpIBER0YUmZjhv2rQJ\nbrcbq1evxuOPP47S0tKw56xevRrHjh2DTqcbskISERFdSGKGc0VFBaZPnw4AmDJlCqqrq8OOV1VV\n4d5772W3JBER0SCJGc5WqxVms1l9bDAYIEkSAKClpQXLli3DggULGMxERESDKCnWQbPZDJvNpj6W\nJAl6vS/PN2zYgM7OTjzyyCNoa2uDw+HAhAkTcM8998T8hfn5GYNQbOrLSKvnzAzTiCvzSCvvSMQ6\nHnqsYzHFDOepU6diy5YtuPvuu1FZWYlJkyapx+bPn4/58+cDAN577z3U1dX1GcwA0Nrac45Fpr7k\n52eMuHq29PSOqDKPxDoeaVjHQ491fH4M5AIoZjjPnDkT27ZtQ0lJCQBgyZIlKCsrg91ux7x584Ke\nywlhdC44MkJEFBAznHU6HRYtWhT0veLi4rDnzZ49e3BLRUREdAHjJiRERESCYTiTEGTe+oKISMVw\nJiIiEgzDmYiISDAMZyIiIsEwnEkMHHImIlIxnImIiATDcCYiIhIMw5mEwF5tIqIAhjMREZFgGM5E\nRESCYTgTEREJhuFMYuCgMxGRiuFMREQkGIYzERGRYBjOJATelYqIKIDhTEREJBiGMxERkWAYziQE\nmb3aREQqhjMREZFgGM5ERESCYTgTEREJhuFMREQkGIYzERGRYBjOREREgmE4kxC4lIqIKIDhTERE\nJBiGMxERkWAYzkRERIJhOJMQeFcqIqIAhjMREZFgGM5ERESCYTiTGNirTUSkYjgTEREJhuFMREQk\nGIYzCYG92kREAQxnIiIiwTCciYiIBMNwJiIiEgzDmYiISDAMZyIiIsEwnImIiATDcCYhyDIXUxER\nKRjOREREgmE4ExERCYbhTEJgpzYRUQDDmYiISDAMZyIiIsEwnImIiATDcCYxcNCZiEjFcCYiIhIM\nw5mIiEgwDGcSAnu1iYgCGM5ERESCSYp1UJIkLFy4EDU1NTAajVi8eDGKiorU42VlZVixYgUMBgMm\nTpyIhQsXQqfTDXmhiYiIElnMlvOmTZvgdruxevVqPP744ygtLVWPORwO/OlPf8LKlSuxatUqWK1W\nbNmyZcgLTERElOhihnNFRQWmT58OAJgyZQqqq6vVYykpKVizZg1SUlIAAB6PB6mpqUNYVEpovCsV\nEZEqZre21WqF2WxWHxsMBkiSBL1eD51Oh9zcXADAypUr0dvbi1tvvbXPX5ifn3GORaZ4jLR6Nmek\njrgyj7TyjkSs46HHOhZTzHA2m82w2WzqYyWYtY+feeYZnDp1Cs8//3xcv7C1tWeARaV45ednjLh6\ntvY4RlSZR2IdjzSs46HHOj4/BnIBFLNbe+rUqdi6dSsAoLKyEpMmTQo6vmDBArhcLixbtkzt3iYa\nCHZqExEFxGw5z5w5E9u2bUNJSQkAYMmSJSgrK4PdbsfkyZOxbt06XH/99XjggQcAAD/60Y9wxx13\nDH2piYiIEljMcNbpdFi0aFHQ94qLi9WvDx8+PDSlIiIiuoBxExISAidrExEFMJyJiIgEw3AmIiIS\nDMOZiIhIMAxnIiIiwTCciYiIBMNwJiIiEgzDmYQgcy0VEZGK4UxERCQYhjMREZFgGM5ERESCYTiT\nEDjiTEQUwHAmIiISDMOZiIhIMAznC4jb48V/P/clPvzqxHAXJRz7tYmIVAznC0hTmx09djfeFzGc\niWhE+mjHSfzipe1wuDzDXZSEwnAmIqIBW/dFHdq6HThS3zXcRUkoDGcSAnu1iUa22tPdw12EhMJw\nJiKiAdFuu3vqbM8wliTxMJwTyPHGbny881TU4zrdeSwMESU8uzMwzmztdQ9jSRJP0nAXgAbPb9/c\nCwC4+7bxYA4T0VDr7HGqX2uDms4dW84JyOOVhrsI/ce7UhGNOF3acHYwnAcTwzkB6dhuJqLzoEfT\nlW13eHjr10HEcE5AHFsmovPBpglnSZbhcHmHsTSJheGcIEb6ZAxeb4vF4fLgq6ozkCS+MxSdct7J\nyUgBwK7twcRwThCtXb3q1zo2nekcrdxwFK9/fBgfxZj9T2Tzh3F+tsn/eGQ3EkTCcE5A0aL5Qgxt\nu8OD0jf34tDJjuEuyohS32IFAJw8YxnmkpDIlDAu8IczW86Dh+GcgKKFcLRoPt7YDXuCXvEeqGtH\nTWM3fr+6criLMiQ8Xgkb9zTAYnMN6uuakn2rLHu5PIZisPUqLedUAFxONZgYzglI3493tbHFit++\nuRdL3qoYugLFYagmeY7IZWX9sGlPI1ZtOoa/lB0a1NdNMvgu5bwccx6wD786gTc/PZrQ4/Y2hxtJ\nBh1yMlLVxzQ4uAnJhSRC01kZqz7dajvPhTk/En1lR2Orr/v5THv87197twOZ6UYYkwxRn6MEiv4C\nHAoZDJIsq3d/m3JZHq4ZP2qYSzQ0bL1upJuMSE/1RQm7tQcPW84JKNEDqT/kBJ8HrixdSU2J7zr7\nTLsNT7y0HS9/cDDm85T+Br2e4TwQ3VZXxK8TjbXXjfRUI9IYzoOO4ZyAooZzhO932xP3xAEg4ddo\nKd32SXGOZZzwT/Dad6wt5vOUzSTYcB6YDotD/XqkL3OMxuOVYHN4kJlmRFqqEQDDeTCxWzsByZAj\nTv7S5pQsy9i6vwkr1h89X8UaFgmezf3m7OcmEeyFGZj2CyCclR6BbHOK2q1tcybm3zoc2HJORFFO\nqNqt9WoauvDXBA/mC4EyJhxv973LE98EOWULWG7HODDa2fOJGs5dNt++2lnmZJhS2K092BjOCSjS\n6dTp8qJDs0l9jz36CUOW5fM+yzmecHF7JNQ1WfoVGEMVLnaHGz0iDQnE+WfG20ut1JtnkGYad1gc\nCT1rOZQ2kG0JGs5KyzkrPQWpyQbodTqG8yBiOCcIbQZFCqTfrNiD5/5WFfM1JP/PvfTBQfz4mc/h\ndPu6QLcdOIMvKk8PXmEHaPknh/GbFXuwv7Y97p8Zqobfvy7ehP9+7quhefF+UMaEB/vPdPsvzjxR\nWtqb9jTgk13x7R52+FQnHn9xOz7cdmLQyieC5k47Vn56VP0/0bJpQipRlxd1W30X+9nmZOh0OqSl\nJiXs3zocGM4Joq+W5+m2vpfa/OvTW7CvphV7jrQACNwO7rWPDgvRBb7zYDMAoKG5J+6fGap1ukqr\nWRKk2zfeHoJ4W8Ier+957ig9KG9vOoa1W2rj6mGpONoKANhQ3hDX7x5udocHbk/fY/PP/a0KWypO\nY8Ou+rBj2tay052Ya+27lJaz2bevdlpqEjchGUQM50QxSBmh3Uv5f17Zifo4g1CS5fMWVLHW54by\nDnH3/FC/fixvfHKkz1nXoeIdrlBazO4ILWftbnLaPd2jcfmDLjlJrNON2yOFXbRae934+Qtf4aX3\nYy81AwLjyl1WZ9gxpVvbmKRX//5E09bte+9zM33hnJ6aNOBu7araNpw8y61itcT6b6EBC56JDfx1\n/RF8WdXU79cJXTrz+seH4/q5J17cjl++srNfvyuotdePXO/P0ttzbTlv3d+Eo/WdUY8rLczB5nB5\n8Pm+01HDVPLPtldEKoUkyWEt6njLq7SYI4Vzl2bdbjwtJWU72f5cu+0/3oaHF2/EvprW+H+on97e\nVIOn/rILxxq71O9VHmuDyyOh8nhbnxebyUbfRWKklrFv5yw9zCYjXP5u79qmbvXrRNDQYkOyUa/e\n9CLdZITbI/V7RYDD5cGza6vw6zf2oCmOHr4LBcM5QWhPwjaHG19UNmH5x0eiPv+L/ZGDO3Rf7vpm\na5+/2+3xorPHiebOvltR0fQr4vqx+PZcJjR5vBLe+OQInn57X9TnROv2PVfvbD6OFRuO4r0v6yIe\nD22hhOaIV5Lwv1/eEVb2eFv6nhjhrB1XjOf+vcpr9GdC2Nb9TWjpsOPjIbwr1heVvv+BgycCN0Xp\ntgVawW3djrCf0TL6ewIcrvALFGuvG2ZTEpKNBjjdEqpPtGPxir14Y330/8mRxOOVcKbdhrF56eqK\ngaz0ZADBdRiPk2cCvXM7D50dvEKOcAznBKE9OcdzEtSekLQG8oFQtinsr0bNlqH9mVXdn40xzqXb\nOdpNH7xS4DW9Q9RybvDfFarudOSuvtCZ4qH112N3o93iQE1DV9CxeC8m1HCO8HztZCeHs+9wHkhr\nsa8x78GkDWHtKob2Prrslfc+0ufE1uuB2WRESpIeLrcXxxu7AQTmTYwExxu78fbGmohLwZo77PBK\nMi7JN6vfy/aPPXf1c0e0E5o7n51ttw+wtImH4ZwgtC2cc1myUuM/iUQSLUC1V7798avXd2teO/6f\n08H39z62bBve2xq5Zak4l27t3iitQm1dKyF28qwF/1K6GaVv7o34/P4uTVNCKSU58vh6rKVwQHBg\naNc2x3sx4fb4l1JFajn3alvOfXdrOwYQzpKktLYjH69rsqCqH7P2Q2nfQ+2aZO1FT18ho9Rx6OfE\nK0mwOz1ITzUiOdkAl1sasuGPofT6x4exaW8jlkcY2mrw7+muDedMf8u5v3dIO3E2cP4405FY4dzZ\n4+zXvvdaDOcEsXrzMfXrwyH3Lh6sca6KKON/A1lLHDqe159X0Ol0ONthR2ePE3/ffjLmc+MJ54Mn\nOoK2W1S4NCfd/cfb8PpHhyHJcsRw/vUbewAELm5sDrcaXP/zyg489ZddQT9Tezr6RZDvOb5yG5P0\nqG3qDmu9hJ4AQ98Cbbe3Nky1FwmyLMMrSeHvhSyrPQ5eSQ7qKQh97dc+Oox3+7hAUpYa9WfCoNV/\nK8Jo4f+bFXvw7Nr9A17Hrg1hi10bzoG6ijTRSyHLMnpdShmD/7+U+lFazpIsqxPnRsqNRFo67Tjr\nD8pDpzrDPgPKjXIuKQiEc6Bbu3/h3NxhR7JRj3FjMtDcYU+Y9fAWmwv/8+cd+OWru/p+cgQM5wSh\nvavUh1tr1a/rm3vw70u/GJTf8erfg29L+Pm+09hzpCXsn8nmcKOqti3oxLll32ksXb1PfW7YP2A/\nu7XjbYn21VJs6erF0jWV+OVfwv+BnJpZtn/6WxW+OnAGjS3WoFaQxxseXrIs4z+f/RJPvrwDLrcX\nHRbfeLxSH29trMHilXux92hL1HIpLdauHicWr9iL36zYE3S81xU65hz8d2pbztpuaG3Z3R4Jz6yq\nxJKVwa19rxS8MM/jCX9/tcq2n4QkyWhsiTw/QZkg5HR74w5T5XdYImz0or04Guj9prWtYm0g92gu\nZDpjhLPLLakf2dALCOVCKt0/5gwALf75GDLCJ+mJ6JimB83p8obNPVGGXS7JT1e/p4RzrIuaULIs\no6WzFwXZabh4VBo8XlmdBT7SHahrh8sjwTDAm8cwnBOQthtz4fLyIXldAFix4ShefL8aoTH5x3f2\n49m1VVi6phJt/hbDyg1HcfBkJ7qsTpxutYYFfX9bztpw7+xxoqo28pKiSMG5sbxBHefq9LeYI80w\ndUWYhbtweTn+9vlx9bHHK+HZd/YHPafXPw5rsbuDglE58X9VdQYAUNsUGGt7d2stfrNij3rR4fGX\nW3lOS8hku9DyhrZKta25aC3n8iMtqGnoQm2TBZ2a3eNCL3xClwLZIiyXWbP5OBa8vhv7joX3rigt\nZ1mOPMFMoQ1a5WuXO3z2r/bkbemjez8a7aQli82lBmav06PuotbVEz1ktBdHoS1nm7/Vn24yquHc\n0eP7nMly5M+VaE76u5q/Ne0SAL6g0V58nW61IsucjIy0ZPV7yphzf+7CZbG74XR7UZBjUmd9t/Yx\nEW+otXX39nvGeSTKsMuif7lxQD/PcE5A/R3zGQjt1X9oS6DOHyiHTnbiuXXBu5LpdDr84Z39KD8S\n3Gpcv6u+Xy0KbYAseG0Xnl1bFbQMw+n2ovJYW1i39tkOO1Z9dgwvf1Dtf53ovzPaRhQ7NJN6PF4J\nB08GL7Xq7AmcXA7UBcZFOy2+k71SpiSD799PkmSUbT+FOk1IRtuZ60y7DZ09zrBdqUL/Du1xpSX3\n8c5TQfX+2keBscSahsByorDXCg2fCLtAbdzj22Dk+XUHwt5HbVki7aYFAHuOtODRP27F3qOtkGU5\naIlW6J3TtC3dgX7WtV2vXinw+3qdHhTkpkGv08Ucc9ZeSDhdXhw+2YGfPvclGlussPrrx5xqVNd2\na8s8XLtodfY41V29otl9uBkvf1CNmoYu6HU63H1TEXQA3v/yBBa8vhtrPz8Ou8ONdoszaLwZ8K13\n1ut0ONMR/xhrS6ev67wgx4S8LF84tw9jONc2deMXL+3A8k/iW0Jqc7gj9uJJsoyDJzqQl5WKi0al\nDagsDGcaEG3oxcrU0BafxytFPKE6XF61q6wvja3WoCVCSkvu//5lFz7wzxx/6f1qPLeuCrsOBc+O\nVcKhtcsR9DiSeFo4ESdMaVqWb3wSWDoTbZMF7YQRpdUR6QYVTrcXv3x1Fxa8titCOIe0djXHlTD4\n2+e1iEa71je0dRs6oauvjSZCN/bQhntouS12F55ZtQ8vve+7WNpc0QiHyxv0meqxuWBzuNUxUO34\ne18T46Kx+Os5LyvV99jfeu51epCemoQsc3LM7tnQ1vIzqythsbuxaW+j2lORbjIixRg+oS9SzwPg\nu8hdv6s+5sY/HRZHzHX30bjcXvz6r+VY+EZ51PevtasXL39wELsPt6ChxYriizKQm5mKSUXZ6nM+\n29OIE/4JoIUh4WxMMuDivHQ0tFjjHjdWzg8F2Sb1vRiubm1ZlvGGf/np7sMtfd6wpMPiwC9e8s0n\nCX3PunqcsDs9KL4oM2x5arwYzjQg0WaHh7aako2GoMlPHm/4BKRIr9lldWL34cjLTrZURN/n+4Ov\nTuBMu03tUgo9iWoD1+7wxByzjBXcikjrqKO1jI5qWqe+svhe/2xH4GTUbXNCkuSI64eVELc5PGHH\ne+xu7D3agm0HzkCW5aCNMeLZtUl7ERUa9OHdtrFPWo2tVtgdbmw7cAaSLIe0nAOvfabdhne/qMXh\nU53qsIZOFz6ObLG7UPpmBf7PKztxpt0WdNKMZ7Z4JErLudA/oclic/ln1cswpSQh25yCLqszam9O\nrGV2ajinGiPOto9Wf42tNryz5TgWLi8P+h/xShKsvW64PV48/uJ2PP32vqDlR/HYfbgF3VYXuq2u\nqJsTaS/QAOCa8aMAAP8551osfOgG3H1TEVweCZ/6t2EdqxlvVowbY4bLLeE//vAFvoyyl4KWGs45\n2nAenpbz6VZb0IXlnhhzQgCgbMcp9Do9aO7sxe9XV6JZM9Nc2fNhdK5pwOVhONOAaNefalvRDz+9\nJeh5kiRjsWbCkdcrR21p6/W+seSfLP0CP39hG17+4GDU9dixxGrxaluUb248GhQ8sizD6fLG3IAj\nVKQurWhhGDpDW3metqXQZXUF3QtYS+kWB4BNexrDji97rxqvfXQYlcfagkLL5vCEjb2H0p4Qw8M5\n+O+xOTzqLQK1CnJ8J6KuHhde+fshvPbRYWze2xj0fiut6JZOO3756i5s3X8m6DV0Op0afKn+YOu2\nudSTZl2tcllOAAAVZUlEQVSTJSicK4+34d9+/zn2H+/fNqZKOCtdsz12t/p7feGcDI9XjtrKjbb5\nytl2e6Bb25Sk/g1a0S7etGO62vW+azYfx3/96UscqAv8L8S7ra7ilOb5h09Fbnkrk77uurEQ104Y\nhduuvQiArz6KRmfg2gm+sFaGakK7tYFAoLs8EpZ/ciRoLkMkLV2BcM7xd4sPVzhX+OdLzLr1UgDB\n/6/dNlfQucPtkbD9wBnkZaXi/pkTYe11Y83mwFyUZn93/eicgXVpAwxnGoB3t9bhD6sr1cexZk6H\nbu9YGeMkqtfpcLrNFtTSauvuhccrYfHKPVF/LlSs1pT2tXcebA66U9Jnexvxkz98gaffqgAQ3xK0\nSCfpaC2j0G4y5STd1hU4Gb21sQb/++UdQc8bnev7B+/o40SnaGi1hvQQuNVJatF0WBxqKzH0oiR0\nzNna61Zn5mqNzfO1pOxOD474A6AupIWn1H+0nee8Xkn9zFyc5zv5a1cinG61warpyt57tBVuj4TV\nmhNjqLLtJ/HbN/cGvffdNicMeh0uykvzP3apvzctxYDsDN/kpk/LG7Dkzb1hLVXlvcs2B9dDp9Wp\nTggzm4xBFzFKKzpa4Ctrh4FAaAGBC7G/bzupfq+pLfJ64IqaVix4bbe66YlCCYus9GQcbeiK+D9b\n39wDHYDv3VaMn35/CnIzU4OOF1+Uqc481ut0uDgvPHhuvHI0Hiu5DlP8QV5dF3stektnLwx6HXIz\nUmHQ65GbmRI2FNYXi80VVF8DIcsyyo+0wKDX4c4bCmFM0qvDbG1dvXjixe1Bu7udOGOByyNhymV5\n+ObUsZh4SRYqj7epn5MWf28Yw5nOm5ZOO8q2n0S95ir/TD929Ym1Jnb97vqwIPB4ZXx14Axqo+yU\nFUm0IKpp6Ap7He2Y5Uc7fFtFKjOknXG0nCONS0UKCp3Ot5uWdgjA5u9W/6wivBWsdZEazvG1KLqs\nrqAgsjk8MffATk9Ngssj4Wh9F55dux8d/ha62WQEELgAcbq9+OM7+2HtdSMvOzXsdZSLiF6HBwaD\n7ySuBJVCHfOPMq5o7Q0MNVwy2hfOdU2BoFm/ux4HT4b3pljtrohd0F1WJ97dWofjjd2o1NwkpNvq\nQmZ6MrLSfSHcY3cFZlmnGtWZx2XbT+JYYzd2HWqGxe5SJzApn5uLRgV37XZbXeqs/HSTMajlPMZ/\noo42iU3bLdrqDyjtemxt6zfaxhYf7zyFxlYrVmwI3ia0pbMXGWlGTJ2YD6fLi1Nng1vesiyjvtmK\ngtw0pCaH94oAviGqSy/KAOD7PEe7Ac3Vl+binunjAQSvSIikpdOO/GwT9P7QvzgvHRabq8/xXm25\nf7+6Ek++vCNskml/nGruwelWG667PA9mkxGX5Kejqc0Gj1fCh9tOwuOVsPNgs3rTF2Xcf1JhNnQ6\nHe66qQhAYC8I5WKogN3adL48+ef+3dyiP3YdasYHIff8fWtjDU72c3wt2nhg6VsVWL87/PZ+itDu\nxnhaztY4JyTlZaVCRvAaULvDjQ3+8kRaC3ndZXm4f+ZEdWxP260dyWVjswAAuw6dxZZ9gXF5u8ON\n3hjjzkWjfSfc363ah6radqzZ4ru4UAJKCec9R1rULs1Lx2SiqCC4W1Pp1rY73TDofaeWRn9rMMkf\n1k6XF9Un2lFdF3m4orHVimfX+mb4T/D/PSdCdqCLNHHQ5vBEnF2tHRbZdagZlcfb8Lu3K9DW7UBm\nejIy03wXIBa7Ww3CjLRk5Pj/dkVLZy9eePcAnvzzThw62aEG7MUh4eyVZDU401ONMGmCrtB/sRFp\nwxsgeJay0nqMdOGr1+mCwvlshx2/eGk7vqo6o4bu2Y7ewLI8r4S2LgdG56ThynE5vnoJucBp73bA\n7vSEvaeh/unWYgDALZPHxHze2Px0JCfp1ZUbkVjsLtgcHvVzAwTqM/QGGNFWTtQ2WdTP2AcD3EZY\nlmW8/6XvZ2+7xteVX1iQAY9XRkOLFTs1k0r3+sNXmT8y0T9Z7spxOTDodernurmzF6aUJGT4L3AH\nguFMQok0xhw6LtmX0A064hW6hCieOy71tUOZYpS/i/D5dQfU79kcHlTUtEEH4PmfTg/7mSsvzcG3\npl2ido1GO6krsjNSkJykD+s5OHiyE4veCF7v/o2vjVW/LhodfEJW1vfmZCjh7KsHZQzuynE5uOvG\nQvzygWn4j3smqz832r9O1ebwqBc2yphjYUGG+jf8Yc3+qOOeWqOyTGp4AsAjs66K+fw9R1qweOWe\noBn62nHDA3XtePmDahyp951Ys9KTg7acVFq8GWlGdaKY4nB9p9pVvPqz4zjtD4SLNF27Y/w9B6db\nbcgyJ8OYpA9qOY/zXwS1+d/HNz89iide3K62xtstDrWbXOmmbQ7ZzjI12YBJRdlotzjV9+XdrXVo\n63bg9Y8Pq/M/PF5JHQ5otzggyTIKcky48tIcJBl02H24JainQekJC/0shLp2wij87t9vwf13TIz5\nvCSDHpeOycDpNmvUi+Wj/vdBuQgDfC1nAGjyX3x4vBKe+1sV/u33X+DZtfvDLqB3HvTdKMOg16Gp\nzRbxf6Szxxl1/kiP3YXfr65EVW07rro0B9f4u+PHjfG9V5+WN8DjldTx9uq6Dni8Eo6f7sbYvHRk\n+td5pyYnYWJhNk4196Ctuxdn2+0Ym5c+4JnaQB/hLEkSFixYgJKSEsyfPx/19cGtjs2bN2Pu3Lko\nKSnB2rVrB1wIosE0GBNKXG5vv8e+orn9a2PVzSi0XZOdPU40tlr9J/HwrsQ0fygr4Rxr33MASDHq\n1bFSRaRJOwDwjzcXqV8roaFQup5zMpKDHje0WmHQ6/DT709BeqoRxiRD0O/LyzYhyaDHscausKVg\nk4tzAQQm3cQj3WRUT9YAkJ9twuTxuWHPG+1vea367BhqT1uwalONeqy2yQJjkh7fv30CvJIcNBY/\nNi8d6SYj9Dodum3OoJZz4Whz0DIo7XBLY6tVXduu7daeWBhYclTgv1BJ1Yw5j8pMRXpqEjosvv2W\nN1ecRrvFgc0Vp9Hr9MDm8OCSfDPMJmMgnDWzmQFfb4ZSJ2fa7fBKEg6eCB7XLb4oE4BvXPR4Y7fa\nCh2dY0J6qq9ru6nNFtQNrEwwC/0sRJKXbYq657vW+LFZkOXAhiahlG2Gr/K35oHADHClt+y9rXXq\nPJWq2nb815++xKLl5ahrskCWZVTVtsOUYsDcb0zwPSdkjHvHwbN4/MVteOq1XWEXOk63V71QvGb8\nKPz4u1er26sqZVIu9GZcdzFyMlJwtL4TJ8/0wOWW1FazQvlsfrKrHpIsq0MAAxUznDdt2gS3243V\nq1fj8ccfR2lpqXrM7XajtLQUy5cvx8qVK7FmzRq0tw98I3qiwbJ+V/Su63j96vXd6pX9uZp/16SY\nXeQ5/lb13TcVBX0/LdV3Yk+LMDM6EoNej4y04G60b99UGPG52olKl2laLlq5Gb5ylW0/hdc+OoTG\nFhvGjEpTb5UIACbNSdpsMsKYpAtruU+ZMArXXuZrecQ7d8CUYsDlhdlBFxe5mSlqV7u2/JP9M4QV\nFrsbn+6ux6t/P4SGFisuHZOBb37tElySn47UZAPunzkRkwqzccf1hdDrdMjLTkVzR6+6xjwjzRfY\n/+/hG/Ho7Mm4+erR6mv/r3++BpM0IVysOQFfMS48nLXlNJuMGJObhuYOO77U9AZtrz6rBtiY3DQU\n5JjQ1tWLHrtLDc1/++7VmDYpH/dML1YD9HhjN06c6Qmr79uu8XU5r9hwFL99c6/aW1PgH/OePX08\njEl6vLWxBu3dvbA73OoEvcI4wjleEy72XSQcj7CPfH1zD7ZXn0VGmjEoxMaNzkBmejIqatpw8qwF\nn5Y3IC8rFct+9nXcc1sxLs5Lx6nmHrzy4UE0tFjR1u3A1ZfmYtrEfAC+3hNFW3cvln98BLLsGyZ4\n5e+H1FtdVh5rw0vvV+NUcw9uu+Yi/PT716qtYF9dmYL+lyYVZmNSUTYsdje27PPNEbmiKHBRAQCT\ni32fQ2Wpp3KRNFAx/+srKiowfbqvu23KlCmorq5Wj9XW1qKoqAgZGb6KnTZtGsrLy/Htb3/7nApE\nJIK+7k1d+uhteHLZV3G/3v0zJ+Kp13ZHPPbobF/X8JxvTMDMGwrx8xe2AfCNWQLhM4K1br5qtDom\nJsty0F2cSr51OW68cjTW72pQx+UU2sDPyzbh6uLcoCEFY5Iet0weo94OdNsBX/dh6MYTmZpZ26nJ\nhqCgKCwwIzM9GXO/MQFJSfGPoCUn6fHMT/4BWeaU4BsrmJPVtbDZ5mS1u/SqS3Pw2d7gSXXaSXnX\nThiFlGQDFjx4AzxeCanJSeq2lIBvnLPyeBs27W2EThfons7LNiEv2xQ0W7ywwIz7Z07Eb1buwezp\n45GanITf/eQWdFtdQS1tZRzfrBlzTDcl4drL8lDbZMH63fUw6HX4h2vGYOv+M3hmlW9TnQljs2B3\nelDXZMF/P+f7fI3JTUPxRZl4dPY1AAJj0+99WaeOK//z18fj3a11MOh1uOmqMVj5aaD3QKGsuR2d\nm4Y5Xx+P1ZuP48FffxpUv5Fm4Q/U5YXZSDLo8NGOk3C4PGhosSIvMxVJSXp8vu80PF4ZD/3jler8\nBMC3nPKGSQX4rKJRvZnM/TMnwpSShO/eVozv3laM1Z8dw6flDerWxFMn5SMv24QJYzNx+FQnahq6\nkJORgr+uPwKPV8Ijs67Cgbp27DzUjEf/uDWoi/uKomzMv2tSWPezTqfD16dcjI92nMJ9d1yOtFQj\nrizKwc6DzdhxsBl6nS7oIg3w7TOeZU5WL/IuHXNuFzoxw9lqtcJsDvxzGAwGSJIEvV4Pq9WqBjMA\npKeno6dnYLcOpHM38/pCdQvFgRiTm6buwCQis8kYNoNzVGYqfvnANNQ1WfC1y/OwYXcD3tkSeUnN\njOsuRlpqEmZeX4hscwoOn+zAH97ZD68k48G7r0BNQxe2V4ff6L0gxxSxezsnMwX/dOul6pjzg3df\ngTc+OaLW4503FOLT8gZ87fI8AMDYfDMenT0Zy96rRkG2CVcV5+LzfafxD9eMUVuHep1ObRkqfx8A\ndc9hwHc1rizXuHJcDn783avR0eNETUMXsszJGJWZihNnLLhj2iW48wZfq/nXD9+IfTWt2FDegJqG\nLhQVmKHT6fDo7GvUMbwZUy4OCucnSr6mBqFWUUjLSru3sk6nw3/cMxl//vAgxoxKw79/72q121e7\nvO0Hd1yOdV/U4dZrxiDDZERDixX7jrUhJyMFmenJuO9bl6u9Bsp4ZGZ6Mgx6Pe66sQhdVhduvXoM\n/rh2P3qdHlxRlINrJ4xCVW077rmtGLVNFhj0Otw6eQx6nR518lKSQa9umap15bgctev08kuyw9Zw\nXz0+V71IUd6LF376dXUSX16Wb+tJ7Ulf2VXLbDJi/l2TcPhUJwpyTJh+7UX4aPtJuDwSrr+iAN+/\n/TIcbehGc4cdxiQ9JhZmw5ikD/osXndZXlB5RmWlBl2UpaUk4a4bC3F1cS7ys01IS01SP38TLs5U\n60M7jj7zhkLYnR5U1XXAYnOiw+LELVfFnuTVX5lpybjpytHYVn0Wn+wM7s3KSk/GA3dNwtf8LV6t\nWf9wKZrabWhqs+E7t4zDlJC//3u3FaP8SAs6e5zIy0rFDVcU+P6m6wvx8gcHUepfCgkAVxfn4qar\nR2PqpHy4PRKqT3RgcnEurhyXgzGj0jDlsryodwqbPX08Zt5QqLaob756NNZsPg6704NvThsbdGEK\n+D7/t1w1But31yMnI0VdvTBQOjnGhsalpaWYMmUK7r77bgDAjBkz8MUXvjscHT16FEuXLsUrr7wC\nAFiyZAmmTZuGO++885wKREREdKGL2dc0depUbN26FQBQWVmJSZMmqcfGjx+PU6dOobu7Gy6XC+Xl\n5bjuuuuGtrREREQXgJgtZ1mWsXDhQhw9ehSAr3V88OBB2O12zJs3D1u2bMGyZcsgSRLmzp2LH/zg\nB+et4ERERIkqZjgTERHR+cdNSIiIiATDcCYiIhIMw5mIiEgwDGciIiLBDEk4c0/uoddXHZeVlWHe\nvHm477778Ktf/Sri7fQotr7qWPHUU09h6dKl57l0iaOveq6qqsL999+PH/zgB/jZz34GlyvyLRcp\nur7qeOPGjZgzZw7mzp2LVatWDVMpE8P+/fsxf/78sO/3O/fkIbBhwwb5ySeflGVZlisrK+Wf/OQn\n6jGXyyXPnDlTtlgsssvlkufMmSO3tbUNRTESWqw67u3tle+44w7Z4XDIsizLP//5z+XPPvtsWMo5\nksWqY8WqVavke++9V166dOn5Ll7CiFXPkiTJ3/ve9+T6+npZlmV5zZo1cm1t7bCUcyTr67N8++23\ny93d3UHnZ+q/V155RZ41a5Z87733Bn1/ILk3JC3nePfkNhqN6p7c1D+x6jglJQVr1qxBSopvK0iP\nx4PU1PCtGCm2WHWsHK+qqsK9997LnolzEKueT5w4gezsbCxfvhzz58+HxWLB+PHjh6uoI1Zfn2Wj\n0QiLxQKn0wlZls/pVocXsnHjxuGFF14IOx8MJPeGJJyj7cmtHOOe3OcuVh3rdDrk5vpuX7Zy5Ur0\n9vbi1ltvHZZyjmSx6rilpQXLli3DggULGMznKFY9d3Z2Yt++ffjhD3+I5cuXY8eOHdi5c+dwFXXE\nilXHAPDQQw9hzpw5mDVrFm6//fag51L87rzzThgM4bfTHEjuDUk4m81m2GyBO7koN8sAgIyMjKBj\nNpsNWVmRb1lH0cWqY+Xx008/jR07duD5558fjiKOeLHqeMOGDejs7MQjjzyCV199FWVlZXj//feH\nq6gjWqx6zs7ORlFREcaPH4+kpCRMnz49rNVHfYtVx01NTXjrrbewefNmbN68Ge3t7Vi/fv1wFTUh\nDST3hiScuSf30ItVxwCwYMECuFwuLFu2TO3epv6JVcfz58/Hu+++i5UrV+LHP/4xZs2ahXvuuWe4\nijqixarnwsJC2O12dQLT3r17cfnllw9LOUeyWHXsdDqh1+uRnJwMvV6P3Nxc9mYOsoHkXnx3ce+n\nmTNnYtu2bSgpKQHg25O7rKxM3ZP7ySefxMMPP6zuyV1QUDAUxUhosep48uTJWLduHa6//no88MAD\nAIAf/ehHuOOOO4azyCNOX59jLY7RDVxf9bx48WI89thjkGUZU6dOxYwZM4a5xCNPX3U8e/ZslJSU\nICUlBePGjcPs2bOHucQjm3I+OJfc497aREREguEmJERERIJhOBMREQmG4UxERCQYhjMREZFgGM5E\nRESCYTgTEREJhuFMREQkmP8PEhGz76MLs7AAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t, y, dy = create_data(100, period=0.3)\n", "period = 1. / freq_grid(t, nyquist_factor=10)\n", "\n", "model = LombScargle().fit(t, y, dy)\n", "power = model.periodogram(period)\n", "plt.plot(period, power)\n", "plt.xlim(0, 1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a data sampling rate of approximately $1$ time unit, we easily find a period of $0.3$ time units. The averaged Nyquist limit clearly does not apply for irregularly-spaced data!\n", "Nevertheless, short of a full analysis of the temporal window function, it remains a useful milepost in estimating the upper limit of frequency." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Scaling with $N$\n", "With these rules in mind, we see that the size of the frequency grid is approximately\n", "$$\n", "N_f = \\frac{f_{max}}{\\delta f} \\propto \\frac{N/(2T)}{1/T} \\propto N\n", "$$\n", "So for $N$ data points, we will require some multiple of $N$ frequencies (with a constant of proportionality typically on order 10) to suitably explore the frequency space.\n", "This is the source of the $N^2$ scaling of the typical periodogram: finding periods in $N$ datapoints requires a grid of $\\sim 10N$ frequencies, and $O[N^2]$ operations.\n", "When $N$ gets very, very large, this becomes a problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fast Periodograms with ``LombScargleFast``\n", "Finally we get to the meat of this discussion.\n", "\n", "In a [1989 paper](http://adsabs.harvard.edu/full/1989ApJ...338..277P), Press and Rybicki proposed a clever method whereby a Fast Fourier Transform is used on a grid *extirpolated* from the original data, such that this problem can be solved in $O[N\\log N]$ time. The ``gatspy`` package contains a pure-Python implementation of this algorithm, and we'll explore it here.\n", "If you're interested in seeing how the algorithm works in Python, check out the code in [the gatspy source](https://github.com/astroML/gatspy/blob/master/gatspy/periodic/lomb_scargle_fast.py).\n", "It's far more readible and understandable than the Fortran source presented in Press *et al.*\n", "\n", "For convenience, the implementation has a ``periodogram_auto`` method which automatically selects a frequency/period range based on an oversampling factor and a nyquist factor:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function periodogram_auto in module gatspy.periodic.modeler:\n", "\n", "periodogram_auto(self, oversampling=5, nyquist_factor=3, return_periods=True)\n", " Compute the periodogram on an automatically-determined grid\n", " \n", " This function uses heuristic arguments to choose a suitable frequency\n", " grid for the data. Note that depending on the data window function,\n", " the model may be sensitive to periodicity at higher frequencies than\n", " this function returns!\n", " \n", " The final number of frequencies will be\n", " Nf = oversampling * nyquist_factor * len(t) / 2\n", " \n", " Parameters\n", " ----------\n", " oversampling : float\n", " the number of samples per approximate peak width\n", " nyquist_factor : float\n", " the highest frequency, in units of the nyquist frequency for points\n", " spread uniformly through the data range.\n", " \n", " Returns\n", " -------\n", " period : ndarray\n", " the grid of periods\n", " power : ndarray\n", " the power at each frequency\n", "\n" ] } ], "source": [ "from gatspy.periodic import LombScargleFast\n", "help(LombScargleFast.periodogram_auto)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFVCAYAAADc5IdQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0lOW9L/DvO/dJJveEcA0EEPCKBm+1UiwVCru2WyrS\noKbaerrO2muf3Z6erd22PWVh9/YQt9u6uytuq+2iLVqwFq9oQREERRAkJBjul0Ag4ZJ7MtfMzPue\nP96Zd26ZmQCZzJPX72ctV5KZXJ6MZL7ze57f87ySoigKiIiIKGsM2R4AERHRFx3DmIiIKMsYxkRE\nRFnGMCYiIsoyhjEREVGWMYyJiIiybFBh3NDQgJqamoTbN2/ejMWLF6O6uhqvvvrqkA+OiIjoi8CU\n7hNefPFFvPXWW8jNzY253e/3o7a2FuvWrYPNZsPSpUsxd+5clJSUZGywREREepS2Mp44cSKeffZZ\nxJ8Ncvz4cVRUVCAvLw9msxmzZs3C7t27MzZQIiIivUobxvPnz4fRaEy43el0Ii8vT/s4NzcXfX19\nQzs6IiKiL4BLbuDKy8uDy+XSPna5XCgoKEj5NTx5k4iIKFHaNeNkJk+ejFOnTqGnpwd2ux27d+/G\nww8/nPJrJElCWxur50wrK8vj45xhfIwzj49x5vExHh5lZXlpP2fQYSxJEgBg/fr1cLvdWLJkCR57\n7DE8/PDDkGUZixcvxqhRoy59tERERF9Q0nBftYmvwjKPr3Yzj49x5vExzjw+xsNjMJUxD/0gIiLK\nMoYxERFRljGMib4AuJOBSGwMYyKdc3r8ePjJLXjz46ZsD4WIkmAYE+nc0TPdAMAwJhIYw5hI54JB\nTlETiY5hTKRzMteLiYTHMCbSOVbGROJjGBPpXECWsz0EIkqDYUykc7LMyphIdAxjIp1jGBOJj2FM\npHMBhjGR8BjGRDoXCHDNmEh0DGMinetnGBMJj2FMpHP9gWC2h0BEaTCMiXSu38/KmEh0DGMinfOH\nKmOjQcrySIgoGYYxkc6F14wZxkTiYhgT6Vx4mtpoZBgTiYphTKRz4QYug8QwJhIVw5hI5/yhytjA\naWoiYTGMiXQuXBnzWEwicTGMiXQu3MDF6xoTiYthTKRz4WlqXteYSFwMYyKd84WmqYOcpiYSFsOY\nSOe0ylhWoHCqmkhIDGMinYu+UATXjYnExDAm0jFZVhAIRsKY68ZEYmIYE+mYP+7yiVw3JhITw5hI\nx+Ivn8gwJhITw5hIx+Ivn8gwJhITw5hIx+IrY57CRSQmhjGRjiWsGQflJJ9JRNnEMCbSsYRpam5t\nIhISw5hIxxIauLi1iUhIDGMiHYtv2GIDF5GYGMZEOhbfsMUGLiIxMYyJdCx8/KVBkgAAAZkNXEQi\nYhgT6Vi4Ejab1T91rhkTiYlhTKRj4Vlps1H9U+c0NZGYGMZEOqZVxqZQZcwwJhISw5hIx8Jrxiaj\numYc5JoxkZAYxkQ6FqmMjQBYGROJimFMpGMJ09Rs4CISEsOYSMfC09ThBi5WxkRiYhgT6ZjWTW1i\nNzWRyBjGRDoWP03NQz+IxMQwJtKx+GlqVsZEYmIYE+kY9xkTjQwMYyId0/YZs5uaSGgMYyId0ypj\ndlMTCY1hTKRj8d3UPIGLSEwMYyIdU7hmTDQiMIyJdCwos5uaaCRIGcayLGPZsmWorq5GTU0Nmpub\nY+5///33cc8992Dx4sVYs2ZNRgdKRBdP29rEyphIaKZUd27atAl+vx9r165FQ0MDamtr8dxzz2n3\nr1ixAm+88Qbsdju+8Y1v4K677kJeXl7GB01EgxO5ahO7qYlEljKM6+rqMHv2bADAzJkz0djYGHO/\n2WxGb28vDAYDFEWBJEmZGykRXTQl1K/FyphIbCnD2Ol0wuFwaB8bjUbIsgyDQf3D/t73vod77rkH\ndrsd8+fPj/ncZMrKWDkPBz7OmTcSHmOrzQwAKC7KAQBYrKYRMe6wkTTWkYqPsRhShrHD4YDL5dI+\njg7i1tZWvPzyy9i8eTPsdjseffRRbNiwAQsWLEj5A9va+oZg2JRKWVkeH+cMGymPsdPlAwB4Pf0A\nAJfLNyLGDYycx3gk42M8PAbzgidlA1dVVRW2bdsGAKivr8f06dO1+3w+HwwGAywWCwwGA4qLi9HX\nx/+pRCKJP5s6wGlqIiGlrIznzZuH7du3o7q6GoDasLV+/Xq43W4sWbIEixYtQnV1NaxWKyZOnIhF\nixYNy6CJaHAS9hmzgYtISCnDWJIkPP744zG3VVZWau8/9NBDeOihhzIyMCK6fPHd1OGPiUgsPPSD\nSMfk+G7qII/DJBIRw5hIx4I8DpNoRGAYE+mYwhO4iEYEhjGRjsV3UzOMicTEMCbSsfCFIUxcMyYS\nGsOYSMfChXCkmzqLgyGipBjGRDoWrowNkgRJiqwhE5FYGMZEOhZeMzYY1EDmPmMiMTGMiXQstjKW\nwCwmEhPDmEjHIpWxBIMUCWciEgvDmEjHtMrYIEEysDImEhXDmEjHoqepDRLPpiYSFcOYSMdkRQ1i\nQH3LbmoiMTGMiXRMVhQYQn/lkiRxnzGRoBjGRDomy0pUZcx9xkSiYhgT6ZhaGathLEkSu6mJBMUw\nJtKxmMqY3dREwmIYE+mYrCCqMmY3NZGoGMZEOqZWxur77KYmEhfDmEjHZEWBFFMZZ3lARDQghjGR\njsmyAqMhss+YDVxEYmIYE+mYrEQauCROUxMJi2FMpGPx3dQsjInExDAm0jFZgbZmzEM/iMTFMCbS\nsehuavU4TIYxkYgYxkQ6pihxDVzMYiIhMYyJdCwoRzdwAQrTmEhIDGMiHYveZ8wGLiJxMYyJdEyW\nEVsZc82YSEgMYyIdU6KuZ2xgAxeRsBjGRDomywqMMYd+sDomEhHDmEinZEWBgshVm8JbnBjFROJh\nGBPpVPgcaimqMo6+nYjEwTAm0qnwdLQhqptavT1rQyKiJBjGRDoly+rb6G5qAGziIhIQw5hIp8Kh\nG30CF8AGLiIRMYyJdCqorRmrHxu0NeNsjYiIkmEYE+mUHLdmzGlqInExjIl0KnwOtXY9Y05TEwmL\nYUykU+EdTFplHHrLnU1E4mEYE+mUnFAZq7ezMiYSD8OYSKeC2pqx+rGBh34QCYthTKRT8WvGksRD\nP4hExTAm0qn4bmoDu6mJhMUwJtKp+DVjycBuaiJRMYyJdCq+mzpSGWdpQESUFMOYSKcSu6lZGROJ\nimFMpFNyXDc1L6FIJC6GMZFOJa+MszYkIkqCYUykUzybmmjkYBgT6VRCZWxgZUwkKoYxkU5pYczK\nmEh4DGMindK2NsVfz5hhTCQcU6o7ZVnG8uXLceTIEZjNZjzxxBOoqKjQ7t+3bx+efPJJKIqC8vJy\nPPnkk7BYLBkfNBGll7hmHJqmlrM2JCJKImVlvGnTJvj9fqxduxaPPPIIamtrtfsURcGyZctQW1uL\nP//5z/jSl76EM2fOZHzARDQ4iWvGodtZGRMJJ2VlXFdXh9mzZwMAZs6cicbGRu2+pqYmFBYWYtWq\nVTh69CjmzJmDyZMnZ3a0RDRoCZUxeOgHkahSVsZOpxMOh0P72Gg0QpbVOa6uri7s3bsXDzzwAFat\nWoUdO3Zg586dmR0tEQ1asm5qHvpBJJ6UlbHD4YDL5dI+lmUZhtBcV2FhISoqKrRqePbs2WhsbMSt\nt96a8geWleVd7phpEPg4Z57oj7GjtRcAkJ9vQ1lZHhwOq/pxgV34sYeNlHGOZHyMxZAyjKuqqrBl\nyxYsXLgQ9fX1mD59unbfhAkT4Ha70dzcjIqKCuzZsweLFy9O+wPb2vouf9SUUllZHh/nDBsJj3F3\ntwcA4Hb50NbWB4+7HwDQ1eUWfuzAyHiMRzo+xsNjMC94UobxvHnzsH37dlRXVwMAVqxYgfXr18Pt\ndmPJkiV44okn8M///M9QFAVVVVWYM2fO0IyciC5beM04fOlEbZqas9REwkkZxpIk4fHHH4+5rbKy\nUnv/1ltvxauvvpqZkRHRZQmvDRslHvpBJDoe+kGkU/EncPESikTiYhgT6ZS2tUmKPfRD5qEfRMJh\nGBPpVHhtWAr9lYePxWRlTCQehjGRTsXvM5Z4NjWRsBjGRDqlKLyEItFIwTAm0iltmprd1ETCYxgT\n6VSkm1r9mN3UROJiGBPpVPw0tVYZs5uaSDgMYyKdir9qEytjInExjIl0KuGqTeymJhIWw5hIp8IN\nXOH9xeH9xsxiIvEwjIl0KtlxmKyMicTDMCbSKe2qTfHT1LxsE5FwGMZEOhXfwCVpx2Fma0RElAzD\nmEinFO3QD/Utp6mJxMUwJtKphLOpeRwmkbAYxkQ6FX8JRQOPwyQSFsOYSKeU0ElbPPSDSHwMYyKd\nilTG6scSu6mJhMUwJtKpxOMw1dtZGBOJh2FMpFMJDVzspiYSFsOYSKe0Qz/ClbEhHMZZGxIRJcEw\nJtKp8KUSI2vG6ls2cBGJh2FMpFPx1zPmoR9E4mIYE+lU/NnU4bfhLU9EJA6GMZFOaZdQ1NaMw7ez\nMiYSDcOYSKcUOW6fMThNTSQqhjGRTiXsM+bZ1ETCYhgT6VTiPuPQ7UxjIuEwjIl0SlszjuumVrjR\nmEg4DGMinYpMU6sfRyrjLA2IiJJiGBPpVHiaOrylKXICF9OYSDQMYyKdUhIuFMFLKBKJimFMpFPh\n6ejQ7HRkmpqHfhAJh2FMpFOyokCSoqapw5UxWBkTiYZhTKRTiqxoAQxEXUKRHVxEwmEYE+mUrCja\nejHAQz+IRMYwJtIpWUZcZRy6nWlMJByGMZFOqZVx5ONIN3WWBkRESTGMiXRKVmLXjA1aNzXTmEg0\nDGMinZJlRWvaAqKuZ8zSmEg4DGMinZIVxDZwhbupmcVEwmEYE+mUurUp8nF4/ZiVMZF4GMZEOqUe\n+jHAPmOGMZFwGMZEOqUkNHCxm5pIVAxjIp1S14wjH0vspiYSFsOYSKfkAY7DlMA1YyIRMYyJdCr+\nOExADWQWxkTiYRgT6VR8ZQyo09Zs4CISD8OYSKdkBTHd1ID6MaepicTDMCbSqfizqQG1o1qWszMe\nIkqOYUykU/HXMwbUaWpWxkTiYRgT6VT8cZgAIEHimjGRgFKGsSzLWLZsGaqrq1FTU4Pm5uYBP+8X\nv/gFnn766YwMkIgujaIoiCuMYTBIPPSDSEApw3jTpk3w+/1Yu3YtHnnkEdTW1iZ8ztq1a3H06NGE\nRhEiyq6Buqklid3URCJKGcZ1dXWYPXs2AGDmzJlobGxMuH/fvn34zne+w3UoIoEoigIFSFwz5j5j\nIiGlDGOn0wmHw6F9bDQaIYdaMS9cuICVK1di2bJlDGIiwYSr38RDP9TGLiISiynVnQ6HAy6XS/tY\nlmUYQnslNm7ciK6uLvzgBz9Ae3s7vF4vpkyZgrvvvjvlDywryxuCYVM6fJwzT+TH2B8IAgCsVlPM\nOE0mIySDJPTYo42UcY5kfIzFkDKMq6qqsGXLFixcuBD19fWYPn26dl9NTQ1qamoAAK+//jpOnDiR\nNogBoK2t7zKHTOmUleXxcc4w0R9jn18N40AgGDNORVYQkGWhxx4m+mOsB3yMh8dgXvCkDON58+Zh\n+/btqK6uBgCsWLEC69evh9vtxpIlS2I+lw1cROIIX5lpoH3GoaKZiASSMowlScLjjz8ec1tlZWXC\n5y1atGhoR0VElyXcx5HYTS1BVngEF5FoeOgHkQ6Fe7TiG7iMBokNXEQCYhgT6VC4mzrh0A9JQpBh\nTCQchjGRDilJ1ox5PWMiMTGMiXQo1TS1zDQmEg7DmEiHIt3UsbcbDDwOk0hEDGMiHZKTdFMbWBkT\nCYlhTKRDWgOXYYCzqRnGRMJhGBPpUNJDPyQJCsDz5IkEwzAm0qFkDVzhj7luTCQWhjGRDilJG7hC\nYcypaiKhMIyJdChZA5dRC+NhHxIRpcAwJtKh8Cx0/AVcwuHMU7iIxMIwJtIhrTKO+wsPZzPXjInE\nwjAm0qFk3dRGNnARCYlhTKRDkco4STc1p6mJhMIwJtKhcNgmWzNmGBOJhWFMpEPaPmNubSIaERjG\nRDqUdJpa4poxkYgYxkQ6lOx6xpETuIZ9SESUAsOYSIfSNXBxnzGRWBjGRDoka4d+xN4ezmaFYUwk\nFIYxkQ6ln6ZmGBOJhGFMpEPJzqbmcZhEYmIYE+lQskso8gQuIjExjIl0SOYlFIlGFIYxkQ6FK18p\nLo0lnsBFJCSGMZEOpb9QxLAPiYhSYBgT6VDyBq7Q/UxjIqEwjIl0SNEauGJv59YmIjExjIl0KGll\nzBO4iITEMCbSISXNJRR5AheRWBjGRDqUbJ8xp6mJxMQwJtKhpPuMeQIXkZAYxkQ6lG7NmJUxkVgY\nxkQ6lOzQDyNP4CISEsOYSIeSHfohafuMh3tERJQKw5hIh+Qk+4x5oQgiMTGMiXQo6fWMeTY1kZAY\nxkQ6FK58jdzaRDQiMIyJdCi8dSlhnzErYyIhMYyJdChpGLObmkhIDGMiHQqHrSmug8vASygSCYlh\nTKRD6aapg9zbRCQUhjGRDiWfplbfsjImEgvDmEiH5FDla4o/gYtXbSISEsOYSIeCwdQNXLxQBJFY\nGMZEOhRMss84fH1j7jMmEgvDmEiHwt3U8WHMC0UQiYlhTKRD6aapWRkTiYVhTKRDwSSVceQErmEf\nEhGlwDAm0qHI2dRJDv3gNDWRUBjGRDoUDKqlb+KhH+pbTlMTiYVhTKRDSaepubWJSEgMYyIdkmUF\nEtjARTRSmFLdKcsyli9fjiNHjsBsNuOJJ55ARUWFdv/69evxpz/9CUajEdOmTcPy5cu1fYxElD1B\nWUkIYoAncBGJKmVlvGnTJvj9fqxduxaPPPIIamtrtfu8Xi9+/etfY/Xq1VizZg2cTie2bNmS8QET\nUXpBWYHRmBjG2jQ1K2MioaQM47q6OsyePRsAMHPmTDQ2Nmr3Wa1WvPLKK7BarQCAQCAAm82WwaES\n0WDJspKwXgxEncDFyphIKCmnqZ1OJxwOh/ax0WiELMswGAyQJAnFxcUAgNWrV8Pj8eC2225L+wPL\nyvIuc8g0GHycM0/ox9ggwWQ0JIzRbLOoby0msccfMhLGONLxMRZDyjB2OBxwuVzax+Egjv74qaee\nwqlTp/Cb3/xmUD+wra3vEodKg1VWlsfHOcNEf4z7+4OQJClhjE6PHwDg8fiFHj8g/mOsB3yMh8dg\nXvCknKauqqrCtm3bAAD19fWYPn16zP3Lli1Df38/Vq5cqU1XE1H2JZum1vYZc5qaSCgpK+N58+Zh\n+/btqK6uBgCsWLEC69evh9vtxjXXXIN169bhxhtvxHe/+10AwIMPPog777wz86MmopSCsqwdfRmN\nW5uIxJQyjCVJwuOPPx5zW2Vlpfb+wYMHMzMqIrosQVmBxZw48WVgAxeRkHjoB5EOBZNNU7MyJhIS\nw5hIh+Qkh37wQhFEYmIYE+lQ0sqY09REQmIYE+lQsjAG1ItH8AQuIrEwjIl0KNk0NaA2ZsryMA+I\niFJiGBPpjKIoocp44D9vo0FiAxeRYBjGRDoTDtpk09QGA9eMiUTDMCbSmXDQJpumNkisjIlEwzAm\n0pmgnK4yllgZEwmGYUykM2nDWGIYE4mGYUykM8F009Rs4CISDsOYSGdkVsZEIw7DmEhngsFBdFMz\ni4mEwjAm0pnw6VrJp6kNrIyJBMMwJtKZyDT1wH/eBimyrkxEYmAYZ1Fbtwf/sXYvWttd2R4K6Ugw\nqJ51mWprk8IGLiKhMIyzaO0HR3HgZBdWvXsw20MhHUm3tckoSayMiQTDMM6i8BNigE+MNITkNGvG\nErc2EQmHYSwCPi/SEErXTW008KpNRKJhGGdR+KlSYRrTENKmqY3cZ0w0UjCMs0iSQk+WfF6kIaSd\nwCUlC2N1KptNXETiYBgLgE+JNJTSnsAVup1ZTCQOhnEWaYUxnxRpCEWmqZPsMw6FMZu4iMTBMM6w\nQFDm+hwNq2CoOyvpNHU4jPnvkkgYDOMM+8dntuFfnt+R5rP4pEhDZzAXigB4CheRSBjGGeYPyOjo\n9QIAPmpoxf96Zht6XP0AIg1cfEqkoZSumzoc0gxjInEwjIfRqr8dgtsXQMOxdgCRNWOmMQ2ldNcz\nNpvUP3t/gJuNiUTBMM4Cq9kIIHqfMdHQ0aapk6wZW0L//voDwWEbExGlxjAeJn/ZfEx7X3uODE9T\ns6uVhlC6aWqrKRTGflbGRKJgGGdQdMhu2NWccP/AT5VElyfdNLXFrP7ZszImEgfDOIPSNcgkmUUk\nuizprmesTVOzMiYSBsM4gwbbrcpZahpK6S6haAk1cPX71cq4td2F5at24cW39w/PAIkoAcM4g8JX\nz0mHWUxDSTv0I+k0dbiBS4bHF8CKl/ag+bwTO/afH7YxElEshnEGBdNcp04agvMw9xxuw4VuzyV/\nPelPukM/oivjtm4PXN6Adp+vn+vIRNnAMM6gtGvGobeXGsUXuj1Y+frn+PkLOy/xO8Tqdffj6bV7\ncby1Z0i+36XYd7wDK1/7nHtgL0Paaeqoytjp8cfc1x46oIaIhhfDOIOSTVNrhfBlpnF3n0/9OUN0\nktK6D49j/8kurHr30JB8v4ulKAr+89UG7DnShuYLfVkZw0h19Ew39h1XD5MJ/7sbTGUcDuOCXAsA\noKOHYUyUDQzjDEo7TR16q1xiGocbcIZKZ6gqCh9KMty8UVOknqipU0pNURSseKkO//nqPiiKol2N\nKe2asT+IPrcaxpNG5wGAdnQrEQ0vhnEGJatYI+F7eXubMnW2cJLn8IyLnjKNnz6l5KJ7Bvo8/qhp\n6mRbm8L7jCPT1BPDYczKmCgrTNkegJ4lDeO4gvlS+7eGPIyzfOEKhvGlOXSqS3u/s9c7iGnqyD5j\nr0+djZg0Oh8AK2OibGFlnEHJ1ox73epVm3bsP3dZ3/9ympw+O3QBjz73iTY1LYLwlCnAML4YB6PC\nuKPHl76bOuoErj6P+m9xwigHDJLEypgoSxjGGZSscn1l8zH0OH3axwNVxrKs4J0dJ3GmzZn0+19O\nGP/2rf3o6PXivd2ntdu0NewslcbOUDAAgMvDNePBajrbq73f2edFIN0+Y1PimnF+rgVFeRZWxkRZ\nwjDOoFQNXGc73FEfJabfZ4cvYN3WE/jLlmMJ94X5g4nfPxCU8cr7h9HS7ko5tnD1FL6cXuwospPG\nzqjK2NsfCeOP9rXiqTV74R6Cpi631482He3LVhQF3c7Ii5jOXq/2OOXYBl6F0ipjv7pmbLMYYTYZ\nUJJvQ3efD4EB/l0RUWYxjDMo1QlcspI6+sIVS+OJzqTfY6DKeGt9K17acAir3j044Nccbu7CE6s/\n035mdBUcDD0JZ6sy7ouamvZGdYqv2XQUB0914W+fnrqs7+/y+vFPv/4Iy1ft0k3geHxB+ANypAGr\n1we3V30cc5OGsVoZ+wLq1iaH3QwAKCmwQQHQ2ecb8OuIKHMYxhm0vfFs0vt8UWFzqeHnH+CqO92h\n6e+TZwfep/vW9pM43hKZ1vT1B+H2+vHUmr041NwNAFkLquh14uiToMIPT3eKkHhvVzN++tsdaGzq\nSPo59UfboShqgA1l4MiygjNtzgH/f2Raj0v9PSaMcsBokNDZ64XTG4DFZIDZNPAWNbMpUhn3uf3I\ny4mEMcCOaqJsYBhniNvrx/bPkzdoJTt2UJYV+PqDg7qi00CVsXbCZpKpZldcY5TPH8SR0z0xTUDZ\nuppPzDR1zIsV9XfpS9LUpSgKXtlyDOe7PNhS15L0+5/vikxPD+VU9YqX9mDZ73fhjY+bhux7AsDB\nk5341z9+htXvHU76OT2hKeoihxXF+VZ09Hrh9vqTTlEDgEGSYDYZ0OfuRyAow2FXD/woyWcYE2UL\nwzhD3L7U65sboxqnor38/hH8w6+2oiuucpNlBQdOdmprvcDAa8bhNqxk1Xb81/j8Qbi8cQEdCGJr\nfQv2HG5L+TsMtT6PHxIAq8WovVgJyjL8oRcHfe7I2qjT48fL7x3Bh3tb0Nru0n7f/Sc7k1b20QHc\nPkRh3Ovux/FWdabhQFNXms++OL975yCazvZiS11LzP/3aD0u9TEpcFhQnGdDj7MfvW4/ckNTz8lY\nTAZ09qr/xqKnqQFubyLKBoZxhqTrdD51LjKNrEQl55a9amV36nzsNPOaTUfxH2vr8emByJV1/ANU\nsOlO/fLEvUjw+YPaVqswl8ePP244jJWvfx4ztnidvV5sqTsTExTbPz+LAyeTr3MD6u+7o/FcwrYq\np0cNkRyrSQvjXpdfq/Gjtz69t7sZH9SdwZ82HkbdkciLhn6/jL992oxfvVKf8LteiKmMhyZwopcD\nTl9wxjSeXY6uPl/MC7LzXe4BP08L41wLChxqhevrDyLXmvoIAYvZqC2VaNPUrIyJsoZhnCEXs092\noLgLRIW5rz+I3YfUED7UHKm+oqvccGj29ycPY0VR4IzbMtTvD6LPFTvWQFTjWW/oyf5shyvh+M1l\nv9+F1e8dwcHQmLr6fPj9OwfxH2vrYxrU4tUfbceL6w/g2dc+j7nd6e6Hw26G1WzUpqmjq+FwGAdl\nGR/vi6zHv/tpMwDghitKAQCvbzuBxqbOmBcugFoZh9dLh6r6O3lOrYrHlORAVhQ0tfam+YpBft/Q\ndqVwE9aZtoG748NrxgW5VhQ6rNrtaSvjqCNPtco4P7OVsc8fxEcNrQx7ogEwjDMkPmhS6XH2Y2fc\nASC9UVXg06/Ua++HL7u492gbPjt0Qbs9HKC+FOdV9/vlhCncgSrjaH0eP3buP4efv/gp1m09EXNf\neCr+XGibVuOJSPPU+U43FEVJqE4B4MgZtVHsZNTsgBx6oeDIMcdMU0evcfv8QfT7g/hwbyu6nf24\nurJYvT30ubOml8X8nP1NkQrd7Q3A6fFj6rgCAIjZ5x32p42HsfaDoylnA+KFZzjuuH4cALU6HgpN\noZCffd1YAMCZqO/b3uPRXqSE14zzHZHKGEi+rSnMErWlLfx1FrMR+TnmjITlwZOd+NkLO7Hqb4dQ\n+3JdwjIM0RcdwzhDoqdUB+OFtw/gXGdkKjI6LI619MBoVP9Xhae/f7Pu85jr0AaCMjZ82oyPP0/e\nwT1Qte4A20K5AAAWvklEQVTtTx3Gf/jbIbyzQ91StLU+0hwVvc4cnv491hK59GJ7jxfv7jyFH/3X\nRzFT8oC6/SaexxeArCjIs5thC02h7jp4Hrvj1q23NbTi5fePwGox4oH502CzqBWeySjhxumjYj73\nXNTUbni9eExJDvJzzAlhcLbDhQ/3tuC93aexO+pFTjy3N4DGEx3agS5nO9zItZlw5cQiAEBrR+r9\n3YPVFJr+vv26MQCgHf7S5+7HT/57B55eq75Ai56mLsyNqoxt6SrjyJ/+6OIc7f2SAhs6+7wpZzYu\nlscXwAtvH0Cvqx/XTi5BR68Xz/ylQfjLZCqKgvYej262wZHYeDZ1BlzqE1lz1DqxK+6Ai3B4OD1+\nPPlyXcLX+oNyygNCwl8bb6Bp6mgnoqZdo08Ui26GuhAKvegwbuv2aJX0ps9OY/FXpyI/xwxJkmLW\nivv9QVjMRq2T2mE3a2vQz7+5X/u8/Bwzet1+vBnqWP7hPdehvCgH48pycbylFwW5FljMRnytajw+\nPXgeJqOE9h4vvP0BrP3gqPbiaFShHYUOK853eaAoijbTsOtgJIDrjrTh5ivLAajV/t5j7Vjy1anw\nB2T8y/M74PEF4A4oqJpSjLZuDyaNzkN5cQ4MkoTW9oHXdtPZd7wDr287ge/MnYorJhTgeEsPyotz\nMLY0F3k5Zq3i3rhLbfxrvuDE+S43epz9sFuNsJqNMZVxsj3GYZaobU/lRVFhnG9D09k+9Dj7UZRn\nHehLL9q7O0+hx9WPv7+9Et/68iSs+tshfLzvLD6sb8G8Gydc9veXZUWbpXGkmZ4fDF9/EH/dehx1\nR9rQ1efD2NJc/NO3r0V51IsWGjkURcGp831oOtuH5vPqf+c7PcixmVDosGr9FuVFORhdkoPRxTko\nybclPcEuUxjGGTBQY9VgRIdPMr3u/oRKE0jcsgSo1bLJGKmAnN6BK2NZSV4ZR5NDl+czSFJMM9SF\nbg+cHj/OdrhhDVW19cfatfu3N57D9sZzmFs1Dg/Mnx4zDbq98Rx6Xf2YPqEQAODIMQ/4YmZ0SS56\n3d1weQMYW5qrVaJTxhbgeEsvrhivfv3SO6/AvV+dguff3I/6Y+3Y1nAW2xoiswVlRXYU5lnRfMEJ\nb38Q9lCj07HQ1LnZZMCh5m4oioK2bg9+9ZcGAEChw4pAQNam3bftPYOK0hwEZQWji3NgNhkwqsge\n6uyOhPxAWtpdeHfHKUwanYd5N01AUJaxeuMhdPT68NSavfjBN6+Ctz+IWyrU32l8mQMHT3XB5fVj\ny94z2vfZdfACOnu92lpxQdSacU6ayjj6MpnhBi4gtqP6YsL4XKcbTrcfY0tzY6bI27o92LjrNIry\nrFhwSwUkScK9d0zBZ4cuYP0nJzH7ujGwWS7taUhRFHzSeA5//fC4NkNw85Wj8L2FV8JqubTLgDo9\nfvznqw040doLh92MGRWFONTcjV/+cTf+4e5rcE1lySV934sRlGW093jh8QXg6w+GXpBZ0n/hCNfV\n58Pxlh7YrEbk2szIsZlQWmBLevWxwThyuht/3Xocx85ECgWjQcKoIju8/UE0ne0d8Nhik9GA0cV2\njC6OBHR5cQ5GFdrhsJtT/n1fKoZxBqRat71c4YaqeAOt87W2u9BwrB1XjC9E5Zj8AQPb5Q0kVOHJ\nKIq6RlnosODo6dgq+HioKr79ujH4YM+ZAU8O21rfinvvmKo9cQLA6o3qHtrKMeoJUvk5loSrWgHq\nVOqR02pgji3N1W5fNHsycm0mzJ01HoB6HrPFYERpKFS2NbTGfJ8rxheiIfRCoavPB7vVBEVRcPJc\nH8oKbZgytgA7D5xHa4cbuw9GGsDe+OgEJEhw2M0oK7ThQFOnNp7RJTnauM51utHV50NxqBkqXr8/\niKfX7kW3sx879p/DtVNK0Hy+Dx29Pu2FzAtvHwAATI8L4x2N5+DxBXHzlaNQd6QNmz47DbcvgJlT\n1YAojK6M7an/tE1G9cmkKM8a88QSnrI+da4PU8cV4IM9Z/D+bjVMq792hXbSV7S3tjfhjY/UGYv8\nXAt+en+VVkW+uuUYAkEZ994xRXsBkJdjwddvrsCbHzdhw6fNuHv25JRj7erz4aOGVpzvcsNoNKA4\nzwqzyYCdB86jpc0Fi8mA66eWoqPXi10HL+Bchxs/XjIz5sXJYLR1e/Drv+5Da7sLt10zGg8tnAGT\n0YBPGs/ijxsO49nXPsfPHpiFivLEx2AotHd78NG+s/hoX2vMEaeSBEyfUIibryrH7deOiXmBrQcn\nz/Wqy0MHLyQEY0m+DfNvnoCvXDf2ol5gnWlzYt2Hx9FwXO1juX5qKW6YVoqJ5XkYW5qrPYZqr4of\n3X0+nOt0q/91uHE29Hagxkm71YRRRXaMKrSrb4vsKC/KQVmhHXk55kv+/5PyL1aWZSxfvhxHjhyB\n2WzGE088gYqKCu3+zZs347nnnoPJZMI999yDe++995IGoTepwnjS6LyYxqWLlazxZf0nJxNuW75q\nt/a+3WrC/JuSTwkW5VljvneycXb0eLHhU3VbEQBcXVmM/U2d+OywOs17TWUxPjt8QWsseuz+KjQc\na0f9sXac7XBj79GB9y6H10jHleYO+Eq1ckwetqlFKsaWRKYLrRYjvvnlyoTPH1VkB6C+IHHYzbjr\ntkmYOq4gFKbqfWc73Bhbmov2Hi9c3gCumlSMGROLsPPAeRw61YVdBy/AYjJgbtV4bNjVDAUKvn7z\nBG28m0OPwZiSXG2MdUfacKylBzfHhbGiqFOpW+vVJ9qSfBs6er14e/tJ9IdO7vpZzSw89/rn2uEk\n0yeo1f/4Uer3D297u25KCVweP/afVLvYp4Sa0nKsJpiMBgSCcto14/CSRXFc9Rtuivv8RAeqppXh\n1S3HEJQVXOj24JlXG/B/a2ahNPT4AcA7O07ijY+aUFpgw5UTi/DRvrN4au1e/HjJ9ejq8+Kzw22Y\nMjYft1xVHvNz5t80AdsaWvHOjlOomlY2YMD1OH3444ZD+Kjh7ICzJUaDhJuvHIV775iKkgIbAkEZ\nL79/BFvrW/Gb1z7HT5beENM1nsqBk514/s39cHr8mH/TBCyZOxWG0IuU264ZA6vZhJWvf47/WrcP\nv3jwJhTkXnylKivqgT7e/iC8/QF4fEF4+gM41+HGpwfPa9Wb3WrErVeXIz/HAqNRwpHmbhwK/ffe\nrtO4784rcM3kzFfo6QSCMk609uJshwtnO9zo7POhrMCGcWW5GFfqwLiy3KTBJCsKGo61471dp3E4\n6kX2l64u15Ydup392HukDWs2HcVbHzdhbtV4fLVqXMyugXjtPR68+VETPmk8BwXAtAmFWHzHFK1x\nM55BkpCfY0F+jiXh32D43PdzHS6c7XTjQpdH/a/bg5Y214AzlABgsxjhsJvhsJthCjVKPvPjO9I8\nmmnCeNOmTfD7/Vi7di0aGhpQW1uL5557DgDg9/tRW1uLdevWwWazYenSpZg7dy5KSrL/jyTb/uX5\nHUnvmzw2/7LCOJkjUdMwA/H4AvhgjxoeEmK3U82aVgajUcKugxcwqsiOBxfMQFmhDT/578Tf4/+9\ntEd7/67bJsJsNGB/U6d22ljl2HyMLcnVwnjKuHxMm1CIiaPz8Pyb+7UGs4nleQl7qQFgXJljwKMq\nJ4yK/KFEV8bJTBjl0N6fPDY/5oXIxNAf3c7959DV54Un1I09dVwBZoSmv9/b3Yy2bi9unF6GhbdW\n4NR5tXK+c9YEbdzN553a1wHA9Ar1aw83d2trzgDQ0ubEs6834nyoQc9hN2PZQzei9uU67TKa5UV2\njC/Lxf/45lX4pPEcxpfmatPE4d8lfHGRqeML4fIEtDAO/3xJklDosKC9x5s2jMMvvArjwri0wI5x\npbk4eKoLr2w+iv6AjIcWzkC/P4g/bzqKJ/+8F//nOzNRVmjHhk+b8dq2EyjOt+InS29AaaE6rffq\nh8fxr3/YDX9QhiQB1XdekTCtZ7ea8NDCGXjmLw343foDeKT6BuSHAi4QlPHBnjNY/8lJuLwBjCnJ\nwfybJuDqScUIygrae7xwevy4urI4Zo3YZDTgu1+fjn6/jB37z+F37xzED+66StvOJisKjp7uxpk2\nFzr71Os+y4qCQ6e6cabNCaNBwoMLpmNOqDM+2qzpZVj0lcl4fdsJPL22Hv+46Jqka8iyoqClzYXD\nzV04da4Pp9ucaOv2wOsLJr0EiwRgRkUhvnTNaNw8ozyhCuzs9eKdnafw4d4W/OovDbi6shgLbqnA\nVROLtMfW5w/i6JluHDzZpZ7E5gvA61OXYvJzzcjLsYSOSTXAZFTfFhXmINgfQEmBDWWFduTaTGmn\nYNt7PNha34qPGlpjdn3Es1mMmFFRhKsrizGuNBeKoiCoKDjX4camPWe0pa5rKosx/2b1/2/8z+5z\n92NzXQs+2HMGb39yEu/uPIVZ08swt2o8rhhfAEmSEAiqx7pu3NWMzXVnEAgqGF+Wi8V3TMG1k0su\neUpZkiQU5VlRlGfFlZOKY+6TFQXdfT4tnC90edDeoy7XOd1+9Hn8aGl3IRCUIWFwPz9lGNfV1WH2\n7NkAgJkzZ6KxsVG77/jx46ioqEBenvrENmvWLOzevRsLFiy4qF/4i+B7C2egNPQPPdURmWHfvG0S\n3h6g0r1c4YsC/OS+G/DSxsMYW+bAh3tbcMtV5dgfOqjD4wto67HRxpXmxlwJ6ltfnoS7Z0/G4ah9\nz2WFNuTnWFBeZMfBU10wGQ3aes9Vk4ohScCBUIBcf0WpFmq3XzsGO/afw8yppSjKs8asZwJqB/T4\nskgADyaMx8eFcbQJoTDec6QNe6IODLluSgnKCmwYW5qL1tDvetOV5cjLseDRpTdon1c5Jg8GCZAV\ndWo4HCKTRufBYjZgy94W7G/qRF6uGbOmjcKGXc3odfVjyrh8+AMyvv93VyIvx4I7bhiHNZuOAgBu\nmFYGSZIwZWwBpoyNfRU/tiQXkqQuExQ4LCgrsOG6KSVY88FRWC1GjIt6bAq0ME49TR0O44HWha+d\nUoINnzZj18ELGF+Wi9uvHQODQYKnP4jXt53A//3dp8ixmuDyBlCUFwliAFh460SUFNiw6t1DKC/K\nQc38aQm/j/ZzJpfgjhvG4cO9Lfj5iztx44xR8AdkfH6iA31u9d/qfXdegTtuGBdTYaVqpJIkCQ8t\nnIH2Hg8+O3QBZ9tdmDtrPNq6PPjs8AW0D7CcYzJKuG5KiTZ7ksxdX5qIHqcPm+tasPwPu/H3X67E\n5LH5yMsxo6PXi9Z2N46c7sbh5q6YpR9LqJ8gx2qCzWqC3WqC3WKEzWKCzWpEfq4FM6eUplyjL863\noWb+dMyZORavbD6G/U2d2N/UiTElObBZjHD7gujo8SZ0fof/3QyWw27GpNF5mDQmT1urzrOb0efx\n43Cz+rsda+mBoqhNgl+rGo+K0Q6MKclFcZ4Vbd0etLa7cPqCEwebu1EfmhUb6DGffd0YzLtpAsaX\nOQYYiSovx4K/v70SC26pwI7Gc/ig7gx2HbyAXQcvwGYxwh+QY2bSSvJtWPSVStx61eiMNmAZJAnF\n+TYU59u0F/CXK+VfrNPphMMReaCMRiNkWYbBYIDT6dSCGAByc3PR1zf0Fd9IVDkmP+Yas3m5Fi3g\nPL4A3v9M7Yj9ec0sPLF6T8LXL/rKZDQ2dWhTt2GDDelbrh4NKdRBGHupRnVNd3yZA489MAuKomDe\njeMxujhHW8cdWxJ5Yl/xP29V31EAm9WEH//mY3V8syux8NaJAIArJhRqv+/XqtR126/fUoH9Jzux\n9M5p2vdy2M24+cpyfHrgPKxmI+ZWjcP6T04iKCv4xpcmYsncqVrjz3VTSnDrVeVYcEsFcqwm5NhM\nMdON0d2/yeTazLhqUhG6+ny4/doxMfcV5FoSKvPKMXnak/z8mybgD387hNJQ6MUzm4xYcud0rH3/\nMG6cEdlOZTIa8PWbKvD2JyfVV8vdHu2iHPfPm4avhda1tf8X147B0dPdKHBY8c3bJiX9XSxmIxbf\nMQWfH+/ATVeWQ5IklBfn4KYZo1Ccb41pcBlTnIvT553aC4Rk7rljCtZsOqrtY442MxTGZYU2/Gjx\nTO1J7Zu3TUJZoQ0f7m1FW7cHt1xVjr+7dWLC+vjNV5bj2sklsJqNaZ8QH5g3DWNLcrBu2wlsrVfX\n9/NzLZh/0wQ8+M1r4HNf/H5ks8mA/33vTPx163FsqWvR+hIsZgNuv24Mrp5UjJJ8G0wmCbKsvtiz\npzmxDFCD/oH503HF+EL8ccOhpLsXSvJtuH5qKaZXFGHKuHyUF+UMWTBUlOfh0aU3oOlsLzbuasae\nw20wGCTYrSaMK83FlZOKcPWkYq2RzmIyaNsXnW4//AEZ/qCMQOitPceCsxecaO/2oL3Hi5Z2Jxqb\nOtHYNPApepKkNk3OuX4sbpoxKmEZoDjfps0QAeo6+P6Tnejs9cFokCAZJNgtRtx0ZflFTfVbzUbc\nccM4zLl+LI6c7sbmuha0drhgMxthMRvVKnxiEe64flzMZWFHEklJccJBbW0tZs6ciYULFwIA5syZ\ng61btwIADh8+jKeffhovvPACAGDFihWYNWsW5s+fPwzDJiIi0o+ULyGqqqqwbds2AEB9fT2mT5+u\n3Td58mScOnUKPT096O/vx+7du3H99ddndrREREQ6lLIyVhQFy5cvx+HD6jTPihUrsH//frjdbixZ\nsgRbtmzBypUrIcsyFi9ejPvuu2/YBk5ERKQXKcOYiIiIMm9krnQTERHpCMOYiIgoyxjGREREWcYw\nJiIiyrJhCWNZlrFs2TJUV1ejpqYGzc3Nw/Fjv5AaGhpQU1OT7WHolt/vx6OPPor7778f9957LzZv\n3pztIelOMBjET3/6UyxduhT33Xcfjh49mu0h6VZHRwfmzJmDpqambA9FlxYtWoSamhrU1NTgZz/7\nWcrPHZarNqU645qGzosvvoi33noLubnpj4ukS/P222+juLgYTz31FHp6enD33Xdj7ty52R6WrmzZ\nsgUGgwFr1qzBrl278Mwzz/D5IgP8fj+WLVsGu92e/pPpovl86ulxq1evHtTnD0tlnOqMaxo6EydO\nxLPPPgvuVsucBQsW4Ic//CEAdcbHaLy06+ZScnfeeSd++ctfAgBaWlpQUJD8vGi6dP/+7/+OpUuX\noqysLNtD0aVDhw7B4/Hg4YcfxoMPPoiGhoaUnz8sYZzsjGsaWvPnz2c4ZFhOTg5yc3PhdDrxox/9\nCD/+8Y+zPSRdMhqNeOyxx/Bv//ZvuOuuu7I9HN157bXXUFxcjNtvvx0A+AI+A+x2Ox5++GH8/ve/\nx+OPP45HHnkkZe4NSxg7HA64XJEr/oQvNkE0Ep09exYPPvgg7r77bnzjG9/I9nB0q7a2Fhs3bsQv\nfvELeL2JV1uiS/faa6/hk08+QU1NDQ4dOoTHHnsM7e2JV1eiSzdp0iR861vf0t4vLCxEW9vA13MH\nhimMU51xTTSStLe34/vf/z4effRRfPvb3872cHTpjTfewG9/+1sAgM1mgyRJfPE+xF566SWsXr0a\nq1evxowZM/Dkk0+itLQ028PSlddeew21tbUAgPPnz8PpdKZcEhiWBq558+Zh+/btqK6uBqCecU2Z\nc6kX06b0nn/+efT19WHlypVYuXIlAOB3v/sdrNbk16Kli7NgwQI89thjeOCBBxAIBPDzn/8cFsvg\nL7dHJILFixfjpz/9Ke6//34Aau6lelHJs6mJiIiyjHM/REREWcYwJiIiyjKGMRERUZYxjImIiLKM\nYUxERJRlDGMiIqIsYxgTERFl2f8HdSDPoyAUx/gAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from gatspy.periodic import LombScargleFast\n", "\n", "t, y, dy = create_data(100)\n", "model = LombScargleFast().fit(t, y, dy)\n", "period, power = model.periodogram_auto()\n", "plt.plot(period, power)\n", "plt.xlim(0, 5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, to illustrate the different computational scalings, we'll evaluate the computational time for a number of inputs, using ``LombScargleAstroML`` (a fast implementation of the $O[N^2]$ algorithm) and ``LombScargleFast``, which is the fast FFT-based implementation:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from time import time\n", "from gatspy.periodic import LombScargleAstroML, LombScargleFast\n", " \n", "\n", "def get_time(N, Model):\n", " t, y, dy = create_data(N)\n", " \n", " model = Model().fit(t, y, dy)\n", " t0 = time()\n", " model.periodogram_auto()\n", " t1 = time()\n", " result = t1 - t0\n", " \n", " # for fast operations, we should do several and take the median\n", " if result < 0.1:\n", " N = min(50, 0.5 / result)\n", " times = []\n", " for i in range(5):\n", " t0 = time()\n", " model.periodogram_auto()\n", " t1 = time()\n", " times.append(t1 - t0)\n", " result = np.median(times)\n", " return result\n", "\n", "N_obs = list(map(int, 10 ** np.linspace(1, 4, 5)))\n", "times1 = [get_time(N, LombScargleAstroML) for N in N_obs]\n", "times2 = [get_time(N, LombScargleFast) for N in N_obs]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFsCAYAAAAtwdttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VPW9//HXzGSy7ytkTyBkARIIIgoEBAQJRK+tWvEq\nIoi2vdra/lzaW6tSdyu29Wq59rIooJW6K5tQIApENoEkhCSE7Hsm+zJJZj2/P4IpVBaBJJNMPs/H\nw0eczHI+CSfznu/3nPP5qhRFURBCCCHEsKC2dQFCCCGEGDgS/EIIIcQwIsEvhBBCDCMS/EIIIcQw\nIsEvhBBCDCMS/EIIIcQw4mDrAs7nwIEDbNu2ja6uLpYvX05cXJytSxJCCCHswqAM/u7ubp577jny\n8vLIyMiQ4BdCCCH6yKCc6p81axadnZ1s2LCBH/3oR7YuRwghhLAbAx78WVlZLF68GACr1crTTz/N\nokWLWLx4MeXl5QA0NTXx3HPP8cgjj+Dr6zvQJQohhBB2a0CDf/Xq1fz+97/HZDIBsGvXLkwmE5s2\nbeKxxx7j5ZdfBuCVV16hsbGR1157jR07dgxkiUIIIYRdG9Bj/BEREbz55ps88cQTABw9epSUlBQA\nkpKSyMnJAXqCXwghhBB9b0BH/PPmzUOj0fTe1uv1uLu7997WaDRYrdaBLEkIIYQYVmx6cp+7uzt6\nvb73ttVqRa2+vJJkcUEhhBDih7Pp5XzJycmkp6eTmppKZmYmsbGxl/0aKpWK+vr2fqhODFcBAR6y\nT4k+JfuU6GsBAR5X/FybBL9KpQJg7ty5ZGRksGjRIgBeeuklW5QjhBBCDBsqxQ7myuWTtOhLMjoT\nfU32KdHXrmbEPygb+AghhBCif0jwCyGEEMOIBL8QQggxjEjwCyGEEMOIBL8QQggxjEjw94Njx75l\n/vwb0Onqer/3v//7Btu3b7ngc7Zv38L+/XuveJsPP/wg5eWlV/z8tWv/xmeffXzFz7+Q4uJCsrKO\nX/QxH3/8AQCHDh3giy8+7fMahBBC/IsEfz/Rah158cU/9N7+rnfBhaSmpjF9+owr3l7P6198G5d+\nft9LT99NSUnxRR+zYcNaAKZMuZ5bbpFlmIUQoj/ZtHPfQPhgTyFH8nV9+pqT4wL5yezRF7xfpVKR\nnHwNoPDxxx9w220/Oef+t956k1On8mhtbWX06Bh+97tnWLv2b/j5+VNRUc7o0TGkpqbR2NjAE0/8\nmrVrN/LWW2+SnZ2J1Wrlzjv/k1mzbjzPlhW2bdtMRsZejEYjjY0N3HHHXezb9zXFxUU8/PAjTJ8+\nkyVL7iIsLIy6ulpGjx7DE0882fsKx48fZePGt3F0dESnq+M//uM2jh07QmHhae64YxG33no7x48f\nZfXq/0WtVhMSEsrjj/+OnTu3c+BABgaDgerqSu6+ewmTJ09h+/YtODo6EhsbR21tDZ9++hFmsxmV\nSsWLL77KZ599TFtbG6+99goJCWMpKyvlZz97mPfff5c9e3ai0TiQlDSRn//8F6xd+zdqa2tobm6i\ntraWX/7y/3Httdf10b+qEEIMDzLi7wff9UR69NHf8sEHf6eqqrL3vs5OPZ6envz5z39lzZoN5Obm\n0NBQ3zvivvnmW/nyy60A7NixjYULb+HAgQxqaqpZtWoNr7/+v2zYsI6Ojo4Lbr+rq4tXX32du+9e\nwqeffsSLL77KE0/8jq1bNwNQW1vNr3/9BKtXb6C9vZ29e7865/n19TpeeOFVHn30v1m/fi1PPfUc\nK1f+D59//gkAr7zyAi++uJI33/w/AgIC2b59CyqVCr1ezx//+GdefvlPvPvuO/j7B7Bgwc3ceefd\nxMePpbKygldf/QurVq0hMjKKQ4cOsmTJ/Xh6evLoo7/p3X5RUSHp6bt46623eeutdVRWlvPNN/tR\nqVQ4OjqycuX/8Mgjj/KPf/z96v+xhBBimLH7Ef9PZo++6Oi8P3l6evHLXz7K888/w/jxSQA4OjrR\n1NTEihVP4uLiSmdnJ2azGeiZKYiMjMJisVBbW8uePbt4/fVVfPbZx5w6lc8vfvFTgDP31zB6dMz3\ntqlSqYiJ6VnzwM3NncjIKAA8PDwwGo0AREZG4+fnD0BiYhIVFWXnvEZ09Cg0Gg3u7u6EhITi4OCA\nu3vP85ubm2lqauSpp3qC2mAwMHnyFEJDw4iJGQNAQEBg77YURen9IOTt7cPzz6/AxcWF8vIyxo1L\nPO/vrby8lLFjx/eu5JiUNJGSkiKA3m0EBgZhNBp+8L+FEEKIHjLi72fTpqUQHh7Re2LfwYPfUF9f\nx4oVL/Dgg/+F0WjoDcbvvi5ceAurVr1OVFQ0bm7uREREkZw8iTfe+Bt//vNfmTXrRoKDQy64zUsd\nr6+sLO+dMThxIpvo6H//YHTh53t7exMYGMgrr/yJN974G/fccx/XXHPtBber0WhQFIWOjg7Wrfs/\nnn32JX7zm9/j5OTU+5jvmkZ/9/NHRESSm5uDxWJBURQyM48TFhZxydqEEMLeNbV1s3Hnqat6DQn+\nfqBSqc4JwUceebQ36BISxlJdXcUvf/kzXn99JQkJ42hoqO99HsCsWTdy+PBB0tJuBWD69Bm4uLjy\n0EMP8OCDS1CpVLi6ul50+2d//df3e746Ojrx/PNP8+CD9xEUNIKpU6df8Hn//v8qlYpHHnmUxx57\nhJ//fBlffPEJkZHR59lez//Hxsbx8ccfUFCQz/jxSfz0p0v53e8eIywsvPfnjoyM4rnnnup9/ejo\n0cyefSM///n9PPjgEkaODGbGjBsuWpsQQti7g7m1PL32MOnHqq7qdWSRnmHo3nvvZMOGf9i6jEFL\nFlQRfU32KXE1OrpMvLvzFIfzdDhq1SyaHcMd8+Ku+PXs/hi/+D4ZKQshxNBwsqSJtVtzaekwMirE\nk+VpCQT5XHjG94eQ4B+G1q/fZOsShBBCXITBZOGjr4rYfbQSjVrFj2dEk3pdOBr11R+hl+AXQggh\nBpGSmjZWb86ltqmTkX6uPHjzWCJGePTZ60vwCyGEEIOAxWpl6zdlbP6mFItVYe41Ydw2MxpHraZP\ntyPBL4QQQthYbVMnqzfnUlLTho+HE/cvjCch0rdftiXBL4QQQtiIoih8dbyKf+wpxGi2cv3YIO6e\nOwZXZ22/bVOCvx/U1FSzZMldxMb+63KLSZMmc999y7nhhut6u/hBTxe9cePGs2XL5xiNRkpLixkz\npud5zzzzPP7+Ab2vuWLFk/ztb29fcV23334z77//CVpt/+1Q36mrq6Ww8DTTpqVc8DGff/4JCxfe\nQklJMRkZe7nvvuX9XpcQQgwWze0G3t6eR05xE27ODixbGM+18UH9vl0J/n4SFRXNG2/87Xvf9/Ly\nOu/3b7ppAbW1NTzzzO/Oe39fGMjL+I4ePUJ5edlFg//dd98hNTWNmJgxva14hRBiODiSr2PDl/no\nu82Mi/Jl6YJ4fDycLv3EPmD3wf9J4RaO60706WtODBzPj0en9elrwr9a1l7Kww8/SExMLMXFRbi6\nupCYOJHDhw/Q0dHOn/70V/bt+4pDh76hpaWV1tYWli17sLfzHfSMxl999UUMBgNOTk488cSTWCwW\nnn76vwkKGkFtbQ1z5syjpKSIgoJTXH/9NH7604coKirk9ddXoigKXl5e/Pd/P82pU/m8994GHB21\nVFdXMWfOPO655z7effcdDAYD48Yl4ubmxjvvrMFqtdLV1cUzzzxPVtYxGhsbWbHiSe64YxGfffYx\nf/jDi+zcuZ0PP3wfrdaR0NAwnnjiyfOu/Jea2ve/fyGE6G+d3Sbe/WcBB0/W4eig5p55Y5g1MWRA\nB2Z2H/y2Ulpa3LuoDsAzz7yAv78/bW1t53z/4Yd/fc4hgR9CpVKRkDCWRx55lEcf/SUuLs78+c9/\n5YUXVpCZeRSVSoXVqvD666tobGzgpz9d2jvyVhSFv/71dW6/fRHXXTeVb789zFtvvcmDD/4XNTXV\nvP76Krq7u7njjlv47LMvcXJy4vbbb+anP32IV155niefXEFERCRbtnzOe+9tYPLkKdTV1bJhwyaM\nRiO33jqfe+9dxuLFSykvL2P69Bl8+ulHPPXUc/j7+7Nx49ukp+/i3nuXsX79Ov7whxc5cSILgLa2\nVtat+z/efvvvuLi48MYbf+Lzzz/B1dUVvV7Pn/70BpWVFfzmN7+W4BdCDDl5pU2s2ZpHc7uBqJGe\nLE+LZ6Sf24DXYffB/+PRaf0yOr+UyMjzT/V7enr2yVT+dx8W3N3de3vln70C36RJkwHw8/PH3d2D\n1taW3ucWFxeycePbvPfeehRF6T3mHxwcgqurGxqNA76+fnh49Fw3+t0H0bKyElaufAkAs9lMWFg4\nAKNGjUKtVuPs7Ny7JsHZq/L5+/vzl7+8iqurK/X1OhITJ5z3Z6quriIqKhoXFxcAkpKSOXz4IGPH\njjvvyn9CCDEUGE0WPv66mH9+W4FapeLW6VEsnBrRJ814roTdB7/9uvi0UH5+LnAbTU2NdHd34+3t\n03tfREQkd921mHHjEikuLiQ3N6fnFS8x1RQeHslTTz1LYGAQmZnHaG1tvWAtarUaq9UKwB//+CIf\nfPA5Li4uvPDCit7v98xMWHqfM3JkMCUlJXR3d+Ps7Mzx40cJD4/4QbUJIcRgVFbbzuotuVQ36Bnh\n68oDNycQNdLTpjVJ8PeTCwfVxQPsYgF3OeFXWVnBI4/8F52dHTz22G9Rq9VAz+p3Dz30K1aufBmj\n0YDBYOBXv3r8PK///f9/7LH/5rnnnsZisaBWq/ntb5+ivl533udFR49iw4Z1xMbGMW9eKg89tBx/\n/wDCwyNpbGwAIClpIo899ghLlz6ASqXCy8ub++9/kF/84qeo1WpCQ8P4+c9/we7dOy9SmxBCDD4W\nq5XtB8v5fH8JFqvCnORQbp81Cqc+bsZzJWR1Pju0ffsWWlpauOuue2xdypAkK6mJvib71PBS19zJ\nmi25FFW14e3uyLKF8YyL8uvTbQQEXHkLXxnx2ymZGRdCiIGlKApfZ1Xzj92FGEwWro0P5J55sbi7\n9H/vlMshwW+H5Ix3IYQYWK0dBt7enk92USOuTg48eEsC1yWMsHVZ5yXBL4QQQlyFo6d0rP/yFB1d\nJhIifVi2IB5fT2dbl3VBEvxCCCHEFejsNvP+rgIycmrROqj5zxtjmD0pFPUgP9YqwS+EEEJcplPl\nzazZkkdjWzcRIzx4IC2BYP+Bb8ZzJST4hRBCiB/IZLbw6d4SdhwuR6VScfPUSG6eFomDxjbNeK6E\nBL8QQgjxA5TX9TTjqarXE+jjwgNpCYwK8bJ1WZdNgl8IIYS4CKtV4cvD5Xy6txiLVWHWxBB+Mms0\nTo62b8ZzJST4hRBCiAuob+lizZZcTle24uXmyNIF8SSO6ttmPANNgl8IIYT4N4qisD+7hr/vPo3B\naGFSbAD33hSLh6ujrUu7ahL8QgghxFna9Ebe2Z5PZmEDLk4alqfFc/3YEXazWJgEvxBCCHHG8dP1\nvLM9n/ZOE3Hh3ty/MAE/r8HbjOdKSPALIYQY9roMZjbtPs2+7BocNGoWzR7NjZPDBn0znisxqIP/\nwIEDbN26leeff97WpQghhLBTBRUtrNmSS0NrN+GB7jxwcwIhAe62LqvfDNrgLy8vJz8/H4PBYOtS\nhBBC2CGzxcpn+0rYfrAMVLDw+gj+Y3rUkGrGcyUG7U8XHh7O0qVLbV2GEEIIO1RZ38Fz679l28Ey\n/L2d+e3dydw2c5Tdhz7YKPizsrJYvHgxAFarlaeffppFixaxePFiysvLbVGSEEKIYcCqKHx5qJxn\n3zlCha6DGUnB/GHZtcSEetu6tAEz4FP9q1ev5osvvsDNrWcxg127dmEymdi0aRNZWVm8/PLLrFq1\naqDLEkIIYecaWrtYuyWPUxUteLpquW9BPBNG+9u6rAE34CP+iIgI3nzzTRRFAeDo0aOkpKQAkJSU\nRE5OzjmPf/XVVwe6RCGEEHZEURQyTtTwzLrDnKpoYWKMP88unzIsQx9sMOKfN28elZWVvbf1ej3u\n7v86e1Kj0WC1WlGr7f84ixBCiP7V3mlkw5enOFpQj7OjhmUL4pk23n6a8VwJm5/V7+7ujl6v7719\nJaEfEODR12WJYU72KdHXZJ8aeEdya/mfDzJpaTcwNtqPX9+VTJCvq63LsjmbB39ycjLp6emkpqaS\nmZlJbGzsZb9GfX17P1QmhquAAA/Zp0Sfkn1qYHUbzfxjTyFfZ1bjoFHxk1mjmTc5DLXFYjf/Dlfz\nQdJmwf/dNMvcuXPJyMhg0aJFALz00ku2KkkIIcQQV1jVyprNuehauggN6GnGExZov814roRK+e4s\nuyHMXj7BicFBRmeir8k+1f/MFitfZJSw9UAZKDB/Sji3pkSjdbDP88WG5IhfCCGE6AtVDXrWbM6l\nrK4dfy9n7l8YT2y4j63LGrQk+IUQQgxJVkVh17eVfPRVEWaLlenjR3LXjTG4OEm0XYz8doQQQgw5\nTW3drN2aR15ZMx6uWu6bP5aJYwJsXdaQIMEvhBBiyFAUhYO5dby7s4Aug5kJo/1ZkhqHl5ujrUsb\nMiT4hRBCDAkdXSY27jjFkXwdTloN96XGkZI4clg347kSEvxCCCEGvZziRtZuy6O1w8joEC+Wp8UT\n6CPNeK6EBL8QQohBy2C08MFXhaQfq0KjVnHbzGhSp0SgVsso/0pJ8AshhBiUiqp7mvHUNXcR4u/G\nAzcnEB4krY+vlgS/EEKIQcVssbLlm1K2fFOGoijMmxzGbTOj0TpobF2aXZDgF0IIMWjUNOpZvTmX\n0tp2fD2duH9hAvER0oynL0nwCyGEsDlFUdhzrIoP0wsxmq1MHTeC/7xxDK7OElN9TX6jQgghbKq5\n3cC6bXmcLGnC3UXL8rQErokLtHVZdkuCXwghhM0czqtj445T6LvNJI7y477UOLzdnWxdll2T4BdC\nCDHg9N0m3t1ZwKHcOhy1au69KZaZE4KlGc8AkOAXQggxoE6WNrFuax7N7QZGBXuyPC2BIF9pxjNQ\nJPiFEEIMCKPJwkdfFbHraCUatYofpUSx4PoINGq1rUsbViT4hRBC9LvS2jZWb86lprGTkX6uPHBz\nApEjPG1d1rAkwS+EEKLfWKxWth4oY3NGKRarwo3XhHL7zFE4aqUZj61I8AshhOgXdU2drN6SS3F1\nGz4eTixbGM/YSF9blzXsSfALIYToU4qi8FVmNf/Ycxqjycp1CUHcPW8Mbs5aW5cmkOAXQgjRh1o6\nDLy9LZ8TxY24OTuwbEE818YH2boscRYJfiGEEH3i23wdG3acoqPLxNgoX5YtiMfHQ5rxDDYS/EII\nIa5KZ7eZ9/5ZwIGTtTg6qLl77hhmJ4dIM55BSoJfCCHEFcsra2bt1lya2gxEjfRgeVoCI/3cbF2W\nuAgJfiGEEJfNZLbw8dfF7DxSgVql4j+mR7Hw+ggcNNKMZ7CT4BdCCHFZymrbWbMll6oGPUG+rjyQ\nlkB0sDTjGSok+IUQQvwgVqvC9kNlfLavBItVYXZyCHfMGo2TNOMZUiT4hRBCXJKuuZM1W/IorGrF\ny92R+xfEMy7az9ZliSsgwS+EEOKCFEVhb1Y1m3YXYjBZmBwXyOKbYnF3kWY8Q5UEvxBCiPNq1Rt5\nZ1seWUWNuDg58ODNCUxJCJLL9IY4CX4hhBDfc6ygnne259PRZSI+wof7F8bj6+ls67JEH5DgF0II\n0avLYOb9XafZf6IGrYOau+bEMOeaUNQyyrcbEvxCCCEAOFXezNqteTS0dhMR5MEDNycQ7C/NeOyN\nBL8QQgxzJrOVT/cVs+NQOaggbWokt0yLlGY8dkqCXwghhrEKXQerN5+ksl5PoI8Ly9MSGB3iZeuy\nRD+S4BdCiGHIalXYcaScT/cWY7Yo3DAhmJ/MHo2zo8SCvZN/YSGEGGYaWrpYsyWXgspWvNwcWbog\njsRR/rYuSwwQCX4hhBgmFEVhf3YNf999GoPRwqTYAO69KRYPV0dblyYGkAS/EEIMA616I+u355NZ\n2ICLk4b7F8YzddwIacYzDEnwCyGEnTtWUM/6L/Np7+xpxrNsQTx+XtKMZ7iS4BdCCDvV2W3m/V0F\nZOTUSjMe0WtQBv+xY8f44IMPAHjyySfx8PCwcUVCCDG05JU1s25rLo1tBiJGePBAmjTjET0GZfB/\n+OGHPPvss2RnZ7Nt2zbuvPNOW5ckhBBDgtFk4ZO9xew8UoFapeKWaZGkTZVmPOJfBmXwWywWHB0d\nCQgI4ODBg7YuRwghhoSy2nZWb8mlukFPkK8rD6QlEB3saeuyxCAz4MGflZXFypUr2bhxI1arlRUr\nVlBQUIBWq+WFF14gPDwcZ2dnjEYjOp0Of3+5tlQIIS7GYrWy7UAZX2SUYrEqzJkUyu03jMJJq7F1\naWIQGtDgX716NV988QVubj3HmXbt2oXJZGLTpk1kZWXx8ssvs2rVKu68806eeeYZzGYzzz777ECW\nKIQQQ0ptUydrtuRSXN2Gj4cTyxbEMzbK19ZliUFsQIM/IiKCN998kyeeeAKAo0ePkpKSAkBSUhI5\nOTkAjB07lpdeemkgSxNCiCFFURT2HKviw/RCjGYr140N4u65Y3Bz1tq6NDHIDWjwz5s3j8rKyt7b\ner0ed3f33tsajQar1YpafXknoQQEyFn/om/JPiX6Wl/uU42tXfxl03EyC+rxcNXy6/9MZnpSSJ+9\nvrBvNj25z93dHb1e33v7SkIfoL6+vS/LEsNcQICH7FOiT/XlPnUwt5Z3dxTQaTAzPtqPpQvi8HZ3\nkn12mLmaD5I2Df7k5GTS09NJTU0lMzOT2NhYW5YjhBCDVkeXiXd3nuJwng4nrYZ758cyMylYWu6K\ny2aT4P9uR507dy4ZGRksWrQIQI7rCyHEeZwobmTdtjxaO4yMDvFieVo8gT6uti5LDFEqRVEUWxdx\ntWSKS/QlmeoXfe1K96luo5kP0ov46ngVGrWKW1OiSJ0SgVoto/zhbshO9QshhDi/wqpW1mzORdfS\nRUiAGw+kJRAeJCediqsnwS+EEIOI2WLl8/0lbDtYBgrMnxLOj1Ki0TpIy13RNyT4hRBikKis72DN\n5lzKdR34ezmzPC2BMWHeti5L2BkJfiGEsDGrVWHnkQo+2VuE2aIwI2kkd86OwcVJ3qJF35O9Sggh\nbKi+pYu1W/MoqGjB01XLfanxTIiRNUpE/5HgF0IIG1AUhf3ZNfx992kMRguTxgSweH4snq6Oti5N\n2DkJfiGEGGCteiPrt+eTWdiAi5OG+xfGM3XcCGnGIwaEBL8QQgygYwX1rP8yn/ZOE/ERPixbEI+f\nl7OtyxLDiAS/EEIMgM5uM+/vKiAjpxatg5q75sQw55pQ1DLKFwNMgl8IIfpZdmE9f3rvKI1tBiJG\nePBAWgLB/m62LksMUxL8QgjRT4wmC5/sLWbnkQrUKhW3TIskbWokDhppxiNsR4JfCCH6QVltO6u3\n5FLdoCckwI2lqfFEB3vauiwhJPiFEKIvWaxWth0o44uMUixWhTmTQvnZ7Um0t3bZujQhAAl+IYTo\nM7VNnazenEtJTRs+Hk4sWxDP2ChfnB0dkPUexWAhwS+EEFdJURT2HKviw/RCjGYr140N4u65Y3Bz\n1tq6NCG+R4JfCCGuQnO7gXVbczlZ2oybswP3pyUwOS7Q1mUJcUES/EIIcYUO5tby7o4COg1mxkf7\nsXRBHN7uTrYuS4iLkuAXQojL1NFl4t2dpzicp8NJq+He+bHMTAqWlrtiSJDgF0KIy5Bd1Mjb2/No\n7TAyOsSL5WnxBPq42rosIX4wCX4hhPgBuo1mPkgv4qvjVWjUKm6bGU3qlAjUahnli6FFgl8IIS6h\nsKqVNZtz0bV0ERLgxgNpCYQHedi6LCGuiAS/EEJcgNli5fP9JWw7WAYKzJ8Szo9SotE6SMtdMXRJ\n8AshxHlU1newZnMu5boO/L2cWZ6WwJgwb1uXJcRVk+AXQoizWK0KO49U8MneIswWhRlJI7lzdgwu\nTvJ2KWxLURTK2ivIbzrN4oBbr/h1ZE8WQogz6lu6WLs1j4KKFjxdtdyXGs+EGH9blyWGMatipbSt\nguO6bI7rTtBsaAFg8WQJfiGEuGKKorA/u4a/7z6NwWhh0pgAFs+PxdPV0daliWHIqlgpaS3neH1P\n2LcYWgFw1jhz7YhkkgMTr+r1JfiFEMNaq97I+u35ZBY24OKk4f6F8UwdN0Ka8YgBZVWsFLeWcUyX\nTabuBK3GNgBcHFy4bsQ1TAwcT6xvDFr11ce2BL8QYtg6eqqeDTvyae80ER/hw7IF8fh5Odu6LDFM\nWBUrhS0lHNedILP+BG3GnjUc3RxcuX7kZCYGJhLrMwqHPgj7s0nwCyGGnc5uM+/vKiAjpxatg5q7\n5sQw55pQ1DLKF/3MYrX0hH19T9i3GzsAcNO6MnXktSQHJjLGZxQatabfapDgF0IMK3llzazbmktj\nm4GIER48kJZAsL+brcsSdsxitXC6pZjjumwy63PoMOkBcNe6MT14ChMDE4nxju7XsD+bBL8QYlgw\nmix8sreYnUcqUKtU3DItkrSpkThopBmP6HsWq4WC5iKO6bLJashBb+oEwMPRnZSQ65kYMJ7R3lED\nFvZnk+AXQti90to21mzJo7pBT5CvKw+kJRAd7GnrsoSdMVvNnGou4rgum+z6k+jNPWHv6ejBjJCp\nJAeOZ5R3FGqVbT9sSvALIeyWxWpl64EyNmeUYrEqzJkUyu03jMJJO/CjLGGfzFYz+U2nOa47QVbD\nSbrMXQB4OXoyM3QayYGJRHtF2Dzsz3bR4G9qauK9995jz549lJaWolariYiIYM6cOdx11134+voO\nVJ1CCHFZaps6Wb05l5KaNnw8nFi2IJ6xUfKeJa6eyWIiv/k0x3TZnGjIpcvcDYC3kxfXjZzExIBE\norzCB1XYn+2Cwf/ee++xc+dO5s2bx0svvURISAgODg5UVlZy6NAhHn74YebPn8+99947kPUKIcRF\nKYrCnmNVfJheiNFs5bqxQdw9dwxuzlpblyaGMJPFRG5TAcfPhH23xQCAj5N376V3kZ5hgzbsz3bB\n4A8KCmKyDRuWAAAgAElEQVT9+vXf+35MTAwxMTHcc8897Nixo1+LE0KIy9HcbmDd1lxOljbj5uzA\n/WkJTI4LtHVZYogyWkzkNuZzvP4EJxpyMViMAPg5+zAtZArJgYlEeIQNuWZPKkVRlIs9wGw28/XX\nXzNnzhyamprYs2cPt91226D6Qevr221dgrAjAQEesk8NMYqicCivjnd3FNBpMDM+2o+lC+Lwdney\ndWmA7FNDidFiJKcxn0zdCU405mE8E/b+zr5MDExkYuB4wj1CbZ6BAQEeV/zcS57c99RTT2GxWJgz\nZw4ABw4cIDs7m2efffaKNyqEEH2lo8vExh2nOJKvw0mr4d75scxMCrb5G7MYOgwWIzkNeRzXZXOy\nMR+j1QRAgItfb9iHuYfYzT51yeA/ceIEW7ZsAcDX15fXXnuNm2++ud8LE0KIS8kuauTt7Xm0dhgZ\nHeLF8rR4An1cbV2WGAK6zd3kNOafCftTmM6EfaCrP8kBiUwMTCTEfaTdhP3ZLhn8iqJQV1dHUFAQ\nAA0NDajVg//kBSGE/eo2mvkgvYivjlehUau4bWY0qVMiUKvt701a9J0uczcnGnLJ1J0gt+kUJqsZ\ngCDXQJIDxzMxMJFgN/tfoOmSwf+zn/2MH//4xyQnJwOQlZXFk08+2e+FHThwgK1bt/L888/3+7aE\nEENHYVUrazbnomvpIiTAjQfSEggPuvLjncK+dZm7yK7P5Xh9NnmNBZgVCwAj3YKYGHAm7N1H2LjK\ngXXJ4L/55pu59tpryczMxMHBgaeeeorAwP49S7a8vJz8/HwMBkO/bkcIMXSYLVY+31/CtoNloMD8\nKeH8KCUarYPMQIpzdZo6yW7I5bgum7ym01jOhH2w2wiSzxyzH+EWZOMqbeeSwW80Gvnkk08oKSnh\n97//PRs2bODBBx/E0dGx34oKDw9n6dKlPP744/22DSHE0FFZ38HqzblU6Drw93JmeVoCY8K8bV2W\nGET0pk6y609yrD6bU02FvWEf4j6yJ+wDxhPkJpd2wg8I/j/84Q/4+vpy8uRJNBoNZWVlPPnkk7z6\n6quXtaGsrCxWrlzJxo0bsVqtrFixgoKCArRaLS+88ALh4eH85S9/oby8nBUrVuDpKX20hRjurFaF\nnUcq+GRvEWaLwoykkdw5OwYXJ+k2LqDDqCerIYfjuhOcai7EqlgBCPMIOTONP55A1wAbVzn4XPKv\n5+TJk3z22Wfs27cPNzc3/vjHP5KWlnZZG1m9ejVffPEFbm49S1/u2rULk8nEpk2byMrK4uWXX2bV\nqlX86le/urKfQghhd+pbuli7NY+CihY8XbXclxrPhBh/W5clbKzd2EFWfU/YF7QU9YZ9uEcoyYGJ\nTAgYT4Crn42rHNwuGfxqtRqj0dh7u7m5+bLP6o+IiODNN9/kiSeeAODo0aOkpKQAkJSURE5Oznmf\nd7mzCkKIoU9RFPZn1/D33acxGC1MGhPA4vmxeLr23+FFMbi1GdvJqs/hmO4Ep5uLUOjpOxfhGdY7\nje/nIusw/FCXDP57772XpUuX0tDQwPPPP8+uXbt46KGHLmsj8+bNo7Kysve2Xq/H3d2997ZGo8Fq\ntV7xZYJX08FIiPORfco2mtu7+euHWRw6WYurswO/vmsisyYNvZao5yP71OVp6WrlUGUmByuPkVt/\nmu+azI7xi+a6sIlMCZ1IgJuM7K/EJYP/1ltvZezYsRw6dAir1cpbb71FXFzcVW3U3d0dvV7fe/tq\nQh+kZa/oW9Je1TaOnqpn/Zf5dHSZiI/wYdmCePy8nGlo6LB1aVdN9qkfpsXQSqYuh+P12RS1lPaO\n7KO9IpkYOJ6JAePxcT5zUmcn1HcO399pv7bsbW5uRqfTcc899/DWW2+xatUqfvnLXzJ69Ogr3mhy\ncjLp6emkpqaSmZlJbGzsFb+WEGJo6+w28/6uAjJyatE6qLlrTgxzrglFbQejfHFpzd0tZNbncEyX\nTUlrGQoKKlREe0X2HLMPHIe3k5ety7Qrlwz+Rx99lFmzZqFSqdixYwdLlizhmWee4b333rvsjX03\nXTd37lwyMjJYtGgRAC+99NJlv5YQYujLK2tm3dZcGtsMRIzw4IG0BIL93WxdluhnTd3NZOpOcEx3\ngpK2MgBUqBjtHcXEwESSAsZK2PejS67Od9ttt/Hxxx/z3HPPER4ezpIlS/jxj3/MJ598MlA1XpJM\noYm+JNOy/c9osvDJ3mJ2HqlArVKRNjWCtKmROGjssxmP7FPQ2NXE8foTHNedoLStHOgJ+xifUUwM\nGE9SwDi8nOQ8iB+qX6f6FUUhJyeHXbt2sXHjRvLy8rBYLFe8QSHE8FZa28bqzbnUNHYS5OvKA2kJ\nRAdL3w571NDVyHFdT9iXtVcAoFapifOJYULgeCYEjMPD0f0SryL62iWD//HHH+ePf/wjS5cuJTw8\nnEWLFvHb3/52IGoTQtgRi9XK1gNlbM4oxWJVmDMplNtvGIWTVmPr0kQf0nU29Ezj12dT0V4F9IR9\nvO8YJgaOJ8l/HO6OcjjHli441a/T6S7Zk/+HPGYgDPcpNNG3ZFq279U2dbJ6cy4lNW34eDixbEE8\nY6OGz3XX9r5P1XXWnxnZZ1PZUQ2cGdn7xjAxIJHEgATctRL2falfpvr/9Kc/ERQUxK233kpUVNQ5\n9xUVFfHRRx9RX1/PypUrr3jjQgj7pigKe45V8WF6IUazlevGBnH33DG4OWttXZq4SrV6Hcd12Ryv\nP0FVRw0AGpWGcX5xTAhMJMk/AVetq42rFOdz0ZP70tPTWbt2LaWlpQQGBqLRaKitrSU8PJz777+f\n2bNnD2StF2TPn6TFwLP30dlAaWrr5u1teZwsbcbN2YF758cxOc72M4S2YC/7VI2+jmO6bI7rsqnR\n1wHgoNIQ7zeGiQGJjPdPwFXrYuMqh4erGfFf8qx+gJaWFsrLy1Gr1YSGhuLtPbhWxbKHPygxeNjL\nm7StKIrCodw63t1ZQKfBzPhoP5YuiMPb3cnWpdnMUN2nFEWhWl/bO41f26kDwEHtQIJvLBMDxzPe\nPx4XBwn7gdavZ/UDeHt7D7qwF0IMPh1dJjbuOMWRfB1OWg33zo9lZlKwXbTcHS4URaGqo+bMpXfZ\n1HXWA6BVO5AUMI7kgPGM84/H2cHZxpWKKyVrWwoh+kR2USNvb8+jtcPI6BAvlqfFE+gjx3iHAkVR\nqOyo5pgum0zdCXRdDQBo1dre5W3H+sXj7DB8Z23siQS/EOKqdBvNfJBexFfHq9CoVdw2M5rUKRGo\n1TLKH8wURaGivarnmH39CRq6GgFwVGt7VrwLTGSsXxxOGlkV0d5cMvh/8Ytf8MYbb5zzvSVLlrB+\n/fp+K0oIMTQUVrayZksuupYuQgLceCAtgfAg6b42mLUZ29lbeYDDtcdo7G4CwEnjyKTAJJIDE0nw\ni8VRwt6uXTD4H3roIfLy8tDpdOecvW+xWBg5cuSAFCeEGJzMFiuf7y9h28EyUGD+lHB+lBKN1sE+\nW+7ag+qOWvZU7ONI7THMigVnjRPXBE0gOTCReN9YHDVyieVwccHgf/nll2ltbeX555/nqaee6l0L\n2cHBAX9//wErUAgxuJTUtLFuWx5V9Xr8vZxZnpbAmDA5+XcwUhSFvKYC9lTsI6+pAIBAF39mhaUw\nZeQkmcYfpn7Q5XyD3VC8TEYMXkP10qv+ZjJb+Gx/CV8eKkdR4IYJwdwxazQuTnKq0KUM9D5lspg4\nUnecPRX7eq+3j/GOZk74DMb6xaFWyczMUNfvl/MJIYa3wspW1m3Lo7apE38vZ5amxhEfOXxa7g4V\n7cYO9lUdYG/lAdpNHahVaiYHJTM7fDrhHqG2Lk8MEhL8QogLMpgsfPJ1Mbu+7VlZ7cZJodw2cxRO\njrKwzmBSo68jvWIfh2qPYbaacXFwYW74DdwQNk3WtRffI8EvhDiv/LJm3tmej66liyBfV5amxsmx\n/EFEURRONReyu2IvuY2nAPB38WNW2HSuG3GNXHMvLkiCXwhxji6DmY++LiL9WBUqVc8Z+7dOj8JR\nls8dFExWM9/WZbKnfC/V+loARnlFMSc8hfH+CXL8XlySBL8QoldOSSPrt+fT2GYg2N+NZQviiQ72\ntHVZAugw6tlffZCvK7+hzdiOWqXmmqAJzA5LIcIzzNbliSFEgl8IQWe3iX/sKWRfdg1qlYq0qZHc\nPDVSrssfBOr0OvZU7udQzVFMVhPOGmfmhM/ghtBp+Dr72Lo8MQRJ8AsxzGUVNrBhxyma2w2EBbqz\nbEE8ESOk+54tKYrC6ZYidpfvI6cxDwA/Zx9mhaVw/chrZIEccVUk+IUYpjq6TLy/q4ADJ+vQqFXc\nmhLFgusicNDIKN9WzFYzR+uy2FOxj8qOagCivSKYHTaDRP8ENGo5z0JcPQl+IYaho6d0bNxZQJve\nSOQID5YtjCc0wN3WZQ1belMn+6sO8nVlBq3GdlSomBiYyJywFKK8ImxdnrAzEvxCDCNteiPv/bOA\nI/k6HDRq7rhhFPOuDUOjllG+Leg660mvyOBgzRGMVhPOGidmh6VwQ+g0/FykQZLoHxL8QgwDiqJw\nOE/He/8soKPLxKgQT5YtiGekn5utSxt2FEWhsKWEPRX7ONGQi4KCj5M3aWHTmRo8GRcHF1uXKOyc\nBL8Qdq6lw8DGHac4froBRwc1i+bEcOOkUNRqla1LG1YsVgvHdNnsqdhLeXsVABGeYcwJm8GEgHFy\n/F4MGAl+IeyUoih8k1PLpt2n0XebiQ3z5r4FcQT5uNq6tGGl09TJ53kH2HpqDy2GVlSomBAwjjnh\nM4jyjEClkg9gYmBJ8Athh5rautmw4xTZRY04aTXcM28MN0wMQS0hM2AauhrZU7GfAzVHMFqMOGoc\nuSF0GjeETifA1c/W5YlhTIJfCDuiKAp7s6r5IL2QLoOFsZE+LJkfh7+3HDceCIqiUNxaxp6KvWTV\nn0RBwdvJi5+MW0iS5wRctfLvIGxPgl8IO9HQ0sXb2/PJK2vGxUnDfalxpCSOlKnkAWCxWsisP8Hu\nin2UtfWsZBjuEcKcsBlMDExkRJA39fXtNq5SiB4S/EIMcVZFIf1YFR99VYTBZCFxlB/33hSLr6d0\nd+tvXeYuMqoP81VFBs2GFlSoSPQfy+ywFEZ7R8mHLjEoSfALMYTVNXfy9rZ8CipacHN2YPFN8Vw/\ndoQETj9r7Griq8oMvqk+TLfFgKNay4yQqcwKm0aga4CtyxPioiT4hRiCrFaFf35bwad7izGarSSP\nCWDxvDF4ucsa7P2ppLWM3RX7yNSdQEHBy9GTmyJmMy1kCm5auVpCDA0S/EIMMdUNet7elkdRdRvu\nLlqWLYxnclygjPL7icVqIavhJHvK91HSVgZAqHswc8JnkByYiINa3kbF0CJ7rBBDhMVq5ctD5Xy+\nvwSzReHa+ED+c+4YPF0dbV2aXeo2d/NNzRG+qthPY3czAOP945kdlkKM9yj5oCWGLAl+IYaACl0H\n67blUVbbjpebI4tviiV5jBxL7g9N3c18VZlBRtVhui3daNVapodcx+zQ6QS5Bdq6PCGumgS/EIOY\n2WJl64EytnxTisWqMG3cCO6cE4O7i9bWpdmdsrYKdpfv5Xj9CayKFU9HD+ZGzGR68HW4O8qaBsJ+\nSPALMUiV1raxbms+lfUd+Hg4sWR+LImj/G1dll2xKlayG3LZU76XotZSAILdRjAnfAaTgiagleP3\nwg7JXi3EIGMyW/gio5TtB8uxKgozkoL5yazRuDrLn2tf6TYbOFjzLekV+2jobgIgwS+WOWEziPUZ\nLcfvhV2TdxIhBpGiqlbWbcujprETfy9nlqTGMTZS1mXvK83dLXxd+Q37qw/RZe7CQe3AtOBrmRWW\nwki3IFuXJ8SAkOAXYhAwmCx8tq+YnUcqUBSYkxzKbTdE4+wof6J9oby9kj3l+ziqy8KqWPHQurMw\nai4pIdfj4ehu6/KEGFCD7l3lwIEDbNu2ja6uLpYvX05cXJytSxKiXxVUtLBuWx665i4CfVxYmhpH\nbLiPrcsa8qyKlZyGPPZU7ON0SzEAI92CmB02g8lBE9Bq5ARJMTwNuuDv7u7mueeeIy8vj4yMDAl+\nYbe6jWY+/qqY3ccqUQHzJofxoxnROGk1ti5tSDNYjByq+Zb0iv3ouhoAiPcdw+ywFOJ9x8jxezHs\nDbrgnzVrFp2dnWzYsIHHH3/c1uUI0S9yS5t4Z3s+Da3djPRzZdmCeEaFeNm6rCGtxdDK3soD7K86\niN7ciYNKw/UjJzM7LIVg9xG2Lk+IQWNAgj8rK4uVK1eyceNGrFYrK1asoKCgAK1WywsvvEB4eDh/\n+ctfKC8v58knn2TlypU88sgj+PrKSU3CvnR2m/nwq0K+zqxGrVKx8PoIbpkWidZBRvlXqqK9mvSK\nfXxbl4lFseCudSM18kZmhF6Pp6OHrcsTYtDp9+BfvXo1X3zxBW5uPQ0wdu3ahclkYtOmTWRlZfHy\nyy+zatUqfvWrXwHwm9/8hubmZl577TVuvPFGbrrppv4uUYgBkV3UyPov82luNxAa4MayhfFEjvC0\ndVlDklWxktt4it0V+yhoLgQgyDWQOWEpTB6RjKMcvxfigvo9+CMiInjzzTd54oknADh69CgpKSkA\nJCUlkZOTc87jX3nllf4uSYgBpe82sWnXaTJyatGoVfzH9CgWXh+Bg0Zt69KGHKPFyKHaY6RX7KOu\nsx6AWJ/RzA5LIcEvFrVKfqdCXEq/B/+8efOorKzsva3X63F3/9flMxqNBqvVilp95X+wAQEynSf6\nVl/tUwdzalj1URbN7QZGhXrxyJ0TiQqWY/mXq6WrlS8Lv+afhXtpN+rRqDXMjLyOhWPmEOkTauvy\nfhB5nxKDxYCf3Ofu7o5er++9fbWhD1Bf3361ZQnRKyDA46r3qfZOI+/9s4DDeTocNCpumxnN/Cnh\naNRq2V8vQ1VHDXsq9vFt7XHMigU3B1fmR8xmRuhUvJw8wTw0/v77Yp8S4mxX80FywIM/OTmZ9PR0\nUlNTyczMJDY2dqBLEKLfKIrCkXwd7/2zgPZOE6OCPVm6IJ5gf1nk5YdSFIXcpgL2lO8lv/k0AIGu\n/swOS2HKiEk4amQZYiGuxoAF/3fXzs6dO5eMjAwWLVoEwEsvvTRQJQjRr1o7DLy7s4CjBfVoHdTc\nOXs0c68JQ62W68Z/CJPFxOG6Y+yp2E+tvg6AGO9o5oTPYKxfnBy/F6KPqBRFUWxdxNWSKTTRly53\nWlZRFA6erOPvuwrQd5sZE+rF0gXxBPm69mOV9qPd2MHeym/YW3WADpMetUrNpMAJzA6fTrjH0Dh+\nfyky1S/62pCa6hfCnjS3G9jwZT5ZRY04aTXcPXcMs5JDUEt3uEuq7qglvWI/h+uOYbaacXFwYV7E\nLGaGTsXbSU6AFKK/SPALcQUURWFfdg3/2HOaLoOF+Agf7kuNI8DbxdalDWqKopDffJo95fvIbToF\ngL+LH7PCpnPdiGtwdnCycYVC2D8JfiEuU0NrF+u353OytBlnRw1L5scyIylYesBfhMlq5tva4+yp\n2Ee1vhaAUV5RzAlPYbx/ghy/F2IASfAL8QNZFYWvj1fxwVdFGIwWxkf7sWR+LL6ezrYubdBqN3aw\nv+ogX1d9Q7uxA7VKzTVBE5gdlkKEZ5ityxNiWJLgF+IH0DV38s72fPLLW3B1cuD+hfFMHTdCRvkX\nUKvXsadiH4drj2KymnFxcObG8JncEDoNH2dvW5cnxLAmwS/ERVitCruOVvLJ10UYzVYmxviz+KZY\nvN3lWPS/UxSFguYi9lTsJacxHwA/Zx9mhaVw/chrcHaQmREhBgMJfiEuoKZRz9vb8imsasXdRcvS\nBfFcGx8oo/x/Y7aaOVqXxe6KvVR11AAQ7RXB7LAZJAWMleP3QgwyEvxC/BuLxcq2g2V8tq8Es8XK\n5LhA7p47Bk836Rh3tg6Tnv1Vh9hbmUGrsR0VKpIDE5kdlkKUV4StyxNCXIAEvxBnqazv4MX3jlFY\n0YKnmyOL541hUmygrcsaFDqMekrayihpLaektYyStjJMVjPOGidmh6VwQ+g0/Fx8bV2mEOISJPiF\nAMxnRvmbM0qxWBWuHzuCu26Mwd1leK7rblWsVHfUnhP0uq6Gcx4z0i2I60dOZmrwZFwcpH+BEEOF\nBL8Y9spq21m3LY8KXQfe7o784s6JRAUMr0V1Okx6Ss8EfHFbOWVt5Rgsxt77nTXOxPnEEOUV0fOf\nZxiuWmlJLMRQJMEvhi2T2crmb0rYdqAcq6KQkjiSO2ePJiLM1677qlsVKzX6Oopby3qn7HWd547m\ng1wDifIKJ9qzJ+hHuAXKSXpC2AkJfjEsFVe3sW5bHtUNevw8nbgvNZ6xUfZ5fLrT1ElJ25nj8q3l\nlLaV020x9N7vrHE6M5oPJ8orgkjPcNxkNC+E3ZLgF8OK0WThs/0l7DhcjqLArOQQbp85Chcn+/hT\nsCpWavW6M1P2PUFf16k75zFBrgFM8IzoDfqRbkEymhdiGLGPdzshfoDTlS2s25ZPXVMnAd7OLE2N\nJy7Cx9ZlXZWe0XzFmdF8GaVtFXRbunvvd9I4Eusz+sxx+Z6gl9G8EMObBL+wewajhY/3FrH720oA\n5l4Txo9nROPkqLFxZZendzR/1pn2tf82mg909SfJcyxRXhFEy2heCHEeEvzCruWVNfPO9jzqW7oZ\n4evKsgXxjA4dGmu9d5q6KP3u2Hxbz7H5LvO5o/kxPqOJPjOSj/QKx107vK5GEEJcPgl+YZe6DGY+\n/KqIr45XoVJB6nXh3Do9Cq3D4BzlWxUrdZ31vVP2xW3l1Ol1KCi9jwlw8WO8fwLRXhFEeUYQ7D5C\nRvNCiMsmwS/sTk5xI+u/zKexzUBIgBvLFsQTNdLT1mWdo8vcRWlrxZkT8HqOzXeZu3rvd1RrGe0d\n1TtlH+kZjoejuw0rFkLYCwl+YTc6u01s2lPI/uwaNGoVt0yLZOH1kWgdbDsqtipWdJ31FJ/V6rb2\n30bz/i5+jPOLJ/rMmfbBbiPQqAfn7IQQYmiT4Bd2IbOwgQ1f5tPSYSQ8yJ1lC+IJD/KwSS1d5m7K\n2ioobi3tvW6+8wKj+e/OtJfRvBBioEjwiyGto8vE33cVcPBkHQ4aFT+aEU3qlHAcNAMzylcUpWc0\n39sgp4wafd25o3lnX8b6xZ1pdxtOiNtIGc0LIWxGgl8MWd/m63h35ynaOk1EjfRk2YI4QgL6d+Tc\nbe6mtK2i53K6tjJKW8vRmzt779eqtYzyjiTqTKvbKK9wPB1tM/MghBDnI8Evhpw2vZF3d57i21P1\naB3U/GTWaOZODkWj7ttRvqIo6LoaekfyJW3lVHfUnjOa93P2Jd5vTM9JeJ4RhLjLaF4IMbhJ8Ish\nQ1EUDuXW8fddp+noMjE61ItlC+IZ4ds3nei6zQbK2irYV19LTnUBJW3l6E1nj+YdiPaK7Lmcziuc\nSM8IvJxkNC+EGFok+MWQ0NxuYOOOU2QWNuCoVXPXjTHMmRSKWqW6otdTFIX6rgZKWst7L6n7/mje\np3cp2mivntG8g1r+ZIQQQ5u8i4lBTVEUMk7Usmn3aToNZuLCvblvQTyB3i6X9TrdZgPl7RW9l9SV\ntpXTYdL33u+gdug9Jj8hLA4/AvFyGlzX/gshRF+Q4BeDVmNrN+t35JNT3ISzo4Z7b4plxoTgS47y\ne0bzjb3H5Utay6jW12JVrL2P8XHyZlJgUm/Yh7oH947mAwI8qK9v79efTQghbEWCXww6iqLwdWY1\nH6QX0m20MC7KlyXz4/Dzcj7v4w0WI+VnzrT/btr+30fzkZ5h55xp7+00NPr1CyFEX5PgF4OKrqWL\n9dvzyStrxsXJgWUL4pk2fgSqM6N8RVFo7G6iuLWs95K6qo6a743mkwMTzzTIiSDUIxitHJsXQghA\ngl8MElZFYc/RSj76ugijycqE0f4svikWN1cVhS0l5yxF227q6H2eg0pDhEcYUWda3UZ7RchoXggh\nLkKCX9hcbVMnb2/L43RlC64eJmbOdELtUcz/5adT+W+jeW8nLyYGJvYuRRvqESKjeSGEuAzyjils\npttk5OMjx9hXdBJcmvG4pg2zupv9bUDbd6P50DPH5Xv62vs4e9u6bCGEGNIk+EW/sSpW2oztNHY1\n09jdRENXI41dzTR0N1Lb0UiHqQ1UoAnpebyboydRXjFEeYUT7RVJmHswWo3Wtj+EEELYGQl+cVW6\nzN00djXR0N3U87Wr6UzIN9HU3YTJav7+kxSwGp1RDD4EOY9gbkIiCYHR+Dh5957EJ4QQon9I8IuL\nMlvNNHW39Ib52SHf2NV0zgI1Z3N1cGGkWxB+zr74ufhi7nTm5KluKiqtKEYXEqMDuGVaFNHB0iRH\nCCEGkgT/MKcoCm3GDhq7G78X7A1dTbQYWs9pY/sdB7UDfs4+RHiF4X8m3L/76ufsi6vWBUVRyCtr\n5ov9JRRUtgIuJI7yk8AXQggbkuAfBrrN3TR2N58J9kYaupvP+tqEyWr63nNUqPBy8iTaKxJ/l3OD\n3d/FF09HD9Sq86+GpygKJ0ub+GJ/CacrWwFIGuXHLdOjiBopgS+EELYkwW8HLFYLzYaW847YG7ub\nzulidzYXBxdGuAbg5+KHn4sP/s5+vcHu6+xz2ZfJKYpCbmkzn2eUUHgm8CeM9ufmaZES+EIIMUhI\n8A8BiqLQbur43slzjWf+v6m75fzT8SoNvi4+hHmE4O/ih5+zT89XFx/8nX1x1fbNcrb/GuGXUlj1\nr8C/ZXokkSMk8P9/e/caFOV593H8u6wcDHgCARMRSJMnmtTEhNa4EY2aoEBkRfHpSGfCTJrOmE5s\n5U2ijWMdZ0waaewhE8dmZJJ2mgPkMbUIUpsM09SEBWwwgCUtMSYcYlQQPAVQTvf9vFhYNWCKcWHZ\n3d/nlYsL82fnmvlyX3u4RETGkjEX/traWt544w1M0+Tpp58mIiLC0yONiku9XbQNvGju0lnnW9+u\nCGTgqkoAABDGSURBVHz3ENvxAJOCJvKdSXGDnmOfOj6cScETr7kd7w6mafJx/Rn2Oer57MsLANz3\nP1NZkXgrcdN0Tr2IyFg05sLf3d3Npk2bKC0tpaqqiqSkJE+P5BbO7fjz/Vvx/e9nv9jmivy1tuND\nrCFE3RQ56Dn2iJBwIkKmeOR97q7gl9bz2QkFX0TEm4y58CckJFBVVcWrr77K7373O0+PM2ymadLe\n0zFoK37g+fazXeeu+ujZAVaLlfCQycyYMP2KrfjLkb9p3Pgx89520zSp7Q/+5/3BT7gjkhWJ8cRG\nK/giIt5gVMJfU1PDjh07eO211zAMg61bt3L06FECAwN57rnniI2N5cUXX6SxsZHHHnuM2bNnk5ub\ny86dO9m8efNojDgs3X3dV4f9a4Hv7use8vsmBU0gfuIMIkIimDp+ChHjI5jaH/mR3o53B9M0+dfn\nzuDXn1TwRUS82YiHPzc3l8LCQkJDQwEoKSmhp6eH/Px8ampq2L59O7t27SI7OxuAiooKNm3aRGBg\nIJmZmSM93lX6jD7OdZ3vD/rAW94uh/2r7vYhvy/EGkzk+IjLz7H3X7E7Xx0fTpCXfuysM/ht/cH/\nCoDv3RGJXcEXEfFaIx7+uLg4du7cyYYNGwA4fPgwCxcuBGDOnDnU1tZedX+bzYbNZhuRWUzTpKOn\nc9Dnxg88337mGtvxAZYAwkOmMH3KzYOfax8fTui4m8bMdrw7mKbJkc+cwW841R/8mc5P2psRFebh\n6URE5EaMePiXLVvG8ePHXbc7OjoIC7scD6vVimEYBAS4Z7u7u6+btv4Pphm8Ld9G1zW24ycGTSBu\nwgznW93GR7heGR8REs7k4IlYA6xumW8sM02Tms/aKLwi+N/vD36Mgi8i4hNG/cV9YWFhdHRcfgX7\njUb//2r3c6r9NKfbW2nuaOXcpQtD3i9kXDDRYZFEhU0lKjSC6NCpRIVNJTp0KpGhEQSPC/rWM3g7\n0zT58N/N5L1bx7Hj57FYIHHOLWQunUm8n37wTmSknsoQ99KakrFi1MOfkJDAe++9R2pqKtXV1cyc\nOfOGft7bHxcD/dvxwZOZOeV215X6wFZ8REg4YYGhQ2/Hd8OF7i6g64bm8EamaVJ9rJXC0gYam7/C\nAsydFYU9MZ6YSOcV/unTX3l2SA+IjJzgl7+3jBytKXG3G/lDctTCPxDdpUuX4nA4XC/ce/7552/o\n5z778NMYnVYmB0/yi+14dzBNk+pPW9nnqKepuR0LcP+dUdjnxzM9Ulv6IiK+zGKa5uDPevUy+kt6\neEzTpOrTVgqvCP7cO6OwJ97K9Kmhnh5vzNDVmbib1pS4m1dc8YvnGKZJ1dFWihz1NLU4gz/vrmjS\n5scr+CIifkbh92HO4J+m0NHAF/3Bt/UH/xYFX0TELyn8Pmgg+PtKGzh+uh2LBWzfjcY+P56bIxR8\nERF/pvD7EMM0+eiT0xQ66jl+ugOLBR74rvMKX8EXERFQ+H3CQPD3Oer50hX8adgT45kWfpOnxxMR\nkTFE4fdihmlyuP8KfyD482dPI22+gi8iIkNT+L2QYZpU1rVQ5Gjgy9bLwbfPjydawRcRkW+g8HsR\nwzCp/KSFQkcDJ1o7CLBYSOy/wlfwRURkOBR+L2AYJh/WtVDoqOdkW6cz+Hf3B3+Kgi8iIsOn8I9h\nhmHyz7pmihwNruAvuPtm0ubHEaXgi4jIt6Dwj0GGYfLP/zRTVHZF8O+5mbT58URNHu/p8URExIsp\n/GOIYZgc+o/zCv/UmU6sARYW3nMzyxV8ERFxE4V/DOgzDP757xaKyi4H/8E5N7P8gXgiFXwREXEj\nhd+D+gyDQ/9upqiskWZX8G8h7YE4pir4IiIyAhR+D+gzDCo+bmZ/WQPNZy9iDbCw6N5bWG5T8EVE\nZGQp/KNoIPhFZQ209Ad/8b238MgDcUydpOCLiMjIU/hHwZDBv286y21xREwK8fR4IiLiRxT+EdRn\nGJTXOrf0W845g7/kvuk8ouCLiIiHKPwjoLfPoPzjU+wva+D0uUuMs1pYkuC8wg+fqOCLiIjnKPxu\n1NtnUF57iqKyBlrPO4P/UILzCl/BFxGRsUDhd4PePoOyWucVvoIvIiJjmcJ/AwYHP4CHE2JItcUq\n+CIiMiYp/N9Cb5+B418nKS5vvBz878XwiC2OKROCPT2eiIjINSn816G3z6D0XycpLmuk7YIz+Enf\niyFVwRcRES+h8A9Db59B6ZGTFJc30Hahi8BxASR9P4bUeQq+iIh4F4X/G/T0Oq/w/1p+OfhLvz+D\nVFssk8MUfBER8T4K/xB6eg1Kj5yguKKRM/3BXzZ3BqnzYpmk4IuIiBdT+K8wEPz95Y2c/aqLIAVf\nRER8jMKPM/gfHDlB8RXBT75/Binz4pgUGuTp8URERNzGr8Pf09vH+zUn+WvF5eCn3B9L8rxYBV9E\nRHySX4Z/IPjF5Q2ca+8mKDCAlHmxpNwfy0QFX0REfJhfhb+nt4+D1Sf4a0WjK/ip82JJVvBFRMRP\n+EX4u3v6OFjjDP759m6CA62k2vqDf5OCLyIi/sOnw9/dc/kK/3yHM/iP2OJIvn8GExR8ERHxQz4Z\n/u6ePv5RfYIDA8EPUvBFRETAx8Lf1dPHwaovOXCoyRX85Q/EsWyugi8iIgI+Ev6unj7+0R/8C1cE\nP/n+WMLGB3p6PBERkTHD68NfcPAYe0qOcqGzh5AgK2nz41g2V8EXEREZiteH/5XCj/uDH8+yuTMU\nfBERkW/g9eF/8n/nMGv6RAVfRERkGAI8PcBQWltbWb169bDum/pAvKIvIiIyTGMy/K+88grTp0/3\n9BgiIiI+Z8yF/80332TFihUEB+sYXBEREXcblef4a2pq2LFjB6+99hqGYbB161aOHj1KYGAgzz33\nHLGxsbz44os0NjbS1tbGJ598wpEjR3jnnXdITk4ejRFFRET8woiHPzc3l8LCQkJDQwEoKSmhp6eH\n/Px8ampq2L59O7t27SI7O/uq79uwYYOiLyIi4mYjvtUfFxfHzp07MU0TgMOHD7Nw4UIA5syZQ21t\n7ZDf96tf/WqkRxMREfE7Ix7+ZcuWYbVaXbc7OjoICwtz3bZarRiGMdJjiIiICB54H39YWBgdHR2u\n24ZhEBBwY39/REZOuNGxRK6iNSXupjUlY8Wov6o/ISGB999/H4Dq6mpmzpw52iOIiIj4rVG74rdY\nLAAsXboUh8NBZmYmAM8///xojSAiIuL3LObAq+5ERETE5425D/ARERGRkaPwi4iI+BGFX0RExI/4\nZPjLy8vZvHmzp8cQH1FeXs4vfvELnnrqKerq6jw9jviA2tpannnmGX7+85/T1tbm6XHEB1zPqbY+\nF/6mpibq6uro6ury9CjiIy5dusS2bdv48Y9/jMPh8PQ44gO6u7vZtGkTixYtoqqqytPjiJczTfO6\nTrX1ufDHxsbyox/9yNNjiA9ZsmQJnZ2d/OlPf2LVqlWeHkd8QEJCAseOHePVV1/lzjvv9PQ44uXy\n8vKu61Rbrwp/TU0NWVlZgPMT/7Zs2UJmZiZZWVk0NTV5eDrxRsNZU2fOnGHbtm1kZ2cTHh7uyXHF\nCwxnTR05coTZs2eTm5vLH/7wB0+OK2PccNZTeXk5+fn5rlNt/5tR/8jeb2u4p/yJDNdw11ROTg5n\nz57l17/+NUlJSTo1Uq5puGuqs7OTTZs2ERgY6PowM5GvG+56eumll4Dhn2rrNeEfOOVvw4YNwH8/\n5e+FF14Y9RnFuwx3TeXk5HhsRvEuw11TNpsNm83msTnFO1xv94Z7qq3XbPXrlD9xN60pcTetKXGn\nkVpPXhP+rxuJU/7Ev2lNibtpTYk7uWs9ee0K1Cl/4m5aU+JuWlPiTu5aT17zHP8AnfIn7qY1Je6m\nNSXu5O71pNP5RERE/IjXbvWLiIjI9VP4RURE/IjCLyIi4kcUfhERET+i8IuIiPgRhV9ERMSPKPwi\nIiJ+ROEXERHxIwq/iJc4fvw4s2bNoqys7KqvP/TQQ5w4cWLYP8Nut4/EeNflmWee4eTJkwCsXbuW\n06dPe3giEf+h8It4kXHjxrF58+arDurwRocOHXKdKrZ7924iIyM9PJGI/1D4RbxIVFQUCxYsICcn\n57/e9+WXX2b58uXY7XZycnJcoe3o6OCnP/0pK1asYP369bS3twOQk5NDeno6GRkZ7Ny503XfjRs3\nkpGRwcqVKykuLgZg7969ZGVlYbfb2bp1KwsWLKC3txeAo0ePsmLFCgB++9vfsmbNGpKTk8nMzKS1\ntZXdu3fT0tLCE088wblz51w7FoZh8Oyzz5KWlobdbic3Nxdw/pHw+OOPs27dOlJSUli/fj09PT20\nt7ezdu1aMjIyyMjI4O9//7t7H2wRH6Xwi3iZDRs2UFpaOmjL/0oHDx7kvffe4y9/+QsFBQU0NjaS\nl5cHQHNzM0888QSFhYXExMSwa9cuTpw4wQcffMC+ffvIz8+nqamJ7u5ufv/73zN79mz27t3L66+/\nzssvv8wXX3wBQEtLC/v27WPr1q3cc889lJaWAlBcXEx6ejpNTU3U19fz1ltv8c477xAXF0dRURFr\n164lKiqK3bt3M3nyZABM0yQvL4/m5maKiorYs2cP7777LgcPHgSgqqqKLVu2cODAAU6ePElpaSkl\nJSXExMSwd+9eXnjhBSorK0fyYRfxGQq/iJcJCwtj27Zt37jlX1FRQVpaGkFBQVitVlavXk1FRQUW\ni4U77riDu+++G4D09HQqKiqIjo4mODiYH/7wh/zxj38kOzuboKAgysrKyM/PZ+XKlTz66KNcvHiR\nY8eOYbFYuOuuu1xngaenp7t2A/72t7+RlpZGbGwsGzdu5K233mL79u1UV1fT2dl5zd/r0KFDrFq1\nCovFQkhICHa7nfLyctfM0dHRWCwWbrvtNs6fP899991HSUkJ69at46OPPuLJJ5908yMt4psUfhEv\nlJiYSGJiItu3bx/y/03T5MqDN03TdG3FW63Wq75utVqxWq3s2bOH7Oxszp49y5o1a2hoaMA0TXbs\n2EFBQQEFBQXk5eWxYMECAEJCQlw/Z8mSJXz44YdUVlYybdo0oqOjqa2t5fHHHwcgJSWFpKQkvukw\n0K/PbBiGa+agoCDX1y0WC6ZpEhcXx4EDB7Db7VRWVvKDH/xg2I+fiD9T+EW81MaNG3E4HLS0tAz6\nP5vNRnFxMV1dXfT29vLnP/8Zm82GaZrU1dXx6aefAvD2228zf/586urqePTRR5k7dy4bN27k9ttv\np76+HpvNxptvvgk4t/ZXrVrFqVOnBgU8KCiIhQsX8stf/pL09HQAKisrmTdvHmvWrOG2227D4XC4\nXmcwbtw4V9SvnLmgoADDMLh48SL79+93zTyUvLw8XnrpJVJSUtiyZQtnzpxxvV5BRK5N4RfxIhaL\nxfXvgS3/vr6+QfdbvHgxixcvZvXq1aSlpRETE0NWVhYAt956K7/5zW+w2+2cP3+en/zkJ8yaNYt7\n772XtLQ0MjIyiImJYdGiRaxbt45Lly5ht9t57LHHeOqpp5gxY8ZVcwxIT0/n888/Jzk5GYDU1FTq\n6upYuXIl69ev58EHH+T48eOu+dauXeu6bbFYWLNmDdHR0aSnp7Nq1SoefvhhkpKSBv3eA7ftdjv1\n9fXY7XaysrL42c9+RlhYmBseZRHfZjG/ae9NREREfIqu+EVERPyIwi8iIuJHFH4RERE/ovCLiIj4\nEYVfRETEjyj8IiIifkThFxER8SMKv4iIiB/5f/yMM1GciS/pAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.loglog(N_obs, times1, label='Naive Implmentation')\n", "plt.loglog(N_obs, times2, label='FFT Implementation')\n", "plt.xlabel('N observations')\n", "plt.ylabel('t (sec)')\n", "plt.legend(loc='upper left');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For fewer than 100 observations, the naive implementation wins out, but as the number of points grows, we observe the clear trends in scaling: $O[N^2]$ for the Naive method, and $O[N\\log N]$ for the fast method. We could push this plot higher, but the trends are already clear: for $10^5$ points, while the FFT method would complete in a couple seconds, the Naive method would take nearly two hours! Who's got the time for that plot?" ] } ], "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.3.5" } }, "nbformat": 4, "nbformat_minor": 0 }