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.

 
***********************************************************************
* KMPP - K-Means++ - Traditional data clustering with a special initialization
* Public Domain - This program may be used by any person for any purpose.
*
* Origin:
*    Hugo Steinhaus, 1956
*
* Refer to:
*    "kmeans++: the advantages of careful seeding"
*    David Arthur and Sergei Vassilvitskii
*    Proceedings of the eighteenth annual ACM-SIAM symposium 
*      on Discrete algorithms, 2007
*
*____Variable_______I/O_______Description___________________Type_______
*    X(P,N)         In        Data points                   Real
*    P              In        Dimension of the data         Integer
*    N              In        Number of points              Integer
*    K              In        # clusters                    Integer
*    C(P,K)         Out       Center points of clusters     Real
*    Z(N)           Out       What cluster a point is in    Integer
*    WORK(N)        Neither                                 Real
*    IFAULT         Out       Error code                    Integer
************************************************************************
      SUBROUTINE KMPP (X, P, N, K, C, Z, WORK, IFAULT)
 
       IMPLICIT NONE
       INTEGER P, N, K, Z, IFAULT
       REAL X, C, WORK
       DIMENSION X(P,N), C(P,K), Z(N), WORK(N)
 
*               constants
       INTEGER ITER                 ! maximum iterations
       REAL BIG                     ! arbitrary large number
       PARAMETER (ITER = 1000,
     $            BIG = 1E33)
 
*                local variables
       INTEGER 
     $         H,          ! count iterations
     $         I,          ! count points
     $         I1,         ! point marked as initial center
     $         J,          ! count dimensions
     $         L,          ! count clusters
     $         L0,         ! present cluster ID
     $         L1          ! new cluster ID
 
       REAL 
     $      BEST,                 ! shortest distance to a center
     $      D2,                   ! squared distance
     $      TOT,                  ! a total
     $      W                     ! temp scalar
 
       LOGICAL CHANGE             ! whether any points have been reassigned
 
************************************************************************
*           Begin.
************************************************************************
       IFAULT = 0
       IF (K < 1 .OR. K > N) THEN       ! K out of bounds
         IFAULT = 3
         RETURN
       END IF
       DO I = 1, N                       ! clear Z
         Z(I) = 0
       END DO
 
************************************************************************
*        initial centers
************************************************************************
       DO I = 1, N
         WORK(I) = BIG
       END DO
 
       CALL RANDOM_NUMBER (W)
       I1 = MIN(INT(W * FLOAT(N)) + 1, N)  ! choose first center at random
       DO J = 1, P
         C(J,1) = X(J,I1)
       END DO
 
       DO L = 2, K                    ! initialize other centers
         TOT = 0.
         DO I = 1, N                     ! measure from each point
           BEST = WORK(I)
           D2 = 0.                         ! to prior center
           DO J = 1, P
             D2 = D2 + (X(J,I) - C(J,L-1)) **2  ! Squared Euclidean distance
             IF (D2 .GE. BEST) GO TO 10               ! needless to add to D2
           END DO                          ! next J
           IF (D2 < BEST) BEST = D2          ! shortest squared distance 
           WORK(I) = BEST 
  10       TOT = TOT + BEST             ! cumulative squared distance
         END DO                      ! next data point
 
************************************************************************
* Choose center with probability proportional to its squared distance
*     from existing centers.
************************************************************************
         CALL RANDOM_NUMBER (W)
         W = W * TOT    ! uniform at random over cumulative distance
         TOT = 0.
         DO I = 1, N
           I1 = I
           TOT = TOT + WORK(I)
           IF (TOT > W) GO TO 20
         END DO                ! next I
  20     CONTINUE
         DO J = 1, P         ! assign center
           C(J,L) = X(J,I1)
         END DO
       END DO               ! next center to initialize
 
************************************************************************
*                      main loop
************************************************************************
       DO H = 1, ITER
         CHANGE = .FALSE.
 
*             find nearest center for each point 
         DO I = 1, N
           L0 = Z(I)
           L1 = 0
           BEST = BIG
           DO L = 1, K
             D2 = 0.
             DO J = 1, P
               D2 = D2 + (X(J,I) - C(J,L)) **2
               IF (D2 .GE. BEST) GO TO 30
             END DO
  30         CONTINUE
             IF (D2 < BEST) THEN           ! new nearest center
               BEST = D2
               L1 = L
             END IF             
           END DO        ! next L
 
           IF (L0 .NE. L1) THEN
             Z(I) = L1                   !  reassign point 
             CHANGE = .TRUE.
           END IF
         END DO         ! next I
         IF (.NOT. CHANGE) RETURN      ! success
 
************************************************************************
*           find cluster centers
************************************************************************
         DO L = 1, K              ! zero population
           WORK(L) = 0.
         END DO
         DO L = 1, K               ! zero centers
           DO J = 1, P
             C(J,L) = 0.
           END DO
         END DO
 
         DO I = 1, N
           L = Z(I)
           WORK(L) = WORK(L) + 1.             ! count
           DO J = 1, P
             C(J,L) = C(J,L) + X(J,I)         ! add
           END DO
         END DO
 
         DO L = 1, K
           IF (WORK(L) < 0.5) THEN          ! empty cluster check
             IFAULT = 1                     ! fatal error
             RETURN
           END IF
           W = 1. / WORK(L)
           DO J = 1, P
             C(J,L) = C(J,L) * W     ! multiplication is faster than division
           END DO
         END DO
 
       END DO                   ! next H
       IFAULT = 2                ! too many iterations
       RETURN
 
      END  ! of KMPP
 
 
************************************************************************
*             test program (extra credit #1)
************************************************************************
      PROGRAM TPEC1
       IMPLICIT NONE
       INTEGER N, P, K
       REAL TWOPI
       PARAMETER (N = 30 000,
     $            P = 2,
     $            K = 6,
     $            TWOPI = 6.2831853)
       INTEGER I, L, Z(N), IFAULT
       REAL X(P,N), C(P,K), R, THETA, W, WORK(N)
 
*             Begin
       CALL RANDOM_SEED()
       DO I = 1, N                      ! random points over unit circle
         CALL RANDOM_NUMBER (W)
         R = SQRT(W)                      ! radius
         CALL RANDOM_NUMBER (W)
         THETA = W * TWOPI                ! angle
         X(1,I) = R * COS(THETA)          ! Cartesian coordinates
         X(2,I) = R * SIN(THETA)
       END DO
 
*             Call subroutine
       CALL KMPP (X, P, N, K, C, Z, WORK, IFAULT)
       PRINT *, 'kmpp returns with error code ', IFAULT
 
*            Print lists of points in each cluster
       DO L = 1, K
         PRINT *, 'Cluster ', L, ' contains points: '
  10     FORMAT (I6, $)
  20     FORMAT ()
         DO I = 1, N
           IF (Z(I) .EQ. L) PRINT 10, I
         END DO
         PRINT 20
       END DO
 
*         Write CSV file with Y-coordinates in different columns by cluster
       OPEN (UNIT=1, FILE='tpec1.csv', STATUS='NEW', IOSTAT=IFAULT)
       IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble opening file'
  30   FORMAT (F8.4, $)
  40   FORMAT (',', $)
  50   FORMAT (F8.4)
       DO I = 1, N
         WRITE (UNIT=1, FMT=30, IOSTAT=IFAULT) X(1,I)
         IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble writing X-coord'
         DO L = 1, Z(I)                     ! one comma per cluster ID
           WRITE (UNIT=1, FMT=40, IOSTAT=IFAULT)
           IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble writing comma'
         END DO
         WRITE (UNIT=1, FMT=50, IOSTAT=IFAULT) X(2,I)
         IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble writing Y-coord'
       END DO
 
*           Write the centroids in the far column
       DO L = 1, K
         WRITE (UNIT=1, FMT=30, IOSTAT=IFAULT) C(1,L)
         IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble writing X-coord'
         DO I = 1, K+1
           WRITE (UNIT=1, FMT=40, IOSTAT=IFAULT)
           IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble writing comma'
         END DO
         WRITE (UNIT=1, FMT=50, IOSTAT=IFAULT) C(2,L)
         IF (IFAULT .NE. 0) PRINT *, 'tpec1: trouble writing Y-coord'
       END DO
       CLOSE (UNIT=1)
 
      END  ! of test program
 