cds-monte-carlo-methods/Exercise sheet 1/exercise_sheet_01.ipynb

{
 "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, 3.96949304e-01],\n",
       "       [3.20000000e+01, 2.84220425e-01],\n",
       "       [6.40000000e+01, 2.15515371e-01],\n",
       "       [1.28000000e+02, 1.59275704e-01],\n",
       "       [2.56000000e+02, 1.00798933e-01],\n",
       "       [5.12000000e+02, 6.87799562e-02],\n",
       "       [1.02400000e+03, 4.87684636e-02],\n",
       "       [2.04800000e+03, 3.85210905e-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/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+9ElEQVR4nO3dd3gU1f7H8fdJAqEHqQoBAcFQAqH3BJAiIB0UFBUQUZqiXqzXe/VexYoNBRUEsSDiD0U6KgJK70WKdBAQpElogYTk/P7YwA1hN9ksSWaTfF7Ps8+T2Zk5+51zzsx+c3aKsdYiIiIiIv4pwOkARERERMQzJWsiIiIifkzJmoiIiIgfU7ImIiIi4seUrImIiIj4MSVrIiIiIn5MyZqIiIiIH1OyJiIiIuLHlKz5EWNMmDFmvTHmjDHmUafjucwYM9EY87IT5aXXZxtjthhjmju1vo+fmab+kLyunIjZnxhj9hljWvm4brr2efEP19MnRJJKy/E1Pfqd18maMeYeY8waY8xZY8xhY8xcY0zT6/lwucZTwCJrbUFr7ShvV9IBKHXW2mrW2kXeLOuuPtOyfjryqT9cllExq79JVqB+mn0YY4KNMeONMfsT/3ldb4xpl2yZIsaYacaYc4nL3ePNPA+fl2rfyezvBK+SNWPME8C7wCtASaAsMAbonGGRpYExJsjpGNLJzcAWp4NIL9moXZySIf1B7SL+zuk+6vTnyzWCgANAMyAE+BfwjTGmXJJlRgOxuHKU3sCHxphqXsxLM0f6h7U2xReuijkL3JnCMlWARcApXF8unZLM2wcMBzYB0cAUIE/ivGeAqcnKeg8Ylfh3KeBb4BiwF3g0WblPJ5Z7EVdj1gbWA2eA/0v8rJeTrJNaeW7jTJxfBvgucd0TwAeplZmWugIWAPHAhcT6vtXNuk8DhxK3bzvQEvgCSABiEtd7Kknd7k5cdivQNQ3bWgtYl7juFODrZPWYWtnJ2yXF8pJtY2qf7ba+Sb0v7QNapbYNKdTnlfU9taE3dZvO/SG1ukoas7t2SbHv4qbPe6ofb/pqan0nSZxPJsZ5DhiP6wA7N3Gd+cANSZZ9NrGcv4FPubofJ93+1LbVqz4KBAOnE9vmbOLrcju1crN8asckC1RMMj3x8nwvYvZUx57e96m8NPRZn4/nnvposvkDgdm4vniPA38CrT3EmNJ+7On76JrPT2Fb+wEzk3zeLuCbJNMHgJpprNf03Dc8HqO8qIcU+6ybuB8EfgA+xLUf7gCqAsOAPxLbqpun9X19JcbePfHv/LiSsVuTzP8CeC2leT70naT94yBefKe4ORZ51R+uicuLCmkLXCLZjpNkfq7EjvockBu4LTGIsCRBrsK1oxYBtgEDE+fdDJwHCiVOBwKHgYa4Rv3WAv9OLLcCsAe4PUm5G3B9oeRNXGZ/YgfJBXRLbKDLBz5vyvMUZyCwEXgnseHzAE1TK9OHuloEPOhh3TBcB4BSidPlgFuSd4Qky9+ZuC0BQE9cO/dNXmzr5Xp8PDHeHkAcV3/BpFa2u3bxWF6SclNcNqX6JoW+5K6OvNiG5PW5D2jlRRt6rNt07g/etNOV7XDTLqntD277vKf6SUNf9VjvScpegetLqDRwFFcSVQtXorQAeCHJspsTt6kIsNTd9nuxrV730cTl+wE/JpneAUSl0EZuj0mJy7hN1ryI2W0dp/C+T+Wlpc9yHcdzd33UzWePAU7i2t8DgBeA+Sn0w3243489Hfuu+vxUtrUCrkQoALgpsZ0PJZZTAVfSEuBtvabnvpFS3KnVA170WTdxj8L1z1zLxDb/NlkZw4B1btablViH7l6zPH1e4rolcf2DVDlxuhYQk2yZ4cDMlOb50HeS9o+rlkmp/fjfscjr/nBNTKku4BoyPJLC/EjgCBCQ5L3JwItJgrw3ybw3gI+STC8B7k/8uzWwO/HvBsAfyT7rWeDTJOU+kGReFK5s1SQr++U0lOc2TqARrv8Gk/+nl2KZPtTVIjx/OVfEtXO2AnKl1rHcrL8B6OzFtkbh+o81aT0uI+WdNXnZydvFq/JSW9aLNnTbl7ypIzfb4ClZu67+no79IdV65dpkLWm7pFaXbvu8l3Xpsa+mVO9Jyu6dZPpb4MMk048A3ydZdmCSee3dtbkX25qmPg+8DbyV+Hd+XCNrRTy0kcdjUuK0p2QttZjd1nEK7/tUng991qfjubs+6uazl5BkJBfXl6MvyZqnY99Vn+/Fth7ANQrVCxiLK/mpjCuZn5HWfSG99o3U4k6pHvCiz7qJ8xdgeJLpl0iSbOFK4jalZdtTqZdcuEYRP07eVsmWG4DrGOpxng9954GUlvHUfvzvWORzf/DmnLUTQLEUfqMtBRyw1iYkeW8/rsz/siNJ/j4PFEgy/RVwd+Lf9yROg+u/tFLGmFOXX7j+UyiZZN0DyeI4ZBNrxs18b8rzFGcZYL+19hJX86bMpLypK7estbuAx4AXgaPGmK+NMaU8LW+Mud8YsyFJXOFAsSSLeNpWd/W4P41lp9YuV5WXhmVTq29PfekaXmyDJ9fb39NaVkpxeFuvl6Vlf/DU569hjOmdeOHRWWPM3JT6qpf1/leSv2PcTCetz6TbtB9XvSSX2ramtS6rA78l+fuwtfakm+VSOyalJMWYPdVxCnXvU3ketimlPns9x/PU6qc6rpGSy8Jx/dyUVintn8mPXSlt6y9Ac1wJzi+4EoNmia9fIG3H7XTcN7w9rrirB1/6bA1co2SXVXUz/XsqZXjFGBOA62fKWGBokllngULJFi+Ea0QxpXlp5bEuvGm/tH6PJ+VNsrYc13BjFw/z/wTKJFbiZWVxZefe+D+guTEmFOjK/3buA8Bea23hJK+C1tr2SdZN2qEOA6WNMSbJe2WS/O1NeZ4cAMq6SVjTWuZ11ZW19itrbVNcBz4LvH55VtLljDE3A+Nwdeai1trCuH4uSlo3nrirx7JpLDu1dimLe6ktm1p9e+pLV/FiG6y79RJdb39Pr7LSUq+XJT8Ap1SXnvp88nKw1k6y1hZIfLVLfO+avnqd/dKTpPt4WVx1mlxq25rWukyarEUk+Tu51I5J4PqSzJdk+kYvY/Z4PPDwvs/lJZNan72e4zl42PeMMeVxnUO2PcnbtXCNXniS0n7szTqpbevlZC0y8e9fSJasgXf1ms77RnofV5L32SsS486N61SAy2pydbvUwE07GdcdJc56eM11s7zhf+fpdbfWxiWZvQMIMsZUSvJeBK7z9VKa54mnvuOpf3rdfl7uZ9dINVmz1kbjOs9gtDGmizEmnzEmlzGmnTHmDWAlrt9mn0p8vznQEdcJuqmy1h7D9R/Jp7h25m2Js1YBp40xTxtj8hpjAo0x4caYeh6KWo7r54ihxpggY0xnoH6S+WktL6lVuDrxa8aY/MaYPMaYJj6U6XNdGdc9t24zxgTjSp5jErcXXP9dVUiyeH5cneBY4rr9cGX53liO6xzFRxPrsRtX12Nay06tvLQsm2J9p9CXkkttG5LXZ1LX1d/Tsay01Ks7qfVdT30eUq6flPrq9fRLT4YYY0KNMUVwjdRM8WFbva5LY0wxoDj/G9GpzNVfVEmldkwC15fYPYkxtcX1RZ9qzJ7qOIW696k8N9uUYp9Nx+N5cjWA35KNFtXCdV6lJyn2Uy+ktn/+ArTAdX7dQWAxrnO8i+I6QT8t9Zqe+8b1HldS67NJRZCkXYwxhXAlIZuSLXNNO1lr2yX5Jy/5q13y5XFdwFAF6GitjUlW1jlcF0P9N/F41QTX3Sq+SGleCtuV1r7jVfuloT9cw6tbd1hr3waeAJ5PDOYArgzye2ttLNAJaIfrqo8xuM5ZSMuw51e4fsO9MhJirY3H1cFq4rpy6DjwCa6rU93FGIvrZMj+uE5QvBfXUOxFX8pLVvbldSviurrlINDTxxh9ratgXFe2HMc1fF0C15cTwKvA88Y1/DrcWrsVeAvXjvcXrpGApV58RtJ67IvrJNmeuDr65flpKju18tL42d7U9zV9yc3npLYNV9Wnmxivt79fd1lpqVcP66dYl576fOLqHusnkdu+ej39MgVfAT/iOll9D67zvdK6rWmpyxq4zsO6/GVxCLjLGNPAzeemeExKNCwxtlO4zg/+3puY8Xw88FT3vpbnbptS67PXfTx346rRGWNMUVyjkJtTWCe1fpqi1LbVWrsD109sixOnT+Pqg0sTtxe8r9d02zfS6biSUp9NKoJrR9F2WWvPw5WfLcNJeQQ0VcY1cvUwrv5zxPxvBK53ksUG4zrx/yiuc/QGWWu3eDHPnTT1nTS0n1f9wR1jr/ppOnsxxqzEdfLop07HIiLpyxizD9cFGPOdjsVbOiZJVqM+6x+y1eOmjDHNjDE3Jg7f9sGV5c9zOi4RyZl0TJKsRn3WP2W3uzSHAd/guqplN9DDWnvY2ZBEJAfTMUmyGvVZP5StfwYVERERyeqy1c+gIiIiItmNkjURERERP5bdzlmTTFKsWDFbrlw5p8MQEclS1q5de9xaW9zpOCRrUbImaWKM6Qh0rFixImvWrHE6HBGRLMUYk9pj4USuoZ9BJU2stTOttQ+FhHh7L0sRERG5HkrWRERERPyYkjURERERP6Zz1kQkR4iLi+PgwYNcuHDB6VAkB8iTJw+hoaHkypXL6VAkG1CyJiI5wsGDBylYsCDlypXDGON0OJKNWWs5ceIEBw8epHz58k6HI9mAfgYVkRzhwoULFC1aVImaZDhjDEWLFtUorqQbJWuSuWLPw6b/Az3mTBygRE0yi/qapCcla5K51n8J3z0IX3SBv/c5HY2Io1588UVGjhzpcf7333/P1q1bMzEiEfFHStYkc9V7EO54Cw6ugTGNYMVHkBDvdFQifknJmoiAkjXJbAEBroRt8Aq4uQnMexomtIVj252OTCRTjBgxgrCwMFq1asX27a5+P27cOOrVq0dERATdu3fn/PnzLFu2jBkzZvDkk09Ss2ZNdu/e7XY5Ecn+jNW5Q5IGSR43NWDnzp3XV5i1sOkbV8IWew6aPQVNHoNAXeou6W/btm1UqVIFgP/M3MLWP0+na/lVSxXihY7VUlxm7dq19O3bl5UrV3Lp0iVq167NwIED6devH0WLFgXg+eefp2TJkjzyyCP07duXDh060KNHDwBOnDjhdjnxT0n73GXGmLXW2roOhSRZlEbWJE3S9XFTxkBETxiyGip3gAUvw9jm8Of66y9bxA8tXryYrl27ki9fPgoVKkSnTp0A2Lx5M5GRkVSvXp1JkyaxZcsWt+t7u5yIZC+6z5o4r0BxuPNTqN4DZj0B41pC46HQ/FnIldfp6CQbSm0ELCO5u0qwb9++fP/990RERDBx4kQWLVrkdl1vlxOR7EUja+I/Kt8BQ1ZCrd6w9D34sAnsW+p0VCLpJioqimnTphETE8OZM2eYOXMmAGfOnOGmm24iLi6OSZMmXVm+YMGCnDlz5sq0p+VEJHtTsib+JW9h6PQ+3D8dEi7BxPau0bYL6Xt+kYgTateuTc+ePalZsybdu3cnMjISgJdeeokGDRrQunVrKleufGX5Xr168eabb1KrVi12797tcTkRyd50gYH4pG7dunbNmjUZ+yGx52DBCFgxBgqVhg7vwK1tMvYzJdtyd7K3SEbSBQaSXjSyJv4rd35o+wr0/wmCC8BXd8J3D8G5E05HJiIikmmUrIn/K1MPHv4Vmj0Nm7+F0fVh83d6ZJWIiOQIStYkawgKhhbPuZK2wmVgaj/4ujecPux0ZCIiIhlKyZpkLSWrQf/50OZl2P0zjG4Aaz/TKJuIiGRbStYk6wkMgsaPwKBlcGN1mPkofN4JTu51OjIREZF0p2RNsq6it0CfmVxo+zb20HrXg+GXj9aD4UVEJFtRsiZpYozpaIwZGx0d7XQoLgEBvHCwHvcGj+LUjY3gh+dgfBs4us3pyERERNKFkjVJk3R9Nmg6aV21JHtjQ6i5qz9Tyr5Awsm98FEkLHodLsU6HZ6IiMh1UbImWV6rqiX58YlmPNCkAs/uDOP2uDc5VPp2WPQKjG0Gh9Y6HaLIFaNGjaJKlSr07t2bxo0bA3Dq1CnGjBlzZZnk09568cUXGTlyZLrFmh6OHj1K7dq1efbZZ+nWrRsJCQlOhySS5ShZk2yhQHAQ/+5YlelDmhIcUoImO+/h3eL/5dK5v+GTVvDDPyH2vNNhijBmzBjmzJnDpEmTWLZsGZB+yZo/Wr16NXfffTevvvoqJUqU4MQJ3dRaJK2UrEm2Uj00hO8HN+H5O6ow9q8wGp1+ha03dYXlH8CHjWHvYqdDlBxs4MCB7Nmzh06dOvHOO+9QoEABAJ555hl2795NzZo1efLJJ6+ZBvjyyy+pX78+NWvW5OGHHyY+3nUhzYgRIwgLC6NVq1Zs377d7edu3LiRqKgoqlatSkBAAMYYXnjhBZ+2YerUqTRs2JCIiAiaNm3KsWPHAOjatSvPP/88kZGR3HjjjcyfPx9wJWsREREAREdHU7x4cZ8+VyRHs9bqpVeaX3Xq1LH+7sDJc/aBT1fZm5+eZZ9+8wN7YWR1a18oZO2MR62NOeV0eJLJtm7d6nQI1lprb775Znvs2DFrrbX58+e31lq7d+9eW61atSvLJJ/eunWr7dChg42NjbXWWjto0CD72Wef2TVr1tjw8HB77tw5Gx0dbW+55Rb75ptvXvV5MTExNiwszK5cudJaa+3zzz9vhw8fbhMSEq5armnTpjYiIuKa108//XTVcsePH7/y94svvmg/+OADa621FStWvPLZ3377re3bt6+11tpevXrZoUOH2oEDB9r58+f7WGtZk7s+B6yxfnAM1ytrvYKcThZFMkroDfn4pE9d5m0+wgszgplx9kXGlfmJxus+x+z40fVg+LC2Tocpkqqff/6ZtWvXUq9ePQBiYmIoUaIEJ0+epGvXruTLlw+ATp06XbPu/PnzqV27NvXr1wegRo0azJs3D2PMVcstXuzdqPPEiROZMmUKFy9e5MiRI7zyyiucP3+e6OhoHn/8cQAuXbpE4cKFAQgKCuL999/3abtFxEXJmmRrxhjaVb+JJpWK8ea87dy7Mpjm+SMYFTCegpN7QngPaPc65C/mdKgiHllr6dOnD6+++upV77/77rvXJF3Jbd68merVq1+ZXrduHbVr175mucjISM6cOXPN+yNHjqRVq1YAfP7556xatYoFCxZQoEABoqKiqFatGlu2bKFOnToEBgYCsGnTJsLDwwH44osv0raxInINnbMmOUKhPLl4qUs43w5qzOECVan91/NMv6EPdut014Phf5uqR1aJYwoWLHhVopR8umXLlkydOpWjR48CcPLkSfbv309UVBTTpk0jJiaGM2fOMHPmzGvKLlq0KJs2bQJgx44dfPfdd/Tq1eua5RYvXsyGDRuueV1O1AB+++03GjduTIECBfj2229ZtmwZ1atXZ/PmzdSsWfPKcps2baJGjRrXXS8i4qJkTXKU2mVvYOYjTXmibThPH29H50uvcjRXKfi2P0zuBdGHnA5RcqCiRYvSpEkTwsPDefLJJ6+Zrlq1Ki+//DJt2rShRo0atG7dmsOHD1O7dm169uxJzZo16d69O5GRkdeUfffdd3P27FnCw8N56KGHmDx5MkWLFvUpzj59+jBq1CgiIyPZsWMHFSpUIH/+/Pz2229XJWubN2++MrImItfPWI0miA/q1q1r16xZ43QY1+WPE+f55/e/sXTnUf5ZdBH9LnxJQGAuaPNfqN0XAvS/THaybds2qlSp4nQYkoO463PGmLXW2roOhSRZlL6NJMcqWzQfnz9Qn3d61ebDi225LeY19uUJg1mPux4Mf2K30yGKiIgoWZOczRhD55qlmf9EMxrWqUPzo4/zWtBgLh1a77ov29JREH/J6TBFRCQH09WgIkDhfLl5rXsNutUO5blpBZl2tCrjin5FjZ/+BVumQecPoGQ1p8OU9DL3GTjyW/qWeWN1aPda+pYpIoJG1kSuUr98EWY/2pTerRvSI/oRhvMYF47vxX4cBQtfgUsXnQ5RRERyGI2siSQTHBTIoy0r0aHGTfxzWhEa7qnCe4Wn0OyX12HrdOj0AZSp53SYcj00AiYiWYhG1kQ8qFC8AF8NaMDzdzblsdhB9I97itPRJ7HjW8O85yD2nNMhiohIDqBkTdLEGNPRGDM2Ojra6VAyhTGGHnVC+fkfzSkc0YHGp19hetDtsGI0jGkEexY5HaJkMTExMTRr1uzKg9jTywMPPECJEiVSvL/ZqVOn6NGjB5UrV6ZKlSosX74cgAsXLlC/fn0iIiKoVq3aNQ95j4+Pp1atWnTo0CFdY05P8+bNIywsjIoVK/LaaymPnLrbnvfee4/w8HCqVavGu+++e9Xy7ubFxsYSFRXFpUu6AEkynpI1SRNr7Uxr7UMhISFOh5KpiuTPzVt3RTD2wRa8l2cQd138F0fPx8PnnWHGIxBzyukQJYuYMGEC3bp1u/JopvTSt29f5s2bl+Iyw4YNo23btvz+++9s3Ljxyj3AgoODWbBgARs3bmTDhg3MmzePFStWXFnvvffe8+t71MXHxzNkyBDmzp3L1q1bmTx5Mlu3bvW4fPLt2bx5M+PGjWPVqlVs3LiRWbNmsXPnzhTn5c6dm5YtWzJlypQM3z4RJWsiadC4YjHmDoukYYuO3HZuBJ/SmYR1X2JHN4DfZzsdnmQBkyZNonPnzgBMnTqVhg0bEhERQdOmTTl27JjP5UZFRVGkSBGP80+fPs2vv/5K//79AcidO/eVh60bYyhQoAAAcXFxxMXFXXnm6MGDB5k9ezYPPvigV3H06tWLnj170qBBA26++WZmz874/WLVqlVUrFiRChUqkDt3bnr16sX06dPdLutue7Zt20bDhg3Jly8fQUFBNGvWjGnTpqU6r0uXLkyaNCnDt09EyZpIGuXJFcgTbcKY9mhLfio9mM4X/8uuc8Hw9T0kfNMXzh51OkTxU7GxsezZs4dy5coB0KJFC1asWMHGjRtp3bo133zzzVXLR0ZGUrNmzWte8+fPT/Nn79mzh+LFi9OvXz9q1arFgw8+yLlz/zvvMj4+npo1a1KiRAlat25NgwYNAHjsscd44403CPDyiR4bN26kQoUKrFy5kkmTJvGf//wnw7bpskOHDlGmTJkr06GhoRw65P7Rce62Jzw8nF9//ZUTJ05w/vx55syZw4EDB7yat3r1ap/jFvGWrgYV8VGlkgX5akBDlu6qyDPzatDo8Bc8unUa7FxA0B1vEBDRExJHJ0QAjh8/fmU0C2DixIlMmTKFixcvcuTIEV555ZWrll+8eHG6ffalS5dYt24d77//Pg0aNGDYsGG89tprvPTSSwAEBgayYcMGTp06RdeuXdm8eTP79u2jRIkS1KlTh0WLFqX6GTExMRw/fvzKOW9Vq1bl77//vq5tatWqFUeOHLnm/REjRlwZoXT32ETjZt+bNWuW2+2pUqUKTz/9NK1bt6ZAgQJEREQQFBSU6rzAwEBy587NmTNnKFiwYJq2SyQtlKyJXKcmFYvReEgUP2+rzNC5zXk4+l3qfP8wx1ZMoliv0ZjCZZ0OUfxE3rx5uXDhAgCff/45q1atYsGCBRQoUICoqCiqVbv6xsuRkZGcOXPmmnJGjhxJq1at0vTZoaGhhIaGXhkx69Gjh9sT8QsXLkzz5s2ZN28eJ06cYMaMGcyZM4cLFy5w+vRp7r33Xr788ku3n7F582YqVapEnjx5AFi3bh0RERHXtU3ejLiFhoZeGe0C10+dpUqVuma5pUuXetye/v37X/mJ+LnnniM0NPTKeinNu3jx4pXtFckw1lq99Erzq06dOlauFR+fYGes/8OOGvGEPffv4vb8iyXtjlnv2IT4S06HluNt3brV6RCstdaGhobamJgYO3z4cPvuu+9aa62dOnWqDQwMtGfPnr2usvfu3WurVavmcX7Tpk3t77//bq219oUXXrDDhw+31lp79OhR+/fff1trrT1//rxt2rSpnTlz5lXrLly40N5xxx1Xpm+77TZ78ODBq5YZN26cLV26tI2JibFnz561jRs3tkuWLLmubfJGXFycLV++vN2zZ4+9ePGirVGjht28eXOK6yTfnr/++staa+3+/fttWFiYPXnyZKrzjh8/bitXruzxM9z1OWCN9YNjuF5Z66Vz1kTSUUCAoWPNMgx6+k0WtpzOZnMrlVa/wLbXoti0cY3T4YkfaNOmDUuWLKFPnz6MGjWKyMhIduzYQYUKFcifP7/P5d599900atSI7du3Exoayvjx4wFo3749f/75JwDvv/8+vXv3pkaNGmzYsIHnnnsOgMOHD9OiRQtq1KhBvXr1aN26dYq36UhISGDXrl3XXNCwceNGevfuTfPmzalXrx6DBg2iSZMmPm+Tt4KCgvjggw+4/fbbqVKlCnfddddVo5RJ68CT7t27U7VqVTp27Mjo0aO54YYbUp23cOFC2rdvnzEbJZKEsfba3/pFUlO3bl27Zo2Sj9RcjLvEmu8/oPqWN8htY5l+Qx+qdX+O8DJFnQ4tx9m2bZtf3H5i/fr1vP3223zxxRdOh+KzzZs3M2HCBN5+++2r3o+KimLcuHGEhYU5FFnm6tatG6+++qrH7XXX54wxa621dTMjPsk+NLImkoGCcwXR5M7HyPXoav4s3pSepz7BjruNl8d/w86/rj1vR7K/WrVq0aJFi3S/KW5mCg8PvyZRA9i9ezeVKlVyIKLMFxsbS5cuXXJMYirO0sia+EQjaz6wlvMbv8POHk7u2Gg+ju/IH9UGM7RNdcoWzed0dNmev4ysSc6hkTVJL7oaVCSzGEO+mt3h1hZcnP00Q7d8w+7fV/GPzQ9TqW5LHr2tEjeG6KoyERG5mn4GFcls+YoQfOc46P0t5QoFMCXXi4Ste5m2b87l5VlbOXH2otMRioiIH1GyJuKUSq0IHLqCgPoD6BM4jwV5n2HHshlEvbGQt37cTnRMnNMRZjs67UMyi/qapCclayJOCi4I7d+EfvMoUqggn+d+lU8Kf8pnCzYQ9cZCRi/cxfnYS05HmS3kyZOHEydO6EtUMpy1lhMnTuhmuZJudIGB+EQXGGSAuAvwy+uw9D3i8hRhbIHBvHkgjGIFghnS4hbuaVCW4KBAp6PMsuLi4jh48OCVJwiIZKQ8efIQGhpKrly5rnpfFxiIL5SsiU+UrGWgPzfAjKFw5Df+vrkdz164n3n7LaVC8vBoy0r0qBNKUKAGxUWyIiVr4gsd8UX8TamaMGAhtPw3NxxcwEfRg/ixxUFKFAzmme9+o9XbvzB9wyESEvSPlohITqBkTcQfBeaCyH/AwCVQvDK3Ln+KaSFv82X3m8iTK5BhX2+g/ajFLN113OlIRUQkgylZE/FnxW+FfnOh3ZuYP1bQ9McOzGm4jVG9IjgXe4nen6zkwc9Ws+fYWacjFRGRDKJz1sQnOmfNAaf+gJmPwe6foUxDLt7xHp9uz8UHC3ZxIS6e+xuVY1jLSoTky5VqUSLiDJ2zJr7QyJpIVlG4LNz7LXT5EI79TvC4KAYGfM/Cx5twZ90yTFy2l2YjF/LZsn3ExSc4Ha2IiKQTjayJTzSy5rAzf8HcJ2HrdLixOnQezVZbnpdnb2XZ7hNULFGAf95RhRZhJZyOVESS0Mia+EIja5ImxpiOxpix0dHRToeSsxUsCXd9Dnd94Urcxrag6ta3mdQ3gnH31yU+wdLv09X0mbCKnX+dcTpaERG5DhpZE59oZM2PxPwNPzwPG76EohWh0wfElm7A58v38d7POzkfG0/vBmV5rNWtFMmf2+loRXI0jayJL5SsiU+UrPmh3Qtg5jDXhQj1BkCrFzh5KZh35+9g0so/yJc7kGEtK3F/o3LkDtKguogTlKyJL5SsiU+UrPmpi2dhwUuw8mMICYWO70LFVuz86wwvz97GLzuOUb5Yfp5tV5nWVUtijHE6YpEcRcma+EL/XotkJ8EFoN3r8MAPkCsvfNkdpg2kUsE4PnugPp/2q0dggOGhL9bS+5OVbDt82umIRUQkFRpZE59oZC0LiLsAv74JS9+FvDdA+5FQrQtx8Ql8tfIP3pm/g9MxcfSsV4YnWodRvGCw0xGLZHsaWRNfKFkTnyhZy0IOb3I9GP7wRqjS0ZW0FbyR6PNxjFqwk8+W7SNPrkCGtKhIvyblyJMr0OmIRbItJWviCyVr4hMla1lM/CVY/gEsehWCguH2V6BmbzCGPcfO8sqcbczfdpQyRfLybLsqtAu/UeeziWQAJWviCyVr4hMla1nU8V0w4xH4YxlUaA4d34MbygGwZOdxXp69ld+PnKF+uSL8q0NVqoeGOBquSHajZE18oWRNfKJkLQtLSIC1E+CnF8AmQMsXoP4ACAgkPsEyZfUB3vpxOyfPx9KtVihPtQ2jZKE8Tkctki0oWRNfKFkTnyhZywZOHYBZj8OunyC0PnT+AIqHAXD6QhyjF+7i0yX7CAwwDGp+CwMiK5A3t85nE7keStbEF0rWxCdK1rIJa2HTNzDvaYg9B82egiaPQWAuAP44cZ5X525j7uYjlArJw9PtKtMpopTOZxPxkZI18YWSNfGJkrVs5uwxmPsUbPkOSoa7RtlK1boye8WeE7w0aytb/jxNzTKFebNHDSqVLOhgwCJZk5I18YVuiisiUKA43Pkp9PoKzh2HcS3hp39DXAwADSsUZcbQprzRowYH/z5PtzHLWLLzuMNBi4jkDErWROR/Kt8BQ1ZCrd6w9D34sAnsWwpAYIDhrrplmDG0KaVvyEvfT1fx9ao/HA5YRCT7U7ImIlfLWxg6vQ/3T4eESzCxPcx6Ai64Hk1VqnBe/m9gI5pULMYz3/3G6/N+JyFBp1OIiGQUJWsi4l6F5jB4OTQcAmsmwJhGsONHAArmycX4PnW5p0FZPly0m0cmr+dCXLyz8YqIZFNK1kTEs9z5oe0r0P8n10Piv7oTvnsIzp0gKDCAEV3C+Wf7KszZfJi7x63g+NmLTkcsIpLtKFkTkdSVqQcP/wrNnoHN38Lo+rD5WwwwIKoCH/auzbbDp+k6Zim7jp5xOloRkWxFyZqIeCcoGFo860raCpeBqQ/A1/fA6cO0Db+Jrx9qRExsPN3GLGPZbl0pKiKSXpSsiUjalKwG/edDm5dh9wIY3QDWfkbN0BCmDW5CyUJ5uH/8KqauPeh0pCIi2YKSNRFJu8AgaPwIDFoGN1aHmY/C550owxGmDmpMwwpFGf5/G3n7x+3oxtsiItdHyZqI+K7oLdBnJnR4F/7cAGMaE7L+Yz7tU5uedcswasEuHpuyQVeKiohcByVrInJ9AgKgbj8YvAIqNIMf/0muibfzWmQgT7UNY/qGP7lv/EpOnot1OlIRkSxJyZqIpI+Q0nD319B9PPy9D/NxMwYzldG9qrHxYDTdxixl7/FzTkcpIpLlKFkTkfRjDFTvAUNWQbUusOhV7lh2N9O7BHP6wiW6jlnKqr0nnY5SRCRLUbImIukvfzHo/olrpC3mFFVmd2dhjZ+4KZ/l3k9W8v36Q05HKCKSZShZE5GME9YOhqyA2n0IWf8xs4Ke4t4b9/PYlA28N3+nrhQVEfGCkjURyVh5QqDju9BnFoHG8O8TTzP5xsl8Mn8D//i/jcReSnA6QhERv6ZkTUQyR/lI133ZGj9Cw+jZLC/0LNEbZnD/hJWcOq8rRUVEPFGyJiKZJ3c+aPMy5sH5FChcgvG53+Leg/+l3+i57D+hK0VFRNxRsiYima90HXhoETR/jvZBq5lwbggff/A6a/edcDoyERG/o2RNRJwRlBuaP03AwMXkLVmRV+x7nJ7QnZ9WrHM6MhERv6JkTUScVaIKeR7+mXMtXqJR4FYazm3PokmvYxP0iCoREVCyJiL+ICCQ/M0exQxazp/5q9B85yvsGXkbcUd3Oh2ZiIjjlKyJiN8ILnELtw7/mR9v+SfFz+3AjmlMzKK3If6S06GJiDhGyZoAYIypYIwZb4yZ6nQskrOZgADa3PcUv7aZzS8JNci76D/EfnwbHNnsdGgiIo5QspYNGGMmGGOOGmM2J3u/rTFmuzFmlzHmmZTKsNbusdb2z9hIRbzXoUlt8t//NcN5jLNH92I/bgYLRsCli06HJiKSqZSsZQ8TgbZJ3zDGBAKjgXZAVeBuY0xVY0x1Y8ysZK8SmR+ySOoaVyzOwMFPcl/e95ke3wh+fQM+joIDq50OTUQk0yhZywastb8CJ5O9XR/YlThiFgt8DXS21v5mre2Q7HU004MW8VLFEgX4bEg7PrvxWfrGPkX0qb+x41vDvGchVjfSFZHsT8la9lUaOJBk+mDie24ZY4oaYz4CahljnvWwzEPGmDXGmDXHjh1L32hFUlCsQDCTBzSkXMMuNDkzglm528GKMTCmEexZ5HR4IiIZSsla9mXcvGc9LWytPWGtHWitvcVa+6qHZcZaa+taa+sWL1483QIV8UaeXIG82Kkao/pG8WL8A/S+9AKnY4HPO8P0oRBzyukQRUQyhJI1P2GMeSGdizwIlEkyHQr8mc6fIZLpbqtckrmPRRJQvgn1Tv6HuSG9sBu+gtENYNssp8MTEUl3Stb8xwvGmNeNMeOMMYOMMTdcZ3mrgUrGmPLGmNxAL2DG9Ycp4rwSBfPwWb/6PHlHBMOOd6FP4KuczXUDTOkN3/SBszoNU0SyDyVr/sMCF4AfcI2ILTPGRHizojFmMrAcCDPGHDTG9LfWXgKGJpa3DfjGWrslY0IXyXwBAYYHIyvw3eDGHMobRq0jz/JL6CDs9jkwuj5s/Bqsx1/+RUSyDGN1MPMLxpgt1tpqSaZvBT6y1t7mYFge1a1b165Zs8bpMEQAOB97iZdmbWPyqj9oVzKad/J+Qp4ja6Fia+jwDhQuk3ohIpnAGLPWWlvX6Tgka9HImv84boypc3nCWrsD8Luz+I0xHY0xY6Ojo50OReSKfLmDeLVbdT66tw7LzxSj7p/DWVftWez+ZTCmIawaBwkJTocpIuITjaz5icSfPL8G1gK/ATWAgtbaTo4G5oFG1sRfHY6O4fEpG1ix5yT3VTH8235Mrn2LoGxj6PQ+FKvodIiSg2lkTXyhkTU/Ya3dCNQEJie+tRC427GARLKom0LyMunBhjzVNozJ26H54UfZ0+RNOLoFPmwMS97Rg+FFJEvRyJr4RCNrkhVsOHCKYV+v58DJ8zzd9AYGnBlDwO8z4aYI6PQB3FTD6RAlh9HImvhCI2sikm3VLFOY2Y9G0rVWKK8u/pseJwdxvP04OH0YxjaHn/8LcRecDlNEJEVK1kQkWysQHMRbd0Uw6u5a7PzrLC1mhzCn2fdQoycsfgs+joQ/VjodpoiIR0rWRCRH6BRRijnDIrn1xoIM/m4fT8Q+zPm7voG4GJhwO8x5Ci6edTpMEZFrKFmTNNGtOyQrK1MkH1MeasiwlpX4fsMh2s7KzaZO86D+AFg11vVg+F0/Ox2miMhVlKxJmlhrZ1prHwoJCXE6FBGfBAUG8HjrW5nycCPiEyzdxm9idN6Hie87B4KC4ctu8P1gOH/S6VBFRAAlayKSQ9UrV4Q5wyJpG34jb/6wnXt+MPx5908Q+Q/Xo6rGNIStepyuiDhPyZqI5FgheXPx/t21GHlnBL8diqbd6NXMKzkAHloEBUrCN/fBlPvgzF9OhyoiOZiSNRHJ0Ywx9KgTypxHI7m5aD4GfrmOZ5fD+T4/QssXYMcPrgfDr5+kB8OLiCOUrImIAOWK5WfqwMYMan4LX68+QIcxK9lcoT8MWgolqsD0wa7z2f7e73SoIpLDKFkTEUmUOyiAp9tWZtKDDTh/MZ6uY5byybZAEvrMhvYj4cAq1xWjKz/Wg+FFJNMoWZM00a07JCdofEsx5g6LpEVYCV6evY0+E9dwtPJ9MHg5lG0Ic5+CT9vBsR1OhyoiOYCeDSo+0bNBJSew1jJ51QH+O2sLuQMDeP6OqtxZpzRm0xSY9wzEnYdmT0OTYRCYy+lwJQvQs0HFFxpZExHxwBjDPQ3KMndYFJVvKsRT327i3gmr+KNMZxi6GsLaw4KXYFwL+HOD0+GKSDalZE1EJBXli+Xn6wENGdE1nI0Homnz7i98sv4s8T0mQs8v4exRGHcbzH/R9fgqEZF0pGRNRMQLAQGG3g1u5qcnomhasRgvz95Gtw+X8fsNzWDISqh5Nyx5Bz5qCvuXOR2uiGQjStZERNLgppC8jLu/Lu/fXYuDJ8/TYdQS3l58lIt3jIL7vof4WNfFB7OHw8UzTocrItmAkjURkTQyxtAxohTzn2hGp4hSjFqwiztGLWFtUAQMWg4NBsHqT1y3+dg53+lwRSSLU7ImIuKjG/Ln5u2eNZnYrx4xsfH0+Gg5L/6wn3O3vQz9f4Rc+WBSd5g2UA+GFxGfKVkTEblOzcNK8MPjUdzf8GY+W76PNu/8yi8x5WHgYoh6En77P9cjq7ZM0yOrRCTNlKxJmuimuCLuFQgO4j+dw5k6sBF5cgXQZ8IqnvhuG383eMr1YPhCpeH/+sKUe+HMEafDFZEsRDfFFZ/oprginl2Ii2f0wl18uGg3hfPl4oWO1egQXhyzYgwsfAUCg+H2EVDrXjDG6XAlE+mmuOILjayJiKSzPLkC+UebMGY+0pRShfPyyOT1DPhyI0fCH4aBS+HGcJgxFL7oAn/vczpcEfFzStZERDJIlZsK8d2gxvyzfRWW7DpG67d/4avduUm4fybc8TYcXOu6YnTFh5AQ73S4IuKnlKyJiGSgoMAABkRV4IfHoggvHcJz037j7k9Wsbd8LxiyAm5u4nrO6ITb4ejvTocrIn5IyZqISCa4uWh+vhrQgNe7V2fr4dO0ffdXPtpwkUu9pkC3cXBiN3wcCb+8AZdinQ5XRPyIkjURkUxijKFnvbLMf6IZzW4tzmtzf6fLh8vYUux2GLIKKneAhSNcD4Y/tM7pcEXETyhZExHJZCUL5eHj++owpndtjkRfpNMHS3ljyQkudPkEen0F50/AJy3hx3/pwfAiomRNRMQJxhjaV7+J+U9E0bVWacYs2k379xazKrgRDF4Bte6DZaPgw8awb4nT4YqIg5SsiYg4qHC+3Iy8M4LPH6hPbHwCd328nH/9cJAzbd6C+2eATYCJd8Csx+HCaafDFREHKFkTEfEDUbcW54fHonigSXm+XLmfNu/8yoLYyjBoGTQaCmsnwpiGsOMHp0MVkUymZE3SRI+bEsk4+YOD+HfHqnw7qDEF8wTxwMQ1DJ26nR01n4H+P0FwQfjqLvh2AJw74XS4IpJJ9Lgp8YkeNyWSsWIvJTBmkeuRVRcvJRB1a3EebFSKyCOfYxa/BXlCoN0bEN5dj6zKQvS4KfGFkjXxiZI1kcxx8lwsX63cz2fL93PszEVuLVmAx6vHcfvulwk4vB7C2sMdb0GhUk6HKl5Qsia+ULImPlGyJpK5Ll6KZ9bGw4xfspeth09TIl8gb5ZdRtTBsZjAXNDmJajdR6Nsfk7JmvhCyZr4RMmaiDOstSzfc4IJS/Yyf9tRKgYe5aOQz6h4fj2Ui4ROo6BIBafDFA+UrIkvgpwOQEREvGeMofEtxWh8SzH2HDvLp0v30XFtCbok/My/9k8meHQjzG3PE9BoMAQEOh2uiKQDjayJTzSyJuI/Tp2P5atVfzB36VqGXfiQVoHrOR5SnYJ3fUhw6epOhydJaGRNfKFbd4iIZHGF8+VmcPOKfPv0nZzt+iVvFHgSc2ofAeOaseyTf3D0lG6mK5KVaWRNfKKRNRH/Za1l3badXJrzDA3O/swOG8qcCs/TpvUdVC1VyOnwcjSNrIkvlKyJT5SsiWQNf62eRp4fn6RA7HEmxLdjaZmHuS+qCi3CShAQoCtHM5uSNfGFkjXxiZI1kSzkQjQX5/6L4I2fcZCSPBn7IH8VqU+/puXpUTuUvLl1IUJmUbImvlCyJj5RsiaSBe1djJ3xCObvvfwQfDvDo+8kMF8I99QvS5/G5ShZKI/TEWZ7StbEF0rWxCdK1kSyqNjzsOgV7PLRxOUtzrhCjzByfwWCAgwdapSif9PyhJcOcTrKbEvJmvhCyZqkiTGmI9CxYsWKA3bu3Ol0OCLiq0NrYfojcHQL5yp1YUzeB5m44RznYuNpUL4ID0ZWoGVlndeW3pSsiS+UrIlPNLImkg1cioWl78Ivb0BwQc63GsGksw2YuHw/h07FcEeNmxh9T22no8xWlKyJL5SsiU+UrIlkI0e3wfShcGgNVLqdS+1GMu9gEIXz5qZppWJOR5etKFkTX+imuCIiOV2JKtD/R7j9Vdi3mKCPGtPh4lya3lLE6chEBCVrIiICrueINhoMg5ZB6dow+wn4rCOc2O10ZCI5npI1ERH5nyLl4f7p0Ol9OPIbfNgYlr4H8Zecjkwkx1KyJiIiVzMGat8PQ1bCLS3hp3/D+FZwZLPTkYnkSErWRETEvUI3Qa9JcOdEiD4IY5vBghFw6aLTkYnkKErWRETEM2OgWlcYsgrCe8Cvb8BHkXBgldORieQYStZERCR1+YpAt4+h91SIPQfj28DcZ1x/i0iGUrImIiLeq9QaBi+Hev1h5YcwpiHsXuh0VCLZmpI1ERFJmzyF4I63oO8cCMgFX3SB6UMg5pTTkYlkS0rWRETEN+WawKCl0PRx2DAZRjeAbbOcjkok21GyJiIivsuVF1q9CAN+hvzFYUpv+KYPnD3qdGQi2YaSNRERuX6lasFDC+G2f8H2OTC6Pmz8GvT8aZHrpmRNRETSR2AuiBoOA5dAsVth2sMwqQecOuB0ZCJZmpI1SRNjTEdjzNjo6GinQxERf1U8DPrNg3ZvwP7lritGV42DhASnIxPJkpSsSZpYa2daax8KCQlxOhQR8WcBAdDgYddtPkLrwZzhMLE9HN/pdGQiWY6SNRERyTg33Az3TYPOY+DoVviwCSx+Ww+GF0kDJWsiIpKxjIFavWHIari1Dfz8H/jkNji8yenIRLIEJWsiIpI5CpaEnl/CXZ/D6cMwtjn8/F+Iu+B0ZCJ+TcmaiIhkrqqdYchKiOgFi9+CjyPhj5VORyXit5SsiYhI5stXBLqMgXu/c42sTbgd5jwFF886HZmI31GyJiIizqnY0nXFaP2HYNVYGNMIdv3sdFQifkXJmoiIOCu4ALR/Ax6YB0HB8GU3+H4wnD/pdGQifkHJmoiI+IeyDV1PP4j8h+tRVaMbwNbpTkcl4jglayIi4j9y5YGW/4aHFkHBG+Gb+2HKfXDmL6cjE3GMkjUREfE/N9WAAQug1Yuw4wfXg+HXT9KD4SVHUrImIiL+KTAXNH0cBi2FElVg+mDX+Wx/73c6MpFMpWRNRET8W7FK0HcOtB8JB1a5rhhd+bEeDC85hpI1ERHxfwEBUH+A6zYfNzeCuU/Bp23h2HanIxPJcErWREQk6yhcFnpPha4fw/Ed8FFT+HUkxMc5HZlIhlGyJiIiWYsxrkdVDVkFle+ABS/BuBbw5wanIxPJEErWREQkaypQAu6cCD0nwdmjMO42mP8ixMU4HZlIulKyJiIiWVuVDq4Hw9e8B5a84/ppdP8yp6MSSTdK1kREJOvLewN0/gDu+x7iY+HTdjD7H3DxjNORiVw3JWsiIpJ93NICBq+AhoNh9XgY3RB2/uR0VCLXRcmaiIhkL7nzQ9tXof+ProfET+oB3z2sB8NLlqVkTUREsqcy9eHhXyHqKdg8FT6oB5u/0yOrJMtRsiZpYozpaIwZGx0d7XQoIiKpCwqG2/4JD/0CIaEwtR9MuRdOH3Y6MhGvKVmTNLHWzrTWPhQSEuJ0KCIi3rsxHB78GVr/F3bNh9ENYN3nGmWTLEHJmoiI5AyBQdBkGAxaBjdWhxmPwOed4eRepyMTSZGSNRERyVmK3gJ9ZkKHd+DQOviwMSwfAwnxTkcm4paSNRERyXkCAqDuA66b6ZaLhB+ehfFt4Og2pyMTuYaSNRERyblCSsM9U6DbJ3ByD3wUCb+8AZdinY5M5AolayIikrMZAzXuhKGroWonWDgCxjaHQ2udjkwEULImIiLikr8Y9JgAvSZDzEn4pBX8+C+IPe90ZJLDKVkTERFJqnJ717lste+HZaPgoyawd7HTUUkOpmRNREQkuTwh0PE911Wj1sJnHWDmY3BBNwSXzKdkTURExJPyUa77sjUaCus+04PhxRFK1kRERFKSOx/cPgL6z4e8heH0IacjkhwmyOkAREREsoTQOq5njAbmcjoSyWGUrImIiHgrKLfTEUgOpJ9BRURERPyYkjURERERP6ZkTURERMSPKVkTERER8WNK1kRERET8mJI1ERERET+mZE1ERETEjxlrrdMxSBZkjDkG7HczKwRI/vC85O8VA45nUGipcRdfZpTj7fKpLZfSfE/zvGkTcK5dnGqTtKyT3u3ibVtpX/F9OX/dV2621hb3cV3Jqay1eumVbi9gbGrvAWv8Kb7MKMfb5VNbLqX5nuZ50yZOtotTbeJku3jbVtpXMq9N0tJWTraLXjnzpZ9BJb3N9PI9p6RXLGktx9vlU1supfme5qlNrn+d9G6XtLSVU7SvePc5IhlOP4NKpjPGrLHW1nU6Drma2sX/qE38k9pFMptG1sQJY50OQNxSu/gftYl/UrtIptLImoiIiIgf08iaiIiIiB9TsiYiIiLix5SsiYiIiPgxJWviOGNMBWPMeGPMVKdjkf8xxnQxxowzxkw3xrRxOh4BY0wVY8xHxpipxphBTscjLsaY/MaYtcaYDk7HItmTkjXJEMaYCcaYo8aYzcneb2uM2W6M2WWMeQbAWrvHWtvfmUhzljS2y/fW2gFAX6CnA+HmCGlsk23W2oHAXYBuHZFB0tImiZ4GvsncKCUnUbImGWUi0DbpG8aYQGA00A6oCtxtjKma+aHlaBNJe7s8nzhfMsZE0tAmxphOwBLg58wNM0eZiJdtYoxpBWwF/srsICXnULImGcJa+ytwMtnb9YFdiSNpscDXQOdMDy4HS0u7GJfXgbnW2nWZHWtOkdZ9xVo7w1rbGOiduZHmHGlskxZAQ+AeYIAxRt+rku6CnA5AcpTSwIEk0weBBsaYosAIoJYx5llr7auORJdzuW0X4BGgFRBijKlorf3IieByKE/7SnOgGxAMzMn8sHI0t21irR0KYIzpCxy31iY4EJtkc0rWJDMZN+9Za+0JYGBmByNXeGqXUcCozA5GAM9tsghYlLmhSCK3bXLlD2snZl4oktNouFYy00GgTJLpUOBPh2KR/1G7+B+1if9Rm4hjlKxJZloNVDLGlDfG5AZ6ATMcjknULv5IbeJ/1CbiGCVrkiGMMZOB5UCYMeagMaa/tfYSMBT4AdgGfGOt3eJknDmN2sX/qE38j9pE/I0e5C4iIiLixzSyJiIiIuLHlKyJiIiI+DElayIiIiJ+TMmaiIiIiB9TsiYiIiLix5SsiYiIiPgxJWsiIiIifkzJmoiIiIgfU7ImItmGMWaaMeZlY8xiY8wRY0wrp2MSEbleStZEJDsJB05ZayOBwUBvh+MREbluStZEJFswxuQDQoB3Et8KAk45FpCISDpRsiYi2UU1YK21Nj5xugaw2cF4RETShZI1EckuwoENSaZrAJucCUVEJP0oWROR7KI6Vydr4WhkTUSyAWOtdToGEREREfFAI2siIiIifkzJmoiIiIgfU7ImIiIi4seUrImIiIj4MSVrIiIiIn5MyZqIiIiIH1OyJiIiIuLHlKyJiIiI+LH/B405I2YNW02wAAAAAElFTkSuQmCC\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": 10,
   "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 10000 \td = 1 \tV_d = 2.000 \tdirect-sampled V_d = 1.000\n",
      "N = 10000 \td = 2 \tV_d = 3.142 \tdirect-sampled V_d = 0.788\n",
      "N = 10000 \td = 3 \tV_d = 4.189 \tdirect-sampled V_d = 0.530\n",
      "N = 10000 \td = 4 \tV_d = 4.935 \tdirect-sampled V_d = 0.304\n",
      "N = 10000 \td = 5 \tV_d = 5.264 \tdirect-sampled V_d = 0.165\n",
      "N = 10000 \td = 6 \tV_d = 5.168 \tdirect-sampled V_d = 0.083\n",
      "N = 10000 \td = 7 \tV_d = 4.725 \tdirect-sampled V_d = 0.039\n"
     ]
    }
   ],
   "source": [
    "import scipy.special\n",
    "\n",
    "def volume_d_ball(d):\n",
    "    \"\"\"Returns the volume of a d-dimensional unit ball.\"\"\"\n",
    "    return np.pi**(d/2)/scipy.special.gamma(d/2 + 1)\n",
    "\n",
    "N = 10000\n",
    "for d in range(1, 8):\n",
    "    print(\"N = {}\".format(N),\n",
    "          \"\\td = {}\".format(d),\n",
    "          \"\\tV_d = {:.3f}\".format(volume_d_ball(d)),\n",
    "          \"\\tdirect-sampled V_d = {:5.3f}\".format(number_of_hits_in_d_dimensions(N, d)/N)\n",
    "         )\n",
    "\n",
    "# TODO: Results not yet as expected."
   ]
  },
  {
   "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
}