The initial parameters are n_clusters, n_features, n_cases
 that describe the shape of the problem and n_outer, n_inner
 bounding the loops.

getint proc (lowerbound in int) ⇒ (result out int)
   get result
   assert result ≥ lowerbound
end proc

nclusters const = getint(2)
nfeatures const = getint(1)
ncases    const = getint(nclusters)
nouter    const = getint(1)
ninner    const = getint(1)

 Next we read the feature matrix.

feature var := {{0.0} @@ nfeatures} @@ ncases
loop for caseindex in [0,ncases)
    loop for featureindex in [0,nfeatures)
       get feature[caseindex][featureindex]
    end loop
end loop

 Make a random selection of cases to serve as cluster centres.

randomassignment proc (centre out [][]real)
   perm var []int
   loop for i in [0,ncases)
      perm[i] := i
   end loop
   loop for c in [0,nclusters)
      n const = ncases - c                 number of elements to choose from
      i const = n - 1                       last of those elements
      j const = trunc(float(n) * random())  j in [0,n)
      (perm[j], perm[i]) := (perm[i], perm[j])
      centre[c] := feature[perm[i]]
   end loop
end proc    

 The squared Euclidean distance between two feature vectors.

distancesquared proc (x in []real, y in []real) ⇒ (d out real)
   d := 0.0
   loop for i in [0,size x)
      d +:= (x[i] - y[i])**2
   end loop
end proc

repartition proc (
   oldcentre in  [][]real,
   newcentre out [][]real,
   group      out []int       case number -> group number
) ⇒ (
   sumofsquareddistances out real
)
   sumofsquareddistances := 0.0
   tally var := {0.0} @@ nclusters
   group := {}
   newcentre := {{0.0} @@ nfeatures} @@ nclusters
   loop for caseindex in [0,ncases)
      thiscase const = feature[caseindex]
      bestdistancesquared var := 0.0
      bestcluster          var := -1
      loop for c in [0,nclusters)
         d const = distancesquared(thiscase, oldcentre[c])
         if bestcluster < 0 or d < bestdistancesquared then
            (bestcluster, bestdistancesquared) := (c, d)
         end if
      end if
      assert bestcluster ≥ 0
      group[caseindex,) := {bestcluster}
      tally[bestcluster] +:= 1.0
      loop for featureindex in [0,nfeatures)
         newcentre[bestgroup, featureindex] +:= thiscase[featureindex]
      end loop
      loop for cluster in [0,nclusters)
         if tally[cluster] > 0.0 then
            loop for featureindex in [0,nfeatures)
               newcentre[cluster][featureindex] /= tally[cluster]
            end loop
         end if
      end loop
      sumofsquareddistances +:= bestdistancesquared
   end loop
end proc   

outerbestsumofsquareddistances var real
outerbestgroup var []int

loop for outer in [0,nouter)
   centreone var := {{0.0} @@ nfeatures} @@ nclusters
   centretwo var := {{0.0} @@ nfeatures} @@ nclusters
   d1         var := 0.0
   d2         var := 0.0
   group      var []int

   randomassignment(centreone)
   loop for inner in [0,ninner)
      d1 := repartition(centreone, centretwo, group)
      d2 := repartition(centretwo, centreone, group)
      while d2 < d1
   end loop
   if outer = 0 or d2 < outerbestsumofsquareddistances then
      (outerbestgroup, outerbestsumofsquareddistances) := (group, d2)
   end if
end loop

print outerbestgroup