#!/usr/bin/env python3
"""Generate/check the wholly fictional Chapter 11 unsupervised-learning case.

Only Python's standard library is used.  The eight records are designed teaching
points, not observations of people or a natural population.  Running without an
explicit mode is an error; ``--check`` never writes.
"""

import argparse
import csv
import io
import json
import math
from pathlib import Path

CSV_NAME = "ch11-fictional-clustering.csv"
JSON_NAME = "ch11-case-results.json"
FEATURES = ("feature_x", "feature_y", "feature_z")
ROWS = [
    {"point_id": "P01", "feature_x": -5, "feature_y": -1, "feature_z": -1},
    {"point_id": "P02", "feature_x": -5, "feature_y": 1, "feature_z": 1},
    {"point_id": "P03", "feature_x": -3, "feature_y": -1, "feature_z": -1},
    {"point_id": "P04", "feature_x": -3, "feature_y": 1, "feature_z": 1},
    {"point_id": "P05", "feature_x": 3, "feature_y": -1, "feature_z": -1},
    {"point_id": "P06", "feature_x": 3, "feature_y": 1, "feature_z": 1},
    {"point_id": "P07", "feature_x": 5, "feature_y": -1, "feature_z": -1},
    {"point_id": "P08", "feature_x": 5, "feature_y": 1, "feature_z": 1},
]
ROUND_DIGITS = 6


def vector(row):
    return [float(row[name]) for name in FEATURES]


def squared_distance(left, right):
    return sum((a - b) ** 2 for a, b in zip(left, right))


def rounded(value):
    value = round(value, ROUND_DIGITS)
    return 0.0 if value == 0 else value


def cluster_objective(points, assignments, centroids):
    return sum(squared_distance(point, centroids[label]) for point, label in zip(points, assignments))


def kmeans(rows, initial_point_ids, transformed_points=None):
    """Run deterministic Lloyd iterations, resolving equal distances by cluster index."""
    points = transformed_points or [vector(row) for row in rows]
    by_id = {row["point_id"]: point for row, point in zip(rows, points)}
    centroids = [by_id[point_id][:] for point_id in initial_point_ids]
    trace = [{
        "stage": "initial_centroids",
        "centroids": [[rounded(value) for value in center] for center in centroids],
    }]
    previous = None
    for iteration in range(1, 100):
        assignments = [min(range(len(centroids)), key=lambda j: (squared_distance(point, centroids[j]), j)) for point in points]
        if any(j not in assignments for j in range(len(centroids))):
            raise ValueError("empty cluster in fixed teaching initialization")
        new_centroids = []
        for j in range(len(centroids)):
            members = [point for point, label in zip(points, assignments) if label == j]
            new_centroids.append([sum(values) / len(members) for values in zip(*members)])
        objective = cluster_objective(points, assignments, new_centroids)
        trace.append({
            "stage": f"iteration_{iteration}",
            "clusters": [[row["point_id"] for row, label in zip(rows, assignments) if label == j] for j in range(len(centroids))],
            "centroids": [[rounded(value) for value in center] for center in new_centroids],
            "objective_sse": rounded(objective),
        })
        if assignments == previous:
            break
        previous = assignments
        centroids = new_centroids
    else:
        raise AssertionError("fixed k-means run did not converge")
    return {"initial_point_ids": list(initial_point_ids), "trace": trace, "final": trace[-1]}


def partition_objective(points, clusters):
    assignment_by_index = {}
    centroids = []
    for label, indexes in enumerate(clusters):
        members = [points[index] for index in indexes]
        centroids.append([sum(values) / len(members) for values in zip(*members)])
        for index in indexes:
            assignment_by_index[index] = label
    assignments = [assignment_by_index[index] for index in range(len(points))]
    return rounded(cluster_objective(points, assignments, centroids))


def euclidean(left, right):
    return math.sqrt(squared_distance(left, right))


def silhouette(points, clusters):
    """Return mean silhouette; singleton clusters use the conventional value zero."""
    labels = {}
    for label, members in enumerate(clusters):
        for index in members:
            labels[index] = label
    values = []
    for index, point in enumerate(points):
        own = clusters[labels[index]]
        if len(own) == 1:
            values.append(0.0)
            continue
        a = sum(euclidean(point, points[other]) for other in own if other != index) / (len(own) - 1)
        b = min(
            sum(euclidean(point, points[other]) for other in members) / len(members)
            for label, members in enumerate(clusters) if label != labels[index]
        )
        values.append((b - a) / max(a, b) if max(a, b) else 0.0)
    return [rounded(value) for value in values], rounded(sum(values) / len(values))


def build_artifacts():
    points = [vector(row) for row in ROWS]
    vertical = [[0, 1, 2, 3], [4, 5, 6, 7]]
    horizontal = [[0, 2, 4, 6], [1, 3, 5, 7]]
    three_clusters = [[0, 2], [1, 3], [4, 5, 6, 7]]
    four_clusters = [[0, 2], [1, 3], [4, 6], [5, 7]]

    good_run = kmeans(ROWS, ("P01", "P05"))
    poor_run = kmeans(ROWS, ("P01", "P02"))
    standard_deviations = [math.sqrt(17), 1.0, 1.0]
    standardized = [[value / scale for value, scale in zip(point, standard_deviations)] for point in points]
    vertical_silhouettes, vertical_mean = silhouette(points, vertical)
    horizontal_silhouettes, horizontal_mean = silhouette(points, horizontal)

    pca_rows = []
    root_two = math.sqrt(2)
    for row in ROWS:
        pca_rows.append({
            "point_id": row["point_id"],
            "pc1": rounded(row["feature_x"]),
            "pc2": rounded((row["feature_y"] + row["feature_z"]) / root_two),
            "pc3": rounded((row["feature_y"] - row["feature_z"]) / root_two),
        })

    results = {
        "schema_version": 1,
        "title": "Fictional Chapter 11 eight-point unsupervised-learning case",
        "fictional": True,
        "provenance": "Eight designed points fixed in tools/ch11-generate-case.py; not measurements of people or a natural population.",
        "n_points": len(ROWS),
        "distance": "Euclidean distance on the named feature space; k-means objective is within-cluster sum of squared Euclidean distances (SSE).",
        "kmeans_initialization_comparison": {
            "same_data_and_k": 2,
            "left_right_initialization": good_run,
            "upper_lower_initialization": poor_run,
            "lesson": "Lloyd's algorithm reaches different local solutions from different initial centroids; initialization is part of the method specification.",
        },
        "scale_sensitivity": {
            "raw_population_variances": {"feature_x": 17.0, "feature_y": 1.0, "feature_z": 1.0},
            "standardization": "Subtract zero means and divide each feature by its population standard deviation; no labels are used.",
            "raw_sse": {"left_right_partition": partition_objective(points, vertical), "upper_lower_partition": partition_objective(points, horizontal)},
            "standardized_sse": {"left_right_partition": partition_objective(standardized, vertical), "upper_lower_partition": partition_objective(standardized, horizontal)},
            "lesson": "Changing units/scale reverses which of these two partitions has lower SSE.",
        },
        "k_sensitivity": {
            "strategy": "Compare predeclared k values on the same raw features; SSE cannot increase as k grows, so inspect both fit and the substantive usefulness/stability of groups.",
            "candidate_partition_sse": [
                {"k": 1, "sse": partition_objective(points, [list(range(8))])},
                {"k": 2, "sse": partition_objective(points, vertical)},
                {"k": 3, "sse": partition_objective(points, three_clusters)},
                {"k": 4, "sse": partition_objective(points, four_clusters)},
            ],
        },
        "pca": {
            "centering": "All three feature means are zero; covariance divides by n=8 (population convention for this designed finite teaching set).",
            "covariance_matrix": [[17.0, 0.0, 0.0], [0.0, 1.0, 1.0], [0.0, 1.0, 1.0]],
            "eigenvalues_descending": [17.0, 2.0, 0.0],
            "explained_variance_ratio": [rounded(17 / 19), rounded(2 / 19), 0.0],
            "component_vectors": [[1.0, 0.0, 0.0], [0.0, rounded(1 / root_two), rounded(1 / root_two)], [0.0, rounded(1 / root_two), rounded(-1 / root_two)]],
            "projections": pca_rows,
            "sign_note": "An eigenvector and its negation are equivalent; the fixed signs above make the artifact reproducible.",
        },
        "silhouette": {
            "strategy": "For each point, a is mean Euclidean distance to other points in its cluster and b is the smallest mean distance to another cluster; s=(b-a)/max(a,b). A singleton is assigned 0. Report the point values and unweighted mean; do not treat the score alone as proof of natural categories.",
            "left_right_partition": {"point_values": dict(zip((row["point_id"] for row in ROWS), vertical_silhouettes)), "mean": vertical_mean},
            "upper_lower_partition": {"point_values": dict(zip((row["point_id"] for row in ROWS), horizontal_silhouettes)), "mean": horizontal_mean},
        },
        "rounding": f"Computed numeric outputs are rounded to {ROUND_DIGITS} decimal places; CSV contains the exact integer teaching coordinates.",
    }

    csv_buffer = io.StringIO(newline="")
    writer = csv.DictWriter(csv_buffer, fieldnames=ROWS[0].keys(), lineterminator="\r\n")
    writer.writeheader()
    writer.writerows(ROWS)
    artifacts = {
        CSV_NAME: csv_buffer.getvalue().encode("utf-8"),
        JSON_NAME: (json.dumps(results, ensure_ascii=False, indent=2) + "\n").encode("utf-8"),
    }
    return ROWS, results, artifacts


def validate_invariants(rows, results):
    assert len(rows) == 8 and len({row["point_id"] for row in rows}) == 8
    comparison = results["kmeans_initialization_comparison"]
    assert comparison["left_right_initialization"]["final"]["objective_sse"] == 24.0
    assert comparison["upper_lower_initialization"]["final"]["objective_sse"] == 136.0
    scale = results["scale_sensitivity"]
    assert scale["raw_sse"] == {"left_right_partition": 24.0, "upper_lower_partition": 136.0}
    assert scale["standardized_sse"] == {"left_right_partition": 16.470588, "upper_lower_partition": 8.0}
    assert [item["sse"] for item in results["k_sensitivity"]["candidate_partition_sse"]] == [152.0, 24.0, 16.0, 8.0]
    assert results["pca"]["eigenvalues_descending"] == [17.0, 2.0, 0.0]
    assert all(item["pc3"] == 0.0 for item in results["pca"]["projections"])
    assert results["silhouette"]["left_right_partition"]["mean"] > results["silhouette"]["upper_lower_partition"]["mean"]


def parse_args(argv=None):
    parser = argparse.ArgumentParser(description="Generate or read-only check the fictional Chapter 11 eight-point case.")
    modes = parser.add_mutually_exclusive_group(required=True)
    modes.add_argument("--check", action="store_true", help="compare generated bytes with files without writing")
    modes.add_argument("--write", action="store_true", help="explicitly overwrite generated files")
    parser.add_argument("--output-dir", type=Path, help="override output directory (mainly for isolated tests)")
    return parser.parse_args(argv)


def main(argv=None):
    args = parse_args(argv)
    rows, results, artifacts = build_artifacts()
    validate_invariants(rows, results)
    script_path = Path(__file__).resolve()
    root = script_path.parents[1]
    default_dir = script_path.parent if script_path.parent.name == "data" and script_path.parent.parent.name == "public" else root / "public" / "data"
    output_dir = args.output_dir or default_dir
    if args.check:
        mismatches = []
        for name, expected in artifacts.items():
            path = output_dir / name
            if not path.is_file() or path.read_bytes() != expected:
                mismatches.append(str(path))
        if mismatches:
            raise SystemExit("CHECK FAILED: generated bytes differ: " + ", ".join(mismatches))
        print("CHECK PASSED: 8 points; k-means SSE 24/136; PCA eigenvalues 17/2/0; CSV/JSON bytes match")
        return 0
    output_dir.mkdir(parents=True, exist_ok=True)
    for name, data in artifacts.items():
        (output_dir / name).write_bytes(data)
    print(json.dumps(results, ensure_ascii=False, indent=2))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())