{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise sheet\n",
    "\n",
    "Some general remarks about the exercises:\n",
    "* For your convenience functions from the lecture are included below. Feel free to reuse them without copying to the exercise solution box.\n",
    "* For each part of the exercise a solution box has been added, but you may insert additional boxes. Do not hesitate to add Markdown boxes for textual or LaTeX answers (via `Cell > Cell Type > Markdown`). But make sure to replace any part that says `YOUR CODE HERE` or `YOUR ANSWER HERE` and remove the `raise NotImplementedError()`.\n",
    "* Please make your code readable by humans (and not just by the Python interpreter): choose informative function and variable names and use consistent formatting. Feel free to check the [PEP 8 Style Guide for Python](https://www.python.org/dev/peps/pep-0008/) for the widely adopted coding conventions or [this guide for explanation](https://realpython.com/python-pep8/).\n",
    "* Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting `Kernel > Restart & Run All` in the jupyter menu.\n",
    "* Each sheet has 100 points worth of exercises. Note that only the grades of sheets number 2, 4, 6, 8 count towards the course examination. Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are just for feedback.\n",
    "\n",
    "Please fill in your name here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "NAME = \"Kees van Kempen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "42a51f89f2fde456dc33fa324c98aed9",
     "grade": false,
     "grade_id": "cell-f051f1e1b3b3e7be",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__Exercise sheet 1__\n",
    "\n",
    "Code from the lecture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "b99ee0ea4e023a0a82dc338192363b57",
     "grade": false,
     "grade_id": "cell-9b66485da15cd77c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "rng = np.random.default_rng()\n",
    "%matplotlib inline\n",
    "\n",
    "def random_in_square():\n",
    "    \"\"\"Returns a random position in the square [-1,1)x[-1,1).\"\"\"\n",
    "    return rng.uniform(-1,1,2)  \n",
    "\n",
    "def is_in_circle(x):\n",
    "    return np.dot(x,x) < 1\n",
    "\n",
    "def simulate_number_of_hits(N):\n",
    "    \"\"\"Simulates number of hits in case of N trials in the pebble game.\"\"\"\n",
    "    number_hits = 0\n",
    "    for i in range(N):\n",
    "        position = random_in_square()\n",
    "        if is_in_circle(position):\n",
    "            number_hits += 1\n",
    "    return number_hits\n",
    "\n",
    "def random_in_disk():\n",
    "    \"\"\"Returns a uniform point in the unit disk via rejection.\"\"\"\n",
    "    position = random_in_square()\n",
    "    while not is_in_circle(position):\n",
    "        position = random_in_square()\n",
    "    return position\n",
    "\n",
    "def is_in_square(x):\n",
    "    \"\"\"Returns True if x is in the square (-1,1)^2.\"\"\"\n",
    "    return np.abs(x[0]) < 1 and np.abs(x[1]) < 1\n",
    "    \n",
    "def sample_next_position_naively(position,delta):\n",
    "    \"\"\"Keep trying a throw until it ends up in the square.\"\"\"\n",
    "    while True:\n",
    "        next_position = position + delta*random_in_disk()\n",
    "        if is_in_square(next_position):\n",
    "            return next_position\n",
    "    \n",
    "def naive_markov_pebble(start,delta,N):\n",
    "    \"\"\"Simulates the number of hits in the naive Markov-chain version \n",
    "    of the pebble game.\"\"\"\n",
    "    number_hits = 0\n",
    "    position = start\n",
    "    for i in range(N):\n",
    "        position = sample_next_position_naively(position,delta)\n",
    "        if is_in_circle(position):\n",
    "            number_hits += 1\n",
    "    return number_hits\n",
    "\n",
    "def naive_markov_pebble_generator(start,delta,N):\n",
    "    \"\"\"Same as naive_markov_pebble but only yields the positions.\"\"\"\n",
    "    position = start\n",
    "    for i in range(N):\n",
    "        position = sample_next_position_naively(position,delta)\n",
    "        yield position\n",
    "\n",
    "def sample_next_position(position,delta):\n",
    "    \"\"\"Attempt a throw and reject when outside the square.\"\"\"\n",
    "    next_position = position + delta*random_in_disk()\n",
    "    if is_in_square(next_position):\n",
    "        return next_position  # accept!\n",
    "    else:\n",
    "        return position  # reject!\n",
    "\n",
    "def markov_pebble(start,delta,N):\n",
    "    \"\"\"Simulates the number of hits in the proper Markov-chain version of the pebble game.\"\"\"\n",
    "    number_hits = 0\n",
    "    position = start\n",
    "    for i in range(N):\n",
    "        position = sample_next_position(position,delta)\n",
    "        if is_in_circle(position):\n",
    "            number_hits += 1\n",
    "    return number_hits\n",
    "\n",
    "def markov_pebble_generator(start,delta,N):\n",
    "    \"\"\"Same as markov_pebble but only yields the positions.\"\"\"\n",
    "    position = start\n",
    "    for i in range(N):\n",
    "        position = sample_next_position(position,delta)\n",
    "        yield position"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "48539a8a121b164023c052bb11a71ca3",
     "grade": false,
     "grade_id": "cell-99420e260f8510fe",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Empirical convergence rate in the pebble game\n",
    "\n",
    "__(a)__ Write a function `pi_stddev` that estimates the standard deviation $\\sigma$ of the $\\pi$-estimate using $n$ trials by running the direct-sampling pebble game $m$ times. Store this data for $n=2^4,2^5,\\ldots,2^{14}$ and $m=200$ in an array. *Hint*: you may use the NumPy function [`np.std`](https://numpy.org/doc/stable/reference/generated/numpy.std.html). __(12 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "6896cc4d2b06bc6cd9db88bbd93fb488",
     "grade": false,
     "grade_id": "cell-0e3910cb420ef878",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.60000000e+01, 4.26468053e-01],\n",
       "       [3.20000000e+01, 2.84594009e-01],\n",
       "       [6.40000000e+01, 1.98253852e-01],\n",
       "       [1.28000000e+02, 1.49882034e-01],\n",
       "       [2.56000000e+02, 1.01855227e-01],\n",
       "       [5.12000000e+02, 6.41574064e-02],\n",
       "       [1.02400000e+03, 4.84787461e-02],\n",
       "       [2.04800000e+03, 3.37458855e-02]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pi_stddev(n,m):\n",
    "    \"\"\"Estimate the standard deviation in the pi estimate in the case of n trials,\n",
    "    based on m runs of the direct-sampling pebble game.\"\"\"\n",
    "    \n",
    "    return np.std([simulate_number_of_hits(n)/n*4 for i in range(m)])\n",
    "    \n",
    "stddev_data = np.array([[2**k,pi_stddev(2**k,200)] for k in range(4,12)])\n",
    "stddev_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "9eec8b91bf40df9b97932746f8682097",
     "grade": true,
     "grade_id": "cell-90b5022604409d60",
     "locked": true,
     "points": 12,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# If done correctly, your code should pass the following tests\n",
    "from nose.tools import assert_almost_equal\n",
    "assert_almost_equal(pi_stddev(2**3,1000),0.582,delta=0.03)\n",
    "assert_almost_equal(pi_stddev(2**4,1000),0.411,delta=0.03)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "bea3c12c013efeb34a936fe5be28149c",
     "grade": false,
     "grade_id": "cell-7103a31a2719c5a0",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(b)__ Write a function `fit_power_law` that takes an array of $(n,\\sigma)$ pairs and determines best-fit parameters $a,p$ for the curve $\\sigma = a n^p$. This is best done by fitting a straight line to the data on a log-log scale ($\\log\\sigma=\\log a+p\\log n$). *Hint*: use [`curve_fit`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) from SciPy. __(12 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "7ef0a128432a026f4caa85746fd729b6",
     "grade": false,
     "grade_id": "cell-8656545dbb5ec948",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import scipy.optimize\n",
    "\n",
    "def fit_power_law(stddev_data):\n",
    "    \"\"\"Compute the best fit parameters a and p.\"\"\"\n",
    "    \n",
    "    # Define the to be fitted function. We use the log, due to the calculation of the uncertainties.\n",
    "    log_f = lambda n, a, p: np.log(a) + p*np.log(n)\n",
    "    \n",
    "    n, 𝜎 = stddev_data.T\n",
    "    a_fit, p_fit = scipy.optimize.curve_fit(log_f, n, np.log(𝜎))[0]\n",
    "    \n",
    "    return a_fit, p_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "1983b7014d8bac9a28f6fd36fddb3589",
     "grade": true,
     "grade_id": "cell-311260afe5628598",
     "locked": true,
     "points": 12,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "from nose.tools import assert_almost_equal\n",
    "assert_almost_equal(fit_power_law(np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]]))[0],3.36,delta=0.05)\n",
    "assert_almost_equal(fit_power_law(np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]]))[1],-0.31,delta=0.03)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "63cac93ff9a1399db5b3eea45c7c6ff4",
     "grade": false,
     "grade_id": "cell-5c9e6c8197ea6ab5",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(c)__ Plot the data against the fitted curve $\\sigma = a n^p$ in a log-log plot. Don't forget to properly label your axes. *Hint*: use [`loglog`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.loglog.html) from matplotlib. __(12 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "df285cc68ad2e9230aee61ae26cda962",
     "grade": true,
     "grade_id": "cell-1755ac9d1fb8402f",
     "locked": false,
     "points": 12,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAEcCAYAAACYg/MAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA89UlEQVR4nO3dd3hUVf7H8fdJoTfpVZCydELviaCAgIACFhAVFEWavaC77uq6Cr9dUVABESxIExUUYUWliSIdBKRJFSEUaRJaQtr5/TGBHcJMMjMkuZPk83qeeZ7MLWe+95xz73xz5hZjrUVEREREglOI0wGIiIiIiHdK1kRERESCmJI1ERERkSCmZE1EREQkiClZExEREQliStZEREREgpiSNREREZEgpmRNREREJIgpWQsixpiaxpiNxpizxpjHnI7nEmPMFGPMq06Ul1GfbYzZZoxp59T6AX6mX/0hdV05EXMwMcbsN8Z0CHDdDO3zEhyupU+IuPPn+JoR/c7nZM0Yc48xZr0x5pwx5ogx5htjTNtr+XC5ynPAMmttYWvt276upANQ+qy1da21y3xZ1lN9+rN+BgqoP1ySWTGrv0l2oH6acxhj8hpjPjDG/J7yz+tGY0yXVMsUN8Z8aYw5n7LcPb7M8/J56fadrP5O8ClZM8Y8BYwFRgJlgOuBCcBtmRaZH4wxYU7HkEEqA9ucDiKj5KB2cUqm9Ae1iwQ7p/uo058vVwkDDgI3AkWBvwOfGWOquC0zHojHlaP0A941xtT1YZ7fHOkf1to0X7gq5hxwZxrL1AaWAadxfbn0cJu3H3gG+AWIAT4F8qXMex6Ynaqst4C3U/4uD8wBjgO/AY+lKndESrkXcTVmY2AjcBb4POWzXnVbJ73yPMaZMr8S8EXKuieBcemV6U9dAUuBJCAupb7/4mHdEcChlO3bCdwMTAOSgdiU9Z5zq9u9KctuB3r6sa2NgJ9T1v0UmJWqHtMrO3W7pFleqm1M77M91jfp96X9QIf0tiGN+ry8vrc29KVuM7g/pFdX7jF7apc0+y4e+ry3+vGlr6bXd9zifDYlzvPAB7gOsN+krLMYuM5t2RdSyvkT+Igr+7H79qe3rT71USAvcCalbc6lvC61UwcPy6d3TLJAdbf3Uy7N9yFmb3XsbXpA5fnRZwM+nnvro6nmDwa+xvXFewI4DHT0EmNa+7G376OrPj+NbX0AmO/2eXuAz9zeHwQa+lmvGblveD1G+VAPafZZD3E/BHwHvItrP9wF1AEeBw6ktFUvb+sH+kqJvXfK3wVxJWN/cZs/Dfi/tOYF0Hfc+0c0PnyneDgW+dQfrorLhwrpDCSSasdxmx+e0lH/CuQBbkoJoqZbkGtx7ajFgR3A4JR5lYELQJGU96HAEaAlrlG/DcA/UsqtCuwDbnErdxOuL5T8Kcv8ntJBwoFeKQ106cDnS3ne4gwFNgNjUho+H9A2vTIDqKtlwENe1q2J6wBQPuV9FaBa6o7gtvydKdsSAtyNa+cu58O2XqrHJ1PivQNI4MovmPTK9tQuXstzKzfNZdOqb9LoS57qyIdtSF2f+4EOPrSh17rN4P7gSztd3g4P7ZLe/uCxz3urHz/6qtd6dyt7Na4voQrAMVxJVCNcidJS4CW3ZbembFNxYIWn7fdhW33uoynLPwAsdHu/C4hKo408HpNSlvGYrPkQs8c6TmN6QOX502e5huO5pz7q4bMnAKdw7e8hwEvA4jT64X4878fejn1XfH4621oVVyIUApRLaedDKeVUxZW0hPharxm5b6QVd3r1gA991kPcb+P6Z+7mlDafk6qMx4GfPaz335Q69PT6r7fPS1m3DK5/kGqlvG8ExKZa5hlgflrzAug77v3jimXSaj/+dyzyuT9cFVO6C7iGDI+mMT8SOAqEuE37BHjZLch73eb9B5jo9v4n4P6UvzsCe1P+bgEcSPVZLwAfuZX7oNu8KFzZqklV9qt+lOcxTqAVrv8GU/+nl2aZAdTVMrx/OVfHtXN2AMLT61ge1t8E3ObDtkbh+o/VvR5XkvbOmrrs1O3iU3npLetDG3rsS77UkYdt8JasXVN/z8D+kG69cnWy5t4u6dWlxz7vY1167atp1btb2f3c3s8B3nV7/ygw123ZwW7zunpqcx+21a8+D7wJvJHyd0FcI2vFvbSR12NSyntvyVp6MXus4zSmB1ReAH02oOO5pz7q4bN/wm0kF9eXYyDJmrdj3xWf78O2HsQ1CtUHmIQr+amFK5mf5+++kFH7Rnpxp1UP+NBnPcT5A/CM2/t/4ZZs4UrifvFn29Opl3Bco4jvpW6rVMs9jOsY6nVeAH3nwbSW8dZ+/O9YFHB/8OWctZNAyTR+oy0PHLTWJrtN+x1X5n/JUbe/LwCF3N7PBPqm/H1Pyntw/ZdW3hhz+tIL138KZdzWPZgqjkM2pWY8zPelPG9xVgJ+t9YmciVfynTnS115ZK3dAzwBvAwcM8bMMsaU97a8MeZ+Y8wmt7jqASXdFvG2rZ7q8Xc/y06vXa4oz49l06tvb33pKj5sgzfX2t/9LSutOHyt10v82R+89fmrGGP6pVx4dM4Y801afdXHev/D7e9YD+/d69N9m37HVS+ppbet/tZlfWCL299HrLWnPCyX3jEpLWnG7K2O06j7gMrzsk1p9dlrOZ6nVz/1cY2UXFIP189N/kpr/0x97EprW38A2uFKcH7AlRjcmPL6Afw7bmfgvuHrccVTPQTSZxvgGiW7pI6H97+mU4ZPjDEhuH6mjAeGu806BxRJtXgRXCOKac3zl9e68KX9/P0ed+dLsrYK13Dj7V7mHwYqpVTiJdfjys598TnQzhhTEejJ/3bug8Bv1tpibq/C1tqubuu6d6gjQAVjjHGbVsntb1/K8+YgcL2HhNXfMq+prqy1M621bXEd+Czw70uz3JczxlQGJuPqzCWstcVw/VzkXjfeeKrH6/0sO712uR7P0ls2vfr21peu4MM2WE/rpbjW/p5RZflTr5ekPgCnVZfe+nzqcrDWzrDWFkp5dUmZdlVfvcZ+6Y37Pn49rjpNLb1t9bcu3ZO1CLe/U0vvmASuL8kCbu/L+hiz1+OBl+kBl5dKen32Wo7n4GXfM8bcgOscsp1ukxvhGr3wJq392Jd10tvWS8laZMrfP5AqWQPf6jWD942MPq6k7rOXpcSdB9epAJc05Mp2aYCHdjKuO0qc8/L6xsPyhv+dp9fbWpvgNnsXEGaMqeE2LQLX+XppzfPGW9/x1j99bj8f97OrpJusWWtjcJ1nMN4Yc7sxpoAxJtwY08UY8x9gDa7fZp9Lmd4O6I7rBN10WWuP4/qP5CNcO/OOlFlrgTPGmBHGmPzGmFBjTD1jTDMvRa3C9XPEcGNMmDHmNqC523x/y3O3Flcn/j9jTEFjTD5jTJsAygy4rozrnls3GWPy4kqeY1O2F1z/XVV1W7wgrk5wPGXdB3Bl+b5YhescxcdS6rEXV9ajv2WnV54/y6ZZ32n0pdTS24bU9enumvp7BpblT716kl7f9dbnIe36SauvXku/9GaYMaaiMaY4rpGaTwPYVp/r0hhTEijF/0Z0anHlF5W79I5J4PoSuyclps64vujTjdlbHadR9wGV52Gb0uyzGXg8T60BsCXVaFEjXOdVepNmP/VBevvnD0B7XOfXRQPLcZ3jXQLXCfr+1GtG7hvXelxJr8+6i8CtXYwxRXAlIb+kWuaqdrLWdnH7Jy/1q0vq5XFdwFAb6G6tjU1V1nlcF0O9knK8aoPrbhXT0pqXxnb523d8aj8/+sNVfLp1h7X2TeAp4MWUYA7iyiDnWmvjgR5AF1xXfUzAdc6CP8OeM3H9hnt5JMRam4SrgzXEdeXQCeB9XFeneooxHtfJkANxnaB4L66h2IuBlJeq7EvrVsd1dUs0cHeAMQZaV3lxXdlyAtfwdWlcX04Ao4AXjWv49Rlr7XbgDVw73h+4RgJW+PAZ7vU4ANdJsnfj6uiX5vtVdnrl+fnZvtT3VX3Jw+ektw1X1KeHGK+1v19zWf7Uq5f106xLb30+ZXWv9ZPCY1+9ln6ZhpnAQlwnq+/Ddb6Xv9vqT102wHUe1qUvi0PAXcaYFh4+N81jUorHU2I7jev84Lm+xIz344G3ug+0PE/blF6fvebjuQdXjM4YY0rgGoXcmsY66fXTNKW3rdbaXbh+Ylue8v4Mrj64ImV7wfd6zbB9I4OOK2n1WXcRXD2KtsdaewEu/2xZj7RHQNNlXCNXj+DqP0fN/0bg+rktNhTXif/HcJ2jN8Rau82HeZ741Xf8aD+f+oMnxl7x03TOYoxZg+vk0Y+cjkVEMpYxZj+uCzAWOx2Lr3RMkuxGfTY45KjHTRljbjTGlE0Zvu2PK8v/1um4RCR30jFJshv12eCU0+7SXBP4DNdVLXuBO6y1R5wNSURyMR2TJLtRnw1COfpnUBEREZHsLkf9DCoiIiKS0yhZExEREQliOe2cNckiJUuWtFWqVHE6DBGRbGXDhg0nrLWlnI5DshclaxKQKlWqsH79eqfDEBHJVowx6T0WTuQq+hlUREREJIgpWRO/GGO6G2MmxcTEOB2KiIhIrqBkTfxirZ1vrR1UtKivT4kRERGRa6Fz1kQk10hISCA6Opq4uDinQ5EcLl++fFSsWJHw8HCnQ5EcQMmaiOQa0dHRFC5cmCpVqmCMcTocyaGstZw8eZLo6GhuuOEGp8ORHEA/g4pIrhEXF0eJEiWUqEmmMsZQokQJjeBKhlGyJlkr/gL88jnoMWfiECVqkhXUzyQjKVmTrLVxOnzxEEy7Hf7c73Q0Io57+eWXGT16tNf5c+fOZfv27VkYkYgEGyVrkrWaPQS3vgHR62FCK1g9EZKTnI5KJGgpWRMRJWuStUJCXAnb0NVQuQ18OwI+7AzHdzodmUiWee2116hZsyYdOnRg505X3588eTLNmjUjIiKC3r17c+HCBVauXMm8efN49tlnadiwIXv37vW4nIjkbMbq3CEJQNOmTe01P27KWvjlM1fCFn8ebnwO2jwBobrUXTLHjh07qF27NgD/nL+N7YfPZGj5dcoX4aXuddNcZsOGDQwYMIA1a9aQmJhI48aNGTx4MA888AAlSpQA4MUXX6RMmTI8+uijDBgwgG7dunHHHXcAcPLkSY/LSfBx72+XGGM2WGubOhSSZFMaWRO/ZOgTDIyBiLth2Dqo1Q2WvgqT2sHhjddetkiQWr58OT179qRAgQIUKVKEHj16ALB161YiIyOpX78+M2bMYNu2bR7X93U5Eck5dJ818Yu1dj4wv2nTpg9nWKGFSsGdH0H9O+C/T8Hkm6H1cGj3AoTnz7CPEXGX3ghYZvJ0peCAAQOYO3cuERERTJkyhWXLlnlc19flRCTn0MiaBI9at8KwNdCoH6x4C95tA/tXOB2VSIaKioriyy+/JDY2lrNnzzJ//nwAzp49S7ly5UhISGDGjBmXly9cuDBnz569/N7bciKScylZk+CSvxj0eAfu/wqSE2FKV9doW1zGnlsk4pTGjRtz991307BhQ3r37k1kZCQA//rXv2jRogUdO3akVq1al5fv06cPr7/+Oo0aNWLv3r1elxORnEsXGEhAMuQCg/TEn4elr8HqCVCkAnQbA3/plLmfKTmapxO+RTKLLjCQjKKRNQleeQpC55EwcBHkLQQz74QvBsH5k05HJiIikmWUrEnwq9QMHvkRbhwBW+fA+Oaw9Qs9skpERHIFJWuSPYTlhfZ/dSVtxSrB7AdgVj84cwT9lC8iIjmZkjXJXsrUhYGLodOrsHcJjG/Bj7Pe4NnPNnH+YqLT0YmIiGQ4JWuS/YSGQetHYchKKFufG3f+i9u3DmXQW7PZeigDbtYrIiISRJSsSfZVohr0nw/dxtIi7+98cOFR5k38Gx8u36OfRkVEJMdQsibZW0gINH2AsOFrCal6I38NnUajRXfz4qTPOXHuotPRiYiIXDMla5IzFK1Anvs+x/Z6n9p5T/LS4cHMefNRVvx62OnIREREromSNfFLhj7IPaMZg2lwJ/meWE9sje48kvwpJWZ2YsrnX5CQlOx0dCKXvf3229SuXZt+/frRunVrAE6fPs2ECRMuL5P6va9efvllRo8enWGxZoRjx47RuHFjXnjhBXr16kVysvZHEX8oWRO/WGvnW2sHFS1a1OlQvCtYkqL3fszFO2dSLm8c9219kP+OHsiBo8edjkwEgAkTJrBgwQJmzJjBypUrgYxL1oLRunXr6Nu3L6NGjaJ06dKcPKkbW4v4Q8ma5Fh5695K0ac3cKjqXfSM/QLebcNPi750OizJ5QYPHsy+ffvo0aMHY8aMoVChQgA8//zz7N27l4YNG/Lss89e9R5g+vTpNG/enIYNG/LII4+QlJQEwGuvvUbNmjXp0KEDO3fu9Pi5mzdvJioqijp16hASEoIxhpdeeimgbZg9ezYtW7YkIiKCtm3bcvy46x+hnj178uKLLxIZGUnZsmVZvHgx4ErWIiIiAIiJiaFUqVIBfa5IrmWt1Usvv19NmjSx2cmxzQvt4VdqWvtSEbtyTD979vRJp0MSB2zfvt3pEKy11lauXNkeP37cWmttwYIFrbXW/vbbb7Zu3bqXl0n9fvv27bZbt242Pj7eWmvtkCFD7Mcff2zXr19v69WrZ8+fP29jYmJstWrV7Ouvv37F58XGxtqaNWvaNWvWWGutffHFF+0zzzxjk5OTr1iubdu2NiIi4qrXokWLrljuxIkTl/9++eWX7bhx46y11lavXv3yZ8+ZM8cOGDDAWmttnz597PDhw+3gwYPt4sWLA6y17MdTfwPW2yA4huuVvV5hTieLIlmhVIOOJP5lPeunjqD5oRmcGruSEx1fp0rr3k6HJuKTJUuWsGHDBpo1awZAbGwspUuX5tSpU/Ts2ZMCBQoA0KNHj6vWXbx4MY0bN6Z58+YANGjQgG+//RZjzBXLLV++3KdYpkyZwqeffsrFixc5evQoI0eO5MKFC8TExPDkk08CkJiYSLFixQAICwvjnXfeCWi7RQQla5J7hOUrRNNB49mytjf5v3mC6gsfZM/mz6h67zhCCutnGQlu1lr69+/PqFGjrpg+duzYq5Ku1LZu3Ur9+vUvv//5559p3LjxVctFRkZy9uzZq6aPHj2aDh06ADB16lTWrl3L0qVLKVSoEFFRUdStW5dt27bRpEkTQkNDAfjll1+oV68eANOmTfNvY0XkCjpnTXKd+s1vouRTK5l3XX+uP7qIc2OacGbtJ3owvDiqcOHCVyRKqd/ffPPNzJ49m2PHjgFw6tQpfv/9d6Kiovjyyy+JjY3l7NmzzJ8//6qyS5QowS+//ALArl27+OKLL+jTp89Vyy1fvpxNmzZd9bqUqAFs2bKF1q1bU6hQIebMmcPKlSupX78+W7dupWHDhpeX++WXX2jQoME114uIKFmTXKpY4UJ0f+wtvmv7OfuSSlFkwWBOTu4JMYecDk1yqRIlStCmTRvq1avHs88+e9X7OnXq8Oqrr9KpUycaNGhAx44dOXLkCI0bN+buu++mYcOG9O7dm8jIyKvK7tu3L+fOnaNevXoMGjSITz75hBIlSgQUZ//+/Xn77beJjIxk165dVK1alYIFC7Jly5YrkrWtW7deHlkTkWtjrEYTJABNmza169evdzqMDLHz8Gm+n/oK/WOnERIaRugtrxLW7AHX0xEkR9mxYwe1a9d2OgzJJTz1N2PMBmttU4dCkmxK30aS69UsX4z+T73O+NrTWJ9wA2HfPEXsB7fCyb1OhyYiIqJkTQQgf55QnunTmbN3zeFlBpMYvZGk8a1gxduQlOh0eCIikovpalARN53rl6NBpZd5ZkZb7jg6ho6L/k7S1i8IvX08lKnrdHiSkb55Ho5uydgyy9aHLv+XsWWKSK6nkTWRVMoXy8+EId3YHjWR4QmPcebIPuzEKPh+JCRedDo8ERHJZTSyJuJBaIjh8Y5/YV2NJ+k7swmDYydz+w//xm7/CtNjHFRq5nSIcq00AiYi2YRG1kTS0KxKcT594la+q/kKA+Kf5dSpE9gPOsK3f4X4806HJyIiuYCSNfGLMaa7MWZSTEyM06FkmaIFwpnQrzGdbrufTvH/YbbpBKvHw4RWsG+Z0+FJNhQbG8uNN954+UHsGeXBBx+kdOnSXu9vtnPnTho2bHj5VaRIEcaOHQvAmDFjqFu3LvXq1aNv377ExcURFxdH8+bNiYiIoG7dugE/+D2rfPvtt9SsWZPq1avzf//neeS0SpUq1K9fn4YNG9K06f/uoHHw4EHat29P7dq1qVu3Lm+99Vaa0+Pj44mKiiIxURcgSRZw+uGkemXPV3Z7kHtG2XX0jL1lzA/2zudftydG1rX2pSLWfjXc2gt/Oh2a+CBYHuQ+btw4O3bs2Awv94cffrAbNmy44gHw3iQmJtoyZcrY/fv32+joaFulShV74cIFa621d955p/3oo49scnKyPXv2rLXW2vj4eNu8eXO7atWqDI87IyQmJtqqVavavXv32osXL9oGDRrYbdu2XbVc5cqV7fHjx6+afvjwYbthwwZrrbVnzpyxNWrUsNu2bfM63VrXQ+ynT5/uNSY9yF2vjHppZE3EDzXKFGbusDbUatmZ1jH/Yna+3tiN02F8C/j1a6fDk2xixowZ3HbbbQDMnj2bli1bEhERQdu2bTl+/HjA5UZFRVG8eHGfll2yZAnVqlWjcuXKgOvB67GxsSQmJnLhwgXKly+PMYZChQoBkJCQQEJCQrrPIe3Tpw933303LVq0oHLlynz9ddbsF2vXrqV69epUrVqVPHny0KdPH7766iuf1y9Xrtzl56UWLlyY2rVrc+jQIa/TAW6//XZmzJiR8RsjkoqSNRE/5QsP5ZXb6vHOfa14Lb4PvRP+xQlbGGbdA58PgHPHnA5Rglh8fDz79u2jSpUqALRv357Vq1ezefNmOnbsyGeffXbF8pGRkVf8dHnptXjx4muKY9asWfTt2xeAChUq8Mwzz3D99ddTrlw5ihYtSqdOnQBISkqiYcOGlC5dmo4dO9KiRYs0y928eTNVq1ZlzZo1zJgxg3/+859Zsj2HDh2iUqVKl99XrFjxclLlzhhDp06daNKkCZMmTfJY1v79+9m4ceNV25p6er169Vi3bt01xS3iC10NKhKgTnXL0rjydbw8ryQtf6nE34st5L4dnxGybxl0/jc0uAvSGYWQ3OfEiRMUK1bs8vspU6bw6aefcvHiRY4ePcrIkSOvWH758uUZHkN8fDzz5s1j1KhRAPz555989dVX/PbbbxQrVow777yT6dOnc++99xIaGsqmTZs4ffo0PXv2TPOZn7GxsZw4ceLyuW116tThzz//vKbt6dChA0ePHr1q+muvvXZ5dBJcp/Sk5mkUcMWKFZQvX55jx47RsWNHatWqRVRU1OX5586do3fv3owdO5YiRYqkOT00NJQ8efJw9uxZChcu7Nd2ifhDyZrINShZKC/j7mnMwojyvDi3IDPjIvigwFQqfjkItnwO3cZAsUrpFyS5Rv78+YmLiwNg6tSprF27lqVLl1KoUCGioqKoW/fKmy9HRkZy9uzZq8oZPXo0HTp0CCiGb775hsaNG1OmTBkAFi9ezA033ECpUqUA6NWrFytXruTee++9vE6xYsVo164d3377rddkbevWrdSoUYN8+fIB8PPPPxMREXFN2+PriFvFihU5ePDg5ffR0dGUL1/+quUuTStdujQ9e/Zk7dq1l5O1hIQEevfuTb9+/ejVq9fldbxNB7h48eLl7RXJLErWRDJAp7plaVG1BKMWlCZqXTmeLLKMoftnEDqhJXT8JzR5UA+GFwCuu+46kpKSiIuLY8uWLbRu3ZpChQoxZ84cVq5cSf369a9YPjNG1j755JPLP4ECXH/99axevZoLFy6QP39+lixZQtOmTTl+/Djh4eEUK1aM2NhYFi9ezIgRIy6vd/PNNzN16lQqVKgAuH4CPXDgAHFxcSQlJfHSSy/xn//8J9O3B6BZs2bs3r2b3377jQoVKjBr1ixmzpx5xTLnz58nOTmZwoULc/78eRYuXMg//vEPwDUyN3DgQGrXrs1TTz11eR1v0wFOnjxJqVKlCA8Pz5RtErlE3x4iGaRo/nD+r3cDpj3Uis/DunHj+VHsyVMLvn4aptwKJ/Y4HaIEiU6dOvHTTz/Rv39/3n77bSIjI9m1axdVq1alYMGCAZfbt29fWrVqxc6dO6lYsSIffPABAF27duXw4cMAXLhwgUWLFl0xQtSiRQvuuOMOGjduTP369UlOTmbQoEEcOXKE9u3b06BBA5o1a0bHjh3p1q0bAMnJyezZs+eKCxo2b95Mv379aNeuHc2aNWPIkCG0adMm4O3xR1hYGOPGjeOWW26hdu3a3HXXXZdHKS9t/x9//EHbtm2JiIigefPm3HrrrXTu3Blw/Tw6bdo0li5devk8ugULFnidDvD999/TtWvXLNk+yd2Mp9/5RdLTtGlTu379eqfDCFoX4hN5c+EuPlyxjwcLrmSEmUa4jYd2L0Cr4RCqQW0n7Nixg9q1azsdBhs3buTNN99k2rRpTocSsK1bt/Lhhx/y5ptvXp4WFRXF5MmTqVmzpoORZZ1evXoxatQor9vrqb8ZYzZYa5t6XEHEC42siWSCAnnCeLFbHeYMacMPBW6h9dlRbM7XDBa/BO/flPEPEJdspVGjRrRv3z7Db4qblerVq3dFogawd+9eatSo4VBEWSs+Pp7bb7891ySm4iyNrElANLLmu4uJSUz4fi8Tlu2mZ94N/Ct8CnkSYjBtnoCoZyFcJydnlWAZWZPcQSNrklE0siaSyfKGhfJkx78w/9FIdha/iRYxI1lVoD0sHw3vRcKBNU6HKCIiQUzJmkgWqVW2CF8MbcOwrs15MGYgg+1fOX/uLPbDW+CbEXDxnNMhiohIEFKyJpKFQkMMD0dV5dvHozhdIYrmp19lYcHusGYivNsK9i51OsQcT6d+SFZQP5OMpGRNxAFVShZk5kMtebFXc545dy/3JL3M6fgQmNYT5g6D2D/TL0T8li9fPk6ePKkvUslU1lpOnjypm+VKhtEFBhIQXWCQcY7GxPHi3C0s3xHNq9ct4I64OZiCJaHraKjTw+nwcpSEhASio6MvP0FAJLPky5ePihUrXnXDXF1gIIFQsiYBUbKWsay1zP/lCC/P20aluF1MLjaF0ud3Qe0erqStcBmnQxSRDKBkTQKhn0FFgoAxhh4R5Vn81I3cUL81rU++yId57yd513cwvjlsmgn6x0pEJFdSsiYSRIoXzMPYPo2YNKAlk7mdTrGvcjDsepg7BKb3htMHnA5RRESymJI1kSB0U60yLHwyihbNWhF14jnGhj9M0u+rYHxLWDMJkpOdDlFERLKIkjWRIFU4Xziv9azPJ4NaMzfPrUSdG8WuvHXhm2fhoy5wYrfTIYqISBZQsiYS5FpWLcG3T0TR7cYWdD75OP8MHU7CHzvg3Taw/A1ISnA6RBERyURK1kSygXzhobzQpTZzh7VlVeFbaH1mFBvytYAlr8Dk9nBks9MhiohIJlGyJn4xxnQ3xkyKiYlxOpRcqUHFYswb3pb7OzbnvrPDGJzwBGeOR2MntYfFL0OC7h8mIpLT6D5rEhDdZ815J89dZNLyfcxduY1nmMqdoT8QX6wqeXpOgMqtnA5PRDzQfdYkEBpZE8mmShTKywtdavP1iO7sbvVvHkz6K3/8eRY+6kzMnMfh4lmnQxQRkQygkTUJiEbWgs/xsxf56PstlF3/Ovea7zidpzTxncdQtsmtTocmIik0siaB0MiaSA5RqnBenuvRlC7PTmVq7Yn8GR9K2fn3sG7MXRyIjnY6PBERCZCSNZEcplThvAzo04fCT6xiebkBNDy9mAKTWzH1/bf4/eR5p8MTERE/KVkTyaFKX1eMyEfe4ux9C4kvWI77o//Bjrdu55WZSzhw8oLT4YmIiI90zpoEROesZTNJiZxbNpa8P/2bC8nhjEy6FyL6MfzmGlQqXsDp6ERyDZ2zJoHQyJpIbhAaRqGbnyF82CryV6zPv8Peo8eWYdw3+jOen/MLB09ppE1EJFgpWRPJTUpWJ8/Ab+DWN2iVdx8L842g4Kb3uXn0Ul744hei/1TSJiISbJSsieQ2ISHQ7CFChq0hT9VI/h76MYuvG8WmDWtpP3oZL3yxhUOnY52OUkREUihZE8mtilWCfp9Dz0lcn3yIBfle4N1K3/PVhv20e/17/valkjYRkWCgCwwkILrAIIc5dxy+eQ62fUFCyTpMLv4UY7YVBODuZpUY2q465YvldzhIkexPFxhIIDSyJiJQqBTc+RH0mUl43CmG7n6En1v+xD2NS/PpuoO0e30Zf5+7lSMxGmkTEclqGlmTgGhkLQeLPQ2L/g4/T4Xi1Th20+uM3V2az9cfxGDo07wSj91cg5KF8jodqUi2o5E1CYRG1kTkSvmLQY934P6vIDmR0rN7MTL8I5Y92oTeTSoyc80BbhnzI4u2/+F0pCIiuYKSNRHxrGo7GLoKWg6D9R9SYeZNjKp3hAWPR1KmSD4enrqeEbN/4dzFRKcjFRHJ0ZSsiYh3eQpC55EwcBHkLQQz7+QvK55m7gO1GNquGp9vOEiXt35k3f5TTkcqIpJjKVkTkfRVagaP/Ag3Pg9b55BnYkueq7iNzwa1xGC4671V/PvbX4lPTHY6UhGRHEfJmoj4JiwvtH/BlbQVqwSzH6Tp6uF882B1+jSrxLvL9nLb+BXsPHrW6UhFRHIUJWsi4p8ydWHgYuj0KuxdSsHJbRhVeSPv39eE42fj6P7OT7y/fB/JybrSXEQkIyhZExH/hYZB60dhyEooWx/mP0aH9Q+zaEAl2tUsxatf7+Ce91frWaMiIhlAyZqIBK5ENeg/H7qNhcObuG5KO96rvorXe9dl66EzdBm7nDkbotH9HEVEAqdkTUSuTUgINH0Ahq6GqjdiFr7InZseZPG9pahdrghPf76ZIdN/5tT5eKcjFRHJlpSsiUjGKFoB+s6C3h/An/sp+0knZtVcxt9uqcbSX4/RacyPfP/rMaejFBHJdpSsiUjGMQbq3wHD1kLd2wn54f94eMcDfHdXAUoWysMDU9bx1y+3cF430hUR8ZmSNRHJeAVLQu/3XSNtsae5Ye7t/LfWtwxrU55P1h7g1reX8/OBP52OUkQkW1CyJiKZp2YXGLYaGvcnbPV4nt03gAXdLQlJljveXckbC3eSkKQb6YqIpEXJmohkrnxFoftY6P9fwFB7YT++r/Ul90QU452le+g5YQV7julGuiIi3ihZE5GscUOk675srR8lz+bpvHpoIHNuiuHw6ThuffsnPlrxm26kKyLigZI1Eck6eQq4nnzw0GLIX5wmK4ewsvp0Ot8Qxj/nb+e+D9dwJCbW6ShFRIKKkjURyXoVmsCgZdDur+Tb/TVjTwxiZosDbDzwJ7eM+ZGvNh3SjXRFRFIoWRMRZ4TlgXYjYPByTPGqtN78POurvk+LEnE8PmsTwz/ZyOkLupGuiIiSNRFxVunaMHAh3DKKAodWMunsMKbU38rCrYe5ZeyP/LDruNMRiog4SsmaAGCMqWqM+cAYM9vpWCQXCgmFVkNhyEpMhca02z2STZXfoXb4cfp/uJZ/fLWV2Pgkp6MUEXGEkrUcwBjzoTHmmDFma6rpnY0xO40xe4wxz6dVhrV2n7V2YOZGKpKO4jfA/V9Bj3coeGoHH118gsnVVjBj1T5ufXs5mw+edjpCEZEsp2QtZ5gCdHafYIwJBcYDXYA6QF9jTB1jTH1jzH9TvUpnfcgiXhgDje+HYWsw1W6m46HxbCr/H8pf3Evvd1cya+0BpyMUEclSYU4HINfOWvujMaZKqsnNgT3W2n0AxphZwG3W2lFAtywOUcR/RcpBnxmw7UsKL3iWaUnPMbfk3Yz4Ip5df5zjr11rERaq/zdFJOfTkS7nqgAcdHsfnTLNI2NMCWPMRKCRMeYFL8sMMsasN8asP35cJ31LFjAG6vWC4esw9e6g55kZ/FTsZTauXMjAj9dzJi7B6QhFRDKdkrWcy3iY5vXGVdbak9bawdbaaimjb56WmWStbWqtbVqqVKkMC1QkXQWKQ6/3oN9sSudJ4Iu8L9PutzfpO24J+0+cdzo6EZFMpWQt54oGKrm9rwgcdigWkYxRoyMMXYVpNpAHQr/hvXOP8tq4iazcc8LpyEREMo2StZxrHVDDGHODMSYP0AeY53BMItcuXxG49Q0YsIAyxQoxmVeI/vghPl2+xenIREQyhZK1IGGMeeka1v0EWAXUNMZEG2MGWmsTgeHAd8AO4DNr7baMiVYkCFRpQ/iwlVxs+Ri9Q3+g3eJuzPx4AolJyU5HJiKSoYyevxccjDHJwOtAceBnYJa19k9no7qaMaY70L169eoP79692+lwRABIit7IiRkPUyZ2N2vyR1H7wYkUKeX1ehoRxxhjNlhrmzodh2QvGlkLHhaIwzUSVglYaYyJcDakq1lr51trBxUtWtTpUEQuC63YiDLPrGJLzcdodGEldnwL/vhpCuifURHJATSyFiSMMdustXXd3v8FmGitvcnBsLxq2rSpXb9+vdNhiFxly6a1JM8dRgS7OFXuRorfPR6KVUp/RZEsoJE1CYRG1oLHCWNMk0tvrLW7AN0fQ8RP9Rs2p/jwpUzIP4i8h1eT8E5z7NrJkKxz2UQke1KyFjweA6YbY6YbY0YYY2YAvzkdlEh2VKlkYe5/YhSvVPqQVfHVMAueIXlKVzixx+nQRET8pmQtSFhrNwMNgU9SJn0P9HUsIJFsrlDeMEY9eCurWk/mmYRHuHBwC/bd1vDTGEhKdDo8ERGf6Zw1CYjOWZPs5MuN0Yye8yMj837MjUmroVwE9BgH5Ro4HZrkMjpnTQKhkTXxizGmuzFmUkxMjNOhiPisZ6OKvDOoC0+bZ3nCPsXFPw/BpHaw5BVIiHM6PBGRNClZE7/o1h2SXTW+/jrmDW/DruI30TJmJLvKdoXlb8B7kXBgjdPhiYh4pWRNRHKN8sXyM3tIK1rUqU6n3/ryfuXR2PgL8OEtsOA5uHjO6RBFRK6iZE1EcpUCecKY0K8xj91UnVd3luf+fG8R2+hBWDsJJrSCPUucDlFE5ApK1kQk1wkJMTzVqSbv9G3E2sMJdNjRjd9vmw1heWF6L5g7FC6ccjpMERFAyZqI5GLdI8rz2SOtSEhKpuuXiSxp9wVEPg2bZ8GElrB9ntMhiogoWROR3C2iUjHmDW9LtdKFeGjmFt4N7Ycd9D0UKgOf3Qef3gdn/3A6TBHJxZSsiV906w7JicoWzceng1pxa/1y/PvbX3n6R0vcgEVw80uw6zsY3xw2zdSD4UXEEboprgREN8WVnMhayztL9/Dmol00ur4Y793XhNIXD8K8R+HAKqh2E3QbC9dVdjpUyaZ0U1wJhEbWRERSGGN47OYavNuvMb8eOcvt41awLb40DFgAXUfDwbWuK0bXvKcHw4tIllGyJiKSSpf65fh8cCsscMe7q/h2+x/Q/GEYugoqt4JvnoOPusDxXU6HKiK5gJI1EREP6lUoylfD21CzbGEGT/+ZF+du4agpDf1mQ8/34MROmNgGfhwNSQlOhysiOZjOWZOA6Jw1yS3iEpIYuWAHM9ccICTEcE/z6xnSrhplQs7Agmdh+1woWx9uG+96QLxIGnTOmgRCyZoERMma5DYHT11g/Pd7mL0h+nLSNrRdNUofWgRfPw3nT0Cbx+DG5yE8n9PhSpBSsiaBULImAVGyJrnVgZMXGPf9bub8fIiwEEO/FpUZ0rI4pVb+CzZOhxLVocc417ltIqkoWZNAKFmTgChZk9zu95PneWfpHr7ceIjwUMO9LSozvEo0xRY/DacPQLOHocNLkLew06FKEFGyJoFQsiZ+McZ0B7pXr1794d27dzsdjojj9p+4lLRFkycshAebleJR8xn5N0yCohVd92Wr0cHpMCVIKFmTQChZk4BoZE3kSr+dOM87S3Yzd9Mh8oaF8kK9M9zzx+uEndoFEX3hlpFQoLjTYYrDlKxJIJSsSUCUrIl4tvf4Od5Zspt5mw9TODyZCZWW0vrIVEz+61w31q1zGxjjdJjiECVrEgjdZ01EJANVK1WIsX0asfDJG2lXpyL37utAr8TXOEIJ+Lw/fHovnD3qdJgiko0oWRMRyQTVSxfirT6NWPRkFBVrtyDy1N8YndyPhJ2LSB7X3HXlqH7ZEBEfKFkTEclE1UsX5p2+jVjwRHt+q/UQnS6OZENcefhqGAlTboM/9zsdoogEOZ2zJgHROWsigdl59CzvLN5J0R0zeCHsE/KGWhLb/Z38bYdASKjT4Ukm0zlrEgglaxIQJWsi1+bXo2eY+s0KOuz7P24K3cThwvUpfNdECleq53RokomUrEkg9DOoiIgDapUtwsgHulB28Dw+KP0C+c78Rt4PbmTlhyM4c/6C0+GJSBBRsiYi4qA6FYoycOjzHLv/RzYVjKT1gYkceb0ln8ydx9m4BKfDE5EgoJ9BxS96goFI5vp95ecUWTKCIomnmGp6cLHtczzQrjZ5w3Q+W06gn0ElEErWJCA6Z00kE8XFcOrLERTf+Qn7ksvyTqHHuPfue2hS+TqnI5NrpGRNAqGfQUVEgk2+ohTvOxHun0e5IuGMufBXtr//EKPmruVCfKLT0YlIFlOyJiISrKreSP7H1hLffCj9QpfSf+PdvPTGGFbsOeF0ZCKShZSsiYgEszwFyNN1FCEPLaLYdSV5/eKrHPv4fv4560diYnUBgkhuoGRNRCQ7qNiUAsNXkBj5HD3C1jB8Rz/+PXokC7cecToyEclkStZERLKLsDyE3fw3Qgf/SP7SNzAy6U3Mp/3428ffceLcRaejE5FMomRNRCS7KVOXAkO+J6njq7QL38rz+/oz4Y2/M/fnaHSFv0jOo2RNRCQ7CgkltM2jhA9fTWiFRvzDvkfpL+/k+ffncfh0rNPRiUgGUrImIpKdFa9KgYcXkNztLZrmOcDL0Q8xfcxzzFy9j+RkjbKJ5ARK1kREsjtjCGk6gDyPr8NWbcdzZip1FtzBcxNmsf/EeaejE5FrpGRNRCSnKFKeAvd/hu39AbXynWLUiWHMf/tR3l/2K4lJyU5HJyIBUrImfjHGdDfGTIqJiXE6FBHxxBhM/TvI9/gGEmvdzqMhc4hc2pvn3/6IX4+ecTo6EQmAng0qAdGzQUWyB7vzW+LmPk6e2GNMSerMhTYv8EiHeuQJ0//qTtCzQSUQ2ltFRHIwU7Mz+R9fR0LE/QwMXUCPlb3525jxbDp42unQRMRHStZERHK6fEXI1/MtGPA1pYrk5/XzL7Jj0gBGf7WW2Pgkp6MTkXQoWRMRyS2qtKXAY6uJbzGcu0N/4N6f7+KVN0azau9JpyMTkTQoWRMRyU3yFCBPl9cIeXgJha8rzaiLIzkx5R5e/exHzsTpwfAiwUjJmohIblShMQWH/0RC1F/pEraBodv68sbr/2LJ9qNORyYiqShZExHJrcLyEH7TCMKG/ETe0jX4Z9Jb8EkfXpr2HSf1YHiRoKFkTUQktytdi4JDlpDYaRRR4Tt4dk9/Jr7xN+ZvinY6MhFByZqIiACEhBLWeijhw1djKjbhb3YyDZbcCyf3Oh2ZSK6nZE1ERP6n+A0UfOi/JHd/h+vj98K7rWHFW5CU6HRkIrmWkjUREbmSMYQ0uR8zbC1U7wCL/gHv3wxHtzgdmUiupGRNREQ8K1IO7p4Od06BM4dgUjtY+iok6uIDkaykZE1ERLwzBur2hGFrof6d8OPrMDESDq51OjKRXEPJmoiIpK9Aceg5EfrNhvjz8EEn+OZ5198ikqmUrImIiO9qdIRhq6HZQ7DmXZjQEvZ+73RUIjmakjXxizGmuzFmUkxMjNOhiIhT8haGW0fDA99ASDhMux2+Ggaxp52OTCRHUrImfrHWzrfWDipatKjToYiI0yq3hiEroO2TsOkTGN8CdvzX6ahEchwlayIiErjw/NDhZXh4KRQqBZ/2g8/6w7ljTkcmkmMoWRMRkWtXviE8/D3c9HfYuQDGNXONtlnrdGQi2Z6SNRERyRih4RD1DAxeAaVqwtzBMOMOOH3Q6chEsjUlayIikrFK/QUe+Ba6/Ad+X+W6YnTtZEhOdjoykWxJyZqIiGS8kBBo8QgMXQUVm8GCZ2BKVzix2+nIRLIdJWsiIpJ5rqsM930Jt02AY9vh3Taw/E1ISnA6MpFsQ8maiIhkLmOgUT8Ytg7+0gmW/BMm3wRHNjsdmUi2oGRNRESyRuEyrgfD3zUVzh6FSe1hySuQEOd0ZCJBTcmaiIhkrTq3wbA1ENEHlr8BE9vCgdVORyUStJSsiYhI1itQHG6fAPd+AYkX4cPOsOA5uHjO6chEgo6SNRERcU71m11XjDYfBGsnwYRWsGeJ01GJBBUlayIi4qy8haDrf+DBbyEsL0zvBXOHwoVTTkcmEhSUrImISHC4viUM/gkin4bNs1wPht/+ldNRiThOyZqIiASP8Hxw8z9g0DIoXBY+ux8+vdd19ahILqVkTUREgk+5BvDwUujwMuxaCOObw8YZejC85EpK1kREJDiFhkPbJ2HICihdB74aCtN6wp+/Ox2ZSJZSsiYiIsGtZA0YsAC6joboda4rRte8B8lJTkcmkiWUrImISPALCYHmD7tu81G5FXzzHHzUBY7vdDoykUynZE1ERLKPYtdDv9nQ8z04scv19IMfX9eD4SVHU7ImIiLZizGuR1UNWws1u8LSV13PGT28yenIRDKFkjXxizGmuzFmUkxMjNOhiEhuV6g03PWx6+Hw54/B5Jtg0UuQEOt0ZCIZSsma+MVaO99aO6ho0aJOhyIi4lK7u+vB8A3vgRVjXT+N/r7S6ahEMoySNRERyf7yXwe3jYP75kJSvOvig6+fhrgzTkcmcs2UrImISM5RrT0MXQ0th8K6D1y3+di9yOmoRK6JkjUREclZ8hSEzqNg4ELX3zPugC8e0YPhJdtSsiYiIjlTpeYweDlEPQdbZ8O4ZrD1Cz2ySrIdJWsiIpJzheWFm/4Gg36AohVh9gOuB8OfOeJ0ZCI+U7ImIiI5X9l68NAS6PgK7FkM41vAz1M1yibZgpI1ERHJHULDoM3jMGSlK3mb9yhMvQ1O/eZ0ZCJpUrImIiK5S4lq0P+/cOubcOhneLc1rJqgB8NL0FKyJiIiuU9ICDQbCMNWQ5VI+O4F+KATHNvhdGQiV1GyJiIiuVfRinDPp9DrfTi1DyZGwg//gcR4pyMTuUzJmoiI5G7GQIM7Yfg6qNMDvn8NJrWDQxucjkwEULImIiLiUrAk3PEh9PkEYk/B+x1g4YsQf8HpyCSXU7ImIiLirlZX14PhG98PK9+BiW3gt+VORyW5mJI1ERGR1PIVhe5vQf/5rnuxfdwN5j8BcTFORya5kJI1ERERb26Ict2XrdVw+PljGN9SD4aXLKdkTUREJC15CsAtr8HAxZC/GJw55HREksuEOR2AiIhItlCxiesZo6HhTkciuYySNREREV+F5XE6AsmF9DOoiIiISBBTsiYiIiISxJSsiYiIiAQxJWsiIiIiQUzJmoiIiEgQU7ImIiIiEsSUrImIiIgEMWOtdToGyYaMMceB3z3MKgqkfnhe6mklgROZFFp6PMWXFeX4unx6y6U139s8X9oEnGsXp9rEn3Uyul18bSvtK4EvF6z7SmVrbakA15Xcylqrl14Z9gImpTcNWB9M8WVFOb4un95yac33Ns+XNnGyXZxqEyfbxde20r6SdW3iT1s52S565c6XfgaVjDbfx2lOyahY/C3H1+XTWy6t+d7mqU2ufZ2Mbhd/2sop2ld8+xyRTKefQSXLGWPWW2ubOh2HXEntEnzUJsFJ7SJZTSNr4oRJTgcgHqldgo/aJDipXSRLaWRNREREJIhpZE1EREQkiClZExEREQliStZEREREgpiSNXGcMaaqMeYDY8xsp2OR/zHG3G6MmWyM+coY08npeASMMbWNMRONMbONMUOcjkdcjDEFjTEbjDHdnI5FciYla5IpjDEfGmOOGWO2ppre2Riz0xizxxjzPIC1dp+1dqAzkeYufrbLXGvtw8AA4G4Hws0V/GyTHdbawcBdgG4dkUn8aZMUI4DPsjZKyU2UrElmmQJ0dp9gjAkFxgNdgDpAX2NMnawPLVebgv/t8mLKfMkcU/CjTYwxPYCfgCVZG2auMgUf28QY0wHYDvyR1UFK7qFkTTKFtfZH4FSqyc2BPSkjafHALOC2LA8uF/OnXYzLv4FvrLU/Z3WsuYW/+4q1dp61tjXQL2sjzT38bJP2QEvgHuBhY4y+VyXDhTkdgOQqFYCDbu+jgRbGmBLAa0AjY8wL1tpRjkSXe3lsF+BRoANQ1BhT3Vo70Yngcilv+0o7oBeQF1iQ9WHlah7bxFo7HMAYMwA4Ya1NdiA2yeGUrElWMh6mWWvtSWBwVgcjl3lrl7eBt7M6GAG8t8kyYFnWhiIpPLbJ5T+snZJ1oUhuo+FayUrRQCW39xWBww7FIv+jdgk+apPgozYRxyhZk6y0DqhhjLnBGJMH6APMczgmUbsEI7VJ8FGbiGOUrEmmMMZ8AqwCahpjoo0xA621icBw4DtgB/CZtXabk3HmNmqX4KM2CT5qEwk2epC7iIiISBDTyJqIiIhIEFOyJiIiIhLElKyJiIiIBDElayIiIiJBTMmaiIiISBBTsiYiIiISxJSsiYiIiAQxJWsiIiIiQUzJmojkGMaYL40xrxpjlhtjjhpjOjgdk4jItVKyJiI5ST3gtLU2EhgK9HM4HhGRa6ZkTURyBGNMAaAoMCZlUhhw2rGAREQyiJI1Eckp6gIbrLVJKe8bAFsdjEdEJEMoWRORnKIesMntfQPgF2dCERHJOErWRCSnqM+VyVo9NLImIjmAsdY6HYOIiIiIeKGRNREREZEgpmRNREREJIgpWRMREREJYkrWRERERIKYkjURERGRIKZkTURERCSIKVkTERERCWJK1kRERESC2P8DUvcZ5B6G/WIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a_fit, p_fit = fit_power_law(stddev_data)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Convergence of standard deviation of direct-sampled $\\pi$ guesses over $n$ throws among $m = 200$ trials\")\n",
    "plt.loglog(*stddev_data.T, label=\"data\")\n",
    "n_range = np.logspace(1, 4, 10)\n",
    "plt.loglog(n_range, a_fit*n_range**p_fit,\n",
    "           label=\"fitted $\\sigma = an^p$\\n($a = {:.3f}, p = {:.3f})$\".format(a_fit, p_fit))\n",
    "plt.xlabel(\"$n$\")\n",
    "plt.ylabel(\"$\\sigma$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "efacdbc9ad60113c9fd39e149f09e955",
     "grade": false,
     "grade_id": "cell-8be1d8340bd0c065",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Volume of a unit ball in other dimensions\n",
    "\n",
    "In this exercise you will adapt the direct-sampling Monte Carlo method of the pebble game to higher dimensions to\n",
    "estimate the $d$-dimensional volume of the $d$-dimensional ball of radius $1$.\n",
    "\n",
    "__(a)__ Adapt the direct-sampling Monte Carlo method to the pebble game in a $d$-dimensional hypercube $(-1,1)^d$ with an inscribed $d$-dimensional unit ball. __(12 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "71be71d20de0f6ee438696d00c9306e8",
     "grade": false,
     "grade_id": "cell-120de066cebda898",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def number_of_hits_in_d_dimensions(N,d):\n",
    "    \"\"\"Simulates number of hits in case of N trials in the d-dimensional direct-sampling pebble game.\"\"\"\n",
    "\n",
    "    def random_in_d_dimensional_box():\n",
    "        \"\"\"Returns a random position in the box [-1,1)^d.\"\"\"\n",
    "        return rng.uniform(-1,1,d)\n",
    "        \n",
    "    number_hits = 0\n",
    "    for i in range(N):\n",
    "        position = random_in_d_dimensional_box()\n",
    "        if is_in_circle(position):\n",
    "            number_hits += 1\n",
    "    return number_hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "521b242a339c4a797d346f192ac1fe2d",
     "grade": true,
     "grade_id": "cell-e8e97677ad30041e",
     "locked": true,
     "points": 12,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "from nose.tools import assert_almost_equal\n",
    "assert_almost_equal(number_of_hits_in_d_dimensions(100,1),100,delta=1)\n",
    "assert_almost_equal(number_of_hits_in_d_dimensions(2000,3),1045,delta=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "fb64e8b2afc2e451908b8743d83a7ea6",
     "grade": false,
     "grade_id": "cell-2b9ff58ee25d374f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(b)__ Compare your estimates for the volume $V_d$ of the $d$-dimensional unit ball for $N=10000$ trials and $d=1,2,\\ldots,7$ to the [exact formula](https://en.wikipedia.org/wiki/Volume_of_an_n-ball) $V_d = \\frac{\\pi^{d/2}}{\\Gamma(\\tfrac{d}{2}+1)}$ in a plot. *Hint*: the Gamma function is available in `scipy.special.gamma`. __(12 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "68784cc26ebb01389c1ad1675f0e5ef8",
     "grade": true,
     "grade_id": "cell-82ea663199a69ee2",
     "locked": false,
     "points": 12,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "0eff565d5e6284f83691c38b751d3ce8",
     "grade": false,
     "grade_id": "cell-beac1a536ba9895a",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Efficiency of the Metropolis algorithm and the 1/2-rule\n",
    "\n",
    "In the lecture we mentioned the \"1/2 rule\" for Metropolis algorithms: if moves are rarely rejected then typically you do not explore the domain efficiently, while if most moves are rejected you are wasting efforts proposing moves. A good rule of thumb then is to aim for a rejection rate of $1/2$. In this exercise you are asked to test whether this rule of thumb makes any sense for the Markov-chain pebble game by varying the throwing range $\\delta$. \n",
    "\n",
    "__(a)__ Estimate the mean square deviation $\\mathbb{E}[(\\tfrac{4\\text{hits}}{\\text{trials}} - \\pi)^2 ]$ from $\\pi$ for different values of $\\delta$ ranging between $0$ and $3$, but fixed number of trials (say, 2000). For a decent estimate of the mean you will need at least $100$ repetitions. __(14 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "836a5749c3515998c0bf7a7b524886b4",
     "grade": true,
     "grade_id": "cell-e3762269c892c3df",
     "locked": false,
     "points": 14,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "1819e4933a0845a4463f559ed8613bb2",
     "grade": false,
     "grade_id": "cell-2b9b2164827f2cb9",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(b)__ Measure the rejection rate for the same simulations as in (a). __(14 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "ec4f5855aa1a96d41a81fe64542ff847",
     "grade": true,
     "grade_id": "cell-9726c572b9859251",
     "locked": false,
     "points": 14,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "bfc0b7ec0fcf893a10443ba881f99ff5",
     "grade": false,
     "grade_id": "cell-6b4d8983ff5d3d37",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(c)__ Plot both the mean square deviation and the rejection rate as function of $\\delta$. How well does the 1/2 rule apply in this situation? __(12 pts)__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "9ef86bd5fb28c1a681ea44586c49477a",
     "grade": true,
     "grade_id": "cell-0677a8cb2d72d30f",
     "locked": false,
     "points": 12,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  }
 ],
 "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.9.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}