# MSC-P-031 customer segmentation. MIT License.
# Same data, same declared construction and same printed contract as the Python reference.
args <- commandArgs(trailingOnly = TRUE)
path <- if (length(args)) args[[1]] else file.path("public", "datasets", "msc-p031-segmentation-panel.csv")
raw_lines <- readLines(path, warn = FALSE, encoding = "UTF-8")
if (raw_lines[[1]] != "customer_id,recency_days,frequency_12m,avg_basket_eur,category_breadth,digital_share,next_quarter_revenue_eur") stop("exact schema required")
if (any(lengths(regmatches(raw_lines[-1], gregexpr(",", raw_lines[-1], fixed = TRUE))) != 6L)) stop("exactly seven cells required per row")
d <- read.csv(path, colClasses = "character", check.names = FALSE)
required <- c("customer_id", "recency_days", "frequency_12m", "avg_basket_eur", "category_breadth", "digital_share", "next_quarter_revenue_eur")
if (!identical(names(d), required)) stop("exact schema required")
n <- nrow(d)
if (n < 100L) stop("too few customers for the declared segmentation")
if (!identical(d$customer_id, sprintf("S%04d", seq_len(n)))) stop("customers must be ordered S0001, S0002, ... without gaps")
for (name in c("recency_days", "frequency_12m", "avg_basket_eur", "category_breadth")) {
  d[[name]] <- suppressWarnings(as.numeric(d[[name]]))
  if (any(!is.finite(d[[name]])) || any(d[[name]] <= 0)) stop("recency, frequency, basket and breadth must be finite and strictly positive")
}
d$digital_share <- suppressWarnings(as.numeric(d$digital_share))
if (any(!is.finite(d$digital_share)) || any(d$digital_share <= 0) || any(d$digital_share >= 1)) stop("the digital share must lie strictly between zero and one")
d$next_quarter_revenue_eur <- suppressWarnings(as.numeric(d$next_quarter_revenue_eur))
if (any(!is.finite(d$next_quarter_revenue_eur)) || any(d$next_quarter_revenue_eur < 0)) stop("the outcome must be finite and non-negative")

SEGMENTS <- 4L; MAX_ITERATIONS <- 100L; STEP <- 3L
SILHOUETTE_MIN <- 0.25; SHARE_MIN <- 0.05; STABILITY_MIN <- 0.60; USEFULNESS_MIN <- 1.50
X <- cbind(log(d$recency_days), d$frequency_12m, log(d$avg_basket_eur), d$category_breadth, d$digital_share)
feature_names <- c("log_recency", "frequency", "log_basket", "breadth", "digital")
centers <- colMeans(X); deviations <- apply(X, 2, sd)
if (any(deviations <= 0)) stop("a declared feature has no variation")
Z <- scale(X, center = centers, scale = deviations)
start_centroids <- function(points) {
  order_index <- order(rowSums(points), seq_len(nrow(points)))
  positions <- round(c(0.125, 0.375, 0.625, 0.875) * (nrow(points) - 1)) + 1
  points[order_index[positions], , drop = FALSE]
}
assign_labels <- function(points, centroids) {
  apply(points, 1, function(row) which.min(rowSums((centroids - matrix(row, nrow = nrow(centroids), ncol = length(row), byrow = TRUE))^2)))
}
cluster <- function(points) {
  centroids <- start_centroids(points)
  labels <- rep(0L, nrow(points))
  for (step in seq_len(MAX_ITERATIONS)) {
    new_labels <- assign_labels(points, centroids)
    changed <- !identical(new_labels, labels)
    labels <- new_labels
    for (k in seq_len(SEGMENTS)) {
      members <- points[labels == k, , drop = FALSE]
      if (nrow(members) == 0L) stop("the declared segmentation collapsed to fewer segments")
      centroids[k, ] <- colMeans(members)
    }
    if (!changed) break
  }
  list(labels = labels, centroids = centroids)
}
fit <- cluster(Z); labels <- fit$labels; centroids <- fit$centroids
sizes <- sapply(seq_len(SEGMENTS), function(k) sum(labels == k))
shares <- sizes / n
sample_index <- seq(1L, n, by = STEP)
silhouette_values <- sapply(sample_index, function(i) {
  own <- labels[i]
  distances <- sqrt(rowSums((Z - matrix(Z[i, ], nrow = n, ncol = ncol(Z), byrow = TRUE))^2))
  inside <- sum(distances[labels == own]) / max(sum(labels == own) - 1L, 1L)
  outside <- min(sapply(setdiff(seq_len(SEGMENTS), own), function(k) mean(distances[labels == k])))
  (outside - inside) / max(inside, outside)
})
separation <- mean(silhouette_values)
odd <- cluster(Z[seq(1L, n, by = 2L), , drop = FALSE])
even <- cluster(Z[seq(2L, n, by = 2L), , drop = FALSE])
first <- assign_labels(Z, odd$centroids); second <- assign_labels(Z, even$centroids)
table_counts <- table(first, second)
choose2 <- function(v) v * (v - 1) / 2
index_sum <- sum(choose2(table_counts))
row_sum <- sum(choose2(rowSums(table_counts))); column_sum <- sum(choose2(colSums(table_counts)))
expected <- row_sum * column_sum / choose2(n); maximum <- (row_sum + column_sum) / 2
stability <- (index_sum - expected) / (maximum - expected)
means <- sapply(seq_len(SEGMENTS), function(k) mean(d$next_quarter_revenue_eur[labels == k]))
ratio <- max(means) / min(means)
grand <- mean(d$next_quarter_revenue_eur)
explained <- sum(sizes * (means - grand)^2) / sum((d$next_quarter_revenue_eur - grand)^2)
flags <- c(
  if (separation >= SILHOUETTE_MIN) "PASS" else "FAIL",
  if (min(shares) >= SHARE_MIN) "PASS" else "FAIL",
  if (stability >= STABILITY_MIN) "PASS" else "FAIL",
  if (ratio >= USEFULNESS_MIN) "PASS" else "FAIL"
)
verdict <- if (all(flags == "PASS")) "SEGMENTATION_READABLE_FOR_ACTION" else "DIAGNOSTIC_BLOCKS_SEGMENTATION_READING"
cat(sprintf("customers=%d\nsegments=%d\n", n, SEGMENTS))
cat(paste0("segment_sizes=", paste(sizes, collapse = ","), "\n"))
cat(paste0("segment_shares=", paste(sprintf("%.6f", shares), collapse = ","), "\n"))
cat(sprintf("smallest_share=%.6f\nbalance_flag=%s\nmean_silhouette=%.6f\nseparation_flag=%s\n",
            min(shares), flags[2], separation, flags[1]))
cat(sprintf("stability_adjusted_rand=%.6f\nstability_flag=%s\n", stability, flags[3]))
cat(paste0("segment_outcome_means=", paste(sprintf("%.6f", means), collapse = ","), "\n"))
cat(sprintf("outcome_ratio_high_low=%.6f\noutcome_variance_explained=%.6f\nusefulness_flag=%s\n", ratio, explained, flags[4]))
for (k in seq_len(SEGMENTS)) {
  cat(paste0("centroid_", k - 1L, "=", paste(sprintf("%s:%.6f", feature_names, centroids[k, ]), collapse = ","), "\n"))
}
cat(sprintf("verdict=%s\n", verdict))