hamiltonian-mechanics-torch.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://colindcarroll.com/2019/04/11/hamiltonian-monte-carlo-from-scratch/\n",
    "# hamiltonian monte carlo, implementation 1:\n",
    "import torch\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 466,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://www.mcmchandbook.net/HandbookChapter5.pdf | 5.2.3.3 The Leapfrog Method\n",
    "def leapfrog(q, p, Vdq, eps):\n",
    "    p -= 0.5 * eps * Vdq(q)\n",
    "    \n",
    "    q.detach_()\n",
    "    q += eps * p\n",
    "    \n",
    "    p -= 0.5 * eps * Vdq(q)\n",
    "    return [q, p]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 467,
   "metadata": {},
   "outputs": [],
   "source": [
    "# kinetic energy\n",
    "def K(pi, p, q):\n",
    "    return 0.5 * (p.transpose(-1, 0)  @ p)\n",
    "\n",
    "# potential energy at point p for distribution pi\n",
    "def V(pi, q):\n",
    "    return pi(q)\n",
    "\n",
    "# closure\n",
    "def mkPiVdq(pi):\n",
    "    def vdq(q):\n",
    "        q.detach_()\n",
    "        q.requires_grad_()\n",
    "        vcur = V(pi, q)\n",
    "        vcur.backward()\n",
    "        return q.grad\n",
    "    return vdq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 470,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "acceptance ratio: 0.65\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fda60513668>]"
      ]
     },
     "execution_count": 470,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt0nGd9J/Dvb2Y0uszoZs1YlizJkmzZWMSOnWiTACWkQNjcTYEucYEDuzmkWQgLB9jdUNq0S/9oKee0pSUUHJayBEg2wAJuMAQaaJImmMT3xHZs3WzrZmt0l0a3ufz2j5nXFopljaV35r3M93NOEs3o0by/ieWvHj3P8z6PqCqIiMhdPFYXQERE5mO4ExG5EMOdiMiFGO5ERC7EcCciciGGOxGRCzHciYhciOFORORCDHciIhfyWXXhUCikjY2NVl2eiMiRDh48OKSq4eXaWRbujY2NOHDggFWXJyJyJBE5m0k7DssQEbkQw52IyIUY7kRELsRwJyJyIYY7EZELMdyJiFyI4U5E5ELLhruIfFNEBkXk1SU+LyLyDyLSISLHROQ688skIqKrkclNTN8C8BUA317i87cDaEn/cyOAf0r/l4gcRv6XvO45/XOes+xEy/bcVfU5ACNXaLILwLc1ZT+AChGpMatAIsqNywX7lZ4nezNjzH09gJ4Fj3vTzxE5UjKp+KufncRfPnUCiSR7reRMOd1bRkTuB3A/ADQ0NOTy0kQZSSYVn//xK3j8pVR/ZXwmhr9573Z4POy9krOYEe59AOoXPK5LP/c6qroHwB4AaGtrY5eIbEVV8Wc/eRWPv9SDj//+Rvg8Hnz5mXZ4RfBX79nGgCdHMSPc9wJ4UESeQGoidVxVB0x4XaKc+uLPT+G7vz2HB962EZ991xYAQFIV//irDgSLfPizu1otrpAoc5kshXwcwG8AbBGRXhG5T0QeEJEH0k32AegC0AHgUQAfy1q1RFkyNRfHP7/QjV07avE/b9sCEYGI4NO3bsb7rq/DY785i/GZmNVlZtVSq2K4WsaZlu25q+ruZT6vAD5uWkVEFvjXExcwF0/iAzdugMil4RcRwQdubMAPDvbiF8fP4w/b6q/wKs7HIHcP3qFKBOCpY/1YV1aEtg2Vr/vcjvoK1FUW46ljHG0k52C4U94bn4nh2dMR3LGt5rKTpiKCO7fX4IWOIYxG5y2okOjqMdwp7/3i+HnEEoq7r1363ru7t9cinlT8/Pj5HFZGtHIMd8p7Tx0bQF1lMXbUVyzZ5o21ZWisKsFTx/pzWBnRyjHcKa+NRufxQscQ7txe8zsTqYuJCO7aXovfdA4jMjmXwwqJVobhTnnt58fPI55U3L29dtm2d11bg6QCP3+VE6tkfwx3yms/PTaAxqoSvLG2bNm2W6pLsWltEP/CVTPkAAx3yltz8QRe6h7Bra3VVxySMYgIbm2txqGzo5iZT+SgQqKVY7hT3jrRP4H5RBLXNbx+bftSrm+oRDypeLV/PIuVEa0ew53y1uFzYwCAnVcR7jsaKtJfO5qVmojMwnCnvHWkZww15UVYV16U8deEgoWoX1OMIz1jWayMaPUY7pS3DveMXnFt+1J21lde7PUT2RXDnfLS0NQcekZmsLPh6sN9R30FBsZncX58NguVEZmD4U556cgKxtsNxg+EIz0cdyf7YrhTXjrcMwqvR3BNbflVf21rbRn8Xg+HZsjWGO6Ulw6fG8PWmlIU+71X/bWFPi9aa8sY7mRrDHfKO4mk4mjP2IomUw07GypwrG8MsUTSxMqIzMNwp7zTMTiF6HwCO+uvfrzdsLOhErOxJE6dnzSxMiLzMNwp7xg3IK1kpYxhZ7rXf5jr3cmmGO6Udw6fG0N5cQGaQoEVv0ZdZTFCQT/vVCXbYrhT3jmSHm/PZLOwpYgIdtRX8k5Vsi2GO+WV2VgC7YOTuLbu6pdALnZtXTm6IlFMzcVNqIzIXAx3yisdg1NIKrBl3fL7ty9ny7pSAED7BU6qkv0w3CmvnE4H8ZZ1wVW/1qVwn1r1axGZjeFOeeXUhUkUeAUbqlY+mWqoryxBUYEHp9hzJxtiuFNeOX1+EhvDQRR4V/+t7/EIWtaWXvxtgMhOGO6UV05fmMLm6lLTXm9zdSlvZCJbYrhT3picjaFvbObiWLkZtqwLYnByDmPT86a9JpEZGO6UN9oHUxOfZvfcgdRvBER2wnCnvHE6PXyyuXr1K2UMRrhzUpXsJqNwF5HbROSUiHSIyEOX+XyDiPxaRA6LyDERucP8UolW59SFSRQVeFBfWWLaa9aUF6G00HfxBweRXSwb7iLiBfAIgNsBtALYLSKti5r9KYAnVXUngHsBfNXsQolW6/SFSWyuLoXHs/JtBxYTEWxeV8qeO9lOJj33GwB0qGqXqs4DeALArkVtFIBxy185gH7zSiQyh9krZQybq0vRfmESqmr6axOtVCbhvh5Az4LHvennFvoLAB8UkV4A+wB8wpTqiEwyEp1HZHIOW7IQ7luqgxidjiEyNWf6axOtlFkTqrsBfEtV6wDcAeAxEXnda4vI/SJyQEQORCIRky5NtDzjRqMWEydTDRdXzJznihmyj0zCvQ9A/YLHdennFroPwJMAoKq/AVAEILT4hVR1j6q2qWpbOBxeWcVEK3BpT5ksDMus44oZsp9Mwv1lAC0i0iQifqQmTPcuanMOwDsAQES2IhXu7JqTbZw6P4nSIh/WlRWZ/tqhYCGqAn6umCFbWTbcVTUO4EEATwM4idSqmOMi8gURuSfd7DMAPioiRwE8DuAjytklspH2C1PYUl26qgM6rmRzdSlODzLcyT58mTRS1X1ITZQufO7hBR+fAPAWc0sjMoeq4tSFSdyxrSZr19hcHcQPDvZCVbP2A4ToavAOVXK9yNQcxmdiaFlr/mSqoaW6FNH5BAbGZ7N2DaKrwXAn1+uKRAEAG7MY7s3hwO9ci8hqDHdyve6hVOA2h1Z/QMdSmkPB9LW4HJLsgeFOrtcVmYLf50FtRXHWrlFdVoiA34tO9tzJJhju5HpdkSiaqgLwmrinzGIigqZwAF1DDHeyB4Y7uV7XUPTimHg2NYeC6IpwWIbsgeFOrhZLJHFuZDon4d4UCqBvbAazsUTWr0W0HIY7udq5kWkkkoqmUPZWyhiawwGoAmeHp7N+LaLlMNzJ1YylibnouW8MB9PX5NAMWY/hTq5mLE3cmIOee2N6qSUnVckOGO7kal2RKKoCfpSXFGT9WsFCH6rLCnkjE9kCw51crSsSRVMWb15arDkURBdvZCIbYLiTq3UNTeVkvN3QFA6gKxLlkXtkOYY7udb4TAxDU/NoDmd/vN3QHApgfCaG0elYzq5JdDkMd3KtXOwpsxhXzJBdMNzJtYyAzeWwDHeHJLtguJNrdQ9F4fUIGtbkLtzXVxSjwCtcDkmWY7iTa3VFoqivLIbfl7tvc5/Xgw1VAQ7LkOUY7uRanZGpnC6DNDSHuDskWY/hTq6UTCrODEdzulLG0BQO4OxwFIkkl0OSdRju5EoDE7OYjSUt6blvDAURSyh6R7mBGFmH4U6udMaCZZAGY4+Zbg7NkIUY7uRKRrA25XAZpMH4beEMw50sxHAnV+oeiqKowIPq0qKcXzsU9CNY6GPPnSzFcCdXOjMURWNVAJ4snpu6FBFBY6gE3Ty0gyzEcCdX6h5OhbtVGqsCHJYhSzHcyXXiiSR6RqYvTmxaoSkUQO/oNObjSctqoPzGcCfX6R+bRSyhlqyUMTSFAkgq0MPlkGQRhju5jnFYhpU994vLIbmBGFmE4U6uY4x1N4ZKLKuhKT3ef2aY4U7WyCjcReQ2ETklIh0i8tASbf6TiJwQkeMi8j1zyyTK3JnhaQT8XoSDhZbVUBnwo7y4gMshyTK+5RqIiBfAIwBuBdAL4GUR2auqJxa0aQHwOQBvUdVREVmbrYKJltM9FEVjKACR3C+DXKgxFGDPnSyTSc/9BgAdqtqlqvMAngCwa1GbjwJ4RFVHAUBVB80tkyhzZ4Zzeyj2UppDAZwZ4oQqWSOTcF8PoGfB4970cwttBrBZRF4Qkf0icptZBRJdjfl4ahmkHcK9sSqAvrEZzMYSVpdCecisCVUfgBYAtwDYDeBREalY3EhE7heRAyJyIBKJmHRpokt6RqeRVFh6A5PBmNA9yztVyQKZhHsfgPoFj+vSzy3UC2CvqsZUtRvAaaTC/neo6h5VbVPVtnA4vNKaiZZ0aaWM9eHexN0hyUKZhPvLAFpEpElE/ADuBbB3UZsfI9Vrh4iEkBqm6TKxTqKMXNwN0gbhbvyA4aQqWWHZcFfVOIAHATwN4CSAJ1X1uIh8QUTuSTd7GsCwiJwA8GsA/11Vh7NVNNFSzgxHUV5cgMqSAqtLQVlRAUJBP/eYIUssuxQSAFR1H4B9i557eMHHCuDT6X+ILGOXZZCGxiqep0rW4B2q5CpnhqbRVGXdnamLNYa4OyRZg+FOrjEbS6B/fMYWk6mGplAAg5NziM7FrS6F8gzDnVzj3Mg0VO0xmWrgihmyCsOdXKMrYhyKHbS4kkuauGKGLMJwJ9fotsFukIsZN1Nx61/KNYY7uUb30BTCpYUoLbJ+GaSh2O9FbXkRh2Uo5xju5BpdEXtsGLZYUziAToY75RjDnVyjeyhq6dF6S2kKBdAdmULqdhCi3GC4kyuMT8cwHJ23Zc+9ORTExGwcI9F5q0uhPMJwJ1foHrbPnjKLNYW5HJJyj+FOrtCdPhS7OWyfZZAGY6iI2xBQLjHcyRW6I1F4BGhYY59lkIb1FcUo8Ap77pRTDHdyhc6hKOrXlMDvs9+3tM/rQcOaEnRFpqwuhfKI/f4mEK1At02XQRqaw0H23CmnGO7keKqK7iGbh3sogDPD00gkuRyScoPhTo53YWIOM7GELde4G5pCAczHk+gfm7G6FMoTDHdyvC4br5QxcHdIyjWGOzmenc5NXQrXulOuMdzJ8boiURQVeLCurMjqUpYUDhaitNDHFTOUMwx3crzuoSgaqwLweOxxburliAiawjxPlXKH4U6O1z0URXPYvkMyhqZQgMMylDMMd3K0WCKJcyPTth5vNzSFAugbm8FsLGF1KZQHGO7kaD0jqbXjTTY6Wm8pTaEAVFNnvRJlG8OdHK3TODfVAcMyG9NLNTsHOalK2cdwJ0frSAflprX277kbP4A6GO6UAwx3crT2wUmsLS1EmY3OTV1Kid+H9RXFaGe4Uw4w3MnROgen0FJt/167oaU6yJ475QTDnRxLVdExOIVNNt52YLFN4SA6I1PcQIyyjuFOjjUwPovofAKbqkutLiVjm9YGMRdPom+UG4hRdjHcybEuTqY6qOduDCF1RCYtroTcLqNwF5HbROSUiHSIyENXaPdeEVERaTOvRKLLM8LdSWPum8Kp3zLaL3DcnbJr2XAXES+ARwDcDqAVwG4Rab1Mu1IAnwTwW7OLJLqc9sEpVJQUoCrgt7qUjJWXFCBcWshJVcq6THruNwDoUNUuVZ0H8ASAXZdp95cAvghg1sT6iJbUmZ5MFbHvhmGXsykc5HJIyrpMwn09gJ4Fj3vTz10kItcBqFfVn5pYG9EVtQ9OOmpIxrBpbRCdg1NQ5YoZyp5VT6iKiAfA3wL4TAZt7xeRAyJyIBKJrPbSlMeGp+YwOh27eEu/k7RUBzE5F8eFiTmrSyEXyyTc+wDUL3hcl37OUArgGgD/JiJnANwEYO/lJlVVdY+qtqlqWzgcXnnVlPcuTaY6ZxmkwVjdw3F3yqZMwv1lAC0i0iQifgD3AthrfFJVx1U1pKqNqtoIYD+Ae1T1QFYqJgIujlk7YU+ZxTYZyyEHuRySsmfZcFfVOIAHATwN4CSAJ1X1uIh8QUTuyXaBRJfTMTiFgN+L2nL7Hq23lHCwEGVFPk6qUlb5MmmkqvsA7Fv03MNLtL1l9WURXVnH4BQ2rnXeShkgdeTeprXcY4ayi3eokiN1DE45ckjG0LK2lOFOWcVwJ8eZnI3h/MSso8N909oghqPzGInOW10KuRTDnRzHiXvKLGb8YGLvnbKF4U6OY+zL4sRlkAbj5qvTF7hihrKD4U6Oc2JgAiV+LzasKbG6lBVbX1GM0iIfTg5MWF0KuRTDnRzn5MAE3rCuFB6P81bKGEQEW2vKGO6UNQx3chRVxYmBCWytKbO6lFVrrSnDa+cnkeSpTJQFDHdylL6xGUzOxl0R7ltrSjE9n8DZkWmrSyEXYriTo5wcSE1AuiPcU++BQzOUDQx3cpSTAxMQAd6wzrkrZQybq0vhEYY7ZQfDnRzlRP8EGqsCCBRmtHOGrRUVeLExHMSJfoY7mY/hTo5y8vwEttY4v9du4IoZyhaGOznG1FwcZ4ensXWd88fbDVtrytA/PouxaW5DQOZiuJNjnDqf6uG6YTLVYPwWYkwUE5mF4U6OYYxNt9a6J9yN93KCQzNkMoY7OcaJgUmUFxegxoEHdCxlbWkRQkE/x93JdAx3coyTA6nJVCce0HElnFSlbGC4kyMkkopT5yddNd5u2FpThvYLU4glklaXQi7CcCdHODMcxUwsgVYXhntrTRnmE0l0Rri3O5mH4U6OYAxbuLXnDvBOVTIXw50c4VjvOPxez8VDLtxkYziA4gIvjvaMW10KuQjDnRzh0NlRXLO+DIU+r9WlmM7n9WBbXTkO94xZXQq5CMOdbG8+nsQrfeO4rqHS6lKy5rqGSpzoH8dsLGF1KeQSDHeyvZMDE5iLJ7HTxeG+s6ECsYTieD+HZsgcDHeyvcPnRgGkAtCtjPd26CyHZsgcDHeyvcM9Y1hXVoTaimKrS8mataVFqKssxuGeUatLIZdguJPtHTo36upeu2FnQyUOn2PPnczBcCdbi0zOoWdkxtWTqYbrGiowMD6LgfEZq0shF2C4k63lw3i7wZgwZu+dzMBwJ1s73DOGAq/gmvXlVpeSda01ZfD7PDh0luPutHoZhbuI3CYip0SkQ0QeusznPy0iJ0TkmIg8IyIbzC+V8tGhs6NorSlDUYH7bl5azO/zYNt63sxE5lg23EXEC+ARALcDaAWwW0RaFzU7DKBNVbcD+AGAvzG7UMo/8UQSx3rHXb2+fbGd9RV4pW8c83HuEEmrk0nP/QYAHaraparzAJ4AsGthA1X9tapOpx/uB1BnbpmUj147P4mZWCIvxtsN122oxHw8yZOZaNUyCff1AHoWPO5NP7eU+wD8bDVFEQGXJlPzYaWM4dLNTBx3p9UxdUJVRD4IoA3Al5b4/P0ickBEDkQiETMvTS70YucwaspTN/fki5ryYtSvKcaLncNWl0IOl0m49wGoX/C4Lv3c7xCRdwL4PIB7VHXuci+kqntUtU1V28Lh8ErqpTwRTyTx7x1DuLkl7Lpj9ZZzc0sYv+kc4rg7rUom4f4ygBYRaRIRP4B7Aexd2EBEdgL4OlLBPmh+mZRvjvaOY3I2jrduDlldSs69tSWM6Hzi4rAU0UosG+6qGgfwIICnAZwE8KSqHheRL4jIPelmXwIQBPB9ETkiInuXeDmijDzfHoEI8JaN+Rfub95UBa9H8Hz7kNWlkIP5MmmkqvsA7Fv03MMLPn6nyXVRnnvudATb6ypQGfBbXUrOlRUVYGd9BZ5rj+Cz/3GL1eWQQ/EOVbKd8ZkYjvSM4eaW/Ou1G97aEsYrfeMYic5bXQo5FMOdbOfFjiEkFbh5c/5Out+8OQRV4N87ODRDK8NwJ9t5rn0IwUIfdtTnz81Li22vq0BZkQ/Pn+aSYVoZhjvZiqriudMRvGljFQq8+fvt6fUIfq8lhOfbh6CqVpdDDpS/f3vIlrqHougbm8nrIRnDzS1hnJ+YRfvglNWlkAMx3MlWjOV/+TyZavi99P+D5zg0QyvAcCdb2ffKAJrDAWyoClhdiuXqKkuwpboUP3v1vNWlkAMx3Mk2+sdm8NvuEey69kr70uWXe3bU4uDZUfSMTC/fmGgBhjvZxt6j/QCAd++stbgS+9i1I/X/4idHXredE9EVMdzJNn58uA/XNVRwSGaBusoS3NC0Bj863MdVM3RVGO5kCycHJvDa+Um8eyeHZBZ794716IxEcbyfB3hQ5hjuZAs/PtIHn0dw57Yaq0uxnTu31cDv9eBHhzk0Q5ljuJPlkknF3iP9uHlzGFXBQqvLsZ3ykgLcsiWMvUf7kUhyaIYyw3Any/22ewQD47MckrmCP9i5HpHJObzYyb1mKDMMd7Lc9w/2IOD34tat1VaXYlu//4a1KC3y4fsHeq0uhRyC4U6W6h2dxt4j/Xjf9XUo9nutLse2igq8uPc/1OOnrwzg7HDU6nLIARjuZKmvP9sFEeCP37bR6lJs76NvbYbXI/jas51Wl0IOwHAnywxOzOL/HujB+66vQ21FsdXl2N7asiK8v60ePzjYi/6xGavLIZtjuJNlHn2+C/FEEg+w156xP35bM1SBPc91WV0K2RzDnSwxEp3Hd/afw64d63lH6lWoqyzBe65bj8dfOofI5JzV5ZCNMdzJEt94vguz8QQ+dgt77Vfrv96yCbFEEnue49g7LY3hTjl3cmACjz7fhV3X1qKlutTqchynKRTAe6+rwz+/cAav9o1bXQ7ZFMOdcmo+nsSnnzyK8mI//vzuN1pdjmP96Z2tqAr68eknj2AunrC6HLIhhjvl1D/+qh0nBybw1+/ZhsqA3+pyHKu8pABffO92nL4whb/7ZbvV5ZANMdwpZ472jOGr/9aJ911fh3e28m7U1bply1rsvqEee57rxMGzI1aXQzbDcKec6B2dxse+ewjVpYV4+O5Wq8txjc/f2YraimI8+L3DODfM05roEoY7ZV3f2Ax2P7ofk7MxfP1DbSgrKrC6JNcIFvqw50NtmIklsPvR/TyOjy5iuFNW9Y/NYPee/RibjuGx+27Etrpyq0tyndbaMnznvhsxNRfHvXv2o3eUAU8Md8qiZ09H8N5/ehGj0Xk8dt+NuLa+wuqSXOua9eX4zn03YnI2hvd89UX8+rVBq0siizHcyXQTszE89MNj+PA3X0Kg0IfH778JOxjsWbetrhxP3P8mVJQU4D9/62V89vtHMT4Ts7ossogvk0YichuALwPwAviGqv71os8XAvg2gOsBDAN4v6qeMbdUsrsLE7P47v6z+N5L5zASnccDb9uIT72zBUUF3Mo3V1pry/Avn/g9/MMz7fjas1341WuD+KMbGvCBmxpQU87N2fKJLHeiuoh4AZwGcCuAXgAvA9itqicWtPkYgO2q+oCI3AvgD1T1/Vd63ba2Nj1w4MBq6ycLqSrOjUzjhY5hPHt6EM+cHERCFW/fshafeEcLe+sWe6V3HF9+ph3PvHYBHhG8/Q1rccuWMN68MYTGqhKIiNUl0gqIyEFVbVuuXSY99xsAdKhqV/qFnwCwC8CJBW12AfiL9Mc/APAVERFd7icH2VYyqZiOJRCdi2NqLo7hqXkMTc1hcGIWZ4an0RmZQvuFKZyfmAUAVJcV4iNvbsSH3rSBG4HZxLa6cnzjw23oGZnGY/vPYu+RfvzyxAUAqT+vzdWlaA4F0BQKYG1ZEaoCflQFC1Fa5EOg0IeSAi88Hv4AcKpMwn09gJ4Fj3sB3LhUG1WNi8g4gCoAph/4+OTLPXj0+fze7jTTn5iqeqmtXvo6VUVSAYUimQQSSUU8qUgkk5iPJzGfSCKWWPoqAb8XzeEgbmpeg+s3VOJNG0PYGA6wJ2hT9WtK8Cd3bMXnbn8DuoeieKFzGIfOjqIzMoUfHurD1Fx8ya8t8Ar8Xg/8Pg+8Hg98HoHXI/B4AIHAI7j45y4X/5X6T6bfD/n4XfPf3tGCu6+tzeo1MhpzN4uI3A/gfgBoaGhY0WtUlBSgpTpoZlmOJJn+lZBLf3lELn2VRwCPCEQEPo/A40n91+/zoMDrQaHPg0ChF4FCH4KFPqwJ+BEKFqIq6Ec4WMggdyARQXM4iOZwEB+6aQOA1A/6ofRvZcNT8xiOzmFqLo7oXBzRuQTmE+kf+PEk4klFMt0RSHUQUp0EINVxMH5R14v/Wp5m3FVxl/Li7N/rkUm49wGoX/C4Lv3c5dr0iogPQDlSE6u/Q1X3ANgDpMbcV1Lwu964Du9647qVfCkRLSIiCJcWIlxaaHUpZLJMlkK+DKBFRJpExA/gXgB7F7XZC+DD6Y/fB+BXHG8nIrLOsj339Bj6gwCeRmop5DdV9biIfAHAAVXdC+B/A3hMRDoAjCD1A4CIiCyS0Zi7qu4DsG/Rcw8v+HgWwB+aWxoREa0U71AlInIhhjsRkQsx3ImIXIjhTkTkQgx3IiIXWnbjsKxdWCQC4KwlF1+dELKwrYID5OP75nvOH0563xtUNbxcI8vC3alE5EAmO7K5TT6+b77n/OHG981hGSIiF2K4ExG5EMP96u2xugCL5OP75nvOH6573xxzJyJyIfbciYhciOG+CiLyGRFREQlZXUu2iciXROQ1ETkmIj8SEVcfkCoit4nIKRHpEJGHrK4n20SkXkR+LSInROS4iHzS6ppyRUS8InJYRJ6yuhYzMdxXSETqAbwLwDmra8mRXwK4RlW3I3Vg+ucsridr0ofCPwLgdgCtAHaLSKu1VWVdHMBnVLUVwE0APp4H79nwSQAnrS7CbAz3lfs7AP8DmR9p6miq+gtVNQ7b3I/UiVxudfFQeFWdB2AcCu9aqjqgqofSH08iFXbrra0q+0SkDsCdAL5hdS1mY7ivgIjsAtCnqketrsUi/wXAz6wuIosudyi864POICKNAHYC+K21leTE3yPVSUtaXYjZcnpAtpOIyL8CuNxhrZ8H8CdIDcm4ypXes6r+JN3m80j9Cv/dXNZGuSEiQQA/BPApVZ2wup5sEpG7AAyq6kERucXqeszGcF+Cqr7zcs+LyDYATQCOigiQGp44JCI3qOr5HJZouqXes0FEPgLgLgDvcPkZuZkcCu86IlKAVLB/V1X/n9X15MBbANwjIncAKAJQJiLfUdUPWlyXKbjOfZVE5AyANlV1yqZDKyIitwH4WwBvU9WI1fVkk4j4kJo0fgdSof4ygD9S1eOWFpZFkuqp/B8AI6r6KavrybV0z/2zqnqX1bWYhWPulKmvACgF8EsROSIiX7O6oGxJTxwbh8KfBPCkm4M97S0APgTg7ek/3yM61I9yAAAAPElEQVTpHi05FHvuREQuxJ47EZELMdyJiFyI4U5E5EIMdyIiF2K4ExG5EMOdiMiFGO5ERC7EcCcicqH/D5JGMqHWx49eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# pi(q): probability distribution over |R^2\n",
    "# piH(q, p) = pi_pq(p|q) pi(p)\n",
    "# H = - log piH(q, p) = - log pi_pq(p|q) - log pi(q)\n",
    "# H =                   - K              - V\n",
    "# dq/dt =  dH/dp\n",
    "# dp/dt = -dH/dq \n",
    "def hmc(pi, q, p):\n",
    "    accept = 0\n",
    "    momentum = st.norm(0, 1)\n",
    "    qs = []\n",
    "    for _ in range(1000):\n",
    "        # H = K(pi, p, q) + V(pi, q)\n",
    "        prob = -np.log(pi(q).detach().numpy()) - np.sum(momentum.logpdf(p.detach().numpy()))\n",
    "        [q_new, p_new] = leapfrog(q, p, mkPiVdq(pi), 1.0)\n",
    "        \n",
    "        # for time reversal\n",
    "        p_new = -p_new\n",
    "        prob2 = -np.log(pi(q_new).detach().numpy()) - np.sum(momentum.logpdf(p_new.detach().numpy()))\n",
    "        \n",
    "        if np.random.rand() < prob2 / prob:\n",
    "            q = q_new\n",
    "            p = p_new\n",
    "            accept += 1\n",
    "        qs.append(q.detach().numpy())\n",
    "    print(\"acceptance ratio: %4.2f\" % (accept / len(qs)))\n",
    "    return np.array(qs)\n",
    "\n",
    "\n",
    "def prob(x):\n",
    "    return torch.pow(5.0, -(x[0] * x[0]))\n",
    "\n",
    "points = hmc(prob, torch.tensor([2.0]), torch.tensor([0.1]))\n",
    "plt.scatter(points, [1 for _ in range(len(points))], c='g')\n",
    "xs = np.arange(-5, 5, 0.1)\n",
    "ys = [prob(torch.tensor([x])) for x in xs]\n",
    "\n",
    "plt.plot(xs, ys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 481,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([1., 1.]) tensor([0., 1.])\n",
      "tensor([1., 1.]) tensor([-0.5000, -0.5000])\n",
      "tensor([1.0625, 1.5625]) tensor([-0.2976, -0.4376])\n",
      "tensor([1.0625, 1.5625]) tensor([-0.2976, -0.4376])\n",
      "tensor([1.1994, 2.2344]) tensor([-0.1865, -0.3474])\n",
      "tensor([1.1994, 2.2344]) tensor([-0.1865, -0.3474])\n",
      "tensor([1.3829, 2.9932]) tensor([-0.1272, -0.2753])\n",
      "tensor([1.3829, 2.9932]) tensor([-0.1272, -0.2753])\n",
      "tensor([1.5982, 3.8208]) tensor([-0.0932, -0.2228])\n",
      "tensor([1.5982, 3.8208]) tensor([-0.0932, -0.2228])\n",
      "tensor([1.8369, 4.7041]) tensor([-0.0720, -0.1845])\n",
      "tensor([1.8369, 4.7041]) tensor([-0.0720, -0.1845])\n",
      "tensor([2.0935, 5.6335]) tensor([-0.0580, -0.1560])\n",
      "tensor([2.0935, 5.6335]) tensor([-0.0580, -0.1560])\n",
      "tensor([2.3646, 6.6019]) tensor([-0.0481, -0.1342])\n",
      "tensor([2.3646, 6.6019]) tensor([-0.0481, -0.1342])\n",
      "tensor([2.6477, 7.6038]) tensor([-0.0408, -0.1173])\n",
      "tensor([2.6477, 7.6038]) tensor([-0.0408, -0.1173])\n",
      "tensor([2.9411, 8.6351]) tensor([-0.0353, -0.1038])\n",
      "tensor([2.9411, 8.6351]) tensor([-0.0353, -0.1038])\n",
      "tensor([3.2433, 9.6923]) tensor([-0.0310, -0.0928])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEKlJREFUeJzt3X2MZXV9x/H3pyxggkTQXWFZFoF0U8TWVrxBUNPQgg1sDOsTCfxRwYdMbSXVpElDS6KJ/1TbxCYGKt0gERoDWBRZ61oEwVDTgAxkeVhWZCE27LrKCA1KtNjVb/+4P+xkuLMzu+fMvXfx/Upu5jz85ny/nBn2c8/DPZOqQpKk35p0A5Kk6WAgSJIAA0GS1BgIkiTAQJAkNQaCJAnoIRCSrE9yZ5JHkmxP8pERY5LkM0l2JnkwyWld60qS+rWqh23sBf6qqu5PciRwX5LbquqReWPOAza015uAz7avkqQp0fkIoar2VNX9bfqnwA5g3YJhm4Drauhu4Kgka7vWliT1p48jhF9LciLwBuCeBavWAU/Om9/Vlu0ZsY0ZYAbgiCOOeOMpp5zSZ4uS9JJ23333/biq1hzI9/YWCEleDnwJ+GhV/eRAt1NVm4HNAIPBoGZnZ3vqUJJe+pL814F+by93GSU5lGEYfKGqvjxiyG5g/bz549sySdKU6OMuowCfA3ZU1acXGbYFeG+72+gM4NmqetHpIknS5PRxyugtwJ8CDyXZ1pb9LXACQFVdBWwFNgI7gZ8B7+uhriSpR50Doaq+DWSJMQV8uGstSdLK8ZPKkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDW9BEKSa5I8leThRdafleTZJNva62N91JUk9afz31RuPg9cAVy3jzH/UVVv76meJKlnvRwhVNVdwDN9bEuSNBnjvIZwZpIHknw9yevGWFeStAx9nTJayv3Aa6rquSQbga8AG0YNTDIDzACccMIJY2pPkjSWI4Sq+klVPdemtwKHJlm9yNjNVTWoqsGaNWvG0Z4kiTEFQpJjk6RNn97qPj2O2pKk5enllFGS64GzgNVJdgEfBw4FqKqrgPcAf55kL/Bz4MKqqj5qS5L60UsgVNVFS6y/guFtqZKkKeUnlSVJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEtBTICS5JslTSR5eZH2SfCbJziQPJjmtj7qSpP70dYTweeDcfaw/D9jQXjPAZ3uqK0nqSS+BUFV3Ac/sY8gm4Loauhs4KsnaPmpLkvoxrmsI64An583vasteJMlMktkks3Nzc2NpTpI0hReVq2pzVQ2qarBmzZpJtyNJvzHGFQi7gfXz5o9vyyRJU2JcgbAFeG+72+gM4Nmq2jOm2pKkZVjVx0aSXA+cBaxOsgv4OHAoQFVdBWwFNgI7gZ8B7+ujriSpP70EQlVdtMT6Aj7cRy1J0sqYuovKkqTJMBAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJElAT4GQ5NwkjybZmeSyEesvSTKXZFt7fbCPupKk/nT+m8pJDgGuBN4G7ALuTbKlqh5ZMPTGqrq0az1J0sro4wjhdGBnVT1RVb8AbgA29bBdSdIY9REI64An583vassWeneSB5PclGT9YhtLMpNkNsns3NxcD+1JkpZjXBeVvwqcWFWvB24Drl1sYFVtrqpBVQ3WrFkzpvYkSX0Ewm5g/jv+49uyX6uqp6vq+TZ7NfDGHupKknrURyDcC2xIclKSw4ALgS3zByRZO2/2fGBHD3UlST3qfJdRVe1NcilwK3AIcE1VbU/yCWC2qrYAf5nkfGAv8AxwSde6kqR+paom3cOiBoNBzc7OTroNSTpoJLmvqgYH8r1+UlmSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkppdASHJukkeT7Exy2Yj1hye5sa2/J8mJfdSVJPWncyAkOQS4EjgPOBW4KMmpC4Z9APjvqvpt4B+BT3WtK0nqVx9HCKcDO6vqiar6BXADsGnBmE3AtW36JuDsJOmhtiSpJ30EwjrgyXnzu9qykWOqai/wLPCqURtLMpNkNsns3NxcD+1JkpZj6i4qV9XmqhpU1WDNmjWTbkeSfmP0EQi7gfXz5o9vy0aOSbIKeAXwdA+1JUk96SMQ7gU2JDkpyWHAhcCWBWO2ABe36fcAd1RV9VBbktSTVV03UFV7k1wK3AocAlxTVduTfAKYraotwOeAf0myE3iGYWhIkqZI50AAqKqtwNYFyz42b/p/gAv6qCVJWhlTd1FZkjQZBoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAnoGAhJXpnktiSPta9HLzLul0m2tdeWLjUlSSuj6xHCZcA3q2oD8M02P8rPq+oP2uv8jjUlSSugayBsAq5t09cC7+i4PUnShHQNhGOqak+b/iFwzCLjXpZkNsndSfYZGklm2tjZubm5ju1JkpZr1VIDktwOHDti1eXzZ6qqktQim3lNVe1OcjJwR5KHqurxUQOrajOwGWAwGCy2PUlSz5YMhKo6Z7F1SX6UZG1V7UmyFnhqkW3sbl+fSPIt4A3AyECQJE1G11NGW4CL2/TFwC0LByQ5OsnhbXo18BbgkY51JUk96xoInwTeluQx4Jw2T5JBkqvbmNcCs0keAO4EPllVBoIkTZklTxntS1U9DZw9Yvks8ME2/Z/A73WpI0laeX5SWZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqSmUyAkuSDJ9iS/SjLYx7hzkzyaZGeSy7rUlCStjK5HCA8D7wLuWmxAkkOAK4HzgFOBi5Kc2rGuJKlnq7p8c1XtAEiyr2GnAzur6ok29gZgE/BIl9qSpH6N4xrCOuDJefO72rKRkswkmU0yOzc3t+LNSZKGljxCSHI7cOyIVZdX1S19N1RVm4HNAIPBoPreviRptCUDoarO6VhjN7B+3vzxbZkkaYqM45TRvcCGJCclOQy4ENgyhrqSpP3Q9bbTdybZBZwJfC3JrW35cUm2AlTVXuBS4FZgB/DFqtrerW1JUt+63mV0M3DziOU/ADbOm98KbO1SS5K0svyksiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJElN17+pfEGS7Ul+lWSwj3HfT/JQkm1JZrvUlCStjE5/Uxl4GHgX8M/LGPtHVfXjjvUkSSukUyBU1Q6AJP10I0mamHFdQyjgG0nuSzIzppqSpP2w5BFCktuBY0esuryqbllmnbdW1e4krwZuS/LdqrprkXozwAzACSecsMzNS5K6WjIQquqcrkWqanf7+lSSm4HTgZGBUFWbgc0Ag8GgutaWJC3Pip8ySnJEkiNfmAb+hOHFaEnSFOl62+k7k+wCzgS+luTWtvy4JFvbsGOAbyd5APgO8LWq+vcudSVJ/et6l9HNwM0jlv8A2NimnwB+v0sdSdLK85PKkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJ6BgISf4hyXeTPJjk5iRHLTLu3CSPJtmZ5LIuNSVJK6PrEcJtwO9W1euB7wF/s3BAkkOAK4HzgFOBi5Kc2rGuJKlnnQKhqr5RVXvb7N3A8SOGnQ7srKonquoXwA3Api51JUn9W9Xjtt4P3Dhi+TrgyXnzu4A3LbaRJDPATJt9PsnDvXW4MlYDP550E8tgn/2yz37ZZ39+50C/cclASHI7cOyIVZdX1S1tzOXAXuALB9rIC6pqM7C5bXe2qgZdt7mSDoYewT77Zp/9ss/+JJk90O9dMhCq6pwlil8CvB04u6pqxJDdwPp588e3ZZKkKdL1LqNzgb8Gzq+qny0y7F5gQ5KTkhwGXAhs6VJXktS/rncZXQEcCdyWZFuSqwCSHJdkK0C76HwpcCuwA/hiVW1f5vY3d+xvHA6GHsE++2af/bLP/hxwjxl9lkeS9JvGTypLkgADQZLUTFUgHAyPwkhyQZLtSX6VZNHbz5J8P8lD7drKAd8GdqD2o8+JPlYkySuT3Jbksfb16EXG/bLty21JxnZTwlL7J8nhSW5s6+9JcuK4elvQx1J9XpJkbt4+/OAEerwmyVOLfbYoQ59p/w0PJjlt3D22Ppbq86wkz87blx+bQI/rk9yZ5JH2//lHRozZ//1ZVVPzAv4EWNWmPwV8asSYQ4DHgZOBw4AHgFPH2ONrGX7w41vAYB/jvg+snuC+XLLPSe/L1sPfA5e16ctG/czbuucmsA+X3D/AXwBXtekLgRuntM9LgCvG3duCHv4QOA14eJH1G4GvAwHOAO6Z0j7PAv5twvtyLXBamz6S4aODFv7M93t/TtURQh0Ej8Koqh1V9ei46h2oZfY5DY8V2QRc26avBd4x5vr7spz9M7//m4Czk2SMPcJ0/ByXVFV3Ac/sY8gm4Loauhs4Ksna8XT3/5bR58RV1Z6qur9N/5ThHZzrFgzb7/05VYGwwPsZpttCox6FsXBHTIMCvpHkvvY4jmk0DfvymKra06Z/CByzyLiXJZlNcneScYXGcvbPr8e0NzPPAq8aS3cjemgW+zm+u506uCnJ+hHrJ20afh+X68wkDyT5epLXTbKRdpryDcA9C1bt9/7s81lGyzLuR2EciOX0uAxvrardSV7N8HMa323vPHrTU58rbl99zp+pqkqy2H3Qr2n782TgjiQPVdXjfff6EvZV4Pqqej7JnzE8qvnjCfd0sLqf4e/jc0k2Al8BNkyikSQvB74EfLSqftJ1e2MPhDoIHoWxVI/L3Mbu9vWpJDczPKzvNRB66HMsjxXZV59JfpRkbVXtaYezTy2yjRf25xNJvsXwHdFKB8Jy9s8LY3YlWQW8Anh6hftaaMk+q2p+T1czvHYzbQ6Kx9zM/4e3qrYm+ackq6tqrA+9S3IowzD4QlV9ecSQ/d6fU3XKKC+RR2EkOSLJkS9MM7xYPo1PbZ2GfbkFuLhNXwy86MgmydFJDm/Tq4G3AI+Mobfl7J/5/b8HuGORNzIrack+F5w7Pp/hOedpswV4b7s75gzg2XmnE6dGkmNfuE6U5HSG/46O9U1Aq/85YEdVfXqRYfu/Pyd5pXzElfOdDM95bWuvF+7eOA7YuuDq+fcYvkO8fMw9vpPhubjngR8Bty7skeHdHg+01/Zx97jcPie9L1v9VwHfBB4Dbgde2ZYPgKvb9JuBh9r+fAj4wBj7e9H+AT7B8E0LwMuAf22/u98BTh73Plxmn3/XfhcfAO4ETplAj9cDe4D/bb+bHwA+BHyorQ/DP6b1ePs5L3oX34T7vHTevrwbePMEenwrw+uUD87793Jj1/3poyskScCUnTKSJE2OgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDX/B82a8jdTShV1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def gravity_simulation():\n",
    "    def gravitation_potential(q): return -1.0 / torch.norm(q)\n",
    "    def gravitational_force(q):\n",
    "        qnorm = q.norm()\n",
    "        qdir = q / qnorm\n",
    "        f =  qdir / qnorm\n",
    "        print(q, f)\n",
    "        return f\n",
    "    def gravitation_potential_dq(q):\n",
    "        q.detach_()\n",
    "        q.requires_grad_()\n",
    "        vcur = gravitation_potential(q)\n",
    "        vcur.backward()\n",
    "        print(q, q.grad)    \n",
    "        return q.grad\n",
    "\n",
    "    q = torch.tensor([1.0, 1.0])\n",
    "    p = torch.tensor([0.0, 1.0])\n",
    "    qs = [q.detach().numpy()]\n",
    "    print(q, p)\n",
    "    for i in range(10):\n",
    "        [q, p] = leapfrog(q, p, gravitational_force, 1.0 / 2.0)\n",
    "\n",
    "        qs.append(q.detach().numpy())\n",
    "    qs = np.array(qs)\n",
    "    plt.scatter(qs[:, 0], qs[:, 1])\n",
    "    plt.xlim(-2, 2)\n",
    "    plt.ylim(-2, 2)\n",
    "\n",
    "gravity_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}