diff --git a/qualtran/bloqs/chemistry/thc/walk_operator.py b/qualtran/bloqs/chemistry/thc/walk_operator.py
index 30de4297b0..d88827bb38 100644
--- a/qualtran/bloqs/chemistry/thc/walk_operator.py
+++ b/qualtran/bloqs/chemistry/thc/walk_operator.py
@@ -12,7 +12,10 @@
# See the License for the specific language governing permissions and
# limitations under the License.
"""Function for building a walk operator for the THC hamiltonian."""
+import attrs
+import numpy as np
from numpy.typing import NDArray
+from openfermion.resource_estimates.utils import QI
from qualtran.bloqs.block_encoding.lcu_block_encoding import SelectBlockEncoding
from qualtran.bloqs.chemistry.thc import PrepareTHC, SelectTHC
@@ -51,3 +54,59 @@ def get_walk_operator_for_thc_ham(
block_encoding = SelectBlockEncoding(select=sel, prepare=prep)
walk_op = QubitizationWalkOperator(block_encoding=block_encoding)
return walk_op
+
+
+def get_reiher_thc_walk_operator(
+ num_bits_theta: int = 16, num_bits_state_prep: int = 10
+) -> QubitizationWalkOperator:
+ """Build the THC walk operator for the Reiher hamiltoninan
+
+ Note currently we spoof the Hamiltonian by over writing prepare's 1-norm
+ value with the correct value and use random THC factors for expediency.
+
+ Parameters are taken from openfermion compute_cost_thc_test.py.
+
+ Args:
+ num_bits_theta: the number of bits of precision for the givens rotations
+ num_bits_state_prep: The number of bits of precision for the preparation
+ of the LCU coefficients using alias sampling.
+
+ Returns:
+ walk_op: A constructed Reiher Hamiltonian walk operator.
+ """
+ # Let's just generate some random coefficients for the moment with parameters
+ # corresponding to the FeMoCo model complex.
+ num_spin_orb = 108
+ num_mu = 350
+ num_bits_theta = 16
+ num_bits_state_prep = 10
+ tpq = np.random.normal(0, 1, size=(num_spin_orb // 2, num_spin_orb // 2))
+ zeta = np.random.normal(0, 1, size=(num_mu, num_mu))
+ zeta = 0.5 * (zeta + zeta.T)
+ eta = np.random.normal(0, 1, size=(num_mu, num_spin_orb // 2))
+ eri_thc = np.einsum("Pp,Pr,Qq,Qs,PQ->prqs", eta, eta, eta, eta, zeta, optimize=True)
+ # In practice one typically uses the exact ERI tensor instead of that from
+ # THC, but that's a minor detail.
+ tpq_prime = (
+ tpq
+ - 0.5 * np.einsum("illj->ij", eri_thc, optimize=True)
+ + np.einsum("llij->ij", eri_thc, optimize=True)
+ )
+ t_l = np.linalg.eigvalsh(tpq_prime)
+ qroam_blocking_factor = QI(num_mu + num_spin_orb // 2)[0]
+ walk_op = get_walk_operator_for_thc_ham(
+ t_l,
+ eta,
+ zeta,
+ num_bits_state_prep=num_bits_state_prep,
+ num_bits_theta=num_bits_theta,
+ kr1=qroam_blocking_factor,
+ kr2=qroam_blocking_factor,
+ )
+ # GIANT HACK: overwrite the lambda value directly
+ # TODO: maybe parse THC hamiltonian files from openfermion directly
+ block_encoding = attrs.evolve(
+ walk_op.block_encoding, prepare=attrs.evolve(walk_op.prepare, sum_of_l1_coeffs=306.3)
+ )
+ walk_op = attrs.evolve(walk_op, block_encoding=block_encoding)
+ return walk_op
diff --git a/qualtran/bloqs/phase_estimation/phase_estimation_of_quantum_walk.ipynb b/qualtran/bloqs/phase_estimation/phase_estimation_of_quantum_walk.ipynb
index 6937510fa2..87c8a58a49 100644
--- a/qualtran/bloqs/phase_estimation/phase_estimation_of_quantum_walk.ipynb
+++ b/qualtran/bloqs/phase_estimation/phase_estimation_of_quantum_walk.ipynb
@@ -17,7 +17,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
@@ -41,7 +41,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
@@ -104,9 +104,43 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 13,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "m_0: ──────────chi_m[1]───@───────────────────────────────────────────────────────────────qft^-1───\n",
+ " │ │ │\n",
+ "m_1: ──────────chi_m[2]───┼───@(0)─────────@(0)───────────────────────────────────────────#2───────\n",
+ " │ │ │ │ │\n",
+ "m_2: ──────────chi_m[3]───┼───┼────────────┼──────@(0)─────────@(0)───────────────────────#3───────\n",
+ " │ │ │ │ │ │ │\n",
+ "m_3: ──────────chi_m[4]───┼───┼────────────┼──────┼────────────┼──────@(0)─────────@(0)───#4───────\n",
+ " │ │ │ │ │ │ │\n",
+ "selection0: ──────────────W───R_L────W^2───R_L────R_L────W^4───R_L────R_L────W^8───R_L─────────────\n",
+ " │ │ │ │ │ │ │ │ │ │\n",
+ "selection1: ──────────────W───R_L────W^2───R_L────R_L────W^4───R_L────R_L────W^8───R_L─────────────\n",
+ " │ │ │ │ │ │ │ │ │ │\n",
+ "selection2: ──────────────W───R_L────W^2───R_L────R_L────W^4───R_L────R_L────W^8───R_L─────────────\n",
+ " │ │ │ │ │ │ │ │ │ │\n",
+ "selection3: ──────────────W───R_L────W^2───R_L────R_L────W^4───R_L────R_L────W^8───R_L─────────────\n",
+ " │ │ │ │\n",
+ "target0: ─────────────────W──────────W^2─────────────────W^4─────────────────W^8───────────────────\n",
+ " │ │ │ │\n",
+ "target1: ─────────────────W──────────W^2─────────────────W^4─────────────────W^8───────────────────\n",
+ " │ │ │ │\n",
+ "target2: ─────────────────W──────────W^2─────────────────W^4─────────────────W^8───────────────────\n",
+ " │ │ │ │\n",
+ "target3: ─────────────────W──────────W^2─────────────────W^4─────────────────W^8───────────────────\n",
+ " │ │ │ │\n",
+ "target4: ─────────────────W──────────W^2─────────────────W^4─────────────────W^8───────────────────\n",
+ " │ │ │ │\n",
+ "target5: ─────────────────W──────────W^2─────────────────W^4─────────────────W^8───────────────────\n"
+ ]
+ }
+ ],
"source": [
"num_sites: int = 6\n",
"eps: float = 1e-2\n",
@@ -128,9 +162,20 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 14,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "GateCounts(t=0, toffoli=0, cswap=168, and_bloq=1474, clifford=11825, rotation=84, measurement=1474)"
+ ]
+ },
+ "execution_count": 14,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"from qualtran import Bloq\n",
"from qualtran.resource_counting import get_cost_value, QECGatesCost, GateCounts\n",
@@ -164,9 +209,19 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 15,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 800 ms, sys: 12.3 ms, total: 813 ms\n",
+ "Wall time: 813 ms\n",
+ "cswap: 295362, and_bloq: 21321353, clifford: 124158198, rotation: 65636, measurement: 21321353\n"
+ ]
+ }
+ ],
"source": [
"from qualtran.cirq_interop.t_complexity_protocol import t_complexity_compat\n",
"\n",
@@ -191,9 +246,19 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 16,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 7.77 s, sys: 92.5 ms, total: 7.87 s\n",
+ "Wall time: 7.94 s\n",
+ "cswap: 134219344, and_bloq: 4391438157, clifford: 20696810140, rotation: 33555056, measurement: 4391438157\n"
+ ]
+ }
+ ],
"source": [
"x_dim, y_dim = 20, 20\n",
"t = 20\n",
@@ -210,9 +275,19 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 17,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 7.76 s, sys: 118 ms, total: 7.88 s\n",
+ "Wall time: 7.96 s\n",
+ "t: 19, cswap: 134219344, and_bloq: 4391438389, clifford: 20696810716, rotation: 33555359, measurement: 4391438389\n"
+ ]
+ }
+ ],
"source": [
"# Or, we can just use the included bloq example directly\n",
"\n",
@@ -231,9 +306,678 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 18,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "T-count: 3.54214e+06\n",
+ "Rotations: 131428\n",
+ "Cliffords: 1.07402e+07\n",
+ "\n"
+ ]
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/fionn/.venvs/qualtran/lib/python3.10/site-packages/cotengra/hyperoptimizers/hyper.py:33: UserWarning: Couldn't import `kahypar` - skipping from default hyper optimizer and using basic `labels` method instead.\n",
+ " warnings.warn(\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
"source": [
"from qualtran.bloqs.phase_estimation.qubitization_qpe import _qubitization_qpe_hubbard_model_small\n",
"from qualtran.drawing import show_flame_graph\n",
@@ -242,6 +986,55 @@
"print(qpe_small.t_complexity())\n",
"show_flame_graph(qpe_small)"
]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Resource Estimates for THC Hamiltonian\n",
+ "\n",
+ "Phase estimation of the Reiher THC hamiltonian from [https://arxiv.org/abs/2011.03494](https://arxiv.org/abs/2011.03494)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 38,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 1.57 s, sys: 668 ms, total: 2.24 s\n",
+ "Wall time: 453 ms\n",
+ "result=11290\n",
+ "CPU times: user 7.9 s, sys: 125 ms, total: 8.03 s\n",
+ "Wall time: 8.07 s\n",
+ "result=2.26283042e+09\n"
+ ]
+ }
+ ],
+ "source": [
+ "from qualtran.bloqs.chemistry.thc.walk_operator import get_reiher_thc_walk_operator\n",
+ "from qualtran.bloqs.phase_estimation import LPResourceState, QubitizationQPE\n",
+ "from qualtran.resource_counting.generalizers import generalize_cswap_approx\n",
+ "\n",
+ "walk_op = get_reiher_thc_walk_operator()\n",
+ "%time result = get_cost_value(walk_op, QECGatesCost(), generalizer=generalize_cswap_approx).total_t_and_ccz_count(ts_per_rotation=0)['n_ccz']\n",
+ "print(f\"{result=}\")\n",
+ "algo_eps = 0.0016\n",
+ "qpe_eps = algo_eps / (walk_op.block_encoding.alpha * np.sqrt(2))\n",
+ "qpe = QubitizationQPE(\n",
+ " walk, LPResourceState.from_standard_deviation_eps(qpe_eps)\n",
+ ")\n",
+ "%time result = get_cost_value(qpe, QECGatesCost(), generalizer=generalize_cswap_approx).total_t_and_ccz_count(ts_per_rotation=0)['n_ccz']\n",
+ "print(f\"{result=:10.8e}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": []
}
],
"metadata": {
@@ -260,7 +1053,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.11.8"
+ "version": "3.10.4"
}
},
"nbformat": 4,