Notebooks/MarkovChainMonteCarloSampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PYTHON VERSION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6.5 |Anaconda, Inc.| (default, Mar 29 2018, 13:32:41) [MSC v.1900 64 bit (AMD64)]\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print(sys.version)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rejection sampling in large dimensions.\n",
    "\n",
    "Already knowing how to uniformly select a sample from a surface $D$ in the plane, the area of any disk can easily be calculated. This method can be extended in an obvious way to more dimensions. The idea is that if a point is randomly taken inside the cube with center 0 and edge length 2, then the probability of this point being inside the unit ball is \n",
    "\n",
    "$$p(d)=\\frac{\\omega(d)}{2^d}.$$\n",
    "\n",
    "So if a large number $N$ of points is taken inside the cube, by the law of large numbers, the percentage of them that will \"fall\" into the ball is a good estimate of $p(d)$. This is exactly what the *vol* function defined in the routine below does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Monte Carlo estimate of ω( 3 ) is  : 4.20320 \n"
     ]
    }
   ],
   "source": [
    "from numpy import random\n",
    "\n",
    "def vol(N,d):\n",
    "    nhits = 0\n",
    "    for i in range(N):\n",
    "        x = random.uniform(-1,1,d)\n",
    "        if sum(x ** 2) < 1: \n",
    "            nhits += 1\n",
    "    return 2 ** d * nhits / N\n",
    "\n",
    "d = 3\n",
    "N = 10_000       \n",
    "print (\"The Monte Carlo estimate of ω(\", d,\") is  : %.5f \" % vol (N,d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function returns the volume of the unit ball in $d$ dimensions using the analytic formula: \n",
    "$$\\omega(d) =\\frac{\\pi^{d/2}}{\\Gamma\\big(\\frac{d}{2}+1\\big)}$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The actual value of ω( 3 ) is  : 4.18879 \n"
     ]
    }
   ],
   "source": [
    "from math import gamma, pi \n",
    "\n",
    "def Vol1(d):\n",
    "    x = d/2\n",
    "    return pi ** x / gamma(x + 1)\n",
    "\n",
    "print(\"The actual value of ω(\", d,\") is  : %.5f \" % Vol1(d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sampling with MCMC\n",
    "The rejection method is not practical for sampling from a large-dimensional sphere. This problem can be solved using ideas from Markov Chain Monte Carlo (MCMC).\n",
    "\n",
    "In particular, for a irreducible Markov chain with symmetric transition probabilities, the only invariant distribution is the uniform distribution in state space. Indeed, because the transition probabilities are symmetric the transition probability matrix will be doubly stochastic. Therefore, in order to randomly select a point from a space $D$, a Markovian chain that evolves in $D$ can be considered, with symmetric transition probabilities and enough steps for it to come close to its equilibrium state. \n",
    "\n",
    "One way to achieve this is as follows. Let $X_n= x\\in D$. Displacement by $z$ where the step distribution $z$ is symmetric around 0 can have two outcomes. The fisrt outcome is making the move, if the point $x+z$ remains inside $D$, i.e. defining $X_{n+1}=x+z$. The second outcome is remaining in $x$, i.e. defining $X_{n+1}=x$. \n",
    "\n",
    "The following algorithm implements this idea to take 1,000 samples of a random variable with a uniform distribution on the unit disk\n",
    "\n",
    "$$\n",
    "D_2=\\{(x,y)\\in\\mathbb{R}^2: x^2+y^2 < 1\\}.\n",
    "$$\n",
    "\n",
    "The chain starts from zero and at each step tries a parallel shift to one of the axes, \n",
    "with a uniform distribution in `(-delta, +delta)`. \n",
    "If the displacement would drive the chain out of the disk, it is not realized and the chain remains in the position it was in at the previous time. The chain takes $N=100$ steps before sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random as r\n",
    "\n",
    "delta = 1.0\n",
    "N = 100 # number of steps to approach equilibrium\n",
    "samples = 1000 \n",
    "point_x = []\n",
    "point_y = []\n",
    "\n",
    "for _ in range(samples):\n",
    "    x = [0,0]  ## start at the centre of th disc. This variable will keep the position of the chain\n",
    "    R_sq = 0.0  ## this variable keeps the squared distance from 0. It saves some computations to keep it\n",
    "    \n",
    "    for _ in range(N):\n",
    "        k = r.choice([0,1])  ## choose a jump direction at random\n",
    "        z = r.uniform(-delta,delta) ## choose a jump size uniformly in (-delta,delta)\n",
    "        x_prop_k = x[k] + z   ## propose a jump by z in the direction k\n",
    "        R_sqprop = R_sq - x[k]**2+ x_prop_k**2 ## compute the squared distance from 0 after the proposed jump \n",
    "        if R_sqprop < 1.0: \n",
    "            R_sq = R_sqprop\n",
    "            x[k]= x_prop_k   ## if the proposed jump leads to a point in the disc, then jump\n",
    "   \n",
    "    point_x.append(x[0])\n",
    "    point_y.append(x[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next up taking a sample from the cylinder:\n",
    "\n",
    "$$\n",
    "C_3=\\{\\,(x_1,x_2,x_3)\\in\\mathbb{R}^{3}:(x_1,x_2)\\in D_2,\\, |x_3| <1\\}.\n",
    "$$ \n",
    "\n",
    "This could be done by selecting $(x_1,x_2)$ uniformly on the disk $D_2$ and then the coordinate $x_3$ uniformly on (-1,+1). The probability that a randomly selected point from cylinder $C_3$ falls inside the ball $D_3=\\{(x_1,x_2,x_3)\\in\\mathbb{R}^3: x_1^2+x_2^2+x_3^2<1\\}$ is then \n",
    "\n",
    "$$ \\frac{|D_3|}{|C_3|}= \\frac{|D_3|}{2|D_2|}$$\n",
    "\n",
    "since the volume of the cylinder $|C_2|$ is the base surface area ($|D_2||$) times its height (2). Using samples from cylinder $C_3$ and by counting how many of them fell into $D_2$ an estimation for the ratio $|D_3|/|D_2|$ can be done.\n",
    "\n",
    "This procedure can be generalized to $d$ dimensions without problem. It is not difficult to modify the above algorithm to take random points in the $d$-dimensional unit ball\n",
    "\n",
    "$$ D_d = \\{(x_1,\\ldots,x_{d})\\in\\mathbb{R}^{d}:\\ x_1^2+\\cdots+x_d^2 <1\\}.$$ \n",
    "\n",
    "The chain by construction remains inside the $D_d$ ball, it is imprtant to let the chain take enough steps to get close to equilibrium. \n",
    "Then just by taking random points from the cylinder \n",
    "\n",
    "$$\n",
    "C_{d+1}=\\{\\,(x_1,\\ldots,x_{d+1})\\in\\mathbb{R}^{d+1}:(x_1,\\ldots,x_d)\\in D_d,\\, |x_{d+1} | <1\\}\n",
    "$$\n",
    "\n",
    "and by choosing their first $d$ coordinates uniformly in $D_d$ and the last one uniformly in (-1,1) the volume can be estimated.\n",
    "The percentage of these points falling inside the ball $D_{d+1}$ can give us an estimate of the ratio $|D_{d+1}|/|D_d|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application 1\n",
    "\n",
    "i) Write a routine that prints for dimensions $d=2,3,...$ the Monte Carlo estimated volume of the unit ball in $d$ dimensions, its actual value and the relative error of the estimate. Your routine should select $N=10^6$ points in the corresponding cube of each dimension, and should reach up to that dimension value where none of the $N$ selected points fell on the ball. \n",
    "\n",
    "ii) Plot the graph $d- \\log p(d)$ for the theoretical value of the probability $p(d)$ (notice that the relative error of the estimate in (i) gets larger as the dimensions increase).\n",
    "\n",
    "*In large dimensions, the fraction of the unit cube occupied by the ball inscribed in it is very small*. That is, in large dimensions, most of the points of the unit cube are located at its \"corners\". This is perhaps counter intuitive given the results from 2 and 3 dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Monte Carlo estimate of ω( 2 ) is  : 3.14336 \n",
      "The actual value of ω( 2 ) is  : 3.14159 \n",
      "The the realtive error for d =  2  is  : 0.03499 \n",
      "The Monte Carlo estimate of ω( 3 ) is  : 4.18478 \n",
      "The actual value of ω( 3 ) is  : 4.18879 \n",
      "The the realtive error for d =  3  is  : 0.03442 \n",
      "The Monte Carlo estimate of ω( 4 ) is  : 4.91784 \n",
      "The actual value of ω( 4 ) is  : 4.93480 \n",
      "The the realtive error for d =  4  is  : 0.08759 \n",
      "The Monte Carlo estimate of ω( 5 ) is  : 5.25264 \n",
      "The actual value of ω( 5 ) is  : 5.26379 \n",
      "The the realtive error for d =  5  is  : 0.08475 \n",
      "The Monte Carlo estimate of ω( 6 ) is  : 5.17613 \n",
      "The actual value of ω( 6 ) is  : 5.16771 \n",
      "The the realtive error for d =  6  is  : 0.27802 \n",
      "The Monte Carlo estimate of ω( 7 ) is  : 4.68672 \n",
      "The actual value of ω( 7 ) is  : 4.72477 \n",
      "The the realtive error for d =  7  is  : 0.38948 \n",
      "The Monte Carlo estimate of ω( 8 ) is  : 4.05350 \n",
      "The actual value of ω( 8 ) is  : 4.05871 \n",
      "The the realtive error for d =  8  is  : 0.52569 \n",
      "The Monte Carlo estimate of ω( 9 ) is  : 3.29779 \n",
      "The actual value of ω( 9 ) is  : 3.29851 \n",
      "The the realtive error for d =  9  is  : 0.84441 \n",
      "The Monte Carlo estimate of ω( 10 ) is  : 2.57331 \n",
      "The actual value of ω( 10 ) is  : 2.55016 \n",
      "The the realtive error for d =  10  is  : 2.15249 \n",
      "The Monte Carlo estimate of ω( 11 ) is  : 1.88211 \n",
      "The actual value of ω( 11 ) is  : 1.88410 \n",
      "The the realtive error for d =  11  is  : 3.59004 \n",
      "The Monte Carlo estimate of ω( 12 ) is  : 1.27795 \n",
      "The actual value of ω( 12 ) is  : 1.33526 \n",
      "The the realtive error for d =  12  is  : 10.12544 \n",
      "The Monte Carlo estimate of ω( 13 ) is  : 0.85197 \n",
      "The actual value of ω( 13 ) is  : 0.91063 \n",
      "The the realtive error for d =  13  is  : 1.65460 \n",
      "The Monte Carlo estimate of ω( 14 ) is  : 0.60621 \n",
      "The actual value of ω( 14 ) is  : 0.59926 \n",
      "The the realtive error for d =  14  is  : 6.62670 \n",
      "The Monte Carlo estimate of ω( 15 ) is  : 0.32768 \n",
      "The actual value of ω( 15 ) is  : 0.38144 \n",
      "The the realtive error for d =  15  is  : 3.08636 \n",
      "The Monte Carlo estimate of ω( 16 ) is  : 0.39322 \n",
      "The actual value of ω( 16 ) is  : 0.23533 \n",
      "The the realtive error for d =  16  is  : 16.45456 \n",
      "The Monte Carlo estimate of ω( 17 ) is  : 0.13107 \n",
      "The actual value of ω( 17 ) is  : 0.14098 \n",
      "The the realtive error for d =  17  is  : 100.00000 \n"
     ]
    }
   ],
   "source": [
    "# Application 1\n",
    "# (i)\n",
    "# Subroutine includes vol function, for the approximation of the volume of the unit sphere.\n",
    "\n",
    "d = 2\n",
    "N = 1000000  \n",
    "while(vol(N,d)>0): # if vol(N,d)>0 then nhits=0\n",
    "    print (\"The Monte Carlo estimate of ω(\", d,\") is  : %.5f \" % vol (N,d))\n",
    "    print(\"The actual value of ω(\", d,\") is  : %.5f \" % Vol1(d))\n",
    "    print(\"The the realtive error for d = \", d,\" is  : %.5f \" % (100*abs(1-vol(N,d)/Vol1(d))))\n",
    "    d += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOX1+PHPyWQPCVsgCWsCBElYBA0oAhIUFbWCawXrvlV/oq22blXRarXW5WtbiwsudSviVpUqiopGcIMAyhJIWAJI2HcIGELI+f0xEx1jApNkZu7M5Lxfr3ll5s6de88lTzhz7/Pc84iqYowxxkQ5HYAxxpjQYAnBGGMMYAnBGGOMhyUEY4wxgCUEY4wxHpYQjDHGAJYQQpKIvCAif2nkZzNFREUk2t9xGRMMTWn/pmksIRhjjAEsIRhjjPGwhBACRGSAiMwXkT0i8hoQ78dtdxCRqSKyXURWiMhVXu8liMiLIrJDRJaKyC0iUub1/moRuV1ElnjW+beI+C02Y8C/7d/rkunVIrJeRDaIyB/8GG5Es4TgMBGJBd4BXgbaAG8A5/hxF68CZUAH4FzgARE50fPe3UAm0A04Cbiwjs//BjgF6A70BO70Y2ymmQtg+x8BZAMnA7eJyEg/bDPiWUJw3rFADPB3VT2gqm8Chf7YsIh0BoYCt6pqhap+BzwLXORZ5dfAA6q6Q1XLgH/WsZl/qepaVd0O3A+M80dsxngEqv3/WVX3quoi4N9Yu/WJJQTndQDW6c+rDK6pb2URKRKRcs9jmA/b3q6qe2ptu6PX+2u93vN+XteyNZ7PGOMvgWr/1m4bwYYmOm8D0FFExOuPoguwsq6VVbV3A7a9HmgjIsleSaELsM5r352AJZ7XnevYhveyLp5tGuMvgWr/nYFir+1Zu/WBnSE472ugCrhBRKJF5GxgkD82rKprga+Av4pIvIj0A64A/uNZ5XXgdhFpLSIdgfF1bOY6EekkIm2APwGv+SM2YzwC1f7vEpFEEekNXIa1W59YQnCYqlYCZwOXAjuA84H/+nEX43B3HK8H3gbuVtWPPe/di7vDeRXwCfAmsL/W5ycDHwGlnofdMGT8JoDt/3NgBTADeERVP/LDNiOe2AQ5poaIXAuMVdXhntergStV9RNHAzPGRyKSifsLToyqVjkbTfixM4RmTEQyRGSIiESJyBHAH3CfRRhjmiHrVG7eYoGngSxgJzAFeMLRiIwxjrFLRsYYYwC7ZGSMMcYjrC4ZpaamamZmptNh/MzevXtJSkpyOowmsWNwmzdv3lZVbeenkBok1Np2JLQJiIzjCGbbDquEkJmZydy5c50O42cKCgrIz893OowmsWNwE5F675ANtFBr25HQJiAyjiOYbdsuGRljjAEcTggiMkpESjxlmW9zMhZjfHG4NisiXUTkMxH5VkQWishpTsRpTGM4lhBExAVMBE4FcoFxIpLrVDzGHI6PbfZO4HVVHQCMxYbxmjDi5BnCIGCFqpZ6bl+fAoxxMB5jDseXNqtAiud5S6yomgkjTnYqd+TnJWrLgGNqryQiVwNXA6SlpVFQUBCU4HxVXl4ecjE1lB2Dz3xps/cAH4nI9UASYBOzmLDhZEKQOpb94i45VZ0ETALIy8vTUBsxYKMYQkOQjsGXNjsOeEFVHxWRwcDLItJHVat/sbEQ/rITCV8SIDKOI5jH4GRCKOPntfY7YafXJrT50mavAEYBqOrXnjmoU4HNtTcWyl92IuFLAkTGcQTzGJzsQygEskUkyzOv6lhgamM29MzMUr4p3ebX4Iypgy9t9nvgRAARycE9YfyWxuxs1da9PDy9mE27K5oQsjG+cywheErTjgemA0txj8woauh2Kg4cZPKc7xk76Rvuemcxe/dbxVsTGPW1WRG5V0RGe1b7A3CViCwAXgUu1UYWDNuyZz8TP1vJ0g27/RG+MYfl6J3KqjoNmNaUbcTHuHj/hqE8PL2EF75azafFm/nbOf0Ymp3qpyiN+UldbVZVJ3g9XwIM8ce+urVzlyso3bKX/CP8sUVjDi0i7lROjI3m7jN688ZvBxMXHcWFz83m9v8uZHfFAadDM6bR2ibFkhwfzaqte50OxTQTEZEQauRltmHa74Zx9fHdeK1wLac8NpOCkl/05RkTFkSEbu1aWEIwQRNRCQHcl5D+dFoOb117HC3iorn034X88Y0F7NpnZwsm/HRLTaJ0S7nTYZhmIuISQo0BXVrz3g1DuW5Ed97+dh0nPfY5nyzZ5HRYxjRIVmoS63dV8EPlQadDMc1AxCYEgLhoFzef0ot3/t8Q2iTFcuVLc5nw7mKqDv7iHiFjQlJNx7JdNjLBENEJoUbfTi2ZOn4oVw3L4qWv13DVS3Mpt+GpJgxkpVpCMMHTLBICQGx0FHecnsv9Z/Vh5vKtnPfU12zY9YPTYRlzSD8lBOtHMIHXbBJCjd8c05XnLx3I2u37OHPilxSt3+V0SMbUKzE2moyW8ZRusTMEE3jNLiEADO/ZjjeuGUyUCOc99TWfFdvQVBO6slKTKLVLRiYImmVCAMjJSOGd64aQlZrEFS8W8vLXq50OyZg6dWvnHnrayAoYxvis2SYEgLSUeF7/7WBO6NWeu94t4r73lnCw2v7oTGjJSm3B7ooqtu+tdDoUE+GadUIASIqL5umL8rj0uEye+2IV17wyj32VNgLJhA4bemqCpdknBABXlHDP6N7cfUYuM5ZuYuykb9iyZ7/TYRkDuO9WBqxj2QScJQQvlw3JYtJFeSzfVM5vnv2GbeWWFIzzOrZKIMYl1rFsAs4SQi0jc9N47pI81mzbx4XPzWHnPrtua5wV7Yqia9skuxfBBJwlhDoc1yOVZy7OY+WWci58bja7frDCeMZZWalJdsnIBJwlhHoc37MdT194NCUb93Dx83PYY3MrGAd1a5fEmm37bBScCShLCIcwold7nvjN0RSt28Wl/y60+kcGERklIiUiskJEbqvj/cdE5DvPY5mI7PTHfrulJlF5sJp1O6zcigkcSwiHcVJuGo+PG8B3a3dy+QuFNiS1GRMRFzAROBXIBcaJSK73Oqp6o6r2V9X+wOPAf/2x76zUFgCUWj+CCSBLCD44tW8Gj53fn7mrt3Pli3OtNn3zNQhYoaqlqloJTAHGHGL9ccCr/tix3YtggiHa6QDCxegjO3CwupqbXl/A1S/P5ZmL84iPcTkdlgmujsBar9dlwDF1rSgiXYEs4NP6NiYiVwNXA6SlpVFQUFDvjlWVhGiYtWAZWQfWNDzyBiovLz9kPOEiEo4jmMdgCaEBzhrQiQMHlVveXMi1r8zjqYuOdjokE1xSx7L6ennHAm+qar2nk6o6CZgEkJeXp/n5+Yfcec+iL6iMjSE/v84c5FcFBQUcLp5wEAnHEcxjsEtGDfTrvM48cFZfPivZwvjJ31Jloz6akzKgs9frTsD6etYdi58uF9Xo1q6FXTIyAWUJoREuOKYLfx7dm4+XbOLFokqrQtl8FALZIpIlIrG4/9OfWnslETkCaA187c+dZ6UmsW7nD9aHZQLGEkIjXXJcJjec0INZ66p4omCl0+GYIFDVKmA8MB1YCryuqkUicq+IjPZadRwwRf38TaGmY3n1NjtLMIFhfQhNcONJPSlcuoqHp5fQqXUCY/p3dDokE2CqOg2YVmvZhFqv7wnEvrO8itzlZKQEYhemmbOE0AQiwuV94zgYH8/NbyykQ6sEBma2cTosE6FsfmUTaHbJqIliooRJFx1Np9YJXPXSXOv0MwHz4/zK1sZMgDiSEETkPBEpEpFqEclzIgZ/apUYy78vG0iUCJf9e47NbGUCxorcmUBy6gxhMXA2MNOh/ftd17ZJPHPx0azfVcHVL82l4oCNBDH+Z/Mrm0ByJCGo6lJVLXFi34F0dNc2PPbr/sxds4Ob31xItd2jYPysZn7lHfus+q7xv5DvVG7I7f1OqH1beRJwXs8Y3liwHt2zhXN7xjoWm6/s9v7w8dN0muW0SbIBDMa/ApYQROQTIL2Ot+5Q1Xd93U5Db+8PtrpuKx8+XIl+exGvzlnL8QNy+PXAznV/OETY7f3ho+ZehNKte8mzEW3GzwKWEFR1ZKC2HepEhHvH9GHdzgr+9PYiOrRKYGh2qtNhmQjw4/zK1rFsAsCGnQZIjCuKiRcMoEf7Flz7yjxWbLax46bpbH5lE0hODTs9S0TKgMHA+yIy3Yk4Ai05PobnLx1IbHQU17wyz2ZcM36RlZpk97uYgHBqlNHbqtpJVeNUNU1VT3EijmDo0CqBxy8YQOmWcm55c4ENFzRN1i01idU2v7IJALtkFATHdU/ltlN7MW3RRp6ZVep0OCbMdWuXRGVVNet32vzKxr8sIQTJVcO6cXrfDB78oJivVm51OhwTxmrmV165xfoRjH9ZQggSEeFv5/ajW7sWXD/5W/t2ZxrN5lc2gWIJIYhaxEXz9EVHs7+qmmv/M5/9VVbewjRc26RYkuOjLSEYv7OEEGTd27XgkfOOZMHanfz5f0ucDseEIRGhmxW5MwFgCcEBo/qkc21+dybP/p7XC9c6HY5pABEZJSIlIrJCRG6rZ51fi8gST0XfyYGIw+ZXNoFgCcEhfzz5CIb2SOXOdxezsGyn0+EYH4iIC5gInArkAuNEJLfWOtnA7cAQVe0N/D4QsdTMr2xVdY0/WUJwiCtK+Oe4AbRrEce1r8y3ORTCwyBghaqWqmolMAUYU2udq4CJqroDQFU3ByIQ61g2gRDy1U4jWZukWJ688CjOfeprbnj1W168fBCuKHE6LFO/joD3Nb4y4Jha6/QEEJEvARdwj6p+WNfGmlLJd9tu95nBe5/PYVO6//+MI6V6bCQcRzCPwRKCw/p1asV9Y3pz61uLePSjEm4Z1cvpkEz96srWtW8XjgaygXygEzBLRPqo6i+uCzalku++yiru/mo6Ce27kp+f7fPnfBUp1WMj4TiCeQx2ySgEnD+wC+MGdeaJgpUUlATkCoPxjzLAu5Z5J2B9Heu8q6oHVHUVUII7QfiVza9sAsESQoi4+4ze9EpP5o9vLGDzngqnwzF1KwSyRSRLRGKBscDUWuu8A4wAEJFU3JeQAlKvxIrcGX+zhBAi4mNcPD5uAOX7q/jD6wts+s0QpKpVwHhgOrAUeF1Vi0TkXhEZ7VltOrBNRJYAnwE3q+q2QMTjnl95rxVMNH5jCSGEZKclc/cZvZm1fKsVwQtRqjpNVXuqandVvd+zbIKqTvU8V1W9SVVzVbWvqk4JVCxZqS3Y9cMBm1/Z+I0lhBAzdmBnTuubzsPTS1iw1u5PMPWrmV/ZJssx/mIJIcSICH89qx9pKfHcMOVb9lTYtz9Tt5p7EVZaCQvjJ5YQQlDLxBj+MbY/a7fv4653Fts1YlOnmvmVrWPZ+IslhBCVl9mG34/syTvfree/89c5HY4JQTXzK5favAjGTywhhLDrRvTgmKw23PXuYvujN3WyoafGnywhhDBXlPD3sf2JjY7ihinfUllV7XRIJsTY/MrGnywhhLiMlgk8dE4/Fq/bzUMfFjsdjgkxNr+y8SdLCGHg5N7pXDy4K89+sYrPrLSF8VIzv7KVsDD+YAkhTPzptBx3aYvXrbSF+UnN0FPrYzL+YAkhTNSUtthbaaUtzE9sfmXjT5YQwkh2WjJ3/SqXWcu38srsNU6HY0JAzfzKlhCMP1hCCDMXDOrC8J7teGDaUvtPwADu+ZVL7W5l4weWEMKMiPC3c/oR64riD69/Z8MNjc2vbPzGEkIYSm8Zz31n9mH+9zuZNNOqojZ3Pdq7Rxqt2Gwdy6ZpHEkIIvKwiBSLyEIReVtEWjkRRzgbfWQHTu2TzmMfL6N4426nwzEO6pWeDMDSDdYOTNM4dYbwMdBHVfsBy4DbHYojbIkIfzmzDykJ0dz02gK7i7kZ69o2iYQYF0s37HE6FBPmHEkIqvqRZ/YpgG9wz01rGqhtizj+enY/lmzYzeOfLnc6nGZBREaJSImIrBCR2+p4/1IR2SIi33keVwY6JleUcER6sp0hmCaLdjoA4HLgtfreFJGrgasB0tLSKCgoCFJYvikvL3c0phhgaMdoJn62gtb7yujWytXgbTh9DP4QjGMQERcwETgJKAMKRWSqqi6pteprqjo+oMHUkpORzAeLN6KqiEgwd20iSMASgoh8AqTX8dYdqvquZ507gCrgP/VtR1UnAZMA8vLyND8/3//BNkFBQQFOx3TUsQcY9dhMXlnpYtoNw4iPaVhSCIVjaKogHcMgYIWqlgKIyBRgDFA7IQRdTkYKr85Zy6bd+0lvGe90OCZMBSwhqOrIQ70vIpcAvwJOVJsBpklS4mN4+Lwj+c2zs3nowxImnJHrdEiRqiOw1ut1GXBMHeudIyLH4+4fu1FV19axjl/Pfit2uIecvvbRFxzZrul/1pFw1giRcRzBPAZHLhmJyCjgVmC4qu5zIoZIM6RHKpcM7srzX67ipNw0Bndv63RIkaiuazG1v8z8D3hVVfeLyDXAi8AJdW3Mn2e/R1Uc4IHZHxGdmkl+fo9Gb6dGJJw1QmQcRzCPwalRRv8CkoGPPR1vTzkUR0S57dQcslKT+OMbCyjfX3X4D5iGKgM6e73uBKz3XkFVt6nqfs/LZ4CjgxFYSnwMnVonULzRRhqZxnNqlFEPVe2sqv09j2uciCPSJMS6eOS8I9mw6wf+8p7jl7UjUSGQLSJZIhILjAWmeq8gIhleL0cDS4MVXE5Gio00Mk1idypHmKO7tua3w7szpXAtnxZvcjqciOIZKj0emI77P/rXVbVIRO4VkdGe1W4QkSIRWQDcAFwarPhy0pMp3VJuJSxMo1lCiEC/H5lNr/RkbntrEbt+OOB0OBFFVaepak9V7a6q93uWTVDVqZ7nt6tqb1U9UlVHqGrQprnLyUihWmHZJrtsZBrHEkIEiot28fC5R7K1fD8PfhC0KxbGYTkZKQAU2x3LppEsIUSovp1actXx3Xh1zlq+WrnV6XBMEHRpk0hirIsl1o9gGskSQgS7cWRPMtsmcvt/F/FDpV1XjnRRVsLCNJElhAgWH+Pir2f3Y822fTz2yTKnwzFBkJORQvHGPdi9nqYxLCFEuMHd2zJuUBeenVXKwrKdTodjAiwnI4VdPxxgw64Kp0MxYcgSQjNw+2m9aJccxy1vLrQy2REux+ZGME1gCaEZSImP4S9n9qV44x6e/nyl0+GYAOrlGWlkCcE0hiWEZuKk3DR+1S+Dxz9dwYrNNiwxUrWIi6ZLm0SWWgkL0wiWEJqRe0b3JjHOxa1vLaK62jodI1UvG2lkGskSQjOS2iKOCb/KZd6aHbz8zRqnwzEBkpORwuqte22osWkwSwjNzFkDOnJ8z3b87cNiynZY5fFIZCUsTGNZQmhmRIQHzuoDwJ/eXmzj1SNQrnUsm0ayhNAMdWqdyK2jejFz2Rbe/nad0+EYP+vUOoGkWJclBNNglhCaqYuO7crRXVtz73tL2LXfzhIiSVSU0CsjhaVW5M40kE8JQUS6icj/RGSriGwWkXdFpFuggzOBExUl/O2cvuzbf5DJS/cf/gMmrORkJLN04267JGgaxNczhMnA60A60AF4A3g1UEGZ4OjRPpnrRvRg9saDzFy2xelwjB/1Sk9hT0UV63b+4HQoJoz4mhBEVV9W1SrP4xV+Obm4CUPX5HcjPVG4693FNtOWD0RklIiUiMgKEbntEOudKyIqInnBjK9Gzo8dy3bZyPjO14TwmYjcJiKZItJVRG4B3heRNiLSJpABmsCKi3Zxce841mzbxxMFVtbiUETEBUwETgVygXEiklvHesm4p8+cHdwIf9LLU9Oo2DqWTQNE+7je+Z6fv621/HLcZwrWnxDGctu6OLN/B54qWMmY/h3o3q6F0yGFqkHAClUtBRCRKcAYYEmt9e4DHgL+GNzwfpIUF03Xtoks3WgJwfjOp4SgqlmBDsQ4647Tc5lRvJm73lnMf648BhFxOqRQ1BFY6/W6DDjGewURGQB0VtX3ROSQCUFErgauBkhLS6OgoMCvwaZG72feyk2N2m55ebnf43FCJBxHMI/Bp4QgImfXsXgXsEhVN/s3JOOEdslx3DqqF3e+s5h3v1vPmQM6Oh1SwFVUVPDEE0/wxRdf1CTA9iISr6r1TSZQV5b8sS9NRKKAx4BLfdm/qk4CJgHk5eVpfn5+A6I/vAVVy/n7jGUMOm4oibG+XgxwKygowN/xOCESjiOYx+BrH8IVwLPAbzyPZ4CbgC9F5KIAxWaC7IJBXejfuRV/eX8Ju/YdcDqcgLv44ospKiri+uuvZ/z48QAJwMuH+EgZ0NnrdSdgvdfrZKAPUCAiq4FjgalOdSz3ykhGFUqs8qnxka9fG6qBHFXdBCAiacCTuE+XZ3LoPyITJqKihPvP6sMZj3/BQ9OLuf+svk6HFFAlJSUsWLDAe9EaoOchPlIIZItIFrAOGAtcUPOmqu4CUmtei0gB8EdVnevHsH2W6zXSaECX1k6EYMKMr2cImTXJwGMz0FNVtwOR/1WyGendoSWXDcli8pzvmf/9DqfDCagBAwbwzTffeC9KAr6sb31VrQLGA9OBpcDrqlokIveKyOiABtsInVonkBwXbSUsjM98PUOYJSLv4b4hDeBcYKaIJAE2UW+EufGknry/cAN3vL2Y/40fQrQrMiuczJ49m5deeokuXbrULOoFVIvIIkBVtV/tz6jqNGBarWUT6tq+qub7N+KGERF6ZSRTbCONjI98TQjXAWcDQ3F3rL0IvKXu++JHBCg245AWcdHcMzqXa16ZzwtfrebKYZE5qvjDDz/82evMzMxFwBnORBMYvdJTeOfbdaiqjRwzh+XrsFMVkS+AStyjKuZoE4qkiMh9uMdvV+O+/HSpqq4/9KdMMJ3SO50TerXn/z5exml9M+jQKsHpkPyua9eutRdVqmpEzRyUk5HCy9+soWzHD3Ruk+h0OCbE+Vrc7tfAHNyXin4NzBaRc5uw34dVtZ+q9gfeA+o85TbOERH+PLo31arc+7/a912ZcJGT4b5jeYn1Ixgf+Hpx+A5goKpeoqoX475j867G7lRVvVtnElYXKSR1bpPIDSdm82HRRmYs3XT4D5iQc0R6MiJQbDWNjA987UOIqnUD2jaaOJeCiNwPXIz7Brd6+yECfTdnU0X6nZA9q5UOLYRbXpvHA8MSiHOF5nXoSPg9BEJibDSZbZNspJHxia8J4UMRmc5PJa/Pp9ZIi9pE5BPc5bJru0NV31XVO4A7ROR23EP57q5rO4G+m7OpmsOdkMmZ2zh/0jcsOJDBLSf2Cl5gDRAJv4dAyclIpmi9JQRzeL52Kt8sIucAQ3CPMpqkqm8f5jMjfYxhMvA+9SQE47xjurXl7KM68uysVZyX15ms1CSnQzIN0Cs9hWmLNrJ3fxVJcQ0rYWGaF58v+6jqW6p6k6reeLhkcDgiku31cjRQ3JTtmcC77dRexEZH8ef/FdksXGGmZm6EYithYQ7jkAlBRPaIyO46HntEpCnnoA+KyGIRWQicDPyuCdsyQdA+OZ7fj8ymoGQLM5ZaPcNwUjPSyPoRzOEc8vxRVZMDsVNVPScQ2zWBdclxmbxWuJZ731vC0OxU4mNcTodkfNCxVQLJ8VbCwhxeZNYkMAER44rintG9+X77Pp6ZWep0OMZHIkJOeopdMjKHZQnBNMiQHqmc3jeDiQUrKNuxz+lwjI9yMpIp3rCb6mrr/zH1s4RgGuxPp+cAcP/7Sx2OxPgqJyOFvZUHWWtJ3ByCJQTTYB1bJTB+RA8+WLyRL5ZvdToc44NeP86NYP0Ipn6WEEyjXDmsG13bJnL31MVUVlU7HY45jCPSkokS92Q5xtTHEoJplPgYFxN+lcvKLXt58avVTodjDiMh1kVmqpWwMIdmCcE02ok5aZzQqz3/mLGczbvrm5c+sojIKBEpEZEVInJbHe9fIyKLROQ7EflCRHKdiLMuOekpLLXJcswhWEIwTTLhV7lUVlXz4AeRf7O5iLiAicCpQC4wro7/8Ceral9PafeHgP8Lcpj16tOxJWu3/8CWPfudDsWEKEsIpkkyU5O46vgs/vvtOgpXb3c6nEAbBKxQ1VJVrQSm4J7o6UehXNp9SI+2AHy5wgYCmLpZpSvTZNeN6MF/56/j7neL+N/1Q3FFhWaJbD/oCKz1el0GHFN7JRG5DrgJiAVOqG9jwS7tXq1Kixh4Y9ZiWu1afsh1I6WceCQcRzCPwRKCabLE2GjuPD2X6ybPZ/Kc77no2F9MTRkp6sp0vzgDUNWJwEQRuQC4E7ikro05Udo9f8N85qzazvDhww85x3KklBOPhOMI5jHYJSPjF6f1TWdwt7Y8Mr2E7XsrnQ4nUMqAzl6vOwGHmgt8CnBmQCNqoOOz27F5z36WbSp3OhQTgiwhGL8QEf48pjfl+6t49KMSp8MJlEIgW0SyRCQWGAtM9V6hVmn304FDX5sJsqHZqQDMWr7F4UhMKLKEYPymZ1oyFx3blVfnfE9JBBZSU9Uq3LP7TQeWAq+rapGI3Csioz2rjReRIhH5Dnc/Qp2Xi5zSoVUC3dslMdPuMDd1sIRg/Op3J2aTHB/Dfe8ticiJdFR1mqr2VNXuqnq/Z9kEVZ3qef47Ve2tqv1VdYSqFjkb8S8Ny27H7NJtVBw46HQoJsRYQjB+1Toplt+PzOaLFVttIp0QdXzPVPZXVTNvzQ6nQzEhxhKC8bsLj+1K93ZJ3D9tqdU5CkHHZLUlxiXMtH4EU4slBON3Ma4o7vxVLqu27uWlr1c7HY6pJSkumqO7tmbWMutHMD9nCcEExIgj2jO8Zzv+MWM528qtVEKoGZbdjiUbdlsZC/MzlhBMwNx5eg77Kg/y2CfLnA7F1DLMM/zUylgYb5YQTMBke4ahTp4dmcNQw1nvDi1pnRhj/QjmZywhmICK9GGo4coVJQzpkcoXy7fa78X8yBKCCSgbhhq6rIyFqc0Sggk4G4YamqyMhanNEoIJOBuGGpo6tEqgR/sWVsbC/MgSggkKG4Yamob2SLUyFuZHlhBM0Ngw1NBjZSyMN0sIJmiy05K58JguNgw1hFgZC+PN0YQgIn8UERWRVCfjMMHz+5E9bRhqCLEyFsabYwlBRDoDJwHfOxWls0RQAAAShklEQVSDCT7vYaifFtsw1FBgZSxMDSfPEB4DbqGOOWlNZLvw2K50S03igWlLqTpow1Cddnx2O8DKWBiIdmKnntml1qnqgkNN9O1Z92rgaoC0tDQKCgoCH2ADlJeXh1xMDeXEMZzeuYrHv93Pff+ZwYguMU3eXrCOQURGAf8AXMCzqvpgrfdvAq4EqoAtwOWquibggTVB7w4pP5axOHNAR6fDMQ4KWEIQkU+A9DreugP4E3CyL9tR1UnAJIC8vDzNz8/3V4h+UVBQQKjF1FBOHMNwVb7Z8TXvf7+Xm88fSou4pjXFYByDiLiAibgvdZYBhSIyVVWXeK32LZCnqvtE5FrgIeD8gAbWRFG1ylgc7kuaiVwBu2SkqiNVtU/tB1AKZAELRGQ10AmYLyJ1JQ8ToUSEP52Ww9bySiZ9vtLpcHw1CFihqqWqWglMAcZ4r6Cqn6nqPs/Lb3C375BnZSwMONCHoKqLVLW9qmaqaibub1pHqerGYMdinDWgS2vOOLIDk2aVsnFXhdPh+KIjsNbrdZlnWX2uAD4IaER+YmUsDDjUh2BMjVtOOYLpizfyfx+X8NC5RzodzuHUdS2lzkERInIhkAcMr3djIdY/1iFJeGf2Mnoc/D4i+sbA+vgayvGE4DlLMM1U5zaJXHJcV579YhWXDckiJyPF6ZAOpQzo7PW6E7C+9koiMhJ3X9lwVa13LGeo9Y+dsqeIybO/59ghw/jmy1lh3zcG1sfXUHansnHc+BHZpMTH8MC0pU6HcjiFQLaIZIlILDAWmOq9gogMAJ4GRqtqWN1oMSzbXcZi7morY9FcWUIwjmuZGMP1J/Rg1vKtfL4sdK9hq2oVMB6YDiwFXlfVIhG51zOUGuBhoAXwhoh8JyJT69lcyKkpYzFrRej+DkxgOX7JyBiAiwZ35cWvV/PXaUsZ2iMVV1RoDn1U1WnAtFrLJng9Hxn0oPzEu4zF4JDvzjGBYGcIJiTERbu4dVQvijfu4a35ZU6H02zVlLHYtd8KCDRHlhBMyDi9bwb9O7fi0Y9K2FdZ5XQ4zVJNGYuibTY/QnNkCcGEDBHhztNz2LR7P8/NWuV0OM1STRmLoq2WEJojSwgmpORltmFU73Se/Hwlm/eExc1qESUqShiW3Y4FW6r4odKSQnNjCcGEnFtP7UVlVTV//2S506E0Sxce25XyA/DGvLWHX9lEFEsIJuRkpSZx4bFdea1wLcs32cxqwTYwszU9WkXx9OelHLDy5M2KJQQTkm44MZvEGBcPflDsdCjNjohwercY1u38gfcXbnA6HBNElhBMSGqTFMv/G9GDGcWb+XrlNqfDaXaObOeiZ1oLnixYaVOdNiOWEEzIumxIJhkt43nww2L7TynIokS4Znh3Sjbt4bOSsKrAYZrAEoIJWfExLm4c2ZMFa3cyvciqowfbGUd2oGOrBJ4sCJv5KkwTWUIwIe3sozrSo30LHppeYvMvB1mMK4qrhmVRuHoHhau3Ox2OCQJLCCakRbuiuPmUIyjdspc35llJi2A7f2AX2iTF8pSdJTQLlhBMyDs5N42jurTi758ss5ulgiwh1sWlx2Uyo3gzxRt3Ox2OCTBLCCbkiQi3neouafHvr6ykRbBdPLgribEunv681OlQTIBZQjBhYVBWG07s1Z4nC1ayc1+l0+E0K60SY7lgUBemLljP2u37nA7HBJAlBBM2bh51BOX7q3jCrmcH3ZXDuhEl8OwsO0uIZJYQTNjolZ7CWQM68sJXq1m/8wenw2lW0lvGc/aATkwpXMvW8nqniTZhzhKCCSs3ndQTFB77eJkj+xeRUSJSIiIrROS2Ot4/XkTmi0iViJzrRIyBcvXwblQerObFr1Y7HYoJEEsIJqx0ap3IRYO78tb8MpYFufCdiLiAicCpQC4wTkRya632PXApMDmowQVB93YtGNU7nRe/Wk35fpvAKBJZQjBh57oRPUiKjeahD0uCvetBwApVLVXVSmAKMMZ7BVVdraoLgYi8i+6a4d3ZXVHFq7O/dzoUEwDRTgdgTEO1SYrlt8O78chHy5i7ejt5mW2CteuOgPckAWXAMY3dmIhcDVwNkJaWRkFBQZOC86fy8vJ648ltG8XEGcVkVq0hJkqCG1gDHeo4wkUwj8ESgglLlw/N4sWv1/DgB8W8cc3gYO22rv/9Gl11T1UnAZMA8vLyND8/v7Gb8ruCggLqiye641YufG4225O7c/7ALsENrIEOdRzhIpjHYJeMTFhKjI3mdydmM3fNDmYsDVo1zjKgs9frTsD6YO08VAzp0Za+HVvy9OelHKy2KrSRxBKCCVvnD+xMVmoSD00vpjo45bELgWwRyRKRWGAsMDUYOw4lIsK1+d0p3bqXj6wKbUSxhGDCVowrij+efATLNpXz5brAj3pR1SpgPDAdWAq8rqpFInKviIwGEJGBIlIGnAc8LSJFAQ/MAaf0TicrNYm/flBs9yVEEEcSgojcIyLrROQ7z+M0J+Iw4e+0vun069SSt1ccoOJA4Avfqeo0Ve2pqt1V9X7PsgmqOtXzvFBVO6lqkqq2VdXeAQ/KAa4o4ZHzjmTT7gqueKGQvTYMNSI4eYbwmKr29zymORiHCWMiwq2jerG9QplsQyGD6uiurfnXBUexaN0urps8nwM2X0XYs0tGJuwN6ZFKTpsoJn62wr6pBtlJuWncf1ZfCkq2cPt/F9lUp2HOyWGn40XkYmAu8AdV3VHXSqE8VhtsnHOoOK3TQR5dWM1dr3zK6O6xTofTrIwb1IWNuyr4x4zlpKXEcfMpvZwOyTRSwBKCiHwCpNfx1h3Ak8B9uMdw3wc8Clxe13ZCeaw22DjnkFFQwMgDLfh41TYmjDuOVomWFILp9yOz2byngomfrSQ9JZ6LBmc6HZJphIBdMlLVkarap47Hu6q6SVUPqmo18AzukgDGNMkfTu5J+f4qnp5pJZqDTUS4b0wfRua0Z8LUIj5cvMHpkEwjODXKKMPr5VnAYifiMJElJyOFM/p14IUvV7N5T4XT4TQ70a4oHh93FAM6t+KGKd8xZ9V2p0MyDeRUp/JDIrJIRBYCI4AbHYrDRJgbT+pJ5cFqnvjMJtFxQkKsi+cuGUin1glc+WJh0CvSmqZxJCGo6kWq2ldV+6nqaFW180vjF1mpSfw6rxP/mb2Gsh023aMTWifF8uJlg4iLcXHJ83NsMqMwYsNOTcS5/oRsRIR/fLLc6VCarc5tEnnxskGUV1RxyfNz2LXvgNMhGR9YQjARp0OrBC461j2Jzsot5U6H02zldkjh6YuOZvW2vZz2z1m8Na/MiuGFOEsIJiJdm9+d+BgX/+fQVJvG7bgeqbxyxTG0TorhD28s4LR/zOLT4k12A1uIsoRgIlJqiziuGJrF+ws3sHjdLqfDadaO6daWqdcN5fFxA6ioOsjlL8zl/Ke/Yd6aOu9FNQ6yhGAi1pXDutEyIYZHPwr6VJumlqgo4YwjO/DxjcO5b0xvSrfu5Zwnv+Kql+ay3EYihQxLCCZitUyI4Zrh3fmsZAtzV9uY+FAQGx3FRYMz+fzmfG46qSdfr9zGKX+fyS1vLrDRSCHAEoKJaJcc15V2yXE8NL3ErluHkKS4aG44MZvPb87n0uOyeOfb9eQ/UsA1L8/j+S9WsXjdLuuAdoDNqWwiWmJsNNef0IMJ7xYxa/lWju/ZzumQjJe2LeKYcEYulw3J5MnPVzJz2RY+9MzC1iIumqO7tmZQVhsGZrahX6eWxMe4HI44sllCMBFv7MAuPP15KQ9PL2FYdioi0uhticgo4B+AC3hWVR+s9X4c8BJwNLANOF9VVzd6h81E5zaJPHBWXwDW7/yBwtXbmbNqO4Wrt/PwdHcfUGx0FP07tWJgVmsy2yaRlhJP+5Q40pLjaZUY06Tfq3GzhGAiXmx0FL8fmc3Nby5ketFGRvXJOPyH6iAiLmAicBJQBhSKyFRVXeK12hXADlXtISJjgb8B5zfxEJqVDq0SGNO/I2P6dwRgx95KCldv/zFJPPV56S8uJ8W6otzJISWe9snun+2S49jw/QE2F64lPtZFQoyLxFgX8THu5wmxntfRLlwuITpKcEW5fzbX5GIJwTQLZw3oyFOfr+TRj5ZxUm46rqhG/cEPAlaoaimAiEwBxgDeCWEMcI/n+ZvAv0RE1DowGq11Uiwn907n5N7uavoVBw6yefd+Nu2pYNPuCjbt3s/mPRXuZbsrWL65nC9WbGVPhWeypKULG7zPKIHoqChcniRR83A3G/dPEYgSQXBXexXPMsH93L3mT2qSzM9aXh3NsPaiH/btI3H+54eNefSRHbjhxGyfj7EulhBMsxDtiuLmU47g0+LN7KusIjk+pjGb6Qis9XpdBhxT3zqqWiUiu4C2wNbaGwvlyZ/CZdKkFp5H9wQgAUireSeOyoOxbN+9l5j4RPYfhMqD+uPPymp+en4QDipUq3p+el5Xu3+qKgdVqVZFAVX3RC4/e+5J99V4FvDjj5/x/lrgyzcEBaoSqomOOvwIrB3rV1NQsM6HrdbPEoJpNkb1yWj05SKPuk4rav9d+7KOe2EIT/4UEZMmERnHEcxjsGGnxviuDOjs9boTsL6+dUQkGmgJ2E0QJixYQjDGd4VAtohkiUgsMBaYWmudqcAlnufnAp9a/4EJF3bJyBgfefoExgPTcQ87fV5Vi0TkXmCuqk4FngNeFpEVuM8MxjoXsTENYwnBmAZQ1WnAtFrLJng9rwDOC3ZcxviDXTIyxhgDWEIwxhjjYQnBGGMMYAnBGGOMh4TTiDgR2QKscTqOWlKp4y7UMGPH4NZVVR0phxqCbTsS2gRExnEErW2HVUIIRSIyV1XznI6jKewYTG2R8u8ZCccRzGOwS0bGGGMASwjGGGM8LCE03SSnA/ADOwZTW6T8e0bCcQTtGKwPwRhjDGBnCMYYYzwsIRhjjAEsITSaiKwWkUUi8p2IzHU6Hl+JyPMisllEFnstayMiH4vIcs/P1k7GeDj1HMM9IrLO8/v4TkROczLGcBaObTsS2jU437YtITTNCFXtH2bjnF8ARtVadhswQ1WzgRme16HsBX55DACPeX4f/T1VSU3jhVvbfoHwb9fgcNu2hNDMqOpMfjmD1xjgRc/zF4EzgxpUA9VzDKYZi4R2Dc63bUsIjafARyIyzzNZejhLU9UNAJ6f7R2Op7HGi8hCz2l3yF8eCGGR0rYjpV1DkNq2JYTGG6KqRwGnAteJyPFOB9TMPQl0B/oDG4BHnQ0nrFnbDi1Ba9uWEBpJVdd7fm4G3gYGORtRk2wSkQwAz8/NDsfTYKq6SVUPqmo18Azh/ftwVAS17bBv1xDctm0JoRFEJElEkmueAycDiw/9qZDmPTH8JcC7DsbSKDV/+B5nEd6/D8dEWNsO+3YNwW3bdqdyI4hIN9zfnMA9L/VkVb3fwZB8JiKvAvm4S+puAu4G3gFeB7oA3wPnqWrIdtrWcwz5uE+pFVgN/Lbm+rHxXbi27Uho1+B827aEYIwxBrBLRsYYYzwsIRhjjAEsIRhjjPGwhGCMMQawhGCMMcYj2ukAjJuI3AOUAynATFX9JMj7Hw3kquqDwdyviWzWrsOLDTsNETV/OKr6iNOxGOMv1q7Di10ycpCI3CEiJSLyCXCEZ9kLInKu5/lqEXlARL4WkbkicpSITBeRlSJyjdd2bhaRQk/xqz97lmWKyFIReUZEikTkIxFJ8Lx3g4gs8aw/xbPsUhH5l+d5VxGZ4Xl/hoh08YrtnyLylYiU1sRpjDdr1+HLEoJDRORoYCwwADgbGFjPqmtVdTAwC3et9HOBY4F7Pds5GcjGXd+kP3C0VzGybGCiqvYGdgLneJbfBgxQ1X7Aj3+AXv4FvOR5/z/AP73eywCGAr8C7DTc/Iy16/BmCcE5w4C3VXWfqu7GXXelLjXLFwGzVXWPqm4BKkSkFe5aMycD3wLzgV64/2AAVqnqd57n84BMz/OFwH9E5EKgqo59DgYme56/jPsPpcY7qlqtqkuANJ+P1jQX1q7DmHUqO8uXDpz9np/VXs9rXkcDAvxVVZ/2/pCIZNZa/yCQ4Hl+OnA8MBq4S0R6NyBO723K4cM3zZC16zBlZwjOmQmcJSIJnuqSZzRyO9OBy0WkBYCIdBSReicCEZEooLOqfgbcArQCWtRa7Svcp/0AvwG+aGRspvmxdh3G7AzBIao6X0ReA74D1uC+ltqY7XwkIjnA1yIC7iF+F+L+5lQXF/CKiLTE/U3oMVXd6flsjRuA50XkZmALcFljYjPNj7Xr8GbDTo0xxgB2ycgYY4yHJQRjjDGAJQRjjDEelhCMMcYAlhCMMcZ4WEIwxhgDWEIwxhjj8f8Bu56pUrlb+/4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# (ii)\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "P = []\n",
    "D = []\n",
    "for d in range(2,17): # for d = 17 (dimension) nhits = 0 in the previous block of code\n",
    "    p = Vol1(d)/2**d\n",
    "    P.append(p)\n",
    "    D.append(d)\n",
    "\n",
    "logP = np.log10(P)\n",
    "\n",
    "plt.subplot(1, 2, 1) \n",
    "plt.plot(D,logP)\n",
    "plt.ylabel('logp')\n",
    "plt.xlabel('dimension')\n",
    "plt.title('d - logp')\n",
    "plt.grid(True)  \n",
    "\n",
    "plt.subplot(1, 2, 2) \n",
    "plt.plot(D,P)\n",
    "plt.ylabel('p')\n",
    "plt.xlabel('dimension')\n",
    "plt.title('d - p')\n",
    "plt.grid(True)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application 2\n",
    "\n",
    "(i) Construct a scatter plot showing the 1000 samples taken from disk $D_2$ with the above MCMC algorithm. Then do the same for delta = 0.01 and for delta = 20. Explain in a text cell why the resulting graph 'breaks down' when delta is \"too small\" or \"too large\"?\n",
    "\n",
    "(ii) Write a routine that prints for dimensions $d=2,3,...,100$ the MCMC-estimated volume of the unit ball in $d$ dimensions, its true value and the relative error of the estimate. You can start from the obvious relation $|D_1|=2$ and work inductively to compute ω(2), ω(3),... etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Application 2\n",
    "###############################################################################################################\n",
    "def mcmcplot(samples,delta):\n",
    "    N = 100         ## number of steps to approach equilibrium\n",
    "    point_x = []\n",
    "    point_y = []\n",
    "\n",
    "    for i in range(samples):\n",
    "        x = [0,0]   ## start at the centre of th disc. This variable will keep the position of the chain\n",
    "        R_sq = 0.0  ## this variable keeps the squared distance from 0. It saves some computations to keep it\n",
    "        for j in range(N):\n",
    "            k = r.choice([0,1])                    ## choose a jump direction at random\n",
    "            z = r.uniform(-delta,delta)            ## choose a jump size uniformly in (-delta,delta)\n",
    "            x_prop_k = x[k] + z                    ## propose a jump by z in the direction k\n",
    "            R_sqprop = R_sq - x[k]**2+ x_prop_k**2 ## compute the squared distance from 0 after the proposed jump \n",
    "            if R_sqprop < 1.0: \n",
    "                R_sq = R_sqprop\n",
    "                x[k]= x_prop_k                     ## if the proposed jump leads to a point in the disc, then jump\n",
    "        point_x.append(x[0])\n",
    "        point_y.append(x[1])\n",
    "    ## Plotting code:\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    plt.scatter(point_x, point_y, color = 'green', s = 1) # parameter s determines the size of each dot in the scatter plot\n",
    "\n",
    "    ax.set_xlim([-1, 1])\n",
    "    ax.set_ylim([-1, 1])\n",
    "    ax.set_aspect('equal')  # set aspect ratio 1:1\n",
    "\n",
    "    plt.show()\n",
    "###############################################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztfX+MXVd95+c4jo1mulvsvIZNQtxgwtql/zDgJkCkbmBnFTqRHKIm2iDomtlUESVesbUYN4hd3pjVqmBLoW2makhbJo6CQgAtu249lYOxTUSEE5w6wQnYZOKwG28i8h4BJN5oYwfO/nHvfT5zfe89v77nx31zPtLTe+/ec8/5nl+f+z3f7/nBOOdISEhIoMKq0AIkJCSMFhKpJCQkkCKRSkJCAikSqSQkJJAikUpCQgIpEqkkJCSQgoRUGGNfYoy9whh7puY+Y4z9FWNskTH2fcbYO4V72xhjz+WfbRTyJCQkhAOVpnI/gA803P8DAG/LP3cA+BsAYIytB9AFcC2AawB0GWPriGRKSEgIABJS4Zw/CuDVhiA3AXiAZzgK4I2MscsA3ADgm5zzVznnPwPwTTSTU0JCQuRY7SmdKwC8KPw/k1+ru34BGGN3INNyMD4+/q7Nmze7kTSBBK//+nX0l/rojHWwetXq4bVXBq8AAC4dv3R4PSE+PPnkk33O+W+ZPOurVlnFNd5w/cKLnN8H4D4A2LJlCz927BiddAkk6C/1MX98HtMT05g/Po+dB3fiE5OfwMx1MwCAPY/twc6DOwEAfzr5p8PrCfGBMfa/TZ/1RSpnAFwp/H8zgJfy69eXrh/xJFNCCSIpdMY62s8XRAIA0xPTy76L34OzA4Atv54wYuCck3wAXAXgmZp7NwL4J2SaybsBPJFfXw/gBQDr8s8LANbL0nrXu97FE+ix+zu7OWbBd39nt9HzvUGP7/7Obt4b9IglS/ANAMe4IReQaCqMsYeQaRwdxtgZZB6di3PSuhfAAoApAIsAlgBM5/deZYz9NwDfy6P6LOe8yeCb4BBV2oUOOmOdNKRJAOMt3Pog2VTigu2wKZY0Es6DMfYk53yLybNpRm2CEvpLfex5bA/6S/1lv4HztpT54/PO0veRRgINkk9vxODqjS4aYQEMf89cN2M9bFJBUxpJi4kLiVRGDGLnp7Rv1HlzXHToqjjL9poq9/Xg7ADja8aNZUnkRINEKiMGV1pDuVNXzT2hIjEVYqxyXw/ODSqfUyULV4S80pBIZcTg2wPjgsRU4hTDFHkubDyDs4PhbF5AnSx8DONWBEx90SE/K3GeSpoDooaquTaysktleyFgMU8leX9aAhXvR9krEytcyjk9MY3dk7uXaRuFJlM39EmeJVqk4U9LoKKax2YTqLNluJLT1NAqK9tkwNWEqYoT8rMShz8qiE2NF4ciomxVclLIbrvMwHe8MQOhp+knxIHYpsmLGkBZOynLSaG9uDK0JgOuJkzZKORnJWsqPrQRF2lQG0tj08pGDUiG2pWB/lIf276xjcSo2GQsLbSGbd/YRmZMLRtLy+nLjKl1MtqUQ1sM221DIpUWYf74PBYWFzB19ZS1Kt7UKacnpjF19RQWFhcwf3wep/qncOOXb8Sp/qna+HQ7qC0pVHl5dDFKHrWoYKrihPyM8vCnSa2nVPnr4iqun+ydHN6fenCKYxZ86sGp2vh0jZm+8mkqQ4GVaKTl3G74E5wgTD6jTCqhG3FV+id7J/nUg1P8ZO9k7XMykqrz+qjKEAor1XZjQyrJ+xMZdD0NunMoZOGr0t/U2YT9H94vja9qvc3g3AC7vr1reL1t0+Vj86i1Acmm4hH9pT5mD89i9shs7RhdZtAso8ku0F/qY/bILGYPn09PZkegMpgOXcQcQ9vH9MQ0ur/fxeDcoNFG0RnrDN3QrmwZrm0lK9oWY6rihPy0dfhTqPVVE8Jkz9QNBZriKKcnC28ClWGPTn5UwsU8Uc5X/K6B0DYVZAeAnUK2B+1dFfe/AOCp/PMjAD8X7v1KuLdPJb22kkpv0OPdQ13ePdwddgxZw7PpQL1Bj3cPd3n3UDcam4BqflTIUpWYbeQwfbbttpigpALgIgDPA9gIYA2ApwG8vSH8fwLwJeH/L3XTbCuplOGy4bW9UTdBzBuVRqBTXkWaUw9OjWT5cm5HKhQ2lWsALHLOT3POzwL4CrJjTuvwIQAPEaTbeujaL3SgOw+krTYAivkqwPLykpVFeR5PwnJQeH+qji69tiogY+y3AbwFwCHh8hsYY8cAvA7gc5zz/0kgUzBQrWi1jUfXgxLbCucmFLI+8vwjuG7Dddh+zXbtMirKd+umrdh3ah+2btoKoHqdUhmdsQ723rx3WD8JJZiqOMUHwK0A/k74/0cA7qkJ+2flewAuz783AvgxgLfWPHsHgGMAjm3YsIFc3aMClYHRt6HPp53Adh2QOCHPtIzEIYzupk4rAQg8T6XuSNMq3AbgTvEC5/yl/Ps0Y+wIgAlkNhqUwi07S9laakdo0hB0tAHfczVU5mNQaTOyeGQbWReawtzjc8ZHqBZHsC6dW8LvXf57lZs6JRjClI34eQ1iNYDTyIY1haH2dyvCbUKmiTDh2joAa/PfHQDPocHIW3zaaqi1XYmr676lBsVK47IHrPJ+7rHqHu4mt28gIAKX8hQyV/HzAD6dX/ssgK1CmFlkNhPxufcCOJET0QkAt6uk11ZS0UW50dd1glg6h4ocsjAU7mJVxOhyjwXBScX3JzSpuGrsMs0ktKYig45dRlywKEJlnRElYiHk2JBIxTMoG6KLORc+YLslZN1cjxgN1CsRiVQ8QfaWNYGOuh9TB6jq/DqEIHpwbD0vMZXLqCCRiie4mFavM2ToHupGo8nYaiom4esQm4Y3CiSXSMUTQm3qU8TZPdxtbWN1qYG4Xseji9hIzgSJVBzA1UK10Op9qCGWSUfzscbGBQH4Nja7QCIVB3D1tgn9FtNx6VLClEyr7C6h5ZKhrgzbNCxKpOIAvtzGvuFKU7HNVwh3ues6Lhv0Q79QdJBIpYUITS4U0J392pTnEB3OZZpVWlZT/k3vuUIilQghawjlBt1GkinygFnw7iG5EbmpE4fQSFxqKuIJBDrzdqrKJgTh2pBK2vi6ATbbD8gWzZUXDFIs1usv9TH3xBzAge3X6m8HoIvpiWkMzg2U02taJOlyEV9d2bpKUzyfae/Ne5XqoalsinIenM329o3+kHhTNgr58aWp2LwhQszZEDUH12tnxKFP7NqVzlwgHx42E/jWVpCGP24Q45BENvYWO7rLhlgmMGr4LvsYjahiGfguj0QqKwg67koXywrEuF1qKi46uS9jqKmHrXw9JNElUlkBEAmiaj8Sn4Y+H29NCte3aielzo/pXKAm433SVFYIqVBUdG/QvGlRgcJlW0zVLzc+2cZHPjuNDGV5TORT8SSV3dx16fgmXRNvVLKprBBSUaloVdeyrPEXiwlFe0koNdn2DVqW10R+lQ6o4uZuiismm5pvrSWRiiFsK0dljYfYYarsEFVaRlUno7AHuPCCmMy3odBUmkAVHxVZ65CWiuw+XiKJVAyh23lVni9DjE/VY+LqraQir26DDa01yWBTflQLA1VtKE1hRSRNJWJSqaocnU6hW7llTcXVW6npzSjz2OjkyfRN63NYYUNysWoqPhCcVCA/S/mjAHo4f2byHwv3tiHbRf85ANtU0rMlFZWhhAs3bBmu3kpN8VJqEqZx+Xwb28QTsoObpk0lc1BSgcJZyjmpzFU8ux7Z8R7rkR3XcRrAOlmatqTiYhhQBd9eGZV4KdN02fCrbDUqXjKXcvuES8JWQWhSeQ+AA8L/TwH4VClMHal8CMAXhf9fBPAhWZouNRWdMDKo2lBcoA0dh3P9iWCFW90Gsdl9qtBkzPbx0ghNKrfgwmNP50phPgrgZQDfB/B1AFfm1z8J4L8I4f4rgE/WpBPs2FObN7LO25VyjN2GjsO5upy9QY9P7p3kmAWfOTBj1cFiIlxVWcRyqiszynzZkArFKmVWca18LOk/AHiIc/4aY+xjAPYCeL/is9nFgMeemq4g7ox1MPu+Wat06tKWraD2fWyqikxVUJWzM9bBQ7c8hPnj8xicGwzLBEBt3fhenWyCsozlMhQPkgfOl9Pg7ACDc8tXLVMdS2sNUzbi5zUI6fCnFP4iAL/gHoc/tgxOOZ6XpaOqqVBrIpTDPdfakWhMb/JmxaSRiGgaypTLUMcdTbk3LgIPf6RnKQO4TPh9M4Cj+e/1AF5AZqRdl/9eL0tTl1SKCjDdQLk3cL9Xqm6DoO4wVIZpW/LWeb4tQ7wymuRusqU0hZPFq4ugpJKl33yWMoA/B/BsTjiHAWwWnv2PyFzRiwCmVdIz0VREUlA1ehWwJSUViDuF+YJuOVCmVwXdThGrJiJDCK+fLoKTiu+PifdHLHAVo5f4nI/NiFQ1FcqGo5J3yslZJulRo21EFEreRCqaEDsw9dvTNSjlMcm7TfqiLSRUx46tPgvU2e2o5NUlp0QqmtCpqNjebD7lodRU6jRFnbRtoEpooeq7bk6TTB5VecWtM1SQSEUTVA3Hl1fIBLGRoUgkOnYslfkrKnYh1fhCaTJVbcmknOqG0eLWGSpIpBII5beLT5KRTfqq81aFIhvddCknhXHuzmblsjx1NbreoMc337O50uCfhj+EpGJT6Sqqp0gidSqsCzQ1uCZvVUibgosOSK2p6MJleeqWV6GNXP2XV1vPVUmk0gBdtU+Eris5Jk1Fds/EWGpLCm0kNCqbhg/o2k2akEilATMHZobrRapgOozQiSdGmKjWthMAVcqIuvPrXi/Dp1HZNh5KAkyk0gCZpkI1dyKkgc/FG5jzanuFywmA5TQpntO9XoYKKVERj0nedepfJ/5EKg1wrb7aDCdkcVI3FBs5fGliqumUw7nSVOrgwptlsnZHp/6TpkJEKjZw9XaRwVVDUYFLwzZlPCHtM5y78RKZ1LuriYSJVCxg6zGw1YSqGofK29WV5mDTWak6uk87hi9Q25Nsh28yJFKxQFPlUDRcWeWLtgqdLS4pOnBV/mzewLEYLClBJQulZtUbqHsZTdNNpGIB1w24PE6uswWoqLE2HbgqPEVDd03KlPGYQLeM6oiacqqBj+FxIhVCmBoK61BuAKFsAVWGRYrxuMrwUfWkQBXZbWDSwXSfqZKVus5dOx84T6RSC5PCFYcjFEY1maZCJbcMYpy+iK1Iq3ymsWk8KuVhaiNreq58r8lLQzGk1M1XGUMyP2xO5olUamDSeXoDtUleppqKK7l1YNvIdTu5qupvKpcqYZoSTvlesQH35N5JK7l1YDLkKeZombSjRCo1oGikNjAdS7tupLbx6zRwV2HrnrOt86ohYTnOYpb2nfvvtO68uvL5akeJVBQQwtgXyn4ig61cVMMRm7AUz1VBpWzKmkCxM2CMW2CUoVpWwUkF8mNPdwD4AbJzf74F4LeFe7/C+eNQ96mkV0UqssJy0cF9GMxsZQglly5sPFuu5NAJG9MLxHSYJyL0bvoqx56+D8BY/vtPADws3PulbppVpCIrLIo3Zvm6qlHXJSgas29DYhXEfMTUQVXhy06lEs7UriQiNKnonvszAeAx4T8JqfhQgcvXVY26VaC029jGQ9GJKYdUMWpSTfBZBzrDMxt5QpOK9NjTUvg5LD/q9HVkx5keBfDBhue8HXvapKmYbPlXharGEapjudZUdOOPlVTq5PKpLZoa/3XTC00qt1aQyj01YT+Sk8da4drl+fdGAD8G8FZZmiG3k1RtQCb2ltiHAJQEqhPe19BCVa7ycNc3CfrQLEOTitLwB8AkgB8CuLQhrvsB3CJLUyQV3xWqmp5JxVNoKi7Lw7Qxq9qo6q7b2gioCNpmuFsXX6g6jl1TWQ35sacTuTH3baXr6wqtBUAHwHNlI2/VRySVGN/onOut6Sk/E9o+UoXegP5QNQqtz5eNwUVctnWlc36VLmJwKcuOPT0I4Cdl1zGA9wI4kRPRCQC3q6QXUlPRhc5blkLtNyEz23w0ydGUfgxep5CwlV08Kpf6ZRKcVHx/Yjmiw7bjyEjEpqFQNzLdDkBlexoFuMrjSGsqvj8mB7SrWtZ13vIUbtQmS75NQ/GlKdjGHevwlRKUNh1fBJxIRYK6Si2TSLGaVjTGUbpKdWTzgRg6dAzE5oK8KYzuZbi0l5XlS6RiUGicX7jr2swjM3zqwSn+3f/z3WWE4rLjmWx2TIVRH3qo1p2LYaaJAdlkGoIpxLiqZE2kYogLNJWK1abUHa/JbhJDJ9d1/6reDwHdo04pDljTiavcmW3bgs4zsrQSqRDBZKhja8CUvTFs5daF6Zs3hqFUGboyycKrEKvOPJYmTcWkPHWekeVlRZOKrzekjZorwoWNppDBdtvGJhlcaSox2FRUw6sSK8UCU9eaigwrmlR8vSGpNBUXKGSw3bYxBGLUcOrQxiGgKWxIhWXPtwtbtmzhx44dAwD0l/qYPz6P6YlpdMY6gSULizaWRQiZ21hOvsEYe5JzvsXk2VXUwvhGZ6yDmetmAAB7HtuD/lLfKr7+Up8knhBpF2XRpo4SQub54/PYeXAn5o/Pe0szBFTbkxiu+I1VWG2arvGDsaFoKACGJBMynralHTOoNYvpiell36MK1fYkhgOQ/R7HJcYJm46bQn5sNmmSufxcTFpqkkOMO+SclZhha3exrcc22kp6A7UFoL3B8lndQ8/TKjzFV6qhVhflCW9NDdWlEbEq7tiNlqE6l226tuVKVS8+y69p6kJTuAJYyd4fXVRpKtSeHZXnqsLE7GHinP48YF950ilr1edN4PKl0TTnpSntOs09kYolqCs79viqQE2EMoTWykKk75JIbSfulZ9PpGIJ6sqmjs90ibvptG0dmYqDtXQ7Z2jtS/ZmV33WJh+UZWBrUyw/n0hlxCF2eJ3OTzltuxx28z2bs6M/H5iMYmhmC9Wy6g2WT8W30XhCaEtlm2Ld7N9EKhoI/YY0kcP0zegqr0XD3HzP5pHxVNV53upsE0VnjEVT0U2zyG8dqSVS0YDK28FHZYe2KdggFmKmhFgfopu1vHLddNgUY5k1yRScVCA/9nQtgIfz+48DuEq496n8+ikAN6ikZ6upiP77qoKNxTDa5vTahPJcjaL+MZudk6xq4Gwq47a9RIKSCtSOPf04gHvz37chP/YUwNvz8GuR7cb/PICLZGmakkpR6eLbp6pRUG8aTQEbwyLn/Pxh4oe6tXGGQuhyr5rToXpgl1iGTcQRS1mrIjSpSM/9AXAAwHvy36sB9AGwclgxXNPHdEZtUeni20ccSxerfLuHu1oVYAqdzlRusLpvvqq8+Xx7NtVPQXiTeyeDvM1tO3xoUhRlsE17aGtZjRPckBMo1v5cAeBF4f8ZANfWheGcv84Y+wWAS/LrR0vPXlGVCGPsDmRHn2LDhg0X3FdZ5yCu+SjWkMwfn8fC4gKuv+p6oFiwrbBwm2I9SiHzkR8fwcLigrLsVd8ybL9mO8YvHl8W3sUamLpyaawfln1t6mzCmovWYOumrWTyqEBclGpSr+W8FYvydNuGbtr9pT7mHp8blt+ub+8aymCKHQd2ZG3xN3GlcSSmbFR8oHDsKYBnAbxZ+P88MlL5awAfEa7/PYA/lKVJeUC7qWel7AUwQQxvOGrUaT9NE+eatvP0DRPtrc5DpJsP3Ta1zPZDsEEX5zSaykgMf0IYPbuHumTqeujxtutJWHXxl13TocuBc5qyMPUIlefAqKSjYvuRyVN1Hy049vROLDfUfjX//btYbqg9DQNDraoFvqkQTTqDzDOgkm5Z/hBwnX6T9lJMopt6cMpJ2mJa1IRlG2dVubiQU1a/VfeDkkqWvvTY0zcA+Boy1/ETADYKz346f+4UgD9QSa9JU1HpIFVhVK+V09MBdXxUDdDWs6Qbv4iTvZN88oFJPnNgxukw0AVx2sbpSzOTDbOj01RCfJpcynUVJbOd6GgqOumK92X7W+iqqU3eKht7jW/NSbQl+FrF6zLOGIZxVXLo1GsiFQUMhywlgxZVA5C5o001qKb7VXNPymFNOqrvzhKrwdpWI6XYVd8G5faiU6+JVBRQFJx4tKnqcEkFTR1cTJ/SoNYU3qajxmr78Q3T/OoaXGVxmRKtKolEaVPx/bGdpq96VnL5uabhkU08KjL7eoM3dQjbIaKuHKG1FhcdWje8D1tQ0lQsSaWuEJuuc063/aOvZ0SoakndQ10++cCk1ipcV9qLL89ILDAx4rs08CdSIUJTB6F6S9s+Y/K8jj2nbAuQPWvb0XUI3qTj+QBF+iabnlMSejkPiVSIQDmEoW7oQ5uQ5qzT3uBCr1OVbaZ7uMtnDsxcMJHK9ZBHp2PIhgihDKMmnbucF5N1Zy7rIZGKZ1B4cprQ9JbWnXBXNX4u20x0ZXX5hjSNw6U7WiV9HQM753qePFMZdDTcsmE/HdGhUdA24XTC23QWSntCnaVffKu7yL+N3LEZs23jVqlPF2mI11RfBMNh8L/AizyRilpB24TzBaoO4mrookOqJofG2761deHaGO3D5mOrqZSfSZpKTaGqGr8Kz4fupjxlxHbCoE/PTF0Yk8mGvve1cWWMDmlA7g16fOaRGT65d3LYHnXkSTYVbu/ytdFqisoqVi27XhwnQ3l8HFJTqTOqNnVgVYKnAmXnF+MKqQEPhzFCexTlkWk2iVS4vapvY38pKmvmwAyJpuJifF0g1O52dZrKBV4Qz0Mfasg6LjXq2n33cJffuf/OWk1F9hJOpBIY1AZhl2N4351WlpcLvCCeSY8avg3Gqhq6imE4aSoWCDXWVSWLmD0Z1OlRe0HaCpV8q3ruVMmnKc1EKprwOdY1scDbwsZg3MZO3WRDcpUWdRqq2oXKfBwboimQSKWm0Orgs+NUvV1UYCOjeKSlLmJzr6tALGPXK4Opy6eJEHW1C10Zk6YiIRWKytbtyKpekKrGLnvWtEFx7l5ToSZlW22OWlNpakuFEbRpsy3KtEzyQ1E/iVQ4TUHqEpPNG0HFtdr0TCjPgkx2E7nE+GLQlEwI31VaoRCMVACsB/BNAM/l3+sqwrwDwHfzYzq+D+DfC/fuB/ACgKfyzztU0qXw/qiOO3XjsElf5xkfna8uDerhQQi7kw3aIGMdVGUPSSq7kZ+dDOAuAJ+vCPOvAbwt/305gJcBvJGfJ5VbdNOlIBUqt22oLRBDaiq+5WqKU1YPIYZzMUO13YcklVMALst/XwbglMIzTwsk44RUVLQQW1cmtXFwVOF6WCqrB5VORKX16WhcKkTpguTaoKn8vPT/Z5Lw1wD4IYBV/DypnMqHRV8AsLbh2TsAHANwbMOGDY0FUtVIZA1Ht2GJb0jf08pdgrpBm3ZY1Q7qSlMxeU7HNqRClCFfVE5JBcBBAM9UfG7SIZVCkwHw7tI1huwwsb0APqMiNIWmovKMCuq8O20FdYM2LdfQHUvFGF1egd0GTUVMo+llGP3wB8C/BPDPAG5tiOt6AP+okq6PafqqFSuq36OkqYS2FYU8rkPmNi7qnOr84nLatg4EFRR5iG7yG4A9JUPt7oowawB8C8B/rrhXEBID8BcAPqeSrgqp2FaEjbt4FODCzuAzXRvIZHBZ5yZDdxPErKlckhPGc/n3+vz6FgB/l//+CIBzgtt46DoGcAjAiXw49SCA31BJV4VUbDWIUGRBna7NsI5CDp0OUW7oIeuAcoKbTj6oNBVbr1cwUgn1UdVUTGayhgb1Wyn0m1+nvMuy+pC9Tj7KtMW4XNn2mtI0CZNIpQZVFRS6k8kQQlOJhWht3fxNcdWhrj1QyFJlH5IRjM9hZ9JUDEilCraFHUsHpEQo+wkVbDqnqsfGpIyqnhGHeFXHrcTSvhKpEKOuMVS5Eavg0o7hotGFsJ9Qoqm+dPKk4kbWedFUaTvdw93htqO6x62oykOBRCoaMO24Om5E085lOw4OjZg0FdfxiGFV60R048ocCDJZTNuBah4TqSgWWp3xVgW2FnyqNKi8A7Ggyu5QF8ZH/lRtKTrGVzFuVa+SjDRMy0SVjBKpKBZacc/VRDVflv0yqLUXl524HLdYJ1TT1m3kL6elasi1TbcMV3WQNBUHmoorDYJz+vVF1PKpwuUQqxy3qKnUTcbSzZ+N/DpenzqCdKkFm0I3jUQqjmHiSdC9H9MQxqemIoKKDHy95esI0tWQxAa6aSRScQzXGg7nbnYTC7l+pgmq3hIdhOiYVO2iN3B/gFrSVCIiFRvtQ6ehU75dC7f35AOTTm1IpnBBANTaSZ3NRNb5VeSoyr+Ll4pNWdiQymqMMPpLfcwfn8f0xDQ6Yx2jOOaPz2PnwZ0AgJnrZpTv95f66A16mNw4ia2btkrT6Yx1KuM3As++Jt40gTWr1mBhcQHbvrENe2/ea1wOlJiemF72TQHS8kN1vXbGOhhfM46dB3di/OLxxvYwODfA+MXjlW2vKv+UZSJrs85hykYhP6qaCgX7m2oq4pwEW3XZRmYbN3ooxGBfqrPRyAy5xf+q2bIhZDcF0vCnGiE6sXitrCr7GPeryhczqAzjpmFl8ohx6bicKWTxhUQqgeGiE/gElVzU8VQZmlU6dBUo3cxiXCZ5DvFySYbalpEKJVmEIB7bRj5U+RXWRVU9V5fXQi5xaYRph3ZZR7px+3B/l5Fcyi0jFc7bvSjP1gVd1fl1nqvLaxVZ2Zy+qAOddGzqrMhb93DXQEp1+NRURtr74wOFh2lwdoBdj+4CYGdxp/aMqHjACs/Jnsf2GHkNRJl1vEtNeRXlBjD0pMwfn8fC4gKuv+p6zHTsPBtNZbPjwA4sLC4AAPZ/eL9WPrS8jrz07QjU3rFGmLJRyA+lpmKrYQzf0qVl7LFMQAs1V8YWJgZQqjQ4tzuP2keZu64rxHzsaR7uVzi/P+0+4fpbADyeP/8wgDUq6VKSCpU9oc6l7MKdS2lLcN04Y+00LtPwYfh2PUwOSSrSY0/ze7+suf5VALflv+8F8Ccq6cakqcjidaGpUDYo141T1W5S5+GpC7MS0FR2KmViU24hSUX13J8LSAXZsRx9AKvz/+8BcEAl3RgNtTqoq2zV6769TTbpqXp4xI5TvhZqfk9oUA3NTcotJKlKqoY6AAAMqklEQVQonVAI4HVkR5YeBfDB/FoHwKIQ5koAzzSkpXzsqQ80uQVNOpLJdR0ZbSAO5XwMFdquqcQib7SaCgiOPQVwef69EcCPAbwVwG9VkMoJFaFj0FTEjq77ZrXVVHRltNkLtUhf1TYUS4cqwydRiSQc2lBviuiHP6Vn7gdwS9uHP73B+a0Byw0nho4lykexDkU1T7EOVWTETz2kLEi4beuuCoQkFZVjT9cBWJv/7uSenrfn/79WMtR+XCXdGEiFcz9uT1OIsvn0APnwJpmcHijTVEy1S1l6SVPRJxWVY0/fi+xo06fz79uF5zcCeALAYk4wa1XSjYVUZO7kUG7kIjzVJssxoZDVd/lSlFEMLxtVBCOVUJ+YNmmiekaGqkZN1RHa1NhNNRWKdG09YLpro1zL1YREKkSoqiCdYYRLyGSrCmuqelPls01k5RJDg7nm2iiduKm1zEQqRCgMmt1D3eE1sWOYVKDLjmUz/GqSi6qhtmlY5RIh2oAtbEglLSgUwUrfJZgs9nO5tV/dIjEVOZvkslnUKC6m04mHYutPCriQQ2cxn276XhcKqsKUjUJ+Qg1/qOKkhk4arr0SpuUVi1ajKwd1/cZSDkjDH3cIbahViUunIbputKZ5j8X+oisHdXmquv9du6oTqTiCygzXqvkPTROeXDRaE02FOmxbYZvH3sDNmT0ye5nrSXWJVByhrkOL1+tmatatk9GdZKXS6FWeUb2mI6sJYiMqijxST5xrijNpKi0nFVNNRUV9pezIdcQmxqF6TSarLSn4njtjWx8UacTmNVRBIpUIEHJuhytNpQq+DNdiOiayF/ddTDjTRWiCMEEilQhANVSI3eahor1Rp2OiZRX3XUw4WwlIpBIBqDpVVWeRGe1Cux9dy2KjqSQyMUMiFQ/w1Uir0pEZ7WLoOCYGZsr0EmiRSMUDQmoFbexArssrFi3Npm5cGMKpkEjFA3QbQCyNQweq2oZKvqldnr41IVXYkJuJrcgXEqkEQlMDUGkcsXSMAiqu6d6geXKfq04RwhDuOr6kqUT0CUUqOm9LlcOoYnkrFVDJXyFz3eQ+V52CyusUW5nHikQqnqDTIKu2USjDx8Qr6jRjeZMWaDJidw93+cwjM8s2dAolv5hubGVYBRtSSVsfaEBrSwDJNgpAtmx9emIac0/MARzYfu127eX2ulsrlMO3fal9XZ3MH5/Hrm/vGv4fXzOOmetmnG1DIAsrljsAZ9thRAFTNsrITH7sKYD34fyRp08B+H84f/bP/QBeEO69QyXdWGwqTVA9i7d405qq5LaaStuGA6r5rdNUdKBTNrKwK0lTsSUVpWNPhfDrAbwKYIyfJ5VbdNNtA6moNsii8Zs2/JDDmXJH8bF3rMymI5OzfL1phTHF7OY2EEgVQpKK1rk/yE4Z/LLwf2RJxVdjCqlpiGnbalyqkHmfZHJWXac4cE037dhhQyq2NpU3cc5fzodRLzPGLpWEvw3A3aVr/50x9hlkR3zcxTl/zVKmKODL9mCz9SN12oNzA4Bn/2U2BtNtGztjHey9ee/wWRM5xeuDs4PM7sXd2DlC1k8wyFgHBMee8vOaTA/AxaVrDMBaAHsBfKbh+ajOUm4zfGhRqgv+YnmDt3WY4gpow/AHwCcA3Ndw/3oA/6iSbhuGPzHDpENTu65duMJjmtjWdoQkFemxp0LYowDeV7pWEBID8BcAPqeSbiIVNVAaD0NrFirpU8sYOs8hEZJUpMee5v+vAvB/AawqPX8I2VGozwB4EMBvqKSbSEUNlFPbXezDKkuzbnavL09L0lQCkEqoTyIVOSiJwOcbu+jITTu26cwJSTCDDamsMjHuJtihv9THnsf2oL/UN7qvgvnj89j16C6MXzx+gXdFN/7piWnsntztxYMx98Qcdh7ciaVzS7VpyuQpZq/OH5+vvE9Rvj7ibCvSNP0AkE2tpzjVsM6V2V/qY9s3tmFhcUE5fq9T83n2NbZ6rDZNmTwyN66LUyNdnkTZOpiqOCE/bR/+UHtGdFA3IzWWZfjUWwn4CKcbtg1AsqmsTJg05Lo1SVV2Ct/eD9uOmbw/dEikMmJQ7VxNy/7rntd5xvfbt0q2kNpCkwdq1JFIZcSgsxixqpE3PR9zx6iSzZW2YOIdW0maSyKVEUPIlcch4/aZnskCyJgJmRqJVBK8oeiMTat6YzBwyvazqdJUVhJpyGBDKmmeSoIU4hyMYo5Isaq3ai7I3OPZXJO5x+ekcVfNKSnP+TCZ17PjwA4sLC5gx4Edlc90xjrYfu12jF883ihLgj7SPJWERlTNaym2oRxfM149F0RhK80CVXNKynM+TOb13H3D3cu+q1B+bkVuU+ACpipOyE8a/tBDZvSl2GnNVBZX83rScKceSDaVBNsOUqy16R7uksSrsgCQSvYEetiQShr+jAisp4nz0ncO0yn6OrvHpynuo4VEKiMCW3vA9mu319tIiOSpi1tHdtNtKBP8gWWaTruwZcsWfuzYsdBiRI1R7XyzR2ax69u70P03XcxePxtanJEFY+xJzvkWk2eTS3lEYeoePdU/hRu/fCOOvng0zqX8NcO0hHiQhj8jCtPhUDG/4/TPTuPkT08CiMvOQT1MS6BH0lRGFIWBVXfoc/cNd2Pq6il84YYvYOrqKWzdtBVAPBsbmeYrwR8SqSQsw6bOJuz/8H6ceOUEFhYXsO/UPgBuZpuqzKZNaB+shj+MsVsBzAL4HQDXcM4rraeMsQ8A+EsAFyHbEPtz+fW3APgKsuNQ/xnAH3HOz9rIlECD8vDJxWxTldm0Ce2DlfeHMfY7AH4N4IsAPllFKoyxiwD8CMC/A3AGwPcAfIhz/gPG2FcB/A/O+VcYY/cCeJpz/jeydJP3Z3TRX+pj7ok5gGf2kzTMCYNg3h/O+Q8556ckwa4BsMg5P51rIV8BcBNjjAF4P4Cv5+H2AvigjTwJ7UdnrIPxi8ex69FdaWFfS+HD+3MFgBeF/2cAXIvszKCfc85fF65fURcJY+wOZEefAsBrjLFnHMgaGh0Ao2pMUM/bKqzGOC7Z+dmdP935652vyx8IilGts02mD0pJhTF2EMC/qrj1ac75/1JIo2qtKm+4XgnO+X0A7stlOmaqmsWMUc0XMLp5G+V8mT4rJRXO+aRp5DnOALhS+P9mAC8hY/c3MsZW59pKcT0hIaHF8OFS/h6AtzHG3sIYWwPgNgD78pWQhwHckofbBkBF80lISIgYVqTCGLuZMXYGwHsA7GeMHcivX84YWwCAXAvZDuAAgB8C+Crn/Nk8ij8DsIMxtojMxvL3iknfZyN3xBjVfAGjm7eUrxJauaAwISEhXqQZtQkJCaRIpJKQkECKVpAKY+xWxtizjLFfM8Zq3XeMsQ8wxk4xxhYZY3f5lNEEjLH1jLFvMsaey7/X1YT7FWPsqfyzz7ecqpCVP2NsLWPs4fz+44yxq/xLaQaFvH2UMdYT6umPQ8ipA8bYlxhjr9TN+WIZ/irP8/cZY+9Uith0H0qfH2RrizYBOAJgS02YiwA8D2AjgDUAngbw9tCyS/K1G8Bd+e+7AHy+JtwvQ8uqkBdp+QP4OIB789+3AXg4tNyEefsogLnQsmrm6/cBvBPAMzX3pwD8E7I5Ze8G8LhKvK3QVLjFcgD30lnhJmTLE4D2L1NQKX8xv18H8G/z5Rqxo41tSwrO+aMAXm0IchOAB3iGo8jmlV0mi7cVpKKIquUAtdP+I8GbOOcvA0D+fWlNuDcwxo4xxo4yxmIlHpXyH4bh2VSDXyCbShA7VNvWH+bDhK8zxq6suN82GPWpaHZ+c7gcICia8qURzQbO+UuMsY0ADjHGTnDOn6eRkAwq5R9lHSlARe5/APAQ5/w1xtjHkGlk73cumVsY1Vc0pMLdLQcIiqZ8McZ+whi7jHP+cq5WvlITx0v592nG2BEAE8jG+DFBpfyLMGcYY6sB/Caa1e9YIM0b5/ynwt+/BfB5D3K5hlGfGqXhT+VygMAyybAP2fIEoGaZAmNsHWNsbf67A+A6AD/wJqE6VMpfzO8tAA7x3CIYOaR5K9katiKbPd527APwH3Iv0LsB/KIYrjcitAVa0Up9MzLWfA3ATwAcyK9fDmChZK3+EbK3+KdDy62Qr0sAfAvAc/n3+vz6FmQ75AHAewGcQOZxOAHg9tByN+TngvIH8FkAW/PfbwDwNQCLAJ4AsDG0zIR5+3MAz+b1dBjA5tAyK+TpIQAvAziX96/bAXwMwMfy+wzAX+d5PoEaz2v5k6bpJyQkkGKUhj8JCQkRIJFKQkICKRKpJCQkkCKRSkJCAikSqSQkJJAikUpCQgIpEqkkJCSQ4v8D/4d2r1PsEKUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFKZJREFUeJzt3XtwlfWdx/H3lyBgWBUk1FIxisVB3XYrbbx0mWqreKPT2G5ta2Z21IyO067szqxT1LbbBbTOeunojoO7Sq3pZTpoy06ncYqrRERmrKGGCYqXRgO4yEKBALo2IVyS7/7x/E76EHI5yfnl3Pp5zZw55/yey/k9Ifnwe57zPM/X3B0RkVjGFboDIlJeFCoiEpVCRUSiUqiISFQKFRGJSqEiIlFFCRUze8LMdpvZ64NMNzN72Mzazew1M/t0atoNZvZOeNwQoz8iUjixRio/Aa4aYvrVwFnhcQvwnwBmdjKwGLgQuABYbGZTI/VJRAogSqi4+zpg3xCzXAP8zBPNwBQzmwFcCax2933uvh9YzdDhJCJFbnyePudU4L3U++2hbbD2Y5jZLSSjHCZPnvyZs88+e2x6KiJs2LChw92nj2bZfIWKDdDmQ7Qf2+i+HFgOUFNT4y0tLfF6JyJHMbP/Ge2y+fr2ZztwWur9TGDHEO0iUqLyFSqNwPXhW6CLgA/cfSfwLHCFmU0NB2ivCG0iUqKi7P6Y2Qrg80CVmW0n+UbnOAB3fxRYBSwA2oEuoD5M22dmdwOvhFXd5e5DHfAVkSIXJVTcvW6Y6Q7cOsi0J4AnYvRDRApPZ9SKSFQKFRGJSqEiIlEpVEQkKoWKiESlUBGRqBQqIhKVQkVEolKoiEhUChURiUqhIiJRKVREJCqFiohEpVARkagUKiISlUJFRKJSqIhIVAoVEYkqVtnTq8ysLZQ1vXOA6Q+Z2cbweNvM3k9N60lNa4zRHxEpnJzvUWtmFcAjwOUkJTdeMbNGd38zM4+7/3Nq/n8E5qZWccDdz8u1HyJSHGKMVC4A2t19i7sfAp4kKXM6mDpgRYTPFZEiFCNURlK69HRgFrAm1TzJzFrMrNnMvhyhPyJSQDFKdGRduhS4Dljp7j2ptmp332FmZwJrzGyTu28+5kNStZSrq6tz7bOIjJEYI5WRlC69jn67Pu6+IzxvAdZy9PGW9HzL3b3G3WumTx9V3WgRyYMYofIKcJaZzTKzCSTBccy3OGY2B5gKvJxqm2pmE8PrKmAe8Gb/ZUWkdOS8++PuR8xsIUkN5ArgCXd/w8zuAlrcPRMwdcCToVphxjnAY2bWSxJw96a/NRKR0mNH/42XhpqaGm9paSl0N0TKlpltcPea0SyrM2pFJCqFiohEpVARkagUKiISlUJFRKJSqIhIVAoVEYlKoSIiUSlURCQqhYqIRKVQEZGoFCoiEpVCRUSiUqiISFQKFRGJSqEiIlEpVEQkKoWKiESlUBGRqPJVS/lGM9uTqpl8c2raDWb2TnjcEKM/IlI4eamlHDzl7gv7LXsysBioISlAtiEsuz/XfolIYRSilnLalcBqd98XgmQ1cFWEPolIgeSzlvJXzew1M1tpZpmKhiOpw3xLqLncsmfPngjdFpGxECNUsqml/DRwhrv/DdAE/HQEyyaNKnsqUhLyUkvZ3fe6+8Hw9kfAZ7JdVkRKS15qKZvZjNTbWuCt8PpZ4IpQU3kqcEVoE5ESla9ayv9kZrXAEWAfcGNYdp+Z3U0STAB3ufu+XPskIoWjWsoicgzVUhaRoqFQEZGoFCoiEpVCRUSiUqiISFQKFRGJSqEiIlEpVEQkKoWKiESlUBGRqBQqIhKVQkVEolKoiEhUChURiUqhIiJRKVREJCqFiohEpVARkajyVfb0NjN7M9T9ed7MTk9N60mVQ23sv6yIlJZ8lT1tBWrcvcvMvgXcD3wjTDvg7ufl2g8RKQ55KXvq7i+4e1d420xS30dEylA+y55m3AQ8k3o/KZQzbTazLw+2kMqeipSGnHd/GEHpUjP7e6AGuCTVXO3uO8zsTGCNmW1y983HrNB9ObAckhIduXdbRMZCXsqeApjZfOB7QG2qBCruviM8bwHWAnMj9ElECiRfZU/nAo+RBMruVPtUM5sYXlcB84D0AV4RKTH5Knv6APBXwK/MDGCbu9cC5wCPmVkvScDd2+9bIxEpMSp7KiLHUNlTESkaChURiUqhIiJRKVREJCqFiohEpVCRguro6uCBlx6go6uj0F2RSBQqUlANrQ3c3nQ7Da0Nhe6KRKJQkbwYaETS0dVB5+FOFl+8mPq59QXsncSkUJG8WLZ+Gbc33U7dyrq+YGlobWDpi0uZPGEyVZVVBe6hxKJQkagGO0bSdSS5nU7T1qa+XZ36ufXcP/9+jVLKjEJFohpoRAJQeVwlAPNnze8LkarKKurn1tPQ2qADtWUkxv1UpIx1dHXQ0NpA/dz6rHZRug7/eUSybP2y5G47DnWfrKPrUBetu1rZ27WXqsoqOro6qFtZR9PWJjoPdbLkC0vGdmMkLxQqMqTMtzMAi+YtGnLe5veaebz1cQAmjpvIoxseZVfnLgCe2/Icb+99m70H9lL/m3r2H9jPZbMuo2lrEwBr313LkheWsPDChTq+UuIUKjKkzK7KcMc92jramP/z+XQe7gTgYO/BvkABeHn7ywBMrJhI685Wunu6OdJ7hMWXLOalbS/RtLWJF7e9yOQJk4cNLyluOqYiQ6qqrGLRvEUDjh7SB2UXPrOQzsOdjLOhf6UO9hyku6ebcTaOb9V8i7Vb19LR2cHpJ57OrTW36qBtGdBIRUYsc5yl81AnS9ctZfmG5UyeMBmAXu+lggow6PGeQdfR6718d813OdjTd2dR2va2adenDGikIiPWd5zFYNaUWbTvb+fVXa/2Te+hZ8hAyUgHCsCcqjk6Zb8MaKQiI5bZRflc9ed4eP3DUdZ5/ozzaeto45FXHgGGPygsxStfZU8nmtlTYfp6MzsjNe07ob3NzK6M0R8ZW5njLHevu5v93fuZVDGJyorKnNa58087adraxMwTZvJoy6M0v9ccqbeSbzmHSqrs6dXAuUCdmZ3bb7abgP3uPht4CLgvLHsuyd33/xq4CviPsD4pYpkDtN+/+PvMOmkW3T3ddPV0Db/gELZ/uL3vecv7W6h9sla7QSUqL2VPw/ufhtcrgcssua3+NcCT7n7Q3bcC7WF9UsQyx1TuXnc31SdVR133zBNmMu34aezp2qMrl0tUvsqe9s3j7keAD4BpWS4LqOxpMamfW8+C2QtY1b6Krfu38pHKj0RZb4VV8MPLf8gfFv5B1wSVsBihkk3Z08Hmybpkqrsvd/cad6+ZPn36CLsoMVVVVvHglQ8y7fhpbPtwG7u7dg+/UBZ6vIclLy456twY3cSp9OSr7GnfPGY2HjgJ2JflslKEGtsa2Xtgb9/74yuOH/W6DOOUylOoPqmahmuO3uXRTZxKT17Knob3N4TX1wJrPKli1ghcF74dmgWcBfw+Qp9kDHV0ddB5qJOLqy/uazvQc2DU63OcXV27+MT0TzB72uyjRie6PULpyTlUwjGSTNnTt4BfZsqemlltmO3HwDQzawduA+4My74B/JKkfvJ/A7e6Z3HWlBRUQ2sDS9ctZcL4CSz67CJOmHBCzuscb+NZ1b6KhtaGo0YnQ10mIMUpyslv7r4KWNWv7V9Tr7uBrw2y7D3APTH6IflRP7ee5zY/R9OWJnD4+NSPs3HXxlGvr8IqOOJHmDppKrVzaplWOa3vc6T06DR9GbGqyirmnTYPSO6bcuLEE3NaX+aU/v3d+2lsa9TopMTpNH0ZlYUXLuy7AdNFMy/i3affZXfXbrp7uke8rlMnn0rF+Aq+dNaXNDopAwoVGZWqyiqWfH4JAF/8xRfZ9uE2IPkmxwc+K+AY0yun86mPfirZjQJOn3K6RidlQKEiOXvwygc51HOIOVVz+OOHf+TptqcxMw72Hn0V8qSKSUeNZI70HmHZ1ctYsWkFmI6hlAuFiuRsTtUcVl+/mgdeeoBHXnmExRcvBoPV7av53f/+jikTp/D+wffp7unmuHHHcbj3MONsXN8xFN2btrzoQK1EUzunlgWzF9B1uIulLy7l8tmXc//8+2m+uZnFlyxm/qz5HO49zPTK6fR6LwtmL9DopAxppCLRNLY1sqp9Fed/7HwWX7IYnL678C+8YCE4zKueR90n6mhsa8z6Dv1SWjRSkWgyZ78uvDAJkKXrliZlOvjzCXOTj5vMnKo5+sq4jGmkItFkzi/p6Orgpe0vJY3hktFs78ovpU+hItE1tDbQtKWJBbMXJLs9/DlwpPwpVCS69KhEuzh/eRQqEp1GJX/ZdKBWRKJSqIhIVAoVEYlKoSIiUSlURCQqhYqIRJVTqJjZyWa22szeCc9TB5jnPDN72czeMLPXzOwbqWk/MbOtZrYxPM7LpT8iUni5jlTuBJ5397OA58P7/rqA6909U9r0381sSmr6Inc/LzxGf6NTESkKuYZKupzpT4Ev95/B3d9293fC6x3AbkDVwETKVK6hcoq77wQIz0PWvzSzC4AJwOZU8z1ht+ghM5s4xLIqeypSAoYNFTNrMrPXB3j0L8I+3HpmAD8H6t29NzR/BzgbOB84GbhjsOVV9lSkNAx77Y+7zx9smpntMrMZ7r4zhMaARXXN7ETgt8C/uHtzat07w8uDZtYAfHtEvReRopPr7k+6nOkNwG/6zxBKof4a+Jm7/6rftBnh2UiOx7yeY39EpMByDZV7gcvN7B3g8vAeM6sxs8fDPF8HLgZuHOCr41+Y2SZgE1AF/CDH/ohIgVlSJ7201NTUeEtLS6G7IVK2zGyDu9eMZlmdUSsiUSlURCQqhYqIRKVQEZGoFCoiEpVCRUSiUqiISFQKFRGJSqEiIlEpVEQkKoWKiESlUBGRqBQqIhKVQkVEolKoiEhUChURiUqhIiJRKVREJKoxL3sa5utJ3Z+2MdU+y8zWh+WfCjfJFpESlo+ypwAHUqVNa1Pt9wEPheX3Azfl2B8RKbAxL3s6mFCW41Jg5WiWF5HilK+yp5NCydJmM8sExzTgfXc/Et5vB04d7INU9lSkNAxbodDMmoCPDjDpeyP4nGp332FmZwJrQq2f/xtgvkHrhbj7cmA5JCU6RvDZIpJHeSl76u47wvMWM1sLzAX+C5hiZuPDaGUmsGMU2yAiRSQfZU+nmtnE8LoKmAe86UkVsxeAa4daXkRKSz7Knp4DtJjZqyQhcq+7vxmm3QHcZmbtJMdYfpxjf0SkwFT2VESOobKnIlI0FCoiEpVCRUSiUqiISFQKFRGJSqEiIlEpVEQkKoWKiESlUBGRqBQqIhKVQkVEolKoiEhUChURiUqhIiJRKVREJCqFiohEpVARkagUKiIS1ZiXPTWzL6RKnm40s+5M7R8z+4mZbU1NOy+X/ohI4Y152VN3fyFT8pSkImEX8FxqlkWpkqgbc+yPiBRYvsueXgs84+5dOX6uiBSpfJU9zbgOWNGv7R4ze83MHsrUBxKR0pWvsqeECoafBJ5NNX8H+CMwgaSk6R3AXYMsfwtwC0B1dfVIPlpE8igvZU+DrwO/dvfDqXXvDC8PmlkD8O0h+qFayiIlYMzLnqbU0W/XJwQRZmYkx2Nez7E/IlJg+Sh7ipmdAZwGvNhv+V+Y2SZgE1AF/CDH/ohIgQ27+zMUd98LXDZAewtwc+r9u8CpA8x3aS6fLyLFR2fUikhUChURiUqhIiJRKVREJCqFiohEpVARkagUKiISlUJFRKJSqIhIVAoVEYlKoSIiUSlURCQqhYqIRKVQEZGoFCoiEpVCRUSiUqiISFQKFRGJSqEiIlHlWkv5a2b2hpn1mlnNEPNdZWZtZtZuZnem2meZ2fpQi/kpM5uQS39EpPByHam8DvwdsG6wGcysAngEuBo4F6gzs3PD5PuAh0It5v3ATTn2R0QKLKdQcfe33L1tmNkuANrdfYu7HwKeBK4JtX4uBVaG+bKpxSwiRS6nEh1ZOhV4L/V+O3AhMA14392PpNqPKeORkS57SlLRsBwLj1UBHYXuxBgp120r1+2aM9oFc6ql7O5DVSTsW8UAbT5E+4DSZU/NrMXdBz2GU6rKdbugfLetnLdrtMvmVEs5S9tJqhNmzAR2kKT7FDMbH0YrmXYRKWH5+Er5FeCs8E3PBOA6oNHdHXgBuDbMN1wtZhEpAbl+pfwVM9sOfBb4rZk9G9o/ZmarAMIoZCHwLPAW8Et3fyOs4g7gNjNrJznG8uMsP3p5Lv0uYuW6XVC+26bt6seSAYOISBw6o1ZEolKoiEhUJREquV4OUKzM7GQzWx0uU1htZlMHma/HzDaGR2O++5mt4X7+ZjYxXI7RHi7POCP/vRydLLbtRjPbk/p3urkQ/RwJM3vCzHYPds6XJR4O2/yamX06qxW7e9E/gHNITsZZC9QMMk8FsBk4E5gAvAqcW+i+D7Nd9wN3htd3AvcNMt+fCt3XLLZl2J8/8A/Ao+H1dcBThe53xG27EVhW6L6OcLsuBj4NvD7I9AXAMyTnlF0ErM9mvSUxUvEcLgcY+97l5BqSyxOg9C9TyObnn97elcBl4XKNYleKv1vDcvd1wL4hZrkG+JknmknOK5sx3HpLIlSyNNDlAIOe9l8kTnH3nQDh+SODzDfJzFrMrNnMijV4svn5983jyakGH5CcSlDssv3d+mrYTVhpZqcNML3UjOpvKh/X/mRlDC8HKKihtmsEq6l29x1mdiawxsw2ufvmOD2MJpuff1H+G2Uhm34/Daxw94Nm9k2SEdmlY96zsTWqf6+iCRUfu8sBCmqo7TKzXWY2w913hmHl7kHWsSM8bzGztcBckn38YpLNzz8zz3YzGw+cxNDD72Ix7La5+97U2x+R3Naj1I3qb6qcdn8GvBygwH0aTiPJ5QkwyGUKZjbVzCaG11XAPODNvPUwe9n8/NPbey2wxsMRwSI37Lb1O9ZQS3L2eKlrBK4P3wJdBHyQ2V0fUqGPQGd5lPorJKl5ENgFPBvaPwas6ne0+m2S/8W/V+h+Z7Fd04DngXfC88mhvQZ4PLz+W2ATyTcOm4CbCt3vIbbnmJ8/cBdQG15PAn4FtAO/B84sdJ8jbtu/AW+Ef6cXgLML3ecstmkFsBM4HP6+bgK+CXwzTDeSG6xtDr97A37z2v+h0/RFJKpy2v0RkSKgUBGRqBQqIhKVQkVEolKoiEhUChURiUqhIiJR/T9WQN0hLQ0y6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX+QHcV94D8Nkoh3nYuFVnFAsQxY3O7ZUPHaAgwcLjm3Z8TmstgJrkO2Y2XPLuIkou6iu9XZ8ZGVQlVi71bJuWhTcTjHayh8MolTXJSwePFaEghlQci3tgVh1ywSCIGD9wlMivdylmT6/ngzT7Oz86Nnpmem56k/Va/eezM93T093d/59vfbP4SUEovFYtHFeWVnwGKxtBdWqFgsFq1YoWKxWLRihYrFYtGKFSoWi0UrVqhYLBataBEqQoivCCF+JIR4MuS8EEL8qRBiXgjxfSHEezznNgshnnE+m3Xkx2KxlIcuTeWrwMaI8zcBlzuf24A/BxBCXAgMA9cAVwPDQoiVmvJksVhKQItQkVI+ArwSEeRm4B7Z5DHgLUKIi4AbgW9JKV+RUr4KfIto4WSxWAxnWUHprAFe8Pw/4RwLO74EIcRtNLUcOjs739vT05NPTi3aOfPGGWqNGl0dXSw7r6gqZ8nCd77znZqUcnWaa4t6wiLgmIw4vvSglHcBdwGsX79eHj58WF/uLLmy7aFtjE6P8tFrP8rIB0fKzo5FASHE82mvLcr7cwJ4m+f/LwIvRRy3GEatUWP04Ci1Ri3ROYCZl2cWfVvam6KEyh7gE44X6H3Aa1LKHwKTwAeFECsdA+0HnWMWwxifGWfb1DbGZ8YTnQO4c8Od9Kzq4c4NdwYKoDihZKkWWro/QojdwAagSwhxgqZHZzmAlPJLwATQD8wDDWDQOfeKEOJO4Aknqj+UUkYZfC0FU2vUGJ8ZZ6B7AIDB3sElYdxjQecADhw/wOzJWQ4cP8CB4wfYNrUNgKHrh4CzQsl7zFJdtAgVKeWmmPMS+N2Qc18BvqIjHxb9qDT4ro6uJedcYTTYO8hA9wD7n9vPQPcAqzpWATDQPcDowVEGewdjhZKlWlhTvCWStA3eK4wAJuYn2HDJBoauH2Lo+iFGD44uElZWQ2kfrFCxRBKkhajgFUYnGydbmkrQeUt7Yef+VBzTjJxufqCpgXR1dLH7yG4m5ifYfWR3K5wrrLo6usrKqiUnrFCpOH7PS9lCJsgTdOKfTwAw/8p8YXkruxzOZWz3p+L4uxFleFK8Rtmgbs3Dzz8MwDef/SZfe/Jr2vPmTd/VfHSUQ1C8FgWklJX7vPe975WWYBbqC3Lk0RG5UF8oLM2RR0ck25Ejj44Enp8+Pi17dvXIB3/wYKK8qd5LUPo6yiHuvtoZ4LBM2T6tptJmpDWsZsGvnXjHtuyZ28MNa2/gspWXcenKS9l4ufp8UVVtI0g70lEO1picDiEruEWHnftjNq67uH9dPxPzE/Ss6mH25Cz96/p54GMPKMdjux/lIYT4jpRyfZprrabS5uTRMOPidN/sA90DbLhkA1f+/JX83uTvccf770iUThlalyU71vvT5sTNy4kizIPijTMojCsMuru6Gbp+iMdefIzZk7N889lvZr6ftGTxBllPUjKsptLmZLELhNk0vHEq2T2k7zshOrStLN4gOzcpIWktvGV+zkXvT5w3Iw+vj0qcusJEocO7kyUPZXjUyoYM3p/SBUSaz7koVOLcmya6P93GOLswm6lRBjVqE++3ncgiVGz3pyLEdWPycn9m6Xq43Yb7n76f6Renef615xnrH0uchyCDbZnuXuuVisYaaitC3FyZqPNJDY3e8KqG3qA0BnsHGekb4eX6ywBMzutZf6vsRp3F+H0uYDUVA8i7kSQ1NHrDq2oEUWn80S//Ebd/83Z23bQrNm2VsijbcGoHxcWQtt9U5qfdbCp52weKMGpG2T367+1Xvj+VstBlq7GEgzXUmk0Znpsi0wmLN03jT5JHa6zNj9KFCs0NwOZorkH7mYDzXwS+63x+APzYc+6nnnN7VNKrmlAxpfLnlY+4eHVpPrrCZBWu54KLuVShApwPPAtcBqwAvge8MyL87cBXPP9fT5pm1YSKKZWwaE3FJY0wU7lmdmFW9t/bL2cXZhPFk1W4eq835dnqpmyhci0w6fn/WeCzEeH/Afj3nv9tL1RUaNfKKaVa4/ejUh6urab/3v5E8ejUVEzRQnVTtlC5Bfiy5/9vAGMhYd8O/BA433PsDHAYeAz4kEqa7ShU8qycZdtshvcNS7Yjh/cNa81jGmGlG1PsZbrJIlR0jFNR3roUuBX4hpTyp55ja2VzivVHgT8RQrwjMBEhbhNCHBZCHF5YWMiWYwNxx3Tk4aZ0XbCb79+c66S4sPEbjdONRd9Jro2iu6ubBz72AN1d3eky7CHtWB4gcvzQuTimRcc4lSRbl96Kb/8fKeVLzvdRIcR+oJemjQZfuEV7KWfOtWHkOc1/sHeQ/c/tZ2J+gvGZ8VzT8X4Xda0KceNfsozlSbqAVNuTVsVxPzQF01HgUs4aat8VEK4beA5nYSjn2ErgAud3F/AMEUZe99OO3Z+8KUMNX6gvyOG9w/LaL18r2Y4cmhwqLG0/uj1UVe3WqIIBLuV+mq7iZ4HPOcf+EBjwhNkOfN533XXAEUcQHQE+qZKeFSpL0VHJdTcUtyG7nyibSt6YKARMzJNL6UKl6E9VhEqRlUaHoTfJaFaVsSGzC7NyeO+wHJocksN7h+VCfcHIhqQ7T1kW7DYFK1QMpchKU5SmonJPUWGSlElRAkj3c1KNz0QB62KFiqGYujBQ3vmKCmPiMHzd5WGysFDFCpU2RFdXJG3cutLKcq33ujIbalTaJndhsmCFShuStCuSpNGlaaBBA9iKtB2U2Xij0jbRQK4DK1QqQB4VxzuiNO9GN7zXESp7h1vHVNPUMfLVVE1FByZqO1aolIwuA2dSipzYlsV2UFajMVEDCMLEfFqhUjJ52j+iMLEyBqF7foxqeJU5Rzoo8tkW9cyzCBW7Rq0GVObtxK0xm4Y84syDuHwmnR+jHD5kvyHdm4PlvWGbrrQKI600KvNjmqZiOiZoNLrczEnCh4XT0R3T5ZkKy0uVNZXSBUSajxUqySjbEDh9fFquHlltjDFSh4tbV5maIPCDsEIlR/Ic3FRUhSq74vbs6pFsR64eWd3KQ9l5SkORhvGysUIlR4LeSO6x/nv7M1WqsjWIopg+Pi17dvXI6ePTrWNVvPd2FyRerFDJkTBNJcm2E0niLuLaomnXoeztjBUqJZBXozB17EcWQ2tRrt12pCzhm0WoWJdySvJy56q6DFXc2Dpdp1H5isuzynKSlmCKWgpUJ3bbU8NQXX5QZflJnduDRuVrsHeQ+uk69VN1ao3aEkHbsaxj0bdFHf9SoIO9g+ZvDp9WxSnzU5ZNpWp2gCLzm/ekO5PQeT8q41HycGXHgbWp5If3IZrmsdA5oCxr2lHpBXl/qozOehAWl8qguDyfsRUqGkj6xjCBoIrn5tE1juYlAJM0LHecSs+unlzyUjRFayph5PmSK12oEL+X8m8CC5zdM/lTnnObaa6i/wywWSW9PISKaVqICkEVz72P4b3Dhc9aDqPdNBVTaFtNBYW9lB2hsmTXQuBCmtt7XEhzu46jwMq4NHULlYX6ghzeNyyHHjq7QHNVybOiFaGpqaRRtsZYdvpFkEWo6HApXw3MSymPSilPAV8Hbla89kbgW1LKV6SUrwLfoqn1FMr4zDg7Ht7BUz96ih2P7Ih05+qe4aqbOFd3lvynnSGbJE03jbFDY6HXlD1Tt+z0ozChfupwKa8BXvD8PwFcExDu14UQ76e5P9DvSSlfCLl2TVAiQojbgNsA1q5dqyHbZ3HdpAPdA2y4ZEOkO1enm7YMsuTfW06jB0eV3ZpJ0nTTqJ+qh15T9q5/ZacfhRH1M62K436Aj7B0g/ZdvjCrOLsT4aeBvc7vIeB/eMLdAfzXuDTTdH90qawmbAqeBR3lkNT+lMRTpDOfZbNQb+7QOLyvuC61rnKjZJvKtcCk5/9ngc9GhD8feM35vQn4C8+5vwA2xaWZRqgknQRo6pwfE8ia/7RGcTddd03eIssvzT2791k1B4CU5QuV2L2UgYs8vz8MPOb8vhA4RtNIu9L5fWFcmmk1lSQCIa/ZybMLsy0Xa1g+klTgKgqotNqet/yLbqhpBGGQplKV51WqUGmmH72XMvDHwFOOwNkH9Hiu/U80XdHzwKBKemm9P0nU7qC3oY4K4TaInl09ofEkqcBVdIXH5Tlu7EYWTSXuGWYZN6JC2ORK04RN6UKl6I8Ol3JYP9996FnXSglLb/r4dOxbughNpchKnNSm0hprs0//WJs4gaZDSEfdr3+rk1a9c46rpFvEs7NCJQX+yuMdNJaHem2aRpElP2mFRFLDbpKGpprfuK5XUm0oagCiv2654bwjtN26lkSAFlGXrFBJQVTDyOLhyVt91kWW/EQ1miC85VmGvcibP9UGmSWcqtBNa6OzmoqhQiWKLC7TrJ4NU4ROFEk1Fa8dIa+3bFQe/NqBitDPGl/SfKmELxIrVDST9MGGqbdJ4jWte6QTrx1BV6Pxx5PVM5eXcTzpczWlHlihEkIeUj+oa6SaTlSF0ZlXU952Lnnkx1+WXvuEDjtRGWGDwqter7uMrVAJIQ+p71bc/nv7E19bVGMvo4tRNEV0HdIID90D83TbgFSxQiWErJUs6K0xNDkk++7uM3qYfl6N3xTVvCjSdHN0ew7TaipZ64AVKjnhH6h0rjUqP2kraprrTNCKTNBUVAmzMaWtq1aoaMT7cMIGKpmg/lcJ1Qquw4tWFGXWhaC0g2xMVlMxRKioenIs6qiO+8lS9kUbNMsUeipjY7JihYpGrCDRj+pmYlnKPqtBM60Qy9rdybtrmLZMrVBpQ1QqQ1UE4NDkkGQ7cmhyKLc8Z9VU0moeWYcJZB1fE5de2vuyQiWGtAPSyiSuMmQdl1EkXqFSRLehSMNw1gGN/ueYRTjq7BZZoRJD1MPVVcnjHp5uG4HON1zaPKhShKbixRQjb1IB4dqe0nbjdK5KaIVKDFE+/DwMdzrV0DCq1DiL3qDdVO0zjiQvCquplCxU/KRpLEk0kSKs82Hx6UxHV1x23x81wsrb1UAe/MGDSss2JHmZhaVphUpC3DEo7jJ/SQxqKoKoCCNrWH5MUf299N3TJ9mO7Lunr+ysBGKCZhOVB7dLtHpkdeIpIv667iesvlih4ilA1crhLUyVgUO6K17Wxl+EpqILr00lDXnfU5bxMXnkwY+qppImbmM1FeK3Pd0K/CPwfeDbwNs9537K2e1Q96ikFyZUknQ7ouwqeb7t/eMbdA7rNlGgSJmfVqaLuK5rWNgk55LkQTdp4i57NX2VbU8/AHQ4v38buM9z7vWkaQYJlTA1L42XJM8H7K+0Ot+SJnZ9vOThttV9varXLU9PYhHE3WfZQiXpvj+9wEHPfy1CJeyBqmyJUSRRmkpeXSL3XNEbW/nz5XqBin4OWco1ynMYF9YUgrTyuPV/yxYqt7B0h8Ilm7F7zo+xeFfCM8Bh4DHgQxHX3eaEO7x27drIgvOSRlMpgjQeoiyV1k0v7TKZWcK4aQ/vTbY6vur9pik31biL8ORFoSMtbxvwL+1pqqYSu+2p59zHHeFxgefYxc73ZcBzwDvi0kzi/anC20OVrG/cIE0lKh8q6amESetSVr3fNJ45Va1JxU1bZHc5DQt1z6r9ioK9bKGi1P0B+oCngZ+PiOurwC1xaZ4Lc3+C0F15vZUtiWcgaZi0q+Wp3m+SkaRJtCZVD1ue3iNd8SWNp2yhsoz4bU97HWPu5b7jKzm7cXsX8IzfyBv0ySJUTNVcykB31zCsbPMe/JZGU1GxiajGm8R7pIu867EJLuW4bU+ngJf9rmPgOuCII4iOAJ9USS+LUEn70E0XRmnyp/uewspWlwqft4FUR7emqHqStEwrpamU8SlDUzHV4OtigjtTtbuQhiLuz/QXh5ekeU1aflaoFECc/aFsim4QQW7KPLxWaeNIkr+qkVaDGt433Np7KQ4rVDSgMh6hqMpZhUbgffPFvQXL0KKS5K9qpL2fJNdZoaIBf4H7F71OShLBkNZAmDQdXfjd00VoKmny2A6aipt37yDJtPeT5DlZoaIBf4FnXQMkiWDIYiAsWwuw6CXsBVNE19v7XK1QcdD5RsoaV9jYiaA3aJZJhVnz6dc6ikjTEk7YC6aIiadWUwnApDeoins1z/yqNnw3DyorzJdhcG0HkpRjEeWjUu+yCJVltBGDvYOLvtNSa9QYnxlnsHeQro6uVOHC8hJ0PGt+gxg7NMaOh3dQP11n+4btoXke7B2kfqoOYmk+xmfG2Ta1jfrpOp3LO6mfrrPj4R0ADF0/pJwXb5punEnjKBvVOhGE/56jyqCro6t1rNaoMXZoDCRsuWZL4nTD0NVOQkkrjcr8/FLvL+WqCqpqECZpRn7CDM1pjMDeuNKUu9c+VVVNJcuzTqvxRWmRecO51v1Z071GW0EHVRbVh64656QMV7TOgWhZ85/Vk2YCpXnZEowt0ck5J1Ty1lRUyUujSTPQa3jvsOy7u89IzSnPBmmS8GwnzjmhksalnMdgNtXrk6adVAh51WTvNIIsY2WqgopBXIWFurkjpst4NlaoKOD196cZZKabqLEpaTWVsKU004yVyYM8hLguTSWofphC0fV0ob4g+VlekFaoxBeU/01UZp81bJCTzoqj0rDyGAcRRtZ7NMEFXwZF523k0RHJRUhphUo8WTw9eVFkow6iyPvPaucwueGXjc6ysZpKQIGovp1NqKi6PCNp76Ps+4/DLZ++u/uMzWNZeJ+dysshybPOYlM5L5/RL+XhDiwanxlXCuMONtI1sCgxwvedEpX7DkL3/dcaNUYPjlJr1JSOx+KUy9SxqcT31u54n/lg7yAjfSORA9rS1pHEpJVGZX50aipZwuhAVzpZ4tGp5QS9MbNsk1KU3ct0jS2INAb9IjSV0gVEmk8REwrLtrWoEHW/ScoiyhMVRVDXLeha10Des6unlEYb5h3zUoXnXSSlCxXitz29ALjPOf84cInn3Ged43PAjSrp6ZhQGNdwdI7x8J/XJfyi7jeJrSatJ0p1eYgHf/CgXD2yWj74gwdj85IE1XL0juMJGzldRU0lT0oVKqhte/o7wJec37fibHsKvNMJfwHN1fifBc6PS1P3Bu1ZiYvLf15X2lH3m2U9GNVyVA3njvTtu7svcV6iUC1HHeN4wuItSxDlnXYWoaJjlvLVwLyU8iiAEOLrwM00N2R3uRnY7vz+BjAmhBDO8a9LKX8CHBNCzDvxTUcleOaNM4weHGWge4DdR3bTON0AAY1TDeZemWOge4D5k/Pc+je3IqXkvlvu431ve1/r+lqjRv10naFrh6ifrreMh95ZqHO1ObY8uIXuC7vpWN5Bx/IONl25iT1zexjoHmDP3J5W2IHuAR569iEWGgvUGrVFRs9ao8ZCfYG+S/u4Ye0NjB4c5Ya1N9B3WR/P//h5tu/fzqYrNrXi8+bDnyc/3hmtblpu+E1XbOKJF59g4zs2tsrKzfvuI7uZf2WeifkJfq3n1/jUez7FHfvuoPcXetn2b7e1Zi8HpT1Xm2Pr5FZ23riT7q7u0HBeGj9pLPr25nfk0RFmXp5h7KYxuru6l8wGrjVqjBwc4dCJQ1zx1ivo6uhiy9VbWnmE5mzboOvGDo3RONVoGXs7lncsqQPD7x9edP1A9wC7Du3iwWce5KZ1N7H9A9tb9+UNs/vJ3Rw8fpCpY1M8dPQhrv/F6xfNJK41aow9PgYCNl2xid1P7gYJm648+3vLNVs42TjJb/3dbyGRfKHvCxw4foAb1t7AnY/cyR3vv4MDxw8EPoOb7r2JY68dY6G+wMgHR5aU6aPHH+WMPMOblr2Ju371LgBu+7vb+Jcz/8Kblr2Jj1/5cT6373NsuWoLY0+Mcc+H7mH9mvXN5/FPM7CMCwIfpgI6hMoa4AXP/xPANWFhpJRnhBCvAauc44/5rl0TlIgQ4jaaW5/ylre+hW1T29j/3H4m5ieWhN06uZWjrx7l+deeB2Dwbwd5esvTrfPjM+PseHgH/ev6mZifoHN5J8Ci6ehbJ7cydXSKqaNTreueeOkJJuYneOjZh5g6NkX9VJ3tH9jOnrk9TB2bYurYFKs7Vi9qaOMz44xOjwKw4vwVTMxP0L+uf1HcT7z4xKL7cPPhz1Mc3in1QCvOifmJVln5y+wvv/uXHHzhILMnZ5v572zmf/P9m1vhvGlvndzaOv7Axx5QWspg+qXpRd/e/Lpls3Vya2B84zPjjP5DM8zDxx8GoHN5Z8tj5aY5enB0yXXuMg1e3Gvd8yN9I3R1dLWu95bPnx3+M97+lre30nDz5g3Ts6qn9Sw7V3QuCrvjkWb63ufr1iGAzhWd7H9uf+u+Bv92kNmTs/Ss6mH25CxHXz3K7MnZwGdw7LVjAMy8PBNapt7wAI8cf+TsMzkxzek3TrP94e1IJJ/4P59g6Lqhs9f+HG9bUniK6BAqQc5QqRhG5drmQSnvAu4CePd73i1/v+/3Gege4KqLr1qiqey8cSfHXj3Gpr/ZxJuXv5nxmxe70Nw33ED3ABsu2RC4tsnOG3dy6o1TSzSVDZdsYKG+wNSxKRqnGy0toH66DpIla4YM9g62zrnXD3QPcNWaq2icatCxooON79jYys+qjlWL8lE/VW9pU3Fu36B1Mtx79H5fdfFVoZrKQPdAS6D0XdpH/dTitHfeuHPRt8raHEPXDTHyDyMMXbdY6Az2DrJQX2Dm5ZnQ+AZ7B1loLCzSVILSCrqufrq+RFMJCuf9Huge4NKVl7Y0laC6MdA9wLtWv4vHX3ycK996Jb/a/at0LOtYErZ+qt6sm8C7Vr+rVYeuWnNVq64MdA9QP1WP1VS87LxxJ6+ffh0hBWM3jQWWqVdTccv29VOvx2oqC/UFZv5phqnXpl4gLWn7Te4HhW1PgUngWuf3MqBG81EvCusNF/VZ071G2TiX55Bu1xjqTyPtvJOgvBbtlfDOgVHdbziOtAP8yrBZ+I23YV6jMOOvHx3PL69yiIqXkm0qTwCXCyEuBV6kaYj9qC/MHmAzTVvJLcBeKaUUQuwB/rcQYidwMXA5cCguwRf/+UXGZ8YjuwR5rm7lqt21Ro3OFZ1L0vDbOuKIymvuq3Q5eO0F3vQ6ly+9v7BrQ+0qCQb45blKnMrqbf4upNuFcbtNLlEr5nnT8pdnGvJYLa/WqIV2cTOTVhp5P8Rve/ozwF/TdB0fAi7zXPs557o54CaV9FQ0lSqi2+uShCxv1Lhr03rldN9n0qHsKuNbsqSlSp7PO2xWNmWPUyn6U8YOhV7yUkdVK2IeXaI8R+SqrpCXNR8mxW36uJe4/Fmhopm4As/LzlGmppIn7oja/nv7Y8MG3VsV7tfkEblpys8KFc3oVOctUk4fn5Y9u3rk9PHp2LBBZV+2EC86Lt3ppSm/c1qo5PEw08Z5LgubqHvPOn2i7O6mToqYoqEjzXNaqKStGKYZO/OiKEEXde+mCtsy8lXEFA0dnNNCJU3hLtSjFzluJ02lKEFXxliKMvOVFtPyE8Y5LVTSEOdOy6r9lLU0ZFXylIQkkwaD7q/slfWqihUqCYmrIGkrkFdYmdYVMLFrpoJqeYXdX5bZ2irxJ8mjalj/WJm860tQGuekUDHxzaGiFSTxLBXlncg7nSJc5WHX6ro3XYZolbDeMEW8DILSOCeFSlXdvknGwBSlXehMJ4tLOI/7LetNnyWs1VQM1VTaQd1PMhJVV5p5GEaL0FTCyNs7mLeGVBbnpFCJowpeg4V69NySMgRjHmkWJRyDiCvjMLJqV1V9qblYoZIDRYy7cNMwyaibh+0lr+UkVYny9mXVNIrUVIqsD1ao5IAuw1xcGmlnwZZB2vsemhySbEcOTQ7llLNoFurh45KqpFEE5TUvQWOFikZU3M15C4IiPBZFxmeCfSEvjaRIgvKQl1C0QkUjcQ9J90PMs6KY8hYu06aiC1PK0o+JmoqOld/airiV1nSvxBa0qpeuNIpaNS4O/2LZVULnCm55kHSVwUJIK43K/JS9SJNOTFCr8ybJ0gc60Nl9jBod3c5gN2gvn7QbkOvaID31BuiaiEr/wPEDzJ6c5cDxA4XkRddG5OMz460tVUzTUIwmrTRqCjMuBL4FPON8rwwI826aC14/BXwf+I+ec18FjgHfdT7vVknXRE0la587y6AxHelnpQgXvGq8/uNVNjCXBWUZaoERnL2Tgc8AXwgI86+By53fFwM/BN4izwqVW5Kma6JQWagvyOF9w3J4bzqvUJbh7W76ZTaAMtL3l0+YZ65IV2wYZT+fpJQpVOaAi5zfFwFzCtd8zyNkChMqRTzULNpCGk3F1Ddw0vh13Yd3MKF3VnKRrtgw4tIzTeiUKVR+7Pv/akz4q4GngfPkWaEy53SLvghcEHHtbcBh4PDatWsTF1IRlajoipH2nrzXFTHfJi4NnYMJ3dG7ceunmKaplN199ZOrUAGmgCcDPjcnESquJgO8z3dMABcAdwN/oJLpvDQV094WcSQZuOXtFnivy6Myh2kQYUJGZ7kX8QzzyG8ei2hlyafx3R/gXwH/F/hIRFwbgL9XSTcvm0pYA8uzm1Fkt0znvSXBP/jNtLeylOUZxfMsiyxxlylURn2G2pGAMCuAbwP/JeCcK5AE8CfA51XSzUuohFWspOp82HVpw2SliKkFUaQtvzyJ06aSXKszHzqpqqayyhEYzzjfFzrH1wNfdn5/HDjtcRu3XMfAXuCI0526F3izSrpBQqXIh6NaAVVcniY0MB1E3YeJ92iioDOJ0oRKWZ8goVKkSp21O+Sum5o0r6Z5bbxUofx1x1EEZeXTChVZjUriNrzhvcOp8pp3w9XtEs+LVjnuG16i8VV114AwyrI/ZREqbTOhUOfEqlqjxtjjYyBgy9VbMg+hd/FO8Esa51xtjoeefYiha4dyGzKeZQJiVPnP1ebYOrmVnTfupLurO1Me4Wz+6qfqrcmYANumtrH/uf2tyYth+XEnCaZ5DkWTdVJoKfeaVhqOgQ4oAAAQSklEQVSV+SlqL+Wkb4g839ZJNjk3jbzyHmSbUtFUdL/9y+iOJXUWJNWOsZqKXgZ7B6mfqtM43aB+qk6tUVOS8kHLGOhi5407F31Xibzy7teO3N9DXWePed/UQC7LGOh47ipxeO8lLrx/yYb66XpudXMJaaVRmZ+kmkraN0nSN1oV7DploOoZymNcj/cZ5mWfKEpTSTISOqt3C2uojSZtZTqXhESe9xpV/l6jq8raJVkEfdWfZ5IuXtZ7tUIlBlNdsSZV8jy9DFGLNLXc7M6ex2H7W/vDm1BmZVGER8gKlRzwV96oUalxDznsbWnScPUsjTXuWpUtOsoe9RtGUD0oW6ilWfO3yO6PNdSG4DeEjc+Ms+ORHQB0Lu9cZOyKc/t54wJav5O4C/N2DWZxyccZDXsv6mXq2BS9F/WG3kdXRxedKzrZNrVtSfmWRa1RY/P9mxe5qPM0xrtpxj3nPXN7mJifYMMlGxYZpaPiyjvfi0grjcr8mKKpqEp/77VpB2eZpNX4iSsHb/cn6j5MGynr5tXbJctbU1F5zkndySqGXT/Y7k85JGnoWY3FVR4p6u3+5OHh8aJT+JbR1dGZZpa4rFApiSSehaLc2lnTy4OkOxTqmi5gqoG+ClihkoC0jT/puICsecmqoeTdXUrSoJIaFnVpM3mXgcld0qxYoZKAuIoQdj6thycqbNRs5awV1oS+f5qwOtNPIszSlJfVVKxQkVJm11R0zC1xzw9NDoVW+qQVtugKXuRArKRx+se+lCX4ykKHG9wKlQLRYZ1X0VTyyFceJE1Xp+0pLKybJ+/SCHHXt5PW4X8mQV6sOKxQKRDd1vksewV50bkJuo6GHRYmbfczSdg0Ltd2IkhTSbp9qxUqOZPnW0xXxdbZQHQ3tiTjJXQLtDyvz4O88lSZcSoobHvqhPspZ9en3eM5finwuHP9fcAKlXR1CJUkhZznG01XJdKpyuuu2KrxFW17MRFTtKcyhUrstqfOuddDjv8VcKvz+0vAb6ukq0OoJHl4Ot6uWSp30mtNqZhJKctL5KJDAGV9zrq6w1kpU6io7vuzRKjQ3JajBixz/l8LTKqkm7emkkcjztJgdBtDTSUP4axbI42Lr8jnnCdlChWlHQqBMzS3LH0M+JBzrAuY94R5G/BkRFqZtj1NQh6NuEhNpej4kqaXJf2sxtkk+UySvo58mPQyyFWooGHbU+Bi5/sy4DngHcDqAKFyRCXTeRtqgx5u2aqxToq2EfmH6WdJP86NnPdb3pRnmDfGd39813wVuKXs7k9SdFRaU9Tbor1Z/vVUyjbItuv4FJ2UKVRUtj1dCVzg/O5yPD3vdP7/tc9Q+zsq6ZYhVMrSVKpW6YPyq3MMjQ68gi+roK/a81GlTKGisu3pdTS3Nv2e8/1Jz/WXAYeAeUfAXKCSrimzlIvAFO2mndCpqbTr88kiVETz+mqxfv16efjwYW3xmby5lMl5s7Tv8xFCfEdKuT7NtefpzkwVcZfaG58ZLzsrS3CXeWynCmsCtUaN0YOj1Bq1TPEU9Xx05bcIrFChuUbsSN9I4FqxSR5mXNgsFaNKlaoKmPwiCcLN7+b7N5tfB9L2m8r8FGVTSToRS+dkuaTXtqvBMC9MmWOTJN6kkwKzgJ1QmA9uQ1Zd+DjrcP4scbt5jZvun5Q8hVU7CsI8vUlFllcWoWK36IjA3VMZAScbJxmfGad+qt7aqsO/1UHcNhdx56O2UYi71u261U/p3TM3z60dCt02QhNxhtkk264ExZulfhlDWmlU5qdIl7JXW8lDE/Bi2qhd/wQ3ncPto67XMYEzL/JyIYdpmrrLXBVs9yc/3IeYZgHqKqv3QX14f4PKu4HlYZvKwkI9v10UVacfFHXvVqhoRtebUsXOYarg8WporkD1C9a83qI6yj+Pci1DmFlNpU2EStxMWO/asiqGtagFmE0dkem9L9U86rqXJBqBrgmGWQWVqS+HtFihopm4ijq8dzhRg6t6ZVTNY1obiR+3TLMIsrRpphWIpr4c0mKFSkEEVdQqCIWyCWpwcYI2q6aSlLyMzkXHoQsrVAyhiEphUsVTJSjPRbzZTS6rsspElSxCxY5TUURl4lgR4y6qNraj1qixfd92Jp+d5Ia1N7TKLu14jiSYWFat8Sin6+x4ePF4lCLKpBDSSqMyP6Yu0lSECmzy2zcIr32kZ1dPoWmbWFZBdjkTwXZ/8qeoCmqSChyFtzzi7CO/+/e/K9f9z3Vy+vh0CTk1i6iyml2YlX339MmhyaHShY0VKm1EFuGl2tB14BV+pgtCEzUWPwv1Bdmzq0fZ65U3WYSKtalUhKQ2HSCzPSEqzaD+f5gtYK42x9bJrey8cSfdXd2p8pIFE20rfsZnxpk9Ocu6C9fx4e4PV9uuklYalfmpiqaSZjRk2Fs/qU1Hx9vZn2baOP0LX6vkPytpy6IsrcY0bQqTtz0FPsDZLU+/C/w/zu7981XgmOfcu1XSrYpQSTNvQ9cAsqwEjRVR7eb48+rfoiMMnd2otHFFXWfKmJgiKFOoKG176gl/IfAK0CHPCpVbkqZbBaGyUF+6haUplSmtxqSa/7Qajk5NK+31UdfpEnqm26CkLFeoJNr3h+Yug1/z/G9boWJyxcnbPe7fkiNNXCaWn9VUihEqSjsUes7vBf6D5/9XHcH0feCLtNEWHUVWnKRzi/L2EvkFgjuhcnjvsHIcVWh47UyuQgUN257Ks5rMArDcd0wAFwB3A38QcX1heymXQZZGFPVWj3vj56ERLLGpPOTYVB6KtqlUlXYUgJXo/gD/Gbgr4vwG4O9V0q2CphKHvyJmadxJNZUk59Ok6cddKmJ433CiNKqC++y8axlXnTKFSuy2p56wjwEf8B1zBZIA/gT4vEq6VREqSYx+pr7twvKVRAhOH5+WPbt6Kj+iNsozV+RK90VQplCJ3fbU+X8J8CJwnu/6vTS3Qn0SuBd4s0q6VREqRbgn8ybsHpLk321w/ff2p8pDXuNXktIOz1OV0oRKWZ+qCJUqVrQ0A/biyLpBuwnjV6Ss5vNMixUqbUCasRx5UIThNu646vU682JZjBUqbYCqMTPv8Rt5NLowl7KJY1FUaXfhlEWo2AmFpiB93yFkWchHZVKi6oZVKnG5NE43Fn27VHlRoipMUiwLK1QMYcs1W+hc0RnbwLyNPknDBr0NQUdcldlxz0etUaN+us7w+4crKRDz5ryyM2Bp4jYwFeHg4jbs8ZlxpfCDvYOM9I1oaQhJ4upY0bHoWwe1Ro3Rg6PUGrVU57MwPjPOjod30LmiM9HzOlewmkqFSdp90KkZJIlry9Vb6Fwer4UlwRWo9dP1VtzeBp5n96TK3bYisEKlwqg07KRdpLzj0UXchvR5NvyqdtsKI62Ft8xPO3p/0hLnhciyhEFcPKqkmVAoZXnbm1qs9+ecxq/m+zWKoDd2mq5Bljd/mPcnDpV8Wq3BPKxQqTj+xq7SENMIiCyNt2N5x6JvVUyxXZjW9TOetCpOmR/b/QlH5+xnXagO0ze1K2NCGRYNtvtjcfFrFFFv+1qjxtihMZCw6cpN7Jnbk8vbePeTu5mYn+CqNVexfcP20HAqWlYZWoMpGlNVsEKlDfE3vLAG6o63AHjipSeYmJ8Achghqmm0cK1RY/P9m/PLZwjWbpMMK1TaEFVD7GDvIPXT9ZamsuGSDbm8jTeu28h9T93HxnUbI8PFNd7xmXEm5ifoX9dvtQaDEc3uU7VYv369PHz4cNnZMBbTDIu/8rVfaQmDBz72QOp4TLuvdkYI8R0p5fo019ph+m1ImiH/eXL71bezumM1t199e6bh86bdlyUYK1QsLfKaL7Pr0C4WGgvsOrQr8XylPPNlyQdrU7G0yDpfJqx7svPGna3vVR2rgGSelKB82a6QuVihYmmR1XUaJpS6u7oX2VKSCixdo4ItBZF2gItj4P0I8BTwBrA+ItxGmtt5zOOsvu8cvxR4nObC2fcBK1TStYPf8ifNQLQi94E2daBcu0CGwW9ZbSpPAr8GPBIWQAhxPvBnwE3AO4FNQoh3Oqe/AHxRSnk58CrwyYz5sWgije3DnWs0PjO+yP7hxjV2aEybbcQabc0lU/dHSvk0gBAiKtjVwLyU8qgT9uvAzUKIp4FfBj7qhLsb2A78eZY8WfSQtisU1C2JW6bA0l5oGacihNgP/Dcp5ZLBI0KIW4CNUspPOf9/A7iGpgB5TEq5zjn+NuBBKeUVIWncRnPrU4AraGpJ7UYXUG0Xx3kso5NV1DnJG5zxnOniPH4ccq7KVP+ZBdMtpfzZNBfGaipCiCngFwJOfU5K+bcKaQSpMTLieCBSyruAu5w8HZYpB+aYTLveFzj39tP2u7d2fWZCiNSjS2OFipSyL23kDieAt3n+/yLwEk3p/hYhxDIp5RnPcYvFUmGKGPz2BHC5EOJSIcQK4FZgj2Nh3gfc4oTbDKhoPhaLxWAyCRUhxIeFECeAa4EHhBCTzvGLhRATAI4WsgWYBJ4G/kpK+ZQTxX8Htgoh5mnuy/yXiknflSXfBtOu9wXte2/2vnxUckKhxWIxFzv3x2KxaMUKFYvFopVKCBUhxEeEEE8JId4QQoS674QQG4UQc0KIeSHEZ4rMYxqEEBcKIb4lhHjG+V4ZEu6nQojvOp89RedTlbjyF0JcIIS4zzn/uBDikuJzmQ6Fe/tNIcSC5zl9qox8JkEI8RUhxI+EEIFjvkSTP3Xu+ftCiPcoRZx2fH+RH+DfAN3AfkLmGAHnA88ClwErgO8B7yw77zH3NYIzFwr4DPCFkHCvl51XhXuJLX/gd4AvOb9vBe4rO98a7+03gbGy85rwvt4PvAd4MuR8P/AgzTFl7wMeV4m3EpqKlPJpKeVcTLDWdAAp5Sng68DN+ecuEzfTnJ6A8/2hEvOSFZXy997vN4B/J2LmeBhCFetWLFLKR4BXIoLcDNwjmzxGc1zZRXHxVkKoKLIGeMHz/4RzzGTeKqX8IYDz/fMh4X5GCHFYCPGYEMJUwaNS/q0wsjnU4DWaQwlMR7Vu/brTTfiGM+2k6qRqU8asp5LjdIBSibqvBNGslVK+JIS4DNgrhDgipXxWTw61oVL+Rj4jBVTy/XfAbinlT4QQn6apkf1y7jnLl1TPyxihIvObDlAqUfclhHhZCHGRlPKHjlr5o5A4XnK+jzqTN3tp9vFNQqX83TAnhBDLgJ8jWv02hdh7k1Ke9Pz9XzSX9ag6qdpUO3V/AqcDlJynOPbQnJ4AIdMUhBArhRAXOL+7gOuBfywsh+qolL/3fm8B9krHImg4sffmszUM0Bw9XnX2AJ9wvEDvA15zu+uRlG2BVrRSf5im1PwJ8DIw6Ry/GJjwWat/QPMt/rmy861wX6uAb9Nc+e7bwIXO8fXAl53f1wFHaHocjgCfLDvfEfezpPyBPwQGnN8/A/w1zRUADwGXlZ1njff2xzRXQfwezTltPWXnWeGedgM/BE477euTwKeBTzvnBc0F1p516l7o6o7ejx2mb7FYtNJO3R+LxWIAVqhYLBatWKFisVi0YoWKxWLRihUqFotFK1aoWCwWrVihYrFYtPL/AaCZNLCAwxsdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "########################\n",
    "samples = 1000\n",
    "delta   = 1.0\n",
    "mcmcplot(samples,delta)\n",
    "########################\n",
    "samples = 1000\n",
    "delta   = 0.01\n",
    "mcmcplot(samples,delta)\n",
    "########################\n",
    "samples = 1000\n",
    "delta   = 20.0\n",
    "mcmcplot(samples,delta)\n",
    "########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### (i)\n",
    "\n",
    "The graphs above have been plotted for `samples=1000` and for `delta=1.0, 0.01 & 20.0` respectively. Notice that for `delta=1.0` the points are distributed uniformly over the entire area of the unit disk, while for the other two values (`delta=0.01` and `delta=20.0`) the 1000 points (samples) show peculiar distributions within the same area. In particular, for `delta = 0.01`, the points are all concentrated near the beginning of the axes (0,0), which is due to the jumps z according to the formula: `z = r. uniform(-delta,delta)` and `x_prop_k = x[k] + z`, i.e. the variable z is a variable which takes values (randomly w.r.t the uniform distribution) within the interval (-0.01,0.01). Therefore, with such small jumps the points appear concentrated in a neighborhood of (0,0). As for the case where `delta = 20.0`, in the corresponding graph the experimental points show a preference on the `x=0` and `y=0` axes over the other points on the unit disk. That is, for large values of delta, the experimental points seem to cluster on the `x=0` and `y=0` axes, which may be due to the interval (-20.0,20.0) in which the variable z takes values. All the numbers within this interval have the same probability of being chosen according to the uniform distribution, in this case only the numbers that satisfy the condition: `R_sqprop < 1.0`, get accepted, where `R_sqprop` is the square of the distance of the (new) point after the jump z from the origin of the axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii)\n",
    "import numpy as np\n",
    "import random \n",
    "\n",
    "###########################################################################################\n",
    "def mcmc(samples,delta,d): ## d = dim\n",
    "    N       = 100          ## number of steps to approach equilibrium\n",
    "    nhits   = 0\n",
    "\n",
    "    x = np.zeros(d)   ## start at the centre of the disc. This variable will keep the position of the chain\n",
    "    R_sq = 0.0        ## this variable keeps the squared distance from 0. It saves some computations to keep it\n",
    "    for i in range(samples):\n",
    "        for j in range(N):\n",
    "            k = random.randint(0,d-1)               ## choose a jump direction at random\n",
    "            z = random.uniform(-delta,delta)        ## choose a jump size uniformly in (-delta,delta)\n",
    "            x_prop_k = x[k] + z                     ## propose a jump by z in the direction k\n",
    "            R_sqprop = R_sq - x[k]**2+ x_prop_k**2  ## compute the squared distance from 0 after the proposed jump \n",
    "            if(R_sqprop < 1.0): \n",
    "                R_sq = R_sqprop\n",
    "                x[k]= x_prop_k                      ## if the proposed jump leads to a point in the disk, then jump\n",
    "                nhits+=1                            ## checking if R_sqprop < 1.0 after equilibrium. \n",
    "    return(nhits/(N*samples))                       ## probability of the speciment to hit the untitary sphere of d-dimensions\n",
    "###########################################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  1 is:  0.75002\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  2 is:  0.68749\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  3 is:  0.63839\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  4 is:  0.60342\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  5 is:  0.57152\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  6 is:  0.5401\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  7 is:  0.51492\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  8 is:  0.49318\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  9 is:  0.47387\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  10 is:  0.45509\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  11 is:  0.443\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  12 is:  0.42763\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  13 is:  0.41308\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  14 is:  0.39851\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  15 is:  0.38947\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  16 is:  0.38037\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  17 is:  0.36777\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  18 is:  0.35836\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  19 is:  0.34988\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  20 is:  0.34401\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  21 is:  0.33303\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  22 is:  0.32853\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  23 is:  0.31933\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  24 is:  0.316\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  25 is:  0.31116\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  26 is:  0.30684\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  27 is:  0.29728\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  28 is:  0.29192\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  29 is:  0.28855\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  30 is:  0.28411\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  31 is:  0.28009\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  32 is:  0.27248\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  33 is:  0.27209\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  34 is:  0.27079\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  35 is:  0.26322\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  36 is:  0.26325\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  37 is:  0.25675\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  38 is:  0.25324\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  39 is:  0.25281\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  40 is:  0.24789\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  41 is:  0.24551\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  42 is:  0.24242\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  43 is:  0.23932\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  44 is:  0.23903\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  45 is:  0.23263\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  46 is:  0.23505\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  47 is:  0.22993\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  48 is:  0.22887\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  49 is:  0.22226\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  50 is:  0.22541\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  51 is:  0.21825\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  52 is:  0.2174\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  53 is:  0.21717\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  54 is:  0.21344\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  55 is:  0.21105\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  56 is:  0.20839\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  57 is:  0.20993\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  58 is:  0.20506\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  59 is:  0.20653\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  60 is:  0.20088\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  61 is:  0.2025\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  62 is:  0.19708\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  63 is:  0.19562\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  64 is:  0.19691\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  65 is:  0.19607\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  66 is:  0.19519\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  67 is:  0.19115\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  68 is:  0.19377\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  69 is:  0.1918\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  70 is:  0.19018\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  71 is:  0.18655\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  72 is:  0.18688\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  73 is:  0.18692\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  74 is:  0.18441\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  75 is:  0.18212\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  76 is:  0.18242\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  77 is:  0.17997\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  78 is:  0.17751\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  79 is:  0.17603\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  80 is:  0.17706\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  81 is:  0.17512\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  82 is:  0.17787\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  83 is:  0.17192\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  84 is:  0.17005\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  85 is:  0.17255\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  86 is:  0.17142\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  87 is:  0.16883\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  88 is:  0.16843\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  89 is:  0.168\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  90 is:  0.16462\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  91 is:  0.16407\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  92 is:  0.16393\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  93 is:  0.16359\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  94 is:  0.16106\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  95 is:  0.16037\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  96 is:  0.16048\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  97 is:  0.16279\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  98 is:  0.15951\n",
      "The probability of the speciment to hit the untitary sphere of d-dimensions =  99 is:  0.15956\n"
     ]
    }
   ],
   "source": [
    "for d in range(1,100):\n",
    "    print(\"The probability of the speciment to hit the untitary sphere of d-dimensions = \", d, \"is: \", mcmc(1000,1.0,d))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}