Concept · intermediate
k-means clustering
Group points into k clusters around moving centroids.
k-means groups points into clusters. Each cluster has a centroid, and the algorithm alternates between assigning points to their nearest centroid and moving each centroid to the mean of its points, until nothing changes. Unlike k-nearest neighbors, it has no labeled examples to learn from. You hand it raw points and a number , and it finds the groups on its own.
The dataset
data.csv holds 150 synthetic points, x and y, generated once by
make_data.py (fixed seed 20261101) and committed. It draws three round
blobs of 50 points, centered at (2, 2), (6, 3), and (3.5, 6), with a
standard deviation of 1.3. The centers are only 3.9 to 4.3 apart, so the
blobs overlap.
The first three rows come one from each blob, because k-means here starts from them. The generator also rejects any draw where a point is ever within of being equally close to two centroids, since rounding could then send it different ways in different languages. The first draw passed.
The method
With :
- Start with the first three points as the centroids.
- Assign every point to its nearest centroid by Euclidean distance.
- If no point changed cluster since the last assignment, stop.
- Update each centroid to the mean of its assigned points, then go back to step 2.
iterations counts the assignment steps, including the final one where
nothing changed. inertia is the sum of squared distances from each point to
its centroid:
Why does it stop? Both steps can only lower the inertia. Moving a point to its nearest centroid can’t increase its squared distance, and the mean is the position with the smallest total squared distance to a group. There are finitely many ways to split the points, so the falling inertia must run out.
On this data:
- Pass 1 assigns all 150 points: 45, 51, and 54 per cluster. The update moves the centroids from (1.25, 1.66), (6.42, 3.18), and (2.41, 5.11) to (1.555, 1.913), (6.203, 3.113), and (3.011, 5.454). Inertia falls from 467.03 to 431.41.
- Pass 2 moves 4 points, and no centroid then moves more than 0.16.
- Passes 3 and 4 each move 2 points from cluster 3 to cluster 1. Cluster 2 no longer changes, so its centroid stays put.
- Pass 5 changes nothing, and the loop stops.
The first step does most of the work. After that, only points near the borders switch sides, and the centroids drift by hundredths.
Reading the results
The final centroids are (1.700, 2.138), (6.205, 3.043), and (3.131, 5.750), holding 52, 50, and 48 points, with an inertia of 423.63. Each lies 0.21 to 0.45 from one of the true centers.
Two things cause the gap. Each blob is only a sample of 50 points, so its own mean isn’t exactly its center. And the blobs overlap: 15 of the 150 points end up in a cluster other than the blob they came from. k-means can’t know where a point came from. It only minimizes inertia, and by that measure it beats the truth: grouping points by their real blobs gives an inertia of 487.47.
The answer also depends on the starting centroids, and the algorithm never chooses for you. We fixed both, so every language has one result to match.
What to watch for across languages
None of the from-scratch implementations calls a ready-made k-means. R has
one, kmeans(), and the R tab’s Built-in switch shows it. Its defaults
aren’t this page’s method, though: it starts from random points and uses
the Hartigan-Wong algorithm, which moves single points between clusters.
Two runs with different seeds found inertias of 422.90 and 422.17, both a
little lower than our 423.63. Passing the first three points and
algorithm = "Lloyd" reproduces our result exactly.
The from-scratch versions use general tools: a distance formula (Python’s
math.dist, Go’s math.Hypot, Fortran’s norm2, written out by hand
elsewhere), sums, counts, and picking the smallest value. Watch how each
holds 2-D points: tuples, nested arrays, a matrix, two-element arrays in PHP and array references in Perl, two-element lists in Tcl, structs in Go, C, C++, and Swift, records in Java, Pascal, OCaml, and Ada, tables with x and y fields in Lua, a defstruct in Common Lisp, an X-Y term in Prolog,
fixed-size arrays in Rust, columns of a 2-D array in Fortran, separate x
and y arrays in awk, BASIC, and PostgreSQL, or a table of records in COBOL. Watch
how each notices that nothing changed: twenty-four compare the whole assignment
with the last one, while awk, BASIC, COBOL, C, and Pascal set a flag as they go. And
watch the four languages with no loops: Haskell’s, Prolog’s, and Elixir’s run
call themselves with the moved centroids, and PostgreSQL repeats the two steps with
WITH RECURSIVE. OCaml has loops, but its run calls itself the same way.
Data preview
| x | y |
|---|---|
| 1.25 | 1.66 |
| 6.42 | 3.18 |
| 2.41 | 5.11 |
| 1.22 | 3.21 |
| 4.95 | 3.83 |
| 2.06 | 1.96 |
| 3.32 | 4.29 |
| 7.44 | 6.26 |
First 8 of 150 rows · Download data.csv
Expected results
Every implementation must print these, to within a relative tolerance of 10⁻⁴.
iterations | 5 |
centroid_1_x | 1.700192 |
centroid_1_y | 2.138462 |
centroid_2_x | 6.204800 |
centroid_2_y | 3.043400 |
centroid_3_x | 3.130833 |
centroid_3_y | 5.749792 |
size_1 | 52 |
size_2 | 50 |
size_3 | 48 |
inertia | 423.627610 |
Implementations
"""k-means clustering (Lloyd's algorithm), using only the standard library.
Starts from the first K points as centroids, then alternates two steps until
no point changes cluster: assign each point to its nearest centroid, and move
each centroid to the mean of its points.
"""
import csv
import math
import sys
K = 3
MAX_ITERATIONS = 100
def nearest(point, centroids):
# min() returns the first of equal distances, so ties go to the lower index.
return min(range(len(centroids)), key=lambda j: math.dist(point, centroids[j]))
def main():
with open(sys.argv[1], newline="") as f:
points = [(float(row["x"]), float(row["y"])) for row in csv.DictReader(f)]
centroids = points[:K]
assignment = None
iterations = 0
while iterations < MAX_ITERATIONS:
iterations += 1
new_assignment = [nearest(p, centroids) for p in points]
if new_assignment == assignment:
break
assignment = new_assignment
for j in range(K):
members = [p for p, a in zip(points, assignment) if a == j]
if members: # an empty cluster keeps its old centroid
centroids[j] = (
sum(x for x, _ in members) / len(members),
sum(y for _, y in members) / len(members),
)
sizes = [assignment.count(j) for j in range(K)]
inertia = sum(math.dist(p, centroids[a]) ** 2 for p, a in zip(points, assignment))
print(f"iterations {iterations}")
for j, (x, y) in enumerate(centroids, start=1):
print(f"centroid_{j}_x {x:.6f}")
print(f"centroid_{j}_y {y:.6f}")
for j, size in enumerate(sizes, start=1):
print(f"size_{j} {size}")
print(f"inertia {inertia:.6f}")
if __name__ == "__main__":
main()Python keeps the whole loop short. math.dist gives the Euclidean distance,
and min(range(K), key=...) picks the nearest centroid; min returns the
first of equal values, so a tie would go to the lower index. Points are
plain tuples, and comparing two lists with == tells you when the
assignments stopped changing.
The one trap is an empty cluster, which has no mean. Our data never produces one, but the code keeps the old centroid if it does.
There was no built-in to skip this time: the standard library has no clustering at all, so in real code you’d reach for a third-party package.
-- k-means clustering (Lloyd's algorithm), using only Ada's standard
-- library.
--
-- Starts from the first K points as centroids, then alternates two steps
-- until no point changes cluster: assign each point to its nearest
-- centroid, and move each centroid to the mean of its points. Ada has no
-- k-means routine, so there is no built-in to avoid.
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO;
with Ada.Command_Line;
with Ada.Strings.Fixed; use Ada.Strings.Fixed;
with Ada.Containers.Vectors;
procedure Main is
K : constant := 3;
Max_Iterations : constant := 100;
-- A record is a value: assigning one copies both fields.
type Point is record
X, Y : Long_Float;
end record;
package Point_Vectors is new Ada.Containers.Vectors (Positive, Point);
subtype Cluster is Positive range 1 .. K;
type Point_Array is array (Positive range <>) of Point;
type Centroid_Array is array (Cluster) of Point;
type Assignment_Array is array (Positive range <>) of Natural;
-- The text of field number Col (counting from 1) in a CSV line.
function Field (Line : String; Col : Positive) return String is
Start : Positive := Line'First;
begin
for I in 1 .. Col - 1 loop
Start := Index (Line (Start .. Line'Last), ",") + 1;
end loop;
declare
Stop : constant Natural := Index (Line (Start .. Line'Last), ",");
begin
return (if Stop = 0 then Line (Start .. Line'Last)
else Line (Start .. Stop - 1));
end;
end Field;
function Read_Points (Path : String) return Point_Array is
File : File_Type;
Points : Point_Vectors.Vector;
X_Col, Y_Col : Natural := 0;
begin
Open (File, In_File, Path);
declare
Header : constant String := Get_Line (File);
begin
for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
if Field (Header, C) = "x" then
X_Col := C;
elsif Field (Header, C) = "y" then
Y_Col := C;
end if;
end loop;
end;
if X_Col = 0 or else Y_Col = 0 then
raise Constraint_Error with "no x or y column";
end if;
while not End_Of_File (File) loop
declare
Line : constant String := Get_Line (File);
begin
if Line /= "" then
-- Point'(...) names the aggregate's type; a bare one could
-- also be a vector, which makes Append ambiguous.
Points.Append (Point'(X => Long_Float'Value (Field (Line, X_Col)),
Y => Long_Float'Value (Field (Line, Y_Col))));
end if;
end;
end loop;
Close (File);
return Result : Point_Array (1 .. Natural (Points.Length)) do
for I in Result'Range loop
Result (I) := Points (I);
end loop;
end return;
end Read_Points;
-- Squared Euclidean distance. It ranks centroids in the same order as
-- the distance itself, so the nearest one needs no Sqrt.
function Squared_Distance (A, B : Point) return Long_Float is
((A.X - B.X) ** 2 + (A.Y - B.Y) ** 2);
-- The closest centroid. Only a strictly smaller distance replaces the
-- best so far, so a tie goes to the lower number.
function Nearest (P : Point; Centroids : Centroid_Array) return Cluster is
Best : Cluster := 1;
begin
for J in 2 .. K loop
if Squared_Distance (P, Centroids (J))
< Squared_Distance (P, Centroids (Best))
then
Best := J;
end if;
end loop;
return Best;
end Nearest;
procedure Say (Name : String; Value : Long_Float) is
begin
Put (Name & " ");
Put (Value, Fore => 1, Aft => 6, Exp => 0); -- six decimals, no E
New_Line;
end Say;
Points : constant Point_Array := Read_Points (Ada.Command_Line.Argument (1));
N : constant Natural := Points'Length;
Centroids : Centroid_Array;
-- 0 means "no cluster yet", so the first pass always counts as a change.
Assignment : Assignment_Array (1 .. N) := [others => 0];
Iterations : Natural := 0;
Inertia : Long_Float := 0.0;
Sizes : array (Cluster) of Natural;
begin
if N < K then
raise Constraint_Error with "need at least 3 points";
end if;
for J in Cluster loop
Centroids (J) := Points (J); -- a copy, so moving it leaves the point
end loop;
while Iterations < Max_Iterations loop
Iterations := Iterations + 1;
declare
New_Assignment : Assignment_Array (1 .. N);
Sums : array (Cluster) of Point := [others => (0.0, 0.0)];
begin
for I in Points'Range loop
New_Assignment (I) := Nearest (Points (I), Centroids);
end loop;
-- "=" on two arrays compares every element.
exit when New_Assignment = Assignment;
Assignment := New_Assignment;
-- Sum each cluster's points in one pass, then divide by its size.
Sizes := [others => 0];
for I in Points'Range loop
Sums (Assignment (I)).X := Sums (Assignment (I)).X + Points (I).X;
Sums (Assignment (I)).Y := Sums (Assignment (I)).Y + Points (I).Y;
Sizes (Assignment (I)) := Sizes (Assignment (I)) + 1;
end loop;
for J in Cluster loop
if Sizes (J) > 0 then -- an empty cluster keeps its old centroid
Centroids (J) := (Sums (J).X / Long_Float (Sizes (J)),
Sums (J).Y / Long_Float (Sizes (J)));
end if;
end loop;
end;
end loop;
Sizes := [others => 0];
for I in Points'Range loop
Sizes (Assignment (I)) := Sizes (Assignment (I)) + 1;
Inertia := Inertia + Squared_Distance (Points (I), Centroids (Assignment (I)));
end loop;
-- 'Image puts a space before a non-negative number; Trim removes it.
Put_Line ("iterations" & Natural'Image (Iterations));
for J in Cluster loop
declare
Num : constant String := Trim (Cluster'Image (J), Ada.Strings.Left);
begin
Say ("centroid_" & Num & "_x", Centroids (J).X);
Say ("centroid_" & Num & "_y", Centroids (J).Y);
end;
end loop;
for J in Cluster loop
Put_Line ("size_" & Trim (Cluster'Image (J), Ada.Strings.Left)
& Natural'Image (Sizes (J)));
end loop;
Say ("inertia", Inertia);
end Main;Each point is a Point record with X and Y fields, and the centroids
are a fixed array indexed by Cluster, a subtype covering 1 to 3. We start
from the first three points, assign each point to its nearest centroid, and
move each centroid to its cluster’s mean. Ada’s = compares two arrays
element by element, so one exit when New_Assignment = Assignment notices
that nothing changed.
Ada’s standard library has no k-means and no distance function, so we write
out the squared distance. One quirk: Append needs Point'(...), because
a bare aggregate could also be a vector, and the call becomes ambiguous.
⍝ k-means clustering (Lloyd's algorithm), from scratch in GNU APL.
⍝ Start from the first K points, then alternate two steps until no point
⍝ changes cluster: assign each point to its nearest centroid, and move each
⍝ centroid to the mean of its points. APL reads right to left, with no
⍝ operator precedence: 2×3+4 is 2×(3+4). ⍝ starts a comment; ← assigns.
∇NAME SAY N
⍝ Print "name value" for a whole number. ⍕ formats a number as text.
⎕←NAME,' ',⍕N
∇
∇NAME SAY6 X;T
⍝ Print "name value" with six decimal places.
T←0 6⍕|X ⍝ |X is the size of X; 0 6⍕ gives 6 decimals
T←(∨\T≠' ')/T ⍝ drop the leading blank ⍕ adds
T←((T[1]='.')/'0'),T ⍝ ⍕ writes .94 for 0.94, so put the 0 back
⎕←NAME,' ',((X<0)/'-'),T ⍝ sign as -, since APL's own minus is ¯
∇
K←3 ⋄ MAX_ITERATIONS←100
⍝ Read the file. ⎕ARG is the command line; ↑¯1↑ takes its last item, the
⍝ data path. ⎕FIO[26] reads the whole file as one string.
TEXT←⎕FIO[26]↑¯1↑⎕ARG
LINES←(TEXT≠⎕UCS 10)⊂TEXT ⍝ ⊂ splits at the 0s of the mask: one item per line
⍝ Split each (¨) line at commas; ⊃ stacks them into a table.
ROWS←⊃{(⍵≠',')⊂⍵}¨LINES
⍝ A column by header name (⍳), minus the header (1↓), as numbers (⍎¨).
⍝ , turns a one-letter name like 'x' into a string; ⍉ makes one row per point.
COLUMN←{⍎¨1↓ROWS[;ROWS[1;]⍳⊂,⍵]}
P←⍉⊃COLUMN¨'x' 'y'
∇A←NEAREST C;D
⍝ Distance from every point (rows) to every centroid (columns): ∘.- subtracts
⍝ each centroid's coordinate from each point's, making a table.
D←(((P[;1]∘.-C[;1])*2)+(P[;2]∘.-C[;2])*2)*0.5
⍝ Per row: where (⍳) is its minimum (⌊/)? Ties go to the first.
A←{⍵⍳⌊/⍵}¨⊂[2]D
∇
∇C←C MOVE A;B;N;E
B←A∘.=⍳K ⍝ table with a 1 where point i is in cluster j
N←+⌿B ⍝ cluster sizes: +⌿ sums down each column
⍝ Clusters with points; an empty one keeps its centroid.
E←(N>0)/⍳K
⍝ (⍉B)+.×P is a matrix product: it sums each cluster's points. Dividing by
⍝ the cluster's size, repeated for x and y (∘.×1 1), gives the mean.
C[E;]←((⍉B)+.×P)[E;]÷N[E]∘.×1 1
∇
∇KMEANS;NEW
C←K↑[1]P ⋄ A←0 ⋄ ITERATIONS←0 ⍝ the first K rows are the starting centroids
NEXT: ITERATIONS←ITERATIONS+1
NEW←NEAREST C
⍝ → jumps; /0 leaves the function once nothing changes.
→(NEW≡A)/0
A←NEW
C←C MOVE A
→(ITERATIONS<MAX_ITERATIONS)/NEXT
∇
KMEANS
'iterations' SAY ITERATIONS
⍝ SAY¨ prints one line per name and value. ∘., joins every prefix with x and y;
⍝ , flattens the table row by row, the same order as ,C.
(,({'centroid_',(⍕⍵),'_'}¨⍳K)∘.,'xy') SAY6¨ ,C
({'size_',⍕⍵}¨⍳K) SAY¨ +⌿A∘.=⍳K
⍝ C[A;] is each point's own centroid; , flattens.
'inertia' SAY6 +/,(P-C[A;])*2
)OFFAssignment is one table: ∘.- subtracts every centroid from every point, so
squaring, adding, and taking *0.5 gives a 150 by 3 table of distances.
For each row, ⍵⍳⌊/⍵ finds the first minimum, so a tie goes to the lower
centroid. To move the centroids, A∘.=⍳K marks each point’s cluster with a
1, and a matrix product, (⍉B)+.×P, sums each cluster’s points in one step.
GNU APL’s {} functions fit on one line, so the loop is a ∇ function that
jumps to a label with →. One trap: 0÷0 is 1 in APL, not an error, so we
divide only for clusters that have points.
# k-means clustering (Lloyd's algorithm), in POSIX awk.
#
# Starts from the first K points as centroids, then alternates two steps until
# no point changes cluster: assign each point to its nearest centroid, and move
# each centroid to the mean of its points. Points are stored as they are read;
# the iterations happen in END.
BEGIN { FS = ","; K = 3; MAX_ITERATIONS = 100 }
# The header line: find the x and y columns by name.
NR == 1 {
for (i = 1; i <= NF; i++) {
if ($i == "x") xcol = i
if ($i == "y") ycol = i
}
if (!xcol || !ycol) {
print "missing column x or y" > "/dev/stderr"
failed = 1
exit 1
}
next
}
{
n++
px[n] = $xcol + 0
py[n] = $ycol + 0
}
# Squared distance from point i to centroid j.
function dist2(i, j) {
return (px[i] - cx[j]) ^ 2 + (py[i] - cy[j]) ^ 2
}
# Index of the centroid nearest to point i. The strict "<" means a tie goes
# to the lower index.
function nearest(i, j, best) {
best = 1
for (j = 2; j <= K; j++)
if (dist2(i, j) < dist2(i, best)) best = j
return best
}
END {
if (failed) exit 1
for (j = 1; j <= K; j++) { cx[j] = px[j]; cy[j] = py[j] }
# assign[i] starts empty, so the first pass always counts as a change.
# The pass that changes nothing is counted too, then ends the loop.
iterations = 0
while (iterations < MAX_ITERATIONS) {
iterations++
changed = 0
for (i = 1; i <= n; i++) {
a = nearest(i)
if (assign[i] != a) { assign[i] = a; changed = 1 }
}
if (!changed) break
# Move each centroid to the mean of its members.
for (j = 1; j <= K; j++) { sx[j] = 0; sy[j] = 0; size[j] = 0 }
for (i = 1; i <= n; i++) {
a = assign[i]
sx[a] += px[i]; sy[a] += py[i]; size[a]++
}
for (j = 1; j <= K; j++)
if (size[j] > 0) { # an empty cluster keeps its centroid
cx[j] = sx[j] / size[j]
cy[j] = sy[j] / size[j]
}
}
# Sizes and inertia from the final assignment.
for (j = 1; j <= K; j++) size[j] = 0
for (i = 1; i <= n; i++) {
size[assign[i]]++
inertia += dist2(i, assign[i])
}
printf "iterations %d\n", iterations
for (j = 1; j <= K; j++) {
printf "centroid_%d_x %.6f\n", j, cx[j]
printf "centroid_%d_y %.6f\n", j, cy[j]
}
for (j = 1; j <= K; j++)
printf "size_%d %d\n", j, size[j]
printf "inertia %.6f\n", inertia
}The arithmetic is easy in awk. Assigning points and averaging clusters are
short loops, and the “did anything change?” check happens inside the loop
that makes the assignment. assign[i] starts out empty, so the first pass
always counts as a change with no special case.
Structure is where awk pushes back. POSIX awk has no arrays of arrays, so
points and centroids are pairs of parallel arrays, px, py, cx, and
cy, which the helper functions dist2 and nearest share as globals. And
since every point is needed before the first pass, all the real work sits
in END. awk has no clustering built in.
10 REM k-means clustering (Lloyd's algorithm), in classic BASIC.
20 REM
30 REM Starts from the first K points as centroids, then alternates two
40 REM steps until no point changes cluster: assign each point to its
50 REM nearest centroid, and move each centroid to the mean of its points.
60 REM
70 REM run.sh passes the data path in the DATA environment variable,
80 REM because bwbasic can't read command-line arguments.
100 K = 3 : MAXIT = 100
110 DIM X(1000), Y(1000), A(1000), FL$(20)
120 DIM CX(3), CY(3), SX(3), SY(3), SZ(3)
130 GOSUB 1000 : REM read the points into X(1) to X(N) and Y(1) to Y(N)
140 FOR J = 1 TO K
150 CX(J) = X(J) : CY(J) = Y(J) : REM copy the first K points
160 NEXT J
170 REM 0 means "no cluster yet", so the first pass always counts as a change.
180 FOR I = 1 TO N
190 A(I) = 0
200 NEXT I
210 IT = 0
220 IF IT >= MAXIT THEN 300
230 IT = IT + 1
240 GOSUB 2000 : REM assign each point; CH = 1 if any point moved
250 IF CH = 0 THEN 300
260 GOSUB 3000 : REM move each centroid to the mean of its points
270 GOTO 220
300 REM Count each cluster and add up the squared distances (inertia).
310 FOR J = 1 TO K
320 SZ(J) = 0
330 NEXT J
340 W = 0
350 FOR I = 1 TO N
360 J = A(I) : SZ(J) = SZ(J) + 1
370 W = W + (X(I) - CX(J)) ^ 2 + (Y(I) - CY(J)) ^ 2
380 NEXT I
390 PRINT "iterations"; IT
400 FOR J = 1 TO K
410 J$ = MID$(STR$(J), 2) : REM STR$ puts a space before the number
420 PRINT "centroid_"; J$; "_x ";
430 PRINT USING "#.######"; CX(J)
440 PRINT "centroid_"; J$; "_y ";
450 PRINT USING "#.######"; CY(J)
460 NEXT J
470 FOR J = 1 TO K
480 PRINT "size_"; MID$(STR$(J), 2); SZ(J)
490 NEXT J
500 PRINT "inertia ";
510 PRINT USING "#.######"; W
520 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read the x and y columns into X(1) to X(N) and Y(1) to Y(N).
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find which fields of the header are x and y
1050 FOR C = 1 TO NF
1060 IF FL$(C) = "x" THEN CXI = C
1070 IF FL$(C) = "y" THEN CYI = C
1080 NEXT C
1090 N = 0
1100 IF EOF(1) THEN 1150
1110 LINE INPUT #1, L$
1120 GOSUB 4000
1130 N = N + 1 : X(N) = VAL(FL$(CXI)) : Y(N) = VAL(FL$(CYI))
1140 GOTO 1100
1150 CLOSE #1
1160 RETURN
2000 REM Assign each point to its nearest centroid. Squared distance ranks
2010 REM centroids the same way as distance, so no SQR is needed, and only
2020 REM a strictly smaller distance wins, so a tie goes to the lower index.
2030 CH = 0
2040 FOR I = 1 TO N
2050 B = 1 : BD = (X(I) - CX(1)) ^ 2 + (Y(I) - CY(1)) ^ 2
2060 FOR J = 2 TO K
2070 D = (X(I) - CX(J)) ^ 2 + (Y(I) - CY(J)) ^ 2
2080 IF D < BD THEN B = J : BD = D
2090 NEXT J
2100 IF B <> A(I) THEN A(I) = B : CH = 1
2110 NEXT I
2120 RETURN
3000 REM Sum each cluster's points in one pass, then divide by its size.
3010 FOR J = 1 TO K
3020 SX(J) = 0 : SY(J) = 0 : SZ(J) = 0
3030 NEXT J
3040 FOR I = 1 TO N
3050 J = A(I)
3060 SX(J) = SX(J) + X(I) : SY(J) = SY(J) + Y(I) : SZ(J) = SZ(J) + 1
3070 NEXT I
3080 FOR J = 1 TO K
3090 REM An empty cluster keeps its old centroid.
3100 IF SZ(J) > 0 THEN CX(J) = SX(J) / SZ(J) : CY(J) = SY(J) / SZ(J)
3110 NEXT J
3120 RETURN
4000 REM Split L$ at its commas into FL$(1) to FL$(NF).
4010 NF = 0 : P = 1
4020 Q = INSTR(P, L$, ",")
4030 NF = NF + 1
4040 IF Q = 0 THEN FL$(NF) = MID$(L$, P) : RETURN
4050 FL$(NF) = MID$(L$, P, Q - P) : P = Q + 1
4060 GOTO 4020The points live in two parallel arrays, X and Y, and the centroids in
CX and CY, all given a fixed size with DIM. Each pass of the loop runs
two subroutines: GOSUB 2000 assigns every point to its nearest centroid
and sets CH to 1 if any point moved, and GOSUB 3000 moves each centroid
to the mean of its points. Squared distances rank centroids the same way as
distances, so we never call SQR, and only a strictly smaller distance
wins, so ties go to the lower cluster. STR$(1) gives " 1", with a
leading space, so MID$ trims it before building names like centroid_1_x.
BASIC has no k-means function to avoid.
/* k-means clustering (Lloyd's algorithm), using only the C standard library.
*
* Starts from the first K points as centroids, then alternates two steps until
* no point changes cluster: assign each point to its nearest centroid, and move
* each centroid to the mean of its points. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define K 3
#define MAX_ITERATIONS 100
#define MAX_LINE 1024
/* A struct groups the two coordinates, so a point can be copied with =. */
typedef struct { double x, y; } Point;
/* Return the 0-based position of name in a header line, or -1. strtok
* replaces each comma with '\0' as it goes, so it works on a copy. */
static int column_index(const char *header, const char *name) {
char copy[MAX_LINE];
strcpy(copy, header);
int index = 0;
for (char *field = strtok(copy, ",\r\n"); field; field = strtok(NULL, ",\r\n")) {
if (strcmp(field, name) == 0) return index;
index++;
}
return -1;
}
/* Read the x and y columns as points. The array doubles in size whenever it
* fills, and the count comes back through *count. */
static Point *read_points(const char *path, size_t *count) {
FILE *file = fopen(path, "r");
if (!file) { perror(path); exit(1); }
char line[MAX_LINE];
if (!fgets(line, sizeof line, file)) { fprintf(stderr, "empty file\n"); exit(1); }
int x_col = column_index(line, "x"), y_col = column_index(line, "y");
if (x_col < 0 || y_col < 0) { fprintf(stderr, "no x or y column\n"); exit(1); }
size_t n = 0, capacity = 64;
Point *points = malloc(capacity * sizeof *points);
while (fgets(line, sizeof line, file)) {
Point p = {0, 0};
int index = 0;
for (char *field = strtok(line, ",\r\n"); field; field = strtok(NULL, ",\r\n")) {
if (index == x_col) p.x = strtod(field, NULL);
if (index == y_col) p.y = strtod(field, NULL);
index++;
}
if (index == 0) continue; /* a blank line */
if (n == capacity) {
capacity *= 2;
points = realloc(points, capacity * sizeof *points);
}
points[n++] = p;
}
fclose(file);
*count = n;
return points;
}
/* Squared Euclidean distance. It ranks centroids in the same order as the
* distance itself, so the nearest one needs no sqrt. */
static double squared_distance(Point a, Point b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
/* The index of the closest centroid. Only a strictly smaller distance
* replaces the best so far, so a tie goes to the lower index. */
static int nearest(Point p, const Point centroids[K]) {
int best = 0;
for (int j = 1; j < K; j++)
if (squared_distance(p, centroids[j]) < squared_distance(p, centroids[best])) best = j;
return best;
}
int main(int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }
size_t n;
Point *points = read_points(argv[1], &n);
if (n < K) { fprintf(stderr, "need at least %d points\n", K); return 1; }
Point centroids[K];
for (int j = 0; j < K; j++) centroids[j] = points[j]; /* copies, not views */
/* -1 means "no cluster yet", so the first pass always counts as a change. */
int *assignment = malloc(n * sizeof *assignment);
for (size_t i = 0; i < n; i++) assignment[i] = -1;
int iterations = 0;
while (iterations < MAX_ITERATIONS) {
iterations++;
int changed = 0; /* C treats 0 as false and anything else as true */
for (size_t i = 0; i < n; i++) {
int j = nearest(points[i], centroids);
if (j != assignment[i]) { assignment[i] = j; changed = 1; }
}
if (!changed) break;
/* Sum each cluster's points in one pass, then divide by its size. */
Point sums[K] = {{0, 0}}; /* elements left out of an initializer start at 0 */
int sizes[K] = {0};
for (size_t i = 0; i < n; i++) {
sums[assignment[i]].x += points[i].x;
sums[assignment[i]].y += points[i].y;
sizes[assignment[i]]++;
}
for (int j = 0; j < K; j++)
if (sizes[j] > 0) /* an empty cluster keeps its old centroid */
centroids[j] = (Point){sums[j].x / sizes[j], sums[j].y / sizes[j]};
}
int sizes[K] = {0};
double inertia = 0;
for (size_t i = 0; i < n; i++) {
sizes[assignment[i]]++;
inertia += squared_distance(points[i], centroids[assignment[i]]);
}
/* Clusters are numbered from 0 inside C and from 1 in the output. */
printf("iterations %d\n", iterations);
for (int j = 0; j < K; j++) {
printf("centroid_%d_x %.6f\n", j + 1, centroids[j].x);
printf("centroid_%d_y %.6f\n", j + 1, centroids[j].y);
}
for (int j = 0; j < K; j++) printf("size_%d %d\n", j + 1, sizes[j]);
printf("inertia %.6f\n", inertia);
free(assignment);
free(points);
return 0;
}Each point is a small struct with x and y fields, so the data is one
array of Point that read_points grows with realloc. Assigning a struct
copies it, so centroids[j] = points[j] can’t overwrite a data point when a
centroid moves.
We notice that nothing changed with a flag: every assignment starts at −1,
and the pass sets changed whenever a point’s nearest centroid differs from
its stored one. nearest compares squared distances, which rank centroids
the same way, so it needs no sqrt, and only a strictly smaller distance
wins, so ties go to the lower index. C has no k-means function to avoid.
// k-means clustering (Lloyd's algorithm), using only the C++ standard library.
//
// Starts from the first K points as centroids, then alternates two steps until
// no point changes cluster: assign each point to its nearest centroid, and move
// each centroid to the mean of its points.
#include <algorithm>
#include <cstddef>
#include <fstream>
#include <print>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
constexpr std::size_t K = 3;
constexpr int MAX_ITERATIONS = 100;
// A struct groups the two coordinates, so a point can be copied with =.
struct Point {
double x;
double y;
};
// Split one line at every comma.
std::vector<std::string> split(const std::string& line) {
std::vector<std::string> fields;
std::stringstream stream(line);
std::string field;
while (std::getline(stream, field, ',')) {
fields.push_back(field);
}
return fields;
}
// The x and y columns of a CSV file, found by name, as points.
std::vector<Point> read_points(const std::string& path) {
std::ifstream file(path);
std::string line;
std::getline(file, line);
std::vector<std::string> header = split(line);
auto column = [&header](const std::string& name) {
auto found = std::find(header.begin(), header.end(), name);
if (found == header.end()) {
throw std::runtime_error("no column " + name);
}
return static_cast<std::size_t>(found - header.begin());
};
std::size_t x_col = column("x");
std::size_t y_col = column("y");
std::vector<Point> points;
while (std::getline(file, line)) {
if (!line.empty()) {
std::vector<std::string> fields = split(line);
points.push_back({std::stod(fields.at(x_col)), std::stod(fields.at(y_col))});
}
}
return points;
}
// Squared Euclidean distance. It ranks centroids in the same order as the
// distance itself, so the nearest one needs no sqrt.
double squared_distance(Point a, Point b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
// The index of the closest centroid. std::min_element returns the first of
// equal smallest values, so a tie goes to the lower index.
std::size_t nearest(Point p, const std::vector<Point>& centroids) {
auto closer = [p](Point a, Point b) {
return squared_distance(p, a) < squared_distance(p, b);
};
auto best = std::min_element(centroids.begin(), centroids.end(), closer);
return best - centroids.begin();
}
int main(int argc, char* argv[]) {
if (argc != 2) {
std::println(stderr, "usage: program data.csv");
return 1;
}
std::vector<Point> points = read_points(argv[1]);
if (points.size() < K) {
std::println(stderr, "need at least {} points", K);
return 1;
}
// This constructor copies the first K points, so moving a centroid
// never changes a data point.
std::vector<Point> centroids(points.begin(), points.begin() + K);
// Empty until the first pass, so that pass never matches it.
std::vector<std::size_t> assignment;
int iterations = 0;
while (iterations < MAX_ITERATIONS) {
iterations++;
std::vector<std::size_t> new_assignment;
for (Point p : points) {
new_assignment.push_back(nearest(p, centroids));
}
// == on two vectors compares their lengths and then every element.
if (new_assignment == assignment) {
break; // this unchanged pass still counts
}
assignment = new_assignment;
// Sum each cluster's points in one pass, then divide by its size.
std::vector<Point> sums(K, Point{0.0, 0.0});
std::vector<int> sizes(K, 0);
for (std::size_t i = 0; i < points.size(); i++) {
sums[assignment[i]].x += points[i].x;
sums[assignment[i]].y += points[i].y;
sizes[assignment[i]]++;
}
for (std::size_t j = 0; j < K; j++) {
if (sizes[j] > 0) { // an empty cluster keeps its old centroid
centroids[j] = {sums[j].x / sizes[j], sums[j].y / sizes[j]};
}
}
}
std::vector<int> sizes(K, 0);
double inertia = 0.0;
for (std::size_t i = 0; i < points.size(); i++) {
sizes[assignment[i]]++;
inertia += squared_distance(points[i], centroids[assignment[i]]);
}
// Clusters are numbered from 0 inside C++ and from 1 in the output.
std::println("iterations {}", iterations);
for (std::size_t j = 0; j < K; j++) {
std::println("centroid_{}_x {:.6f}", j + 1, centroids[j].x);
std::println("centroid_{}_y {:.6f}", j + 1, centroids[j].y);
}
for (std::size_t j = 0; j < K; j++) {
std::println("size_{} {}", j + 1, sizes[j]);
}
std::println("inertia {:.6f}", inertia);
}Each point is a small struct with x and y fields, as in C, but the
points, centroids, and assignment all live in a std::vector, which grows
as you push_back, so there’s no realloc.
Where C sets a flag as it goes, we notice that nothing changed by comparing
the whole assignment with the last one: == on two vectors checks their
lengths and then every element, and the first pass compares against an
empty vector, so it never matches. nearest hands squared distances to
std::min_element, which returns the first of equal smallest values, so
ties go to the lower index. C++ has no k-means function to avoid.
*> k-means clustering (Lloyd's algorithm), in COBOL.
*>
*> Starts from the first K = 3 points as centroids, then alternates
*> two steps until no point changes cluster: assign each point to
*> its nearest centroid, and move each centroid to the mean of its
*> points.
IDENTIFICATION DIVISION.
PROGRAM-ID. K-MEANS.
ENVIRONMENT DIVISION.
INPUT-OUTPUT SECTION.
FILE-CONTROL.
SELECT DATA-FILE ASSIGN TO DYNAMIC WS-FILE-NAME
ORGANIZATION IS LINE SEQUENTIAL.
DATA DIVISION.
FILE SECTION.
FD DATA-FILE.
01 DATA-LINE PIC X(200).
WORKING-STORAGE SECTION.
01 WS-FILE-NAME PIC X(256).
01 WS-K PIC 9 VALUE 3.
01 WS-MAX-ITERATIONS PIC 999 VALUE 100.
01 WS-EOF PIC X VALUE "N".
88 END-OF-FILE VALUE "Y".
*> One CSV line split into fields, and where the columns are.
01 WS-FIELDS.
05 WS-FIELD PIC X(40) OCCURS 20 TIMES.
01 WS-FIELD-COUNT PIC 99.
01 WS-POINTER PIC 999.
01 WS-LINE-LENGTH PIC 999.
01 WS-X-COLUMN PIC 99 VALUE 0.
01 WS-Y-COLUMN PIC 99 VALUE 0.
*> The points, each with the cluster it is currently assigned to.
01 WS-N PIC 9(4) VALUE 0.
01 WS-POINTS.
05 WS-POINT OCCURS 5000 TIMES.
10 WS-X PIC S9(9)V9(6).
10 WS-Y PIC S9(9)V9(6).
10 WS-CLUSTER PIC 9 VALUE 0.
*> The centroids, with the sums and sizes used to move them.
01 WS-CENTROIDS.
05 WS-CENTROID OCCURS 9 TIMES.
10 WS-CX PIC S9(9)V9(12).
10 WS-CY PIC S9(9)V9(12).
10 WS-SUM-X PIC S9(12)V9(6).
10 WS-SUM-Y PIC S9(12)V9(6).
10 WS-SIZE PIC 9(4).
01 I PIC 9(4).
01 J PIC 99.
01 WS-NEAREST PIC 9.
01 WS-DIST PIC S9(12)V9(12).
01 WS-BEST-DIST PIC S9(12)V9(12).
01 WS-ITERATIONS PIC 999 VALUE 0.
01 WS-CHANGED PIC X.
88 ASSIGNMENT-CHANGED VALUE "Y".
01 WS-INERTIA PIC S9(12)V9(12) VALUE 0.
*> Output: a name, then the value edited to six decimal places.
01 WS-NAME PIC X(20).
01 WS-RESULT PIC S9(12)V9(12).
01 WS-RESULT-OUT PIC -(12)9.9(6).
01 WS-COUNT-OUT PIC Z(8)9.
01 WS-INDEX-OUT PIC Z9.
PROCEDURE DIVISION.
MAIN.
ACCEPT WS-FILE-NAME FROM ARGUMENT-VALUE
OPEN INPUT DATA-FILE
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
END-READ
PERFORM SPLIT-LINE
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-FIELD-COUNT
IF WS-FIELD(I) = "x"
MOVE I TO WS-X-COLUMN
END-IF
IF WS-FIELD(I) = "y"
MOVE I TO WS-Y-COLUMN
END-IF
END-PERFORM
IF WS-X-COLUMN = 0 OR WS-Y-COLUMN = 0
DISPLAY "missing column x or y" UPON SYSERR
MOVE 1 TO RETURN-CODE
STOP RUN
END-IF
PERFORM UNTIL END-OF-FILE
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
NOT AT END PERFORM ADD-POINT
END-READ
END-PERFORM
CLOSE DATA-FILE
*> The first K points are the starting centroids.
PERFORM VARYING J FROM 1 BY 1 UNTIL J > WS-K
MOVE WS-X(J) TO WS-CX(J)
MOVE WS-Y(J) TO WS-CY(J)
END-PERFORM
*> Every assignment pass counts, including the final one that
*> changes nothing.
MOVE "Y" TO WS-CHANGED
PERFORM UNTIL NOT ASSIGNMENT-CHANGED
OR WS-ITERATIONS >= WS-MAX-ITERATIONS
ADD 1 TO WS-ITERATIONS
PERFORM ASSIGN-POINTS
IF ASSIGNMENT-CHANGED
PERFORM MOVE-CENTROIDS
END-IF
END-PERFORM
*> Inertia: the sum of squared distances to assigned centroids.
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
MOVE WS-CLUSTER(I) TO J
PERFORM DISTANCE
ADD WS-DIST TO WS-INERTIA
END-PERFORM
MOVE WS-ITERATIONS TO WS-COUNT-OUT
DISPLAY "iterations " FUNCTION TRIM(WS-COUNT-OUT)
PERFORM VARYING J FROM 1 BY 1 UNTIL J > WS-K
MOVE J TO WS-INDEX-OUT
MOVE SPACES TO WS-NAME
STRING "centroid_" FUNCTION TRIM(WS-INDEX-OUT) "_x"
DELIMITED BY SIZE INTO WS-NAME
MOVE WS-CX(J) TO WS-RESULT
PERFORM PRINT-RESULT
MOVE SPACES TO WS-NAME
STRING "centroid_" FUNCTION TRIM(WS-INDEX-OUT) "_y"
DELIMITED BY SIZE INTO WS-NAME
MOVE WS-CY(J) TO WS-RESULT
PERFORM PRINT-RESULT
END-PERFORM
PERFORM VARYING J FROM 1 BY 1 UNTIL J > WS-K
MOVE J TO WS-INDEX-OUT
MOVE WS-SIZE(J) TO WS-COUNT-OUT
DISPLAY "size_" FUNCTION TRIM(WS-INDEX-OUT) " "
FUNCTION TRIM(WS-COUNT-OUT)
END-PERFORM
MOVE "inertia" TO WS-NAME
MOVE WS-INERTIA TO WS-RESULT
PERFORM PRINT-RESULT
STOP RUN.
*> Split DATA-LINE on commas into WS-FIELD(1), WS-FIELD(2), ...
SPLIT-LINE.
MOVE SPACES TO WS-FIELDS
MOVE 0 TO WS-FIELD-COUNT
MOVE 1 TO WS-POINTER
MOVE FUNCTION LENGTH(FUNCTION TRIM(DATA-LINE TRAILING))
TO WS-LINE-LENGTH
PERFORM UNTIL WS-POINTER > WS-LINE-LENGTH
ADD 1 TO WS-FIELD-COUNT
UNSTRING DATA-LINE DELIMITED BY ","
INTO WS-FIELD(WS-FIELD-COUNT)
WITH POINTER WS-POINTER
END-UNSTRING
END-PERFORM.
*> Store one data line's point.
ADD-POINT.
PERFORM SPLIT-LINE
ADD 1 TO WS-N
COMPUTE WS-X(WS-N) = FUNCTION NUMVAL(WS-FIELD(WS-X-COLUMN))
COMPUTE WS-Y(WS-N) = FUNCTION NUMVAL(WS-FIELD(WS-Y-COLUMN)).
*> Squared distance from point I to centroid J, into WS-DIST.
*> Squared distance picks the same nearest centroid as distance.
DISTANCE.
COMPUTE WS-DIST = (WS-X(I) - WS-CX(J)) ** 2
+ (WS-Y(I) - WS-CY(J)) ** 2.
*> Assign every point to its nearest centroid; note any change.
ASSIGN-POINTS.
MOVE "N" TO WS-CHANGED
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
MOVE 1 TO WS-NEAREST
MOVE 1 TO J
PERFORM DISTANCE
MOVE WS-DIST TO WS-BEST-DIST
*> A strict "<" sends a tie to the lower centroid index.
PERFORM VARYING J FROM 2 BY 1 UNTIL J > WS-K
PERFORM DISTANCE
IF WS-DIST < WS-BEST-DIST
MOVE J TO WS-NEAREST
MOVE WS-DIST TO WS-BEST-DIST
END-IF
END-PERFORM
IF WS-NEAREST NOT = WS-CLUSTER(I)
MOVE WS-NEAREST TO WS-CLUSTER(I)
MOVE "Y" TO WS-CHANGED
END-IF
END-PERFORM.
*> Move each centroid to the mean of its points. An empty cluster
*> keeps its old centroid.
MOVE-CENTROIDS.
PERFORM VARYING J FROM 1 BY 1 UNTIL J > WS-K
MOVE 0 TO WS-SUM-X(J) WS-SUM-Y(J) WS-SIZE(J)
END-PERFORM
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
MOVE WS-CLUSTER(I) TO J
ADD WS-X(I) TO WS-SUM-X(J)
ADD WS-Y(I) TO WS-SUM-Y(J)
ADD 1 TO WS-SIZE(J)
END-PERFORM
PERFORM VARYING J FROM 1 BY 1 UNTIL J > WS-K
IF WS-SIZE(J) > 0
COMPUTE WS-CX(J) ROUNDED = WS-SUM-X(J) / WS-SIZE(J)
COMPUTE WS-CY(J) ROUNDED = WS-SUM-Y(J) / WS-SIZE(J)
END-IF
END-PERFORM.
*> Print "name value", rounded to six decimal places.
PRINT-RESULT.
COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
DISPLAY FUNCTION TRIM(WS-NAME) " "
FUNCTION TRIM(WS-RESULT-OUT).The algorithm maps neatly onto COBOL paragraphs: ASSIGN-POINTS,
MOVE-CENTROIDS, and a DISTANCE helper, repeated with PERFORM UNTIL
until nothing changes. Each point’s row in the table stores its cluster
number, so a change is noticed the moment it happens. Numbers are
fixed-point decimals, with centroids held to 12 places in PIC S9(9)V9(12).
Paragraphs take no arguments, though. DISTANCE works on whatever the
shared I and J hold, so every caller sets them first. Printing names
like centroid_2_x takes a STRING statement and an edited picture to
drop the leading zero. The result is 227 lines, about four times the
Python’s 55. GnuCOBOL has MEAN and VARIANCE functions, but nothing for
clustering.
;;;; k-means clustering (Lloyd's algorithm), using only standard Common Lisp.
;;;;
;;;; Starts from the first k points as centroids, then alternates two steps
;;;; until no point changes cluster: assign each point to its nearest
;;;; centroid, and move each centroid to the mean of its points.
(defparameter *k* 3)
(defparameter *max-iterations* 100)
;; Read 1.25 as a double, not a single float with about 7 digits.
;; eval-when makes the setting apply while the compiler reads this file, too,
;; not only once the program runs.
(eval-when (:compile-toplevel :load-toplevel :execute)
(setf *read-default-float-format* 'double-float))
;; defstruct defines a record type: make-point builds one, and point-x and
;; point-y read its fields.
(defstruct point x y)
(defun split-on-comma (line)
"The fields of LINE, split at every comma."
(loop for start = 0 then (1+ end)
for end = (position #\, line :start start)
collect (subseq line start end)
while end))
(defun parse-number (text)
"The number written in TEXT, read by the Lisp reader. Binding *read-eval*
to nil stops the reader from running code written as #.(...)."
(let ((*read-eval* nil))
(coerce (read-from-string text) 'double-float)))
(defun read-points (path)
"The x and y columns of the CSV file at PATH, as a list of points."
(with-open-file (in path)
(let* ((header (split-on-comma (read-line in)))
(xi (or (position "x" header :test #'string=) (error "no column x")))
(yi (or (position "y" header :test #'string=) (error "no column y"))))
(loop for line = (read-line in nil)
while line
unless (string= line "")
collect (let ((fields (split-on-comma line)))
(make-point :x (parse-number (nth xi fields))
:y (parse-number (nth yi fields))))))))
(defun squared-distance (p q)
"Squared Euclidean distance. It ranks centroids in the same order as the
distance itself, so finding the nearest one needs no sqrt."
(+ (expt (- (point-x p) (point-x q)) 2)
(expt (- (point-y p) (point-y q)) 2)))
(defun nearest (point centroids)
"The index of the centroid closest to POINT, counting from 0. Only a
strictly smaller distance replaces the best so far, so a tie goes to the
lower index."
(loop with best = 0
with best-distance = nil
for centroid in centroids
for j from 0
for d = (squared-distance point centroid)
when (or (null best-distance) (< d best-distance))
do (setf best j best-distance d)
finally (return best)))
(defun members (j points assignment)
"The points assigned to cluster J."
(loop for p in points
for a in assignment
when (= a j) collect p))
(defun mean (numbers)
(/ (reduce #'+ numbers) (length numbers)))
(defun update (points assignment centroids)
"Each centroid moved to the mean of its points. An empty cluster keeps its
old centroid."
(loop for centroid in centroids
for j from 0
for group = (members j points assignment)
collect (if (null group)
centroid
(make-point :x (mean (mapcar #'point-x group))
:y (mean (mapcar #'point-y group))))))
(let* ((points (read-points (second sb-ext:*posix-argv*)))
;; subseq copies the first k points into a new list.
(centroids (subseq points 0 *k*))
(assignment nil)
(iterations 0))
;; Before the first pass the assignment is nil, the empty list, which no
;; real assignment equals. equal compares two lists element by element;
;; eq would only ask whether they are the same object in memory.
(loop while (< iterations *max-iterations*)
do (incf iterations)
(let ((new (mapcar (lambda (p) (nearest p centroids)) points)))
(when (equal new assignment)
(return))
(setf assignment new
centroids (update points assignment centroids))))
(let ((inertia (loop for p in points
for a in assignment
sum (squared-distance p (nth a centroids)))))
(format t "iterations ~d~%" iterations)
;; Clusters are numbered from 0 inside the program and from 1 in the output.
(loop for c in centroids
for j from 1
do (format t "centroid_~d_x ~,6f~%" j (point-x c))
(format t "centroid_~d_y ~,6f~%" j (point-y c)))
(loop for j from 0 below *k*
do (format t "size_~d ~d~%" (1+ j) (count j assignment)))
(format t "inertia ~,6f~%" inertia)))Each point is a defstruct record, so (point-x p) reads its x. Assignment
is a mapcar of nearest over the points, and nearest keeps a centroid
only if it is strictly closer, so a tie would go to the lower index. To
notice that nothing changed, we compare the whole new assignment with the
last one using equal, which checks two lists element by element; eq
would only ask whether they are the same list in memory.
Common Lisp has no clustering, mean, or distance function to skip, so we write the squared distance and the mean ourselves. An empty cluster would keep its old centroid, though this data never produces one.
# k-means clustering (Lloyd's algorithm), using only Elixir's standard
# library.
#
# Starts from the first k points as centroids, then alternates two steps
# until no point changes cluster: assign each point to its nearest centroid,
# and move each centroid to the mean of its points. Elixir has no while
# loop, so the repetition is a function that calls itself.
defmodule KMeans do
@k 3
@max_iterations 100
# The x and y columns as a list of {x, y} tuples. The header gives their
# positions.
def read_points(path) do
[header | rows] = path |> File.read!() |> String.split("\n", trim: true)
names = String.split(header, ",")
xi = Enum.find_index(names, &(&1 == "x"))
yi = Enum.find_index(names, &(&1 == "y"))
Enum.map(rows, fn row ->
fields = String.split(row, ",")
# Float.parse returns {number, rest}; matching "" rejects trailing junk.
{x, ""} = fields |> Enum.at(xi) |> Float.parse()
{y, ""} = fields |> Enum.at(yi) |> Float.parse()
{x, y}
end)
end
# Squared Euclidean distance. It ranks centroids in the same order as the
# distance itself, so finding the nearest one needs no sqrt.
def squared_distance({ax, ay}, {bx, by}), do: (ax - bx) ** 2 + (ay - by) ** 2
# The index of the closest centroid. Enum.min_by keeps the first of equal
# distances, so a tie goes to the lower index.
def nearest(point, centroids) do
centroids
|> Enum.with_index()
|> Enum.min_by(fn {c, _j} -> squared_distance(point, c) end)
|> elem(1)
end
# Move each centroid to the mean of its points. An empty cluster keeps its
# old centroid.
def update(points, assignment, centroids) do
centroids
|> Enum.with_index()
|> Enum.map(fn {c, j} ->
members = for {p, ^j} <- Enum.zip(points, assignment), do: p
case members do
[] -> c
_ ->
n = length(members)
{members |> Enum.map(&elem(&1, 0)) |> Enum.sum() |> Kernel./(n),
members |> Enum.map(&elem(&1, 1)) |> Enum.sum() |> Kernel./(n)}
end
end)
end
# One assignment step, then either stop or update and go again. previous
# is nil before the first pass, so the first comparison always counts as a
# change. == compares two lists element by element. Returns the number of
# assignment steps, the final assignment, and the centroids.
def run(points, iterations, previous, centroids) do
step = iterations + 1
assignment = Enum.map(points, &nearest(&1, centroids))
cond do
assignment == previous -> {step, assignment, centroids}
step == @max_iterations -> {step, assignment, update(points, assignment, centroids)}
true -> run(points, step, assignment, update(points, assignment, centroids))
end
end
def k, do: @k
end
[path] = System.argv()
points = KMeans.read_points(path)
{iterations, assignment, centroids} = KMeans.run(points, 0, nil, Enum.take(points, KMeans.k()))
six = fn v -> :erlang.float_to_binary(v, decimals: 6) end
IO.puts("iterations #{iterations}")
# Clusters are numbered from 0 inside the program and from 1 in the output.
centroids
|> Enum.with_index(1)
|> Enum.each(fn {{x, y}, j} ->
IO.puts("centroid_#{j}_x #{six.(x)}")
IO.puts("centroid_#{j}_y #{six.(y)}")
end)
for j <- 0..(KMeans.k() - 1) do
IO.puts("size_#{j + 1} #{Enum.count(assignment, &(&1 == j))}")
end
inertia =
Enum.zip(points, assignment)
|> Enum.map(fn {p, a} -> KMeans.squared_distance(p, Enum.at(centroids, a)) end)
|> Enum.sum()
IO.puts("inertia #{six.(inertia)}")We hold each point as a tuple, {x, y}, and the centroids as a list of
them. Elixir has no while loop, so run calls itself with the moved
centroids until a fresh assignment list == the previous one.
Elixir and Erlang ship no k-means, and Hex packages such as Nx aren’t
installed. The steps use general tools: Enum.min_by picks the nearest
centroid and keeps the first of equal distances, so ties go to the lower
index. Enum.sum and / give each mean. A comprehension, for {p, ^j} <- ...,
gathers a cluster’s members; the ^ pins j so only matching pairs pass.
Distances stay squared, which ranks centroids the same way without a square
root.
! k-means clustering (Lloyd's algorithm), using only standard Fortran.
!
! Starts from the first K points as centroids, then alternates two steps until
! no point changes cluster: assign each point to its nearest centroid, and move
! each centroid to the mean of its points.
program k_means
implicit none
integer, parameter :: dp = kind(1.0d0) ! double precision
integer, parameter :: k = 3, max_iterations = 100
real(dp), allocatable :: points(:, :)
real(dp) :: centroids(2, k), inertia
integer, allocatable :: assignment(:), new_assignment(:)
integer :: i, j, n, iterations, members
character(len=20) :: name
! Fortran stores arrays column by column, so each point is a column:
! points(:, i) is the x and y of point i.
points = read_points()
n = size(points, 2)
centroids = points(:, 1:k) ! assignment copies; it is not a view
allocate(assignment(n), new_assignment(n))
assignment = 0 ! no cluster yet, so the first pass always differs
iterations = 0
do while (iterations < max_iterations)
iterations = iterations + 1
do i = 1, n
new_assignment(i) = nearest_centroid(points(:, i))
end do
if (all(new_assignment == assignment)) exit
assignment = new_assignment
do j = 1, k
members = count(assignment == j)
if (members > 0) then ! an empty cluster keeps its old centroid
! sum with a mask adds only the elements where the mask is true.
centroids(1, j) = sum(points(1, :), mask=assignment == j) / members
centroids(2, j) = sum(points(2, :), mask=assignment == j) / members
end if
end do
end do
inertia = 0
do i = 1, n
inertia = inertia + sum((points(:, i) - centroids(:, assignment(i)))**2)
end do
print "(a, 1x, i0)", "iterations", iterations
do j = 1, k
write(name, "(a, i0)") "centroid_", j ! builds the name in a string
call say(trim(name) // "_x", centroids(1, j))
call say(trim(name) // "_y", centroids(2, j))
end do
do j = 1, k
write(name, "(a, i0)") "size_", j
print "(a, 1x, i0)", trim(name), count(assignment == j)
end do
call say("inertia", inertia)
contains
! The index of the centroid closest to a point. minloc returns the first
! of equal minimums, so a tie goes to the lower index.
integer function nearest_centroid(point)
real(dp), intent(in) :: point(2)
real(dp) :: distance(k)
integer :: j
do j = 1, k
distance(j) = norm2(point - centroids(:, j)) ! Euclidean distance
end do
nearest_centroid = minloc(distance, dim=1)
end function nearest_centroid
! Every row of the file as a column of a 2-by-n array.
function read_points() result(points)
real(dp), allocatable :: points(:, :)
real(dp), allocatable :: xs(:), ys(:)
character(len=256) :: path, line, text
integer :: unit, status, x_col, y_col
real(dp) :: x, y
call get_command_argument(1, path)
open(newunit=unit, file=path, status="old", action="read")
read(unit, "(a)") line
x_col = field_index(line, "x")
y_col = field_index(line, "y")
allocate(xs(0), ys(0))
do
read(unit, "(a)", iostat=status) line
if (status /= 0) exit
text = field(line, x_col)
read(text, *) x ! an internal read parses the number
text = field(line, y_col)
read(text, *) y
xs = [xs, x] ! grows the array by one
ys = [ys, y]
end do
close(unit)
allocate(points(2, size(xs)))
points(1, :) = xs ! the first row holds every x
points(2, :) = ys
end function read_points
! The position of a name in a comma-separated header line.
integer function field_index(header, name)
character(len=*), intent(in) :: header, name
integer :: i, fields
fields = 1
do i = 1, len_trim(header)
if (header(i:i) == ",") fields = fields + 1
end do
do field_index = 1, fields
if (field(header, field_index) == name) return
end do
error stop "column not found: " // name
end function field_index
! The k-th comma-separated field of a line.
function field(line, k) result(text)
character(len=*), intent(in) :: line
integer, intent(in) :: k
character(len=:), allocatable :: text
integer :: start, i, comma
start = 1
do i = 1, k - 1
start = start + index(line(start:), ",")
end do
comma = index(line(start:), ",")
if (comma == 0) then
text = trim(line(start:))
else
text = line(start:start + comma - 2)
end if
end function field
! Print "name value" with six decimals. The f0.6 format drops the zero
! before the point (0.5 prints as .500000), so we put it back.
subroutine say(name, x)
character(len=*), intent(in) :: name
real(dp), intent(in) :: x
character(len=40) :: text
write(text, "(f0.6)") x
if (text(1:1) == ".") text = "0" // text
if (text(1:2) == "-.") text = "-0" // text(2:)
print "(a, 1x, a)", name, trim(text)
end subroutine say
end program k_meansFortran stores arrays column by column, so each point is a column of a
2-by-150 array and points(:, i) is one point. centroids = points(:, 1:k)
copies the values, so moving a centroid can’t overwrite a data point.
minloc returns the first of equal distances, so a tie goes to the
lower-numbered centroid, and clusters are numbered from 1, matching the output.
A new centroid is sum(points(1, :), mask=assignment == j) over the member
count, skipped when the count is zero. The helper is nearest_centroid
because nearest is already an intrinsic: it returns the next representable
number. gfortran has no k-means intrinsic, so there was nothing to skip.
// k-means clustering (Lloyd's algorithm), using only the standard library.
//
// Starts from the first K points as centroids, then alternates two steps until
// no point changes cluster: assign each point to its nearest centroid, and move
// each centroid to the mean of its points.
package main
import (
"encoding/csv"
"fmt"
"math"
"os"
"slices"
"strconv"
)
const (
k = 3
maxIterations = 100
)
type point struct{ x, y float64 }
func readPoints(path string) ([]point, error) {
f, err := os.Open(path)
if err != nil {
return nil, err
}
defer f.Close()
records, err := csv.NewReader(f).ReadAll()
if err != nil {
return nil, err
}
xi, yi := slices.Index(records[0], "x"), slices.Index(records[0], "y")
if xi < 0 || yi < 0 {
return nil, fmt.Errorf("need columns x and y")
}
points := make([]point, 0, len(records)-1)
for _, record := range records[1:] {
x, err := strconv.ParseFloat(record[xi], 64)
if err != nil {
return nil, err
}
y, err := strconv.ParseFloat(record[yi], 64)
if err != nil {
return nil, err
}
points = append(points, point{x, y})
}
return points, nil
}
func distance(a, b point) float64 {
return math.Hypot(a.x-b.x, a.y-b.y)
}
// nearest returns the index of the closest centroid. The comparison is strict,
// so a tie goes to the lower index.
func nearest(p point, centroids []point) int {
best := 0
for j := 1; j < len(centroids); j++ {
if distance(p, centroids[j]) < distance(p, centroids[best]) {
best = j
}
}
return best
}
func main() {
points, err := readPoints(os.Args[1])
if err != nil {
fmt.Fprintln(os.Stderr, err)
os.Exit(1)
}
// Clone, so moving a centroid never overwrites one of the points.
centroids := slices.Clone(points[:k])
var assignment []int // nil until the first pass
iterations := 0
for iterations < maxIterations {
iterations++
newAssignment := make([]int, len(points))
for i, p := range points {
newAssignment[i] = nearest(p, centroids)
}
if slices.Equal(newAssignment, assignment) {
break
}
assignment = newAssignment
var sums [k]point
var counts [k]int
for i, p := range points {
j := assignment[i]
sums[j].x += p.x
sums[j].y += p.y
counts[j]++
}
for j := range k {
if counts[j] > 0 { // an empty cluster keeps its old centroid
n := float64(counts[j])
centroids[j] = point{sums[j].x / n, sums[j].y / n}
}
}
}
var sizes [k]int
inertia := 0.0
for i, p := range points {
j := assignment[i]
sizes[j]++
d := distance(p, centroids[j])
inertia += d * d
}
fmt.Printf("iterations %d\n", iterations)
for j, c := range centroids {
fmt.Printf("centroid_%d_x %.6f\n", j+1, c.x)
fmt.Printf("centroid_%d_y %.6f\n", j+1, c.y)
}
for j, size := range sizes {
fmt.Printf("size_%d %d\n", j+1, size)
}
fmt.Printf("inertia %.6f\n", inertia)
}A point struct with x and y fields keeps the loop readable, and
math.Hypot gives the distance. The nearest centroid comes from a loop with
a strict <, so a tie would go to the lower index. slices.Equal tells us
when the assignments stopped changing; it treats the nil slice we start
with as unequal to any 150 assignments, so the first pass always continues.
The trap is aliasing. points[:3] isn’t a copy: it shares memory with
points, so moving a centroid would overwrite a data point. slices.Clone
makes a real copy. As in Python, an empty cluster keeps its old centroid,
and there is no built-in clustering to skip.
-- k-means clustering (Lloyd's algorithm), using only Haskell's base library.
--
-- Starts from the first k points as centroids, then alternates two steps until
-- no point changes cluster: assign each point to its nearest centroid, and move
-- each centroid to the mean of its points. Haskell has no while loop, so the
-- repetition is a function that calls itself with the new centroids.
import Data.List (elemIndex, minimumBy)
import Data.Ord (comparing)
import System.Environment (getArgs)
import Text.Printf (printf)
-- A point is a pair of coordinates, (x, y).
type Point = (Double, Double)
k :: Int
k = 3
maxIterations :: Int
maxIterations = 100
-- Split a line at every comma. break returns the text before the first comma
-- and the rest, which still starts with the comma.
splitOnComma :: String -> [String]
splitOnComma line = case break (== ',') line of
(field, ',' : rest) -> field : splitOnComma rest
(field, _) -> [field]
-- The x and y columns as points. The header gives their positions.
readPoints :: String -> [Point]
readPoints text = case lines text of
header : rows
| Just xi <- elemIndex "x" (splitOnComma header)
, Just yi <- elemIndex "y" (splitOnComma header) ->
[ (read (fields !! xi), read (fields !! yi))
| row <- rows, not (null row), let fields = splitOnComma row ]
_ -> error "no x or y column"
-- Squared Euclidean distance. It ranks centroids in the same order as the
-- distance itself, so finding the nearest one needs no sqrt.
squaredDistance :: Point -> Point -> Double
squaredDistance (x1, y1) (x2, y2) = (x1 - x2) ^ (2 :: Int) + (y1 - y2) ^ (2 :: Int)
-- The index of the closest centroid. zip [0 ..] numbers the centroids, and
-- minimumBy keeps the first of equal distances, so a tie goes to the lower index.
nearest :: [Point] -> Point -> Int
nearest centroids p =
fst (minimumBy (comparing snd) [(j, squaredDistance p c) | (j, c) <- zip [0 ..] centroids])
-- Move each centroid to the mean of its points. An empty cluster keeps its
-- old centroid.
update :: [Point] -> [Int] -> [Point] -> [Point]
update points assignment centroids =
[ if null members then c else (mean (map fst members), mean (map snd members))
| (j, c) <- zip [0 ..] centroids
, let members = [p | (p, a) <- zip points assignment, a == j] ]
where
mean values = sum values / fromIntegral (length values)
-- One assignment step, then either stop or update and go again. The previous
-- assignment is Maybe a list: Nothing before the first pass, so the first
-- comparison always counts as a change. Returns the number of assignment
-- steps, the final assignment, and the final centroids.
run :: [Point] -> Int -> Maybe [Int] -> [Point] -> (Int, [Int], [Point])
run points iterations previous centroids
| Just assignment == previous = (step, assignment, centroids)
| step == maxIterations = (step, assignment, moved)
| otherwise = run points step (Just assignment) moved
where
step = iterations + 1
assignment = map (nearest centroids) points
moved = update points assignment centroids
main :: IO ()
main = do
[path] <- getArgs
text <- readFile path
let points = readPoints text
(iterations, assignment, centroids) = run points 0 Nothing (take k points)
sizes = [length (filter (== j) assignment) | j <- [0 .. k - 1]]
inertia = sum [squaredDistance p (centroids !! a) | (p, a) <- zip points assignment]
printf "iterations %d\n" iterations
-- Clusters are numbered from 0 inside the program and from 1 in the output.
mapM_ (\(j, (x, y)) -> printf "centroid_%d_x %.6f\ncentroid_%d_y %.6f\n" j x j y)
(zip [1 :: Int ..] centroids)
mapM_ (\(j, size) -> printf "size_%d %d\n" j size) (zip [1 :: Int ..] sizes)
printf "inertia %.6f\n" inertiaPoints are pairs, (Double, Double), and the assignment is a plain list of
cluster numbers. With no loops, run repeats by calling itself with the
moved centroids. Its previous assignment is a Maybe [Int]: Nothing
before the first pass, then Just the last list, so “nothing changed” is
one == between whole lists.
nearest numbers the centroids with zip [0 ..] and takes minimumBy on
squared distance. minimumBy keeps the first of equal values, so a tie goes
to the lower index; maximumBy keeps the last, so it would need care. The
exponents are typed (2 :: Int), or -Wall warns that it guessed.
// k-means clustering (Lloyd's algorithm), using only the Java standard
// library.
//
// Starts from the first K points as centroids, then alternates two steps
// until no point changes cluster: assign each point to its nearest centroid,
// and move each centroid to the mean of its points.
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
class Main {
static final int K = 3;
static final int MAX_ITERATIONS = 100;
// A record is a small immutable class: Java writes the constructor, the
// x() and y() accessors, and equals for us.
record Point(double x, double y) {}
// The named column of a CSV file, as an array of doubles. Reading a file
// can fail with an IOException, and Java makes us say so with "throws".
static double[] readColumn(Path path, String name) throws IOException {
List<String> lines = Files.readAllLines(path);
int index = Arrays.asList(lines.get(0).split(",")).indexOf(name);
if (index < 0) {
throw new IllegalArgumentException("no column " + name);
}
double[] values = new double[lines.size() - 1];
for (int row = 1; row < lines.size(); row++) {
values[row - 1] = Double.parseDouble(lines.get(row).split(",")[index]);
}
return values;
}
static double distance(Point a, Point b) {
return Math.hypot(a.x() - b.x(), a.y() - b.y());
}
// The index of the closest centroid. The comparison is strict, so a tie
// goes to the lower index.
static int nearest(Point p, Point[] centroids) {
int best = 0;
for (int j = 1; j < centroids.length; j++) {
if (distance(p, centroids[j]) < distance(p, centroids[best])) {
best = j;
}
}
return best;
}
public static void main(String[] args) throws IOException {
Path path = Path.of(args[0]);
double[] xs = readColumn(path, "x");
double[] ys = readColumn(path, "y");
var points = new Point[xs.length];
for (int i = 0; i < points.length; i++) {
points[i] = new Point(xs[i], ys[i]);
}
// Records can't be changed, so moving a centroid replaces it in this
// array and never touches the points.
Point[] centroids = Arrays.copyOf(points, K);
int[] assignment = null; // nothing assigned before the first pass
int iterations = 0;
while (iterations < MAX_ITERATIONS) {
iterations++;
var newAssignment = new int[points.length];
for (int i = 0; i < points.length; i++) {
newAssignment[i] = nearest(points[i], centroids);
}
// Arrays.equals compares the contents; == would only ask whether
// both names point at the same array. Against null it is false.
if (Arrays.equals(newAssignment, assignment)) {
break;
}
assignment = newAssignment;
var sumX = new double[K];
var sumY = new double[K];
var counts = new int[K];
for (int i = 0; i < points.length; i++) {
int j = assignment[i];
sumX[j] += points[i].x();
sumY[j] += points[i].y();
counts[j]++;
}
for (int j = 0; j < K; j++) {
if (counts[j] > 0) { // an empty cluster keeps its old centroid
centroids[j] = new Point(sumX[j] / counts[j], sumY[j] / counts[j]);
}
}
}
var sizes = new int[K];
double inertia = 0;
for (int i = 0; i < points.length; i++) {
int j = assignment[i];
sizes[j]++;
double d = distance(points[i], centroids[j]);
inertia += d * d;
}
// printf follows the computer's language settings unless told
// otherwise; Locale.ROOT keeps the decimal point a point.
Locale root = Locale.ROOT;
System.out.printf(root, "iterations %d%n", iterations);
for (int j = 0; j < K; j++) {
System.out.printf(root, "centroid_%d_x %.6f%n", j + 1, centroids[j].x());
System.out.printf(root, "centroid_%d_y %.6f%n", j + 1, centroids[j].y());
}
for (int j = 0; j < K; j++) {
System.out.printf(root, "size_%d %d%n", j + 1, sizes[j]);
}
System.out.printf(root, "inertia %.6f%n", inertia);
}
}A one-line Point record holds each pair, and Java writes the constructor,
the x() and y() accessors, and equals for us. Math.hypot gives the
distance, and a loop with a strict < picks the nearest centroid, so a tie
would go to the lower index. Each pass builds a fresh int[] of
assignments, and Arrays.equals compares it with the last one. Don’t use
== or .equals on arrays: both ask whether the two names point at the same
array, so equal assignments would never match. Against the null
we start with, Arrays.equals is false, so the first pass always continues.
Records can’t change, so a moving centroid is a new Point, and there is no
built-in clustering to skip.
// k-means clustering (Lloyd's algorithm), using only Node's built-in modules.
// Starts from the first K points as centroids, then alternates two steps until
// no point changes cluster: assign each point to its nearest centroid, and move
// each centroid to the mean of its points.
import { readFileSync } from "node:fs";
const K = 3;
const MAX_ITERATIONS = 100;
// Node has no CSV parser. No field here is quoted, so we split on commas and
// find the x and y columns by their header names.
function readPoints(path) {
const [header, ...lines] = readFileSync(path, "utf8").trim().split("\n");
const names = header.split(",");
const ix = names.indexOf("x");
const iy = names.indexOf("y");
return lines.map((line) => {
const fields = line.split(",");
return [Number(fields[ix]), Number(fields[iy])];
});
}
function distance([x1, y1], [x2, y2]) {
return Math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2);
}
function nearest(point, centroids) {
// Only a strictly smaller distance replaces the best so far, so ties go to
// the lower index.
let best = 0;
for (let j = 1; j < centroids.length; j++) {
if (distance(point, centroids[j]) < distance(point, centroids[best])) best = j;
}
return best;
}
// == on two arrays compares identity, not contents, so compare element by element.
const sameAssignment = (a, b) => b !== null && a.every((cluster, i) => cluster === b[i]);
const points = readPoints(process.argv[2]);
const centroids = points.slice(0, K).map(([x, y]) => [x, y]);
let assignment = null;
let iterations = 0;
while (iterations < MAX_ITERATIONS) {
iterations++;
const newAssignment = points.map((p) => nearest(p, centroids));
if (sameAssignment(newAssignment, assignment)) break;
assignment = newAssignment;
for (let j = 0; j < K; j++) {
const members = points.filter((_, i) => assignment[i] === j);
if (members.length > 0) {
// An empty cluster keeps its old centroid.
const sumX = members.reduce((acc, [x]) => acc + x, 0);
const sumY = members.reduce((acc, [, y]) => acc + y, 0);
centroids[j] = [sumX / members.length, sumY / members.length];
}
}
}
const sizes = centroids.map((_, j) => assignment.filter((a) => a === j).length);
const inertia = points.reduce((acc, p, i) => acc + distance(p, centroids[assignment[i]]) ** 2, 0);
console.log(`iterations ${iterations}`);
centroids.forEach(([x, y], j) => {
console.log(`centroid_${j + 1}_x ${x.toFixed(6)}`);
console.log(`centroid_${j + 1}_y ${y.toFixed(6)}`);
});
sizes.forEach((size, j) => console.log(`size_${j + 1} ${size}`));
console.log(`inertia ${inertia.toFixed(6)}`);Points are two-element arrays, and destructuring in the parameter list,
function distance([x1, y1], [x2, y2]), unpacks them where they arrive.
nearest is a plain loop, because arrays have no min with a key;
replacing the best only on a strictly smaller distance sends a tie to the
lower index.
The trap is the stopping test. == on two arrays compares identity, not
contents: [1, 2] == [1, 2] is false. Written the Python way, the test
never sees the assignments settle, and the loop runs all 100 iterations
before it stops. sameAssignment compares element by element instead. The
standard library has no clustering to skip.
# k-means clustering (Lloyd's algorithm), using only Base Julia.
#
# Starts from the first K points as centroids, then alternates two steps until
# no point changes cluster: assign each point to its nearest centroid, and move
# each centroid to the mean of its points.
using Printf
const K = 3
const MAX_ITERATIONS = 100
distance(a, b) = sqrt(sum((a .- b) .^ 2))
# argmin returns the first of equal distances, so ties go to the lower index.
nearest(point, centroids) = argmin([distance(point, c) for c in centroids])
function main(path)
lines = Iterators.drop(eachline(path), 1) # skip the header
points = [parse.(Float64, split(line, ",")) for line in lines]
centroids = points[1:K] # a slice is a new vector, so points is never changed
assignment = Int[]
iterations = 0
while iterations < MAX_ITERATIONS
iterations += 1
new_assignment = [nearest(p, centroids) for p in points]
new_assignment == assignment && break
assignment = new_assignment
for j in 1:K
members = points[assignment .== j]
if !isempty(members) # an empty cluster keeps its old centroid
centroids[j] = sum(members) / length(members)
end
end
end
inertia = sum(distance(p, centroids[a])^2 for (p, a) in zip(points, assignment))
println("iterations $iterations")
for (j, (x, y)) in enumerate(centroids)
@printf("centroid_%d_x %.6f\n", j, x)
@printf("centroid_%d_y %.6f\n", j, y)
end
for j in 1:K
println("size_$j $(count(==(j), assignment))")
end
@printf("inertia %.6f\n", inertia)
end
main(ARGS[1])Points are vectors, and Julia’s arithmetic works on whole vectors: sum
over a cluster’s members adds them coordinate by coordinate, so the new
centroid is just sum(members) / length(members). Picking the members is a
mask, points[assignment .== j], which keeps the points whose assignment
equals j.
argmin returns the first of equal values, so a point tied between two
centroids goes to the lower-numbered one, as in Python. Clusters are
numbered from 1, matching the output lines. Neither Base nor Statistics
has a k-means function, so there was nothing to skip.
-- k-means clustering (Lloyd's algorithm), using only Lua's standard library.
--
-- Lua has no clustering, so both steps are written out. Starting from the
-- first K points as centroids, we assign each point to its nearest centroid,
-- move each centroid to the mean of its points, and stop when no point
-- changes cluster.
local K = 3
local MAX_ITERATIONS = 100
-- Each point is a small table with named fields: {x = 1.25, y = 1.66}.
local function read_points(path)
local f = assert(io.open(path, "r"))
local header, xcol, ycol = f:read("l"), nil, nil
local i = 0
for field in header:gmatch("[^,]+") do
i = i + 1
if field == "x" then xcol = i elseif field == "y" then ycol = i end
end
local points = {}
for line in f:lines() do
local fields = {}
for field in line:gmatch("[^,]+") do fields[#fields + 1] = field end
points[#points + 1] = { x = tonumber(fields[xcol]), y = tonumber(fields[ycol]) }
end
f:close()
return points
end
local function distance(p, q)
return math.sqrt((p.x - q.x) ^ 2 + (p.y - q.y) ^ 2)
end
-- Returns the index of the nearest centroid. The test is strictly less
-- than, so a tie keeps the lower index.
local function nearest(point, centroids)
local best, best_dist = 1, distance(point, centroids[1])
for j = 2, #centroids do
local d = distance(point, centroids[j])
if d < best_dist then best, best_dist = j, d end
end
return best
end
-- == on two tables asks whether they are the same table, not whether they
-- hold the same values, so we compare the assignments element by element.
local function same(a, b)
if #a ~= #b then return false end
for i = 1, #a do
if a[i] ~= b[i] then return false end
end
return true
end
local points = read_points(arg[1])
local centroids = {}
for j = 1, K do
centroids[j] = { x = points[j].x, y = points[j].y }
end
local assignment = {}
local iterations = 0
while iterations < MAX_ITERATIONS do
iterations = iterations + 1
local new_assignment = {}
for i, p in ipairs(points) do
new_assignment[i] = nearest(p, centroids)
end
if same(new_assignment, assignment) then break end
assignment = new_assignment
for j = 1, K do
local sx, sy, n = 0, 0, 0
for i, p in ipairs(points) do
if assignment[i] == j then
sx, sy, n = sx + p.x, sy + p.y, n + 1
end
end
if n > 0 then -- an empty cluster keeps its old centroid
centroids[j] = { x = sx / n, y = sy / n }
end
end
end
local sizes, inertia = {}, 0
for j = 1, K do sizes[j] = 0 end
for i, p in ipairs(points) do
local a = assignment[i]
sizes[a] = sizes[a] + 1
inertia = inertia + distance(p, centroids[a]) ^ 2
end
print("iterations " .. iterations)
for j, c in ipairs(centroids) do
print(string.format("centroid_%d_x %.6f", j, c.x))
print(string.format("centroid_%d_y %.6f", j, c.y))
end
for j = 1, K do
print(string.format("size_%d %d", j, sizes[j]))
end
print(string.format("inertia %.6f", inertia))Each point is a small table with named fields, {x = 1.25, y = 1.66}, so
we write p.x and p.y. Lua has no distance function, so distance
writes out the square root of the summed squares, and nearest replaces
its best only on a strictly smaller distance, so a tie keeps the lower
index.
To notice that nothing changed, we compare the whole assignment with the
last one. Lua’s == can’t do that: on two tables it asks whether they are
the same table, so {1, 2} == {1, 2} is false. A short same function
compares them element by element instead. Lua has no clustering or
statistics library to avoid.
(* k-means clustering (Lloyd's algorithm), using only OCaml's standard
library.
Starts from the first k points as centroids, then alternates two steps
until no point changes cluster: assign each point to its nearest centroid,
and move each centroid to the mean of its points. The repetition is a
function that calls itself, the usual way to loop in OCaml. *)
(* A record with named fields; p.x reads one. *)
type point = { x : float; y : float }
let k = 3
let max_iterations = 100
(* The x and y columns as an array of points. The header gives their
positions. *)
let read_points path =
In_channel.with_open_text path (fun ic ->
match In_channel.input_lines ic with
| [] -> failwith "empty file"
| header :: rows ->
let names = String.split_on_char ',' header in
let index name =
match List.find_index (String.equal name) names with
| Some i -> i
| None -> failwith ("column not found: " ^ name)
in
let xi = index "x" and yi = index "y" in
rows
|> List.filter (fun row -> row <> "")
|> List.map (fun row ->
let fields = String.split_on_char ',' row in
{ x = float_of_string (List.nth fields xi);
y = float_of_string (List.nth fields yi) })
|> Array.of_list)
(* Squared Euclidean distance. It ranks centroids in the same order as the
distance itself, so finding the nearest one needs no sqrt. *)
let squared_distance a b = ((a.x -. b.x) ** 2.) +. ((a.y -. b.y) ** 2.)
(* The index of the closest centroid. Only a strictly smaller distance
replaces the best, so a tie goes to the lower index. *)
let nearest centroids p =
let best = ref 0 in (* ref makes a mutable cell: !best reads, := writes *)
Array.iteri
(fun j c ->
if squared_distance p c < squared_distance p centroids.(!best) then best := j)
centroids;
!best
(* Move each centroid to the mean of its points. An empty cluster keeps its
old centroid. *)
let update points assignment centroids =
Array.mapi
(fun j c ->
let sum_x = ref 0. and sum_y = ref 0. and count = ref 0 in
Array.iteri
(fun i p ->
if assignment.(i) = j then begin
sum_x := !sum_x +. p.x;
sum_y := !sum_y +. p.y;
incr count
end)
points;
if !count = 0 then c
else
let n = float_of_int !count in
{ x = !sum_x /. n; y = !sum_y /. n })
centroids
(* One assignment step, then either stop or update and go again. previous is
an option: None before the first pass, so the first comparison always
counts as a change. = compares two arrays element by element. Returns the
number of assignment steps, the final assignment, and the centroids. *)
let rec run points iterations previous centroids =
let step = iterations + 1 in
let assignment = Array.map (nearest centroids) points in
if previous = Some assignment then (step, assignment, centroids)
else
let moved = update points assignment centroids in
if step = max_iterations then (step, assignment, moved)
else run points step (Some assignment) moved
let () =
let points = read_points Sys.argv.(1) in
let iterations, assignment, centroids =
run points 0 None (Array.sub points 0 k)
in
Printf.printf "iterations %d\n" iterations;
(* Clusters are numbered from 0 inside the program and from 1 in the output. *)
Array.iteri
(fun j c -> Printf.printf "centroid_%d_x %.6f\ncentroid_%d_y %.6f\n" (j + 1) c.x (j + 1) c.y)
centroids;
for j = 0 to k - 1 do
let size = Array.fold_left (fun n a -> if a = j then n + 1 else n) 0 assignment in
Printf.printf "size_%d %d\n" (j + 1) size
done;
let inertia = ref 0. in
Array.iteri
(fun i p -> inertia := !inertia +. squared_distance p centroids.(assignment.(i)))
points;
Printf.printf "inertia %.6f\n" !inertiaEach point is a record, { x : float; y : float }, and the assignment is
an array of cluster numbers. OCaml has while loops, but the usual way to
repeat is recursion, so run calls itself with the moved centroids. Its
previous assignment is an option: None before the first pass, then
Some the last array. OCaml’s = compares arrays element by element, so
“nothing changed” is one previous = Some assignment.
The standard library has no k-means, so the means and the
nearest-centroid search are written out. hypot exists, but squared
distance ranks centroids the same way without a square root. A tie goes
to the lower index, since only a strictly smaller distance replaces the
best.
% k-means clustering (Lloyd's algorithm), using only core GNU Octave.
%
% Starts from the first K points as centroids, then alternates two steps until
% no point changes cluster: assign each point to its nearest centroid, and move
% each centroid to the mean of its points. Each point is a row of a 150-by-2
% matrix, so both steps work on every point at once.
% The named column of a CSV file, as a column vector of numbers.
function values = read_column(path, name)
lines = strsplit(strtrim(fileread(path)), "\n");
header = strsplit(strtrim(lines{1}), ",");
index = find(strcmp(header, name));
if isempty(index)
error("no column %s", name);
end
values = zeros(numel(lines) - 1, 1);
for row = 2:numel(lines)
fields = strsplit(strtrim(lines{row}), ",");
values(row - 1) = str2double(fields{index});
end
end
K = 3;
MAX_ITERATIONS = 100;
args = argv();
points = [read_column(args{1}, "x"), read_column(args{1}, "y")];
n = rows(points);
centroids = points(1:K, :); % a copy: changing it never changes points
assignment = [];
distances = zeros(n, K);
iterations = 0;
while iterations < MAX_ITERATIONS
iterations += 1;
% points - centroids(j, :) subtracts the 1-by-2 row from all 150 rows
% (broadcasting), and sum(..., 2) adds across each row, giving every
% point's squared distance to centroid j.
for j = 1:K
distances(:, j) = sum((points - centroids(j, :)) .^ 2, 2);
end
% The smallest squared distance is also the smallest distance. min with
% [] and 2 works along each row; its second output is the column of the
% minimum, the first one on a tie, so ties go to the lower index.
[~, new_assignment] = min(distances, [], 2);
if isequal(new_assignment, assignment)
break;
end
assignment = new_assignment;
for j = 1:K
members = points(assignment == j, :); % the rows assigned to j
if rows(members) > 0 % an empty cluster keeps its old centroid
centroids(j, :) = sum(members, 1) / rows(members);
end
end
end
% Pair each point with its own centroid's row, then add every squared difference.
inertia = sum(sum((points - centroids(assignment, :)) .^ 2));
printf("iterations %d\n", iterations);
for j = 1:K
printf("centroid_%d_x %.6f\n", j, centroids(j, 1));
printf("centroid_%d_y %.6f\n", j, centroids(j, 2));
end
for j = 1:K
printf("size_%d %d\n", j, sum(assignment == j));
end
printf("inertia %.6f\n", inertia);The points are rows of a 150-by-2 matrix, so each step works on all of them
at once. points - centroids(j, :) subtracts one centroid from every row,
and sum(..., 2) adds across each row to give squared distances. Then
min(distances, [], 2) returns each row’s nearest column, the first one on a
tie, so ties go to the lower-numbered centroid.
isequal compares the new assignment with the last, and it checks shape as
well as values. A new centroid is sum(members, 1) / rows(members), written
out instead of calling mean. Core Octave has no kmeans, and the
statistics package isn’t installed.
{ k-means clustering (Lloyd's algorithm), using only Free Pascal's run-time
library.
Starts from the first K points as centroids, then alternates two steps
until no point changes cluster: assign each point to its nearest centroid,
and move each centroid to the mean of its points. None of the installed
units has a k-means routine, so there is no built-in to avoid. }
program KMeans;
{$mode objfpc}{$H+} // 32-bit Integer, and strings longer than 255 characters
uses
SysUtils;
const
K = 3;
MaxIterations = 100;
type
// A record is a value: assigning one copies both fields.
TPoint = record
X, Y: Double;
end;
TPointArray = array of TPoint;
TCentroids = array[0..K - 1] of TPoint; // a fixed-size array, indexed 0..2
function ReadPoints(const Path: String): TPointArray;
var
F: TextFile;
Line: String;
Fields: TStringArray;
XCol, YCol, I: Integer;
begin
Result := nil;
AssignFile(F, Path);
Reset(F);
ReadLn(F, Line);
Fields := Line.Split(',');
XCol := -1;
YCol := -1;
for I := 0 to High(Fields) do
begin
if Fields[I] = 'x' then XCol := I;
if Fields[I] = 'y' then YCol := I;
end;
if (XCol < 0) or (YCol < 0) then
raise Exception.Create('no x or y column');
while not Eof(F) do
begin
ReadLn(F, Line);
if Line = '' then
Continue;
Fields := Line.Split(',');
SetLength(Result, Length(Result) + 1); // grow by one
Result[High(Result)].X := StrToFloat(Fields[XCol]);
Result[High(Result)].Y := StrToFloat(Fields[YCol]);
end;
CloseFile(F);
end;
{ Squared Euclidean distance. It ranks centroids in the same order as the
distance itself, so the nearest one needs no Sqrt. }
function SquaredDistance(const A, B: TPoint): Double;
begin
Result := Sqr(A.X - B.X) + Sqr(A.Y - B.Y);
end;
{ The index of the closest centroid. Only a strictly smaller distance
replaces the best so far, so a tie goes to the lower index. }
function Nearest(const P: TPoint; const Centroids: TCentroids): Integer;
var
J: Integer;
begin
Result := 0;
for J := 1 to K - 1 do
if SquaredDistance(P, Centroids[J]) < SquaredDistance(P, Centroids[Result]) then
Result := J;
end;
procedure Say(const Name: String; Value: Double);
begin
WriteLn(Name, ' ', Value:0:6); // width 0, six decimals
end;
var
Points: TPointArray;
Centroids, Sums: TCentroids;
Assignment: array of Integer;
Sizes: array[0..K - 1] of Integer;
N, I, J, Iterations: Integer;
Changed: Boolean;
Inertia: Double;
begin
Points := ReadPoints(ParamStr(1));
N := Length(Points);
if N < K then
raise Exception.Create('need at least 3 points');
for J := 0 to K - 1 do
Centroids[J] := Points[J]; // a copy, so moving a centroid leaves the point alone
// -1 means "no cluster yet", so the first pass always counts as a change.
SetLength(Assignment, N);
for I := 0 to N - 1 do
Assignment[I] := -1;
Iterations := 0;
while Iterations < MaxIterations do
begin
Inc(Iterations);
Changed := False;
for I := 0 to N - 1 do
begin
J := Nearest(Points[I], Centroids);
if J <> Assignment[I] then
begin
Assignment[I] := J;
Changed := True;
end;
end;
if not Changed then
Break;
// Sum each cluster's points in one pass, then divide by its size.
for J := 0 to K - 1 do
begin
Sums[J].X := 0;
Sums[J].Y := 0;
Sizes[J] := 0;
end;
for I := 0 to N - 1 do
begin
J := Assignment[I];
Sums[J].X := Sums[J].X + Points[I].X;
Sums[J].Y := Sums[J].Y + Points[I].Y;
Inc(Sizes[J]);
end;
for J := 0 to K - 1 do
if Sizes[J] > 0 then // an empty cluster keeps its old centroid
begin
Centroids[J].X := Sums[J].X / Sizes[J];
Centroids[J].Y := Sums[J].Y / Sizes[J];
end;
end;
for J := 0 to K - 1 do
Sizes[J] := 0;
Inertia := 0;
for I := 0 to N - 1 do
begin
Inc(Sizes[Assignment[I]]);
Inertia := Inertia + SquaredDistance(Points[I], Centroids[Assignment[I]]);
end;
// Clusters are numbered from 0 inside the program and from 1 in the output.
WriteLn('iterations ', Iterations);
for J := 0 to K - 1 do
begin
Say('centroid_' + IntToStr(J + 1) + '_x', Centroids[J].X);
Say('centroid_' + IntToStr(J + 1) + '_y', Centroids[J].Y);
end;
for J := 0 to K - 1 do
WriteLn('size_', J + 1, ' ', Sizes[J]);
Say('inertia', Inertia);
end.Each point is a TPoint record with X and Y fields, held in one
dynamic array, and the three centroids sit in a fixed-size
array[0..K - 1]. A record is a value, so Centroids[J] := Points[J]
copies it, and moving a centroid can’t change a data point.
We notice that nothing changed with a flag: every assignment starts at −1,
and the pass sets Changed whenever a point’s nearest centroid differs
from its stored one. Nearest compares squared distances, which rank
centroids the same way, so it needs no Sqrt, and only a strictly smaller
distance wins, so ties go to the lower index. None of the installed units
has a k-means routine to avoid.
# k-means clustering (Lloyd's algorithm), using only Perl and its core
# modules.
#
# Starts from the first K points as centroids, then alternates two steps
# until no point changes cluster: assign each point to its nearest centroid,
# and move each centroid to the mean of its points.
use strict;
use warnings;
use List::Util qw(sum);
my $K = 3;
my $MAX_ITERATIONS = 100;
# Each point is a reference to a two-element array, [x, y]; $p->[0] is x.
sub read_points {
my ($path) = @_;
open my $fh, '<', $path or die "can't open $path: $!";
chomp(my $header = <$fh>);
my @fields = split /,/, $header;
my ($ix) = grep { $fields[$_] eq 'x' } 0 .. $#fields;
my ($iy) = grep { $fields[$_] eq 'y' } 0 .. $#fields;
die "no x and y columns" unless defined $ix && defined $iy;
my @points;
while (my $line = <$fh>) {
chomp $line;
my @row = split /,/, $line;
push @points, [$row[$ix], $row[$iy]];
}
close $fh;
return @points;
}
sub squared_distance {
my ($p, $q) = @_;
return ($p->[0] - $q->[0]) ** 2 + ($p->[1] - $q->[1]) ** 2;
}
# Only a strictly smaller distance replaces the best, so ties go to the
# lower index.
sub nearest {
my ($point, @centroids) = @_;
my $best = 0;
for my $j (1 .. $#centroids) {
$best = $j
if squared_distance($point, $centroids[$j])
< squared_distance($point, $centroids[$best]);
}
return $best;
}
my @points = read_points($ARGV[0]);
# Copy the first K points so moving a centroid never changes a point.
my @centroids = map { [@$_] } @points[0 .. $K - 1];
my @assignment;
my $iterations = 0;
while ($iterations < $MAX_ITERATIONS) {
$iterations++;
my @new_assignment = map { nearest($_, @centroids) } @points;
# Arrays have no ==, so join each into a string to compare them whole.
last if "@new_assignment" eq "@assignment"; # this pass still counts
@assignment = @new_assignment;
for my $j (0 .. $K - 1) {
my @members = map { $points[$_] }
grep { $assignment[$_] == $j } 0 .. $#points;
next unless @members; # an empty cluster keeps its old centroid
$centroids[$j] = [
sum(map { $_->[0] } @members) / @members,
sum(map { $_->[1] } @members) / @members,
];
}
}
my @sizes = map { my $j = $_; scalar grep { $_ == $j } @assignment } 0 .. $K - 1;
my $inertia = sum(map { squared_distance($points[$_], $centroids[$assignment[$_]]) }
0 .. $#points);
print "iterations $iterations\n";
for my $j (0 .. $K - 1) {
printf "centroid_%d_x %.6f\n", $j + 1, $centroids[$j][0];
printf "centroid_%d_y %.6f\n", $j + 1, $centroids[$j][1];
}
printf "size_%d %d\n", $_ + 1, $sizes[$_] for 0 .. $K - 1;
printf "inertia %.6f\n", $inertia;Each point is a reference to a two-element array, so $p->[0] is its x,
and the centroids start as copies of the first three points, so moving a
centroid never moves a point. A loop keeps the first nearest centroid, and
only a strictly smaller distance replaces it, so ties go to the lower
index.
Perl arrays have no ==, so we compare whole assignments as strings:
putting an array inside quotes, as in "@assignment", joins its cluster
numbers with spaces, and eq compares the two strings. sum comes from List::Util, a core module. Perl has no clustering
built in, so there’s nothing to skip.
<?php
// k-means clustering (Lloyd's algorithm), using only PHP and the extensions
// built into it.
//
// Starts from the first K points as centroids, then alternates two steps
// until no point changes cluster: assign each point to its nearest centroid,
// and move each centroid to the mean of its points.
declare(strict_types=1);
const K = 3;
const MAX_ITERATIONS = 100;
// Each point is a two-element array [x, y].
function read_points(string $path): array
{
$lines = file($path, FILE_IGNORE_NEW_LINES | FILE_SKIP_EMPTY_LINES);
array_shift($lines); // drop the header row "x,y"
$points = [];
foreach ($lines as $line) {
[$x, $y] = explode(',', $line); // list-style destructuring
$points[] = [(float) $x, (float) $y];
}
return $points;
}
function squared_distance(array $p, array $q): float
{
return ($p[0] - $q[0]) ** 2 + ($p[1] - $q[1]) ** 2;
}
// Only a strictly smaller distance replaces the best so far, so a tie goes
// to the lower index. Squared distances rank the same as distances.
function nearest(array $point, array $centroids): int
{
$best = 0;
for ($j = 1; $j < count($centroids); $j++) {
if (squared_distance($point, $centroids[$j]) < squared_distance($point, $centroids[$best])) {
$best = $j;
}
}
return $best;
}
$points = read_points($argv[1]);
// Assigning an array copies it, so moving a centroid never moves a point.
$centroids = array_slice($points, 0, K);
$assignment = null;
$iterations = 0;
while ($iterations < MAX_ITERATIONS) {
$iterations++;
$new = array_map(fn($p) => nearest($p, $centroids), $points);
// === is true only when both arrays hold the same keys and values, in
// the same order, with the same types.
if ($new === $assignment) {
break;
}
$assignment = $new;
for ($j = 0; $j < K; $j++) {
$members = array_values(array_filter(
$points,
fn($i) => $assignment[$i] === $j,
ARRAY_FILTER_USE_KEY // the callback gets each point's index
));
$n = count($members);
if ($n > 0) { // an empty cluster keeps its old centroid
$centroids[$j] = [
array_sum(array_column($members, 0)) / $n,
array_sum(array_column($members, 1)) / $n,
];
}
}
}
$sizes = array_fill(0, K, 0);
$inertia = 0.0;
foreach ($points as $i => $p) {
$sizes[$assignment[$i]]++;
$inertia += squared_distance($p, $centroids[$assignment[$i]]);
}
echo "iterations $iterations\n";
foreach ($centroids as $j => [$x, $y]) {
printf("centroid_%d_x %.6f\n", $j + 1, $x);
printf("centroid_%d_y %.6f\n", $j + 1, $y);
}
foreach ($sizes as $j => $size) {
printf("size_%d %d\n", $j + 1, $size);
}
printf("inertia %.6f\n", $inertia);Each point is a two-element array, so $p[0] is its x. PHP arrays are
values, not references: assigning one makes a copy, so the centroids start
as copies of the first three points, and moving a centroid never moves a
point. A loop keeps the first nearest centroid, and only a strictly smaller
distance replaces it, so ties go to the lower index.
In PHP, === compares two whole arrays, key by key and
value by value, so one comparison tells you the assignments stopped
changing. array_column pulls out every x or every y for the means. PHP has
no clustering built in, so there’s nothing to skip.
-- k-means clustering (Lloyd's algorithm, k = 3), in PostgreSQL.
--
-- Starts from the first 3 points as centroids, then alternates two steps until
-- no point changes cluster: assign each point to its nearest centroid, and move
-- each centroid to the mean of its points. The loop is a WITH RECURSIVE query
-- whose state is one row of arrays: the centroids and every point's cluster.
WITH RECURSIVE points AS (
-- A table has no row order; run.sh numbers rows in file order as row_id.
SELECT row_id AS id, x, y FROM data
),
steps (iteration, cx, cy, assignment, done) AS (
-- Iteration 0: the first 3 points are the centroids; nothing is assigned yet.
SELECT 0,
array_agg(x ORDER BY id),
array_agg(y ORDER BY id),
NULL::int[],
false
FROM points
WHERE id <= 3
UNION ALL
SELECT s.iteration + 1, moved.cx, moved.cy, assigned.assignment,
assigned.assignment IS NOT DISTINCT FROM s.assignment
FROM steps s,
LATERAL (
-- Assign: each point's nearest centroid (ties go to the lower index).
SELECT array_agg(nearest.j ORDER BY p.id) AS assignment
FROM points p,
LATERAL (
SELECT j FROM generate_series(1, 3) AS j
ORDER BY sqrt((p.x - s.cx[j]) ^ 2 + (p.y - s.cy[j]) ^ 2), j
LIMIT 1
) AS nearest
) AS assigned,
LATERAL (
-- Update: each centroid moves to the mean of its points; an empty
-- cluster keeps its old centroid.
SELECT array_agg(coalesce(sum_x / n, s.cx[j]) ORDER BY j) AS cx,
array_agg(coalesce(sum_y / n, s.cy[j]) ORDER BY j) AS cy
FROM (
SELECT j, sum(p.x) AS sum_x, sum(p.y) AS sum_y, nullif(count(p.id), 0) AS n
FROM generate_series(1, 3) AS j
LEFT JOIN points p ON assigned.assignment[p.id] = j
GROUP BY j
) AS clusters
) AS moved
-- Stop once an assignment pass changed nothing (or after 100 passes).
WHERE NOT s.done AND s.iteration < 100
),
final AS (
SELECT * FROM steps ORDER BY iteration DESC LIMIT 1
),
members AS (
-- Each point with its final cluster.
SELECT p.x, p.y, a.j, f.cx[a.j] AS cx, f.cy[a.j] AS cy
FROM final f,
unnest(f.assignment) WITH ORDINALITY AS a (j, id)
JOIN points p ON p.id = a.id
),
sizes AS (
-- count(m.j) is 0 for an empty cluster.
SELECT j, count(m.j) AS size
FROM generate_series(1, 3) AS j
LEFT JOIN members m USING (j)
GROUP BY j
),
inertia AS (
SELECT sum((x - cx) ^ 2 + (y - cy) ^ 2) AS inertia FROM members
)
SELECT name, value
FROM (
SELECT 0 AS position, 'iterations' AS name, iteration::text AS value FROM final
UNION ALL
SELECT 2 * j - 1, format('centroid_%s_x', j), round(f.cx[j]::numeric, 6)::text
FROM final f, generate_series(1, 3) AS j
UNION ALL
SELECT 2 * j, format('centroid_%s_y', j), round(f.cy[j]::numeric, 6)::text
FROM final f, generate_series(1, 3) AS j
UNION ALL
SELECT 6 + j, format('size_%s', j), size::text FROM sizes
UNION ALL
SELECT 10, 'inertia', round(inertia::numeric, 6)::text FROM inertia
) AS results
ORDER BY position;Each half of an iteration is natural SQL. “Nearest centroid” is ORDER BY
distance LIMIT 1 over three candidates, and “move to the mean” is a
GROUP BY on the cluster with sum and count.
The loop is the hard part. SQL has none, so we use WITH RECURSIVE, where
each step is one row holding the whole state as arrays: centroid x values,
centroid y values, and every point’s cluster. PostgreSQL rejects aggregates
in the recursive step, so the assign and update steps sit in LATERAL
subqueries. “The first three points” also needs file order, which a table
doesn’t have; it comes from the row_id column the harness adds while
loading. Core PostgreSQL has no k-means function to skip.
% k-means clustering (Lloyd's algorithm), using SWI-Prolog and the
% libraries that come with it.
%
% Starts from the first K points as centroids, then alternates two steps
% until no point changes cluster: assign each point to its nearest
% centroid, and move each centroid to the mean of its points. Prolog has
% no loops, so the repetition is a rule, run/7, that calls itself with the
% moved centroids.
:- initialization(main, main).
k(3).
max_iterations(100).
% read_column(+Path, +Name, -Values): Values are the numbers in the column
% called Name. findall collects every Value that satisfies the goals in
% the middle, one per line.
read_column(Path, Name, Values) :-
read_file_to_string(Path, Text, []),
split_string(Text, "\n", "", [HeaderLine|Lines]),
split_string(HeaderLine, ",", "", Header),
nth0(Index, Header, Name),
findall(Value,
( member(Line, Lines),
Line \== "",
split_string(Line, ",", "", Fields),
nth0(Index, Fields, Field),
number_string(Value, Field) ),
Values).
% A point is the term X-Y. The - only joins the two numbers into a pair;
% it subtracts nothing, because only "is" does arithmetic.
% squared_distance(+P, +Q, -D): D is the squared Euclidean distance. It
% ranks centroids in the same order as the distance, so no sqrt is needed.
squared_distance(X1-Y1, X2-Y2, D) :-
D is (X1 - X2) ** 2 + (Y1 - Y2) ** 2.
% nearest(+Centroids, +Point, -J): J is the number of the closest
% centroid, counting from 1. aggregate_all(min(D, I), Goal, min(_, J))
% finds the smallest D over every solution of Goal and the I that went
% with it. nth1 tries the centroids in order, and when several share the
% smallest D, aggregate_all keeps the first: a tie goes to the lower index.
nearest(Centroids, Point, J) :-
aggregate_all(min(D, I),
( nth1(I, Centroids, C),
squared_distance(Point, C, D) ),
min(_, J)).
% mean(+Values, -Mean)
mean(Values, Mean) :-
sum_list(Values, Sum),
length(Values, N),
Mean is Sum / N.
% update(+Points, +Assignment, +Centroids, -Moved): Moved holds each
% centroid moved to the mean of its points. pairs_keys_values joins each
% point to its cluster number, as J-(X-Y).
update(Points, Assignment, Centroids, Moved) :-
pairs_keys_values(Labeled, Assignment, Points),
findall(New,
( nth1(J, Centroids, Old),
findall(P, member(J-P, Labeled), Members),
move(Members, Old, New) ),
Moved).
% move(+Members, +Old, -New): an empty cluster keeps its old centroid.
move([], Old, Old).
move([P|Ps], _, X-Y) :-
pairs_keys_values([P|Ps], Xs, Ys),
mean(Xs, X),
mean(Ys, Y).
% run(+Points, +Done, +Previous, +Centroids, -Iterations, -Assignment,
% -Final): one assignment step, then either stop or update and go
% again. Done counts the steps so far. Previous is the last assignment, or
% the atom none before the first pass, so the first comparison always
% counts as a change. == asks whether two terms are identical, so it
% compares the whole assignment list at once.
run(Points, Done, Previous, Centroids, Iterations, Assignment, Final) :-
Step is Done + 1,
maplist(nearest(Centroids), Points, New),
( New == Previous
-> Iterations = Step, Assignment = New, Final = Centroids
; update(Points, New, Centroids, Moved),
( max_iterations(Step)
-> Iterations = Step, Assignment = New, Final = Moved
; run(Points, Step, New, Moved, Iterations, Assignment, Final)
)
).
% point_inertia(+Centroids, +Point, +J, -D): the squared distance from
% Point to its centroid, number J.
point_inertia(Centroids, Point, J, D) :-
nth1(J, Centroids, C),
squared_distance(Point, C, D).
main :-
current_prolog_flag(argv, [Path|_]),
read_column(Path, "x", Xs),
read_column(Path, "y", Ys),
pairs_keys_values(Points, Xs, Ys),
k(K),
length(Start, K),
append(Start, _, Points),
run(Points, 0, none, Start, Iterations, Assignment, Centroids),
maplist(point_inertia(Centroids), Points, Assignment, Squares),
sum_list(Squares, Inertia),
format("iterations ~d~n", [Iterations]),
forall(nth1(J, Centroids, X-Y),
format("centroid_~d_x ~6f~ncentroid_~d_y ~6f~n", [J, X, J, Y])),
% aggregate_all(count, ...) counts the solutions of a goal: here, the
% positions in Assignment that hold J.
forall(between(1, K, J),
( aggregate_all(count, member(J, Assignment), Size),
format("size_~d ~d~n", [J, Size]) )),
format("inertia ~6f~n", [Inertia]).Each point is the term X-Y. The - only pairs the two numbers, since
only is does arithmetic, and pairs_keys_values splits points back
into their x and y lists. Prolog has no loops, so run repeats by calling
itself with the moved centroids. Its previous assignment starts as the
atom none, and “nothing changed” is one == between whole lists.
nearest asks aggregate_all(min(D, I), ...) for the smallest squared
distance and the centroid number that went with it. When several
solutions share the minimum, it keeps the first, so a tie goes to the
lower index. There’s no clustering, mean, or distance built in to skip.
# k-means clustering (Lloyd's algorithm), using only base R.
#
# Starts from the first k points as centroids, then alternates two steps until
# no point changes cluster: assign each point to its nearest centroid, and move
# each centroid to the mean of its points. kmeans() is built in; it is avoided
# here so the loop stays visible.
k <- 3
max_iterations <- 100
args <- commandArgs(trailingOnly = TRUE)
points <- as.matrix(read.csv(args[1])[, c("x", "y")])
# One centroid per row; rows 1..k of the data are the starting centroids.
centroids <- points[1:k, , drop = FALSE]
nearest <- function(point) {
distance <- (centroids[, 1] - point[1])^2 + (centroids[, 2] - point[2])^2
which.min(distance) # first minimum, so ties go to the lower index
}
assignment <- NULL
iterations <- 0
while (iterations < max_iterations) {
iterations <- iterations + 1
new_assignment <- apply(points, 1, nearest)
if (identical(new_assignment, assignment)) break
assignment <- new_assignment
for (j in 1:k) {
members <- points[assignment == j, , drop = FALSE]
if (nrow(members) > 0) { # an empty cluster keeps its old centroid
centroids[j, ] <- colSums(members) / nrow(members)
}
}
}
sizes <- tabulate(assignment, nbins = k)
# centroids[assignment, ] lines up each point with its own centroid.
inertia <- sum((points - centroids[assignment, , drop = FALSE])^2)
cat(sprintf("iterations %d\n", iterations))
for (j in 1:k) {
cat(sprintf("centroid_%d_x %.6f\n", j, centroids[j, 1]))
cat(sprintf("centroid_%d_y %.6f\n", j, centroids[j, 2]))
}
cat(sprintf("size_%d %d\n", 1:k, sizes), sep = "")
cat(sprintf("inertia %.6f\n", inertia))Matrices fit k-means well. The points are a two-column matrix,
colSums(members) / nrow(members) moves a centroid, and
centroids[assignment, ] lines each point up with its own centroid, so the
inertia is one sum(). which.min() returns the first minimum, which is the
lower-index tie rule, and identical() checks whether anything changed.
The trap is that subsetting a matrix down to one row silently turns it into
a plain vector unless you pass drop = FALSE.
kmeans() is built in, but its default algorithm is Hartigan-Wong, not the
loop on this page. kmeans(points, centers = points[1:3, ], algorithm = "Lloyd") reproduces our result exactly: 5 iterations, sizes 52, 50, and 48,
and inertia 423.627610.
# k-means clustering with kmeans(), R's built-in function.
#
# The built-in version of ../r/main.R. kmeans() defaults to the Hartigan-Wong
# algorithm and to random starting centroids, so we ask for Lloyd's algorithm
# and pass the first k points as the centers, to match the from-scratch steps.
k <- 3
args <- commandArgs(trailingOnly = TRUE)
points <- as.matrix(read.csv(args[1])[, c("x", "y")])
fit <- kmeans(points, centers = points[1:k, ], algorithm = "Lloyd", iter.max = 100)
# fit$iter counts the assignment passes, including the last one that changed
# nothing, and tot.withinss is the inertia.
cat(sprintf("iterations %d\n", fit$iter))
for (j in 1:k) {
cat(sprintf("centroid_%d_x %.6f\n", j, fit$centers[j, 1]))
cat(sprintf("centroid_%d_y %.6f\n", j, fit$centers[j, 2]))
}
cat(sprintf("size_%d %d\n", 1:k, fit$size), sep = "")
cat(sprintf("inertia %.6f\n", fit$tot.withinss))kmeans() is part of base R’s stats package, but its defaults don’t run
the method on this page. It picks random starting centroids, and it uses
the Hartigan-Wong algorithm, which moves single points between clusters
whenever that lowers the total.
So we pass centers = points[1:k, ] to start from the first three points
and algorithm = "Lloyd" for the assign-then-update loop, and the results
match exactly. fit$iter counts the passes the same way, including the
last one, and tot.withinss, the total within-cluster sum of squares, is
the inertia.
# k-means clustering (Lloyd's algorithm), using only the standard library.
# Starts from the first K points as centroids, then alternates two steps until
# no point changes cluster: assign each point to its nearest centroid, and move
# each centroid to the mean of its points.
K = 3
MAX_ITERATIONS = 100
# The csv library is a bundled gem since Ruby 3.4, so split lines by hand.
# Fine for this data: no quoted fields or embedded commas.
def read_points(path)
header, *rows = File.readlines(path, chomp: true).map { |line| line.split(",") }
ix = header.index("x") or abort "no column x"
iy = header.index("y") or abort "no column y"
rows.map { |row| [Float(row[ix]), Float(row[iy])] }
end
def squared_distance(a, b)
a.zip(b).sum { |p, q| (p - q)**2 }
end
# min_by returns the first of equal values, so ties go to the lower index.
def nearest(point, centroids)
(0...centroids.length).min_by { |j| squared_distance(point, centroids[j]) }
end
points = read_points(ARGV[0])
centroids = points.first(K)
assignment = nil
iterations = 0
while iterations < MAX_ITERATIONS
iterations += 1
new_assignment = points.map { |p| nearest(p, centroids) }
break if new_assignment == assignment # this unchanged pass still counts
assignment = new_assignment
K.times do |j|
members = points.select.with_index { |_, i| assignment[i] == j }
next if members.empty? # an empty cluster keeps its old centroid
centroids[j] = members.transpose.map { |coords| coords.sum / members.length }
end
end
sizes = (0...K).map { |j| assignment.count(j) }
inertia = points.zip(assignment).sum { |p, a| squared_distance(p, centroids[a]) }
puts "iterations #{iterations}"
centroids.each.with_index(1) do |(x, y), j|
puts format("centroid_%d_x %.6f", j, x)
puts format("centroid_%d_y %.6f", j, y)
end
sizes.each.with_index(1) { |size, j| puts "size_#{j} #{size}" }
puts format("inertia %.6f", inertia)Ruby reads close to the description of the algorithm. Points are
two-element arrays, min_by picks the nearest centroid and returns the
first of equal values, and transpose turns a cluster’s points into an x
column and a y column, so the new centroid is one map. Array ==
compares element by element, which makes “did anything change?” a single
comparison.
Ruby pushed back only on reading the file. csv is a bundled gem rather
than part of the standard library, so we split each line on commas by
hand, which is fine for this file. The standard library has no clustering
to skip.
//! k-means clustering (Lloyd's algorithm), using only the standard library.
//!
//! Starts from the first K points as centroids, then alternates two steps until
//! no point changes cluster: assign each point to its nearest centroid, and move
//! each centroid to the mean of its points.
use std::env;
use std::fs;
const K: usize = 3;
const MAX_ITERATIONS: usize = 100;
type Point = [f64; 2];
/// Reads the x and y columns by header name. There is no CSV parser in the
/// standard library; the file has no quoted fields, so commas suffice.
fn read_points(path: &str) -> Vec<Point> {
let text = fs::read_to_string(path).expect("cannot read data file");
let mut lines = text.lines();
let header: Vec<&str> = lines.next().expect("empty file").split(',').collect();
let column = |name| {
header
.iter()
.position(|&h| h == name)
.expect("column not found")
};
let (ix, iy) = (column("x"), column("y"));
lines
.map(|line| {
let fields: Vec<&str> = line.split(',').collect();
[fields[ix].parse().unwrap(), fields[iy].parse().unwrap()]
})
.collect()
}
fn squared_distance(a: &Point, b: &Point) -> f64 {
(a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2)
}
fn nearest(point: &Point, centroids: &[Point]) -> usize {
// A strict < keeps the first of equal distances, so ties go to the lower index.
let mut best = 0;
for j in 1..centroids.len() {
if squared_distance(point, ¢roids[j]) < squared_distance(point, ¢roids[best]) {
best = j;
}
}
best
}
fn main() {
let path = env::args().nth(1).expect("usage: main <data.csv>");
let points = read_points(&path);
let mut centroids: Vec<Point> = points[..K].to_vec();
// Option<Vec> stands for "no assignment yet" before the first pass.
let mut assignment: Option<Vec<usize>> = None;
let mut iterations = 0;
while iterations < MAX_ITERATIONS {
iterations += 1;
let new_assignment: Vec<usize> = points.iter().map(|p| nearest(p, ¢roids)).collect();
if assignment.as_ref() == Some(&new_assignment) {
break; // this unchanged pass still counts
}
for (j, centroid) in centroids.iter_mut().enumerate() {
let members: Vec<&Point> = points
.iter()
.zip(&new_assignment)
.filter(|&(_, &a)| a == j)
.map(|(p, _)| p)
.collect();
if !members.is_empty() {
// An empty cluster keeps its old centroid.
let n = members.len() as f64;
*centroid = [
members.iter().map(|p| p[0]).sum::<f64>() / n,
members.iter().map(|p| p[1]).sum::<f64>() / n,
];
}
}
assignment = Some(new_assignment);
}
let assignment = assignment.expect("at least one pass");
let inertia: f64 = points
.iter()
.zip(&assignment)
.map(|(p, &a)| squared_distance(p, ¢roids[a]))
.sum();
println!("iterations {iterations}");
for (j, [x, y]) in centroids.iter().enumerate() {
println!("centroid_{}_x {x:.6}", j + 1);
println!("centroid_{}_y {y:.6}", j + 1);
}
for j in 0..K {
let size = assignment.iter().filter(|&&a| a == j).count();
println!("size_{} {size}", j + 1);
}
println!("inertia {inertia:.6}");
}Rust gives each piece a precise type. A point is [f64; 2], a fixed-size
array, and “no assignment yet” is None rather than a sentinel value, so
the first pass can’t be mistaken for a repeat. Two assignments are compared
with == on vectors, which checks every element.
Nearest-centroid search is a plain loop. Only a strictly smaller distance replaces the best, so a tie goes to the lower index. We compare squared distances, which needs no square root, and inertia uses them directly. An empty cluster keeps its old centroid. The standard library has no clustering, so there was no built-in to skip.
// k-means clustering (Lloyd's algorithm), using only Swift's standard library
// and Foundation, which ships with it.
//
// Starts from the first k points as centroids, then alternates two steps until
// no point changes cluster: assign each point to its nearest centroid, and
// move each centroid to the mean of its points.
import Foundation
let k = 3
let maxIterations = 100
// A struct is a value type: copying a Point copies both numbers.
struct Point {
var x: Double
var y: Double
}
// The named column of a CSV file, as numbers. A function that can fail is
// marked throws, and every call to it is marked try.
func readColumn(_ path: String, _ name: String) throws -> [Double] {
let text = try String(contentsOfFile: path, encoding: .utf8)
// split skips empty pieces, so a blank last line doesn't become a row.
var lines = text.split(separator: "\n")
let header = lines.removeFirst().split(separator: ",").map(String.init)
guard let index = header.firstIndex(of: name) else {
fatalError("no column \(name)")
}
// Double(_:) gives nil for text that isn't a number, and ! stops the
// program there instead of carrying on with a missing value.
return lines.map { line in
Double(line.split(separator: ",", omittingEmptySubsequences: false)[index])!
}
}
func squaredDistance(_ a: Point, _ b: Point) -> Double {
(a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)
}
// The index of the nearest centroid. A strict < keeps the first of equal
// distances, so a tie goes to the lower index.
func nearest(_ point: Point, _ centroids: [Point]) -> Int {
var best = 0
for j in 1..<centroids.count
where squaredDistance(point, centroids[j]) < squaredDistance(point, centroids[best]) {
best = j
}
return best
}
let path = CommandLine.arguments[1]
let xs = try readColumn(path, "x")
let ys = try readColumn(path, "y")
// zip pairs the two columns row by row.
let points = zip(xs, ys).map { Point(x: $0, y: $1) }
var centroids = Array(points.prefix(k))
// An empty array stands for "no assignment yet"; a real one has 150 entries.
var assignment: [Int] = []
var iterations = 0
while iterations < maxIterations {
iterations += 1
let newAssignment = points.map { nearest($0, centroids) }
// == on arrays compares every element, so this asks whether any point moved.
if newAssignment == assignment {
break // this unchanged pass still counts
}
for j in 0..<k {
let members = zip(points, newAssignment).filter { $0.1 == j }.map { $0.0 }
// An empty cluster keeps its old centroid.
if !members.isEmpty {
let n = Double(members.count)
centroids[j] = Point(x: members.map(\.x).reduce(0, +) / n,
y: members.map(\.y).reduce(0, +) / n)
}
}
assignment = newAssignment
}
let inertia = zip(points, assignment)
.map { squaredDistance($0.0, centroids[$0.1]) }
.reduce(0, +)
func show(_ name: String, _ value: Double) {
print(name, String(format: "%.6f", value))
}
print("iterations", iterations)
for (j, c) in centroids.enumerated() {
show("centroid_\(j + 1)_x", c.x)
show("centroid_\(j + 1)_y", c.y)
}
for j in 0..<k {
print("size_\(j + 1)", assignment.filter { $0 == j }.count)
}
show("inertia", inertia)A point is a small Point struct with x and y fields, and zip
pairs the two columns into points. Swift arrays compare with ==
element by element, so we notice that nothing changed by comparing the
whole new assignment with the last one. Before the first pass the last
assignment is an empty array, which never equals one with 150 entries.
Nearest-centroid search is a plain loop. Only a strictly smaller distance replaces the best, so a tie goes to the lower index. We compare squared distances, which needs no square root, and inertia uses them directly. An empty cluster keeps its old centroid. Swift has no clustering, so there was no built-in to skip.
# k-means clustering (Lloyd's algorithm), using only Tcl's built-in commands.
#
# Tcl has no clustering (nor any matrix or statistics package; those live in
# Tcllib, which isn't installed), so both steps are written out. Starting
# from the first K points as centroids, we assign each point to its nearest
# centroid, move each centroid to the mean of its points, and stop when no
# point changes cluster.
set K 3
set MAX_ITERATIONS 100
# Each point is a two-element list {x y}; lassign unpacks it.
proc read_points {path} {
set f [open $path r]
set header [split [gets $f] ,]
set xcol [lsearch -exact $header x]
set ycol [lsearch -exact $header y]
set points {}
while {[gets $f line] >= 0} {
set fields [split $line ,]
lappend points [list [lindex $fields $xcol] [lindex $fields $ycol]]
}
close $f
return $points
}
proc distance {p q} {
lassign $p px py
lassign $q qx qy
return [expr {hypot($px - $qx, $py - $qy)}]
}
# Returns the index of the nearest centroid. The test is strictly less
# than, so a tie keeps the lower index.
proc nearest {point centroids} {
set best 0
set best_dist [distance $point [lindex $centroids 0]]
for {set j 1} {$j < [llength $centroids]} {incr j} {
set d [distance $point [lindex $centroids $j]]
if {$d < $best_dist} {
set best $j
set best_dist $d
}
}
return $best
}
proc say {name value} {
puts "$name [format %.6f $value]"
}
set points [read_points [lindex $argv 0]]
set centroids [lrange $points 0 $K-1]
set assignment {}
set iterations 0
while {$iterations < $MAX_ITERATIONS} {
incr iterations
set new_assignment {}
foreach p $points {
lappend new_assignment [nearest $p $centroids]
}
# A list is a string, so eq compares the whole assignment at once.
if {$new_assignment eq $assignment} break
set assignment $new_assignment
for {set j 0} {$j < $K} {incr j} {
set sx 0.0
set sy 0.0
set n 0
# foreach can walk two lists side by side.
foreach p $points a $assignment {
if {$a == $j} {
lassign $p x y
set sx [expr {$sx + $x}]
set sy [expr {$sy + $y}]
incr n
}
}
if {$n > 0} { ;# an empty cluster keeps its old centroid
lset centroids $j [list [expr {$sx / $n}] [expr {$sy / $n}]]
}
}
}
set inertia 0.0
foreach p $points a $assignment {
set inertia [expr {$inertia + [distance $p [lindex $centroids $a]] ** 2}]
}
puts "iterations $iterations"
set j 0
foreach c $centroids {
incr j
lassign $c x y
say centroid_${j}_x $x
say centroid_${j}_y $y
}
for {set j 0} {$j < $K} {incr j} {
# lsearch -all returns every matching index; count them.
puts "size_[expr {$j + 1}] [llength [lsearch -all -exact $assignment $j]]"
}
say inertia $inertiaEach point is a two-element list, and lassign unpacks it into x and y.
hypot in expr gives the distance, and nearest replaces its best only
on a strictly smaller distance, so a tie keeps the lower index. Because
every Tcl list is also a string, $new_assignment eq $assignment compares
the whole assignment in one test, which is how we notice nothing changed.
foreach p $points a $assignment walks the points and their clusters side
by side, and lsearch -all returns every index in a cluster, so its
llength is the cluster’s size. Tcl has no clustering or statistics
commands to avoid, and Tcllib, with its matrix and statistics packages,
isn’t installed.