diff --git a/Exercise sheet 7/exercise_sheet_07.ipynb b/Exercise sheet 7/exercise_sheet_07.ipynb new file mode 100644 index 0000000..15caf5d --- /dev/null +++ b/Exercise sheet 7/exercise_sheet_07.ipynb @@ -0,0 +1,599 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e301721f", + "metadata": {}, + "source": [ + "# Exercise sheet\n", + "\n", + "Some general remarks about the exercises:\n", + "* For your convenience functions from the lecture are included below. Feel free to reuse them without copying to the exercise solution box.\n", + "* For each part of the exercise a solution box has been added, but you may insert additional boxes. Do not hesitate to add Markdown boxes for textual or LaTeX answers (via `Cell > Cell Type > Markdown`). But make sure to replace any part that says `YOUR CODE HERE` or `YOUR ANSWER HERE` and remove the `raise NotImplementedError()`.\n", + "* Please make your code readable by humans (and not just by the Python interpreter): choose informative function and variable names and use consistent formatting. Feel free to check the [PEP 8 Style Guide for Python](https://www.python.org/dev/peps/pep-0008/) for the widely adopted coding conventions or [this guide for explanation](https://realpython.com/python-pep8/).\n", + "* Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting `Kernel > Restart & Run All` in the jupyter menu.\n", + "* For some exercises test cases have been provided in a separate cell in the form of `assert` statements. When run, a successful test will give no output, whereas a failed test will display an error message.\n", + "* Each sheet has 100 points worth of exercises. Note that only the grades of sheets number 2, 4, 6, 8 count towards the course examination. Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are just for feedback.\n", + "\n", + "Please fill in your name here:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ce67703", + "metadata": {}, + "outputs": [], + "source": [ + "NAME = \"\"\n", + "NAMES_OF_COLLABORATORS = \"\"" + ] + }, + { + "cell_type": "markdown", + "id": "d32aa225", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "ca31ce2e", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "2139ee04c015c84926865a605fd18730", + "grade": false, + "grade_id": "cell-58dadbd252381fce", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "**Exercise sheet 7**\n", + "\n", + "Code from the lectures:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d10dcc0d", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "b22c1a51a55c2aa684aa56d0a75c8141", + "grade": false, + "grade_id": "cell-d0c2906276dc11be", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "rng = np.random.default_rng() \n", + "import matplotlib.pylab as plt\n", + "%matplotlib inline\n", + "\n", + "def potential_v(x,lamb):\n", + " '''Compute the potential function V(x).'''\n", + " return lamb*(x*x-1)*(x*x-1)+x*x\n", + "\n", + "def neighbor_sum(phi,s):\n", + " '''Compute the sum of the state phi on all 8 neighbors of the site s.'''\n", + " w = len(phi)\n", + " return (phi[(s[0]+1)%w,s[1],s[2],s[3]] + phi[(s[0]-1)%w,s[1],s[2],s[3]] +\n", + " phi[s[0],(s[1]+1)%w,s[2],s[3]] + phi[s[0],(s[1]-1)%w,s[2],s[3]] +\n", + " phi[s[0],s[1],(s[2]+1)%w,s[3]] + phi[s[0],s[1],(s[2]-1)%w,s[3]] +\n", + " phi[s[0],s[1],s[2],(s[3]+1)%w] + phi[s[0],s[1],s[2],(s[3]-1)%w] )\n", + "\n", + "def scalar_action_diff(phi,site,newphi,lamb,kappa):\n", + " '''Compute the change in the action when phi[site] is changed to newphi.'''\n", + " return (2 * kappa * neighbor_sum(phi,site) * (phi[site] - newphi) +\n", + " potential_v(newphi,lamb) - potential_v(phi[site],lamb) )\n", + "\n", + "def scalar_MH_step(phi,lamb,kappa,delta):\n", + " '''Perform Metropolis-Hastings update on state phi with range delta.'''\n", + " site = tuple(rng.integers(0,len(phi),4))\n", + " newphi = phi[site] + rng.uniform(-delta,delta)\n", + " deltaS = scalar_action_diff(phi,site,newphi,lamb,kappa)\n", + " if deltaS < 0 or rng.uniform() < np.exp(-deltaS):\n", + " phi[site] = newphi\n", + " return True\n", + " return False\n", + "\n", + "def run_scalar_MH(phi,lamb,kappa,delta,n):\n", + " '''Perform n Metropolis-Hastings updates on state phi and return number of accepted transtions.'''\n", + " total_accept = 0\n", + " for _ in range(n):\n", + " total_accept += scalar_MH_step(phi,lamb,kappa,delta)\n", + " return total_accept\n", + "\n", + "def batch_estimate(data,observable,k):\n", + " '''Devide data into k batches and apply the function observable to each.\n", + " Returns the mean and standard error.'''\n", + " batches = np.reshape(data,(k,-1))\n", + " values = np.apply_along_axis(observable, 1, batches)\n", + " return np.mean(values), np.std(values)/np.sqrt(k-1)" + ] + }, + { + "cell_type": "markdown", + "id": "d0a863dd", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "79366cc7caf65397bcf57caebfe394c5", + "grade": false, + "grade_id": "cell-5dae26c757131cdf", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "## Lattice scalar field & heatbath algorithm\n", + "\n", + "**(100 points)**\n", + "\n", + "**Goal**: Implement and test the heatbath algorithm for the lattice scalar field simulation.\n", + "\n", + "As in the lecture we consider the scalar field $\\varphi : \\mathbb{Z}_w^4 \\to \\mathbb{R}$ on the $w^4$ lattice with periodic boundary conditions (here we use the notation $\\mathbb{Z}_w = \\{0,1,\\ldots,w-1\\}$ for the integers modulo $w$) and lattice action \n", + "\n", + "$$\n", + "S_L[\\varphi] = \\sum_{x\\in\\Lambda}\\left[ V(\\varphi(x)) - 2\\kappa \\sum_{\\hat{\\mu}}\\varphi(x)\\varphi(x+\\hat{\\mu})\\right],\\qquad V(\\varphi) = \\lambda (\\varphi^2-1)^2+\\varphi^2,\n", + "$$\n", + "\n", + "where the second sum is over the 4 unit vectors $\\hat{\\mu} = (1,0,0,0), (0,1,0,0), (0,0,1,0), (0,0,0,1)$.\n", + "\n", + "In the lecture notes the Metropolis-Hastings algorithm for local updates to $\\phi$ has been implemented to sample from the desired distribution $\\pi(\\varphi) = \\frac{1}{Z} e^{-S_L[\\varphi]}$. As in the case of the Ising model this local updating suffers from **critical slowing down** when close to the symmetry-breaking transition. So one should always investigate alternative Monte Carlo updates. One of these is the **heatbath algorithm** that we have already discussed in the setting of the harmonic oscillator. It entails selecting a uniform site $x_0\\in \\mathbb{Z}_w^4$ and sampling a new random value $Y$ for $\\varphi(x_0)$ with density proportional\n", + "to $f_Y(y)\\propto\\pi(\\varphi(x)|_{\\varphi(x_0)=y})$ on $\\mathbb{R}$. To sample such a $y$ we will make use of **acceptance-rejection sampling** (recall the lecture of week 3)." + ] + }, + { + "cell_type": "markdown", + "id": "96f10149", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "1f6b27e4b5c5e89df7aa07edc57f6471", + "grade": false, + "grade_id": "cell-4ac4b782063a8df9", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(a)__ Show with a calculation, to be included below in markdown and LaTeX, that for any choice of $c>0$ we have\n", + "\n", + "$$\n", + "f_Y(y)\\propto \\exp\\left( - c (y - s/c)^2 - \\lambda (y^2 - v)^2 \\right),\n", + "$$\n", + "\n", + "where $s$ and $v$ are given by\n", + "\n", + "$$\n", + "s := \\kappa \\sum_{\\hat{\\mu}} (\\varphi(x_0+\\hat{\\mu}) + \\varphi(x_0-\\hat{\\mu})), \\qquad v := 1 + \\frac{c-1}{2\\lambda}.\n", + "$$\n", + "\n", + "**(15 pts)**" + ] + }, + { + "cell_type": "markdown", + "id": "41c68351", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "e4ef65e5cbb01bcc92764b99989dff55", + "grade": true, + "grade_id": "cell-f61eb44774462387", + "locked": false, + "points": 15, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "source": [ + "YOUR ANSWER HERE" + ] + }, + { + "cell_type": "markdown", + "id": "26ea9ea4", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "14b441ebb964df8db04811878804faf4", + "grade": false, + "grade_id": "cell-182ca6c04fd76105", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(b)__ Implement acceptance/rejection sampling of the random variable $Y$ with pdf $f_Y(y)\\propto\\exp\\left( - c (y - s/c)^2 - \\lambda (y^2 - v)^2 \\right)$. You may treat $c, s, \\lambda$ as given parameters for this. Test your sampling by producing a histogram for $\\lambda = 3$, $s = 0.3$ and three different values for $c$ (e.g. 0.5, 1.0, 2.0). _Hint:_ use proposal distribution $f(y) = \\sqrt{\\frac{c}{\\pi}}\\exp\\left( - c (y - s/c)^2\\right)$ and appropriate acceptance probability $A(y)$. **(25 pts)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1aff44bf", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "243cd5ce582605801df65116cb68a5c7", + "grade": false, + "grade_id": "cell-a7996b41ba76514e", + "locked": false, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def sample_Y(s,lamb,c):\n", + " '''Sample Y width density f_Y(y).'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02a240a7", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "c1f71a59b63866ec594ac54065efe09b", + "grade": true, + "grade_id": "cell-183431967b5b0320", + "locked": true, + "points": 15, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "from nose.tools import assert_almost_equal\n", + "assert_almost_equal(np.mean([sample_Y(0.8,1.5,2.0) for _ in range(800)]),0.71,delta=0.06)\n", + "assert_almost_equal(np.mean([sample_Y(-0.7,0.5,0.6) for _ in range(800)]),-0.60,delta=0.06)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b948366e", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "831381086e376d13aaae36346439fb1d", + "grade": true, + "grade_id": "cell-35ac6348edb77983", + "locked": false, + "points": 10, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Plotting\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "id": "068b2d7d", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "e5fb80691276aacf937880fd6f3ca538", + "grade": false, + "grade_id": "cell-757abc5290a5d8cd", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(c)__ We still have a free parameter $c>0$ that we can choose as we see fit. We would like to choose $c$ such that the acceptance probability $\\mathbb{P}(\\text{accept})=\\int_{-\\infty}^{\\infty} \\mathrm{d}y f(y) A(y)$ is maximized. Show by a calculation, to be included in markdown and LaTeX below, that $c$ must then be a positive root of the polynomial\n", + "\n", + "$$\n", + "-c^3 + (1 - 2 \\lambda)c^2 + \\lambda c + 2 s^2 \\lambda = 0.\n", + "$$\n", + "\n", + "_Hint:_ evaluate $\\frac{d}{d c}f(y) A(y)$ and set it equal to $0$. **(10 pts)**" + ] + }, + { + "cell_type": "markdown", + "id": "abfbbc40", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "b97bd524247c5525a006a1a4ec006fc9", + "grade": true, + "grade_id": "cell-8fec7758f04ca5ae", + "locked": false, + "points": 10, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "source": [ + "YOUR ANSWER HERE" + ] + }, + { + "cell_type": "markdown", + "id": "95847e73", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "b090f3a14794e0f04ac15a701d8777bd", + "grade": false, + "grade_id": "cell-e3e8d6226e51229b", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "In fact, it can be easily shown that the polynomial above has only a single positive root when $\\lambda>0$, so this uniquely determines an optimal $c$. Computing it accurately can be costly, but it is possible to obtain a good approximation based on a Taylor series estimate, which yields:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76c7f67e", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "0a559e8802ae7f9ec73e690e845ff95d", + "grade": false, + "grade_id": "cell-e9fe1618aa90e9e6", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "outputs": [], + "source": [ + "def approx_optimal_c(s,lamb):\n", + " '''Compute a value of c that is close to a positive root of the polynomial above.'''\n", + " u = np.sqrt(1+4*lamb*lamb)\n", + " return ((3 + 3*u*(1-2*lamb)+4*lamb*(3*lamb-1) \n", + " + np.sqrt(16*s*s*(1+3*u-2*lamb)*lamb +(1+u-2*u*lamb+4*lamb*lamb)**2)) /\n", + " (2+6*u-4*lamb))\n", + "\n", + "def sample_Y_optimal(s,lamb):\n", + " c = approx_optimal_c(s,lamb)\n", + " return sample_Y(s,lamb,c)\n", + " \n", + "approx_optimal_c(0.3,3)" + ] + }, + { + "cell_type": "markdown", + "id": "b1882f8b", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "5fc984f5a02b48d3599cf61bbd09167d", + "grade": false, + "grade_id": "cell-3562ebb7562ac692", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(d)__ Implement the MCMC update using the heatbath with (approximately) optimal $c$. **(10 pts)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08b69f13", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "af478c34871b3055d1d61c8a652f9d32", + "grade": true, + "grade_id": "cell-2045e775f42306f1", + "locked": false, + "points": 10, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def heatbath_update(phi,lamb,kappa):\n", + " '''Perform a random heatbath update on the state phi (no return value necessary).'''\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "def run_scalar_heatbath(phi,lamb,kappa,n):\n", + " '''Perform n heatbath updates on state phi.'''\n", + " for _ in range(n):\n", + " heatbath_update(phi,lamb,kappa)" + ] + }, + { + "cell_type": "markdown", + "id": "465f6b26", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "748558ea83d69a646873ac1cc9b3682c", + "grade": false, + "grade_id": "cell-8c54007d89b7ec4b", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "__(e)__ Gather $800$ samples of $|m| = |w^{-4} \\sum_{x\\in\\mathbb{Z}_w^4} \\varphi(x)|$ for $w=3$, $\\lambda = 1.5$ and $\\kappa = 0.08, 0.09, \\ldots, 0.18$ for the heatbath algorithm as well as for the Metropolis-Hastings algorithm (width $\\delta=1.5$) from the lecture. Store your data in the HDF5-file `scalarfield.hdf5`. You may use $500$ sweeps for equilibration $1$ sweep in between measurements. **(15 pts)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7013334", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "d09ba2c288f6d97cdb635656d064a586", + "grade": true, + "grade_id": "cell-3ed06bbba8ae4e6e", + "locked": false, + "points": 15, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Gathering and storing data\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "id": "e3df064d", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "b3cc3deba1dc548f2ce06ad7871ef11b", + "grade": false, + "grade_id": "cell-195d69ed7b2975ce", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "**(f)** Read your data from `scalarfield.hdf5` and :\n", + "* plot your estimates for $\\mathbb{E}[|m|]$ with error bars obtained from the heatbath and Metropolis algorithm for comparison;\n", + "* plot only the errors in your estimate so you can compare their magnitude for both algorithms;\n", + "* plot the estimated autocorrelation time in units of sweeps for both algorithms. \n", + "\n", + "Where do you see an improvement? **(25 pts)**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "845380df", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "0bfbc9b1733ef2ebe02375d39c880a8e", + "grade": true, + "grade_id": "cell-064aaf697d09bad1", + "locked": false, + "points": 25, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Reading data and plotting\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c977d3ba", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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 +}