restart:

line:= proc(Pt1::[realcons,realcons], Pt2::[realcons,realcons])
global Lines;
     Lines[Pt1,Pt2]:= ()
end proc:

Point:= proc(x::realcons, y::realcons)
global Points;
     Points[x,y]:= ()
end proc:

Rand:= RandomTools:-Generate(float(range= 0..1, method= uniform), makeproc):

# Firstly, set the size limit for the grid, which shall define how many

# points there are in both the vertical and horizontal direction

size:= 40;

# Iteration through x, using the line function to create a number of

# horizontal lines up to the limit point

for x to size do  line([0, size], [x, x])  end do:

# Now iterating through the diagonals, using two sets of

# lines, both starting from the line going through the origin. The first

# iterates through the lines that pass through positive x values at y=0,

# the second iterates through the lines passing through negative x values

# at y=0

# d must be increased by 2 each time to ensure that equilateral triangles

# are formed.

for d from 0 by 2 to size do
     #Upward diagonals

     line([d, size], [0, size-d]); line([0, size-d], [d, size]);

     #Downward diagonals:

     line([d, 0], [0, d]); line([size, d], [d, size])

end do:

# The program now iterates through potential y values for the potential

# points. Then it tests the value of y to see what orientation the lattice

# is at that point, then randomly selects x values to enter, iterating

# between all the possible x values in the orientation before moving on to

# the next y value.

for y from 0.5 by 0.5 while y<size do
     (i,s):= `if`(frem(y,1)=.5, [0,1], `if`(frem(y,2)=1, [0,2], [1,2]))[];
     seq(`if`(Rand() > 0.5, Point(x,y), [][]), x= i..size-1, s)

end do:

plots:-display(
     [plot([indices(Lines)], thickness= 0, color= gray),
      plot(
          [indices(Points)],
          style= point, symbol= soliddiamond, symbolsize= 1, color= red
      )], gridlines= false
);

restart:

FuzzyCMeans:= proc(
     X::Matrix,   #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:= (x-> LinearAlgebra:-Norm(x,2))}
)
option `Written by Carl J Love, 2016-Apr-11`;
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 clusrters and number of menas
     m:= fuzziness,    #fuzziness factor
     n:= op([1,1], X), #number of points
     d:= op([1,2], X), #dimension
     N:= Norm,
     #U[i,j] is the degree (on a 0..1 scale) to which point i belongs to cluster j.
     U:= Matrix((n,c), datatype= float[8]),
     #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 are the cluster centers.
     e:= 2/(m-1),
     i, j, k, jj, #loop indices
     UM # U^~m
;
     for jj to max_iters while LA:-Norm(U-U1, infinity) >= epsilon do
          U:= copy(U1);
          UM:= U^~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;
          U1:= Matrix(
               (n,c),
               (i,j)-> 1/add((N(X[i,..]-C[j,..])/N(X[i,..]-C[k,..]))^e, k= 1..c),
               datatype= float[8]
          )
     end do;
     if jj > max_iters then error "Didn't converge" end if;
     userinfo(1, FuzzyCMeans, "Used", jj, "iterations.");
     (C,U1)                
end proc:

#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= (x-> LinearAlgebra:-Norm(x,2))
);

FuzzyCMeans: Used 9 iterations.

#Plot. First separate points into clusters:
(CL1,CL2):= selectremove(i-> U[i,1] > U[i,2], [$1..nops(X)]):
plots:-display(
     [plot(
          map(CL-> XM[CL,..], [CL1,CL2]), color= [red, green],
          style= point, symbol= diagonalcross, symbolsize= 16
      ),
      plot(
          `[]`~(convert(C, listlist)), color= [red, green],
          style= point, symbol= circle, symbolsize= 32
      )
     ], gridlines= false
 );

J:= [0,4,1,2,3]:

G[0]:= add(a[k]*v[k], k= 1..4);

for j from 2 to nops(J) do
     k:= J[j];
     F[k]:= expand(t[k]*G[J[j-1]]);
     for ii to 4 do for jj to 4 do
          F[k]:= expand(algsubs(t[ii]*v[jj]= M[ii,jj], F[k]))
     od od;
     G[k]:= expand(subs(epsilon= E[k], F[k]))
od;