#!/usr/bin/env python3 """ qaoa_circuit.py — QAOA Pipeline for E=mc² Finsler Route Complete QAOA pipeline: 1. QUBO encoding of E=mc² path cost through Hachimoji state space 2. QUBO → Ising → Pauli string conversion 3. QAOA circuit construction (8 qubits, p=2 layers) via Cirq 4. Classical optimization of γ, β parameters 5. Measurement → bitstring → Hachimoji state decoding 6. Classical validation via HiGHS solver Research-Stack Trace: chaos_game → Finsler metric → QUBO → QAOA → Hachimoji state """ from __future__ import annotations import json import math import random import time import hashlib from collections import Counter from dataclasses import dataclass, field from pathlib import Path from typing import Any, Optional from datetime import datetime, timezone # --------------------------------------------------------------------------- # 0. Optional backend imports # --------------------------------------------------------------------------- try: import cirq _HAS_CIRQ = True except ImportError: _HAS_CIRQ = False print("[WARN] Cirq not installed — using fallback simulator") try: import highspy import numpy as np _HAS_HIGHS = True except ImportError: _HAS_HIGHS = False try: import scipy.optimize as _sp_opt _HAS_SCIPY = True except ImportError: _HAS_SCIPY = False try: import numpy as np _HAS_NUMPY = True except ImportError: _HAS_NUMPY = False # ========================================================================= # I. Data Models (from qaoa_adapter.py) # ========================================================================= @dataclass class QUBO: """Quadratic Unconstrained Binary Optimization. Minimize x^T Q x where x_i in {0, 1}. matrix[(i,j)] = coefficient for x_i * x_j (upper triangular). """ n: int matrix: dict[tuple[int, int], float] = field(default_factory=dict) offset: float = 0.0 def __post_init__(self) -> None: keys = list(self.matrix.keys()) for i, j in keys: if i > j: self.matrix[(j, i)] = self.matrix.pop((i, j)) def energy(self, x: list[int]) -> float: e = self.offset for (i, j), qij in self.matrix.items(): e += qij * x[i] * x[j] return e @dataclass class Ising: """Ising Hamiltonian: H = sum h_i s_i + sum J_ij s_i s_j + offset. s_i in {+1, -1}. J uses upper triangular (i < j). """ n: int h: list[float] = field(default_factory=list) J: dict[tuple[int, int], float] = field(default_factory=dict) offset: float = 0.0 def energy(self, s: list[int]) -> float: if len(s) != self.n: raise ValueError(f"expected {self.n} spins, got {len(s)}") e = self.offset for i in range(self.n): e += self.h[i] * s[i] for (i, j), Jij in self.J.items(): e += Jij * s[i] * s[j] return e @dataclass class PauliSum: """Pauli string representation of an Ising Hamiltonian. terms: list of (pauli_string, coefficient) e.g. ("ZZIIIIII", 0.5), ("ZIIIIIII", -1.0) offset: constant (identity) term """ n: int terms: list[tuple[str, float]] = field(default_factory=list) offset: float = 0.0 # ========================================================================= # II. Hachimoji State Constants # ========================================================================= # 8 Hachimoji states — indices 0-7 HACHIMOJI_STATES: list[str] = ["A", "T", "G", "C", "B", "S", "P", "Z"] HACHIMOJI_INDEX: dict[str, int] = {s: i for i, s in enumerate(HACHIMOJI_STATES)} # Semantic descriptions per state HACHIMOJI_DESC: dict[str, str] = { "A": "trivial — lowest energy basin", "T": "room — moderate energy", "G": "tight — low energy", "C": "marginal — moderate-high energy", "B": "collision — highest energy (avoid)", "S": "symmetric — high energy", "P": "potential — moderate energy", "Z": "zero — low-moderate energy", } # Greek mapping (from qaoa_adapter.py) GREEK_STATES: list[str] = ["\u03A6", "\u039B", "\u03A1", "\u039A", "\u03A9", "\u03A3", "\u03A0", "\u0396"] # ========================================================================= # III. QUBO → Ising → Pauli Conversion (from qaoa_adapter.py) # ========================================================================= def qubo_to_ising(qubo: QUBO) -> Ising: """Convert QUBO to Ising Hamiltonian. x = (1 + s) / 2 maps binary x in {0,1} to spin s in {+1,-1}. """ n = qubo.n a: dict[int, float] = {} b: dict[tuple[int, int], float] = {} for (i, j), qij in qubo.matrix.items(): if i == j: a[i] = a.get(i, 0.0) + qij else: key = (min(i, j), max(i, j)) b[key] = b.get(key, 0.0) + qij offset = qubo.offset for i in range(n): offset += 0.5 * a.get(i, 0.0) for (i, j), bij in b.items(): offset += 0.25 * bij h = [0.0] * n for i in range(n): hi = 0.5 * a.get(i, 0.0) for (j1, j2), bij in b.items(): if j1 == i: hi += 0.25 * bij if j2 == i: hi += 0.25 * bij h[i] = hi J: dict[tuple[int, int], float] = {} for (i, j), bij in b.items(): J[(i, j)] = 0.25 * bij return Ising(n=n, h=h, J=J, offset=offset) def ising_to_pauli(ising: Ising) -> PauliSum: """Convert Ising Hamiltonian to Pauli strings. Mapping: s_i → Z_i, s_i s_j → Z_i Z_j, offset → I (constant). """ terms: list[tuple[str, float]] = [] for i in range(ising.n): if abs(ising.h[i]) > 1e-15: ps = ["I"] * ising.n ps[i] = "Z" terms.append(("".join(ps), ising.h[i])) for (i, j), Jij in ising.J.items(): if abs(Jij) > 1e-15: ps = ["I"] * ising.n ps[i] = "Z" ps[j] = "Z" terms.append(("".join(ps), Jij)) return PauliSum(n=ising.n, terms=terms, offset=ising.offset) # ========================================================================= # IV. QAOA Circuit Construction (Cirq) # ========================================================================= def build_qaoa_circuit( pauli: PauliSum, p_layers: int, gamma: list[float], beta: list[float], ) -> Any: """Build a QAOA circuit from Pauli Hamiltonian. Structure: |+>^⊗n → [e^{-iγ_k H_C} e^{-iβ_k H_M}]_{k=1..p} → measure Returns cirq.Circuit if cirq available, else dict description. """ if not _HAS_CIRQ: return { "type": "qaoa_circuit_description", "n_qubits": pauli.n, "p_layers": p_layers, "gamma": gamma, "beta": beta, "cost_hamiltonian": [{"pauli": t[0], "coeff": t[1]} for t in pauli.terms], "offset": pauli.offset, "note": "Install cirq for runnable circuits", } n = pauli.n qubits = cirq.LineQubit.range(n) circuit = cirq.Circuit() # Initial Hadamards: prepare |+>^⊗n circuit.append(cirq.H.on_each(*qubits)) for layer in range(p_layers): g = gamma[layer] if layer < len(gamma) else gamma[-1] b = beta[layer] if layer < len(beta) else beta[-1] # Cost Hamiltonian: e^{-iγ H_C} for ps_str, coeff in pauli.terms: angle = 2.0 * g * coeff if abs(angle) < 1e-15: continue z_positions = [i for i, c in enumerate(ps_str) if c == "Z"] if len(z_positions) == 1: # Single Z: Rz(2γh_i) rotation circuit.append(cirq.rz(angle)(qubits[z_positions[0]])) elif len(z_positions) == 2: # ZZ interaction: use CZ + Rz decomposition i, j = z_positions # CZ**(angle/π) gives controlled phase of angle circuit.append(cirq.CZ(qubits[i], qubits[j]) ** (angle / math.pi)) circuit.append(cirq.rz(angle)(qubits[i])) circuit.append(cirq.rz(angle)(qubits[j])) elif len(z_positions) > 2: # Higher-order: decompose into individual rotations for idx in z_positions: circuit.append(cirq.rz(angle / len(z_positions))(qubits[idx])) # Mixer Hamiltonian: e^{-iβ H_M} where H_M = sum X_i circuit.append(cirq.rx(2.0 * b).on_each(*qubits)) # Measure in computational basis circuit.append(cirq.measure(*qubits, key="result")) return circuit def compute_circuit_depth(pauli: PauliSum, p_layers: int) -> int: """Compute theoretical circuit depth for the QAOA circuit. Gate model (one moment = gates that can run in parallel): - Initial H: 1 moment (all H gates parallel) - Per layer: * Single-Z rotations: 1 moment (all parallel on distinct qubits) * ZZ interactions: depends on overlap Non-overlapping pairs run in parallel; overlapping pairs sequential. Each ZZ needs: CZ + Rz(i) + Rz(j) = 3 gates, but CZ and Rz can partially overlap → ~2 moments per ZZ group * Rx mixer: 1 moment (all parallel) - Measurement: 1 moment """ n = pauli.n # Count single-Z and ZZ terms single_z = [t for t in pauli.terms if sum(1 for c in t[0] if c == "Z") == 1] zz_terms = [t for t in pauli.terms if sum(1 for c in t[0] if c == "Z") == 2] # Single-Z rotations: all run in parallel (different qubits) = 1 moment single_z_depth = 1 if single_z else 0 # ZZ interactions: compute overlapping groups zz_pairs = [] for ps_str, coeff in zz_terms: z_pos = [i for i, c in enumerate(ps_str) if c == "Z"] if len(z_pos) == 2: zz_pairs.append(tuple(z_pos)) if zz_pairs: # Greedy coloring: each "color" = parallel group # Two pairs can run in parallel if they share no qubits groups = [] for pair in zz_pairs: placed = False for group in groups: # Check if pair conflicts with any pair in group if not any(set(pair) & set(existing) for existing in group): group.append(pair) placed = True break if not placed: groups.append([pair]) # Each group: CZ + 2*Rz ≈ 2-3 moments zz_depth = len(groups) * 2 else: zz_depth = 0 per_layer = single_z_depth + zz_depth + 1 # +1 for Rx mixer total = 1 + p_layers * per_layer + 1 # H init + layers + measure return total def circuit_to_json(circuit: Any) -> dict: """Serialize a Cirq circuit to JSON-friendly dict.""" if _HAS_CIRQ and hasattr(circuit, "to_json"): return {"type": "cirq.Circuit", "moments": str(circuit)} elif isinstance(circuit, dict): return circuit else: return {"type": "circuit", "repr": str(circuit)} # ========================================================================= # V. Fallback Circuit Simulator (when Cirq unavailable) # ========================================================================= class FallbackSimulator: """Minimal statevector simulator for QAOA (no Cirq needed). Simulates the QAOA circuit by applying: - H gates (all qubits) - Rz(θ) rotations - CZ(θ) controlled phases - Rx(θ) rotations Uses explicit statevector evolution (2^n amplitudes). """ def __init__(self, n: int): self.n = n self.N = 1 << n # 2^n states def _apply_h(self, state: list[complex], q: int) -> list[complex]: """Apply Hadamard gate to qubit q.""" new_state = [0j] * self.N inv_sqrt2 = 1.0 / math.sqrt(2.0) for i in range(self.N): bit = (i >> q) & 1 j = i ^ (1 << q) # flip qubit q if bit == 0: new_state[i] += inv_sqrt2 * (state[i] + state[j]) new_state[j] += inv_sqrt2 * (state[i] - state[j]) # bit==1 handled by the bit==0 case of its partner return new_state def _apply_rz(self, state: list[complex], q: int, angle: float) -> list[complex]: """Apply Rz(angle) to qubit q.""" new_state = list(state) half_angle = angle / 2.0 for i in range(self.N): bit = (i >> q) & 1 if bit == 1: new_state[i] *= complex(math.cos(half_angle), -math.sin(half_angle)) else: new_state[i] *= complex(math.cos(half_angle), math.sin(half_angle)) return new_state def _apply_cz(self, state: list[complex], q1: int, q2: int, angle: float) -> list[complex]: """Apply CZ**(angle/π) = diag(1,1,1,e^{i*angle}).""" new_state = list(state) for i in range(self.N): b1 = (i >> q1) & 1 b2 = (i >> q2) & 1 if b1 == 1 and b2 == 1: new_state[i] *= complex(math.cos(angle), math.sin(angle)) return new_state def _apply_rx(self, state: list[complex], q: int, angle: float) -> list[complex]: """Apply Rx(angle) to qubit q.""" new_state = [0j] * self.N c = math.cos(angle / 2.0) s = -1j * math.sin(angle / 2.0) for i in range(self.N): bit = (i >> q) & 1 j = i ^ (1 << q) if bit == 0: new_state[i] += c * state[i] + s * state[j] else: # bit==1: handled when we see the flipped partner new_state[i] += s.conjugate() * state[j] + c * state[i] return new_state def run(self, pauli: PauliSum, gamma: list[float], beta: list[float], shots: int = 1000) -> dict[str, int]: """Run QAOA and return measurement counts.""" # Initialize |0...0> state = [0j] * self.N state[0] = 1.0 # Apply H^⊗n to get |+...+ for q in range(self.n): state = self._apply_h(state, q) p = len(gamma) for layer in range(p): g = gamma[layer] b = beta[layer] # Cost layer for ps_str, coeff in pauli.terms: angle = 2.0 * g * coeff if abs(angle) < 1e-15: continue z_positions = [i for i, c in enumerate(ps_str) if c == "Z"] if len(z_positions) == 1: state = self._apply_rz(state, z_positions[0], angle) elif len(z_positions) == 2: i, j = z_positions state = self._apply_cz(state, i, j, angle) state = self._apply_rz(state, i, angle) state = self._apply_rz(state, j, angle) elif len(z_positions) > 2: for idx in z_positions: state = self._apply_rz(state, idx, angle / len(z_positions)) # Mixer layer for q in range(self.n): state = self._apply_rx(state, q, 2.0 * b) # Compute probabilities probs = [abs(a)**2 for a in state] total_prob = sum(probs) if total_prob > 0: probs = [p / total_prob for p in probs] # Sample shots counts: dict[str, int] = {} for _ in range(shots): r = random.random() cumsum = 0.0 for i in range(self.N): cumsum += probs[i] if r <= cumsum: bits = format(i, f"0{self.n}b") counts[bits] = counts.get(bits, 0) + 1 break return counts # ========================================================================= # VI. Classical Optimization Loop # ========================================================================= def qaoa_expectation( params: list[float], pauli: PauliSum, qubo: QUBO, p_layers: int, shots: int = 1000, ) -> float: """Compute <ψ|H_C|ψ> for given γ, β parameters. params = [γ_1, ..., γ_p, β_1, ..., β_p] Returns the QUBO energy of the most probable measurement outcome. """ gamma = params[:p_layers] beta = params[p_layers:] if _HAS_CIRQ: circuit = build_qaoa_circuit(pauli, p_layers, gamma, beta) simulator = cirq.Simulator() samples = simulator.run(circuit, repetitions=shots) counts_int = samples.histogram(key="result") # Convert int keys to bitstring keys counts: dict[str, int] = {} for val, cnt in counts_int.items(): bits = format(val, f"0{qubo.n}b") counts[bits] = cnt else: sim = FallbackSimulator(qubo.n) counts = sim.run(pauli, gamma, beta, shots=shots) # Compute expectation value: weighted average energy total = sum(counts.values()) expectation = 0.0 for bits, count in counts.items(): x = [int(b) for b in bits] energy = qubo.energy(x) expectation += (count / total) * energy return expectation def optimize_qaoa_params( pauli: PauliSum, qubo: QUBO, p_layers: int = 2, shots: int = 2000, maxiter: int = 100, seed: int = 42, ) -> dict: """Optimize γ, β to minimize <ψ|H_C|ψ>. Uses COBYLA (gradient-free) via scipy if available, otherwise simple grid search + random restart. """ random.seed(seed) t0 = time.time() if _HAS_SCIPY: # Use COBYLA for gradient-free optimization def objective(params): return qaoa_expectation(params, pauli, qubo, p_layers, shots) # Initial guess: spread γ across [0.1, 1.0], β across [0.1, 0.5] x0 = [] for k in range(p_layers): x0.append(0.5 + 0.3 * (k / max(p_layers, 1))) # γ for k in range(p_layers): x0.append(0.3 - 0.1 * (k / max(p_layers, 1))) # β # Bounds-like constraints via penalties not needed for COBYLA result = _sp_opt.minimize( objective, x0, method="COBYLA", options={"maxiter": maxiter, "rhobeg": 0.1}, ) opt_gamma = result.x[:p_layers].tolist() opt_beta = result.x[p_layers:].tolist() best_energy = result.fun n_iters = result.nfev status = result.message else: # Fallback: grid search + random perturbation best_energy = float("inf") opt_gamma = [0.5] * p_layers opt_beta = [0.3] * p_layers gamma_grid = [0.1, 0.3, 0.5, 0.7, 1.0, 1.5] beta_grid = [0.1, 0.2, 0.3, 0.4, 0.5] n_iters = 0 for g1 in gamma_grid: for g2 in gamma_grid: for b1 in beta_grid: for b2 in beta_grid: params = [g1, g2, b1, b2] if p_layers == 2 else [g1, b1] energy = qaoa_expectation( params, pauli, qubo, p_layers, shots ) n_iters += 1 if energy < best_energy: best_energy = energy opt_gamma = params[:p_layers] opt_beta = params[p_layers:] # Random refinement for _ in range(50): rg = [random.uniform(0.05, 2.0) for _ in range(p_layers)] rb = [random.uniform(0.05, 0.8) for _ in range(p_layers)] energy = qaoa_expectation(rg + rb, pauli, qubo, p_layers, shots) n_iters += 1 if energy < best_energy: best_energy = energy opt_gamma = rg opt_beta = rb status = "grid_search_fallback" runtime = time.time() - t0 return { "gamma": opt_gamma, "beta": opt_beta, "best_energy": best_energy, "iterations": n_iters, "runtime_s": round(runtime, 4), "status": status, } # ========================================================================= # VII. Measurement → Bitstring → Hachimoji State # ========================================================================= def measurements_to_solution(counts: dict[str, int], n: int) -> list[int]: """Extract most probable bitstring from measurement counts.""" if not counts: return [0] * n best_bitstring = max(counts, key=counts.get) if len(best_bitstring) < n: best_bitstring = best_bitstring.zfill(n) return [int(b) for b in best_bitstring[-n:]] def bitstring_to_hachimoji(bits: list[int]) -> dict: """Map a binary solution to Hachimoji state interpretation. For this QUBO (selecting the "center" state), the optimal is typically a SINGLE 1 (one state selected) due to the diagonal dominance. The state with x_i=1 is the chosen Hachimoji center. """ active_states = [] for i, b in enumerate(bits): if b == 1: active_states.append({ "index": i, "state": HACHIMOJI_STATES[i], "greek": GREEK_STATES[i], "description": HACHIMOJI_DESC[HACHIMOJI_STATES[i]], }) # Determine primary state (if exactly one active, that's it; # if multiple, pick the one with lowest diagonal energy; # if none, all states are "off" → report A as default) if len(active_states) == 1: primary = active_states[0] elif len(active_states) > 1: # In multi-active case, prefer lower-energy states # Sort by diagonal QUBO energy (lower = better) diag_energies = { 0: -1.0, 1: 0.5, 2: 0.3, 3: 0.8, 4: 2.0, 5: 1.5, 6: 1.0, 7: 0.7, } active_states.sort(key=lambda s: diag_energies.get(s["index"], 0)) primary = active_states[0] else: primary = { "index": 0, "state": "A", "greek": "\u03A6", "description": HACHIMOJI_DESC["A"], } active_states = [primary] return { "primary_state": primary, "active_states": active_states, "bitstring": "".join(str(b) for b in bits), "hamming_weight": sum(bits), } # ========================================================================= # VIII. Classical HiGHS Solver Validation # ========================================================================= def solve_qubo_highs(qubo: QUBO, time_limit: float = 10.0) -> dict: """Solve QUBO using HiGHS MIP solver (or SA fallback). Wraps qubo_highs.py solve_qubo_highs if available, otherwise uses built-in simulated annealing. """ Q_dict = dict(qubo.matrix) n = qubo.n t0 = time.time() if _HAS_HIGHS: try: return _solve_highs_mip(Q_dict, n, time_limit) except Exception as exc: print(f"[WARN] HiGHS failed: {exc}, falling back to SA") return _solve_sa(Q_dict, n, time_limit) else: print("[INFO] highspy not installed, using SA") return _solve_sa(Q_dict, n, time_limit) def _solve_highs_mip(Q: dict, n: int, time_limit: float) -> dict: """Solve QUBO via HiGHS MIP linearization.""" # QUBO: min x^T Q x = sum_i Q[i,i]*x_i + sum_{i= x[i]+x[j]-1 col_entries: dict[int, list[tuple[int, float]]] = {} row_idx = 0 for (i, j), yi in y_index_map.items(): # y <= x[i] if yi not in col_entries: col_entries[yi] = [] col_entries[yi].append((row_idx, 1.0)) if i not in col_entries: col_entries[i] = [] col_entries[i].append((row_idx, -1.0)) row_idx += 1 # y <= x[j] if yi not in col_entries: col_entries[yi] = [] col_entries[yi].append((row_idx, 1.0)) if j not in col_entries: col_entries[j] = [] col_entries[j].append((row_idx, -1.0)) row_idx += 1 # y >= x[i] + x[j] - 1 => -y + x[i] + x[j] <= 1 if yi not in col_entries: col_entries[yi] = [] col_entries[yi].append((row_idx, -1.0)) if i not in col_entries: col_entries[i] = [] col_entries[i].append((row_idx, 1.0)) if j not in col_entries: col_entries[j] = [] col_entries[j].append((row_idx, 1.0)) row_idx += 1 starts = [] indices_list = [] values_list = [] nnz = 0 for col in range(num_vars): starts.append(nnz) if col in col_entries: for ri, val in col_entries[col]: indices_list.append(ri) values_list.append(val) nnz += 1 starts.append(nnz) lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise lp.a_matrix_.start_ = np.array(starts, dtype=np.int32) lp.a_matrix_.index_ = np.array(indices_list, dtype=np.int32) lp.a_matrix_.value_ = np.array(values_list) # Row bounds row_lower = [] row_upper = [] for _ in range(len(y_index_map)): row_lower.append(-1e30) row_upper.append(0.0) row_lower.append(-1e30) row_upper.append(0.0) row_lower.append(-1e30) row_upper.append(1.0) lp.row_lower_ = np.array(row_lower) lp.row_upper_ = np.array(row_upper) h = highspy.Highs() h.setOptionValue("time_limit", time_limit) h.setOptionValue("output_flag", False) h.passModel(model) h.run() sol = h.getSolution() x_vals = sol.col_value solution = [int(round(x_vals[i])) for i in range(n)] objective = sum(Q.get((i, j), 0.0) * solution[i] * solution[j] for i in range(n) for j in range(n) if (i, j) in Q) status_val = h.getInfoValue("primal_solution_status")[1] status_map = {0: "unknown", 1: "infeasible", 2: "feasible", 3: "optimal"} status_str = status_map.get(status_val, f"status_{status_val}") return { "solution": solution, "objective": objective, "status": status_str, "runtime_s": round(time.time() - t0, 4), } def _solve_sa(Q: dict, n: int, time_limit: float, seed: int = 42) -> dict: """Simulated annealing fallback.""" random.seed(seed) t0 = time.time() x = [random.randint(0, 1) for _ in range(n)] def energy(xv): return sum(Q.get((i, j), 0.0) * xv[i] * xv[j] for i in range(n) for j in range(n) if (i, j) in Q) current_energy = energy(x) best_x = list(x) best_energy = current_energy T = 10.0 iterations = 0 while time.time() - t0 < time_limit and T > 1e-8: i = random.randint(0, n - 1) x[i] = 1 - x[i] new_energy = energy(x) delta = new_energy - current_energy if delta < 0 or random.random() < math.exp(-delta / max(T, 1e-30)): current_energy = new_energy if current_energy < best_energy: best_x = list(x) best_energy = current_energy else: x[i] = 1 - x[i] T *= 0.9995 iterations += 1 return { "solution": best_x, "objective": best_energy, "status": "SA_heuristic", "iterations": iterations, "runtime_s": round(time.time() - t0, 4), } # ========================================================================= # IX. Brute-Force Exact Solver (for small n, validates everything) # ========================================================================= def solve_qubo_brute_force(qubo: QUBO) -> dict: """Brute-force over all 2^n assignments. Exact for n <= 20.""" n = qubo.n best_energy = float("inf") best_solution = [0] * n all_energies: dict[str, float] = {} for val in range(1 << n): x = [(val >> i) & 1 for i in range(n)] e = qubo.energy(x) bits = "".join(str(b) for b in x) all_energies[bits] = e if e < best_energy: best_energy = e best_solution = list(x) # Count degeneracy degeneracy = sum(1 for e in all_energies.values() if abs(e - best_energy) < 1e-10) return { "solution": best_solution, "objective": best_energy, "degeneracy": degeneracy, "status": "exact", "total_evaluated": 2 ** n, } # ========================================================================= # X. Main Pipeline: QAOA for E=mc² # ========================================================================= def build_emc2_qubo() -> QUBO: """Build the E=mc² QUBO encoding the Finsler path cost. Variable x_i = 1 if the equation's center is in Hachimoji state i. Diagonal Q_ii = self-energy (Finsler α component). Off-diagonal Q_ij = interaction energy (Finsler β drift). """ Q = { (0, 0): -1.0, # A (trivial) — lowest self-energy (1, 1): 0.5, # T (room) (2, 2): 0.3, # G (tight) (3, 3): 0.8, # C (marginal) (4, 4): 2.0, # B (collision) — highest self-energy (5, 5): 1.5, # S (symmetric) (6, 6): 1.0, # P (potential) (7, 7): 0.7, # Z (zero) # Off-diagonal: coupling between states (0, 1): -0.3, # A favors T (0, 2): -0.2, # A favors G (4, 5): 0.5, # B repels S (4, 6): 0.5, # B repels P } return QUBO(n=8, matrix=Q, offset=0.0) def run_qaoa_pipeline( p_layers: int = 2, shots: int = 2000, optimize_maxiter: int = 100, seed: int = 42, ) -> dict: """Run the complete QAOA pipeline for E=mc². Returns comprehensive result dict with all stages documented. """ print("=" * 70) print("QAOA Pipeline: E=mc² → Hachimoji State") print("=" * 70) overall_t0 = time.time() # --- Stage 1: Build QUBO --- print("\n[Stage 1] Building QUBO for E=mc²...") qubo = build_emc2_qubo() print(f" QUBO: n={qubo.n} variables, {len(qubo.matrix)} terms") print(f" Hachimoji states: {HACHIMOJI_STATES}") for (i, j), val in sorted(qubo.matrix.items()): si = HACHIMOJI_STATES[i] sj = HACHIMOJI_STATES[j] if j != i else "" label = f"{si}{sj}" if i != j else f"{si} (diag)" print(f" Q[{label}] = {val:+.2f}") # --- Stage 2: QUBO → Ising --- print("\n[Stage 2] Converting QUBO → Ising Hamiltonian...") ising = qubo_to_ising(qubo) print(f" Ising: n={ising.n}, offset={ising.offset:.4f}") h_str = ", ".join(f"h[{HACHIMOJI_STATES[i]}]={ising.h[i]:+.4f}" for i in range(ising.n)) print(f" Linear terms: {h_str}") j_strs = [f"J[{HACHIMOJI_STATES[i]}{HACHIMOJI_STATES[j]}]={Jij:+.4f}" for (i, j), Jij in ising.J.items()] if j_strs: print(f" Coupling terms: {', '.join(j_strs)}") # --- Stage 3: Ising → Pauli --- print("\n[Stage 3] Converting Ising → Pauli strings...") pauli = ising_to_pauli(ising) print(f" Pauli sum: {len(pauli.terms)} non-zero terms, offset={pauli.offset:.4f}") for ps, coeff in pauli.terms: z_pos = [i for i, c in enumerate(ps) if c == "Z"] label = " * ".join(f"Z_{HACHIMOJI_STATES[i]}" for i in z_pos) print(f" {label}: {coeff:+.4f}") # --- Stage 4: Optimize γ, β --- print(f"\n[Stage 4] Optimizing QAOA parameters (p={p_layers})...") print(f" Classical optimizer: {'COBYLA (scipy)' if _HAS_SCIPY else 'grid search'}") print(f" Shots per evaluation: {shots}") print(f" Max iterations: {optimize_maxiter}") opt_result = optimize_qaoa_params( pauli, qubo, p_layers=p_layers, shots=shots, maxiter=optimize_maxiter, seed=seed ) opt_gamma = opt_result["gamma"] opt_beta = opt_result["beta"] print(f" Optimal γ = {[round(g, 4) for g in opt_gamma]}") print(f" Optimal β = {[round(b, 4) for b in opt_beta]}") print(f" Best <ψ|H_C|ψ> = {opt_result['best_energy']:.4f}") print(f" Optimization iterations: {opt_result['iterations']}") print(f" Optimization runtime: {opt_result['runtime_s']:.2f}s") # --- Stage 5: Build final circuit and measure --- print(f"\n[Stage 5] Running QAOA circuit (shots={shots})...") # Compute circuit depth (works for both Cirq and fallback) circuit_depth = compute_circuit_depth(pauli, p_layers) if _HAS_CIRQ: circuit = build_qaoa_circuit(pauli, p_layers, opt_gamma, opt_beta) simulator = cirq.Simulator() samples = simulator.run(circuit, repetitions=shots) counts_int = samples.histogram(key="result") counts: dict[str, int] = {} for val, cnt in counts_int.items(): bits = format(val, f"0{qubo.n}b") counts[bits] = cnt circuit_moments = str(circuit) else: sim = FallbackSimulator(qubo.n) counts = sim.run(pauli, opt_gamma, opt_beta, shots=shots) circuit = build_qaoa_circuit(pauli, p_layers, opt_gamma, opt_beta) circuit_moments = str(circuit) if isinstance(circuit, dict) else str(circuit) # Top 5 measurement outcomes sorted_counts = sorted(counts.items(), key=lambda kv: -kv[1])[:5] print(f" Top 5 measurement outcomes:") for bits, cnt in sorted_counts: x = [int(b) for b in bits] e = qubo.energy(x) active = [HACHIMOJI_STATES[i] for i, b in enumerate(x) if b == 1] states_str = ",".join(active) if active else "(none)" print(f" |{bits}> (states={states_str}): count={cnt}, energy={e:.4f}") # Best solution best_bits = measurements_to_solution(counts, qubo.n) best_energy = qubo.energy(best_bits) best_bitstring = "".join(str(b) for b in best_bits) # --- Stage 6: Decode to Hachimoji state --- print(f"\n[Stage 6] Decoding measurement to Hachimoji state...") hachi = bitstring_to_hachimoji(best_bits) primary = hachi["primary_state"] print(f" Most probable bitstring: |{hachi['bitstring']}>") print(f" Hamming weight: {hachi['hamming_weight']}") print(f" Primary Hachimoji state: {primary['state']} ({primary['greek']})") print(f" Description: {primary['description']}") if len(hachi["active_states"]) > 1: print(f" All active states:") for s in hachi["active_states"]: print(f" - {s['state']} ({s['greek']}): {s['description']}") # --- Stage 7: Classical validation --- print(f"\n[Stage 7] Classical validation...") # Exact brute-force (n=8 → 256 evaluations) print(f" Brute-force exact solution:") exact = solve_qubo_brute_force(qubo) exact_bits = exact["solution"] exact_hachi = bitstring_to_hachimoji(exact_bits) print(f" Exact optimum: |{''.join(str(b) for b in exact_bits)}>") print(f" Exact energy: {exact['objective']:.4f}") print(f" Exact state: {exact_hachi['primary_state']['state']}") print(f" Degeneracy: {exact['degeneracy']} solutions") # HiGHS / SA solver print(f" HiGHS/SA solver:") classical = solve_qubo_highs(qubo, time_limit=5.0) classical_bits = classical["solution"] classical_hachi = bitstring_to_hachimoji(classical_bits) print(f" Solver solution: |{''.join(str(b) for b in classical_bits)}>") print(f" Solver energy: {classical['objective']:.4f}") print(f" Solver state: {classical_hachi['primary_state']['state']}") print(f" Solver status: {classical['status']}") # --- Summary --- total_runtime = time.time() - overall_t0 print(f"\n{'=' * 70}") print(f"RESULT SUMMARY") print(f"{'=' * 70}") print(f" Equation: E = mc²") print(f" QUBO variables: 8 (Hachimoji states)") print(f" QAOA depth (p): {p_layers}") print(f" Circuit depth: {circuit_depth}") print(f" Optimal γ: {[round(g, 4) for g in opt_gamma]}") print(f" Optimal β: {[round(b, 4) for b in opt_beta]}") print(f"") print(f" QAOA best state: {primary['state']} (Greek: {primary['greek']})") print(f" QAOA best energy: {best_energy:.4f}") print(f" Exact best state: {exact_hachi['primary_state']['state']}") print(f" Exact best energy: {exact['objective']:.4f}") print(f" Classical state: {classical_hachi['primary_state']['state']}") print(f" Classical energy: {classical['objective']:.4f}") print(f"") print(f" Agreement QAOA==Exact: {primary['state'] == exact_hachi['primary_state']['state']}") print(f" Agreement Classical==Exact: {classical_hachi['primary_state']['state'] == exact_hachi['primary_state']['state']}") print(f" Total pipeline runtime: {total_runtime:.2f}s") print(f"{'=' * 70}") # Build comprehensive receipt receipt = { "schema": "rrc_qaoa_emc2_v1", "computed_at": datetime.now(timezone.utc).isoformat(), "equation": "E = mc^2", "qaoa_config": { "n_qubits": 8, "p_layers": p_layers, "shots": shots, "optimize_maxiter": optimize_maxiter, "backend": "cirq" if _HAS_CIRQ else "fallback_statevector", }, "qubo": { "n": qubo.n, "matrix": {f"({i},{j})": v for (i, j), v in qubo.matrix.items()}, "offset": qubo.offset, }, "ising": { "h": {HACHIMOJI_STATES[i]: ising.h[i] for i in range(ising.n)}, "J": {f"({HACHIMOJI_STATES[i]},{HACHIMOJI_STATES[j]})": v for (i, j), v in ising.J.items()}, "offset": ising.offset, }, "pauli_terms": [ {"pauli": t[0], "coefficient": t[1], "description": " * ".join( f"Z_{HACHIMOJI_STATES[i]}" for i, c in enumerate(t[0]) if c == "Z" )} for t in pauli.terms ], "circuit": { "depth": circuit_depth, "n_qubits": qubo.n, "p_layers": p_layers, "gamma": opt_gamma, "beta": opt_beta, }, "optimization": { "gamma": opt_gamma, "beta": opt_beta, "best_expectation": opt_result["best_energy"], "iterations": opt_result["iterations"], "runtime_s": opt_result["runtime_s"], "status": opt_result["status"], }, "measurement": { "top_outcomes": [ { "bitstring": bits, "count": cnt, "probability": round(cnt / shots, 4), "energy": round(qubo.energy([int(b) for b in bits]), 4), "active_states": [HACHIMOJI_STATES[i] for i, b in enumerate(bits) if b == "1"], } for bits, cnt in sorted_counts ], "most_probable": { "bitstring": best_bitstring, "energy": round(best_energy, 4), }, }, "hachimoji_result": { "primary_state": primary["state"], "primary_greek": primary["greek"], "description": primary["description"], "active_states": [s["state"] for s in hachi["active_states"]], "hamming_weight": hachi["hamming_weight"], }, "validation": { "exact": { "solution": exact_bits, "objective": exact["objective"], "degeneracy": exact["degeneracy"], "primary_state": exact_hachi["primary_state"]["state"], "status": exact["status"], }, "classical_solver": { "solution": classical_bits, "objective": classical["objective"], "primary_state": classical_hachi["primary_state"]["state"], "status": classical["status"], "runtime_s": classical.get("runtime_s", 0), }, }, "agreement": { "qaoa_matches_exact": primary["state"] == exact_hachi["primary_state"]["state"], "classical_matches_exact": classical_hachi["primary_state"]["state"] == exact_hachi["primary_state"]["state"], }, "runtime": { "optimization_s": opt_result["runtime_s"], "total_s": round(total_runtime, 4), }, } return receipt # ========================================================================= # XI. CLI # ========================================================================= def main(): import argparse parser = argparse.ArgumentParser(description="QAOA Pipeline for E=mc² Finsler Route") parser.add_argument("--p-layers", type=int, default=2, help="QAOA depth (default: 2)") parser.add_argument("--shots", type=int, default=2000, help="Measurement shots (default: 2000)") parser.add_argument("--maxiter", type=int, default=100, help="Optimizer max iterations (default: 100)") parser.add_argument("--seed", type=int, default=42, help="Random seed (default: 42)") parser.add_argument("--output", "-o", type=str, default="/mnt/agents/output/e2e/qaoa_receipt.json", help="Output receipt path") parser.add_argument("--quiet", action="store_true", help="Suppress console output") args = parser.parse_args() if args.quiet: import sys old_stdout = sys.stdout sys.stdout = open("/dev/null", "w") try: receipt = run_qaoa_pipeline( p_layers=args.p_layers, shots=args.shots, optimize_maxiter=args.maxiter, seed=args.seed, ) finally: if args.quiet: sys.stdout.close() sys.stdout = old_stdout # Write receipt output_path = Path(args.output) output_path.parent.mkdir(parents=True, exist_ok=True) with open(output_path, "w") as f: json.dump(receipt, f, indent=2, default=str) # Print key results (even in quiet mode, print these) if not args.quiet: print(f"\nReceipt written to: {output_path}") # Always print the concise result line primary = receipt["hachimoji_result"]["primary_state"] depth = receipt["circuit"]["depth"] print(f"\nRESULT: Hachimoji state = {primary}, circuit depth = {depth}") return receipt if __name__ == "__main__": main()