restart:

FuzzyCMeans:= proc(
     X::Matrix(datatype= float[8]),   #The rows of X are the data points.
     #number of clusters and number of means:
     {clusters::And(posint, satisfies(x-> x > 1)):= 2},
     {fuzziness::satisfies(x-> x > 1):= 2},   #fuzziness factor
     {epsilon::positive:= 1e-4},   #convergence criterion
     {max_iters::posint:= 999},   #max number of iterations
     {Norm::procedure:= proc(V) local x; sqrt(add(x^2, x= V)) end proc}
)
option `Written by Carl J Love, 2016-Apr-13`;
description
    "Fuzzy C-Means clustering algorithm.
     See http://home.deib.polimi.it/matteucc/Clustering/tutorial_html/cmeans.html
     and https://en.wikipedia.org/wiki/Fuzzy_clustering
    "
;
uses LA= LinearAlgebra;
local
     c:= clusters,     #number of clusters and number of means
     m:= fuzziness,    #fuzziness factor
     n:= op([1,1], X), #number of points
     d:= op([1,2], X), #dimension
     N:= Norm,         #d-dimensional vector norm
     #U[i,j] will be the degree (on a 0..1 scale) to which point i belongs to cluster j.
     U:= Matrix((n,c), datatype= float[8]),
     UM:= Matrix((n,c), datatype= float[8]), # U^~m
     #U1 is just an update copy of U.
     U1:= LA:-RandomMatrix((n,c), generator= 0..1, datatype= float[8]),
     C:= Matrix((c,d), datatype= float[8]), #The rows of C will be the cluster centers.
     e:= 2/(m-1), #fuzzy exponent used to update U
     err:= epsilon,  #max matrix difference on current iteration
     i, j, k, jj #loop indices
;
     for jj to max_iters while err >= epsilon do
          ArrayTools:-Copy(U1, U);
          ArrayTools:-Copy(U, UM);
          map[evalhf, inplace](`^`, UM, m);
          for j to c do
               C[j,..]:= add(UM[i,j]*X[i,..], i= 1..n) / add(UM[i,j], i= 1..n)
          end do;
          for i to n do  for j to c do
               U1[i,j]:= 1/add((N(X[i,..]-C[j,..])/N(X[i,..]-C[k,..]))^e, k= 1..c)
          end do  end do;
          LA:-MatrixAdd(U,U1,1,-1,inplace); #U:= U-U1
          err:= max(map[evalhf, inplace](abs, U))
     end do;
     if jj > max_iters then  
          error "Didn't converge; error = %1. Increase max_iters or epsilon.", err
     end if;
     userinfo(1, FuzzyCMeans, "Used", jj-1, "iterations.");      
     (C,U1)                
end proc:

AssignToClusters:= (U::Matrix)->
     map(proc(i) local p; member(max(U[i,..]), U[i,..], p); p end, [$1..op([1,1], U)])
:
     

PlotClusters2D:= proc(X::Matrix, U::Matrix)
local A:= AssignToClusters(U), k, J:= [$1..op([1,1], U)], c:= op([1,2], U);
     plot(
          [seq(X[select(i-> A[i]=k, J), ..], k= 1..c)],
          color= [seq(COLOR(HUE, k/c), k= 1..c)],
          style= point, symbol= diagonalcross, symbolsize= 16, _rest
     )
end proc:

PlotClusterMeans2D:= proc(C::Matrix)
local c:= op([1,1], C), k;
     plot(
          `[]`~(convert(C, listlist)),
          color= [seq(COLOR(HUE, k/c), k= 1..c)],
          style= point, symbol= circle, symbolsize= 32, _rest
     )
end proc:

PlotClustersAndMeans2D:= (X::Matrix, C::Matrix, U::Matrix)->
     plots:-display(
          [PlotClusters2D(X,U), PlotClusterMeans2D(C)],
          gridlines= false, _rest
     )
:
 

#Test data:
X:=
     [[5.37,19.50],[5.73,19.44],[4.66,20.03],[5.66,20.38],[4.22,21.97],
      [5.08,20.67],[5.08,19.08],[4.54,20.06],[4.35,19.82],[5.19,20.72],
      [4.48,19.95],[5.76,19.81],[4.15,18.68],[6.37,20.60],[5.58,21.13],
      [5.76,19.91],[5.85,19.02],[5.00,19.71],[5.42,20.31],[4.77,21.45],
      [8.61,19.48],[10.70,20.31],[10.99,20.28],[11.68,21.28],[9.12,20.77],
      [10.30,20.07],[10.40,21.62],[10.95,20.34],[9.79,20.29],[9.69,20.86],
      [10.02,21.45],[11.05,20.19],[9.20,17.96],[10.49,19.88],[9.61,19.49],
      [10.33,19.59],[9.29,20.94],[10.17,19.64],[10.97,20.32],[10.08,19.16]
     ]
:
XM:= Matrix(X, datatype= float[8]):

infolevel[FuzzyCMeans]:= 1:

(C,U):= FuzzyCMeans(
     XM,
     #All options are set to their default values here:
     clusters= 2, fuzziness= 2, epsilon= 1e-4, max_iters= 999,
     Norm= proc(V) local x; sqrt(add(x^2, x= V)) end proc
);

FuzzyCMeans: Used 8 iterations.

PlotClustersAndMeans2D(XM, C, U);

(C,U):= CodeTools:-Usage(FuzzyCMeans(XM, clusters= 3));

FuzzyCMeans: Used 60 iterations.

memory used=60.99MiB, alloc change=0 bytes, cpu time=530.00ms, real time=534.00ms

PlotClustersAndMeans2D(XM, C, U);

(C,U):= CodeTools:-Usage(FuzzyCMeans(XM, clusters= 3));

FuzzyCMeans: Used 47 iterations.
memory used=47.79MiB, alloc change=0 bytes, cpu time=483.00ms, real time=481.00ms

PlotClustersAndMeans2D(XM, C, U);