denoising_learned_prior.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Denoising with a learned DNN prior\n",
    "\n",
    "We consider the classical image or data denoising problem, where the goal is to remove zero-mean white Gaussian noise from a given image or data point. In more detail, our goal is to obtain an estimate of a vector $y_0 \\in \\mathbb R^n$ from the noisy observation\n",
    "$$\n",
    "y = y_0 + \\eta, \n",
    "$$\n",
    "where $\\eta$ is zero-mean Gaussian noise with covariance matrix $\\sigma^2/n I$, and $y_0$ lies in the range of the generator, i.e., $y_0=G(x_0)$.\n",
    "\n",
    "We consider the following two-step denoising algorithm:\n",
    "1. Obtaine an estimate $\\hat x$ of the latent representation by minimizing the empirical loss\n",
    "$$\n",
    "f(x) = \\|G(x) - y_0\\|_2^2\n",
    "$$\n",
    "using gradient descent.\n",
    "2. Obtain an estimate of the image as $\\hat y = G(\\hat x)$.\n",
    "\n",
    "We obtain G as the encoder of a learned autoencoder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.utils as utils\n",
    "from torch.autograd import Variable\n",
    "import torchvision.datasets as dset\n",
    "import torchvision.transforms as transforms\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import torch.nn.functional as F\n",
    "\n",
    "import random\n",
    "import numpy as np\n",
    "import collections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mnist_train = dset.MNIST(\"./\", train=True, transform=transforms.ToTensor(), target_transform=None, download=True)\n",
    "mnist_test  = dset.MNIST(\"./\", train=False, transform=transforms.ToTensor(), target_transform=None, download=True)\n",
    "\n",
    "batch_size = 1\n",
    "train_set = [ex for ex in torch.utils.data.DataLoader(dataset=mnist_train,batch_size=batch_size,shuffle=True)]\n",
    "test_set = [ex for ex in torch.utils.data.DataLoader(dataset=mnist_test,batch_size=batch_size,shuffle=True)]\n",
    "\n",
    "# construct training and test set only consisting of twos\n",
    "def extract_nu(dset,nu):\n",
    "    eset = []\n",
    "    for image,label in train_set:\n",
    "        if label.numpy() == nu:\n",
    "            eset.append((image,label))\n",
    "    return eset\n",
    "train_twos = extract_nu(train_set,2)\n",
    "test_twos = extract_nu(test_set,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specification of the autoencoder\n",
    "\n",
    "Next, we specify the encoder and decoder parameterized by lengthbottleneck:\n",
    "- Encoder: 784 neurons -> 400 neurons -> lengthbottleneck neurons\n",
    "- Decoder: lengthbottleneck neurons -> 400 neurons -> 784 neurons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Encoder(nn.Module):\n",
    "    def __init__(self,lenbottleneck = 20):\n",
    "        super(Encoder,self).__init__()\n",
    "        self.fc1 = nn.Linear(28*28, 20*20)\n",
    "        self.fc2 = nn.Linear(20*20,lenbottleneck)                \n",
    "    def forward(self,x):\n",
    "        x = x.view(-1, 28*28)\n",
    "        x = F.relu(self.fc1(x))\n",
    "        x = F.relu(self.fc2(x))\n",
    "        return x\n",
    "    \n",
    "class Decoder(nn.Module):\n",
    "    def __init__(self,lenbottleneck = 20):\n",
    "        super(Decoder,self).__init__()\n",
    "        self.fc1 = nn.Linear(lenbottleneck,20*20)\n",
    "        self.fc2 = nn.Linear(20*20,28*28)\n",
    "    def forward(self,x):\n",
    "        x = F.relu(self.fc1(x))\n",
    "        x = F.relu(self.fc2(x))\n",
    "        x = x.view(28,28)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(encoder,decoder,trainset,learning_rate = 0.001,epoch = 10):\n",
    "    parameters = list(encoder.parameters()) + list(decoder.parameters())\n",
    "    loss_func = nn.MSELoss() # mean square loss after the decoder \n",
    "    print(\"learning_rate: \",learning_rate)\n",
    "    optimizer = torch.optim.SGD(parameters, lr=learning_rate)\n",
    "    ctr = 0\n",
    "    for i in range(epoch):\n",
    "        print(\"at epoch \", i + 1 , \"/\", epoch)\n",
    "        for (image,label) in trainset:\n",
    "            ctr += 1\n",
    "            image = Variable(image)\n",
    "            optimizer.zero_grad()\n",
    "            output = encoder(image)\n",
    "            output = decoder(output)\n",
    "            loss = loss_func(output,image)\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "    print(\"trained in \",ctr,\" iterations\")\n",
    "    return encoder, decoder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training autoencoders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "learning_rate:  1.0\n",
      "at epoch  1 / 10\n",
      "at epoch  2 / 10\n",
      "at epoch  3 / 10\n",
      "at epoch  4 / 10\n",
      "at epoch  5 / 10\n",
      "at epoch  6 / 10\n",
      "at epoch  7 / 10\n",
      "at epoch  8 / 10\n",
      "at epoch  9 / 10\n",
      "at epoch  10 / 10\n",
      "trained in  100000  iterations\n"
     ]
    }
   ],
   "source": [
    "# train autoencoder for all digits\n",
    "K = 10\n",
    "encoder10 = Encoder(K)\n",
    "decoder10 = Decoder(K)\n",
    "encoder,decoder = train(encoder10,decoder10,train_set[:10000],1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the weight distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4000,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VNX5x/HPk4SEsEiAhEUggAooKmqNS7VVAwrBBQqIAuUnrvxc61atVooW7aLWn/sWBVFRUMFaqgh1QbEVEBRFAdEAApEl7DuBJM/vjzu00zBJJsmdObM879crr8zcOXPvl5vJw825954jqooxxpjEkuI6gDHGGP9ZcTfGmARkxd0YYxKQFXdjjElAVtyNMSYBWXE3xpgEZMXdGGMSkBV3Y4xJQFbcjTEmAaW52nB2drZ26tTJ1eaNMSYuff755xtVNaemds6Ke6dOnZg/f76rzRtjTFwSkZXhtLNuGWOMSUBW3I0xJgFZcTfGmARkxd0YYxKQFXdjjElAVtyNMSYBWXE3xpgEZMXdGGMSkBV3Y4xJQFbcjTEmAVlxN8aYBGTF3RhjEpAVd2OMSUA1FncRGSciJSLyTTVtzhKRL0VkkYh87G9EY4wxtRXOkft4oKCqF0UkC3gK6KeqRwOD/YlmjDGmrmos7qo6C9hcTZNhwJuquirQvsSnbMYYY+rIjz73rkBzEflIRD4XkUt8WKcxxph68GMmpjTgRKAXkAnMFpE5qvpd5YYiMhIYCZCbm+vDpo0xxoTix5F7MTBdVXep6kZgFnBcqIaqWqiqeaqal5NT4xSAxhhj6siP4v434OcikiYijYBTgCU+rNcYY0wd1dgtIyITgbOAbBEpBu4GGgCo6jOqukREpgMLgQrgeVWt8rJJY4wxkVdjcVfVoWG0eRB40JdExhhj6s3uUDXGmARkxd0YYxKQFXdjKtu3D5YuhfXrXScxps78uM7dmJiRn//fz2fODONN770HTzwBq1bBxImwYwcMHQqbN8OgQTB6NHToEJG8xkSKHbmb5DZrFlx+uVfEn3sOOnWCk06CoiL47js49FDYvRvWroWvv3ad1piw2ZG7SU4LFkBJCfTuDd98A82aHdwmOxt+/3vv8QcfwPDh8MtfwoMPgkh08xpTS3bkbpJLRQXcfz/06QPbt3tFOlRhr6xXL1iyBD780OumMSbG2ZG7SS6PPAKTJ8O8edCxY+3em5UF//gHfPFFZLIZ4yMr7ia5XHKJ17XSunWd3p4/OBvozfk3FJJRvpspHW4K76StMVFm3TImOSxdCv37Q4sWdS7swT5rUcCFxY9w3prnfAhnjP+suJvEt307/OIXcP75kOLPR76kYS63Hvc+I374PUyf7ss6jfGTFXeT2FTh0kvhjDPgqqt8XfWaRkdw8/EzvXUbE2OsuJvEJuIV9ccei8jqf2zUBYqLvStwjIkhVtxNwuqy4wsYOxb69oWMjMhtqFUr7z+PTz6J3DaMqSUr7iYxqXJd0U3ede0Rlj8gi7ubP8bKgpH0PrP0oCEQjHHBirtJSGdsmEKTsm3e0AJRMCt7IJ83P4fWe1dGZXvG1KTG4i4i40SkRESqnV1JRE4SkXIRudC/eMbUTY9tn/DkEQ9Damp0NijC410eY23DzjQvXRedbRpTjXCO3McDBdU1EJFU4H5ghg+ZjKmXBhWlPNHlURY07xn1bedveJ0xiwZFpTvImOrUWNxVdRawuYZmNwBTgBI/QhlTV81L1/HS3G40LNvpZPsftBqKUAGvvOJk+8YcUO8+dxFpBwwAnql/HGPq54oVo/io1WD2pjVxsn2VFMZ1vg/uuw/Ky51kMAb8GVvmEeA3qlouNQyDKiIjgZEAubm5PmzamCALFvDTTW9zySlLncb4IqsnPDchev39xoTgR3HPAyYFCns2cK6IlKnqW5UbqmohUAiQl5enPmzbmP9o1Yo/dH+FXWnVD+Eb6lJFXwf/EoETToCnn4aRI63IGyfqXdxVtfOBxyIyHng7VGE3JqK+/RYaNOCL5r1cJ/GkpsJLL3kDlV18ses0JgnVWNxFZCJwFpAtIsXA3UADAFW1fnYTG26/Hc49Fzi8Tm/3/WhexJvU47bbYPBg3wYsMyZcNRZ3VR0a7spU9dJ6pTGmLpYsgblz4bXX4DX/VlvvO00LCuCee2D2bDj9dD8iGRM2m6zDxL9HHoHrr4fMTNdJ/psIfPwxNGzoOolJQlbcTfz74x8hLUY/yg0bwlNPQW6uN568MVFiHYEmvr36KmzZEt4k1660awe/+503trwxUWLF3cSv7dvhhhti96j9gH79vBua3n/fdRKTRKy4m/hVWAh9+kCnTq6TVE8Ebr0VvvvOdRKTRGL8kMeYKqjCpEnwXGxOUH3wlTYj/L1RypgaWHE38UkEPv0U0tOrbRZTE2e89BIsX+5dHmlMhFm3jIk/qnDllbBjh+sktXPqqd6QBKWlrpOYJGBH7iYuBB+BH7VtLr/99mPaP9fCXaC66NoVjjsO3ngDhg93ncYkOCvuJu5csLaQt9uO5LWe1Y9CGpN+/WtYaVPxmciz4m7iimgFubuXUHjYn11HqZvevb3vZWWQlnbQOQE76Wr8Yn3uJq6opHD9T2azNb2V6yh1N2ECXH216xQmwVlxN/FDlbsXXUR26Y+uk9RPnz4wZQpsrmn2SmPqzoq7iRtH7fiMLju+YFN6W9dR6icnxxtn5oUXXCcxCcyKu4kb5615jnfaXolKAnxsb7kF2rRxncIksAT4LTFJQZWs/SVMb3Op6yT+OOEEGDaMhmU7XScxCarG4i4i40SkRES+qeL1X4rIwsDXpyJynP8xTdITYdSxU9mSkUBHu88/z03fX+c6hUlQ4VwKOR54AnipitdXAGeq6hYR6Ys3AfYp/sQzyeigIQNU+dPX5/F4l8dZk1m3afRiUv/+nH7NbWSW7WBPWlPXaUyCqfHIXVVnAVWe1lfVT1V1S+DpHKC9T9mMAaDrjs/puHsJaxt2rrlxPGnViq+yzuTMDZNdJzEJyO+bmK4A3vV5nSbJXbC2MHFOpFYyqcNtpFfsdR3DJCDfiruI5OMV959V02YkMBIgNzfXr02bJDAjUU6kVvJN1s9Iq9hnXTPGd74cColID+B5oL+qbqqqnaoWqmqequbl5OT4sWmT6FR5qFshGzPauU4SMVesGMVFqx9yHcMkmHoXdxHJBd4E/kdVbaoZ46tRS4Zx7NZPXMeIqA9aDaXP+hcRrXAdxSSQcC6FnAjMBrqJSLGIXCEiV4vIgcExRgMtgadE5EsRmR/BvCaJZO0r4ZRN7/J9kxNcR4mooibHszv1EHpsneU6ikkgNfa5q+rQGl6/ErjSt0TGBPQsmcSn2RewN62J6yiRJcKzh93P9gYtXScxCcSG/DVOVTcNXoqWM61t4hw3VPdvndeygCb7t8CePZCZGb1QJmEl3rVlJjGoMrnDzXyVdabrJFFz29IrvUm/jfGBFXcTky77YTQXrHnWdYyoeq/1cBg/3nUMkyCsuJuYI1pB73UvseiQn7qOElVzWp4H33xj0/AZX1hxNzGnx9ZZ7ErLYnmTHq6jRFVZSjoUFkKDBq6jmARgxd3EnMzynbze4VbXMdwYNAgq7Hp3U392tYyJLarMyT7fdQp3VOG002DaNDjmGNdpTByzI3cTU/JLXuOm765xHcOZ/J7CJC5mwnkTyc+v/vJJY6pjR+4masIpVOesn8BHrS6KfJgY9kGrYYxZNJCxne8DEddxTJyy4m4ioi5HnIfs38Sx2z7h3u4T/Q8UR4qaHM+ruXeSpvspk3TXcUycsm4ZEzOy9pXwRodbbehbEd5uexVZ+0pcJzFxzIq7iRmrGh3JS51Gu44RE9rs/YFnPz+RlIoy11FMnLLibmJCy9I1jJt3rHe1iGFdZmfWNezMiVs/cB3FxCkr7iYmnFXyOkub5tkJxCAftB7G2etfcR3DxCkr7iYm9CqZyIetqx1dOul82GoIs1sm8TX/pl6suBvn0ir28UPj7nye1ct1lJiyNb0Vn2QPgHXrXEcxcSicmZjGiUiJiHxTxesiIo+JSJGILBSRn/gf0ySyMmnAA0e+QEWKXZlbWc+SSfC//+s6holD4Ry5jwcKqnm9L9Al8DUSeLr+sUwyeXzB6bTfbdPvhvJpdj/46CPYutV1FBNnaizuqjoL2FxNk/7AS+qZA2SJSFu/AprE1mnnN7QqXc2PmUe4jhKTdqU1g1694K23XEcxccaPv4PbAauDnhcHlq31Yd0mwfUqmciHrYagYqd/qnT77ZBi+8fUjh/FPdS1ayEvVhaRkXhdN+Tm5vqwaRPvSlMy+Sg7uceSqdGpp8K2bbBzJzRJ8MnCjW/8OBwoBjoEPW8PrAnVUFULVTVPVfNycnJ82LSJa6pM6DSKZU2Oc50k9l13HUyY4DqFiSN+FPepwCWBq2ZOBbapqnXJmBr97/Lb6bPuRdcx4sPgwTZ5tqmVcC6FnAjMBrqJSLGIXCEiV4vI1YEm04DlQBHwHHBtxNKahJGi5fRa/ypLmp7sOkp8KCiAhQthTcg/io05SI197qpa7W2DqqrAdb4lMnGnLsP7Hrvtn2xvkM2qxkf5HygRZWTA03aVsQmfnYI3TmSW7WBK+xtdx4gvF18MqamuU5g4YbcEmuhL9nlS66q8HHr0gDlzoHNn12lMjLMjdxN1J255nzuWjHAdI/6kpsKAAfD6666TmDhgxd1EXc+SSXb5Y10NGQKvveY6hYkDVtxNVKVV7OP0jW8xM8duXKqTn/8cLrnE66IxphpW3E1Utdi3jg9aD2Njw/auo8Sn1FS48UbYXN1wT8bYCVUTZSUNc3m8y+OuY8SVypeadts+n2d2XwKLF9vMVaZKduRuoiajfDfPz+tBg/K9rqPEtaVN82DPHu+mJmOqYMXdRM0pm6axJb01+1Mbuo4S30S8E6s2HIGphnXLmKjpWTKJma2GuI6RGC67DBYtcp3CxDAr7iY6VClNzfTmBDX1160bHH64NxRws2au05gYZMXdRIcIfzrqZdcpEsuDD8L69fDII66TmBhkfe4mKkYtHsbR2z51HSOxDBzo3a1q17ybEKy4m4g7ZP8mTtn0Dssa93AdJbF06wZt2sCsWa6TmBhkxd1E3JkbJvNZi77sTbMp4nx3773QvLnrFCYGWZ+7ibj08j1Mb3uZ6xiJ6bzzYMcOKCuDNPt1Nv8R1pG7iBSIyFIRKRKRO0K8nisiM0VkgYgsFJFz/Y9q4tWUDjcxr0Uf1zESV79+8N57rlOYGBPONHupwJNAX6A7MFREuldqNgp4XVVPAIYAT/kd1MSngcWPceHqh13HSGwDB8LEia5TmBgTzpH7yUCRqi5X1X3AJKB/pTYKHBJ43AywiR4NAL3XvcTyxse6jpFw8vP/8zVw4mB2TPy7NySBMQHhFPd2wOqg58WBZcHuAYaLSDHehNk3hFqRiIwUkfkiMn/Dhg11iGviSbvd35NTWsyXzeswyaoJ25aMNozrfB/s3u06iokh4RT3UMPOaaXnQ4HxqtoeOBd4WUQOWreqFqpqnqrm5eTk1D6tiSs5pcVMbn8TFWLzfkbaW+2ugxS7+M38Rzin14uBDkHP23Nwt8sVQAGAqs4WkYZANlDiR0gTn75snm9H7VGSUb4bDjsMli+3SyMNEN6R+zygi4h0FpF0vBOmUyu1WQX0AhCRo4CGgPW7JLEjdizgnm8udB0jaZSmNoJevWDKFNdRTIyosbirahlwPTADWIJ3VcwiERkjIv0CzW4FrhKRr4CJwKWqWrnrxiSRXiUTKW7U1XWM5PLLX8Irr7hOYWJEWHc9qOo0vBOlwctGBz1eDJzubzQTr0Qr6Fkykd8c+67rKMnl3HNh2TJQtRmajA0/YPyXtX8DX2adxQ9NjnEdJblkZMCtt8KWLa6TmBhgxd34bkt6axve15W5cw+edNUkJSvuxlcNyvfy7PwTvas3TPSdfLJ35P71166TGMesuBtfnbbp7+xMy/Ku3jDRl5ICw4bZiVVjo0Iaf52z/mX+0eYS1zGS2+WXw8KFrlMYx6y4G9+IVlCa0ohZ2QNdR0luXbt6NzRt3QpZWa7TGEesW8b4RhHuPXoSe9Kauo5inngCbr/ddQrjkBV345sHFvah065FrmMYgEGDvLtVS0tdJzGOWHE3vuiweymH7fqaVZndXEcxAB06QI8eMG1azW1NQrI+d+OLc9a9zAethlGRYh+pmDF6NDRu7DqFccR+E40vVFKY0WaE6xgmWH4+7NzpzbHa1M6DJBvrljH1p8oLncewvEkP10lMZTffDIWFrlMYB6y4m3q76fvrOGODDTUbk4YPhxde8AYTM0nFumVMrQUPXZJevoc3SiYxoeNd7gKZqp1xBuzdC/Pnw0knuU5josiKu/m3UONNzZxZ/XvO2Pgm3zY9mY0ZlafVNTFBxLvmvUUL10lMlFm3jKmXVC1jcvubXMcw1SkogMxM7wjeJI2wiruIFIjIUhEpEpE7qmhzkYgsFpFFIvKqvzFNLErRcma0GcG8lgWuo5iaXHYZvPWW6xQmimos7iKSCjwJ9AW6A0NFpHulNl2AO4HTVfVowA7lksDlK0YxaPUjrmOYcIwY4Z1YNUkjnCP3k4EiVV2uqvuASUD/Sm2uAp5U1S0Aqlrib0wTa1IqyihYN555Lfq4jmLCMWAAzJsHq1e7TmKiJJwTqu2A4E9EMXBKpTZdAUTkX0AqcI+qTq+8IhEZCYwEyM3NrUteEyNO3TyNNQ0PZ1Xjo1xHMQHVnhDPzISXX/a+m6QQTnEPNdNu5Ytm04AuwFlAe+ATETlGVbf+15tUC4FCgLy8PLvwNg5UNWNbuaTxam7I0y8mVp13HixaBC1b2gTaSSCcbplioEPQ8/bAmhBt/qaq+1V1BbAUr9ibBNSgopTPWhQwJ/t811FMbajCRRfBJ5+4TmKiIJziPg/oIiKdRSQdGAJMrdTmLSAfQESy8bpplvsZ1MSOi1Y/xBUr7KaluCMCV10FzzzjOomJghqLu6qWAdcDM4AlwOuqukhExohIv0CzGcAmEVkMzARuU9VNkQpt3BGtoO/asXxisy3FpxEjvGGA1693ncREWFh3qKrqNGBapWWjgx4rcEvgyySw47d+xN7UJixtmuc6ignDwedMmjPzzTehSRMXcUwU2R2qplb2pjZmbOf77IRcPMvPhzlzoLzcdRITQVbcTdgalu2kqPFxzM6+wHUUUx8icOed8O67rpOYCLLibsI26MfHuGbZr13HMH649lp4+mnXKUwEWXE3YUmt2E+/NU/zTtsrXUcxfrj4Ypg7F1ascJ3ERIgVdxOWn218i7UNO7Os6fGuoxg/ZGZ6t6+2b+86iYkQK+4mLMWNulJ42P2uYxg/HXssTJ4MpaWuk5gIsMk6TI2yS39kQ3o7tqdnu45i/DZunHfn6rBhrpMYn9mRu6nRZStGc8Fam2Q5IV1zDTz1lOsUJgKsuJtqHbJvIz/f+CZvt73KdRQTCf36wQ8/wOLFrpMYn1m3jKnWeWuf55/ZA9iWnuM6iomEtDT47DNo29Z1EuMzK+6mWh+3GkyF/YGX2A491Lvm/dxzoWNH12mMT+y31lTp8B1fUkEK6zI7u45iIm3ZMnjEpkxMJFbcTZVuKPoVR+6Y5zqGiYYbb4QXX4QtW1wnMT6x4m5COmbbv8gpLWaWDe2bHDp08E6uzpjhOonxifW5m5AGr36ISbm/oSLFPiJJY+xYSE11ncL4JKwjdxEpEJGlIlIkIlVOnCkiF4qIiogN9h3nHupayPQ2l7qOYaIpNRVeeAFefdV1EuODGou7iKQCTwJ9ge7AUBHpHqJdU+BXwFy/Q5rounD1wzTQfexPyXAdxURbx45w771QUeE6iamncI7cTwaKVHW5qu4DJgH9Q7S7F3gA2OtjPuOD/PyDv6rSYde3DFv1J3alHhK9gCZ25OdDo0bwzjuuk5h6CqdDtR2wOuh5MXBKcAMROQHooKpvi4gN+B0Hqirww1b9mb+2u4G9aTYNW1ISgd/+FpYuhQtsUpZ4Fk5xDzWfmv77RZEU4GHg0hpXJDISGAmQm5sbXkITNenlezhqx1yePOJh11GMS4MGed/Ly+0EaxwLp1umGOgQ9Lw9sCboeVPgGOAjEfkBOBWYGuqkqqoWqmqequbl5Njt7LFmX2oml520iJ0NmruOYlybOxd69fJGjDRxKZziPg/oIiKdRSQdGAJMPfCiqm5T1WxV7aSqnYA5QD9VnR+RxCYiWpSu5b6v+6Mh/1AzSScvDzZssHlW41iN3TKqWiYi1wMzgFRgnKouEpExwHxVnVr9Gkw8GFz8f6xv2NHrczUJLdT5lpkzKy1ITYUxY2DUKCgogBS73zHehHWHiqpOA6ZVWja6irZn1T+WiabsvcX0XTuOK/O+ch3FxJKBA2HlSti3Dxo2dJ3G1JL9d2xoU7qSlzv+jo0NbT5NE0QEbrkFvvvOO7lq4ooV9yTXdP9mFh3yU6Z0uMl1FBOrrr/e7lqNQ1bck9yoxcM4e/0rrmOYWHbvvXDPPV73jIkbVtyT2Imb3+PQvcuY2epi11FMLDvzTDj8cG9IYBM3bMi/JCVawdXLbuO5zn+iLCXddRzjWI1X0IwdCy1aRC2PqT87ck9SivDUEf/HrJxBrqOYeNChA6xaBU895TqJCZMV9ySUXr6Hc9eOZUFWvl3XbsLXvDncfTcsWeI6iQmDFfckNKj4UU7d/I4VdlM7bdrA734H111nwxLEASvuSabtnuVcvPovPHP4X1xHMfHo2muhWTNYu9Z1ElMDO6GaZHqtf5WJub9hTebhrqOYeJSWBn/9K5SVwZ49kJnpOpGpghX3JJJasZ8JHe9CsD+pTT3dcw/s3Qt/sb8AY5UV9wRT1SQczfet57EFP+Pqn8xjV4Os6IYyiedXv4Kjj4ZLL4VjjnGdxoRgfe5J4rqim/hn9gAr7MYfrVp5R+/XXmsnV2OUHbkngVM3vs2R2z/jwZPGuo5iEsnVV8Oxx9pVVzHKjtyTwLImxzHm6NcpTW3kOopJJKmpcMYZ3o1N821unlhjxT3BDSp+lNKUTL5reqLrKCZRZWfDsGGwa5frJCZIWMVdRApEZKmIFInIHSFev0VEFovIQhH5QEQ6+h81ueXnH/xVkz5rx9Pvx6fZn5IR+YAmeV10EZx2Gtx8s+skJkiNxV1EUoEngb5Ad2CoiHSv1GwBkKeqPYDJwAN+BzW103nn11y9/DbuPmYKe9Kauo5jEt3jj3vXvJeVuU5iAsI5oXoyUKSqywFEZBLQH1h8oIGqBo8fNwcY7mdIE1p1R++5u7/lySMe4YfGR0cvkEleTZvCo4/CsmXelHzt2rlOlPTCKe7tgNVBz4uBU6ppfwVgU6a7ospxWz/m41aDXScxyei112D6dHj/fUi3oaRdCqfPPdR1TiEvbBWR4UAe8GAVr48UkfkiMn/Dhg3hpzRhG/DjE1y77BZSK/a7jmKS0R13eCdYr7rKrn93LJziXgx0CHreHlhTuZGInA3cBfRT1dJQK1LVQlXNU9W8nJycuuQ11Th626f8z8p7+f3Rb1Ce0sB1HJOMUlJgwgRYtw5WrHCdJqmFU9znAV1EpLOIpANDgKnBDUTkBOBZvMJe4n9ME46BxY/x5yNftEHBjFuNGnldM506wZdfuk6TtGrsc1fVMhG5HpgBpALjVHWRiIwB5qvqVLxumCbAG+LdrbZKVftFMHdCqXxi9L+mNwtDmz0r2J+Swb3dJ9rdgsY34VxuW+VnVQRWroTevWHKFPj5z33NZmoW1vADqjoNmFZp2eigx2f7nCuphfNLdUCL0rX85auzebnTaGa0GRG5UMbUVqdOXhfN4MHwySfQpYvrREnF7lCNY032b+GBhX2Y3uYyK+wmNvXuDffdB88/7zpJ0rGBw+JYl50L+KxFXyZ0vMt1FGOqduWV3pUzX33lddf06OE6UVKwI/c41HzfevquHceC5j0pPPx+62c3sU/Eu8Gpd2/47DPXaZKCFfc4k7trCU988VNySotdRzGmdgYO9Lpnzj8fZs1ynSbhWbdMHDls50Ie/OocCg9/wPrYTXw6/3yYNAkaN/a6auyvzoix4h4nGpbvojizC6OP+SuLmp3mOo4xddezp/d91CgoLYU//hEa2E13frNumRiXUb6bG77/FX/4+gL2pWZaYTeJ4+abYdEi79rfYutm9JsV9xjWZccXPPv5iTTbv5G7j57iOo4x/mrZEt5+G847D8baFJB+s26ZGJRSUYagZJTv5sWOdzOz9RDXkYyJjJQUuPNO7/Hs2TBxIowZA1k2kXt9WXH3SW3uKq3OUdvmcEPRjbzd9iqmHXqlPys1Jh507Qr79sFRR3n98CNGeMXf1IkV9xiRouX8bvFQum+fzfhOv2d6m0tdRzKmRqEOamo7NtK/tWwJzzzj3fR0//0wZIhX3DNsmsi6sOLuWIvStRyz/VNm5QzivdbD+fOR4ylNbeQ6ljF1Vu+Cn5cHb7zhPX7oIbj1Vl9yJRsr7o503fE5v/jxCX628S3+dui1zMoZxKfZNpCmSUx1LvhW2OvMinsUZZTv5pht/+LzFudw0ubprGp0FMNPfpDt6dmuoxljEowV9yg4adN0+q57gZM2z+CrrDNZ0Lwnr9hgX8aYCEqa4l7fCTHClVJRxhG7vuLYbf+kw+5veaTr0xxStpn5zc/h8S6PsSW9dWQ2bIwxQcIq7iJSADyKNxPT86r650qvZwAvAScCm4CLVfUHf6PGjgP/UYhW0HrvSjrtWsTPdD//zBnAvYsGcOie5XzT7HTmNe8NqnzQepjbwMaYpFNjcReRVOBJ4By8ybLnichUVV0c1OwKYIuqHiEiQ4D7gYsjETiSqjzps2cPzJsHq1bBypWctOlE5rUs4NU5h5FCOT80OprPWhQAMProN21yamOMc+EcuZ8MFKnqcgARmQT0B4KLe3/gnsDjycATIiKqqj5mrb3ycu97SgotS9eSUb6bzIpd7E1pRH5+F47b+jGH7llG47JtZFTsgY6/5acb/86g4kfJ2r+BrP0l3HjCa6zPyGXUkjtY37Aj6zM6sq9lJgCXnrz4oMsWrbAbU3f6MvKbAAAHAUlEQVS+Xjef5MIp7u2A1UHPi4FTqmoTmFB7G9AS2OhHyFp55x329LuI9IpShAoe7DaWf7S5hNc+z6M0pRF7Uxsxr0Ufnj38QY7Z9i/a7fmeXanN2JmWBaqsbNydibm/YWuDHLamt2JLg1ZUpKRxw08+PWhTdj26MXXn113dJrRwinuoAZcrH5GH0wYRGQmMDDzdKSJLw9h+KNmE+x/H0sth6eXkBC/b9TWs/gsfVW678vd1jFOHXNEXq9ksV+0kXa56DvmeiPurYziNwinuxUCHoOftgTVVtCkWkTSgGbC58opUtRAoDCdYdURkvqrm1Xc9fovVXBC72SxX7Viu2knmXOGMyjMP6CIinUUkHRgCTK3UZipwYGqgC4EPnfe3G2NMEqvxyD3Qh349MAPvUshxqrpIRMYA81V1KjAWeFlEivCO2G2MWmOMcSis69xVdRowrdKy0UGP9wKD/Y1WrXp37URIrOaC2M1muWrHctVO0uYS6z0xxpjEYyPhG2NMAorZ4i4ig0VkkYhUiEiVZ5VFpEBElopIkYjcEbS8s4jMFZHvReS1wMlgP3K1EJH3Aut9T0Sah2iTLyJfBn3tFZFfBF4bLyIrgl47Plq5Au3Kg7Y9NWi5y/11vIjMDvy8F4rIxUGv+bq/qvq8BL2eEfj3FwX2R6eg1+4MLF8qIn3qk6MOuW4RkcWB/fOBiHQMei3kzzRKuS4VkQ1B278y6LURgZ/79yIyovJ7I5zr4aBM34nI1qDXIrm/xolIiYh8U8XrIiKPBXIvFJGfBL3m7/5S1Zj8Ao4CugEfAXlVtEkFlgGHAenAV0D3wGuvA0MCj58BrvEp1wPAHYHHdwD319C+Bd5J5kaB5+OBCyOwv8LKBeysYrmz/QV0BboEHh8KrAWy/N5f1X1egtpcCzwTeDwEeC3wuHugfQbQObCe1Cjmyg/6DF1zIFd1P9Mo5boUeCLEe1sAywPfmwceN49Wrkrtb8C7ECSi+yuw7jOAnwDfVPH6ucC7ePcGnQrMjdT+itkjd1Vdoqo13eT076ERVHUfMAnoLyIC9MQbCgHgReAXPkXrH1hfuOu9EHhXVXf7tP2q1DbXv7neX6r6nap+H3i8BiiB/77vzCchPy/V5J0M9Arsn/7AJFUtVdUVQFFgfVHJpaozgz5Dc/DuN4m0cPZXVfoA76nqZlXdArwHFDjKNRSY6NO2q6Wqswhxj0+Q/sBL6pkDZIlIWyKwv2K2uIcp1NAI7fCGPtiqqmWVlvuhtaquBQh8b1VD+yEc/MH6Q+BPsofFG1Ezmrkaish8EZlzoKuIGNpfInIy3tHYsqDFfu2vqj4vIdsE9seBoTTCeW8kcwW7Au/o74BQP9No5hoU+PlMFpEDNzzGxP4KdF91Bj4MWhyp/RWOqrL7vr+cjucuIu8DbUK8dJeq/i2cVYRYptUsr3eucNcRWE9b4Fi8ewQOuBNYh1fACoHfAGOimCtXVdeIyGHAhyLyNbA9RDtX++tlYISqVgQW13l/hdpEiGXhDqVRr89UDcJet4gMB/KAM4MWH/QzVdVlod4fgVx/ByaqaqmIXI33V0/PMN8byVwHDAEmq2p50LJI7a9wRO3z5bS4q+rZ9VxFVUMjbMT7cyctcPQVasiEOuUSkfUi0lZV1waKUUk1q7oI+Kuq7g9a99rAw1IReQH4dTRzBbo9UNXlIvIRcAIwBcf7S0QOAd4BRgX+XD2w7jrvrxDqM5RGOO+NZC5E5Gy8/zDPVNXSA8ur+Jn6UaxqzKWqm4KePoc33PeB955V6b0f+ZAprFxBhgDXBS+I4P4KR1XZfd9f8d4tE3JoBPXOUMzE6+8Gb2iEcP4SCEfwUAs1rfegvr5AgTvQz/0LIORZ9UjkEpHmB7o1RCQbOB1Y7Hp/BX52f8Xri3yj0mt+7q/6DKUxFRgi3tU0nYEuwGf1yFKrXCJyAvAs0E9VS4KWh/yZRjFX26Cn/YAlgcczgN6BfM2B3vz3X7ARzRXI1g3v5OTsoGWR3F/hmApcErhq5lRgW+AAxv/9FamzxvX9Agbg/W9WCqwHZgSWHwpMC2p3LvAd3v+8dwUtPwzvl68IeAPI8ClXS+AD4PvA9xaB5Xl4s1QdaNcJ+BFIqfT+D4Gv8YrUBKBJtHIBpwW2/VXg+xWxsL+A4cB+4Mugr+Mjsb9CfV7wunn6BR43DPz7iwL747Cg994VeN9SoK/Pn/eacr0f+D04sH+m1vQzjVKuPwGLAtufCRwZ9N7LA/uxCLgsmrkCz+8B/lzpfZHeXxPxrvbaj1e/rgCuBq4OvC54kx8tC2w/L+i9vu4vu0PVGGMSULx3yxhjjAnBirsxxiQgK+7GGJOArLgbY0wCsuJujDEJyIq7McYkICvuxhiTgKy4G2NMAvp/bNa+AfIQQEMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117b21978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = decoder.fc1.weight.data.numpy()\n",
    "a = np.ndarray.flatten(a)\n",
    "print(a.shape)\n",
    "import matplotlib.mlab as mlab\n",
    "n, bins, patches = plt.hist(a, 50, normed=1, facecolor='blue', alpha=0.75)\n",
    "\n",
    "# add a 'best fit' line\n",
    "mu, sigma = 0, 0.26\n",
    "y = mlab.normpdf( bins, mu, sigma)\n",
    "l = plt.plot(bins, y, 'r--', linewidth=1)\n",
    "\n",
    "plt.axis([-1.1, 1.1, 0, 1.75])\n",
    "\n",
    "#plt.hist(a, normed=True, bins=1000)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "learning_rate:  1.0\n",
      "at epoch  1 / 10\n",
      "at epoch  2 / 10\n",
      "at epoch  3 / 10\n",
      "at epoch  4 / 10\n",
      "at epoch  5 / 10\n",
      "at epoch  6 / 10\n",
      "at epoch  7 / 10\n",
      "at epoch  8 / 10\n",
      "at epoch  9 / 10\n",
      "at epoch  10 / 10\n",
      "trained in  100000  iterations\n"
     ]
    }
   ],
   "source": [
    "# train autoencoder for all digits\n",
    "K = 20\n",
    "encoder20 = Encoder(K)\n",
    "decoder20 = Decoder(K)\n",
    "encoder20,decoder20 = train(encoder20,decoder20,train_set[:10000],1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# given two lists of images as np-arrays, plot them as a row\n",
    "def plot_images(top,bottom):\n",
    "    fig, axes = plt.subplots(nrows=2, ncols=10, sharex=True, sharey=True, figsize=(20,4))\n",
    "    for images, row in zip([top, bottom], axes):\n",
    "        for img, ax in zip(images, row):\n",
    "            ax.imshow(img, cmap='Greys_r')\n",
    "            ax.get_xaxis().set_visible(False)\n",
    "            ax.get_yaxis().set_visible(False)\n",
    "    fig.tight_layout(pad=0.1)\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABawAAAEsCAYAAAAvofT2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XmczmUb//9zso99GfsytsiWyJLIUpGUJSWyxV1KSWnBTVmLVlp0p1QqpaS6kRLJVpE1W4juLHMjy9gNxjK/P37fvt/7PI+D6+NzXdfMZ2Zez/+Oo2M+11nzmWs5ux7nOyYlJcUAAAAAAAAAAJDWrkjrBQAAAAAAAAAAYAwb1gAAAAAAAACAgGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACgQ1rAAAAAAAAAEAgsGENAAAAAAAAAAgENqwBAAAAAAAAAIHAhjUAAAAAAAAAIBDYsAYAAAAAAAAABELWyxkuUqRISnx8fJSWgvRu9erVB1NSUuIu9s+5f3Ax3DsIB/cPwsH9g3Bw/yAc3D8IB/cPwsH9g3Bw/yAcoe6fv13WhnV8fLxZtWqV/1UhQ4uJidl5qX/O/YOL4d5BOLh/EA7uH4SD+wfh4P5BOLh/EA7uH4SD+wfhCHX//I0jQQAAAAAAAAAAgXBZ37D+XzExMZFcB9KplJQUXz/H/QNjuH8QHu4fhMPP/cO9A2N47kF4uH8QDu4fhIP7B+Hg/kE4/Nw/fMMaAAAAAAAAABAIbFgDAAAAAAAAAAKBDWsAAAAAAAAAQCCwYQ0AAAAAAAAACAQ2rAEAAAAAAAAAgcCGNQAAAAAAAAAgENiwBgAAAAAAAAAEAhvWAAAAAAAAAIBAyJrWCwAAAEDay5Ili1X/+uuvYqZQoUKiV7p06aitCQAAAEDmwzesAQAAAAAAAACBwIY1AAAAAAAAACAQ2LAGAAAAAAAAAAQCG9YAAAAAAAAAgEAgdBEAAADm8ccft+oaNWqImQsXLqTWcgAAAABkUnzDGgAAAAAAAAAQCGxYAwAAAAAAAAACgQ1rAAAAAAAAAEAgsGENAAAAAAAAAAgEQhfTsSuukP+/YfHixaLXvHlzqz537lzU1gQACL4nn3xS9Fq0aGHVrVu3FjMnTpwQvZdeesmqR40aFebqkFY6d+4ccuaPP/5IhZUgo5o0aZLo3XfffVY9e/ZsMXP77bdHbU0AMh7tfc6wYcNEL2/evFb92GOPiZkyZcqI3pIlS6x61qxZl7tEAGmoQYMGotetWzdf1+rRo4dV58uXT8wsXLhQ9ObPn2/VY8aM8fX4GRnfsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCJxhfRHdu3e36m+//VbMJCYmptZyVNo5OI0bNxa9ZcuWWXW9evWitiZEVq5cuURv8+bNole2bFmrPn78uJjp06ePVU+bNi3M1SFotm/fbtXr168XM+3atUut5SCN1KhRw6p//PFHMaOdrRYTE2PVKSkpYiZ37tyiN3z4cKt+6qmnxMxNN91k1cuXLxczSHsVK1YMObNhw4ZUWAkygsmTJ4ue+/7aGPlcc+utt4qZwYMHi97zzz8fxuqQUbhnDj/99NNipnDhwqKXlJRk1S+//LKYcV/fEFzFixe3au13577PMUbmbtSsWVPM9O7dW/Qefvhhq7777rvFDOda+1egQAGr1vY47rzzTtHLmtXe3rrrrrvEzIIFC0SvVatWIdek3T87d+606kWLFokZbR/JvTdOnz4d8vGhc/PajJHvIxo1aiRmGjZsKHra79gP7TNUs2bNQvbcz0vG6J+rVq9e7Xtt6Q3fsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACgdDFixgzZoxVawfx9+zZU/TOnz8ftTXdfPPNVq2FD2iuvfbaaCwHqWDmzJmiV6ZMGdFzD/bPkyePmGnQoIFVew1ddO8f7X767LPPRO/IkSOerp+ZZMmSRfS+++470fvpp5+seuTIkZ6uv2TJEqvWghtKliwpenv27PF0fQSPFg7k3gf58+eP6hrcgBIvwYxaqBpSV69evUQvb968IX/OfX4C/nbllVdadadOncTMFVeE/q6MFnqkhVAj83njjTdE76GHHrJq7f7RArDce+qZZ54RM/Xr17fq1q1be1onUt/HH39s1dpzhva6N2XKFKuOj48XM2fOnBG9bt26WbUW9knoon/Lli2z6ipVqnj6Ofdv/fDhw2KmXr16onfq1KlL1hfjfubWXvd69Ogheu57Ke255eTJk57WkNm4n4sHDRokZrJnz55ay4koLZhxzpw5ole0aNFUWE0w8A1rAAAAAAAAAEAgsGENAAAAAAAAAAgENqwBAAAAAAAAAIHAGdbGmP79+4teXFycVXft2lXMjB8/XvRWr14duYU5PvnkE6v2cg6gMcZs2LAhGsvBJbjnpmnnBs+bN0/0ihUrZtWxsbFiRjuHz4uEhISQM+45gMYYM27cOKvOli2bmNHW6f4c9P9O2jnT7hnAXs+w7tevn1WvXbtWzLRt21b0Jk6c6On6CJ5vv/1W9KJ9ZrUfbg7EjBkzxEz79u1TazmZjnaWp3Zeq/u+Qnu9cc9IB/7WvXt3q/Z77nRSUpLovf/++76uhfRrwoQJote3b1/R086sjpRGjRpF7drw74YbbgjZ27Ztm5hxz6vW7NixQ/Qefvhh0UtOTrbqRx99VMy4Z2ZPnjw55OPj/zd48GCrvvHGG8VMuXLlRM/Nf/L62nH11Vdb9bp16zz9nOvtt98Wvfvvv1/03GwpLZts7ty5vtaQ0Q0dOtSqve6JuQ4ePCh67tnpGjfrzhhjtm/fbtXamesvvvii6Ln3gaZw4cKi98gjj1i1lu+QUfANawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBwIY1AAAAAAAAACAQMl3oYvny5UVPOwA9R44cVp2YmChmfv/998gtzHHLLbeInnbguuvAgQOiV7t27YisCTrtnpo/f75Vx8fH+7q2FoZXrVo10bv++uutWgsMcIM+br75ZjGj/S1oIYsu7e8D0vnz50XvzJkzoqcF6Xlx/Phxqz527JiYKVSokK9rI+29++67oleiRImIXf/kyZNW/dNPP4mZhg0bip6XkEc3FEt7jbv22mtFb9WqVSGvjdCKFi0qetprlxuyqD2H/PXXX5FbGKA4cuSI6O3cuTMNVoJo0d7LduzY0aq1wPtoBixq8uTJY9VPPPGEmHnllVdSazn4Pz777DPRy5rV3tZwQ8kibdiwYVbdu3dvMdO/f3+rJnTROzc80a0jzW/IoksL29Ps37/fqglYjCx3T+ymm24SM7t37xa9Q4cOReTx3d+vMcY0bdpU9FavXm3V1atXFzPa654bnE7oIgAAAAAAAAAAUcaGNQAAAAAAAAAgENiwBgAAAAAAAAAEAhvWAAAAAAAAAIBAyPChi27Y3cKFC8WMG7BojDGnT5+26rZt24qZEydOhLe4/5EzZ06rnjJlipi54orQ/3+he/fuonfhwgX/C4OlQYMGoqcF5BUoUCDktbT7xw09fO6558SMFlSzZs0aq964caOYWbFihVWXLl1azGh/C66xY8eK3rRp00L+HIypWLGi6HkJtfRLC4Zt0aKF6I0ZMyZqa4B/bsjivffeK2a8vC5oZs2aJXpuQFFCQoKY0UIX58yZY9VeQhizZ88uegsWLBC95s2bW7UbTgL/3IBFzUcffSR6e/fujcZykAF06tTJ18+571NHjRolZrTQIfe986lTp3w9PqJLCwy+5557RM/9fUbSvn37RG/Pnj1Wfc0114gZ977T3kMRuhhdxYsXFz0tQHzr1q1WPX/+/KityRgZdK4F3rdq1cqq69evL2bcz2dIX9zntyZNmqTRSjIPN/xWew7+5ZdfrHrDhg1RXZMX2v6Ql89MGi0wMqPiG9YAAAAAAAAAgEBgwxoAAAAAAAAAEAhsWAMAAAAAAAAAAiFDnWFdvnx50XPPuyxYsKCna73++utWvXTpUv8L82D9+vVWXaRIkZA/89dff4medgYoIqd69eqi5+We0s4ZGj58uOi9//77Ia+lnR+6c+dOq37wwQdDXkc7D9I9u90YY1577TWr/uSTTzz9HKRnnnlG9M6fPy96X3/9dUQeL1++fKJXuHDhiFwbkdWrVy/Rc8/49Hte9W+//SZ62vOPdma1yz0TzhhjWrdubdXumdbGeDujLU+ePKJ39913WzVnWPtz/fXXi572GuCea62dgw9cTIUKFXz9nHuW8DvvvCNmtPON161bZ9VaHsMHH3zga03wb/v27VZdpkwZMeP39cwv7XPVv//9b6vWzrB23XjjjaJXu3Zt0Vu7du1lrA6X4n4mN0bPwdByf1KTlmnUsmVLq9ZygLR7CsHkZr0YY0zXrl2tWntvlZycLHqPPfZY5BaWybh7E9rrifYckZqKFSsmegMHDhQ9LVPMC20/JqPiG9YAAAAAAAAAgEBgwxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAIKTb0MX27duL3meffSZ6OXLkCHmt2bNni97o0aP9LcyDJ598UvS8BNWcPXvWqq+77rqQMwhPiRIlrHrcuHGefu7YsWNWXbNmTTFz5MiRkNfRAh2XLFkiepUqVbJqNzhLo933I0eOFD1CzvxzQ35uuOEGMbN161bRc4OkvHJDJ8qWLStmtDArpL0JEyaInhYy5sWWLVusukGDBmImKSnJ17U1bhDjlClTxEy/fv18XbtWrVq+fg62pk2bip6X1wkgyNz3Pu+9956YcUPTtddc+Ddp0iTRc0MWUztgUXPu3DnRW7NmjVUfPXpUzLiBwVqQ1x133CF6hC76V6BAAau+8847xYz7OcsYPfA5NWlrcsXHx0d/IfAlW7ZsVv3WW2+JGS0g3Q1Z1J5rOnbsKHrffPPN5S4RFzF+/PhUfbxBgwaJ3qOPPmrV2h6ktq/jxfHjx0Vv8uTJvq6VHqX9OwgAAAAAAAAAAAwb1gAAAAAAAACAgGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACId2ELlasWNGqtWAVLwGLWijhhx9+KHq5c+e26hMnToS8tqZZs2ai9+yzz4pelixZQl5rxowZVr1jxw5fa4J3hQsXtmo3kOFi3EAzLYBBO3j/xhtvtOrnnntOzLghQ5pFixaJ3ubNm63abwgavHODzvLmzStm5s6dG7HHc8NcqlevLma0AFCkLi1g0cvrl1fu61UkAxa9GDhwoOi5IcF169b1dC33ORH+aEHVmuTkZKueNm2ar8fr3bu36A0bNkz03LAiNwTNGBmo7XdNSN/Onz8vej/++KNVL1y4UMz8+eefUVtTZuOGShmjh5CldciiGwRsjDF9+vQRvY0bN1p1t27dxIwWlo3oGjt2bMiZp59+WvQOHjwYjeV49sEHH4ieFkqK1OWGwBpjzIsvvih6LVu2tGqvAXluwKoWyPf99997uhaCR3v/PGrUKNHzukfkR2xsrOi5z4FvvPGGmMko73/4hjUAAAAAAAAAIBDYsAYAAAAAAAAABAIb1gAAAAAAAACAQAjkGdYPPfSQ6LlnxRQqVMjXtbXzZaZPny567pmfe/fuDbkmY4xZtWqVVc+aNUvMeDmr1D1H0hj9TCREV/78+a3ay1njxhhTrFgxq+7Zs6eY0c6nds84ds/3NMaYkydPip57Zuxbb73laZ2IrptuusmqteefAQMGiJ57n+XLl0/MtGvXTvS8nE997733it6RI0es+qmnngp5HXjXv39/q+7bt6+Y0f7W06vTp0+LXmJiYkSuXblyZdHbtm1bRK6dkRQoUMCqtecQTUpKilV7PRPUzeYYMmSIp59zaWdNulkA2bNnFzNTpkzx9XhIP7QMGvfeQHTFxcWJnt/zqmfOnCl6r732mlVr58xee+21ord7926rbtOmjZg5fPjw5S4RqSBnzpyid/vtt4f8OTfbIKh+/vlnq9Ze4xBZ7rn6Y8aMETPu5/Rw3HbbbVa9Z8+eiF0baU9773HgwAHRK1myZNTWoO0/uZkS2ucsLTcvPeIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACgQ1rAAAAAAAAAEAgBCJ08ZFHHrHql19+WcxoITvRFBsba9UVK1YUM++9957oucFZXkP6XC1atBC97du3+7oW/HPDMrTQFi/BDa+//rqvx9fCyzp06CB68+fP93V9RFfRokWteuvWrWJm5MiRotejRw+r1p7/tABHL883WkCSFuSJyHFDdiIZsOgGBBujv4amV+79qoUxjRs3LrWWk2506tTJqr2EPRtjzNKlS0PO9OnTR/Tc8BfNmTNnRG/lypVWfeWVV4oZ93n0+eefFzOELgLRp/19eqG9x3j88cdFz/2c07p1azFTv3590Tt06JBVew1YdD/baYGOrvPnz4verl27PD0eJC0Q2A0vW79+vZjxGggcNBkpYDuo6tSpY9WRDFjUfPfdd1bdsWNHMUM4ePr1zTffiJ57jxljTJMmTSLyeA8//LDoNWvWLOTPlS9fPiKPH0R8wxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIhECELroBUX4DFjds2CB6bhjd3LlzxUy9evVEzz243A1BM8aYrFmj959PCxBavHixVSckJIiZESNGiN6FCxcitq7M7sEHHxS9GTNmROz6bshi8+bNxczy5csj9niILvf5pn379mKmePHiouc+t0ybNk3MFChQQPT27t1r1dr9qt0/2vMG/Clbtqzode/ePSLX3rNnj+gNGDBA9KZPnx6Rx/NLCwe57rrrfF3Lff0iYDG6tHAr1/Dhw0Uvd+7cVp2cnCxmSpUqJXpuWFrNmjXFzLp166w6Li5OzGhBbCtWrBA9RE7fvn1FTwv1RcbRtm1bXz+3atUq0fMSJJ+YmCh6c+bM8bUGjfv5LzY2NuTPaM9t7777bsTWlNncf//9IWeWLVuWCiuJjt9//92q3RBuRJ6X94kNGzYUvalTp1p106ZNxYwWtlejRg2r1t57aM8RTz31VMh1Ipj2798vel9++WVErq29D3ffBxtjTM6cOa26Z8+eYubnn3+26vfffz/M1aUN3lkCAAAAAAAAAAKBDWsAAAAAAAAAQCCwYQ0AAAAAAAAACIRAnGHt5YzllJQUq/7www/FzMMPPyx6SUlJIa+tnWvt0s7PmjRpUsif88s9Q1vr/fnnn2Jm1KhRoscZ1t7kypVL9IYOHWrVffr0ETPuvXmxnhe7du2yas6rTt/c5wjt7MXdu3eLXkxMjFVrf8PaPXbbbbdZtXaG9cGDB/XFIiK0szqLFCni61q//fabVWtn7p08edLXtSPJPbN65syZYiZv3ryptBqEY8KECVb9yiuviBntfj5z5oxVa+e2u+dVa9wcB42WH3L11VeLHmdYR1ft2rVFz33tAoLskUceueyfOXz4cBRWknl16dIl5MzgwYNTYSXR4b4v9/Iah/C45+P7+Ts3xpjx48eLXqFChUTvvffes2rtrH8tb8Z9vXzyyScvd4nIgLZt2yZ62ln/bt6dliEybNgwq+YMawAAAAAAAAAAwsCGNQAAAAAAAAAgENiwBgAAAAAAAAAEAhvWAAAAAAAAAIBACEToYt++fa16yJAhYuaBBx6w6oULF0Z1TSVKlLBqLXjIr4kTJ4qeG1q1adMmMfPRRx9Z9Q8//CBmzp07F97iMrEFCxaIXv369VN1DW6w5vDhw8XMyJEjU2s5CJMbjJiQkBDVx8uePXvImeTk5KiuIbOLi4sTPb8hrKNHj7bq1A5YjI2NFT0thKZz585WHcmAxcWLF0fsWrBpAXnTpk2z6vPnz4uZbNmyiZ4bGPzFF1/4WpOXcK2zZ8+K3po1a3w9HoDMQXs9y5cv32VfR/sMB/+034vr6NGjqbCS8PXq1Uv03P2EV199NbWWgyjQwqM7dOhg1f369RMzr7/+uug99thjVr106VIx89VXX13uEpEBrVy50tfPuc8/Wmj6wYMHfV07NfENawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBwIY1AAAAAAAAACAQAhG6+MEHH1yyTgtTp061aj/BHMbI0Cxj9NA8LdgI0dOnTx/Rq1u3rq9raQFQL7zwglXff//9YqZo0aKilzWr/SfpBjIYY8wbb7wheloIBDKftWvXhpwpUKBAKqwEkVCjRg2r/vzzzyN27bJly4re0KFDrbpBgwZiplatWhFbg+vUqVOip4XhIrR9+/ZZ9YULF8TMFVfI7yzUqVMnamvSwjjd17iePXuGvM7+/ftFb/Xq1f4XhjTnvseuVq2amNHCyAGNdv+4gbLGGFO9evWQ1zp9+rRVf/nll/4XhgylYcOGVq19PnMDiQcNGhTVNSHtvfnmm6J38803i97tt99u1aNGjRIzhC5GToUKFUQvISFB9LR9nfTKDUkfOHCgmNF6QcM3rAEAAAAAAAAAgcCGNQAAAAAAAAAgENiwBgAAAAAAAAAEQiDOsE5r7rmdxhhzww03+LqWexbOuHHjxAznVae9AQMGiF6WLFl8XevBBx8UvcmTJ1v1v/71LzHz3XffiV7NmjWtOn/+/GJGuzdnzJgRcp2AMZz1Gm1nzpwRvezZs/u6lvvcUr9+fV/X0Vx//fWilzt37ohdP5QTJ06IXvPmzUWP+9WfmTNnWvWQIUPEzPPPPx+xxytVqpRVJyYmihntNdZLPkhKSopVazkgiC73HERj5PuVcLj3QYsWLcQMZ1inX7ly5RI97Z7ye3aoez7sSy+9JGa8nFedlJQkeo8++qhVb968+TJXh0uJiYlJ6yV4on32mjNnjlVrn+9btmxp1RnpfFzotHs6T548IX9OyxWBfx07drTqTz/9VMxo542/+uqrVr1nzx4xk17+js+dO2fV6TWDgb8MAAAAAAAAAEAgsGENAAAAAAAAAAgENqwBAAAAAAAAAIHAhjUAAAAAAAAAIBAyXehikSJFRE8LI/Jy8L0bsGiMMfXq1bPqI0eOXMbqEC1lypSx6pIlS4oZLSTBPVT/7bffFjNuwKJGCzMrXbp0yDXs27dPzBCwiHDUrVs3rZeQod1yyy2i5wbz5MiRw9O13NcrN7wnqNyQPGNkyOzw4cPFDAGL0aMFy7Rp00b0GjVqZNVew4jduYIFC17G6v6fY8eOiZ4bsvjuu+/6uja8i42NteoOHTqImYYNG0bs8ZKTk61ae++DtLd06VLRa9asWcif0wKDX3vtNdFzAw61YKucOXOK3ltvvWXVFSpUCLkmjRtWa4wx7733nq9rwRvt/UJaq1y5suj961//Ej33eatJkyZiZuvWrZFbGNKFUaNGiZ4WKu5Kr4F4QVW+fHmrvnDhgphxX3O03rJly8TMkiVLRM8Na4yk4sWL+/o59995w4YNkVhOquMb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACIdOdYT19+nTRc8/q8+quu+4SPc7dC6Ybb7zRqvPkySNmtHPU5s+fb9UTJ04UM3379g15rREjRogZ7YzPxMREq37llVfEDHAx7pn5SUlJabSSzGvRokWi5/4da7kJ6YX2POned2PGjBEzPJelrZMnT4pe06ZNRW/SpElWXbt2bTGjnYPvPtf88MMPntY1bdo0q164cKGY2bt3r6drIXLcsxjvu+++qD7eli1brFp7r460p/1da88jWiaM68EHHxS9WrVqWbX7ntgYY0qVKiV6fs+sPnz4sFWPHTvW13Xgn/u3b4wx8fHxVq3lLXzzzTcRW8Pjjz9u1f369RMzBQoUED33XOKNGzdGbE0IJjcTyxj5euXmmV2Mewb6ggUL/C8Mwssvv2zVWj7d888/H/I61113nafeoEGDLmN1qSN79uxWre1HDRw4MJVW4x/fsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACIcOHLjZr1syqtXAQL7QgoJUrV/q6FlLfwYMHff1cnTp1rPqpp54SMz179hQ9LZjMdeHCBdFLSEiwajcwALiU06dPW/XZs2fTaCX4X08//bRVa88PTzzxhOjlzJkzamvSuOtywxSNIVAxo7v//vvTegkIADf8zq9Tp06J3ttvvy16bsgjgum5554TvcaNG4tey5YtrdpLCKMxxjRq1Mjfwhza+2s3YNEYQvOCQAsqu+WWW6z6iy++EDPffvut6Lm/v7Zt24qZ/Pnzi54b8qh9ZmzRooXorV27VvSQuq666irR+/rrr636xx9/FDPLli0TvXXr1ll1p06dxEyvXr1ETwvkdGnB1+7zz6pVq0JeB/4tXrxY9LR9FjcQOE+ePFFbU7S5wZ7z5s1Lo5WEh29YAwAAAAAAAAACgQ1rAAAAAAAAAEAgsGENAAAAAAAAAAgENqwBAAAAAAAAAIGQ4UMX7777bqtesWKFp587ceLEJa9jjB7qgWDKlSuXVZ85c0bM5MiRQ/SKFi1q1T169PD1+Nq9Mn36dNG75557fF0fMEaGybj3PYLhmWee8dQbPXq0VUcymPGbb74RPTdceNy4cb6uDSBz+uOPP6y6bt26Yub48eOptRykgtatW4ueG3CYL18+MeM1iNGL8+fPW/WkSZPEzEMPPRSxx0Pk/P7776I3fPhwqx45cqSY6dChg6eeF+7jPf/882KGEPNgqlatmui5IZoVKlQQMz179hQ9N7je7/vrXbt2id6oUaNEj5DF1LV8+XJPvYEDB1r1kCFDxEyxYsVEz72ntNe9SHL3NN3PcMYY889//jOqa0gtfMMaAAAAAAAAABAIbFgDAAAAAAAAAAKBDWsAAAAAAAAAQCBk+DOs+/btm9ZLQAC45/p8/vnnYqZLly6ilzWrvz+Rjz/+2KpffvllMbNhwwZf1wYu5sCBA1adkJCQRitBJLjnWmvnXANApDVs2DCtl4B0rGDBglY9YsQIMdOiRQvRa9y4sVUfPHhQzHz66aeiN2HCBKvetm2bl2UiAJKTk0XPze9wa+BvX375pehy796TAAAgAElEQVT16tXLqrX7p2zZsqLnnlmdmJgoZmbPni167hnr+/fvFzNJSUmih/RhzJgxnuYeffTRKK8k8+Ib1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACgQ1rAAAAAAAAAEAgZPjQRcAYGdxy7733ihmtB6Qn586ds+pKlSql0UoAAAD00EWtBwDhmjJlyiVrAOkL37AGAAAAAAAAAAQCG9YAAAAAAAAAgEBgwxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBwIY1AAAAAAAAACAQ2LAGAAAAAAAAAAQCG9YAAAAAAAAAgEBgwxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAIGT1+4MpKSmRXAcyGe4fhIP7B+Hg/oFf3DsIB/cPwsH9g3Bw/yAc3D8IB/cP/OIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCDGX8/X8mJiYA8aYndFbDtK5cikpKXEX+4fcP7gE7h2Eg/sH4eD+QTi4fxAO7h+Eg/sH4eD+QTi4fxCOS94/f7usDWsAAAAAAAAAAKKFI0EAAAAAAAAAAIHAhjUAAAAAAAAAIBDYsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACIevlDBcpUiQlPj4+SktBerd69eqDKSkpcRf759w/uBjuHYSD+wfh4P5BOLh/EA7uH4SD+wfh4P5BOLh/EI5Q98/fLmvDOj4+3qxatcr/qpChxcTE7LzUP+f+wcVw7yAc3D8IB/cPwsH9g3Bw/yAc3D8IB/cPwsH9g3CEun/+dlkb1s4D+P1RZCApKSm+fo77B8Zw/yA83D8Ih5/7h3sHxvDcg/Bw/yAc3D8IB/cPwsH9g3D4uX84wxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBwIY1AAAAAAAAACAQ2LAGAAAAAAAAAAQCG9YAAAAAAAAAgEDImtYLAAAAQNorWrSoVe/fv1/MjB07VvRSUlKs+pdffhEzs2bNCnN1AAAAADILvmENAAAAAAAAAAgENqwBAAAAAAAAAIHAhjUAAAAAAAAAIBDYsAYAAAAAAAAABAKhiwAAwJMsWbKI3vnz59NgJYiGNm3aWHX//v3FTOnSpUUvLi4uamsCAAAAkPnwDWsAAAAAAAAAQCCwYQ0AAAAAAAAACAQ2rAEAAAAAAAAAgcCGNQAAAAAAAAAgEAhdTEdiY2OtumPHjmJmypQpIX8uKSkpsgsDAKSJmJgY0UtJSbFq7bVi8ODBonf11Vdb9eHDh8XM0aNHRe+HH36w6lmzZomZOXPmiB4io1SpUqK3e/fukD/3/PPPi17jxo2tukaNGmJmz549ojd69GirfuaZZ0I+PjK+Jk2aiF7fvn1Fr2HDhlb90UcfiZkRI0ZEbF2AMcZky5bNqrNnzy5mTp48mVrLQSro16+f6Llhw7/88ouYWb16tejNnj07cgtD1LivL8YYU7hwYdFLTk626oSEBDGzdetW0btw4UIYq0NaGj58uOgVLVrUquvWrStmjh07Jnrbt2+36mLFiokZ7fOR+37n3Llz+mIzMb5hDQAAAAAAAAAIBDasAQAAAAAAAACBwIY1AAAAAAAAACAQOMPaI+0cPu3Mop9//jlqa7jvvvus+oUXXhAz7733nugtX77cqocOHSpmlixZEubqEAm5c+e26ttuu03MaGfsnT592qrdc7iMMWbHjh1WvX79ejHjnn2LYGjQoIHoFSlSRPQKFSpk1drv8+OPP47cwpDmbrrpJtHr1KmTVXft2lXM5MqVK+S13XPcjJGZCMbI8wHLlSsnZgoWLGjVU6dODfn48MbLedWaLFmyiJ575l7WrPJt4qFDh0Rv2LBhvtaAjO2JJ54QvXbt2oX8uUcffVT0Dh48KHoTJkzwtzCkWyVLlhS9W2+91aq13Abt59zPPlu2bBEz7tnF2vnGCIb4+Hirfuqpp8RM+/btQ16natWqoqftA7jnIGtn7/O5KnJ69+4teh06dBC95s2bW7X72TrS3P2gnTt3ipkPP/xQ9EaOHBm1NUHS9u20HKBoatmypeg9++yzVr1p0yYxo2V/bNu2LXILCzi+YQ0AAAAAAAAACAQ2rAEAAAAAAAAAgcCGNQAAAAAAAAAgENiwBgAAAAAAAAAEAqGLxpguXbqInnuwvxtqZ4wesBjN0EU3KCJnzpyefi5//vxWrf27IPV1795d9NxwNC0gQAvw+OOPP6xaC2Y8evSoVVeuXFnMfPvtt6LnhrNpj+8GXBhjzJdffil6kLTQvPr161t1vnz5xExiYqLoVahQwarLli0rZrR7igC89EEL2rzxxhtF795777VqLTjPL+1apUqVsuqKFSuKGTfESAv8mzJlSpirw+U4e/as6LnBZGfOnBEze/fuFT2CpWCMMeXLl7fqatWq+bqO9j5Vez1DxpItWzardgOEjdHDp9ywPfdzz8W472+198VuKJcWzHjkyBFPj4foatOmjVXXqlVLzPz111+il5SUZNV//vmnmNHe+9x2221W7QafG2PM+PHj9cUipAEDBlj1fffdJ2aqVKkieu77y3PnzomZU6dOiV5ycrJVu5+bjTEmLi5O9PLmzWvV7uugMcYMHTpU9Nz3yj169BAz8GbIkCEhe34DFhMSEkTPvVeMMaZEiRJWrX3e1gLv3Z9za2OMmT17tug1atTIqrV9gYyCb1gDAAAAAAAAAAKBDWsAAAAAAAAAQCCwYQ0AAAAAAAAACIRMd4Z10aJFRW/SpEmilzt3bqvetGmTmPnwww8jtzDH9ddfL3pXXXWVr2uVKVPGqlesWOHrOvDOPT/rySefFDPuObPGGFO6dGmrPnHihJjZunWr6DVr1uySj2+MMQcOHLDqK66Q/7+qVatWolegQAGr1s7v084L5AxrqWHDhqLXsWNH0XPPuFq/fr2Y0c7YW7t2rVXfeeedYubKK68MuU4E09133y16999/v+j5PbPaPdPPPdfRGP2sTvf5xssZf+3atRMzq1atEr3Nmzfri0XYbrnlFtGLjY21au0sYfd5Bvibm82RJ08eTz+3a9cuq9byYLZt2yZ6xYoVs+p9+/Z5ejykPe0sz549e1r1HXfcIWbKlSsnem6+Q44cOcSMdma/+95Ve3903XXXWfWyZcvEzOrVq0UP0aWdb96iRQur1l6rli5dKnobN260au0+0HKj3Awa7T23e5+PGTNGzEB3+PBhq9beIx47dkz0Fi5caNXaa8fKlStFz70PNA0aNBC9m2++2ar/+c9/ihnt+a5evXohHw/eaK8L7l6eRjvTfuTIkVb91VdfiRkt38XdY3TPuDdGf25p3bq1VWv/LtrPTZ482arbtm0rZjIKvmENAAAAAAAAAAgENqwBAAAAAAAAAIHAhjUAAAAAAAAAIBDYsAYAAAAAAAAABEKmC138/vvvRU87lP38+fNW/dFHH4mZL774InILc/Tu3Vv0ChUq5OtaWuAfIqdChQqi161bN6vWwkGqVq0a8tpuUKIxxpw8eVL03Pu1RIkSYqZgwYJW7YaTGKPfY26gWnJyspjp2rWr6LlhS1qAZGZTrVo1T7333nvPqrXnHy8qVaokem6okTHGjBgxwtf1EV3ua5MWiuqG5GnOnTsnenv37hW93bt3W/X+/fvFjBbWWqpUKauOj48PuaYbbrhB9LRQSe7NyNCCeLVAPDdIRgutcgONgL+54apubYwxKSkpoucGH2nBdlroWfHixa1ae65LTEzUF4tUo7237Ny5s+i5QbBuaLwx+j2lPb+5smXLJnru66d2nauvvtqq3RBGYwhdjLaKFSuK3j333CN6MTExVv3uu++KmXXr1oV8vA0bNoie+xnKGPlZ6/rrrxczffr0sWot9HHRokUh15QZzZo1y6pnzpwpZtxgxmhbvnx5yF6zZs3ETOPGjUXPDYfV7rHU/vdLL9y9Fze42Rj5e9H+zgYPHhyxNR09etSqx48f7+nnnn76aasePXq0p5+rXbu2t4VlAHzDGgAAAAAAAAAQCGxYAwAAAAAAAAACgQ1rAAAAAAAAAEAgsGENAAAAAAAAAAiEDB+66AY21ahRQ8xcuHBB9H7++WernjJlSkTX5XLDtLTwqezZs1u1FkrjBhQYIwPc4J8WkHfnnXeKnhvKUrZsWTGjBZr98MMPVj158mQxM3/+fNHr27evVT/88MNiplixYlZ97NgxMaOF0pw+fdqqtYBFDSGLkvZco4VAfPbZZxF5PC10MVeuXBG5NqLv2WeftWrt9csNXDVGBn9oEhISRG/JkiVWvWrVKjGjPW+4r1eVK1cWM27Anxs2Y4wxDRo0ED33tXHu3LliBqH16NFD9LTfgXs/rVmzRsxo4dWAMcZUr17dqrVgz7Nnz4qeG3K2ZcsWMaMFvrphfgcPHvS0TkRXu3btrLphw4ZiRgser1WrllW7gb7GyGA9Y+TzVlJSkpjZuHGj6Lnvzd33ydrjaUGQWqjkoUOHRA/+DBkyRPSqVKkietOnT7dqN8w1HFr43eeff27VXbp0ETM1a9a06scee0zMELqoSy9/Q25InhYqrn3+O3LkiFUTsOjdn3/+adUHDhwQMx9//LFVT5gwIapr8ssNldU+Z2mBwHv27InamoKGb1gDAAAAAAAAAAKBDWsAAAAAAAAAQCCwYQ0AAAAAAAAACIQMdYb1bbfdJnqPP/64VWtn4O3cuVP03n77bauO5DkxxYsXF70XXnjBqrWzhN1z/1577TUxM23atDBXh//lnhmknVddvnx50YuNjbXqb775RsysXr1a9NwzrH/99VcxU7RoUdG7+eabrTp//vxixj0DvWDBgmLmjz/+EL1GjRpZtXZ+n3YeLqTSpUuLnnamdHJyckQez/2dG6OfnY60p50x3KxZM6vWzhzWzoN1aecATp06VfTefPPNkNfSxMfHW7X7/GeMMS1atLDquLg4MeOeX2oMZ1ZHyrXXXit62rmr7nukTz75JGprQsZz1VVXhZxJTEwUvTNnzli19tqlnUe7Y8eOkNeOFG1NkXqtTs+6desmevfdd59Va2eZa+8l3feuWbPKj6naa97x48etWjvP9OTJk6LnftZKSUkRM+7zZIcOHcSM9l49UlkkmZF7drmWi6Fldezdu9eq9+3bF9mFOXbt2mXV7mc4Y4wpWbKkVV999dViRssn0c5cR9qbNGmS6LnPdxrtfOpXXnklImuCMe+8805aL8ETd1/SGLnPor1eavuQM2bMiNzCAo5vWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBkG5DF7UAoUGDBole3rx5rfr8+fNi5pdffhE9LZDKjxIlSojelClTRK969epWrQW5LF682KoJWIw+N/jj1ltvFTNVq1YVvXXr1ln1G2+8IWY2b94sem4oQ4MGDcRMy5YtRe+aa66xai3QzLV06VLRW7RoUcifI2DRu4oVK1p1vXr1xMxvv/0Wscdzw+9atWolZrRATqQu7XfQr18/0XNfP7TA1ZiYGNFzQxa//PJLMbNgwYKQ6/Rq5syZVq2Fk1WoUMGqtdBX7s3ocZ+LjNGDXdywspUrV0ZtTUjf3PcdxsjXIE1SUpLo5c6dO+S1582bJ3qRClnMmTOn6J0+fdqqCVg0pmzZsqLXtm1b0XN/f9prgvbf3HXs2DHRc8PQtTntd6U93pYtW6xau3/dzwHlypUTM7Vr1xY9Qhf9a968uVXny5dPzGif57Xg+NT0wQcfiF6TJk2sWgtfd0PXjCF0MS3cddddVj106FAxo4VmurRg2LFjx4oeodYZy6hRo6xaC+jVXmO09+KuU6dOiZ62f5lR8Q1rAAAAAAAAAEAgsGENAAAAAAAAAAgENqwBAAAAAAAAAIHAhjUAAAAAAAAAIBDSTejisGHDrPqBBx4QMyVLlhQ995Dy//73v2Lm559/Fr26detatRZYtGPHDtErXry4Vd9+++1ipk6dOqJ3xRX2/zvYuXOnmLn33ntFD9Hlhnu64SvG6EFhbijMmTNnPD3eHXfcYdXt2rUTM02bNhW9MmXKWLUbumaMMRMnTrRqLbSTkI/Icv/+GzduLGbWr18fscdznxdz5MghZvbv3x+xx4M37muTFryrBXJ6oQVxLFmyxKq/++47MaOFvvrlPi/u3btXzLjhfW4IozF6WOzatWutWgu3QmjXXXedp7msWe23hVp4kBf33HOP6HXu3Fn03GDR3bt3ixk31HPy5Mm+1oTI0gKgvXDDpY0xplChQlatvYfRwhojxQ1YhO4f//iH6Gnh4O57D+09sNY7d+6cVR89elTMuJ+XjDHm4MGDl6yNkYGyxsj377ly5RIz7uu3tm4tFBD+ua9X2uf75cuXi14k30/7sWnTJtGbP3++Vbdv317MeAkghX9aUKIWqOju/XgJETbGmK1bt1r1wIEDxYz7PgbpR82aNUXvnXfeET13z8h9Px0O7TnQ/Ry1aNGiiD1e0PANawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBEIgzrLNnz27Vo0aNEjOPP/64VWfLls3XY7lnSRljTIcOHURvxIgRVn38+HExo51V5Z5jpp1XrZ2R5l5fO1+0SJEiVs1ZtNFXrFgxqz5//rynn3PPyxoyZIiYOXnypOi594t2Ll6ePHlEzz237c033xQz7777rr5YRM2VV15p1dqZ0kWLFhW9a665xqpjY2PFjHamVqNGjazaPRfUGP2+Q3S5Z4u7v6dwaGc2umdY//TTTxF7PM2vv/5q1e7zpjHybEDt3tSkpKRYtfu3oT0+5H9vr2es7tu3z9fjue/R3PdQxhiTN2/ekNfR7h33nOvKlSuLGe01FtHl90xpN+PDGJn/op1lTsZG6nP/1sqVKydmtL9Z971OcnKymDly5Ijouc8/u3btEjNbtmwRPffzmPY+Weu5evfuLXoxMTFWrX2G27NnT8hrQ9eqVSvRa9GihVW7+wTGGLNixQrR8/v6FU1XXXWVVWv3oZvVAe+0e6Nbt25W3b9/fzFTvXp10XPP0Nc+82sZY24uDedVp2/uPdWrVy8x4+7JGSM/r2j3j/b+x4vcuXOLXps2baz6/fff93Xt9IBvWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBEIjQxSpVqlh1165dxYzfkMUDBw5YddOmTcVMtWrVQl5HC4jyEjTihnVcTNas9q9CCydq2LChVZ85c0bMNG7cWPTcwC0tKGfv3r2e1pnZuAfvJyYmihkthMYN1bj11lvFjPZ7cA/s14JcZsyYIXpaoATSXqVKlaxaC01o27at6LmhfCdOnBAzWs99TtLCXbSAEkSOFgjXrFkzq86ZM6eva7uvZ8YYM3fuXNH7/vvvrVoLvIomLYzEDaryGjxSsmRJqyY01Bst8NkLN+xOC8SqX7++6PXp08eqtfdsR48eFT33eUx7L1KwYEGrbteunZjRXptfeeUV0UPkaH+LbmiV+97WGBmwaIwxZ8+etWotoBipr0yZMlatfYbSflfua866devEjBaQ535e0Wa+/fZb0XOf77TPXtrzjxvoeujQITHj0u5pLz8HXfPmzUXPve+2b98uZlatWhW1NUWSGy6qfb7X/v3gjRZA736ucj+LGSNfc4wx5tSpU1ZdoEABMaO9v23durVVa4GybjCjMfpzEtLehQsXrHrevHliRgu8r1ChglVr75+1PSN3b0ALa8yfP7/otW/f3qrfeustMdO3b1/RS4/4hjUAAAAAAAAAIBDYsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCIEIXXRDwGJjY31dRwuAcUPs3EAYY4w5duyY6LmH8e/fv1/MaIFmbkCU12ApN5BKC6Fwgz60cC8tDMT9b/DFF194WhPk789vYJ0Wsnb69GnR27Fjh1VrQS4ELKYf7vPG4cOHxYwWwOLeZ+59YYweMuaGH7mBtsbIUAhjZMjj0qVLxQx0brhTkyZNxEy9evUi8lgzZ84UPe139dtvv0Xk8a655hrR04JqKlasaNXaf4NSpUr5WoMbNOm+VkJ34403+vo5N9jlpptu8nXtFStWiJ4WLPXDDz9Y9bXXXitmOnToYNVaoKMbzIjocwPrvNJCq9ywTS2YEanPDVjVgp807meRzZs3i5n58+eLnhtupYWXaRYuXBhyRnv/7gZg+X2d0j57wZs2bdqInhtCtmvXLjETqfc5kaQFqrmhgNq/i/a+Ct6sXr1a9Nz3ylqIXenSpUXP3Q/S/q6116/4+HirfuaZZ8TMI488Inrjx4+36okTJ4oZpD53r/C7777zdZ1hw4aJXvny5UM+XseOHcXMqFGjRM/dP3jwwQfFzO+//27Vr776qr7YgOMb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACIRCHbp06dcqqtXN7CxUqZNU///yzmHnttddE75dffrHq//73v2LGPSvLGGNOnDhh1Z06dRIzAwYMEL0yZcqInks7Rzs5OdmqtfO43bPVtP9Ob7zxhuhpZyXDG/fMu8KFC/u6zoULF0RP+x27fwtz58719XgIhm7duln1HXfcIWa8nDuvnW+3aNEi0WvXrp1Va2cna+epa89J8Mb9XVWtWjXkjFf/+c9/rPrZZ58VMzt37vR17SuukP+/2r1ftNcOL2cB1q1bV8zExcVd7hKNMcasX7/eqtetW+frOpnNvn37rFp7DklISBA997w79znFGP11cMGCBVatnbe3YcMGfbH/wz2H3xhj7rzzTqvWzjH3e54y/NNyXNzcBjfXxRj9Nc/NAvByJjGiz/1daa8bGvc9xfHjx8WM9v726NGjl7G6y+N+zjJGnj3r1po9e/aI3tSpU/0uK1Nxz3M2Rn/Od1+vVq1aJWa0s6DTWu/evUXP/XdesmSJmPnrr7+itqbMaPLkyZesw6G9fg0ZMsSqn3zySTFTtGhR0Xvuuees+vvvvxcz7ucApG9alotLO2e6QYMGote5c+eQ13L3KrXsiI0bN4a8TlrjG9YAAAAAAAAAgEBgwxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAIAQidHHTpk1W/dZbb4mZP//806pnzJgRscd3AxY1bdu2FT0t0MyLLVu2iJ4bBnnkyBExM3v2bKv+5JNPfD0+dFq4VJ06daxaC6fT7h832E4LL8uSJYvolS9f3qq1MBKkX1999VVUrz9lyhSr1oIbtNBFLRQU/mjBKm7AqhswZoweovvpp59atd+ARS0crXr16iF/Tnv+0V73brnlFquuXbv2Zazu/9ECt7TwPoTmhlS54ZXG6PeFG5iphZBpYVerV6+2ai8Bixo3YNEYGdy3bNkyMXPw4EFfjwdvihQpInrXXHON6Lnhqtprixbc54aO+b1/EFluiKZbG2NM/vz5Rc9LcHQ0AxY1sbGxote9e3er1u5N99/l448/FjM8/3jjBtkbo3+u2rZtm1WvXLkyamsKR40aNay6U6dOYsZ9Dpw3b15U14To0u5X933qsWPHxMz48eNFr1ChQlat7X+1bNnycpeIDKhLly6i16ZNG6vOmzevmHFD0k+dOhXZhaUSvmENAAAAAAAAAAgENqwBAAAAAAAAAIHAhjUAAAAAAAAAIBDYsAYAAAAAAAAABEIgQhddWqCFFvQRTcOHD7fqrl27+rqOFrA4ZswY0Vu8eLFVawFcWnAf/HGDpYwxZvTo0aJXokQJq3YDEozRD7BfunSpVWthe8OGDRO9bNmyWbUWzAhcTMGCBa1aC/fTaM9T8McNTjXGmL1791q1+3syxpjExETR27NnT0TWpAVxaAFUFStWtOpWrVqJmapVq4peuXLlrNrL85b2ejZhwgTR497058svv7TqkSNHihntXvVCC24tVaqUVVeuXFnMaCGPgwYNsmrtnjtz5oxVb968WczMmTNHXywi4rXXXhO9+vXri54bSKX9/dasWVP0kpOTw1gdosUNOPX6nsJ9L6u95qW2sWPHip6XgGA3ZNZ9zoJ32nsaN5Rac/jw4Wgs57JogZwvvPCCVWuh224g8bhx4yK7MASOFnj/4IMPil6VKlWs2n0PjugrW7as6GnB4kG0Zs0aq27atKmYcV+zixUrJmb+85//RHZhUcA3rAEAAAAAAAAAgcCGNQAAAAAAAAAgENiwBgAAAAAAAAAEQiDPsE7t86rvvPNO0Rs6dKiva82dO9eqJ0+eLGamTZvm69rwr0yZMlY9YsQIMaOdHeWe87pkyRIxM2nSJNFzf8fatbt06SJ67jmgJUuWFDPAxbhnceXIkUPMaGcBnjhxImprymzOnj0rerly5bJq7Uxp7efuuOMOq9bOy9fOiC1cuLBV33DDDWJGO0/YPdO4QYMGYkY7x9+L8+fPW/XMmTPFzNSpU31dG1JSUpJVa2fUaWfZaWebu+Li4kSvSZMmVl2rVi0xo50VGB8fH/Lx3PeE27ZtEzOLFi0KeR14554zXaBAATGjPfe454tXr17d0+O5z5HuGcjG6M+RiK4NGzZY9Y4dO8RMpUqVRM897zd79uwRXVco2pn9vXr1Cvlz7uuUMcYMGTIkImuC/l7TzSgwRr6H0d7LRpP2WqWd4+++t9q4caOY+ec//xm5hSHdOnLkSMgZLSOgWrVqVr1p06aIrSmjc99HvPzyy2KmQoUKordy5Uqr1vZ53GyitJA7d+6QMxkl/45vWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIBDasAQAAAAAAAACBkOqhizlz5hS91DwQvOBW4eMAAAm0SURBVEePHqI3duxY0dMCX1x//fWX6H399ddWTcBiMNx3331W3bp1azGTJUsW0du3b59VP/XUU2JmxYoVIR//yiuvFL1y5cqJnhssQhgeLocbVHPs2DExc/LkydRaTqb066+/it7VV19t1VqAmRZ254auuuErxughxW6YnvbcpoUYZc1qvyXwEuihOX78uOhNmTLFqp955hkxc+jQIV+Ph9B69uwpep999pnoueGJGi1Azb3HU1JSxMyFCxdE78CBA1a9bt06MTN//nyr1oJzEFnt27e36qNHj4oZN1jPGGPq1atn1VWqVPH0eG7YHQGLweSGuXrlBvoao4f6rl+/3qq1oGH3ucYYYx577DGrbtmypZjRwo5dY8aMET3CgCNHe+3QgsBdiYmJ0VjO/3XzzTdb9dChQ8WMFiTsrl3bT1i6dGmYq0N6c9ddd4medv+4du/eLXp//PFHRNaUGZUoUcKqb7vtNjGjhS66c0888YSYmThxougNGjTocpeo0l4vtfBfNxzbizp16oheeniO4hvWAAAAAAAAAIBAYMMaAAAAAAAAABAIbFgDAAAAAAAAAAIh1c+wTs3zqo0xpk2bNlb95JNPihn3nFCvPvzwQ9F78803fV0L0eWeo6id6arJly+fVbtnBF/sWu55evfcc4+Y0c7Tc88HdM+UBS7l3LlzVh0TEyNmcuXKlVrLyZSWLVsmeu7Z07feequY0Z5H3PMe3bOpjTGmSJEiIX9Oo50n7J5Je+bMGTGjnU3snjv6xRdfiJmXXnop5JoQPXv27BG9G264QfTefvttq+7evbuY0c4X3rRpk1Vv375dzBw5ckT03JwI7Qz4WbNmiR4iRzs/v0uXLlatnS9funRp0StVqpRVe30P8+mnn3qaQ9py/86NkeeWGyPPDq1UqZKYGT58uOi574G1PKGqVauKnpYT48Unn3xi1cOGDfN1HXijnUWtZa1Ur17dqt3P8sbIz2ca7b6rW7eu6DVt2tSqtfdau3btEr1Ro0ZZNa9VwVC8eHGr1jLHIsk9u3j06NFiRnsuc/Nepk+fLmaSk5PDXF3m5f7NbtiwQcyULVtW9Nz3LdpzTd++fUXvpptusmotg0H7XOVyXz+NMaZMmTKip2URudx8kG3btoX8mSDiG9YAAAAAAAAAgEBgwxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAIGT4RLeKFStatRak4MWBAwdE76uvvvJ1LaS+o0ePWrUWHKYF1LkH2g8ePFjM9OnTR/Tcg/fz5MnjaZ0HDx606qVLl3r6OcAYY66//nqrjouLEzPacxkiZ8aMGaLnBipqYWVuSIwxMphVC8zUwhpPnjxp1Tlz5hQzbtiLMfL5RwvYmjdvnui5QUMJCQliBunDAw88cMn6YtzATi3UE8HkBt0ZY8zy5cut+pZbbhEz2nOWF1OnThW9cePG+boWUpcbUmiMMfXr1xc99zVH++yl/Vzu3LkveZ1wfP3116LXrVu3iF0f/mjvM2rVqmXVPXv2FDOtWrUKee2rrrpK9NwQMmNkkPCiRYvEjBakp4W4IXrc9xnGGNOvXz/Ra9u2rVVr73d/++030StUqJBVlytXTszUrFlT9LRAPC9efvllq3799dd9XQfe3HvvvaL38MMPi979999v1VpgphbEWLlyZat2P8NFmxau/o9//MOq586dm1rLiSi+YQ0AAAAAAAAACAQ2rAEAAAAAAAAAgcCGNQAAAAAAAAAgENiwBgAAAAAAAAAEQoYPXdy9e7dVaweSJycni97+/fut+t///reYWbFiRZirQ2pxD5nXQhOuu+460XPDy7Sfi42NFT03rFGj3T9PP/20VX///fchrwP8zQ3lO3XqlJg5fPiw6DVq1MiqCfv0Twubmz59ulXv27dPzLRv31703NDg2rVrixkthMYNUdN+51o45OTJk0OuE9AQspixuKGLd9xxh5jRwsvc3nfffSdmunbtGubqkFY2btwoelqIkxs0rgUsFi5cOGLrcsOkP/vsMzHTv3//iD0eIue5554TPff32aVLFzGjBd2dPn3aqrdt2yZmFi5cKHrLli2z6m+++UZfLNJUgQIFRK98+fKiV7VqVasuVaqUmNHec0fKiRMnRG/AgAGi9+6770ZtDZCOHDkietrzj9vr3LmzmLnmmmtE7+6777ZqbS8oe/bsoue+b8qSJYuY0T7Pf/7551bdq1cvMZOSkiJ66RHfsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCBn+DGv3/Dz3XCxjjPnrr79Eb+vWrVFbE1LfV199ZdWVKlUSM7///rvoValSxaq1s4BiYmJEz71/FixYIGY++ugjfbGAT1u2bLHqJUuWiBnt3NGEhISorQmS9nvRegCQFiZOnGjVa9asETOlS5cWPfe9FjK+F198UfS2b99u1T169BAzZcuWFT03E2bXrl1iZuXKlaK3ePFiq54zZ46+WKQL77zzziVrZE6HDh0SvbVr14peXFycVderV0/MaOdh58uXz6q1vaB58+aJ3qJFi6yaM9AzFi0TQesNGjTI1/Xd/Cntfj169KjorVu3ztfjpUd8wxoAAAAAAAAAEAhsWAMAAAAAAAAAAoENawAAAAAAAABAILBhDQAAAAAAAAAIhAwfunjy5EmrJtgKxughMZHkBsckJSVF9fEAY4z5+OOPL1kbY0yWLFlETwtiBADAGGNWrFjhqQcYY8z06dMvWQNAJHz44YeeekBQnTp1yqrZq5T4hjUAAAAAAAAAIBDYsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCBk+dBFIC4QsIqgIWAQAAAAAAEHGN6wBAAAAAAAAAIHAhjUAAAAAAAAAIBDYsAYAAAAAAAAABAIb1gAAAAAAAACAQGDDGgAAAAAAAAAQCGxYAwAAAAAAAAACgQ1rAAAAAAAAAEAgsGENAAAAAAAAAAiErH5/MCUlJZLrQCbD/YNwcP8gHNw/8It7B+Hg/kE4uH8QDu4fhIP7B+Hg/oFffMMaAAAAAAAAABAIbFgDAAAAAAAAAAIh5nK+nh8TE3PAGLMzestBOlcuJSUl7mL/kPsHl8C9g3Bw/yAc3D8IB/cPwsH9g3Bw/yAc3D8IB/cPwnHJ++dvl7VhDQAAAAAAAABAtHAkCAAAAAAAAAAgENiwBgAAAAAAAAAEAhvWAPD/tWPHAgAAAACD/K2nsaMwAgAAAGBBWAMAAAAAsCCsAQAAAABYENYAAAAAACwIawAAAAAAFoQ1AAAAAAALwhoAAAAAgIUAJ8jfKLP7geUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b68e630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Check output of autoencoder: print one input/output pairs\n",
    "in_imgs = []\n",
    "out_imgs = []\n",
    "numprint = 10\n",
    "for i,(img,label) in enumerate(test_set):\n",
    "    if i >= numprint:\n",
    "        break\n",
    "    out_img = decoder20(encoder20(Variable(img)))\n",
    "    in_imgs += [img[0][0].numpy()] # img is 1x1x28x28 tensor\n",
    "    out_img.data.clamp_(0, 1)\n",
    "    out_imgs += [out_img.data.numpy()]\n",
    "\n",
    "fig = plot_images(in_imgs,out_imgs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Denoising with trained decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to estimate the latent representation\n",
    "def invert_prior(net, noisy_image,latent_param,learning_rate=0.1, num_steps=500):\n",
    "    optimizer = torch.optim.SGD([latent_param], lr=learning_rate) \n",
    "    iter_idx = 0\n",
    "    while iter_idx <= num_steps:\n",
    "        # correct the values of updated input signal\n",
    "        optimizer.zero_grad()\n",
    "        output = net(latent_param)\n",
    "        loss = F.mse_loss(output, noisy_image, size_average=False)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        iter_idx += 1    \n",
    "    return latent_param\n",
    "\n",
    "# denoise by recovering estimating a latent representation and passing that through the decoder\n",
    "def denoise(net,noisy_image,K,learning_rate=0.1):\n",
    "    rec_rep = Variable(torch.randn(K), requires_grad=True)\n",
    "    rec_rep = invert_prior(net, noisy_image, rec_rep,learning_rate)\n",
    "    if(sum(rec_rep.data.numpy()*rec_rep.data.numpy()) > 10000):\n",
    "        print(\"stepsize to large?\")\n",
    "    return net(rec_rep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize denoising performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABawAAAEsCAYAAAAvofT2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXnclnPa/494MPadwZDdWJItZYsUlVIqslRklyWMLesYYdAI2SKiqOxEZEtaLMkypLIvyeSHMRkzYxv6/fW8Xs/3c7y5rzlv88w183ze/x2H477v6zrP73aeeh3vJgsWLAhjjDHGGGOMMcYYY4wx5l/NQv/qD2CMMcYYY4wxxhhjjDHGRPiFtTHGGGOMMcYYY4wxxpg6wS+sjTHGGGOMMcYYY4wxxtQFfmFtjDHGGGOMMcYYY4wxpi7wC2tjjDHGGGOMMcYYY4wxdYFfWBtjjDHGGGOMMcYYY4ypC/zC2hhjjDHGGGOMMcYYY0xd4BfWxhhjjDHGGGOMMcYYY+oCv7A2xhhjjDHGGGOMMcYYUxf81z9SvNJKKy1Ye+21/0kfxfy788ILL3y6YMGClX/ov3v8mB/CY8c0Bo8f0xg8fkxj8PgxjcHjxzQGjx/TGDx+TGPw+DGNoaHx89/8Qy+s11577Xj++eerfyrzH02TJk3e/7H/7vFjfgiPHdMYPH5MY/D4MY3B48c0Bo8f0xg8fkxj8PgxjcHjxzSGhsbPf+OWIMYYY4wxxhhjjDHGGGPqgn/oX1j/T5o0afJTfg7zb8qCBQsq/dzxxx+fcp9//nkR/9d/5eH597//PeVWWGGFIl588cVTzfvvl/8D56uvvko1a621VsqtvvrqRfzmm2+mmtmzZ6fcwgsvXMQrrbRSqll66aUb/HvvvvtuqnnllVdSbv3112/wd6+66qpF/Mknn6Sab7/9NuVatWpVxOPHj081r7/+esrtscceRTx//vxUc8stt6RcLQwYMCDl3nnnnSL+7rvvUs3XX3+dcssss0wRL7bYYqlmlVVWKeK//e1vqWbMmDEpt8EGGxTxs88+m2q6du2acksttdSPfsaIiCWWWCLlllxyySL+/e9/n2pefvnlBj8Dre86pnW+/tDPff/990VcyxyOyPeTxs+kSZNSrhZOPPHElNP58NZbb6Wa7t27p9w333xTxMsuu2yqWWih8v8NL7/88qnmxRdfTDmda/vss0+q0XseEfHSSy8V8dtvv51qevfunXL6ryA222yzVDNjxoyU03WD7rHuFTQ/6bvomljL9Y2I2GijjYr4ySefTDVDhgxJuYa4/fbbU+4Pf/hDEevYjYho06ZNyuk9p++ve8c666yTam688caU+/DDD4t4tdVWSzVNmzZNuT/96U9FTGNVf3dEvi9dunRJNW+88UbK6d5Me4n+PRq7tJ+tt956RTxhwoQGa+gz6RiMiLjkkktSrhb23HPPlNPx2759+1Tz2muvpZyuyTNnzkw1m2++eRHTvJszZ07K6XX59NNPU826666bcvPmzSvitm3bpho6M+k6+tlnn6WaFi1apJzuqbQ+rLHGGimn6HWKiLj00ksb/PsfffRRynXu3LnBmr59+zb4mYhTTz015f7f//t/RUznTT2L0M+tueaaqUbXqD59+qSaKVOmpJyu9zSv6XP++c9/LmI6s/3yl79MOR3XDz30UKqhM5OelbfaaqtUQ+udouM3IuKPf/xjEdPn/tnPfpZyu+22WxHfdtttqebyyy9v8DMRd911V8rp+YvOVXTe089O1/zYY48t4hdeeCHV0Jqk94XWEVqT9Mx7xRVXpJr99tsv5fSMTXvxBx980ODnXGSRRVKNntt03kVETJ48OeXatWtXxPSMSmeYX/ziF0VM61///v1TrhbOPPPMlNNxTddg7ty5Kaf7Me17I0aMKOJBgwalGjpX6DMana+33XbblNPnZJovdB/0XPbqq6+mGpr/ut7R2fnBBx8sYh0XERH33ntvyu2yyy4N/u7dd9895fSsoe8XIiIOPPDAlKuFyy67LOV0f1x00UVTzdZbb51yY8eOLWI6Y+tzI5156Ryj7yZ+/vOfpxo9J0ZEvPfee0VM76N0fY/Ia8tyyy2Xaui9jv4cnbunTp1axHfccUeqobmgzwI0Dugsp+9+6PmEzuIN4X9hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXVBZumhMY/j4449TbsUVVyxikntR4/taZF5/+ctfipgaxZO4S+VZW2yxRaohgYY22id5IokFvvzyywY/0yGHHJJyjz76aBGTZEPFMSrD+6HPpM3xSX5AUiq9xyRIqcrTTz+dcs2bNy9iunYkLdDPqdKEiCzpJKkICTVUQHXVVVelmquvvjrlVCB29NFHp5ovvvgi5VQUSPNsp512avDvkeho+vTpRUzi1McffzzlVLZEwokdd9wx5VQOufLKK6eaquhaE5HnNokK6Zrr3D7iiCNSjco5SNZBuWbNmhUxiWhV+BeRpTAkvLnppptSTq8Licg23njjlNM1l6S2KmgiAQzJnlSiNmrUqFRDY0rvZ1VBsHLnnXemnK6JNHZI1KPXgO6T7lU6DyN4rVNB6MiRI2v6OZW2/Pa3v001HTp0SDkVd5GElsSIOsZoL9l1112LmK4vfRdd/1SeG8HSRd2HScRblZYtW6acnmtoD6Lxo5+T0LWcJJp07VQwSMLgWbNmNfhzNKZp/9Q1Q4VRESzC1TVLRZQRed+n89gDDzyQcnqOI6kw7Qkq9SURUlVIknXBBRcUMY0LWjf0d5HQW+W1JCCknErINtxww1RD11PXRJIg0vlP1w0SuNG+r+Plr3/9a6rRczHtw7QmqbSYpM0ko9T7QPOzKiSDU7lWrZLQfv36FTFdc10jDj/88FRDUi6VIKpIOoLPIiq63H777VMNPQfoPkBC0LXXXjvlnnvuuSImWeNTTz1VxLRmqHwvIs8FEsruv//+KadznUSXVaE1WM/qtMfQuqxyWpp7ej548803U80xxxyTcnpuIakcnUFV2kmSa7p/ugY+88wzqYaey1UES3NPhaOtWrVKNfRzW265ZREPHTo01Rx22GEpp9LOxx57LNVURWWYEfk8qbL5CH7+U4E1iW917Ov9jeDnTz3P0vykfVbXJFozSICuewzte3SvdF+ldwVdu3YtYpJTUk7F19dee22qufDCC1NOxfTdunVLNVXwv7A2xhhjjDHGGGOMMcYYUxf4hbUxxhhjjDHGGGOMMcaYusAvrI0xxhhjjDHGGGOMMcbUBe5hbf4lUN/ImTNnFjH1q37jjTdSTvunUp8o7auovWEjuJ+e9n+kHj577bVXymm/1J133jnVfPvttymnfRypJ5z2Q4vI/ZWoF2GnTp2KmPrFUv9Q7a1Efa4///zzlNOemz9lH1DqR6Z9Bak3IPWZ1utJPdr22WefIqbenXp9I3J/u9GjR6ca7b0dkfv3UW/o9ddfP+W0dxz1NaN+v9ozbN111001J5xwQhFTLznqO689trRfYkTukx6R5//8+fNTTVVofgwfPryIN91001RDfZB1TFE/vTPPPLOItZ9zBPdo0/7U2u8uIvdJj8jjgPrfUl9M7Qmn9zyC107tG/mb3/wm1WjvaepxfPHFF6fca6+9VsTHHXdcqhk2bFjK6b5A63sVtLddRO5l16dPn1RDfQB1HenRo0eq0e9P6xp5G7T/HPkXqJee9l6kvYv25ltuuaWIaczRWjdt2rQi1r6WEfneUU9O6vmn+zf1Ln3xxRdTTq+L3oPGQH28V1111SJWd0YE7zna55l6a+o1nzp1aqqh89gdd9xRxLTf0D3WPsHUr5r6SOq9op7S9Dn1bPe73/0u1XTs2LGIdW5E5H7VEfm+aE/iCL4G2223XRFT/8uq0FlgyJAhRXzPPfekGvIR6LmNelFrP3XqRb3mmmumnJ4l6dxKY0rXG/WzRHBfWV239t5771Tz0EMPpZz21t1ggw1SjY5pcnzQuqzjhfZqcgto32XqvV8VWid1L6azAfXtbt26dRHT84P2p6V+1XT20XPVgAEDUo32io/IPfN79eqVaubMmZNy2nuW1nx1IUXksUE9+/V6kveIes/qOZXWRJrXkydPLmK6BlXZd999U073K1oj6L5rr2R69jrooIOKmOYCzetHHnmkiGnc0zqizyJ05qZnGO1VfPLJJ6caclOo/4Se2XQc0FmA+nFr33c6u9M4f/jhh4tY97PGQN9Pe6DfeuutDX6miOyGuOaaa1KNehpor6K+2tqLntZ8Oj9rv2hyG9AZUP1W6m2JiDj22GNTTl1k6hWIiBg7dmwR6/NaRH5XEZH3IfIB0NzTfe+nenb3v7A2xhhjjDHGGGOMMcYYUxf4hbUxxhhjjDHGGGOMMcaYusAvrI0xxhhjjDHGGGOMMcbUBX5hbYwxxhhjjDHGGGOMMaYusHTR/Et49913U04b76uMIIJlcCq9UXFhRJaXkZxo2WWXTTkVKbRq1SrVUDN+FROSEEFlABFZfnLRRRelGkLlQBpHREycOLGIVSoQEbHSSiulnIqqSPBFgjqFREdVIVGWyjZJsEXSPBUy0XXR8UKyBfo5HRskUiBRoYopVK4VwbIMlSnQfCEJhEpSaEyrMOTZZ59NNSR80LmgkrWIiG+++SblVFRFMqSqkORHZTK6HkWwOEqlNyQEVWEIrW1//etfU06lMCrHjIgYOHBgyqnUjKQiJJ6sRZ5DkjGdDyTPUcmsrq0RfF90nNPn3m+//VJOxWMkwakCjXudZyphjOD1Qcc0rZEqHDziiCNSDe1nuufQekhjXKWSJCuidUV/lwqV6DNFRGyyySZFTEK+559/vohpvyE50u23317EtP6qoCoii2mXW265VFMVEvdpju6LXqeILO/aaaedUs2nn35axDTvaL3faKONGqz54IMPUk6FnHr2+iF03NFcoDGsslyVp0XkMwtJ5S6//PKUu/TSS4tYx+EPfU7dr+l+VoX2QRW6kmCW5qyu9yRm1XPxlClTUs1SSy2VciphVTldBK9bKlQlQTHJ71QMRmLW999/P+VU8kiCTP2cdIY544wzUm727NlFfNNNN6Wajz/+uMEcicqqQtJOFZfSHFKBW0Qe+/R81LNnzyImGSZJ1vQsQLI/ks/pmUllnFQTkeWbJEYjeZjKx7t165ZqdO8ngSWda3RvonMOiVLbtWtXxDTuq0LXXPdHFfZG8Ppz//33FzGtbTofr7766lRDz616HiBhJUlft9lmmyJW+WcEC1b1rEx7Bd0rld3Rs+1ee+1VxLTv0vXV+UlnUJVTRmSpJJ3lqkJnfH1eIPn66NGjU0731dNOOy3VqICe9hyVx0bkPYbE7vQsdN999xWxrg8RfA6+4YYbipj2L1p/dK9QwWJEfi9A52B6r/P0008XMZ2b6Plg0qRJRUzPxFXwv7A2xhhjjDHGGGOMMcYYUxf4hbUxxhhjjDHGGGOMMcaYusAvrI0xxhhjjDHGGGOMMcbUBX5hbYwxxhhjjDHGGGOMMaYusHTx3xhqgK7NziMi2rRpU8TU8P1/m8UWWyzlmjVrVsQkTyRZhTadJwHD559/XsQk36PfrQIoktmQyEBlSHTNSSyg0gkVBkSweEi/M30XlUmQCIjEHypbUqlSBEvdVMShconGoCK/iCy46dq1a6ohAYIKJkjyo9eTalTaEpEFVyTiINGISjVeffXVVKMyrYgsACRJIN0rlU7MnDkz1ajAg8STJKTTuUYCoVVWWSXldP7TelcVkqjp79d1M4Ln/0MPPVTEvXv3TjUqUlHJUQSLXPR3H3LIIanmqquuSjkVi+jcj+BroNItkrd26NAh5a644ooiHjduXKo59thji3jNNddMNeecc07KqcSIxE4kIG3fvn0R09w79NBDU64hSCCiYiCaY7QHqCiH9iVda2hfJFGYrkckSyPZqQp+vv/++1RD11vFPDS+SOKp0lmSJw4dOrSIacyTSEvHAAkdad3WMwV9pqqofDQiy9jovtB6pN+HziIqeqJ5R/NcZYkkhd1yyy1TTuU9dNYj8dGBBx5YxCr5jWDZnY47EsbpdSHBdufOnVNOrxXN4UUWWSTl9L7Qd6nK3LlzU07nqJ6FIlh6qPd0zpw5qeaJJ54oYpJT0llAzye0b9AZWM8id911V6ohadS9995bxCQsV5FoRETTpk2LeOrUqalG1zL6vjQ2F1988SKme0DzUc+uJGasCu0D+uxFe6U+Q0VE7LjjjkVMAnGdVyQkpnVE5zFdc/pMep6n50ES3usaOG/evFTTvHnzlDvooIOKWEWmEXmNouezESNGpNwpp5zSYM1ZZ52VciqDozFdFVq79fxOYj1aO1XYRjLB6dOnFzE9P5BIWEWwRx55ZKqh8aPPTCQEpTMYPTMpJC/U5w6SMj/44INFTM8PtD/rXkyfe7311ks53QfovlSF9md9p0Hn2T322CPldD7SnJ0wYUIR9+/fP9WQyFPP8D169Eg19F1UVNi9e/dUc8EFF6ScvrvbfvvtU83dd9+dcvqspWt5RBY63nzzzammS5cuKafiST2HR7AMXMcw7eEkkG0I/wtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXeAe1j9Anz59ilh7iUZE/PGPf/zf+jjIxIkTU057ikVEPPPMM0VM/Y/+t6F+VtqXeIUVVkg11BfzpZdeKuJZs2alGu2RRn05qUeu9guk3lzUW017RVF/qa233jrltOfU6quvnmq0F2xE7t1EvZW09yD1tdXe2xG5by/1/KQ+jtqnjXpzVUV7v0bkXtDU76lly5Ypp33MtI9vRO77p329I7hXnvakpR5bu+yyS8pp/6z58+enGuq52bdv3yKmuUD9Z7WHGN3PTz75pIip7yn1etVrRX0kqS+79s/6KfuA6poYkecM9W2j/o/aD5F63nXq1KmItd9mBN9j7Rf/5JNPphrqW6s9/Xr16pVqqG+b3ncaY7reRuS+x9p3OiJi7733/tG/FRFx3nnnpZyuidSXnXq2ah/ZIUOGpJoq0N/abbfdiph65FL/c+1PTWNc+5e+8cYbqWaNNdZIOV2naU+gfUl759G4vOOOO1JO5/4LL7yQarbbbruU0z2d1sPrrruuiKnXN/W21HV00003TTW6rkXkvsu19KesFVp/9b7TvaLzkI7x9ddfP9XoHCYfAvXP1z2W/v7777+fcnqO23///VMNOQvUkUD3U/eEiDyGaZ7pWYA8A9rrMiL3WKW5Tz1d1SVBZ62q0HOHriO1rtvav1R7Lkfk9XfMmDGphjwcei6ePHlyqtF9IyLi0UcfLeKOHTummmWWWSbljjjiiCKm85iu0xHZrdKzZ89Uo2OT+nGPGjUq5XRtobMlnad1PtK+sO+++6ZcLdC90mtAvaipz7SOO/Lb6LMQ7fvvvvtuyuk6Sc9CdJbUuUbnI5pDupbSGKPnT+0PO3z48FSjZ2XqnUzr5FtvvVXEuh5F8D6raxKt3VVRt0FEvsd07chloL2S6X5qP2M6I9IervsQPZuQ0+jUU08t4lrO8xG5LzE931NvZr1X5Gei/VKh9zr6DE696Wmf1TMgeTeqQmdHXePJ4UHPY7pu0HXSPvA0Vn71q1+lnK75dO/ouugaQfuQnmcj8ppLfa7JgaD7FY07ffdDLiTyl+n8pHWE1iRdA+lMXwX/C2tjjDHGGGOMMcYYY4wxdYFfWBtjjDHGGGOMMcYYY4ypC/zC2hhjjDHGGGOMMcYYY0xd4BfWxhhjjDHGGGOMMcYYY+oCSxd/gAsvvLCI27dvn2oOOuiglFNByU+JCkNIsEj8lMKXnwpqHq888sgjKXfaaaelnAqYnn322VTz5ZdfFjHJy0gIpZ+ThAj091Qc065du1RDQrx11lmniEkSQxIRlS2RIEDHyxNPPJFqqKm+So2aN2+eam6//faUU7EmCWeqoiKyiHwNDj300FQzaNCglNO5TUIWFcWQFIJEono/SY521FFHpdy4ceOK+Oabb041JOKZNGlSES+//PKp5umnn045lfyQKEI/O8krCK2jMfanP/0p5V588cUiJpFdVUhKdc899xQxyYkOOOCAlFM5xx577JFqdL0hadOuu+6acipIo3lNoiNdp0iQSZIWveYq/YmIOPDAA1NO60hqq/ITEsCQQGjDDTcsYhIfkYBq+vTpRUxymSq0atUq5XQMkPiExM0qkqK1R+cwCUpJmKKCOpI6qVQ4IuLTTz8tYhWzRUScccYZKXfYYYcVsQpbIvIci6hNtqwyJhIHUk7ld7RmkgRM56fKixoDCYNVSEXnDNrzdO6TNFnHPX0Xkp6pDIqkWR9++GHKqZiaxj2tkbp/0jmHzhD6uWhdUcmiCtYiWOaloieVzEXwnq57Fc2zqvziF79IudGjRxcx7d8kHVP52zvvvJNq9LvQ96Xvp4LpNm3apJorr7wy5VT8qOLUCB4bm2yySRHTGqzno4gs6iIx9ssvv1zEdDag+aH7Lp29SOCm16V169appir0PKhrtQosI1iwWssaodeOnoXofKK/mwR5Cy2U/82d7nN0Lh8/fnzK6bpM+wmdv3Qs9OnTJ9XotatFQBiRxzmdD0jgqOcvXZMbA50j9HlFn7cjeI1Yc801i5iuuUrlSMxI67lecxLk6fkgIl9PElMvtdRSKaf3ne4Vrd369zbbbLNUo5J2lSJG8FldfzfdA5IU65pE54Oq0Pqj94/2Kjp367mX1hY989H3nTJlSsrp8xFJX/WsHJElurRu0pqkUuRzzjkn1ejzWUReF0m83bZt2yKmeUafU9c7ehagdUvPEd27d081VfC/sDbGGGOMMcYYY4wxxhhTF/iFtTHGGGOMMcYYY4wxxpi6wC+sjTHGGGOMMcYYY4wxxtQF7mEdEf3790857fnUq1evVHPZZZelHPXg/KkYNWpUEVP/LmLGjBn/jI/TKLRnbkRtPbqvuOKKlKult632p6Zrt/POO6fcgw8+WMTUL4z6cWs/RO3HFhFx8MEH84f9H1CvM2LmzJlFTH1ttVcf9aWi3lxdunQpYuqLTDntV0q9V6vy+OOPp5z2m6RxQb2LtV+W9lCMyD2n9Fr+0N/TnlPaSyqC+5Vqn0jqnffaa6+lnPYrve6661IN9SB/4403ipjGnfZq33vvvVPNH/7wh5TbYYcdipj6nlJvNbpWPxXXXHNNyp1//vlFTP3mdD2IyP26qI+s9u+kXrra6zEi93t77rnnUk3Xrl1TTsfw22+/nWqor2GPHj2KmPpGUh9O7elOPeiOO+64IqbettrnOiKvbfoZI7gn5cUXX1zEq622WqqpAq0P6rPQnrIREZtuumnKLbPMMkVMvWd1vlBPRVp7dF2h/tEDBgxIOe1xTH1faf1t1qxZEZ966qmphnrZ6b5LvR+1Rzj9Hhrjeq/Gjh2bak4++eSU0x7d1Le9KnrPI/LaSuvvL3/5y5TTHqC0l+haQ2cvOt9qX0k6M9Gc0hz1rFxrrbVSbuLEiUXcsmXLVPPYY4+l3Pz584uY/BK6n955552pRl0oEbl/MvXBp/upY5H6PldFe/pHROy3335FTL08b7rpppTT70d7nvbEpB7hX331VcrpmeWGG25INdQ//swzzyxi6vdJZx91FtA1p/6w9913XxFTj3ntnU7nAPo5XUupDz25bH7zm98U8eWXX55qqkJnDz23UU/gE044IeX0fKDrdETuOUxr8Lbbbptyeu3o3EjjQPfCbt26pRrqxTphwoQGa+he6d/T825EnleDBw9ONfr8EpHn8SGHHJJqaGzo/qh+iYg8z2pFe+1G5OdN2itojVCPEnnB9HxN/YbJMaZjk94d0H659tprFzGNO+qjrWOf1i393RH5vpMTS3uZ33333Q3WROS5Tn39yd2l/YynTp2aag4//PCUqwV6RlR/Dz1jkHdMn+fJnaPPOeShINShRmNFHVUR+fmP/Ex0zfUdET1nkKvnoosuKuJp06alGnWq0XpL5y19x6nnjAh+lj766KOLmN5fVMH/wtoYY4wxxhhjjDHGGGNMXeAX1sYYY4wxxhhjjDHGGGPqAr+wNsYYY4wxxhhjjDHGGFMX+IW1McYYY4wxxhhjjDHGmLrg/5x0kZqkX3LJJSmn0hsSAZG85qeiQ4cOKUdyP+WTTz5JORWG1AMkmLjrrruKmKRgJA3Q5vQkfFCJAEkJSWj2yiuvFDGJpOjnVMRx0kknpRqVBUXkcfb111+nmldffTXlVGxEQsX777+/iFu0aJFqSF6hDfpJujF79uyUU9EIjc2qkPBBfz8JGKj5/5dfflnEJJBUGRqtIyTZ+PTTTxv8+yS40fmhUooIHncqGlKpSQTLpVQUQ2I3lUKoUDKCxRQ6hnVORbB8ScUm9F2qQhIaFT+qdCOCxTEqUiGh1+qrr17EJK1TAVZEXkdIukHiM5WCDhs2LNXQNVAh0gYbbJBqSFqiaxmtW3o9v/vuu1RD4tu33nqriFXCE8EiMv1Mt912W6o54IADUq4h6DOq0JAEiySM07WG1izdE2qRu0ZkCRgJHWne6WegfZj2CZXijBw5MtXstddeKad7c5MmTVKNCmFI2ElnNJX7kbCTZDMqzqG1vSokZdVzDc0xkve0bt26iElsp8JIWv9JhKTrGt0XkhBqjs5jNBf0DEPCVxWxReR9gu7VuHHjipjEaHQmVLkeia0op/Oa5HBVoXVb5Wx0vqa5p88+JPbUsx2tPyTR1HtMYk86Z2iOzpsrrLBCym2//fZFfPXVV6eanj17ptzmm29exHSW1TMwCadpfurZjtbgjh07ppye31V83hhIPLnHHnsU8ahRo1LNM888k3J6H0gqqbJlkrDSuFMJ4rfffptqSCCuguvbb7891ZBIT59X9Lkggse5SkhJ+qo5qtHxG5GfafQZOYJFqUOHDi3iI488MtVUhT77TjvtVMQke1966aVTTtdhWoNVMEjr7TnnnJNyW265ZRHTuKf96/TTTy/iESNGpBoSOPbu3buIScLaqlWrlFt11VXbpRQbAAAgAElEQVSLmASHOg5oLiy88MIpp2cyOsvRmqRndZL4VoWeV2bMmFHEJJ2mdz96dmzatGmq0ftAz5+0tumzLF27TTbZJOV0/yDBos7PiCz6HjNmTKohKamexV9++eVUo+8raJ7RuVT3Rjq3qfgyIq+59HzWtm3blGsI/wtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXeAX1sYYY4wxxhhjjDHGGGPqgv946aKKYyZOnJhqqMG7yrVIevGXv/ylcR/uf6CSnVtuuSXVkLxG6dOnT8qRyOlfDclPtNE/yR3efvvtlFt//fWLmEQuJCRQSAaw5557FjE1pldJQ0QWdsyZMyfVkMBNm9Wr6C4iS8gi8mcnQZLKFVRKEcHCIhVzkfCKpE0qkSPZU1VeeOGFlFMBggoPa4XEJirG2GqrrVINyeAUlWBERPTr1y/lVORCAhgVAUXksUGyp27duqXcIYccUsQkwdExRvITkkrqtSNhEs0rHWcknqwK3SuVwnTq1CnVkLhPxVwkl1KBWK3rtK5bdM9VShORhRYknBk9enTKPfXUU0X8wAMPpJoTTjgh5VSqodckIn92Eqgdc8wxKaeyzSlTpqQaktDoZyI5URVImKnrHwmZSa6qazB9Rv1u7733XqqhsarQHkTSRZVhTps2LdXQ59T7qTLkCBaoqcSS5KMqHVJBTUTEBRdckHIq/6S9k+Sq+v1IgDNw4MCUqwW65rp3zZo1K9WQVFLXYBU/ReR9l85Hek6OyHsHrVkXX3xxyulYpLHy2WefpZxKxWuVI+n5ncadStZU1hTB0jOdCySoojOT7p90z6tCZw+9BnTN6UyoUqWNN9441bzxxhtFTM899AylZ3ySkNH5ViXxJOsmqa0+G5BYc968eSmnoi4S3qtEjqSohAoASSRKgkM9/+222241/b1aWG655VJu8uTJRbzffvulGpJ56dpC51SVvVMN3RedQyT7JPn5008/XcQq047gZwN9HltvvfVSDQmX9UxPZ2edVyrFjuB1RM9fJJ4jofZJJ51UxPSsWRUVLEZkUSDdFzoTqVyU1mUVgurciOC1Rddzuuf0XVSoSNeO/p6e35s1a5ZqdC5EZJEnyf10zuy4446pZubMmSmn45XeqdB6ruOVzh5VIfG1SjNJyEfnXpXZ0z3+/PPPi3jSpEmpRs8eEVk0ecQRR6Qaeu+hY4Pui0qZI/KZr0WLFqmGzj96ricBu6439OxO68+jjz5axBdddFGqob1/3333LWJ6L1AF/wtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXfAf1cOa+iZpr1vqHUMMGTKkiLUv1k+N9uujXnnKRx99lHLa77NeoV7UrVq1KmLqZ0V96bTPjvadjoh49tlni5h6UFEvJ+2/Sz3+qB+Z9mSjvm3UE077lVJv3z/96U8N/j3qWau9danvK/Ur1X5H1Atxgw02SDntWUa9uapCv0vnOvUQo7652td+r732SjXaO2r48OGphsaG9lqkHnjUA0p7CVOvNe3/FhGxxhprFDH14aMen/o5aWzoWKSegtTbUXuhUv816mmqPX9pLlSFrov2IqR+0dTP7oADDihimp86FqmvGfXo1r7h1DNNe3BG5B6m2oc0ImLRRRdNOe0BR+P1hhtuSDmdV9QTV3utde/ePdU89thjKad9h8nl8MEHH6Scruc07qpA11u/v+5JEdx3WXupU59/7VtJ/XAvvfTSlOvYsWMRz58/P9VQ3zqdw9Q3k8aT7gvaFzUionPnzimnaxb13R86dGgRUy9Y6uOt/ampx7P2OIzI84D6F1ZFzyIR+axDfVBp3GuPburXqnsz3XP6TCeeeGIR07ynea49R6l/ofbsjsh976mHNfkPdt999yKmXpd63tPe1BHcG1/7qa+55pqphvZB9S/Q+jRgwICUqwVaN3VfoD6d2hs6Is/Z6dOnpxo9L9D5gfrO69ign6Mzvu7zzZs3TzXUj1971j788MOphp4fZs+eXcT0HKk9nelZQdeoiIiPP/64iA899NBUQ+cF7fX/U/pf9DNF5Oup/ocI/px6b6hHrp5haB+iOavrHc2zn//85ymn95P2Ba2JyL2EL7nkklRD/dR1PaXxqusN9W6nPU2fLamX8EEHHZRyet5UN1JjGDlyZMrpWn3UUUelGjo7q6OFzqT6LETPmvQspM8PNFbI59C+ffsiprMO+Z822mijIqY1gvZLfW599dVXU03r1q2LmM7JtM/qnKX3ELp/RuSxOG7cuFRDLptaoPOdPs9TL2p6ZtN+/HfeeWeq0fMIrW30vk33Qj0PRfA+q2OD9iE6B+u5Xs/FEdx7upbzur7P1PdaEbxGqONi/PjxqYbmwoQJE4qYnnfpPV1D+F9YG2OMMcYYY4wxxhhjjKkL/MLaGGOMMcYYY4wxxhhjTF3gF9bGGGOMMcYYY4wxxhhj6gK/sDbGGGOMMcYYY4wxxhhTF/zbShdJjHbbbbelHAnNFGooP3DgwGofrAZOPvnklCMJoKIitO22267BmnplqaWWSjkVwKisIyJLqiKynEOFNxFZSLDiiiummu+//z7lVACzyy67pJoZM2ak3BJLLFHEJCggcaAKJUgm0a9fv5RT4d/zzz+falQmQ4IUlb1EZCkEySToeur3I7FTVUha0qFDhwZ/jqRCKlJSEUhE/n4kb6XrqUIWlRhEsMBRxzR9broP+tn79u2bakhwo1Iokl6omJHkFSTU0XlFQjoVnUTkuUDSnaqQ+EPFeSrajGDpjX4uEhWq0GKZZZZJNSTWfO6554pYxVIRLAtTicgDDzyQami9a9u2bRFPnTo11dx3330p16dPnyLWNSMiokuXLkWsMp2ILD6KyOsWyU3XX3/9lNN7RaLfKtD8UbkNrUUqAIzI94lEvCrFIcEinVdUekvzdfXVV085FSjROnP55Zen3FZbbVXEtNaRyFRF0S1atEg1un+SIIbOeiq7IpGfigsjsohR5T6NgeRzKtyh+dOtW7eUU0ET3U8VhtPfJ8mjXnOamzSmdJ6RkI/GlK5HJBUncZYKfejZ4NRTTy3iAw88MNXQ2q7nLxVrRbCITUVP9F2qQmdgPTvPmjUr1agIKSKPF1p/VIRG35ekpPqZ6MzWo0ePlNPzAUmkaP9WcWCbNm1SDe0TL774YhHT9dX5SZJokq/rnFERXATLVLt27VrEtB5Uhc4ehx12WBGTKEyvU0QWENOerms+PaPSvFJJ38EHH5xq6Cyp6xuNOxoHOs5J8kbrq+4xJPHVz0Tfl/a0QYMGFbGesyL4PcC1115bxCS5rQpJLHUPJREbSTv1mYnGjz4H0L2bNm1aymkdCXPpWU+FdHQ+IPGj7nMkuaV1Y968eUVM53Kdj/QsRM8G+sxG8ju6drpG0BpcFZJYtmvXroivuOKKVEPSchWp07OBvt9TSWFEPrtG5Pmoz0YR/MyvdR9++GGqoTWwZcuWRXzOOeekGpKSqpSexp3O/yeffDLV0BlQ1xYSw9L5Ts/LNDar4H9hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXfBvI11UAcKNN96YamoRLFKz8xEjRqTckksuWcQkiakFkvSdf/75KUfN/xWVXZH86d8FElP84he/KGISuKkcLiKiWbNmRUySBIXEBiRb0Wb8Sy+9dKpRSUNEliTQ+FEJWURujk+/m4ReKk4gYYjK/UhqSVIIFauQCJLmlUpTSMJVFZVaRuTxQnIOkmWoDISkQp07dy5iklPq943IAjG65pTT70fCEpICquCGBDskX9JxTQIhlQORbLR169Ypp5ILEteooCki4ssvvyxikudUhdZOlduRJFQlmhFZRPHxxx+nGl3fSZRz//33p9zRRx9dxCSNo2uu6xvNWV03I/K41rkfETF+/PiUUwmfCkQislSSBB4kF1UZyMSJE1MNXTuV7JDYswoq4o3I5wUaJypYjMj7AsmKdB0jYVMtQiz9jBEsotQ9gMRWJLZT6O+RXHqNNdYoYhIv6fXcb7/9Ug3JX1RoRJ975MiRKafSM5KlVWXBggUpp2JRkspNmjQp5XTfJWne/vvvX8Qk3CHx+FlnnVXEF110UapRmXVEvsd07Ujqq3sOyYNuueWWlBswYEARq2QyIu+fjz32WKqhOavrNt2XVq1apdy9995bxHq2bQx0dlaJ+WabbZZqSEo/dOjQIiaZl+6Va6+9dqpRwWtEFnLS/k0CX13faG2jOaT7IM1rWltmz55dxCSo03PjmDFjUs0pp5yScr179y5iknvRmqhn9ZNPPjnVVIXGvn4uGud0ztDzAj3vqhSMBG56DyLyfdA5FcHCaRWT0Zg+++yzU+6oo44qYpKz0XOcPnvpuTUij1faw48//viU69mzZxGTpFSlxRERF1544Y9+xsZA+7reU3qWJnGoriWnn356qtHxQ+9+VJgekfeYYcOGpRqV8Ubw2qLQNdcz98UXX5xq6B2VPlfRdVK5IL0XoHOpzqG//vWvqYb2cD1L9erVK9VURSWlEXnvPfHEE1MNnd/1OuhYiYjYe++9i5ie0++8886U22effYqYntN1fkZkWSN9bnpfojLKjh07pho6o+g5n8aPrtMkdKTneRVG0ucmVI5Nz4NV8L+wNsYYY4wxxhhjjDHGGFMX+IW1McYYY4wxxhhjjDHGmLrAL6yNMcYYY4wxxhhjjDHG1AV12cNa+3ZGRJx33nlFvMIKK1T63YssskjKUf8a7eGjvRDpM0VEPP/880VMvTVr6bX9zTffpNxpp53W4M/9u0DfT/vsUL9h6sGkfbCon572LaIetrX01qW+ZjQ29txzzyKm3tvamzUi94lbccUVUw199oceeqiIqdfQ8ssvn3KK9h2MiPjoo4+KmHpAUf9Q7SVH/ZeqQj3fH3/88SKmHsQ777xzyumaQH2ftb8UzWHqub7yyisXMfVqpv5yeq+WW265VEPjR3sR0pimfmRad8kll6SaPn36FPEf/vCHVEO9trVHG/VMo2unvQdfffXVVFMV6i9Xy+/XfmgRuW8szQ+d17vttluqod7XV155ZRFTv1RCe9jTmNb7EpF7AdKaQT0wtZcbjQ1dA6lf89y5c1NOe9LeeuutqYbm1aabblrEkydPTjUnnXRSyjUEra3a9/D9999PNdS7T3sX07qmPemotx19f+0jR/1Fad7pd6E9l/wAuqdTv2jtiReR+/tSj1NdM2gPpL6OOsapzyP1Jf3jH/9YxOqkiIgYPHhwytXCvvvum3JTp04tYurxTPNFz5t0Phk7dmwRk/9h1113TTntdX3wwQenGuqrrT3QyYtBvXx1faC9mnqVPvfcc0VM1069FPR76JlC1yM6H1FfZO0rSc8B1BO8FqiXpo7Pbt26pRrt0RsR8eCDDxaxnvUi8p5D/gVdayPyHKK95K677mrw52htoz7T1157bRHT2Zm+n0LPiLpu0TOq9p2OiLjmmmuK+He/+12qoTV4jz32KGLq410Vulf6GcgDRP4e9TLQuVj3OHoOoOdB7e1NfdnpnKH3j+4n7ft6dqU1ivYP3Y9pv9Q+/jSmqY+39rqmsyXN9VGjRhWxum0ag64ZEfl8Sa4Vun/6c+RX0H3hhBNOSDV0JtSe0kOGDEk1dOZWyL2kzoCI7Fdo3rx5qiH/gJ5v6e/pfNQ5FZGfNSNyr30am/SeRc8I9PxZFVqDdW2ZOXNmqqFrrs/OdEZRV47uExF8pp41a1YR05lX51lExDbbbFPE9H6InBa6BpPzhq6drufkmNAcjQM6m+v6o/tZBL+H0DMKnbeq4H9hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXVAX0sXjjjuuiElMoQK3fzbadF4lVhERN954Y8qpTIEkSrVAEpx333230u+qR0j8o2IBlW5EcPN/veYkAFShBgl2VJQYkUUN2oQ+ImKLLbZIOW2OT8IrEj6o5Ie+76mnnppyKi0gyZqKMEi2MH78+JRT2RPJn1SUQ6j8oDGsueaaKaffR2VlESw/0nFHY4NEBgqJMPSa07WjtUWFCyT7mzFjRsrp5yRJDMkcVFxFwhCVMpDUkq6Tyk/od3/55Zcpp9IkEjpWZejQoSmnAo1zzjkn1dA4V4nQKquskmp0/1JZWgRLPfRa0fykufDEE08UMa1bKiKLyOvGzTffnGpIHKNiLhI76VraqVOnVKPCkogsAfz1r3+dakg217dv3yI+++yzU00V6LupDIlkp7T+vfXWW0VMsisdTyQcIymXCuJITETCpjlz5hQxfV8VBUVkKTRJz+gMo0JbkgftuOOORaxjIiLLTyPymkF7AgmpVMD6q1/9KtVURedmRBbutGjRItWQKEfvH8ko9TxEa8H06dNTTvclkoLRuq3SQ5JI0Xlevx+tkbRm6b3af//9U41eT7qWJDTT+/LOO++kGloP9exKAtKqkLB3jTXWKGISCJPoSfeTCy64INWo5IyeaQ488MCU07MWzVnKLb744kVMaw1J3XTc0X5Kn1MFubRGHHbYYUVM4mGSVuk5Q0W4P/T3VFBJz0JVWXLJJVNO12o6p9K6pcJaPetF5LFJc3+ttdZKOR1n9LtpHKgE7O233041KvSOqG38kDBcrxWJA88999wiJhEtnS1VPEvzmuajnglJAl0VOg/oPkCSNXqW1fMOCTlVtkdnHboG+plUbBzBEnM9k5Egk8Sz++yzTxHffvvtqUafsyLy+KH3AiphnTJlSk2fqUOHDkV87733phoSDqrE96WXXko1LVu2TLlaoOs5bNiwIqb9kvZsPffTeyU9P+v7xoiI0aNHp5xKQkmuTOcBfYahcz+NYRXI0ni94oorUk6/H52xdW8iQXn79u1TTr8z3Rf6frpOkjCzCv4X1sYYY4wxxhhjjDHGGGPqAr+wNsYYY4wxxhhjjDHGGFMX+IW1McYYY4wxxhhjjDHGmLrAL6yNMcYYY4wxxhhjjDHG1AV1IV1UyWJVwSJJyB5//PEifuSRR1INCW5UoEayDhIL/FTccsstKTdp0qQi/uCDD1KNyh0iWFbxr4YEDCoqfPTRR1MNyTm0Wf3GG2+calS8scgii6Sa1157LeVUvKECkYgst4qI+Pbbb4uYxgoJJlQUQdeAJH0q1aC/165duyImIZ+KBiKy+O33v/99qiHZ08SJE4tYJWiNga6dzm2Sac2aNSvlzjzzzCImAZXOIRIOklxKpR4kGWratGnK6fpDskgSYaicR6UtESxJuOeee4qYxE4qblDRZwQLAFVmM3ny5FRDUiGdCyS8qQpJbx588MEiPv3001MNiUpVdkLiPJUZfvPNN6mGxpSKsjp37pxqnnnmmZS74YYbiphEfbSO6PwnoVi/fv1STvciGps6/0nKR4Ik3Z+POeaYVEPrlt4/kk1VQaU8EVkqotKsCBbGqdSNZFC1iJfuvPPOlNN9afvtt081JAHTz6T7a0QeXxFZnjhhwoRU8/rrr6ecrkcqWIzIn33fffdNNbQv6R7wxRdfpJrVVlst5fS70NmyKnT20HMGzXO6nroP0rqm8h4ah7QHfffdd0U8b968VNOsWbOUUzkSndmuv/76lGvbtm0Rk1iKrp2uD7T+6n2nsw+Nc92vVQgYwXuerod0bqzKHXfckXK6btx1112phu6fisJIbDVw4MAipn1Dzw8R+V6NHDky1dCZW59hSKKusr+IiO7duxcxCcZuu+22lNP5T9JkFbrSPkVrqf4uEjPSuqzzcdq0aammKjT2dXyef/75qYZk2SrvojOwnv9oDSZRqj7DDB8+PNXoXhWRRbs09+hco89sdEaj5w7Nkbzs0ksvLWI6j9F30XchJGakM4OKk0kkXxX6XbrH6FyMYEGd7jH0rkK/H+3zundEZPlcz549Uw0JT/W70HlTr29ExFZbbVXEe++9d6qhs6M+55Cw+8YbbyxiemdF50K9VyTMpDODCohJrlwVetYbPHhwEavoN4Kf3VWwSlJ6XUtJ1Nq1a9eU0/VnoYXyv++ldwz6c7TW0PuvlVZaqYhVmBnB64+KPGnd0vuuEs+IiEGDBqWcjmnad0lGqe89evXqlWroWaAh/C+sjTHGGGOMMcYYY4wxxtQFfmFtjDHGGGOMMcYYY4wxpi7wC2tjjDHGGGOMMcYYY4wxdUFd9LCupcey9pIbMWJEqqFeQ9T/TKG+1gr1CR02bFiDP1cV7WFLOeoJdd5556VcPfawpr5Jiy22WBFTv0Lqw6c9pxZeeOFUo713qJ8p9QTX3pLUw4d6wn388cc/+nsiuCev3lP6vtTHTHuiUS8n7cFEvcjoGuh90H6tERG77757yumcpZ6CN910U8rVwq233ppy2seMrpP2VorI4+e+++5LNdqPjHpAETfffHMRU99e6qm18sorFzH1z6I+eNq3jXrJ0T3WvpTUR01/juan9iuMyD1TqW8l9XHUNaKWtbxWqE+u9hmke0X9X7VnKvXr0h641OuR+pNpr7Mjjjgi1VD/5nPOOaeIr7322lRDfZZr8QFQb7VRo0Y1+HO6vlMvS+rxqXNW+9ZFcN9anTPah7Qq2msuIt8n9TFERMydOzfltIcs9VPWHnHaazMiYujQoSmnDgHqHUhzWHN0L9UNEpHnC91L6hOsv//5559PNS1btixiOrORt0HPAtTfj/pT6+fUHvSNQft2RuS9mHqS03qvY4POjdrTdb/99ks1dB47++yzi/iss85KNXTNt9hiiyKmuUBzWPcOWuuuu+66lNN5Tn0lr7766iKmNZPOjTpeaO4//fTTDf4uOqv3798/5WphiSWWSDnd9w8++OBUQ/unnkGpH6X21h0yZEiqIeeOjnOdwxG8x6q3hfqGUw/QJ598soi33nrrVENrsM5tdVJE5Hn18MMPpxrtex+R91haf+j5Qft70jn5sssuS7laoP7b2uf5oIMOSjXkPtFzODk+dP+mM5SekyPyGYauE7lWtA87/Rw5AnS/onMq9afW/ZjO1+pcoHHfqVOnlFN/hbqKIrLnJCKiTZs2RUzPg1XReRYRscsuuxQxnTfpmUJ7OtP11ec48iaQQ0nXXJpD1N+8FucOPY/pc+Txxx+fai6//PKU69GjRxHvueeeqUbHC/Utp+dd3dfpXRCt+Xo9u3TpkmqqQuuPegvoHQPNR+2rT2NDe0jTPae9Uft20/pH612TJk2KmBwXdAbUuUAuBXq/p8+t2tM6Io8X9QZG8LOAvsN58cUXUw2d1/U8S3t/FfwvrI0xxhhjjDHGGGOMMcbUBX5hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQvqQrrYr1+/Ij7jjDNSzZFHHlnEKhT6qVH5kTY2bwwkSFJpwaxZs1LNyJEji3jChAmphsRd9QgJv7Tpu8ozIliusOuuuxYxiU1UVEOyIJKtaMN8kmKRiFFFHCQMIHHMWmutVcTUrP7uu+9OuZdeeqmITznllFTz4IMPFvGzzz6bajbddNMGPxN9FxUURORrR3KQqqhsISKLUY8++uhUQyIwlYhsueWWqUbnlUqOIiJ+/vOfp5zKMqZPn55qSOChsgwSppFoTmVPJDkjcYKKzyZNmpRq9DMsvfTSqUZlQRF5nKt8L4KFVypNUZlXYyDJhspPHnvssVSz0047pZzWkQhIhTO0ZpCEUIVe22yzTap59dVXU073DxL60FqqYlQSfJG8Zu211y5ikgrpPCNRq8pNI7LgiiQ4JFpTodhPJR4i4Y2u9yTpo/VPx/24ceNSjQqiSDhIa7mK5Ug0TGNOofVez2wRWZxFcki6djpWaP1dddVVi5jkgjQPdB2lfYoEP/oZ5syZk2qqoiKtiLz2kISMxN877LBDEZP4TQXBDzzwQKqh85juZ2PHjk01KveKyHsc7bk0FlWORD9H65+uRyo9isj7G4nYVHQckaVDtG/QmVvHJ8mIq7LVVlul3EUXXVTEJ5xwQqr52c9+lnIqcdNxGJGvFUk099prrwZ/juRMJArTveSrr75KNSo4i4g46qijinjEiBGpRp8VIrIQneaZPouQMJ3Wd33Wa9euXaqhfUnPOvfff3+qqQoJc1VERgIuEq/p3kDnBR37JCEjyaOu57S+q9w+Is89Om+SvHDRRRctYnqOpH1W5WgkF9TxSp+b5MY6Z2+88cZUQyI2fTdA+3VVNtxww5TTPYX2ZxVPRkQcfvjhRUzvNPT+kVjvtttuSzl9RqT7QnN9o402KmJa8+nsrLK73r17pxoSjvbt27eI6byl6xa9O3jhhRdSbsCAAUVMc5++i0rpadzTc04t0N6rz8Ak2tV3gBF5HpM8WvdwlVBH8Hqu51e65nSe1b9HAtKuXbumnI5z2r9IGKkiaD0rR+T3AHTeo/mhou1hw4alGjo36Z6t75Cq4n9hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXVAX0kUVzJBw5n+b0aNHFzFJlGph4MCBKfeb3/wm5b777rtKv//fFZIfqTihbdu2qWaTTTZJOW28TzK6WqSWv/zlL1NOm+qTiEOb3kfk5vTdunVLNXTPVQBFjf5JUKeiCBIz6nchyQcJblRG8tRTT6Uaku688sorRbzBBhukmqqQSEWFbS+//HKqIRGYSjxo/Og9JlkZXTu95iScoPupUgQaKyRUVInRmWeemWro2s2ePbvBGhUAkkTllltuSTkdm61atUo1JK/Qa38qIbQAACAASURBVL7CCiukmqqQVEjlKltvvXWqobG/5557FrFKNyKyZIwEWDQO9HeRdEzXtogsKCHJLInPdN8jwSH9Lh37tF+qrIdkoyQo0c+kIpmIiEGDBqXcySefXMQku6vCvHnzUk6FOyRxIaGkiqTonvz+978vYpLRkSxXJbvXX399g38/Io8nmiskD1I55vbbb59qSHalwmASo6n07PLLL081KiCMiHjrrbeKmIQ7JH7TtYZEU1Uh+ZxK8kjYdPzxx6ecCppUmkq/i84rM2bMSDk9s9A9Jxmwrg8tWrRINSSx1HGue1IEC4PXW2+9IqZ9WNcaOuvR3/v+++9/NI7gvVmlujSvq0LnmnPPPbeISYREn1NlSLQm1yLN0zNiRJbs6j4ZwQJQFYeSzJrGsP4uEiTT/L/zzjuLWCWlEXn9IUEV7ek6h0hGPHz48JTTe0zSRbqetUBri543dU2OYDmtiibpM6mwjcYKnf90vJJckM4ntUD3SqV1JG+lfVbnOp0b9XmB1j9dxyKy/PzUU09NNbR/de/evYhpbaP9uRZIMqsCviFDhqQaOrvqMwQ9e+ncUylrRMQpp5yScro30tmdJKj6c/SugH5On+1UAvtDuc0337yIaZ4NHjy4iPX+RvAz2xVXXFHENH47dOiQcrp2k8ywKiRrVPklSYppzkyZMqWI6eyo5xEav7Ru6ZpLMsOmTZumnO5fdO6n66n3ht6R0dqpIk0am/o5SbR5zz33pFz//v2LWM9oEbze6fsnFWpXxf/C2hhjjDHGGGOMMcYYY0xd4BfWxhhjjDHGGGOMMcYYY+oCv7A2xhhjjDHGGGOMMcYYUxfURQ/rfzXU57V169aVftcHH3xQxNp7KOL/Xr9qQvv8RERst912Ray9q34op72/6Ppqr07qeUU9ZLVf6dJLL51qqG/vY489VsTUO0p7+0bkfsrUD5d6eGlvPuq/pj2RqPct9VieO3duEdPnpt7eBx54YBHPnz8/1VRl8cUXTzntPUb90Og+aO+vMWPGpBrtX7XGGmukmnfeeSfldEzRz1Evce1jRr0X6RpovyyaC9RPSvtSUr857fVIvTR1DkfkOUT99GgO/TPHz7Rp01JOPyf1NaN7fMABBxSx9hSLiOjatWsR67WMiOjRo0fKaQ9D7TUZwX0xR44cWcS6HkVE7LjjjimnvdWoXyn1f9WennR9tY/aKquskmo++uijlNM+atprMoL7644fP76Iqf/tYYcdlnINQddE/z71h6S+ddpXW3suR0S8/vrrRax9QyN4fdC1burUqamGevDpeki99KiXnc4DomXLlim3zjrrFPHbb7+darR3n/Z9jOC+q9rvm9Y17bEaEbH//vsXsfo1GgP1qNRzBfVkX3/99VNO9zxaQ/SsRf1btZdxRO6NSj0y6Xyie5y6DyJ4DdFeiNqfMoIdJjqGaZ49+uijRUxnA1pD9DvTz1HvWd131Y/QGKinvK6b1NuX7vHuu+9exDoXI/Jnp16X1Gda9xfqr04/p3OWnhVqOUfpuTWC1zvtWX3rrbemGnUykMeB7rGeF2jtpnV57NixRUxzvyrTp09PuS233LKIe/XqlWr0LBKR5yP1WG3evHkRa8/nCH6uU28AeQvoLKv3oZb1NiKPRbrHdFYfN25cg39v1KhRRUxrKa1bq622WhGTM4D2QvXb0JpRFbrHTzzxRBHr547gM4rO45kzZ6YaPd/SfKHe+/p8RH3nab3TcwT1BH7kkUdSTsdw+/btUw25Xa666qoiPuaYY1KN+p+oDzNdc3Uu0Bij9UDHD/WdrgrdY+2r/+6776Yaulc6H8kdU8u7H1rbdL+kfV49MRF53SL30hJLLJFy+lyse0AEPzNdcsklRaxnnYiIbbfdtoh1TY5gP4eeFe+6665UQ/4B7YtO79aq4H9hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXfB/Trq40korpdwZZ5yRctTEX1HBYkREixYtivinFIX9J0EyOBXxaLP8CJaPqEyLpBf6cyT+IJFUly5dipiEDyQyaNu2bRE///zzqWbllVdOORU3kCiC5EBPPfVUEe+0004N/hyJeUgKodIbFRBGsPBKJXYkhajKcsstl3IqiqDxQ2IBnaMq4YnI6wHJ0b766quUU0kT3ZdaZI00DujvTZw4sYh32GGHVEOSlvfee6+IVaISkSWdJOgkCc6SSy5ZxCR/IlmjCjQ+/fTTVFMVEnio5JGuAclrdF6R1EzXDZKUksBjv/32K2ISO5GITMf+XnvtlWro/unnInkg7Xt633Uc0mfq2LFjqlHpbETEK6+8UsQkZtx1111TTteIfv36pZoqkKRPRbRz5sxJNSQrUoEZjUtdC0i6q3MsIq9ZJH955plnUq5Vq1ZFrALPCJZy6TpKgjFax1TGROuazh/ah0k4qGOO5k/fvn1TTvcS+i5VIfHbQw89VMQkATrppJNSTgU/NDdUmqlSu4gsb4zI6/3HH3+cakhMq2KgrbbaKtWQBFalRiQ4o/3l9NNPL2I64+taQHsQjQ0VT+q5LiLi22+/TTnd96mmKiQuVRkSSblobdF7RTIxldfefPPNqUbFsBFZttmmTZtUQ/ugzj1djyIibrzxxpRTsTFJz2j/VuHzcccdl2pUmkd7gF6niPz9SLZHAq7WrVsXMZ3ZqkLPDzo/SIZJe6yOa1q3apH4qvQxIs9jWg9oH1BBHYnWSdysc0GFyBEsNFPR+eDBg1ONrgf6MxE893T9UZlZBK9benYmWWxVZsyYkXL6XLXZZpulGjpr6PsSElrrswnNFxK/6bM0nQ9oHdE5S2c5Wkd0DNOZjGR3ep4jCaLKcQla8/UZnM78JBfV9YAkt1Whd2Ivv/xyEdOZgd5X6Lseem+nY2Po0KGpht7FqGxTx2FExKGHHppy9913XxGrNDSCRcL6HmmLLbZINSR+1HWSBJn69/R6R/DZSp/LaRyed955KTdw4MAG/x49kzaE/4W1McYYY4wxxhhjjDHGmLrAL6yNMcYYY4wxxhhjjDHG1AV+YW2MMcYYY4wxxhhjjDGmLvg/18P6zjvvTDnqiVQL++yzT8pRjyCTadKkScppL7mDDz441VAfo2222aaIqUea9u+jnoKdOnVKOe1HRL+bejq//fbbRdynT59UQz2n9Oeo17b25orIfduoD6j2yqKeqtSfVvvZUY9e6lmm/de++OKLVFMV7TdHv1+vSQT3TtceUNozNyL329XeYBHcF1P7UtG4p2u+/PLLN/iZ6O8tWLCgiKl/FvWC1j6D9Lu1bzetm3Pnzk25xRdfvIipFxj1ENO+wHfffXeqqQr1QNd+wtTzt5Yey9Qnfd111y1i6kNM/eO1J+3MmTNTjd7ziNzr7P3330811BNy2rRpRUzjnPZQnUN07XSdpv6Ts2fPTrkjjzyyiLVfawT3Bdb5T3OoCtTPT+/n8OHDUw2tPdqrj8aO3l8aOw8//HDKaf/xc845J9XQ31OXhK5FP/Rz2iuUzkLkFdB+59rTOiKic+fORUx7Zy3rGvWQ1T0+Io/VWbNmpZqq0O864ogjinjHHXdMNdQnWPuJ0v6ivWep3yf1N1b/AfWGprVcfR20rtH6q/ODerrSvr///vsXsfZAjsh9Oql3Mq1Zeq30fBbBPWT1/EW926tC/S+1Nyr5F2jduu2224qYzrfqf6H+reT00L/3+OOPp5qePXumnI4D6q2p/aoj8r2hvshvvvlmyml/+gsvvDDV6P5Ja+Kvf/3rlNM9iK6BrrcRuecn9Xyn3vu1QD+nZ7mxY8emGurJq+cT6huuz3V0tqRzhj6jkeeE5rGORXpmpHVS9xR6PqJ9R38/XSd9hqLxe/HFF6ecegNontF6oF4n8h9UhdYyfSalff6oo45Kueuuu66IaY3ScyqdPWgP17FJvZrpGVF7AOv4jWC3i64jdL6m3um6llDPbPWYnH/++amG+qLrvk6OC3oXoteAnDQ012uB+pTrfKRnS9qbdK3WZ7iI3Heenompn7JeK3IM0ZlT1ynqx00MGjSoiNWfFsGeNe15Ts4Hfb7WZ6oIdtVMmDChiOkZhvYK9QLqs3xV/C+sjTHGGGOMMcYYY4wxxtQFfmFtjDHGGGOMMcYYY4wxpi7wC2tjjDHGGGOMMcYYY4wxdYFfWBtjjDHGGGOMMcYYY4ypC/7jpYu77LJLEe+8886Vfg81nddm56Z2VCIQkaU+d9xxR6rZc889U06FHSTPUQEgCRhIqKNCBJLDkWRDxSYkHCR5ogrNSEBFQoJapHnPPPNMEZO8a8MNN0y5d999t4hJrEdCJhUEtGnTJtVUhaSSKgSge0xjQ4U6JOZSeSIJV8eMGZNyuv4MGDAg1Rx++OEppyJGEuuRoG7rrbcuYppnNM5VBkIyCZWt0PVt0aJFyj377LNFTOIjErK89tprRfxTiRsisggyIs8rktnQ/H/yySeLWCVnEflekaxstdVWS7nNN9+8iEl0cumll6bc4MGDi1jXv4iI+fPnp5yud3/+859TDa0RKt5o165dqvnwww8b/Ey0Jo0fP76IaQ+n+aH7AgmLqvD111+nnApaVNoXwXI2FW6RoEXX0dGjR6canfcRWWSz7777phpaH3TfpXl3zTXXpJzeuylTpqQauna77757EZPUSUWbJAmkOaXiI5XBRET0798/5WiM/1Q89NBDKadzg6RgtGaplGuLLbZINS+99FKDv4f2PL0P+hl/6HcdeuihRTxixIhUQ/uLjnM6V5HwWec+CR21ZtKkSalmk002STnd455++ulUQ2cBvce0L1aFzoQqnxsyZEiqufrqq1NO52zv3r1TjcrgSIZOUvFtt922iGnfpz1I98+FFsr/rorOMHrGpj1BBVwReQzT31MJIokgac/T/ZPEgSRhVdElybaqQtIzFcbR8263bt1STgV4JNpV6SqtGSTp03WZJIinnHJKyumZqUOHDqmGhJy0Pyu63kbkdYrOTJojsSfJ9nRNUvlwBMvIVSJHYveq6PWNyPOB7hWdW1SYq8+oEVm2R3Ohe/fuKdeyZcsipvtCY1HP0/oZI3iPWWqppYqYroHOs4h8tiGp+PHHH1/EJHhV2XFEvlZ0rlAhaET+7CQbrYpep4j8fETjld79qACYBJl77LFHEdP5kmSNekahPe7+++9POd13SAKtn4n+Hknp6RroGYXmkD6PqOQ7gvdUlX+T5Jb2ikMOOeRHP2NV/C+sjTHGGGOMMcYYY4wxxtQFfmFtjDHGGGOMMcYYY4wxpi7wC2tjjDHGGGOMMcYYY4wxdYFfWBtjjDHGGGOMMcYYY4ypC/7jpYsqGnruuedq+rm//OUvP/p7IrJwwtSOClkictP5NdZYI9VQg3cVJOm9i8iyHpKDkDhQpTQkAlp22WVTTiV9JAsiCdg777xTxCTvIpmVilTod6uEjO4B/W4VBJBckMQCO+20UxGr3LAxkDxCBQwqkIuImDlzZsqp/I1EZHo/DzjggFQzaNCglFOp0HHHHZdqXn/99ZRbd911i5gEpM2bN085ve8k6CShmEq3VFgSkYUhJCci0ZGKhlZcccVUQyIpFdn9lJC05MgjjyxiugYkCzv99NOLmISgTz31VBGT5IykF5988kkRk3RRBYsReT7OmDEj1dBeqNItlb1ERHTq1Cnl+vXrV8T6fSPyd1ZJagQLK1UcQ9+FRC4qi6V5XQXaA5o1a1bEKjGN4H1CxZ4kJlOB0gYbbJBqSNZ4/vnnFzEJU+j+ah3JoGj91b1ZxWERvOfccMMNRUziN123zz333FRDa4gKjNZbb71Uo/cuIgtQdV9uDK1bt065pZdeuohpjJEcSaVRKqiKyOu07i0REU888UTKqZiwVatWqYaup94rEh3rWYQ+J+0TJJbSsU/XSSVZJBIlKZeK2EjMSPu3Xhc6M5H4qBZIVKhr6bRp01JN165dU06fYXTNjMjrDd2Xm266KeVUtERywR122CHl7r333iImES/tzSpUvPXWW1ONCl4jsiCT5l6XLl2KuBZhekSW6tJaSmc0nddNmzZNNVWhM77OvV/96lephp4pVEJNNSoK03N6BAvEb7vttiImqeVFF12Ucjo+6dlk1VVXTbn77ruviOn5nsaiiljpOadv375FrELxCJZh6vghcSrtJ/oZ6FxeFZ2fEfnZgCR9eo4h6Prq8zw93+vzWUSWDZPUnH7uggsuKOJXXnkl1eh+HZEFz3SdaH1VYeVmm22WanTfoz2Vnpf0mYbOViQu1Xm81lprpZqqkLx13LhxRUzr5CWXXJJyenajZy8VperfiuC1Rff1hx9+ONXQeVb3Bt07IvjZS88D9F1I8qh7A41zPZvrWI3gM5Gu1bSHk4hWrwudD3bbbbeUawj/C2tjjDHGGGOMMcYYY4wxdYFfWBtjjDHGGGOMMcYYY4ypC/zC2hhjjDHGGGOMMcYYY0xd8B/fw1p7a5r6QPuSRuR+OdRnmnpXDhs2rIipP6z21KKeaT179ky5r7/+uoi//fbbVEP9wcaOHVvE1OeHeh8us8wyRTx9+vRUQ729te8e9UvV/mBUo73IIiJGjRpVxNRP78MPP0w57bVI/S6rQj1F9T5QrzPK7bjjjkVMPe+0v7n25Yvg3rp6X2hMU7837XG15ZZbppr33nsv5bTHZi39LiMiVl999SKePHlyqtl0002LmHruUQ+6qVOn/ujvieBrrteO+p5WRftq0++nPpXUi177TE+YMCHV6H2huUd9zbTf+GWXXZZqTjrppJTT+6d9FiNy/8mIiDPOOKOItZdmBPf4vPvuu4tY+wdGROyzzz5FTP1E6R7r3KMevDSvdH7o36/KKqusknLaF1T74kfwHqC9CGltPfbYY4uYekHS39t+++2L+KCDDko1NJ7UPaB9tiN4n1Dfwssvv5xqaM3q3LlzEdN3OfDAA4uY9m/qL697DvWefPbZZ1NO+xxqL+PGQPNO1z/6nNRDUe/V7bffnmq0/yTdAzozaR9LWrPoXuneTGsBXQPti669PSO4d6eum7qfR0S8+uqrRUz9E2n9XWeddYqYrlOLFi1STq/d/vvvn2qqQmukrkn0OWnff/zxx4v4tNNOSzXaJ5T2t1qcNLoeRXDfee2rr+eHiNx3OiLPUVoPrr322pTTnp86ViIijj766CK+/vrrUw39Pe1PreMpIp+9IvL3ozH9U6K9km+55ZZUoy6biDxnaF5rb2+659Qz+7rrriticnzQmNbnSHpmI/+KrlOjR49ONeSN0fMQ9bnXMXzCCSekmhEjRqRcLVCvbT1r0NisSseOHVNu+PDhRXzppZemGvp++rvo2VZ7QW+99dapRtexiIiTTz65iGnvIPeHnp3JC7Hoooum3JVXXlnEQ4cOTTV03tL5QO4Pfb7u0aNHqqE+0+PHjy/io446KtUcfvjhKffkk08WMa23/fv3T7la6NOnT8pp33DaK26++eaU0+exXr16pRo9v9LeT2uL9g3Xs1ZExAMPPJBy+m6L7gt5U7S/+QcffJBq6DlZnR20X86dO7eIqc/1FltskXLPPPNMEdO7Cuq5rvsHneWq4H9hbYwxxhhjjDHGGGOMMaYu8AtrY4wxxhhjjDHGGGOMMXWBX1gbY4wxxhhjjDHGGGOMqQv8wtoYY4wxxhhjjDHGGGNMXfAfL1009ck333yTcioYI7nMkCFDUk7FJiTWW2qppYqYmuyTUEOlENrgPoIFgPq7XnjhhVTTtm3blFPxBzW0VxFkRBYpqcAoIksJqfE+fT8Vm+j1jmD5kspHVOTQGEh89sQTTxQxyeHo+6nsgGSYb775ZhGvuuqqqeaVV15JuQMOOKCIScxIohqVj9A9p8+g1+W3v/1tqlEpTUTEwgsvXMQkitDPSWP6iy++SLmvvvqqiEkkSqIaHcPt27dPNVUhYaRKzVQEEhFx4403ppxKt2jOqhBKZRYREb179045lT3tvffeqeb+++9PubvuuquISexE41xlQHSdSJ6jUo/BgwenGl3P77jjjlRD64gKTmltI/GQCvf0M1Zl0qRJKaciUZIHbbzxximncj29bxFZQkb38qOPPko5Xd8fffTRVLPiiiumnAqF6O+R/EV/ju7TBhtskHK6d9H6pH+PBEokKmzXrl0Rn3322TX9nK5jJGCtCgmbXnvttSKm9ZfGlMo+SdqpayvNMZLpaF0ta3tEltb913/lxwzKqchTz2wRPId0/SOJpoq6aA2rRUZM45f2fV2nL7/88lRzxBFHpFwt0DjQz/n++++nGrp/J554YhGTeFIFkioZjmChrZ5l6VxOc2/mzJlFTIJ2OsfpeYi+L8l5P/vssyI+9NBDU40K+FRSGBGxyy67pJzueXPmzGnw70fk+0ljsyokXdVz43HHHZdqSDCm+55K3iKykJP2KhLUqYBez6gRLE9UyVqrVq1SzVVXXZVy+jw2ZsyYVENj6oYbbihiOv/p+k4SVhJPvvjii0VM84XONSpd1eeQxkBzSNekK664ItXQZ582bVoRb7TRRqlGzzE0Vo488siU03Guz78RWcwYkfedTp06pRqSt7Zu3bqIVRQdwUJXlS7qe5CIfL6l8UvPgy1btixiOjc9+OCDKadrNT13VIXOd7o/qmg9gvdeFWLS+FEhJ52R6Bys0kN6z0P3SuXj9KxHa6k+f9LzIEkl77vvviKuRS6q58aI/LkjsqBSBboReYxFRNx6661FTM8CVfC/sDbGGGOMMcYYY4wxxhhTF/iFtTHGGGOMMcYYY4wxxpi6wC+sjTHGGGOMMcYYY4wxxtQFfmFtjDHGGGOMMcYYY4wxpi6wdNH8S6Am7CokIbGJCmAiItZdd90iJinMp59+WsQkPCQpw7LLLlvEKomIiFh++eVTTsUx9F1mz56dcipuUalABAugtthiiyJWAWFElkKQPIyEPipJWX/99VMNSRI+/PDDIv6pGu9HcPN/vQ8kdiJhiEqp9HNHROy6665F/Mgjj6Sa+fPnp5xKNkhgVov47PPPP081NDZUknDllVemGhJTqBjn22+/TTUqjCRpJ8khNde5c+dU8/LLL6ecCr3GjRuXag477LCUq4UOHTqknN4/FXVF8H1X+SSNA5W8ksxQxZ4REYcffngRL7fccqlmnXXWSbntttuuiEmU2rNnz5TTNWnQoEGphiQ0P/vZz4qY5tmoUaOKmOQrJI5RQa9KlSJY6NO9e/civu6661JNFXbeeeeU0/lCIiS6BypZ3HzzzVONCtToe5AoTNdkvUcRETvttFPKqTiLZLl073S+0H2idUzrSOKismVaQ0gao7I0GrskNNPzAs37qpAwTiWLJFKl76zi0nPPPTfV0L1S6FylZxY900SwHFLPTCRLev755xv8ORIUz5o1K+VUBkxCbRX2qoQoojaZNImGaU0eOXJkEdO5qip//vOfU073dDrLEnqNa5E70xmGBNdaR/JEGlMqo6T7QqJLncckBKXvd9FFFxXxpZdemmpUxEZnqN/97ncpt/vuuxfx4osvnmpo79LPrgLLxqAS34i85t99992phkS7em6kPUb3QpIS0r6nEtbzzz8/1ehzT0R+zqHnrObNm6ecyu+GDx+eamiu6/dTWWREvnb9+vVLNSTDXXvttYuY1mkS0Os6SWJ3FRLXCj1/6jMhrQf0jKjjhdZXFQDutddeNX0mPYNef/31qYbmgp7f6dmE1lcVcpJAm+S0eq1Ihq7CUXqWpmunz4gPPfRQqiG5n66T9LtJAFgLdP7p2LFjEeu6EpHX6Yh8r2hf0DMK7Z+Ezo8ePXqkGjoXdu3atYgnT57c4GeKyOOa1ls6F2od/T19jtR9KYJl1XqOoXMTSSV1Db733ntTTRX8L6yNMcYYY4wxxhhjjDHG1AV+YW2MMcYYY4wxxhhjjDGmLvALa2OMMcYYY4wxxhhjjDF1gXtYm38J2ms3IuKBBx4oYuq7td5666Wc9vh85513Us2WW25ZxNSfjPrKan9N7ZcdwT39/v73vxcx9e587733Uk77WVHPqY033jjltNcZ9YIdO3ZsEVMvsM8++yzltIcg9eGje6X9pagvVVWoz6n2gqbenXSvtB/ZDjvskGrefffdIm7atGmq0T7XEblflvamjsh9riNynz/qCbfCCiuk3NVXX13E1BuQ7rH2F6a+o3vvvXcRU68s7UMaka+nXssI7vGuvdNr7ctZC9TX8KSTTipi6vFO9/25554rYvou2iuY7jn13NQef9pn8YfQ/uLUW/L7779POR1n++67b6qhua5r7tChQ1NNixYtipj6uNFnGjZsWBF369Yt1VCfZe19SPO6CrT+qX+B+pHTur3hhhsWMc1p7Vd4wgknpJqPP/445bR/IY2vMWPGpFyzZs2KmNZtWo9ee+21IqbrTb0Q+/btW8TUd197flIPR+qZ3b9//yI+66yzUg311lSob3lV5s2bl3J6PqG5Qb1KtScm9fLUHpw0f6mvpI4pWnvI36F7rPb2jOBeiPfcc08RL7bYYqmG9n39/TTPdNxRP0rdbyLyGYbm/vbbb59yEydOLGLyPVSlVatWKaf3nTwuNGd1LNA+oeOAzuADBw5MOR0HtH/T79K+1nTmph7S2suTHCY0h7TXLPWU1vlI80x7b0dk18rFF1+carRPaEQ+a1177bWppiq0N2t/fFojqPe0PnvR/qV+B+pD//TTT6ec9rVt06ZNqqFnPe0B3KVLl1RDPay1j7d6DCIiTjzxxJTT9YZ6wer+Qc8B5FLQdZmur/Y3pp+j3r5VIXfOQguV//aRzsk0znUeDRgwINVMmTKliGkdob1fn4nVhRIR8be//a3B30Xjh36Xjmtap+le6f6sHquI3DudepnTs62Os6OOOirVqN8hIj8fUR/vqtAerk4LehdDrRtbrAAAIABJREFUfZf1nQadObWG3o3QM7H2Yab1gPp467Md7UPkHdI9Rj93BO+Fuv7Q2ea8884r4j59+qQaOqfp84lekwg+M+j+rA6IqvhfWBtjjDHGGGOMMcYYY4ypC/zC2hhjjDHGGGOMMcYYY0xd4BfWxhhjjDHGGGOMMcYYY+oCv7A2xhhjjDHGGGOMMcYYUxc0+UdEVttss80CbY5uzH/TpEmTFxYsWLDND/13jx/zQ3jsmMbg8WMag8ePaQweP6YxePyYxuDxYxqDx49pDB4/pjE0NH7+G/8La2OMMcYYY4wxxhhjjDF1gV9YG2OMMcYYY4wxxhhjjKkL/qGWIE2aNPkkIt7/530c829O0wULFqz8Q//R48f8CB47pjF4/JjG4PFjGoPHj2kMHj+mMXj8mMbg8WMag8ePaQw/On7+m3/ohbUxxhhjjDHGGGOMMcYY88/CLUGMMcYYY4wxxhhjjDHG1AV+YW2MMcYYY4wxxhhjjDGmLvALa2OMMcYYY4wxxhhjjDF1gV9YG2OMMcYYY4wxxhhjjKkL/MLaGGOMMcYYY4wxxhhjTF3w/9u78/Cqqqvx4zvMYSaBMIcxgIgQBhEQQVGcSmXUh8HH1gfBoaKiaNVSWnmsCiIWrQYUsVqFVhEVxSIgxSoIQiQo8xwmGSIJUwhT8vvj/fXtu/dack/OHXKSfD//rdV1z91tNvecu3uftTiwBgAAAAAAAAAEQrnCFNeuXbugadOmUVoKirv09PSsgoKCOj/3n7N/8HPYOwgH+wfhYP8gHOwfhIP9g3CwfxAO9g/Cwf5BOELtn/8o1IF106ZNzZo1a/yvCiVaXFxc5sX+c/YPfg57B+Fg/yAc7B+Eg/2DcLB/EA72D8LB/kE42D8IR6j98x+FOrB23sDvS1GCFBQU+Hod+wfGsH8QHvYPwuFn/7B3YAyfPeEoU0Z2I8zPzy+ClRQd9g/Cwf5BONg/CAf7B+Hws3/oYQ0AAAAAAAAACAQOrAEAAAAAAAAAgcCBNQAAAAAAAAAgEDiwBgAAAAAAAAAEgu+hiwAAAACiq2zZsiJ34cKFIlhJ+ErbgEUAAAD4wy+sAQAAAAAAAACBwIE1AAAAAAAAACAQOLAGAAAAAAAAAAQCPawBAABgEhISrPjo0aOipmvXriJXuXJlK960aZOoOXToUJirK72Ka79qAAAAwC9+YQ0AAAAAAAAACAQOrAEAAAAAAAAAgcCBNQAAAAAAAAAgEDiwBgAAAAAAAAAEAkMXAQAAYC699FIrnjBhgqjJzc0Vuf79+0dtTQAAAABKH35hDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQOLAGAAAAAAAAAAQCQxeLsSuuuELkVq1aVQQrAQAUJ6NGjRK5ESNGWPGJEydETXx8vMh98sknVjxt2rQwV4fCaNGihcidPn3airOyskTNyy+/LHJXXXVVyGuvXbtW5Nq2bWvFGzdu1BeLUm/w4MEiN3LkSCvesWOHqHn00UdFLi8vL3ILQ7HVsGFDKz516pSoSUpKErlf/vKXVpyTkyNqtm3bZsWrV68WNe7nLYLr9ttvFzn382f37t2iZsWKFSL3+uuvR2xd8CchIUHkjh49asV9+vQRNXXr1hW5xYsXW7H23ISSRTtLc597U1JSRE1iYqLIufeYOnXqiJqbbrpJ5M6ePRtynaUdv7AGAAAAAAAAAAQCB9YAAAAAAAAAgEDgwBoAAAAAAAAAEAj0sPaoXr16nuoOHjwYtTVcc801VvzKK6+ImnLl5J90w4YNVnznnXeKGq1vG4LpjjvuEDm3h5fWe3bZsmVWrPWIROz16NFD5M6cOWPF6enpsVoOSoBhw4ZZ8bhx40RN8+bNRa5ChQpWXLlyZVGTn58vcpdccokV33LLLaImLS3NiufOnStq4I/fz/Jz586JXI0aNazY3RPGGFOxYkWRo2c1NPfcc4/IuZ8FGu0ZJjs7W+QmTJjgb2EotgYNGiRyvXv3tuKhQ4eKGu2zLCMjw4rdXtiaH374QeS0vuyIPfe7+pQpU0SNO6vDq1//+tci5+4p7Xs5ImfTpk0i16pVK5ErU8bf7zHdXvTHjh0TNW5Pe2OMmTlzphW//fbbvt4f0bVnzx6Ra9SokcjFxcVF5P3Onz8vctrcDXdPzZo1S9RMmjQpImsqrviFNQAAAAAAAAAgEDiwBgAAAAAAAAAEAgfWAAAAAAAAAIBA4MAaAAAAAAAAABAIDF00xiQnJ4vcQw89ZMUXLlwQNd99953IzZkzJ3ILc7jDEt1BVz/n8OHDVnzq1KmIrQmRNX/+fCu+/PLLRU2VKlVELisry4q1v/GQIUOseOnSpaJm3rx5Ile1alUr1oYRuANvjDFm2rRpIgdJG6RQq1YtK27Xrl2sloNiZvTo0SLn7qmaNWtG7P20YTYNGjSw4jp16oialJQUK3aH+xljzBtvvBHm6lAY2r3E/VtqtMFk7vCps2fP+l8YSoxmzZr5ep028DWSn2MIJvf5UhsY7H4/M0beT9zB1cboA4Pr169vxfHx8aKmbt26VqwNb0Qw1K5d24pzc3NFzeLFi0Vu+/btVqwNEe7YsaPI/eY3v7Fi7d745JNP6ouFRTvTcL+nukM1vTpw4IDIHT9+POTrtOeY9u3bi9z06dOt+MorrxQ1d999d8j3g3fu2cSaNWtETevWraP2/tpARff7kfb5U716dZFzB4c+99xzokbb+2PHjg25zpKCX1gDAAAAAAAAAAKBA2sAAAAAAAAAQCBwYA0AAAAAAAAACAR6WBtjZs2aJXJu35uFCxeKmk2bNkVtTR06dBA5rSeSF40aNbLic+fO+boO/Bs+fLjITZ06VeTcXnkarW9S48aNrXj//v2ixu3N17RpU1Hz1FNPiVxSUpIVb9iwQdR46QUGYyZPnixyPXv2FLl9+/bFYjkoZu677z6Re/rpp0WuqHu9li9fXuTcz6jrrrtO1MyePVvkTp8+HbmFFaFq1aqJ3IkTJ4pgJf917bXX+nqd1h/W7R2q9YxE6dOlSxdfr9P+baxfvz7c5SDghg0bZsX9+vUTNdochZycnJDX1p6L3edp93PMGGPKli1rxdqMmEGDBomcNhMG0eV+d545c6ao+fbbbyP2fo8//rgV33vvvaJm8+bNVvz2229H7P1LEu1M5dNPP7XiNm3aiBrt++f48eOteO3atWGu7r8ee+wxkXPvc7fddpuo0fphjxkzJmLrKm1OnjxpxXv27BE1XnpYa/eOGTNmWPFrr70marTZUu6aWrZsKWpSU1NFbujQoVbszrEyxpj7779f5OhhDQAAAAAAAABAjHFgDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQSt3QxT/+8Y8ipw2FcQc9LV++XNRkZGREbF2uIUOGiFyzZs18XevRRx8Ndzm4iMTERJFzh6306tXL17Vzc3NFThsc4zb61wYWrVu3zopTUlJETfv27UXOHaCmDUHr37+/yEHq1KmTp7olS5ZEeSUoDpo0aWLF7mAOY/ThHNHkftYYY0zFihWtWBu66Grbtq3IaQO23n///UKsLriKesCipkqVKiFrCgoKRG7RokUid+jQoYisCcXb4MGDrbh79+6eXucOA58/f76oqVy5sv+FIXC6desmcgMHDrRibQhiuXLyq6s70Ez73EpOThY5ty4hIUHUxMXFWbF2f3MHMyL6Ro8eLXLu3zOSAxY17vPJqFGjRM2kSZOsmKGL3rlDxTMzM4toJf81efJkkXOHxfbs2VPUDB8+XOSys7OteMKECWGurvR48sknrTg9PV3UuN9N0tLSRM3f//53X++/c+fOkDVLly71dC33XMfd98bo9z33Hrpy5UpP71cc8QtrAAAAAAAAAEAgcGANAAAAAAAAAAgEDqwBAAAAAAAAAIHAgTUAAAAAAAAAIBBK/NDFyy67zIq1IQ3x8fEit337div+4YcfIrswhzv86N577xU17uCPvXv3iprPP/9c5ErK0KoguPrqq0Vu+vTpIte6deuQ17pw4YLIvffee1Y8fvx4UaM1+u/Ro4cVjx07VtS4QxZPnTolarR/C+6gIwYs+nf48GGRW79+vcg9/PDDsVgOAs7dB+79LBznz58Xue+//96K//GPf4gabQ/ffffdVqwNF3WHYuXl5Yma3r17ixz3r8jQBv96Gbp48OBBkXPvU8bo9zOUPs8884wVa88Umh07dlixOyTaGGM+++wz/wtDTDVo0MCK3aFkxhjTsWNHkevatasVV61aVdRUqFBB5I4fP27F+fn5okZ7dnYHjbvfszTaAN1oD/cr7R555BGR074nv/POO7FYzv9yP7eysrJEzaWXXmrF7mekMXKAHP5HEIYsenHfffdZcf369T297vrrr7dihi565z4P7Nu3T9Q88cQTsVpOWFJTU329zr3vlWT8whoAAAAAAAAAEAgcWAMAAAAAAAAAAoEDawAAAAAAAABAIJT4HtZz5861Yq2v0JkzZ0RuwYIFVpyenh7ZhTlefPFFK9b6trm++uorkRs1alTE1gRj2rRpY8VpaWmiplmzZiJXUFBgxatXrxY17t/cGGO++eYbK/bav2vgwIFWnJSUJGoSEhKs2O1NbYzemy8xMdGKR4wYIWreffddT+ss7dq1aydyWi/x7OzsWCwHATJo0CCRc3sfantF6+dZpoz9/0Vr97hp06aJ3B/+8IeQ69S467ruuutEjdtDuU6dOqLmqquu8vX+CK1v374iV7FixZCv+/Of/yxy2qwMLypVqiRybt9yrYfswoULfb0fYs99XnB71xuj9zt3ewB/9913ombr1q1hrg7RcOWVV4rcc889Z8VNmjQRNY0bN47YGtx+/Fo/05o1a4qcu1+90K7jpfc1/LvxxhtFrkWLFiKnPWPHkjbzoXPnzlZ8zTXXxGo5CJP2XVp7/vHbg/iDDz7w9ToYk5GRUdRLCMn9t2+MMa+++qqnOpf2LLVx40Z/CyuG+IU1AAAAAAAAACAQOLAGAAAAAAAAAAQCB9YAAAAAAAAAgEDgwBoAAAAAAAAAEAglaujizJkzRa5Vq1YhX7dp0yaRcweGRNKUKVNE7rbbbrNibRhSTk6OFb/55puRXRgEd5hMSkqKqClbtqzI7dy504pvvvlmUfPTTz+FfH93UJkx+sCOa6+91oq1QRHx8fFWvHLlSlGza9eukGtiwKJ39957rxUnJyeLmhkzZsRqOcYYY8aMGSNys2fPtmIvexP+aZ8Zw4cPF7lq1apZsTbsqXz58iJ3+vRpK3aHDxtjzKxZs0Ku06v3338/5JrcAUnacCRtgCQio3v37r5e5/fzXvus69evn8i5+0K7B7n3Lnd/o2i4g5yNkZ9t2r9p9/nIGDnwOTc3N8zVIRq0e9e4ceNErmfPnrFYzv9yhx7WrVtX1Ghrz8rKsuLatWuHfC/t+frkyZMhXwf/mjdvLnLHjx8XOW2YfSx99NFHIud+Zzt//nysloNCmjRpkhWPHj1a1GjP4V5ow/aef/55X9dC0atXr57IpaenW3GDBg0i9n6LFi2K2LWKI35hDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQOLAGAAAAAAAAAARCsRm6eOONN1rxU089JWq6du0a8jrHjh0TuQULFoic2yhdGwSUn58f8v3uuecekRs4cKDI1ahRw4q1oULPPvusFS9ZsiTk+yM8TZo08fW6gwcPWrE7PM0YYxo3bixyv/rVr6zYHaZojDGJiYki5w4/0oYa3X///VackZEhanbv3i1y8M8dFKMN/dmzZ0/E3s8dBqQNWNSG++Xl5Vnx66+/HrE1QdKG+mpD8WrVqmXF7vA5Y4wpKCgQub1791rx2rVrRc2RI0dCrtOrNm3aWHGlSpVEjTvcSqupUqWKyH3xxRdWrH0mIjR3gPDPcZ9r9u/f7+l1VatWteKpU6eKmr59+4pcuXL2Y+jWrVtFjTvQ7K9//aunNSG63M8nY+SzjntvMcaYffv2iZy7fw4dOhTm6hAN/fv3F7kuXbpE7Prnzp2zYu27kHbPW7hwoRUfPXpU1LifNcYYc9VVV1mxl6GL2lDq7OzskK8radzvR5mZmRG7tvs5og0v077DLF++PGJr8GPDhg0il5aWZsU//vhjrJaD/y81NVXkZs6cKXKdO3f2dX33O/+wYcNEzbJly3xdG8E0ffp0kYvUkEVtMGvTpk1Fzv0epT1vlRT8whoAAAAAAAAAEAgcWAMAAAAAAAAAAoEDawAAAAAAAABAIASyh/VDDz0kchMmTLBirXeexu11tmjRIlGj9V8bMmSIFWu9fVevXi1y9evXt+Jf/OIXoqZevXoi5/ar0fpgaX2JEV1uH0WtZ5rW86p9+/ZW/P7774sarYe12z9P68N34sQJkXP7sE+bNk3UfPXVVyKH6IqLi7Piw4cPixrtb+zSPjMGDRokcgMGDLBirWes5sCBA57q4I/bJ1Prw1y3bl2R03qeu3Jzc0XO7aen3Tu01/nl7s9GjRqJmh49elhx5cqVRY3239ftxw1/tJ7hGq2HtBdPPvmkFWvPPtoa3Gef8uXLi5obbrjBit2e1sYYM3fuXE/rROR06tRJ5Nzn4uTkZFFz/PhxkXOf1UtyL8biTOsD3bBhQ1/X+vjjj0UuPT3div/2t7+JGr+zVtxZL8YYc+edd1rxlClTQl7H7bdujDEXLlzwtabiLJI9q13us+uZM2dEzSOPPBK19/fr22+/FTn3/lUa90pRu+OOO0SuY8eOvq6l9dV/7bXXrJh+1SXfAw88IHLus2nLli1FjTY/qHr16lZcoUIFUaOdGb3zzjtWPHToUFGj9cMujviFNQAAAAAAAAAgEDiwBgAAAAAAAAAEAgfWAAAAAAAAAIBA4MAaAAAAAAAAABAIgRi62LVrVyt+8MEHRY3XIYuuzZs3W7E2AOaKK64IeR2t2flll10mcjVq1LDiKlWqhLy2McacO3fOirUhaDk5OVbsDng0xpirr75a5ObMmeNpDZDcATDa/tGGLrpDWbTBnhp3GMfKlStFTVpamsj985//9HR9xNaWLVus+OzZs6JmzJgxIjds2DArjo+PFzUVK1YUOW04oys/P1/k9u/fH/J18EYbltGtWzcr1gZxeBmwqA0iW79+vch98sknVhzJgavawCl3L/bu3VvUuEMlteF6Gi97GpI76NOrhQsXWrF2z9MGu7jDy7TPGW2Apvtcow2Iad26tRVff/31ouazzz4TuUgOFoWk3Tc2btxoxa1atRI12jOTOxhW+wx57733CrtEhMkdBK4NU9X+rbtDNN3PFWOMGTt2rMjt2rWrsEv0TBtwvW3btkJfx+v3Ovh31113WfG///1vUbNixYpYLScsM2bMsGLtPAHRpT17fPrppyKXlJRkxe751M/RnpNQsu3Zs0fkbr/9dl/XcvedNtB63LhxIteuXTsr1p6Dtefl4ohfWAMAAAAAAAAAAoEDawAAAAAAAABAIHBgDQAAAAAAAAAIBA6sAQAAAAAAAACBEIihi40aNbLimjVrhnyNNkxnzZo1IpeVlWXF7nBDY/QBCO5AqH379okabfCGO1jKK3eomjYI0l1DkyZNRM2pU6dEzh0EeezYMT9LLJXcgTOXX355xK6tDeBzh4h8//33ooYBi8XHwIEDrbhhw4aiRhti5/5bX758uajRhtF16NDBit39+3Pvl5GRIXLwJjEx0YoHDx4sam699VYr9jpw0LVu3TqR0+5727dvt+KTJ0/6ej93oIcxxrRp00bkBg0aZMXa56Q2jNILd1hO5cqVRQ3D9aQBAwaErNGGpVWrVs2Kr7vuOlHjDlg0Rv47WLt2rajRBh998803VnzTTTeJmtTUVCvW/h2wB2LPHRJtjBz8nZCQIGrcv7kx8u83evRoUcPQxdjr06ePFbdt21bUuAMWjTEmOzvbiidMmCBqojlg0Svte1QomZmZUVhJ6eU+HxljTN++fa34L3/5S6yWE3HuuUN6enoRraT0evHFF0XurbfeEjl3WHWXLl1Ezc033yxy7nOLO/jcGGOeffZZkSsug0MRXYcPH7ZibUjxli1bRG7OnDlW7N6vjTFmypQpVqwNbywO+IU1AAAAAAAAACAQOLAGAAAAAAAAAAQCB9YAAAAAAAAAgEAIRA/rSpUqWfGJEydEjdvj+fe//72omTp1asTW5PbR1no2aj32ypYtG/La2n+/vLw8K9b6TLt9Qb/88ktRQ4+/yOrRo4cV9+rVK2LXPnDggMi5f/dNmzZF7P0Qe6+88ooVv/vuu6JGy/nl9q968803Rc22bdsi9n6QMwK6du0qapKSkqxY6/mpce8VX3/9taiZO3euyK1cuTLktcuUkf9/tbv2Vq1aiZr27duLnHtvcvsge6XNmHD7cdOr2ButP7Vr9+7dIuf2V//d734napo1ayZy7ufKpEmTRM28efNErnPnzlas3WPdeSGXXnqpqOnWrZvIefl3AP+Sk5NFzu0x7z7bGmNMgwYNRM79t7948eIwV4dI6NSpkxVrMwzKlZNfJXNycqxY62kfBO5/Py+2bt0ahZWUXjfccIPIuXvqyJEjsVpOxLlzTWbOnFlEK8H/dfToUZH7+OOPLxobY8zLL78c8nX9+vUTNW5fdmOMeeKJJ6xY67UNGKPPfHB7X2tnkO6+c+eMGGPMsmXLwlpbLPALawAAAAAAAABAIHBgDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQAjF0cfbs2VackpIiaubMmWPF0R564Q4MGTlypKjRho+4Lly4IHKbN28WOXcwjTao5sMPP7RiBixG1vDhw0VuyJAhVuxlqKYmKytL5LShWO4wqYULF/p6PwTDggULYvp+7gAGd9gfIm/nzp1WrA0Ui4uLs2LtvqB9RmzZssWKP/vsM1Gzbt26kGts2LChyF1yySUi17FjRytu3bq1qElNTRW5evXqhVyDF+5/X2P0oTcIzR1+pw171oZjXn/99VasfYZkZ2eL3Pr16614+fLlntbZokULK9aGlrqDubVrHzp0yNP7IXK0Zyb3s8C9JxljTPXq1UXO/RwtDkOASoNGjRpZcdWqVT29zv0eN2LECFETyYHTfmkDXEN57LHHorCS0mvAgAEha5o3bx6DlYRPu19qw7JRfGn3tO7du1uxdkalnW1NnDjRis+cOSNqXn311cIuEaVErVq1Qta430lr1qwZreVEFb+wBgAAAAAAAAAEAgfWAAAAAAAAAIBA4MAaAAAAAAAAABAIHFgDAAAAAAAAAAIhEEMXXU8//bTIaUOqoskdluZnMIcxxqxevVrkZs2aJXJHjhyx4l27dokaL8O14M2VV14pcpMnTxa5OnXqhLzW6dOnRe6jjz6yYm343htvvCFyBQUFVux336F0cgetucPKjDEmOTk5VssplcqXLy9yx44ds2JtcJV2j3MHsLhDEY0xpmLFiiJ3/vx5K77ssstEjTYcsmfPnlbcoUMHUVOhQgWRK1Om8P/ft3aPW7JkichpQ4oR2r/+9S8rvu+++0SNNnTRdfbsWZHTBji6nzXuECJjjGnWrJnIjR8/3ooTEhJEjTsEW3uu0vYTIuf2228XOW0Aq8vrkD6GigfT9u3bI3Kd6dOni5w2YGzu3LkReT+NO+DVGGOaNm0a8nXbtm2zYu5JkZWYmBiypl27djFYSfi0e9yqVauKYCUoSoMHDxa5L774QuRq165txX369BE1DF2EV7m5uSLn7rHWrVvHajkRxS+sAQAAAAAAAACBwIE1AAAAAAAAACAQOLAGAAAAAAAAAARCIHtYx7pf9Z133ilyffv2Dfm6c+fOidybb75pxZ9//rmomT9/vsi5fU+1vsjwr27dulY8ceJEUdOwYcOQ11m6dKnIjRkzRuQ2btxoxVrP7KysLJGrV6+eFReXvm0Ihp07d4as0T634F+vXr2sWOvxXKVKFSuuWbOmqHH71xtjTPv27a1Y6/Xo9sc2xpgaNWpYcZMmTURN2bJlQ+a0Huh+nTx50oq///57UbNixYqIvV9pt3jxYivOzMwUNdq+cGn7Wdu/rVq1suIXXnhB1Gj9PePi4kKuwb2far0gEV1bt24VuUOHDolc/fr1rVh7nv/ggw9E7k9/+lMYq0O0uM8UP/30k6jx0oNY62XepUsXkfvmm2+seP/+/SGvrbn11ltFbvTo0SLnZU7NSy+95GsN8ObUqVMi5z4zefk7RZs2C8SdK/L222/HajkIsB9++EHktN733bp1s2Jttoz7nKY9yyH2hgwZInLu98E1a9aImkh+RrhnhdpMI7evtfb8XhzwC2sAAAAAAAAAQCBwYA0AAAAAAAAACAQOrAEAAAAAAAAAgcCBNQAAAAAAAAAgEAI5dDGa+vTpI3IPP/ywyLlDEDVbtmwRuXnz5lmxNnRRc/78eU918OeWW26x4t69e4ua/Px8kdu3b58Va4Ncjh49GvL9u3fvLnJehjzu2LEjZA3wH5UrVw5ZE+uhtiXd9u3brfjAgQOixh02V6aMt/+vuHbt2heNg8IdGLlnzx5R8+WXX1rxsmXLRM2CBQsiuq7SzB3GOWrUKFEzYcIEkevZs6cVa58X8fHxIpecnGzF2lBPbcCie31tD0yaNMmKtWF/iCx3mFjjxo1FjfY3dof+aENhtc8HBNOSJUus+ODBg6LGy9BFzc033yxy7hCp2bNni5qkpCSRGzBggBX3799f1KSkpIRc06pVq0Tu3XffDfk6+JeRkSFy7qB67XkhmiZOnChydevWFTntHADo27evyKWmpoqce9Z09uxZUaPdQxFbTZs2Fblhw4aJXPXq1a3Y/W5kjDEtWrQQuUWLFlnx3r17RY17jmWMvKdpz2TuuYB2b3ziiSdELmj4hTUAAAAAAAAAIBA4sAYAAADaB13bAAAIOklEQVQAAAAABAIH1gAAAAAAAACAQIhqD2ut31Osew+6vUMfffRRUdOuXbuQ18nKyhI5t6+iMd57VrsqVKhgxVofI/jXsmVLK9b6/Gi9hvLy8qy4atWqokbrYe3us8cff9zTOt01fPfdd55eBxhjzIkTJ0LWeOnPD+/cntXr1q0TNe7/5u7nvTF6X+s6deqEubrw5OTkiNzx48dFLjMz04oXLlwoap555pnILQyFtnjxYk859941ZMgQUZObmytymzZtsuJy5eTjZUJCgsj9+OOPVvzyyy+Lmq1bt4ocouvpp5+2YvffuDHG1KhRI+R1KlWqJHLubBAEl3t/mzFjhqjRemtee+21VqzNiNG+I1599dVW3LlzZ1FTq1YtkevRo4cVe50T4c4PcucQGWNMdna2p2vBH61HuNvD2v0OF2nuPId77rlH1Gh9Ze++++6orQn+uX17tWeWSHJnob3wwgueXueuS+uhrz2HI7Zq1qwpcto9pnnz5lbsnkEaY8yRI0dEzv2O2KRJE1HTunVrkfMyC81VXGeI8AtrAAAAAAAAAEAgcGANAAAAAAAAAAgEDqwBAAAAAAAAAIHAgTUAAAAAAAAAIBCiOnTR74BFbTCPNtjOi/bt21uxNiRGGyLlWrNmjci98847vtakYchidO3YscOKtSEq2iCXxMREK05LSxM1p0+fFrnBgwcXdonGGGM2bNhgxa+++qqv66B06tixY8gabeBoSkqKFW/bti1iaypt3nvvPZFzB1dpA2DcITHGyKFUjRo1EjXa3/PUqVNWrN33zp07J3InT5604o0bN4qa9PR0kXMH7n399deiBsXD888/f9HYK21IjTaA7+DBg76uj+hq0KCBFd90002iJi4uLuR1qlevLnLa3kDx8PHHH4tcfHy8yLn3vN69e4sabfhwmzZtrLhevXqFXWKhTJ061YonT54c1feD9NVXX4mcO4y3V69eoubDDz8UOXeQsDao7IorrvCUcw0YMCBkDWKvbdu2Ivfcc89Z8eHDh0XNypUrRS4pKcmKO3XqJGrcgbLG+L+nzZ4924p/+9vf+roOoisjI0Pk3HuHMcbcddddVtyvXz9R06JFC5Fzv4O7Z0/h2LVrlxUPHTo0YteOJX5hDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQOLAGAAAAAAAAAARCVIcu+uVlwKI2iEMb3nPkyBEr1oaJucMdjDEmJyfHipcsWRJyTQiu7du3W/GWLVtETbdu3UTOHRKlDVuoWLGirzWtXr1a5EaMGOHrWoAxctCRVx06dLBihi76t3PnzpC5zz//XNR4+fxp1qyZqNEGV7kDZN1BMsbIQYnGGDNv3jwr1vaTe28ENOyT4s0dSJWamurrOnl5eSI3d+5cX9fSVKtWzYpPnDgRsWtD2rNnj8i99dZbIteyZUsrPn/+vKjRhgG7+0wbDlyunPzqeuzYMSvWvg+++OKLIvfaa6+JHGJr/fr1IucOfZ05c6ao6dGjh8g1btzYirU9pu1F96xg8ODBombFihUih6Ln/s2NMaagoMCK77jjDlEzcuTIqK1J8+CDD4rcSy+9FNM1IHK0YbHudy9tGKf2LOW+Lj8/X9Roz1KHDh2yYm0o8tixY0WuOOIX1gAAAAAAAACAQODAGgAAAAAAAAAQCBxYAwAAAAAAAAACIZA9rL3Q+rto3D58mZmZomb//v0RWROCa+nSpVa8ceNGUZOVlSVyHTt2tOKyZcuKGq2fnnv9jz76SNRo/fSAcGzevNmKtX6TGrdXMqLL7TtmjN57DACKwvjx461Y6yVcpoz8zcuzzz5rxdF+vqZnddHT7mdubvny5Z6u1bVrVytu06aNqKlatarIud/tvv76a1Hj9rlG8XHXXXd5qmvVqpUVa58P2n71eqaA4NG+z7tzeLp37y5qtM+R+Ph4K969e7eo0fZPWlqaFWt9/VHyzZ8//6KxMfpzkzvDKC4uTtS4ffaNkfvc7d1ekvALawAAAAAAAABAIHBgDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQiu3QxcOHD3uqcwcpMGARxhgzcuTIqF5/0KBBVjxv3ryovh+gue2220Ru1apVRbASAEBx8dNPP1nxAw88UEQrQWny7bffXjQGLmbr1q1FvQTE2N69e0Vu3LhxF42BoqINeF2xYkURrKR44RfWAAAAAAAAAIBA4MAaAAAAAAAAABAIHFgDAAAAAAAAAAKBA2sAAAAAAAAAQCAU26GLQJAxZBFBwIBFAAAQSqVKlUQuLy+vCFYCAADwP/iFNQAAAAAAAAAgEDiwBgAAAAAAAAAEAgfWAAAAAAAAAIBAoIc1AAAAAJRS9KsGCqdGjRpWfOzYsSJaCUqChg0bitz+/fuLYCVAsPALawAAAAAAAABAIHBgDQAAAAAAAAAIBA6sAQAAAAAAAACBwIE1AAAAAAAAACAQfA9dLCgoiOQ6UMqwfxAO9g/Cwf6BX+wdhIP9g3CwfxAO9g/Cwf5BONg/8ItfWAMAAAAAAAAAAoEDawAAAAAAAABAIMQV5uf5cXFxR4wxmdFbDoq5JgUFBXV+7j9k/+Ai2DsIB/sH4WD/IBzsH4SD/YNwsH8QDvYPwsH+QTguun/+o1AH1gAAAAAAAAAARAstQQAAAAAAAAAAgcCBNQAAAAAAAAAgEDiwBgAAAAAAAAAEAgfWAAAAAAAAAIBA4MAaAAAAAAAAABAIHFgDAAAAAAAAAAKBA2sAAAAAAAAAQCBwYA0AAAAAAAAACAQOrAEAAAAAAAAAgfD/APaaMgs+MuIUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c03f4e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "K = 20\n",
    "Sigmas = [0.3*i for i in range(10)] # noise variances\n",
    "noisy_imgs = []\n",
    "rec_imgs = []\n",
    "for (img,label),sigma in zip(test_set[:len(Sigmas)],Sigmas):\n",
    "    img = img[0][0]\n",
    "    noise = np.sqrt(sigma)*torch.norm(img)/np.sqrt(28*28)*torch.randn(28, 28)\n",
    "    noisy_img = img + noise\n",
    "    rec_img = denoise(decoder20,Variable(noisy_img),K,0.1)\n",
    "    rec_img.data.clamp_(0, 1)\n",
    "    noisy_imgs += [noisy_img.numpy()]\n",
    "    rec_imgs += [rec_img.data.numpy()]\n",
    "\n",
    "# plot and save to file\n",
    "fig = plot_images(noisy_imgs,rec_imgs)\n",
    "fig.savefig(\"denoising_ex_0.1.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate plots for paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 782\n",
    "def estimate_MSE(decoder,dataset,Sigmas,K,numit=1):\n",
    "    MSEs = np.zeros(len(Sigmas))\n",
    "    for i,sigma in enumerate(Sigmas):\n",
    "        for (img,label) in test_twos[:numit]:\n",
    "            img = img[0]\n",
    "            noise = np.sqrt(sigma)*torch.norm(img)/np.sqrt(N)*torch.randn(1, 28, 28)\n",
    "            noisy_img = img + noise\n",
    "            rec_img = denoise(decoder,Variable(noisy_img),K)\n",
    "            rec_img.data.clamp_(0, 1)\n",
    "            MSErec = torch.norm( rec_img.data.numpy() - img[0] )**2 / torch.norm(img)**2\n",
    "            MSEs[i] += MSErec/numit\n",
    "    return MSEs\n",
    "\n",
    "numit=100\n",
    "Sigmas = [0.3*i for i in range(22)] # noise variances\n",
    "\n",
    "K = 10\n",
    "MSEs10 = estimate_MSE(decoder10,test_set,Sigmas,K,numit)\n",
    "\n",
    "K = 20\n",
    "MSEs20 = estimate_MSE(decoder20,test_set,Sigmas,K,numit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Slope of MSE10: 0.03319207866007233\n",
      "Slope of MSE20: 0.05402469204793347\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:2: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \n",
      "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  import sys\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4VEXbwOHfJNn0XkkhECB0UCACr0pTEITQFGxURbGABTu2z1eliiIKCIgoVhQLLyTSm0oPSA0tJCQESO89uzvfHydggJRNSM/c15XL7J6Zs3MAn3N2yjNCSomiKIrSOJjVdgMURVGUmqOCvqIoSiOigr6iKEojooK+oihKI6KCvqIoSiOigr6iKEojooK+oihKI2JS0BdCDBJCnBZCRAghXi/h+EQhRKIQ4nDRz+PFjk0QQpwt+plQlY1XFEVRKkaUtzhLCGEOnAEGALHAAeBhKWV4sTITgSAp5dTr6roCYUAQIIGDQDcpZWoVXoOiKIpiIgsTynQHIqSUkQBCiFXAcCC8zFqagcBmKWVKUd3NwCDgx9IquLu7y+bNm5twakVRFOWKgwcPJkkpPcorZ0rQ9wUuFHsdC/Qoodz9QojeaN8KpkkpL5RS17esD2vevDlhYWEmNEtRFEW5QggRbUo5U/r0RQnvXd8ntA5oLqXsDGwBVlagLkKIyUKIMCFEWGJioglNUhRFUSrDlKAfCzQt9toPuFS8gJQyWUqZX/TyC6CbqXWL6i+TUgZJKYM8PMr9dqIoiqJUkilB/wAQKIQIEEJYAg8Ba4sXEEJ4F3s5DDhZ9PtG4B4hhIsQwgW4p+g9RVEUpRaU26cvpdQLIaaiBWtzYIWU8oQQ4j0gTEq5FnhOCDEM0AMpwMSiuilCiPfRbhwA710Z1FUURVFqXrlTNmtaUFCQVAO5iqIoFSOEOCilDCqvnFqRqyiK0oiooK8oitKIqKCvKIpynYPRqXy7N5rTcZkYjXWrC/xmmbI4S1EUpdHIKdDz5LcHScrSZqG72lnSvbkrPVq40rOFG228HDAzK2kJUuUYjAb2xe0jNDIUnZmOd29/t8rOXRIV9BVFUYpZuTuapKx8Fj7ShdwCA3sjU9gXlcyGE3EAONvqim4CbvRs4Uq7Jo4VvglIKQlPCSfkXAgbzm8gKTcJe509Q1sOrY5LuoYK+oqiKEXScwtZsvMcfdt4ENzZB4DRQdr60tjUHPZFprA3Mpl9USlsCo8HwNHagu4B2g2gR4Ab7X0cMS/lJnAh8wJ/RP5BSGQI5zPOY2FmQW/f3gS3DKa3X2+szK2q/RpV0FcURSmy/K9I0nMLefmeNjcc83Oxxa+bLfd38wPgUlou+6KSr94ItpzUbgIOVha083akrbcDbZo44OtqJDp/D1ti1nM48TAA3by6MaHDBAY0G4CTlVPNXSAq6CuKogCQlJXPl39HMaSTNx19yw/EPs42jOzix8gu2k0gLj2PfVHJHDifwonLyfx2eieGCwcxtz+DEEbMCr1pqhtFT8+76ebTkrYuDtjr7Kr7sm6ggr6iKAqwePs58goNTBvQulL1PRx0eHhGI7NCuZS9BWGbg5e1B52dR+Eme5KQ4sapyxl8fSqdFfIQAFYWZrT2cqBtE+1bQSdfJ3q0cKvKy7qBCvqKojR6F9Ny+W5vNPd39aOVp73J9a4MyIZGhrI+aj1JuUk46BwYFDCI4BbBdPPqhpm4dmZ8XqGBiIQsTsVlcupyBqfiMtl+OpHVB2Pp4u/M78/cUdWXdw0V9BVFafQ+23oWieT5/oEmlb8yIBsaFUpUehQ6Mx29/XozpMWQcgdkrXXmdPR1uqELKSkxgTSjzU1dhylU0FcUpVGLTMxi9cFYxvVshp+LbanlUvNS2Xh+I6GRoVcHZIO8ghjffnyVDMi6b3sJ96wEeGwDiKpbB3A9FfQVRWnU5m85i6W5Gc/0a3nDsVx9Ljsu7CA0MpRdF3ehl3paObfi+a7PMyRgCN723iWcsRLSLsCpELj9uWoN+KCCvqIojVj4pQzWHbnEM31b4ulgDVy7QnZL9BZy9Dl42noyrv04hrQYQmuX1oiqDswHlmv/vW1S1Z63BCroK4rSaH206TSO1hZM7tWCE8knrhmQtdfZM7D5wKsDsuZm5tXTiIIcOLQS2g4BZ//q+YxiVNBXFKVROhidyrZzJ+nX7SLjNy8iKj3q6grZIS2G0KdpnxpZIcux1ZCbCj2eqv7PQgV9RVEamdS8VDZEbWD+3p+wb3WOA+naCtlx7cdxT7N7anaFrJSwfxl4dYRm1TtV8wqTgr4QYhCwAG27xOVSytmllBsFrAZuk1KGCSGao+2Xe7qoyF4pZc3czhRFUYqUNCBrKPCij+943u47Bh97n9ppWPQuiD8Owz6r9gHcK8oN+kIIc2ARMACIBQ4IIdZKKcOvK+cAPAfsu+4U56SUt1ZRexVFUUxS2oDs2PZj2bzfl4wMDz6Z1Bcri2rqqzfFviVg4wKdRtfYR5rypN8diJBSRgIIIVYBw4Hw68q9D8wFXq7SFiqKopiotJTFxQdkt5xM5LOYg8wd1bp2A35aDJwKhTueB131L8q6wpSg7wtcKPY6FuhRvIAQogvQVEoZIoS4PugHCCH+ATKAt6SUf91MgxVFaVzyCg3ojRJ7q9LDVWkpi68fkDUYJR9tOk0LDzvu6+JbU5dQsgPLAQFB1T9NszhTgn5JHU1X9w8TQpgB84GJJZS7DPhLKZOFEN2ANUKIDlLKjGs+QIjJwGQAf//qn7KkKEr9cORCGlN/PERGrp63hrRjVDe/q3Pkr6yQDYkM4UjiEUAbkB3fYXypA7Jrj1zkTHwWix7pioV5Le4WW5ADB69M02xaox9tStCPBYq3yg+4VOy1A9AR2FH0l9EEWCuEGCalDAPyAaSUB4UQ54DWQFjxD5BSLgOWAQQFBTWsDSkVpQHJyCvEYJC42FlW6+dIKflq13lmrT+Jp4M1rTzteeWXo6w5EsW93VPYm7DlhhWygwMGlzkgW6A3Mn/zWTr4OHJvxybV2v5yHfsZ8tJumKYppaz6hV/XMSXoHwAChRABwEXgIeCRKwellOmA+5XXQogdwMtFs3c8gBQppUEI0QIIBCKrsP2KotSAAr2Rb/acZ8HWs+TrjUzt14on+7Solj7x9JxCXvnlCJvC4+nfzos5ozpwKvUfPjuwiSNpuzh6KB97czfGtBvL0JbBJq+Q/TnsAjEpOXz16G1VusdthUkJ+5aCVydodjsAhvR0EhcuQhYU4P3fd6v148sN+lJKvRBiKrARbcrmCinlCSHEe0CYlHJtGdV7A+8JIfSAAXhKSplSFQ1XFKVm7DidwHsh4UQmZtOntQf2VhZ8vPkMa/65yAcjOnJ7K/fyT2KiwxfSmPrDIeLSc3mivw4Lx62MCn3l6oDsvc0HEhXdhrDjLvyd7srQpj4mBfy8QgOfbj3Lbc1d6Nvao8raWynn/4aEcBi2EGk0krb6FxIXLMCQlobzQw9W+9O+kLJu9aYEBQXJsLCw8gsqilKtopKy+SAknK2nEmjuZsvbwe25q60nQgh2nE7gnf+dICYlhxG3+vDmkPZ4OFR+9aqUkhW7zjNny984uh/H1es4l3NiShyQlVKy9sgl/rsunMy8Qp7u24op/VqW+a1j2Z/nmPnHKX5+8j90D3CtdDurxKoxEL2bnN7fEjf3Y/JPnsQ2KAivN9/Aul27Sp9WCHFQShlUbjkV9BVFKS4zr5CF2yJYsSsKKwtznr2rFRPvaH5DUM0rNLBoewRLdp7DWmfOa4Pa8kh3/wp3nUSnJjB17decy9mJuW0MoA3IDmkxpMwVsinZBbwfEs7v/1yklac9c+7vTLdmLiVeT6+527nFz5mVj3WvUNuqXGo0hbO6Eh/bjcxDF7Dw9sbr1VdwGDTopp/uVdBXFKVCjEbJr4dimbPhNElZ+Yzq5serg9pczT5ZmoiELN5ec5w9kcnc2tSZGSM70sGn7FQGV1bI/nBiDYeT9oIw4qbzZ2zHkQxuUfaA7PW2n07grd+Pcyk9lwn/ac4rA9tgV2x65/zNZ1iw9Szrpt5JJ7+a3YS8OGNuLslvjCV50wmwsMLt8Sdwe3wSZjZVM0dfBX1FUUx2KCaV/649wZHYdLr4O/Pu0A7c0tTZ5PpSStYcvsgHISdJzSlg4u0BvHhP62vm1pe0QlYWOmKVH8Q7/cYxrF23Sj/tZuXrmbfxNCv3nMfHyYYZIzvSt40nKdkF9JqzjT5tPFg8plulzn2zpJRkbthA/Ny56C/H4dDRBa9Pf0HnU7WpH1TQVxSlXPEZecxZf4rf/rmIp4MVr9/blhG3+lZ6dkt6TiFzN57ih/0xeDlY805wO5r5pBIS+e8KWTsLO2z1XYmObstdzXoyb3QXnGx1pZ80Mw4Kc8C1RbmffzA6hdd+PUZEQhYju/hirTPjpwMX2DStN608HSp1TTcj79Qp4mfMJOfAAaz8PfFqFY7da/+7OmunKqmgryhKqfIKDazYFcXCbRHoDZJJvQKY0q9VmateK2LD6RO8v/17Us32Ym6ViIWwoLdfb9o79mXlFlsSMw1Mv7cdj97RvOyneylhaW9IiYInd4LbjbtbXS9fb2DRtggW7ziH3igZ1c2PeaNvqZLrMpU+NZXEBQtI+3k15o6OeDz/HM5JnyB0Onjyr2pJrmZq0FeplRWlkUnKymf0kj1EJWUzoL0Xbw5uR3N3u5s+7w0rZG3A37oDF2J6k5PZCXOrlny45RJNnHT88lQP07qPLuyDuKOAgNUTYNIW0JU9xmBlYc6L97RhcGdvvtsbzbN3mbbZeVWQhYWk/riKxIULMWZn4zJmDB5Tp2CecgRWnoLhi2osm2ZpVNBXlEbms61nry5S6tfG86bOVdYesldWyF5Oz+W9deH89s9FBnbwYu6oW3CyKaM7p7j9y8DKCYZ+Ar88Chte1343QdsmjnwwotNNXF3FZO/eTdzMmRREnMPu9v/gNX06VoFFN5z1S8HWDTqOqrH2lEYFfUVpRGKSc/hhfwwPBDWtdMCv6B6y3k42fD62GxfTcvFxsjZ9sDbjMoT/D7o/CR3vg8uHYdcCaH4ndKr94HlFwYULxM+eQ9bWreiaNsVv4WfY3333v9eZeh5O/wF3Tiv3W0pNUEFfURqR+VvOYCYEz99dsS6P0lIWDwoYxJCAISbtIevrXMGpiQe/BqPh383C73obLuyHtc9Bk87g0bpi56tixuxskpYuI+Wrr0Cnw2PaNFwnTsDM6rpFarWUTbM0KugrSiNx8nIGaw5fZHLvFjRxMu2J8/qUxTozHb18exHcMpjefr2rbw9ZfQEc/AoCB/w7eGuug1ErYMmdWv/+41vB0rZ6Pr8MUkoy1q0jYd5H6BMScBw2FM+XXkLn5XVj4YJsOPQNtBsKTrWcyrmICvqK0kjM23gaeysLnu5T9gyY0lIWT+gwgQHNBtTMHrIn10JWPHSffO37jj5w3zL4bhSsf0UbGK1BuceOEz9jBrmHD2PdsSO+Cz7BtkuX0isc/Qny0mts03NTqKCvKI1A2PkUtp5K4JWBbXC2vTEtcmkDsi90fYHBAYPxtveu2QbvX6bNy295943HWvWHXi/BX/O0zcRvfeTGMlVMn5REwvz5pP/2O+ZubnjPmIHTyBEIszJy8l/JptmkM/j3rPY2mkoFfUVp4KSUzNlwCg8HKx69o/nV9/VGPfsv7yc0quQB2TaubWqnwZcOa1M1B86C0oJq3+lamZAXwacLeFY+UVmJshKhIBNp70fKt9+RtHgxxoICXB99FPdnnsbc3r78c0TthMRTMHxxrU/TLE4FfUVp4HacTuTA+VTeH9ERG505J5JPXDMg66BzYFDAoKt7yJqJWtxRCmD/F6CzLfsJ3twC7l8OS3rBzxPgiW1gZUIgLo+UcPh7WP8aWecLiT/ehIKUQuz79Mbz9dexCggw/Vz7lhVN07z/5ttVhVTQV5QGzGjUnvL9PHLIsl7PsDWhVwdke/tpKYurdUC2onJS4Nhq6DIGbMpZvOXQBO7/Ar4ZAaEvwcglN/dEnZ0E654nf/964k81JzsyF0unQpr2Tsa+9Z8Q0wo8JoC9CVNdr0zT7PVSnZimWZwK+orSQKXmpfLh3z8RY7UOc9sYFh+BIK+gmh2QrahDK8GQD7c9YVr5Fn2h7+uwYxY0vwO6jq/c557ZhGH1FJLCCkg5442ZtRmer76K6yMPI85vgwNfwPYPYOcc6DBCa1/T7qXfZPZ/AcLs3+mmdYgK+orSgBQfkP374i4MUo+VlQ9TurzAkBa1MCBbEUYDHPgSmvcCr/am1+v9CsTsgT9eAZ+u0KSj6XULspEb3iL9lx9JOOaKIU+H030j8Jw2DQv3oh3B2g7WfpIitDn3h7/Xvo006Qzdn9BW2RafOpqfBYe+hfbDtdlGdYxJCdeEEIOABWjbJS6XUs4updwoYDVwW9Gm6AghpgOT0LZLfE5KubGsz1IJ1xSlYkpaIetl60UzqzvZftCPLx8ZdtPpFmrEqVBY9Qg88C20H1axulmJ2vx9K3uYvAOsTMioefEgOYsnEb89k7wUS2xu6YzXW29j06mcm0Z+lrax+f7lkHACrJ2hy1jtqd61hXbjCn0RHttYo7N2qizLphDCHDgDDABi0TZKf1hKGX5dOQcgFLAEphZtjN4e+BHoDvgAW4DWUkpDaZ+ngr6ilK+kFbIOOgfuaX4PQ1oMoZ3LLfT98E9aeNjx0+Se1brnapVZOQySz8HzR7SB2oo6/zesHAodRsL9X5be9WLQU7jufRKWrCTjvA0Wbk54vv4mjsHBFftzkhKid2tdPyfXad9UWvWH5AiwdoTJO2t01k5VZtnsDkRIKSOLTrwKGA6EX1fufWAu8HKx94YDq6SU+UCUECKi6Hx7TPhcRalRF1JymLBiP0/1ackDtzWt7eaUqKQVsiUNyC7aHkFSVj5Lx3WtHwE/8bQ2xfHudyoX8EHLyXPXW7D1PW3+fgn96cZLJ0l5ZyJJe9MAW9wmjcf9mWcxs6tEllEhtHGE5ndoeYIOfq2tIs6KhxE3OahcjUz50/UFLhR7HQv0KF5ACNEFaCqlDBFCvHxd3b3X1a0ba5EV5TqLtkcQmZTNa78dxdbKnODOdaM/tqQVskFeQUzsMJH+zfrfMCCbml3Akh3n6N/Oi27NankTcFPt/wLMLaHrhJs7zx3TIHqPlo3Ttxv43AqANBrJ+uIt4r/4hcIsc+xv64DXzAVYNq2im7ujN/Sbrs3WiTsGvl2r5rzVwJSgX9Lt6mqfkBDCDJgPTKxo3WLnmAxMBvD39zehSYpStS6k5PDLwVgeDGpKZFIW0346jL2VBX1rqS+8pBWygS6BJq2QXbLzHFkFel4ZWEuLqyoqLwOO/KjNZ7dzv7lzmZnByKWwtJeWn+fJP8k/c4b4154hOyoLK3cb/Ge9j92ACo4ZmMrCEvxqZ1tGU5kS9GOB4rdDP+BSsdcOQEdgR9HXyCbAWiHEMBPqAiClXAYsA61PvwLtV5QqsWh7BGZmgmkDWmNjac7Dy/by1HcH+XZSD25rXjNPyyWtkPWy9WJch3EMCTBthezl9Fy+3n2ekV18adOk5rcHrJQjP0JB1o15dirLzg1GrcCwdAiJk+4h9VA6ZjqJ15i+uLz2GcLyxjQUjYkpQf8AECiECAAuAg8BV5fKSSnTgau3ZyHEDuDlooHcXOAHIcTHaAO5gcD+qmu+oty8K0/5Y3s2u5p98ptJ3XlgyR4e++oAP07uSUff6pnTLqUkPDmckMgQ1ketJzkvGQedA/cG3MuQFkMqvEL2061nMUrJtP61m3bYZEajlmfHN6jKukSkwUDa7vMkbvTHkJ2Oc0drPD5YgkWbupP/pjaVG/SllHohxFRgI9qUzRVSyhNCiPeAMCnl2jLqnhBC/Iw26KsHppQ1c0dRasOVp/ynimWfdLe34tvHezD6891MWLGfn5/6Dy09qmCZf5ELmRcIjQwlNPLaFbLBLYLp5derUitkzyVm8XNYLON6NqOpa82nHK6UyO3abJeRy6rkdNn79xM/cxb5p05he9tteD06BOveI8Cijqw4rgPUxuhKo3YhJYd+83Ywtmcz3h3W4YbjkYlZjF6yBysLM1Y/fXvFNwIpJiUvhY3nNxIaGXrNgGxwi+ASB2Qrasr3h9h+OoE/X+2Hu309CXI/PAQXw2DaiZsKzIUXLxL/4TwyN2zAwscbr1dfxWHgwPoxc6mKqI3RFcUEJT3lF9fCw55vJnXnoWV7Gbd8Hz8/9Z8KBdRcfS7bY7YTGhXK7ou7qy1l8dHYNEKPXea5uwPrT8BPPQ9nNkDvlysd8I25uSQv/5Lk5ctBCNyfnYrbY49hZlP5m3NDp4K+0miV1Jdfkg4+Tnw18TbGfrmP8V/u58fJPcvc2PvKgGxIZAhbY7ZWakC2oj7ceBoXWx1P9KpAFsjadmC5lp+m26MVriqlJHP9euI/nIf+8mUcB9+L58svo/OpG9Ns6zIV9JVGq7yn/OKCmruyZGw3nvgmjElfH+DbST2wsfx3T9jKDMjm6w2sO3KZP45dJtDTnn5tPenWzAWdecVSG++OSOKvs0m8NaQdDtal34zqlIIcLT9NJbYRzAsPJ27mTHLDDmLVrh2+c+dge9tt1dTQhkcFfaVRMvUpv7i+bTz55MEuPPvjIZ787iDLxwcRn3uxwgOyCZl5fLc3hh/2RZOUVYCvsw1/nU1k6Z+ROFhb0Lu1B3e18aRvGw/cyumqkVIyZ+NpvJ2sGduzWaX+LGrF8V8gL61C0zT1KSkkfrKAtNWrMXd2psl//4vzqPsR5mVvyK5cSwV9pVFauM30p/zihnT2Jj7bn1l//kS/72eRQQRgWsrio7FpfLXrPCFHL6E3Su5q48mjdwRwRys3sgsM/H02kW2nEth+OpHQo5cRAm5t6sxdbTzp19aTDj6ONwxMbjwRz5ELacy9vzPWunoS/KTUNhjx7ADNbi+/eGEhqT/+SOLCRRizs3EZNxaPKVMwd6qDqaHrARX0lUYnJjmHXw9V7Cn/ygrZkMgQdl/cjXUTPal5XrR3eIhPhz6Kj0PJfcmFBiMbjsfx1a4oDsWkYWdpzpgezZh4e3Oau/+b78XeyoJBHb0Z1NEbo1Fy4lIG204lsO10Ah9vOcNHm8/g5WhFv6IbwJ2t3LGyMGPeptO09LDjvq71KLtJzF6IPwZDF5SbnyZr1y7iZ86i4Nw57O64A6/pr2PVqlUNNbRhUkFfaXRM7csvb4VsSBh8ti2Cr5zSeGOw9zVP4anZBfywP4bv9kZzOT2PZm62vBPcntFBfuX2u5uZCTr5OdHJz4nn+weSmJnPjtMJbD+dQOjRy6w6cAFLczMCveyJSMhiydiuWFRwHKBW7V8K1k7QaXSpRQpiYoifPYesbdvQ+fvjt3gR9v36NaopmNVFBX2lUSnvKb+0lMUlDci2HiDJyC3ki7+icLLRMfWuQE7HZfLVrih+/+ci+Xojd7Ry4/3hHenX1hNzs8oFLA8HK0YHNWV0UFMKDUYOnE9h+6kEtp1KoHdrDwZ2aHJTfyY1KuOSloa4x1NgeWNmS0NWNslLl5Ly9dcInQ6Pl17EdcIEzBp56oSqpIK+0qiU9pRfmRWyQgj+b2gHMvL0zNt0ho0n4jl2MR0rCzPu6+rLxNsDqjz/jc7cjNtbunN7S3feHFKB3aXqioNfa3nnr0t7LI1GMtatI2HeR+gTE3EaPhyPF19E51UPNn+pZ1TQVxqN65/yS0tZXJE9ZM3MBHNHdSZfb+DIhXReHdSGh2/zx8VOPZneQF8AYV9B4D3aDlNFco8dI/6DGeQeOYJ1p074ffYpNrfeWosNbdhU0FcajUXbIzAzL6Rtq3NM3fr51ZTFN7tCVmduxuIxdTudbp0Q/j/ITrg6TVOfmEjCx/NJ//13zN3d8Z45E6cRwxFm9Wh8oh5SQV9p8AxGAyFn/mTtpa+xaxXOjAN51b5CVinB/mXg2hLZ9E5SvvySpMWfYywowO3xSbg99RTm9lWX0E4pnQr6SoNU0oCsub01A5oNZHTb4RVOWaxUUk4KHP8VjqyCi2FkNnmShOEjKIiOxr5vX7xefw3L5s1ru5WNigr6SoNyZQ/Z0KhQotKj0JnpCPK8g4sRTXmwwwDe79ultpvY8Onz4ewmLdCf2QjGQvIt2hJ/oifZq9ZhGRBA0y+WYd+rV223tFFSQV+p964MyIZGhnI48TAA3by6Mb79eAY0G8DMddGQc5Gp/drVcksbMCkhNkzbBev4r1qKBXsvDJ0fI+mQJGXVBsysjXi+9hquYx5p9LtX1SYV9JV6qaQ9ZEsakC0+Y8fL0bTVt0oFpJ6Hoz9rT/Up58DCBtoFIzs+QPqRNBI+/BRDSgrOo+7H44UXsHBzq+0WN3oq6Cv1hsFoYF/cPkIj/10h62nrWeaA7MLtZzEzEzzdt2I5dpQy5KXDiTVaoI/Zrb3XvBf0ehHaDSMnPIL4N2aSd+IENl274rV0KTYdb9ygRqkdJgV9IcQgYAHadonLpZSzrzv+FDAFMABZwGQpZbgQojlwEjhdVHSvlPKpqmm60hiUtkJ2UMAgglsElzkgqz3lX2ScesqvOrsXwrb3QZ8H7q3h7neg0wPg3JTCuDgS3nqfjJAQLLy88PnwQxyDh6jUCXVMuUFfCGEOLAIGALHAASHEWilleLFiP0gplxSVHwZ8DAwqOnZOSqlWWigVcmVANiQy5OoK2V6+vQhuGUxvv94m7SG7cPtZzNVTftVJOgtb/g8C+sBdb4JPVxACY34+KUuWkLR0GRgMuD39FO5PPIGZbT3Zp7eRMeVJvzsQIaWMBBBCrAKGo212DoCUMqNYeTugbm28q9QLJQ3IeuraM6Lp89zf5l46efuYnL9GPeVXgw3TQWcLI5eCvYe2e9XmzSTMmUthbCwOAwbg+dqrWPr51XZLlTKYEvR9gQvFXscCPa4vJISYArwIWAJ3FTsUIIT4B8gA3pJS/lX55ioNTUkDss0dW2KfPZykuHZE6104d1Ty7abDWFkcpU0TB9o1caSdtwPtvB1p6+1Y4taF6im/ip3ZBBGb4Z4ZYO9B3pkzxM+cRc7evVgFBuL/9VfycHcAAAAgAElEQVTY9exZ261UTGBK0C/p0eqGJ3kp5SJgkRDiEeAtYAJwGfCXUiYLIboBa4QQHa77ZoAQYjIwGcDf37+Cl6DUN2WlLO7ichdv/ZxCZm4h304Moou/M+cSsjl5OUP7ictg88l4fgr79znE19nm6k2gnbcjrnaW6im/KukLYON0cAvE0OZBEt//gNRVqzCzt8frrbdweehBhIWaE1JfmPI3FQs0LfbaD7hURvlVwOcAUsp8IL/o94NCiHNAayCseAUp5TJgGUBQUJDqGmqATElZfOpyFuNX7McoJT8+0ZNOflrCs/Y+jrT3cbzmXAmZ+YRfuRFczuTkZW3TEWPRvx5LCzP1lF9V9i1BJkaQ5vY8iUOGYcjIwOWhB3F/9lksXFxqu3VKBZkS9A8AgUKIAOAi8BDwSPECQohAKeXZopdDgLNF73sAKVJKgxCiBRAIRFZV45W6r6QB2ZJSFh+MTuHRrw5gZ2XBt5N60sqz9DwsQgi8HK3xcrSmX5t/U+/mFRo4E6/dAJo42ain/KqQlUD2Tx8Rf7gF+Qmrse3eHa8338C6jcpXVF+VG/SllHohxFRgI9qUzRVSyhNCiPeAMCnlWmCqEKI/UAikonXtAPQG3hNC6NGmcz4lpUypjgtR6o7SUhZP7DCR/s3635CyeOeZRJ78NgxvJxu+ndQdP5fKzfqw1pnT2c+Zzn7ON30NChRevEj8Cw+TecwWXRNHfD+ZjcPAe9QUzHpOSFm3elOCgoJkWFhY+QWVOqW0FbLBLYLLTFkcevQyL/z0D4GeDnwzqTvu9uVPxVSqlzE3l+QvviB5+XIwFOA2sANuM7/HzFp9c6rLhBAHpZRB5ZVToy9KpZW3h2x5KYtX7Y/hjd+P0a2ZC8sn3FbiLByl5kgpyfjjDxI+nIc+Lg7HtjZ43pKLbvpXoAJ+g6GCvlIhFdlDtixLd55j1vpT9G3jwedjumFjaV4DrVdKkxceTtzMmeSGHcSqfTt8nx6M7ckPYPgisHYs/wRKvaGCfiOQrzfw55kkQo9eIiEzn67+LgQ1d6FrMxccrU17uq7MHrIlkVIyd+NpPt9xjuDO3nz8wK1YWqi89rVFn5JC4icLSFu9GnNnZ5q891+cgwciPu8BPl3glkfKP4lSr6ig30AV6I3sOpdEyJHLbAqPIzNPj7OtDj8XGz7feQ7DdokQ0MbLgduauxLU3IWg5q74OttcPUdFB2TLYzBK3v7fcX7YF8OYHv68N7yjyStslaolCwtJ/eEHEhcuwpibi+v4cbhPmYK5oyNsfQ8yL8MD34DaurDBUUG/AdEbjOyJTCbkyGU2nIgjPbcQB2sLBnZoQnBnb+5o5Y7O3IycAj2HY9I4cD6VsOgUfjsUy7d7owHwdjbD3+88edb7ic75B4M0EOgSeFN7yIJ2E3rx58OEHL3MM31b8srANmoWSC3J+nsX8bNmUXDuHHZ33onX9Nexalm0piElSkuq1vlBaNq9dhuqVAsV9Os5g1GyPyqFkKOX2HA8juTsAuytLBjQ3oshnbzp1dodK4tr+8ttLS24vZU7t7dyByCvsIDfTu7kfxHrOJ25h3BDHsY0JwrT78QqrxsOhg5k27mwV2/AzS4Rd3sr3O0tcbGzRGde/pNgboGBp78/yI7TiUy/ty1P9lGLpmpDQXQ08XPmkrVtGzp/f/wWL8a+X99rb76b3gIzC+j/bm01U6lmKujXQ0aj5GBMKiFHLvHH8TgSM/Ox0ZlzdztPgjv70LeNB9a6sgdGpZSEJ4cTEhnC+qj1JOcl46BzYETgEAYHDMbLsj2Hir4NHDyfysebz5R4HmdbHe72VrjZWWr/tf/3v252VrjY6pi36TRh0anMuq8TD3dXaTZqmiErm+SlS0j5eiVCp8PjpRdxnTABs+t3r4rcAadC4K63wdGnVtqqVD8V9OuozLxCLqXlcTEth4tpeVxMzeVSWi4X03I5n5RNcnYBVhZm3NXWkyGdvbmrrSe2luX/dZY0INvHrw9DWgy5YUC2mZs9I7toGROz8vUkZOSRnF1AclY+SVkFJGXlk5xVQHJ2PkmZBZyMyyA5q4D03MJrPlNnLlj4cFeGdK5c15BSOdJoJH3tWhI/+hh9YiJOw4fj8dKL6Dw9byxs0MP618G5Gfxnas03VqkxKujXkkKDkWMX068J5pfScoktep2Rp7+mvM5c4O1kg6+zDf3aetIr0J2723lhb1X+X2FVDMjaW1lg72FPC4/yr61AbyQlu+imkF2Ar7M1rTwdyq+oVJnco0eJmzGDvCNHse7cGb+Fn2Fzyy2lVwhbAYkn4cHvQKfm5DdkKujXggK9kQkr9rMnMvnqe47WFvg42+DnYkP3AG0WjY+zDb4uWqD3sLfCrAIzXa6skA2JDGH3xd3opf7qgOyQFkNoYtekOi4N0JKdNXGypomTCh41rTAhgcSP55O+Zg3mHu54z5qF0/BhiLJm4eSkwPYZ2uYobYNrrrFKrVBBv4ZJKXnz92PsiUzmrSHt6BXogY+zNQ4mzpcvy82ukFXqL2NBAanffEPS4s+RhYW4PfE4bk8+hbm9XfmVt30A+ZkwaDaoGVUNngr6NWzpn5GsPhjLc3cH8nivFjd9vtIGZCu6Qlapn6SUZG3fQfyc2RRGx2B/1114vfYqls2amXaCuONw8Cu47XHwal+9jVXqBBX0a9CG43HM2XCK4M7eTOsfeFPnqsiArNIw5UdGEj9zFtl//41lixY0/eIL7HvdafoJpIQNr4O1E/SdXn0NVeoUFfRryLHYdF746R9uberMvNG3VGphUlWvkFXqJ0NGBkmLFpPy/feY2djgNf11XB55BKGrYBdh+P/g/F8weB7YulZPY5U6RwX9GnA5PZdJKw/gZmfFsnFB5c6hL660Adlp3aYxOGBwtQ7IKnWLNBhI+/VXEj9ZgCE1FedRo/B44Xks3NwqfrLCXNj0Nnh2gG6PVn1jlTpLBf1qlp2vZ9LXYeQUGPj16R54OJTf7aIGZJXr5Rw8SNyMGeSHn8SmWze8vliGTYcOlT/h7s8gPQYmrANzFQYaE/W3XY0MRskLPx3mVFwGX068jTZNSp+rrgZklZIUXr5MwryPyAgNxaJJE3w+mofj4ME3l7cocif89TG0GwYBvauusUq9YFLQF0IMAhagbZe4XEo5+7rjTwFT0LZEzAImSynDi45NByYVHXtOSrmx6ppft83ZcIrN4fG8O7T9NXu5FqcGZJWSGPPySF6xguQvloPRiPszT+P2+OOY2VZuK0ntpAb480PYMRvcA7UpmkqjU27QF0KYA4uAAUAscEAIsfZKUC/yg5RySVH5YcDHwCAhRHu0jdQ7AD7AFiFEaymloYqvo85ZtT+GZX9GMv4/zZh4R8A1x1LyUth4fiOhkaFqQFa5hpSSzE2bSZg7l8KLF3G45x48X30VSz/fmztxVgL89oSWX6fzgzDkY7AqffN5peEy5Um/OxAhpYwEEEKsAoYDV4O+lDKjWHk74MrGu8OBVVLKfCBKCBFRdL49VdD2Omt3RBJvrTlO79YevBOszX3O1eeyPWY7oVGhakBWKVHe6TPEz5xJzr59WAUG4v/119j17HHzJ476C36dBHnpMOwz6DJOLcJqxEwJ+r7AhWKvY4Eb/iUKIaYALwKWwF3F6u69ru5NPrLUbecSs3jqu4MEuNvxyUOd2B+3l5DIELbGbFUDskqJ9KmpJH32GamrfsLcwQGvd97G5YEHEBY3OeRmNMJfH8GOmeDaAsb+Bk06Vk2jlXrLlH9VJT0SyBvekHIRsEgI8QjwFjDB1LpCiMnAZAB///qbejc1u4BHv96Puc1FegZd5r51/6cGZJVSSb2e1J9+IunTzzBkZuLy0EN4PPcs5s7ON3/y7CStO+fcNug4CoZ+AlYq6Z1iWtCPBZoWe+0HXCqj/Crg84rUlVIuA5YBBAUF3XBTqA/OpUbz+K/LSHbeg5llIqHRldtDVqlFhXmwawF0GQNOftX6Udl79xE/cyb5Z85g27UzXm++iXWHzlVz8vO7tO6cnBQI/gS6TVTdOcpVpgT9A0CgECIAuIg2MHvNbslCiEAp5dmil0OAK7+vBX4QQnyMNpAbCOyviobXBTcMyOqghV0nJnaeogZk66P9y7SukMjtMDEUzExfRGeqgtiLJMyZQ+bmzeh8ffF9dTwO0bMRv++EI30g8B5oPbByNx2jEXbN1xKouQTA4z+DdxXdSJQGo9ygL6XUCyGmAhvRpmyukFKeEEK8B4RJKdcCU4UQ/YFCIBWta4eicj+jDfrqgSn1feZOSQOyrjp/8hMGMa7jCN6+947abqJSGXnp8PfH4OgLMXtg72K4/dkqO70xJ4ekL74g5csVYG6Ox/PP4To6GLMVfaBJJ2h2O5zZoP2EAl4dteAfOBD8gsq/AWUnw++TIWILdLgPhi4Aa8cqa7/ScAgp61ZvSlBQkAwLC6vtZlzjygrZ6wdkB7cYjLO+B+/9nsawW3xY8NCtarPv+mrbDPhzLkzeqc1lP7sZnvoLPG5usF1KSUboHyTMm4c+Lg7H4GA8X34JnZcX/Piw1uf+5J/g2VZLgJZ0pij4b9JuPtIANq4QOEC7CbS8G2yu6/OP2QurH4WcJBg0C4Imqe6cRkgIcVBKGVReObUitxSmrJANPRrHy78foau/M3NHdVYBv77KSoQ9i6D9CPC5FYLnw6Ie8PuTMGlLpdMU5J44QfyMmeQeOoR1+/b4fvwRtl27agcPfQNn1sPAmVrABy1Qe7TRfu54HnJTtZvCmU1wdhMc/QmEOfj/B1rfo30LOLMBtr4Hzv7w+BbwLmN3LEVBPenfoKQVstcPyBqMkrkbT7F0ZyRBzVxYOq4bbvZqoLbeWv+61p8/ZZ+2UhXgxBpYPQH6vQV9XqnQ6fTJySR+8glpv/yKuYsLHtNewPm++xDmRV00KVGw5E7w6QLj10JZu1pdYTTAxYP/fguIP/bvsfbDtfn31moMqTFTT/oVkJqXyobzG25YITuhwwQGNBtwzYBsek4hz676hz/PJDKmhz//N7QDlhZqCma9lRYDYV9qM3bci+1x0GEEnBwFO2dr3SomDIjKwkJSvv+epEWLMebm4jphAu5TnsHcodhUSaMB1jwNwgxGfG5awAetT79pd+3n7ncgPVbrgrJ21Prw1bdMxUSNNuiXNCDbyrkVL3R9gcEBg/G2976hzpn4TJ74JoxLabnMHNmJR3rU3zUFSpEdswEBfV6/8djgD7V882uehie2gUXp3+ay/vqb+FmzKIiMxK5XL7ymv45VixJ2Rtv9mdZXP2IJODe98bipnPwgSKVEViquQQV9KWWZ/eqlDciaskJ2w/HLvPjzEeysLFg1uSfdmqlNJ+q9hFNw5Efo+Qw4lbBQ3NYVhn4KPz4IO+doT9jXKTh/nvjZc8jasQNdM3/8Pl+Mfd++Jf87jDuubUDebijc8lA1XJCilK/BBP28QgMjF+9mSKcmPNKjGa52lsDNpyw2GiWfbDnDp9siuKWpM0vHdqOJk3VNXZZSnbZ/ADo7uPPF0su0GQRdxsLf86HNYG36JGDIyiJ5yRKSV36DmU6H58sv4TJ+PGaWliWfR58Pv00Ga2cIXqC6Y5Ra02CCfkp2Ae72lszbdIbPtkUw8FYdTXzC2RO/udQB2fJk5BUybdVhtp5KYHQ3P94f0bFCu14pddjFg3ByHfR9A+zK2Xlq4Ew4twN+fwo5eSfpoZtImP8xhsQknEaOxPPFaVh4eJR9ju0zIOEEPPJz+Z+nKNWowQR9H2cbFowJ5NujJ/n19Fq2Z5+Fs2Av2/BIy5d45raROFVgdsO5xCye+CaMmOQc3hvegXE9m6kpmQ3Jlv+CrTv855nyy1o7wfCF5H4ymrihA8iLScXmllvwWrwYm06dyq9/fhfs+hS6TtAGhRWlFjWYoH8h4wLD1gy7mrL4Ad9nyUzqxG9h2Sw9lc/2sMM8dkcAI7r4lvu0vvVkPC+sOozOwozvHu9BzxbqyaxBidwBUTth4CyTkpAVxieQuGwD6Zs9sLBOwuflyTg+9iLClJk3eRmw5ilwaaZ9Y1CUWtZg5ulLKVlxfAV3+t55zYBsvt5AyJHLfPl3FOGXM3C1s2RMD3/G9WyGp+O1ffNGo2TR9gg+3nKG9t6OLBsfhK+zzU1fk1KHSAnL74bMeHj2IOhKH58xFhSQ8vVKkpcsQRYW4jp+LG58j7lOwtO7TduE5H9T4PAP8Oh68O9ZhReiKNcydZ5+gwn65ZFSsi8qhS//jmLLyXgszARDb/Fh0p0BdPBxIitfz8s/H2HDiThGdvFl1n2dVP99Q3RyHfw0FoYthK7jSiwipSRr+3biZ8+hMCYG+7vvxuu1V7H099dSHqwYpE2XDJ5f9medCoVVj8Cd06D/u1V+KYpSnAr6ZTiflM3Xu8/zc9gFcgoM9GzhSkp2AREJWbwxuB2T7gxQ/fcNkdEAn98O0ghP7ykxvUJ+RATxs2aTvWsXli1b4vXGdOzvuC6J3qa3tPn2Y3+DVneX/FlZibC4Jzh4F83xL2VWj6JUERX0TZCeW8hPB2JYuTua7AI9Cx/uyp2B7jXy2UotOPyDttDqgW+01AXFGDIySFy4kNTvf8DM1haPZ6fi8vDDCJ3uxvMU5sHS3lCQpXXzXJ8ATUpYNQYiNmsJ3LzaV+NFKYpGBf0KMBglBXojNpaqO6fB0ufDZ93Azh2e2H51nrw0GEj75VcSP/kEQ1oazqNH4/HC81i4lrP47uIhWN5f22R85OfXHvvnO60v/54PqjQ9s6KUReXeqQBzM6ECfkMX9hWkX4Bhn14N+DlhYcTNmEn+yZPYBHWjyZtvYt2unWnn8+0KvV7S0jG3GwptB2vvp0ZrCdya3Qk9p1TTxShK5amgrzR8+VlajvzmvaBFPwovXybhww/J+GM9Ft7e+H78EQ733lvxcZzer2jpkdc9D017aN08a57Wjo2sQDI1RalBKugrDd/ezyEnCeOd00levJjkL5aDlLg/8wxuTzyOmU0lp+VaWMLIpbC0D/zxEvh0hehdMHyxlt9eUeogk4K+EGIQsABtu8TlUsrZ1x1/EXgcbUvEROAxKWV00TEDcCX5d4yUclgVtV1RypeTgtz1KZmG20mY/A6Fly7hMGgQXq+8jM63hCRrFeXVAfpN1zYyCf8ftA2GWx8pv56i1JJyg74QwhxYBAwAYoEDQoi1UsrwYsX+AYKklDlCiKeBucCDRcdypZS3VnG7FcUkeT/9H/EbdOQknMeqTRv8V67Erkf3qv2Q25+H0xsg9by2N62a7qvUYaY86XcHIqSUkQBCiFXAcLTNzgGQUm4vVn4vMLYqG6koFaVPTSVx3mzSftuMubUdTf7vDZxHj0ZYVEOPprkFTFgHhTlaOmZFqcNM+T/AF7hQ7HUs0KOM8pOA9cVeWwshwtC6fmZLKddUuJWKYiKp15O66icSP/sMY2YGLoG5eHz6C+bNO1bvB+usy0zpoCh1hSlBv6TvqiVO7hdCjAWCgD7F3vaXUl4SQrQAtgkhjkkpz11XbzIwGcDfXw2AKZWTvWcP8TNnkn82AtugW/Dy2Ir1gIlQ3QFfUeoRU+aUxQLF93XzAy5dX0gI0R94Exgmpcy/8r6U8lLRfyOBHUCX6+tKKZdJKYOklEEe5eUlV5TrFMTGEvvss8Q8+hjG3Dx8P/sU/+E2WLtbaNMqFUW5ypSgfwAIFEIECCEsgYeAtcULCCG6AEvRAn5CsfddhBBWRb+7A3dQbCxAUW6GMTubhE8+IXLwELL+3oXHCy/QIjQEx46eiBO/Qc+nwd6ztpupKHVKud07Ukq9EGIqsBFtyuYKKeUJIcR7QJiUci3wIWAPrC5a4HJlamY7YKkQwoh2g5l93awfRakwKSUZISEkzPsIfXw8jkOH4vnyS+i8vLREZ78/rW1LePtztd1URalzTJrKIKX8A/jjuvfeKfZ7/1Lq7QZM2FpIUUyTe/wE8TNmkPvPP1h36IDv/PnYdi3qMcxKgJVDtVQID/94YyI0RVHUilylftAnJ5Mwfz7pv/6Guasr3jM+wGnkyH93r8qM0wJ+eiyM+RkCetdugxWljlJBX6nTZEEBKd99T9LixRjz8nCdOBH3Z57G3KHYNocZl+DrYC3wj/kFmt9R+gkVpZFTQV+ps7L+/JP4WbMpiIrCrncvvF6fjlWLgGsLpV3QnvCzk2Dc7+Bf1hISRVFU0FfqnPyoKBJmzyFr504smzen6dIl2Pfpc2PB1GhYGQy5aTB+DfiVm0pcURo9FfSVOsOQlUXS4s9J+fZbzCwt8XzlFVzHjUVYlrDVYEqU9oSfnwHj/6flt1cUpVwq6Cu1ThqNpP++hoT58zEkJeF03314TnsBi9IW6iWf0wJ+YQ6MXws+Kp+fophKBX2lVuUePkzcjJnkHTuGzS234PX5Ymw6lTHLN+msFvANBVqSsyZqRrCiVIQK+kqtKIxPIOGjeWSsXYeFpyc+c+fgGBz87xTMkiScgm+GgTTChBC14biiVIIK+kqNMubnk/L1SpKWLoXCQtwmT8b9ycmY2dmVXTE+XAv4wgwmhoJHm5ppsKI0MCroKzVCSknWtm3Ez55D4YUL2N99N16vvYqlKVlV447BN8PB3FLr0nEPrP4GK0oDpYK+Uu3yIyKInzmL7N27sWzVkqZfLsf+DhMXUF0+ogV8na0W8N1aVm9jFaWBU0FfqTaG9HQSFy0i9fsfMLOzw+uNN3B5+CGETmfaCS4egm9HgJWjFvBdA8qvoyhKmVTQV6qcNBhIW/0LiQsWYEhLw/mBB/B4/jksXE3YStBogOhdcGw1HPsV7Ny1gO/SrPobriiNgAr6SpXKOXCAuJmzyD95EtugILzefAPrdu3KriQlXD4Mx36B479C5mWwtId2Q+Hut8HJr2YaryiNgAr6SpUovHSJ+A8/JHP9Biy8vfGd/zEOgwZRtL9CyZLOaoH+2GpIOacN1AbeAx3vh9aDwNK25i5AURoJFfSVm2LMzSX5yxUkL18OUuI+ZQpuj0/CzMam5AoZl+D4b1qgv3wYEBDQC+58QXuyt3Gp0fYrSmOjgr5SKVJKMjdsIP7DD9FfuozDvYPweuUVdD4+NxbOSYGTa7Wn+vN/AxJ8usDAmdDhPnD0rvH2K0pjZVLQF0IMAhagbZe4XEo5+7rjLwKPA3ogEXhMShlddGwC8FZR0Q+klCurqO1KLck7dYr4GTPJOXAAq7Zt8Zk9G7vu3W8sWJgHm9+GsK/AWAhugdB3OnQapaZeKkotKTfoCyHMgUXAACAWOCCEWHvdXrf/AEFSyhwhxNPAXOBBIYQr8H9AECCBg0V1U6v6QpTqp09NJXHBAtJ+Xo25oyNN3n0X59GjEObmNxZOiYLVE7R59t0ehW4TwfsWKKuPX1GUamfKk353IEJKGQkghFgFDAeuBn0p5fZi5fcCY4t+HwhsllKmFNXdDAwCfrz5pis1RRYWkvrjKhIXLsSYnY3LmDF4TJ2CuZNTyRVOhsCaZ0AAD/0IbQfXaHsVRSmdKUHfF7hQ7HUsUNb2RJOA9WXU9a1IA5Xalb17N3EzZ1IQcQ672/+D1/TpWAWWkgbBUAhb3oU9C7U++9Er1fx6RaljTAn6JX0flyUWFGIsWlfOlW2OTKorhJgMTAbwNyUXi1LtCi5cIH72HLK2bkXXtCl+ixZif9ddpU/BTI+F1Y9C7H7oPhnu+QAsrGq20YqilMuUoB8LNC322g+4dH0hIUR/4E2gj5Qyv1jdvtfV3XF9XSnlMmAZQFBQUIk3FKVmGLOzSVq6jJSvvgKdDo9p03CdOAEzqzIC+Nkt8NsTWo77USu0efaKotRJpgT9A0CgECIAuAg8BDxSvIAQoguwFBgkpUwodmgjMFMIcWXy9T3A9JtutVLlpJRkrFtHwryP0Cck4DR8GB4vvojOy6v0SkYD7JgFf84Dz/bwwDfg3qrmGq0oSoWVG/SllHohxFS0AG4OrJBSnhBCvAeESSnXAh8C9sDqoq//MVLKYVLKFCHE+2g3DoD3rgzqKnVH7rFjxM+YSe7hw1h37Ijvgk+w7dKl7EqZ8fDrJDj/F3QZC/d+qFbQKko9IKSsW70pQUFBMiwsrLab0SjoExNJmP8J6b/9hrm7O57TpuE0ckTZu1cBRP2lBfy8DBjyEXQZUzMNVhSlVEKIg1LKoPLKqRW5jZAsKCDl2+9IWrwYY0EBrpMew/3ppzG3ty+7otEIf38M22eAa0sYt0ZtWago9YwK+o1M5o4dJMyaTUF0NPZ9+uD5+mtYBZiQpz4nBX6bDBGbtYHaoQvAyqH6G6woSpVSQb+RyI+MIn72LLL//AvL5s1pumwp9r17l18xNxX++R72LIKcJK07J2iSWlmrKPWUCvoNnCEzk6TFn5Py7beYWVvj+eqruI4dg7C0LLtifDjsXwpHf4bCHPD/Dzz8g7boSlGUeksF/QZKGo2k//47CR/Px5CSgtP99+H5wgtYuLuXXsmgh9N/wP5l2qwcC2voNFpbbOXdueYaryhKtVFBvwHKOfQP8TNmkHfiBDZduuC1dCk2HTuUXiE7GQ6thANfQkYsOPlD//9C1/Fga8IWh4qi1Bsq6DcghfHxJMz7iIx167Dw9MTnw7k4BgeXnjrh0mHtqf7YL2DIh4DecO8caHMvmJWQOVNRlHpPBf0GwJifT8pXX5O0bBno9bg99STuTzyBmZ3djYUNhRD+Py3YX9gHOjttcVX3J8CznL1sFUWp91TQr8eklGRu2ULCnLkUxsbiMKA/nq++iqXn/7d379FVlWcex78/kkDIBRISEsJNASlaL4MYqS2txbZcKla6KlZHqzh21bGtjC7LtGjHcdAiKrRVO+CIgNjWyrJUWKyqoEvb1elMbQmIQwkXDbUYIMmBgCYEcuOZP/YBY4CVcyBk55zzfNY661z25Tybyy877373+8412L4AAAx9SURBVOZDfRVE9kBdFdRXB891VfC3PwTL8ofBpLkw+gbonRf2oTjnuoiHfqI5uA8iW2jcvJGqxato2LKLXkW9GXpdCdmF/wu/uAya6o/fLq0X5BbDwNFBl8tzvgQd3XnrnEs6HvrdVUsjRLZB9Wao2Rw8V5fTWltNZFMu+9/NpkeGUXzpYfLH9ER9e0HuhXDOhCDcc0sgpxhyBwSPzDzvW++c89APnRkc2Ak15VD916B/fE057H0HrDVYJ60XVvAJDkRGEXk9i9aGRvKmTqb/nXeRPmCIh7lzLmYe+mHZtCK4mFqzBRo//OjzvLOg+Hw496rgufh8DlbUUj33ERq3bSPr0ksp/uG9ZJ57bni1O+cSlod+VzMLBi17/YFgDPqLrgsGLSu+APqfC5l9jq3avGsX1XPmU7dmDekDSxj02GPkTpp48i6YzjnXAQ/9rnSkFdbcEwxvcOG1MHUhpB8/HMKRQ4fY9/Ri9i1ZAhKFM+6g4JvfpEdmZghFO+eSiYd+V2lpDEapLF8Fn74DJjx4XO8ZM6PulVeonjeflj176HPllymaOZOMgQNDKto5l2w89LvC4Q9g+Y3BeDYTfwSfmXH8KuXlVD30EIfK1tPrvPMYNO9Rsko7nA/BOefiElNHbUmTJW2T9K6kWSdYfrmkDZJaJE1rt6xV0sboY3VnFZ4w6qrgmSmw80/wtaePC/yW2lr2/Pv9/O2aaTRV7GDA7NkMW/FrD3zn3BnR4Zm+pDRgATABqATWSVptZuVtVtsJ3ALMPMEuDpnZ6E6oNfHsfQd+8TVo2Ac3vADnfPHYImtuZv/zzxP5zwUcaWig3803Ufid75DWt2+IBTvnkl0szTtjgXfNbAeApOXAVOBY6JvZe9FlR85AjYmpsgyeuxbUA275LQwac2xR/R//h+q5c2mqqCB73DiK772HXiNGhFiscy5VxBL6g4D327yvBD4Vx3dkSioDWoCHzWxVHNsmpu2vwq+nQ04RfONFKAgCvWnnTqoffoT6N94gY+hQBi9cQM4VV3gXTOdcl4kl9E+USBbHdww1s92ShgNvSNpkZhUf+wLpNuA2gKFDh8ax627oredg9QwYcAHcuAJyimitP8i+p56idtkylJFB/+/dTb/p0+nR0exVzjnXyWIJ/UpgSJv3g4HdsX6Bme2OPu+Q9HvgYqCi3TqLgEUApaWl8fxA6T7a3nQ1fDxc90ssI5sPVq0i8uOf0BKJ0HfqVPrffTcZxUVhV+ucS1GxhP46YKSkYcAu4Hrghlh2LikfaDCzRkmFwDjg0VMttttqe9PVBdPgq09yqHwrVXPmcPjt/yPzoosY/LMn6D06Na9nO+e6jw5D38xaJN0BrAXSgKVmtlnSA0CZma2WdCmwEsgHviJptpmdD5wHPBW9wNuDoE2//CRflZja3nR12XdpvvhOIvf9Bx+sXEla/0JK5s6l79SrkQ9j7JzrBmTWvVpTSktLraysLOwyOtZYDxWvw58WwvtvcuSK2ex/J4e9C5/kSHMzBdNvpuD220nLyQm7UudcCpC03sw6vMHH78iNR30Etr8CW1+Cit9BayOWmU/9WTOpfugVmv++k5zx4yme9QN6nn122NU659xxPPQ7sq8iCPltL8PONwGDvkOh9FYas8dQ/fNXObjsV/QcPpwhTy8i53OfC7ti55w7KQ/99sxg91tB0G99CSJbgs8HXAjjZ8GoK2nNPpu9CxZS+9yD9MjMpGjWD+h3440oIyPc2p1zrgMe+gDNh4OxcY6e0X+4C5QGZ30GLnkYRl0J+Wdhra0cePFFIj/9Nq3795M37Rr633UX6QUFYR+Bc87FJPVCv6khmG92z0bYvRH2vB2czR9pgfTewfg4X7gPPjEJsvod26xhwwaqfzSHw+Xl9B4zhuKnF9H7/PNDPBDnnItfcod+Yz1UbQqCfc/RgN8KFh0iKKsASkbDyAkwZCwM+zz0zPrYLpqrqqiZN58PX3qJ9OJiBs6fT58pV/rQCc65hJQ8od/SCJXrPjp737MxGOXy6IgR2UUwcHQw92zJPwSv+ww66aTiRxobqV26lL2LnobWVgq+fTuF3/oWPbKyTri+c84lguQJ/UP7YdmU4HXuwCDYL7gmeC4ZDX1KYtqNmVH32mvUPPIozbt2kTtxIkXf/1d6Dh58Bot3zrmukTyhnzsgGNFywIXB6Jan4PC27VTPnUvDm2/Sa+RIhi57huzLLuvkQp1zLjzJE/rwsUlK4tF64ACRJ37G/uXL6ZGbS/F9/0b+ddeh9OT643HOuZRONWtpYf8LL7D38Sdorasj//rrKZxxB+n5+WGX5pxzZ0TKhv7BP/+F6jlzaNy+nayxYyn+4b1kjhoVdlnOOXdGpVzoN1XuombePOrWriVj4EAGPf44uRMneBdM51xKSJnQP9LQwL7Fi9m3ZClIFP7LDApuvZUemZlhl+acc10m6UPfzPjw5ZepmTeflqoq+kyZQtHM75FRElsXTuecSyZJHfqHy8upmvMQh9avp9cnz2PQj+eTdcklYZflnHOhScrQb6mtJfLTxziwYgVpeXkMeGA2eddcg9LSwi7NOedCFdMcfpImS9om6V1Js06w/HJJGyS1SJrWbtl0Se9EH9M7q/ATseZmap99lopJkzmwciX9br6JEWvXkP/1r3vgO+ccMZzpS0oDFgATgEpgnaTV7ea63QncAsxst20/4H6glGAQnPXRbfd3Tvkfaaqs5P3b/pmmHTvI/uxnKb5nFr1GjOjsr3HOuYQWS/POWOBdM9sBIGk5MBU4Fvpm9l502ZF2204CXjOz2ujy14DJwPOnXXk7GUVF9BwyhKKZM8m5Yrx3wXTOuROIJfQHAe+3eV8JfCrG/Z9o20ExbhsX9ezJkKf+60zs2jnnkkYsbfonOmW2GPcf07aSbpNUJqksEonEuGvnnHPxiiX0K4Ehbd4PBnbHuP+YtjWzRWZWamal/fv3j3HXzjnn4hVL6K8DRkoaJqkncD2wOsb9rwUmSsqXlA9MjH7mnHMuBB2Gvpm1AHcQhPUW4AUz2yzpAUlXA0i6VFIlcC3wlKTN0W1rgQcJfnCsAx44elHXOedc15NZrM3zXaO0tNTKysrCLsM55xKKpPVmVtrRejHdnOWccy45eOg751wK8dB3zrkU0u3a9CVFgL+fxi4Kgb2dVE4YvP7wJfoxeP3hC+MYzjKzDvu8d7vQP12SymK5mNFdef3hS/Rj8PrD152PwZt3nHMuhXjoO+dcCknG0F8UdgGnyesPX6Ifg9cfvm57DEnXpu+cc+7kkvFM3znn3EkkTeh3NKVjdydpqaQaSX8Nu5ZTIWmIpN9J2iJps6Q7w64pHpIyJf1F0tvR+meHXdOpkJQm6S1Jvw27llMh6T1JmyRtlJRw47FIypO0QtLW6P+FT4ddU3tJ0bwTndJxO22mdAT+sd2Ujt2apMuBeuDnZnZB2PXES1IJUGJmGyTlAuuBrybK34GCqdayzaxeUgbwR+BOM3sz5NLiIulugulJ+5jZVWHXEy9J7wGlZpaQ/fQlPQv8t5ktjo5KnGVmB8Kuq61kOdM/NqWjmTUBR6d0TBhm9gcgYUcgNbM9ZrYh+rqOYETWMzJL2plggfro24zoI6HOiCQNBqYAi8OuJRVJ6gNcDiwBMLOm7hb4kDyh32XTMrqOSTobuBj4c7iVxCfaNLIRqCGY2zmh6gceA74PtJ+rOpEY8Kqk9ZJuC7uYOA0HIsAz0Sa2xZKywy6qvWQJ/dOZ0tF1Ikk5wG+Au8zsw7DriYeZtZrZaIIZ3sZKSphmNklXATVmtj7sWk7TODMbA3wZ+G602TNRpANjgCfN7GLgINDtri8mS+ifzpSOrpNE28J/AzxnZi+GXc+piv5K/ntgcsilxGMccHW0TXw58AVJvwy3pPiZ2e7ocw2wkqDpNlFUApVtfkNcQfBDoFtJltA/nSkdXSeIXghdAmwxs5+EXU+8JPWXlBd93Rv4ErA13KpiZ2b3mNlgMzub4N//G2b2jZDLiouk7GgnAKLNIhOBhOnNZmZVwPuSRkU/+iLQ7ToypIddQGcwsxZJR6d0TAOWmtnmkMuKi6TngfFAYXTqyfvNbEm4VcVlHHATsCnaLg5wr5m9HGJN8SgBno32BOtBMC1oQnZ7TGDFwMrg/IF04FdmtibckuI2A3guevK5A/inkOs5TlJ02XTOORebZGnecc45FwMPfeecSyEe+s45l0I89J1zLoV46DvnXArx0HfOuRTioe+ccynEQ98551LI/wOci2zENOGryQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c7667f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# least squares fit\n",
    "a,b = np.linalg.lstsq(np.vstack([Sigmas, np.ones(len(Sigmas))]).T,MSEs10)[0]\n",
    "print(\"Slope of MSE10:\", a)\n",
    "MSEs10fit = b + np.array(Sigmas)*a\n",
    "\n",
    "# least squares fit\n",
    "a,b = np.linalg.lstsq(np.vstack([Sigmas, np.ones(len(Sigmas))]).T,MSEs20)[0]\n",
    "print(\"Slope of MSE20:\", a)\n",
    "MSEs20fit = b + np.array(Sigmas)*a\n",
    "\n",
    "\n",
    "plt.plot(Sigmas,MSEs10)\n",
    "plt.plot(Sigmas,MSEs20)\n",
    "plt.plot(Sigmas,MSEs10fit)\n",
    "plt.plot(Sigmas,MSEs20fit)\n",
    "plt.show()\n",
    "\n",
    "mtx = np.array([Sigmas, MSEs10, MSEs20, MSEs10fit, MSEs20fit])\n",
    "np.savetxt(\"./denoise_learned.dat\", mtx.T , delimiter='\\t')\n"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}