# Deterministic reference engine for the 18 MSC validation dossiers. MIT; R 4.5+.

engine_version <- "msc-repro-v2"
args <- commandArgs(trailingOnly = TRUE)
if (length(args) < 1L || length(args) > 2L) stop("usage: Rscript msc-validation-reference-v1.R PAGE_ID [DATASET]")
page_id <- args[[1]]
script_arg <- commandArgs(trailingOnly = FALSE)
file_arg <- sub("^--file=", "", script_arg[grep("^--file=", script_arg)])
script_dir <- dirname(normalizePath(file_arg))
data_path <- if (length(args) == 2L) args[[2]] else file.path(dirname(script_dir), "datasets", "msc-validation-inputs-v1.csv")

fields <- c(
  "bundle_version", "page_id", "dataset_id", "record_id", "entity_id", "time_index",
  "arm_or_group", "variable", "value_number", "value_text", "value_type", "unit", "role",
  "data_status", "license"
)
required <- list(
  "MSC-P-003" = c("estimate_pp", "standard_error_pp", "z_critical"),
  "MSC-P-005" = c("loading_1", "loading_2", "loading_3", "loading_4"),
  "MSC-P-006" = c("construct_defined", "content_review_complete", "cognitive_pretest_complete", "independent_confirmation_declared", "reliability_model_declared", "invariance_scope_declared"),
  "MSC-P-008" = c("trust_satisfaction_correlation", "separation_threshold", "sample_size"),
  "MSC-P-010" = c("local_decision_value_eur", "experiment_cost_eur", "operational_risk_eur", "anticipated_gain_pp", "minimum_useful_gain_pp", "action_reversible", "assignment_feasible", "value_horizon_days", "valuation_scenario_prespecified", "risk_probability_model_declared"),
  "MSC-P-012" = c("baseline_rate", "target_rate", "alpha_two_sided", "power", "attrition_rate", "design_effect"),
  "MSC-P-014" = c("treated_geo_post_outcome", "counterfactual_post_outcome", "standard_error", "z_critical"),
  "MSC-P-015" = c("treated_pre", "treated_post", "control_pre", "control_post", "standard_error", "z_critical"),
  "MSC-P-017" = c("assigned_treatment", "assigned_control", "observed_treatment", "observed_control", "exposed_treatment", "exposed_control"),
  "MSC-P-019" = c("residual_mean", "residual_sd", "maximum_absolute_standardized_residual", "maximum_leverage", "maximum_cooks_distance", "sample_size", "parameter_count", "cluster_count", "holdout_rmse"),
  "MSC-P-023" = c("baseline_price_eur", "baseline_volume", "variable_cost_eur", "fixed_cost_change_eur", "elasticity", "price_ratio"),
  "MSC-P-028" = c("predicted_transactions_period_1", "predicted_transactions_period_2", "predicted_transactions_period_3", "mean_transaction_value_eur", "contribution_margin_rate", "discount_rate_per_period", "frequency_value_correlation", "independence_diagnostic_completed", "monetary_model_calibrated"),
  "MSC-P-029" = c(paste0("predicted_risk_", 1:8), paste0("outcome_", 1:8), "temporal_holdout_declared"),
  "MSC-P-031" = c("within_distance_1", "nearest_distance_1", "within_distance_2", "nearest_distance_2", "within_distance_3", "nearest_distance_3", "cluster_count", "initialization_seed"),
  "MSC-P-032" = c(paste0("jaccard_", 1:5), "dissolution_threshold", "resample_count"),
  "MSC-P-033" = c(paste0("train_y_", 1:5), paste0("forecast_abs_error_", 1:4), "seasonal_period", "covered_observations", "total_observations", "nominal_coverage"),
  "MSC-P-035" = c("current_spend", "previous_adstock", "adstock_decay", "saturation_alpha", "saturation_beta", "saturation_gamma"),
  "MSC-P-039" = c("mean_a", "mean_b", "sd_a", "sd_b", "n_a", "n_b", "confidence_level")
)
integer_keys <- c(
  "sample_size", "construct_defined", "content_review_complete", "cognitive_pretest_complete",
  "independent_confirmation_declared", "reliability_model_declared", "invariance_scope_declared",
  "action_reversible", "assignment_feasible", "value_horizon_days", "valuation_scenario_prespecified", "risk_probability_model_declared", "independence_diagnostic_completed", "monetary_model_calibrated", "assigned_treatment", "assigned_control",
  "observed_treatment", "observed_control", "exposed_treatment", "exposed_control", "cluster_count",
  "initialization_seed", "jaccard_intersection", "jaccard_union", "resample_count",
  "covered_observations", "total_observations", "n_a", "n_b", "parameter_count",
  "temporal_holdout_declared", "seasonal_period", paste0("outcome_", 1:8)
)
expected_units <- list(
  "MSC-P-003" = c("percentage_point", "percentage_point", "dimensionless"),
  "MSC-P-005" = rep("dimensionless", 4), "MSC-P-006" = rep("binary", 6),
  "MSC-P-008" = c("dimensionless", "dimensionless", "respondent"),
  "MSC-P-010" = c("EUR", "EUR", "EUR", "percentage_point", "percentage_point", "binary", "binary", "day", "binary", "binary"),
  "MSC-P-012" = c("probability", "probability", "probability", "probability", "probability", "dimensionless"),
  "MSC-P-014" = c("index", "index", "index", "dimensionless"),
  "MSC-P-015" = c("index", "index", "index", "index", "index", "dimensionless"),
  "MSC-P-017" = rep("unit", 6),
  "MSC-P-019" = c("outcome", "outcome", "dimensionless", "dimensionless", "dimensionless", "observation", "parameter", "cluster", "outcome"),
  "MSC-P-023" = c("EUR_per_unit", "unit", "EUR_per_unit", "EUR", "dimensionless", "ratio"),
  "MSC-P-028" = c(rep("transaction", 3), "EUR_per_transaction", "proportion", "rate per period", "dimensionless", "binary", "binary"),
  "MSC-P-029" = c(rep("probability", 8), rep("binary", 9)),
  "MSC-P-031" = c(rep("distance", 6), "cluster", "seed"),
  "MSC-P-032" = c(rep("dimensionless", 6), "resample"),
  "MSC-P-033" = c(rep("outcome", 9), "period", "observation", "observation", "probability"),
  "MSC-P-035" = c("spend", "spend", "dimensionless", "response", "dimensionless", "spend"),
  "MSC-P-039" = c("outcome", "outcome", "outcome", "outcome", "unit", "unit", "probability")
)
if (is.null(required[[page_id]])) stop("no explicit branch for page_id")

header <- strsplit(readLines(data_path, n = 1L, warn = FALSE), ",", fixed = TRUE)[[1]]
if (!identical(header, fields)) stop("dataset header must exactly match the bundle contract")
all_rows <- read.csv(data_path, stringsAsFactors = FALSE, check.names = FALSE, colClasses = "character")
keys <- paste(all_rows$page_id, all_rows$dataset_id, all_rows$record_id, all_rows$variable, sep = "|")
if (anyDuplicated(keys)) stop("duplicate bundle key")
if (any(all_rows$bundle_version != "msc-inputs-v1") || any(all_rows$license != "CC0-1.0")) stop("bundle version and CC0 license are mandatory")
populated <- nzchar(trimws(all_rows$value_number)) + nzchar(trimws(all_rows$value_text))
if (any(populated != 1L)) stop("exactly one value column must be populated")
rows <- all_rows[all_rows$page_id == page_id, , drop = FALSE]
if (!nrow(rows)) stop("no atomic slice")
if (any(rows$record_id != "protocol") || any(rows$role != "input")) stop("unexpected record or role in atomic slice")
if (any(!rows$data_status %in% c("synthetic", "parameter"))) stop("data_status must be synthetic or parameter")
if (!setequal(rows$variable, required[[page_id]]) || length(rows$variable) != length(required[[page_id]])) stop("slice variables must exactly match branch contract")
unit_contract <- setNames(expected_units[[page_id]], required[[page_id]])
if (any(rows$unit != unit_contract[rows$variable])) stop("input unit does not match the branch contract")
expected_type <- ifelse(rows$variable %in% integer_keys, "integer", "number")
if (any(rows$value_type != expected_type)) stop("input value_type does not match the branch contract")
if (any(!nzchar(trimws(rows$value_number)))) stop("v1 reference branches require numeric inputs")
numeric_values <- suppressWarnings(as.numeric(rows$value_number))
if (any(!is.finite(numeric_values))) stop("numeric inputs must be finite")
names(numeric_values) <- rows$variable
for (index in seq_along(numeric_values)) if (expected_type[[index]] == "integer" && numeric_values[[index]] != floor(numeric_values[[index]])) stop(paste(rows$variable[[index]], "must be an integer"))
x <- numeric_values
data_status <- paste(sort(unique(rows$data_status)), collapse = "+")

sha256_text <- function(value) {
  temp <- tempfile(fileext = ".txt")
  on.exit(unlink(temp), add = TRUE)
  con <- file(temp, open = "wb")
  writeChar(value, con, eos = NULL, useBytes = TRUE)
  close(con)
  if (nzchar(Sys.which("sha256sum"))) {
    # GNU coreutils escapes the line with a leading backslash whenever the
    # hashed path contains one, as every Windows temporary path does, so the
    # digest is read by pattern rather than by position.
    line <- system2("sha256sum", shQuote(temp), stdout = TRUE)[[1]]
    digest <- regmatches(line, regexpr("[0-9a-f]{64}", line))
    if (length(digest)) return(digest[[1]])
  }
  if (nzchar(Sys.which("certutil"))) {
    lines <- system2("certutil", c("-hashfile", shQuote(temp), "SHA256"), stdout = TRUE)
    candidates <- gsub(" ", "", lines[grepl("^[0-9A-Fa-f ]{64,}$", lines)])
    if (length(candidates)) return(tolower(candidates[[1]]))
  }
  stop("sha256sum or certutil is required")
}

rows <- rows[order(rows$dataset_id, rows$record_id, rows$variable), , drop = FALSE]
canonical <- paste(apply(rows[, fields, drop = FALSE], 1L, function(row) paste(sprintf("%s=%s", fields, trimws(row)), collapse = "|")), collapse = "\n")
slice_hash <- sha256_text(canonical)

metrics <- NULL
diagnostic <- NULL
if (page_id == "MSC-P-003") {
  if (x[["standard_error_pp"]] <= 0 || x[["z_critical"]] <= 0) stop("confidence interval requires positive SE and critical value")
  margin <- x[["z_critical"]] * x[["standard_error_pp"]]
  metrics <- c(estimate_pp = x[["estimate_pp"]], ci_lower_pp = x[["estimate_pp"]] - margin, ci_upper_pp = x[["estimate_pp"]] + margin)
  diagnostic <- "parameter_confidence_interval_not_prediction_or_posterior"
} else if (page_id == "MSC-P-005") {
  loadings <- x[paste0("loading_", 1:4)]
  if (any(abs(loadings) > 1)) stop("standardized loadings must be within [-1, 1]")
  metrics <- c(minimum_loading = min(loadings), maximum_loading = max(loadings), ave = mean(loadings^2))
  diagnostic <- "convergent_evidence_does_not_establish_global_validity"
} else if (page_id == "MSC-P-006") {
  if (any(!x %in% c(0, 1))) stop("scale-development gates must be binary")
  metrics <- c(completed_gates = sum(x), required_gates = length(x))
  diagnostic <- "all_prespecified_development_gates_required"
} else if (page_id == "MSC-P-008") {
  r <- x[["trust_satisfaction_correlation"]]; threshold <- x[["separation_threshold"]]; n <- x[["sample_size"]]
  if (!(r > -1 && r < 1 && threshold > 0 && threshold < 1 && n > 3)) stop("invalid correlation screening inputs")
  fisher_se <- 1 / sqrt(n - 3)
  metrics <- c(correlation = r, correlation_ci_lower = tanh(atanh(r) - 1.95996398454 * fisher_se), correlation_ci_upper = tanh(atanh(r) + 1.95996398454 * fisher_se), threshold_excess = r - threshold, screening_signal = as.numeric(r > threshold))
  diagnostic <- "correlation_screen_only_not_htmt_or_cfa_discriminant_validity"
} else if (page_id == "MSC-P-010") {
  if (min(x[c("local_decision_value_eur", "experiment_cost_eur", "operational_risk_eur", "anticipated_gain_pp", "minimum_useful_gain_pp")]) < 0 || x[["value_horizon_days"]] <= 0) stop("decision values, risks, costs and gains must be nonnegative")
  if (any(!x[c("action_reversible", "assignment_feasible", "valuation_scenario_prespecified", "risk_probability_model_declared")] %in% c(0, 1))) stop("decision screen flags must be binary")
  net <- x[["local_decision_value_eur"]] - x[["experiment_cost_eur"]] - x[["operational_risk_eur"]]
  candidate <- as.numeric(net > 0 && x[["anticipated_gain_pp"]] >= x[["minimum_useful_gain_pp"]] && x[["action_reversible"]] == 1 && x[["assignment_feasible"]] == 1 && x[["valuation_scenario_prespecified"]] == 1 && x[["risk_probability_model_declared"]] == 1)
  metrics <- c(local_net_screen_eur = net, anticipated_gain_pp = x[["anticipated_gain_pp"]], minimum_useful_gain_pp = x[["minimum_useful_gain_pp"]], passes_local_screen = candidate, value_horizon_days = x[["value_horizon_days"]])
  diagnostic <- "local_nonuniversal_governance_screen_not_voi_or_effect_estimate"
} else if (page_id == "MSC-P-012") {
  p0 <- x[["baseline_rate"]]; p1 <- x[["target_rate"]]
  if (!(p0 > 0 && p0 < 1 && p1 > 0 && p1 < 1 && p0 != p1 && x[["alpha_two_sided"]] > 0 && x[["alpha_two_sided"]] < 1 && x[["power"]] > 0.5 && x[["power"]] < 1 && x[["attrition_rate"]] >= 0 && x[["attrition_rate"]] < 1 && x[["design_effect"]] >= 1)) stop("invalid sample-size inputs")
  z_alpha <- qnorm(1 - x[["alpha_two_sided"]] / 2)
  z_power <- qnorm(x[["power"]])
  pbar <- (p0 + p1) / 2
  numerator <- z_alpha * sqrt(2 * pbar * (1 - pbar)) + z_power * sqrt(p0 * (1 - p0) + p1 * (1 - p1))
  raw <- numerator^2 / (p1 - p0)^2
  metrics <- c(raw_n_per_arm = raw, planned_n_per_arm = ceiling(raw * x[["design_effect"]] / (1 - x[["attrition_rate"]])), z_alpha_two_sided = z_alpha, z_power = z_power)
  diagnostic <- "normal_approximation_equal_allocation_quantiles_computed_from_alpha_and_power_round_up"
} else if (page_id == "MSC-P-014") {
  if (x[["standard_error"]] <= 0 || x[["z_critical"]] <= 0) stop("geo illustration requires positive supplied SE and critical value")
  effect <- x[["treated_geo_post_outcome"]] - x[["counterfactual_post_outcome"]]
  margin <- x[["z_critical"]] * x[["standard_error"]]
  metrics <- c(conditional_incremental_outcome = effect, ci_lower = effect - margin, ci_upper = effect + margin)
  diagnostic <- "illustration_with_supplied_counterfactual_and_se_not_geo_experiment_fit"
} else if (page_id == "MSC-P-015") {
  if (x[["standard_error"]] <= 0 || x[["z_critical"]] <= 0) stop("DiD illustration requires positive supplied SE and critical value")
  effect <- (x[["treated_post"]] - x[["treated_pre"]]) - (x[["control_post"]] - x[["control_pre"]])
  margin <- x[["z_critical"]] * x[["standard_error"]]
  metrics <- c(att_did = effect, ci_lower = effect - margin, ci_upper = effect + margin)
  diagnostic <- "supplied_cluster_robust_se_parallel_trends_and_no_anticipation_required"
} else if (page_id == "MSC-P-017") {
  if (min(x[c("assigned_treatment", "assigned_control")]) <= 0) stop("assigned counts must be positive")
  if (!(x[["observed_treatment"]] >= 0 && x[["observed_treatment"]] <= x[["assigned_treatment"]] && x[["observed_control"]] >= 0 && x[["observed_control"]] <= x[["assigned_control"]] && x[["exposed_treatment"]] >= 0 && x[["exposed_treatment"]] <= x[["assigned_treatment"]] && x[["exposed_control"]] >= 0 && x[["exposed_control"]] <= x[["assigned_control"]])) stop("observed and exposed counts must lie within assigned counts")
  attr_t <- 1 - x[["observed_treatment"]] / x[["assigned_treatment"]]
  attr_c <- 1 - x[["observed_control"]] / x[["assigned_control"]]
  metrics <- c(attrition_treatment = attr_t, attrition_control = attr_c, differential_attrition = attr_t - attr_c, treatment_exposure = x[["exposed_treatment"]] / x[["assigned_treatment"]], control_contamination = x[["exposed_control"]] / x[["assigned_control"]])
  diagnostic <- "diagnostics_do_not_repair_bias"
} else if (page_id == "MSC-P-019") {
  if (x[["residual_sd"]] <= 0 || x[["sample_size"]] <= x[["parameter_count"]] + 1 || x[["parameter_count"]] < 1 || x[["cluster_count"]] < 2 || min(x[c("maximum_absolute_standardized_residual", "maximum_leverage", "maximum_cooks_distance", "holdout_rmse")]) < 0 || x[["maximum_leverage"]] > 1) stop("invalid regression diagnostic inputs")
  metrics <- c(residual_mean = x[["residual_mean"]], residual_sd = x[["residual_sd"]], maximum_absolute_standardized_residual = x[["maximum_absolute_standardized_residual"]], maximum_leverage = x[["maximum_leverage"]], maximum_cooks_distance = x[["maximum_cooks_distance"]], cluster_count = x[["cluster_count"]], holdout_rmse = x[["holdout_rmse"]])
  diagnostic <- "descriptive_inspection_only_no_universal_cutoff_influence_conclusion_requires_case_deletion_refit_sensitivity"
} else if (page_id == "MSC-P-023") {
  if (x[["price_ratio"]] <= 0 || x[["baseline_price_eur"]] <= 0 || x[["baseline_volume"]] <= 0 || x[["variable_cost_eur"]] < 0 || x[["variable_cost_eur"]] >= x[["baseline_price_eur"]]) stop("invalid price scenario inputs")
  new_price <- x[["baseline_price_eur"]] * x[["price_ratio"]]
  new_volume <- x[["baseline_volume"]] * x[["price_ratio"]]^x[["elasticity"]]
  baseline_contribution <- (x[["baseline_price_eur"]] - x[["variable_cost_eur"]]) * x[["baseline_volume"]]
  scenario_contribution <- (new_price - x[["variable_cost_eur"]]) * new_volume - x[["fixed_cost_change_eur"]]
  metrics <- c(new_price_eur = new_price, conditional_volume = new_volume, conditional_volume_change = new_volume / x[["baseline_volume"]] - 1, baseline_contribution_eur = baseline_contribution, scenario_contribution_eur = scenario_contribution, incremental_contribution_eur = scenario_contribution - baseline_contribution)
  diagnostic <- "conditional_scenario_using_supplied_elasticity_not_demand_estimation_or_causal_effect"
} else if (page_id == "MSC-P-028") {
  if (min(x[c("predicted_transactions_period_1", "predicted_transactions_period_2", "predicted_transactions_period_3", "mean_transaction_value_eur")]) < 0 || !(x[["contribution_margin_rate"]] >= 0 && x[["contribution_margin_rate"]] <= 1 && x[["discount_rate_per_period"]] >= 0 && x[["discount_rate_per_period"]] < 1 && x[["frequency_value_correlation"]] >= -1 && x[["frequency_value_correlation"]] <= 1) || x[["independence_diagnostic_completed"]] != 1 || x[["monetary_model_calibrated"]] != 1) stop("invalid CLV aggregation inputs")
  contribution <- x[["mean_transaction_value_eur"]] * x[["contribution_margin_rate"]]
  flows <- sapply(1:3, function(index) x[[paste0("predicted_transactions_period_", index)]] * contribution / ((1 + x[["discount_rate_per_period"]])^index))
  metrics <- c(mean_contribution_eur = contribution, discounted_period_1_eur = flows[[1]], discounted_period_2_eur = flows[[2]], discounted_period_3_eur = flows[[3]], expected_discounted_clv_eur = sum(flows), frequency_value_correlation = x[["frequency_value_correlation"]])
  diagnostic <- "aggregation_only_after_declared_independence_diagnostic_and_monetary_calibration_low_correlation_is_not_proof"
} else if (page_id == "MSC-P-029") {
  risks <- x[paste0("predicted_risk_", 1:8)]; outcomes <- x[paste0("outcome_", 1:8)]
  if (any(risks <= 0 | risks >= 1) || any(!outcomes %in% c(0, 1)) || x[["temporal_holdout_declared"]] != 1) stop("holdout risks must be open probabilities with binary outcomes")
  brier <- mean((risks - outcomes)^2); log_loss <- -mean(outcomes * log(risks) + (1 - outcomes) * log1p(-risks))
  positives <- sum(outcomes); if (positives == 0 || positives == length(outcomes)) stop("average precision requires both outcome classes")
  true_positives <- 0; predicted_positives <- 0; average_precision <- 0
  for (threshold in sort(unique(risks), decreasing = TRUE)) { group <- outcomes[risks == threshold]; group_positives <- sum(group); true_positives <- true_positives + group_positives; predicted_positives <- predicted_positives + length(group); average_precision <- average_precision + (group_positives / positives) * (true_positives / predicted_positives) }
  gap <- mean(risks) - mean(outcomes)
  metrics <- c(sample_size = length(risks), mean_predicted_risk = mean(risks), observed_event_rate = mean(outcomes), mean_predicted_minus_observed_gap = gap, brier_score = brier, log_loss = log_loss, threshold_grouped_average_precision = average_precision)
  diagnostic <- "average_precision_is_threshold_grouped_not_generic_pr_auc_tiny_synthetic_holdout_not_deployment_validation_or_causal_effect"
} else if (page_id == "MSC-P-031") {
  silhouettes <- numeric(3)
  for (index in 1:3) { a <- x[[paste0("within_distance_", index)]]; b <- x[[paste0("nearest_distance_", index)]]; if (a < 0 || b < 0 || max(a, b) == 0) stop("per-observation silhouette distances must be nonnegative and nondegenerate"); silhouettes[[index]] <- (b - a) / max(a, b) }
  if (x[["cluster_count"]] < 2) stop("at least two clusters are required")
  metrics <- c(cluster_count_in_parent_partition = x[["cluster_count"]], sampled_observation_count = length(silhouettes), sample_mean_silhouette = mean(silhouettes), sample_minimum_silhouette = min(silhouettes), sample_maximum_silhouette = max(silhouettes), initialization_seed = x[["initialization_seed"]])
  diagnostic <- "three_sampled_points_from_larger_partition_not_full_partition_mean_stability_and_actionability_still_required"
} else if (page_id == "MSC-P-032") {
  values <- x[paste0("jaccard_", 1:5)]
  if (any(values < 0 | values > 1) || x[["dissolution_threshold"]] < 0 || x[["dissolution_threshold"]] > 1 || x[["resample_count"]] != length(values)) stop("Jaccard series must match declared resample count and [0,1] domain")
  metrics <- c(mean_jaccard = mean(values), minimum_jaccard = min(values), maximum_jaccard = max(values), dissolved_resample_share = mean(values < x[["dissolution_threshold"]]), resample_count = x[["resample_count"]])
  diagnostic <- "labels_aligned_before_clusterwise_resample_distribution_summary"
} else if (page_id == "MSC-P-033") {
  train <- x[paste0("train_y_", 1:5)]; errors <- x[paste0("forecast_abs_error_", 1:4)]; period <- as.integer(x[["seasonal_period"]])
  if (period < 1 || period >= length(train) || any(errors < 0) || x[["total_observations"]] != length(errors) || x[["covered_observations"]] < 0 || x[["covered_observations"]] > x[["total_observations"]] || x[["nominal_coverage"]] <= 0 || x[["nominal_coverage"]] >= 1) stop("invalid MASE or coverage inputs")
  naive_mae <- mean(abs(train[(period + 1):length(train)] - train[1:(length(train) - period)])); if (naive_mae <= 0) stop("in-sample seasonal naive MAE must be positive")
  forecast_mae <- mean(errors); mase <- forecast_mae / naive_mae
  coverage <- x[["covered_observations"]] / x[["total_observations"]]
  metrics <- c(forecast_mae = forecast_mae, in_sample_seasonal_naive_mae = naive_mae, mase = mase, empirical_coverage = coverage, coverage_gap = coverage - x[["nominal_coverage"]])
  diagnostic <- "mase_scaled_by_declared_in_sample_seasonal_naive_error_rolling_origin_still_required"
} else if (page_id == "MSC-P-035") {
  decay <- x[["adstock_decay"]]
  if (!(decay >= 0 && decay < 1) || min(x[c("current_spend", "previous_adstock")]) < 0 || min(x[c("saturation_alpha", "saturation_beta", "saturation_gamma")]) <= 0) stop("invalid adstock or saturation parameters")
  adstock_level <- x[["current_spend"]] + decay * x[["previous_adstock"]]
  half_life <- if (decay == 0) 0 else log(0.5) / log(decay)
  level_beta <- adstock_level^x[["saturation_beta"]]
  response <- x[["saturation_alpha"]] * level_beta / (x[["saturation_gamma"]]^x[["saturation_beta"]] + level_beta)
  metrics <- c(adstock_level = adstock_level, geometric_half_life_periods = half_life, saturated_response = response)
  diagnostic <- "single_recurrence_and_shape_illustration_transform_order_initialization_and_identification_required"
} else if (page_id == "MSC-P-039") {
  if (min(x[c("sd_a", "sd_b", "n_a", "n_b")]) <= 0 || min(x[c("n_a", "n_b")]) <= 1) stop("invalid Welch inputs")
  va <- x[["sd_a"]]^2 / x[["n_a"]]; vb <- x[["sd_b"]]^2 / x[["n_b"]]
  se <- sqrt(va + vb)
  df <- (va + vb)^2 / (va^2 / (x[["n_a"]] - 1) + vb^2 / (x[["n_b"]] - 1))
  if (x[["confidence_level"]] <= 0 || x[["confidence_level"]] >= 1) stop("confidence level must be in (0,1)")
  tcrit <- qt((1 + x[["confidence_level"]]) / 2, df)
  effect <- x[["mean_a"]] - x[["mean_b"]]
  metrics <- c(mean_difference = effect, standard_error = se, welch_df = df, student_t_critical = tcrit, ci_lower = effect - tcrit * se, ci_upper = effect + tcrit * se)
  diagnostic <- "exact_student_quantile_estimand_and_sampling_design_choose_test_not_normality_alone"
}

cat(sprintf("page_id=%s\n", page_id))
cat(sprintf("branch_id=%s-v1\n", tolower(page_id)))
cat(sprintf("engine_version=%s\n", engine_version))
cat("bundle_contract_complete=true\n")
cat(sprintf("slice_sha256=%s\n", slice_hash))
for (key in sort(names(metrics))) cat(sprintf("metric.%s=%.6f\n", key, metrics[[key]]))
cat(sprintf("diagnostic=%s\n", diagnostic))
cat(sprintf("data_status=%s\n", data_status))