K-means++ clustering
K-means
 This data was partitioned into 7 clusters using the K-means algorithm.
The task is to implement the K-means++ algorithm. Produce a function which takes two arguments: the number of clusters K, and the dataset to classify. K is a positive integer and the dataset is a list of points in the Cartesian plane. The output is a list of clusters (related sets of points, according to the algorithm).

For extra credit (in order):

 Provide a function to exercise your code, which generates a list of random points.
 Provide a visualization of your results, including centroids (see example image). 
 Generalize the function to polar coordinates (in radians).
 Generalize the function to points in an arbitrary N space (i.e. instead of x,y pairs, points are an N-tuples representing coordinates in ℝN). 
If this is different or more difficult than the [naive] solution for ℝ2, discuss what had to change to support N dimensions.
Extra credit is only awarded if the examples given demonstrate the feature in question. To earn credit for 1. and 2., visualize 6 clusters of 30,000 points in ℝ2. It is not necessary to provide visualization for spaces higher than ℝ2 but the examples should demonstrate features 3. and 4. if the solution has them.

include c:\cxpl\codes;  \intrinsic 'code' declarations
 
def     N = 30000;              \number of points
def     K = 6;                  \number of clusters
int     Px(N), Py(N), Pc(N),    \coordinates of points and their cluster
        Cx(K), Cy(K);           \coordinates of centroid of cluster
 
 
func Centroid;  \Find new centroids of points grouped with current centroids
int  Change, Cx0(K), Cy0(K), C, Count, I;
[Change:= false;
for C:= 0 to K-1 do                       \for each centroid...
        [Cx0(C):= Cx(C);  Cy0(C):= Cy(C); \save current centroid
        Cx(C):= 0;  Cx(C):= 0;  Count:= 0;\find new centroid
        for I:= 0 to N-1 do               \for all points
            if Pc(I) = C then             \ grouped with current centroid...
                [Cx(C):= Cx(C) + Px(I);
                 Cy(C):= Cy(C) + Py(I);
                 Count:= Count+1;
                ];
        Cx(C):= Cx(C)/Count;  Cy(C):= Cy(C)/Count;
        if Cx(C)#Cx0(C) or Cy(C)#Cy0(C) then Change:= true;
        ];
return Change;
];
 
 
proc Voronoi;                   \Group points with their nearest centroid
int  D2, MinD2, I, C;           \distance squared, minimum distance squared
[for I:= 0 to N-1 do            \for each point...
        [MinD2:= -1>>1;         \find closest centroid
        for C:= 0 to K-1 do
                [D2:= sq(Px(I)-Cx(C)) + sq(Py(I)-Cy(C));
                if D2 < MinD2 then
                        [MinD2:= D2;  Pc(I):= C];  \update closest centroid
                ];
        ];
];
 
 
proc KMeans;                    \Group points into K clusters
int  Change, I;
repeat  Voronoi;
        Change:= Centroid;
        SetVid($101);           \show result on 640x480x8 screen
        for I:= 0 to N-1 do Point(Px(I), Py(I), Pc(I)+1);
        for I:= 0 to K-1 do Point(Cx(I), Cy(I), \bright white\ $F);
until   Change = false;
 
 
proc Random(X, Y);              \Return random X,Y biased for polar coordinates
int  X, Y;
real A, D;
[D:= float(Ran(240));                   \distance: 0..239
A:= float(Ran(314159*2)) / 10000.0;     \angle:    0..2pi
X(0):= fix(D*Cos(A)) + 320;             \rectangular coords centered on screen
Y(0):= fix(D*Sin(A)) + 240;
];
 
 
int  I;
[for I:= 0 to N-1 do Random(@Px(I), @Py(I));    \random set of points
 for I:= 0 to K-1 do Random(@Cx(I), @Cy(I));    \random set of cluster centroids
KMeans;
I:= ChIn(1);                    \wait for keystroke
SetVid($03);                    \restore normal text screen
]