From 7cebbb082057c3247dd5d097badd282cb2d2d87e Mon Sep 17 00:00:00 2001
From: jaunatisblue <jaunatisblue@gmail.com>
Date: Sat, 9 May 2026 15:21:56 +0800
Subject: [PATCH] =?UTF-8?q?=E6=9B=B4=E6=96=B0baseline?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 baseline_mps_expectation.py | 39 ++++++++++++++++++++++++++++++++-----
 1 file changed, 34 insertions(+), 5 deletions(-)

diff --git a/baseline_mps_expectation.py b/baseline_mps_expectation.py
index 8a7e68e..a72a3f2 100644
--- a/baseline_mps_expectation.py
+++ b/baseline_mps_expectation.py
@@ -39,6 +39,24 @@ def build_observable(nqubits):
     return hamiltonians.SymbolicHamiltonian(form=form)
 
 
+def exact_expectation(circuit, nqubits):
+    import numpy as np
+
+    state = circuit().state(numpy=True).reshape(-1)
+    probabilities = np.abs(state) ** 2
+    indices = np.arange(state.size)
+
+    value = 0.0
+    for qubit in range(nqubits - 1):
+        left = (indices >> (nqubits - 1 - qubit)) & 1
+        right = (indices >> (nqubits - 2 - qubit)) & 1
+        value += 0.5 * np.sum(probabilities * (1 - 2 * left) * (1 - 2 * right))
+
+    flip_q0 = 1 << (nqubits - 1)
+    value += 0.25 * np.vdot(state[indices ^ flip_q0], state).real
+    return float(value)
+
+
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument("--nqubits", type=int, default=20)
@@ -47,7 +65,9 @@ def main():
     parser.add_argument("--seed", type=int, default=42)
     parser.add_argument("--cut-ratio", type=float, default=1e-12)
     parser.add_argument("--svd-control", default="V")
-    parser.add_argument("--no-exact", action="store_true")
+    parser.add_argument("--exact", action="store_true")
+    parser.add_argument("--exact-max-qubits", type=int, default=24)
+    parser.add_argument("--preprocess", action="store_true")
     args = parser.parse_args()
     logging.getLogger("qibo.config").setLevel(logging.ERROR)
     logging.getLogger("qtealeaves").setLevel(logging.ERROR)
@@ -55,10 +75,17 @@ def main():
     circuit = build_circuit(args.nqubits, args.nlayers, args.seed)
     observable = build_observable(args.nqubits)
     exact = None
-    if not args.no_exact:
-        exact = float(observable.expectation_from_state(circuit().state()).real)
+    if args.exact:
+        if args.nqubits > args.exact_max_qubits:
+            raise ValueError(
+                f"--exact is limited to {args.exact_max_qubits} qubits by default."
+            )
+        exact = exact_expectation(circuit, args.nqubits)
 
-    print(f"nqubits={args.nqubits} nlayers={args.nlayers} seed={args.seed}")
+    print(
+        f"nqubits={args.nqubits} nlayers={args.nlayers} "
+        f"seed={args.seed} preprocess={args.preprocess}"
+    )
     if exact is not None:
         print(f"exact={exact:.16e}")
     print("bond_dim expval abs_error rel_error seconds")
@@ -72,7 +99,9 @@ def main():
             svd_control=args.svd_control,
         )
         start = time.perf_counter()
-        value = float(backend.expectation(circuit, observable, preprocess=False).real)
+        value = float(
+            backend.expectation(circuit, observable, preprocess=args.preprocess).real
+        )
         elapsed = time.perf_counter() - start
         abs_error = float("nan") if exact is None else abs(value - exact)
         rel_error = float("nan") if exact is None else abs_error / max(abs(exact), 1e-15)