941 lines
54 KiB
Plaintext
941 lines
54 KiB
Plaintext
{
|
|
"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": 1,
|
|
"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": 2,
|
|
"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": 3,
|
|
"id": "071e5617",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABIgElEQVR4nO3dd1gU59oG8HsbLIqAVBHsDRsaG0SioigqRY1IIEcTezcxOcYkxmPKMeZoYosxdo2JJSBYIoINpNhQQcWKiooBlS5SZJct8/3hxyaEIsjuzszu87surkt2h9lnFfee9523CBiGYUAIIYQYCSHbBRBCCCH6RMFHCCHEqFDwEUIIMSoUfIQQQowKBR8hhBCjQsFHCCHEqFDwEUIIMSoUfIQQQowKBR8hhBCjQsFHCCHEqFDwEUIIMSoUfIQQQowKBR8hhBCjQsFHCCHEqFDwEUIIMSoUfIQQQowKBR8hhBCjQsFHCCHEqFDwEUIIMSoUfIQQQowKBR8hhBCjQsFHCCHEqIjZLoAQwp68EjnCkzORmlWEIpkSFlIxXJpZILC3M2zMTdkujxCdEDAMw7BdBCFEv1IyCvFzXBri7+YCAORKteY5qVgIBoBnJzvMGdQePVpYsVMkITpCwccDdFVOtGl3YjqWRaVCplShtv/9AgEgFYuw2McFE9xb660+QnSNgo/D6KqcaNvL0LuNMoX61Qf/PzOJEIt9OlP4EYNBwcdRdFVOtC0loxDBWxNRplBVejwvYiVk6SlQK2QQNW4KC/cANOkxvNIxZhIRQme4w9XZSo8VE6IbFHwcRFflRBdm7ErCydvZVS6kynMfQdK0OQRiCRT5Gcjauwj2gV/DtFl7zTECATC8iwM2Teij56oJ0T4a1ckxKRmFWBaVWmPold6KR+HZ36EqyoWocVPY+H4EaYtuKFOosSwqFa7OVnRVTqrIK5Ej/m5utb0HJnat/vadAAIIoHz2tFLwMQwQeycX+SVyuq9MeI+Cj2N+jkuDTKmq9rmyh1fwLG4n7EZ/BpPmHaEqKaj0vEypwoa4NLoqJ1WEJ2fW+nz+8Q0ovR4DRimHiUM7mLWr+jskABB+ORMzB7bTUZWE6AcFH4fUdlUOAM/P7IGlx7swdXIBAIib2FZ6nq7K2cMwDNRqNVQqFZRKpebrn99X95g+jkmSdIXcxLnG+m2Gz4H1sJmQP06F7M/rEIgkVY6RKdVIfVqsy79GQvSCgo9DarsqZ9QqyJ+mway9Gx5vmg5GVY5GHdxhNXgKhJK/Qk6fV+VqtVqnH9xcCY26HiMQCCAWizVfIpGo0vfVPaaLY0xMTGBmZlbp+ZvPrAFZ7f+eAqEI0hZdUXozFsVXomDRZ1SVY4pkCh39NhGiPxR8HJKaVVRpysLfqUoLAbUSL+6chcOEFRAIRcjd/y2enwtF00Hva46TKdXYsf8YzmxN1PmHPcMwr/xA1teHvVQqrdN5dFWfSCSCUMjdFQAfhF7B/atP6nawWg3ls6fVPmUhrdoSJIRvKPg4pEimrPE5wf+36pr09ofY3Prln/uOqRJ8AGBp5wjvPt46+XD/+/dc/qAnlbk0s4CpOKvKhZWqtBCyRykwa98PArEJZOlXUXo7Hrb+C6ucQyIEOjVroq+SCdEZCj4OsZDW/M8hkppD9I97ejXp2qEt3gvqqaWqiCEY19sZa6LvVn1CIEDxlaPIP74BYNQQW9qjqdd0NOroXuVQhUKBn/49Huqp72PChAmwsLDQQ+WEaB8FH4fUdFVewbz7UBQnH4FZ296ASIzipD/QqH3fSsdIxUK4ONJVOanM1twUgzraVZnHJ2pkiWbjl7/y5wUCwNvVGUEjvsPGjRuxePFiBAUFYfbs2ejRo4cOKydE+6ivikPG9a551B0AWHoEw8SxAx5vmYknW2fBxKEdLPsHVTqGATCuV+3nIcZprmd7SMWi1/pZqViEuZ7tMWTIEISFheHmzZtwcnKCn58f+vfvj127dkEme8XoGUI4glZu4ZiaVteoC1pdg7yKtlcFUiqViIyMxMaNG5GcnIyJEydi1qxZaN++fdUTEcIR1OLjmIZelc/xpA8cUrMJ7q2x2KczzCQiCAS1HysQvFyjs7al8MRiMUaPHo1jx44hMTERQqEQ/fv3h7e3Nw4ePAilsuYBW4SwhVp8HERrdRJdu5ZZiA1xaYi9kwsBXk6DqVCx88fgTnaY49m+3kvgyWQyhIeHY+PGjXj06BGmT5+OadOmwcnJSavvgZDXRcHHUbQ7A9GH/BI5wi9nIvVpMc4nX4WlmQRjvdwxrpd29nq8du0aNm3ahJCQEAwePBizZ8/GkCFDaCoMYRUFH4fp8qqckH9atWoVMjIysHbtWq2fu7i4GHv27MHGjRtRVlaGWbNmYeLEibCxsdH6axHyKhR8PPD3q/ITcafRvVM7DO7VSWtX5YQAQGhoKMLDwxEWFqaz12AYBufPn8fGjRsRERGB0aNHY/bs2XBzc4PgVTcdCdES6m/gARtzU8wc2A5rgnqi3eOT+FdbJWYObEehR7TKyckJjx8/1ulrCAQCzfSHtLQ0dOvWDePHj0evXr2wZcsWlJSU6PT1CQEo+HhHKpXSfCmiE05OTsjMrH37Im2ytbXFwoULce/ePSxfvhxRUVFo2bIl5s2bh5s3b+qtDmJ8KPh4hoKP6Erz5s2RlZUFtbruo4m1QSgUYvjw4Th06BBSUlJgY2MDb29vDBw4EL///jvkcrle6yGGj4KPZyj4iK6YmprCysoKOTk5rNXQokULfPPNN0hPT8f8+fOxfft2tGzZEosWLcLDhw9Zq4sYFgo+nqHgI7qk7+7OmkgkEgQEBCA6OhoJCQmQy+Xo27cvfHx8EBERAZVKxXaJhMco+HiGgo/okrOzs84HuNRXp06dsHr1amRkZOCdd97Bt99+i7Zt22LZsmXIyspiuzzCQxR8PEPBR3RJHyM7X5eZmRkmTZqECxcu4ODBg0hPT0fnzp0RFBSEuLg40MwsUlcUfDxDwUd0iStdna/Sq1cvbN26Fenp6RgwYADmzp2Lrl27Yt26dSgsLGS7PMJxFHw8Y2pqSsFHdIaLXZ21sbS0xLx583Djxg1s2rQJ58+fR5s2bTBt2jQkJyezXR7hKAo+nqEWH9ElLnd11kYgEGimP6SmpqJdu3YICAhAv379sGPHDrx48YLtEgmHUPDxjFQqpXlNRGf40tVZGwcHByxatAj379/HV199hQMHDqBly5b46KOPkJqaynZ5hAMo+HiGWnxEl/jW1VkbkUgEX19fHDlyBElJSWjcuDE8PT01u8grFAq2SyQsoeDjGQo+oksWFhZgGAZFRUVsl6JVrVu3xrJly/Dnn39i5syZ+Pnnn9GqVSssWbIEGRkZbJdH9IyCj2co+IguCQQC3t7nqwsTExPN9Ifo6Gg8f/4cPXv21Owir+/l2gg7KPh4hoKP6JqzszPv7/PVRZcuXbBu3Tr8+eef8Pf3x+LFi9GhQwd8//33yM3NZbs8okMUfDxDwUd0zZBbfNVp3Lgxpk2bhqSkJPz++++4ffs2OnTogPHjx+PMmTM0Md4AUfDxDAUf0TVjC74KAoEA/fr1wy+//IIHDx6gT58+mDp1Knr06IENGzYY3H1PY0bBxzMUfETXDGFKQ0NZW1vj448/RmpqKtauXYvY2Fi0bt0as2bNQkpKCtvlkQai4OMZCj6ia4Y0paGhBAKBZvrDzZs34eTkBD8/P80u8vR/kZ8o+HiGgo/omrF2db6Ko6MjlixZgocPH+Kzzz7Dnj170LJlSyxcuBBpaWlsl0fqgYKPZyj4iK5RV2ftxGKxZvrD+fPnIRAI0L9/f3h7e+PgwYNQKpVsl0heQcDQkCVeKS0thb29PUpLS9kuhRgolUoFMzMzlJSUwMTEhO1yeEEmk2H//v3YsGEDHj16hOnTp2P69Olo3rw526WRalCLj2cqdmeg6xWiKyKRCA4ODnj69CnbpfCGVCrF+PHjcfbsWURFRSE7OxvdunXT7CJPE+O5hYKPZ8RiMYRCIXWnEJ2i7s7X5+rqqmn5DRs2DAsWLICLiwtWr16NgoICtssjoODjJbrPR3SNRnY2XJMmTTBr1ixcvXoVO3fuxJUrV9CuXTvNLvLUa8MeCj4eouAjukYjO7WnYvDLrl27cO/ePXTt2hXjx49H7969sXXrVrpfzwIKPh6i4CO6Rl2dumFra4uFCxfi7t27WL58OaKiotCyZUvMmzcPN2/eZLs8o0HBx0MUfETXqKtTt4RCoWb6Q0pKCmxsbODt7a3ZRZ42m9YtCj4eouAjukZdnfrj7OyMb775Bunp6Zg/fz62b9+Oli1bYtGiRUhPT2e7PIMkZrsAUn8UfETXqKtT/yQSCQICAhAQEIC7d+9i06ZN6NOnD9zc3DB79myMHDkSIpHotc6dVyJHeHImUrOKUCRTwkIqhkszCwT2doaNuamW3wn30QR2HhowYAD+97//4a233mK7FGKgysrK0LRpU5SVlUEgELBdjtEqKytDaGgoNm7ciOzsbMyYMQNTp06Fg4NDnX4+JaMQP8elIf7uy/0F5cq/5hNKxUIwADw72WHOoPbo0cJKB++Am6irk4cqJrEToitmZmZo3Lgx8vLy2C7FqJmZmWmmPxw4cADp6elwcXFBcHAw4uLiap0SsTsxHcFbE3HydjbkSnWl0AMA2f8/duJWNoK3JmJ3YrqO3w13UPDxEHV1En2g+3zc0qtXL2zZsgXp6el46623MHfuXHTt2hXr1q1DYWFhpWN3J6ZjWdRtlClUeFWfHsMAZQoVlkXdNprwo+DjIQo+og/Ozs50n4+DLC0tMW/ePNy4cQObNm3C+fPn0aZNG0ybNg3JyclIySjEsqhUlCkqt/BUZcXI2f8t/lwVgMwNk1F6M67S82UKNZZFpeJaZqH+3gxLKPh4iIKP6AO1+LhNIBBopj+kpqaiXbt2CAgIQODX21CmqLqkYcGJjRCIJHD+YDds/T9B/okNKM99VOkYmVKFDXGGv8USBR8PUfARfaDg4w8HBwcsWrQIF1JuQe3QGUDlAUnqchle3DkHq4ETIDQxg7RFVzRq74bSm7GVjmMYIPZOLvJLDHseIQUfD1HwEX2grk7+OXj1KUSiqh/ryoLHEAiFkFg7aR6T2LeB4h8tPuBlZIZfNux/dwo+HqLgI/pALT7+Sc0qqjJ6EwDUijIITBtVekxo2gjq8rIqx8qUaqQ+LdZZjVxAwcdDFHxEHyj4+KdIVv12ZUKJGRh55ZBj5C8gNDGr4TwKrdfGJRR8PETBR/SBujr5x0Ja/WJcYmsnMGoVFAV/XciU5zyExK5VDeeR6KQ+rqDg4yEKPqIPTZs2RXl5OW2bwyMd7c0hFlSduCc0kaJRpzdReHoP1OUyyDJv4UXaBTTuOrjKsVKxEC6OTfRRLmso+HiIgo/og0AgQPPmzam7kweKi4vx448/4vtZb0OpUlV7jLX3HDDKcmT+NB55h3+AjfccmFTT4mMAjOvlrOOK2UWLVPMQBR/Rl4ruzo4dO7JdCqnG48eP8dNPP2Hbtm0YPHgwQn7dhh33xDh5O7vKii0isyawD/hPrecTCIDBnewMfuFqavHxEAUf0Rca4MJN165dw8SJE9G9e3e8ePECFy9eRFhYGNzd3THXsz2k4tfbxUEqFmGOZ3stV8s9FHw8RMFH9IWCjzsYhsHx48fh7e2NESNGwMXFBWlpaVi3bh3atm2rOa5HCyss9nGBmaR+H+9mEiEW+7jA1dlKy5VzD3V18hAFH9EXZ2dn3L17l+0yjJpcLsfvv/+O1atXg2EYLFiwAO+++y5MTWvujpzg3hoAsCwqFTJl7QtVCwQvW3qLfVw0P2foKPh4iIKP6IuTkxNiY2NffSDRumfPnmHz5s346aef0KVLF/zwww/w9vau8/6IE9xbw9XZChvi0hB7JxcCvJycXkEsUEMkEmNwJzvM8WxvFC29ChR8PETBR/SFujr17+HDh1i7di127doFPz8/REVFoUePHq91LldnK2ya0Af5JXKEX85E6tNiFMkUePrnA5TnPETo/z42+IEs1aHg4yEKPqIvzs7OFHx6cuHCBaxatQoxMTGYNm0arl27Bmdn7UwrsDE3xcyB7TTfJycLMXHiUtiYf66V8/MNDW7hIalUCrncsFdPJ9zQrFkz5ObmQqmsfiks0jBqtRp//PEHBgwYgHfeeQf9+/dHeno6VqxYobXQq46rqyvS09NRVFSks9fgMmrx8RC1+Ii+iMVi2NraIisrS6cfxMbmxYsX+O2337B69WpYWlrik08+QUBAAMRi/XwkSyQS9OzZE5cuXYKXl5deXpNLqMXHQxR8RJ/oPp/25OTk4KuvvkLr1q0RFRWFbdu24eLFiwgKCtJb6FVwd3fHhQsX9PqaXEHBx0MUfESfaLHqhktNTcWMGTPQqVMnZGVlISEhAYcPH8bAgQPrPEpT29zc3JCYmMjKa7ONgo+HKPiIPlGL7/UwDIP4+Hj4+/tj4MCBaN68Oe7cuYPNmzfDxcWF7fI0LT6mtkl+Boru8fGQqakpBR/RGwq++lEqlQgPD8eqVavw/PlzLFiwAPv27YOZWfV737HF2dkZIpEI6enpaNOmDdvl6BUFHw+JxWKo1WoolUq93xcgxsfZ2Rk3btxguwzOKy4uxrZt2/Djjz+iZcuWWLJkCfz8/CAUcrNjTSAQaFp9xhZ83PwXIbUSCAQ0pYHoDbX4apeZmYlPP/0UrVu3RmJiIvbt24eEhASMGjWKs6FXwVjv83H7X4XUiO7zEX2h4Kve1atX8d5778HV1RXl5eVISkpCaGgo+vXrx3ZpdWasIzsp+HiKgo/oi5OTEzIzM41yEMQ/MQyDY8eOYejQofD19UW3bt1w//59rF27lpfdhb1798a1a9eMrveIbhDxFAUf0Rdzc3OYmpri2bNnsLa2ZrscVsjlcuzduxerV6+GQCDAJ598guDgYJiYmLBdWoOYm5ujQ4cOSElJ4VVLtaGoxcdTFHxEn4y1u7OgoADfffcd2rRpg5CQEKxevRopKSl4//33eR96Fdzc3Iyuu5OCj6co+Ig+Gdsk9gcPHuCDDz5A+/btcffuXRw7dgzHjx/HsGHDWJtwrivu7u5GN8CFgo+nKPiIPhlLiy8xMRGBgYHo168fzM3NcePGDezcuROurq5sl6Yzxtjio3t8PEXBR/TJkINPpVLh8OHDWLVqFR4/foyPP/4Yv/zyC8zNzdkuTS9cXFyQl5eH3Nxc2NnZsV2OXlDw8RQFH9EnZ2dnJCcns12GVr148QK//vorVq9ejaZNm2LhwoV4++23jW5RCKFQiL59++LixYvw9fVluxy9oK5OnqLgI/pUMaXBEGRnZ+PLL79E69atcfz4cezYsQMXLlxAYGCg0YVeBWO7z0fBx1MUfESfDKGr8/bt25g+fTpcXFyQm5uLM2fO4NChQxgwYIDBDVipL2O7z2eclzcGgIKP6JOzszMvg69ih4SVK1fi0qVLmDt3Lu7evWs097Lqys3NDRcvXoRareb8MmvaQMHHU7RWJ9EnW1tblJSUoKysjHO7DFRHoVAgPDwcK1euRGlpKRYsWICwsDBe1M4GOzs72NjY4M6dO+jcuTPb5egcBR9PUYuP6JNAIICjoyOePHmCdu3asV1OjYqKirBt2zasXbsWbdu2xTfffAMfHx+jaMU0VMV9PmMIPvpt4CkKPqJvXO7uzMjIwMKFC9GmTRtcunQJBw4cQFxcHKe3BeIaY7rPR78RPEXBR/SNiyM7r1y5ggkTJqBHjx5QqVS4fPkyfv/9d/Tp04ft0njHmEZ2UvDxFAUf0TeujOxkGAZHjx6Fl5cX/P394erqigcPHmD16tVo1aoV2+XxVo8ePXDv3j2UlpayXYrO0T0+nqLgI/rm7OyMR48esfb6crkce/bswapVqyCRSLBgwQIEBQUZzGLRbDM1NUX37t2RlJSEQYMGsV2OTlGLj6co+Ii+sdXVmZ+fj2XLlqF169YICwvDunXrcOXKFbz33nsUelpmLBvTUvDxFAUf0Td9d3Xev38f8+bNQ4cOHXD//n2cPHlS08Vp7BPOdcXNzc0o7vNR8PEUBR/RN32N6jx//jwCAgLg7u4OS0tL3Lx5Ezt27EC3bt10/trGrmKAC8MwbJeiU3SPj6dMTU0p+IheOTo6IisrCyqVCiKRSKvnVqlU+OOPP7By5UpkZWXh448/xq+//mo0OyRwRevWraFUKpGZmYkWLVqwXY7OUPDxFLX4iL6ZmJigadOmyMnJgaOjo1bOWVpaip07d2LNmjWwtbXFJ598grffflvrwUrqRiAQaO7zUfARzqHgI/qWVyKHtUcQPj14E5LGmbCQiuHSzAKBvZ1hY25ar3NlZWVh/fr12Lx5MwYMGIBff/0V/fv3p3t3HFBxn2/cuHFsl6IzFHw8RcFH9CUloxA/x6Uh/m4uyjt44XSGHEAOAEAqzsKa6Lvw7GSHOYPao0cLq1rPdevWLaxevRoHDhzAu+++i3PnzqFDhw66fxOkztzd3fH111+zXYZOUfDxFAUf0YfdielYFpUKmVIFhgEgrPyRIVOqAQAnbmUj4W4eFvu4YIJ760rHMAyD2NhYrFq1CsnJyZg3bx7u3r0LW1tbPb0LUh99+/bFlStXoFAoIJFI2C5HJyj4eIqCj+jay9C7jTKF+pXHMgxQplBhWdRtAMAE99ZQKBTYt28fVq1aBZlMhgULFmD//v2QSqW6Lp00gIWFBVq3bo3r16+jV69ebJejExR8PEXBR3QpJaMQy6JSqw09RcFjPNk+D41dPGDr/0ml58oUanwbeRvX4iIRsmEF2rdvj6VLl2LkyJG0WDSPVNznM9Tgo99EnqLgI7r0c1waZEpVtc8VnNgEU8ea78uVlSsR/VSIQ4cO4dSpU/D19aXQ4xlDX8GFfht5ioKP6EpeiRzxd3NR3Rzm0lvxEEobQ9qqR40/LxAKUWrRGq06dtVhlUSXDH2LIgo+nqLgI7oSnlz9epxq+QsUnt6DpkOmvvIcAgDhl7m1hRGpu65du+LJkyd49uwZ26XoBAUfT0kkEiiVSqjVrx54QEh9pGYVQa6s+ntVmLAL5j28Ibawe+U5ZEo1Up8W66I8ogcikQi9e/fGxYsX2S5FJ2hwC08JBAJIpVLI5XKYmZmxXQ7hEbVajYKCAuTk5CA7Oxs5OTmV/pxo0hNoUnlfu/LsB5A9SoHj5B/r/DpFMoWWKyf6VLFu5/Dhw9kuReso+HisoruTgo/IZLIqAVZdqOXk5CAvLw8WFhawt7eHvb09HBwcNH/u3bs38osdkZz/j/P/eR3K59nI3DAZAMCUywBGjad582sMQwupYc4BMxZubm7YsmUL22XoBAUfj9F9PsPFMAyePXtWa4D9/c8ymazaIHNycsIbb7xR6XFbW9ta97Fj4u/jRvTdSt2d5j2Ho3HngZrviy4egPJ5NqyHz632HGKBGq0s6eOFz9zc3DB16lQwDGNwS8nRbyaPUfDxi1wuR25ubp2CLDc3F40bN64SZg4ODujRo0eVgLO0tNTah9O43s5YE3230mNCiRSQ/DXxXCCRQiA2gaiRZbXnUKnU+HriCJxy64WgoCCMGjUKTZo00Up9RD8cHR1hbm6OtLQ0g1tWjoKPxyj42MUwDJ4/f16nIMvOzsaLFy9gZ2dXJbQcHBzQvXv3So/b2dnB1LR+Cz9ri625KQZ1tMPJ29nVTmkAAKsB42v8eYEAGN7dCd/fu4U//vgDe/fuxZw5c+Dt7Y2goCD4+vpS9zxPVNznM7TgEzCGvuOgAcorkSM8OROrd4TAtbcbnB1sXnuVfFKZQqHQhNWrgiw3NxdSqbTaIKvuz02bNuVNl1FKRiGCtp6HrA7Llf2TmUSE0BnucHW20jyWn5+PgwcPIiQkBElJSfDz80NwcDC8vb1r7XYl7Fq9ejUePHiA9evXs12KVlHw8cjfV8kHUOkejFQsBAPUeZV8Y8EwDIqKiuoUZDk5OSguLoatrW2dgsze3t5g151Uq9Xwmvk1Hln3glpY90EqZhIhFvt0rrJQ9d9lZWUhPDwcoaGhuHXrFsaMGYPg4GAMHjwYYjF1QnHJuXPn8OGHHyIpKYntUrSKgo8nqqySXwOBAJCKRdWukm8oFAoF8vLy6jzwQyKR1CnIHBwc0LRpU1peC8Dnn3+O06dPY+rynfjh5H2d/d5lZGRg3759CAkJwZ9//omAgAAEBwfjrbfeon8HDigrK4OtrS3y8vIMqnuago8H6rNKfoW6XHlzBcMwKCkpqXOr7Pnz57CxsalTkNnZ2aFRo0Zsv0Ve2bhxI9auXYuzZ8/C1tYW1zILsSEuDbF3ciHAX1sRAX/1NAzuZIc5nu0rdW/W1/379xEaGoqQkBAUFBQgMDAQwcHB6NevH2+6iA1R3759sXbtWnh4eLBditZQ8HFcSkYhgrcmokzx14LBf66qvDMyoyxHkzd8YO09q9Lj1d1r0RelUon8/Pw6t8oEAkGlkYu1hZq1tTVEIpHe35MxiIiIwIwZM3DmzBm0a9eu0nP5JXKEX85E6tNiFMkUsJBK4OLYBON6af/e8q1btzQhWF5ejuDgYAQHB8PV1ZVCUM/mzZuHtm3b4t///jfbpWgNBR/HzdiVVOvoOnW5DJk/TYB94NeQtuxW6TmBABjexQGbJvTRSi2lpaV1DrJnz56hadOmdWqV2dvbo3Hjxlqpkby+S5cuwcfHB5GRkejXrx/b5QB42Rtw9epVTQhKpVJNCLq4uLBdnlHYvXs3Dh8+jH379rFditZQ8HFYXokcHitOVbtuYoWS6zF4fmYvms/aVu2VsKlYiHOfDan2ilylUqGgoKDOw/HVanWdW2U2NjY0UIFHHjx4gLfeegubNm3CqFGj2C6nWgzD4MKFCwgJCcG+fftgb2+PoKAgBAUFoW3btmyXZ7Du3buHoUOH4tGjR2yXojUUfBy2Kf4+1vxjBY1/ytr7BaQtutY4r0osUKO3SRaaP79VJcgKCgpgaWlZ5+H45ubm1M1kgPLz8+Hh4YEPPvgAc+dWvxIL16hUKpw5cwYhISHYv38/2rRpg+DgYAQGBsLZ2Znt8gwKwzCwtbXFjRs34OjoyHY5WkHBx2EfhV7BoatPanxe+TwHjzdNQ/OZWyCxalbjca2ZbIyyL6wSZLa2tpBIaD1FYyaTyTB06FD0798f33//PdvlvBalUolTp04hJCQEhw4dQrdu3RAcHIxx48bB3t6e7fIMgo+PD2bMmIExY8awXYpWUPBx2JRfL+FUak6NzxeeDYEs/SqajV9e63m8XOyxfWJfbZdHeE6tViM4OBhCoRB79+41iOkDcrkcJ06cQEhICCIjI9G3b18EBwfj7bffhrW1Ndvl8dY333yDsrIyLF9e+2cNX/D/N92AWUhrv0dWeuMUzLsNqcN5qFVHqvr000+RlZWFnTt3GkToAYCpqSn8/f2xZ88ePHnyBDNmzMDRo0fRpk0b+Pn5Yffu3Sgupn0C68vd3d2gdmQ3jN92A+XSzAKm4ur/iWSZt6EqyUcjl7dqPYdULISLIy0OTCr76aefEBkZiUOHDhns6jONGjVCYGAgwsPDkZGRgeDgYISEhMDZ2Rnjxo1DWFgYXrx4wXaZvNCvXz8kJSVBpVK9+mAeoODjsHG9a75JX3ojBo069ofQtPbJ2QyAcb3oZj/5y6FDh/C///0PUVFRRtP9Z2FhgQkTJuDIkSN4+PAhRo4ciS1btqB58+YYP348IiIiIJfL2S6Ts5o2bQonJyfcvHmT7VK0goKPwypWya9uIKXNiHmw9V9Q688LBC9X1KCFq0mFxMRETJ8+HYcPH0abNm3YLocV1tbWmDp1Kk6ePIk7d+7Aw8MDP/zwAxwdHTFlyhScOHECSqWS7TI5x83NDYmJiWyXoRUUfBw317M9pOLXW6VEKhZhjmd7LVdE+CotLQ1vv/02du7ciT59tLOoAd85ODhgzpw5SEhIwLVr19C9e3csWbIEzZs3x+zZsxEfH28w3XsNZUj3+Sj4OK5HCyssHNYOUJXX6+cYhRwTXZuwslwZ4Z68vDyMHDkSX3/9NXx9fdkuh5OcnZ3x8ccf48KFC0hMTETLli0xf/58tGzZEh999BESExNhzIPgDanFR9MZOI5hGEyZMgWPJC3wpNmbkCvVdVolf6RjGX7/7xxER0eja9eu+iuYcE5ZWRm8vLzg6emJ7777ju1yeOf27duaJdPkcrlmtZiePXsa1YIOCoUCTZs2xZMnT2BhYcF2OQ1CwcdxW7ZswY8//ogLFy7gQaGyXqvk7927FwsXLsSpU6fQqVMn1t4DYY9KpUJgYCDMzMywa9cug5m2wAaGYXDt2jWEhIQgNDQUEokEwcHBCAoKQpcuXdguTy8GDBiAr7/+Gl5eXmyX0iAUfBxWsWjwmTNnKgVXfVbJ37lzJ5YsWYK4uLgqq+0Tw/fRRx8hJSUFx44dg6kpDXLSFoZhcOnSJU0I2tjYaELQkP+fLVy4EFZWVli8eDHbpTQIBR9H5eXloU+fPli1ahUCAgIadK4tW7bgu+++Q3x8PFq1aqWlCgnXrV27Flu3bsXZs2dhZWXFdjkGS61W4+zZswgJCUF4eDhatWqFoKAgvPPOO2jRogXb5WlVeHg4fvvtNxw+fJjtUhqEgo+DVCoVfHx84Orqih9++EEr51y/fj3WrFmD+Ph4WsTXCOzfvx/z58/HuXPn0LJlS7bLMRpKpRKxsbEIDQ3FwYMH0aVLFwQFBSEwMBAODg5sl9dgmZmZ6NWrF7Kzs3l9f5OCj4O+/PJLJCQkIDo6Wqtb+6xatQqbN29GfHy8wayyTqo6d+4cxowZg+PHj+ONN95guxyjVV5ejhMnTiA0NBQRERHo06cPgoODMXbsWF4vHODk5IQzZ87weh4o3enmmMjISOzYsQMhISFa389uwYIFmDRpEry8vJCTU/Pi14S/7t69i7Fjx+K3336j0GOZiYkJ/Pz8sGvXLjx9+hSzZ8/G8ePH0aZNG/j6+mLXrl0oKipiu8x6c3d35/20Bgo+Dnnw4AGmTJmC0NBQNGtW8zZDDfHFF18gMDAQQ4cORX5+vk5eg7AjJycHI0eOxLJlyzBixAi2yyF/Y2ZmhoCAAISFhSEzMxPjx49HWFgYWrRogbFjx2Lfvn28WTfUzc2N9xPZqauTI8rKytC/f39MmjQJ8+fP1+lrMQyDRYsW4eTJk4iJiaGBDwagtLQUQ4YMwfDhw/Hf//6X7XJIHT179gyHDh1CSEgILly4gJEjRyI4OBgjRozg7CjchIQEfPrpp7xu9VHwcQDDMJg6dSpevHiB33//XS83jRmGwb///W+cO3cOJ0+e5P2EVGOmUqkwduxYWFlZYefOnbwedGDMcnJysH//foSGhuLatWsYPXo0goKC4OXlxakNo0tLS2Fvb4+CggLOhvOrUPBxwNatW7F27VpcuHAB5ubmentdhmEwb948zTwvfb420Q6GYfDBBx8gNTUVUVFRMDExYbskogWPHz9GWFgYQkNDkZaWhoCAAAQHB2PAgAEQiV5v7V5teuONN7Bp0ya4ubmxXcproeBjWVJSEkaOHInTp0/DxcVF76+vVqsxc+ZMpKWlITIyEo0a1b7NEeGWlStX4rfffsPp06dhaWnJdjlEBx4+fIh9+/YhJCQE2dnZCAwMRHBwMNzd3Vlr3c+aNQtdunTBhx9+yMrrNxQNbmFRfn4+xo0bh40bN7ISegAgFAqxefNmtGjRAmPGjIFMJmOlDlJ/oaGh+PHHHxEVFUWhZ8DatGmDzz77DFeuXEFsbCxsbGwwZcoUtGnTBp9++ikuX76s98Wz+T6yk1p8LFGpVPD19UW3bt2wcuVKtsuBSqXChAkTUFRUhAMHDvC2795YnD59GgEBAYiOjoarqyvb5RA9YxgG169f1yyeLRKJEBQUhODgYL0sSn/79m34+fnh/v37On8tXaDgY8lXX32FuLg4xMTEaH2+3utSKBQIDg6GSqVCWFgYp26ok7/cvn0bnp6e2L17N4YNG8Z2OYRlDMMgKSkJoaGhCA0NhZWVlWbd0PbtdbMfp1qthrW1Ne7duwc7OzudvIYuUfCxICoqCjNmzEBSUpLO5uu9rvLycowbNw5SqRR79+7lTCiTl7KysvDmm2/iq6++wqRJk9guh3CMWq3GuXPnNOuGOjs7Izg4GO+8847Wl64b4jMa3UdNh6CpE4pkSlhIxXBpZoHA3lUXy+caCj49e/jwIdzd3bF//3689dZbbJdTLblcjtGjR8PW1ha//vorJ0aREaCkpASenp4YNWoUvvzyS7bLIRynVCoRHx+PkJAQHDhwAC4uLggODkZgYGCDLrhTMgrxc1waom89hQAMVPjr86FiezTPTnaYM6g9erSwavgb0QEKPj0qKyuDh4cHJk6cqPNJ6g1VVlYGPz8/tGrVCtu2baN93FimVCoxZswYODg4YNu2bTRXj9RLeXk5oqOjERISgoiICPTq1QtBQUEICAiAjY1Nnc+zOzEdy6JSIVOq6rQh9mIfF0xwb93wN6BlFHx6NHXqVJSWluptknpDlZaWYuTIkejatSs2bNjAi5oNEcMwmD17NtLT0xEREUH3XkmDlJWV4ejRowgNDcWxY8fg4eGBoKAgjBkzptbRwS9D7zbKFOoaj/knM4kQi306cy78KPj0ZNu2bVizZo3eJ6k3VHFxMby9vdGvXz+sXbuWwo8Fy5cvR2hoKBISEtCkSRO2yyEGpKSkBBEREQgNDUVsbCyGDBmC4OBg+Pn5oXHjxprjUjIKEbw1EWUKleaxouQIlF6PQXluOhp3HgRbv4+rfQ0ziQihM9zh6myl67dTZxR8epCcnIyRI0ciISGBtfl6DVFYWIihQ4diyJAhWLFiBYWfHu3duxeLFi3C+fPn0bx5c7bLIQassLBQs27o+fPnK60b+mHYDZy8nV2pe/PFnXOAQICyh5fBKMprDD6BABjexQGbJvTR0zt5NQo+HcvPz0efPn3www8/YNy4cWyX89oKCgowZMgQ+Pv7Y+nSpWyXYxRiY2MRFBSE2NhYvczNIqRCbm4uDhw4gJCQEKSkPoDlxJ/ACKof5PYsYRdURXk1Bh8AmIqFOPfZEM6M9qQRCzqkUqkwfvx4BAQE8Dr0AMDa2honT57EwYMH8e2337JdjsG7efMmgoODERoaSqFH9M7Ozg4zZ85EbGwsPt9yCMIG9vIIAIRfztROcVpAk7R06L///S/KysqwfPlytkvRCjs7O0RHR8PT0xOmpqZYuHAh2yUZpCdPnsDX1xerV6/G4MGD2S6HGLknLwBVA9tIMqUaqU+LtVRRw1Hw6UhUVBS2b9+OpKQkg5oE3qxZM8TExGDQoEEwMTHh/LQMvikuLoavry9mzJiB8ePHs10OISiSKbV0HoVWzqMNhvOJzCEPHz7E5MmTsX//fs6tzKINTk5OOHXqlCb8Zs+ezXZJBkGhUCAwMBD9+vXDokWL2C6HEACAhVQ7MWEh5c40HAo+LSsrK0NAQAAWLVrE2ZVZtKFly5aIiYmBp6cnTExMMHXqVLZL4rWKuXoikQg///wzjZwlnOHSzAKm4izIlZXn7zFqFVDxxajBKMsBoQgCYdVBMFKxEC6O3JmKQ8GnZfPmzUPHjh2Noguwbdu2iImJweDBg2FiYoL33nuP7ZJ4a9myZbhy5Qri4+MNqmuc8N+43s5YE323yuPPz4bg+dnfNd+X3oyFpce7sBpQtYueATCul7Muy6wX+h+mRdu2bcP58+dx8eJFo7li79ChA06ePAkvLy+YmJggKCiI7ZJ457fffsP27dtx/vx5Xi1uQIyDrbkpBnW0qzKPz2rA+GpD7p8EAmBwJzvOTGUAKPi0Jjk5GYsWLcLp06eN7sOrc+fOOH78OIYNGwYTExO8/fbbbJfEGzExMVi4cCHi4uIM8n4wMQxzPdvj9L28Siu31JVULMIcT91sj/S6aB6fFnBhJ3W2de/eHUePHsWsWbNw5MgRtsvhhevXr+Pdd99FWFgYOnfuzHY5hNSoRwsrLPZxgVRcv8h4uVanC6eWKwMo+BqsYufysWPH8n6SekO98cYbiIiIwJQpU3DixAm2y+G0zMxM+Pr6Yt26dRg4cCDb5RDySuPdWsH+8WmIGBVedSdHIHi5RicXF6gGKPgabOnSpSgtLTWYSeoN1a9fPxw6dAgTJkxAbGws2+VwUlFREXx9fTFv3jwEBwezXQ4hdRIaGoq8xEMImeGO4V0cYCoWVmkBSsVCmIqFGN7FAaEz3DkZegCt1dkgR48exbRp05CUlARHR0e2y+GU+Ph4BAYG4sCBAwY9raO+FAoFfHx80LFjR6xfv95oBkERfsvKykKPHj1w5MgR9O3bFwCQXyJH+OVMpD4tRpFMAQupBC6OTTCuF+3AbrAqdlIPDw/HgAED2C6Hk6Kjo/Gvf/0LERERcHNzY7sc1jEMg8mTJ6OgoAAHDhygaQuEFxiGwdtvv42uXbti2bJlbJejFdTV+RpkMhnGjRuHzz//nEKvFkOHDsXOnTsxatQoJCcns10O67755hvcunULv//+O4Ue4Y09e/bg/v37+PLLL9kuRWuoxfcapk2bhuLiYoSEhFBXVR388ccfmDlzJo4fP44ePXqwXQ4rduzYgWXLluHcuXNwcHBguxxC6uTJkyfo2bMnjh49it69e7NdjtbQZWc9bd++HefOnTOqSeoNNXr0aJSXl2PEiBGIjo42um12jh8/ji+++ALx8fEUeoQ3GIbBzJkzMWvWLIMKPYCCr16Sk5Px+eefG+Uk9YYKDAyEQqGAt7c3Tp06hU6dOrFdkl5cvXoV7733Hg4ePGg075kYht9++w1//vkn9u/fz3YpWkfBV0cFBQVGP0m9of71r3+hvLwcQ4cORVxcHNq1a8d2STr1559/ws/PDxs2bICHhwfb5RBSZ48fP8bChQtx4sQJmJiYsF2O1lHw1YFaraZJ6loyadIklJeXw8vLC/Hx8WjVqhXbJelEYWEhfHx88O9//5t+ZwivMAyD6dOnY+7cuejZsyfb5egEBV8dLF26FCUlJTRJXUtmzJiB8vJyDBkyBPHx8XB25s6q7dpQXl6OsWPHwsvLCx9//DHb5RBSL7/88guePn2KL774gu1SdIZGdb7CsWPHMHXqVJqkrgOrVq3C5s2bER8fbzB/twzD4P3330dpaSnCwsIgElXdm4wQrsrIyECvXr0QExMDV1dXtsvRGWrx1SI9PR0TJ05EeHi4wXwwc8mCBQsgl8vh5eWFuLg42Nvbs11Sgy1ZsgRpaWmIiYmh0CO8wjAMpk2bhvnz5xt06AEUfDWiSer68cUXX0Aul2Po0KGIjY2FjY0N2yW9ti1btiA0NBTnzp1Do0aN2C6HkHrZtm0b8vPz8dlnn7Fdis5RV2cNpk+fjqKiIpqkrgcMw2DRokU4efIkYmJiYGVlxXZJ9RYVFYWpU6fi9OnTaN+eW3uPEfIqjx49Qp8+fRAbG4tu3bqxXY7OUYuvGjt27MCZM2dokrqeCAQC/O9//4NcLsfw4cNx8uRJWFhYsF1WnSUnJ2PSpEk4fPgwhR7hHYZhMHXqVCxYsMAoQg+gFl8Vly9fxvDhw5GQkECbg+oZwzCYN28eUlJScOzYMV4sEpCeng4PDw+sX7+edp4nvLRp0yb88ssvOHv2rNGsIUvB9zcFBQXo3bs3vv/+ewQGBrJdjlFSq9WYOXMm0tLSEBkZyel7Zc+ePYOHhwdmzZqFDz/8kO1yCKm3hw8fom/fvkhISECXLl3YLkdvKPj+n1qthp+fH1xcXLB69Wq2yzFqarUakyZNQlZWFg4fPgypVMp2SVXI5XJ4e3ujT58+WLVqFdvlEFJvarUaXl5e8PHxwcKFC9kuR68o+P7fN998g+joaJw6dQoSiYTtcoyeUqnEhAkTUFxcjAMHDsDUlDsbW6rVaowfPx5KpRKhoaEQCml3L8I/69evx549e3DmzBmjm3pDwQeapM5VCoUCQUFBUKvVCAsL48wFyeeff44zZ84gOjqak61RQl7l/v37cHNzw9mzZ41y8XSjv1StmKQeEhJCoccxEokEISEhUKlUmhYW2zZu3IiDBw/ijz/+oNAjvKRWqzF58mR88cUXRhl6gJEHX8Uk9c8++4wmqXOUiYkJwsLCUFRUhEmTJkGlUrFWy+HDh7F06VIcPXqU1xPtiXH76aefoFarMX/+fLZLYY1Rd3VOnz4dhYWF2LdvH83X47iysjL4+fmhVatW2LZtm97vq126dAk+Pj6IjIxEv3799PrahGjLvXv38Oabb+L8+fPo0KED2+WwxmhbfBWT1Hfs2EGhxwNmZmY4fPgw0tLSMHfuXOjzeu3BgwcYPXo0tm/fTqFHeEulUmHy5MlYsmSJUYceYKTBd/nyZXz22Wc4cOAAmjRpwnY5pI4aN26MyMhIXL16FR999JFewi8/Px8jR47Ef/7zH4waNUrnr0eIrvz4448QCoX44IMP2C6FdUbX1VlQUIA+ffpgxYoVNEmdpwoLCzF06FAMGTIEK1as0FmLvaysDMOGDYOHhwdWrFihk9cgRB/u3LkDDw8PXLhwAe3atWO7HNYZVfCp1Wr4+/ujU6dONEmd5woKCjBkyBD4+/tj6dKlWj+/Wq1GUFAQxGIx9uzZQ3P1CG+pVCq89dZbmDBhAubOnct2OZxgHAuz/b9vv/0WRUVFdPVuAKytrXHy5EkMHjwYpqam+M9//qPV8y9cuBA5OTk4ceIEhR7htdWrV0MqlWL27Nlsl8IZRhN8x44dw+bNm5GUlMSZidCkYezs7BAdHQ1PT0+YmppqbdmldevWISoqCmfPnuXUijGE1NetW7ewYsUKXLp0iS7g/sYogi89PR2TJk1CWFgYTVI3MM2aNUNMTAwGDRoEExOTBs9NOnjwIFasWIGzZ8/C2tpaS1USon9KpRKTJk3Ct99+izZt2rBdDqcYfPBVTFL/9NNPaZK6gXJycsKpU6cwaNAgmJqaYtasWa91nsTERMyYMQPHjh1D69attVskIXq2cuVKWFpaYubMmWyXwjkGP7hlxowZePbsGU1SNwIPHjyAp6cnvv76a0yZMqXSc3klcoQnZyI1qwhFMiUspGK4NLNAYG9n2JibIi0tDQMGDMD27dvh4+PD0jsgRDtu3LiBwYMHIykpCa1atWK7HM4x6BbfL7/8gtOnT9NO6kaibdu2iImJweDBgyGRSPDee+8hJaMQP8elIf5uLgBArlRrjpeKs7Am+i7ebG2J05v/g2+++YZCj/CeQqHApEmT8N1331Ho1cBgg+/KlSv49NNPkZCQQJPUjUiHDh1w4sQJDB06FMlFjXA8uzFkShWq69eQ/X8Ixt3Lh8jrIzRyddVztYRo34oVK2BjY4Np06axXQpnGWRX57Nnz9C7d28sX74c77zzDtvlEBas2H8OG85nQSD5a1Rm1p7PIX9yBwLhy73HRE1s4DRjs+Z5M4kQi306Y4J7a32XS4hWXLt2DV5eXrh8+TJatGjBdjmcZXAtPrVajQkTJmD06NEUekYqJaMQO1OKKoVeBWvvWWjSY3i1P1emUGNZVCpcna3g6myl4yoJ0S6FQoGJEydixYoVFHqvYHATO5YtW4aioiJ8//33bJdCWPJzXBpkytfbvkimVGFDXJqWKyJE97777js4Ojpi8uTJbJfCeQbV4jt+/Dg2bdpEk9SNWF6JHPF3c6u9pwcAhXG/ojDuV0isnWA18D1IW1W+r8cwQOydXOSXyGFjTpPXCT9cvXoVP//8M65cuUID+eqAN8H3quHojx49wsSJE7Fv3z6apG7EwpMza3yu6eDJkNi0gEAkQentBOTsXwrHyesgaVr590UAIPxyJmYOpMV8CfeVl5dj4sSJ+OGHH+Dk5MR2ObzA+eCry3D0Ae1tkPTbMixcuBADBw5kq1SiYyqVCsXFxSgqKqryVfF4RJ415Grban/etHknzZ/Nu3uh9FY8yu4nQdLHv9JxMqUaqU+LdfpeCNGWb7/9Fi1btsT777/Pdim8weng252YjmVRqa8cjh59OxvCPlNg37+Hniskr8IwDORyeY1B9aqvvx9XVlYGc3NzWFhY1PhVbmZX9+IEAgDV94kWyRTa+QsgRIeSk5OxadMmXL16lbo464Gzwfcy9G6jTKF+9cECIdQCIb47ehsCAWg4uhao1WqUlJTUOaBqCzMAsLS0rBJSTZo0qfS9g4NDjYHWpEkTNG7c+JUL7X4UegWZV59UfT+yEsif3IG0ZXdAKELp7QTIM27A2mt6tecxNzG4cV/EwMjlckycOBFr1qxB8+bN2S6HVzgZfCkZhVgWlVol9JSF2cg/sQHlj1MBsQSNO3mg6dAZmnlZNBz9ZX9/XcLqVceUlpaiUaNGNYZUxZednR3atWtX63H63OHApZkFTMVZlbrEAYBRq1CYsBuKgkxAIITExhl2Y/8DiY1z1ZOoFAjZtAqP95fDz88PPj4+sLe319M7IKRu/vvf/6J9+/b417/+xXYpvMPJCewzdiXh5O3sKt2b2fu+gqiRFWxGzIVaVors0P/AvMdwWPQZpTlGIACGd3HApgl99Fz162MYBi9evHjtkPr7l1KprLUr8FWtroovc3NziEQitv9q6i2vRA6PFaeqBF99mIqFODzVFefjTiIiIgLR0dHo2rUr/Pz84O/vj27dulG3EmHVpUuX4Ofnh5SUFDRr1oztcniHcy2+2oajK59nw6K3HwRiE4jMTWDWpjcUeX9WOkafw9GVSmW1oVTf+1fFxcUwMTF5ZUhZWlqiRYsWtR4jlUqN+kPZ1twUgzraVXvhVBcCATC4kx06tXZCp0mTMGnSJMjlcsTHxyMiIgL+/i8HwlSEYMVegIToi0wmw8SJE/Hjjz9S6L0mzrX4NsXfx5rou9VesRdfiYI88zasR8yFWlaCnNAvYTVgAhp16l/pOKlYiI+Hdax2ODrDMJDJZA0aZFHxJZPJamwx1fR4dcc1adKE5h1qUUpGIYK3JqJMUf9J7GYSEUJnuNfYVc4wDG7evIkjR44gIiICN27cgJeXF/z8/ODr6wsHB4cGVk9I7T7//HOkpaUhLCzMqC9yG4JzwfdR6BUcqmZwAgAo8jKQF7ES5TkPAUaNxt28YOP7UbX/+LbFD+Dw8Fi1ISYSiV4rpP55bOPGjekXj6PqNTjq/73OWp25ubk4evQoIiIicPLkSXTq1An+/v7w9/eHq6sr/X4QrUpMTMSYMWNw7do1uu/cAJwLvim/XsKp1JwqjzOMGo83TkWTniNg0W8s1Ioy5Ef+CImNE5oOnlLleJcmCsxzFVcbXNQ1ZRxeNR2mgkAASMUiLPZxadCI4PLyciQkJCAiIgIRERFQKpWaLtHBgwdDKpW+9rkJKSsrwxtvvIGlS5ciMDCQ7XJ4jXPBV1OLT/XiOTLXjUeLj0IhlDYGALy4ex6FCbvQfNqGKse/3dMJa4J66rpcwnHXMguxIS4NsXdyIcBfcz+Bl13iDF7e05vj2V6rI4EZhsHt27c1XaIpKSkYMmSIpkuUVhci9fXJJ58gIyMDoaGhbJfCe5wb3FLTcHRRI0uILR1QfCUKFm5jwZSXoeR6DCT2baqcQyoWwsWR9uAjgKuzFTZN6IP8EjnCL2ci9WkximQKWEglcHFsgnG9nHUyCEogEKBLly7o0qULPv30U+Tl5eHYsWOIiIjAwoUL0aFDB01rsGfPntQlSmp17tw57NmzB9evX2e7FIPAuRZfbcPRy7MfoCB6CxQ5DwGhCNKW3WHtPRuixlaVjjMVC3HusyG0yDDhJIVCgdOnT2u6RGUymSYEhwwZAjMzM7ZLJBzy4sUL9OzZE8uXL8fYsWPZLscgcC74gJrn8dUFH+fxEePFMAzu3LmDiIgIHDlyBFeuXIGnpyf8/f3h6+tLK3IQfPzxx8jOzsbevXvZLsVgcDL4dDkcnRAuKygo0HSJHj9+HG3bttW0Bnv16kVdokbm9OnTCAoKwvXr12FjY8N2OQaDk8EH6G84OiFcpVAocPbsWU2XaGlpKXx9feHv7w8vLy80atSI7RKJDpWWlqJHjx5YtWoVRo8ezXY5BoWzwQfofzg6IVx29+5dTZdocnIyBg4cqOkSdXauZs1Rwmsffvghnj17hl27drFdisHhdPAB7A1HJ4TLnj17huPHjyMiIgLHjh1Dq1atNF2ivXv3fuUuFoTb4uLiMH78eFy/fh3W1tZsl2NwOB98FfQ9HJ0QvlAqlTh37pymNVhYWKjpEh06dCgaN27MdomkHkpKSuDq6op169bBz8+P7XIMEm+CjxBSN2lpaZoQvHjxIgYMGAB/f3/4+fmhRYsWbJdHXmHu3Ll48eIFfvnlF7ZLMVgUfIQYsMLCQhw/fhxHjhzB0aNH4ezsrAnBvn37Upcox8TExGDSpEm4fv06rKys2C7HYFHwEWIklEolEhMTNaNE8/PzNV2iw4YNg7m5OdslGrWioiK4urpi48aNGDlyJNvlGDQKPkKM1P379zVriV64cAEeHh6a1mCrVq3YLs/ozJw5EyqVCtu2bWO7FINHwUcIwfPnz3HixAkcOXIEUVFRcHR01IRgv379IBKJ2C7RoJ04cQLTp0/HtWvXYGlpyXY5Bo+CjxBSiUqlQmJioqY1mJOTAx8fH/j7+8Pb2xtNmtAC8Nr0/PlzdO/eHdu2bYO3tzfb5RgFCj5CSK0ePnyoCcHz58/jzTff1Gy227p1a7bL471p06ZBJBJh8+bNbJdiNCj4CCF1VlxcrOkSjYyMhL29vaZL1N3dnbpE6+no0aOYPXs2rl+/Ti1pPaLgI4S8FpVKhYsXL2pag0+fPoWPjw/8/PwwfPhwWFhYsF0ipxUWFqJ79+7YuXMnvLy82C7HqFDwEUK04tGjR5oQPHv2LNzc3DRdom3btmW7PM6ZPHkyzMzMsGHDBrZLMToUfIQQrSspKcHJkyc1XaLW1taaLtE333wTYrGY7RJZFRkZiQ8++ADXrl2j+ZMsoOAjhOiUWq3GpUuXNK3BzMxMjBw5En5+fhgxYoTRDd9/9uwZunfvjt27d8PT05PtcowSBR8hRK8yMjI0IXj69Gn069dPs7NE+/bt2S5P595//31YWVlh3bp1bJditCj4CCGsKS0tRXR0NCIiIhAZGQlLS0tNl6iHh4fBdYn+8ccfWLBgAVJSUmjXDBZR8BFCOEGtViM5OVnTGnz06BFGjBih6RJt2rQp2yU2SH5+Prp3747Q0FAMGDCA7XKMGgUfIYSTMjMzERkZiYiICCQkJKB3796aLtGOHTuyXV69jR8/Hvb29lizZg3bpRg9Cj5CCOe9ePECMTExmn0Gzc3NNVMlPDw8IJFI2C6xVgcOHMDnn3+Oq1evolGjRmyXY/Qo+AghvKJWq3HlyhVNl+iDBw8wfPhw+Pv7Y8SIEbC2tma7xEpyc3Ph6uqK8PBweHh4sF0OAQUfIYTnnjx5oukSjYuLwxtvvKHpEu3UqRMEAgGr9QUFBaFFixZYuXIlq3WQv1DwEUIMRllZGU6dOqXpEjUzM9OE4IABA3TSJZpXIkd4ciZSs4pQJFPCQiqGSzMLBPZ2xqmjh7FkyRJcuXIFZmZmWn9t8noo+AghBolhGFy9elUTgvfu3YO3tzf8/f0xcuRI2NjYNOj8KRmF+DkuDfF3cwEAcqVa85xULISaYVB2PwnL3x+C8SPfatBrEe2i4COEGIWnT58iKioKEREROHXqFHr06KFpDXbu3LleXaK7E9OxLCoVMqUKtX6CMgzMTMRY7OOCCe6tG/weiHZQ8BFCjI5MJkNsbKymNSiRSDQhOHDgQJiYmNT4sy9D7zbKFOoaj/knM4kQi306U/hxBAUfIcSoMQyDa9euaUIwNTUVw4YN03SJ2tnZaY5NyShE8NZElClUf/28UoH8ExsgS78KtawEYitHNB30Psza9an0OmYSEUJnuMPV2Upfb43UgIKPEEL+Jjs7G5GRkThy5AhiYmLQrVs3zTJq667IcPJ2dqXuTXW5DEUX9sO8+1CILO1Qdj8JeYd/QPMp6yG2ctAcJxAAw7s4YNOEPtW8KtEnCj5CCKmBTCZDfHw8IiIiEHEiFoIxywDRq0eGPtk+D5Ye76KxS+V5e6ZiIc59NgQ25qa6KpnUgZDtAgghhKukUimGDx+O9evXY9HWP2Bi8urQU5U+g6LgMUzsWlZ5TgAg/HKmDiol9UHBRwghdZCaVYxyVe3HMCol8g6vhHl3L0hsWlR5XqZUI/VpsY4qJHVFwUcIIXVQJFPW+jzDqJF3ZBUgEsN62KxazqPQdmmknij4CCGkDiykNe8NyDAM8qPWQVVaCLu3v4BAVPOxFlJuL6htDCj4CCGkDlyaWcBUXP1HZsHxn6HIz4D9uC8hlNQ8cEUqFsLFsYmuSiR1RKM6CSGkDvJK5PBYcarS0mQAoHyeg8cbpwAiCQRCkeZx6xFzYd51cKVjaVQnN9TcHieEEKJha26KQR3tqszjE1vao9XnR1758wIBMLiTHYUeB1BXJyGE1NFcz/aQikWvPrAaUrEIczzba7ki8joo+AghpI56tLDCYh8XmEnq99H5cq1OF1qujCOoq5MQQuqhYqHpuuzOIBC8bOnR7gzcQoNbCCHkNVzLLMSGuDTE3smFAC8np1eQioVg8PKe3hzP9tTS4xgKPkIIaYD8EjnCL2ci9WkximQKWEglcHFsgnG9nGkgC0dR8BFCCDEqNLiFEEKIUaHgI4QQYlQo+AghhBgVCj5CCCFGhYKPEEKIUaHgI4QQYlQo+AghhBgVCj5CCCFGhYKPEEKIUaHgI4QQYlQo+AghhBgVCj5CCCFGhYKPEEKIUaHgI4QQYlQo+AghhBgVCj5CCCFGhYKPEEKIUaHgI4QQYlQo+AghhBgVCj5CCCFGhYKPEEKIUfk/rSfBGNhI43EAAAAASUVORK5CYII=\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": 33,
|
|
"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",
|
|
" '''Pick a random node y from the neighbors of x in graph with uniform\n",
|
|
" probability.'''\n",
|
|
" y_list = np.fromiter(graph.neighbors(x), dtype=int)\n",
|
|
" return np.random.choice(y_list)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 34,
|
|
"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",
|
|
" n = example_graph.number_of_nodes()\n",
|
|
" Q = np.zeros((n, n))\n",
|
|
" for x in range(n): \n",
|
|
" for k in example_graph.neighbors(x):\n",
|
|
" Q[x,k] = round(1/example_graph.degree(x),2)\n",
|
|
" return Q\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
|
|
}
|