# Reproduction code for the Monte Carlo study in
# "Degree Distributions in Psychological Networks"
#
# This single script generates the numerical results used in the simulation
# section, including the conditional time to complete failure required for
# Table 3. It also writes the manuscript figures based on the simulations.
#
# Required packages: igraph and ggplot2.
# Run from the project directory with:
#   Rscript Manuscript_Simulation_Reproduction_CompleteFailure_20260817.R

library(igraph)
library(parallel)
library(ggplot2)

command_line <- commandArgs(trailingOnly = FALSE)
script_argument <- grep("^--file=", command_line, value = TRUE)
script_path <- if (length(script_argument)) {
  normalizePath(
    sub("^--file=", "", script_argument[1L]),
    winslash = "/",
    mustWork = TRUE
  )
} else {
  normalizePath(
    "Manuscript_Simulation_Reproduction_CompleteFailure_20260817.R",
    winslash = "/",
    mustWork = TRUE
  )
}
project_directory <- dirname(script_path)

# =============================================================================
# 1. USER CONFIGURATION
# =============================================================================

# Enter 1 for the primary manuscript design.
# Enter 2 for the otherwise identical beta = .05 Watts-Strogatz sensitivity design.
STUDY_DESIGN <- 1L

# Enter 0 to run all four threshold conditions. Enter 1, 2, 3, or 4 to run
# only Sensitive, Moderate, Elevated, or Resilient thresholds, respectively.
# Thus, entering 1 runs the first condition of the study.
THRESHOLD_CONDITION <- 0L

# Leave as NA to use the planned number of replications. Set to a small integer
# (for example, 10) for a quick code check.
REPLICATIONS_OVERRIDE <- NA_integer_

# Results are independent of worker count because every replication receives
# its own deterministic seed.
N_WORKERS <- 3L

# Leave as NULL to use the study-specific output folder defined below.
OUTPUT_DIRECTORY_OVERRIDE <- NULL

# The replication-level file is always saved because it is the auditable source
# for every reported mean, standard deviation, percentage, and paired contrast.
MAKE_FIGURES <- TRUE
MAKE_REPRESENTATIVE_TOPOLOGY_FIGURE <- TRUE

study_conditions <- list(
  `1` = list(
    label = "Primary manuscript study: beta = .10",
    n_nodes = 20L,
    edge_counts = c(Sparse = 20L, Denser = 40L),
    ws_rewiring = 0.10,
    n_replications = 1000L,
    base_seed = 3887891L,
    output_directory =
      "Output_Manuscript_Reproduction_CompleteFailure_Study1_20260817"
  ),
  `2` = list(
    label = "Watts-Strogatz sensitivity study: beta = .05",
    n_nodes = 20L,
    edge_counts = c(Sparse = 20L, Denser = 40L),
    ws_rewiring = 0.05,
    n_replications = 1000L,
    base_seed = 3887891L,
    output_directory =
      "Output_Manuscript_Reproduction_CompleteFailure_Study2_20260817"
  )
)

study_key <- as.character(as.integer(STUDY_DESIGN))
if (!study_key %in% names(study_conditions)) {
  stop("STUDY_DESIGN must be 1 or 2.")
}
design <- study_conditions[[study_key]]

if (!is.na(REPLICATIONS_OVERRIDE)) {
  design$n_replications <- as.integer(REPLICATIONS_OVERRIDE)
}
if (design$n_replications < 1L) stop("At least one replication is required.")

output_directory_name <- if (THRESHOLD_CONDITION == 0L) {
  design$output_directory
} else {
  paste0(design$output_directory, "_Threshold", THRESHOLD_CONDITION)
}
if (!is.na(REPLICATIONS_OVERRIDE)) {
  output_directory_name <- paste0(
    output_directory_name,
    "_Smoke_n",
    design$n_replications
  )
}

output_directory <- if (is.null(OUTPUT_DIRECTORY_OVERRIDE)) {
  file.path(project_directory, output_directory_name)
} else {
  requested_output <- as.character(OUTPUT_DIRECTORY_OVERRIDE)
  if (grepl("^([A-Za-z]:[/\\\\]|/)", requested_output)) {
    requested_output
  } else {
    file.path(project_directory, requested_output)
  }
}

metrics <- c("degree", "strength", "closeness", "betweenness")
target_order <- c("Low", "Average", "High")
graph_keys <- c("S_ER", "S_BA", "S_WS", "D_ER", "D_BA", "D_WS")
threshold_ranges <- list(
  Sensitive = c(0.10, 0.20),
  Moderate = c(0.20, 0.40),
  Elevated = c(0.40, 0.60),
  Resilient = c(0.60, 0.80)
)
if (!THRESHOLD_CONDITION %in% 0:4) {
  stop("THRESHOLD_CONDITION must be 0, 1, 2, 3, or 4.")
}
active_threshold_ranges <- if (THRESHOLD_CONDITION == 0L) {
  threshold_ranges
} else {
  threshold_ranges[as.integer(THRESHOLD_CONDITION)]
}
threshold_codes <- vapply(
  active_threshold_ranges,
  function(x) sprintf("%.2f-%.2f", x[1L], x[2L]),
  character(1L)
)
threshold_display <- setNames(
  sprintf("%s (%s)", names(active_threshold_ranges), threshold_codes),
  threshold_codes
)

network_names <- c(
  ER = "Erdos-Renyi",
  BA = "Barabasi-Albert",
  WS = "Watts-Strogatz"
)
density_names <- c(S = "Sparse", D = "Denser")

# The simulated partial correlations are positive and uniformly distributed.
weight_lower <- 0.10
weight_upper <- 0.20
minimum_precision_eigenvalue <- 0.05
exclude_isolates_as_targets <- TRUE

# Explicit version and RNG settings reproduce the archived R 4.4.2 stream.
RNGversion("4.4.2")
RNGkind("Mersenne-Twister", "Inversion", "Rejection")

expected_versions <- c(R = "4.4.2", igraph = "2.1.1")
if (as.character(getRversion()) != expected_versions[["R"]]) {
  warning("Exact replay was validated with R 4.4.2; current R is ", getRversion(), ".")
}
if (as.character(packageVersion("igraph")) != expected_versions[["igraph"]]) {
  warning(
    "Exact replay was validated with igraph 2.1.1; current igraph is ",
    packageVersion("igraph"),
    "."
  )
}

stopifnot(
  design$n_nodes == 20L,
  identical(names(design$edge_counts), c("Sparse", "Denser")),
  all(design$edge_counts >= design$n_nodes - 1L),
  all(design$edge_counts <= choose(design$n_nodes, 2L)),
  N_WORKERS >= 1L
)

dir.create(output_directory, recursive = TRUE, showWarnings = FALSE)

# =============================================================================
# 2. GRAPH GENERATION AND EDGE WEIGHTS
# =============================================================================

density_from_key <- function(key) {
  density_names[[strsplit(key, "_", fixed = TRUE)[[1L]][1L]]]
}

network_from_key <- function(key) {
  network_names[[strsplit(key, "_", fixed = TRUE)[[1L]][2L]]]
}

adjust_edge_count <- function(graph, target_edges) {
  graph <- simplify(graph, remove.multiple = TRUE, remove.loops = TRUE)
  target_edges <- as.integer(target_edges)

  while (gsize(graph) > target_edges) {
    removable <- setdiff(seq_len(gsize(graph)), as.integer(bridges(graph)))
    if (!length(removable)) {
      stop("No non-bridge edge was available for edge-count matching.")
    }
    graph <- delete_edges(
      graph,
      removable[sample.int(length(removable), 1L)]
    )
  }

  while (gsize(graph) < target_edges) {
    missing_edges <- as_edgelist(
      complementer(graph, loops = FALSE),
      names = FALSE
    )
    selected <- missing_edges[sample.int(nrow(missing_edges), 1L), ]
    graph <- add_edges(graph, as.integer(selected))
  }

  graph
}

generate_topologies <- function(n_nodes, target_edges, ws_rewiring) {
  target_edges <- as.integer(target_edges)

  # Barabasi-Albert graphs are generated first to preserve the study's random
  # number sequence, then matched to the requested edge count.
  edges_per_new_vertex <- max(1L, as.integer(round(target_edges / n_nodes)))
  graph_ba <- sample_pa(
    n = n_nodes,
    power = 1,
    m = edges_per_new_vertex,
    directed = FALSE
  )
  graph_ba <- simplify(graph_ba, remove.multiple = TRUE, remove.loops = TRUE)
  ba_edges_before_matching <- gsize(graph_ba)
  graph_ba <- adjust_edge_count(graph_ba, target_edges)
  graph_ba <- set_graph_attr(
    graph_ba,
    "edges_before_matching",
    value = ba_edges_before_matching
  )

  graph_er <- sample_gnm(
    n = n_nodes,
    m = target_edges,
    directed = FALSE,
    loops = FALSE
  )
  graph_er <- set_graph_attr(
    graph_er,
    "edges_before_matching",
    value = gsize(graph_er)
  )

  lattice_neighbors <- max(
    1L,
    as.integer(round(target_edges / n_nodes))
  )
  graph_ws <- sample_smallworld(
    dim = 1,
    size = n_nodes,
    nei = lattice_neighbors,
    p = ws_rewiring,
    loops = FALSE,
    multiple = FALSE
  )
  ws_edges_before_matching <- gsize(graph_ws)
  if (ws_edges_before_matching != target_edges) {
    stop("The configured Watts-Strogatz graph did not have the target edge count.")
  }
  graph_ws <- set_graph_attr(
    graph_ws,
    "edges_before_matching",
    value = ws_edges_before_matching
  )

  stopifnot(
    gsize(graph_er) == target_edges,
    gsize(graph_ba) == target_edges,
    gsize(graph_ws) == target_edges
  )

  # The return order matches the original analysis after all three graph
  # generators have consumed their random numbers.
  list(ER = graph_er, BA = graph_ba, WS = graph_ws)
}

add_partial_correlation_weights <- function(
    graph,
    lower = 0.10,
    upper = 0.20,
    minimum_eigenvalue = 0.05
) {
  edge_list <- as_edgelist(graph, names = FALSE)
  n_edges <- nrow(edge_list)
  magnitudes <- runif(n_edges, lower, upper)

  off_diagonal <- matrix(0, vcount(graph), vcount(graph))
  if (n_edges) {
    off_diagonal[edge_list] <- -magnitudes
    off_diagonal[cbind(edge_list[, 2L], edge_list[, 1L])] <- -magnitudes
  }

  minimum_off_diagonal_eigenvalue <- min(
    eigen(off_diagonal, symmetric = TRUE, only.values = TRUE)$values
  )
  scale_factor <- 1
  if (1 + minimum_off_diagonal_eigenvalue < minimum_eigenvalue) {
    scale_factor <-
      (1 - minimum_eigenvalue) / abs(minimum_off_diagonal_eigenvalue)
  }

  precision <- diag(vcount(graph)) + scale_factor * off_diagonal
  diagonal_scale <- 1 / sqrt(diag(precision))
  partial_correlation <-
    -outer(diagonal_scale, diagonal_scale) * precision
  diag(partial_correlation) <- 1

  if (n_edges) {
    E(graph)$pcor <- partial_correlation[edge_list]
    E(graph)$weight <- abs(E(graph)$pcor)
  }
  graph <- set_graph_attr(graph, "pcor_scale_factor", value = scale_factor)
  graph
}

prepare_graphs <- function(configuration) {
  density_codes <- c(Sparse = "S", Denser = "D")
  graphs <- unlist(lapply(names(configuration$edge_counts), function(density) {
    generated <- generate_topologies(
      configuration$n_nodes,
      configuration$edge_counts[[density]],
      configuration$ws_rewiring
    )
    names(generated) <- paste(density_codes[[density]], names(generated), sep = "_")
    generated
  }), recursive = FALSE)

  graphs <- lapply(graphs, add_partial_correlation_weights,
    lower = configuration$weight_lower,
    upper = configuration$weight_upper,
    minimum_eigenvalue = configuration$minimum_precision_eigenvalue
  )
  graphs[configuration$graph_keys]
}

# =============================================================================
# 3. CENTRALITY, TARGET SELECTION, AND CASCADE PROCESS
# =============================================================================

compute_centrality <- function(graph, metric) {
  coupling <- abs(E(graph)$pcor)
  distance <- 1 / pmax(coupling, .Machine$double.eps)

  value <- switch(
    metric,
    degree = degree(graph, mode = "all", loops = FALSE),
    strength = strength(
      graph,
      mode = "all",
      loops = FALSE,
      weights = coupling
    ),
    closeness = harmonic_centrality(
      graph,
      mode = "all",
      weights = distance,
      normalized = TRUE
    ),
    betweenness = betweenness(
      graph,
      directed = FALSE,
      weights = distance,
      normalized = TRUE
    ),
    stop("Unknown centrality metric.")
  )

  value[!is.finite(value)] <- 0
  as.numeric(value)
}

select_targets <- function(graph, centrality, exclude_isolates = TRUE) {
  eligible <- seq_len(vcount(graph))
  if (exclude_isolates) {
    nonisolates <- which(degree(graph, loops = FALSE) > 0)
    if (length(nonisolates)) eligible <- nonisolates
  }

  values <- centrality[eligible]
  if (!length(values) || any(!is.finite(values))) {
    stop("No valid target candidates were available.")
  }

  select_one <- function(candidates) {
    candidates[sample.int(length(candidates), 1L)]
  }

  high_candidates <- eligible[values == max(values)]
  low_candidates <- eligible[values == min(values)]
  distance_from_mean <- abs(values - mean(values))
  average_candidates <- eligible[distance_from_mean == min(distance_from_mean)]

  # This draw order is retained for exact reproduction of the original study.
  high_target <- select_one(high_candidates)
  low_target <- select_one(low_candidates)
  average_target <- select_one(average_candidates)

  list(
    targets = c(
      High = high_target,
      Average = average_target,
      Low = low_target
    ),
    tie_counts = c(
      High = length(high_candidates),
      Average = length(average_candidates),
      Low = length(low_candidates)
    ),
    centrality_sd = stats::sd(values),
    eligible_n = length(eligible)
  )
}

run_cascade <- function(graph, target, thresholds) {
  adjacency <- as.matrix(
    as_adjacency_matrix(graph, attr = "weight", sparse = FALSE)
  )
  incident_weight <- rowSums(adjacency)
  state <- integer(vcount(graph))
  state[target] <- 1L
  state_changing_steps <- 0L

  for (step in seq_len(vcount(graph))) {
    active_incident_weight <- as.numeric(adjacency %*% state)
    active_fraction <- numeric(vcount(graph))
    defined <- incident_weight > 0
    active_fraction[defined] <-
      active_incident_weight[defined] / incident_weight[defined]

    updated_state <- pmax(
      state,
      as.integer(active_fraction >= thresholds)
    )
    if (identical(updated_state, state)) break
    state <- updated_state
    state_changing_steps <- state_changing_steps + 1L
  }

  c(
    final_active_count = sum(state),
    final_active = mean(state),
    steps_to_stopping_point = state_changing_steps
  )
}

topology_diagnostics <- function(graph, replication, graph_key) {
  graph_degree <- degree(graph, loops = FALSE)
  graph_components <- components(graph, mode = "weak")
  edges_before <- as.integer(graph_attr(graph, "edges_before_matching"))

  data.frame(
    replication = replication,
    graph_key = graph_key,
    density_condition = density_from_key(graph_key),
    network = network_from_key(graph_key),
    n_nodes = vcount(graph),
    n_edges = gsize(graph),
    density = edge_density(graph, loops = FALSE),
    edges_before_matching = edges_before,
    edges_added = max(0L, gsize(graph) - edges_before),
    edges_removed = max(0L, edges_before - gsize(graph)),
    mean_degree = mean(graph_degree),
    degree_variance = stats::var(graph_degree),
    maximum_degree = max(graph_degree),
    global_clustering = transitivity(
      graph,
      type = "global",
      isolates = "zero"
    ),
    mean_path_length = mean_distance(
      graph,
      directed = FALSE,
      unconnected = TRUE,
      weights = NA
    ),
    n_components = graph_components$no,
    n_isolates = sum(graph_degree == 0),
    largest_component = max(graph_components$csize),
    connected = as.integer(graph_components$no == 1L),
    stringsAsFactors = FALSE
  )
}

weight_diagnostics <- function(graph, replication, graph_key) {
  weights <- E(graph)$weight
  data.frame(
    replication = replication,
    graph_key = graph_key,
    density_condition = density_from_key(graph_key),
    network = network_from_key(graph_key),
    n_weights = length(weights),
    weight_sum = sum(weights),
    weight_sum_squares = sum(weights^2),
    minimum_weight = min(weights),
    maximum_weight = max(weights),
    within_graph_sd = stats::sd(weights),
    scale_factor = as.numeric(graph_attr(graph, "pcor_scale_factor")),
    stringsAsFactors = FALSE
  )
}

simulate_replication <- function(replication, threshold_limits, threshold_code, cfg) {
  set.seed(cfg$base_seed + as.integer(replication))
  graphs <- prepare_graphs(cfg)
  metrics_for_run <- cfg$metrics

  cascade_rows <- list()
  topology_rows <- list()
  target_rows <- list()
  agreement_rows <- list()
  weight_rows <- list()

  for (graph_key in names(graphs)) {
    graph <- graphs[[graph_key]]
    thresholds <- runif(
      vcount(graph),
      min = threshold_limits[1L],
      max = threshold_limits[2L]
    )

    topology_rows[[graph_key]] <-
      topology_diagnostics(graph, replication, graph_key)
    weight_rows[[graph_key]] <-
      weight_diagnostics(graph, replication, graph_key)

    centralities <- setNames(
      lapply(metrics_for_run, function(metric) compute_centrality(graph, metric)),
      metrics_for_run
    )

    # Every metric begins target tie-breaking from the same random-number state.
    target_rng_state <- .Random.seed
    target_information <- setNames(lapply(metrics_for_run, function(metric) {
      assign(".Random.seed", target_rng_state, envir = .GlobalEnv)
      select_targets(
        graph,
        centralities[[metric]],
        exclude_isolates = cfg$exclude_isolates_as_targets
      )
    }), metrics_for_run)

    for (metric in metrics_for_run) {
      target_info <- target_information[[metric]]
      targets <- target_info$targets

      target_rows[[paste(graph_key, metric, sep = "__")]] <- data.frame(
        replication = replication,
        graph_key = graph_key,
        density_condition = density_from_key(graph_key),
        network = network_from_key(graph_key),
        metric = metric,
        high_target = unname(targets["High"]),
        average_target = unname(targets["Average"]),
        low_target = unname(targets["Low"]),
        high_ties = unname(target_info$tie_counts["High"]),
        average_ties = unname(target_info$tie_counts["Average"]),
        low_ties = unname(target_info$tie_counts["Low"]),
        high_average_same = as.integer(targets["High"] == targets["Average"]),
        high_low_same = as.integer(targets["High"] == targets["Low"]),
        average_low_same = as.integer(targets["Average"] == targets["Low"]),
        centrality_sd = target_info$centrality_sd,
        eligible_n = target_info$eligible_n,
        stringsAsFactors = FALSE
      )

      # Retain the original cascade evaluation order.
      for (target_class in c("High", "Average", "Low")) {
        result <- run_cascade(
          graph,
          unname(targets[target_class]),
          thresholds
        )
        cascade_rows[[paste(graph_key, metric, target_class, sep = "__")]] <-
          data.frame(
            replication = replication,
            density_condition = density_from_key(graph_key),
            network = network_from_key(graph_key),
            metric = metric,
            target = target_class,
            target_vertex = as.integer(unname(targets[target_class])),
            threshold = threshold_code,
            final_active_count = as.integer(result["final_active_count"]),
            final_active = unname(result["final_active"]),
            steps_to_stopping_point = as.integer(
              result["steps_to_stopping_point"]
            ),
            complete_failure = as.integer(
              result["final_active_count"] == cfg$n_nodes
            ),
            stringsAsFactors = FALSE
          )
      }
    }

    metric_pairs <- combn(metrics_for_run, 2L, simplify = FALSE)
    for (pair in metric_pairs) {
      first <- target_information[[pair[1L]]]$targets
      second <- target_information[[pair[2L]]]$targets
      agreement_rows[[paste(graph_key, pair, collapse = "__")]] <- data.frame(
        replication = replication,
        graph_key = graph_key,
        density_condition = density_from_key(graph_key),
        network = network_from_key(graph_key),
        metric_1 = pair[1L],
        metric_2 = pair[2L],
        high_agreement = as.integer(first["High"] == second["High"]),
        average_agreement = as.integer(
          first["Average"] == second["Average"]
        ),
        low_agreement = as.integer(first["Low"] == second["Low"]),
        stringsAsFactors = FALSE
      )
    }
  }

  list(
    cascades = do.call(rbind, cascade_rows),
    topology = do.call(rbind, topology_rows),
    targets = do.call(rbind, target_rows),
    target_agreement = do.call(rbind, agreement_rows),
    weights = do.call(rbind, weight_rows)
  )
}

# =============================================================================
# 4. RUN THE MONTE CARLO STUDY
# =============================================================================

detected_workers <- suppressWarnings(detectCores(logical = TRUE))
if (!is.finite(detected_workers)) detected_workers <- 1L
workers <- max(
  1L,
  min(
    as.integer(N_WORKERS),
    design$n_replications,
    max(1L, detected_workers - 1L)
  )
)

configuration_for_workers <- c(
  design,
  list(
    metrics = metrics,
    graph_keys = graph_keys,
    weight_lower = weight_lower,
    weight_upper = weight_upper,
    minimum_precision_eigenvalue = minimum_precision_eigenvalue,
    exclude_isolates_as_targets = exclude_isolates_as_targets
  )
)

cat(sprintf("Study design: %d\n", STUDY_DESIGN))
cat(sprintf("Threshold condition: %d\n", THRESHOLD_CONDITION))
cat(sprintf("Design: %s\n", design$label))
cat(sprintf("Replications per threshold range: %d\n", design$n_replications))
cat(sprintf("Workers: %d\n", workers))
cat(sprintf("Output directory: %s\n", output_directory))

cluster <- makeCluster(workers)
invisible(clusterEvalQ(cluster, {
  RNGversion("4.4.2")
  RNGkind("Mersenne-Twister", "Inversion", "Rejection")
  library(igraph)
}))
clusterExport(
  cluster,
  c(
    "metrics", "graph_keys", "network_names", "density_names",
    "weight_lower", "weight_upper", "minimum_precision_eigenvalue",
    "exclude_isolates_as_targets", "density_from_key", "network_from_key",
    "adjust_edge_count", "generate_topologies",
    "add_partial_correlation_weights", "prepare_graphs",
    "compute_centrality", "select_targets", "run_cascade",
    "topology_diagnostics", "weight_diagnostics", "simulate_replication"
  ),
  envir = environment()
)

cascade_by_threshold <- vector("list", length(active_threshold_ranges))
topology_raw <- NULL
target_raw <- NULL
target_agreement_raw <- NULL
weight_raw <- NULL

tryCatch(
  {
    for (threshold_index in seq_along(active_threshold_ranges)) {
      threshold_name <- names(active_threshold_ranges)[threshold_index]
      threshold_limits <- active_threshold_ranges[[threshold_index]]
      threshold_code <- threshold_codes[threshold_index]

      cat(sprintf(
        "\nRunning %s thresholds (%s) ...\n",
        threshold_name,
        threshold_code
      ))

      replications <- parLapply(
        cluster,
        seq_len(design$n_replications),
        simulate_replication,
        threshold_limits = threshold_limits,
        threshold_code = threshold_code,
        cfg = configuration_for_workers
      )

      cascade_by_threshold[[threshold_index]] <- do.call(
        rbind,
        lapply(replications, `[[`, "cascades")
      )

      # The same seeded graph and target realizations are reused across threshold
      # ranges, so diagnostics are retained once.
      if (threshold_index == 1L) {
        topology_raw <- do.call(rbind, lapply(replications, `[[`, "topology"))
        target_raw <- do.call(rbind, lapply(replications, `[[`, "targets"))
        target_agreement_raw <- do.call(
          rbind,
          lapply(replications, `[[`, "target_agreement")
        )
        weight_raw <- do.call(rbind, lapply(replications, `[[`, "weights"))
      }

      cat(sprintf(
        "Completed %s: %s cascade results.\n",
        threshold_name,
        format(nrow(cascade_by_threshold[[threshold_index]]), big.mark = ",")
      ))
    }
  },
  finally = {
    try(stopCluster(cluster), silent = TRUE)
  }
)
cluster <- NULL

cascade_raw <- do.call(rbind, cascade_by_threshold)
rownames(cascade_raw) <- NULL
rownames(topology_raw) <- NULL
rownames(target_raw) <- NULL
rownames(target_agreement_raw) <- NULL
rownames(weight_raw) <- NULL

expected_cascade_rows <-
  design$n_replications *
  length(active_threshold_ranges) *
  length(graph_keys) *
  length(metrics) *
  length(target_order)

stopifnot(
  nrow(cascade_raw) == expected_cascade_rows,
  nrow(topology_raw) == design$n_replications * length(graph_keys),
  all(cascade_raw$final_active >= 1 / design$n_nodes),
  all(cascade_raw$final_active <= 1),
  all(cascade_raw$final_active_count == design$n_nodes * cascade_raw$final_active),
  all(cascade_raw$steps_to_stopping_point >= 0L),
  all(cascade_raw$steps_to_stopping_point <= design$n_nodes - 1L),
  all(
    cascade_raw$complete_failure ==
      as.integer(cascade_raw$final_active_count == design$n_nodes)
  )
)

# =============================================================================
# 5. CELL SUMMARIES FOR TABLES AND FIGURES
# =============================================================================

safe_sd <- function(x) {
  if (length(x) >= 2L) stats::sd(x) else NA_real_
}

safe_se <- function(x) {
  if (length(x) >= 2L) stats::sd(x) / sqrt(length(x)) else NA_real_
}

split_by_columns <- function(data, columns) {
  split(
    data,
    interaction(data[columns], drop = TRUE, lex.order = TRUE)
  )
}

cascade_summary <- do.call(rbind, lapply(
  split_by_columns(
    cascade_raw,
    c("density_condition", "network", "metric", "target", "threshold")
  ),
  function(cell) {
    final_active <- cell$final_active
    all_steps <- cell$steps_to_stopping_point
    complete <- cell$final_active_count == design$n_nodes
    complete_steps <- all_steps[complete]
    n_total <- length(final_active)
    n_complete <- length(complete_steps)
    final_se <- safe_se(final_active)
    all_duration_se <- safe_se(all_steps)
    complete_duration_se <- safe_se(complete_steps)

    data.frame(
      density_condition = cell$density_condition[1L],
      network = cell$network[1L],
      metric = cell$metric[1L],
      target = cell$target[1L],
      threshold = cell$threshold[1L],
      mean_final_active = mean(final_active),
      sd_final_active = safe_sd(final_active),
      se_final_active = final_se,
      ci95_lo_final_active = max(0, mean(final_active) - 1.96 * final_se),
      ci95_hi_final_active = min(1, mean(final_active) + 1.96 * final_se),
      mean_steps_unconditional = mean(all_steps),
      sd_steps_unconditional = safe_sd(all_steps),
      se_steps_unconditional = all_duration_se,
      n_complete_failures = n_complete,
      percent_complete_failures = 100 * n_complete / n_total,
      mean_steps_to_complete_failure = if (n_complete) {
        mean(complete_steps)
      } else {
        NA_real_
      },
      sd_steps_to_complete_failure = safe_sd(complete_steps),
      se_steps_to_complete_failure = complete_duration_se,
      ci95_lo_steps_to_complete_failure = if (n_complete >= 2L) {
        max(0, mean(complete_steps) - 1.96 * complete_duration_se)
      } else {
        NA_real_
      },
      ci95_hi_steps_to_complete_failure = if (n_complete >= 2L) {
        mean(complete_steps) + 1.96 * complete_duration_se
      } else {
        NA_real_
      },
      n_replications = n_total,
      stringsAsFactors = FALSE
    )
  }
))
rownames(cascade_summary) <- NULL

high_results <- cascade_raw[cascade_raw$target == "High", ]
low_results <- cascade_raw[cascade_raw$target == "Low", ]
paired_results <- merge(
  high_results[
    , c(
      "replication", "density_condition", "network", "metric", "threshold",
      "final_active"
    )
  ],
  low_results[
    , c(
      "replication", "density_condition", "network", "metric", "threshold",
      "final_active"
    )
  ],
  by = c(
    "replication", "density_condition", "network", "metric", "threshold"
  ),
  suffixes = c("_high", "_low"),
  sort = FALSE
)
paired_results$high_low_difference <-
  paired_results$final_active_high - paired_results$final_active_low

contrast_summary <- do.call(rbind, lapply(
  split_by_columns(
    paired_results,
    c("density_condition", "network", "metric", "threshold")
  ),
  function(cell) {
    difference <- cell$high_low_difference
    difference_se <- safe_se(difference)
    data.frame(
      density_condition = cell$density_condition[1L],
      network = cell$network[1L],
      metric = cell$metric[1L],
      threshold = cell$threshold[1L],
      mean_high_low_difference = mean(difference),
      median_high_low_difference = median(difference),
      sd_high_low_difference = safe_sd(difference),
      se_high_low_difference = difference_se,
      ci95_lo_high_low_difference = mean(difference) - 1.96 * difference_se,
      ci95_hi_high_low_difference = mean(difference) + 1.96 * difference_se,
      percent_high_greater = 100 * mean(difference > 0),
      percent_equal = 100 * mean(difference == 0),
      percent_high_lower = 100 * mean(difference < 0),
      n_replications = length(difference),
      stringsAsFactors = FALSE
    )
  }
))
rownames(contrast_summary) <- NULL

# Table 2 shading: maximum unrounded mean change across graph families within
# each metric-by-density-by-threshold comparison. Exact ties are all flagged.
contrast_summary$is_table2_max <- FALSE
table2_groups <- split_by_columns(
  contrast_summary,
  c("density_condition", "metric", "threshold")
)
for (group in table2_groups) {
  group_key <- with(
    group,
    paste(density_condition, metric, threshold, sep = "__")
  )[1L]
  full_key <- with(
    contrast_summary,
    paste(density_condition, metric, threshold, sep = "__")
  )
  candidates <- which(full_key == group_key)
  contrast_summary$is_table2_max[candidates] <-
    contrast_summary$mean_high_low_difference[candidates] ==
      max(contrast_summary$mean_high_low_difference[candidates])
}

# Derived changes reported in the Results: denser minus sparse Delta_F.
sparse_contrast <- contrast_summary[
  contrast_summary$density_condition == "Sparse",
  c("network", "metric", "threshold", "mean_high_low_difference")
]
names(sparse_contrast)[4L] <- "mean_high_low_difference_sparse"
denser_contrast <- contrast_summary[
  contrast_summary$density_condition == "Denser",
  c("network", "metric", "threshold", "mean_high_low_difference")
]
names(denser_contrast)[4L] <- "mean_high_low_difference_denser"
density_change_summary <- merge(
  sparse_contrast,
  denser_contrast,
  by = c("network", "metric", "threshold"),
  sort = FALSE
)
density_change_summary$denser_minus_sparse_difference <-
  density_change_summary$mean_high_low_difference_denser -
    density_change_summary$mean_high_low_difference_sparse

# Derived most-minus-average and average-minus-least final-cascade gaps.
target_mean_columns <- c(
  "density_condition", "network", "metric", "threshold",
  "target", "mean_final_active"
)
target_means <- cascade_summary[, target_mean_columns]
high_means <- target_means[target_means$target == "High", -5L]
average_means <- target_means[target_means$target == "Average", -5L]
low_means <- target_means[target_means$target == "Low", -5L]
names(high_means)[5L] <- "mean_final_active_most"
names(average_means)[5L] <- "mean_final_active_average"
names(low_means)[5L] <- "mean_final_active_least"
target_position_gaps <- Reduce(
  function(x, y) merge(
    x,
    y,
    by = c("density_condition", "network", "metric", "threshold"),
    sort = FALSE
  ),
  list(high_means, average_means, low_means)
)
target_position_gaps$most_minus_average <-
  target_position_gaps$mean_final_active_most -
    target_position_gaps$mean_final_active_average
target_position_gaps$average_minus_least <-
  target_position_gaps$mean_final_active_average -
    target_position_gaps$mean_final_active_least

topology_summary <- do.call(rbind, lapply(
  split_by_columns(topology_raw, c("density_condition", "network")),
  function(cell) data.frame(
    density_condition = cell$density_condition[1L],
    network = cell$network[1L],
    n_replications = nrow(cell),
    mean_edges = mean(cell$n_edges),
    mean_density = mean(cell$density),
    mean_degree = mean(cell$mean_degree),
    mean_degree_variance = mean(cell$degree_variance),
    sd_degree_variance = safe_sd(cell$degree_variance),
    mean_maximum_degree = mean(cell$maximum_degree),
    mean_global_clustering = mean(cell$global_clustering),
    mean_path_length = mean(cell$mean_path_length),
    percent_connected = 100 * mean(cell$connected),
    mean_components = mean(cell$n_components),
    mean_isolates = mean(cell$n_isolates),
    mean_largest_component = mean(cell$largest_component),
    stringsAsFactors = FALSE
  )
))
rownames(topology_summary) <- NULL

target_summary <- do.call(rbind, lapply(
  split_by_columns(target_raw, c("density_condition", "network", "metric")),
  function(cell) data.frame(
    density_condition = cell$density_condition[1L],
    network = cell$network[1L],
    metric = cell$metric[1L],
    n_replications = nrow(cell),
    mean_high_ties = mean(cell$high_ties),
    mean_average_ties = mean(cell$average_ties),
    mean_low_ties = mean(cell$low_ties),
    percent_high_tied = 100 * mean(cell$high_ties > 1L),
    percent_average_tied = 100 * mean(cell$average_ties > 1L),
    percent_low_tied = 100 * mean(cell$low_ties > 1L),
    percent_high_average_same = 100 * mean(cell$high_average_same),
    percent_high_low_same = 100 * mean(cell$high_low_same),
    percent_average_low_same = 100 * mean(cell$average_low_same),
    mean_centrality_sd = mean(cell$centrality_sd, na.rm = TRUE),
    stringsAsFactors = FALSE
  )
))
rownames(target_summary) <- NULL

target_agreement_summary <- do.call(rbind, lapply(
  split_by_columns(
    target_agreement_raw,
    c("density_condition", "network", "metric_1", "metric_2")
  ),
  function(cell) data.frame(
    density_condition = cell$density_condition[1L],
    network = cell$network[1L],
    metric_1 = cell$metric_1[1L],
    metric_2 = cell$metric_2[1L],
    n_replications = nrow(cell),
    percent_high_agreement = 100 * mean(cell$high_agreement),
    percent_average_agreement = 100 * mean(cell$average_agreement),
    percent_low_agreement = 100 * mean(cell$low_agreement),
    stringsAsFactors = FALSE
  )
))
rownames(target_agreement_summary) <- NULL

weight_summary <- do.call(rbind, lapply(
  split_by_columns(weight_raw, c("density_condition", "network")),
  function(cell) {
    total_n <- sum(cell$n_weights)
    total_sum <- sum(cell$weight_sum)
    total_sum_squares <- sum(cell$weight_sum_squares)
    pooled_mean <- total_sum / total_n
    pooled_variance <-
      (total_sum_squares - total_sum^2 / total_n) / (total_n - 1L)

    data.frame(
      density_condition = cell$density_condition[1L],
      network = cell$network[1L],
      n_weights = total_n,
      mean_weight = pooled_mean,
      sd_weight = sqrt(max(0, pooled_variance)),
      minimum_weight = min(cell$minimum_weight),
      maximum_weight = max(cell$maximum_weight),
      mean_within_graph_sd = mean(cell$within_graph_sd),
      percent_scaled = 100 * mean(cell$scale_factor < 1),
      mean_scale_factor = mean(cell$scale_factor),
      minimum_scale_factor = min(cell$scale_factor),
      stringsAsFactors = FALSE
    )
  }
))
rownames(weight_summary) <- NULL

# Table 3 intentionally contains only the least and most central targets.
table3_complete_failure_duration <- cascade_summary[
  cascade_summary$target %in% c("Low", "High"),
  c(
    "density_condition", "network", "metric", "target", "threshold",
    "mean_steps_to_complete_failure", "sd_steps_to_complete_failure",
    "n_complete_failures", "percent_complete_failures", "n_replications"
  )
]

names(table3_complete_failure_duration)[
  names(table3_complete_failure_duration) == "target"
] <- "target_centrality"

# Table 3 shading: minimum unrounded conditional mean across graph families.
# Cells without a complete cascade are undefined and cannot be flagged.
table3_complete_failure_duration$is_table3_fastest <- FALSE
table3_groups <- split_by_columns(
  table3_complete_failure_duration,
  c("density_condition", "metric", "target_centrality", "threshold")
)
for (group in table3_groups) {
  group_key <- with(
    group,
    paste(density_condition, metric, target_centrality, threshold, sep = "__")
  )[1L]
  full_key <- with(
    table3_complete_failure_duration,
    paste(density_condition, metric, target_centrality, threshold, sep = "__")
  )
  candidates <- which(full_key == group_key)
  estimable <- candidates[
    !is.na(
      table3_complete_failure_duration$mean_steps_to_complete_failure[candidates]
    )
  ]
  if (length(estimable)) {
    fastest <- min(
      table3_complete_failure_duration$mean_steps_to_complete_failure[estimable]
    )
    table3_complete_failure_duration$is_table3_fastest[estimable] <-
      table3_complete_failure_duration$mean_steps_to_complete_failure[estimable] ==
        fastest
  }
}

# Density comparison for the conditional Table 3 outcome. These means can be
# based on different completed-cascade subsets, so counts and percentages are
# retained for both densities.
table3_density_fields <- c(
  "network", "metric", "target_centrality", "threshold",
  "mean_steps_to_complete_failure", "n_complete_failures",
  "percent_complete_failures"
)
table3_sparse <- table3_complete_failure_duration[
  table3_complete_failure_duration$density_condition == "Sparse",
  table3_density_fields
]
table3_denser <- table3_complete_failure_duration[
  table3_complete_failure_duration$density_condition == "Denser",
  table3_density_fields
]
names(table3_sparse)[5:7] <- paste0(names(table3_sparse)[5:7], "_sparse")
names(table3_denser)[5:7] <- paste0(names(table3_denser)[5:7], "_denser")
conditional_duration_density_comparison <- merge(
  table3_sparse,
  table3_denser,
  by = c("network", "metric", "target_centrality", "threshold"),
  sort = FALSE
)
conditional_duration_density_comparison$sparse_minus_denser_steps <-
  conditional_duration_density_comparison$mean_steps_to_complete_failure_sparse -
    conditional_duration_density_comparison$mean_steps_to_complete_failure_denser

# =============================================================================
# 6. WRITE AUDITABLE DATA FILES
# =============================================================================

write.csv(
  cascade_raw,
  file.path(output_directory, "cascade_replication_level.csv"),
  row.names = FALSE
)

write.csv(
  cascade_summary,
  file.path(output_directory, "cascade_summary_all_conditions.csv"),
  row.names = FALSE
)
write.csv(
  paired_results,
  file.path(output_directory, "high_low_difference_replication_level.csv"),
  row.names = FALSE
)
write.csv(
  contrast_summary,
  file.path(output_directory, "table2_change_in_failure_rate.csv"),
  row.names = FALSE
)
write.csv(
  table3_complete_failure_duration,
  file.path(output_directory, "table3_steps_to_complete_failure.csv"),
  row.names = FALSE
)
write.csv(
  density_change_summary,
  file.path(output_directory, "density_change_in_high_low_difference.csv"),
  row.names = FALSE
)
write.csv(
  target_position_gaps,
  file.path(output_directory, "target_position_failure_rate_gaps.csv"),
  row.names = FALSE
)
write.csv(
  conditional_duration_density_comparison,
  file.path(output_directory, "conditional_duration_density_comparison.csv"),
  row.names = FALSE
)
write.csv(
  topology_raw,
  file.path(output_directory, "topology_diagnostics_replication_level.csv"),
  row.names = FALSE
)
write.csv(
  topology_summary,
  file.path(output_directory, "topology_diagnostics_summary.csv"),
  row.names = FALSE
)
write.csv(
  target_raw,
  file.path(output_directory, "target_selection_replication_level.csv"),
  row.names = FALSE
)
write.csv(
  target_summary,
  file.path(output_directory, "target_selection_summary.csv"),
  row.names = FALSE
)
write.csv(
  target_agreement_raw,
  file.path(output_directory, "target_agreement_replication_level.csv"),
  row.names = FALSE
)
write.csv(
  target_agreement_summary,
  file.path(output_directory, "target_agreement_summary.csv"),
  row.names = FALSE
)
write.csv(
  weight_summary,
  file.path(output_directory, "edge_weight_diagnostics.csv"),
  row.names = FALSE
)
write.csv(
  weight_raw,
  file.path(output_directory, "edge_weight_diagnostics_replication_level.csv"),
  row.names = FALSE
)

saveRDS(
  list(
    configuration = configuration_for_workers,
    cascade_summary = cascade_summary,
    table2 = contrast_summary,
    table3 = table3_complete_failure_duration,
    density_change_summary = density_change_summary,
    target_position_gaps = target_position_gaps,
    conditional_duration_density_comparison =
      conditional_duration_density_comparison,
    topology_summary = topology_summary,
    target_summary = target_summary,
    target_agreement_summary = target_agreement_summary,
    edge_weight_summary = weight_summary
  ),
  file.path(output_directory, "manuscript_simulation_results.rds")
)

# =============================================================================
# 7. MANUSCRIPT FIGURES 5-9
# =============================================================================

if (MAKE_FIGURES) {
  plot_summary <- cascade_summary
  plot_summary$target <- factor(
    plot_summary$target,
    levels = target_order
  )
  plot_summary$network <- factor(
    plot_summary$network,
    levels = c("Erdos-Renyi", "Barabasi-Albert", "Watts-Strogatz")
  )
  plot_summary$density_label <- factor(
    plot_summary$density_condition,
    levels = c("Sparse", "Denser"),
    labels = c(
      "Sparse (20 edges; ~10% density)",
      "Denser (40 edges; ~20% density)"
    )
  )
  plot_summary$threshold_label <- factor(
    plot_summary$threshold,
    levels = threshold_codes,
    labels = unname(threshold_display[threshold_codes])
  )

  metric_titles <- c(
    degree = "degree",
    strength = "strength",
    closeness = "harmonic closeness",
    betweenness = "betweenness"
  )
  metric_figure_numbers <- c(
    degree = "05",
    strength = "06",
    closeness = "07",
    betweenness = "08"
  )

  for (metric in metrics) {
    metric_data <- droplevels(plot_summary[plot_summary$metric == metric, ])
    plot_metric <- ggplot(
      metric_data,
      aes(
        target,
        mean_final_active,
        shape = network,
        color = network,
        group = network
      )
    ) +
      geom_errorbar(
        aes(
          ymin = ci95_lo_final_active,
          ymax = ci95_hi_final_active
        ),
        width = 0.15,
        linewidth = 0.4,
        position = position_dodge(width = 0.4)
      ) +
      geom_line(
        linewidth = 0.6,
        position = position_dodge(width = 0.4)
      ) +
      geom_point(
        size = 3,
        position = position_dodge(width = 0.4)
      ) +
      scale_y_continuous(
        "Mean final active proportion",
        limits = c(0, 1)
      ) +
      scale_x_discrete(NULL) +
      scale_color_brewer("Graph family", palette = "Set1") +
      scale_shape_manual(
        "Graph family",
        values = c(
          "Erdos-Renyi" = 16,
          "Barabasi-Albert" = 17,
          "Watts-Strogatz" = 15
        )
      ) +
      facet_grid(density_label ~ threshold_label) +
      labs(title = paste("Cascade size by", metric_titles[[metric]], "target selection")) +
      theme_minimal(base_size = 13) +
      theme(
        strip.text = element_text(face = "bold"),
        plot.title = element_text(face = "bold", hjust = 0.5),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "bottom"
      )

    metric_file <- sprintf(
      "figure%s_failure_rate_%s.png",
      metric_figure_numbers[[metric]],
      metric
    )
    ggsave(
      file.path(output_directory, metric_file),
      plot_metric,
      width = 14,
      height = 10,
      dpi = 300,
      bg = "white"
    )
  }

  contrast_plot_data <- contrast_summary
  contrast_plot_data$threshold_label <- factor(
    contrast_plot_data$threshold,
    levels = threshold_codes,
    labels = unname(threshold_display[threshold_codes])
  )
  contrast_plot_data$density_label <- factor(
    contrast_plot_data$density_condition,
    levels = c("Sparse", "Denser"),
    labels = c(
      "Sparse (20 edges; ~10% density)",
      "Denser (40 edges; ~20% density)"
    )
  )
  contrast_plot_data$metric_label <- factor(
    contrast_plot_data$metric,
    levels = metrics,
    labels = c("Degree", "Strength", "Harmonic closeness", "Betweenness")
  )
  contrast_plot_data$network <- factor(
    contrast_plot_data$network,
    levels = c("Erdos-Renyi", "Barabasi-Albert", "Watts-Strogatz")
  )

  plot_difference <- ggplot(
    contrast_plot_data,
    aes(
      threshold_label,
      mean_high_low_difference,
      color = network,
      group = network
    )
  ) +
    geom_hline(yintercept = 0, color = "grey65", linewidth = 0.35) +
    geom_errorbar(
      aes(
        ymin = ci95_lo_high_low_difference,
        ymax = ci95_hi_high_low_difference
      ),
      width = 0.10,
      linewidth = 0.4,
      position = position_dodge(width = 0.25)
    ) +
    geom_line(
      linewidth = 0.7,
      position = position_dodge(width = 0.25)
    ) +
    geom_point(
      size = 2.4,
      position = position_dodge(width = 0.25)
    ) +
    scale_color_brewer("Graph family", palette = "Set1") +
    scale_y_continuous(
      expression("Mean change in failure rate (" * Delta[F] * ")")
    ) +
    scale_x_discrete(NULL) +
    facet_grid(density_label ~ metric_label) +
    labs(
      title = expression(
        "Change in failure rate (" * Delta[F] * ") by density"
      )
    ) +
    theme_minimal(base_size = 12) +
    theme(
      strip.text = element_text(face = "bold"),
      plot.title = element_text(face = "bold", hjust = 0.5),
      panel.grid.minor = element_blank(),
      axis.text.x = element_text(angle = 30, hjust = 1),
      legend.position = "bottom"
    )

  ggsave(
    file.path(output_directory, "figure09_change_in_failure_rate.png"),
    plot_difference,
    width = 15,
    height = 8,
    dpi = 300,
    bg = "white"
  )
}

# =============================================================================
# 8. REPRESENTATIVE TOPOLOGY FIGURE
# =============================================================================

if (MAKE_REPRESENTATIVE_TOPOLOGY_FIGURE) {
  representative_features <- c(
    "degree_variance",
    "maximum_degree",
    "global_clustering",
    "mean_path_length",
    "n_components",
    "n_isolates",
    "largest_component"
  )

  choose_representative <- function(cell) {
    modal_connected <- as.integer(mean(cell$connected) >= 0.5)
    candidates <- cell[cell$connected == modal_connected, , drop = FALSE]
    values <- candidates[, representative_features, drop = FALSE]
    centers <- vapply(values, median, numeric(1L), na.rm = TRUE)
    scales <- vapply(values, stats::mad, numeric(1L), na.rm = TRUE)
    fallback <- vapply(values, stats::sd, numeric(1L), na.rm = TRUE)
    invalid_scale <- !is.finite(scales) | scales <= 0
    scales[invalid_scale] <- fallback[invalid_scale]
    scales[!is.finite(scales) | scales <= 0] <- 1
    standardized <- sweep(as.matrix(values), 2L, centers, "-")
    standardized <- sweep(standardized, 2L, scales, "/")
    candidates$representativeness_score <- rowSums(standardized^2)
    candidates$modal_connected <- modal_connected
    candidates <- candidates[order(
      candidates$representativeness_score,
      candidates$replication
    ), , drop = FALSE]
    candidates[1L, , drop = FALSE]
  }

  selected_representatives <- do.call(rbind, lapply(
    split_by_columns(topology_raw, c("density_condition", "network")),
    choose_representative
  ))
  rownames(selected_representatives) <- NULL

  graph_diagnostics_for_check <- function(graph) {
    graph_degree <- degree(graph)
    graph_components <- components(graph)
    c(
      n_edges = gsize(graph),
      degree_variance = stats::var(graph_degree),
      maximum_degree = max(graph_degree),
      global_clustering = transitivity(
        graph,
        type = "global",
        isolates = "zero"
      ),
      mean_path_length = mean_distance(
        graph,
        directed = FALSE,
        unconnected = TRUE,
        weights = NA
      ),
      n_components = graph_components$no,
      n_isolates = sum(graph_degree == 0),
      largest_component = max(graph_components$csize),
      connected = as.integer(graph_components$no == 1L)
    )
  }

  regenerate_topologies <- function(replication) {
    set.seed(design$base_seed + as.integer(replication))
    list(
      Sparse = generate_topologies(
        design$n_nodes,
        design$edge_counts[["Sparse"]],
        design$ws_rewiring
      ),
      Denser = generate_topologies(
        design$n_nodes,
        design$edge_counts[["Denser"]],
        design$ws_rewiring
      )
    )
  }

  network_codes <- c(
    "Erdos-Renyi" = "ER",
    "Barabasi-Albert" = "BA",
    "Watts-Strogatz" = "WS"
  )
  representative_graphs <- list()

  for (row_index in seq_len(nrow(selected_representatives))) {
    selected_row <- selected_representatives[row_index, , drop = FALSE]
    regenerated <- regenerate_topologies(selected_row$replication)
    network_code <- network_codes[[selected_row$network]]
    graph <- regenerated[[selected_row$density_condition]][[network_code]]
    observed <- graph_diagnostics_for_check(graph)
    expected <- unlist(selected_row[names(observed)], use.names = TRUE)
    if (!isTRUE(all.equal(
      as.numeric(observed),
      as.numeric(expected),
      tolerance = 1e-12
    ))) {
      stop("A representative graph did not match its stored diagnostics.")
    }
    representative_graphs[[paste(
      selected_row$density_condition,
      network_code,
      sep = "_"
    )]] <- graph
  }

  write.csv(
    selected_representatives,
    file.path(output_directory, "representative_topology_selection.csv"),
    row.names = FALSE
  )

  topology_colors <- c(
    "Barabasi-Albert" = "#D95F59",
    "Watts-Strogatz" = "#4C9F70",
    "Erdos-Renyi" = "#4C78A8"
  )

  graph_layout <- function(graph, network, replication) {
    if (network == "Watts-Strogatz") {
      return(layout_in_circle(graph, order = seq_len(vcount(graph))))
    }
    set.seed(900000L + as.integer(replication))
    layout_with_fr(graph, niter = 1500L, grid = "nogrid")
  }

  plot_representative_graph <- function(graph, selected_row) {
    network <- selected_row$network
    density <- selected_row$density_condition
    coordinates <- graph_layout(
      graph,
      network,
      selected_row$replication
    )
    coordinates <- norm_coords(
      coordinates,
      xmin = -1,
      xmax = 1,
      ymin = -1,
      ymax = 1
    )
    graph_degree <- degree(graph)

    if (network == "Watts-Strogatz") {
      edge_list <- as_edgelist(graph, names = FALSE)
      ordinary_distance <- abs(edge_list[, 1L] - edge_list[, 2L])
      circular_distance <- pmin(
        ordinary_distance,
        vcount(graph) - ordinary_distance
      )
      lattice_radius <- if (density == "Sparse") 1L else 2L
      long_range <- circular_distance > lattice_radius
      edge_color <- ifelse(long_range, "#E6862A", "#A7AFB8A6")
      edge_width <- ifelse(long_range, 2.25, 1.15)
    } else {
      edge_color <- rep("#9DA7B3A6", gsize(graph))
      edge_width <- rep(1.15, gsize(graph))
    }

    plot(
      graph,
      layout = coordinates,
      vertex.size = 9 + 2 * sqrt(graph_degree),
      vertex.color = adjustcolor(topology_colors[[network]], alpha.f = 0.92),
      vertex.frame.color = "white",
      vertex.label = NA,
      edge.color = edge_color,
      edge.width = edge_width,
      edge.curved = 0,
      margin = 0.17,
      rescale = FALSE,
      asp = 1
    )
    title(main = network, cex.main = 1.15, font.main = 2, col.main = "#17324D")
    mtext(
      sprintf(
        "rep %d | m=%d | C=%.3f | L=%.3f",
        selected_row$replication,
        selected_row$n_edges,
        selected_row$global_clustering,
        selected_row$mean_path_length
      ),
      side = 1,
      line = 0.35,
      cex = 0.70,
      col = "#5F6B76"
    )
  }

  draw_topology_panel <- function() {
    network_order <- c("Barabasi-Albert", "Watts-Strogatz", "Erdos-Renyi")
    density_order <- c("Sparse", "Denser")
    layout(matrix(seq_len(6L), nrow = 2L, byrow = TRUE))
    par(
      mar = c(2.5, 2.2, 3.0, 1.2),
      oma = c(4.2, 4.4, 3.3, 1.0),
      xpd = NA,
      bg = "white",
      family = "sans"
    )

    for (density in density_order) {
      for (network in network_order) {
        selected_row <- selected_representatives[
          selected_representatives$density_condition == density &
            selected_representatives$network == network,
          , drop = FALSE
        ]
        graph_key <- paste(density, network_codes[[network]], sep = "_")
        plot_representative_graph(
          representative_graphs[[graph_key]],
          selected_row
        )
        if (network == network_order[1L]) {
          density_label <- if (density == "Sparse") {
            "20 edges\ndensity 0.105"
          } else {
            "40 edges\ndensity 0.211"
          }
          mtext(
            density_label,
            side = 2,
            line = 1.2,
            cex = 0.87,
            font = 2,
            col = "#17324D"
          )
        }
      }
    }

    mtext(
      "Representative simulated graph topologies",
      side = 3,
      outer = TRUE,
      line = 1.45,
      cex = 1.45,
      font = 2,
      col = "#17324D"
    )
    mtext(
      paste(
        "Node size is proportional to degree. Each graph is nearest its cell",
        "median across recorded topology diagnostics."
      ),
      side = 1,
      outer = TRUE,
      line = 2.0,
      cex = 0.82,
      col = "#5F6B76"
    )
  }

  png(
    file.path(output_directory, "figure04_representative_topologies.png"),
    width = 3900,
    height = 2600,
    res = 320,
    bg = "white",
    type = "cairo-png"
  )
  draw_topology_panel()
  dev.off()
}

# =============================================================================
# 9. REPRODUCIBILITY RECORD
# =============================================================================

configuration_lines <- c(
  sprintf("Study design: %d", STUDY_DESIGN),
  sprintf("Threshold condition: %d", THRESHOLD_CONDITION),
  sprintf("Study label: %s", design$label),
  sprintf("Completed: %s", format(Sys.time(), tz = "America/Los_Angeles")),
  sprintf("Base seed: %d", design$base_seed),
  sprintf("Replications per threshold range: %d", design$n_replications),
  sprintf("Workers used: %d", workers),
  sprintf("Nodes: %d", design$n_nodes),
  sprintf(
    "Edges: Sparse=%d; Denser=%d",
    design$edge_counts[["Sparse"]],
    design$edge_counts[["Denser"]]
  ),
  sprintf("Watts-Strogatz rewiring probability: %.2f", design$ws_rewiring),
  sprintf("Metrics: %s", paste(metrics, collapse = ", ")),
  sprintf("Threshold ranges: %s", paste(threshold_codes, collapse = ", ")),
  sprintf("Validated R version: %s", expected_versions[["R"]]),
  sprintf("Validated igraph version: %s", expected_versions[["igraph"]]),
  sprintf("RNG kind: %s", paste(RNGkind(), collapse = ", ")),
  sprintf("Review script MD5: %s", unname(tools::md5sum(script_path))),
  paste0(
    "Primary-study archived backbone SHA256: ",
    "E9CC833FB8D35F399BBA578C9B33D812EDE6AC108D29DF9EA544C927EDF59D74"
  ),
  paste0(
    "Primary-study archived beta=.10 runner SHA256: ",
    "9E5EA48D317CE8F801BF602A7EA11CF4C4C8E72F32EE6E7E1D755FA9E7BED52D"
  ),
  paste0(
    "Primary-study archived result RDS SHA256: ",
    "22E42EB105DDB4709793D854F4C085D6A50E0770570B5604E0A98CCD536B4EF6"
  ),
  "",
  capture.output(sessionInfo())
)
writeLines(
  configuration_lines,
  file.path(output_directory, "reproducibility_record.txt")
)

file.copy(
  script_path,
  file.path(output_directory, basename(script_path)),
  overwrite = TRUE
)

generated_output_names <- c(
  "cascade_replication_level.csv",
  "cascade_summary_all_conditions.csv",
  "high_low_difference_replication_level.csv",
  "table2_change_in_failure_rate.csv",
  "table3_steps_to_complete_failure.csv",
  "density_change_in_high_low_difference.csv",
  "target_position_failure_rate_gaps.csv",
  "conditional_duration_density_comparison.csv",
  "topology_diagnostics_replication_level.csv",
  "topology_diagnostics_summary.csv",
  "target_selection_replication_level.csv",
  "target_selection_summary.csv",
  "target_agreement_replication_level.csv",
  "target_agreement_summary.csv",
  "edge_weight_diagnostics.csv",
  "edge_weight_diagnostics_replication_level.csv",
  "manuscript_simulation_results.rds",
  "reproducibility_record.txt",
  basename(script_path)
)
if (MAKE_FIGURES) {
  generated_output_names <- c(
    generated_output_names,
    sprintf("figure%02d_failure_rate_%s.png", 5:8, metrics),
    "figure09_change_in_failure_rate.png"
  )
}
if (MAKE_REPRESENTATIVE_TOPOLOGY_FIGURE) {
  generated_output_names <- c(
    generated_output_names,
    "representative_topology_selection.csv",
    "figure04_representative_topologies.png"
  )
}
output_files <- file.path(output_directory, generated_output_names)
if (any(!file.exists(output_files))) {
  stop(
    "Expected output files were not generated: ",
    paste(basename(output_files[!file.exists(output_files)]), collapse = ", ")
  )
}
output_hashes <- data.frame(
  file = basename(output_files),
  md5 = unname(tools::md5sum(output_files)),
  stringsAsFactors = FALSE
)
write.csv(
  output_hashes,
  file.path(output_directory, "output_file_md5.csv"),
  row.names = FALSE
)

cat("\nSimulation and all requested summaries completed successfully.\n")
cat(sprintf("Primary output: %s\n", normalizePath(output_directory)))