From 6180d6c1bba00c92fdad7789e73a4e0b478cf389 Mon Sep 17 00:00:00 2001 From: Anton Date: Fri, 11 Sep 2026 09:36:28 +0200 Subject: [PATCH 1/2] This is a single-qubit research demo. The objective of the demo is to make a link between p-bit probabilistic computing, Riemannian geometry, and the quantum Bloch sphere. It explores this combination for a potential future quantum simulator. A pure single-qubit state is treated first as a point on the Riemannian manifold CP^1, represented by the Bloch sphere S^2 with the Fubini-Study metric. Three independent p-bits statistically realize the Bloch coordinates, gates act as CP^1 isometries, and one additional p-bit samples projective measurements. --- examples/quantum_manifold_pbit.py | 576 ++++++++++++++++++++++++++++++ 1 file changed, 576 insertions(+) create mode 100644 examples/quantum_manifold_pbit.py diff --git a/examples/quantum_manifold_pbit.py b/examples/quantum_manifold_pbit.py new file mode 100644 index 0000000..dc1560c --- /dev/null +++ b/examples/quantum_manifold_pbit.py @@ -0,0 +1,576 @@ +""" + +This is a single-qubit research demo. The objective of the demo is to make +a link between p-bit probabilistic computing, Riemannian geometry, and the +quantum Bloch sphere. It explores this combination for a potential future +quantum simulator. + +A pure single-qubit state is treated first as a point on the Riemannian +manifold CP^1, represented by the Bloch sphere S^2 with the Fubini-Study +metric. Three independent p-bits statistically realize the Bloch +coordinates, gates act as CP^1 isometries, and one additional p-bit samples +projective measurements. + +So the the qubit lives on CP^1, p-bits realize that state statistically, +gates move it geometrically, and measurements are sampled probabilistically. + +GATE_SEQUENCE represents a small quantum circuit test. + +ManifoldPBitQubit # single-qubit simulator object +├── CP1BlochManifold # geometry +└── PBitCP1Backend # realization using p-bits + +The demo reports: + * the raw p-bit Bloch-vector norm before projection back to CP^1; + * an effective pre-projection purity indicator (1 + |raw_means|^2) / 2; + * Fubini-Study reconstruction error and target-state fidelity; + * tangent-space error decomposition into local polar/azimuthal components; + * mean and standard deviation over multiple independent p-bit trials; + * Born-rule measurement counts and projective-collapse behavior. + +The raw norm/purity values are diagnostics of the stochastic state estimate; +finite sampling and finite p-bit bias can move the raw estimate away from the +unit sphere, so they should not by themselves be interpreted as physical +mixed-state tomography. + +Results: + The quantum manifold p-bit representation is numerically stable and reproduces + single-qubit circuit behavior well. + +The question is how far we can go from here. + +Could coupled p-bits provide a stochastic representation of joint probability +distributions and correlations associated with points on higher-dimensional +quantum-state manifolds? + +""" + +import numpy as np + +from p_kit.psl import PCircuit +from p_kit.solver.csd_solver import CaSuDaSolver + + +# --------------------------------------------------------------------------- +# p-bit numerical parameters +# --------------------------------------------------------------------------- + +BIAS_CLIP = 0.995 + +STATE_NT = 12_000 +STATE_BURN_IN = 2_000 +STATE_SHOTS = 32 + +MEASURE_NT = 1_000 + +DT = 0.02 + +N_TRIALS = 8 # multi-seed statistics; each trial re-runs the whole + # 4-gate circuit, so cost scales linearly with this. +BASE_SEED = 42 + +GATE_SEQUENCE = [ # quantum circuit + ("h", ()), + ("t", ()), + ("ry", (np.pi / 3.0,)), + ("rz", (-np.pi / 5.0,)), +] + + +# --------------------------------------------------------------------------- +# Riemannian manifold: CP^1 represented by the unit Bloch sphere +# --------------------------------------------------------------------------- + +class CP1BlochManifold: + """ + Pure single-qubit state manifold CP^1 represented as S^2. + Riemannian geometry is implemented cleanly with NumPy here. + """ + + X_AXIS = np.array([1.0, 0.0, 0.0]) + Y_AXIS = np.array([0.0, 1.0, 0.0]) + Z_AXIS = np.array([0.0, 0.0, 1.0]) + + @staticmethod + def point(r): + r = np.asarray(r, dtype=float) + norm = np.linalg.norm(r) + if norm == 0.0: + raise ValueError("The zero vector cannot define a pure qubit state.") + return r / norm + + @classmethod + def from_angles(cls, theta, phi): + return cls.point( + [ + np.sin(theta) * np.cos(phi), + np.sin(theta) * np.sin(phi), + np.cos(theta), + ] + ) + + @classmethod + def angles(cls, point): + r = cls.point(point) + theta = float(np.arccos(np.clip(r[2], -1.0, 1.0))) + phi = float(np.arctan2(r[1], r[0])) + return theta, phi + + @classmethod + def tangent(cls, point, ambient_vector): + p = cls.point(point) + v = np.asarray(ambient_vector, dtype=float) + return v - np.dot(v, p) * p + + @classmethod + def exp_map(cls, point, tangent_vector): + p = cls.point(point) + v = cls.tangent(p, tangent_vector) + theta = np.linalg.norm(v) + if theta < 1e-15: + return p.copy() + return cls.point(np.cos(theta) * p + np.sin(theta) * (v / theta)) + + @classmethod + def log_map(cls, point_a, point_b): + p = cls.point(point_a) + q = cls.point(point_b) + dot = float(np.clip(np.dot(p, q), -1.0, 1.0)) + angle = float(np.arccos(dot)) + if angle < 1e-15: + return np.zeros(3) + if np.pi - angle < 1e-12: + raise ValueError("Log map is not unique for antipodal states.") + direction = q - dot * p + direction /= np.linalg.norm(direction) + return angle * direction + + @classmethod + def geodesic(cls, point_a, point_b, t): + if not 0.0 <= t <= 1.0: + raise ValueError("t must lie in [0, 1].") + p = cls.point(point_a) + return cls.exp_map(p, t * cls.log_map(p, point_b)) + + @classmethod + def distance(cls, point_a, point_b): + p = cls.point(point_a) + q = cls.point(point_b) + dot = float(np.clip(np.dot(p, q), -1.0, 1.0)) + return 0.5 * float(np.arccos(dot)) + + @classmethod + def fidelity(cls, point_a, point_b): + p = cls.point(point_a) + q = cls.point(point_b) + return float(np.clip(0.5 * (1.0 + np.dot(p, q)), 0.0, 1.0)) + + @classmethod + def rotate(cls, point, axis, angle): + r = cls.point(point) + n = cls.point(axis) + c = np.cos(angle) + s = np.sin(angle) + rotated = c * r + s * np.cross(n, r) + (1.0 - c) * np.dot(n, r) * n + return cls.point(rotated) + + @classmethod + def born_expectation(cls, point, axis): + r = cls.point(point) + n = cls.point(axis) + return float(np.clip(np.dot(r, n), -1.0, 1.0)) + + @classmethod + def born_probability(cls, point, axis): + expectation = cls.born_expectation(point, axis) + return 0.5 * (1.0 + expectation) + +def tangent_frame(theta, phi): + """Orthonormal tangent basis (e_theta, e_phi) at Bloch angles (theta, phi).""" + e_theta = np.array( + [ + np.cos(theta) * np.cos(phi), + np.cos(theta) * np.sin(phi), + -np.sin(theta), + ] + ) + e_phi = np.array([-np.sin(phi), np.cos(phi), 0.0]) + return e_theta, e_phi + +def tangent_decompose(manifold, base_point, other_point): + """ + Decompose the Riemannian reconstruction error into the local + (e_theta, e_phi) tangent basis. + + Returns (d_theta, d_phi, total_angle), where total_angle is the + Bloch-sphere tangent norm and equals 2 * d_FS. + + At the Bloch poles, spherical tangent components are undefined, + so d_theta and d_phi are returned as NaN. + """ + tangent = manifold.log_map(base_point, other_point) + norm = np.linalg.norm(tangent) + + if norm < 1e-12: + return 0.0, 0.0, 0.0 + + theta, phi = manifold.angles(base_point) + + # Spherical-coordinate tangent directions are undefined at the poles. + if abs(np.sin(theta)) < 1e-10: + return np.nan, np.nan, float(norm) + + e_theta, e_phi = tangent_frame(theta, phi) + + return ( + float(np.dot(tangent, e_theta)), + float(np.dot(tangent, e_phi)), + float(norm), + ) + +# --------------------------------------------------------------------------- +# p-kit realization and measurement backend +# --------------------------------------------------------------------------- + +class PBitCP1Backend: + """Realize CP^1 points using p-bit expectation values and sample measurements.""" + + def __init__( + self, + bias_clip=BIAS_CLIP, + state_nt=STATE_NT, + state_burn_in=STATE_BURN_IN, + state_shots=STATE_SHOTS, + measure_nt=MEASURE_NT, + dt=DT, + seed=1234, + ): + self.bias_clip = float(bias_clip) + self.state_nt = int(state_nt) + self.state_burn_in = int(state_burn_in) + self.state_shots = int(state_shots) + self.measure_nt = int(measure_nt) + self.dt = float(dt) + self.rng = np.random.default_rng(seed) + + def _next_seed(self): + return int(self.rng.integers(0, 2**31 - 1)) + + def _build_independent_circuit(self, target_means): + target_means = np.asarray(target_means, dtype=float).reshape(-1) + circuit = PCircuit(len(target_means)) + circuit.J = np.zeros((len(target_means), len(target_means))) + clipped = np.clip(target_means, -self.bias_clip, self.bias_clip) + circuit.h = np.arctanh(clipped) + return circuit + + def realize(self, manifold, target_point): + """ + Realize a pure-state manifold point with three p-bits. + + Returns (raw_means, reconstructed): raw_means is the direct p-bit + expectation estimate (may have norm < 1 -- see module docstring), + reconstructed is its projection onto CP^1. + """ + target = manifold.point(target_point) + circuit = self._build_independent_circuit(target) + + solver = CaSuDaSolver( + Nt=self.state_nt, dt=self.dt, i0=1.0, expected_mean=0.0, + seed=self._next_seed(), + ) + samples = solver.solve(circuit, n_shots=self.state_shots) + + assert samples.ndim == 3, ( + f"expected (Nt, n_shots, n_pbits) samples, got shape {samples.shape}" + ) + raw_means = np.mean(samples[self.state_burn_in:, :, :], axis=(0, 1)) + reconstructed = manifold.point(raw_means) + return raw_means, reconstructed + + def sample_measurement(self, manifold, point, axis, shots=1024): + """Sample a projective measurement using p-bit solver shots.""" + shots = int(shots) + if shots <= 0: + raise ValueError("shots must be positive.") + + expectation = manifold.born_expectation(point, axis) + circuit = self._build_independent_circuit([expectation]) + + solver = CaSuDaSolver( + Nt=self.measure_nt, + dt=self.dt, + i0=1.0, + expected_mean=0.0, + seed=self._next_seed(), + ) + + # CaSuDaSolver returns a 3-D trajectory for n_shots > 1: + # (Nt, n_shots, n_pbits). + # Use at least 2 solver shots to avoid the special n_shots == 1 return. + solver_shots = max(shots, 2) + samples = solver.solve(circuit, n_shots=solver_shots) + + assert samples.ndim == 3, ( + f"expected (Nt, n_shots, n_pbits), got {samples.shape}" + ) + + # One terminal p-bit value per logical measurement shot. + terminal = np.asarray(samples[-1, :shots, 0]) + + plus = int(np.sum(terminal > 0)) + + return { + "+": plus, + "-": shots - plus, + "shots": shots, + "empirical_p_plus": plus / shots, + "target_p_plus": 0.5 * (1.0 + expectation), + "expectation": expectation, + } + + def sample_one(self, manifold, point, axis): + result = self.sample_measurement(manifold, point, axis, shots=1) + return +1 if result["+"] == 1 else -1 + + +# --------------------------------------------------------------------------- +# Single-qubit manifold simulator +# --------------------------------------------------------------------------- + +class ManifoldPBitQubit: + H_AXIS = np.array([1.0, 0.0, 1.0]) / np.sqrt(2.0) + + def __init__(self, backend=None): + self.M = CP1BlochManifold() + self.backend = backend if backend is not None else PBitCP1Backend() + self.history = [] + self.reset() + + def reset(self): + north = self.M.Z_AXIS.copy() + self.ideal_state = north.copy() + self.raw_means, self.state = self.backend.realize(self.M, north) + self.history = [] + self._record("reset |0>") + return self + + def prepare(self, theta, phi): + target = self.M.from_angles(theta, phi) + self.ideal_state = target.copy() + self.raw_means, self.state = self.backend.realize(self.M, target) + self.history = [] + self._record(f"prepare(theta={theta:.4f}, phi={phi:.4f})") + return self + + def _record(self, gate): + raw_norm = float(np.linalg.norm(self.raw_means)) + d_theta, d_phi, _tangent_norm = tangent_decompose( + self.M, self.ideal_state, self.state + ) + self.history.append( + { + "gate": gate, + "ideal": self.ideal_state.copy(), + "state": self.state.copy(), + "raw": self.raw_means.copy(), + "raw_norm": raw_norm, + "purity": 0.5 * (1.0 + raw_norm ** 2), + "fidelity": self.M.fidelity(self.ideal_state, self.state), + "fs_error": self.M.distance(self.ideal_state, self.state), + "tangent_theta": d_theta, + "tangent_phi": d_phi, + } + ) + + def _rotation(self, axis, angle, name): + self.ideal_state = self.M.rotate(self.ideal_state, axis, angle) + target_from_pbits = self.M.rotate(self.state, axis, angle) + self.raw_means, self.state = self.backend.realize(self.M, target_from_pbits) + self._record(name) + return self + + def rx(self, angle): + return self._rotation(self.M.X_AXIS, angle, f"Rx({angle:.4f})") + + def ry(self, angle): + return self._rotation(self.M.Y_AXIS, angle, f"Ry({angle:.4f})") + + def rz(self, angle): + return self._rotation(self.M.Z_AXIS, angle, f"Rz({angle:.4f})") + + def x(self): + return self._rotation(self.M.X_AXIS, np.pi, "X") + + def y(self): + return self._rotation(self.M.Y_AXIS, np.pi, "Y") + + def z(self): + return self._rotation(self.M.Z_AXIS, np.pi, "Z") + + def h(self): + return self._rotation(self.H_AXIS, np.pi, "H") + + def s(self): + return self._rotation(self.M.Z_AXIS, np.pi / 2.0, "S") + + def t(self): + return self._rotation(self.M.Z_AXIS, np.pi / 4.0, "T") + + def bloch_vector(self): + return self.state.copy() + + def ideal_bloch_vector(self): + return self.ideal_state.copy() + + def fidelity_to_ideal(self): + return self.M.fidelity(self.ideal_state, self.state) + + def fs_error(self): + return self.M.distance(self.ideal_state, self.state) + + def sample_counts(self, axis="Z", shots=1024): + axis_vector = self._axis_vector(axis) + return self.backend.sample_measurement(self.M, self.state, axis_vector, shots=shots) + + def measure_once(self, axis="Z", collapse=True): + axis_vector = self._axis_vector(axis) + outcome = self.backend.sample_one(self.M, self.state, axis_vector) + if collapse: + collapsed = outcome * self.M.point(axis_vector) + self.ideal_state = collapsed.copy() + self.raw_means, self.state = self.backend.realize(self.M, collapsed) + self._record(f"measure {axis} -> {outcome:+d}") + return outcome + + def _axis_vector(self, axis): + if isinstance(axis, str): + name = axis.upper() + if name == "X": + return self.M.X_AXIS + if name == "Y": + return self.M.Y_AXIS + if name == "Z": + return self.M.Z_AXIS + raise ValueError("axis string must be 'X', 'Y', or 'Z'.") + return self.M.point(axis) + + def print_history(self): + print("gate/state history (single trial)") + print("==================================") + header = ( + f"{'operation':18s} {'fidelity':>10s} {'d_FS':>10s} " + f"{'purity':>8s} {'|raw|':>8s} {'d_theta':>9s} {'d_phi':>9s}" + ) + print(header) + print("-" * len(header)) + for row in self.history: + print( + f"{row['gate']:18s} " + f"{row['fidelity']:10.6f} " + f"{row['fs_error']:10.6f} " + f"{row['purity']:8.5f} " + f"{row['raw_norm']:8.5f} " + f"{row['tangent_theta']:9.5f} " + f"{row['tangent_phi']:9.5f}" + ) + print( + " (purity = effective pre-projection indicator " + "(1+|raw_means|^2)/2; d_theta/d_phi = tangent components of the " + "ideal -> reconstructed error, in the local polar/azimuthal frame)" + ) + +def run_single_qubit_trial(seed): + q = ManifoldPBitQubit(backend=PBitCP1Backend(seed=seed)) + for name, args in GATE_SEQUENCE: + getattr(q, name)(*args) + return q + + +def summarize_trials(qubits): + n_steps = len(qubits[0].history) + print(f"aggregate statistics over {len(qubits)} independent p-bit trials") + print("=" * 78) + header = ( + f"{'step':18s} {'fidelity':>12s} {'d_FS':>12s} " + f"{'purity':>12s} {'|d_theta|':>12s} {'|d_phi|':>12s}" + ) + print(header) + print("-" * len(header)) + for i in range(n_steps): + gate_name = qubits[0].history[i]["gate"] + fid = np.array([q.history[i]["fidelity"] for q in qubits]) + dfs = np.array([q.history[i]["fs_error"] for q in qubits]) + pur = np.array([q.history[i]["purity"] for q in qubits]) + dth = np.array([abs(q.history[i]["tangent_theta"]) for q in qubits]) + dph = np.array([abs(q.history[i]["tangent_phi"]) for q in qubits]) + + print( + f"{gate_name:18s} {fid.mean():12.6f} {dfs.mean():12.6f} " + f"{pur.mean():12.6f} {dth.mean():12.6f} {dph.mean():12.6f}" + ) + print( + f"{' (std)':18s} {fid.std():12.6f} {dfs.std():12.6f} " + f"{pur.std():12.6f} {dth.std():12.6f} {dph.std():12.6f}" + ) + print() + + +def main(n_trials=N_TRIALS): + print("Single-qubit manifold-first p-bit simulator") + print("=====================================================") + print(f"Running {n_trials} independent trials of H, T, Ry(pi/3), Rz(-pi/5)...") + print() + + trials = [run_single_qubit_trial(seed=BASE_SEED + i) for i in range(n_trials)] + + trials[0].print_history() + print() + summarize_trials(trials) + + q = trials[0] + print("final state (trial 0)") + print("======================") + print("ideal Bloch vector :", q.ideal_bloch_vector()) + print("p-bit Bloch vector :", q.bloch_vector()) + print(f"final fidelity : {q.fidelity_to_ideal():.6f}") + print(f"final d_FS error : {q.fs_error():.6f}") + print() + + shots = 2048 + print(f"measurement counts ({shots} shots, trial 0 final state)") + print("=========================================================") + print("axis ideal P(+) manifold P(+) p-bit P(+) counts") + print("-" * 66) + for axis_name, axis_vector in [ + ("X", q.M.X_AXIS), ("Y", q.M.Y_AXIS), ("Z", q.M.Z_AXIS), + ("XY", np.array([1.0, 1.0, 0.0])), + ]: + ideal_p = q.M.born_probability(q.ideal_bloch_vector(), axis_vector) + manifold_p = q.M.born_probability(q.bloch_vector(), axis_vector) + result = q.sample_counts(axis_vector, shots=shots) + print( + f"{axis_name:4s} {ideal_p:10.4f} {manifold_p:13.4f} " + f"{result['empirical_p_plus']:10.4f} " + f"+:{result['+']:4d} -:{result['-']:4d}" + ) + + print() + print("projective-collapse demonstration") + print("==================================") + before = q.bloch_vector() + outcome = q.measure_once("Z", collapse=True) + after = q.bloch_vector() + print("state before Z measurement :", before) + print(f"observed outcome : {outcome:+d}") + print("post-measurement p-bit state:", after) + print(f"fidelity to collapsed pole : {q.fidelity_to_ideal():.6f}") + repeat = q.sample_counts("Z", shots=shots) + print( + f"repeat Z measurement : +:{repeat['+']} -:{repeat['-']} " + f"(P(+)= {repeat['empirical_p_plus']:.4f})" + ) + print() + +if __name__ == "__main__": + main() From 838f3c2108cd41aacbe668b799af01ba32106de4 Mon Sep 17 00:00:00 2001 From: Anton Date: Tue, 15 Sep 2026 11:03:25 +0200 Subject: [PATCH 2/2] Moved demo 1 to qunatum_folder Added demo 2: calibrate and execute a quantum circuit with a p-bit/Riemannian runtime backend. --- .../demo_001_quantum_manifold_pbit.py} | 34 +- ..._quantum_pbit_runtime_noise_calibration.py | 1180 +++++++++++++++++ 2 files changed, 1200 insertions(+), 14 deletions(-) rename examples/{quantum_manifold_pbit.py => quantum_sim/demo_001_quantum_manifold_pbit.py} (94%) create mode 100644 examples/quantum_sim/demo_002_quantum_pbit_runtime_noise_calibration.py diff --git a/examples/quantum_manifold_pbit.py b/examples/quantum_sim/demo_001_quantum_manifold_pbit.py similarity index 94% rename from examples/quantum_manifold_pbit.py rename to examples/quantum_sim/demo_001_quantum_manifold_pbit.py index dc1560c..ac9e8ed 100644 --- a/examples/quantum_manifold_pbit.py +++ b/examples/quantum_sim/demo_001_quantum_manifold_pbit.py @@ -1,8 +1,7 @@ """ - This is a single-qubit research demo. The objective of the demo is to make a link between p-bit probabilistic computing, Riemannian geometry, and the -quantum Bloch sphere. It explores this combination for a potential future +quantum Bloch sphere. It explores this combination for a potential future quantum simulator. A pure single-qubit state is treated first as a point on the Riemannian @@ -11,40 +10,47 @@ coordinates, gates act as CP^1 isometries, and one additional p-bit samples projective measurements. -So the the qubit lives on CP^1, p-bits realize that state statistically, +So the qubit lives on CP^1, p-bits realize that state statistically, gates move it geometrically, and measurements are sampled probabilistically. -GATE_SEQUENCE represents a small quantum circuit test. +GATE_SEQUENCE represents a small quantum circuit test. ManifoldPBitQubit # single-qubit simulator object ├── CP1BlochManifold # geometry └── PBitCP1Backend # realization using p-bits The demo reports: - * the raw p-bit Bloch-vector norm before projection back to CP^1; - * an effective pre-projection purity indicator (1 + |raw_means|^2) / 2; - * Fubini-Study reconstruction error and target-state fidelity; - * tangent-space error decomposition into local polar/azimuthal components; - * mean and standard deviation over multiple independent p-bit trials; - * Born-rule measurement counts and projective-collapse behavior. + +* the raw p-bit Bloch-vector norm before projection back to CP^1; +* an effective pre-projection purity indicator (1 + |raw_means|^2) / 2; +* Fubini-Study reconstruction error and target-state fidelity; +* tangent-space error decomposition into local polar/azimuthal components; +* mean and standard deviation over multiple independent p-bit trials; +* Born-rule measurement counts and projective-collapse behavior. The raw norm/purity values are diagnostics of the stochastic state estimate; finite sampling and finite p-bit bias can move the raw estimate away from the unit sphere, so they should not by themselves be interpreted as physical mixed-state tomography. +The final test shows that after applying different gates to the qubit, the +p-bit representation reproduces the stochastic measurement distribution +expected for the resulting single-qubit state. Each gate changes the quantum +state, and therefore the expected measurement probabilities, while the same +p-bit encoding and sampling mechanism is used throughout. + Results: - The quantum manifold p-bit representation is numerically stable and reproduces - single-qubit circuit behavior well. - +The quantum manifold p-bit representation is numerically stable and reproduces +single-qubit circuit behavior well. + The question is how far we can go from here. Could coupled p-bits provide a stochastic representation of joint probability distributions and correlations associated with points on higher-dimensional quantum-state manifolds? - """ + import numpy as np from p_kit.psl import PCircuit diff --git a/examples/quantum_sim/demo_002_quantum_pbit_runtime_noise_calibration.py b/examples/quantum_sim/demo_002_quantum_pbit_runtime_noise_calibration.py new file mode 100644 index 0000000..edb7531 --- /dev/null +++ b/examples/quantum_sim/demo_002_quantum_pbit_runtime_noise_calibration.py @@ -0,0 +1,1180 @@ +""" +Demo 2: calibrate and execute a quantum circuit with a p-bit/Riemannian runtime backend. + +## Purpose + +Demo 1 establishes the basic single-qubit representation: + +Demo 2 follows a more familiar quantum-computing workflow and adds an explicit +noise model and calibration loop: + +## Workflow + +1. Define an intended quantum circuit. +2. Define backend imperfections using the NoiseModel class. +3. Execute the circuit with PBitRuntimeBackend. +4. Realize the output state statistically with p-bits. +5. Estimate the resulting quantum state on CP^1 using a Riemannian mean. +6. Compare the produced state with the ideal target using Fubini-Study distance. +7. Calibrate the programmable gate parameters. +8. Execute the corrected circuit again. +9. Compare state fidelity and Born measurement probabilities before/after calibration. + +## Main components + +QuantumCircuit +Describes the intended quantum computation and programmable gate parameters. + +NoiseModel +Explicitly describes backend imperfections. In this demo it contains +reproducible coherent gate-angle errors. The class is separate from the +quantum circuit so that the intended computation remains ideal while the +runtime backend models imperfect execution. + +PBitRuntimeBackend +Executes the circuit using the geometric single-qubit model and realizes +the resulting state statistically with p-bits. + +Quantum-state manifold +A pure single-qubit state is represented as a point on CP^1, equivalently +the Bloch sphere S^2. The manifold constrains the stochastic representation +to valid pure quantum states. + +Riemannian state estimator +Combines repeated noisy p-bit state realizations using an intrinsic +Frechet/Karcher mean on the quantum-state manifold. + +Fubini-Study / geometric loss +Measures the distance between the estimated runtime state and the desired +ideal target state. + +Calibration +Adjusts the programmable gate commands so that execution on the imperfect +stochastic runtime backend produces a state closer to the desired target. + +Corrected QuantumCircuit +The calibrated circuit is specialized to the imperfect runtime backend. +Its command parameters therefore do not necessarily equal the ideal +parameters. + +## Interpretation + +The p-bits provide the stochastic realization of the quantum state. + +CP^1 geometry provides the state constraint, the Riemannian estimator, and the +Fubini-Study calibration objective. + +The NoiseModel class makes runtime imperfections explicit rather than hiding +them inside the circuit or the backend implementation. + +The main question explored by this demo is: + +Can a noisy stochastic p-bit quantum runtime be calibrated using +quantum-state Riemannian geometry so that it still produces the +expected quantum-state result? + +## IMPORTANT + +The gate imperfections used in this file are deliberately simulated so that +the calibration step has a reproducible error to correct. They are not meant +to model a particular physical quantum processor or p-bit device. + +The same NoiseModel/runtime interface could later be replaced or populated +with experimentally measured distortions from a physical Probana/p-kit +backend. + +This is a stochastic research demo, not a physical quantum computer. +""" + +from dataclasses import dataclass +import copy +import numpy as np + +from p_kit.psl import PCircuit +from p_kit.solver.csd_solver import CaSuDaSolver + + +# ============================================================================ +# Demo configuration +# ============================================================================ + +TARGET_THETA = 0.95 +TARGET_PHI = 0.70 + +# p-bit state realization +STATE_NT = 2_000 +STATE_BURN_IN = 400 +STATE_SHOTS = 4 +N_REALIZATIONS = 7 +DT = 0.02 +BIAS_CLIP = 0.995 + +# Measurement +MEASUREMENT_SHOTS = 2048 +MEASUREMENT_NT = 250 + +# Reproducibility +BASE_REALIZATION_SEED = 1000 +MEASUREMENT_SEED = 9000 + +# Calibration search +INITIAL_CALIBRATION_STEP = 0.08 +MIN_CALIBRATION_STEP = 0.00125 +MAX_CALIBRATION_ROUNDS = 18 + +# -------------------------------------------------------------------------- +# Explicit demo noise parameters. +# +# These are deliberately injected coherent gate-angle errors. They make the +# noisy-backend calibration reproducible and visible: +# +# actual Rz = gain * commanded Rz + offset +# actual Ry = gain * commanded Ry + offset +# +# They are not claimed to model a particular physical device. +# -------------------------------------------------------------------------- + +RZ_NOISE_OFFSET = 0.045 +RY_NOISE_GAIN = 0.960 +RY_NOISE_OFFSET = -0.025 + + +# ============================================================================ +# CP^1 / Bloch-sphere geometry +# ============================================================================ + +class CP1BlochManifold: + """Pure single-qubit state manifold CP^1 represented by Bloch S^2.""" + + X_AXIS = np.array([1.0, 0.0, 0.0]) + Y_AXIS = np.array([0.0, 1.0, 0.0]) + Z_AXIS = np.array([0.0, 0.0, 1.0]) + + @staticmethod + def point(vector): + vector = np.asarray(vector, dtype=float) + norm = np.linalg.norm(vector) + + if norm == 0.0: + raise ValueError("The zero vector cannot define a pure qubit state.") + + return vector / norm + + @classmethod + def tangent(cls, point, vector): + point = cls.point(point) + vector = np.asarray(vector, dtype=float) + return vector - np.dot(vector, point) * point + + @classmethod + def exp_map(cls, point, tangent_vector): + """ + Sphere exponential map. + + The sphere geodesic angle is twice the Fubini-Study distance, but the + same intrinsic mean point is obtained because that constant metric + scaling does not change the minimizer. + """ + point = cls.point(point) + tangent_vector = cls.tangent(point, tangent_vector) + angle = np.linalg.norm(tangent_vector) + + if angle < 1e-15: + return point.copy() + + return cls.point( + np.cos(angle) * point + + np.sin(angle) * tangent_vector / angle + ) + + @classmethod + def log_map(cls, point_a, point_b): + """Sphere logarithmic map used by the intrinsic mean.""" + point_a = cls.point(point_a) + point_b = cls.point(point_b) + + dot = float(np.clip(np.dot(point_a, point_b), -1.0, 1.0)) + angle = float(np.arccos(dot)) + + if angle < 1e-15: + return np.zeros(3) + + if np.pi - angle < 1e-12: + raise ValueError("Log map is not unique for antipodal points.") + + direction = point_b - dot * point_a + direction /= np.linalg.norm(direction) + + return angle * direction + + @classmethod + def fs_distance(cls, point_a, point_b): + """Fubini-Study distance for pure single-qubit states.""" + point_a = cls.point(point_a) + point_b = cls.point(point_b) + + dot = float(np.clip(np.dot(point_a, point_b), -1.0, 1.0)) + + return 0.5 * float(np.arccos(dot)) + + @classmethod + def fidelity(cls, point_a, point_b): + """Pure-state fidelity.""" + point_a = cls.point(point_a) + point_b = cls.point(point_b) + + return float( + np.clip( + 0.5 * (1.0 + np.dot(point_a, point_b)), + 0.0, + 1.0, + ) + ) + + @classmethod + def rotate(cls, point, axis, angle): + """ + Single-qubit unitary action represented as a Bloch-sphere rotation. + """ + point = cls.point(point) + axis = cls.point(axis) + + c = np.cos(angle) + s = np.sin(angle) + + return cls.point( + c * point + + s * np.cross(axis, point) + + (1.0 - c) * np.dot(axis, point) * axis + ) + + @classmethod + def born_expectation(cls, point, axis): + point = cls.point(point) + axis = cls.point(axis) + + return float(np.clip(np.dot(point, axis), -1.0, 1.0)) + + @classmethod + def born_probability(cls, point, axis): + return 0.5 * (1.0 + cls.born_expectation(point, axis)) + + +M = CP1BlochManifold() +H_AXIS = M.point([1.0, 0.0, 1.0]) + + +# ============================================================================ +# Geometry-aware averaging +# ============================================================================ + +def projected_euclidean_mean(points): + """Mean in R^3 followed by projection back to CP^1.""" + points = np.asarray(points, dtype=float) + return M.point(np.mean(points, axis=0)) + + +def riemannian_mean(points, max_iter=64, tol=1e-10): + """ + Intrinsic Frechet/Karcher mean on CP^1. + + This is the geometry-aware state estimate used by the calibration loss. + """ + points = np.asarray(points, dtype=float) + mean = projected_euclidean_mean(points) + + for _ in range(max_iter): + logs = np.asarray([M.log_map(mean, point) for point in points]) + update = np.mean(logs, axis=0) + + if np.linalg.norm(update) < tol: + break + + mean = M.exp_map(mean, update) + + return M.point(mean) + + +# ============================================================================ +# Minimal single-qubit circuit interface +# ============================================================================ + +@dataclass +class Gate: + name: str + angle: float | None = None + label: str | None = None + + +class QuantumCircuit: + """ + Minimal one-qubit circuit used only by this research demo. + + The interface intentionally resembles a normal quantum-circuit workflow. + """ + + def __init__(self): + self.gates = [] + + def copy(self): + return copy.deepcopy(self) + + def h(self): + self.gates.append(Gate("h")) + return self + + def x(self): + self.gates.append(Gate("x")) + return self + + def y(self): + self.gates.append(Gate("y")) + return self + + def z(self): + self.gates.append(Gate("z")) + return self + + def t(self): + self.gates.append(Gate("t")) + return self + + def rx(self, angle, label=None): + self.gates.append(Gate("rx", float(angle), label)) + return self + + def ry(self, angle, label=None): + self.gates.append(Gate("ry", float(angle), label)) + return self + + def rz(self, angle, label=None): + self.gates.append(Gate("rz", float(angle), label)) + return self + + def set_parameter(self, label, value): + found = False + + for gate in self.gates: + if gate.label == label: + gate.angle = float(value) + found = True + + if not found: + raise KeyError(f"No circuit parameter labelled {label!r}.") + + return self + + def get_parameter(self, label): + for gate in self.gates: + if gate.label == label: + return float(gate.angle) + + raise KeyError(f"No circuit parameter labelled {label!r}.") + + def __str__(self): + parts = ["|0>"] + + for gate in self.gates: + if gate.angle is None: + parts.append(f"-- {gate.name.upper()}") + elif gate.label is None: + parts.append(f"-- {gate.name.upper()}({gate.angle:.4f})") + else: + parts.append( + f"-- {gate.name.upper()}({gate.label}={gate.angle:.4f})" + ) + + return " ".join(parts) + + +# ============================================================================ +# Exact geometric circuit execution +# ============================================================================ + +def apply_gate_ideal(state, gate): + """Apply one ideal gate to a Bloch point.""" + if gate.name == "h": + return M.rotate(state, H_AXIS, np.pi) + + if gate.name == "x": + return M.rotate(state, M.X_AXIS, np.pi) + + if gate.name == "y": + return M.rotate(state, M.Y_AXIS, np.pi) + + if gate.name == "z": + return M.rotate(state, M.Z_AXIS, np.pi) + + if gate.name == "t": + return M.rotate(state, M.Z_AXIS, np.pi / 4.0) + + if gate.name == "rx": + return M.rotate(state, M.X_AXIS, gate.angle) + + if gate.name == "ry": + return M.rotate(state, M.Y_AXIS, gate.angle) + + if gate.name == "rz": + return M.rotate(state, M.Z_AXIS, gate.angle) + + raise ValueError(f"Unsupported gate: {gate.name!r}") + + +def execute_ideal(circuit): + """Execute a circuit exactly on CP^1.""" + state = M.Z_AXIS.copy() + + for gate in circuit.gates: + state = apply_gate_ideal(state, gate) + + return M.point(state) + + +# ============================================================================ +# Explicit noise model +# ============================================================================ + +@dataclass(frozen=True) +class GateAngleError: + """ + Coherent single-gate angle error. + + effective_angle = gain * commanded_angle + offset + """ + gain: float = 1.0 + offset: float = 0.0 + + def apply(self, angle): + return self.gain * float(angle) + self.offset + + +class NoiseModel: + """ + Minimal Qiskit-Aer like noise container for this research demo. + + The noise model belongs to the runtime backend, not to the ideal QuantumCircuit. + That keeps the intended circuit separate from imperfections of the device + or simulator used to execute it. + """ + + def __init__(self): + self.gate_angle_errors = {} + + def add_gate_angle_error(self, gate_name, gain=1.0, offset=0.0): + gate_name = str(gate_name).lower() + + if gate_name not in {"rx", "ry", "rz"}: + raise ValueError( + "Gate-angle noise is supported only for rx, ry and rz." + ) + + self.gate_angle_errors[gate_name] = GateAngleError( + gain=float(gain), + offset=float(offset), + ) + return self + + def apply(self, gate): + distorted = copy.deepcopy(gate) + error = self.gate_angle_errors.get(gate.name) + + if error is not None: + distorted.angle = error.apply(gate.angle) + + return distorted + + def describe(self): + if not self.gate_angle_errors: + return ["ideal noise model"] + + lines = [] + + for gate_name in sorted(self.gate_angle_errors): + error = self.gate_angle_errors[gate_name] + lines.append( + f"{gate_name.upper()}: actual angle = " + f"{error.gain:.6f} * commanded angle " + f"{error.offset:+.6f} rad" + ) + + return lines + + +# ============================================================================ +# p-bit / Riemannian runtime backend +# ============================================================================ + +class PBitRuntimeBackend: + """ + Stochastic p-bit single-qubit runtime backend. + + Gate evolution is represented geometrically on CP^1. The resulting state + is realized statistically by three p-bits whose means represent the Bloch + coordinates. + + An explicit NoiseModel is attached to the backend so noisy execution and + calibration are visible and reproducible. + """ + + def __init__( + self, + noise_model=None, + state_nt=STATE_NT, + state_burn_in=STATE_BURN_IN, + state_shots=STATE_SHOTS, + dt=DT, + bias_clip=BIAS_CLIP, + ): + self.noise_model = noise_model if noise_model is not None else NoiseModel() + + self.state_nt = int(state_nt) + self.state_burn_in = int(state_burn_in) + self.state_shots = int(state_shots) + self.dt = float(dt) + self.bias_clip = float(bias_clip) + + if self.state_shots <= 1: + raise ValueError("state_shots must be > 1.") + + def _distort_gate(self, gate): + """Apply the explicit runtime-backend noise model to one gate.""" + return self.noise_model.apply(gate) + + def effective_geometric_state(self, circuit): + """ + Deterministic geometric state after applying the backend distortion. + + This is not used as the calibration observation. Calibration uses + stochastic p-bit realizations below. + """ + state = M.Z_AXIS.copy() + + for gate in circuit.gates: + effective_gate = self._distort_gate(gate) + state = apply_gate_ideal(state, effective_gate) + + return M.point(state) + + def _realize_state_once(self, target_state, seed): + """ + Realize one Bloch point statistically using three independent p-bits. + """ + target_state = M.point(target_state) + + circuit = PCircuit(3) + circuit.J = np.zeros((3, 3)) + + clipped = np.clip( + target_state, + -self.bias_clip, + self.bias_clip, + ) + circuit.h = np.arctanh(clipped) + + solver = CaSuDaSolver( + Nt=self.state_nt, + dt=self.dt, + i0=1.0, + expected_mean=0.0, + seed=int(seed), + ) + + # CaSuDaSolver multi-shot path: + # samples.shape == (Nt, n_shots, n_pbits) + samples = solver.solve( + circuit, + n_shots=self.state_shots, + ) + + assert samples.ndim == 3, ( + f"expected (Nt, n_shots, n_pbits), got {samples.shape}" + ) + + raw_mean = np.mean( + samples[self.state_burn_in:, :, :], + axis=(0, 1), + ) + + return raw_mean, M.point(raw_mean) + + def run( + self, + circuit, + n_realizations=N_REALIZATIONS, + seed_base=BASE_REALIZATION_SEED, + ): + """ + Execute a circuit and estimate its output state from p-bit realizations. + """ + effective_state = self.effective_geometric_state(circuit) + + raw_vectors = [] + points = [] + + for index in range(int(n_realizations)): + raw, point = self._realize_state_once( + effective_state, + seed=int(seed_base) + index, + ) + + raw_vectors.append(raw) + points.append(point) + + raw_vectors = np.asarray(raw_vectors) + points = np.asarray(points) + + intrinsic_mean = riemannian_mean(points) + + return { + "effective_state": effective_state, + "raw_vectors": raw_vectors, + "points": points, + "state": intrinsic_mean, + "raw_radius_mean": float( + np.mean(np.linalg.norm(raw_vectors, axis=1)) + ), + "dispersion_fs": float( + np.mean( + [ + M.fs_distance(intrinsic_mean, point) + for point in points + ] + ) + ), + } + + def sample_measurement( + self, + state, + axis, + shots=MEASUREMENT_SHOTS, + seed=MEASUREMENT_SEED, + ): + """ + Sample a projective measurement with one p-bit. + + The p-bit expectation is set to the Born expectation r dot n. + """ + shots = int(shots) + + if shots <= 0: + raise ValueError("shots must be positive.") + + expectation = M.born_expectation(state, axis) + + circuit = PCircuit(1) + circuit.J = np.zeros((1, 1)) + + clipped = float( + np.clip( + expectation, + -self.bias_clip, + self.bias_clip, + ) + ) + circuit.h = np.array([np.arctanh(clipped)]) + + solver = CaSuDaSolver( + Nt=MEASUREMENT_NT, + dt=self.dt, + i0=1.0, + expected_mean=0.0, + seed=int(seed), + ) + + solver_shots = max(shots, 2) + samples = solver.solve( + circuit, + n_shots=solver_shots, + ) + + assert samples.ndim == 3, ( + f"expected (Nt, n_shots, n_pbits), got {samples.shape}" + ) + + terminal = np.asarray(samples[-1, :shots, 0]) + plus = int(np.sum(terminal > 0)) + minus = shots - plus + + return { + "+": plus, + "-": minus, + "shots": shots, + "empirical_p_plus": plus / shots, + "target_p_plus": M.born_probability(state, axis), + } + + +# ============================================================================ +# Calibration +# ============================================================================ + +class CalibrationEvaluator: + """ + Evaluate candidate commanded parameters using the stochastic backend. + + Candidate evaluations are cached. Common seeds are used for all + candidates so the optimization is not comparing unrelated random draws. + """ + + def __init__( + self, + backend, + template_circuit, + target_state, + theta_label="theta", + phi_label="phi", + ): + self.backend = backend + self.template_circuit = template_circuit.copy() + self.target_state = M.point(target_state) + + self.theta_label = theta_label + self.phi_label = phi_label + + self.cache = {} + self.execution_count = 0 + + @staticmethod + def _key(theta, phi): + return ( + round(float(theta), 12), + round(float(phi), 12), + ) + + def evaluate(self, theta, phi): + key = self._key(theta, phi) + + if key in self.cache: + return self.cache[key] + + candidate = self.template_circuit.copy() + candidate.set_parameter(self.theta_label, theta) + candidate.set_parameter(self.phi_label, phi) + + execution = self.backend.run( + candidate, + n_realizations=N_REALIZATIONS, + seed_base=BASE_REALIZATION_SEED, + ) + + state = execution["state"] + + result = { + "theta": float(theta), + "phi": float(phi), + "circuit": candidate, + "execution": execution, + "distance": M.fs_distance(state, self.target_state), + "fidelity": M.fidelity(state, self.target_state), + } + + self.cache[key] = result + self.execution_count += 1 + + return result + + +def calibrate_two_parameters( + evaluator, + theta_initial, + phi_initial, +): + """ + Simple derivative-free coordinate search. + + The objective is the Fubini-Study distance between the target ideal state + and the Riemannian mean of repeated p-bit executions. + """ + theta = float(theta_initial) + phi = float(phi_initial) + + step_theta = INITIAL_CALIBRATION_STEP + step_phi = INITIAL_CALIBRATION_STEP + + history = [] + + current = evaluator.evaluate(theta, phi) + + for round_index in range(MAX_CALIBRATION_ROUNDS): + candidates = [ + (theta, phi), + (theta + step_theta, phi), + (theta - step_theta, phi), + (theta, phi + step_phi), + (theta, phi - step_phi), + ] + + results = [ + evaluator.evaluate(candidate_theta, candidate_phi) + for candidate_theta, candidate_phi in candidates + ] + + best = min(results, key=lambda result: result["distance"]) + + if best["distance"] + 1e-15 < current["distance"]: + theta = best["theta"] + phi = best["phi"] + current = best + else: + step_theta *= 0.5 + step_phi *= 0.5 + + history.append( + { + "round": round_index, + "theta": theta, + "phi": phi, + "distance": current["distance"], + "fidelity": current["fidelity"], + "step_theta": step_theta, + "step_phi": step_phi, + } + ) + + if max(step_theta, step_phi) < MIN_CALIBRATION_STEP: + break + + return { + "theta": theta, + "phi": phi, + "result": current, + "history": history, + } + + +# ============================================================================ +# Reporting helpers +# ============================================================================ + +def print_state_comparison( + label, + state, + target_state, +): + print(label) + print("-" * len(label)) + print(f"state = {np.array2string(state, precision=7)}") + print( + f"d_FS = {M.fs_distance(state, target_state):.8f}" + ) + print( + f"fidelity = {M.fidelity(state, target_state):.9f}" + ) + print() + + +def print_calibration_history(history): + print("CALIBRATION") + print("===========") + print( + f"{'round':>5s} " + f"{'theta_cmd':>11s} " + f"{'phi_cmd':>11s} " + f"{'d_FS':>11s} " + f"{'fidelity':>12s} " + f"{'step':>9s}" + ) + print("-" * 68) + + for row in history: + print( + f"{row['round']:5d} " + f"{row['theta']:11.6f} " + f"{row['phi']:11.6f} " + f"{row['distance']:11.7f} " + f"{row['fidelity']:12.9f} " + f"{max(row['step_theta'], row['step_phi']):9.6f}" + ) + + print() + + +def measurement_table( + backend, + target_state, + before_state, + after_state, +): + print("MEASUREMENT CHECK") + print("=================") + print( + f"{'axis':>4s} " + f"{'ideal P(+)':>12s} " + f"{'before P(+)':>13s} " + f"{'after P(+)':>12s} " + f"{'after sampled':>14s}" + ) + print("-" * 63) + + axes = ( + ("X", M.X_AXIS), + ("Y", M.Y_AXIS), + ("Z", M.Z_AXIS), + ) + + for axis_index, (name, axis) in enumerate(axes): + ideal_p = M.born_probability(target_state, axis) + before_p = M.born_probability(before_state, axis) + after_p = M.born_probability(after_state, axis) + + sampled = backend.sample_measurement( + after_state, + axis, + shots=MEASUREMENT_SHOTS, + seed=MEASUREMENT_SEED + axis_index, + ) + + print( + f"{name:>4s} " + f"{ideal_p:12.6f} " + f"{before_p:13.6f} " + f"{after_p:12.6f} " + f"{sampled['empirical_p_plus']:14.6f}" + ) + + print() + + +# ============================================================================ +# Main demo +# ============================================================================ + +def main(): + print("Demo 2: p-bit quantum-circuit calibration + execution") + print("======================================================") + print() + + # ---------------------------------------------------------------------- + # 1. Define the intended quantum circuit in the normal forward direction. + # ---------------------------------------------------------------------- + + intended_circuit = QuantumCircuit() + intended_circuit.h() + intended_circuit.rz(TARGET_THETA, label="theta") + intended_circuit.ry(TARGET_PHI, label="phi") + intended_circuit.t() + + print("INTENDED CIRCUIT") + print("================") + print(intended_circuit) + print() + + # ---------------------------------------------------------------------- + # 2. Compute the ideal target state. + # ---------------------------------------------------------------------- + + target_state = execute_ideal(intended_circuit) + + print("Ideal target state:") + print(np.array2string(target_state, precision=7)) + print() + + # ---------------------------------------------------------------------- + # 3. Build an explicit noise model and attach it to the runtime backend. + # + # This mirrors the usual simulator pattern: + # + # circuit -> backend(noise_model=...) -> result + # + # The ideal circuit itself remains unchanged. + # ---------------------------------------------------------------------- + + noise_model = NoiseModel() + noise_model.add_gate_angle_error( + "rz", + offset=RZ_NOISE_OFFSET, + ) + noise_model.add_gate_angle_error( + "ry", + gain=RY_NOISE_GAIN, + offset=RY_NOISE_OFFSET, + ) + + backend = PBitRuntimeBackend( + noise_model=noise_model, + ) + + print("NOISE MODEL") + print("===========") + for line in noise_model.describe(): + print(line) + print( + "(These coherent gate errors are deliberately injected for the demo.)" + ) + print() + + # ---------------------------------------------------------------------- + # 4. Execute the intended circuit before calibration. + # ---------------------------------------------------------------------- + + before = backend.run( + intended_circuit, + n_realizations=N_REALIZATIONS, + seed_base=BASE_REALIZATION_SEED, + ) + + before_state = before["state"] + + print_state_comparison( + "BEFORE CALIBRATION", + before_state, + target_state, + ) + + print( + f"mean p-bit realization dispersion d_FS = " + f"{before['dispersion_fs']:.8f}" + ) + print( + f"mean raw Bloch radius = " + f"{before['raw_radius_mean']:.8f}" + ) + print() + + # ---------------------------------------------------------------------- + # 5. Calibrate the programmable Rz/Ry parameters. + # ---------------------------------------------------------------------- + + evaluator = CalibrationEvaluator( + backend=backend, + template_circuit=intended_circuit, + target_state=target_state, + theta_label="theta", + phi_label="phi", + ) + + calibration = calibrate_two_parameters( + evaluator, + theta_initial=TARGET_THETA, + phi_initial=TARGET_PHI, + ) + + print_calibration_history(calibration["history"]) + + calibrated_theta = calibration["theta"] + calibrated_phi = calibration["phi"] + + calibrated_circuit = intended_circuit.copy() + calibrated_circuit.set_parameter("theta", calibrated_theta) + calibrated_circuit.set_parameter("phi", calibrated_phi) + + # ---------------------------------------------------------------------- + # 6. Re-execute the calibrated circuit. + # ---------------------------------------------------------------------- + + after = backend.run( + calibrated_circuit, + n_realizations=N_REALIZATIONS, + seed_base=BASE_REALIZATION_SEED, + ) + + after_state = after["state"] + + print("CALIBRATED CIRCUIT") + print("==================") + print(calibrated_circuit) + print() + + print_state_comparison( + "AFTER CALIBRATION", + after_state, + target_state, + ) + + # ---------------------------------------------------------------------- + # 7. Summary. + # ---------------------------------------------------------------------- + + before_distance = M.fs_distance(before_state, target_state) + after_distance = M.fs_distance(after_state, target_state) + + before_fidelity = M.fidelity(before_state, target_state) + after_fidelity = M.fidelity(after_state, target_state) + + print("FINAL COMPARISON") + print("================") + print( + f"intended theta = {TARGET_THETA:.6f}" + ) + print( + f"calibrated commanded theta = {calibrated_theta:.6f}" + ) + print( + f"intended phi = {TARGET_PHI:.6f}" + ) + print( + f"calibrated commanded phi = {calibrated_phi:.6f}" + ) + print() + + print( + f"before d_FS = {before_distance:.8f}" + ) + print( + f"after d_FS = {after_distance:.8f}" + ) + print( + f"before fidelity = {before_fidelity:.9f}" + ) + print( + f"after fidelity = {after_fidelity:.9f}" + ) + + if after_distance > 0.0: + print( + f"FS-error reduction factor = " + f"{before_distance / after_distance:.3f}x" + ) + + print( + f"candidate runtime-backend executions = " + f"{evaluator.execution_count}" + ) + print() + + # The ideal programmed circuit with the calibrated parameters is not + # expected to equal the target. The calibration is compensating the + # imperfect backend, not recovering the original ideal parameters. + ideal_calibrated_state = execute_ideal(calibrated_circuit) + + print( + "Ideal simulator d_FS at the calibrated commanded parameters = " + f"{M.fs_distance(ideal_calibrated_state, target_state):.8f}" + ) + print( + "This can remain non-zero because the calibrated commands are chosen " + "for the imperfect p-bit backend." + ) + print() + + # ---------------------------------------------------------------------- + # 8. Measurement verification. + # ---------------------------------------------------------------------- + + measurement_table( + backend, + target_state, + before_state, + after_state, + ) + + print("INTERPRETATION") + print("==============") + print( + "The circuit is first executed normally on a stochastic p-bit backend. " + "The resulting CP^1 state is estimated with a Riemannian mean. " + "Fubini-Study distance to the ideal target is then used as the " + "calibration objective. The corrected circuit is finally executed " + "again on the same backend." + ) + print() + print( + "The key result is not that the calibrated command angles reproduce " + "the ideal command angles. The goal is that the calibrated p-bit " + "backend output moves closer to the desired quantum state." + ) + + +if __name__ == "__main__": + main()