Repository navigation
Expand file tree
/
Copy pathbase_code_doing_stat_praparation_no_amp.py
More file actions
305 lines (247 loc) · 9.88 KB
/
Copy pathbase_code_doing_stat_praparation_no_amp.py
File metadata and controls
305 lines (247 loc) · 9.88 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
# %%
# ============================================================
# BASELINE: ALL WALKS (uniform superposition) +
# VALIDITY ORACLE (parity-based) +
# FULL COST PHASE ORACLE
# Then: statevector analysis + classical comparison table
#
# No amplitude amplification, no QSVT.
# ============================================================
import numpy as np
import pandas as pd
import pennylane as qml
import src.classical_funcs as cf
# -------------------------
# Problem + cost matrix
# -------------------------
np.random.seed(42)
n = 5
start_node = n - 1
cost_matrix_raw = cf.generate_cost_matrix(n)
all_walks = cf.generate_all_walks(n, start_node=start_node)
all_costs_raw = cf.find_all_cost(cost_matrix_raw, all_walks)
C_max = float(np.max(all_costs_raw))
cost_matrix = cost_matrix_raw / C_max # entries are phases in radians in [0, ~1]
# -------------------------
# Encoding (time-register)
# -------------------------
T_steps = n - 1
n_qubits_step = int(np.ceil(np.log2(n - 1)))
def twire(t, q):
return f"t{t}_q{q}"
def pos_wires_at_t(t):
return [twire(t, q) for q in range(n_qubits_step)]
time_wires = [twire(t, q) for t in range(T_steps) for q in range(n_qubits_step)]
state_wires = [f"s{i}" for i in range(n - 1)] # parity flags
good_wire = "good"
phase_anc = "phase"
# Optional work wires for big MCX; not required here (controls are small), but kept consistent
m = len(time_wires)
work_wires = [f"w{k}" for k in range(max(0, m - 2))]
wires = time_wires + state_wires + [good_wire, phase_anc] + work_wires
dev = qml.device("lightning.qubit", wires=wires)
# -------------------------
# Init: uniform over ALL time-register bitstrings
# -------------------------
def init_pos():
for w in time_wires:
qml.Hadamard(wires=w)
# -------------------------
# Validity compute (your parity-based oracle)
# Writes parity flags into state_wires and ANDs into good_wire
# -------------------------
def check_pos_compute():
for t in range(T_steps):
pos = pos_wires_at_t(t)
for i in range(n - 1):
bits = format(i, f"0{n_qubits_step}b")
for w, b in zip(pos, bits):
if b == "0":
qml.PauliX(wires=w)
qml.MultiControlledX(wires=pos + [state_wires[i]])
for w, b in zip(pos, bits):
if b == "0":
qml.PauliX(wires=w)
qml.MultiControlledX(wires=state_wires + [good_wire])
# -------------------------
# FULL cost oracle (Start + transitions + return)
# Uses phase kickback via toggling phase_anc
# -------------------------
def apply_transition_phase(t):
pos_t = pos_wires_at_t(t)
pos_tp1 = pos_wires_at_t(t + 1)
for i in range(n - 1):
for j in range(n - 1):
bits_i = format(i, f"0{n_qubits_step}b")
bits_j = format(j, f"0{n_qubits_step}b")
for w, b in zip(pos_t, bits_i):
if b == "0":
qml.PauliX(wires=w)
for w, b in zip(pos_tp1, bits_j):
if b == "0":
qml.PauliX(wires=w)
controls = pos_t + pos_tp1
qml.MultiControlledX(wires=controls + [phase_anc])
qml.PhaseShift(float(cost_matrix[i, j]), wires=phase_anc)
qml.MultiControlledX(wires=controls + [phase_anc])
for w, b in zip(pos_tp1, bits_j):
if b == "0":
qml.PauliX(wires=w)
for w, b in zip(pos_t, bits_i):
if b == "0":
qml.PauliX(wires=w)
def cost_oracle_full():
# Start -> first
pos0 = pos_wires_at_t(0)
for j in range(n - 1):
bits = format(j, f"0{n_qubits_step}b")
for w, b in zip(pos0, bits):
if b == "0":
qml.PauliX(wires=w)
qml.MultiControlledX(wires=pos0 + [phase_anc])
qml.PhaseShift(float(cost_matrix[start_node, j]), wires=phase_anc)
qml.MultiControlledX(wires=pos0 + [phase_anc])
for w, b in zip(pos0, bits):
if b == "0":
qml.PauliX(wires=w)
# Inner transitions
for t in range(T_steps - 1):
apply_transition_phase(t)
# Last -> start
pos_last = pos_wires_at_t(T_steps - 1)
for i in range(n - 1):
bits = format(i, f"0{n_qubits_step}b")
for w, b in zip(pos_last, bits):
if b == "0":
qml.PauliX(wires=w)
qml.MultiControlledX(wires=pos_last + [phase_anc])
qml.PhaseShift(float(cost_matrix[i, start_node]), wires=phase_anc)
qml.MultiControlledX(wires=pos_last + [phase_anc])
for w, b in zip(pos_last, bits):
if b == "0":
qml.PauliX(wires=w)
# -------------------------
# QNode: Uniform -> validity compute -> cost phase
# (No Grover, no QSVT)
# -------------------------
@qml.qnode(dev)
def circuit_baseline():
init_pos()
check_pos_compute() # marks allowed walks into good_wire (and sets state_wires)
cost_oracle_full() # encodes loss as phase
return qml.state()
state = circuit_baseline()
# -------------------------
# Decode helpers (robust to wire order)
# -------------------------
def bits_by_wire_from_index(idx):
bitstr = format(idx, f"0{len(wires)}b")
# PennyLane uses wires list ordering; map bitstring to wires robustly:
# The least-significant bit corresponds to the last wire in 'wires'.
return bitstr, {w: int(b) for w, b in zip(reversed(wires), reversed(bitstr))}
def decode_path_from_bits(bits_by_wire):
time_bits = "".join(str(bits_by_wire[w]) for w in time_wires)
path = []
for t in range(T_steps):
chunk = time_bits[t * n_qubits_step : (t + 1) * n_qubits_step]
path.append(int(chunk, 2))
return path
def is_valid_classical(path):
# valid iff it's a permutation of 0..n-2
return sorted(path) == list(range(n - 1))
def classical_cost_norm(path):
# cost in normalized matrix
cost = float(cost_matrix[start_node, path[0]])
for t in range(T_steps - 1):
cost += float(cost_matrix[path[t], path[t+1]])
cost += float(cost_matrix[path[-1], start_node])
return float(cost)
def wrap_to_pi(x):
return (x + np.pi) % (2*np.pi) - np.pi
# -------------------------
# Build table + phase check
# We use one valid state as reference to compare relative phases:
# dphi(pi) should match dcost(pi) (mod 2pi) for valid paths.
# -------------------------
def build_table_with_phase_check(statevec, amp_eps=1e-12):
ref_phase = None
ref_cost = None
# find a reference VALID state with good=1 (and phase_anc=0 ideally)
for idx, amp in enumerate(statevec):
if abs(amp) < amp_eps:
continue
_, bw = bits_by_wire_from_index(idx)
path = decode_path_from_bits(bw)
good_q = int(bw[good_wire])
valid_c = is_valid_classical(path)
phase_bit = int(bw[phase_anc])
if good_q == 1 and valid_c and phase_bit == 0:
ref_phase = float(np.angle(amp))
ref_cost = float(classical_cost_norm(path))
break
if ref_phase is None:
raise RuntimeError("No valid+good reference state found (try lowering amp_eps).")
rows = []
for idx, amp in enumerate(statevec):
if abs(amp) < amp_eps:
continue
bitstr, bw = bits_by_wire_from_index(idx)
path = decode_path_from_bits(bw)
good_q = int(bw[good_wire])
valid_c = is_valid_classical(path)
phase_bit = int(bw[phase_anc])
phase = float(np.angle(amp))
dphi = float(wrap_to_pi(phase - ref_phase))
if valid_c:
cost = float(classical_cost_norm(path))
dcost = float(cost - ref_cost)
diff_wrapped = float(abs(wrap_to_pi(dphi - dcost)))
else:
cost = np.nan
dcost = np.nan
diff_wrapped = np.nan
rows.append({
"idx": idx,
# "bitstr": bitstr,
"path": path,
"good_quantum": good_q,
"is_valid_classical": int(valid_c),
"phase_anc_bit": phase_bit,
"amp_abs": round(float(abs(amp)), 4),
"prob": round(float(abs(amp)**2), 4),
"phase": round(float(phase), 4),
"dphi": round(float(dphi), 4),
"cost_norm": round(float(cost), 4) if not np.isnan(cost) else np.nan,
"dcost": round(float(dcost), 4) if not np.isnan(dcost) else np.nan,
"diff_wrapped": round(float(diff_wrapped), 4) if not np.isnan(diff_wrapped) else np.nan,
})
return pd.DataFrame(rows)
df = build_table_with_phase_check(state)
# -------------------------
# Report + save
# -------------------------
# Compare quantum good flag vs classical validity
mismatch = df[df["good_quantum"] != df["is_valid_classical"]]
print("mismatches (good_quantum != is_valid_classical):", len(mismatch))
# Probability mass by classical validity (should be near 1 on all states because uniform;
# but note: we're summing only populated basis states; should sum to 1)
p_valid = float(df.loc[df["is_valid_classical"] == 1, "prob"].sum())
p_invalid = float(df.loc[df["is_valid_classical"] == 0, "prob"].sum())
print("Prob mass valid:", p_valid)
print("Prob mass invalid:", p_invalid)
print("Total prob mass (rows summed):", p_valid + p_invalid)
# Phase consistency on valid states
max_diff = float(np.nanmax(df.loc[df["is_valid_classical"] == 1, "diff_wrapped"].values))
print("max |wrap(dphi - dcost)| over valid states:", max_diff)
# Sort for convenience: valid & good on top, then by cost
df_sorted = df.sort_values(
["is_valid_classical", "good_quantum", "cost_norm"],
ascending=[False, False, True],
).reset_index(drop=True)
out_path = "tsp_baseline_validity_cost_phase.csv"
df_sorted.to_csv(out_path, index=False)
print("Saved:", out_path)
# optional preview
print("\nTop-20 valid states (by increasing cost_norm):")
print(df_sorted[df_sorted["is_valid_classical"] == 1][["path","good_quantum","cost_norm","phase","diff_wrapped"]].to_string(index=False))
# %%