{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Filter Example\n", "\n", "This example demonstrates the connection between MKS and signal\n", "processing for a 1D filter. It shows that the filter is in fact the\n", "same as the influence coefficients and, thus, applying the `predict`\n", "method provided by the `MKSLocalizationnModel` is in essence just applying a filter." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "if __import__('pyfftw'):\n", " import pyfftw\n", "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we construct a filter, $F$, such that\n", "\n", "$$F\\left(x\\right) = e^{-|x|} \\cos{\\left(2\\pi x\\right)} $$\n", "\n", "We want to show that, if $F$ is used to generate sample calibration\n", "data for the MKS, then the calculated influence coefficients are in\n", "fact just $F$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xmcm9V56PHfo12zejyeGRt7xhi84BBsszkkEBgCYbtJ\nTLgkAbJBW0ra0OQ2aS+QNsHuTbNv7eW2CYFSkiYhadIWspSYQIYASQj76h28DfYsnn20S+f+IWnQ\naKRXr0aaTXq+H/xBenfb8qMzz3nOOWKMQSmlVHVxzPUDKKWUmn0a/JVSqgpp8FdKqSqkwV8ppaqQ\nBn+llKpCGvyVUqoKlSX4i8idItIjIs9bHPOPIrJHRJ4VkU3luK9SSqnpKVfL/y7g4nw7ReRS4ERj\nzBrgBuCbZbqvUkqpaShL8DfGPAoMWhyyBfhO6tjHgUYRaSvHvZVSShVvtnL+y4FDGe+7U9uUUkrN\nAe3wVUqpKuSapft0A+0Z71ektk0hIjrZkFJKFckYI8UcX86Wv6R+5XIf8CEAETkLGDLG9OS7kDFG\nf5Xh16233jrnz1Apv46//UqWvfvUOX+OSvqln8/y/ZqOsrT8ReT7QCfQLCIHgVsBTzKOm9uNMb8Q\nkctEZC8wDlxXjvsqNVsiiRjheGSuH0OpsilL8DfGXGPjmBvLcS+llFKl0w7fCtbZ2TnXj1BR/Ota\n5/oRKop+PueWBv8Kpv+4yiMYCyMIzjWLp51fVVPp53NuafBXqoDxaJAmXz1OcRCOR+f6cZQqi9kq\n9VRqwRqLBKl1+zAYxqNBfC7PXD+SUiXTlr9SBYzHQtS5/dS5/YxHQ3P9OEqVhbb8lSpgPBKk1u3H\nYBiLBuf6cZQqCw3+ShUwFp2c9lGqEmjwV6qA8WhoouWvaR9VKTT4K1XAeDRIndsPmvZRFUSDv1IF\nJFv+Pm35q4qiwV+pApI5f7/m/FVF0eCvVAHJtE8NYBiLaPBXlUHr/JUqYDwaos7jp8bt07SPqhja\n8leqgIm0jzH0BYbm+nGUKgsN/koVMJ6e3sEktNpHVQwN/koVEIxF8Lu8xBMJAjFN+6jKUJacv4hc\nIiI7RWS3iNyUY3+DiNwnIs+KyAsicm057qvUbAgnonidbvwuD+GYzuqpKkPJwV9EHMBtwMXAycDV\nInJS1mEfBV4yxmwCzge+KiL6U4daECLxZPD3uTyEdClHVSHK0fLfDOwxxhwwxkSBe4AtWccYoD71\nuh44ZoyJleHeSs24cDyK1+nB5/QQimnwV5WhHK3v5cChjPeHSX4hZLoNuE9EXgPqgPeV4b5KzYpI\nPIrHmfynoi1/VSlmK/VyMfCMMeZtInIi8ICIbDDGjOU6eOvWrROvOzs7dbk3NaeSwd+NQxyEYuG5\nfhyl6Orqoqurq6RrlCP4dwMdGe9XpLZlug74PIAxZp+IvAqcBDyZ64KZwV+puRZKpX1cEte0j5oX\nshvF27ZtK/oa5cj5PwGsFpGVIuIBrgLuyzrmAHAhgIi0AWuBV8pwb6VmXLLD16UdvqqilNzyN8bE\nReRGYDvJL5M7jTE7ROSG5G5zO/BZ4F9F5PnUaf/bGDNQ6r2Vmg2ReBSPww2gLX9VMcqS8zfG3A+s\ny9r2rYzXR0jm/ZVacMLxKF6XB0E7fFXl0Fp7pSwkTIJoIobHkfynYowhlojjcjjn+MmUKo3O6qmU\nhUg8htfpRkQQkWTeX1M/qgJo8FfKQjhV5pnmc3oIxbXcUy18GvyVspDs7H09O+pzebXlryqCBn+l\nLITjEbwuz8R7TfuoSqHBXykL4YwyT0infTT4q4VPg79SFpIdvplpHw9BbfmrCqDBXykL4XgErzMj\n7aMze6oKocFfKQuR7Gofl1b7qMqgwV8pC+F4bGI6ZwCv00M4rqt5qYVPg79SFsLxCL6MtI/X6Sai\nwV9VAA3+SlnITvt4nC4N/qoiaPBXykIkEZsU/L1ODyEN/qoCaPBXykI4FsE7qeWvaR9VGTT4K2Uh\nnIhmtfw17aMqgwZ/pSxE4lF8mS1/h5tIPDaHT6RUeZQl+IvIJSKyU0R2i8hNeY7pFJFnRORFEfl1\nOe6r1EzLntXT43RrqaeqCCUv5iIiDuA24ALgNeAJEbnXGLMz45hG4P8BFxljukVkSan3VWo2TK32\ncTMQGpnDJ1KqPMrR8t8M7DHGHDDGRIF7gC1Zx1wD/MQY0w1gjOkvw32VmnGhWHRSh69PO3xVhShH\n8F8OHMp4fzi1LdNaYLGI/FpEnhCRD5bhvkrNuEgiOrXaJ6E5f7XwzdYavi7gNOBtQC3wOxH5nTFm\nb66Dt27dOvG6s7OTzs7OWXhEpabKlfbRnL+aa11dXXR1dZV0jXIE/26gI+P9itS2TIeBfmNMCAiJ\nyG+AjUDB4K/UXArHc7T8NfirOZbdKN62bVvR1yhH2ucJYLWIrBQRD3AVcF/WMfcC54iIU0RqgDcB\nO8pwb6VmVHa1j9b5q0pRcsvfGBMXkRuB7SS/TO40xuwQkRuSu83txpidIvJL4HkgDtxujHm51Hsr\nNdMiWSt5eRya9lGVoSw5f2PM/cC6rG3fynr/FeAr5bifUrMlHI/ic2naR1UeHeGrlIXslr9XO3xV\nhdDgr5SFSDyK15U1n7+WeqoKoMFfKQuheBSP4/XsqKZ9VKXQ4K+UhUiOUs9wXBdwVwufBn+lLGSn\nfZLBX9M+auHT4K+UhXBW2kfr/FWl0OCvlIXs6R28To8Gf1URNPgrZSEUj+BzTk77aPBXlUCDv1IW\nIvHJC7h7HC6iiTjGmDl8KqVKp8FfKQuRRBSP8/Wcv4jgdjh1oJda8DT4K5VHPBEnnkjgdkyeBUVT\nP6oSaPBXKo9IPIbX6UZEJm1PLuiiwV8tbBr8lcojnJXySfNqrb+qABr8lcojHIvgzaj0SfNq2kdV\nAA3+SuURSUyu9EnTKR5UJdDgr1Qe2Us4piU7fDXtoxa2sgR/EblERHaKyG4RucniuDNFJCoiV5Tj\nvkrNpHA8kjv4O1xa6qkWvJKDv4g4gNuAi4GTgatF5KQ8x30B+GWp91RqNmRP7ZCmOX9VCcrR8t8M\n7DHGHDDGRIF7gC05jvsL4MdAbxnuqdSMS5d6ZtP5fVQlKEfwXw4cynh/OLVtgogcB1xujPlnYHLR\ntFLzVDgeydPh69I6f7XglWUBdxu+AWT2BVh+AWzdunXidWdnJ52dnTPyUEpZCedJ+3icbkIxDf5q\n7nR1ddHV1VXSNcoR/LuBjoz3K1LbMp0B3CPJoZJLgEtFJGqMuS/XBTODv1JzJRyP4stX7aMtfzWH\nshvF27ZtK/oa5Qj+TwCrRWQlcAS4Crg68wBjzAnp1yJyF/DTfIFfqZkQjcdw5xita6XcHb7TeQal\nZkrJOX9jTBy4EdgOvATcY4zZISI3iMif5jql1HsqVYyj4wMcf8d7uHfvI0WdF4lH8Thyl3oWW+f/\nvR3bOf6O9zAYGi3qPKVmSlnq/I0x9xtj1hlj1hhjvpDa9i1jzO05jv0jY8x/lOO+Stnx81d+S53b\nz493dxV1XigexevKN8K3uJb/v+/+NfWeGrbv/0NR5yk1U3SEr6p4T/fu5i9Pfy9P9uwsahGWfC1/\nr9NT1PQO8USc5/v2ceOmK3iqd7ft85SaSRr8VcXbcWw/b12xEZfDRW9g0PZ5kTzTOxSb8z86PkCT\nt54zl65nx7H9ts9TaiZp8FcVzRjDwdFeVtYv5YTG43h1+Ijtc/PP7VNczn//yFFWNixlZcNSDo3q\nGEc1P2jwVxWtLzhEjctLncfPstpmjgYGbJ+br9qn2Jz/wdEeOhraaK1ZxGgkQCAasn2uUjNFg7+q\naAdGeuhoWApAW+1ijo7bD/5hiw7fYur89w8nW/4OcbCivoWDoz22z1VqpmjwVxXt4MhROupbAVhW\nu5ieIlr+4XgUb84O3+Jy/q+N9bG8bgkAHQ1LOTiiwV/NPQ3+qqJ1j/WzvK4FgLaa4lr+VoO8ikn7\nHAuNsMTfCEB7XYvm/dW8oMFfVbSB0AjNqcDbVruYnmLTPnkXc7Ef/PuDwyzxLwKgtXYxvYEh2+cq\nNVM0+KuKNhAaZbGvHoCltYs5Gjhm+9y8E7s5im/5N/saAGj1L6IvaL/cVKmZosFfVbTB0AiLU4G3\n2dfIQBHTK1jX+dsr9TTGMBAcZrE/FfxrmujTlr+aBzT4q4o2kBH86z1+AtEQUZuBO5n28UzZnqzz\nt9fyH40EcDtd+F1eIBn8e4oYaKbUTNHgrypaZtrHIQ4avLWMRMZtnZvs8J06C2cxpZ7JlE/jxPvW\nmib6gtryV3NPg7+qaJktf4BF3joGw2O2zo0kSq/26Q8O0+x//f7NvgYGQiPEE3Fb5ys1UzT4q4oV\njkcJx6PUe2omtjV5621PqxyK5Uv72K/2GQyN0pTx5eN2umjw1BbV96DUTNDgrypWMvDWk1xALmmR\nr56hsL3AG0lE8eZK+xQxn/9oZJyGjC8fgMW+egZtPoNSM6UswV9ELhGRnSKyW0RuyrH/GhF5LvXr\nURE5pRz3VcpKMuVTP2nbIm8dQ3bTPhZz+xTT4dvgqZ20bbGvgYHgiK3zlZopJQd/EXEAtwEXAycD\nV4vISVmHvQKca4zZCHwW+Hap91WqkIHQyKSUC0CTt46hkL3gn6/ax+t0E7bZ4TsSCUxp+Tdpy1/N\nA+Vo+W8G9hhjDhhjosA9wJbMA4wxvzfGDKfe/h5YXob7KmUps9InbZG3iLSPVbVPES3/+uzg761n\nIKQtfzW3yhH8lwOHMt4fxjq4/wnw32W4r1KWsit9ABq9tbbTPuF4FF/eDl97Of+RyDj1OdI+upav\nmmtTmzUzSETOB64DzpnN+6rqNJgj+Nd5ahiLBm2dn396BxfRRIyESeAQ6/bTaJ60z7HgcJ4zlJod\n5Qj+3UBHxvsVqW2TiMgG4HbgEmOM5RDHrVu3Trzu7Oyks7OzDI+pqs1AaJSVDW2TttW7/YxFCgf/\nWKoO3+VwTtknIhMVPz7X1J8MMuVK+yz21bNn8HDBZ1Aqn66uLrq6ukq6RjmC/xPAahFZCRwBrgKu\nzjxARDqAnwAfNMbsK3TBzOCv1HQNhEY4tXXNpG12W/7heCTnvD5p6bx/oeCvHb5qJmQ3irdt21b0\nNUoO/saYuIjcCGwn2YdwpzFmh4jckNxtbgc+DSwG/kmSRddRY8zmUu+tlJXB0OjUtI/bz1gkUPDc\nfNM5p9ldyjFfzl87fNVcK0vO3xhzP7Aua9u3Ml5fD1xfjnspZVey1HNytU+9x8+ojZZ/JB7Lme9P\ns1vxk7PlX8QoY6Vmio7wVRVrIFfL31PDuI2cf6G0j9fpslXrPxoZn1rqqWkfNQ9o8FcVK9cI3zq3\nn9Fo4bRPvtG9aV6np2DLP56IE4xFqPP4J21v9NYxEh7Xyd3UnNLgrypSMBomYRLUuHyTtte6fQRj\nERImYXl+OB7Dk2Px9jQ7tf6j0SB1bt+UclCXw0m9p4Zhm1NLKzUTNPireW08GuT9P/87/vrh/4cx\nxvZ5A+HkbJqZk7pBck7/GpeX8WjI8vxwPILPZRH8Ha6CHb7JMs/anPuafPVFzewZT8T5iwe/zrX3\nf45QLGL7PKXy0eCv5rW7X7ofl8PJ0727eejQ07bPG8yR8kmr89QwWqDiJxKPWrb8vTY6fEdyzOiZ\nVmyn7wMHnmTP0GEC0RA/2v2Q7fOUykeDv5q3jDHc9dIv+Oszr+aDb7iEn+57zPa5uaZ2SKtz+wrW\n+pej2mc0PHWAV1qTr57BIso9f7T7Ia49+VI+cfr7+JcXfm77PKXy0eCv5q2Xjr2K1+nmjUtO4Pz2\nU+k69EzBXH1arkqftHpPTcFRvrYGeRWo9hmNBqjLE/yLmd8nGo/xaPfzXLLqLDYvW89wZJxXh4/Y\nOlepfDT4q3nrwYNPcUHHGQCsbFhKrdtne1qEXJU+aXXuGsYKVPzkm9cnzeN0ES7Q4RuIhqnN6nBO\nK6bc8+WB/bTXt7LIW4dDHLyt/TQeOviUrXOVykeDv5q3nuzZxVnL3jDxfkPLal7of8XWuQNB67TP\naIGWfyQexWdV6ukonPYJxELUuvME/yKmdX7y6E5Ob3t9iYy3HPdGHj/6sq1zlcpHg7+al4wxPNe7\nh40tqye2nbLkBPvB3yrnb2N+n8It/8LBfzwazB/8ffY7fF869ioblpw48f70tnU81bPb1rlK5aPB\nX81L3WN9OB1OltU2T2x745ITeLG/4LyAQO6FXNKSOX8b1T4lBv9ANIw/T9qnmJz/7sHDrG1aMfF+\nZcNSovEo3WP9ts5XKhcN/mpeeq5vHxtbVk+q01/b1F5Uzj97Cce05Chf65Z/yM7EbgU6fMejFmkf\nmzl/Ywx7hw6zOiP4iwint63jyaM7C56vVD4a/NW89FzfHja2nDhpW1tNE8FYhOFw4ZGxVtU+dR5/\nwfl9IgWCv9fpKjjCt1DO307LvycwiMfhmvJ7Ob1tHU/3aupHTZ8GfzWjwvHotGawfLZ376R8PyRb\nvKsal/Hq8GsFz7ca5FVrY36fgmkfOx2+0RA1Fi1/OyN89wweZk1Gqz9tY8sanu/bW/D8bAOhER0h\nrAAN/moG/ebws5z+b3/MWd+/gS8/8X3b5yVMghf697GpZc2Ufcngb13jbozJOZ1zWp3bTyAatryG\nnfn8C3f4hqbMLZTW5KtnKDxWcMqKPUOHWL1oavDf0HIiL/a/OrHiWCHGGP7ud//KWd+/gbf84CM8\nfkSrhaqdBn81I/YPH+HGB7/OHW+/id9e/c/cu+9R7n/1cVvnvjp8hEZvHYv9U9M2doL/eDSEy+HC\n7/Lm3F/r9jFehmqfQnP7WKV9vE43HqerYNXRq8NHOHHR8inbG721tNU22e4D+cmeh3n48DM88f47\n+Hrnx7h++xd5TTuMq1pZgr+IXCIiO0Vkt4jclOeYfxSRPSLyrIhsKsd91fxkjOGTD9/Gx067krOO\nO5lmfyOfP+cj/J/f321rGuPn+vayISvfn3ZC43G8UiDtYzXAC9LB33pit0ii0JTOpaV9wF7e/9Bo\nL+31rTn3bWpZw3N9eyzPh+QI4a8+eQ+fP+cGGr21nNe+iWtPvoybH/lmUZPlqcpScvAXEQdwG3Ax\ncDJwtYiclHXMpcCJxpg1wA3AN0u9r5p5w+Fxdg4cKNhKzrb9wB8YCo9x3cmXTWw7Z/kpLPbV88v9\nfyh4/rO9e3OmfABWNR5XMOdvVeMPUOOy0fKPWQ/y8jhdRBLWHb5W1T5gr+LnsEXw39iymuf6Cpe+\n3rvvUdrrW9mcMWDuxlOv4MDIUX518MmC52cKREPsHjzEiI1OdzW/lWMZx83AHmPMAQARuQfYAmTW\noW0BvgNgjHlcRBpFpM0Y01OG+6sye2X4Nf7+99/hke7nWFbbTE9gkA+94WI+cfpVBRcsj8Zj/P3j\n32Xbm/8Ip8M5sV1E+NDJl3DPrge57IQ3W17jub69XHz8NTn3rWxo48Co9cemUPCvdfsJxEpr+Xuc\nbsIFOk4DsRA1eVJPkOr0DeYf5WuM4fBYHyvytfxbV/Ofe39j+QwAP9nTxYdOvmTSNo/Tza1vvo7P\n/PZOzluxyfL3Csmg/4U/fI8f7nqQ1pom+oJDnN9+Grds/gAdDW0Fn0HNP+UI/suBQxnvD5P8QrA6\npju1TYN/BmMM49EQw+ExBsNjQHLN2SZvPQ3e3PPCZ4ol4uwZPMzzfXt5vn8frw4fwet0s7apnbOX\nn8LZx50yKSDnOv+OF37Kbc/+Bx/d+G7+79v+FzVuH72BQW555Ft84Bd/x12XfCrvTJUA39v5AMfV\nLqGz/dQp+y5bdRafeewO+gJDtNQsyvsMLx/bnzft0+JfRDAWZiwSnLJCVprVAC+wmfaJR/E683/R\neRz26vxr3LmfEVJpH4uW/3BkHGMMjXnWBDi5eRW7Bg9adk73B4d4tncP/3LRLVP2va3jdO568Rf8\ny4u/4CMbt+R9jmAszIfu/3uW+Br57dX/TLO/kbFIkDte/BmX/sdf8WcbL+eGDVtwO/OHk1giziOH\nn+Ox115g71A34XiEJb5G1jS1s3npeja2rs7bR5MpEA3x2vgxIvEoCWPwuTws9tVPzHtUrHgijogg\nyJS1HypdWRZwL7dN37kOYOIvI/1XImT85WTvy/iLSx9nuU9yXNNiX75nmXzt7GMnX1sE3A4XHqdr\nYq740UiA0WiAkfA4w5Fx3A4Xi7x1LPLWAcJYNMBAaASP083xDctY1biM5XVLWOJvxOVwEYiG2D98\nhJ2DB9lx7ADL6prZsORENrScyAUdpxOKRdgxcIDPP/5v9AQGuHJtJ+9bd8GUTsRnendz8yPfotFT\ny8/f/SVWNiyd2Nda08S3L/rffOrR2/nAL/4P37vsMzkD73B4nK8/9UN+8D+25vyHVOv2c9Hxm/mv\nvY9w/YZ3TtkPsGvwIMvqmvN+wYgI7XWtHBrtZX3zypzHWA3wSj5H4eAfikfxWAQzr42VvAIF0j6F\nRvmmUz75glKN28eqhmXsOLafTa2502Q/f+X3vK3jdPzu3IH11jdfx+X3fYor157HEv/UL+RIPMqf\nPvBl2mqa+MfzPz7ReKjz+Plfp72HK1afy82PfJN79z3Kl8/98ynPsW+om3t2PshP9nRxXN0SLug4\ng/etexs+p4f+0DAv9b/KZx+/m50DB1mzaAWbWtewvG4JDd5awrEo/cEh9o8c5eBID4fH+ghEQyyt\nXZz6ohBC8TCDoVHGokEaPLWpXzUTn59YIk40ESMcjxKKRwjFIoRi4YnXiVSfhyH5f0FwSPKLwJH6\nQkh/OTiy/p/+e3GIA4E5+RIppc+mHMG/G+jIeL8itS37mPYCx0x4+65W0r+nt7z1bN5y7jkTfznw\n+m/YvL7h9X1M3mdy7cs+v8A+sj4g1tee/D772tFEjEg8lqoUMdR7aqhPfWAbvXU5W3DGGPqDw+wf\nOcIrw0c4On6MAyM9RBMxalw+3rBkFe9ecy4nN6/K+RPCZSe8mU+ecRW7Bw/xw50PcsV9f0N7fRub\nWlbjdDh5pnc33WP93HTmNbxn7fk5P7wOcfC5c/6Um3/zTT58/2f57qWfntKZ+bWnfshFK8/kDc3H\nTzk/bcuJ5/D1p36UN/g/3bOb01rX5j0foL2hjUOjPZbBP1elUFqt20/AVsu/QM6/YLVPuGDaxzr4\n50/5pKXz/vmC/0/3PZr3zxpgddMKrlh9Hl964gd86dw/m7QvnojzFw99A5c4+Ebnx3L+1NjR0Mb3\nLvsM/7HnYa69/3Oc0HgcG1pOJJqI8UzvHrrH+njP2vP54Tu2saapfcr571l7PpD86eLF/ld4rm8v\nR8cHODjai8/pZrG/kUtXnUVHfRvt9a0s8Tfm/HxG4zGGwmOMRMYZiYxPLNbjdrhwOVz4nG58Lg8+\nlxef05N87fRM+mnFGEPCJDAw6f8YSJDAGLL2G5L/GRLGYDAYk/zXPtPh/7HfPMpvH0mubyECz03j\nGuUI/k8Aq0VkJXAEuAq4OuuY+4CPAj8UkbOAIat8/5f//otleKzKISK01CyipWYRZy5dP+3rrG1q\n59NvvpabN3+A3x95mR0D+0mYBOef/j7ectwbC+Z9HeLgC+d+hE903ca193+Ouy/9m4kf1X9z+Fl+\ntu8xtl/5NctrvHX5Rj7263/g8GhvzsD2VM8uTm9bZ3mNjvpWDlrk/QdCo7xxyZK8+71ON9FEjFgi\njitPGqxQ2qfQAu7pfVZ/pk2+evYN5W0DcWi0lxV1LXn3A2xsXc2zvXuAS6bsOzo+wI6BA3S2n2Z5\njU+c8T7O++GNvHft+ZyxNFmrEUvE+eTDtzEUHuPuS/7GMqUjIvzPtZ2848SzeeTwc+wePITL4eSd\nJ5zNaW1rC36uAPwuL2cuXT/tz7fb6Zr4NzJdIoJT0p+H/OnR+eCKS9/FFZe+a+L9Vz73paKvUXLw\nN8bEReRGYDvJ6qE7jTE7ROSG5G5zuzHmFyJymYjsBcaB60q9r5o+t9PFW1ds4K0rNhR9rkMcfPW8\nj/LxX/8j197/OT591ofZN9TNpx+7g3++8K9o9jcWvPelx7+Jn+57jD/b9O4p+5/q2cUNG/LnnwHa\n69s4ONKbd/9ggQ5fEZlI/TTm6UsJF0j7FKrzD8TClikfSE/rnL/lfyjPF2SmTS1ruPul/8657+ev\n/JYLV55h+RMMwCJvHV8+96P8yfYv8rlz/pTl9S18+YnvkzCGuy6+pWAnf5rX6ebClWdw4cozbB2v\n5lZZcv7GmPuBdVnbvpX1/sZy3EvNPafDyTfO/xjfeOpHfPTBr9FS08SdF99su9X2zhPP4fOPf3dK\n8B8IjtAXHJo0g2UuHfWtPH7kpbz7C1X7ANS4/YxHg3mDf6HpHbxOt2Wp53g0VLADs1Dap3usjzdl\nlGfmctLiDvaPHM05puC+fY/xsdOutDw/7aLjz6TG/Zd89cl7GA6Pcfnqc/nIxi22Wu1qYZqXHb5q\n/nM5nPzVmVfzV2dmZ/gKe8txb+S18X5eHT7CqsZlE9sfP/oyp7eus6xIgmSe+dBo/pb/saD1IC9I\ndvpa5f1D8Yh1tU+BnH/AYi7/tMW+Bstqn2TL3zrt43G6Obl5FU/17J70k9zh0V72DXfz1uX2f7o7\nZ/kGzinieLWw6fQOata5HE4uW/XmKQuydx16hvPaCw/+XpHK+eerdOgPDdOSo3IlU63Lx7hFrX8k\nHitpbp9kZ2+BtI+Nln97nXXaB+C8FZt4+PAzk7b9ZM/DvPOEs7XlrvLS4K/mxLtOPJv7MoJ/wiR4\n6NDTOccHZFvkrcMpjpyBMxqPMRoJpEpl8ytU7hmJR/E4LEo9HdY5/0Kje8F6eoeR8DiReCzv5HSZ\nOttP5VcHnny9Us0Yfry7iyvXdhY8V1UvDf5qTmxeup7hyFiqUgV+99pLNHrqWJNjBstc2uvbclb8\nHAuN0OStL5g6qikwuVs4HsFr0dFZqOU/Hg1azusDyS+gaCKWc4rl5MjeFls146e2riEcj04scfnY\nay/gEClYMquqmwZ/NSecDid/tvFyvvTE9zHG8M3n7+Wa9RfaHiCTLPecmvfvDw6zpEDFEaSndS6Q\n9nEUSPuO1KqMAAAR60lEQVRYdPgGooXTPiJCs6+BYzkWcrdT45/mEAdXn3Qh33zuv4gn4nzlyR/w\n55uuqLoRq6o4GvzVnPnA+osYiQR413/dTPdoLx9Yf5Htc9sb2jicM/gP5Rypms0q7WOMsVHqWaDD\n12I650wtNYvoCwxO2X54rJf2AjX+ma4/5Z0827eHLfd+CrfDxZVrzrN9rqpOWu2j5ozH6eZ7l32G\nBw8+aWtysUwd9W3sGjw4ZXt/cJgWGy3/Grc/71z60UQMp8NhmTpKT+9gjMnZwh63Ue0DyakzenME\nf6upnHPxu73857s+x2OvvcDFx28umPZSSlv+ak41emu5Ys15BQeHZWuvb+VQjoFe/cEhW9eyavlH\n4rGJuZfycYgDp8NBNE/qJxAN4y+Q9oHkusQ9uVr+NgZ4TblW7WKuWHMetRaTySmVpsFfLUgdDbk7\nfK1mDM1kVecfjkdsjWpNruObO/jbqfYBaKlpoi8wNGV7sS1/pYqlwV8tSCvqWuge60tOvJWhPzTM\nEp+Nlr/Fgi6hApO6pVnl/YM2c/6tNU30BAambNfgr2aaBn+1INW4fdR7aqakTPoDwyyx0fKvcecf\n5BWKhS1H96Z5nfnn9B8vsIRjWluOnP9oJEA4Hi04RYVSpdDgrxas9vo2Do1MTv30BYdslXpaTesc\njkfxuey0/PPX+o9HQwVLPSF3h+/h0T7LefyVKgcN/mrBylXr/9p4P8vr8k/nnGbV4RuOR/DZaPlb\nBf9ALERNngVUMuXq8LUzm6dSpdLgrxas9vrkoi5pgWiIYDRMs52cv9ufP/jHrOfyT/M43ITzdPgm\nV/EqXHXTUtNEf3B4Ut9FsTX+Sk2HBn+1YHU0tE6a3fO1sX6W1TXbSpfUWkzvEIpH8NpI+3idLsLx\n3Iu426328Trd1Lp9DGSM8j0w0qOdvWrGafBXC1Zyfp/Xg3/3WD/LbbaYC6V9bLX8C6V9bCxIDskB\na5mL07wy1M2JNuc4Umq6Sgr+ItIkIttFZJeI/FJEpvy8LSIrROQhEXlJRF4QkY+Vck+l0jrqWyd1\n+HaP9dnK90OBln8sWkTOP3+df43NwVYrG5ZyYOToxPu9Q92sblpu61ylpqvUlv/NwK+MMeuAh4Bb\nchwTAz5hjDkZeDPwURE5qcT7KsVxdUvoCQwSTQXgYtIltS5/3lLPsO20T/5pnYM2lnFMW9mwlP2p\n4B+MhekNDNJR32brXKWmq9TgvwW4O/X6buDy7AOMMUeNMc+mXo8BOwBt1qiSeZxuWmua6B7rA2D3\n4EHWNnXYOtfn8hCJx4gn4lP2FVrFKy25lKNVqae9tM/xja+3/F8dPkJ7fWveheWVKpdSg3+rMaYH\nkkEesGx2icjxwCbg8RLvqxQAb2g+fmIe+92Dh1i32F7wFxFqXF4CsfCUfeFYFJ+NEb5ep4dwjrn4\n0/0AdieqW9WwjL1D3QA837eXU5acYOs8pUpRcFZPEXkAyPwZVAAD/G2Ow3Ovq5e8Th3wY+DjqZ8A\n8tq6devE687OTjo7Ows9pqpSp7au4ZnePVzQcTpHAwMc37DU9rnpTt96T82k7Xbn9vG5PIRyVPvY\nrfRJO6XlBHYNHCQcj/J0z25Oa1tn+1xVnbq6uujq6irpGgWDvzHm7fn2iUiPiLQZY3pEZCmQc1Vt\nEXGRDPzfNcbcW+iemcFfKSunta7jK0/+gBf7X2Htovai0iXJWv+pnb6hmL20j8/pybkKVyAWtjW1\nQ+ZzrGpcxov9r/BU7y6uWZ/3n5xSwNRG8bZt24q+Rqlpn/uAa1OvPwzkC+z/ArxsjPmHEu+n1CRn\nLF3H7sGD/NuOB3jrig1FnZuv3DNsc2K3fC3/QDRoO9+ftnnper778i/pCwzxRk37qFlQavD/IvB2\nEdkFXAB8AUBElonIz1KvzwbeD7xNRJ4RkadF5JIS76sUAH6Xlw+94VJ+sqeLa04qrsWcr9zTdtrH\n6SEYndpnMG5zdG+m96+/iH/f/Wuu3/Au7exVs6KklbyMMQPAhTm2HwHekXr9GKCfZjVjbt78fj55\nxvuKWgkMkqt55Wr5J6d0tpfzHwpP7b6yO6NnpvXNK3n1T35U9O9BqenSEb5qwRORaQXN5Jz+OYJ/\nLGIr7eN3eXPm/MejIWptzOiZTQO/mk0a/FXVqnX7COQY6GV3eod8OX+7C7koNZc0+KuqlT/nH8Vv\nM+efr+VfbNpHqdmmwV9VrZq81T52W/5eQjkGidldyEWpuaTBX1WtujyrednN+ecd5KVpH7UAaPBX\nVavW7WcsT9rHzsRueQd5RYMa/NW8p8FfVS2rtI/PWXiQVv5BXuGi6/yVmm0a/FXVyjfCNxSzN6Wz\nz2XR4VvkCF+lZpsGf1W1al2+3Dl/u4O8nB6CuYJ/LKjVPmre0+Cvqlat28d4LPf0Dv4Sqn0CRc7q\nqdRc0OCvqla+6R3CMXsdvn6n1ZTOmvNX85sGf1W1at1T0z7GGNsreeWb3iFQxBKOSs0VDf6qatXl\naPnHEnEEbM2s6XUl1/A1ZvIaRuPRIH7t8FXznAZ/VbVq3b4pdf5hm529AA5x4HG6pqR+Apr2UQuA\nBn9VtXwuD+F4dNIi7oFYiBq3/VZ7roFexS7jqNRc0OCvqpZDHPhdk8s1A0XOy5M90MsYQzAW1jp/\nNe+VFPxFpElEtovILhH5pYg0WhzrSK3idV8p91SqnLLX8U2uv1tEyz9roFcoHsHtcOHU1bjUPFdq\ny/9m4FfGmHXAQ8AtFsd+HHi5xPspVVY1Lh/jGXP6B6Ih/MW0/LPSPgGdzlktEKUG/y3A3anXdwOX\n5zpIRFYAlwF3lHg/pcoqu9wzUGTKxufyEoq/PtBL8/1qoSg1+LcaY3oAjDFHgdY8x30d+GvA5Nmv\n1JzInt8nGCuu5Z49xYMGf7VQFFzAXUQeANoyN5EM4n+b4/ApwV1E/gfQY4x5VkQ6U+db2rp168Tr\nzs5OOjs7C52i1LTUZef8o+GiavSTOf/XW/6BmC7komZeV1cXXV1dJV2jYPA3xrw93z4R6RGRNmNM\nj4gsBXpzHHY28C4RuQzwA/Ui8h1jzIfyXTcz+Cs1k2rcPsYmtfyLTftMzvnrEo5qNmQ3irdt21b0\nNUpN+9wHXJt6/WHg3uwDjDGfMsZ0GGNOAK4CHrIK/ErNpmTOP7Pap7jgXePyEYhl5vx1IRe1MJQa\n/L8IvF1EdgEXAF8AEJFlIvKzUh9OqZlWmzXFQyBaXMs/u89gJBKg3lNT1mdUaiYUTPtYMcYMABfm\n2H4EeEeO7Q8DD5dyT6XKqcY1OXgHYiH8RbTcs/sMRiPjNHhqy/qMSs0EHeGrqlp2y73Yln+dZ/I6\nwKPa8lcLhAZ/VdVq3T4CsaxSzyKqdWrdfsYimcE/SIMGf7UAaPBXVS3Z8p/+9A650j7a8lcLgQZ/\nVdWyV/MqdmK3WvfktM9IJECdW4O/mv80+KuqVpfV8i+2Tj97eojRSEDTPmpB0OCvqlp2qedoJECj\n1361Tl2Olr+mfdRCoMFfVbWpdfrFlWpOrfYZp6GILw+l5ooGf1XVarKC/3BknMYign/2egBa6qkW\nCg3+qqrVul7P+Rtjig7edTlKPeu1w1ctABr8VVVr8NYyHB7DGEMgFsLjcON22h/4XpfRZxCOR4km\nYjq3j1oQSpreQamFzu/y4nK4GI+Gkvn+IvP1PpeHmIkTiUcZCI2w2FePSMFZy5Wacxr8VdVr9jdw\nLDRMMBamsch8vYjQ5K1nIDTKseAIi315l7FWal7R4K+qXrOvgWPBEWKJGPXTmJSt2d/AQGiEgdAw\ni331M/CESpWfBn9V9Rb7ksHbYKZVprk49eUxEBql2a8tf7UwaPBXVW+xv4FjoRFc4iiqzDOt2dfI\nQGiYY0Ft+auFQ4O/qnrNvkYGgsO4nS4WeeuKPz/15XEsNEKz5vzVAlFSqaeINInIdhHZJSK/FJGc\nn3wRaRSRfxeRHSLykoi8qZT7KlVOi33JDtvewCCtNU3TOL+BgeBIstrH3zADT6hU+ZVa538z8Ctj\nzDrgIeCWPMf9A/ALY8x6YCOwo8T7KlU2LTVNHA0M0BsYom0awb/Zl2z59wYGafEvmoEnVKr8Sg3+\nW4C7U6/vBi7PPkBEGoC3GmPuAjDGxIwxIyXeV6myaa9vpXu0j57AAC01xQfvZn8j/cEhDo320lHf\nOgNPqFT5lZrzbzXG9AAYY46KSK5P/iqgX0TuItnqfxL4uDEmmONYpWZdR30rB0Z7cDucHN+4rOjz\nVzYs5dXhIxwe62OFBn+1QBQM/iLyANCWuQkwwN/mONzkucdpwEeNMU+KyDdIpotuzXfPrVu3Trzu\n7Oyks7Oz0GMqNW3H1S2hPzhELBGno76t8AlZTmg8jh0DB1jib6RJq33ULOjq6qKrq6uka4gxueK1\nzZNFdgCdxpgeEVkK/DqV1888pg34nTHmhNT7c4CbjDHvzHNNU8ozKTUdF//kk4xFgjx29T9N6/yT\n//WDnLLkRO55x9byPphSNogIxpii5hUpNfh/ERgwxnxRRG4CmowxN+c47mHgemPMbhG5FagxxtyU\n55oa/NWsC0bDhOKRabfcB4Ij1Lh9+FyeMj+ZUoXNRfBfDPwIaAcOAO81xgyJyDLg28aYd6SO2wjc\nAbiBV4DrjDHDea6pwV8ppYow68F/JmjwV0qp4kwn+Ot8/kopVYU0+CulVBXS4K+UUlVIg79SSlUh\nDf5KKVWFNPgrpVQV0uCvlFJVSIO/UkpVIQ3+SilVhTT4K6VUFdLgr5RSVUiDv1JKVSEN/kopVYU0\n+CulVBXS4K+UUlWopOAvIk0isl1EdonIL0WkMc9xfykiL4rI8yLyPRHR5Y6UUmoOldryvxn4lTFm\nHfAQcEv2ASJyHPAXwGnGmA0kF3S/qsT7KhtKXeBZTaZ/nuWlf55zq9TgvwW4O/X6buDyPMc5gVoR\ncQE1wGsl3lfZoP+4ykv/PMtL/zznVqnBv9UY0wNgjDkKtGYfYIx5DfgqcBDoBoaMMb8q8b5KKaVK\n4Cp0gIg8ALRlbgIM8Lc5Dp+y+K6ILCL5E8JKYBj4sYhcY4z5/rSeWCmlVMlKWsBdRHYAncaYHhFZ\nCvzaGLM+65grgYuNMden3n8QeJMx5sY819TV25VSqkjFLuBesOVfwH3AtcAXgQ8D9+Y45iBwloj4\ngDBwAfBEvgsW+xtQSilVvFJb/ouBHwHtwAHgvcaYIRFZBnzbGPOO1HG3kqzwiQLPAH9ijImW+vBK\nKaWmp6Tgr5RSamGaFyN8ReTK1CCwuIiclrXvFhHZIyI7ROSiuXrGhUpEbhWRwyLydOrXJXP9TAuN\niFwiIjtFZLeI3DTXz7PQich+EXlORJ4RkT/M9fMsNCJyp4j0iMjzGdtsDbjNNC+CP/AC8G7g4cyN\nIrIeeC+wHrgU+CcR0T6B4n3NGHNa6tf9c/0wC4mIOIDbgIuBk4GrReSkuX2qBS9BslDkVGPM5rl+\nmAXoLpKfx0wFB9xmmxfB3xizyxizh2QZaaYtwD3GmJgxZj+wB9APS/H0C3P6NgN7jDEHUv1U95D8\nXKrpE+ZJ7FmIjDGPAoNZm+0OuJ0w3/8ClgOHMt53p7ap4twoIs+KyB12fhxUk2R/Bg+jn8FSGeAB\nEXlCRK6f64epEAUH3GYrtdTTNovBYn9jjPnpbD1HJbL6swX+Cfg7Y4wRkc8CXwP+ePafUqkJZxtj\njohIC8kvgR2p1qwqn4KVPLMW/I0xb5/Gad0ky0jTVqS2qQxF/Nl+G9Av2uJ0Ax0Z7/UzWCJjzJHU\n//tE5D9JptY0+JemR0TaMgbc9hY6YT6mfTLz0/cBV4mIR0RWAasBrQ4oQuqDkHYF8OJcPcsC9QSw\nWkRWpqYiv4rk51JNg4jUiEhd6nUtcBH6mZwOYWqsvDb1Ot+A20lmreVvRUQuB/4vsAT4mYg8a4y5\n1Bjzsoj8CHiZ5ACxPzc6MKFYXxKRTSQrLPYDN8zt4ywsxpi4iNwIbCfZWLrTGLNjjh9rIWsD/jM1\njYsL+J4xZvscP9OCIiLfBzqBZhE5CNwKfAH4dxH5I1IDbgteR2OpUkpVn/mY9lFKKTXDNPgrpVQV\n0uCvlFJVSIO/UkpVIQ3+SilVhTT4K6VUFdLgr5RSVUiDv1JKVaH/D/ypppTsDNN1AAAAAElFTkSu\nQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x0 = -10.\n", "x1 = 10.\n", "x = np.linspace(x0, x1, 1000)\n", "def F(x):\n", " return np.exp(-abs(x)) * np.cos(2 * np.pi * x)\n", "p = plt.plot(x, F(x), color='#1a9850')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we generate the sample data `(X, y)` using\n", "`scipy.ndimage.convolve`. This performs the convolution\n", "\n", "$$ p\\left[ s \\right] = \\sum_r F\\left[r\\right] X\\left[r - s\\right] $$\n", "\n", "for each sample." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import scipy.ndimage\n", "\n", "\n", "n_space = 101\n", "n_sample = 50\n", "np.random.seed(201)\n", "x = np.linspace(x0, x1, n_space)\n", "X = np.random.random((n_sample, n_space))\n", "y = np.array([scipy.ndimage.convolve(xx, F(x), mode='wrap') for xx in X])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this problem, a basis is unnecessary, as no discretization is\n", "required in order to reproduce the convolution with the MKS localization. Using\n", "the `ContinuousIndicatorBasis` with `n_states=2` is the equivalent of a\n", "non-discretized convolution in space." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pymks import MKSLocalizationModel\n", "from pymks import PrimitiveBasis\n", "\n", "p_basis = PrimitiveBasis(n_states=2, domain=[0, 1])\n", "model = MKSLocalizationModel(basis=p_basis)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the model using the data generated by $F$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.fit(X, y)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To check for internal consistency, we can compare the predicted\n", "output with the original for a few values" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-0.41059557 0.20004566 0.61200171 0.5878077 ]\n", "[-0.41059557 0.20004566 0.61200171 0.5878077 ]\n" ] } ], "source": [ "y_pred = model.predict(X)\n", "print y[0, :4]\n", "print y_pred[0, :4]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a slight linear manipulation of the coefficients, they agree perfectly with the shape of the filter, $F$. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD7CAYAAACCEpQdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecnVWd+PHPub1MrymTTHonBFIgFEnoIF1EQKUposLq\nT1yluRvcXVzZVRcVCyJKkY6ggdAhhSSUBFJJb5MpyfR6ezm/P2YymXJr5k7ulO/79ZoX9ynneU4u\nd773zHnO+R6ltUYIIcTwYkh3BYQQQhx/EvyFEGIYkuAvhBDDkAR/IYQYhiT4CyHEMCTBXwghhiFT\nuivQk1JKxp4KIUSStNYqmfMHZMtfay0/KfhZsmRJ2uswlH7k/ZT3c6D+HIsBGfyFEEL0Lwn+Qggx\nDEnwH8IWLVqU7ioMKfJ+ppa8n+mljrW/qL8opfRAq5MQQgxkSin0UHjgK4QQfTFu3DiUUkPuZ9y4\ncSl7j6TlL4QYcjpawumuRspF+3dJy18IIURCJPgLIcQwJMFfCCGGIQn+QggxDEnwF0KIYWjAJXYT\nQojhYv369SxZsoSGhgZuvfVWwuEwZWVlPProoxw+fLhf7y3BXwgh0mTevHnY7XZuuOEGbrnlls79\nTqez3+8t3T5CCJFGK1euZOHChQA0NDQAMGbMmH6/r0zyEkIMOYlM8hr9yJV9vk/lba/0qfymTZs4\n44wzaG5uxmAw8Mgjj3DbbbdFPT+Vk7yk20cIMSz1NXCnwooVKygtLeWpp57ivffe49JLLz1u95bg\nL4QQabJ8+XJuueUWbrzxRmbMmJHS3D3xSLePEGLIGQy5fcLhMAUFBXz00UdMmTIloTIDLrePUuox\npVS1UmpzjHN+o5TarZTaqJSak4r7CiHEYLRp0ybuuecevF4vy5cv7/dhnZGkpOWvlDoDaAOe1FrP\njnD8IuAOrfUXlVKnAL/WWp8a5VrS8hdC9MlgaPkfiwHX8tdarwYaY5xyOfBkx7kfA9lKqeJU3FsI\nIUTyjtc4/9FAeZftyo59QgwK/oAff8Cf7moIkTIDcrTP/fff3/l60aJFstanSLtJc2aQnZfDlg/W\np7sqQrBixQpWrFjRp2ukbLSPUqoUeDVKn/8fgeVa6+c7tncAZ2mtqyOcK33+YsBRSmFx2PC5POmu\nikiA9PnHl8puH9XxE8lS4AYApdSpQFOkwC/EQJVbUkR2cX66qyFEyqSk20cp9QywCMhXSh0ElgAW\nQGut/6S1fl0pdbFSag/gAm5OxX2FOF5OuPIMWsvr0l0NIVImJcFfa319AufckYp7CZEO4xbPxm6w\npLsaQqTMgHzgK8RAE7YawGxOdzWESBkJ/kIkwBPyYzQY010NIVJGgr8QCfAGfRiVLH8hhg75NAuR\ngMr1u2ltaU13NYRIGQn+QiRgx2Pv09LQlO5qCJEyktJZiAQopbCNyMZzSL4ABgOZ5BWftPyFSJC/\nyZ3uKogUuv/++1FK9frpml4m1vnRzkvUoUOHeOCBB1i2bBl33XUXZWVltLW1UV19fOa/SstfiAQo\npTDYzIQ8ktxtMBjoLX+3283ixYt5/fXXyc/PZ926dTz44IN89atf5ZJLLsEcZVixtPyFOI5CoRDA\ngA4mYnB5/vnnmTt3Lvn57SlDioqK2Lx5M1rrqIE/1ST4CxGHy+/BfvJodCAkXwAiJfx+P5MnT+7c\ndrlcGI1GrrrqquNWBwn+QsQRVGHG/8tiANw+yeop+u66666jvr6eN954g6VLl1JVVcVJJ53E448/\njsdzfD5j0ucvRBxVbXVc+spdNFbXs+GOv5Fty0h3lUQcA73P/1hJn78Qx5E35MdmspBZlIMvHEh3\ndYRICUnvIEQc3qAfu8lKSIfxBmW0jxgaJPgLEYcn6MNmshAMh/AEfemujhApId0+QsRRWVVF6/ZD\n2EwWvCFp+YuhQYK/EHFs2biZPa983B78pdtHDBES/IWIo83dhsVqYeP/LePDD9akuzpCpIT0+QsR\nR5vLhcVmxeNy09gkid0Gg9LSUpRKauTjoFBaWpqya0nwFyIOl9uN1WYjFArT5naluzoiAQcOHEh3\nFQY8Cf5CxOH2uLDZbYSCIdpcbemujhApIX3+QsSRNbKA0pmTsNqtuNyS1lkMDdLyFyKOyaefgFEZ\nqdxZhtsj3T5iaJDgL0Qc3qCffHs2F9z6JaxmS7qrI0RKSLePEHEcye2Tm5OLtsivjBga5JMsRBye\noA+7yYrdZJX0DmLIkG4fIeLwBv3YjBaCppDM8BVDhrT8hYjjwIYduOqbsZkseCS3jxgiJPgLEceG\nZ9+nas9BbEbJ7SOGDgn+QsQR8PrJdGSw45PNvPvfT6e7OkKkhAR/IeII+AJkOTOwGE246pvTXR0h\nUkIe+AoRR9AfIDMjE1/AT8AnyziKoUGCvxBxhHwBcpxZ+P0Bgj7p8xdDgwR/IeLIOGEUhXkF+AJ+\ngtLyF0OEBH8h4ij42lxGFY+k1dVGyB9Md3WESAl54CtEHN5ge3qH8WNKGb/konRXR4iUkOAvRAzB\ncIiQDmMxmHBa7QSd8isjhgb5JAsRw5HUDkopbCYLvmAArXW6qyVEn0nwFyIGb8iHzdSextmgDFiM\nJryS4kEMARL8hYihrrEB76ZDnduS4kEMFRL8hYjhwMEyKp9f17ltM1mk5S+GBAn+QsTQ7GrBZDF3\nbu/+j2Xs3LMrjTUSIjUk+AsRQ6urDZP1aPDXviBNLZLfRwx+EvyFiKHF1YbZenTdXpPFTEtbWxpr\nJERqSPAXIobWtjbMVmvnttlmocXVmsYaCZEaEvyFiMGRl8nYuVM6t81WK60S/MUQILl9hIhh7KxJ\nzM87p3PbYrPS6pJuHzH4SfAXIgZv0I/ddLTb57wfX8/sifPSWCMhUkO6fYSIwRvyYTMefeCbnZNN\nyJjGCgmRIhL8hYjB05HR8wi7yYon6EtjjYRIDQn+QsTgDfqwden2kfQOYqiQ4C9EDHs376C5orZz\n22aS4C+GBgn+QsSw/tVV7Fu3rXNbcvuIoSIlwV8pdaFSaodSapdS6q4Ix89SSjUppT7r+PlJKu4r\nRH/zer04HI7O7fVLV/L6w8+msUZCpEafh3oqpQzAw8A5QBWwTin1T631jh6nrtJaX9bX+wlxPPm9\nPjIczs5tIwZaG1vSWCMhUiMVLf8FwG6tdZnWOgA8B1we4TyVgnsJcVy1B/+Mzu0MpxOfx5vGGgmR\nGqkI/qOB8i7bFR37elqolNqolFqmlJqRgvsK0e/8Xh8ZzqMtf6fDic8nQz3F4He8Zvh+CozVWruV\nUhcB/wCmRDv5/vvv73y9aNEiFi1a1N/1EyKi3FljGDtmTOd2piMDv1eCv0ivFStWsGLFij5dIxXB\nvxIY22W7pGNfJ611W5fXbyilfq+UytNaN0S6YNfgL0Q6jb1yHtOmTevczsrIIOCV0T4ivXo2in/6\n058mfY1UdPusAyYppUqVUhbgWmBp1xOUUsVdXi8AVLTAL8RA4g11z+0z75T5nHT3FWmskRCp0eeW\nv9Y6pJS6A3ib9i+Tx7TW25VSt7Uf1n8CrlZKfQcIAB7gK329rxDHg7dHeodMeyYha4wCQgwSKenz\n11q/CUztse+RLq9/B/wuFfcS4njyBH09cvtY8MgMXzEEyAxfIWLwBv3YjF1y+5isMsNXDAkS/IWI\nIhwOU796N1bj0QXc23P7yGgfMfhJ8BciCpfXTeNTn2E2Hu0dlayeYqhQWut016EbpZQeaHUSw1NZ\ndSUTxo8n5D4a7MPhMOZsO56GVixmS4zSQhw/Sim01kllUZCWvxBRNLc1Y7B0HxNhMBjQ/hBNLsnv\nIwY3Cf5CRNHU2oLR0ntAnMFiorGlOQ01EiJ1JPgLEUWzqxWj1dxrv8FiorlNWv5icJPgL0QURpuZ\nkaf0TkFltJhokuAvBjkJ/kJEkTeykNlfP7vXfpPVLC1/MehJ8BciCm+oe2qHI07/r+spmTI+DTUS\nInUk+AsRRc+8PkdkZmcRJJyGGgmROhL8xZB3sKWaRze/mnQ5T9CHzdg7+NtNVjzHMMv34Q1/p8bd\nmHQ5IfqDBH8x5G2o2c3zO99LulzPdM5H2EyWY8rv8+SWN9haty/pckL0Bwn+YsircTdS2VabdLm9\nO3bTsKuy136b0Zp0fh9/wM/6bz5OVXPy9RCiP0jwF0NejbuRihfWU1l3OKlyn674kN0rN/bafywt\n/017t6McZuoCrUmVE6K/SPAXQ96zD/yRtjd38emOzUmVc3s82Oz2Xvs/fGwZb/ztlaSutWnXVox5\ndmpcsoCdGBgk+Ishr7b8MCjYsntbUuW8Hg+OCMHfiKK1Jblx/tv37kKF4cPXVyRVToj+IsFfDHlt\n9c3kThnNrn17kirniRL87Q4HbrcnqWvtPbAPp83BxpdXJVVOiP4iwV8Med7GNqbOmUHZwbKkyvk8\nXhyOjF77HXY7Xk9ywf9wTTVT58zAXS8zg8XAIMFfDGmNrc2EfUHmnTyPqoqqpMrmTy1h/JSJvfY7\nHA68Xm9S15p6w1n84L4f4292Ew7LBDGRfhL8xZC27cAuLNkOzj5rEflnTEqq7MQLTubkU+f12p/h\ncOL1JBf8K9tqmTFyIgazkX2HDiZVVoj+0DtZuRBDiCnfyaKHbmHejDkE92QlVTZaeodzLr8Q/wm5\nSV2rsq2W0RkF2PIy2HZgF5NGj0uqvBCpJi1/MaTVeZspKRzJCEcete4mguFQwmU9IX/E9A6ZzgxC\nvdP8R9XicxEMh8mxZjD9itPwm6XbR6SfBH8xpNW4Gyly5GI2miiwZ1OdxDh7b9AXMb1Dsrl9qlx1\njM4oQCnFwivPQef0vqYQx5sEfzGkVbsbKHK0d9GMzihMKs1DtG4fm9GCN5j4DN89NeUUm7IBKHLk\nSnI3MSBI8BdDWo27keLO4F9ARRLBv3z1NoKe3kHeZrImld7hxaefY8fjKwAocuZS7ZLgL9JPgr8Y\n0qpdR1v+TevLePXlpQmX3f/UGvyu3qN6bCZLUondyg6WUTKmBIBiRx7VbknxINJPgr8Y0l793u9p\nPlADgLE1yMYP1yVcNuwPkpPRe4RQXWU16/71mYSvU1VRxYRxEwAolm4fMUBI8BdDmru2halj2gPv\n5AmTqKmqTrhsOBAi25nZa3+G3UGgNfEZvg2Hapg+cTIA9oCRz55Nfm0BIVJNgr8YstxeDyGXjykd\nwf+EydNpPlyfUNlwOIwOhCK2/LMzsgn7Ex8y2lrTxOzJMwEozsin/JVPZZavSDsJ/mLI2la2G1OW\nDYu5fcTOyVNn46ltTijwtnndYFCYTb0H9OdmZhP2BxOqgz8YIGzQzOkI/iPzi0ApDjXUJPEvESL1\nJPiLIWvngd048o5224wpHAlAee2huGXdfi85X4icDiLD5oCwxh+IP+Kn1tPE7F9+hQy7s3OfNdfJ\n5/t2xi0rRH+S4C+GrL0VZWTmH03DYDAYmPWDi6nzx8+sabAYKb3p9MjHDAaU2UizK/6qXJVttYzO\nLOy2LyM/i51le+OWFaI/SW4fMWSNXziTa8fd3m3f7DPm0RBui1s22gSvI6b835cwWOL/+lS21TE6\no3vwzy7MZ3/FgbhlhehP0vIXQ1aNu5ERWQXd9o3OLKSyrS5uWW+UvD5HOLMy8YYDca9zJKFbV6de\nthjnmPy4ZYXoTxL8xZBV7WqgyNk9++bojEKqEpjl6wv6scZo+duMZnwJzPKN1PI/bfGZmEdnxy0r\nRH+S4C+GrOouqR2OGJ1RQEVr/OAfr+XfPss3fvDfvnsnhebucwUkv48YCCT4iyHrSEbPrtqTu8Xv\n9ik7eJDGTdEXXUk0udu79z+J71D3B8NFjlwOS4oHkWYS/MWQdbi5jmJHXrd9Dr+RVT97Lm7ZbVs+\nZ9+y9VGPW42WhLp9PHUtnDx1Vrd9I5x50vIXaSfBXwxJ/oCfjd96kmyjo9v+SSNKafqsDLc3dnoG\nl9uFxRq922fdgy+zevmqmNcor6lChzWlxSXd9ku3jxgIJPiLAa2prYUxp06nobUpqXK7Kw5gdFhw\n2rsHf7vVhinDyp7K/THLuzxuLDZb1OMGFC2tsecLbC/bgzXHicHQ/dcs2+Lk8DPrqW9O7gugvKYK\nR15WQpPLhIhHgr8Y0B5/7XkqPt7Bw8//Jaly2w/sxp7XOykbgMlmpbY5dp+7y+3Caou+4pbZasXt\nif3XQ0NzI2ZH72sopfBuOsTW/Ttilu/pqfdfxtPYyovvv5pUOSEikeAvBrR/vLEUc66D5196Ialy\nu8r2kJmfE/GYyWamsTn2XxIejydm8LdarbS5XTGv4Q37yRlXFPGYIy+THQf2xCzf0/48F4Xnz+C5\n115OqpwQkUjwFwPavkMH+eHP7sN/9ggCocSSqQEcqDhIblHkiVQWh42GOF0u+WNHMHHOtKjHrXYr\nbrc75jXGTJvAF/71yxGPZRfmsfdg7K6nrjwBHx9UbuaWq77KxyvXJFxOiGgk+IsBq8Hbgv3amfz0\n1h8xZdJkPjq0LeGyNQ11FI8cEfHYgm9fwsiJY2OWH79gBvMv+ELU41arDU+8h8YBLw5z5OcGBUWF\nHKwsj1m+qxUVGzixcBLfv/qb1O2qTPp5gRA9SW4fMWCtrdrKgpEzsBjNXDT+VN7Y/yFnlsxOqOyM\nK09jcu6YiMdKpo1HOXunau7KFwqQY8uIevzi734Fu9Ue8xrugAenKfI5xSNHcOhQ/OyiR7y+/yMu\nGn8qI/OLWPTzG9jYsJdzsuclXF6InqTlLwasVRWbOHN0e7C/aPypvHngY8I6sUVQaj1NFDoi9/k7\nzTZcgd5r83YVb4ZvZkYmQRV7QRdXwIszSsv/rPMXM2LB5Jjlj/CHArx/8FMuHHcKAJeccT5rD29N\nqKwQ0UjwFwPW6spNnDn6RAAm5owmx5rBp9W7Eipb426iyB49+LvjBf+gL2ZWT5vRgi8UO7FbrOB/\n8pyTMY3PjXisp7+++QIZe72McLZPWDuz5EQ+qNycUFkhopHgLwakspbDuAM+puUd7Zu/cNwpvLDu\nzYTK13maKHREDq4Okw1XIHZ/vTcUwBort4/RHDe9Q2VFJaFWX8RjRfYcaj2JzV34/SN/pMR1dDGY\nOYWTKG+tpt7TnFB5ISKR4C8GpMdefZaprlyUUp37Jnpy+e03lyS0DGONu4nCKC1/RwLdPns+2hJz\nvd9EEruteOpVtr7zccRjBY7Egr8/4GfbyvV874bbOveZjSZOHTmT1dL6F30gwV/0q7Vb19PmiT0e\nPpKn//Q4eU3dxyNceeaFoDUvr3wjZtkWbxvBtuhdLjtXbmD5M6/FvMbmf6zm0J4Yid1MVrxxcvu4\nXC6yMiNPNCuwZdPgbSEUjv3c4K/LnseWm8FZcxZ2239myYm8v2ddzLKReIN+KhPIaiqGPgn+ot+s\n3vwJp58wn2sf+n8JP6gFCIaClH22k29c9bVu+w0GA5PmzWTZyrdilt+wYwuHHniv218N3a7f5qEq\nRmAHCPr8ZDiij/bZuOJj3vzPJ2Jew+P2kJ2ZFfGY2Wgi0+Kg0Rd7Kci3Vr3LjAVzeu0v9WfxxxsS\n+yvoiGA4xG3v/C/nPnEHW/clN7tYDD0S/EW/8Pi8XPblK7n8/92AcXoBD37ydMJlX1q+DEuWg3lT\new/rHDuulP37Y0+O2n1wH47cyC1ugOzMbDyu2BO0Aj4/mc7owd9useFtjX0Nr9tDdlb0RVva/rmd\nT7fF7ropO1DGhAnje+1fPOd0dDDM8o1rY5Y/QmvNv6/5M75QgHE7FededTHBJCbNiaEnJcFfKXWh\nUmqHUmqXUuquKOf8Rim1Wym1USnVuykjhpRLvnMttgwHL/3iMf56wT28tu9D/rbt7YTKvrDsZWac\nemLEY5PHT6TyYEXM8vsrysgqiNzfD+3B3+uO/cA36AuQ6XBGPZ7hcBLwx+728Xu85GZFr4dnbx1b\nd8aeuFZ05lQuvPiiXvsNBgMT58/kqX8+H7P8EY9s/iefHN7On877Ea/88gmC/gBX33lzQmXF0NTn\n4K+UMgAPAxcAM4HrlFLTepxzETBRaz0ZuA34Y1/vK/rfa2vf5a2NHyRd7rcv/oUPXnmbd15ahslo\nIs+exd8u/jd+sf5ZVpRviFv+49Ufcv4550U8dvIJcwhET7kDQPmhSnILo6+Rm5uVjc8T+4Fv0B+M\n2fLPcGYQ9MUe6mnJz2BEUeRZxgDZBXmUxZnl6x1tZcHMkyMeW7RoEWtWrY5ZHmDZ3rX8ectrPHHR\nT8iyOrFZrLz+0lJe++tLPPXmS3HLd6W15tPqnexpqkyqnBh4UtHyXwDs1lqXaa0DwHPA5T3OuRx4\nEkBr/TGQrZQqTsG9RT9YuuZtJnxhNldedCl3r/ojX3z5R7y4a3lCK1c1elv5n1ce4YHf/oKZ46d0\n7h+fPZIHv/Btfrk+9kIqwXAIPbeQWy6/PuLx8888m5zrY//hWHXoEEVFkROqAeRkZeN3xw7+ufNL\nGVkcPXBnJRD8J3x7ESfOOiHq8fzCAioPRQ+iWmvKW2soySyMePzrl17DgQ074/b7X3/xl/jXcVd0\nW0h+wYyT+NF//4Rbb/wGh+prYpaH9txCz2x/h/P/fiffe//XXPny3cy5ehHvrk++cSAGhlSkdxgN\ndG2+VND+hRDrnMqOfdWRLvjZ7vbZi4aOB3ZOh5Os7KxuD/C01rhcLpqamwmHw2it0R2/BM6MDHLy\nclGAQRkwKIVC4fV6cLcd7acNBIMEw0GsNhuZ2ZmEtSasw+iOe7taXbQ1t2CzWrFb7NisVkxGE1ar\nBZvVRkiH8YcC+EIBAuEg1bU1lFdU4PV6UAYDdqsNh9VOcVERY0eWkGG2YzF2TysQDoeprDvMgUPl\nlNdUcajmMOZcJ/njR2A3WckwO8iyOBjhzKMkswiTwRjxf0Kb38Pn9fvZWrePbfUH+Lx+Pxv+8g7h\nsiZGTShl2rSp/PCm21kw46SI5QGWb1jLt+68nQOf7eDSW67h46XLycvK4f3yz3h86xv8cv1zvHTp\nf1KSGT2w3rf6T9x8+7f40enf6HXs7DFzuXnJ9/lo6mecOiNya3ZL3V6mnzOPSaPHRTxe7Mil2deG\nJ+jDbor8J0Cb182saTOi1nHWzFlM+OZZUY8DFFx1IqOLR0U9nuXMIOSPHfzdQR/2KCOOAEaMKOZw\ndcRfAQDqvc1YjWYyLY6Ix0+bNQ/76BzW7dnMKVMifyEuXfM23vpWrl54Ya9j/337fbz2xjK+/fB9\n/HPJo1HrsbepkmteW8IJBRP4ySk3cmbJbOpbm7jp09u58OzzmHLaiTz+qz/G/Gx9sOljfvPUI+zc\nsZOqfQdpPlSPszCbW/5wD7OLJzGzYAIz88dFzYXk8XnZuH8bZFlo9Xto9bvwhvxkWpzYw0bqy6op\nKRrJ+FFjGZlX1GsNBX8oQFvAQ5vfQ4urhV27duP1+9DhMDabnQybg6zMTMaVjsNqMGMxmjkSboLB\nEE0tzXh9Xjx+L36/n2A4hMFooHjUyM44Y1QGlFKEgiFaGpsxGY/+roZ1GIPBSG5+brc4g9b4/X7q\nampBKZRSGAztMctsMlNYdPSLX2vd/m8JBKivqyPcsX2sBmRun/mzT+LIv0uZDBSePY0RVx8NFhqN\nQtGwajeH/v7Z0S8FpVAK8s6YzKir54LWaHTnm12/YhfVf99w5CIogwKDouALUyi5Zj7Gji8KAK2h\nesV2Kv7+KToUQgfD6FAYrSHrrIkUXDMHg1JYjGYsBhMWo5mmtfuofHUjBpMR0ISDIULBMLkLJ5Bz\nyXTaAh4MKIwGIwalaF6xl5qn12OwGDFn2rFk2HFkOZl6zjzmX3YWnqCPVr+bFr+LqrZ69v7jYzxr\nD5JVlIfZYsZkbv/fZ5k9gtDcAqbllTIzfxxziibz1ennY5hzCys/XcunWzayccNGTltwKkXTxvLA\nrx7kqlMuINNix6AM1HmauPelh/jLd3/O+V+/gtUvv01x7tEP3Xml8zmvdD5/3vIq1y67n39c/jMK\nIoyhX7p3DVvr9/PWl34Z8f+r2WiisNHMg396iFceejLiOR9Wfc7CkbMiHoP2X7JRGYVUttYyKbck\n4jnTv/YFrp56TtRrFOXkYyiJ/kAY2odEWo3R8/9MmTyVaQ9cGfMa7oAn6nBTgFEjR7H98+h9/uWt\nNYyJ8UULcMvDd7MzeJhTohx/6LHfM/+CMzEZI/+qr/n7O5z30g94/+CnnD12bq/j63Zs4txrvshv\nn3yEG+Z8sXN/YVYey373LAeXVHHzPbdz2oJTOeWyxTzxv39k0uhxhHWYFr+bleUbeXr727z2/d8z\ndfYMTp57Mt+66ZssmDGHHQf3okdk8HndPl7YtZydjQcZm1lMicpl9a//TjAQJOAP0Hy4Hl9dKxkT\nizjnwZvJtDjIsjixGs20BtxUllXwwc+ex9fqJtDqRQdCGKwmzCMyKbnvPELhEGE0GWY7GWY74Ro3\nOx56E4Op/T0JB4OEgiEsBZlM+Ndz8YUC+EPBzjjjq2yi8r/fRZmMKKOh/cegsI/KZcbdFwMQCocJ\nowmHw7gqGtj78zfRYQ1dBpzZRuUw6d6LOhujKIUCvFVN7HrwTY4EPa01aLCNymbC3Re0n9tBKYW3\nsold/7UMHQwTZUBbQlIR/CuBrikSSzr29TxnTJxzOoU8sVtUnW4A/pzYqQDcCPw1ifOvB/6U5PkP\nRz+stcYXChDWYcJa4/qSG+cjNjLs0R8sdtVwWRNrNq9j54E9uD0evH4vwWCQc89YzOITF2Lu+Qte\nDHMmzoRr2jfrmxtZ8sj/8ErVh/z86ZfwBH1kWhyEdZgvTzmbXTt3MmFUadT7f/OESzlUW82scxbw\nyT9WMLboaMu4xt3Iv615lMcvvDdqixzgm1+7ifvuvhceinz8w6qtXDft3Jjvw9isIspba6IG/1h5\nfaB9kpcnGHnm7RG+UCBmegen2YbfGH2Mvta6PatnjPdi0aJF7DREH3Nf3lrLmMzYvaOnjZrF8vIN\n3DCjd8s+HA7z4evL+ctT0YekZlmd/HLRHXx/+W949+r/I9d29Etx58G9nHXOIi6+9opugb+rsUWj\neO+xV9hKkqhDAAAU90lEQVR87zZu+uF3OOtnXydn/jha/W4cZhsnFEzga9PP58ldP+n1fvb8S8Ef\nCrCrsZxth/cz5norNqsNh83OzIlTOX32ArJiDL3lm7/pfNnibqOxtRmr2UJWZiYGZcBqNHcf+vuD\n6JeK6N+TPP+eJM//cZLn3919M9qw5pi01n36AYzAHqAUsAAbgek9zrkYWNbx+lTgoxjX0+L4CISC\nut7drBu9rQmXCYVCev7VZ+v86SX6z0uf0VV11TocDusbXv8v/fOP/xa3vM/v06Zsu35n3aqI9Zn6\nl+t1vbs55jV+tOJ3+vGtb0Q9vuBvt+qy5sPR6xD069I/XR31eDAU1CWPXKXD4XDUc1x+j57w52ui\nHvcEfHpcjHtorfXWun367Be+F/X4jQ98X1/177fGvMbBlmp94hM3RazrM2+/om3FWToUCsW8htZa\n/2T1o/r2d3/Vfs3qSv3IK0/qrAlFevHNl8ct21W9p1nXe5p1IBRMqpzom464mVTs7vMDX611CLgD\neBv4HHhOa71dKXWbUupbHee8DuxXSu0BHgG+29f7ir4zGYzk2bPIscZoUfVgMBhY+9xbnLhwHj+8\n805GjxqFtTCTz9dt5Adzr4lb3mK2cNJ5p/HLx3r/ibS5di8lGYXk2SNPjDoi02ti/YZPIx7TWre3\n/KOkdgA6n7v4oyRmO5LRM1Zrymo04wsGOvthe2pytWA4FHseQJE9l1p39BQPm9dtwB6I/IzniDGZ\nRdhMFvY09R7++uLK1zjz8vN69X9Hcu+Cr7Oxdjfz/vZNxk+dxF0/voszLjybd/+c3KphebYs8mxZ\nUZ9NiYEjJX3+Wus3gak99j3SY/uOVNxLpJ/JaOK9x14B2nPPLP9sLdPGTe71MDua2274Bt/7zu3o\n3+tuAfbBX/8Ca8gLkRe/6tTweTlvLH0VvnZfr2MtfjcWoxm7OfZ4UIfJiivgjVjnxpZmPB8dhN7P\nrDsZDUZMBiP+cDDis4H9FWVU/HZVzO6CPFsmzX4XgVCwd5cdcLi8iisv7TlwrreFI2exturzbusX\nhHWYiukGnrwz8vOXnuxmK0uv+DlNvjZKr38EowTvIU9m+Io+sZgtXHDKIkqLRydc5saLv8y475zF\njoaybvs/eHs5J4+PPkrniFmTZ1BfFXl4YkXjYbK88ds0h/6whv3lZRGPVdUcpublTXGvESu5W11z\nAyZb9GcG0P4FkmvNpN7bEvF406F6TpoWf/GambbRPPNi9yG0n1bvwmm2MS0v+jOcnvJsWUzIHiWB\nf5iQ4C+OO5PRxLVnX87SvUfXonV7PVR/foAbL702bvn5M06krTryMoYr16xi32/fj3sNX2Uzh+si\nD7NsbmvBaIn/BbL/x0spjzJJq6G5CYs9zmw0oNCRQ12E7J7hcBhvXQvzp0ee6dzVrKxS3vvV893G\n+y/du5rLJp4Rt6wYviT4i7S4bNIZvLR7BVtq9wLwwvtLcRTnRB3f39WMcVMIefzUNTf0Onagspys\ngviLpFjslqjr4La42zAlEPwBmtoit9obm5swJxD8a1/dworVq3rt3162B4PV1G3IbTSnzpyLyWbm\nzY+XA/Bp9c6O4H963LJi+BqQ4/zF0De7YCK3n3glN735M2YVTGDHM6uYMT9+Kxfa/3KwFWSxbvtG\nLjr17G7HyqsqyC8siFLyKIvDRmNL5Ietba42TNbYXTYARouZVndbxGNNLU3YHNHH+B/hO9zC59s+\nhyu6729WXubf+6W45Y+YNHcG9z34U25v/SH5N8zjnlO+xsScxLvixPAjLX+RFkopbpp1MWuu+wOL\nx5zElpXr+fIViQe7CWeeEDEv/eHD1RQVx54YBWC122hqidxqb3G1YbLGf3htsphocUVOyawtBgon\nxA++hUWFVB7uvZB7bbCZ2XMTz394+aWXsePjzZy76GxWfuU3XBtnroQQ0vIXaWUzWbhp1sXcVB59\n1axIvnj7dZjyew8Jra2pYeHChRFK9Livw05za+RlEJ25mYxdGP/Bs8lqptUVueU/acEszhlvj3uN\n4uJiDh7s/dwgkdm9XT3w3Xt54Lv3Jny+ENLyF4PSmMwiylt6P7ANGMJMHD8hbvkzb7mUSQtmRjw2\ncuIY5nwpdu4fALPVQqs78iplroA3ap6arkpGjaahtq7X/vaEbokHfyGSJcFfDEpjMosob+s93HP8\nzWdw0QW9Ux30VDKhFGNm5ODsDfpjpnY44tz/uIHpCyI/p/AEoy8j2dW4krE01/d+8FzRWptUy1+I\nZEnwF4PSmMxiKiL0+dd6miMmnevJGWMRd28odlK3Ixx2JwEdeTUsV8CL0xy/22fhvFPbkxD2UN5W\nw1gJ/qIfSfAXg9KYzPbkbl2FdZh6TzMF9uhLJx7hNNtwRwn+vmB7eod4bCZL1EXc4yV1O2LSqFIC\nY7unbA6Ggnxyz7MUW+MPWRXiWEnwF4NSgT2bmg93U9dydKx/k7cNp9mWUKvdabbjDkZv+dsSCNx2\nY/QZvhX7ygi1xc4cCpBjzcAV8OLrkmdo24FdhOrdsbNYCtFHEvzFoKSUovXVHazdsr5zX7yEbl05\nYnT77Ni4jdpdsZdXhPaWvy9Ky3/1o0vZt2FH3GsYlIF8ezZ1nqMjj9Zt20RGsbT6Rf+S4C8GrdxR\nBWzaubVze2/VQZyuxPKa713/Ocv/EDlj5ablH7Hvk9gLqwNYY7T8vW4vOVnxu58Aiuw51LqPPvT9\nfPd28kfLKqeif0nwF4NWcclIdu7Z3bn99ltvseupxNaUVUFN7d7I6wl5PF4c9vgPa9c89TpL//Bs\nxGM+j5fczMSCf4Ejh5ouqZ1379vD6DEyO1f0Lwn+YtAqLS1l3/59nduVhyrJL4qfCwcgLzsXnyfK\nA1+fF7s98rq5XVlMZjxuT8RjAY+P3OzEum4qlm3kzVeXdW4fLDvI+PHjEyorxLGS4C8GrelTprFz\n07bObJaHDx+mOIHUDgB5WTn43VGCv9eHM4GWv81mw+eLfI2Ax09BgsEfV4Ad29ufD4TCIXIum8nV\nV8dZ1ECIPpLgLwate2/+PoVzJ/CHDe1993U1tYwaOSpOqXZ5ObkEvZH76/1eH05H/HWVHXYn3ih/\nPVhGZVGUl9hfISOKi6mpbp+t/NBnL1I0ooiLZkpGTtG/JPiLQSvD7uSd37/II1tf5bPqXTTVNTJu\n9NiEyhZm5xH0Rl7GsWj2OMZPjp8iwulw4PNGHs5ZdMfpjB2ZWL/9mFElNNTWs7pyM09vf5vfnXOn\nLKgi+p0EfzGojcks4n++8F2++94vCToNTJ0wObFyI0oY+b0zIx4rvWAOM2efEPcaDrsdv6938A+F\nQ3iDfuwJzBUAGDd6LPU1dXz//V/z0OLvU+SQYZ6i/0lWTzHoXTj+FNZWbeGxr9Ww8KT5CZXJsjsJ\nj3YS1mEMqnsbyBcMJDTDd/HF51IxLtxrv6cj8Pe8bjSTxk7A1dDCV6adwxdKElvTQIi+kuAvhoT7\nTr2RQkduQqkdoH39XKvJjDfo75V9s32Gb/zg77Q6CRDqtd8V8CSU1O2I02bN5fp/uYU7534l4TJC\n9JV0+4ghwWo08y8nfSmpvnKnKfIsX2/QhzWh3D7miLl93EFfUsE/w+7kyX/7DSbp5xfHkQR/MWw5\nzXZcgd7j9H2hQEItf1uUGb7VDbWEKiOv8CXEQCHBXwxbTrMNV4TkbhXvbSXsizwSqKtoWT03bdrE\n7sdWpqSOQvQXCf5i2Nr9p+Vs+OyzXvsrnvsEAr0f5PZkM1nxRWj5N7U0Y3EkNtJHiHSR4C+GLV9d\nG4dreq8Gpv0hcjLiPzhurmtkw7883Wt/U2szNkf8GcJCpJMEfzFsWR02Gluauu0Lh8PoYIjsjMy4\n5bOdmYQ8kVv+EvzFQCfBXwxbdoed5taWbvta3G1gNGAyxh8FnZORRTjQe6hna2srDkf8xHBCpJME\nfzFs2Z0OWnoE/6a2ZgzmxIZcOm0OCGv8ge6tf3OmjVETE0szIUS6yCQvMWw5nU5aW7sPyQwRpuDs\naQmVNxgMKLORxrZminOPJnGbfu58zAb51RIDm3xCxbC1+LpLUMbuK3+ZHFYmXH9awtcwmI20tLV1\nC/6ugJdRGfkpq6cQ/UG6fcSwNXpMCYbsHqkdgv6EFoA/4oSHr8ee032hdXfAg9MkD3zFwCbBXwxb\nDrMVd7B7Vs5EZ/d2XsNq77WIuyvgTSq9gxDpIMFfDFtOkw13j9w+3qAvqeBvN1nwhbrPBpbgLwYD\nCf5i2IqU28cb8ieU1O0Iq9HcK79P5c79hD3x00MIkU4S/MWw5TD3zuq5b+8+6j7dn/A1IuX3Wf/w\na9SWHUpJHYXoLxL8xbB1eG85a37xUrd92zZtZf97GxO+hs3YO79PwOMnP9HF24VIEwn+YtiymSw0\n7avuts/ldmGxJZ6Ube3PX+T919/uti/ok+AvBj4J/mLYysvMIejt3mp3ez1YbYk/rDWZTLS5Xd32\nBb0BCnLyUlJHIfqLBH8xbBXm5BPydn8w63a7sSUR/K02K25P94fGYW+AwpyClNRRiP4iwV8MW4W5\n+YR7Bn+PB1sS3T4WmxW329257fF7sU4uIMMmid3EwCbpHcSwlePMQofCuL0eHLb2GbnFk0soDpck\nfA2bzYbbczT4+3SQCT8+H4NB2lViYJPgL4Ytg8FA6ZKL8OkgR9rp4xbMIM+WlfA1bHZ7t24fV8CL\nQyZ4iUFAgr8Y1vInjMQT8pFL++It3qA/qRm+l3z7Wrzhow+N3QGPzO4Vg4L8bSqGNYe5e4qH9hm+\niSd2c1jt+MPBzm1XwIvDJMFfDHwS/MWw5jDbcAWPBn9f0I8tifQOPWf4Sl4fMVhI8BfDmtPUPcWD\nN+THZkp8tI/daMEXPDpiqKKqEl9FU4wSQgwMEvzFsObs0e2za9VGmg7VJVze2qPl/8nqD9n+4pqU\n1lGI/iDBXwxrW558n5XvLO/c3r70Q+oqDydc3ma04Ony5dHc2oLdIQu5iIFPgr8Y1oLNHqoPHw32\nAV+ATEdGjBLdbVn7KW/95K+d261trTiczpTWUYj+IMFfDGvZuTlUHqrq3A76/EkFf6fNTqBLfqDq\n6mpy8ySpmxj4JPiLYe38xeeyfs3Hndshf5AsZ2bC5TOdmQR9R4P/xrXrueTci1JaRyH6gwR/Max9\n66obqN1RTl1zAwAhf4DsjCSCvyODoL99nH+tq5FwvpUbL/5yv9RViFTqU/BXSuUqpd5WSu1USr2l\nlMqOct4BpdQmpdQGpdQnfbmnEKk0Mr+I/Ekj+ctrzwKQfdoEivITz8iZlZFJyN8+1PPDw9u4csmt\nSXUbCZEufW353w28q7WeCrwP3BPlvDCwSGt9ktZ6QR/vKURK3fm7/6BtXPvErJwrZlKYV5hw2Sxn\nJqGOlv/K8g2cNWZOv9RRiFTra/C/HHii4/UTwBVRzlMpuJcQ/eLiGWewqqJ96UZfKJDUDN/xY0uZ\n9Isr0FqzsmIjZ5VI8BeDQ18DcpHWuhpAa30YKIpyngbeUUqtU0rd2sd7CpFSs/LHU+dtpry1Bq01\nZmPi+Q4dJhs+HWR3UwVGZWBC9qh+rKkQqRP3U66Uegco7rqL9mD+kwin6yiXOV1rfUgpVUj7l8B2\nrfXqaPe8//77O18vWrSIRYsWxaumEMfMaDBy5ugTeevAx1iTaPUDWE1mfMEAK8o3cFbJHJRS/VRL\nIY5asWIFK1as6NM1lNbR4nUChZXaTntffrVSagSwXGs9PU6ZJUCr1vpXUY7rvtRJiGPx/I73eGn3\nCnY0HGTLjU/EL9DF+Ee/TMaaBu78xu3cfOZV/VRDIaJTSqG1Tqrl0ddun6XATR2vbwT+GaFSDqVU\nRsdrJ3A+sLWP9xUipc4sOZEP1q6m9YN9SZe1aCOfP7+KxZPn90PNhOgffW355wEvAGOAMuAarXWT\nUmok8KjW+hKl1HjgFdq7hEzA01rrn8e4prT8RVooowHCmmQ/fxPvuoCGV7bQuKsq/slC9INjafn3\nKfj3Bwn+Il3Gz53Ogc92JB38lVLMvfhM1i9b1U81EyK2dHT7CDFkfOmaq4+57HVfvjaFNRGi/0nL\nX4guqhtrKc5NfJLXsZYRIpWk20cIIYYh6fYRQgiREAn+Q1hfJ4GI7uT9TC15P9NLgv8QJr9cqSXv\nZ2rJ+5leEvyFEGIYkuAvhBDD0IAc7ZPuOgghxGAz6Id6CiGE6H/S7SOEEMOQBH8hhBiGBkTwV0pd\nrZTaqpQKKaVO7nHsHqXUbqXUdqXU+emq42CllFqilKpQSn3W8XNhuus02CilLlRK7VBK7VJK3ZXu\n+gx2SqkDSqlNSqkNSqlP0l2fwUYp9ZhSqloptbnLvlyl1NtKqZ1KqbeUUtnxrjMggj+wBbgSWNl1\np1JqOnANMB24CPi9kqWSjsWvtNYnd/y8me7KDCZKKQPwMHABMBO4Tik1Lb21GvTCtC8CdZLWekG6\nKzMI/ZX2z2NXdwPvaq2nAu8D98S7yIAI/lrrnVrr3bQvEdnV5cBzWuug1voAsBuQD0vy5Avz2C0A\ndmuty7TWAeA52j+X4tgpBkjsGYw6lsBt7LH7cuDIEnRPAFfEu85A/x8wGijvsl3ZsU8k5w6l1Eal\n1J8T+XNQdNPzM1iBfAb7StO+lvc6pdSt6a7MEFGkta4G0FofBoriFYi7gHuqxFgI/j6t9avHqx5D\nUaz3Fvg98B9aa62U+i/gV8A3jn8theh0utb6kFKqkPYvge0drVmROnHH8B+34K+1Pu8YilXSvkTk\nESUd+0QXSby3jwLyRZucSmBsl235DPaR1vpQx39rlVKv0N61JsG/b6qVUsVa62ql1AigJl6Bgdjt\n07V/eilwrVLK0rEW8CRARgckoeODcMRVwNZ01WWQWgdMUkqVKqUswLW0fy7FMVBKOZRSGR2vncD5\nyGfyWCh6x8qbOl7fCPwz3gWOW8s/FqXUFcBvgQLgNaXURq31RVrrbUqpF4BtQAD4rqz0krT/UUrN\noX2ExQHgtvRWZ3DRWoeUUncAb9PeWHpMa709zdUazIqBVzrSuJiAp7XWb6e5ToOKUuoZYBGQr5Q6\nCCwBfg68qJS6BSijfZRk7OtILBVCiOFnIHb7CCGE6GcS/IUQYhiS4C+EEMOQBH8hhBiGJPgLIcQw\nJMFfCCGGIQn+QggxDEnwF0KIYej/Ayt1e+thAVeNAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, F(x), label=r'$F$', color='#1a9850')\n", "plt.plot(x, -model.coef_[:,0] + model.coef_[:, 1], \n", " 'k--', label=r'$\\alpha$')\n", "l = plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some manipulation of the coefficients is required to reproduce the filter. Remember the convolution for the MKS is\n", "\n", "$$ p \\left[s\\right] = \\sum_{l=0}^{L-1} \\sum_{r=0}^{S - 1} \\alpha[l, r] m[l, s - r] $$\n", "\n", "However, when the primitive basis is selected, the `MKSLocalizationModel` solves a modified form of this. There are always redundant coefficients since\n", "\n", "$$ \\sum\\limits_{l=0}^{L-1} m[l, s] = 1 $$\n", "\n", "Thus, the regression in Fourier space must be done with categorical variables, and the regression takes the following form:\n", "\n", "\n", "$$ \\begin{split}\n", "p [s] & = \\sum_{l=0}^{L - 1} \\sum_{r=0}^{S - 1} \\alpha[l, r] m[l, s -r] \\\\\n", "P [k] & = \\sum_{l=0}^{L - 1} \\beta[l, k] M[l, k] \\\\\n", "&= \\beta[0, k] M[0, k] + \\beta[1, k] M[1, k]\n", "\\end{split}\n", "$$\n", "\n", "where\n", "\n", "$$\\beta[0, k] = \\begin{cases}\n", "\\langle F(x) \\rangle ,& \\text{if } k = 0\\\\\n", "0, & \\text{otherwise}\n", "\\end{cases} $$\n", "\n", "This removes the redundancies from the regression, and we can reproduce the filter." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 0 }