{ "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": {}, "outputs": [], "source": [ "#PYTEST_VALIDATE_IGNORE_OUTPUT\n", "import pymks\n", "\n", "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "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": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD7CAYAAABpJS8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYW+V9L/DvKx2ts2k2LzM2tsc2i806HnYoFIZAWpJm\nMeW2TZO0t0zztL2hSRvyJL0kXdKmcG/Tm6TNjd1AQm6SJy6kKQ0UauxQmgQwNg4Y7IDBYxvvnk0z\nmhnteu8fZ5kj6ZyjM9bMSBp9P8/jZ6RzdKRjefzTT7/3d95XSClBRESLn6fSJ0BERAuDAZ+IqE4w\n4BMR1QkGfCKiOsGAT0RUJxjwiYjqBAM+EVGdYMAnIqoTDPhERHVCqfQJmHV0dMjVq1dX+jSIiGrG\nyy+/PCyl7HTz2KoK+KtXr8aePXsqfRpERDVDCHHU7WNZ0iEiqhMM+EREdWJOAr4Qotdh32YhRL8Q\n4r65eC0iIjo3ZQd8IUQ/gEdt9vUCgJRyB4Co0wcDERHNr7IDvhbMB2123w0gqt0eBNBf7usREdG5\nme8afgTAqOl++zy/HhER2eCgLRFRnZjvgB8F0KbdjgAYKXyAEGJACLFHCLFnaGhonk+HaPa27vs3\nfP75hyt9GkRlm5eAL4SIaDe3AejRbvcA2FH4WCnlVilln5Syr7PT1cViRAvqL174Jr7x2o+QzmYq\nfSpEZZmLLp3NAPq0n7qdACCl3Ks9ph9AVL9PVIvGU5OVPgWispQ9tYKU8jEAjxVs22S6vbXc1yCq\nlGwua9yOJibREYo4PJqounHQlsjBeGrKuD2WZIZPtY0Bn8jBeHIm4EeTsQqeCVH5GPCJHEylE8bt\n6XSygmdCVD4GfCIH8UzC8jZRLWLAJ3JgzvDNt4lqEQM+kYNpc0knw5IO1TYGfCIHUxlzhh+v4JkQ\nlY8Bn8iBeaCWg7ZU6xjwiRxMa1l9oy+EaQ7aUo1jwCdyoNft20PNefV8olrEgE/kYCqdQEgJoMnf\nwC4dqnkM+EQOpjMJhJUAwkqAJR2qeQz4RA6m0gk0+EJo8AVZ0qGax4BP5CCeTiDsCyCkBNmHTzWP\nAZ/IwVQmgZASRIMvyBo+1TwGfCIHyUwaQcWPoOJHghk+1biyF0DRVrqKAuiVUj7osL+Hi6FQrUnl\n0mhRGhH0+pHMpit9OkRlKSvDF0L0AoCUcgeAqH6/YP+gtn+wcD9RtUtm0/B7FQQY8GkRKLekczfU\n7B0ABgH0WzzmAe1nD9e0pVqTzmbg8yoIKj6kc5m8JQ+Jak25AT8CYNR0v928Uwvwg0KIsYLHEdWE\nVC6NgMeHgNcPAMzyqabN66CtECIC9RvAFwH8kxCix+IxA0KIPUKIPUNDQ/N5OkSzlsxm4Pf6EPD6\nAACJbKrCZ0R07soN+FEAbdrtCICRgv0DAL6oDebeA2Bz4RNIKbdKKfuklH2dnZ1lng7R3Epl0/B7\nfQgqzPCp9pUb8LcB0LP2HgA7ACOzzyOlfAwz9X6impDOZeDzKEZJJ5Fhhk+1q6y2TCnlXiFEnxCi\nH0DUNCi7E8AmKeWDQoj7hBCDANrYlkm1JpVNI2Aq6SRZ0qEaVnYfvlUQl1JuMt0u6s0nqgVSyry2\nTIAlHaptvNKWyEZGa8H0mWv4LOlQDWPAJ7KRyqnZfMCjIGh06TDDp9rFgE9kQy/f+L0+BBQO2lLt\nY8AnspHKZgAgrw+fg7ZUyxjwiWykc1rA9ygIctCWFgEGfCIbViUdZvhUyxjwiWykzAFfH7RlDZ9q\nGAM+kQ094PtY0qFFggGfyIZeww8ww6dFggGfyMZMDV+B1+OFV3iM3nyiWsSAT2TDXMPXf+qtmkS1\niAGfyIa5LRNQSzv6NqJaxIBPZCNpuvAKUAdvOWhLtYwBn8iGVUknzZIO1TAGfCIbhQHf51E4aEs1\njQGfyEaqqIavGB8CRLWo7AVQhBCboS5d2Gu12IkQohfaMojaModENaEow2eXDtW4sjJ8LZhDSrkD\nQFS/X+AzWqDvsdlPVJVSpj58QM30WdKhWlZuSeduzCxMPgig37xTy/53A+pSh6Y1b4mqXmFbpt+r\ncNCWalq5AT8CYNR0v71g/5UA2oUQvUKI+8p8LaIFlcym4RUeeD1eAGpph22ZVMsWYtB2RM/stYw/\njxBiQAixRwixZ2hoaAFOh8idVDZj1O8BtUuHF15RLSs34EcBtGm3IwBGCvaPQC316I+9svAJpJRb\npZR9Usq+zs7OMk+HaO6ksmlj0jRAvdKWXTpUy8oN+NugdeBoP3cAgBAiom17zLQ/Aq2eT1QLUrkM\nfJ6ZRja/12e0ahLVorICvqlU0w8gahqU3antH4TavbMZQDvbMqmWpLLpopIO2zKplpXdhy+l3Gqx\nbZPFfgZ7qilqwM/P8FnDp1rGK22JbKSyaaMlE9D68FnDpxrGgE9kI5XL79JhWybVOgZ8Ihtsy6TF\nhgGfyEbhoG3A60Mml0VO5ip4VkTnjgGfyEYql8mr4fu0AVx26lCtYsAnslGY4fuNgM86PtUmBnwi\nG0VtmR41+LOOT7WKAZ/IRiqXRsBTnOGzU4dqFQM+kY1UNmPU7YGZhVCY4VOtYsAnslFUw/ewhk+1\njQGfyIbapZN/4ZW+nagWMeAT2UhnM3mDtj5m+FTjGPCJLEgpkSxqy9QyfPbhU41iwCeykMllISFt\navgM+FSbGPCJLKRyatkm4Cnu0tH3EdWasgO+EGKzEKK/1CLlXMScaomexee3Zaq308zwqUaVFfCF\nEL0AIKXcAXVlq16bx/UDuK2c1yJaSPrArFUNnxdeUa0qN8O/G+ri5IC6WHl/mc9HVBX01kvzlbZ6\nlw4vvKJaVW7AjwAYNd1vL3yAEKJX+wZAVDP0LN7qSlu2ZVKtWohB27YFeA2iOaXX6fPnw9e6dJjh\nU40qN+BHMRPQIwBGzDvdZPdCiAEhxB4hxJ6hoaEyT4dobuidOP68kg4zfKpt5Qb8bQB6tNs9AHYA\ngBAiom/TungGALRZDepKKbdKKfuklH2dnZ1lng7R3NCDesBiPnzW8KlWlRXwpZR7AaMLJ6rfB7BT\n2/+YlPIxbVvE4imIFsT2I7vxiWe/6rrDxrItc5ZTK0gp8VcvfAvfPvD0LM+WaH4opR/iTEq51WLb\nJovHFD2OaCFIKfF72/8WWZnD9d2XYPP5N5c8xmjLNF145fV44RUe1x8arw69ja/vexwA8N61NyAS\naJz9yRPNIV5pS4ve0YnTyGoLj+86td/VMcaVtqaSDqBm/G4vvHrx1AHj9u7Tv3B1DNF8YsCnRe/A\nyBEAQIMviINjx10dk7To0gHUvny3XTqD4yeNDwz9HIgqiQGfFr2jE6cBAP3n9eHw+ElXx6Qt+vD1\n+25r+MdjZ3Fh2yosCbfinYkzszhjovnBgE+L3vHJIbT4G3BB23kYSUy4qsFbXWkLqBm/29kyj08O\nYUVTJ7obO3F8ki3HVHkM+LToHY8NobupE0vCrQCA4eloiSNmrrQtLOn4PYqrtkwpJU7EhrCicQlW\nNHbieIwBnyqPAZ8WveOTZ7GisROdIbUz+Mz0WMlj0hZtmYD6AeDmG8JIYhyJbAormjqxoqkTJyeH\nkNMGjokqhQGfFr0Tk8PobpzJ8IfipQP+zHz4xSUdNxn+0PQ4AKAzFMGKxk6kchkMufhmQTSfGPBp\nUUtm04ilptEZjqAzrGb4Z10E3pnpkQsGbT3uBm3HkjEAQFuw2fRBMz6rcyeaa2VfeEVUzaIJNfBG\nAk3oCLYAgKtMO5lNwys88Hq8edsDXsVVW+ZoYgIA0BpsgqI9x2iCAZ8qiwGfFrVochIAEAk0wudV\n0BpswlC8dMBPZzNF9XtAnUBtOpMoebwe8NuCzcbVuiPxidmcOtGcY8CnRU0P+K1BdVqDFn8jJlJT\nJY9L5TJF9XtAu9I2WTrDH9O+WbQGm4yy0AgzfKowBnxa1IzAG2gCALQEGjCufQg4SWbTRS2ZgLqo\nuZtFzEcTMTT4ggh4ffBpc/Aww6dK46AtLWpRbfA0EjQHfBcZvk3A97m88GosETM+ZDzCg9ZgE0Y4\naEsVxoBPi5q5hg+oJZ1xFyWddC5jrGFr5ve6u/BqLDGBVu1DBgDag80YSTDDp8piwKdFbSwxCcXj\nRaMvBABoDjRgwmWGH7AYtPW7bMucSE2j2d9g3G8LtbCGTxXHgE+L2lgyhkigEUIIAEBEq+FLKR2P\ns6vh+73uZsucSsfR5A8b99uDzazhU8WVPWgrhNgMdW3bXinlgxb7B7Sba6WUny739YhmI6oFfF2z\nvwGpXAaJbAohJWB7nF1Jx+2FV7HUNBr9IeN+R6jFaNUkqpSyMnx9jVptofJo4Zq12tKHO7QVr3q0\n+0QLJpqYRCQwU0tv0YJ/qYFbtaRj0aXj9blaAGUyHTfKSADQHmxBNDnpevEUovlQbknnbqjZPQAM\nAigM6D2mbYOYWfCcaEFEk5N5g6czAd+5NTOZzVh36XjUK22dSkJSSjXD982UdPRzGE+Vbgklmi/l\nBvwIgFHT/XbzTinlVtOat70A9pT5ekSzEk1O5pV0mrSsezIddzzOvi1TLfM4deoksilkZQ5NppKO\nHvCjCQZ8qpwFGbTVSj17pZR7LfYNCCH2CCH2DA1xznCaW2OJ/Bp+gy8IQB1UdWJXw9fLPE69+JMp\n9bkbTYO2+jmMubjoi2i+lBvwowDatNsRACM2j+u3G7DVvgX0SSn7Ojs7yzwdohnJbBrTmYRx0RUA\nNGhZ91TaeT4cuxq+/iHgdLVtLD0NYObbBDAT8PULwYgqodyAvw0zdfkeADsAQAgR0R8ghBjQu3c4\naEsLSa/Tt5oGbRsUNcMvVdJR2zKtLrxSPwScBl9nMnyrgM8MnyqnrICvl2i0QB41lWx2mrY/IIQ4\nJIQoveoE0RwqvMoWmAnCpTJ82yttjQzfIeBrHybmQdsIa/hUBcruwzcNypq3bdJ+7gDQWu5rEJ2L\nmRkrTTV8RQ/4pQdtrUo6fqOGb1/SmUxpJR1TDb/ZH4ZHeIyFUYgqgVfa0qJlTJxmKukEFT88wuOi\nhm/TlqmVeZwGbWNawDeXdDzCg5ZAg7EgC1ElMODToqWXT8x9+EIINPiCRp3dipTSfmoFbY58p7bM\nmFHSCeVtjwQaWcOnimLAp5qQkzk88NJ38cVd/w+ZXNbVMWMWNXwAaPCFHFetysocJKTtbJmAOqhr\nRy8XNZlq+Op5NM0q4D/02hP47E+2lCw/EbnFBVCoJvzwrf/CV37+GABgVfMy/OZFt5U8ZiwRg1d4\nijLtUhm+Xp93ast0zPBT0/AKD4KKP297a6DR9RTJr5x9C597/iEA6tXBn77qt1wdR+SEGT7VhG/t\nfwprI91YF+nGv7z1nKtjogUzZeoafSFMOWT4evbu1JbpPGirzpRZ+LpqScddDf/bB55Gsz+M67su\nwXd+sZ1z8NCcYMCnqjcan8DeswfxgfU34T091+PFUweMDhwnhfPo6MK+IKZcZPjW0yO7GLRNT6Oh\n4FsFoLZmumnLlFLipydew40rLsNHN74bown1709ULgZ8qnovntoPALi+6xJctXwDJCT2DR0qeZw6\nj05xwC+V4evlGr9lH37pQdvpdNKYwsGsNdCE8dRUyTGId2JncGJyCNd1XYLruy+FR3jwX8dfcTyG\nyA0GfKp6z598HSElgMs71+GyznUAgFeG3ip5XOE8OrpSNfykQ4Y/05ZpX9KJZxKWc+1HtOsBSq24\ndWDkCADg8s51aAk04NKOHuw6dcDxGCI3GPCp6u06fQBXLrsQPq+ClkAD1rQsx+vDgyWPiyYnjSBr\nVqpLRy/XWAV8Y/I0hww/nkkibBXwjQnUnMtRB8eOAQDWt64AAFzWuQ6vDQ8iJ3OOxxGVwoBPVS2Z\nTePg2DFc2rHO2LY+sgJvR0+UPDaajOXNo6MrleE7lXSMydMcMvzpdBIhi5KO2/l03ho7jhWNncY4\nwGWd6zCZjmNw/JTjcUSlMOBTVXtr7DgyuSw2dqwxtq2LrMCR8VOOtfBUNo2pdP5MmbpGLcO3y5id\nSjp+F22Z05kEQgUtmcDMFb+lAv6bY+/g/NaVxv1LtTLWvqG3HY8jKoUBn6ra/pHDAIAN7auNbWsj\nXUjlMjgWO2t7nNXEaTp9QHU6nbQ81rlLp3RbZjyTQlixyPC18pLT9ArZXBaHoifyAv761hUIKn68\n6mKgmsgJAz5Vtf3DhxFSAljTvMzY1tPSBQAYHD9pe1zUmBrZuoYP2E+gZgR8q5KOi7bM6UzCsobf\n6iLDPzk1gmQ2jbWRbmOb4vFiQ9tq7HcxbkHkhAGfFszXX/1XXPDwb+LWR+/FIRc1eAA4MHoEF7Wt\ngtfjNbatbF4KADges18hLWrMlGnRh68F4+mMdYZv1PAdSjrOGX7Sskun2R+GgHC8huCE9nda0ZS/\nGNDGjjXYP3LYcS1d3XA8irt+dD96vvHruP9n30DW5VQUtPgx4NOC+MHB/8RfvfgILutch7PTUfzO\nf3yx5NWjUkocGD6cV84BgKXhVvg8Ck5MuinpWF94BcC2U8fpSlshBHwexbaGn5M5JDIp4zXMvB6v\nOmOmY4Y/DADoaujI276xfQ0mUtM4Pll6GdBP/uc/4OUzb+LmlVfg4defxJde3lbyGKoPDPg076bS\ncfzFC99E39IL8b1f/Tz+zy9/HIeiJ/D9N3c6Hndychjjqam8AVtAnWq4u7EDxxwyfD2Ltqrh6/X1\nc6nhq9sV27bMuPatwSrD18/HMeBPagG/sTjgAyjZjrrn9BvY+c7L+FTfb+Dh2z+DD66/Gf/4yg/x\n9thxx+OoPpQd8IUQm4UQ/UKI+85lP9WWnMxhPDk1qzLBt/Y/hZHEBO6/5iNQPF7csrIXF3f04NsH\nnnYsURgDtm2ri/ataFqC4w6DtnYzZQJA2KeVdGzmxHdqywTU1ky7ko4e8K1q+Pr5OM2nc2JyGJFA\nY9E3hIvaVsEjPNg/fNj2WADYsu/fEAk04sMb7wAAfO7ajyCk+PHXu77teJxZTuYwloi5Kh9RbSkr\n4AshegFjZauoft/tfqod2VwWW/Y9jiu/O4AN3/oQLnnkI/jqz39QMvDHUtP42is/xC0re9G37EIA\nalnkty66DQdGjjheMbt/5DAEBC5qX1W0b2XTEsfyRjSpzpRpXnVKF1JKlXTsa/iAevGVXTlK/9Zg\n1YcPlJ4i+eTkcFF2rz5fAGtbuowPQSsTySk8c3Q37jr/l42B6Y5QBB+77H3YfnQ3Xj7zpu2xum1v\n/hjXfO9juPiRD2Pjt34bX3jxEU7PvIiUOz3y3QCe0W4PAugHsHcW++kcZHNZjCenkMpl0OQPIawE\ni2ZmtHJ6ahQvn3kTe8+8ieHEONqDLdi09ALcsrIXIZ91Rgqo3TAf//GX8fOzB3Fj92W455I7sevU\nAfztS9/BvqG38bVb/8ToXin0jdeeQDQ5iT+98jfytr9v7Y343M8ewpODL+CKJedbHrt/+DB6Il2W\nE5F1N3bi7PQYEplU0TTEgD6PTvFMmUDpDL9UScfnUZDMWWf4+oeIXUmnNdiEIxOnLfcBag2/sH6v\n29ixBrtPv2F77LPH9iKdy+BXe67N2/57l9yJh19/En/70nfwz3f+pe3vyt/s+jb+8ZUfonfJ+fjd\ni38Frw4dwtdffRxPHd6Fr9xyLzYtvcD2tTO5LH56Yh9ePLUfZ6bG4PMqWN28DBd39ODyznVoDjTY\nHqubSscxnU4iJ3MI+4Jo9IVc/V6Te+UG/AiAUdP99lnunxMfeeqvcUKrfeq/HwJCu6/9RP52y31i\n5lFW2/OPLXgd4/UKXsf0C1t4DoXPqz9U8SgIeBQEFD9yUmIiOYXx1CTGEpOIJmMYT05BYubrdtDr\nx5qWLqyLdGNtpAvLGzrQGmxCMpvC8dgQ9o8cxt6zB436sN+joDPcipH4OLbsexzN/jA+uP5mfHTj\nu7FOu5wfUL/aP3rwP3H/z/4JPo+Cr936Sbx37Q0QQuBjl70P/7TvR/jzFx7G555/CH9zw0DRf85o\nchJb9z2O21dfZcyBo2sONOC6rovx1OFd+LOrP2z5H/v1kcO2HwYrm5YAAE5MDuW1MOrGEjHLi64A\nUw3fpksnpQXzgMcm4HsV2wy/3JLOyclhXLn0Qst9G9vX4F/f/glGExNoCzYX7f+PIy+hI9SC3oL3\nrMEXwsev2IzPPf8QfnJiH35pxWVFxz702hP4x1d+iN/ecDv++vp7jK6oj2y4A/c++2W8//HP4o97\nfx3/44oP5n24n50ew3cObMd333gGp6dGoHi8WBpuQzyTxKg297+AwPrWFbhiyXpcseR8tAabMJVO\n4NjEGRyNncHRidM4Mn6qaK2AoNeP9lAL2oJNaAk0oiXQgIDHh1Qug0wui0QmhXgmiXgmiUR25nZW\n5uCBgEd4IISABwJCiLzbHiEg4IEQ6piQwMxPYdpnNt8Vru2b/w4eMb/DqhVfAEUIMQBgAADOO++8\nc3qOrsYOKB6vEQT12qP+7zNzf+ZfzLhVsE//RzXuF2zP31bwGKvnlO6e1/ycmVwWqVwaiUwaQqhB\notnfgPOalqI12IS2YDNag03weRTEUtMYikdxKHoCrw0fwpOHXyi6gnRl0xL0Lb0QvZecj01LL8DG\njjUIeH1IZdN46fQv8P03duK7v9iOb+7/d1zfdQmuWb4REhJPH9mFAyNHcM3yDfjqLZ8oKjXcc+l7\ncHZ6DF979YdYF+nGf7/kzrz9X9n7GGKpOP60Lz+7192++mp89qdb8Fb0eN6FRgAwnpzCsdhZfOii\nd1ke262dy4nJYcuAr2f4VmbaMu0yfDWY231r0d87K3pJx6pLB1D/LfUxEHOrqXpsAtHkpGVJB5gZ\nuD0wcgQ3dF9acM5p/PjYXtzZc13R8wLAhzbcji37/g0PvPQd3Nh9ad4H7DNHd+Pzzz+MO1ZfnRfs\nAeCaro14ZvPf43/+7Bv4u5e/j++/uRPvX3cjmvxhvHL2Lex852WkchnctOJyfOH638PNK64wvi2O\nJWLYN3QIe88exN6zB7H96G5se/PHxnN7hAfLG9qxqnkpbl99FVY1L0OTvwFCqO/FcHwcQ/EoxhJq\nknNmahSpXBo+jwKfR0FQ8SOoBNARjiCk+BFSAgh6/fAKDyTUhEVKCQmJnFT/6LclJKSUyMnczH1t\nm3qMenxhImJO+sohIefsuWaj3IAfBdCm3Y4AGJnlfkgptwLYCgB9fX3n9Bn6xRt//1wOW5RS2TSG\n4+MYS8YQ8PqxNNxqWccG1JLFDd2X4obuS/H5a38H33vjGfzLW8/hSy9vg4TExvY1+PIv34v3r7vR\nMogAwGeu/hAGx0/iz1/4JlY2LcW7Vl8JQK2/f/P1J3H3BbcUtVXq3rX6Snz2p1vw9OFdRQH/gFar\nvrijx/JYPSie0r61FIomJ7Es3Ga5L+ziSluP8ECx+Tv7PG66dIrLTIB6ta2ExERquugaAbsOHd3G\njtUA1Pe2MOC/cHI/Yqlp3L76KstjA14f/mTT3fjkc/+AJwafx3vWXg8AeGP0KP5w55dwccca/MMt\nn7D8d24ONOArt9yLX1t3A/7vq/+Kr+97HJlcFisaO/GhDbfjoxvfbfmh2xpswk0rL8dNKy8HoCYz\nx2JnMZVOIKj40d3YYVs2o/lRbsDfBqBPu90DYAcACCEiUsqo3X6aP36vD12NHbZBw05nOIJ7e+/C\nvb13IZ5JQkBY1sYLeYQHX73lj7H5R/fjD3b+HR648WO4sH0VBrY/iLZgMz5z9Ydsj13e0I7LO9dh\n+9Hd+Hjv5rx9r2vdKHpWW2hZg1od1PvWC0UTMVzYZv2NUfF44fcoiNtm+GnLHnyd3+uz7cOfNko6\n9oO2gPXiLHY9+LqOUATLwm3YP3ykaN/TR3YhrASLPgjMPnj+zXjo9Sfxqf/6GlqDTZAS+KMffwkN\nvhAevv0zjuM4AHDreZtw63mbkMqmkc5lLMdWnAghcJ520RxVRlkFIynlXgAQQvQDiOr3AewssZ+q\nWEgJuAr2urAviEfu+DNsaF+Njz/7ZbzrsU9iLBHDlts+hY5QxPHY21ZfhZ+fPYiz02N5218ZegvL\nGtrRGbY+PuD1YUm41Ri7KTSWjFledGU+Z9sMP5exrd8D6hiIbVtm2nnQdmbGzOI6vv536S64ytZs\nQ8ca7B/J78XPyRy2H92Nm1debvu6gPpB9607PotIoBF3P/F5/LcnP4+QEsCjd/7lrBIEv9c362BP\n1aHsGr5WkinctslpPy0+neEIfvCeL2D70d0YTUzg9tVXYUm4teRx71p1Jf7X7u/hmaO78VtavV5K\niRdP7sc1XRsdj+1qaDfKIGbGTJk2NXxADcj2Nfy0bf0eUANeLDVtuc8YtLXJlvWs3mqpw5OTwxAQ\ntqUoALi4fQ2eO/ZzTKXjRtDdN3QIp6dGcPtq+29Tuq7GDjyz+e/x74dfhFd48CtrrrEdb6DFh1fa\n0pzxeRX8as+1+O0Nt7sK9oB6QdHKpiXYfmS3se3oxGmcnh7F1cs3OB7b1dhhGfBHtats20PFnSy6\nBqcMP5txrC071fBLl3Ts58Q/OTmMJeGI44fNtV0XIytzeNG0AtYTg89D8Xhx63mbbI8za/KHcfcF\nt2Dz+Tcz2NcZBnyqKCEE3rXqSvzkxKtG1qyvYXvNcucMf3lDB05MDhddEaq3A1q1LurCStA+w8+l\njbVrrfi9PqRtu3TU57QriekzZlqtenVyyvqiK7Orll2EoNeP546pa9xKKfGjQz/DL3VfZjlRHJEZ\nAz5V3PvW/RKS2TR+8NZzAICnDu/C8oZ2rLPo/DDrburAdCaB8VT+GrEjcTXgtzsFfF/Avg8/mykx\naKsYV+MWimeSCCp+235q/QIkq5LOCZurbM2Cih/XLN+IZ4/thZQSe868geOTQ7iz5zrH44gABnyq\nAlcsWY9LOnrwyP6ncCI2hOeOv4L3rr2h5EUoejdLYVlnNDEOAGgPtdgeG1aCxgBrIbVLx7mk49Sl\nY1fOAdSB02Z/uGjQVkqpTqtg06Fj9is912Bw/CReOv0LbHn1cbT4G3DnWgZ8Ko0BnypOCIE/uPwD\nODh2DO8oO91HAAAQP0lEQVT+lz+FAPARbfIvJ13GxVf5c+roNfw2hxJHyBfAlF3Az6VtJ04D1JKO\n0+RpTp0ygNqaOVZQw48mJxHPJF11y3xg3U1oCzbjd//ji3jqyC7cc+l72TVDrjDgU1V4T891+Nil\nvwafR8EXbrgHq0wrXNnRg+PJyfzr+Ubi4xAQzm2ZStChpJNGwCHD9ztm+NarXZm1BosnUDN68F0E\n/JAvgC23fQpN/jDu7LkOf3T5B0oeQwRUwdQKRICa5d9/7Udx/7UfdX3MklAEisdbVNIZSUygJdBg\ne6UsoE6v4DS1QqPfPmNWM3ybGn7aTYbfWLSu7YmY1oPvsh/+uq6L8eJvbnH1WCIdM3yqWV6PF8vC\nbUVX244mJhzr94B64VXc6cIrr/2FZ2pbpn1Jx64HX2e1CEqpq2yJ5gIDPtU0tRc/v4Y/kphw7NAB\n1JJOIpuynM8/mUk5dukEvD5kctmiSeqA0oO2gBrwC9syT00Ow+dRbK8sJpoLDPhU09SAn1/DH42P\nO/bgAzNXwsYzqaJ9yWzaOcP36guZF5d1ptMJBEuVdIJNGE9O5X1gHIudRXdjx7xPj0v1jb9dVNO6\nGztwamokL3i6yfCdVr0qNWjr0zp4rAZu3ZZ0cjKHCdP0DMcnh9DdaD+HDtFcYMCnmtbV0IF0LoOh\n6SgAIJ3NYCQ+UbI04rTqVSKbcgz4+j6r1kw3JZ1ObUK5IdOEccdjZ41FXYjmCwM+1TSjNVMb9Dwb\nj0JCYnmJwU+nVa/czKUDwHI+HTd9+Mu1qZ1PTamLwSUyKZyZHsMKBnyaZwz4VNO6tDKIXsfXF0RZ\n1mA/4yRgWvXKIsNPlsjw/UYNPz/Dz8mcWtIpEfD1czs9pZ6z/mG1wmFaZKK5wIBPNa2rUc2W9att\nT0+rWfPyRuflkxt81hl+JpdFVuYcB2317L9wXduENgBcaiGRpVrAP6UF/OOxswCAlY3M8Gl+MeBT\nTWsNNCGo+I2Lr/SseXnYOeDr0wIXrnqlZ+2lrrQFiks68RJTI+tCSgCRQCNOayWdY3rAZ0mH5hkD\nPtU0IQRWNy/DkYnTANQZJ4Nef8mpgkNGSSc/w09k1SzdsUvHpqQzMxe+c4YP6FM7q99KDo+fQsDr\nK1mGIipX2VMrCCE2Q12svFdK+aDF/gHt5lop5afLfT2iQusjK/Ha8CEAwGD0JNa0LIcQwvEYPcMv\nnEBN7613GrTVPwwK2zJn5sIvHfB7IstxYOQIAODg2DGsjXTbLhRPNFfKyvCFEL0AIKXcASCq3zft\n7wewQ1vmsEe7TzSnzm9diaMTZxBPJ/F29ATWlphHHzB36eQH/KSLDF9fHCWRyc/wSy1vaLYusgLv\nTJxBKpvGwbFjOL91ZcljiMpVbknnbqjZPQAMAigM6D2mbYPafaI5tb51BSQkXh85jGOxMyUXTgFM\nXToFg7ZuavgBRd2nfzjoSi1vaLYu0o2szGH/yGEci53FusiKkscQlavckk4EwKjpft5IWcEC5r0A\nthU+gVbyGQCA8847r8zToXrUu/QCAMB3f7EdWZnDehfZss+rwOdRihZBSRoB375LR9+XLKzha8/l\nJsNfrwX4R/Y/DQDYpP0diObTggzaaqWevVLKvYX7pJRbpZR9Usq+zk72IdPsdTd2YE3Lcjx68FkA\n7oOnOkVyfoafdJHh6+vV2mX4pS68AoAL21YhrATx6MFn4fcouHLpha7OmagcJQO+EGLA4o9epokC\n0FsLIgBGrJ8F/Rywpfn07tXXAAAu61znur0x5AsWXXilB3znQVs14CcKJl7Tvy24Ken4vAruWHM1\nAODmlVeU7N0nmgslSzoFZZlC2wD0abd7AOwAACFEREoZ1W4P6N07Qoh+bYCXaE59YtOvqytAzWJt\nV8sMP1N60Dbo1Wv45z5oCwB/fu3vYG2kG79xwa2uz5moHGWVdPQSjZbxR00lm52m7Q8IIQ4JIcZs\nnoaobGFfEB/v3Yyelq5ZHVOY4esLm5RT0nGT4QPqIut/3HuXceUt0Xwruw/f6huAlHKT9nMHgNZy\nX4NoPlgtc6i3WgaU0oO2hXPpz/Th2x9LVEm80pbqVlgJGmUYnZ7h69MnWFE8XniFxzLDDyp+LmJC\nVYu/mVS3rEo6SRcZPqBm8UVtmZmE63IOUSUw4FPdCvuKB23d1PDV/f7iLh0XUyMTVRIDPtWtsGKV\n4Zfu0tH3F5V00kljjh6iasSAT3UrrAQwVTho6+JKW8C6pBPPJFxddEVUKQz4VLfCviASmVTeAuip\nbBoe4YFSYubKoEVJZzrNkg5VNwZ8qlt6cDZ36iSz6ZLlHMC6pBPPJBFiSYeqGAM+1S09OJsXQUm5\nDPhBxW+Uf3Rqlw4zfKpeDPhUt2amSJ6p45dawFxn1aUznUmyhk9VjQGf6lbYIsNPZjMlB2wBuy4d\nZvhU3RjwqW7ZZfh+b+kZR4JKwHLyNLZlUjVjwKe6ZSxzWFDDd5oaWRfw+vJKOtlcFslsmlfaUlVj\nwKe6pWfjcVOGH8+kXNXhC0s6+kRqbqdGJqoEBnyqW1br2sZdDrwWlnT0slCQNXyqYmUHfCHEZiFE\nvxDivhKPc9xPtNBm2jLNGb7LgF9Q0jHWs2XApypWVsDX1qrV572P6vctHtcP4LZyXotorjXoAT8z\n+4Af8PqRzmWQzWW14/SSDmv4VL3KzfDvhrquLQAMAuh3eCxRVTFKOunCko6LtkxFHdhNZTPqc2SY\n4VP1KzfgRwCMmu63Fz5ACNHLdWypGvm9PigeL6bOqaSjrXqVVT8s9A8NXnhF1WwhBm0dF+wUQgwI\nIfYIIfYMDQ0twOkQzQgrgbwunYTrLh19XVt14FZ/DpZ0qJqVvMJECDFgsXlQr9tjJqBHAIwUHFsy\nu9fWxN0KAH19fdLNSRPNlbASNLp0cjKHRNZlwNdKOvr8+TMLmDPDp+pVMuBbLVJusg1An3a7B8AO\nABBCRKSUUQA9QogeqB8KbdoHwN4yz5lozoR8AaPDRu+6cRPwQwUtnTNdOszwqXqVVdLRg7fWhRM1\nBfOd2v7HpJSPadsi5bwW0XwwZ/j6NMluAn6Dkt/SOZmOAwAa/eH5OE2iOVF60pASrL4BSCk3WTzG\n6ZsCUUWETRn+TMAv3aVjTLymHTOZUgN+A2v4VMV4pS3VtXPO8H3FGX5ICZRcKYuokhjwqa6Zu3Ti\ns6jh6xn+lFbKmUxNo9EXmqezJJobDPhU10K+oNFDr2f4bubD0Wv4+iLok+k4Gv0M+FTdGPCproWV\ngHGV7GxKOuGCkk4sHWeGT1WPAZ/qWtiU4evB21XAL5iWYSoVZ4cOVT0GfKprDUoQ8UwSUkrEtHp8\ns4vA7fV4EVT8Rg0/lmYNn6ofAz7VtbAvCAmJRCaFWGoKAFzX4huUoFHDn0rH0cSAT1WOAZ/qml6a\nmcrEEdN66Zt87kozDb6QMfFaLBVHAwdtqcox4FNdm1kEJYnJ1DSCih8+F4uYA/kXbU1x0JZqAAM+\n1TW9Xj+RmkYsNY1mf4PrY9WLthKIZ5JIZtNoCTTO12kSzQkGfKprrYEmAMBYIjbr1spGfwixVBxj\niVjecxFVKwZ8qmutQS3gJ2OIpabQNIvWykigCdFkDGPJWN5zEVUrBnyqaxGtDDOWiCGWis8q4LcG\nGjGWiM1k+Az4VOUY8KmuFWX4syjptAabMJ6cwkh8Qr3Pkg5VOQZ8qmt+rw8NvqBRw2+axaBta7AJ\nEhJHJ04b94mqGQM+1b1WrRYfTUyiOTCLgK9l9IPjJ9X7DPhU5Rjwqe61BptwYnIY05kEOkMtszoO\nAA6Pn0JYCSLg9c3XKRLNibJXvBJCbIa6mHmvlPJBi/29UNe7hWm5Q6Kq0RFqwc/PvgUA6Ay5X4lT\nD/gHx95Be6h5Xs6NaC6VleFrwRxSyh0Aovr9Ap/RAn2PzX6iiupq6EA0OQkAaJ9Fht8WVIP8RGoa\n3Y0d83JuRHOp3JLO3VCzewAYBNBv3qll/7sBQEr5oGmRc6KqsdwUrFc1L3V/XEM7PEL9L9TVwIBP\n1a/cgB8BMGq6316w/0oA7UKIXiHEfVZPIIQYEELsEULsGRoaKvN0iGbv/NaVxu2VTe4Dvt/rM6Zm\nWG96DqJqtRCDtiN6Zq9l/HmklFullH1Syr7Ozs4FOB2ifJuWXgAA6GrsQFDxz+rYa5ZvBABcvXzD\nnJ8X0VwrOWgrhBiw2Dyo1+0BtGnbIgBGCh43ArXUA+2xVwLgwC1VlWUNbXji/Q+c04VT//umP8SH\nN9zBgE81oWTAl1Juddi9DUCfdrsHwA4AEEJEpJRRqMFdz+oj0Or5RNXmiiXnn9NxrcEm3LTy8jk+\nG6L5UVZJx1Sq6QcQNQ3K7tT2D0Lt3tkMoJ1tmURElVN2H77VNwAp5SaL/Qz2REQVxCttiYjqBAM+\nEVGdYMAnIqoTDPhERHWCAZ+IqE4w4BMR1Qkhpaz0ORiEEEMAjp7j4R0AhufwdOYKz2t2eF6zw/Oa\nncV4XquklK7mpamqgF8OIcQeKWVf6UcuLJ7X7PC8ZofnNTv1fl4s6RAR1QkGfCKiOrGYAr7TJG+V\nxPOaHZ7X7PC8Zqeuz2vR1PCJiKqZEKLXvOqfi/XAHfefi5rO8AvXyBVCbBZC9DusruW4nwio3t8j\nbXW4ASHEAzb7H9Aft8Dn5fi6lXi/tFX2pBDikPZni8VjFuz90mYUftR8foD9euAu1wuftZoN+NXy\nBjqcX1X9ws3mdes5oFXb75HpdfsB7NBmn+3R7hcaEEIcwsyiQwvF9nUr9X4BaJNSCinlWgB3AbD6\nnVqw90v7+5tfx3E9cBf7z0nNBvxqeQMdVNUvnNvXZUCrut8jXY/ptQa1+4XukVKu1f7tFpLT61bk\n/So4lz5tbY5ClXq/gNLrgZfaf05qNuBbqMgbaKfKf+Gq7j8oqiegVdXvkU5b+1kf2OsFsMfiYT0V\nKlk6vW5F3i+dljj8s83uSr1fFbOYAn5VqtJfuKr7D1rlAa1qaN+49poH/3RSyge1D8N2m29I86JS\nr+vSbdpyq0UqfN6l1gMvtf+clL3i1XwpsXi6lYq8gS7cZnfO+si7EOI2IUT/QmX6lXpdN0oFNO0x\n83ne1fp7pOuXUn66cKP2/2VUW0Z0BNbfkOaci9et9PtlWZKs1PtlUmo9cMv95aragF9i8XQrC/4G\nuvxQWvBfOKfzquR/UJfvV6UDWkX+I7ohhBgwfej1a/+e+nntwczYxloARU0C88Tydavk/Sr6HanU\n+6W1WPYJITZLKR+TUu4VQvQJ6/XANznsL4+Usib/ANgMYAzAZtO2Aah14AHTtped9s/zOfYAeKZg\nW0T72Wu6vQVqr+1CnJPl6xac14B2+76FOi/938d0u79S71e1/R7p74f2+35I+9lvc16bAdy3UOdl\n97qVfr+01+0BsKVgW8Xfr0r+4YVX80jLMD4tpfx907aXpbbIu561AuiRc3RhhcvzKnpdi/Ma1PYv\nzBWAM222o1C/Ydwl1Qy24u8X0WLBgE9EVCfYpUNEVCcY8ImI6gQDPhFRnWDAJyKqEwz4RER1ggGf\niKhOMOATEdWJ/w8ijirZIDfy9wAAAABJRU5ErkJggg==\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": {}, "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": {}, "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": {}, "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": {}, "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": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD7CAYAAABpJS8eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd828X9+PHXSdbwXnHsbMfOBhJIHPYI4BAIZRQCoRva\nkhY66Uih/bZ0U6DQFvi1kEJLS8soYW/isEfIJJA9nOXYiae8tKX7/SHJS8uJbMuO389H00h3n9Pn\njSK/fbrPfe6U1hohhBDHPkOyAxBCCDEwJOELIcQwIQlfCCGGCUn4QggxTEjCF0KIYUISvhBCDBOS\n8IUQYpiQhC+EEMOEJHwhhBgmUpIdQFcjRozQxcXFyQ5DCCGGjHXr1tVrrQt6c+ygSvjFxcWsXbs2\n2WEIIcSQoZTa19tjZUhHCCGGCUn4QggxTEjCF0KIYaJPxvCVUrO11uuj1C0CbMBsrfUdfXE+IYTo\nSx6Ph6qqKpxOZ7JDicpqtTJ27FhMJtNRv0bCCV8pVQ48AJRGqJsNoLWuUEqVxPrFIIQQyVJVVUVm\nZibFxcUopZIdThitNQ0NDVRVVTFx4sSjfp2Eh3S01hVAZZTqxQR69wSPKU/0fEII0decTif5+fmD\nMtkDKKXIz89P+BtIf4/h5wCNXZ7n9/P5hBDiqAzWZB/SF/HJRVsh4nh972oe2fJassMQw8CyZcvI\nzc1l+fLlLF++nKuuuoply5b12ev3941XNiAv+DgHaOh5gFJqCbAEYPz48f0cjhBH7m8rHuNA62G+\nNGNBskMRx7iSkhKuvvpqFi1aBEB5eTkVFRV99vr90sNXSuUEHz4BlAQflwBhkWutl2mty7TWZQUF\nvbo7WIgB9fIPH2DdDx9LdhhiGFi/fj1z5szpeJ6Tk0NJSUmMFkemL2bpLALKlFKLtNbLg8UrgTla\n6/VKqbLgTB6bzNARQ5G71ZHsEMQA+sX7D7GlYU+fvuaM/In8+oyvxT1uzZo13HLLLQAsX76cRYsW\nMXv27D6Loy9m6SzXWud2SfZored0ebxMa12hte67gSghhDgGVVRUsHbtWq666qo+7dmHDKrF04QQ\nItl60xPvDzabjby8PJYsWUJ5eXm/JHyZpSNEL/n9/mSHII5ha9eupbw8cKtSfyR7kIQvRFyZ00YB\n0GJvS3Ik4lhVWVnJ7bffTmNjI5WV0e5jTZwM6QgRx5hvnkFDbT0+g052KOIYVVJSwooVK/r9PNLD\nFyKOtgMNqFQTXu1LdihCJER6+ELEUXXXm6g0E7VfrqcwPS9+AyEGKenhCxGD2+MGQNs9VFVXJTka\nIRIjCV+IGFrt7R2P27o8FmIokoQvRAzN7S0dj1tllo4Y4iThCxFD16mY7XZ7EiMRInGS8IWIwZqR\nSuZnpgHSwxdDn8zSESIGZU4h44LJpM4dy7S5M5MdjhAJkYQvRAy19XU4Nx/GXJqPMhuTHY4YIPPm\nzQsru/rqq7nxxhux2+0sXLgwrP7aa6/l2muvpb6+vmM9+5C33nor7jkrKiqw2Ww0NjZSVlZGSUkJ\nOTk5cdsdCUn4QsSwY/t2mh5YjWXGSLZP2gKTz052SOIYZLPZeOCBB3jyyScBuOqqq1i8eHHYL45E\nScIXIobQtEzXllo2r/sEPpvkgMSAiNUjT0tLi1k/YsSIXvXou1q2bBnz58/veF5ZWSmrZQox0Nq7\nzL13OGQjFNF/ysrKuj3vy41PQqSHL0QM7Y7OqZhOpzOJkYhj2aJFi1i+PLCHVGNjY8detqHlkvuK\nJHwhYujWw5eEL/pJSUkJS5cu7Xje14k+JOEhHaXUIqVUuVJqaZz6JYmeS4iBNv20E8n/0VmgwCUJ\nXwxxCSV8pdRsAK11BWALPe9RXxmsr+xZL8RgZ8lOxzIpn8I/XMjp116c7HCESEiiPfzFgC34uBKI\n9D3k9uDfJVrr9QmeT4gBtXvrDuyr9pOdl4M2q2SHI0RCEk34OUBjl+f5XSuDCb5SKdXU4zghhoT1\nb67C9vB6HG/u4dOK1ckOR/QjrQf3jmZ9EV+/TstUSuUQ+AZwG/B3pVTYxFKl1BKl1Fql1Nq6urr+\nDEeII+Z0OkFBwzs72fH2hmSHI/qJ1WqloaFh0CZ9rTUNDQ1YrdaEXifRWTo2ILQFUA7Q0KN+CXCb\n1tqmlKoEFgF3dD1Aa70MWAZQVlY2ON9tMWy5nE6UOYUUswmPy53scEQ/GTt2LFVVVQzmTqfVamXs\n2LEJvUaiCf8JIHS3QAlQAYGevdba1vVArfVymakjhhqXy4XBZAwkfLcn2eGIfmIymZg4cWKyw+h3\nCSV8rfV6pVSZUqocsHW5KLsSmKO1vkMptTTYu88L9uaFGDJczmDCt0gPXwx9Cd94FSmJa63ndHl8\nR896IYaKOV+eT8q8cdQ9+THONtkARQxtcqetEDGYctPIVYWc+PMvsae1JtnhCJEQWTxNiBh2vrsR\n25p9ZKSl41H+ZIcjREKkhy9EDFuf/xCfx8uYjBHsXr8OPpfsiIQ4etLDFyIGj9tDisVMzSeV1L+x\nLdnhCJEQSfhCxOBzezCZTVgsFvweX7LDESIhkvCFiMHr8mCymLFYrWi3F79fxvHF0CUJX4gYfG4v\nZos5cEu7BqfbleyQhDhqkvCFiGHSzy7ivG8tItVqBaOiub012SEJcdRklo4QMehsM/kFI5h73VWs\nmtqMOc2S7JCEOGrSwxcihpoXNlK9cTdWUyDRO32yno4YuiThCxGF3++n6ZlP2bd+O/s37aLpX+vZ\nX7U/2WEJcdQk4QsRhd3lAB1Ylra5thHHh/s5XFub7LCEOGqS8IWIwtbWAkCq1Up6ahoArfa2ZIYk\nREIk4QsRRUt7ILmnpqaRnpYOQJu9PZkhCZEQSfhCRNFqD0zBTEu1khlM+O0OWSJZDF0yLVOIKHKL\nCii6ayHz5i/AUOdEpZlweWUTFDF0JZzwlVKLCOxtOzvSZidKqdkEtj9Ea7080fMJMVA82osh3Ux2\nRialY6cx6u6LOeHMsvgNhRikEhrSCSZztNYVgC30vIdbgom+JEq9EIPSnso9tDy9mbr9NViMJgBc\nPunhi6Er0TH8xQR69wCVQHnXymDvfw0EtjrssuetEIPe/v37aHt9J82HG3G2Omhctpo1b36Q7LCE\nOGqJJvwcoLHL8/we9XOBfKXUbKXU0gTPJcSAagteoE1PS8NsMOJcX82BPfuSHJUQR28gZuk0hHr2\nwR5/N0qpJUqptUqptXV1dQMQjhC9024PJPzMtAyy07MAcDidyQxJiIQkmvBtQF7wcQ7Q0KO+gcBQ\nT+jYuT1fQGu9TGtdprUuKygoSDAcIfqO3RGYc5+Rlt4xLdPpkIQvhq5EE/4TBGfgBP+uAFBK5QTL\nlnepzyE4ni/EUNDucACBHr4pxQQpBpxOR5KjEuLoJZTwuwzVlAO2LhdlVwbrKwnM3lkE5Mu0TDGU\nnHLJPEb9v0uZVFIKgHlkJgarKclRCXH0Ep6Hr7VeFqFsToR6SfZiSHH7vCijgVSTFYBpt13OyRNP\nTXJUQhw9WVpBiCjWv70K22MbMWoFgMVoxuWV9fDF0CUJX4godn+yDfvbe0gzB3r4+x96j/cffjHJ\nUQlx9GQtHSGicDqdqBQDBkOgX9ReWcdhmaQjhjDp4QsRhcvpQpmMHc9TzCY8LhnSEUOXJHwhonC6\nXBjMXRN+Ch63rKUjhi5J+EJE4dc+UlItHc9NFgte6eGLIUzG8IWI4ozvXkFeY+cCrzljR+Cvlz6S\nGLok4QsRhcvnxmI0dzw/5zuL2NywJ4kRCZEYSfhCRLH+sZW4PG4ILvlnMcqOV2Jok4QvRBQH1+3s\nmJIJsPHJt9jyzir4QhKDEiIBMiApRBRetweTpXNIx9HQir2yPokRCZEYSfhCROFzezGZOxO+xWrB\n7/ElMSIhEiMJX4goevbwrVYreP14fd4kRiXE0ZOEL0QUBquJjJzMjuep1lQAWuxtyQpJiITIRVsh\noii9dSEXlpze8XzkmCLMk/JpdznIy8yJ0VKIwUkSvhBRuHweLMbODU/O+kw5L2Xvwtzl7lshhpKE\nh3SUUouUUuVKqaVxjotZL8Rgc+Cet9i6onNXzlDyd/lkeQUxNCWU8JVSswG01hUEtjKcHeW4cmB+\nIucSYiC5PW4cG2tormnoKNvywcfU/rKCbTu2JzEyIY5eoj38xYAt+LgSKE/w9YQYFJrbWwFItVo7\nC70+vIfaaGpuSlJUQiQm0YSfAzR2eZ7f8wCl1OzgNwAhhoxQwrd2SfjpaekAtNntSYlJiEQNxLTM\nvAE4hxB9KjT1Mi01raMsIzWY8B0yLVMMTYkmfBudCT0HaOha2ZvevVJqiVJqrVJqbV1dXYLhCNE3\nnG4nxpHp5ObldpRlBHv47dLDF0NUogn/CaAk+LgEqABQSoUmKZcEZ/EsAfIiXdTVWi/TWpdprcsK\nCgoSDEeIvpE/aiSFv57POQs7L0uNLBiJ5biRWDJSkxiZEEcvoYSvtV4PHbNwbKHnwMpg/XKt9fJg\nmdypIoaM0NTLruvhTyotJf87p1Mya1qywhIiIQmP4Qd76BVa62VdyuZEOKa0yy8EIQbU4qVf5+7H\nHuj18Z9s3Ej9Xe+yf+uujjJrMPk7fb1fE//lD1dyzpc/g9/v732wQvQTWUtHDAv/u/Mhfvj5b/b6\n+Nq6Otw7G/C5OhdKs7e0cfiW11jx5Iu9fp0vXnkN7zzyElv27jiieIXoD5LwxTGvY5cq1fs27Y52\nADK6zNJJt6bha3LQYmvu/es0B2b0VNbs7/3JhegnkvDFMa/R1UrqyWOxjMiMf3CQ3eEAIDO9s012\n8LHD6ez162SMCszysRZk9bqNEP1FFk8Tx7y9h6vwNTlw17Xi9/u7bVsYjT049TIzPaOjLNVsBQXO\nI0j4xV86HRqqsRtk/R2RfNLDF8e8d957F/fOBrI/P4tGZ0uv2qSkWTCNyyY3K7ujzGAwoExGXEeQ\n8J3jrWiHlw0fbzjiuIXoa5LwxTHvwKGDAFimFVDrsMU5OuD4c8oo+Nm5jB89rlt51twJ5BUX9eo1\nWuxtHFq7m6YH1/DOSyuPLGgh+oEkfHHMqzlUA4BjQzVrP+1dT9vpdQFgSTF1Ky/9xjymXFDWq9f4\nZNcWGv+6CoCGOtn8XCSfJHxxzKutrQWg9ZktvPfee71q8+aTr1D3+7dIwdit3JJi6pz1E8eWys5l\nlG0NjTGOFGJgyEVbccxrqG8gJTsVb7OD/VUHetWmvuYwnqpmLCZzt/LNv3yO2jFFcP5NcV9j5/49\nAKTkpNLa0PupnEL0F+nhi2Ne8cLZzPrWhRgzrRyqqelVG6fThTIZwyu8Gmdb7xZP23dgHwAjJ43F\nYZMVNkXyScIXxzw9Jp3pZ5xEal4GDYd6tyKry+HEYA5P+CkWEx5X74Z0Dh48iDIZWPSTr1H4/bOP\nKGYh+oMM6Yhj3q5Vmyg8YS5ZI/NoruvdblUulwuDKfzHI8VswtvLhF+ycA62EhNTJk3G27AKu8dJ\nmskav6EQ/UR6+OKY5vf72fXnFVSuWM+lP/ky4398fq/aZY7OJ3vqqLByk9mE1927m6jsWTB17ky8\nte20Vexi18G9RxK6EH1OEr44ph2sPwRePyNHjmTShFJsBicenzduu+OvPouZNy0MK59w8nTy5kzs\n1bk/XbEaQ7WdtoMNtCzfxMZtm444fiH6kiR8cUzbUVUJwOii0fhqWml+djPbDuyO287pdWM2mMLK\nZ3/2HEZ99sS47f1+P7uWvUn1O1soDt68tbdaFlATySUJXxzT9lQFZsqMGz0GT307ba/uYP2WjXHb\nvXPH/9h0zyth5ZYUM05P/CGdA3U1aLeP0WPGUDq2GICqmuojC16IPiYJXxzT9h4MzLufOHo8UyZM\nAmDHvvg9/NbDTbibHWHlb9/7FJ/c+J+47T/dvRWA4rHjmTw2sAtozeFDvY5biP4gs3TEMW30rFLy\nbzqD2cfNwuUILHq2d/++uO18bg9pmelh5SaTCb83/u5V2/fsBGBycSn52bkoSwq1h2uPMHoh+lbC\nCV8ptQiwAbO11ndEqF8SfFiqtf5JoucT4ki4LBrL1ALG5hdhQIFRUXWwKm47n9uLyWIOK7dYrWiP\nL277XcG7bKdPnAzASX9czKyJJxxh9EL0rYSGdJRSswG01hWALfS8S305ENrvtiT4XIgBs/6D1Rg3\nN2E2mkgxpmDOSae+Lv7NVz6PF5M5POFbrRbwa5xuV8z2JWedQMH/ncuJU44HYPSYMTTr8CEiIQZS\nomP4iwn07gEqgZ4JvaRLWWXwuRAD5oOnKqh/qvMi7fn3Xc9xNy6I2y59ehFjpheHlVstgRunWu2x\nl0po0nYKSseQlRrYQKV9fRXrnnjjCCIXou8lOqSTA3RdBjC/a2WwZx8yG3giwfMJcURaGppIy+nc\npnBMXhE7muJPjxx17amcVnJ6WPmkmdNJnz8Jj449rPP+S29gcLd2PG/YsJe978kmKCK5BuSibXCo\nZ73Wen2EuiXAEoDx48cPRDhiGLE3tTJy4piO560bDrDx5VcD301jcPncmI3h8/BnnjqbbM9qDObY\nPzobn3qH1C4XffMLRuBtdeL1eUkxylwJkRyJDunYgLzg4xygIcpx5dEu2Gqtl2mty7TWZQUFBQmG\nI0R3rmY7uSPyOp47q200v7Obw02xx/Erf/Asq/77Wli5CQN+hwe7O/Z4vL2hhfzCER3PC0cWgl+z\n91D8C8ZC9JdEE/4TdI7LlwAVAEqpnNABSqklodk7ctFWDCSn24WvzUXXjsSEsYFvkZ/u3ha1ndvj\nxtfsBJ8Oq1u34gMO3fQS27dvj9Cys72n2UHhqM6tEEcXBdbl2Xmg8oj/O4ToKwkl/NAQTTCR27oM\n2azsUn67Umq3Uqp3yxQK0Uds7jZG/qqchV+4oqOsZEIxAFsrd0Rt12pvB8BqsYTVpaWmAdAWPCaS\nrft3gV8zdszYjrLxwcd7anq3AYsQ/SHhwcQeF2ZDZXOCf1cAuYmeQ4ij0ehqJaUwg0njOhc7mz5x\nKgC7D+yJ2q7FHrjYarWGL2WcZk0FYs/S2Rzc2nDiuAkdZeedcy6j7r2EsTMnHcF/gRB9S5ZWEMes\nj7d8QlvFLgztnatjnlAyDWVJoa45+h6zoWSempoaVpeeFrgQa3c6o7bPmzSKorsWctGFF3WUjcoa\ngTIZaXDIVocieSThi2PW6tVraFm+CaOzcyx+VP5IJv/1KkoXzInazm+A1FPGMa6kOKwuPfhLoD3G\nkM5huw1Dupni/NEdZbmWTFqe3MQbL4RfCBZioEjCF0PC4aY68qaN4aEXHut1m+pDgdUpp4zrvn59\nQVoOdQ5bpCYAZORkkXvdHOacPjesbsKECWRcPJX8MYVR279b8SYtz2xmRGp2R5nRYMSxuopPV4XN\nTI5q1eZ1ZJcW8e7Gj3rdRohYJOGLIeEfzz5K0/Zqfnzj93rdpra2FgyK8SPHdCs/+PgaPnz45ajt\nnL7A8seR5uGPGzOWrEumUzAhfDeskK0fbaT9zUqMhu574lpz0rDV937uwv/d/mtaKg/z2z/f3us2\nQsQiCV8MCSveDSxL0NbWjtbh0yUjqa+rx5SdisHQ/WPesuMQNZ9Gnx65bs0aqr/zPJ9+sC6szogB\nn81Bc0v0sXiHw4HBEj4fIj03i9bG3o/hpy8ILLy2Y9fOXrcRIhZJ+GJI2L5tO4ZsK/m/Oo8dTb2b\n2tjc0IQ1O3yJY7PVgtsRffGzdns7ePykRthwvKXBxuGbX+Ot51+P2t7pcGC0hH87yMrLwW5rjdAi\nQgweB5/Y95F6yjiqPt2F3x9/SWYh4pGELwY9l8+D5bqZfPnhWzBYUnjzQO/WpDn+poso/821YeWW\nVAtelztquzaHHeicc99VVlpgMTSHM/qdtk67k5QICT9/5Ah8/vhLKwM89NLjNFZs56TZsyHTzMf7\nt/aqnRCxSMIXg96ndbtx+T1cPesCMtc386df3Nardk1+O2NHjwkrt6Ra8Tijb1NotwcSfkZaeMLP\nDCZ8Z4xpmW6XG5M1fGnla5ZeT8Gvy3F6o/+yCfn3I/+m9dkt3Puz2xj5f+ex0ym7ZYnEScIXg96f\n/3YvTQ+t5cT8UvLbLOx4eQ21TfVx223/9zs0bw5fuyZnZD4p2eHDNSH2YO89lNy7yuzo4UdP+Kf+\nbBHn3nldWHlBWmDFkQZnS+zAgc3vr2fMiZM4cfRUciwZfFSzJW4bIeKRhC8GvXdeewOq2hiVVcBV\nl1wBXj/LnnkkZpuG5iZsr23DtrMmrO7Cb1/D2FuiL+uUO7qAtLMnMrJgZFid1WwBg8Llin4NwOF1\nkWEJ/3ZgqzxM4/0fse7T2ENSb274AOfhFuaVn4dBGTBVVPP3r/0mZhshekMSvhjUvD4v1ZsqKTlx\nGgBfvfRzKLORZ196Pma7nVWBWTiFI8Pny6emmLF7nFFn+4ybUUrO52cxZtToiPUjrj6RiafMiHru\nzf99m70rNoaVmzzg/LiG7ZWxZ908+L9/A3D91V8GYGLBWNr21LF1n8zWEYmRhC8GtddXv42v3c1Z\nZ54JBC6ajjlxEpvfj30D057qwEyeMYXh8+W3rVzH4T++TZsj8t2yDrcT7ddYIszDByi88HhGHT8x\nYh3AoXe3Ub85fKP04jGBtXUOBm8Ii2bD9k9JHZ3D2bNOBeCS8y8E4NFXnorZToh4JOGLQW35ay8A\nsOjCyzrK5l+8AF2QyqaDu6K22xdM+ONGhV+0dTS14t7ZQENL5JugXvnXU9Tc+ByeaFM36x3U19RG\nPbfP5cGaFr4Oz6Rgwq85HP0CrNPrxnXxWL73n993lC067xKUycjKt9+M2k6I3pCELwbMfcv/SWZx\nAc+/H30Oe08HHHWkTyvk3BM7txv8+fduJv/bp7GheXfUdjX1hwEoHh2+i1pGeuDCa0NL5OUVnMHx\n+ez0rIj1u+54nbceeCbquX0ub8SF18YXjgGDojbGJuobanfg9Lo5r6SsoywrLYP8KaPZuu7TqO16\nqrM3MfrM6XzrDxH3HRLDlCR8MSDqbA38YMl3aNtXz0+W301bnB2jQtpmZfK5v9zU7W7ZCVmFmA0p\nHGg9HLXdjPlzGXXfpZw0fWZYXWZGIOE3RUv4TicYVOACbQRGcwoed+RpnX6/H+32kpYefsNXijEF\n69gcPCr6XPwnlz9J/R/fJd/X/RfGOZ+9ED09F7sn+uygEJ/fx7ff+DMtZjd//ekdrN4ie+mKAEn4\nYkBcfuMX8DS1862//YL26en86O3/F3eJhAO2w+xrPkRZ0bRu5QZloOHu93jkF/dFbdvoaMFsNpNj\njTC1MtjDt7VGnh7pcrlQKdF/NIzmFDxRbtxqcbShrClkZISfF+C0O7/I1MVnRn3tLVu24t7VwKSi\n4m7lN17/DdIXTmF9bfSNW0Kuv+cW3t33Md///vdRJiNXXHuN3KkrgD5I+EqpRUqpcqXU0qOpF8e+\nR1c8wwePv8bMS87kvm/+iptP/iJPPPoYX7n12zHb/fXRBzn0g5cY0RJ+E5PFaKL+YPSx8LeWv4Lr\npZ0opcLqCkcWYRqXjYfIPW2304XBZIxYB5BiNkXv4Rth1J8+w6Vfj7xLel5qNo0x5uHXHKwmJcva\ncUdvyOzCqfidXio+fj9qW4Db/nUP/7zpTko2+Pntpd9h0Xe/wsE1O1h6z69jthPDQ0IJXyk1Gzp2\ntrKFnve2Xgwtmyq3ccLFp3Pitxfyv+1v4PF547bx+Lzcs2E5mTNG8/zfHwfgGzMvJXObnUd+/zdW\nrn8vatt33nsXvH4WzDk7rC63sIDW2uhLHG97/2PaPj4Yse7k006m4GfnMmZy+Pg+wJiTJlGw4Lio\nr51iSsEbJeE7PIHx/9SUyMNBB55bx3t/eDLqa9cfqiVtRHZYeY4lg5Y/fsCDv/pL1LZ1tgZ+/u2l\nZEwYwXN3B6Z2/uf3fyO7tJA//+J29tVGfj+68vv9vLr7Q65+4Rd88cmf8/qad+K2EUNHoj38xUDo\np64S6Hk3S7x6MYD8fj8fbFrL0nt+zU/v/z3bGvf1Kmm32Nu4/HtfYuaME9j06iqMxdnc9Na9zPy/\ny7jhtthf3B7a9CIHsu08/vxyJhQGZswYDUaefvgJ8Gn+8vD9Udtu27CZnEmjyEgNHw8fNWYUrqZ2\nvFHib21qIS078rBKWnBRNIc38iyccafPoPjKk6PGddwVZ1J8WVnEusq9e2hctpqqTZEvKHsb7DRt\nir74W/PhRnIK8yLWTTx+Moe27o06PHPvEw/ha3Pxm9t/T35WYGdRs8nMQw/9g9zPz+LeLdEvNAPc\n/u/7yC4p5JqffZO9LYd4+e9PsuDUeZx2zQXsOxz/l4XX56WyuZqXKz/ktif/yiur3oj67yOSI9E9\nbXOArnvF5R9hfZ+47e2Hefa+/3QWKIUCpp19EhPnTKetsYX3//MyWmu0JjB2rDXTzy9j7PGlNB9u\nYNWjr6MMCoVCGQwYjQZmXXAao6cV01hVy+pn3sDn8+H3+fH7/fj8PmZddjYFE0dRV3mQjc+/h9Fg\nxGQ2YTKZSDGZOPmSc8gfU8jB7XtZ//r7OJ0uXC4nHpcHj9vN7OsWkDtqBIc/2cOuNzaQkZlBVnY2\nOdnZmM1mzvjM+aRnZrDug9WsefN9GhsaaWm00W5rxdliZ87vriI/fwQNH+zCubueGTNmcPKsMkrG\njCczLYPS0lL2tRzinr/dx3uvv0XN1n34WgMX/cyT8/mXWoPJkEJqRQ2fm38FP73uu6QYu38k7n7s\nfn72g5txHmpm3KnT+Of/+zvnnXQGFfvX8rWvf537K17G7/fzwM/+GPbvsmHXJm758c0suP5KFpac\n1q2ubOpMskpG8t6KtyL+m7Y52mnaVc2piyL3EcaNHQdePzsOVDKjeEpYvcPWRt60yHPlmw41UPf7\nt3hHncqCJaeEn7ulBaM7+vWFktOOZ19L5OGk6toanOursTdFXhUzf0Q+vnYXLo8biyl8qMpQlE7J\nCVMjtj31tFPZ/OpHvL1xFeeedHpY/dPPPYMh1cSSy77YrfzKcxay2VzLA588x+emzWdOUfjr//Wp\nf3Lzdd+T++JAAAAfOElEQVTFOjKLb5x5Jbdd81N2nbmHa759Hav+t4JJr0zihp/dxJ9/9NtuF9D9\nfj//ePExfnfHbdQcrCb/5nMAqLvtLTz7bBhSTYyYMpYTyk7k4gsu5OqFV+Dz+9i6bStVh2v4aONa\nNm3eRNXeA+TPLiZv3mTanXa2/uoF0rMzycrLJjc/jxEFBZx8zmkUz5hMe2sb7z67AltLM802G62t\nbbhdLiaddxJFJ5XibGhl3b9ex2Q2Y7VYsFgtmC0WZp57MqOmjKehuo5VT1fgdrvxer34vL7Az/Ml\nZ1I4ZTxN+w6x4dl3MBqMGIwGjEYjSinmXHoO+eMKqdm5nw0vv4ff50ej0f7AZ+WUxfPJGT2Cg5sr\n2VKxJpCHlEIpUEpx+hcuIiMvi70fb2frW+uhyzWssivmce81P8Wg+veyasKbmCdKKbUEWAIwfnzk\nr9jxvLbjQ9a+/G7gSZc3sdLUSE7aAdyHWtj/zBvB84X+Dw5k28m1HsSxr5F9rwXHRv2BXwZaw66s\nFjJ843Bsr+XQ8ncD/4AGBYbA382TLGQZx9K8+QD7V36E9mu014/2+cHnZ2NOPZYpI7CvqcL22HoM\nJiMqxYjRbMRgSiGndj9WcyvVO7ew7/0N+BwetKuzR/SaZScpBem0vb6T1he2YcqyYs1OJy07kxFj\nCplSNJF27aK2+hC7X/mIdU+/RceCAykGRt93KQBNFeugqo2pp89i7ilzuejscoqKRtFgdrBh/1bu\nWPljbv3f+/zh1t9w5sXnU1RYyKlnncEmQw2PPPwAaLj9kftY+sVvdcQ2f8Jcdry0muJTZrDs1rs5\neeYcvnbJ57r9u1y95Ms0v7uTG++8JOK/2ynnn8WKB59i+/7dTB1f2q3u2XdeQXv9nH3GWRHbzimb\nwzNnFXOorZFI97y6Wx3k5uVGbJtqsuDZb6P2cORZPhV/eIzmww1w7QMR6121LTRURb55KnQhODMj\nM2L9yIKRoKGyeh/TJ0zuVtfucZD5tdksOuVLEdtecu5FPHTrX3h6xYthCd+v/VTV1VB62vEdG613\n9YM5i3nwgWXMv6ecQx/vxtzll80Hm9by3etuILUom+0bNjNuZOAO4+kTJrPxhfd48o0XWPLNb3Dv\nT25jxbp3+cPvfs+WzVt47/W3+GDF29h21mBMNzNn4Vl8/+wbmZ5fTM2sA7z41mt8+NEqdm3cxsp/\nPMP7Gz7ij02vAFDz3RfQ7sA1FGUyklaUQ+5xY5mYNQpl9nMgN4t2WyuNBw6zo8WBdnl5e/8GMhdO\nxdtgp/YPgam9ymzEmGrCaDbhK8mkbaKFloN17F67Bb/Xi9/jR3t8aI+Pdd59pLWNx13ZSP1/30Wl\nGFFGA8oYSMwtk6xkmQ/RsqmK/a99EPh5DuYDgN1FdlJbC7FvqKb2uVUEMzkE/kfD8VZS7fnY1lVy\n6KU1gTc3mEsAqo8zYh6dTfOHu6h/MXgndjAXHSpV3LP4lsAL9adAr/fo/gC3A+XBx4uApUdS3/PP\nnDlz9LHI7/f3+liHy6l3Vu3Rmyq36wO2w/pga51udrRqn88Xs53P59Nrtm3Ud/7nr/qHf7pVf+fO\nn+pHNr+mV+xdo+vammK2bXfY9Q/u/oXOLi3UgAZ0zldm64l/v0r/4d1/66bW5qhtdx/cq62FWTol\n26pXb/24o/yeJx7UgJ5//ZVR2z6+8jltOb5Q/+mVh8Pqbnv57zrjoil6/c5PI7bdWLtLj77/cv3q\nnlVhdQ6XUxuyLPqiGxZHjRnQi5d+PWJ94UklOmfK6KhxT19wsjbnZ0Ssu/2R+zSgH3z+0Yj1P/rz\nLzWgX3x/RVjdzsYDevT9l+undrwVsa3L7dIGa4qeffk5YXXrDm3Xo++/XC/f9mbUuH9w9y80oL/y\ni293lDU7WnVGcYE2pJr0ijXvRG3r8Xr0dbd+T0+683I9+v7Ldc5X52hAW4uy9Rd/eoOubaqP2lZr\nrQ811uqnV7+uH9n8mn5820q99N5f69/+40/6nY9XaZfbFbOt1lo3tDTpnYf26pq2Bl3dUqd37K/U\nrfa2uO1C4v38DGXAWt3LnJ1oD/8JIDSYWQJUACilcrTWtmj1w02kmSLRWM0WJo0pPuJzGAwGyqbO\npGxq+LzzeNKsqdx106+466Zf0eZoZ2/NAdwpmrEjihiZFrmXHFIyegJPPfM0l5x3IZ/98bX84Nab\nGUcOP7rhe1gKMnn8rgejtr363Eu4Y+nFbDaE97T3m5qZ8vmzOGnS8RHbjkrPR/s1++qrobh7XYvX\nTtEdF7H4zGsjts0Ljm+3B5dB7snn9mCKsJ59iNlixu+JPDbdHOzh52SFX3gFKJ1YgmlCDo2O8Jk6\nL654lcM/e53WR78Ak8Pbmk1mTv3+Z9EF4St9vlb5EUZl4PziyNcWAO783q088s9/8chdy5hcXIon\n38Q6334sCyfxi9O/SHlZ5G9TELiH4B+//DPNrnYOttWRdpWZnL9kkJMR+ea0ngpzC/js3PkdzxdP\nPa9X7ULyMnPIy8zpLMgccUTte+56Nlwl9C5ordcDKKXKAVvoObAyTr0YpDJS0zm+ZBqzx0+Pm+xD\nFp52Pg+98DiTv3AGd697gq/e+l3cje389Pe/7P5D2oNSivMnlLFy0we09ljX5u333mFmVnHUtnnW\nTA7d9CL/veehsLqGYDLNs0YeVslKywAF9vbIa+l43V5STNETvsVixe+JPKXTq32BnbmyI793Z59x\nFgW3zCOnOHxRt+27duBrsFNSOC7quS+74nKqMtrDbsC67Us/xvpKNTmWyBeqIZD0/nH/g/gdHv7v\nqzdx17J7aXA2c+93fs2Pv3Bj1HZdZVvSmZFfTHHe6F4nezF4JDyGr7VeFqFsTqx6cey5tvxKruVK\n2twOnp38Knuu3MPPv/qDuO0Kawzs+eFzPFD4b370hRuAwPTPLb96nuk/zIdLI7czGoyYs9I4VB2+\n/PEHa1bR+LePaJveCKXhbQ0GA2kzikgdETlheeP08C0WMzpKwj9+3lyKbr+QaZMjX3jNtwbOGWlN\n/L0HAguuzZocfUro9IxxtK05wIszVnL1mRcDsGrzOlp3H+aCyy6K2i7kM6eX89cn/0mrvY0rz7uY\n0tHFcduIY4d8zxF9KsOcyhfP+2yvkj3AdQsXo8xGnni6c276/1YElj6+aN78aM0C5yrIofFQ+Lo0\n23Zsx7mxhgxj9E1OJi+9kOM+c1rEuqILjmdqefShEYvVivb6I06PDE31jDYPP9uUTt3v3uS5fz4e\nVlcd5aarrqZljqXpwTU8+mRn+/seDXzL+ebnvhq1XVc3XPEVln7xW5LshyFJ+CKp8jJzGHPSZD55\nZ01HAn3rvbfBqLhy3sWx2xaOoLWuOaz8UHA1ytLg6pSRpKZYsEeZh593zmSmnzsnYh1A2QVnknv9\nXFy+8Juv3n1+BQ33fYhFRf7ybDVb8NXbqdoTPhc/cNNV7GGSyWMnYi3KZsOadR1lK199ndRROTHH\n4IUASfhiEFiwcAHu+jau+cn1ON0utm74lJySorhjxEVjRuFqbAvraR+uCyxdXDI6esLf9afXqbgr\nvJcN0Hq4CW2PfsNQ8dQSUueMweMPP+bArn24ttZGXXgNwJyVhq0xfGlmy8Q8Sk4/IWq7jvOfMJnq\nLXvw+/08tv5VDq3fzUnzwu8nEKInSfgi6X51482MnFnMM48tZ/7/vk/DzmqmnBh9R6mQU889k8yF\nUznU2n1/24a6eozp5phJ19vsDMy1j2D7/z3HqodfidrW3tCKc0stLfbwm6scdgcGc+xLY6k56bQ0\nhX8zSb1oEgu/cXXMtgAnn3IK3mYn8+5fwo/WPMD4S2Zz509/G7edEJLwRdKNGVFEzYbdPPbCcrQB\n8r9/Btde/7W47c479zwyF06lwd098fpMkDE+9rQ9k9WMO8oGJ9rrx2IJvws2ZNsHH9N4zwdU14Rf\nMHbY7RgssRN+Rm4Wdlv3mNtcdpocrYxKj38z+pUXBq5kV1cd5M6zb6Ty2dWcfnz0aw5ChCT9Tlsh\nIDBzZtFJ87ls5rmsnbeNU0bF7+EXpuXisznYUbOHEwo6p+NM/tJZTPSHLz3QldlqobUhvJft9wfu\nzLRYo1/wTQ3WtTrawuqcDicpMWb4ABTPnIrd0n2WzwefrKXmO8+z78+z4aQrY7a/9IwL+OdLj7Pg\nlHMZlR++0boQ0UjCF4OKyZjCaaMj32wVdmy7n8M3v8ZzrSVcObNzzZ1GZwsTs8L3su3Kkmql0Rk+\nw8fuCmzMYrVGHw5KS00DoC3CjVvm7FQyxsb+dnHRV69k/ycmtNYdN+Vt2b0dfJqJo6LPwe/q2oWR\nl18WIhZJ+GLImjZhEhgV+w90n/Gy5rdPoeadDguitx17XAmtKeGbmDS3B4ZarJboPfy04PaFbfbw\nHv4JXytnYpzdvHKtWXj8XlrcdrItgZVAt+/ZCcBxpdNiNRUiITKGL4asFGMKlryMbjdfeX1eWjdV\n42uJsgF50Bmfu5CRX54bXmFUZF8zkxNOjb51Q5o12MN3hPfwHV5Xx/LL0VS+/ymHfvwK67du7Cjb\nd2A/ALMmRb/pSohEScIXQ1rGiGwautx8daC2GvyakQUFMdulmayR18M3GUmfV8Lk46P3tE+YeQJ5\n3zmNMaXh0z5X3/UsW//zdsxz56Zn4W91sadqX0dZ9cGDGDOtslyB6FeS8MWQlls0gra6zp2vdlXt\nBaBoZPhaNV2t+d9K9n7v2bANOlraW/BUNeO1R/+GUFBQgPW4QixZ4csQ23YdwnE4/GJwVxNGB8bp\nD9R0biqSOX0UxQtPitlOiERJwhdD2pmfvYDMS6Z3bIi+tzowNDKmaHTMdgat8Le7aWrtnpx37tpF\n3W/fZMuHH0dt67W7cWyopupAVXid040lNfaQTnHwwuzBw51DUZY5Yzjr2s/EbCdEoiThiyHtzHln\nYSob1bExeJvXiak4l8kTSmK2y0gLXCxtaOl+x2ubPbCCZmicPpLmxiaaHljNxo/WhdX53F5SU8N7\n/l1NHhvYietQbeey0AeqD1CU2rvVSYU4WpLwxZCWZ8zAvaeJbVWBPWRHzhhPwc3ncPJJES7IdpGR\nEVigrLGl+0bo7cGEn54WPeFnWgO/LOyO8Nk4fpeXtPTobQEKcvLJOGU8qaMCS0fX2xrZ8f2n2LQ8\n+obuQvQFSfhiSLMfaKD+9repeCuwhWWjMzCtMrQMcTSZ6YGEb+sxpNMeTOKhbwCx2jp6JHy314N5\nQg6F42LfAwAw/TsLKDo1sMvJx7s2AzBxXPS1f4ToC5LwxZB2wqTpAOzcE+jhP/3X/9D4x/ewpkRf\nGgFgQkkxqaeMA7OxW3m7IzSkE31YJis9sLGK09V9ExK338uIH5/NeVfHH4vPt2bR0B74dvHJzi0A\nTCmJsM2VEH1IEr4Y0mYUT8FSkMnT/3qcprZmavYeQDfHnoMPcOJJJ5F73RxyR3W/K3bs5GJyvjKb\n0tIIO6cEhRK+w9k94du9gefR1sLvavPdr/DiD+7H7/fzxz/+EWVJYcEp8+K2EyIRkvDFkJZiTOF3\nf7odR7WNy7/1BVoam7FmRx+OCQkl5Z5r4mcV5pF22nhGFkRfoybdkkr+D89k1oXdN1DZtWc3tb9a\nydb34u/kmZaWisPWxq2P/Ima9bv4/A++zuRxE+O2EyIRCSd8pdQipVS5UmpplPolwT+3J3ouISL5\n4Rdu4PiFp7Hq/Q9pqmsgPTf+zUvVlQeo/s7zVLzwarfyQzU1uHc3YPDpqG0NBgPZ00aTNqL7RuV1\njQ14a1pR3uhtQ3Lz8/C0OHjKv4EzfvU5Hv7NvXHbCJGohBK+Umo2gNa6ArCFnnepLwcqgvvalgSf\nC9HnXnvkWY77+aV42p1k52bHPT47LQM8flraui9TvGrFu9Tf+S6OlvBlE7pyrjvI7o1bu5XZ2gJT\nQ0MXdWMZUVCAdvlwOBz8+3t/JMUoy1qJ/pdoD38xEJrXVgn0TOglXcoqg8+F6HOj80Zy+7wbMZfm\nUzoz/gJkuZmBXwptbd0XQHM4AzNv4iXt2sfWsf7l7tMoba2BH4WczPjfMHIzAuf/DMdTkh37JjEh\n+kqi3YocoLHL8267NwR79iGzgSd6voBSagmwBGD8+PEJhiOGs/kT5vLYo48xu3BK3GPzs/OAznn3\nIaELsbE2EgcwmFJwu7qvthn6tpDdi/VwfvP9n2FJs3LX9T+Le6wQfWVALtoGh3rWa63DrmZprZdp\nrcu01mUFcRa8EiKeyyadybjM+JuC5IV6+O3dE77LFbiIGy/hGyMkfHNGKpYZIykqLIp7/gmFY/j7\nz+/GbIo9fVSIvhS3hx/sgfdUGRq3B/KCZTlA5E1CoVxr/ZOjC1GIvmcxW8g8t5SRU8Z2K3c5nZBi\nwGCI3RcymlPwuLsn/OKZU8j/7ulMKZX59GJwipvwewzL9PQEENpMswSoAFBK5WitbcHHS7TWdwQf\nlwd/UQiRdOO+fBrjgzduhcy44BR2ZbfEbZtiNuHp0cMPLbecGuemLyGSJaEhndAQTXD2ja3LkM3K\nLuW3K6V2K6WaoryMEElhNZho7bFrVfaEAkbMjT+3YO5NlzHrhgu7lb3y32c5/NPXUJ740zKFSIaE\n54JF+gagtZ4T/LsCkCUAxaC09eZnaJ66gb9ccFNH2YFte3Dtr4/bNn98Ie2e7nfaNtU34Gt0kB28\nE1eIwUbutBXDVorFhMvRPWmvfux19j0Yf9XKuvV72PvmJ93K7HY7ymyUOfVi0JKEL4Ytk9WM29l9\naQWPy0OKOX7Crlyxgd1PrepW5rA7MPSirRDJIglfDFumVAtuR4+E73ZjNMVP2maLGZ/H163M4XBg\ntJj6NEYh+pJ0R8SwZbFaaW/sPiMn0MOPn7RNZjN+T/f9cHNLRzGC8E1RhBgsJOGLYWvqvJPYWbWn\nW5nX48FkiT+t0mK14Hd3T/ilnykj0x59WWUhkk0Svhi2TlxwOvX7uyf3qdfNI88af5aN2WzB7/V3\nK7N7XaSlxN7AXIhkkoQvhi2jT9HS0H1P29TiPEbnjI3SotP5X72MvbNUt7I3b/kXGXlZcNnv+jRO\nIfqKXLQVw9aH/36F3T98tltZzeqd2HbWxG2bl5ePyrXi9XdeuHW1tGPwx2gkRJJJwhfDVnpaGvg1\nbY7OBdR2/+Mdtr60KkargENb99L60rZud+p6nR4sqTKkIwYvSfhi2EpPD2yF2NDSOazj9/iwWOLv\nSXtg825aX9hGY0vniiE+twdrWvTNz4VINkn4YtjKCG5y0tTaPeGbrfETvtUSSOwtXXr4PqeX1FRJ\n+GLwkoQvhq3QrlYNXXrpupc9/LRgTz60nr7WmrRTx1E6a3qsZkIklczSEcPWjFnHk/nZ47BkBJK3\n2+MGv8ZqjT8On2oJHBMaw3f5PGR/bhannHx2/wUsRIIk4Ytha+r0aWQumIwlKzCW7/H7GHHLOZx5\nzvy4bdNT0wBodwQ2O7d7nGi/lrXwxaAmQzpi2DL6FN66dhqaA0M6HnyYJ+RSNDr+FoWnn3c2RXct\nZPzUwNr5lfv3UnPjc6x+7q3+DFmIhCSc8JVSi5RS5UqppXGOi1kvxECr3r2P2p+vYNXbgeWQG22N\ntL+7l8b9h+O2zUhLx5Buxktg4n1jcKZPRnDmjxCDUUIJP7g5eWijE1voeYTjyoH435OFGEA5GYGN\nzFvbA+Pw1YcP0fzfj9nzyY64bW2HG2h5ZjM7twWObW5tBiBTNj8Rg1iiPfzFBDYyB6gEyhN8PSEG\nTF52YDO2lrZWoHM8PjQ+H0t7Uyttr+2kctduAJqCPfzsTEn4YvBKNOHnAI1dnuf3PEApNVs2LheD\nUW5moIffFuzht9oDUyx7M5c+Iy3wS8Ee/CXREnyN7IysPo9TiL4yEBdt82JVKqWWKKXWKqXW1tXV\nDUA4QgTkZwV6+O3BufRtjkDSTrPGT/iZaYE5/A5nYP37rJG5pJ9fyoTxE/ojVCH6RNxpmUqpJRGK\nK0Pj9nQm9BygoUfbuL374CboywDKysp0b4IWoi9kp2eSfc0sJp48A4B2e6C3npEa/8JrZlrgGHsw\n4Y8oHkX2VScwcUJx/wQrRB+Im/CDCTmaJ4Cy4OMSoAJAKZWjtbYBJUqpEgK/FPKCvwDWJxizEH3C\nYDBQUD6dvEmjASidOZWCX57PrLKT4rbt6OE7Agm/ub0Vv9OL1Sjz8MXgldCQTih5B2fh2Lok85XB\n+uVa6+XBspxEziVEfzDUOqjeXxV4bE7BVJTZMbYfS9GIQkb9v0s565qLAHj10Wc59P0X8dpdcVoK\nkTwJj+FrrZdprSu6fhPQWs+JcEyp9O7FYLPvT2/w9kPPA7Bz+07aKnbhaLXHbWdJMaGMBtx+D9A5\nHJQXvC4gxGAkd9qKYS3FYsJlDwzLbN24iZblm3C2tsdpBUaDkdYnN7F2xQdA8MKvQZFukdUyxeAl\na+mIYS2zIJf6qsCdtaHx+Izg+Hw8be/tZU/BNgD279mHOTcdg0H6UGLwkk+nGNZOPGUObfvr2VOz\nH6fTCUBWLxO+wWTE5XTh9/vZ//EOik+c0p+hCpEwSfhiWLuk/ELQ8MjLyzvn1Kf3PuG73W52NVWR\n/pmpXPmFxf0ZqhAJk4QvhrXPLbiCghtOwznOitMVmGGTlda75REMphTcLhfr6naQflYxX7388/0Z\nqhAJkzF8MazlZmRzxoJ5fNK+j1Ovns/HE5qxmuPveAWQkmrGj+bpV58n06MozRnTz9EKkRjp4Yth\nb4ou4P3/vozN3UbaiPhz8EPOvOvLnPbDK3jh9w/jem4nSql+jFKIxEkPXwx7eS0mbE9t4n/2f+DJ\nTIGv9a6dxWhme+VOXHWtnHLdqf0bpBB9QHr4Ytj74sJFoODgK5/QtGJbr9tVPruadb9+CoArF1za\nX+EJ0Wck4Ythb8yIIrImFgJgSDH2ul3j5gP4atsxpJq47KwF/RWeEH1GEr4QwIy5MwGOaBzeZA4s\nlFZ0/ETMJlk0TQx+kvCFAOafG9iszVnX0us2JrMJZU3haz/9dn+FJUSfkoQvBHDD4uswT8onZ/Lo\nXrcxWywY0kxceqps1yyGBpmlIwQwKqeA0365mHxz76dlapcXX4uLmQWl/RiZEH1HEr4QQX+98McY\nVe+/9P7rbw+xYdsnmI2mfoxKiL4jCV+IoBn5xUd0/Nxps5g7bVb/BCNEP5AxfCGEGCYk4QshxDCR\n8JCOUmoRYANma63viFA/m8AG53TZ31YIIcQAS6iHH0zmaK0rAFvoeQ+3BBN9SZR6IYQQAyDRIZ3F\nBHr3AJVAedfKYO9/DYDW+g7ZxFwIIZIn0YSfAzR2eZ7fo34ukK+Umq2UWhrpBZRSS5RSa5VSa+vq\n6hIMRwghRDQDcdG2IdSzD/b4u9FaL9Nal2mtywoKCgYgHCGEGJ7iXrRVSi2JUFwZGrcH8oJlOUBD\nj+MaCAz1EDx2LiAXboUQIgniJnyt9bIY1U8AZcHHJUAFgFIqR2ttI5DcQ736HILj+dGsW7euXim1\nL15MUYwA6o+ybX+SuI6MxHVkJK4jcyzGNaG3ByY0LVNrvV4pVaaUKgdsXS7KrgTmaK0rlVK24FBO\nfqRpmz1e76jHdJRSa7XWZfGPHFgS15GRuI6MxHVkhntcCc/Dj/QNQGs9J0K9DOUIIUQSyZ22Qggx\nTBxLCT/WtYZkkriOjMR1ZCSuIzOs41Ja64E4jxBCDGtKqdldbz7txbI0MeuPxpDu4fdcqkEptUgp\nVR7jJq+Y9ULA4P0cBW9SXKKUuj1K/e2h4wY4rpjnTcb7FbzZUyuldgf/PBDhmAF7v4ITW57sGh9E\nX5aml8vWHLEhm/AHyxsYI75B9YE7kvMO54Q22D5HXc5bDlQEJ0GUBJ/3tEQptZvOe18GStTzJuv9\nAvK01kprXQpcBUT6TA3Y+xX87+96npjL0vSi/qgM2YQ/WN7AGAbVB66355WENug+RyElXc5VGXze\n0/Va69Lgv91AinXepLxfPWIp01pH+swk6/2C+MvSxKs/KkM24UeQlDcwmkH+gRt0P6AMnoQ2qD5H\nIcElSEIX9mYDayMcVpKkIctY503K+xUS7Dj8L0p1st6vpDmWEv6gNEg/cIPuB3SQJ7RBI/iNa32k\nlWeDK9JWEFiwcKB+USftvL00P3jXf5gkxx1vWZp49Udl0O5pG2cNn0iS8gb2wvxoMYeuvCul5iul\nygeqp5+s8/ZGvIQWPKY/4x6sn6OQcq31T3oWBn9eGoN7TzQQ+RtSn+vFeZP9fkUckkzW+9VFvGVp\nItYnatAm/Dhr+EQy4G9gL38pDfgHLlZcyfwB7eX7leyElpQfxN5QSi3p8kuvPPjvGYprLZ3XNkqB\nsEkC/STieQfJ+xX2GUnW+xWcYlmmlFqktV7ei2VpotUnRms9JP8QWJStCVjUpWwJgXHgJV3K1sWq\n7+cYS4AVPcpygn/P7vL4AQJzbQcipojn7RHXkuDjpQMVV+jfp8vj8mS9X4PtcxR6P4Kf993Bv8uj\nxLUIWDpQcUU7b7Lfr+B5S4AHepQl/f1K5h+58aofBXsYP9Faf6NL2TodXGso1GsFSnQf3VjRy7jC\nzhshrspg/cDcAdg5zbaRwDeMq3SgB5v090uIY4UkfCGEGCZklo4QQgwTkvCFEGKYkIQvhBDDhCR8\nIYQYJiThCyHEMCEJXwghhglJ+EIIMUz8f4nzIzu0GqzpAAAAAElFTkSuQmCC\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." ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 1 }