{
 "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.10028581e-01],\n",
       "       [3.20000000e+01, 3.00020182e-01],\n",
       "       [6.40000000e+01, 1.96900792e-01],\n",
       "       [1.28000000e+02, 1.47571487e-01],\n",
       "       [2.56000000e+02, 1.03273922e-01],\n",
       "       [5.12000000e+02, 6.88982665e-02],\n",
       "       [1.02400000e+03, 4.86827171e-02],\n",
       "       [2.04800000e+03, 3.75673313e-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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8DElEQVR4nO3dd3gU1f7H8fdJI3SkWQDp0iH0mgAKCkgHBUQRBZGmohex/LwXG+JVREVEBRFUUEFUFCleadKkKiCCIKBcQJDeA2nn98cG7hJ2k90lyWySz+t59nkyOzNnv3POmdlvzszOGGstIiIiIhKcQpwOQERERES8U7ImIiIiEsSUrImIiIgEMSVrIiIiIkFMyZqIiIhIEFOyJiIiIhLElKyJiIiIBDElayIiIiJBTMlaEDHGVDLG/GyMOW2MedjpeC4yxkw1xrzoRHnp9dnGmF+NMS2cWj/Az/SrP6SsKydiDibGmD+NMa0CXDdd+7wEh6vpEyLu/Dm+pke/8zlZM8bcZYxZb4w5Y4w5YIyZb4xpdjUfLlcYASy11ua31o7zdSUdgNJmra1mrV3qy7Ke6tOf9dNRQP3hooyKWf1NsgL10+zDGJPLGDPZGLMn+Z/Xn40xbVMsU9gY85Ux5mzycnf5Ms/L56XZdzL7O8GnZM0Y8xjwBvAScC1wIzAB6JRhkfnBGBPmdAzppDTwq9NBpJds1C5OyZD+oHaRYOd0H3X68+UKYcBeoDlQEPgnMNMYU8ZtmbeBOFw5Sm/gHWNMNR/m+c2R/mGtTfWFq2LOAHekskwVYClwAteXS0e3eX8Cw4HNwElgBhCZPO9JYFaKst4ExiX/fQPwBXAY+AN4OEW5TySXewFXY9YBfgZOA58nf9aLbuukVZ7HOJPnlwK+TF73KDA+rTL9qStgMZAInE+u75s8rPsEsD95+7YDtwAfA0lAbPJ6I9zqdlfysluBLn5sa23gp+R1ZwCfpajHtMpO2S6plpdiG9P6bI/1Tdp96U+gVVrbkEp9XlrfWxv6Urfp3B/Sqiv3mD21S6p9Fw993lv9+NJX0+o7bnE+nhznWWAyrgPs/OR1FgLXuC37VHI5x4EpXN6P3bc/rW31qY8CuYBTyW1zJvl1sZ1aeVg+rWOSBSq4TU+9ON+HmL3Vsbf3AyrPjz4b8PHcWx9NMX8gMBfXF+8R4C+gtZcYU9uPvX0fXfH5qWzrfcAct8/bCcx0m94LRPlZr+m5b3g9RvlQD6n2WQ9x9we+A97BtR/uAKoCjwD/TW6rrt7WD/SVHHu35L/z4krGbnKb/zHwcmrzAug77v1jHz58p3g4FvnUH66Iy4cKaQMkkGLHcZsfntxRnwYigJuTg6jkFuRaXDtqYWAbMDB5XmngHFAgeToUOAA0wjXqtwH4V3K55YDdwG1u5W7E9YWSO3mZPckdJBzomtxAFw98vpTnLc5QYBPwenLDRwLN0iozgLpaCvT3sm4lXAeAG5KnywDlU3YEt+XvSN6WEKAHrp37eh+29WI9Ppocb3cgnsu/YNIq21O7eC3PrdxUl02tvkmlL3mqIx+2IWV9/gm08qENvdZtOvcHX9rp0nZ4aJe09gePfd5b/fjRV73Wu1vZq3F9CZUADuFKomrjSpQWAyPdlt2SvE2FgZWett+HbfW5jyYvfx/wH7fpHUBMKm3k8ZiUvIzHZM2HmD3WcSrvB1SeP32Wqziee+qjHj57AnAM1/4eAowEFqbSD//E837s7dh32eensa3lcCVCIcD1ye28P7mccriSlhBf6zU9943U4k6rHvChz3qIexyuf+ZuSW7zL1KU8Qjwk4f1vk2uQ0+vb719XvK61+L6B6ly8nRtIDbFMsOBOanNC6DvuPePy5ZJrf3437HI5/5wRUxpLuAaMjyYyvxo4CAQ4vbep8CzbkHe7TbvFeBdt+kVQJ/kv1sDu5L/bgj8N8VnPQVMcSv3frd5MbiyVZOi7Bf9KM9jnEBjXP8NpvxPL9UyA6irpXj/cq6Aa+dsBYSn1bE8rL8R6OTDtsbg+o/VvR5XkfrOmrLslO3iU3lpLetDG3rsS77UkYdt8JasXVV/T8f+kGa9cmWy5t4uadWlxz7vY1167aup1btb2b3dpr8A3nGbfgiY7bbsQLd57Ty1uQ/b6lefB8YCryX/nRfXyFphL23k9ZiUPO0tWUsrZo91nMr7AZUXQJ8N6HjuqY96+OwVuI3k4vpyDCRZ83bsu+zzfdjWvbhGoXoCE3ElP5VxJfPf+LsvpNe+kVbcqdUDPvRZD3H+AAx3m34Bt2QLVxK32Z9tT6NewnGNIr6Xsq1SLPcArmOo13kB9J37U1vGW/vxv2NRwP3Bl2vWjgJFUzlHewOw11qb5PbeHlyZ/0UH3f4+B+Rzm/4E6JX8913J0+D6L+0GY8yJiy9c/ylc67bu3hRx7LfJNeNhvi/leYuzFLDHWpvA5Xwp050vdeWRtXYnMAx4FjhkjPnMGHODt+WNMX2MMRvd4qoOFHVbxNu2eqrHPX6WnVa7XFaeH8umVd/e+tIVfNgGb662v/tbVmpx+FqvF/mzP3jr81cwxvRO/uHRGWPM/NT6qo/1/rfb37Eept3r032b9uCql5TS2lZ/67IG8Ivb3westcc8LJfWMSk1qcbsrY5TqfuAyvOyTan12as5nqdVPzVwjZRcVB3X6SZ/pbZ/pjx2pbatPwAtcCU4P+BKDJonv34A/47b6bhv+Hpc8VQPgfTZmrhGyS6q6mH6tzTK8IkxJgTXaco4YKjbrDNAgRSLF8A1opjaPH95rQtf2s/f73F3viRrP+IabuzsZf5fQKnkSrzoRlzZuS8+B1oYY0oCXfjfzr0X+MNaW8jtld9a285tXfcOdQAoYYwxbu+Vcvvbl/K82Qvc6CFh9bfMq6ora+0n1tpmuA58Fvj3xVnuyxljSgOTcHXmItbaQrhOF7nXjTee6vFGP8tOq11uxLO0lk2rvr31pcv4sA3W03rJrra/p1dZ/tTrRSkPwKnVpbc+n7IcrLXTrbX5kl9tk9+7oq9eZb/0xn0fvxFXnaaU1rb6W5fuyVott79TSuuYBK4vyTxu09f5GLPX44GX9wMuL4W0+uzVHM/By75njCmL6xqy7W5v18Y1euFNavuxL+ukta0Xk7Xo5L9/IEWyBr7VazrvG+l9XEnZZy9JjjsC16UAF0VxebvUxEM7GdcdJc54ec33sLzhf9fpdbPWxrvN3gGEGWMqur1XC9f1eqnN88Zb3/HWP31uPx/3syukmaxZa0/ius7gbWNMZ2NMHmNMuDGmrTHmFWANrnOzI5LfbwF0wHWBbpqstYdx/UcyBdfOvC151lrglDHmCWNMbmNMqDGmujGmvpeifsR1OmKoMSbMGNMJaOA239/y3K3F1YlfNsbkNcZEGmOaBlBmwHVlXPfcutkYkwtX8hybvL3g+u+qnNvieXF1gsPJ696HK8v3xY+4rlF8OLkeu3J5Pfpbdlrl+bNsqvWdSl9KKa1tSFmf7q6qv6djWf7Uqydp9V1vfR5Sr5/U+urV9EtvhhhjShpjCuMaqZkRwLb6XJfGmKJAMf43olOZy7+o3KV1TALXl9hdyTG1wfVFn2bM3uo4lboPqDwP25Rqn03H43lKNYFfUowW1cZ1XaU3qfZTH6S1f/4AtMR1fd0+YDmua7yL4LpA3596Tc9942qPK2n1WXe1cGsXY0wBXEnI5hTLXNFO1tq2bv/kpXy1Tbk8rh8wVAE6WGtjU5R1FtePoZ5PPl41xXW3io9Tm5fKdvnbd3xqPz/6wxV8unWHtXYs8BjwTHIwe3FlkLOttXFAR6Atrl99TMB1zYI/w56f4DqHe2kkxFqbiKuDReH65dAR4H1cv071FGMcrosh++G6QPFuXEOxFwIpL0XZF9etgOvXLfuAHgHGGGhd5cL1y5YjuIavi+P6cgIYDTxjXMOvw621W4HXcO14f+MaCVjpw2e412NfXBfJ9sDV0S/O96vstMrz87N9qe8r+pKHz0lrGy6rTw8xXm1/v+qy/KlXL+unWpfe+nzy6l7rJ5nHvno1/TIVnwD/wXWx+m5c13v5u63+1GVNXNdhXfyy2A/caYxp6OFzUz0mJXskObYTuK4Pnu1LzHg/Hnir+0DL87RNafXZqz6ee3DZ6IwxpgiuUcgtqayTVj9NVVrbaq3dgesU2/Lk6VO4+uDK5O0F3+s13faNdDqupNZn3dXiylG0ndbac3DptGV1Uh8BTZNxjVw9iKv/HDT/G4Hr7bbYYFwX/h/CdY3eIGvtrz7M88SvvuNH+/nUHzwx9rJT09mLMWYNrotHpzgdi4ikL2PMn7h+gLHQ6Vh8pWOSZDXqs8EhWz1uyhjT3BhzXfLw7b24svwFTsclIjmTjkmS1ajPBqfsdpfmSsBMXL9q2QV0t9YecDYkEcnBdEySrEZ9Nghl69OgIiIiIlldtjoNKiIiIpLdKFkTERERCWLZ7Zo1ySRFixa1ZcqUcToMEZEsZcOGDUestcWcjkOyFiVr4hdjTAegQ4UKFVi/fr3T4YiIZCnGmLQeCydyBZ0GFb9Ya+dYawcULOjrvSxFRETkaihZExEREQliStZEREREgpiuWRORHCE+Pp59+/Zx/vx5p0ORHCAyMpKSJUsSHh7udCiSDShZE5EcYd++feTPn58yZcpgjHE6HMnGrLUcPXqUffv2UbZsWafDkWxAp0FFJEc4f/48RYoUUaImGc4YQ5EiRTSKK+lGyZpkrrhzsPlz0GPOxAFK1CSzqK9JelKyJpnr52nwZX/4uDMc/9PpaEQc9eyzzzJmzBiv82fPns3WrVszMSIRCUZK1iRz1e8Pt78G+9bDhMaw+l1ISnQ6KpGgpGRNREDJmmS2kBBXwjZ4NZRuCguegA/awOHtTkcmkilGjRpFpUqVaNWqFdu3u/r9pEmTqF+/PrVq1aJbt26cO3eOVatW8c033/D4448TFRXFrl27PC4nItmfsbp2SPzg9ripB37//ferK8xa2DzTlbDFnYXmI6DpMAjVT90l/W3bto0qVaoA8NycX9n616l0Lb/qDQUY2aFaqsts2LCBvn37smbNGhISEqhTpw4DBw7kvvvuo0iRIgA888wzXHvttTz00EP07duX9u3b0717dwCOHj3qcTkJTu597iJjzAZrbT2HQpIsSiNr4pd0fdyUMVCrBwxZB5Xbw+IXYWIL+Ovnqy9bJAgtX76cLl26kCdPHgoUKEDHjh0B2LJlC9HR0dSoUYPp06fz66+/elzf1+VEJHvRfdbEefmKwR1ToEZ3+PYxmHQLNBkKLZ6C8NxORyfZUFojYBnJ068E+/bty+zZs6lVqxZTp05l6dKlHtf1dTkRyV40sibBo/LtMGQN1O4NK9+Ed5rCnyudjkok3cTExPDVV18RGxvL6dOnmTNnDgCnT5/m+uuvJz4+nunTp19aPn/+/Jw+ffrStLflRCR7U7ImwSV3Iej4FvT5GpISYGo712jb+fS9vkjECXXq1KFHjx5ERUXRrVs3oqOjAXjhhRdo2LAhrVu3pnLlypeW79mzJ6+++iq1a9dm165dXpcTkexNPzCQgNSrV8+uX78+Yz8k7iwsHgWrJ0CBEtD+dbjp1oz9TMm2PF3sLZKR9AMDSS8aWZPgFZEX2rwE/b6HXPngkzvgywFw9qjTkYmIiGQaJWsS/ErVhweXQfMnYMsX8HYD2PKlHlklIiI5gpI1yRrCckHLp11JW6FSMOs++Kw3nDrgdGQiIiIZSsmaZC3XVoN+C+HWF2HXIni7IWz4UKNsIiKSbSlZk6wnNAyaPASDVsF1NWDOw/wx9hZOHbjKJyqIiIgEISVrknUVKQ/3zmFemScpemor4e81YcNnL5KYkOB0ZCIiIulGyZr4xRjTwRgz8eTJk06H4hISQru+T7H/rqVszRVF3d9e5ffRTfjlp9VORyYiIpIulKyJX9L12aDpqHKlytR54js21BvDtYkHqPR1O+a+NYy/jgZJUikiIhIgJWuSbZiQEOq2f4DIR9azq1grbj86hTPjmvLZV7M5H5/odHgiAIwbN44qVarQu3dvmjRpAsCJEyeYMGHCpWVSTvvq2WefZcyYMekWa3o4dOgQderU4amnnqJr164kJSU5HZJIlqNkTbKd3NdcS5WhMznS4UOKhcdyx8a+fPnyfSz4eTd6Yoc4bcKECcybN4/p06ezatUqIP2StWC0bt06evXqxejRoylevDhHj+qm1iL+UrIm2VbRup25ZvhPHL6pJ3clfk3lr9rwwvj32HZAzxkVZwwcOJDdu3fTsWNHXn/9dfLlywfAk08+ya5du4iKiuLxxx+/Yhpg2rRpNGjQgKioKB588EESE12jxaNGjaJSpUq0atWK7du3e/zcTZs2ERMTQ9WqVQkJCcEYw8iRIwPahlmzZtGoUSNq1apFs2bNOHz4MABdunThmWeeITo6muuuu46FCxcCrmStVq1aAJw8eZJixYoF9LkiOZq1Vi+9/H7VrVvXZiUJO5faU/+uZu3IAnb6M53t87NW2aNnLjgdlmSirVu3Oh2Ctdba0qVL28OHD1trrc2bN6+11to//vjDVqtW7dIyKae3bt1q27dvb+Pi4qy11g4aNMh++OGHdv369bZ69er27Nmz9uTJk7Z8+fL21VdfvezzYmNjbaVKleyaNWustdY+88wzdvjw4TYpKemy5Zo1a2Zr1ap1xev777+/bLkjR45c+vvZZ5+148ePt9ZaW6FChUuf/cUXX9i+fftaa63t2bOnHTp0qB04cKBduHBhgLWWNXnqc8B6GwTHcL2y1ivM6WRRJDOElm9O/mFrOf/9i/Rc9w5/b97IyF8GULd1T+5uVJqwUA0yS/BatGgRGzZsoH79+gDExsZSvHhxjh07RpcuXciTJw8AHTt2vGLdhQsXUqdOHRo0aABAzZo1WbBgAcaYy5Zbvny5T7FMnTqVGTNmcOHCBQ4ePMhLL73EuXPnOHnyJI8++igACQkJFCpUCICwsDDeeuutgLZbRFyUrEnOEZGHyNtfgqhuFP5iMG8de5mvFyyj1+pBPNKxCc0qFnU6QhGPrLXce++9jB49+rL333jjjSuSrpS2bNlCjRo1Lk3/9NNP1KlT54rloqOjOX369BXvjxkzhlatWgHw0UcfsXbtWhYvXky+fPmIiYmhWrVq/Prrr9StW5fQ0FAANm/eTPXq1QH4+OOP/dtYEbmChhMk5ylRl1yDl2NbPEWHsHVMOj2EmVNeZ8CH6/jv0XNORyc5UP78+S9LlFJO33LLLcyaNYtDhw4BcOzYMfbs2UNMTAxfffUVsbGxnD59mjlz5lxRdpEiRdi8eTMAO3bs4Msvv6Rnz55XLLd8+XI2btx4xetiogbwyy+/0KRJE/Lly8cXX3zBqlWrqFGjBlu2bCEqKurScps3b6ZmzZpXXS8i4qJkTXKmsAhMiycJGbScAjdUZFzEeHrtHkHv17/k1e9+4+wFPQVBMk+RIkVo2rQp1atX5/HHH79iumrVqrz44ovceuut1KxZk9atW3PgwAHq1KlDjx49iIqKolu3bkRHR19Rdq9evThz5gzVq1dnwIABfPrppxQpUiSgOO+9917GjRtHdHQ0O3bsoFy5cuTNm5dffvnlsmRty5Ytl0bWROTqGWt1KwPxX7169ez69eudDiN9JCXCmnexi17gQpLh+Qu9WJynDU+0q0rnqBJpnmaSrGHbtm1UqVLF6TAkB/HU54wxG6y19RwKSbIojayJhIRC4yGYwT8SWboeL4VPZiLP88bM7+j2zio27zvhdIQiIpKDKVkTuahwWejzDXQYR42QPSzK/RQxRz6j69vLePzzTRw6fd7pCEVEJAfSr0FF3BkDde/FVGxN2Nx/MGz7R9xZeD0PbLyXm7cc5OFbKtC3SVkiwvR/TpY2/0k4+Ev6lnldDWj7cvqWKSKCRtZEPCtwA/T8BLpP4QZ7iG9z/R+jCs1hzLxfuO2NZSz+7W+nIxQRkRxCI2si3hgD1btCuRaYBU/SafPHtCr+IyPiB3D/1LO0qFSMf7avSvli+ZyOVPylETARyUI0siaSljyFoetEuOtz8hLL+NgnmFNxLlv/PMBtry/jvR92oV9Vi4hIRlGyJn4xxnQwxkw8efKk06FkvptuhcGrMfXup8be6fxY6BmGltnP6Pm/MXj6T5zRvdnEB7GxsTRv3vzSg9jTy/3330/x4sVTvb/ZiRMn6N69O5UrV6ZKlSr8+OOPqc47f/48DRo0oFatWlSrVi3gh79nhgULFlCpUiUqVKjAyy97HzktU6YMNWrUICoqinr16vm8fmJiIrVr16Z9+/YAxMXFERMTQ0KC9nvJeErWxC/W2jnW2gEFCxZ0OhRnRBaA9mOh7zxCQ8MZ9tdwFpT9nB+37qbT+BXsPHTl43pE3H3wwQd07dr10qOZ0kvfvn1ZsGBBqss88sgjtGnTht9++41NmzZddg8wT/Ny5crF4sWL2bRpExs3bmTBggWsXr06XeNOD4mJiQwZMoT58+ezdetWPv30U7Zu3ep1+SVLlrBx40Yu3ivSl/XffPPNy+orIiKCW265hRkzZmTMRom4UbImEogyTWHQSmj6CJUPfs3agk8TdXYlncavZN4vB5yOToLY9OnT6dSpEwCzZs2iUaNG1KpVi2bNmnH48OGAy42JiaFw4cJe5586dYply5bRr18/wJVsXHzYurd5xhjy5XNdkxkfH098fHyaN4nu2bMnPXr0oGHDhpQuXZq5c+cGvE2+Wrt2LRUqVKBcuXJERETQs2dPvv7663Rbf9++fcydO5f+/ftftl7nzp2ZPn16um2HiDdK1kQCFZ4bWj8P/RcRUaA4ryW9wnu5x/Ov6UsYPW8bCYlJTkcoQSYuLo7du3dTpkwZAFq2bMnq1avZtGkTrVu3ZubMmZctHx0dTVRU1BWvhQsX+v3Zu3fvplixYtx3333Url2b/v37c/bs2TTnJSYmEhUVRfHixWndujUNGzZM9XM2bdpEuXLlWLNmDdOnT+e5557LsG26aP/+/ZQqVerSdMmSJdm/f7/HZY0x3HrrrdStW5eJEyf6tP6wYcN45ZVXCAm5/CuzevXqrFu3LuC4RXylX4OKXK0SdWDAUljxBk2XvcKyvJt4esXd3LO3I2/1rkPRfLmcjlCCxJEjRy6NZgFMnTqVGTNmcOHCBQ4ePMhLL7102fLLly9Pt89OSEjgp59+4q233qJhw4Y88sgjvPzyy7zwwgupzgsNDWXjxo2cOHGCLl26pPrcz9jYWI4cOXLp2raqVaty/Pjxq9qmVq1acfDgwSveHzVq1KURSk8/8PE2Arhy5UpuuOEGDh06ROvWralcuXKq63/77bcUL16cunXrsnTp0suWCQ0NJSIigtOnT5M/f36/tkvEH0rWRNJDaDg0fxxTpQN5vnmIN/ZN4If9q+g/bhAj776N2jde43SEEgRy587N+fOuJ2F89NFHrF27lsWLF5MvXz5iYmKoVq3aZctHR0dz+vSV10GOGTOGVq1a+fXZJUuWpGTJkpdGxrp3737pQvrU5l1UqFAhWrRowYIFC7wma1u2bKFixYpERkYC8NNPP1GrVq2r2iZfRtxKlizJ3r17L03v27ePG264weOyF98vXrw4Xbp0Ye3atTRt2tTr+itXruSbb75h3rx5nD9/nlOnTnH33Xczbdo0AC5cuHBpe0UyjLVWL738ftWtW9eKF4kJ1v74jk184Vp7dmRx+69nhtmPV+22SUlJTkeWo23dutXpEKy11pYsWdLGxsba4cOH2zfeeMNaa+2sWbNsaGioPXPmzFWV/ccff9hq1ap5nd+sWTP722+/WWutHTlypB0+fHiq8w4dOmSPHz9urbX23LlztlmzZnbOnDnWWmtvvvlmu2/fvsvKnzRpki1RooSNjY21Z86csU2aNLErVqy4qm3yRXx8vC1btqzdvXu3vXDhgq1Zs6bdsmXLFcudOXPGnjp16tLfjRs3tvPnz/d5/SVLltjbb7/90vSRI0ds5cqVvcblqc8B620QHMP1ylovXbMmkt5CQqHRQEKGrCa8dEOeC/2AivN78vK0uZyPT9/bNUjWc+utt7JixQruvfdexo0bR3R0NDt27KBcuXLkzZs34HJ79epF48aN2b59OyVLlmTy5MkAtGvXjr/++guAt956i969e1OzZk02btzI008/fWl9T/MOHDhAy5YtqVmzJvXr16d169a0b9+epKQkdu7cecUPGjZt2kTv3r1p0aIF9evXZ9CgQTRt2jTgbfJVWFgY48eP57bbbqNKlSrceeedl41SXqyDv//+m2bNmlGrVi0aNGjA7bffTps2bdJc35slS5bQrl27jNw0EQCMtbqZp/ivXr169uLP3iUV1pL08zTi5j4FCef5JE9vWvd7gVJFCzgdWY6zbdu2y2694JSff/6ZsWPH8vHHHzsdSsC2bNnCBx98wNixYy97PyYmhkmTJlGpUiWHIstcXbt2ZfTo0V6311OfM8ZssNbW87iCiBcaWRPJSMYQUuceIoet51SpFtwfO5XT42NYt/oHpyMTh9SuXZuWLVum+01xM1P16tWvSNQAdu3aRcWKFR2IKPPFxcXRuXPnHJOYirOUrIlkhvzXUbzf5xxuO5HrzXGi5ndh7eRHSYqLdToyccD999+f7jfFDQb79++/4vYW2VVERAR9+vRxOgzJIXLGXiUSDIyhWMMeRD6yno2FWtFg7wccfLUBZ35f6XRkIiISxJSsiWSy3IWKUW/YDBbWmYCNO0ee6bdzbNajcOGM06GJiEgQUrIm4gBjDK069ubQPUuYFdKGwls+4OybDWDXYqdDy9b0gyrJLOprkp6UrIk4qHaFG2nx6FT+WWQMB88kwcddSPxqMMQeT3tl8UtkZCRHjx7Vl6hkOGstR48e1c1yJd3o1h0SEN26I33FJyYxZu4mCqwdy8CwbyFvUUJvfw2qdnQ6tGwjPj6effv2XXqCgEhGioyMpGTJkoSHh1/2vm7dIYFQsiYBUbKWMeZs+oupX8zmpZCJVOIPqNIR2o2B/Nc6HZqIpAMlaxIInQYVCSIdat3A6CH3MDTvGF5N6EHC9gXYtxvAxk9A/1iJiORIStZEgsxN1+bni4ea8/tNA7gtdhS7KQmzB8G0bnDiv06HJyIimUzJmkgQKhAZznv31KX7bbdw68knGR85gKQ9P8LbjWDNREhKcjpEERHJJErWRIKUMYZBLcrz4f2NmRzXmtviX+VI4dow/3GY0haO/O50iCIikgmUrIkEuWYVi/Ltw9HkLlaGensGMbf8v7CHf4N3msLy1yAx3ukQRUQkAylZE8kCShTKzcwHG9OrwY0M+bUyQ695l7jyt8Ki52FSSziwyekQRUQkgyhZE78YYzoYYyaePHnS6VBynMjwUEZ3rcm/u9Xg+73Q9I++LKs9Fnv6b5jYEhY+C/G6h5iISHaj+6xJQHSfNWdt2X+Skd/8yoY9x6lZxPJ20S8otedLKFIBOo6H0o2dDlFEPNB91iQQGlkTyYKqlyjIrIGNmdSnHrGh+Yne3p1/FXjRdXf+KW1g7nC4cNrpMEVEJB0oWRPJoowxtK56LQuGxfBK95p8f6EqdY4+z3f5OmPXvQ8TGsPOhU6HKSIiV0nJmkgWFxpiuLNeKZYMb8GwdlGMONubbhdGcjDWuG6k+9VAOHfM6TBFRCRAStZEsonI8FAGxJRn2YiWNGrelltjX2RCYmcSN80kaXx9+HW20yGKiEgAlKyJZDMFc4czok1lvn/8NvbWHk6n+BfZdjY/fH4vCZ/2htMHnQ5RRET8oF+DSkD0a9CsY9fhM7z+3VZKbPuAR8O/wITlIqTNS4TXvQeMcTo8kRxFvwaVQGhkTSSbK18sH+PvbkDbgS/zVPF32BhXgvBvH+LQhLYkHf3D6fBERCQNStZEcoioUoUYO6gb53t/w/g8g8lz6Gfi3mrI79+8ik1McDo8ERHxQqdBJSA6DZq1JSVZFq5eT/6FI2ic9BM7wquQ2PEtqtSo73RoItmaToNKIDSyJpIDhYQYbm1Sn7pPLWRZ9VEUj99LuVlt+Hrco+w6eNzp8ERExI2SNZEcLCI8lJjuQwl/eD1/FmtJp2MfEDchhremfc7Bk3rOqIhIMFCyJiLkLXw9lYbO4lTnDymZ6xyDfh/AnNf689q8jZyMjXc6PBGRHE3JmohcUiCqM/kf28D5aj15IGQOXVf3YNi/32bisl2cj090OjwRkRxJyZqIXC53IfLd+Q70+ZoSBcOZwkhy/2cEt786j5nr9pKQmOR0hCIiOYqSNRHxrFwLIh5aA42GcHfYIj5LeJR5X31ImzeXs3LnEaejExHJMZSsiYh3EXmhzUuYft9TtEgRpka8yoizrzH0/e95+qtfOH1e17OJiGQ0JWsikrZS9TEPLoPmT9LarmRFvqc4tX4GbV5fxvLfDzsdnYhItqZkTUR8E5YLWj6FeXAZeYuXYXz4W7yS+DL/mPwdT325WaNsIiIZRMmaiPjn2mrQbyHc+iJN2MyyvE9gN3zEbWN/4IcdGmUTEUlvStZExH+hYdDkIcygVUSWjOLl8EmMT3yOf06ZwxOzNnNKo2wiIulGyZqIBK5Iebh3DrR/g9phf7Ao95MU2PgubccuYcn2Q05HJyKSLShZE5GrExIC9e7DDF5DeIWW/F/YdN5P+D9GT/2Sxz/fpCcgiIhcJSVrIpI+CpaAXp9Bt8lUjjzK/Mj/o8SmN2k3diGLf/vb6ehERLIsJWsikn6MgRrdMUPWElq9C8PCvmBa4hO8+eFn/GPmJk6e0yibiIi/lKyJSPrLWxS6vQ+9PqNM3ji+yvUsVX75N+3H/oeFWzXKJiLiDyVrIpJxKrXFDFlDSN176R86lxlJjzF52oc8NmMjJ87FOR2diEiWoGRNRDJWZEHo8Abc+y3XF8zNpxGjqL/leTqPnc/3GmUTEUmTkjURyRxlozGDVkGTh+gZtoRZicP4bNp7DPvsZ46f1SibiIg3StZEJPNE5IFbX8T0X0jhYtcxOeI1bvn1Ke4YO4fvfj3odHQiIkFJyZqIZL4SdQkZ8AO0eJr24ev5ImkYc6eP4+FPfuKYRtlERC6jZE1EnBEWAS2ewAxcTv4bbmJcxNt0+e0x7h77JfN/OeB0dCIiQUPJmog4q3gVQvr9B24bTfOI35iVNIzln73K0OnrOXrmgtPRiYg4TsmaiDgvJBQaDyZk8I9ElmnAS+GTuWf7UPqNncHczRplE5GcTcmaiASPwmUJ6fM1dHyLepH7mZH0DzbNeI6h09ZyRKNsIpJDKVkTAIwx5Ywxk40xs5yORXI4Y6BOH0KHriX8plY8Hf4pD+54kCFjP+LbzX85HZ2ISKZTspYNGGM+MMYcMsZsSfF+G2PMdmPMTmPMk6mVYa3dba3tl7GRivihwPWE9PoEuk+hSp6TTE96gl0zn+aZWes5H5/odHQiIplGyVr2MBVo4/6GMSYUeBtoC1QFehljqhpjahhjvk3xKp75IYv4wBio3pWwh9cTUqM7j4R9RZ/NfXjyzcnsPnzG6ehERDKFkrVswFq7DDiW4u0GwM7kEbM44DOgk7X2F2tt+xSvQ758jjFmgDFmvTFm/eHDh9N5K0RSkacwId0mQu9Z3JgvibFnRrBi/APM3bDT6chERDKckrXsqwSw1216X/J7Hhljihhj3gVqG2Oe8rSMtXaitbaetbZesWLF0jdaEV9UbE3kw2uJrXUvfcw8anzdhskfTdFpURHJ1pSsZV/Gw3vW28LW2qPW2oHW2vLW2tEZGJfI1YksQN4ub5LQ51vy5Ymk3+5hLH21J3v268cHIpI9KVkLEsaYkelc5D6glNt0SUDfZpJthJWLpvA/1rGnygO0iltI7olNWLfgI6fDEhFJd0rWgsdIY8y/jTGTjDGDjDHXXGV564CKxpiyxpgIoCfwzdWHKRJEwnNTuscYjvWaz9nwa6i/+iG2vNGF88d1I10RyT6UrAUPC5wHvsM1IrbKGFPLlxWNMZ8CPwKVjDH7jDH9rLUJwNDk8rYBM621v2ZM6CLOKl6pESVHrGZpyYFUPL6MuHH1ObLyQ7Bez/yLiGQZxupgFhSMMb9aa6u5Td8EvGutvdnBsK5gjOkAdKhQocIDv//+u9PhiFxh1eqV5FkwjCh2cOjaGIr3mgCFSqW9okgmMMZssNbWczoOyVo0shY8jhhj6l6csNbuAILuJ5fW2jnW2gEFCxZ0OhQRj5o0akrRhxYzMd9A8h5cw4VxDYhfPRGSkpwOTUQkIErWgsfDwDRjzDRjzBPGmOnAH04HJZIVlSySn77DRjO55qesiS9P+ILHOf9+Gzii+7KJSNajZC1IWGs3AVHAp8lvLQF6ORaQSBYXERbCw91uIbbH5zzDYOL+2kLihMaw4nVITHA6PBERn+maNQlIvXr17Pr1650OQ8Qne4+d45lpC+l5eBxtQ9eRdF0tQjqNh+trOh2a5DC6Zk0CoZE1Ecn2ShXOw8TBt7Om/psMjBvGyb//i53YAhY9D/HnnQ5PRCRVStZEJEfIFRbKsx2r0anXQNonvcY3thksfw3ei4b/rnE6PBERr5SsiV+MMR2MMRNPnjzpdCgiAWlb43o+ebgNk4o8Tp+4Jzh56hT2g9tg3gi4cMbp8ERErqBkTfyiW3dIdlC6SF6+GNSEMg070uTUKOblbo9dOxEmNIadi5wOT0TkMkrWRCRHyhUWyvOdqvPvu5rwxLl7uI/nOJsUCtO6wuzBcO6Y0yGKiABK1kQkh2tf8wbmPNSMQ9fUoc7hkay8/l7sps9gQiPYqsfpiojzlKyJSI5XtmhevhzchDsaVaD3H7cx/Jo3ictdHGbeAzPugdN/Ox2iiORgStZERIDI8FBe7FyDcb1qs+BIMZoeeZpdtYbDju/g7Qaw8RM9GF5EHKFkTUTETcdartOiRQrm45Y1dZhU/WNsscowe5Drerbje5wOUURyGCVr4hfdukNygnLF8jF7SFN6NbiRUWsS6BH3T060fAn2rnX9YnTNe3owvIhkGj1uSgKix01JTjH75/08/dUvhBjDCy0K0HnfGMyuhVCqEXR8C4rd5HSIkoXocVMSCI2siYikonPtEix4JIaaJQvy6HfH6HH2Hxxq9SYc2Q7vNoVlYyAx3ukwRSQbU7ImIpKGG4vkYXr/hrzSrSbbDp6m2YJrmVJ7Bkk3tYPFL8CklnBgk9Nhikg2pWRNRMQHxhjurF+KRY815+ZKxXlu8RE6/N2fPa0mwplDMLElLHxWD4YXkXSnZE1ExA/FC0Ty7j11ead3Hf4+dYGb5+XnjcrTSKjZE1a87jo1uudHp8MUkWxEyZqISADa1rieRY81p1udEryx4jCtd93JtlYfQWIcTGkDc4fDhdNOhyki2YCSNRGRABXME84r3WsxrV9DEpKSaPttGM+VfJ+4eg/Cuvddt/n4faHTYYpIFqdkTUTkKjWrWJTvhsXQv1lZPtxwhJjNt7Hu5s8gPA9M7wZfDdSD4UUkYErWxC+6Ka6IZ3kiwnimfVW+HNyUgrnDuWNeIo9e8xbnGj0Gv3zuemTVr7P1yCoR8ZtuiisB0U1xRbyLS0jinaW7GL/kd/LlCmNMTAg3b38Bc2AjVG4Pt78G+a9zOkxxgG6KK4HQyJqISDqLCAvhkVYVmftwNGWK5qXfggv0C3+ZE83+CTsXukbZfp6mUTYR8YmSNRGRDHLTtfmZNbAJIztUZfWfJ2myrDpfNpyJLV4Nvh4CH3eG4386HaaIBDklayIiGSg0xHBf07J8NyyGuqWv4bFFZ7jj/P9xKGY07Nvg+sXo6ncgKdHpUEUkSClZExHJBKUK5+Gj+xvw2h21+P3wOZotKsuUqE9JKt0UFjwJH7SBQ785HaaIBCElayIimcQYQ7e6JVn4WHNaV7uW55adot3hh/hvizfg6E54Lxp+eFUPhheRyyhZExHJZMXy5+Ltu+owqU89jsfG0+K74rxe6WMSKrWHJS/CxBbw189OhykiQULJmoiIQ1pXvZbvH2tOj/o38ubqk9z8Zx+2NX8Pzh2FSTfD9/+C+FinwxQRhylZExFxUIHIcEZ3rcGnDzQixEDb7/IzstQHxNXsDSvfhHeawp8rnQ5TRBykZE1EJAg0Ll+EBcNieLB5OaZtPEGzrZ1ZGzMVbCJMbQffPgbnTzkdpog4QMma+EWPmxLJOJHhoTzVtgqzBzelSL5c3PmfCB4qNIHjtQbAhikwoRHs+I/TYYpIJtPjpiQgetyUSMaKT0xi4rLdjF+8k/MJiQwsf5xHzr5J5PEdUONOaPMy5C3idJjiJz1uSgKhZE0ComRNJHMcPXOBKSv/5MNVf3L+wnlevfZ7Op7+jJDIgtDuFajWFYxxOkzxkZI1CYSSNQmIkjWRzHUyNp6Pf/yTySv+oHjsLt7O9wEV4ndgK7XF3P46FLje6RDFB0rWJBBK1iQgStZEnHH2QgKfrv0v7//wO7fHfs2I8FmEhEcQ1mYUpk4fjbIFOSVrEgj9wEBEJAvJmyuM/tHlWPpEK8p2eIJ7Il5n3flSmDkPc2RCGxKP7HY6RBFJZ0rWRESyoMjwUO5uVJrpT/Tir04zeS3XECIObSJ+fEM2z3yR+Hg9skoku9BpUAmIToOKBJfEJMuSdRvJ9/0IGiWsZaupyO6mL9O6RUtyhYU6HZ4k02lQCYSSNQmIkjWR4GSTktjynyncuOZZciedZWpYd8Ji/kHPxuXJExHmdHg5npI1CYROg4qIZCMmJIQabfpRYPjPnCjbngGJM2iyqBsDXp7E20t2cuq8To+KZDVK1kREsiGTtyjF+34Ed82kXL4EPkp6mvBF/+KWl+fz2n+2c+xsnNMhioiPlKyJiGRnN91GxMNrCanXlwFhc5kXNoL1S7+m2b8XM2ruVg6dOu90hCKSBl2zJn4xxnQAOlSoUOGB33//3elwRMQff66Abx6CY7tZVagDg/7uRGxoPnrUK8WDzctR8po8TkeY7emaNQmEkjUJiH5gIJJFxZ2DpaPhx/Ek5CnO9KLDeHFnaayFLrVLMKhFecoVy+d0lNmWkjUJhJI1CYiSNZEsbv9P8PVQOPQr5yp1YXxEfyb/fJr4xCTa1bieIS0rUOX6Ak5Hme0oWZNA6Jo1EZGcqEQdGLAUWj5Dnp1zGbHzHtZ1PM6A6HIs+e0Qbd9czsivtzgdpYigZE1EJOcKi4Dmj8ODy6FIeQrMG8STJ55l1ZDKDGtVkdo3XuN0hCKCkjURESleGe7/Dtq8DH8so+DkZgwrtILOta53OjIRQcmaiIgAhIRCo0EwaJXrFOm3j8KHHeDoLqcjE8nxlKyJiMj/FC4Lfb6GjuPh4C/wThNYOQ4SE5yOTCTHUrImIiKXMwbq3AND1kCFVvD9P2FyKzioHxyIOEHJmoiIeFbgeugxDe6YCif3wcTmsHgUJFxwOjKRHEXJmoiIeGcMVOsCQ9ZCjTtg2SvwXgzsXed0ZCI5hpI1ERFJW57C0OVd6D0LLpyBya1hwVMQd9bpyESyPSVrIiLiu4qtYchqqN8fVk+ACY1h91KnoxLJ1pSsiYiIf3Llh9vHwH3zITQcPurkenRV7AmnIxPJlpSsiYhIYEo3gYErodmjsPETeLshbPvW6ahEsh0layIiErjwSGj1LDywGPIVgxm9Yea9cOaQ05GJZBtK1kRE5OrdEAUPLIGb/wnb58HbDWDTZ2Ct05GJZHlK1sQvxpgOxpiJJ0+edDoUEQk2oeEQM9x1arToTfDVgzD9Djix1+nIRLI0JWviF2vtHGvtgIIFCzodiogEq2I3wX0LoO0rsGcVTGgEaydBUpLTkYlkSUrWREQk/YWEQMMHYfCPULI+zBsOU2+HIzudjkwky1GyJiIiGeea0nDPV9BpAhz61fVg+BWv68HwIn5QsiYiIhnLGKjdG4asg5tug4XPwvs3w4HNTkcmkiUoWRMRkcyR/1ro8THc+RGcOgATW8Ci5yH+vNORiQQ1JWsiIpK5qnaCIWugVk9Y/hq8Fw3/XeN0VCJBS8maiIhkvjyFofMEuPtL18jaB7fBvBGuh8SLyGWUrImIiHMq3OL6xWiDAbB2ouvB8DsXOR2VSFBRsiYiIs7KlQ/avQL3L4CwXDCtK8weDOeOOR2ZSFBQsiYiIsHhxkYwcAVE/8P1qKq3G8LWr52OSsRxStZERCR4hEfCLf+CAUsh/3Uwsw/MuAdO/+10ZCKOUbImIiLB5/qa8MBiaPUs7PjO9WD4n6frwfCSIylZExGR4BQaDs0ehUEroXgV+Hqw63q243ucjkwkUylZExGR4Fa0IvSdB+3GwN61rl+MrnlPD4aXHEPJmoiIBL+QEGjwgOs2H6Ubw/wRMKUNHN7udGQiGU7JmoiIZB2FboTes6DLe3BkB7zbDJaNgcR4pyMTyTBK1kREJGsxxvWoqiFrofLtsPgFmNQS/trodGQiGULJmoiIZE35isMdU6HHdDhzCCbdDAufhfhYpyMTSVdK1kREJGur0t71YPiou2DF665To3tWOR2VSLpRsiYiIllf7mug03i4ZzYkxsGUtjD3H3DhtNORiVw1JWsiIpJ9lG8Jg1dDo8GwbjK83Qh+/97pqESuipI1ERHJXiLyQpvR0O8/rofET+8OXz6oB8NLlqVkTfxijOlgjJl48uRJp0MREUldqQbw4DKIGQFbZsH4+rDlSz2ySrIcJWviF2vtHGvtgIIFCzodiohI2sJywc3/BwN+gIIlYdZ9MONuOHXA6chEfKZkTUREsr/rqkP/RdD6edi5EN5uCD99pFE2yRKUrImISM4QGgZNH4FBq+C6GvDNQ/BRJzj2h9ORiaRKyZqIiOQsRcrDvXOg/euw/yd4pwn8OAGSEp2OTMQjJWsiIpLzhIRAvftdN9MtEw3fPQWTb4VD25yOTOQKStZERCTnKlgC7poBXd+HY7vh3Wj44RVIiHM6MpFLlKyJiEjOZgzUvAOGroOqHWHJKJjYAvZvcDoyEUDJmoiIiEveotD9A+j5KcQeg/dbwX/+CXHnnI5McjglayIiIu4qt3Ndy1anD6waB+82hT+WOx2V5GBK1kRERFKKLAgd3nT9atRa+LA9zBkG5/X0Fsl8StZERES8KRvjui9b46Hw04d6MLw4QsmaiIhIaiLywG2joN9CyF0ITu13OiLJYcKcDkBERCRLKFnX9YzR0HCnI5EcRsmaiIiIr8IinI5AciCdBhUREREJYkrWRERERIKYkjURERGRIKZkTURERCSIKVkTERERCWJK1kRERESCmJI1ERERkSBmrLVOxyBZkDHmMLDHw6yCQMqH56V8ryhwJINCS4un+DKjHF+XT2u51OZ7m+dLm4Bz7eJUm/izTnq3i69tpX0l8OWCdV8pba0tFuC6klNZa/XSK91ewMS03gPWB1N8mVGOr8untVxq873N86VNnGwXp9rEyXbxta20r2Rem/jTVk62i14586XToJLe5vj4nlPSKxZ/y/F1+bSWS22+t3lqk6tfJ73bxZ+2cor2Fd8+RyTD6TSoZDpjzHprbT2n45DLqV2Cj9okOKldJLNpZE2cMNHpAMQjtUvwUZsEJ7WLZCqNrImIiIgEMY2siYiIiAQxJWsiIiIiQUzJmoiIiEgQU7ImjjPGlDPGTDbGzHI6FvkfY0xnY8wkY8zXxphbnY5HwBhTxRjzrjFmljFmkNPxiIsxJq8xZoMxpr3TsUj2pGRNMoQx5gNjzCFjzJYU77cxxmw3xuw0xjwJYK3dba3t50ykOYuf7TLbWvsA0Bfo4UC4OYKfbbLNWjsQuBPQrSMyiD9tkuwJYGbmRik5iZI1yShTgTbubxhjQoG3gbZAVaCXMaZq5oeWo03F/3Z5Jnm+ZIyp+NEmxpiOwApgUeaGmaNMxcc2Mca0ArYCf2d2kJJzKFmTDGGtXQYcS/F2A2Bn8khaHPAZ0CnTg8vB/GkX4/JvYL619qfMjjWn8HdfsdZ+Y61tAvTO3EhzDj/bpCXQCLgLeMAYo+9VSXdhTgcgOUoJYK/b9D6goTGmCDAKqG2MecpaO9qR6HIuj+0CPAS0AgoaYypYa991Irgcytu+0gLoCuQC5mV+WDmaxzax1g4FMMb0BY5Ya5MciE2yOSVrkpmMh/estfYoMDCzg5FLvLXLOGBcZgcjgPc2WQoszdxQJJnHNrn0h7VTMy8UyWk0XCuZaR9Qym26JPCXQ7HI/6hdgo/aJPioTcQxStYkM60DKhpjyhpjIoCewDcOxyRql2CkNgk+ahNxjJI1yRDGmE+BH4FKxph9xph+1toEYCjwHbANmGmt/dXJOHMatUvwUZsEH7WJBBs9yF1EREQkiGlkTURERCSIKVkTERERCWJK1kRERESCmJI1ERERkSCmZE1EREQkiClZExEREQliStZEREREgpiSNREREZEgpmRNRLINY8xXxpgXjTHLjTEHjTGtnI5JRORqKVkTkeykOnDCWhsNDAZ6OxyPiMhVU7ImItmCMSYPUBB4PfmtMOCEYwGJiKQTJWsikl1UAzZYaxOTp2sCWxyMR0QkXShZE5Hsojqw0W26JrDZmVBERNKPkjURyS5qcHmyVh2NrIlINmCstU7HICIiIiJeaGRNREREJIgpWRMREREJYkrWRERERIKYkjURERGRIKZkTURERCSIKVkTERERCWJK1kRERESCmJI1ERERkSD2/5/Iks43N8WrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "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": null,
   "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",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    return number_hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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
}