<jupyter_start><jupyter_text>**Qiskit Global Hackathon**  Quantum Anomaly Detection with Quantum Density Estimation  Data LoadWe will work with the following datasets: *   SpamBase: dataset from UCI Machine Learning Repository, includes 3485 emails classified as spam or non-spam. This dataset has 58 attributes. We treat spam as outliers. The anomaly ratio is 0.2*   Cardiotocography: data is also from UCI Machine Learning Repository and is related to heart diseases. This dataset contains 22 attributes. It describes 3classes: normal, suspect, and pathological. Normal patients are treated as inliers and the remaining as outliers. The anomaly ratio is 0.2  Imports and Libraries<jupyter_code>import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
!pip install qiskit
!pip install pylatexenc
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from qiskit import Aer, execute<jupyter_output><empty_output><jupyter_text>Load from Drive<jupyter_code>#Loading SpamBase dataset file
!gdown --id 1Hdq-Svp31YkTj7wFZY3D2fsnBM1tbr_t
#Loading Cardiotocography dataset file
!gdown --id 1sVygYNPrTQI9d-8kVnx7w5nGWOSosveX
spambase = np.load("SpamBase.npy")

cardio = np.load("Cardiotocography.npy")<jupyter_output><empty_output><jupyter_text>PreprocessingThis function extracts and separates the features and the labels from the whole dataset. We will apply it to both datasets. Normal data will be labeled as '1' and anomalies will be labeled as '0'.<jupyter_code>from sklearn.preprocessing import MinMaxScaler

def preprocessing(data):
    features = data[:,:-1]
    labels = data[:,-1]

    scaler = MinMaxScaler()
    scaler.fit(features)
    features = scaler.transform(features)

    labels = labels[:,np.newaxis]
    labels[labels == -1] = 0

    return features, labels
spam_X, spam_y = preprocessing(spambase)
cardio_X, cardio_y = preprocessing(cardio)
spam_X.shape
cardio_X.shape<jupyter_output><empty_output><jupyter_text>Random Fourier FeaturesTo perform dimensionality reduction, we will use Random Fourier Features, an approach stablished by <NAME>, 2007. RFF is a type of quantum feature map from data to quantum states, in such a way that it can fit more precisely the underlying probability distribution of data.<jupyter_code>from sklearn.kernel_approximation import RBFSampler
"""
  Code from https://arxiv.org/abs/2004.01227
"""
class QFeatureMap: 

  def get_dim(self, num_features):
    pass

  def batch2wf(self, X):
    pass

  def batch2dm(self, X):
    pass

  def batch2dm(self, X):
    psi = self.batch2wf(X)
    rho = np.einsum('...i,...j', psi, np.conj(psi))
    return rho

class QFeatureMap_rff(QFeatureMap):

  def __init__(self, rbf_sampler):
    self.rbf_sampler = rbf_sampler
    self.weights = np.array(rbf_sampler.random_weights_)
    self.offset = np.array(rbf_sampler.random_offset_)
    self.dim = rbf_sampler.get_params()['n_components']

  def get_dim(self, num_features):
    return self.dim

  def batch2wf(self, X):
    vals = np.dot(X, self.weights) + self.offset
    vals = np.cos(vals)
    vals *= np.sqrt(2.) / np.sqrt(self.dim)
    norms = np.linalg.norm(vals, axis=1)
    psi = vals / norms[:, np.newaxis]
    return psi<jupyter_output><empty_output><jupyter_text>Cardio 4x4 Mixed StateFirst, we will work with Cardiotocography dataset.The number of RFF we will use depends on the size of the density matrix we will build. In this case, we will use a 4x4 density matrix, so we will reduce our samples to arrays of length 4. Transforming each sample will allow us to convert those values into two-qubit states.<jupyter_code># Create the RandomFourierFeature map
dim = 4
rm_state = 2
sigma = 1.5  #dataset dependant
gamma = 1 / (2 * sigma**2)

feature_map_fourier = RBFSampler(gamma=gamma, n_components=dim, random_state=rm_state)
X_feat_train = feature_map_fourier.fit(cardio_X)

rffmap = QFeatureMap_rff(rbf_sampler=feature_map_fourier)

Crff = rffmap.batch2wf(cardio_X)

print(Crff.shape)<jupyter_output>(2060, 4)<jupyter_text>Once our samples have been reduced, it's important to normalize them. Also we will split the dataset into three groups, a training partition with only normal data (70%), and a test partition with a group of normal data (10%) and a group of outlier data (20%).<jupyter_code>Cones = []
Ctrain = []
Ctest_zeros = []
Ctest_ones = []

for i in range(Crff.shape[0]):
    if (cardio_y[i] == 1):
        Cones.append(Crff[i] / np.linalg.norm(Crff[i]))
    elif (cardio_y[i] == 0):
        Ctest_zeros.append(Crff[i] / np.linalg.norm(Crff[i]))

Ctest_ones = Cones[:len(Ctest_zeros)//2]
Ctrain = Cones[len(Ctest_zeros)//2:]

print(len(Ctrain))
print(len(Ctest_zeros))
print(len(Ctest_ones))
print(Ctrain[0])
print(sum(Ctrain[0]**2))<jupyter_output>[ 0.6607845   0.60963091 -0.29528374 -0.32329786]
1.0000000000000002<jupyter_text>Now, we will define the training state using the Mixed State approach. We will use the density matrix corresponding to the training state, $\rho_1$. So we have$$\rho_1 = \frac{1}{m} \sum |\phi_1 \rangle \langle \phi_1 |$$being $|\phi_1 \rangle$ the states labeled with '1' (aka training normal data), and $m$ the amount of samples. To construct the density matrices, only train samples will be used.<jupyter_code>import copy
def round_lambda(lambda_param_2):
  lambda_param = copy.copy(lambda_param_2)
  for i in range(4):
    if lambda_param[i] < 0.00000001:
      lambda_param[i] = 0.00000001
  lambda_param[3] = 1 - lambda_param[0] - lambda_param[1] - lambda_param[2]
  return lambda_param
Z_train = np.outer(np.abs(Ctrain[0]), np.abs(Ctrain[0]))
for i in range(1, len(Ctrain)):
    Z_train += np.outer(np.abs(Ctrain[i]), np.abs(Ctrain[i]))

Z_train *= 1/len(Ctrain)
print(Z_train)
lambda_P1_temp, U_train = np.linalg.eigh(Z_train)
lambda_train = round_lambda(lambda_P1_temp)
print(lambda_train)
print(U_train)
lambda_train = lambda_train / np.linalg.norm(lambda_train)
print(lambda_train)
print((sum(lambda_train**2)))<jupyter_output>[0.02166764 0.0333373  0.07158875 0.99664145]
0.9999999999999999<jupyter_text>The quantum circuit we propose will be based in a Qiskit gate called `Isometry`, that allows to prepare arbitrary unitary rotations.We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|00 \rangle$ obtained. This value will be proportional to the density value estimated by the circuit.<jupyter_code>backend2 = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 4x4
Cpred_ones = []

for i in range(len(Ctest_ones)):

    qc = QuantumCircuit(4, 2)

    qc.initialize(Ctest_ones[i], [0, 1])
    qc.initialize(lambda_train, [2, 3])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.cnot(3, 1)
    qc.cnot(2, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Cpred_ones.append((counts['00'])/5000)

print(Cpred_ones)
# Mixed state case
Cpred_zeros = []

for i in range(len(Ctest_zeros)):

    qc = QuantumCircuit(4, 2)
    
    qc.initialize(Ctest_zeros[i], [0, 1])
    qc.initialize(lambda_train, [2, 3])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.cnot(3, 1)
    qc.cnot(2, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Cpred_zeros.append((counts['00'])/5000)

print(Cpred_zeros)<jupyter_output>[0.0684, 0.1196, 0.0228, 0.0174, 0.0092, 0.5514, 0.789, 0.0748, 0.0082, 0.0876, 1.1682, 0.0102, 1.1918, 0.0184, 0.04, 0.1338, 0.3542, 0.2672, 0.248, 0.28, 0.011, 0.0152, 0.0174, 0.5672, 0.2106, 0.0248, 0.009, 0.021, 0.0268, 0.041, 0.2644, 0.3378, 0.016, 0.7752, 1.2374, 0.5204, 0.0172, 0.0104, 0.0708, 0.037, 0.1128, 0.0136, 0.0252, 0.115, 0.0616, 0.6658, 0.0126, 0.2878, 0.0214, 0.0672, 0.0874, 0.0114, 0.0582, 0.3758, 0.009, 0.2992, 0.0364, 0.1518, 0.0722, 0.3478, 0.111, 0.0648, 0.5902, 0.0242, 1.0012, 0.0164, 0.0106, 0.2602, 0.0152, 0.0772, 1.1338, 0.593, 0.0168, 0.8484, 0.0268, 0.0232, 0.0114, 0.0176, 0.4136, 0.7822, 0.3392, 0.0552, 1.6286, 0.7188, 0.1488, 0.0286, 0.0144, 0.0146, 0.84, 0.3134, 0.014, 0.015, 0.0148, 0.064, 0.06, 0.117, 0.056, 0.0116, 0.1024, 1.4278, 0.0128, 0.2292, 0.026, 0.0262, 0.069, 0.356, 0.0104, 0.0924, 0.0272, 0.015, 0.0238, 0.0434, 0.474, 0.0202, 0.0856, 0.0202, 0.0542, 0.063, 0.3098, 0.0218, 0.0302, 0.028, 0.405, 0.0528, 0.0622, 0.1258, 0.0154, 0.0336, 0.0142, [...]<jupyter_text>Once we have obtained the predicted density values, we need to determine a threshold $t$ such that if a sample has a density value lower than $t$, it is classified as an outlier. For this, we will define a validation partition from our test groups, and we'll try to find the best possible value for $t$ using only validation samples (based on the best possible value of weighted F1-score). After that, we'll use this threshold to compare with all the remaining test samples.<jupyter_code>Cpred_val = np.concatenate((Cpred_ones[:100], Cpred_zeros[:100]))
print(Cpred_val.shape)

Cy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Cy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Cy_val, Cpred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Cpred_val > t).astype(int)
    metric = f1_score(Cy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.3758<jupyter_text>We have obtained that our best threshold is at $t = 0.3758$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Cpr_zeros = (np.array(Cpred_zeros) > T).astype(int)

Cpr_ones = (np.array(Cpred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(106), np.zeros(312))), 
        np.concatenate((Cpr_ones[100:], Cpr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.8431    0.8782    0.8603       312
         1.0     0.5914    0.5189    0.5528       106

    accuracy                         0.7871       418
   macro avg     0.7172    0.6985    0.7065       418
weighted avg     0.7793    0.7871    0.7823       418<jupyter_text>Cardio 4x4 Pure StateNow, we will define another way to represent the training density matrix: the Pure State approach. First, we will define the state  $|\psi_1 \rangle$. $$ |\psi_1 \rangle = \frac {\sum |\phi_1 \rangle} {\left \| \sum |\phi_1 \rangle \right \|}$$being $|\phi_1 \rangle$ the states labeled with '1'. To construct this state, only train samples will be used.Once obtained this state, it represents all the information from training states. So we build the density matrix as$$\rho_1 = |\psi_1 \rangle \langle ψ_1 |$$<jupyter_code># Given the eigenvalues this function return a Unitary gate which converts the |0> -> |lambda>
def create_U_train(x_train_param, seed=0):
  np.random.seed(seed)
  size_U = len(x_train_param)
  U_train = np.zeros((size_U, size_U))
  x_1 = x_train_param
  U_train[:, 0] = x_1
  for i in range(1, size_U):
    x_i = np.random.randn(size_U)
    for j in range(0, i):
      x_i -= x_i.dot(U_train[:, j]) * U_train[:, j]
    x_i = x_i/np.linalg.norm(x_i)
    U_train[:, i] = x_i
  
  return U_train
phi_train = np.sum(np.abs(Ctrain), axis=0)
phi_train = phi_train / np.linalg.norm(phi_train)

U_train = create_U_train(phi_train)

print(U_train)<jupyter_output>[[ 0.53308287  0.2309175   0.69538315 -0.42301541]
 [ 0.46304288 -0.88279704 -0.0144998   0.07778448]
 [ 0.26452084  0.21352962  0.27364652  0.89975075]
 [ 0.65684295  0.34892921 -0.66434149 -0.0738656 ]]<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|00 \rangle$ obtained.<jupyter_code>backend = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 4x4
Cpred_ones = []

for i in range(len(Ctest_ones)):

    qc = QuantumCircuit(2, 2)
    qc.initialize(Ctest_ones[i], [0, 1])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Cpred_ones.append((counts['00'])/5000)
    except:
        Cpred_ones.append(0)

print(Cpred_ones)
# Pure state case
Cpred_zeros = []

for i in range(len(Ctest_zeros)):

    qc = QuantumCircuit(2, 2)
    qc.initialize(Ctest_zeros[i], [0, 1])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Cpred_zeros.append((counts['00'])/5000)
    except:
        Cpred_zeros.append(0)

print(Cpred_zeros)<jupyter_output>[0.0598, 0.1026, 0.0114, 0.0046, 0.0002, 0.556, 0.8118, 0.0616, 0, 0.0732, 1.194, 0, 1.2268, 0.0088, 0.0272, 0.127, 0.3356, 0.2658, 0.2252, 0.2706, 0.0002, 0.0048, 0.0056, 0.5674, 0.1846, 0.008, 0, 0.0084, 0.013, 0.0212, 0.2624, 0.3156, 0.0086, 0.8002, 1.223, 0.5228, 0.006, 0, 0.0554, 0.0276, 0.1162, 0, 0.0146, 0.1088, 0.038, 0.6774, 0, 0.29, 0.011, 0.0558, 0.0818, 0, 0.045, 0.3692, 0.0004, 0.315, 0.0282, 0.1598, 0.0712, 0.3218, 0.086, 0.0512, 0.5896, 0.0152, 0.9942, 0, 0.0004, 0.261, 0.0052, 0.0766, 1.1412, 0.62, 0, 0.842, 0.0176, 0.0182, 0.0004, 0.0142, 0.3906, 0.7982, 0.3298, 0.0412, 1.6536, 0.727, 0.1456, 0.0172, 0, 0.0072, 0.8646, 0.319, 0.0042, 0.0018, 0.0024, 0.0574, 0.0498, 0.098, 0.041, 0.0002, 0.0984, 1.4236, 0.0002, 0.232, 0.0134, 0.0162, 0.0604, 0.353, 0, 0.082, 0.0126, 0.0056, 0.012, 0.0444, 0.452, 0.0042, 0.074, 0.0086, 0.047, 0.066, 0.311, 0.012, 0.0176, 0.0048, 0.4014, 0.034, 0.0452, 0.128, 0.0084, 0.0138, 0.0014, 0.1082, 0.8632, 0.0014, 0.2948, 0.0034, 0.0724, 0.044, 0[...]<jupyter_text>Once we have obtained the predicted density values, we need to determine a threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier. For this, we will define a validation partition from our test groups, and we'll try to find the best possible value for $t$ (based on weighted F1-score) using only validation samples. After that, we'll use this threshold to compare with all the remaining test samples.<jupyter_code>Cpred_val = np.concatenate((Cpred_ones[:100], Cpred_zeros[:100]))
print(Cpred_val.shape)

Cy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Cy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Cy_val, Cpred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Cpred_val > t).astype(int)
    metric = f1_score(Cy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.3692<jupyter_text>We have obtained that our best threshold is at $t = 0.3692$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Cpr_zeros = (np.array(Cpred_zeros) > T).astype(int)

Cpr_ones = (np.array(Cpred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(106), np.zeros(312))), 
        np.concatenate((Cpr_ones[100:], Cpr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.8431    0.8782    0.8603       312
         1.0     0.5914    0.5189    0.5528       106

    accuracy                         0.7871       418
   macro avg     0.7172    0.6985    0.7065       418
weighted avg     0.7793    0.7871    0.7823       418<jupyter_text>Cardio 8x8 Mixed StateNow, we will use a 8x8 density matrix, so we will call a new RFF sampler to reduce our samples to arrays of length 8. Transforming each sample will allow us to convert those values into three-qubit states.<jupyter_code># Create the RandomFourierFeature map
dim = 8
rm_state = 2
sigma = 1.5  #dataset dependant
gamma = 1 / (2 * sigma**2)

feature_map_fourier = RBFSampler(gamma=gamma, n_components=dim, random_state=rm_state)
X_feat_train = feature_map_fourier.fit(cardio_X)

rffmap = QFeatureMap_rff(rbf_sampler=feature_map_fourier)

Crff = rffmap.batch2wf(cardio_X)

print(Crff.shape)<jupyter_output>(2060, 8)<jupyter_text>Once our samples have been reduced, it's important to normalize them. Also we will split the dataset into three groups, a training partition with only normal data (70%), and a test partition with a group of normal data (10%) and a group of outlier data (20%).<jupyter_code>Cones = []
Ctrain = []
Ctest_zeros = []
Ctest_ones = []

for i in range(Crff.shape[0]):
    if (cardio_y[i] == 1):
        Cones.append(Crff[i] / np.linalg.norm(Crff[i]))
    elif (cardio_y[i] == 0):
        Ctest_zeros.append(Crff[i] / np.linalg.norm(Crff[i]))

Ctest_ones = Cones[:len(Ctest_zeros)//2]
Ctrain = Cones[len(Ctest_zeros)//2:]

print(len(Ctrain))
print(len(Ctest_zeros))
print(len(Ctest_ones))
print(Ctrain[0])
print(sum(Ctrain[0]**2))<jupyter_output>[ 0.45302627  0.48066672  0.39341633  0.02643167 -0.45585263  0.42733187
  0.13250014  0.01677011]
0.9999999999999999<jupyter_text>Then we build our traning matrix $ρ_1$:<jupyter_code>Z_train = np.outer(np.abs(Ctrain[0]), np.abs(Ctrain[0]))
for i in range(1, len(Ctrain)):
    Z_train += np.outer(np.abs(Ctrain[i]), np.abs(Ctrain[i]))

Z_train *= 1/len(Ctrain)
print(Z_train)
lambda_P1_temp, U_train = np.linalg.eigh(Z_train)
lambda_train = round_lambda(lambda_P1_temp)
print(lambda_train)
print(U_train)
lambda_train = lambda_train / np.linalg.norm(lambda_train)
print(lambda_train)
print((sum(lambda_train**2)))<jupyter_output>[0.00184545 0.00270234 0.00418252 0.72947213 0.00877911 0.01471787
 0.01731873 0.68355584]
1.0<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|000 \rangle$ obtained. This value will be proportional to the density estimated by the circuit.<jupyter_code>backend2 = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 8x8
Cpred_ones = []

for i in range(len(Ctest_ones)):

    qc = QuantumCircuit(6, 3)

    qc.initialize(Ctest_ones[i], [0, 1, 2])
    qc.initialize(lambda_train, [3, 4, 5])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.cnot(5, 2)
    qc.cnot(4, 1)
    qc.cnot(3, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)

    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Cpred_ones.append((counts['000'])/5000)

print(Cpred_ones)
# Mixed state case
Cpred_zeros = []

for i in range(len(Ctest_zeros)):

    qc = QuantumCircuit(6, 3)

    qc.initialize(Ctest_zeros[i], [0, 1, 2])
    qc.initialize(lambda_train, [3, 4, 5])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.cnot(5, 2)
    qc.cnot(4, 1)
    qc.cnot(3, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)
    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Cpred_zeros.append((counts['000'])/5000)

print(Cpred_zeros)<jupyter_output>[1.2566, 0.8392, 0.823, 0.868, 0.692, 0.2758, 0.4006, 0.4426, 0.7928, 0.2826, 0.417, 1.2288, 0.1878, 0.6434, 0.674, 0.4046, 0.406, 0.3468, 0.2662, 0.338, 0.751, 0.8, 0.7308, 0.4616, 0.413, 0.3488, 0.7376, 0.796, 0.3672, 0.3688, 0.56, 0.504, 0.4114, 0.1402, 0.2516, 0.1738, 0.8018, 0.4604, 0.3956, 0.387, 0.631, 0.6238, 0.6196, 0.5906, 0.3412, 0.5356, 0.6682, 0.5734, 0.866, 0.6056, 0.5622, 0.8014, 1.2236, 0.483, 0.7098, 0.1732, 0.593, 0.2206, 0.8026, 0.4662, 0.7884, 0.3162, 0.4838, 0.667, 0.3612, 0.304, 0.7814, 0.4824, 0.5766, 0.617, 0.4032, 0.4908, 0.859, 0.4014, 0.732, 0.9402, 0.5302, 0.6352, 0.2528, 0.4224, 0.5464, 0.7854, 0.3786, 0.4862, 0.4392, 0.7524, 0.6824, 0.7852, 0.4416, 0.4924, 0.9034, 0.8302, 0.8522, 0.779, 1.1684, 0.3416, 0.991, 0.728, 0.7612, 0.421, 0.796, 0.622, 0.6702, 0.3504, 0.5222, 0.4788, 0.6306, 0.541, 0.7512, 1.072, 0.4098, 0.6652, 0.2902, 0.8412, 0.716, 0.9194, 1.1188, 0.6078, 0.3336, 0.8592, 0.3268, 0.3844, 0.8128, 0.651, 0.3386, 0.3244, 0.8476, 0.1844, 0.9834, 0.4[...]<jupyter_text>Once we have obtained the predicted density values, we need to determine the best possible threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier.<jupyter_code>Cpred_val = np.concatenate((Cpred_ones[:100], Cpred_zeros[:100]))
print(Cpred_val.shape)

Cy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Cy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Cy_val, Cpred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Cpred_val > t).astype(int)
    metric = f1_score(Cy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.504<jupyter_text>We have obtained that our best threshold is at $t = 0.504$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Cpr_zeros = (np.array(Cpred_zeros) > T).astype(int)

Cpr_ones = (np.array(Cpred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(106), np.zeros(312))), 
        np.concatenate((Cpr_ones[100:], Cpr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.7627    0.4327    0.5521       312
         1.0     0.2656    0.6038    0.3689       106

    accuracy                         0.4761       418
   macro avg     0.5141    0.5182    0.4605       418
weighted avg     0.6366    0.4761    0.5057       418<jupyter_text>Cardio 8x8 Pure StateNow, we'll build a new training density matrix of size 8x8 using the Pure State approach.<jupyter_code>phi_train = np.sum(np.abs(Ctrain), axis=0)
phi_train = phi_train / np.linalg.norm(phi_train)

U_train = create_U_train(phi_train)

print(U_train)<jupyter_output>[[ 0.28684681  0.36163466 -0.64907043  0.3039476   0.29823376 -0.3628901
   0.15672877 -0.16752971]
 [ 0.43304774 -0.19150853  0.0425077   0.01946241 -0.66638594 -0.35659966
  -0.20048597 -0.40271614]
 [ 0.44413513 -0.01120385 -0.3346563   0.02332315 -0.35547746  0.59271958
   0.24672536  0.38926279]
 [ 0.30137368  0.50622778  0.63860468  0.26581483  0.02231948  0.06976633
   0.39028257 -0.1294136 ]
 [ 0.42656204  0.29150598  0.0298679  -0.72831994  0.27817334  0.09661608
  -0.32529554 -0.09588955]
 [ 0.39638647 -0.61313144  0.11383773  0.26809529  0.49785862  0.24052065
  -0.04173887 -0.27323611]
 [ 0.13336378  0.21047524  0.11205656  0.46037306  0.02482021  0.01655667
  -0.74680723  0.39345019]
 [ 0.29196429 -0.26635886  0.17516306 -0.14753714  0.1194335  -0.5634755
   0.23825873  0.6346844 ]]<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|000 \rangle$ obtained.<jupyter_code>backend = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 8x8
Cpred_ones = []

for i in range(len(Ctest_ones)):

    qc = QuantumCircuit(3, 3)
    qc.initialize(Ctest_ones[i], [0, 1, 2])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Cpred_ones.append((counts['000'])/5000)
    except:
        Cpred_ones.append(0)

print(Cpred_ones)
# Pure state case
Cpred_zeros = []

for i in range(len(Ctest_zeros)):

    qc = QuantumCircuit(3, 3)
    qc.initialize(Ctest_zeros[i], [0, 1, 2])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Cpred_zeros.append((counts['000'])/5000)
    except:
        Cpred_zeros.append(0)

print(Cpred_zeros)<jupyter_output>[1.9668, 1.5754, 1.1734, 1.5554, 1.0302, 0.4776, 0.8382, 0.0148, 1.6616, 0.5102, 0.222, 2.0816, 0.2726, 0.9778, 1.1416, 0.0964, 0.0302, 0.5978, 0.0792, 0.7052, 1.4324, 1.0488, 1.2854, 0.901, 0.776, 0.6954, 1.4032, 1.3468, 0.4798, 0.7234, 1.0234, 0.952, 0.0068, 0.0656, 0.4422, 0.1028, 1.4892, 0.51, 0.6576, 0.7458, 1.1822, 0.9592, 1.1216, 1.0438, 0.5476, 1.132, 1.0804, 1.0708, 1.3552, 1.0736, 1.0506, 1.4146, 1.4974, 0.9296, 1.3532, 0.172, 1.105, 0.3474, 1.526, 0.9158, 1.3762, 0.5742, 1.0202, 1.0788, 0.1654, 0.0868, 1.5812, 0.972, 1.0764, 1.1818, 0.2366, 1.0306, 1.6814, 0.7872, 1.3044, 1.5446, 1.1526, 1.2164, 0.4124, 0.8394, 1.0416, 1.0838, 0.8034, 0.8376, 0.53, 1.3444, 1.3006, 1.3008, 0.876, 0.9544, 1.9064, 1.6878, 1.691, 1.4648, 1.7734, 0.6882, 1.8432, 1.1332, 1.4366, 0.8836, 0.9822, 1.285, 1.085, 0.6948, 0.931, 0.3902, 1.2404, 1.0554, 1.1688, 2.22, 0.0062, 1.1246, 0.5546, 1.2656, 1.4906, 1.9544, 1.7356, 1.0424, 0.3922, 1.6586, 0.7158, 0.8556, 1.6272, 1.3982, 0.2464, 0.6006, 1.7632, 0.0[...]<jupyter_text>Once we have obtained the predicted density values, we need to determine the best possible threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier.<jupyter_code>Cpred_val = np.concatenate((Cpred_ones[:100], Cpred_zeros[:100]))
print(Cpred_val.shape)

Cy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Cy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Cy_val, Cpred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Cpred_val > t).astype(int)
    metric = f1_score(Cy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>1.128<jupyter_text>We have obtained that our best threshold is at $t = 1.128$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Cpr_zeros = (np.array(Cpred_zeros) > T).astype(int)

Cpr_ones = (np.array(Cpred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(106), np.zeros(312))), 
        np.concatenate((Cpr_ones[100:], Cpr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.6815    0.5897    0.6323       312
         1.0     0.1351    0.1887    0.1575       106

    accuracy                         0.4880       418
   macro avg     0.4083    0.3892    0.3949       418
weighted avg     0.5429    0.4880    0.5119       418<jupyter_text>SpamBase 4x4 Mixed StateNow, we will work with the SpamBase dataset.We will do the same work done previously, using both Pure and Mixed States and with both 4x4 and 8x8 density matrices.Now we transform our data samples into arrays of length 4 and splitting them into training and test samples.<jupyter_code># Create the RandomFourierFeature map
dim = 4
rm_state = 2
sigma = 1  #dataset dependant
gamma = 1 / (2 * sigma**2)

feature_map_fourier = RBFSampler(gamma=gamma, n_components=dim, random_state=rm_state)
X_feat_train = feature_map_fourier.fit(spam_X)

rffmap = QFeatureMap_rff(rbf_sampler=feature_map_fourier)

Srff = rffmap.batch2wf(spam_X)

print(Srff.shape)
Sones = []
Strain = []
Stest_zeros = []
Stest_ones = []

for i in range(Srff.shape[0]):
    if (spam_y[i] == 1):
        Sones.append(Srff[i] / np.linalg.norm(Srff[i]))
    elif (spam_y[i] == 0):
        Stest_zeros.append(Srff[i] / np.linalg.norm(Srff[i]))

Stest_ones = Sones[:len(Stest_zeros)//2]
Strain = Sones[len(Stest_zeros)//2:]

print(len(Strain))
print(len(Stest_zeros))
print(len(Stest_ones))
print(Strain[0])
print(sum(Strain[0]**2))<jupyter_output>[ 0.35145319 -0.48431271 -0.78982361  0.13453817]
1.0<jupyter_text>We build our training density matrix with mixed states:<jupyter_code>Z_train = np.outer(np.abs(Strain[0]), np.abs(Strain[0]))
for i in range(1, len(Strain)):
    Z_train += np.outer(np.abs(Strain[i]), np.abs(Strain[i]))

Z_train *= 1/len(Strain)
print(Z_train)
lambda_P1_temp, U_train = np.linalg.eigh(Z_train)
lambda_train = round_lambda(lambda_P1_temp)
print(lambda_train)
print(U_train)
lambda_train = lambda_train / np.linalg.norm(lambda_train)
print(lambda_train)
print((sum(lambda_train**2)))<jupyter_output>[0.03324847 0.06247062 0.07618934 0.99457888]
1.0<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|00 \rangle$ obtained. This value will be proportional to the density estimated by the circuit.<jupyter_code>backend2 = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 4x4
Spred_ones = []

for i in range(len(Stest_ones)):

    qc = QuantumCircuit(4, 2)

    qc.initialize(Stest_ones[i], [0, 1])
    qc.initialize(lambda_train, [2, 3])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.cnot(3, 1)
    qc.cnot(2, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Spred_ones.append((counts['00'])/5000)

print(Spred_ones)
# Mixed state case
Spred_zeros = []

for i in range(len(Stest_zeros)):

    qc = QuantumCircuit(4, 2)
    
    qc.initialize(Stest_zeros[i], [0, 1])
    qc.initialize(lambda_train, [2, 3])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.cnot(3, 1)
    qc.cnot(2, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Spred_zeros.append((counts['00'])/5000)

print(Spred_zeros)<jupyter_output>[0.2156, 0.3606, 0.0436, 0.1496, 0.2066, 0.2036, 0.0158, 3.1044, 0.1336, 0.3454, 0.3726, 0.3356, 0.2196, 0.4036, 0.0176, 0.7834, 0.0586, 0.5296, 0.44, 0.0194, 0.4462, 0.0468, 0.0192, 0.1288, 0.6674, 0.0754, 0.7868, 0.5258, 0.0568, 0.064, 0.8208, 0.2622, 0.0216, 0.3928, 0.4178, 0.291, 0.3442, 0.6596, 0.0646, 0.0186, 1.8216, 0.7712, 0.0188, 0.0132, 1.0964, 0.392, 1.0282, 1.3628, 0.0494, 0.0918, 0.0474, 0.0542, 0.1338, 0.7372, 0.4646, 0.4912, 0.0204, 1.3506, 0.1338, 1.1238, 0.0084, 0.3286, 0.0878, 0.0334, 0.5476, 0.3642, 0.9862, 0.2754, 0.1064, 0.083, 0.2706, 0.3434, 0.932, 0.4802, 0.8014, 0.3824, 0.5466, 0.6966, 0.9994, 1.233, 0.885, 0.932, 0.0476, 0.5266, 0.1754, 0.07, 0.1778, 0.0358, 0.1724, 0.42, 0.034, 0.0142, 0.3766, 0.0072, 0.1166, 0.026, 0.5254, 0.0242, 0.9968, 1.101, 0.077, 0.1932, 0.2958, 0.3598, 0.3786, 0.0328, 0.068, 0.4258, 1.2466, 0.5896, 0.4102, 0.4264, 0.0924, 0.0162, 0.6212, 0.047, 0.0212, 0.0316, 0.0086, 0.1494, 0.5156, 0.3574, 0.3154, 0.0262, 0.102, 0.3564, 0.3788, 0.02[...]<jupyter_text>We need to determine a threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier.<jupyter_code>Spred_val = np.concatenate((Spred_ones[:100], Spred_zeros[:100]))
print(Spred_val.shape)

Sy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Sy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Sy_val, Spred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Spred_val > t).astype(int)
    metric = f1_score(Sy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.542<jupyter_text>We have obtained that our best threshold is at $t = 0.542$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Spr_zeros = (np.array(Spred_zeros) > T).astype(int)

Spr_ones = (np.array(Spred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(248), np.zeros(597))), 
        np.concatenate((Spr_ones[100:], Spr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.8476    0.7454    0.7932       597
         1.0     0.5250    0.6774    0.5915       248

    accuracy                         0.7254       845
   macro avg     0.6863    0.7114    0.6924       845
weighted avg     0.7529    0.7254    0.7340       845<jupyter_text>SpamBase 4x4 Pure StateWe build our training density matrix using Pure States.<jupyter_code>phi_train = np.sum(np.abs(Strain), axis=0)
phi_train = phi_train / np.linalg.norm(phi_train)

U_train = create_U_train(phi_train)

print(U_train)<jupyter_output>[[ 0.37950892  0.47081098  0.63630529 -0.47898392]
 [ 0.24651561 -0.12507789 -0.6121573  -0.74084345]
 [ 0.77958893 -0.56123548  0.13372207  0.24366829]
 [ 0.43294818  0.66910932 -0.44999715  0.40292823]]<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|00 \rangle$ obtained.<jupyter_code>backend = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 4x4
Spred_ones = []

for i in range(len(Stest_ones)):

    qc = QuantumCircuit(2, 2)
    qc.initialize(Stest_ones[i], [0, 1])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Spred_ones.append((counts['00'])/5000)
    except:
        Spred_ones.append(0)

print(Spred_ones)
# Pure state case
Spred_zeros = []

for i in range(len(Stest_zeros)):

    qc = QuantumCircuit(2, 2)
    qc.initialize(Stest_zeros[i], [0, 1])
    qc.isometry(U_train.T, [], [0, 1]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Spred_zeros.append((counts['00'])/5000)
    except:
        Spred_zeros.append(0)

print(Spred_zeros)<jupyter_output>[0.1994, 0.3466, 0.039, 0.1376, 0.2048, 0.1986, 0, 3.1662, 0.1142, 0.3322, 0.3788, 0.3246, 0.2044, 0.3908, 0.0108, 0.7788, 0.0512, 0.5198, 0.4514, 0.0054, 0.4592, 0.0468, 0.0078, 0.1198, 0.6582, 0.0654, 0.7822, 0.5348, 0.0454, 0.0518, 0.8216, 0.2322, 0.0078, 0.383, 0.4154, 0.2672, 0.3358, 0.6744, 0.0556, 0.0108, 1.8162, 0.7796, 0.0116, 0.0024, 1.0976, 0.3788, 1.0236, 1.3672, 0.0376, 0.0674, 0.0382, 0.0462, 0.114, 0.7316, 0.4502, 0.4874, 0.0038, 1.3668, 0.1206, 1.1388, 0.0008, 0.3284, 0.0844, 0.0236, 0.5538, 0.348, 0.9748, 0.2622, 0.0938, 0.0684, 0.2568, 0.3398, 0.9408, 0.468, 0.7904, 0.3754, 0.5616, 0.6876, 1.0208, 1.2268, 0.8996, 0.9942, 0.0446, 0.4912, 0.1658, 0.0508, 0.1764, 0.025, 0.1712, 0.3882, 0.0278, 0.0058, 0.3816, 0, 0.1056, 0.0056, 0.505, 0.004, 0.9772, 1.1118, 0.0572, 0.19, 0.2826, 0.353, 0.3726, 0.0214, 0.0558, 0.4222, 1.244, 0.593, 0.4036, 0.4236, 0.0794, 0.0062, 0.6442, 0.0396, 0.0096, 0.0204, 0.0014, 0.1398, 0.5142, 0.3516, 0.2936, 0.011, 0.082, 0.3366, 0.3892, 0.0172, [...]<jupyter_text>We need to determine a threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier.<jupyter_code>Spred_val = np.concatenate((Spred_ones[:100], Spred_zeros[:100]))
print(Spred_val.shape)

Sy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Sy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Sy_val, Spred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Spred_val > t).astype(int)
    metric = f1_score(Sy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.5198<jupyter_text>We have obtained that our best threshold is at $t = 0.5198$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Spr_zeros = (np.array(Spred_zeros) > T).astype(int)

Spr_ones = (np.array(Spred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(248), np.zeros(597))), 
        np.concatenate((Spr_ones[100:], Spr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.8462    0.7370    0.7878       597
         1.0     0.5169    0.6774    0.5864       248

    accuracy                         0.7195       845
   macro avg     0.6815    0.7072    0.6871       845
weighted avg     0.7495    0.7195    0.7287       845<jupyter_text>SpamBase 8x8 Mixed StateWe use a new RFF sampler to reduce the data to arrays of lenght 8. Then we split our data as before.<jupyter_code># Create the RandomFourierFeature map
dim = 8
rm_state = 2
sigma = 1  #dataset dependant
gamma = 1 / (2 * sigma**2)

feature_map_fourier = RBFSampler(gamma=gamma, n_components=dim, random_state=rm_state)
X_feat_train = feature_map_fourier.fit(spam_X)

rffmap = QFeatureMap_rff(rbf_sampler=feature_map_fourier)

Srff = rffmap.batch2wf(spam_X)

print(Srff.shape)
Sones = []
Strain = []
Stest_zeros = []
Stest_ones = []

for i in range(Srff.shape[0]):
    if (spam_y[i] == 1):
        Sones.append(Srff[i] / np.linalg.norm(Srff[i]))
    elif (spam_y[i] == 0):
        Stest_zeros.append(Srff[i] / np.linalg.norm(Srff[i]))

Stest_ones = Sones[:len(Stest_zeros)//2]
Strain = Sones[len(Stest_zeros)//2:]

print(len(Strain))
print(len(Stest_zeros))
print(len(Stest_ones))
print(Strain[0])
print(sum(Strain[0]**2))<jupyter_output>[-0.34027012 -0.40001974 -0.05008387 -0.25346561  0.53922711  0.19674345
  0.2569611  -0.51180502]
1.0<jupyter_text>We build our training density matrix using Mixed States.<jupyter_code>Z_train = np.outer(np.abs(Strain[0]), np.abs(Strain[0]))
for i in range(1, len(Strain)):
    Z_train += np.outer(np.abs(Strain[i]), np.abs(Strain[i]))

Z_train *= 1/len(Strain)
print(Z_train)
lambda_P1_temp, U_train = np.linalg.eigh(Z_train)
lambda_train = round_lambda(lambda_P1_temp)
print(lambda_train)
print(U_train)
lambda_train = lambda_train / np.linalg.norm(lambda_train)
print(lambda_train)
print((sum(lambda_train**2)))<jupyter_output>[0.00412226 0.00503751 0.00747034 0.7330724  0.01024042 0.01468074
 0.0190798  0.67957505]
1.0<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|000 \rangle$ obtained. This value will be proportional to the density estimated by the circuit.<jupyter_code>backend2 = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 8x8
Spred_ones = []

for i in range(len(Stest_ones)):

    qc = QuantumCircuit(6, 3)

    qc.initialize(Stest_ones[i], [0, 1, 2])
    qc.initialize(lambda_train, [3, 4, 5])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.cnot(5, 2)
    qc.cnot(4, 1)
    qc.cnot(3, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)

    counts = execute(qc, backend2, shots=20000).result().get_counts()
    Spred_ones.append((counts['000'])/5000)

print(Spred_ones)
# Mixed state case
Spred_zeros = []

for i in range(len(Stest_zeros)):

    qc = QuantumCircuit(6, 3)

    qc.initialize(Stest_zeros[i], [0, 1, 2])
    qc.initialize(lambda_train, [3, 4, 5])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.cnot(5, 2)
    qc.cnot(4, 1)
    qc.cnot(3, 0)

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)
    counts = execute(qc, backend2, shots=20000).result().get_counts()
    try:
        Spred_zeros.append((counts['000'])/5000)
    except:
        Spred_zeros.append(0)

print(Spred_zeros)<jupyter_output>[0.0616, 0.02, 0.0596, 0.1126, 0.1856, 0.0788, 0.8884, 0.144, 0.043, 0.0278, 0.232, 0.0136, 0.014, 0.2666, 0.1744, 0.1472, 0.4238, 0.0456, 0.2086, 0.4572, 0.0384, 0.3562, 0.21, 0.1102, 0.1384, 0.201, 0.0196, 0.0486, 0.3134, 0.2972, 0.0154, 0.0126, 0.8714, 0.0606, 0.0458, 0.0186, 0.4808, 0.1086, 0.1412, 0.3704, 0.0092, 0.261, 0.0056, 0.116, 0.0346, 0.0004, 0.649, 0.057, 0.016, 0.6096, 0.0406, 0.1456, 0.119, 0.3538, 0.0214, 0.1672, 0.0782, 0.0222, 0.3208, 0.2108, 0.3198, 0.4384, 0.0044, 0.005, 0.0086, 0.6886, 0.2306, 0.0202, 0.1014, 0.0908, 0.1164, 0.0602, 0.1032, 0.1862, 0.0034, 0.0544, 0.0064, 0.0208, 0.8406, 0.029, 0.1106, 0.4398, 0.096, 0.036, 0.2546, 0.0826, 0.0166, 0.223, 0.2656, 0.1222, 0.1074, 0.1356, 0.044, 0.1108, 0.0002, 0.0316, 0.0894, 0.5982, 0.7574, 0.1582, 0.1666, 0.102, 0.0444, 0.1542, 0.1554, 0.0612, 0.346, 0.1168, 0.1824, 0.02, 0.0638, 0.0582, 0.0248, 0.0448, 0.0314, 0.095, 0.4242, 0.0614, 0.0744, 0.0632, 0.0806, 0.115, 0.0314, 0.0638, 0.1242, 0.4748, 0.006, 0.022, 0.22[...]<jupyter_text>We need to determine a threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier.<jupyter_code>Spred_val = np.concatenate((Spred_ones[:100], Spred_zeros[:100]))
print(Spred_val.shape)

Sy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Sy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Sy_val, Spred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Spred_val > t).astype(int)
    metric = f1_score(Sy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.1222<jupyter_text>We have obtained that our best threshold is at $t = 0.1222$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Spr_zeros = (np.array(Spred_zeros) > T).astype(int)

Spr_ones = (np.array(Spred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(248), np.zeros(597))), 
        np.concatenate((Spr_ones[100:], Spr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.7620    0.5846    0.6616       597
         1.0     0.3592    0.5605    0.4378       248

    accuracy                         0.5775       845
   macro avg     0.5606    0.5725    0.5497       845
weighted avg     0.6438    0.5775    0.5959       845<jupyter_text>SpamBase 8x8 Pure StateWe build our training density matrix using Pure States.<jupyter_code>phi_train = np.sum(np.abs(Strain), axis=0)
phi_train = phi_train / np.linalg.norm(phi_train)

U_train = create_U_train(phi_train)

print(U_train)<jupyter_output>[[ 2.42533676e-01  3.92465306e-01 -6.70307607e-01  2.91091225e-01
   2.67895241e-01 -4.20590156e-01 -6.43884790e-02 -1.70937208e-02]
 [ 3.96977281e-01 -1.22555550e-01  3.64045805e-02  2.60429951e-02
  -6.95863502e-01 -2.91149292e-01 -5.06351380e-01 -3.63535790e-05]
 [ 2.06555201e-01  1.73205406e-01 -2.28882531e-01  7.00712557e-02
  -2.51969816e-01  5.40217153e-01  1.42867681e-01 -7.03066929e-01]
 [ 2.22028386e-01  5.51803480e-01  6.94026834e-01  2.93379612e-01
   1.23616632e-01 -1.80470956e-01  3.93990677e-02 -1.70494422e-01]
 [ 4.33649605e-01  3.06048557e-01 -5.08853786e-03 -7.41394147e-01
   2.34198061e-01  2.05822626e-01 -2.12803105e-01  1.61537237e-01]
 [ 3.73541690e-01 -5.31680987e-01  1.05928774e-01  2.79501890e-01
   5.36466451e-01  1.51244397e-01 -3.80668763e-01 -1.81272982e-01]
 [ 3.71143602e-01  6.25581780e-02 -4.59237958e-02  3.98627595e-01
  -1.48527718e-01  4.65894247e-01  2.29219143e-01  6.36919617e-01]
 [ 4.75701431e-01 -3.40881226e-01  4.41153768e-02 -1.92134164e-01
  -[...]<jupyter_text>We will need a different circuit for each sample $|x \rangle$ in the test partition. To obtain a more accurate predicted value, we will run the circuit many shots, and then we will store the amount of $|00 \rangle$ obtained.<jupyter_code>backend = Aer.get_backend('qasm_simulator')
np.random.seed(1234)
## Density Estimation Approach matrices 8x8
Spred_ones = []

for i in range(len(Stest_ones)):

    qc = QuantumCircuit(3, 3)
    qc.initialize(Stest_ones[i], [0, 1, 2])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Spred_ones.append((counts['000'])/5000)
    except:
        Spred_ones.append(0)

print(Spred_ones)
# Pure state case
Spred_zeros = []

for i in range(len(Stest_zeros)):

    qc = QuantumCircuit(3, 3)
    qc.initialize(Stest_zeros[i], [0, 1, 2])
    qc.isometry(U_train.T, [], [0, 1, 2]) # ArbRot as a isometry

    qc.measure(0, 0)
    qc.measure(1, 1)
    qc.measure(2, 2)

    counts = execute(qc, backend, shots=20000).result().get_counts()
    try:
        Spred_zeros.append((counts['000'])/5000)
    except:
        Spred_zeros.append(0)

print(Spred_zeros)<jupyter_output>[0.0636, 0.023, 0.036, 0.0396, 0.0112, 0.041, 0.0174, 0.1252, 0.0438, 0.0006, 0.0216, 0.0002, 0.003, 0.2998, 0.0648, 0.0756, 0.2608, 0.0772, 0.0364, 0.0424, 0.0684, 0.0008, 0.0062, 0.0604, 0.086, 0.0948, 0.0302, 0.0756, 0.007, 0, 0.0002, 0.025, 0.001, 0.1122, 0.0026, 0.0274, 0.1362, 0.0386, 0.0016, 0.057, 0.007, 0.0724, 0.0094, 0.0014, 0, 0, 0.1166, 0.0278, 0.015, 0.1096, 0.0026, 0.0008, 0.0004, 0.1416, 0.002, 0, 0.0194, 0.001, 0.055, 0.447, 0.0544, 0.1132, 0.0042, 0.0004, 0, 0.0724, 0.0346, 0.009, 0.001, 0.067, 0.1576, 0.0238, 0.2234, 0.0322, 0.0076, 0.111, 0.0118, 0.0376, 0.0478, 0.0006, 0.0262, 0.0474, 0.0154, 0.0616, 0.0252, 0.0334, 0.0014, 0.0342, 0.0194, 0.0012, 0.026, 0.0138, 0.004, 0.2384, 0, 0.0604, 0.1152, 0.285, 0.0242, 0.2776, 0.0604, 0.171, 0, 0.0454, 0.053, 0.079, 0.0182, 0.0002, 0.291, 0.0052, 0.1194, 0.112, 0.0028, 0.0824, 0.0208, 0.0724, 0.5996, 0.0584, 0.0632, 0.0004, 0.0036, 0.0126, 0.0038, 0.0278, 0.0212, 0.99, 0.001, 0.0092, 0.1216, 0.0172, 0.0038, 0.0948, 0.0184, [...]<jupyter_text>Once we have obtained the predicted density values, we need to determine a threshold $t$ such that if a sample has a density lower than $t$, it is classified as an outlier. For this, we will define a validation partition from our test groups, and we'll try to find the best possible value for $t$ using only validation samples. After that, we'll use this threshold to compare with all the remaining test samples.<jupyter_code>Spred_val = np.concatenate((Spred_ones[:100], Spred_zeros[:100]))
print(Spred_val.shape)

Sy_val = np.concatenate((np.ones(100), np.zeros(100)))
print(Sy_val.shape)
from sklearn.metrics import roc_curve, auc

fpr, tpr, thresholds = roc_curve(Sy_val, Spred_val, pos_label=1)

plt.figure()
lw = 2
plt.plot(fpr, tpr, color="darkorange",
    lw=lw, label="ROC curve")
plt.plot([0, 1], [0, 1], color="navy", lw=lw, linestyle="--")
print(auc(fpr, tpr))
from sklearn.metrics import f1_score

m, T = 0,0
for t in thresholds:
    preds = (Spred_val > t).astype(int)
    metric = f1_score(Sy_val, preds, average='weighted')
    if (metric > m):    
        m = metric
        T = t

print(T)<jupyter_output>0.0616<jupyter_text>We have obtained that our best threshold is at $t = 0.0616$. So we will use this threshold to compare all test samples and to calculate some metrics for our estimator.<jupyter_code>Spr_zeros = (np.array(Spred_zeros) > T).astype(int)

Spr_ones = (np.array(Spred_ones) > T).astype(int)
from sklearn.metrics import classification_report

print(classification_report(
        np.concatenate((np.ones(248), np.zeros(597))), 
        np.concatenate((Spr_ones[100:], Spr_zeros[100:])), digits=4))<jupyter_output>precision    recall  f1-score   support

         0.0     0.6988    0.6566    0.6770       597
         1.0     0.2782    0.3185    0.2970       248

    accuracy                         0.5574       845
   macro avg     0.4885    0.4876    0.4870       845
weighted avg     0.5753    0.5574    0.5655       845
{'fraction_non_alphanumeric': 0.14771920572784616, 'fraction_numerical': 0.17652069853761435, 'mean_word_length': 3.672479150871873, 'pattern_counts': {'":': 0, '<': 133, '<?xml version=': 0, '>': 167, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 5, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '42', 'dir': 'jupyter-structured-clean-dedup', 'id': '474994', 'license': 'Apache-2.0', 'n_tokens_mistral': 26194, 'n_tokens_neox': 20560, 'n_words': 4928, 'path': 'AnomalyDetection_with_RFF.ipynb', 'repo_name': 'diegour1/QuantumAnomalyDetection', 'size': '49313'}