{ "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 `LocalizationRegressor` is in essence just applying a filter." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.ndimage\n", "from sklearn.pipeline import Pipeline\n", "\n", "from pymks import (\n", " PrimitiveTransformer,\n", " LocalizationRegressor,\n", " coeff_to_real\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Construct a filter, $F$, such that\n", "\n", "$$F\\left(x\\right) = e^{-|x|} \\cos{\\left(2\\pi x\\right)} $$\n", "\n", "Try 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": 7, "metadata": {}, "outputs": [], "source": [ "def filter_(x):\n", " return np.exp(-abs(x)) * np.cos(2 * np.pi * x)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5xdZXno8d+z77PnnsxMksmFCQkhEO4MoUhRrCASsUoPFLDo8dRKi9a21rbH4l1brReqpVp7OC1S0IoKejSCyB1vCISrSUgCIQm5Z5LMZK77/p4/1lp79uxZa++1Z/Zcdub5fj75MHvtd8162YRnnnnWs95XjDEopZQ6/gVmegJKKaWmhwZ8pZSaIzTgK6XUHKEBXyml5ggN+EopNUeEZnoCpbS1tZmurq6ZnoZSStWUZ5555rAxpr34+KwO+F1dXWzYsGGmp6GUUjVFRHa5HdeSjlJKzREa8JVSao7QgK+UUnOEBnyllJojqhbwReQcEXmxxPvrRGSjiGwVkZuqdV2llFL+VCXgi8jNwINe309E6oFvAJcAa4DLReScalxbKaWUP1UJ+MaYDwPnlhiyFnjWGHPAGJMB7gbWVePaSiml/JmuGn4ncKjgdQ+w0G2giNwgIhtEZENPT8+0TE6pSvxy74vcs+2xmZ6GUhWbzgevckWvI26DjDG3ArcCdHd362L9ata55iefBOBtKy4kEgzP8GyU8m+6MvwDQFvB63b7mFI161hyaKanoFRFpizgi0iziCyzXz4JnCciHSISAq4CHp6qays1HfqSAzM9BaUqUq0unc8APwZW2PX3NwBXAncAGGMGgQ8CjwKbgQeNMY9X49pKTaeRTDL/dV9ycAZnolTlqlLDN8Z8AvhE0eHHgdsLxqwH1lfjekrNlMIg35vQgK9qiz5pq1QFhtKJ/NfDmUSJkUrNPhrwlarAcEHAH0qPzOBMlKqcBnylKlAY8IfTyRIjlZp9NOArVYGhgjLOkJZ0VI3RgK9UBQoz/JG0BnxVWzTgK1WBsTV8DfiqtmjAV6oCw3YffiwY0S4dVXM04CtVAaczp62uWTN8VXM04CtVgeFMEkGYF2vKZ/tK1QoN+EpVYCidoD4coz5Sx7D24asaowFfqQqMpBPEwzHioaj24auaowFfqQoMZRLEQzHqw3Xah69qjgZ8pSownE4QD0ftDF8DvqotGvCVqsBIJkUsGKEuFB2zVLJStUADvlIVSOfSRINhYqEIyWx6pqejVEWqtQHKOhHZKCJbReQmjzH/0x6zTUTuFpGGalxbqemUzGaIBMNEg2ESmRTG6LbLqnZMOuCLSD3wDeASYA1wuYicUzRmAfBJ4AJjzCrgENYOWErVlFQ2bQf8CAZDOpeZ6Skp5Vs1Mvy1wLPGmAPGmAxwN7CuaEwEqAecrP4AkKrCtZWaVlbADxENhQG0rKNqSjW2OOzEytgdPcBJhQOMMbtF5CvASyJyN7AA+EO3byYiNwA3ACxbtsxtiFIzJp3LEAmEiQWjACSzKRqJz/CslPKnWjdtc0WvI4UvRKQZ+H3gAuBnwInA77l9I2PMrcaYbmNMd3t7e5Wmp1R1JO0MPxa0MvxERn9RVbWjGhn+AaCt4HW7fazQpcBLxpiXsLL8QeADwL1VuL5S0yZfww9ZOU1CSzqqhlQjw38SOE9EOkQkBFwFPCwizSLi1GReBS4SkXn2625gSxWurdS0SuUyRAIhokGnhq8Zvqodkw74xphBrI6bR4HNwIPGmMeBK4E77DHPAl8DfiMiLwGnAp+e7LWVmm5puy0zFrQy/GRGM3xVO6pR0sEYsx5YX3TsduD2gte3ALdU43pKzQRjjF3DHy3paIavaok+aauUT5lcFoPJP3gFetNW1RYN+Er55DxkZdXwnQxfSzqqdmjAV8onJ7hHgiFiWtJRNUgDvlI+pfIBPzzah68ZvqohGvCV8illl3SigXC+pKM1fFVLNOAr5ZOT4YeD2oevapMGfKV8SmXtm7Zj2jK1pKNqhwZ8pXxK5ewafiBMJBBCEJJa0lE1RAO+Uj6lCrp0RIRoMKwZvqopGvCV8qmwpGP9M5S/katULdCAr5RP+Qw/YK1IEg6E8seUqgUa8JXyyanhR/MZflgDvqopGvCV8skp6YSDVoavJR1VazTgK+VT4ZO2YHXraIavaokGfKV8KnzSFqxMP60ZvqohGvCV8ilZ8KQtWLV8zfBVLalKwBeRdSKyUUS2ishNHmPCInKziLwiIrtFpLUa11ZquqSL2jKtLh3N8FXtmHTAF5F64BvAJcAa4HIROcdl6L8BQ8BJwDKgb7LXVmo65bt0AgVdOjnN8FXtqMYWh2uBZ40xBwBE5G5gHfCsM0BEFgKvA043xpgqXFOpaVf4pC1Y/fiDqeGZnJJSFalGSacTOFTwugdYWDTmNMAAj9hln2/bvxmMIyI3iMgGEdnQ09NThekpVR2pbIaABAgGgoCV4Se1pKNqSLVu2uaKXkeKXncA24DLgFOBg8An3b6RMeZWY0y3Maa7vb29StNTavJS2XQ+uwcr4Ke1pKNqSDUC/gGgreB1u32sUC8wZIxJGmOywP8DTqnCtZWaNqlcJl+/B71pq2pPNQL+k8B5ItIhIiHgKuBhEWkWkWX2mF8BrxeRLvv15fZ5StWMZDad79ABiOqTtqrGTDrgG2MGgQ8CjwKbgQeNMY8DVwJ32GP6gfcCPxKRzVglni9N9tpKTad0NpPvwQddPE3Vnmp06WCMWQ+sLzp2O3B7weuHgDOrcT2lZkIqlyZSUNLRxdNUrdEnbZXyKZVNEx1z01aXVlC1RQO+Uj4ls2nCBTX8SMDa8UofLVG1QgO+Uj6lc5n85icwuqZOJpedqSkpVREN+Er5ZJV0Crt0rK91eQVVKzTgK+VTMpsZ05YZtrN97cVXtUIDvlI+pXOZfJCH0VUztVNH1QoN+Er5lCp68Mr5Wjt1VK3QgK+UT8U1fOcGblIzfFUjNOAr5VMym3Yt6WiGr2qFBnylfErnMmNWy9SbtqrWaMBXyievtkwt6ahaoQFfKZ9SxW2ZdravJR1VKzTgK+VTqrgtM1/S0Qxf1QYN+Er5kMllyZmca1um1vBVrdCAr5QPThY/pi3TLuno0gqqVmjAV8oH58asW5dOWjN8VSOqEvBFZJ2IbBSRrSJyU5mxfysiG6txXaWmi3NjNhxwWTxNa/iqRkw64ItIPfAN4BJgDXC5iJzjMfZC4J2TvaZS0220pOOylo526agaUY0Mfy3wrDHmgDEmA9wNrCseJCJtwFeAP6vCNZWaVqMlHbfVMjXDV7WhGgG/EzhU8LoHWFg4QEQE+C/g74CDpb6ZiNwgIhtEZENPT08VpqfU5DmdOIVtmVFdWkHVmGrdtM0VvY4Uvf4Q8GtjzGPlvpEx5lZjTLcxpru9vb1K01NqcpygHnXJ8PVJW1UrqhHwDwBtBa/b7WOFlgPvFpEtwMPASSLyiypcW6kJ2d63lzs230/OFOcq7txKOqFAEEEqKunc9+oTPLl/c2WTVapKQuWHlPUk8J8i0gEcBa4CPiYizUCzMeY1Y8wHncEi0gX8xBhzURWurdSEvOunn2VX/0EWxFu5rOv8suNTLgFfRIgEQ75LOrsHDvG+B78IwCt/fBd14egEZq7UxE06wzfGDAIfBB4FNgMPGmMeB64E7pjs91eq2o6O9LOr37qV9Pju532dM9qWOTZHigTDJH324f9q74v5r58++JKvc5SqpqrU8I0x640xa4wxq4wxn7GP3W6Mudhl7E5jzGnVuK5SE7HpyM781y8d3eXrHLcnbcH6AZD2WdLZ1rsn//XGwzt8naNUNemTtmrOefXYPgAuXnI2O47t93WO25O21uuw7z783QOHWNmymLa6ZnbYc1BqOmnAV3PO3sFDhAMhuheeTM9IH4lMquw5bjV8sFbM9Lt42p6BQyxt7GBJQzt7BrXlWE0/Dfhqztkz0ENnw3wW1s8H4PBIX9lznBp+ZFwNP+S7S2fPYA+LG9pZ3NjOngEN+Gr6acBXc44TeNvrWgA4NFw+4Ds3Zsdl+MGwry6d4XSCo4l+O8PvYN/gYYwxE5i9UhOnAV/NOXsHeljS0EFHvBWAnpHesueULumUz/APDVvXWFA/jyWN7SSyKQ6PHKt06kpNigZ8NacYY+gZ6WNBfSvtcf8ZvndJJ+zrSdujiQEA5sWaWBifB8DB4fI/aJSqpmo8eKVUzRhIDZM1OVpjjbTFmgHo8VXScc/wwwF/D171Jq2A3xptJBO2xh/RDF9NMw34ak7pSw4C0BJtJBwM0Rpr5JDPkk5AAoQCwTHHo8Eww5lE2fOPJvoBaI01krWXcziS0ICvppcGfDWnOAG/NdoAQEu0gWPJobLnpbLpcT34AOFgyFdJp7egpOOs33NkpN/3vJWqBg34ak5xAm9rrBGwAn6/j4CfzmXG1e8BIoGwry0Ojyb6CUqApkgcgyEoAc3w1bTTgK/mlD67lt5iZ/hNkXqOpQbLnpfMpsfV78HK8P1sYt6bGKA11oiIIAjzYk0c1QxfTTPt0lFzSm9BDR+gOdqQL/OUkspmXAO+3ydtjyYGaLWvCTC/ronDmuGraaYBX80pTnBvjtbn/+mnpJPKpd1LOsGwrz78/tQQTfY1AebHmrWGr6adBnw1p/QmBmgI1+Wz9eZIPcdSQ2WfevXK8P22ZQ6lR2iKxPOv59c1aw1fTTsN+GpO6UsO5Ov3YJV0Mrls2dbKlEcNP+pztcyB1DAN4YKArzV8NQOqEvBFZJ2IbBSRrSJyk8v7MRF5SES2i8g2tzFKTYe+xCAtsdFaelPEKrOUa830KumEg/5q+IOpERoidfnX8+uaOJYa8tXho1S1TDrgi0g98A3gEmANcLmInOMy9AvGmBXAmcA1InLWZK+tVKX6koNjbp46tfxyAT+dzbj24UeCYXImRzaXLXn+QHqYhvBowG+NNQGjT+AqNR2qkeGvBZ41xhwwxmSAu4F1hQOMMQljzIP21yPAK8CCKlxbqYr0FpV0nAx/MD1c8jyvtkwn6y+V5edMjqF0gsaCGr7zQ8dPh5BS1VKNgN8JHCp43QMs9BosIguA38Ha/Nzt/RtEZIOIbOjp0TXDVXX1JQZoiY0G/Hg4Blgll1JSuQyRgPtNW+d9L0Np6/5AYUmn1Z5DX0IzfDV9qnXTNlf0OuI2SESiwPeBjxpjXFesMsbcaozpNsZ0t7e3V2l6SlkrZfYlB8dk+E6ZZWiCN22dY6VaMwdSw2OuBaPPAWiGr6ZTNQL+AaCt4HW7fWwMEYkA9wA/NcbcXoXrKlWRwpUyHfV2hj9ULsP3rOGXz/AH09b3LuzScX7o9GrAV9OoGgH/SeA8EekQkRBwFfCwiDSLyDIAEYkD64FfGGM+X4VrKlWxvqKnbKGCDD+Xdi3pODX8dIkM3ykXNUYKM3y7pKM3bdU0mnTAN8YMAh8EHgU2Aw8aYx4HrgTusIetBS4G/peIbLH/aOBX06p4pUyooIbvuVqmXdIpmeFbJZ3Cm7aNkThBCeQXc1NqOlRl8TRjzHqsDL7w2O3A7fbXjwHRalxLqYkqXjgNrAenQoFg/saql7THk7bRYPkunYHU+JKOiPhex0epatEnbdWcUbw0MliBtz4UYyhdOsNPej14lW/LLFXSsW/aFpR0nHn0JTTgq+mjAV/VrEdee4af7viN7/HFK2U66iN1ZTN8z9Uy7WOl1tMZvWk7NuC3VJjhP31gC9/f9mh+AxWlKqXr4aua9NT+zbzrp/8AwG2XfYTLus4ve07xSpmO+lAsH5TdZHJZciZX8sGrUrteDZQI+Id8bmS+f+gIV6//OOlchoHUMH982lt9nadUIc3wVU36v7/9CY2ROC3RBv77pYd8ndOXGKA+HBsXuBsidQyXyPCd9W68tjiEMhl+aphYKJIf62iJNvrO8G/feB85k2NBvJX/+O16zfLVhGjAVzUnkUnxyGvPcPWqN/KHJ/8ej+95vmwNHqxlFVqLyjkA8TIZftLe0cq1LdPng1fF2T1UVtJ5bM/zrF14Ch9Zez27+g+y8fAOX+cpVUgDvqo5zx3aRiKb4qLFZ/C7i88gncvwYs/2sucVr5TpaChTw3eC+UTX0hlOJ/IPeBVqjTUwkBouu2Jmb2KATYd3cOHiM7h46dkA/HzP8yXPUcqNBnxVc57YtwlBWLvoVM5qXwnACz2vlD2veFkFR7kundGA775aJpQu6YxkksRD4wO+c/O43J66m47swGA4d8HJdMRbWT1vGU/s21TyHKXcaMBXNeeJ/RtZ07aclmgD8+uaWdrYwQt+MvyilTIdVpdOiYBvB/OJtmUOZ5LUhcY/huK0h5ZrzXy5dw8AJ7cuBeDsjlU83/Ny2V26lCqmAV/VFGMMvz38Kud0rMofW9W6NB8US+lNDIzpwXdYGf7ESjpRH0/aDmcSrgHf73o6W3tfoyXaQEe8FYAz21fSlxxk98ChkucpVUwDvqopuwcOMZAa5tT5XfljK1uWsOPYvpKdK24rZToaInUksikyHpuYOPX5qMeettYY7wx/JJPKL+FQyO96Oi/37uak1iWICGAFfPBXxlKqkAZ8VVM2H9kJwJr5y/PHVrR0ksim2Dt42PO8wfQIWZMb99AVjK6n45XlO8E87FLS8VPDH04niJfI8Et16hhj2Nq7m1UtS/PHVs9bRiQQ0oCvKqYBX9WUTUd2IAir5y3LH1vRvBiA7X17Pc/LL6vgluHbLZNerZn5Gr5rhh+0xkyght/io4bfmxigNzHAytYl+WORYJjV805g0xFtzVSV0YCvZsxv9m3ij+79DB9+7Gtl95R1bDqygxNbOseUSJY1Wbtllqpp51fKdKnhOx00I2UyfLeAHwwECUqgZFtmIpN0Lek0ReIIUnJf272D1q5vyxrH7gi6pm05mw7v8HXjNp3N8MWn/5ur1n+c7297tOx4dfzSgK9mxAs9r3DdvZ9i05EdfH/bo7zvgS/4Cl6bj+wcU84BWBBvJRQIsqdkwB+/UqYjHray7+FM0vVcJ+C71fDB+kFQMsNPJ6gLjd8ELiABmqP1Jbc5dMpUnQ1tY46vmb+cI4l+Dgwf9TzX8YWnv82/PPt9dhzbz189egvf2/pI2XPU8UkDvpp2OZPjbx7/Om3xFh75w3/hsxe+j1/t+y3373Td5jjvWHKI3QOHxtywBSvLXtzQxp5B7z2Qe+2yiduDV3V2hu+1vIJT0nGr4YPVrulVw8/ksqRymfw1irXGSi+vsG/ICfjzxxw/rc36obepzBO3ewd6uPXFH3Pd6kv4zXX/zgWL1vCpX9/G4RHXHUbVcU4Dvpp29776BJuP7OTv117PvFgT159yKUsbO7ht470lz3spf8O2a9x7ixvaS5Z0ektl+CEnwy9X0vEI+MGwZ1vmiP1bg/NbRLFy6+nsGzxMNBhmfqx5zPFT5nUBlK3j/8fGnyAifOjcawgHQ/zTRX/GcCbJF5/+Tsnz1PGpKgFfRNaJyEYR2SoiN010jKotGw5s4U8f/BJv+O6fc+1PPsXP97xQ9pxsLsvNG+5iVetS3r7idwErQ3/n6kv59b6NJW+8OsHt1KKSDsDSxo4yJR1naWS3ko6d4U+wpBMOhjxLOs5vDW43bcG6iVwu4C+qn59vyXQ0ROroalpUMuBnclm+t/UR1i2/gMV2SWhl6xL+6JRLuWvLQ+w4tt/zXMdLR3Zx40M387t3vZ933vtpfvLqr/WBrxo26YAvIvXAN4BLgDXA5SJyTqVj1Mw4OHSU+159gv/adD/rt/+qZGujYyST5BO/+k/e8aOb+PW+jaxqXcqOY/u47t5P8ZVnvlcyIPxo+y95uW8Pf33uNQTtDheAa07+PQThx9t/6XnuxsM7aKtrZoH9AFKhpY0dHBzu9Vym2GulTBjdyHziJZ2w53o4+Qzfo6RjLaDmXcPfN3R4XP3esaatq2RJ56kDL9GXHORtJ75uzPG/POdqIsEwN2+4y/NcgPtefYIrfvh3PL77OU6ZdwK7+g/wpw9+iffc/zlfJaGjiX4e3PU039r8AN/d8jC/2beJfp8359XUqMZ6+GuBZ40xBwBE5G5gHfBshWOq5uYNd3HQXmfcyYucDEkQ1+Nj3pOxY/A47nqu1/Vk7Lix37f09UOBENFgmGgwTM4YBlLDHEsNWS17yYF8614ql6YpEqc11sSJzZ2saOlkRfNiOhvaaI01ksqm2TPQw+YjO3n20DaeObg13wVSaO3CU7j+lDdzxYoLx2W1W4++xo0PfZmtvbt5z5rLuen8d1EfriOZTfN3P/83vrzhOzRH613Xa8/ksvzzM9/llHldvPXEC8a8t6B+HucuOJn7dz7Fh869Zty5YGX4a+YvH5ftglXSAdgzcIgVLYvHvW89dDW+fg+jwbhchu/2w8I6HsqvqFlsOB/wPUo6scaS+9ruHTzM6zpPc33vtPkncu+rTzCQGh6zX67j/h1PEg2GecPSs8Yc74i38t7TruDrz/+A95915bh7ImB1UH3g4X/mtLYT+eZb/p62uhayuSy3bbyPzz91J2/6/l/x5Td8gEtPOG/MecYYfrH3Re7cfD8/2/kU2aKH4QThpNYlnN1xEmd3rOLsjlXMizUynEmyu/8guwYOsqv/ALv6D/Ba/0GG0gmyJkc8FKWtrpm2uhZaY400R+tpiTYQDUZI5dJkclkSmRSJTIqRTJKRTJJE1vo6m8sREEFECCCIBBCBAAHrmAiC9b4w+to5R+x5u/29q7bCVOkfL3zfuCW1J6sa360TKPxdugc4aQJjABCRG4AbAJYtW+Y2pKxf7/strx7bn880jf0xjr5mzPHC94pf58/1OG4do2hs/oDr9d3mMG5uLtcpVB+O0RptpDVm/VnW2EEkGGYoPULPSB8/3fEbjib6Xc8Fq+vjnI5V/MnpV3DugpNZ0tDOweFefrH3Bb6z5SH+4tF/4dNPfJPrVl/KBZ1rAPjZzqf4zpaHaI7W8+11n8iv3AhWueOf3/Dn9CeH+OSvb2N50yLeuGzsL3H/veVBdhzbz22XfYSAjP/l8i1da/mHJ+9g70APixvbx7yXyqbZ1rubi4uCl8MJ+PsGj7gG/N6E+zo6MFpf927L9H7SFqzM3yvDL1fSaYk20J8aJpPLEir4jQes8tfBoaN01ntl+FZp66UjO1m76NQx7xljeGDXU1y0+EzqXZZmvvGsd3Dn5vv54tP/ze1vGVth3d63l/c+8E8sa1rAHZd/LN/KGgwEed8Zb+OiJWfw5w9/lffc/zne0nU+V696I03Rel7seYXvbn2Ebb27aY018r7T38alJ5xHV/MiEpkkrx7bxws923nu0DYe3LWB73p0C8VCEU5oXMCypoU0ReoJiDCUTnBk5BibjuygNzlAf3Jo3A8TgFgwQiwUIRaKUheKEAtGCQYCGGPIGQMYciaHwWoeyBmDwWDsf+bscfljztf2+LEJW/UYzLiE7zOvey/uf+Mmrlo/Poo/+fE9aP7GYIy5FbgVoLu7e0LFwnt+/x8nctqslcllSWXT+XJFYyQ+Lji4OTrSz/Zj+zg4fJTexADRYJgF8XmsnreMBfXzxo1fUD+PM9pX8P4zr+QXe1/gmxvv4+vP/4CvPX8PAKFAkGtPfhN/030d7fGWcecHA0G+9qYP8Y4f3cSfPfRlfvj2z+Wzx6Mj/dy84S7OX3gqbz5hret8L1t+Pv/w5B3cv/NJ3nv6FWPe29a7h3QuM64l0+HUqPd5dOr0JQdde/BhNBh7ZfjJshl+2LMPfySTAnDtw4fRewr9ySHm1TWNee/gcC9Zk/Mu6difxcYjO8YF/M1Hd7J74BB/cfZVnte98awr+aenvsXTB7Zw3sLVgPU5vef+zxGQAHde/nHXz2z1vBO49w++yNeeu4fbNt47prvqjLYVfPWNf8HbTryQWFEralfzIn5v2bmA9QNpV/8BXuh5hcH0CLFghKWNHSxtWsCCeKtrQlDIGMNgeoRUNk04ECIcCBENhcueN9dVI+AfAAr/RrbbxyodozyEAkFCgaBn0PAyr65pXBDxQ0R4/ZKzeP2Ss+hLDhZ0xyynqWh7wGL14Tr+6y0f5W0//N+8+6f/wJ2Xf5zlzYv44CNfpT85xGd/9088fzU+sbmTk1qW8MCup8cFfOfm5Jq2E13PXeQE/CH3exB9yQFW15/g+l5AAsRCEc8afjqXQRCCHsEkEgiR8ijpjGTK3LS1A2pvcmDcf6t9Hj34jgXxVubHmlzr+D/b8RSCjCu5FHrvaW/lmxvv5S8f/Sr/7+2fJyRB3vOzz7F74BDfu+LT+Qfa3ESDYT7cfS3vP+tKNh3ewXAmyarWpSx0SSTciAhdzYvoal7ka7zb+W5lLFVaNQL+k8B/ikgHcBS4CviYiDQDzcaY17zGVOHaaoq1RBu4wKOG7KWzoY071n2M6+/7LG/5wYdpCNfRlxzki6+/0TNDd1zWtZZ/f/FHHEsOjdl79sWeV4iHYixvWuh6XjQYpq2uOR8ki3ktnOaIh2IMlWjLjAbDnj+oIsEwiaz7bwcj5Wr4JdbTcW6gL/YI+CJiPXHr0qnzs11Pce6Ck11/E3PEwzH+480f4eqffJwLv/N+ALImy7+96a/H/cbgpS4Updv+7UDNfpP+/ccYMwh8EHgU2Aw8aIx5HLgSuKPMGHWcWjN/OT/7Hzdzw+m/z6UnnMd33vop/uiUN5c979ITziOTy/Lo7rH385/c/xLdC08e09lTbHFDO/sGj4w7PrpSpntJB6yAPJz2vmnr1YMPTlumRw3fDvh1niUdez0dl4C/33noyqOGD9aN261HXxvTnbRn4BAbD7/KW7rcS2eFzlmwivuu/BJ/cNLr+R+r3sC9V36JdUU31NXxoyo1fGPMemB90bHbgdtLjVHHt454Kx/9nXdXdM45C1bRVtfMAzuf4h0rLwKsG65bju7iihXXlTy3s34+r7j08Ts3RefVlQj44Vi+/FIsadeJvUQDJQJ+/qat6y0rWmJ2hu/SqbNv8DAN4bqSZbTzFq7m3174Ic8c3Jrv5rlvx28A656IHyfPW8YXXn+jr7GqtukdDjWrBCTApSecxyO7n81nrU8f2ILB8DsLS5cZOhva2Dd0eFzH1eoQt80AABRsSURBVJHEMYBxT6sWKpXhp3MZzxu2AOFguEQNv3wfPuC6gNq+Qe8efMcFnacRCgR5fPfoHrfrt/+K09pO5MTmzpLnqrlHA76ada448XUMpIa5385UH9r1NPFQjLM6XDt58zob2hlKJ+hPDY85fnTEak+d59GlA1bJxXtphYxnSyZYffiebZmZJAEJeJ7fHLGyd7clkvcNHfGs3zsaI3HO7TiZx/c8B1grhj57aNu4h62UAg34ahZ6/ZIzOaFpAf+16X5GMknu3fEEl3WtHdfmV6wz35o59sat8zzC/LoyGX6JtsxSGX4kUHpphbpQxPOGbzAQpDlS73HTtqdk/d5x6Qnd/Pbwq2w9+hrf3HgfQQlw5crXlz1PzT0a8NWsE5AAf3zaW3nywGau/cmn6EsO8q5TLyt7Xme9taJk8dPD+YAf825RjYdiJdsyS9Xww4FQycXTvMo5DrcVM5PZNIdHjuXbTUu55uQ3EQtG+ItH/4VvbryXt6+8aNyDa0qBBnw1S71nzTou7DydDQe3cN3qSzjfR5ugE+SKe/GPJJySTomAH/bO8FPlMvwS6+F77XZVyG09nf1levALzatr4rMX/gmbDu9gQf08PnXB/yp7jpqbqrtQg1JVEgoEueuKT7F/6Gg+cy+no66FUCA4rqRzZKSfWChS8sG1+nBdyR2voiXaMiNB7/Xwrd2uygX88Rn+viGrvdRPSQfgnadcyiUndNMcbSh5v0HNbZrhq1krIAEWN7T5XrQqGAiyID5vXC/+0UR/yeweJlvDt5ZWcFsldDid9Nz8xNESa8hv0OJwylKLG/0FfLDaYDXYq1I04KvjSmfD/HHr6RxJ9Jes34PVpZPOZVy7bVK5DNGg9w3jcDCEwZDJZce9N5xx396wkFtJx/ktZZHP326U8kMDvjquLG5oH1fD7030l+zBh8Jdr8Zn+clMquSTtk5W7VbW8XPTtiXawLHkENmCHxh7B3poq2suW/9XqhIa8NVxpbO+jf2DR8gVLJ17ZKS/ZA8+FK6JP76On8qlS2f4dgeP2+Yrw5lk2UXvWqKNGMyY5wf2DPawtLGj5HlKVUoDvjqudDa0kcplOGI/bGWMoWekj7YSi4jB6Jr4bq2ZyUy6zINX3hm+04dfynx7lczCXaR2DxzKr/GvVLVowFfHFefJVOem57HUECOZZNlauJPhj7iUdBLZVOmAb2f4bq2Zfko6ztz2Dx0FrI029mqGr6aABnx1XCl+2vaA3d64sFzAL5Hhp7Ll1tKxA75HDb9cHd6ZmzPXnuE+ktk0SzTgqyrTgK+OK07funPjdr/dolkuw3daJ4fcSjrlMnz7veIMP5vLksymy/bhO5uG7LcD/u4BazdQzfBVtWnAV8eV1lgjsVAk34t/YNgqkywqsxNTPsMvKulkclmyJlfypq1T0ilu6XS2NyzXh18XitISbeCAXdLZY5ejNOCratOAr44rIsLShg529Vs7aDoZfke8teR5Xl06TtbuK8MvKukMl9nesNCi+rZ8GWqnPfcletNWVdmkAr6IdIvIcyKyTURuEXHf9FNEvicir9rj/lX8Pjqp1ASc1LqEl/v2ALDj2D6WNLSXrMFDQR9+0Zr4SR8BP+xx0za/Fn6Zkg5AV9MCdvbvB+Dl3t0sbeyoeA9jpcqZbIb/beB6Y8wqrI3J3+Ex7k5gBXAKsBL4/UleVylPJ7UuZeex/SSzaV7p28vKliVlz3GCa/GuV07AL7d4GjBu1yvnh0e5Lh2AFS1L2Nl/gHQ2w7be3ZzUurTsOUpVasIBX0SWA8PGmE32obuAdW5jjTHrjSWLtaet+07USlXBya1LyZoc2/v28krfXla0LC57jneGb9XhS2X4Mfs9Z6xjtKRTug8fYGXrYjK5LDv697O9by+rfPyQUqpSk8nwO4FDBa97KBPIRSQOvB14rMSYG0Rkg4hs6Onp8RqmlKfT21cA8KNXfslwJuEr4AcDQWLByIRq+M7GLImikk4lGf5Ke473bHucZDbNGfa/g1LVVHZ5ZBF5CHBbsu8DQK7omGcqY9ftbwPuNMZs9RpnjLkVuBWgu7t7/PKDSpWxvGkRi+rn87Xn7wHgnDJbIzrqwtFxffijNXzvLN15L1HU4TNawy8f8E+Z10UkEMrP+XWdp/uas1KVKJvhG2MuMcacVfwHOMDYHwTt9rFx7GD/f4A+Y8ynqzFxpbyICG9cejYAC+KtnDq/y9d58VBsXFvmaA2//OJpxWvpOPcD4j66dGKhCBctOROA09tOpL3MUhBKTcSEN0AxxmwXkWYRWW2M2QJcC9wP+dLNIntMECuzHwA+WI1JK1XOh869hsMjx3jnKZcSDAR9nWOtiV95hh+zA/r4Gr71w6POZ7fNR89/Nwb4y7Ov8jVeqUpNdser64G7RKQeeAD4ln18LXA70AUsBd4FbANesjsynzLGvHuS11bKU2dDG998y00VnRMPxxgpumnrp4bvvJfIFAX8tP+btgAnz1vGnZd/zPd8larUpAK+MeYp4CyX449hBXuMMTvRB7xUDXDP8Mt36XiVdJwM389NW6WmgwZipWx14di4tsxExs7wS2TpIkIsGBlX0hnJJAlIQLcdVLOGBnylbG4Zfipn37QNlP5lOBoMu5Z04qGo7z15lZpqGvCVsrl26fjI8J33x/Xh+1gaWanppAFfKVvcpQ/fyfCjgdJlGdeSTjrhax0dpaaLBnylbPFQbPyDV3aZJhoqHfDdSjp+drtSajppwFfKFg9HSeUyZHLZ/LGkveRxqT58sEo6bl06WtJRs4kGfKVszkYlhVl+MpMiIAFCZR7eigbDroun6RLHajbRgK+UrT7sbIIyeuM2mU37aquMBSMuXTqa4avZRQO+Urb8EskFrZkpnwHfvaST8LWOjlLTRQO+Ujan/FL48FW5DcwdMZeSzkgmpSUdNatowFfK5nTUFO56lcxmym6PCE5Jp2i1zLRm+Gp20YCvlM3pmZ9Ihm89eDV+tUyt4avZRAO+UrZ8l86YDD9dtiUTxnfppLMZ0rmM76WRlZoOGvCVso1m+KMBfyST9PW0bCwYyS/DAIUrZWqGr2YPDfhK2eKh8W2ZIz7LMk5JxxiTP6/weyo1G0wq4ItIt4g8JyLbROQWESn5/UTkahEZmMw1lZoqbm2ZvgN+MEzO5PJP6TrfQ9fSUbPJZDP8bwPXG2NWYe1p+w6vgSKyEvgQoGvFqlnJrS3Tb8CP2atpOr34zvfQm7ZqNplwwBeR5cCwMWaTfeguYJ3H2BjWlofvnej1lJpqoUCQSCA0pi1zJJPymeFbAT+RTdrnORuYa0lHzR6TyfA7gUMFr3uAhR5jbwH+3Rjz0iSup9SUixftepXwm+E72xzaN27zNXwt6ahZpOyetiLyENDm8tYHgFzRsXH9ayLyB4AYY75V/J7H9W4AbgBYtmyZn1OUqpq6ol2vrJJO+bbMmP1DwenFHy3paIavZo+yAd8Yc4nbcRFZwdgfBO3AAZehK4E3isgW+3Xc/voMY0yqeLAx5lbgVoDu7m5Tbn5KVVN9OMaQ3ZaZzWVJZtO+b9qCVQKC0Zu2fn5YKDVdJlzSMcZsB5pFZLV96FrgYQARids/EDDGfNEYs9IYs9oYsxqr7r/aLdgrNdMKtzl0grefgO+stOnU7gfTIwA0RuJTMU2lJmSyXTrXA3eJyMvAUcAp26zFDv5K1ZJ4OMqIneE7dXg/ZRnn5qzz28Fgygr4DeG6qZimUhNStqRTijHmKeAsl+OPAV0e5zRM5ppKTaV4KMaRRD9QGPDLl2Xq7cDuPKU7mB4hIAFty1Szij5pq1SBuoKNzEcDfvmg7XTjDDklndQwDeEYIvrYiZo9NOArVWBsDd9/wHdq+EMFGX6D1u/VLKMBX6kC8YK2zIoCfmhsSWcgNaL1ezXraMBXqkA8HGMkPTbDj/nqw48gCEN2d85QWgO+mn004CtVIB6KkcimyOayDKT8t1aKCPFwNF/SGUgN0xDRgK9mFw34ShVwbr6OZFIMpIYB/7309eG6fP3fyvC1hq9mFw34ShUo3PVqMG0HfJ+BOx6KjdbwtaSjZiEN+EoVKNz1qj81jCD5DpxyrGUZrDLQYGqYRi3pqFlGA75SBZxsfiA1nA/afnvpnXV4cibHYDqhyyqoWUcDvlIFWmONAPQmB+lPDdMYqfd9rpPhH0sOkTO5/PdSarbQgK9UgXzATwxUXJZpjjbQlxyiN2nt4tkS1YCvZhcN+EoVaI06Gf6AneH7L8u0RhvpSw7Qm7ACvmb4arbRgK9UgZaotbZfb2KAwfSI7w4d59y+5CBHRo4Boz88lJotNOArVSAcDNEYidObGKA/NVRZhm9n9Dv69495rdRsoQFfqSIt0QarpJOcWMB/tW/fmNdKzRYa8JUq0hpt5MjIMY4mBmira67oPIBXj+0jIAGatC1TzTKTCvgi0i0iz4nINhG5RUQ8v5+IfMQet0tEzpzMdZWaSvPrmni5bw8GQ1tdi+/znIz+lb69tEQbCHj/76DUjJjs38hvA9cbY1ZhbWL+DrdBIvJRYBVwmjHmBOC3k7yuUlNmUf189g0eBmB+JRm+HfAPDfeyuKFtSuam1GRMOOCLyHKsDck32YfuAta5jIsAfwL8ubNxuTEmN9HrKjXVOhva818vKfi6nAXxeQhifw8N+Gr2mUyG3wkcKnjdAyx0GbcMCAI/EJEtIrJeRBZ4fVMRuUFENojIhp6enklMT6mJObF5Uf7rrma3v9LuYqEIHXGrBLS8aVGZ0UpNv7IBX0QeEpHni//Ybxdn6m47RXQA+4CrjTGrgYeBf/W6njHmVmNMtzGmu73df3alVLWc2b4y/3WlvfTL7R8WZ3acVNU5KVUNoXIDjDGXuB0XkRVA4e+t7cABl6G9QNoYM2C//gHw3grnqdS06WpexBcuupFFDfMr3oT88xf9GT9+5Zdc3nX+FM1OqYkrG/C9GGO2i0iziKw2xmwBrgXuBxCROLDIGLMd2Aq0i8h5xpingcuBJ6swd6WmzPWnvnlC561qXcrfnHddlWejVHVMtkvneuAuEXkZOAp8yz6+Fqt049ygvQ74uohsBi4D/m6S11VKKVWhCWf4AMaYp4CzXI4/BnQVvH4O64eAUkqpGaJPhiil1ByhAV8ppeYIDfhKKTVHaMBXSqk5QgO+UkrNERrwlVJqjhBjzEzPwZOI9AC7Jnh6G3C4itOpFp1XZXReldF5VeZ4ndcJxphxa9PM6oA/GSKywRjTPdPzKKbzqozOqzI6r8rMtXlpSUcppeYIDfhKKTVHHM8B/9aZnoAHnVdldF6V0XlVZk7N67it4SullBrreM7wlVJKFdCAr5RSc0TNB3wR6RKRvUXH6kTkOyKyTUR+bW+47naur3GTmFubvY9v4Z8dHmNvF5G9BeO+Uc25TPR6ItItIs/Zn9EtIjKlf2dE5Hsi8qp9vX8Vjy2npuvzEpF1IrJRRLaKyE0THVPlOcXsrUe325+T17weE5GdBZ/Rx6ZhbmWvOQOf15lF/w++IiKPTWTuVZrPOSLyYsHr+SJyv/3f8n4Rmedxnq9xJRljavYP8CGsjdQHi45/Avi8/fVlwI89zvc1rorzfTPwQ4/3bgeumsbPztf1sHYsW2N//R3gD6Z4Xm8DBGvj+58Cb5+pzwuox3rwbyHW3hG/AM6pdMwUzCsGXGp/XQe8AJzlMu4xoHu6/k75ueZMfF4uc7gB+MpMfF7AzcARYGPBsduAP7W//lPgFo9zfY0r9aemM3xjzFeMMR0ub70JuMse8zNgrUem6HdctXwc+NwUfv+qsn/jGTbGbLIP3QWsm8prGmPWG0sW2IwVGGbKWuBZY8wBY0wGuJvx//5+xlSVMSZhjHnQ/noEeAVYMJXXrKJp/7wKiUgI+Gvgy9N1zULGmA8D5xYdfhPwXfvrUv+P+R3nqaYDfgmdWJm/ox+YP4lxkyYibwRGjLWvrxsDfFVEXhaRO0WkcSrmUeH1ij+fHqYpANv7Ir8dK+tyMx2fl59//xn7jABEZAHwO7jvE22Au+3SyS12sJtq5a45o58X8C7g58aYvS7vzcTnBTDfGNMHYIw5BniVavyO8zTrA75dq3ze5U9nmVNzRa8jkxw32TmWy+5vNMYsAU7F2h/445XOo8J5+b3epD+fCueF/VvWbcCdxpitHt+m6p+XBz///lX/jPwQkSjwfeCjTiAocrkxpgs4G1iEVcqYan6uOVOfVxBrP+0veAyZic8LrB80hbw+D7/jPE3XT7AJM8ZcMoHTDmAtPnTQft2ClUlMdNyk5igiFwJRY+316/U9EvY/0yJyD/C3lc6j0nn5uJ7z+Tja7WNTNi872P8foM8Y8+kS36Pqn5cLP//+U/IZlSMiEeAe4KfGmNvdxhR8RsMish44b6rn5eOaM/J52a7FKidtd3tzJj4vW6+INBhjBkWkGSuBmcw4T7M+w5+gh7H+4yIilwGbjDFp+/UKu1xQclyVjcvuRSQuIisKXr9ZbMDVwG+mYB6F13e9XuG87P8xmkVktX3atVif2VTNKYh1MzYF3Fj03kx8Xk8C54lIh/3r/VXAwyLSLCLLSo2Zgrnk2X9/1wO/MMZ8vuB4fl5idfJcbH8dBq5k6v9OuV5zpj8vez4B4O+BWfN5FXgEuMb+esz/YyKy2p5PyXG+TeUd6an+gxUUNgBZ+5/X2cfjwPeAl4EngJUF5+wELi43ropzPA943uX4xcDOgtc/wOpe2Ap8E4hP8Wfnej2Xea0Fnrc/o68DwSmcUxfWr/tbCv7cMZOfF1bX0CZgG/AJ+9h7gMdKjZni/3YXA8miz+nzhfPC6t75ObDD/oy+DASmeF6u15zpz8u+5tXAj4qOTfvnBXwGeBEYwYpZb8D6LecB+/N4AGgvGG+ALvtrz3F+/+jSCkopNUccryUdpZRSRTTgK6XUHKEBXyml5ggN+EopNUdowFdKqTlCA75SSs0RGvCVUmqO+P+TSKNXWC5P8AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = np.linspace(-10.0, 10.0, 1000)\n", "y = filter_(x)\n", "\n", "plt.plot(x, y, color='#1a9850');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, generate the sample data 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": 21, "metadata": {}, "outputs": [], "source": [ "np.random.seed(201)\n", "\n", "def convolve(x_sample):\n", " return scipy.ndimage.convolve(\n", " x_sample,\n", " filter_(np.linspace(-10.0, 10.0, len(x_sample))),\n", " mode=\"wrap\"\n", " )\n", "\n", "x_data = np.random.random((50, 101))\n", "\n", "y_data = np.array([convolve(x_sample) for x_sample in x_data])" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(50, 101) (50, 101)\n" ] } ], "source": [ "print(x_data.shape, y_data.shape)" ] }, { "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 `PrimitiveTransformer` with `n_state=2` is the equivalent of a\n", "non-discretized convolution in space." ] }, { "cell_type": "code", "execution_count": 33, "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", "(101, 2)\n" ] } ], "source": [ "pipeline = Pipeline(steps=[\n", " ('discretize', PrimitiveTransformer(n_state=2)),\n", " ('regressor', LocalizationRegressor())\n", "])\n", "\n", "pipeline.fit(x_data, y_data)\n", "y_predict = pipeline.predict(x_data)\n", "fcoeff = pipeline.steps[1][1].coeff\n", "\n", "print(y_predict[0, :4].compute())\n", "print(y_data[0, :4])\n", "print(fcoeff.shape)" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3ybxf3A8c9py3vv7E0SspwwA2GTsAKElrJKS6Eto7TQ0jbtr4PuQWmBlpCyW0bbsJJSZkjYJDh77+U43tvW1v3+kOzYlmRLsRzHyff9evmFn3vu9JwV9NXpdM/3lNYaIYQQxz9Df3dACCHE0SEBXwghThAS8IUQ4gQhAV8IIU4QEvCFEOIEYervDnQnKytLDx06tL+7IYQQA8qqVauqtdbZXcuP6YA/dOhQSkpK+rsbQggxoCil9oUrlykdIYQ4QUjAF0KIE4QEfCGEOEEc03P4QghxNHk8HkpLS3E6nf3dlajYbDaKioowm81R1ZeAL4QQQaWlpSQnJzN06FCUUv3dnW5prampqaG0tJRhw4ZF1SZuUzpKqalKqfXdnJ+jlNqolNqmlJofr+sKIUS8OJ1OMjMzj/lgD6CUIjMzM6ZPI3EJ+EqpB4B3Ij2eUioReBQ4HxgPzFZKTY3HtYUQIp4GQrBvE2tf4xLwtdb3AtO6qTIDWK21Ltdae4FFwJx4XFuIo2173QE+LdvY390QImZHa5VOAVDZ4bgKyAtXUSl1m1KqRClVUlVVdVQ6J0Qs/rzq39z3waP93Q0hYnY0l2X6uxxbwlXSWi/UWhdrrYuzs0PuDBai371098N8dtuT/d0NIWJ2tFbplANZHY6zg2VCDDiVm8PetS5E3BiNRiZOnNh+/OqrrxKPvGJ9FvCVUqlAqtZ6P7ACeEIplQPUAvOAH/fVtYUQYiCz2+2sXbs27o8br1U69wOLgRHB+fezgSuBZwG01s3AXcAyYDPwjtb6/XhcWwghRHTiMsLXWv8E+EmX4veBpzvUWQIsicf1hDgWOFxO7FZbf3dD9JGffPwEm2v2xPUxT8ocxv1n3NJjPYfDweTJkwEYNmwYr7zySlyuL3faChEjU7INb5OThpYmCfiiT/TVlI4EfCFiNPoHszm4cTc+1XXhmTieRDMSH2gkW6YQscq0kXDmEDAb+7snQsREAr4QMdr8jReoeegTGpob+rsrQsREAr4QR8C9rZodO3b2dzfEcaq5ublPHlcCvhAxcHvc7b83t/bNi1KIviIBX4gYNLW2tP/e7Gjtx54IETsJ+ELEoK7DvH1zS0s3NYU49kjAFyIWJoV1Yi4gUzpi4JF1+ELEwGSzkH7zNJwbyikaO7y/uyNETGSEL0QMmhwtaK8f+7RCUnLT+7s7QsREAr4QMdi4YSMV33+Tpv9u5dDBQ/3dHSFiIgFfiBi0zds3v7WDT99a3r+dESJGEvCFiEFz6+GlmE6nsx97Io5nv/zlL5k4cSJTpkzho48+4rLLLovL48qXtkLEoNlxeGWOwyEB/3g3a9askLIvfOEL3H777bS2tjJnzpyQ8zfffDM333wz1dXVzJs3r9O55cuX93jNzz77jJdeeok1a9awaNEi5s2bx89//vMj/RM6kRG+EDFodTjaf3e6JOCL+Pv000+55JJLMJlMXHzxxVRWVnLppZfG5bFlhC9EDPJHDCLlqvE0vrxJpnROAN2NyBMSEro9n5WVFdWIPhyr1dr+38LCQgoLC4/ocbqK1xaHc5RSG5VS25RS8yPU+XKwznal1CKlVFI8ri3E0ZQ1JJ+kC0eRcfupjJ99Sn93RxyHiouL+fjjjwFYvHgxZWVlVFVVxeWxex3wlVKJwKPA+cB4YLZSamqXOrnAT4HTtNajgUoCe9wKMaDU1NTgrWgmbfIgUgbn9Hd3xHFo5syZjB8/njlz5vC3v/2NZ599lquuuorW1t7nborHlM4MYLXWuhxAKbUImAOs7lDHAiQCSUATUA64EWKAWfbSG1T++V1SvnshB9y74Jz+7pE4Hj3wwAOdjq+//vq4PG48An4BgRF7mypgVMcKWusDSqkHgS3BN4Rc4AvhHkwpdRtwG8DgwYPj0D0h4qdt3r7sxRLcg0vhjn7ukBAxiNcqna6be1o6HiilUoHLgdOAt4DhwLnhHkhrvVBrXay1Ls7Ozo5T94SID4fTCQaFxW7F45IPqWJgiccIvxzI6nCcHSzr6AJgi9Z6C4FRfjOBsdHrcbi+EEeNy+lEmY2YLCa8Lk9/d0f0Aa01Sqn+7kZUtNYx1Y/HCH8FMF0plaOUMgHzgKVKqVSlVNuczG5gplIqI3hcDGyNw7WFOKpcThcGixGTxYLHLQH/eGOz2aipqYk5kPYHrTU1NTXYbLao2/R6hK+1blZK3QUsA8zAP7XW7yulbgZuBmZprVcrpR4BPlNK+YC1BOfphRhIRp07hdpcjW9dFc21son58aaoqIjS0tK4LYPsazabjaKioqjrx+XGK631EmBJl7Kngac7HD8EPBSP6wnRX7JOGkxBlouCKXZ21uzv7+6IODObzQwbNqy/u9FnJLWCEDGo2HsQX3kThaOHYBki+fDFwCIBX4gYfPLYErb+9R3qdh2i4rMd/d0dIWIiAV+IGHicbkwWM1vfLuHAEx/3d3eEiIkEfCFi4HW7MVkt2GxW/B5ff3dHiJhIwBciBl63F7PVgtVmA68fn1+Cvhg4JOALEQOf24PFamlf+9zQ0tTPPRIiepIPX4gYDLruVCYUjcJ0MJC5sLGliYzktH7ulRDRkYAvRAwSJuYxfPAYBk1KZ6lpJ2a7tb+7JETUZEpHiBjUbjxAS0U9uXm5WAan4VPH/i34QrSRgC9EDPY/uIz1iz+ivqyalg/2UF5V0d9dEiJqEvCFiJLf70d7fFhtNsp27KPh+XXs3be3v7slRNQk4AsRpWZn4Itau91OQkICAE2tLf3ZJSFiIgFfiCg1NDcCYLfZSLIHAn6LBHwxgEjAFyJKjS3NANhtdpISkgBocfR+Y2khjhZZlilElKxJNjLuPI3ps04nyZQIQEurBHwxcEjAFyJaZiO2CbkUDRnE6KQCcn5+PiedNrm/eyVE1OIypaOUmqOU2qiU2qaUmh+hjlkp9YBSaqdS6oBSSpKJiwGlsqYKx5oyWmubSE1MwZSbhMEmYyYxcPQ64CulEoFHgfOB8cBspdTUMFX/BrQAo4DBQH1vry3E0bR92zbqHlvJ/q278Hu8NL+zgy3rNvV3t4SIWjxG+DOA1Vrrcq21F1gEzOlYQSmVB5wO/Ex3EIdrC3HUNAfn65MSkjD4FY0vbWLT5+v6uVdCRC8eAb8AqOxwXAXkdakzAdDAe8Fpn+eCnwyEGDCag6t0khISSE1MBsDhcPRnl4SISbyWZfq7HFu6HOcA24GLgJOACuCn4R5IKXWbUqpEKVUyUHaOFyeGtiWYSfYkEqx2UOB0Ovu5V0JELx4BvxzI6nCcHSzrqA5o0Vq7tNY+4FVgXLgH01ov1FoXa62Ls7Oz49A9IeKjNRjwkxMSMRgMKLMRlwR8MYDEI+CvAKYrpXKUUiZgHrBUKZWqlBocrPMxcJZSamjweHawnRADxthTJ5F575kMGzoUIBjwXf3bKSFi0Os1ZVrrZqXUXcAywAz8U2v9vlLqZuBmYJbWulEpdQvwmlLKDHwK3N7bawtxNFlTE7COyiItORWA8b+5itNHFvdzr4SIXlwWEWutlwBLupQ9DTzd4fhdYFI8ridEf9i9ZSetKw5guF4BkJSVit9q7OdeCRE9yaUjRJRWLf+U+qdWYVaBIF+9dBsb35aZSTFwyG2CQkSpbUVOsj2worhi2WZasw71Z5eEiImM8IWIksvpRJmNGAyBl43JYsHrdvdzr4SIngR8IaLkcrpQ5sNz9iarGY/L0489EiI2EvCFiJLT5cLQIeCbLWa8EvDFACJz+EJEqfiG89FnHM4aYrZa8Lkl4IuBQwK+EFEypyeSoQ8H/It+cAPranb1Y4+EiI1M6QgRpe0fr6Pu873tx0lJSbiVr/86JESMJOALEaXNiz/hwH9Xtx/v/XgT+16Udfhi4JCAL0SUvC4PJou5/fjQxt3ULt3Wjz0SIjYS8IWIktftwWw9nPnbZrOhPT5kLx8xUEjAFyJKXpcHs6VzwMevcbgkRbIYGCTgCxEln8eLpdMI3w5AQ0tTf3VJiJjIskwhojR6/hxOK5zQfmy32cCgaGxpIj8zpx97JkR0JOALESV/ipmM7Mz24zk3XsWyIeUkpaf0Y6+EiJ5M6QgRpUNL1lG2fnf7sc0UmN5x+iSBmhgYJOALEaWaV9azr2Rr+3Hp1j3UPbuaPXv39GOvhIheXAK+UmqOUmqjUmqbUmp+D3W/p5TaGI/rCnG0ON0u8OvAypygpup6HJ/s51B5eT/2TIjo9TrgK6USgUeB84HxwGyl1NQIdc8AruvtNYU42hqDK3HsHQJ+YkJgI5Tm1pZ+6ZMQsYrHCH8GsFprXa619gKLgDldKymlsoAHgW/E4ZpCHFUN7QHf3l6WnJAEQIujtV/6JESs4hHwC4DKDsdVQF7HCkopBTwD3AdUdPdgSqnblFIlSqmSqqqqOHRPiN5rbG0GICEhob0sKfh7i4zwxQARry9t/V2OLV2OvwN8orVe3tMDaa0Xaq2LtdbF2dnZceqeEL2Tnp1B3gNzOOeKi9rLUhKTUQlm3D5vP/ZMiOjFI+CXA1kdjrODZR0NA25SSm0FlgKjlFIfxuHaQhwVbu3FkGghJSm5veykk04i/0+XMOmcGf3YMyGiF4+AvwKYrpTKUUqZgHnAUqVUqlJqMIDW+i6t9Rit9VjgPGCH1npmHK4txBH520tP8W5J9GOOvfv20fjKJqr2lbWXWY2BD7IuX/S7Xu0u28f8v/06+o4KEUe9Dvha62bgLmAZsBl4R2v9PnAl8GxvH1+IvnDHvK9ywfSzoq6/f98+mt/aQV15TXuZ1+mmduFKVi79KOrHueTqK/jNHT9iw87NMfVXiHiIyxy+1nqJ1nq81nq01vr+YNnTWutZYeru1VpPCHkQIY6SI0ln3BxciZNkP/ylbYLFhnN1GQd27o36cUp3BOruPnQg5j4I0Vtyp6044dS7mrFPL8Kandxz5aAWR2AlTnLi4TYpwWWZDmf06ZEtqYE3DJ9dXnri6JPkaeKEU9lciyHVhrOuDL/fj8HQc/BtdTgASE5Mai8zGU1gMuCKIeAXXnQynvIMVKo19o4L0UsyzBAnnPc+WE7LuzuxTSmgurU+qjatrYEpneTg3bVtDGYjTqcr6mtbTx9EwvRCdpfui77DQsSJBHxxwtl/qBSApItGUe1qjKrN9EvOIv+Ryxk5YmSnckt2MgZbdB+UfX4fhyrLqfjxOyx55j+xdVqIOJApHXHCOVgeWFrpa3Cydd9OTsoc2mMbt9+LMhmwmztPxUz89dXMKIhuDcL2A7s5eO9/Aaipqo6t00LEgYzwxQmnvCKQ3aP24U95839vRNVm9fufUf/COgy+zuVWoznqdfgbdx1OrVxfUxddZ4WIIwn44oRTXVWFsgY+3B44GN3yyJ0bttH6/p6QEf7epz/m0ydfj+oxtu/bBYDBbqa5tiGGHgsRHxLwxQknffwgBl8+FWOylUNlh6Jq43I6waCwWjoH/ObdlVRui+5NY9e+vQBkDM/HUd8cU5+FiAcJ+OKEkzp9CMU3Xog9I5maiujm0p1OF8psDCk3Wcx43NFtcVhaVgoKLrjlKtLmjo+pz0LEgwR8ccIpqzhEuimJlOx0Gqtqo2rjcjkxWMIHfK8rujn8jPFF5F09lTPOmYlxWh6tnujX7wsRD7JKR5xwPr/nObwXnc5ZN13CyoNbomrj8/sx2rpm/QaLzRL1fLx5dBaTimaR4DHi3l3L/uoyxuYPj6nvQvSGjPDFCaXZ0YK/xU12Tg7TTp2Ba1gCnijy2Z/5rSuZ/tANIeVpBdlYokzRsGvHTtL9dso27qH69x+wYu2qmPsvRG9IwBcnlO2lewDIy80lwWXEsf4QuytLe2zn8rqxGUNH+OfccTXDv31eVNde+eMX2fLP9xlaMAiAfVGuEBIiXiTgixPKrgOBgF+UV0DN9lJq//YZn639vMd2q198j30vhdazGs24o1iH39jajK/JRX5BPiOKhgJQWn4wts4L0Usyhy9OKHvK9gMwtHAwKcHMlzv27uqxXemqHfh9vpDydS9/wM63lsFN3bfftDtw09WgwkGMHjQCgEMVXTeGE6JvScAXJxR7bhrJc09i4rgJJBptAOw5sLfHdl6PB4stNMNla00jjp1VPbbfvGcHAMMGDSEjJQ2DzURVZc/thIgnCfjihGLOSSL54tGMGzwSu9ECRsXBgz1PrfjcHswpSSHlVpsN7fX3mGZ55/7dAIwdNgqA0befx7BxJx/hXyHEkYnLHL5Sao5SaqNSaptSan6Y8zal1LtKqV1Kqe3h6ghxNOzatxtjvYcksx2T0YQlLZGKQxU9tvO6vZgt5pBymy3wKaHZ2dpt+5ThOaReP5mpYwNBftTMyfjz7EfwFwhx5Hod8JVSicCjwPnAeGC2UmpqmKq/01qPACYBX1RKTe7ttYWI1et/e5Hy3y9HKQXA9O9ewfBrTumxncFmJiE1dIRvswameRqau0+z7E+3kn72SAZnFwBgrHCy/bP1sXZfiF6Jxwh/BrBaa12utfYCi4A5HStorZ1a63eCvzuAnUBuHK4tREwaauqwpx0O3OOmTqQ5pec9bkf/9FIu/l7oOvycwnwsIzJwerpPr7Bh3XqS61T7G83u/5aw8ZG3Yuy9EL0Tj4BfAFR2OK4C8iJVVkrlAqcCKyKcv00pVaKUKqmqki+1RHw11zWSlJHSfmwod7DtjZ6XZbp9HqzG0CmdmbPPIet7Z5GQkhim1WHv/XURB5/+tP04KzsLb5MTbxQ3fQkRL/Fah+/vchx6hwqglLIC/wF+pLUOu7ec1nqh1rpYa12cnZ0dp+4JEeCsbyE1M739uLxkJxXPrKSirvvBxd6Hl7H5rZUh5Zbgm4DT1/0Iv7m6gbSczPbj3Nxc8Gv2lMnNV+LoiUfALweyOhxnB8s6UUpZgJeAN7TWT8fhukLExO/342lsJTPrcOAdXDQYgA0dNicJp3lNKfUHKkPKt3y2nsqfLWXL1sjt/X4/ztpmsnNz2ssK8/IB2FG6O6a/QYjeiEfAXwFMV0rlKKVMwDxgqVIqVSk1GEAplQAsAT7UWv8mDtcUImZN7lZSr5/MGbPPaS8bOXQYAFuC6+TDcbld4NdYgytyOtIeL97yJhoaIm+Gfqi2Eu3yUlhY0F42KL8IgN0H98f8dwhxpHod8LXWzcBdwDJgM/CO1vp94Erg2WC1GcAs4CtKqa3BHwn84qiqdTWRcOpgpk0vbi8bOySwLn7X/sgj7abWFgDsYQJ+gj0BgGZH5GWZ63duBmBw0aD2stNmnELmvWeSObIgUjMh4i4uN15prZcQGMF3LHsaeDr4+3Ig9DZFIY6inQf34t5VQ9Ksw/8rnjzyJAD2HYg80m5sbQIOr7nvKDEY8Fu6CfiGNDsZd5zKObMOf7IYnjsY66gsnKbQdA1C9BVJniZOGB98+AHVf/iQ5oM17WX5mTkMvn82oy+ZEbGdw+3CmJNIWkZ6yLmk9oDfErF9o8GJbWIeE0aMbS9LsybhXFVGyYqeVwgJES8S8MWA5Pa4OfP62Sxd/VHUbQ4eKgNgxKBhncoHjxpGPY6I7ZIzUsm9/wLOueKikHPZWTlYT8rBlhx5Wea6DetwrCkjx57WXmY0GGn813o+Xvxu1P0/VFPJ9KvPZUcw46cQsZKALwakF955lY+ff5Mbrwu9GSqStuyUo4s6B/ymFfvY9E7okss2ruCSS2uYfPgjR44g81unM2Ly2JBzbT55fTl1f/+cBFPnKSFraiJ11dFtsQjw/T/8lJKXl3HPL34YdRshOpKALwak/773JgCV+w9F3aa6sgpDgoUke+fR+MF3N7DtjcgBf9v27VT/6SN2rQ3dDtFqCKzDd3WTE9/R2orBYgxJrpaYnkxTlNsjAjiHBnLvbNu5Peo2QnQkAV8MSJ9/GgjQabfPYG9DdEG/tqYWa2pCSLnFZsXjcEVsV11TjXt7NZ7W0DotDU1UzH+Ld//zesT2jlYHBmvo+oiUjDQcdc1R9d3n97HZXo1tagH71kdeQipEdyTgiwFHa439mnGc+bsbsY7K4v3SdVG1Gzq3mCnfuDik3Gq34XFFHqG3tAZW4LQtwewo0ZaAr9ZBY0PkkbrT4cAYJtNmRlYm7sbus2y2eWvth1Rs2ceIyWNx1zazZtemqNoJ0ZEEfDHg7G4oo9Hu5RuX3UD6Ic3fH18YVTt/QQInnRqapNVqt+F1Rk6N0BJMfZxoD01nnJwQSMTmcDojtnc5nJhsofP/c265hqz5s3B6u0/LAPCXx/5K9e8/4P9uuYe8By9hn45+7l+INhLwxYDz1KvP0/zeLqZkjMSyro5PHl1CczfLItvs/ngDVIWuxrEn2PG5Iicxaw2O8LvO/QMk2QKjfmc3AX/CTecx5e5LQspHDhmGKTuRGkfP8/grl39C2sh8rj71YpKTk1l5aHOPbYToSgK+GHBeefElWt7YztjsoVwx5zK0y8tTS17sto3b42b/w8vZsyw0B/3Ft3+RQb+cHbGtwWbGNCiVjLTQdfgGgwFlMuB2Rf4OwJybRN7owSHlvhoHzW/vYPPubd32fefBvdTvKGP6rNMxGYxkbnbxzM8f6baNEOFIwBcDzu61WymcOAKDwcBtV90ERsW/l7zcbZs9hw6AhpwwGVjT09JwWzVah8+LP+7USeT86BxGjRoV9nzy9MGkD4m8vcOejzbQsCl0G0V3bTONL29izYbuv4N4bNEzoOG6udcAkNZiouy9jeyr6HlrRiE6koAvBpSNu7firGhg+qmBO2PzM3PIPmkIaz+MvKwSYHfZvkD9nNCtGg5t2kPjK5toaGkK29blDXyhG24dPsCwb8xi7IXTI1572/MfseuNVSHlQ/IDuXXabgiL5I2338SYaOG6C64EYPY5F4KG5998qdt2QnQlAV8MKM+/GRjJX3b+4U3Vzjh3Ji3VDWwri3wH6p6yQK6cwrzQZGUV2w/Q/NYOqupqQs4BfLDkHap+sxxXS/i7ca1Gc7fr8H0uN7YwX/iOKBwCwKGqyHvq+rUfw5UjufrBu7CYA284X7rwSjAo3l62NGI7IcKRgC/6zUfrVzL5ipls2hP9jUQrN61BWU1cPetwwP+/e+eT94c5rG2MnPGytDwwih5SMCjkXFJiYKVNbVNd2LZVhyrw7KsnwRKaPA1g809f443f/CPitX0uLzZ7aNvBuYVgUHS3s9vW2v3U+Jq5+qzDf292WiZpI/LY8PmaiO268vi8nHXzZSx8NXI/xfFPAr7oF36/n7lfmse6xR9x/S/vxOePLmuk5byhXPHPezutmJlcNBqrycyuhshz2kWTR5F575lMOenkkHPJbQG/MXxOe2fwC9m2JZghvH6cLZHX0/tdXhISQ9fwm4wmTMlWaqvDf7IAeH3ZWzS8tJFBxsxO5SefWUyr0RvVkk6AX6/8B6v3b+b267/GweqQ/YnECUICvugX3/j196jZfIAJ82ZSNSORP5Z0v8oGoMXjYGP1bk4b0jloG5QB58vbWLzghYhtXVaNdVQWRRmhX662Bfz65sbwbZ1OMChslvAZvo0WEx5X+MDr9XnRHh/2hNCAD3DK769n8q2hSdnarPxsBS3v7KQgOatT+Q/+bz7pd5zChupdEdu2efbj11i4fjFnXX4BPoeby2+9tsc24vgkAV8cdRt3b+WJ3zxC1vhBrHnxPb409nweeP1J/vDCo922e+a//6biwQ/Jd4QGT8/+enZ9Hvnu09WffA7rqtv3oO0oNTkVgPrG8OvhXS4XymyM+NgmsxmvO/wcvsvnIftn53HeF0PX4QPk5efToCN/OjhQWoqymhiUnd+pvDg3kKytpLz7rRmXrv6Ir1xwDamf1vPqXQ8zY975rH7tfR5f/Fy37cTxKS4BXyk1Rym1USm1TSk1/0jriIHF7XHzzV9/j2t+dTvbaqPfqu++f/0RjIrnnvonJqOJ+0+/Bfe/tvKDW7/Nxt2RA9j/lr6Fe2sVM0dPCzmXkZtFU1XkG5g+e2UptS+HrsEHOPXM08h/5HLGTJsQ9nxiThqJYyIvuzRZIwd8h8+FOS85sGl5GI0l+1n7/HsRH7uyrBxbZlJI4rXshDQcf1/LH+7+WcS2bo+bK6++CmVQPH3fn7AYzSx+7HksGUncfftdNLZGl8enoqWW33/+PH9Z+S9qGsJ/zyEGhl4HfKVUIvAocD4wHpitlJoaax3RP/x+P59vXcfHGz7Hr/1Rt3vs1X+QNWYQC370R5bXbODc/9zNvMfv5T/vLem23bv7SliXUcNv33yCC6efBUCC2cYzf38Sf4ub3z7+l4ht132+GntBOkPzQr94zSvMx1XbjN8f/m9oqmvAlhp+Dj7ZnogyGXD4wk/LTLpyJqO/G3naZdD0sWROHRb2XHl1Fc3v7qQuzAboADXr91H6+tqIj11XUUNKdugNXwBZqRmUbtgZ8W/+x5sv0bS3ilt/dDczTpoCQG56Nj974FeoMek8uqb7exc+37qOaVfO4uTvXcrDa17iJw/8ivzhg7j3wZ9GvGZXWmsqWmpZvXcTXl/ku5nF0RGPLQ5nAKu11uUASqlFwBxgdYx14uaN7Z+w9H9vtx8rpQAYNm4URSOG4GxpZdX7n6G17vQzfMIY8ocW0tzQxJoPVgCBjSoMBgNGo5HRJ48jtzCfxroG1q9Yjdfrxefz4fP78Pp8jCueSGZBDnWVNWz6bC0moxGrxYrdasNisTBu8gTSMtKpqqhk85oNtDhaaXW04nC6cDhbmXLh6aRnZVK9v5w9a7eSlpJKTmY22RlZ2CwWRo0dg91mZ+vObZSsXsXB8jIqq6qoqa6hoa6Oc+6cR15OLrVbD+KvauXMaadwfvFZZKYeDhgOj4vn3nqJ1956nXUlqzm0ZR/eBgemwhSG3T+HMemDSd3u5Ntzb6F4TOgXnBt3b+WqW69nx3urMWckcs+ffsL8r9/Dk5te5wHiLPwAACAASURBVFd3/5iXVjzI3/7zFN+86sshbWub6vnmH3/AyDPGcvcpX+h07oqZF2HPT2Ppm+/Ar0P/Tf1+P4e27GXMGaG5cAAGFQ0Cr59tB3YxbkjoDVItdY2k5WWFaQmuxlbqX1jH2tRVXDL8tJDzTp8Hqyl0KqjNlHlns732QNhz+w7so3HRRg7N3A9hbubNzM7C1+LG6XaF/Y7A2dJKwbDQNziA6TOms+O91azavoHpYyeFnP/HoufBoPjBV7/VqfyHX/4W+ws9LNzyOtdOvJAhKaH3JqzcvIYzTjsdb4uLs756Oc988RHezl3K90ru5U/33M8zf3+SJx57nCtmhr4RHqgs457f/5iS9WtIuv5kap2N1PzlYzx768keM4iTpkzkvJmz+NJFVzK8YEh7u1ang+VrPuHD1Z+xadsW0oblkjd1BC0eJysX/Jf0jAxysrLJz82jKL+ASeNPJicvF5fLyeYNm6iqq6Wytoq6+jocLidDJo4if9Rg3M0OVr/9CXabnQSbnQS7HbvNzshxo8nKz6GxoZG1n5bg8rhxupx4vT68Pi/jpk8kt6iAxpp6Nn6yGpPJhNFoxGQyYVAGxk6dQHp2JtXllWxetR6/3x/40X601kw6o5jUzHSqSsvZvm4LikAsavuZfEYxiakplO09wK6N2zrd+Dd+xiRuOvVyjIbIU4lHpGvQi/UHuB54rMPxdcAjsdbpcO42oAQoGTx4sD4Sp//9axoI+Um+8iRdsGCuzvnFBWHPp157si5YMFdn//icsOfTbp6qCxbM1ZnfnRn2fPrXZ+iCBXN1xp2nhT2feffpumDBXJ1+6/Sw57PuO0sXLJirU6+fHPZ89s/O0wUL5uqUq8d3KjfYTNqanawnPXajHvH4F3XShaMOn1dog92sTSk2PfXZr+rBC6/WtmmFGtD2/DQ95oJiff0Pv6F/8MTv9I8+XKgv//f3tLKZNAalR583Vf/qqT/rF999VW/cv10/ueF1nTpzhMZk0OffcqWurq/t9LzvOrhX2wvStDHZqj/dWBLy73LhrVdrQC9c8lz4f7cvXaQxGXRZdUXIuWWrP9aA/spP7w7b9udP/kmbClP0fz97N+x5c3qCHj/71LDn1uzYqAF980/vCnv+pCtO15lTh4Y9p7XWdy59UJ/63G1hz/3zzUUa0Pc/8UDY8zfM/6YG9Mbd20LOeXxeXfTYVfrXnzwbtu1zb7+sAX3fQ/eHPZ88MldnTxwS9tzBpipdeN+5+qSrzgg5V11fq5OH5WiD3ayXfPxOp3Mut0vf+ot7tDHJqjEoPefb1+n9jRX6w60l+s//WqhPu/ZCbbCZAq+nUXn6nmUP67+vX6zv/ON8PfmKmTpleI7GoDSgLWOy9NRnv6onPXOzNmUktJe3v15njdCTn/2KnvHEV7Qh0RLyeki6ZIwuWDBX5/72orCvl5Srxwdezz85N/zr/YbJumDBXJ31/bPDv56/Vhx4vX/7jLDnM+44NfB6/saM8K/3e8/UBQvm6rSbp4Z/Pc+fFXi9f2lS6GPfdZp2eFxh/+2iAZToMPE1LpuYA10/34W7JTGaOmitFwILAYqLi8Pf696DZ6/6KTsnXh+4aHCaQmtNRnYW6RnpuN1u9p27FwBD8N3WoAxk5WaTlpaG2+2m9PIDaMDn86K1xuv3kZ2XS2paKi3NLZReth+zyYzZGHjXN5tM5OXlkZiYSGtrK4duOYTH68XpduFwOmh1ORk6chhJKck0nF1PzbWVJNoTSLQnkJKYTHJCElnpmSijgbqr6in/YTmV9TVU1lRRWVuN2+1m4mlTsSba8RS34LvDzfCiIYwsHEpKl+WCTde18NG6FXyw6lPWbFhHXW0dXp+XswdNJjchg6K/fJlTh05kROHQsM/fJ+O+wD2/mM/K15bxo6XfBiD1ukkknjWMs2+dy7cevILzps0MaTe8YAiLX1vMRWedxwWXXMyutVvISQuMqj9ct4J3nnqFkbMmc+ul14W97vXzruWzxct48cMlfGfuLZ3OrTqwGcvYbOZeMCds2ysuvYzH3B9gzk0OOef3+/E0OsjIzAjbNjMl8AmouTn8nHZjeS3uushfrH7w50Vs+WgNXPdYyLn6psDKn7Sk1LBt83IDo+tdZXsZP2x0p3OVrXX4tZ9BqeHn/y8/80KU2cgHH38Id3U+t7+xgqS7Z/CtUVeEbVuQlMUMdxGvvPwsv3rqz/zoK4F/Z7/fzylXnEPT3kp+9eSfufT08zu1s5gtLPzxA9z3lbu48rbrWGHcz6nPf53WlQeof3IVGBSjzp7M/fN/yrXnd7j2xMvg3sCvVfU1/GfpYrbVHcA4KB2zwUTSJZVYTGYmjhvP6VNnMGvq6eSmd0iD8dUncbpd7Dq4j11le9l38ADmtAQGjxyK1+1hXfaFZKVlkpOZTW56FmlJKaSnpZOalILT46LspjIaW5ppcbTS1NpCi7OV/EEFZGZn4WxxUHrBPuw2OzaLDavZjNFgJCcnG3tiIi0tzRy8sQyvx4vH58Xj9aC1pmjIIJKSk2m9qJnyG8sxGgwYlAFlMGBQioKiAqw2G40XN1L5jYr2kb8/OJIfNGQwNruNxgsbqP5W4F4MgwrMsucV5GMxxis8dxDuXSCWH+A84D8dju8G7o+1TrifadOmHfE7nOi9/RUH9YJXntX/t+B3+sG3n9Hv7P1c+/3+Htv97h+PaBR6zKWn6Pf2rdLVrfW6aMYYrawmvWr7hojtXB63HvfEdfru9/4Scm7+h4/pUU9cq70+b9i2FS21umDBXP3UhtdDzjU4m3XOLy7Qv333ybBtWxytGtCzv/GFsOdzpwzX6WMKIvZ72lWztDHJGvbcr576swb00//7d9jzf/7XQg3ov770VMi5F5cv1vbiQv3k0kURrz3mi2foqd+9PKT88fVLdMGCuXp3fVnEtk2tzTqhMF1bc1L09gO79TtbPtXfX/43nXzpWD332zdFbNfRmort+tlNb+pnP1usH/rX43rdzk1RtRN9iz4c4a8AnlBK5QC1wDzgx0qpVCBVa70/Up04XFv0oUE5BXx97o0xt7vvhjvYU1HKa2oDN7zxC5rf3kHjym1c9Z0vM3VU+JUwABaTmXOHFPPuvhK8Pi+mDiOclQc2MTl7VMQ5zUxbCjUPfsy/Nlq5+eHOnwJqnY2YshMZVjQkbNsEmx2Mqj0NcldelwdTmA1M2litVvye8DeONTYH8vOkp4Qf4Z979rnk//VyBk0K/d5h9fq1OEoOkmYKv4Yf4Pq7v8aTG1/H5fNg7bDk9Dff/zn2XDvDvp4fsW2SPZFfPfA7vnPtbYweNJzkK8aRPHsMt333Tv5w1u0R23U0OWcUk3OCfT8lqiaiH/U64Gutm5VSdwHLADPwT631+0qpm4GbgVmR6vT22uLY9ei9v+F3rhY2VO/mee+/2ZQ8nGd+9dce241oTmXB917imcyLueWywNRPTUMd7970CJd963q4LHw7o8EItU72bN0Zcm7jji00vbUD/+TIOeuNdgtub/illV63h4QIK3wALFYrOkLAn3BWMbm/vpBxY8Jvcp6bnIEyGqhxht70tXv/XgCmjI78Jjktdwx/ff9FPtq+ivPGnQpAadUhdr29ijOvDd3dq6tvf/FWVqz5nJraWq649DK+cP4VZCek9dhODExxmSTSWi8BlnQpexp4urs64viWYk3kjMKJnHHHxKjbzDvtYu6pdfDsoufbA/6/ly4Gr5/p48Ov0GmTlJVKTXloXprVa9fS9MomvF+PPA9/8t+u4/TBoev7AVJG55OVmRn2HIDNZgO/DrvSxmcEY0YCqQmh3y0ApFmSqH9hHcvq3+bm8Z2X8ZQeOICyGBmcUxjx2kWkUfHDt1hYW8R5vw8E/IdffBx8mhvmfSliu45e+G10O4aJgU/utBXHlEE5BeRMGELJe5+2l739QSAr5LXB9MCRpOdl01wVmg+nrDKQOG1ofugmJG3sJhsOb/hNTIbedDqn3BR5tDxq0jgSzx+B0xPafv3K1TS9vhWzDj8VZTGZcZYcZNPK0LX4lYcqsGUmh9x01dHJw8dhyUxi1cqS9rLXFr+GMdHCl+d8IWI7cWKSgC+OObMuOo/W0loeWfQUAKs/X4UtN5VRg8Lf3NQmryAPV03ozVcVlYGbnrprf+jl1ax84Z2w51xed6f58a4mnVFM6ryJ+MPE9E0r1tK0ZCuJETJtAliS7dTXhr5RuZSXjKGha+S7GjRhBAc2BnLqbCrdyfYP1zHytIkRc/+IE1cfrPsRond+eNt3eOXxF7j7q99kQ2IlZZv2MGLG+B7bjZs4ns/Xraa0tpzBWYfz3ldVVaHMRrJTI0/L1G84QKs9fJri9d//D+rsfXDOt8KeN2PE7/Tg9LihS1xvbWlBmQztuezDsacl0RQm4Od97TTOLAy9+a2raTOK2fX+Ou585Xf8t/JzrEPT+c4dd/fYTpx4ZIQvjjmTR45n3449fPMv83lnfwkJF47k0nlze2w395qryLzzNBrp/OVsbU0tphRbt1MjZpsVtzP8lI67tgXtjpy+ecXryyn/9uvs2h2audLR6kBZux9XJaWn0Frfebctr99HZWsd+YmR36TazJ51AQD/+vh1rhx9NrtKNh/R6ipx/JMRvjgm5Wfm8MhXfsLB5mpemfE+Xx4feZPxNgVJgZu8ypqrmZA1vL188u2zyanuPnWT1W4NCbpttNeP1RZ5eiTBHlg22RQmGZmj1YHRGnk6CCAjJ4uq0s456jfu2UrF75dT+8PiQGKSbtxw0dW8ed/X+PK865g9/ZzuK4sTmgR8cUwrTMrizilXR1XX5jJQMf8tnj80iAt/djhKNvpaKcjtfi7carPhdYUuy/T7/WiPD6st8hx8QnD7whZH6Cogp9OJqYeA/8Uf3sqC9a+htW7P+7Ru+ybcu2pJMYRujdiVyWjixd/9vcd6QkjAF8eNEfmD8dU72bt3b6fydf98j5MnnQzhU9IDkJSeGnbqpcUZCOJWazcjfFvkEf5p37mS0rrIe9YCZNhT8Pp9NLpbSbUGdvLaunsHAGOHj+6uqRAxkTl8cdwwm8xYMhIpP1jWqfzg6+uo2Rg+m2WbS++9gaE/DV166fC4sM8oYsjoyCt8EhMCQbrFEbrJuUt5SU1L6fba1VtKqf3bZ2zYsbm9bNe+wIbsU0b1/GW1ENGSgC+OK0nZadRWVLcf1zc3ol1eMrPDp0ZuYzdZaA2zDt9gMZL+1WKmn3dmxLZDhw4h6ZIxZOSHXmPTyx9x6L3NYVp1uIbTh3N9Odv2HL5LuO2mq3C5/4U4UhLwxXElIzeLpsrDSxx3HNwLQG525B2rALa/v47yhz/E4eq8wscV3LSju3X4Q4cOJeWycWQWhl5j/3sbKC8JTffQ0ZCCQFA/UF7aXuazG0mbUNTtyiIhYiX/N4njyoTTp2CZkNO+mcTeg/sAKMzr/kvb5opaXBsqqGnsvIXf1u1bOXTXEkre+ihiWxMGfA1OGppC8+F4Xe5uv/AFGF4QuAO4rOLwSp3cy0/m/J/d1G07IWIlAV8cVy659kqSrhlPnTOwxLK08hAoKMqLnI8GICkpkByttkvAb25pQXt8WLrZ8aq6rJKK77/JB/97N+Sc1+XBntD9SptRRYElpOWVh7/cPdRSQ35Sz2vwhYiFBHxxXMlPzEL7NfsbAsFzyLQx5P/1CmadeVa37ZISAwG/rqnzRuhtq3Tall6GbRv80rbVGfqlrd/lxdZNW4DM1HTMucm4CEwfuT1u1nz7Bfa/ta7bdkLESgK+OK407qng0J2LWfzfQGLWGkcjyqDITgq/EXib5PYRfucUB80tgaWWbUsvw2nbcczpCE2/rL0+EhIit20z+YFrGTM3sJ/uxr3b8VW3tC/RFCJeJOCL48q4IaPAr9m1bzcAyxa/ReOL60mzRs5nD5CVkYUxOxGXr/PNV83Bm6naRvHhJLcF/C5f+Lp9HvIeuozL7wy/pWNHmbYU6oI58Tfs3ALAiCHdJ4sTIlYS8MVx5aQhozAkWnjjtf/i9XnZtHIdztVl7XuFRnL6WWeQ+4sLKBw1tFN5anY6CTOHUlhQEL4hhwO+w9k54Du8bpRSJFp6vlu2dFEJy/7wIgALngncNXv65B5yKggRIwn44rhiMVv46vfvoGrDPm69/x4aauqwpvQ8pZJgDqyk6ZoTP2/4INKun8yIESMjtrVZrKTNm8iw4s67WpVVllP/7GpK14cmVevKU9NC1cb9PPvmIj779ztMmXsW502NvPZfiCPRq4CvlCpWSq1RSm1XSj2kVPhhlFLq30qp3cF6D6u2hCFC9IHHfvRHcicP55+PPEFNZRX2tO6ncwDqy2uo/tNHfLq88/LLVrcD7ddYjN1nIcm6eDx54ztPwZRVHaL1k/00lFdHaHVYemYGnkYHj1ctpXBeMa8//p8e2wgRq96O8J8DbtBajwaygUg5bP8BjADGASOBy3t5XSEiMhgMLH7hJYrmX4CzqZXk9PAbiHdkNZhwb6+m7ODBTuXvLnqdQ7e/Rk1ZZbftVa0zZHvF+uC6/JSk8NsbdpSVlYV2+9hRX8pzD/yd/MycHtsIEasjDvhKqWFAq9Z6U7DoRWBOuLpa6yU6wAdsBnrexkeIXpgxdjI/veg2MCoycrpPqwCQnhzYuLulpaVTedudt8mJ3X9K2Pubt1n2+KudyuobA0s8k6MI+OnJgTel0xyDOHtQ93v3CnGkepMtswDoOOypoodArpRKAK6gm7yFSqnbgNsABg+OvAepED25afzF7F/0EBcMKe6xblZqBgDNrV0CfnCpZZK9+yWSBosJt6vz/H9Dc2CEn5bUffI0gNuvv4XVa9fw4E3ze6wrxJHqMeArpd4Fwg2R7gD8Xcoi7uMWnLd/EviH1npbpHpa64XAQoDi4mLdU/+EiMSgDPzktJujqpsaHIV3HeG3LbVM7SFoG80mPO7OSzodbgfKZiItpecppVPHT2Pj/z7tsZ4QvdFjwNdanx+uXCk1gs5vBNlAeYS6CngMqNda//wI+ilEnzIZTViHZ2JN6byE0hlcapnUzY1XACaLCY/L3als/BnTyP/zpUyd0v1uW0IcLUc8h6+13gWkKqXa1qJdCyyFwNRN8A0BpZQReBpwA9/sVW+F6EOjfnIJE+fO7FRWNGE4qReN6TFrZbgRfqsn8GZhN3efPE2Io6W3O17dALyolEoE3gb+GSyfQSDIDwUGATcC24EtwRWZK7XWkgpQHFPsJmvIOvxBxWPJTymL0OKwcVefgTZ0Xm382dIPqX3mc/zzPNDzylAh+lyvAr7WeiUQsqRAa72cQLBHa70XucFLDAB7H1mGs3ArnPOt9rLGpkZM7p5vGxl+xkRqnZ3TI+/dvhtnyUGSJSeOOEbInrZCBLmqmqjxdR6bvPPnf7GzZBN8vfu2zsomqqs7fxJoaWkBdTi5mhD9TQK+EEEWmwWPo/OUjtvlwmjp+WWy9pl3OLhtL3z3cJmjtRVlMcmuVeKYIf8nChFksdtwO7sEfLc7qoBvtlrwub2dylpbWzFaZUwljh0S8IUIstiseJ2dl1Z63R6M5p6DtsVqwe/pHPANVjPW7J7vshXiaJHhhxBBeaMG00znEb7H5cZkiby9YRur1Ybf4+tUNu1rF5FWPyGufRSiNyTgCxF01s2X0rp3RaeywnMnYDFEEfBtVnSXgO/wurCbrHHtoxC9IQFfiKAEs41WT+cRfu7MMeQmZPTYtviiM9loq0RrTVv27xUL/ovRYoar+qS7QsRM5vCFCFr90nL2/GAxfv/hFFENlXXoVk83rQKGnTQK+7RC3P7D8/iVm/bTtL/nXPhCHC0S8IUI0m4fvupWGlqb2svW/OwlVjy6pMe2LdUNuLZU0tTa3F7mdbmx2mRKRxw7JOALEZSYEEiQVtNQ117m93ixWHsO2ps/WEXNXz6hsubwJihepwdbQs/72QpxtEjAFyIoOSlwR2xdU0N7md/jw2KNmPW7nd0aSJDWMZ++z+3BZpfEaeLYIQFfiKDkxMCa+dpOI3wf1ihG+HZ7YCTf2GFKx5SdSFa+bFUojh2ySkeIoKIhg7BNysdvPFymPT6stp5H6Qn2wHRQi6MVAL/2k3nfWVw4VZboiGOHjPCFCJo2YzoZ3zyFjMLAqNzr85JyzQROPmt6j20TbIERfnNwhO/0Bu7YTZBc+OIYIgFfiKC2m6TacuJ7/D6SzhnBmCnje2w7cepkMu48jdwhhQCUVZZT9dv32fz+qr7rsBAxkoAvRNChPQco/+7/WPbmuwA0OVvwlDbgbXH10BLyc3OxTcjFkhQY6VfX1+DZW4e32dmnfRYiFhLwhQhKsiXgb3ZTX18PQFnFIap+uYw173zSY1tXixPH2jLKyg4CUNsYeIzkJEmeJo4dvQr4SqlipdQapdR2pdRDSqluH08pdY1Sqqm7OkL0l/TkVACaWgJLK5uD/22bn+9O7aEq6hasZN3nawCobw7sftW28keIY0FvR/jPATdorUcD2cDcSBWVUiOB7wA97xcnRD/ISE4HoLklMCZpcgS+gLVHEfCTEgPbGLYGV+nUB0f4qSkS8MWx44gDvlJqGNCqtd4ULHoRmBOhro3Apua3HOn1hOhrmSlpALS0BoJ2201UbUsuu5NsCwR8hzM4Z282YB6eQW52bh/0VIgj05sRfgFQ2eG4CsiLUPchYIHWektPD6qUuk0pVaKUKqmqquqpuhBxk2RPJOHUwWQMCQTpFmcg8EczpZOcGLhL1+F0ADB4wiiy7zuLSSdP6qPeChG7Hm+8Ukq9C2SFOXUH4O9SFnIPulLqKkBprf8ZTYe01guBhQDFxcU6mjZCxMug285k2OiJAGQX5ZF20xTGTTipx3Ypwbl6Z3CE37a0027qOS2DEEdLjwFfa31+uHKl1Ag6vxFkA+Vhqo4EzlFKbQ0eJwR/P1lr7Q5TX4h+k2Cy4vAEgnZCWjIJpw+hsLCwx3bpKWlk3nsmE88/BYD3Xn2Tyr8sxTmnFVL7tMtCRO2Ip3S01ruAVKXU2GDRtcBSAKVUQvANAa3177XWI7XWY7XWYwnM+4+VYC+ORVvnv8riXz4JQFV1Fe7dtWi3r4dWYDdbsY7KIiEzJdC2ohJvWRPJ9sQ+7a8QsejtKp0bgBeVUjuAWqBt2mYGweAvxEBiNBhwOgIj/M8//Izq339AdVlFj+2UUnhXlbNrwzbg8Be/mSnpfddZIWLUq+RpWuuVwOQw5cuBoRHaJPXmmkL0JXtqEg2VtcDhJZZJUY7Sq59bxer6RPgqHCorw2A3Y7PIBiji2CF32grRwfhpk2jYXUFFXVX7Esu2FTg9MZiNuFyBL2t3rt5M7klD+qyfQhwJCfhCdDD7vAvBr3nuzZfbl1hGOw9vMJtwu9xUNdehhyZz9pyw6x2E6DcS8IXo4MbZ15B0wSjKjI3tSyxTokyPYAwG/FVV20m7YQrf+eZdfdlVIWImG6AI0UF+Zg4zv3EFe00NjD1rKumtW6MP+BYTHpebpZs/w6JMTMoZ1ce9FSI2MsIXooupGSP55ONPMGcmklQ8CKs5upunpt5zGVNuu4iF3/oNrsfXYTWa+7inQsRGRvhCdGHc3sih3y/j1UYzvkRP1O0yh+TT1NhI454Kpp9zeh/2UIgjIyN8Ibq4cc41AGxY8DYVCz6Oul3Nmj18smAJaJhz3kV91T0hjpgEfCG6GFk0lMSiDAAMlug/BO9+YzWtn+wDg+KG2Vf3VfeEOGIS8IUIY3TxBACUMfqXiNkamOtPG5FHTlq4fINC9C8J+EKEce7Z5wDgrm6Ouo0lGPAvue0LfdInIXpLAr4QYdwy70ZUgpnEosyo25gtFgypNm6ad30f9kyIIyerdIQIY1zRCE7+8RUkKVvUbTwtLvwNTmbkj+vDnglx5CTgCxHBr794D7HswPOjH/yQ1095izSr5AcUxyal9bG7qVRxcbEuKSnp724IIcSAopRapbUu7louc/hCCHGCkIAvhBAniF4FfKVUsVJqjVJqu1LqIaVUxMdTSv0gWG+fUmpSb64rhBAidr0d4T8H3KC1Hk1gE/O54SoppX4EjAYmaK2HABt6eV0hhBAxOuKAr5QaRmBD8k3BoheBOWHqWYCvAXe2bVyutfYf6XWFEEIcmd6M8AuAyg7HVUBemHqDASPwslJqq1JqiVIqN9KDKqVuU0qVKKVKqqqqetE9IYQQHfUY8JVS7yql1nb9CZ7uOlIPlzg8BygDrtFajwWWAg9Hup7WeqHWulhrXZydnR3lnyGEEKInPd54pbUOuzGnUmoE0DFDVDZQHqZqHeDRWjcFj18Gbomxn0IIIXrpiO+01VrvUkqlKqXGaq23AtcCbwIopRKAfK31LmAbkK2Umq61/hyYDayI5hqrVq2qVkrtO8IuZgHVR9i2L0m/YiP9io30KzbHa7+GhCvs1Z22SqkZwEIgEXgb+JbW2qeUmgU8rbUeGqw3BXgMSAK2Al/TWtce8YWj61tJuDvN+pv0KzbSr9hIv2JzovWrV7l0tNYrgclhypcDQzscrwFm9OZaQgghekfutBVCiBPE8RzwF/Z3ByKQfsVG+hUb6VdsTqh+HdPZMoUQQsTP8TzCF0II0YEEfCGEOEEM+ICvlBqqlDrYpcyulHohmJ3zk2Den3Bto6rXi75lBdNJdPzZE6Hu00qpgx3qPRrPvhzp9WLJiBqnfv1bKbU7eL2HlVKqN/2PQ3/mKKU2KqW2KaXmH2mdOPfJFrwDflfweYrUr+VKqb0dnqMfH4W+9XjNfni+JnV5De5USi0/kr7HqT9TlVLrOxxnKqXeDP5bvqmUyojQLqp63dJaD9gf4DsE8vk0dyn/CfCb4O8XAYsjtI+qXhz7eyHwSoRzTwPzjuJzF9X1CNw4Nz74+wvAVX3cr8sARSD/0hvAFf31fBG4v2QfgRxRJuBDYGqsdfqgXzbgn8fPYgAAA/BJREFUguDvdmAdMDlMveVA8dH6fyqaa/bH8xWmD7cBD/bH8wU8ANQAGzuUPQl8Pfj714GHIrSNql53PwN6hK+1flBrnRPm1HkEsneitX4LmBFhpBhtvXj5P+DXffj4caWizIgaT1rrJTrAB2wmfEK+o2UGsFprXa619gKLCP37o6kTV1prp9b6neDvDmAnEDEh4THmqD9fHSmlTMA9wB+P1jU70lrfC0zrUnwe8K/g7929xqKtF9GADvjd6JrJsxHI7EW9XlNKnQM4dCC9RDga+LNSaodS6h9KqeS+6EeM14s2I2rcBdNzXEFg1BXO0Xi+ovn7++05AlCBzLOnEj5diQYWBadOHgoGu77W0zX79fkCbgQ+0FofDHOuP54vgEytdT2A1roBiDRVE229iI75gK8iZOtUShX00DSaTJ6x1OttH3sa3X9Ta10EnATUBuv3Sg/9ivZ6vX5+YuwXwU9ZTwL/0Fpvi/AwcX++Iojm74/7cxQNpZQV+A/wo7ZA0MVsHUhvMgXIJzCV0deiuWZ/PV9G4D7gdxGq9MfzBYE3mo4iPR/R1ovoaL2DHTEdIVtnD8oJJB+qCB6nERhJHGm9XvVRKXUGYNWBlBORHsMZ/K9HKfUS8L1Y+xFrv6K4Xtvz0yZSRtS49SsY7B8D6rXWP+/mMeL+fIURzd/fJ89RT1RgY6GXgP9v5+xZ4oiiMPwcxWYJWJkyLGixf2Bt3SqphSyuXVp/gk0KG5u01kpKi0CwCAiKmCIKKSSwkA+C+w9sJEGDXIt7lJvNuFnM3szivA8MM3PnDOflnZnD3JnLfRdC2CqKSTz6YWY7QDO3riFyluKX0yF+TvpedLAMv5wzM3sUQjg3s2niC8y/xN3J2L/h35M94sXFzJ4B3RDCL9+f9c8FA+NGzB9v92ZWszjF9M3+U3OANnCUQUeavzBfqssfjGkza/hpHaJnuTRNEn/GXgIrfcfK8OsYaJrZY+/ePwf2LM4S+2RQTAYtt/j9uwO8DyGsJ+23uiyO5Gn59hSwSP57qjBn2X65nglgFRgbvxL2gSXf/u0ZM7OG6xkYNzQ5/0jnXohF4SNw5etlb68B28A34AMwl5zTA1p/ixuhxiZwUtDeAnrJ/hvi6IUvwCZQy+xdYb4CXfPAiXu0AUxm1FQndvc/J8vrMv0ijhrqAl+Bl972AjgYFJP52rWAiz6f1lNdxNE7h8Cpe/QKmMisqzBn2X55zjbwtq/tv/sFrAGfgJ/EmrVA7OXsuh+7wEwSH4C6b98ZN+yiqRWEEKIiPNRPOkIIIfpQwRdCiIqggi+EEBVBBV8IISqCCr4QQlQEFXwhhKgIKvhCCFERrgE2w38DXiajgwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = np.linspace(-10.0, 10.0, 101)\n", "coeff = coeff_to_real(fcoeff, (101,))\n", "plt.plot(x, filter_(x), label=r'$F$', color='#1a9850')\n", "plt.plot(x, -coeff.real[:, 0] + coeff.real[:, 1], \n", " 'k--', label=r'$\\alpha$')\n", "l = plt.legend()\n", "plt.show()" ] }, { "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.8.5" } }, "nbformat": 4, "nbformat_minor": 1 }