Golf the K-means algorithm K-means is a standard unsupervised clustering algorithm, which, given a set of "points" and a number of clusters K, will assign each "point" to one of K clusters.

Pseudo-Code of K-means

Note that there are many variants of K-means. You have to implement the algorithm I am describing below. You may have some variation on the algorithm or use built-ins as long as you would get the same result as this algorithm given the same initial points. 

In this challenge, all inputs will be points on the 2D plane (each point is represented by its coordinates in x and y).

Inputs: K, the number of clusters
        P, the set of points

Choose K points of P uniformly at random
Each chosen point is the initial centroid of its cluster

Loop:
     For each point in P:
         Assign to the cluster whose centroid is the nearest (Euclidean distance)
         In case of a tie, any of the tied cluster can be chosen

     Recompute the centroid of each cluster:
         Its x coordinate is the average of all x's of the points in the cluster
         Its y coordinate is the average of all y's of the points in the cluster

Until the clusters don't change from one iteration to the next

Output: the set of clusters    


Inputs and Outputs


You may take K and P through STDIN, or as a function argument, etc.
P and the points in P can be represented using any structure that is natural for set/lists in your language of choice.
K is a strictly positive integer.
You may assume that inputs are valid.
There will always be at least K points in P.
You may output the clusters to STDOUT, return them from a function, etc.
The ordering of the clusters and the ordering inside the clusters is unimportant.
-You can either return groups of points to represent clusters, or each point labeled with an identifier for the cluster (e.g. an integer).


Test cases

Since the resulting clusters depend on which points were initially chosen, you may not all get the same results (or the same result each time you run your code).

Therefore, only take the output as an example output.

Input:
  K = 1
  P = [[1,2.5]]
Output:
  [[[1,2.5]]]

Input:
  K = 3
  P = [[4,8], [15,16], [23,42], [-13.37,-12.1], [666,-666]]
Output:
  [[[666,-666]],[[-13.37,-12.1],[4,8]],[[15,16],[23,42]]]

Input:
  K = 2
  P = [[1,1], [1,1], [1,1]]
Output:
  [[[1,1]],[[1,1],[1,1]]]


Scoring

This is code-golf, so the shortest answer in bytes wins.
C++, 479 474 bytes



Only ~20x as much as Matlab!

Golfed

#define V vector&lt;P&gt;
#define f float
struct P{f x,y,i=0;f d(P&amp;p){return(p.x-x)*(p.x-x)+(p.y-y)*(p.y-y);}f n(P&amp;p){return i?x/=i,y/=i,d(p):x=p.x,y=p.y,0;}f a(P&amp;p){x+=p.x,y+=p.y,i++;}};P z;int l(P a,P b){return a.d(z)&lt;b.d(z);}f m(f k,V&amp;p){f s=p.size(),i,j=0,t=1;V c(k),n=c,d;for(random_shuffle(p.begin(),p.end());j&lt;k;c[j].i=j++)c[j]=p[j];for(;t;c=n,n=V(k)){for(i=0;i&lt;s;i++)d=c,z=p[i],sort(d.begin(),d.end(),l),j=d[0].i,p[i].i=j,n[j].a(p[i]);for(j=t=0;j&lt;k;j++)t+=n[j].n(c[j]);}}


The input/output to the algorithm is a set of points (struct P) with x and y; and the output is the same set, with their i's set to indicate the index of the output cluster that the point ends in. 

That extra i is also used to identify the clusters. In the main loop, the closest centroid to each point is found by sorting a copy of the current centroids by proximity to that point.

This handles degenerate cases (empty clusters) by keeping the corresponding centroids' previous position (see definition of P::n, which also returns distance-to-previous-centroid). A few chars could be saved by assuming that these will not crop up.

Ungolfed, with main

#include &lt;cstdio&gt;
#include &lt;ctime&gt;
#include &lt;cstdlib&gt;
#include &lt;vector&gt;
#include &lt;algorithm&gt;
using namespace std;

#define V vector&lt;P&gt;
#define f float
struct P{
    f x,y,i=0;
    f d(P&amp;p){return(p.x-x)*(p.x-x)+(p.y-y)*(p.y-y);} // distance squared
    f n(P&amp;p){return i?x/=i,y/=i,d(p):x=p.x,y=p.y,0;} // normalize-or-reset
    f a(P&amp;p){x+=p.x,y+=p.y,i++;}                     // add coordinates
};
P z;int l(P a,P b){return a.d(z)&lt;b.d(z);}            // closer-to-z comparator 
f m(f k,V&amp;p){
    f s=p.size(),i,j=0,t=1;V c(k),n=c,d;
    for(random_shuffle(p.begin(),p.end());j&lt;k;c[j].i=j++)
        c[j]=p[j];                                // initial random assignment
    for(;t;c=n,n=V(k)){                           
        for(i=0;i&lt;s;i++)                          // assign to clusters
            d=c,z=p[i],sort(d.begin(),d.end(),l),
            j=d[0].i,p[i].i=j,n[j].a(p[i]);       // and add those coords
        for(j=t=0;j&lt;k;j++)t+=n[j].n(c[j]);        // normalize &amp; count changes
    }        
}

int main(int argc, char **argv) {
    srand((unsigned long)time(0));

    int k;
    V p;
    sscanf(argv[1], "%d", &amp;k);
    printf("Input:\n");
    for (int i=2,j=0; i&lt;argc; i+=2, j++) {
        P n;
        sscanf(argv[i], "%f", &amp;(n.x));
        sscanf(argv[i+1], "%f", &amp;(n.y));
        p.push_back(n);
        printf("%d : %f,%f\n", j, p[j].x, p[j].y);
    }

    m(k,p);
    printf("Clusters:\n");
    for (int q=0; q&lt;k; q++) {
        printf("%d\n", q);
        for (unsigned int i=0; i&lt;p.size(); i++) {
            if (p[i].i == q) printf("\t%f,%f (%d)\n", p[i].x, p[i].y, i);
        }
    }
    return 0;
}