- To perform hamiltonian monte carlo, we use the hamiltonian and its derivatives to provide a momentum to our proposal distribution --- That is, when we choose a new point from the current point, our probability distribution for the new point is influenced by our current momentum.

- For some integral necessary within this scheme, Euler integration doesn't cut it since the error diverges to infinity

- Hence, we need an integrator that guarantees that the energy of out system is conserved. Enter the leapfrog integrator. This integrator is also
*time reversible*-- We can run it forward for`n`

steps, and then run it backward for`n`

steps to arrive at the same state. Now I finally know how Braid was implemented, something that bugged the hell out of 9th grade me when I tried to implement Braid-like physics in my engine!

- The actual derivation of the integrator uses Lie algebras, Sympletic geometry, and other diffgeo ideas, which is great, because it gives me motivation to study differential geometry
`:)`

- Original paper: Construction of higher order sympletic integrators

```
import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg
# dq/dt = dH/dp | dp/dt = -dH/dq (a = -del V)
def leapfroge(dhdp, dhdq, q, p, dt):
p += -dhdq(q, p) * 0.5 * dt # halfstep momentum
q += dhdp(q, p) * dt # fullstep position
p += -dhdq(q, p) * 0.5 * dt # halfstep momentum
return (q, p)
def euler(dhdp, dhdq, q, p, dt):
pnew = p + -dhdq(q, p) * dt
qnew = q + dhdp(q, p) * dt
def planet(integrator, n, dt):
STRENGTH = 0.5
q = np.array([0.0, 1.0]); p = np.array([-1.0, 0.0])
# H = STRENGTH * |q| (potential) + p^2/2 (kinetic)
def H(qcur, pcur): return STRENGTH * np.linalg.norm(q) + np.dot(p, p) / 2
def dhdp(qcur, pcur): return p
def dhdq(qcur, pcur): return STRENGTH * 2 * q / np.linalg.norm(q)
qs = []
for i in range(n):
(q, p) = integrator(dhdp, dhdq, q, p, dt)
qs.append(q.copy())
return np.asarray(qs)
NITERS = 15
TIMESTEP = 1
plt.rcParams.update({'font.size': 22, 'font.family':'monospace'})
fig, ax = plt.subplots()
planet_leapfrog = planet(leapfroge, NITERS, TIMESTEP)
ax.plot(planet_leapfrog[:, 0], planet_leapfrog[:, 1], label='leapfrog',
linewidth=3, color='#00ACC1')
planet_euler = planet(euler, NITERS, TIMESTEP)
ax.plot(planet_euler[:, 0], planet_euler[:, 1], label='euler',
linewidth=3, color='#D81B60')
legend = plt.legend(frameon=False)
ax.set_title("leapfrog v/s euler: NITERS=%s dt=%s" % (NITERS, TIMESTEP))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
plt.show()
plt.savefig("leapfrog-vs-euler.png")
```