<jupyter_start><jupyter_text>Author: Jiaying (<NAME>This implementation is based on pyQuil version 0.2.4. It is exactly the same experiment with the QISKit implementation, thus there is no detailed explanation here.The input in this implementation is  $A=\begin{pmatrix}2 & 1\\ 1 & 2\end{pmatrix}$ and $\vec{b}=\begin{pmatrix}0\\ 1\end{pmatrix}$<jupyter_code>from pyquil import Program, get_qc
from pyquil.gates import *
from math import pi
import numpy as np<jupyter_output><empty_output><jupyter_text>I find that there is a quantum gate RY in pyQuil, but it is only for 1 qubit. What I need is the controlled-RY, which can input two qubits as the control qubit and target qubit. Thus we can use the below cell to define the quantum gate by providing the unitary matrix (reference: http://docs.rigetti.com/en/stable/basics.htmldefining-new-gates).Another way to solve this problem is by using the `CONTROLLED` keyword, which the SDK tools will happily consume, but pyQuil is still learning to recognize: https://github.com/rigetti/pyquil/pull/808.<jupyter_code># First we define the controlled RY gate from a matrix
from pyquil.parameters import Parameter, quil_exp, quil_sin, quil_cos
from pyquil.quilbase import DefGate

# Define the new gate controlled-RY from a matrix
theta = Parameter('theta')
cry = np.array([
    [1, 0, 0, 0],
    [0, 1, 0, 0],
    [0, 0, quil_cos(theta / 2), -1 * quil_sin(theta / 2)],
    [0, 0, quil_sin(theta / 2), quil_cos(theta / 2)]
])

gate_definition = DefGate('CRY', cry, [theta])
CRY = gate_definition.get_constructor()<jupyter_output><empty_output><jupyter_text>The quantum circuit is:<jupyter_code>p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)<jupyter_output><empty_output><jupyter_text>Wavefunction Simulation (The two states that we need to pay attention to is the `0001` and `1001`):<jupyter_code>from pyquil.api import WavefunctionSimulator
wfn = WavefunctionSimulator().wavefunction(p)
print(wfn)<jupyter_output>(-0.0191221955+0j)|0000> + (0.0961337685+0j)|0001> + (0.9760625312+0j)|1000> + (0.1941509088+0j)|1001><jupyter_text>Run the circuit on a quantum virtual machine:<jupyter_code>#### 1024 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)

# run the program on a QVM
n_trials = 1024
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
#print(result[0])
#print(result[1])
#print(result[2])
#print(result[3])
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25
#### 2048 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)

# run the program on a QVM
n_trials = 2048
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25
#### 4096 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1) 

# run the program on a QVM
n_trials = 4096
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25
#### 8192 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)

# run the program on a QVM
n_trials = 8192
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25
#### 10240 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)

# run the program on a QVM
n_trials = 10240
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25
#### 15360 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)

# run the program on a QVM
n_trials = 15360
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25
#### 20480 trials ####
p = Program() 
p += gate_definition

p += X(3) # Initilize input b

### Phase estimation ###
## Superposition
p += H(1)
p += H(2)
## Controlled-U
p += PHASE(pi/2, 2)
p += PHASE(pi, 1)
p += CNOT(2, 3)
## QFT
p += SWAP(1, 2)
p += H(2)
p += CPHASE(-pi/2, 1, 2)
p += H(1)

### Controlled rotation Ry ###
p += X(1)
p += CRY(pi/8)(1, 0)
p += CRY(pi/16)(2, 0)

### Uncomputation ###
p += X(1)
p += H(1)
p += CPHASE(pi/2, 1, 2)
p += H(2)
p += SWAP(1, 2)
p += CNOT(2, 3)
p += PHASE(-pi, 1)
p += PHASE(-pi/2, 2)
p += H(2)
p += H(1)

# run the program on a QVM
n_trials = 20480
qc = get_qc('9q-square-qvm')
result = qc.run_and_measure(p, trials=n_trials)
x1_square = 0
x2_square = 0
for i in range(n_trials):
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 0:
        x1_square += 1
    if result[0][i] == 1 and result[1][i] == 0 and result[2][i] == 0 and result[3][i] == 1:
        x2_square += 1
p = x1_square/x2_square
print(x1_square)
print(x2_square)
print(p)  # Ideally, p=0.5²=0.25<jupyter_output>210
774
0.2713178294573643
{'fraction_non_alphanumeric': 0.21175177375836915, 'fraction_numerical': 0.07324872589187569, 'mean_word_length': 3.7589158345221114, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 24, 'https://': 1, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '253334', 'license': 'MIT', 'n_tokens_mistral': 4883, 'n_tokens_neox': 4536, 'n_words': 1248, 'path': 'Pyquil_HHL_2by2.ipynb', 'repo_name': 'yangjy0826/pyquil', 'size': '10013'}