B-plot writes when the script runsDictionary of Applied Machine Learning · $k$-means
Numerical companion to the entry $k$-means: it recomputes what the entry states and prints one line per check
Numerical companion to the glossary entry 'kmeans' ($k$-means). A photograph of the Oetscher massif (assets/oetscher.jpg, the same photograph and the same patches the 'gmm' entry uses) is cut into square patches. Each patch is a data point whose feature vector holds two numbers, how green and how blue the patch is on average, so the clusters and their centroids can be drawn in the plane.
Run it without installing anything:uv run https://dictionaryofml.org/terms/kmeans.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download kmeans.py into a folder and run uv run kmeans.py there. With NumPy and Matplotlib already installed, python3 kmeans.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download kmeans.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""k-means with three clusters on the greenness and blueness of image
patches: vegetation, mountain and haze, and sky, each patch assigned to
exactly one of them.
Purpose
-------
Numerical companion to the glossary entry 'kmeans' ($k$-means). A
photograph of the Oetscher massif (assets/oetscher.jpg, the same
photograph and the same patches the 'gmm' entry uses) is cut into
square patches. Each patch is a data point whose feature vector holds
two numbers, how green and how blue the patch is on average, so the
clusters and their centroids can be drawn in the plane.
The demo checks the entry's claims: Lloyd's algorithm never increases
the clustering error and reaches a fixed point after finitely many
iterations, where each centroid is the mean of the data points assigned
to it; the three clusters are the vegetation, the mountain with the
haze, and the sky; the assignment is a hard clustering, so the
per-cluster copies of the photograph add back up to the photograph; the
assignment coincides with the EM algorithm for a GMM whose cluster
probabilities are equal and whose covariance matrices are a small
multiple of the identity matrix; the clustering error decreases
with the number of clusters; and, for the k-medians variant that
averages distances instead of squared distances, the centroid step
computes a geometric median (Weiszfeld's iteration), which a few
far-away patches barely move while they shift the mean.
Deterministic: the centroids are initialized by splitting the patches
into k equal parts ordered by blueness (no randomness). Self-contained:
numpy + matplotlib only.
Blocks
------
[B-patches] Cut the photograph into 128x128 patches and measure how
green and how blue each one is; check the patch count.
[B-lloyd] Run Lloyd's algorithm with k = 3: assign every patch to its
nearest centroid, then move each centroid to the mean of
its patches. Check that the clustering error never
increases, that the assignments stop changing after
finitely many iterations, and that at the fixed point each
centroid is the mean of its patches.
[B-clusters] The three clusters are the vegetation, the mountain with
the haze, and the sky. Check each against the average
color of its patches and against where they sit in the
photograph.
[B-hard] Hard clustering: every patch belongs to exactly one
cluster. Check that the membership indicators of a patch
sum to one and that no patch is shared.
[B-images] The photograph with the patch grid drawn on it, the
photograph at patch resolution, and one copy per cluster
that keeps the patches assigned to that cluster and blacks
out the rest. Because every patch is assigned to exactly
one cluster, the three copies add back up to the
photograph.
[B-gmm] The EM algorithm for a GMM with equal cluster probabilities
and covariance matrices that are a small multiple of the
identity matrix assigns the patches like k-means. Check
the agreement.
[B-k] The effect of the number of clusters: k = 2, 3, 4, each
patch painted with the average color of its cluster. Check
that the clustering error decreases with k.
[B-medians] k-medians: the loss is the Euclidean distance itself
instead of its square. The assignment step is unchanged,
the centroid step becomes the geometric median of the
assigned patches, computed by Weiszfeld's fixed-point
iteration. Check that the average distance never
increases, that at the fixed point the unit vectors from
each centroid to its patches sum to zero, that the
assignment mostly agrees with k-means, and that ten
far-away patches added to the vegetation cluster move its
mean much more than its geometric median.
[B-plot] Write the patch scatter, the three centroids and the
boundaries between the clusters for the entry's figure,
plus the preview. Each boundary is the perpendicular
bisector of the line connecting two centroids: check that
it crosses that line at its midpoint at a right angle.
Outputs
-------
kmeans_patches.csv : green, blue -- every patch
kmeans_points.csv : x1, x2 -- every patch (the figure's scatter)
kmeans_centroids.csv : x1, x2 -- the three centroids
kmeans_boundary1.csv, kmeans_boundary2.csv, kmeans_boundary3.csv :
x1, x2 -- the two end points of each boundary
segment between two clusters (empty when the two
regions do not touch inside the figure's box)
kmeans_rightangle1.csv, kmeans_rightangle3.csv :
x1, x2 -- three corners of a small square marking
the right angle where the boundary crosses the
line connecting the two centroids, at its midpoint
kmeans_oetscher_raster.png : the photograph with the patch grid
kmeans_oetscher_original.png : the photograph at patch resolution
kmeans_oetscher_vegetation.png : the patches assigned to the vegetation
kmeans_oetscher_mountain.png : the patches assigned to the mountain
kmeans_oetscher_sky.png : the patches assigned to the sky
kmeans_oetscher_k2.png, _k3.png, _k4.png : each patch painted with the
average color of its cluster for k = 2, 3, 4
kmeans.png : preview (checking only) -- the patches, the
centroids and the cluster boundaries
"""
from pathlib import Path
import urllib.request
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.image import imread
OUT_DIR = Path(__file__).parent
ASSET_URL = "https://dictionaryofml.org/demoassets/"
def demo_asset(name):
"""Path to the input file `name`, fetched once if it is not already here.
Inside the repository the committed copy is used (pythondemos/ first,
then assets/), so a run there reproduces the committed outputs exactly.
Run from anywhere else -- from the copy published on dictionaryofml.org,
say -- the file is downloaded and cached beside the script, so the demo
is runnable without the repository.
"""
for cand in (OUT_DIR / name, OUT_DIR.parent / "assets" / name):
if cand.exists():
return cand
cached = OUT_DIR / name
print(f" fetching {name} from {ASSET_URL}{name}")
with urllib.request.urlopen(ASSET_URL + name, timeout=180) as resp:
cached.write_bytes(resp.read())
return cached
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
/Users/junga1/dictionaryappliedml/pythondemos/kmeans.py:500: DeprecationWarning: Arrays of 2-dimensional vectors are deprecated. Use arrays of 3-dimensional vectors instead. (deprecated in NumPy 2.0) abs(np.cross(direction, mid - seg[0])) < 1e-9 photograph cut into 27 x 36 = 972 patches of 128 x 128 pixels
Cut the photograph into 128x128 patches and measure how green and how blue each one is; check the patch count.
PATCH = 128
photo = imread(demo_asset("oetscher.jpg")) / 255.0
rows = photo.shape[0] // PATCH
cols = photo.shape[1] // PATCH
tiles = photo[:rows * PATCH, :cols * PATCH].reshape(
rows, PATCH, cols, PATCH, 3).mean(axis=(1, 3))
rgb = tiles.reshape(-1, 3)
green = rgb[:, 1] # average greenness, on a 0 to 1 scale
blue = rgb[:, 2] # average blueness, on the same scale
X = np.column_stack([green, blue])
with open(OUT_DIR / "kmeans_patches.csv", "w") as f:
f.write("green,blue\n")
for a, b in X:
f.write(f"{a:.4f},{b:.4f}\n")
print(f" photograph cut into {rows} x {cols} = {len(X)} patches "
f"of {PATCH} x {PATCH} pixels")
check("[B-patches] every patch has a two-number feature vector",
X.shape == (rows * cols, 2))
[ok] [B-patches] every patch has a two-number feature vector Lloyd's algorithm: 12 iterations, clustering error 0.0269 -> 0.0179
Run Lloyd's algorithm with k = 3: assign every patch to its nearest centroid, then move each centroid to the mean of its patches. Check that the clustering error never increases, that the assignments stop changing after finitely many iterations, and that at the fixed point each centroid is the mean of its patches.
NRCLUSTER = 3
def clustering_error(A, cents):
"""Average squared Euclidean distance to the nearest centroid."""
d2 = ((A[:, None, :] - cents) ** 2).sum(axis=2)
return float(d2.min(axis=1).mean())
def assign(A, cents):
"""Index of the nearest centroid for each row of A."""
return ((A[:, None, :] - cents) ** 2).sum(axis=2).argmin(axis=1)
def lloyd(A, cents):
"""Lloyd's algorithm from the given centroids; returns the trace too."""
errors = [clustering_error(A, cents)]
lab = assign(A, cents)
for _ in range(300):
new = np.stack([A[lab == c].mean(axis=0) for c in range(len(cents))])
new_lab = assign(A, new)
errors.append(clustering_error(A, new))
if np.array_equal(new_lab, lab) and np.allclose(new, cents):
return new_lab, new, errors
lab, cents = new_lab, new
return lab, cents, errors
def start_from_blueness(A, k):
"""k initial centroids: means of k equal parts ordered by blueness."""
parts = np.array_split(np.argsort(A[:, 1]), k)
return np.stack([A[g].mean(axis=0) for g in parts])
start = start_from_blueness(X, NRCLUSTER)
label, centroids, errors = lloyd(X, start.copy())
order = np.argsort(centroids[:, 1]) # greenest first, sky last
centroids = centroids[order]
label = np.argsort(order)[label]
print(f" Lloyd's algorithm: {len(errors) - 1} iterations, clustering error "
f"{errors[0]:.4f} -> {errors[-1]:.4f}")
check("[B-lloyd] no iteration increases the clustering error",
bool(np.all(np.diff(errors) <= 1e-12)))
check("[B-lloyd] the assignments stop changing after finitely many "
"iterations", len(errors) - 1 < 300)
check("[B-lloyd] at the fixed point each centroid is the mean of its patches",
all(np.allclose(centroids[c], X[label == c].mean(axis=0))
for c in range(NRCLUSTER)))
check("[B-lloyd] one more iteration changes nothing",
np.array_equal(assign(X, centroids), label))
[ok] [B-lloyd] no iteration increases the clustering error [ok] [B-lloyd] the assignments stop changing after finitely many iterations [ok] [B-lloyd] at the fixed point each centroid is the mean of its patches [ok] [B-lloyd] one more iteration changes nothing vegetation 566 patches, centroid greenness 0.287, blueness 0.143, average color RGB [0.23 0.29 0.14], 21% of its patches in the upper half mountain and haze 202 patches, centroid greenness 0.441, blueness 0.563, average color RGB [0.27 0.44 0.56], 89% of its patches in the upper half bright sky 204 patches, centroid greenness 0.672, blueness 0.943, average color RGB [0.36 0.67 0.94], 100% of its patches in the upper half
The three clusters are the vegetation, the mountain with the haze, and the sky. Check each against the average color of its patches and against where they sit in the photograph.
NAMES = ("vegetation", "mountain and haze", "bright sky")
where = np.argwhere(np.ones((rows, cols), dtype=bool)) # (row, col) per patch
for c in range(NRCLUSTER):
sel = label == c
print(f" {NAMES[c]:<18} {int(sel.sum()):3d} patches, centroid greenness "
f"{centroids[c, 0]:.3f}, blueness {centroids[c, 1]:.3f}, average "
f"color RGB {np.round(rgb[sel].mean(axis=0), 2)}, "
f"{(where[sel][:, 0] < rows / 2).mean():.0%} of its patches in the "
f"upper half")
check("[B-clusters] every average lies on the 0 to 1 scale",
float(X.min()) >= 0.0 and float(X.max()) <= 1.0)
check("[B-clusters] the sky is the bluest cluster",
centroids[2, 1] == centroids[:, 1].max())
check("[B-clusters] the sky sits in the upper half of the photograph",
float((where[label == 2][:, 0] < rows / 2).mean()) > 0.8)
check("[B-clusters] the vegetation is the least blue and sits mostly in the "
"lower half", centroids[0, 1] == centroids[:, 1].min()
and float((where[label == 0][:, 0] >= rows / 2).mean()) > 0.7)
[ok] [B-clusters] every average lies on the 0 to 1 scale [ok] [B-clusters] the sky is the bluest cluster [ok] [B-clusters] the sky sits in the upper half of the photograph [ok] [B-clusters] the vegetation is the least blue and sits mostly in the lower half 972 of 972 patches are assigned to exactly one cluster
Hard clustering: every patch belongs to exactly one cluster. Check that the membership indicators of a patch sum to one and that no patch is shared.
member = np.stack([(label == c).astype(float) for c in range(NRCLUSTER)])
print(f" {len(X)} of {len(X)} patches are assigned to exactly one cluster")
check("[B-hard] the membership indicators of each patch sum to one",
np.allclose(member.sum(axis=0), 1.0))
check("[B-hard] no patch is shared between clusters",
bool(np.all(member.max(axis=0) == 1.0)))
[ok] [B-hard] the membership indicators of each patch sum to one [ok] [B-hard] no patch is shared between clusters
The photograph with the patch grid drawn on it, the photograph at patch resolution, and one copy per cluster that keeps the patches assigned to that cluster and blacks out the rest. Because every patch is assigned to exactly one cluster, the three copies add back up to the photograph.
def shrink_to_width(img, width_px):
"""Area-average an image down to width_px, keeping its aspect ratio.
Each output pixel is the mean of the input it covers, which is the right
filter for a photograph and needs nothing beyond numpy.
"""
if img.shape[1] <= width_px:
return np.asarray(img, dtype=float)
height_px = round(img.shape[0] * width_px / img.shape[1])
out = np.asarray(img, dtype=float)
for axis, n_out in ((1, width_px), (0, height_px)):
n_in = out.shape[axis]
edges = np.linspace(0.0, n_in, n_out + 1)
cum = np.cumsum(out, axis=axis)
cum = np.concatenate(
[np.zeros_like(np.take(cum, [0], axis=axis)), cum], axis=axis)
grid = np.arange(n_in + 1)
cum = np.apply_along_axis(
lambda col, g=grid, e=edges: np.interp(e, g, col), axis, cum)
shape = [1] * out.ndim
shape[axis] = n_out
out = np.diff(cum, axis=axis) / np.diff(edges).reshape(shape)
return out
def save_image(arr, name, zoom=6):
"""Write an RGB array as a PNG, enlarged so the patches stay visible."""
img = np.clip(arr, 0.0, 1.0).repeat(zoom, axis=0).repeat(zoom, axis=1)
plt.imsave(OUT_DIR / name, img)
# The raster is printed at 0.5\textwidth, and the book's textwidth is 5.40 in,
# so it is 2.70 in across on the page; at 300 ppi that is 810 px against the
# 1152 px of the full-detail view, which went into the PDF at over 400 ppi.
# The photograph is shrunk first and the patch rules drawn afterwards, so every
# rule stays exactly one white pixel; shrinking the ruled image would blur the
# rules, which are what shows the reader the extent of one data point.
SHRINK, TARGET_DPI = 4, 300
RASTER_WIDTH_IN = 2.70
view = photo[:rows * PATCH, :cols * PATCH:, :][::SHRINK, ::SHRINK]
cell = -(-round(RASTER_WIDTH_IN * TARGET_DPI) // cols) # ceil: at least 300 ppi
view = shrink_to_width(view, cols * cell)
view[::cell, :, :] = 1.0 # horizontal rules
view[:, ::cell, :] = 1.0 # vertical rules
view[-1, :, :] = 1.0
view[:, -1, :] = 1.0
save_image(view, "kmeans_oetscher_raster.png", zoom=1)
check("[B-images] the raster has one cell per patch",
view.shape[:2] == (rows * cell, cols * cell))
save_image(tiles, "kmeans_oetscher_original.png")
COPIES = ("kmeans_oetscher_vegetation.png", "kmeans_oetscher_mountain.png",
"kmeans_oetscher_sky.png")
kept = [tiles * member[c].reshape(rows, cols, 1) for c in range(NRCLUSTER)]
for fname, arr in zip(COPIES, kept):
save_image(arr, fname)
check("[B-images] one copy of the photograph per cluster",
all((OUT_DIR / f).exists() for f in COPIES))
check("[B-images] the three copies add back up to the photograph",
np.allclose(sum(kept), tiles))
[ok] [B-images] the raster has one cell per patch [ok] [B-images] one copy of the photograph per cluster [ok] [B-images] the three copies add back up to the photograph constrained GMM, variance 1e-4: assignment agrees with k-means on 99.4% of patches, 100% of posteriors above 0.99
The EM algorithm for a GMM with equal cluster probabilities and covariance matrices that are a small multiple of the identity matrix assigns the patches like k-means. Check the agreement.
def spherical_em(A, cents, var):
"""EM with equal cluster probabilities and covariance matrices var*I."""
post = None
for _ in range(300):
d2 = ((A[:, None, :] - cents) ** 2).sum(axis=2)
logp = -d2 / (2.0 * var)
post = np.exp(logp - logp.max(axis=1, keepdims=True))
post /= post.sum(axis=1, keepdims=True)
new = (post.T @ A) / post.sum(axis=0)[:, None]
if np.allclose(new, cents):
break
cents = new
return post.argmax(axis=1), cents, post
lab_em, cent_em, post = spherical_em(X, start.copy(), 1e-4)
lab_em = np.argsort(order)[lab_em]
agree = float((lab_em == label).mean())
print(f" constrained GMM, variance 1e-4: assignment agrees with k-means on "
f"{agree:.1%} of patches, {(post.max(axis=1) > 0.99).mean():.0%} of "
f"posteriors above 0.99")
check("[B-gmm] the constrained GMM assigns the patches like k-means",
agree > 0.99)
check("[B-gmm] its centroids agree with the k-means centroids",
np.abs(cent_em[order] - centroids).max() < 0.01)
[ok] [B-gmm] the constrained GMM assigns the patches like k-means [ok] [B-gmm] its centroids agree with the k-means centroids k = 2: clustering error 0.0338, cluster sizes [634, 338] k = 3: clustering error 0.0179, cluster sizes [566, 202, 204] k = 4: clustering error 0.0110, cluster sizes [360, 208, 203, 201]
The effect of the number of clusters: k = 2, 3, 4, each patch painted with the average color of its cluster. Check that the clustering error decreases with k.
errors_by_k = {}
for k in (2, 3, 4):
lab_k, cent_k, err_k = lloyd(X, start_from_blueness(X, k))
errors_by_k[k] = err_k[-1]
palette = np.stack([rgb[lab_k == c].mean(axis=0) for c in range(k)])
save_image(palette[lab_k].reshape(rows, cols, 3), f"kmeans_oetscher_k{k}.png")
print(f" k = {k}: clustering error {err_k[-1]:.4f}, "
f"cluster sizes {np.bincount(lab_k).tolist()}")
check("[B-k] the clustering error decreases with the number of clusters",
errors_by_k[2] > errors_by_k[3] > errors_by_k[4])
check("[B-k] k = 3 reproduces the clustering of the entry",
abs(errors_by_k[3] - errors[-1]) < 1e-12)
[ok] [B-k] the clustering error decreases with the number of clusters [ok] [B-k] k = 3 reproduces the clustering of the entry k-medians from the k-means centroids: 5 iterations, average distance 0.1163 -> 0.1150, assignment agrees with k-means on 98.5% of patches
k-medians: the loss is the Euclidean distance itself instead of its square. The assignment step is unchanged, the centroid step becomes the geometric median of the assigned patches, computed by Weiszfeld's fixed-point iteration. Check that the average distance never increases, that at the fixed point the unit vectors from each centroid to its patches sum to zero, that the assignment mostly agrees with k-means, and that ten far-away patches added to the vegetation cluster move its mean much more than its geometric median.
def geometric_median(A, z):
"""Weiszfeld's fixed-point iteration for the geometric median of A's rows."""
for _ in range(1000):
d = np.maximum(np.linalg.norm(A - z, axis=1), 1e-12)
new = (A / d[:, None]).sum(axis=0) / (1.0 / d).sum()
if np.linalg.norm(new - z) < 1e-12:
return new
z = new
return z
def median_error(A, cents):
"""Average Euclidean distance to the nearest centroid."""
return float(np.linalg.norm(A[:, None, :] - cents, axis=2).min(axis=1).mean())
medians = centroids.copy()
lab_med = assign(X, medians)
med_errors = [median_error(X, medians)]
for _ in range(300):
new = np.stack([geometric_median(X[lab_med == c], medians[c])
for c in range(NRCLUSTER)])
new_lab = assign(X, new)
med_errors.append(median_error(X, new))
if np.array_equal(new_lab, lab_med) and np.allclose(new, medians, atol=1e-10):
break
lab_med, medians = new_lab, new
agree_med = float((lab_med == label).mean())
print(f" k-medians from the k-means centroids: {len(med_errors) - 1} iterations, "
f"average distance {med_errors[0]:.4f} -> {med_errors[-1]:.4f}, assignment "
f"agrees with k-means on {agree_med:.1%} of patches")
check("[B-medians] the assignment step is unchanged: the nearest centroid is "
"the same under the distance and the squared distance",
np.array_equal(assign(X, medians),
np.linalg.norm(X[:, None, :] - medians, axis=2).argmin(axis=1)))
check("[B-medians] no iteration increases the average distance",
bool(np.all(np.diff(med_errors) <= 1e-12)))
check("[B-medians] at the fixed point the unit vectors from each centroid to "
"its patches sum to zero",
all(np.linalg.norm(((X[lab_med == c] - medians[c])
/ np.linalg.norm(X[lab_med == c] - medians[c],
axis=1)[:, None]).sum(axis=0)) < 1e-6
for c in range(NRCLUSTER)))
check("[B-medians] k-medians assigns most patches like k-means", agree_med > 0.95)
vegetation = X[label == 0]
spiked = np.vstack([vegetation, np.tile([2.0, 2.0], (10, 1))])
mean_shift = float(np.linalg.norm(spiked.mean(axis=0) - vegetation.mean(axis=0)))
median_shift = float(np.linalg.norm(
geometric_median(spiked, vegetation.mean(axis=0))
- geometric_median(vegetation, vegetation.mean(axis=0))))
print(f" ten far-away patches added to the {len(vegetation)} vegetation patches: "
f"the mean moves by {mean_shift:.3f}, the geometric median by "
f"{median_shift:.3f}")
check("[B-medians] the far-away patches move the mean at least five times as "
"far as the geometric median", mean_shift > 5.0 * median_shift)
[ok] [B-medians] the assignment step is unchanged: the nearest centroid is the same under the distance and the squared distance [ok] [B-medians] no iteration increases the average distance [ok] [B-medians] at the fixed point the unit vectors from each centroid to its patches sum to zero [ok] [B-medians] k-medians assigns most patches like k-means ten far-away patches added to the 566 vegetation patches: the mean moves by 0.044, the geometric median by 0.004 [ok] [B-medians] the far-away patches move the mean at least five times as far as the geometric median
Write the patch scatter, the three centroids and the boundaries between the clusters for the entry's figure, plus the preview. Each boundary is the perpendicular bisector of the line connecting two centroids: check that it crosses that line at its midpoint at a right angle.
BOX = (0.0, 1.12, 0.0, 1.05) # xmin, xmax, ymin, ymax of the figure
def boundary(cents, i, j, box):
"""End points of the boundary between clusters i and j inside the box.
The boundary is the part of the bisector of the two centroids on
which no third centroid is closer, clipped to the box.
"""
mid = (cents[i] + cents[j]) / 2.0
normal = cents[j] - cents[i]
direction = np.array([-normal[1], normal[0]])
direction /= np.linalg.norm(direction)
lo, hi = -np.inf, np.inf
# closer to i (and j) than to every other centroid: linear in t
for k in range(len(cents)):
if k in (i, j):
continue
a = 2.0 * direction @ (cents[k] - cents[i])
b = cents[k] @ cents[k] - cents[i] @ cents[i] \
- 2.0 * mid @ (cents[k] - cents[i])
if a > 0:
hi = min(hi, b / a)
elif a < 0:
lo = max(lo, b / a)
elif b < 0:
return None
# inside the box: four more linear constraints
for coord, (lo_c, hi_c) in enumerate(((box[0], box[1]), (box[2], box[3]))):
d = direction[coord]
for bound, sign in ((lo_c, -1.0), (hi_c, 1.0)):
a = sign * d
b = sign * (bound - mid[coord])
if a > 0:
hi = min(hi, b / a)
elif a < 0:
lo = max(lo, b / a)
elif b < 0:
return None
if lo >= hi:
return None
return np.stack([mid + lo * direction, mid + hi * direction])
with open(OUT_DIR / "kmeans_points.csv", "w") as f:
f.write("x1,x2\n")
for a, b in X:
f.write(f"{a:.4f},{b:.4f}\n")
with open(OUT_DIR / "kmeans_centroids.csv", "w") as f:
f.write("x1,x2\n")
for a, b in centroids:
f.write(f"{a:.4f},{b:.4f}\n")
segments = []
for n, (i, j) in enumerate(((0, 1), (0, 2), (1, 2)), start=1):
seg = boundary(centroids, i, j, BOX)
segments.append(seg)
with open(OUT_DIR / f"kmeans_boundary{n}.csv", "w") as f:
f.write("x1,x2\n")
if seg is not None:
for a, b in seg:
f.write(f"{a:.4f},{b:.4f}\n")
check("[B-plot] every boundary segment lies inside the box",
all(seg is None or (BOX[0] - 1e-9 <= seg[:, 0].min()
and seg[:, 0].max() <= BOX[1] + 1e-9
and BOX[2] - 1e-9 <= seg[:, 1].min()
and seg[:, 1].max() <= BOX[3] + 1e-9)
for seg in segments))
check("[B-plot] the boundaries separate the clusters",
np.array_equal(assign(X, centroids), label))
SQUARE = 0.03 # side of the right-angle marker, data units
right_angles = {}
for n, (i, j) in ((1, (0, 1)), (3, (1, 2))):
mid = (centroids[i] + centroids[j]) / 2.0
along = centroids[j] - centroids[i]
along /= np.linalg.norm(along)
across = np.array([-along[1], along[0]])
corners = np.stack([mid + SQUARE * along, mid + SQUARE * (along + across),
mid + SQUARE * across])
right_angles[n] = corners
with open(OUT_DIR / f"kmeans_rightangle{n}.csv", "w") as f:
f.write("x1,x2\n")
for a, b in corners:
f.write(f"{a:.4f},{b:.4f}\n")
seg = segments[n - 1]
direction = seg[1] - seg[0]
check(f"[B-plot] boundary {n} crosses the connecting line at its midpoint "
"at a right angle",
abs(np.cross(direction, mid - seg[0])) < 1e-9
and abs(direction @ along) < 1e-9)
fig, ax = plt.subplots(figsize=(6.2, 5.0))
ax.plot(X[:, 0], X[:, 1], "o", color="0.6", markersize=2.2,
linestyle="none", label="patch")
ax.plot(centroids[:, 0], centroids[:, 1], "x", color="black", markersize=9,
markeredgewidth=2, linestyle="none", label="cluster centroid")
for n, seg in enumerate(segments):
if seg is not None:
ax.plot(seg[:, 0], seg[:, 1], color="black", linewidth=1.4,
linestyle="--", label="cluster boundary" if n == 0 else None)
ax.plot(centroids[:, 0], centroids[:, 1], color="black", linewidth=0.8,
label="line connecting two centroids")
for corners in right_angles.values():
ax.plot(corners[:, 0], corners[:, 1], color="black", linewidth=0.8)
ax.set_aspect("equal")
for c in range(NRCLUSTER):
ax.annotate(NAMES[c], centroids[c], xytext=(6, 6),
textcoords="offset points", fontsize="small")
ax.set_xlim(BOX[0], BOX[1])
ax.set_ylim(BOX[2], BOX[3])
ax.set_xlabel("average greenness of the patch")
ax.set_ylabel("average blueness of the patch")
ax.set_title("Patches of the Oetscher photograph and k-means with k = 3")
ax.legend(frameon=False, loc="lower right", fontsize="small")
fig.tight_layout()
fig.savefig(OUT_DIR / "kmeans.png", dpi=150)
plt.close(fig)
check("[B-plot] the scatter, the centroids and the boundaries were written",
(OUT_DIR / "kmeans_points.csv").exists()
and (OUT_DIR / "kmeans_boundary3.csv").exists())
passed = sum(1 for _, ok in report if ok)
print(f"\n{passed}/{len(report)} checks pass")
[ok] [B-plot] every boundary segment lies inside the box [ok] [B-plot] the boundaries separate the clusters [ok] [B-plot] boundary 1 crosses the connecting line at its midpoint at a right angle [ok] [B-plot] boundary 3 crosses the connecting line at its midpoint at a right angle [ok] [B-plot] the scatter, the centroids and the boundaries were written 28/28 checks pass
B-plot writes when the script runs