import numpy as np

# Importing standard Qiskit libraries
from qiskit import QuantumCircuit, transpile, Aer, IBMQ
from qiskit.tools.jupyter import *
from qiskit.visualization import *
from ibm_quantum_widgets import *
from qiskit.providers.aer import QasmSimulator

# Loading your IBM Quantum account(s)
provider = IBMQ.load_account()

!pip install mlxtend
!pip install lightgbm
!pip install git+https://github.com/qiskit-community/Quantum-Challenge-Grader.git


import pandas as pd

import warnings
warnings.filterwarnings('ignore')

import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
from mlxtend.plotting import plot_decision_regions
import seaborn as sns
import networkx as nx
from IPython.display import clear_output
import plotly.graph_objects as go
from plotly.subplots import make_subplots



from sklearn import datasets
from sklearn.datasets import make_hastie_10_2
from sklearn.ensemble import GradientBoostingClassifier, AdaBoostClassifier, RandomForestClassifier
from sklearn import metrics
from sklearn.metrics import confusion_matrix, accuracy_score, precision_score, recall_score,  f1_score, classification_report
from sklearn.preprocessing import LabelEncoder, StandardScaler, MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn import svm
from sklearn.svm import SVC
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.neural_network import MLPClassifier

import lightgbm as lgbm

from qiskit import Aer, execute
from qiskit.opflow import Z, I, PauliExpectation, CircuitSampler, StateFn
from qiskit.utils import QuantumInstance
from qiskit.circuit import QuantumCircuit, Parameter, ParameterVector
from qiskit.circuit.library import HGate, RXGate, RYGate, RZGate, CXGate, CRXGate, CRZGate 
from qiskit.circuit.library import PauliFeatureMap, ZFeatureMap, ZZFeatureMap
from qiskit.circuit.library import TwoLocal, NLocal, RealAmplitudes, EfficientSU2
from qiskit.algorithms.optimizers import GradientDescent, COBYLA, SLSQP, ADAM
from qiskit.visualization import plot_state_city
from qiskit_machine_learning.kernels import QuantumKernel
from qiskit_machine_learning.algorithms import QSVC
from qc_grader.utilities.graph_util import display_maxcut_widget, QAOA_widget
from qiskit.providers.aer.noise import NoiseModel
from qiskit.providers.aer.noise.errors import pauli_error, depolarizing_error
from qiskit.utils.mitigation import complete_meas_cal, CompleteMeasFitter
from qiskit.opflow.gradients import Gradient, NaturalGradient, QFI, Hessian
#from qiskit.algorithms.MulticlassExtension import AllPairs





iris = datasets.load_iris()

print(iris)

iris.feature_names

iris.target_names

iris_df_data = pd.DataFrame(iris.data, columns = iris.feature_names)
iris_df_target = pd.DataFrame(iris.target, columns = ['Iris Label'])




iris_df_data  # printing iris.data

iris_df_target  # printing iris.target; 0 : Setosa, 1: Versicolour, 2 : Virginica

iris_df_data.describe()  # standard decription of data

iris_df_data.info()   # standard info about data 

iris_df_data.columns  # listing out column names of data

iris_df_data['sepal length (cm)'].hist()  # histogram plotting of data

iris_df_data['sepal width (cm)'].hist()

iris_df_data['petal length (cm)'].hist()

iris_df_data['petal width (cm)'].hist()

iris_df_target['Iris Label'].hist()  # This means all flowers have equal number of samples

# scatterplot of data

irisSB = sns.load_dataset('iris')
sns.set_style("darkgrid")


sns.FacetGrid(irisSB, hue = "species", height = 6).map(plt.scatter, 'sepal_length', 'sepal_width').add_legend()

sns.FacetGrid(irisSB, hue = "species", height = 6).map(plt.scatter, 'petal_length', 'petal_width').add_legend()

sns.FacetGrid(irisSB, hue = "species", height = 6).map(plt.scatter, 'sepal_length', 'petal_length').add_legend()

sns.FacetGrid(irisSB, hue = "species", height = 6).map(plt.scatter, 'sepal_width', 'petal_width').add_legend()

# Pair Plot

sns.pairplot(irisSB, hue = 'species')

# Distplot for data

plot = sns.FacetGrid(irisSB, hue = "species")
plot.map(sns.distplot, "sepal_length").add_legend()
  
plot = sns.FacetGrid(irisSB, hue = "species")
plot.map(sns.distplot, "sepal_width").add_legend()
  
plot = sns.FacetGrid(irisSB, hue = "species")
plot.map(sns.distplot, "petal_length").add_legend()
  
plot = sns.FacetGrid(irisSB, hue = "species")
plot.map(sns.distplot, "petal_width").add_legend()
  
plt.show()

# Boxplot graph of features VS flower types

def graph(y):
    sns.boxplot(x = "species", y = y, data = irisSB)
  
plt.figure(figsize=(10,10))
      
# Adding the subplot at the specified
# grid position
plt.subplot(221)
graph('sepal_length')
  
plt.subplot(222)
graph('sepal_width')
  
plt.subplot(223)
graph('petal_length')
  
plt.subplot(224)
graph('petal_width')
  
plt.show()

# Outliers for the features

sns.boxplot(x = 'sepal_length', data = irisSB)

sns.boxplot(x = 'sepal_width', data = irisSB)

## To remove Sepal Width Outliers using IQR(Interquartile Range)
## Code is contained in Single Hashes!! 

## IQR
#Q1 = np.percentile(irisSB['sepal_width'], 25, interpolation = 'midpoint')
  
#Q3 = np.percentile(irisSB['sepal_width'], 75, interpolation = 'midpoint')

#IQR = Q3 - Q1
  
#print("Old Shape: ", irisSB.shape)
  
## Upper bound
#upper = np.where(irisSB['sepal_width'] >= (Q3+1.5*IQR))
  
## Lower bound
#lower = np.where(irisSB['sepal_width'] <= (Q1-1.5*IQR))
  
## Removing the Outliers
#irisSB.drop(upper[0], inplace = True)
#irisSB.drop(lower[0], inplace = True)
  
#print("New Shape: ", irisSB.shape)
  
#sns.boxplot(x = 'sepal_width', data = irisSB)

sns.boxplot(x = 'petal_length', data = irisSB)

sns.boxplot(x = 'petal_width', data = irisSB)

# Correlation Matrix

iris_df_data.corr()  # Petal Length and Petal Width have high correlation

# Correlation Matrix Heatmap

corr = iris_df_data.corr()
fig, ax = plt.subplots(figsize=(5,4))
sns.heatmap(corr, annot=True, ax=ax, cmap = 'coolwarm')  # cmap = 'darkgrid' for light colours


sepal_length = iris['data'][:, 0]
sepal_width = iris['data'][:, 1]
petal_length = iris['data'][:, 2]
petal_width = iris['data'][:, 3]
target = iris['target']


#X = np.column_stack((sepal_length, sepal_width, petal_length, petal_width))
X = np.column_stack((sepal_length, sepal_width, petal_length, petal_width))
x = np.column_stack((sepal_length, sepal_width, petal_length))  # correlation of petal width and petal length is high, hence reducing variable
Y = target

x_train, x_test, y_train, y_test = train_test_split(X, Y, test_size = 0.2, random_state = 345)

# Reduce dimensions
n_dim = 4
pca = PCA(n_components = n_dim).fit(x_train)
x_train = pca.transform(x_train)
x_test = pca.transform(x_test)

# Normalise
std_scale = StandardScaler().fit(x_train)
x_train = std_scale.transform(x_train)
x_test = std_scale.transform(x_test)

# Scale
samples = np.append(x_train, x_test, axis = 0)
minmax_scale = MinMaxScaler((-1, 1)).fit(samples)
x_train = minmax_scale.transform(x_train)
x_test = minmax_scale.transform(x_test)

C = 1  # Regularlization Parameter
clf_CSVC = make_pipeline(SVC(C = C, decision_function_shape = 'ovo'))   # default kernel is ' RBF '
model_CSVC = clf_CSVC.fit(x_train, y_train)  # CSVC Model Creation

model_CSVC.score(x_test, y_test)   # CSVC Model score against test data

# Plotting Decision regions
# The dimension should be 2 for plotting purposes
# Also choose 2 features of the dataset


#plot_decision_regions(X = x_test, y = y_test, clf = model_CSVC, legend = 1)

# Split dataset
sample_train, sample_test, label_train, label_test = train_test_split(
     x, Y, test_size = 0.2, random_state = 345)

# Reduce dimensions
n_dim = 3
pca = PCA(n_components = n_dim).fit(sample_train)
sample_train = pca.transform(sample_train)
sample_test = pca.transform(sample_test)

# Normalise
std_scale = StandardScaler().fit(sample_train)
sample_train = std_scale.transform(sample_train)
sample_test = std_scale.transform(sample_test)

# Scale
samples = np.append(sample_train, sample_test, axis=0)
minmax_scale = MinMaxScaler((-1, 1)).fit(samples)
sample_train = minmax_scale.transform(sample_train)
sample_test = minmax_scale.transform(sample_test)

# Select
train_size = 120
sample_train = sample_train[:train_size]
label_train = label_train[:train_size]

test_size = 30
sample_test = sample_test[:test_size]
label_test = label_test[:test_size]

# Backend and Quantum Instance

backend = Aer.get_backend('statevector_simulator')
q_instance = QuantumInstance(backend, shots = 8192, seed_simulator = 345, seed_transpiler = 345)

# Feature Map to convert classical data to quantum data, basically a parameterized Quantum Circuit
# Can be used to map seemingly unseparable data to separable form, with respect to hyperplanes, in some higher dimension
# This is achieved by e.g. rotation gate parameters depend on data features, etc.

dim = 3  # Feature Dimension = number of qubits
feature_map = ZZFeatureMap(feature_dimension = dim, reps = 1, entanglement = 'linear', insert_barriers = True, parameter_prefix = 'x').decompose() 
feature_map.draw('mpl', style = 'iqx')

# Quantum Kernel

q_kernel = QuantumKernel(feature_map = feature_map, quantum_instance = q_instance)

q_circuit = q_kernel.construct_circuit(sample_train[0], sample_train[1])
q_circuit.decompose().decompose().draw(output = 'mpl')

# Noise Model : Measurement error


def get_noise(pm,pg):
    error_meas = pauli_error([('X',pm), ('I', 1 - pm)])
    error_gate1 = depolarizing_error(pg,1)
    error_gate2 = error_gate1.tensor(error_gate1)
    
    noise_model = NoiseModel()
    noise_model.add_all_qubit_quantum_error(error_meas, "measure") # measurement error is applied to measurements
    noise_model.add_all_qubit_quantum_error(error_gate1, ["h"]) # single qubit gate error is applied to H gates
    noise_model.add_all_qubit_quantum_error(error_gate1, ["p"]) # single qubit gate error is applied to P gates
    noise_model.add_all_qubit_quantum_error(error_gate2, ["cx"]) # two qubit gate error is applied to CX gates
    
    return noise_model


noise_model = get_noise(0.01,0.01)



# Without Noise


job = execute(q_circuit, backend = backend, shots = 10000, seed_simulator = 345, seed_transpiler = 345)
counts = job.result().get_counts(q_circuit)

plot_histogram(counts)


# With Noise


Noisy_job = execute(q_circuit, backend = backend, shots = 10000, seed_simulator = 345, seed_transpiler = 345, noise_model = noise_model)
Noisy_counts = Noisy_job.result().get_counts(q_circuit)

plot_histogram(Noisy_counts)


# Error mitigation
meas_calibs, state_labels = complete_meas_cal(qr = 3, circlabel =' mcal')


cal_results = execute(meas_calibs, backend = backend, shots = 10000, seed_simulator = 345, seed_transpiler = 345, noise_model=noise_model).result()

meas_fitter = CompleteMeasFitter(cal_results, state_labels, circlabel='mcal')
meas_filter = meas_fitter.filter

mitigated_results = meas_filter.apply(Noisy_job.result())
mitigated_counts = mitigated_results.get_counts(q_circuit)

plot_histogram(mitigated_counts)



plot_histogram([counts, Noisy_counts, mitigated_counts], legend = ['Without Noisy','Noisy', 'Mitigated'])

# As we can see from the plot, the noise will not affect the score a lot, as the plots are similar in nature

# Kernel Matrix Evaluation

matrix_train = q_kernel.evaluate(x_vec = sample_train)
matrix_test = q_kernel.evaluate(x_vec = sample_test, y_vec = sample_train)

fig, axs = plt.subplots(1, 2, figsize=(10, 5))

axs[0].imshow(np.asmatrix(matrix_train), interpolation = 'nearest', origin = 'upper', cmap = 'Blues')
axs[0].set_title("training kernel matrix")

axs[1].imshow(np.asmatrix(matrix_test), interpolation = 'nearest', origin = 'upper', cmap = 'Reds')
axs[1].set_title("testing kernel matrix")

plt.show()

qpc_svc = make_pipeline(SVC(kernel='precomputed', C = 100))   # QSVC precomputed
qpc_svc.fit(matrix_train, label_train)
qpc_score = qpc_svc.score(matrix_test, label_test)

print(f'Precomputed kernel classification test score: {qpc_score}')   # best possible score we got after after trying many combinations

""" NOTE : On using the multiclass extension function = " AllPairs " on a previous version of Qiskit, we get score of 1!
                     Hope they introduce it back! """

""" This work of ours is a result of 'Qiskit Fall Fest Kolkata Chapter Hackathon', hosted by Mr. Ritajit Majumdar and Ms. Debasmita Bhowmik 
 from Indian Statistical Institute, in collaboration with IBM Quantum"""

# Authors
# Vismay Joshi, National Institute Of Technology, Agartala
# Tanjin Adnan Abir, Jahangirnagar University
# Abhishek Rawat, TCG CREST
# Suprabhat Sinha, D Y Patil International University
# Rajat Lakhera, Sri Dev Suman Uttarakhand University

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright