mirror of
https://github.com/ImperialCollegeLondon/ReCoDE_MCMCFF.git
synced 2025-06-26 08:51:16 +02:00
423 lines
44 KiB
Plaintext
423 lines
44 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5ac56056-ca33-4f13-8e36-564b94144c1e",
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"source": [
|
|
"<h1 align=\"center\">Markov Chain Monte Carlo for fun and profit</h1>\n",
|
|
"<h1 align=\"center\"> 🎲 ⛓️ 👉 🧪 </h1>"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"id": "eb5d773e-4cc0-48ae-bb71-7ece7ab5f936",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from numba import jit\n",
|
|
"\n",
|
|
"# This loads some custom styles for matplotlib\n",
|
|
"import json, matplotlib\n",
|
|
"\n",
|
|
"with open(\"assets/matplotlibrc.json\") as f:\n",
|
|
" matplotlib.rcParams.update(json.load(f))\n",
|
|
"\n",
|
|
"np.random.seed(\n",
|
|
" 42\n",
|
|
") # This makes our random numbers reproducable when the notebook is rerun in order\n",
|
|
"\n",
|
|
"\n",
|
|
"def show_state(state, ax=None):\n",
|
|
" if ax is None:\n",
|
|
" f, ax = plt.subplots()\n",
|
|
" ax.matshow(state, cmap=\"Greys\", vmin=-1, vmax=1)\n",
|
|
" ax.set(xticks=[], yticks=[])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "337f1de8-d743-441f-bc15-387bcfff558d",
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"source": [
|
|
"# Doing Monte Carlo!\n",
|
|
"\n",
|
|
"Now that we can evaluate the energy of a state there isn't that much more work to do Markov Chain Monte Carlo on it. I won't go into the details of how MCMC works but put very simply:\n",
|
|
"\n",
|
|
"We want to calculate thermal averages about a physical system. For example, is this bag of H20 molecules solid or liquid at T = -20C? Our Ising model is much simpler so the equivalent question would be what's the average color of this system at some T?\n",
|
|
"\n",
|
|
"It turns out that this question is pretty hard to answer using maths, it can be done for the 2D Ising model but for anything more complicated it's pretty much impossible. This is where MCMC comes in. MCMC is a numerical method that gives us a rule to probalistically jump from one state of the system to another. \n",
|
|
"\n",
|
|
"If we perform many such jumps many times we get a (Markov) chain of states. The great thing about this chain is that if we average a measurement over it, such as looking at the average proportion of white pixels, the answer we get will be close to the real answer for this system and will converge closer and closer to the true answer as we extend the chain. \n",
|
|
"\n",
|
|
"I've written a very basic MCMC sampler for the 2D Ising model below. It needs:\n",
|
|
"- an initial start to start the chain\n",
|
|
"- to know how many steps to take\n",
|
|
"- the temperature we want to simulate at\n",
|
|
"- a way to measure the energy of a state, which we wrote in a previous chapter\n",
|
|
"\n",
|
|
"It then loops over:\n",
|
|
"- modify the state a little, here we just flip one bond\n",
|
|
"- accepting $p = \\exp(-\\Delta_E / T)$ based on how the energy changed\n",
|
|
"- if we rejected, change the state back to how it was"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 24,
|
|
"id": "2586a542-35f2-419e-9aa2-2bb9e9ab74b9",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABuEAAAHoCAYAAABXfkJFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAABDrAAAQ6wFQlOh8AABIuElEQVR4nO3df5xlaV0f+M+XaaZLRgcGWxyEEklYd+KYLBslEuMiZhBF3SVAcHdMNtGgLqImSsjKjzUgiQgY0Y2IZGVdjKu2xrhMEjFBGiGx/RHYdSJOGMENA0UAnVpm+DFQ3QM8+0fdkra6695b95xT59yq9/v16lf1Pfc+z/me73me59zqb597q7UWAAAAAAAAoD/3GTsAAAAAAAAAOG4U4QAAAAAAAKBninAAAAAAAADQM0U4AAAAAAAA6JkiHAAAAAAAAPRMEQ4AAAAAAAB6pggHAAAAAAAAPVOEAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwrI2quqOqWlU9tqf+XjDr7wUrtH31rO039hTL3rF9Xh/9AQAAAAAA41KEg32q6rGzgtgbx46lL1X1ebNjumOAvltVtb77BQAAAACAdXZq7ABgRC9PcjbJ9gptn5PkxUne12tEAAAAAADAsaAIx4nVWtvOagW4tNbeFwU4AAAAAADgAD6OkrV36fezVdUjqupnq+oPq+pCVd1eVd9TVZeN9St9J9zsIyh/bfbwy/c+anH/x1Me9J1wVfUZVfWtVfWaqvqDqvpoVX2kqn6nqp5XVZ/W43E/oKpeVFW3zfazU1Xvqao3VtVzLo01yTtnDx+275juuOR1n1VVf6eq/nVVvXPW3wer6req6tur6qp9+3/BpR9Dua/fyz6esqq+pKrOzmK8WFV3VtW/qKov6ysnAAAAAAAwFe6E4zh5ZJL/Nbt3t/1akgcl+W+y+7GRD03ynUv08a+T7CT5qiR/OHu85/Yl2v9XSf5Jkj9K8vtJ3pLkM5N8SZJ/mOS/q6ovb63tLNHXgarqfknOJ/mC2b5en+SeJA+ebXt0kh+YvfzXk3x6kqfMXvOLl3R16Z2AX5XkR5K8J8k7kvxWkuuT/MVZ/F9ZVU9qre0V2G5N8lNJ/ubs8U/NiffvJvnB2cP/J8lvZvecfG2Sr62qp7fWfmLZ4wcAAAAAgKlThOM4+TtJvi/JC1trn0ySqnpMdgtyz6iql7bWtuZ10Fp7cVX9VnYLUre31r7xkDHckeSmJG/ci2EWxwOS/FySr57F+ZJD9rvfX81use2Xk/yV1trHL9nXVUm+fO9xa+1VVfX67Bbhtucc0/+d5NGttd++dGNVPTjJa5M8McnXJ/n5Wb+vSfKaqvqbs8dX7LeqnpDkHyV5b5InX9p/Vf2lWd8/VlVvaq29fcnjBwAAAACASfNxlBwnb07yfZcWv1pr/zbJv8nuWP+KoQNorb2ntfaGS2OYbb87yd+ePfyrPezqs2c/X39pAW62r0+01t5w2A5ba2/bX4CbbX9fkv959nCV2F8w+/nN+/tvrZ1P8g+S3DfJ/7RC3wAAAAAAMEnuhOM4ee0lH5V4qduTPCHJ5xxFEFVVSf5Sksdk9yMXPy1Jzf4kyef3sJs3z35+T1VtJ/lXs0JfJ1V1Kslfzu5HUF6fZCO7cX/G7CWHir2qziT5C0k+lOR1B7zsTbOff/Gw8QIAAAAAwFQpwnGcvPuA7R+a/dwYOoCq+uwkv5TkS+e87Nqu+2mtvbGqXprkWUl+Okmrqtuz+/1v/7y19m8O22dVfX6S1yT5M3NedtjYH35Ju4/v1icP9FmH7BsAAAAAACZLEY7j5JOLXzK4V2W3AHc+ux/D+B+S3N1au7eqrk5yoa8dtda+p6pemd3vavuy7N599y1JvqWqXpfka/d/VOUCv5jdAty/SPLSJG9L8sHW2idmBbrfz6fu5lvWVbOfH8xugW+e7UP2DQAAAAAAk6UIBz2pqmuSfE2STyT5uit8POQj+t5na+2dSX5k9idV9WVJfi7J45P8rST/2zL9VNUNSf5skj9K8uTW2if2vWTV2LdmP+9trX3jin0AAAAAAMDauc/YAcAEXZz9PGyR+v7ZnVMfPuD72f5al6CW0Vr79SSvnj38ry55atExPXD2871XKMAl82O/N/nj75PbH89/TvLWJGeq6rFz+gAAAAAAgGNFEQ4u959nPx9xpcLSHH+Y5K4kD6iqb7j0iar66iTP7Cm+VNWTquoxVXWffds/LcnjZg/fdclTd2a3EPfZVXXdFbp8R3Y/zvMLq+ox+/r8piQ3zwlnL18HfZfc985+/p9V9fgrHMtVVfWXq+rRc/YBAAAAAABrRREO9mmtvSvJ7yT57CS/W1U/XVWvqqq/t6DdJ5J8/+zhz1TVb1TVz1bVbyf5lSQv6zHML0/ypiTvr6p/U1X/Z1X9yyTvSfLoJLcn+SeXxHZvkl/O7p1wv1NVPzM7phfPnr8zyStmz/9aVb1hFvtbk/xkkhfPieX/mv08V1VnZ/2+6pJ935Lk7ya5Psm/qarfr6p/Mev/Ddn9LrhzSR7ZOSsAAAAAADARvhMOruzJSV6S3WLXzUmuym7R6wfnNWqt/VBV3ZHkWUluTPKFSX4vyV9vrf1MVT23p/henWQnyZfN9nEmyd1J/iC73wn3v7fWPryvzbck+UCSr0ry9dmd/+9K8uzZ838nye8m+bYkfyG7HzP5fyf5e9kt6j3ngFiel6QleVJ283bf2fZv3ntBa+1lVXUuyXcmeWySr0zy8STvS/Jvk/zLJL90iOMHAAAAAIBJq9ba2DEAAAAAAADAseLjKAEAAAAAAKBninAAAAAAAADQs5WKcFX1iKp6ZVXdWlUfr6rfW7JdVdWzq+rdVfWxqvrNqnr0KjEAAAAAAADAVK16J9yNSb42yR8k+Y+HaPc9Sb4vyQ8n+bok70vyuqr6UyvGAQAAAAAAAJNTrbXDN6q6T2vtk7O/vzrJF7fWvnBBm40kf5jkx1prz51tuzrJ25O8trX2jEMHAgAAAAAAABO00p1wewW4Q/rSJNcm+YVL+rmY5JeSfM0qcQAA0J+qaiv8efXYcc9TVd96SawvHzseAKZh3a95VfWCBbHujB0jAEfnpF/XquobqurfVdUHq+ojVfWWqvr2qlr1kwChN6eOcF83zH7evm/725J8blV9WmvtY0cYDwAAf9JPXWHb9Um+Ksk9SX7xCs//+qARdVBVD0vyj5K0JDVyOABMy3G55v2HJLdeYfu9RxwHAOM6sde1qvqxJM9IspPk3Oy1NyV5eZKbquqvrnhTEfTiKItw1yW50FrbX7W+K7v/KHJdkisW4arq2uzeRZck+YIv+IKrv/u7v3vz4Q9/+B1JPjFMuAAM7Kokn5XkrTfddNOFsYMBktbaN+7fVlWPze4vbttXen6qqqqS/O/Z/eSHf5rkb44bEQBTcoyuea9prb1g7CAAGNdJva5V1VOyW4B7f5LHtNbeMdv+2Ul+LcmTknxnkv+1/1BhOUdZhOvimUmev/fg4x//eB7+8IePGA4APXpUkreMHQRw7Dw9u//78W8n+cyRYwEAAKB/z5n9/J69AlyStNb+sKq+Lckbkzy7qn7U3XCM5SiLcHclOV1VG/vuhrsuux8RdNecti9L8qq9BzfddNNDkvzWQx/60NzvfvcbJNjj4OLFi7ntttty44035uqrrx47nMmSp+XI03LkaXkf/ehH8573vCdJ7hw7FuB4qaqHJ3lpdj9e5eW55D9zAQAAsP6q6qFJvijJxST/bP/zrbU3VdV/TvKQJI9O8htHGyHsOsoi3N53wf2X2f1s1z03JHn3vO+Da619KMmH9h6fO3cuSXK/+90vn/7pn95/pMfEhQsX0lrLNddck9OnT48dzmTJ03LkaTnytBIfKwz0ZvYxlD+Z3fe5T2uttd1NAHAs/fmqekl2/4PzB5L8dpJfbq1dHDcsAFjJYa5r//Xs521zagtvzm4R7r+OIhwjOcoi3G9kt5D21MyKcFV13yRPTvLaI4wDAIAjVlVthWZvaq099pBtviPJY5M8u7X29hX2CQCdHOE1L0n+29mfS72nqv56a+1NK/QHAH/ChK9re99X9a45/b1732vhyK1UhKuq+yX5mtnDhyW5tqr+6uzxm1prd1bVuSQPa609IklaaztV9QNJXlBVdyZ5a3a/NPEzk/yjLgcBAMDk/dQKbW5f/JJPqao/neTF2f2eSe8vARjL4Ne8JP9vdr8H51eSvDPJ1Un+bHY/gvnLk7y2qv5ia+13V4gFAC411eva3kfk3TOn34/Mfn7GIeOB3qx6J9yDcvnnrO49/orsfuHhVVfo/yVJKsmzknxWkluTfFVr7T+tGAcAAGugtfaNQ/Z/ycdQ3je7H0Ppo24BGMXQ17zZPn76Cpt/LcmvVdUvJnlKkhcl+bqhYwHgeHNdg27us0qj1todrbU64M8bZ695bGvt8/a1a621H2itbbbWNlprj26t/Wb3wwAA4IT720kek+QH/K9/AE64F85+fuXsa0AAYJ0ddF3bu8vtmjlt9+6W+3DvUcGSjvI74QAAOKGq6tUrNLu9tfbiJV/7pNnPr6yqL9/33OftvaaqvjDJR1pr/gclAIM4gmvewr5mP69OcibJ+3rqF4ATaMLXtTtmPx82p+3mvtfCkVOEAwDgKPzNFdq8Kbvf8XYYf3HOc58z+/PBFWIBgGUd1TXvIJ95yd8/cuCrAGA5U72u/c7s541V9WmttY9doe2j9r0WjpwiHAAAg2ut1cD9P/ag56rqBdn9Qu8fa619x5BxAMDQ17wlfP3s5++31nz8FgCdTPW61lrbqqr/J8mfT/LUJP/00kazT0h5aJL3J/GVWIxmpe+EAwAAAODoVdXnVtU3VNXpfdurqv7HJD8w2/TDRx8dABxOx+va3nMvqapHXNL2QUleMXv44tbaJ/uOG5blTjgAAACA9fHAJD+T5JWzOwDem+QzktyY5OGz17y8tfZPRooPAA5j5etaa+0Xq+rHk3xbkrdW1euT3JvkpiTXJnlNkpcPfgQwhyIcAAAAwPrYSvKD2f2em0ck+QvZ/aSj9yf5+ST/W2vtDeOFBwCH0um61lp7RlX9epJvT/LlSa5KcnuSn0zy4+6CY2yKcAAAHKi19sYkY38HQCettRckecHIYQAwcetyzWut/X9J/uex4wBg2k7Sda219rNJfrafiKBfvhMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADo2amxAwAADnbu3LnTSf5skjuTfGLkcAA4vKuSfFaSt950000Xxg5mylzzANaea94SXO8AjoWlr3mKcAAwbX82yZvHDgKAzh6V5C1jBzFxrnkAx4Nr3nyudwDHx8JrniIcAEzbnUnyxV/8xTl9+nSvHW9ubvba39i2trbmPj/veOe1nWqe1jFm1tOiucWnXGnuPfCBD8wrXvGKZLaeM1ena96q6/yitl10mT/HbS0/SWvJFM/dScr/WLqc93U8P655nUzyejfUOFw0N6b4e81YuZhnrDyNkYsh3z+tmsex3tNNcSwOZd3er1+4cCFvectbkiWueYpwADBtn0iS06dPZ2Njo9eOt7e3e+1vbIvyM+9457Wdap7WMWbWU99rz3G2YO75uKnFOl3zVl3nF7Xtosv8OW5r+UlaS6Z47k5S/sfS5byv4/lxzetkkte7ocbhorkxxd9rxsrFPGPlaYxcDPn+adU8jvWebopjcSjr+H59ZuE17z5d9wAAAAAAAAD8SYpwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD07NTYAQAA9KGqRmk7lnWMeQittbnPy1N3i3K46BwMtd9VGTMn01jndd5+h5o7U3XccjHFtWKMPE4xD8k4uehyfTlu84NuhpxX88bTFK+VQ+kyX4eck2Ocny7Hs45r16oxd/mdaMhcTPEcrJrHLmN8jByfOXMmZ8+eXaoPd8IBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9OzV2AAAAsKqqGjsE5pji+ZliTPxJm5ub2d7eHjuMwU1xLLbW5j4/L+ZFbVftt2vfq+530T7nPT/WuR3r/EzRULkYyhhjvOt+r9R2Z2cn58+fX7nPk+ag612XtbiLoa4BXeJdx7nR5XinuI6v2zWgyzieYp6GmgNTvBZ2MeQ4PShXh7nmuRMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU2MHAABwXLXW5j5fVUcUyfLmxTzFeBmfcQFHbx2vL13MO55FuZinS9uTlOPjZqjxtMiQfU9pnyxva2srGxsbl20fcj6OMdeHXGvH+N2lS79TvH6Ptf5PcX0a633BULnocr2b4vF0Mfa/c7gTDgAAAAAAAHqmCAcAAAAAAAA9U4QDAAAAAACAninCAQAAAAAAQM8U4QAAAAAAAKBninAAAAAAAADQM0U4AAAAAAAA6NmpsQMAABbb3NzM9vb2Zdtbayv3uahtVa3cN7vWMYddYl51PA65z3U8B3DSbW1tZWNj47Lt5vOnyMWnyMVy5Gk5i/LU5b33qqZ67sbIBYuN8V5+kSHn1by+pzhGh8pFl36HGjNjrV2LcrhqzEOOpzFyddLm5dgxuRMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU2MHAACsrqrGDoGBtNbmPj/UuV+031WNNVbnHY/5Mz7nh8Posj4ZT/1wDoA9V5rTZ86cydmzZ0eIZj1tbm5me3v70O3Gev80VN9d+l01F0P+rjXFa+UYOR7SUL+zdul3qJiG0uXcDdW2y7xc1Hbs8+NOOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGenxg4AAIDLVdXk9ttaO8JIljNWnuiH88dR6bJ+LRqn8/petzHeJd4hrxHHKcecXOv2HovxTHU9HSquKa7jUzzWdYtprPdPxy3Hi0wxj/N0yfFQx3oU12B3wgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAenZq7AAAAKastTb3+ao6okj+pEVxAVc21TnNp2xubmZ7e/uy7WOte132O6/tkGNtjP0u6rdLHqc4L1c9nuN23ruwHnPSbW1tZWNj47LtQ66nXfodak6u29rVRZccj5WLVfc75Pu2eTENNY7XcV4ONbfGek83b79jjLczZ87k7NmzS/XhTjgAAAAAAADomSIcAAAAAAAA9GylIlxV3VBVv1pV91TV+6vqpVV19RLtPrOqXllV7561/b2qevoqMQAAAAAAAMBUHfo74arquiRvSPKOJE9O8pAkL0tyvyTfsaD5P0tyQ5LnJnl3kq9J8uNV9YnW2k8cNhYAAAAAAACYokMX4ZI8Pcm1SZ7UWvtAklTVqSSvqKoXtdbee6VGVXV9kq9I8k2ttVfPNr+hqh6V5H9IoggHAAAAAADAsbDKx1E+Icnr9wpwM78w6+vxc9rdd/bzg/u2fzBJrRAHAAAAAAAATNIqd8LdkOQnL93QWru7qt43e+6KWmtbVfW6JM+tqt9PspXdgt7jk/y1FeIAABhc1TT/r9C8uFprK7VbZF6/XfuGo2Kcrq9F527RGjWUdRtTXfLU5VhPUp7G6DdZ/b3BorZdDPWe5CSZ6vuvK8W1s7OT8+fPjxDNetrc3Mz29navfXYZL4varjqfpzqGxzDUe5ku/Q7ZdgxDvVcZcm6t2m+Xvoecl0ONi6FyfBRWKcJdl+TuK2y/K8kDF7R9cpKfT3Lb7PEnknxna+2fz2tUVddm9yMwkyRPe9rTrr/55ptz8eLFXLhwYdm4T5y93MjRfPK0HHlajjwt7+LFi2OHAAAAAAAwmFWKcCup3VLl/5Hkv0jyDUnel+Qrk/xIVd3VWjs7p/kzkzx/78Ett9ySm2++Obfddtvkq5xTcOutt44dwlqQp+XI03LkabGqyjXXXDN2GAAAAAAAg1ilCHdXkvtfYft1ST5whe17vjbJU5P8udbaW2fb3lhVD0ryQ0nmFeFeluRVew+e+MQnXp/kzTfeeKN/wJ3jwoULufXWW/PIRz4yp0+fHjucyZKn5cjTcuRpeffcc0/uuOOOscMAAAAAABjEKkW427Pvu9+q6v5JHjx77iBfkN2Pn/y9fdt/J8k3V9X9WmsfvVLD1tqHknxo7/G5c+eSJFdffbV/5F7C6dOn5WkJ8rQceVqOPC127733jh0CAAAAAMBg7rNCm19J8riqesAl256a5JNJXjen3buSXJXkz+3b/kVJ/uigAhwAAAAAAACsm1WKcK9M8uEkr6mqx1fVNyX5wSSvbK29d+9FVXWuqv7gknavTfLuJL9YVX+9qm6qqpck+cYkP7ryEQAAAAAAAMDEHPrjKFtrd1XVTdktnL0muwW5VyV53r6XXnVp/621D8/afX+SlyR5QJJ3JnlmkpevEDsAAAAAAABM0irfCZfW2tuSPG7Bax57hW1/kOS/X2WfAHCSbW1tZWNj47LtVTVCNEzdUOOiS7+ttUH6BdhjLfkUuVjOvGvTWOaduy7xLhoTrtO7FuV4jFycpPyznLHe68+bH13WkKHWn6HWzC7rxJBthzJGTGOtxUONxS5zYIocz66dnZ2cP39+qX2s8nGUAAAAAAAAwByKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9OzV2AADAYpubm9ne3h47DI5Qa23ltlXVYyT9mGJMi8w7B+t4PPMMNd4W9Xvc8kg/tra2srGxcdn2qY6XodaKLvNynqHyOFS8yXTP/UG6rItd8thlLK6a46mu86vmYqx4122MM6wh19MuVl3bpng86zjXx2q7ar9D/j7b5Xq3br/jdbnOjvUedYrvNcc+t+6EAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAenZq7AAAALhcVY0dAidIl/HWWhukX1gXQ43zef3Om3d8ypB5WvW8D7Xedu17VVNd56ca11Gb4phhXFO8toxxHZ1iv8lw56DLeR9qzAwVU5f9jmWMeTnF95JDjsUhxuqZM2dy9uzZpfpwJxwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGenxg4AAACmqKrGDgFOpM3NzWxvb1+2vbW2cp/mcz+6nAO6O27jeNF4Om7HOwY5nLatra1sbGxctn0dz9tQMY913Vm3612X/C9qO0Yu1nEODOW45aLLeJrXtkueFsV00PM7Ozs5f/78UvtwJxwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGenxg4AAFg/rbW5z1fVEUUCTNmitWIe6wj7GRPLGypX8/qd4nyf4pjxHupTTtKxrqMuc5phzTs3i+ZVl/PaZb+r9rvIqtelRfGOtT4NleNV97nIUO8Lulwrh7rODhnTqv12sW5zdpEu52fIdXOPO+EAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICenRo7AABg/VTV2CEAE2E9YCpaawc+t2icdmnLLnlazjrmyfw4mbqsm3TTJbeL2naZs/P6Hmo8DDUOhxy/Q62LY+V4qH7Het82xto15LVy3dbiKb6nOIocuhMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU2MHAAAAcBy01i7btrOzk/Pnz48QzclTVYP0e6XzuqyhYoKjZBzD0Vo05+Zdl4acr6v2veg6OsYaM1aOp5iLRVZ9H7SOeZrX76KYurxf7GKKY2aesd6vj50nd8IBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9OzV2AAAAwLhaawc+V1VHGAmsbt44HssUYwJgGjY3N7O9vX2k++xyXZr3nrDLe8kubVd9nzrW9XmK76uHGhNd2naJaVHbofY71rkd6ve4MfpdZKz99sGdcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9U4QDAAAAAACAnp0aOwAAAGA5rbUDn6uqlfvt0vYkmZd/+rO1tZWNjY2xw5g0c3Y5i+bsFPM41Do/lHXM8brpcu2R//XU5ZyPNSfn9TvW8Yw1d07SnB3rd5Oh3pN3GU9Dnfd1nAPr5ijGmzvhAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9U4QDAAAAAACAnp0aOwAAYLGtra1sbGwcul1VDRDNydJam/u8HHOUjLdpu9L5OXPmTM6ePTtCNOtrc3Mz29vbh263aL0eylDzcqzjOU7GWjPnnbvjto6v4/Gs2/mZYkz046Df8Rad8y7Xhy5txxiLQ+ZiqH675GnV9WmK7xm6/B49ZP5XzfE6rsVTfB80T5d4u5z3VcfqYX7PcyccAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9W6kIV1U3VNWvVtU9VfX+qnppVV29ZNuHVNVPVdWdVfWxqnpbVf21VeIAAAAAAACAKTp12AZVdV2SNyR5R5InJ3lIkpcluV+S71jQ9sFJfjPJ7yf51iQfSnJjktOHjQMAAAAAAACm6tBFuCRPT3Jtkie11j6QJFV1KskrqupFrbX3zmn70iRbSb66tfaJ2bZzK8QAAAAAAAAAk7VKEe4JSV6/V4Cb+YUkr0zy+CSvvlKjqro2ydcn+VuXFOAAgA6qau7zrbWV27JLnoBlXWnN3dnZyfnz50eI5uSZt17Pux52NWTfq+7zJF27pvhep8t+nbvlyDHHwebmZra3ty/bPtbc6LLfoa6FY1xjk+HyONb5mafLvyl06XcsU43rIGONmaHeV0/xvdlRrDOrfCfcDUluv3RDa+3uJO+bPXeQP5/k6iT3VtWbqure2ffJvaSq7rtCHAAAAAAAADBJq9wJd12Su6+w/a4kD5zT7vrZz1cl+YkkL0jyF5K8MMknkzznoIazu+iu3Xv8tKc97fqbb745Fy9ezIULFw4T+4mylxs5mk+eliNPy5Gn5V28eHHsEAAAAAAABrNKEW5Ve3fdvb619ndnf/+1qvqMJM+qqhe21j52QNtnJnn+3oNbbrklN998c2677bbRbkteJ7feeuvYIawFeVqOPC1HnharqlxzzTVjhwEAAAAAMIhVinB3Jbn/FbZfl+QDV9h+abskecO+7eeSPC/JI5K89YC2L8vuHXRJkic+8YnXJ3nzjTfe6B9w57hw4UJuvfXWPPKRj8zp06fHDmey5Gk58rQceVrePffckzvuuGPsMAAAAAAABrFKEe727Pvut6q6f5IHZ993xe3zHxf0u3HQE621DyX50N7jc+fOJUmuvvpq/8i9hNOnT8vTEuRpOfK0HHla7N577x07BAAAAACAwdxn8Usu8ytJHldVD7hk21Oz+71urzuoUWvtXdm90+1x+576yiQfy+IiHQAAAAAAAKyFVe6Ee2WS70zymqp6UZKHJPnBJK9srb1370VVdS7Jw1prj7ik7fOS3FJVP5Lkl5M8Ksmzkry0tXbPaocAAMff5uZmtre3D92uqgaIBg5v3vf4GqcAh9Nl3ezyvepTXK/Hisl1bdei8TQvF8ctT11ysY6udDxnzpzJ2bNnR4hmPW1tbWVj48APBjvQvLE05Jxcdd0b8rrTpe8xDJWLsfIw1PuRIc/7UNfvMfpdpMuYGer8THEOrHo8Ozs7OX/+/FL7OPSdcK21u5LclOTjSV6T5MXZ/b62Z+576VXZV+Rrrf3LJDdn9264f5XkW5M8P8n3HjYOAAAAAAAAmKpV7oRLa+1tufxjJfe/5rEHbP/5JD+/yn4BAAAAAABgHazynXAAAAAAAADAHIpwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAenZq7AAAANZZa+3A56rqCCOZNrngOFhlHJ85cyZnz54dIJrja2trKxsbG4duN+/8LDp389byLvsdap9d2nbhmvcpJ+14DyIP4xtj/aIfm5ub2d7ePnS7sdbiVfvuEtM6ju8xjneK73Omat7xjDW3hjo/XY5nXttF8Y7x3niMte8wv+e5Ew4AAAAAAAB6pggHAAAAAAAAPVOEAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0LNTYwcAALDOqmrsEDhBWmsHPmcsfoo8nUxDnfehxkyXfucd65BO0vwZKsfHLYdjjcUupjin51mU4+M2pk6Sra2tbGxsXLZ9yHPaZc7Oi2uofhflYt5+x1rHu7wfWfXcT3Et7hLTFI9nkaHmRxddYurSdt3O31FcZ90JBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9OzU2AEAAIttbW1lY2Pjsu1VNUI0cGWttQOfM1aBMc1bnxJr1NRN8foyxTEzRp4W9bto7gGHN2/edbneLWo7xhrTZQ2Z4jo9lC5r8VBj5rhdH7rkqUsuuozjLjmeYkxD9XsUa4U74QAAAAAAAKBninAAAAAAAADQM0U4AAAAAAAA6JkiHAAAAAAAAPRMEQ4AAAAAAAB6pggHAAAAAAAAPTs1dgAAwGKbm5vZ3t4eO4zBtdbmPl9VRxQJq3B+hifHu6wVx9tB1zzn/VNO0rEmJ+94VzXFPM2LadGcHkqX/XY5nlXPzxTPK/1Y9Xo3lEVjbV5c857r0u9YuhzPUHN2qDwNdX6GWmsX9T3W+8UpjpmhdMlxl+tolxyPzZ1wAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICenRo7AABgmlprBz5XVYPsc6h+FxnjWKdqXi4WmZcrOaZPxgwwNOtMP7q8rxjKqu9Xulq1b2Px+Nra2srGxsZl26d6zoeaO136XTVXU5zryerHs6jdGGvxWDF1mT9jjMVFfY91PKv2u6jvKa5vQ57bPe6EAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAenZq7AAAgGmqqrFDODJTPNbW2tzn1zFmupuX4ymOCVjF1tZWNjY2Dt1uqPmxjuvxuuly/ThJ+Xed7cdQeVw0Fp0/jsKQ16xVx/CQY3+MeTXWNWuotmO9z5liTEPts8uYmeI1a6yYhrJqLnZ2dnL+/Pml9uFOOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAz06NHQAAwHHVWpv7fFWt9NyQxtrvGLqcn7FMMSbo2+bmZra3ty/bvmjODuW4zbt5eZzqtWesc3+SzDsH8j+847bOMF1drgHrtk50mVdDHc86/v4xxXPbJaahxvFQ526ofsc6r1McT0cx79wJBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9OzU2AEAABxXVTV2CMzh/MA0bW1tZWNjY+wwjq15a19rbeW2rLdF557F5JApGGud7nJtGcOQMQ2Vi3lt1/G8D5WnRbmY4nhct/mzyBjndsjxdlDbM2fO5OzZs3Pb7nEnHAAAAAAAAPRMEQ4AAAAAAAB6pggHAAAAAAAAPVOEAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ6fGDgAAgH611g58rqqOMJJPmWJMwHqxVgxvqjmealws1uXczXvvMJYuMRnH7Le5uZnt7e0j3edQ82qo8T1Wv1PM03HLxRT7HeqaNeR5X3W/U7yeDTkWh2i7s7OT8+fPL9WHO+EAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAerZSEa6qbqiqX62qe6rq/VX10qq6+pB9fFdVtar6V6vEAAAAAAAAAFN16rANquq6JG9I8o4kT07ykCQvS3K/JN+xZB/XJ3l+kj867P4BAAAAAABg6g5dhEvy9CTXJnlSa+0DSVJVp5K8oqpe1Fp77xJ9vDTJv0jysBX2DwDAHFU1dggAK9vc3Mz29vbYYbAmWmsrt12362WXeLvkaShTjGkdzcvjuo1xuls0r4ZaR07SWrzIGHNyrPV0qP0Oeb0bY7wNOT/mPT/F6+xQ8Q55XvvI4yofR/mEJK/fK8DN/MKsr8cvalxVX5bkryR59gr7BgAAAAAAgMlbpQh3Q5LbL93QWrs7yftmzx2oqq5K8vIk399ae98K+wYAAAAAAIDJW+XjKK9LcvcVtt+V5IEL2j4jyTVJfvgwO6yqa7P7EZhJkqc97WnX33zzzbl48WIuXLhwmK5OlL3cyNF88rQceVqOPC3v4sWLY4cAAAAAADCYVYpwK6mqByV5YZK/0Vo77L+8PjPJ8/ce3HLLLbn55ptz2223TfKzTafm1ltvHTuEtSBPy5Gn5cjTYlWVa665ZuwwAAAAAAAGsUoR7q4k97/C9uuSfOAK2/e8MMnvJvl3VfWAS/Z/avb4I621jx/Q9mVJXrX34IlPfOL1Sd584403+gfcOS5cuJBbb701j3zkI3P69Omxw5kseVqOPC1HnpZ3zz335I477hg7DAAAAACAQaxShLs9+777rarun+TB2fddcfvckOQx2S3i7XdXkick+ddXatha+1CSD+09PnfuXJLk6quv9o/cSzh9+rQ8LUGeliNPy5Gnxe69996xQwAAAAAAGMwqRbhfSfLcqnpAa+3u2banJvlkktfNafddSR6wb9uPJPlYkudk9y45AACOoaoaOwQAThjXnuWNkasu+1zHryZZt5gXxWt+dbe1tZWNjY3Ltg+Z23nndaj9Lup3qJi6zLl5+x1rbnTJ01i5GMpYuZhnimviFMdxF+v8vmGVItwrk3xnktdU1YuSPCTJDyZ5ZWvtvXsvqqpzSR7WWntEkrTWbt3fUVXdnd2PoXzjCnEAAAAAAADAJN3nsA1aa3cluSnJx5O8JsmLs/t9bc/c99KrslqRDwAAAAAAANbaSkWy1trbkjxuwWseu0Q/C18DAAAAAAAA6+bQd8IBAAAAAAAA8ynCAQAAAAAAQM8U4QAAAAAAAKBnK30nHAAAi7XW5j5fVUcUyfLmxTzFeIHjZ2trKxsbG5dttwZxJUONi3W8hs+zbvEmi88Bw1rHMbNuNjc3s729feh2Q82Nsda9MfpddKxTXH+6HE8Xq/5+2CWmsdrOs2icdvk9eqy2UzNkvGPNnz3uhAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHp2auwAAACOq6oaO4RDW8eYj5PW2tznnR9OMvODo2Q8jW/eOVi0HqzabxddYpqiVY5nZ2cn58+fHyCa42lraysbGxuXbV80RoeaG13M2++i4xmj7VTX+C65mGeMde+45bjL3JrivByy7arGGjNd1plluRMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU2MHAACMo7U29/mqOqJIgD2L5t28eWvOAjCEoa49Xd6LHrdr3hSPp8t7EpazubmZ7e3tQ7frkvt553VRv0OtBVMc//MM+Xv0ULkY49yNtUYsOp4pjreh8jjkNXpVQ8U05Hk9qO8zZ87k7NmzS/XhTjgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9U4QDAAAAAACAninCAQAAAAAAQM8U4QAAAAAAAKBnp8YOAABYbGtrKxsbG5dtr6qV++zSFjhYa23ltovmpXkLB5s397rMnUVzeox5OcWY1tFQY6aLLteQMRiLn3KSjjW58vGeOXMmZ8+eHSGa9TTG73hDrTFD9TvUvBpyrR3qd4EpXrPmGfK911Bt58XcZW6NNS+77HPVsdil7ZDjeIgc7+zs5Pz580u91p1wAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICenRo7AABgsc3NzWxvb48dBpworbUDn6uqI4wEWMZQ83KK832KMa2jMfI479qyyJDxdolrjH4Z3irnbmdnJ+fPnx8gGqZg3hrU5X3zvLaLxuFQMXUx1Lq36rEuajuUsdb/sd4jTXG/Xc7BGG2H3OcQ5+fMmTM5e/bsUq91JxwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU2MHAACMo7U29/mqOqJIYP0smj/HibWCo7a5uZnt7e2xw+jFUGuFebe+pnru5sV1kq55Y5mX47HGzKL9GhfDGSu3UxxrQ+Viqjme4rxaNaYu46lLHrq0HetaONR+hzoHU+y3y34X6SMud8IBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADomSIcAAAAAAAA9EwRDgAAAAAAAHp2auwAAIDFtra2srGxcdn2qhohGjgZ5s2v1toRRjIu68zyrjQudnZ2cv78+RGiOXnmzcuxxnGX/c47nrHWIOvByeS8D2+sHJ+k9zNTtLm5me3t7cu2dzkvi9oOdV3q0m5eTIviXTWmLv0OmcOhfv/o0nbV453q+jLGettln0ONt0X9jpGnIefH2L/buxMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU2MHAACMo6rGDgHWlvnDlVxpXJw5cyZnz54dIRq4srHWr9bagc9ZU4GTaGtrKxsbG4du12XNnLcWHzdd8jSv7ZA5HONa2eV4hop3rBzPM+R7lS7HOy+uoc7Bon5XbbsoD0ONi6N4H+pOOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGenxg4AAFhsc3Mz29vbY4fxx1prK7etqh4jAYBdJ+n6so7Huo4xM6xF7yfXbcx0OZ55bbvkoct7dqary3hZdRwuMtR8HWudGGruLIp3qP0Otcasuk+WN8XzM1TbIa93Q4zHnZ2dnD9/fqnXuhMOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANCzU6s0qqobkvxoki9N8uEk/zTJ/9JauzinzYOTfHeSxyf500k+mOTfJnlOa+1dq8QBAIyjqsYOAdZWa+3A58yt8c07P4s4f+M6bnNrjJhP2vg/acc7Nccth12O57jlgm4WjYcua1eX/a5qqHiH7HteLrrsc8hcDGWK422M6/dUz90Yca3jWjHUnF7WoYtwVXVdkjckeUeSJyd5SJKXJblfku+Y0/SLZq//ySS/leRMku9N8u+r6gtba3ceNhYAAAAAAACYolXuhHt6kmuTPKm19oEkqapTSV5RVS9qrb33gHa/nuSG1trH9zZU1W8keXeSv5Hkh1aIBQAAAAAAACZnle+Ee0KS1+8V4GZ+YdbX4w9q1Fq7+9IC3Gzbe5LcmeRzVogDAAAAAAAAJmmVItwNSW6/dENr7e4k75s9t7Sq+vwkD0rythXiAAAAAAAAgEla5eMor0ty9xW235Xkgct2UrvfhvePk7w3yc8teO212f0IzCTJ0572tOtvvvnmXLx4MRcuXFh2lyfOXm7kaD55Wo48LUeelnfx4sWxQwAAAAAAGMwqRbi+vCDJTUm+urV2z4LXPjPJ8/ce3HLLLbn55ptz2223pbU2YIjHw6233jp2CGtBnpYjT8uRp8WqKtdcc83YYQAAAAAADGKVItxdSe5/he3XJfnAFbZfpqq+JcnfT/K01tq5JZq8LMmr9h488YlPvD7Jm2+88Ub/gDvHhQsXcuutt+aRj3xkTp8+PXY4kyVPy5Gn5cjT8u65557ccccdY4cBAAAAADCIVYpwt2ffd79V1f2TPDj7vivuSqrqSUl+PMnfb6395DI7bK19KMmH9h6fO7dbt7v66qv9I/cSTp8+LU9LkKflyNNy5Gmxe++9d+wQAEax+6nsjMUnaRxf5ta4hpxbq57bIWOa1/dQY7HL8ZgfJ9Oi8+6auJ7mnddF53SMtaDLODxpY/gkrdVDnbsuORzrOjvGOB5rrK1jjg/a75kzZ3L27Nml+rjPCvv9lSSPq6oHXLLtqUk+meR18xpW1WOz+/1vP9Fa+wcr7BsAAAAAAAAmb5Ui3CuTfDjJa6rq8VX1TUl+MMkrW2vv3XtRVZ2rqj+45PGfSfKaJO9I8tNV9ehL/vzpTkcBAAAAAAAAE3Loj6Nsrd1VVTcl+dHsFtU+nN3va3vevpdeta//L8nud8ndP8n5fa/9qSTfeNhYAAAAAAAAYIpW+U64tNbeluRxC17z2H2PX53k1avsDwAAAAAAANbJKh9HCQAAAAAAAMyhCAcAAAAAAAA9W+njKAEAgP611gbru6oG63udLMrDkOeA5WxtbWVjY2PsMP7YcZs7J2mMn7T5fpLG6joe63E7nivFfObMmZw9e3aEaNbT5uZmtre3e+2zy1jqsibO2++ifrvE3GW/Q+yzq6HWiaFiHmu8dTHUfrucu3WbA4uMkeMuVn2/uLOzk/Pnzy+1D3fCAQAAAAAAQM8U4QAAAAAAAKBninAAAAAAAADQM0U4AAAAAAAA6JkiHAAAAAAAAPRMEQ4AAAAAAAB6dmrsAAAAAPZrrR34XFUNtt95fc+LiZNrrLE6j7G6nLHOz6qmGO+isdYl5ikebxerHo/5fHxtbW1lY2Pjsu2LxkqXMTHGvOpyPEPmYox+p2isHA71nnuo69KQ17sp7reLMXI8RtszZ87k7Nmzc9vucSccAAAAAAAA9EwRDgAAAAAAAHqmCAcAAAAAAAA9U4QDAAAAAACAninCAQAAAAAAQM8U4QAAAAAAAKBninAAAAAAAADQs1NjBwAAAOyqqrnPt9ZW7nte20X7HcNQMXXJ4Sp97+zs5Pz584Ptk+ENOWZWNVRMU1wLjpsu526K52eKMXUxxfnO8bW5uZnt7e1e+xxyTs7re4rvM1eNd1HbsYx1/Zjie/Iu53YMxy2mLmOiy1qxzu9/3QkHAAAAAAAAPVOEAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOjZqbEDAAAW29raysbGxmXbq2qwfbbWDnxuyP2uG3liXZyk8ThvXg7pSjk+c+ZMzp49O0I0HJUx5tZJms+LnKTr8KK17bgd7xgW5bDLeBvr2jSPMXP8dFknhhrD6zj2h4q5y5wbar6OtQ5M8XjGWuPXLRdd2s6Laarvcw6Ka2dnJ+fPn1+qD3fCAQAAAAAAQM8U4QAAAAAAAKBninAAAAAAAADQM0U4AAAAAAAA6JkiHAAAAAAAAPRMEQ4AAAAAAAB6pggHAAAAAAAAPTs1dgAAwOpaa3Ofr6qV++7SdgxD5mKMfjmZxhrH9ONK529nZyfnz58fIRqOyrx5a84Ob91yvChe4wlOjq2trWxsbIwdxh9b9D50nqHWpy4xdWk7lLFyPNS1ZR3Pzxj/LrDoWMfIxaI8dDme4/Z+5aDjOXPmTM6ePbtUH+6EAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD07NTYAQAAi21ubmZ7e3vsMCatqsYOAZbSWjvwOeO4H/PyOC//Q+z3zJkzOXv27GD7ZNq6jDfrwcnkvK+vIa8vcKku73OGeh/apd+hrpXm5HLGyv+6jbej/h1iWavm4rjNjyHHWx/cCQcAAAAAAAA9U4QDAAAAAACAninCAQAAAAAAQM8U4QAAAAAAAKBninAAAAAAAADQM0U4AAAAAAAA6JkiHAAAAAAAAPTs1NgBAADAummtHfhcVR1hJEzJlcbFzs5Ozp8/P0I062tzczPb29uXbZ8374ZkTjMVi+bAcRurrrXdydO0HXS966LLOe9ynZ233y5r16K2q64TQx3rImOt40Ploss6PdQaP9R4G+t4jpux8jTEuT3M73nuhAMAAAAAAICeKcIBAAAAAABAzxThAAAAAAAAoGeKcAAAAAAAANAzRTgAAAAAAADo2amxAwAA5roqSR74wAeOHQdwiZ2dnQOfO3PmzGBt6W5e/odw4cKFvb9edaQ7Xk9zr3lHfe72mJdMxaI5cNzG6hSvl2OtQ6s66jxdsn675s03yd/xuozveWOty9o1Vkyr9rvIFNfxoda1IX8nGmNcTPV3vFX3O+T1bKgczzPknD6o78P8nlettcPGNbpz5849NMnW53/+5+fTP/3Txw5nsi5cuJDf/u3fzpd8yZfk9OnTY4czWfK0HHlajjwt7yMf+Uje/va3J8nmTTfd9J6x45mqc+fOfXGSN48dBwCdPeqmm256y9hBTJlrHsCx4Zo3h+sdwLGy8JrnTjgAmLa3JnlUkjuTfGLkWAA4vKuSfFZ213Pmc80DWG+uectxvQNYf0tf8xThAGDCbrrppgtJ/C9SgPX2rrEDWAeueQDHgmveAq53AMfGUte8+wwdBQAAAAAAAJw0inAAAAAAAADQM0U4AAAAAAAA6JkiHAAAAAAAAPRMEQ4AAAAAAAB6pggHAAAAAAAAPVOEAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAerZSEa6qbqiqX62qe6rq/VX10qq6eol2VVXPrqp3V9XHquo3q+rRq8QAAAAAAAAAU3XoIlxVXZfkDUmuTvLkJM9N8q1JXrZE8+9J8n1JfjjJ1yV5X5LXVdWfOmwcAAAAAAAAMFWnVmjz9CTXJnlSa+0DSVJVp5K8oqpe1Fp775UaVdVGkuck+aHW2g/Ptv27JG9P8qwkz1ghFgAAAAAAAJicVT6O8glJXr9XgJv5hVlfj5/T7kuzW7z7hb0NrbWLSX4pydesEAcAAAAAAABM0ipFuBuS3H7phtba3dn9aMkbFrTL/rZJ3pbkc6vq01aIBQAAAAAAACZnlY+jvC7J3VfYfleSBy5od6G1tnOFdjV7/mNXalhV12b3Lrokybd927c95ClPeUo++tGPHiLsk+fixYupqtxzzz259957xw5nsuRpOfK0HHla3iVr+FVjxgEAAAAAMIRVinBjeGaS5+89OHfuXJ7ylKfkPe95z4ghrYdrrrkmd9xxx9hhTJ48LUeeliNPh/POd77z85K8a+w4AAAAAAD6tEoR7q4k97/C9uuSfOAK2y9td7qqNvbdDXddkjZ7/iAvS/KqvQcf/vCHP/cZz3jG+cc97nGPfvKTn/yflw/9ZPm5n/u562+55ZY3P/GJT3zUzTff/P6x45kqeVqOPC1Hnpb3S7/0Sw95/etf/1unTp3a+uZv/uaxwwEAAAAA6FW11g7XoOrfJvn/WmtPumTb/bNbRPtbrbVXH9DuLyc5l+SRrbX/cMn2H0rylNba5x0ihocm2Uqy2VpzO9wB5Gk58rQceVqOPC1PrgAAAACA4+w+K7T5lSSPq6oHXLLtqUk+meR1c9r9RpIPzV6bJKmq+yZ5cpLXrhAHAAAAAAAATNIqRbhXJvlwktdU1eOr6puS/GCSV7bW3rv3oqo6V1V/sPd49hGUP5DkWVX1d2Z3xv1cks9M8o+6HAQAAAAAAABMyaG/E661dldV3ZTkR5O8JrsFuVcled6+l151hf5fkqSSPCvJZyW5NclXtdb+0yHD+FCS75v95GDytBx5Wo48LUeelidXAAAAAMCxdejvhAMAAAAAAADmW+XjKAEAAAAAAIA5FOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOjZ5IpwVXVDVf1qVd1TVe+vqpdW1dVLtKuqenZVvbuqPlZVv1lVjz6KmMewSp6q6sGz191aVR+uqvdU1c9W1cOOKu6jtup42tfHd1VVq6p/NVScY+uSp6p6SFX9VFXdOZt7b6uqvzZ0zGPosD59ZlW9crY+3VNVv1dVTz+KmMdQVY+YHe+tVfXxqvq9JdudqHUcAAAAADjeTo0dwKWq6rokb0jyjiRPTvKQJC9Lcr8k37Gg+fck+b4kz07yu0m+PcnrquqRrbX/NFjQI+iQpy+avf4nk/xWkjNJvjfJv6+qL2yt3Tlk3Eet43ja6+P6JM9P8kcDhTm6Lnmqqgcn+c0kv5/kW5N8KMmNSU4PGPIoOo6nf5bkhiTPTfLuJF+T5Mer6hOttZ8YLOjx3Jjka5P8dnb/s8ey/+HjxKzjAAAAAMDxV621sWP4Y1X1nCTPS/K5rbUPzLZ9a5JXzLa994B2G0n+MMmPtdaeO9t2dZK3J3lta+0ZRxH/UemQpwck+Uhr7eOXbHtodosCf6+19kNDx36UVs3Tvj7+aZKW5GHZzd3XDRjyKLrkqap+OsmfSvKY1tonjiLesXSYd9cneV+Sb2qtvfqS7W9K8vHW2k1Dx37Uquo+rbVPzv7+6iRf3Fr7wgVtTtQ6DgAAAAAcf1P7OMonJHn93j9wz/xCduN8/Jx2X5rk2tlrkySttYtJfim7d5wcNyvlqbV296UFuNm29yS5M8nnDBHoyFYdT0mSqvqyJH8lu3flHGcr5amqrk3y9UlecdwLcDOrjqf7zn5+cN/2Dyap/sKbjr0C3CGdtHUcAAAAADjmplaEuyHJ7ZduaK3dnd27SG5Y0C772yZ5W5LPrapP6yvAiVg1T5epqs9P8qDs5uq4WTlPVXVVkpcn+f7W2vuGCnAiVs3Tn09ydZJ7q+pNVXXv7HvSXlJV953Tbl2tlKfW2laS1yV5blV9QVV9RlV9fXYLdz82XLhr56St4wAAAADAMTe1Itx1Se6+wva7kjxwQbsLrbWdK7Sr2fPHyap5+hOqqpL84yTvTfJzvUQ2LV3y9Iwk1yT54Z5jmqJV83T97Oerkrwlu0WlH07yXUle2F94k9FlPD05ux+1eFt2vzfvZ5N8d2vtn/cZ4Jo7aes4AAAAAHDMnRo7AEb1giQ3Jfnq1to9I8cyGVX1oOwWkf7G7OPwuLK9Iv7rW2t/d/b3X6uqz0jyrKp6YWvtYyPFNhmzYvf/keS/SPIN2b1z7iuT/EhV3dVaOztmfAAAAAAADGNqRbi7ktz/CtuvS/KBK2y/tN3pqtrYdxfFdUna7PnjZNU8/bGq+pYkfz/J01pr53qMbUpWzdMLk/xukn9XVQ+YbTuV5NTs8Uf2f7femusy75LkDfu2n0vyvCSPSPLWztFNx6p5+tokT03y51pre/l446zY+0NJFOF2nbR1HAAAAAA45qb2cZS3Z993K1XV/ZM8OJd/T9D+dknyX+7bfkOSdx/Du3FWzdPea5+U5MeT/P3W2k8OEuE0rJqnG5I8Jrv/6L/35y8l+arZ3x83RLAjWjVP/3FBvxsd45qaVfP0BUk+keT39m3/nSSfU1X36zPINXbS1nEAAAAA4JibWhHuV5I87pK7j5LdO0g+meR1c9r9Rna/Z+mpexuq6r7Z/R6m1/Yf5uhWzVOq6rHZ/f63n2it/YOB4puKVfP0XUm+Yt+f/5Dkt2Z///cDxDqmlfLUWntXdu9021+U/MokH8viIt26WXU8vSvJVUn+3L7tX5Tkj1prH+0zyDV20tZxAAAAAOCYq9ba2DH8saq6LsltSd6e5EVJHpLkZUl+prX2HZe87lySh7XWHnHJtmdn9zvOvie7hYFnJHl8kke21v7TUR3DUVg1T1X1Z5L8ZpKtJP9TdosHe+5srf2/R3MER6PLeLpCX2/M7sdQft2gQY+g47z7b5PckuQfJ/nlJI/K7jx8aWvtfzmqYzgKHebdZ2T3LriLSb4vu98J9/gkz0ry/NbaPzzK4zgKs7v7vmb28NuT/Okkz5w9flNr7c6Tvo4DAAAAAMffpL4TrrV2V1XdlORHk7wmyYeTvCq73y91qatyeewvSVLZ/Yftz0pya5KvOo7/cNshT1+S3e+0un+S8/te+1NJvnGAcEfTcTydGF3y1Fr7l1V1c5LvTfJt2S0wPT/JiwcO+8itmqfW2odn7b4/u+vUA5K8M7tFqZcPHvg4HpTkn+3btvf4K5K8MSd8HQcAAAAAjr9J3QkHAAAAAAAAx8HUvhMOAAAAAAAA1p4iHAAAAAAAAPRMEQ4AAAAAAAB6pggHAAAAAAAAPVOEAwAAAAAAgJ4pwgEAAAAAAEDPFOEAAAAAAACgZ4pwAAAAAAAA0DNFOAAAAAAAAOiZIhwAAAAAAAD0TBEOAAAAAAAAeqYIBwAAAAAAAD1ThAMAAAAAAICe/f8QZXXNrVdSbAAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 2200x550 with 4 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"@jit(nopython=True, nogil=True)\n",
|
|
"def energy(state):\n",
|
|
" E = 0\n",
|
|
" N, M = state.shape\n",
|
|
" for i in range(N):\n",
|
|
" for j in range(M):\n",
|
|
" # handle the north and south neighbours\n",
|
|
" if 0 <= (i + 1) < N:\n",
|
|
" E -= state[i, j] * state[i + 1, j]\n",
|
|
"\n",
|
|
" # handle the east and west neighbours\n",
|
|
" if 0 <= (j + 1) < M:\n",
|
|
" E -= state[i, j] * state[i, j + 1]\n",
|
|
"\n",
|
|
" return 2 * E / (N * M)\n",
|
|
"\n",
|
|
"\n",
|
|
"# While writing numba it's useful to keep the list of supported numpy functions open:\n",
|
|
"# https://numba.pydata.org/numba-doc/dev/reference/numpysupported.html\n",
|
|
"@jit(nopython=True, nogil=True)\n",
|
|
"def mcmc(initial_state, steps, T, energy=energy):\n",
|
|
" N, M = initial_state.shape\n",
|
|
" assert N == M\n",
|
|
"\n",
|
|
" current_state = initial_state.copy()\n",
|
|
" E = N**2 * energy(current_state)\n",
|
|
" for i in range(steps):\n",
|
|
" i, j = np.random.randint(N), np.random.randint(N)\n",
|
|
"\n",
|
|
" # modify the state a little, here we just flip a random pixel\n",
|
|
" current_state[i, j] *= -1\n",
|
|
" new_E = N**2 * energy(current_state)\n",
|
|
"\n",
|
|
" if (new_E < E) or np.exp(-(new_E - E) / T) > np.random.random():\n",
|
|
" E = new_E\n",
|
|
" else:\n",
|
|
" current_state[i, j] *= -1 # reject the change we made\n",
|
|
"\n",
|
|
" return current_state\n",
|
|
"\n",
|
|
"\n",
|
|
"Ts = [4, 5, 50]\n",
|
|
"\n",
|
|
"ncols = 1 + len(Ts)\n",
|
|
"f, axes = plt.subplots(ncols=ncols, figsize=(5 * ncols, 5))\n",
|
|
"\n",
|
|
"initial_state = np.ones(shape=(50, 50))\n",
|
|
"axes[0].set(title=\"Initial state\")\n",
|
|
"show_state(initial_state, ax=ax)\n",
|
|
"\n",
|
|
"for T, ax in zip(Ts, axes[1:]):\n",
|
|
" # initial_state = rng.choice([1,-1], size = (50,50))\n",
|
|
"\n",
|
|
" final_state = mcmc(initial_state, steps=100_000, T=T)\n",
|
|
" show_state(final_state, ax=ax)\n",
|
|
" ax.set(title=f\"T = {T}\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5d1874d4-4585-49ed-bc6f-b11c22231669",
|
|
"metadata": {},
|
|
"source": [
|
|
"These images give a flavour of why physicists find this model useful, it gives window into how thermal noise and spontaneous order interact. At low temperatures the energy cost of being different from your neighbours is the most important thing, while at high temperatures, it doesn't matter and you really just do your own thing.\n",
|
|
"\n",
|
|
"There's a special point somewhere in the middle called the critical point $T_c$ where all sorts of cool things happen, but my favourite is that for large system sizes you get a kind of fractal behaviour which I will demonstrate more once we've sped this code up and can simulate larger systems in a reasonable time. You can kinda see it for 50x50 systesm at T = 5 but not really clearly."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "5f728039-a975-4083-b68e-a13b4f2d1f87",
|
|
"metadata": {},
|
|
"source": [
|
|
"The code we have so far is really just a sketch of a solution. So this is a good time to step back and think about what are aims are and how this software will fulfil them. I see three broad areas on which it needs improvement:\n",
|
|
"\n",
|
|
"**Functionality**\n",
|
|
"Right now we can't really do much except print nice pictures of states, but (within the fiction of this project) we really want to be able to do science! So we need to think about what measurements and observations we might want to make and how that might affect the structure of our code.\n",
|
|
"\n",
|
|
"**Testing**\n",
|
|
"I've already missed at least one devastating bug in this code, and there are almost certainly more! Before we start adding too much new code we should think about how to increase our confidence that the individual components are working correctly. It's very easy to build a huge project out of hundreds of functions, realise there's a bug and then struggle to find the source of that bug. If we test our components individually and thoroughly, we can avoid some of that pain.\n",
|
|
"\n",
|
|
"**Performance**\n",
|
|
"Performance only matters in so far as it limits what we can do. And there is a real danger that trying to optimise for performance too early or in the wrong places will just lead to complexity that makes the code harder to read, harder to write and more likely to contain bugs. However I do want to show you the fractal states at the critical point, and I can't currently generate those images in a reasonable time, so some optimisation will happen!"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "486f066c-f027-44e8-8937-8636a52f32fb",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Functionality\n",
|
|
"\n",
|
|
"The main thing we want to be able to do is to take measurements, the code as I have writting it doesn't really allow that because it only returns the final state in the chain. Let's say we have a measurement called `average_color(state)` that we want to average over the whole chain. We could just stick that inside our definition of `mcmc` but we know that we will likely make other measurements too and we don't want to keep writing new versions of our core functionality!\n",
|
|
"\n",
|
|
"## Exercise 1\n",
|
|
"Have a think about how you would implement this and what options you have."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "c28b0a86-28f8-426f-9013-70e962f02256",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Solution 1\n",
|
|
"So I chatted with my mentors on this project on how to best do this and we came up with a few ideas:\n",
|
|
"\n",
|
|
"### Just save all the states and return them\n",
|
|
"\n",
|
|
"The problem with this is the states are very big and we don't want to waste all that memory. For an NxN state that uses 8 bit integers (the smallest we can use in numpy) 1000 samples would already use 2.5Gb of memory! We will see later that we'd really like to be able to go a bit bigger than 50x50 and 1000 samples!\n",
|
|
"\n",
|
|
"### Pass in a function to make measurements\n",
|
|
"```python\n",
|
|
"\n",
|
|
"def mcmc(initial_state, steps, T, measurement, energy=energy):\n",
|
|
" ...\n",
|
|
"\n",
|
|
" current_state = initial_state.copy()\n",
|
|
" E = N**2 * energy(current_state)\n",
|
|
" for i in range(steps):\n",
|
|
" measurements[i] = measurement(state)\n",
|
|
" ...\n",
|
|
"\n",
|
|
" return measurements\n",
|
|
"```\n",
|
|
"\n",
|
|
"This could work but it limits how we can store measurements and what shape and type they can be. What if we want to store our measurements in a numpy array? Or what if your measurement itself is a vector or and object that can't easily be stored in a numpy array? We would have to think carefully about what functionality we want."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "c7c9575f-2450-4298-a507-90f0c1b9b284",
|
|
"metadata": {
|
|
"tags": []
|
|
},
|
|
"source": [
|
|
"### Use Inheritance\n",
|
|
"```python\n",
|
|
"# This class would define the basic functionality of performing MCMC\n",
|
|
"class MCMCSampler(object):\n",
|
|
" def run(self, initial_state, steps, T):\n",
|
|
" ...\n",
|
|
" for i in range(steps):\n",
|
|
" self.measurement(state)\n",
|
|
"\n",
|
|
" \n",
|
|
"# This class would inherit from it and just implement the measurement\n",
|
|
"class AverageColorSampler(MCMCSampler):\n",
|
|
" measurements = np.zeros(10)\n",
|
|
" index = 0\n",
|
|
" \n",
|
|
" def measurement(self, state):\n",
|
|
" self.measurements[self.index] = some_function(state)\n",
|
|
" self.index += 1\n",
|
|
" \n",
|
|
"color_sampler = AverageColorSampler(...)\n",
|
|
"measurements = color_sampler.run(...)\n",
|
|
"```\n",
|
|
"\n",
|
|
"This would definitely work but I personally am not a huge fan of object oriented programming so I'm gonna skip this option!"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7d05d25d-c9ba-406d-9977-0ca4aeb430a7",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Use a generator\n",
|
|
"This is the approach I ended up settling on, we will use [python generator function](https://peps.python.org/pep-0255/) "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 44,
|
|
"id": "f73d6335-6514-45b1-9128-d72122d8b0b7",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"array([[[-9.86830992e+148, -9.86830992e+148, -9.86830992e+148, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" ...,\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000]],\n",
|
|
"\n",
|
|
" [[-9.86830992e+148, -9.86830992e+148, -9.86830992e+148, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" ...,\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000]],\n",
|
|
"\n",
|
|
" [[-9.86830992e+148, -9.86830992e+148, -9.86830992e+148, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" ...,\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000]],\n",
|
|
"\n",
|
|
" ...,\n",
|
|
"\n",
|
|
" [[-9.86830992e+148, -9.86830992e+148, -9.86830992e+148, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" ...,\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000]],\n",
|
|
"\n",
|
|
" [[-9.86830992e+148, -9.86830992e+148, -9.86830992e+148, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" ...,\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000]],\n",
|
|
"\n",
|
|
" [[-9.86830992e+148, -9.86830992e+148, -9.86830992e+148, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" ...,\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000],\n",
|
|
" [ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000, ...,\n",
|
|
" 1.00000000e+000, 1.00000000e+000, 1.00000000e+000]]])"
|
|
]
|
|
},
|
|
"execution_count": 44,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"@jit(nopython=True, nogil=True)\n",
|
|
"def mcmc(initial_state, steps, T, energy=energy):\n",
|
|
" N, M = initial_state.shape\n",
|
|
" assert N == M\n",
|
|
"\n",
|
|
" current_state = initial_state.copy()\n",
|
|
" E = N**2 * energy(current_state)\n",
|
|
" for i in range(steps):\n",
|
|
" i, j = np.random.randint(N), np.random.randint(N)\n",
|
|
"\n",
|
|
" # modify the state a little, here we just flip a random pixel\n",
|
|
" current_state[i, j] *= -1\n",
|
|
" new_E = N**2 * energy(current_state)\n",
|
|
"\n",
|
|
" if (new_E < E) or np.exp(-(new_E - E) / T) > np.random.random():\n",
|
|
" E = new_E\n",
|
|
" else:\n",
|
|
" current_state[i, j] *= -1 # reject the change we made\n",
|
|
"\n",
|
|
" yield current_state # give the state out to the enclosing function but don't actually return\n",
|
|
"\n",
|
|
" return # this signals that we're done\n",
|
|
"\n",
|
|
"\n",
|
|
"initial_state = np.ones(shape=(50, 50))\n",
|
|
"np.array([s for s in mcmc(initial_state, steps=10, T=5)])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"id": "193b778f-5913-48f1-9df6-304ab50ceb4e",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python [conda env:recode]",
|
|
"language": "python",
|
|
"name": "conda-env-recode-py"
|
|
},
|
|
"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
|
|
}
|