``` 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
```