import numpy as np
rng=np.random.default_rng(420)
# Ongeëtiketteerde, gestandaardiseerde modelsamples.
X=np.vstack([rng.normal(c,.45,(70,2)) for c in [(-1,-1),(1,0),(0,1.5)]])
def kmeans(X,k,seed):
    r=np.random.default_rng(seed)
    c=X[r.choice(len(X),k,replace=False)].copy()
    for it in range(100):
        lab=((X[:,None,:]-c[None,:,:])**2).sum(axis=2).argmin(axis=1)
        nieuw=np.array([X[lab==j].mean(axis=0) if np.any(lab==j)
                        else X[r.integers(len(X))] for j in range(k)])
        if np.max(np.abs(nieuw-c))<1e-8: c=nieuw; break
        c=nieuw
    lab=((X[:,None,:]-c[None,:,:])**2).sum(axis=2).argmin(axis=1)
    sse=np.sum((X-c[lab])**2)
    return sse,c,lab,it+1
sse,c,lab,it=min([kmeans(X,3,s) for s in range(8)],key=lambda v:v[0])
print("k3 SSE",round(sse,3),"iteraties",it)
print("clusteromvang",np.bincount(lab).tolist())
print("centra",np.round(c,3))
print("SSE per k",[(k,round(min(kmeans(X,k,s)[0] for s in range(8)),2)) for k in (1,2,3,4)])