{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bayesian signal de-blending with Gaussian Random Fields"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Florent Leclercq,
\n",
"Imperial Centre for Inference and Cosmology, Imperial College London,
\n",
"florent.leclercq@polytechnique.org"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.linalg\n",
"import math\n",
"from matplotlib import pyplot as plt\n",
"np.random.seed(123456)\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"N=1000\n",
"epsilon=0.000000000001"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup signal covariance"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# A Lorentzian distribution\n",
"gamma=2.5\n",
"signalcovar1=np.array([1./(np.pi*gamma*(1+x*x/(gamma*gamma))) for x in range(N)])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt8FfWd//HXJ+fkSiAECHcQEFTAK8ao9daKAtoKvdgW2t3a6q7tbv1t3bb7W3tv7WXX3uzNX3ft1tbaWnVdrWi1qNit1VYlyDUgEBAhXAOEEBJy//z+OBOMIeEEkjmTnLyfj8d5nDkz38n5ZHKSd2a+M/M1d0dEROR4MqIuQERE+j6FhYiIJKWwEBGRpBQWIiKSlMJCRESSUliIiEhSCgsREUlKYSEiIkkpLEREJKl41AX0lhEjRvikSZOiLkNEpF9Zvnz5PncvStYubcJi0qRJlJaWRl2GiEi/YmZvdKedDkOJiEhSCgsREUlKYSEiIkkpLEREJCmFhYiIJKWwEBGRpBQWIiKS1IAPi8MNzXz7D6+xbX9d1KWIiPRZAz4sahua+eVftvLNJ9dFXYqISJ814MNi1JAcPvmOqSwp28OL5fuiLkdEpE8a8GEBcNOlk5kwLJfbH19Hc0tr1OWIiPQ5CgsgJzPGF66dzoY9Nfz2lW1RlyMi0ucoLAJzZ47m4inD+f4zGzlY1xh1OSIifYrCImBmfPm6GVQfaeIHz26KuhwRkT5FYdHO9DFD+NCFE7nvpTfYtKcm6nJERPoMhUUHn776dAZlxbj9iXW4e9TliIj0CaGGhZnNM7MNZlZuZrd1svxyM3vVzJrN7Pp28881s7+aWZmZrTazD4ZZZ3vDBmXxz1efxp837WPp+r2pelsRkT4ttLAwsxhwF3ANMANYZGYzOjTbBnwUuL/D/DrgI+4+E5gH/MDMhoZVa0d/c9EpTB2Zzzd+v46G5pZUva2ISJ8V5p5FCVDu7lvcvRF4AFjQvoG7b3X31UBrh/kb3X1TML0T2AskHSO2t2TGMvjSu2awdX8dv3xxa6reVkSkzwozLMYB29u9rgjmnRAzKwGygM29VFe3XHFaEbPPGMmPnyunsqYhlW8tItLnhBkW1sm8E+oxNrMxwH3Ax9z9mEurzexmMys1s9LKysqTLLNrX3jndBqaW/jukg29/rVFRPqTMMOiApjQ7vV4YGd3VzazIcDvgS+6+0udtXH3u9292N2Li4p6/yjVlKJ8PnbJZB5avp01FdW9/vVFRPqLMMNiGTDNzCabWRawEFjcnRWD9o8Cv3L3/w6xxqRuuXIqwwdl8bXHy3QqrYgMWKGFhbs3A7cAS4D1wEPuXmZmt5vZfAAzu8DMKoD3A/9pZmXB6h8ALgc+amYrg8e5YdV6PENyMvnsnNMpfaOKx1fviqIEEZHIWbr8t1xcXOylpaWhfO2WVmf+T17gQG0jz33m7eRmxUJ5HxGRVDOz5e5enKydruDuhliG8ZXrZrKrup7/fD6lJ2WJiPQJCotuKpk8jHedPYb/+NNmdhw8EnU5IiIppbA4AZ+7djru8O9PvRZ1KSIiKaWwOAHjhubyiStO5fFVO3nl9QNRlyMikjIKixP0iStOZUxBDl97vIyW1vQ4OUBEJBmFxQnKzYpx2zVnULbzEA8v3558BRGRNKCwOAnzzxlL8SmFfGfJBmrqm6IuR0QkdAqLk2CWOJV2f20jP3muPOpyRERCp7A4SWeNL+D954/nnhdf5/V9tVGXIyISKoVFD3x27ulkx2N88/froi5FRCRUCoseGDk4h1uunMqz6/fyYvm+qMsREQmNwqKHPnbJJIbmZfLw8oqoSxERCY3Cooey4zFmnzGKpev30NRyzPhMIiJpQWHRC+bOHMWh+mZe2rI/6lJEREKhsOgFl59WRG5mjCVlu6MuRUQkFAqLXpCTGeOK04p4umwPrboFiIikIYVFL5l75ij21jSwquJg1KWIiPQ6hUUvufL0UcQzjCVle6IuRUSk1ykseklBXiYXTRnO02W7SZehakVE2igsetHcmaPYsq+W8r2Hoy5FRKRXKSx60dUzRgPorCgRSTsKi140uiCHcycM5el16rcQkfSisOhlc2eOZnVFNTsPHom6FBGRXhNqWJjZPDPbYGblZnZbJ8svN7NXzazZzK7vsOwGM9sUPG4Is87eNGfmKACe1qEoEUkjoYWFmcWAu4BrgBnAIjOb0aHZNuCjwP0d1h0GfAW4ECgBvmJmhWHV2ptOLcpn6sh8nUIrImklzD2LEqDc3be4eyPwALCgfQN33+ruq4GOd+CbCzzj7gfcvQp4BpgXYq29au7MUbyy9QBVtY1RlyIi0ivCDItxwPZ2ryuCeb22rpndbGalZlZaWVl50oX2trkzR9PS6jy7XnsXIpIewgwL62Red69W69a67n63uxe7e3FRUdEJFRems8YVMLYgR2dFiUjaCDMsKoAJ7V6PB3amYN3ImRlzZo7m+Y2V1DU2R12OiEiPhRkWy4BpZjbZzLKAhcDibq67BJhjZoVBx/acYF6/MWfGKBqaW3l+Y985PCYicrJCCwt3bwZuIfFHfj3wkLuXmdntZjYfwMwuMLMK4P3Af5pZWbDuAeDrJAJnGXB7MK/fKJk8jKF5mTorSkTSQjzML+7uTwJPdpj35XbTy0gcYups3XuAe8KsL0zxWAazzxjFM+t209TSSmZM1z+KSP+lv2Ah0nCrIpIuFBYhahtu9WkdihKRfk5hEaKjw62u263hVkWkX1NYhGzOzFHsOaThVkWkf1NYhGz2GRpuVUT6P4VFyDTcqoikA4VFCrQNt7q5UsOtikj/pLBIgTeHW9WhKBHpnxQWKTC6IIdzJgzV2Nwi0m8pLFJk7sxRGm5VRPothUWKzJ2ZOBSl4VZFpD9SWKSIhlsVkf5MYZFCGm5VRPorhUUKtQ23uvS1vVGXIiJyQhQWKXTWuALGFOTorCgR6XcUFilkZsyZMUrDrYpIv6OwSLG5M0druFUR6XcUFimm4VZFpD9SWKRY23CrS9fvoamlNepyRES6RWERgbbhVl/eciDqUkREukVhEYHLphWRk5mhs6JEpN9QWEQgN0vDrYpI/xJqWJjZPDPbYGblZnZbJ8uzzezBYPnLZjYpmJ9pZvea2RozW29mnwuzzijMnTlaw62KSL8RWliYWQy4C7gGmAEsMrMZHZrdBFS5+1TgTuCOYP77gWx3Pws4H/h4W5Cki7bhVp9ep7OiRKTvC3PPogQod/ct7t4IPAAs6NBmAXBvMP0wMNvMDHBgkJnFgVygETgUYq0p1zbcqvotRKQ/CDMsxgHb272uCOZ12sbdm4FqYDiJ4KgFdgHbgO+6e9qdOjR35ii2VNZSvrcm6lJERI4rzLCwTuZ17M3tqk0J0AKMBSYDnzGzKce8gdnNZlZqZqWVlf3vimgNtyoi/UWYYVEBTGj3ejyws6s2wSGnAuAA8CHgD+7e5O57gReB4o5v4O53u3uxuxcXFRWF8C2ES8Otikh/EWZYLAOmmdlkM8sCFgKLO7RZDNwQTF8PPOfuTuLQ05WWMAi4CHgtxFojo+FWRaQ/CC0sgj6IW4AlwHrgIXcvM7PbzWx+0OznwHAzKwc+DbSdXnsXkA+sJRE6v3D31WHVGqW24VafWqu9CxHpuyzxj3z/V1xc7KWlpVGXcVLm/+QFGppa+cOtl5E4GUxEJDXMbLm7H3OYvyNdwd0HLCqZyIY9Nby6TRfoiUjfpLDoA647ZyyDsmL89pVtUZciItIphUUfkJ8dZ/6543hi9U6qjzRFXY6IyDEUFn3Eh0omUt/UymMrd0RdiojIMRQWfcRZ4ws4c9wQ7n95G+ly0oGIpA+FRR+yqGQir+2uYeV2dXSLSN+isOhD5p8zljx1dItIH5Q0LMxsvJl91sweM7NlZva8mf0/M3unmSlsetHgnEzmnzOWx1ft4lC9OrpFpO847h97M/sFcA+JW4TfASwC/hF4FpgHvGBml4dd5ECyqGQiR5paeGxlx9toiYhEJ55k+ffcfW0n89cCjwT3fJrY+2UNXGePL2DGmERH999cOFFXdItIn3DcPYu2oDCz8zsuM7Pr3L3R3cvDKm4gMjMWXTiR9bsOsbqiOupyRESA7ndw/8zMzmp7YWaLgC+GU5IsOHcsuZnq6BaRvqO7YXE9cK+ZTTezvyfRbzEnvLIGtiE5mVx3zhgWr9pJjTq6RaQP6FZYuPsWEuNR/A+J4Jjj7jpGEqJFJROpa2xh8Sp1dItI9I7bwW1ma3jrUKjDgBjwspnh7meHWdxAdu6EoZwxejC/fWUbH77wlKjLEZEBLtnZUO9KSRVyDDPjQxdO5MuPlbGmopqzxhdEXZKIDGDJDkPtd/c3unoAmFl+CuockBacO46czAzuV0e3iEQsWVg8ZmbfM7PLg7GwATCzKWZ2k5ktIXFxnoSgIDeTd509lsUrd3C4oTnqckRkAEt2ncVsYCnwcaDMzKrNbD/wa2A0cIO7Pxx+mQPXopKJ1Da28Lg6ukUkQsn6LACeAta4+/awi5FjzZo4lNNHJTq6F5XoYnkRiUbSU2c9MbjC71JQi3TCzFhUMoHVFdWs3aGzlUUkGt29KO8lM7sg1EqkS+85bzzZ8Qxd0S0ikeluWLwD+KuZbTaz1Wa2xsxWh1mYvKkgL5N3nj2Gx1bupFYd3SISge6GxTXAqcCVwHUkrr+4LtlKZjbPzDaYWbmZ3dbJ8mwzezBY/rKZTWq37Gwz+6uZlQXhlNPNWtPSh0omcrihmSdWq6NbRFKvu7f7aLuu4giJK7rbHl0ysxhwF4mgmQEsMrMZHZrdBFS5+1TgThJjZmBmcRJnXH3C3WcCbwcG9E2Szj+lkGkj87n/FZ1nICKp162wMLP5ZrYJeB34E7CVxFlSx1MClLv7FndvBB4AFnRoswC4N5h+GJhtiQEc5gCr3X0VgLvvd/eW7tSarhId3RNZtf0gZTvV0S0iqdXdw1BfBy4CNrr7ZGA28GKSdcYB7f8NrgjmddrG3ZuBamA4cBrgZrbEzF41s//bzTrT2ntnjSMrnsED2rsQkRTrblg0uft+IMPMMtz9j8C5SdbpbIi3joeuumoTBy4FPhw8v8fMZh/zBmY3m1mpmZVWVlYm/Sb6u6F5WbzzrDH8bsUO6hrV0S0iqdPdsDgY3APqz8BvzOyHQLK/VhXAhHavxwMde2ePtgn6KQqAA8H8P7n7PnevA54EZnV8A3e/292L3b24qKiom99K/7aoZCI1Dc08sXpX1KWIyABy3LAws5+Y2SUk+hbqgFuBPwCbSX421DJgmplNDsbqXggs7tBmMXBDMH098FxwEeAS4GwzywtC5ApgXfe/rfR1waRCTi0apGsuRCSlku1ZbAK+C5QB/wac6e73uvuPgsNSXQr6IG4h8Yd/PfCQu5eZ2e1mNj9o9nNguJmVA58GbgvWrQK+TyJwVgKvuvvvT+o7TDNtHd0rth1k/a5DUZcjIgOEJf6RT9LI7BQSewYLgRzgfuBBd98YbnndV1xc7KWlpVGXkRJVtY1c+K2lLCqZwNcWnBl1OSLSj5nZcncvTtbuRK6zuMPdzwM+BLyXxN6CRKBwUBbXnDWaR1bs4EjjgD6jWERSpLvXWWSa2XVm9hsS11dsBN4XamVyXItKJlJT38zv16ijW0TCl6yD+2ozu4fE2Uk3kzgr6VR3/6C76060Ebpw8jCmjFBHt4ikRrI9i88DfwWmu/t17v4bd69NQV2SRFtH9/I3qtiwuybqckQkzSUbKe8d7v4zdz+QqoKk+953/niyYrp1uYiEr7sX5UkfNGxQFnPPHM0jr1ZQ36SObhEJj8Kin1tUMoFD9c0sXqlbl4tIeBQW/dzFU4Zz9vgCvvP0BqqPDOi7uItIiBQW/ZyZ8c13n8X+ww187+kNUZcjImlKYZEGzhpfwEcunsR9L73Byu0Hoy5HRNKQwiJNfGbOaRTlZ/OFR9fQ3NIadTkikmYUFmlicE4mX7luJmU7D/Grv74RdTkikmYUFmnk2rNGc8VpRXzv6Q3srq6PuhwRSSMKizRiZnx9wZk0tzq3P1EWdTkikkYUFmlm4vA8/mn2NJ5cs5s/vrY36nJEJE0oLNLQ3182hakj8/nSY2t1C3MR6RUKizSUFc/gG+8+k4qqI/z4uU1RlyMiaUBhkaYumjKc980az93Pb2HjHt2VVkR6RmGRxj5/7RkMyo7zxUfX0p3hc0VEuqKwSGPD87P53DVn8MrWAzy8vCLqckSkH1NYpLkPFE+g+JRCvvXkeqpqG6MuR0T6KYVFmsvIML7xnjOpqW/m355aH3U5ItJPKSwGgDNGD+GmyybzUGkFr7yuQQ9F5MSFGhZmNs/MNphZuZnd1snybDN7MFj+splN6rB8opkdNrPPhlnnQPCp2dMYNzSXL/5uDY3NutGgiJyY0MLCzGLAXcA1wAxgkZnN6NDsJqDK3acCdwJ3dFh+J/BUWDUOJHlZcb42fyYb9xzm5y+8HnU5ItLPhLlnUQKUu/sWd28EHgAWdGizALg3mH4YmG1mBmBm7wa2ALrJUS+5asYo5swYxQ+XbmT7gbqoyxGRfiTMsBgHbG/3uiKY12kbd28GqoHhZjYI+Ffga8d7AzO72cxKzay0srKy1wpPZ1+dP5MMM76yuEzXXohIt4UZFtbJvI5/nbpq8zXgTnc/fLw3cPe73b3Y3YuLiopOssyBZezQXD599Wk899pelpTtjrocEeknwgyLCmBCu9fjgZ1dtTGzOFAAHAAuBL5tZluBW4HPm9ktIdY6oHz0bZOYPmYIX128jsMNzVGXIyL9QJhhsQyYZmaTzSwLWAgs7tBmMXBDMH098JwnXObuk9x9EvAD4Fvu/pMQax1Q4rEMvvmeM9lTU8+dz2yMuhwR6QdCC4ugD+IWYAmwHnjI3cvM7HYzmx80+zmJPopy4NPAMafXSjhmTSxkUclEfvHi66zdUR11OSLSx1m6dHIWFxd7aWlp1GX0K9V1TVz5vf9l/LA8HvmHtxHL6KwLSUTSmZktd/fiZO10BfcAVpCXyRffNZ1V2w/y0/8tj7ocEenD4lEXINF697njeO61Sr779EZys+LcdOnkqEsSkT5IYTHAmRnf/8A5NLe08vUn1hHPMG5426SoyxKRPkaHoYTMWAY/WnQec2aM4iuLy7jvpTeiLklE+hiFhQCJwPjJh2Zx1fSRfOl3a7n/5W1RlyQifYjCQo7Kimdw14dn8Y7Ti/j8o2t4aNn25CuJyICgsJC3yI7H+OnfnM/lpxXxr4+s1nCsIgIoLKQTOZkx7v7b87nk1BH8y8OreHSFAkNkoFNYSKdyMmP87CPFXDR5OJ95aBWLV3W8rZeIDCQKC+lSblaMn3+0mOJJw/jnB1fy+9W7oi5JRCKisJDjysuK84uPXsB5E4byTw+s4A9rFRgiA5HCQpIalB3nlzeWcM74Am65fwVPaxwMkQFHYSHdkh8ExsxxBXzy/ldZun5P1CWJSAopLKTbhuRk8qsbS5g+Zgj/8OtX+d8Ne6MuSURSRGEhJ6QgN5P7bryQaaPyufm+5fx5k8Y+FxkIFBZywgryMvn1TRdyalE+f3dvKS+W74u6JBEJmcJCTkrhoCx+83cXMnnEIG66d5k6vUXSnMJCTtqwQVn8+u8Sexg337ecTz2wggO1jVGXJSIhUFhIj4zIz+bRf7yEW6+axpNrdnH19//EE6t3ki7D9YpIgsJCeiwrnsGtV53G4//nUsYV5nLL/Sv4+H3L2XuoPurSRKSXKCyk15wxegiP/MPb+Nw1Z/CnjZVc9f0/8d+l27WXIZIGFBbSq+KxDD5+xak89anLOH30YP7l4dXc8ItlVFTVRV2aiPRAqGFhZvPMbIOZlZvZbZ0szzazB4PlL5vZpGD+1Wa23MzWBM9Xhlmn9L4pRfk8ePPF3L5gJqVbDzD3zue5769baW3VXoZIfxRaWJhZDLgLuAaYASwysxkdmt0EVLn7VOBO4I5g/j7gOnc/C7gBuC+sOiU8GRnGRy6exJJbL2fWKYV86bEyFv7sJV7fVxt1aSJygsLcsygByt19i7s3Ag8ACzq0WQDcG0w/DMw2M3P3Fe7eNoBCGZBjZtkh1iohmjAsj1/dWMK3rz+b9bsOMe8Hz3P385tp0V6GSL8RZliMA9oP4lwRzOu0jbs3A9XA8A5t3gescPeGkOqUFDAzPlA8gWc/fQWXn1bEt558jff+9C9s3FMTdWki0g1hhoV1Mq/jv5LHbWNmM0kcmvp4p29gdrOZlZpZaWWl7lHUH4waksPdf3s+P150HtsP1PHOH/2ZHz67icbm1qhLE5HjCDMsKoAJ7V6PBzqOzXm0jZnFgQLgQPB6PPAo8BF339zZG7j73e5e7O7FRUVFvVy+hMXMuO6csTzzz5dzzZljuPPZjVz27ee464/lugJcpI8KMyyWAdPMbLKZZQELgcUd2iwm0YENcD3wnLu7mQ0Ffg98zt1fDLFGidDw/Gx+tOg8fn3ThZw+egjfWbKBi/9tKbf9z2o27NbhKZG+xMK8YMrMrgV+AMSAe9z9m2Z2O1Dq7ovNLIfEmU7nkdijWOjuW8zsi8DngE3tvtwcd+9yAIXi4mIvLS0N7XuR8G3aU8Mv/rKVR16toL6plUunjuDGSyfx9tNGkpHR2RFLEekpM1vu7sVJ26XL1bUKi/RRVdvIb5dt41d/eYPdh+qZPGIQH7tkEu+bNZ5B2fGoyxNJKwoL6feaWlp5au1u7nnhdVZuP8jgnDgLL5jARy6exIRheVGXJ5IWFBaSVl7dVsUvXtzKk2t24e7MnTmaGy+dTPEphZjpEJXIyVJYSFraefAI9730Bve/vI3qI02cNa6AGy+dxDVnjiEnMxZ1eSL9jsJC0lpdYzOPrtjBPS+8zubKWvKyYlw2bQRXTR/FlWeMZHi+LvgX6Q6FhQwIra3OXzbv5w9lu1i6fi+7qusxg1kTC7lq+iiumj6SqSPzdahKpAsKCxlw3J2ynYd4dv0elq7fy5od1QCcMjyP2WeM4qoZI7lg0jAyY7ozv0gbhYUMeLuqj7B0/V6Wrt/Di5v309jcypCcOG8/fSRXzRjFFacVUZCbGXWZIpFSWIi0U9vQzJ837WPp+j0899pe9tc2Es8wSiYP46rpo7hk6gimjswnpov/ZIBRWIh0oaXVWbm9imfX7+XZdXvYtPcwAIOyYpwzYSjnTRzKuRMKOXfCUIoGq6Nc0pvCQqSbtu2vo/SNA6zYdpCV2w+yftchmoOxNiYMy+XcCYWcF4TIjLFDyI7rFF1JH90NC907QQa8icPzmDg8j/fOGg/AkcYW1u6sZsW2KlZuP0jp1gM8vipxw+SsWAYzxg4J9j6GMmtiIeMLc3W2laQ97VmIdMPu6npWbq9ixbaDrNh2kNU7DlLflBiDY0R+FtPHDGHqyPzEoyjxrGs9pD/QnoVILxpdkMO8gjHMO3MMkLhv1YbdNazYfpAV26rYuKeGB17ZzpGmlqPrFOZlBgEy+M0gGZnP2IIc7YlIv6M9C5Fe0trq7Kw+Qvnew299VB7mYF3T0XZ5WTFOLcpn2sh8Tg0CZNLwQYwrzCVfd9WVFNOehUiKZWQY4wvzGF+Yx9tPH3l0vruzv7aR8r2H2bT3MJuDEPnL5v08smLHW75GQW4m44bmMr4wl3GFuW9OD81jfGEuQ/MytVcikVBYiITMzBiRn82I/GwumjL8LcsO1Texee9htlcdYUfVESqq6thx8Aiv76vlhfJ91DW2vKV9XlaMcUPbB0ke4wpzGT0kh5GDsykanK0xPyQU+lSJRGhITibnTSzkvImFxyxzdw7WNbHjYCJEKqqOsONgIlR2HDzCyu0H33J4q01eVoyiwdkU5WczckjiuWjwm4+Rg3MoGpzN8EFZxHXrE+kmhYVIH2VmFA7KonBQFmeOK+i0zeGGZnZUHWHPoXoqaxqoPNyQeK5pYG9NPRt21/BCzT4O1Td38vVhWF4WRYOzGTYoi8K8LIbmZVKYl3jPwmD66Ly8LAbnxDXE7QClsBDpx/Kz45w+ejCnjx583Hb1TS3sO9zA3po3w6QtXPYeauBAbQPrdx2iqq6R6iNNtHZx3kuGwdC3BEgmQ/OyKMjNZEhOJoNz4gzOiTMkNzHdNq/tWXsy/ZfCQmQAyMmMHe18T6a11TlU30RVXRNVdY0crGukqrZt+q3POw7WU7bzENVHmo7pX+lMbmaMIblxBncIkcE5cfKy4gzKjjMoK5Z4zo4xqG1e+/lZiWUKntRSWIjIW2RkWLD3kMVkBnV7veaWVg43NHPoSDOH6puoqX/zuaa+iUNHEs/t5x+sa2TbgToONzRT29DcrcBpkxXPID87Tl5WIlRysmLkZmaQmxkjLytOTmaM3KzE69yseOI5M4PcrPavE22y4zFyMmPkZLZNJ54zY6azzwIKCxHpFfFYxtGQOVmtrU5dUwt1Dc3UNrZQG4RIbWMztQ0t1DU2c7ghsfxwYzN1DS1HQ+ZIU+Kx73AjR5qOcKSxhfqmlqPLTkaG8ZbwyMnMICczRnY8g+zgue11Vjwj8RxLTCdexxLTb5mXcbR9Vix2dH5mzMiOZ5AZe/ORFcsgM25kxjKIZ0QbXAoLEekzMjKM/Ox4r1+c6O40NLe+GSrtgqSusZn6plYamltoCJ7rj/Nc39RCQ3Piuaa+mcqaBhqbW2lobqWxpZXG5tbgdUuXfT8nw4w3AyRmbwZKPIMzxxXw40Xn9d6bdSLUsDCzecAPgRjwX+7+7x2WZwO/As4H9gMfdPetwbLPATcBLcA/ufuSMGsVkfRlZsFhptTeMbi5pWOABKHSLlgamltoammlsdlpamk9+mhscRqbg9fNb847ujz4Gk0tzsRhuaF/L6GFhZnFgLuAq4EKYJmZLXb3de2a3QRUuftUM1sI3AF80MxmAAuBmcBY4FkzO83dT25fUkQkAvFYBvFYBj04MtdnhHk6QQlQ7u5b3L0ReABY0KHNAuDeYPphYLYlDsotAB5w9wZ3fx0oD76eiIhEIMywGAdsb/e6IpjXaRt3bwaqgeHdXBczu9nMSs2stLKyshdLFxGR9sIMi8667Tt293TVpjvr4u7uMrLoAAAGHklEQVR3u3uxuxcXFRWdRIkiItIdYYZFBTCh3evxwM6u2phZHCgADnRzXRERSZEww2IZMM3MJptZFokO68Ud2iwGbgimrwee88QAG4uBhWaWbWaTgWnAKyHWKiIixxHa2VDu3mxmtwBLSJw6e4+7l5nZ7UCpuy8Gfg7cZ2blJPYoFgbrlpnZQ8A6oBn4pM6EEhGJjkbKExEZwLo7Up7uxCUiIkmlzZ6FmVUCb/TgS4wA9vVSOWFQfT2j+npG9fVMX67vFHdPejpp2oRFT5lZaXd2xaKi+npG9fWM6uuZvl5fd+gwlIiIJKWwEBGRpBQWb7o76gKSUH09o/p6RvX1TF+vLyn1WYiISFLasxARkaQGVFiY2Twz22Bm5WZ2WyfLs83swWD5y2Y2KYW1TTCzP5rZejMrM7NPddLm7WZWbWYrg8eXU1Vfuxq2mtma4P2PuQrSEn4UbMPVZjYrhbWd3m7brDSzQ2Z2a4c2Kd2GZnaPme01s7Xt5g0zs2fMbFPwXNjFujcEbTaZ2Q2dtQmpvu+Y2WvBz+9RMxvaxbrH/SyEWN9XzWxHu5/htV2se9zf9xDre7BdbVvNbGUX64a+/XqVuw+IB4lbjmwGpgBZwCpgRoc2/wj8RzC9EHgwhfWNAWYF04OBjZ3U93bgiYi341ZgxHGWXws8ReLOwRcBL0f4895N4hzyyLYhcDkwC1jbbt63gduC6duAOzpZbxiwJXguDKYLU1TfHCAeTN/RWX3d+SyEWN9Xgc924+d/3N/3sOrrsPx7wJej2n69+RhIexY9GYwpdO6+y91fDaZrgPV0MoZHP7AA+JUnvAQMNbMxEdQxG9js7j25ULPH3P15Evc9a6/95+xe4N2drDoXeMbdD7h7FfAMMC8V9bn7054YXwbgJRJ3fY5EF9uvO7rz+95jx6sv+NvxAeC3vf2+URhIYdGTwZhSKjj8dR7wcieLLzazVWb2lJnNTGlhCQ48bWbLzezmTpZ3a+CqFFhI17+kUW/DUe6+CxL/JAAjO2nTV7bjjST2FDuT7LMQpluCw2T3dHEYry9sv8uAPe6+qYvlUW6/EzaQwqIngzGljJnlA/8D3OruhzosfpXEYZVzgB8Dv0tlbYFL3H0WcA3wSTO7vMPyvrANs4D5wH93srgvbMPu6Avb8Qsk7vr8my6aJPsshOWnwKnAucAuEod6Oop8+wGLOP5eRVTb76QMpLDoyWBMKWFmmSSC4jfu/kjH5e5+yN0PB9NPAplmNiJV9QXvuzN43gs8yrFjo/eFgauuAV519z0dF/SFbQjsaTs0Fzzv7aRNpNsx6FB/F/BhDw6wd9SNz0Io3H2Pu7e4eyvwsy7eN+rtFwfeCzzYVZuott/JGkhh0ZPBmEIXHN/8ObDe3b/fRZvRbX0oZlZC4ue3PxX1Be85yMwGt02T6Ahd26HZYuAjwVlRFwHVbYdcUqjL/+ii3oaB9p+zG4DHOmmzBJhjZoXBYZY5wbzQmdk84F+B+e5e10Wb7nwWwqqvfR/Ye7p43+78vofpKuA1d6/obGGU2++kRd3DnsoHiTN1NpI4S+ILwbzbSfxSAOSQOHRRTmJkvikprO1SErvJq4GVweNa4BPAJ4I2twBlJM7seAl4W4q335TgvVcFdbRtw/Y1GnBXsI3XAMUprjGPxB//gnbzItuGJEJrF9BE4r/dm0j0gy0FNgXPw4K2xcB/tVv3xuCzWA58LIX1lZM43t/2OWw7Q3As8OTxPgspqu++4LO1mkQAjOlYX/D6mN/3VNQXzP9l22euXduUb7/efOgKbhERSWogHYYSEZGTpLAQEZGkFBYiIpKUwkJERJJSWIiISFIKC5EQmdmk9nckFemvFBYiIpKUwkIkRcxsipmtMLMLoq5F5EQpLERSwMxOJ3Hfr4+5+7Ko6xE5UfGoCxAZAIpI3P/pfe5eFnUxIidDexYi4asmca+lS6IuRORkac9CJHyNJEbDW2Jmh939/qgLEjlRCguRFHD3WjN7F/CMmdW6e2e3JRfps3TXWRERSUp9FiIikpTCQkREklJYiIhIUgoLERFJSmEhIiJJKSxERCQphYWIiCSlsBARkaT+P40M8AtGP+/qAAAAAElFTkSuQmCC\n",
"text/plain": [
"