import numpy as np
import matplotlib.pyplot as plt
np.random.seed(26)
X = np.random.normal(0,1.8,[180,2])
X[60:120] += 8
X[120:,0] += 16
plt.axes(aspect=1)
plt.scatter(X[:,0],X[:,1],c='m',edgecolor='k')
plt.show()

n_krachuk = 3 # จำนวนกระจุก
n_thamsam = 100 # จำนวนทำซ้ำสูงสุด
tol = 0.0001 # ค่าความเปลี่ยนแปลงสูงสุดที่ยอมให้หยุดได้
sumlueak = np.random.choice(len(X),n_krachuk,replace=0)
X_cen = X[sumlueak] # จุดเซนทรอยด์ตั้งต้น เลือกแบบสุ่ม
# วนซ้ำเพื่อปรับเซนทรอยด์
for i in range(n_thamsam):
    raya2 = ((X_cen[None]-X[:,None])**2).sum(2) # วัดระยะห่างจากจุดถึงเซนทรอยด์
    klum = raya2.argmin(1) # ตัดสินกลุ่มของจุดโดยเลือกเซนทรอยด์ที่ใกล้สุด
    X_cen_mai = np.empty_like(X_cen) # จุดเซนทรอยด์ใหม่
    # วนซ้ำเพื่อหาตำแหน่งเซนทรอยด์ใหม่
    for j in range(n_krachuk):
        if(len(X[klum==j])): # ถ้ามีสมาชิกในกลุ่ม
            X_cen_mai[j] = X[klum==j].mean(0) # กำหนดเซนทรอยด์ใหม่เป็นตำแหน่งเฉลี่ยของทุกจุดในกลุ่ม
        else: # ถ้าในกลุ่มว่างเปล่าก็ให้สุ่มเซนทรอยด์ใหม่
            X_cen_mai[j] = X[np.random.randint(len(X))]
    if(np.allclose(X_cen,X_cen_mai,atol=tol)): # ถ้าความเปลี่ยนแปลงน้อยกว่าค่าที่กำหนดก็ให้หยุด
        X_cen = X_cen_mai
        break
    X_cen = X_cen_mai # ย้ายจุดเซนทรอยด์ไปยังตำแหน่งใหม่

raya2 = ((X_cen[None]-X[:,None])**2).sum(2)
z = raya2.argmin(1)

plt.axes(aspect=1)
plt.scatter(X[:,0],X[:,1],c=z,edgecolor='k',cmap='rainbow')
plt.scatter(X_cen[:,0],X_cen[:,1],300,'#EEAA55',marker='*',edgecolor='#9999DD',lw=2)
plt.show()

sse = 0
for i in range(n_krachuk):
    sse += np.sum(raya2[z==i,i])
print(sse) # ได้ 1106.13721312

sse = (raya2*(z[:,None]==np.arange(n_krachuk))).sum()

np.random.seed(4)
X = np.random.normal(0,3.,[480,2])**3/60
X[120:240] += 5
X[360:480,0] += 12
X[360:480,1] += 7
X[240:360,0] += 16
X[240:360,1] += 2
plt.axes(aspect=1)
plt.scatter(X[:,0],X[:,1],c='g',alpha=0.6,edgecolor='k')
plt.show()

n_krachuk = 4
n_thamsam = 50
tol = 0.00001
n_sumchut = 10 # จำนวนครั้งที่สุ่มจุดใหม่
sse_noisut = np.inf # ค่า sse น้อยสุด
for s in range(n_sumchut):
    sumlueak = np.random.choice(len(X),n_krachuk,replace=0)
    X_cen = X[sumlueak]
    for i in range(n_thamsam):
        raya2 = ((X_cen[None]-X[:,None])**2).sum(2)
        klum = raya2.argmin(1)
        X_cen_mai = np.empty_like(X_cen)
        for j in range(n_krachuk):
            if(len(X[klum==j])):
                X_cen_mai[j] = X[klum==j].mean(0)
            else:
                X_cen_mai[j] = X[np.random.randint(len(X))]
        if(np.allclose(X_cen,X_cen_mai,atol=tol)):
            X_cen = X_cen_mai
            break
        X_cen = X_cen_mai

    raya2 = ((X_cen[None]-X[:,None])**2).sum(2)
    klum = raya2.argmin(1)

    # หา sse ในแต่ละรอบ ถ้าได้น้อยกว่าเดิมก็เก็บค่านั้นและตำแหน่งเซนทรอยด์ไว้
    sse = (raya2*(klum[:,None]==np.arange(n_krachuk))).sum()
    print(sse)
    if(sse_noisut>sse):
        sse_noisut = sse
        X_cen_disut = X_cen

raya2 = ((X_cen_disut[None]-X[:,None])**2).sum(2)
z = raya2.argmin(1)

plt.axes(aspect=1)
plt.scatter(X[:,0],X[:,1],c=z,edgecolor='k',cmap='rainbow')
plt.scatter(X_cen_disut[:,0],X_cen_disut[:,1],300,'#EEAA55',alpha=0.8,marker='*',edgecolor='#9999DD',lw=2)
plt.show()

3696.79020957
3814.36451799
1432.78366256
1432.75873231
3853.71442022
1432.78366256
1432.78366256
3756.04317335
4290.05406862
3756.04317335

class Kmeans:
    def __init__(self,n_krachuk,n_sumchut=10):
        self.n_krachuk = n_krachuk
        self.n_sumchut = n_sumchut
    
    def rianru(self,X,n_thamsam=100,tol=0.0001):
        sse_noisut = np.inf
        for s in range(self.n_sumchut):
            sumlueak = np.random.choice(len(X),self.n_krachuk,replace=0)
            X_cen = X[sumlueak]
            for i in range(n_thamsam):
                raya2 = ((X_cen[None]-X[:,None])**2).sum(2)
                klum = raya2.argmin(1)
                X_cen_mai = np.empty_like(X_cen)
                for j in range(self.n_krachuk):
                    if(len(X[klum==j])):
                        X_cen_mai[j] = X[klum==j].mean(0)
                    else:
                        X_cen_mai[j] = X[np.random.randint(len(X))]
                if(np.allclose(X_cen,X_cen_mai,atol=tol)):
                    break
                X_cen = X_cen_mai
            
            X_cen = X_cen_mai
            raya2 = ((X_cen[None]-X[:,None])**2).sum(2)
            klum = raya2.argmin(1)
            sse = (raya2*(klum[:,None]==np.arange(self.n_krachuk))).sum()
            if(sse_noisut>sse):
                sse_noisut = sse
                self.X_cen = X_cen
    
    def thamnai(self,X):
        raya2 = ((self.X_cen[None]-X[:,None])**2).sum(2)
        return raya2.argmin(1)

np.random.seed(26)
X = np.random.normal(0,1.8,[180,2])
X[60:120] += 8
X[120:,0] += 16
km = Kmeans(n_krachuk=3)
km.rianru(X)
z = km.thamnai(X)
plt.axes(aspect=1)
plt.scatter(X[:,0],X[:,1],c=z,edgecolor='k',cmap='rainbow')
plt.scatter(km.X_cen[:,0],km.X_cen[:,1],300,'#EEAA55',marker='*',edgecolor='#9999DD',lw=2)
plt.show()

def plotkmeans(n,X):
    km = Kmeans(n)
    km.rianru(X)
    z = km.thamnai(X)
    plt.figure()
    plt.axes(aspect=1)
    plt.scatter(X[:,0],X[:,1],c=z,edgecolor='k',cmap='rainbow')
    plt.scatter(km.X_cen[:,0],km.X_cen[:,1],300,'#EEAA55',marker='*',edgecolor='#9999DD',lw=2)
    plt.show()

X = np.random.uniform(-0.5,0.5,[840,2])
plotkmeans(7,X)

X = np.random.normal(0,1,[840,2])
plotkmeans(7,X)

from sklearn import datasets
X,_ = datasets.make_moons(n_samples=180,noise=0.05)
plotkmeans(2,X)