import numpy as np
import math
# Importing standard Qiskit libraries
from qiskit import QuantumCircuit, transpile, Aer, IBMQ
from qiskit.tools.jupyter import *
from qiskit.quantum_info import *
from qiskit.visualization import *
from ibm_quantum_widgets import *
from qiskit.providers.aer import QasmSimulator
from qiskit.extensions import Initialize
from qiskit.providers.ibmq.managed import IBMQJobManager
# Loading your IBM Quantum account(s)
provider = IBMQ.load_account()
import matplotlib.pyplot as plt

backend=provider.get_backend('ibmq_qasm_simulator')

data=np.load("data.npy")
print("No of datapoints:", data.shape[1])
print("Dimension:", data.shape[0])

def compute_mean(x,d,k):                #helperfunction to compute cluster means(used in all the follwing subparts)
    l=[]
    for i in range(k):
        a=np.mean(x,axis=1,where=x[d-1]==i+1)[:d-1] 
        a[np.isnan(a)]=0                                    #handling empty slice and true dvide warnings
        l.append(a)                             

    return np.array(l)

def quant_clustering(data,n=400,k=4,max_iter=20):
    dim=data.shape[0]
    means=np.load("means.npy")
    done=False
    ite=0
    z=list(np.random.choice([1,2,3,4],n))
    
    error=[]
    X_clus=np.vstack((data,z))
    dim_new=X_clus.shape[0]
    err= np.sum(np.linalg.norm(X_clus[:dim_new-1].transpose()-means[(X_clus[dim_new-1]-1).astype("int64")],axis=1)**2,axis=0)
    error.append(err)
    
    while (not done) and ite<=max_iter:
        done=True
        circ=[]
        for i in range(n):
            for j in range(k):
                circuit=QuantumCircuit(12,1)
                dt=np.hstack((data[:,i],np.zeros(2**9-len(data[:,i]))))
                mean=np.hstack((means[j],np.zeros(2**9-len(means[j]))))
                
                a_norm=(dt/np.linalg.norm(dt)) if np.linalg.norm(dt)!=0 else np.array([0]*2**9)
                b_norm=(mean/np.linalg.norm(mean)) if np.linalg.norm(mean)!=0 else np.array([0]*2**9)
                
                c=np.hstack((a_norm,b_norm))/ 2**0.5
                c_norm=c/np.linalg.norm(c)
                
                circuit.initialize(c_norm,range(10))
                Z=(np.linalg.norm(dt)**2 + np.linalg.norm(mean)**2)
                circuit.initialize([Z**(-0.5) * np.linalg.norm(dt),-Z**(-0.5) * np.linalg.norm(mean)],11)
                circuit.h(11)
                circuit.cswap(11,10,9)
                circuit.h(11)
                circuit.measure(11,0)
                circ.append(circuit)
    
        circ = transpile(circ, backend=backend)
        job_manager = IBMQJobManager()
        job = job_manager.run(circ, backend=backend,name="clustering_ite_{}".format(ite),shots=8192)
        result= job.results()
        prob=np.array([result.get_counts(i)["0"]/((result.get_counts(i)["0"] if "0" in result.get_counts(i) else 0)
                       +(result.get_counts(i)["1"] if "1" in result.get_counts(i) else 0)) for i in range(n*k)])
        distances=list(4*Z*abs(prob-0.5))
        
        for i in range(n):
            candidate=np.argmin(distances[k*i:k*i+k])+1
            if z[i]!=candidate:
                done=False
                z[i]=candidate
                
        err=0
        for i in range(n):
            err+=(distances[4*i+(z[i]-1)])
        error.append(err)
            
        means=compute_mean(np.vstack((data,z)),len(data)+1,k)
        ite+=1
        
        if ite==20 and done==False:
            print("Not Converged")
            print("Final_value:", np.round(error[-1],2))
        elif done==True:
            print("Convergence:", np.round(error[-1],2))
        
    plt.plot(error)
    plt.grid(True)
    plt.xlabel("Iteration number")
    plt.ylabel("Objective function")
    plt.title("Quantum")
    plt.show()
    

    return z

z1=quant_clustering(data)
np.save("cluster1.npy",np.array(z1))

from collections import Counter
def plot_result(z,title):
    plt.title(title)
    c=Counter(z).most_common()
    cx=[x for x,y in c]
    cy=[y for x,y in c]
    plt.xticks(np.arange(0,5,1))
    plt.yticks(np.arange(0,100,5))
    plt.grid(True)
    plt.xlabel("cluster number")
    plt.ylabel("count")
    plt.bar(cx,cy)
    plt.show()
    
plot_result(z1[:100],"covid")
plot_result(z1[100:200],"normal")
plot_result(z1[200:300],"opacity")
plot_result(z1[300:],"viral")

cluster_stat=list(enumerate(z1))
cluster_stat.sort(key= lambda x: x[1])
cluster1=[]
cluster2=[]
cluster3=[]
cluster4=[]
a=[0]*4
b=[0]*4
c=[0]*4
d=[0]*4
for i,j in cluster_stat:
    if j==1:
        if i<100:
            a[0]+=1
        elif i>=100 and i<200:
            a[1]+=1
        elif i>=200 and i<300:
            a[2]+=1
        else:
            a[3]+=1
    if j==2:
        if i<100:
            b[0]+=1
        elif i>=100 and i<200:
            b[1]+=1
        elif i>=200 and i<300:
            b[2]+=1
        else:
            b[3]+=1
    if j==3:
        if i<100:
            c[0]+=1
        elif i>=100 and i<200:
            c[1]+=1
        elif i>=200 and i<300:
            c[2]+=1 
        else:
            c[3]+=1
            
    if j==4:
        if i<100:
            d[0]+=1
    
        elif i>=100 and i<200:
            d[1]+=1

        elif i>=200 and i<300:
            d[2]+=1
        else:
            d[3]+=1
            
def plot_result2(a,title):
    plt.title(title)
    plt.grid(True)
    plt.xticks(np.arange(0,5,1))
    plt.yticks(np.arange(0,100,5))
    plt.bar(range(4),a)
    plt.xlabel("type")
    plt.ylabel("count")
    plt.show()
    
plot_result2(a,"cluster 1")
plot_result2(b,"cluster 2")
plot_result2(c,"cluster 3")
plot_result2(d,"cluster 4")

means=np.load("means.npy")
circuit=QuantumCircuit(12,1)
dt=np.hstack((data[:,0],np.zeros(2**9-len(data[:,0]))))
mean=np.hstack((means[1],np.zeros(2**9-len(means[1]))))
                
a_norm=(dt/np.linalg.norm(dt)) if np.linalg.norm(dt)!=0 else np.array([0]*2**9)
b_norm=(mean/np.linalg.norm(mean)) if np.linalg.norm(mean)!=0 else np.array([0]*2**9)
                
c=np.hstack((a_norm,b_norm))/ 2**0.5
c_norm=c/np.linalg.norm(c) 
                
circuit.initialize(c_norm,range(10))
Z=(np.linalg.norm(dt)**2 + np.linalg.norm(mean)**2)
circuit.initialize([Z**(-0.5) * np.linalg.norm(dt),-Z**(-0.5) * np.linalg.norm(mean)],10)
circuit.h(11)
circuit.cswap(11,10,9)
circuit.h(11)
circuit.measure(11,0)
circuit.draw("latex")

circuit = transpile(circuit, backend=backend)
job_manager = IBMQJobManager()
job = job_manager.run([circuit], backend=backend,name="clustering_ite_{}".format(1),shots=8192)
result=job.results()

