* MSC-P-032 segmentation stability. MIT License.
BEGIN PROGRAM Python3.
import csv, math
required = ["customer_id","recency_days","frequency_12m","avg_basket_eur"]
path = "public/datasets/msc-p032-stability-panel.csv"
with open(path, newline="", encoding="utf-8") as handle:
    reader = csv.DictReader(handle)
    if reader.fieldnames != required: raise ValueError("exact schema required")
    rows = list(reader)
if len(rows) < 100: raise ValueError("too few customers for the declared stability analysis")
if any(set(r) != set(required) or any(r[k] is None for k in required) for r in rows): raise ValueError("exactly four cells required per row")
if [r["customer_id"].strip() for r in rows] != ["T%04d" % i for i in range(1, len(rows) + 1)]: raise ValueError("customers must be ordered T0001, T0002, ... without gaps")
X = []
for r in rows:
    v = [float(r["recency_days"]), float(r["frequency_12m"]), float(r["avg_basket_eur"])]
    if any(not math.isfinite(x) or x <= 0 for x in v): raise ValueError("recency, frequency and basket must be finite and strictly positive")
    X.append([math.log(v[0]), v[1], math.log(v[2])])
CANDIDATES, PROPOSED, B, ITERS = (2, 3, 4), 4, 40, 100
RECOVERY_MIN, MARGIN_MIN = 0.75, 0.10
SEED, MULT, INC, MOD = 20260912, 1103515245, 12345, 2147483648
n = len(X)
means = [sum(c) / n for c in zip(*X)]
devs = [math.sqrt(sum((x - m) ** 2 for x in c) / (n - 1)) for c, m in zip(zip(*X), means)]
if any(dv <= 0 for dv in devs): raise ValueError("a declared feature has no variation")
Z = [[(x - m) / dv for x, m, dv in zip(row, means, devs)] for row in X]
dist = lambda a, b: sum((x - y) ** 2 for x, y in zip(a, b))
class Stream:
    def __init__(self, seed): self.state = seed
    def draw(self, size):
        self.state = (MULT * self.state + INC) % MOD
        return self.state % size
def start(points, k):
    order = sorted(range(len(points)), key=lambda i: (sum(points[i]), i))
    return [list(points[order[int(round((j + 0.5) / k * (len(order) - 1)))]]) for j in range(k)]
def assign(points, centroids):
    return [min(range(len(centroids)), key=lambda k: (dist(p, centroids[k]), k)) for p in points]
def cluster(points, k):
    centroids = start(points, k); labels = [-1] * len(points)
    for _ in range(ITERS):
        new = assign(points, centroids)
        changed = new != labels; labels = new
        for j in range(k):
            members = [p for p, l in zip(points, labels) if l == j]
            if not members: raise ValueError("the declared segmentation collapsed to fewer segments")
            centroids[j] = [sum(c) / len(members) for c in zip(*members)]
        if not changed: break
    return labels, centroids
def jaccard(a, b):
    u = len(a | b)
    return len(a & b) / u if u else 0.0
def recovery(points, k):
    labels, _ = cluster(points, k)
    original = [{i for i, l in enumerate(labels) if l == j} for j in range(k)]
    totals = [0.0] * k
    stream = Stream(SEED)
    for _ in range(B):
        sample = [points[stream.draw(len(points))] for _ in range(len(points))]
        _, centroids = cluster(sample, k)
        rebuilt_labels = assign(points, centroids)
        rebuilt = [{i for i, l in enumerate(rebuilt_labels) if l == j} for j in range(k)]
        for j in range(k):
            totals[j] += max(jaccard(original[j], candidate) for candidate in rebuilt)
    return [t / B for t in totals]
def scrambled(points):
    stream = Stream(SEED); columns = []
    for column in zip(*points):
        values = list(column)
        for i in range(len(values) - 1, 0, -1):
            s = stream.draw(i + 1)
            values[i], values[s] = values[s], values[i]
        columns.append(values)
    return [list(row) for row in zip(*columns)]
reference = scrambled(Z)
print("customers=%d" % n)
print("resamples=%d" % B)
print("proposed_segments=%d" % PROPOSED)
flags = {}
for k in CANDIDATES:
    observed = recovery(Z, k); null = recovery(reference, k)
    weakest = min(observed); null_weakest = min(null); margin = weakest - null_weakest
    flags[k] = "STABLE" if weakest >= RECOVERY_MIN and margin >= MARGIN_MIN else "NOT_STABLE"
    print("recovery_k%d=%s" % (k, ",".join("%.6f" % v for v in observed)))
    print("weakest_k%d=%.6f" % (k, weakest))
    print("null_weakest_k%d=%.6f" % (k, null_weakest))
    print("margin_k%d=%.6f" % (k, margin))
    print("flag_k%d=%s" % (k, flags[k]))
print("verdict=%s" % ("PROPOSED_SEGMENTATION_STABLE" if flags[PROPOSED] == "STABLE" else "PROPOSED_SEGMENTATION_NOT_STABLE"))
END PROGRAM.