In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from qiskit import QuantumCircuit
from qiskit.circuit.library import TwoLocal #based
from qiskit.circuit.library import EfficientSU2 #I used this, for this example is the best
from qiskit.circuit.library.n_local import n_local #In qiskit v3 it replaces TwoLocal or EfficientSU2
from qiskit.quantum_info import Statevector, SparsePauliOp
from qiskit_algorithms.optimizers import SPSA
from scipy.optimize import minimize #minimize or SPSA
from matplotlib.lines import Line2D
import warnings
warnings.filterwarnings('ignore', category=DeprecationWarning) #For TwoLocal and EfficientSU2
np.random.seed(42)
N = 6 # number of transmitters
positions = np.random.rand(N, 2) * 80 # 80×80 meters
freq_count = 3
num_qubits = N * freq_count # 18 qubits (one-hot encoding)
lambda_penalty = 30.0 # penalty for invalid assignment
interf_range = 30.0 # maximum interference range
p = 3 # Increased depth for better results
n_iters = 180
print(f"Number of transmitters: {N}")
print(f"Number of frequencies: {freq_count}")
print(f"Total qubits: {num_qubits}")
print(f"\nTransmitter positions:")
for i, pos in enumerate(positions):
print(f" Transmitter {i}: x={pos[0]:.2f}m, y={pos[1]:.2f}m")
Number of transmitters: 6 Number of frequencies: 3 Total qubits: 18 Transmitter positions: Transmitter 0: x=29.96m, y=76.06m Transmitter 1: x=58.56m, y=47.89m Transmitter 2: x=12.48m, y=12.48m Transmitter 3: x=4.65m, y=69.29m Transmitter 4: x=48.09m, y=56.65m Transmitter 5: x=1.65m, y=77.59m
In [2]:
def draw_transmitter(ax, x, y, color, label, size=8):
"""Draw transmitter as triangle with antenna"""
# Antenna (line)
ax.plot([x, x], [y, y + size*0.8], color='black', linewidth=3, zorder=5)
# Base (triangle)
triangle = patches.RegularPolygon((x, y), 3, radius=size*0.6,
orientation=0,
facecolor=color,
edgecolor='black',
linewidth=2,
zorder=4)
ax.add_patch(triangle)
# Label
ax.text(x, y-size*0.3, label, fontsize=12, color='blue',
ha='center', va='center', weight='bold', zorder=6)
fig, ax = plt.subplots(figsize=(10, 10))
for i, (x, y) in enumerate(positions):
draw_transmitter(ax, x, y, '#D3D3D3', str(i)) # Pastel gray
circle = plt.Circle((x, y), 50, color='#B4E7F5', alpha=0.3, zorder=1) # Pastel blue
ax.add_patch(circle)
ax.set_title('Transmitter Positions (before frequency assignment)', fontsize=14, weight='bold')
ax.set_xlabel('Distance [m]', fontsize=12)
ax.set_ylabel('Distance [m]', fontsize=12)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.legend(['Transmitter', 'Interference zone (r=50m)'], loc='upper right')
plt.tight_layout()
plt.show()
distance_matrix = np.zeros((N, N))
for i in range(N):
for j in range(N):
distance_matrix[i, j] = np.linalg.norm(positions[i] - positions[j])
print("Distances between transmitters [m]:")
print(" ", end="")
for j in range(N):
print(f"T{j:2d} ", end="")
print()
for i in range(N):
print(f"T{i:2d}: ", end="")
for j in range(N):
if i == j:
print(" - ", end="")
else:
print(f"{distance_matrix[i,j]:5.1f} ", end="")
print()
Distances between transmitters [m]:
T 0 T 1 T 2 T 3 T 4 T 5
T 0: - 40.1 65.9 26.2 26.6 28.4
T 1: 40.1 - 58.1 58.0 13.6 64.2
T 2: 65.9 58.1 - 57.4 56.7 66.0
T 3: 26.2 58.0 57.4 - 45.2 8.8
T 4: 26.6 13.6 56.7 45.2 - 50.9
T 5: 28.4 64.2 66.0 8.8 50.9 -
In [3]:
pauli_terms = []
print(f"Penalty coefficient λ = {lambda_penalty}")
print(f"Interference range = {interf_range}m")
# 1. ONE-HOT CONSTRAINT: (Σ_k Z_k - 1)^2
# For one-hot encoding (one |1⟩): we want Σ Z_k = +1
# Since Z|0⟩ = +1, Z|1⟩ = -1:
# - {|1⟩,|0⟩,|0⟩}: (-1)+(+1)+(+1) = +1 ✓
# - {|0⟩,|1⟩,|0⟩}: (+1)+(-1)+(+1) = +1 ✓
# - {|0⟩,|0⟩,|1⟩}: (+1)+(+1)+(-1) = +1 ✓
#
# Expand (Z_0 + Z_1 + Z_2 - 1)^2:
# = Z_0^2 + Z_1^2 + Z_2^2 + 2*Z_0*Z_1 + 2*Z_0*Z_2 + 2*Z_1*Z_2 - 2*Z_0 - 2*Z_1 - 2*Z_2 + 1
# Since Z^2 = I:
# = 3*I + 2*Z_0*Z_1 + 2*Z_0*Z_2 + 2*Z_1*Z_2 - 2*Z_0 - 2*Z_1 - 2*Z_2 + 1
# = 4*I + 2*Z_0*Z_1 + 2*Z_0*Z_2 + 2*Z_1*Z_2 - 2*Z_0 - 2*Z_1 - 2*Z_2
print("Exactly one bit |1⟩ per transmitter")
print("Method: (Σ Z - 1)^2 with CORRECT signs")
constraint_count = 0
for i in range(N):
base = i * freq_count
# Quadratic terms Z_i * Z_j: +2*lambda (cross terms penalize multiple |1⟩)
for k1 in range(freq_count):
for k2 in range(k1+1, freq_count):
label = ['I'] * num_qubits
label[base + k1] = 'Z'
label[base + k2] = 'Z'
pauli_string = ''.join(label)
pauli_terms.append((pauli_string, lambda_penalty * 2.0))
constraint_count += 1
# Linear terms Z_i: -2*lambda (prefer |1⟩ since Z|1⟩=-1)
for k in range(freq_count):
label = ['I'] * num_qubits
label[base + k] = 'Z'
pauli_string = ''.join(label)
pauli_terms.append((pauli_string, lambda_penalty * (-2.0))) # NEGATIVE!
constraint_count += 1
# Constant term: +4*lambda (offset)
label = ['I'] * num_qubits
pauli_string = ''.join(label)
pauli_terms.append((pauli_string, lambda_penalty * 4.0))
constraint_count += 1
if i == 0: # Example for first transmitter
print(f" Transmitter {i}: 3 cross(+) + 3 linear(-) + 1 offset = 7 terms")
print(f"Total constraint terms: {constraint_count}")
#INTERFERENCE TERMS
interference_count = 0
interference_details = []
for i in range(N):
for j in range(i+1, N):
dist = distance_matrix[i, j]
if dist >= interf_range:
continue
loss_db = 20 * (interf_range / dist)
for f in range(freq_count):
p1 = i * freq_count + f
p2 = j * freq_count + f
label = ['I'] * num_qubits
label[p1] = 'Z'
label[p2] = 'Z'
pauli_string = ''.join(label)
coef = loss_db / 4.0
pauli_terms.append((pauli_string, coef))
interference_count += 1
interference_details.append({
'tx1': i, 'tx2': j, 'freq': f,
'dist': dist, 'loss': loss_db, 'coef': coef
})
print(f"Total interference terms: {interference_count}")
print(f"\nExamples:")
for idx, detail in enumerate(interference_details):
print(f" Transmitter {detail['tx1']} <-> {detail['tx2']}, "
f"f={detail['freq']}: dist={detail['dist']:.1f}m, "
f"loss={detail['loss']:.1f}dB, coef={detail['coef']:.2f}")
# Create Hamiltonian
H = SparsePauliOp.from_list(pauli_terms)
print(f"\nHamiltonian created: {len(pauli_terms)} Pauli terms")
Penalty coefficient λ = 30.0 Interference range = 30.0m Exactly one bit |1⟩ per transmitter Method: (Σ Z - 1)^2 with CORRECT signs Transmitter 0: 3 cross(+) + 3 linear(-) + 1 offset = 7 terms Total constraint terms: 42 Total interference terms: 15 Examples: Transmitter 0 <-> 3, f=0: dist=26.2m, loss=22.9dB, coef=5.72 Transmitter 0 <-> 3, f=1: dist=26.2m, loss=22.9dB, coef=5.72 Transmitter 0 <-> 3, f=2: dist=26.2m, loss=22.9dB, coef=5.72 Transmitter 0 <-> 4, f=0: dist=26.6m, loss=22.6dB, coef=5.65 Transmitter 0 <-> 4, f=1: dist=26.6m, loss=22.6dB, coef=5.65 Transmitter 0 <-> 4, f=2: dist=26.6m, loss=22.6dB, coef=5.65 Transmitter 0 <-> 5, f=0: dist=28.4m, loss=21.2dB, coef=5.29 Transmitter 0 <-> 5, f=1: dist=28.4m, loss=21.2dB, coef=5.29 Transmitter 0 <-> 5, f=2: dist=28.4m, loss=21.2dB, coef=5.29 Transmitter 1 <-> 4, f=0: dist=13.6m, loss=44.0dB, coef=10.99 Transmitter 1 <-> 4, f=1: dist=13.6m, loss=44.0dB, coef=10.99 Transmitter 1 <-> 4, f=2: dist=13.6m, loss=44.0dB, coef=10.99 Transmitter 3 <-> 5, f=0: dist=8.8m, loss=68.0dB, coef=17.00 Transmitter 3 <-> 5, f=1: dist=8.8m, loss=68.0dB, coef=17.00 Transmitter 3 <-> 5, f=2: dist=8.8m, loss=68.0dB, coef=17.00 Hamiltonian created: 57 Pauli terms
In [4]:
print(f"Total: {constraint_count} terms")
print("\nFirst transmitter (T0) - all 7 terms:")
constraint_idx = 0
for i in range(1): # Only first transmitter for detailed view
base = i * freq_count
print(f"\nTransmitter {i} constraint terms:")
# Cross terms
for k1 in range(freq_count):
for k2 in range(k1+1, freq_count):
pauli_word, coef = pauli_terms[constraint_idx]
highlighted = [f"q{base+k1}:Z", f"q{base+k2}:Z"]
print(f" [{constraint_idx:2d}] {pauli_word} | coef={coef:6.1f} | Cross: {highlighted}")
constraint_idx += 1
# Linear terms
for k in range(freq_count):
pauli_word, coef = pauli_terms[constraint_idx]
highlighted = [f"q{base+k}:Z"]
print(f" [{constraint_idx:2d}] {pauli_word} | coef={coef:6.1f} | Linear: {highlighted}")
constraint_idx += 1
# Constant term
pauli_word, coef = pauli_terms[constraint_idx]
print(f" [{constraint_idx:2d}] {pauli_word} | coef={coef:6.1f} | Offset (Identity)")
constraint_idx += 1
# Display ALL interference terms
print(f"Total: {interference_count} terms")
print("\nAll interference terms (full listing):")
interf_start_idx = constraint_count
for idx, detail in enumerate(interference_details):
pauli_idx = interf_start_idx + idx
pauli_word, coef = pauli_terms[pauli_idx]
q1 = detail['tx1'] * freq_count + detail['freq']
q2 = detail['tx2'] * freq_count + detail['freq']
print(f" [{pauli_idx:2d}] T{detail['tx1']} <-> T{detail['tx2']} f{detail['freq']}: "
f"q{q1}:Z x q{q2}:Z | "
f"dist={detail['dist']:4.1f}m, loss={detail['loss']:5.1f}dB, coef={coef:5.2f}")
if idx < len(interference_details) - 1:
# Show the actual Pauli word every 5 terms for verification
if idx % 5 == 4:
print(f" Pauli word: {pauli_word}")
# Matrix visualization of Pauli operators
print("Legend: Z=Pauli-Z, I=Identity\n")
# Create a visual matrix for first few terms
print("Constraint terms (T0 example):")
print(" ", end="")
for q in range(num_qubits):
print(f"q{q:2d}", end=" ")
print()
for idx in range(min(7, len(pauli_terms))):
pauli_word, coef = pauli_terms[idx]
print(f"[{idx:2d}] ", end="")
for q, char in enumerate(pauli_word):
if char == 'Z':
print(f" Z ", end=" ")
else:
print(f" · ", end=" ")
print(f"| coef={coef:6.1f}")
print("\nInterference terms (first 10):")
print(" ", end="")
for q in range(num_qubits):
print(f"q{q:2d}", end=" ")
print()
for idx in range(constraint_count, min(constraint_count + 10, len(pauli_terms))):
pauli_word, coef = pauli_terms[idx]
print(f"[{idx:2d}] ", end="")
for q, char in enumerate(pauli_word):
if char == 'Z':
print(f" Z ", end=" ")
else:
print(f" · ", end=" ")
print(f"| coef={coef:6.1f}")
# Summary statistics
z_counts = []
for pauli_word, _ in pauli_terms:
z_count = pauli_word.count('Z')
z_counts.append(z_count)
print(f"Total Pauli terms: {len(pauli_terms)}")
print(f"Average Z operators per term: {np.mean(z_counts):.2f}")
print(f"Terms with 1 Z: {z_counts.count(1)} (linear)")
print(f"Terms with 2 Z: {z_counts.count(2)} (quadratic)")
print(f"Terms with 0 Z: {z_counts.count(0)} (constant)")
# Export all Pauli words to file
with open('pauli_words_complete.txt', 'w') as f:
f.write("CONSTRAINT TERMS:\n")
for idx in range(constraint_count):
pauli_word, coef = pauli_terms[idx]
f.write(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f}\n")
print(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f}\n")
f.write("\n\nINTERFERENCE TERMS:\n")
for idx in range(constraint_count, len(pauli_terms)):
pauli_word, coef = pauli_terms[idx]
detail_idx = idx - constraint_count
if detail_idx < len(interference_details):
detail = interference_details[detail_idx]
f.write(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f} | "
f"T{detail['tx1']}<->T{detail['tx2']} f{detail['freq']}\n")
print(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f} | "
f"T{detail['tx1']}<->T{detail['tx2']} f{detail['freq']}\n")
else:
f.write(f"[{idx:3d}] {pauli_word} | coef={coef:8.2f}\n")
Total: 42 terms
First transmitter (T0) - all 7 terms:
Transmitter 0 constraint terms:
[ 0] ZZIIIIIIIIIIIIIIII | coef= 60.0 | Cross: ['q0:Z', 'q1:Z']
[ 1] ZIZIIIIIIIIIIIIIII | coef= 60.0 | Cross: ['q0:Z', 'q2:Z']
[ 2] IZZIIIIIIIIIIIIIII | coef= 60.0 | Cross: ['q1:Z', 'q2:Z']
[ 3] ZIIIIIIIIIIIIIIIII | coef= -60.0 | Linear: ['q0:Z']
[ 4] IZIIIIIIIIIIIIIIII | coef= -60.0 | Linear: ['q1:Z']
[ 5] IIZIIIIIIIIIIIIIII | coef= -60.0 | Linear: ['q2:Z']
[ 6] IIIIIIIIIIIIIIIIII | coef= 120.0 | Offset (Identity)
Total: 15 terms
All interference terms (full listing):
[42] T0 <-> T3 f0: q0:Z x q9:Z | dist=26.2m, loss= 22.9dB, coef= 5.72
[43] T0 <-> T3 f1: q1:Z x q10:Z | dist=26.2m, loss= 22.9dB, coef= 5.72
[44] T0 <-> T3 f2: q2:Z x q11:Z | dist=26.2m, loss= 22.9dB, coef= 5.72
[45] T0 <-> T4 f0: q0:Z x q12:Z | dist=26.6m, loss= 22.6dB, coef= 5.65
[46] T0 <-> T4 f1: q1:Z x q13:Z | dist=26.6m, loss= 22.6dB, coef= 5.65
Pauli word: IZIIIIIIIIIIIZIIII
[47] T0 <-> T4 f2: q2:Z x q14:Z | dist=26.6m, loss= 22.6dB, coef= 5.65
[48] T0 <-> T5 f0: q0:Z x q15:Z | dist=28.4m, loss= 21.2dB, coef= 5.29
[49] T0 <-> T5 f1: q1:Z x q16:Z | dist=28.4m, loss= 21.2dB, coef= 5.29
[50] T0 <-> T5 f2: q2:Z x q17:Z | dist=28.4m, loss= 21.2dB, coef= 5.29
[51] T1 <-> T4 f0: q3:Z x q12:Z | dist=13.6m, loss= 44.0dB, coef=10.99
Pauli word: IIIZIIIIIIIIZIIIII
[52] T1 <-> T4 f1: q4:Z x q13:Z | dist=13.6m, loss= 44.0dB, coef=10.99
[53] T1 <-> T4 f2: q5:Z x q14:Z | dist=13.6m, loss= 44.0dB, coef=10.99
[54] T3 <-> T5 f0: q9:Z x q15:Z | dist= 8.8m, loss= 68.0dB, coef=17.00
[55] T3 <-> T5 f1: q10:Z x q16:Z | dist= 8.8m, loss= 68.0dB, coef=17.00
[56] T3 <-> T5 f2: q11:Z x q17:Z | dist= 8.8m, loss= 68.0dB, coef=17.00
Legend: Z=Pauli-Z, I=Identity
Constraint terms (T0 example):
q 0 q 1 q 2 q 3 q 4 q 5 q 6 q 7 q 8 q 9 q10 q11 q12 q13 q14 q15 q16 q17
[ 0] Z Z · · · · · · · · · · · · · · · · | coef= 60.0
[ 1] Z · Z · · · · · · · · · · · · · · · | coef= 60.0
[ 2] · Z Z · · · · · · · · · · · · · · · | coef= 60.0
[ 3] Z · · · · · · · · · · · · · · · · · | coef= -60.0
[ 4] · Z · · · · · · · · · · · · · · · · | coef= -60.0
[ 5] · · Z · · · · · · · · · · · · · · · | coef= -60.0
[ 6] · · · · · · · · · · · · · · · · · · | coef= 120.0
Interference terms (first 10):
q 0 q 1 q 2 q 3 q 4 q 5 q 6 q 7 q 8 q 9 q10 q11 q12 q13 q14 q15 q16 q17
[42] Z · · · · · · · · Z · · · · · · · · | coef= 5.7
[43] · Z · · · · · · · · Z · · · · · · · | coef= 5.7
[44] · · Z · · · · · · · · Z · · · · · · | coef= 5.7
[45] Z · · · · · · · · · · · Z · · · · · | coef= 5.6
[46] · Z · · · · · · · · · · · Z · · · · | coef= 5.6
[47] · · Z · · · · · · · · · · · Z · · · | coef= 5.6
[48] Z · · · · · · · · · · · · · · Z · · | coef= 5.3
[49] · Z · · · · · · · · · · · · · · Z · | coef= 5.3
[50] · · Z · · · · · · · · · · · · · · Z | coef= 5.3
[51] · · · Z · · · · · · · · Z · · · · · | coef= 11.0
Total Pauli terms: 57
Average Z operators per term: 1.47
Terms with 1 Z: 18 (linear)
Terms with 2 Z: 33 (quadratic)
Terms with 0 Z: 6 (constant)
[ 0] ZZIIIIIIIIIIIIIIII | coef= 60.00
[ 1] ZIZIIIIIIIIIIIIIII | coef= 60.00
[ 2] IZZIIIIIIIIIIIIIII | coef= 60.00
[ 3] ZIIIIIIIIIIIIIIIII | coef= -60.00
[ 4] IZIIIIIIIIIIIIIIII | coef= -60.00
[ 5] IIZIIIIIIIIIIIIIII | coef= -60.00
[ 6] IIIIIIIIIIIIIIIIII | coef= 120.00
[ 7] IIIZZIIIIIIIIIIIII | coef= 60.00
[ 8] IIIZIZIIIIIIIIIIII | coef= 60.00
[ 9] IIIIZZIIIIIIIIIIII | coef= 60.00
[ 10] IIIZIIIIIIIIIIIIII | coef= -60.00
[ 11] IIIIZIIIIIIIIIIIII | coef= -60.00
[ 12] IIIIIZIIIIIIIIIIII | coef= -60.00
[ 13] IIIIIIIIIIIIIIIIII | coef= 120.00
[ 14] IIIIIIZZIIIIIIIIII | coef= 60.00
[ 15] IIIIIIZIZIIIIIIIII | coef= 60.00
[ 16] IIIIIIIZZIIIIIIIII | coef= 60.00
[ 17] IIIIIIZIIIIIIIIIII | coef= -60.00
[ 18] IIIIIIIZIIIIIIIIII | coef= -60.00
[ 19] IIIIIIIIZIIIIIIIII | coef= -60.00
[ 20] IIIIIIIIIIIIIIIIII | coef= 120.00
[ 21] IIIIIIIIIZZIIIIIII | coef= 60.00
[ 22] IIIIIIIIIZIZIIIIII | coef= 60.00
[ 23] IIIIIIIIIIZZIIIIII | coef= 60.00
[ 24] IIIIIIIIIZIIIIIIII | coef= -60.00
[ 25] IIIIIIIIIIZIIIIIII | coef= -60.00
[ 26] IIIIIIIIIIIZIIIIII | coef= -60.00
[ 27] IIIIIIIIIIIIIIIIII | coef= 120.00
[ 28] IIIIIIIIIIIIZZIIII | coef= 60.00
[ 29] IIIIIIIIIIIIZIZIII | coef= 60.00
[ 30] IIIIIIIIIIIIIZZIII | coef= 60.00
[ 31] IIIIIIIIIIIIZIIIII | coef= -60.00
[ 32] IIIIIIIIIIIIIZIIII | coef= -60.00
[ 33] IIIIIIIIIIIIIIZIII | coef= -60.00
[ 34] IIIIIIIIIIIIIIIIII | coef= 120.00
[ 35] IIIIIIIIIIIIIIIZZI | coef= 60.00
[ 36] IIIIIIIIIIIIIIIZIZ | coef= 60.00
[ 37] IIIIIIIIIIIIIIIIZZ | coef= 60.00
[ 38] IIIIIIIIIIIIIIIZII | coef= -60.00
[ 39] IIIIIIIIIIIIIIIIZI | coef= -60.00
[ 40] IIIIIIIIIIIIIIIIIZ | coef= -60.00
[ 41] IIIIIIIIIIIIIIIIII | coef= 120.00
[ 42] ZIIIIIIIIZIIIIIIII | coef= 5.72 | T0<->T3 f0
[ 43] IZIIIIIIIIZIIIIIII | coef= 5.72 | T0<->T3 f1
[ 44] IIZIIIIIIIIZIIIIII | coef= 5.72 | T0<->T3 f2
[ 45] ZIIIIIIIIIIIZIIIII | coef= 5.65 | T0<->T4 f0
[ 46] IZIIIIIIIIIIIZIIII | coef= 5.65 | T0<->T4 f1
[ 47] IIZIIIIIIIIIIIZIII | coef= 5.65 | T0<->T4 f2
[ 48] ZIIIIIIIIIIIIIIZII | coef= 5.29 | T0<->T5 f0
[ 49] IZIIIIIIIIIIIIIIZI | coef= 5.29 | T0<->T5 f1
[ 50] IIZIIIIIIIIIIIIIIZ | coef= 5.29 | T0<->T5 f2
[ 51] IIIZIIIIIIIIZIIIII | coef= 10.99 | T1<->T4 f0
[ 52] IIIIZIIIIIIIIZIIII | coef= 10.99 | T1<->T4 f1
[ 53] IIIIIZIIIIIIIIZIII | coef= 10.99 | T1<->T4 f2
[ 54] IIIIIIIIIZIIIIIZII | coef= 17.00 | T3<->T5 f0
[ 55] IIIIIIIIIIZIIIIIZI | coef= 17.00 | T3<->T5 f1
[ 56] IIIIIIIIIIIZIIIIIZ | coef= 17.00 | T3<->T5 f2
In [10]:
ansatz = TwoLocal(num_qubits, 'ry', 'cz', entanglement='linear', reps=p)
#ansatz = EfficientSU2(num_qubits, entanglement='linear', reps=p)
#ansatz = n_local(num_qubits=num_qubits, rotation_blocks='ry', entanglement_blocks='cz', entanglement='linear', reps=p)
print(f"QAOA depth: p={p}")
print(f"Number of parameters: {ansatz.num_parameters}")
def qaoa_circuit(params):
qc = QuantumCircuit(num_qubits)
qc.h(range(num_qubits))
bound = ansatz.assign_parameters(params)
qc.compose(bound, inplace=True)
return qc
iteration_data = []
def cost_function(params):
qc = qaoa_circuit(params)
state = Statevector(qc)
energy = state.expectation_value(H).real
iteration_data.append({
'params': params.copy(),
'energy': energy,
'state': state
})
if len(iteration_data) % 10 == 0:
print(f" Iteration {len(iteration_data):3d}: energy = {energy:8.2f}")
return energy
initial_params = np.random.uniform(0, 2*np.pi, ansatz.num_parameters)
#result = minimize(cost_function, initial_params, method='COBYLA',
# options={'maxiter': n_iters, 'disp': False})
#optimizer = SPSA(maxiter=n_iters, learning_rate=0.01, perturbation=0.1)
optimizer = SPSA(maxiter=n_iters)
result = optimizer.minimize(fun=cost_function, x0=initial_params)
optimal_params = result.x
final_energy = result.fun
print(f" Number of iterations: {len(iteration_data)}")
print(f" Best energy: {final_energy:.2f}")
print(f" Initial energy: {iteration_data[0]['energy']:.2f}")
print(f" Improvement: {iteration_data[0]['energy'] - final_energy:.2f}")
QAOA depth: p=3 Number of parameters: 72 Iteration 10: energy = 489.32 Iteration 20: energy = 595.88 Iteration 30: energy = 608.86 Iteration 40: energy = 588.00 Iteration 50: energy = 638.37 Iteration 60: energy = 499.08 Iteration 70: energy = 347.95 Iteration 80: energy = 328.27 Iteration 90: energy = 371.17 Iteration 100: energy = 373.88 Iteration 110: energy = 277.79 Iteration 120: energy = 264.09 Iteration 130: energy = 304.44 Iteration 140: energy = 263.08 Iteration 150: energy = 239.07 Iteration 160: energy = 226.18 Iteration 170: energy = 230.14 Iteration 180: energy = 188.94 Iteration 190: energy = 226.49 Iteration 200: energy = 223.56 Iteration 210: energy = 190.62 Iteration 220: energy = 158.19 Iteration 230: energy = 163.64 Iteration 240: energy = 155.37 Iteration 250: energy = 142.08 Iteration 260: energy = 147.37 Iteration 270: energy = 118.18 Iteration 280: energy = 121.78 Iteration 290: energy = 130.40 Iteration 300: energy = 131.84 Iteration 310: energy = 116.72 Iteration 320: energy = 113.51 Iteration 330: energy = 133.30 Iteration 340: energy = 113.35 Iteration 350: energy = 90.81 Iteration 360: energy = 94.29 Iteration 370: energy = 83.99 Iteration 380: energy = 75.64 Iteration 390: energy = 70.22 Iteration 400: energy = 84.29 Iteration 410: energy = 55.27 Number of iterations: 411 Best energy: 38.72 Initial energy: 622.66 Improvement: 583.94
In [11]:
energies = [d['energy'] for d in iteration_data]
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(energies, linewidth=2, color='#A8DADC') # blue
plt.axhline(y=final_energy, color='#E63946', linestyle='--', linewidth=2,
label=f'Minimum: {final_energy:.2f}') # red
plt.xlabel('Iteration', fontsize=12)
plt.ylabel('Energy', fontsize=12)
plt.title('Optimization Convergence', fontsize=14, weight='bold')
plt.grid(True, alpha=0.3)
plt.legend()
plt.subplot(1, 2, 2)
improvements = np.abs(np.array(energies) - final_energy) + 1e-10
plt.semilogy(improvements, linewidth=2, color='#90BE6D') # green
plt.xlabel('Iteration', fontsize=12)
plt.ylabel('Distance from minimum', fontsize=12)
plt.title('Convergence (log scale)', fontsize=14, weight='bold')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
final_qc = qaoa_circuit(optimal_params)
final_state = Statevector(final_qc)
probs = final_state.probabilities()
top_k = 10
top_indices = np.argsort(probs)[-top_k:][::-1]
print(f"\nTop {top_k} most probable states:")
for rank, idx in enumerate(top_indices, 1):
binary = format(idx, f'0{num_qubits}b')
prob = probs[idx]
assignment_str = []
for i in range(N):
chunk = binary[i*3:(i+1)*3]
assignment_str.append(f"T{i}:{chunk}")
print(f" {rank:2d}. |{binary}⟩ P={prob:.4f} ({prob*100:.2f}%)")
print(f" → {' '.join(assignment_str)}")
# Probability distribution visualization
plt.figure(figsize=(14, 5))
plt.subplot(1, 2, 1)
top_50 = np.argsort(probs)[-50:][::-1]
plt.bar(range(len(top_50)), probs[top_50], color='#A8DADC', edgecolor='#457B9D', linewidth=1.5) # blue
plt.xlabel('State rank', fontsize=12)
plt.ylabel('Probability', fontsize=12)
plt.title('Top 50 most probable states', fontsize=14, weight='bold')
plt.grid(True, alpha=0.3, axis='y')
plt.subplot(1, 2, 2)
plt.hist(probs[probs > 1e-6], bins=50, color='#F4A261', edgecolor='#E76F51', linewidth=1.5) # orange
plt.xlabel('Probability', fontsize=12)
plt.ylabel('Number of states', fontsize=12)
plt.title('Probability distribution', fontsize=14, weight='bold')
plt.yscale('log')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Top 10 most probable states:
1. |100100100001010010⟩ P=0.2588 (25.88%)
→ T0:100 T1:100 T2:100 T3:001 T4:010 T5:010
2. |100100010001010010⟩ P=0.2023 (20.23%)
→ T0:100 T1:100 T2:010 T3:001 T4:010 T5:010
3. |100000100001010010⟩ P=0.0576 (5.76%)
→ T0:100 T1:000 T2:100 T3:001 T4:010 T5:010
4. |100100100001010001⟩ P=0.0531 (5.31%)
→ T0:100 T1:100 T2:100 T3:001 T4:010 T5:001
5. |100000010001010010⟩ P=0.0450 (4.50%)
→ T0:100 T1:000 T2:010 T3:001 T4:010 T5:010
6. |100100010001010001⟩ P=0.0415 (4.15%)
→ T0:100 T1:100 T2:010 T3:001 T4:010 T5:001
7. |101100100001010010⟩ P=0.0291 (2.91%)
→ T0:101 T1:100 T2:100 T3:001 T4:010 T5:010
8. |100100100001100010⟩ P=0.0259 (2.59%)
→ T0:100 T1:100 T2:100 T3:001 T4:100 T5:010
9. |101100010001010010⟩ P=0.0227 (2.27%)
→ T0:101 T1:100 T2:010 T3:001 T4:010 T5:010
10. |100100010001100010⟩ P=0.0202 (2.02%)
→ T0:100 T1:100 T2:010 T3:001 T4:100 T5:010
In [12]:
best_idx = np.argmax(probs)
best_bin = format(best_idx, f'0{num_qubits}b')
print(f"\nMost probable state: |{best_bin}⟩")
print(f"Probability: {probs[best_idx]:.4f} ({probs[best_idx]*100:.2f}%)")
assignment = []
print("\nDecoding one-hot encoding:")
for i in range(N):
chunk = best_bin[i*3:(i+1)*3]
print(f" Transmitter {i}: {chunk}", end="")
freq = chunk.find('1')
if freq == -1:
print(f" INVALID (no bit=1), fallback → f0")
freq = 0
else:
ones_count = chunk.count('1')
if ones_count > 1:
print(f" INVALID ({ones_count} bits=1), fallback → f{freq}")
else:
print(f" frequency {freq}")
assignment.append(freq)
colors = ['#FFB3BA', '#BAFFC9', '#BAE1FF'] # red, green, blue
freq_names = ['f0 (red)', 'f1 (green)', 'f2 (blue)']
print("\nAssigned frequencies:")
for i, f in enumerate(assignment):
print(f" Transmitter {i} → {freq_names[f]}")
# Check collisions
print("\nInterference collision check:")
collision_count = 0
for i in range(N):
for j in range(i+1, N):
if assignment[i] == assignment[j]:
dist = distance_matrix[i, j]
if dist < interf_range:
collision_count += 1
loss = 20 * (interf_range / dist)
print(f"COLLISION: Transmitter {i} <-> {j} on {freq_names[assignment[i]]}, "
f"dist={dist:.1f}m, loss={loss:.1f}dB")
if collision_count == 0:
print("No interference collisions!")
else:
print(f"Total collisions: {collision_count}")
Most probable state: |100100100001010010⟩ Probability: 0.2588 (25.88%) Decoding one-hot encoding: Transmitter 0: 100 frequency 0 Transmitter 1: 100 frequency 0 Transmitter 2: 100 frequency 0 Transmitter 3: 001 frequency 2 Transmitter 4: 010 frequency 1 Transmitter 5: 010 frequency 1 Assigned frequencies: Transmitter 0 → f0 (red) Transmitter 1 → f0 (red) Transmitter 2 → f0 (red) Transmitter 3 → f2 (blue) Transmitter 4 → f1 (green) Transmitter 5 → f1 (green) Interference collision check: No interference collisions!
In [13]:
freq_colors = [colors[f] for f in assignment]
fig, ax = plt.subplots(figsize=(12, 10))
# Interference lines
for i in range(N):
for j in range(i+1, N):
dist = distance_matrix[i, j]
if dist < interf_range and assignment[i] == assignment[j]:
# Collision - red line
ax.plot([positions[i,0], positions[j,0]],
[positions[i,1], positions[j,1]],
color='#FF6B6B', linewidth=3, alpha=0.7, zorder=2)
elif dist < interf_range:
# In range but different frequencies - green line
ax.plot([positions[i,0], positions[j,0]],
[positions[i,1], positions[j,1]],
color='#95E1D3', linewidth=1.5, linestyle='--', alpha=0.5, zorder=1)
# Transmitters
for i, (x, y) in enumerate(positions):
draw_transmitter(ax, x, y, freq_colors[i], str(i), size=6)
circle = plt.Circle((x, y), interf_range,
color=freq_colors[i], alpha=0.3, zorder=0)
ax.add_patch(circle)
legend_elements = [
Line2D([0], [0], marker='^', color='w', label=f'Frequency {i} ({freq_names[i].split()[0]})',
markerfacecolor=colors[i], markersize=15, markeredgecolor='black', markeredgewidth=2)
for i in range(freq_count)
]
legend_elements.extend([
Line2D([0], [0], color='#FF6B6B', linewidth=3, label='Interference (collision)'),
Line2D([0], [0], color='#95E1D3', linewidth=2, linestyle='--', label='In range (OK)')
])
ax.legend(handles=legend_elements, loc='upper left', fontsize=11)
ax.set_title('Result: Frequency Assignment', fontsize=16, weight='bold')
ax.set_xlabel('Distance [m]', fontsize=12)
ax.set_ylabel('Distance [m]', fontsize=12)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
