{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Implementing the NFFT" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import division\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to solve the following:\n", "\n", "$$\n", "\\hat{f}_k = \\sum_{j=0}^{M-1} f_j e^{2\\pi i k x_j}, \n", "$$\n", "\n", "for complex values $\\{f_j\\}$ at points $\\{x_j\\}$ satisfying $-1/2 \\le x_j < 1/2$,\n", "and for integer wavenumbers $k$ in the range $-N/2 \\le k < N$.\n", "\n", "A straightforward implementation of this sum would require $\\mathcal{O}[MN]$ operations, but the nonequispaced fast Fourier transform (NFFT) allows this to be computed in $\\mathcal{O}[M\\log N]$.\n", "In this post, we'll work on writing a pure Python implementation of this NFFT algorithm." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Straightforward implementation\n", "\n", "As our first step, we'll define a simple reference implementation of the nonuniform direct Fourier transform (*NDFT*):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def ndft(x, f, N):\n", " \"\"\"nonuniform discrete Fourier transform\"\"\"\n", " k = -(N // 2) + np.arange(N)\n", " return np.dot(f, np.exp(2j * np.pi * k * x[:, np.newaxis]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try evaluating this on some sinusoidal data, with a frequency of 10 cycles per unit time:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOWh//HPkz0TIJCEJSRAwiICQWUVRStVUbQq1qVq\ne11arb9bu7tVa1vp4r0tvdr29t7aa62Ver1V64pWEbXizo7sIGEPJEASyDZZZ57fH+dMMkASskxm\nwsz3/XrNK2fOOTPz5GQy3/M8z3meMdZaREQktsVFugAiIhJ5CgMREVEYiIiIwkBERFAYiIgICgMR\nEUFhICIiKAxERASFgYiIAAmRLkBHZWVl2by8vEgXQ0TkpLJq1apSa+3AE+130oRBXl4eK1eujHQx\nREROKsaY3R3ZT81EIiKiMBAREYWBiIhwEvUZiIh0RmNjI0VFRdTV1UW6KGGRkpJCbm4uiYmJXXq8\nwkBEolJRURF9+/YlLy8PY0yki9OjrLWUlZVRVFREfn5+l55DzUQiEpXq6urIzMyM+iAAMMaQmZnZ\nrVqQwkBEolYsBEFAd39XhYGIRF7xWtizLNKliGkKAxGJvLfnwet3RboUvcott9zC888/H7bXUxiI\nSORVH3JuUcpai9/vj3Qx2qUwEJHI85aCtwysjXRJQmbXrl2MHTuWm266iYKCAp566inOOussJk+e\nzLXXXkt1dTUAP/vZz5g2bRoFBQXcfvvt2AgdA11aKiKRZS3UlIK/EeqrIKVfyF/ip69uZNP+ypA+\n5/ih/Xjw8gnt7rNt2zYWLFjA6NGjueqqq3j77bdJS0vjV7/6FY888gg/+clP+Na3vsVPfvITAG68\n8UZee+01Lr/88pCWtSMUBiISWfWVThCAU0PogTCIlBEjRjBjxgxee+01Nm3axMyZMwFoaGjgrLPO\nAuDdd99l/vz5eL1eysvLmTBhgsJARGJQTWnLsrccMkaG/CVOdAbfU9LS0gCnz2D27Nn87W9/O2p7\nXV0dd9xxBytXrmTYsGHMmzcvYiOm1WcgIpHlLWtZDg6GKDJjxgw++ugjCgsLAaipqeGzzz5r/uDP\nysqiuro6rFcPHUs1AxGJrKNqBmVt73cSGzhwIE8++SQ33HAD9fX1APziF7/glFNO4etf/zoFBQUM\nGTKEadOmRayMJlI91501depUqy+3EYlCq/8KC7/tLM/+Gcz8bkiedvPmzYwbNy4kz3WyaO13Nsas\nstZOPdFj1UwkIpEVqA2YuKitGZwM1EwkIpFVUwoJqZDaX2EQQQoDEYksbxmkZUFKf6hRGESKwkBE\nIqumFDyZkJKumkEEqc9ARCLLW+rUDNKynGWJCIWBiERWTRl4spzagWoGEaMwEJHICtQMPJlQVwG+\nxkiXKGTOPvvsSBehwxQGIhI5DV5o9DpB4Ml01nnLI1umEPr4448jXYQOUxiISOQEmoUCNYPgdVGg\nT58+ACxZsoTzzjuPuXPnMnLkSO677z6efvpppk+fzsSJE9m+fTsAr776KmeeeSaTJk3iwgsv5MCB\nAwAcOnSI2bNnM2HCBG677TZGjBhBaWlo+1d0NZGIRE6gw9iTCcl9j14XSm/cByXrQ/ucQybCJb/s\n8O5r165l8+bNZGRkMHLkSG677TaWL1/O7373O37/+9/z29/+lnPOOYelS5dijOHxxx9n/vz5PPzw\nw/z0pz/l/PPP5/7772fRokX8+c9/Du3vQgjDwBgTD6wE9llrLzPGZADPAnnALuBL1trD7r73A7cC\nPuA71to3Q1UOETmJBMYVeLIg2TmLjqaaQbBp06aRnZ0NwKhRo7jooosAmDhxIu+++y4ARUVFXHfd\ndRQXF9PQ0EB+fj4AH374IS+99BIAc+bMYcCAASEvXyhrBt8FNgOBycjvA96x1v7SGHOfe/8Hxpjx\nwPXABGAo8LYx5hRrrS+EZRGRk0GgFpCWBUluGPTEzKWdOIPvKcnJyc3LcXFxzffj4uJoamoC4Nvf\n/jZ33nknV1xxBUuWLGHevHlhK19I+gyMMbnAF4DHg1bPBRa4ywuAK4PWP2OtrbfW7gQKgemhKIeI\nnGRqgpqJPBnOchR1IHdWRUUFOTk5ACxYsKB5/cyZM3nuuecAWLx4MYcPHw75a4eqA/m3wL1A8Dc+\nD7bWFrvLJcBgdzkH2Bu0X5G77jjGmNuNMSuNMSsPHYreL8sWiVneUohLdEYfxydCcnpMDzybN28e\n1157LVOmTCErK6t5/YMPPsjixYspKCjg73//O0OGDKFv374hfe1uNxMZYy4DDlprVxljZrW2j7XW\nGmM6PVe2tfYx4DFwprDuVkFFpPcJTEVhjHM/LboGngW+9H7WrFnMmjWref2SJUual4O3zZ07l7lz\n5x73POnp6bz55pskJCTwySefsGLFiqOanUIhFH0GM4ErjDGXAilAP2PM/wIHjDHZ1tpiY0w2cNDd\nfx8wLOjxue46EYk1gUnqAjQKuVV79uzhS1/6En6/n6SkJP70pz+F/DW6HQbW2vuB+wHcmsHd1tp/\nMcb8GrgZ+KX78xX3IQuB/zPGPILTgTwGWN7dcojISchb1jK+AJyriiqKIleeXmrMmDGsWbOmR1+j\nJwed/RKYbYzZBlzo3sdauxF4DtgELAK+qSuJRGJUoJkoIMQ1g5PlmxxDobu/a0gHnVlrlwBL3OUy\n4II29nsIeCiUry0iJ6HAvEQBaZnOOmtb+hG6KCUlhbKyMjIzMzHdfK7ezlpLWVkZKSkpXX4OjUAW\nkcjwNToT03mO6TPwNUBDdcuI5C7Kzc2lqKiIWLkSMSUlhdzc3C4/XmEgIpHRPC/RMc1EgW3dDIPE\nxMTmEbxyYpqoTkQio3nAWXDNwF3W11+GncJARCIjeCqKgCicufRkoTAQkchorWYQaDKK4VHIkaIw\nEJHICMxBdOylpaCaQQQoDEQkMrylgGmZoA4guZ8zV1FPzFwq7VIYiEhk1JRC6gCIi29ZZ4ympIgQ\nhYGIRMaxA84CPJkxPY11pCgMRCQyasqO7jwOCIxClrBSGIhIZHhLjx5wFqBmoohQGIhIZNSUtl4z\n8GSpAzkCFAYiEn5+P9SWt91nUHcEfE3hL1cMUxiISPjVHQHrP3qMQUAgIGrViRxOCgMRCb/WRh8H\nBMYdqN8grBQGIhJ+zfMStdGBDAqDMFMYiEj4tVszyDp6HwkLhYGIhF9rM5YGqGYQEQoDEQm/wPcV\ntNaBrDCICIWBiISft9SZlC4h+fhtCUnONoVBWCkMRCT8akpbrxUEeDLVZxBmCgMRCT9v2YnDQDWD\nsFIYiEj4tTVjaYDCIOwUBiISfm3NWBqQlqUwCDOFgYiEl7Vtz1ga4MlwwsDa8JUrxikMRCS86qvA\n19B+zcCTBU110FATvnLFOIWBiIRXewPOAjTWIOwUBiISXs0Dzk7QZwD6xrMwUhiISHgFPuBPdGkp\n6LuQw6jbYWCMGWaMedcYs8kYs9EY8113fYYx5i1jzDb354Cgx9xvjCk0xmw1xlzc3TKIyEkk0PTT\nbgeymonCLRQ1gybgLmvteGAG8E1jzHjgPuAda+0Y4B33Pu6264EJwBzgD8aY+BCUQ0ROBu3NWBoQ\nCAONQg6bboeBtbbYWrvaXa4CNgM5wFxggbvbAuBKd3ku8Iy1tt5auxMoBKZ3txwicpLwlkJCCiSl\ntb1PSjrEJahmEEYh7TMwxuQBk4BlwGBrbbG7qQQY7C7nAHuDHlbkrhORWBAYcGZM2/sY445CVs0g\nXEIWBsaYPsALwPestZXB26y1Fuj06BFjzO3GmJXGmJWHDh0KUUlFJKJONOAswJOpDuQwCkkYGGMS\ncYLgaWvti+7qA8aYbHd7NnDQXb8PGBb08Fx33XGstY9Za6daa6cOHDgwFEUVkUirKW2/vyBA8xOF\nVSiuJjLAn4HN1tpHgjYtBG52l28GXglaf70xJtkYkw+MAZZ3txwicpI40SR1AZrGOqwSQvAcM4Eb\ngfXGmE/ddT8Efgk8Z4y5FdgNfAnAWrvRGPMcsAnnSqRvWmt9ISiHiJwMak4wfXWAagZh1e0wsNZ+\nCLTVE3RBG495CHiou68tIieZxlporOlYGKRlQe1h8PsgTlef9zSNQBaR8GkecNbBZiKsEwjS4xQG\nIhI+HRlwFqCBZ2GlMBCR8OnIjKUBmpIirBQGIhI+HZmxNKB55lKFQTgoDEQkfJprBh28mij4MdKj\nFAYiEj41pWDiIaX/ifdVM1FYKQxEJHy8pc6HfHvzEgUkJENS35amJelRCgMRCR9vecc6jwM8GaoZ\nhInCQETCp6a0YwPOAtKy1GcQJgoDEQmfjs5LFKApKcJGYSAi4dPRGUsDPFmaxjpMFAYiEh6+Rqg7\n0vk+A41ADguFgYiER+AMv7N9Bk210ODtmTJJM4WBiIRHoCO4M2GggWdhozAQkfCo6cS8RAEaeBY2\nCgMRCQ9vJ+YlCgjsq4FnPU5hICLh0ZnvMghQzSBsFAYiEh6BZqLUjI4/Jk1hEC4KAxEJD28ppA6A\n+E58225yujOxnTqQe5zCQETCo7MDzgDi4jQKOUwUBiISHt6yzvUXBHgyNfAsDBQGIhIenZ2kLsCT\nqSkpwkBhICLh4e1iGKRlqs8gDBQGIiH05sYSFm0oiXQxeh+/v/PfZRCgPoOw6ES3voi0xVrLo+9t\nZ/6irQDcO2cs3zhvFKYj3+gVC+qOgPV1vgMZnMfUHga/D+LiQ182ARQGIt3m91t+8Y/NPPHRTq44\nfSjGwPxFWymtauBHXxhHXJwCoUsDzgI8mWD9UHukZdyBhJzCQE4qfr/lg8JSBngSOS23A1+q3sMa\nmvzc8/xaXvl0P1+dmcePvzAegMy0ZJ74aCdlNfX8+prTSUqIfIvsuqIjHPY2cu7orPAHVE0XJqkL\nCASIt0xh0IMUBnJSqKlv4vlVRTz58S52ltaQEGf46dwJfOXMESF/LWtth5p3auqb+MbTq3n/s0Pc\nc/FY7pjV0iz048vGkdU3ifmLtnLY28ijX5lMWvKJ/906+tqd9fSy3Tz4ykaa/Jb8rDRuOTuPa6bk\ndqhMIeHtwiR1AZ6MoOc4JWRF6jWqSiBtYMSbwBQGvVhReTWZaSmkhusfthfaW+5lwce7eHbFXqrq\nmzhjWH9+c93pLPx0Pw+8tIEtxVX85PLxJMZ3/szbWsve8lo2FVeyqbiSzcWVbNpfyaHqes4ZncUl\nBUOYPX4w/T1Jxz22vKaBrz65gvVFR/jV1RO5btrwo7YbY7hj1mgy05K4/8X1fPnxZfzllmlkpB3/\nXEe8Dby16QBvbCjhw8JSBvZJZvzQfozL7sd49zYsI7VLIdHo8/OzVzfx1NLdzBo7kLlnDGXBx7t5\ncOFG/uPNrVw3bRg3n53HsAxPp5+7U5prBl1sJoIe60Sua/RRWl1P7oAePgbHqj4Er98Fm16BpD6Q\nMxlyp0HudMid2rXg7IaIfcoYY+YAvwPigcettb+MVFl6HV8ja/7+S0Zv/m8qjIftg2YyfPpl9Bs/\nu+UsqT0NNXBgk9PpNurzEJ/Y82UOIWsty3aW85ePdvLWpgPEGcMlE7P56sw8Jg8fAMAVp+cw/80t\n/M97O9h2sIo/fGVKqx+0x9p2oIpnVuxlXdERthRXUVXfBECcgfysNCYN788ATxL/3HKQf245SEKc\n4ezRWVxaMISLJgwhIy2JfUdqufHPy9h3uJY//ssULpowpM3Xu27acAZ4kvj239ZwzR8/5qlbzySn\nfyrlNQ0s3ljC6xtK+LiwlCa/Jad/KtdNHcZhbwObiyt5Z/MB/NZ5nr7JCZya3ZeJOf25Yfowxgzu\ne8LftbymgTueXsXSHeX8v/NGcu+sHOJ99Xxx0kxW7znMXz7axZMf7+KJj3Yye/xgvjoznzPzM0Jb\nM6kpgx3vwrpnnftdGmcQmLk0tJeXllbX89Qnu3lq6W4qahv5xZUF3DB9+IkfGAobX4J/3AX1VTDz\nu86X9xQthw9/63S0A2SMdMNhGpzxZUhK69EiGWttj75Aqy9qTDzwGTAbKAJWADdYaze19ZipU6fa\nlStXhrQc1fVNPLdiL0P7p/D5UweRnBD5KxX825dQ/vfvklW3i7UpU2mM8zCmZiXpxosfQ8Og00gZ\nOxtGne+8SaoPwIENULIBDqx3fpbvANy/a/owOPs7MPlGbEIKu8q8LN9ZxvKdhyk67CUh3hAfF0dC\nnHFuQfdPy03nhunDSUkM33Gprm/i+89+ylubDjDAk8iXzxzOjTPyGJKe0ur+L60p4gcvrGdwv2Qe\nv2kaYz3VsOxRyBwNp38Z4hOaw+Wx93fwzy0HSUqI47ScdMZlu2ffQ/sxdnBfUpNafk9rLeuKKnhj\nQwmvry9mT7mX+DjDjJEZbD9YQ01DE3++eRrT8zs26dryneXcumAFaUkJjBqUxtId5fj8luEZHi6d\nmM0lBUM4LTe95YPY76fu4HaKt62mas9azKHN9K/aRmpTJQuaZrNt1M3cMmtCmx/eW0uquO2vKzhQ\nWc/8L57KlY2LYMm/Q10FZJ8OYy6CMRdR0mcCTy3fy/8t28NhbyOzxw/mN9edQZ+u1kZ9jVC0Agrf\nge3vwP5PAQsp/eG06+DS+cc9pK7RxzPL97C2qIImv8Xn99PkszT5nVucr44n913O6wO/TtX07zA9\nP5O8TE+XQ2v7oWr+/OFOXlhVRH2TnwvHDaa+yccH20r51/NGce/FY3uuT6WmzKkNbHwJss+AKx+F\nweNbtjd4Yf8aKFqBb+9ymnYtI6HhCHE/LMIkpnbpJY0xq6y1U0+4X4TC4CxgnrX2Yvf+/QDW2n9v\n6zGhDANrLa+tK+YX/9jEgcp6ANJTE7n89GyumpzLpGH9O/ZGqymFra87/wBYsO6NoJ+pGZB3DqTn\ntP9cFfvwLfoh8ZtfZrd/EO+Puosb/uV2EuLj2H6ggjffeoPGrW9ztlnL5LhC4vFjTRzG+lueY0A+\nDCmAwRNhSAE+v5+6935D2oFVVMYPYIG9lMe8n6cKDxlpSYwe2Ae/tTQe8w/o81vqG33sr6gjOz2F\n71wwhmum5Ha4KWZzcSUfFZZy8YQhnWp+KDrs5bYFK9l2sJp7Lx7LzWfndSiIPt17hG8t+JirG17h\n24kLSfA5X5FoB57K8tHf56HPcli3r5LMtCRuOiuPG88a0aFaRIC1lk3FlbyxvoTXNxRjLfzhK5MZ\nl92vw88BznG54+nVGHACYOIQxg/wY8p3wuGdToiX74SDm+DQVmgM+qrHAXkwaDwNjY0k7XiLg2Tw\n68Zr2DbkMm793BguKRhCgvv3WbyxhO8/+ylpyQn83wW1jF71Czi0xTmBGHG280G9d5lzhU7qABh9\nIQ35F/C3sjH87N2DjBmYxuP/MpHcPnHQVA9Ndc6tsRYaqp2z2foqJ1jqq6C+kurKw5Tt30lu5Wri\nG6qcyeVyp8KoC2D0BTB00nFt4o0+P8+vKuI/39lGcUUdQ9NTSE6MJ/6YE5PEOMNTB69iuZ3AN+v+\nlWo8ZPVJZnr+AKblZTAtL4Nx2f2Ib+cD3FrLyt2Heez9Hby9+QCJ8XFcPTmX287NZ9TAPjT5/Mx7\ndSP/u3QPX5iYzcNfOt1575UWwuIfOX+fETNh5HmQd27HaujH2vQKvHanc9xm3Qczv3fcpH3WWtbs\nPcKLq4t4dW0xFbUNnNqnlr9+93IG9W39hOhEensYXAPMsdbe5t6/ETjTWvutth4TqjAoPFjNvIUb\n+bCwlIKcfvz0iglU1/t4cXURb24soa7RT35WGldNyuHKSTltf5htfNmp5nV0ZGTmaMj/HOS7b6bA\nVRFN9fDJf2Pf+zUNTU38d9MV9L/gbr4669TjAqm0up6/frKblz7eyLj6tczutxdvajbb4/PZbkZQ\n6U+h0eenoclPfZOfI94GahqamG62cGfqa8zwr6EhoS/Vp3+VAZ//DqbPwHaL/HFhKb9evJU1e46Q\nl+nh+7NP4fLThrZ61lThbWTh2n08t7KI9fsqAPAkxXP3Rc6HenyccQKyoQbqK6G+2v2QcT5stu0/\nxP+8s4k4fz23njmUsSOGOscrdUD7x9Va2PIPmhb9kISK3bzpm0rJmT9iWEMhp6z/D3JtCaviT+PA\njB9x/qwLO1fLaayDI7udD+jgD2tvKcQltHGLd37GJx59P3Az8U5tLvB8tYePfs0+g2HgqTBovHPG\nOGgCDBwLyX1a9tmzFP+iHxK3fxXb40Ywr+4Gdqafya3n5FNV18Qjb33Gxdlefjfg76TseNM5Sbj4\n32DsJRB4T9Uehu3vwra3oPAtqDkEGHzxScT76jt+jAC/SaDCn0q57cMy/6ns7n8W+dMv4ZKp40j3\nHN9E6fdbXl23n9+89Rm7yrxMGt6fey4ay9mj22kj/8ddsOJxfEn92Dzsep6J+wLv7vWz70gtAGlJ\n8fT3JJGcEEdSQhyJ8c7PJPdnaXU9G/dXMsCTyI1n5XHTWSPI6pN81EtYa3n8g5382xubOTM3hSfy\nl+BZ+SgkpjrBtne5E4YY56Qr/zznNuIsSG6n2c5bDq/fDRtecGplVz4Kgycctcveci8vr9nHi2v2\nsbO0hpTEOC6eMISrJucyc1Rmc9B3RVSEgTHmduB2gOHDh0/ZvXt3l1/T29DE7/9ZyOMf7CA1MZ57\nLh7Ll88ccdTZRFVdI29sKOHF1UUs3eHMhXJmfgZ3XzyWaXnumcCx1bzLHoF+OYABE+f+s5mWf7qK\nItj5vnPb/ZH7ZgKGTHTONArfhrJC3o8/k5/Wf4W7r7uISyZmt/u71Db4eH51ES+tLsLnt86bPvAP\nEPRP0DclgckjnLOnof1Tnernh7+BTQshIQXO/xGc3Wb+As4/yDubD/Ifi7eypaSKU4f05e6LxnLB\nuEFYC5/sKOPZFXtZtLGEhiYfFw88zNeydzKuYQMlxUX4aivJTKgjM6GB+MYq52y0o0w8DDsTxsx2\nmjUGT2g5rgAHt8Ci+5w26YGn0jD737hndQavfLofgBnD+/Dg0GWcuvUPmNojcPoNzu98bC2trgIO\nboYDG52fh7Y4H/qV+2hubgNI6gsZ+dBnkPN7+JucgVD+ppabr+no+34f+BuPvu/JdJ4nY6TzQd28\nnNfxdmFrYeNL2LfnYY7s5tOkydxX9SX22oH8dug7XFjxPCY+CT53N8y4AxKS234uvx+KP3Wadeqr\nKK+P4/m1hyitM1w6KY8z8rOdxyekOKGU3A+S+7G9Ko77/7GL5UVeZo0dxD0Xj2XFznKeXVnE5uJK\nkhLimDNhCNdNG8ZZIzMxhjbfSx2qie9fAx88Aptfdcoy5RZKCm5jWVkqa/YcobK2kQb3ZKjB56fR\n5yepoYKcukIG2VJGnH4eF517TvsXZFjLmkV/YcjSn5Ntyqkcey39Lv8352/ua4R9q93/5/ec2pWv\nwQn5lP4c1SIQ3ErQVOcsn/cDOOd7R/XhrdhVzsOLtzZ/3swYmcFVk3O5pGAIfVNC09fX28MgbM1E\n1lre3HiAn7+2iX1Harl6ci73X3rqcWcFx9pb7uWVT/fx9LI9FFfUcc2UXB4ctY2+b/+g3Wpeu3yN\nzht653uw4z3Yu5xaTzZ3Vn2ZFQmT+dNNU5k0/ARnwqFw6DNY+C3nw/T+PR16iN9veW19MY8s3squ\nMi+n5aZTVt1A9ZFDzE7ZxA0ZhZxWv4rEmmLnARmjsOk5FNcmsrKkiSO+FE4ZkcPUscNJSE2HpL74\nE5J5YV0ZL6wrZfTQTO699HT69e0LiSlQWeycsW5bDMVrnefsOxTGXAijL4TdH8PyPzkfTp9/AKbe\n2tw/8Oq6YnL6pzJlhHssa4/ABw/Dsj86gT3tNueM/eBmp6O9sqjlF03q65yJZ44K+rAe6XxgezKP\nDqPeoKkeVjwO783H1lXQkJROcoMbfBc8CP3aP7Foy+GaBr7hdj5/8/OjuGt2Szt6fZOP/353O48u\nKaRPcgIPXj6BuWcMbf5At9aycX8lz63cy8tr9lFZ10RO/1Qy+ySxrqiCvEwPd140lssmZnetbf7Q\nVqejdd2zzt/z9OudfjFjoGR9UB/aBjfQg/TLdS6qGHU+jJx1dHPPoa3w+j2w8z28GeP5VsWXWek7\nhf+5cSpnjWql47ux1gmEnR9AbTktJ4HuTxPnLMcnwGnXO7UJV1l1Pf/+xhaeX1VEdnoKXzlzOHPP\naKcloht6exgk4HQgXwDsw+lA/rK1dmNbj+lKGDT5/OyafzYD63YTFxdHSlKi2+4d/EeLc5pwhk1r\nuawrqPnE29DEn95cxcjlD3J5/CeU9RvHgBseJy67oM3X9fstS3eW8dLqfXx2oIqUxHhSk+JJTXRu\nKe4yTQ38dcU+RmT24S+3TOv5y/uCvTcf3n0IflzaqauNGn1+XlhVxMF//heX+t9jVMNnGPzOl5CM\nPM9pHx51AfQf1vyYsup6fv7aJl7+dD9jBvXhl1efxrjsvtz57FoWbSzhhunD+OkVBW0PzKoqcWpQ\n2xY7TRv1lc7fbcot8PkfdXwg0uHd8M+fw/q/Q1yi86E/aJzTJBNolkkf1vs+8DvCW+4EXvkOOOdO\n5/3cTQ1Nfh5cuIG/Ld/LnAlDeOS609lcXMV9L6xj28FqrjxjKD++bDyZ7ZxY1TX6WLzpAM+t2Etx\nRS1fP3ckV3ei/6ldh3fDx7+HNU85Z98BJt752w4ucPvQCqDfUNiz1Kn97Hgf6isA41zOOep854N9\n2R+dmtn5P4apX2PvkXq++uQKdpfV8BW3FaG20Uddg4/aRvfW4KOu0ceYwX25anIOM/Iz2w04v9/y\ntxV7mL9oKzX1TXz9cyP59vmj8ST13IWdvToMAIwxlwK/xbm09Alr7UPt7d/VmsE//3Qvg+MqGTek\nD3HHVuOwTrX+4CYoWedU48GprgeCISkN3n4QW3uE59Ju4IFDFzJhWBYPXVlAQU76Ua9VeLCal9YU\n8fKa/ew7Ukuf5AQmDe9PQ5OfusbgN5Bz39vQxLljBvKf109qtW21R6143GmHvWsr9G370shW+Rrh\n5wOdM+aJ1zgf/jlTTlhLenfLQR54aT3FlXUMTU+luKKWB74wnq/NzOv4lSGBq1U8WTCwiwOQasog\npd9Jd8kLjTetAAAO5klEQVRtJFhreeKjXTz0j01kp6eyv6KW7H4pPPTFiXz+1EGRLp6j+iCse87p\nXxpS4PS5tNcs5muC/ath+z+dW9EK5/Ng8o1ObSro+v4KbyPfeWYNH2w7hCcpwT2xi2s5sUuMJykh\njjV7jlBd79SArpw0lC9OymX0oD5HveyGfRU88PIG1u49woyRGfx8bkGHLhHurl4fBp3VE5eWHqWx\n1rkMrmiFc73v3hVQ7c4+OWQiXPlH7OAJvPLpfn7xj82U19Rz44wRfO2cfJZsPcSLq4tYW1RBnIHP\nnTKQqybnMnvc4KMuVzxWT4027ZCNL8Pfb4ZvfHxcZ9YJVZXAw2PhCw87TS6dUF3fxK8XbeGNDSX8\n6urTes8HirTr3S0H+cEL67ikYAj3zDm165ee9ka1R5yrooJqs8c60f9qbYOPtzYf4MXVRbz/2SH8\nFk7PTeeqybnMGjuQJz7cyVNLd5ORlsyPvjDuqGa1nqYw6C5rnc7f8h3O5XhBZ5EVtY385q3P+Osn\nu5oHBY3L7sfVk3O44vShDOrXtUvAwmrnB7DgMrhpodO80xklG+CPM+HaBTDhyp4pn8hJ6mBlHQvX\n7ueF1fvYXFwJOIMabzorj+/PPoX01PDWSDsaBlEU7yFmjHOm0MrZQnpqIvOumMA1U3J5f9shZp0y\niPFDO3fNecR1Z4h/d+aZEYlyg/qlcNu5I7nt3JFs2l/Jks8O8rkxA49rVu5tFAbdUJCT3uv/wG0K\nngmys7ozA6VIDBk/tN9Jc6IY+Xl1JTICg7m6VDNwH9OVScdEpFdSGMSq+ERnoEyXw8B0bUi+iPRK\nCoNY5sns2kyQNaVOzUJfQSgSNRQGsSwtq+sdyOo8FokqCoNY5snsYgdymfoLRKKMwiCWeTK6UTPQ\nlUQi0URhEMs8WU77f2cHHtaUqmYgEmUUBrHMk+lMr1xf1fHH+P3ODI0aYyASVRQGsawrA89qDztz\n+asDWSSqKAxiWVempNCAM5GopDCIZZ4u1Aya5yVSM5FINFEYxLLACOLODDxrnpdINQORaKIwiGVd\naibSjKUi0UhhEMuS+0J8UssHfEfUBPoM1EwkEk0UBrHMmM6PQvaWQnK/9r9WUEROOgqDWOfJajnb\n74iaUs1WKhKFFAaxrrNTUng1+lgkGikMYl1nZy71lqnzWCQKKQxinSez8x3IqhmIRB2FQazzZEJd\nBfgaT7yvtZqxVCRKKQxiXfNYg/IT71tfBb4G1QxEopDCINZ1ZuCZBpyJRC2FQaxrnrm0A/0GNZqk\nTiRaKQxiXVdqBhp9LBJ1FAaxrjMzlwb2UQeySNRRGMS65plLOxAGmrFUJGopDGJdfCIkp3e8mSgh\nBZLSer5cIhJW3QoDY8yvjTFbjDHrjDEvGWP6B2273xhTaIzZaoy5OGj9FGPMenfbfxpjTHfKICGQ\n1sGBZ4EBZ/qTiUSd7tYM3gIKrLWnAZ8B9wMYY8YD1wMTgDnAH4wx8e5jHgW+Doxxb3O6WQbpro7O\nXKoBZyJRq1thYK1dbK1tcu8uBXLd5bnAM9baemvtTqAQmG6MyQb6WWuXWmst8Ffgyu6UQUKgozOX\n1miSOpFoFco+g68Bb7jLOcDeoG1F7rocd/nY9RJJnaoZKAxEolHCiXYwxrwNDGll0wPW2lfcfR4A\nmoCnQ1k4Y8ztwO0Aw4cPD+VTS7A0Nwysbb8/oKZMYwxEotQJw8Bae2F7240xtwCXARe4TT8A+4Bh\nQbvluuv20dKUFLy+rdd+DHgMYOrUqbat/aSbPJngq4eGauerMFvTWAuNNQoDkSjV3auJ5gD3AldY\na71BmxYC1xtjko0x+TgdxcuttcVApTFmhnsV0U3AK90pg4RAR0YhNw84UzORSDQ6Yc3gBP4LSAbe\ncq8QXWqt/Vdr7UZjzHPAJpzmo29aa33uY+4AngRScfoY3jjuWSW8Ap3CNWUwIK/1fTTgTCSqdSsM\nrLWj29n2EPBQK+tXAgXdeV0JsQ7VDDRjqUg00whkaRk70N7AM81YKhLVFAbSyZqBOpBFopHCQCC5\nH8Qlth8GNaVg4iGlf9v7iMhJS2EgztgCT2ZLJ3FrvKXOPpqXSCQqKQzE4cls/3uQveXqPBaJYgoD\ncZxo5tKaUg04E4liCgNxnGh+Is1LJBLVFAbi8GS132egGUtFoprCQByeTKg7Ar6m47f5Gp1tqhmI\nRC2FgTgC/QG1h4/fFuhYVp+BSNRSGIijvVHImopCJOopDMTR3ijkwDrVDESilsJAHM0zl7ZSM9CM\npSJRT2Egjo7UDNRMJBK1FAbiaC8MAjWD1IzwlUdEwkphII6EJGfCulZrBqWQOgDiu/tdSCLSWykM\npIUno+2agfoLRKKawkBatDUK2Vum/gKRKKcwkBZtzU+kSepEop7CQFqkZbV9NZHCQCSqKQykRaDP\nwNqWdX6/molEYoDCQFp4sqCpDhpqWtbVHQHrUweySJRTGEiL1sYaaMCZSExQGEiL1sKgeSoK9RmI\nRDOFgbQInP0fVTPQjKUisUBhIC3arRkoDESimcJAWgTCIHjgmVfNRCKxQGEgLVLSIS7hmGaickjq\nA4kpkSuXiPQ4hYG0MMYdhRxUM9DoY5GYEJIwMMbcZYyxxpisoHX3G2MKjTFbjTEXB62fYoxZ7277\nT2OMCUUZJEQ8mS3feQxOMKjzWCTqdTsMjDHDgIuAPUHrxgPXAxOAOcAfjDHx7uZHga8DY9zbnO6W\nQULo2PmJNGOpSEwIRc3gN8C9QNAcBswFnrHW1ltrdwKFwHRjTDbQz1q71Fprgb8CV4agDBIqnsxj\nOpA1FYVILOhWGBhj5gL7rLVrj9mUA+wNul/krstxl49dL71FcM3AWvUZiMSIE351lTHmbWBIK5se\nAH6I00TUI4wxtwO3AwwfPrynXkaCpWVB7WHw+6DRC7561QxEYsAJw8Bae2Fr640xE4F8YK3bB5wL\nrDbGTAf2AcOCds911+1zl49d39ZrPwY8BjB16lTb1n4SQp5MwDqB0FDtrlMYiES7LjcTWWvXW2sH\nWWvzrLV5OE0+k621JcBC4HpjTLIxJh+no3i5tbYYqDTGzHCvIroJeKX7v4aETPDAs5qyo9eJSNTq\nkW84t9ZuNMY8B2wCmoBvWmt97uY7gCeBVOAN9ya9RfCUFIGagZqJRKJeyMLArR0E338IeKiV/VYC\nBaF6XQmx4DCorzp6nYhErR6pGchJrHnm0tKWMFDNQCTqKQzkaME1g7pKiE925iYSkaimMJCjJSRD\nUl+n87i+0qkVaMYQkainMJDjeTLcmkGF+gtEYoTCQI4XmLm0rlL9BSIxQlNYy/HSspyagVdTUYjE\nCoWBHC8wjXVNmUYfi8QINRPJ8TyZUFUC/kZIU81AJBaoZiDH82Q6QQCqGYjECIWBHC+401gdyCIx\nQWEgxwvuNFbNQCQmKAzkeMFhoJqBSExQGMjxgmsDurRUJCYoDOR4ngznp4mHlP6RLYuIhIXCQI6X\n0t8JAk8GxOktIhIL9J8ux4uLc4JAncciMUNhIK3zZKnzWCSGaASytO68e/Q9BiIxRGEgrSu4OtIl\nEJEwUjORiIgoDERERGEgIiIoDEREBIWBiIigMBARERQGIiKCwkBERABjrY10GTrEGHMI2N3Fh2cB\npSEsTiipbF2jsnWNytY1J3PZRlhrB57oSU6aMOgOY8xKa+3USJejNSpb16hsXaOydU0slE3NRCIi\nojAQEZHYCYPHIl2AdqhsXaOydY3K1jVRX7aY6DMQEZH2xUrNQERE2hG1YWCM+bUxZosxZp0x5iVj\nTP+gbfcbYwqNMVuNMRdHoGzXGmM2GmP8xpipQevzjDG1xphP3dsfw1229srnbovosTumLPOMMfuC\njtelES7PHPe4FBpj7otkWVpjjNlljFnvHquVES7LE8aYg8aYDUHrMowxbxljtrk/B/SiskX8vWaM\nGWaMedcYs8n9//yuuz40x81aG5U34CIgwV3+FfArd3k8sBZIBvKB7UB8mMs2DhgLLAGmBq3PAzb0\ngmPXVvkifuyOKec84O5IHy+3LPHu8RgJJLnHaXyky3VMGXcBWZEuh1uWzwGTg9/vwHzgPnf5vsD/\nbC8pW8Tfa0A2MNld7gt85v5PhuS4RW3NwFq72Frb5N5dCuS6y3OBZ6y19dbanUAhMD3MZdtsrd0a\nztfsjHbKF/Fj14tNBwqttTustQ3AMzjHS1phrX0fKD9m9Vxggbu8ALgyrIVytVG2iLPWFltrV7vL\nVcBmIIcQHbeoDYNjfA14w13OAfYGbSty1/UW+W419D1jzLmRLswxeuOx+7bbFPhEpJoVXL3x2BzL\nAm8bY1YZY26PdGFaMdhaW+wulwCDI1mYVvSW9xrGmDxgErCMEB23k/o7kI0xbwNDWtn0gLX2FXef\nB4Am4OneVrZWFAPDrbVlxpgpwMvGmAnW2speUr6wa6+cwKPAz3E+5H4OPIwT/NK6c6y1+4wxg4C3\njDFb3LPgXsdaa40xvelSx17zXjPG9AFeAL5nra00xjRv685xO6nDwFp7YXvbjTG3AJcBF1i3QQ3Y\nBwwL2i3XXRfWsrXxmHqg3l1eZYzZDpwChLyzryvlI0zHLlhHy2mM+RPwWk+W5QTCfmw6y1q7z/15\n0BjzEk7TVm8KgwPGmGxrbbExJhs4GOkCBVhrDwSWI/leM8Yk4gTB09baF93VITluUdtMZIyZA9wL\nXGGt9QZtWghcb4xJNsbkA2OA5ZEo47GMMQONMfHu8kicsu2IbKmO0quOnfvGD/gisKGtfcNgBTDG\nGJNvjEkCrsc5Xr2CMSbNGNM3sIxzgUUkj1drFgI3u8s3A72phhrx95pxqgB/BjZbax8J2hSa4xbJ\n3vEe7nkvxGnD/dS9/TFo2wM4V35sBS6JQNm+iNOmXA8cAN50118NbHTLuxq4PELHrtXy9YZjd0w5\nnwLWA+vcf4jsCJfnUpwrPLbjNLdFrCytlG0kzhVOa933WETLB/wNp1m00X2v3QpkAu8A24C3gYxe\nVLaIv9eAc3CaqdYFfa5dGqrjphHIIiISvc1EIiLScQoDERFRGIiIiMJARERQGIiICAoDERFBYSAi\nIigMREQE+P96CXCbPgiqJgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = -0.5 + np.random.rand(1000)\n", "f = np.sin(10 * 2 * np.pi * x)\n", "\n", "k = -20 + np.arange(40)\n", "f_k = ndft(x, f, len(k))\n", "\n", "plt.plot(k, f_k.real, label='real')\n", "plt.plot(k, f_k.imag, label='imag')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the NFFT shows strong features at wave numbers $\\pm 10$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fast transform: initial implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Kernel function" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# equations C.1 from https://www-user.tu-chemnitz.de/~potts/paper/nfft3.pdf\n", "\n", "def phi(x, n, m, sigma):\n", " b = (2 * sigma * m) / ((2 * sigma - 1) * np.pi)\n", " return np.exp(-(n * x) ** 2 / b) / np.sqrt(np.pi * b)\n", "\n", "def phi_hat(k, n, m, sigma):\n", " b = (2 * sigma * m) / ((2 * sigma - 1) * np.pi)\n", " return np.exp(-b * (np.pi * k / n) ** 2)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from numpy.fft import fft, fftshift, ifftshift\n", "\n", "N = 1000\n", "sigma = 1\n", "n = N * sigma\n", "m = 20\n", "\n", "# compute phi(x)\n", "x = np.linspace(-0.5, 0.5, N, endpoint=False)\n", "f = phi(x, n, m, sigma)\n", "\n", "# compute phi_hat(k)\n", "k = -(N // 2) + np.arange(N)\n", "f_hat = phi_hat(k, n, m, sigma)\n", "\n", "# compute the FFT of phi(x)\n", "f_fft = fftshift(fft(ifftshift(f)))\n", "\n", "# assure they match\n", "np.allclose(f_fft, f_hat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The expanded algorithm" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "\n", "\n", "def nfft1(x, f, N, sigma=2):\n", " \"\"\"Alg 3 from https://www-user.tu-chemnitz.de/~potts/paper/nfft3.pdf\"\"\"\n", " n = N * sigma # size of oversampled grid\n", " m = 20 # magic number: we'll set this more carefully later\n", " \n", " # 1. Express f(x) in terms of basis functions phi\n", " shift_to_range = lambda x: -0.5 + (x + 0.5) % 1\n", " x_grid = np.linspace(-0.5, 0.5, n, endpoint=False)\n", " g = np.dot(f, phi(shift_to_range(x[:, None] - x_grid), n, m, sigma))\n", " \n", " # 2. Compute the Fourier transform of g on the oversampled grid\n", " k = -(N // 2) + np.arange(N)\n", " g_k = np.dot(g, np.exp(2j * np.pi * k * x_grid[:, None]))\n", " \n", " # 3. Divide by the Fourier transform of the convolution kernel\n", " f_k = g_k / phi_hat(k, n, m, sigma)\n", " \n", " return f_k" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = -0.5 + np.random.rand(1000)\n", "f = np.sin(10 * 2 * np.pi * x)\n", "N = 100\n", "\n", "np.allclose(ndft(x, f, N),\n", " nfft1(x, f, N))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Speedup #1: using an FFT\n", "\n", "We'll replace this slow sum\n", "\n", "``` python\n", "g_k = np.dot(g, np.exp(2j * np.pi * k * x_grid[:, None]))\n", "```\n", "\n", "With the FFT-based version of the sum\n", "\n", "``` python\n", "g_k_n = fftshift(ifft(ifftshift(g)))\n", "g_k = n * g_k_n[(n - N) // 2: (n + N) // 2]\n", "```" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "from numpy.fft import fft, ifft, fftshift, ifftshift\n", "\n", "\n", "def nfft2(x, f, N, sigma=2):\n", " \"\"\"Alg 3 from https://www-user.tu-chemnitz.de/~potts/paper/nfft3.pdf\"\"\"\n", " n = N * sigma # size of oversampled grid\n", " m = 20 # magic number: we'll set this more carefully later\n", " \n", " # 1. Express f(x) in terms of basis functions phi\n", " shift_to_range = lambda x: -0.5 + (x + 0.5) % 1\n", " x_grid = np.linspace(-0.5, 0.5, n, endpoint=False)\n", " g = np.dot(f, phi(shift_to_range(x[:, None] - x_grid), n, m, sigma))\n", " \n", " # 2. Compute the Fourier transform of g on the oversampled grid\n", " k = -(N // 2) + np.arange(N)\n", " g_k_n = fftshift(ifft(ifftshift(g)))\n", " g_k = n * g_k_n[(n - N) // 2: (n + N) // 2]\n", " \n", " # 3. Divide by the Fourier transform of the convolution kernel\n", " f_k = g_k / phi_hat(k, n, m, sigma)\n", " \n", " return f_k" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = -0.5 + np.random.rand(1000)\n", "f = np.sin(10 * 2 * np.pi * x)\n", "N = 100\n", "\n", "np.allclose(ndft(x, f, N),\n", " nfft2(x, f, N))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Speedup #2: Truncating the basis function sum\n", "\n", "The expression of the inputs in terms of the basis functions $\\phi$ takes the form of a matrix-vector product:\n", "\n", "``` python\n", "np.dot(f, phi(shift_to_range(x[:, None] - x_grid), n, m, sigma))\n", "```\n", "\n", "The NFFT algorithm is designed so that this matrix will be sparse, as we can see by visualizing it:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8lNX1/99nZhIgIQkJgSwTdoKAYt0AcQVBwGpdvrb+\nrFtt/RbRutXautRW/dpWahVFqyAuba0oUq1KrRoWZVFAdlkFwp7JHiAJCSSZmfP743kmDDEkgQzJ\nZHLfrxevzPPc+zxzE2bOc+85536OqCoGg8FgiFwcrT0Ag8FgMJxcjKE3GAyGCMcYeoPBYIhwjKE3\nGAyGCMcYeoPBYIhwjKE3GAyGCKfFDb2IjBeRLSKSLSIPtfT7GwwGQ3tDWjKPXkScwFbgUiAHWAH8\nWFU3tdggDAaDoZ3R0jP6YUC2qu5Q1WpgJnBVC4/BYDAY2hWuFn4/N7A36DgHGF63k4hMACYAOHGe\nHUN8y4zOYDCcFMThgCgXvpgoauKUlLgyYh1V7D6UBKUuXBU+pKoG9fugGU6GcvYXq2q3E71+3KhY\nLdnna1LfVeuqslR1/Im+V0vS0oa+SajqdGA6QLwk6XAZ3cojMhhaBunQAWe3ZCqGpOMZ6eLyMSv4\nU8oSiv3VTNx+Hbvn9KbXa9vwHyhFa6pbe7hNR4Fq+9+BI6cz6usrJ/428/S93Sd+NZTs87E8q2eT\n+jrTtiU31kdExgNTACfwmqpOqtN+I/Ag1m9dDtyhqt80dK2IJAHvAr2BXcB1qrq/oXG0tKH3AD2C\njjPscwYDEhWNo0sC3gFuci+MIWPsbqb2e5cUZzTfm3Ev7oVeYtfn4SssQqurIQJ1mrSqCm+Ohw45\nHvp+CpuBaxhmt+aSQS5Nm28aTgQF/PhDci87JvkSQTFJEZldJya5E7hYVfeLyGVYE9zhjVz7EDBf\nVSfZCS0PYT0sjklLG/oVQKaI9MEy8NcDN7TwGOpHBHFF4eyayI6J/Th77Cae7/FfAO7PuYzlc07D\nvbCKjps9+Ir3ta3ZVBtBa6rxFRUhRUW4vwKdBBO5AIC+LAXA25oDNEQ8ilKjIXuU1sYkAUQkEJOs\nNfSquiSo/zKOLHIauvYqYKTd7x/AAsLJ0KuqV0TuArKwliNvqOrGlhzDMVFFa6rx5hfQ8/ECih6H\nGznfbiyjF9b/hzE0YYgIrvQ0Dp2aTs7IKEaNWcvT6Z9T6ffxi13X8G1WJu4vD+FcstE8oA2NEqoZ\nPU2MSQZxG/BpE65NUdU8+3U+kNLYQFrcR6+qnwCftPT7GiIYVbyeXKI8ufSZA7segesYYTcW0YMi\nq9vJHofDifOUvpQMTab40sM8NfzfXNe5lP9WduS+5dfTZW4nur2/Ed/BCvAbB0w4oii+prsEk0Vk\nZdDxdDu+eNyIyCgsQ3/B8VynqioijQ44LIOxBkMokahoyyV3ez++N/ZbXuj5H6IQ7veMY+mc03Av\nrKbjJg/e/ILm+f39Pnybt9Fl8za6vAmv04fX7aa+rAUw/vU2gL/pU4JiVT2ngfYmxSRF5HTgNeAy\nVS1pwrUFIpKmqnkikgYUNjZQY+jrQVwuHAnx+Pq5yT8/jrhx+bwycAYDoqL5c8mpvPHFSAY9l4s/\nvxB/VVVEBgUjiVqX3BMF7H8Cbq51yZXTy/j+DUEo4Avd2q/RmKSI9AT+DdysqlubeO1s4CfAJPvn\nR40NxBj6elCvF1/JPijZR+py4Dm4v9YVAJksM4bBYGgOIrj69OLA2SnkjlbuvzCLO7vsZFW1j4nr\nb0KzupL29mb85eWot2W/bccxo2+QY8UkRWSi3T4N+D3QFXhZRAC8qnpOI/HMScAsEbkN2A1c19hY\nWlQC4UQwefRtH2d8PNozneJhiZSOreCFc95lfEwVM8sTeWTZNQycVI7m5OOvqDS+a0OzmKfvrWrE\nndIg3/tetH76SaPp8QC4M/Ka9V4tiZnRG046vrIy2FBG0gZIegOeYxDP2W2ZrI4ov7WzSwLaK53C\n4V04PK6Ml854h5Gd/LxemsofvrqCQU/vRz35+CsrjcsvDFE0lK6bsMEY+jDGGR+P9k6naHgilePK\nefHMdxjdycffy7rzxJIrGTTJGI1ww3egFA6UkvwNMB2e4nSestsGsDKiHmoRiYIvAr9KxtCHMb6y\nMlhXRtd10PVVeJohPG23hY3REMHZrzelZ3av9bXenbibZYd93LH+RshKIu2dzfhKy4xbxhD2WDtj\nIw9j6OvicOIc2I/iYV3ZN/YQzwx9j6tjD/JhRWceWP4jkuZ0JHlFCb5vtxvDBaCKL3snnbN3MuBf\n8DGJfEwiAN35FjAphYa2hOBrjthOmGIqTNXF74OCYpI2HiR2eQx/zh7HxupDXB5Tys+GLKH4bD8H\nB3TBER0FEnkfCIOhPaNAjUqT/rUlzIy+HmpTK1cAz9dNrfwaiMzlncHQLBxOHJ06Iu5U9p/TjTt+\n/x63xBcz/5CTu1bfQOycznT7ej+yKzdsdwdbefRty4g3BWPoDQZDw4gg0dE4uyVTeVo6notdrL75\nOQ74vdy584dsz+pL+uJDRG314N93AP/W7cRv3c6MtzOYYWt09WQ90DYmSP42NltvCm3K0IvLBWcM\nJP+8eKLHFTF98FsMiY5i8v5Mpi4Yg/tziF+Tj3fXHpOFYjCECtVa+eToHA99PoNrHz7Xbswng3wg\nMmIxZkbfygRkCTznxxMzvoBXBs1gUFQUfy4ZxGsLR1pGfm0+/rxm6pUYDIZ2iyL4IjB02WYMfUCW\nIHXKEpgCv+bc2raA39zIErRzAi6G5K62iyGKMZeu4em0BZT7vYx+5Te4Fx1xMRjJYkN9GNdNO0Jc\nLhyJiXgz08m9MBb32D1M6z+TFGc0jxUOZ8Wj5xC7Pg9/UbERNgsXAi4GTy7Rnlz6ZMH2R+Bae1LQ\nw64pEAkuBsPJQRGq1dnawwg5xtAfA/V6j1Q7WgL8+Ui1I1A6sMKsIAwhJ1AztvI0q2bs+DEr+VPq\n4iOBz8/60uuNbLMiOUlYG6aM68ZgaBgRHB064EjpxsEhaXhGOdh4/Yvs9lZz+9YbKM5yk/ZVBa5t\nOVaB6xZWJgx36gY9t3BkRRIIfJoVycnFBGMNkYMIjk6dKL9sCLmXwF0Xz+XexGzWV9fw8403U5OV\nTOrSMlj77fEZY1X8hw/j372Xjrv30u9juPJXQwGIZjfp7AaM+8QQnqgKPjUzekOkoIq/spLY978m\n833IIp4szgIgia3A1rat4ReU+10xxMr9/uamKRT7q5m4/Tp2z+lN+uJKorZ4rJWFcYMYbPxmRh8Z\nSFQ0jqQu1Axwk3thJzLHbeevfd6ni8PFQ3kXs/HRIcRszMNXVIxWV5tAa1skKPe7Q46Hvp/CNQ8N\nsxtzySAXMCsLw9FYwdjIM4uR9xs1Aa2pxldQiKOgkIzFcOhPcFttoLWKaFaaQGuoEMGV4abitDQ8\nI11cOeZrnkxZRrG/mgnZ/4+cOb3o+fo2M6s2hAUmGGuIOCQqGt/wweReGEOvsbuY1m8WyY5oHik4\nj//OHYp7oZeOn69r3qpGFe/eHDrszaHvp7ABuIbAzNqDG4+ZVRvCCp/JozdEElpTjePLtWR8Cb6n\n4Oe1qxovfe2i2cZpFZmIy4WjSwK+/m7yzu9M0rhcpp8yg16uaP5UfAYLHj2fzuvy8BcUtat9IqHe\nGSsi44EpWHVfX1PVSXXaBwJ/A84Cfquqz9jnTwHeDeraF/i9qj4vIo8DPweK7LZHVPWThsZhDH0r\nEviyeTMzyLsglu7jcpie+TYZrg48UXg2Sx4dTuw6synLEHrU68VXXALFJaQtA56Fuzm/tr0jy9ut\n+9IfoqwbEXECLwGXAjnAChGZraqbgrrtA+4Brg6+VlW3AGcE3ccDfBDU5bnAQ6EpGEPfigS+bFJc\nQvpS4C9wZ+2sGrMpyxD5iCCuKJzJSRwe7MZzUQeGjd3AlIzP+MWey1kzdxDuhYfpsNmDr2T/SY/j\nWKJmIZvRDwOyVXUHgIjMBK4Cag29qhYChSJyeQP3GQ1sV9XdJzoQY+gNTaZWFmKAm9wLYmplIV4u\nuZAP552Le4GX2PV5+AqLTLaSoWmoojXVePPyceXl02s+FDwG13MesJ+etmxFS014FKGm6RIIySKy\nMuh4uqpODzp2A3uDjnOA4ScwrOuBd+qcu1tEbgFWAr9S1f0N3cAYekOTOUoW4iuCZCG01qdvViBh\nhsOJM74z2ttNwYgEvOMOMPX0Gaw/3IM/L/4+6fOdJKwuQHML8B861O4fzqocz4apYlU952SOR0Si\ngSuBh4NOTwWexFqAPAk8C/ysofu0D0MfvC3/9DRG/uErHk1ex07vYSZsuZF9WemkfXUQZ7bHbMs3\nRBZ+H74DpbC2lG5rganwf/bGuAGsAMxegqORUG6Y8gA9go4z7HPHw2XAalUtCJwIfi0irwIfN3aT\n9mHo62zLX/afKK7gbAA6sIs0dgENfOBtuQBJT6H0zO7kjlayfzCNVdU+bl93M2QlkbK0FMeuXPzl\n5eZBYTC0UZTjmtE3xgogU0T6YBn464EbjvMeP6aO20ZE0lQ1zz68BitruUHajqG3l6D5Px6Md9wB\npp/+Fud2dPLSgR48s3j80UvQysrQvrctF0D2Tjpn72TAv+D7E61ZUXe+tbpgZkYGQyQQqmCsqnpF\n5C4gCyu98g1V3SgiE+32aSKSiuVnjwf8InIfMFhVy0QkFitj5/Y6t35aRM7AMju76mn/Dm3H0NtL\n0G5Tl8JUeMyekYNZghoaICirY/vEvpxz6Sae7/FfAO7bezmr5gwmfVEVHTd78Oab6mTtHUVCWnjE\nzm//pM65aUGv88EurPvdayuArvWcv/l4x9F2DL2hcQKxiG7JVAxJI2eUi3U/nkKBr5qJ2dfjmdOT\n9MUVuLbl4t+/v324mIKyOno9lk/RY3Bjbb74gRbP6jCENwrUGK0bQ1gTiEXYkgP9PoFrfh2QG8jB\nTQ5gVj6GlsHZNQlfPzd5F8QRPzafVwe+Rf+oDvypeAhvfnERA5/PxZ9fGGabAcXo0RvCH3G5cMTF\n4e/nJv+8eP7+y+cYEh3F5P2ZTF04Bvd8iF+Tb21tN+l0hpOIr2QflOwjbTkwGe7jvNq2/iwLy1WU\nErqdseGEMfQRhnq9+Pbvh5X7SVkJD75wZH9GWyiiLi4XjoR4fP3c5J8fR+dx+UwfOIMPys7kjS9G\nkvGFn7i1+WE4EzRECmZGH4SI9ADeBFKwHoTTVXWKiCRhifH0xooIXxfYtSUiDwO3YXkP7lHVrGaN\n3hBxqNdbOxNMXQ48B/czAoBMlgHh/aBqCRxxcUiPNIqHdqXk0sM8O/xfXB17kNkVMdy/4joG/KES\n3ZuHv6IS/MZRdzyoipnR18GLtfV2tYjEAatEZC5wKzBfVSeJyEPAQ8CDIjIYK4/0VCAdmCciA1S1\nzX4SawuYDHTjubATA8du46XeHxDjcDLsH/fjXlhDp425+IpLjCSAIWT4y8thUzmJmyDxHzCV/ky1\n2/qx1sRgmoEVjG2yBEKb4YQNvZ2wn2e/LheRzVjaDlcBI+1u/wAWAA/a52eqahWwU0SysUR/lp7o\nGFqb4AImPRZCxR/gVluUrLeRBDCcDAKFXIak4bn4SCGXAl81d2w3hVyaj6kZe0xEpDdwJvA1kBK0\naysfy7UD1kNgWdBlOfa5+u43AZgA0JGYUAzR0NapZ3fy/RdmcWeXndYO5WfuJWVJO9idHFzI5RNT\nyCXUWMFY46P/DiLSGXgfuM/ezVXbpqoqIsftr7AV4KYDxEtS+Pk7RHD16UXpmSnkXgJ3XTyXexOz\nWVvtZeLGm6jJSib9rc2RbXBamnp2J39MIh+TCEB3lpjdyYaQEMrCI+FCswy9iERhGfkZqvpv+3RB\nQItBRNKAQvt8KAR+wgNVvDt2EbtjF5nvQxbxZNlCUUlsBbYag3McBPKtg7NsBkRF8+eSU2szbWI/\nW2eybAwnnVDvjA0XmpN1I8DrwGZVnRzUNBv4CTDJ/vlR0Pm3RWQyVjA2E1h+ou9viByOlWUDRzJt\n/K00ttZCoqJxdEmg5hQ3uRfG0GfcTl7u+y+SHC7OfOuXZCzwErMhF19RsQn0hxhTHPxozgduBtaL\nyFr73CNYBn6WiNwG7AauA7DFfGZhVVfxAr9oyxk3BsPJRGuq8RUV4SgqIuNLqAmq6dsk7X9bBNDf\nJ4OCEQn4x+7nlSAhwLcfv9zo0NeDKtT4jaGvRVW/hGPuLBh9jGv+CPzxRN/TYAgVR+0gHhFP1Lhi\nlp05s3YHcfrnkLC6De8gDujQryml+xrg5aOFADuzzLgX68Fy3RhDbzBEBHV3EPMifJ+zQIRTYjYg\n6SnsH55G3ujuPHrBf7gtIZ/Tlt1Ih6x4ui87gOzOxVd20GxIikDMzlhD20QEiY7GmdyVytPS8Vwc\nxZhL1/B02gJu2n4VW+f0w73oEFFbPfj3HWjf+deq+CsqYNsO4rbtIG4mzCKVWaTiZiPQ/uIF7QmT\nXtmKOGJjEXcq+8/pxh2/f49b4ouZf8jJXatvIDarM92W70d25eI7WGFmWPWhilZV4fXkEu3JpU8W\nbH8EruVcoIAeWJXJzF/OYDCum1bDX1EBW7cTv3U7M97OYIat09+T9VZ7aw6uOdibgByp3b+Tk7++\nuoZbJ/+S1KVlOLZ7TE6+wdBChLBmbNjQJgx9xOFw4jgtk6LhiVSOK+fFM99hdCcffy/rzhNLruRf\n88aSuLKIlG1LUVUz0zYYWggr68Zo3RhCgd+Hf923dF0HXV+FpxnC03bTAFYCbdyNEogJdO9GxZA0\nhv5hJU90//roSldfVuLa6mk/la4MbYJQb5gSkfHAFKyasa+p6qQ67QOBvwFnAb9V1WeC2nYB5Vjm\nwKuq59jnj6kQfCyMoTeEnkBMwNZkWfdJsB5LG6505XDiiI05SiL46WHvc23nMiMRHEGEynUjIk7g\nJawC3znAChGZraqbgrrtA+4Brj7GbUapanGdcw9Rj0JwQ2Mxht4QGoLjDWel4Bmt3HvRHO7usoOh\nq36MzulKytJSHDvasOiY3/cdieDp9LVEmTASwZFAiLNuhgHZqroDQERmYqn41hp6VS0ECkXk8uO4\n77EUgo+JMfQtQcCVkZTI4cFuPBdHs/Snz1KDcs+eH/DNnIG4Fx2mw2YPvpL9bTO90RYd89saQAPe\ng0/pwqecRTe2WF1og7N4Q7vjOLJukkVkZdDxdFuQMYAb2Bt0nAMMp+koVt0OH/BK0L2PpRB8TIyh\nbwlsV4avqJgOa6vpVe7mB8Nu4pWBM/hn77n86X/y+UfSRWTE9yLumw6mTJ7B0EqoCt6mG/rigN/8\nJHGBqnpEpDswV0S+VdVFwR2aqhAcmYY+aIPQoVPTyRkZxagxa3k6/XMq/T5Gv/Ib3F8eIurblt0g\nFFwmL3Z8/eJdbdChYTBEFCF03TRLsVdVPfbPQhH5AMsVtIhjKwQfk8g09EEbhKI8ufSZA7segets\nw9qDJUB4uxEkKhrfuafiuagTA8ZuZ2qf94lzuPhN3kjmzT0T98IaOixYb5QLDYYQEmIf/QogU0T6\nYBn464EbmnKhiMQCDrt6XywwFvg/u/lYCsHHJDINfQSgNdU4Fq+hx2I49McjJQrhMH1s9cI2Z96D\ns1aGdWXf2ENMHjoLPw4eWP4jkuZ0JHlFiclaMbQqoTL0quoVkbuALKz0yjdsFd+Jdvs0EUkFVgLx\ngF9E7gMGA8nAB3YhJxfwtqp+Zt+6XoXghogsQ+9w4uwcS9G1p3Lg0kM8P2wml8ccZtbBBB7++n9I\nntuRriuK8W3ZYYxIa1A3a+Xv8BIDAOjHGiC8V1mGyCfUefSq+gnwSZ1z04Je54O91f9oyoDvHeOe\nJRxDIfhYRJah9/vwlZWR9LelJP0NXmAgL9hN/Y0hadvUkYvwjFbuvmhebQnH/332vrafvmkIC4wE\nQnvDXiHQI42iYUlHrRL6zf9p7QpBc/KNq+FkUzd9s04JR1Mz1hAKVMFrCo+0M+wVAhvLSNrIUasE\ns0IwGE6AwMospRubH0jljovncX/iNtZX1zBh001UZ3UjdYkl5Ofb3+Cu/pOGkSk2GAxtC3tVqr3T\nKRqeSMXYg/z1rLfx1CTy2JKrSZ3nInFlEerJx3/o8MlflQZWZjt3k/mL3cwjjnn2qiyRbcC2Vl2Z\nmeLghtahjr7KvrGHeWbYv7g69iD9v/ipyVQxNExgVbqu7CgRPYgQAb2TgBpDb2hx6tFXmUp/pmIy\nVQwngO06Kf/+kKP0iFZV+5i4/iY0qyspy0rhmy3tNqBtgrGG9k2Qf7XszFQ8l8AdI+dxaeymo/2r\nO3Pxl5a1W0MR1tiuk9j3vj5Kjwg4SpOovaJqfPSG9kxQBlKhnYH0Ym0G0kSS53YkfUURmpOPz7iQ\nDG0WwWeybgztFpOBFJ6I4OjQAUe3ZCqGpJEzysW6H085usjL4gpc23JNkZcmYnz0BoOh5TiOjBmv\nJ48Oe3Po9wlc8+sIKPLSSoRY6yZsMIbeYGgi4nLh6JKAd0AGeRfEkjpuL9P6v0OGqwO/LxzK+/NG\ncMoLe/EXFYdGZtpkzLQ8Gpkage3C0IvLhSMhHl8/N/nnx/H2vc8yICqaP5ecyhtfjCTjCz9xa/ON\nDryhQdTrxVdcghSXkL4EeBrurBWbg34sNTLTEYDJummjBOvApy6H+58zOvBtmoBLo2c6RecmcnDs\nQb694J/MKO/K75ZeTfe5USStLEZz8lpmE5AhYlATjDUYwoSAS2NDGV03QNfXYBxnAJDJKsC4NJqL\nuFw44uLw906nYEQCjNvHK6f/k7Ojnbx8oA8zHx9PwppCNLcA/6FDEbUKjqBfpRZj6A3tBomKxpHU\nhZoBbnIv6ES/8Tt4uc97dHG4eCjvYjb+7nRiNuTiKypu9wVd1Ou1tGb276f7GuBl+B1Da9s783XE\nPkxN1k17owFp3GGrr8eflWykcdsQWlONr6AQR0EhGYuh6im4rdbHXkU0K4wLr52jagx9+yNYGnfX\nXgZ9Hs+c3iP454jxTH3gRc4928lLB3rwzKLLSJ8vEbuUNRjaEya9sj3j9x21lH3s5bNrmwawHDB+\nYYOhQUSQ6Gic3ZKpPC2dU59cx6S0hRzwe7lz5w/Z/llf0r88RNRWD/59B9Ca6lYZZijnaCIyHpiC\nVUrwNVWdVKd9IPA34Czgt6r6jH2+B/AmkIKV3j9dVafYbY8DPweK7Ns8YleyOibG0LckgQ96cley\n7+zFRWPW8Yx7HpV+H3fvvppNWQNwLz5E1Let+0E3GE4KqmhVFd4cD9E5HrZ9Btdyrt2YTwb5QOtO\nmBTBH6KsGxFxAi8BlwI5wAoRma2qm4K67QPuAa6uc7kX+JWqrhaROGCViMwNuva5wEOhKRhD35IE\nPuieXHr/Npc9v4XrCKR6FtODYsCsDEKKw4kzvjPa243ncWXq6TM4v6ODaQfc/Hnx90mf7yRhdYHl\ncqusbO3RGsKAEE7ohwHZqroDQERmAlcBtYZeVQuBQhG5/KgxqOYBefbrchHZDLiDrz0emm3o7afW\nSsCjqleISBLwLtAb2AVcp6r77b4PA7dh2bJ7VDWrue//HRxOnIP6UzQsidKxlTw39N2jC4TP6UjS\nhxuMdnt7we/Dd6AU1paSdjX8n63UCDCAFYB5sBqCCG0w1g3sDTrOAYYf701EpDdwJvB10Om7ReQW\nLNv7q4CNPRahmNHfC2wG4u3jh4D5qjpJRB6yjx8UkcHA9cCpQDowT0QGqGpov2d+H76NW74jvAVH\nxLf8IX1Dw4lwVJ72476jcrQnLx5ngtuG1qPpH7VkEVkZdDxdVaeHcigi0hl4H7hPVcvs01OBJ7FG\n+iTwLPCzhu7TLEMvIhnA5cAfgfvt01cBI+3X/wAWAA/a52eqahWwU0SysZY2S5szBkPoCGi5+Pq7\nyTu/M13G5jH/tPf4Q/HpvPX5hWR87qfzujz8BUXNloo4Kk/7qqNztE1w29CaHMeMvlhVz2mg3QP0\nCDrOsM81CRGJwjLyM1T130fGpwVBfV4FPm7sXs2d0T8P/AaICzqXYvuXAPKxosZgLWOWBfXLsc99\nBxGZAEwA6EjMkSBm925svbsHV49ZxhPdvz5aivXLSlxbPUaKtRkEtFwoLiFtGfAsXIGVXdTfSEUY\n2gEK+P0hc92sADJFpA+Wgb8euKEpF4qIAK8Dm1V1cp22tCAbew2wobH7nbChF5ErgEJVXSUiI+vr\no6oqIsc97bOXP9MB4iVJa4OYe3Po+5sc1gHXYKRYwx673m3J/5x2VK3bDys688DyH5E0tyPJy0vw\nfbvdxEsM4YECIfLRq6pXRO4CsrDSK99Q1Y0iMtFunyYiqVh+9njALyL3AYOB04GbgfUista+ZSCN\n8mkROcMe7S7g9sbG0pwZ/fnAlSLyfaAjEC8ibwEFgSeOiKQBhXb/Zi1jDG0Qu95t4j+WHlXrFky9\n2zaFCI6YGMSdyv6zu5E/xsvOy15jwSEHd665gU5z4uj+9QFkdy6+soNt/qEdynCQbZg/qXNuWtDr\nfCxbWJcvoX4ZTVW9+XjHccKGXlUfBh4GsGf0D6jqTSLyF+AnwCT750f2JbOBt0VkMlYwNhNsZ2y4\nYEselF1xOrmX+Hngok/5RZe9LDvs4/Z1N+GYk2hJHuzMsTI5DIb2gCr+igrYup34rduJf+eIiFwP\n22sQUQkOERj3Pxl59JOAWSJyG7AbuA7AXrLMwsoD9QK/CHnGTXOxJQ86z1rGgFkwm67MpisAqWy2\nunCSZqFGetdgCAPEaN0cC1VdgJVdg6qWAKOP0e+PWBk6bRd7GVv6gyHkjfbx8AX/ZUJCLosOw51r\nb6RDVjzdlx3Av37r8RljI71rMIQHZkZvCCxj42YuI24mvE933qc7AG42AhG2jI1kRBBXFM7kJA4P\ncuP+4zYmZ3wKwH17L2fVnMGkL6qi46YcfCX7jSRFe0BBQ5d1EzYYQ29ov6iiNdV48/Jx5eVTMAJu\n5Hy78QA9WQJEZkqpIzbWCq6e0438S2t4YsRH3BJfzPxDTu5afQO9n6xBduXiO1jRDt2ExtAb2giB\nOrmemwZxlQ/AAAAgAElEQVTSaVwh0we/xalR0Tyz7xReWXgJ7vkQvzYf7649ZtdpO+So4OrbMIMM\nZtjJHz1Z375XpRH4dTCGPkIJ1MlNnbIEpsCva1UCIdOWzIjEmWpYIIKjQwccqd0pPyOVnEuEn41c\nyINdNzJy/Y84+FkqqUvKcW734C8tMxv8wg1j6A3tmVp9mr7pFJybgIwrYdqQt/i6sj+TF4/DPU9I\nWF2AP7+wfevTqOI/fBj/rj102rWHzA9hMR1ZzNl0Zged2QGY4HpYEsINU+GEMfSGJlOrT7NqP91X\nAS8d0agJ6NOYuWkDOJw4+/eu1y/+i1U30HlOZ7rP2thO/eLhQyTOT4yhNxhCiKNjR8tl871UPKMc\n3DJqEY8krye7porbt9xIp8diSczaSsKsMmZ4j/jFe7EeMLP8sMBk3RgM4YtERePsmkjVIDeeizpy\n5tjNvNDzY5wI9+aMx/NoJh03efDt249WV5+UqVuwy6b/R7CE6FphuE7sBIwxD3eOX50r/DGG3hAx\naE013vwCnPkF9PwCSp4ITpcsw8Uq41o6FrYAnWSkUnJOMr/63dtc17mUzyo7cM/K/0fCnFi6fb0P\n9uZFtmtJMcFYw0miTtFkz0gX48es5Dfdv2Dijh+xM6sP6YsridriwX+g1GzcMYQeW4COzeV02byN\n1//Zh9ftpj6sA9rLSkRMMNZwkqhTNLnPZ7AFuI0LgDwyrNKR7eSL1g4J6Bz1TmfX76L461lvM7qT\njzfLknls6VWkzo0icWUR6sk3OkctgZnRGyKaQP53SjcOnp5GzigHm/7fi+z2VjNhy42UzEkn/cuD\nOLPtlYXJ/w4NAZ2jdWX0/BE8zRCetpsGGJ2jlicCd4sZQ284QiD/e/deOu7eS///wJX3W+mTHdhF\nOrsAY3QaJeCKS0rk8GA3OSOjOe/SDUx2Z3HHniv4Jmsg7sWH6bDZYzR0wg2TR28wtDwBKYecnwwk\nfmw+rw58i/5RHfhT8RDe/OIi3F/4ifsmH+/uveGTAB1wxdkaOr3nQ+7v4HrOA/ZFtIZOJGCybgxt\nFnG5cCQm4s1MJ/fCWNxj9zCt/0xeLrmQD+edi3uBl9j1efiLiptd+DuUBKQc0iYvgclwH+fVtpk6\ntoaTQnh89EOKMfTtBPV68RUVIUVFuJcAf4aJXAAofVkKGINpoNbtVHPBaXgujmbEWMvlVINyz54f\n1LqdXF9uMC6nJiAi44EpWDVjX1PVSXXaBwJ/A84CfquqzzR2rYgkAe8CvbFqxl6nqvsbGocx9AZD\nEwm4kXz93ORdEFevK2ng87mW1k8YrYqOC9vt5Jq/il7zIff3AZcTBLud2uBv1mRC5boRESfwEnAp\nkAOsEJHZqropqNs+4B7g6uO49iFgvqpOEpGH7OMHGxqLMfSGiEdcLhxdEth76ykkj/PwyoC36eWK\n5smis3jn8/PJ+NxP5/V5ePfkNGicA24kSvaRtpx6XUlmVdTGUUIpgTAMyFbVHQAiMhO4CqucqvV2\nqoVAoYhcfhzXXgWMtPv9A6u6nzH0hrZPrXJmPzf5I+KJGlfMq6f+kyHRUUzen8lHj40hYXU+/oKi\n7yhnqteLr7iE9GeWwDNwd+1uWePnN9RD02f0ySKyMuh4uqpODzp2A3uDjnOA4U28d0PXpqhqnv06\nH0hp7GbG0BtaBYmKxpHUhZqBbuKezOGl3h8Q43Dym9xLWDD3DNwLa+i0MRdfcQlaXX1EOXPlflJW\nAi/Cg0HfmVi+NsbaEBKOw3VTrKrnnMShNIqqqkjjIzaG3tAqaE01voJCHAWFVFwEt3KB3XKI3uEY\nHBbB1acXpWel4LlEuffiOdzdZQdrq71M2HAT/qxkUpaWwjdbzEaytk7oAhAeoEfQcYZ9rrnXFohI\nmqrmiUgaUNjYzYyhN7Q4AZ+5d0AGeefH8uFdT5Ph6sDvC4fy/vxzyfjCR+y6MEv1VMW7YxexO3Yx\n4D34lC58ylkAJLMV2BrRAcp2Rej+I1cAmSLSB8tIXw/cEIJrZwM/ASbZPz9q7GbG0BtanIDPXIpL\nSF8Cd/7lgtq2fkE+c3G5cHbpgr+fm31PVDF98Fu1PvmpC8dYdW/X2H75yspW+m0MkYRo6LJuVNUr\nIncBWVgpkm+o6kYRmWi3TxORVGAlEA/4ReQ+YLCqltV3rX3rScAsEbkN2A1c19hYjKE3hC3BfvnE\ny4/2yZ9I3VuJisbRJYGaU9zkXhhDn3E7ebnvv0h2RPNQ/vmsf/QMYjbk4isqPml69YY2QAgLj6jq\nJ8Andc5NC3qdD3b1mSZca58vAUYfzziMoW8qAYXBnukUD0vkyUfeYHxMFTPLE3lk2TV0nxtN0opi\nNLcAf0WlURgMQ7SmGl9REY6iIjK+hJqn4Oe1sYEaolkRXnEBQ6tgJBDaMwGFwQ1lJG2A594YxHN2\nUyargbYv9nVU8e8RCTjGFbP8rJm8eKAvUxaOxT1fSFhjin8bIpwI/FgbQ2+o5TvFv/8K37cDjiEp\n/u1w4kyIx9/bepD4x+7nldPfYtXh3jyz6DLS5wsJawqtVZF5kBhagxD66MMJY+gNLYffZz1I9u+n\n+xrgZXjMrqcaeJC09VXRUQTr+w9JwzPKwcbrXyTHW8WEbTdQmJVB2pcVuLblGH3/cMIY+jAkSPv7\n0KluPBdHs+zWZ6lBuWv3layfcwruRYfp8K3R/ja0MHX0/ft9DFf+ytL3d7GHdPYAEfZwiwDEFB4J\nQ4K0v6Py8uk9D67/XUB/pKTFtb8DwleeWwYSNy6fVwbOYEBUNH8uOZU3vhhJxudhqJ9uMBgimrZv\n6MOMgPBV6nNL4Dm4nxG1bZlGV6VtUyfGwLh9TB0yg3M7Onlxfy8mLx7HoGeKTIyhrROB/23G0Bsi\nExEcMTFIegoHzupO3mgfO694la8O+7lj3Y24srqQsrQU2eXBV3awaemwDcQYwIozGDdMG8cEYw2G\nNoQq/ooK2LaDuG07iHsXxnEGAGlsBiKyBnTj2AHiinGn4xnl4JZRi3gkeT3ZNVX8/NubKM9KJfWr\ncmT15vYbHDaG3mA4eQQ0cHz93eSd35mkcblMP2UGvVzRDH73bjK+8NN5XZ4leRAuGjhtDTtA3Omj\n5fT/CJYQzRX2qiSWHcSyw+rWmmNsbSLwl293hl6iovGdeyq5F3ai37gdvNznPbo4XDyUdzFz5p2F\ne6GXDl+sM1vgW4GABg7FJaQtA549oh1vdOMNLYFgsm6+g4h0AV4DTsN6Dv4M2MIx6hmKyMPAbVgZ\nZfeoalZz3v9E0JpqHIvXkLEYqv4Et9Vuga+ijy2P21rmPVij3XNhJ+ZPeLpBjXbzIDIYQozx0dfL\nFOAzVf2hiEQDMcAj1FPPUEQGY0ltngqkA/NEZICqmviVTbBGe4+FcOsfwlyjvaUQwdmvN6VnpZA7\n2seDF37CxC4elh32MWHdTbiyupD6zqamB1UNhoYwhv4IIpIAXATcCqCq1UC1iByrnuFVwExVrQJ2\nikg2Vl3EpSc6BgjSNs+0tM2Tx3n4bPB7PFF4NrM+Pw/3F1Y9UOPXbcOo4sveSefsnQyYBR/QjQ/o\nBhwJrBrzbggZEWgimjOj7wMUAX8Tke8Bq4B7OXY9QzfYjlaLHPvcdxCRCcAEgI7ENDiIo7TNlwLP\nwJVYuw/7Gb9uRCAuF/6hp9Y+yOsr7h0zb515kBtCgnHdfPfas4C7VfVrEZmC5aappan1DOtiF9id\nDhAvSRH4Zz8ORHB06oSkp1iui0v8PHDRp/yiy16+t/zHOOYkkrK0FMfOnIh1XajXiyz9pvZBXl9x\n7wiMnxlaiwi0OM0x9DlAjqp+bR+/h2Xoj1XPsDn1E9svqlb1pCDXxWy6MpuupNpuC8W4LgxhRqDG\n7pkp5F4Cd108l3sTs1lfXcPPN95MTVYy6W9txl9eHl75+mqybo5CVfNFZK+InKKqW7Aqnmyy/9VX\nz3A28LaITMYKxmaCLVkY4dTVeZexJaw4+x1L532RrfO+2ui8GyKIoBq7me9DFvFk2ZLXSXaN3bCd\nnETg16+5WTd3AzPsjJsdwE8BB/XUM7RrJc7CehB4gV+0l4ybk67zbjC0ELWifTcNJGZ8Aa8MmsGg\nqCj+UjKY1xaOxP05xK/Nx59XgP/w4dYe7gkRSh+9iIzHyk50Aq+p6qQ67WK3fx+oBG5V1dUicgpW\nmnqAvsDvVfV5EXkc+DlWjBTgEbvs4DFplqFX1bXAOfU01VvPUFX/CPyxOe8ZChyxsUhGGvuGJlM4\npoYnR3zIjXElzKmM4u5VP6bvk9XInlx8Bysi0udtMJwotaJ9U5bAFPg159a2nUgd37AkRIZeRJzA\nS8ClWK7uFSIyW1U3BXW7DMu7kQkMB6YCw20vyRlB9/EAHwRd95yqPtPUsbS7nbGApYGyJZuELdkk\nvAVv0oM37fBBb9aZwF591CmikXOJgx9f8hW/67aa8Zt+SHGWm7SvTBENQxtHCaXrZhiQrao7AERk\nJlaaebChvwp4U1UVWCYiXQIxzqA+o4Htqrr7RAfSLg294QSoU0Sj/8ewAidXMpRodpOO9RmM1PWP\nREXjHzYYz0Ux9Bq7i2n9ZpHsiOaRgvP477yhuBd4iV2fi9eTa2IsbRjhuFw3ySKyMuh4up0xGMAN\n7A06zsGatdNIHzcQbOivB96pc93dInILsBL4VUB94FhEnqEXwZWawuFBbnIv6sDZYzfxfI//AnDf\n3stZNWcwfV/Zga94H+qtMV9KQ5PQmmrkq7VkfAW+p+DntdIZXvq2513LEchxGPpiVa3PdR0y7Pjn\nlcDDQaenAk9irT2eBJ7Fkp85JpFn6FXx5uXjysun5+dQ9DjcWJt3fYCeLDFfSEPosPc5OFK7U3pW\nCj98Iou7u+xgbbWXn6+/GZ3T1drnsCM3/FIJDfUTurlfU1LKG+tzGbBaVQtqhxf0WkReBT5ubCCR\nZ+jbGfWlbk4b8hZnRzt58UBf3ntsnEndPJnY+xz8dirhp+914VM7o6obW6wuRK5LKyIJ3VdkBZAp\nIn2wjPf1wA11+swG7rL998OB0jr++R9Tx21Tx4d/DbChsYEYQ9/GqS9183e2BARALF+bFYwhdIgg\n0dHUXHAanoujGTF2A5PdWfhQ7tlzBWvmDMK96DAdNnvw5hc0fr9wI4TqlarqFZG7gCys9Mo37DTz\niXb7NOATrNTKbKz0yp8GrheRWKyMndvr3PppETnDGi276mn/DqJhPsOLlyQdLvVma54cRHC506kY\nko5npIvLx6zgTylLKPZXM3H7dezJ6k364kocyzehNdUtNy6DwdAo8/S9Vc3xm8d076EDrru/SX2/\neen+Zr1XS9J+Z/QOJ87OsWjPdIqHJVI6toIXznmX8TFV9J13Jt3musj8ezFb/tCJaypG2Pn0ubjJ\nBSJy85zBYMBIIEQWfh++sjLYUEbSBkh6A55jEM8BmawGjF81nHDGx9f7UJ5Znsgjy65h4KRyNCcf\nf0Wl2eRmaBZGvdJgaCWO9VAG68FsTHsYErzJ7vQ0ckYd2WS321vNtS/+Ovw22YV2w1TY0OYMvTM+\nHnqkUTQsidKxlTw39F0ujznMrIMJPPz1/5A8pyNJH24wMzuDobWpu8nuP0c22QGkswQIw5WzMfSt\nj6+sDDaWkbQRkv4GLzCQF+y2/qwBWkGb3OHEGd8Z7e2m4NwEasaVsn7420wvTeepxZeTNt9Jl9WF\naG6BJTkc5gFwg6G9cpw7Y9sMbc7QhyV+H74DpbC2lG5rgWkwztIjYgArgDCctRgMDeCIi0N6plM8\nNIl9Yw8xeegsroyt5P2D8fxm+bWc8qcKdG9eRK6cxR95lt4Y+pYkaBflpl93596L5nB3lx2sqvYx\ncf1NaFZXUpaZXZSG1sdfXg4bt5C4ERL/Di8xgJfstv6sidyJi/HRG5pN0C7KAXfs4lNO8i7K4D0B\nF7u4bMxKJqV+RbG/mjt3/IidWX3o9do2KxBm9gQYDIBx3RhambpyB6/96nnOiHY1WKnKm+OhQ46H\nvp/CFuAahtl3yyODvMidmRkMJ4ox9IbWpK7cwSN/HVbbZipVhT/icsH3TiH/vAScY4uZftpbnBHt\nYsr+/ry4aAzu+ULcp+uNJlErY2b0hsjE4cTRqSPiTmX/Od3Iv7SGJ0Z8xC3xxQz66mZiszrTbfl+\nZJeputUc1OuFVRtJWQW8CI/w3Qd1BG7KbHsYQ2+ISPw+q+rW1u3Eb91O/NswgwxmkEFP1ltdWnmI\nJ4pERePokkDNKW6iniw8umDI3KG4F1oFQ3xFxWhVVWsP19DaqJFAMBjaHFpTja+oCEdREb5RpmDI\ncVGnfKRnlIPrLlnCY91X8f3N11KYlUHal2G2s7WZmDz6VkKioqgeMxTPSBfjx6zkN92/4M6dP2T7\nZ31JX1xJ1LZc/PsOmKwRgyHU1NnZ2u9jWIWDKxmKiz2kswcIkz0itnwyh0NwrwiMj4S9odeaGqI/\nW0Gfz6yskdu4AMgng3wgTD5kBkOABmbBYavvEgmohsz1Zmb0hvaFCK7ePSk7IxXPJXD7yM95IGkL\nG2uqmbDpJg5ldSdtSTms2WwMVoAGZsEQxvouBguzYcrQ7lDFu3M3MTt3k/kBfE4sn9sbvBLIJoHs\nSPxOGE4WgdVOWgplZ6Ry2RML+HXXTWyuqeH2zTdS+VkKaUvKcezw4C8ta7XJgwnGtnccThyxMUiP\nNLY+GsMzQ9/j6tiDfFjRmQeW/4ikOR1JXlFiaYCUl7f2aA2G8CKw2rEnDws/6MRCzgYgnu3Esz0s\n6uuG0tCLyHhgClYpwddUdVKddrHbv49VSvBWVV1tt+0CyrH+JN5ANSsRSQLeBXpjlRK8TlX3NzQO\nY+iPB7/PMuCbyul3A0ylP1Ptpn62cmZrf0gNrUMgjXP3/2bSa+yuY6Zxej25ERnsixiUkP3/iIgT\neAmr7msOsEJEZqvqpqBulwGZ9r/hwFT7Z4BRqlpc59YPAfNVdZKIPGQfP9jQWIyhNxhCQCCNM+Op\nInxPmTTOZhOID51pxYfuGDmP+xO3sb66hgmbbqI6qxupS8pg7bchd/GEMBg7DMhW1R0AIjITuAoI\nNvRXAW+qVbx7mYh0EZE0Vc1r4L5XASPt1/8AFhBphl6ionEkdaFmgJvcCzrRb/wOXu7zHl0cLh7K\nu5g5888i8697rA0w1dVm9mQwtEWC40P/hnnEMc+ODyWyDdh28uJDobuxG9gbdJzD0bP1Y/VxA3n2\nSOaJiA94RVWn231Sgh4E+UBKYwNpc4Zea6rxFRTiKCgkYzFUPRVIuQSoog9LzczJENE4OnbEkdqd\n8u+l4hnl4JZRi3gkeT3ZNVXcvuVGOj0Wh3NbTqsGNNsqx7lhKllEVgYdTw8yxqHgAlX1iEh3YK6I\nfKuqi4I7qKqKND7iNmfoDYaThXTogDO5K4cGpeEZGc1Fl67jGfc8qtTPXbuvZN+jvejwrQdfyf5W\n3aDnP3wY/649dNq1h/4fwRKiucIOanZiJ2BiRSeM6vEUHikOBEiPgQfoEXScYZ9rUh9VDfwsFJEP\nsFxBi4CCgHtHRNKAwsYGagy9wWCjVVV4PblEeXLpPQ/2PArXMcJuLcFJSUSvFoN1gTr8oYBpff9F\nF4eLR/Iv5LN55+Be4CVmQzvQBQqd62YFkCkifbCM9/XADXX6zAbusv33w4FS24DHAg5VLbdfjwX+\nL+ianwCT7J8fNTYQY+gNBgNwtC5Qzchgl2gNfU5GQNmWLXAmJXJ4sJuckdGcd+kGJruzqEG5+NVf\n4158mA6bW3YVFapgrKp6ReQuIAsrvfINVd0oIhPt9mnAJ1ipldlY6ZU/tS9PAT6wsi9xAW+r6md2\n2yRglojcBuwGrmv8dwrzYGW8JOlwGd3aw2i3BNwZlael4xkZxZgxa3g6bQEH/F7u2nkt27L60ev1\nbKM3ZAgL5ul7qxpxpzRIXEKGnnXBPU3qu+iTB5v1Xi2JmdEbGiTgzoj25NInC7Y/DNdyrt1aQAYF\nxh/cGIEdoand2fzLdH42agEPdt3I1ppqJnx7IwezUkn9qhzndg++kn2tPVpDeM99T4hmGXoR+SXw\nv1h/mvVYy44YjrFrS0QeBm7DihXdo6pZzXl/g6FNENgRumsPmffuYTEdWWwHTzuzg87sANpQAFUE\nV2oKhwe78VwczYixR9wt9+z5Ad/MGYh70WFcX25ok6s8I2oWhIi4gXuAwap6SERmYQUbBlPPri0R\nGWy3nwqkY+WHDlDVBj/fmadXMPk/S2tnPu43vzVpYwZDa6KKNy8fV14+veZD7u/hes6zG/fR0xZu\nO+n2UgRXrx6Un5FKzigHPxu1gHlDQnDbpmfdtBma67pxAZ1EpAZrJp8LPEz9u7auAmaqahWwU0Sy\nsdKFljb0BtvWxXJ/7xG1M582M+sxGOojyI1TXx789VMeqHXjmAlNI6jitdNMMz+ExXQMwT0xrptg\n7ET+Z4A9wCFgjqrOEZFj7dpyA8uCbhHYAfYdRGQCMAGgIzEnOkRDKBDBERNj1ZM9uxv5Y7w8dv5s\nbo0vZMEhB3euuYFe/+dF9uThKzto6sk2RpAbp748+FQjY9yqWBumIs/SN8d1k4g1S+8DHAD+JSI3\nBfdp6q6tuti7y6aDlXVzomM0hADVo+vJvgPvkM47pAPQgw2hrScrgqNTJyQ9hdIzu5M7Wrn/wiyG\nx2Rz+7qbYU4SKUtKcezKxV9ebma8htBjZIqPYgywU1WLAETk38B5HHvXVlN2iRnq4nDijO+Mv08G\nBSMS8I/dzyunv8W5HZ30mT2B9M8dJKwuQHML8B861Pa1fVTxV1ZC9k46Z+9kwL/gYxL5mKF051ur\nC2bGazh5mBn90ewBzhWRGCzXzWhgJVBB/bu2ZgNvi8hkrGBsJrC8Ge/fPvD78B0ohTWldF8DvAyP\n2cv8AfafrzWNnkRF4zvvVDwXdmLwuK282OtDYhxOfpN7CQvmnoF7YQ3RC9cbgTlD28D46I9GVb8W\nkfeA1Vgb5tZguVs6U8+uLXtH2CwsiU4v8IvGMm4M4Y/WVONYuIYeC6H8D3Br7W7KQ/S24+wR+L1p\n8zjj46FHGkXDkyi9tJLnhr7L5TGHmXUwgYeWXcvAp8rRnHz8FZXtLO5yXFo3bQazM9ZgaGkcTpyd\nY9He6RQNS+Tg2IO8cNZMxsbU0CfrNlLmRZG0shjNycN/6HA7M7TNo7k7Y+Pj3DrszDub1Hf+4kfN\nzliDwXAM/D58ZWWwroyu66Dra/Asp/IsMIBVQPPcceJy4UiIx9/XTd75ccSMK+C1QW8xICqav5QM\n5rUFIxk0OQ9/fiH+qirjUgtGTc1YQ3smMAvtmU7R8ETKLq3ghXPeZXxMFX3n/ozuc6NJWlFsBYXb\n3XI/vFCv15JSKNlH6grgebi/VoUTMvk6olU4m00EPviMoW8r2Nk32ttNwbkJVI8tY8O5M5hems5T\nX15O2nwnXVYVWoa2sjL0H9bALHRDGV03QNfX4TkG8RyQyWrAZMIYIoTIs/PG0LcZAtk3a0vpthaY\nBuM4A4ABrACMoTWEEIcTR2wMkpFKyTnJFF1axaRz3+e6zqX8t7Ijj//xp3Rbvg/25uE7WBFRKzjx\nR57vxhh6g6G9E7RJ7cDZKeSN9rHzilf56rCfO9bdiCsLTpl8mL/vOoPX7d3PSSyNzImFYjZMGQyG\nMMYu5FE1cgiei6MYc+kankr7gkq/j1/suoZvszJxf3mIqG89+AqCqs8FbVKLy95J3LtHVotpbAYi\n0vbVi6Bmw5TBYAhjVNGqKqKzVlq1Ax4JLoVYRA+KgNC6+I7K8Dkvjk7jCpk++C1OjYom89934P4c\n4tfm4y8oajs7t9vCGI8TY+jbA3XytivGHeSvZ73N6E4++mTdRurcKBJXFqGefJO3bTguvpPhMwV+\nbRemyeRrIMTlB1sCY+gNbZJ68rafZghPE5q8bUOY0Ugt1nv2/IDiR3u3eC3WNkGIffQiMh6YglUz\n9jVVnVSnXez272PVjL1VVVeLSA/gTSz1XwWmq+oU+5rHgZ+DvUSDR1T1k4bGYQy9wRBp2C6cQHGQ\n3vMh93dHFwdxsg+vw4mjU0ccfXuy75xkCi+t4ckRH3JjXAkDv7yZznM6023ZfmRPbsRl1jREqLJu\nRMQJvARciiXLvkJEZqvqpqBul2HpfmUCw4Gp9k8v8Cvb6McBq0RkbtC1z6nqM00dizH04UCgGEVa\nCmVnpOK5BP734gVcEf8Nt2++kcrPUkhbUo5jhylGYQghfp8lQb0lm4Qt2STMgDfpwZv0oBfrrS6t\nPMSWR0PpuhkGZKvqDgARmYkl7R5s6K8C3lRLi2aZiHQJqP8CeQCqWi4im7Hqd2ziBDCGPhwIFKPY\nuZuYnbvJ/AAW0omFnEs824lnu5HmjSTsmbRkpLFvaDKFY47MpE9ZfAtxc2Lp9nX7m0mHBcrxGPpk\nEVkZdDzdrqURwA3sDTrOwZqt00gfN7aRBxCR3sCZYAc9LO4WkVuwFIN/FajLfSyMoT8WgcpK6Skc\nOLs7eWO8PHr+x9yWkM+CQw5+9fTtdF92ANmdayorGY6PujPpt47MpHuzzurSykNs1zT9j198skXN\nRKQz8D5wn6qW2aenAk9iPZaeBJ4FftbQfYyhPxaBykrbdhC3bQdxM2EWqcwiFYBklpovo6HtUnci\nM9rHwxf8lwkJuSw6DL/88x3tdiITwjz6phRbOmYfEYnCMvIzVPXfgQ6qWlA7VpFXgY8bG4gx9AZD\nCJCoaJxdE9lxez/OHLuZF3p+jBPh3pzxLJ9zGu5FVXTc5MGbXxAe6Xv1TGTepzvv0x1o5xOZ0P3/\nrAAyRaQPlvG+HrihTp/ZwF22/344UGpX5xPgdWCzqk4OviDIhw9wDbChsYEYQ28whACtqcabX0DP\nJwooeQJu5Hy7pYxedsFvE0LnqNKYeY/5a8tivnSgB88suuzo0piVlS0/PlXwheYRp6peEbkLyMJK\nr2CQS3MAAAuYSURBVHzDLsA00W6fBnyClVqZjZVe+VP78vOBm4H1IrLWPhdIo3xaRM7Act3sAm5v\nbCym8Ighcgi4I9yp7D+7G/ljvDx2/mxujS9kwSEHD0y63Q5y5rU7d0R7obmFRxI6pup5PW5pUt/P\nsv9iCo8YDC1OwB2xdTvxW7f///bOPTiq+orjn7OQQMAQAoGQLBSiFS21ijgCWgQ7VUR0Bu1MLZ0y\naoeKr6qdjo74aNXW97R2bEel2GKrMqCDOnWmo4GIQFERtbzCO2CsLEhAAkQQ8jr94/42rCG7CSGb\n/e3mfGbu7N3fvXf5cuCee3+Pcw595sE8iplHMQD9u8hwRKhXL0KFAzh4brBU9+aLy/h1/lbW1dYx\nY8M0aksHUPzSRhpramypbkt4/vLbHszRZyIxRUK23Z/dVCBkfk0+9664+ptFQmpqUq3W6GAaDx8+\ntlT3dSgjlzJGAZDPVmCrLdWNhwIZWDM2Mxx9jGPbOzqf398753jH9smXXafYcUyRkJKpxwqEgBUJ\nMYzEKGjm9fsyw9HHOLZ+5fCnOebYkkJscrQx+RyaGCRHi9Tl88D7VzGorPux5GjJqHJlGMlG6bDJ\nWJ/IDEdvdA7Nk6M9HyRHAxhOECBoD9QTQ7KyCfXrS92ZYSIX5fDOjCfpFerGnZFLWFZ2NuGldeSs\n30nD3i/R2lp7eHYGGWhjc/SGkUK0rpaG3VWEdlcxZClc//A4d+QQw/gAaHlZpmRlE+qbR90ZYXZe\n1IuhEyuZddqrFISyOeflOwgvraf3up007NlrD4gTJQNtZY7eMNIQraulYc8eQnv2MHg5NDwGNxA8\nJE5t6QERW1z+gjz+ctczfL9niFn7wzyxfLIrLu/Wr6dLgZCk0KFJzbzBHL1hdAWaFZf/3XOjmg4l\ntbh8TG78r78bJjIhm3GXruMP4YXc/L8rWVN6JuFlR+ixyZPc+ApYcXDDCGgqIVdSzL6HjjJ7xMt8\nLzuLp6pP57kllwQl5Fa5EnKpiHA0/CAmN37Wri8YVgY7mnLj7+NbPkYN2xu9YQTElpDLvwLujsm+\n2qYScjFJtTbe1Z97LgoSar13pJGb1kwjqzSPwhUHkMqIRbEanUjHpUDwCXP0RmqISao1fMb2byTU\nKna1FZJ6u4kQyskhNGggB84tJPJDpWLKLFbX1jOjfBqNpQUUfnCA0PadFkHalVBQW0dvtEh0HLKg\nP19/p4jIxdmsvP4pjmojt1ROYcPC4RT/x6NxSCN40Bw+TOP2Snpvr2T4azD5lmDcuoAtwBYr9tJV\nschYo0Wi45CRnWRFdjKsDK65/wJ3cC9D2At4Ng5pGElGuncnlJtL47Bidl+YBxP38dezX+K87G48\nu7+E+Q9OIm9VlX8rfXzR0YGYo89gJCubhjEjjltnfe/uC/n3ovMJL62n5+K1ts7aSApaX09DdTVU\nVzNwFfAM/Ibzm46fwof+9ZhUbdWNkV5oXS2h5auPW2cN9U1rrc29G74Sys1FBg9i33n9qZpYy6Nj\n3mBqbjVvH+7B7R//hLyFvRmwYH3H19XNwJcec/RG0oiG9382/ducdtl2ni1ZQN9Qd2bumsDCslGE\nl9TRa/0u6iM7M/LmMk6Oxpoa2FhD3sat5L0MLzCUFxgKQImrrdvxPQJFG7zrZ5w05uiNpBEN7x/8\naBVHH4XpTT2Ko5QkCO83PEeE7sO+1ZTv/sYJi7mz32bW19UyY8M0vi4dSNF7B2H1pvRbrWRpio1O\nIRqqPrSYzx8I8ezIuYzvCbMPFPPY8itcqLqbwDp0KNVqja6IKvUx+e4X05vFLt99HhXkUZHeQ4Jd\ncXmliMwBrgSqVPUs19YPeAUYRlCz8BpVrXbH7gGmE/SqblfVUtd+HvAPIIegTuId6nsdw1QQDVXf\nf4Dwj+ARRvKIO5TUUHXDSCUihHr0IDRoIDXnDCLygxCbfvwMFXVHuWHTNA4uHETR8hq6bYvQeOBg\n0noKCmgHvtGLyCTgaYKasX9T1cebHRd3fDJBzdjrVfW/ia5N5H/j6mjN14rIeOAr4MUYR/8ksE9V\nHxeRmUC+qt4tIiOAecBooBgoA4araoOIrARuBz4kcPR/VtW3WjOU1Yw1jCQSdbBFhRwcGQzF/GLC\nEq7ss4YbN/6Mw28XUvR+DaHtyXWwHcXJ1oztI/10bPeJbTp3Uf0rCf8sEekGbAEuBXYAHwE/VdUN\nMedMBm4jcPRjgKdVdUyia+P530RaW32jV9VlIjKsWfMU4GK3/09gCXC3a5+vqkeBT0WkAhgtIpVA\nH1Vd4f5yLwJXAa06eiMzacrDPjxMzsO7jp+oXVpPr3JLs5t0VGk8cuRY6cE3YCk5LGUsfdhGH7Z1\nucCxDpyMHQ1UqOp2ABGZT+AjN8ScM4XgJVqBFSLSV0SKCN7W410bz//Gpb1j9IWqusvtfwEUuv0w\nsCLmvB2urc7tN29vERGZAcxwX4+W6YLydupMFQXgoqTSg87XW0vwP+cLYDxuLQXAXGAu21r/hXSz\nMaSf5nTTC3DGyVxcQ3VpmS4oaOPpPUXk45jvs1V1dsz3MPB5zPcdEJMUKv454Vaujed/43LSk7Gq\nqiLSoa9bzlizAUTk45PpiqWCdNOcbnrBNHcG6aYXAs0nc72qTuooLZ1BW/1vqJ2/v9t1L3CfVa49\nAgyJOW+wa4u4/ebthmEYmUo8f9iWcxJdG8//xqW9jv5N4Dq3fx3wr5j2qSLSQ0RKgNOBla6bcVBE\nxrpZ5mtjrjEMw8hEPgJOF5ESEckGphL4yFjeBK6VgLHAAecvE10bz//GpS3LK+cRDPwXiMgO4AHg\nceBVEZkOfAZcA6Cq60XkVYIJg3rgVlWNzmzcwrHllW/R9onY2a2f4h3ppjnd9IJp7gzSTS94pFlV\n60Xkl0ApwRLJOc5H3uSOzyJYgTgZqCBYXvnzRNe6n27R/yai1eWVhmEYRnrT3qEbwzAMI00wR28Y\nhpHheO3oRWSSiGwWkQoXAeYdIlIpIutEZHV0aZeI9BORRSKy1X3mp1jjHBGpEpHymLa4GkXkHmfz\nzSJymUeaHxSRiLP1ahdV6IVmERkiIu+KyAYRWS8id7h2b+2cQLOXdhaRniKyUkTWOL0PuXZvbewN\nqurlRjABsQ04FcgG1gAjUq2rBZ2VQEGztieBmW5/JvBEijWOB0YB5a1pBEY4W/cASty/QTdPND8I\n3NnCuSnXDBQBo9x+LkH4+gif7ZxAs5d2BgQ4xe1nEaRTGeuzjX3ZfH6jbwofVtVaIBoCnA5MIQhN\nxn1elUItqOoyYF+z5ngam9JYqOqnBKsBRneK0BjiaI5HyjWr6i51yahUtQbYSBDd6K2dE2iOR0o1\na8BX7muW2xSPbewLPjv6eKHBvqFAmYh84lI3QDtClFNAojQWPtv9NhFZ64Z2ol10rzRLkBvqXII3\nzrSwczPN4KmdRaSbiKwmCBJapKppY+NU4rOjTxfGqepI4HLgVgmyfTahQR/S6zWs6aDR8RzBUN5I\nYBfwx9TKOR4ROQV4DfiVqh6MPearnVvQ7K2dVbXB3W+DCRImntXsuJc2TjU+O/q2hA+nHFWNuM8q\n4A2CruEJhyingBNNY5FyVHW3u9Ebgec51g33QrOIZBE4zLmq+rpr9trOLWn23c4AqrofeBeYhOc2\n9gGfHX1bwodTioj0FpHc6D4wESinHSHKKeCE0likQN9xRG9mx9UEtgYPNIuIAH8HNqrqUzGHvLVz\nPM2+2llEBohIX7efQ5CrfRMe29gbUj0bnGgjCA3eQjBbfl+q9bSg71SCWf01wPqoRqA/8A6wlaD4\nSr8U65xH0AWPpouenkgjcJ+z+Wbgco80vwSsA9YS3MRFvmgGxhEMGawFVrttss92TqDZSzsDZwOr\nnK5y4Leu3Vsb+7JZCgTDMIwMx+ehG8MwDKMDMEdvGIaR4ZijNwzDyHDM0RuGYWQ45ugNwzAyHHP0\nhmEYGY45esMwjAzn/2LiFEj5vBpuAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sigma = 3\n", "n = sigma * N\n", "m = 20\n", "x_grid = np.linspace(-0.5, 0.5, n, endpoint=False)\n", "shift_to_range = lambda x: -0.5 + (x + 0.5) % 1\n", "\n", "mat = phi(shift_to_range(x[:, None] - x_grid), n, m, sigma)\n", "plt.imshow(mat, aspect='auto')\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By design, each row of the matrix contains just a single nonzero clump of entries, of width approximately $2m$.\n", "We could make this basis function sum *far* more efficient if we construct this as a sparse rather than a dense matrix.\n", "Because each row has the same number of nonzero entries (that is, $2m$), this can be done quite efficiently using scipy's compressed sparse row (CSR) sparse matrix format.\n", "We will build the matrix from an array of values, an array of column indices, and an array of indices specifying the beginning and end of each sequential row within these arrays:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8lNX1/99nZhIgIQkJgSwTdoKAYt0AcQVBwGpdvrb+\nrFtt/RbRutXautRW/dpWahVFqyAuba0oUq1KrRoWZVFAdlkFwp7JHiAJCSSZmfP743kmDDEkgQzJ\nZHLfrxevzPPc+zxzE2bOc+85536OqCoGg8FgiFwcrT0Ag8FgMJxcjKE3GAyGCMcYeoPBYIhwjKE3\nGAyGCMcYeoPBYIhwjKE3GAyGCKfFDb2IjBeRLSKSLSIPtfT7GwwGQ3tDWjKPXkScwFbgUiAHWAH8\nWFU3tdggDAaDoZ3R0jP6YUC2qu5Q1WpgJnBVC4/BYDAY2hWuFn4/N7A36DgHGF63k4hMACYAOHGe\nHUN8y4zOYDCcFMThgCgXvpgoauKUlLgyYh1V7D6UBKUuXBU+pKoG9fugGU6GcvYXq2q3E71+3KhY\nLdnna1LfVeuqslR1/Im+V0vS0oa+SajqdGA6QLwk6XAZ3cojMhhaBunQAWe3ZCqGpOMZ6eLyMSv4\nU8oSiv3VTNx+Hbvn9KbXa9vwHyhFa6pbe7hNR4Fq+9+BI6cz6usrJ/428/S93Sd+NZTs87E8q2eT\n+jrTtiU31kdExgNTACfwmqpOqtN+I/Ag1m9dDtyhqt80dK2IJAHvAr2BXcB1qrq/oXG0tKH3AD2C\njjPscwYDEhWNo0sC3gFuci+MIWPsbqb2e5cUZzTfm3Ev7oVeYtfn4SssQqurIQJ1mrSqCm+Ohw45\nHvp+CpuBaxhmt+aSQS5Nm28aTgQF/PhDci87JvkSQTFJEZldJya5E7hYVfeLyGVYE9zhjVz7EDBf\nVSfZCS0PYT0sjklLG/oVQKaI9MEy8NcDN7TwGOpHBHFF4eyayI6J/Th77Cae7/FfAO7PuYzlc07D\nvbCKjps9+Ir3ta3ZVBtBa6rxFRUhRUW4vwKdBBO5AIC+LAXA25oDNEQ8ilKjIXuU1sYkAUQkEJOs\nNfSquiSo/zKOLHIauvYqYKTd7x/AAsLJ0KuqV0TuArKwliNvqOrGlhzDMVFFa6rx5hfQ8/ECih6H\nGznfbiyjF9b/hzE0YYgIrvQ0Dp2aTs7IKEaNWcvT6Z9T6ffxi13X8G1WJu4vD+FcstE8oA2NEqoZ\nPU2MSQZxG/BpE65NUdU8+3U+kNLYQFrcR6+qnwCftPT7GiIYVbyeXKI8ufSZA7segesYYTcW0YMi\nq9vJHofDifOUvpQMTab40sM8NfzfXNe5lP9WduS+5dfTZW4nur2/Ed/BCvAbB0w4oii+prsEk0Vk\nZdDxdDu+eNyIyCgsQ3/B8VynqioijQ44LIOxBkMokahoyyV3ez++N/ZbXuj5H6IQ7veMY+mc03Av\nrKbjJg/e/ILm+f39Pnybt9Fl8za6vAmv04fX7aa+rAUw/vU2gL/pU4JiVT2ngfYmxSRF5HTgNeAy\nVS1pwrUFIpKmqnkikgYUNjZQY+jrQVwuHAnx+Pq5yT8/jrhx+bwycAYDoqL5c8mpvPHFSAY9l4s/\nvxB/VVVEBgUjiVqX3BMF7H8Cbq51yZXTy/j+DUEo4Avd2q/RmKSI9AT+DdysqlubeO1s4CfAJPvn\nR40NxBj6elCvF1/JPijZR+py4Dm4v9YVAJksM4bBYGgOIrj69OLA2SnkjlbuvzCLO7vsZFW1j4nr\nb0KzupL29mb85eWot2W/bccxo2+QY8UkRWSi3T4N+D3QFXhZRAC8qnpOI/HMScAsEbkN2A1c19hY\nWlQC4UQwefRtH2d8PNozneJhiZSOreCFc95lfEwVM8sTeWTZNQycVI7m5OOvqDS+a0OzmKfvrWrE\nndIg3/tetH76SaPp8QC4M/Ka9V4tiZnRG046vrIy2FBG0gZIegOeYxDP2W2ZrI4ov7WzSwLaK53C\n4V04PK6Ml854h5Gd/LxemsofvrqCQU/vRz35+CsrjcsvDFE0lK6bsMEY+jDGGR+P9k6naHgilePK\nefHMdxjdycffy7rzxJIrGTTJGI1ww3egFA6UkvwNMB2e4nSestsGsDKiHmoRiYIvAr9KxtCHMb6y\nMlhXRtd10PVVeJohPG23hY3REMHZrzelZ3av9bXenbibZYd93LH+RshKIu2dzfhKy4xbxhD2WDtj\nIw9j6OvicOIc2I/iYV3ZN/YQzwx9j6tjD/JhRWceWP4jkuZ0JHlFCb5vtxvDBaCKL3snnbN3MuBf\n8DGJfEwiAN35FjAphYa2hOBrjthOmGIqTNXF74OCYpI2HiR2eQx/zh7HxupDXB5Tys+GLKH4bD8H\nB3TBER0FEnkfCIOhPaNAjUqT/rUlzIy+HmpTK1cAz9dNrfwaiMzlncHQLBxOHJ06Iu5U9p/TjTt+\n/x63xBcz/5CTu1bfQOycznT7ej+yKzdsdwdbefRty4g3BWPoDQZDw4gg0dE4uyVTeVo6notdrL75\nOQ74vdy584dsz+pL+uJDRG314N93AP/W7cRv3c6MtzOYYWt09WQ90DYmSP42NltvCm3K0IvLBWcM\nJP+8eKLHFTF98FsMiY5i8v5Mpi4Yg/tziF+Tj3fXHpOFYjCECtVa+eToHA99PoNrHz7Xbswng3wg\nMmIxZkbfygRkCTznxxMzvoBXBs1gUFQUfy4ZxGsLR1pGfm0+/rxm6pUYDIZ2iyL4IjB02WYMfUCW\nIHXKEpgCv+bc2raA39zIErRzAi6G5K62iyGKMZeu4em0BZT7vYx+5Te4Fx1xMRjJYkN9GNdNO0Jc\nLhyJiXgz08m9MBb32D1M6z+TFGc0jxUOZ8Wj5xC7Pg9/UbERNgsXAi4GTy7Rnlz6ZMH2R+Bae1LQ\nw64pEAkuBsPJQRGq1dnawwg5xtAfA/V6j1Q7WgL8+Ui1I1A6sMKsIAwhJ1AztvI0q2bs+DEr+VPq\n4iOBz8/60uuNbLMiOUlYG6aM68ZgaBgRHB064EjpxsEhaXhGOdh4/Yvs9lZz+9YbKM5yk/ZVBa5t\nOVaB6xZWJgx36gY9t3BkRRIIfJoVycnFBGMNkYMIjk6dKL9sCLmXwF0Xz+XexGzWV9fw8403U5OV\nTOrSMlj77fEZY1X8hw/j372Xjrv30u9juPJXQwGIZjfp7AaM+8QQnqgKPjUzekOkoIq/spLY978m\n833IIp4szgIgia3A1rat4ReU+10xxMr9/uamKRT7q5m4/Tp2z+lN+uJKorZ4rJWFcYMYbPxmRh8Z\nSFQ0jqQu1Axwk3thJzLHbeevfd6ni8PFQ3kXs/HRIcRszMNXVIxWV5tAa1skKPe7Q46Hvp/CNQ8N\nsxtzySAXMCsLw9FYwdjIM4uR9xs1Aa2pxldQiKOgkIzFcOhPcFttoLWKaFaaQGuoEMGV4abitDQ8\nI11cOeZrnkxZRrG/mgnZ/4+cOb3o+fo2M6s2hAUmGGuIOCQqGt/wweReGEOvsbuY1m8WyY5oHik4\nj//OHYp7oZeOn69r3qpGFe/eHDrszaHvp7ABuIbAzNqDG4+ZVRvCCp/JozdEElpTjePLtWR8Cb6n\n4Oe1qxovfe2i2cZpFZmIy4WjSwK+/m7yzu9M0rhcpp8yg16uaP5UfAYLHj2fzuvy8BcUtat9IqHe\nGSsi44EpWHVfX1PVSXXaBwJ/A84Cfquqz9jnTwHeDeraF/i9qj4vIo8DPweK7LZHVPWThsZhDH0r\nEviyeTMzyLsglu7jcpie+TYZrg48UXg2Sx4dTuw6synLEHrU68VXXALFJaQtA56Fuzm/tr0jy9ut\n+9IfoqwbEXECLwGXAjnAChGZraqbgrrtA+4Brg6+VlW3AGcE3ccDfBDU5bnAQ6EpGEPfigS+bFJc\nQvpS4C9wZ+2sGrMpyxD5iCCuKJzJSRwe7MZzUQeGjd3AlIzP+MWey1kzdxDuhYfpsNmDr2T/SY/j\nWKJmIZvRDwOyVXUHgIjMBK4Cag29qhYChSJyeQP3GQ1sV9XdJzoQY+gNTaZWFmKAm9wLYmplIV4u\nuZAP552Le4GX2PV5+AqLTLaSoWmoojXVePPyceXl02s+FDwG13MesJ+etmxFS014FKGm6RIIySKy\nMuh4uqpODzp2A3uDjnOA4ScwrOuBd+qcu1tEbgFWAr9S1f0N3cAYekOTOUoW4iuCZCG01qdvViBh\nhsOJM74z2ttNwYgEvOMOMPX0Gaw/3IM/L/4+6fOdJKwuQHML8B861O4fzqocz4apYlU952SOR0Si\ngSuBh4NOTwWexFqAPAk8C/ysofu0D0MfvC3/9DRG/uErHk1ex07vYSZsuZF9WemkfXUQZ7bHbMs3\nRBZ+H74DpbC2lG5rganwf/bGuAGsAMxegqORUG6Y8gA9go4z7HPHw2XAalUtCJwIfi0irwIfN3aT\n9mHo62zLX/afKK7gbAA6sIs0dgENfOBtuQBJT6H0zO7kjlayfzCNVdU+bl93M2QlkbK0FMeuXPzl\n5eZBYTC0UZTjmtE3xgogU0T6YBn464EbjvMeP6aO20ZE0lQ1zz68BitruUHajqG3l6D5Px6Md9wB\npp/+Fud2dPLSgR48s3j80UvQysrQvrctF0D2Tjpn72TAv+D7E61ZUXe+tbpgZkYGQyQQqmCsqnpF\n5C4gCyu98g1V3SgiE+32aSKSiuVnjwf8InIfMFhVy0QkFitj5/Y6t35aRM7AMju76mn/Dm3H0NtL\n0G5Tl8JUeMyekYNZghoaICirY/vEvpxz6Sae7/FfAO7bezmr5gwmfVEVHTd78Oab6mTtHUVCWnjE\nzm//pM65aUGv88EurPvdayuArvWcv/l4x9F2DL2hcQKxiG7JVAxJI2eUi3U/nkKBr5qJ2dfjmdOT\n9MUVuLbl4t+/v324mIKyOno9lk/RY3Bjbb74gRbP6jCENwrUGK0bQ1gTiEXYkgP9PoFrfh2QG8jB\nTQ5gVj6GlsHZNQlfPzd5F8QRPzafVwe+Rf+oDvypeAhvfnERA5/PxZ9fGGabAcXo0RvCH3G5cMTF\n4e/nJv+8eP7+y+cYEh3F5P2ZTF04Bvd8iF+Tb21tN+l0hpOIr2QflOwjbTkwGe7jvNq2/iwLy1WU\nErqdseGEMfQRhnq9+Pbvh5X7SVkJD75wZH9GWyiiLi4XjoR4fP3c5J8fR+dx+UwfOIMPys7kjS9G\nkvGFn7i1+WE4EzRECmZGH4SI9ADeBFKwHoTTVXWKiCRhifH0xooIXxfYtSUiDwO3YXkP7lHVrGaN\n3hBxqNdbOxNMXQ48B/czAoBMlgHh/aBqCRxxcUiPNIqHdqXk0sM8O/xfXB17kNkVMdy/4joG/KES\n3ZuHv6IS/MZRdzyoipnR18GLtfV2tYjEAatEZC5wKzBfVSeJyEPAQ8CDIjIYK4/0VCAdmCciA1S1\nzX4SawuYDHTjubATA8du46XeHxDjcDLsH/fjXlhDp425+IpLjCSAIWT4y8thUzmJmyDxHzCV/ky1\n2/qx1sRgmoEVjG2yBEKb4YQNvZ2wn2e/LheRzVjaDlcBI+1u/wAWAA/a52eqahWwU0SysUR/lp7o\nGFqb4AImPRZCxR/gVluUrLeRBDCcDAKFXIak4bn4SCGXAl81d2w3hVyaj6kZe0xEpDdwJvA1kBK0\naysfy7UD1kNgWdBlOfa5+u43AZgA0JGYUAzR0NapZ3fy/RdmcWeXndYO5WfuJWVJO9idHFzI5RNT\nyCXUWMFY46P/DiLSGXgfuM/ezVXbpqoqIsftr7AV4KYDxEtS+Pk7RHD16UXpmSnkXgJ3XTyXexOz\nWVvtZeLGm6jJSib9rc2RbXBamnp2J39MIh+TCEB3lpjdyYaQEMrCI+FCswy9iERhGfkZqvpv+3RB\nQItBRNKAQvt8KAR+wgNVvDt2EbtjF5nvQxbxZNlCUUlsBbYag3McBPKtg7NsBkRF8+eSU2szbWI/\nW2eybAwnnVDvjA0XmpN1I8DrwGZVnRzUNBv4CTDJ/vlR0Pm3RWQyVjA2E1h+ou9viByOlWUDRzJt\n/K00ttZCoqJxdEmg5hQ3uRfG0GfcTl7u+y+SHC7OfOuXZCzwErMhF19RsQn0hxhTHPxozgduBtaL\nyFr73CNYBn6WiNwG7AauA7DFfGZhVVfxAr9oyxk3BsPJRGuq8RUV4SgqIuNLqAmq6dsk7X9bBNDf\nJ4OCEQn4x+7nlSAhwLcfv9zo0NeDKtT4jaGvRVW/hGPuLBh9jGv+CPzxRN/TYAgVR+0gHhFP1Lhi\nlp05s3YHcfrnkLC6De8gDujQryml+xrg5aOFADuzzLgX68Fy3RhDbzBEBHV3EPMifJ+zQIRTYjYg\n6SnsH55G3ujuPHrBf7gtIZ/Tlt1Ih6x4ui87gOzOxVd20GxIikDMzlhD20QEiY7GmdyVytPS8Vwc\nxZhL1/B02gJu2n4VW+f0w73oEFFbPfj3HWjf+deq+CsqYNsO4rbtIG4mzCKVWaTiZiPQ/uIF7QmT\nXtmKOGJjEXcq+8/pxh2/f49b4ouZf8jJXatvIDarM92W70d25eI7WGFmWPWhilZV4fXkEu3JpU8W\nbH8EruVcoIAeWJXJzF/OYDCum1bDX1EBW7cTv3U7M97OYIat09+T9VZ7aw6uOdibgByp3b+Tk7++\nuoZbJ/+S1KVlOLZ7TE6+wdBChLBmbNjQJgx9xOFw4jgtk6LhiVSOK+fFM99hdCcffy/rzhNLruRf\n88aSuLKIlG1LUVUz0zYYWggr68Zo3RhCgd+Hf923dF0HXV+FpxnC03bTAFYCbdyNEogJdO9GxZA0\nhv5hJU90//roSldfVuLa6mk/la4MbYJQb5gSkfHAFKyasa+p6qQ67QOBvwFnAb9V1WeC2nYB5Vjm\nwKuq59jnj6kQfCyMoTeEnkBMwNZkWfdJsB5LG6505XDiiI05SiL46WHvc23nMiMRHEGEynUjIk7g\nJawC3znAChGZraqbgrrtA+4Brj7GbUapanGdcw9Rj0JwQ2Mxht4QGoLjDWel4Bmt3HvRHO7usoOh\nq36MzulKytJSHDvasOiY3/cdieDp9LVEmTASwZFAiLNuhgHZqroDQERmYqn41hp6VS0ECkXk8uO4\n77EUgo+JMfQtQcCVkZTI4cFuPBdHs/Snz1KDcs+eH/DNnIG4Fx2mw2YPvpL9bTO90RYd89saQAPe\ng0/pwqecRTe2WF1og7N4Q7vjOLJukkVkZdDxdFuQMYAb2Bt0nAMMp+koVt0OH/BK0L2PpRB8TIyh\nbwlsV4avqJgOa6vpVe7mB8Nu4pWBM/hn77n86X/y+UfSRWTE9yLumw6mTJ7B0EqoCt6mG/rigN/8\nJHGBqnpEpDswV0S+VdVFwR2aqhAcmYY+aIPQoVPTyRkZxagxa3k6/XMq/T5Gv/Ib3F8eIurblt0g\nFFwmL3Z8/eJdbdChYTBEFCF03TRLsVdVPfbPQhH5AMsVtIhjKwQfk8g09EEbhKI8ufSZA7segets\nw9qDJUB4uxEkKhrfuafiuagTA8ZuZ2qf94lzuPhN3kjmzT0T98IaOixYb5QLDYYQEmIf/QogU0T6\nYBn464EbmnKhiMQCDrt6XywwFvg/u/lYCsHHJDINfQSgNdU4Fq+hx2I49McjJQrhMH1s9cI2Z96D\ns1aGdWXf2ENMHjoLPw4eWP4jkuZ0JHlFiclaMbQqoTL0quoVkbuALKz0yjdsFd+Jdvs0EUkFVgLx\ngF9E7gMGA8nAB3YhJxfwtqp+Zt+6XoXghogsQ+9w4uwcS9G1p3Lg0kM8P2wml8ccZtbBBB7++n9I\nntuRriuK8W3ZYYxIa1A3a+Xv8BIDAOjHGiC8V1mGyCfUefSq+gnwSZ1z04Je54O91f9oyoDvHeOe\nJRxDIfhYRJah9/vwlZWR9LelJP0NXmAgL9hN/Y0hadvUkYvwjFbuvmhebQnH/332vrafvmkIC4wE\nQnvDXiHQI42iYUlHrRL6zf9p7QpBc/KNq+FkUzd9s04JR1Mz1hAKVMFrCo+0M+wVAhvLSNrIUasE\ns0IwGE6AwMospRubH0jljovncX/iNtZX1zBh001UZ3UjdYkl5Ofb3+Cu/pOGkSk2GAxtC3tVqr3T\nKRqeSMXYg/z1rLfx1CTy2JKrSZ3nInFlEerJx3/o8MlflQZWZjt3k/mL3cwjjnn2qiyRbcC2Vl2Z\nmeLghtahjr7KvrGHeWbYv7g69iD9v/ipyVQxNExgVbqu7CgRPYgQAb2TgBpDb2hx6tFXmUp/pmIy\nVQwngO06Kf/+kKP0iFZV+5i4/iY0qyspy0rhmy3tNqBtgrGG9k2Qf7XszFQ8l8AdI+dxaeymo/2r\nO3Pxl5a1W0MR1tiuk9j3vj5Kjwg4SpOovaJqfPSG9kxQBlKhnYH0Ym0G0kSS53YkfUURmpOPz7iQ\nDG0WwWeybgztFpOBFJ6I4OjQAUe3ZCqGpJEzysW6H085usjL4gpc23JNkZcmYnz0BoOh5TiOjBmv\nJ48Oe3Po9wlc8+sIKPLSSoRY6yZsMIbeYGgi4nLh6JKAd0AGeRfEkjpuL9P6v0OGqwO/LxzK+/NG\ncMoLe/EXFYdGZtpkzLQ8Gpkage3C0IvLhSMhHl8/N/nnx/H2vc8yICqaP5ecyhtfjCTjCz9xa/ON\nDryhQdTrxVdcghSXkL4EeBrurBWbg34sNTLTEYDJummjBOvApy6H+58zOvBtmoBLo2c6RecmcnDs\nQb694J/MKO/K75ZeTfe5USStLEZz8lpmE5AhYlATjDUYwoSAS2NDGV03QNfXYBxnAJDJKsC4NJqL\nuFw44uLw906nYEQCjNvHK6f/k7Ojnbx8oA8zHx9PwppCNLcA/6FDEbUKjqBfpRZj6A3tBomKxpHU\nhZoBbnIv6ES/8Tt4uc97dHG4eCjvYjb+7nRiNuTiKypu9wVd1Ou1tGb276f7GuBl+B1Da9s783XE\nPkxN1k17owFp3GGrr8eflWykcdsQWlONr6AQR0EhGYuh6im4rdbHXkU0K4wLr52jagx9+yNYGnfX\nXgZ9Hs+c3iP454jxTH3gRc4928lLB3rwzKLLSJ8vEbuUNRjaEya9sj3j9x21lH3s5bNrmwawHDB+\nYYOhQUSQ6Gic3ZKpPC2dU59cx6S0hRzwe7lz5w/Z/llf0r88RNRWD/59B9Ca6lYZZijnaCIyHpiC\nVUrwNVWdVKd9IPA34Czgt6r6jH2+B/AmkIKV3j9dVafYbY8DPweK7Ns8YleyOibG0LckgQ96cley\n7+zFRWPW8Yx7HpV+H3fvvppNWQNwLz5E1Let+0E3GE4KqmhVFd4cD9E5HrZ9Btdyrt2YTwb5QOtO\nmBTBH6KsGxFxAi8BlwI5wAoRma2qm4K67QPuAa6uc7kX+JWqrhaROGCViMwNuva5wEOhKRhD35IE\nPuieXHr/Npc9v4XrCKR6FtODYsCsDEKKw4kzvjPa243ncWXq6TM4v6ODaQfc/Hnx90mf7yRhdYHl\ncqusbO3RGsKAEE7ohwHZqroDQERmAlcBtYZeVQuBQhG5/KgxqOYBefbrchHZDLiDrz0emm3o7afW\nSsCjqleISBLwLtAb2AVcp6r77b4PA7dh2bJ7VDWrue//HRxOnIP6UzQsidKxlTw39N2jC4TP6UjS\nhxuMdnt7we/Dd6AU1paSdjX8n63UCDCAFYB5sBqCCG0w1g3sDTrOAYYf701EpDdwJvB10Om7ReQW\nLNv7q4CNPRahmNHfC2wG4u3jh4D5qjpJRB6yjx8UkcHA9cCpQDowT0QGqGpov2d+H76NW74jvAVH\nxLf8IX1Dw4lwVJ72476jcrQnLx5ngtuG1qPpH7VkEVkZdDxdVaeHcigi0hl4H7hPVcvs01OBJ7FG\n+iTwLPCzhu7TLEMvIhnA5cAfgfvt01cBI+3X/wAWAA/a52eqahWwU0SysZY2S5szBkPoCGi5+Pq7\nyTu/M13G5jH/tPf4Q/HpvPX5hWR87qfzujz8BUXNloo4Kk/7qqNztE1w29CaHMeMvlhVz2mg3QP0\nCDrOsM81CRGJwjLyM1T130fGpwVBfV4FPm7sXs2d0T8P/AaICzqXYvuXAPKxosZgLWOWBfXLsc99\nBxGZAEwA6EjMkSBm925svbsHV49ZxhPdvz5aivXLSlxbPUaKtRkEtFwoLiFtGfAsXIGVXdTfSEUY\n2gEK+P0hc92sADJFpA+Wgb8euKEpF4qIAK8Dm1V1cp22tCAbew2wobH7nbChF5ErgEJVXSUiI+vr\no6oqIsc97bOXP9MB4iVJa4OYe3Po+5sc1gHXYKRYwx673m3J/5x2VK3bDys688DyH5E0tyPJy0vw\nfbvdxEsM4YECIfLRq6pXRO4CsrDSK99Q1Y0iMtFunyYiqVh+9njALyL3AYOB04GbgfUista+ZSCN\n8mkROcMe7S7g9sbG0pwZ/fnAlSLyfaAjEC8ibwEFgSeOiKQBhXb/Zi1jDG0Qu95t4j+WHlXrFky9\n2zaFCI6YGMSdyv6zu5E/xsvOy15jwSEHd665gU5z4uj+9QFkdy6+soNt/qEdynCQbZg/qXNuWtDr\nfCxbWJcvoX4ZTVW9+XjHccKGXlUfBh4GsGf0D6jqTSLyF+AnwCT750f2JbOBt0VkMlYwNhNsZ2y4\nYEselF1xOrmX+Hngok/5RZe9LDvs4/Z1N+GYk2hJHuzMsTI5DIb2gCr+igrYup34rduJf+eIiFwP\n22sQUQkOERj3Pxl59JOAWSJyG7AbuA7AXrLMwsoD9QK/CHnGTXOxJQ86z1rGgFkwm67MpisAqWy2\nunCSZqFGetdgCAPEaN0cC1VdgJVdg6qWAKOP0e+PWBk6bRd7GVv6gyHkjfbx8AX/ZUJCLosOw51r\nb6RDVjzdlx3Av37r8RljI71rMIQHZkZvCCxj42YuI24mvE933qc7AG42AhG2jI1kRBBXFM7kJA4P\ncuP+4zYmZ3wKwH17L2fVnMGkL6qi46YcfCX7jSRFe0BBQ5d1EzYYQ29ov6iiNdV48/Jx5eVTMAJu\n5Hy78QA9WQJEZkqpIzbWCq6e0438S2t4YsRH3BJfzPxDTu5afQO9n6xBduXiO1jRDt2ExtAb2giB\nOrmemwZxlQ/AAAAgAElEQVTSaVwh0we/xalR0Tyz7xReWXgJ7vkQvzYf7649ZtdpO+So4OrbMIMM\nZtjJHz1Z375XpRH4dTCGPkIJ1MlNnbIEpsCva1UCIdOWzIjEmWpYIIKjQwccqd0pPyOVnEuEn41c\nyINdNzJy/Y84+FkqqUvKcW734C8tMxv8wg1j6A3tmVp9mr7pFJybgIwrYdqQt/i6sj+TF4/DPU9I\nWF2AP7+wfevTqOI/fBj/rj102rWHzA9hMR1ZzNl0Zged2QGY4HpYEsINU+GEMfSGJlOrT7NqP91X\nAS8d0agJ6NOYuWkDOJw4+/eu1y/+i1U30HlOZ7rP2thO/eLhQyTOT4yhNxhCiKNjR8tl871UPKMc\n3DJqEY8krye7porbt9xIp8diSczaSsKsMmZ4j/jFe7EeMLP8sMBk3RgM4YtERePsmkjVIDeeizpy\n5tjNvNDzY5wI9+aMx/NoJh03efDt249WV5+UqVuwy6b/R7CE6FphuE7sBIwxD3eOX50r/DGG3hAx\naE013vwCnPkF9PwCSp4ITpcsw8Uq41o6FrYAnWSkUnJOMr/63dtc17mUzyo7cM/K/0fCnFi6fb0P\n9uZFtmtJMcFYw0miTtFkz0gX48es5Dfdv2Dijh+xM6sP6YsridriwX+g1GzcMYQeW4COzeV02byN\n1//Zh9ftpj6sA9rLSkRMMNZwkqhTNLnPZ7AFuI0LgDwyrNKR7eSL1g4J6Bz1TmfX76L461lvM7qT\njzfLknls6VWkzo0icWUR6sk3OkctgZnRGyKaQP53SjcOnp5GzigHm/7fi+z2VjNhy42UzEkn/cuD\nOLPtlYXJ/w4NAZ2jdWX0/BE8zRCetpsGGJ2jlicCd4sZQ284QiD/e/deOu7eS///wJX3W+mTHdhF\nOrsAY3QaJeCKS0rk8GA3OSOjOe/SDUx2Z3HHniv4Jmsg7sWH6bDZYzR0wg2TR28wtDwBKYecnwwk\nfmw+rw58i/5RHfhT8RDe/OIi3F/4ifsmH+/uveGTAB1wxdkaOr3nQ+7v4HrOA/ZFtIZOJGCybgxt\nFnG5cCQm4s1MJ/fCWNxj9zCt/0xeLrmQD+edi3uBl9j1efiLiptd+DuUBKQc0iYvgclwH+fVtpk6\ntoaTQnh89EOKMfTtBPV68RUVIUVFuJcAf4aJXAAofVkKGINpoNbtVHPBaXgujmbEWMvlVINyz54f\n1LqdXF9uMC6nJiAi44EpWDVjX1PVSXXaBwJ/A84CfquqzzR2rYgkAe8CvbFqxl6nqvsbGocx9AZD\nEwm4kXz93ORdEFevK2ng87mW1k8YrYqOC9vt5Jq/il7zIff3AZcTBLud2uBv1mRC5boRESfwEnAp\nkAOsEJHZqropqNs+4B7g6uO49iFgvqpOEpGH7OMHGxqLMfSGiEdcLhxdEth76ykkj/PwyoC36eWK\n5smis3jn8/PJ+NxP5/V5ePfkNGicA24kSvaRtpx6XUlmVdTGUUIpgTAMyFbVHQAiMhO4CqucqvV2\nqoVAoYhcfhzXXgWMtPv9A6u6nzH0hrZPrXJmPzf5I+KJGlfMq6f+kyHRUUzen8lHj40hYXU+/oKi\n7yhnqteLr7iE9GeWwDNwd+1uWePnN9RD02f0ySKyMuh4uqpODzp2A3uDjnOA4U28d0PXpqhqnv06\nH0hp7GbG0BtaBYmKxpHUhZqBbuKezOGl3h8Q43Dym9xLWDD3DNwLa+i0MRdfcQlaXX1EOXPlflJW\nAi/Cg0HfmVi+NsbaEBKOw3VTrKrnnMShNIqqqkjjIzaG3tAqaE01voJCHAWFVFwEt3KB3XKI3uEY\nHBbB1acXpWel4LlEuffiOdzdZQdrq71M2HAT/qxkUpaWwjdbzEaytk7oAhAeoEfQcYZ9rrnXFohI\nmqrmiUgaUNjYzYyhN7Q4AZ+5d0AGeefH8uFdT5Ph6sDvC4fy/vxzyfjCR+y6MEv1VMW7YxexO3Yx\n4D34lC58ylkAJLMV2BrRAcp2Rej+I1cAmSLSB8tIXw/cEIJrZwM/ASbZPz9q7GbG0BtanIDPXIpL\nSF8Cd/7lgtq2fkE+c3G5cHbpgr+fm31PVDF98Fu1PvmpC8dYdW/X2H75yspW+m0MkYRo6LJuVNUr\nIncBWVgpkm+o6kYRmWi3TxORVGAlEA/4ReQ+YLCqltV3rX3rScAsEbkN2A1c19hYjKE3hC3BfvnE\ny4/2yZ9I3VuJisbRJYGaU9zkXhhDn3E7ebnvv0h2RPNQ/vmsf/QMYjbk4isqPml69YY2QAgLj6jq\nJ8Andc5NC3qdD3b1mSZca58vAUYfzziMoW8qAYXBnukUD0vkyUfeYHxMFTPLE3lk2TV0nxtN0opi\nNLcAf0WlURgMQ7SmGl9REY6iIjK+hJqn4Oe1sYEaolkRXnEBQ6tgJBDaMwGFwQ1lJG2A594YxHN2\nUyargbYv9nVU8e8RCTjGFbP8rJm8eKAvUxaOxT1fSFhjin8bIpwI/FgbQ2+o5TvFv/8K37cDjiEp\n/u1w4kyIx9/bepD4x+7nldPfYtXh3jyz6DLS5wsJawqtVZF5kBhagxD66MMJY+gNLYffZz1I9u+n\n+xrgZXjMrqcaeJC09VXRUQTr+w9JwzPKwcbrXyTHW8WEbTdQmJVB2pcVuLblGH3/cMIY+jAkSPv7\n0KluPBdHs+zWZ6lBuWv3layfcwruRYfp8K3R/ja0MHX0/ft9DFf+ytL3d7GHdPYAEfZwiwDEFB4J\nQ4K0v6Py8uk9D67/XUB/pKTFtb8DwleeWwYSNy6fVwbOYEBUNH8uOZU3vhhJxudhqJ9uMBgimrZv\n6MOMgPBV6nNL4Dm4nxG1bZlGV6VtUyfGwLh9TB0yg3M7Onlxfy8mLx7HoGeKTIyhrROB/23G0Bsi\nExEcMTFIegoHzupO3mgfO694la8O+7lj3Y24srqQsrQU2eXBV3awaemwDcQYwIozGDdMG8cEYw2G\nNoQq/ooK2LaDuG07iHsXxnEGAGlsBiKyBnTj2AHiinGn4xnl4JZRi3gkeT3ZNVX8/NubKM9KJfWr\ncmT15vYbHDaG3mA4eQQ0cHz93eSd35mkcblMP2UGvVzRDH73bjK+8NN5XZ4leRAuGjhtDTtA3Omj\n5fT/CJYQzRX2qiSWHcSyw+rWmmNsbSLwl293hl6iovGdeyq5F3ai37gdvNznPbo4XDyUdzFz5p2F\ne6GXDl+sM1vgW4GABg7FJaQtA549oh1vdOMNLYFgsm6+g4h0AV4DTsN6Dv4M2MIx6hmKyMPAbVgZ\nZfeoalZz3v9E0JpqHIvXkLEYqv4Et9Vuga+ijy2P21rmPVij3XNhJ+ZPeLpBjXbzIDIYQozx0dfL\nFOAzVf2hiEQDMcAj1FPPUEQGY0ltngqkA/NEZICqmviVTbBGe4+FcOsfwlyjvaUQwdmvN6VnpZA7\n2seDF37CxC4elh32MWHdTbiyupD6zqamB1UNhoYwhv4IIpIAXATcCqCq1UC1iByrnuFVwExVrQJ2\nikg2Vl3EpSc6BgjSNs+0tM2Tx3n4bPB7PFF4NrM+Pw/3F1Y9UOPXbcOo4sveSefsnQyYBR/QjQ/o\nBhwJrBrzbggZEWgimjOj7wMUAX8Tke8Bq4B7OXY9QzfYjlaLHPvcdxCRCcAEgI7ENDiIo7TNlwLP\nwJVYuw/7Gb9uRCAuF/6hp9Y+yOsr7h0zb515kBtCgnHdfPfas4C7VfVrEZmC5aappan1DOtiF9id\nDhAvSRH4Zz8ORHB06oSkp1iui0v8PHDRp/yiy16+t/zHOOYkkrK0FMfOnIh1XajXiyz9pvZBXl9x\n7wiMnxlaiwi0OM0x9DlAjqp+bR+/h2Xoj1XPsDn1E9svqlb1pCDXxWy6MpuupNpuC8W4LgxhRqDG\n7pkp5F4Cd108l3sTs1lfXcPPN95MTVYy6W9txl9eHl75+mqybo5CVfNFZK+InKKqW7Aqnmyy/9VX\nz3A28LaITMYKxmaCLVkY4dTVeZexJaw4+x1L532RrfO+2ui8GyKIoBq7me9DFvFk2ZLXSXaN3bCd\nnETg16+5WTd3AzPsjJsdwE8BB/XUM7RrJc7CehB4gV+0l4ybk67zbjC0ELWifTcNJGZ8Aa8MmsGg\nqCj+UjKY1xaOxP05xK/Nx59XgP/w4dYe7gkRSh+9iIzHyk50Aq+p6qQ67WK3fx+oBG5V1dUicgpW\nmnqAvsDvVfV5EXkc+DlWjBTgEbvs4DFplqFX1bXAOfU01VvPUFX/CPyxOe8ZChyxsUhGGvuGJlM4\npoYnR3zIjXElzKmM4u5VP6bvk9XInlx8Bysi0udtMJwotaJ9U5bAFPg159a2nUgd37AkRIZeRJzA\nS8ClWK7uFSIyW1U3BXW7DMu7kQkMB6YCw20vyRlB9/EAHwRd95yqPtPUsbS7nbGApYGyJZuELdkk\nvAVv0oM37fBBb9aZwF591CmikXOJgx9f8hW/67aa8Zt+SHGWm7SvTBENQxtHCaXrZhiQrao7AERk\nJlaaebChvwp4U1UVWCYiXQIxzqA+o4Htqrr7RAfSLg294QSoU0Sj/8ewAidXMpRodpOO9RmM1PWP\nREXjHzYYz0Ux9Bq7i2n9ZpHsiOaRgvP477yhuBd4iV2fi9eTa2IsbRjhuFw3ySKyMuh4up0xGMAN\n7A06zsGatdNIHzcQbOivB96pc93dInILsBL4VUB94FhEnqEXwZWawuFBbnIv6sDZYzfxfI//AnDf\n3stZNWcwfV/Zga94H+qtMV9KQ5PQmmrkq7VkfAW+p+DntdIZXvq2513LEchxGPpiVa3PdR0y7Pjn\nlcDDQaenAk9irT2eBJ7Fkp85JpFn6FXx5uXjysun5+dQ9DjcWJt3fYCeLDFfSEPosPc5OFK7U3pW\nCj98Iou7u+xgbbWXn6+/GZ3T1drnsCM3/FIJDfUTurlfU1LKG+tzGbBaVQtqhxf0WkReBT5ubCCR\nZ+jbGfWlbk4b8hZnRzt58UBf3ntsnEndPJnY+xz8dirhp+914VM7o6obW6wuRK5LKyIJ3VdkBZAp\nIn2wjPf1wA11+swG7rL998OB0jr++R9Tx21Tx4d/DbChsYEYQ9/GqS9183e2BARALF+bFYwhdIgg\n0dHUXHAanoujGTF2A5PdWfhQ7tlzBWvmDMK96DAdNnvw5hc0fr9wI4TqlarqFZG7gCys9Mo37DTz\niXb7NOATrNTKbKz0yp8GrheRWKyMndvr3PppETnDGi276mn/DqJhPsOLlyQdLvVma54cRHC506kY\nko5npIvLx6zgTylLKPZXM3H7dezJ6k364kocyzehNdUtNy6DwdAo8/S9Vc3xm8d076EDrru/SX2/\neen+Zr1XS9J+Z/QOJ87OsWjPdIqHJVI6toIXznmX8TFV9J13Jt3musj8ezFb/tCJaypG2Pn0ubjJ\nBSJy85zBYMBIIEQWfh++sjLYUEbSBkh6A55jEM8BmawGjF81nHDGx9f7UJ5Znsgjy65h4KRyNCcf\nf0Wl2eRmaBZGvdJgaCWO9VAG68FsTHsYErzJ7vQ0ckYd2WS321vNtS/+Ovw22YV2w1TY0OYMvTM+\nHnqkUTQsidKxlTw39F0ujznMrIMJPPz1/5A8pyNJH24wMzuDobWpu8nuP0c22QGkswQIw5WzMfSt\nj6+sDDaWkbQRkv4GLzCQF+y2/qwBWkGb3OHEGd8Z7e2m4NwEasaVsn7420wvTeepxZeTNt9Jl9WF\naG6BJTkc5gFwg6G9cpw7Y9sMbc7QhyV+H74DpbC2lG5rgWkwztIjYgArgDCctRgMDeCIi0N6plM8\nNIl9Yw8xeegsroyt5P2D8fxm+bWc8qcKdG9eRK6cxR95lt4Y+pYkaBflpl93596L5nB3lx2sqvYx\ncf1NaFZXUpaZXZSG1sdfXg4bt5C4ERL/Di8xgJfstv6sidyJi/HRG5pN0C7KAXfs4lNO8i7K4D0B\nF7u4bMxKJqV+RbG/mjt3/IidWX3o9do2KxBm9gQYDIBx3RhambpyB6/96nnOiHY1WKnKm+OhQ46H\nvp/CFuAahtl3yyODvMidmRkMJ4ox9IbWpK7cwSN/HVbbZipVhT/icsH3TiH/vAScY4uZftpbnBHt\nYsr+/ry4aAzu+ULcp+uNJlErY2b0hsjE4cTRqSPiTmX/Od3Iv7SGJ0Z8xC3xxQz66mZiszrTbfl+\nZJeputUc1OuFVRtJWQW8CI/w3Qd1BG7KbHsYQ2+ISPw+q+rW1u3Eb91O/NswgwxmkEFP1ltdWnmI\nJ4pERePokkDNKW6iniw8umDI3KG4F1oFQ3xFxWhVVWsP19DaqJFAMBjaHFpTja+oCEdREb5RpmDI\ncVGnfKRnlIPrLlnCY91X8f3N11KYlUHal2G2s7WZmDz6VkKioqgeMxTPSBfjx6zkN92/4M6dP2T7\nZ31JX1xJ1LZc/PsOmKwRgyHU1NnZ2u9jWIWDKxmKiz2kswcIkz0itnwyh0NwrwiMj4S9odeaGqI/\nW0Gfz6yskdu4AMgng3wgTD5kBkOABmbBYavvEgmohsz1Zmb0hvaFCK7ePSk7IxXPJXD7yM95IGkL\nG2uqmbDpJg5ldSdtSTms2WwMVoAGZsEQxvouBguzYcrQ7lDFu3M3MTt3k/kBfE4sn9sbvBLIJoHs\nSPxOGE4WgdVOWgplZ6Ry2RML+HXXTWyuqeH2zTdS+VkKaUvKcezw4C8ta7XJgwnGtnccThyxMUiP\nNLY+GsMzQ9/j6tiDfFjRmQeW/4ikOR1JXlFiaYCUl7f2aA2G8CKw2rEnDws/6MRCzgYgnu3Esz0s\n6uuG0tCLyHhgClYpwddUdVKddrHbv49VSvBWVV1tt+0CyrH+JN5ANSsRSQLeBXpjlRK8TlX3NzQO\nY+iPB7/PMuCbyul3A0ylP1Ptpn62cmZrf0gNrUMgjXP3/2bSa+yuY6Zxej25ERnsixiUkP3/iIgT\neAmr7msOsEJEZqvqpqBulwGZ9r/hwFT7Z4BRqlpc59YPAfNVdZKIPGQfP9jQWIyhNxhCQCCNM+Op\nInxPmTTOZhOID51pxYfuGDmP+xO3sb66hgmbbqI6qxupS8pg7bchd/GEMBg7DMhW1R0AIjITuAoI\nNvRXAW+qVbx7mYh0EZE0Vc1r4L5XASPt1/8AFhBphl6ionEkdaFmgJvcCzrRb/wOXu7zHl0cLh7K\nu5g5888i8697rA0w1dVm9mQwtEWC40P/hnnEMc+ODyWyDdh28uJDobuxG9gbdJzD0bP1Y/VxA3n2\nSOaJiA94RVWn231Sgh4E+UBKYwNpc4Zea6rxFRTiKCgkYzFUPRVIuQSoog9LzczJENE4OnbEkdqd\n8u+l4hnl4JZRi3gkeT3ZNVXcvuVGOj0Wh3NbTqsGNNsqx7lhKllEVgYdTw8yxqHgAlX1iEh3YK6I\nfKuqi4I7qKqKND7iNmfoDYaThXTogDO5K4cGpeEZGc1Fl67jGfc8qtTPXbuvZN+jvejwrQdfyf5W\n3aDnP3wY/649dNq1h/4fwRKiucIOanZiJ2BiRSeM6vEUHikOBEiPgQfoEXScYZ9rUh9VDfwsFJEP\nsFxBi4CCgHtHRNKAwsYGagy9wWCjVVV4PblEeXLpPQ/2PArXMcJuLcFJSUSvFoN1gTr8oYBpff9F\nF4eLR/Iv5LN55+Be4CVmQzvQBQqd62YFkCkifbCM9/XADXX6zAbusv33w4FS24DHAg5VLbdfjwX+\nL+ianwCT7J8fNTYQY+gNBgNwtC5Qzchgl2gNfU5GQNmWLXAmJXJ4sJuckdGcd+kGJruzqEG5+NVf\n4158mA6bW3YVFapgrKp6ReQuIAsrvfINVd0oIhPt9mnAJ1ipldlY6ZU/tS9PAT6wsi9xAW+r6md2\n2yRglojcBuwGrmv8dwrzYGW8JOlwGd3aw2i3BNwZlael4xkZxZgxa3g6bQEH/F7u2nkt27L60ev1\nbKM3ZAgL5ul7qxpxpzRIXEKGnnXBPU3qu+iTB5v1Xi2JmdEbGiTgzoj25NInC7Y/DNdyrt1aQAYF\nxh/cGIEdoand2fzLdH42agEPdt3I1ppqJnx7IwezUkn9qhzndg++kn2tPVpDeM99T4hmGXoR+SXw\nv1h/mvVYy44YjrFrS0QeBm7DihXdo6pZzXl/g6FNENgRumsPmffuYTEdWWwHTzuzg87sANpQAFUE\nV2oKhwe78VwczYixR9wt9+z5Ad/MGYh70WFcX25ok6s8I2oWhIi4gXuAwap6SERmYQUbBlPPri0R\nGWy3nwqkY+WHDlDVBj/fmadXMPk/S2tnPu43vzVpYwZDa6KKNy8fV14+veZD7u/hes6zG/fR0xZu\nO+n2UgRXrx6Un5FKzigHPxu1gHlDQnDbpmfdtBma67pxAZ1EpAZrJp8LPEz9u7auAmaqahWwU0Sy\nsdKFljb0BtvWxXJ/7xG1M582M+sxGOojyI1TXx789VMeqHXjmAlNI6jitdNMMz+ExXQMwT0xrptg\n7ET+Z4A9wCFgjqrOEZFj7dpyA8uCbhHYAfYdRGQCMAGgIzEnOkRDKBDBERNj1ZM9uxv5Y7w8dv5s\nbo0vZMEhB3euuYFe/+dF9uThKzto6sk2RpAbp748+FQjY9yqWBumIs/SN8d1k4g1S+8DHAD+JSI3\nBfdp6q6tuti7y6aDlXVzomM0hADVo+vJvgPvkM47pAPQgw2hrScrgqNTJyQ9hdIzu5M7Wrn/wiyG\nx2Rz+7qbYU4SKUtKcezKxV9ebma8htBjZIqPYgywU1WLAETk38B5HHvXVlN2iRnq4nDijO+Mv08G\nBSMS8I/dzyunv8W5HZ30mT2B9M8dJKwuQHML8B861Pa1fVTxV1ZC9k46Z+9kwL/gYxL5mKF051ur\nC2bGazh5mBn90ewBzhWRGCzXzWhgJVBB/bu2ZgNvi8hkrGBsJrC8Ge/fPvD78B0ohTWldF8DvAyP\n2cv8AfafrzWNnkRF4zvvVDwXdmLwuK282OtDYhxOfpN7CQvmnoF7YQ3RC9cbgTlD28D46I9GVb8W\nkfeA1Vgb5tZguVs6U8+uLXtH2CwsiU4v8IvGMm4M4Y/WVONYuIYeC6H8D3Br7W7KQ/S24+wR+L1p\n8zjj46FHGkXDkyi9tJLnhr7L5TGHmXUwgYeWXcvAp8rRnHz8FZXtLO5yXFo3bQazM9ZgaGkcTpyd\nY9He6RQNS+Tg2IO8cNZMxsbU0CfrNlLmRZG0shjNycN/6HA7M7TNo7k7Y+Pj3DrszDub1Hf+4kfN\nzliDwXAM/D58ZWWwroyu66Dra/Asp/IsMIBVQPPcceJy4UiIx9/XTd75ccSMK+C1QW8xICqav5QM\n5rUFIxk0OQ9/fiH+qirjUgtGTc1YQ3smMAvtmU7R8ETKLq3ghXPeZXxMFX3n/ozuc6NJWlFsBYXb\n3XI/vFCv15JSKNlH6grgebi/VoUTMvk6olU4m00EPviMoW8r2Nk32ttNwbkJVI8tY8O5M5hems5T\nX15O2nwnXVYVWoa2sjL0H9bALHRDGV03QNfX4TkG8RyQyWrAZMIYIoTIs/PG0LcZAtk3a0vpthaY\nBuM4A4ABrACMoTWEEIcTR2wMkpFKyTnJFF1axaRz3+e6zqX8t7Ijj//xp3Rbvg/25uE7WBFRKzjx\nR57vxhh6g6G9E7RJ7cDZKeSN9rHzilf56rCfO9bdiCsLTpl8mL/vOoPX7d3PSSyNzImFYjZMGQyG\nMMYu5FE1cgiei6MYc+kankr7gkq/j1/suoZvszJxf3mIqG89+AqCqs8FbVKLy95J3LtHVotpbAYi\n0vbVi6Bmw5TBYAhjVNGqKqKzVlq1Ax4JLoVYRA+KgNC6+I7K8Dkvjk7jCpk++C1OjYom89934P4c\n4tfm4y8oajs7t9vCGI8TY+jbA3XytivGHeSvZ73N6E4++mTdRurcKBJXFqGefJO3bTguvpPhMwV+\nbRemyeRrIMTlB1sCY+gNbZJ68rafZghPE5q8bUOY0Ugt1nv2/IDiR3u3eC3WNkGIffQiMh6YglUz\n9jVVnVSnXez272PVjL1VVVeLSA/gTSz1XwWmq+oU+5rHgZ+DvUSDR1T1k4bGYQy9wRBp2C6cQHGQ\n3vMh93dHFwdxsg+vw4mjU0ccfXuy75xkCi+t4ckRH3JjXAkDv7yZznM6023ZfmRPbsRl1jREqLJu\nRMQJvARciiXLvkJEZqvqpqBul2HpfmUCw4Gp9k8v8Cvb6McBq0RkbtC1z6nqM00dizH04UCgGEVa\nCmVnpOK5BP734gVcEf8Nt2++kcrPUkhbUo5jhylGYQghfp8lQb0lm4Qt2STMgDfpwZv0oBfrrS6t\nPMSWR0PpuhkGZKvqDgARmYkl7R5s6K8C3lRLi2aZiHQJqP8CeQCqWi4im7Hqd2ziBDCGPhwIFKPY\nuZuYnbvJ/AAW0omFnEs824lnu5HmjSTsmbRkpLFvaDKFY47MpE9ZfAtxc2Lp9nX7m0mHBcrxGPpk\nEVkZdDzdrqURwA3sDTrOwZqt00gfN7aRBxCR3sCZYAc9LO4WkVuwFIN/FajLfSyMoT8WgcpK6Skc\nOLs7eWO8PHr+x9yWkM+CQw5+9fTtdF92ANmdayorGY6PujPpt47MpHuzzurSykNs1zT9j198skXN\nRKQz8D5wn6qW2aenAk9iPZaeBJ4FftbQfYyhPxaBykrbdhC3bQdxM2EWqcwiFYBklpovo6HtUnci\nM9rHwxf8lwkJuSw6DL/88x3tdiITwjz6phRbOmYfEYnCMvIzVPXfgQ6qWlA7VpFXgY8bG4gx9AZD\nCJCoaJxdE9lxez/OHLuZF3p+jBPh3pzxLJ9zGu5FVXTc5MGbXxAe6Xv1TGTepzvv0x1o5xOZ0P3/\nrAAyRaQPlvG+HrihTp/ZwF22/344UGpX5xPgdWCzqk4OviDIhw9wDbChsYEYQ28whACtqcabX0DP\nJwooeQJu5Hy7pYxedsFvE0LnqNKYeY/5a8tivnSgB88suuzo0piVlS0/PlXwheYRp6peEbkLyMJK\nr2CQS3MAAAuYSURBVHzDLsA00W6fBnyClVqZjZVe+VP78vOBm4H1IrLWPhdIo3xaRM7Act3sAm5v\nbCym8Ighcgi4I9yp7D+7G/ljvDx2/mxujS9kwSEHD0y63Q5y5rU7d0R7obmFRxI6pup5PW5pUt/P\nsv9iCo8YDC1OwB2xdTvxW7f///bOPTiq+orjn7OQQMAQAoGQLBSiFS21ijgCWgQ7VUR0Bu1MLZ0y\naoeKr6qdjo74aNXW97R2bEel2GKrMqCDOnWmo4GIQFERtbzCO2CsLEhAAkQQ8jr94/42rCG7CSGb\n/e3mfGbu7N3fvXf5cuCee3+Pcw595sE8iplHMQD9u8hwRKhXL0KFAzh4brBU9+aLy/h1/lbW1dYx\nY8M0aksHUPzSRhpramypbkt4/vLbHszRZyIxRUK23Z/dVCBkfk0+9664+ptFQmpqUq3W6GAaDx8+\ntlT3dSgjlzJGAZDPVmCrLdWNhwIZWDM2Mxx9jGPbOzqf398753jH9smXXafYcUyRkJKpxwqEgBUJ\nMYzEKGjm9fsyw9HHOLZ+5fCnOebYkkJscrQx+RyaGCRHi9Tl88D7VzGorPux5GjJqHJlGMlG6bDJ\nWJ/IDEdvdA7Nk6M9HyRHAxhOECBoD9QTQ7KyCfXrS92ZYSIX5fDOjCfpFerGnZFLWFZ2NuGldeSs\n30nD3i/R2lp7eHYGGWhjc/SGkUK0rpaG3VWEdlcxZClc//A4d+QQw/gAaHlZpmRlE+qbR90ZYXZe\n1IuhEyuZddqrFISyOeflOwgvraf3up007NlrD4gTJQNtZY7eMNIQraulYc8eQnv2MHg5NDwGNxA8\nJE5t6QERW1z+gjz+ctczfL9niFn7wzyxfLIrLu/Wr6dLgZCk0KFJzbzBHL1hdAWaFZf/3XOjmg4l\ntbh8TG78r78bJjIhm3GXruMP4YXc/L8rWVN6JuFlR+ixyZPc+ApYcXDDCGgqIVdSzL6HjjJ7xMt8\nLzuLp6pP57kllwQl5Fa5EnKpiHA0/CAmN37Wri8YVgY7mnLj7+NbPkYN2xu9YQTElpDLvwLujsm+\n2qYScjFJtTbe1Z97LgoSar13pJGb1kwjqzSPwhUHkMqIRbEanUjHpUDwCXP0RmqISao1fMb2byTU\nKna1FZJ6u4kQyskhNGggB84tJPJDpWLKLFbX1jOjfBqNpQUUfnCA0PadFkHalVBQW0dvtEh0HLKg\nP19/p4jIxdmsvP4pjmojt1ROYcPC4RT/x6NxSCN40Bw+TOP2Snpvr2T4azD5lmDcuoAtwBYr9tJV\nschYo0Wi45CRnWRFdjKsDK65/wJ3cC9D2At4Ng5pGElGuncnlJtL47Bidl+YBxP38dezX+K87G48\nu7+E+Q9OIm9VlX8rfXzR0YGYo89gJCubhjEjjltnfe/uC/n3ovMJL62n5+K1ts7aSApaX09DdTVU\nVzNwFfAM/Ibzm46fwof+9ZhUbdWNkV5oXS2h5auPW2cN9U1rrc29G74Sys1FBg9i33n9qZpYy6Nj\n3mBqbjVvH+7B7R//hLyFvRmwYH3H19XNwJcec/RG0oiG9382/ducdtl2ni1ZQN9Qd2bumsDCslGE\nl9TRa/0u6iM7M/LmMk6Oxpoa2FhD3sat5L0MLzCUFxgKQImrrdvxPQJFG7zrZ5w05uiNpBEN7x/8\naBVHH4XpTT2Ko5QkCO83PEeE7sO+1ZTv/sYJi7mz32bW19UyY8M0vi4dSNF7B2H1pvRbrWRpio1O\nIRqqPrSYzx8I8ezIuYzvCbMPFPPY8itcqLqbwDp0KNVqja6IKvUx+e4X05vFLt99HhXkUZHeQ4Jd\ncXmliMwBrgSqVPUs19YPeAUYRlCz8BpVrXbH7gGmE/SqblfVUtd+HvAPIIegTuId6nsdw1QQDVXf\nf4Dwj+ARRvKIO5TUUHXDSCUihHr0IDRoIDXnDCLygxCbfvwMFXVHuWHTNA4uHETR8hq6bYvQeOBg\n0noKCmgHvtGLyCTgaYKasX9T1cebHRd3fDJBzdjrVfW/ia5N5H/j6mjN14rIeOAr4MUYR/8ksE9V\nHxeRmUC+qt4tIiOAecBooBgoA4araoOIrARuBz4kcPR/VtW3WjOU1Yw1jCQSdbBFhRwcGQzF/GLC\nEq7ss4YbN/6Mw28XUvR+DaHtyXWwHcXJ1oztI/10bPeJbTp3Uf0rCf8sEekGbAEuBXYAHwE/VdUN\nMedMBm4jcPRjgKdVdUyia+P530RaW32jV9VlIjKsWfMU4GK3/09gCXC3a5+vqkeBT0WkAhgtIpVA\nH1Vd4f5yLwJXAa06eiMzacrDPjxMzsO7jp+oXVpPr3JLs5t0VGk8cuRY6cE3YCk5LGUsfdhGH7Z1\nucCxDpyMHQ1UqOp2ABGZT+AjN8ScM4XgJVqBFSLSV0SKCN7W410bz//Gpb1j9IWqusvtfwEUuv0w\nsCLmvB2urc7tN29vERGZAcxwX4+W6YLydupMFQXgoqTSg87XW0vwP+cLYDxuLQXAXGAu21r/hXSz\nMaSf5nTTC3DGyVxcQ3VpmS4oaOPpPUXk45jvs1V1dsz3MPB5zPcdEJMUKv454Vaujed/43LSk7Gq\nqiLSoa9bzlizAUTk45PpiqWCdNOcbnrBNHcG6aYXAs0nc72qTuooLZ1BW/1vqJ2/v9t1L3CfVa49\nAgyJOW+wa4u4/ebthmEYmUo8f9iWcxJdG8//xqW9jv5N4Dq3fx3wr5j2qSLSQ0RKgNOBla6bcVBE\nxrpZ5mtjrjEMw8hEPgJOF5ESEckGphL4yFjeBK6VgLHAAecvE10bz//GpS3LK+cRDPwXiMgO4AHg\nceBVEZkOfAZcA6Cq60XkVYIJg3rgVlWNzmzcwrHllW/R9onY2a2f4h3ppjnd9IJp7gzSTS94pFlV\n60Xkl0ApwRLJOc5H3uSOzyJYgTgZqCBYXvnzRNe6n27R/yai1eWVhmEYRnrT3qEbwzAMI00wR28Y\nhpHheO3oRWSSiGwWkQoXAeYdIlIpIutEZHV0aZeI9BORRSKy1X3mp1jjHBGpEpHymLa4GkXkHmfz\nzSJymUeaHxSRiLP1ahdV6IVmERkiIu+KyAYRWS8id7h2b+2cQLOXdhaRniKyUkTWOL0PuXZvbewN\nqurlRjABsQ04FcgG1gAjUq2rBZ2VQEGztieBmW5/JvBEijWOB0YB5a1pBEY4W/cASty/QTdPND8I\n3NnCuSnXDBQBo9x+LkH4+gif7ZxAs5d2BgQ4xe1nEaRTGeuzjX3ZfH6jbwofVtVaIBoCnA5MIQhN\nxn1elUItqOoyYF+z5ngam9JYqOqnBKsBRneK0BjiaI5HyjWr6i51yahUtQbYSBDd6K2dE2iOR0o1\na8BX7muW2xSPbewLPjv6eKHBvqFAmYh84lI3QDtClFNAojQWPtv9NhFZ64Z2ol10rzRLkBvqXII3\nzrSwczPN4KmdRaSbiKwmCBJapKppY+NU4rOjTxfGqepI4HLgVgmyfTahQR/S6zWs6aDR8RzBUN5I\nYBfwx9TKOR4ROQV4DfiVqh6MPearnVvQ7K2dVbXB3W+DCRImntXsuJc2TjU+O/q2hA+nHFWNuM8q\n4A2CruEJhyingBNNY5FyVHW3u9Ebgec51g33QrOIZBE4zLmq+rpr9trOLWn23c4AqrofeBeYhOc2\n9gGfHX1bwodTioj0FpHc6D4wESinHSHKKeCE0likQN9xRG9mx9UEtgYPNIuIAH8HNqrqUzGHvLVz\nPM2+2llEBohIX7efQ5CrfRMe29gbUj0bnGgjCA3eQjBbfl+q9bSg71SCWf01wPqoRqA/8A6wlaD4\nSr8U65xH0AWPpouenkgjcJ+z+Wbgco80vwSsA9YS3MRFvmgGxhEMGawFVrttss92TqDZSzsDZwOr\nnK5y4Leu3Vsb+7JZCgTDMIwMx+ehG8MwDKMDMEdvGIaR4ZijNwzDyHDM0RuGYWQ45ugNwzAyHHP0\nhmEYGY45esMwjAzn/2LiFEj5vBpuAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from scipy.sparse import csr_matrix\n", "\n", "col_ind = np.floor(n * x[:, np.newaxis]).astype(int) + np.arange(-m, m)\n", "vals = phi(shift_to_range(x[:, None] - col_ind / n), n, m, sigma)\n", "col_ind = (col_ind + n // 2) % n\n", "row_ptr = np.arange(len(x) + 1) * col_ind.shape[1]\n", "spmat = csr_matrix((vals.ravel(), col_ind.ravel(), row_ptr), shape=(len(x), n))\n", "\n", "plt.imshow(spmat.toarray(), aspect='auto')\n", "plt.colorbar()\n", "np.allclose(spmat.toarray(), mat)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "from numpy.fft import fft, ifft, fftshift, ifftshift\n", "\n", "\n", "def nfft3(x, f, N, sigma=2):\n", " \"\"\"Alg 3 from https://www-user.tu-chemnitz.de/~potts/paper/nfft3.pdf\"\"\"\n", " n = N * sigma # size of oversampled grid\n", " m = 20 # magic number: we'll set this more carefully later\n", " \n", " # 1. Express f(x) in terms of basis functions phi\n", " shift_to_range = lambda x: -0.5 + (x + 0.5) % 1\n", " col_ind = np.floor(n * x[:, np.newaxis]).astype(int) + np.arange(-m, m)\n", " vals = phi(shift_to_range(x[:, None] - col_ind / n), n, m, sigma)\n", " col_ind = (col_ind + n // 2) % n\n", " row_ptr = np.arange(len(x) + 1) * col_ind.shape[1]\n", " mat = csr_matrix((vals.ravel(), col_ind.ravel(), row_ptr), shape=(len(x), n))\n", " g = mat.T.dot(f)\n", " \n", " # 2. Compute the Fourier transform of g on the oversampled grid\n", " k = -(N // 2) + np.arange(N)\n", " g_k_n = fftshift(ifft(ifftshift(g)))\n", " g_k = n * g_k_n[(n - N) // 2: (n + N) // 2]\n", " \n", " # 3. Divide by the Fourier transform of the convolution kernel\n", " f_k = g_k / phi_hat(k, n, m, sigma)\n", " \n", " return f_k" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = -0.5 + np.random.rand(1000)\n", "f = np.sin(10 * 2 * np.pi * x)\n", "N = 100\n", "\n", "np.allclose(ndft(x, f, N),\n", " nfft3(x, f, N))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Choosing ``m``\n", "\n", "Finally, we should be a bit more careful about choosing a suitable value of ``m`` for our problem.\n", "The paper offers a way to estimate ``m`` from a desired error tolerance for the result of the computation." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def C_phi(m, sigma):\n", " return 4 * np.exp(-m * np.pi * (1 - 1. / (2 * sigma - 1)))\n", "\n", "def m_from_C_phi(C, sigma):\n", " return np.ceil(-np.log(0.25 * C) / (np.pi * (1 - 1 / (2 * sigma - 1))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's add this to our nfft function" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "from numpy.fft import fft, ifft, fftshift, ifftshift\n", "\n", "\n", "def nfft(x, f, N, sigma=2, tol=1E-8):\n", " \"\"\"Alg 3 from https://www-user.tu-chemnitz.de/~potts/paper/nfft3.pdf\"\"\"\n", " n = N * sigma # size of oversampled grid\n", " m = m_from_C_phi(tol / N, sigma)\n", " \n", " # 1. Express f(x) in terms of basis functions phi\n", " shift_to_range = lambda x: -0.5 + (x + 0.5) % 1\n", " col_ind = np.floor(n * x[:, np.newaxis]).astype(int) + np.arange(-m, m)\n", " vals = phi(shift_to_range(x[:, None] - col_ind / n), n, m, sigma)\n", " col_ind = (col_ind + n // 2) % n\n", " indptr = np.arange(len(x) + 1) * col_ind.shape[1]\n", " mat = csr_matrix((vals.ravel(), col_ind.ravel(), indptr), shape=(len(x), n))\n", " g = mat.T.dot(f)\n", " \n", " # 2. Compute the Fourier transform of g on the oversampled grid\n", " k = -(N // 2) + np.arange(N)\n", " g_k_n = fftshift(ifft(ifftshift(g)))\n", " g_k = n * g_k_n[(n - N) // 2: (n + N) // 2]\n", " \n", " # 3. Divide by the Fourier transform of the convolution kernel\n", " f_k = g_k / phi_hat(k, n, m, sigma)\n", " \n", " return f_k" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = -0.5 + np.random.rand(1000)\n", "f = np.sin(10 * 2 * np.pi * x)\n", "N = 100\n", "\n", "np.allclose(ndft(x, f, N),\n", " nfft(x, f, N))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from pynfft import NFFT\n", "\n", "def cnfft(x, f, N):\n", " \"\"\"Compute the nfft with pynfft\"\"\"\n", " plan = NFFT(N, len(x))\n", " plan.x = x\n", " plan.precompute()\n", " plan.f = f\n", " # need to return a copy because of a\n", " # reference counting bug in pynfft\n", " return plan.adjoint().copy()" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(cnfft(x, f, N),\n", " nfft(x, f, N))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "fast nfft:\n", "10 loops, best of 3: 23.4 ms per loop\n", "\n", "wrapped C-nfft/pynfft package:\n", "10 loops, best of 3: 40 ms per loop\n" ] } ], "source": [ "x = -0.5 + np.random.rand(10000)\n", "f = np.sin(10 * 2 * np.pi * x)\n", "N = 10000\n", "\n", "#print(\"direct ndft:\")\n", "#%timeit ndft(x, f, N)\n", "#print()\n", "print(\"fast nfft:\")\n", "%timeit nfft(x, f, N)\n", "print()\n", "print(\"wrapped C-nfft/pynfft package:\")\n", "%timeit cnfft(x, f, N)" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Fast Template Periodogram", "language": "python", "name": "ftp" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" } }, "nbformat": 4, "nbformat_minor": 2 }