Convolution of functions with Tensorflow

tf-function-convolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiable convolution of two probability distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We draw samples from a uniform distribution, smear them with a Gaussian. \n",
    "\n",
    "Then we perform a maximum likelihood fit to a probability density function which is a convolution of a uniform distribution with a Gaussian where we fit the resolution parameter $\\sigma$ of the Gaussian to check whether we get back the size of the Gaussian smearing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "session = tf.InteractiveSession()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a grid of points on which we will work: mainly plotting but also calculating the convolution of the functions later on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also define the bin centers (at which the functions to be convoluted are evaluated). We make sure that these are centered at zero (which we need for at least for one of the convoluted functions) and that we have an odd number of bins. In this example we use the same grid for both functions in the convolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbins = 10001\n",
    "xmin = -5\n",
    "xmax = +5\n",
    "\n",
    "bin_width = (xmax - xmin) / float(nbins)\n",
    "\n",
    "# bin boundaries\n",
    "xboundaries = np.linspace(xmin, xmax, nbins + 1, dtype = 'float32')\n",
    "\n",
    "# bin centers\n",
    "xpoints = xboundaries[:-1] + 0.5 * bin_width"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the uniform distribution function\n",
    "\n",
    "see also https://www.tensorflow.org/api_docs/python/tf/distributions/Uniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "uniform = tf.distributions.Uniform(0.,2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAF3CAYAAACfa4MXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFI1JREFUeJzt3W+spGd93+Hvj10vRrYhIqYL9a6ylmJV2gDBcDBVjMwR0MpukP0iiWRHabFKta2UbUmTqjKl8gu3UhoikVaq27JKkFAa4tC0VbftIjcpmSKgdtfGxtHaLJwsBq9Ja4P5t078Z83dF3uIDssedjxn5jz37HNdkqXzzDya+c0tjz+emXPuqdZaAIB+vGToAQCA7yfOANAZcQaAzogzAHRGnAGgM+IMAJ0RZwDojDgDQGfEGQA6I84A0JmdQ93x5Zdf3vbt2zfU3Q/i6aefziWXXDL0GEvNGm6dNdw6azgfY1vH+++//2uttVdNc+5gcd63b1/uu+++oe5+EJPJJKurq0OPsdSs4dZZw62zhvMxtnWsqi9Pe663tQGgM+IMAJ0RZwDojDgDQGfEGQA6I84A0BlxBoDOiDMAdEacAaAzU8W5qq6vquNVtVZVt53j+lur6smqenD9n78z/1EBYBzOu31nVe1IcmeSv5bkZJKjVXW4tfbwWaf+Xmvt4AJmBIBRmWZv7WuSrLXWTiRJVd2V5KYkZ8cZ6NzzL3w3D3/9hbTjTww9ylL7k2+8kNWhh+CCNk2cr0jy2Ibjk0neco7zfqaqrkvyhST/sLX22NknVNWBJAeSZPfu3ZlMJi964GV26tSp0T3mebOGW/O/v3o6H3ro2eTo0aFHWXqX7vpEXn2JX9vZCs/nzc3rW6n+a5Lfba09W1V/N8lHkrz97JNaa4eSHEqSlZWVNqZvI0nG9w0si2ANt+ar934leeiP81vvXskrL9k19DhL6d4vPZV/8fHP53VXvyk/8ZdfMfQ4S83zeXPTxPnxJHs3HO9Zv+wvtNa+vuHwN5N8YOujAYvy2itekd0vv3joMZbSE995dugRGIFp3pM5muSqqrqyqnYluTnJ4Y0nVNVrNhzemOSR+Y0IAONy3lfOrbXTVXUwyd1JdiT5cGvtWFXdkeS+1trhJP+gqm5McjrJU0luXeDMAHBBm+oz59bakSRHzrrs9g0/vy/J++Y7GgCMk181BIDOiDMAdEacAaAz4gwAnRFnGJGWNvQIF4xmKVkgcYYRqqEHWGLWju0gzgDQGXEGgM6IMwB0RpwBoDPiDACdEWcA6Iw4A0BnxBlGxMYZsBzEGQA6I84wRra5mlmVxWPxxBkAOiPOANAZcQaAzogzAHRGnAGgM+IMAJ0RZwDojDjDiNggbH7stsYiiTOMUNmFZGZWju0gzgDQGXEGgM6IMwB0RpwBoDPiDACdEWcA6Iw4A0BnxBnGxM4Zc9Ns6cICiTMAdEacYYTKNlczs3ZsB3EGgM6IMwB0RpwBoDPiDACdEWcA6Iw4A0BnxBkAOiPOMCL2tJofm62xSOIMI2QfjdnZhITtIM4A0BlxBoDOiDMAdEacAaAz4gwAnRFnAOjMVHGuquur6nhVrVXVbT/kvJ+pqlZVK/MbEQDG5bxxrqodSe5MckOS/Uluqar95zjvsiTvTXLvvIcE5sPGGfNjKVmkaV45X5NkrbV2orX2XJK7ktx0jvP+WZJfS/LMHOcDgNGZJs5XJHlsw/HJ9cv+QlW9Mcne1tp/n+NswIKUba5mVvZXYxvs3OoNVNVLknwwya1TnHsgyYEk2b17dyaTyVbvfqmcOnVqdI953qzh1nzxy88nST796U/nsl0iM4uHnjydJPns/ffnm3+yY+Bplpvn8+amifPjSfZuON6zftn3XJbktUkm6/83/uokh6vqxtbafRtvqLV2KMmhJFlZWWmrq6uzT76EJpNJxvaY580abs2XP/No8sixXHvttXnlJbuGHmcptc8/kdx/NG9805vyhr0/MvQ4S83zeXPTvK19NMlVVXVlVe1KcnOSw9+7srX2rdba5a21fa21fUnuSfIDYQYApnPeOLfWTic5mOTuJI8k+Vhr7VhV3VFVNy56QAAYm6k+c26tHUly5KzLbt/k3NWtjwUA42WHMADojDgDQGfEGUak2SJsbqwliyTOMEL+wnkLLB7bQJwBoDPiDACdEWcA6Iw4A0BnxBkAOiPOANAZcQaAzogzjIhtM+bHWrJI4gwAnRFnGKGyy9XMLB3bQZwBoDPiDACdEWcA6Iw4A0BnxBkAOiPOANAZcQaAzogzjEizrdXcWEsWSZwBoDPiDCNU9rmaWdlejW0gzgDQGXEGgM6IMwB0RpwBoDPiDACdEWcA6Iw4w4jYN2OerCaLI84A0BlxhjGyj8bMLB3bQZwBoDPiDACdEWcA6Iw4A0BnxBkAOiPOANAZcQaAzogzjEhrdrWaF0vJIokzAHRGnGGEyjZXM7N2bAdxBoDOiDMAdEacAaAz4gwAnRFnAOiMOANAZ8QZYAb2IGGRpopzVV1fVceraq2qbjvH9X+vqv64qh6sqk9V1f75jwoA43DeOFfVjiR3Jrkhyf4kt5wjvh9trb2utfaGJB9I8sG5TwrMjX00ZldWj20wzSvna5KstdZOtNaeS3JXkps2ntBa+/aGw0viHR8AmNnOKc65IsljG45PJnnL2SdV1S8m+eUku5K8fS7TAcAITRPnqbTW7kxyZ1X9fJJ/muTdZ59TVQeSHEiS3bt3ZzKZzOvul8KpU6dG95jnzRpuzdqXnk+SfOpTn8rLdnp7dhbHvvZCkuSBBx7I04/uGHia5eb5vLlp4vx4kr0bjvesX7aZu5L823Nd0Vo7lORQkqysrLTV1dXpprxATCaTjO0xz5s13Jq1HSeS44/krW99ay67+KKhx1lKO7/4teS+e3P11VfnzfteOfQ4S83zeXPTfOZ8NMlVVXVlVe1KcnOSwxtPqKqrNhz+dJIvzm9EABiX875ybq2drqqDSe5OsiPJh1trx6rqjiT3tdYOJzlYVe9M8nySb+Qcb2kDANOZ6jPn1tqRJEfOuuz2DT+/d85zAcBo2SEMRqT5I8e5sZYskjgDQGfEGUaoyp9RzcrSsR3EGQA6I84A0BlxBoDOiDMAdEacAaAz4gwAnRFnAOiMOMOItNjWal6aLcJYIHGGEbKPxuysHdtBnAGgM+IMAJ0RZwDojDgDQGfEGQA6I84A0BlxBoDOiDOMiH0z5sdSskjiDACdEWcYobLN1eysHdtAnAGgM+IMAJ0RZwDojDgDQGfEGQA6I84A0BlxBoDOiDOMiF2t5sduayySOMMIlZ00Zmbt2A7iDACdEWcA6Iw4A0BnxBkAOiPOANAZcQaAzogzAHRGnGFEbJwxP82WLiyQOANAZ8QZRqhscjUza8d2EGcA6Iw4A0BnxBkAOiPOANAZcQaAzogzAHRGnAGgM+IMI2JXqzmylCyQOAO8CPYgYTtMFeequr6qjlfVWlXddo7rf7mqHq6qh6rqf1bVj81/VAAYh/PGuap2JLkzyQ1J9ie5par2n3XaA0lWWmuvT/L7ST4w70EBYCymeeV8TZK11tqJ1tpzSe5KctPGE1prf9Ra+7P1w3uS7JnvmAAwHtPE+Yokj204Prl+2Wbek+TjWxkKAMZs5zxvrKp+IclKkrdtcv2BJAeSZPfu3ZlMJvO8++6dOnVqdI953qzh1pw48VyS5JOf/GR27fCrTbM4/tQLSZIHP/e5PHdyx8DTLDfP581NE+fHk+zdcLxn/bLvU1XvTPL+JG9rrT17rhtqrR1KcihJVlZW2urq6oudd6lNJpOM7THPmzXcmoezlnzheK677rpcfJGwzOJlJ76e/J978oaf/Mn81I9fPvQ4S83zeXPTvK19NMlVVXVlVe1KcnOSwxtPqKqrk3woyY2ttSfmPyYAjMd549xaO53kYJK7kzyS5GOttWNVdUdV3bh+2q8nuTTJf6iqB6vq8CY3Bwyo2ThjbiwlizTVZ86ttSNJjpx12e0bfn7nnOcCgNGyQxiMUPldsJmVxWMbiDMAdEacAaAz4gwAnRFnAOiMOANAZ8QZADojzgDQGXEGmIHd1lgkcQaAzogzjFDFLlezskEY20GcAaAz4gwAnRFnAOiMOANAZ8QZADojzgDQGXGGEWl2zpibFmvJ4ogzAHRGnGGEbKQxO0vHdhBnAOiMOANAZ8QZADojzgDQGXEGgM6IMwB0RpwBoDPiDCNig7D5sZYskjgDQGfEGUbILlezs7sa20GcAaAz4gwAnRFnAOiMOANAZ8QZADojzgDQGXGGEbFvxvxYSxZJnAGgM+IMI1R20tgCa8fiiTMAdEacAaAz4gwAnRFnAOiMOANAZ8QZADojzgDQGXGGEWm2tZqbZjFZIHEGgM6IM4yQPa5mZ3M1toM4A0BnxBkAOjNVnKvq+qo6XlVrVXXbOa6/rqo+W1Wnq+pn5z8mAIzHeeNcVTuS3JnkhiT7k9xSVfvPOu0rSW5N8tF5DwgAY7NzinOuSbLWWjuRJFV1V5Kbkjz8vRNaa4+uX/fdBcwIAKMyzdvaVyR5bMPxyfXLAIAFmOaV89xU1YEkB5Jk9+7dmUwm23n3gzt16tToHvO8WcOt+dKjzyVJJv9rkpf4m6CZrH3zhSTJQw89lPzptv4n9ILj+by5af7NejzJ3g3He9Yve9Faa4eSHEqSlZWVtrq6OsvNLK3JZJKxPeZ5s4Zb8+DpLyRrX8zq21bzkpeI8yxe/pVvJPd8Jq97/euz+lf+0tDjLDXP581N87b20SRXVdWVVbUryc1JDi92LGCRvGienaVjO5w3zq2100kOJrk7ySNJPtZaO1ZVd1TVjUlSVW+uqpNJfi7Jh6rq2CKHBoAL2VQfmLTWjiQ5ctZlt2/4+WjOvN0NAGyRHcIAoDPiDACdEWcA6Iw4A0BnxBkAOiPOMCKtDT3BBcRaskDiDACdEWcYobJF2MysHdtBnAGgM+IMAJ0RZwDojDgDQGfEGQA6I84A0BlxBoDOiDOMiE2t5qdZTRZInAFeBFuQsB3EGQA6I84A0BlxBoDOiDMAdEacAaAz4gwAnRFnAOiMOMOYNBtnzIulZJHEGQA6I84wMna42pqygGwDcQaAzogzAHRGnAGgM+IMAJ0RZwDojDgDQGfEGQA6I84wIja1mh87hLFI4gzwIpRtXNgG4gwAnRFnAOiMOANAZ8QZADojzgDQGXEGgM6IMwB0RpxhRGycMT+WkkUSZwDojDjDyJQNrrbE+rEdxBkAOiPOANAZcQaAzogzAHRmqjhX1fVVdbyq1qrqtnNc/9Kq+r316++tqn3zHhQAxuK8ca6qHUnuTHJDkv1Jbqmq/Wed9p4k32it/XiS30jya/MeFADGYppXztckWWutnWitPZfkriQ3nXXOTUk+sv7z7yd5R5U/OACAWeyc4pwrkjy24fhkkrdsdk5r7XRVfSvJjyb52jyGPJ8/fPj/5c7J2nbc1ZZ8+1t/nn/18KeHHmOpWcOt+eo3/3zoES4Yv3rkkfybJfjvTs96fz5f+tKd+e33nJ277TFNnOemqg4kOZAku3fvzmQymcvtPvLk6Tz/9Om53NYiXVQv5PmnvzP0GEvNGm7Nqy5KfuI1bW7PvTF69nTLGy9vea49k+effmbocZZa78/nP3s2gz1Xponz40n2bjjes37Zuc45WVU7k7wiydfPvqHW2qEkh5JkZWWlra6uzjDyD1pN8vfnckuLNZlMMq/HPFbWcOus4da9dKc1nAf/Lm5ums+cjya5qqqurKpdSW5Ocviscw4neff6zz+b5BOt2WIfAGZx3lfO658hH0xyd5IdST7cWjtWVXckua+1djjJbyX57apaS/JUzgQcAJjBVJ85t9aOJDly1mW3b/j5mSQ/N9/RAGCc7BAGAJ0RZwDojDgDQGfEGQA6I84A0BlxBoDOiDMAdEacAaAz4gwAnRFnAOhMDfX9FFX1ZJIvD3Lnw7k82/Qd1xcwa7h11nDrrOF8jG0df6y19qppThwszmNUVfe11laGnmOZWcOts4ZbZw3nwzpuztvaANAZcQaAzojz9jo09AAXAGu4ddZw66zhfFjHTfjMGQA645UzAHRGnAdSVb9SVa2qLh96lmVTVb9eVZ+vqoeq6j9X1Y8MPdOyqKrrq+p4Va1V1W1Dz7NsqmpvVf1RVT1cVceq6r1Dz7SsqmpHVT1QVf9t6Fl6JM4DqKq9Sf56kq8MPcuS+oMkr22tvT7JF5K8b+B5lkJV7UhyZ5IbkuxPcktV7R92qqVzOsmvtNb2J/mrSX7RGs7svUkeGXqIXonzMH4jyT9O4gP/GbTW/kdr7fT64T1J9gw5zxK5Jslaa+1Ea+25JHcluWngmZZKa+1PW2ufXf/5OzkTlyuGnWr5VNWeJD+d5DeHnqVX4rzNquqmJI+31j439CwXiL+d5ONDD7Ekrkjy2IbjkxGWmVXVviRXJ7l32EmW0r/MmRco3x16kF7tHHqAC1FV/WGSV5/jqvcn+Sc585Y2P8QPW8PW2n9ZP+f9OfM24+9s52xQVZcm+Y9Jfqm19u2h51kmVfWuJE+01u6vqtWh5+mVOC9Aa+2d57q8ql6X5Mokn6uq5MzbsZ+tqmtaa/93G0fs3mZr+D1VdWuSdyV5R/P3gNN6PMneDcd71i/jRaiqi3ImzL/TWvtPQ8+zhK5NcmNV/Y0kFyd5eVX9+9baLww8V1f8nfOAqurRJCuttTFt/L5lVXV9kg8meVtr7cmh51kWVbUzZ36B7h05E+WjSX6+tXZs0MGWSJ35v+qPJHmqtfZLQ8+z7NZfOf+j1tq7hp6lNz5zZhn96ySXJfmDqnqwqv7d0AMtg/VfojuY5O6c+UWmjwnzi3Ztkr+Z5O3r/+49uP4KEObKK2cA6IxXzgDQGXEGgM6IMwB0RpwBoDPiDACdEWcA6Iw4A0BnxBlGoqrevP4d2BdX1SXr30f82qHnAn6QTUhgRKrqn+fMfsYvS3KytfarA48EnIM4w4hU1a6c2VP7mSQ/1Vp7YeCRgHPwtjaMy48muTRn9ia/eOBZgE145QwjUlWHk9yVM19d+prW2sGBRwLOwfc5w0hU1d9K8nxr7aNVtSPJZ6rq7a21Tww9G/D9vHIGgM74zBkAOiPOANAZcQaAzogzAHRGnAGgM+IMAJ0RZwDojDgDQGf+P/1qtNAzT0rIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110ea2f90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.plot(xpoints, uniform.prob(xpoints).eval())\n",
    "plt.xlabel('x')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### sample points from the uniform distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "True sigma value used for smearing and number of points to draw:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_sigma = 0.1\n",
    "n_points = 10000\n",
    "\n",
    "sampled_points = uniform.sample(n_points).eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "apply Gaussian smearing to each of the points (shift each point with a random number drawn from a Gaussian distribution of width defined above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampled_points += tf.random_normal((n_points,), \n",
    "                                   mean = 0, \n",
    "                                   stddev = true_sigma).eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAF3CAYAAABqlQinAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAE1xJREFUeJzt3W+MZfV93/HPN6wdEDRe26RTtKCupSBVKDSxuyW0rtQxtBEGK/AgsZy6MThI+wSrtkIVb5IHVatIxaocUqeV21WIvG5JCUpigQz9QzGjKqrsGGwMtUnqjYMLK2xqGxOvLbfd5NcHc2zNDjvMzM4d7nf2vl4SmnvOPffe3/3pLu855957psYYAQB6+oF5DwAA2JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY/vmPYAkufjii8fBgwfnPYxX1Le//e1ceOGF8x7GnmYOd84c7pw5nI1Fm8fHHnvsa2OMH97Kti1CffDgwTz66KPzHsYramVlJcvLy/Mexp5mDnfOHO6cOZyNRZvHqvryVrd16BsAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGmvx17MAztbBIw+ctvz0HTfMaSSwO+xRA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCN+XoWsGXrvwqV+DoU7DZ71ADQmFADQGNCDQCNeY8aWChOOcpeY48aABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMd+jBs4pvifNucYeNQA0JtQA0JhQA0Bj3qMGdpX3jGFn7FEDQGNCDQCNCTUANLalUFfV01X1ZFU9XlWPTuteV1UPVdUXp5+vndZXVX2oqo5X1RNV9abdfAIAcC7bzofJ3jLG+Nqa5SNJHh5j3FFVR6bl9yd5a5LLp/9+IsmHp58Am/LhMzjdTg5935jk2HT5WJKb1qz/6Fj1yST7q+qSHTwOACysGmNsvlHVnyZ5IclI8m/HGEer6ptjjP3T9ZXkhTHG/qr6eJI7xhh/MF33cJL3jzEeXXefh5McTpKlpaW/cc8998zyebV38uTJXHTRRfMexp5mDnduu3P45IkXX7LuygOv2dZtdnv79dbffrv3vxmvw9lYtHl8y1ve8tgY49BWtt3qoe+/M8Y4UVV/OclDVfVHa68cY4yq2rz4p9/maJKjSXLo0KGxvLy8nZvveSsrK1m05zxr5nDntjuHt6w7LJ0kT7/z5W+//ja7vf1662+/3fvfjNfhbJjHjW0p1GOME9PP56vqY0muSvLVqrpkjPHcdGj7+WnzE0kuW3PzS6d1wAJY/x4zsDObhrqqLkzyA2OMb02XfzLJP0tyf5Kbk9wx/bxvusn9Sd5TVfdk9UNkL44xntuNwQPzJ8ywu7ayR72U5GOrb0NnX5LfHmP8p6r6dJJ7q+rWJF9O8vZp+weTXJ/keJLvJHn3zEcNvCIWIcI+ZU53m4Z6jPGlJD92hvVfT3LtGdaPJLfNZHQAsOCcmQwAGhNqAGjMn7kEzmmL8D475zZ71ADQmFADQGMOfQOt+foUi84eNQA0JtQA0JhD3wBrONRON/aoAaAxe9TA9/nOMfRjjxoAGhNqAGjMoW/gFeXDWrA9Qg3sKd5HZ9E49A0AjdmjBubKHjK8PHvUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjfkeNcDLcMpT5s0eNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANDYlkNdVedV1Wer6uPT8huq6lNVdbyqfqeqXj2t/8Fp+fh0/cHdGToAnPu2s0f93iRPrVn+QJI7xxg/kuSFJLdO629N8sK0/s5pOwDgLGwp1FV1aZIbkvzmtFxJrknyu9Mmx5LcNF2+cVrOdP210/YAwDZtdY/615P8YpK/mJZfn+SbY4xT0/KzSQ5Mlw8keSZJputfnLYHALZp32YbVNXbkjw/xnisqpZn9cBVdTjJ4SRZWlrKysrKrO56Tzh58uTCPedZM4c7t34Ob7/y1MYbkyQvec15Hc6GedzYpqFO8uYkP1VV1yc5P8kPJfmXSfZX1b5pr/nSJCem7U8kuSzJs1W1L8lrknx9/Z2OMY4mOZokhw4dGsvLyzt8KnvLyspKFu05z5o53Ln1c3jLkQfmN5g94ul3Lp+27HU4G+ZxY5se+h5j/NIY49IxxsEk70jyiTHGO5M8kuSnp81uTnLfdPn+aTnT9Z8YY4yZjhoAFsROvkf9/iS/UFXHs/oe9F3T+ruSvH5a/wtJjuxsiACwuLZy6Pv7xhgrSVamy19KctUZtvlukp+ZwdgAYOE5MxkANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANDYtv7MJXBuefLEi7nlyAPzHgbwMuxRA0Bj9qgBtuHguiMQH7nuwjmNhEVhjxoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaGzTUFfV+VX1h1X1uar6fFX902n9G6rqU1V1vKp+p6pePa3/wWn5+HT9wd19CgBw7trKHvX/SXLNGOPHkvx4kuuq6uokH0hy5xjjR5K8kOTWaftbk7wwrb9z2g4AOAubhnqsOjktvmr6byS5JsnvTuuPJblpunzjtJzp+murqmY2YgBYIFt6j7qqzquqx5M8n+ShJH+S5JtjjFPTJs8mOTBdPpDkmSSZrn8xyetnOWgAWBT7trLRGOPPk/x4Ve1P8rEkf22nD1xVh5McTpKlpaWsrKzs9C73lJMnTy7cc541c7hzSxckt195avMN2ZDX4WyYx41tKdTfM8b4ZlU9kuRvJdlfVfumveZLk5yYNjuR5LIkz1bVviSvSfL1M9zX0SRHk+TQoUNjeXn5rJ/EXrSyspJFe86zZg537jfuvi8ffHJb/xtgnY9cd6HX4Qz497yxrXzq+4enPelU1QVJ/n6Sp5I8kuSnp81uTnLfdPn+aTnT9Z8YY4xZDhoAFsVWfpW+JMmxqjovq2G/d4zx8ar6QpJ7qupXk3w2yV3T9ncl+XdVdTzJN5K8YxfGDQALYdNQjzGeSPLGM6z/UpKrzrD+u0l+ZiajA4AF58xkANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNbRrqqrqsqh6pqi9U1eer6r3T+tdV1UNV9cXp52un9VVVH6qq41X1RFW9abefBACcq7ayR30qye1jjCuSXJ3ktqq6IsmRJA+PMS5P8vC0nCRvTXL59N/hJB+e+agBYEFsGuoxxnNjjM9Ml7+V5KkkB5LcmOTYtNmxJDdNl29M8tGx6pNJ9lfVJTMfOQAsgG29R11VB5O8McmnkiyNMZ6brvpKkqXp8oEkz6y52bPTOgBgm/ZtdcOquijJ7yV53xjjz6rq+9eNMUZVje08cFUdzuqh8SwtLWVlZWU7N9/zTp48uXDPedbM4c4tXZDcfuWpeQ9jT/M6nA3zuLEthbqqXpXVSN89xvj9afVXq+qSMcZz06Ht56f1J5Jctubml07rTjPGOJrkaJIcOnRoLC8vn90z2KNWVlayaM951szhzv3G3fflg09u+fd1zuAj113odTgD/j1vbCuf+q4kdyV5aozxa2uuuj/JzdPlm5Pct2b9u6ZPf1+d5MU1h8gBgG3Yyq/Sb07yc0merKrHp3W/nOSOJPdW1a1Jvpzk7dN1Dya5PsnxJN9J8u6ZjhgAFsimoR5j/EGS2uDqa8+w/Uhy2w7HBQDEmckAoDWhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoLF98x4A8Mo5eOSB05Zvv3JOAwG2zB41ADQm1ADQmFADQGNCDQCNCTUANCbUANCYr2cB7MCTJ17MLWu+9vb0HTfMcTSci+xRA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjW0a6qr6rap6vqr+x5p1r6uqh6rqi9PP107rq6o+VFXHq+qJqnrTbg4eAM51W9mj/kiS69atO5Lk4THG5UkenpaT5K1JLp/+O5zkw7MZJgAspk1DPcb4b0m+sW71jUmOTZePJblpzfqPjlWfTLK/qi6Z1WABYNGc7XvUS2OM56bLX0myNF0+kOSZNds9O60DAM7Cvp3ewRhjVNXY7u2q6nBWD49naWkpKysrOx3KnnLy5MmFe86zZg637/YrT522vHTBS9exPevn0Gvy7Pj3vLGzDfVXq+qSMcZz06Ht56f1J5Jctma7S6d1LzHGOJrkaJIcOnRoLC8vn+VQ9qaVlZUs2nOeNXO4fbcceeC05duvPJUPPrnj39cX2vo5fPqdy/MbzB7m3/PGzvbQ9/1Jbp4u35zkvjXr3zV9+vvqJC+uOUQOAGzTpr9KV9V/SLKc5OKqejbJP0lyR5J7q+rWJF9O8vZp8weTXJ/keJLvJHn3LowZABbGpqEeY/zsBldde4ZtR5LbdjooAGCVM5MBQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADQm1ADQmFADQGNCDQCNCTUANCbUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADS2b94DAHbPwSMPzHsIwA7ZowaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGjM96gBZuhM311/+o4b5jASzhX2qAGgMaEGgMaEGgAaE2oAaEyoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqAGhMqAGgMaEGgMaEGgAa25VQV9V1VfXHVXW8qo7sxmMAwCKYeair6rwk/zrJW5NckeRnq+qKWT8OACyCfbtwn1clOT7G+FKSVNU9SW5M8oVdeCxgjYNHHpj3EIAZ241QH0jyzJrlZ5P8xC48DsCesP4XqKfvuGFOI2Ev2o1Qb0lVHU5yeFo8WVV/PK+xzMnFSb4270HsceZwh/6ROdyxs5nD+sAuDWZvW7TX4l/d6oa7EeoTSS5bs3zptO40Y4yjSY7uwuPvCVX16Bjj0LzHsZeZw50zhztnDmfDPG5sNz71/ekkl1fVG6rq1UnekeT+XXgcADjnzXyPeoxxqqrek+Q/JzkvyW+NMT4/68cBgEWwK+9RjzEeTPLgbtz3OWRhD/vPkDncOXO4c+ZwNszjBmqMMe8xAAAbcApRAGhMqBuoqturalTVxfMey15TVf+iqv6oqp6oqo9V1f55j2mvcKrfnamqy6rqkar6QlV9vqreO+8x7VVVdV5VfbaqPj7vsXQk1HNWVZcl+ckk/2veY9mjHkryo2OMv57kfyb5pTmPZ09wqt+ZOJXk9jHGFUmuTnKbOTxr703y1LwH0ZVQz9+dSX4xiQ8LnIUxxn8ZY5yaFj+Z1e/ts7nvn+p3jPF/k3zvVL9s0RjjuTHGZ6bL38pqaA7Md1R7T1VdmuSGJL8577F0JdRzVFU3JjkxxvjcvMdyjvj5JP9x3oPYI850ql+ROUtVdTDJG5N8ar4j2ZN+Pas7K38x74F0NbdTiC6KqvqvSf7KGa76lSS/nNXD3ryMl5vDMcZ90za/ktVDkXe/kmODqrooye8led8Y48/mPZ69pKreluT5McZjVbU87/F0JdS7bIzx9860vqquTPKGJJ+rqmT1kO1nquqqMcZXXsEhtrfRHH5PVd2S5G1Jrh2+b7hVWzrVLy+vql6V1UjfPcb4/XmPZw96c5Kfqqrrk5yf5Ieq6t+PMf7hnMfViu9RN1FVTyc5NMZYpJPS71hVXZfk15L83THG/573ePaKqtqX1Q/fXZvVQH86yT9wFsGtq9XfsI8l+cYY433zHs9eN+1R/+MxxtvmPZZuvEfNXvevkvylJA9V1eNV9W/mPaC9YPoA3vdO9ftUkntFetvenOTnklwzvfYen/YMYabsUQNAY/aoAaAxoQaAxoQaABoTagBoTKgBoDGhBoDGhBoAGhNqWEBV9Tenv+F9flVdOP095R+d97iAl3LCE1hQVfWrWT2/8gVJnh1j/PM5Dwk4A6GGBVVVr87qOb6/m+RvjzH+fM5DAs7AoW9YXK9PclFWz5V+/pzHAmzAHjUsqKq6P8k9Wf1zq5eMMd4z5yEBZ+DvUcMCqqp3Jfl/Y4zfrqrzkvz3qrpmjPGJeY8NOJ09agBozHvUANCYUANAY0INAI0JNQA0JtQA0JhQA0BjQg0AjQk1ADT2/wGSAJG75zMK7QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111070450>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "plt.hist(sampled_points, bins = np.linspace(xmin, xmax, 101))\n",
    "plt.xlabel('x')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution of the uniform distribution with itself\n",
    "\n",
    "Tensorflow supports convolutions between tensors but not (continuous) functions or distributions. \n",
    "We therefore need to project our functions into a tensor, keeping the dependencies on the distribution parameters such that we can take derivatives with respect to them. The tensor represents the function values on grid defined above.\n",
    "\n",
    "We first look at the convolution of the Uniform distribution with itself e.g. to check that we got the function range correctly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expected analytic result (obtained with `SymPy`): \n",
    "\n",
    "![analytic convolution](https://gist.githubusercontent.com/andreh7/bd76330507675e6010bddea5005a7a37/raw/6a1eca9e0e5f3f37489fbbda7c20c9c08ef02ca7/uniform-analytic-convolution.png \"analytic convolution\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "define the placeholder object for our x values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = tf.placeholder(dtype = 'float32', \n",
    "\n",
    "                   # leave the exact number of entries \n",
    "                   # in the tensor open\n",
    "                   shape = (None,)\n",
    "                  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Padding of the convolution inputs: \n",
    "\n",
    "For the definition of the dimensions of the two input tensors, see e.g. [this link](https://github.com/tensorflow/tensorflow/blob/8753e2ebde6c58b56675cc19ab7ff83072824a62/tensorflow/python/ops/nn_ops.py#L2403-L2411) ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that what the convolution operator in Tensorflow actually calculates for input tensors $a$ and $b$ is (with `padding = 'SAME'`):\n",
    "\n",
    "$$\n",
    "  z_k = \\sum_{j = 0}^{|b|-1}  a_{k + j - \\left \\lfloor{(|b| - 1) / 2}\\right \\rfloor} \\cdot b_j\n",
    "$$\n",
    "\n",
    "(where the summation is over valid indices only) which is a [cross-correlation](https://en.wikipedia.org/wiki/Cross-correlation) while a [convolution](https://en.wikipedia.org/wiki/Convolution#Discrete_convolution) would be similar to (note the sign of $j$ in the second term):\n",
    "\n",
    "$$\n",
    "  z_k = \\sum_{j = 0}^{|b|-1}  a_{k + j - \\left \\lfloor{(|b| - 1) / 2}\\right \\rfloor} \\cdot b_{|b| - 1 - j}\n",
    "$$\n",
    "\n",
    "i.e. with the tensor $b$ mirrored.\n",
    "\n",
    "To reflect this, we invert the argument $x$ to the second uniform distribution in the following lines of code:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# note that we reshape the tensors into NWC (minibatch, width, channels) format.\n",
    "\n",
    "conv_values = tf.nn.conv1d(\n",
    "\n",
    "    # dimensions of this input (for NWC format) is [ batch, in_width, in_channels]\n",
    "    tf.reshape(uniform.prob(x), (1,-1,1)), \n",
    "    \n",
    "    # dimensions of the filter are [ filter_width, in_channels, out_channels ]\n",
    "    # note the minus sign of x\n",
    "    tf.reshape(uniform.prob(-x), (-1,1,1)), \n",
    "        \n",
    "    stride = 1,\n",
    "    padding = 'SAME',\n",
    "    data_format = 'NWC'\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the convolution result we get the same shape as predicted from the analytical calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAF3CAYAAABqlQinAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xlwm/l93/HPlwQP8RIPQJRW9wF6vV7vZXl3tRJg1XZSO/V0/Uev9IjbemanU2eaTNJJnXY602umaTtTN5lm0m7jTjdN2jSTY+xxnbgb21xCWq12V3t5LxOkbq4OgOJ9E/j1Dz5YU1pRBAmAz/MA79eMRsADEPjiEakPn9/v+X0fc84JAAAEU53fBQAAgLUR1AAABBhBDQBAgBHUAAAEGEENAECAEdQAAAQYQQ0AQIAR1AAABBhBDQBAgBHUAAAEWMTvAiQpGo26AwcO+F3GlpqZmVFra6vfZYQa+7B07MPSsQ/Lo9b247lz57LOuVgxzw1EUB84cECvvvqq32Vsqf7+fp08edLvMkKNfVg69mHp2IflUWv70cwuFftchr4BAAgwghoAgAAjqAEACDCCGgCAACOoAQAIMIIaAIAAI6gBAAgwghoAgAAjqAEACLCigtrMLprZj8zsDTN71dvWbWbPm1na+7vL225m9htmNmRmb5nZY5X8AAAAVLONHFH/BefcI865o979r0v6vnMuLun73n1J+qKkuPfnGUm/Va5iAQCoNaUMfT8t6Tnv9nOSvrxq+++4FS9J6jSzXSW8DwAEknNOV6byfpeBKldsUDtJ/8/MzpnZM962XufcNe/2dUm93u3dkq6s+tqr3jYAqCq/e/ay/vnpOZ0ZHvW7FFSxYq+edcI5N2JmOyQ9b2bvr37QOefMzG3kjb3Af0aSent71d/fv5EvD73p6ema+8zlxj4sHfuwNH90bl6S9L9+cE4LVxp9ribc+F5cW1FB7Zwb8f6+aWZ/IulxSTfMbJdz7po3tH3Te/qIpL2rvnyPt+3O13xW0rOSdPToUVdLlzeTau+SbpXAPiwd+3DznHP6ldPfl5TTyFKrTp487ndJocb34trWHfo2s1Yzay/clvTTkt6W9G1JX/Ge9hVJ3/Juf1vSz3lnfz8paWLVEDkAVIXBG9O6ObWgribTG1fGNTG35HdJqFLFzFH3SjplZm9KelnS/3XO/ZmkX5P0U2aWlvR5774kfVfSeUlDkv6bpH9Y9qoBwGepdEaS9Nc/1qi8k84MZ32uCNVq3aFv59x5SQ/fZfuopM/dZbuT9LWyVAcAAfXCYEZHdrTp6M682n4c0QuDWX3hQRa4oPzoTAYAGzS/lNPLF24pGY8pUmd66nCPBgYzWjlOAcqLoAaADXrl4i0tLOeV6ItKkhJ9MY2Mz+ni6KzPlaEaEdQAsEGpdFaN9XV64mC3JCkZj3rbM36WhSpFUAPABg0MZnT0QJdaGldO89nf06p93S0aGOSEMpQfQQ0AG3Bzcl7vX59SIh67bXsiHtWZ4ayWcrQURXkR1ACwAaeGVo6aE95wd0EiHtPMYk6vXRrzoyxUMYIaADYglc6qp7VRD+zquG37U0d6VF9nSqUZ/kZ5EdQAUKR83imVzupEPKq6OrvtsY7mBj26t5MTylB2BDUAFOn961PKTi98ZH66IBGP6a2RCY3NLG5xZahmBDUAFGnAO1pO3jE/XZDoi8o56TTtRFFGBDUAFCmVzuj+ne3a0dF818cf2r1dHc0RDQwy/I3yIagBoAhzizm9cmHsI2d7rxapr9OJeFSpdJZ2oigbghoAinD2wqgWc/k156cLEvGYrk3MazgzvUWVodoR1ABQhFQ6q8ZInR732oau5cSRlSNuupShXAhqAChCKp3REwe71dxQf8/n7e1u0aFoK8u0UDYENQCs4/rEvAZvTN9zfnq1RDyql87f0sJyrsKVoRYQ1ACwjsLR8Xrz0wWJeExzSzmdu0g7UZSOoAaAdaTSWcXam3T/zvainn/scI8a6k0DtBNFGRDUAHAP+bzTqaGsEvGozGz9L5DU2hTRY/u6mKdGWRDUAHAP73wwqVszi0oWOexdkOyL6Z0PJpWdXqhQZagVBDUA3EOhbejxI8WdSFZQOPHs9BDD3ygNQQ0A95BKZ/TArg7F2ps29HWfuG+7uloaWE+NkhHUALCGmYVlnbs0pkTfxo6mJam+znT8SFSpdIZ2oigJQQ0Aazh7YVRLObfh+emCZF9MN6cW9OMbU2WuDLWEoAaANQwMZtXcUKdP7e/a1NcX5qlTDH+jBAQ1AKxhpW1oz7ptQ9eya/s2xXe0fXhCGrAZBDUA3MXI+JyGMzNFtw1dSyIe08sXbml+iXai2ByCGgDuIjW4chT8mb7NzU8XJPqiWljO6+ULt8pRFmoQQQ0Ad5FKZ7Wzo1lHdrSV9DpPHuxRY30dXcqwaQQ1ANwht4m2oWvZ1livTx/sUoq+39gkghoA7vCjkQlNzC0pUeKwd0EiHtP716d0c3K+LK+H2kJQA8AdUoMZmUknNtg2dC0fLtPiqBqbQFADwB1S6awevG+7ulsby/J6H9/ZoWhbI/PU2BSCGgBWmZpf0muXx0pelrVaXZ3pxJGoUums8nnaiWJjCGoAWOWl87e0nHdKbLJt6FqSfTGNzizq3WuTZX1dVD+CGgBWSaUzamms33Tb0LUU5ruZp8ZGEdQAsMrAYEbHDvWoMVLe/x53dDTr/p3tzFNjwwhqAPBcHp3VxdHZss5Pr5bsi+nVi2OaXVyuyOujOhHUAOBJDa0c7ZZr/fSdEvGoFnN5nT1PO1EUj6AGAE9qMKvdndt0KNpakdf/9IFuNUXquJoWNoSgBgBJy7m8Tg+Xp23oWpob6vXEoR5OKMOGENQAIOnNqxOaml8u+7KsOyXjUQ3dnNYH43MVfR9UD4IaALSyLMtMOn6kp6LvU/hF4BRH1SgSQQ0AWlnf/NCeTnW2lKdt6Fr6etu0o72JeWoUjaAGUPMm5pb0xpVxJSu0LGs1M1MiHtOpoaxytBNFEQhqADXvzPBKaCYrtCzrTsm+qMZnl/T2yMSWvB/CjaAGUPMG0lm1NUX0yN7OLXm/n7QTZfgb6yOoAdQ059xK29DDPWqo35r/EnvamvTg7g4NcEIZikBQA6hpl0ZndXVsbkvmp1dLxGN67dKYphdoJ4p7I6gB1LTC8HOl10/fKRGPajnvdGZ4dEvfF+FTdFCbWb2ZvW5m3/HuHzSzs2Y2ZGb/x8wave1N3v0h7/EDlSkdAEo3kM5qb/c27e9p2dL3/dT+Lm1rqGeeGuvayBH1L0h6b9X9fyfpG865I5LGJH3V2/5VSWPe9m94zwOAwFnK5XVmeFSJeKxibUPX0hSp17HDtBPF+ooKajPbI+kvSfpt775J+qykP/Se8pykL3u3n/buy3v8c7bVPwEAUIQ3roxremF5y+enCxLxqC5kZ3Tl1qwv749wKPaI+j9J+hVJee9+j6Rx51zhLIirknZ7t3dLuiJJ3uMT3vMBIFBSgxnV15mOHfYrqFfmxTmqxr1E1nuCmX1J0k3n3DkzO1muNzazZyQ9I0m9vb3q7+8v10uHwvT0dM195nJjH5au1vfhd87N6WCH6fWzpzf9GqXsQ+ecuptNf/ziu7pv7vyma6gGtf69eC/rBrWk45L+spn9jKRmSR2Sfl1Sp5lFvKPmPZJGvOePSNor6aqZRSRtl/SR0xqdc89KelaSjh496k6ePFniRwmX/v5+1dpnLjf2YelqeR+Ozy7q4vee1z/6XFwnT/Zt+nVK3Yc/NfqWvvv2NZ1IJBXZonXcQVTL34vrWfe7wjn3q865Pc65A5L+hqQfOOf+lqQfSvor3tO+Iulb3u1ve/flPf4D5xwNbQEEyumhUeXd1i/LulOyL6ap+WW9eZV2ori7Un59+yeSfsnMhrQyB/1Nb/s3JfV4239J0tdLKxEAyi+Vzqi9OaKH92z3tY7jR3pkRjtRrK2Yoe8POef6JfV7t89Levwuz5mX9FfLUBsAVIRzTql0VscPR30fbu5sadRDezqVSmf1i5/f/BA8qlftTogAqFnnszMaGZ9Tos+fs73vlIxH9caVcU3MLfldCgKIoAZQc1KDK8PMSZ/npwsS8ZhyeaczwyzTwkcR1ABqTiqd1YGeFu3t3tq2oWt5dF+n2poiXE0Ld0VQA6gpi8t5nTk/6vvZ3qs11Nfp2OEeDQxmxCIZ3ImgBlBTXrs8ptnFnBI+tQ1dSzIe1dWxOV0apZ0obkdQA6gpA4MZRepMxw4Hq7PxT9qJskwLtyOoAdSUVDqrx/Z1qb25we9SbrO/p0V7u7cxT42PIKgB1IzR6QW9/cFE4Ia9JcnMlIjHdGZ4VEu5/PpfgJpBUAOoGaeHR+WclOgLzolkqyXjUU0vLOv1y+N+l4IAIagB1IzUYEbbtzXok7v9bRu6lmOHo6qvM+apcRuCGkBNKLQNPXFkJQyDaPu2Bj2yt5N5atyGoAZQE4ZuTuv65Hwg56dXS8SjeuvquMZnF/0uBQFBUAOoCYWj1BOBD+qYnJNODXFUjRUENYCakEpndCjWqj1dwWgbupaH92xXe3NEqUGCGisIagBVb34pp5fOjwbmIhz3Eqmv04kjUaXStBPFCoIaQNU7d2lM80t5JQNyWcv1JOIxfTAxr+HMjN+lIAAIagBVbyCdUUO96YmDwWobupbCCW8s04JEUAOoAanBrD61v0utTRG/SynK3u4WHYy2KsUyLYigBlDlMlMLevfaZKAua1mMRDyqM8OjWljO+V0KfEZQA6hqp71lTmE4kWy1RDymuaWczl0a87sU+IygBlDVBtIZdbU06BP3dfhdyoYcO9yjSJ0x/A2CGkD1+rBtaDymuoC2DV1LW1NEj+3v4oQyENQAqtePb0wpM7UQ+Laha0nGo3p7ZFKj0wt+lwIfEdQAqtbA4MrRaNjmpwsKJ8DRTrS2EdQAqlYqnVVfb5t2bm/2u5RNeXD3dnW2NGiAdqI1jaAGUJXml3I6e+FW6JZlrVZfZ7QTBUENoDq9fOGWFpfzoZ2fLkjGY7o5taDBG9N+lwKfENQAqlIqnVFjfV1o2oau5QTtRGseQQ2gKqXSWX36YJe2Ndb7XUpJ7uvcpiM72j68njZqD0ENoOrcnJzX+9enQj0/vVoiHtXZ86OaX6KdaC0iqAFUnUI3r7DPTxck4zEtLOf1ysVbfpcCHxDUAKpOKp1RtK1RH98Zrraha3niULca6+toJ1qjCGoAVSWfdzo1lFUihG1D19LSGNHRA10fNnBBbSGoAVSVd69NKju9WDXD3gWJeEzvX5/Szcl5v0vBFiOoAVSVwvDwiSPVFtQrn4d2orWHoAZQVVLpjO7f2a4dHeFsG7qWB3Z1qKe1keHvGkRQA6gas4vLevXimJJ91bEsa7W6OtOJeFSnhrLK52knWksIagBV4+yFW1rMhb9t6FqS8Ziy04t67/qk36VgCxHUAKpGajCrpkidPn2g2+9SKiLxYTtR5qlrCUENoGqk0hk9frBbzQ3hbhu6lh0dzbp/Zzt9v2sMQQ2gKlybmFP65rSSVdI2dC2JeFSvXBjT7OKy36VgixDUAKpCYTi4Gk8kWy0Rj2kxl9fZC7QTrRUENYCqMDCY0Y72JvX1tvldSkU9frBbTZE6pQaZp64VBDWA0MutahtqVh1tQ9fS3FCvxw92M09dQwhqAKH3zgcTGp9dUrKvOpdl3SkZjyl9c1rXJub8LgVbgKAGEHqF+enjVdY2dC2JPpZp1RKCGkDoDQxm9In7OhRta/K7lC3xsd52xdqbaCdaIwhqAKE2vbCs1y6PKVHly7JWMzMlvHaiOdqJVj2CGkConT0/qqWcU7JK24au5TN9MY3PLumdDyb8LgUVRlADCLVUOqttDfX61IEuv0vZUoX5eOapq9+6QW1mzWb2spm9aWbvmNm/9LYfNLOzZjZkZv/HzBq97U3e/SHv8QOV/QgAatlAOqMnD3WrKVKdbUPXEm1r0ifu62CeugYUc0S9IOmzzrmHJT0i6Qtm9qSkfyfpG865I5LGJH3Ve/5XJY1527/hPQ8Ayu7q2KzOZ2Zqan56tUQ8pnOXxjS9QDvRarZuULsV097dBu+Pk/RZSX/obX9O0pe920979+U9/jmr9g4EAHzxk7ahtTU/XZCMR7Wcd3ppeNTvUlBBRc1Rm1m9mb0h6aak5yUNSxp3zhV+jbsqabd3e7ekK5LkPT4hqaecRQOAtHK1rF3bm3U4Vt1tQ9fyqQNd2tZQT5eyKhcp5knOuZykR8ysU9KfSLq/1Dc2s2ckPSNJvb296u/vL/UlQ2V6errmPnO5sQ9LF+Z9mHdO/e/N6lO9Eb3wwgu+1eH3Pox3St9787L+wvZwn1Tm934MsqKCusA5N25mP5R0TFKnmUW8o+Y9kka8p41I2ivpqplFJG2X9JFxGefcs5KelaSjR4+6kydPbvpDhFF/f79q7TOXG/uwdGHeh69fHtPs917UX0t+Uicfvs+3Ovzeh+cjF/SvvvOuDj/0uPZ2t/hWR6n83o9BVsxZ3zHvSFpmtk3ST0l6T9IPJf0V72lfkfQt7/a3vfvyHv+Bc44V+QDKKpXOyqx22oaupTA/f2oo3EfUWFsxc9S7JP3QzN6S9Iqk551z35H0TyT9kpkNaWUO+pve878pqcfb/kuSvl7+sgHUulQ6o0/u3q7u1ka/S/HV4Vibdm1vZplWFVt36Ns595akR++y/bykx++yfV7SXy1LdQBwF1PzS3rt8rj+wWcO+V2K7wrtRP/s7etazuUVqaePVbXhXxRA6JwZHlUu75Ss0fXTd0r2xTQ5v6y3RmgnWo0IagChM5DOqLWxXo/uq622oWs5fjgqMyk1yDx1NSKoAYROKp3VscM9aozwX5gkdbU26qHd21lPXaX4LgcQKpdGZ3RpdLZm24auJRGP6fUr45qcX/K7FJQZQQ0gVAptQxM1dlnL9STiUeXyTi8O0U602hDUAEIllc5od+c2HYy2+l1KoDy6r0utjbQTrUYENYDQWM7l9eLQqJJ9UXGtn9s1Rup07HCU61NXIYIaQGi8eXVcUwvLzE+vIdkX1eVbs7o0OuN3KSgjghpAaAwMZlVn0lOHuSDf3RR+gRngqLqqENQAQiOVzuihPZ3qbKnttqFrOdDToj1d22gnWmUIagChMDG3pDeujCvZx7D3WlbaicZ0ZnhUS7m83+WgTAhqAKHw4lBWeSclWZZ1T5/pi2p6YVlvXBn3uxSUCUENIBQG0lm1N0X08N5Ov0sJtGOHo6ozKcXwd9UgqAEEnnNOA4MZHTvcowauDnVP27c16JG9nZxQVkX4jgcQeBdHZzUyPqcE89NFScRjeuvquMZnF/0uBWVAUAMIvEK3Leani5PsiyrvpNO0E60KBDWAwBsYzGpfd4v299A2tBgP7+lUe3OEdqJVgqAGEGhLubzODGe5CMcGROrrdNxrJ+qc87sclIigBhBor18e18xijrahG5Toi2pkfE7ns7QTDTuCGkCgpdIZ1deZnjpC29CNSHq/2LBMK/wIagCBNjCY0aN7O9XR3OB3KaGyt7tFB3paWKZVBQhqAIE1NrOot0YmGPbepEI70YXlnN+loAQENYDAOj2clXMr863YuGRfTHNLOb12iXaiYUZQAwis1GBWHc0RPbR7u9+lhNKTh7oVqTOWaYUcQQ0gkJxzSqUzOn4kqghtQzelvblBj+3rUop56lDjux9AIA1nZvTBxDzz0yVKxKN6+4MJjU4v+F0KNomgBhBIheFaGp2UJtEXk3PSqSGOqsOKoAYQSKl0Vgejrdrb3eJ3KaH2yd3b1dnSwPB3iBHUAAJnYTmnM8OjXISjDOrrTMePRJVKZ2gnGlIENYDAOXdpTHNLtA0tl2Q8qhuTC0rfnPa7FGwCQQ0gcFLprCJ1picP0za0HE54v/AM0E40lAhqAIGTSmf02P4utTVF/C6lKuzu3KbDsVbaiYYUQQ0gUEanF/T2yCTz02WWiMd09vyo5pdoJxo2BDWAQCksI2J+urySfVEtLOf16sUxv0vBBhHUAAIllc6qs6VBD9I2tKyePNSjhnraiYYRQQ0gMFa3Da2vM7/LqSotjREd3d/NPHUIEdQAAiN9c1o3JheYn66QRF9U712b1M2peb9LwQYQ1AACo7B8iPnpykh6+/UUR9WhQlADCIyBdFZHdrTpvs5tfpdSlR7Y1aGe1kbaiYYMQQ0gEOaXcjp7fpSLcFRQXZ3pRDyqVDqrfJ52omFBUAMIhFcvjmlhOf/h8CwqIxGPKTu9oPevT/ldCopEUAMIhFQ6o4Z60xOHuv0upaoVRixYphUeBDWAQBhIZ3V0f7daGmkbWkm9Hc36WG+7Bgjq0CCoAfju5tS83rs2qUQf89NbIRGP6pULY5pbpJ1oGBDUAHx32msbyvz01kj0xbSYy+vshVG/S0ERCGoAvksNZtXT2qgHdnX4XUpNeOJgtxojdSzTCgmCGoCvnHMaSGd1Ih5VHW1Dt0RzQ72eONjNCWUhQVAD8NX716eUnV6gG9kWS8SjGrwxresTtBMNOoIagK9+0jaUE8m2UuEXI87+Dj6CGoCvUumsPtbbrt6OZr9LqSn372xXrL2JeeoQWDeozWyvmf3QzN41s3fM7Be87d1m9ryZpb2/u7ztZma/YWZDZvaWmT1W6Q8BIJzmFnN6+eItjqZ9YGZKxKM6lc7QTjTgijmiXpb0y865ByQ9KelrZvaApK9L+r5zLi7p+959SfqipLj35xlJv1X2qgFUhZcv3tLicl6JPuan/ZCMxzQ2u6R3Ppj0uxTcw7pB7Zy75px7zbs9Jek9SbslPS3pOe9pz0n6snf7aUm/41a8JKnTzHaVvXIAoZcazKgxUqfHD9A21A/Hj6yMZDBPHWwbmqM2swOSHpV0VlKvc+6a99B1Sb3e7d2Srqz6sqveNgC4TSqd1eMHurWtsd7vUmpSrL1JD+zq+PCEPgRT0U11zaxN0h9J+kXn3KTZT9Y7OuecmW1oksPMntHK0Lh6e3vV39+/kS8Pvenp6Zr7zOXGPiydn/twbD6vH9+Y08Od86H+dwz79+GB5kV97+KS/vTPf6htEf/WsYd9P1ZSUUFtZg1aCenfc879sbf5hpntcs5d84a2b3rbRyTtXfXle7xtt3HOPSvpWUk6evSoO3ny5OY+QUj19/er1j5zubEPS+fnPvzDc1clvam/94Un9fEQdyQL+/dhw56svvvbZ9W4+wGd/Hjv+l9QIWHfj5VUzFnfJumbkt5zzv3HVQ99W9JXvNtfkfStVdt/zjv7+0lJE6uGyAFA0splFqNtTbp/Z7vfpdS0owe61NxAO9EgK+aI+rikvyPpR2b2hrftn0r6NUl/YGZflXRJ0l/zHvuupJ+RNCRpVtLfK2vFAEIvn3dKpbM62RfT6mk0bL2mSL2ePNTDCWUBtm5QO+dOSVrrJ+lzd3m+k/S1EusCUMXevTapWzOLXNYyIBLxmP71d97V1bFZ7elq8bsc3IHOZAC2XOHorbA8CP5Keg1nGP4OJoIawJZLDWb18V0d2tFO29AgOLKjTTs7mrmaVkAR1AC21Ozisl69dOvDozj47yftRLPK0U40cAhqAFvq7PlbWso5LmsZMMm+mCbnl/XW1XG/S8EdCGoAW2ognVFzQ52OHujyuxSscvxIVGbMUwcRQQ1gS6XSWT1xsEfNDbQNDZLu1kZ9cvd22okGEEENYMt8MD6noZvTXNYyoBLxqF6/Mq7J+SW/S8EqBDWALVM4qzjJZS0DKRGPKZd3OjM86ncpWIWgBrBlBtJZ9XY0Kb6jze9ScBeP7etSa2M9y7QChqAGsCVyeafTQ1kl4rQNDarGSJ2OHe7hhLKAIagBbIm3RyY0PrvE/HTAJeIxXRqd1aXRGb9LgYegBrAlCsOpJ2gbGmiFX6QGOKoODIIawJYYSGf14O4O9bQ1+V0K7uFgtFW7O7cpxTKtwCCoAVTc9MKyXrs0RjeyEDAzJfuiOjM8qqVc3u9yIIIawBZ4aXhUy3mnJEEdCsl4TFMLy3rzCu1Eg4CgBlBxqXRGLY31emx/p9+loAhPHY6qzpinDgqCGkDFDaSzevJQj5oitA0Ng+0tDXp4byftRAOCoAZQUVduzepCdoZlWSGTiMf01tVxjc8u+l1KzSOoAVRUoXkGJ5KFSzIeVd5JL9JO1HcENYCKSqUzum97sw7HWv0uBRvw8N5OtTdFaCcaAAQ1gIpZzuVpGxpSDfV1eupIjwYGs3LO+V1OTSOoAVTMWyMTmpxfVqKP+ekwSsRjGhmf04Us7UT9RFADqJjUYFZm0vHDBHUYFda9c/a3vwhqABWTSmf00O7t6mpt9LsUbMK+nhbt72nhalo+I6gBVMTk/JJevzKuZB9ne4dZIh7VmfOjWlymnahfCGoAFXFmeFS5vGNZVsgl4zHNLub02uUxv0upWQQ1gIoYGMyotbFej+6jbWiYHTvco/o6Y5mWjwhqABWRSmd17HBUDfX8NxNm7c0NemxfpwYGmaf2Cz9BAMru0uiMLt+aVZJlWVUhEY/p7Q8mNDq94HcpNYmgBlB2A7QNrSqJeFTOSadpJ+oLghpA2aUGM9rTtU0Helr8LgVl8NCeTm3f1qAU66l9QVADKKulXF5nhkdpG1pF6utMJ45ElUrTTtQPBDWAsnrzyrimFpaV5LKWVSURj+r65LyGbk77XUrNIagBlNVAOqs6k546QlBXkxPeL14vMPy95QhqAGWVSmf0yN6VOU1Ujz1dLToUa6WdqA8IagBlMzG7pDevjHO2d5VKxmM6e2FU80s5v0upKQQ1gLI5PZxV3on101Uq2RfV/FJe5y7RTnQrEdQAyiaVzqi9KaKH99A2tBo9cbBHDfWmAdqJbimCGkBZOOc0MJjVU0d6FKFtaFVqbYroU/u7aCe6xfhpAlAWF7IzGhmfY366yiXiMb13bVI3p+b9LqVmENQAyqJwNnCSoK5qhX/f00McVW8VghpAWaTSGe3vadE+2oZWtU/c16Hu1kalGP7eMgQ1gJItLq+0DeVouvrVee1EB2gnumUIagAle/3ymGYWc0rQNrQmJOJRZacX9P71Kb9LqQkENYCSDaQzqq8oLTjCAAAN7ElEQVQzHTvc43cp2AKFEwYHaCe6JQhqACVLpbN6bF+n2ptpG1oLdm5vVl9vG+1EtwhBDaAkt2YW9aORCZZl1ZhEPKaXL97S3CLtRCuNoAZQktNDWTkn5qdrTCIe1eJyXi9fvOV3KVWPoAZQklQ6o47miB6ibWhNeeJgjxojdUoxT11xBDWATXPOKZXO6kQ8qvo687scbKFtjfV6/EA3fb+3AEENYNOGM9O6NjHP/HSNSsSjGrwxresTtBOtpHWD2sz+u5ndNLO3V23rNrPnzSzt/d3lbTcz+w0zGzKzt8zssUoWD8BfhYszMD9dmwq/oKU4qq6oYo6o/4ekL9yx7euSvu+ci0v6vndfkr4oKe79eUbSb5WnTABBlEpndCjWqj1dtA2tRffvbFe0rYllWhW2blA75wYk3Xla39OSnvNuPyfpy6u2/45b8ZKkTjPbVa5iAQTHwnJOL52/RdvQGlZXZ0rGozo1lFU+TzvRStnsHHWvc+6ad/u6pF7v9m5JV1Y976q3DUCVOXdxTHNLtA2tdYm+qG7NLOrda5N+l1K1IqW+gHPOmdmGf5Uys2e0Mjyu3t5e9ff3l1pKqExPT9fcZy439mHpStmHf/DjRdWbtDTyrvpvvFfewkKk1r8PbSEvSfoff3ZWXzrcuOnXqfX9eC+bDeobZrbLOXfNG9q+6W0fkbR31fP2eNs+wjn3rKRnJeno0aPu5MmTmywlnPr7+1Vrn7nc2IelK2Uf/oe3Ujp6oF1f+Pyx8hYVMnwfSv/1/ZSu5iI6eXLz3wvsx7Vtduj725K+4t3+iqRvrdr+c97Z309Kmlg1RA6gSmSnF/TOB5NK9jE/DSkZj+rcpTHNLCz7XUpVKmZ51v+WdEbSx8zsqpl9VdKvSfopM0tL+rx3X5K+K+m8pCFJ/03SP6xI1QB8dXqIZVn4iUQ8pqWc09kLo36XUpXWHfp2zv3sGg997i7PdZK+VmpRAIJtYDCrrpYGfeK+7X6XggA4eqBLzQ11GhjM6rP3967/BdgQOpMB2JCVtqEZnYjHaBsKSVJzQ72eONhD45MKIagBbMjgjWndnFpg2Bu3ScSjGs7MaGR8zu9Sqg5BDWBDCkdNBDVWK5xYyNW0yo+gBrAhLwxmFN/Rpl3bt/ldCgIkvqNNOzuaaSdaAQQ1gKLNL+X08oVbXC0LH2FmSnjtRHO0Ey0rghpA0V65eEsLy3kl+hj2xkcl+mKamFvSj0Ym/C6lqhDUAIqWSmfVWF+nJw52+10KAujEkajMpAHmqcuKoAZQtIHBjI4e6FJLY8mXCUAV6m5t1IP3bWeZVpkR1ACKcnNyXu9fn2J+GveUiEf12uVxTc0v+V1K1SCoARTllNc2NMn8NO4hEY8pl3c6M0w70XIhqAEUJZXOKtrWqI/v7PC7FATYp/Z3qaWxnmVaZURQA1hXPu+USmd14khUdbQNxT00Rup07BDtRMuJoAawrveuTyo7vcD8NIqSiEd1cXRWl0dn/S6lKhDUANZVGMakbSiKkfDaiQ5wVF0WBDWAdaXSGd2/s107Opr9LgUhcCjaqt2d2xj+LhOCGsA9zS3m9MqFMY6mUbRCO9EXh0a1nMv7XU7oEdQA7unshVEt5vLMT2NDkn0xTS0s682r436XEnoENYB7SqWzaorU6XHahmIDnjrcozqTXhhkmVapCGoA95RKZ/T4wW41N9T7XQpCpLOlUQ/t6WSeugwIagBruj4xr8Eb00oy7I1NSMajevPKuCZmaSdaCoIawJoKR0Nc1hKbkeiLKe+kF4cZ/i4FQQ1gTQPprGLtTfpYb7vfpSCEHtnbqfamiAZoJ1oSghrAXeXzTqfSGSXiUZnRNhQb11Bfp2OHezQwmJFzzu9yQougBnBX73wwqbHZJeanUZJEX0wj43O6kJ3xu5TQIqgB3FWh/ePxI8xPY/OSXqMcrqa1eQQ1gLtKpTN6YFeHYu1NfpeCENvf06p93S0s0yoBQQ3gI2YWlnXu0hhne6MsEvGozgyPanGZdqKbQVAD+IizF0a1lHPMT6Mskn0xzSzm9PrlMb9LCSWCGsBHDAxm1dxQp6MHuvwuBVXg2OEe1dcZl73cJIIawEek0hk9eahHTRHahqJ0Hc0NenRvJyeUbRJBDeA2I+NzGs7McLUslFUiHtOPRiZ0a2bR71JCh6AGcJvU4MrwZJLrT6OMEn1ROSedHuKoeqMIagC3SaWz2tnRrCM72vwuBVXk4T2d6miOsExrEwhqAB/K5Z1ODWVpG4qyq68znYhHlUpnaSe6QQQ1gA/9aGRCE3NLSvQxP43yS8RjujYxr6Gb036XEioENYAPpQYzMpNO0DYUFVD4vuJqWhtDUAP4UCqd1YP3bVd3a6PfpaAK7e1u0aFoK/PUG0RQA5AkTc0v6bXLY0rSNhQVlIhH9dL5US0s5/wuJTQIagCSpJfO39Jy3rF+GhWV7ItpfimvcxdpJ1osghqApJVuZC2N9XpsH21DUTlPHupRQ73pBYa/i0ZQA5AkDQxmdOxQjxoj/LeAymltiuixfV1KDXJCWbH4iQSgy6Ozujg6qwTdyLAFkn0xvXttUpmpBb9LCQWCGoBSQyvDkKyfxlYo/EJIO9HiENQAlBrManfnNh2KtvpdCmrAJ+7brq6WBi57WSSCGqhxy7m8Tg/TNhRbZ6WdaIx2okUiqIEa9+bVCU3NL7MsC1sqEY8qM7Wg969P+V1K4BHUQI1LpTOqM+n4kR6/S0ENKcxT06VsfQQ1UONS6awe2tOpzhbahmLr7Nq+TfEdbUrR93tdBDVQw2aWnN64Mq4ky7Lgg0Q8prMXbml+iXai90JQAzXsvdGccnnHsiz4ItkX1eJyXi9fuOV3KYFWkaA2sy+Y2Y/NbMjMvl6J9wBQurdHc2priuiRvZ1+l4Ia9MTBHjXW12lgkHnqeyl7UJtZvaTflPRFSQ9I+lkze6Dc7wOgNM45vZ3N6djhHjXUM7iGrbetsV6fPtjFPPU6KvHT+bikIefceefcoqTfl/R0Bd4HQAkujc4qO+eYn4avEvGYfnxjSmPzeb9LCaxIBV5zt6Qrq+5flfREBd7nrv783Rv6zf6hrXq7TZucmNOvv3va7zJCjX1YmonZJUli/TR8lYhH9Wt/Kv37V+b1O+eD+/Pc1hTR//zqlkXZbSoR1EUxs2ckPSNJvb296u/vL8vrvpdZ1tLMclleq5IaLKelGRb6l4J9WJoWSSd3OV340cu6SEeyTZueni7b/1+1KO+cTu6N6PrUUqB/nmcX5Nu/s5W7fZuZHZP0L5xzf9G7/6uS5Jz7t2t9zdGjR92rr75a1jqCrr+/XydPnvS7jFBjH5aOfVg69mF51Np+NLNzzrmjxTy3EnPUr0iKm9lBM2uU9DckfbsC7wMAQNUr+9C3c27ZzH5e0vck1Uv67865d8r9PgAA1IKKzFE7574r6buVeG0AAGoJiycBAAgwghoAgAAjqAEACDCCGgCAACOoAQAIMIIaAIAAI6gBAAgwghoAgAAjqAEACDCCGgCAACv71bM2VYRZRtIlv+vYYlFJWb+LCDn2YenYh6VjH5ZHre3H/c65oi4GH4igrkVm9mqxlzjD3bEPS8c+LB37sDzYj2tj6BsAgAAjqAEACDCC2j/P+l1AFWAflo59WDr2YXmwH9fAHDUAAAHGETUAAAFGUAeAmf2ymTkzi/pdS9iY2X8ws/fN7C0z+xMz6/S7prAwsy+Y2Y/NbMjMvu53PWFjZnvN7Idm9q6ZvWNmv+B3TWFlZvVm9rqZfcfvWoKIoPaZme2V9NOSLvtdS0g9L+lB59xDkgYl/arP9YSCmdVL+k1JX5T0gKSfNbMH/K0qdJYl/bJz7gFJT0r6Gvtw035B0nt+FxFUBLX/viHpVyRxssAmOOf+n3Nu2bv7kqQ9ftYTIo9LGnLOnXfOLUr6fUlP+1xTqDjnrjnnXvNuT2klaHb7W1X4mNkeSX9J0m/7XUtQEdQ+MrOnJY045970u5Yq8fcl/anfRYTEbklXVt2/KkJm08zsgKRHJZ31t5JQ+k9aOVjJ+11IUEX8LqDamdmfS9p5l4f+maR/qpVhb9zDvfahc+5b3nP+mVaGIn9vK2sDzKxN0h9J+kXn3KTf9YSJmX1J0k3n3DkzO+l3PUFFUFeYc+7zd9tuZp+UdFDSm2YmrQzZvmZmjzvnrm9hiYG31j4sMLO/K+lLkj7nWG9YrBFJe1fd3+NtwwaYWYNWQvr3nHN/7Hc9IXRc0l82s5+R1Cypw8x+1zn3t32uK1BYRx0QZnZR0lHnXC01pS+ZmX1B0n+U9BnnXMbvesLCzCJaOfnuc1oJ6Fck/U3n3Du+FhYitvIb9nOSbjnnftHvesLOO6L+x865L/ldS9AwR42w+8+S2iU9b2ZvmNl/8bugMPBOwPt5Sd/TyklQf0BIb9hxSX9H0me97703vCNDoKw4ogYAIMA4ogYAIMAIagAAAoygBgAgwAhqAAACjKAGACDACGoAAAKMoAYAIMAIaqAGmdmnvWt4N5tZq3c95Qf9rgvAR9HwBKhRZvZvtNJfeZukq865f+tzSQDugqAGapSZNWqlx/e8pKecczmfSwJwFwx9A7WrR1KbVnqlN/tcC4A1cEQN1Cgz+7ak39fK5VZ3Oed+3ueSANwF16MGapCZ/ZykJefc/zKzekkvmtlnnXM/8Ls2ALfjiBoAgABjjhoAgAAjqAEACDCCGgCAACOoAQAIMIIaAIAAI6gBAAgwghoAgAAjqAEACLD/D/oF2xzKdXVUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1113c0e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "session.run(tf.global_variables_initializer())\n",
    "\n",
    "plt.figure(figsize=(8,6))\n",
    "\n",
    "plt.plot(xpoints, \n",
    "         session.run(conv_values,\n",
    "            feed_dict = {x : xpoints,\n",
    "                    # sigma: 0.3,\n",
    "                        }\n",
    "                    ).reshape(-1))\n",
    "plt.xlabel('x')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution of uniform distribution with the Gaussian distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "build the Gaussian distribution with a yet unknown parameter sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = tf.Variable(0.2, dtype = 'float32')\n",
    "\n",
    "gauss = tf.distributions.Normal(loc=0., scale = sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution of uniform and normal distribution\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "session.run(tf.global_variables_initializer())\n",
    "\n",
    "# we need NWC data format\n",
    "\n",
    "conv_values = tf.nn.conv1d(\n",
    "    tf.reshape(uniform.prob(x), (1,-1,1)), \n",
    "             \n",
    "    # indices of the filter are (delta_i, input_filter, output_filter)\n",
    "    tf.reshape(gauss.prob(-x), (-1,1,1)), \n",
    "\n",
    "    stride = 1,\n",
    "    padding = 'SAME',\n",
    "    data_format = 'NWC'\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot the convolution result for the initial value of the sigma parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAF3CAYAAABqlQinAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xlwm/d95/HPFxdJ8AJJUdRpS5ZcO7ITOwnjOI3b0nGPOE3jbKfJNNsmTuuOZztOm27Tadx2Z3Z3ptNru02attutNknX6SZ10iO1N5ukcWyziXP4ik/5kizbuiiT4iUSIAkQ+O0feEBJFikCxAM8AJ73a8Yj4MEj8ovHpD74Hc/vZ845AQCAxhQJugAAALA2ghoAgAZGUAMA0MAIagAAGhhBDQBAAyOoAQBoYAQ1AAANjKAGAKCBEdQAADQwghoAgAYWC7oASdq0aZPbtWtX0GXUVTqdVmdnZ9BlNDWuYfW4htXjGvojbNfx0UcfPeWcGyzn3IYI6l27dumRRx4Juoy6Gh0d1cjISNBlNDWuYfW4htXjGvojbNfRzF4p91y6vgEAaGAENQAADYygBgCggRHUAAA0MIIaAIAGRlADANDACGoAABoYQQ0AQAMjqAEAaGBlBbWZvWxmT5nZ42b2iHes38zuMbOD3p993nEzs0+Z2SEze9LM3lTLNwAAQCurpEV9vXPuaufcsPf8dkn3OuculXSv91ySbpR0qfffrZL+2q9iAQAIm2rW+r5J0oj3+A5Jo5I+7h3/nHPOSfq+maXMbKtzbqyaQgE0p3zB6dT8kuaXlrWQzWshl9dy3snJyTnJOangnJwk5/1pksxMEZMiZivPzXseMcnMOybvHDvzZ2cipoGuhLraYjKzQN8/UC0r5uk6J5m9JGlakpP0N865/WY245xLea+bpGnnXMrMviLpj5xzD3iv3Svp4865R17zNW9VscWtoaGhN995551+vq+GNz8/r66urqDLaGpcw+rV4hqmc04Pn1zWM5N5vTRb0NSiU379f2ZqIh6RtnVFtLsnoqs3R3XlpqhiEX+Dm59Df4TtOl5//fWPntVDfUHltqivc84dN7PNku4xs+fOftE558ysol9F59x+SfslaXh42IVp1xQpfDvF1ALXsHp+XsOFbF5/df8hffqBw1rMFbS1t13X7O3TxQNJbU11qLstpo5EVMlEVNGIrbSUI5FSi/lMC3mlde2kQqnF7c60uAvOqeDOnOPkVCicOZ5eWtZkeknjp5f0/KtzevTojEaPLWl7qkO//c4f0nuu2uZbS5ufQ39wHddWVlA75457f46b2ZclXSPp1VKXtpltlTTunX5c0s6z/voO7xiAFnViZkG/cscjembstH7mqm269Ucu0ZXbexqm2zm7XND9z4/rL+47qI/e+bgeOHhKf/Czr1c8yo0vaHzr/pSaWaeZdZceS/pJSU9LulvSzd5pN0u6y3t8t6QPebO/r5U0y/g00Lqm0ln9wqcf1JGpjD774WH9xQfeqNfv6G2YkJakRCyin7pii+667Tr92jv26h8ePabf+ocnVCgE1CcPVKCcFvWQpC97v3QxSV9wzn3dzB6W9CUzu0XSK5Le753/VUnvknRIUkbSL/leNYCG4JzTr//9Yzoxs6DP/8pbNbyrP+iSLigaMX3sJy9TWyyiP/3GC3rDjpRuuW530GUBF7RuUDvnDku6apXjk5JuWOW4k3SbL9UBaGh//9BRPXDolH7/vVc2fEif7bbr9+rxo7P64689pxsu36xdmzqDLglYEwM0ADYkvbSs//6N5/XW3f36hbdeFHQ5FTEz/cHPXqlY1PSHX3s26HKACyKoAWzI5773iibTWX38xssbajy6XJu72/UffmyP/vXAq3rmxOmgywHWRFADqFguX9Bnv/OSfvSHBvWmi/qCLmfDbn7bLnXEo/rsd14KuhRgTQQ1gIrd99y4JuaW9MFrLw66lKr0JuN63/AO3f34Cc1kskGXA6yKoAZQsS8+fFRDPW26/rLBoEup2vuHdyqbL+irT50MuhRgVQQ1gIqcXszp2wcndNPV2xVrgQVDrtjWo0sGO/Uvj7MuExpT8/+WAair0ecnlMs7/dQVQ0GX4gsz03uu2qaHXprSqfmloMsBzkNQA6jINw6c1KauNl29s3knkb3WOy7fLEn61gsTAVcCnI+gBlC25XxB//b8hG64fLOiPu9CFaQrt/VqU1dCo88T1Gg8BDWAsj0zdlpzS8v64b0DQZfiq0jE9KM/NKhvHZxg/W80HIIaQNkePDwlSbr2ktYKakl6+55Nmsnk9ML4XNClAOcgqAGU7cGXJrV7U6eGetqDLsV3w7uKY+6PvDwdcCXAuQhqAGXJF5wefGlKb93dPJtvVOKi/qQ2dbXp0VcIajQWghpAWQ5PzGtucVlvvrh1Znufzcz0ll19evjlqaBLAc5BUAMoy1PHZyVJb9iRCriS2nnzxX06Nr2giTnup0bjIKgBlOWp47Nqj0e0Z7B1926+YluvJOnAidmAKwHOIKgBlOXp47Pat7WnJZYNXcu+bT2SpANse4kG0rq/cQB8Uyg4HThxWq/f3ht0KTXV2xHXzv4O9qdGQyGoAazr8Km0Mtm8rmzxoJakK7b20vWNhkJQA1jX8yeLi4C8bmtPwJXU3r5tPXp5MqO5xVzQpQCSCGoAZTg0Pi8zae/mrqBLqbnSh5GD4/MBVwIUEdQA1nVwfE47+5Jqj0eDLqXmSh9GXiSo0SAIagDrOjQ+H4rWtCTt7OtQIhrRoQmCGo2BoAZwQfmC0+FTaV0akqCORSPatSmpF8fTQZcCSCKoAazj6FRG2eWC9oQkqCVpz2CXDtOiRoMgqAFcUGlSVVha1FJxnPoV7wMKEDSCGsAFvei1LMPWos4XnF6ZpPsbwSOoAVzQK5MZDXQm1NMeD7qUurnEW8/8xQmCGsEjqAFc0NGpjHb2J4Muo64u7i8G9bHpTMCVAAQ1gHUcmcroopAFdW8yrp72mI5MEdQIHkENYE25fEHHZxZCF9SSdNFAkqBGQyCoAaxpbGZR+YILZ1D3J3VkkqBG8AhqAGt6Zao4meqigfAF9c7+pI5NLyhfcEGXgpAjqAGsqdT1G9YWdTZf0KunF4MuBSFHUANY05GpjBLRiIZ62oMupe529hU/nDBOjaAR1ADWdGQyox19HYpGLOhS6q7Ui0BQI2gENYA1HZte0I4QdntL0rZUhyImHSOoETCCGsCaxmYXtK03fN3ekpSIRTTY3aYTs4xRI1gENYBVLS3ndWo+qy0hDWpJ2trbobHZhaDLQMgR1ABW9erskiRpW29HwJUEZ1uqXWO0qBEwghrAqk54LcmtqZC3qGcW5Rz3UiM4BDWAVZ30WpJbQ9yi3trbroVcXrMLuaBLQYgR1ABWtdKiDvkYtSS6vxEoghrAqsZmFtXTHlNnWyzoUgJT6vZnQhmCRFADWNXY7KK2pcLb7S2dmUh3YoYWNYJDUANY1djsQqi7vSVpsLtNsYitjNcDQSCoAaxqbHZRW0I8kUySohHTUE/7yng9EASCGsB5FnN5TaWzoV2V7Gxbe9s1Rtc3AkRQAzhPqas3zKuSlQz1tuvVOYIawSGoAZxnfK64KhlBLW3ubtPE6aWgy0CIlR3UZhY1s8fM7Cve891m9qCZHTKzL5pZwjve5j0/5L2+qzalA6iVca8FOdjdFnAlwRvsbtPc0rIWsvmgS0FIVdKi/qikZ896/seSPuGc2ytpWtIt3vFbJE17xz/hnQegiUx4LerN3bSoS9dgnO5vBKSsoDazHZJ+WtKnvecm6R2S/tE75Q5J7/Ue3+Q9l/f6Dd75AJrExNySYhFTqiMedCmB2+z1KpQ+vAD1Vm6L+pOSfltSwXs+IGnGObfsPT8mabv3eLuko5LkvT7rnQ+gSUzMLWlTV5siET5jl7r/xwlqBGTdtQHN7N2Sxp1zj5rZiF/f2MxulXSrJA0NDWl0dNSvL90U5ufnQ/ee/cY1rN5a1/C5VxbVLsf1lXQ6W9w56zs/eFrJyefPe52fQ39wHddWziK+b5f0HjN7l6R2ST2S/lxSysxiXqt5h6Tj3vnHJe2UdMzMYpJ6JU2+9os65/ZL2i9Jw8PDbmRkpMq30lxGR0cVtvfsN65h9da6hn/yxLe1Z1O7RkbeUv+iGkyh4PSbo19TastOjYxcft7r/Bz6g+u4tnW7vp1zv+Oc2+Gc2yXp5yXd55z7BUn3S/o577SbJd3lPb7bey7v9fscm7kCTWVifokZ355IxLSpq03j3KKFgFRzH/XHJf2mmR1ScQz6M97xz0ga8I7/pqTbqysRQD3lC06TBPU5Nve0MUaNwFS0f51zblTSqPf4sKRrVjlnUdL7fKgNQAAm00sqOO6hPtvm7jYdZxlRBISVyQCc48w91AR1yWB3G7dnITAENYBzlAKJFvUZg93tmkwvaTlfWP9kwGcENYBzlMZiB7tYlaxkc3ebnJMm09mgS0EIEdQAzkGL+nwri54w8xsBIKgBnGNibkndbTF1JKJBl9IwNnUVg3oyTVCj/ghqAOfgHurzbepKSJIm5+n6Rv0R1ADOUVrnG2f0d3pBTYsaASCoAZxjcn5JA14LEkVdbTElYhEmkyEQBDWAc0ylsystSBSZmQY6E3R9IxAENYAV+YLTzEJOAwT1eQa6EpqiRY0AENQAVkxnsnJOtKhX0d/Zpsl5xqhRfwQ1gBWlFmM/k8nOs6kzwRg1AkFQA1hRGoOl6/t8/YxRIyAENYAVKy1qgvo8A11tWsjllckuB10KQoagBrBiyrtPmBb1+UrXhFY16o2gBrCiNAbbR1Cfp3RvOePUqDeCGsCKqXRWPe0xxaP80/BapeGAKVYnQ53x2whgxWQ6qwFmfK+qtKzqKbq+UWcENYAVU/OsSraWMy1qghr1RVADWMHyoWtLJqJqj0dY9AR1R1ADWDGZzjLjew3F9b7bmEyGuiOoAUiSCgWn6UyWnbMuYKCLRU9QfwQ1AEnS6cWc8gWn/k4mk61loDPBntSoO4IagKQz9wfT9b22vs6EptO5oMtAyBDUACSxfGg5+pIJTWfo+kZ9EdQAJJ1ZGpOgXlt/Z0KZbF6LuXzQpSBECGoAks60qJlMtrZUMi5JmsnQ/Y36IagBSDqzNCYt6rX1JYvXhu5v1BNBDUBScTJZV1tMbbFo0KU0rFKLmqBGPRHUACQVu3NLQYTVlVrUdH2jnghqAJKKrcRSEGF1rPeNIBDUACTRoi7HmclkBDXqh6AGIKkYPila1BfUFosqmYhqmq5v1BFBDUCSNLOQUx8t6nWx6AnqjaAGoHzBaXYhp1QHQb2evs44k8lQVwQ1AM0t5uSc6PouQ18ywWQy1BVBDWBlzJXJZOtLJRNMJkNdEdQAVoKH27PW15eMM5kMdUVQA1gZc+2lRb2uVDKh04s5LecLQZeCkCCoAazMYqZFvb7+ZFzOSbMLtKpRHwQ1gJUWNbO+19fXWdqYg6BGfRDUADSTycpM6iGo15VaWe+bCWWoD4IagGYWcurtiCsasaBLaXh9Kzto0aJGfRDUADSdYbGTcrEnNeqNoAbAOt8VWBmjZtET1AlBDYCdsyrQmYgqHjW6vlE3BDUAzSywF3W5zIzVyVBXBDUAzaSLk8lQnuLqZAQ16oOgBkIuly9obmmZFnUF+pIJTafp+kZ9ENRAyJVW2GKMunx9yYRmFmhRoz7WDWozazezh8zsCTM7YGb/1Tu+28weNLNDZvZFM0t4x9u854e813fV9i0AqEZprJWgLl9vR5wlRFE35bSolyS9wzl3laSrJb3TzK6V9MeSPuGc2ytpWtIt3vm3SJr2jn/COw9AgyotH0rXd/lSyfjKdQNqbd2gdkXz3tO495+T9A5J/+gdv0PSe73HN3nP5b1+g5mx3BHQoNiLunK9ybiWlgtazOWDLgUhUNYYtZlFzexxSeOS7pH0oqQZ59yyd8oxSdu9x9slHZUk7/VZSQN+Fg3AP+xFXblUR2m9b1rVqL1YOSc55/KSrjazlKQvS7q82m9sZrdKulWShoaGNDo6Wu2XbCrz8/Ohe89+4xpWb35+Xo+89Jwk6alHH9SLcTq/ynHsZLGN8s1vfVcpy/Bz6AN+n9dWVlCXOOdmzOx+SW+TlDKzmNdq3iHpuHfacUk7JR0zs5ikXkmTq3yt/ZL2S9Lw8LAbGRnZ8JtoRqOjowrbe/Yb17B6o6Oj2hTfoujBw7rxx0fEKFV54odO6X88/qAuveIqLRx5ip9DH/D7vLZyZn0Pei1pmVmHpJ+Q9Kyk+yX9nHfazZLu8h7f7T2X9/p9zjnnZ9EA/FPakIOQLl9pcZgZZn6jDsppUW+VdIeZRVUM9i85575iZs9IutPMfl/SY5I+453/GUl/Z2aHJE1J+vka1A3AJ8UNOZhIVonS9ZrN5LQ54FrQ+tYNaufck5LeuMrxw5KuWeX4oqT3+VIdgJorbsjBRLJKnGlRZwlq1BwrkwEhN53JqY8WdUW62mKKRoxFT1AXBDUQcrPsRV0xM1Oqg0VPUB8ENRBypclkqExvMs5kMtQFQQ2EWDbvtJDLq6+TFnWlejvimqVFjTogqIEQS+eKd06yF3XlUmzMgTohqIEQK22pzPKhlUux1SXqhKAGQowW9cb1MpkMdUJQAyE27wU1C55UrrcjrrnFZRVYeBE1RlADIZYmqDesdM1oVKPWCGogxEpj1NxHXblSUJd6JYBaIaiBEEvnnKIRU2ciGnQpTae0J3WaoEaNEdRAiKVzjp2zNqjHm4BHUKPWCGogxNI5p17Gpzek1PWdZowaNUZQAyFWalGjcila1KgTghoIsfkc91BvVC9BjTohqIEQy+QcM743KBaNqKstxqxv1BxBDYTYfM7Roq5Cb0ec+6hRcwQ1EFLL+YIWllnspBqpZJwWNWqOoAZC6vTisiTGqKuRSsYZo0bNEdRASJW2aKRFvXGpjgRBjZojqIGQmskUt2gsrbCFyvV0xLmPGjVHUAMhNeO1qFnwZONKXd+OHbRQQwQ1EFKz3nRlxqg3LtURV95JmWw+6FLQwghqIKRWxqgJ6g0rje+XeieAWiCogZCaoUVdtdK1m+VmatQQQQ2E1MxCVh2x4gpb2JhebyLezEI24ErQyvgNBUJqNpNTMsb2ltUodX3TokYtEdRASM0u5NSVIKirwRg16oGgBkJqZiGnToanq7IyRk1Qo4YIaiCkZjJZdcZpUVejIx5VzM5MzANqgaAGQmp2IadOxqirYmbqTJhmmUyGGiKogRByzhWDmhZ11TrjtKhRWwQ1EEKZbF65vFMny3xXrTNmBDVqiqAGQqg0S5kWdfU648ZkMtQUQQ2EUGnnLMaoq0dQo9YIaiCEZmlR+6YrfuaDD1ALBDUQQqWVtFjwpHrJuCmdzSuXLwRdCloUQQ2E0Jkx6oALaQGlXgm6v1ErBDUQQqVZynR9V6/Lu4Z0f6NWCGoghGYXckpEI0rwL0DVSr0S3KKFWuHXFAih2YWsepNxmdGirlZpnJ+gRq0Q1EAIzWRySnUwQO2H0vABO2ihVghqIIRmMrmVLRpRnU7GqFFjBDUQQrMLuZUtGlGdZEyKRlhGFLVDUAMhVAxqFvr2g5mptyOuaVrUqBGCGgihmUyWrm8fpZJxxqhRMwQ1EDLZ5YLS2TyTyXyU6oivrPYG+I2gBkKmtIJWLy1q36SSCbq+UTMENRAyK0FNi9o3qWScyWSoGYIaCJnZhWLLL5VkMplfUh0Jbs9Czawb1Ga208zuN7NnzOyAmX3UO95vZveY2UHvzz7vuJnZp8zskJk9aWZvqvWbAFC+UsuPMWr/9CXjSmfzyi6zgxb8V06LelnSx5xz+yRdK+k2M9sn6XZJ9zrnLpV0r/dckm6UdKn3362S/tr3qgFsWCmo6fr2T2kGPTtooRbWDWrn3Jhz7gfe4zlJz0raLukmSXd4p90h6b3e45skfc4VfV9Sysy2+l45gA0phQm3Z/mn1xtGoPsbtVDRGLWZ7ZL0RkkPShpyzo15L52UNOQ93i7p6Fl/7Zh3DEADmFnIyUzqbieo/dLnfejhXmrUQqzcE82sS9I/SfoN59zps3fdcc45M3OVfGMzu1XFrnENDQ1pdHS0kr/e9Obn50P3nv3GNdyYAweXlIxJ3/7Wv3ENfTA/P69TzzwpSfr2gz9Q+uWy/1nFWfhZXFtZP1FmFlcxpD/vnPtn7/CrZrbVOTfmdW2Pe8ePS9p51l/f4R07h3Nuv6T9kjQ8POxGRkY29g6a1OjoqML2nv3GNdyYL598TAPzMxoZGeEa+mB0dFRXveEa/Zfv3a8dey7TyPDO9f8SzsPP4trKmfVtkj4j6Vnn3J+d9dLdkm72Ht8s6a6zjn/Im/19raTZs7rIAQRsdoGds/xWup6MUaMWymlRv13SByU9ZWaPe8d+V9IfSfqSmd0i6RVJ7/de+6qkd0k6JCkj6Zd8rRhAVWYy7Jzlt662mGLsoIUaWTeonXMPSLI1Xr5hlfOdpNuqrAtAjcxkstrR1xF0GS3FzNiYAzXDymRAyExncurvZFUyv/V2xOn6Rk0Q1ECI5AtOpxdzLB9aA33JBF3fqAmCGgiR2YWcnDtz3y/8k0rGNU1QowYIaiBESlsx9tGi9l1vR0KzdH2jBghqIERKY6jcnuW/PiaToUYIaiBEptPFIKFF7b9UMq5MNq+l5XzQpaDFENRAiEx5LWpmffuvNEFvlnFq+IygBkKEru/aKV1TJpTBbwQ1ECLTmZxiEVNXGxtH+K2PrS5RIwQ1ECIzmaxSyYTO3v0O/igty8qEMviNoAZCZDqd4x7qGmFjDtQKQQ2EyHQmy4zvGjnT9U2LGv4iqIEQmcmwxWWtJBNRxaPGZDL4jqAGQmQqk+XWrBop7qCV0OwCXd/wF0ENhIRzbmUyGWoj1RGn6xu+I6iBkEhn88rlHZPJaqi4MQctaviLoAZCYjrNhhy1lmKrS9QAQQ2ERClAmExWO3R9oxYIaiAkVra4ZDJZzfR1Juj6hu8IaiAkzuxFTYu6VvqSCS0tF5TJLgddCloIQQ2EBGPUtTfg9VZMztOqhn8IaiAkSgtxlNakhv9K96hPpQlq+IegBkJiJpNVT3tMsSi/9rXSR1CjBviNBUJiOpNjIlmNDRDUqAGCGgiJaVYlq7n+LoIa/iOogZCYybDFZa11t8UUj5omCWr4iKAGQmIqzYYctWZm6ksmNJVeCroUtBCCGggB55wm00srY6ionf7OBF3f8BVBDYRAJpvXYq6gga62oEtpeQNdBDX8RVADIVAKDrq+a6+/s42ghq8IaiAETs0Xx0w3dRHUtTbQmWAyGXxFUAMhcKZFTdd3rfUlE5pbXFZ2uRB0KWgRBDUQAqW1p5lMVnule6nZRQt+IaiBECh1xQ7Q9V1zrE4GvxHUQAhMzi+pIx5VMhELupSWx8Yc8BtBDYQAi53UT+k6M6EMfiGogRA4lc4y47tOVlrU86xOBn8Q1EAITKWXaFHXSV8yITO6vuEfghoIgan5LKuS1Uk0Ykp1xDXFrG/4hKAGWpxzTqfSWW7NqiPW+4afCGqgxaWzeWWXC9yaVUf9nYmVe9eBahHUQIub9CY1sSpZ/dCihp8IaqDFsdhJ/Q12t62srw5Ui6AGWhzLh9bfYFe7pjM51vuGLwhqoMVNpYstO2Z9189gd/FaT6ZpVaN6BDXQ4k7Roq67UlCPnyaoUT2CGmhxE3NL6mmPqT0eDbqU0CgF9cQcQY3qEdRAixufW1wJDtTH5lJQM6EMPiCogRY3fnpJm7vbgy4jVEoz7GlRww8ENdDiJuaXaFHXWVssqlQyTlDDFwQ10MKcc16LmqCut8GuNoIavlg3qM3ss2Y2bmZPn3Ws38zuMbOD3p993nEzs0+Z2SEze9LM3lTL4gFcWDqb10IuT4s6AIPdbYxRwxfltKj/t6R3vubY7ZLudc5dKule77kk3SjpUu+/WyX9tT9lAtiI8dOLkqTNPQR1vQ1206KGP9YNaufctyRNvebwTZLu8B7fIem9Zx3/nCv6vqSUmW31q1gAlSkFxWAXk8nqrdT17ZwLuhQ0uY2OUQ8558a8xyclDXmPt0s6etZ5x7xjAAIw7gU1Ler6G+xu00Iur3Q2H3QpaHKxar+Ac86ZWcUfGc3sVhW7xzU0NKTR0dFqS2kq8/PzoXvPfuMaru+7L+ckSS888YhOJOy817mG1VvrGp46Xrz2X/nmt7Slk3m76+FncW0bDepXzWyrc27M69oe944fl7TzrPN2eMfO45zbL2m/JA0PD7uRkZENltKcRkdHFbb37Deu4fq+/7XnFD94WD/9EyMyOz+ouYbVW+saRg9O6H899ZAu2Xe1rtndX//Cmgw/i2vb6Me8uyXd7D2+WdJdZx3/kDf7+1pJs2d1kQOos/G5RQ12ta0a0qgtlhGFX9ZtUZvZ30sakbTJzI5J+s+S/kjSl8zsFkmvSHq/d/pXJb1L0iFJGUm/VIOaAZRpYm5Jgz1MJAvCoLdb2avezHtgo9YNaufcB9Z46YZVznWSbqu2KAD+mJhb0s7+ZNBlhFJ/Z0KJWISgRtWY4QC0sJOnF1mVLCBmpq297RqbJahRHYIaaFEL2bxmMjltS3UEXUpobelp10mCGlUiqIEWNTa7IEna2ssYdVC29rbrhPf/AdgoghpoUaUu1629tKiDsqW3Q6+eXlShwOpk2DiCGmhRZ4KaFnVQtqXalcs7TaazQZeCJkZQAy1qbKbY5bqFoA7MFu/WOMapUQ2CGmhRY6cXNdCZUHs8GnQpoVUadhhjnBpVIKiBFjU2s0BrOmCl688tWqgGQQ20qLHZRSaSBWygM6FENEJQoyoENdCiikFNizpIkYhpqLdNJ+n6RhUIaqAFZbLLml3IaWuKoA7a1p4OWtSoCkENtKATM8Vg2EbXd+C2ptp1fIYWNTaOoAZa0NHpjCRpRx9BHbSL+pMam11ULl8IuhQ0KYIaaEFHp4pBfRE7ZwVuZ39S+YLT2Azd39gYghpoQUcmM2qLRTTIzlmB29lX/LB0xPvwBFSKoAZa0NHpjC7qT8rMgi4l9C4aKAZ1aTgCqBRBDbSgI1ML2km3d0PY0tOueNRoUWPDCGpl8w2oAAAKmUlEQVSgxTjndHQqw/h0g4hGTNtTHQQ1NoygBlrMdCan+aVlWtQNZGd/cmWCH1ApghpoMaWW205uzWoYFxHUqAJBDbSYlaCmRd0wLupPajqT0+nFXNCloAkR1ECLeWkiLTNp10Bn0KXAc7E38/vlU+mAK0EzIqiBFvPixLy2pzrUkWAf6kaxd3OXpOL/G6BSBDXQYg6Nz2vPYFfQZeAsF/V3KhoxHRonqFE5ghpoIYWC0+FTBHWjScQiunggqRfH6fpG5QhqoIWcmF3QYq6gPZsZn240ewe7dIiub2wAQQ20kBcnii22vbSoG86ezV16+VSaXbRQMYIaaCEvemOgezYT1I1m72CXlguOFcpQMYIaaCHPn5xTXzKugc5E0KXgNUozvw++Svc3KkNQAy3kwNisrtjWy65ZDejSoS5FTHp27HTQpaDJENRAi8jlC3rh5Lz2besJuhSsIpmIac9glw6cmA26FDQZghpoES9OzCubL+gKgrphXbm9V08fp0WNyhDUQIs44AXAvq0EdaO6YluPTp5e1MTcUtCloIkQ1ECLOHDitNrjEV3CrVkN64ptvZJE9zcqQlADLeKxo9O6cluvohEmkjWqK7YXezuePk5Qo3wENdACFnN5PX18Vm/e1Rd0KbiAnva49m7u0qOvTAddCpoIQQ20gCePzSqXdxq+uD/oUrCOa3b365GXp5UvuKBLQZMgqIEW8MgrU5KkN19Mi7rRvXV3v+aWlrmfGmUjqIEW8NBLU9oz2Kl+ViRreG/ZVez1eOilqYArQbMgqIEmt5jL63svTupHLh0MuhSUYVuqQzv7O/TdF08FXQqaBEENNLkHX5rS0nJBI5cR1M3i+ss264FDp7SYywddCpoAQQ00udHnx9UWi+jaSwaCLgVluuF1Q1rMFWhVoywENdDECgWnrz99Utft3aT2eDToclCmay/pV2ciqnueeTXoUtAECGqgiT388pTGZhf1nqu3BV0KKtAWi+r6yzfra0+f1NIy3d+4MIIaaGJffuy4OuJR/cS+oaBLQYXeN7xTM5mcvvnMeNCloMER1ECTmk5n9S+PH9fPXLVVyUQs6HJQoev2btK23nbd+fCRoEtBgyOogSb1hYeOaDFX0C9ftzvoUrAB0YjpF992sb598JSePDYTdDloYAQ10ISm01n9zb+9qOsvG9TlW9jWsll96G27lErG9clvHgy6FDQwghpoQn/6jec1v7Ss2298XdCloApdbTH96o/t0X3PjevrT58Muhw0KIIaaDJff/qkPv/gEd1y3W5dtqU76HJQpV++brf2be3Rf/qXp3R8ZiHoctCAahLUZvZOM3vezA6Z2e21+B5AGD1w8JQ+eudjesOOXv3WT10WdDnwQTwa0ac+cLWWlgv68Gcf0tgsYY1z+R7UZhaV9FeSbpS0T9IHzGyf398HCJOFbF6f/OYLuvlvH9LFA0n97YfforYYC5y0ir2bu7X/g8Mam13UTX/5HX3tqTE5xzaYKKrFPR3XSDrknDssSWZ2p6SbJD1Tg+8FtKTFXF4nZxd14MRpfe/wKf2/J8c0ncnpZ67apj/4d1equz0edInw2dv2DOgff/Vt+o9ffEK/+vkfaM9gp979hm26Zne/dm/q1FBPu6IRC7pMBKAWQb1d0tGznh+T9NYafJ9V3fvsq/rL+w+dd3ytD6drfmZd4y9c6DPu2t/j/Bfm5hbU9eS3K/w6a33fyj55V1Knn/VcsMoNfI+FhQV1PHx/meev9fUrfM8VNnIqvRbOSXOLOaWzZ1arao9H9OOvG9LNP7xrZYtEtKbLt/To/37k7brr8RP6wkNH9Bf3HVTB+2GJRkxdbTF1JqJKtsUUNZOZFDFTJFL808xkkqwJ8/z07IL+/JnvBF3GmrraYvq7W+oWZecIbJUEM7tV0q2SNDQ0pNHRUV++7oGJZeXSy2t804oOr3n8Ql+n3N+P7mhe0Vx67e9d4S9aTb/OBb7GWl9/I/9OVPoeliMFxWJL5X/vNWtd/YVK34Nf16Kjz9STiKu3zbSjO6Id3RHFI6eVfvlJjb5c4Rdbx/z8vG+/e2FVi2s4IOnXXiel9yb10mxBE5mCJhedFpadFpeXtZTPyUlyheKHvsKy96er/MNko4hbXrn0XNBlrCmzpMB+V2oR1Mcl7Tzr+Q7v2Dmcc/sl7Zek4eFhNzIy4ss3H5H06758pdoaHR2VX+85rLiG1eMaVo9r6A+u49pqMev7YUmXmtluM0tI+nlJd9fg+wAA0PJ8b1E755bN7COS/lVSVNJnnXMH/P4+AACEQU3GqJ1zX5X01Vp8bQAAwoSVyQAAaGAENQAADYygBgCggRHUAAA0MIIaAIAGRlADANDACGoAABoYQQ0AQAMjqAEAaGAENQAADcwq3cu4JkWYTUh6Jeg66myTpFNBF9HkuIbV4xpWj2voj7Bdx4udc4PlnNgQQR1GZvaIc2446DqaGdewelzD6nEN/cF1XBtd3wAANDCCGgCABkZQB2d/0AW0AK5h9biG1eMa+oPruAbGqAEAaGC0qAEAaGAEdQMws4+ZmTOzTUHX0mzM7L+Z2XNm9qSZfdnMUkHX1CzM7J1m9ryZHTKz24Oup9mY2U4zu9/MnjGzA2b20aBralZmFjWzx8zsK0HX0ogI6oCZ2U5JPynpSNC1NKl7JF3pnHuDpBck/U7A9TQFM4tK+itJN0raJ+kDZrYv2KqazrKkjznn9km6VtJtXMMN+6ikZ4MuolER1MH7hKTflsRkgQ1wzn3DObfsPf2+pB1B1tNErpF0yDl32DmXlXSnpJsCrqmpOOfGnHM/8B7PqRg024OtqvmY2Q5JPy3p00HX0qgI6gCZ2U2Sjjvnngi6lhbxy5K+FnQRTWK7pKNnPT8mQmbDzGyXpDdKejDYSprSJ1VsrBSCLqRRxYIuoNWZ2TclbVnlpd+T9LsqdnvjAi50DZ1zd3nn/J6KXZGfr2dtgJl1SfonSb/hnDsddD3NxMzeLWncOfeomY0EXU+jIqhrzDn346sdN7PXS9ot6Qkzk4pdtj8ws2uccyfrWGLDW+salpjZhyW9W9INjvsNy3Vc0s6znu/wjqECZhZXMaQ/75z756DraUJvl/QeM3uXpHZJPWb2f5xzvxhwXQ2F+6gbhJm9LGnYORemRemrZmbvlPRnkn7MOTcRdD3NwsxiKk6+u0HFgH5Y0r93zh0ItLAmYsVP2HdImnLO/UbQ9TQ7r0X9W865dwddS6NhjBrN7i8ldUu6x8weN7P/GXRBzcCbgPcRSf+q4iSoLxHSFXu7pA9Keof3s/e41zIEfEWLGgCABkaLGgCABkZQAwDQwAhqAAAaGEENAEADI6gBAGhgBDUAAA2MoAYAoIER1EAImdlbvD28282s09tP+cqg6wJwPhY8AULKzH5fxfWVOyQdc879YcAlAVgFQQ2ElJklVFzje1HSDzvn8gGXBGAVdH0D4TUgqUvFtdLbA64FwBpoUQMhZWZ3S7pTxe1WtzrnPhJwSQBWwX7UQAiZ2Yck5ZxzXzCzqKTvmtk7nHP3BV0bgHPRogYAoIExRg0AQAMjqAEAaGAENQAADYygBgCggRHUAAA0MIIaAIAGRlADANDACGoAABrY/wcVPOseVqUFaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111004290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "\n",
    "plt.plot(xpoints, \n",
    "         session.run(conv_values,\n",
    "            feed_dict = {x : xpoints, }\n",
    "                    ).reshape(-1))\n",
    "plt.xlabel('x')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### build a log likelihood function\n",
    "\n",
    "Note that we treat the convolution result as a piecewise constant function to evaluate the log likelihood.\n",
    "\n",
    "TODO: We should also ensure that the function is properly normalized.\n",
    "\n",
    "Since the final binning is fixed, we can bin the sampled values into an array before performing the maximum likelihood estination. We use [np.histogram](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.histogram.html) for this purpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_contents = np.histogram(sampled_points, bins = xboundaries)[0]\n",
    "\n",
    "assert(len(bin_contents) == len(xpoints))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that our PDF is zero in some places which causes problems when taking the logarithm. We threshold the PDF with some small value (e.g. $10^{-9}$) to protect against taking the logarithm of zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "define the negative log likelihood function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_contents_ph = tf.placeholder(dtype = 'float32', \n",
    "\n",
    "                   # leave the exact number of entries in \n",
    "                   # the tensor open\n",
    "                   shape = (None,))\n",
    "\n",
    "nll = - tf.reduce_sum( \n",
    "    bin_contents_ph * tf.log(\n",
    "        \n",
    "        # protect against bins with zero probability: set\n",
    "        # them to a small positive value\n",
    "        tf.maximum(1e-9,\n",
    "                   tf.reshape(\n",
    "                    conv_values, (-1,))\n",
    "                  )\n",
    "    \n",
    "    ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### plot the negative log likelihood function for different values of sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_grid = np.linspace(0.05, 0.2, 100 + 1)\n",
    "\n",
    "nll_values = [ session.run(nll,\n",
    "           feed_dict = { x : xpoints, \n",
    "                        sigma: sigma_val, \n",
    "                        bin_contents_ph : bin_contents\n",
    "                       })\n",
    "           for sigma_val in sigma_grid\n",
    "         ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAHlCAYAAAB258/mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8XXX9x/HXJ3snTdokTZM2XeneKRtt2bsgw6KyRFGW+kNEESsgFAERAQcKUgUVGQotUGaBsmlpoXvvpiNNd5M2+/v7495gKE1ym+Tec2/yfj4e95Hk3JH3edTiu+d8hznnEBEREWlJlNcBREREJDKoNIiIiEhAVBpEREQkICoNIiIiEhCVBhEREQmISoOIiIgERKVBREREAqLSICIiIgFRaRAREZGAxHgdINx07drVFRYWeh2j1SoqKkhOTvY6RrvTeUUWnVdk0XlFlmCc19y5c7c757q19DqVhoMUFhYyZ84cr2O02syZMxk3bpzXMdqdziuy6Lwii84rsgTjvMxsfSCv0+0JERERCYhKg4iIiAREpUFEREQCotIgIiIiAVFpEBERkYCoNIiIiEhAVBpEREQkICoNIiIiEhCVBhEREQmISoOIiIgERKVBREREAqLSICIiIgFRaRAREZGAqDSIiIhIQFQaREREJCAqDSIiIhIQlYYgcs6xt7KGA9V1XkcRERFpM5WGIFpdVs7w217n9SVbvY4iIiLSZioNQZSdlgBA6d5Kj5OIiIi0nUpDEKXGx5AcF82WPSoNIiIS+VQagsjMyElP0JUGERHpEFQagiw3LYGtutIgIiIdgEpDkOWmJVC6t8rrGCIiIm2m0hBkDbcn6uud11FERETaRKUhyLqnJ1Bb79heoasNIiIS2VQagiynYdrlHpUGERGJbCoNQZbrLw1bNYNCREQinEpDkOWmqzSIiEjHoNIQZF1T4omOMko17VJERCKcSkOQRUcZ2anxutIgIiIRz7PSYGbXm9kyM1tsZvc2On6zma0ys+Vmdmqj4z80s0X+1/+o0fHbzGyTmc3zP85o6bNCLUcLPImISAcQ48UvNbPxwARghHOuysyy/ccHAxOBIUAeMMPMioBBwHeBI4Bq4FUze8k5t8r/kb9zzt130O845Gc550K+T3VuWgKryspD/WtFRETalVdXGq4G7nbOVQE457b5j08AnnLOVTnn1gKr8BWFQcAs59x+51wt8A7wtRZ+R1OfFXK56Qka0yAiIhHPkysNQBFwvJlNBiqBG51znwA9gI8bva7Ef2wRMNnMsoADwBnAnEavu87MLvUf+7Fzblczn/UlZnYVcBVATk4OM2fObPMJNlaxvZp9VbW8OuNtEmKsXT/7YOXl5e2ePxzovCKLziuy6Lwii5fnFbTSYGYzgNxDPHWL//dmAkcBY4FnzKxPU5/lnFtqZvcArwMVwDyg4TbDw8AdgPN//S3w7cPJ6px7BHgEoLi42I0bN+5w3t6i3embeHbFPPqPGEvfbint+tkHmzlzJu2dPxzovCKLziuy6Lwii5fnFbTS4Jw7qannzOxq4DnnnANmm1k90BXYBBQ0emm+/xjOuceAx/zvvwvflQOcc6WNPvdR4CX/j01+Vqg1rAq5dU9l0EuDiIhIsHg1pmEqMB7AP9AxDtgOvABMNLN4M+sN9Adm+1/XMFiyJ77xDE/6f+7e6HPPw3crg+Y+K9Q+X+BJ4xpERCSCeTWmYQowxcwW4ZsNcZn/qsNiM3sGWALUAtc2mu3wX/+Yhhr/8d3+4/ea2Uh8tyfWAd8DcM4191khpaWkRUSkI/CkNDjnqoFvNfHcZGDyIY4f38TrL2nm9xzys0ItMS6a9MRYSlUaREQkgmlFyBDJ1QJPIiIS4VQaQiQnPUG3J0REJKKpNIRIblq8rjSIiEhEU2kIkdy0BLaXV1FbV+91FBERkVZRaQiR3PRE6h2UlVd5HUVERKRVVBpCJDc9HtBaDSIiErlUGkKk8aqQIiIikUilIUS0wJOIiEQ6lYYQyUyOIy46SqVBREQilkpDiJgZOenxlOr2hIiIRCiVhhDKTdMCTyIiErlUGkIoJy2B0r2acikiIpFJpSGEctMS2LLnAL4NPUVERCKLSkMI5aYnUFlTz94DtV5HEREROWwqDSGUm65plyIiErlUGkJIazWIiEgkU2kIof+tCnnA4yQiIiKHT6UhhHLSEjCDLVqrQUREIpBKQwjFxUSRnRrP5t260iAiIpFHpSHE8jIS2bxbVxpERCTyqDSEmK806EqDiIhEHpWGEOuRkcim3VrgSUREIo9KQ4jlpSdQVVvPzopqr6OIiIgcFpWGEMvLSATQuAYREYk4Kg0h1lAaNmlcg4iIRBiVhhDr8fmVBpUGERGJLCoNIZaRFEtibLRKg4iIRByVhhAzM7pnJLBZS0mLiEiEUWnwgG/apQZCiohIZFFp8EBeuhZ4EhGRyKPS4IG8jETK9lVRVVvndRQREZGAqTR4IC+jYYts3aIQEZHIodLggR5a4ElERCKQSoMH8rRWg4iIRCCVBg/kpvtuT6g0iIhIJFFp8EBCbDRdU+K1VoOIiEQUlQaP9MhI0FoNIiISUVQaPJKXobUaREQksqg0eKShNDjnvI4iIiISEJUGj+RlJLK/uo49B2q8jiIiIhIQlQaP9PAv8LRJtyhERCRCqDR4JE8LPImISIRRafCIFngSEZFIo9LgkazkOOJiolQaREQkYqg0eMTMyEtP0JgGERGJGCoNHtJaDSIiEklUGjzkKw0aCCkiIpFBpcFDeRmJlO6rpKau3usoIiIiLVJp8FCPjAScg617dLVBRETCn0qDhzTtUkREIolKg4c+Lw3aIltERCKASoOH8tJ9pWHTLpUGEREJTG29dxsdqjR4KDEumq4pcZSoNIiISADmb9zNT989wKJNezz5/Z6WBjO73syWmdliM7vXfyzLzN42s3Iz+8NBrx9jZgvNbJWZPWRm5j+eaWZvmNlK/9cu/uPmf90qM1tgZqNDf5bNy++SxMZd+72OISIiYa6+3nHbi4upqYdeWUmeZPCsNJjZeGACMMI5NwS4z/9UJTAJuPEQb3sY+C7Q3/84zX/8Z8Cbzrn+wJv+nwFOb/Taq/zvDysFmUls3KkrDSIi0ryp8zbx2YbdXFgUS2pCrCcZvLzScDVwt3OuCsA5t83/tcI59z6+8vA5M+sOpDnnPnbOOeAJ4Fz/0xOAx/3fP37Q8Secz8dAhv9zwkZBF9+qkHUe3qMSEZHwVl5Vy92vLGNEfjrH9ojxLId3vxmKgOPNbDK+gnCjc+6TZl7fAyhp9HOJ/xhAjnNui//7rUBOo/dsPMR7tjQ6hpldhe9KBDk5OcycOfOwT6a19pfVUFvveP61t+ma2PYOV15eHtL8oaLziiw6r8ii8wp//1lRzbZ9NVw1xNhfccCz8wpqaTCzGUDuIZ66xf+7M4GjgLHAM2bWx38VodWcc87MDusznHOPAI8AFBcXu3HjxrUlwmGJWbmdvy+eRY+iERzdN6vNnzdz5kxCmT9UdF6RRecVWXRe4W39jgpef+NdvjaqB985d6Sn5xXU0uCcO6mp58zsauA5f0mYbWb1QFegrIm3bALyG/2c7z8GUGpm3Z1zW/y3H7Y1ek9BE+8JC/ldfNMuN+7az9G0vTSIiEjHcuf0pcREGz89faDXUTwd0zAVGA9gZkVAHLC9qRf7bz/sNbOj/LMmLgWm+Z9+AbjM//1lBx2/1D+L4ihgT6PbGGEhLyMRMzTtUkREvuS9lWW8saSU607oR05agtdxPB3TMAWYYmaLgGrgsoZbE2a2DkgD4szsXOAU59wS4Brg70Ai8Ir/AXA3vtsbVwLrgYv8x18GzgBWAfuBK4J/WocnLiaK7mkJlOzUtEsREfmfmrp6bn9xCb2ykrjyuN5exwE8LA3OuWrgW008V9jE8TnA0EMc3wGceIjjDri2TUFDID9TazWIiMgXPfHRelZtK+evlxYTHxPtdRxAK0KGhYIuWqtBRET+Z3t5FQ+8sYKvFnXjxEHZXsf5nEpDGCjITKR0XyVVtXVeRxERkTBw32vLOVBTx6SzBuNf/DgsqDSEgYIuSTinjatERAQWluzh6TkbueLYQvplp3gd5wtUGsJAQaZvDfGNKg0iIp2ac45bX1hEVnIc15/Y3+s4X6LSEAYKMv1rNWgGhYhIpzZ13iY+3bCbm04bSJpH+0s0R6UhDOSkJhAXHaUZFCIinVh5VS2/ftm3v8QFo/NbfoMHvFynQfyiooweXRIp0QwKEZFO6/dvrmTbvioeubSYqKjwGfzYmK40hIn8Lom60iAi0kmtLitnygdruXBMPiMLMryO0ySVhjCR3yVJYxpERDoh5xy3v7iEhJhobjrN+/0lmqPSECYKMhPZtb+G8qpar6OIiEgIvbGklHdXlPGjk4volhrvdZxmqTSEiYIuvmmXJbpFISLSaVTW1HHH9CX0z07h0qN7eR2nRSoNYeLztRo0GFJEpNN45N01bNx5gNvPGUJsdPj/X3L4J+wkCrporQYRkc6kZNd+/jRzFWcMy+WYfl29jhMQlYYwkZkcR1JctGZQiIh0Ene9vBSAW84c7HGSwKk0hAkz026XIiKdxAertvPywq1cM64fPTISvY4TMJWGMFKQmaiBkCIiHVxNXT23vbCYgsxErvpKH6/jHBaVhjDSsFaDc87rKCIiEiRPfLSeldvKmXTmYBJio72Oc1hUGsJIQWYSFdV17Npf43UUEREJgrJ9VTzwxgqO79+VkwfneB3nsKk0hJF8zaAQEenQfvPaMg7U1HHr2UMwC8/9JZqj0hBGGhZ40gwKEZGO57MNu3hmTglXHtebftkpXsdpFZWGMFKQ2XClQTMoREQ6krp6xy+nLSY7NZ7rT+zvdZxWU2kII6kJsXRJimXDzgqvo4iISDt6Zs5GFm7awy1nDiIlPsbrOK2m0hBmemUls36Hbk+IiHQUu/dXc++ryziiMJNzRuR5HadNVBrCTO+uKg0iIh3Jb19fwZ4DNdw+ITIHPzam0hBmemUlsXnPASpr6ryOIiIibbR48x7+NWs9lxzVi0Hd07yO02YqDWGmMCsZ57RFtohIpHPOceu0xWQkxXHDyQO8jtMuVBrCTK8s37TLddtVGkREItlzn25izvpd/PS0AaQnxXodp12oNISZwqxkANbt0AwKEZFItbeyhl+/soyRBRlcOKbA6zjtJnLnfXRQGUmxpCXEaDCkiEgEe+CNleyoqGLK5cVERUX24MfGdKUhzJgZhV2TdaVBRCRCLdu6l8c/WsfFR/RkeH6G13HalUpDGOqVpdIgIhKJGgY/pibE8JNTOsbgx8ZUGsJQ76wkNu06QHVtvddRRETkMLwwfzOz1u7kJ6cOoEtynNdx2p1KQxjqlZVMvaZdiohElPKqWu56eSnDeqQzcWxPr+MEhUpDGCrs6pt2qcGQIiKR46E3V1K6t4pfTRhCdAca/NiYSkMY6qVplyIiEWVl6T6mvL+WiWMLGNWzi9dxgkalIQxlJceREq9plyIikcA537bXyfEx3HTaQK/jBJVKQxgyM3plJelKg4hIBHhpwRY+WrODG08dQGYHHPzYmEpDmCrMSmbddpUGEZFwVlFVy+TpSxnaI41vHNExBz82ptIQpgq7JlGy6wA1dZp2KSISrh56ayVb91byqwlDO+zgx8ZUGsJUr6xkausdm3cf8DqKiIgcwqpt+3jsvbVcVJzP6A48+LExlYYw9b+NqzQYUkQk3DjnmDTVN/jxpx188GNjKg1hqjCrYa0GjWsQEQk3L/oHP/7k1AFkpcR7HSdkVBrCVLfUeBJjo1m3XVcaRETCyb7KGu58aQnD89O5uBMMfmxMW2OHqYZpl7rSICISXh6csZKy8ioevbS4Uwx+bExXGsJYoXa7FBEJK8u37uNvH65j4tiejCjoWNteB0KlIYz16prExp0HqKt3XkcREen0nHNMmraItIQYbjq14217HQiVhjDWOyuZ6rp6TbsUEQkDz3+2idlrd/LT0wZ2yG2vA6HSEMYaNq7SHhQiIt7ac6CGu15eyqieGVxUXOB1HM+oNISxhi2yNa5BRMRb97++nJ0V1dwxYShRnWzwY2MqDWEsJzWB+Jgo7UEhIuKhRZv28I+P13PJUb0Y2iPd6zieUmkIY1FRRu+uyaxVaRAR8UR9veMXUxeRmRzPDad0zsGPjak0hLk+3ZJZo9IgIuKJp+dsZN7G3fz8jIGkJ8Z6HcdznpUGM7vezJaZ2WIzu9d/LMvM3jazcjP7w0GvH2NmC81slZk9ZGbmP36bmW0ys3n+xxmN3nOz//XLzezU0J5h++jbLYUNO/dTVVvndRQRkU5lZ0U197y6jCN6Z3LeqB5exwkLnqwIaWbjgQnACOdclZll+5+qBCYBQ/2Pxh4GvgvMAl4GTgNe8T/3O+fcfQf9jsHARGAIkAfMMLMi51xE/b9v324p1NU7NuzYT/+cVK/jiIh0Gve8sozyylruPHco/n+ndnpeXWm4GrjbOVcF4Jzb5v9a4Zx7H195+JyZdQfSnHMfO+cc8ARwbgu/YwLwlHOuyjm3FlgFHNHO5xF0fbr5pl2uLtMtChGRUJmzbidPz9nIlcf3pkj/YPucV6WhCDjezGaZ2TtmNraF1/cAShr9XOI/1uA6M1tgZlPMrEuj92xs5j0RoU+3FABWl5V7nEREpHOoravnF1MXkZeewA9O6O91nLAStNsTZjYDyD3EU7f4f28mcBQwFnjGzPr4ryIcroeBOwDn//pb4NuHmfUq4CqAnJwcZs6c2YoYwdMl3vhw4WqGWEmLry0vLw+7/O1B5xVZdF6RRef1Ra+tq2HZ1mquHxXPJx+93/7B2sjLP6+glQbn3ElNPWdmVwPP+UvCbDOrB7oCZU28ZROQ3+jnfP8xnHOljT73UeClRu8pONR7DpH1EeARgOLiYjdu3Lgmz8sLg1Z+TEV1HePGHdvia2fOnEm45W8POq/IovOKLDqv/9m6p5Jr35rJ+AHduOGisWE5lsHLPy+vbk9MBcYDmFkREAdsb+rFzrktwF4zO8o/a+JSYJr//d0bvfQ8YJH/+xeAiWYWb2a9gf7A7PY+kVDo2y2F1WXltO5CjIiIBOqOl5ZQW++4/RwNfjwUT2ZPAFOAKWa2CKgGLmu4NWFm64A0IM7MzgVOcc4tAa4B/g4k4ps10TBz4l4zG4nv9sQ64HsAzrnFZvYMsASoBa6NtJkTDfp2S2ZfZS1l5VVkpyZ4HUdEpEN6Z0UZ0xdu4YaTi+iZleR1nLDUZGkwsxuae6Nz7v7W/lLnXDXwrSaeK2zi+By+PA0T59wlzfyeycDk1qUMHw2DIdeUVag0iIgEQWVNHbdOW0Tvrsl876t9vI4Ttpq7PZHqfxTjmyLZw//4PjA6+NGkQd9szaAQEQmmh2euZt2O/dwxYSjxMdFexwlbTV5pcM7dDmBm7wKjnXP7/D/fBkwPSToBoHtaAomx0azeprUaRETa29rtFTw8czVnj8jjuP5dvY4T1gIZCJmDb9xBg2r/MQmRqCjz70GhKw0iIu3JOccvpy0iPiaKSWcO8jpO2AtkIOQT+KZFPg8YvpUW/x7MUPJlfbqlMG/jLq9jiIh0KNMXbuG9ldu57ezBZKdpzFhLWrzS4B9MeAWwC9gBXOGc+3Wwg8kX9e2WTMmuA1TWROQEEBGRsLOvsoY7XlrC0B5pXHJ0oddxIkKg6zTUAfWNHhJifbul4Jzv3puIiLTd/W+sYNu+Ku48dxjRUVqTIRAtlgYz+yHwL3wrNmYD/zSz64MdTL6oYeOqNdq4SkSkzRZt2sPjH67jW0f2YmRBhtdxIkYgYxquBI50zlUAmNk9wEfA74MZTL6oT1dNuxQRaQ919Y5bnl9IZnI8N546wOs4ESWQ2xOG7/ZEgzr/MQmhxLhoemQkqjSIiLTRk7PWM79kD5POGkR6YqzXcSJKIFca/gbMOmj2xGNBTSWH1Kdbsm5PiIi0wbZ9ldz76nKO7ZfFOSPyvI4TcQKZPXE/vtkTO/FtKnWFc+6BYAeTL9PGVSIibXPnS0upqq3njgnakKo1Dmf2hPM/NHvCI32zU9hfXcfWvZVeRxERiTjvrSzjhfmb+f64vp/v6SOHR7MnIkjfrppBISLSGpU1dUya6tuQ6ppxfb2OE7E0eyKCNN646th+Wh9dRCRQf3p7Fet27OefVx5JQqw2pGotzZ6IINmp8aTEx7B6m2ZQiIgEatW2ch5+ZzXnjtSGVG11uLMnAM5Fsyc8YWb07ZbMKk27FBEJiHO+NRkSY6O55czBXseJeC2WBufc/Wb2DnCs/9AVzrnPghtLmlKUk8rby8u8jiEiEhH+++kmZq3dyV3nDaNbarzXcSJeoLMn5gH/AaYCO8ysZ/AiSXOKclLZXl7Fzorqll8sItKJlVc77np5KWN6dWHi2AKv43QILV5p8M+UuBUo5X/jGRwwPLjR5FCKclMBWFG6j6P6ZHmcRkQkfD29vJq9B+qYfN5QorQhVbsIZEzDD4EBzrkdwQ4jLSvK8c2gWKnSICLSpI/X7OC9TbVcPa4vA3PTvI7TYQRye2IjsCfYQSQwuWkJpCbEsLx0n9dRRETCUlVtHbc8v5BuicYPTujvdZwOpckrDWZ2g//bNcBMM5sOVDU8719eWkLMzCjKSWVFqWZQiIgcyp9nrmF1WQU3jIknMU5rMrSn5m5PpPq/bvA/4vwP8VhRTiqvLNqCc05rp4uINLKmrJw/vr2Ks0fkMbybLpK3tyZLg3Pu9lAGkcAV5aTw79k1lJVXkZ2a4HUcEZGw4FuTYRHxsVFMOmsQS+Z+7HWkDqe52xMPOOd+ZGYv4pst8QXOuXOCmkyaNCDHP4Nia7lKg4iI338/3cRHa3Yw+byhZKcmsMTrQB1Qc7cn/uH/el8ogkjg+uf8b9qllkQVEYGdFdVMnr6EMb26cPFYLSUULM3dnpjr//pO6OJIILqmxNElKZaV2zSDQkQE4M7pS9hXWctd5w3TmgxB1NztiYUc4rYE/sWdnHNa3MkjDTMolm9VaRAR+WDVdp77dBPXju/LgNzUlt8grdbc7YmzQpZCDltRTipTP9ukGRQi0qlV1vjWZCjMSuJ6rckQdE0u7uScW9/w8B/q7/9+G7AzJOmkSUW5qeyrqmXLnkqvo4iIeOaPb69i3Y79TD5vGAmxWpMh2FpcEdLMvotvs6q/+A/l49u4SjxUlO1bTnqFVoYUkU5qRek+/vzOar42qgfH9tOg8FAIZBnpa/Fti70XwDm3EsgOZihpWVGjGRQiIp1Nfb3j5ucWkhwfwy1nDvI6TqcRSGmocs59vg+zmcVw6AGSEkJdkuPolhqv5aRFpFN6cvYG5q7fxS/OHExWSrzXcTqNQErDO2b2cyDRzE4GngVeDG4sCcSAnFRdaRCRTqd0byX3vLKMY/pmcf7oHl7H6VQCKQ0/A8qAhcD3gJedc7cENZUEpH9OCitLy6mv14UfEek8bn9xMdV19dx13jDNHgux5qZcNhjlnHsUeLThgJmd5Zx7KXixJBADclI5UFNHya4D9MxK8jqOiEjQvbGklJcXbuUnpw6gsGuy13E6nUCuNDxqZkMbfjCzi4FJwYskgeqvwZAi0omUV9Xyy2mLGJCTylVf6eN1nE4pkNJwAfCEmQ30T7+8BjgluLEkEEU5vmmXy1UaRKQTuO+15WzdW8mvzx9GbHQg//cl7a3F2xPOuTVmNhHf2gwbgFOccweCnkxalJoQS156gq40iEiHN3f9Lh7/aB2XHV3I6J5dvI7TaR3O3hOZQDQwy8zQ3hPhYUCu9qAQkY6turaem59bQPe0BG48dYDXcTo17T0R4QbnpfHeyu1U1dYRH6MlVEWk4/nzO6tZUVrOY5cVkxIfyPh9CZbmbgrt8u81sa+Jh4SBQd3TqK13rNQiTyLSAa3aVs4f3lrFWcO7c+KgHK/jdHrNVbYn8V1tmIvvNkXjybAO0NDVMDC4exoAS7bsZWiPdI/TiIi0H99S0QtIjIvm1rOHeB1HaKY0OOfO8n/tHbo4crh6ZSWTFBfN0i17vY4iItKunpy9gU/W7eLe84fTLVVLRYeD5gZCjm7ujc65T9s/jhyu6ChjQG4qSzarNIhIx7FlzwHufmUZx/bL4sLifK/jiF9ztyd+28xzDjihnbNIKw3unsYL8zfjnJaTFpHI55xj0tRF1NbX8+vzhmup6DDS3O2J8aEMIq03qHsa/5q1gU27tXyGiES+6Qu3MGPpNm45Y5CWyA8zWlKrAxic5x8MqVsUIhLhdu+v5rYXFjOsRzpXHFvodRw5iEpDBzAwNxUz3wwKEZFIduf0pezeX8M95w8nRktFhx39iXQASXEx9M5K1gwKEYlo76/czn/mlvC9r/b5/AqqhJcWl9ZqYhbFHmC9c662/SNJawzKS2NByW4oUA8Ukcizv7qWnz23gD5dk7n+hP5ex5EmBLIe55+A0cACfAs8DQUWA+lmdrVz7vUg5pMADe6exvQFW9hfo0FDIhJ57nttBSW7DvDs948mIVZL4oerQP5ZuhkY5Zwrds6NAUYBa4CTgXuDGU4C17Ay5MZ99R4nERE5PHPX7+JvH67l0qN7MbYw0+s40oxASkORc25xww/OuSXAQOfcmrb8YjO73syWmdliM7vXfyzLzN42s3Iz+8NBr59sZhvNrPyg4/Fm9rSZrTKzWWZW2Oi5m/3Hl5vZqW3JG+4a7v9tUGkQkQhSVVvHT//r28HyptMGeh1HWhDI7YnFZvYw8JT/568DS8wsHqhpzS81s/HABGCEc67KzLL9T1UCk/DdAhl60NteBP4ArDzo+JX4NtfqZ2YTgXuAr5vZYGAiMATIA2aYWZFzrq41mcNddmo8mclxutIgIhHlj2+tYtW2cv52xVjtYBkBArnScDmwCviR/7HGf6wGaO0CUFcDdzvnqgCcc9v8Xyucc+/jKw9f4Jz72Dm35RCfNQF43P/9f4ATzbd82ATgKedclXNurf8cjmhl3rBnZgzunsaGvSoNIhIZlm7Zy59mruZro3owfkB2y28Qz7VY65xzB8zs98Dr+JaPXu6ca7jC0Nr9mIuA481sMr6CcKNN3AwQAAAgAElEQVRz7pNWflYPYKM/a62Z7QGy/Mc/bvS6Ev+xLzGzq4CrAHJycpg5c2Yro3grpbaKkn11vPnW20RHdaxlV8vLyyP2z6U5Oq/IovNqP3X1jjs+riQxxnFCl11B+f3682p/gUy5HIfvX/Lr8M2eKDCzy5xz77bwvhlA7iGeusX/ezOBo4CxwDNm1sd5tHmCc+4R4BGA4uJiN27cOC9itNmu9BJeXTefgiHFFOWkeh2nXc2cOZNI/XNpjs4rsui82s9f3lnNur3L+MM3RnHW8Lyg/A79ebW/QG4g/RY4xTm3HMDMioB/A2Oae5Nz7qSmnjOzq4Hn/CVhtpnVA12BskCDN7IJKABKzCwGSAd2NDreIN9/rMMa3D0d8C0n3dFKg4h0HGvKyrn/jRWcMjiHM4d19zqOHIZAxjTENhQGAOfcCiC2jb93Kv7xEP4SEgdsb+VnvQBc5v/+AuAtfxl5AZjon13RG+gPzG5T6jDXp1syMYZWhhSRsFVf7/jZfxcSFxPFnecO1Q6WESaQKw1zzOyvwD/9P38TmNPG3zsFmGJmi4Bq4LKGWxNmtg5IA+LM7Fx8VzmW+KdlfgNIMrMS4K/OuduAx4B/mNkqYCe+GRM45xab2TPAEqAWuLajzpxoEBsdRY/UKBZr4yoRCVP/mrWe2et2cu/5w8lOS/A6jhymQErD1cC1wA/8P7+Hb5XIVnPOVQPfauK5wiaO3wTcdIjjlcCFTbxnMjC51UEjUK+0KBZs3oNzTg1eRMLKpt0HuPuVZRzfvysXFud7HUdaIZDZE1XA/f6HhLneaVG8W1JNya4DFGRqSWkRCQ/OOW5+biEOuOu8YfpHTYRqsjSY2UJ8UywPyTk3PCiJpE16p/uGqcwv2a3SICJh4z9zS3h3RRm3nzNE/22KYM1daTgrZCmk3eSnRhEXHcXCkj1Bm8YkInI4SvdWcsdLSziiMJNLjurldRxpgyZLg3NufSiDSPuIiTIGdU9lQcker6OIiOCc45bnF1JVW889FwwnqoMtPNfZBDLlUiLMsPx0Fm3aQ329J2tliYh87oX5m5mxdBs/OXUAvbsmex1H2kiloQMa3iODfVW1rNtR4XUUEenEyvZVcesLixnVM4Mrju3tdRxpBwGVBjNLNLMBwQ4j7WNYvm9lSN2iEBEv3frCIvZX1fGbC4Z3uP1wOqsWS4OZnQ3MA171/zzSzF4IdjBpvf7ZKSTERqk0iIhnXlqwmZcXbuVHJ/enX7aWte8oArnScBu+LaV3Azjn5gG6zhTGYqKjGJKXzsJNu72OIiKd0PbyKn45bTEj8tO56vg+XseRdhRIaahxzh38T1aNsAtzw3qks2jTXuo0GFJEQuzWaYspr6zlNxeOICZaQ+c6kkD+NBeb2TeAaDPrb2a/Bz4Mci5po+H56RyoqWN1WbnXUUSkE5m+YAvTF27hhyf11267HVAgpeF6YAhQBTwJ7AF+FMxQ0nbD/YMh52/ULQoRCY0d5VVMmraI4fnpfO8rui3REQVSGgY6525xzo31P37h3yRKwljvrikkx0WzcJMGQ4pIaPzSf1viPt2W6LAC+VP9rZktNbM7zGxo0BNJu4iOMob2SNcMChEJCd2W6BxaLA3OufHAeKAM+IuZLTSzXwQ9mbTZ8Px0lmzZS01dvddRRKQD2+6/LTFCtyU6vICuHznntjrnHgK+j2/Nhl8GNZW0i2H5GVTX1rN86z6vo4hIB/bLaYt0W6KTCGRxp0Fmdpt/q+yGmRP5QU8mbTa8h28wpMY1iEiwNF7Eqb9uS3R4gVTCKfgWdjrVOTfOOfewc25bkHNJO+iVlURaQozGNYhIUJTtq2LS1EWMKMjQIk6dRJNbYzdwzh0diiDS/syM4fkZLCjRtEsRaV/OOSZNXURFVR33XTBctyU6iSb/lM3sGf/XhWa2oNFjoZktCF1EaYuRBRks27qP/dW1XkcRkQ7khfmbeXXxVm44pUi3JTqR5q40/ND/9axQBJHgGN0rg7p6x4KSPRzVJ8vrOCLSAZTureSX0xYzumcG39VtiU6lySsNzrkt/m+vcc6tb/wArglNPGmrUQVdAJi7fpfHSUSkI3DOcfNzC6mqreO+C0doy+tOJpCbUCcf4tjp7R1EgqNLchx9uibz2QaVBhFpu2fnlvDWsm3cdOpA+nRL8TqOhFiTtyfM7Gp8VxT6HDSGIRX4INjBpP2M6tmFt5dvwzmHmf5VICKts3n3Ae54cQlH9s7k8mMKvY4jHmjuSsOTwNnAC/6vDY8xzrlvhSCbtJPRvTLYWVHN+h37vY4iIhHKOcdP/7uAOuf4zQUjiNJtiU6puTENe5xz65xzF/vHMRwAHJBiZj1DllDabHRP37iGT3WLQkRa6Z8fr+e9ldv5+RmD6JmV5HUc8UggK0KebWYrgbXAO8A64JUg55J2VJSTSkp8jEqDiLTKuu0V3PXyMr5S1I1vHql/M3ZmgQyEvBM4CljhnOsNnAh8HNRU0q6io4yRBRl8ul6LPInI4amrd/z42fnERhv3nj9c46I6uUBKQ41zbgcQZWZRzrm3geIg55J2NrpnBsu27qWiSos8iUjgHnl3DXPX7+JXE4aSm57gdRzxWCClYbeZpQDvAv8ysweBiuDGkvY2qlcX6h3M15LSIhKgZVv38rs3VnD60FwmjMzzOo6EgUBKwwR8gyD/D3gVWI1vFoVEkNH+RZ4+1SJPIhKA6tp6/u/p+aQlxnDnuUN1W0KAwDasanxV4fEgZpEgSk+KpW+3ZD7doCsNItKyB99cwdIte3n00mKyUuK9jiNhIpDZE/vMbO9Bj41m9ryZadHxCDK6Zxc+27AL55zXUUQkjM1dv5OHZ67mouJ8Th6c43UcCSOB3J54APgJ0APIB27Et/DTU8CU4EWT9jamVxd27a9h7XYNSRGRQ6uoquWGZ+aTl5HIpLMGex1HwkwgpeEc59xfnHP7nHN7nXOPAKc6554GugQ5n7Sj0b0aFnnSLQoRObS7Xl7Khp37ue/CEaQmxHodR8JMIKVhv5ldZGZR/sdFQKX/OV3njiD9uqWQmqBFnkTk0N5evo1/zdrAd47rzVF9sryOI2EokNLwTeASYBtQ6v/+W2aWCFwXxGzSzqI+X+RJpUFEvmhXRTU//c8C+men8ONTBngdR8JUILMn1tD0FMv32zeOBNvYwkx+N2MFe/bXkJ6kS48i4tuM6papC9m1v5opl48lITba60gSpgKZPVFkZm+a2SL/z8PN7BfBjybBcETvTJyDT9bt9DqKiISJDzfX8vLCrfzfyUUM7ZHudRwJY4HcnngUuBmoAXDOLQAmBjOUBM/IggzioqOYtXaH11FEJAyU7NrPP5dWM7awC9/7Sl+v40iYC6Q0JDnnZh90TBsYRKiE2GhGFmQwe62uNIh0dnX1jh8/Mx/n4P6LRhIdpVUfpXmBlIbtZtYX/0wJM7sA2BLUVBJUR/bJZNHmvZRr8yqRTu2x99cwa+1OvjkojoLMJK/jSAQIpDRcC/wFGGhmm4AfAVcHNZUE1RG9M6mrd8zVLAqRTmvJ5r3c99oKThuSy3E9WhwTLwIEUBqcc2uccycB3YCBzrnjnHPrgp5MgmZMry7ERBmz1mhcg0hnVFlTx4+e/oyMpFju+towbUYlAWuxXppZPHA+UAjENPyPyzn3q6Amk6BJiothaI90jWsQ6aTufmUZK0rLefzbR5CZHOd1HIkggdyemIZve+xaoKLRQyLYkX0ymV+ymwPVdV5HEZEQemdFGX//cB2XH1PIV4u6eR1HIkwgN7LynXOnBT2JhNSRvTP5yztr+GzjLo7p29XrOCISAjsrqrnx2fn0z07hZ6cP9DqORKBArjR8aGbDgp5EQqq4MJMog1lrdItCpDNwznHzcwvYvb+aByaO1KqP0iqBXGk4DrjczNYCVYABzjk3PKjJJKjSEmIZnJemRZ5EOomnP9nIa4tLufn0gQzJ06qP0jqBlIbTg55CPHFEYRb/mrWeqto64mP0rw6Rjmp1WTm3v7iEY/pm8d3j+3gdRyJYIFMu1x/qEYpwElxH9smkqraeBSV7vI4iIkFSXVvPj56aR3xsFPdfNJIorfoobRDImAbpoMYWZgJo6qVIB3b/GytYuGkPd39tOLnpCV7HkQin0tCJZSbHMSAnlY+1yJNIh/Th6u385d3VXHxEAacNzfU6jnQAnpUGM7vezJaZ2WIzu9d/LMvM3jazcjP7w0Gvn2xmG82s/KDjl5tZmZnN8z++0+i5y8xspf9xWWjOLLIc2SeTuet3UV1b73UUEWlHu/dX8+Nn5tM7K5lJZw32Oo50EJ6UBjMbj2/BqBHOuSHAff6nKoFJwI2HeNuLwBFNfOTTzrmR/sdf/b8jE7gVONL/vlvNrEs7nkaHcEzfruyvrmPext1eRxGRduKc42f/Xcj28ioenDiKpDjtLSHtw6srDVcDdzvnqgCcc9v8Xyucc+/jKw9f4Jz72Dl3OLtrngq84Zzb6ZzbBbwBaJGqgxzdN4sog/dXbfc6ioi0k3/P3siri7fyk1MHMCxf0yul/ZhzLvS/1GwevuWpT8NXEG50zn3S6PnLgWLn3HWHeG+5cy7loNf+GigDVgD/55zbaGY3AgnOuTv9r5sEHHDO3XeIz7wKuAogJydnzFNPPdVepxpy5eXlpKSktPzCRn710QGiDH5xVGKQUrVda84rEui8IksknNfm8npu+/AA/btE8ePiBKIC2IwqEs6rNXRegRs/fvxc51xxS68L2jUrM5sBHGrkzS3+35sJHAWMBZ4xsz6udQ3mReDfzrkqM/se8DhwwuF8gHPuEeARgOLiYjdu3LhWxAgPM2fO5HDzz6lazsPvrGb0UceSlhAbnGBt1JrzigQ6r8gS7udVWVPHeX/6kJTEOv72vePJTgtstkS4n1dr6bzaX9BuTzjnTnLODT3EYxpQAjznfGYD9UCrNkBwzu1ouM0B/BUY4/9+E1DQ6KX5/mNykOP6d6Wu3vHxas2iEIlk9766nKVb9nLfhcMDLgwih8OrMQ1TgfEAZlYExAGtuqluZt0b/XgOsNT//WvAKWbWxT8A8hT/MTnIqJ4ZJMZG84HGNYhErLeXbWPKB2u5/JhCThiY43Uc6aC8GlI7BZhiZouAauCyhlsTZrYOSAPizOxc4BTn3BL/tMxvAElmVgL81Tl3G/ADMzsH39bdO4HLAZxzO83sDqBhrMSvnHNaxegQ4mOiObJPJu+pNIhEpNK9lfz42fkM6p6m3SslqDwpDc65auBbTTxX2MTxm4CbDnH8ZuDmJt4zBV9BkRYc168rd05fyubdB8jLCN8BkSLyRXX1jv97eh4Hquv4/cWjtHulBJVWhBTAN64BNPVSJNL8+Z3VfLh6B7efM4R+2R1vpoCEF5UGAWBATipdU+I1rkEkgsxdv4v731jB2SPyuLA43+s40gmoNAgAZsZx/bL4YNV2vFi7Q0QOz54DNfzwqc/Iy0hg8nlDsQDWYxBpK5UG+dyx/bqyvbyaZVv3eR1FRJrhWyZ6AVv3VPLQxFFhu76KdDwqDfK5hnENukUhEt7++fF6Xlm0lZtOG8ContpSR0JHpUE+1z09kb7dknlvpUqDSLhavHkPd0xfyrgB3fjOcX28jiOdjEqDfMHx/bsxa+0OKmvqvI4iIgepqKrl+ic/o0tSLL+9cARRURrHIKGl0iBf8NWiblTW1DNrrdbBEgk3k6YuYt2OCh74+iiyUuK9jiOdkEqDfMHRfbNIiI3i7WXbvI4iIo08O2cjz322ietP6M/RfbO8jiOdlEqDfEFCbDTH9u3KW8u2aeqlSJhYUbqPSdMWcVSfTH5wYn+v40gnptIgXzJ+YDYbdu5ndVmF11FEOr391bVc+69PSYmP4aGJo4jWOAbxkEqDfMn4gdkAvLWs1OMkIvLLaYtZVVbOA18fpe2uxXMqDfIlPTISGZibylsa1yDiqf/MLeE/c0u4fny/z9dREfGSSoMc0viB2cxZt4u9lTVeRxHplFaW7mPS1EUc2TuTH55U5HUcEUClQZpwwsBsausd763QQk8ioba/upar//UpyfHRPHSxxjFI+FBpkEMaVZBBRlKsblGIhJhzjlueX8TqsnIenDiKHI1jkDCi0iCHFBMdxVeLujFz+Tbq6zX1UiRUnvpkI89/tokfnVjEsf00jkHCi0qDNOmEgdnsqKhmfslur6OIdAqLNu3h1hcWc3z/rlx/Qj+v44h8iUqDNOmrRd2IMrQ6pEgI7K2s4donPyUzKY4Hvj5S+0pIWFJpkCZlJMUxumcX3lqu0iASTM45bnxmPiW7DvD7b2hfCQlfKg3SrBMGZbNo01627DngdRSRDuuRd9fw+pJSbj59IGMLM72OI9IklQZp1imDcwF4fbFWhxQJho9W7+CeV5dxxrBcrjyut9dxRJql0iDN6pedQv/sFF5ZtMXrKCIdTuneSq7/92cUdk3mnvOHY6ZxDBLeVBqkRacNzWX22p3sKK/yOopIh1FTV891T35KRVUtf/7WGFITYr2OJNIilQZp0alDcql3MGOpblGItJe7X1nGJ+t2cff5wyjKSfU6jkhAVBqkRUPy0ijITOSVRVu9jiLSIbw4fzOPvb+Wy48pZMLIHl7HEQmYSoO0yMw4bUguH6zarg2sRNpoRek+fvrfBYzp1YWfnzHI6zgih0WlQQJy2tDu1NQ5LfQk0gb7Kmv4/j/mkhQXw5++OZq4GP0nWCKL/hcrARlVkEF2ajyvLNQtCpHWcM5x47PzWb9zP3/4hjaiksik0iABiYoyTh2Sy8wV2zhQXed1HJGI8+d31vDaYt8CTkf1yfI6jkirqDRIwE4fmktlTT3vrNAtCpHD8e6KMn7z2jLOHN5dCzhJRFNpkIAd0TuTjKRYXtUsCpGAbdixn+v//Rn9s1P5zQVawEkim0qDBCwmOoqTB+Xw5tJtVNXqFoVISw5U1/G9f87FOcdfLhlDUlyM15FE2kSlQQ7LGcO7s6+qlndXbPc6ikhYc87xs+cWsGzrXh68eBSFXZO9jiTSZioNcliO69eVzOQ4ps7b5HUUkbD22PtrmTZvMzeeMoDxA7K9jiPSLlQa5LDERkdx5rDuzFhSSnlVrddxRMLSeyvLuOvlpZw2JJdrxvX1Oo5Iu1FpkMN27qg8qmrreU0DIkW+ZP2OCq570jfw8bcXjdDAR+lQVBrksI3u2YX8LolMm7/Z6ygiYaW8qpbvPjEHM3j00mKS4zXwUToWlQY5bGbGhJF5vL+yjLJ92i5bBKC+3nHD0/NYXVbBH78xmp5ZSV5HEml3Kg3SKueO7EG9g5cW6GqDCMCDb67k9SWl/PyMQRzbr6vXcUSCQqVBWqV/TiqDuqcxbZ5Kg8j0BVt48M2VXDAmn28fW+h1HJGgUWmQVjt3ZB7zNu5m3fYKr6OIeGbRpj38+Nl5jOnVhcnnDdXAR+nQVBqk1c4ekYcZvKABkdJJbdtXyXefmENmUhx//tYY4mOivY4kElQqDdJqeRmJHFGYydR5m3DOeR1HJKQqa+r43j/msnt/DY9eVky31HivI4kEnUqDtMmEkT1YU1bBgpI9XkcRCRnnHD9/biGfbdjN774+giF56V5HEgkJlQZpkzOHdychNoqn52z0OopIyPxp5mqe+2wTN5xcxGlDu3sdRyRkVBqkTdITYzljWHdenLeZ/dVaVlo6vlcXbeE3ry1nwsg8rj+hn9dxREJKpUHa7OvFBeyrquXlhVpWWjq2dXvq+L+n5zOqZwb3nD9cMyWk01FpkDY7oncmvbsm88wnukUhHVfp3koe+LSKLkmxPHJJMQmxmikhnY9Kg7SZmXFRcQGz1+1kdVm513FE2t3+6lqufPwTDtQ6Hrt8rGZKSKel0iDt4vwxPYiOMl1tkA6nrt7xg3/PY8nmvVw9Ip5B3dO8jiTiGZUGaRfZqQmcODCb/35aQk1dvddxRNrN5OlLmbG0lFvPHsLIbO1aKZ2bSoO0m6+PLWB7eTVvLt3mdRSRdvH4h+uY8sFarji2kMuOKfQ6jojnPCsNZna9mS0zs8Vmdq//WJaZvW1m5Wb2h0avTTKz6Y1ef3ej5+LN7GkzW2Vms8yssNFzN/uPLzezU0N5fp3RV4u6kZMWz9OfbPA6ikibvbWslNtfXMxJg7L5xZmDvY4jEhY8KQ1mNh6YAIxwzg0B7vM/VQlMAm48xNvuc84NBEYBx5rZ6f7jVwK7nHP9gN8B9/h/x2BgIjAEOA34k5lpuHMQxURHccGYfN5ZUcaWPQe8jiPSagtL9nDdk58xOC+NByeOIjpKUytFwLsrDVcDdzvnqgCcc9v8Xyucc+/jKw+fc87td8697f++GvgUyPc/PQF43P/9f4ATzTd5egLwlHOuyjm3FlgFHBHc05KJY3vigH99rKsNEpk27tzPFX//hC5JcUy5bCzJ8RrHINLAvNhoyMzmAdPwXQGoBG50zn3S6PnLgWLn3HWHeG8GvtJwknNujZktAk5zzpX4n18NHAncBnzsnPun//hjwCvOuf8c4jOvAq4CyMnJGfPUU0+149mGVnl5OSkpKZ5mePDTSlbtruO3X00iLrp9/oUWDucVDDqv8FJe7Zg86wB7qhy/OCqRvJQv/rsqUs+rJTqvyBKM8xo/fvxc51xxS68LWoU2sxlA7iGeusX/ezOBo4CxwDNm1se10GDMLAb4N/CQc25Ne2V1zj0CPAJQXFzsxo0b114fHXIzZ87E6/yx+dv55l9nsTejPxeMyW/5DQEIh/MKBp1X+KiqreOSx2azo7KSf1x5JEf2yfrSayLxvAKh84osXp5X0EqDc+6kpp4zs6uB5/wlYbaZ1QNdgbIWPvYRYKVz7oFGxzYBBUCJv1SkAzsaHW+Q7z8mQXZM3yz6Z6fw9w/Xcv7oHlpqV8Jefb3jx8/MZ/banTx08ahDFgYR8W5Mw1RgPICZFQFxwPbm3mBmd+IrBD866KkXgMv8318AvOUvIy8AE/2zK3oD/YHZ7XYG0iQz47JjClm0aS+fbtjldRyRZjnnuHP6Ul5asIWbTx/IOSPyvI4kEra8Kg1TgD7+8QhPAZc13Jows3XA/cDlZlZiZoPNLB/fbY3BwKdmNs/MvuP/rMeALDNbBdwA/AzAObcYeAZYArwKXOucqwvZGXZyXxvdg9SEGP72wTqvo4g069H31ny+FsNVX+njdRyRsObJsGD/DIhvNfFcYRNvO+Q1budcJXBhE89NBia3IqK0UVJcDBPHFjDlg3Vs2XOA7umJXkcS+ZLnPyvhrpeXcebw7kw6c7BupYm0QCtCStBcenQh9c5p+qWEpXdXlPGTZxdwdJ8s7r9oBFFai0GkRSoNEjQFmUmcODCHf8/eQGWN7gxJ+Ji/cTff/+dc+mWn8JdLxxAfo3XfRAKh0iBBdcWxheyoqGbqZ5q4IuFh1bZyLv/bbLJS4nji20eQlhDrdSSRiKHSIEF1TN8shvZI4y/vrqGuPvQLiYk0tmXPAS59bBbRUVH849tHkp2W4HUkkYii0iBBZWZcM64fa7dX8MqiLV7HkU5sV0U1lzw2m32Vtfz9irEUdk32OpJIxFFpkKA7dUgufbom86e3V+PFsuUi5VW1XPH3T9iwcz+PXlbM0B7pXkcSiUgqDRJ00VHG97/alyVb9vLOipYW/RRpX5U1dVz1xBwWbtrDHy4exVFa7VGk1VQaJCTOHdWD7ukJ/Gnmaq+jSCdSU1fPdU9+yoerd3DfhcM5ZcihtsMRkUCpNEhIxMVE8d3j+zB77U7mrNvpdRzpBOrrHTc+O58ZS7dxx4QhnDeqfTZPE+nMVBokZCYeUUCXpFhdbZCgc87xi2mLmDZvMzedNoBLji70OpJIh6DSICGTFBfDFcf25q1l21iyea/XcaSDcs5xx0tLeXLWBq4e15drxvXzOpJIh6HSICF12dGFpMbH8LsZK7yOIh2Qc457X1v++QZUN506wOtIIh2KSoOEVHpSLFd9pQ9vLCnlM22bLe3soTdX8fDM1XzjyJ788ixtQCXS3lQaJOSuOK43Wclx3Pf6cq+jSAfy53dW87sZK7hgTD53ThiqwiASBCoNEnIp8TFcM74fH6zawQertnsdRzqAv763hrtfWcbZI/K45/zh2rFSJEhUGsQT3zyyJ3npCfzmteVaJVLa5K/vreHO6Us5c1h3fnfRCKJVGESCRqVBPJEQG80PT+rPvI27mbF0m9dxJEI1FIYzhuXywMSRxETrP2kiwaS/YeKZ80fn07trMve9tpx67YAph2nK+2u5c/pSTh+ay4MTRxGrwiASdPpbJp6JiY7ihpOLWF66j6nzNnkdRyLIo++u4VcvLeG0Ibk8dLEKg0io6G+aeOrMYd0ZkZ/OPa8uo6Kq1us4EgH++PYqJr/sG8Pw+2+oMIiEkv62iaeiooxbzxlC6d4q/vj2Kq/jSBhzzvHAjBX85rXlnDsyjwcnjlRhEAkx/Y0Tz43u2YWvjerBX99by/odFV7HkTDknOO+15fzwIyVXDAmn99epEGPIl7Q3zoJCz89fSAx0cad05d6HUXCTH294/YXl/DHt1dz8RE9uff84ZpWKeIRlQYJCzlpCVx3Qj/eWFLKeyvLvI4jYaK2rp4b/zOfv3+4ju8c15u7zhuqhZtEPKTSIGHjyuN60ysridtfXEJNXb3XccRjVbV1XPOvT3nu003ccHIRt5w5SEtDi3hMpUHCRnxMNL84czCrtpXz+IfrvI4jHqqoquU7j8/h9SWl/PKswfzgxP4qDCJhQKVBwspJg7I5cWA2v319BRt27Pc6jnhgR3kV33j0Yz5YtZ17LxjOt4/r7XUkEfFTaZCwYmbced5QoqOMnz23QPtSdDIlu/Zz4Z8/YtnWffzlkmIuKi7wOpKINKLSIGGne3oiN58xkA9X7+CZORu9jiMhsmzrXs5/+EO2l1fxz+8cycmDc7yOJK509rYAABMoSURBVCIHUWmQsHTx2J4c2TuTO6cvpXRvpddxJMg+XrODi/78EQDPfv8YxhZmepxIRA5FpUHCUlSUcff5w6mureeW5xfpNkUHNm3eJi55bBbZaQn89+pjGJCb6nUkEWmCSoOErd5dk7nh5CJmLC1l1tY6r+NIO3PO8ce3V/HDp+YxumcX/vv9Y8jvkuR1LBFphkqDhLUrj+vNyIIMHl9cxcadmk3RUdTW1fPz5xfxm9eWM2FkHk9ceQTpSbFexxKRFqg0SFiLiY7ioYmjAPjhU59Rq0WfIt6eAzVc8fdP+PfsDVw3vh8PfH0k8THRXscSkQCoNEjY65mVxGWD4/l0w24e/P/27jzOqrr+4/jrMzvbsG/DjoAJssmiKKhoKmZiEoqKBampidmP0h6Z9tDKitT6pfLItNA0fwVprrmgIbiQCLLIDsMOA7ILDMswM3x+f9xDjysNzGHuvXO4w/v5eNwHZ537+XDPnPnc7/d7zplSGHU4koA12/Zy5e+nM2PVdh4a3oM7LzlVN20SSSNZUQcgEsZZBVlsz27K+KkrOPuUJgw4pXHUIclx+mjldm59bjYZBv9301n076ArJETSjVoaJG389IpudGhSh/+ZNJcdew9GHY6E5O48N2Mt35jwMU3r5fLKmIEqGETSlIoGSRu1c7J47Nre7NxbqvENaaKkrJy7X1zAvS8vZGDnJrx429m0bawrJETSlYoGSSvdCurz869144PCbYx7c2nU4cgxbN59gGuenMHEWeu5fXAnJozqR36erpAQSWca0yBpZ0S/tizZtIc/fbiaU1vU4yo9n+CEs3xnOXc99iF7S8p4fOQZXNq9ZdQhiUgSqGiQtHTvZaexYksx97y0kI5N69KnXcOoQxJi4xeeeH8VD848QNtGtXnuxjN1h0eRGkTdE5KWsjIzGH9db1o2yOOWv8xm0679UYd00tu1r5RvP/sJ495cyhnNMnntuwNVMIjUMCoaJG01qJ3Dn77ZlwOl5Xzr6Vns2lcadUgnrXnrP+eyxz7gveVbue/yrozplUs9jV8QqXFUNEha69y8Hk98ow+rtu7lhmdmse9gWdQhnVTKD8WeHzH88X/jDn+/ZQDfOqeDbtgkUkOpaJC0d06nJjx6bS/mrtvJrc/N4WCZLsWsDhs/3891f5zBQ5OXMeT0FrzxvUH0bquxJSI1mYoGqRGGnN6SccN68P7yrYz9+zzKD+lR2qn0z/kbufSRD1hQtIuHr+rJY9f2pn4tdUeI1HS6ekJqjKv7tWHX/lJ+8cYS6uRk8qthPcjMUDN5Mu3Ye5CfvLKQ1+dvomebBjwyohftm9SJOiwRqSYqGqRG+fa5HSkuKeORKYXsPVjO/17di5wsNaglw9uLPuPHLy1g1/5S7rrkVG45tyNZmfq/FTmZqGiQGmfsRV2ok5vJL99YGtxcqA+1cvTo5araVlzCz/+5mFfmbaRry3z+cuOZnNYyP+qwRCQCKhqkRrr53FOol5fNj19awKinZzJhVF9dAnic3J1/zCnigdcXs7ekjO9d2Jkxgzup5UbkJKaiQWqsa/u3pW5uFmMnzeOqP3zEH7/ZlzaN9LCkMNZs28s9Ly9g+ort9G3XkF8N607n5rpRk8jJTkWD1GiX9yygfq1sxvx1DkPHf8jvR/ZhwCmNow7rhLW3pIzxU1cw4YPV5GZl8MDXTue6/m3J0IBSESHCSy7N7LtmttTMFpnZg8GyxmY21cyKzWx83La1zez1uO3Hxa0bbWZbzWxe8Lopbt0oMysMXqOqN0M5UZzbpSmv3j6QxnVzuX7Cx/x5+mrcdUlmPHfnlXlFXPCbaTw+bSVf7dmSKT84j+vPaqeCQUT+I5KWBjMbDFwB9HT3EjNrFqw6APwEOD14xXvY3aeaWQ4wxcwudfc3g3WT3P32I96jEXAf0BdwYLaZveruO1OUlpzAOjSpw0u3nc3YSfO4/7XFzC/axU+HdtM4B+DjVdsZ99ZS5q77nO6t6vP7kX30ADARqVBU3RPfAca5ewmAu28J/t0LfGhmneI3dvd9wNRg+qCZzQFaV/IelwDvuPsOADN7BxgC/C2ZiUj6qJeXzZPf6MsjUwp57N1CZq7ewW+v7kX/Do2iDi0SSzbt5sG3ljJ12Vaa5+fy669356o+bdSyICJHFVX3RBdgkJl9bGbvmVm/sDuaWQPgcmBK3OKvm9l8M3vBzNoEy1oB6+O22RAsk5NYRoYx9qIuPH/rADLMGPHkR/zqjSWUlJVHHVq1Kdy8h+9NnMtXHv2A2Wt38qNLv8R7dw1mRD+NXRCRY7NU9e2a2b+AFhWsugf4BbGWgzuAfsAkoKMHwZjZaKBvBV0OWcBrwGR3/12wrDFQHHRz3AKMcPcLzOxOIM/dHwi2+wmw390friDWm4GbAZo3b95n4sSJCecfleLiYurWrRt1GEmXirwOlDkTlx1k2voyCuoaI7+US7cm1Xs/h+r8vNbuLue1laXM3lxOTiZc2DabyzpmUyc7+YWCjsP0orzSSyryGjx48Gx371vZdikrGo75pmZvAb9296nB/ErgLHffGsyPpuKi4SliBcIdR/m5mcAOd69vZtcC57v7LcG6J4Bp7n7M7om+ffv6J598kliCEZo2bRrnn39+1GEkXSrzmrpsC/e9soh1O/ZxSbfm3HtZ12q7NDPVn5e7M33FdiZ8uIqpy7ZSLzeLUWe354aBHWhUJydl76vjML0or/SSirzMLFTRENWYhpeBwcBUM+sC5ADbjrWDmT0A1AduOmJ5S3ffFMwOBZYE05OBX5rZ4RFdFwN3Jyd8qUkGn9qMAWMbM+HD1Yx/dwUXLnuPGwd24KaBHWhcNzfq8KrkQGk5L80t4unpq1m+uZgmdXMY++UujD6nvR4sJSJVFlXR8BTwlJktBA4Co+K6JtYA+UCOmX2N2B/73cS6NZYCc8wMYLy7/wm4w8yGAmXADmA0gLvvMLOfA7OC9/zZ4UGRIkfKy85kzOBOfP2M1ox7cwmPT1vJ09NXc02/tnz73I60alAr6hAr5e58umEXz3+yntc+3cjuA2V0bZnPQ8N7cHnPAvKydSttEUlMJEWDux8Erj/KuvZH2a3Cjld3v5ujtCC4+1PEChSRUFrUz+N31/Tm9gs68fi0VTw3Yy3PzVjL0J4FjOjXhn7tG51wgwVXbS3mrUWf8eKcIlZsKSYvO4Mh3VpwTf+2nNmhEUGRLSKSMN0RUqQCnZrV4zdX9+T7F3fhj++v4oXZG3hxbhFtGtViWO/WDDujFe0aR/NI6PJDzsKiXbyzeDOTF31G4ZZiAPq0a8i4Yd35So+W5Ov+EyKSAioaRI6hVYNa3D+0Gz8cciqTF33GP2YX8ei7hTwypZCOTetwbuemnNelKWd2bETtnNT8OpWWH6JwczEfr97Ov1duZ8aq7ew5UEZmhtG/fSNGntmWi7u1oCANulBEJL2paBAJoXZOFlf2bs2VvVuz8fP9vLnwM95fvpWJs9bx53+vISvD6NSsLl0L8unaMp/TWuZT0KAWzfNzQxcTew6UUvT5fop27mf9jn0s/WwPizbuZtnmPRwsOwRA20a1uax7Swac0phBnZum9AoIEZEjqWgQOU4FDWpx48AO3DiwAwdKy5m1ZgczVm1n0cbdfFC4jRfnFH1h+3q5WTTNzyUvK5OcrAxyszLIzsxg38EyNu/YBzPeZff+UvaUlH1hvwa1s+lWkM/os9vTrSCfPu0a0rqhntIpItFR0SCSgLzsTAZ1bsqgzk3/s2zLngMUbi5m8+4DbN5dwubdB9haXEJJaTklZYcoKTvE3oNl1MnJokWdDNq3aky9vCxa1M+jVYNatGpYi9YNatG0Xq4GMYrICUVFg0iSNauXR7N6eaG2jd2kpWeKIxIRSY7IHo0tIiIi6UVFg4iIiISiokFERERCUdEgIiIioahoEBERkVBUNIiIiEgoKhpEREQkFBUNIiIiEoqKBhEREQlFRYOIiIiEoqJBREREQlHRICIiIqGoaBAREZFQVDSIiIhIKCoaREREJBQVDSIiIhKKigYREREJRUWDiIiIhGLuHnUMJxQz2wqsjTqOBDQBtkUdRAoor/SivNKL8kovqcirnbs3rWwjFQ01jJl94u59o44j2ZRXelFe6UV5pZco81L3hIiIiISiokFERERCUdFQ8zwZdQAporzSi/JKL8orvUSWl8Y0iIiISChqaRAREZFQVDSIiIhIKCoa0oSZDTGzZWa2wsx+VMH6XDObFKz/2Mzax63rYWYfmdkiM1tgZnnVGfuxVDUvM8s2s2eCfJaY2d3VHfuxhMjrXDObY2ZlZjb8iHWjzKwweI2qvqgrV9W8zKxX3DE438xGVG/klUvkMwvW55vZBjMbXz0Rh5PgsdjWzN4OfscWx59XopZgXg8Gx+ISM3vUzKz6Ij+2EHl9P/gs5pvZFDNrF7cu9ecOd9frBH8BmcBKoCOQA3wKdD1im9uAPwTT1wCTguksYD7QM5hvDGRGnVMS8roOmBhM1wbWAO2jzuk48moP9ACeBYbHLW8ErAr+bRhMN4w6pyTk1QXoHEwXAJuABlHnlIzc4tY/AvwVGB91PsnKC5gGXBRM1wVqR51TEo7Fs4Hpwc/IBD4Czo86p+PIa/DhzwH4Ttw5sVrOHWppSA/9gRXuvsrdDwITgSuO2OYK4Jlg+gXgwqB6vhiY7+6fArj7dncvr6a4K5NIXg7UMbMsoBZwENhdPWFXqtK83H2Nu88HDh2x7yXAO+6+w913Au8AQ6oj6BCqnJe7L3f3wmB6I7AFqPTuc9Uokc8MM+sDNAfero5gj0OV8zKzrkCWu78TbFfs7vuqKe7KJPJ5OZBH7I9yLpANbE59yKGEyWtq3OcwA2gdTFfLuUNFQ3poBayPm98QLKtwG3cvA3YRa1XoAriZTQ6a6n5YDfGGlUheLwB7iX1jXQc87O47Uh1wSGHySsW+qZaU2MysP7ET9sokxZUMVc7NzDKA3wB3piCuRCXymXUBPjezF81srpk9ZGaZSY+waqqcl7t/BEwldu7YBEx29yVJj7BqjjevG4E3q7hvlahoqPmygIHAyODfK83swmhDSor+QDmxpu4OwA/MrGO0IUllzKwl8BfgW+7+X9/Y09RtwBvuviHqQJIsCxhErBjqR6zJfHSUASWDmXUCTiP2Db0VcIGZDYo2quNnZtcDfYGHqvN9VTSkhyKgTdx862BZhdsETfb1ge3Eqs333X1b0KT1BnBGyiMOJ5G8rgPecvdSd99CrI/yRLnHfJi8UrFvqiUUm5nlA68D97j7jCTHlqhEchsA3G5ma4CHgW+a2bjkhldlieS1AZgXNJWXAS+TXueOo7kSmBF0txQT+6Y+IMnxVVWovMzsy8A9wFB3LzmefROloiE9zAI6m1kHM8shNiDw1SO2eRU4PFp2OPCux0bHTAa6m1nt4I/uecDiaoq7MonktQ64AMDM6gBnAUurJerKhcnraCYDF5tZQzNrSGxMyuQUxXm8qpxXsP1LwLPu/kIKY6yqKufm7iPdva27tyf2rfxZd/+vUe8RSeRYnAU0MLPDY08uIL3OHUezDjjPzLLMLJvYOfFE6Z6oNC8z6w08Qaxg2BK3qnrOHVGPFtUr9KjarwDLifUD3xMs+1lw4EBsYM/zwApgJtAxbt/rgUXAQuDBqHNJRl7ERnI/H+S1GLgr6lyOM69+xL7J7SXWcrIobt8bgnxXEGvGjzyfRPMKjsFSYF7cq1fU+STrM4v7GaM5ga6eSMKxeBGxq68WAH8GcqLOJwnHYiaxP7pLgnPHb6PO5Tjz+hexgZuHf49ejds35ecO3UZaREREQlH3hIiIiISiokFERERCUdEgIiIioahoEBERkVBUNIiIiEgoKhpEREQkFBUNIiIiEoqKBhE5IZlZcdQxiMgXqWgQERGRUFQ0iEjCzKynmb1vZovN7JCZuZn9LG79ODMbEzd/v5ndGUy/bGazzWyRmd1cwc9ub2YL4+bvNLP7g+nrzWymmc0zsydOoEc3i9RIKhpEJCFmlgdMAu50967AL4g97fG+uM0mAVfHzV8dLAO4wd37EHtK6R1m1jjk+54GjADOcfdexB6VPjKRXETk2LKiDkBE0t6XgTnuPjOYnw8M8bgH27j7XDNrZmYFQFNgp7uvD1bfYWZXBtNtgM7EHjBUmQuBPsAsMwOoBWw55h4ikhAVDSKSqNOJPQXxsDOAORVs9zyxx5u3IGhlMLPziRUdA9x9n5lNI/Zk03hlfLFV9PB6A55x97sTjF9EQlL3hIgkajvQA8DMugDDgIkVbDcJuIZY4fB8sKw+sVaHfWb2JeCsCvbbDDQzs8Zmlgt8NVg+BRhuZs2C925kZu2SlJOIVEAtDSKSqL8BQ4PBituAa939v7oX3H2RmdUDitx9U7D4LeBWM1sCLANmVLBfaTCociZQBCwNli82s3uBt80sAygFxgBrk56hiABgcd2OIiIiIkel7gkREREJRUWDiIiIhKKiQUREREJR0SAiIiKhqGgQERGRUFQ0iIiISCgqGkRERCSU/wePLxQpu8Ez6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112f80a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.plot(sigma_grid, nll_values)\n",
    "plt.xlabel('$\\sigma$ value')\n",
    "plt.ylabel('negative log likelihood')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### getting the location of the minimum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "minimum NLL (-61257.703125) is at 0.099500\n"
     ]
    }
   ],
   "source": [
    "index = np.argmin(nll_values)\n",
    "print(\"minimum NLL (%f) is at %f\" % (nll_values[index], sigma_grid[index]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the minimum value of sigma using an optimizer "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate the derivative of the negative log likelihood with respect to sigma:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "nll_derivatives = [\n",
    "    session.run(tf.gradients(nll, sigma),\n",
    "           feed_dict = { x : xpoints, \n",
    "                        sigma: sigma_val, \n",
    "                        bin_contents_ph : bin_contents\n",
    "                       })[0]\n",
    "    \n",
    "    for sigma_val in sigma_grid\n",
    "]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For illustration purposes we draw arrows corresponding to the derivatives at selected points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "arrow_points = np.arange(7) * 0.02 + 0.06\n",
    "\n",
    "arrow_x_size = 0.05\n",
    "\n",
    "arrow_points_y = []\n",
    "arrow_points_derivative = []\n",
    "\n",
    "for sigma_val in arrow_points:\n",
    "    # (re)calculate the function value and derivative\n",
    "    # at some points\n",
    "    nll_val, nll_derivative = session.run([\n",
    "        nll, \n",
    "        tf.gradients(nll, sigma)[0]\n",
    "        ],\n",
    "               \n",
    "        feed_dict = { x : xpoints, \n",
    "                      sigma: sigma_val, \n",
    "                      bin_contents_ph : bin_contents\n",
    "                      })\n",
    "\n",
    "    \n",
    "    arrow_points_y.append(nll_val)\n",
    "    arrow_points_derivative.append(nll_derivative)\n",
    "    \n",
    "arrow_points_y = np.array(arrow_points_y)\n",
    "arrow_points_derivative = np.array(arrow_points_derivative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAHjCAYAAACgviz7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0VVXexvHvTu8JJBACiYQiSkdqCOMrjF1nxI4d7DKOOiJW7AyOYlfGNiN2BWwoFixILPTeCR1CTYGQXu9+/8glEzRACEnOvcnzWesucur9bUF5PLscY61FRERE5Eh8nC5AREREvINCg4iIiNSIQoOIiIjUiEKDiIiI1IhCg4iIiNSIQoOIiIjUiEKDiIiI1IhCg4iIiNSIQoOIiIjUiJ/TBXiamJgYm5iY6HQZtZafn09oaKjTZdQ5tcu7qF3eRe3yLvXRrkWLFmVaa1sc6TyFht9JTExk4cKFTpdRaykpKQwePNjpMuqc2uVd1C7vonZ5l/polzFma03OU/eEiIiI1IhCg4iIiNSIQoOIiIjUiEKDiIiI1IhCg4iIiNSIQoOIiIjUiEKDiIiI1IhCg4iIiNSIQoOIiIjUiEKDiIiI1IhCg4iIiNSIQoOIiIjUiEKDiIiI1IhCQz3KKSrli6U7SNtb4HQpIiIix0yhoR7tLyjljklLmb0x0+lSREREjplCQz2KiwzC18ewfV+h06WIiIgcM4WGeuTn60OriCB1T4iISKOg0FDPEpoHk6YnDSIi0ggoNNSz+GYhbN+nJw0iIuL9FBrqWUKzEPbkFFNUWu50KSIiIsdEoaGeJTQPBmBntrooRETEuyk01LP4ZiEAGtcgIiJeT6Ghnh140qAZFCIi4u0UGupZy/Ag/H21VoOIiHg/x0KDMeY2Y8xaY8wqY8z4KvvvN8ZsMMakGmPOrLL/DmPMSvf5/6iy/1FjzA5jzFL355wj3ash+foY2kQFk6YZFCIi4uX8nPhSY8wQYCjQ01pbbIxp6d7fBbgM6Aq0Bn40xnQCOgM3Av2BEmC6MeYra+0G9y2ft9Y+87vvqPZe1toGn8aQ0DyE7eqeEBERL+fUk4aRwJPW2mIAa226e/9QYJK1tthauxnYQEVQ6AzMs9YWWGvLgJ+BC4/wHYe6V4OLbxas7gkREfF6jjxpADoBJxtjxgFFwGhr7QKgDTC3ynnb3ftWAuOMMdFAIXAOsLDKeX83xlzj3neXtXbfYe71B8aYm4CbAGJjY0lJSTnmBlZVuq+ErPxSpv84kyA/U6f3/r28vLw6r98TqF3eRe3yLmqXd3GyXfUWGowxPwKtqjk0xv29zYEkoB8wxRjT/lD3stauMcY8BXwP5ANLgQPdDK8CYwHr/vVZ4LqjqdVa+wbwBkDfvn3t4MGDj+byI8pptpNP1i+hffe+dIoNr9N7/15KSgp1Xb8nULu8i9rlXdQu7+Jku+otNFhrTzvUMWPMSOAza60F5htjXEAMsANIqHJqvHsf1to3gTfd1z9BxZMDrLV7qtz3P8BX7s1D3quhxTf737TL+g4NIiIi9cWpMQ1TgSEA7oGOAUAm8CVwmTEm0BjTDjgemO8+78BgyeOoGM/woXs7rsp9L6CiK4PD3auhJRxY4EmDIUVExIs5NaZhIjDRGLOSitkQw91PHVYZY6YAq4Ey4NYqsx0+dY9pKHXvz3bvH2+M6UVF98QW4GYAa+3h7tWgYsICCPL30WBIERHxao6EBmttCXDVIY6NA8ZVs//kQ5x/9WG+p9p7NTRjDPHNQrRWg4iIeDWtCNlAEpoFk7ZXTxpERMR7KTQ0kPhmIWzXkwYREfFiCg0NJKF5MDlFZewvLHW6FBERkVpRaGgg8ZpBISIiXk6hoYEcmHapGRQiIuKtFBoaSELzigWeNK5BRES8lUJDA4kM9ics0E9PGkRExGspNDSQirUagjWmQUREvJZCQwNKaK4FnkRExHspNDSg+GbBbN9XSMWK2SIiIt5FoaEBJTQLoaCknL35JU6XIiIictQUGhpQ5SuyNRhSRES8kEJDA0pofmCtBo1rEBER76PQ0IAOhIatWQoNIiLifRQaGlBYoB8xYQFsU2gQEREvpNDQwNpGh7IlK9/pMkRERI6aQkMDaxsdwjYt8CQiIl5IoaGBJUaHsmt/EUWl5U6XIiIiclQUGhpY2+iKwZB62iAiIt5GoaGBJUaHArAlU+MaRETEuyg0NLADoUHTLkVExNsoNDSwyBB/okL8NYNCRES8jkKDA9pGh+pJg4iIeB2FBgckRofoSYOIiHgdhQYHtI0OZWd2IcVlmnYpIiLeQ6HBAYnRIbgsbNfbLkVExIsoNDigbeUMCnVRiIiI91BocECie4GnLZkaDCkiIt5DocEBzUMDCA/006qQIiLiVRQaHGCMoW2MZlCIiIh3UWhwiNZqEBERb6PQ4JDE6BDS9hZQVu5yuhQREZEaUWhwSNvoUMpclp3ZRU6XIiIiUiMKDQ6pfNulxjWIiIiXUGhwyIFpl1qrQUREvIVCg0NahAcS7O/LFg2GFBERL6HQ4BBjDG2jQ/SkQUREvIZCg4MSo0P1pEFERLyGQoOD2kaHsG1vAS6XdboUERGRI1JocFDb6FBKylzsztG0SxER8XwKDQ6qfHGVxjWIiIgXUGhwUNuYA6/I1rgGERHxfAoNDoqLCCLQz4dNGXlOlyIiIl4gu6CED9YUk1tU6sj3OxoajDG3GWPWGmNWGWPGu/dFG2NmGmPyjDETfnd+H2PMCmPMBmPMS8YY497f3BjzgzFmvfvXZu79xn3eBmPMcmNM74Zv5aH5+BjaxYSyKUPdEyIicnjWWsZ8vpKftpU59oTasdBgjBkCDAV6Wmu7As+4DxUBDwGjq7nsVeBG4Hj35yz3/vuAGdba44EZ7m2As6uce5P7eo/SoUUYmzIVGkRE5PCmLt3B1yt2cUFHf7q1iXSkBiefNIwEnrTWFgNYa9Pdv+Zba3+jIjxUMsbEARHW2rnWWgu8C5zvPjwUeMf98zu/2/+urTAXiHLfx2O0bxHKtr0FlJTpbZciIlK9HdmFPDx1FX3bNuOc9v6O1eHn2DdDJ+BkY8w4KgLCaGvtgsOc3wbYXmV7u3sfQKy1dpf7591AbJVr0qq5ZleVfRhjbqLiSQSxsbGkpKQcdWNqqzizjHKX5ZPpKbQOO/YMl5eX16D1NxS1y7uoXd5F7fJsLmsZv6CI0jIXw9oWUZBf4Fi76jU0GGN+BFpVc2iM+7ubA0lAP2CKMaa9+ylCrVlrrTHmqO5hrX0DeAOgb9++dvDgwcdSwlFpvj2bN5bPonliFwZ3q+4f1dFJSUmhIetvKGqXd1G7vIva5dn+88sm1u5dw/iLe3BJ3wRH21WvocFae9qhjhljRgKfuUPCfGOMC4gBMg5xyQ4gvsp2vHsfwB5jTJy1dpe7+yG9yjUJh7jGI7RvEQbApkzNoBARkYOt3Z3D09+lckaXWC7pE3/kC+qZk2MapgJDAIwxnYAAIPNQJ7u7H3KMMUnuWRPXAF+4D38JDHf/PPx3+69xz6JIAvZX6cbwCGGBfsRGBLIxXYMhRUTkf4pKy/nHpKVEBPvzrwu7454w6CgnxzRMBCYaY1YCJcDwA10TxpgtQAQQYIw5HzjDWrsa+BvwNhAMfOv+ADxJRffG9cBW4FL3/m+Ac4ANQAFwbf036+hVzKDQkwYREfmfZ75LZe3uXN66th/RYYFOlwM4GBqstSXAVYc4lniI/QuBbtXszwJOrWa/BW49pkIbQPsWoXy5dCfWWo9IkiIi4qxZGzL572+buTqpLUNOaOl0OZW0IqQH6NAijJyiMrLyS5wuRUREHJZdUMJdU5bRvkUoD5zT2elyDqLQ4AEODIbcmK4uChGRpsxay5ipK8nMK+bFYScRHODrdEkHUWjwAB1aVLy4SitDiog0bZ8v2cHXy3dx5+md6B7vzKqPh6PQ4AFaRwYT5O+jJw0iIk1Y2t4CHv5iFf0Sm3HLKR2cLqdaCg0eoOLFVXoHhYhIU1VW7uLOyUsxwHOX9sLXxzMHxSs0eIj2LULZqFdki4g0Sa+mbGTh1n08fn5XEpqHOF3OISk0eIgOLcJI21tAcVm506WIiEgDWpqWzQsz1nNez9ac36vNkS9wkEKDh+jQIhSXxbF3pIuISMPLLy7jH5OW0CoiiLHnd/P4tXoUGjxEhwPvoFAXhYhIkzH2q9Vs3VvAs5f2JDLYuVde15RCg4doF1Mx7XJjhgZDiog0BdNX7mbSgjRuOaUDSe2jnS6nRhQaPERooB9xkUEaDCki0gTsySni/s+W061NBHee1snpcmpMocGDVMyg0JMGEZHGzOWyjP54GYWl5bww7CQC/Lznr2LvqbQJ6NAijE0Zebhf9ikiIo3Q27O38Ov6TB48twsdW4Y5Xc5RUWjwIO1jQsktKiMjr9jpUkREpB6k7s7lyelrOa1zS64ccJzT5Rw1hQYP0qHlgRdXqYtCRKSxKSot545JS4gI8uPJi3p4/PTK6ig0eJADb7vclKnBkCIijc1T09eydncuT1/Sk5iwQKfLqRWFBg8SFxFEsL8vG/TiKhGRRiUlNZ23Zm1hRHIiQ05o6XQ5tabQ4EF8fAwdW4axfo9Cg4hIY5GVV8zoj5dzQmw49519otPlHBOFBg/TKTacdXtynS5DRETqgLWWez5ZTk5RKS9e3osgf1+nSzomCg0eplNsGOm5xWQXlDhdioiIHKP3521jxtp07jvrRE5sFeF0OcdMocHDdGoVDsA6dVGIiHi19Xty+edXq/m/Ti0YkZzodDl1QqHBw3SKPRAa1EUhIuKtikrLue2jJYQF+vHMJT3w8fG+6ZXV8XO6ADlY68ggwgL9WK/QICLitcZPT2Xt7lwmjuhLy/Agp8upM3rS4GGMMRwfG0aqQoOIiFdKSU1n4qzNDB/Ylj+fGOt0OXVKocEDdWoZrmmXIiJeKDOvmNEfL6Nds0DuP6ez0+XUOYUGD3R8bBhZ+SVk6h0UIiJew1rLnR8tIiunkKs6lHr99MrqKDR4oBNaaTCkiIi3GTvpF37duI+S+ZMYfv7pTpdTLxQaPNCBGRTqohAR8Xwul4v7n5rAm4v2UrBxAcP6xOHn1zjnGTTOVnm5luGBRAT56UmDiIiHy8zM5Jprr2dpyzPxDc4n65sXGD77Z6fLqjd60uCBjDGc0ErLSYuIeLJffvmFXr16MbeoNQEt2pL5zfN069iWnj17Ol1avVFo8FDHx4azbk8e1lqnSxERkSrKy8sZO3YsQ4YMYW9QGyL6/IWcBVMp2ryYa665xuny6pVCg4fq1DKM/YWlpOdqBoWIiKfYtWsXZ5xxBg8//DAmJIroc+6gZM8m9v38Nj4+PlxxxRVOl1ivFBo8VCfNoBAR8SjfffcdPXv25KeffgIM0eeOwvgFkvHleCgv4/TTTycuLs7pMuuVQoOH+t87KDSDQkTESaWlpdx3332cddZZZGRkABAx4EKCE3uxb8YblO3dDtDouyZAsyc8VkxYINGhAazbrScNIiJOycrKYujQocyaNatyX0BcJ6JOvpr8tb+St/x7AMLCwjj//POdKrPB6EmDBzs+Nox16QoNIiJOiY6OZvr06cydO5fXX3+dyJhWxPz1bsrz9pI1fULleRdffDEhISEOVtowFBo8WKfYindQaAaFiIhzwsLCGDBgAIWFhfgNuAK/yJZkTnsaW5xfeU5T6JoAhQaP1ik2nLziMnbuL3K6FBGRJm3BggU8PPErwrr9mf2zPiKyNIugoIpXXickJHDKKac4XGHDUGjwYP8bDKkuChERp2RnZ3PJdbcS+eebKNq2gpy5H/PBBx8wcuRIAK6++mp8fJrGX6dNo5VeqlNsGADrFRpERBxhrWXEdTdQ3PsKrKuMzK+e4aEHx3Daaadxzz33EBwczNVXX+10mQ1GocGDRYUE0DI8kLWaQSEi4oiXX36Zn/dFEBh3PFnfvsTJfbrz8MMPA9CqVSvefPNNTjzxRIerbDgKDR6uc1wEa3YpNIiINLSFCxcyZsIHRA64iNzFXxG+fxMffvghvr6+ledcfvnlDlbY8BwLDcaY24wxa40xq4wx4937oo0xM40xecaYCb87v48xZoUxZoMx5iVjjHHvf9QYs8MYs9T9OafKNfe7z081xpzZsC2sG11aR7AhPZeSMpfTpYiINBnZ2dlccs31RJ15GyUZW8hOeYv333+/0a/4eCSOLO5kjBkCDAV6WmuLjTEt3YeKgIeAbu5PVa8CNwLzgG+As4Bv3ceet9Y+87vv6AJcBnQFWgM/GmM6WWvL66FJ9aZzXASl5ZYN6Xl0aR3hdDkiIo2etZbrrr+egh6XEBgQTOakMYy57x5OP/10p0tznFNPGkYCT1priwGstenuX/Ottb9RER4qGWPigAhr7VxbsWjBu8CRlt4aCkyy1hZbazcDG4D+ddyOetclriIorN6V43AlIiJNw4QJE5ixwxCceBL7ZvyH5K7teOSRR5wuyyM4tYx0J+BkY8w4KgLCaGvtgsOc3wbYXmV7u3vfAX83xlwDLATustbucx+fe5hrKhljbgJuAoiNjSUlJeXoWlOPXNYS4APfz19FTO6GI56fl5fnUfXXFbXLu6hd3kXt+p/U1FTue+Z1YoY9Qf7aX/HbOo+/P/gffv311/opshac/P2qt9BgjPkRaFXNoTHu720OJAH9gCnGmPa2dksfvgqMBaz712eB647mBtbaN4A3APr27WsHDx5cizLqT+fVs8j19WXw4KQjnpuSkoKn1V8X1C7vonZ5F7WrQnZ2NtfefCvNzh5FeW4me7/7N99+8SlnnHFG/RVZC07+ftVbaLDWnnaoY8aYkcBn7pAw3xjjAmKAjENcsgOIr7Id796HtXZPlfv+B/iqyjUJ1V3jbbrERfDNil1Ya3GP/xQRkTpkreX6G24gr/NfCYlowe4P7uH+u+7wuMDgNKfGNEwFhgAYYzoBAUDmoU621u4CcowxSe5ZE9cAX7ivrzqU9QJgpfvnL4HLjDGBxph2wPHA/LpuSEPoEhfO/sJSdmk5aRGRevHKK6/w/YY8Qk88mexf3iOpYyyPPvqo02V5HKfGNEwEJhpjVgIlwPADXRPGmC1ABBBgjDkfOMNauxr4G/A2EEzFrIkDMyfGG2N6UdE9sQW4GcBau8oYMwVYDZQBt3rbzIkDDsyaWL0zh9ZRwQ5XIyLSuCxatIh7xj1P9OXjKdy8mMBNv/DR0iX4+Tn1V6TncuSfiLW2BLjqEMcSD7F/IX+chom19pDrd1prxwHjalel5zihVQTGVMygOK1LrNPliIg0Gvv37+eSy68k8uxR2OICsr5+nm8+m0Tr1q2dLs0jaUVILxAW6Efb5iGs0bRLEZE6Y63lxhtvZH/70wlo0ZbMr5/l3jtGahzDYSg0eIkurSO0VoOISB169dVX+XrFbsJ7ncX+OR/TLyGcxx57zOmyPJpCg5fo3CqCrVkF5BWXOV2KiIjXW7x4MaMffZLos/5O8Y61+K2ZzkcffaRxDEeg0OAlDgyGXKunDSIixyQnJ4dLhl1O5Nl3Yq0l48vxvPfu27RpU+36f1KFQoOXqJxBodAgIlJrB8Yx7Es4mcC4TmR9+yL33Ho9Z511ltOleQWFBi/RKiKIqBB/DYYUETkGr732GtMWbSai/wXkLPqKPi19efzxx50uy2uo88ZLGGPoEhfB6p0KDSIitbFkyRJGPTiWmKuepWTPRnyWT+WjRQs0juEo6EmDF+kcF8Ha3bmUlbucLkVExKtUjGO4jMiz7sD4+pPxxVO89/ZE4uPjj3yxVFJo8CJd4iIoLnOxJSvf6VJERLyGtZabbrqJrFb9CUroxt7v/s3om6/h7LPPdro0r6PQ4EUODIZcpS4KEZEae/311/liXioRAy8lb/n39GpWytixY50uyyspNHiRDi3C8Pc1rNmV63QpIiJeYenSpYx64FFi/nIXpVlpsOhjJk2apHEMtaTQ4EUC/HzoFBvOyh37nS5FRMTj5efnc8mlw4g483ZMQDCZXzzFu2/9V+MYjoFCg5fpER/Jih37cb8UVEREqmGt5dlnnyWjZR+C2vZk7w+vMer6yznnnHOcLs2rKTR4me5tothfWEra3kKnSxER8VhvvPEGszdmEpl8GXkrf6JHeIHGMdQBdep4mR7xkQAs35HNcdEhDlcjIuJ5rLV8PeNXWvz1Hsr27cQu+IjJC+fh7+/vdGleT08avEyn2HACfH1YsV3jGkREquOyEHb63/EPDWfvtKd5d+J/SEhIcLqsRkFPGrxMgJ8PnePCWbY92+lSREQ80oSfNjB7YxbXdQ9hxPyZHHfccU6X1GjoSYMX6h4fycodObhcGgwpIlLV7A2ZvDBjHRee1IaT2/gpMNQxhQYv1CM+irziMjZrZUgRkUrpuUXcMXkp7WNCGXt+N4wxTpfU6Cg0eKEDgyE1rkFEpEK5y3LHR0vJLSrllSv7EBqo3vf6oNDghTq2CCPI34flCg0iIgC8OGM9czZl8fjQbpzQKtzpchothQYv5OfrQ9fWkazYocGQIiK/rc/k5Z/Wc1HveC7tq1kS9UmhwUt1b1MxGFKvyRaRpiw9p4h/TF5CxxZhjD2/q9PlNHoKDV6qR3wkhaXlbMzQYEgRaZrKyl3c9tES8ovLeeXK3oQEaBxDfVNo8FI94qMAWK71GkSkiXr+x3XM27yXcRd04/hYjWNoCAoNXqp9TCihAb6s0BsvRaQJmpmazr9nbuSyfglc2FtvrWwoCg1eysfH0K1NpGZQiEiTsyO7kDsnL6VzXASPnqdxDA1JocGL9YiPZPWuHEo1GFJEmoiSMhd//3AxZeWWV67sTZC/r9MlNSkKDV6se3wUJWUuUnfnOl2KiEiDeGr6WpZsy2b8xT1oFxPqdDlNjkKDF+vRxr0ypMY1iEgTMH3lLt78bTMjkhM5p3uc0+U0SQoNXqxtdAgRQX6aQSEijd7mzHzu/ng5PROieOCczk6X02QpNHgxYww9E6JYmqYnDSLSeBWVljPy/UX4+hpeubI3AX76q8sp+ifv5U5KiCJ1dw55xWVOlyIiUi8emrqS1D25vDCsF22igp0up0lTaPByvds2w2VheZq6KESk8ZmyII2PF23ntiEdGXxCS6fLafIUGrzcSQnNAFi8bZ/DlYiI1K1VO/fz0Bcr+VPHGO44rZPT5QgKDV4vMsSfji3DWLRVoUFEGo/9haWMfH8xzUICeOGyXvj6GKdLEhQaGoXex0WxJC0ba63TpYiIHDOXy3LXlGXszC7k31eeRExYoNMliZtCQyPQ+7hmZBeUsilTb7wUEe/3+i+b+HHNHh44pzN92jZ3uhypQqGhEejd1j2uQV0UIuLl5mzM4unv1nJujziuHZTodDnyOwoNjUDHFmGEB/mxeJtmUIiI99qTU8RtHy2mXUwoT13UA2M0jsHTKDQ0Aj4+hpOOa8YSzaAQES9VWu7i1g8Wk19czqtX9SEs0M/pkqQaCg2NRO/jokjdk0thmQZDioj3eeKbNSzcuo+nLu5Bp9hwp8uRQ3AsNBhjbjPGrDXGrDLGjHfvizbGzDTG5BljJvzu/HHGmDRjTN7v9gcaYyYbYzYYY+YZYxKrHLvfvT/VGHNmQ7TLKb2Pa4a1sDFbr8kWEe/y5bKdvDVrC9cOSuS8nq2dLkcOw5HQYIwZAgwFelpruwLPuA8VAQ8Bo6u5bBrQv5r91wP7rLUdgeeBp9zf0QW4DOgKnAW8YoxptC9e37dxCWD59+RveO+995wuR0SkRtbtyeXeT5bTt20zvYjKCzj1pGEk8KS1thjAWpvu/jXfWvsbFeHhINbaudbaXdXcayjwjvvnT4BTTcXomaHAJGttsbV2M7CB6kNHo/DI/fdQkrGNvSaCGTNmOF2OiMgR5RaVcst7iwgN9OPfV/bG31c95p7Oqd+hTsDJ7u6En40x/Y7hXm2ANABrbRmwH4iuut9tu3tfozRo0CCKd6whoPWJzJo9x+lyREQOy1rL3R8vZ+veAiZccRKxEUFOlyQ1UG/DU40xPwKtqjk0xv29zYEkoB8wxRjT3jq0pKEx5ibgJoDY2FhSUlKcKOOYREVFUbxzPuG9zmLL3gKmTp1KVFSU02XVmby8PK/8fTkStcu7qF1155tNJUxfV8plJwRQtG0FKdvq/jv0+1X36i00WGtPO9QxY8xI4DN3SJhvjHEBMUBGLb5qB5AAbDfG+AGRQFaV/QfEu/dVV+sbwBsAffv2tYMHD65FGc5q27Yt4197F4DA1p3x8fHBG9txKCkpKY2qPQeoXd5F7aobszdk8sl38zi3exz/uuKkeluPQb9fdc+p7ompwBAAY0wnIADIrOW9vgSGu3++GPjJHUa+BC5zz65oBxwPzD+mqj1YYmIiMQHllBfmEtjmRGbPnu10SSIif7Azu5DbPlpCTKCLR87uoAWcvIxToWEi0N4YsxKYBAw/0DVhjNkCPAeMMMZsd8+CwBgz3hizHQhx73/Ufa83gWhjzAZgFHAfgLV2FTAFWA1MB2611pY3VAMbmjGG5OSBFO9cS1B8F4UGEfE4xWXl3Pj2XPbl5BG3cRotm0c6XZIcJUeW3LLWlgBXHeJY4iH23wPcU83+IuCSQ1wzDhhX60K9THJyMj9+PIeQDv1YuDKVkpISAgICnC5LRARrLVc+8zmrsoNJ//Jp3n/rmSNfJB5H81sakeTkZIrSVgJgWh7P0qVLHa5IRATS0tJIvuouFmYHs3/uJ3QJL2HgwIFOlyW1oNDQiPTu3RubtRVXSRGB8V3VRSEijrLW8sYbb9Bj8F/ZGXcyhVuWkP3Lu9x+++0ay+ClFBoakcDAQDp17EDxzjUEHdddoUFEHLN582ZOO+00Rv7jbkLPvIPy/Gwyv3yali1iGDZsmNPlSS0pNDQy3bp1ozhtFf4t2jJrwRIcWvpCRJool8vFyy+/TLdu3fhpZgox592Db0gUGZ+Pw1WYw80330xgYKDTZUotKTQ0Ml27dqVo2wqM8WGfX3PS0tKOfJGISB1Yv349p5xyCrfffjsFBQVEnTKc4MReZH3/CiV7NuLn58ctt9zidJlyDBQaGpmuXbtSvGsdtqyEoAR1UYghs6E/AAAgAElEQVRI/SsvL+fZZ5+lR48e/PbbbwCEnHgykQMuInfJN+Sv+BGAiy++mNat9RZLb6bQ0Mg0b96c9m0TKN6ZSmCCBkOKSP3at28fp556KqNHj6aoqOJdg/4t2hF99h0UbV/F3h/fqDz39ttvd6pMqSMKDY3QgamXAbEd+G3eIqfLEZFGrFmzZkyfPp1ff/2VJ554gvDoWFpcOAZXcT4ZU/8FrjIA+vbtS1JSksPVyrFSaGiEkpOTKU5bifHxJTWrhPz8fKdLEpFGLCgoiD/96U8EBAYRfOqt+IVFVwx8zM+uPEfTLBsHhYZGKDk5meKda7HlZQS06cqCBQucLklEGrmvv/6acdOWE9yuN3t/eAWbubnyWMuWLbn00ksdrE7qikJDI9StWzdCA/0p3r1e4xpEpN6tXLmSax58kYgBF5O75BsKV/3EtGnT6NChAwC33HKLplk2EgoNjZCvry9JSUkUp60kMO54fpvTaF/uKSIOS09P5y9X3UzokJsrBz6+9NJLnHnmmdx+++2aZtnIKDQ0UgcGQxpffxZuycLlcjldkog0MsXFxZx3yRWUJl2HqziPjKn/4taRN/O3v/0NgBEjRnDDDTcQFxfncKVSVxQaGqnk5GSKt6/GusopiTqOdevWOV2SiDQi1lpuvOlmNrc6xT3w8QlOG9SfF154ofKciIiIg7bF+yk0NFJJSUlQWkTJnk1a5ElE6txTTz3FtG0+BLfrTdb3r9A+0ofJkyfj5+d30Hkay9C4KDQ0UpGRkXTr1o2ibcsJbH0iv8ya63RJItJIfP755/zzve8qVnxc/BWBOxYxbdo0oqKinC5N6plCQyOWnJxM0ZalGD9/Zm9Id7ocEWkElixZwog7H6L5WbdRtG0FOT+/xaeffkrHjh2dLk0agEJDI1Y5rqGslAzfaLKyspwuSUS82K5duzjv0isJP2c0roL9ZEz9F6/+ewKDBw92ujRpIAoNjVhycjK2rJii7asJTjyJuXPVRSEitVNYWMjQCy6kdMAIfILDSf9sLP8YeQM33HCD06VJA1JoaMQ6dOhAixYtKNq6lICW7ZgxSytDisjRs9Zy7XXXsTGqD0HxXcn65kXO6NeF8ePHO12aNDCFhkbMGOMe17AEgF/X7XG4IhHxRmPHjuXr1FzCe53N/jlTaO+3jw8//BBfX1+nS5MGVuvQYIz5R10WIvUjOTmZkj2bKC/MZVtJKKWlpU6XJCJeZMqUKfzrzU9pftpNFKyfh/+a6UybNo2IiAinSxMHHMuThlF1VoXUm+TkZLAuirYuwz++G8uWLXO6JBHxEgsWLOC62+4m5vz7KN27nZzvX+Lzzz8jMTHR6dLEIccSGvSOUy/Qp08f/P39KdqyBL+IFnz1y0KnSxIRL7B9+3aGXnQpEX+5FzBkfDqW/7zyMoMGDXK6NHHQsYQGW2dVSL0JDg6md+/eFG1ZCkDK2l0OVyQini4/P5+/njeUsn5X4R8dT+YXT3L3yGu5+uqrnS5NHHbY0GCMyTXG5FTzyQXaNFCNcoySk5Mp27+H0n272FygJV1F5NBcLhfXXHMNW8K7EXJ8Evtm/IezTmrHP//5T6dLEw9w2NBgrQ231kZU8wm31mrYrJdITk4GoGjrUmyL49myLc3hikTEUz300EN8t3YvkQMvJXfpt7R3bee9997Dx0eT7eTYZk9sq8tCpP5UhoYtS/EJDOHjGfMdrkhEPNH777/PM299QvTZt1O0bQUBK6Yy7csvCQsLc7o08RAaCNkEtG7dmrZt21K0dTnWuvhp1Q6nSxIRDzN79mxuvOMeWlz4IGV5e8n59jm++PwzEhISnC5NPIgGQjYRycnJuIpyKdm9gQ15/k6XIyIeZPfu3Zx/0aVE/fVefPyDyPj0cd567WX69+/vdGniYfwOd9AYc6i1GAyg51VeJDk5mY8++oiizUsISLqY3XtzaNVci7OINHW5ubk8MOZB7ICrCYhtT8anY3ng1msZNmyY06WJBzrSk4bwQ3zCgBfrtzSpSwfGNRRuWoDx8eW9H7Reg0hT53K5uOKKK9jbOonQE//EvpkTOa9vex555BGnSxMPddgnDdbaxxqqEKlfPXr0ICQkhIKd6ygv2M8Pq/Zxt9NFiYijfHx8SDzlYlZkxpC77DtOYAdvvfUzxmjImlTvSN0TDx/msLXWjq3jeqSe+Pn5MWDAAGbOnEnh5sVsPCGJcpfF10f/cRBpqpZs28d32S1p7Z9Hzo5f+PK7bwkODna6LPFgR+qeyK/mA3A9cG891iX14MDyr4UbF1DuF8yytH0OVyQiTtmRXciN7y6iVUQQ958cw/Kli4mLi3O6LPFwR1rc6dkDH+ANIBi4FpgEtG+A+qQOVa7XsHkx1lXOp3NSHa5IRJyQV1zG9W8voLi0nDeH9yU8wGjxJqmRI/4pMcY0N8b8E1hORXdGb2vtvdba9HqvTupUUlISAK6iPIp3rOWnNfotFGlqyl2W2z9awvr0PP59ZW+Ojw13uiTxIkd698TTwAIgF+hurX3UWqtn2l6qWbNmdOnSBajoothV7MeenCKHqxKRhvTEN2v4aW06j57Xlf/r1MLpcsTLHOlJw11Aa+BBYGfVF1YZY3Lqvzypa5VTLzcuAGDmWj1tEGkqPpi3lTd/28y1gxK5Oqmt0+WIFzrSmAYfa21wNS+uCrfWamUgL3QgNJRmbqUsJ4PvVmhJaZGm4Jd1GTz8xSqGnNCCB8/t4nQ54qUOO+VSGp8DoQEqnjbMjmpBcVk5gX56aalIY5W6O5dbP1jM8S3DePmK3ppqLbWm4bJNTKdOnWjevDlQERqKXTB/816HqxKR+pKRW8x1by8gOMCXiSP6ERao/1eU2lNoaGKMMf+berl1OcZVxk8a1yDSKBWWlHPDuwvZm1/Cm8P70TpKCzfJsXEsNBhjbjPGrDXGrDLGjHfvizbGzDTG5BljJvzu/HHGmDRjTN7v9o8wxmQYY5a6PzdUOTbcGLPe/RneMC3zfAdCgy0rpjhtBTPW7MFavbRUpDFxuSx3fbyU5duzefGyXnSPj3S6JGkEHHlOZYwZAgwFelpri40xLd2HioCHgG7uT1XTgAnA+mpuOdla+/fffUdz4BGgLxWv8V5kjPlSU0YPHteQu3Y229qexPr0PDppvrZIo/HUd2v5ZsVuHjy3M2d0beV0OdJIOPWkYSTwpLW2GODAQlHW2nxr7W9UhIeDWGvnWmt3HcV3nAn8YK3d6w4KPwBnHXvp3q9fv374+lYMfCxcPxewfLtit7NFiUideX/uVl7/eRNXJ7Xl+j+1c7ocaUScGhHTCTjZGDOOioAw2lq74Bjud5Ex5v+AdcCd1to0oA2QVuWc7e59f2CMuQm4CSA2NpaUlJRjKMVZeXl5Naq/Y8eOpKamUp6/j6C8nXwyz5eefp47/bKm7fI2apd38YZ2Lcso44VFxfRs4cvgiAx+/vnnI17jDe2qDbWr7tVbaDDG/AhU90xsjPt7mwNJQD9gijGmva1dx/o04CN3N8fNwDvAn4/mBtbaN6h4twZ9+/a1gwcPrkUZniElJYWa1H/mmWeSmlrx7oniDfNIC2tDu+79aBsdWs8V1k5N2+Vt1C7v4untWrljP6//NIcurSP48OaBhNZwpoSnt6u21K66V2/dE9ba06y13ar5fEHF//V/ZivMB1xATC2/J+tANwfwX6CP++cdQEKVU+Pd+4SDxzWkzfkKgO9WqYtCxFvtzC7k+ncWEBXsz8QR/WocGESOhlNjGqYCQwCMMZ2AACCzNjcyxlR9l+t5wBr3z98BZxhjmhljmgFnuPcJB4eG8px04kNcTF+p0CDijXKKSrn2rQUUFJcz8dp+xEYEOV2SNFJOhYaJQHtjzEoqXrM9/EDXhDFmC/AcMMIYs90Y08W9f7wxZjsQ4t7/qPtet7unbS4DbgdGAFhr9wJjqXjh1gLgcfc+ARISEoiPj6/cjszdwuJt2ezerxdYiXiTkjIXI99fxMaMPF67ug8nttIK/1J/HAkN1toSa+1V7u6K3tban6ocS7TWNrfWhllr4621q93773Fv+7h/fdS9/35rbVdrbU9r7RBr7doq95pore3o/rzV4A31cFWfNqQv/h6A71fraYOIt7DWct9ny5m1IYunLurBoI616uUVqTGtCNmEVQ0Ny3/7ng4xIeqiEPEiz/+4ns8W72DU6Z24qE/8kS8QOUYKDU1Y1dBQUlJC92Yu5m3ey978EgerEpGamDR/Gy/NWM+lfeO57c8dnS5HmgiFhiasV69eBAf/by36wIzVlLssP67Z42BVInIkP63dw5ipKzmlUwvGXdAdY/TWSmkYCg1NmL+/P/3796/cXjd/JvHNgvlOXRQiHmtZWja3frCEznHhvHJlb/x99Z9xaTj609bEVe2imDN7Nmd2bcWv6zPZX1jqYFUiUp0tmflc9/YCosMCtBaDOEKhoYmrGhr27NlDnxhLSblLTxtEPExWXjEj3pqPy1reua4/LcO1FoM0PIWGJi4pKemg7az1S0iMDuGLZVo8U8RT5BeXcd3bC9idU8R/h/ejQ4swp0uSJkqhoYmLiYnhhBNOqNyeM2c25/Vqw+yNWaTnaKEnEaeVlrv42weLWbFjPxMu702fts2cLkmaMIUGOaiLYvbs2Qzt1RprYdryo3kTuYjUNWst9326gp/XZfDEBd05rUus0yVJE6fQIAeFhhUrVtAi0EX3NpF8sVRdFCJOevq7VD5dvJ07T+vEZf2Pc7ocEYUGOTg0uFwu5s2bx9BerVm+fT+bMvIcrEyk6Xpr1mZeSdnIFQOO4/ZTtXiTeAaFBuHEE08kKiqqcnv27Nn8pUdrjIEvl+10sDKRpumLpTt4bNpqzuway+PnddXiTeIxFBoEHx8fBg4cWLk9e/ZsWkUGkdQumi+W7sT9AlIRaQA/r8vgrinLGNCuOS9edhJ+WrxJPIj+NApwcBfF3LlzKS8vZ2iv1mzOzGfFjv0OVibSdCxNy2bk+4tIbBbIvcmRBPn7Ol2SyEEUGgQ4ODTk5OSwevVqzu4Wh7+v4Yul6qIQqW8b0vMY/uY8KMpl23v30uX49k6XJPIHCg0CQP/+/fHx+d8fh9mzZxMZ4s/gE1oybdlOyl3qohCpL2u27eG8575nX1Ym6/9zB4/ceydBQVrxUTyPQoMAEBYWRs+ePSu3Z8+eDcCFJ7UhPbeYX9ZnOFWaSKNVUFDAo08+yxlPTCOvqJTdkx+iR/s4Lr/8cqdLE6mWQoNU+v0iTwCndo6leWgAUxakOVWWSKNTWlrKa6+9RsfO3Xl9rR8mtBnpnzxGacZmnnnmmYOe+ol4Ev3JlEpVQ8OGDRtIT08nwM+H83u14cc1e8jKK3awOhHv53K5+OCDDzjxxBP522134PrTTQS0SCTj8yco3rGGc889lyFDhjhdpsghKTRIpaqhAWDOnDkADOuXQGm55fMlWiFSpDastUybNo1evXpx1VVXsWnLNmLOv5/A+C5kfv0cRZsX4+Pjw/jx450uVeSwFBqkUtu2bYmLi6vcPtBFcUKrcHomRDFlYZrWbBA5SikpKQwaNIjzzjuPFStWgPEh5tw7CenQj73fvULBml8AuP766+nSpYvD1YocnkKDVDLGVDuuAWBY3wTW7cljaVq2E6WJeKWJEydyzjnnVD61A2h++khCu5zCvplvkbdsOgChoaE89thjTpUpUmMKDXKQQYMGVf68YMECSkpKAPhrzziC/X2ZslADIkVq6rrrriM9PZ0PP/yQnj17EnXKcMJPOpv9c6aQM//TyvNGjx590FM+EU+l0CAHqfqkobi4mCVLlgAQHuTPOd3jmLZsFwUlZU6VJ+J1wsLCaNOmDTsiuxGZdAm5i78m+5d3K4+3atWK0aNHO1ihSM0pNMhBTjrpJAIDAyu3D+qi6JdAXnEZ36zY7URpIl5p5syZXHjfS4QOvJy8VTPZ+8NrBx1/7LHHCAsLc6g6kaOj0CAHCQgIoF+/fpXbVUNDv8RmtI8J1ZoNIjX0448/ctE9zxL2f9dSkDqbrK+f5x//uIOYmBgAOnfuzHXXXedwlSI1p9Agf/D7wZAHZkwYY7ikbwLzt+xlQ3quU+WJeIXp06dz8V1PEv7nmynctJCMaeO55+7RPPfccwwdOhSA8ePH4+fn53ClIjWn0CB/UDU07Ny5k23btlVuX9I3ngBfH96ds9WJ0kS8wldffcWlo8YReebtFG9bScbnTzDmvnt58sknMcZw4YUXMnjwYM4991ynSxU5KgoN8gcDBw48aLtqF0VMWCB/6RnHp4u2k1tU2tCliXi8qVOncvld42h27l2U7FpP+mdjeeTBBxg7dizGGABOPfVUXnzxxcptEW+h0CB/0LJlSzp27Fi5XTU0AAwfmEh+STmfLtre0KWJeLRPPvmEq0aPo/lf76EkfQt7Pn6EsQ+P4ZFHHjkoIAQGBtKjRw8HKxWpHYUGqdahFnkC6JkQRa+EKN6dsxWXXpktAsDkyZMZfvc4mg+9n9Ks7aRPeYgnH3+YMWPGOF2aSJ1RaJBqVQ0Ny5YtIy8v76DjI5IT2ZSZz68bMhu6NBGP88EHHzDi7rFEX/AgZdm72TP5QZ4e9xj33HOP06WJ1CmFBqlW1dBQXl7OggULDjp+Tvc4YsICeWf2lgauTMSzvPPOO1x/3xO0uOgRynMz2DN5DM8/OZZRo0Y5XZpInVNokGp16dKFiIiIyu3fd1EE+PlwRf8EZqamszUrv6HLE/EIb775JjePeZKWFz9Ked5e9nz0ABOefoLbb7/d6dJE6oVCg1TL19eXpKSkyu1Zs2b94Zwrk9riawzvafqlNEGvvfYaf3v4aVpe8jjl+dmkTx7Dq88/xciRI50uTaTeKDTIIVXtopgzZw4ul+ug47ERQZzVrRWTF6bpfRTSpEyYMIE7Hn+e2Esfp7wgm/RJD/DGi09z4403Ol2aSL1SaJBDqhoasrOzWbt27R/OuXZQIrlFZUzW0tLSRDz//PPc9cTLxA4bS3lBDhmTH2Tiv5/j2muvdbo0kXqn0CCHNGDAgIPmlv9+XANAn7bN6ZfYjP/+upnSctcfjos0JpMmTeK+Z14ndtg/KS/YT+aUB3nn1Re4+uqrnS5NpEEoNMghRURE0L1798rt6kIDwC2ndGBHdiHTlu1sqNJEGtwTTzzBW9N+JvbSsZTnZ5M55SHef+NlLr/8cqdLE2kwCg1yWIdb5OmAISe05ITYcF77eaMWe5JGyeVyMWvd7ooxDPn7yPz4QSZNfJVLLrnE6dJEGpRCgxxW1dCQmppKZuYfF3Py8THcfEp71u3JY2ZqekOWJ9Ig5m3ex+bjziHAVUTWJw/z8Tv/4YILLnC6LJEGp9Agh1U1NADMnTu32vP+2rM1baKCeTVlY0OUJdJgfl2fwbVvz6dNVDBPntGKX6Z/yXnnned0WSKOcCw0GGNuM8asNcasMsaMd++LNsbMNMbkGWMmVDk3xBjzdZXzn6xyLNAYM9kYs8EYM88Yk1jl2P3u/anGmDMbsn2NRfv27YmNja3cPlQXhb+vDzec3I6FW/exYMvehipPpF79tHYP17+zkHYxYUy6KYmYUH/69+/vdFkijnEkNBhjhgBDgZ7W2q7AM+5DRcBDwOhqLnvGWnsicBIwyBhztnv/9cA+a21H4HngKfd3dAEuA7oCZwGvGGN866lJjZYxpkbjGgCG9UugWYg/r+lpgzQC01fu5ub3FnFCbDgf3TiA6LBAp0sScZxTTxpGAk9aa4sBrLXp7l/zrbW/UREeKllrC6y1M90/lwCLgXj34aHAO+6fPwFONRXzBIcCk6y1xdbazcAGQP+LUAtVQ8P8+fMpLS2t9ryQAD+GJycyY206a3blNFR5InXu8yXbufXDxXRvE8kHNw4gKiTA6ZJEPIJToaETcLK7O+FnY0y/ml5ojIkC/grMcO9qA6QBWGvLgP1AdNX9btvd++QoVQ0NhYWFLFu27JDnjkhOJDzQjxd+XNcQpYnUuffnbmXUlGUMaNec964fQESQv9MliXgMv/q6sTHmR6BVNYfGuL+3OZAE9AOmGGPaW2sPO1/PGOMHfAS8ZK3dVIe13gTcBBAbG0tKSkpd3brB5eXl1Xn9JSUl+Pv7Vz5heOedd/7wquyqTk0wTF21h7e/mEFiZN30CNVHuzyB2uVZvt1cyuTUEnq28GVE+0IWzPntoOPe2q4jUbu8i6PtstY2+AeYDgypsr0RaFFlewQwoZrrJlIRGKru+w4Y6P7ZD8gEDHA/cH915x3u06dPH+vNZs6cWS/3HThwoAUsYIcNG3bYc/cXltgej35nR0ycV2ffX1/tcpra5RlcLpd99vtU2/ber+zfPlhkS8rKqz3P29pVU2qXd6mPdgELbQ3+/naqe2IqMATAGNMJCKDiL/tDMsb8E4gE/vG7Q18Cw90/Xwz85P4H8CVwmXt2RTvgeGB+nbWgianpYEiAiCB/bj6lPTNTM1i0dV99lyZyTFwuy2PTVvPSjPVc0ieely47CX9fzUYXqY5T/2ZMBNobY1YCk4Dh7r/oMcZsAZ4DRhhjthtjuhhj4qno1ugCLDbGLDXG3OC+15tAtDFmAzAKuA/AWrsKmAKspuLJxq3W2vIGa2EjUzU0pKWlkZZ2+BdUjUhOJCYsgOd+SK3v0kRqrbTcxV0fL+Pt2Vu44U/tGH9xD3x9zJEvFGmi6m1Mw+HYihkQVx3iWOIhLqv232RrbRFQ7Vqu1tpxwLhalCi/M3DgwIO258yZQ0JCwiHPDwnw45ZTOvDPr9cwe2MmyR1i6rtEkaNSVFrOrR8sZsbadO4+8wT+NrjDQS9oE5E/0jM4qZG4uDjatWtXuX2kLgqAq5LaEhsRyHPfrzswrkTEI+wvLOWaifP5KTWdsed349YhHRUYRGpAoUFq7GjGNQAE+fvy9z8fz8Kt+/ROCvEY6TlFDHt9Dku27eOFYb24Oqmt0yWJeA2FBqmxqqFhyZIlFBQUHPGaYX0TaBcTyriv11Ba7qrP8kSOaHNmPhe9Npttewt4c3g/hvbS0i0iR0OhQWqsamgoKytj4cKFR7wmwM+HB87pzMaMfD6ct60+yxM5rBXb93Pxq7PJLy7noxuT+L9OLZwuScTrKDRIjXXr1o2wsLDK7Zp0UQCc1rklyR2ief7HdewvqH4JapH69PO6DC57Yw5B/r58fMtAeiZEOV2SiFdSaJAa8/PzY8CAAZXbNQ0NxhgePLcL+wtLeXHG+voqT6RaHy9M47q3F3BcdCif/S2ZDi3CjnyRiFRLoUGOyu8HQ9Z0VkSX1hFc1i+Bd+dsYVPGoZegFqkr1lpemrGeuz9ZTnKHaKbcnERsRJDTZYl4NYUGOSpVQ0NWVhbr1tX8xVSjTj+BIH9fnvhmTX2UJlKptNzFA5+v4Lkf1nFh7za8Obwf4XrxlMgxU2iQo5KUlHTQdk27KABahAdy65CO/LgmnV/WZdR1aSIA5BaVct3bC/hofhp/H9KRZy/pSYCf/lMnUhf0b5IclaioKLp27Vq5fTShAeDaQYm0iwnloS9WUlSqVb2lbu3MLuSS1+YwZ2MW4y/qwegzT9CiTSJ1SKFBjtrRLvJUVZC/L+PO78bWrAJe0qBIqUMrtu/n/H/PYse+Qt6+tj+X9jv0MuciUjsKDXLUBg0aVPnz6tWr2bfv6N5kmdwxhot6x/PGL5tI3Z1b1+VJEzR95W4ufX0O/r4+fPq3ZP50vN51IlIfFBrkqFV90gAwd+7co77HmHM7Ex7kxwOfr8Dl0nsppHastfx75gZueX8RJ8aFM/XWQXSKDXe6LJFGS6FBjlrHjh2Jifnf/8kdbRcFQPPQAB48twuLtu7jw/laKVKOXlFpOaOmLOPp71IZ2qs1H92YRIvwQKfLEmnUFBrkqBljjmlcwwEX9m7DoI7RPPXtWvbkFNVVedIEpOcWceV/5/H5kh2MPqMTLwzrRZC/r9NliTR6Cg1SK1VDw7x58ygrKzvqexhj+Of53Skpd3Hfp8v1+mypkeXbsznv5Vms3pnDK1f25u9/Pl4zJEQaiEKD1ErV0JCfn8+KFStqdZ92MaHcf/aJzEzN4AO90EqO4PMl27nktTn4+hg+GTmQc7rHOV2SSJOi0CC10rdvX/z8/Cq3a9tFAXDNwEROPj6GcV+v0RLTUq2ychdPfLOGOycvo1dCFF/+fRBdW0c6XZZIk6PQILUSHBxM7969K7ePJTT4+Bievrhi1b47pyyjtNxVFyVKI5GZV8w1E+fzxi+buDqpLe/fMIDoMA14FHGCQoPUWl0MhjygVWQQT1zQnWVp2Uz4acOxliaNxJJt+/jry7+xaOs+nr64B2PP74a/r/6zJeIU/dsntVY1NGzZsoWdO3ce0/3O7RHHhSe1YcLMDSzednQLRknjYq3lg3lbGfb6XHCVc3uXUi7uE+90WSJNnkKD1NrAgQMP2p4zZ84x3/PRoV2Jiwzi7x8sJiuv+JjvJ95nb04+l780nTGfr6R812qWP3MV/Tq20gwJEQ+g0CC1Fh8fz3HHHVe5faxdFAARQf68dlUfMvNLuGPSUsq1WmSTsHXrVl577TXOvGQEPUZ/wJwdpWT/+j6b37mX9ye+Tv/+/Z0uUURQaJBjVJfjGg7o1iaSfw7txm8bMnnuh9Q6uad4luLiYn744QdGjRpF586dSUxMZPSEj1lz3F8xASGkT36I/bMn8fT4p7jgggucLldE3PyOfIrIoSUnJzNp0iQAFi1aRFFREUFBQcd830v7JbB42z7+PXMjPeOjCDjmO4rTNm/ezLfffsu3337LTz/9REFBAQAmIJjoc+4krPupFG1bQeaX4ynP38ctt9zCqFGjHK5aRKpSaJBjUvVJQ2lpKYsWLTroLZjH4tHzurJqZw53TVnGecFrKS0t5fTTT6+Te0vDmjVrFhdffDG7d+8+aH9AbAdizrsXv89IaPQAABpISURBVKhYsmd9yP5Zk8C6OOuss3j55Zc1jkHEw6h7Qo5Jjx49CAkJqdyuqy4KgCB/X165sjelJcW8vSmIYVdfx7p16+rs/tJwBg0aRGpqKqNGjcLX1xcwhPc7n1ZXP4Px82fPpDHs/+1DsC569OjB5MmTD1o8TEQ8g0KDHJP/b+/e46Oo7/2Pvz7kAhIIJNxv5SJCRRQql+KxKFAQCiikgHpaKvhDpYpgD/UCtZ4q9agPRX/aA1RRUaBQvPw05GCVooAUFFCQS4IoBBXBAAFCIECAJN/zxw78FkzMkOxmCLyfj8c+mJ3vzM7n89hl952Z2dm4uLjTTlKLZGgA2LxmBV/PnkhsYn3ie49l4KAUDhw4ENFtSMVITEwkJSWFpm3aU/+mSST3uo2jmZ+S9co4jn2bDkCjRo1YsGABiYmJAVcrIsVRaJByO/NkyEj+8FTv3r0Z3u8q9r7zDNWaXkZOm4HcdPPNZfqBLAlObm4ud911F31vn0hhnwlUbXwp+96bQvbb/0VR/iEAEhISWLBgAc2aNQu4WhEpiUKDlFt4aNizZw/btm2L2GObGdOmTePiuFxylrxMwqXXsPpYY+6///6IbUOia/78+bTr0InXt19EvUETKMjZSdar48hb/96pZapUqcK8efNOuzS5iJx7FBqk3M68yNOKFSsi+vjx8fE88sgj1N61hkNrF1Drp0N4adlWZsyYEdHtSGRlZWUxbNgwbh7/KDbgIRLaXcuB5XPIeeMh/nz/WIYOHXpq2WeffZaBAwcGWK2I+KHQIOWWnJzMpZdeeup+pM9rAEhKSiItLY1jH8/hyBcfkdx7NOOnvR3xgCLl55zjnXfeoV2HTiw50oQGwx6mKP8wu2b/np/EZbFxw3omTJhA8+bNARg3bhxjx44NuGoR8UOhQSIiGhd5OlOHDh2YPWsm2WlPcmTLKmr3/i3DJjzL9u3bo7I9KZtbb72VqW9/SMKwx0m4rBe5H73G0dQ/Me3RCSxevJhLLrkEgCZNmnD99dfzzDPPBFyxiPil0CARER4a0tPTyc3Njcp2UlJSmPTwf5I9/3GOZH5C/NUj6XvnI+Tl5UVle3J29hzM59DlN1J/6J8oys9j1+zf06dhPp9nbGTUqFGnXXehV69ezJ071/sKpohUBgoNEhHhocE5x6pVq6K2rT/+8Y/cOOSXZL/9GEe/WsuR9oMZMPZRioqKorZN+WFFRY65q7bz82c+ZGOOUWfnR1RZ9CSvPz+Z119/nYYNG35vnQ4dOlCjRo0AqhWRslJokIho06YNycnJp+5H6xAFhL5R8corr3Blh8vJfutR8r9exzcNujP0oekR/bqn+JO+M5chz3/EH97eSPvGtVj4u2t4+MZuZGzcwKBBg4IuT0QiSKFBIqJKlSqnfYsimqEBoHr16qSmplK/ThJ73pzE4U1LWVvYjOHPzNcvY1aQ3CMneCg1neunLOfb/UeYPKwDc2//KS3rJlCzZk1q1aoVdIkiEmEKDRIx4YcoVq5cSWFhYVS316xZM1JTU6kaF8Pe/3ma3NVvsSI7juHTFpN/IrrbvpAVFBYxd9V2ej69lDmrvmHEVS344Pc9GNqpqX4rQuQ8p9AgERMeGg4dOkRGRkbUt9mtWzdefPFFwHFgyQz2L36Zj3fkM2zav9h9MD/q27/QLPsymwF/Wc4f3t5I63o1WDC2Ow/fcBm1LooLujQRqQAKDRIxXbp0Oe1M+GgfojjpN7/5Dffddx8Ahz55m+zUx0nfkcPAv/yL1V/tr5Aazndf7j7EyFdWc8uM1Rw9Uchff30lr43uRrvG+o0IkQuJQoNETEJCAh07djx1Pzw0rFq1KqqHKx5//HEGDBgAwJEvVrDz1d+Rd2Afv3pxJTOWf6UTJMvo2/1HGP/6Ovo+u4w13+TwYP9LWTT+Gn5xeSMdihC5ACk0SEQVd5GnWbNm0a9fv6h+Hz8mJoa5c+eeujLlib3b+eK/R9Es9iCTFmxizNy17D98PGrbP99kHzrGw2kZ9Hp6Ke9syOL27q1Ydl9Pbr+mFVVjdV0FkQuVQoOU22effUZWVhZwemjIzMxk9OjRjBgxgurVq0e9jsTERNLS0khKSgLAHT/Csv8azuCWxqJNu7nu/37IPzN2Rb2Oymz3wXwe+Z8Muj+5mNkrv2Fop6Ysva8Hf+h/KUkJ8UGXJyIBCyw0mNlYM9tsZhlm9qQ3r46ZLTGzPDObErZsdTN7J2z5J8LGRppZtpmt8263hY2NMLMt3m1ExXZ44cjNzaVx48a0adOGuXPnnjY2ffp0gFMf5NHWunVr3nzzzVN7NYqKCnnl/n/nuf6NqV+zGnfMXsP419aRe+REhdRTWezIOcJDqel0f3IJsz7+hgGXN2bRf1zD47+8gka1Lgq6PBE5R8QGsVEz6wkMAjo4546ZWX1vKB94CGjv3cJNds4tMbN44AMz+4Vz7l1v7DXn3N1nbCMZ+BPQGXDAGjNLc87lRKmtC1aPHj0YPHgwqampbNmypdhlwi/8FG29evXiueee4+67Qy+JgwcPcs/IYSz/6GPmrNvP1CVbWbYlm/v6tmVop2bEVLlwj82v+/YAL/1rG++m76KKwdBOTbnz2tb8qE709wyJSOUT1J6GO4EnnHPHAJxze7x/DzvnlhMKD6c4544455Z408eBtUDTUrbRF1jknNvvBYVFQL/ItiEnPfXUU8TFlfy1u4oMDQB33XUXo0ePPnV/69at/Prfb2Zsj5bMH3M1zesk8MD/28gNU5ZfcN+wOFFYxDsbshj6148YPHUFH36ZzaifteTD+3ry+C+vUGAQkRJZEGeVm9k6YD6hD/F84F7n3Cdh4yOBzmfuPfDGahMKDb2dc9u8ZR8HsoEvgf9wzn1rZvcC1Zxzj3rrPQQcdc5NLuYx7wDuAGjQoEGnefPmRbLdCpWXlxfY9fynTZvGG2+8UexYv379eOCBB8r82GXpq6CggHvvvZf169efmpeSksK4ceNCv4+xq5DXvzjO/nxH5wYxDGodT7OaFZujK/L52nu0iA+/LWDZzgJyjznqXWRc1zyOnzWN5aLYyO5tCfJ1GE3qq3JRX/717NlzjXOuc6kLOueicgPeB9KLuQ3y/v1vwICuwFd4AcZbdyQwpZjHjAXeBX4XNq8OUNWbHg0s9qbvBf4YttxDhMLJD9bdqVMnV5ktWbIksG3v37/f1alTxxE6HHTabfz48eV67LL2lZ2d7Vq2bHlaLS+88MKp8a1fbXdDH57hLvvP91zzBxa4Ua9+4j7bnlOuWs9GtJ+vI8cKXOpnO9wtL69yLSYscC0nLHCjXl3tFn++2xUUFkVtu0G+DqNJfVUu6ss/4FPn47M9auc0OOd6lzRmZncCb3mFrjazIqAuob0FP2Q6sMU592zYdvaFjb8EPOlN7wR6hI01BZb6rV/OXlJSEg8//DBjx44tdiwIdevWJS0tjauuuurUz2ePGTOGtm3bcu211/LSC9NY+be/sW7Tl/xt1Q5mrPiKwVNXcHXrOvyqa3P6tGtAfGzl+pLR8YIiVm7bx/x13/FeehaHjxfSpPZFjO3Zmpu7/ojGtXVio4iUTSAnQgKpQE9giZm1AeKBvT+0gpk9CtQCbjtjfiPnXJZ39wbgc296IfCYmZ38tLoOmBiZ8qUko0ePZurUqWzevPm0+RV9TkO49u3bM2fOHAYPHoxzjoKCAoYMGcLSpUt54YUXyMnJYdn773FPSgqjurfkbyu/YdZHXzNm7lrqJMQzpFNTbuzclNb1awbWQ2nyjhWw7MtsFmbsYvHmPRzKL6Bm1VgGXtGYlCub0LVFMlUu4BM+RSQyggoNM4AZZpYOHAdGeHsdMLOvgUQg3swGE/qwPwg8CGwG1npXopvinHsJGGdmNwAFwH5ChzZwzu03sz8DJ8+VmOScu7DOeAtAXFwckydPZuDAgafNDzI0ANxwww089thjTJwYyo379u2jW7duHD58GIApU6aQkpJCjaqx/Pbai0MXM9qSzbzV25mx/CumL9tG6/o16NOuAde1a0CHprUD/RA+UVjE+m8PsHzrXj7auo+123MoKHIkVY+j32UNue6yhnS/pC7V4nQhJhGJnEBCgwt9A2J4CWMtSlit2Hdo59xEStiD4JybQSigSAXq378/ffr0YdGiRafmBXV4ItwDDzxAeno6c+bMATgVGAAWL17Mpk2baNeuHQAxVYyebevTs2199hzK5x8bslj0+W6mL9vGX5dmUrdGPF1aJNOpeRJdWiTTrnEicTHROYxRVOTYeeAoG3bksn7HAdZtP8DGnbkcPVGIGVzepBa3dW9Fj7b16Nw8idgo1SEiEtSeBjmPmRlPP/00HTt2pKioCAh+T8OxY8dYs2YNbdu2xcyK/S2KadOmMWXKlO/Nr1+zGiOvbsnIq1ty4Mhxlnyxhw+/yObTb3J4Nz10hcn42Cq0qpvAxfVrcHG9GlxcL4F6NapSp0ZVkhPiSaoed9qH+a5du5g4cSIHcnPpetXPOHy8kL15x9h9MJ89B4+RlZvPtr15ZGbnkbnnMEe9n/qOj6lCu8aJ3NSlGT9tmcxVF9ehdnVdqVFEKoZCg0TF5Zdfzm233VbhV4QsTmZmJn379iUzM/MHl5s5cyaPPfYYiYkl/3Jj7erxpPykKSk/CV0mZPfBfD79Oof1Ow6QuSeP9J25vLsxi6JivskcW8WIj61CXEwVjh7JI792f6x+Ndr9aWGx22pS+yJa1Uvg5q7JtK5fg/aNa3Fpo8RKd2KmiJw/FBokaiZNmsTf//53Dh06FOiehosvvpg1a9Ywfvx4Zswo+WhVXl4es2bNOnUlST8aJFZjwBWNGHBFo1Pz8k8UsiPnCHvzjrMv7zj7Dh9j/+HjHC8o4kRhEScKHRvSM1j6wSLciXzuGfNbWv2oCXVrVKVBYlXq16xGvZpVdT6CiJxzFBokaho0aMCDDz7IxIkTqVWrVqC11KpVi5dffpkhQ4Zw++2389133xW73NSpUxkzZky5fva5WlwMrevXpHX9kpdx17fj/S5JZGZm8tubri7ztkREKpL2c0pU3XPPPXTs2DGqP4t9Nvr3709GRgYjR44sdnzz5s0sXrw46nWYGX369OHHP/5x1LclIhIpCg0SVdWqVeP5558PuozT1K5dm1deeYUFCxbQqFGj741PnTo1gKpERM59Cg0SdV27dg26hGINGDCAjIwMbrnlltPmz58/n+3btwdUlYjIuUuhQS5oSUlJzJw5k7S0NBo2bAhAUVHRObd3RETkXKDQIAJcf/31ZGRkMHx46JpjL774Ivn5+aWsJSJyYVFoEPEkJycze/ZsUlNTiYmJKfFnvkVELlQKDSJnGDRoEBkZGeTk5ARdiojIOUWhQaQYderUYdy4cUGXISJyTlFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMSXQEODmY01s81mlmFmT3rz6pjZEjPLM7MpZyz/npmt95Z/3sxivPnJZrbIzLZ4/yZ5883M/mJmW81sg5ldWfFdioiInB8CCw1m1hMYBHRwzl0GTPaG8oGHgHuLWe1G51wHoD1QDxjmzZ8AfOCcuwT4wLsP8AvgEu92B/DXKLQiIiJyQQhyT8OdwBPOuWMAzrk93r+HnXPLCYWH0zjnDnqTsUA84Lz7g4CZ3vRMYHDY/FkuZCVQ28waRaMZERGR850550pfKhobNlsHzAf6EQoI9zrnPgkbHwl0ds7dfcZ6C4GuwLvAb5xzhWZ2wDlX2xs3IMc5V9vMFhAKJsu9sQ+AB5xzn57xmHcQ2hNBgwYNOs2bNy8qPVeEvLw8atSoEXQZEae+Khf1Vbmor8olGn317NlzjXOuc2nLxUZ0q2cws/eBhsUMPehtOxnoBnQBXjezVq6UFOOc62tm1YA5QC9g0RnjzszOKgk556YD0wE6d+7sevTocTarn1OWLl1KZa6/JOqrclFflYv6qlyC7CuqocE517ukMTO7E3jLCwmrzawIqAtk+3jcfDObT+jwwyJgt5k1cs5leYcf9niL7gSaha3a1JsnIiIiZynIcxpSgZ4AZtaG0DkKe0ta2MxqnDwfwcxigQHAZm84DRjhTY8gdNjj5PxbvG9RdANynXNZkW5ERETkQhDVPQ2lmAHMMLN04Dgw4uShCTP7GkgE4s1sMHAdsA9IM7OqhMLOEuB577GeIHR4YxTwDXCjN/8fQH9gK3AEuLUC+hIRETkvBRYanHPHgeEljLUoYbUuJSy/D/h5MfMdMKaMJYqIiEgYXRFSREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxxZxzQddwTjGzbOCboOsoh7rA3qCLiAL1Vbmor8pFfVUu0eiruXOuXmkLKTScZ8zsU+dc56DriDT1Vbmor8pFfVUuQfalwxMiIiLii0KDiIiI+KLQcP6ZHnQBUaK+Khf1Vbmor8olsL50ToOIiIj4oj0NIiIi4otCg4iIiPii0FBJmFk/M/vCzLaa2YRixqua2Wve+CozaxE2doWZfWxmGWa20cyqVWTtP6SsfZlZnJnN9Pr53MwmVnTtP8RHX9eY2VozKzCzoWeMjTCzLd5tRMVVXbqy9mVmHcNegxvM7KaKrbx05XnOvPFEM9thZlMqpmJ/yvla/JGZ/dP7P7Yp/H0laOXs60nvtfi5mf3FzKziKv9hPvoa7z0XG8zsAzNrHjYW/fcO55xu5/gNiAEygVZAPLAeaHfGMncBz3vTNwOvedOxwAagg3e/DhATdE8R6OtXwDxvujrwNdAi6J7Ooq8WwBXALGBo2PxkYJv3b5I3nRR0TxHoqw1wiTfdGMgCagfdUyR6Cxt/DpgLTAm6n0j1BSwF+njTNYDqQfcUgdfivwErvMeIAT4GegTd01n01fPk8wDcGfaeWCHvHdrTUDl0BbY657Y5544D84BBZywzCJjpTb8J/NxLz9cBG5xz6wGcc/ucc4UVVHdpytOXAxLMLBa4CDgOHKyYsktVal/Oua+dcxuAojPW7Qsscs7td87lAIuAfhVRtA9l7ss596Vzbos3/R2wByj16nMVqDzPGWbWCWgA/LMiij0LZe7LzNoBsc65Rd5yec65IxVUd2nK83w5oBqhD+WqQBywO/ol++KnryVhz8NKoKk3XSHvHQoNlUMT4Nuw+zu8ecUu45wrAHIJ7VVoAzgzW+jtqru/Aur1qzx9vQkcJvQX63ZgsnNuf7QL9slPX9FYN9oiUpuZdSX0hp0Zoboiocy9mVkV4Gng3ijUVV7lec7aAAfM7C0z+8zMnjKzmIhXWDZl7ss59zGwhNB7Rxaw0Dn3ecQrLJuz7WsU8G4Z1y0ThYbzXyzwM+DX3r8pZvbzYEuKiK5AIaFd3S2B35tZq2BLktKYWSNgNnCrc+57f7FXUncB/3DO7Qi6kAiLBboTCkNdCO0yHxlkQZFgZq2BSwn9hd4E6GVm3YOt6uyZ2XCgM/BURW5XoaFy2Ak0C7vf1JtX7DLeLvtawD5CaXOZc26vt0vrH8CVUa/Yn/L09SvgPefcCefcHkLHKM+Va8z76Ssa60ZbuWozs0TgHeBB59zKCNdWXuXp7SrgbjP7GpgM3GJmT0S2vDIrT187gHXervICIJXK9d5RkhRgpXe4JY/QX+pXRbi+svLVl5n1Bh4EbnDOHTubdctLoaFy+AS4xMxamlk8oRMC085YJg04ebbsUGCxC50dsxC43Myqex+61wKbKqju0pSnr+1ALwAzSwC6AZsrpOrS+emrJAuB68wsycySCJ2TsjBKdZ6tMvflLf82MMs592YUayyrMvfmnPu1c+5HzrkWhP4qn+Wc+95Z7wEpz2vxE6C2mZ0896QXleu9oyTbgWvNLNbM4gi9J54rhydK7cvMfgK8QCgw7Akbqpj3jqDPFtXN91m1/YEvCR0HftCbN8l74UDoxJ43gK3AaqBV2LrDgQwgHXgy6F4i0RehM7nf8PraBNwXdC9n2VcXQn/JHSa05yQjbN3/4/W7ldBu/MD7KW9f3mvwBLAu7NYx6H4i9ZyFPcZIzqFvT0TgtdiH0LevNgKvAvFB9xOB12IMoQ/dz733jmeC7uUs+3qf0ImbJ/8fpYWtG/X3Dl1GWkRERHzR4QkRERHxRaFBREREfFFoEBEREV8UGkRERMQXhQYRERHxRaFBREREfFFoEBEREV8UGkQkUGbWwcyWmdkmMysyM2dmk4KuS0S+Txd3EpHAmFk1Qle1u8U5t9rM/kzoKqD3O705iZxztKdBRILUG1jrnFvt3d8AJCswiJybFBpEJEjtCf2uwUlXAmsDqkVEShEbdAEickHbx///tdI2wC+Bfwu0IhEpkc5pEJHAmFkN4O9AS2AvMN45pz0NIucohQYRERHxRec0iIiIiC8KDSIiIuKLQoOIiIj4otAgIiIivig0iIiIiC8KDSIiIuKLQoOIiIj48r90u9Oe+ny/GgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11f4da390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (8,8))\n",
    "plt.plot(sigma_grid, nll_values)\n",
    "\n",
    "ax = plt.gca()\n",
    "\n",
    "# and the derivative vectors\n",
    "plt.quiver(\n",
    "    \n",
    "    # starting points\n",
    "    arrow_points, arrow_points_y,\n",
    "    \n",
    "    # arrow directions\n",
    "    np.ones(len(arrow_points)) * arrow_x_size, arrow_points_derivative * arrow_x_size,\n",
    "    \n",
    "    angles = 'xy',\n",
    "    \n",
    "    zorder = 2,\n",
    "    \n",
    "    # smaller values imply longer arrows !\n",
    "    scale = 1500,\n",
    "    scale_units = None,\n",
    "    # units = 'xy',\n",
    "    # width = 0.01,\n",
    ")\n",
    "\n",
    "ymin, ymax = plt.ylim()\n",
    "\n",
    "plt.ylim(ymin - 0.2 * (ymax - ymin), ymax)\n",
    "\n",
    "plt.xlabel('$\\sigma$')\n",
    "plt.ylabel('NLL')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can thus also use a gradient based optimization algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find the minimum NLL using an optimizer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "set the initial guess for $\\sigma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "session.run(sigma.assign(0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the optimizer object. This uses `scipy.optimize.minimize` to drive the optimization but uses tensorflow to evaluate \n",
    "\n",
    "Since we specify bounds on the variable but we do not have constraints, `scipy.optimize.minimize` should use the `L-BFGS-B` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/tensorflow/contrib/learn/python/learn/datasets/base.py:198: retry (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use the retry module or similar alternatives.\n"
     ]
    }
   ],
   "source": [
    "optimizer = tf.contrib.opt.ScipyOptimizerInterface(nll, \n",
    "                    # the variables to minimize\n",
    "                    var_list = [ sigma ], \n",
    "                                                   \n",
    "                    # avoid going to negative values of sigma\n",
    "                    var_to_bounds = { sigma: (0.01, 1)},\n",
    "                                                   \n",
    "                    method = 'L-BFGS-B',\n",
    "                    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a callback function to monitor the progress of the optimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_callback(nll, sigma):\n",
    "    print \"neg. log likelihood is\",nll,\"for sigma=\",sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the optimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "neg. log likelihood is -60997.5 for sigma= 0.2\n",
      "neg. log likelihood is -56434.52 for sigma= 0.01\n",
      "neg. log likelihood is -61202.332 for sigma= 0.13707024\n",
      "neg. log likelihood is -58122.992 for sigma= 0.016932994\n",
      "neg. log likelihood is -61257.17 for sigma= 0.10237107\n",
      "neg. log likelihood is -61257.414 for sigma= 0.096986935\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.703 for sigma= 0.09923459\n",
      "neg. log likelihood is -61257.703 for sigma= 0.09934764\n",
      "neg. log likelihood is -61257.707 for sigma= 0.099354655\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "neg. log likelihood is -61257.707 for sigma= 0.09935468\n",
      "INFO:tensorflow:Optimization terminated with:\n",
      "  Message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n",
      "  Objective function value: -61257.707031\n",
      "  Number of iterations: 5\n",
      "  Number of functions evaluations: 24\n"
     ]
    }
   ],
   "source": [
    "optimizer.minimize(session,\n",
    "                   feed_dict = { x : xpoints, \n",
    "                                 bin_contents_ph : bin_contents\n",
    "                               },\n",
    "                   # arguments to call the loss callback function with\n",
    "                   fetches = [ nll, sigma ],\n",
    "                   \n",
    "                   loss_callback = loss_callback,\n",
    "                  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the minimization successfully completes, the $\\sigma$ variable object contains the optimized value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maximum likelihood estimate for sigma is 0.09935468\n"
     ]
    }
   ],
   "source": [
    "print \"maximum likelihood estimate for sigma is\",sigma.eval()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

uniform-analytic-convolution.png