TestWithGaussian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook just contains quick tests of our code snippets in a si;ple case with Gaussian correlated noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.daf.persistence import Butler\n",
    "import lsst.geom as geom\n",
    "import lsst.afw.math as afwMath\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.geom as afwGeom\n",
    "import lsst.meas.algorithms as measAlg\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "# Morgadditions\n",
    "import time\n",
    "import sys\n",
    "sys.path.append('/project/mschmitz/repos/RaulModelPrototyping')\n",
    "import helperFuncs as hf\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = 8,8\n",
    "cmap = 'PuOr' # colormap; diverging is best for looking at noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adapting Arun's local cov code\n",
    "Select how far to compute the covariance for (in both directions, for both axes):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "howFar = 3\n",
    "nbNoiseReals = 10000\n",
    "printCheckpoints = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make dummy covariance matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataSize = 5,5\n",
    "nbPix = np.prod(dataSize)\n",
    "mean = np.zeros(nbPix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fdaf0f42890>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdgAAAHLCAYAAABrm4msAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X24nXV95/v3N0+EQCQJISGQEAKJQDQQIEIVH1BQU9uKnaqFPogOXFgrPTO1MyPWOdpxxjn0nLbOdMY6TSsVa0d0nFEzRyoiiHpaRYKmPITGRAQSEwh5gIQ8772/54+90msbs7Ozv/u+s/dO3q/rWtda6173Z/3urL2yf/v7ux9+kZlIkqRmjRnuDZAk6VhkBytJUgvsYCVJaoEdrCRJLbCDlSSpBXawkiS1wA5WknRMiIjbImJTRDzSz+sREX8aEWsj4qGIuKTPa9dHxJrO7fomtscOVpJ0rPgUsPQwr/88sKBzuwn4BEBETAM+DFwOXAZ8OCKmDnVj7GAlSceEzPwWsPUwq1wDfDp7fReYEhGzgDcCd2fm1szcBtzN4TvqI2IHK0k6XpwJrOvzfH1nWX/Lh2TcUN9AkqS+JseM7GZf4++7m+cfBfb0WbQsM5cN4i3iEMvyMMuHxA5WktSobvYxn1c2/r4P85U9mblkCG+xHpjT5/lsYENn+ZUHLb9vCO0ADhFLklowZsyYxm8NWA68o3M08c8Bz2fmRuAu4A0RMbVzcNMbOsuGxApWknRMiIjP0luJTo+I9fQeGTweIDP/G3An8CZgLbALeFfnta0R8e+BBzpv9ZHMPNzBUke2PU5XJ0lq0qSYmueNvbLx913Z/aUHhzhEfFQ5RCxJUgscIpYkNSqCpvaZ/rTu5t+yTXawkqSGBRFjh3sjhp1DxJIktcAKVpLUuFaGiEcZPwFJklpgBStJalgQYf1mBytJalTvUcQe5OSfGJIktcAKVpLUMIeIwQpWkqRWWMFKkhrnPlg7WElSwyIcIgaHiCVJaoUVrCSpcV7JyQpWkqRWWMFKkhrmbDpgBStJUiusYCVJjWptwvVRxg5WktQ4T9NxiFiSpFZYwUqSGhZeyQkrWEmSWmEFK0lqnPtg7WAlSQ2LcIgYHCKWJKkVVrCSpMZ5HqwVrCRJrbCClSQ1zPlgwQ5WktSw3kslepCTf2JIktQCK1hJUsMcIgYrWEmSWmEFK0lqnPtgj3IHOy4m5AQmDTp3weI59Uaj9k+Mrp219vbvquWAnUwr5SadWPwi9+yt5YA9+8eXcifG9lJux/6TS7kxY6KUA5g0sZaLPVtKuU07Tqo1CJw2vfY974naP3IMXaUcwJ59tZ/JxAk9tQazu5YDYn/x98DYE2q5fS/UckDPpNMHnXnqqafYvHlz/T9JPyIcIoaj3MFOYBLzeeWgc9+597+U2+wZP6WUm7Dle7UGN62s5YAV3deVche+dHIpN27Xk6UcwKMbZ5Zyi8bfVcp9e+PgvzcAJ540oZQDWHx+LTdu9WdKuT+772W1BoHfetdppdwLE+aXcpPiuVIO4LEnan8QXnDWnlJuzP7nSzmAMU/fXwtOObeWW//tWg7Ye/Etg85cccUV5fY0MIeIJUmN80pOQzzIKSKWRsTqiFgbEYP/80mSpGNUuYKNiLHAx4HXA+uBByJieWauamrjJEmjUdDbRRzfhlLBXgaszczHM3MfcAdwTTObJUnS6DaUfbBnAuv6PF8PXD60zZEkjXa9l0p0H+xQOthDHdqdP7NSxE3ATQDjOXEIzUmSRosxYxs/+2fUGcqfGOuBvieozgY2HLxSZi7LzCWZuWQc9VMmJEkaTYZSwT4ALIiIecBPgGuBX2tkqyRJo1cEMdYh4nIHm5ldEXEzcBcwFrgtMx9tbMskSRrFhnShicy8E7izoW2RJB0DAvfBgldykiQ1LYavg42IpcB/pndk9S8z89aDXv8Y8NrO00nAjMyc0nmtG3i489pTmfnmoWzLUe1gL1g8p3Rd4SXTfqfc5hcfv72Um3JK7Rqdz+y6uJQDuGTauoFXOoRHf1S78P65819cyvVmi8FVz5Ril7/ijFLuyfX1yRfIHbXY7FeVcje+96WlHMCYbT8o5Z7dUbu+7xOr6xelv2hJ7Wc57sd3lHKbzryhlANg5jml2Iwdd5dyOffqUk69juQCSJn5u33W/x2g7y/t3Zm5uKntsYKVJDVr+A5y+qcLIPVuRhy4AFJ/Vxi8DvhwWxvjYV6SpNFiekSs6HO76aDXD3UBpDMP9UYRMReYB9zbZ/HEzvt+NyLeMtSNtYKVJDWqxYOcNmfmkgGaPtjPXACp41rgC5k/NWHwWZm5ISLOAe6NiIcz80fVjbWDlSQ1K2DM8AwRH9EFkDquBd7bd0FmbujcPx4R99G7f7bcwTpELEk6VvzTBZAiYgK9nejyg1eKiPOAqcB3+iybGhEndB5PB66g/323R8QKVpLUqOE6D7a/CyBFxEeAFZl5oLO9DrgjM/sOH18A/HlE9NBbfN461OlX7WAlSceMQ10AKTM/dNDzPzhE7u+BRU1uix2sJKlZEYRXcnIfrCRJbbCClSQ1a/iOIh5R7GAlSY3yYv+9/BNDkqQWWMFKkpoV4RAxVrCSJLXi6FawMY6e8VMGHatOOQfwy+dcX8o9uPlPSrnJM8eWcgA7x55fyi3a9cel3Kon3l3KAcw/Z3Ipt/Wc3y7ldm3ZW8qdOu2EUg7g6Z21n+XJJ51WynXt6inlADa/UJvq7oxZE0u5eZfuLOUAfrilNkXezmm1/8tj9nYPvFI/tj9f/N6dXvu/zJja1JMjkafpOEQsSWpYDOOE6yOJQ8SSJLXAClaS1CzPgwWsYCVJaoUVrCSpYV6LGOxgJUkNC4eIAYeIJUlqhRWsJKlZnqYDWMFKktQKK1hJUsO8FjFYwUqS1AorWElSoyK8FjHYwUqSWuBBTke5g42unUzY8r1B56acckW5zeqsOJdOf18pt/Krs0o5gJ+c9dFSbv55v17K7fphbaYQgOe2n1jKrX5kUyn3qulfK+Ue3PtLpRzAvHMHP/MTwBM/fr6UW3jeyaUcwPR/+FAtePLbS7EtJ76s1h5w/vjaz/Ifd7+qlHvou0+UcgBv+pWXlHJrflLb+zZ+Qn2v3exJ5ahaYgUrSWqWE64DHuQkSVIrrGAlSY1yPthedrCSpMaFQ8QOEUuS1AYrWElSsxwiBqxgJUlqhRWsJKlRQTBmjBWsHawkqVkBY8Y4QOonIElSC6xgJUmNCjzICaxgJUlqhRWsJKlZEe6D5Wh3sPt3waaVg449s+vicpOTZ44t5aqz4ixeurGUA1j53b8s5f72sTeXcksvW1/KAcS+7lJu2sVzS7m9J76rlLv0qS+WcgCbun6xlFu8//Zagw9sruWALz377lLuZQvnlXLTT6r/6tg25o2l3Plr/rSWu3JRKQewvWdhKTdnTm1mpL17a/+vNDJZwUqSGudpOnawkqSGBZ6mAx7kJElSK6xgJUnNCoeIwQpWkqRWWMFKkhoVBGOcD9YOVpLUMIeIAYeIJUlqhRWsJKlxnqZjBStJOoZExNKIWB0RayPilkO8/s6IeDYiVnZuN/Z57fqIWNO5XT/UbbGClSQ1KmJ4JlyPiLHAx4HXA+uBByJieWauOmjVz2XmzQdlpwEfBpYACTzYyW6rbo8drCSpccM0RHwZsDYzHweIiDuAa4CDO9hDeSNwd2Zu7WTvBpYCn61ujEPEkqRjxZnAuj7P13eWHexXIuKhiPhCRMwZZPaI2cFKkpoVEGOj8RswPSJW9Lnd9LMt/4w86Pn/Bs7OzAuBrwMHpr86kuygHNUh4p1MY0X3dYPOXTJt3cAr9dfm2PNLuZ+c9dFSrjrlHMDin/v7WptPXlbK7T9lcSkHsO7Z2v6VCbtqf9Odse+RUu7+Z19RygEsml3777Fn0e+Wcut+srOUA3jNhRNLuV27u0q5iZvuK+UAVm+9qJTbePINpdyMk04s5QCmdtWmdHx2d226y+kTt5ZyAN2cXs6OIpszc8lhXl8PzOnzfDawoe8Kmbmlz9O/AP6wT/bKg7L3VTcUrGAlSQ0Leidcb/p2BB4AFkTEvIiYAFwLLP+pbYvo+9fPm4HHOo/vAt4QEVMjYirwhs6yMg9ykiQdEzKzKyJuprdjHAvclpmPRsRHgBWZuRz4PyLizUAXsBV4Zye7NSL+Pb2dNMBHDhzwVGUHK0lq3HBdKjEz7wTuPGjZh/o8/gDwgX6ytwG3NbUtQ+pgI+IJYAfQDXQNMDYuSToORHglJ2imgn1tZm5u4H0kSTpmOEQsSWqcs+kM/SjiBL4WEQ8e4nwkACLipgPnLD333JD2F0uSNGoMtYK9IjM3RMQM4O6I+MfM/FbfFTJzGbAM4IILLhzSSbuSpJEvwgnXYYgdbGZu6Nxviogv0nsdyG8dPiVJOtY5RDyEIeKIOCkiJh94TO9JubXL7UiSdIwZSgU7E/hiRBx4n/+emV9tZKskSaOXp+kAQ+hgO9MB1S4qKknSMc7TdCRJjeq9FrH7YI9qBzvpxLFc+NLJg849+qPx5TYX7frjUm7+eb9eyv3tY28u5aA+K87iuf+plFux7c9KOYC5p+wp5cZuX13K7T/l5aXcvAUTSjmA8T07SrlxL/y4lJt/Uk8pB5Bdk0q5F7rnlnJru15WygFcOOU7teC2tbXcpsH/zvknPd2l2LT5td8fY7fVZu8B6J44gmbTcYgYcDYdSZJa4RCxJKlRAYwZ6xCxFawkSS2wgpUkNSvCfbBYwUqS1AorWElS4zxNxw5WktSwwA4WHCKWJKkVVrCSpGYFhBWsFawkSW2wgpUkNc59sHawkqSGBUFnKtPjmkPEkiS14OhWsD17GbfryUHHzp3/4nKTq554dym364d7S7mll9Vnw9h/yuJSrjorzpKpv13KAazY+l9Kub2nX13K7d5Tm9Wku7urlAPIqM3Ec9f9tdlbXnFlbWYbgBc9/NFSbtai95Ryu+PUUg6Adc+UYvte/M5S7u6v1v9PXv0L80u5MfufK+WejkWlHMC0crIF4RAxWMFKktQK98FKkhrnaTp2sJKkhoVDxIBDxJIktcIKVpLUOIeIrWAlSWqFFawkqWHhPlisYCVJaoUVrCSpWc6mA9jBSpIaFsAYr0XsELEkSW2wgpUkNc4hYitYSZJaYQUrSWqWl0oEjnIHu2f/eB7dOHPQuXNrM0YBMP+c2tRhz20/sZSLfbVp1QDWPVv7Qs49ZU8pV51yDmDJtN8p5T73j39Vyr1o8vhSridLMQDG7qpNc7Z04apS7unds0s5gIkX/5tSbsIz3yzlxs96YykHsH32taXc+id3lnKvX1r/XPfu6ynldu85qZSb9cztpRzA3qk3lrNNC8IhYhwiliSpFQ4RS5Ka5RAxYAUrSVIrrGAlSY1zH6wdrCSpYV7JqZdDxJKkY0ZELI2I1RGxNiJuOcTr74uIVRHxUETcExFz+7zWHRErO7flQ90WK1hJUrOG6WL/ETEW+DjwemA98EBELM/MvufO/QBYkpm7IuI9wP8N/Grntd2Zubip7bGClSQdKy4D1mbm45m5D7gDuKbvCpn5jczc1Xn6XaB+ovQArGAlSQ0btgnXzwTW9Xm+Hrj8MOvfAPxtn+cTI2IF0AXcmplfGsrG2MFKkhrX0hDx9E4HeMCyzFzWt9lDZA55PbeI+A1gCfCaPovPyswNEXEOcG9EPJyZP6purB2sJGm02JyZSw7z+npgTp/ns4ENB68UEVcDHwRek5l7DyzPzA2d+8cj4j7gYqDcwboPVpLUqOhcyanp2xF4AFgQEfMiYgJwLfBTRwNHxMXAnwNvzsxNfZZPjYgTOo+nA1cAtQuLd1jBSpKOCZnZFRE3A3cBY4HbMvPRiPgIsCIzlwP/D3Ay8D+i91zdpzLzzcAFwJ9HRA+9xeetBx19PGhHtYM9MbazaPxdgw+ueqbc5tZzfruUW/3IpoFXOoRpF88deKV+TNhVG1AYu311Kbf39KtLOajPivOr57+rlFv56FtLuV1z317KATy+YVYpN3PGvFLu9M2F/xsdq7a+spRbOGlsKTdx1cdLOYBd895dyj27YXspN+O0wc/gdcCTT2wt5arbOufct5VyAEOYdKwVw3Ulp8y8E7jzoGUf6vP4kL/4MvPvgUVNbotDxJIktcAhYklS47xUoh2sJKlhEU64Dg4RS5LUCitYSVLjnHDdClaSpFZYwUqSmjVMs+mMNHawkqTGeRCxQ8SSJLXCClaS1Kig91Sd450VrCRJLbCClSQ1zmOc7GAlSU3zKGLgKHewO/afzLc3Dn7Wj8tfcUa5zV1b9g680iG8avrXSrm9J9ZmiwE4Y98jpdz+U15eyu3e013KAbxo8vhSrjorzuKXfKGU+/T331DKAbx47XtKuRNe/r5S7tN312d9ue43isENO0qxnjmvLTYI03Y9UMq95sR7aw1uv7iWA6Zv+2otOKs2o1LXWe+otQfU/zerLVawkqTGeYyTBzlJktQKK1hJUqOC8DQdjqCCjYjbImJTRDzSZ9m0iLg7ItZ07qe2u5mSJI0uRzJE/Clg6UHLbgHuycwFwD2d55IkAb2n6TR9G20G7GAz81vA1oMWXwPc3nl8O/CWhrdLkjSKRUTjt9GmepDTzMzcCNC5n9HcJkmSNPq1fpBTRNwE3AQwc0b9fFZJ0igRnqYD9Qr2mYiYBdC539Tfipm5LDOXZOaSU6ZMKzYnSdLoUu1glwPXdx5fD3y5mc2RJI12B2bTOd73wQ44RBwRnwWuBKZHxHrgw8CtwOcj4gbgKeBtbW6kJGl0CS9jNHAHm5nX9fPSVQ1viyRJxwyv5CRJatyYUTik27Sj2sGOGROceNKEQeeeXL+r3Oap004o5R7c+0ul3KVPfbGUA7j/2VeUcvMWDP4zBeju7irlAHqylts19+2lXHVWnHdcclMpB/D9jR8u5XZPfkkpd+bcp0o5gPEb7inltpxa+57v399TygGMO6H2i3f6rO2l3JNjX13KAcyd/FApt++l7y3lxm9bWcoBdJ96ajmrdljBSpKa5Wk6gB2sJKlhXuy/l8d5SZLUAitYSVLjLGCtYCVJaoUVrCSpce6DtYKVJKkVVrCSpMaNxgnSm2YHK0lqVHgeLOAQsSRJrbCClSQ1LhwjtoKVJKkNVrCSpMa5D9YOVpLUsAinq4Oj3MFOmgiLzy8Ec0e5zad3ji3l5p07pZTb1PWLpRzAotm1H8f4ntrnk1Gb5g5g7K71pdzjG2aVci9e+55SrjrlHMAls/5dKfd3Oz5fyl01tZYD+Pb6Xyvlzjm5u5Q7c3dtejyAVXtfWcqtWLOglDvvouLcisD3+Oel3ITVL5RyC8+/tJTTyGQFK0lqnAWsBzlJktQKK1hJUuO8FrEVrCSpBQeu5tTk7cjajaURsToi1kbELYd4/YSI+Fzn9fsj4uw+r32gs3x1RLxxqJ+BHawk6ZgQEWOBjwM/DywErouIhQetdgOwLTPnAx8D/rCTXQhcC7wEWAr8Wef9yuxgJUkNCyKavx2By4C1mfl4Zu4D7gCuOWida4DbO4+/AFwVvW9+DXBHZu7NzB8DazvvV2YHK0kaLaZHxIo+t5sOev1MYF2f5+s7yw65TmZ2Ac8Dpx5hdlA8yEmS1KiIJKJ+/vFhbM7MJYdr+hDLDt6Q/tY5kuygWMFKko4V64E5fZ7PBjb0t05EjANOAbYeYXZQ7GAlSc1KILubvw3sAWBBRMyLiAn0HrS0/KB1lgPXdx6/Fbg3M7Oz/NrOUcbzgAXA94byMThELElqWBLZc/RbzeyKiJuBu4CxwG2Z+WhEfARYkZnLgU8Cfx0Ra+mtXK/tZB+NiM8Dq4Au4L2ZR9ar98cOVpJ0zMjMO4E7D1r2oT6P9wBv6yf7UeCjTW2LHawkqXlDK/6OCUe1g409Wxi3+jODzuXsV5XbPPmk00q5J378fCm3eP/tA6/Ujz2LfreUG/fCj0u5u+6fXMoBLF24qpSbOWNeKXfCy99Xyu2e/JJSDuqz4lwx+e2l3Pefrs3eA/CqdV8q5fadckMp1z1xUSkHsPD7tRmOznjNH5VyTz+9q5QDuHT+9lLuyRfOKOW+ee+6gVfqx6uvOrucVTusYCVJDUsYhn2wI40drCSpceEQsafpSJLUBitYSVKzMj3ICStYSZJaYQUrSWqeBznZwUqSmuYQMThELElSK6xgJUmNG45rEY80VrCSJLXAClaS1LAE3AdrBStJUgusYCVJzUo8TYej3MFu2nESf3bfywadu/G9Ly232bWr9kNeeN7JtQYf2FzLAet+srOUm39S7d/4iivnlnIAT++eXcqdvvmuUu7Td88s5c6c+1QpB3DV1NpsOtVZcS45vTbLDMDK77+xFhx7Qim2/8Q5tfaAnstr/84p2/6ulNv5ostLOYBxu9eXcuds+FwpN/fl7yjloHeG8JHD03TAIWJJklrhELEkqXGepmMFK0lSK6xgJUkNcx8s2MFKkhpnBwsOEUuS1AorWElSoyI9yAmsYCVJaoUVrCSpYe6DBTtYSVIbHCJ2iFiSpDZYwUqSGuYQMVjBSpLUCitYSVLjAvfBHtUO9rTp4/itd5026NyYbT8ot7n5hdpUd9P/4UOl3JeefXcpB/CaCyeWctk1qZR70cMfLeUAJl78b0q5VVtfWcpd9xulGOM33FMLAt9e/2ul3KvWfamUK085Byy+pDYN4Kce/Gel3IUbf7+UA7hz5++VcgsuuqSUO2/rn5ZyANvOf18pt3bcS0q5s7unlHIAk8tJtcUKVpLUMPfBgh2sJKlpiafpcAQHOUXEbRGxKSIe6bPsDyLiJxGxsnN7U7ubKUnS6HIkFeyngP8KfPqg5R/LzD9qfIskSaOcQ8RwBBVsZn4L2HoUtkWSpGPGUM6DvTkiHuoMIU9tbIskSaNeZE/jt9Gm2sF+AjgXWAxsBP64vxUj4qaIWBERKzZvea7YnCRp9OgMETd9G2VKHWxmPpOZ3ZnZA/wFcNlh1l2WmUsyc8n0U+vneEmSNJqUTtOJiFmZubHz9JeBRw63viTpOJJAz+irOJs2YAcbEZ8FrgSmR8R64MPAlRGxmN6P8QmgfvkiSZKOQQN2sJl53SEWf7KFbZEkHRMSsmu4N2LYeSUnSVLD0iFinK5OkqRWHNUKticm8sKE+YPOPbtjT7nNM2bVZqjh5LeXYi9bOK/WHrBrd21I5YXuuaXcrEXvKeUAJjzzzVJu4aSxtQY37CjFtpz6S7X2gHNOrv0Fvu+UG2oNjj2hlqM+K847L60dPvG/n/zrUg7ghIeeLuXOPW17Kdc9/dpSDmDSxFoNsnDR9FLuhJ0/KuUAuhj879ZWjcLTappmBStJOuZFxLSIuDsi1nTuf+YCSRGxOCK+ExGPdi6k9Kt9XvtURPy4zzX4Fw/Uph2sJKlZ2dkH2/RtaG4B7snMBcA9necH2wW8IzNfAiwF/lNE9L2Aw7/OzMWd28qBGrSDlSQdD64Bbu88vh14y8ErZOYPM3NN5/EGYBNwWrVBO1hJUsM6p+k0fRuamQcukNS5n3G4lSPiMmAC0HfH+Ec7Q8cfi4gBD5rwNB1JUvPaOU1nekSs6PN8WWYuO/AkIr4OnH6I3AcH00hEzAL+Gri+c0lggA8AT9Pb6S4D3g985HDvYwcrSRotNmfmkv5ezMyr+3stIp45cJnfTge6qZ/1XgR8Bfi3mfndPu994PLAeyPir4B/NdDGOkQsSWrWyDzIaTlwfefx9cCXD14hIiYAXwQ+nZn/46DXZnXug979twNeg98OVpJ0PLgVeH1ErAFe33lORCyJiL/srPN24NXAOw9xOs7fRMTDwMPAdOA/DNSgQ8SSpIbliLvQRGZuAa46xPIVwI2dx58BPtNP/nWDbdMOVpLUvB4v9u8QsSRJLbCClSQ1K0feEPFwsIKVJKkFR7WCHUMXk+K5QeeeWP1Cuc15l+4s5bac+LJSbvpJ9Y904qb7Srm1XbVt3R2nlnIA42e9sZSbuOrjpVzPnNeWcvv39wy8Uj/O3H1PKdc9cVEpt//EOaUcwIUbf7+Uq86K80tzf7OUA/jers+XcuPW1LZ108zrSjmAfH5/KTdzw221BidMruWArskjaTYd54MFh4glSS1Ih4gdIpYkqQ1WsJKkZmV6mg5WsJIktcIKVpLUMA9yAitYSZJaYQUrSWqeRxHbwUqSGpYOEYNDxJIktcIKVpLUsIT0NB0rWEmSWmAFK0lqnvtg7WAlSU1zujpwiFiSpFYc1Qp2z77gsSfGDjp30ZIzym3+cMueUu788V8r5baNqU3jBrB660Wl3IVTvlNrcN0ztRywffa1pdyuee8u5abteqCUG3dClHIAq/a+spRb+P0Pl3I9l9dyAHfu/L1S7oSHni7lqlPOAVw26e2l3FfWf6aUO6P7kVIO4HuP1373TL209j2f8OB/LOVGnMQhYqxgJUlqhftgJUkNczYdsIOVJDXOg5zAIWJJklphBStJapYHOQFWsJIktcIKVpLUMPfBgh2sJKlxSXoUsUPEkiS1wQpWktSoTEgPcrKClSSpDVawkqSGpRUsVrCSJLXiqFawEyf0cMFZg5/dZtyP7yi3uXPa9aXcP+5+VSl3/po/LeUANp58Qy24bW0ptu/F76y1B6x/cmcp9+yG7aXca068t5SbPqvWHsCKNQtKuTNe80el3JRtf1fKASy46JJS7tzTap/PuDV/XcpBfVacX5j9G6Xc/7uuvq379m4p5VY9tq2Uu/DFby3lRpy0ggWHiCVJLchuO1iHiCVJaoEVrCSpWQ4RA1awkiS1wgpWktSo9DQdwA5WktSCHjtYh4glSWqDFawkqVke5ARYwUqS1AorWElS46xg7WAlSQ3LdMJ1cIhYknQciIhpEXF3RKzp3E/tZ73uiFjZuS3vs3xeRNzfyX8uIiYM1KYdrCSpYb0HOTV9G6JbgHsycwFwT+f5oezOzMWd25v7LP9D4GOd/DZgwNlZju4QcXYzZv/zg45tOrM4ywwwZm/th/LQd58o5c6/clEpBzDjpBNrwU2TS7G7v7q+1h7w+qWzS7kZp82sNbj94lLsybGvrrUHnHdRlnJPP73lqtEwAAAMxElEQVSrlNv5ostLOYDzttZmceqefm0pt2nmdaUcwBndj5Ry1VlxfnHOb5ZyAH/+7f9ayi06a0cp9+SOOaUcwLzd3xp0ZkzXC+X2RqFrgCs7j28H7gPefyTBiAjgdcCv9cn/AfCJw+XcBytJalaOyIOcZmbmRoDM3BgRM/pZb2JErAC6gFsz80vAqcBzmXlgx/J64MyBGrSDlSSNFtM7nd8ByzJz2YEnEfF14PRD5D44iDbOyswNEXEOcG9EPAwcauLkAYe47GAlSQ3LtuaD3ZyZS/ptNfPq/l6LiGciYlanep0FbOrnPTZ07h+PiPuAi4H/CUyJiHGdKnY2sGGgjR3wIKeImBMR34iIxyLi0Yj4F53lR3REliTp+JI5Ig9yWg5c33l8PfDlg1eIiKkRcULn8XTgCmBVZibwDeCth8sf7EiOIu4Cfi8zLwB+DnhvRCzkyI/IkiRpuN0KvD4i1gCv7zwnIpZExF921rkAWBER/0Bvh3prZq7qvPZ+4H0RsZbefbKfHKjBAYeIOzuFD+wY3hERj9G7c7d8RJYk6dg20mbTycwtwFWHWL4CuLHz+O+BQ54KkpmPA5cNps1BnQcbEWfTOx59PwcdkQX0d0SWJEnHnSM+yCkiTqZ3R++/zMztvacFHVHuJuAmgLNmz6psoyRpVHE2HTjCDjYixtPbuf5NZv6vzuIjPSJrGbAM4NLFC2tn7kuSRo+ReR7sUXckRxEHvTtzH8vMP+nz0oBHZEmSdLw6kgr2CuA3gYcjYmVn2e/TewTW5yPiBuAp4G3tbKIkaTRJh4iBIzuK+P8D+tvh+jNHZEmSJK/kJElqWlrBgh2sJKkFTrh+lDvY2L+TMU/fP/jgzHPKbW5/fm8p96ZfeUmtvZ6FpRzA1K7i9HHFvxSv/oX5tfaAvft6Srknn9hayk3f9tVSbu7kh0o5gO/xz0u5S+cf6rrgAxu3uz594Lbz31fKTZpYmxI6n99fygF87/EzSrl9e7eUctUp5wDe/aqbS7kHn/m/Srl5Y1YMvFI/9s4Y/NSMPeNOLrengVnBSpKala1d7H9Uqf35KkmSDssKVpLUqCRH3LWIh4MVrCRJLbCClSQ1ztN07GAlSU3zPFjAIWJJklphBStJapwVrBWsJEmtsIKVJDUqM+npcfpvO1hJUuPsXx0iliSpFVawkqTGpSXsUe5gx54AU84ddGzGjrvLTZ56+vml3Jqf1Ir7OXPqs1M8u3tWKTdt/q+XcmP2P1fKAezec1Ip9+yG2kwzzJpXiu176Xtr7QETVr9Qyj35Qm22mHM2fK6UA1g7rjb708JF00u5mRtuK+UApl767lJu1WPbSrlFZ+0o5aA+K86lMz9Qyn3y/k+UcgAXTRn8bEORTinXJitYSVKjMqGnNqPlMcUOVpLUuJ50iNiDnCRJaoEVrCSpcQ4RW8FKktQKK1hJUqMyPU0HrGAlSWqFFawkqXEWsHawkqQWeLF/h4glSWqFFawkqVGJQ8RgBStJUiusYCVJzcr0NB2Odge77wVY/+1Bx3Lu1fU2x4wvxcZPqBX3e/d2l3IA0yduLeXGbltfyj0di0o5gFnP3F7KzTn3baVc11nvKOXGb1tZygEsPP/SUu6b964r5ea+vPZvBDi7e0opd8LOH9UanDC5lgMmPPgfS7kLX/zWUu7JHXNKOYB5Y1aUctVZcW64/D2lHMAPNvzbQWeiZ3+5vYF4kJNDxJIktcIhYklSozzIqZcVrCRJLbCClSQ1K90HC3awkqQWpNPVOUQsSVIbrGAlSY3qPcjJIWIrWEmSWmAFK0lqVkKP+2CtYCVJx76ImBYRd0fEms791EOs89qIWNnntici3tJ57VMR8eM+ry0eqE07WElSo5Le03Savg3RLcA9mbkAuKfz/Ke3O/Mbmbk4MxcDrwN2AV/rs8q/PvB6Zg54HVY7WElS43qy+dsQXQMcuIj67cBbBlj/rcDfZuauaoN2sJKk0WJ6RKzoc7tpENmZmbkRoHM/Y4D1rwU+e9Cyj0bEQxHxsYg4YaAGPchJktSs9qar25yZS/p7MSK+Dpx+iJc+OJhGImIWsAi4q8/iDwBPAxOAZcD7gY8c7n2OagfbM+l09l78M8PeI9LsSUe/ze5Dfi+OIDexlptWSvXaO/XGUm5+sb3qJIDdp55aTNa9+qqzS7muIbRZnTyuq/gT6Zpc/UkeffN2f6uc3Tvj1aXcRVO2lHKVKecOuPiM/zDozFo2ltsbiTKz37lNI+KZiJiVmRs7Heimw7zV24EvZuY/zed3oPoF9kbEXwH/aqDtcYhYktSoEXqQ03Lg+s7j64EvH2bd6zhoeLjTKRMRQe/+20cGatAhYklS40bgtf5vBT4fETcATwFvA4iIJcBvZeaNnednA3OAbx6U/5uIOA0IYCXwWwM1aAcrSTrmZeYW4KpDLF8B3Njn+RPAmYdY73WDbdMOVpLULKerA9wHK0lSK6xgJUmNSpwPFuxgJUktcLo6h4glSWqFFawkqVlOVwdYwUqS1AorWElSo5IReaGJo84KVpKkFljBSpIaZwV7lDvYH/zgB5snTZr0ZD8vTwc2H83tGWX8fPrnZ9M/P5vDO94/n7ltvKlDxL2Oagebmaf191pErDjcPH/HOz+f/vnZ9M/P5vD8fNQmh4glSc1KK1jwICdJkloxkirYZcO9ASOcn0///Gz652dzeH4+LbGCHUEdbGb6RT8MP5/++dn0z8/m8Px82tF7kFMM92YMO4eIJUlqwYjoYCNiaUSsjoi1EXHLcG/PSBIRT0TEwxGxMiJWDPf2DLeIuC0iNkXEI32WTYuIuyNiTed+6nBu43Dp57P5g4j4Sef7szIi3jSc2zhcImJORHwjIh6LiEcj4l90lvvdaUPnWsRN30abYe9gI2Is8HHg54GFwHURsXB4t2rEeW1mLvZ0AgA+BSw9aNktwD2ZuQC4p/P8ePQpfvazAfhY5/uzODPvPMrbNFJ0Ab+XmRcAPwe8t/N7xu+OWjPsHSxwGbA2Mx/PzH3AHcA1w7xNGqEy81vA1oMWXwPc3nl8O/CWo7pRI0Q/n42AzNyYmd/vPN4BPAacid+dVhy40ETTt9FmJHSwZwLr+jxf31mmXgl8LSIejIibhntjRqiZmbkRen+RAjOGeXtGmpsj4qHOEPJxPwQaEWcDFwP343enNXawI6ODPdShZqPwo2zNFZl5Cb1D6O+NiFcP9wZpVPkEcC6wGNgI/PHwbs7wioiTgf8J/MvM3D7c26Nj20joYNcDc/o8nw1sGKZtGXEyc0PnfhPwRXqH1PXTnomIWQCd+03DvD0jRmY+k5ndmdkD/AXH8fcnIsbT27n+TWb+r85ivzstcIi410joYB8AFkTEvIiYAFwLLB/mbRoRIuKkiJh84DHwBuCRw6eOS8uB6zuPrwe+PIzbMqIc6Dw6fpnj9PsTEQF8EngsM/+kz0t+d9SaYb/QRGZ2RcTNwF3AWOC2zHx0mDdrpJgJfLH3dwPjgP+emV8d3k0aXhHxWeBKYHpErAc+DNwKfD4ibgCeAt42fFs4fPr5bK6MiMX0FhVPAO8etg0cXlcAvwk8HBErO8t+H787rRmNFWfTItNPQZLUnLOnjMn/85UTGn/fG7+y98HRdLrisFewkqRjzCjdZ9o0O1hJUqMS6LaDHREHOUmSdMyxgpUkNc4hYitYSZJaYQUrSWpUMjpnv2maHawkqVkeRQw4RCxJUiusYCVJjTpwLeLjnRWsJEktsIKVJDXOCtYOVpLUgp481FTfxxeHiCVJaoEVrCSpUR7k1MsKVpKkFljBSpKa5YUmACtYSZJaYQUrSWqU1yLuZQcrSWqcQ8QOEUuSjgMR8baIeDQieiJiyWHWWxoRqyNibUTc0mf5vIi4PyLWRMTnImLCQG3awUqSGteTzd+G6BHgnwHf6m+FiBgLfBz4eWAhcF1ELOy8/IfAxzJzAbANuGGgBu1gJUnHvMx8LDNXD7DaZcDazHw8M/cBdwDXREQArwO+0FnvduAtA7XpPlhJUqNy9J6mcyawrs/z9cDlwKnAc5nZ1Wf5mQO9mR2sJKlRG3fmXf/u7/ZMb+GtJ0bEij7Pl2XmsgNPIuLrwOmHyH0wM798BO9/qAso52GWH5YdrCSpUZm5dJjavXqIb7EemNPn+WxgA7AZmBIR4zpV7IHlh+U+WEmSej0ALOgcMTwBuBZYnpkJfAN4a2e964EBK2I7WEnSMS8ifjki1gMvB74SEXd1lp8REXcCdKrTm4G7gMeAz2fmo523eD/wvohYS+8+2U8O2GZvxyxJkppkBStJUgvsYCVJaoEdrCRJLbCDlSSpBXawkiS1wA5WkqQW2MFKktQCO1hJklrw/wP4fKAwlAPfbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "truCov = np.eye(nbPix, nbPix)\n",
    "np.random.seed(121)\n",
    "for i in range(nbPix):\n",
    "    for j in range(i+1,nbPix):\n",
    "        truCov[i,j] = (np.random.rand(1) - .5)\n",
    "        truCov[j,i] = truCov[i,j]\n",
    "truCov = truCov.T.dot(truCov)\n",
    "# I do want some uncorrelated entries:\n",
    "truCov[-1,:-5] = 0\n",
    "truCov[:-5,-1] = 0\n",
    "# unit variance\n",
    "truCov /= np.max(truCov)\n",
    "\n",
    "for j in range(nbPix):\n",
    "    truCov[j,j] = 1\n",
    "\n",
    "plt.imshow(truCov, cmap=cmap, interpolation='Nearest', vmin=-1, vmax=1)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5, 5, 7, 7) 1225\n"
     ]
    }
   ],
   "source": [
    "covariances = np.zeros((*(dataSize),\n",
    "                    2*howFar+1, 2*howFar+1)) \n",
    "print(covariances.shape, covariances.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make nosie realizations and compute covariance estimator on each:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " > Noise realization 0 generated in 0.0007195472717285156s\n",
      " > Cross terms computed in 0.001476287841796875s\n",
      " > Noise realization 500 generated in 0.0004246234893798828s\n",
      " > Cross terms computed in 0.0014433860778808594s\n",
      " > Noise realization 1000 generated in 0.00040411949157714844s\n",
      " > Cross terms computed in 0.001466989517211914s\n",
      " > Noise realization 1500 generated in 0.0004036426544189453s\n",
      " > Cross terms computed in 0.001434326171875s\n",
      " > Noise realization 2000 generated in 0.00040650367736816406s\n",
      " > Cross terms computed in 0.0014178752899169922s\n",
      " > Noise realization 2500 generated in 0.0004303455352783203s\n",
      " > Cross terms computed in 0.0014009475708007812s\n",
      " > Noise realization 3000 generated in 0.0004742145538330078s\n",
      " > Cross terms computed in 0.0015568733215332031s\n",
      " > Noise realization 3500 generated in 0.00040602684020996094s\n",
      " > Cross terms computed in 0.0014119148254394531s\n",
      " > Noise realization 4000 generated in 0.0004649162292480469s\n",
      " > Cross terms computed in 0.0015225410461425781s\n",
      " > Noise realization 4500 generated in 0.0004050731658935547s\n",
      " > Cross terms computed in 0.003044605255126953s\n",
      " > Noise realization 5000 generated in 0.00040602684020996094s\n",
      " > Cross terms computed in 0.0014095306396484375s\n",
      " > Noise realization 5500 generated in 0.0004277229309082031s\n",
      " > Cross terms computed in 0.0014951229095458984s\n",
      " > Noise realization 6000 generated in 0.0004062652587890625s\n",
      " > Cross terms computed in 0.0015759468078613281s\n",
      " > Noise realization 6500 generated in 0.0004048347473144531s\n",
      " > Cross terms computed in 0.0014445781707763672s\n",
      " > Noise realization 7000 generated in 0.00046181678771972656s\n",
      " > Cross terms computed in 0.001489400863647461s\n",
      " > Noise realization 7500 generated in 0.0004050731658935547s\n",
      " > Cross terms computed in 0.0014145374298095703s\n",
      " > Noise realization 8000 generated in 0.00040459632873535156s\n",
      " > Cross terms computed in 0.0014297962188720703s\n",
      " > Noise realization 8500 generated in 0.0004150867462158203s\n",
      " > Cross terms computed in 0.0014147758483886719s\n",
      " > Noise realization 9000 generated in 0.0004134178161621094s\n",
      " > Cross terms computed in 0.0014371871948242188s\n",
      " > Noise realization 9500 generated in 0.0004055500030517578s\n",
      " > Cross terms computed in 0.0014088153839111328s\n",
      "All done in 16.100027799606323s\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "for j in range(nbNoiseReals):\n",
    "    startReal = time.time()\n",
    "    noiseReal = np.random.multivariate_normal(mean, truCov, 1)[0].reshape(dataSize)\n",
    "    startCovComp = time.time()\n",
    "    if printCheckpoints and not j%500:\n",
    "        print(' > Noise realization {} generated in {}s'.format(\n",
    "              j, startCovComp-startReal))\n",
    "    for xidx,dx in enumerate(range(-howFar,howFar+1)):\n",
    "        noiseShiftx = np.roll(noiseReal, shift=dx, axis=1)\n",
    "        for yidx,dy in enumerate(range(-howFar,howFar+1)):\n",
    "            noiseShiftxy = np.roll(noiseShiftx, shift=dy, axis=0)\n",
    "            covariances[:,:,xidx,yidx] += noiseShiftxy * noiseReal\n",
    "    if printCheckpoints and not j%500:\n",
    "        print(' > Cross terms computed in {}s'.format(\n",
    "              time.time()-startCovComp))\n",
    "covariances /= (nbNoiseReals - 1)\n",
    "print('All done in {}s'.format(time.time()-start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test our estimator\n",
    "For each shift considered, we can access the corresponding entry of the covariance for every single pixel from the last two axes.\n",
    "\n",
    "For instance, the diagonal corresponds to a shift of $0,0$, ie:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fdaf0e18c50>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAAHVCAYAAABFZXqOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADXhJREFUeJzt3N+LpQd9x/HvNzu7GLHgRXMh2dAIFWkQmsAQhPQqeLH+QG8T0CthbypEEEQv/QfEG28WDRYUg6AXEiwSMEECNrqJUUxXIYjFRWFbRDQU3O7utxc7lLTdOmf7Oc8+e868XjAwZ/bh8OHZw7znmXlmemYKAPj/uWftAQCwy4QUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAAQOFnnSPjNn6q1LPDUALO5q/Xtdm6u9ybGLhPRMvbX+uv5uiaeGxXT7Bs1S7rnHuV3K9evX1p6wl16vFzc+1qsbAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAENgopN19rrt/0d2vd/dnlh4FALvi2JB296mq+mJVvb+qHqqqJ7v7oaWHAcAu2OSK9NGqen1mfjkzV6vqmar6yLKzAGA3bBLS+6vq1296fPnoY/9Nd5/v7ovdffFaXd3WPgC4q20S0r7Fx+Z/fWDmwswczszhQZ3JlwHADtgkpJer6oE3PT5bVb9ZZg4A7JZNQvqjqnpXd7+zu89U1RNV9e1lZwHAbjg47oCZudbdn6iq71bVqap6emZeW3wZAOyAY0NaVTUz36mq7yy8BQB2jr9sBAABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgcrD2A29Ptax92T/eptSfsre4ba0/YT7P5oT4rA0BASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgCBY0Pa3U9395Xu/tmdGAQAu2STK9KvVNW5hXcAwE46NqQz8/2q+t0d2AIAO+dgW0/U3eer6nxV1em6d1tPCwB3ta3dbDQzF2bmcGYOD+rMtp4WAO5q7toFgICQAkBgk19/+XpV/aCq3t3dl7v748vPAoDdcOzNRjPz5J0YAgC7yLd2ASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAIGDtQdwe2ZurD0Bbtv16/+x9oS95XPC+lyRAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgCBY0Pa3Q909/Pdfam7X+vup+7EMADYBQcbHHOtqj41M690919U1cvd/dzM/PPC2wDgrnfsFenM/HZmXjl6/49Vdamq7l96GADsgk2uSP9Ldz9YVY9U1Uu3+LfzVXW+qup03buFaQBw99v4ZqPufltVfbOqPjkzf/if/z4zF2bmcGYOD+rMNjcCwF1ro5B29+m6GdGvzcy3lp0EALtjk7t2u6q+XFWXZubzy08CgN2xyRXpY1X1sap6vLtfPXr7wMK7AGAnHHuz0cy8WFV9B7YAwM7xl40AICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAgYOlnrhbo5cwc2PtCXDbvG7ZZ2oHAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASODWl3v6W7f9jdP+nu17r7c3diGADsgoMNjvlTVT0+M2909+mqerG7/3Fm/mnhbQBw1zs2pDMzVfXG0cPTR2+z5CgA2BUb/Yy0u09196tVdaWqnpuZl25xzPnuvtjdF6/V1W3vBIC70kYhnZnrM/NwVZ2tqke7+z23OObCzBzOzOFBndn2TgC4K93WXbsz8/uqeqGqzi2yBgB2zCZ37d7X3W8/ev/eqnpfVf186WEAsAs2uWv3HVX1D919qm6G9xsz8+yyswBgN2xy1+5Pq+qRO7AFAHaOv2wEAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACBws9cQzN5Z66hPt1KnF/stOvOvXr609YW953S7H63Z9rkgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgMDGIe3uU9394+5+dslBALBLbueK9KmqurTUEADYRRuFtLvPVtUHq+pLy84BgN2y6RXpF6rq01V14/86oLvPd/fF7r54ra5uZRwA3O2ODWl3f6iqrszMy3/uuJm5MDOHM3N4UGe2NhAA7mabXJE+VlUf7u5fVdUzVfV4d3910VUAsCOODenMfHZmzs7Mg1X1RFV9b2Y+uvgyANgBfo8UAAIHt3PwzLxQVS8ssgQAdpArUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkBASAEgIKQAEBBSAAgIKQAEhBQAAkIKAAEhBYCAkAJAQEgBICCkABAQUgAICCkABIQUAAJCCgABIQWAgJACQEBIASAgpAAQEFIACAgpAASEFAACQgoAASEFgICQAkCgZ2b7T9r9r1X1L1t/4mX8ZVX929oj9pRzuxzndjnO7XJ26dz+1czct8mBi4R0l3T3xZk5XHvHPnJul+PcLse5Xc6+nlvf2gWAgJACQEBIqy6sPWCPObfLcW6X49wuZy/P7Yn/GSkAJFyRAkBASAEgcKJD2t3nuvsX3f16d39m7T37oruf7u4r3f2ztbfsm+5+oLuf7+5L3f1adz+19qZ90d1v6e4fdvdPjs7t59betE+6+1R3/7i7n117y7ad2JB296mq+mJVvb+qHqqqJ7v7oXVX7Y2vVNW5tUfsqWtV9amZ+Zuqem9V/b3X7db8qaoen5m/raqHq+pcd7935U375KmqurT2iCWc2JBW1aNV9frM/HJmrlbVM1X1kZU37YWZ+X5V/W7tHftoZn47M68cvf/HuvmJ6f51V+2HuemNo4enj97cjbkF3X22qj5YVV9ae8sSTnJI76+qX7/p8eXyCYkd0t0PVtUjVfXSukv2x9G3H1+tqitV9dzMOLfb8YWq+nRV3Vh7yBJOckj7Fh/z1Sc7obvfVlXfrKpPzswf1t6zL2bm+sw8XFVnq+rR7n7P2pt2XXd/qKquzMzLa29ZykkO6eWqeuBNj89W1W9W2gIb6+7TdTOiX5uZb629Zx/NzO+r6oXys/5teKyqPtzdv6qbP0J7vLu/uu6k7TrJIf1RVb2ru9/Z3Weq6omq+vbKm+DP6u6uqi9X1aWZ+fzae/ZJd9/X3W8/ev/eqnpfVf183VW7b2Y+OzNnZ+bBuvl59nsz89GVZ23ViQ3pzFyrqk9U1Xfr5g0b35iZ19ZdtR+6++tV9YOqend3X+7uj6+9aY88VlUfq5tf1b969PaBtUftiXdU1fPd/dO6+YX2czOzd7+qwfb5E4EAEDixV6QAsA1CCgABIQWAgJACQEBIASAgpAAQEFIACPwn5ajmumbHk5cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(covariances[:,:,3,3], cmap=cmap, interpolation='Nearest', vmin=-1, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conversely, we can access the covariance entries with all considered pixels for one given pixel; for instance that last pixel with several 0 covariance counterparts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAHLCAYAAAB1UrkfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+w5XV93/Hna1cWYn4UdFERloAJEyVRcbOD6WgrIhjMdFzSmAhNlTQwTDLQ1ObHiLWVBusMppnYdIa02QqRJI5ojMZtxKGIUhsNhtVsQEDCBhvZWRJc0Bgjgrv77h/nu+nx5t7de+953z333H0+Zr5zz/f3+8Dcfd/3+/v5fr+pKiRJUo910w5AkqS1xMQqSVIjE6skSY1MrJIkNTKxSpLUyMQqSVIjE6skaaYluSHJI0k+t8D6JPmvSXYluSvJ5rF1lyR5YJgu6YjHxCpJmnXvAi44xPpXAWcM0+XAfwNI8jTgauDFwNnA1UlOmDQYE6skaaZV1SeAxw6xyVbgt2vkDuD4JCcBPwzcWlWPVdWXgVs5dIJeFBOrJGmtOxl4aGx+97BsoeUTecqkB5AkCeA784zaz5Otx3ycv7kH+MbYom1VtW2Jh8k8y+oQyydiYpUktdjPk3wvL2095t18+BtVtWXCw+wGNo3NnwLsGZafM2f57ROey1awJKnPunXrWqcm24HXD6ODfwj4m6p6GLgFeGWSE4ZBS68clk3EilWSNNOSvIdR5bkxyW5GI32PAaiq/w7cDPwIsAv4OvCvhnWPJXkrcOdwqGuq6lCDoBYXj6+NkyR1eGpOqO9bf07rMXfu/4PPNLSCjyhbwZIkNbIVLElqkdB5XXRkf+/hjgQTqySpSUjWTzuIqbMVLElSIytWSVKb9lbwDPK/gCRJjaxYJUlNQmK9ZmKVJLUYjQp28JJ/WkiS1MiKVZLUxFYwWLFKktTKilWS1MZrrCZWSVKTxFYw2AqWJKmVFaskqY1PXrJilSSplRWrJKmJb7cBK1ZJklpZsUqSWqzIi85nkIlVktTG221sBUuS1MqKVZLUJD55CStWSZJaWbFKktp4jdXEKklqktgKBlvBkiS1smKVJLXxPlYrVkmSWlmxSpKa+D5WMLFKkpqMHmno4CX/tJAkqZEVqySpia1gsGKVJKmVFaskqY3XWI9AYt24cWOdeuqpK32aIyLUtENoU2TaIbRJ7Zt2CH2ytv7WXTu/MWvn9/8vv/hF9u59dEX+AUhsBcMRSKynnnoqn/zkJ1f6NEdEDjw57RDa7OOYaYfQZsM39047hDYHNjxt2iG0qqyd6mWt/P6/5J+cM+0Q1ry19eexJGmqfPKSg5ckSWplxSpJahKyhtr/y2XFKklSIytWSVKL0SMNrddMrJKkNuvWr51b+ZbLPy0kSWpkxSpJ6pGQ9dZr/heQJKmRFaskqUXwGiuYWCVJXTKdxJrkAuDXgfXAO6vq2jnr3wG8fJh9KvCMqjp+WLcfuHtY98WqevWk8ZhYJUkzK6MnUlwHnA/sBu5Msr2q7j24TVX927Ht/zXworFDPF5VZ3XGZGKVJPWYzuCls4FdVfXgKITcBGwF7l1g+4uBq1cyIAcvSZJWs41JdoxNl89ZfzLw0Nj87mHZP5Dku4HTgY+NLT5uOO4dSS7sCNiKVZLUYoUGL+2tqi2HOe1cC7089yLg/VW1f2zZqVW1J8lzgI8lubuq/mK5wYKJVZLUJbDuyLeCdwObxuZPAfYssO1FwBXjC6pqz/DzwSS3M7r+OlFitRUsSZpldwJnJDk9yQZGyXP73I2SfB9wAvDHY8tOSHLs8Hkj8BIWvja7aFaskqQW07iPtar2JbkSuIXR7TY3VNU9Sa4BdlTVwSR7MXBTVY23iZ8H/GaSA4wKzWvHRxMvl4lVkjTTqupm4OY5y94yZ/4/zrPfp4Dnd8djYpUk9UiIT15a2jXWJFuT3JVk5zA8+aUrFZgkSbNoqRXrbcD2qqokLwDeBzy3PyxJ0syZzqjgVWdJibWqvjY2++0sfK+QJOko40P4R5b8p0WSH03yeeDDwE8vsM3lB5+SsXfv3kljlCRpZiw5sVbVB6vqucCFwFsX2GZbVW2pqi0bN26cNEZJ0ixIWLd+Xes0iw4bdZIrhsFKO5M8++DyqvoE8D3DTbWSJIlFXGOtqusYvZKHJN+bJMPgpc3ABuDRFY5RkjQjvN1m6aOCfwx4fZJvAo8Dr53zFAtJ0lEqU3rR+Wqz1FHBbwfevkKxSJI083zykiSph/exAr7dRpKkVlaskqQmPisYTKySpCaxFQzYCpYkqZUVqySph7fbAFaskiS1smKVJDWJ11ixYpUkqZUVqySpReKzgsHEKklq5OAlW8GSJLWyYpUk9YiDl8CKVZKkVlaskqQWvo91xMQqSWoTW8G2giVJ6mTFKknqYSsYsGKVJKnVilesAUKt9GmOjHXHTDuCNutq7fxVuf/YE6cdQps187syWFPfZ838/q/c734I69atnX9blstWsCSpR2DdOhuh/heQJKmRFaskqUVw8BJYsUqS1MqKVZLUI/EaK1askiS1smKVJLXxdhsTqySpSfB2G7AVLElSKytWSVKP2AoGK1ZJklpZsUqSWoSwzvexmlglSU1sBQO2giVJamXFKklq4+02VqySpBmX5IIk9yfZleSqedb/VJIvJdk5TJeNrbskyQPDdElHPFaskqQWyZF/0XmS9cB1wPnAbuDOJNur6t45m763qq6cs+/TgKuBLUABnxn2/fIkMZlYJUltptAKPhvYVVUPAiS5CdgKzE2s8/lh4NaqemzY91bgAuA9kwRkK1iSNMtOBh4am989LJvrx5LcleT9STYtcd8lMbFKknoEsj6tE7AxyY6x6fJ/eNZ/oObM/0/gtKp6AfBR4MYl7LtktoIlSavZ3qracoj1u4FNY/OnAHvGN6iqR8dm/wfw9rF9z5mz7+3LDfQgK1ZJUoswetF557QIdwJnJDk9yQbgImD7t8SVnDQ2+2rgvuHzLcArk5yQ5ATglcOyiSypYk3yk8Abh9mvAT9bVX82aRCSJC1HVe1LciWjhLgeuKGq7klyDbCjqrYDP5fk1cA+4DHgp4Z9H0vyVkbJGeCagwOZJrHUVvAXgJdV1ZeTvArYBrx40iAkSWvDNB5pWFU3AzfPWfaWsc9vAt60wL43ADd0xrOkxFpVnxqbvYNRP1qSJBKfvASTXWO9FPjIfCuSXH5wBNeX9u6d4BSSJM2WZY0KTvJyRon1pfOtr6ptjNrE/ODmzRMPXZYkzQbfbrOIijXJFWPPV3x2khcA7wS2zhnCLEnSUe+wFWtVXcfoOYwkORX4APC6qvrzFY5NkjRDEl90DktvBb8FeDrwG0kA9h3mxl1J0lHEVvDSRwVfBlx22A0lSTpK+UhDSVIPb7cBfKShJEmtrFglSS1Gzwr2GquJVZLUw1YwYCtYkqRWVqySpBYB1q23FWzFKklSIytWSVKPxGusWLFKktTKilWS1MbbbUyskqQmwcQKtoIlSWplxSpJ6hGIFasVqyRJnaxYJUltvMZqYpUkNQkhMbHaCpYkqZEVqySpR2wFgxWrJEmtVrxiLaBYG3/BVE07gj7r2D/tENoc+8j/mXYIbR4/8Zxph9BqLV1uWyu//yv9NbzdxlawJKlJbAUDtoIlSWplxSpJamMr2IpVkqRWVqySpCbxGitWrJIktbJilST18O02gIlVktQkwLq1dPPyMtkKliSpkRWrJKmNrWArVkmSWlmxSpJ6+EhDwMQqSWoSYisYW8GSJLWyYpUk9bAVDFixSpLUyopVktTGa6wmVklSE5+8NGIrWJI005JckOT+JLuSXDXP+p9Pcm+Su5LcluS7x9btT7JzmLZ3xGPFKknqMYWH8CdZD1wHnA/sBu5Msr2q7h3b7E+BLVX19SQ/C/wK8Nph3eNVdVZnTFaskqRZdjawq6oerKongZuAreMbVNXHq+rrw+wdwCkrGZAVqySpyVRedH4y8NDY/G7gxYfY/lLgI2PzxyXZAewDrq2qP5g0oCUl1iTPBX4L2Ay8uap+ddIAJElrxwq0gjcOie+gbVW1bfyU8+xT8x0oyb8EtgAvG1t8alXtSfIc4GNJ7q6qv5gk4KVWrI8BPwdcOMlJJUlapL1VteUQ63cDm8bmTwH2zN0oyXnAm4GXVdUTB5dX1Z7h54NJbgdeBEyUWJd0jbWqHqmqO4FvTnJSSdLak+HJS53TItwJnJHk9CQbgIuAbxndm+RFwG8Cr66qR8aWn5Dk2OHzRuAlwPigp2VZkWusSS4HLgfYtGnTYbaWJGl5qmpfkiuBW4D1wA1VdU+Sa4AdVbUd+M/AdwC/l9F9tl+sqlcDzwN+M8kBRoXmtXNGEy/LiiTWof+9DWDz5s3z9rolSWvPNJ68VFU3AzfPWfaWsc/nLbDfp4Dnd8dz2FZwkivGbp59dncAkiStJYetWKvqOkY330qSdEg+0nDpt9s8C9gBfBdwIMkbgDOr6qsrEZwkaXYkvugclphYq+qvWOEnVkiSNMt88pIkqY0vOvdZwZIktbJilST1mMLbbVYjE6skqY2Dgm0FS5LUyopVktQijG65OdpZsUqS1MiKVZLUxrFLJlZJUhdHBQO2giVJamXFKklq49glK1ZJklpZsUqSWoR4uw1WrJIktbJilSS1cVCwiVWS1MhWsK1gSZJaWbFKknrE223AilWSpFYrXrEGSO1f6dMcGVk/7QjaFGvnu7zwtOumHUKbO/7unGmH0GrN/O7Dmvr9Xym+3WbEVrAkqU3sg9oKliSpkxWrJKnNOlvBVqySJHWyYpUk9fB2G8DEKklq4kP4R2wFS5LUyIpVktTGgtWKVZKkVlaskqQ2XmO1YpUkqZUVqySpjS86N7FKkprE+1gBW8GSJLWyYpUktYm9YCtWSZI6WbFKktp4jdXEKklqkvjaOLAVLElSKytWSVIbC1YrVkmSWlmxSpLa+KzgZVasSS5Icn+SXUmu6g5KkjSbDj59qWta3DkPnZOSHJvkvcP6Tyc5bWzdm4bl9yf54Y7/BktOrEnWA9cBrwLOBC5OcmZHMJIkLcUic9KlwJer6nuBdwBvH/Y9E7gI+H7gAuA3huNNZDkV69nArqp6sKqeBG4Ctk4aiCRp1oWkd1qExeSkrcCNw+f3A6/I6OBbgZuq6omq+gKwazjeRJaTWE8GHhqb3z0s+3tJLk+yI8mOL+3dO0l8kqSj28aD+WSYLp+z/rA5aXybqtoH/A3w9EXuu2TLGbw0358Q9S0zVduAbQA/uHlzzbO9JGmNSYqk/Z/8vVW15VCnnWfZ3CAW2mYx+y7ZcirW3cCmsflTgD2TBiJJ0jIsJif9/TZJngL8I+CxRe67ZMtJrHcCZyQ5PckGRhd+t08aiCRpxhVQ+3unw1tMTtoOXDJ8fg3wsaqqYflFw6jh04EzgD+Z9D/DklvBVbUvyZXALcB64IaqumfSQCRJs65IHTiyZ1wgJyW5BthRVduB64HfSbKLUaV60bDvPUneB9wL7AOuqFpcNj+UZT0goqpuBm6e9OSSJE1qvpxUVW8Z+/wN4McX2PdtwNs64/HJS5KkPpMXfDPPZwVLktTIilWS1KTgCF9jXY1MrJKkNrEVbCtYkqROVqySpB5VDl7CilWSpFZWrJKkPg5eMrFKkrrYCgZbwZIktbJilSS1OdLPCl6NrFglSWpkxSpJalKA11itWCVJamTFKknqUXi7DSZWSVIbb7cBW8GSJLWyYpUktfF2GytWSZJaWbEuQahph6B5fPbhq6cdQpsnpx2AFrRWfv+zokf3GiuYWCVJbUysYCtYkqRWVqySpBYpBy+BFaskSa2sWCVJTbzGCiZWSVInW8G2giVJ6mTFKklqYisYrFglSWplxSpJahO8xmrFKklSIytWSVITr7GCiVWS1KXwdhtsBUuS1MqKVZLUxFYwWLFKktTKilWS1Ma325hYJUltbAWDrWBJklpZsUqSehRwwIrVilWSpEZWrJKkJgW1b9pBTN2SE2uSG4B/BjxSVT/QH5IkaTaVrWCW1wp+F3BBcxySJK0JS65Yq+oTSU7rD0WSNPO83WZlBi8luTzJjiQ7vrR370qcQpKkQ0rytCS3Jnlg+HnCPNucleSPk9yT5K4krx1b964kX0iyc5jOWsx5VySxVtW2qtpSVVtO3LhxJU4hSVptarjG2jlN5irgtqo6A7htmJ/r68Drq+r7GV3m/C9Jjh9b/0tVddYw7VzMSb3dRpK0Vm0Fbhw+3whcOHeDqvrzqnpg+LwHeAQ4cZKTmlglSU2G2206p8k8s6oeBhh+PuNQGyc5G9gA/MXY4rcNLeJ3JDl2MSddzu027wHOATYm2Q1cXVXXL/U4kqQ1qP92m41JdozNb6uqbQdnknwUeNY8+715KSdJchLwO8AlVX//JoE3AX/FKNluA94IXHO4Yy1nVPDFS91HkqRl2ltVWxZaWVXnLbQuyV8nOamqHh4S5yMLbPddwIeBf19Vd4wd++Hh4xNJfgv4xcUEbCtYktRj9Q1e2g5cMny+BPjQ3A2SbAA+CPx2Vf3enHUnDT/D6Prs5xZzUhOrJGmtuhY4P8kDwPnDPEm2JHnnsM1PAP8U+Kl5bqt5d5K7gbuBjcB/WsxJfVawJKnJ6nofa1U9CrxinuU7gMuGz78L/O4C+5+7nPOaWCVJfQ74EH5bwZIkNbJilST1qNXVCp4WK1ZJkhpZsUqSmvg+VjCxSpIala1gW8GSJHWyYpUk9ajydhusWCVJamXFKklq4uAlsGKVJKmVFaskqY+jgk2skqQmZSsYbAVLktTKilWS1KSgvN3GilWSpEZWrJKkPl5jXfnEWkBl/Uqf5ojIgSenHUKbfRwz7RDarDvumdMOoc061tY/Smvldx/W0u9/reyxHRVsK1iSpE62giVJPQpbwVixSpLUyopVktTEt9uAiVWS1MbBS2ArWJKkVlaskqQeDl4CrFglSWplxSpJauI1VjCxSpLaFOWoYFvBkiR1smKVJLWognLwkhWrJEmdrFglSU3KihUrVkmSWlmxSpJ6lBUrmFglSY1qv4nVVrAkSY2sWCVJPWwFA1askiS1smKVJLUob7cBlpFYk2wCfht4FnAA2FZVv94dmCRp9hwwsS6rYt0H/EJVfTbJdwKfSXJrVd3bHJskSTNnyYm1qh4GHh4+/22S+4CTAROrJB3NHLwETDh4KclpwIuAT89ZfnmSHUl27N27d5JTSJI0U5Y9eCnJdwC/D7yhqr46vq6qtgHbADZv3lwTRShJmhlWrMtMrEmOYZRU311VH+gNSZI0i6p80TksoxWcJMD1wH1V9Wv9IUmSNLkkT0tya5IHhp8nLLDd/iQ7h2n72PLTk3x62P+9STYs5rzLucb6EuB1wLljgfzIMo4jSVpTRoOXOqcJXQXcVlVnALcN8/N5vKrOGqZXjy1/O/COYf8vA5cu5qRLTqxV9UdVlap6wVggNy/1OJIkrbCtwI3D5xuBCxe749CdPRd4/1L398lLkqQeteoGLz1zuEWUqno4yTMW2O64JDsYPafh2qr6A+DpwFeq6uBF492Mbi09LBOrJGk12zgkvYO2DXeeAJDko4yeBDjXm5dwjlOrak+S5wAfS3I38NV5tlvUXS4mVklSk1qJ97HuraotC56x6ryF1iX56yQnDdXqScAjCxxjz/DzwSS3M3o+w+8Dxyd5ylC1ngLsWUzAvt1GktSiatUNXtoOXDJ8vgT40NwNkpyQ5Njh80ZGA3TvraoCPg685lD7z8fEKklaq64Fzk/yAHD+ME+SLUneOWzzPGBHkj9jlEivHXv2/RuBn0+yi9E11+sXc1JbwZKkNqvp7TZV9SjwinmW7wAuGz5/Cnj+Avs/CJy91PNasUqS1MiKVZLUxLfbgIlVktRl9d3HOhW2giVJamTFKklqUbaCAStWSZJaWbFKknqUFSuYWCVJjXzRua1gSZJaWbFKknrUijyEf+ZYsUqS1MiKVZLUoqhV9azgabFilSSpkRWrJKmNt9scgcQaINRKn+bIWHfMtCNos64y7RDa7D/2xGmH0GbN/K4M1tT3WTO//yv4u+99rICtYEmSWtkKliS1sWK1YpUkqZUVqySpRVVx4MAauq6+TCZWSVIb86qtYEmSWlmxSpLalCWrFaskSZ2sWCVJLargwIFpRzF9JlZJUpsDZSvYVrAkSY2sWCVJbWwFW7FKktTKilWS1KLK223AilWSpFZWrJKkNhasJlZJUiMfwm8rWJKkVlaskqQWha1gsGKVJKmVFaskqUeVt9uwjMSa5DjgE8Cxw/7vr6qruwOTJM0eBy8tr2J9Aji3qr6W5Bjgj5J8pKruaI5NkqSZs+TEWlUFfG2YPWaY/BNFko5yDl4aWdbgpSTrk+wEHgFurapPz1l/eZIdSXZ8ae/ejjglSZoJy0qsVbW/qs4CTgHOTvIDc9Zvq6otVbXlxI0bO+KUJK12NbrG2jnNoolGBVfVV5LcDlwAfK4lIknSzCpfG7f0ijXJiUmOHz5/G3Ae8PnuwCRJmkXLqVhPAm5Msp5RYn5fVf1hb1iSpFkzGrw0m+3bTssZFXwX8KIViEWSpJnnk5ckST0KDniN1WcFS5LWpiRPS3JrkgeGnyfMs83Lk+wcm76R5MJh3buSfGFs3VmLOa+JVZLUolh1t9tcBdxWVWcAtw3z3xpz1cer6qzhFtJzga8D/2tsk186uL6qdi7mpCZWSVKbA9U7TWgrcOPw+UbgwsNs/xrgI1X19UlOamKVJK1mGw8+yW+YLl/Cvs+sqocBhp/POMz2FwHvmbPsbUnuSvKOJMcu5qQOXpIk9ViZ18btraotC61M8lHgWfOsevNSTpLkJOD5wC1ji98E/BWwAdgGvBG45nDHMrFKkmZWVZ230Lokf53kpKp6eEicjxziUD8BfLCqvjl27IeHj08k+S3gFxcTk61gSVKLVTh4aTtwyfD5EuBDh9j2Yua0gYdkTJIwuj67qEf3WrFKktqssufmXwu8L8mlwBeBHwdIsgX4maq6bJg/DdgE/O85+787yYlAgJ3AzyzmpCZWSdKaVFWPAq+YZ/kO4LKx+f8LnDzPducu57wmVklSj+G1cUc7r7FKktTIilWS1KLwfaxgYpUkNfK1cbaCJUlqZcUqSerha+MAK1ZJklpZsUqSWhSr7gERU2HFKklSIytWSVIbK9YjkFg/+6d/uvfbnvrtf7nS5wE2AnuPwHmOBL/L6uR3WZ38Lkvz3St1YFvBIyueWKvqxJU+B0CSHYd6Z98s8busTn6X1cnvotXGVrAkqUdZsYKDlyRJarWWKtZt0w6gkd9ldfK7rE5+l1XEihVSPtdRktTgOSesq7e9fEPrMf/FB5/4zKxdd7YVLElSo5lPrEm2Jrkryc4kO5K8dNoxLVeSnxy+y11JPpXkhdOOabmSPDfJHyd5IskvTjueSSS5IMn9SXYluWra8UwiyQ1JHknyuWnHMokkm5J8PMl9Se5J8m+mHdNyJTkuyZ8k+bPhu/zytGNatuFZwZ3TLJr5xArcBrywqs4Cfhp455TjmcQXgJdV1QuAtzLb11seA34O+NVpBzKJJOuB64BXAWcCFyc5c7pRTeRdwAXTDqLBPuAXqup5wA8BV8zw/5cngHOr6oXAWcAFSX5oyjFpAjOfWKvqa/X/LxR/O6N7lGdSVX2qqr48zN4BnDLNeCZRVY9U1Z3AN6cdy4TOBnZV1YNV9SRwE7B1yjEtW1V9gtEfPTOtqh6uqs8On/8WuA84ebpRLU+NfG2YPWaYZvLfsYMPiOicZtHMJ1aAJD+a5PPAhxlVrWvBpcBHph2EOBl4aGx+NzP6D/haleQ04EXAp6cbyfIlWZ9kJ/AIcGtVzex3MbGukcRaVR+squcCFzJqoc60JC9nlFjfOO1YROZZNqO/7mtPku8Afh94Q1V9ddrxLFdV7R8uZ50CnJ3kB6Ydk5ZvJhNrkiuGwUo7kzz74PKhzfU9STZOMbwlmftdkryA0XXirVX16LTjW4qF/r/MuN3AprH5U4A9U4pFY5IcwyipvruqPjDteDpU1VeA25nR6+C2gkdmMrFW1XVVddbwF95TkwQgyWZgAzAzCWnOd3kK8AHgdVX151MObcnGv0tVrZXkcydwRpLTk2wALgK2Tzmmo97wO389cF9V/dq045lEkhOTHD98/jbgPODz041Kk1gLT176MeD1Sb4JPA68dmww06x5C/B04DeGvxX2zdqN0QcleRawA/gu4ECSNwBnzlq7rqr2JbkSuAVYD9xQVfdMOaxlS/Ie4BxgY5LdwNVVdf10o1qWlwCvA+4erk0C/LuqunmKMS3XScCNwwj0dcD7quoPpxzTss1qldnJJy9Jklqcdvy6+g8v7X3y0mUfnr0nL62FilWStBrM8HXRTiZWSVKLAvabWGdz8JIkSauVFaskqY2tYCtWSZJaWbFKkloUs/tGmk4mVklSD0cFA7aCJUlqZcUqSWpx8FnBRzsrVkmSGlmxSpLaWLGaWCVJjQ7UfK8wPrrYCpYkqZEVqySphYOXRqxYJUlqZMUqSerhAyIAK1ZJklpZsUqSWvis4BETqySpja1gW8GSpDUqyY8nuSfJgSRbDrHdBUnuT7IryVVjy09P8ukkDyR5b5INizmviVWS1OZA9U4T+hzwz4FPLLRBkvXAdcCrgDOBi5OcOax+O/COqjoD+DJw6WJOamKVJK1JVXVfVd1/mM3OBnZV1YNV9SRwE7A1SYBzgfcP290IXLiY83qNVZLUombzdpuTgYfG5ncDLwaeDnylqvaNLT95MQc0sUqSWjz8d3XLL3/yGxubD3tckh1j89uqatvBmSQfBZ41z35vrqoPLeL48z3cuA6x/LBMrJKkFlV1wRTOed6Eh9gNbBqbPwXYA+wFjk/ylKFqPbj8sLzGKkk6mt0JnDGMAN4AXARsr6oCPg68ZtjuEmAxFbCJVZK0NiX50SS7gX8MfDjJLcPyZye5GWCoRq8EbgHuA95XVfcMh3gj8PNJdjG65nr9os47SsqSJKmDFaskSY1MrJIkNTKxSpLUyMQqSVIjE6skSY1MrJIkNTKxSpLUyMQqSVKj/wcHSWngAAAABElEQVRFUYfuGOfKfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(covariances[-1,-1,:,:], cmap=cmap, interpolation='Nearest', vmin=-1, vmax=1)\n",
    "plt.colorbar()\n",
    "\n",
    "# show shifts in the labels\n",
    "plt.xticks(range(2*howFar + 1), range(-howFar,howFar+1))\n",
    "plt.yticks(range(2*howFar + 1), range(-howFar,howFar+1))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.99857504 -0.13869636 -0.06391662  0.00206449]\n",
      "[ 1.00000000e+00 -1.45608499e-01 -5.98825052e-02  8.19815336e-04\n",
      " -1.93509663e-01]\n"
     ]
    }
   ],
   "source": [
    "print(covariances[-1,-1,3:,3])\n",
    "print(truCov[-1,-5:][::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \"Bookkeeping kind of thing\"\n",
    "What Arun wanted to avoid but I don't know, I actually like to look at things that way (probably out of habit because of my OT cost matrices I suppose)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a plot with the coordinates of each pixel as x-axis label of the covariance matrix, and y-axis deltas for a given pixel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['0,0' '0,1' '0,2' '0,3' '0,4']\n",
      " ['1,0' '1,1' '1,2' '1,3' '1,4']\n",
      " ['2,0' '2,1' '2,2' '2,3' '2,4']\n",
      " ['3,0' '3,1' '3,2' '3,3' '3,4']\n",
      " ['4,0' '4,1' '4,2' '4,3' '4,4']]\n"
     ]
    }
   ],
   "source": [
    "pixCoords = []\n",
    "for i in range(dataSize[0]):\n",
    "    rowCoords = []\n",
    "    for j in range(dataSize[0]):\n",
    "        rowCoords += [str(i)+\",\"+str(j)]\n",
    "    pixCoords += [rowCoords]\n",
    "pixCoords = np.array(pixCoords)\n",
    "print(pixCoords)\n",
    "pixCoords = pixCoords.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, pick the pixel we'll look at:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixIdx = 2,3 #4,4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mapping from covariance entry to \"shift\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['2,3', '2,2', '2,1', '2,0', '2,-1', '1,3', '1,2', '1,1', '1,0', '1,-1', '0,3', '0,2', '0,1', '0,0', '0,-1', '-1,3', '-1,2', '-1,1', '-1,0', '-1,-1', '-2,3', '-2,2', '-2,1', '-2,0', '-2,-1']\n"
     ]
    }
   ],
   "source": [
    "ylabs = [str(pixIdx[0] - int(coord[0])) +\n",
    "         str(\",\") +\n",
    "         str(pixIdx[1] - int(coord[-1]))\n",
    "        for coord in pixCoords]\n",
    "print(ylabs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqsAAAKkCAYAAAAnTXYkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xm83VV97//XOyFhlikIgYRBQTECgkSch1unaBHsrVWwKlpsHKDVi729VvtTq+Xaaq+1rfjTXKViq+BUe9FiKVe0VguUiFEmKSmCxIABgkyBkOFz/9g7dnM4GfbZ65yzT87r+XjsR873+137/V1nJ4SVz1lrfVNVSJIkScNoxmR3QJIkSdocB6uSJEkaWg5WJUmSNLQcrEqSJGloOViVJEnS0HKwKkmSpKHlYFWSJEm/lOScJKuSXL2Z60nyl0mWJ/lRkif3XDs1yQ3d16kt+uNgVZIkSb0+AyzawvWXAId3X4uB/x8gyd7Ae4GnAscD702y16CdcbAqSZKkX6qq7wCrt9DkJOCz1XEZsGeSucCLgYuranVV3QVczJYHvdtkh0EDJEmStG12z6NrAw9Nah8e4O5rgAd7Ti2pqiV9RBwI3NJzvKJ7bnPnB+JgVZIkaYJs4CEO41mT2oer+IcHq2rhABEZ5Vxt4fxAnAYgSZKkfqwA5vcczwNWbuH8QBysSpIkqR8XAK/r7grwNODuqroVuAh4UZK9ugurXtQ9NxCnAUiSJE2gGTMmuVa4ccuXk5wHPA+Yk2QFnRX+swCq6hPAhcBLgeXAGuAN3Wurk3wAuKIb9f6q2tJCrW2SqoGnEkiSJGkb7JI963EznjOpffjhxq99f8A5qxPKyqokSdKECcnMye7ElOKcVUmSJA0tB6uSJEkaWk4DkCRJmiDJECyw2jC5t++XlVVJkiQNLSurkiRJE8YFVv2ysipJkqSh5WBVkiRJQ8tpAJIkSRNo0hdYTTF+WpIkSRpaVlYlSZImTEisFfbDT0uSJElDy8GqJEmShpbTACRJkiZI5wlW7rPaDyurkiRJGlpWViVJkiaMC6z65aclSZKkoeVgVZIkSUPLaQCSJEkTyAVW/bGyKkmSpKFlZVWSJGmCJC6w6pefliRJkoaWg1VJkiQNLacBSJIkTaAZM6wV9sNPS5IkSUPLwaokSZKGltMAJEmSJkxI3Ge1H1ZWJUmSNLSsrEqSJE2QxAVW/fLTkiRJ0tBysCpJkqSh5TQASZKkCeTjVvvjpyVJkqShZWVVkiRpwoQZM9y6qh9WViVJkjS0HKxKkiRpaDkNQJIkaQK5wKo/flqSJEkaWg5WJU05SZ6XZMUWrn8iyf/Xc/yWJD9Pcl+SfcahP/clecyAGYckqSSb/YlXkg8mefsg99nGvpyY5Pzxvo80HSWdBVaT+ZpqHKxKEyTJq5Ms7Q5sbk3yjSTPmux+bbK1AeAYM9+X5G9bZm6LqnpzVX2g24dZwEeAF1XVblV1Z3dQeFjD++1WVTe2yhtNkn2B1wGf7B4/LcnFSVYnuT3Jl5LM7SPvW9333ZPkh0lO2nStqi4AjkxydPNvRJL65GBVmgBJzgQ+CvxPYD/gIODjwElbet+w2VLVb4jtB+wEXDPZHRnQ64ELq+qB7vFewBLgEOBg4F7gr/vIexswt6oeBSwG/nbEYPe87nlJmlQOVqVxlmQP4P3A6VX1d1V1f1Wtq6qvVdV/77bZMclHk6zsvj6aZMfuteuSnNCTt0OSO5I8eTP3OyHJsiS/SPKvvdWxJDcl+b0kP0pyd5IvJNkpya7AN4ADupXf+5Ic0K2MfjnJ3ya5B3hnkjW9P0pPcly3QjdrRD8WAe8CXtXN+2H3/AFJLuhWBJcn+e0tfHYvTXJtknuT/CzJ7424/o4kq7qV6jf0nP9Mkj9O8jjg+u7pXyS5JMl3usc/7PbrVZu7/4i8T3Qrmfcm+eckB/dcrySHJZnd/ex/p3t+ZpLvJXlP93hGkncm+Y8kdyb5YpK9t3b/rpcA/7zpoKq+UVVfqqp7qmoN8DHgmduYRVX9qKrWbzoEZgHze5p8G/jVbc2TtO1mzJgxqa+pZur1WJp6nk6nsvfVLbR5N/A04BjgScDxwB92r50HnNLT9sXAHVV15ciQ7gD2HOBNwD50fmR8waaBb9crgUXAocDRwOur6n46g6GV3R9p71ZVK7vtTwK+DOwJ/C86g5hX9uS9Bji/qtb19qWq/pFOJfkL3bwn9Xw/K4ADgFcA/zPJ8zfzuXwaeFNV7Q4cCVzSc21/YA/gQOA04Owke43ow78DT+we7llVv1JVz+keP6nbry9s5t4j/SbwAWAOsAz43MgGVfUQnc/j/UmeALwTmAmc1W3yu8DLged2v/+7gLO38f5H8Z8D79E8hz6rx0m+nuRB4HI6v69Ley5fBxyS5FH9ZEpSaw5WpfG3D53B5fottPlN4P1Vtaqqbgf+CHht99rngROT7NI9fnX33Gh+G/hkVV1eVRuq6lxgLZ2B8CZ/WVUrq2o18DU6A+QtubSq/r6qNnZ/BH0unQEZSWbSGUj/zVYy6LafDzwL+B9V9WBVLQM+1fO9jrQOWJDkUVV114gB+jo6n9m6qroQuA94/Lb0Y4z+oaq+U1Vr6fzj4und7+dhqupq4I/p/OPk94DXVtWG7uU3Ae+uqhXdnPcBr9jG6RV70vlR/yN0q+fvAf57P99QVZ0A7A68FLioqjb2XN50rz37yZS0NSGZMamvqWbq9Viaeu4E5mxlQHIAcHPP8c3dc1TVcjpVrpd1B6wnsvnB6sHAO7pTAH6R5Bd0frR7QE+b23q+XgPstpX+3zLi+P/QGUA+BnghcHdV/dtWMjY5AFhdVb2DrpvpVEdH8+t0BlI3d3/0/vSea3eO+AfAtnwvg/jl51BV9wGrefjn2utcOnNJL6yqG3rOHwx8tef35jpgA515tVtzF52B5cN0F4p9A3hbVf3LNuQ8THew/w3gxUlO7Lm06V6/6DdTklpysCqNv0uBB+n8+HdzVtIZyGxyUPfcJpumApwEXNsdwI7mFuCsqtqz57VLVZ23Df2sbTlfVQ8CX6RTDX4tW66qjsxcCeydpHfQdRDws1HfXHVFVZ0EPBr4++59J8svq6hJdgP25uG/R70+DnydzgCwd8eHW4CXjPj92amqRv3+R/gR8LjeE915s/8X+EBVbVN1ewt2AB7bc/wE4KaqumfAXEkaiINVaZxV1d10fkR7dpKXJ9klyawkL0nyoW6z84A/TLJvkjnd9r1bPp0PvAh4C5uvqgL8b+DNSZ6ajl2T/OqIweHm/BzYJ50FYVvzWTqr008c0c/RMg9J9+dOVXUL8K/AB7sLu46mM9/0EfM/u4uVfjPJHt35sPfQqUK28HPgYfuidhdJPW8L73lpkmclmU1n7url3e9nZL9fCxxH5/P5XeDc7uAW4BPAWZsWZ3V/v7d1R4gL6cx13XSfA+nM4T27qj4xSj9en+Sm0YKSHNH987dz98/ia+jMef3nnmbPpVOxldRQgvus9snBqjQBquojwJl0Fk3dTqfCdgadaiF05jgupVM9uwq4sntu0/tvpVOhfQaw2QVBVbWUzrzVj9H5sfFyOoOmbenjj+kMmm/s/ph6cz/ipqq+B2wErqyqm7YQ+6Xur3cm2TTf9BQ6PyJfSWde53ur6uLNvP+1wE3dnQjeTHeubAPvozOI/EWSVyaZR2fO61VbeM/ngffS+fH/cXQqyw+T5CA6W5S9rqruq6rP0/l9/fNuk78ALgD+Kcm9wGXAU7exz5+lM2DeuXv8RjoD7vfmP3dwuK+n/Xzge5vJCp3PYBWdP49vA141Yk7wKXT3dJWkyZSqzf3kT5I2L8klwOer6lOT3ZdBdSuLT6yqP9jM9c8AK6rqD0e7PlGS/E9gVVV9dBva/hOdeazXjeE+L6OzMOyVW20sqS977nhAPeeAN05qH7520we+X1ULJ7UTfZiKG3xLmmRJngI8mSn2UIPNqaoJf8rWWFTVu/po+6IB7vM1OjtFSNKkcxqApL4kOZfOop63j1jVL0lSc1ZWJfWlqk6d7D5MtKp6/WT3QdL2YyoucppMVlYlSZI0tIaysrpDZtdsdtl6wz484ZhHPGimjW168Eyfkevvb54JwLo1zSPvZ1sfa96fXXYeh391blzbPPLBdbOaZwLsnPZbW967bnz2y58xI80zd9mpeSQAefDO5pmr7t21eSbAvnPa/92yMe0/2Bls6cFsY/fgQ+3/XO00e+PWG41FtdpR7T9l3Tj9f2Dmjltv06+H7tt6mz5t3GX/5pkAP/jBD+6oqn3HJVzjZigHq7PZhcN41tYb9uHSS/6qad4mG2e1fxLh7Du39WFAfVq1rHnk0g2nbL3RGBx95LZsC9qfHdbcvPVGfbrm1m158FD/jpp1UfPMf7m17X9Tm+y86+zmmccc0TwSgB2ub7+O6uPffkrzTIA3v6H9/0/vm31Y88xdMj4PuLrupvb/YH3CQQ82zwSYse7u9pm3Xd48E4A9H7v1Nv1a0feD07Zq7bHvbJ4JsMsuu7T/H0GfkkzJR55OJj8tSZIkDa2hrKxKkiRtr2bMsFbYDz8tSZIkDS0Hq5IkSRpaTgOQJEmaMCFxn9V+DFRZTTI/ybeSXJfkmiRvG6XNm5NclWRZku8mWTDIPSVJkjR9DFpZXQ+8o6quTLI78P0kF1fVtT1tPl9VnwBIciLwEWDRgPeVJEmachIXWPVroE+rqm6tqiu7X98LXAccOKJN7+7muwI1yD0lSZI0fTSbs5rkEOBY4BE7GSc5HTgTmA38Sqt7SpIkafvWZLCaZDfgK8DbR1RSAaiqs4Gzk7wa+EPg1FEyFgOLAWaxc4tuSZIkDZ0ZM9s/Tnh7NvCkiSSz6AxUP1dVf7eV5ucDLx/tQlUtqaqFVbVwB9o/vlGSJElTz0CV1SQBPg1cV1Uf6Tl/BkBVfSzJ4VV1Q/fSrwI3PDJJkiRpGkjITBdY9WPQaQDPBF4LXJVkWffcu4AjgO91j89I8gJgHXAXo0wBkCRJkkYz0GC1qr4LPGLiRZK30llQRVU9Yu9VSZIkaVuMyxOsquqE8ciVJEmayoILrPrlpAlJkiQNrXGprEqSJGkUsbLaLyurkiRJGlpDWVl9wjHzufSSv2qauXDv32mat8lXbzy3eeaeezyzeSbAz9cc2zzzyXvf0jwT4Jr/mNU887GHPW4cMptHdlz78+aRT33GAc0zAW5esaZ9aN3bPhOoec9unvnG049sngkw464fNM+8/d4Hm2fedP19zTMBnrSw/Z/XHX5yfvNMgFUHntY+dL/HtM8EHn3vxc0z6+AXNM/U5EuyCPgLYCbwqar6kxHX/xz4L93DXYBHV9We3WsbgKu6135aVScO0pehHKxKkiRtl6bAPqtJZgJnAy8EVgBXJLmgqq7d1Kaq/ltP+98BeitiD1TVMa36M9yfliRJkiba8cDyqrqxqh6i8wTSk7bQ/hTgvPHqjINVSZKk6WVOkqU9r8Ujrh8I9M7zW9E99whJDgYOBS7pOb1TN/eyJC8ftLNOA5AkSZogQ7LP6h1VtXAL10frYG2m7cnAl6tqQ8+5g6pqZZLHAJckuaqq/mOsnbWyKkmSpF4rgPk9x/OAlZtpezIjpgBU1crurzcC3+bh81n7ZmVVkiRpogRmDPkCK+AK4PAkhwI/ozMgffXIRkkeD+wFXNpzbi9gTVWtTTIHeCbwoUE642BVkiRJv1RV65OcAVxEZ+uqc6rqmiTvB5ZW1QXdpqcA51dV7xSBJwCfTLKRzk/w/6R3F4GxcLAqSZKkh6mqC4ELR5x7z4jj943yvn8FjmrZl4EGq0nmA58F9gc2Akuq6i9GtDkTeCOwHrgd+K2qunmQ+0qSJE1FQ7LAakoZdNLEeuAdVfUE4GnA6UkWjGjzA2BhVR0NfJkB5y1IkiRp+hioslpVtwK3dr++N8l1dPbh6n3Cwbd63nIZ8JpB7ilJkjRlJcTKal+aLUdLcgidrQku30Kz04BvbOb9izdtTnvHHatbdUuSJElTWJPBapLdgK8Ab6+qezbT5jXAQuDDo12vqiVVtbCqFs6Zs3eLbkmSJGmKG3g3gCSz6AxUP1dVf7eZNi8A3g08t6rWDnpPSZKkKWlq7LM6VAbdDSDAp4HrquojPefPAKiqjyU5FvgksKiqVg1yP0mSJE0vg1ZWnwm8FrgqybLuuXcBRwDf6x5/GNgN+FJnbMtPq+rEAe8rSZI05bh1Vf8G3Q3gu3Q+94dJ8lbgzG6bFwxyD0mSJE1f4/IEq6o6YTxyJUmSNL34uFVJkqSJkrjAqk9+WpIkSRpaVlYlSZImkE+w6s9wDlazAxtn7dk08qs3nts0b5Nfe8ypzTO/f8dHtt5oDHbfb2bzzPtnHtE8E+CoNf+reea1N72peeZhj9m9eSbA6se8tXnmmjvHZ4vjffbesXnmbfe3/7MKsNuu+zbPXL9mY/NMgDvuO7J55gFzd2qeeehx9zfPBPj3Ox9snnn/3u3/vgaYsXZD88x77h6n/173H4e/s2fMap8p9XAagCRJkobWcFZWJUmStkOJ+6z2y8qqJEmShpaVVUmSpIkS3LqqT35akiRJGloOViVJkjS0nAYgSZI0YeI+q30aqLKaZH6SbyW5Lsk1Sd42SpvnJLkyyfokrxjkfpIkSZpeBp0GsB54R1U9AXgacHqSBSPa/BR4PfD5Ae8lSZKkaWagaQBVdStwa/fre5NcBxwIXNvT5iaAJOPzmBdJkqQpIu4G0Ldmn1aSQ4BjgcvH+P7FSZYmWXr7HXe06pYkSZKmsCYLrJLsBnwFeHtV3TOWjKpaAiwBOO7Jx1aLfkmSJA0Vn2DVt4Erq0lm0Rmofq6q/m7wLkmSJEkdg+4GEODTwHVV9ZGe82ckOWPQzkmSJGl6G3QawDOB1wJXJVnWPfcu4AjgewBJngJ8FdgLeFmSP6qqJw54X0mSpCkoLrDq06C7AXwXeMTEiyRvBc7strkCmDfIfSRJkjQ9jcsTrKrqhPHIlSRJmsoSfIJVn6xDS5IkaWg5WJUkSdLQGpdpAJIkSRqd+6z2ZygHq1l/P7Pv/LemmXvu8cymeZt8/46PbL1Rn46bc2bzTIBl/zi3eebPDjqreSbAYY//zeaZa/59bfPMX9yzc/NMgOuvXtU889lz/ql5JsD3176seeahj92zeSbATT+5u3nmgsfv1jwTYM4P39M+dLdXNo+8c+enNM8EOGJW+z+vP37g2c0zAX502U3NM1/66+Ozac4NP2v/A9VZs9tnztuleaSmsKEcrEqSJG2X4tZV/fLTkiRJ0tBysCpJkqSh5TQASZKkCZK4wKpfVlYlSZI0tKysSpIkTaC4wKovflqSJEkaWg5WJUmSNLQGmgaQZD7wWWB/YCOwpKr+YkSbHbttjgPuBF5VVTcNcl9JkqQpyQVWfRu0sroeeEdVPQF4GnB6kgUj2pwG3FVVhwF/DvzpgPeUJEnSNDFQZbWqbgVu7X59b5LrgAOBa3uanQS8r/v1l4GPJUlV1SD3liRJmmpCmDHDymo/ms1ZTXIIcCxw+YhLBwK3AFTVeuBuYJ9R3r84ydIkS29f3f753ZIkSZp6mgxWk+wGfAV4e1XdM/LyKG95RFW1qpZU1cKqWrjv3nu06JYkSZKmuIH3WU0yi85A9XNV9XejNFkBzAdWJNkB2ANYPeh9JUmSppzAjBluxtSPgT6tJAE+DVxXVR/pOX9GkjO6hxcAp3a/fgVwifNVJUmStC0Graw+E3gtcFWSZd1z7wKOAL7XPf408DdJltOpqJ484D0lSZI0TQy6G8B3GWVOapK3Amd22zwI/MYg95EkSdoeBPdZ7dfAc1ZHU1UnjEeuJEmSppdxGaxKkiRpFIkLrPrkpyVJkqSh5WBVkiRJQ2s4pwGsWwOrlm29XR9+vubYpnmb7L7fzOaZy/5xbvNMgGMW3do8c9lln2qeCfCN605snrno+BXNM/PQhuaZAHsfe3DzzLU7v6F5JsBxP/1q88xV68dn2vsx685tH3rFHe0zgb+//U3NM5+y4NDmmXN2HZ//jdw148XNM4+44S+bZwIc8byjmmfes3FB80yA+fN3a565du34/D24PfNxq/2xsipJkqShNZyVVUmSpO1Q8AlW/fLTkiRJ0tBysCpJkqSh5TQASZKkiRIXWPXLyqokSZKGlpVVSZKkCRLCjJnWCvvR96eVZH6SbyW5Lsk1Sd62De85K8ktSe4bWzclSZI0HY1laL8eeEdVPQF4GnB6kq3tXvw14Pgx3EuSJEnTWN/TAKrqVuDW7tf3JrkOOBC4dgvvuQwgcUKxJEmaxlxg1beBJk0kOQQ4Frh80I4kWZxkaZKlt//i/kHjJEmStB0Y8wKrJLsBXwHeXlX3DNqRqloCLAFY+IR5NWieJEnSMPIJVv0Z06eVZBadgernqurvRlybmWRZ9/X+Fp2UJEnS9NR3ZTWdiaefBq6rqo/0nD8DoKo+BhzTrIeSJEmatsZSWX0m8FrgV3oqqC8FjgDuHO0NST6UZAWwS5IVSd435h5LkiRNUUmYMWNyX1NN34PVqvpuVaWqjq6qY7qvC4FD6EwNGO09v19V86pqRvfX9w3WbUmSJI2XJIuSXJ9keZJ3jnL99Ulu7ylcvrHn2qlJbui+Th20L82eYFVVJ7TKkiRJ0uRIMhM4G3ghsAK4IskFVTVym9IvVNUZI967N/BeYCFQwPe7771rrP3xcauSJEkTaArsBnA8sLyqbgRIcj5wElvYU7/Hi4GLq2p1970XA4uA88bamaH/tCRJkjShDgRu6Tle0T030q8n+VGSLyeZ3+d7t5mVVUmSpIkSyMxJX+Q0J8nSnuMl3f3uNxmtgyP3wP8acF5VrU3yZuBc4Fe28b19GcrB6v3szdINpzTNfPLet2y90RjcP/OI5pk/O+is5pkAyy77VPPMY572r80zAZbdfHzzzHV7tN9R7Zbbx+cvnNlr2v/Q44CHrm6eCXD57c9onnnUvPH5q+nBo/5b88xbfjY+T9x77tE7Nc9c88D65pk7rfp280yA61c/qXnmrbud1jwT4NG77tw8c6/1K5pnAtz+wNzmmXN2Wt08cwP7N8/Uw9xRVQu3cH0FML/neB6wsrdBVfXuAPW/gT/tee/zRrz322PtKDgNQJIkSQ93BXB4kkOTzAZOBi7obZCk918+JwLXdb++CHhRkr2S7AW8qHtuzIaysipJkrQ9Chn6BVZVtb77sKeLgJnAOVV1TffJpEur6gLgd5OcCKwHVgOv7753dZIP0BnwArx/02KrsXKwKkmSpIfp7qF/4Yhz7+n5+g+AP9jMe88BzmnVFwerkiRJE2gqPkVqMg13HVqSJEnTmoNVSZIkDS2nAUiSJE2QZEo8wWqoDPxpJTknyaoko27imOTNSa5KsizJd5MsGPSekiRJmh5aDO0/Q+eZr5vz+ao6qqqOAT4EfKTBPSVJkqakGTMyqa+pZuDBalV9h87+Wpu7fk/P4a4M+MgtSZIkTR8TMmc1yenAmcBsOs+NlSRJkrZqQgarVXU2cHaSVwN/CJw6sk2SxcBigP33P3AiuiVJkjShkjBjpgus+jHRn9b5wMtHu1BVS6pqYVUt3HPPvSe4W5IkSRpG41JZ7T5Plqr6WJLDq+qG7qVfBW7Y/DslSZK2b1NxkdNkGniwmuQ84HnAnCQrgPcCRwDf6zY5I8kLgHXAXYwyBUCSJEkazcCD1ao6ZeS5JF+ns6CKqnrboPeQJEnS9DQu0wCq6oTxyJUkSZrSfIJV3/y0JEmSNLQmZOsqSZIkQZiaT5GaTFZWJUmSNLQcrEqSJGloDeU0gF12nsnRR+7eNPOa/5jVNG+To9b8r+aZhz3+N5tnAnzjuhObZy67+fjmmQDHHPzR5plL7/p488yD93iweSbAzHuub565bo+nN88EOPTw2c0zZ228t3kmwA73/aR55mG7bmyeCVDrd2meed+Gg5tnLl//lOaZAEfveWn70LuWt88EWNX2/1cAbNzQPhPY+7D2/3+ZedeK5pkbdtq/eebQcIFV3/y0JEmSNLQcrEqSJGloDeU0AEmSpO1RgBkz3Q2gH1ZWJUmSNLSsrEqSJE2UxAVWffLTkiRJ0tBysCpJkqShNfA0gCTnACcAq6rqyFGunwm8EVgP3A78VlXdPOh9JUmSpiIft9qfFpXVzwCLtnD9B8DCqjoa+DLwoQb3lCRJ0jQwcGW1qr6T5JAtXP9Wz+FlwGsGvackSdJUFKys9mui56yeBnxjtAtJFidZmmTp7XfcMcHdkiRJ0jCasMFqktcAC4EPj3a9qpZU1cKqWrjvnDkT1S1JkiQNsQnZZzXJC4B3A8+tqrUTcU9JkqShE4jTAPoyLoPVJGcAVNXHkhwLfBJYVFWrxuN+kiRJ2j612LrqPOB5wJwkK4D3AkcA3+s2+TCwG/ClJAA/raoTB72vJEnSVOQCq/602A3glJHnknwdOLN7/QWD3kOSJEnT07hMA6iqE8YjV5IkSdPLhCywkiRJEoTQnRapbTTR+6xKkiRJ28zKqiRJ0kSJC6z6NZyD1Y1r2WHNzU0jH3vY45rmbXLtTW9qnrnm38dnK9pFx69onrluj2OaZwIsvevjzTMX7vXW5plLV/9V80yAtfu3X5f4wIMbmmcCbNiwvnlmZXbzTICLLt+9eeYznndw80yAR111VvPMuUe9pXnmA9mneSYAt/y8eeRDj3t980yAi/+x/d+tL/jVw5pnAsxY94vmmbflqOaZezdP1FTmNABJkiQNreGsrEqSJG2nfIJVf6ysSpIkaWhZWZUkSZogcYFV36ysSpIkaWg5WJUkSdLQchqAJEnSBHKBVX8GrqwmOSfJqiRXb+b6c5JcmWR9klcMej9JkiRNHy2mAXwGWLSF6z8FXg98vsG9JEmSNI0MPA2gqr6T5JAtXL8JIMnGQe8lSZI0tcXdAPrkAitJkiQNraFZYJVkMbAY4KB5cye5N5IkSeMgLrDq19BUVqtqSVUtrKqFc/bZc7K7I0mSpCEwLoPVJGckOWM8siVJkjR9tNi66jzgUuDxSVYkOQ04Arize/0pSVYAvwF8Msk1g95TkiRpKgowI5nU11QJs5dQAAAgAElEQVTTYjeAU0aeS/J14Mzu9SuAeYPeR5IkSdPPuCywqqoTxiNXkiRpqnOBVX+GZoGVJEmSNJKDVUmSJA2todlnVZIkabsXfIJVn4ZysPrgullcc+t+TTMfe1jTuF867DG7N8/8xT07N88EyEMbmmfecvv4/Ad38B4PNs9cuvqvmmcu3Pt3mmcCfOHHf90881G7z2qeCbCx2mfOXLOifSiwaMG1zTNve2B81o/udOzvN8+c/fN/bp45a+6Lm2cC3DPv5OaZK26+v3kmwAsXtf8zsPah8XlC+QMP7to8c+7Pz22euXavNzbP1NQ1lINVSZKk7VGIC6z65JxVSZIkDS0Hq5IkSRpaTgOQJEmaKC6w6puVVUmSJA0tK6uSJEkTyAVW/bGyKkmSpKHlYFWSJElDa+DBapJzkqxKcvVmru+Y5AtJlie5PMkhg95TkiRpKgowI5nU11TTorL6GWDRFq6fBtxVVYcBfw78aYN7SpIkaRoYeLBaVd8BVm+hyUnApmexfRl4fjIFh/WSJEmacBOxG8CBwC0AVbU+yd3APsAdvY2SLAYWA8yde+AEdEuSJGmCxd0A+jURC6xG+x2pR5yoWlJVC6tq4V577TMB3ZIkSdJokixKcn13zdE7R7l+ZpJrk/woyTeTHNxzbUOSZd3XBYP2ZSIqqyuA+cCKJDsAe7DlaQOSJEnbqQz9E6ySzATOBl5IZxx3RZILquranmY/ABZW1ZokbwE+BLyqe+2BqjqmVX/GpbKa5IwkZ3QPLwBO7X79CuCSqnpEZVWSJElD4XhgeVXdWFUPAefTWYP0S1X1rapa0z28DJg3Xp1psXXVecClwOOTrEhyGnAEcGe3yaeBfZIsB84EHlFKliRJ0oSZk2Rpz2vxiOu/XG/UtaJ7bnNOA77Rc7xTN/eyJC8ftLMDTwOoqlNGnkvydToDU6rqQeA3Br2PJEnS9mAIFljdUVULt3B9m9YbASR5DbAQeG7P6YOqamWSxwCXJLmqqv5jrJ0dl2kAVXVCt2wsSZKkqWXTeqNN5gErRzZK8gLg3cCJVbV20/mqWtn99Ubg28Cxg3RmIhZYSZIkCUgY+gVWwBXA4UkOBX4GnAy8urdBkmOBTwKLqmpVz/m9gDVVtTbJHOCZdBZfjZmDVUmSJP1Sd1/8M4CLgJnAOVV1TZL3A0ur6gLgw8BuwJe6z3r6aVWdCDwB+GSSjXR+gv8nI3YR6JuDVUmSJD1MVV0IXDji3Ht6vn7BZt73r8BRLfsylIPVnXMPR826qG3otT9vm9e1+jFvbZ55/dWrtt5oDPY+9uCtN+rT7DXj81yJmfdc3zxz7f6j/nc1kC/8+K+bZwK86og3NM9cds0rmmcCrDn4lc0zb1w5t3kmwH6PPrR55v53NP67quva1c9qnrlgl5nNM3e69uzmmQBrDn1T88zbV97TPBPg0fvu1zzz5pvGZzvy8fgM5j+2/Rrqw5onDpchWGA1pUzEE6wkSZKkMRnKyqokSdL2akasrPbDyqokSZKGloNVSZIkDS2nAUiSJE2QJC6w6pOVVUmSJA0tK6uSJEkTaAo8wWqo9F1ZTXJOklVJrt7G9mcluSXJff13T5IkSdPZWKYBfAZY1Ef7rwHHj+E+kiRJmub6ngZQVd9Jckgf7S+DzoRiSZKkaS0+wapfQ7PAKsniJEuTLL199b2T3R1JkiQNgaEZrFbVkqpaWFUL991798nujiRJ0rhIJvc11TQdrCaZmWRZ9/X+ltmSJEmafppsXZXkDICq+hhwTItMSZIkaSxbV50HXAo8PsmKJKcBRwB3bqb9h5KsAHbptn/fIB2WJEmaqkL3KVaT+JpqxrIbwCkjzyX5OnDmZtr/PvD7/XdNkiRJ012TaQBVdUKLHEmSJKmXj1uVJEmaQG6z2p+h2bpKkiRJGsnKqiRJ0kTxCVZ9G8rB6r3rduNfbn1W08ynPuOApnmbrLlzbfPMZ8/5p+aZAGt3fkPzzAMeurp5JsC6PZ7ePPOBBzc0z3zU7rOaZwIsu+YVzTOPeeKXm2cCfPbKFzXPfNzytzTPBNjx6aOuAx3IZy/er3kmwCmvGYfQle2fDrhx/n9pngmw95ormmc+d+dLmmcCcM+xzSPn3PWPzTMBmHto88j1B72ueWb7v601lTkNQJIkSUNrKCurkiRJ26spuNXppLKyKkmSpKFlZVWSJGmChKn5FKnJZGVVkiRJQ8vBqiRJkoaW0wAkSZImkNus9sfKqiRJkobWwIPVJIuSXJ9keZJ3jnL9A0l+lGRZkn9KMj6780uSJE0BSSb1NdUMNFhNMhM4G3gJsAA4JcmCEc0+XFVHV9UxwNeB9wxyT0mSJE0fg1ZWjweWV9WNVfUQcD5wUm+Dqrqn53BXoAa8pyRJkqaJQRdYHQjc0nO8AnjqyEZJzgJeB9wNjPog6SSLgcUA+z3amQKSJGk7FJ9g1a9BK6ujfdyPqJxW1buraj7wOeCM0YKqaklVLayqhXvsufeA3ZIkSdL2YNDB6gpgfs/xPGDlFtp/Hvj1Ae8pSZI0JQUXWPVr0MHqFcDhSQ5NMhs4GbggyQeT/BpAksN72p8I/HjAe0qSJGmaGGjOalWtT3IGcBEwEzinqq5JchRwQbfZnyR5PLARuBl48yD3lCRJ0vQx8BOsqupC4MIRp2dV1aXd6/7YX5IkqSs+kqkv4/JxVdWLxyNXkiRJ08vAlVVJkiRtuxlTcJHTZLIQLUmSpKHlYFWSJElDayinAcyYEXbedXbTzJtXrGmat8k+e+/YPPP7a1/WPBPguJ9+tXnm5bc/o3kmwKGHt/39B9iwYX3zzI3j9PDgNQe/snnmZ698UfNMgNc9eXHzzCtvfW/zTIAHdn9i88wDD/5p80yAWSu/2Tzzzn3a/92ybt3G5pkAO+zY/sekc+bes/VGY3DzzOc0zzx49x81zwR46MjTm2fOumtZ88wN++zTPHNo+ASrvllZlSRJ0tBysCpJkqShNZTTACRJkrZHYWo+8nQyWVmVJEnS0LKyKkmSNIEsrPbHyqokSZKGloNVSZIkDa2BB6tJFiW5PsnyJO8c5fqHk/w4yY+SfDXJnoPeU5IkaapKMqmvqWagwWqSmcDZwEuABcApSRaMaHYxcGRVHQ38O/AHg9xTkiRJ08egC6yOB5ZX1Y0ASc4HTgKu3dSgqv6pp/1lwCsGvKckSdKUNWPqFTcn1aDTAA4Ebuk5XtE9tzm/BXxjtAtJFidZmmTpL36xesBuSZIkaXsw6GB1tH8bjPq09CTvBtYDnxvtelUtqaqFVbVwzz33HrBbkiRJ2h4MOg1gBTC/53gesHJkoySnAicAz6+qUQezkiRJ27vEfVb7NWhl9Qrg8CSHJpkNnAxckOSDSX4NOrsFAP8DOLGq1gx4P0mSJE0jA1VWq2p9kjOAi4CZwDlVdU2So4ALus0+BuwIXNzdLuGyqnrzIPeVJEmaquIKq74M/LjVqroQuHDE6VlVdWn3+mGD3kOSJEnT07g8waqqXjweuZIkSZpeBq6sSpIkadu5wKo/41JZlSRJklqwsipJkjRBEphhabUvQzlY3WUnOOaIxqF1b+PAjtvun9k889DH7tk8E2DV+hOaZx41b3z+CM3a2P73qzK7eebMNSuaZwLcuHJu88zHLX9L80yAK299b/PMJ8/9o+aZAN+794vNM5+/V/tMgH9Z8ermmY/ZbUPzzAMf+GbzTIBr1z6reebSGw5vngnw+Ce13z783/it5pkAs6+/r3nmgiOOa54p9XIagCRJkobWUFZWJUmStlfOAuiPlVVJkiQNLSurkiRJEyiWVvtiZVWSJElDy8GqJEmShpbTACRJkiaQswD6M3BlNcmiJNcnWZ7knaNc/40k1yTZmGThoPeTJEnS+NqG8d2OSb7QvX55kkN6rv1B9/z1SV48aF8GGqwmmQmcDbwEWACckmTBiGZXA/8V+M4g95IkSdL428bx3WnAXVV1GPDnwJ9237sAOBl4IrAI+Hg3b8wGraweDyyvqhur6iHgfOCk3gZVdV1VXT/gfSRJkrYDIZnc1zbY6viue3xu9+svA89PJ/wk4PyqWltVPwGWd/PGbNDB6oHALT3HK7rnJEmSNDVty/jul22qaj1wN7DPNr63L4MusBpteD6mhyQnWQwsBjhonuNdSZK0/UmKZExDpZbmJFnac7ykqpb0HG/L+G5zbZqNDTcZdLC6ApjfczwPWDmWoO6HtATguGOPnvTfRUmSpO3UHVW1pUXv2zK+29RmRZIdgD2A1dv43r4MOg3gCuDwJIcmmU1nQu0FST6Y5NcGzJYkSdLEG3V8N6LNBcCp3a9fAVxSVdU9f3J3t4BDgcOBfxukMwNVVqtqfZIzgIuAmcA5VXVNkqO6naU7aP0rYF/gH5Isq6qBtzGQJEmacgqoDZPdiy3awvju/cDSqroA+DTwN0mW06montx97zVJvghcC6wHTq8a7Bse+KEAVXUhcOGI07Oq6tLu9a8CXx30PpIkSZoYo43vquo9PV8/CPzGZt57FnBWq76MyxOsrJxKkiSNpkhtnOxOTCkDP8FKkiRJGi8OViVJkjS0xmUagCRJkjZjyBdYDRsrq5IkSRpaQ1lZzYN3ssP1f9s0s+Y9u2neJrvtum/zzJt+cnfzTIBj1p279UZ9evCo/9Y8E2CH+37SPPOiy3dvnrlowbXNMwH2e/ShzTN3fPqZzTMBHtj9ic0zv3fvF5tnAjxz91c2z7zytj9qngnw7Fv+vnnmQ3uc1jxzw05HNc8EWHDle5tnHvDcP2ueCXDbbWuaZx532D3NMwFuvu+A5pn/fMktW2/Up+c8/5DmmcOjwAVWfbGyKkmSpKHlYFWSJElDayinAUiSJG2v4gKrvlhZlSRJ0tCysipJkjRRqty6qk9WViVJkjS0HKxKkiRpaA08WE2yKMn1SZYneeco1/dOcnGSG7q/7jXoPSVJkqas2ji5rylmoMFqkpnA2cBLgAXAKUkWjGj2TuCbVXU48M3usSRJkrRVg1ZWjweWV9WNVfUQcD5w0og2JwGbHp10LvDyAe8pSZKkaWLQ3QAOBHqfs7YCeOqINvtV1a0AVXVrkkePFpRkMbAY4KD99xiwW5IkScPI3QD6NWhlNaOcq7EEVdWSqlpYVQv33XPXAbslSZKk7cGgldUVwPye43nAyhFtfp5kbreqOhdYNeA9JUmSpqxMwUVOk2nQyuoVwOFJDk0yGzgZuCDJB5P8WrfNBcCp3a9PBf7PgPeUJEnSNDHQYLWq1gNnABcB1wFfrKprgKOA27rN/gR4YZIbgBd2jyVJkqStGvhxq1V1IXDhiNOzqurS7vU7gecPeh9JkqSprwAXWPVjXJ5gVVUvHo9cSZIkTS8DV1YlSZK0jYop+RSpyTQulVVJkiSpBQerkiRJGlpDOQ1g1b278vFvP6Vp5htPP7Jp3ibr17Qv5S94/G7NMwG44o7mkbf87P7mmQCH7dr+c33G8w5unnnbA/OaZwLsf8dFzTM/e/F+zTMBDjz4p80zn7/XF5tnAlx52x81z3zy/u9tngmw7MpxmPo/c8fmket2nr/1RmOw8antP9c97/pe80yA+x818sGNg9vhgRXNMwEes/ILzTMPfvrrmmeub544THyCVb+srEqSJGloDWVlVZIkaXvlE6z6Y2VVkiRJQ8vBqiRJkoaW0wAkSZImjAus+mVlVZIkSUPLyqokSdKEsbLarzFVVpMsSnJ9kuVJ3rkN7c9JsirJ1WO5nyRJkqanvgerSWYCZwMvARYApyRZsJW3fQZY1HfvJEmSNK2NZRrA8cDyqroRIMn5wEnAtZt7Q1V9J8khY+mgJEnS9iLlPqv9Gss0gAOBW3qOV3TPDSTJ4iRLkyy97/67B42TJEnSdmAsldWMcq4G7UhVLQGWABw073ED50mSJA0fF1j1ayyV1RXA/J7jecDKTQdJ5idZ1n29edAOSpIkafoaS2X1CuDwJIcCPwNOBl6d5IPAv1XVV4FjGvZRkiRJ01TfldWqWg+cAVwEXAd8saquAY4CbhvtPUnOAy4FHp9kRZLTxt5lSZKkKaw2Tu5rihnTQwGq6kLgwhGnZ1XVpZtpf8pY7iNJkqTprdnjVqvqxa2yJEmSJPBxq5IkSRPI3QD61ayyKkmSJLVmZVWSJGkCham3yGkyWVmVJEnS0BrKyuq+c3bgzW/Yt2nmjLt+0DRvkzvuO7J55pwfvqd5JsDf3/6m5pnPPXqn5pkAtX6X5pmPuuqs5pk7Hfv7zTMBrl39rOaZp7ymeSQAs1Z+s3nmv6x4dfNMgGff8vfNM5ddOT5rS4958kXNMz/z/f/aPPPoW9/VPBPgwvvf0Tzz8Cc9uXkmwONX/2XzzLuOOLN5JsDyHZ7YPPOQDXs2z9y9eaKmsqEcrEqSJG2fXGDVL6cBSJIkaWhZWZUkSZooxZR8itRksrIqSZKkoeVgVZIkSUPLaQCSJEkTxgVW/bKyKkmSpKE15sFqkiOSXJpkbZLf20K7f0zywyTXJPlEkpljvackSdJUl9o4qa+pZpDK6mrgd4E/20q7V1bVk4AjgX2B3xjgnpIkSZpGxjxYrapVVXUFsG4r7e7pfrkDMJvOpg2SJEnSVk3IAqskFwHHA98AvryZNouBxQAHzZs7Ed2SJEmaYC6w6teELLCqqhcDc4EdgV/ZTJslVbWwqhbO2af9c4YlSZI09fRVWU1yOvDb3cOXVtXKbX1vVT2Y5ALgJODifu4rSZK0XShgo5XVfvRVWa2qs6vqmO5rswPVJN9McmCS3ZLM7Z7bAXgp8OPBuixJkqTpYsxzVpPsDywFHgVsTPJ2YAFwH3AYnd0CHgVckGRHYCZwCfCJQTstSZKk6WHMg9Wqug2YN/J8kiOBr1TVA8ADwFPG3j1JkqTtSUGtn+xOTCnNdwOoqquBM1vnSpIkafqZkK2rJEmSBFAusOrThGxdJUmSJI2Fg1VJkiQNraGcBrAxO3Hf7MOaZt5+74NN8zY5YO5O7UN3e2X7TOApCw5tnrnmgfGZJH7fhoObZ8496i3NM2f//J+bZwIs2GVm+9CV97bPBO7c52XNMx+z2/j8iOyhPU5rHzpzx/aZwGe+/1+bZ77+uDc1z/zazX/TPBNgxx/d1jzzsfves/VGY7BhzsnNM3fZaXxqSQuOmtM8c8f7/6N55nrajgGGjk+w6ouVVUmSJA0tB6uSJEkaWkM5DUCSJGm7VO4G0C8rq5IkSdomSfZOcnGSG7q/7jVKm2OSXJrkmiQ/SvKqnmufSfKTJMu6r2O2dk8Hq5IkSROm+wSryXwN5p3AN6vqcOCb3eOR1gCvq6onAouAjybZs+f6f6+qY7qvZVu7oYNVSZIkbauTgHO7X58LvHxkg6r696q6ofv1SmAVsO9Yb+hgVZIkaXqZk2Rpz2txH+/dr6puBej++ugtNU5yPDAb6N3j7Kzu9IA/T7LV/f/GvMAqyRHAXwNPBt5dVX82SptdgC8BjwU2AF+rqtHKxZIkSdPD5C+wuqOqFm7uYpL/C+w/yqV393OTJHOBvwFOraqN3dN/ANxGZwC7BPgfwPu3lDPIbgCrgd9llPLvCH9WVd9KMhv4ZpKXVNU3BrivJEmSxklVvWBz15L8PMncqrq1OxhdtZl2jwL+AfjDqrqsJ/vW7pdrk/w18Htb68+YpwFU1aqqugJYt4U2a6rqW92vHwKuBOaN9Z6SJElT2qatqybzNZgLgFO7X58K/J+RDboFyq8Cn62qL424Nrf7a+gUPK/e2g0nbM5qdxXYy+isHBvt+uJNcyfuuOOOieqWJEmStt2fAC9McgPwwu4xSRYm+VS3zSuB5wCvH2WLqs8luQq4CpgD/PHWbjghDwVIsgNwHvCXVXXjaG2qagmduQsce+yTayL6JUmSpG1XVXcCzx/l/FLgjd2v/xb42828/1f6vWdfg9UkpwO/3T18aXc7gm2xBLihqj7az/0kSZK2LwU16QusppS+pgFU1dk9m7hudqCa5JtJDux+/cfAHsDbB+uqJEmSpptBtq7aH1gKPArYmOTtwALgPuAwYHWSeXS2OfgxcGVnLi0fq6pPjZ4qSZK0nds48FOkppUxD1ar6jZGWdmf5EjgK1X1ALACyNi7J0mSpOms+QKrqroaOLN1riRJkqafCdkNQJIkSXT2WXWBVV8mbJ9VSZIkqV9WViVJkiZMtXiK1LQylIPVGaxnl/yiaeZN19/XNG+TQ4+7v3nmnTs/pXkmwJxd2/9277Tq280zAZavb/8ZPJB9mmfOmvvi5pkAO117dvPMjfP/S/NMgHXrNjbPPPCBUR90N7ANOx3VPHPdzvObZwIcfeu7mmd+7ea/aZ75soNf2zwT4N/WfLF55g43tP/+AVbtd0rzzLp7s08yH8h+K89pHzp79+aR63c/rHmmpi6nAUiSJGloDWVlVZIkaXtVLrDqi5VVSZIkDS0Hq5IkSRpaTgOQJEmaKFU+brVPVlYlSZI0tKysSpIkTRj3We3XmCurSY5IcmmStUl+bwvtzkpyS5Lx2ehUkiRJ261BpgGsBn4X+LOttPsacPwA95EkSdI0NeZpAFW1CliV5Fe30u4ygCRjvZUkSdL2w31W++ICK0mSJA2toVlglWQxsBjgoPkHTnJvJEmSxkG5wKpffVVWk5yeZFn3dUDLjlTVkqpaWFUL5+yzd8toSZIkTVF9DVar6uz/1979B9tV1Xcff38DBKSCgYQfIQmCBEtTfoRyq3Zs4UFRQfsItujgQ5FOg2kdHOtQH8HB+thHW6l2ascBf2Swj7RV0QqOVIuIkVZbfmgQhAhKKD8khh8KUX6FQHK/zx9nxx7DvTl7n7POvfvc+37NnLlnn7P2Z697V7LPuuuuvXZmLq8eGyYrFxGrI8LhUUmSJA2k72kAEbE/sAbYExiPiLcDy4DHgaV0VgsgIj4I/C9g94hYD1ycme8dsN6SJEkjKCG9g1UTg6wG8ACwePvXI+Jw4LLM3FSVeyfwzr5rKEmSpFmr+AVWmbkWOKd0riRJ0ozgBVaNuHSVJEmSWsvOqiRJklqrNeusSpIkzXzpHawacmRVkiRJrdXKkdWnng5uv2enoplHjRW9h8Ev3PHwU8UzD9vla8UzATbOeVXxzB8+clTxTIAj511XPvS+B4tHPrr4tOKZAE8e/MfFM/d+8jvFMwF23jWKZ962+beLZwIs++7/KZ45/uLymQD/+sSfFc/c9ZYHimd++8nPF88EeNHubyie+ZX1/1Q8E+CArWuLZ377ruF8Zu11TPlzy9wb/6p45oyWeIFVQ46sSpIkqbXsrEqSJKm1WjkNQJIkaWZKGPcOVk04sipJkqTWcmRVkiRpyrh0VVOOrEqSJKm17KxKkiSptfrurEbEYRFxXURsjoh37KDcMRFxa0TcGREfiYjyizJKkiSNgm3rrE7nY8QMMrL6CPA24G96lPsYsBI4tHqcOMAxJUmSNIv03VnNzIcy8zvAM5OViYiFwJ6ZeV1mJvAPwCn9HlOSJEmzy7BXA1gErO/aXl+99iwRsZLOCCwLF05YRJIkacS5GkBTw77AaqL5qTlRwcxclZljmTm2117zh1wtSZIkjYJGI6sRcTbw5mrz1Zm5occu64HFXduLgV77SJIkzVBJegerRhqNrGbmRZm5vHpM2umMiNURsSgz7wcei4iXVKsAvAn40oB1liRJ0izR95zViNgfWAPsCYxHxNuBZcDjwFI6qwUAvAX4FPAc4MrqIUmSJPXUd2c1Mx/gl//ED0BEHA5clpmbqnJrgMP7rqEkSdIMkQk5gmudTqfiqwFk5lrgnNK5kiRJmn2GvXSVJEmSfiEdWW1o2EtXSZIkSX2zsypJkqTWauU0gN3mjvNrBz5VNHPnuy8tmrfNE3ufWTzzB5t+p3gmwGHrPlI88/7nriieCcDGO4tHPv3CPyyeuf7eJ4pnAvxkw6PFM497zjeKZwIsWFi+rmvWHVo8E+CA4/6meOa8jf9ZPBPg0KN+o3jmIfuUb6ud1/1j8UyAr6z/p+KZr1n8B8UzAb58X/mfwdObHy6eCXDb7RuLZx75wlOLZ85o6TSAphxZlSRJUmu1cmRVkiRppsqtjqw24ciqJEmSWsvOqiRJklrLaQCSJElTxQusGnNkVZIkSa3lyKokSdIUSe9g1VijkdWIOCwirouIzRHxjpr7zI+IayLi8Yi4sL9qSpIkaTZqOrL6CPA24JQG+zwF/DlwePWQJEmSamnUWc3Mh4CHIuI1DfZ5AviPiFjatHKSJEkzzbjTABrxAitJkiS1VmsusIqIlcBKgAMXHzDNtZEkSRoCl65qrOfIakScHRE3V4+evciIeF1X+bG6FcnMVZk5lpljCxbsXXc3SZIkzWA9R1Yz8yLgol7lImI18KbM/CLwxQJ1kyRJ0izXaBpAROwPrAH2BMYj4u3AMuBxYCmd1QIm2u+eap+5EXEK8MrMvG2AekuSJI0kpwE003Q1gAeAxdu/HhGHA5dl5qZJ9juor9pJkiRpVitygVVmrgXOKZElSZIkbdOa1QAkSZJmuswkx7dMdzVGiuusSpIkqbUcWZUkSZoyrrPalCOrkiRJqiUi9o6IqyNiXfV1r0nKbe1ad/+KrtcPjogbqv0/FxFzex2znSOruZU5z/y8aORDi1YUzdtmzubyvx3dcv09xTMBDvsfRxTP3PdXnlM8E4CH9igeefVX1xfPfMWJz1oco4h999mvfOijR5fPBO7d6djimb96VBbPBHjggSeLZz6x54uLZwL86iMfKZ65dcFpxTMf2u+NxTMBDti6tnjml+/7x+KZAL+75IzimZ/41oXFMwGOOPCx4pn3PrakeObBm75ZPFPFnAeszswLIuK8avvcCcptyszlE7z+18CHM/PSiPg4sAL42I4O6MiqJEnSVMnOOqvT+RjQycAl1fNLgFPq7hgRAbwM+EKT/e2sSlwu0NkAABhjSURBVJIkqa79MvN+gOrrvpOU2y0i1kTE9dUNoQDmAz/LzG3LIawHFvU6YDunAUiSJM1ISW6d9gusFkTEmq7tVZm5attGRHwd2H+C/c5vcIwDM3NDRLwA+EZE3Ao8OkG5nvO+7KxKkiTNLj/NzLHJ3szMEyZ7LyIejIiFmXl/RCwEHpokY0P19a6I+DfgaOAyYF5E7FyNri4GNvSqrNMAJEmSVNcVwJnV8zOBL21fICL2iohdq+cLgJcCt2VmAtcAp+5o/+05sipJkjRFOnewmvZpAIO4APh8RKwAfgS8HiAixoA/ycyzgF8DPhER43QGRi/IzNuq/c8FLo2I9wM3AZ/sdUA7q5IkSaolMx8GXj7B62uAs6rn1wITrpeZmXcBL2pyzL47qxFxOv+9rtbjwFsy83sTlPsqsLA61reAszNzpH+lkCRJ6tf4aI+sTrlB5qzeDRyXmUcC7wNWTVLuDZl5FHA4sA/VcLEkSZLUS98jq9UQ7zbX07mia6Jy25Yp2BmYS40lCiRJkiQoN2d1BXDlZG9GxFV05idcyX/ftWD7MiuBlQAHLl5YqFqSJEltMvIXWE25gZeuiojj6XRWJ7ovLACZ+So681Z3pXObrYnKrMrMscwcWzB/3qDVkiRJ0gzQaGQ1Is4G3lxtvhpYAFwMnFRdHTapzHwqIq6gc0/Zq/uoqyRJ0mhLHFltqNHIamZelJnLM3M5nY7u5cAZmXlHd7mIWB0RiyLiudXdDYiInel0cH9QqO6SJEma4QaZs/oeYD7w0YgA2JKZYxExB1gKPALsCVxR3cVgJ+AbwMcHq7IkSZJmi0FWAziLavHX7SwDLsvMTcAm4Df7PYYkSdJMkl5g1VjxO1hl5lrgnNK5kiRJmn0GXg1AkiRJGpbiI6uSJEmaRDoNoClHViVJktRarRxZjWeeYM4DN5QN3e8FZfMqj/58c/HMV//+rxfPBHh0fFnxzL22rC+eCcAQfus84TVLi2dufnq8eCbAvfc8UjxzwcavFs8EeP4etxTP/DZ/VDwT4Jilj/Yu1NDOm4bzf2DjYeWn/u++W/nxifz5M8UzAb591wHFM5/evMPlwPv2iW9dWDzzj3/nrcUzAW588APFMw+es6Z45uZ9jy2e2SY5vmW6qzBSHFmVJElSa9lZlSRJUmu1chqAJEnSjJRJbvUCqyYcWZUkSVJrObIqSZI0RZJk3KWrGnFkVZIkSa1lZ1WSJEmt1fc0gIg4HTi32nwceEtmfm+7MrsD/wwcAmwF/iUzz+v3mJIkSaPOO1g1M8jI6t3AcZl5JPA+YNUk5f4mMw8DjgZeGhEnDXBMSZIkzSJ9j6xm5rVdm9cDiyco8yRwTfX86Yj47kTlJEmSZoVMR1YbKjVndQVw5Y4KRMQ84H8Cqyd5f2VErImINT/Z+HihakmSJGmUDdxZjYjj6XRWz91BmZ2BzwIfycy7JiqTmasycywzx/bZ67mDVkuSJEkzQKNpABFxNvDmavPVwALgYuCkzHx4B7uuAtZl5t/1VUtJkqQZwmkAzTQaWc3MizJzeWYup9PRvRw4IzPv6C4XEasjYlH1/P3A84C3F6qzJEmSZolB7mD1HmA+8NGIANiSmWMRMQdYCjwSEYuB84EfAN+tyl2YmRcPVm1JkqTRk5mMj+d0V2OkDLIawFnAWRO8tQy4LDM3AeuB6PcYkiRJmt0GGVmdUGauBc4pnStJkqTZp3hnVZIkSZNzFkAzpdZZlSRJkopzZFWSJGkKpUOrjTiyKkmSpNZq58jqTrvCvEOKRu772NVF87aZv/9hxTPX/Xg4v0MsWVL+zmA/2bSweCbA3ktPL54555mfFc/c9NSvFM8E+MmGR8uHLjy4fCbw9OFnF8+c+8Ph3HL53scPKJ75gg2fK54JcOfOv148c9kRC4pn7rfh74tnAux1zB8Xz7zt9o3FMwGOOPCx4pk3PviB4pkAx+z3ruKZn7zhY8Uzj5q3o/sMabZpZ2dVkiRpBsqE8fHprsVocRqAJEmSWsvOqiRJklrLaQCSJElTaDxdDaAJR1YlSZLUWo6sSpIkTSEvsGrGkVVJkiS1Vt+d1Yg4PSJuqR7XRsRRk5T7y4i4LyKGs3CiJEmSZqxBpgHcDRyXmRsj4iRgFfDiCcr9C3AhsG6AY0mSJI28TG+32lTfndXMvLZr83pg8STlrgeIiH4PJUmSpFmq1AVWK4ArBwmIiJXASoADF+1bok6SJEmt48BqMwNfYBURx9PprJ47SE5mrsrMscwc22fv5w1aLUmSJM0AjTqrEXF2RNxcPQ6IiCOBi4GTM/Ph4VRRkiRJs1WjaQCZeRFwEUBEHAhcDpyRmXd0l4uI1cCbMvPHpSoqSZI0E4w7D6CRQaYBvAeYD3y0GmldAxARc4ClwCPV9gcjYj2we0Ssj4j3DlhnSZIkzRKDrAZwFnDWBG8tAy7LzE1VuXcC7+z3OJIkSTNF4gVWTRW/3WpmrgXOKZ0rSZKk2cfbrUqSJKm1io+sSpIkaRKZ3sGqIUdWJUmS1FrtHFl9+nFY/62ikfn8E4rm/cKcXYpH7jJ3OL9DbN68tXjmgt0eKZ4JsNPG9cUzH4gjimcufPCS4pkASw55ffHMLQe+qXgmwC4bby6eueywY4pnAvz7N+4rnvn83xrOz/WgrfOKZ+76xH8Vz2TuHuUzgbk3/lXxzCNfeGrxTIB7H1tSPPPgOWuKZwJ88oaPFc9c8eK3FM+8acO7i2e2iUtXNePIqiRJklrLzqokSZJaq53TACRJkmYg11ltzpFVSZIktZYjq5IkSVMlvcCqKUdWJUmS1Fp2ViVJktRafXdWI+L0iLilelwbEUdNUu6YiLg1Iu6MiI9ERPRfXUmSpNGW49P7GDWDjKzeDRyXmUcC7wNWTVLuY8BK4NDqceIAx5QkSdIs0ndnNTOvzcyN1eb1wOLty0TEQmDPzLwuMxP4B+CUfo8pSZKk2aXUagArgCsneH0R0H3fzPXVa88SESvpjMBy4H7PK1QtSZKk9uiss+pqAE0MfIFVRBxPp7N67kRvT/DahC2Umasycywzx/aZt/ug1ZIkSdIM0GhkNSLOBt5cbb4aWABcDJyUmQ9PsMt6fnl6wGJgQx/1lCRJGn0J4yN4kdN0ajSympkXZebyzFxOp6N7OXBGZt7RXS4iVkfEosy8H3gsIl5SrQLwJuBLpSovSZKkmW2QOavvAeYDH61Wo9qSmWMRMQdYCjxSlXsL8CngOXTmtU40t1WSJEl6lr47q5l5FnDWBG8tAy7LzE1VuTXA4f0eR5IkaaZIvN1qU8XvYJWZazPznNK5kiRJml4RsXdEXB0R66qve01Q5viIuLnr8VREnFK996mIuLvrveW9juntViVJkqbQeE7vY0DnAasz81BgdbX9SzLzmq5rnF4GPAl8ravI/972fmbe3OuAdlYlSZJU18nAJdXzS+h9s6dTgSsz88l+D2hnVZIkaXZZEBFruh4rG+y7X7XaE9XXfXuUPw347Hav/WVE3BIRH46IXXsdsNQdrCRJktRLJjn9F1j9NDPHJnszIr4O7D/BW+c3OUhELASOAK7qevldwAPAXGAVnZtK/d8d5bSyszq++/5sPvpZUyBmjcUjdAOvrRP+Wy6Qu1v53L2LJ8LmvSZaEGNwS4eQuXUImQBb588fUnJ5x778oOKZW4onduwxhMwtQ/iXtWWPYfxrHS0Hb/pm8czN+x5bPBPgqHkT3b9nMDdteHfxzKMPeH/xTNWXmSdM9l5EPBgRCzPz/qoz+tAOot4AfDEzn+nKvr96ujki/h/wjl71cRqAJEnSFNm2dNV0PgZ0BXBm9fxMdnyzpzey3RSAqoNLdbOoU4C1vQ5oZ1WSJEl1XQC8IiLWAa+otomIsYi4eFuhiDgIWAL8+3b7fzoibgVuBRYAPYfRWzkNQJIkSe2TmQ8DL5/g9TV03SwqM+8BFk1Q7mVNj2lnVZIkaQpN//VVo8VpAJIkSWotR1YlSZKmSlLiIqdZpdHIakScXi3iektEXBsRR9XYZ35EXBMRj0fEhf1XVZIkSbNN05HVu4HjMnNjRJxEZzHXF/fY5yngz4HDq4ckSZJUS6POamZe27V5PbC4xj5PAP8REa4cLUmSZrUEcny6azFaBrnAagVwZamKSJIkSdvr6wKriDieTmf1t0tVJCJWAisBlixZUipWkiSpVcbTC6ya6DmyGhFnR8TN1eOAiDgSuBg4uVoYdvvyr+sqP1a3Ipm5KjPHMnNswYIFzb4LSZIkzUg9O6uZeVFmLs/M5XRGYi8HzsjMO7rLRcTqiFiUmV/cVr66m4EkSZLUl6bTAN4DzAc+GhEAWzJzLCLmAEuBRybaKSLuAfYE5kbEKcArM/O2vmstSZI0ihLGvcCqkaarAZxF131fuywDLsvMTZPsd1DzqkmSJGm2K3IHq8xcC5xTIkuSJEnaxtutSpIkTZEEvNtqM4OssypJkiQNlSOrkiRJU8iR1WYcWZUkSVJrtXJk9aabbvrp7rvvfm/N4guAnxauwjAyh5VrXUenrrP9+x9W7myv62z//oeVO9vrOlO//+cXPramQCs7q5m5T92yEbEmM2vfKWu6MoeVa11Hp66z/fsfVu5sr+ts//6HlTvb6zrbv/9h8gKr5pwGIEmSpNZq5ciqJEnSjJSOrDY1E0ZWV41I5rByrevo1HW2f//Dyp3tdZ3t3/+wcmd7XWf7968WiUy795IkSVPhkL3m5AUv32Va6/CGy56+cZTm+ToNQJIkaQo5DaCZmTANQJIkSTPUyHRWI+LEiPhhRNwZEedN8P7eEXF1RKyrvu5VIPP1EfH9iBiPiNrD5TVyPxQRP4iIWyLiixExr0Dm+6q8myPiaxFxQIm6dpV7R0RkRCzokff3EfFQRKyd5P1dI+Jz1fFuiIiDatazV+6xEfHdiNgSEacWyjwnIm6rfq6rI6LW+nw1cv8kIm6t2uo/ImLZoJld5U6t2qnnv9eIWBIR10TE7dW/8z+doEzj9qqZ26i9amY2aq+amf20Vc/crrK12isidouIb0fE96rMv5igTD/nwDq5jc6DNTP7OQfWyW10HqyT2VW21jmwq/xOEXFTRHx5gvcat1XN3H4/s3aU2bitamT29XnVK7erTKO2mmqdpatiWh+jZiQ6qxGxE3ARcBKwDHjjBB8a5wGrM/NQYHW1PWjmWuD3gG8WruvVwOGZeSRwB/CuApkfyswjM3M58GXgPYXqSkQsAV4B/KhXJvAp4MQdvL8C2JiZS4EPA39dI7NO7o+APwQ+UzOvTuZNwFjVTl8APlgo9zOZeUTVVh8E/rZAJhGxB/A24IZ61WQL8GeZ+WvAS4CzJ2j/ftqrTm7T9qqT2bS96mT201Z1cpu212bgZZl5FLAcODEiXrJdmUbnwAa5Tc+DdTIbnQMb5DY9D9bJbHoO3OZPgdsnea+ftqqT2/gzq0ZmP23VK7Px51XN3H7bSi03Ep1V4EXAnZl5V2Y+DVwKnLxdmZOBS6rnlwCnDJqZmbdn5g9L1zUzv5aZW6rN64HFBTIf7dr8FTq/vA1c18qHgXfWyczMbwKP7KBIdzt9AXh5RPT8Na9Xbmbek5m3AOO9shpkXpOZT1abddqpbm7jtqrxcwV4H50O1VM1qklm3p+Z362eP0bnA2DRdsUat1ed3KbtVTOzUXvVzOynrer8XKFBe2XH49XmLtVj+7o0PQfWym16HqyZ2fQcWDe3UXvV/LlCg3MgQEQsBl4DXDxJkcZtVSe3n8+sGpmN26pGZj+fV3V+rtCwrTQaRqWzugi4r2t7Pc8++e+XmfdD58MC2LdAZj+a5v4RcGWJzIj4y4i4Dzider+p9syNiNcCP87M79XIq+MXx6xOgD8H5hfKHqYV9G6n2iLi7Ij4LzqdlbcVyDsaWJKZk/5prMf+BwFH8+xRvoHaawe5fauZ2ai9dpQ5SFtNlttPe1V//rwZeAi4OjO3r2vTc2Dd3MYaZtY5B9bObXoe7JXZ5znw7+h0mCb7ZayvtqqR248mmXXbqmdmH59XPXOH8Hk1HAnj49P7GDWj0lmdaCRn0N+ahpHZKDcizqfzJ8NPl8jMzPMzc0mV99YemT1zI2J34Hya/YlmoGO2UUT8ATAGfKhUZmZelJmHAOcC7x4kKyLm0BlN+LM+938ucBnw9u1GPGCA9uqR25c6mU3bq1dmv201WW6/7ZWZW6s/my4GXhQRhzfZfypz62Y2OAfWzm16HtxRZj/nwIj4XeChzLyx7j7Tldsks25b1c1s2k69cof0eaWWGJXO6npgSdf2YmDDdmUejIiFANXXhwpk9qNWbkScCfwucHpmz8Vum9b1M8DvF6jrIcDBwPci4p7q/e9GxP41snseMyJ2Bp5H7z9vT5uIOIHOCfC1mbl5CIe4lJp/AtyBPYDDgX+r2uklwBV1LrCIiF3odKg+nZmXT1Ckr/aqkdtYncym7dWwnrXbqkdu3+0FkJk/A/6NZ89hbnoOrJvbtx1lNjwH1s7tUvc8uKPMfs6BLwVeW5W/FHhZRPzTdmX6aas6uU3VymzYVk3rWbedeuUO4/NqKDoXWE3vY9SMSmf1O8ChEXFwRMwFTqNzcv9ARLyuKnMFcGb1/EzgSwUyh1LXiDiRzijNa7vm2A2aeWhX+dcCPxg0NzNvzcx9M/OgzDyITsflNzLzgTo/iG0i4q0Rse035+52OhX4RtMPqklyi+jOrP5U+wk67dTog79HbndbvQZYN0hmZv48Mxd0tdP1VZ3X9Ng/gE8Ct2fm33a9PlB71cxtpE5m0/aqmdm4rXrl9tNeEbFPVFdhR8RzgBOAHwx4Dqyb20idzD7OgXVzG50He2X2cw7MzHdl5uKq/Gl0/s/8waBtVTO3kTqZTduqZmbjz6teuaU+r9ROI9FZrebKvRW4is7FCp/PzO8DRwDb/iFeALwiItbRuRLwgkEzI+J1EbEe+C3gKxFxVaG6XkhndOXq6Czd8fES339ErI2IW4BX0rliskRdG4mIzwLXAb8aEesjYgVwGPBwVeSTwPyIuBM4h5pXwfbKjYjfrNrq9cAnIuL7Ber6IeC5wD9X7XRFiboCb43O8jI3Vz+DMyeJapLZj5cCZ9AZobi5eryawdurZ24f7VWnrk3bq05m47aqmdvUQuCa6v/3d+jMrfwyA5wD6+b2cR6sU9dG58AGuU3Pg3UySxm0rXrm9vOZVaOu/bRVr8zGn1c1czVDjfTtViPiqsx8Vdszh5U7KnWNznp4v5edFQeKGUaudR2t3FHJHFbuqJwDhpU5rNxRyRxW7qhkDjN3mA6eNyffe+zcaa3DH/7L5pG63epId1YlSZJGiZ3V5kZiGoAkSZJmp52nuwKSJEmzxbbVAFSfI6uSJElqLUdWJUmSppAjq804sipJkqTWsrMqSZKk1nIagCRJ0lQZ0VueTidHViVJktRajqxKkiRNkQS2OrLaiCOrkiRJai07q5IkSWotpwFIkiRNIS+wasaRVUmSJLWWI6uSJElTJIHx8emuxWhxZFWSJEmtZWdVkiRJreU0AEmSpKniHawac2RVkiRJreXIqiRJ0hRJHFltypFVSZIktZadVUmSJLWW0wAkSZKmkNMAmnFkVZIkSa3lyKokSdIUGs+Y7iqMFEdWJUmS1Fp2ViVJktRaTgOQJEmaIq6z2pwjq5IkSWotO6uSJElqLacBSJIkTZV0GkBTjqxKkiSptRxZlSRJmiIJjI9Pdy1GiyOrkiRJai07q5IkSWotpwFIkiRNIS+wasaRVUmSJLWWI6uSJElTyJHVZhxZlSRJUmvZWZUkSVItEfH6iPh+RIxHxNgOyp0YET+MiDsj4ryu1w+OiBsiYl1EfC4i5vY6pp1VSZKkKZLVHaym8zGgtcDvAd+crEBE7ARcBJwELAPeGBHLqrf/GvhwZh4KbARW9DqgnVVJkiTVkpm3Z+YPexR7EXBnZt6VmU8DlwInR0QALwO+UJW7BDil1zG9wEqSJGmK3P9EXvUX//nUgmmuxm4RsaZre1VmriqYvwi4r2t7PfBiYD7ws8zc0vX6ol5hdlYlSZKmSGaeON116CUivg7sP8Fb52fml+pETPBa7uD1HbKzKkmSpF/IzBMGjFgPLOnaXgxsAH4KzIuInavR1W2v75BzViVJklTSd4BDqyv/5wKnAVdkZgLXAKdW5c4Eeo7U2lmVJElSLRHxuohYD/wW8JWIuKp6/YCI+FeAatT0rcBVwO3A5zPz+1XEucA5EXEnnTmsn+x5zE4nV5IkSWofR1YlSZLUWnZWJUmS1Fp2ViVJktRadlYlSZLUWnZWJUmS1Fp2ViVJktRadlYlSZLUWv8fDv6KZR8QFo4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x864 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['figure.figsize'] = 12,12\n",
    "plt.imshow(truCov, cmap=cmap, interpolation='Nearest', vmin=-1, vmax=1)\n",
    "plt.xticks(range(nbPix), pixCoords)\n",
    "plt.yticks(range(nbPix), ylabs)\n",
    "plt.title('Cov entry to shift, pixel {}'.format(pixIdx))\n",
    "plt.colorbar()\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 8,8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And same for those we get from our estimators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAKkCAYAAACplz9YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8ZXdZH/7PMzEBLIRAhktIwkWJAgoGTMFWfl64NehPgnfQSlRotIVa29IK0h8CXgrWX2l9NVangtCKQbygQcMvIhethWAChkCIQEAwQwKYhHAnkMzz+2Ovkc3hzJxZa+asPfvM+/167dfsdX/2Ojsnz3me9V2rujsAALDddq06AAAAjg0STwAAZiHxBABgFhJPAABmIfEEAGAWEk8AAGYh8QQA2KGq6sVV9ZGqescBlldV/XJVXV1VV1TVQ5aWnVtV7xle5x6JeCSeAAA710uSnH2Q5Y9NcsbwOi/Jf0+Sqrpzkp9J8rAkD03yM1V1p8MNRuIJALBDdfefJ7nxIKuck+R/9sIlSU6qqlOS/JMkr+nuG7v7o0lek4MnsIfkyw53BwAAx6I71F371nxupTF8Jh+7Mslnl2bt6e49I3ZxapJrlqb3DvMONP+wSDwBACa4NZ/LffPwlcbw9vzxZ7v7rMPYRW0yrw8y/7BotQMAHLv2Jjl9afq0JNceZP5hkXgCABy7LkzypGF0+zck+Vh3X5fk4iSPqao7DYOKHjPMOyxa7QAAE+3ateIa3r6DL66qC5J8S5LdVbU3i5HqxydJd/9qkouSfFuSq5N8OsmPDMturKqfTXLpsKvndffBBikdkuo+7HY9AMAx58vrpP6qXd+00hjetu9VbznMazxnpeIJADBJpeq4VQexVlzjCQDALCSeAADMQqsdAGCCqqNgcNGtqz38WCqeAADMQsUTAGASg4vGUvEEAGAWEk8AAGah1Q4AMNHKBxetGWcLAIBZqHgCAExSqVLDG8PZAgBgFhJPAABmodUOADDB4slF7uM5hoonAACzUPEEAJjE4KKxnC0AAGYh8QQAYBZa7QAAExlcNI6KJwAAs1DxBACYoMrgorGcLQAAZiHxBABgFlrtAAAT7dqlhjeGswUAwCwkngAAzEKrHQBgkkqV+3iOoeIJAMAsVDwBACaoMrhoLGcLAIBZSDwBAJiFVjsAwEQemTmOswUAwCxUPAEAJqns2uV2SmOoeAIAMAuJJwAAs9BqBwCYyOCicZwtAABmIfHkmFZVP11Vv36I676kqn5uu2M6HFX1yar6ilXHMaequrKqvmXVcWyXMd/RLfZz0O9vVd2lqt5VVbc93GMdQix/WVVfs93Hge1WtRhctMrXupF4clSrqvdX1WeGhOrDVfUbVXX7ifv6lqrauzyvu3+hu59yZKJdve6+fXe/b9VxbJfNkqfu/prufsM2HOtLvi+rMON39BlJfqO7P5skVfVLVfWeqvpEVf11VT3pUHdUVU8YktiPVdVHquqlVXXi0iq/lOR5Rzh+YA1IPFkH39Hdt0/ykCT/MMl/GLuDqtrR1zPv9M/H9qqq2yQ5N8lvLs3+VJLvSHLHYdl/rap/fIi7/D9JvrG775jkK7IYT7D8B8OFSb61qk453NiB9SLxZG109weTvDrJ1yZJVf1IVV01VGTeV1U/tn/d/dWqqvqpqvpQkguGbe8xVE8/WVX3qKrnVNVvLm33O1X1oaFS8+dj2oFV9c+W4nlnVT1kmH//qnpDVd00tIUfN8z/huFYxy3t4zur6orh/UOr6k3DdtdV1X+rqhOW1u2qempVvSfJe5bm3Xd4/+1V9VdV9fGquqaqnrO07b2Hdc+tqr+tquur6llLy48bWrzvHT7PW6rq9GHZ/arqNVV141DV+r6l7b5t+OyfqKoPVtXTD3K+fnQ4Xx+tqour6l7D/KqqFw6Vso9V1RVV9bVVdV6SH0zy74ef36uG9d9fVY8a3j9n+Bn+5hDD26vqq6rqmcP+rqmqxyzFsOl3qKr+wQG+L7uq6hnDebmhql5RVXcetrntcNwbhp/ZpVV1t0P43uz/WZxXVdcOP+t/u7T877+jVfX9Q5wnDtOPHb5Dd9nqZ7OFhyW5qbv/vsLb3T/T3X/d3fu6+81J/neSf3QoO+vua7r7+qVZtya579LyzyZ5S5LHbNwW1s2uXbtW+lo36xcxx6wh8fm2JH81zPpIkv87yYlJfiTJC/cne4O7J7lzknsleVKSxya5dmhH3767r93kMK9OckaSuyZ5a5KXHWJs35vkOcNxTkzyuCQ3VNXxSV6V5E+Gff7LJC+rqq/u7kuyqCo9YmlXP5Dkt4b3tyb510l2Z/E//Ecm+RcbDv34LJKGB2wS1qeGeE5K8u1J/nlVPX7DOg9P8tXDvp9dVfcf5v+bJE/M4nyfmORHk3x6SMheM8R412GdX6kvJOgvSvJj3X2HLP5AeN0Bztfjk/x0ku9KcpcskpoLhsWPSfJNSb5qiP37k9zQ3Xuy+Hn84vDz+47N9p1Fle5/JblTFt+Vi7P4XXdqFu3dX1tad9PvUHd/Kpt/X34ii3P+zUnukeSjSc4f9nVuFtXB05OcnOTHk3zmADFu5luz+O49Jskz9ifTy7r7t5O8KckvV9XJWZzvp3T33x3Cz+ZgHpjkXQdaWFW3y6LbcOWhfpiqenhVfSzJJ5J8d5L/smGVq5J83aHuD9gZJJ6sgz+oqpuS/EWSP0vyC0nS3X/c3e/thT/LIrn7v5a225fkZ7r75u4+pASgu1/c3Z/o7puzSCS/rqrueAibPiWLhOjSIZ6ru/sDSb4hye2TPL+7P9fdr0vyR1kkBcki2XpiklTVHbJI9C4YYnlLd1/S3bd09/uzSJi+ecNx/2N337jZ5+vuN3T324eK1RXDfjdu/9zu/kx3vy3J2/KFROApSf5Dd79r+Dxv6+4bskjS3t/dvzHE9dYkv5fke4btPp/kAVV1Ynd/dFi+mR8bYr+qu2/J4md65lD1/HySOyS5X5Ia1rnuAPvZzP/u7ouH/f5OFont87v780lenuTeVXXScI62+g5tFvezunvv0nfke2pxqcPns0g479vdtw4/v4+PiPu53f2p7n57kt/IF74jGz01iz9W3pDkVd39R8P8rX42B3NSFgnigfxqFt+Piw9hX0mS7v6LodV+WpL/lOT9G1b5xHBcWGOVql0rfa2b9YuYY9Hju/uk7r5Xd/+L/UnW0Ga8ZGgr3pRF0rZ7abu/2z9Q4lAM7eXnD23Uj+cL/6PcfZDN9js9yXs3mX+PJNd0976leR/IovqWLKpT31WLa+y+K8lbh4Q1Q4v4j4ZW6sezSM42xnLNQT7Pw6rq9VX1d0Pl6cc32f5DS+8/nUWSfLDPc68kDxtayTcN5/0Hs6guJ4vK1rcl+UBV/VlVHag1e68srhncv48bk1SSU4fk/L9lUUn8cFXtqS8emLKVDy+9/0yS67v71qXp7P+ch/Ad2izuVy7FfVUWlem7ZVFlvTjJy4eW+S8OFe9Dtfyz/EAW350v0d03ZZFQf22S/3dDbAf72RzMR7NI9r9EVf2n4Vjf1919CPvaGO8Hk/x/WST9y+6Q5Kax+wPWm8STtTQkar+XxejYu3X3SUkuyiJ52W/j/yS3+p/mDyQ5J8mjsmiZ3nv/4Q4hpGuSfOUm869Ncnp98Z+l90zywSTp7ndmkWQ8Nl/cZk+S/57kr5Oc0d0nZtGa3hjLwT7Tb2UxiOP0ofL0q4f4WQ72ea5J8mfDHwL7X7fv7n8+fJ5Lu/ucLFq9f5DkFQfZ/49t2M/tuvuNw35+ubu/PsnXZNFy/3eH8HlHOYTv0GbHuibJYzfEfdvu/mB3f767n9vdD0jyj7OoQB7ySPAskv397pnFd2ezuM/M4tKHC5L88obYDviz2cIVWZznjcd6bhbfzceMrN5u9GX50u/T/bOoogLHEIkn6+qEJLdJ8ndJbqmqx2brgQofTnLyQVrnd0hyc5Ibknx5hpb+Ifr1JE+vqq+vhfsObeM3Z3Gt5b+vquNrcb/J78gXV39+K4trB78pi0rWcjwfT/LJqrpfkkNJIDZ+nhu7+7NV9dAsEtsxn+dnq+qM4fM8aLim8I+SfFVV/dDweY6vqn9YiwFUJ1TVD1bVHYe29sezqAZu5leTPHP/9YdVdcfhOtkM+3vYUC38VJLPLu3nw1mMkj4StvoObfZ9+dUkP19fGAh1l6o6Z3j/rVX1wFoMFvt4Fq33W4dlz6mqN2wRz/9TVV8+nJMfSfLbG1eoxT02fzOLP0J+JMmpVbX/ut8D/mwO4Vz8ZZKTqmp/JT5V9cwsvjOPHi6z2BjL+6vqhzfb2fA9uOfw3blXkp9P8tql5bdJ8vVZXJMKa6sq7uM5ksSTtdTdn8giWXtFFm3CH8iiunewbf46iyrR+4ZW5MZW5v/Movr4wSTvTHLJiHh+J4v/uf5WFteu/UGSO3f357IYaPTYJNcn+ZUkTxpi2e+CJN+S5HX9xSOBnz58rk8k+R/ZJBHZwr9I8ryq+kSSZ+fA1cfN/Odh/T/JIol6UZLbDef9MUmekEVF7kNJXpBFApckP5Tk/cOlAT+e5J9utvPufuWw3cuHdd+RxTlKFgN9/kcWP9cPZPGHwC8Ny16UxTWkN1XVH4z4PJvFcNDv0AG+L/91WOdPhvN6SRaDu5JFS/t3szhfV2VxPfL+OyacnsUthg7mz5JcnUWC9kvd/SebrPMfk+zt7v8+XGP6T5P8XFWdcQg/m4Odi88leUm++Of1C1lUXt9TXxjZ/9NJUou7K5ycA/838oAkb0zyySw+97uS/LOl5Y9L8obefIAfsIPVhEt2ABihqi5P8sgDVA7vneRvkhw/DIhaiVrckul/J3nwZoPVNqz78CRP7e4DDYDa6lhvTvLk7n7HlO3haHHSbe7R33SP1T6D5FXv/9m3dPdZKw1iBDedBthm3X3mqmPYSnf/XRZ3EjiUdf8ii7tMTD3Ww7ZeC9iJtNoBAJiFiifACg33aD3Uuw0AR5l1HOCzSiqeAADMYtsrnrtPvnPf656nbfdhdozP3eJvgbE+/7kD3bGHA/nyWw5433k2c8cjdQenY0etbpzUWvK7f7x3vONt13f3XVYdB+Nse+J5r3uelje94aLtPsyO8bc33G7VIaydD13zsVWHsHYe8pF/veoQ1sq+b7tg65X4IifccuOqQ1grfvePd8ZX7v7AqmOoqrV8bOUqOVsAAMzC4CIAgIl27VLDG8PZAgBgFhJPAABmodUOADBJpcp9PMdQ8QQAYBYqngAAE1QZXDSWswUAwCwkngAAzEKrHQBgol3H1apDWCsqngAAzELFEwBgiqrUcWp4YzhbAADMQuIJAMAstNoBACaoGFw0loonAACzUPEEAJiiVDzHUvEEAGAWEk8AgB2sqs6uqndV1dVV9YxNlr+wqi4fXu+uqpuWlt26tOzCw41Fqx0AYIo1uI9nVR2X5Pwkj06yN8mlVXVhd79z/zrd/a+X1v+XSR68tIvPdPeZRyqeo/tsAQBwOB6a5Orufl93fy7Jy5Occ5D1n5jkgu0KRuIJALC+dlfVZUuv8zYsPzXJNUvTe4d5X6Kq7pXkPkletzT7tsN+L6mqxx9usFrtAAATHCX38by+u886yPLNAuwDrPuEJL/b3bcuzbtnd19bVV+R5HVV9fbufu/UYFU8AQB2rr1JTl+aPi3JtQdY9wnZ0Gbv7muHf9+X5A354us/R1PxBACYopJdR/ngoiSXJjmjqu6T5INZJJc/sHGlqvrqJHdK8qaleXdK8unuvrmqdif5xiS/eDjBSDwBAHao7r6lqp6W5OIkxyV5cXdfWVXPS3JZd++/RdITk7y8u5fb8PdP8mtVtS+LLvnzl0fDTyHxBADYwbr7oiQXbZj37A3Tz9lkuzcmeeCRjEXiCQAwwVEyuGitHPUXJgAAsDOoeAIATFGVUvEcRcUTAIBZSDwBAJiFVjsAwBTrcR/Po8qos1VV51TVFVV1+fDczodvV2AAAOwsYyuer01yYXd3VT0oySuS3O/IhwUAcHRzO6XxRiWe3f3Jpcl/kAM/ZB4AAL7I6AsTquo7q+qvk/xxkh89wDrnDa34y66/4cbDjREAgB1g9OCi7n5lkldW1Tcl+dkkj9pknT1J9iTJ1z/4QaqiAMDOU2Vw0Uhbnq2qeuowmOjyqrrH/vnd/edJvrKqdm9rhAAA7AhbVjy7+/wk5ydJVd23qmoYXPSQJCckuWGbYwQAOCp5ctE4Y1vt353kSVX1+SSfSfL93a2VDgDAlsaOan9BkhdsUywAAOxgnlwEADBBlft4jmUoFgAAs1DxBACYwrPaR3O2AACYhcQTAIBZaLUDAExS7uM5koonAACzkHgCADALrXYAgAnKqPbRnC0AAGah4gkAMIUnF42m4gkAwCwkngAAzEKrHQBgkjK4aCRnCwCAWah4AgBMUBVPLhpJxRMAgFlIPAEAmIVWOwDARO7jOY6KJwAAs1DxBACYotxOaSxnCwCAWUg8AQCYhVY7AMAEVQYXjaXiCQDALFQ8AQAmKoOLRnG2AACYhcQTAIBZaLUDAExhcNFoKp4AAMxCxRMAYIJKZdcuFc8xtj3xvLWPy003n7jdh9kxTj3l+FWHsHbefcWHVh3C2rntPR+26hDWytXXfXrVIaydk070e38Mv/s5Vmi1AwAwC612AIApKtm1Sw1vDGcLAIBZSDwBAJiFVjsAwAQV9/EcS8UTAIBZqHgCAExRZXDRSM4WAACzkHgCADALrXYAgIk8MnMcFU8AAGah4gkAMEHFk4vGcrYAAJiFxBMAgFlotQMATFEGF42l4gkAwCxUPAEAJqhUdh2nhjeGswUAwCwkngAAzEKrHQBgCoOLRlPxBABgFiqeAAATeXLROM4WAACzkHgCADALrXYAgAmqyuCikVQ8AQB2sKo6u6reVVVXV9UzNln+w1X1d1V1+fB6ytKyc6vqPcPr3MONRcUTAGCHqqrjkpyf5NFJ9ia5tKou7O53blj1t7v7aRu2vXOSn0lyVpJO8pZh249OjUfiCQAw0RqMan9okqu7+31JUlUvT3JOko2J52b+SZLXdPeNw7avSXJ2kgumBnPUny0AACY7Nck1S9N7h3kbfXdVXVFVv1tVp4/c9pCpeAIATFFJHbfywUW7q+qypek93b1naXqzAHvD9KuSXNDdN1fVjyd5aZJHHOK2o0g8AQDW1/XdfdZBlu9NcvrS9GlJrl1eobtvWJr8H0lesLTtt2zY9g1TA0202gEAdrJLk5xRVfepqhOSPCHJhcsrVNUpS5OPS3LV8P7iJI+pqjtV1Z2SPGaYN5mKJwDABJU66gcXdfctVfW0LBLG45K8uLuvrKrnJbmsuy9M8hNV9bgktyS5MckPD9veWFU/m0XymiTP2z/QaCqJJwDADtbdFyW5aMO8Zy+9f2aSZx5g2xcnefGRikXiCQAwkScXjTOqPlxVPzgMtb+iqt5YVV+3XYEBALCzjK14/k2Sb+7uj1bVY5PsSfKwIx8WAAA7zajEs7vfuDR5SRbD6gEAjjlVa/HkoqPK4ZytJyd59ZEKBACAnW3S4KKq+tYsEs+HH2D5eUnOS5LTTjt9s1UAANaewUXjbFnxrKqnVtXlw+seVfWgJL+e5JwNd7r/e929p7vP6u6zTj755CMdMwAAa2jLxLO7z+/uM7v7zCwqpL+f5Ie6+93bHh0AADvG2Fb7s5OcnORXqipJbtni+aAAADtSVWXXcQYXjTF2VPtTkjxlm2IBAGAH8+QiAICJDC4aR30YAIBZSDwBAJiFVjsAwBSeXDSaswUAwCxUPAEAJqiUwUUjqXgCADALiScAALPQagcAmMLgotGcLQAAZiHxBABgFlrtAAATVJJdxxnVPoaKJwAAs1DxBACYosrgopGcLQAAZiHxBABgFlrtAAATeWTmOCqeAADMQsUTAGCCiornWCqeAADMQuIJAMAstNoBAKaopLTaR1HxBABgFiqeAAATGVw0joonAACzkHgCADALrXYAgAkqlSqt9jFUPAEAmIWKJwDAFGVw0VgqngAAzELiCQDALLTaAQAm8uSicWZIPCu7jvNDOVR162dXHcLauecZu1cdwtr59L1/ctUhrJWv+NS7Vx3C2rnpuK9cdQhrxe9+jhUqngAAE5TBRaO5xhMAgFlIPAEAmIVWOwDARAYXjaPiCQDALCSeAADMQqsdAGCSMqp9JBVPAABmoeIJADBFGVw0loonAACzkHgCADALrXYAgAkqya7Sah9DxRMAgFmoeAIATGRw0TgqngAAzELiCQDALLTaAQCmqHhy0UgqngAAzELFEwBggkoZXDSSiicAALOQeAIAMAutdgCAKQwuGk3FEwCAWah4AgBMZHDROCqeAADMQuIJAMAstNoBACaoJLtKq30MFU8AAGYh8QQAYBZa7QAAU5RR7WOpeAIA7GBVdXZVvauqrq6qZ2yy/N9U1Tur6oqqem1V3Wtp2a1VdfnwuvBwY1HxBACYpI76JxdV1XFJzk/y6CR7k1xaVRd29zuXVvurJGd196er6p8n+cUk3z8s+0x3n3mk4lHxBADYuR6a5Orufl93fy7Jy5Ocs7xCd7++uz89TF6S5LTtCkbiCQCwvnZX1WVLr/M2LD81yTVL03uHeQfy5CSvXpq+7bDfS6rq8YcbrFY7AMBER8Hgouu7+6yDLN8swN50xap/muSsJN+8NPue3X1tVX1FktdV1du7+71Tgx1V8ayq+1XVm6rq5qp6+tSDAgAwi71JTl+aPi3JtRtXqqpHJXlWksd1983753f3tcO/70vyhiQPPpxgxlY8b0zyE0kOu9QKALDOqnLUDy5KcmmSM6rqPkk+mOQJSX5geYWqenCSX0tydnd/ZGn+nZJ8urtvrqrdSb4xi4FHk41KPIdgPlJV3344BwUAYPt19y1V9bQkFyc5LsmLu/vKqnpeksu6+8Ik/ynJ7ZP8Ti0eAfq33f24JPdP8mtVtS+LLvnzN4yGH21brvEcLmw9L0lOO+30LdYGAGC7dPdFSS7aMO/ZS+8fdYDt3pjkgUcylm1JPLt7T5I9SXLmmQ/Z9AJWAIB1dxQMLlorWw4uqqqnLt2x/h5zBAUAwM6zZcWzu8/P4o73AAAs2VUqnmOMarVX1d2TXJbkxCT7quonkzyguz++HcEBALBzjB3V/qFs42OUAADYuTy5CABggqoyuGgkz2oHAGAWKp4AABOtwZOLjioqngAAzELiCQDALLTaAQCmKE8uGkvFEwCAWah4AgBM5MFF46h4AgAwC4knAACz0GoHAJigsnh6EYdOxRMAgFlIPAEAmIVWOwDARG7jOY6KJwAAs1DxBACYwpOLRlPxBABgFhJPAABmodUOADCR23iOo+IJAMAsVDwBACaolCcXjaTiCQDALCSeAADMQqsdAGAit/EcR8UTAIBZqHgCAExkcNE4Kp4AAMxC4gkAwCy02gEApihPLhpLxRMAgFmoeAIATFAxuGgsFU8AAGax7RXPffs6n/rULdt9mB3jg3s/u+oQ1s4DT7xs1SGsnTNv/2urDmGt/OUnX7bqENbOp27we38Mv/s5Vmi1AwBMVHrHozhdAADMQsUTAGCiXQYXjaLiCQDALCSeAADMQqsdAGAKTy4aTcUTAIBZSDwBAJiFVjsAwASV8sjMkVQ8AQCYhYonAMBECp7jqHgCADALiScAALPQagcAmMjgonFUPAEAmIWKJwDARLsUPEdR8QQAYBYSTwAAZqHVDgAwQZX7eI6l4gkAwCxUPAEAJiqji0ZR8QQAYBYSTwAAZqHVDgAwkcFF46h4AgAwCxVPAIAJqpJdSp6jqHgCADALiScAALPQagcAmEinfRwVTwAAZqHiCQAwUSl5jqLiCQDALCSeAADMQqsdAGAinfZxJlU8q+rsqnpXVV1dVc840kEBAHBkbJW3VdVtquq3h+Vvrqp7Ly175jD/XVX1Tw43ltGJZ1Udl+T8JI9N8oAkT6yqBxxuIAAAHFmHmLc9OclHu/u+SV6Y5AXDtg9I8oQkX5Pk7CS/MuxvsikVz4cmubq739fdn0vy8iTnHE4QAADrp1K12tchOJS87ZwkLx3e/26SR9Zi5+ckeXl339zdf5Pk6mF/k01JPE9Ncs3S9N5hHgAAR5dDydv+fp3uviXJx5KcfIjbjjJlcNFm6XV/0QpV5yU5L0lOPfW0CYcAADi6VXWqeusVt9fuqrpsaXpPd+9Zmt4ybzvIOoey7ShTEs+9SU5fmj4tybVfFNHiA+9Jkgc96MEr/4kAAOxQ13f3WQdZvmXetrTO3qr6siR3THLjIW47ypRW+6VJzqiq+1TVCVlcdHrh4QQBAMC2OJS87cIk5w7vvyfJ67q7h/lPGEa93yfJGUn+8nCCGV3x7O5bquppSS5OclySF3f3lYcTBADA2ukkfeuqozioA+VtVfW8JJd194VJXpTkf1XV1VlUOp8wbHtlVb0iyTuT3JLkqd2H94En3UC+uy9KctHhHBgAgO23Wd7W3c9eev/ZJN97gG1/PsnPH6lYPLkIAGCSTvW+VQexVjyrHQCAWUg8AQCYhVY7AMBUR/ngoqONiicAALNQ8QQAmKQTg4tGUfEEAGAWEk8AAGah1Q4AMFEZXDSKiicAALNQ8QQAmKLb7ZRGUvEEAGAWEk8AAGah1Q4AMJX7eI6i4gkAwCwkngAAzEKrHQBgEqPax1LxBABgFiqeAAATlcFFo6h4AgAwC4knAACz0GoHAJikkxhcNIaKJwAAs1DxBACYouPJRSOpeAIAMAuJJwAAs9BqBwCYxJOLxlLxBABgFiqeAAATeXLROCqeAADMQuIJAMAstNoBACYxuGisbU88jz/ulpxyh5u2+zA7xnV7FaHH2nvyw1cdwtp563V3X3UIa+XWfZ9fdQhr55Q7fGLVIawVv/s5Vqh4AgBMouI5lj+xAACYhcQTAIBZaLUDAExQ7T6eY6l4AgAwCxVPAIBJDC4aS8UTAIBZSDwBAJiFVjsAwFQGF42i4gkAwCwkngAAzEKrHQBgEqPax1LxBABgFiqeAAATVQwuGkPFEwCAWUg8AQCYhVY7AMAkBheNpeIJAMAsVDwBAKboeHLRSCqeAADLMQ7xAAASz0lEQVTMQuIJAMAstNoBACYxuGgsFU8AAGah4gkAMFEZXDSKiicAALOQeAIAMAutdgCASQwuGkvFEwCAWah4AgBM0Un2qXiOoeIJAMAsJJ4AAMxCqx0AYJJO+pZVB7FWVDwBAJiFiicAwCRtcNFIKp4AAMxidOJZVS+uqo9U1Tu2IyAAAHamKRXPlyQ5+wjHAQCwfvrW1b7WzOjEs7v/PMmN2xALAAA72LZc41lV51XVZVV12fXXy1EBANimUe3dvSfJniT5+gc/sLfjGAAAK9VGtY9lVDsAwDGoqu5cVa+pqvcM/95pk3XOrKo3VdWVVXVFVX3/0rKXVNXfVNXlw+vMrY4p8QQAmGR4ctEqX4fnGUle291nJHntML3Rp5M8qbu/JovB5f+lqk5aWv7vuvvM4XX5VgeccjulC5K8KclXV9Xeqnry2H0AALBy5yR56fD+pUkev3GF7n53d79neH9tko8kucvUA46+xrO7nzj1YAAAHFG7q+qypek9w1ibQ3G37r4uSbr7uqq668FWrqqHJjkhyXuXZv98VT07Q8W0u28+2D48MhMAYKrVDy66vrvPOtDCqvrTJHffZNGzxhykqk5J8r+SnNvd+4bZz0zyoSyS0T1JfirJ8w62H4knAMAO1d2POtCyqvpwVZ0yVDtPyaKNvtl6Jyb54yT/obsvWdr3dcPbm6vqN5I8fat4JJ4AAFOs/+2ULkxybpLnD//+4cYVquqEJK9M8j+7+3c2LNuftFYW14du+Th1o9oBAI5Nz0/y6Kp6T5JHD9OpqrOq6teHdb4vyTcl+eFNbpv0sqp6e5K3J9md5Oe2OqCKJwDAMai7b0jyyE3mX5bkKcP730zymwfY/hFjjynxBACYpJNe61b77LTaAQCYhYonAMBU+w776UHHFBVPAABmIfEEAGAWWu0AAFO0wUVjqXgCADALFU8AgEnW/slFs1PxBABgFhJPAABmodUOADBRG1w0ioonAACzkHgCADALrXYAgCm6PTJzJBVPAABmoeIJADCJ+3iOpeIJAMAsJJ4AAMxCqx0AYCr38RxFxRMAgFmoeAIATNEGF42l4gkAwCwkngAAzEKrHQBgkk7ak4vGUPEEAGAWKp4AAFMZXDSKiicAALOYoeJZya7jt/8wO8Tuu99u1SGsnb9970dXHcLa2f2Rn1t1CGtl37ddsOoQ1s4Juz676hDWit/9HCu02gEAJmlPLhpJqx0AgFmoeAIATNExuGgkFU8AAGYh8QQAYBZa7QAAk3Syz5OLxlDxBABgFiqeAACTuJ3SWCqeAADMQuIJAMAstNoBAKZwH8/RVDwBAJiFxBMAgFlotQMATGJU+1gqngAAzELFEwBgkk57ctEoKp4AAMxC4gkAwCy02gEAJuhO2n08R1HxBABgFiqeAACTtIrnSCqeAADMQuIJAMAstNoBAKZorfaxVDwBAJiFiicAwER9q4rnGCqeAADMQuIJAMAstNoBAKYwuGg0FU8AAGah4gkAMEF7ctFoKp4AAMxC4gkAwCy02gEAJtqn1T7K6IpnVZ1eVa+vqquq6sqq+lfbERgAADvLlIrnLUn+bXe/tarukOQtVfWa7n7nEY4NAODo5XZKo42ueHb3dd391uH9J5JcleTUIx0YAAA7y2ENLqqqeyd5cJI3H4lgAADYuSYPLqqq2yf5vSQ/2d0f37DsvCTnJck9T1cMBQB2Jq32cSZVPKvq+CySzpd19+9vXN7de7r7rO4+a/fJdz7cGAEA2AGmjGqvJC9KclV3/+cjHxIAADvRlFb7Nyb5oSRvr6rLh3k/3d0XHbmwAACObt2d3nfLqsNYK6MTz+7+iyS1DbEAALCDeXIRAMAk7uM5lme1AwAcg6rqzlX1mqp6z/DvnQ6w3q1VdfnwunBp/n2q6s3D9r9dVSdsdUyJJwDAsekZSV7b3Wckee0wvZnPdPeZw+txS/NfkOSFw/YfTfLkrQ4o8QQAmKIX9/Fc5eswnZPkpcP7lyZ5/KFuONzl6BFJfnfM9hJPAIBj0926+7pk8Uj0JHc9wHq3rarLquqSqtqfXJ6c5Kbu3j+sf28O4RHqBhcBAEzS6VtXPrhod1VdtjS9p7v37J+oqj9NcvdNtnvWiGPcs7uvraqvSPK6qnp7ko9vsl5vtSOJJwDA+rq+u8860MLuftSBllXVh6vqlO6+rqpOSfKRA+zj2uHf91XVG5I8OIsnWJ5UVV82VD1PS3LtVsFqtQMAHJsuTHLu8P7cJH+4cYWqulNV3WZ4vzuLBwm9s7s7yeuTfM/Btt9IxRMAYILFk4tW3mo/HM9P8oqqenKSv03yvUlSVWcl+fHufkqS+yf5taral0XB8vnd/c5h+59K8vKq+rkkf5XFI9UPSuIJAHAM6u4bkjxyk/mXJXnK8P6NSR54gO3fl+ShY44p8QQAmGjfelc8Z+caTwAAZiHxBABgFlrtAACTrP3gotmpeAIAMAsVTwCAKYZntXPoVDwBAJiFxBMAgFlotQMATNAGF42m4gkAwCwkngAAzEKrHQBgitZqH0vFEwCAWah4AgBM1PtuWXUIa0XFEwCAWUg8AQCYhVY7AMAU3elbDS4aQ8UTAIBZqHgCAEzQ6exzO6VRVDwBAJiFxBMAgFlotQMATOTJReNse+J5ax+Xm24+cbsPs2Ocesrxqw5h7bz7ig+tOoS1c9t7PmzVIayVq6/79KpDWDsnnej3/hh+93OsUPEEAJjCs9pHc40nAACzkHgCADALrXYAgIm02sdR8QQAYBYqngAAE3R39u3rVYexVlQ8AQCYhcQTAIBZaLUDAEyk0z6OiicAALNQ8QQAmKiVPEdR8QQAYBYSTwAAZqHVDgAwQXeyb9+qo1gvKp4AAMxC4gkAwCy02gEAJtrXRrWPoeIJAMAsVDwBACYyuGgcFU8AAGYh8QQAYBZa7QAAE3R7ZOZYKp4AAMxCxRMAYCIFz3FUPAEAmIXEEwCAWWi1AwBMtE+vfRQVTwAAZqHiCQAwQcfgorFUPAEAmIXEEwCAWWi1AwBM0e3JRSOpeAIAMAsVTwCAidxOaRwVTwAAZjE68ayq21bVX1bV26rqyqp67nYEBgDAzjKl1X5zkkd09yer6vgkf1FVr+7uS45wbAAARy338RxvdOLZ3Z3kk8Pk8cPLaQcA4KAmDS6qquOSvCXJfZOc391v3rD8vCTnJclpp51+uDECABx92uCisSYNLuruW7v7zCSnJXloVX3thuV7uvus7j7r5JNPPhJxAgCw5g5rVHt335TkDUnOPiLRAACwY41utVfVXZJ8vrtvqqrbJXlUkhcc8cgAAI5yvW/VEayXKdd4npLkpcN1nruSvKK7/+jIhgUAwE4zZVT7FUkevA2xAACwg3lkJgDABIv7eBrVPoZHZgIAMAsVTwCAKTrZZ3DRKCqeAADMQuIJAMAstNoBACboeGTmWCqeAADHoKq6c1W9pqreM/x7p03W+daqunzp9dmqevyw7CVV9TdLy87c6pgSTwCAifb1al+H6RlJXtvdZyR57TD9Rbr79d19ZnefmeQRST6d5E+WVvl3+5d39+VbHVDiCQBwbDonyUuH9y9N8vgt1v+eJK/u7k9PPaDEEwBgfe2uqsuWXueN2PZu3X1dkgz/3nWL9Z+Q5IIN836+qq6oqhdW1W22OqDBRQAAU3SnVz+46PruPutAC6vqT5PcfZNFzxpzkKo6JckDk1y8NPuZST6U5IQke5L8VJLnHWw/Ek8AgB2qux91oGVV9eGqOqW7rxsSy48cZFffl+SV3f35pX1fN7y9uap+I8nTt4pHqx0AYIL9t1Na5eswXZjk3OH9uUn+8CDrPjEb2uxDspqqqiyuD33HVgeUeAIAHJuen+TRVfWeJI8eplNVZ1XVr+9fqaruneT0JH+2YfuXVdXbk7w9ye4kP7fVAbXaAQCOQd19Q5JHbjL/siRPWZp+f5JTN1nvEWOPKfEEAJho9WOL1otWOwAAs1DxBACYoj2rfSwVTwAAZiHxBABgFlrtAAATdJLet+oo1ouKJwAAs1DxBACYaF8bXDSGiicAALOQeAIAMAutdgCAKTrZZ3DRKCqeAADMQuIJAMAstNoBACboJJ6YOY6KJwAAs1DxBACYSMVzHBVPAABmse0Vz7e97fLr73a3O35gu48zwe4k1686iDXjnI3nnI1zlJ6vZ606gIM5Ss/ZUcv5Gu9oPWf3WnUAjLftiWd332W7jzFFVV3W3WetOo514pyN55yN43yN55yN43yN55wdmMFF42m1AwAwC4OLAACmaBXPsY7liueeVQewhpyz8ZyzcZyv8ZyzcZyv8ZwzjpjqlqoDAIz1lXfa1c9/5PErjeH7fu9zb1mna3C12gEAJtJqH+dYbrUDADCjYzLxrKpzquqKqrq8qi6rqoevOqajWVX94HC+rqiqN1bV1606pqNdVd2vqt5UVTdX1dNXHc86qKqzq+pdVXV1VT1j1fEc7arqxVX1kap6x6pjWQdVdXpVvb6qrqqqK6vqX606pqNZVd22qv6yqt42nK/nrjqmo9Hidkq10te6OSYTzySvTfJ13X1mkh9N8usrjudo9zdJvrm7H5TkZ+NC80NxY5KfSPJLqw5kHVTVcUnOT/LYJA9I8sSqesBqozrqvSTJ2asOYo3ckuTfdvf9k3xDkqf6jh3UzUke0d1fl+TMJGdX1TesOCZ2gGMy8ezuT/YXRlX9gyz+aOEAuvuN3f3RYfKSJKetMp510N0f6e5Lk3x+1bGsiYcmubq739fdn0vy8iTnrDimo1p3/3kWf+BwCLr7uu5+6/D+E0muSnLqaqM6evXCJ4fJ44eX/1dy2I7ZwUVV9Z1J/mOSuyb59hWHs06enOTVqw6CHefUJNcsTe9N8rAVxcIOV1X3TvLgJG9ebSRHt6ET8ZYk901yfnc7Xxt1sm/fqoNYL8dkxTNJuvuV3X2/JI/Pon3MFqrqW7NIPH9q1bGw42x2oZLqCkdcVd0+ye8l+cnu/viq4zmadfetwyVppyV5aFV97apjYv0dM4lnVT11GEx0eVXdY//8oV31lVW1e4XhHXU2nq+qelAW18Ke0903rDq+o9GBvmMckr1JTl+aPi3JtSuKhR2qqo7PIul8WXf//qrjWRfdfVOSN8Q1xV9i/7PaV/laN8dM4tnd53f3mcNfb19eVZUkVfWQJCckkUwt2XC+vizJ7yf5oe5+94pDO2otn7PuljSNc2mSM6rqPlV1QpInJLlwxTGxgwy/81+U5Kru/s+rjudoV1V3qaqThve3S/KoJH+92qjYCY7Vazy/O8mTqurzST6T5PuXBhvxpZ6d5OQkvzLk67es01MSVqGq7p7ksiQnJtlXVT+Z5AFae5vr7luq6mlJLk5yXJIXd/eVKw7rqFZVFyT5liS7q2pvkp/p7hetNqqj2jcm+aEkb6+qy4d5P93dF60wpqPZKUleOlznuSvJK7r7j1YcEzuAR2YCAExwn5N29XO+6YSVxvDDr7p5rR6Zecy02gEAWC2JJwAAszhWr/EEADgs+0e1c+hUPAEAmIWKJwDARCqe46h4AgAwC4knAACz0GoHAJhiTR9buUoqngAAzELFEwBggk5yq4rnKCqeAADMQuIJAMAstNoBACYyuGgcFU8AAGah4gkAMEEn2bdv1VGsFxVPAABmIfEEAGAWWu0AAFN4ctFoKp4AAMxCxRMAYIKOiudYKp4AAMxC4gkAwCy02gEAJtJqH0fFEwCAWah4AgBMtK9r1SGsFRVPAABmIfEEAGAWWu0AABO4j+d4Kp4AAMxC4gkAwCy02gEApmit9rFUPAEAmIWKJwDABJ1k375VR7FeVDwBAJiFxBMAgFlotQMATGRw0TgqngAAzELFEwBgIhXPcVQ8AQCYhcQTAOAYVFXfW1VXVtW+qjrrIOudXVXvqqqrq+oZS/PvU1Vvrqr3VNVvV9UJWx1T4gkAMEEPTy5a5eswvSPJdyX58wOtUFXHJTk/yWOTPCDJE6vqAcPiFyR5YXefkeSjSZ681QElngAAx6Duvqq737XFag9NcnV3v6+7P5fk5UnOqapK8ogkvzus99Ikj9/qmAYXAQBMcN2n+uLn/p/P7l5xGLetqsuWpvd0954juP9Tk1yzNL03ycOSnJzkpu6+ZWn+qVvtTOIJADBBd5+96hi2UlV/muTumyx6Vnf/4aHsYpN5fZD5ByXxBADYobr7UYe5i71JTl+aPi3JtUmuT3JSVX3ZUPXcP/+gXOMJAMCBXJrkjGEE+wlJnpDkwu7uJK9P8j3Deucm2bKCKvEEADgGVdV3VtXeJP8oyR9X1cXD/HtU1UVJMlQzn5bk4iRXJXlFd1857OKnkvybqro6i2s+X7TlMRcJKwAAbC8VTwAAZiHxBABgFhJPAABmIfEEAGAWEk8AAGYh8QQAYBYSTwAAZvH/A+8mY0o0swAPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x864 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams['figure.figsize'] = 12,12\n",
    "plt.imshow(covariances[pixIdx[0],pixIdx[0],:,:], cmap=cmap, interpolation='Nearest', vmin=-1, vmax=1)\n",
    "\n",
    "plt.xticks(range(2*howFar + 1), range(-howFar,howFar+1))\n",
    "plt.yticks(range(2*howFar + 1), range(-howFar,howFar+1))\n",
    "plt.title('Partial covariances estimates, pixel {}'.format(pixIdx))\n",
    "plt.colorbar()\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 8,8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Annoying to match up, but if you do, they _do_ match.\n",
    "\n",
    "We need to figure out a way to make visualization less of a hassle, but at least the estimator does what it's supposed to do."
   ]
  }
 ],
 "metadata": {
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}