#!/usr/bin/env Rscript

# ============================================================================
# Additional file 3
# Reproduce VAP evidence counts and prepare evidence-map input tables
# ============================================================================
# Manuscript:
#   Dataset-structured evidence synthesis for machine-learning prediction
#   models: a methodological framework and worked example
#
# Purpose:
#   1. Read the VAP dataset-structured reconstruction workbook.
#   2. Recalculate the evidence counts reported in Table 3.
#   3. Compare recalculated counts with the workbook's Evidence_Counts sheet.
#   4. Export node, edge, crosswalk, dataset-family, and validation-exercise
#      tables that can be used to reproduce or redesign the evidence map.
#
# Usage:
#   Rscript Additional_file_3_reproduce_evidence_counts_and_prepare_map.R \
#     VAP_dataset_structured_evidence_reconstruction_supplement.xlsx \
#     evidence_map_outputs
#
# Arguments:
#   1. Input workbook path (optional; default shown above)
#   2. Output directory (optional; default: additional_file_3_outputs)
#
# Required package:
#   readxl
#
# Notes:
#   - The workbook stores one row per report in Study_Reconstruction. Two
#     worked-example multiplicity rules are therefore declared explicitly:
#       * Schurink contributes two AUC-type analyses.
#       * Giang contributes two VAP validation tasks.
#   - A report-to-validation-exercise crosswalk is also declared explicitly,
#     because the source workbook intentionally separates report-level and
#     validation-exercise-level sheets.
#   - Temporal same-site validation is not counted as fully external.
#   - The Calvert/Faucher abstract is treated as companion evidence to the
#     Faucher preprint and is not counted as an independent study group.
# ============================================================================

options(stringsAsFactors = FALSE)

required_packages <- c("readxl")
missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0L) {
  stop(
    "Missing required R package(s): ",
    paste(missing_packages, collapse = ", "),
    ". Install with install.packages() before running this script.",
    call. = FALSE
  )
}

args <- commandArgs(trailingOnly = TRUE)
input_file <- if (length(args) >= 1L && nzchar(args[1L])) {
  args[1L]
} else {
  "VAP_dataset_structured_evidence_reconstruction_supplement.xlsx"
}
output_dir <- if (length(args) >= 2L && nzchar(args[2L])) {
  args[2L]
} else {
  "additional_file_3_outputs"
}

if (!file.exists(input_file)) {
  stop("Input workbook not found: ", normalizePath(input_file, mustWork = FALSE), call. = FALSE)
}
if (!dir.exists(output_dir)) {
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
}

clean_text <- function(x) {
  x <- as.character(x)
  x[is.na(x)] <- ""
  trimws(x)
}

read_sheet_text <- function(path, sheet_name) {
  out <- readxl::read_excel(
    path,
    sheet = sheet_name,
    col_types = "text",
    .name_repair = "minimal"
  )
  out <- as.data.frame(out, stringsAsFactors = FALSE, check.names = FALSE)
  out[] <- lapply(out, clean_text)
  out
}

require_columns <- function(data, required, sheet_name) {
  missing <- setdiff(required, names(data))
  if (length(missing) > 0L) {
    stop(
      "Sheet '", sheet_name, "' is missing required column(s): ",
      paste(missing, collapse = ", "),
      call. = FALSE
    )
  }
}

is_yes <- function(x) {
  tolower(clean_text(x)) %in% c("yes", "y", "true", "1")
}

as_integer_safe <- function(x, field_name) {
  value <- suppressWarnings(as.integer(clean_text(x)))
  if (any(is.na(value))) {
    stop("Non-integer value found in field '", field_name, "'.", call. = FALSE)
  }
  value
}

normalise_count <- function(x) {
  x <- ifelse(is.na(x), "", as.character(x))
  x <- gsub("[[:space:]]+", " ", x)
  trimws(x)
}

make_node_id <- function(prefix, label) {
  slug <- tolower(label)
  slug <- gsub("[^a-z0-9]+", "_", slug)
  slug <- gsub("^_+|_+$", "", slug)
  paste0(prefix, ":", slug)
}

count_delimited <- function(x, delimiter = ";") {
  x <- clean_text(x)
  vapply(
    x,
    function(value) {
      if (!nzchar(value)) return(0L)
      parts <- clean_text(strsplit(value, delimiter, fixed = TRUE)[[1L]])
      sum(nzchar(parts))
    },
    integer(1)
  )
}

required_sheets <- c(
  "Study_Reconstruction",
  "Validation_Exercises",
  "Dataset_Map",
  "Evidence_Counts",
  "Synthesis_Decision",
  "Coding_Legend"
)
available_sheets <- readxl::excel_sheets(input_file)
missing_sheets <- setdiff(required_sheets, available_sheets)
if (length(missing_sheets) > 0L) {
  stop(
    "Input workbook is missing required sheet(s): ",
    paste(missing_sheets, collapse = ", "),
    call. = FALSE
  )
}

study <- read_sheet_text(input_file, "Study_Reconstruction")
validation <- read_sheet_text(input_file, "Validation_Exercises")
dataset_map <- read_sheet_text(input_file, "Dataset_Map")
expected_counts <- read_sheet_text(input_file, "Evidence_Counts")
synthesis_decision <- read_sheet_text(input_file, "Synthesis_Decision")
coding_legend <- read_sheet_text(input_file, "Coding_Legend")

require_columns(
  study,
  c(
    "Report_ID", "Report", "Data_source_family", "Analytic_dataset",
    "Validation_design", "Models/pipeline", "Main_performance_summary",
    "Calibration_ready", "Decision_analytic", "Dependence_notes",
    "Synthesis_handling"
  ),
  "Study_Reconstruction"
)
require_columns(
  validation,
  c(
    "VE_ID", "Evidence_unit", "Source_family", "Validation_design",
    "AUROC_eligible", "External_validation_count", "Calibration_ready",
    "Decision_analytic_ready", "Notes"
  ),
  "Validation_Exercises"
)
require_columns(
  dataset_map,
  c(
    "Source_family", "Reports", "Analytic_datasets", "Validation_exercises",
    "Reused_source", "Primary_interpretation"
  ),
  "Dataset_Map"
)
require_columns(
  expected_counts,
  c("Evidence_descriptor", "Conventional_count", "Dataset_structured_count"),
  "Evidence_Counts"
)

# ---------------------------------------------------------------------------
# Worked-example multiplicity and crosswalk rules
# ---------------------------------------------------------------------------
# These rules make task-level counts explicit where the workbook deliberately
# stores only a report-level row.
report_auc_multiplicity <- data.frame(
  Report_ID = sprintf("S%02d", 1:10),
  Reported_AUROC_summaries = c(2L, 0L, 0L, 1L, 1L, 1L, 2L, 1L, 1L, 1L),
  stringsAsFactors = FALSE
)

# Reports that contribute dependent evidence because they contain multiple
# AUC-type analyses/tasks or belong to the reused MIMIC-III/companion family.
dependent_report_groups <- data.frame(
  Report_ID = c("S01", "S07", "S08", "S09", "S10"),
  Study_group = c("Schurink", "Giang", "Liang", "Faucher/Calvert", "Faucher/Calvert"),
  stringsAsFactors = FALSE
)

# Validation exercise to report crosswalk.
ve_report_crosswalk <- data.frame(
  VE_ID = sprintf("VE%02d", 1:11),
  Report_ID = c(
    "S01", "S02", "S03", "S04", "S05", "S06",
    "S07", "S07", "S08", "S09", "S10"
  ),
  stringsAsFactors = FALSE
)

# Source names differ slightly between sheets; map them to the canonical labels
# used in Dataset_Map.
source_alias <- data.frame(
  Validation_source = c(
    "UMCU", "NTDB", "OUTCOMEREA", "Hospital Tacchini",
    "Qatar/Hamad registry", "LifeCenter", "MIMIC-III family"
  ),
  Canonical_source = c(
    "UMCU", "National Trauma Data Bank", "OUTCOMEREA", "Hospital Tacchini",
    "Qatar/Hamad trauma registry", "LifeCenter Hospital", "MIMIC-III family"
  ),
  stringsAsFactors = FALSE
)

if (!all(report_auc_multiplicity$Report_ID %in% study$Report_ID)) {
  stop("The AUC multiplicity table contains Report_ID values absent from the workbook.", call. = FALSE)
}
if (!all(dependent_report_groups$Report_ID %in% study$Report_ID)) {
  stop("The dependent-report table contains Report_ID values absent from the workbook.", call. = FALSE)
}
if (!all(ve_report_crosswalk$VE_ID %in% validation$VE_ID)) {
  stop("The validation crosswalk contains VE_ID values absent from the workbook.", call. = FALSE)
}
if (!all(ve_report_crosswalk$Report_ID %in% study$Report_ID)) {
  stop("The validation crosswalk contains Report_ID values absent from the workbook.", call. = FALSE)
}

# ---------------------------------------------------------------------------
# Reproduce Table 3 evidence counts
# ---------------------------------------------------------------------------
external_count <- sum(
  as_integer_safe(validation$External_validation_count, "External_validation_count")
)
reused_rows <- dataset_map[is_yes(dataset_map$Reused_source), , drop = FALSE]
reused_labels <- sub(
  "[[:space:]]+family$",
  "",
  reused_rows$Source_family,
  ignore.case = TRUE
)

reproduced_counts <- data.frame(
  Evidence_descriptor = c(
    "Included reports",
    "Reported primary model entries",
    "Reported AUROC/c-statistic summaries",
    "Unique data-source/cohort groups",
    "Unique analytic datasets",
    "Unique validation exercises",
    "Clearly fully external validations",
    "Internal, temporal, or unclear validations",
    "Calibration-ready validation exercises",
    "Decision-analytic validation exercises",
    "Reused or overlapping dataset families",
    "Reports contributing dependent estimates"
  ),
  Conventional_count = c(
    as.character(nrow(study)),
    as.character(nrow(study)),
    as.character(sum(report_auc_multiplicity$Reported_AUROC_summaries)),
    "", "", "", "", "", "", "", "", ""
  ),
  Dataset_structured_count = c(
    "", "", "",
    as.character(nrow(dataset_map)),
    as.character(sum(as_integer_safe(dataset_map$Analytic_datasets, "Analytic_datasets"))),
    as.character(nrow(validation)),
    as.character(external_count),
    as.character(nrow(validation) - external_count),
    as.character(sum(is_yes(validation$Calibration_ready))),
    as.character(sum(is_yes(validation$Decision_analytic_ready))),
    paste0(length(reused_labels), ": ", paste(reused_labels, collapse = "; ")),
    paste0(
      nrow(dependent_report_groups),
      " reports / ",
      length(unique(dependent_report_groups$Study_group)),
      " study groups after companion merge"
    )
  ),
  stringsAsFactors = FALSE
)

expected_aligned <- expected_counts[
  match(reproduced_counts$Evidence_descriptor, expected_counts$Evidence_descriptor),
  ,
  drop = FALSE
]
if (any(is.na(match(reproduced_counts$Evidence_descriptor, expected_counts$Evidence_descriptor)))) {
  stop("One or more reproduced evidence descriptors are absent from Evidence_Counts.", call. = FALSE)
}

count_check <- data.frame(
  Evidence_descriptor = reproduced_counts$Evidence_descriptor,
  Reproduced_conventional = reproduced_counts$Conventional_count,
  Workbook_conventional = expected_aligned$Conventional_count,
  Reproduced_dataset_structured = reproduced_counts$Dataset_structured_count,
  Workbook_dataset_structured = expected_aligned$Dataset_structured_count,
  Match = (
    normalise_count(reproduced_counts$Conventional_count) ==
      normalise_count(expected_aligned$Conventional_count)
  ) & (
    normalise_count(reproduced_counts$Dataset_structured_count) ==
      normalise_count(expected_aligned$Dataset_structured_count)
  ),
  stringsAsFactors = FALSE
)

# ---------------------------------------------------------------------------
# Prepare crosswalks and evidence-map nodes
# ---------------------------------------------------------------------------
validation$Report_ID <- ve_report_crosswalk$Report_ID[
  match(validation$VE_ID, ve_report_crosswalk$VE_ID)
]
validation$Report <- study$Report[match(validation$Report_ID, study$Report_ID)]
validation$Canonical_source <- source_alias$Canonical_source[
  match(validation$Source_family, source_alias$Validation_source)
]
if (any(!nzchar(validation$Canonical_source))) {
  stop("At least one validation-exercise source could not be mapped to Dataset_Map.", call. = FALSE)
}

source_nodes <- data.frame(
  node_id = make_node_id("source", dataset_map$Source_family),
  label = dataset_map$Source_family,
  node_type = "Data source/cohort",
  layer = 2L,
  source_family = dataset_map$Source_family,
  validation_design = "",
  reused_source = dataset_map$Reused_source,
  external_validation = "",
  calibration_ready = "",
  decision_analytic_ready = "",
  notes = dataset_map$Primary_interpretation,
  stringsAsFactors = FALSE
)
source_id_lookup <- setNames(source_nodes$node_id, source_nodes$label)

report_reuse <- dataset_map$Reused_source[
  match(study$Data_source_family, dataset_map$Source_family)
]
report_nodes <- data.frame(
  node_id = paste0("report:", study$Report_ID),
  label = study$Report,
  node_type = "Publication/report",
  layer = 1L,
  source_family = study$Data_source_family,
  validation_design = study$Validation_design,
  reused_source = report_reuse,
  external_validation = "",
  calibration_ready = study$Calibration_ready,
  decision_analytic_ready = study$Decision_analytic,
  notes = study$Dependence_notes,
  stringsAsFactors = FALSE
)

validation_nodes <- data.frame(
  node_id = paste0("validation:", validation$VE_ID),
  label = validation$Evidence_unit,
  node_type = "Validation exercise",
  layer = 3L,
  source_family = validation$Canonical_source,
  validation_design = validation$Validation_design,
  reused_source = dataset_map$Reused_source[
    match(validation$Canonical_source, dataset_map$Source_family)
  ],
  external_validation = validation$External_validation_count,
  calibration_ready = validation$Calibration_ready,
  decision_analytic_ready = validation$Decision_analytic_ready,
  notes = validation$Notes,
  stringsAsFactors = FALSE
)

domain_nodes <- data.frame(
  node_id = c(
    "domain:discrimination",
    "domain:other_non_auc",
    "domain:calibration",
    "domain:decision_analytic"
  ),
  label = c(
    "Discrimination (AUROC/c-statistic)",
    "Other or non-AUROC performance",
    "Calibration",
    "Decision-analytic evidence"
  ),
  node_type = "Performance domain",
  layer = 4L,
  source_family = "",
  validation_design = "",
  reused_source = "",
  external_validation = "",
  calibration_ready = "",
  decision_analytic_ready = "",
  notes = c(
    "Ranking performance",
    "Dynamic, accuracy, Gini, RMSE, or other non-AUROC evidence",
    "Absolute risk accuracy",
    "Net benefit or related clinical decision evidence"
  ),
  stringsAsFactors = FALSE
)

nodes <- rbind(report_nodes, source_nodes, validation_nodes, domain_nodes)

# ---------------------------------------------------------------------------
# Prepare evidence-map edges
# ---------------------------------------------------------------------------
report_source_edges <- data.frame(
  from = paste0("report:", study$Report_ID),
  to = unname(source_id_lookup[study$Data_source_family]),
  relation = "uses_data_source",
  weight = 1,
  note = study$Dependence_notes,
  stringsAsFactors = FALSE
)

source_validation_edges <- data.frame(
  from = unname(source_id_lookup[validation$Canonical_source]),
  to = paste0("validation:", validation$VE_ID),
  relation = "contains_validation_exercise",
  weight = 1,
  note = validation$Validation_design,
  stringsAsFactors = FALSE
)

report_validation_edges <- data.frame(
  from = paste0("report:", validation$Report_ID),
  to = paste0("validation:", validation$VE_ID),
  relation = "reports_validation_exercise",
  weight = 1,
  note = validation$Notes,
  stringsAsFactors = FALSE
)

performance_edges <- vector("list", nrow(validation))
for (i in seq_len(nrow(validation))) {
  targets <- if (is_yes(validation$AUROC_eligible[i])) {
    "domain:discrimination"
  } else {
    "domain:other_non_auc"
  }
  if (is_yes(validation$Calibration_ready[i])) {
    targets <- c(targets, "domain:calibration")
  }
  if (is_yes(validation$Decision_analytic_ready[i])) {
    targets <- c(targets, "domain:decision_analytic")
  }
  performance_edges[[i]] <- data.frame(
    from = paste0("validation:", validation$VE_ID[i]),
    to = targets,
    relation = "reports_performance_domain",
    weight = 1,
    note = validation$Notes[i],
    stringsAsFactors = FALSE
  )
}
performance_edges <- do.call(rbind, performance_edges)

companion_edge <- data.frame(
  from = "report:S10",
  to = "report:S09",
  relation = "companion_of",
  weight = 1,
  note = "Calvert/Faucher abstract treated as companion evidence to Faucher preprint",
  stringsAsFactors = FALSE
)

edges <- rbind(
  report_source_edges,
  source_validation_edges,
  report_validation_edges,
  performance_edges,
  companion_edge
)
edges$edge_id <- sprintf("E%03d", seq_len(nrow(edges)))
edges <- edges[, c("edge_id", "from", "to", "relation", "weight", "note")]

# ---------------------------------------------------------------------------
# Additional analysis-ready summaries
# ---------------------------------------------------------------------------
dataset_family_summary <- data.frame(
  source_id = source_nodes$node_id,
  Source_family = dataset_map$Source_family,
  Report_count = count_delimited(dataset_map$Reports, ";"),
  Reports = dataset_map$Reports,
  Analytic_datasets = as_integer_safe(dataset_map$Analytic_datasets, "Analytic_datasets"),
  Validation_exercises = as_integer_safe(dataset_map$Validation_exercises, "Validation_exercises"),
  Reused_source = dataset_map$Reused_source,
  Primary_interpretation = dataset_map$Primary_interpretation,
  stringsAsFactors = FALSE
)

validation_exercise_summary <- data.frame(
  VE_ID = validation$VE_ID,
  Report_ID = validation$Report_ID,
  Report = validation$Report,
  Evidence_unit = validation$Evidence_unit,
  Canonical_source = validation$Canonical_source,
  Validation_design = validation$Validation_design,
  AUROC_eligible = validation$AUROC_eligible,
  External_validation_count = validation$External_validation_count,
  Calibration_ready = validation$Calibration_ready,
  Decision_analytic_ready = validation$Decision_analytic_ready,
  Notes = validation$Notes,
  stringsAsFactors = FALSE
)

# ---------------------------------------------------------------------------
# Write outputs
# ---------------------------------------------------------------------------
write_csv <- function(data, filename) {
  utils::write.csv(
    data,
    file = file.path(output_dir, filename),
    row.names = FALSE,
    na = "",
    fileEncoding = "UTF-8"
  )
}

write_csv(reproduced_counts, "evidence_counts_reproduced.csv")
write_csv(count_check, "evidence_counts_check.csv")
write_csv(nodes, "evidence_map_nodes.csv")
write_csv(edges, "evidence_map_edges.csv")
write_csv(ve_report_crosswalk, "report_validation_crosswalk.csv")
write_csv(dataset_family_summary, "dataset_family_summary.csv")
write_csv(validation_exercise_summary, "validation_exercise_summary.csv")
write_csv(synthesis_decision, "synthesis_decision.csv")
write_csv(coding_legend, "coding_legend.csv")

readme_lines <- c(
  "Additional file 3 outputs",
  "=========================",
  paste0("Generated: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")),
  paste0("Input workbook: ", normalizePath(input_file, winslash = "/", mustWork = TRUE)),
  "",
  "Output files:",
  "- evidence_counts_reproduced.csv: recalculated Table 3 counts.",
  "- evidence_counts_check.csv: comparison with the workbook Evidence_Counts sheet.",
  "- evidence_map_nodes.csv: nodes for report -> data source -> validation exercise -> performance-domain mapping.",
  "- evidence_map_edges.csv: typed edges for Sankey, alluvial, or network visualization.",
  "- report_validation_crosswalk.csv: explicit report-to-validation-exercise mapping.",
  "- dataset_family_summary.csv: one row per source family.",
  "- validation_exercise_summary.csv: one row per validation exercise with canonical source and report link.",
  "- synthesis_decision.csv: worked-example synthesis decisions.",
  "- coding_legend.csv: operational definitions used in the workbook.",
  "",
  paste0("Count check passed: ", ifelse(all(count_check$Match), "YES", "NO")),
  "",
  "Interpretive rules:",
  "- Temporal same-site validation is not counted as fully external.",
  "- Repeated MIMIC-III reports are mapped to one reused source family.",
  "- The Calvert/Faucher abstract is companion evidence to the Faucher preprint.",
  "- Calibration-ready evidence requires extractable or synthesis-ready calibration information.",
  "- Evidence-map edge types can be filtered according to the intended visualization."
)
writeLines(readme_lines, con = file.path(output_dir, "README_outputs.txt"), useBytes = TRUE)

capture.output(
  sessionInfo(),
  file = file.path(output_dir, "R_session_info.txt")
)

if (!all(count_check$Match)) {
  warning(
    "One or more recalculated counts do not match the workbook Evidence_Counts sheet. ",
    "Inspect evidence_counts_check.csv.",
    call. = FALSE
  )
}

cat("Completed successfully.\n")
cat("Input workbook: ", normalizePath(input_file, winslash = "/", mustWork = TRUE), "\n", sep = "")
cat("Output directory: ", normalizePath(output_dir, winslash = "/", mustWork = TRUE), "\n", sep = "")
cat("Evidence-count check passed: ", ifelse(all(count_check$Match), "YES", "NO"), "\n", sep = "")