def leapfrog(p0, q0, N, epsilon, mgradU):
    """
    Leapfrog integration step
    :param p0: initial momentum variables
    :param q0: initial coordinate variables (thetas)
    :param N: number of leapfrog steps to take
    :paran eps: step size for HMC steps
    :param mgradU: callable giving grad(log(posterior)) w/rt theta
    """
    # Initial half-step in momentum
    p = p0 + 0.5*epsilon * mgradU(q0)
    # Full step in coordinates
    q = q0 + epsilon * p
    # If more than one leapfrog step desired, alternate additional steps
    # in p and q until desired trajectory length reached
    for i in range(N-1):
        p = p + epsilon * mgradU(q)
        q = q + epsilon * p
    # Final half-step in momentum, with reflection (for reversibility)
    p = p + 0.5*epsilon * mgradU(q)
    # Return
    return p, q