In [13]:
ansatz = TwoLocal(op.num_qubits, 'ry', 'cx', reps=3, entanglement='circular')
ansatz.draw('mpl')
Out[13]:
import numpy as np
from IPython.display import display, HTML
from docplex.mp.model import Model
from qiskit import Aer
from qiskit.circuit.library import TwoLocal
from qiskit.algorithms import VQE
from qiskit.aqua.algorithms import NumPyMinimumEigensolver
from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SPSA, SLSQP
from qiskit.opflow import Z, X, I, Gradient
from qiskit.optimization.applications.ising import docplex
from qiskit.optimization.algorithms import MinimumEigenOptimizer
from qiskit.optimization.problems import QuadraticProgram
from qiskit.opflow.primitive_ops.pauli_sum_op import PauliSumOp
from qiskit.utils import QuantumInstance
from pandas import DataFrame
from origin import Origin
from doublets import Doublet
from triplets import Triplet
C:\Users\david\.conda\envs\Qiskit\lib\site-packages\qiskit\optimization\__init__.py:92: DeprecationWarning: The package qiskit.optimization is deprecated. It was moved/refactored to qiskit_optimization (pip install qiskit-optimization). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
warn_package('optimization', 'qiskit_optimization', 'qiskit-optimization')
data_set = np.load('simplified_experiment.npy', allow_pickle=True)
hits_per_plane = [len(k) for k in data_set]
for i, hits in enumerate(hits_per_plane):
print(f'Number of hits in plane {i}: {hits}')
# constants
LIMIT_SIGMA_X = 3
LIMIT_SIGMA_Y = 3
# testing with smaller samples
TEST_SAMPLE_SIZE = 4
sub_dict_0 = dict(list(data_set[0].items())[:TEST_SAMPLE_SIZE])
sub_dict_1 = dict(list(data_set[1].items())[:TEST_SAMPLE_SIZE])
sub_dict_2 = dict(list(data_set[2].items())[:TEST_SAMPLE_SIZE])
sub_dict_3 = dict(list(data_set[3].items())[:TEST_SAMPLE_SIZE])
origin_object = Origin(data_set[0])
print(f'Sigma x limit set to: {LIMIT_SIGMA_X}')
print(f'Sigma x limit set to: {LIMIT_SIGMA_Y}')
print(f'Test sample size: {TEST_SAMPLE_SIZE}')
Number of hits in plane 0: 2807 Number of hits in plane 1: 2801 Number of hits in plane 2: 2793 Number of hits in plane 3: 2783 Sigma x limit set to: 3 Sigma x limit set to: 3 Test sample size: 4
# Doublets
doublet_list_01 = []
doublet_list_12 = []
doublet_list_23 = []
# fill plane 0 -> 1 doublets
for p0_key, p0_value in zip(sub_dict_0.keys(), sub_dict_0.values()):
for p1_key, p1_value in zip(sub_dict_1.keys(), sub_dict_1.values()):
doublet = Doublet(p0_key, p0_value, p1_key, p1_value)
if doublet.is_valid_doublet(origin_object, LIMIT_SIGMA_X, LIMIT_SIGMA_Y):
doublet_list_01.append(doublet)
# fill plane 1 -> 2 doublets
for p1_key, p1_value in zip(sub_dict_1.keys(), sub_dict_1.values()):
for p2_key, p2_value in zip(sub_dict_2.keys(), sub_dict_2.values()):
doublet = Doublet(p1_key, p1_value, p2_key, p2_value)
if doublet.is_valid_doublet(origin_object, LIMIT_SIGMA_X, LIMIT_SIGMA_Y):
doublet_list_12.append(doublet)
# fill plane 2 -> 3 doublets
for p2_key, p2_value in zip(sub_dict_2.keys(), sub_dict_2.values()):
for p3_key, p3_value in zip(sub_dict_3.keys(), sub_dict_3.values()):
doublet = Doublet(p2_key, p2_value, p3_key, p3_value)
if doublet.is_valid_doublet(origin_object, LIMIT_SIGMA_X, LIMIT_SIGMA_Y):
doublet_list_23.append(doublet)
print(f'Number of doublets created from plane 0 -> 1: {len(doublet_list_01)}')
print(f'Number of doublets created from plane 1 -> 2: {len(doublet_list_12)}')
print(f'Number of doublets created from plane 2 -> 3: {len(doublet_list_23)}')
count_01 = 0
count_12 = 0
count_23 = 0
for d in doublet_list_01:
if d.is_correct_match():
count_01 += 1
for d in doublet_list_12:
if d.is_correct_match():
count_12 += 1
for d in doublet_list_23:
if d.is_correct_match():
count_23 += 1
print('---------------------------------')
print(f'Fraction of correct doublets recognized for set limits: {np.round(count_01 / len(sub_dict_1), 3)}')
print(f'Fraction of correct doublets recognized for set limits: {np.round(count_12 / len(sub_dict_2), 3)}')
print(f'Fraction of correct doublets recognized for set limits: {np.round(count_23 / len(sub_dict_3), 3)}')
Number of doublets created from plane 0 -> 1: 16 Number of doublets created from plane 1 -> 2: 16 Number of doublets created from plane 2 -> 3: 13 --------------------------------- Fraction of correct doublets recognized for set limits: 1.0 Fraction of correct doublets recognized for set limits: 1.0 Fraction of correct doublets recognized for set limits: 1.0
# Triplets
triplet_list_012 = []
triplet_list_123 = []
for d1 in doublet_list_01:
for d2 in doublet_list_12:
if d1.get_p2_key() == d2.get_p1_key():
triplet = Triplet(d1, d2)
if triplet.is_valid_triplet(origin_object, 1/20 * LIMIT_SIGMA_X, 1/20 * LIMIT_SIGMA_Y):
triplet_list_012.append(triplet)
for d2 in doublet_list_12:
for d3 in doublet_list_23:
if d2.get_p2_key() == d3.get_p1_key():
triplet = Triplet(d2, d3)
if triplet.is_valid_triplet(origin_object, 1/20 * LIMIT_SIGMA_X, 1/20 * LIMIT_SIGMA_Y):
triplet_list_123.append(triplet)
print(f'Number of triplets created from plane 0 -> 1 -> 2: {len(triplet_list_012)}')
print(f'Number of triplets created from plane 1 -> 2 -> 3: {len(triplet_list_123)}')
count_012 = 0
count_123 = 0
correct_state = ''
for i, t in enumerate(triplet_list_012):
if t.is_correct_match():
count_012 += 1
correct_state += '1'
else:
correct_state += '0'
for j, t in enumerate(triplet_list_123):
if t.is_correct_match():
count_123 += 1
correct_state += '1'
else:
correct_state += '0'
print('Correct triplet list state: ', correct_state)
print(f'Fraction of correct triplets recognized for set limits: {np.round(count_012 / len(sub_dict_2), 3)}')
print(f'Fraction of correct triplets recognized for set limits: {np.round(count_123 / len(sub_dict_3), 3)}')
Number of triplets created from plane 0 -> 1 -> 2: 6 Number of triplets created from plane 1 -> 2 -> 3: 4 Correct triplet list state: 1010111111 Fraction of correct triplets recognized for set limits: 1.0 Fraction of correct triplets recognized for set limits: 1.0
The q's are in this case triplets of consecutive hits through the detector (e.g hit on plane 1, 2, 3 or 2, 3, 4) So the objective transforms to:
\begin{equation} \tag{QUBO: triplets} H(T) = \sum_i a_i T_i + \sum_{i<j} b_{ij}T_iT_j \end{equation}How to calculate: Let $P_i$ = {1, 2, 3} and $P_j$ = {2, 3, 4} be the sets of planes where the $T_i$ and $T_j$ come from.
a_i = []
b_ij = []
for t_i in triplet_list_012:
a_i.append(max(t_i.get_sigma_diff_at_origin(origin_object)))
b_fixed_i_variable_j = []
for t_j in triplet_list_123:
b_fixed_i_variable_j.append(t_i.conflict_with_other_triplet(origin_object, t_j))
b_ij.append(b_fixed_i_variable_j)
# build model with docplex
mdl = Model('QUBO')
x_i = [mdl.binary_var() for i in range(len(triplet_list_012))] # creating binary variables
x_j = [mdl.binary_var() for j in range(len(triplet_list_123))] # creating binary variables
# Create the Hamiltonian,
objective = mdl.sum(a_i[i] * x_i[i] for i in range(len(x_i))) + \
mdl.sum(b_ij[i][j] * x_i[i] * x_j[j] for i in range(len(x_i)) for j in range(len(x_j)))
mdl.minimize(objective)
print(mdl.export_as_lp_string())
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: QUBO
Minimize
obj: 0.000784435655 x1 + 0.114230452797 x2 + 0.006672050882 x3
+ 0.114758190618 x4 + 0.004781701073 x5 + 0.017807995856 x6 + [
- 1.939587965365 x1*x7 + 0.628073662404 x2*x8 + 0.995824239618 x2*x9
- 1.986696526418 x3*x8 + 0.956346572371 x4*x7 + 0.614610212867 x4*x8
- 1.990365011654 x5*x9 - 1.924500008322 x6*x10 ]/2
Subject To
Bounds
0 <= x1 <= 1
0 <= x2 <= 1
0 <= x3 <= 1
0 <= x4 <= 1
0 <= x5 <= 1
0 <= x6 <= 1
0 <= x7 <= 1
0 <= x8 <= 1
0 <= x9 <= 1
0 <= x10 <= 1
Binaries
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
End
qubo = QuadraticProgram()
qubo.from_docplex(mdl)
op, offset = qubo.to_ising()
print(f'offset: {offset}')
print(f'operator: \n {op}')
offset: -0.45126943962136257 operator: SummedOp([ 0.24205627784319803 * IIIIIIIIIZ, -0.260102464151179 * IIIIIIIIZI, 0.24500104036116124 * IIIIIIIZII, -0.25374869346397266 * IIIIIIZIII, 0.2464047759200888 * IIIIIZIIII, 0.2316585031120662 * IIIIZIIIII, -0.24244849567057897 * IIIZIIIIIZ, 0.12290517412415872 * IIIZIIIIII, 0.11954332154642025 * IIIZIIZIII, 0.07850920780053372 * IIZIIIIIZI, 0.0930015813932435 * IIZIIIIIII, -0.24833706580219372 * IIZIIIIZII, 0.07682627660841651 * IIZIIIZIII, 0.12447802995229879 * IZIIIIIIZI, 0.124317596504399 * IZIIIIIIII, -0.24879562645669778 * IZIIIZIIII, -0.2405625010402131 * ZIIIZIIIII, 0.2405625010402131 * ZIIIIIIIII ])
C:\Users\david\.conda\envs\Qiskit\lib\site-packages\qiskit\aqua\operators\operator_base.py:46: DeprecationWarning: The package qiskit.aqua.operators is deprecated. It was moved/refactored to qiskit.opflow (pip install qiskit-terra). For more information see <https://github.com/Qiskit/qiskit-aqua/blob/master/README.md#migration-guide>
warn_package('aqua.operators', 'qiskit.opflow', 'qiskit-terra')
# VQE module can't use the docplex to ising mdoels, so here is a quick workaround:
operator_list = []
for operator in op:
aux_list = str(operator).split(' * ')
t = tuple([aux_list[1], aux_list[0]])
operator_list.append(t)
H = PauliSumOp.from_list(operator_list)
# Exact Result
ee = NumPyMinimumEigensolver(H)
result = ee.run()
print('Ising: ', np.real(result['eigenvalue']))
Ising: -3.4392591327910327
seed = 10598 # seeding to take the random out of the calculation ... at least a bit
backend = Aer.get_backend('qasm_simulator')
quantum_instance = QuantumInstance(backend, seed_simulator=seed, seed_transpiler=seed, shots=8000)
# Choose an optimizer
optimizer_type = 'SPSA'
if optimizer_type == 'COBYLA':
optimizer = COBYLA(maxiter=3000)
elif optimizer_type == 'L_BFGS_B':
optimizer = L_BFGS_B(maxfun=5000)
elif optimizer_type == 'SPSA':
optimizer = SPSA(maxiter=1000)
elif optimizer_type == 'SLSQP':
optimizer = SLSQP(maxiter=4000)
ansatz = TwoLocal(op.num_qubits, 'ry', 'cx', reps=3, entanglement='circular')
ansatz.draw('mpl')
intermediate_info = {
'nfev': [],
'parameters': [],
'energy': [],
'stddev': []
}
def callback(nfev, parameters, energy, stddev):
intermediate_info['nfev'].append(nfev)
intermediate_info['parameters'].append(parameters)
intermediate_info['energy'].append(energy)
intermediate_info['stddev'].append(stddev)
# construct VQE
vqe = VQE(ansatz=ansatz,
optimizer=optimizer,
quantum_instance=quantum_instance,
callback=callback)
result_vqe = vqe.compute_minimum_eigenvalue(H)
# print results
print('energy: ', result_vqe.eigenvalue.real)
energy: -3.379160088488672
print(intermediate_info['nfev'][-1])
2051
print(result_vqe.eigenstate)
{'1111110101': 0.9794130895592523, '0100101101': 0.011180339887498949, '1000100011': 0.011180339887498949, '1111111101': 0.0873212459828649, '0111110101': 0.06892024376045111, '0000100010': 0.011180339887498949, '1111010101': 0.02958039891549808, '1000101001': 0.015811388300841896, '1011110101': 0.0806225774829855, '1111110110': 0.06614378277661477, '1101110101': 0.057008771254956896, '1110010101': 0.019364916731037084, '0000101010': 0.027386127875258306, '0010101000': 0.015811388300841896, '0000000101': 0.011180339887498949, '1000111011': 0.011180339887498949, '0100101010': 0.011180339887498949, '1111110001': 0.035355339059327376, '1100110101': 0.027386127875258306, '1111000101': 0.02958039891549808, '0001101010': 0.015811388300841896, '0000101101': 0.011180339887498949, '1111110011': 0.03872983346207417, '0010101001': 0.011180339887498949, '0011100101': 0.011180339887498949, '1110110101': 0.05, '0100101001': 0.015811388300841896, '0011110101': 0.019364916731037084, '0111111101': 0.015811388300841896, '0100010101': 0.011180339887498949, '0101110101': 0.011180339887498949, '1000101111': 0.011180339887498949, '0000010101': 0.027386127875258306, '0001110101': 0.011180339887498949, '1111100101': 0.019364916731037084, '0000101001': 0.015811388300841896, '1011111101': 0.015811388300841896, '1111110111': 0.011180339887498949, '1010101011': 0.011180339887498949, '1111110100': 0.015811388300841896}
The state with the highest porobability has be read backwars --> 1010111111, which is the correct matching!