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.

 
(import (scheme base) ; headers for R7RS Scheme
        (scheme file)
        (scheme inexact)
        (scheme write)
        (srfi 1 lists)
        (srfi 27 random-bits))
 
;; calculate euclidean distance between points, any dimension
(define (euclidean-distance pt1 pt2)
  (sqrt (apply + (map (lambda (x y) (square (- x y))) pt1 pt2))))
 
;; input
;; - K: the target number of clusters K
;; - data: a list of points in the Cartesian plane
;; output
;; - a list of K centres
(define (kmeans++ K data)
  (define (select-uniformly data)
    (let loop ((index (random-integer (length data))) ; uniform selection of index
               (rem data)
               (front '()))
      (if (zero? index)
        (values (car rem) (append (reverse front) (cdr rem)))
        (loop (- index 1) (cdr rem) (cons (car rem) front)))))
  ;
  (define (select-weighted centres data)
    (define (distance-to-nearest datum)
      (apply min (map (lambda (c) (euclidean-distance c datum)) centres)))
    ;
    (let* ((weights (map (lambda (d) (square (distance-to-nearest d))) data))
           (target-weight (* (apply + weights) (random-real))))
      (let loop ((rem data)
                 (front '())
                 (weight-sum 0.0)
                 (wgts weights))
        (if (or (>= weight-sum target-weight) (null? (cdr rem)))
          (values (car rem) (append (reverse front) (cdr rem)))
          (loop (cdr rem)
                (cons (car rem) front)
                (+ weight-sum (car wgts))
                (cdr weights))))))
  ;
  (let-values (((pt rem) (select-uniformly data)))
              (let loop ((centres (list pt))
                         (items rem))
                (if (= (length centres) K)
                  centres
                  (let-values (((pt rem) (select-weighted centres items)))
                              (loop (cons pt centres)
                                    rem))))))
 
;; assign a point into a cluster
;; input: a point and a list of cluster centres
;; output: index of cluster centre
(define (assign-cluster pt centres)
  (let* ((distances (map (lambda (centre) (euclidean-distance centre pt)) centres))
         (smallest (apply min distances)))
    (list-index (lambda (d) (= d smallest)) distances)))
 
;; input
;; - num: the number of clusters K
;; - data: a list of points in the Cartesian plane
;; output
;; - list of K centres
(define (cluster K data)
  (define (centroid-for-cluster i assignments)
    (let* ((cluster (map cadr (filter (lambda (a-d) (= (car a-d) i)) (zip assignments data))))
           (length-cluster (length cluster)))
      ; compute centroid for cluster
      (map (lambda (vals) (/ (apply + vals) length-cluster)) (apply zip cluster))))
  ;
  (define (update-centres assignments)
    (map (lambda (i) (centroid-for-cluster i assignments)) (iota K)))
  ;
  (let ((initial-centres (kmeans++ K data)))
    (let loop ((centres initial-centres)
               (assignments (map (lambda (datum) (assign-cluster datum initial-centres)) data)))
      (let* ((new-centres (update-centres assignments))
             (new-assignments (map (lambda (datum) (assign-cluster datum new-centres)) data)))
        (if (equal? assignments new-assignments)
          new-centres
          (loop new-centres new-assignments))))))
 
;; using eps output, based on that in C - only works for 2D points
(define (save-as-eps filename data clusters K)
  (when (file-exists? filename) (delete-file filename))
  (with-output-to-file
    filename
    (lambda ()
      (let* ((W 400)
             (H 400)
             (colours (make-vector (* 3 K) 0.0))
             (max-x (apply max (map car data)))
             (min-x (apply min (map car data)))
             (max-y (apply max (map cadr data)))
             (min-y (apply min (map cadr data)))
             (scale (min (/ W (- max-x min-x))
                         (/ H (- max-y min-y))))
             (cx (/ (+ max-x min-x) 2))
             (cy (/ (+ max-y min-y) 2)))
 
        ;; set up colours
        (for-each 
          (lambda (i) 
            (vector-set! colours (+ (* i 3) 0) (inexact (/ (modulo (* 3 (+ i 1)) 11) 11)))
            (vector-set! colours (+ (* i 3) 1) (inexact (/ (modulo (* 7 i) 11) 11)))
            (vector-set! colours (+ (* i 3) 2) (inexact (/ (modulo (* 9 i) 11) 11))))
          (iota K))
 
        (display ;; display header
          (string-append 
            "%!PS-Adobe-3.0\n%%BoundingBox: -5 -5 " 
            (number->string (+ 10 W)) " " (number->string (+ 10 H)) "\n"
            "/l {rlineto} def /m {rmoveto} def\n"
            "/c { .25 sub exch .25 sub exch .5 0 360 arc fill } def\n"
            "/s { moveto -2 0 m 2 2 l 2 -2 l -2 -2 l closepath "
            "	gsave 1 setgray fill grestore gsave 3 setlinewidth"
            " 1 setgray stroke grestore 0 setgray stroke }def\n"))
 
        ;; display points
        (for-each ; top loop runs over the clusters
          (lambda (i)
            (display
              (string-append (number->string (vector-ref colours (* i 3)))
                             " "
                             (number->string (vector-ref colours (+ (* i 3) 1)))
                             " "
                             (number->string (vector-ref colours (+ (* i 3) 2)))
                             " setrgbcolor\n"))
            (for-each ;loop over points in cluster
              (lambda (pt)
                (when (= i (assign-cluster pt clusters))
                  (display
                    (string-append (number->string (+ (* (- (car pt) cx) scale) (/ W 2)))
                                   " "
                                   (number->string (+ (* (- (cadr pt) cy) scale) (/ H 2)))
                                   " c\n"))))
              data)
            (let ((center (list-ref clusters i))) ; display cluster centre
              (display 
                (string-append "\n0 setgray "
                               (number->string (+ (* (- (car center) cx) scale) (/ W 2)))
                               " "
                               (number->string (+ (* (- (cadr center) cy) scale) (/ H 2)))
                               " s\n"))))
          (iota K))
        (display "\n%%EOF")))))
 
;; extra credit 1: creates a list of n random points in n-D unit square
(define (make-data num-points num-dimensions)
  (random-source-randomize! default-random-source)
  (map (lambda (i) (list-tabulate num-dimensions (lambda (i) (random-real)))) (iota num-points)))
 
;; extra credit 2, uses eps visualisation to display result
(define (tester-1 num-points K)
  (let ((data (make-data num-points 2)))
    (save-as-eps "clusters-1.eps" data (cluster K data) K)))
 
;; extra credit 3: uses radians instead to make data
(define (tester-2 num-points K radius)
  (random-source-randomize! default-random-source)
  (let ((data (map (lambda (i) 
                     (let ((ang (* (random-real) 2 (* 4 (atan 1))))
                           (rad (* radius (random-real))))
                       (list (* rad (cos ang)) (* rad (sin ang)))))
                   (iota num-points))))
    ;; extra credit 2, uses eps visualisation to display result
    (save-as-eps "clusters-2.eps" data (cluster K data) K)))
 
;; extra credit 4: arbitrary dimensions - already handled, as all points are lists
(define (tester-3 num-points K num-dimensions)
  (display "Results:\n")
  (display (cluster K (make-data num-points num-dimensions))) 
  (newline))
 
(tester-1 30000 6)
(tester-2 30000 6 10)
(tester-3 30000 6 5)
 