docs: def expand_tpm(tpm):
      """Broadcast a state-by-node TPM so that singleton dimensions are expanded
      over the full network.
      """
      unconstrained = np.ones([2] * (tpm.ndim - 1) + [tpm.shape[-1]])
      return tpm * unconstrained