Exercise for qosf

qosf.py

import numpy as np
import matplotlib.pyplot as plt
import time
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Operator

# Define basic quantum gates
def get_X():
    return np.array([[0, 1], [1, 0]])

def get_H():
    return np.array([[1, 1], [1, -1]]) / np.sqrt(2)

def get_I():
    return np.array([[1, 0], [0, 1]])

def get_CNOT():
    return np.array([[1, 0, 0, 0],
                     [0, 1, 0, 0],
                     [0, 0, 0, 1],
                     [0, 0, 1, 0]])

class NaiveStatevectorSimulator:
    def __init__(self, num_qubits):
        self.num_qubits = num_qubits
        # Initialize state to |0...0>
        self.state = np.zeros(2**num_qubits)
        self.state[0] = 1

    def apply_single_qubit_gate(self, gate, target_qubit):
        # Construct the full operator using kronecker products
        operator = np.array([[1]])
        for i in range(self.num_qubits):
            if i == target_qubit:
                operator = np.kron(operator, gate)
            else:
                operator = np.kron(operator, get_I())
        self.state = operator @ self.state

    def apply_CNOT(self, control, target):
        if control >= self.num_qubits or target >= self.num_qubits:
            raise ValueError("Qubit index out of range")

        # Construct CNOT matrix for the specific control and target qubits
        operator = np.eye(2**self.num_qubits)
        for i in range(2**self.num_qubits):
            binary = format(i, f'0{self.num_qubits}b')
            if binary[control] == '1':
                # Flip the target qubit
                new_binary = list(binary)
                new_binary[target] = '1' if binary[target] == '0' else '0'
                new_i = int(''.join(new_binary), 2)
                # Swap rows
                operator[i], operator[new_i] = operator[new_i].copy(), operator[i].copy()

        self.state = operator @ self.state

def measure_runtime(max_qubits):
    qubit_range = range(1, max_qubits + 1)
    times = []

    for n in qubit_range:
        start_time = time.time()

        # Create and run a test circuit
        sim = NaiveStatevectorSimulator(n)

        # Apply some test operations
        sim.apply_single_qubit_gate(get_H(), 0)
        if n > 1:
            sim.apply_CNOT(0, 1)
        sim.apply_single_qubit_gate(get_X(), 0)

        elapsed_time = time.time() - start_time
        times.append(elapsed_time)

        print(f"Completed simulation for {n} qubits in {elapsed_time:.4f} seconds")

    return qubit_range, times

# Run the benchmark
max_qubits = 10  # Adjust this based on your computer's capabilities
qubit_range, times = measure_runtime(max_qubits)

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(qubit_range, times, 'bo-')
plt.yscale('log')
plt.xlabel('Number of Qubits')
plt.ylabel('Runtime (seconds)')
plt.title('Quantum Circuit Simulation Runtime vs Number of Qubits')
plt.grid(True)
plt.show()

# Test the simulator with a simple circuit
test_sim = NaiveStatevectorSimulator(2)
print("\nInitial state:", test_sim.state)

# Apply Hadamard to first qubit
test_sim.apply_single_qubit_gate(get_H(), 0)
print("After H on qubit 0:", test_sim.state)

# Apply CNOT
test_sim.apply_CNOT(0, 1)
print("After CNOT:", test_sim.state)

# Apply X to first qubit
test_sim.apply_single_qubit_gate(get_X(), 0)
print("Final state:", test_sim.state)