Computing_pi_throwing_dart.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import pi,sqrt\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The number _pi_ ($\\pi$)\n",
    "The number $\\pi$ is a mathematical constant. Originally defined as the ratio of a circle's circumference to its diameter, it now has various equivalent definitions and appears in many formulas in all areas of mathematics and physics. It is approximately equal to 3.14159. It has been represented by the Greek letter __\"$\\pi$\"__ since the mid-18th century, though it is also sometimes spelled out as \"pi\". It is also called Archimedes' constant.\n",
    "\n",
    "Being an irrational number, $\\pi$ cannot be expressed as a common fraction (equivalently, its decimal representation never ends and never settles into a permanently repeating pattern). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is the value of $\\pi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value of pi:  3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "print(\"Value of pi: \",pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the logic behind computing pi by throwing dart randomly?\n",
    "Imagine a square dartboard.\n",
    "\n",
    "Then, the dartboard with a circle drawn inside it touching all its sides.\n",
    "\n",
    "And then, you throw darts at it. Randomly. That means some fall inside the circle, some outside. But assume that no dart falls outside the board.\n",
    "\n",
    "![boards](https://raw.githubusercontent.com/tirthajyoti/Stats-Maths-with-Python/master/images/boards.png)\n",
    "\n",
    "At the end of your dart throwing session, you count the fraction of darts that fell inside the circle of the total number of darts thrown. Multiply that number by 4.\n",
    "\n",
    "The resulting number should be pi. Or, a close approximation if you had thrown a lot of darts.\n",
    "\n",
    "The idea is extremely simple. If you throw a large number of darts, then the **probability of a dart falling inside the circle is just the ratio of the area of the circle to that of the area of the square board**. With the help of basic mathematics, you can show that this ratio turns out to be pi/4. So, to get pi, you just multiply that number by 4.\n",
    "\n",
    "The key here is to simulate the throwing of a lot of darts so as to make the fraction equal to the probability, an assertion valid only in the limit of a large number of trials of this random event. This comes from the [law of large number](https://en.wikipedia.org/wiki/Law_of_large_numbers) or the [frequentist definition of probability](https://en.wikipedia.org/wiki/Frequentist_probability).\n",
    "\n",
    "See also the concept of [Buffon's Needle](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Center point and the side of the square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Center point\n",
    "x,y = 0,0\n",
    "# Side of the square\n",
    "a = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to simulate a random throw of a dart aiming at the square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def throw_dart():\n",
    "    \"\"\"\n",
    "    Simulates the randon throw of a dirt. It can land anywhere in the square (uniformly randomly)\n",
    "    \"\"\"\n",
    "    # Random final landing position of the dirt between -a/2 and +a/2 around the center point\n",
    "    position_x = x+a/2*(-1+2*random.random())\n",
    "    position_y = y+a/2*(-1+2*random.random())\n",
    "    \n",
    "    return (position_x,position_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.21862368906636753, 0.04575803244910337)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "throw_dart()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to determine if the dart landed inside the circle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_within_circle(x,y):\n",
    "    \"\"\"\n",
    "    Given the landing coordinate of a dirt, determines if it fell inside the circle\n",
    "    \"\"\"\n",
    "    # Side of the square\n",
    "    a = 2\n",
    "    \n",
    "    distance_from_center = sqrt(x**2+y**2)\n",
    "    \n",
    "    if distance_from_center < a/2:\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_within_circle(1.9,1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_within_circle(1.2,1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_within_circle(0.4,-0.74)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.5129931631226134 0.3782598009424025\n",
      "This one landed inside the circle!\n"
     ]
    }
   ],
   "source": [
    "r1,r2=throw_dart()\n",
    "print(r1,r2)\n",
    "if is_within_circle(r1,r2):\n",
    "    print(\"This one landed inside the circle!\")\n",
    "else:\n",
    "    print(\"This one did not land inside the circle!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now throw a few darts!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_throws = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_inside_circle=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_throws):\n",
    "    r1,r2=throw_dart()\n",
    "    if is_within_circle(r1,r2):\n",
    "        count_inside_circle+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the ratio of `count_inside_circle` and `n_throws`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio = count_inside_circle/n_throws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Is it approximately equal to $\\pi$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.2\n"
     ]
    }
   ],
   "source": [
    "print(4*ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Not exactly. Let's try with a lot more darts!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_throws = 10000\n",
    "count_inside_circle=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_throws):\n",
    "    r1,r2=throw_dart()\n",
    "    if is_within_circle(r1,r2):\n",
    "        count_inside_circle+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio = count_inside_circle/n_throws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.1476\n"
     ]
    }
   ],
   "source": [
    "print(4*ratio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now we are approaching $\\pi$ :-). Let's functionalize this process and run a number of times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_pi_throwing_dart(n_throws):\n",
    "    \"\"\"\n",
    "    Computes pi by throwing a bunch of darts at the square\n",
    "    \"\"\"\n",
    "    n_throws = n_throws\n",
    "    count_inside_circle=0\n",
    "    for i in range(n_throws):\n",
    "        r1,r2=throw_dart()\n",
    "        if is_within_circle(r1,r2):\n",
    "            count_inside_circle+=1\n",
    "            \n",
    "    result = 4*(count_inside_circle/n_throws)\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now let us run this experiment a few times and see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed value of pi by throwing 3 darts is: 2.6666666666666665\n",
      "Computed value of pi by throwing 10 darts is: 4.0\n",
      "Computed value of pi by throwing 31 darts is: 3.3548387096774195\n",
      "Computed value of pi by throwing 100 darts is: 3.16\n",
      "Computed value of pi by throwing 316 darts is: 3.1645569620253164\n",
      "Computed value of pi by throwing 1000 darts is: 3.148\n",
      "Computed value of pi by throwing 3162 darts is: 3.135989879822897\n",
      "Computed value of pi by throwing 10000 darts is: 3.1736\n",
      "Computed value of pi by throwing 31622 darts is: 3.1553981405350706\n",
      "Computed value of pi by throwing 100000 darts is: 3.14972\n",
      "Computed value of pi by throwing 316227 darts is: 3.143653135247781\n",
      "Computed value of pi by throwing 1000000 darts is: 3.144012\n",
      "Computed value of pi by throwing 3162277 darts is: 3.141502151772283\n",
      "Computed value of pi by throwing 10000000 darts is: 3.1412476\n"
     ]
    }
   ],
   "source": [
    "n_exp=[]\n",
    "pi_exp=[]\n",
    "n = [int(10**(0.5*i)) for i in range(1,15)]\n",
    "for i in n:\n",
    "    p = compute_pi_throwing_dart(i)\n",
    "    pi_exp.append(p)\n",
    "    n_exp.append(i)\n",
    "    print(\"Computed value of pi by throwing {} darts is: {}\".format(i,p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAAF1CAYAAACef1IVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8VNX5+PHPkxCWAIIQCXtQKSjqt2wiCMjivuBWWxdcqFWsYqvV1g3rUqto3a0rxUoLWLRaRVT0J5ssgoqKigiICgmLCLIoBAhJnt8f584wmcwkd8hsmTzv12teM3OXc8+9s9znnnPuOaKqGGOMMaZuyUp1BowxxhiTfBYAGGOMMXWQBQDGGGNMHWQBgDHGGFMHWQBgjDHG1EEWABhjjDF1kAUAxpgqicgqEVmV6nzUlIjMFpGk3fcsIneIiIrI4BjW6eStMz5xOTPGsQDA1CkicoiI/F1ElojINhEpEZF1IvKGiPxGRBqmOo/JluwTY13nneBnpzofxtRLdQaMSRYRuQ24HRf4LgT+BWwH8oHBwDjgSqB3irKYro5NdQbi5GIgN4nbexyYDBQmcZvG+GYBgKkTROQW4E6gCPilqr4fYZnTgOuTnbd0p6pfpzoP8aCqST0Rq+omYFMyt2lMLKwKwGQ8EekE3AHsAU6JdPIHUNXXgZMirP8rEZnjVRnsFJHPReRmEWkQYdlV3qOJiDwsIkXeOotF5ExvmXoicouIfCUiu0TkaxG5OkJag73i4jtEpJ+ITPfy8JOIvC0ilUoqRGS8t06nqtILHBev6H+Q915DHrPD9yksrRHeciNEZIhXjfCTiPzoVaccGukYi0gXEXlZRLaIyA4ReU9ETg1NL9J6EdIJ1q+LyCUi8ol3nL8XkX+KSOsI6/iu6hCRBSKyW0Qah02f42332bDp3bzp/46UR+/9iJDtDwo73ndEyEMnEZksIpu878kiL0g1Ji6sBMDUBb8GcoDJqrqkqgVVdXfoexG5B7gZdyX3PK7K4GTgHuBEETleVfeEJZMDvAO0AKYA9YHzgZdF5ATgKuAoYBqwG/gl8HcR2aiqL0TI1lFeHqYDTwCdgbOBY0TkBFWd6+soVLYVVyoyAijwXges8pnGacAZuH15GugGnAIcKSLdvKtgwLW/AObjjssbwGfAQcArwJv7uA9/AE4AXgDeAgbgPu/BInKUqm7cx3RnAH2BgV66iEgu7rOAytUiQ0PWi2Yx7hjfDqwGxofMmx22bAHwAfANMAF3zM4FpojIcao6y/+uGBOFqtrDHhn9wP0pK3BZjOv189YrBFqHTK8HTPXm3RK2zipv+lSgQcj0gd70zcCHQPOQeQcBJcAnYWkN9tZR4OqweWd4078CskKmj/emd4qwP4H07gibPtv9FUQ9DquAVWHTRnhplQLHhs0b4827IcrncGXY9JND9nOEz8/mDm/5EqBH2LyHvXnPxrKfYcsO9dK4P2Taid60/+c9Hxwy7xVvWocIeRwclrYCs6Nst1PIsbg9bF5g+2+m8vdkj8x5WBWAqQvaeM9rYlzvUu/5r6r6XWCiqpbi2gqUA5dFWfdaDSlNUHeV/i2wP3Cjqm4NmfcN7sr4CBHJjpDWSuDJ0AmqOgV4F1caMDDG/YqnyaoaftU71nvuE5ggIh1wJ9WVwDOhC6vqNFzpxr6YoKqfhE27A9gGXBCpmsan94BdVLzSPxYX8Nwe8h4RycIFV1+patE+bi/cauCvoRNU9W1cMNon4hrGxMgCAFMXiPcc661uPb3nmeEzVHUFLqA4UESah83eqpEbzq3znj+KMG8tkA1UqrsG5qpqeYTps73nHhHmJcuiCNMCJ8H9Q6Z1954XRNmXefu4/XfDJ6jqNlxxe0MgYluE6qjqLlwQ0F1EWnqThwIfquoCYAN7g4OeQHMifE9qYLGqlkWYXkTF42rMPrMAwNQFgRNv+xjXa+Y9r48yf33YcgHboixfCsETVMR5uPYD4TZESS9QKhG+/WTaGj7BKyEBF9AEBPIYbV+iTa9OIo/NDFzwOERE9scFWoHSjpnAUBER9gYCVdX/x6rScfWUYv/bJk7si2TqgsDVZaz3swdO1JGuymFv1UK0E3685EeZHshX6PYDV9eRGviGl1Qk04/ec7R9iTa9OrEcm1gFruiPA4bg/i9DA4A84Oe475UC1jDP1CoWAJi64DncLYC/EJFuVS0YVmccqFseHGG5zrgShW9D6/MTZIBXzxwukK/QOvAt3nOHCMtH6+CoDCBK+4N4CeSxX5R9GbCP6Q4KnyAizXBVDruAL/cxXXCNNX/EneCHAjuBBd68QCBwCtAf+ExD7nioRjkVS0eMSQkLAEzGU9VVuIZh9YE3It0/DyAiJ+FuZwv4p/d8q4gcELJcNvAA7vdT4X7wBPkZ7tbBIBE5A3fyWwmE3gb4gfd8edjyRwDXREn/B++5Y41zGoXXOG42rtHiFWF5Owl3lb0vLhKR8DYQd+CK/v+jYbd1xsKrg5+Dy/MvgXmB9FT1W9zdEdfgeheMpf7/ByIHaMYklfUDYOoEVb1HROrhWnB/KCLv4RqwBboCPgZ3ol0Uss57IvI34AZgiYi8BOzA3bZ2OK5q4f4kZP8t4EERORn4lL39AOwCfhPWqG4K7tbA80WkPfA+7sR+hjfvVxHSn4E7wf1PRN7EXemuVtUJcd6PUbi7HZ4UkVPY2w/AL7y8ncHeKgy/pgHzReRFXJuMAd5jFXBTHPI8A9fXQSsq1/HPAH4T8jqWNM8Tkam4BqGlwBxVnVPDvBoTEysBMHWGqv4Fd+J+HHeF+GvgT8CpwNe4W/oGhK1zI64Tn69wfcn/Hve7uRU4XlVLkpD193HF/Q2Aq3EByEzgmPCThtd6/VjgRdy+Xo07yV4APBUl/XG4e/eb4YKdu9h7YosbVV2K61vhFdyti9fi7ns/i73tNH6MuHJ0D+NKR7p76R2C6wvhaFX9vsaZrnhiD7/KD8wrxZUU+HUN8B/c7Xx/xh3voVWuYUwCiKoNAmZMOvK6kJ0F3Kmqd6Q2N4klIpNwQcohqrrcx/J34Epzhqjq7MTmzpjMZCUAxpikEJGsKH30H4vr5napn5O/MSY+rA2AMSZZ6gNFIjILWIYrOj8MOB7Xpe+oFObNmDrHAgBjTLLswQ0YNBQ3qE4ubpCl/wL3RujS1xiTQNYGwBhjjKmDrA2AMcYYUwdZAGCMMcbUQRYApICIdBIRFZHxqc5LuhKREd4xGpHqvNR2IrJKRFalOh/xJCK9ReQdEdnkfU8WpzpPNSUid3j7MjjVeYmHTPzeZRprBGiMqVVEZD/gDdxwvxNwDQm/q3IlE3ciMhsYpKpS3bImPVkAkBprceOUJ3oUudrsFWAh0YfiNXVXH1zXvKNV9Z5UZ8aY2soCgBRQ1T24+6BNFKq6DQuQTGRtved1Kc2FMbWctQFIgWhtAERkvDe9k4hcISKfi8guEdkgImO9YU4jpddeRB4Tka+85TeLyAci8uew5VZ5j/1E5CHv9R6vW9XAMvVE5CoRWSgiP4pIsYh8IiJXRxrG1aurf1lEvhGRnd4680Xkwih5Pcjbl5Xe8pu9/XxaRFqGpVupDUDIPuSKyP0iUigiu730bhSRSsWR4lwjIku947NWRB4XkWax1lN6eZotInnefqz3tv+FiPw6yvGJ2pYhkF7YtGBdsIicLyIfeZ/DOu9za+AtN9TLy48iskVEJoQewwjbaubt91rvOCwVkd9HOmbe8keJyEsi8p2IlIhIkYg8IyJtIyw728tzfRG5TUSWe8dlfJUHdO/6x4rIW973YZeIrBCRe0O/897vQoF/eZOe87ZZbVsR71iqd2z7iMgb3rZURDp5ywzxPtOl3jHdKSJLROR2EWkYIc3Qz+kccb+5Yi/dySLSLkpeenn7+pO3neki0q+mxydk2cBnkeN9Fl976ywTkctDlvutuN/eThFZIyJ3SuShmsPTD3wOg7z3GvKYHWF5X79VCflfFJEuIvKCiHwvIuUS0i5CRH4mIv/2vscl3u/i3yLys7D0rvDSCx8Z81JverFUHP4b7zPcJSKNQqadLiIzZO9vfZ2IvCsiFUborI2sBCA9/Q04EZgK/D9gCG54186EDRoibmjbt4EWuAFJ/ofrYKUbbljUu8LSro8b1KSFl/aPwLdeWjneNk8ElgPP40acGwL8Hdd5y0Vh6T0FLPW2vR5oiRsjfYKIdFXVYBAiIm1wY6zvB7wJvIyrxz3QS/dx9g5NW5UcL+9tcaPBlQJnAvd66d0ZtvwTwJW4K8axuF7nTscVJefgOqiJRXPcqHYlwEveNs8B/iki5ar6r6pWjsHvcAP/vIobSvcE4A9ACxGZAkzG1YWPBY4GLgTyvHXC1Qeme3mf7L3/BfAo0JWwXvjEBTP/AHYDrwFFuNESLwOGiUhfVS2MsJ2XgSNxn8urQLUD8ojIFbjv0Q5cp0Df4wY/utHbVn9V3QpsxX223dk7umGg8Z/fRoD9gJtxgw/9E3e8AgM63YgbTOg99rYx6I/7HQ0WkeO8IYLDXYX7Pr0GvIv7nZwL/FxEuocOSSwiR+M+h/q43+pKb39mE2VI4RiOT7jJXl7exH3HzwHGisge4P+AS4DXcYManQ7cBhQD90XKR4jA5zACKKDi721V2LKx/lYBDsYNgLUCmAQ0whskSkSOxB2/prjjvRT3mQ0HzhCRY1U1MKJnYLCmY3Hf5YDAf2gj3Pdhtpd2M6AnMFdVd3rTRgLP4NqYTMW1N2mFO36/Bp6MeIRqC1W1R5IfuBHQFBgfNn28N70Q6BgyvR7uBKtAn5Dp9XEnbwUuiLCdDmHvV3nLTgcaR1j+Dm/+34HskOnZuHHvFTgjbJ2DI6RTH/fj2wO0C5n+Oy+NayKs0xhoFPJ+hLfsiCj78GbY8q1wf0xbgZyQ6QO95ZcDzcPyGDimq2L47NR7jAs7Rt1wf25Lw5aPuB9h6c2O8jlsAw4Nmd4A+AIowwVKg0LmZQHveOt1j3LM5gENQqa3wI2CqLiRBQPTu+BOiitDPz9v3lBv+6+ETZ/tpfMZkBfD8SzABRk/4gYCCp33pJfm2FiOaZTtDA757K6IssxBeJ2jhU2/y1vv3Cif04/AEWHznvfm/SpkmuCq/iL9jq4Jyd/gGh6fwGfxYdh3/iDvc92C+98I/W02x53cNgL1fB7T2YBWMT/wvfP7W+0UcgzuiZCeAF9684eHzTvXm74MyAqZvhoXMEnItHW4/6cy4K6Q6Wd4afw5ZNpH3vFvFSE/vr/n6fqwKoD09BcNubpS1VLgOe9tn5DlhuF+NK+p6vPhiahqUZT0r1fVHaETvKK/q3GR7h805ErHe3093g8vbBtfR9huCe6qux4u+g63M8I6O9SLun36fejy6oZ+nYIb0rZryHKXeM93a8hVkpfHm2PYXqhi4LqwY7QUVypwqIg03cd0wz2mql+GbGM38ALuZP+Gqr4bMq8cmOi9/XmU9G7WkKtRVd3M3hKi0OqLK3FXbteo6trQBFR1Ju7Ka1iU/fyzqm7ys3OeC3HB2OOqGt4uZjTwE3BReFFtDSxW1WcizVDVb9T7Zw/ziPd8YpQ0H1PVz8OmBa44Q3+vR+O+m3NUdUrY8o/jgrFwNTk+N4V957/BBYHNcSe+tSHztuKucPOAiFUXNeD3txqwgcglA0fjrvYXqOqk0Bmq+gJu37pScUjvmcABwBEAItINaIMrufuYiv9PgdehQ0CDC+wrlRLG+D1PS1YFkJ4WRZgWOJnvHzKtr/c8LYa0d+Gu0sJ1wRXffwXcKpGrhXfi7l4IEpGOuKLIY4GOuGK1UKF/Jq8B9wBPiMiJuKqL+bir5lj6pN6mqisjTI90jHp4z/OobCHuxx2rr1Q10rj1ge03x/0x11Sk70Gg4dtHEeYF/tDbR5hXiivaDjfbe+4RMi1QHz3IK3IN1wpXKtQlQj4+iLB8VXp6z5WKv1V1i4h8AhyD++P/NMa0I4maPxFpjLsSPwu3b01xV50B0U6Mfn+vgX19N2xZVLVMRObhir9D1eT41OT7szrC/H0Ry2814NPQQDVE1GMRMn0A7rs8J2TaCNz/02fsLf6fgbt4uk5EmqrqT9687VT8jkwCHgS+EJEXcJ/dfFXdGCUPtYoFAOkpUn1e4ESVHTKtufe8Fv++j3KyDTQe+xlunPVomgReiMhBuB/L/sBcXF3fNlzRWifc1XfwykRVV4tIH1zR6UnA2d6sIhF5QFUf87kPkY4PRD5GgUZSG8IX9v50/bQ5qMn2ayLSXRClPublRJi3SSPXXwfunw9tTBb4Lvypmvw1iTAt1vvxA9uNdrtnYHrzKPNjFTF/XvuXmbgr9iW4kpaN7L3yu52Q73IYv7/XqN/FKvK2z8dH3Z000fIV6/dnX+3LbyXad2hfjkVoO4CHvec1qrpCRGYAN+AC3Q9xI1O+6ZW4AqCqD4nIJlw7j98D1wIqIu8Cf9K97Q1qJQsAarfAjyuWIrtoV9qBP4RXVPXsKMuEuw53svi1qo4PnSEi57O3+H3vxl2R9rkiUg9XVH0crm3AoyKyQ1Wf9bltvwJX6vnAN2F5zPbyH0sAFaty77nSb01E4nVS8yNPRLIjBAGtvefQE0LgdbMoJR1RxViSE7qt1rj2DeHahC1XU9Hydwbu5P8vVR0ROsNrvFpVUOxXYB/yo8xvHWFaso9POqjuPyrScYIIx0JV14nIctxJvgGuLUig+mUerk3EcbiGyRC5pOXfwL+93+vRuBKiS4G3ReRQr0qjVrI2ALXbQu85UqvvWC3DBRR9vashPzp7zy9HmDeoqhVVtVRVP1LV+4Dzvcln+txuLAJDzA6IMK8viQ+Ct3jPHSLM653gbYeqh/vzCjfYew4dijfwvRqYyAyFbXdw+AzvD7c7rtrqy/D5cbbP3+UYfBwtPS8YjfQdTZfjE0mZl494lXhVJ+qxCJv+cdj0GbjSqitxpQMzAFS1GPddP5aKVQMRqepWVX1TVS/HNdhuQXJ+IwljAUDtNhXX0vZ074q7AolyH3IkXrHX33FR9GOh98GGpNfGa0QTsMp7Hhy23Im428XC1+8jIpGufgLTiv3mNwb/9p5HS8V7yuvj2iMk2iJcKcAFIpIbsv0WuNs9k2lMaGMxLw+3em+fC1nucVzR98Mi0iU8EXH3+sfrj2+it63fiUjnsHl34a7MJkapE46nVd7z4NCJXjVXdbfF+fUe7m6UY0TkjLB5V1O5/h/S5/hEEqg+65ik7c3HHb8BInJO6Azv/TG4WwfD2/sErupvDnsfeH047jbIHwhrRyEiJ3mlleFaec+J+M9KGqsCqMVUtUREfomre3/eu194Ie7+2kNxkW0sn/FduGL53+Jaec/EFY+3wrUN6I9rebzUW/5JXOvx/4rIy96yh+Pq91/E3ZoT6gJglFd/thJ3dXww7m6G3extbR03qvquiIwFRuIa8ryM+0MdhisqXMfeYvq4U9X1IjIJ18/BYhF5A/enfQquoVKPqtaPo/W4OuwlIvIarp73HFzA96SqBhpNoarLRORS3H3yX4jIW7g/1hzcn/1AXP34ITXNlKquEpFrcXeNfCwiL3ppD8I1RlyGa2SaaFNx38nrROQI3NVmR+A0XJ8ANT7JqaqKyG9wt2u+LCKBfgACVWFv4X47oeuky/GJZAbwS+B/IvImrpHwalWdkIiNecfvEtzxe8HrC2MZruX/mbiGtxd7d8SEmoX7jbcCloXd2TID1ybpAOClCFVYk4FdXgPNVbhGoQNxfV18hLulutayAKCWU9VFItIduAlXFXA07oewkhjrLVV1j4icibv1aATuz68J7g/nW+DPuFaxgeU/E5EhwF9xJ7R6uAj6bFx1QngA8B/cSehoXIveRrigYTLwoKouiSW/MbgS90dxBS64+QE31sAtwBoi334VT5fjGn6dj+twpxB4DLgf+FWCtx0QqOu8BzgPd7vXN7gOWf4evrCqThSRT3G3fw7BdUK0AxcwvYRrJBcXqvqkiKwE/ojrnCgX10r8ftz94NEaksWNqu4QkaG44zEY9yf/DS4ofojK3+V93c58r/TkbvZW3b3vbfNEwgIAb52UH58oxuH6KTgP15iuHq6VfEICAABVfd+7M+VW3Pd5GK7/gv/gbm1cHmGdzeJGi+xJ5Tr+93Hf68YR5oH7Xz3RW/cUXHXLalzQ9ZS6bt1rLYm9zY4xmUFc16ErgMmqWqkKxRhjMpm1ATAZT0RaS1gf5159fKDK4ZXk58oYY1LLqgBMXXAtcL64gUrW424jOhbX4ck0XP/qxhhTp1gAYOqCd3ANrU7A3bpTiiv6fwx4ZB/uXTfGmFrP2gAYY4wxdVBGlwDk5eVpp06d4prmjh07aNy4cVzTTEe2n5nF9jOz2H5mlnjv50cffbRJVQ+obrmMDgA6derEokXx7ap59uzZDB48OK5ppiPbz8xi+5lZbD8zS7z3U0R8DeZkdwEYY4wxdZAFAMYYY0wdZAGAMcYYUwdZAGCMMcbUQRYAGGOMMXWQBQDGGGNMHWQBgDHGGFMHpTQAEJFbRERF5PFqljtCRN4VkZ0islZEbhMRSVY+jTHGmEyTsgBARPrixkn/rJrl9sP15b4BOBL4PfAn4LpE57EumjRpEp06dWLo0KF06tSJSZMmpTpLxhhjEiAlAYCINAMmAb8BtlSz+HAgF7hEVZeo6svAfcB1VgoQX5MmTeKyyy5j9erVqCqrV69m5MiRFgQYY0wGSlUJwFjgJVWd6WPZfsBcVd0ZMu1toC3QKQF5q7NGjx7Nrl27KkwrLi5m9OjRKcqRMcaYREn6aIAicjnwW6CfqpZ4Y7QvUdWroyz//4A1qnppyLSOwGrgaFVdELb8SGAkQH5+fq/JkyfHNf/bt2+nSZMmcU0zXQwdOpRI3wcRYeZMP7Fa7ZPJn2co28/MYvuZWeK9n0OGDPlIVXtXt1xSBwMSka7APcBAVS2JYdXws5JEmY6qjsWVMNC7d2+N90ASmTw4RYcOHSgsLKw0vWPHjhm7z5n8eYay/cwstp+ZJVX7mewqgH5AHrBEREpFpBQYBFzlvW8QYZ3vgNZh01p5zxsSl9W659prr600LTc3l7vvvjsFuTHGGJNIyQ4AXgWOALqHPBYBk73XkUoFFgADRaRhyLTjgXXAqkRmtq75+c9/XmnanXfeyfDhw1OQG2OMMYmU1ABAVbd6LfmDD2AHsNl7ryIyRkRmhKz2PFAMjBeRw0XkbOAm4CFNdgOGDLdq1apK05o3b578jBhjjEm4dOwJsA1wcOCNqm7DXfG3xZUWPAE8CDyUktxlsNWrV1eaNnfu3BTkxBhjTKIltRFgJKo6OOz9iAjLfA4ck6Qs1VkWABhjTN2RjiUAJkUiVQF8++23rF27NvmZMcYYk1AWAJig0BKAvLy84GsrBTDGmMxjAYABoLS0lKKiouD7Y489NvjaAgBjjMk8FgAYANatW0dZWRkA+fn59OrVKzjPAgBjjMk8FgAYoGLxf0FBAYcddhhZWe7rsWTJErZsqW7MJmOMMbWJBQAGqNgAsFOnTuTm5tKjRw8AVJX58+enKGfGGGMSwQIAA1QuAQAYOHBgcJpVAxhjTGaxAMAAlUsAwAIAY4zJZBYAGCByCcCAAQOC0xYtWsTOnTuTni9jjDGJYQGAASIHAK1ataJr164A7Nmzh/fffz8leTPGGBN/FgAYysvLIwYAYNUAxhiTqSwAMGzYsIGSEjcSc4sWLWjatGlwngUAxhiTmSwAMBEbAAaEBgALFiygtLQ0SbkyxhiTSBYAmKjF/+ACgnbt2gGwfft2Fi9enNS8GWOMSQwLAEyFEoDwAEBErBrAGGMykAUApkIJQHgVAFg7AGOMyUQWAJgqqwCgYgAwb948VDUp+TLGGJM4FgCYKhsBAhx22GHsv//+AGzcuJHly5cnKWfGGGMSxQKAOk5Vqy0ByMrKon///sH3Vg1gjDG1nwUAddymTZsoLi4GoGnTpjRv3jzictYOwBhjMosFAHVceANAEYm4nAUAxhiTWSwAqOOqK/4P6NWrF40aNQJcm4E1a9YkPG/GGGMSxwKAOq66BoAB9evX56ijjgq+t1IAY4yp3SwAqOP8lgCAVQMYY0wmsQCgjvNbAgCV+wMwxhhTe1kAUMfFUgLQr18/srOzAViyZAlbtmxJaN6MMcYkTlIDABEZJSKficiP3mOBiJxazTonesv9JCKbRGSKiHRJVp4zXSwBQJMmTejRowfg+g+YP39+QvNmjDEmcZJdArAGuBHoCfQGZgKvisj/RVpYRA4EpgBzgR7AcUAj4M2k5DbDbd26lW3btgHQqFEjDjjggGrXsXYAxhiTGZIaAKjqFFWdpqorVXWFqo4GfgL6RVmlF5AD3OytsxgYAxwsInlJynbGCr/6j9YHQCgLAIwxJjOkrA2AiGSLyHlAE+C9KIstAvYAl3nLNwUuAT5U1U1JymrGiqUBYMCAAQOCrxctWsTOnTvjnCtjjDHJUC/ZGxSRI4AFQENgO3CWqn4eaVlVXSUixwP/BZ7ABSyfACdXkf5IYCRAfn4+s2fPjmv+t2/fHvc0U2X69OnB1zk5ORX2q6r97NixI4WFhezZs4dnnnmG7t27JziniZNJn2dVbD8zi+1nZknZfqpqUh9AfaAzrg3AGGATcHiUZVsDK4C/4doAHAPM9h5Z1W2rV69eGm+zZs2Ke5qpct111ymggN5zzz0V5lW1n5dffnlwvb/85S8JzmViZdLnWRXbz8xi+5lZ4r2fwCL1cT5OehWAqpaoq89fpKo3A4uBP0RZfBSwQ1VvUNVPVHUOcCEwCDg6SVnOWPtSBQDWDsAYYzJBOvQDkAU0iDIvFygLmxZ4nw55r9ViuQUwVGgAsGDBAkpLS+OaL2OMMYmX7H4A7hWRgSLSSUSOEJExwGBgkjd/jIjMCFnlDaCniNwuIj8TkZ7Ac0AR8FEy856J9rUEoKCggPbt2wOu7mrx4sVxzpkxxphES/ZVdGtgIrAcmAEcCZysqtO8+W2AgwMLq+pM4ALgDFzjv7dxdwWcpKo7kpjvjLNjxw5++OEHwA3007p1a9/riohVAxhjTC2X7H4ARqhqgao2UNVWqnqcqr4dNr8vLUB5AAAgAElEQVRT2DqTVbWnqjZR1QNUdZiqLk1mvjNRaPF/hw4dyMqK7atgAYAxxtRuVo9eR+1r8X9A+MBAruGpMcaY2iJqPwAi8iKuB76vvddVUVU9N75ZM4m0rw0AA7p168b+++/Pli1b2LhxI8uXL+eQQw6JZxaNMcYkUFUlAAfguuEFaOW9j/ZolcA8mgSoaQlAVlZWhV4BrRrAGGNql6glAKo6JOT14KTkxiRNTUsAwFUDTJ06FXABwOWXXx6XvBljjEk8awNQR8UrAAiwEgBjjKldfAcA3n37z4vIShHZ4T0/H20oX5PealoFANCzZ08aNWoUTG/NmjVxyJkxxphk8BUAiMiZuI53egAvAX/2nnsAi7z5ppbYtWsX3333HQDZ2dm0a9dun9KpX78+ffv2Db63UgBjjKk9/JYA3AdMAbqp6k2q+pCq3gR0A6biBusxtURhYWHwdfv27alXb98HhbRqAGOMqZ38BgAdgHEadrO3934s0D7eGTOJE4/6/wALAIwxpnbyGwAsAg6LMu9w4OP4ZMckQ2j9f00DgL59+5KdnQ3AkiVL2Lx5c43SM8YYkxx+A4DrgKtE5EYR6Soi+3vPNwFXAteKSG7gkbjsmngILQHY1waAAU2aNKFnz57B9/Pnz69ResYYY5LDbwDwAXAQMAZYCmzynu/xpr8P/BTyMGksnlUAYNUAxhhTG/lt/XUpYJ29Z4h43AIYauDAgTz00EOABQDGGFNb+AoAVHV8gvNhkijeJQChXQIvWrSI4uJicnOtJsgYY9KZ9QRYx+zZs4e1a9cG33fo0KHGaebl5XHooYcCUFpayvvvv1/jNI0xxiSWBQB1zJo1aygvLwegbdu2NGjQIC7pWjsAY4ypXSwAqGPiXfwfYAGAMcbULhYA1DHxbgAYEBoALFiwgNLS0rilbYwxJv6iBgAiMlNEDvFeXywiLZOXLZMoiSoBKCgoCLYn2LFjB5988knc0jbGGBN/VZUADASae6+fAw5OfHZMoiWqBACsGsAYY2qTqgKAIuCXInI4IMCBItIt2iM52TU1lagSALAAwBhjapOq+gEYAzwJXIvrBOj5KMuJNz87vlkziZCsAGDevHmoKiIS120YY4yJj6gBgKr+Q0ReA34GzAFG4br/NbVUWVlZhaGA4x0AHHroobRo0YLNmzezadMmli1bFuwfwBhjTHqpsidAVd0AbBCRO4EpqrouOdkyibB+/fpg6/wDDjgg7r31ZWVlMWDAAF577TXAVQNYAGCMMenJ122Aqnqnqq4Tkfoi0ktEjvee6yc6gyZ+EtkAMMDaARhjTO3gux8AEbkB2IAbGfBt73mDiPwpQXkzcZbI+v8ACwCMMaZ28BUAiMi1uEaBzwNDgEO95+eBMSLye5/pjBKRz0TkR++xQEROrWYdEZFrRWSZiOwWkfUicq+f7ZmKkhEA9OzZM1i1sHr1aoqKihKyHWOMMTXjtwRgFHCvqo5S1Tmqutx7HgXcB/gKAIA1wI1AT6A3MBN4VUT+r4p1HgSu8tY7FDgF1yjRxCgZVQA5OTn07ds3+N5KAYwxJj35DQA6ALOizJsNtPeTiKpOUdVpqrpSVVeo6mjgJ6BfpOVFpCvwO+AMb91vVPUTVX3TZ75NiGSUAIBVAxhjTG3gNwAoBE6IMu94b35MRCRbRM4DmgDvRVnsDOAb4CQR+UZEVonIv0SkVazbM8kpAQALAIwxpjYQVa1+IZGrgceAfwIv4RoDtgJ+CYwAfq+qT/raoMgRwAKgIbAdGK6qb0RZ9mkv/U+BP+E6HHrAm91PVcsjrDMSGAmQn5/fa/LkyX6y5dv27dtp0qRJXNNMBlXlpJNOoqSkBIDXX3+dxo0bR12+Jvu5c+dOhg0bRllZGQCvvvoqzZo126e0Eq22fp6xsv3MLLafmSXe+zlkyJCPVLV3tQuqqq8HcDmuDr8cKPOe1wCX+U3DS6c+0BnXBmAMsAk4PMqyY3En/S4h07p4046qblu9evXSeJs1a1bc00yG7777Tr3jps2bN692+ZruZ58+fYLbmzJlSo3SSqTa+nnGyvYzs9h+ZpZ47yewSH2cj33fBqiq/8C1BSjA1dkXAB1UdZzfNLx0StS1AVikqjcDi4E/RFl8PVCqqitCpn0FlAIdY9luXZes4v8AqwYwxpj05jsAAPCCiyJV/cB7rr7+wF8eGkSZNx+oJyKhIxEehOvBcHXkVUwkyWoAGGABgDHGpLeYAoCaEpF7RWSgiHQSkSNEZAwwGJjkzR8jIjNCVpkOfAz8U0R6iEgPXDuE94FFycx7bZfsEoABAwYEX3/00Ufs2LEj4ds0xhjjX1IDAKA1MBFYDswAjgROVtVp3vw2QPBqX10jv9OA73H3/r+Na3dwhkZoAGiiS3YJQMuWLenWzY0SXVpayvvvv5/wbRpjjPGvysGA4k1VR8Q6X1XX4+42MDWQ7AAAXDXA0qVuAMm5c+cydOjQpGzXGGNM9ZJdAmBSJNlVAGDtAIwxJp3FFACIyP5eHf4FIrK/N62hiFggkcZUNWUlAAELFixgz549SdmuMcaY6vkdDChbRP6Gq39/F5gAHOjNfhm4PTHZM/GwefNmtm/fDkCTJk1o0aJFUrbbsWNHOnZ0d2sWFxfzySefJGW7xhhjquf3yv0eXEdAV+Nuw5OQeVOAYXHOl4mj8Kt/Eali6fiyagBjjElPfgOAi4GbVPU5IHx8169xQYFJU6H1/8kq/g+wAMAYY9KT3wCgOe5EH0l9IDs+2TGJEFoCkKwGgAGhAcC8efMoL7e7N40xJh34DQCW4Ebmi+RkXGc9Jk2logFgwKGHHkrLli0B+OGHH1i2bFlSt2+MMSYyvwHAX4ErRWQccBxukJfuInIXcAWujYBJU6m4BTBARCr0CmjVAMYYkx58BQCqOgW4AHfyn4ZrBDgON1TvRar6dqIyaGoulSUAYO0AjDEmHfnuCVBVXwReFJEuQB6wGVgepwGBTAKlshEgWABgjDHpKOYOfFR1haq+p6rL7OSf/n788Ue2bt0KQMOGDcnPz096Hnr06EFubi4AhYWFFBYWJj0PxhhjKvJVAuB1AlQlVb2h5tkx8RZa/N+xY8ek9gEQkJOTQ79+/Zgxww30OHfuXIYPH570fBhjjNnLbxVApMF49gf2A7YBWwALANJQKhsAhho4cKAFAMYYk0Z8BQCqemCk6SJyFDAW+G08M2XiJ9UNAAOsHYAxxqSXGg3io6rvA/cDj8cnOybe0qUEoG/fvtSr5+LNpUuX8sMPP6QsL8YYY+IzHPAPQNc4pGMSIF1KAHJzc+nVq1fw/bx581KWF2OMMf5HA8yN8GguIv2AvwBfJDabZl+lSwAAVg1gjDHpxG8JwHbgp7DHD8B8oDVwVUJyZ2osXaoAwAIAY4xJJ37vArgU1/1vqF3AGuADVd0T11yZuCguLmbjxo0A1KtXjzZt2qQ0P/379w++/vjjj9mxYweNGzdOYY6MMabu8nsXwPgE58MkQHgfANnZqR20sWXLlhx22GF88cUXlJaWsnDhQo499tiU5skYY+qqeDQCNGkqner/A6wawBhj0kPUAEBENorI934fycy08ccCAGOMMdFUVQXwBJXr/U0tkk4NAANCA4CFCxeyZ88ecnJyUpgjY4ypm6IGAKp6RxLzYRIgHUsAOnToQEFBAatXr6a4uJiPP/6Yo446KtXZMsaYOsfaAGSwdCwBAKsGMMaYdOA7ABCRfiIyTkTmiMgH4Y9EZtLsm3QsAQALAIwxJh347QnweGAO0B4YAGzEdQ70c6AlsMRnOqNE5DMR+dF7LBCRU32u+zMR+UlEtvtZvq7bvXs369evByArK4v27dunOEd7hQYA8+bNo7y8PIW5McaYuslvCcBfgEeBwMn6z6o6FOgC7AFm+0xnDXAj0BPoDcwEXhWR/6tqJRGpD0zGBSHGh6KiIlRdG8527dqlVUO7Qw45hLy8PAA2b97Ml19+meIcGWNM3eM3AOgGTAPKcXcGNAZQ1dXAHcBoP4mo6hRVnaaqK1V1haqOxnUr3K+aVe8DPgP+6zO/dV66Fv8DiAgDBgwIvrdqAGOMST6/AcAuIEvdJeV64OCQeT/iqgZiIiLZInIe0AR4r4rlTgVOA34f6zbqsnRtABhg7QCMMSa1/I4F8CluyN93gBnAzSKyFijBVQ987neDInIEsABoiGtHcJaqRlxfRNoA/wDOVtWfRMRP+iOBkQD5+fnMnj3bb9Z82b59e9zTTIR33323wvtY85zo/czNzQ2+nj59esqOaW35PGvK9jOz2H5mlpTtp6pW+wBOAUZ5r9sBH+OqA8qBQqCXn3S89esDnXFtAMYAm4DDoyw7A9feIPB+BLDd77Z69eql8TZr1qy4p5kIF198seKqa3Ts2LExr5/o/dyzZ482btw4mMdVq1YldHvR1JbPs6ZsPzOL7Wdmifd+AovUxznSVxWAqr6pqk94r9cCvXAlAt2Bzqr6UQwBR4m6NgCLVPVmYDHwhyiLDwVuF5FSESkFngUae+9H+t1mXZTuVQD16tWjX7+9TT+sGsAYY5LL722AQySk/N0LMr5S1c9UtSQOeWgQZd4RuCAj8LgN2Om9tgaBVUjnRoAB1g7AGGNSx28bgBnABhF5EXhBVaM22quKiNwLvAEUAU2BC4DBeLcXisgYoI+qHgugqkvC1u8NlIdPNxWVlpayZs2a4PuOHTumMDfRWQBgjDGp4/cugCNwjfFOBOaJSKGI3O+dkGPRGpgILMcFFUcCJ6vqNG9+GyreYWD2wdq1aykrKwOgdevWNGzYMMU5iuyoo44iK8t9Bb/88ks6dOjApEmTUpwrY4ypG/y2AfhCVW9T1UNwnfhMAs4CPhCRlSLyV5/pjFDVAlVtoKqtVPU4VX07bH6nKtYfr6pN/GyrLgut/0/X4n+AV155pcL7NWvWMHLkSAsCjDEmCWIeDEhVF6vqzaraGTgdaATcHPecmX0WWv+fjg0AA0aPHl2pG+Di4mJGj/bVr5Qxxpga8NsGIEhEWgBnA+cCg3CN8p6Pc75MDdSGBoAAhYWFMU03xhgTP37vAthPRC4RkTdxPQE+CmwBzgNaqepFCcyjiVG63wIYEK1xYuvWrZOcE2OMqXv8VgF8DzyN6/lvBO6k/ytV/Z+q7k5U5sy+qS0lAHfffXeFHgEDysrK2LJlSwpyZIwxdYffAOC3QL6qnqmq/1HVHYnMlKmZ2tIIcPjw4YwdO5aCggJEhEBXE99//z3nnnsupaWlKc6hMcZkLr93AYxX1R8TnRlTc+Xl5RQVFQXfp3MAAC4IWLVqFeXl5fz3v3v7dnrnnXe44YYbUpgzY4zJbDHfBWDS23fffUdJieucsWXLljRpUnvumvzFL37BbbfdFnz/8MMPM378+NRlyBhjMpgFABmmtjQAjOb222/nrLPOCr6/4oorWLBgQQpzZIwxmckCgAxTWxoARpOVlcW///1vDj/8cABKSko4++yzWbt2bYpzZowxmcUCgAxT20sAAJo0acKUKVNo0aIF4Ko1zjzzTHbu3JninBljTOawACDD1PYSgICDDjqI//73v2RnZwOwaNEiLr/8ctxQ18YYY2oqagAgIt+KyDd+H8nMtIkuUwIAgKFDh/LII48E30+aNIkHHngghTkyxpjMUVUJwMthjxygGfAB8Lr33AzXnfBLic2m8SsTqgBCjRo1issuuyz4/sYbb2TatGlVrGGMMcaPqGMBqOofA69F5Bbga+DU0E6ARKQJLhiwPgLSgKpmVAkAgIjwxBNP8OWXXzJ//nxUlfPPP5/333+frl27pjp7xhhTa/ltAzAKuD+8B0BV3Q484M03KbZx48ZgQ7lmzZrRvHnzFOcoPurXr8/LL79Mhw4dANi2bRunn346W7duTXHOjDGm9vIbADQD8qPMaw3Unt5mMlimXf2Hys/PZ8qUKTRq1AiAFStWcP7551NWVpbinBljTO3kNwB4DbhfRM4RkQYAItJARH4J3AdMTVQGjX+ZHAAA9OjRg+eeey74/q233uLmm29OYY6MMab28hsAXAnMAV4EikVkK1AMvADM9eabFMu0BoCRnHvuudxyyy3B9/fffz8TJ05MYY6MMaZ2itoIMJSqbgPOEpHDgCNx1QHfAR+q6tIE5s/EINNLAALuuusuPv/8c6ZOdQVPl112GV27duXII49Mcc6MMab2iKkjIFX9whsZ8D5V/Zed/NNLXSgBANdd8MSJE+nWrRsAu3fv5swzz2T9+vUpzpkxxtQevgMAEWklIveJyAwRWe6VBiAi14hIv8Rl0fhVV0oAAPbbbz+mTJnC/vvvD8C6des4++yz2bVrV4pzZowxtYOvAEBE+gBfAb8AVgGdgQbe7DbA9YnInPEvE/sAqE7nzp154YUXyMpyX+OFCxdy5ZVXWnfBxhjjg98SgIeBWUAX4ApAQuZ9APSJc75MjLZu3cqPP7r+mHJzc8nLy0txjpLj+OOP58EHHwy+Hz9+PI8++mgKc2SMMbWD3wCgJ/CkqpYD4ZdXPwCt4porE7Pwq38RqWLpzHLNNdcwYsSI4Pvrr7+ed955J3UZMsaYWsBvALANOCDKvIOADfHJjtlXdaUBYCQiwtNPP03fvn0BKC8v59xzz2XlypUpzpkxxqQvvwHAFOBOETkoZJqKSB7wR+B/cc+ZiUldq/8P16BBA/73v//Rrl07ALZs2cLpp58erBYxxhhTkd8A4CbcgD9LcR0CATwNLAd2Arf5SURERonIZyLyo/dYICKnVrH8YBGZIiLrRaTYW/dSn3muU+p6AADQpk0bXnnlFRo0cO1Tv/zySy688ELKy8tTnDNjjEk/vgIAVd0C9MUN+rMamA58iwsM+qvqTz63twa4EdemoDcwE3hVRP4vyvJHA58D5wCHA08BY0XkAp/bqzPqchVAqCOPPJJnn302+H7q1Kn8+c9/TmGOjDEmPfnqCRBAVUuAZ73HPlHVKWGTRovIlUA/4LMIy98TNukpERmCux3x+X3NRyayEoC9hg8fzqeffsr9998PwD333MMRRxzBeeedl+KcGWNM+vDbD0CZ1xdApHm9RCTmIdlEJFtEzsONJPheDKvuB2yJdXuZzkoAKhozZgwnn3xy8P2ll17Kxx9/nMIcGWNMehE/naaISDnQV1U/iDCvL/CuqjaovGbEtI4AFgANge3AcFV9w+e6pwGv4KodKuXFW2YkMBIgPz+/1+TJk/0k7dv27dtp0iS9Rj8uLi7m1FNdU4qcnBzeeuutYOc4+yod9zNW27dv56qrrqKoqAiAAw44gKeffpoWLVpUWKa276cftp+ZxfYzs8R7P4cMGfKRqvaudkFVjfgAOgLHeI9yXAdAx4Q9TgAmAcuipRMh3fq4ngR7A2OATcDhPtbrj2uIeKXfbfXq1UvjbdasWXFPs6Y+//xzxfXPoJ07d45Lmum4n/ti+fLl2qxZs+Dx6d+/v+7atSs4P1P2szq2n5nF9jOzxHs/gUXq4xxZ1WXir4HZuB4AFdcAb3bY4y3gTODOaiONvQFHiaquVNVFqnozsBj4Q1XriMgAYBpwm6o+5XdbdUVo/b8V/1fUpUsXJk+eHCwRmT9/PldffbV1F5xhJk2aRKdOnRg6dCidOnVi0qRJqc6SMWmvqgDgSeAI4Oe4rn+He+9DH12BFqr6nxrmIWr1gYgcgzv536mqj9RgOxnLGgBW7aSTTuK+++4Lvh83bhxPPPFECnNk4kVVefjhh7n00ktZvXp1cEyM3/zmN/ztb39j+/btqc6iMWkragCgqhvVDf+7BDgQeNl7H/r4SlV3+92YiNwrIgNFpJOIHCEiY4DBuGoERGSMiMwIWX4w7uT/NDBJRFp7j2i9EtZJ1gCwetdffz0XXnhh8P3vf/978vPzk3rFGLhKzcrKsqvUfVReXs6SJUt46qmnOP/88+nQoQPXXXcdJSUlFZbbvXs3N954I02bNqVFixZ0796d008/nVGjRnHffffxn//8h/nz51NUVERpaWmK9saY1PJ7G2BjoHNV/cur6lIf6bQGJnrP23C3/p2sqm9789sAB4csPwLIxfU2+MeQ6auBTv6ynvmsBKB6IsLYsWNZvnw5H374IarK999/D7jjd9lll7Fq1SqGDh3K7t272b17N7t27YrL8+7du9m4cSObNm0KVj2sXr2aSy65hBdeeIETTjiBdu3aBR+tW7cmOzs7lYcrbZSWlvLJJ58wZ84c5s6dy9y5c9m8eXNMaWzZsoUtW7bw6aefRpyfnZ1N27Zt6dixY/DRoUOHCu+bN29eaXyNSZMmMXr0aAoLC+nYsSN33303w4cP3+d9NSbZ/AYAS6g8CFC4av+xVHVELPO991WuY6wEwK9GjRrxyiuv0LFjx0q9A+7atYtbb701qfkpKytj6tSpTJ06tcL0rKwsWrduTbt27Wjfvn2F4CD0kYmto3fu3MkHH3wQPOG/99577Nixo8p1RCRim4569eqRlZVVqXQgXFlZGUVFRRQVFTF//vyIyzRp0qRCULBlyxZee+21YNqrV6/m8ssvZ9u2bVx44YXk5uZSr57vblZ8s6DDxJPfb+iQCNNa4O4COAG4Jm45MjGzEgD/2rVrl/YNAMvLy1m3bh3r1q3jww8/jLrcfvvtVykoCA8Ypk+fzq233pq2J4xt27bx3nvvMXfuXObMmcOHH35Y7Qk7Ly+PY445hoEDB3LMMcfwxRdf8Nvf/pbi4uLgMrm5uYwdO5bzzz+f77//nqKiIgoLC4OP0PcbNlQ/ltn27dv58ssv+fLLL6Mus3PnTkaNGsWoUaMAF4A0atSI3NxcGjVqVOl1VfMiLTd//nweeughdu3aBewNOkpLS7n44osTNgJoKoIOC3SSw1c/AFUmIPJXoKOqXhyfLMVP7969ddGiRXFNc/bs2QwePDiuadbEzp07yc3NBVxR5q5du+Jy5ZFu+xlPnTp1qhA0BeTk5NCzZ08aNGhAw4YN4/p8wgknsG7dukrbbNasGeeeey5r164NPjZu3Jiwfc/Ozua4446jf//+tGzZkry8PPLy8oKvW7ZsScOGDRO2/e+//5558+YFr/AXL15c7VgNHTt2rHDC79q1a1yL43ft2sWaNWuqDBKqK4VItaysLHJzc6t9BIIJv8tOnz6d0aNHs3PnzuC2GjVqxKOPPso555wDEAyode8t29VOq2reK6+8wk033VRpm4899hgXXHAB9erVC5buxFMmBToi4qsfgHgEAMcC/1PVZjVKKAHqQgCwfPlyDjnkEMCd2L799tu4pJtu+xlPkyZNYuTIkRGvGBP1g49lm7t372b9+vXBgGDNmjUVAoTAo7or5X3VuHHjSoFBda9Dg4bQP7W2bdsybNgwysvLmTNnDsuWLat2+4ccckjwZD9w4MCYSrUS8b1VVbZs2VIhKLjlllsijjQZOBEXFxfbIFQJJiLBYCAnJyf4OvzhZ9769etZvHgxZWV7O7XNzs6mT58+HHjggRUCzsDrmk5buXIl8+bNq7DNeP0PJTMAeAj4haqmXdlzIgKA7t2707x58wrTfvWrX3HVVVdRXFzMKaecUmmdESNGMGLECDZt2hSMmkNdeeWVnHvuuRQVFXHRRRdVmn/99dczbNgwli9fzhVXXFFh3ubNm/n8888B6NWrV8R64XvuuYejjz6a9957j1tuuaXS/EceeYTu3bszffp0/vrXvwKwdevW4H4+88wzdO3alalTp/Lggw9WWn/ChAl06NCBF154gaeeqtxNw0svvUReXh7jx49n/Pjxlea/+eab5Obm8uSTT/Liiy9Wmj979mwAHnjgAV5//fUK8xo1asS0adMAuOuuu5gxY0aF+S1btuTll18G4Oabb2bBggUAbNiwgW+//Zbdu3dTUFDA3XffzYcffsjixYsrrN+lSxfGjh0LwMiRI1mxYkWF+d27d+eRR9zdqRdeeCFr1qypML9fv36MGTMGcAMVff755+zevZsGDRpw4IEHcsEFFwQHKzr55JMrXPUAnHbaafzxj679a+iJTVUpLS1l0KBBDBgwgG+++YYHHngg2OiwpKQkqVetjRo1Cp7wdu/2fWMQWVlZdO/enYKCApYvX06zZs2oX79+cH6s373Q7y0k7ru3YcMGVqxYUeEkn5WVRZcuXcjPz0dVadGiBePHj2fnzp3ceeedLFq0iPLycsrKyigvL2f//ffniiuuYOfOnUycOJHCwsLgvPLycpo0aULv3r3ZuXMnH374IWvWrEn76itTcwUFBRXade0LvwGAr7JiEan8y3A9+h0C/AyofFYxSRGoDwRo27ZtxKsSU1l+fj75+fnk5OTwzjvvAFRZ3x4PHTt2pHHjxnFJS0TIycmhffv2nHLKKRQXF/PSSy9VWGbhwoURT8b16tWjbdu27Nmzhw4dOtCoUSPWr1/Pt99+S2lp6T6dZMIDl6ry3bRpU5o1a0azZs148sknGThwYIXgszbIz8+nfv36bNu2jcLCQvbbbz/atGlDfn4+4PYzKysruJ/Nmzev9Nm3b9+eCy5wA5t+/vnnlY57ePA5d+7cSkFH4IrxjTfeoLCwsEKA0a1bN379619TXFzMX/7yF7Zt2xacV1ZWRseOHenZsyfFxcW8/vrrlJSUVAhAcnNz2bRpU4Ur1FCB6i0gYrDZsGFDGjZsiKpW+l8SkWA1Q1lZGVu3bg3Oq6pkKysrq0LVQSYqLCxM2rb8jgUwK8LkXbjhfV9R1TfjnbF4qAtVAKNHj+aee9ygibfddht33um7U8Yqpdt+Jkom7+e+VHUE/qw3bdrEDz/8wKZNm3y99nMv/ezZs+nTpw+NGjWK2z5G2kamfp6Q/HrqdK4uKy8vZ8+ePZSWlkZ8VDUvfP5vfvObiG1vWrZsGSzhg8rtGGoy7cYbb4x4S2sySwB89alfWx91YSyA4cOHB/u5f/bZZ+OWbrrtZ9a0xFwAACAASURBVKJk+n5OnDhRCwoKVES0oKBAJ06cGPdtlJeX67Zt2/Trr7/W1q1bB7+PoY+CgoK4bzeSTP88A5K5n8n4DqV6mxMnTtTc3NwK39nc3NyEbjeR28TnWAAxn1Rx3QIfgFd6kM6PuhAA9O/fP/jlmTFjRtzSTbf9TBTbz/hKxR9pKPs8M4sFOvvGbwDg+z4KETlFRN7DFf1/B+wSkfdE5FS/aZj4q419AIwbNw4RqdRgLt2UlZVxyy23UFBQQHZ2Nr17V1+i5seAAQM47rjj4pJWuhk+fDhjx46loKAAEaGgoCChxcXGxMvw4cNZtWoV5eXlrFq1Kinf2cA2Z86cmbRthvIVAIjIFcBUYDuu059fes/bgde8+SbJSkpKWLt2LeAa1XTo0CFuad92223k5ORUeU/6o48+iohU6skuUzz33HOMGTOG0047jfHjx3PXXXelOku1Qir+SI0xsfPbY8wtwFhVvTJs+tMi8jQwGngmrjkz1Qq9Laht27YVbp+qqeOPP565c+cyefJkfve730VcZuLEieTl5XHSSSfFbbvpZPr06eTl5cV95MCZM2cmrNc2Y4zxy28VQEvgf1HmvYzrFtgkWSKL//v27UuLFi2YOHFixPkrVqxg0aJFnHfeeeTk5MR12+ni+++/p1mz+PdvVb9+/Yw9ZsaY2sNvADALGBRl3iBgTnyyY2KRyEGAcnJy+NWvfsUHH3zAV199VWn+hAkTACp0XPTqq68ybNgw2rdvT/369Wnfvj2jRo3y1TdBtHrxW2+9NWLXxpMnT+aoo44iNzeXZs2aMWzYML744gtf+/bDDz/w29/+ljZt2nDCCSdw6KGH8tBDDwVLU1auXImIMGvWLL7++mtEBBGJGgyB6wSoYcOGFBYWMmzYMJo2bUpeXh6jRo2qdI90JrcBMMbUHn4DgMeAi0TkKRE5UUR6eM9PAxcBD4tIt8Ajcdk1oRLdADBwco80bv3zzz9Ply5d6NOnT3Das88+S/369fnd737H448/zrBhwxg3bhzDhg2La77uvfdezj//fNq3b8/999/PLbfcwmeffUb//v1ZuXJllevu2rWLwYMHM27cOM455xyuvPJKDjzwQK6//nr+8Ic/ANC6dWsmTJjAIYccQn5+PhMmTGDChAn079+/yrTLy8s56aSTaNiwIffddx+nnXYaTz75JOedd17c9t0YY+LGz60CQHnYoyzC+8C0Mj9pJuOR6bcBXnLJJcFbrZ555pm4ph3Yz4MPPlg7d+5cYd68efMU0LvuuqvC9B07dlRK57nnnlNAFy5cGJz2j3/8QwEtKioKTuvfv78ee+yxldYfPXq0ZmdnB99/++23mp2drbfddluF5datW6fNmjXTiy++uMr9evjhhxXQcePGBfezvLxczzrrLBURXbZsWXDZQYMG6cEHH1xlegGB/hguvfTSCtNvuOEGBXT69OnV7msipdP3NpFsPzOL7ee+Ic63AQ4JewyN8H5oyGuTBMm4BfDCCy9k5cqVLFy4MDht4sSJiEil1t2BUQlVlW3btrFp0yYGDBgAwEcffRSX/Lz00kuUlZVx3nnnBXuj27RpEzk5ORx11FHMnDmzyvVff/118vLyGDFiRHCaiPCnP/0JVeWNN96oUf6uvfbaiO/DxzAwxphU83UXgKq+m+iMmNglKwC48847mThxIn379qWkpIQXX3yRAQMGcOCBB1ZYdunSpdx4443MnDmzQjeeQIW+vmsiMBhPt26Ra5qquxNi1apV/OxnPyM7O7vC9EB6NRlNUUTo0qVLhWlt2rShadOmNe7a0xhj4i3mgeNFpB5uIKAKVLU4wuImQcrKyigqKgq+T1QA0LlzZ/r168cLL7zAI488wptvvsnmzZsrjVq4detWBg0aROPGjbnrrrvo3Lkzubm5lJSUcOqpp1Y7NGq02+LCByIJpDNt2rSIjQNrOka43Z5njKkr/I4G2AwYA5yF1w1whMWyI0wzCbJu3brgACytWrVK6AArF130/9u783ArqjPf498fiAhCVECBaAMq0auttooxGqMBEwcwakyc0auJynWMQxxi03Hovoo3orY+Dgkxthq4jYaOA4qCojhEHDAYB5xo40UFURRRRhHe+0fVwWJzhn3OqT2cfX6f56mHU6tWVb3v3odTq1atqjqe008/ncmTJzN27Fg6d+7MEUccsVadqVOnsmDBAu699961BsvNmjWrqH1ssskmzJ07d53ywjPngQMHAsmb9RrqBWjMgAEDePnll1m1atVavQCvv/76muUtFRG89dZb7LjjjmvK5s2bxxdffJH7XRpmZq1V7OnS7cDRwJ3AacDP65msjEp5C2Cho446ivXXX5+bbrqJBx54gIMPPnitd67D12fehWf6V199dVH7GDhwIK+99tpaTx6cM2cO999//1r1Dj/8cDp27Mill15ab69CY08uBDj44IP5+OOPufPOO9eURQSjR49GEgcd1LonW2ffHJadb+12zczyVuwlgB8A/ysi/rOUwVjxyvkOgB49ejBs2DDuvfdegHW6/wH23ntvevTowfDhwznrrLPYYIMNmDhxIgsWLChqHyeffDLXX389+++/PyeddBILFy7klltuYbvttuOll15aU2/gwIGMGjWKCy+8kD322IPDDjuMnj17MmfOHB566CF22WUXbr311gb3M2LECG699VZGjBjBzJkzkcTVV1/NpEmTOPvss9l2222b+el8rVOnTkyfPp0jjjiCIUOG8Pzzz3PHHXcwbNgw3/dvZlWn2B6AOYCv8VeRcr8EqO6g37NnT4YOHbrO8l69ejFp0iT69+/P5ZdfzuWXX06fPn2YNGlSUdvffvvtGTt2LIsXL+a8885j/Pjx3HDDDfU+ZviCCy7gvvvuY8MNN2TUqFGcc8453HXXXey0006ccsopje6nS5cuTJs2jZNOOom7776bm2++mdmzZzN69Giuu+66omJtSIcOHXj44YdZvnw5F110ERMnTuTUU0/lrrvuatV2zcxKoph7BYFhwAtAv2LqV8tUy88BOPnkk9c8A+DGG2/MffvVkmep5ZXn8OHDo3PnzrlsqxT8fdYW51lbKvUcgGJvA5wk6YfAbEnvAuvc0xURu6+zopVMW3wNsJmZVY9i7wIYDZxD0gswG/iylEFZ08o5CNDMzGpPsYMATwZGRsSoUgZjxVm9ejVz5sxZM+8eADMza65iBwEuBVr9LFdJZ0h6WdLn6TRdUqP3R0naUdITkpZJ+kDSJWrnT2uZP38+K1asAJL757t3717hiGzs2LEsX7680mGYmRWt2AbA9cCIHA687wMXAbsCuwGPAfdK2qm+ypK+ATwCzAe+DfwCuAA4r5VxtGnZ6//u/jczs5Yo9hJAL+A7wJuSprHuIMCIiIua2khE3FdQNFLSacCewMv1rDIc6AqcEBHLgFclbQecJ+nadLRju+MBgGZm1lrFNgAOB74COgH71bM8SM7siyapI3AE0A14poFqewJPpQf/OpOBfwMGAC1/c0sb5gGAZmbWWsXeBrhl07WKI2lHYDqwAbAYOCwiXmmgeh+SywZZ8zPL1mkASBoBjADo3bs306ZNyyHqry1evDj3bTbXM8983V5auXJlSeKphjzLwXnWFudZW5xniRXzsIA8J5I3CQ4kGQMwClgA7NBA3SnAHwrK+pP0OOzR1L5q9UFAQ4cOXfMQoHvuuack+6iGPMvBedYW51lbnGfLUOSDgIp+d6qkrSTdIumVdDT+K5JulrRVMxscX0bE7IiYEREXAy8B5zZQ/UOSM/2szdJ/59NOeQyAmZm1VlENAEmDSA7UPyV5GNCd6b8/BWZK2rWVMXRuYNl0YG9JG2TK9gPmAu+2Yp9tVkS4AWBmZq1W7CDA0cBMYGhErHkpkKSuwKR0+b5NbUTSVcCDwHtAd+BYYDBwULp8FLB7RPwgXeX/ApcCt0v638A2wK+Ay9Nujnbnk08+YcmSJQB0796dTTbZpMIRmZlZW1RsA2B34MjswR8gIpamjwku9nVnfYCx6b+LSG79GxoRk9PlfYGtM9tfJGk/4CZgBrAQuAa4tsj91ZzCs/92/kwkMzNroWIbAMuAng0s6wEU9Qi0iDixucsjuUNgn2K23x74FkAzM8tDsYMAHwSukvS9bGE6PwqYmHdgVj9f/zczszwU2wNwHnAf8ISkj0lG4G+WTs8AvyxNeFbIDQAzM8tDsQ8C+gT4nqQDSZ7J3xeYBzwXEVNKGJ8V8CUAMzPLQ7E9AABExMPAwyWKxYrgHgAzM8tDg2MAJPWU9F+SDmikzgFpnc0aqmP5cg+AmZnlobFBgOcAW5E8jrchU4At8RiAsli0aBGLFi0CoEuXLmy66aYVjsjMzNqqxhoARwK/beyBO+my3wGH5h2YrSvb/d+vXz8/A8DMzFqssQZAf2BWEdt4neTVvFZi7v43M7O8NNYAWAZ8o4htdEvrWol5AKCZmeWlsQbAX4FDitjGoWldKzH3AJiZWV4aawDcBJwk6YSGKkj6n8DPgBvzDszW5R4AMzPLS4PPAYiIP0u6HvgPSWeS3P8/BwigH3AAsBtwXUTcU45g2zs3AMzMLC+NPggoIn4paRrJLYHnA53TRSuAvwCHRsQDJY3Q1vAlADMzy0uTTwKMiInAREnr8fUbAT+JiK9KGpmtZcmSJSxYsACATp060bdv3wpHZGZmbVnRjwJOD/jzSxiLNaLwGQAdOhT7IkczM7N1+SjSRvj6v5mZ5ckNgDbCDQAzM8uTGwBthAcAmplZntwAaCPcA2BmZnlyA6CNcA+AmZnlyQ2ANsI9AGZmlic3ANqA5cuXM2/ePAA6dOjA5ptvXuGIzMysrXMDoA1477331vy8xRZb0KlTpwpGY2ZmtcANgDbA3f9mZpY3NwDaAA8ANDOzvLkB0Aa4B8DMzPJW1gaApIslvSDpc0kfS5ooaYci1jtA0nRJX0haIOk+SduUI+Zq4B4AMzPLW7l7AAYDNwPfBfYFvgIeldSjoRUkbQncBzwF7AL8EOgCTCp1sNXCPQBmZpa3ot8GmIeIOCA7L+l4YBGwFzCxgdUGAZ2AiyNiVbreKOAxSb0iYkEJQ64KbgCYmVneKj0GoHsaw8JG6swAVgInS+ooqTtwAvBCezj4r1y5kvfff3/NfL9+/SoYjZmZ1QpFROV2Lt0NfAvYre7svoF6ewN/AnqRNBhmAkMj4qN66o4ARgD07t170Pjx43ONefHixXTr1i3XbTbmww8/5JhjjgGgZ8+eTJgwoSz7LXeeleI8a4vzrC3Os2WGDBnyYkTs1mTFiKjIBFwLzAW2aqJeH+At4DckYwD2AaalU4fG1h00aFDk7fHHH899m03tDwgg9txzz7Lutz1wnrXFedYW59kywIwo4jhc1jEAdSRdBxwNDImId5qofgawJCIuzKx/HPAeyWDCp0sWaBXw9X8zMyuFsjcAJF1PcvAfHBFvFLFKV6Dw8kDdfKXHMJScGwBmZlYK5X4OwE3Az4BjgIWS+qRTt0ydUZKmZlZ7ENhV0qWSviVpV+A/SHoAXixn/JXgZwCYmVkplPsM+nSSkf9TgXmZ6fxMnb7A1nUzEfEYcCxwKMngv8kkdwUcGBFLyhN25bgHwMzMSqHczwFQEXVOrKdsPJDvcP42wj0AZmZWCjV/Db0tW7169VqvAvYzAMzMLC9uAFSxefPmsXLlSgB69erFhhtuWOGIzMysVrgBUMXc/W9mZqXiBkAV8wBAMzMrFTcAqph7AMzMrFTcAKhi7gEwM7NScQOgirkBYGZmpeIGQBXzJQAzMysVNwCqVES4B8DMzErGDYAq9dFHH7F8+XIANt54YzbaaKMKR2RmZrXEDYAq5bN/MzMrJTcAqpQbAGZmVkpuAFQpDwA0M7NScgOgSrkHwMzMSskNgCrlHgAzMyslNwCq0Lhx45gyZcqa+VmzZlUwGjMzq0VuAFSZcePGMWLEiDWvAQa48sorGTduXAWjMjOzWuMGQJUZOXIkS5cuXats2bJljBw5skIRmZlZLXIDoMrMmTOnWeVmZmYt4QZAlenXr1+zys3MzFrCDYAqc8UVV9C1a9e1yrp27coVV1xRoYjMzKwWuQFQZYYPH86YMWPo378/kujfvz9jxoxh+PDhlQ7NzMxqyHqVDsDWNXz4cB/wzcyspNwDYGZm1g65AWBmZtYOuQFgZmbWDpW1ASDpYkkvSPpc0seSJkraoYj1JOkcSW9IWiFpnqSryhGzmZlZLSr3IMDBwM3AC4CAfwUelbR9RHzayHrXAD8CLgBeATYC+pY2VDMzs9pV1gZARByQnZd0PLAI2AuYWN86krYFzgJ2iojXM4tmlipOMzOzWlfpMQDd0xgWNlLnUOAd4EBJ70h6V9IdkjYrS4RmZmY1qNLPAbgeeAmY3kidrYD+wNHAiUAAo4GJkvaMiNXZypJGACPS2cWS3kx/3oiktyGrsCw739DPvYAFTSXWhPpiaUm9YnKqr8x5Os+WcJ7Nq+c8Gy9rS3k2tKxa8+xf1BoRUZEJuBaYC2zVRL0xJAf9bTJl26Rl32nG/sY0VZadb+TnGTnkvk4sLalXTE7O03k6T+fpPFtXr6FlbTHP7FSRSwCSrgOOAfaNiHeaqD4P+Coi3sqUvQ18BTTnDTn1jTEoLJtYxM95KHZ7TdUrJqf6ypxnvpxn8+o5z8bLnGe+8sizoWVtMc81lLYcykbS9STd+YNj7UF9DdXfH5gMDIyI/07LtgZmk/QAPF/KeOuJZ0ZE7FbOfVaC86wtzrO2OM/aUqk8y/0cgJuAn5Gc/S+U1CedumXqjJI0NbPao8Bfgdsk7SJpF+A24DlgRhnDrzOmAvusBOdZW5xnbXGetaUieZa1B0BSQzu7PCIuS+vcTtI7MCCzXl/gBuBAYBnwCHBeRMwvZbxmZma1quyXAMzMzKzyKv0cADMzM6sANwDMzMzaITcAzMzM2iE3AHIiaZikNyW9Len0SsdTKpLul7RQ0oRKx1Iqkv5B0jRJsyT9TdJPKh1TqUh6Js3xVUmXVDqeUpLUIX0baS3/7r4r6WVJL0l6vNLxlIqkAZIeS/+PviapV6VjypukHdPvsW5aJunHue7DgwBbT9J6wOvAvsAnJLcn/iAi5lU0sBKQNAToBpwQEYdXOp5SSO866R0RL6XvnHgR2DYillY4tNxJ+kZEfC6pI/A0cFpEvFTpuEpB0hnA3sB6Nfy7+y6wQ0QsrnQspSTpCeDXEfGkpI2A5RGxotJxlYqk7sDfgf4RsSSv7boHIB+7A7Mi4r30IHEPyeuLa05EPA58Uek4Siki5tUdBCPiI5KXVdXcGQZARHye/rh+OtWktCH3E9rPfeU1S9I/Aisj4kmAiFhUywf/1KHAo3ke/MENAAAk7ZN2bX8gKSSdWE+d0yX9XdJySS9K2juz+JvAe5n594HNSxx2s+WQZ5uQZ56SdgM6sfb3WxXyylPSc8BHJH9gqu7sP6c8rwZ+DawuXLda5JRnAE+klzqGlyXwZsohz28BX0i6T9JMSf9atuCbIee/t0cBd+UdoxsAiW7Aq8DZJA8aWouko0jeXHglsAvwDPCQpLp3EaiebVbjtZXW5tlW5JKnpJ7AncBJUZ3XynLJMyK+Q9Jg3VnSDqUOugValaekfYCIiGfKFnHL5PF97hURg4BDgH+WtGPJo26+1ua5HjAY+AVJ7+sgSYeVPuxmy+vv0MbAHsBDuUfY2jcQ1doELAZOLCh7Dvh9QdnbwKj05+8C92WWXQGcUulc8s4zUzYYmFDpHEqZJ9AZeBI4vtI5lPr7zCy7CDi/0rnknSfwK+AD4F3gQ2AJ8IdK51KG7/Pqwm1U29TC73MPkt6qumWnkTxNtuL5lOL7BH4O3FmKuNwD0ARJ6wODgCkFi6aQHPgBngf+Ucno8S7AYcCD5Yuy9YrMs80rJk9JAm4HHouIP5Y1wJwUmefGSkdPS9oA2B94o5xxtlYxeUbEVRGxeSSPFz8aeCgiTiproK1U5Pe5YTpYDCXvV9kXeK2ccbZWkX+HXgB6SuqZ/l/9PjCrfFG2XjP/3h4FjC9FHG4ANK0X0BEofO/AfKAPQER8BZwLTAVeAW6JiLnlDDIHTeYJIOlR4E/AMEnvS9qzfCHmopg89yL5T/djfX0LTjV2pTammDx7AJMlvUxy58oTEfFA+ULMRVG/tzWgmDx7A09L+hvwLMlZ4wvlCzEXxfy9XUXSW/U48DKwALi7jDHmodi/t72AnUnef5O79Uqx0RpVeA1Y2bKImEj+73euhKby/GF5wymZBvOMiKepncZxY3m+Q3IWUgsa/b1dUyliGjCtDPGUSlPf5z+VPaLSaOrv0BRgp7JGVBpN5bmApGFXErXyR66UFgCrWPdsYjPWbb21Zc7TebZFztN5tkVVkacbAE2IiC9JHgSzX8Gi/UhGbdYE5+k82yLn6TzbomrJ05cAWDNgZmA62wHoJ2ln4NOImANcC/xR0vPAX4BTSe79/20l4m0p5+k8cZ5Vy3k6T8qdZ6Vvj6iGieS2tqhnuj1T53SS24hWkLTc9ql03M7TeTpP51npuJ1n283T7wIwMzNrhzwGwMzMrB1yA8DMzKwdcgPAzMysHXIDwMzMrB1yA8DMzKwdcgPAzMysHXIDwMzMrB1yA8BqiqTLJIWkyfUsmyBpWhljGZzGskO59tkckraT9JSkJWmcAyodU5akMyVV9EElko6UdGI95dMkTahASGa5cQPAatX+kr5d6SCq3NXAxsAhwJ7AvMqGU5WOBE6sdBBmpeAGgNWiT0neEz6y0oGUkqQNWrmJ/wE8EhFTI+LZiFhR5H67tHK/7ZY/O6smbgBYLQrgSuAQSTs2VCm9XLCgnvKQdGZm/l1JoyX9StI8SYskXaPEMEmvSfpC0r2SNqlnV9+U9EDa1T5H0qn17PN7kp6QtFTSJ5J+L6l7ZvmJaVy7p93Py4ALGsltZ0lT0+0tlDROUu902YC0a31r4Nx0u9Ma2VZIOk/Sv0v6GHglLT9I0iOSPpL0uaRnJe1f32csaZd0+VJJMyXtXVCvs6QbJX0m6VNJ1wGd6olly/Rz/jz9zCdKGlhQJySdm35Hn6T7Pz9ddoKkd9L93NZYI0rS7cBPge+n2wxJlxXUOVbS7DSehyRtkVk2IF1nuKQ7JX0GTEyXdUw/mzmSVqS/Q8dm1t03XfebmbLpklZJ2jhT9oqkK9KfN5Z0q6S5kpan2/59Q/mZVfyFCZ485TkBl5G8a7sD8AYwPrNsAjCtsG492wjgzMz8u8D7wJ+BA0l6FgK4juQFHj8BhgMLgd9m1huc1nuPpEFyAPC7tOxHmXp7kbwM5C5gGHA88AEwIVPnxHS9/wbOB4YAuzTwGWwKfAZMB34MHJfG/zKwPtAZ2IOky39c+vP2jXymkda9K81/WFp+JvCLNK/9SN5utgrYq+AzXpru+2fAUODZ9Dvqmql3HbAc+GVa589pzJGp0xl4B3gTOIrk4Pxq+ln1KIj3/fSz3h+4IS37DfA48COSl7CsAH7VSN5bA48Bf00/oz2ALdJl09Lv9RngUOBokve4T8qsPyDz2d2Ufkb7psuuAFYC/5J+fmPSuseky7sAXwJHpfNd0/llwEFpWQ9gNXBgOn8bye/8UcD30+99TKX/T3qq3qniAXjylOdE5qBOctBcBWyTzremATAb6Jgpex74CtgyU/YbYH5mfnC6rTEF238EeDYz/xTweEGdfdN1d8jkEsDZRXwGV5E0AL6RKds9e4DJ5DW6iO0FMLOJOh1IXi8+Gbit4DOOugNfWrZzWlZ34OqZHtguKtjeG6zdADg1/cy3ypRtkR4YLy6I9/GCbc0jaaBlP5O7geeayGut35lM+TRgEbBJpuycdN9d0vkB6fw9Bev2AJYAlxaUTwLezMxPB27M/D4sAMYDV6Vlh5D8fn8jnX8VOKtS//c8tb3JlwCslo0F5gAX57CtaRGxKjM/G3g3Iv5eULappPUL1r2nYP7PwKC0G7gryQC8uyWtVzcBT5OcIQ4qWPfBImLdHZgSEZ/XFUTE8yQH/O8VsX591tmvpC0k3SHpA5ID80qSM+5tCqquJDlg1pmV/lvXXb4jsAFwXybe1dn51O7AXyPinUy990nepV6Y19SCbf0deDH7mZB8X5sX5tUML0TEwsx8XV6F2yz87HYgOaP/U0H5XcA2kjZL558C6i6V7JPOP1FQ9rdMTi8BF0g6XVLhd2C2DjcArGZFxFckZ+XHSerfys19VjD/ZQNlIulmz/qonvn1gF7AJkBH4GaSA2XdtILkGvg/FKw7v4hY+zZQbz7J2WdLrLU9SR2A+4HvApeQXJL4NvAQycE86/P0IAxARHyZ/lhXr0/6b32fU1Zz8ir2+2rNQMr6tkc92yyMuW8D5XXzdeNIngR2SK/5703SAHgK2C0du1BXVudM4F6S7+NNSW9LOrrIXKwdcgPAat1tJAeSi+pZtpyCg3UDg/haa7N65r8i6dL9jKSb+FKSA2jhdFvBusXcFz+vnn0C9Ca5Q6IlCvc7ENiFpMv5DxHxRETMILl23Vwfpv/W9zlllSKvcij87OputyzMpXf6b10ufyFpUA4mGX/wJPAasBj4AbArmQZARHwWEb+IiD7APwHPAeMkbZ9PGlZr3ACwmhbJrW2jgZ/z9ZlXnfeB7pKyXbb7k7/D6pl/MSJWRcQSkkFx20bEjHqmuS3Y33PAAQV3EXyb5Jr00y3MoVDdgX7NrYNpL8teLdjWKySNsUMz2+qQnU89R3LpZMtMvc1JeiHyyqtQa3sJ6vMqycDIIwrKjwTeioiPAdLLC68C55Jc658ZEUGS64UkvUj15h0RL5PcJdKB5HZPs3WsV+kAzMrgd8A/kxwonsiUP0wy+Ow2SdcAW5IMNMvb0PRWrSdI7hjYj7UPbhcCUyWtJhl09gXQDzgIGBkRbzVzf9cC6v6FIAAAAY5JREFUpwGTJf0foBvJwMBXgP9qTSIZb5A0oK6R9GugO3A5yYj8ZomITySNAS6X9BXJWe4padxZt5P05Dwk6RKSg+JlJD0pv2tZGk16AzhU0o9J8p3bwkbZGhHxqaR/B/4lzXcGye/FMOCYgupPAmcAkzNjUJ4ieYjT2xFR13uCpKdJxpu8StLrcArJYMPnWxOv1S73AFjNi4ilJLeZFZYvILmVbAuSa6fHAccW1svBySTdtfeS3IJ2RkTcn4njaZIBXZsCfyS5V/xCktvMirnmv5b0DHIIyVn1f5LcgvYUsF/m+nurpD0rPyG5lDEB+DdgFGs3sJrjQpLLHZeQxDyXpCFTuM8fkhyU/wDcAfw/YHBElOoSwM3AlDS2F4AROW33EpLP6zTgAZLv/7iIGF9Qr66L/8l6ygrP/qeT3C0ygeQOh17A0HSgpNk6lPQomZmZWXviHgAzM7N2yA0AMzOzdsgNADMzs3bIDQAzM7N2yA0AMzOzdsgNADMzs3bIDQAzM7N2yA0AMzOzduj/AwCrBpHEOkaNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,5))\n",
    "plt.title(\"Computing pi with \\nincreasing number of random throws\",fontsize=20)\n",
    "plt.semilogx(n_exp, pi_exp,c='k',marker='o',lw=3)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.xlabel(\"Number of random throws\",fontsize=15)\n",
    "plt.ylabel(\"Computed value of pi\",fontsize=15)\n",
    "plt.hlines(y=3.14159,xmin=1,xmax=1e7,linestyle='--')\n",
    "plt.text(x=10,y=3.05,s=\"Value of pi\",fontsize=17)\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### So, effectively we can average a few experiments and take that value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment number 1 done. Computed value: 3.1411944\n",
      "Experiment number 2 done. Computed value: 3.1404488\n",
      "Experiment number 3 done. Computed value: 3.1419368\n",
      "Experiment number 4 done. Computed value: 3.1427144\n",
      "Experiment number 5 done. Computed value: 3.1412816\n",
      "Experiment number 6 done. Computed value: 3.141308\n",
      "Experiment number 7 done. Computed value: 3.1412728\n",
      "Experiment number 8 done. Computed value: 3.1412608\n",
      "Experiment number 9 done. Computed value: 3.1424392\n",
      "Experiment number 10 done. Computed value: 3.1415536\n",
      "Experiment number 11 done. Computed value: 3.1418888\n",
      "Experiment number 12 done. Computed value: 3.1413784\n",
      "Experiment number 13 done. Computed value: 3.141464\n",
      "Experiment number 14 done. Computed value: 3.1418816\n",
      "Experiment number 15 done. Computed value: 3.141988\n",
      "Experiment number 16 done. Computed value: 3.1420944\n",
      "Experiment number 17 done. Computed value: 3.1415448\n",
      "Experiment number 18 done. Computed value: 3.1422472\n",
      "Experiment number 19 done. Computed value: 3.142232\n",
      "Experiment number 20 done. Computed value: 3.14278\n",
      "---------------------------------------------------------------------------\n",
      "Average value from 20 experiments: 3.1417\n"
     ]
    }
   ],
   "source": [
    "n = 5000000\n",
    "sum=0\n",
    "for i in range(20):\n",
    "    p=compute_pi_throwing_dart(n)\n",
    "    sum+=p\n",
    "    print(\"Experiment number {} done. Computed value: {}\".format(i+1,p))\n",
    "print(\"-\"*75)\n",
    "pi_computed = round(sum/20,4)\n",
    "print(\"Average value from 20 experiments:\",pi_computed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error percentage can be easily computed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error percentage:  -0.0034169423615242416\n"
     ]
    }
   ],
   "source": [
    "error_pct = 100*(pi - pi_computed)/pi\n",
    "print(\"Error percentage: \", error_pct)"
   ]
  }
 ],
 "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.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}