{
 "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
}