{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7711bd97",
   "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",
    "* For some exercises test cases have been provided in a separate cell in the form of `assert` statements. When run, a successful test will give no output, whereas a failed test will display an error message.\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": null,
   "id": "e39c230e",
   "metadata": {},
   "outputs": [],
   "source": [
    "NAME = \"Kees van Kempen\"\n",
    "NAMES_OF_COLLABORATORS = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e7caa84",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c4606a8",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "1ba146ece4effd3e8dcf33fe63ec69dc",
     "grade": false,
     "grade_id": "cell-455926422ebe4f85",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "**Exercise sheet 4**\n",
    "\n",
    "Code from the lectures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6d3bb610",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "ac562d84531ae19e1d920c10a8861b8e",
     "grade": false,
     "grade_id": "cell-6b2d0465c2abaaa6",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "import networkx as nx\n",
    "\n",
    "rng = np.random.default_rng()\n",
    "%matplotlib inline\n",
    "\n",
    "def draw_transition_graph(P):\n",
    "    # construct a directed graph directly from the matrix\n",
    "    graph = nx.DiGraph(P)        \n",
    "    # draw it in such a way that edges in both directions are visible and have appropriate width\n",
    "    nx.draw_networkx(graph,connectionstyle='arc3, rad = 0.15',width=[6*P[u,v] for u,v in graph.edges()])\n",
    "\n",
    "def sample_next(P,current):\n",
    "    return rng.choice(len(P),p=P[current])\n",
    "\n",
    "def sample_chain(P,start,n):\n",
    "    chain = [start]\n",
    "    for _ in range(n):\n",
    "        chain.append(sample_next(P,chain[-1]))\n",
    "    return chain\n",
    "\n",
    "def stationary_distributions(P):\n",
    "    eigenvalues, eigenvectors = np.linalg.eig(np.transpose(P))\n",
    "    # make list of normalized eigenvectors for which the eigenvalue is very close to 1\n",
    "    return [eigenvectors[:,i]/np.sum(eigenvectors[:,i]) for i in range(len(eigenvalues)) \n",
    "              if np.abs(eigenvalues[i]-1) < 1e-10]\n",
    "\n",
    "def markov_sample_mean(P,start,function,n):\n",
    "    total = 0\n",
    "    state = start\n",
    "    for _ in range(n):\n",
    "        state = sample_next(P,state)\n",
    "        total += function[state]\n",
    "    return total/n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f7814d8",
   "metadata": {},
   "source": [
    "## Markov Chain on a graph\n",
    "\n",
    "**(50 points)**\n",
    "\n",
    "The goal of this exercise is to use Metropolis-Hastings to sample a uniform vertex in a (finite, undirected) connected graph $G$. More precisely, the state space $\\Gamma = \\{0,\\ldots,n-1\\}$ is the set of vertices of a graph and the desired probability mass function is $\\pi(x) = 1/n$ for $x\\in\\Gamma$. The set of edges is denoted $E = \\{ \\{x_1,y_1\\}, \\ldots,\\{x_k,y_k\\}\\}$, $x_i,y_i\\in\\Gamma$, and we assume that there are no edges connecting a vertex with itself ($x_i\\neq y_i$) and there is at most one edge between any pair of vertices. The **neighbors** of a vertex $x$ are the vertices $y\\neq x$ such that $\\{x,y\\}\\in E$. The **degree** $d_x$ of a vertex $x$ is its number of neighbours. An example is the following graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "071e5617",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABHp0lEQVR4nO3dd1iT198G8DshgaCAyBIRHBUBUYYbRRG1ddVZ9adV3NaBq4qr4mqtWrXW1jrqbC1o5S1W3FtwISgORBAQB1sEFZkJGc/7h8WKBGQkeTK+n+vyuiwJya1Fbs55zjkPh2EYBoQQQoiO4LIdgBBCCFElKj5CCCE6hYqPEEKITqHiI4QQolOo+AghhOgUKj5CCCE6hYqPEEKITqHiI4QQolOo+AghhOgUKj5CCCE6hYqPEEKITqHiI4QQolOo+AghhOgUKj5CCCE6hYqPEEKITqHiI4QQolOo+AghhOgUKj5CCCE6hYqPEEKITqHiI4QQolOo+AghhOgUKj5CCCE6hcd2AEIAIKdAhODbaYh/noc8oQQmAh6crE0wop0tzI0M2I5HCNEiHIZhGLZDEN0VnZqLbWFJuJyYDQAQSWTvHhPwuGAAeDtawre7PdzsTNkJSQjRKlR8hDWBEc+w5lQ8hBIpKvsq5HAAAU8P/v2d4OPRVGX5CCHaiaY6CSvelt5DFItlH30uwwDFYinWnHoIAFR+hJBaoREfUbno1FyM2h2BYrG0zMelxfl4eeoXCJ/dBdfQBPW7j0fdVt5lnmPI10PQVA+42pqqLjAhRKvQqk6ictvCkiCUSMt9/NW5HeDo8WE7OxAWAxfg5bntKMlOLvMcoUSK7WFJqopKCNFCVHxEpXIKRLicmF3ump6sRIiihHCYevmAq28IgV0r1LHvhMLY0DLPYxggNCEbLwtEKkxNCNEmVHxEpYJvp8n9uORVOjhcLvhmjd59jG/VDOIPRnwAwAEQfEf+6xBCyMfQ4haiUvHP88psWSglExeDY1CnzMe4BnUgKyku91yhRIZ7T1+gxMMO+vr6SstKVIv2chJVoeIjKpUnlMj9OJdvCEZUtuQYURG4+oZyn3/87EXsmeIFfX191K9fv0a/DAzom6k6qHwv53NsvpBIezmJQlHxEZUyEcj/kuOZNQIjk0L8Kv3ddGfJi6fgWzaR+/wRgwfgpwP+KCwsxOvXryv8lZCQUOFjPB6vxqUpEAiU9nekSz62l1P4bwmei8vClcQc2stJFIKKj6iUk7UJDHjPy013cvUFqOPYGblXD8C83xyUvHiCoqRIWPtsLPcaAh4XTg2NweFwYGRkBCMjI9jZ2VUrB8MwKCoqqrQ0Hz16VOFjXC63xqVpaCh/FKtraC8nYQvt4yMqlVMgguf6S3Kv81VlHx8AGPC4CF/ck7XrPgzDoLi4uNLSrOwXgEqL0dTUtMLH6tSpAw6Hw8qfW5Eq2suZc/xHCJ9FQyYWQq9ufZh4DIOxW58yz6G9nKS2qPiISkkkEvRceQjJknrgcKu/qJjDAfo4N8BvPu2VkE41alOaMpms0mKs7FfdunXVpjSnBkTh/MOsctObJdnJ4Ne3AYfHh/hlKp4f/AZWI1bBwNr+3XO04WuAsIumOonKPHr0CD4+PhA0coSg1SiIJNX/mUvA04Ovt/3Hn6jGDA0NYWhoCBsbm2p/rlAorLQYU1JSEB0dLfcxiURS49I0MjJSWGlWtJcTAPTLXNPlgAMOJK8zyxTf+3s5abUnqQkqPqJ0DMNg165dWLZsGVauXImZM2fiQGRyla/vlDLkc+Hf30mnp7gEAgEaNmyIhg0bVvtzRSIRXr9+jdzcXLnFmJ6ejgcPHsh9TCQS1bg0jY2Ny5RmRXs5S708ux2FMRfBSETQb9Achs3Lj+xK93JO82pe7b8HQmiqkyhVVlYWpkyZgoyMDAQGBqJly5bvHqvq3RnAyGCoz6cVfSwqKSmpsDA/9ksoFKJevXrvirDYfSTyzZ0qfT9GJoUoPR7ClBjU8xgOjl75n9GHujfC5pHuSvoTE21GIz6iNEePHsX06dMxefJkHD58uNxmcx+PpnC1NcX2sCSEJmSDg/+WrwP/3Y9PnBwNvyHtqfRYpK+vDysrK1hZWVX7c8VicZnS/P5KDmJeVf45HK4eBHatUBgbivy7p2DSflC55+QJxdXOQghAxUeUID8/H/PmzUNoaCiCg4Ph6elZ4XNdbU3xm097vCwQIfhOGuIz85EnFMNEwIdTQ2MMb2uLi6fy8eMPyzB56GdqsziDVB2fz4elpSUsLS0BAM2f3UXMq4yqfbJMBsnrTLkPmQj4iopIdAyd1UkUKjw8HG3atAHDMLh3716lpfc+cyMDTPNqjs0j3bF3fAdsHumOaV7NYW5kgOHDh0MsFuPIkSNKTk9U4e1ezvLfeqSFuSiMuwxZSTEYmRTFT26j8OFlCJq4lXtu6V5OQmqCrvERhRCLxfj222+xZ88e/PbbbxgyZIhCX//06dPw8/NDTEwM9PT0FPraRLUq2sspLXqD7CPrUPLiKcDIwKtnBeN2A2Hs3rfca7C9l5NoNio+Umvx8fHw8fFBgwYNsHfvXlhbWyv8PRiGQffu3TF58mSMHz9e4a9PVGtqQBTOx2WhJt98GJkMViWZOOLXH7a2tgrPRrQfTXWSGmMYBtu2bUO3bt0wZcoUnDhxQimlBwAcDgdr167FypUrIRLRvfg0Xb/GXDCSkhp9rqE+D20E2XBzc8M333yD3NxcxYYjWo+Kj9RIRkYG+vXrh/379+P69euYPn260heedO3aFa1atcKuXbuU+j5EuY4fP46pw/tgoJ0YhvzqfQsy5HOx7POW2LXOH9HR0Xjx4gUcHBzw008/0Q9EpMqo+Ei1HT58GG3btoWHhweuX78OBwcHlb33mjVrsHbtWhQWFqrsPYliMAyDDRs2YPr06Thx4gR+nfM/+PdvCUO+Hj72MxOH8/aMTv/+Ld9ta7G1tcXevXsRGhqKsLAwODo6IjAwEDJZ1Q9FILqJrvGRKsvLy8OcOXNw/fp1BAYGolOnTqzkGDVqFFxdXbF06VJW3p9Un1AoxNSpUxEbG4ujR4+WuTZ3Py33o3s5ezhawtfbvtJTe65cuYLFixejuLgY69evR+/evWn7C5GLio9UydWrVzFu3Dj07t0bmzZtgpGREWtZHj16hC5duiAhIQFmZmas5SBVk5mZiaFDh6JJkyb4/fffUadOHbnPq2wvZ1VXbzIMgyNHjuCbb76BnZ0d1q9fj3bt2inyj0O0ABUfqVRJSQlWrlyJP/74A7t27cLAgQPZjgQAmDp1KszMzPDDDz+wHYVU4vbt2xg6dCi++uorLFu2TGUjMLFYjL179+K7776Dt7c3vv/+e3zyyScqeW+i/ugaH6lQbGwsOnXqhLi4OERHR6tN6QHAihUrsHv3bmRmyj/Vg7Dv77//Rt++fbF582YsX75cpdOOfD4f06dPR2JiIpycnNChQwfMnTsX2dnZKstA1BcVHylHJpPhl19+gbe3N2bOnImQkJAandGoTLa2tpgwYQK+//57tqOQD8hkMqxcuRILFizA+fPnMWzYMNayGBkZYcWKFXj48CFkMhlatmyJNWvW0OIoHUdTnaSM9PR0TJgwAQUFBQgICIC9vfre+y4nJwdOTk64efMmTWOpicLCQowfPx6ZmZn4559/0KBBA7YjlZGUlAR/f39cu3YNK1euxKRJk8Dj0ZHFuoZGfOSdoKAgtG3bFt27d8fVq1fVuvQAwMLCArNnz8aqVavYjkIApKSkoGvXrjA2NsalS5fUrvQAwN7eHkFBQQgJCcGhQ4fg4uKCkJAQ0M//uoVGfAS5ubmYNWsWoqKiEBAQgA4dOrAdqcry8vLQokULXLx4Ea1bt2Y7js4KDw/H8OHD4efnh/nz52vENgKGYXDmzBksXrwYxsbG2LBhQ5UPVSeajUZ8Oi4sLAxubm6oV68e7ty5o1GlBwAmJiZYvHgxli1bxnYUnbV//34MGTIEe/bsgZ+fn0aUHvD2GLx+/frh7t27mDZtGkaPHo0hQ4bg4cOHbEcjysYQnSQUChk/Pz/GxsaGOXXqFNtxaqW4uJixtbVlbty4wXYUnSKRSBg/Pz+mefPmTGxsLNtxaq24uJjZuHEjY2FhwXz11VdMeno625GIktCITwfFxMSgQ4cOePr0KaKjo9GvXz+2I9WKQCDAihUr4O/vz3YUnfHmzRsMGjQId+/eRWRkJJydndmOVGsCgQALFixAYmIi6tevDxcXF/j7++PNmzdsRyMKRsWnQ2QyGTZt2oSePXti/vz5CA4OhoWFBduxFGLixIlITU3FhQsX2I6i9ZKSktC5c2c0a9YMZ86cgbm5OduRFKp+/fpYv3497t27h4yMDDg4OODnn3+mQ7C1CC1u0REpKSkYP348xGIxAgIC0KxZM7YjKVxQUBA2bdqEyMhIjbnOpGkuXbqE0aNHY+XKlZgxYwbbcVQiJiYGS5YswcOHD/H9999j1KhR4HJpzKDJ6P+elmMYBgcOHED79u3Rp08fXL58WStLDwBGjBgBsViMkJAQtqNope3bt2P06NH466+/dKb0AMDFxQUnT57Evn378Msvv6BDhw40s6DhaMSnxV6/fo0ZM2bg/v37OHDgANq0acN2JKU7deoUFi5ciPv370NPT4/tOFpBLBZj7ty5CAsLw/Hjx9G8eXO2I7GGYRgEBwdj6dKlaNasGdavX68T/660DY34tNSFCxfg6uqKBg0a4Pbt2zrzj7Nfv34wMzNDYGAg21G0wsuXL9GnTx8kJycjIiJCp0sPeLsFYsSIEYiLi8OQIUPQv39/+Pj44NmzZ2xHI9VAxadliouL8fXXX2PixInvpmYMDQ3ZjqUyHA4Ha9euxcqVK2kxQi3FxcWhU6dOaNeuHY4dOwYTExO2I6kNPp8PX19fJCYmwt7eHu3atcO8efOQk5PDdjRSBVR8WuTevXto3749MjMzER0djc8++4ztSKzo1q0bnJ2dsXv3brajaKxTp07B29sby5cvx8aNG2nauALGxsZYtWoV4uLiUFJSAicnJ6xbtw5FRUVsRyOVYXMTIVEMiUTC/PDDD4ylpSUTEBDAyGQytiOx7s6dO4y1tTVTUFDAdhSNIpPJmI0bNzINGzZkwsPD2Y6jcRITE5nhw4czjRo1Ynbv3s2IxWK2IxE5aHGLhnv27BnGjRsHLpeL/fv3o0mTJmxHUhujRo2Cq6srli5dynYUjSASiTBt2jRER0fj6NGjaNy4MduRNFZkZCQWL16M7OxsrFu3DgMHDqQtNmqEik9DMQyDP//8EwsWLMDixYsxb948mo76QGJiIjw9Pd+dxEEq9vz5c3zxxRewsbHB/v37UbduXbYjaTyGYXDq1CksWbIEpqam2LBhAzp37sx2LAIqPo308uVLTJs2DQkJCThw4ABcXV3ZjqS2vvrqK1hYWGDdunVsR1Fbd+/exZAhQzBx4kSsWLGCNmcrmFQqRUBAAFasWIH27dtj3bp1cHR0ZDuWTqOvcA1z9uxZuLm5oUmTJrh16xaV3kesWLECu3btQmZmJttR1NLhw4fRu3dv/Pjjj1i1ahWVnhLo6elhwoQJSEhIgIeHB7p27Yrp06fT1ySLaMSnIYqKirB48WIcPXoUf/zxB3r27Ml2JI0xf/58iEQibNu2je0oaoNhGKxevRp79uxBSEgI2rZty3YknfHq1SusW7cO+/btg6+vLxYuXEhbRVSMfrzTALdv30a7du3w6tUrREdHU+lV0zfffINDhw7hyZMnbEdRC0VFRRg5ciROnz6NmzdvUumpmJmZGTZu3Ig7d+4gJSUFDg4O+PXXX1FSUsJ2NJ1BxafGJBIJ1qxZg379+mHlypU4cOAALdKoAUtLS8yePRurVq1iOwrr0tLS0K1bNwgEAoSGhsLa2prtSDqrSZMm2L9/P86ePYtTp06hZcuWOHToEGQyGdvRtB5Ndaqpx48fY9y4cRAIBPjjjz9gZ2fHdiSNlpeXhxYtWuDixYto3bo123FYERERgWHDhmHu3LlYuHAhLa9XM5cuXcLixYvBMAzWr1+PXr16sR1Ja9GIT80wDIO9e/fCw8MDw4cPx/nz56n0FMDExASLFi3C8uXL2Y7CioCAAAwaNAg7d+7EokWLqPTUUM+ePREZGYmFCxdi6tSp6Nu3L6Kjo9mOpZVoxKdGsrOzMXXqVDx9+hSBgYE6OzJRluLiYjg4OCA4OBidOnViO45KSKVSLF26FMHBwTh27BhatWrFdiRSBSUlJdi1axe+//579O7dG6tXr6bDKRSIRnxq4uTJk3Bzc4ODgwMiIyOp9JTA0NAQK1as0JmTXPLy8jBkyBBERkYiMjKSSk+D6OvrY9asWUhMTETTpk3Rtm1bLFiwAC9fvmQ7mlag4mNZYWEhZsyYgZkzZ+LQoUNYv349DAwM2I6ltSZMmICUlBStv5HokydP0KVLFzRq1Ajnzp2DhYUF25FIDZiYmOC7775DbGwsCgsL4eTkhPXr16O4uJjtaBqNio9FN2/eRJs2bVBUVITo6Gh4eXmxHUnr8fl8rF69GkuXLoW2zvKHhYWhS5cumDFjBnbs2AF9fX22I5Fasra2xo4dO3Dt2jXcvHkTDg4O+P333yGVStmOpplUfSo2YRixWMysWrWKsbKyYv7++2+24+gcqVTKuLm5Mf/88w/bURTut99+Y6ysrJjz58+zHYUoUXh4ONO1a1emVatWzPHjx+mOLNVEi1tU7NGjR/Dx8YGpqSl+//132NjYsB1JJ508eRKLFi3C/fv3teJwb4lEgnnz5uHChQs4duwYWrRowXYkomQMw+DEiRNYsmQJLCwssGHDBp1ZtFVbNNWpIgzDYOfOnejSpQvGjh2LM2fOUOmxqH///qhfvz4OHDjAdpRae/XqFfr27YukpCTcuHGDSk9HcDgcDBw4ENHR0Rg3bhyGDRuG4cOHIzExke1oao+KTwWysrIwaNAg7Nq1C1euXMGsWbNoHxXLOBwO1q1bh5UrV2r0UVHx8fHo1KkT3NzccOLECZiamrIdiagYj8fD5MmTkZiYiPbt28PT0xO+vr54/vw529HUFhWfkh07dgzu7u5wdXXFjRs30LJlS7YjkX9169YNLVu2xK5du9iOUiNnzpyBl5cXli5dik2bNmnFlC2puTp16mDJkiWIj4+HoaEhWrVqhZUrVyI/P5/taOqH3UuM2is/P5+ZPHky06xZM+batWtsxyEVuHPnDmNtbc0UFBSwHaXKZDIZ89NPPzHW1tbM1atX2Y5D1NTTp08ZHx8fxtramtm6dStTUlLCdiS1QSM+Jbhx4wbc3d3BMAyio6Ph6enJdiRSgTZt2sDLywtbtmxhO0qViEQiTJ48Gfv370dERAS6du3KdiSippo2bYqAgACcPn0ax48fh7OzM/7v//5Pa7fxVAet6qxAToEIwbfTEP88D3lCCUwEPDhZm2BEO1uYG8nfYC4Wi/Hdd99h9+7d2LFjB4YOHari1KQmEhMT4enpicTERLW++8WLFy/wxRdfwMrKCn/++SeMjIzYjkQ0yMWLF7Fo0SLo6elhw4YN8Pb2ZjsSa6j4PhCdmottYUm4nJgNABBJ/rtFiIDHBQPA29ESvt3t4WZn+u6x+Ph4+Pj4wMrKCvv27aPbvWiYKVOmwNLSEuvWrWM7ilzR0dEYPHgwxo4di2+//ZbulE5qRCaTISgoCP7+/u9OgXFxcWE7lspR8b0nMOIZ1pyKh1AiRWV/KxwOIODpwb+/E8Z0aoLt27dj1apVWL16NaZNm0YrNjVQamoq3N3dERsbq3Y/tBw5cgRTp07Fr7/+ilGjRrEdh2gBkUiE3377DWvXrkW/fv3w3XffoXHjxmzHUhkqvn+9Lb2HKBZX/SaQAh4Xpk8vQpZ4BYGBgXBwcFBiQqJs8+fPR0lJCbZu3cp2FABv936uWbMGO3fuxJEjR9C+fXu2IxEt8+bNG2zcuBE7duzApEmTsHTp0ipP99fkcpC6oOLD2+nNUbsjUCwue+5d3u3jKIy5iJLsZ6jbsjssBswr97l6kOLvaV3QtikdAqzpsrOz4eTkhKioKDRr1ozVLEVFRZg8eTIeP36MkJAQOuyAKFVGRga+/fZb/PPPP1i0aBFmz54NgUAg97k1vRykTuhCAYBtYUkQSsof9sozMke9LiNh5PpZhZ8r4+hh17VnSkxHVMXS0hKzZs3CypUrWc2Rnp4OLy8vcLlcXL58mUqPKJ2NjQ127tyJq1evIjw8HI6Ojti/f3+5Q7ADI55h1O4InH+YBZFEVqb0AED478fOxWVh1O4IBEY8U+Gfoup0vvhyCkS4nJgt95peHccuqOPQGVxDkwo/n2GA0IRsvCwQKTElURU/Pz+cPXsWsbGxrLz/zZs30alTJwwfPhyBgYEwNDRkJQfRTU5OTjhy5AgOHjyIXbt2oU2bNjh16hQYhnnvclDlayCAt98Xi8VSrDn1UC3LT+enOn+7/BibLySW+8nlfa+vBECalyN3qhN4O7yf95kDpnk1V1ZMokKbNm3CtWvXcOTIEZW+78GDBzF37lzs2bMHgwcPVul7E/IhhmFw7NgxLFmyBPWbuyPH3Qclcu6CVBh3GbnX/4I0Lxt6devD/POvIbD770bahnw9BE31gKutqerCf4TOj/jin+dVWnpVIZTIEJ9JxwJpC19fX9y6dQuRkZEqeT+ZTIalS5fC398fly5dotIjaoHD4WDw4MGIiYmBcadhEMlZ+Ff89C5eh/0Bi/5fw27+32gw5gfwTMuuihZKpNgelqSq2FXCYzsA2/KEEoW8TtT9WPxRfA9NmjRB48aNYWdnRzcA1VCGhoZYsWIF/P39lX6n9vz8fIwdOxavXr3CzZs3YWlpqdT3I6S6coVSPCupCw63fPG9uXYA9Ty/hEEjJwAAz7j8Ir/3Lwepy2pPnS8+E4Fi/gp4MjEuXbqE5ORkJCcnIzMzExYWFu+KsEmTJuV+b2JS8bVDwq6JEydi48aNuHjxInr16qWU93j69CkGDRoEDw8P/N///R/9oETUUvDtNLkfZ2RSiDKTYGjfCem/fQVGWoI6LTxg2mMSuPyyBccBEHwnTW0uB+l88TlZm8CA91zudCcjkwKlvxgZGEkJwNUDh1v2FHwBj4tRn3XDNK8J7z4mlUqRkZHxrghTUlJw//59HD9+HCkpKUhOTgafz6+wFJs0aQIrKys6oYMlfD4f3333HZYuXYqIiAiFH0pw5coV/O9//8PSpUsxe/ZsOvSAqK2KLgdJC3MBmQRFCdfRwGc9OFw9ZB/+Hm/Cg1C/+7gyz1W3y0E6X3zD29li8wX5N258c/0Q3lz/691/F8aGop7nlzDtNqbM8xgAw9valvmYnp4e7OzsYGdnJ/cgYYZh8OrVq3elWFqQkZGR736fl5cHOzu7CkeNNJ2qXCNHjsT69esREhKi0HNX9+zZg6VLlyIwMBC9e/dW2OsSogwVXQ7i/DuqM243EDwjs7e/7zBEbvG9fR2x8kJWk84Xn4WRAbo7WOL8w6xyS3RNu40pV3If4nCAHo6W1Z675nA4MDc3h7m5Odq2bSv3OcXFxWVKMSUl5d10akpKCjIyMmBhYVHpqJGmU2uOy+VizZo1WLRoEQYNGlTr+91JJBIsWLAAp0+fxtWrV+Ho6KigpIQoT0WXg/QERtCTc02v4tfhKypSrel88QHATG97XH2UU+7klqoQ8PTg622vhFRvF1k4OjpW+A3y/enU0oKMiYnByZMn35Ulj8eTW4qlv2/QoAFNp1aif//+WLt2LQ4cOIBx48r/FFtVr1+/xsiRI8HhcBAREaHWd4Eg5H2VXQ4ycvkU+bdPwPCTdoAeD/lRR1HHvkO55wl4XDg1NFZF3CrR+X18pWpyVqchnwv//i3h49FUecFqoXQ69cNRY+nv359OrWjUSNOpb6/HjR8/HgkJCTX6u0hISMCgQYPQr18//Pjjj+Dx6OdNojlyCkTwXH9J/joIqQSvLuxCYdxlcHh81HXqhvo9JoLDK/vvxIDHRfjinmqzqpOK7z01uTuDupZeVb0/nfphKZZOp5qbm1c4ldq4cWPUq1eP7T+G0vXr1w8DBgzAyPFTqnUw77lz5+Dj44O1a9diypQpLCQnpPamBkThXNxzvF2fWT0cDtDHuQF+81GfQ9ap+D5wP+3tAaxnotNgYKCPEul/fz2lB7D2cLSEr7e9Wp1EoCyl06mVjRpLp1MrGjVqw3Rq0IUILNh3HoaftAOHw/nowbwMw2DLli344YcfEBQUBC8vL/bCE1ILMpkMUxavxiW4Arzqz3io48ktVHxyREVFwWfydMz7NQjxmfnIE4phIuDDqaExhrdV/1tuqBLDMHj9+nWFpZiSkoI3b97A1ta2wlGjra0tDAzU9++0dCaguEQMcCou8NKZgMW9W+D6/h8QERGB48ePo2nTpqoLS4gCiUQiTJgwAWlpafD5did+DkvWistBVHxy+Pv7QyaTqe3duDVN6XSqvFJMTk4uM50qbyq1SZMmrE2n1uTaL0daApusSJzZ6g9jY/W5oE9IdeTm5mLo0KEwNzdHYGAgBAKB1lwOouKTo1WrVti3bx86derEdhSdIJVKkZmZKbcU359OragUlTWdWtF9GkuJX6UjY+8s1HXyhMXABWUeM+RzETS1s1pN7xBSVampqejfvz969uyJn376qcxWnvtpudgeloTQhGxw8HZzeilNuRxExfeBxMREeHt7Iy0tTeOvS2mLyqZTS39f2XRq6dmp1Z1OnRoQJXd/Z6msQ8vBSETg1bMqV3zqeEGfkKqIiYnB559/jjlz5sDPz6/CU4VeFogQfCdNIy8H0brqDxw9ehSDBw+m0lMjHA4HZmZmMDMzQ5s2beQ+p7i4GKmpqWVGiWFhYe9+XzqdWtmo8f3p1Mru0wi8vRULV1AXfHMnSHIzyz2ujgfzEvIxoaGhGDlyJH755Rd8+eWXlT7X3MhAbc7erC4qvg+EhIRgxYoVbMcg1WRoaAgHBwc4ODjIffzD6dSUlBQ8ePAAJ0+efDdq5HK570pQbO8NicAB8u7cJRMVIffqATT4cg0Kos9VmEndDuYlpDKHDh3CnDlzEBQUhB49erAdR6mo+N7z/PlzxMXFaf3/dF2kp6cHW1tb2NrawtPTs9zj70+npqSkYMe9YkiF8kf9uVcCYOTWGzyTym8hpG4H8xIiD8Mw2LRpE7Zs2YKLFy/CxcWF7UhKR8X3nuPHj6Nv3746f1KJLvpwOvVo7i3Ex78o97ySrCcQJkej4cRfqvS66nQwLyEfkkqlmD9/Pi5evIjr16/Dzs6O7UgqQcX3npCQkFqdx0i0R0UH8wpTYiB5k4W07RMBAEyJEGBkyMyZK7cM1elgXkLeJxQK4ePjg5cvX+LatWswNTVlO5LKUPH9Kz8/H1evXsVff/318ScTrVfRwbxG7n1Qt+V/p7Dk3fwHkjdZMOszs9xrqNvBvISUevXqFQYPHoxGjRrhzJkzan2AhDLQ0sV/nTlzBp6ennQbHwLg7X0a5eHyBdAzqv/uF4cvAIenD7065TfYy7tPIyFsS05ORteuXdGpUyccPHhQ50oPoBHfOyEhIRgyZAjbMYiaqOw+je+r8H6NMhlcLfi0lYGolXv37mHAgAFYuHAh5s6dy3Yc1tCID0BJSQlOnz6NwYMHsx2FqJGZ3vYQ8Gp281k+j4Pru1dg/vz5KCoqUnAyQqrv/Pnz6N27N37++WedLj2Aig8AEBYWBicnJ1hbW7MdhagRNztT+Pd3giG/ev9MDPlcrBzYGvfDTuD58+dwc3PDlStXlJSSkI8LCAiAj48PDh8+jOHDh7Mdh3U01Qma5iQVKz1gt6YH8x48eBBHjx7Fl19+iS+++ALr1q2DkZGR8oMTgrd79NavX4/ffvsNoaGhcHZ2ZjuSWtD5szplMhns7OwQGhpa4akfhLx/MK+4pAQy7n8/M1blYN7Xr19j3rx5uHz5Mvbs2YNevXqpLjzRSVKpFLNnz8b169dx+vRp2NjYsB1Jbeh88d28eRMTJkxAXFwc21GIBnhZIELXsQvQ/tPB0DeqV+2DeU+dOoVp06bh888/x4YNG2gVMVGKoqIijB49GgUFBfjnn3/o6+wDOl98S5cuBQCsXbuW5SREE0ilUhgbGyMrK6vG99p78+YN/Pz8cP78eezatQt9+vRRcEqiy3JycjBw4EDY29tj7969dBKVHDq/uIWu75HqePr0KSwtLWt1g9l69ephz5492L17N6ZNm4bJkycjNzdXcSGJznry5Ak8PT3h7e2NP//8k0qvAjpdfAkJCXjz5g3at6d7ppGqiY2NRevWrRXyWr1790ZMTAwMDAzQunVrnDhxQiGvS3TT7du30bVrV8ydOxfr1q2r8D56RMeLj+69R6rrwYMHaNWqlcJez9jYGNu3b0dAQADmzp2LsWPH4tWrVwp7faIbzpw5g759+2L79u3w9fVlO47a0+nv+DTNSaorNjZWocVXqkePHrh//z7MzMzQunVrhISEKPw9iHb6/fffMWHCBBw9epS+n1WRzi5uyczMhLOzM7KysmgenFSZm5sb9u3bh3bt2intPa5evYpJkyahffv2+PXXX2FhYaG09yKai2EYrF69Gn/88QdOnz4NR0dHtiNpDJ0d8R0/fhz9+vWj0iNVJpFIkJiYiJYtWyr1fbp164bo6Gg0atQILi4u+Pvvv5X6fkTzSCQSTJs2DUePHkV4eDiVXjXpbPGFhIRg6NChbMcgGiQpKQk2NjaoU6eO0t+rTp06+PHHH3HkyBGsWLECw4cPR1ZWltLfl6i/wsJCDBkyBCkpKQgLC6OjFmtAJ4svLy8P165dQ9++fdmOQjSIIld0VpWHhwfu3r2LFi1awNXVFQcPHoSOXp0gAF68eIEePXrA0tISx48fr9W2Gl2mk8V3+vRpdOvWjb5oSLUoekVnVQkEAqxbtw4nT57EunXrMGTIEGRmZqo8B2FXUlISunTpgr59+2Lfvn3g8/lsR9JYOll8tJqT1ISyVnRWVfv27REVFQVXV1e4ublh//79NPrTETdv3kS3bt2waNEifPfdd7RHr5Z0blWnSCSCtbU1Hj58SHPjpFpatWqFgwcPws3Nje0ouHv3LiZOnAgbGxvs2rULtrZ0p3dtdeLECUyaNAn79u3DgAED2I6jFXRuxBcWFgZnZ2cqPVItJSUlePz4sdqsnmvTpg1u3bqFzp07o02bNti9ezeN/rTQrl278NVXX+HEiRNUegqkc8VH05ykJhITE9GkSRMIBAK2o7zD5/OxfPlyXLp0CTt37kSfPn2QnJzMdiyiAAzDYMWKFdiwYQOuXr2Kjh07sh1Jq+hU8clkMjrdgNQIGys6q8rFxQURERHo2bMn2rdvjx07dkAmk7Edi9SQWCzGpEmTcObMGYSHh8Pe3p7tSFpHp4rv1q1bqF+/Plq0aMF2FKJh2FrRWVU8Hg9LlizB5cuXsX//fvTq1QtPnjxhOxappvz8fAwcOBDZ2dkIDQ2FlZUV25G0kk4VH01zkppie0VnVTk7O+P69esYMGAAOnbsiC1bttDoT0M8f/4c3t7eaNy4MUJCQlC3bl22I2ktKj5CqkCdpzo/pKenBz8/P4SHh+Pvv/9G9+7d8ejRI7ZjkUokJCSgS5cuGDJkCHbu3Akej8d2JK2mM8UXHx+P/Px8pR4uTLSTUChESkqKxk2ROzg44PLlyxgxYgQ6d+6MTZs2QSqVsh2LfCA8PBzdu3fH8uXLsXz5ctqjpwI6U3whISF07z1SI/Hx8fjkk0808kBzLpeLOXPmIDIyEidOnEDXrl3x8OFDtmORfx05cgRDhgzBH3/8gYkTJ7IdR2foTAvQodSkpjRpmrMizZs3x8WLFzF27Fh069YNP/zwAyQSCduxdNq2bdswa9YsnD59ms4NVjGdKL6MjAwkJiaie/fubEchGkhTFrZ8DJfLha+vL6KionDx4kV4eHggJiaG7Vg6RyaTYcmSJdiyZQuuXbtGl19YoBPFd+zYMfTv358OdSU1ou5bGaqradOmOHfuHKZPn46ePXti9erVEIvFbMfSCSUlJRg3bhyuXLmC69evo1mzZmxH0kk6UXy0mpPUhraM+N7H4XAwZcoU3LlzBzdu3EDHjh1x7949tmNptby8PPTv3x+FhYW4cOECLCws2I6ks7T+kOo3b97Azs4O6enpdBsiUm2FhYWwsLBAfn6+1i4xZxgGf/75JxYuXIjp06dj2bJlGrmQR51lZGSgX79+6Nq1K7Zs2QI9PT22I+k0rR/xnT59Gl5eXlR6pEYePnwIBwcHrS094O3ob/z48bh37x7u3buHdu3aISoqiu1YWiMuLg5dunTBl19+ia1bt1LpqQGtLz6a5iS1oY3TnBWxsbHB0aNH8c033+Dzzz/HN998A6FQyHYsjXb16lX06NEDq1evxpIlS2iPnprQ6uITiUQ4c+YMBg4cyHYUoqG0YStDdXA4HIwePRr379/Ho0eP0LZtW0RERLAdSyMFBwdj2LBhOHDgAMaOHct2HPIerS6+0NBQtG7dGg0aNGA7CtFQ2rais6oaNGiA4OBgfPvttxg6dCgWLFiA4uJitmNpjJ9//hlff/01zp07h08//ZTtOOQDWl18NM1JakuXpjrlGTFiBO7fv4/09HS4ubnh2rVrbEdSazKZDH5+fti1axeuX78Od3d3tiMRObR2VadMJkOjRo1w5coVjTtjkaiH/Px8WFtbIy8vjxYk4O3xWjNnzsSIESOwdu1aunvAB0QiEcaPH4+MjAyEhITAzMyM7UikAlo74ouMjIS5uTmVHqmxuLg4ODk5Uen9a+jQoXjw4AFev34NV1dXhIWFsR1JbeTm5qJPnz6QSqU4d+4clZ6a09rio7M5SW3p6vW9ypiZmeHPP//EL7/8grFjx8LX1xf5+flsx2JVamoqunbtCnd3dwQFBUEgELAdiXyEVhYfwzDvTj0npKZ0bUVndQwYMAAxMTEQiURwcXHB+fPn2Y7EipiYGHTp0gUTJ07E5s2b6e4vGkIr/y/Fx8ejuLgYbdu2ZTsK0WC6vrDlY0xNTbF3717s3LkTU6ZMwVdffYU3b96wHUtlQkND8emnn2Ljxo3w8/OjPXoaRCuLr3Q1J30hktqgqc6q6dOnD2JiYqCnpwcXFxecPn2a7UhK99dff2HkyJEICgrCqFGj2I5DqkkrV3V26tQJa9euRa9evdiOQjRUbm4u7Ozs8ObNG5q+qoaLFy/iq6++gpeXFzZv3oz69euzHUmhGIbBpk2bsGXLFpw8eRIuLi5sRyI1oHX/otPT05GUlAQvLy+2oxANFhsbC2dnZyq9aurVqxfu378PY2NjuLi44NixY2xHUhipVIqvv/4a+/fvR3h4OJWeBtO6f9V07z2iCDTNWXNGRkb49ddfcfDgQfj5+WHMmDF4+fIl27Fqpbi4GCNHjsT9+/dx9epV2Nrash2J1ILWFR+d1kIUgVZ01p6Xlxeio6PRoEEDuLi44PDhw2xHqpFXr17hs88+A5/Px5kzZ2Bqasp2JFJLWlV8ubm5uHHjBvr06cN2FKLhaEWnYtSpUwc//fQTgoOD4e/vj//973948eIF27GqLDk5GZ6enujcuTMOHDgAAwMDtiMRBdCq4jt9+jS6d+8OIyMjtqMQDUdTnYrVpUsX3L17F82aNYOrqysOHToEdV9Xd+/ePXh6emLGjBnYuHEjXe/VIlq1qnPkyJHo3bs3Jk+ezHYUosFycnJgb2+P169f05YYJbh58yYmTpwIBwcH7NixA9bW1mxHKuf8+fMYM2YMduzYgWHDhrEdhyiY1vwIIxQKcfbsWbr3Hqm10mlOKj3l6NixI+7cuYNWrVrBzc0NAQEBajX6CwgIgI+PDw4fPkylp6W0pvguXboEV1dXWFlZsR2FaDia5lQ+AwMDfP/99zh9+jR+/PFHDBw4EOnp6axmYhgG69atw/LlyxEWFoZu3bqxmocoj9YUH63mJIpCKzpVp23btrh16xY6dOiANm3aYN++fayM/qRSKWbOnImgoCCEh4ejZcuWKs9AVEcrik8qleLYsWMYPHgw21GIFqAVnaqlr6+PlStX4sKFC9i2bRv69euHlJQUlb1/UVERvvjiCzx69AhXrlyBjY2Nyt6bsEMrii8yMhJWVlZo3rw521GIhmMYhoqPJa6uroiIiICXlxfatWuHnTt3Kn30l5OTg169eqFevXo4efIkTExMlPp+RD1oRfHRNCdRlKysLDAMgwYNGrAdRSfx+XwsXboUYWFh2LdvHz799FM8ffpUKe/15MkTdOnSBT169MD+/fuhr6+vlPch6ofHdoDaKr33XlBQENtRiBYovb5HKzrZ1apVK1y/fh2bN29Gx44dsXLlSvj6+la6ly6nQITg22mIf56HPKEEJgIenKxNMKKdLcyNym48j4qKwqBBg7B8+XLMmDFD2X8comY0fh9fXFwc+vbti+TkZPpmRWpty5YtiI+Px/bt29mOQv6VkJCASZMmgcfjYe/evbC3ty/zeHRqLraFJeFyYjYAQCSRvXtMwOOCAeDtaAnf7vZwszPF6dOnMW7cOOzZs4fWBegojR/x0b33iCI9ePAAbm5ubMcg73F0dMSVK1fw66+/wsPDA/7+/pgzZw709PQQGPEMa07FQyiRQt6P8MJ/S/BcXBauJOagh+lr/LNhHo4dO4bOnTur+E9C1IXGX+Oj63tEkWgrg3rS09PD119/jYiICISEhKBbt274MSQSa049RLFYfum9j2GAYrEUpzINsHAXlZ6u0+ipzrS0NLi5ueH58+d0GyJSawzDoH79+khKSoKFhQXbcUgFZDIZVvyyFwHp5uDw/7t2l7JpeJnnMZISGLfpD7Pe08t83JCvh6CpHnC1NVVFXKKGNHqq89ixY/j888+p9IhCpKenw8DAgEpPzXG5XLywaAPuiyy8/1N7Y7/gd7+XlQiR9qsP6jh1Lff5QokU28OS8JtPexWkJepIo6c6jxw5QtOcRGFomlMz5BSIcDkxG5VNVRUlXIdenXowsCu/H5NhgNCEbLwsECkvJFFrGjPi+3CpsoArQ3SJJTp07cF2NKIlaOO6Zgi+nfbR5xTEXETd1j0rXPTGARB8Jw3TvOjQC12k9sVX2VLluh4j8emvEWWWKhNSUw8ePECnTp3YjkE+Iv55XpnvAx+SvHkBUeoDmPefU+FzhBIZ4jPzlRGPaAC1nuoMjHiGUbsjcP5hFkQSWbkvdobLg0giw7m4LIzaHYHAiGfsBCVagaY6NUOeUFLp4wUPLsHA1hl808rv85cnFCsyFtEgalt8b/fnVG+p8ppTD6n8SI3IZDLExcXRVKcak0gkiIqKQtqTxEqfV/jgEoxa9/zo65kIaFGcrlLLqc7o1FysORWPYvF/IzxGIsbLc9shfHYPMmEBeKYNUb/7OBg2/29lVrFYhjWn4uFqa0pLlUm1pKSkwMTEBKampmxHIf+SSCS4d+8ewsLCEBYWhmvXrsHW1ha2n00Ez7AeJEz563fCtIeQFryUu5rzfQIeF04NjZUVnag5tRzxbQtLglAiLfMxRiYFz9gC1qN/gN28IJh6+SD76HpIcrPKPK90qTIh1UHTnOwrHdH9+OOPGDBgACwsLDBhwgQ8e/YMEyZMQGJiIh48eIDA1bOgp6cn9zUKH1xEHYcu4BrUqfS9GADD29oq4U9BNIHajfjeLVX+YHqTqy+Aabcx7/67jn1H8Oo1gOh5Enim/52k//5S5Q8PpiWkIrSiU/UqGtF5e3tjwoQJ2LdvH6ysrMp9noWRAbo7WOL8w6xy3yfM+8766PtyOEAPR0v6/qDD1K74qrJUGQCkha8hfpUOfcvG5R6jpcqkumJjY+Hl5cV2DK1W06KTZ6a3Pa4+ykGxWPrxJ39AwNODr7f9x59ItJbaFd/HlioDACOVIOfYjzBy6QW+uV25x2mpMqmuBw8ewNfXl+0YWkWRRfchNztT+Pd3+ncBXOXfL95nyOfCv78TrQHQcWpXfB9bqswwMuSc2ATo8WD22fQKn0dLlUlVyWQyxMfHw9nZme0oGk2ZRSePj0dTAKj07gylOJy3Iz3//k7vPo/oLrU7pPrroLsIuZch9zGGYfDy1C+QvMmC1YhV4PIrnqO3Fafjy2ZiODo6wtHREba2tnTrIiLX48eP0bNnTyQnJ7MdRaNUVnTe3t7w8vJSaNFV5H5aLraHJeHiwyyIxWJA779tCqX34+vhaAlfb3sa6REAajjic7I2gQHvudzpzldnt0H8MhUNRn1faenxOUATU33cv38bf//9NxISEpCfn48WLVq8K8L3fxkZGSnzj0TU3IMHD2hFZxWoekRXVa62pvjNpz2+/WETIrM5cGzXE3lCMUwEfDg1NMbwtuXvwE50m9qN+HIKRPBcf6lc8UnevED6jkmAHh8c7n9Lmc36zoRRq7LndRrwuAhf3LPMF/ubN2+QmJiIhISEMr8ePXoEMzMzuYXYuHHjCpdNE+2xdu1a5ObmYsOGDWxHUSvqMqKrqqFDh2LkyJEYNWoU21GImlO74gOAqQFRcpcqVwWHA/RxblDlW47IZDKkpqYiISEB8fHxZUrx5cuXaN68+bsidHJyevf7evXqVT8cUUtjxoxB7969MX78eLajsErTiu59DMOgYcOGiIyMRJMmTdiOQ9ScWhZfdGouRu2OqNFSZUXeZLKwsFDuKDExMRF169aVO0ps1qwZeDy1m0EmlXBzc8O+ffvQrl07tqOolCYX3YeePn0KT09PpKen07V88lFqWXzA+2d1Vnepckulr9piGAYZGRnlRogJCQnIzMxEs2bN5I4Szc3NlZqLVJ9EIoGJiQlycnJQp07lp31oOm0qug/99ddf+Pvvv/HPP/+wHYVoALUtPqC0/DRrqXJxcTGSkpLKFWJ8fDz4fL7cUWLz5s2hr6/Pam5dlZCQgP79++Px48dsR1E4bS66D82ZMwd2dnZYuHAh21GIBlDrOTkfj6ZwtTXF9rAkhCZkg4O3m9NLqeNSZUNDQ7i4uMDFxaXMxxmGQVZWVpkyvHr1KhISEpCamorGjRvLLUUrKyuaulGiBw8eaM1RZfKKzs7ODt7e3pg4cSJ+//13WFpash1TKW7cuIH//e9/bMcgGkKtR3zve1kgQvCdNMRn5mvdUuWSkhI8fvy43AgxISEBUqlUbiG2aNECAoGA7ega77vvvoNQKMTatWvZjlJtlRVd6YhOW4vufcXFxbCwsEBOTg4MDQ3ZjkM0gMYUn67KyckpN22akJCAp0+fwsbGRm4p2tjY0CixikaOHIlBgwZhzJgxH38yy6jo5Lt27RrmzZuHW7dusR2FaAgqPg0lkUjw9OlTudswioqK4ODgUK4QHRwcULduXbajq5XWrVsjMDAQ7u7ubEcph4quajZu3IjU1FRs2bKF7ShEQ1DxaaHc3Fy5o8SkpCRYWlrKHSXa2dmBy1XL2zMqTUlJCerVq4fXr1+rxbQxFV3NfPHFFxgxYgS+/PJLtqMQDUHFp0OkUilSUlLkbsPIzc2Fvb19uW0YDg4OMDExYTu6UsTGxuKLL75AQkICK+9PRVd7DMPAxsYGN27cQNOmTdmOQzQEFR8BAOTn51e4Wb9evXpyR4lNmzbVuCPdcgpECL6dhvjneUh4koLnqc8wbeQAjGin/EVSVHSK9+zZM3Tu3BkZGRl0XZtUGRUfqZRMJkNaWprcqdMXL17gk08+kbtZv379+mxHLyM6NRfbwpJwOTEbAMqcBVu6Lcbb0RK+3e3hZmeqkPekolO+Q4cOISgoCEeOHGE7CtEgVHykxoqKivDo0SO5pSgQCOSOEj/55BPw+fyPv7gCqeogBCo61Zs7dy4aNWqERYsWsR2FaBAqPqJwDMMgMzNTbiGmp6ejSZMmZUaHpb8sLCwUPl2lzKPvqOjY17FjR2zatAndunVjOwrRIFR8RKWEQmGZzfrvL7ThcDhyR4n29vYwMKj+9beKDjt/fmAJRBkJ725vpWdsjkZTd5Z5jrzDzqno1EvpxvXs7GytP2eVKBYVH1ELDMMgOztb7igxOTkZtra2ckvR2tq6wlFiRbe3en5gCeq27gFjtz4V5uFwgM9aWmFqSw4VnZq6fv065s6di6ioKLajEA2j1md1Et3B4XBgZWUFKyurctNWYrEYT548eTdCvHnzJgICApCQkICSkhK5m/XNbJrgcmJ2je7pCAAMA5y9n4brP69Hzy4dtP6sS01048YNdO7cme0YRAPRiI9otFevXskdJT6v7wLjziMBXvm7Xjw/sATinBQAAN+sEUy9xkLQxLXc8wx4HMz/zBHTvJor/c9Bqm/YsGEYNmwYRo8ezXYUomGo+IhWmvvXHRy9nyn3MVFGAvjmduDo8VH48Apenf8NDSduAb9+w3LPHereCJtHuis5LakuhmHQqFEjXL9+Hc2aNWM7DtEwunVGFdEZ+SXSCh8zsHEE16AOODw+jFx6waBRSxQ/ln+dKE8oVlZEUgupqamQSqV0WgupESo+opVMBNW4fM3hAJA/8WEiUO2eQ1I1pdf36LQWUhNUfEQrOVmbwIBX/stbJixA8ZPbYCQlYGRSFMSGQpT6AIbN2pZ7roDHhVNDY1XEJdVEC1tIbdCqTqKVhrezxeYLieU+zsikyL0SCPGrNIDDBd/cFpZfLAPf3Lb8cwEMb1v+44R9N27cwMaNG9mOQTQUFR/RShZGBujuYFluH59enXpoOGHzRz+fkclQX/gcepJiAMo9vJpUj1AoxIMHD9CuXTu2oxANRVOdRGvN9LaHgFezu0cY6vPQtDABLVu2xJ9//gla/Kw+7ty5AycnJ7qpMqkxKj6itdzsTOHf3wmG/Op9mRvyuVj2eUsc2r4BR48exZYtW+Dl5YWYmBglJSXVQdf3SG1R8RGt5uPRFP79W8KQr4ePLQDkcN6e0fn+AdUdO3ZEZGQkxowZg169emH+/PnIy8tTfnBSISo+UltUfETr+Xg0RdBUD/RxbgADHheCD1Z7CnhcGPC46OPcAEFTPcrdlUFPTw/Tp09HbGws3rx5g5YtW+Kvv/6i6U8WMAyDGzduwMPDg+0oRIPRyS1Ep7wsECH4ThriM/ORJxTDRMCHU0NjDG9b9Tuwh4eHw9fXF2ZmZti6dSucnZ2VnJqUSk1NRbt27ZCVlUV7+EiN0apOolPMjQxqffZmly5dEBUVhR07dqB79+6YNGkSli9fDiMjIwWlJBWhjetEEWiqk5Aa4PF4mD17NmJiYpCZmQlnZ2cEBwfT9KeS0fU9oghUfITUgrW1Nf78808EBgZi1apV6Nu3LxITy2+cJ4pBxUcUgYqPEAXw8vLC3bt30adPH3Tp0gXLli1DUVER27G0ikgkQkxMDNq3b892FKLhqPgIURA+n4/58+fj/v37ePz4MZydnXH06FGa/lSQO3fuwNHRkTauk1qj4iNEwWxsbPDXX39h7969WLJkCQYMGIDHjx+zHUvj0TQnURQqPkKUpFevXoiOjoaXlxc6deqEVatWobi4mO1YGouKjygKFR8hSqSvr4/Fixfj7t27ePDgAVq3bo1Tp06xHUsjRURE0MZ1ohC0gZ0QFTp79ixmz54NZ2dn/Pzzz3QH8SpKS0tDmzZt8OLFC9rDR2qNRnyEqFCfPn3erUxs37491qxZA5FIxHYstUcb14kiUfERomIGBgZYtmwZoqKicPPmTbi4uODcuXNsx1JrdH2PKBIVHyEsadq0KY4ePYqffvoJ06dPx4gRI5CWlsZ2LLVE1/eIIlHxEcKyAQMGIDY2Fq1atYK7uzs2bNiAkpIStmOpDZFIhOjoaHTo0IHtKERLUPERogYMDQ2xatUqREREICwsDO7u7ggNDWU7llq4e/cuHBwc6BBwojBUfISoEXt7e5w8eRJr167FxIkTMXr0aGRkZLAdi1V0fY8oGhUfIWqGw+FgyJAhiI2NRbNmzeDq6orNmzdDLBazHY0VdONZomi0j48QNZeQkIBZs2YhKysL27ZtQ7du3diOpFKNGzfGxYsX0aJFC7ajEC1BxUeIBmAYBsHBwZg/fz569uyJDRs2oEGDBmzHUrr09HS4u7vTxnWiUDTVSYgG4HA4GDFiBOLi4mBlZYXWrVtj69atkEgkbEdTqtJpTio9okhUfIRoEGNjY2zcuBFhYWEIDg5Ghw4dcOPGDbZjKQ0tbCHKQMVHiAZq1aoVQkNDsXDhQgwfPhyTJ09GdnY227EUjjauE2Wg4iNEQ3E4HIwePRpxcXEwMTFBq1atsHPnTkilUrajKURJSQnu3buHjh07sh2FaBkqPkI0XL169bB582acP38eAQEB8PDwwK1bt9iOVWt3795FixYtaOM6UTgqPkK0hJubG65cuYJZs2Zh0KBBmD59Ol69esV2rBqj63tEWaj4CNEiXC4X48ePR1xcHHg8HpydnbFv3z7IZDK2o1UbXd8jykL7+AjRYnfu3IGvry+4XC62b98Od3d3tiNVWZMmTXD+/Hk4ODiwHYVoGRrxEaLF2rZti/DwcEyaNAl9+vTBnDlzkJuby3asj8rIyEBhYSGd1kKUgoqPEC3H5XIxZcoUxMXFQSQSwdnZGQEBAVDnyR7auE6UiYqPEB1hbm6OnTt3IiQkBL/88gu6d++OmJgYtmPJRQdTE2Wi4iNEx3Ts2BGRkZH48ssv0atXL/j5+SEvL4/tWGVERETQik6iNFR8hOggPT09zJgxA7GxsXj9+jWcnZ1x6NAhtZj+pI3rRNmo+AjRYZaWlti3bx+CgoLwww8/4NNPP8XDhw9ZzXTv3j00b94cxsbGrOYg2ouKjxACT09PREVFYfDgwfDy8sLixYtRUFDASha6vkeUjYqPEAIA4PF4mDNnDmJiYpCRkQFnZ2ccPnxY5dOfdH2PKBttYCeEyHXlyhX4+vqiUaNG2Lp1q8r21DVt2hRnz56Fo6OjSt6P6B4a8RFC5PLy8sLdu3fRu3dvdO7cGcuWLUNRUZFS3zMzMxP5+fl0WgtRKio+QkiF+Hw+/Pz8EB0djaSkJLRq1QpHjx5V2vTnjRs30KlTJ9q4TpSKio8Q8lGNGjXCoUOHsGfPHixZsgQDBw7EkydPFP4+dH2PqAIVHyGkynr16oXo6Gh069YNHTt2xLfffguhUKiw16dbERFVoMUthJAaSUlJwfz583Hv3j1s2bIF/fv3r9XrlZSUwMzMDBkZGTAxMVFQSkLKoxEfIaRGGjdujODgYGzduhVz587F0KFDkZycXOPXi46OxieffEKlR5SOx3YAQohm69u3L2JiYrBx40a0bdsWfn5+8PPzg4GBQaWfl1MgQvDtNMQ/z0OeUILM5Mdo4D0GLwtEMDeq/HMJqQ2a6iSEKMzTp08xd+5cJCQkYOvWrfjss8/KPSc6NRfbwpJwOTEbACCS/Hd3eB5HBj09HrwdLeHb3R5udqaqik50CBUfIUThjh8/jrlz56J9+/b46aefYGtrCwAIjHiGNafiIZRIUdl3Hg4HEPD04N/fCT4eTVUTmugMusZHCFG4gQMHIjY2Fk5OTnB3d8fGjRvxx7XHWHPqIYrFlZceADAMUCyWYs2phwiMeKaSzER30IiPEKJUSUlJmLJoNZ41HwLw9Ms8JsnNwstz21GSHg/w+Kjr6In6n04Fh6v37jmGfD0ETfWAq62paoMTrUUjPkKIUtnb28Nh6GxwePxyj708tx16dUxhOzsANhN/hTD1AfLvnCzzHKFEiu1hSaqKS3QAFR8hRKlyCkS4nJgNBuWPIZO8yULdll3B4elDz6g+DJu1gzgnpcxzGAYITcjGywKRqiITLUfFRwhRquDbaRU+ZtJ+EArjrkAmFkKSn4PiJ1EwbNa23PM4AILvVPw6hFQH7eMjhChV/PO8MlsW3iewc0HBvbNI/el/ACND3da9YOhQ/sgyoUSG+Mx8ZUclOoJGfIQQpcoTSuR+nGFkyPq/Fajj2AWN/Q7Ddu5ByIQFyA37vYLXESszJtEhVHyEEKUyEcifWJIV50Oalw3jtgPA4fGhZ2gCI9dPUfw4qoLXKb84hpCaoOIjhCiVk7UJDHjlv9Xo1akHXr0GyL97CoxMCpmwAAUxF8G3albuuQIeF04NjVURl+gA2sdHCFGqnAIRPNdfknudryTrCV5d2AXxi6cAVw+Cxi4w6z0DenVNyzzPgMdF+OKedIYnUQha3EIIUSoLIwN0d7DE+YdZ5U5s0W/wCazH/FDp53M4QA9HSyo9ojA01UkIUbqZ3vYQ8PQ+/kQ5BDw9+HrbKzgR0WVUfIQQpXOzM4V/fycY8qv3LceQz4V/fyc6rowoFE11EkJUovQuC3R3BsI2WtxCCFGp+2m52B6WhNCEbHDwdnN6KQGPCwZvr+n5etvTSI8oBRUfIYQVLwtECL6ThvjMfOQJxTAR8OHU0BjD29rSQhaiVFR8hBBCdAotbiGEEKJTqPgIIYToFCo+QgghOoWKjxBCiE6h4iOEEKJTqPgIIYToFCo+QgghOoWKjxBCiE6h4iOEEKJTqPgIIYToFCo+QgghOoWKjxBCiE6h4iOEEKJTqPgIIYToFCo+QgghOoWKjxBCiE6h4iOEEKJTqPgIIYToFCo+QgghOoWKjxBCiE6h4iOEEKJT/h+wExZaoJ2omgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "edges = [(0,1),(1,2),(0,3),(1,4),(2,5),(3,4),(4,5),(3,6),(4,7),(5,8),(6,7),(7,8),(5,7),(0,4)]\n",
    "example_graph = nx.Graph(edges)\n",
    "nx.draw(example_graph,with_labels=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5f082e4e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.  , 0.33, 0.  , 0.33, 0.33, 0.  , 0.  , 0.  , 0.  ],\n",
       "       [0.33, 0.  , 0.33, 0.  , 0.33, 0.  , 0.  , 0.  , 0.  ],\n",
       "       [0.  , 0.5 , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  ],\n",
       "       [0.33, 0.  , 0.  , 0.  , 0.33, 0.  , 0.33, 0.  , 0.  ],\n",
       "       [0.2 , 0.2 , 0.  , 0.2 , 0.  , 0.2 , 0.  , 0.2 , 0.  ],\n",
       "       [0.  , 0.  , 0.25, 0.  , 0.25, 0.  , 0.  , 0.25, 0.25],\n",
       "       [0.  , 0.  , 0.  , 0.5 , 0.  , 0.  , 0.  , 0.5 , 0.  ],\n",
       "       [0.  , 0.  , 0.  , 0.  , 0.25, 0.25, 0.25, 0.  , 0.25],\n",
       "       [0.  , 0.  , 0.  , 0.  , 0.  , 0.5 , 0.  , 0.5 , 0.  ]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#h = [[0 for _ in range (9)] for _ in range(9)]\n",
    "n = example_graph.number_of_nodes()\n",
    "h = np.zeros((n, n))\n",
    "for x in range(n): \n",
    "    for k in example_graph.neighbors(x):\n",
    "        h[x,k] = round(1/example_graph.degree(x),2) \n",
    "\n",
    "h"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "952c97cd",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "a9ae9f2c0eb7e99b85fd90ae657a9202",
     "grade": false,
     "grade_id": "cell-a8413209ddb70ab0",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "A natural proposal transition matrix is $Q(x \\to y) = \\frac{1}{d_x} \\mathbf{1}_{\\{\\{x,y\\}\\in E\\}}$. In other words, when at $x$ the proposed next state is chosen uniformly among its neighbors.\n",
    "\n",
    "__(a)__ Write a function `sample_proposal` that, given a (networkX) graph and node $x$, samples $y$ according to transition matrix $Q(x \\to y)$. _Hint_: a useful Graph member function is [`neighbors`](https://networkx.org/documentation/stable/reference/classes/generated/networkx.Graph.neighbors.html). **(10 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "449da919",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "db0642dd7c9f53225d13fb9674c89a56",
     "grade": false,
     "grade_id": "cell-df840ef2dee49b9f",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def sample_proposal(graph,x):\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "604ce3f4",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "690a774380e6a4383ccb5da2dee8a737",
     "grade": true,
     "grade_id": "cell-c70fec7da167b7be",
     "locked": true,
     "points": 10,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "from nose.tools import assert_almost_equal\n",
    "assert sample_proposal(nx.Graph([(0,1)]),0)==1\n",
    "assert_almost_equal([sample_proposal(example_graph,3) for _ in range(1000)].count(4),333,delta=50)\n",
    "assert_almost_equal([sample_proposal(example_graph,8) for _ in range(1000)].count(5),500,delta=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7a2e658",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "99c1e08eed2976f4625452d9eed6a5d5",
     "grade": false,
     "grade_id": "cell-4c74e7dae57e50e3",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(b)__ Let us consider the Markov chain corresponding to the transition matrix $Q(x \\to y)$. Produce a histogram of the states visited in the first ~20000 steps. Compare this to the exact stationary distribution found by the function `stationary_distributions` from the lecture applied to the transition matrix $Q$. _Hint_: another useful Graph member function is [`degree`](https://networkx.org/documentation/stable/reference/classes/generated/networkx.Graph.degree.html). **(15 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f52fcbd5",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "44377fd9fa8cd0c3bc78ab03753db1ce",
     "grade": false,
     "grade_id": "cell-ca4f5c3685b72c8a",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def chain_Q_histogram(graph,start,k):\n",
    "    '''Produce a histogram (a Numpy array of length equal to the number of \n",
    "    nodes of graph) of the states visited (excluding initial state) by the \n",
    "    Q Markov chain in the first k steps when started at start.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "def transition_matrix_Q(graph):\n",
    "    '''Construct transition matrix Q from graph as two-dimensional Numpy array.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "# Compare histogram and stationary distribution in a plot\n",
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "812cd2ef",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "c3e08e0577c953ca68fe8159cadc0634",
     "grade": true,
     "grade_id": "cell-bdea40714e10d0e8",
     "locked": true,
     "points": 15,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert_almost_equal(transition_matrix_Q(example_graph)[3,4],1/3,delta=1e-9)\n",
    "assert_almost_equal(transition_matrix_Q(example_graph)[3,7],0.0,delta=1e-9)\n",
    "assert_almost_equal(transition_matrix_Q(example_graph)[2,2],0.0,delta=1e-9)\n",
    "assert_almost_equal(np.sum(transition_matrix_Q(example_graph)[7]),1.0,delta=1e-9)\n",
    "assert chain_Q_histogram(nx.Graph([(0,1)]),0,100)[1] == 50\n",
    "assert len(chain_Q_histogram(example_graph,0,100)) == example_graph.number_of_nodes()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ebe4946",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "f4338397c60be0fc0069e30bec121faa",
     "grade": false,
     "grade_id": "cell-5a9debbf863b9e3b",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(c)__ Determine the appropriate Metropolis-Hastings acceptance probability $A(x \\to y)$ for $x\\neq y$\n",
    "and write a function that, given a graph and $x$, samples the next state with $y$ according to the Metropolis-Hastings transition matrix $P(x \\to y)$. **(10 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37804448",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "639522f9ff68136165c0d9daa174b4ba",
     "grade": false,
     "grade_id": "cell-92e3b11a3af7eb99",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def acceptance_probability(graph,x,y):\n",
    "    '''Compute A(x -> y) for the supplied graph (assuming x!=y).'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "def sample_next_state(graph,x):\n",
    "    '''Return next random state y according to MH transition matrix P(x -> y).'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "058d3728",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "6bf922f8787f02c80681390ce8a8f84c",
     "grade": true,
     "grade_id": "cell-fbfc2607999d0c99",
     "locked": true,
     "points": 10,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert_almost_equal(acceptance_probability(example_graph,3,4),0.6,delta=1e-9)\n",
    "assert_almost_equal(acceptance_probability(example_graph,8,7),0.5,delta=1e-9)\n",
    "assert_almost_equal(acceptance_probability(example_graph,7,8),1.0,delta=1e-9)\n",
    "assert_almost_equal(acceptance_probability(nx.Graph([(0,1)]),0,1),1,delta=1e-9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbf205bf",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "da69b710a3b1986ac9f9997235c1a472",
     "grade": false,
     "grade_id": "cell-97f7d289177acf39",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(d)__ Do the same as in part (b) but now for the Markov chain corresponding to $P$. Verify that the histogram of the Markov chain approaches a flat distribution and corroborate this by calculating the explicit matrix $P$ and applying `stationary_distributions` to it. _Hint_: for determining the explicit matrix $P(x\\to y)$, remember that the formula $P(x\\to y) = Q(x\\to y)A(x\\to y)$ only holds for $x\\neq y$. What is $P(x\\to x)$? **(15 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5de68351",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "bd25c9f6b39133973fd2d1594e210b30",
     "grade": false,
     "grade_id": "cell-3b7fde395331916d",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def chain_P_histogram(graph,start,n):\n",
    "    '''Produce a histogram of the states visited (excluding initial state) \n",
    "    by the P Markov chain in the first n steps when started at start.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "def transition_matrix_P(graph):\n",
    "    '''Construct transition matrix Q from graph as numpy array.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "# plotting\n",
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b5cc504",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "1b756f8dff589048e791db63f89d8624",
     "grade": true,
     "grade_id": "cell-8c4b6c60ee96f037",
     "locked": true,
     "points": 10,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert_almost_equal(transition_matrix_P(example_graph)[3,4],1/5,delta=1e-9)\n",
    "assert_almost_equal(transition_matrix_P(example_graph)[3,7],0.0,delta=1e-9)\n",
    "assert_almost_equal(transition_matrix_P(example_graph)[2,2],0.41666666,delta=1e-5)\n",
    "assert_almost_equal(np.sum(transition_matrix_P(example_graph)[7]),1.0,delta=1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4e9d8aa",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "22ff259ba7ad6066040f3417cc8d4483",
     "grade": true,
     "grade_id": "cell-a63274314d1b2a2d",
     "locked": true,
     "points": 5,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert len(chain_P_histogram(example_graph,0,100)) == example_graph.number_of_nodes()\n",
    "assert_almost_equal(chain_P_histogram(example_graph,0,20000)[8],2222,delta=180)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89f65139",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "9b277d17c7b55bb15e5708e6ca90c4f9",
     "grade": false,
     "grade_id": "cell-da50533d4d46c293",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## MCMC simulation of disk model\n",
    "\n",
    "**(50 points)**\n",
    "\n",
    "Recall that in the disk model with we would like to sample the positions $x = (x_1,y_1,\\ldots,x_N,y_N)\\in [0,L)^{2N}$ of $N$ disks of radius $1$ in the torus $[0,L)^2$ with uniform density $\\pi(x) = \\mathbf{1}_{\\{\\text{all pairwise distance }\\geq 2\\}}(x) / Z$, where $Z$ is the unknown partition function of the model. We will assume $L > 2$ and $N\\geq 1$. For the purposes of this simulation we will store the state $x$ in a `np.array` of dimension $(N,2)$ with values in $[0,L)$. Such a configuration can be conveniently plotted using the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "683d2de3",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "4c7cb2d90103cd790145cae15d46e99e",
     "grade": false,
     "grade_id": "cell-d661f575ab9f80ea",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def plot_disk_configuration(positions,L):\n",
    "    fig,ax = plt.subplots()\n",
    "    ax.set_aspect('equal')\n",
    "    ax.set_ylim(0,L)\n",
    "    ax.set_xlim(0,L)\n",
    "    ax.set_yticklabels([])\n",
    "    ax.set_xticklabels([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_xticks([])\n",
    "    for x,y in positions:\n",
    "        # consider all horizontal and vertical copies that may be visible\n",
    "        for x_shift in [z for z in x + [-L,0,L] if -1<z<L+1]:\n",
    "            for y_shift in [z for z in y + [-L,0,L] if -1<z<L+1]:\n",
    "                ax.add_patch(plt.Circle((x_shift,y_shift),1))\n",
    "    plt.show()\n",
    "    \n",
    "# Example with N=3 and L=5\n",
    "positions = np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4]])\n",
    "plot_disk_configuration(positions,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efa2f1c9",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "d29e40f9d0e8d304d266f016150aeab7",
     "grade": false,
     "grade_id": "cell-527056f298dedf7c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(a)__ Write a function `two_disks_overlap` that tests whether disks at position $\\mathbf{x}_1 \\in [0,L)^{2}$ and position $\\mathbf{x}_2 \\in [0,L)^{2}$ overlap and a function `disk_config_valid` that checks whether a full configuration is valid (non-overlapping and non-touching). _Hint:_ The minimal separation in the $x$-direction can be expressed as a function of `x1[0]-x2[0]` and the minimal separation in the y-direction as a function of `x1[1]-x2[1]`. Then use pythagoras. **(15 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c010eb",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "bbf8ae650d9d4a66a3e4ac4722a3f60e",
     "grade": false,
     "grade_id": "cell-1b2a61bf719003e0",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def two_disks_overlap(x1,x2,L):\n",
    "    '''Return True if the disks centered at x1 and x2 (represented as 2-element arrays) overlap in [0,L)^2.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    \n",
    "def disk_config_valid(x,L):\n",
    "    '''Return True if the configuration x (as two-dimensional array) is non-overlapping in [0,L)^2.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0f2024a",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "05fc821fed0690aabc13fd66da800694",
     "grade": true,
     "grade_id": "cell-9f544deda0526691",
     "locked": true,
     "points": 10,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert two_disks_overlap(np.array([1,1]),np.array([1,1]),5)\n",
    "assert two_disks_overlap(np.array([0.6,0.6]),np.array([4.1,0.5]),5)\n",
    "assert two_disks_overlap(np.array([0.3,0.3]),np.array([4.6,4.6]),5)\n",
    "assert not two_disks_overlap(np.array([1,1]),np.array([3.1,1]),7)\n",
    "assert not two_disks_overlap(np.array([1,1]),np.array([1,3.1]),7)\n",
    "assert not two_disks_overlap(np.array([1,1]),np.array([1.01+np.sqrt(2),1.01+np.sqrt(2)]),6)\n",
    "assert two_disks_overlap(np.array([1,1]),np.array([0.99+np.sqrt(2),0.99+np.sqrt(2)]),6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68e57e57",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "8705cb3a08fc636963fe61806348d093",
     "grade": true,
     "grade_id": "cell-699454de327d56d5",
     "locked": true,
     "points": 5,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert disk_config_valid(np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4]]),5)\n",
    "assert not disk_config_valid(np.array([[0.1,0.5],[2.1,1.5],[3.2,3.4],[4.1,2.3]]),5)\n",
    "assert disk_config_valid(np.array([[1,1],[3.1,1],[1,3.1]]),6)\n",
    "assert not disk_config_valid(np.array([[1,1],[3.1,1],[1,3.1],[2.5,2.5]]),6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe7c5736",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "db023e65fc2613d624f29977a82bf1f1",
     "grade": false,
     "grade_id": "cell-b0e00f11b0f85525",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(b)__ Assuming $N \\leq \\lceil \\frac12 L -1 \\rceil^2$ where $\\lceil r\\rceil$ is the smallest integer larger or equal to $r$, write a function `generate_initial_positions` that produces an arbitrary non-overlapping (and non-touching) initial condition given $N$ and $L$. The layout need not be random, any deterministic layout is ok (e.g. grid). **(10 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f86a04ed",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "b17579d74c8ff8572d81f18997b49c90",
     "grade": false,
     "grade_id": "cell-53c1bd894d6fe27d",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def generate_initial_positions(N,L):\n",
    "    '''Return array of positions of N disks in non-overlapping positions.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    \n",
    "plot_disk_configuration(generate_initial_positions(33,14.5),14.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d3dc3db",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "fa73c21133aa2b913b931b266e7782f0",
     "grade": true,
     "grade_id": "cell-e1c80d8b59301b93",
     "locked": true,
     "points": 10,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "for l in [4,9.2,14.5]:\n",
    "    for n in range(1,int(np.ceil(l/2-1)**2)+1):\n",
    "        assert disk_config_valid(generate_initial_positions(n,l),l), \"Failed for n = {}, l = {}\".format(n,l)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb5dbe07",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "b4bb37c92263c95c89055362336f894a",
     "grade": false,
     "grade_id": "cell-1bd76964ec9620ce",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(c)__ Write a function `remains_valid_after_move` that determines whether in a non-overlapping configuration $x$ moving the $i$th disk to `next_position` results in a valid non-overlapping configuration. **(10 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56252046",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "04ff82e1b0e3386f2cb15412695c1864",
     "grade": false,
     "grade_id": "cell-d54b4fa9b2f8eb92",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def remains_valid_after_move(x,i,next_position,L):\n",
    "    '''Returns True if replacing x[i] by next_position would yield a valid configuration,\n",
    "    otherwise False.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "167bdc3e",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "37fe75b415d4e31369d6b6bb1041609b",
     "grade": true,
     "grade_id": "cell-e4902309b9869f3d",
     "locked": true,
     "points": 10,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "assert remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],0,[4.5,0.5],5)\n",
    "assert not remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],1,[4.5,0.5],5)\n",
    "assert not remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],2,[3.2,2.5],5)\n",
    "assert remains_valid_after_move([[0.1,0.5],[2.1,1.5],[3.2,3.4]],2,[3.2,3.8],5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8cd0c44",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "8e0cfe7c5bce12cc34e939c5c151ed65",
     "grade": false,
     "grade_id": "cell-582ca3e158908e6c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "__(d)__ Implement the Metropolis-Hastings transition by selecting a uniformly chosen disk and displacing it by $ (\\delta\\,\\mathcal{N}_1,\\delta\\,\\mathcal{N}_2)$ where $\\delta>0$ is a parameter and $\\mathcal{N}_i$ are independent normal random variables (make sure to keep positions within $[0,L)^2$ by taking the new position modulo $L$). Test run your simulation for $L=11.3$ and $N=20$ and $\\delta = 0.3$ and about $10000$ Markov chain steps and plot the final state. **(15 pts)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccaa8f9f",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "24f09b7137c0184ddb872ef558f756a9",
     "grade": true,
     "grade_id": "cell-e180b99ca699c610",
     "locked": false,
     "points": 15,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def MH_disk_move(x,L,delta):\n",
    "    '''Perform random MH move on configuration x, thus changing the array x (if accepted). \n",
    "    Return True if move was accepted, False otherwise.'''\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "\n",
    "# Test run and plot resulting configuration\n",
    "# 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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}