{ "metadata": { "name": "", "signature": "sha256:bc886572cb397b961f99cb1ee3d91d14ff0e583f6d8ee074a5a19c03e9f6122e" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "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", "collapsed": false, "input": [ "%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()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "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", "collapsed": false, "input": [ "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');" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAFVCAYAAAApGgzgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuQpNV53/HfzOzNoN7RFGmZuJyCcoSfWrlYZQAFtEvQ\nJRHEUrBB60qtyrFugBAgbYx3wZgC4gjFwrUsNistJoVWXlROoKwgEktUkCoWFjCrEtiBgEroiHWV\nVJUqxRmT0c7Aam8znT+6e6anp69vv5dzzvv9VG1VX7a7z5x++33O5TnnHavVagIAAGEYL7oAAABg\ncARuAAACQuAGACAgBG4AAAJC4AYAICAEbgAAArJulBeb2cWS7nHOvaft8ZslXSNptvHQ9c65H47y\nWQAAYITAbWa3Svo3kl7v8PQFkn7LOfdC0vcHAABrjTJUfkTSByWNdXjuQkm3m9kzZnbbCJ8BAABa\nJA7czrmvSjrd5elHJF0v6b2SLjWzDyT9HAAAsGKkOe4e7nfOzUuSmT0haVrSE93+c61Wq42Ndeq4\nAwAQpcRBL/XAbWaTkl4ys7dJOqZ6r/tgr9eMjY1pdnYh7aKgRbVaoY5zQD1njzrOHnWcvWq1kvi1\naQTumiSZ2Yckvck591BjXvspSSck/Q/n3JMpfA4AAKU35snVwWq07rJFCzof1HP2qOPsUcfZq1Yr\niYfK2YAFAICAELgBAAgIgRsAgIAQuAEACAiBGwCAgBC4AQAICIEbAICAELgBAAgIgRsAgIAQuAEA\nCAiBGwCAgBC4AQAICIEbAICAELgBAAgIgRsAgIAQuAEACAiBGwCAgEQZuG954LBueeBw0cUAACB1\nUQZuAABiReAGPMJoEYB+CNwAAASEwA0AQEAI3AAABITADQBAQAjcAAAEJLrAfe+jL+i1+eN6bf64\n7n30haKLAwBAqqIK3Pc++oK+/6O55fvf/9Gcdh+Y0Y//z0KBpQIGQ6MTwCCiCtyvtATtprmFE9r/\n2EsFlAYYXLdG582ff5Z13QBWiSpwA6Hq1uhcOHaygNKUDxvfICRRBe4t506teWyqslG7dmwtoDQA\nAKQvqsC9Z+e0piobl+9PVTZq303bdc7ZlQJLBfTXrdFZOWNDAaUB4LOoArck7dqxVeNj0viY6Gkj\nGN0anesmovuJAhhRdGeFc86uaKqySVOVTfS0ERQanYgR+QPpW1d0AQDUNRudzdsA0El0PW4gFqzr\nBtAJgRvwEJsJ5YcGEkITfOBm/sRvfD/JsJlQPmggIUTBB24ASIoGEkIUZXLa3hu3FV0EYCRbzp1a\n1ROU2EwIQB09biTCEHi22EwoH+y2mC3yB7JB4AY80WwINUeMWNedPRpI2SF/IDsEbsBTbCaUDxpI\n2SB/IDtRznEDwKDY+AahGanHbWYXm9lTHR6/0syeM7PDZnbtKJ/RC/MnAOAn8geykzhwm9mtkh6S\ntLHt8fWS7pP0PknvkvQJM3vLKIXshPmTYs0tHNfcwvGe/4eGFVBe5A9kZ5Qe9xFJH5Q01vb4FklH\nnHNHnXOnJD0r6bIRPqcj5k/8RsMKAPkD2UgcuJ1zX5V0usNTmyUdbbm/IGky6ecgLM1lYjSsAJBg\nmY0sktOOSmr9hiqS1p7F21Srw32pbz+vqhdfnV312FmTm3THxy8e+r3KYph6ueaz35QkHbzj8jXP\n3fngYS3V6rf3P/ay7v7kyoY3ExONAZgxSbW17zs+Phb995Pk77vzwcN6bb4+9dBap836jL3OhpV2\nfVDPa6VVF9Rt+rII3D+QdJ6ZTUl6Q/Vh8r39XjQ7O9wQ6q4d52v3gRnNLZyQVJ8/2XvDtkTvVQbV\namWoellcrEfd9te0D4G/+OqsPvz7T2rXjq065+zK8uu2nNN5569PXX1+1N/PsPUs9a7Te65/p6SV\n76F9rXcZJanjflqPd+o43Trudi4pu1EaMmms465Jkpl9yMyua8xr/46kb0g6LOmgc+4nKXzOGknm\nT9jxazSDDoGTmDI4phWKt/fGbaUO1AjLSD1u59yPJG1r3H6k5fGvS/r6SCUbQKf1l7SW/bFrx1bd\n/fDzy7eRLY59oBzYOQ1DGWZtJokpg2G9K4BhELgxlD07p7VuYmUF4LqJMYbAR8S0AmLGNET6CNwY\nyr2PvqDTiyvp4qcXa6zPTgHrXQEMisCNNXrteNYrkYqd0pJjWiFfJKkiZARurJJ0x7OFYyc7vu70\n4lJmZS0bGkaIRbeGEw2qwQQfuJk/SVe/pUndEqlah89bX7dw7GT6hSwhtpDNB40jhCD4wN2q34+O\nH+XouiVStW9Y3zT5po00rFLQr0HFsT06GkfFoac9nGgCd78fHT/KwQyyNKnTbZY0FYdjOx1shINQ\nRBO4+/3oBv1Rlr3lN8jSpHPOrixnQDcfZ0lTtno1jAg4QLlEE7iRnkGWJjUzoId9HZKhYZS9YUaN\nyt7AR7GiCdz9fnQM5Q4u6dIkljSNpl+iZbeGEcf2cLrlA9A4QiiiCdz9fnT8KBG6bg0jju3B9csH\nYNQoe90aTvNvnCTBckDRBG6p/4+u/XmGuxALAs5g+uUDMGqUrW4Np6Ovn9Cplj0fSLDsLarA3e9H\nx48SseLYRgi6NZxOL3XeB4IEy85Guqwnyou12QjRlnOnVvX4JPIBEJ6oetyjYhMLIB6dpsLIByhW\np0TKbmhQdUfgbmATi3SwBS18N2o+AA385NobTt2MjYkGVQ8E7gY2sViNAIxYjZIPQAN/dK0Np1/6\nh5s7/p/r/tXbci5VWJjjxrLmsCIB2198N9nrVce9Gvj7btqeZbGi0Ww4SdIdH7lIuw/MaG7hhKR6\nMJ+qbNIlv3J2kUX0XnSBu9+Jrfl8c7ireZukFYSOoI4Q7dqxVXc//LwkqXLGhoJLE4boAvcgOg13\nTVU2qnLGei0cOyVpJWkFAJpo4KevtQeOwZRyjrvbcJckNrEA0BVZ6fBBKQN3N+smxtnEAkBP7FKH\nopUycHNRBiBu/ZZsjbJqgl3qULRSznHv2Tm9KpOR+WwgHt2WbO3asXXkQMu1DdLRrdFEguVgStnj\nlhjuAmJ0ywOH1ySPSeXekwHxKW3gZrhrNXaDAoAwlHKoHKtlObQI5G39xPiqS0RK5LD4jiHy4RC4\n25TxAGI3KMRk85kbtFSrkcMSOHZy7I7ADSA6rbtxZdHTJpigSKWd48YKlsfFo9OlLMuIHBbEjMAN\ndoMCBkAC52o0EotD4IYklschfFkGVi7nCZ+Ueo6beaoVrRv909NGaDoF1vGx9K42RQInfBJMj5th\nmfjwnSItnQLrUk1aOHaygNIA2QomcANAUUjgzBf5BL0RuAEEr1tgvfMj70jl/UngzA/5BP0RuIFI\nlLmXkkdgJYEzH73yCVBH4AYiQC8l+8DK2vAVrY3E6/f+FbkqOSNwY9ko1yhGseilEFjz0t5IPLW4\npLmF46k1Eskn6I/ADWSM7PnsULf565bBn1YjkXyC/oII3GWeu4sV32m66KUgJuQT9OZ94GbuLj58\np+mjl4K8dGokph1gmfboLVHgNrNxM3vQzA6b2VNm9o/bnr/ZzL7XeO4pM/vlpAVk7i4+fKfZoJeC\nPLQ3EsfHRIDNWdItT6+StME5t83MLpa0r/FY0wWSfss5xxgokBO2rc0eyZt1rZdNTWtbWQwu6VD5\ndklPSpJz7ruSLmp7/kJJt5vZM2Z22wjl82LujgSYdPnwnSJOrIzIR7ORODE+rqNvnCRXJWdJA/dm\nSfMt9xfNrPW9HpF0vaT3SrrUzD6Q8HOYu4sQ3ykQvvk3TurU4tLyfXJV8pN0qHxeUutZdtw5t9Ry\n/37n3LwkmdkTkqYlPdHrDavV7iftu669RLv/+NvLt3v93yxMTIxJ6l3GEPhU/qK/0yy1/y15Hj+x\nHKv9VKsV3fngYb02f1yStP+xl3X3J+lpp6nfMdQatJvmFk7oC4+/rEN3XTHy55flWE4iaeCekXSl\npK+Y2SWSlrOKzGxS0ktm9jZJx1TvdR/s94azs91baZMbJ5bn7iY3TvT8v1lYXKxJ6l1G31WrFa/K\nX/R3mpX2er730Rf0f+d+Jkn63c8/rT07pzP9/BiO1X6q1Yp+9/NPr1qZ8OKrs/rw7z+pXTu2MnKT\nglHOF0tLtVSOv3uuf6ekeI/lURokSYfKH5d03MxmVE9Mu9nMPmRm1znnjkq6TdJTkp6W9D3n3JOJ\nSwgEimVv2WFlQvHWT6wNH+Sq5CNRj9s5V5N0Q9vDP2x5/hHV57mB0uoVXPbdtD2TzyQxC3nZfOYG\nzS0c11J9kGc5VwXZ834DFgBox8qE4u29cZvu/Mg7Cts7oMyrfbwJ3L5+CWzNiaQILtlhZYIf2OGs\nGN4Ebh8xR4lREFyyxU5xKKukWeW5K2Lurog5yjIpw3xs6w5TBJd0+bhTXHPUsAzHNooTTOAGQuRj\ncAEQNobKe2COEgDgGwJ3D8xRAgB8Q+DugwQYAOiuiAu7lH21jxeBu7nnsI9fAssdEAtfl1wCw2C1\njyeB+8VXZ5dvl/FLAAAMhu1uPc0qZ8kVYsLSoOz4VLfN4dvm7awvKIPy8qLHjWIwdAqkowzDt76c\nL1jt42ngLtuXACBsDN/mh9U+ngTusyY3Ld8u45cAABhc2Vf7eBG47/j4xaX+EhAXX4YUkR+Gb/NV\n9tU+XiSnvfUX3+z1tpA+JcDEjH2eEao9O6e1+8CM5hZOSOLa1MiWFz3udvRYgHSVfcOKPJR9+Bb5\n8TJwA0hPGTKefVD24Vvkh8BdUvTAyoOMZyAuBO4EQh/KpwcGYBg09P1C4C4hemDlQsYzRuFrQ7+I\ni5v4wpvAXeYvAcgSG1ZgFDT0/eNN4EZ+fOyBxTIU5+uV7sh4BuLhXeCO5QTuM996YL4OxQ3r3kdf\n8PZKd2Q85yPGkUMfG/pl51XgjuUEHgKfemCxDMXF8ncArXxr6MOzwM2JLz/0wAAMyqeGPjwL3Cin\nWIbiYvk7gHY09P3iVeAO4cTHHHz6YhmK27NzmivdAcicV4Hb9xM4c/DZiWUojivdAciaV4Fb8vsE\nzhx8dmIZimte6c7Hv6M14zn03f+AMvPisp6tmifw5m0AALDCu8Dtsy3nTq0aKpf8m4MfRmzrTQGg\nDAjcQ9izc1q7D8xobuGEpJU5eACIHQ19f3g3x+07n+fgAQDxo8c9JObgAQBFInADKWNIEe2aGfwc\nG0iDl4Gbg7uc+N7z0dxEqHl7z87pgksEYBjMcQMlwiZCQPgI3ECJsIkQED4CNwAAAfFyjtt3zMUi\nVLFtIgSUET1uoER8v5APgP4I3EDJsIlQvrgUMNJG4AZKJpYrsYWALH5kIdEct5mNS3pA0lZJJyRd\n65z725bnr5R0p6TTkr7knPtiCmUFgKD0yuLnOgdIKmmP+ypJG5xz2yTdJmlf8wkzWy/pPknvk/Qu\nSZ8ws7eMWlAAAJA8cG+X9KQkOee+K+milue2SDrinDvqnDsl6VlJl41USpTOLQ8cXt4mEgjVlnOn\n1jxGFn+6yniuSLocbLOk+Zb7i2Y27pxbajx3tOW5BUmT/d6wWmWuLWsh1fHExJiksMrcFEKZQ65f\nKZxy/+GnL9NHP/MNvXa0vsXsWZObdOiuKwou1WBCqePQj+UkkgbueUmttdQM2lI9aLc+V5G0dqKn\nzewsyRpZqlYrQdXx4mJNUnjHRSj1HGr9SuHUcdOnrj5fdz/8/PLtEMoeUh2HeiyP0tBIOlQ+I+n9\nkmRml0hq3S/xB5LOM7MpM9ug+jD5dxKXEAAU7pAoWfxIW9Ie9+OS3mdmM437HzOzD0l6k3PuITP7\nHUnfUL1hcNA595MUygogJez+B4QrUeB2ztUk3dD28A9bnv+6pK+PUC4AAKKR5jXZ2YAFACAp3OmI\nsiFwAwAQEAI3vMPezgAGUdZzBYEbXmFv52IxVJqNvTduIyEwZWU+VxC44ZVeezujvFp7Vnc+SMMC\n5T5XELgBeK29Z/Xiq7Ol6VkBnRC4IxXqkCd7O6NdmXtW6K7M5woCN7yyZ+e0piobl+9PVTZq303b\n2XEKyFho0xEhnSvSTqIjcMM7u3Zs1fiYND6mUrSe0VuZe1Z5CXU6IoRzRbckuit3/7cLkr4ngRve\nYW9ntGrvWZ01ucnbnlWoQp2OCOFc0a1uJf1F0vckcAMjCjWfoJ3Pa2Jbe1Z3fPzioosDFIrADcD7\nNbGtPau3/uKbiy5OdJiOyE63upX0a0nfk8ANINihUqSD6YjsdEui+9q+X/+fSd+TwB0hn4c8AfiJ\n6YjspJ1ER+COTLchzyP/+6cFlgq+Y6gUTEdkJ+0kukTX44a/vt9lyPOzX/qu9t4Qzl7JnfZ1TvN6\ntlhtz85p7T4w08x2XR7OQ/o4jtNVxnqkxw1AUhhrYpG+WFZFlAk97sisnxjXqcWlVY9NVTYyZ4W+\nmsN5zdu+KWPPCuiEHndkNp+5QeNjK/ebQ57MWWWDREAAeSNwR6hyxgaGPHPQKRHwo5/5hjdrnwHE\niaHyFPmSdLJuYtzrIc9YdFr7/NrR49r/2EskdiFIRZ+7YpZm3dLjBgAgIATuiLTOt86/cbLo4qTK\nx7nkTmufz5rcxPQEuvLxOEb6ss7UJ3BHon2+9dTikuYWjkcx3+rrPtqdtjI8dNcVQU9P7L1xG8Ol\nGfHxOKYhESYCdyQ6zbcu1RTFXtM+76PN2mcMyrfj2MeGBAZD4AZGEML1gEPFxiDZ8q0hkYVYjyEC\ndyRi3ms65r8N5cFxjLQQuCPR7dJxMfQCY/7bUB6+Hcc0JMJF4E6JD0keMc+3xvy3oTx8Oo59a0hg\ncATuFPiS5BHzfGvMfxvKw7fj2KeGRCzy6MQRuFNQhiQPAPHxrSERurw6cQRuIIHWbFXWPgOQ8uvE\nEbhTQJIHkC4fckYQtpiPIQJ3CkjyANLjS84IwlXUMZRXJ47AnRKSPIB0kDOCURV1DOXVieOynilp\nJnk0bxcl5rnWmP82lAfHcdx27diqux9+fvl2FgjcALyy5dypVcOcEjkjWYqxIVHkMZRHJ46h8gjE\nuh8vyomcEYwq9mOIwA0MKeZsVV+QM4JRxXwMEbiBIZDxnA82BsGoYj6GCNzAEMh4BlA0ktNSFGOS\nBwDAL/S4gSGwSx6Aog3d4zazn5P0Z5KqkhYkfcQ59/dt/+d+Sdsbz9ckXeWcmx+9uECx9uyc1u4D\nM5pbOCFpJVsVAJqyHn1N0uO+QdL/cs5dJunLku7o8H8ukHS5c+49zrn3ErSzQ4Zz9tqX28WcrQrA\nf0nmuLdL+sPG7Scl3dn6pJmNSzpP0kNm9vOSDjrn/nSkUqKjbhnOu3ZsjS6Lsl3rlbny5ssuebEr\nW85Ikcd0rGKty56B28yukfTbbQ//naRmD3pB0mTb82dI2i/pvsb7P2Vmf+2ce7nXZ1WrnACH9cqP\nO2c4f+Hxl3XorivWPBdTHU9MjEnK52/q9Fm9Pj+mevZVjHWc5zE9CF/KgbV6Bm7n3EFJB1sfM7PH\nJDW/0Yqkn7a97Jik/c65443//y1Jb5fUM3DPzrIOdmi1zg8vLdXW1Ge1WomqjhcX6398Hn9Tp8/q\n9vmx1bOPYq3jPI/pfmKsY99GNEZpGCWZ456R9P7G7V+V9HTb8ybpWTMbN7P1ki6V9DeJS4iuyHAG\nAD/kufV0ksD9J5J+xcyekXStpH8vSWZ2s5ld6Zx7RfWkte9IekrSocZjpZbFlxr7frwoHvvgA/4Z\nOjnNOfczSf+6w+N/1HL7PtXnuJGxPC4hBwDwBzunBY4M52w1l9s1b+/ZOS3Jn3kyAOXDzmkITl5r\n17mgCPKS1zHN1EccCNwISp7BlAuKIA80EDEsAjeCQjBFbDimsxfbDpMEbkTh6OsnUn9PltsB4ctj\nRCPvhgGBOwextfaK1CmYjo9JlTM2pP5ZZV9ux3GbDxqI2cp6RKOIqQ4Cd8by+FL33ritNFnOnYLp\nVGWT1k1kcyiX9YIizLvmp+wNxNAVMdVB4M4Y81fpyzOYNpfbTVU2lepEynGbr7I2EPMQ44gGgRvB\nKWswRbw4prOT9YhGEQ0DAneGbnngcMch3NBbe4hfjL2UsitzzkKWIxpFTHUQuDO2+cwNzF8hOMy7\nxqXsOQtZj2jkPdVB4M4B81fZKXMvImt5Hrfs6JUtchaylfdUB3uV54D9xLMx/8bJ5X3EpZVexK4d\nW1Op57IHEo5bwE/0uBGsU4tLax6jFwGsRc5CXOhxI0h7b9yma+75lmo5fRaQtSyPsz07p7X7wIzm\nFuo7DDZzFhAmetwIFr0IYHDk2sSDwJ2RNJOmSNzpjMxnxCbL3zprxePBUHkGOi29yGo/7bLbtWOr\n7n74+eXbANBJTFNeBO4MdFp6sVSTxsfHCihN3Mh8zlbrya7ZE0z7BNgcnWre3rNzOtX3B/KQZ8OA\noXKgA9aH56PsG4MASRC4M5Bl0hTz3dkjmKSv23HLxiDA8AjcGSBpKjytgYVgAsBnBO6MsPQC6I8l\nffnae+O2qJK0yorAnRGWXoSLYJIfRqdWkFcRBh+mKwncnuPH3F/avQiCSb4YnSKvAsMhcHus0495\nbuG4TnfYoxvpIpislmUDktEp8iowHAK3x7qtB184drKA0pQLwWRFVr1BH4YcgRARuAO0VBPD5sjN\nqL1Bpnv6yyKvgobR8EKpMwK3R9oPmk4/5ibmwNJDYMkOc7eDIa8CwyBwZ2jUpKn2H3M75sCSaW0g\nEVj6G6U3yNzt4IbJqwilZ4hsELg91/wxIxsElv7oDeaDvAoMisDtueaPef3E2q+KtcXIS9Is+0F7\n62wMghD4Mq3G1cECsfnMDVqq1TS3cELSSq8Ho9ly7tSqoXJpJbDQ61nR6ypsva4atmfntHYfmOG4\nRfC6TasVca6gxx0Q1hanj2Hg7HU6bn3puQCD8mlajcAdEObAskGDKFvtxy0JgdmjYTS4ZqJfSHVG\n4C5Qe3ZzKAdNbGgQ5cunnkuMaBgNb/6Nk33rzKdrGBC4PcAPDYA0WJJev0Y+DaPhneqwjXR7nfk0\nrUbg9kC/HxoZt8l0WuvKyEbxfOq5hIZGfrF8mVYjcKM0OOn5waeeS2gG6U3TMBper+W2rR0AX6bV\nCNwe4IeWD4YQR5PmyI8vPZcY0TAa3uYzNwRVZwRuDwzzQ2PYHD4ZdOqh/bj1pecSmkEb+TSMhhdS\nnRG4PRHSQRMqRjbSxdRD/gZt5NMwGkxrw/Mrf3UkmDojcHuCH1r2ep30GMkYHlMPxejUyOeiI8Pr\n1PCcWziu0x0yzH2TOHCb2dVm9p+6PHedmT1vZt8xsw8kL168yG4uBiMb+SCQpK9ZpzTy09Gp4blU\nkxaOnSygNMNJFLjN7H5JfyBpzXWrzOxsSZ+WtE3SFZI+Z2YbRilkbLoNMYbQ0gtFt4YRJ730MPWA\nMvJhdC5pj3tG0g3qELgl/VNJM865U865eUlHJPFLbtFtiDGEll4ImHvNB9nLCFmnhuf4mFQ5Y6Wf\n6evIaM/AbWbXmNnLbf8udM79eY+XVSQdbbm/IGkyjcICg2DuNT9MPSBUnRqeU5VNWtdY0+1zB6Dn\nZT2dcwclHRzyPedVD95NFUlrz6RtqtXytNLffl5VL746u+qxsyY3aWmppomJsczqojR1PCaptvbh\n8fF63U5M1AeKqOfRVasV/YM3/5wk6aLzf2H58UHq+NC/u2Kkzy2b1jptr99O9T3qcV6GOr7r2ku0\n+4+/vXz7c4eek1T/21/5cecOwBcef1mH7kp+7KYhi+txPyfpP5jZRkmbJG2R9L1+L5qdLb4Vk5dd\nO85fc43ivTds0y0PHNbiYi2TuqhWK6Wp4y3ndL7G9qeuPl+zswtaXKxHdeo5HZ3qkzpOX2udttdv\np/q+5/p3rnlsUGWp48mNE8vXmZ/cOLG6zjo0/iVpaSmdc/QoDaNRloPV1PKnmdnNZnalc+7vJO2X\n9Iykv5R0u3OOyds2nYYYfUh6iAFzrygTX+dhQ+dz8mXiHrdz7tuSvt1y/49abn9R0hdHK1rcmtnN\nzdtI164dW3X3w88v325F4yhbzUDSvL1n53TBJQpfe502j+Fu87C7dmzlvDKiPTun14yM7rtpe8Gl\nqmMDFkSJZV/F8DmhJ1S96pREzGz5mnxJ4AaQGgJJ+qjT4vjaAcgiOQ1AiTD1UJwt53ZOxPSpd4j0\nEbgBpIZAkr5edXrO2RVv52FDEWLDk6FyAKkhoz99/erU13lYZIfAXSCWfyFGBJL09apTX+dhkR2G\nyhEtGkXFYKlj+qhTtCJwAwDQhY8dAIbKAQAICIEbAICAELgBAAgIc9wAEDgf52GRHXrcAAAEZKxW\n63LR0XzVynDt1yKV5fq6RaOes0cdZ486zl61WhlL+lp63AAABITADQBAQAjcAAAEhMANAEBACNwA\nAASEwA0AQEAI3AAABITADQBAQAjcAAAEhMANAEBACNwAAASEwA0AQEAI3AAABITADQBAQAjcAAAE\nhMANAEBACNwAAASEwA0AQEAI3AAABITADQBAQAjcAAAEhMANAEBACNwAAASEwA0AQEAI3AAABITA\nDQBAQAjcAAAEhMANAEBACNwAAARkXdIXmtnVkn7DOfebHZ67X9J2SQuSapKucs7NJy4lAACQlDBw\nNwLz5ZJe6PJfLpB0uXPu/yUtGAAAWCvpUPmMpBskjbU/YWbjks6T9JCZPWtmHxuhfAAAoEXPHreZ\nXSPpt9se/qhz7s/N7N1dXnaGpP2S7mu8/1Nm9tfOuZdHLSwAAGXXM3A75w5KOjjkex6TtN85d1yS\nzOxbkt4uqVfgHqtWK0N+DIZFHeeDes4edZw96thfWWSVm6RnzWzczNZLulTS32TwOQAAlE7irHLV\ns8VrzTtmdrOkI865r5nZlyV9R9IpSYecc6+MVkwAACBJY7Varf//AgAAXmADFgAAAkLgBgAgIARu\nAAACQuAGACAgo2SVj6yxy9oDkrZKOiHpWufc3xZZphg0luF9SdI5kjZK+qykVyQdkrQk6XuSbnLO\nkZk4IjN7i+rLHf+56nV7SNRxqszs9yRdKWm9pC+ovnPjIVHPqWich78o6ZdVr9PrJC2KOk6FmV0s\n6R7n3HsZuNiFAAAC30lEQVTM7K3qUK9mdp2kT0g6Lemzzrkner1n0T3uqyRtcM5tk3SbpH0FlycW\nvylp1jl3maR/KemA6nV7e+OxMUm/XmD5otBoIP1HSW+oXqf3iTpOVWOHxnc2zhHvlvRL4lhO2+WS\nznTOXSrpM5L+QNRxKszsVkkPqd6BkjqcI8zsbEmflrRN0hWSPmdmG3q9b9GBe7ukJyXJOfddSRcV\nW5xofEXSXY3b46qvp7/AOfd047H/LulfFFGwyOyV9CeSftK4Tx2n73JJL5vZf5X0NUl/IelC6jlV\nP5M0aWZjkiYlnRR1nJYjkj6olet6dDpHvEPSjHPuVOMqmkdUH4XuqujAvVlS6+U+FxvDNhiBc+4N\n59zrZlZRPYjfodXf9euq/0CRkJl9VPVRjW82HhrT6ovuUMfpqEq6UNJvSPqkpP8s6jltM5I2SfqB\n6iNI+0Udp8I591XVh7+bWut1QfV63SzpaIfHuyo6SM5Lat0Qd9w5t1RUYWJiZv9I0rckfdk594jq\ncypNFUk/LaRg8fiYpPeZ2VOS/omkh1UPMk3UcTr+XtI3nXOnnXM/lHRcq09q1PPoblW9x2eqH8tf\nVj2foIk6Tk/reXiz6vXaHgcrkuZ6vUnRgXtG0vslycwukfRSscWJg5n9vKRvSrrVOXeo8fALZvau\nxu1flfR0p9diMM65dznn3u2ce4+kFyV9WNKT1HHqnlU9T0Nm9guqX33wL6nnVJ2plZHPOdWTljlf\nZKNTvT4n6Z+Z2UYzm5S0RfXEta4KzSqX9LjqvZaZxn2u3Z2O21XvldxlZs257n8raX8j6eH7kv5L\nUYWLVE3SbtWvQ08dp8Q594SZXWZmz6ne0bhR0o9EPadpr6Q/NbNnVO9p/57qKyWo4/Q0M/LXnCMa\nWeX7JT2j+jF+u3PuZK83Y69yAAACUvRQOQAAGAKBGwCAgBC4AQAICIEbAICAELgBAAgIgRsAgIAQ\nuAEACMj/B3JGOmCPQUEPAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 2 }, { "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 \u2013 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", "collapsed": false, "input": [ "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)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's use the ``gatspy`` tools to plot the periodogram:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "t, y, dy = create_data(100, period=2.5)\n", "freq = freq_grid(t)\n", "print(len(freq))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "750\n" ] } ], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "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);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFVCAYAAADc5IdQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8HOd95/HP7C4WvQMEu0iR4ogqtLpkFVuWJTnyWTk5\nsZPIsS9W4rvESa44d+f4fInvLnfJ+eKX7BTHsS3Zlotsx5IlV5nqhaIkUiwiJZZhbyAIgCB62Tr3\nx+wMdhdYLAgC3MHo+369+CK3YDEccvHd3/P8nmcM27YRERGR0gmV+gBERETe7hTGIiIiJaYwFhER\nKTGFsYiISIkpjEVEREpMYSwiIlJi0wpj0zSvN03z+Unuv9s0zc2mab5imuYnZv/wREREgq9oGJum\n+WngAaA87/4y4IvAHcC7gX9nmuaCuThIERGRIJtOZXwA+A3AyLt/LXDAsqx+y7ISwMvAu2b5+ERE\nRAKvaBhblvUYkJzkoTqgP+v2IFA/S8clIiLythE5h6/tB2qzbtcCvVN9gW3btmHkF9giIiKBVjT4\nziWM9wIXmabZCAzjDFF/YcqjMQy6uwfP4VvKdLS21uo8zzGd47mnczz3dI7Pj9bW2qLPOZswtgFM\n07wXqLEs6wHTNP8MeBJnuPsblmV1zORARURE3s6M83zVJlufwuaePu3OPZ3juadzPPd0js+P1tba\nosPU2vRDRESkxBTGIiIiJaYwFhERKTGFsUjAJZIpPv/wNjbv6Sz1oYhIAeeytElE5oFdh3vZd7yP\nfcf7uG5tW6kPR0QmocpYJOBS6XSpD0FEilAYiwRcKn1ely+KyAwojEUCTmEs4n8KY5GAS6UUxiJ+\npzAWCTjNGYv4n8JYJODSGqYW8T2FsUjAJRXGIr6nMBYJuGRKw9QifqcwFgm4REJhLOJ3CmORgIsn\nFcYifqcwFgm4eCJV6kMQkSIUxiIB51bG4VDR65uLSIkojEUCLpF0KuOQwljEtxTGIgEXT6gyFvE7\nhbFIwGmYWsT/FMYiAadhahH/UxiLBJxbGdvaiEvEtxTGIgHnLm3SBSNE/EthLBJwbmWs6xqL+JfC\nWCTgEm4Y67rGIr6lMBYJuPFhahtbE8civqQwFgm47L2plcUi/qQwFgmwtG17w9SgJi4Rv1IYiwRY\nIu+KTUnNG4v4ksJYJMDywzitcWoRX1IYiwRY/uUT1VEt4k8KY5EAi+dVxlprLOJPCmORAJtQGauB\nS8SXFMYiAZY/Z6zKWMSfFMYiAaY5Y5H5QWEsEmDJvEo4rcpYxJcUxiIBlh++GqYW8SeFsUiAueuK\nDcO5rTAW8SeFsUiAuc3TZRHnra5uahF/UhiLBJh7laaycCaM1cAl4ksKY5EAc4elvcpY22GK+JLC\nWCTA3DljL4xVGYv4ksJYJMDSXmUcBjRnLOJXCmORAHMr40jYaafWOmMRf1IYiwSYO0U83k2tMBbx\nI4WxSIB5DVzqphbxNYWxSIBNnDNWGIv4kcJYJMAmdFOrgUvElxTGIgFmp3MbuFQZi/iTwlgkwFIT\nKmOFsYgfKYxFAsydM466c8Zq4BLxJYWxSIC5hXAk002d1naYIr6kMBYJMDt/b+qUGrhE/EhhLBJg\nE7upVRmL+FFkqgdN0wwBXwHWATHgE5ZlHcx6/IPAZwEb+KZlWV+dw2MVkbM04apNCmMRXypWGd8D\nRC3LuhH4DHB/3uNfBO4AbgL+s2ma9bN/iCIyU+n86xkrjEV8qVgY3wSsB7AsaxNwTd7jCaABqAQM\nnApZRHzCzkwR6xKKIv425TA1UAcMZN1OmaYZsizL7QK5H9gKDAM/tixrIP8F8rW21s7oQOXs6DzP\nvflwjssrygBobKhybpdH5sVxu+bTsc5XOsf+UCyMB4DsfykviE3TXA78KXABMAJ8zzTND1mW9ehU\nL9jdPXgOhyvT0dpaq/M8x+bLOR4ejgEwNhoHYGg4Ni+OG+bPOZ7PdI7Pj+l84Ck2TL0ReD+AaZo3\nADuzHqsAUkAsE9BdOEPWIuIT+TtwJTVnLOJLxSrjx4E7TNPcmLl9n2ma9wI1lmU9YJrmt4FXTNMc\nAw4AD83doYrI2cpfZ5xWGIv40pRhbFmWDXwy7+59WY9/CfjSHByXiMwCdVOLzA/a9EMkwNwrJka0\nA5eIrymMRQJMlbHI/KAwFgmwtHbgEpkXFMYiAZa/N7UauET8SWEsEmCqjEXmB4WxSIBNuJ6xwljE\nlxTGIgHmhm/IcDaPd4etRcRfFMYiAeaGbyhkEAoZKItF/ElhLBJg45WxgWGoMhbxK4WxSIDlVMaG\nga0wFvElhbFIgOVWxoa3I5eI+IvCWCTA3OZpw4BQSMPUIn6lMBYJsHTa9qpiAw1Ti/iVwlgkwNK2\nTSjzLg+FDLTMWMSfFMYiAeZWxuCsNVZlLOJPCmORAEvbNkbICWOngUthLOJHCmORAEunGa+MtemH\niG8pjEUCLG3bhL3KWN3UIn6lMBYJMGfO2PmzNv0Q8S+FsUiA5c4Zo25qEZ9SGIsEWG43taFhahGf\nUhiLBJhtj4exYRjYKo1FfElhLBJgqfR4A5c2/RDxL4WxSIClbXLmjNXAJeJPCmORAMvvplZlLOJP\nCmORALNtm5DWGYv4nsJYJMDSdm43tYapRfxJYSwSYKl0dmVskE6X+IBEZFIKY5EAy92bWg1cIn6l\nMBYJMDvresaGNv0Q8S2FsUiA5e/AZduqjkX8SGEsElBp28Yma5g6s8RJUSziPwpjkYBKZxYVZzdw\nZd8vIv6hMBYJKHc4OpS1HaZzf8kOSUQKUBiLBJS7jGn8QhGZ+5XGIr6jMBYJKDd0s7fDBDVwifiR\nwlgkoNL5w9TenHHJDklEClAYiwRUakIDl3O/rX5qEd9RGIsElO2GsZFfGSuMRfxGYSwSUG7mepWx\nuqlFfEthLBJQ3jpjr4Erc7/SWMR3FMYiATXeTZ07TK0sFvEfhbFIQE3cgSv3fhHxD4WxSEAVWtqk\ndcYi/qMwFgmodF43tbc3tcJYxHcUxiIB5XVTu3PGodz7RcQ/FMYiAeVWxkbmXW5omFrEtxTGIgHl\nDkeHJ2yHqTAW8RuFsUhA5c8Za2mTiH8pjEUCKr+bWpdQFPEvhbFIQHlzxl4DlypjEb9SGIsE1Hg3\ntfO7KmMR/1IYiwSUWxlPaOBSGIv4TmSqB03TDAFfAdYBMeATlmUdzHr8WuB+wADagX9jWVZ87g5X\nRKZr4pxxZpg6XbJDEpECilXG9wBRy7JuBD6DE7wAmKZpAF8HPm5Z1i3As8DKuTpQETk7E7qpvU0/\nVBmL+E2xML4JWA9gWdYm4Jqsx9YAPcCfmab5AtBgWZY1FwcpImfPDV23IjbQph8iflUsjOuAgazb\nqczQNUALcCPwj8DtwHtN03zP7B+iiMxEOjMc7c0Zh9w541IdkYgUMuWcMU4Q12bdDlmW5c449QAH\n3GrYNM31OJXz81O9YGtr7VQPyyzReZ57fj/HNSedz9F1dRW0ttZSU1Oec3s+mC/HOZ/pHPtDsTDe\nCNwNPGKa5g3AzqzHDgE1pmmuyjR13QI8WOwbdncPzvRYZZpaW2t1nufYfDjHfX2jAIwMx+juHmR0\nxOmt7O0b8f2xw/w4x/OdzvH5MZ0PPMXC+HHgDtM0N2Zu32ea5r1AjWVZD5im+QfA9zPNXBsty/rV\nOR2xiMwab844fwcudVOL+M6UYWxZlg18Mu/ufVmPPw9cPwfHJSLnaGI3tRq4RPxKm36IBFT+VZvc\nbmotbRLxH4WxSEAVqow1TC3iPwpjkYByC2Aj8y5354w1TC3iPwpjkYDytsM0tDe1iN8pjEUCqnAD\nV8kOSUQKUBiLBJS705a3HaYuoSjiWwpjkYDyKuPMu9ytkDVnLOI/CmORgLLz5oy16YeIfymMRQIq\n/3rGqoxF/EthLBJQExq41E0t4lsKY5GAchu4MoWxt95YWSziPwpjkYAab+BSZSzidwpjkYCyyVy1\nycifMy7ZIYlIAQpjkYByu6YndlMrjUX8RmEsElDj1zN2bqubWsS/FMYiAZXfTW14c8YlOyQRKUBh\nLBJQEy4U4XVTK41F/EZhLBJQtjtnHMqvjBXGIn6jMBYJqPHKmMzvmTDWOLWI7yiMRQIqfztMt5ta\nhbGI/yiMRQJK22GKzB8KY5GAGl/a5DZwqZtaxK8UxiIB5RbA7pt8fJhaaSziNwpjkYAqdAlFDVOL\n+I/CWCSg3DljI2/TD3fJk4j4h8JYJKC8SyiGcjf9UGUs4j8KY5GAstO564wNNEwt4lcKY5GAmjBn\nHNIlFEX8SmEsElATLxSRuV9pLOI7CmORgPLmjPM2/VAWi/iPwlgkoMaHqZ3bXmWsXT9EfEdhLBJQ\n7uYehpE/Z6wwFvEbhbFIQE2cM1Y3tYhfKYxFAsodjTa8Syg6vyuLRfxHYSwSUGnbxjCyhqlVGYv4\nlsJYJKDstO0FMGQNU6uBS8R3FMYiAZW2ba9pCzRMLeJnCmORgEqnya2MQxqmFvErhbFIQDmV8fht\nbfoh4l8KY5GAStu5c8Yhbfoh4lsKY5GASqdtr2kLsq5nrNJYxHcUxiIBZdvj1TBkL20q0QGJSEEK\nY5GAStu217QF43tUqzIW8R+FsUhApQutM1YYi/iOwlgkoOwJDVzqphbxK4WxSEClbXKWNukSiiL+\npTAWCajJhqkNNGcs4kcKY5GAyt8OE5xAVmEs4j8KY5GAyq+MwRm2VmUs4j8KY5GAStvkbPoBbmWs\nMBbxG4WxSEDl700NTke1hqlF/EdhLBJQtm1PUhk71zkWEX9RGIsEVP4lFMGtjBXGIn4TmepB0zRD\nwFeAdUAM+IRlWQcned7XgR7Lsv7bnByliJw1e5JhasPQph8iflSsMr4HiFqWdSPwGeD+/CeYpvmH\nwGWA3uIiPjJ5N7UqYxE/KhbGNwHrASzL2gRck/2gaZo3AtcBXwOMCV8tIiVh2zY2hYapS3NMIlJY\nsTCuAwaybqcyQ9eYprkI+BzwpyiIRXzFrX4nbvqhBi4RP5pyzhgniGuzbocsy0pn/vwhoAV4AlgI\nVJmmuceyrO9M9YKtrbVTPSyzROd57vn5HCeSKQDKyyM5xxmJhDFChq+PPdt8Oc75TOfYH4qF8Ubg\nbuAR0zRvAHa6D1iW9Y/APwKYpvl7wMXFghigu3tw5kcr09LaWqvzPMf8fo5jCSeMk8lUznHaaZtk\nOu3rY3f5/RwHgc7x+TGdDzzFwvhx4A7TNDdmbt9nmua9QI1lWQ/kPVdjXyI+4V6ZabLtMDNFs4j4\nyJRhbFmWDXwy7+59kzzv27N5UCJybtyG6fwwdrbDTE/yFSJSStr0QySA3AauvCwmHDLUwCXiQwpj\nkQAq1E0dMgxSCmMR31EYiwSQXWDOWNczFvEnhbFIALmBm18Zh7UDl4gvKYxFAmi8mzr3/lBo/DER\n8Q+FsUgAeXPGk22HqTAW8R2FsUgAed3U+Q1cIYWxiB8pjEUCqOCmH4aBjXMhCRHxD4WxSAAVauBy\nb6uJS8RfFMYiAWQX2PTDC2MNVYv4isJYJICmGqZ2Hj/vhyQiU1AYiwRQob2pwxqmFvElhbFIAI1v\nh5l7v5vN2hJTxF8UxiIBVGiYWpWxiD8pjEUCqOCFItTAJeJLCmORAHLD1ijYwKUwFvEThbFIAHnr\njAstbdIwtYivKIxFAmiq6xmDKmMRv1EYiwRQoesZu93VymIRf1EYiwSQG7YTd+By3vJa2iTiLwpj\nkQAqPEzt/G4rjEV8RWEsEkCFh6nVwCXiRwpjkQDyKuMCS5s0TC3iLwpjkQDSJRRF5heFsUgAjW+H\nmXt/WDtwifiSwlgkgNzK18hLY0PrjEV8SWEsEkDFLxRx3g9JRKagMBYJoMINXJnHlcYivqIwFgkg\n22vgyr1fDVwi/qQwFgkgb8640DpjVcYivqIwFgmggpt+qIFLxJcUxiIBpHXGIvOLwlgkgAqtM9YO\nXCL+pDAWCaCC3dSqjEV8SWEsEkCFNv1w1xnb6fN+SCIyBYWxSAAV2vTDvalhahF/URiLBFC60Dpj\nQ8PUIn6kMBYJoEJLm3ShCBF/UhiLBJAauETmF4WxSAC5c8LhsDb9EJkPFMYiAeSGccFNPxTGIr6i\nMBYJIDdsI3kdXOMNXOf9kERkCgpjkQAqVhmn0lpoLOInCmORACocxs7vqoxF/EVhLBJA6UzlG8nf\ngctwd+BSGov4icJYJICKD1MrjEX8RGEsEkDe0qZQ/naYWmcs4kcKY5EAShcI47A2/RDxJYWxSABp\nnbHI/KIwFgmgVGryynh8B67zfkgiMgWFsUgAucPQ4fxNP1QZi/iSwlgkgFIpp/SdMEyduak5YxF/\nURiLBFChbmpdtUnEnxTGIgGUTtsYaJ2xyHwRmepB0zRDwFeAdUAM+IRlWQezHr8X+I9AEngT+GPL\nsvQuFymxVNqeEMQw3sClHbhE/KVYZXwPELUs60bgM8D97gOmaVYC/xu41bKsm4F64ANzdaAiMn2p\ntD1hiBrGh61TGqYW8ZViYXwTsB7AsqxNwDVZj40B77QsayxzOwKMzvoRishZS6dtwuGJYeztwKXK\nWMRXphymBuqAgazbKdM0Q5ZlpTPD0d0Apmn+e6Dasqxnin3D1tbaGR+sTJ/O89zz9TkOGUTCoQnH\nWFYRdX6PRvx9/Bnz4RjnO51jfygWxgNA9r9UyLIsb7uAzJzy3wKrgd+czjfs7h4822OUs9TaWqvz\nPMf8fo7j8RQGE99vQ6MJAEZHE74+fvD/OQ4CnePzYzofeIoNU28E3g9gmuYNwM68x78GlAMfzBqu\nFpESc4apJ769vXXGGqYW8ZVilfHjwB2maW7M3L4v00FdA2wBfh94CXjONE2Av7cs6ydzdbAiMj2p\ndNrrnM6mdcYi/jRlGGfmhT+Zd/e+rD+HZ/2IROScpdI20chklbEauET8SJt+iARQoW5qVcYi/qQw\nFgmggpt+6EIRIr6kMBYJoEKbfmiYWsSfFMYiAZQuEMbg7MKlLBbxF4WxSAAVGqYGZxcuXShCxF8U\nxiIBY9t2Zph68rd3KKQGLhG/URiLBIybs1MOU6syFvEVhbFIwKTSzo61hYapQ4ahyljEZxTGIgHj\nzgcXqoxDqoxFfEdhLBIwRcPYUBiL+I3CWCRg3DAuOEwd0jC1iN8ojEUCJq3KWGTeURiLBEwqVWzO\nGG36IeIzCmORgEnZxYapQ6qMRXxGYSwSMOPD1AU2/TDQDlwiPqMwLiHbtmk/PaxmGplVqZSzzniq\npU22/s+J+IrCuIRe2N7OXz64iSdePVrqQ5EAKdZNHdbe1CK+ozAuoR0HewDYuq+7xEciQeKOtBSq\njA0tbRLxHYWxH+jnosyiYt3Uzt7U5/OIRKQYhXEJuT8qbaWxzKKim35onbGI7yiMS8gwMj8s9XNR\nZlHxTT90CUURv1EY+4B+LMpsms6FIkCBLOInCuMS8gpj/UyUWeSFcbjAOmM3jDVULeIbCmORgHFD\nNmQUnjPOfp6IlJ7C2Bf0Q1FmTypdfNMP53n6fyfiFwrjOfb1n+/i4af2ARCLpzjY3u895jZw6Uei\nzKbxYerCS5tAc8YifqIwnmOv7erk2W0nAPjyYzv56+9uZd/xPmB8aZPSWGZTsaVNZRHnbR9PaLGx\niF8ojM+jXUd6AejuG3Xu0MommQPe0qYCc8bRsjAA8WTqvB2TiExNYVwCkUyXq7fph4YLZRYVG6Yu\nj2TCWJWxiG8ojM+T7z+9z/uzV7AUqFxEzkWxYepoWWaYWpWxiG8ojOdQdsX7zNYTJTwSeTspdj1j\nzRmL+I/CeA4VWzoyPkw998cibx/FduAqd+eME05l3NEzzP/61us88PPd5+cARWQChfEccq+ek88N\nX21NLXPBXWdceJjabeBKMxpL8jff3crRzkFe3XXqvB2jiORSGM+hVJHr1BmzsB/mWDypBjDJUexC\nEVFvmDpFd98ow2NJ77FYXPPIIqWgMJ5DyekOU8/w9YfHEvzxF1/inx5/a4avkCuWSPH9Z/bR1Tsy\nK683E3uO9vKNX+wmmdJ85kwVG6bOroyHRhM5j50eGJvbgxORSSmM51ChYerZ0tHjhOa2fd2z8no/\nffkwz2w5wYO/3DMrrzcTX/jBdja+dYpjnUMlO4b56EB7P28e6gHG/99NpzJ2w7i+OgpAT7/CWKQU\nFMZzqNgw9bmWxrHE7A4pHu8cBCCZLE1VOhYfHy4diSWmeKZks22bv/nuVr70ox3Ytu1tc1l0zjiR\nYnDEOc8XLKwFoEeVsUhJKIznUKFuajuTvkbe7bM126HpHkWhzSLm2tDIeADnD59KYd6ObjjnLVVk\nadP4OuPxYeoVbhirMhYpiUipDyDIinVTu3E80/6r2b4EXqn3yh7MCuDh0eQUz5Rse4/1eX8+MxCb\nxjD1+A5cYzFndGXFwjpAlbFIqagynkOFKmM3RNu7z21eNHEOTU77jvfxd4/sYGTMPxVodjWsynj6\n9hzt9f7cMzA2fj3jaezANTgaB2DZghpChqHKWKREFMZzqNCc8WMvHSKZSnOsa+owPnJqgNFY4Qox\ncQ7D1P/8k7fYebCHJzcfH7+zxJd01DD1zBw6OX5Zzp6BMZJFrmc8XhmPzxnXVUdprI2qMhYpEYXx\nHCo0TN07GONwx4B3e7Jh6l1HzvBXD23h0RcOFnz9QpXxtr1dDI7Epzw294fwZK9RqmXL2cPUY1kf\nQl7ddYq/f2RHToPXTI3GkpwJUODYtk3f0Pi/9ZmBMe8DXFXF5LNQXmWccOaMK6JhyiIhmusq6BuM\naVmZSAkojOfQVNth5jZfTXzeqcyypee3txd8jcQkewu/vreL//HAqzycdWGKbB09w3z3ScvruM2u\nnlLeD+HSpPHQ6HiojGV1ij/81D52HOzhideOndPrj4wl+dSXX+Z/fHNz8U73eWI0liKRTGd1Q8cY\nznyoqS4YxrnrjGsqywBorq/ABs4Mxub+wEUkh8J4Dk01J5xdkc40+iarao9llie9vqdr0q957MVD\nOQEfi6eIxVN87We7sDKNQMk5Xh9dSPYwdfZOUKnMB4fewcIV7bNbT/CX39jE7iNnCj5n+/5u4ok0\nw2NJzgzMXuCkbZuTp4fPadpgpvqHnb/HsgU1hEMGZwbGGB5LEo2EKMsMR+dzLxQRywxT11aNhzGo\no1qkFBTGcyQWT/HdpyavTgHGssJmpsPCk/3wH7884+RfczrvB20skWLvsV427e70PhTM9vrl6coZ\npk5knx/nyAZHJp9Htm2bHzyzn/buYZ7bVngkobN3fAlQ9nKgc/X5723jLx7cxE9fPjxrrwlOY9Zf\nf2cL33vKKvic/swQdUNNOY215fQMjDE8lig4RA0QMgzKIiEGR+IkU2lqKp0NP5rrFMYipaIwniPF\nGpCOdAxOev+zW0/wqS+/PGnjVv9QLGcf6skrsamvPpFfTccSKYbzOqrjiRRb9naxZ4oqcy4MjSQw\ncK4q5FbGqXTaG44fzGvwevjpfbz4RjsdPSPesPuuw2cKznlmB/BshfHASJwD7f3e955ND/5iNwdP\nDvDctvaCy9j6h90wjtJUV0H/UJyBkQTVmaHnQqKRkDc6kD1MDVreJFIKCuM5UmzZ0frN4/Of2QH7\n8NP76B+KYx3vy3n+U68f51Nf3sjuI+PLWCYLY/eHdqFiO78JKhZPMTCcG8bDY0m+8pO3+MIP35jy\nIhSDI3E27e70ghDgjQOnOXRyoODXgPP33bynk76h3KHiwVGnoqssHw/jgeGE93fJbkp76vVjPLv1\nBN9eb7E1azvQWCLFk5uP8aUf7ZjwgaYrqzLOHyGYqaOnxj9UHe8ampUmM3Ca/Hqz5m47C+wX7oZx\nfXWUhhqnwo3FU1SXT72FQLQs7I2AeMPUqoxFSkZhPEfiZzHUO1ncZVdCtm3zy1ePALDVGp8Lzpl3\nzgRisSHm/Io9lkgxkNd5nR3ybjXaO0mX7V899Dpf+9ku9mXmmvuH4/zDozv5P9/ZMmWI7zjYw1d/\nuot/euzN3GMbiVNbFaU8GvGGqbOP1x3GTqXTvLyzw7v/V68dBeCK1S0A/PjFQ7x5qIdNuztzXr+7\nb5RI2PkvP1uB43bFL2yqIm3bHC7yQWS6jmRe123COtE9POnz3Dnj+upyGmrKvfuLVsZl4/PJXmVc\nN7eVcSyRYsPOkwp7kUkojOfI96aYL87XPxTnRN6a4+yAfOAXu715ZTfijp4a5GD7+PpSt3N7qjCO\nJ1LEExOHqQeHCy+DGhxNsGVvF//5nzby+IZDOY/1ZIY5O3qcoHgrc6ECgK7MMHAiOfF43PA+mBVc\nadtmaDRJTVUZFVnD1ENZ5yEWT5FIpnhh+0n6huJcurIJGJ9/v9pszfk+2cPGI2NJhkYTXLS0HmBC\nVQ7wvacsfvTcgbO6JKVbGd965RIAjhcIzbN1+JRzbm5Ztxgg5//H6f5Rb5RgIDNnXFcTpT5TGUPh\nZU0u92IRgPd10bIwdVVlcxKWe46c4bNff41vPbGXzz+8LafqFxGF8Zw5kBWU0/G5b27O+QHVn7V2\n9LVdnV4HrFu1/q+HXud41g/oRDLNhp0ncyrGfJPNY8cSaQYKNEYBfP/pffzi1SMAOc1R2Tt3uY1R\n2X/n031jrN90jD/9uw05xwmTV16jsSRp26a2sozyqDOEunlPJ1us3CtSbdjZwcNP76M8Guajd66h\nIupUeJGwwTUXL8h57qmsoV13jnhRcxW1VWX0DuV+AOnocZq/1m8+NuF7ZhsZS/LWoR7vw8+pMyNU\nV0S45IJGAE6enqUwzvQU3LRuEQAnMp35gyNxPv3Pr3L/v7wBQJ87TF0VpaE6qzKuKFYZj7/1FzZV\neX9urq/gzOBYztTDuRqNJfn6z3czMBzn8gub6RkY40s/2lGS7vOzYds2Pf1jWnct54X2pp4DZ1NZ\nZcvejCI/ON2gHhpN8I1f7p7wtclUmm89sXfK1580jOMpBqaojLO3WsxeG93dN36sXZOEcVffKD96\n/gAAT79+nN9+72qqyiMYhpETxvFEimhZ2FvWVFNZ5g3Rf/Wnu7zn1VaVMTiS4PGXnOr8P/zG5bQ1\nVrGktZrDuneEAAAUiElEQVSD7QPUV0cpLwtz21VL2LS7k0gkxOn+MWLxFD98br/3d1zQUEljTTmd\neQ1c2UvBtlpdXJsJ9l2Hz/DG/tN86D2rSCTTfOarrzISSzKStLlqVRNdvaOsWFhLW1MVIcPgZM/M\nwnjnwR4e33CI37ltNauX1nOwvZ+2pioWZz48uGHs7ph2rHOIzt4R+ofiVETDlEfDOZVxoTXGrmjW\nsqe2xqwwrqvgcMcg/UNxGmvLJ/vSs/bEa0fpH47zr29eya/ftIJv/WovL+/s4IU32rnjmmXn/Prp\ntM1Ipj+gpsjw/HTE4ikeffEg2/Z10zsYY3FLNf/+Ny6nLetDi8wftm1ztHOQwx2DHOt0fnWeGaWq\nIkJDTTn11c6oUltjFQubq1jYVEVzXUXB7WTnisJ4Dsz0E/9ff3dr0ecMjsTZebBnwv3DYxMbh1Lp\ndM6Ve4YnCeOxeHLaVVA6c3m+kGF4w9DgBO/IWIKT3cOURUIkkml2HDjtPf7ymx28/GYHt1+zlI/c\nviYnjF/ZdYrB4ThrljUAUFNVNunxLGqqYnCkn+GxJIuaq1i7whmiXrW4noPtA6xa4gw/f+T2NXz4\n1tV87We7eOPAaV7ccZIX3zjpvU5rYyUNteUc6xpiNJakMtPotD/zQSISDrH3WB+2bdPdN+pVoHU1\nUZLJtPdD/6XtJ1jeUkUqbbOwqYqySIgFjZWc7B7Gtm0Mo/Abuf30ML967SgXLKzljmuWkUqn+e6T\ne+kZiPG339/Ov737EsbiKa5f3oBhGCxtrWHP0V6GxxI8v/2E9zqb93RxZmDMC836rDnjqiKVcXnW\nnLHbwAW5HdVnE8anzowwNJJgcUt1zhB5d98oT24+TmNtOb92/XIMw+DDt65iy94ufvHKEW5Zt4iK\n6Mx+DNm2zStvneLRFw56jWzXrV3AfXetpTw6+RrrYoZGE/zdIzs4dHKAmsoyLl7ewN5jffzVt1/n\nk/dcxmUrm2f0umcjlU5zut/ZSS0WT7G4pZraqmjxL5znegdjHGzvp6I8THVFGVUVEVrqKwpefWw6\n9h3v49EXD3LgxHihEA4ZLGisZCye4nDHwKSbM0XCIRY2VbKwaTyg25qqWNBQSU1l2ZTv75lSGM+B\nuVynW6iKPd0/calOR88I1rE+Llxcx7IFNTnreF3DY8lJg3wytu0MnzfURHP2Q+7qHeVAez82cMu6\nRTy3rX3SDwwvbG/nN9+9KmcI/jvrnTW0Fy52rhpUVxXFnuSzzMLmKvZl3lBLWqq9+++5ZSVV5RHe\ne81SwLk4Qnk0TEsmVF7acTLndS5a2uB9UOgdjFFZHsG2bY50DNDaUMGqxfW8truTkz0jvL5nvAHs\nJxsOYWBQW1VGS30Fuw+fYV+m431hs1MxLW6p5tSZEXoHYzRlmqHyxRMp7v/hdvqG4rzy1inWXdjM\n0c5BegZizpKuRIqv/9wZ+TCXOx9Q3DB+9a1TjMZSXLd2Adv2dfPMluOMxJKsW+0EREN2ZVw59Vs7\nkrlMZmNtec4PFnfI+uipQVYvqefZrSd4estxmmrL+e3bLvJ2+sr2842HeXyDs8a6rjrKf/vdq7wq\n8pHnD5BMpfnwrau8DwC1VVHed91yfvryYdZvOsY9t1w45bH2DsbYsOMknb0jhMMhmmrLKYuEeG13\nJ+3dw0QjIa5Y3ULPwBib93RxqmeET/3WO3I+nExHd98of//oTk6eHubGyxby8bsuJhIO8cpbHXx7\nvcWXH3uTz370apa3TTwHs+F03ygbdnawYefJnC1ODQPMZQ1cd0kbN1++yGtCDIojpwZ46vXjvL6n\na0IwNtdVcOd1y3jXusVn9QHrRPcQP37hIDsyP4euWN3ClWtauKCtlsUt1d45dHpVEvQNxjh1ZsT5\n1TNCR+b3yRonK8sjLGisZEFDpfN7YyVtjVW0NlRSW1U243+fKd+xpmmGgK8A64AY8AnLsg5mPX43\n8JdAEvimZVkPzugoAmaqMI6EQ+c0B9VTYOeoZ7acmHDf576x2ftzQ02U92SajCZTXx31qguA5Qtq\nJr2QRU//GE+/ftwbLr1kRSO7j/SyZa8zz3rZhc1stbq91/rzj1zJ9v2neePAabp6R9m+f/L5WHc5\n1OKW6kk/qa5YVMdLOzq857gqohF+/eaVE56/oLEScOZwqysi3H3TSlYvqaemsozWBuexjp4RFrdU\nc7rf2bXqkhVNXHxBI6/t7mTv0V5e39tFNBLitquWsn7zMWxs7rzWGVY93DHIc9ucc76o2TmelYtq\n2bavmwPt/VyXF8a27QylvrC9nb6hOM115fQMxPjZxiPEM01un/3Y1XzlJ2/RecaZ6zaXOfPQSxc4\nr+/unLZuVTPDowl2ZZa5rc6MClSVR7z/X8XmjN0pi6a86tdtinvzUA9XrWnlkecPkErbdPWO8qVH\ndvAXH7ualsz5A/jlq0d4fMNhmusquGRFIxt2dvCFH27nU791Bb2DY2yxulm1uI7rL2nL+T53XruM\nF99o55evHuWqNa2TBlz/UIxvr9/Lhh0dk46WhEMG161dwIdvXU1zfQXJVDqz9vwk//jYm3z63itz\nusansvvIGb76010MjSa489pl/NZtqwllPqTceNkiyssi/NPjb/IPP97JX/7etdRXn32lmrZtYvEU\nY/EUY/Eko7EUo/Ekp3pG2LSn06veKsvD3HBpG3VVUcJhg33H+tib+fXU5uN85I6LzkuFXkwylebQ\nyQFOnRmho2eYMwMxWuorWNJazZKWGpa0VhcMprRts+PAaZ7afNxbxrm4pZp3XtrmTTv0DcXZvq+b\nHzyzn5+9fJjbrlrKe65akrNqIN/p/lF+uuEwr7x1ChtYs6yBD926ynuP5AsZBnVVUeqqohP+D7r7\nvp/qGabjzAhdvaPOr75R2ruHc5Y1ZquIhqmpLKOmsoxIptfnS5+6tcjZLF4Z3wNELcu60TTN64H7\nM/dhmmYZ8EXgGmAE2Gia5s8sy5p8H8a3ka//bOKcruuWdyzi+Sl2iZqpySrRbH1DcZ4t8H2vvKiF\nskiIzXu6WNBQyUfvXENbUxV//tVXJzz3b743PpR+1/XLqYiG2X2kl5ffdILywsV1LG6p9sJ49dJ6\nzOWNrFxUx9d+touNmQaz5W01HOucGPZLW2sm7bRdvmD8jZIdxoUsW1Dj/XnVknovRJ3v7bzWa7tO\n0Ts4xmimG3v1knouzlSjT285TlfvKNeYrdx1w3KOdg7S2lDJ7Vcv42hmy1H3+N03urncCU/reB/X\nrR0Pn/buIb78+HjI1lSW8bmPX8vnH97Gq7tOAdDWWMnS1mo+8YG1vPLWKZa2VHvDxEtbnb9LR2a/\n8tVLGxgeTU4IY8MwaKiJcrp/rGgYu+e4IS+MW+orWdJSzZ6jvfzLc/uJJ9N8/K6LiSdSfP+Z/fy/\n72/nz377HSxorGT9pmP8+MVDNNWV8+cfuZKWBmdY75EXDvK/H3qdRCqNYcDv3H7RhGG9yvII971/\nLV/60Q4e/MVu/svvXEldJuCSqTTPbj3BL1454k1L3HntMi5d0UQqbXO6f4yh0QSXrmzKmSOOhEP8\nm/eZxBNpXt11igd/uYd/+4FLvObHtG2z/3gfJ7qHOTM4RirlTLvsPdrHie4hwiGD3/s1k3dfMfFD\n69VmKx9814U8/tIh7v/hG/zJBy8rOIectm3au4exjvVy9NQgx7uH6O4bZSyWKrj+3wAuXt7AOy9b\nyHUXt02oAs8MjPHL147ywvZ2vvgvO7h0ZRO/dv1yLrmg0Tu3sUSK/Sf62HOkl56BMUZiScZiKSrL\nI9RVl1FbFc1skxoiEnZ+b2yoIhVP0lxfQWtDJdUVkaJDsKf7R3nxjZNs2HFyyubPimiYi5c3cunK\nJpa0VGPbNinb5lTPCM9sPeH1mly2sok7r3P+ffO/9+BInGe3nuC5be38/JUjPPHaUa42W7ntqqVc\ntLQewzBIptIMjiR4cvMxntt2gmTKZmlrNR+6dRWXX9g84yFlwzBorHV2tnOnxVxp26ZvMOaFc1fv\nKKf7RxkaTTA0kmBwNEH76WGSqTRGoe0Q8xQL45uA9QCWZW0yTfOarMfWAgcsy+oHME3zZeBdwKPT\n+6sGV34n9e/esYa2xkqqK8u8H75TuXndoim7omdqYDhOTWUZn773Sr7zlMWipio27OzghksXertt\njcSSXHbhxE/dbU1VXpgA/PpNK7jnlguxjo03eLXUV1BXFaWtsZI9R3uJhEPefM+lK5swDLwAuWJ1\nixdmN12+kNd2dbJuVTONteU585nghO/S1uqc28UszQpjdwjc5Ybx1n3dORuGrFvVTGtDJYuaq7zg\nu3ZtG7VVUf7rvVd6z1u5qJaQAWnbGXFwQ2TFwlqiZSGe39bOrkNnqK0u4+o1C1i/+RgDw3FWLakj\nkUzz++9fS21VlFuvWMIPnt0PwJVrWjEMg1WL61m1OPdT/OKWagzDmSaor47SWl/ButXN/ODZ/ZRH\nwyzJOjf1XhhP/dZ2w3iyeeHLVzWzftMxNu/pYmlrNTdfvohQyGA0nuLxlw7xFw9uoqo8wvBYksba\ncj5975VetXzXDRfQXF/Bt57YS1tjFR+7c82Ev4/3fS5s5tYrl/DC9nb++wOvcc3FC0gk07x5qIfB\nEecCFh+5/SJuvXJJToU1VSOVYRh8/K6LOd0/ypa9XXScHua2q5fS3TvKFqtr0s1eImGDdaua+cCN\nKwpWUAAfeOcF9A/FeG5bO//zodf51zet5MLFddRWldEzMMbJ0yPsO96Hdaw3Z+onmuknqCqPUFEe\nobI8QmU0TEU0QkV5mLrqKO9Y1TLlHH1TXQUfu9Pk3e9YzL88d4Bdh8+w6/AZFjVXURENMxJLTdr5\n7f6/ma6ayjJWLKxlxaJab666trKMwdEE1jHn73agvR/bdpoE33vVUpYvrGFRczVNteV0941y8vQw\nx7uG2HOsjzcOOKNik53zW9Yt4o5rl3kfNidTWxXlnlsu5K7rL+CVXad4btsJNu/pYvOeLiqiYRLJ\ndM5IWnNdBR9810puuGThnDZghQyDproKmuoquDizkuJcFQvjOiB7F4OUaZohy7LSmceyU2cQKPw/\n+W2ksjySs/tTU125F3CxeMobUv70vVfytz/YPuHr77vrYk50DXEkbxjk/TdcwBOZDS6mct0lC8FO\nc6xziFNncnduuvnyRSxdUMNnP3o1tm1zx7XLWNxS7c1FZ8/H/t9/dwPgrGGurojwqS9vBJwg/sCN\nKwC4aFkDKxfVcbhjgNuvduZt33f9cnYf6eV3br/Ie62ayjKuW9vGpt2dTtfz1Uv5xStHSds2/+qd\nK/jt2y7yGn/WrWrm+kvauOv65VRVRKgqL8sZbszu/i2kuqKMS1Y00jsY4+bLF+U8Vl8d5YK2Wq/C\nBSdg3R/y77tuOQ/9ai8t9RWsWzXxg0lZJMxv3W7yw6etnOVUkXCI9127nJ+/csT5tNw3ysF25+3z\nu3es4b2Z8+P9W6xbxP4TfdRXl3N35nxOprwszIduXcWbB3u4dm0bhmHQ1ljFdWsX0FhbntPgsqip\nmuOdQ94HhEJ+89ZV/OCZ/d465mzvyIRxa0MF/+FD67wfanffuILWhgpe2H6S7r5Rrr+kjfffcMGE\n+fHr1rZx+YXNlJeFi/5A/Ogda1jUXMVjLx7yGu3qqqPcee0yfu/uy4iNnP165LJIiP/04Xfw6IsH\neX5bO9990ulLiJaFuHndIi5d0URzXQWRiEE67Sx3qyyyYxk4Qf/RO00uWtrAt9fv9VYL5Guuq+CK\n1S2YyxtZtaSOtsaqWQuG5W21/Nd7r+RwxwBPbj7GVqubUMigsjzCkpZq1q5o5NIVTV4jXTQSYizu\nbOwzNJIgkUyTSKVJZn6vrIrS0TXE6b5RTveP0X56iLcOn+GtAlu7GobTNPnuKxZz7cULJkwDNNVV\neCNE4MyD7zpyhjMDMcIhAyNkUBkNc+3atrMa6i+PhnnPlUu49YrF7Dvex3Pb2jnZM0xFWZhoWdip\nwi9o5NYrlngjIfONMdUyHNM07wdesyzrkczt45ZlLcv8+XLg85Zl/avM7S8CL1uW9djcH7aIiEhw\nFPsIsRF4P4BpmjcAO7Me2wtcZJpmo2maUZwh6omTjCIiIjKlYpWxwXg3NcB9wNVAjWVZD5im+QHg\nczih/g3Lsv55jo9XREQkcKYMYxEREZl783OmW0REJEAUxiIiIiWmMBYRESkxhbGIiEiJnZcLRRTb\n41pmT2bb0s9blvWeUh9LEGW2gf0mcAFQDvwfy7J+XtqjChbTNMPAA8AawAb+yLKsXVN/lcyEaZoL\ngK3Aey3L2lfq4wka0zS3Mb451iHLsv6g0HPP11WbCu5xLbPHNM1PAx8FJm76LLPld4Fuy7I+Zppm\nI/AGoDCeXR8A0pZl3Wya5ruBv0Y/L2Zd5oPl14CZXYRbpmSaZgXAdAuj8zVMnbPHNc7FJWT2HQB+\nA6a5M7nMxCM4a+vBef9M7/qTMm2WZf0U+MPMzRVAb+Fnyzn4AvDPwOxvhC8A7wCqTNN80jTNZzOF\naEHnK4wn3eP6PH3vt43MVqQKhzlkWdawZVlDpmnW4gTzfy/1MQWRZVkp0zQfAv4B+H6JDydwTNP8\nOM4Iz1OZu/QBfvYNA1+wLOt9wB8BD0+Ve+crEAeA7ItFuhebEJl3TNNcBjwHfMeyrB+W+niCyrKs\nj+PMGz9gmmZlkafL2bkPuMM0zeeBK4Bvm6bZVuRr5OzsAx4GsCxrP9ADLCr05PM1Z7wRuBt4ZJI9\nrkXmjcwPrKeAP7Ys6/lSH08Qmab5MWCpZVn/FxgF0plfMkssy3q3++dMIP+hZVmdJTykILoPp2n5\nT0zTXIwzQlxwSuB8hfHjOJ/CNmZu33eevu/blfY4nTufxblU6OdM03Tnju+yLGvihXJlph4FHjJN\n80WgDPiPlmWd/bUURUrrG8C3TNN8KXP7vqlGhLU3tYiISImpiUpERKTEFMYiIiIlpjAWEREpMYWx\niIhIiSmMRURESkxhLCIiUmIKYxERkRL7/1vRnVLgtVbuAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 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", "collapsed": false, "input": [ "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);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFVCAYAAADVDycqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeUXGdh9/HfzBZt1660K8lFsmXZvu6GuDdKQPELAWIH\nOCeEEwLBgdACIYGY8MJLCCam2BgDtrFxwTZgMLgAtuVuy5as3mXpSqu6q7K9l+nvHzP3zp3Zqatd\n6dnV93OOztHsnZ199s7s/T39+mKxmAAAgDn8x7oAAAAgFeEMAIBhCGcAAAxDOAMAYBjCGQAAwxDO\nAAAYpqBwtizrMsuyXs7w9fdblrXKsqzllmXdMPHFAwDg+JM3nC3L+qqkeyTNSPt6maRbJS2W9HZJ\nn7Isa85kFBIAgONJIS3nZkl/K8mX9vWzJTXbtt1n23ZI0uuS3jbB5QMA4LiTN5xt235MUjjDoTpJ\nfZ7HA5JmTlC5AAA4bpUewff2Sar1PK6V1JPrG2KxWMznS2+AAwAwrRUdfEcSztslnWFZVoOkIcW7\ntH+Q6xt8Pp86OgaO4EeiEE1NtZznScY5nnyc48nHOT46mppq8z8pTTHhHJMky7I+IqnGtu17LMv6\nsqRnFe8ev9e27UNFlwAAAKTwHeW7UsWopU0+asOTj3M8+TjHk49zfHQ0NdUW3a3NJiQAABiGcAYA\nwDCEM46pZZsP6b8fWK1QOHKsiwIAxiCccUzd+9Q27Ts8oB0tffmfDADHCcIZZmD5OwC4CGcYgWwG\ngCTCGUYgnAEgiXAGAMAwhDOMwJ7rAJBEOMMIZDMAJBHOAAAYhnAGAMAwhDOMwJgzACQRzgAAGIZw\nhhFoOANAEuEMI/jYhgQAXIQzzEA2A4CLcAYAwDCEM4xAwxkAkghnmIF0BgAX4QwjMCEMAJIIZxiB\npVQAkEQ4AwBgGMIZRojFjnUJAMAchDMAAIYhnGEExpwBIIlwhhHo1gaAJMIZAADDEM4wAt3aAJBE\nOAMAYBjCGQAAwxDOAAAYhnAGAMAwhDOMwFIqAEginAEAMAzhDCOwlAoAkghnGIFubQBIIpwBADAM\n4Qwj0K0NAEmEMwAAhiGcYQTGnAEgiXAGAMAwhDOMwJgzACQRzjAC3doAkEQ4AwBgGMIZRqBbGwCS\nCGcAAAxDOAMAYBjCGQAAwxDOAAAYhnCGEVhKBQBJhDMAAIYhnGEEllIBQBLhDCPQrQ0ASaW5DlqW\n5Zd0h6QLJAUk3WDb9i7P8esl/ZekmKT7bNu+axLLCgDAcSFfy/k6SeW2bV8p6UZJt6Qdv1XSYklX\nSfp3y7JmTnwRcTygWxsAkvKF81WSlkiSbdsrJV2cdjwkqV5SpSSf4i1oAABwBHJ2a0uqk9TveRyx\nLMtv23Y08fgWSWslDUn6g23b/ekvkK6pqXZcBUVxptp5rq+vmnJlnmrlnYo4x5OPc2ymfOHcL8n7\nzrnBbFnWAkmfl3SKpGFJD1uW9SHbtn+f6wU7OgaOoLgoRFNT7ZQ7zz09w+qozPdxNMdUPMdTDed4\n8nGOj47xVIDydWsvk/ReSbIs63JJmzzHKiRFJAUSgd2ueBc3UDTGnAEgKV9T5XFJiy3LWpZ4/AnL\nsj4iqca27Xssy/qlpOWWZY1Kapb0wOQVFdMZS6kAIClnONu2HZP0mbQv7/Ac/5GkH01CuQAAOG6x\nCQmMQLc2ACQRzjAC3doAkEQ4AwBgGMIZRqBbGwCSCGcAAAxDOMMIjDkDQBLhDACAYQhnAAAMQzjD\nCDFuaAYALsIZAADDEM4AABiGcAYAwDCEM8zAkDMAuAhnAAAMQzgDAGAYwhlGoFcbAJIIZwAADEM4\nAwBgGMIZZqBfGwBchDMAAIYhnAEAMAzhDACAYQhnGIG7UgFAEuEMAIBhCGcAAAxDOMMIMXq1AcBF\nOAMAYBjCGQAAwxDOAAAYhnAGAMAwhDMAAIYhnAEAMAzhDCPEWEsFAC7CGQAAwxDOAAAYhnAGAMAw\nhDMAAIYhnAEAMAzhDACAYQhnGIGVVACQRDgDAGAYwhkAAMMQzjACvdoAkEQ4AwBgGMIZAADDEM4A\nABiGcIYZWEsFAC7CGQAAwxDOAAAYhnCGEejUBoAkwhkAAMMQzgAAGIZwhhHo1gaApNJcBy3L8ku6\nQ9IFkgKSbrBte5fn+CWSbpHkk3RA0sds2w5OXnEBAJj+8rWcr5NUbtv2lZJuVDyIJUmWZfkk3S3p\n47ZtXyPpRUkLJ6ugAAAcL/KF81WSlkiSbdsrJV3sOXampC5JX7Ys6xVJ9bZt25NRSAAAjif5wrlO\nUr/ncSTR1S1JjZKulPQTSe+W9C7Lst458UXEcYFBZwBw5RxzVjyYaz2P/bZtRxP/75LU7LSWLcta\nonjL+uVcL9jUVJvrMCbIVDvPM2dWTrkyT7XyTkWc48nHOTZTvnBeJun9kh61LOtySZs8x3ZLqrEs\na1Fiktg1kn6R7wd2dAyMt6woUFNT7ZQ7z319I1OqzFPxHE81nOPJxzk+OsZTAcoXzo9LWmxZ1rLE\n409YlvURSTW2bd9jWdYnJf06MTlsmW3bzxRdAkBSjH5tAHDlDGfbtmOSPpP25R2e4y9LumwSygUA\nwHGLTUgAADAM4Qwz0KsNAC7CGQAAwxDOAAAYhnAGAMAwhDOMwJAzACQRzgAAGIZwBgDAMIQzjBCj\nXxsAXIQzAACGIZwBADAM4QwAgGEIZxiCQWcAcBDOAAAYhnAGAMAwhDOMwFIqAEginAEAMAzhDACA\nYQhnGIFebQBIIpwBADAM4QwAgGEIZwAADEM4wwwMOgOAi3AGAMAwhDMAAIYhnGGEGP3aAOAinAEA\nMAzhDACAYQhnAAAMQzjDDAw5A4CLcAYAwDCEMwAAhiGcYQR6tQEgiXAGAMAwhDMAAIYhnGGEGP3a\nAOAinAEAMAzhDACAYQhnAAAMQzjDEAw6A4CDcAYAwDCEMwAAhiGcYQSWUgFAEuEMAIBhCGcAAAxD\nOAMAYBjCGQAAwxDOAAAYhnAGAMAwhDOMwEoqAEginAEAMAzhDACAYQhnGCHGFmEA4CrNddCyLL+k\nOyRdICkg6QbbtndleN7dkrps2/7apJQSAIDjSL6W83WSym3bvlLSjZJuSX+CZVmflnSemNMDAMCE\nyBfOV0laIkm2ba+UdLH3oGVZV0q6VNLPJfkmo4AAABxv8oVznaR+z+NIoqtblmWdIOmbkj4vghkA\ngAmTc8xZ8WCu9Tz227YdTfz/Q5IaJT0taZ6kKsuyttm2/WCuF2xqqs11GBNkqp3nutrKKVfmqVbe\nqYhzPPk4x2bKF87LJL1f0qOWZV0uaZNzwLbtn0j6iSRZlvWPks7KF8yS1NExMP7SoiBNTbVT7jz3\nD4xMqTJPxXM81XCOJx/n+OgYTwUoXzg/LmmxZVnLEo8/YVnWRyTV2LZ9T9pzmRCGcWMlFQAk5Qxn\n27Zjkj6T9uUdGZ73y4ksFAAAxzM2IQEAwDCEM4wQY1QEAFyEMwAAhiGcAQAwDOEMAIBhCGeYgSFn\nAHARzgAAGIZwBgDAMIQzjECvNgAkEc4AABiGcAYAwDCEMwAAhiGcYQYGnQHARTgDAGAYwhkAAMMQ\nzjACd6UCgCTCGQAAwxDOAAAYhnCGEWL0agOAi3AGAMAwhDMAAIYhnAEAMAzhDACAYQhnAAAMQzgD\nAGAYwhlGYCkVACQRzgAAGIZwBgDAMIQzAACGIZxhBO5KBQBJhDMAAIYhnAEAMAzhDDPQqw0ALsIZ\nAADDEM4AABiGcIYR6NUGgCTCGQAAwxDOAAAYhnAGAMAwhDMAAIYhnAEAMAzhDACAYQhnGCEWYzEV\nADgIZwAADEM4AwBgGMIZRqBTGwCSCGcAAAxDOAMAYBjCGQAAwxDOMAODzgDgIpwBADAM4QwAgGEI\nZxiBXm0ASCKcAQAwTGmug5Zl+SXdIekCSQFJN9i2vctz/COSvigpLGmzpM/atk0jCACAI5Cv5Xyd\npHLbtq+UdKOkW5wDlmVVSvofSe+wbftqSTMlvW+yCgoAwPEiXzhfJWmJJNm2vVLSxZ5jo5KusG17\nNPG4VNLIhJcQxwfuSgUArpzd2pLqJPV7Hkcsy/Lbth1NdF93SJJlWV+QVG3b9gv5fmBTU+24C4vC\nTbXzXFNbMeXKPNXKOxVxjicf59hM+cK5X5L3nfPbth11HiTGpL8v6XRJHyzkB3Z0DBRbRhSpqal2\nyp3nwYHRKVXmqXiOpxrO8eTjHB8d46kA5evWXibpvZJkWdblkjalHf+5pBmSrvd0bwNFo1MbAJLy\ntZwfl7TYsqxlicefSMzQrpG0RtI/SVoq6SXLsiTpx7ZtPzFZhQUA4HiQM5wT48qfSfvyDs//Sya8\nRAAAHOfYhARGYLI2ACQRzgAAGIZwBgDAMIQzAACGIZwBADAM4QwAgGEIZwAADEM4wwgx1lIBgItw\nBgDAMIQzAACGIZwBADAM4QwjMOIMAEmEMwAAhiGcAQAwDOF8HAmFI/ri7a/pj6/vOdZFGYt+bQBw\nEc7HkYOdwxoYDukJE8MZwJT01Bt79dU7l2s0GD7WRZlWCGcAwLj94dXd6uwb1fb9vce6KNMK4Qwj\n0KsNTG27DvQd6yJMK4QzAGBcvNvu7js8cAxLMv0QztNIc2ufnl6xL+txn+8oFgbAtDccSI4zD46E\njmFJpp/SY10ATJzvPrxWkvSeq08TOQxgsvUMBNz/e4MaR46W8zQUjkSPdRGKx12pgCmn1xvOo4Tz\nRCKcpyEf7WYAR8GApyt7eDTMrV8nEOE8DTG2DOBoGPKEczQW02gwcgxLM70QztPEVJ+MQX3bLKPB\nsF7fdEjRKO8MsnOuOw21MyTRtT2RCOdpoqN3xP2/j6YzjtBDz9q67+lteirH7H9gKBHGTfWVicdT\nu5FgEsJ5GsoWzcdjaA+PhnXzw2v15t7uY12UKWV/+6Akae+h/mNcEpjMCeM5iXCm5TxxCOdpKFsI\nZ4vm5tY+DU/TGu/m3V3a0dqnHz6y4VgXZVKEI1E9v6ZF/UPBCX3dyvL4KssRlscgh6ERp+VcIYnl\nVBOJcJ6G/EW8q63tg/ruw2v1v79aN3kFKsBkTfKcksvKivDCmlb95oWd+sWf35zQ1y0tiVflIow5\nj9sfX9+jh5+zp/W4/dBoSKUlPjXUVriPMTHYhOR4kqHp7IxVH+gYOsqFOTqm+8qO1o549/OhrsLf\nv66+UdVVl6mstCTrc5xA8R+HQyETIRqLuXd/u/D0Rp1/2uxjXKLJMTQSUnVlmaor4lFCt/bEoeU8\nDU33QCpGbJrPA3eWrlTMKKyefahrSF+5c7nuenJrzuc5/Q1+P+E8Hn2DwYz/n24GR0KqrihTFeE8\n4QjnaShrOGf4et/w9L1wSJr2a7ScbvvSAscy9iQmeK3f2Znzec5mEjScx6e7f9T9/1Rf5phNOBLV\n0GhYdVVlqqook0Q4TyS6taehmGIZJ395cyoWi2npxoN6cIl9tIp1TEzzbC5aoMhNIuiFGZ+u4yCc\nnR6B+poZbrf2UGB6/q7HAi3n6SjLBdW7td6Oll79cpoH8/HAGRMutPs+GC5sgpyzBSzbMY6Pd/b8\ndA3n3qH4vtoza8pVOYNu7YlGOE9DmS6ngWBE3Z5N6geGs18wYrHYUZ/lXEi4hMJR7T7YX1RgTFa4\nDI+GNGDSkECBv2ahvdTOeQtP0Ezj7v7RaT1rOZ03kIemaTg7LeeZ1TNUUV4iv89HOE8gwnma8GZQ\npkD6zoNrdPvvN+V8jWji++58cqs+9YNXFAjFu0CXbT6kVzccmLjCjtP9z2zTdx5co427ugr+nslq\n+N1w0wv64u2vT86LF8EZE57oXzOUqJyFs7S0X1jTomdWFrZ72LZ9PfqPO5brj8v2TFj5TNDWM6yH\nnrPdvxOvIU9ITdflRX2D8cp+fU25fD6fqipKp+3veiwQztNEvpbngc78S21u+N7LWr+jQ2u2t0tK\n3g7u3qe2GdEFvmJrmySppW2g4O+ZrHW6Tqs5aki3b6E9BIW2hMOR+PNCWXpQfv3CTj368q6CeljW\n2R2SpGdXtxT0s4+14dGwQuH8Y/O3/36TXl53QM+u3D/mmLe1HAhNz7X2vU7LuSa+r3ZVRSmbkEwg\nwnm6mKCM8O6l/LW7V2h/gUEYjcWOWlDlWp+bLjLJ3fOT/fq5PPDM9ryzrtMVOlzhtJhDGVrO3t3k\nvHu6ZxNMBF15qVmXm1A4OqbSOjgS0pd/+rrufCL3UjMpOa7cOxgYc8zp1i4r9bu//3TT2Rd/72fV\nxcO5uqJ03N3am3Z1au9htor1MuuvBeOWOhNb+uWS7Xpt08GiXyd96cx9T28r6Pu+csdyff3uFUX9\nrJTWXhG5XszS2yNtOS/deFD2/p6sx50W5kQbDYb1yvoDWcM0mpht78hUimg0NqZFXWh5nRZzpnDu\n9azbLaSl5GwnW0zdbWNzpz550/Nav6Oj8G8q0q9f2KFv/GKldrb2ul/bsLNTwXBUG5o781Y2y8vi\nlcRMLeP4zll+1VSWKZjo9t51sM/9/3TQ0j6k8jK/e9OL6soyhcLRolcEjAbDuu3RTfr2A2t0sIAe\nvuMF4TxNeC/CQ6MhvbrhoO5/envW57+6MXNwp+/Lvb9tMO/PDoUj6hkIqK0nfysqm6IirojFt0cy\noSkcieqBZ7bre79en/U52bp9j9TvXmrWg8/aevy13RmPp7dQ0nMkEo3qP+96Y0zZC23ph3OEs3dc\nsZD79zqvUcyEsKUbD6q9e1hPT+JdsV7dEP8b2LoneVOUvqFkK7izb3TM93iVJXoCRoNjKyiDIyHV\nVJaqvKxEgVBUW/Z06aYH1+qBJdn/JqeScCSqQ11DOqmx2l0xMLO6XFLqOSzE3kPJ3rkVbx6euEJO\ncYTzNOG9OBdyEfRekLzG84FwtiksVqtny9BiZlUXszHGkXQ7Z7vpQySafM3IJLWcWxJ3hdp9IHNX\nX/pM8fTzNzAcUlf/qHa09KYcK7Qy4YZzhud7JzuNBvKH83hai/nGvCeSN4S9qxi68nTZO+99ps/J\n0EhYNZVlmlHqVzAUUXNrn6TkvImpoLm1T79+fkfGpWBt3cOKRGM6uanG/Vp9Yuy5t8gd0fZ47nx2\nuGt4nKWdfgjnacLbwjmSJSs7EheRTLIFqLfmW4z/d98qz2sX/n0+xX/ff//ZMj2+NHPL0nEk3doj\nWVqF3nPthNjew/36p5tf0s0Pr834/GKXpjmhNKM88/h6rqVwUmpgeNc2F1qZCIUTS6kytZxHvC3n\n/N3ao+MI52jUaW1nPr77YL82FTFrP533PfSuSfZWevKFjHOO0z8nkWhUw4GwqivKVF5eomAoOmnD\nH5Ppvqe36YW1rbo/w9BWS2JPd2841yVazsXeIW3P4eT141D39ArnnoFAUfveexHO08QjL+10/78t\n7d7FEzXOtS7L+N941hKnj+cV8wo+n0+Hu4fVMxDQn5bvzfncQsJ5657ulO0WHUHPRXdjc6fue2qb\norFYxnD+9gNrJCUrN0OjITe4vnb3G/rGL1amfM+uA9krQfHnxMtdVurXroN9Y1ov6RfA9LfA2+3t\nDVNvJSEWiykSjY59L2Ixt8chEo2l9BSkv/a9T23TY3kqSM5So2ImDA4mbkWYLfy/8+Aa3fboxnGv\nY/eGcP+wN5yT5yrTRC9HLBbTSNApY+rfl3N+nJZzNBZzJ85NlRuJtPcM63AiKN/c1zPmM+DcKOfk\nOclwTnZrFxfObd3DKi/z65R5tWrrHp426+H7h4L62s/f0NfvWZn/yRkQztOE965Sf1y6y/3//rYB\n/cstr07Iz7jnT6m3JXxl/QGt2d4+5o9paDSkTbs6Uy6cL68/oFseWe8+d8wfYJHd2oW2RPO1FNt7\nR3TLbzfo678Y+wcU8Myy/fHvN+n1zYfU2j6Y0goKR8aGVywW0xdue0033vWGgqGIuvvj4/HO+fjV\n8zt000NrtdZuz1oup8XaOxDQTQ+u1XceXJNyfCSYPuac+nt6W87ebmhv2UPhqH7wmw3634dSW/uR\naOrCvHB47Pvr9eflexWNxtTannl+gjNBKBCKFBymzs/oz7DRi7dyNN77TXtbxd5AHvBUZHpyhHMw\nFHU/sukVCKciVZ0Yc5ak9sR8jJjGTtIz0U5PD1ogGBkz98QZdjm5qdr9mhPOuSo16WKxmNp7RjSn\nvkonzq5SOBJzZ4FPdZt3dykYjqpknDePIZynIW835rfuXz0prytJDz5r644ntig9Jn/0u4267dFN\nuuW3G9SZaDE89KytrXt71DsY0IGOwTFBX2zL2RvuPQMBbdqVeUlRpuB8fnWLO87Vk2gxZ5phGsww\nC/db96/W719pdh+HI1Hd9ruNKc8ZSYzD9g+HUoLRufC/vumQJGnXweRY22NLd+k7D65xKx3hRLmd\n57SnTbZLL296q9TbmsvWcl69vV07Wnq162C/ejy7x6VXfNKXAg1lWC7z25ea9c37Vmn9zrG9K07L\nORbLPMHM4Q1a5//B0NjZv96Ld3+e7v1svJOW+oeCbmCOBMLuLmq9A9lDxls5Sm85DyVa/dWVZW44\ndw/EP2exWObPlWn2Jrqa33XRyZLiQeOtfB3oGNTMmnLVVpW7X3PGnIu5C1f/cEiBUERzGirdWd8d\neSbiTbbOvpGiZ5xn4gy7/Pc/XTqu7yecp6Fix3zGw1v7T28J7E4Eypt7e3T7H1J3JfP5fLr1dxu1\nentqq3HJyv1FtSi8AfLNe1fqtkc3pSzDCIQi2rCzc0y39uHuYf3mxZ2668ktidfJ/jOzbUTxhmdS\nTzgS1da9qUutegaSF5fNu5Pjoj398Yu9U6bSkvifXzQa05+X79NuT0hm25nrUNeQegYCY3alSv89\nvMedltzTK/alnPd7n0qOJe5oSS4nGvNa6eGTYReo59fENxj5yR82j3kfvWXJtJuWJK3Z3q7P/Wip\n1todisViKUu00u+c5m3pjvez7u16jUSTP28kENacWVXy+3w5x5y9FYlAMKJte7v1pdtfU2v7oAYT\n56emosxd2+0t87HaRatnIODu6pXNqm1tuuvJLdrR0iu/z6f3XLZAPklPvLZH37xvlR59pVnDoyF1\n9QdSxpul+Hpnv8+nQ92Fj7G298S7zuc0VKpxZjycu45hOO862Kev3vmG7n+msCWkQ6OhjL140VhM\nW/d0q3FmhU6YXTWushDOGBdv6OXK1PQWXzgSzXhBHQ1G3K6yfFo7BlOWCDktuf/7i5V6MjFz/M4n\ntuj2P2zSyjdTZ8c64dDRO5ryOJNCWjgZJ0x5WpYPPJNcOpNtkwXvhBGn1ZHpBhWBUERfv2elvnnv\nygzhnNba9Rx3wuD3r+xSNt61vumt2/QJXfk2mkjf2MMb7unl7h8O6ge/Wa87n4hXll5a16rRYCTl\nMzUwFNTQaMgdA/WOv+ebGJdNf+I8N86siD9OtJ5HAmFVV5RqZk15zu7Z9NbyDx7ZoP7hkF5Y2+r2\nVFRXlmlG2dgJfZl6HqR4JXfJyv05N/7p7h/Nue4+m2Aoom//crW+9cDqrO9fR++I7npyq1Zta1dL\n+6AWnlCrWXUVshbUu895cU2r9iQmgM5PC+ey0hKd2FitlvbBgseNnevDnPpK9704Vt3asVhMDySW\nn67a1p73hiXd/aP66p3x+STp71nvQEDDgbAWnlA3ZnlqoQhnjEu22eHprabyspKUyU/hyNgJSJle\ns3cwoFXbMi87eXld9n2+n3x9jw51DbldSukXUW/gDo+Gc45Z5gpuR6Z11NlaRrandRovS/z1D3cn\nL0Z9QwFFo7GM64edEB8aDY85PjAc0lq7Xcs2H1IsFkvZGKOQXZu8laj0oB/bbZv7otXaMajh0ZCW\nbT6kaCyW1nJOvvahriE99uoubdvX4w5r+Hxjx5H7h4O6+eF1+q+7V+hQ11DKRbOQ2eKZOC3n+YkJ\nTf1DwcSs+pgqZ5SqvmaGegcDWXtzci2zc8O5oizjbPts56+1Y0i/e7lZ37p/dcrfSCQa1eBISKFw\nRP9xx3J979frU5YfFWLVtnb1DQbVNxjMujmRt4ImSeefNluS9IUPXqBvfeISveeyBQqGo3ousQ3r\nSZ7xZscp82oUDEX12Vtf1WtZ9lLwcsO5wRvOx6blfKBjKKViuSbHnBBJ+vMb+zQSCKutZ0Q/fGSD\n2jwzzZ09H+bOqhx3eQhnjIt3/am3Ff3J772c8rxoNKabPBOOIpFY1pa23x8fS/7MLa/qyz9dprue\n3Jp1PXYuuVq83hblw8/bKcETi8UUCEZybsCRLlOXVrYwTJ+h7TzP21LoHQym3AvYy+kWl6QX1rSO\nOf6zx7fo3qe2acPOzpTQGhoNjxl7T+e9II4N59TfZ2g07N4i0GtOQ/xC1DsQ1N1/elP3PrVNL61t\nTXm/nVZ0e8+wvn7PSi3deCjlNXw+nxt8FYlg6xsKuhfN3Qf7U8J5Q3OnPv3DV7SxubhtTJ1wdrpm\nB4ZD7s+Nh3O5wpFY1lZuts1XDncNJ7u1K0vd38ErW+XNO6brXe/725ea9a8/fk2bdyf/FgrdVtex\nz/P8bfsyt7ydSV/XXjpfFyyarasvOEFS/HwsmFurCxbFw9oZqknv1paSgR4MR3X/M9tT5jJk0t6b\nDOeGRLf4sQrndYn5Eu+78lRJqX+vfUPBlGtHKBzV8s2H1DizQh9dfKYGR0L67UvJuShtie76uQ3j\n69KWCGeMw2NLd+vWRza4j3PNnE7f3nFDjouo3+fTgc6hlJZWZ9+IwpGobnpoTdbvS5erNeV97RVb\n21LulPTi2lZ95tZX9b1frZNU2BK0TBfpbC2j9G4y5yLd2Zu8GP3q+R36z7veSHne3FnxP/DuPBc6\nR0vHYFoPQcidpJZNd/+o20pMr5SkjzkPjoTcmbleJzXGW1LDgbC2JwJgd1oLzzn/2Xaei0Si7mfm\nxMb4xd+7EuFAx5AGPV3Za+0OhcJRPeK5MKb78/K9+u7Da1Pe+76hgEr8Pp3QWJV4HHR/btWMEtXX\nxic3Pbe6Rf/78NoxLVXnvauvST0PPYMBd0JYTWVZSiXGaUVnC3xn7bCUDC0pWRH707K97tcOdmZe\nD7xuR4cxWVgUAAAS2klEQVS+ee8qd9MThxMWM6vLZbf0Zvyb3d82IJ+kv7l6ob704Qs1q64i5fjC\nE+rcmcd+n08nNo4NnkvPnqt//7u36MJEkG/ZnXstenvPiEr8Ps2qrVCJ369ZdTPGDIXl0z8UTDlf\n4xGLxbR6e7tK/D791SXzVVbqd4fZOntH9JU7lqfs7rbnUL+C4aguPL1Rf/kXJ+nMk2dqQ3On+zlp\nT/SGEc44atp7hvXn5Xu131PLP1TErj651sQuWbV/TBCEIzG9vvmQdmXZKSuTbEG0o6V3zOt4xyyf\neiO+VaQzQzpQQMs507hUpqDw+eK7aXmHAIYS3eovrhvbCvY6wQ3nwloUvYPBlCAaGg3n3AO7uqJU\nwXBU9v5e3fboRnUnWug1lWWSkhWQQCiiH/1uowZHQmqsrxjzOk4lYmQ0rJKS+EXcCSqHO+afZVxx\ncCQ51HDy3Hg47z6YDJolq/Zr696xvSmDw8GMXdC9gwE9tnS3mlv7tMFzk5C+waDqqss1szoewgPD\nweQs64oyd+bxn5fv1c7WPq18s039w0F3ApPzuTlhdmrXbt9g0J2VX11ZltJynpe4UGebxObtFu1I\nBJR3Pba39ZttY4unV+xTa8egHnw2dZvQ9p4R1VaV6S/ObFIgGNG+w6kt71gspv1tg5ozq0oV5WN7\nRaT4ENWpJ9RKin+es92A5txTZ+m6a06TlLoiIZP2nmE11VfKnwj9Exur1T8UzDve6y33Dx/ZoBvv\nemPMJNNi7Gsb0IGOIb3ljEbVVJbp5KZqHewcUjgS1R+X7VU4EtWKrW3uTV+ccX9rfr18Pp+uvWyB\npOReEE5laA7d2jhabvx5cTe3KMbKN9v0ZNo9f3/1/A7tLXJ8Ldt44M2/Wqclq8be3s+R3t1YSMt5\nsMAJSY0zKxRT6hrQ4dGQnk2UJ9NayLec3qiPLj7THdvzdmtncvpJMyVJK988rJfXJ8flh0dDGskx\n7rxgbvyC+/3frNemXV367cvxyoUTUE44r9ne7nZpnjqvTgvmpHZrOt3aw4GQSvzxS0trojVYmgjr\nQDCiLXu6tGV35uGK1o5B3fZofIb/osTvsydtB7pMEweHRsMZZ1d7h0VWvtmmDc2d+v6v16mzb1R1\n1eWqq4pXQPqHQ24Q1laVqyHxuzvae0b008c268afr9Cbe7vdgD0xLZwj0ZgbnNUVZar0BN38RGUj\n04Y3UuosZaf1mKni6/f5UsL5cPewvnrncr2+6ZAbuoe7R5LL8iJRdfaOam5Dlc4+pSF+XtIqOF19\noxoOhMe8p+nef+VCSdIV583L+byTmqpVXup3V25k0j8c1NBo2P3cSMnzmX4DjGwrJ3Yd7Hc/Y0+O\ncxvhWCymJ16Lf+/V58e78ufPqVU4ElNL+6BWeCaVrk2ErzN/5MzEZLmzT2lQid/nfq7bekZUOaNU\ntYkK7ngQzjBKpjHm9HHJfNI36ChU+hKiQu64lG+HMsfsRBfhT/6w2f3a0GhY63Z0yifpJ1+6Zsz3\nnH1qg9510clu12i2i7qjvnaGykv9Y3oOtu7t0X8/kLre/R1vPcn9/4K5qRdkZ31vQ60TzvHz4IzB\nnX1Kg669dL6+/rGL9NnrznO/b25inerQaNit2DhjjvPn1Lq/w62/3Zh13NNr9sxKNzwl6Z/fd07O\n56/Z3q6bHlqTMkPfO264eXeX7npyi7bvj19YZ1aXp2w56bR4a6vK3Ilijm37e9yu4kdebNaBRCCc\n4OnanZfoOTjQMaSZNeUqK/WntJxPSVSCOhPv48PP2frKHcvd1nhX/6jbTe5007albWdZUV4ia0G9\nuvoD7vvy2NLd6uwb1X1Pb3Pnf4QjUXc4oKt/VNFYTHMaKnX2qQ0qLfFp1bb2lJ4Gpycs/bOQ7oJF\ns/X9f7lCH333mTmfV1ri16nzanWgczBrZdlOvA9OJUyKt5wl6WCi8hGORHX77zfp0z98Vbc9unFM\nBXrF1viNMkr8Ph3sHMr4N9IzEMg6f2RgOKgfPrJBm3Z16ZxTG3R+ojv+lHnx9+q51S0KR6LuePuW\n3d0KR6JqPtCnkxqrVZdY511RXqoz59drX9uAOvtGdLhrWCc1Vo97praUJ5wty/JblnWXZVnLLct6\n2bKsRWnH329Z1qrE8RvGXQpgAk3EhJJgKFL02Fc273zrSe5mFN6uyZ6BgFo7BhMX8bFdiVWJUHbC\nOde+55I0o8zvjpU6Mk3akaT3Xr7A/b8TGg6n67mhtjzlcUvHoEr8Pn3pwxequqJMZaUlKT+vsb5S\npSV+7WztHbMU7LyFsyQlJ90UorqyzL1YS1JTfaXOO23WmOfNTbS8fvPiTu060K/fvLDDPbbrYL/K\nSv368DsXKRKNpYzFn9RYrerKMvl9PvUNBVJazvPn1qQsg/IOt7R2DLpr273d2mfOTy45mpOoqFR4\nxpxn11WouqJU3f3x/ZZfWndAXf2jemndAY0EwhoaDevkphrVVJYlw9kzm1mK92Y45+RQ17Ai0ai2\n7kkd1114Qp2k+Lhoc2uf2wqd21Cp6op41/bBzqGUbmBngln6ZyGTxvrKrHu+e5120kzFYskNTdI5\n2wyfk2jNS8kZ4E5v2eNLd7vzVDbt6tK//vg1/ff9q7X7YL9isZg27epS5YwSfegd8WjalDbG/cbW\nw/qPO5bpG/euHFPRCYQibkXx/NNm61MfONfdXtUpk1PRe/tbTlRD7QzZ+3u099CAgqGo22p2OJ/N\nZ1buVzQWc4cAxitfy/k6SeW2bV8p6UZJtzgHLMsqk3SrpMWS3i7pU5ZlzTmi0gATYMnK7F3Xhfp/\n961ya/ZH6h+utXJ2kTckWtXvuWxByterKuIX9qoMM6MzKfH7VVuV2o32fy6bn/G53olKp3taLl6z\nauPl+vPyfbr3qTfV2j6kebOr3FslSlKl5yJdU1mmslLfmJb7hYtm64LT4y2PQucOVM4o0Rnz61Mq\nF7PqZrhd7d7yn5eYIezoHw7puVX7dc+f3lRL+6BOnVerv3zryTq5qVoV5SX66OIzZc2v17svni+/\nz6fG+gq1dY+4a8xrq+KB/T+fvFSfu/48XX7uXPe1P/+358vyhPBCzwX4rFPGhrO3nDWVZZo3q0pt\n3cN6zdMbtHzLYTfA5s2q0pyGSnX2jmhgOOiG5qc/cK4uspp03TUL3QBtbu3TnkMDY8731efHu5wf\nfNbWdx9e6/bWzEmMeV9/zWkqK/XrV8/vUFffiIZHQ+4EvfkFhHOhFp0YryQ0Z9hHfn/bgJZvOaza\nqrKUEDtlbq3qqsu1bken9h7u13OrW9Q4s0I/+7e36bqrF+rExmrtaxvQ3X/cqpb2QXX2jercU2fp\nojObJMV7TxydfSO6/+ntisXiwwR3/+lN91aXG3Z26s4ntmhf24CuPv8EfenDF7it4Pi5qkz5W7Lm\n18taUK/+4ZBeXh+fI3LWgmSlQpLOWxj/HDpLPZ1K0njl+6u/StISSbJte6VlWRd7jp0tqdm27T5J\nsizrdUlvk/T7IyoRYIB896a++XNX68afvV7w63108Zn6xr2rMh773PXxruEPvmORFl8yX1/+6TJJ\n8TFLaeyMYK/Lz5nrjonFYrGUuzj93bvO0KVnz9WSlS3uuJzDG/iN9ZU6d+GslCGFslK/rjhvnns7\n0GWb492H6RtP1HlmbVeUl6QExfw5NaqrLteH3rFIpaWFj6CVl/r1g89cpZk1M1JvrFBT7q6Fra8p\nd7tLzzm1QS+uTZ1U552Ud8Gi2ZpRXqJvfvwShSNRVZSXuttSSvFxzg3NnXphbat8vmT3dGN9pRrr\nK1Nmi8+fU6OPLj5T33loja6/5jRVlJfq+5+5Qn2DwZSWtjOOX+MZc6yuLNUFpzdq18F+LVm1XyV+\nn646f56WbjykH/wmvqnOopNmajgQ1u6D/fri7fHP17xZVVp4Qp0+d/35kpJj04+/ttsdV/7bt52m\nx5buVonfp8vOmaeHnkv2HjicNbdzZ1Xpg287TY+81KyPf/u5lPObaRb+eJ0xv16lJT499cZejQbD\namkfVGNdhUpL/Xpl/QGFIzF94r1nu/MTpPhyykusOXpxXat7M5mPLj5TlTNK9YGrF+oDVy/UIy/u\n1HOrW9ytif/CalJjfaUWnVSnbft6tKOlVw21M/TLJdsVjkT1z+87R5t3d2nFm2363I+WpnRxn7Wg\nXv9wrTWm+9nn8+ltF56op97Yp4+8+wxVVZTp7AUNWrG1TW9sbZPf50uppEnxfcZn1pS7lbxT5x1Z\nRSdfONdJ8lZ1I5Zl+W3bjiaOeatEA5IyV8Ex6RZfPN/dQnE85s2qcndgMlFNZdmYGZyz6yr09Y9d\npN0H+/XWMxr17KoW/e7lzEtq3v6WE1VVUarFF89Xfc0MbdvbrVt/t1GRaEwff89Z2tHSq+Vbxt7o\nfU5DZcbu7Ya6GXr/lae6Y84ff89ZeuCZ7e55/KtL5uu51S166xmNkqSTmmr0uevP088e36I59ZU6\nZ+EsvbL+gK46f57bOvT7fG7L0Pn9JLl7Dkvx2rizXOPsUxr0qQ+cq+6BgHa09GpmTblm11Voz6F+\nvfuik/VXl8Rbzd/+5KVav6NDz65u0Y6WXi2YUyOfz6fPXX++O4b39gtPTAnnr/zdW90g9FqQ1rLy\n7q3s8/n02evO08//uFXzZlfpX/7mXLfb17u87e/ffYb+8OpuXXn+PNVWlqmlfVDrd3aqoXaG6qrL\n9ZF3neH2GjjjkXXV5Srx+3XtpQvUOxjUlefO048e3aiRQFhnLWjQBYtma9OuLl139ULtOtivEr9P\nV543TyOBsDt5qbTE726Z6nX2KQ1u1+kZJ9ePWcN97mmz3EqK81789EtvcyfxNc6Mbz3pveg7u2rV\nVJbpH661tG1fj+Y0VOqaC07QU8v3KhiO6uKz5ujD7zxddkuf2rqHVVbq15nz61VW6k/5LL7l9MaU\n8syeWZFSKauaUaprL52vcxfOUlN9paoqSt3P36IT69zz4R1HX3zJfA0Hwtq0u1v9QwF19wd0xTm5\nJ3kVq66qXJedPVfLthzWMytSe7NmVpfrY9daemuixev1vqtO1cGuIR3sHNJfX3GKLkz7/f/m6oVa\nvb1dPQMBNc6s0CVnxTtsF188X3c9uVU3J5ZCStK5C2fpsnPn6i+sJoXCUW3Z063zFs7S2ac0aN7s\nKl14emPWO4Vdf81pWnzJfLdFffm5c/Xbl5o1HAjrLy86KaViKsU//1ecM09LVu1XQ+0Md/XCePly\n7WdsWdYtklbYtv1o4nGLbdvzE/8/X9LNtm3/deLxrZJet237sSMqEQAAx7l8fU3LJL1XkizLulyS\n9y4G2yWdYVlWg2VZ5Yp3ab8x9iUAAEAx8rWcfZLukHRB4kufkHSRpBrbtu+xLOt9kr6peMjfa9v2\nnZNcXgAApr2c4QwAAI4+NiEBAMAwhDMAAIYhnAEAMAzhDACAYQrbF7BIlmX5lZzlHZB0g23buzzH\n3y/pG5LCku6zbfsXk1GO6ayAc/wRSV9U/BxvlvRZ27aZ/VeEfOfY87y7JXXZtv21o1zEaaGAz/Il\nim8d7JN0QNLHbNvOfN9FZFTAOb5e0n9Jiil+Tb7rmBR0GrAs6zLF9wB5Z9rXi8q9yWo5syf35Mt1\njisl/Y+kd9i2fbXiO7e975iUcmrLeo4dlmV9WtJ5il/UMD65Pss+SXdL+rht29dIelHSwmNSyqkt\n32fZuSZfJenfLctit8dxsCzrq5LukTQj7etF595khXPKntySMu7Jbdt2SJKzJzeKk+scj0q6wrZt\n5/ZMpZIm5hZLx5dc51iWZV0p6VJJP1e8VYfxyXWez5TUJenLlmW9Iqnetm37qJdw6sv5WZYUklQv\nqVLxzzKVzfFplvS3Gns9KDr3JiucM+7J7TnGntxHLus5tm07Ztt2hyRZlvUFSdW2bb9wDMo41WU9\nx5ZlnaD4BjyfF8F8pHJdLxolXSnpJ5LeLeldlmW9UyhWrnMsxVvSayVtkfQn27YLu30YUiS2r850\nA+uic2+ywrlfkneHfOdmGVK8gN5jtZLy33kd6XKdY+de3D+U9C5JHzzahZsmcp3jDykeHE9L+k9J\nf29Z1seOcvmmi1znuUvxFodt23ZY8dZfeqsP+WU9x5ZlLVC8knmKpFMlzbUs60NHvYTTW9G5N1nh\nzJ7cky/XOZbiXa0zJF3v6d5GcbKeY9u2f2Lb9sWJSR83S/q1bdsPHptiTnm5Psu7JdVYlrUo8fga\nxVt3KE6uc1whKSIpkAjsdsW7uDFxis69Sdm+kz25J1+ucyxpTeLfUs+3/Ni27SeOaiGnuHyfY8/z\n/lGSZdv2fx39Uk59BVwvnAqQT9Iy27b/7diUdOoq4Bz/m6S/V3y+SrOkf070VKBIlmWdqnhl/crE\nqplx5R57awMAYBg2IQEAwDCEMwAAhiGcAQAwDOEMAIBhCGcAAAxDOAMAYBjCGQAAw/x//9/njRa5\nNGcAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 6 }, { "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", "collapsed": false, "input": [ "from gatspy.periodic import LombScargleFast\n", "help(LombScargleFast.periodogram_auto)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Help on function periodogram_auto in module gatspy.periodic.modeler:\n", "\n", "periodogram_auto(self, oversampling=5, nyquist_factor=3)\n", " Compute the score on an automatically-determined grid\n", " \n", " Warning: depending on the data window function, the model may be\n", " sensitive to periodicity at higher frequencies than this function\n", " returns! The final number of frequencies will be\n", " Nf = int(0.5 * oversampling * nyquist_factor * len(self.t))\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" ] } ], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, "input": [ "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);" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAFVCAYAAADc5IdQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHHd95/F3dff03PeM7tOSXZIP+T7wAcZgg1m8gQTY\nOIEFb9jNskk2m+yGENjlybPZLMkSQxISAhhiHHBsMGCMwcjGtmzZsizZOq2rdN/X3GdPT3dX7R/V\nXX3M9PRImlHXlD+v59Ez091z1JRm+tPf3+/7+5XhOA4iIiJSPqFyH4CIiMjbncJYRESkzBTGIiIi\nZaYwFhERKTOFsYiISJkpjEVERMpsUmFsmubNpmmuGef++0zT3Gia5mumaX566g9PREQk+EqGsWma\nnwUeAioL7q8AvgLcDbwL+E+mac6ajoMUEREJsslUxvuBXweMgvtXAvsty+qzLCsBvAq8c4qPT0RE\nJPBKhrFlWT8BkuM81AD05dweABqn6LhERETeNiIX8Ll9QH3O7XqgZ6JPcBzHMYzCAltERCTQSgbf\nhYTxHuBS0zSbgSHcIeovT3g0hkFHx8AFfEuZjPb2ep3naaZzPP10jqefzvHF0d5eX/JjziWMHQDT\nNO8H6izLesg0zT8GnsUd7v6OZVmnzudARURE3s6Mi3zVJkevwqafXu1OP53j6adzPP10ji+O9vb6\nksPU2vRDRESkzBTGIiIiZaYwFnkbuMjTUSJyjhTGIgE3GEvwO3+9hqdePVTuQxGRIhTGIgG373gv\ngMJYxMcUxiIBl0ppiFrE7xTGIgFna75YxPcUxiIBp8pYxP8UxiIBl7Ttch+CiJSgMBYJONtWZSzi\ndwpjkYBTGIv4n8JYJOCSCmMR31MYiwRcMqk5YxG/UxiLBNyowljE9xTGIgE3mkyV+xBEpASFsUjA\njSZUGYv4ncJYJOAS6co4HCp5fXMRKROFsUjAZeaMFcYi/qUwFgm4zDB1OKwwFvErhbFIwGUauEKG\nwljErxTGIgGXSFfGIQ1Ti/iWwlgk4DKVsbbFFPEvhbFIwGUauHRdYxH/UhiLBFxmmFrXNRbxL4Wx\nSMDF08PUKQ1Ti/iWwlgk4LzK2HZwNFQt4ksKY5GAy71QhOaNRfxJYSwSYLbtkExlw1jzxiL+pDAW\nCbBEweUTNW8s4k8KY5EAK7x8osJYxJ8UxiIBVnj5RIWxiD8pjEUCrLAy1i5cIv6kMBYJsDFzxim7\nyEeKSDkpjEUCbMwwtZY2ifiSwlgkwMY0cGlpk4gvKYxFAqywYUsNXCL+pDAWCbDChi01cIn4k8JY\nJMAy21+GDAOApK0GLhE/UhiLBFimEq6ocP/UNWcs4k8KY5EAy4xKV4TdP3UNU4v4k8JYJMC8yjiS\nrowVxiK+pDAWCbDMnHEk7M4ZpzRnLOJLCmORAMtWxmFAlbGIXymMRQJszDC1GrhEfElhLBJgmWHq\nTAOXKmMRf1IYiwSY100dUTe1iJ8pjEUCrHCYWpt+iPiTwlgkwAqHqVUZi/iTwlgkwLTOWGRmUBiL\nBJi3zljd1CK+pjAWCTCvMlY3tYivKYxFAqywm1o7cIn4k8JYJMAczRmLzAgKY5EAS9nqphaZCSIT\nPWiaZgj4OrAKiAOftizrQM7jHwY+DzjAP1uW9Y1pPFYROUfe0iZVxiK+Vqoy/hAQtSzrVuBzwIMF\nj38FuBu4Dfjvpmk2Tv0hisj5yl61Sd3UIn5WKoxvA1YDWJa1Abih4PEE0ARUAwZuhSwiPuGk+7VU\nGYv424TD1EAD0J9zO2WaZsiyrExL5oPAJmAI+LFlWf2FX6BQe3v9eR2onBud5+k3E85xZVUFAC3N\nNQBEKyMz4rgzZtKxzlQ6x/5QKoz7gdz/KS+ITdNcBPw+sBgYBr5vmuZHLMv60URfsKNj4AIOVyaj\nvb1e53mazZRzPDgUB2AkNgrA0FB8Rhw3zJxzPJPpHF8ck3nBU2qYeh3wAQDTNG8Btuc8VgWkgHg6\noM/iDlmLiE8U7k2d1DC1iC+VqoyfBO42TXNd+vYDpmneD9RZlvWQaZqPAK+ZpjkC7Ae+O32HKiLn\nasw6YzVwifjShGFsWZYDfKbg7r05j38V+Oo0HJeITIHCburMbRHxF236IRJgdmE3dUrbYYr4kcJY\nJMBS2g5TZEZQGIsEmKMduERmBIWxSIAVdlMrjEX8SWEsEmCZC0NENGcs4msKY5EAyxTC2W7qMh6M\niBSlMBYJsExlHDIMDCM7hywi/qIwFgmwzJxxKOQGstYZi/iTwlgkwPIrYwNlsYg/KYxFAixbGRuE\njGw4i4i/KIxFAsyrjEMGRkiVsYhfKYxFAix3mDpkaG9qEb9SGIsEmO24QQzuW3VTi/iTwlgkwGzH\nIZT+KzcMQ+uMRXxKYSwSYLbt5FTGWmcs4lcKY5EAcytjN4wNw1A3tYhPKYxFAiyvMlY3tYhvKYxF\nAsx2yKmM1U0t4lcKY5EAcytj9311U4v4l8JYJMBsx8HIq4zLfEAiMi6FsUiA2bZDOJRdZ6wGLhF/\nUhiLBJjtZBu4DA1Ti/iWwlgkwAq7qVUYi/iTwlgkwGwHb85Ym36I+JfCWCTAcrup3e0wFcYifqQw\nFgkwxylo4FIWi/iSwlgkwFJ2bgMXOEpjEV9SGIsEWO46YzVwifiXwlgkwGyb/MpYc8YivqQwFgkw\nJ+d6xiE1cIn4lsJYJMBs2yGct+mHqmMRP1IYiwSU7Tg4ZK/alFnipCgW8R+FsUhAZfahNnIq49z7\nRcQ/FMYiAZUZjg7ldFO795ftkESkCIWxSEDZtvs2t5saUBOXiA8pjEUCKhO6uTtwgRq4RPxIYSwS\nUClvzti9HfLmjMt1RCJSjMJYJKDsgjljDVOL+JfCWCSgMvtQe9cz1jC1iG8pjEUCKrOCyauM02+1\nsknEfxTGIgFlj6mM3ftVGYv4j8JYJKBS3pyxezukTT9EfEthLBJQhXPGhqFNP0T8SmEsElCF3dQh\ndVOL+JbCWCSgCueMjZC6qUX8SmEsElCF3dTZyrhMByQiRSmMRQJqbDe1KmMRv1IYiwSUXdBNrUso\niviXwlgkoIpXxmU7JBEpQmEsElDam1pk5lAYiwTUmMo4pMpYxK8UxiIB5YWxKmMR31MYiwSUt7Sp\n8HrGCmMR34lM9KBpmiHg68AqIA582rKsAzmP3wg8CBjACeDfW5Y1On2HKyKTNXbOOD1MbZftkESk\niFKV8YeAqGVZtwKfww1eAEzTNIBvAZ+yLOsO4AVg6XQdqIicm7Fzxun7VRmL+E6pML4NWA1gWdYG\n4Iacxy4DuoA/Nk3zJaDJsixrOg5SRM7dmMoYbfoh4lelwrgB6M+5nUoPXQO0AbcCXwPeC7zHNM13\nT/0hisj5KNZNrU0/RPxnwjlj3CCuz7kdsiwrM+PUBezPVMOmaa7GrZzXTPQF29vrJ3pYpojO8/Tz\n+zmuO+m+jm5oqKK9vZ66ukr3dmO17489Y6Yc50ymc+wPpcJ4HXAf8IRpmrcA23MeOwjUmaa5LN3U\ndQfw7VLfsKNj4HyPVSapvb1e53mazYRz3NsbA2B4KE5HxwCxYbe3sqdn2PfHDjPjHM90OscXx2Re\n8JQK4yeBu03TXJe+/YBpmvcDdZZlPWSa5u8A/5pu5lpnWdYvL+iIRWTKZOaMM5dO9IapNUot4jsT\nhrFlWQ7wmYK79+Y8vga4eRqOS0QuUGZuOGxo0w8Rv9OmHyIBVbgDly6hKOJfCmORgPKWNhn5m37Y\n2vRDxHcUxiIBlZkbNtJ/5ZltMVUZi/iPwlgkoArXGRvam1rEtxTGIgHlOLqEoshMoTAWCShvmFrd\n1CK+pzAWCahsN7V7W93UIv6lMBYJqMJhaq8yVje1iO8ojEUCqvCqTaqMRfxLYSwSUGOu2qRuahHf\nUhiLBFSmgSuzvjiz3lhZLOI/CmORgCq2HaYqYxH/URiLBJR31abCYWpdtknEdxTGIgFV2MBleNth\nluuIRKQYhbFIQDneph/uWw1Ti/iXwlgkoMbsTa3tMEV8S2EsElCFl1AMaTtMEd9SGIsElJPeaUub\nfoj4n8JYJKCylbF721A3tYhvKYxFAmrsdpju/SqMRfxHYSwSUGMauNRNLeJbCmORgPI2/chUxqFM\nGJftkESkCIWxSEBlLpWYnTN236qBS8R/FMYiAVV4PWNt+iHiXwpjkYAq3Js68zaz5ElE/ENhLBJQ\n3iUUvTnjzP2qjEX8RmEsElCOXbDOGA1Ti/iVwlgkoMasM9be1CK+pTAWCaix64zT9yuNRXxHYSwS\nUN6ccUE3taOFxiK+ozAWCajsMLV7O1sZl+mARKQohbFIQGWGqTNLmrI7cCmNRfxGYSwSUM6YC0Xo\nEooifqUwFgmozHB0enQ6O0ytTT9EfEdhLBJQtuNgGDnD1JnKGFXGIn6jMBYJKMd2vACGnEsoqoNL\nxHcUxiIBZTuON18M2vRDxM8UxiIBZdsUVMbp+5XGIr6jMBYJKLcyzt7OdlOX6YBEpCiFsUhA2U7+\nnHHI66ZWGov4jcJYJKBs2/GatiDnesYqjUV8R2EsElC2Q34DV6abWlks4jsKY5GAcpc2ZW9n5o9V\nGYv4j8JYJKDcTT/GWWesMBbxHYWxSEA5Yxq41E0t4lcKY5GAcueMs7cNdVOL+JbCWCSg7HG2wzTQ\nnLGIHymMRQKqcDtMcANZhbGI/yiMRQKqsDIGd9haDVwi/qMwFgko2yGvmxrc2xqmFvEfhbFIQBXu\nTQ1uR7Vtl+d4RKQ4hbFIQBVezxjcYWpVxiL+ozAWCajC7TABDAzNGYv4UGSiB03TDAFfB1YBceDT\nlmUdGOfjvgV0WZb1Z9NylCJyzhzHoaAwJhQytOmHiA+Vqow/BEQty7oV+BzwYOEHmKb5u8CVgP7E\nRXxkvG5qw1A3tYgflQrj24DVAJZlbQBuyH3QNM1bgZuAbwLGmM8WkbJwHAcHxs4Za52xiC+VCuMG\noD/ndio9dI1pmnOBLwK/j4JYxFcy1e/YTT/cxi4R8ZcJ54xxg7g+53bIsqzMwoiPAG3AM8AcoMY0\nzd2WZf3LRF+wvb1+oodliug8Tz8/n+NEMgVAZWUk7zgjkTBGyPD1seeaKcc5k+kc+0OpMF4H3Ac8\nYZrmLcD2zAOWZX0N+BqAaZqfBFaUCmKAjo6B8z9amZT29nqd52nm93McT7hhnEym8o7TsR2Stu3r\nY8/w+zkOAp3ji2MyL3hKhfGTwN2maa5L337ANM37gTrLsh4q+FiNfYn4RObKTOOtM04XzSLiIxOG\nsWVZDvCZgrv3jvNxj0zlQYnIhcls7DG2m9rAdrQFl4jfaNMPkQDK9GgVNnCFQ4YauER8SGEsEkCZ\nbuoxm34YBimFsYjvKIxFAsgpMmes6xmL+JPCWCSAJhqmtpXGIr6jMBYJoGw3df79oZC2wxTxI4Wx\nSADZRbqpQ6qMRXxJYSwSQF4DV2icvakVxiK+ozAWCaCim34YBg7Zdcgi4g8KY5EAKtbAlbmteWMR\nf1EYiwSQU7SBKx3GGqoW8RWFsUgAFWvgCnthfNEPSUQmoDAWCaDMKLQxzpwxoF24RHxGYSwSQF5l\nXPAXnslmzRmL+IvCWCSAinVTh9XAJeJLCmORAMpWxkW6qTVMLeIrCmORAMqEbbE5Y4WxiL8ojEUC\nyFtnrKVNIjOCwlgkgIoOUxuaMxbxI4WxSAAVu55xdgeui35IIjIBhbFIAJVq4NI6YxF/URiLBJDt\nbfqRf38mmx2FsYivKIxFAqj0MLXCWMRPFMYiAVRsb2pthyniTwpjkQAqdglF7cAl4k8KY5EAsnUJ\nRZEZRWEsEkCZytcoSGNDO3CJ+JLCWCSASl8o4qIfkohMQGEsEkDFG7jSjyuNRXxFYSwSQI7XwJV/\nv5Y2ifiTwlgkgIpWxtqBS8SXFMYiAeSUuISiduAS8ReFsUgAFVtnrGFqEX9SGIsEUNF1xtqBS8SX\nFMYiAVRqzliVsYi/KIxFAqjYph9h7cAl4ksKY5EAKrbph+GtM77YRyQiE1EYiwSQXWSdsS4UIeJP\nCmORACp6PWPtTS3iSwpjkQDKVL5hLW0SmREUxiIBlFm6NGadsSpjEV9SGIsEUNEwVje1iC8pjEUC\nKBO2kYIOrpAuoSjiSwpjkQAqNUyd0tomEV9RGIsEUPFhavetKmMRf1EYiwSQna58I4U7cOmqTSK+\npDAWCaBUauIGLl0oQsRfFMYiAZQqss44c31jrTMW8ReFsUgAZbqpC8NYF4oQ8SeFsUgAlRqmVmUs\n4i8KY5EAShWpjLM7cF30QxKRCSiMRQIouzd1kU0/NEwt4isKY5EASqXc0nfsph/uWw1Ti/iLwlgk\ngIoOU2tpk4gvKYxFAsi2HQzUwCUyU0QmetA0zRDwdWAVEAc+bVnWgZzH7wf+EEgCbwH/xbIs/ZWL\nlFnKdsYEMWgHLhG/KlUZfwiIWpZ1K/A54MHMA6ZpVgN/AdxpWdbtQCPwwek6UBGZvJTtEA6PDWNv\nmFqVsYivlArj24DVAJZlbQBuyHlsBHiHZVkj6dsRIDblRygi58y2nTHzxZCzA5cqYxFfmXCYGmgA\n+nNup0zTDFmWZaeHozsATNP8A6DWsqznS33D9vb68z5YmTyd5+nn63McMoiEQ2OOsaIq6r6NRvx9\n/Gkz4RhnOp1jfygVxv1A7v9UyLIsb7uA9Jzy/wOWA78xmW/Y0TFwrsco56i9vV7neZr5/RyPjqYw\nDGPMMQ7GEgDEYglfHz/4/xwHgc7xxTGZFzylhqnXAR8AME3zFmB7wePfBCqBD+cMV4tImRUbpvbW\nGWuYWsRXSlXGTwJ3m6a5Ln37gXQHdR3wJvAfgLXAi6ZpAvydZVk/na6DFZHJSdm2t/VlLi1tEvGn\nCcM4PS/8mYK79+a8H57yIxKRC5ayHaIVYwe+QmrgEvElbfohEkCpYsPUqoxFfElhLBJAdpFNP3Sh\nCBF/UhiLBFDRyljD1CK+pDAWCaBiYQzuxSO0A5eIvyiMRQKo2DA1uLtw2fa4D4lImSiMRQLGcZx0\nZTz+n3c4ZKiBS8RnFMYiAZMJ2mLD1KGQ5oxF/EZhLBIwmaAtNkwdMlQZi/iNwlgkYFJ2qcrYUGUs\n4jMKY5GAKRnGhsJYxG8UxiIBkyo1TK0GLhHfURiLBIytylhkxlEYiwRMKjWJbmplsYivKIxFAiaz\nu1bxYeqQKmMRn1EYiwRMdph6/D/vkJGdVxYRf1AYl1FHb4y/eXwLJzuHyn0oEiCplLvX5URLmxw1\ncIn4isK4jB5/YR+7Dvfw8DO7y30oEiClljaFDUOVsYjPKIzLKPOEmNQTo0whu8ScsaGlTSK+ozD2\nAz0vyhQq1U0dDumqTSJ+ozAuo8xTpaM0linkDVOHtc5YZKZQGJeRYaSfLPW8KFPI24HLKBbG7lC2\nmrhE/ENh7AN6SpSpVHIHrvT9ymIR/1AYl5FXGOtJUaZQdpi6yDrjdBiriUvEPxTG0yyZsjU/JxdV\nKt2dVXSYOhPG+r0U8Q2F8TT7va+u5U+/sb7ER+lJUabOZC4UAdqFS8RPFMbTLJG06eofAeCVbSf5\n/a+upW9oFMg2cOkpUaZSqW7qTEgrjEX8Q2F8ET38yz0Mx5Ns298JZOeMlcYylUpdz7gi4v7ZJ5Ja\nbCziFwrjMqisCAO564xFpo43TF1kzjia/v0bTaYu2jGJyMQUxhfJD1/c773vPUdmhqnV1SpTqNQw\ndWUkHcYJVcYifqEwnka5Ibt649Exj4//VClyYUoNU0cr3D97VcYi/qEwnkalGmSKjCKKXJBS1zP2\nhqlVGYv4hsJ4Gk22W1Wj1DKVSl1CMZpu4BpNuJXxyc4h/vzhjTz09M6Lc4AiMobCeBplrp5TirJY\nppK36UfRYepMA5dNLJ7kS9/fxNEzg6zfeeaiHaOI5FMYT6NUievUGVOwH+Ymq4OzvbHz/nwJnlKb\nfuRWxh29MYZGkt5j8VHNI4uUg8J4GpWcM06/Pd8oPtsb4x+ffIsvfOv18/wK+fqHR3nw8S0cONk3\nJV/vfGw/0MU//uQtrYG9ACWHqXMq48FYIu+xzvQGNSJycSmMp1GxYWqvEL7ANO4diLvfZ4p2Uvrx\nSwfYebiHh5/ZMyVf71w5jsPfPrGNTXs7OHp2oCzHMFPtO97L9gPuZjKZ37vJVMaZMG6sjQLQ1acw\nFikHhfE0KjlMnX7rnGcaZxpwpkp3uirKbEpysY3kDJHGcoZOZWKO4/Cl72/mb5/YjuM43tWYSs4Z\nJ1IMDLthvGROPYC3dauIXFwK42lUrGLNhu+FrW2arr2FizyHT7vcIdPC4VMpLrdnYCCWyBmmLra0\nKbPOODtMvTgTxqqMRcoiUu4DCLKiYVxQMJ9v/9aUh3GZL1yhMD4/e470eO93949MYpg6u854JO6O\nRiyZ0wCoMhYpF1XG06jYnHH/sHvVpvU7T1/Q17+QJqc395zlT77+mjc07QeZIVNQGJ+L3Tlh3NUX\nL91NnbMD10DM/V1cOKuOkGGoMhYpE4XxNCpWuf7gxf30Dca92+NVxrbt8Iv1hzneMVj0619IGH/z\nZzvp6h/huTeOefd5c9hlKo0H08EAMBTTnPFkHTrV773fPTBCstQ648jYOeOG2ijN9VFVxiJlojCe\nRhM1cJ3qGs65NTb93rTO8uOXD/LDNfvHPJaRSI39+smUzQ9+ZXGic2jCY8tUT5nL6eUfRXnSeDCn\nMh4ZzYbxK9tP8uXHtjA8BU1dwyMJOgK0LttxHHoHsy9iuvtHvPNUUzX+LJRXGSfcOeOqaJiKSIjW\nhip6B+Ikx/m9EpHppTCeRhPtwGU7E0dfpmLZcbC76NcYrzJ+eetJvr96Dw8/s3vcz7GO9vCX33vT\n+565VXAq/SRcrsp4IGdoeiSnU/yx5/ex+0gPv9xw5IK+/tBIgj/4u1f484c3BiZwYvEUiaSdbcDq\njzM84p7H2qJh7FbG8aS7tKmuugKA1sYqHKB7ID7u54nI9FEYT6N1O04VfSyeEzbnG36Jca6605se\n/j58avx1uj9bd5gDJ7LDmvHRFMMjCb782Bb2HO0FKFtQ5c4T5+4ElTk9vROExHMbj/Jn31zPjkNd\nRT9m675OHMcNsKkMHNt2ON4xOO7/x3TrG3J/joWz6giHDLr7RxgcSRKNhKiIjL9ErSKSrYwHhhPU\n12TDGNRRLVIOCuNpMjySYN1bxRu0im07aNsO8dHUpK7oNF5l7O2wWWSoeaigMSqeSLH3WF9eE1C5\nruaTN0yd92LF/VkGijR1OY7DD9bs50xPjDWbTxT9+md6ssPTUzlU/aXvb+KL39nIT189NGVfE2D3\n4W7+4pE3+d5zVtGP6UsPUTfXVdLSUElX/wjDI4miQ9QAIcOgIhJiYHiUZMqmrtrd8KO1QWEsUi4K\n42kyHJ94fvPZnMapXI/+ai+f+crL9BRUbrbtsOtwtzfXC+PPGWfasIpV24WfE0+kGBopCOhkipe3\nnmCT1THhzzDVBmIJDKAyGvZerKRsm0T6xcHAcHZudDCW4NHn9vLSlhOc7Bzyft6dh7uLVva5Adw5\nRWHcPzzKgZPuSMOuQz0lPvrcfPsXuzl0qp81m0/k/b/n6htyz0ljXZSW+ir6BkfpH05Qmx56LiYa\nCdHd7/6O5Q5Tg5Y3iZSDwnialOp0PnI6O4zs5CTnmi1uZXfkTP4w82PP7+NvHt/Khl3ZK+skxqlg\nS+36FSt4kRBPpLylVhlDsQSPrLb4xyffyju2Qt39I6zZfDwvKNa9dYpdh4vPc4P7867fcXrMsqrB\nmBsiNZURL4z7hxJejZ+79Om5N47ywubj/MuzFpv3Zl80jCZsfrnhKF/5wdYxP+vZvMp4agIndzrg\n2NnBvMazC9EzEM97QXamZ3jcj/PCuDZKY51b4cZHU9RWTryFQLQi7E2VeMPUqoxFykZhPE3OZZ3s\neHGXzAnz+GiKN/a4IbznaLb6yq1yM6E5Olo8jB3HYbBgydBoIsXAUP6xJnMaz/rTT/anuobGbL/5\nxe9s5HvP7WV3+ph6BuJ85xe7+ZvHt+Y1qBXauq+Th36+i3/4yVt59w8Oj1JXXUFlRdgbps6thjNh\nnLJtXt2enY9/ZsNRAK69tA2AJ9ceZMeh7rwXLuBWxpn50qmq/g6fdqviua012I7DoZP9JT5jkl83\nvVwp04R1vGP87vjMnHFjbSVNdZXe/SUr45wtT73KuGF6K+N4IsUr204q7EXGoTCeJoVBM5G+wVFe\nL9gApD+nCnzwB1u99zOXXdyyr4M395z17s8EaHyC/apHE/aYIdzxKuNcA7EEr+88zRce2sCPXz6Y\n91hmKP50epnWjoPZ5qkz3cM4jjOmOgXYe9xtFDucMzpgp18o1NVU5A1T585xxxMpRhMpXtpykt7B\nUa5Y2uLen/7Y6832vO+z81C2Qh8eSTIYS7B8fiNA3jrvjH951uLxF/ZNOBpQKDPCcec18wG3Op4K\nh9Ihf8eqeQAcz/m6nX0x70VKZs64oS5bGUPxZU0Z0ZwlbZnPi1aEaaipmJaw3H24m89/63Ue/uUe\n/urRzWOmYUTe7hTG0yR3SHUyvvX0Lk53Z4cic8Ni/4k+wmH3vyoz/P21H7+Vdx3aZMpm9YajvPpW\n8Q7u8ar1kdGJw/i7v9zDL9a7S4pe3pptjsqdZ84M/+4/kb30YmffCM+8foQ//PtX8obkwV1+UygW\nT2I7DvXVFVSlh1A37j7DGwXz1mu3neTRX+2lMhrm4/dcRlXUrfAiYYMbzFl5H3s6Z2g3M188t7WG\nhpqKMWFwqmuIl7ac4Lk3jvFGzoucQsMjSXYc7PI2dDnVNUxtVYSVi5sBONk18fruyTqUHv6+fdVc\nAG/zl4HhUT77T+t58HH3BVruMHVTbU5lXFWqMs7+6c9pqfHeb22sontgZMKRjXMViyf51tO76B8a\n5apLWunqH+GrP9zm+8tkOo5DZ18sMMvgxN+0N/U0ON8nsqM588RDBRtcZMJjMJbgrx/dPOZzEyl7\nwg1CMp9baLxh6lwHc4Zdc3cUy22GOpsOvdww7uiNeZX0828e4yPvXk5DTQWGYeTNFY8mUkQrwl4n\ndV11hTeNpjGbAAAUc0lEQVQH/Y2ndnof11BTQf9wgqfSHcv/9TdWMbu5hvnttRw40U9jbZRoRZj3\nXLeADbvPEAkbdPaNMDKa5PEX9nkvjmY1VdNUV8mZnhiO43gjDRt3ZwN4894Oblo5G3Cr/S37O/nY\nu5eTSNr86TfWE4snGU46XLeshY7eGEvm1DO7pYaQYXCyc/y53VK2H+jiybUH+Xd3LefShY0cONHH\n7JYa5rXVUl9T4VXcz250G/+Onh3kTM8wfYOjVFeGqawI51XGxdYYZ0Rzlj3Nbs4J44YqDp0aoG9w\nlOb6yvE+9Zw98/oR+oZG+bXbl/Jvb1vCw7/cw6vbT/HS1hPcfcPCC/76tu14ozR1JYbnJyM+muJH\nLx9g894OegbizGur5Q9+/Spm57xokZnDcRyOnBng0KkBjp5x/53pjlFTFaGprtLrt5jdXMOc1hrm\ntNTQ2lBVdAe76aIwngbjNVZNRm74FNM/PDqm0oSxS5bArZYj4WwFNDgyfmVsO8Ur41x2+vJ8IcPI\na4Y62xtjMJbgVNcwlemqduv+Tu/xdTtOs27Hae66bj4fv8fMGwZdt+M0/UOjmAubAKirqRj3xcyc\n1lr6h3sZGkkyr63Wq0SXzWvkwIl+Ll3gfv79772Uj757Gd94aidb93eydtsp1m7Ljha0N1fTVF/J\n0bODjIymqE43Ou1PD51XRELsOdqL4zh09Mb4yg+3AdBUV0kyaXvD7mu3HGdRWw0p22FOSw0VkRCz\nmqvTnd3ZkB/Pic4hnll/hCVz6rn7xoWkbJvvPbuHrv44X35sC//xvssZGU1x8yL3Z1rQXsfuIz0M\njSRYs+W493U27j5Ld/+IN1fcmDNnXFOiMs69TGamgQvyO6rPJYxPdw8zOJxgXltt3hB5R2+MZzce\no7m+kvffvAjDMPjonct4c89Zfv7aYe5YNZeq6Pk9DTmOw2s7TvOjlw54IwQ3rZzFA/eupDJ6fpcB\nHYwl+NsntnHwZD911RWsWNTEnqO9/O9H3uAzH7qSK5e2ntfXPRcp26azb4RYPEl8NJV+QRYt/Ykz\nXM9AnAMn+qiqDFNbVUFNVYS2xqqiVx+bjL3HevnRywfYfzxbKIRDBrOaqxkZTXHoVP+42xZHwiHm\ntFQzpyUb0LNbapjVVE1ddcWEf9/nS2E8DSaat71QmYaqQuPN853sHGLb/k4uXdDE0rkN4wb20Ehy\nTBVejOO4c5RNdVH2Hcuvgg+kq+LbV83lhU3Hx9057OWtJ/noncu9J06A7z3rrqFdOtfdQaqhJjrm\nqlbgDqXuPeYG5ry2Wu/+D99xCbVVEe66fgHg7sccDYVpS4fK2m0n877OpQua2JZ+odAzEKe6MoLj\nOBw+PUB7UxXL5jXy+q4znOwa5o3d2Qawn75yEAODuuoK2puq2HWo2zueOa013nGd7h6mZyBOS7oZ\nqtBoIsWDj2+hd3CU9TtPc9WyVo6eGaCrP+69kPnW07sAMAvCeP2O08TiKW5aOYvNezt4/s1jDMeT\nXL3cDYim3Mq4euI/7UjYfTJprq/Me2LJDFkfOT3A8vmNvLDpOL96ww3T33zPpd5OX7l+tu4QP33F\nHbFoqI3yZ799nVdFPrFmP8mUzUfvXOa9AKivifK+mxbx1KuHWL3hKB+645IJj7VnIM4r205ypmeY\ncDhES30lFZEQr+86w4mOIaKRENcsb6Orf4SNu89yumuYP/rY1XkvTiajozfG3/1oOyc7h7j1yjl8\n6t4VRMIhXttxikdWW/zDT97i8x+/nkWzx56DqdDZG+OV7ad4ZfvJvC1ODQPMhU3cdPlsbr9qbt4L\n7CA4fLrfnR7afXZMMLY2VHHPTQt556p55/QC63jHID9+6QDbDrh9LNcsb+Pay9pYPLueeW213jl0\ne1US9A7EOd097P7rGuZU+u14jZPVlRFmNVczq6nafdtczezmGtqbqqmvqTjv/58J/2JN0wwBXwdW\nAXHg05ZlHch5/D7gfwFJ4J8ty/r2eR1FwEwUxkvm1Oc1Lp2rYo0vP3/t8Jj7/vzhN7z3qysj3HNj\n8SHB5vrKvK9d7Di7+kZYvcFdVgRwxdIWdh7q5k3LHea9cmkLb1pnvcaiz/32dWzb38nW/Z2c6hpm\ny77x1y5n5kjnt9WO+0p16dx61rpFKvNas8OFldEw9922dMzHz2quBtwXJHXVFXzw1iUsn9+YDlP3\nsVNdw8xrq6Wzb4ShkSSXL2lhxeJmXt91hj1Heti4+yzRSIi7rlvA6o1HcXB4300LveN9MX0O5rbW\nese4eW8H+0/0cVNBGDuOO5T68lb3iba1oYqu/hGeXneY0fTOXZ//xPV8/cm3vM1JzIVu9b9glvv1\nM8veVi1rZSiWYOdht4t9WbopraYyQiQcIpmyS84ZZ6YsWgqq30xT3FsHu7jusnaeWLOflO1wtjfG\nV5/Yxv/8xPW0pc8fwC/WH+anrxyirbGKlYubeWX7Kb78+Bb+6GPX0DMwwptWB8vmNXDz5bPzvs89\nNy5k7baT/GL9Ea67rH3cgOsbjPPI6j28su3UuKMl4ZDBTStn8dE7l9PaWEUyZfPor/by8taTfO0n\nb/HZ+6/N6xqfyK7D3XzjqZ0MxhLcc+NCPnbXckLpFym3XjmXyooI//jkW/z9j7fzvz55I421516p\n2o67oc/IaIqR0SSxeIrYaJLTXcNs2H3Gq96qK8PccsVsGmqihMMGe4/2sif977mNx/it917KlZdM\nf4VeSjJlc/BkP6e6hjjVNUz3QJz2xirmt9cyv62O+e21RYPJdhy27e/kuY3HsHJeZL/jitnetEPv\n4Chb9nbw2PP7+Nmrh7jrugW8+7r5easGCnX2xXjqlUO8tuM0DnDZwiY+cucyr3GzUMgwaKiJ0lAT\nHfM7mNn3/XTXEKe6hznbE3P/9cY40TE07gglQFU0TF11BXXVFUTSjZJf/aM7S5zN0pXxh4CoZVm3\nmqZ5M/Bg+j5M06wAvgLcAAwD60zT/JllWcW7X94m/vQb64s+dsm8hgsK42L25gzDjCcWT/LCJjc8\nDPKXU11/WTvhsMHG3WeZ1VzNJ9+/gvamKj77T2N/jv/7/U3e+x+8dTEV4RA7D3V7u40tndfAvNZa\nL4yXzW/gsoVNLJ5Tzzee2uk1mC2eXT9mLTXA/Pa6cbeqXDgr+4eSWxkXs3BWnff+JfMa8l6ILE7/\n0b2+8zQ9AyPE0t3Yy+c3siI9/P3cG0fp6B3hBrOde29ZxJEzbuX83usXesd99Myg93kA5iL3c62j\nvd6cM8CJjkH+4ckdnEk36NVVV/DFT93AXz262buM5uzmaha01/Lp+y7ntR2nWdBW6w0TZ36WzMVF\nli9oYiiW9MI48/0Nw6CpLkpn30jJMM688GoqCOO2xmrmt9Wy+0gPP3hxH6NJm0/du4LRRIp/fX4f\nf/2vW/jjf3c17U3VrN5wlJ+sPUhLQyWfvf9a2prcYb0nXjrAX3z3DRIpG8OA33zvpWOG9aorI3zq\n3hV89Yfb+PbPd/E/fvNaGtIBl0zZvLDpOD9/7TBDI0nmttZwz40LuWJJCynbobNvhMFYgiuWtuTN\nEUfCIf79+0xGEzbrd57m27/YzX/84OXecjbbcdh3rJfjHUN0D7jXfbYdhz1HejneMUg4ZPDJ95u8\nK90Zn+t6s50Pv/MSnlx7kAcf38rvffjKonPItuNwomMI62gPR04PcKxjkI7eGCPxVNFLsBjAikVN\nvOPKOdy0YvaYKrC7f4RfvH6El7ac4Cs/3MYVS1t4/82LuHxxs3du44kU+473svtwj7sTWzzJSNyd\nimmoraC+JpreJjVEJOy+bW6qITWapLWxivamamqrIiWHYDv7Yry89SSvbDuZt+qjUFU0zIpFzVyx\ntIX5bbU4jkPKcTjdNczzm457U11XLm3hnpvc/9/C7z0wPMqLm0/wwqbjPP3aYZ55/QjXm+3cdd0C\nLl3QiGEYJFPutq7PbjzKi5uPk0w5LGiv5SN3LuOqS1rPe0jZMAya6ytprq9k5ZKWvMdsx6F3IO6F\n89meGJ197nTd4HCCgViCE51DJFM2BpP7/qXC+DZgNYBlWRtM07wh57GVwH7LsvoATNN8FXgn8KPJ\n/ahvHw/cu4K29C/6RFtkZtx36xKeHqfSvVCZiwJ89reu5fvPWsxrr+OlLSe4+fLZ7Exv1BGLJ735\n2Fzz22rzrgT1b29bwofuuAQrZ91ze1MVDTVRZjdXs/tID5FwyJvvuXxJC4YBu9IBcs2lbV6o3X7V\nXNbvPM3Vy9torq/Mm88EtwN6QXs2gCcTxgsKwjjXwnQYb9rbwaacDUNWLWulvbGKeW21nEz/rDeu\nnE19TZQ/uf9a7+OWzq0nZIDtuEPDmRBZMqeeaEWINVtOsPNQN/W1FVx/2SxWbzxK/9Aoy+Y3kEja\n/IcPrKS+Jsqd187nsef3AXDtZe0YhsGyeY0sm5f/Kn5eay2G4U4TNNZFaW+sYtWyVh57YR+V0TDz\nc85NoxfGE/9pZ8J4vHnhq5a1snrDUTbuPsuC9lpuv2ouoZBBbDTFk2sP8j+/vYGayghDI0ma67NB\nDHDvLYtpbazi4Wf2MLu5hk/cc9mYn8f7Ppe0cue183lpywm+8NDr3LBiFomkzVsHuxgYdn9Xf+u9\nl3LntfPzKqyJGqkMw+BT966gsy/Gm3vOcqpziLuuX0BHT4w3rbN0jjOdEwkbrFrW6o2eFPPBdyym\nbzDOi5tP8OfffYNfu20pl8xroL6mgq7+EU52DrP3WC/W0Z68qZ9oup+gpjJCVWWE6soI1dEwVdEI\nVZVhGmqjXL2sbcI5+paGKj5xj8m7rp7HD17cz85D3ew81M3c1hqqomGG4ym6+kbGdH5nfm8mq666\ngiVz6lkyt96bq66vrmAglsA66v5s+0/04Thuk+B7rlvAojl1zG2tpaW+ko7eGCc7hzh2dpDdR3vZ\nmh4VG++c37FqLnffuJAF7XXjHImrvibKr92+lPffvIj1O07zwubjbNx9lo27z1IVDZNI2nkjaa0N\nVXz4nUu55fI509qAFTIMWhqqaGmo8l7AX6hSYdwA5O5ikDJNM2RZlp1+LLccGwCK/ya/jSyd25B3\njdn62qgXcLF4kl+96XbEfuET1/OX39s05vM//M5L2HGoyxu6zZhsSN98xRyMdAdh/qUa3TndBe11\nfO7j1+M4DnffsIA5LTXePO681uwT+5d+9xb3HQeqKiP80ddedY/vjqXce8tiAC5d2OT9vO+5zp23\nfd/Ni9h5uJv733uZ97Xqqiu4aeVsNuw6Q2VFmLuum8/PXztMynb4N+9YzMfuWu41/qxa1sotl8/m\n/TcvoqYyQk1VJG+4Mbf7t5jaqgouX9JMz0Cc26+am/dYY210TGW+dG699yR/z40L+e4v99CWDr1C\nFZEwH3uvyeO/srhhRXY5VSQc4n03LuLp1w67r5Z7Y95FOX777st4T3pe2/u/uGou+4710lhXyX23\nLin6s0QrwnzkzmW8daCLG1fOxjAMZrfUcOOKWbQ0VOY1uMxtqeXYmUHvBUIxv3HnMh57fp+3jjnX\n1ekwbm+q4g8/crX3pHbfrUtob6ripS0n6eiNcfPls/nALYvHzI/ftHI2V13SSmVFuOQT4sfvvox5\nrTX8eO1BXt7qzu831Ea558aFfPK+K4kPn/t65IpIiP/20av50csHWLP5hNeXEK0IcfuquVyxpIXW\nhioiEQPbdl/sVZfYsQzcoP/4PSaXLmjikdV7iq5eaG2o4prlbZiLmlk2v4HZzTVTFgyLZtfzJ/df\ny6FT/Ty78SibrA5CIYPqygjz22pZuaSZK5a0eI100UjIW744OJwgkbRJpGyS6bfVNVFOnR2kszdG\nZ98IJzoH2XGomx2Hxt9FzzDcpsl3XTOPG1fMGjMN0NJQ5Y0QgTsPvvNwN939ccIhAyNkUB0Nc+PK\n2ec01F9ZEebOa+fzrmvmsfdYLy9uPsHJriGqKsJEK8JuFb64mTuvmZ93WdiZxJhogwPTNB8EXrcs\n64n07WOWZS1Mv38V8FeWZf2b9O2vAK9alvWT6T9sERGR4Cj1EmId8AEA0zRvAbbnPLYHuNQ0zWbT\nNKO4Q9TFJ0tFRERkXKUqY4NsNzXAA8D1QJ1lWQ+ZpvlB4Iu4of4dy7L+aZqPV0REJHAmDGMRERGZ\nfjNzpltERCRAFMYiIiJlpjAWEREpM4WxiIhImV2UC0WU2uNapk5629K/sizr3eU+liBKbwP7z8Bi\noBL4P5ZlPV3eowoW0zTDwEPAZbg7t/5ny7JKX9JMzplpmrOATcB7LMvaW+7jCRrTNDeT3RzroGVZ\nv1PsYy/WVZuK7nEtU8c0zc8CHwcGy30sAfbbQIdlWZ8wTbMZ2AoojKfWBwHbsqzbTdN8F/CX6Pli\nyqVfWH4TGHtpIrlgpmlWAUy2MLpYw9R5e1zjXlxCpt5+4NdhkjuTy/l4AndtPbh/P5O7/qRMmmVZ\nTwG/m765BOgp/tFyAb4M/BNwqtQHynm5GqgxTfNZ0zRfSBeiRV2sMB53j+uL9L3fNtJbkSocppFl\nWUOWZQ2aplmPG8xfKPcxBZFlWSnTNL8L/D3wr2U+nMAxTfNTuCM8z6Xv0gv4qTcEfNmyrPcB/xl4\ndKLcu1iB2A/kXiwyc7EJkRnHNM2FwIvAv1iW9Xi5jyeoLMv6FO688UOmaVaX+HA5Nw8Ad5umuQa4\nBnjENM3ZJT5Hzs1e4FEAy7L2AV3A3GIffLHmjNcB9wFPjLPHtciMkX7Ceg74L5ZlrSn38QSRaZqf\nABZYlvUlIAbY6X8yRSzLelfm/XQg/65lWWfKeEhB9ABu0/LvmaY5D3eEuOiUwMUK4ydxX4WtS99+\n4CJ937cr7XE6fT6Pe6nQL5qmmZk7vteyrLEXypXz9SPgu6ZpvgxUAH9oWda5X0tRpLy+Azxsmuba\n9O0HJhoR1t7UIiIiZaYmKhERkTJTGIuIiJSZwlhERKTMFMYiIiJlpjAWEREpM4WxiIhImSmMRURE\nyuz/A712qdKWJIA2AAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 8 }, { "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", "collapsed": false, "input": [ "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]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [ "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');" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFsCAYAAAAtwdttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8leWd///XOSf7TsgCScgONwnIjkGIyuKGoGNda7Uu\nFe10ujgd+3Vm2p+1nS7TvePYOm3FDe2iVtECglbBBYTIvphwQ3bInpB9OSfnnPv3R0IKCsiS5GR5\nPx8PH8lZ70/Ibd7nuu/r/lw2y7IQERGR0cHu6wJERERk8Cj4RURERhEFv4iIyCii4BcRERlFFPwi\nIiKjiIJfRERkFPHzdQGnYhjGEuA2IAT4mWma+3xckoiIyIgwVEf8waZpPgD8ArjK18WIiIiMFLah\n2sDHMIxQ4HHgYdM0631dj4iIyEgw6If6DcPIAX5imuYiwzDswBPANMAJrDBNs8gwjBjgZ8B3Ffoi\nIiL9Z1AP9RuG8TDwJBDYe9cNQIBpmvOB/wB+2Xv/L4F44L8Nw7hpMGsUEREZyQZ7xF8I3Ag833s7\nF9gAYJpmnmEYc3q/v3uQ6xIRERkVBnXEb5rmq4D7hLvCgZYTbnt6D/+LiIjIAPD15Xwt9IT/cXbT\nNL3n8gaWZVk2m61/qxIRERnazjv4fB38W4DrgJcNw5gHnPP1+jabjbq61n4vTEav2Nhw7VPSr7RP\nSX+LjQ3/7Cedhq+C//g1hKuBKw3D2NJ7+14f1SMiIjIqDNnr+M+BpU/S0p80OpP+pn1K+ltsbPh5\nH+rXRDoREZFRRMEvIiIyiij4RURERhEFv4iIyCii4BcRERlFFPwDYNeuHVxzzUJqa2v67vu//3uc\n9evXnvY169evZfPm9897m1/72gOUl5ee9+ufeur3vPbaK+f9+tMpLi5k797dZ3zOK6+8BEBe3lb+\n9rfV/V6DiIj8g4J/gPj7B/DjH3+/7/ZndRdcunQ5ubmXnff2et7//DsYDlT3w02b3qGkpPiMz1m1\n6ikAcnIu4frrPzcgdYiISA9fd+4bcC9tLGT7wdp+fc+5k+O4dXHmaR+32WzMmjUHsHjllZe46aZb\nT3r8d7/7DaZZQHNzM5mZE/n2tx/lqad+z9ixMRw5Uk5m5kSWLl1OQ0M9Dz/8TZ566nl+97vfsG/f\nHrxeL7fd9gUWLbriFFu2eOONNWzZ8j4ul4uGhnpuueV2PvjgPYqLi/ja1x4kN/dy7r77diZMmEBN\nTTWZmZN4+OHv9L3D7t07ef75ZwgICKC2toZ/+qeb2LVrO4WFh7nlls9zww03s3v3Tp588v+w2+0k\nJibx//7ft3nrrfVs3boFp9NJZeVR7rjjbubOzWH9+rUEBARgGJOprq5i9eq/4na7sdls/PjHP+e1\n116hpaWFX/7yp2RnT6GsrJR//uev8ec/v8DGjW/hcPgxffpMvvKVr/PUU7+nurqKxsZjVFdX841v\n/BsXXzyvn36rIiKjg0b8A+B4U6SHHvoPXnrpT1RUHO17rKOjnYiICH7969+ycuUq8vMPUF9f1zfi\nvu66G9iwYR0Ab775BsuWXc/WrVuoqqrkiSdW8thj/8eqVU/T1tZ22u13dnby858/xh133M3q1X/l\nxz/+OQ8//G3WrVsDQHV1Jd/85sM8+eQqWltbef/9d096fV1dLT/60c956KH/5LnnnuKRR37AL37x\nv7z++qsA/PSnP+LHP/4Fv/nNH4iNjWP9+rXYbDba29v52c9+zU9+8iteeOFZYmJiufba67jttjvI\nyprC0aNH+PnP/4cnnlhJamoaeXnbuPvu+4iIiOChh/69b/tFRYVs2vQ2v/vdM/zud09z9Gg5H364\nGZvNRkBAAL/4xf/y4IMP8eKLf7rwX5aIyCgz4kf8ty7OPOPofCBFRETyjW88xA9/+CgXXTQdgICA\nQI4dO8b3vvcdgoND6OjowO3uWbDQZrORmpqGx+OhurqajRvf5rHHnuC1117BNA/y9a9/GaD38Soy\nMyd+aps2m42JEw0AQkPDSE1NAyA8PByXywVAamo6Y8fGADBt2nSOHCk76T3S0zNwOByEhYWRmJiE\nn58fYWE9r29sbOTYsQYeeaQnqJ1OJ3Pn5pCUNIGJEycBEBsb17cty7L6PghFRY3hhz/8HsHBwZSX\nlzF16rRT/ruVl5cyZcpFOBwOAKZPn0lJSRFA3zbi4uJxuZxn/bsQEZEeGvEPsAULLiU5OaVvYt+2\nbR9SV1fD9773Ix544F9wuZx9wXj867Jl1/PEE4+RlpZOaGgYKSlpzJo1m8cf/z2//vVvWbToChIS\nEk+7zc86X3/0aHnfEYP9+/eRnv7JD0anf31UVBRxcXH89Ke/4vHHf8+dd97DnDkXn3a7DocDy7Jo\na2vj6af/wH/913/z7//+/xEYGNj3nONdo4///CkpqeTnH8Dj8WBZFnv27GbChJTPrE1EZKQrrGjm\nJ3/cdUHvMeJH/L5gs9lOCsEHH3yInTu3A5CdPYXnnlvJN77xz0RHjyU7eyr19XV9rwNYtOgKHnvs\nl/z0p78GIDf3Mnbv3slXv3o/nZ0dXHbZIkJCQs64/RO//uP+nq8BAYH88Iff5dixY0ybNoP583M5\neDD/lK/75Pc2m40HH3yIb33rQSzLS2hoGN/5zveprq76xPZ6vjeMyfz2t/9LSkoqF100nS9/+V7G\njBnDhAnJfT93amoaP/jBI8yZk4PNZiM9PZPFi6/gK1+5D8vyMm3aTC67bCGFhYdOW5uIyEjm9nh5\nfXMJb2wr+8cyd+dJi/SMQnfddRurVr3o6zKGLC2oIv1N+5RciIr6dp5c8zHlNW3ERAaxYnk2C2ZN\nOO+Rj0b8o5BGyiIiQ5/Xsnh7x1H++m4Rbo+X3GnjuX3JRIIDLyy6Ffyj0HPP/cXXJYiIyBkca+ni\nqXUFFJQ1Eh7izz3XTGHmpNh+eW8Fv4iIyBBhWRZ5+TU8/9YhOp1upmeM5Z5rs4gMDei3bSj4RURE\nhoC2zm6ef9Nk+8FaAv0d3LN0MpdOG9/vp2cV/CIiIj52oLiBp98ooKnNRWZiJCuWZxE35vRXb10I\nBb+IiIiPOLs9vLypkI27KnDYbdx0eTpLc1Kw2wduEraCfwBUVVVy9923YxiT++6bPXsu99yzgoUL\n5/V18YOeLnpTp17E2rWv43K5KC0tZtKkntc9+ugPiYmJ7XvP733vO/z+98+cd10333wdf/7zq/j7\n+5/3e5ytmppqCgsPs2DBpad9zuuvv8qyZddTUlLMli3vc889Kwa8LhGRoaKkqoU/rMmn5lgHCTGh\n3L88m5Rx4QO+XQX/AElLS+fxx3//qfsjIyNPef/VV19LdXUVjz767VM+3h8G8zK+nTu3U15edsbg\nf+GFZ1m6dDkTJ07qa8UrIjLSuT1e1m0tY82WUryWxVVzJ3DT5en4+zkGZfsjPvhfLVzL7tr9/fqe\nM+Mu4sbM5f36nvCPlrWf5Wtfe4CJEw2Ki4sICQlm2rSZfPTRVtraWvnVr37LBx+8S17ehzQ1NdPc\n3MSXvvQAl122sO/1NTXV/PznP8bpdBIYGMjDD38Hj8fDd7/7n8THj6O6uoolS66ipKSIQ4dMLrlk\nAV/+8lcpKirkscd+gWVZREZG8p//+V1M8yB//OMqAgL8qaysYMmSq7jzznt44YVncTqdTJ06jdDQ\nUJ59diVer5fOzk4effSH7N27i4aGBr73ve9wyy2f57XXXuH73/8xb721npdf/jP+/gEkJU3g4Ye/\nc8qV/5Yu7f9/fxGRgVZ9rIMn1+RTUtVCdEQg912bRVZq9KDWMOKD31dKS4v7FtUBePTRHxETE0NL\nS8tJ93/ta9886ZTA2bDZbGRnT+HBBx/ioYe+QXBwEL/+9W/50Y++x549O7HZbHi9Fo899gQNDfV8\n+cv39o28Lcvit799jJtv/jzz5s1nx46P+N3vfsMDD/wLVVWVPPbYE3R1dXHLLdfz2msbCAwM5Oab\nr+PLX/4qP/3pD/nOd75HSkoqa9e+zh//uIq5c3Ooqalm1aq/4HK5uOGGa7jrri/xxS/eS3l5Gbm5\nl7F69V955JEfEBMTw/PPP8OmTW9z111f4rnnnub73/8x+/fvBaClpZmnn/4DzzzzJ4KDg3n88V/x\n+uuvEhISQnt7O7/61eMcPXqEf//3byr4RWRYsSyLTbsreGljIS63l0umxHPHlZMICRr4U6+fNOKD\n/8bM5QMyOv8sqamnPtQfERHRL4fyj39YCAsLIzU1HTh5Bb7Zs+cCMHZsDGFh4TQ3N/W9tri4kOef\nf4Y//vE5LMvqO+efkJBISEgoDocf0dFjCQ/vOdd0/AxBWVkJv/jFfwPgdruZMCEZgIyMDOx2O0FB\nQX2L75y4Kl9MTAz/8z8/JyQkhLq6WqZNm3HKn6mysoK0tHSCg4MBmD59Fh99tI0pU6aecuU/EZHh\noLHVyTNvFHCg5BihQX7ctzybuZPjfFbPiA/+kevM5+sPHswHbuLYsQa6urqIihrT91hKSiq33/5F\npk6dRnFxIfn5B3re8TPmACQnp/LII/9FXFw8e/bsorm5+bS12O12vF4vAD/72Y956aXXCQ4O5kc/\n+l7f/T1HJjx9rxk/PoGSkhK6uroICgpi9+6dJCennFVtIiJD0faDtazacJD2LjdT06K599osxoQH\nfvYLB5CCf4CcPqjOHGBnCrhzCb+jR4/w4IP/QkdHG9/61n9gt9uBntX1vvrVf+UXv/gJLpcTp9PJ\nv/7r/zvF+3/6+2996z/5wQ++i8fjwW638x//8Qh1dbWnfF16egarVj2NYUzmqquW8tWvriAmJpbk\n5FQaGuoBmD59Jt/61oPce+/92Gw2IiOjuO++B/j617+M3W4nKWkCX/nK13nnnbfOUJuIyNDT0dXN\nC38/xLaPawjws/PFqyaxcGbikBjEaHW+EWj9+rU0NTVx++13+rqUYUkrqUl/0z41uhSUHmPlugIa\nW52kjY/g/uuyGRfdv814YmPDtTqfnGwIfKgUERlVXN0eXnmvmL/vOILdZuOG3DSWzU/BYbf7urST\nKPhHIM14FxEZXGXVrTy5Np/K+nbGRYdw/3XZpI2P8HVZp6TgFxEROU9er8Ub28p4fXMJHq/FkllJ\n3Lwog0D/wWnGcz4U/CIiIuehtrGDlWsLKKxoJiosgC8ty2Jq2lhfl/WZFPwiIiLnwLIs3t9byV/e\nKcTZ7eHirDjuvMogLHjwm/GcDwW/iIjIWWpud/HsGwXsLWogONCPB67PZl72OF+XdU4U/CIiImdh\n16E6nl1/kLbObrJSxnDfsiyiI4J8XdY5U/CLiIicQafTzZ/fPszm/VX4+9m5fclElsxJwj5Mr5tW\n8IuIiJyGWd7IU+sKqG/uIiU+nBXXZZMYE+rrsi6Igl9EROQTut1eXvugmA155WCD5fNTuH5BGn6O\nodWM53wo+EVERE5wpLaNJ9fkc7SujbioYFZcl01mYqSvy+o3Cn4RERF6mvG8ub2c1e8X4/ZYLJyR\nwK2LMwkKGFlRObJ+GhERkfNQ39TJynUFHDrSRERoAPcuncz0zBhflzUghnTwG4axGLjdNM37fV2L\niIiMPJZlsWV/NX96+xBdLg+zJsVy1zUGESEBvi5twAzZ4DcMIwOYAQy/iyRFRGTIa+lwsWqDya5D\ndQQFOLhvWRbzp47DNkwv0ztbQzb4TdMsAn5lGMbzvq5FRERGlj2F9Ty7/iAt7S4mTYhixbIsYqKC\nfV3WoPBJ8BuGkQP8xDTNRYZh2IEngGmAE1jRG/oiIiL9qsvl5sWNhby3pxI/h41bF2Vy1dwJ2O0j\ne5R/okEPfsMwHgbuBNp677oBCDBNc37vB4Jf9t4nIiLSbwormlm5Jp/apk6SYsN44LpskuLCfF3W\noPPFiL8QuBE4fgg/F9gAYJpmnmEYc058smmaXxzc8kREZCRxe7z8bUsJ67aWgQVL5yVzQ246/n7D\nvxnP+Rj04DdN81XDMFJPuCscaDnhtscwDLtpmt7BrUxEREaaivp2Vq7Jp6ymlZjIIFYsz2bShChf\nl+VTQ2FyXws94X/cOYd+bGz4Zz9J5Bxon5L+pn1qcHm9Fms3F/Psuny63V6uvDiZFf80lZAgf1+X\n5nNDIfi3ANcBLxuGMQ/Yd65vUFfX2u9FyegVGxuufUr6lfapwXWspYun1hVQUNZIeIg//3z9FGZO\niqW9tYv21i5fl9cvLuSDpC+D3+r9uhq40jCMLb237/VRPSIiMoxZlkVefg3Pv3WITqeb6Rljuefa\nLCJDR24znvNhsyzrs581tFn6JC39SaMz6W/apwZeW2c3z79psv1gLYH+Dm6/YiKXThs/YpvxxMaG\nn/cPNhQO9YuIiJy3A8UNPP1GAU1tLjITI1mxPIu4MSG+LmvIUvCLiMiw5Oz28PKmQjbuqsBht3HT\n5ekszUkZVc14zoeCX0REhp2Sqhb+sCafmmMdJMSEcv/ybFLG6cqJs6HgFxGRYcPt8bJuaxlrtpTi\ntSyumjuBmy5Px9/P4evShg0Fv4iIDAvVxzp4ck0+JVUtREcEct+1WWSlRvu6rGFHwS8iIkOaZVls\n2l3BSxsLcbm9XDIlnjuunKRmPOdJwS8iIkNWY6uTZ94o4EDJMUKD/LhveTZzJ8f5uqxhTcEvIiJD\n0vaDtazacJD2LjdT06K599osxoQH+rqsYU/BLyIiQ0pHVzd//Pshtn5cQ4CfnS9eNYmFMxNHbDOe\nwabgFxGRIaOg9BhPvVHAsRYnaeMjuP+6bMZFqxlPf1Lwi4iIz3W7PbzyXjFvbT+C3Wbjhtw0ls1P\nwWG3+7q0EUfBLyIiPlVW3cqTa/OprG9nXHQI91+XTdr4CF+XNWIp+EVExCe8Xov1eWW89kEJHq/F\nkllJ3Lwog0B/NeMZSAp+EREZdLWNHaxcW0BhRTNRYQF8aVkWU9PG+rqsUUHBLyIig8ayLN7fW8lf\n3inE2e1h7uQ4vni1QViwmvEMFgW/iIgMiuZ2F8++UcDeogaCA/144LpscrLjdZneIFPwi4jIgNt1\nqI5n1x+krbObrJQx3Lcsi+iIIF+XNSop+EVEZMB0Ot38+e3DbN5fhb+fnduXTGTJnCTsGuX7jIJf\nREQGhFneyFPrCqhv7iIlPpwV12WTGBPq67JGPQW/iIj0q263l9c+KGZDXjnYYPn8FK5fkIafQ814\nhgIFv4iI9JsjtW08uSafo3VtxEUFs+K6bDITI31dlpxAwS8iIhfM67V4c3s5q98vxu2xWDgjgVsX\nZxIUoJgZavQbERGRC1Lf1MnKdQUcOtJERGgA9y6dzPTMGF+XJaeh4BcRkfNiWRYfHqjmj38/RJfL\nw6xJsdx1jUFESICvS5MzUPCLiMg5a+1wsWqDyc5DdQQFOLhvWRbzp45TM55hQMEvIiLnZG9hPc+s\nP0hLu4tJE6JYsSyLmKhgX5clZ0nBLyIiZ6XL5ebFjYW8t6cSP4eNWxdlctXcCdjtGuUPJwp+ERH5\nTIUVzaxck09tUydJsWE8cF02SXFhvi5LzoOCX0RETsvt8fK3LSWs21oGFiydl8wNuen4+6kZz3Cl\n4BcRkVOqqG9n5Zp8ympaiYkMYsXybCZNiPJ1WXKBFPwiInISr2Xxzo6jvPxuEW6Pl9xp47l9yUSC\nAxUZI4F+iyIi0udYSxdPrSugoKyRsGB/7lk6hVmTYn1dlvQjBb+IiGBZFnn5NTz/1iE6nW6mZ4zl\nnmuziAxVM56RRsEvIjLKtXV28/ybJtsP1hLo7+CepZO5dNp4NeMZoRT8IiKj2IHiBp5+o4CmNheZ\niZGsWJ5F3JgQX5clA0jBLyIyCjm7Pby8qZCNuypw2G3cdHk6S3NS1IxnFFDwi4iMMiVVLfxhTT41\nxzpIiAnl/uXZpIwL93VZMkgU/CIio4Tb42Xth6Ws/bAMr2Vx1dwJ3HR5Ov5+Dl+XJoNIwS8iMgpU\nNbTz5Jp8SqtbiY4I5L5rs8hKjfZ1WeIDCn4RkRHMa1ls3NnTjKfb7WX+1HF84YqJhAT5+7o08REF\nv4jICHWspYun3yggv7SnGc/9y7OZMznO12WJjw3J4DcMYz7wQO/NB03TbPZlPSIiw4llWWzLr+GF\nE5vxLJ1MZFigr0uTIWBIBj9wPz3BnwPcBvzBt+WIiAwPbZ3drHrTZIea8chpDNXgd5im6TIMowpY\n7OtiRESGg31FDTzzRgHN7S4ykyJZsTybuKhgX5clQ8ygB79hGDnAT0zTXGQYhh14ApgGOIEVpmkW\nAR2GYQQACUD1YNcoIjKcdLncvLSxkHf3VOKw27h5YQbXXJysZjxySoMa/IZhPAzcCbT13nUDEGCa\n5vzeDwS/7L3vD8Dve+v78mDWKCIynBRWNLNyTT61TZ0kxYayYnk2yfFqxiOnN9gj/kLgRuD53tu5\nwAYA0zTzDMOY0/v9LuDeQa5NRGTYcHu8vL65hDe2lYEFS3OSueHSdPz97L4uTYa4QQ1+0zRfNQwj\n9YS7woGWE257DMOwm6bpPZf3jY3Vp1vpX9qnpL/15z5VVtXCr/60m+LKZuKjQ/jm7bOYkj62395f\nRjZfT+5roSf8jzvn0Aeoq2vtv4pk1IuNDdc+Jf2qv/Ypr2Xx1kdHePX9Itwei0unjefzSyYSHOin\nfXaUuZAPkr4O/i3AdcDLhmHMA/b5uB4RkSGpvqmTp9YVYB5pIiLEn3uWZjFjYoyvy5JhyFfBb/V+\nXQ1caRjGlt7bOq8vInICy7LYvL+KP799mC6Xh1mTYrnrGoOIkABflybDlM2yrM9+1tBm6RCX9Ccd\n6pf+dr77VEu7i+c2HGT34XqCAhzcceUk5k8dp2Y8Qmxs+HnvBL4+1C8iIqew+3Adz60/SEtHN5OT\no/jSsixiItWMRy6cgl9EZAjpdLr58zuH2byvCj+Hnc8vzuSKuROwa5Qv/UTBLyIyRJjljTy1roD6\n5i6S48O4f3k2ibFhvi5LRhgFv4iIj3W7Pax+v4Q3PyoHGyyfn8L1C9Lwc6gZj/Q/Bb+IiA+V17Ty\n5Np8KuraiRsTzP3Ls8lIjPR1WTKCKfhFRHzA67VYn1fGax+U4PFaLJqZyK2LMgkMcPi6NBnhFPwi\nIoOstrGDlWsLKKxoJjIsgC9dm8VFarkrg0TBLyIySCzL4r29lbz4TiHObg9zJ8fxxasNwoL9fV2a\njCIKfhGRQdDU5uTZ9QfZV9RASKAfD1yfzbzscb4uS0YhBb+IyADbsreS37y8h7bObqakjuHea7OI\njgjydVkySin4RUQGSEdXN3/8+yG2flxDgJ+dO66cxKJZiWrGIz6l4BcRGQD5pcd4al0Bja1OJiVH\ncffVBuPHhvq6LBEFv4hIf3J1e/jru0W8vfModpuNG3LTuOf6qRw71u7r0kQABb+ISL8pqWph5dp8\nqho6GD82hBXLs0kbH4FDHfhkCFHwi4hcILfHyxtby1jzYSker8UVc5K4+fIMAvzVjEeGHgW/iMgF\nqGpoZ+XafEqqWhkTHsh9y7LITo32dVkip6XgFxE5D17LYtOuCl7eVIjL7eWSKfHcceUkQoLUjEeG\nNgW/iMg5OtbSxTNvFPBxaSNhwf6sWJ7NnMlxvi5L5Kwo+EVEzpJlWeTl1/DCW4focLqZljGWe5ZO\nJios0NeliZw1Bb+IyFlo6+zm+TdNth+sJdDfwV3XGFw+PQGbmvHIMKPgFxH5DPuLG3j6jQKa21xk\nJkayYnkWcWNCfF2WyHlR8IuInIbT5eHFTYW8u7sCh93GTZenszQnBbtdo3wZvhT8IiKnUFjRzMq1\n+dQ2dpIYG8r9y7NJjg/3dVkiF0zBLyJyArfHy9+2lLBuaxlYcM3FyXzusjT8/dSMR3yrobORXbV7\nOXjsMP911b+d9/so+EVEelXUtfHk2nzKa9qIiQzivmVZGMljfF2WjGJNzmZ21+5nZ81eSlrKALDb\nLqwFtIJfREY9r2Xx9+1HeOW9YtweL7nTxnP7kokEB+pPpAy+Vlcbu2v3sbN2L0VNpVhY2LAxaUwm\nc+KmMz1u6gW9v/ZqERnV6ps7eXpdAQfLmwgP8eeea6Ywc1Ksr8uSUaa9u4M9dfvZVbMPs7GwL+zT\nI1OZHT+dmXEXERHQP3NMFPwiMipZlsWHB6r509uH6HR6mDkxhruvmUxEaICvS5NRotPdyd66j9nZ\ne97ea3kBSItIZlb8dGbFTSMqMLLft6vgF5FRp6XDxaoNJrsO1REU4OBL12ax4KJxasYjA67L7eRA\nfT47avdS0GDitjwATAhPZHZcT9iPDR7YRZ4U/CIyquw5XM+z6wto6ejGmBDFfcuyiIkK9nVZMoK5\nPC4ONBxkV81eDjQcpNvbDUBC6Dhm947s40IG7/SSgl9ERoVOp5u/vHOYD/ZV4eewc9viTK6cOwG7\nRvkyALq9bvIbTHbV7mVffT4ujwuA+JBYZsVNZ3b8dMaHxvukNgW/iIx4h440sXJtPvXNXSTHhbHi\numySYsN8XZaMMB6vh4ONh9lZs5e9dR/T5ekCICYomllJ05kdN53EsPE+P6Wk4BeREavb7WX1B8W8\nmVcONlh2SQr/lJuGn+PCroMWOc7j9XC4qbg37A/Q7u4AYExgFAsSLmZ2/HSSw5N8HvYnOmPwG4YR\nC3wVuB6YCHiBQuA14P9M06wf8ApFRM5DeU0rK9fmc7SunbioYFYszyYzqf9nSMvo47W8FDWVsLN2\nH3tq99Pa3QZAZEA4C5MWMDt+OqkRyRfcaGegnDb4DcP4KnAj8CpwD1AGdANpwCJgtWEYL5um+b+D\nUKeIyFnxei3W55Xx2gcleLwWC2cmcuuiDIICdIBTzp9lWZS0lLOrZi+7avfR7GoBIMw/lEsTL2F2\n3DQyotKGbNif6Ez/J1SYprnkFPd/3PvfbwzDuGlgyhIROXe1jR2sXFdA4dFmIkMDuPfaLKZljPV1\nWTJMWZZFeetRdtbuZVfNPhqdTQCE+AUzf/xcZsVPZ1JUBg778FrHwWZZ1hmfYBiGH3CtaZp/Mwwj\nhp7D/s+YpnnmFw4eq66u1dc1yAgSGxuO9qnhxbIs3t9byV/eKcTZ7WHO5DjuutogLNjf16UB2qeG\nE8uyqGw8QRwzAAAgAElEQVSvZmfNXnbW7qW+swGAIEcQ02OnMCtuGpOjJ+Jn9+0RpNjY8POeNHA2\nlT8JOIC/ATZgCZADfPl8Nyoi0l+a25w8s/4g+4oaCA704/7rspmXHT+kJlPJ0FfdXtMb9vuo6agF\nIMARwJz4GcyKm0529CT8HUPjg+SFOpvgn2ua5lQA0zTrgDsMw9g/sGWJiHy2HQdrWfWmSVtnN1kp\nY7hvWRbREUG+LkuGibqOBnbW7mVnzR4q26sB8Lf7MSP2ImbHT2fq2MkEOEZeC+ezCX6bYRgJpmlW\nAhiGEQ94BrYsEZHT6+hy88e/H2Lrx9X4+9n5whUTWTw7Sc145DMdX9N+V+1eylsrAPCzObgoJpvZ\ncdO5KCaLIL+R/eHxbIL/R8AuwzC29N7OAR4cuJJ6GIaxGLjdNM37B3pbIjJ8FJQe46k3CjjW4iRt\nfDgrlmczfmyor8uSIazJ2cyu2n3sqtlLSUs50LOmfXa0waz46UyPmUKI/+hp2/yZwW+a5p8Mw3gP\nmEfP5XxfM02zaiCLMgwjA5gBjOyPXSJy1lzdHv76XhFv7ziK3Wbjn3LTWHZJiprxyCmdbk17Y0wm\ns3vXtA/zH50fGD8z+A3DCKTnOn4D+AbwDcMwfmKapmugijJNswj4lWEYzw/UNkRk+CitbuHJNflU\nNXQwLjqE+6/LJm18hK/LkiGmrbudvbUH2Fm7l0ONRSetaT8nfjoz+nFN++HsbA71/xaoA2YDbno6\n+D0FfPFcNmQYRg7wE9M0FxmGYQeeAKYBTmCFaZpFhmH8AMgEvmKaZtO5vL+IjDwer5d1W8tYs6UU\nj9fiitlJ3LQwg0D/4XXdtAwcX61pP5ydTfDPNk1zpmEY15im2WYYxl3AgXPZiGEYDwN3Am29d90A\nBJimOb/3A8EvgRtM03zkXN5XREau6mMdPLkmn5KqFsaEB/KlZVlMSR3YdcpleOhyO9lfn8/O065p\nP52xwWN8XOXQdTbB7zUM48TrGWLo6dl/Lgrpaf97/NB9LrABwDTNPMMw5pzqRaZpntNRBREZ/izL\nYuOuCl7eVIjL7WXelHjuuHISoUEj4xpqOT8nr2lfQLfXDZy4pv104kJifFzl8HA2wf8Y8DYwzjCM\nx4DPAd8/l42YpvmqYRipJ9wVDrSccNtjGIbdNM1z/UAB9HTFEulP2qd8o6G5k8f+spvdh+oID/Hn\nm1+YRe70RF+X1S+0T527bk83e6rz+bB8Bzsq9+N0OwFIDB/HJcmzmZ88m6SI8T6ucvg5m1n9qwzD\n2EnPwjx2YLlpmvsucLst9IT/cecd+oBaYUq/UntV38jLr+GFt0zau9xclD6We6+dTFRY4Ij4XWif\nOntnWtP+8sT5J69p7xy9f/8v5IPk2czqHwskmKb5G8Mwvg08YhjGo6Zp5p/3VmELcB3wsmEY84AL\n/SAhIsNUW2c3L7xl8lFBLQH+du662uDyGQlquTuKnHFN+8SLmR039Na0H87O5lD/n4E1hmFYwM3A\n/wC/Ay47j+0dX9hnNXDlCU2B7j2P9xKRYe5AcQNPv1FAU5uLjMQIVizPJn5MiK/LkkHw2WvazyA1\nYsKwWOZ2uDmb4B9jmubjhmE8DjzXe+j/G+e6IdM0S4H5vd9bwFfO9T1EZGRwujy89G4hm3ZV4LDb\nuPGydJbOS8Zh1x/5kez4mvY7a/awu3Yfza6ew/TDcU374exse/XPpucSvIWGYcw4y9eJiHxKUUUz\nK9fmU9PYSWJMKCuWZ5MyThPfRqrTrWkf6hfC/PEXMzt+OhOj0ofdmvbD2dkE+L8DPwd+2dtk50Pg\n3wa2LBEZadweL3/bUsq6raVgwdUXT+DGy9Lx99Mf/JHGsiwq2qp6wr5230lr2ueMm83s+OlMHjNR\nYe8jNsuyTvmAYRjjP6sn/9k8ZxBYo3VWpwwMzcDufxX17axck09ZTStjI4JYsTwLI3n0NFgZLfvU\n6da0nxaTPeLWtPe12Njw857peKYR/38bhlFBz3n9Qyc+YBhGFvAlYDw9HflERD7Fa1m8vf0If32v\nGLfHy4KLxvGFKyYRHKizhSNFbUc9u2r3srNm76ha0344O+3/faZp3mMYxnLgScMwJgGV9PTqTwKK\ngJ+bprlmcMoUkeGmobmLp9blc7C8ifAQf+6+ZgqzJsX6uizpB1rTfng748du0zTXAmsNw4gGMuhp\n1VtimuaxwShORIYfy7L48EA1f3r7EJ1ODzMyY7hn6WQiQjXqG85Ou6b9WIPZcdOZNsrWtB/Ozup4\nW2/QK+xF5IxaOlw8v8Fk56E6AgMc3Lt0MrnTxqvxyjDV4mplT+1+dtTspbhZa9qPFDrRJiL9Yk9h\nPc+uP0hLu4tJE6K4b1kWsVEaAQ43p1vTPiMqldlxWtN+JFDwi8gF6XS6eXHjYd7fW4Wfw8atizK5\nau4E7HaN8oeLju5O9tZ/zK6avRxs1Jr2I93Z9Op/xTTNmz5x3zumaS4ZuLJEZDg4dKSJlWvzqW/u\nYkJcGPcvzyYpLszXZclZ6HJ3sb++4FNr2ieHJzJLa9qPaKcNfsMwVgMzgATDMEo+8ZrygS5MRIYu\nV7eH1zeXsCGvHGyw7JIUrl+Qhr+fWq0OZZZlcbipmC2VeeytO/CJNe1nMCtumta0HwXONOK/BxgD\n/C/wdeD4cTs3UD2wZYnIUOS1LD7Kr+GV94poaHESGxXEiuXZTEyK8nVpcgZt3e3kVe1kS2UeNR11\nAMSHxDI7bjqz46czLjTexxXKYDrTdfzNQDNw/eCVIyJD1aEjTby48TAlVa34OWxck5PM9QtSCQrQ\nVKGhyLIsippL2Vyxjd11+3F73fjZ/ZgbP5PcxHlkRKbqaotRSv/HisgZ1TR28NdNRew81DNSvDgr\njpsuz9CM/SGqo7uDvOpdbK7Mo7q9BoC4kBhyE+aRM362Lr8TBb+InFpbZzdrtpSycddRPF6LjMQI\nPr94IhmJmt091Bxf7nZzxTZ21e6l2+vGYXMwO246uYnzmBiVrtG99FHwi8hJ3B4vG3ceZc2HpbR3\nuYmJDOKWRZnMMWIVHkNMp7uTj6p3s7liW1+f/NjgsSxIyGHe+DmEB+gKC/k0Bb+IAD2jxp1mHX99\nt4japk6CA/24dVEmS2Ynabb+EGJZFmWtR9hckcfOmj24vN3YbXZmxk0jNyGHSWMysNv0+5LTU/CL\nCCVVLfzlncMcPtqMw27jitlJXJ+bRliwllAdKrrcXWyv2c3mijyOtlUCMDYomtyEHOYlzFE3PTlr\nCn6RUay+uZNX3ytmW37PJLCZE2O4ZVEm46JDfFyZHFfecpTNldvYXrMHl8eF3WZnRuxUchPmYURn\nanQv50zBLzIKdTrdrNtaxlvbj+D2eEmJD+fzSzIxktWpbSjocjvZWbOHzZXb+pa9HRMYxVXJi7gk\nYY7a58oFUfCLjCIer5f391Ty2uYSWju6GRMeyE2XpzNvyjjsmrjnc0daK9lcuY0d1bvp8jixYeOi\nmGxyE3LIHmtodC/9QsEvMgpYlsW+ogZe2lRIVUMHgQEOPndZOlfNnUCgv8PX5Y1qTo+LnTV72Vy5\njbKWIwBEBUayOPky5o+fy5ggdUWU/qXgFxnhymtaeXFjIQVljdhscPmMBG7ITSMyLNDXpY1qlW3V\nbK7cxkfVu+h0d2HDxtSxk8lNnEd2tIHDrg9kMjAU/CIjVGOrk9UfFLNlXxUWMDU9mlsXZZIUq2u7\nfcXl6WZ37T42V26juLkMgMiACBamLmB+wsVEB2mOhQw8Bb/ICON0edjwUTnr88pwdXtJjA3ltsWZ\nTE0b6+vSRq2jLVWsObSRvOqddLg7sWEjO9ogNzGHqWOzNLqXQaXgFxkhvF6LLQeqWP1+MU1tLiJC\nA7h9SRqXTkvAbtfEvcHW7elmd91+NlfkUdTcs7J5eEAYV6csZn7CxcQER/u4QhmtFPwiI8DHpcd4\naWMhR2rbCPCzs3x+KktzkgkO1P/ig62mvZbNlXnkVe+kvbsDgGnxWVwcO4dpMdka3YvP6a+CyDBW\nWd/OS5sK2VfUAMD8qeO48bJ0oiOCfFzZ6NLtdbO37gCbK7ZxuKkYgDD/UK5MXsiChByyU1Kpq2v1\ncZUiPRT8IsNQS7uL1zeX8N6eSryWxeTkKG5bPJGUcWrbOphqO+rZUpnHtqodtHW3AzApKoPcxBym\nx07Fz64/sTL0aK8UGUa63R7e2n6EdVvL6HJ5iI8O4dZFGczIjNHKeYPE7XWzrz6fzRXbMBsLAQj1\nD2HJhMtYkJhDfEisjysUOTMFv8gw4LUsPsqv4ZX3imhocRIW7M8dV2Zw+YwE/Bzq5jYY6jsb2FL5\nEVurttPqagMgMyqN3IR5zIidir9DCxrJ8KDgFxniDh1p4sWNhympasXPYeOai5NZPj+FkCAFzUDz\neD3sbyhgc8U2Dh47jIVFiF8wiybkkpuQw7jQeF+XKHLOFPwiQ1RNYwd/3VTEzkN1AMydHMfNCzOI\njQr2cWUjX0NnIx9WfcTWyo9odvVMykuPTCU3IYeZcdMI0OhehjEFv8gQ09bZzdoPS3ln51E8XouM\nhAhuWzKRzEStyDaQPF4PHzccZHNlHvkNJhYWwX5BXJ60gNyEHBLCxvm6RJF+oeAXGSLcHi8bd1Ww\nZksJ7V1uYiKDuHlhBnMnx2ni3gBq7Griw8qP+LBqO03OZgBSI5LJTchhdvx0AhwBPq5QpH8p+EV8\nzLIsdh2q4+V3i6ht7CQ40I9bF2WyZHYS/n6auDcQvJaX/AaTzZXbOFB/EAuLIEcglyZeQm5CDknh\nCb4uUWTAKPhFfKikqoUX3znMoaPNOOw2lsxO4voFqYSHaJQ5EJqczWyt3M6Wyo9odDYBkByeRG5i\nDrPjZhDkpxULZeRT8Iv4QENzF6+8V8S2/BoAZk6M4ZZFmYyLDvFxZSOP1/JScOwwWyq2sb+hAK/l\nJdARwIKEHHITc0gOT/J1iSKDSsEvMog6nW7WbS3jre1HcHu8pMSHc9viTCanaDnW/tbsbGVr1XY+\nrMyjoasRgAlhCSxInMfc+BkE+amtsYxOCn6RQeDxenl/bxWvfVBMa0c3Y8IDuenydOZNGYddE/f6\njdfyYjYWsrkij331H+O1vATY/Zk/fi65ifNIDk/SREkZ9YZc8BuGsQS4DQgBfmaa5j4flyRy3izL\nYn9xAy9tKqKyvp1Afwefuyydq+ZOINBfq7T1l1ZXG1ures7d13f2LFiUGDae3IQc5o6bSbCfeh+I\nHDfkgh8INk3zAcMwZgBXAQp+GZbKa1p5aVMh+aWN2Gxw2fQEPndpGpFhmkDWHyzL4lBjEZsrt7G3\n7mM8lgd/uz/zxs0hNzGH1Ihkje5FTmHIBb9pmmsNwwgFvgE87Ot6RM5VY6uT1R8Us2VfFRYwNS2a\nWxdnkhQb5uvSRoQ2VzvbqnewpSKP2s56AMaFxnNpwjwuHjeTEH9NkBQ5k0EJfsMwcoCfmKa5yDAM\nO/AEMA1wAitM0ywyDOMHQCbwIPAT4LumadYPRn0i/cHp8rDho3LW55Xh6vaSGBvKbYsymZo+1tel\nDXuWZVHYVMzmyjz21O7HbXnws/sxN34WuYk5ZESmanQvcpYGPPgNw3gYuBNo673rBiDANM35vR8I\nfgncYJrmI73Pfw6IAf7bMIzXTNN8ZaBrFLkQXq/FlgNVrH6/mKY2FxGhAdy+JI3caeNx2NWA50K0\nd3eQV72TzRV51HTUAhAfEktuQg4Xj59NmH+ojysUGX4GY8RfCNwIPN97OxfYAGCaZp5hGHNOfLJp\nmncPQk0i/aKg9BgvbiykvLaNAD87y+ensjQnmeDAIXcWbdiwLIvi5jI2V25jV+0+3F43fjYHc+Jn\nkJuQQ2ZUukb3IhdgwP86mab5qmEYqSfcFQ60nHDbYxiG3TRN7/luIzY2/HxfKnJKn7VPHalp5Zm1\nH7O9twHP4jkT+OLSLGK0ct55a3d18H5pHm8XfcCRlioAxofFsSQjl4VplxAROLznSOjvlAwVvhiW\ntNAT/sddUOgD1NW1XlhFIieIjQ0/7T7V0uHi9c0lvLe7Eq9lYUyI4rYlmaSOi8DqdmtfPEeWZVHa\nUs7mijx21u6l29uNw+ZgVtw0chPmMWlMBjabDWeLRR3D99/2TPuUyPm4kA+Svgj+LcB1wMuGYcxD\nl+vJMNDt9vD3HUdZt7WUTqeH+OgQbl2YwYyJMTrsfB463Z1sr97N5so8Ktp6RvcxQdEsSMzhkvFz\nCQ8Y3qN7kaFsMIPf6v26GrjSMIwtvbfvHcQaRM6JZVnkFdTwyrvFNLR0ERbszxeuSGfhzET8HJq4\ndy4sy6K89SibK7axo2YPLm83dpudGbEXkZuYgzEmE7tN/6YiA81mWdZnP2tos3QITfrT8cOyh482\n8Zd3CimpasHPYeOK2RNYPj+FkCB/X5c4rHS5u9hes4ctFds40lYJwNigMcxPyOGS8XOIDIzwcYUD\nT4f6pb/Fxoaf96FGTT0W+YSq+nZ+v3o/O806AOZOjuOmhRnEaeLeOekZ3eexo2Y3To8Lu83O9Jgp\nLEicR1b0RI3uRXxEwS/Sq72rmzVbStm46yhuj0VGQgS3LZ5IZlKkr0sbNpweFztqdrO5Io/y1qMA\njAmM4orky5mfcDFRgfq3FPE1Bb+Mem6Pl427KlizpYT2Ljdx0SHceGkacyfHaeLeWTraWsnmyjy2\nV++iy+PEho2pY7PITcxhytjJGt2LDCEKfhm1LMti16F6Xn63kNrGToID/bh1USa3XT2Z5qYOX5c3\n5Lk8LnbW7mNLxTZKWsoBiAyIYNGES1mQcDFjgqJ8XKGInIqCX0alkqoWXnznMIeONuOw21gyO4nr\nF6QSHhJAgJbLPSWXp5vy1qMUN5VS3FLK4cYSujxd2LCRPdYgN2EeU8dOxmHXv5/IUKbgl1GlobmL\nV94vYtvHPR33ZmTGcMuiDMaPVc/3T2pyNlPcXEZxcynFzWUcaa3Aa/2j19bYoDEsTJrP/ISLGRsc\n7cNKReRcKPhlVOh0unljWxlvbT9Ct9tLSnw4ty3OZHLKGF+XNiR4vB4q26spai6lpLmM4uYyjnU1\n9j3usDlIDk8iPTKFtMgU0iNTNFFPZJhS8MuI5vF6eX9vFa9/UExLRzdjwgO58bJ0Lpk6DvsonrjX\n0d1JSUtZ74i+jNKWclweV9/jYf6hXBSTTXpkCumRqSSHJxHgUP8CkZFAwS8jkmVZ7C9u4KVNRVTW\ntxPo7+Bzl6Zx1cXJBI6yc/iWZVHXWU9RcxklvYftq9prTnrO+ND43tF8KumRKcQFqxWxyEil4JcR\n50htGy9uPEx+aSM2G1w2PYHPXZpGZFigr0sbFH2T8HpDvqS5jLbu9r7HAxwBTBqT2TuaTyEtIpkQ\n/xAfViwig0nBLyNGU5uT1e8Xs3lfFRYwNS2aWxdlkhQ3shd8OT4J7/i5+SOtFXgsT9/j0UFjmBM9\nse/cfGLoeM28FxnFFPwy7DldHt78qJz1eeU4uz0kxoRy6+JMLkof6+vS+t3xSXjHZ9uXNJfRcMIk\nPLvNzoTwxL5z85qEJyKfpOCXYctrWXy4v5pX3y+iqc1FRIg/ty3J5NJp43HYR0anuJ5JeOV95+ZL\nW8pxnjAJL9Q/hItiskiPSCUtMoWUiCQCHAE+rFhEhjoFvwxLBaXHeHFjIeW1bfj72Vk+P4WlOSkE\nBw7fXfr4JLzjM+1LeifhWfxjBc1xofGkR/Qcsk+PStUkPBE5Z8P3r6SMSlUN7by0sZC9RQ0AXDJl\nHDddnk50RJCPKzt33Z5uylqP9p2bL24uPXkSnt2fiVHpJ107r0l4InKhFPwyLLR0uHh9cwnv7a7E\na1kYE6K4bUkmqeOGz1ruzc6Wk87Nl39iEt6YwChmx03vOzefGKZJeCLS/xT8MqR1uz28veMoa7eW\n0un0ED8mmFsXZTJj4tA+xO21vFS0Vfedmy9uLqOh61jf43abnQlhiSeN5rWojYgMBgW/DEmWZfFR\nQS1/fbeIhpYuQoP8uP2KiSyamYifY+hN3Ot0d1LSXN53br6kpezkSXh+IUwdm9U3216T8ETEVxT8\nMuQUHm3mLxsPU1zZgp/DxtUXT2D5/FRCg4ZGy9ieSXgNvefm/9EJ76RJeCFxfZ3wMiJTiAuJHdJH\nKERk9FDwy5BR29jBX98tYodZB8CcyXHcvDCDuKhgn9bV7emmvLWi79x8cXMZrd1tfY8H2P3JjErr\nOzefFplCqCbhicgQpeAXn2vv6mbNllLe2XkUj9ciIyGC2xZPJDPJN41nmjqb2VP7cd9EPE3CE5GR\nRMEvPuP2eNm0q4K/bSmhvctNTGQQNy/MYO7kuEE7LO61vFS2VZ+07vwnJ+ElhSWc1AlPk/BEZDhT\n8ItP5Jce4/m3DlFzrIPgQAe3LMrgitlJ+PsN7Mi5091JafORvpAvbSmny+PsezzEL5hZ46eSFNyz\n9nxyxAQCNQlPREYQBb8MqpYOFy++U8jWj6ux2WDxrESuz00jIqT/w9WyLOo7j/WG/Kkn4cWHxDHr\nxOVoQ2KIj4ukrq613+sRERkKFPwyKCzL4sMD1by4sZC2zm5SxoVzzzWTSRkX3m/b6PZ0c6Stouew\nfVPppybh+X9iEl5qZDJh/qH9tn0RkeFAwS8DruZYB6veNCkoayTQ38Hnl0xkyezEC15Ip9nZSknL\nP0L+SOtR3CdMwosKjGRW3LS+oE8KS9AkPBEZ9RT8MmDcHi8b8sr525ZS3B4v0zLGcudVk4iJPPfL\n87yWl6r2Gop6Q76kuZT6T03CG98X8umRqZqEJyJyCgp+GRCFFc08t+EgFXXtRIYG8IUrJzHHOPsm\nNp3uLkpbyvtG86eahDdl7OS+oE/RJDwRkbOi4Jd+1dHl5pX3inh3dwUWsHBGAjcvzCDkDF33LMui\noetYz2i+paflbWVb9Scm4cUys28039MJz24beq17RUSGOgW/9AvLsthp1vHHtw/R3OZi/NgQvnj1\nJJLHB9HpbqOxrYtOdxdd7i66PE463V20d3dQ3nqU4uZSWl2fnoR3fPGatIgUwgI0CU9EpD8o+OVT\nvJYXp8dFl7s3rD3OE74/HuDOvu9bujoorW2kzdmBLc1NVLBFh93N44UuKPzs7WkSnojI4FHwjyCW\nZeHydtPp7jwpmLvcvcHt6fpHgH/y8RO+d3qcJx1mPyv+YPezE+IfRKh/MEF+QQQ5Agn2CybIL5Ag\nRxDBfkEE+QX2fO29nRA2jjGBUVrARkRkkCj4hwDLsnB73WcZzKcahTv7DqF7Le85b99us/eF8djg\nMX3f94R0cG+AB50Q5kG0tHp5a1sVFTUugv0Cufkyg8unTVCAi4gMcQr+C+T2uulyO/sOg58umP9x\nfvtUo3DnSYvAnC0btr4wjgqM7PneL5BgR09InxzgJwd30Amjbn+731kHtrPbw+ubS3jroyN4LT/m\nTUnk84snEhGqGfUiIsPBqA1+j9eD83ggnxjMJ4Txid+fMrg9XXR73ee1/SBHIEF+QYQHhBMXHNMb\n2kEE995/yuB2BBHsF9gX2oGOgEEdYR8obmDVmyb1zV3ERAZx19UGU9PHDtr2RUTkwg374C9sKOVo\nQz1dvee1jx8u/9T3nxiRu7zd57W9AEcAwY5AQvyDiQ4e0zu6DuwN7aBPjbp7Avzkc9yBjsBhdSla\nS7uLv7xzmG35NdhtNpbmJHN9bhqB/pqAJyIy3Az74P/22z89q+f52/36wjcyMOKEID5FUH9q1N0z\n4g50BI6q2eaWZbF5XxUvbSqkvctN2vhw7r5mMsnx/ddfX0REBtewD/6bpyzD1enpC/HjQf3JQ+V+\n9mH/ow6q6mMdrNpwkIPlTQQGOLj9ioksmZWE3a7JeyIiw9mwT8Nbpy7XEqr9yO3x8sa2MtZ+WIbb\n42VGZgx3XjWJ6IggX5cmIiL9YNgHv/SfQ0eaeG7DQaoaOogKC+COKycxa9LZ99cXEZGhb8gFv2EY\ns4GvATbgYdM0a31c0ojX0dXNy+8W8d6eSmzAolmJ3HRZBiFBQ273EBGRCzQU/7IHAv8KXAVcArzu\n23JGLsuy2GHW8ae/H6K53UVibCh3XzOZzMRIX5cmIiIDZMgFv2maHxqGcQnwLeBWX9czUtU3d/LC\nW4fYV9SAn8POTZenc/XFyfg5hs9lhiIicu4GJfgNw8gBfmKa5iLDMOzAE8A0wAmsME2zyDCM/wIm\nAr8CdgBLgUeBBwejxtHC4/Xyzo6jrP6gBGe3h6yUMdx1tUF8dIivSxMRkUEw4MFvGMbDwJ3A8XVX\nbwACTNOc3/uB4JfADaZpfrf3+YuApwEX8PuBrm80Katu5dn1BymraSUs2J87r5rE/KnjNHlPRGQU\nGYwRfyFwI/B87+1cYAOAaZp5hmHMOfHJpmluAjYNQl2jhtPl4bXNxby1/QiWBfOnjuO2xZmEh6i/\nvojIaDPgwW+a5quGYaSecFc40HLCbY9hGHbTNM99WTn5TPuK6nn+zUM0tHQRFxXMF68xmJIa7euy\nRETER3wxua+FnvA/7oJDPzZWLWQ/qbGliyf///buPbjK+s7j+DsJIVySIJeA3MFAvqFaREW8FCm3\nQLBrh5baiiLR7rZruzvb3bF1dtxept1Zt+tWZ2d2tjPdznQLkVrX21JhCSAgykUsAkWsfGO4iFyE\nUJAQCCGXZ/94nnRTFEngJOecPJ/XP3CePOfkew6/4fP8nuc539/SXby24xBZmRncPWMsXykx9ddv\nI40pSTSNKUkVyQj+jcBdwLNmdiuw80pfUJ37/l9zEPDa7w7z7Lo9nK1vpHBIPmWlxQwbmEvNh2eT\nXV5aKCjI05iShNKYkkS7kgPJzgz+IPrzRaDEzDZGjx/sxBq6tMPHz7C4YjeVB0/Ro3sWC2YVMXXC\nUPXXFxGRP8oIguDSe6W2IO5H0g2NzSzfvJ/lm9+jqTngpqIC7i0pom9eTrJLS0uanUmiaUxJohUU\n5O70rcEAAAv/SURBVF32jC7lGvhI+/iBkyyqcD44cZa+eTksKCnihqKCZJclIiIpSsGfpmrrGnh2\nXRWv7TxCBjDjpmF8cco19MzRP6mIiFycUiLNBEHAG+8c4+mXK6k528CwglzK5hiFQ9RfX0RELk3B\nn0aqP6yjfJWza+8JsrtlcvfUQkpuHq7++iIi0mYK/jTQ1NzM6t8e5H9e28v5xmauHdWX+2cbA/uq\nv76IiLSPgj/F7TtSw6IVuzlwrJbcntmUzSnm1k8NUn99ERG5LAr+FHXufCMvvrqPl98M++tP/vRg\nvjx9DLk9s5NdmoiIpDEFfwra8e5xnlrtnKipZ1DfniwsLWbcyL7JLktERLoABX8K+bC2nl+trmSr\nV5OVmcGf3T6Ku24fSXY39dcXEZHEUPCngOYgYP2Owzz3ShV19U2MGdqHslJjaEFusksTEZEuRsGf\nZIeqa1lU4VQdOkXPnG4snG1MmTCETN28JyIiHUDBnyQNjU28tOk9Vrwe9tefWDyQe2eO5apc9dcX\nEZGOo+BPgnfeO8niit0cPVlHv/wcFpQYE8YOSHZZIiISAwr+TlRb18Aza99l41sfkJEBJROH84Up\no+nRXf8MIiLSOZQ4nSAIAl5/+yhPr3mX2roGRgzMpWxOMaMH5ye7NBERiRkFfwc7dvIs5Sudt/ef\npHt2Jl+eNoaSm4eRlan++iIi0vkU/B2ksamZVb99n6Ub9tHQ2Mx11/Tj/llGwVU9k12aiIjEmIK/\nA+w5fIpFK5yD1bXk98rmq3eOY9K4geqvLyIiSafgT6C6+kZeeHUva988SABMuX4wX5qq/voiIpI6\nFPwJsq2ymiWrKzl5up6r+/WirNSwEeqvLyIiqUXBf4VOnq5nyepKtlVW0y0rg89/ZhSfu20U2d10\n856IiKQeBf9lam4OWLf9EM+v38O5800UDevDwtJihgzonezSRERELkrBfxkOHqtlUcVu9hyuoVdO\nNx6YU8zk8YPVX19ERFKegr8dzjc08dKm/VRsOUBTc8CkcQOZP2MsfdRfX0RE0oSCv43e3n+C8grn\n2Id19M/vwf2zjfGF/ZNdloiISLso+C+h5ux5nllTxea3w/76sycNZ+7ka8jpnpXs0kRERNpNwX8R\nQRCwadcHPLO2itq6BkZenccDpcWMvDov2aWJiIhcNgX/xzh64iyLVzrvvHeSnOws7pk+hhkT1V9f\nRETSn4K/lcamZiq2HOA3G/fT2NTM+ML+LJhVxIA+6q8vIiJdg4I/UnXoFItW7ObQ8TPk9+7OfSVF\nTLQC9dcXEZEuJfbBf/ZcI8+v38Mr2w8RAFMnDOFLUwvp1UP99UVEpOuJbfAHQcCbXs2Slys5VXue\nwf17UVZaTNHwq5JdmoiISIeJZfCfqDnHU6sq2VF1nG5ZGcy9YzRzbhmp/voiItLlxSr4m5sD1mw7\nyAuv7qX+fBM2/CoWlhqD+6u/voiIxENsgv/A0dMsqtjNviOn6d2jG/dG/fV1856IiMRJlw/++oYm\nlm7Yx6o33qc5CLj12kHcM30s+b27J7s0ERGRTtelg3/X3j+weKVz/NQ5BvTpwcLZxnXXqL++iIjE\nV5cM/poz5/n1mnd5/fdHyczIYM4tI/j85NHkZKu/voiIxFuXCv4gCNiw8wj/va6KM+caGT04j7LS\nYkYMUn99ERER6ELBf+QPZyhf6ew+8CE53bOYP3MsM24cRmambt4TERFpkZLBb2aDgGXufvOl9m1o\nbOY3G/exbNN+GpsCJowZwIJZRfTL79EJlYqIiKSXlAx+4DvA/rbs+K0nX+H9o6fpk9udBSVF3Fik\n/voiIiIXk3LBb2bfAJ4CHm7L/gePnWbajUOZN6WQXj1S7u2IiIiklIwgCDr8l5jZLcCP3X2amWUC\nPwXGA/XAX7j7HjP7ETAWGAhUAtOBR939+U967Zoz54P6s/Ud+wYkVgoK8qiuPp3sMqQL0ZiSRCso\nyLvsU9sdPkU2s0eABUBttGku0N3db48OCJ4A5rr79y943uJLhT5Afu/uVCv4RURE2qQzVqWpAr4I\ntBydTAYqANx9CzDx457k7gs7oTYREZFY6fDgd/cXgMZWm/KAmlaPm6LT/yIiItLBknE3XA1h+LfI\ndPfmK3i9jIICNeiRxNKYkkTTmJJUkYyZ9kbgTgAzuxXYmYQaREREYqkzZ/wtXx94ESgxs43R4wc7\nsQYREZFY65Sv84mIiEhq0E11IiIiMaLgFxERiREFv4iISIx0yeb2ZjYdmO/uX0t2LZL+zGwG8BWg\nF/C4u+ubKHJFzOwm4K8JG5s94u7HklySpLn2rGrb5Wb8ZlYITAC0Lq8kSk93/zrwE2BWsouRLiEH\n+FtgOXBbkmuRNGdmGbRjVdsuF/zuvsfdn0x2HdJ1uPsyM+sN/A3wyySXI12Au28CPgV8G9iR5HIk\n/T1EuKrtubbsnFan+tuyyl9SC5S008aVIwcAjwPfd/fjSSxX0kAbx9TNwFZgDvAD4FtJK1hSWhtz\nb2a0bZKZzbvUAndpM+OPVvn7OeEpMmi1yh/w94Sr/Im0WTvG1BPAIOCfzWxepxcqaaMdYyoX+AXw\nr8CSzq5T0kNbx5O7z3P3bwBb2rKqbTrN+FtW+SuPHv/JKn9m9ier/Ln7/Z1bnqShNo0pdy9LTnmS\nhto6ptYB65JSoaST9uZem1a1TZsZv1b5k0TTmJJE05iSROqo8ZTOAzDRq/yJaExJomlMSSIlZDyl\nc/BrlT9JNI0pSTSNKUmkhIyndLrG30Kr/EmiaUxJomlMSSIldDxpdT4REZEYSedT/SIiItJOCn4R\nEZEYUfCLiIjEiIJfREQkRhT8IiIiMaLgFxERiREFv4iISIwo+EVERGJEwS+SJsxslJk1m9nMC7bv\nN7MR7XiNtzqmwrYzs/8ys+HR35eb2dXJrkkkLhT8IumlAfi5meW22paO7TenEv3/4+6fc/cPkluO\nSHykY69+kTg7DKwCngD+8pN2NLNHgfuApug5j0Q/yjWzF4BCoBL4c3evMbOfADOj/Ze6+4+iA4z/\nAK4FsoB/cfdfm9kDQBnQH9gAfAEY7u6NZnYdsMTdrzezfwKmA/2A44Rriz8IDAGWm9kUYBswBTgI\n/Fu0fwCUu/vjZjYVeBQ4A4wD3gLuBXoCTwODovf1Q3d/qZ2fp0jsaMYvkn6+Dcy+8JR/a2Z2J3AX\ncCNwAzAGeCj68TDgMXe/HtgHfDe6VFDq7hOA24ExZpYDfBfY6u4Tgc8C/2Bmo6PXGQpMcPdvAluA\n2dH2+UC5mRUC5u63ubsBVcB97v5jwgOYO939BGHIZ0T1DQU+DUwC5kXvA+A24K8Ig39E9LvmAvui\n2hYAd7T3gxSJIwW/SJpx99PA1/joKf/WpgG/cvd6d28CfgHMIAzZt9x9a7RfebT9EFBnZhuAvwO+\n5+71hGcAHjKz7cB6oBfh7D8AtrVaC7wcuCf6+93R794DPGxmXzezJwjDu/cnvLVpwC/dPXD3OmBJ\nq5p3ufthdw+Ad4C+wCZgrpm9CEwG/vHSn56IKPhF0pC7rwZWA09eZJdMwll068ctl/YaL9jeGB0c\n3AJ8j/D0/WYzGxv9/D53v8HdbwA+A6yMnlvX6nWWAZ81szuA9939sJndRHiJAeBZwiVFW9fUnprP\ntdoeAJnuXgUUEx4g3AG88QmvLSIRBb9I+noYmEV4vfxCa4H5ZtbDzLoRXldfSxis15vZtdF+XwVW\nm9l4whn9q+7+HeD3gEXP+SaAmQ0GtgPDuSDAo7MDFYTX6MujzVOAV9z9Pwln6bMI7xOA8OAj+2Nq\nLjOzTDPrRXgdv6XmjzCzhwiv6z9HeBlgoJnlX+SzEpGIgl8kvfzxDv5Wp/w/cpOuuy8nnIVvBXYR\nXsv/95YfA4+Z2U7Cm+4ec/edwGZgl5m9Ge3/v8APgZ7RVwDXAI+4+96ojgu/TVBOOAN/Lnr8DOFB\nxvZo2wqg5f6AZYQ3941q9b5+RniD3+8Ib/hb6u5LL3zfrR4vASx6H+uBH7h7zcd/bCLSIiMI0vGb\nQCIiInI5NOMXERGJEQW/iIhIjCj4RUREYkTBLyIiEiMKfhERkRhR8IuIiMSIgl9ERCRGFPwiIiIx\n8n+CsGDjIc8YgAAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 10 }, { "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": {} } ] }