docs: def jsd(p1, p2):
      """Compute Jensen-Shannon divergence between p1 and p2.
      It assumes that p1 and p2 are already normalized that each of them sums to 1.
      """
      m = (p1 + p2) / 2
      return (kld(p1, m) + kld(p2, m)) / 2