<jupyter_start><jupyter_text>Quantum Circuit Simulator From Scratch<NAME>, March 2021.In this notebook I aim to implement very simple simulator of Quantum Circuits using only NumPy, which supports all basic functionality (qubits, gates and measurements).I don't intend it to be used by other people, because there are already hunderds of such simulators. Goal is to check how easy it is to implement such a simulator from scratch. I just told someone that it's very easy and even trivial, but I realized that I haven't done that myself.Also I don't intend to make it efficient, because anyway simulation is exponetially hard. So every operation will have $O(2^N)$ complexity.**Model**Our simulator will simulate a register with $N$ qubits. It will be explicitly represented in memory by a state vector, i.e. normalized complex vector with $2^N$ elements.The most confusing part for me, when working with qubits, is to have consistent "[endianness](https://en.wikipedia.org/wiki/Endianness)" - that is, whether our register is big-endian or little-endian. In this project it will be always little-endian: qubit with index 0 is encoded by least significant bit in state index. So, for $N=2$, vector $(a_0, a_1, a_2, a_3)^T$ encodes state $a_0 |0 \rangle  |0 \rangle + a_0 |1 \rangle  |0 \rangle + a_2 |0 \rangle  |1 \rangle + a_3 |1 \rangle  |1 \rangle$, and qubits in tensor product are numbered from left to right, starting from 0.**Operations**There are two fundamental operations which quantum computer must perform: apllying gate ad measurement.**Gate.** To apply gate, specified by matrix and list of qubits, we will first represent this gate as $2^N \times 2^N$ unitary matrix acting on entire state vector and then just multipy it by the vector.**Measurement.** To apply measurement on given qubit (with index i) we will first respresent state as $| \psi \rangle = | 0 \rangle_i \otimes |\psi_0 \rangle + | 1 \rangle_i \otimes |\psi_1 \rangle$ (this is splitting $2^N$-dmesional vector into two $2^{N-1}$-dimesnional vectors. Then define probabilities $p_0 = \langle \psi_0 | \psi_0 \rangle, p_1 = \langle \psi_1 | \psi_1 \rangle$. Then with probability $p_0$ replace state with $  | 0 \rangle_i \otimes |\psi_0 \rangle $, and otherwise - with $| 1 \rangle_i \otimes |\psi_1 \rangle $. Then normalize. Which option we chose determines result of the measurement.<jupyter_code>import numpy as np

class QubitRegister:
    """Resgiter of qubits of fixed size."""
    
    def __init__(self, n):
        self.n = n
        # Initalize 'all zeros' state.
        self.state = np.zeros(2**self.n, dtype=np.complex128)
        self.state[0] = 1
    
    def gate(self, A, qubit_ids):
        """Apply given gate to given set of qubits."""
        k = len(qubit_ids)
        assert A.shape == (2**k, 2**k), "Matrix has wrong shape."
        assert np.allclose(np.eye(2**k), A @ A.conj().T), "Matrix is not unitary."
        for id in qubit_ids:
            assert id >= 0  and id < self.n
        # Construct full matrix, acting on entire state space.
        U = np.zeros((2**self.n, 2**self.n), dtype = np.complex128)
        mask = sum(2**i for i in range(self.n) if i not in qubit_ids)
        def extract_idx(idx):
            ans = 0
            for i in range(k):
                ans += ((idx>>qubit_ids[i])%2)<<i
            return ans
        for i in range(2**self.n):
            for j in range(2**self.n):
                if (i & mask) == (j & mask):
                    U[i][j] = A[extract_idx(i)][extract_idx(j)]             
        # Apply the matrix.
        self.state = U @ self.state        
    
    def measure(self, qubit_id):
        """Measures given qubit in computational basis."""
        # Split into parts corresponding to basis states.
        p = np.zeros((2, 2**self.n), dtype = np.complex128)
        for i in range(2**self.n):
            p[(i >> qubit_id) % 2][i] = self.state[i]
        assert np.allclose(p[0]+p[1], self.state)
        # Compute probability of outcomes.
        prob0 = np.dot(p[0].conj(), p[0])
        prob1 = np.dot(p[1].conj(), p[1])
        assert np.allclose(prob0 + prob1, 1)
        # Choose the outcome.
        result = 0 if (np.random.random() < prob0) else 1
        # Collapse wavefunction and normalize it.
        self.state = p[result]
        self.state /= np.linalg.norm(self.state)
        return result
    
    # Here implementation of the simulator ends.
    # Below are just definitions of common gates, for convenience.
    def X(self, qubit_id):
        self.gate(np.array([[0, 1], [1, 0]]), [qubit_id])
        
    def Y(self, qubit_id):
        self.gate(np.array([[0, -1j], [1j, 0]]), [qubit_id])

    def Z(self, qubit_id):
        self.gate(np.array([[1, 0], [0, -1]]), [qubit_id])
        
    def Ry(self, theta, qubit_id):
        U = np.array([[np.cos(theta/2), -np.sin(theta/2)], [np.sin(theta/2), np.cos(theta/2)]])
        self.gate(U, [qubit_id])
        
    def H(self, qubit_id):
        """Hadamard gate."""
        self.gate(np.array([[1, 1], [1, -1]] / np.sqrt(2)), [qubit_id])
        
    def CNOT(self, control, target):
        U = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])
        self.gate(U, [control, target])
        
    def QFT(self, qubit_ids):
        """Quantum Fourier transform."""
        k = len(qubit_ids)
        w = np.exp(2j*np.pi/(2**k))
        U = np.array([[w**(i*j) for i in range(2**k)] for j in range(2**k)]) / np.sqrt(2**k)
        self.gate(U, qubit_ids)<jupyter_output><empty_output><jupyter_text>That's it! Strictly speaking, gates definitions are not necessary part of the simulator. So implementation of simulator itself takes about 50 lines of code.Now let's verify it works correctly.  Demo: preparing Bell state.<jupyter_code>reg = QubitRegister(2)
reg.H(0)
reg.CNOT(0, 1)
print(reg.state)<jupyter_output>[0.70710678+0.j 0.        +0.j 0.        +0.j 0.70710678+0.j]<jupyter_text>Demo: quantum teleportationThis is the circuit for quantum teleportation:![Teleportation](pic/teleportation.jpg)Let's teleport state $| \psi \rangle = \sin(\frac{\theta}{2}) | 0 \rangle + \sin(\frac{\theta}{2}) | 1 \rangle$.Let's fix value of $\theta$ and then many times prepare qubit 0 in state $| \psi \rangle$, apply teleportation circuit and measure qubit 2. Qubit 2 must be in state $| \psi \rangle$ before this measurement, so probability of measurement outcome 1 is $\sin^2(\theta/2)$.<jupyter_code>def teleport_circuit(reg):
    # Prepare shared pair of qubits (1 and 2).
    reg.H(2)
    reg.CNOT(2, 1)
    
    # Alice entangles her qubit from entangled pair (1) with qubit to teleport.
    reg.CNOT(0, 1)
    reg.H(0)
    
    # Alice makes 2 measurements
    m0 = reg.measure(0)
    m1 = reg.measure(1)

    # Bob uses results of measurement to rotate hois qubit.
    if m1==1: reg.X(2)
    if m0==1: reg.Z(2)
    
    
def test_teleport(theta):
    times = 1000 # Increase to improve accuracy. 
    measure1_count = 0
    for _ in range(times):
        reg = QubitRegister(3)
        reg.Ry(theta, 0)    
        teleport_circuit(reg)
        measure1_count += reg.measure(2)
    
    estimated_prob = measure1_count/times
    expected_prob = np.sin(0.5*theta)**2
    print('Estimated probability: %.06f' % estimated_prob)
    print('Expected probability: %.06f' % expected_prob)
    print('Error: %.06f' % np.abs(expected_prob - estimated_prob))
    
test_teleport(1.234)<jupyter_output>Estimated probability: 0.328000
Expected probability: 0.334767
Error: 0.006767<jupyter_text>Demo: Deutsch-Jozsa algorithmSee https://en.wikipedia.org/wiki/Deutsch%E2%80%93Jozsa_algorithm.We will consider function $f$ on N=5 qubits which is either balanced or constant. It will be implemented by an oracle acting on another qubit.Our task is to determine whether function is constant or balanced.<jupyter_code>N = 5

def deutsch_jozsa(oracle):
    """`oracle` is oracle for function f which is either constant or balanced."""
    # Prepare state.
    reg = QubitRegister(N+1)
    reg.X(N)
    for i in range(N+1):
        reg.H(i)

    # Act with oracle.
    oracle(reg)

    # Apply Hadamard gates.
    for i in range(N):
        reg.H(i)

    # Measurement.
    result = sum([2**i * reg.measure(i) for i in range(N)])
    if result == 0:
        return "Constant"
    else:
        return "Balanced"
    
# Generates oracle for random balanced function.
def random_balanced_function():
    U = np.zeros((2**(N+1), 2**(N+1)), dtype=np.complex128)
    Psi_1 = set(np.random.choice(range(2**N), 2**(N-1), replace=False))
    for psi in range(2**N):
        f = 1 if psi in Psi_1 else 0 # value of function on this state.
        for i in range(2):
            U[psi + (2**N) * i][psi + (2**N) * ((i + f) % 2)] = 1.0 
    
    def oracle(reg): 
        reg.gate(U, list(range(N+1)))
    return oracle
    
# Generates oracle for random constant function.
# There are only 2 such functions: no-op, and X on target qubit.
def random_constant_function():
    if np.random.random() < 0.5:
        def oracle(reg): 
            pass
        return oracle
    else:
        def oracle(reg): 
            reg.X(N)
        return oracle
    
for i in range(10):
    oracle = random_constant_function()
    result = deutsch_jozsa(oracle)
    assert(result == 'Constant')
print('Constant OK.')
    
for i in range(100):
    oracle = random_balanced_function()
    result = deutsch_jozsa(oracle)
    assert(result == 'Balanced')    
print('Balanced OK.')
a = np.array([0, 1j])
print(np.dot(a.conj(), a))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12884333821376281, 'fraction_numerical': 0.028864254340096214, 'mean_word_length': 2.840562248995984, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '88878', 'license': 'MIT', 'n_tokens_mistral': 3311, 'n_tokens_neox': 3079, 'n_words': 1042, 'path': 'SimpleSimulator.ipynb', 'repo_name': 'fedimser/quant_comp', 'size': '9572'}