import numpy as np ; import matplotlib.pyplot as plt
# import direct connaissant les entêtes (N°dep,longitude,latitude,Nbre habitants)
ardeche = np.loadtxt('ardeche_sup_1000.csv', delimiter=';', skiprows=1, dtype=str)
ardeche = np.char.replace(ardeche, ',', '.'); ardeche=ardeche.astype(float)
# initialise la classe à 0 pour chaque ville et initialise arbitrairement 3 centroïde
colk=np.ones((len(ardeche[:,1]),1))
Cent=np.array([[4.2,44.4],[4.5,44.8],[4.8,45.2]])

def affichage(Dep,Clas,CdG,n): # Affichage des points et centroïdes
    plt.figure(n)
    plt.scatter(CdG[:,0],CdG[:,1],c="red",s=150);plt.scatter(Dep[:,1],Dep[:,2],c=Clas);
    plt.xlabel("Longitude en °");plt.ylabel("Latitude en °");
    plt.title("Villes de plus de 1000 habitants");plt.grid(True);plt.show()

def classe(Dep,Clas,CdG): # Attribution de sa classe à chaque ville
    for i in range(0,len(Dep[:,1]-1)):
        Clas[i]=1
        d1=np.sqrt((Dep[i,1]-CdG[0,0])**2+(Dep[i,2]-CdG[0,1])**2)
        d2=np.sqrt((Dep[i,1]-CdG[1,0])**2+(Dep[i,2]-CdG[1,1])**2)
        d3=np.sqrt((Dep[i,1]-CdG[2,0])**2+(Dep[i,2]-CdG[2,1])**2)
        if d2<=d1 :
            Clas[i]=2
        if d3<=d2 :
            Clas[i]=3
    return(Clas)

def nouvcent(Dep,Clas): # calcul du centroïde d'une classe (barycentre)
    id1=0;id2=0;id3=0;
    X1,Y1,X2,Y2,X3,Y3=0,0,0,0,0,0;CdG1=[X1,Y1];CdG2=[X2,Y2];CdG1=[X3,Y3]
    for i in range(0,len(ardeche[:,1])-1):
        if Clas[i]==1:
            X1=X1+Dep[i,1];Y1=Y1+Dep[i,2];id1=id1+1
        elif Clas[i]==2:
            X2=X2+Dep[i,1];Y2=Y2+Dep[i,2];id2=id2+1
        else:
            X3=X3+Dep[i,1];Y3=Y3+Dep[i,2];id3=id3+1
    X1=X1/id1;Y1=Y1/id1;X2=X2/id2;Y2=Y2/id2;X3=X3/id3;Y3=Y3/id3 # pb de divpar 0 si cluster vide
    CdG=np.array([[X1,Y1],[X2,Y2],[X3,Y3]])
    return(CdG)

affichage(ardeche,colk,Cent,1)# affiche carte et centroïdes initiaux
colk=classe(ardeche,colk,Cent) # calcul de la classe de chaque ville
affichage(ardeche,colk,Cent,2) # affiche classe de chaque ville et de son centroïde

i=0;stop=0
while i<50 and stop!=1: # calcul et affichage classe et centroïde tant qu'il bouge
    Cent1=nouvcent(ardeche,colk)
    if np.array_equal(Cent1,Cent):
        stop=1
    Cent=Cent1 ; affichage(ardeche,colk,Cent,3+i) ; colk=classe(ardeche,colk,Cent) ; i+=1