In [29]:
from qiskit.ml.datasets import *
from qiskit import QuantumCircuit
from qiskit.aqua.components.optimizers import COBYLA, ADAM, SPSA, SLSQP, POWELL, L_BFGS_B, TNC, AQGD
from qiskit.circuit.library import ZZFeatureMap, RealAmplitudes
from qiskit.quantum_info import Statevector

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.preprocessing import normalize
from sklearn.model_selection import train_test_split
from sklearn.utils import shuffle

import warnings
warnings.filterwarnings("ignore")

%matplotlib inline
In [30]:
# constants
n = 4
RANDOM_STATE = 42
LR = 1e-3
class_labels = ['yes', 'no']
In [31]:
def normalizeData(DATA_PATH = "../../Data/Processed/data.csv"):
    """
    Normalizes the data
    """
    # Reads the data
    data = pd.read_csv(DATA_PATH)
    data = shuffle(data, random_state=RANDOM_STATE)
    X, Y = data[['sex', 'cp', 'exang', 'oldpeak']].values, data['num'].values
    # normalize the data
    X = normalize(X)
    X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.3, random_state=RANDOM_STATE)
    return X_train, X_test, Y_train, Y_test
In [32]:
X_train, X_test, Y_train, Y_test = normalizeData()
In [33]:
sv = Statevector.from_label('0' * n)
feature_map = ZZFeatureMap(n, reps=1)
var_form = RealAmplitudes(n, reps=1)
circuit = feature_map.combine(var_form)
circuit.draw(output='mpl')
Out[33]:
In [34]:
def get_data_dict(params, x):
    parameters = {}
    for i, p in enumerate(feature_map.ordered_parameters):
        parameters[p] = x[i]
    for i, p in enumerate(var_form.ordered_parameters):
        parameters[p] = params[i]
    return parameters
In [35]:
X_train[0]
Out[35]:
array([0.4472136 , 0.89442719, 0.        , 0.        ])
In [36]:
data = X_train[0]
params = np.array([0.1, 1.2, 0.02, 0.1, 0.1, 1.2, 0.02, 0.1])
circ_ = circuit.assign_parameters(get_data_dict(params, data))
circ_.draw(plot_barriers=True)
Out[36]:
In [37]:
def assign_label(bit_string, class_labels):
    hamming_weight = sum([int(k) for k in list(bit_string)])
    is_odd_parity = hamming_weight & 1
    if is_odd_parity:
        return class_labels[1]
    else:
        return class_labels[0]
In [38]:
def return_probabilities(counts, class_labels):
    shots = sum(counts.values())
    result = {class_labels[0]: 0,
    class_labels[1]: 0}
    for key, item in counts.items():
        label = assign_label(key, class_labels)
        result[label] += counts[key]/shots
    return result
In [39]:
return_probabilities({'00' : 10, '01': 10, '11': 20}, class_labels)
Out[39]:
{'yes': 0.75, 'no': 0.25}
In [40]:
def classify(x_list, params, class_labels):
    qc_list = []
    for x in x_list:
        circ_ = circuit.assign_parameters(get_data_dict(params, x))
        qc = sv.evolve(circ_)
        qc_list += [qc]
        probs = []
    for qc in qc_list:
        counts = qc.to_counts()
        prob = return_probabilities(counts, class_labels)
        probs += [prob]
    return probs
In [41]:
x = np.asarray([X_train[0]])
classify(x, params=np.array([0.8, -0.5, 1.5, 0,5, 0.8, -0.5, 1.5, 0,5]), class_labels=class_labels)
Out[41]:
[{'yes': 0.6338620140853923, 'no': 0.3661379859146076}]
In [42]:
def cost_estimate_sigmoid(probs, expected_label): # probability of labels vs actual labels
    p = probs.get(expected_label)
    sig = None
    if np.isclose(p, 0.0):
        sig = 1
    elif np.isclose(p, 1.0):
        sig = 0
    else:
        denominator = np.sqrt(2*p*(1-p))
        x = np.sqrt(200)*(0.5-p)/denominator
        sig = 1/(1+np.exp(-x))
    return sig
In [43]:
def mse_cost(probs, expected_label):
    p = probs.get(expected_label)
    actual, pred = np.array(1), np.array(p)
    return np.square(np.subtract(actual,pred)).mean()
In [44]:
def rmse_cost(probs, expected_label): 
    p = probs.get(expected_label)    
    actual, pred = np.array(1), np.array(pred)
    return np.sqrt(np.square(np.subtract(actual,pred)).mean())
In [45]:
cost_list = []
def cost_function(X, Y, class_labels, params, shots=100, print_value=False):
    # map training input to list of labels and list of samples
    cost = 0
    training_labels = []
    training_samples = []
    for sample in X:
        training_samples += [sample]
    for label in Y:
        if label == 0:
            training_labels += [class_labels[0]]
        elif label == 1:
            training_labels += [class_labels[1]]
    probs = classify(training_samples, params, class_labels)
    # evaluate costs for all classified samples
    for i, prob in enumerate(probs):
        cost += mse_cost(prob, training_labels[i])
    cost /= len(training_samples)
    # print resulting objective function
    if print_value:
        print('%.4f' % cost)
    # return objective value
    cost_list.append(cost)
    return cost
In [46]:
cost_function(X_train, Y_train, class_labels, params)
Out[46]:
0.22354962821107346
In [47]:
def compare_optim(optims):
    for opt in optims:
        if opt == 'SPSA':
            optimizer = SPSA(max_trials = 10)
        elif opt == 'SLSQP':
            optimizer = SLSQP(maxiter = 10)
        elif opt == 'TNC':
            optimizer = TNC(maxiter = 10)
        elif opt == 'POWELL':
            optimizer = POWELL(maxiter = 10)
        elif opt == 'ADAM':
            optimizer = ADAM(maxiter = 10)
        elif opt == 'L_BFGS_B':
            optimizer = L_BFGS_B(maxfun = 1000 ,  maxiter = 10)
        elif opt == 'COBYLA':
            optimizer = COBYLA(maxiter = 10)
        elif opt == "AQGD":
            optimizer = AQGD(maxiter=10)
        cost_list = []
        # define objective function for training
        objective_function = lambda params: cost_function(X_train, Y_train, class_labels, params, print_value=False)
        # randomly initialize the parameters
        np.random.seed(RANDOM_STATE)
        init_params = 2*np.pi*np.random.rand(n*(1)*2)
        # train classifier
        opt_params, value, _ = optimizer.optimize(len(init_params), objective_function, initial_point=init_params)
        # print results
        print()
        print('opt_params:', opt_params)
        print('opt_value: ', value)
In [48]:
compare_optim(["COBYLA", "ADAM", "SPSA", "SLSQP", "POWELL", "L_BFGS_B", "TNC", "AQGD"])
opt_params: [2.35330497 6.97351416 4.59925358 4.76148219 0.98029403 1.98014248
 0.3649501  5.44234523]
opt_value:  0.20637154014580933

opt_params: [2.34335531 5.98351175 4.58925731 3.77148558 0.97016639 0.99014585
 0.35494922 5.45235819]
opt_value:  0.28483386150197443

opt_params: [ 2.04927472  6.15709295  4.97282334  3.95333517  1.30242588  0.97419756
 -0.39591884  5.77216108]
opt_value:  0.22048075626519356

opt_params: [ 2.67580429  6.73085048  3.67672534  3.86981027  0.77933055  1.82929249
 -0.770624    5.59825639]
opt_value:  0.17596546056746962

opt_params: [ 0.12309248  7.91780579  8.04906378  3.16355487 -0.49750456 -0.64397726
  0.18959737  4.29371839]
opt_value:  0.1650836282262177

opt_params: [ 2.90497708  6.41074425  3.89803076  3.09054985  0.22344208  1.9978181
 -0.64476327  5.36431209]
opt_value:  0.16971839038652609

opt_params: [2.17425836 6.41728759 4.09463744 4.72008948 0.90205469 1.14542438
 0.24674838 5.83797257]
opt_value:  0.21844982063207888
------------------------------------------------------------
IndexError                 Traceback (most recent call last)
<ipython-input-48-cbe20ebd7e7c> in <module>
----> 1 compare_optim(["COBYLA", "ADAM", "SPSA", "SLSQP", "POWELL", "L_BFGS_B", "TNC", "AQGD"])

<ipython-input-47-b07473820ca6> in compare_optim(optims)
     24         init_params = 2*np.pi*np.random.rand(n*(1)*2)
     25         # train classifier
---> 26         opt_params, value, _ = optimizer.optimize(len(init_params), objective_function, initial_point=init_params)
     27         # print results
     28         print()

~/Documents/VENVS/qosf/lib/python3.6/site-packages/qiskit/aqua/components/optimizers/aqgd.py in optimize(self, num_vars, objective_function, gradient_function, variable_bounds, initial_point)
    316                 # Calculate objective function and estimate of analytical gradient
    317                 objval, gradient = \
--> 318                     self._compute_objective_fn_and_gradient(params, objective_function)
    319 
    320                 logger.info(" Iter: %4d | Obj: %11.6f | Grad Norm: %f",

~/Documents/VENVS/qosf/lib/python3.6/site-packages/qiskit/aqua/components/optimizers/aqgd.py in _compute_objective_fn_and_gradient(self, params, obj)
    149 
    150         # return the objective function value
--> 151         obj_value = values[0]
    152 
    153         # return the gradient values

IndexError: too many indices for array
In [56]:
cost_list = []
optimizer = L_BFGS_B(maxiter=100)
# optimizer = ADAM(maxiter=100, lr=LR)

# define objective function for training
objective_function = lambda params: cost_function(X_train, Y_train, class_labels, params, print_value=True)
# randomly initialize the parameters
np.random.seed(RANDOM_STATE)
init_params = 2*np.pi*np.random.rand(n*(1)*2)
# train classifier
opt_params, value, _ = optimizer.optimize(len(init_params), objective_function, initial_point=init_params)
# print results
print()
print('opt_params:', opt_params)
print('opt_value: ', value)
0.2878
0.2878
0.2878
0.2878
0.2878
0.2878
0.2878
0.2878
0.2878
0.2047
0.2047
0.2047
0.2047
0.2047
0.2047
0.2047
0.2047
0.2047
0.1954
0.1954
0.1954
0.1954
0.1954
0.1954
0.1954
0.1954
0.1954
0.2012
0.2012
0.2012
0.2012
0.2012
0.2012
0.2012
0.2012
0.2012
0.1792
0.1792
0.1792
0.1792
0.1792
0.1792
0.1792
0.1792
0.1792
0.1744
0.1744
0.1744
0.1744
0.1744
0.1744
0.1744
0.1744
0.1744
0.1733
0.1733
0.1733
0.1733
0.1733
0.1733
0.1733
0.1733
0.1733
0.1718
0.1718
0.1718
0.1718
0.1718
0.1718
0.1718
0.1718
0.1718
0.1715
0.1715
0.1715
0.1715
0.1715
0.1715
0.1715
0.1715
0.1715
0.1710
0.1710
0.1710
0.1710
0.1710
0.1710
0.1710
0.1710
0.1710
0.1704
0.1704
0.1704
0.1704
0.1704
0.1704
0.1704
0.1704
0.1704
0.2180
0.2180
0.2180
0.2180
0.2180
0.2180
0.2180
0.2180
0.2180
0.1697
0.1697
0.1697
0.1697
0.1697
0.1697
0.1697
0.1697
0.1697
0.1700
0.1700
0.1700
0.1700
0.1700
0.1700
0.1700
0.1700
0.1700
0.1694
0.1694
0.1694
0.1694
0.1694
0.1694
0.1694
0.1694
0.1694
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1680
0.1680
0.1680
0.1680
0.1680
0.1680
0.1680
0.1680
0.1680
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1672
0.1672
0.1672
0.1672
0.1672
0.1672
0.1672
0.1672
0.1672
0.1671
0.1671
0.1671
0.1671
0.1671
0.1671
0.1671
0.1671
0.1671
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1673
0.1669
0.1669
0.1669
0.1669
0.1669
0.1669
0.1669
0.1669
0.1669
0.1665
0.1665
0.1665
0.1665
0.1665
0.1665
0.1665
0.1665
0.1665
0.1662
0.1662
0.1662
0.1662
0.1662
0.1662
0.1662
0.1662
0.1662
0.1661
0.1661
0.1661
0.1661
0.1661
0.1661
0.1661
0.1661
0.1661
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1660
0.1659
0.1659
0.1659
0.1659
0.1659
0.1659
0.1659
0.1659
0.1659
0.1658
0.1658
0.1658
0.1658
0.1658
0.1658
0.1658
0.1658
0.1658
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1690
0.1657
0.1657
0.1657
0.1657
0.1657
0.1657
0.1657
0.1657
0.1657
0.1656
0.1656
0.1656
0.1656
0.1656
0.1656
0.1656
0.1656
0.1656
0.1654
0.1654
0.1654
0.1654
0.1654
0.1654
0.1654
0.1654
0.1654
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651
0.1651

opt_params: [3.26333934 7.90762422 4.89434665 3.15547364 0.49794644 2.49425875
 0.18009411 5.13728541]
opt_value:  0.16508072929082526
In [59]:
fig = plt.figure()
plt.plot(range(0,378,1), cost_list)
plt.xlabel('Steps')
plt.ylabel('Cost value')
plt.title("L_BFGS_B Cost value against steps")
plt.show()
fig.savefig('../../Output/Figures/costvssteps.jpeg')
In [21]:
def test_model(X, Y, class_labels, params):
    accuracy = 0
    training_labels = []
    training_samples = []
    for sample in X:
        training_samples += [sample]
    probs = classify(training_samples, params, class_labels)
    for i, prob in enumerate(probs):
        if (prob.get('yes') >= prob.get('no')) and (Y_test[i] == 0):
            accuracy += 1
        elif (prob.get('no') >= prob.get('yes')) and (Y_test[i] == 1):
            accuracy += 1
    accuracy /= len(Y_test)
    print("Test accuracy: {}\n".format(accuracy))
In [22]:
test_model(X_test, Y_test, class_labels, opt_params)
Test accuracy: 0.7613636363636364

In [ ]: