airtool.qmd

---
title: "AIR -- This is prototype software, use without warranty."
format: 
  dashboard:
    theme: spacelab
    orientation: columns
    # embed-resources: true
    header-includes:
      - '<meta name="color-scheme" content="light">'

resources:
  - readme_md_files
server: shiny
---

```{=html}
<style>
  .plot-container, .shiny-text-output, .panel, .card {
    border: none;
    box-shadow: none;
  }
  .custom-text {
    font-size: 20px;
    color: green;
    font-weight: bold;
  }

  /* FILE-INPUT */
  
  /* Browse button base look  (unchanged colour) */
  .shiny-input-container .btn-file{
    background:#1976d2 !important;
    color:#fff        !important;
    border:0          !important;
    border-radius:4px !important;
    box-shadow:none   !important;
    padding:8px 16px;
    text-transform:uppercase;
    transition:box-shadow .2s cubic-bezier(.4,0,.2,1);
  }
  
  /* Hover / focus — add the same elevation you see on other buttons */
  .shiny-input-container .btn-file:hover,
  .shiny-input-container .btn-file:focus{
    background:#1976d2 !important;           /* stay blue            */
    box-shadow:0 4px 6px rgba(0,0,0,.25);     /* subtle, matches rest */
    outline:0;
  }
  
  /* Active click (tiny dip) */
  .shiny-input-container .btn-file:active{
    background:#1565c0 !important;
    box-shadow:0 2px 4px rgba(0,0,0,.25);
  }
  
  /* Progress bar wrapper — push down & fill width */
  .shiny-input-container .shiny-file-input-progress{
    display:block;
    width:100%;
    margin:10px 0 0 0;                       /* 10 px top margin fixes overlap */
  }
  
  /* Progress bar track & fill */
  .shiny-input-container .progress{
    height:8px;                               /* less “squished”     */
    margin:0;
    border-radius:4px;
  }
  .shiny-input-container .progress-bar{
    background:#1976d2;
    border-radius:4px;
  }

  /* progress bar separation */
  /* The bar itself (not the wrapper) needs the margin */
  .progress.shiny-file-input-progress {
    margin-top: 10px;            /* pushes bar below Browse button */
  }
  
  /* let the shadow show through */
  .shiny-input-container .btn-file {
    overflow: visible !important;  /* un-clip shadow                     */
  }
  
  /* optional: trim the corner radius a bit more */
  .shiny-input-container .btn-file,
  .progress.shiny-file-input-progress {
    border-radius: 2px !important; /* crisper corners                    */
  }
  
  /* for Browse buttons & progress bar */
  .shiny-input-container .btn-file          { box-shadow:none; }            /* reset */
  .shiny-input-container .btn-file:hover,
  .shiny-input-container .btn-file:focus    { box-shadow:0 4px 8px rgba(0,0,0,.30)!important; }
  
  .progress.shiny-file-input-progress,
  .progress.shiny-file-input-progress .progress-bar{
    height:12px;
    border-radius:4px;
    background:#1976d2;
    background-image:none!important;        /* kill default stripes/dots */
  }
  
  /* optional tighter corners */
  .shiny-input-container .btn-file          { border-radius:2px; }

</style>

<style> .main-container { max-width: unset; } </style>
<div id="style-test" style="background-color: white; color: black; position: absolute; left: -9999px;"></div>

<!-- Hidden element for detecting overridden styles -->
<div id="style-test" style="background-color: white; color: black; display: none;"></div>

<script>
document.addEventListener("DOMContentLoaded", () => {
  const testElem = document.getElementById("style-test");
  if (!testElem) return;

  const computedBg = window.getComputedStyle(testElem).backgroundColor;
  const computedColor = window.getComputedStyle(testElem).color;

  if (computedBg !== "rgb(255, 255, 255)" || computedColor !== "rgb(0, 0, 0)") {
    alert("It appears a dark mode extension might be overriding the styles. Please disable it for the best experience.");
  }
});
</script>
```

```{r initial-setup}
#| context: setup
#| echo: false
#| include: false

# Choose conservative defaults for general users
DEFAULT_XMS_GB <- 1   # 1g initial heap
DEFAULT_XMX_GB <- 4   # 4g max heap

# Hard safety ceilings so some lunatic doesn't ask for 5000g on a ThinkPad
MAX_ALLOWED_GB <- 256  # tune this to whatever you consider sane

read_heap_gb <- function(var, default_val) {
  raw <- Sys.getenv(var, unset = "")
  if (!nzchar(raw)) {
    return(default_val)
  }

  # numeric conversion with basic sanity checks
  val <- suppressWarnings(as.numeric(raw))

  if (is.na(val)) {
    return(default_val)
  }

  # must be positive
  if (val <= 0) {
    return(default_val)
  }

  # cap outrageous requests
  val <- min(val, MAX_ALLOWED_GB)

  # floor to integer GB (JVM wants clean "16g", not "16.2g")
  floor(val)
}

calc_heap_settings <- function() {
  xms_gb <- read_heap_gb("AIR_JAVA_XMS_GB", DEFAULT_XMS_GB)
  xmx_gb <- read_heap_gb("AIR_JAVA_XMX_GB", DEFAULT_XMX_GB)

  # ensure logical ordering: Xms <= Xmx
  if (xms_gb > xmx_gb) {
    # if user fat-fingered it, bump max up to match min
    xmx_gb <- xms_gb
  }

  list(
    xms_gb = xms_gb,
    xmx_gb = xmx_gb
  )
}

heap <- calc_heap_settings()

params <- c(
  paste0("-Xms", heap$xms_gb, "g"),
  paste0("-Xmx", heap$xmx_gb, "g")
)

options(shiny.maxRequestSize = Inf)

suppressMessages({
  library(AIPW)
  library(base64enc)
  library(DiagrammeR)
  library(DiagrammeRsvg)  # SVG export
  library(doParallel)
  library(dplyr)
  library(e1071)
  library(earth)
  library(foreach)
  library(ggplot2)
  library(hal9001)
  library(hash)
  library(here)
  library(jsonlite)
  library(logger)
  library(nnet)
  library(rJava)
  library(randomForest)
  library(readr)
  library(rpart)
  library(rsvg)
  library(scales)
  library(sets)
  library(shiny)
  library(shinyjs)
  library(shinyWidgets)
  library(sl3)
  library(tidyr)
  library(tmle3)
  library(xgboost)
  library(xml2)
})

## ── session-wide defaults ───────────────────────────────────────────────
AIRHome <- here::here()
setwd(AIRHome)

log_dir  <- file.path(AIRHome, "logs")
dir.create(log_dir, recursive = TRUE, showWarnings = FALSE)


# ── guard: run heavy log setup only once per R-session ───────────────────
if (exists(".airtool_log_init", envir = .GlobalEnv)) {
  log_file <- get(".airtool_log_file", envir = .GlobalEnv)
} else {
  log_file <- file.path(
    log_dir,
    strftime(Sys.time(), "%Y-%m-%d_%H-%M-%S"),
    "airtool.log"
  )
  dir.create(dirname(log_file), recursive = TRUE, showWarnings = FALSE)
  assign(".airtool_log_init", TRUE,      .GlobalEnv)
  assign(".airtool_log_file", log_file,  .GlobalEnv)
}

## ── 1. logger setup  ----------------------------------------------------
logger::log_threshold("DEBUG")
layout_plain <- function(level, msg, namespace, .logcall, .topcall, .topenv) {
  sprintf("%s [%s] %s", level, format(Sys.time(), "%Y-%m-%d %H:%M:%S"), as.character(msg))
}
logger::log_layout(layout_plain)
logger::log_appender(logger::appender_file(log_file))

for (nm in c("debug", "info", "warn", "error")) {
  f <- get(paste0("log_", nm), envir = asNamespace("logger"))
  assign(paste0("log_", nm), f, .GlobalEnv)
}
log_info("AIR Tool startup – log initialised.")

## ── 2. Logging connections: anchor in globalenv! ------------------------
if (!exists(".airtool_log_con", envir = .GlobalEnv)) {
  log_con <- file(log_file, open = "a")
  assign(".airtool_log_con", log_con, .GlobalEnv)
} else {
  log_con <- get(".airtool_log_con", envir = .GlobalEnv)
}

## Quarto/knitr needs stdout for a while; only divert messages
sink(log_con, append = TRUE, type = "message")

## ── 3. Clean-up: only close what we opened, and ONLY at process end -----
.on_session_end <- function(...) {
  # Remove all sinks attached to this connection
  if (sink.number(type = "message") > 0)
    sink(type = "message")
  if (exists(".airtool_log_con", envir = .GlobalEnv)) {
    con <- get(".airtool_log_con", .GlobalEnv)
    if (isOpen(con)) close(con)
    rm(".airtool_log_con", envir = .GlobalEnv)
  }
}

reg.finalizer(.GlobalEnv, .on_session_end, onexit = TRUE)
if ("shiny" %in% loadedNamespaces()) {
  shiny::onStop(.on_session_end)
} else {
  setHook(packageEvent("shiny", "onLoad"),
          function(...) shiny::onStop(.on_session_end),
          action = "append")
}

## --- JVM streams, stack trace, log pruning as before ---
redirect_java_streams <- function(target) {
  fos <- rJava::.jnew("java/io/FileOutputStream", target, TRUE)
  os  <- rJava::.jcast(fos, "java/io/OutputStream")
  ps  <- rJava::.jnew("java/io/PrintStream", os)
  rJava::.jcall("java/lang/System", "V", "setOut", ps)
  rJava::.jcall("java/lang/System", "V", "setErr", ps)
}

options(
  shiny.fullstacktrace = TRUE,
  error = function() {
    logger::log_error(paste(capture.output(traceback(20L)), collapse = "\n"))
    invokeRestart("abort")
  }
)

time_threshold <- Sys.time() - 30 * 24 * 60 * 60
old_logs <- list.files(log_dir, full.names = TRUE, recursive = TRUE)
old_logs <- old_logs[file.info(old_logs)$mtime < time_threshold]
if (length(old_logs)) file.remove(old_logs)

theme_set(ggplot2::theme_gray(base_family = "DejaVu Sans"))  # ggplot, grid, etc.

#options(warn = 2)
#options(shiny.trace = TRUE)
options(
  shiny.fullstacktrace = TRUE,
  shiny.error = function() {
    cat("\n=== SHINY ERROR =========================================\n")
    print(sys.calls())
    cat("---------------------------------------------------------\n")
    # re-throw so Shiny still shows the UI banner
    geterrmessage()
  }
)

set.seed(123)

## Also tell the R graphics device so axis labels use it too
grDevices::png(type = "cairo")            # ensure Pango/Cairo backend
par(family = "DejaVu Sans")               # default plotting family

# Constants for Java JDK URLs
TETRAD_PATH <- Sys.getenv("TETRAD_PATH")

if (TETRAD_PATH == "") {
  stop("The TETRAD_PATH environment variable is not set. Please set it to the path of the Tetrad jar.")
}

# ––– shield any .jcall from killing R –––––––––––––––––––––––––––––––––
safe_jcall <- function(...) {
  tryCatch(
    .jcall(...),
    Throwable = function(e) {
      msg <- .jcall("java/lang/Throwable", "S", "toString", e)
      log_debug("Java exception trapped: {msg}", msg = msg)
      NULL                                   # propagate safe NULL
    }
  )
}

# make it visible to Reference-class methods & workers
assign("safe_jcall", safe_jcall, envir = .GlobalEnv)

# ----- Initialize Java and Check Version -----
initialize_java <- function() {
  .jinit(parameters = params)
  .jaddClassPath(get_tetrad_path())

  java_version <- safe_jcall("java/lang/System", "S", "getProperty", "java.version")
  print(paste("Java version:", java_version))
}

# ── helper: always return a usable TETRAD class-path ────────────────────
get_tetrad_path <- function() {
  if (exists("TETRAD_PATH", envir = .GlobalEnv) &&
      nzchar(get("TETRAD_PATH", envir = .GlobalEnv))) {
    return(get("TETRAD_PATH", envir = .GlobalEnv))
  }
  ## fall back to the environment variable if the symbol vanished
  Sys.getenv("TETRAD_PATH", unset = "")
}

# make it visible everywhere
assign("get_tetrad_path", get_tetrad_path, envir = .GlobalEnv)

initialize_java()

redirect_java_streams(log_file)

# ── one-stop Java ⇢ DOT helper ───────────────────────────────────────
graph_to_dot <- function(graph) {
  dot <- safe_jcall(
    "edu/cmu/tetrad/graph/GraphSaveLoadUtils",
    "Ljava/lang/String;",
    "graphToDot",
    graph
  )

  if (is.null(dot) || !nzchar(dot)) {
    log_warn("graphToDot() returned empty result.")
    return("")                          # harmless for grViz()
  }
  dot
}

# make it visible to Reference-class methods & parallel workers
assign("graph_to_dot", graph_to_dot, envir = .GlobalEnv)

```

```{r session-setup}
#| context: server
#| echo: false

# --- Buffers --- #
results_buffer <- reactiveVal(data.frame())
results_out_buffer <- reactiveVal(data.frame())
datafile_buffer <- reactiveVal(NULL)
superlearner_output_buffer <- reactiveVal(list())
graphtxt_buffer <- reactiveVal(NULL)
dotfile_buffer <- reactiveVal(NULL)
rv <- reactiveValues()
locked <- reactiveVal(FALSE)
step1lock <- reactiveVal(FALSE)

# ── NEW: keep TMLE rows, ML rows, and the 1-row summary apart ────────────────
tmle_buffer        <- reactiveVal(data.frame())   # 1 row per adjustment set
ml_buffer          <- reactiveVal(data.frame())   # 1 row per ML algorithm
combined_buffer    <- reactiveVal(data.frame())   # exactly 1 row (ribbon plot)
classifier_buffer  <- reactiveVal(NULL)            # all 5 classifier ATEs ("No" mode only)

# create a shared reactiveValues object to pass variables around
rv$dot <- NULL  # define once in app init

# --- BEGIN INLINE: scripts/AIR_functions.R ---
getLocalTags <- function() {
  if (!isLocal()) {
    return(NULL)
  }
  
  htmltools::tagList(
    htmltools::tags$script(paste0(
      "$(function() {",
      "  $(document).on('shiny:disconnected', function(event) {",
      "    $('#ss-connect-dialog').show();",
      "    $('#ss-overlay').show();",
      "  })",
      "});"
    )),
    htmltools::tags$div(
      id="ss-connect-dialog", style="display: none;",
      htmltools::tags$a(id="ss-reload-link", href="#", onclick="window.location.reload(true);")
    ),
    htmltools::tags$div(id="ss-overlay", style="display: none;")
  )
}

isLocal <- function() {
  Sys.getenv("SHINY_PORT", "") == ""
}

disconnectMessage <- function(
    text = "An error occurred. Please refresh the page and try again.",
    refresh = "Refresh",
    width = 450,
    top = 50,
    size = 22,
    background = "white",
    colour = "#444444",
    overlayColour = "black",
    overlayOpacity = 0.6,
    refreshColour = "#337ab7",
    css = ""
) {
  
  checkmate::assert_string(text, min.chars = 1)
  checkmate::assert_string(refresh)
  checkmate::assert_numeric(size, lower = 0)
  checkmate::assert_string(background)
  checkmate::assert_string(colour)
  checkmate::assert_string(overlayColour)
  checkmate::assert_number(overlayOpacity, lower = 0, upper = 1)
  checkmate::assert_string(refreshColour)
  checkmate::assert_string(css)
  
  if (width == "full") {
    width <- "100%"
  } else if (is.numeric(width) && width >= 0) {
    width <- paste0(width, "px")
  } else {
    stop("disconnectMessage: 'width' must be either an integer, or the string \"full\".", call. = FALSE)
  }
  
  if (top == "center") {
    top <- "50%"
    ytransform <- "-50%"
  } else if (is.numeric(top) && top >= 0) {
    top <- paste0(top, "px")
    ytransform <- "0"
  } else {
    stop("disconnectMessage: 'top' must be either an integer, or the string \"center\".", call. = FALSE)
  }
  
  htmltools::tagList(
    getLocalTags(),
    htmltools::tags$head(
      htmltools::tags$style(
        glue::glue(
          .open = "{{", .close = "}}",
          
          "#shiny-disconnected-overlay { display: none !important; }",
          
          "#ss-overlay {
             background-color: {{overlayColour}} !important;
             opacity: {{overlayOpacity}} !important;
             position: fixed !important;
             top: 0 !important;
             left: 0 !important;
             bottom: 0 !important;
             right: 0 !important;
             z-index: 99998 !important;
             overflow: hidden !important;
             cursor: not-allowed !important;
          }",
          
          "#ss-connect-dialog {
             background: {{background}} !important;
             color: {{colour}} !important;
             width: {{width}} !important;
             transform: translateX(-50%) translateY({{ytransform}}) !important;
             font-size: {{size}}px !important;
             top: {{top}} !important;
             position: fixed !important;
             bottom: auto !important;
             left: 50% !important;
             padding: 0.8em 1.5em !important;
             text-align: center !important;
             height: auto !important;
             opacity: 1 !important;
             z-index: 99999 !important;
             border-radius: 3px !important;
             box-shadow: rgba(0, 0, 0, 0.3) 3px 3px 10px !important;
          }",
          
          "#ss-connect-dialog::before { content: '{{text}}' }",
          
          "#ss-connect-dialog label { display: none !important; }",
          
          "#ss-connect-dialog a {
             display: {{ if (refresh == '') 'none' else 'block' }} !important;
             color: {{refreshColour}} !important;
             font-size: 0 !important;
             margin-top: {{size}}px !important;
             font-weight: normal !important;
          }",
          
          "#ss-connect-dialog a::before {
            content: '{{refresh}}';
            font-size: {{size}}px;
          }",
          
          "#ss-connect-dialog { {{ htmltools::HTML(css) }} }"
        )
      )
    )
  )
}

#' Show a nice message when a shiny app disconnects or errors
#'
#' This function is a version of disconnectMessage() with a pre-set combination
#' of parameters that results in a large centered message.
#' @export
disconnectMessage2 <- function() {
  disconnectMessage(
    text = "Your session has timed out.",
    refresh = "",
    size = 70,
    colour = "white",
    background = "rgba(64, 64, 64, 0.9)",
    width = "full",
    top = "center",
    overlayColour = "#999",
    overlayOpacity = 0.7,
    css = "padding: 15px !important; box-shadow: none !important;"
  )
}

# ── tiny validation helpers ──────────────────────────────────────────────
check_numeric <- function(x, min = -Inf, max = Inf, id = "") {
  shiny::validate(
    shiny::need(is.numeric(x) && length(x) == 1 && is.finite(x),
                sprintf("[%s] must be a single numeric value", id)),
    shiny::need(x >= min && x <= max,
                sprintf("[%s] must be between %g and %g", id, min, max))
  )
  x
}

check_column_exists <- function(df, col, id = "") {
  shiny::validate(
    shiny::need(!is.null(df),
                sprintf("Data not loaded, can’t check %s yet", id)),
    shiny::need(col %in% names(df),
                sprintf("Column “%s” not found (%s)", col, id))
  )
  col
}

escape_dot <- function(x) gsub('[^\\w\\s]', '_', x, perl = TRUE)  # DOT-safe

# keeps thresholds inside the data range and numeric
safe_thresh <- function(val, vec, id) {
  rng <- range(vec, na.rm = TRUE)
  check_numeric(val, rng[1], rng[2], id)
}

fix_knowledge <- function(df){
  # Store original column names
  original_colnames <- colnames(df)
  
  # Detect numeric vs non-numeric columns
  numeric_cols <- sapply(df, function(col) all(!is.na(suppressWarnings(as.numeric(as.character(col))))))
  # check if column header is missing and data conform to expectations. If so, process and return
  if (any(!is.na(suppressWarnings(as.numeric(original_colnames))))) {
    
    # Confirm exactly one numeric and one character-type column exist
    if (sum(numeric_cols) == 1 && sum(!numeric_cols) == 1) {
      new_colnames <- c("level", "variable")
      new_colnames_ordered <- rep(NA, length(df))
      new_colnames_ordered[numeric_cols] <- "level"
      new_colnames_ordered[!numeric_cols] <- "variable"
      
      # Move original column names to first row
      df <- rbind(setNames(as.list(original_colnames), names(df)), df)
      
      # Now assign the new column names
      colnames(df) <- new_colnames_ordered
    } else {
      return("Unable to read knowledge file data. Please make sure file contains a header with the following column names: level, variable. 'variable' should contain the name of each variable used, and 'level' should be a numeric value to represent an estimated causal hierarchy (see readme file for a detailed description).")
    }  
  } else if (sum(numeric_cols) == 1 && sum(!numeric_cols) == 1) {
      colnames(df)[numeric_cols] <- "level"
      colnames(df)[!numeric_cols] <- "variable"
    } else {
      return("Unable to read knowledge file data. Please make sure file contains a header with the following column names: level, variable. 'variable' should contain the name of each variable used, and 'level' should be a numeric value to represent an estimated causal hierarchy (see readme file for a detailed description).")
    }
  
  return(df)
}

# change color of nodes in graph
change_node_color <- function(dot_code, node, color) {
  node <- escape_dot(node)
  # Remove any accidental extra quotes from the color string
  color <- trimws(gsub("['\"]", "", color))
  node_definition <- paste0("\"", node, "\" [style=filled, fillcolor=\"", color, "\"];")
  dot_code <- sub("digraph g \\{", paste0("digraph g {\n  ", node_definition), dot_code)
  dot_code <- gsub("\'", "\"", dot_code)
  return(dot_code)
}


paint_nodes <- function(dot, nodes, col) {
  nodes <- nodes[nzchar(nodes)]
  for (n in nodes) {
    dot <- change_node_color(dot, n, col)
  }
  dot
}
AIR_getGraph <- function(data, knowledge, graphtxt_buffer = NULL) {
  headers_string <- "PD\tfrac_ind\tfrac_dep\tunif\t \tBIC\t \t#edges\tn_tests_ind\tn_tests_dep"
  log_debug(headers_string)


r_df_to_tetrad <- function(df) {
  require(rJava)

  log_debug("r_df_to_tetrad(): START")

  n_rows <- nrow(df)
  n_cols <- ncol(df)
  log_debug("Rows: {r}, Cols: {c}", r = n_rows, c = n_cols)

  # ------------------ VARS LIST ------------------
  log_debug("Constructing vars ArrayList...")
  vars <- .jnew("java.util.ArrayList")
  log_debug("vars class = {cls}", cls = .jclass(vars))

  log_debug("Adding ContinuousVariables...")
  for (name in names(df)) {
    log_debug("  Creating variable for column '{nm}'", nm = name)
    v <- .jnew("edu.cmu.tetrad.data.ContinuousVariable", as.character(name))
    log_debug("   OK: created; class={cls}", cls = .jclass(v))

    vars$add(v)
    log_debug("   OK: added to list")
  }

  log_debug("Casting vars → java.util.List...")
  vars_list <- .jcast(vars, "java.util.List")
  log_debug("vars_list class = {cls}", cls = .jclass(vars_list))

  # ------------------ DATA BOX ------------------
  log_debug("Constructing DoubleDataBox...")
  dbox <- .jnew(
    "edu.cmu.tetrad.data.DoubleDataBox",
    as.integer(n_rows),
    as.integer(n_cols)
  )
  log_debug("dbox class = {cls}", cls = .jclass(dbox))

  log_debug("Setting DoubleDataBox values...")
  for (j in seq_len(n_cols)) {
    col <- as.numeric(df[[j]])
    log_debug("  Column {j}: length={len}", j = j, len = length(col))

    for (i in seq_len(n_rows)) {
      # Create a java.lang.Double
      dbl <- .jnew("java.lang.Double", col[i])
      dbox$set(
        as.integer(i - 1L),
        as.integer(j - 1L),
        dbl
      )
    }
  }
  log_debug("Finished populating DoubleDataBox.")

  log_debug("Casting dbox → DataBox...")
  dbox_cast <- .jcast(dbox, "edu.cmu.tetrad.data.DataBox")
  log_debug("dbox_cast class = {cls}", cls = .jclass(dbox_cast))

  # ------------------ DATASET ------------------
  log_debug("Constructing BoxDataSet...")
  dataset <- .jnew(
    "edu.cmu.tetrad.data.BoxDataSet",
    dbox_cast,
    vars_list
  )
  log_debug("dataset created; class = {cls}", cls = .jclass(dataset))

  log_debug("Casting dataset → DataSet interface...")
  dataset_cast <- .jcast(dataset, "edu.cmu.tetrad.data.DataSet")
  log_debug("dataset_cast class = {cls}", cls = .jclass(dataset_cast))

  log_debug("r_df_to_tetrad(): SUCCESS")
  return(dataset_cast)
}

np_dataset_to_df <- function(j_np_dataset, orig_names = NULL) {
  require(rJava)

  log_debug("np_dataset_to_df(): START")

  # Get the underlying Matrix: dataSet.getDoubleData()
  j_mat <- .jcall(
    j_np_dataset,
    "Ledu/cmu/tetrad/util/Matrix;",
    "getDoubleData"
  )
  log_debug("np_dataset_to_df(): got Matrix; class={cls}", cls = .jclass(j_mat))

  n_rows <- .jcall(j_mat, "I", "getNumRows")
  n_cols <- .jcall(j_mat, "I", "getNumColumns")
  log_debug("np_dataset_to_df(): rows={r}, cols={c}", r = n_rows, c = n_cols)

  # Pull values via Matrix.get(int row, int col)
  mat <- matrix(NA_real_, n_rows, n_cols)

  for (i in seq_len(n_rows)) {
    for (j in seq_len(n_cols)) {
      mat[i, j] <- .jcall(
        j_mat,
        "D",
        "get",
        as.integer(i - 1L),
        as.integer(j - 1L)
      )
    }
  }

  df <- as.data.frame(mat)

  if (!is.null(orig_names) && length(orig_names) == n_cols) {
    colnames(df) <- orig_names
  }

  log_debug("np_dataset_to_df(): DONE; returning data.frame with names = {n}",
            n = paste(colnames(df), collapse = ", "))

  df
}




nonparanormal_transform <- function(j_dataset) {

  j_dataset <- .jcast(j_dataset, "edu.cmu.tetrad.data.DataSet")

  .jcall(
    "edu/cmu/tetrad/data/DataTransforms",
    returnSig = "Ledu/cmu/tetrad/data/DataSet;",
    method = "getNonparanormalTransformed",
    j_dataset,
    sig = "(Ledu/cmu/tetrad/data/DataSet;)Ledu/cmu/tetrad/data/DataSet;"
  )
}

# 1. Build Java dataset and NP transform
isbinary <- apply(data, MARGIN = 2, FUN = function(x) {length(unique(x))==2 })
log_debug("nvars: {ib}", ib = paste0(isbinary, collapse = ", "))

data_binary <- data[, isbinary, drop = FALSE]
data_continuous <- data[, !isbinary, drop = FALSE]

log_debug("Discovered {db} binary and {dc} continuous columns...", db = ncol(data_binary), dc = ncol(data_continuous))  

j_dataset    <- r_df_to_tetrad(data_continuous)
log_debug("Class after java transform: {cls}", cls = .jclass(j_dataset))

j_np_dataset <- nonparanormal_transform(j_dataset)
log_debug("Class after NP transform: {cls}", cls = .jclass(j_np_dataset))

# 2. Convert NP dataset back to an R data.frame
#np_df <- np_dataset_to_df(j_np_dataset, orig_names = names(data))
np_df <- np_dataset_to_df(j_np_dataset, orig_names = names(data_continuous))
np_df <- cbind(data_binary, np_df)
log_debug("After NP back to R: class={cls}", cls = paste(class(np_df), collapse = ", "))
# optional:
# log_debug(paste(capture.output(dplyr::glimpse(np_df)), collapse = "\n")) 
  MC_passing_cpdag_already_found <- FALSE
  best_cpdag_seen_so_far <- NULL
  best_cpdag_seen_so_far_num_edges <- Inf
  best_cpdag_seen_so_far_params <- "<unset>"
  cpdag_graph_when_PD_is_1 <- NULL
  last_ts <- NULL

  suppress_output <- function(expr) {
    sink("/dev/null")
    on.exit(sink(), add = TRUE)
    force(expr)
  }

  for (i in seq(0, 15)) {
    pd <- 0.5 + (i * 0.1)

    tryCatch({
      ts <- TetradSearch$new(np_df)
      last_ts <- ts

      # Validate and apply knowledge tiers
      if (!is.null(knowledge) && nrow(knowledge) > 0) {
        if (!all(c("level", "variable") %in% colnames(knowledge))) {
          log_debug("Knowledge file is missing required columns 'level' and 'variable'; skipping tier assignment.")
        } else {
          # Add knowledge to specific tiers
          for (j in seq_len(nrow(knowledge))) {
            level <- knowledge$level[j]
            variable <- knowledge$variable[j]
            
            if (is.numeric(level) && is.character(variable)) {
              tryCatch({
                ts$add_to_tier(level, variable)
              }, error = function(e) {
                log_debug(sprintf("Failed to add variable '%s' to tier %s: %s", 
                                as.character(variable), as.character(level), conditionMessage(e)))
              })
            } else {
              log_debug(sprintf("Invalid knowledge row at index %d: level=%s, variable=%s",
                              j, deparse(level), deparse(variable)))
            }
          }
        }
      }
      
      ts$use_sem_bic(penalty_discount = pd)

      suppress_output(ts$run_boss())
      g2 <- ts$get_java()

      bic <- tryCatch(g2$getAttribute("BIC"), error = function(e) NA_real_)
      num_edges <- tryCatch(g2$getNumEdges(), error = function(e) NA_integer_)

      suppress_output(ts$use_fisher_z(use_for_mc = TRUE))
      ret <- suppress_output(ts$markov_check(g2))

      cpdag_graph_when_PD_is_1 <- g2

      result_line <- tryCatch({
        sprintf("%.1f\t%.4f  \t%.4f   \t%.4f  \t%.2f  \t%.0f  \t%.0f  \t\t%.0f", pd, ret$frac_dep_ind, ret$frac_dep_dep, ret$ad_ind, bic, num_edges, ret$num_tests_ind, ret$num_tests_dep)
      }, error = function(e) {
        "<incomplete>"
      })

      if (!is.null(ret) && is.list(ret) && !is.null(ret$ad_ind) && ret$ad_ind > 0.1) {
        log_debug(result_line)

        if (!MC_passing_cpdag_already_found || (num_edges < best_cpdag_seen_so_far_num_edges)) {
          best_cpdag_seen_so_far <- g2
          best_cpdag_seen_so_far_num_edges <- num_edges
          best_cpdag_seen_so_far_params <- result_line
          MC_passing_cpdag_already_found <- TRUE
        }
      }

      # If PD==1.0 and still no valid CPDAG found, keep something usable
      if (pd == 1.0 && is.null(best_cpdag_seen_so_far)) {
        best_cpdag_seen_so_far <- g2
        best_cpdag_seen_so_far_num_edges <- num_edges
        best_cpdag_seen_so_far_params <- result_line
      }

    }, error = function(e) {
      log_debug(sprintf("Failed iteration for PD=%.1f: %s", pd, conditionMessage(e)))
    })
  }

  # Bail out early if absolutely nothing passed any threshold
  if (!MC_passing_cpdag_already_found && is.null(cpdag_graph_when_PD_is_1)) {
    stop(structure(
      list(
        message = "No valid graph found; TetradSearch may have failed",
        call = match.call()
      ),
      class = c("simpleError", "error", "condition")
    ))
  }

  log_debug("The best cpdag (the one with fewest edges among those for which unif > 0.1) has these attributes:")
  log_debug(headers_string)
  log_debug(best_cpdag_seen_so_far_params)

  graph <- if (!is.null(best_cpdag_seen_so_far)) {
    # A valid CPDAG passed Markov Check
    best_cpdag_seen_so_far
  } else {
    # Nothing passed, but fallback to something rather than crashing
    log_debug("Falling back to cpdag at PD=1.0, as no unif > 0.1 CPDAG was found")
    cpdag_graph_when_PD_is_1
  }

  graph_str <- tryCatch(
    graph_to_dot(graph)
  )

  if (!is.null(graphtxt_buffer)) {
    if (is.null(graph_str) || !nzchar(graph_str)) {
      log_debug("graph_str is NULL or empty, not updating graphtxt_buffer.")
      graphtxt_buffer(NULL)
    } else {
      graphtxt_buffer(
        tryCatch(
          gsub("(?s)Graph Attributes:.*", "", graph_str, perl = TRUE),
          error = function(e) {
            log_debug(sprintf("Regex failed on graph_str: %s", conditionMessage(e)))
            graph_str  # fallback to raw string
          }
        )
      )
    }
  }

  return(list(graph, last_ts, MC_passing_cpdag_already_found, best_cpdag_seen_so_far))
}

AIR_getAdjSets <- function(ts, tv, ov, MC_passing_cpdag_already_found, best_cpdag_seen_so_far){
  TREATMENT_NAME = tv
  RESPONSE_NAME = ov
  MAX_NUM_SETS = 3
  MAX_DISTANCE_FROM_POINT = 4
  MAX_PATH_LENGTH = 4
  NEAR_TREATMENT = 1
  NEAR_RESPONSE = 2
  
  
  log_debug("Identification parameters:")
  log_debug("maximum number of adjustment sets = {p}", p = MAX_NUM_SETS)
  log_debug("maximum distance from target endpoint (TREATMENT or RESPONSE) = {p}", p = MAX_DISTANCE_FROM_POINT)
  log_debug("maximum path length = {p}", p = MAX_PATH_LENGTH)

  len_adj_sets_t <- 0
  len_adj_sets_r <- 0        
  union_of_two_lists <- NULL

  count_undirected_from_dot <- function(dot_text) {
    lines <- strsplit(dot_text, "\n")[[1]]
    undirected <- grep("arrowtail=none[^]]*arrowhead=none", lines, value = TRUE)
    length(undirected)
  }
  undirected_edges <- count_undirected_from_dot(rv$dot)
  log_debug("Number of undirected edges found: {ue}", ue = undirected_edges)

    #if (MC_passing_cpdag_already_found == TRUE) {
    log_debug("Searching for adjustment set(s) on the *** Treatment *** side:")

    # Z1
    adj_sets_treatment = ts$get_adjustment_sets(best_cpdag_seen_so_far, TREATMENT_NAME, RESPONSE_NAME,
                                                MAX_NUM_SETS,
                                                MAX_DISTANCE_FROM_POINT,
                                                NEAR_TREATMENT,
                                                MAX_PATH_LENGTH)
    
    ts$print_adjustment_sets(adj_sets_treatment)
    
    log_debug("Searching for adjustment sets on the *** Response *** side:")
    # Z2
    adj_sets_response  = ts$get_adjustment_sets(best_cpdag_seen_so_far, TREATMENT_NAME, RESPONSE_NAME,
                                                MAX_NUM_SETS,
                                                MAX_DISTANCE_FROM_POINT,
                                                NEAR_RESPONSE,
                                                MAX_PATH_LENGTH)
    
    ts$print_adjustment_sets(adj_sets_response)

    log_debug("adj_sets_treatment type: {type_t}, length: {len_t}",
          type_t = typeof(adj_sets_treatment), len_t = length(adj_sets_treatment))
    log_debug("adj_sets_response type: {type_r}, length: {len_r}",
          type_r = typeof(adj_sets_response), len_r = length(adj_sets_response))
  ## Initialize flag variables to indicate to AIR Step 3 that no/only one adjustment set is yet found
  flag_no_adjustment_set_found = FALSE
  flag_only_one_adjustment_set_found = FALSE
  
  # Determine the union and differences of the two adjustment sets
  union_of_two_lists = union(adj_sets_treatment, adj_sets_response)
  near_treatment_not_near_response = setdiff(adj_sets_treatment, adj_sets_response)
  near_response_not_near_treatment = setdiff(adj_sets_response, adj_sets_treatment)

  for (i in seq_along(adj_sets_treatment)) {
       set <- adj_sets_treatment[[i]]
      
       if (length(set) != 0) {
        len_adj_sets_t <- len_adj_sets_t + 1
       }
   }
   log_debug("Length of non-zero adjustment sets for treatment {set}", set = len_adj_sets_t)

    

  for (i in seq_along(adj_sets_response)) {
       set <- adj_sets_response[[i]]
      
       if (length(set) != 0) {
        len_adj_sets_r <- len_adj_sets_r + 1
       }
   }
   log_debug("Length of non-zero adjustment sets for response {set}", set = len_adj_sets_r)  

  
  log_debug("Total number of adjustment sets encountered (ignoring duplicates) = {msg}", msg = length(union_of_two_lists))
  log_debug("Size of Treatment - Response adjustment sets = {msg}", msg = length(near_treatment_not_near_response))
  log_debug("Size of Response - Treatment adjustment sets = {msg}", msg = length(near_response_not_near_treatment))


   #} 
   if (MC_passing_cpdag_already_found == FALSE) {
    #union_of_two_lists <- NULL
    log_debug("MC_passing_cpdag_already_found = FALSE")
    log_debug("*** WARNING: Markov Condition not satisfied for current graph! Proceed with caution. ***")
  } 
 
  
  log_debug("union of two lists = {union}", union = paste(capture.output(str(union_of_two_lists)), collapse = "\n"))
  log_debug("target adjustment sets: {tar}", tar = len_adj_sets_t)
  log_debug("response adjustment sets: {res}", res = len_adj_sets_r)
  ## Now consider all cases where either set (or both sets) of adjustment sets is (are) empty.
  # if both empty, we need to set the corresponding flag:
  if (length(union_of_two_lists) == 0 | (len_adj_sets_t + len_adj_sets_r == 0)) {
    # Then no adjustment sets and we must be working with a cpdag (at least one undirected
    #   edge) rather than a DAG. (For a DAG, there's always an adjustment set--namely the
    #   set of parents of the treatment variable, which can be empty, but that's still
    #   a valid adjustment set.) There are multiple solutions here, but get the end-user
    #   involved.
    if (undirected_edges == 0){
      log_debug("No adjustment set found. The structure of this data is not amenable to causal inference as no adjustment sets can be found.")
  
      flag_no_adjustment_set_found = TRUE    
      return(list("error","structure"))
    } else if (undirected_edges > 0){
      log_debug("No adjustment set found. Revise knowledge (see AIR job aid) so search result has no undirected edges")
  
      flag_no_adjustment_set_found = TRUE    
      return(list("error","edges"))
    }

  } else if (length(union_of_two_lists)==1) {
    # Only one adjustment set is found altogether, and so we set the corresponding flag:
    flag_only_one_adjustment_set_found = TRUE
    return_first_adj_set  = union_of_two_lists[[1]]
    return_second_adj_set = union_of_two_lists[[1]]  # same adjustment set
  } else {
    # We have at least two distinct adjustment sets; we prefer ones only from each side if practical, so, test alternatives first.
    log_debug("At least two adjustment sets found.")
    
    # if no adjustment set found near response that was not already near treatment:
    if (length(near_response_not_near_treatment) == 0) {
      log_debug("In this case, there is no adjustment set near response that is not also near treatment.")
      return_first_adj_set  = adj_sets_treatment[[1]]
      return_second_adj_set = adj_sets_treatment[[2]]
    } else if (length(near_treatment_not_near_response) == 0) {
      log_debug("In this case, there is no adjustment set near treatment that is not also near response.")
      return_first_adj_set  = adj_sets_response[[1]]
      return_second_adj_set = adj_sets_response[[2]]
    } else {
      # At least one adjustment set is near treatment but not response, and vice versa.
      #   This is the ideal case to obtain greater diversity of adjustment sets, which 
      #   is also why we might want to set max_num_sets higher.
      #   Return two distinct adjustment sets, one from each side.
      log_debug("In this case, we have found at least one adjustment set near treatment but not near response; and vice versa.")
      return_first_adj_set  = near_treatment_not_near_response[[1]]
      return_second_adj_set = near_response_not_near_treatment[[1]]
    }
  }
  
  #if (MC_passing_cpdag_already_found == TRUE) {
   log_debug("Summary of results:")
   log_debug("First adjustment set to return: {msg}", msg = return_first_adj_set)
   log_debug("Second adjustment set to return: {msg}", msg = return_second_adj_set)
   log_debug("Flag status for no adjustment set found: {msg}", msg = flag_no_adjustment_set_found)
   log_debug("Flag status for only one adjustment set found: {msg}", msg = flag_only_one_adjustment_set_found)
  
   ### return the two adjustment sets plus the two flags instead.
   return(list(return_first_adj_set, return_second_adj_set))
  #} else {return(NULL)}
}

scale_ <- function(x){
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

runSuperLearner <- function(settings,
                            AIRHome,
                            data_df,
                            tv_dir,
                            ov_dir,
                            tv_threshold,
                            ov_threshold,
                            log_file){ 
  log_debug("Started super-learner with {line}", line = paste(c(settings, AIRHome, tv_threshold, ov_threshold), collapse = ', '))
  Z_level <- settings$Z_level
  doc_title <- settings$doc_title
  
  log_debug("Reading in data")
  confounders <- as.character(strsplit(x = settings$confounders, split = " ")[[1]])
  treatment <- as.character(settings$varName)
  outcome =  as.character(df_vars[df_vars$var == "OV",]$val)

  df <- data_df
  if (isTRUE(tv_dir)){ df[[treatment]] <- ifelse(df[[treatment]] >= tv_threshold, 1, 0)} else { df[[treatment]] <- ifelse(df[[treatment]] <= tv_threshold, 1, 0)}
  if (isTRUE(ov_dir)){ df[[outcome  ]] <- ifelse(df[[outcome  ]] >= ov_threshold, 1, 0)} else { df[[outcome  ]] <- ifelse(df[[outcome  ]] <= ov_threshold, 1, 0)}

  ## ────────────────────────────────────────────────────────────────
  ##  Coerce the *outcome* to numeric if it’s still a factor
  ##  (TMLE + sl3 learners expect a numeric Y for a “continuous”
  ##   outcome; leaving it as factor triggers the warning you saw)
  ## ────────────────────────────────────────────────────────────────
  if (is.factor(df[[outcome]])) {
    log_debug(sprintf("Outcome %s is a factor – coercing to numeric {0,1}.", outcome))

    ## factor levels are usually 1 / 2 after the binarisation above,
    ## so subtract 1 to make them {0,1}
    df[[outcome]] <- as.numeric(df[[outcome]]) - 1
  }

  #### TMLE -------------------------------------
  ##### Define Superlearner -------------------
  # sl3_list_learners("binomial") 
  
  log_debug("Building learner list")
  lrnr_mean <- sl3::make_learner(sl3::Lrnr_mean)
  lrnr_glm <- sl3::make_learner(sl3::Lrnr_glm)
  lrnr_hal <- sl3::make_learner(sl3::Lrnr_hal9001)
  lrnr_nnet <- sl3::make_learner(sl3::Lrnr_nnet)

  ## sl3 ≥ 0.3.0:  Lrnr_ranger            (C++ backend, seed-respecting)
  ## sl3 ≤ 0.2.x:  Lrnr_randomForest      (R backend)
  lrnr_rf  <- if ("Lrnr_ranger" %in% getNamespaceExports("sl3"))
                sl3::make_learner(sl3::Lrnr_ranger)
              else
                sl3::make_learner(sl3::Lrnr_randomForest)

  lrnr_glmnet <- sl3::make_learner(sl3::Lrnr_glmnet)
  lrnr_xgboost <- sl3::make_learner(sl3::Lrnr_xgboost, max_depth = 4, eta = 0.01, nrounds = 100)
  lrnr_earth <- sl3::make_learner(sl3::Lrnr_earth)
  if (length(confounders) > 1) {
    sl_ <- sl3::make_learner(sl3::Stack, unlist(list(lrnr_mean, 
                                           lrnr_glm,
                                           lrnr_hal,
                                           lrnr_rf,
                                           lrnr_glmnet,
                                           lrnr_xgboost,
                                           lrnr_earth,
                                           lrnr_nnet),
                                      recursive = TRUE))
  } else {
    sl_ <- sl3::make_learner(sl3::Stack, unlist(list(lrnr_mean, 
                                           lrnr_glm,
                                           lrnr_hal,
                                           lrnr_rf,
                                           lrnr_xgboost,
                                           lrnr_earth,
                                           lrnr_nnet), 
                                      recursive = TRUE))
  }
  # DEFINE SL_Y AND SL_A 
  # We only need one, because they're the same
  ##### Define Formulae --------------------------
  Q_learner <- sl3::Lrnr_sl$new(learners = sl_, 
                           metalearner = sl3::Lrnr_nnls$new(convex = T)) # output model
  g_learner <- sl3::Lrnr_sl$new(learners = sl_, 
                           metalearner = sl3::Lrnr_nnls$new(convex = T)) # treatment model
  learner_list <- list(Y = Q_learner,
                       A = g_learner)
  
  # PREPARE THE THINGS WE WANT TO FEED IN TO TMLE3
  ate_spec <- tmle3::tmle_ATE(treatment_level = 1, control_level = 0)
  

  
  ##### Nodes ------------------
  nodes_ <- list(W = confounders, # covariates
                 A = treatment,
                 # Z = mediators, # unnecessary unless doing mediation analysis
                 Y = outcome)
  
  ##### RUN TMLE3 -------------------------------
  set.seed(123)
  ### this is where the parallel is breaking
  log_debug("Starting TMLE")
  tryCatch({
    tmle_fit_ <- tmle3::tmle3(tmle_spec = ate_spec,
                 data = df,
                 node_list = nodes_,
                 learner_list = learner_list)
  }, error = function(e) {
    log_debug("Error in tmle3 call: {msg}", msg = conditionMessage(e))
    stop(e)
  })
  log_debug("Pulling out TMLE scores")
  tmle_task <- ate_spec$make_tmle_task(df, nodes_)
  
  initial_likelihood <- ate_spec$make_initial_likelihood(
    tmle_task,
    learner_list
  )
  
  ## save propensity score for diagnosis
  propensity_score <- initial_likelihood$get_likelihoods(tmle_task)$A
  propensity_score <- propensity_score * df[,..treatment] + (1 - propensity_score) * (1 - df[,..treatment])
  
  plap_ <- tibble(exposure = df[,..treatment] |> pull(),
                  pscore = propensity_score |> pull())
  
  plap_$sw <- plap_$exposure * (mean(plap_$exposure)/propensity_score) + (1 - plap_$exposure) * ((1 - mean(plap_$exposure)) / (1 - propensity_score))

  ## inside runSuperLearner() – after tmle_fit_ is available
  prefix <- tolower(Z_level)    # "z1"  or  "z2"

  ## Extract the three CI numbers *as vectors*
  est  <- as.numeric(tmle_fit_$summary[[8]])
  lci  <- as.numeric(tmle_fit_$summary[[9]])
  uci  <- as.numeric(tmle_fit_$summary[[10]])

  tmle_row <- tibble::tibble(
    Treatment = treatment,
    Group     = Z_level,
  
    ## fill the columns that belong to *this* adjustment set …
    z1_Mean = if (prefix == "z1") est else NA_real_,
    z1_LCI  = if (prefix == "z1") lci else NA_real_,
    z1_UCI  = if (prefix == "z1") uci else NA_real_,
  
    z2_Mean = if (prefix == "z2") est else NA_real_,
    z2_LCI  = if (prefix == "z2") lci else NA_real_,
    z2_UCI  = if (prefix == "z2") uci else NA_real_
  )

  # save outcome predictions for diagnosis
  # initial_likelihood_preds was formerly labeled outcome_preds
  initial_likelihood_preds <- initial_likelihood$get_likelihoods(tmle_task,"Y")
  # define and fit likelihood
  factor_list <- list(
    tmle3::define_lf(LF_emp, "W"),
    tmle3::define_lf(LF_fit, "A", sl_),
    tmle3::define_lf(LF_fit, "Y", sl_, type = "mean")
  )
  likelihood_def <- tmle3::Likelihood$new(factor_list)
  likelihood <- likelihood_def$train(tmle_task)
  likelihood_values <- rowMeans(likelihood$get_likelihoods(tmle_task,"Y"))
  
  # print("super learner coefficients for PS model")
  g_fit <- tmle_fit_$likelihood$factor_list[["A"]]$learner
  # g_fit$fit_object$full_fit$learner_fits$Lrnr_nnls_TRUE
  
  # print("super learner coefficients for outcome model")
  Q_fit <- tmle_fit_$likelihood$factor_list[["Y"]]$learner
  # Q_fit$fit_object$full_fit$learner_fits$Lrnr_nnls_TRUE
  
  ## generate counterfactuals
  ### counterfactual where all treatments set to 1
  intervention1 <- tmle3::define_lf(LF_static, "A", value = 1)
  
  cf_likelihood1 <- tmle3::make_CF_Likelihood(likelihood, intervention1)
  
  cf_likelihood_values1 <- cf_likelihood1$get_likelihoods(tmle_task, "A")
  
  # We can then use this to construct a counterfactual likelihood:
  ### counterfactual where all treatments set to 0
  # set values
  intervention0 <- tmle3::define_lf(LF_static, "A", value = 0)
  # generate counterfactual likelihood object
  cf_likelihood0 <- tmle3::make_CF_Likelihood(likelihood, intervention0)
  # get likelihoods from object
  cf_likelihood_values0 <- cf_likelihood0$get_likelihoods(tmle_task, "A")
  # We see that the likelihood values for the A node are all either 0 or 1, as would be expected from an indicator likelihood function. In addition, the likelihood values for the non-intervention nodes have not changed.
  log_debug("Building Output")
  ## output individual row values
  # df_out <- df[,c(nodes_$A, nodes_$Y, nodes_$W)]
  df_out <- df |> select(nodes_$A, nodes_$Y, nodes_$W)
  df_out$exposure <- plap_$exposure
  df_out$rownum <- rownames(df_out)
  df_out$pscore <- plap_$pscore
  df_out$sw <- plap_$sw
  df_out$tmle_est <- tmle_fit_$estimates[[1]]$IC
  df_out$initial_likelihood_preds <- initial_likelihood_preds
  df_out$likelihood_values <- likelihood_values
  df_out$counterfactual_0 <- cf_likelihood_values0
  df_out$counterfactual_1 <- cf_likelihood_values1
  df_out$g_fit_pred <- g_fit$predict() 
  df_out$Q_fit_pred <- Q_fit$predict()
  log_debug("Finished SuperLearner – returning results")

  # -------- FINAL RETURN (nothing after this executes) ------------
  list(
    tmle = tmle_row,                                # one‐row summary
    diag = list(doc_title = doc_title, data = df_out)  # diagnostics
  )
}

processResults <- function(settings,
                           AIRHome,
                           tv_dir,
                           ov_dir,
                           tv_threshold,
                           ov_threshold,
                           model_in, model_yn,
                           model_ate,
                           log_file,
                           move_results = FALSE) {
  log_debug("Started processResults()")
  log_debug("Settings: {line}",
    line = paste(c(settings, AIRHome, tv_dir, ov_dir, tv_threshold, ov_threshold, model_yn,  model_ate), collapse = ", "))

  treatment <- as.character(settings$varName)
  outcome <- as.character(df_vars[df_vars$var == "OV", ]$val)
  confounders <- as.character(unique(Zvars$Z))
  doc_title <- settings$doc_title

  df <- req(datafile_buffer())
  set.seed(123)

  # Bin conversion logic
  if (isTRUE(tv_dir)){ df[[treatment]] <- ifelse(df[[treatment]] >= tv_threshold, 1, 0)} else { df[[treatment]] <- ifelse(df[[treatment]] <= tv_threshold, 1, 0)}
  if (isTRUE(ov_dir)){ df[[outcome  ]] <- ifelse(df[[outcome  ]] >= ov_threshold, 1, 0)} else { df[[outcome  ]] <- ifelse(df[[outcome  ]] <= ov_threshold, 1, 0)}


  # Split into train/test
  test_size <- floor(0.3 * nrow(df))
  samp <- sample(nrow(df), test_size)
  ## ML learners downstream need a *numeric* 0/1 response.  Coerce once here
  ## and keep the original column name.
  y_train <- df[-samp, outcome, drop = FALSE] |>
             mutate(!!outcome := as.numeric(.data[[outcome]]))
  x_train <- df[-samp, setdiff(names(df), outcome), drop = FALSE]
  y_test  <- df[samp,  outcome, drop = FALSE] |>
             mutate(!!outcome := as.numeric(.data[[outcome]]))
  x_test <- df[samp, setdiff(names(df), outcome), drop = FALSE]
  train <- cbind(label = y_train[[1]], x_train)
  test <- cbind(label = y_test[[1]], x_test)
  colnames(train)[1] <- "label"
  colnames(test)[1] <- "label"
  xtest_0 <- mutate(x_test, !!treatment := 0)
  xtest_1 <- mutate(x_test, !!treatment := 1)

  log_info("Data loaded and split.")

  if (model_yn == "No") {
    log_info("Running new ML models")

    ## --- helper: a “safe” GLM that gracefully falls back to lm() -------------
    safe_glm <- function(formula, data, ...) {
      tryCatch(
        {
          suppressWarnings(                     # keep harmless warnings quiet
            glm(formula, data = data, family = binomial, ...)
          )
        },
        warning = function(w) {
          log_info(sprintf("⚠️  glm(): %s – falling back to linear regression.", conditionMessage(w)))
          lm(formula, data = data)              # fall back → numeric preds
        },
        error = function(e) {
          log_info(sprintf("❌  glm() failed: %s – using lm() instead.", conditionMessage(e)))
          lm(formula, data = data)
        }
      )
    }

    ## all four learners now see a numeric 0/1 outcome
    model_lm  <- safe_glm(label ~ ., data = train)
    ## --- helper: run rpart in "class" mode when the response is 0/1 ----------
    safe_rpart <- function(formula, data, ...) {
      y <- data[[all.vars(formula)[1L]]]
      is_binary <- is.numeric(y) && setequal(unique(y), c(0, 1))
      method    <- if (is_binary) "class" else "anova"

      tryCatch(
        {
          suppressWarnings(
            rpart::rpart(formula, data = data, method = method, ...)
          )
        },
        warning = function(w) {
          log_info(sprintf("⚠️  rpart(): %s – continuing with produced tree.", conditionMessage(w)))
          invokeRestart("muffleWarning")        # keep going
        }
      )
    }

    ## --- helper: pick the right SVM mode & silence benign warnings --------
    safe_svm <- function(formula, data, ...) {
      y <- data[[all.vars(formula)[1L]]]
      is_binary <- is.numeric(y) && setequal(unique(y), c(0, 1))
      svm_type  <- if (is_binary) "C-classification" else "eps-regression"

      tryCatch(
        {
          suppressWarnings(
            e1071::svm(formula,
                       data = data,
                       type = svm_type,
                       ...)
          )
        },
        warning = function(w) {
          log_info(sprintf("⚠️  svm(): %s – continuing with produced model.", conditionMessage(w)))
          invokeRestart("muffleWarning")
        }
      )
    }

    model_dt  <- safe_rpart(label ~ ., data = train)               # tree
    model_svm <- safe_svm(label ~ ., data = train)                 # SVM
    model_rf  <- randomForest(label ~ ., data = train)

    pred_lm0 <- predict(model_lm, xtest_0)
    pred_lm1 <- predict(model_lm, xtest_1)
    pred_dt0 <- predict(model_dt, xtest_0)
    pred_dt1 <- predict(model_dt, xtest_1)

    ## SVM: use probabilities when we trained a classifier (otherwise numeric)
    if (model_svm$type == "C-classification") {
      pred_svm0 <- attr(predict(model_svm, xtest_0, probability = TRUE),
                        "probabilities")[, "yes"]
      pred_svm1 <- attr(predict(model_svm, xtest_1, probability = TRUE),
                        "probabilities")[, "yes"]
    } else {
      pred_svm0 <- as.numeric(predict(model_svm, xtest_0))
      pred_svm1 <- as.numeric(predict(model_svm, xtest_1))
    }
    pred_rf0 <- predict(model_rf, xtest_0)
    pred_rf1 <- predict(model_rf, xtest_1)

    lm_ate <- mean(pred_lm1) - mean(pred_lm0)
    dt_ate <- mean(pred_dt1) - mean(pred_dt0)
    svm_ate <- mean(as.numeric(pred_svm1)) - mean(as.numeric(pred_svm0))
    rf_ate <- mean(as.numeric(pred_rf1)) - mean(as.numeric(pred_rf0))

    covariates <- setdiff(names(train), "label")
    ## choose the proper outcome_type ------------------------------------------
    out_type <- if (is.numeric(train$label) &&
                    setequal(unique(train$label), c(0, 1)))
                  "binomial" else "continuous"

    task <- sl3::make_sl3_Task(data         = train,
                               outcome      = "label",
                               covariates   = covariates,
                               outcome_type = out_type)
    learners <- list(
      sl3::Lrnr_glm$new(),
      sl3::Lrnr_glmnet$new(),
      sl3::Lrnr_xgboost$new(),
      sl3::Lrnr_earth$new(),
      sl3::Lrnr_nnet$new(),
      sl3::Lrnr_svm$new()
    )
    sl <- sl3::Lrnr_sl$new(
      learners     = sl3::Stack$new(learners),
      metalearner  = sl3::Lrnr_nnls$new()
    )
    sl_fit <- sl$train(task)

    ## -- helper: always hand rowMeans() a 2-D object ------------------------
    make_pred <- function(task) {
      p <- sl_fit$predict(task)
      if (is.null(dim(p)))               # <- only one learner returned
        p <- matrix(p, ncol = 1)         #    promote to 1-col matrix
      rowMeans(p, na.rm = TRUE)
    }

    pred_sl0 <- make_pred(sl3::make_sl3_Task(data = xtest_0,
                                             covariates = covariates))
    pred_sl1 <- make_pred(sl3::make_sl3_Task(data = xtest_1,
                                             covariates = covariates))
    sl_ate <- mean(pred_sl1) - mean(pred_sl0)

    log_info("Model training complete")

    ## five learners → PICK ONE that represents “the classifier”
    ate_df_all <- data.frame(
      algorithm = c("Logistic Regression", "Decision Tree", "Support Vector Machine",
                    "Random Forest", "ML Model Summary"),
      flag      = c(lm_ate, dt_ate, svm_ate, rf_ate, sl_ate),
      Treatment = outcome
    )

    ## keep only the stacked-SL row so we have ONE reference ATE
    ate_df_ref <- dplyr::filter(ate_df_all,
                                algorithm == "ML Model Summary")

    classifier_buffer(ate_df_all)   # all 5 ATEs — used for spread display in interpretation
    ml_buffer(ate_df_ref)           # stacked-SL summary is the reference ATE

  } else {
    log_info(sprintf("Using existing model config: %s", model_yn))

    if (model_yn == "Yes") {
      if (is.null(model_in)) {
        log_info("⚠️  No model file was supplied – skipping custom-model ATE.")
        return(invisible(NULL))            # bail out, nothing downstream to do
      }
  
      # pred_m0 <- predict(model_in, xtest_0)
      # pred_m1 <- predict(model_in, xtest_1)
      # m_ate   <- mean(as.numeric(pred_m1)) - mean(as.numeric(pred_m0))
      pred_m0 <- tryCatch(predict(model_in, xtest_0, type = "prob"), error = function(e) NULL)
      pred_m1 <- tryCatch(predict(model_in, xtest_1, type = "prob"), error = function(e) NULL)

      if (!is.null(pred_m0) && is.matrix(pred_m0)) {
        # use the probability of the "positive" class (last column by convention)
        pred_m0 <- pred_m0[, ncol(pred_m0)]
        pred_m1 <- pred_m1[, ncol(pred_m1)]
        log_debug("Using predicted probabilities from classification model.")
      } else {
        # fallback for regression-style models
        pred_m0 <- predict(model_in, xtest_0)
        pred_m1 <- predict(model_in, xtest_1)
      }
      m_ate <- mean(as.numeric(pred_m1)) - mean(as.numeric(pred_m0))
      log_debug("ATE: mean(pred_m1)={p1}, mean(pred_m0)={p0}, diff={d}",
            p1 = mean(as.numeric(pred_m1)),
            p0 = mean(as.numeric(pred_m0)),
            d  = m_ate)
      } else if (model_yn == "ATE") {
        m_ate <- model_ate
      } else if (model_yn == "null"){
        m_ate <- 0
      }

    ate_df <- data.frame(
      algorithm = "Existing Model",
      flag      = m_ate,
      Treatment = outcome
    )

    ml_buffer(ate_df)
  }

  # ---------------------------------------------------------------------------
  # small helpers
  # ---------------------------------------------------------------------------
  safe_drop <- function(.data, cols) {
    dplyr::select(.data, -dplyr::any_of(cols))
  }
  
  force_numeric_cols <- function(.data, cols) {
    missing <- setdiff(cols, names(.data))
    if (length(missing) > 0) .data[missing] <- NA_real_
    dplyr::mutate(.data, dplyr::across(all_of(cols), as.numeric))
  }
  
  # ---------------------------------------------------------------------------
  # 1. pull whatever the workers have written so far --------------------------
  # ---------------------------------------------------------------------------
  results <- results_buffer() |> tibble::as_tibble()
  
  ## ── NEW: bail out if *all* learners failed ─────────────────────────────
  if (nrow(results) == 0 ||
      purrr::every(
        results[c("z1_Mean","z1_LCI","z1_UCI",
                  "z2_Mean","z2_LCI","z2_UCI")],
        ~ all(is.na(.x))
      )) {
    log_info("❌ All TMLE calls returned NA — nothing to summarise.")
    shiny::showNotification(
      "All adjustment-set learners failed – no ATEs to display.",
      type = "error", duration = 8
    )
    return(invisible(NULL))   # abort before any plotting math
  }
  
  # ---------------------------------------------------------------------------
  # 2. schema guards – be SURE the numeric columns are really numeric ----------
  # ---------------------------------------------------------------------------
  if (!"Group" %in% names(results)) results <- dplyr::mutate(results, Group = "z")
  
  # keep Treatment as character; coerce everything else (except Group)
  results <- results |>
    dplyr::mutate(
      Treatment = as.character(Treatment),
      dplyr::across(-c(Group, Treatment), ~ as.numeric(unlist(.x)))
    )
  
  log_info(
    "results snapshot >>>\n",
    paste(
      utils::capture.output({
        print(utils::head(results, 6))
        str(results)
      }),
      collapse = "\n"
    ),
    "\n<<< end snapshot"
  )
  
  # ---------------------------------------------------------------------------
  # 3. (A) build the *causal* results frame – **ML stays outside**
  # ---------------------------------------------------------------------------
  causal_df <- results |>                       # **only TMLE rows here**
    dplyr::mutate(
      Treatment = outcome %||% Treatment        # ensure present / non-NA
    ) |>
    ## rename ONLY *after* the sig-test so those columns still exist
    dplyr::rename(
      z1_ATE     = z1_Mean,
      z1_ATE_LCI = z1_LCI,
      z1_ATE_UCI = z1_UCI,
      z2_ATE     = z2_Mean,
      z2_ATE_LCI = z2_LCI,
      z2_ATE_UCI = z2_UCI
    )

  # ---------------------------------------------------------------------------
  # 3. (B) bind in the *ML* ATEs and derive significance flags
  #       – this is now a clean join of two equally-shaped tables
  # ---------------------------------------------------------------------------
  ml_df <- ml_buffer()                           # <- comes from runSuperLearner()
  
  final_results <- dplyr::full_join(causal_df, ml_df,
                                    by = "Treatment") |>
    dplyr::mutate(
      z1_sig = ifelse(
        !is.na(z1_ATE_LCI) & !is.na(z1_ATE_UCI) &
        (flag < z1_ATE_LCI | flag > z1_ATE_UCI),
        1, 0
      ),
      z2_sig = ifelse(
        !is.na(z2_ATE_LCI) & !is.na(z2_ATE_UCI) &
        (flag < z2_ATE_LCI | flag > z2_ATE_UCI),
        1, 0
      ),
      significance = pmax(z1_sig, z2_sig, na.rm = TRUE)
    )

  ## ── NEW: extra guard & verbose logging ──────────────────────────────
  if (!"Treatment" %in% names(final_results)) {
    log_info("⚠️  Treatment column absent – creating it from settings$varName.")
    final_results <- dplyr::mutate(final_results,
                                   Treatment = settings$varName)
  } else {
    final_results <- dplyr::mutate(final_results,
                                   Treatment = tidyr::replace_na(Treatment,
                                                                 settings$varName))
  }

  # crash loudly *here* (where the problem is easy to see) if anything
  # critical is still NA/NULL
  req_cols <- c("Treatment", "z1_ATE_LCI", "z1_ATE_UCI",
                "z2_ATE_LCI", "z2_ATE_UCI")
  missing_vals <- purrr::keep(final_results[req_cols], ~ all(is.na(.x)))
  if (length(missing_vals) > 0) {
    stop("processResults(): critical columns still NA – see log for details")
  }

  # log a concise snapshot for post-mortem debugging
  log_info(
    "final_results snapshot →\n",
    paste(
      utils::capture.output({
        print(utils::head(final_results, 6))
        str(final_results)
      }),
      collapse = "\n"
    )
  )
  
  # ---------------------------------------------------------------------------
  # 4. provide a one-row "combined" summary the UI helpers expect -------------
  # ---------------------------------------------------------------------------
  combined_row <- final_results |>
    summarise(across(
      c(z1_ATE, z1_ATE_LCI, z1_ATE_UCI,
        z2_ATE, z2_ATE_LCI, z2_ATE_UCI,
        flag, z1_sig, z2_sig, significance),
      ~ first(na.omit(.x))
    )) |>
    mutate(Group = "combined",
           Treatment = first(final_results$Treatment),
           algorithm = first(final_results$algorithm))
  
  # ---------------------------------------------------------------------------
  # 5. stash both flavours for later use --------------------------------------
  #    • the detailed two-row frame (if you need it somewhere else)
  #    • the one-row “combined” frame the plots/text expect
  # ---------------------------------------------------------------------------
  log_debug("results_out_buffer: {cr}", cr = jsonlite::toJSON(combined_row))
  results_out_buffer(combined_row)
  assign("results_detailed", final_results, envir = .GlobalEnv)   # <<- optional
  
  log_info("Results processed and stored.")
}

## Create four simple methods to help parse a line

# ── DOT-parsing helpers (quoted edges) ──────────────────────────────────
# Examples of the lines we want to catch:
#   "mission_urgency" -> "region_sensitivity";
#   "A" -> "B";

## 1. does this line specify an edge?  ───────────────────────────────────
is_edge_line <- function(line) {
  grepl('^\\s*"[^"]+"\\s*->\\s*"[^"]+"', line)
}

## 2. split the line into the two node names  ────────────────────────────
split_line <- function(input_string) {
  # grab everything that is inside matching double quotes
  unlist(regmatches(input_string,
                    gregexpr('"([^"]+)"', input_string, perl = TRUE))) |>
    gsub('"', '', x = _ , fixed = TRUE)          # strip the quotes
}

## 3. first node (source)  ───────────────────────────────────────────────
first_node <- function(line) {
  split_line(line)[1]
}

## 4. second node (target) ───────────────────────────────────────────────
second_node <- function(line) {
  split_line(line)[2]
}

## The key methods defined here are descendants (and its dual: ancestors)
descendants <- function(node, children) {
  checked_so_far <- set()
  seen_not_checked <- set(node)
  while (!set_is_empty(seen_not_checked)) {
    check_these <- seen_not_checked
    for (n in check_these) {
      for (m in children[[n]]) {
        if ( !(m %e% checked_so_far)) {
          seen_not_checked <- seen_not_checked | set(m)
        }
      }
      seen_not_checked <- seen_not_checked - set(n)
      checked_so_far <- checked_so_far | set(n)
    }
  }
  return(checked_so_far)
}

ancestors <- function(node, children){
  checked_so_far <- set()
  seen_not_checked <- set(node)
  while (!set_is_empty(seen_not_checked)) {
    check_these <- seen_not_checked
    for (n in check_these) {
      for (m in parents[[n]]) {
        if ( !(m %e% checked_so_far)) {
          seen_not_checked <- seen_not_checked | set(m)
        }
      }
      seen_not_checked <- seen_not_checked - set(n)
      checked_so_far <- checked_so_far | set(n)
    }
  }
  return(checked_so_far)}

get_orphans <- function(graphtxt_buffer) {
  buf <- graphtxt_buffer()
  if (is.null(buf) || !nzchar(buf)) return(character(0))

  lines <- strsplit(buf, "\n", fixed = TRUE)[[1]]
  edge_flags <- vapply(lines, is_edge_line, logical(1))
  edge_lines <- lines[edge_flags]

  if (!length(edge_lines)) return(character(0))

  parents  <- vapply(edge_lines, first_node,  character(1), USE.NAMES = FALSE)
  children <- vapply(edge_lines, second_node, character(1), USE.NAMES = FALSE)

  nodes <- unique(c(parents, children))

  # in-degree via tabulation of "child" occurrences
  indeg <- table(children)
  has_parent <- names(indeg)

  orphans <- setdiff(nodes, has_parent)
  return(sort(orphans))
} 

get_spinsters <- function(graphtxt_buffer) {
  buf <- graphtxt_buffer()
  if (is.null(buf) || !nzchar(buf)) return(character(0))

  lines <- strsplit(buf, "\n", fixed = TRUE)[[1]]
  edge_flags <- vapply(lines, is_edge_line, logical(1))
  edge_lines <- lines[edge_flags]
  if (!length(edge_lines)) return(character(0))

  parents  <- vapply(edge_lines, first_node,  character(1), USE.NAMES = FALSE)
  children <- vapply(edge_lines, second_node, character(1), USE.NAMES = FALSE)

  nodes <- unique(c(parents, children))

  outdeg <- table(parents)           # how many times each node appears as a parent
  has_child <- names(outdeg)

  spinsters <- setdiff(nodes, has_child)  # nodes that never appear as a parent
  sort(spinsters)
}

get_X_descendents <- function(TV, graphtxt_buffer) {
  buf <- graphtxt_buffer()
  if (is.null(buf) || !nzchar(buf))           # <<< guard added
    return(character(0))

  lines <- strsplit(buf, "\n")[[1]]
  
  nodes <- set()
  for (line in lines) {
    if (is_edge_line(line)) {
      if (!(first_node(line)  %e% nodes)) {
        nodes <- nodes | set(first_node(line))
      }
      if (!(second_node(line) %e% nodes)) {
        nodes <- nodes | set(second_node(line))
      }
    }
  }

  ## Create children dictionaries
  children <- hash()
  for (n in nodes) {
    children_of_n <- set()
    for (line in lines) {
      if (is_edge_line(line)) {
        if ((first_node(line)  == n) && (!(second_node(line) %e% children_of_n))) {
          children_of_n <- children_of_n | set(second_node(line))
        }
      }
    }
    children[[n]] <- children_of_n
  }
  return(unlist(descendants(TV, children)))
}

get_X_ancestors <- function(TV, graphtxt_buffer) {
  buf <- graphtxt_buffer()
  if (is.null(buf) || !nzchar(buf))           # <<< guard added
    return(character(0))

  lines <- strsplit(buf, "\n")[[1]]
  
  nodes <- set()
  for (line in lines) {
    if (is_edge_line(line)) {
      if (!(first_node(line)  %e% nodes)) {
        nodes <- nodes | set(first_node(line))
      }
      if (!(second_node(line) %e% nodes)) {
        nodes <- nodes | set(second_node(line))
      }
    }
  }

  ## Create children dictionaries
  children <- hash()
  parents <- hash()
  for (n in nodes) {
    children_of_n <- set()
    for (line in lines) {
      if (is_edge_line(line)) {
        if ((first_node(line)  == n) && (!(second_node(line) %e% children_of_n))) {
          children_of_n <- children_of_n | set(second_node(line))
          parents_of_n <- parents_of_n | set(first_node(line))
        }
      }
    }
    children[[n]] <- children_of_n
    parents[[n]] <- parents_of_n
  }
  return(unlist(ancestors(TV, parents)))
}

get_ribbon_plot <- function(AIRHome) {
  dfr <- results_out_buffer()

  if (nrow(dfr) == 0) {
    return(tags$em("No results – model fitting failed. Check the log."))
  }

  nz_max <- function(x) if (all(is.na(x))) NA_real_ else max(x, na.rm = TRUE)
  nz_min <- function(x) if (all(is.na(x))) NA_real_ else min(x, na.rm = TRUE)

  if (nrow(dfr) == 0 || all(is.na(dfr$z1_ATE))) {
    return(tags$div("No causal estimates available for this run."))
  }

  # --- CI bounds ---
  lo_all     <- nz_min(c(dfr$z1_ATE_LCI, dfr$z2_ATE_LCI))
  hi_all     <- nz_max(c(dfr$z1_ATE_UCI, dfr$z2_ATE_UCI))
  lo_overlap <- nz_max(c(dfr$z1_ATE_LCI, dfr$z2_ATE_LCI))
  hi_overlap <- nz_min(c(dfr$z1_ATE_UCI, dfr$z2_ATE_UCI))
  overlap_exists <- hi_overlap >= lo_overlap && !any(is.na(c(lo_overlap, hi_overlap)))

  # --- Handle no overlap ---
  if (!overlap_exists) {
    sendSweetAlert(
      session = shiny::getDefaultReactiveDomain(),
      title = "Consensus could not be reached for this scenario",
      text = "The causal estimates derived from your two models (Z₁ and Z₂) do not agree within their confidence ranges.
This suggests that the data or model assumptions may not support a consistent causal interpretation for this configuration.",
      type = "error",
      btn_labels = "Continue",
      closeOnClickOutside = FALSE
    )
  }

  dfr0 <- dfr[1, ]

  p <- ggplot(dfr0, aes(x = Treatment)) +
    ## background
    geom_linerange(aes(ymin = -1.05, ymax = 1.05),
                   lwd = 6,
                   col = "black",
                   alpha = 1,
                   stat = "unique",
                   lineend = "round",
                   position = position_nudge(x = 0)) +

    ## +/- annotations
    annotate("text",
             x = 1, y = -1,
             label = "-", hjust = 2.5, vjust = 0.25,
             size = 5, color = "white") +
    annotate("text",
             x = 1, y = 1,
             label = "+", hjust = -0.35, vjust = 0.4,
             size = 5, color = "white") +

    ## red baseline ribbon (everything outside CIs)
    geom_linerange(aes(ymin = -1, ymax = 1),
                   lwd = 3.5,
                   col = "#C00000",
                   alpha = 1,
                   stat = "unique",
                   lineend = "round",
                   position = position_nudge(x = 0)) +

    ## yellow ribbons for Z1 and Z2 CIs
    geom_linerange(aes(ymin = z1_ATE_LCI, ymax = z1_ATE_UCI),
                   lwd = 3.5,
                   col = "#FCB514",
                   alpha = 1,
                   stat = "unique",
                   position = position_nudge(x = 0)) +
    geom_linerange(aes(ymin = z2_ATE_LCI, ymax = z2_ATE_UCI),
                   lwd = 3.5,
                   col = "#FCB514",
                   alpha = 1,
                   stat = "unique",
                   position = position_nudge(x = 0)) +

    ## conditional green overlap (only if overlap exists)
    {if (overlap_exists)
      geom_linerange(aes(ymin = lo_overlap, ymax = hi_overlap),
                     lwd = 3.5,
                     col = "#378855",
                     alpha = 1,
                     stat = "unique",
                     position = position_nudge(x = 0))
    } +

    ## Z1 range + mean
    geom_linerange(aes(ymin = z1_ATE_LCI, ymax = z1_ATE_UCI),
                   lwd = 3.5, col = "#9394A2",
                   alpha = 1, stat = "unique",
                   lineend = "round", position = position_nudge(x = 1)) +
    geom_point(aes(y = z1_ATE),
               col = "white", cex = 3, pch = 1, stroke = 1.25,
               position = position_nudge(x = 1)) +

    ## Z2 range + mean
    geom_linerange(aes(ymin = z2_ATE_LCI, ymax = z2_ATE_UCI),
                   lwd = 3.5, col = "#D4C7C7",
                   alpha = 1, stat = "unique",
                   lineend = "round", position = position_nudge(x = 0.5)) +
    geom_point(aes(y = z2_ATE),
               col = "white", cex = 3, pch = 1, stroke = 1.25,
               position = position_nudge(x = 0.5)) +

    ## horizontal zero line
    geom_segment(aes(x = 0.6, xend = 1.35, y = 0, yend = 0), lwd = 1.2) +

    coord_flip(clip = "off") +
    scale_x_discrete(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), limits = c(-1.1, 1.1)) +
    labs(y = "", x = "") +
    theme_void(base_size = 10) +
    theme(
      panel.background = element_rect(fill = "transparent", color = NA),
      plot.background  = element_rect(fill = "transparent", color = NA),
      aspect.ratio     = 0.2
    )
  if (rv$model_yn != "null") {
     p +
      ## arrow for algorithm flag
    geom_segment(data = dfr,
                 aes(x = 2.5, xend = 1.25, y = flag, yend = flag, color = algorithm),
                 arrow = arrow(length = unit(0.5, "cm")),
                 lwd = 1.2,
                 color = "#0F9ED5") +
    geom_point(data = dfr,
               aes(x = 2.5, y = flag, shape = algorithm),
               size = 3, color = "#0F9ED5")
  } else {
    p
  }

}


get_figure_caption <- function(AIRHome, df_vars) {
  caption <- paste0("Risk Difference: This chart represents the difference in outcomes resulting from a change in your experimental variable,",df_vars[1,][[2]],". The x-axis ranges from negative to positive effect, where the treatment, ", df_vars[2,][[2]]," either increases the likelihood of the outcome or decreases it, respectively. The midpoint corresponds to 'no significant effect.")
  return(caption)
}

get_ui_interpretation <- function(AIRHome, df_vars, Zvars, model_yn, n_obs) {
  dfr <- results_out_buffer()

  if (nrow(dfr) == 0) {
    return(tags$em("No results – model fitting failed. Check the log."))
  }

  # ── helpers ────────────────────────────────────────────────────────────────
  nz_max <- function(x) if (all(is.na(x))) NA_real_ else max(x, na.rm = TRUE)
  nz_min <- function(x) if (all(is.na(x))) NA_real_ else min(x, na.rm = TRUE)

  fmt_pct_signed <- function(x) {
    v <- round(x * 100, 0)
    if (v > 0) paste0("+", v, "%") else paste0(v, "%")
  }
  fmt_range <- function(lo, hi) paste0(fmt_pct_signed(lo), " to ", fmt_pct_signed(hi))
  fmt_range_abs <- function(lo, hi) {
    paste0(round(min(abs(lo), abs(hi)) * 100, 0), "\u2013",
           round(max(abs(lo), abs(hi)) * 100, 0), "%")
  }

  x_var <- df_vars[1, ][[2]]
  y_var <- df_vars[2, ][[2]]

  # ── CI overlap ─────────────────────────────────────────────────────────────
  lo_overlap <- nz_max(c(dfr$z1_ATE_LCI, dfr$z2_ATE_LCI))
  hi_overlap <- nz_min(c(dfr$z1_ATE_UCI, dfr$z2_ATE_UCI))
  overlap_exists <- !is.na(lo_overlap) && !is.na(hi_overlap) && hi_overlap >= lo_overlap

  # ── CI width note ──────────────────────────────────────────────────────────
  z1_width   <- dfr$z1_ATE_UCI[1] - dfr$z1_ATE_LCI[1]
  z2_width   <- dfr$z2_ATE_UCI[1] - dfr$z2_ATE_LCI[1]
  max_width  <- max(z1_width, z2_width, na.rm = TRUE)
  ci_width_note <- if (max_width < 0.10) {
    "The uncertainty range is narrow, so this is a fairly precise estimate \u2014 not just a directional guess."
  } else if (max_width < 0.25) {
    "The uncertainty range is moderate; the direction of the effect is reliable but the exact size is approximate."
  } else {
    paste0(
      "Note: the uncertainty range is wide (spanning more than ",
      round(max_width * 100, 0),
      " percentage points). The data supports only a rough directional conclusion here, not a precise effect size. ",
      "Collecting more data or refining the causal graph may narrow this range."
    )
  }

  # ── Z1/Z2 point estimate proximity ─────────────────────────────────────────
  z1_pt <- dfr$z1_ATE[1]
  z2_pt <- dfr$z2_ATE[1]
  pt_gap <- abs(z1_pt - z2_pt)
  proximity_note <- if (pt_gap < 0.05) {
    paste0(
      "Both independent analyses give nearly identical estimates ",
      "(Analysis A: ", fmt_range(dfr$z1_ATE_LCI[1], dfr$z1_ATE_UCI[1]),
      "; Analysis B: ", fmt_range(dfr$z2_ATE_LCI[1], dfr$z2_ATE_UCI[1]),
      "), which adds confidence to this finding."
    )
  } else if (pt_gap < 0.10) {
    paste0(
      "The two analyses give similar estimates ",
      "(Analysis A: ", fmt_range(dfr$z1_ATE_LCI[1], dfr$z1_ATE_UCI[1]),
      "; Analysis B: ", fmt_range(dfr$z2_ATE_LCI[1], dfr$z2_ATE_UCI[1]), ")."
    )
  } else {
    paste0(
      "The two analyses agree on direction but differ in magnitude ",
      "(Analysis A: ", fmt_range(dfr$z1_ATE_LCI[1], dfr$z1_ATE_UCI[1]),
      "; Analysis B: ", fmt_range(dfr$z2_ATE_LCI[1], dfr$z2_ATE_UCI[1]),
      "). Treat the exact size with some caution."
    )
  }

  # ── Effect direction (from consensus range, or full range if no overlap) ───
  if (overlap_exists) {
    flagdir <- if (lo_overlap > 0) "positive" else if (hi_overlap < 0) "negative" else "none"
  } else {
    all_bounds <- c(dfr$z1_ATE_LCI, dfr$z1_ATE_UCI, dfr$z2_ATE_LCI, dfr$z2_ATE_UCI)
    flagdir <- if (max(all_bounds, na.rm = TRUE) < 0) "negative" else if (min(all_bounds, na.rm = TRUE) > 0) "positive" else "none"
  }

  # ── Sample size caveat ─────────────────────────────────────────────────────
  sample_caveat <- if (!is.null(n_obs) && !is.na(n_obs) && n_obs < 500) {
    "Note: this dataset is relatively small, which means there may not have been enough data to detect a small-to-moderate effect even where the analyses agree. Treat directional conclusions with caution."
  } else {
    "These estimates are only as reliable as the variables included in your causal graph. If an important variable is missing, the true effect may differ."
  }

  # ══════════════════════════════════════════════════════════════════════════
  # NULL MODE — causal effect only, no model comparison
  # ══════════════════════════════════════════════════════════════════════════
  if (isTRUE(model_yn == "null")) {

    if (!overlap_exists) {
      return(paste0(
        "<strong>The two causal analyses gave conflicting answers.</strong><br><br>",
        "AIR estimates the effect of <em>", x_var, "</em> on <em>", y_var,
        "</em> using two different sets of background variables. ",
        "Here those two approaches disagree completely \u2014 ",
        "Analysis A estimates the effect at ", fmt_range(dfr$z1_ATE_LCI[1], dfr$z1_ATE_UCI[1]),
        ", while Analysis B estimates it at ", fmt_range(dfr$z2_ATE_LCI[1], dfr$z2_ATE_UCI[1]), ". ",
        "When the ranges don\u2019t even overlap, no single reliable conclusion can be drawn.<br><br>",
        "This often happens when the dataset is small, the causal graph includes uncertain relationships, ",
        "or a key variable is missing.<br><br>",
        "<strong>What to do:</strong> Review your causal graph for missing variables or incorrect edges. ",
        "Consider collecting more data before drawing conclusions."
      ))
    }

    direction_phrase <- if (flagdir == "positive") {
      paste0("is associated with a <strong>", fmt_range_abs(lo_overlap, hi_overlap),
             " increase</strong> in the probability of <em>", y_var, "</em>")
    } else if (flagdir == "negative") {
      paste0("is associated with a <strong>", fmt_range_abs(lo_overlap, hi_overlap),
             " decrease</strong> in the probability of <em>", y_var, "</em>")
    } else {
      paste0("does <strong>not meaningfully change</strong> the probability of <em>", y_var,
             "</em> \u2014 the effect is too small to distinguish from random variation")
    }

    consensus_text <- if (flagdir != "none") {
      paste0(
        "The range where both analyses agree \u2014 the most defensible single estimate \u2014 is ",
        "<strong>", fmt_range(lo_overlap, hi_overlap), "</strong>. "
      )
    } else {
      paste0(
        "The range where both analyses agree spans roughly ",
        "<strong>", fmt_range(lo_overlap, hi_overlap), "</strong>, ",
        "which AIR treats as no meaningful effect. "
      )
    }

    return(paste0(
      "<strong>Causal Effect Summary</strong><br><br>",
      "Based on your data and causal graph, changing <em>", x_var,
      "</em> from its baseline range to its experimental range ", direction_phrase, ".<br><br>",
      consensus_text, proximity_note, "<br><br>",
      ci_width_note, "<br><br>",
      "<strong>Caveat:</strong> ", sample_caveat
    ))
  }

  # ══════════════════════════════════════════════════════════════════════════
  # MODEL COMPARISON MODES: Yes / ATE / No
  # ══════════════════════════════════════════════════════════════════════════

  # ── State 1: No consensus ─────────────────────────────────────────────────
  if (!overlap_exists) {
    return(paste0(
      "<strong>The two causal analyses gave conflicting answers.</strong><br><br>",
      "AIR computes the effect of <em>", x_var, "</em> on <em>", y_var,
      "</em> using two different sets of background variables. ",
      "Here those two approaches disagree completely \u2014 ",
      "Analysis A estimates the effect at ", fmt_range(dfr$z1_ATE_LCI[1], dfr$z1_ATE_UCI[1]),
      ", while Analysis B estimates it at ", fmt_range(dfr$z2_ATE_LCI[1], dfr$z2_ATE_UCI[1]), ". ",
      "When the ranges don\u2019t even overlap, no single reliable conclusion can be drawn.<br><br>",
      "Wide uncertainty ranges often accompany this result and can indicate a small dataset, ",
      "uncertain relationships in the causal graph, or a missing variable driving the disagreement.<br><br>",
      "<strong>What to do:</strong> Review your causal graph for missing variables or incorrect edges. ",
      "Consider collecting more data before drawing conclusions."
    ))
  }

  # ── State 2: Internally inconsistent ──────────────────────────────────────
  if ((all(dfr$z1_ATE_UCI < dfr$z2_ATE) & all(dfr$z2_ATE_LCI > dfr$z1_ATE)) |
      (all(dfr$z1_ATE_LCI > dfr$z2_ATE) & all(dfr$z2_ATE_UCI < dfr$z1_ATE))) {
    dir_word <- if (flagdir == "positive") "a positive" else if (flagdir == "negative") "a negative" else "a possible"
    return(paste0(
      "<strong>The two analyses agree on direction but disagree on size.</strong><br><br>",
      "Both analyses find ", dir_word, " effect of <em>", x_var, "</em> on <em>", y_var,
      "</em>, but they differ meaningfully on how large that effect is. ",
      "Analysis A estimates <strong>", fmt_pct_signed(z1_pt),
      "</strong> and Analysis B estimates <strong>", fmt_pct_signed(z2_pt),
      "</strong> \u2014 a ", round(pt_gap * 100, 0), " percentage-point spread. ",
      "Their confidence ranges do overlap (roughly ", fmt_range(lo_overlap, hi_overlap),
      "), so a causal effect in that region is plausible, but the disagreement means ",
      "estimates should be treated as rough approximations rather than precise values.<br><br>",
      "This usually reflects limited data or structural uncertainty in the causal graph, ",
      "not necessarily a problem with your model."
    ))
  }

  # ── State 3: Standard interpretation ──────────────────────────────────────
  flag <- dfr$flag[1]

  inZ1 <- any(flag >= dfr$z1_ATE_LCI & flag <= dfr$z1_ATE_UCI, na.rm = TRUE)
  inZ2 <- any(flag >= dfr$z2_ATE_LCI & flag <= dfr$z2_ATE_UCI, na.rm = TRUE)
  bias_status <- if (inZ1 & inZ2) "unbiased" else if (inZ1 | inZ2) "partial" else "biased"

  # Degree of divergence from nearest CI boundary
  divergence <- if (!inZ1 & !inZ2) {
    if (flag > max(dfr$z1_ATE_UCI[1], dfr$z2_ATE_UCI[1])) {
      flag - max(dfr$z1_ATE_UCI[1], dfr$z2_ATE_UCI[1])
    } else {
      min(dfr$z1_ATE_LCI[1], dfr$z2_ATE_LCI[1]) - flag
    }
  } else if (!inZ1) {
    if (flag > dfr$z1_ATE_UCI[1]) flag - dfr$z1_ATE_UCI[1] else dfr$z1_ATE_LCI[1] - flag
  } else if (!inZ2) {
    if (flag > dfr$z2_ATE_UCI[1]) flag - dfr$z2_ATE_UCI[1] else dfr$z2_ATE_LCI[1] - flag
  } else {
    0
  }
  divergence_pp <- round(abs(divergence) * 100, 0)

  model_dir <- if (flag > hi_overlap) "overestimating" else if (flag < lo_overlap) "underestimating" else "accurate"

  # Bias source explanation with structural roles
  Z1_vars <- Zvars[Zvars$grp == "Z1", ]$Z
  Z2_vars <- Zvars[Zvars$grp == "Z2", ]$Z

  build_source <- function(vars, role) {
    if (length(vars) == 0 || all(is.na(vars)) || all(!nzchar(vars))) return(NULL)
    paste0("<em>", paste(vars, collapse = ", "), "</em> (", role, ")")
  }

  z1_role <- paste0(
    "direct background ", if (length(Z1_vars) == 1) "cause" else "causes",
    " of <em>", x_var, "</em>"
  )
  z2_role <- paste0(
    if (length(Z2_vars) == 1) "a variable" else "variables",
    " that influence the pathway between <em>", x_var, "</em> and <em>", y_var, "</em>"
  )

  bias_sources <- Filter(Negate(is.null), list(
    if (!inZ1) build_source(Z1_vars, z1_role) else NULL,
    if (!inZ2) build_source(Z2_vars, z2_role) else NULL
  ))

  bias_source_text <- if (bias_status != "unbiased" && length(bias_sources) > 0) {
    paste0(
      "<strong>Likely source of bias:</strong> ",
      paste(bias_sources, collapse = "; and "),
      " \u2014 highlighted in red on the causal graph. ",
      "Review these variables and the edges connected to them."
    )
  } else ""

  # Classifier spread (No mode only)
  classifier_text <- ""
  if (isTRUE(model_yn == "No")) {
    clf <- classifier_buffer()
    if (!is.null(clf) && nrow(clf) > 1) {
      clf_flags <- clf$flag[clf$algorithm != "ML Model Summary"]
      clf_flags <- clf_flags[!is.na(clf_flags)]
      if (length(clf_flags) > 1) {
        clf_spread <- max(clf_flags) - min(clf_flags)
        classifier_text <- if (clf_spread < 0.05) {
          paste0(
            "The five classifiers AIR trained internally gave consistent predictions ",
            "(range: ", fmt_range_abs(min(clf_flags), max(clf_flags)), "), ",
            "which reinforces this conclusion."
          )
        } else if (clf_spread < 0.15) {
          paste0(
            "The five classifiers AIR trained internally were moderately consistent ",
            "(range: ", fmt_range_abs(min(clf_flags), max(clf_flags)), ")."
          )
        } else {
          paste0(
            "Note: the five classifiers AIR trained internally gave a wide spread of predictions ",
            "(range: ", fmt_range_abs(min(clf_flags), max(clf_flags)), "). ",
            "This spread suggests the model may be sensitive to the choice of algorithm. ",
            "Interpret the specific effect size with caution."
          )
        }
      }
    }
  }

  # Headline
  headline <- if (bias_status == "unbiased") {
    "<strong>No bias detected.</strong>"
  } else if (bias_status == "partial") {
    "<strong>Possible bias detected.</strong>"
  } else {
    "<strong>Bias detected.</strong>"
  }

  # What the causal analysis found
  direction_phrase <- if (flagdir == "positive") {
    paste0("a <strong>", fmt_range_abs(lo_overlap, hi_overlap),
           " increase</strong> in the probability of <em>", y_var, "</em>")
  } else if (flagdir == "negative") {
    paste0("a <strong>", fmt_range_abs(lo_overlap, hi_overlap),
           " decrease</strong> in the probability of <em>", y_var, "</em>")
  } else {
    paste0("<strong>no meaningful change</strong> in the probability of <em>", y_var,
           "</em> \u2014 the effect is too small to distinguish from random variation")
  }

  found_text <- paste0(
    "<strong>What the causal analysis found:</strong> When <em>", x_var,
    "</em> moves from its baseline range to its experimental range, the expected outcome is ",
    direction_phrase, ". ", proximity_note
  )

  # Model comparison paragraph
  model_pct      <- fmt_pct_signed(flag)
  consensus_range <- fmt_range(lo_overlap, hi_overlap)

  model_para <- if (bias_status == "unbiased") {
    paste0(
      "Your model\u2019s predicted effect (", model_pct,
      ") falls within the range both analyses expect (", consensus_range,
      "). The model is performing as expected."
    )
  } else if (flagdir == "none") {
    paste0(
      "Your model predicts an effect of ", model_pct, " on <em>", y_var,
      "</em>, but the causal analysis finds no real effect ",
      "(both analyses agree the effect is near zero: ", consensus_range, "). ",
      "The relationship your model learned is likely spurious."
    )
  } else {
    over_under_word  <- if (model_dir == "overestimating") "overstating" else "understating"
    divergence_clause <- if (divergence_pp > 0) {
      paste0(" \u2014 your model is ", over_under_word, " the effect by roughly <strong>",
             divergence_pp, " percentage points</strong>")
    } else ""
    warning_clause <- if (bias_status == "partial") {
      " Your model\u2019s prediction falls within one analysis range but outside the other, which is a warning sign rather than a definitive verdict."
    } else {
      " Your model\u2019s prediction falls outside both analysis ranges."
    }
    paste0(
      "Your model predicts an effect of ", model_pct, " on <em>", y_var,
      "</em>, but the causal analysis puts the best estimate at ", consensus_range,
      divergence_clause, ".", warning_clause
    )
  }

  # Assemble all parts
  parts <- c(
    headline,
    "<br><br>",
    model_para,
    "<br><br>",
    found_text,
    if (nzchar(bias_source_text)) c("<br><br>", bias_source_text) else NULL,
    if (nzchar(classifier_text))  c("<br><br>", classifier_text)  else NULL,
    "<br><br>",
    ci_width_note,
    if (bias_status != "biased") c("<br><br>", "<strong>Caveat:</strong> ", sample_caveat) else NULL
  )

  paste0(parts, collapse = "")
}


# ──────────────────────────────────────────────────────────────
#  Histogram for the experimental (X) variable
# ──────────────────────────────────────────────────────────────
get_histogram_x <- function(df, xvar, dir, threshold) {

  v       <- df[[xvar]]
  if (isTRUE(dir)){ treated <- ifelse(v >= threshold, 1, 0)} else { treated <- ifelse(v <= threshold, 1, 0)}

  plot_df <- data.frame(
    x         = v,
    Treatment = factor(ifelse(treated, "Treated", "Untreated"),
                       levels = c("Untreated", "Treated"))   # lock legend order
  )

  colors <- c(Treated   = "#5C9AFF",
              Untreated = "#EAE1D7")

  untreated_label <- paste0("Baseline") # (< ", signif(threshold, 4), ")")
  treated_label   <- paste0("Experimental") #  (≥ ",  signif(threshold, 4), ")")

  bw <- diff(range(v, na.rm = TRUE)) / 30
  if (bw <= 0 || is.na(bw) || is.infinite(bw)) bw <- 1        # <- arrow fix

  ggplot(plot_df, aes(x = x, fill = Treatment, colour = Treatment)) +
    geom_rug(sides = "b") +
    geom_histogram(binwidth = bw, colour = "black") +
    geom_vline(xintercept = threshold,
               colour = "grey20", linetype = "dashed", linewidth = 1) +
    scale_fill_manual(values = colors,
                      breaks = c("Untreated", "Treated"),
                      labels = c(untreated_label, treated_label)) +
    scale_colour_manual(values = colors,
                        breaks = c("Untreated", "Treated"),
                        labels = NULL) +
    labs(title = paste("Distribution of", xvar),
         x = NULL, fill = NULL) +
    guides(colour = "none") +
    theme_minimal(base_size = 10) +
    theme(text             = element_text(colour = "#666666", face = "bold"),
          panel.background = element_rect(fill = "transparent", colour = NA),
          plot.background  = element_rect(fill = "transparent", colour = NA),
          legend.position  = "bottom",
          legend.direction = "horizontal")
}

# ──────────────────────────────────────────────────────────────
#  Histogram for the outcome (Y) variable
# ──────────────────────────────────────────────────────────────
get_histogram_y <- function(df, yvar, dir, threshold) {

  v       <- df[[yvar]]
  if (isTRUE(dir)){ success <- ifelse(v >= threshold, 1, 0)} else { success <- ifelse(v <= threshold, 1, 0)}

  plot_df <- data.frame(
    x      = v,
    Result = factor(ifelse(success, "Success", "Fail"),
                    levels = c("Fail", "Success"))
  )

  colors <- c(Success = "#5C6CFF",
              Fail    = "#EAE1D7")

  fail_label    <- paste0("Fail") # (< ",  signif(threshold, 4), ")")
  success_label <- paste0("Success") # (≥ ", signif(threshold, 4), ")")

  bw <- diff(range(v, na.rm = TRUE)) / 30
  if (bw <= 0 || is.na(bw) || is.infinite(bw)) bw <- 1

  ggplot(plot_df, aes(x = x, fill = Result, colour = Result)) +
    geom_rug(sides = "b") +
    geom_histogram(binwidth = bw, colour = "black") +
    geom_vline(xintercept = threshold,
               colour = "grey20", linetype = "dashed", linewidth = 1) +
    scale_fill_manual(values = colors,
                      breaks = c("Fail", "Success"),
                      labels = c(fail_label, success_label)) +
    scale_colour_manual(values = colors,
                        breaks = c("Fail", "Success"),
                        labels = NULL) +
    labs(title = paste("Distribution of", yvar),
         x = NULL, fill = NULL) +
    guides(colour = "none") +
    theme_minimal(base_size = 10) +
    theme(text             = element_text(colour = "#666666", face = "bold"),
          panel.background = element_rect(fill = "transparent", colour = NA),
          plot.background  = element_rect(fill = "transparent", colour = NA),
          legend.position  = "bottom",
          legend.direction = "horizontal")
}

get_updated_graph <- function(dot, xvar, yvar, Zvars = NULL) {
  lines <- strsplit(dot, "\n", fixed = TRUE)[[1]]

  make_decl <- function(node, colour) {
    sprintf('  "%s" [style=filled, fillcolor="%s"];', node, colour)
  }

  # ---- build lists safely ----
  Z1 <- if (!is.null(Zvars)) Zvars$Z[Zvars$grp == "Z1"] else character(0)
  Z2 <- if (!is.null(Zvars)) Zvars$Z[Zvars$grp == "Z2"] else character(0)

  nodes_by_color <- list(
   # "#FFC107" = c(xvar, yvar),
    "#5C9AFF" = xvar,
    "#5C6CFF" = yvar,
    "#9394A2" = Z1,
    "#D4C7C7" = Z2
  )

  # flatten + clean
  nodes_vec <- unique(na.omit(unlist(nodes_by_color)))
  nodes_vec <- nodes_vec[nzchar(nodes_vec)]
  if (!length(nodes_vec)) return(dot)

  # escape for regex
  esc <- function(s) gsub("([][{}()+*^$|?.\\-])", "\\\\\\1", s)
  pat <- sprintf('^\\s*"(%s)"\\s*\\[', paste(esc(nodes_vec), collapse = "|"))

  # drop existing declarations for these nodes
  keep <- !grepl(pat, lines, perl = TRUE)
  lines <- lines[keep]

  # synthesize new declarations
  new_lines <- c(
    if (nzchar(xvar)) sapply(xvar, make_decl, colour = "#5C9AFF"),
    if (nzchar(yvar)) sapply(yvar, make_decl, colour = "#5C6CFF"),
    if (length(Z1))   sapply(Z1,   make_decl, colour = "#9394A2"),
    if (length(Z2))   sapply(Z2,   make_decl, colour = "#D4C7C7")
  )

  # insert after the digraph line (more robust than hard-coded 2L)
  insert_at <- which(grepl("^\\s*digraph\\b", lines))[1]
  if (is.na(insert_at)) insert_at <- 0L
  lines <- append(lines, new_lines, after = insert_at)

  dot <- paste(lines, collapse = "\n")
  dot <- sub("(digraph[^\\{]*\\{)", "\\1\n  graph [tooltip=\"Updated Causal graph\"]", 
           dot, perl = TRUE)
  return(dot)
}

get_final_graph <- function(dot, xvar, yvar, Zvars, dfr = results_out_buffer()) {
  # Highlight treatment and outcome nodes
  dot <- change_node_color(dot, xvar, "'#5C9AFF'")
  dot <- change_node_color(dot, yvar, "'#5C6CFF'")

  # Extract Z1 and Z2 sets
  Z1 <- Zvars[Zvars$grp == "Z1", ]$Z
  Z2 <- Zvars[Zvars$grp == "Z2", ]$Z
  log_debug("get final graph vars")
  log_debug("dfr$flag: {msg}", msg = dfr$flag)
  log_debug("dfr$z1_ATE_UCI: {msg}", msg = dfr$z1_ATE_UCI)
  log_debug("dfr$z2_ATE_UCI: {msg}", msg = dfr$z2_ATE_UCI)
  log_debug("dfr$z1_ATE_LCI: {msg}", msg = dfr$z1_ATE_LCI)
  log_debug("dfr$z2_ATE_LCI: {msg}", msg = dfr$z2_ATE_LCI)
  log_debug("Z1: {msg}", msg = Z1)
  log_debug("Z2: {msg}", msg = Z2)

  # Ribbon plot flag logic to color additional nodes
# Determine whether flag is inside each interval
in_z1 <- with(dfr, flag >= z1_ATE_LCI & flag <= z1_ATE_UCI)
in_z2 <- with(dfr, flag >= z2_ATE_LCI & flag <= z2_ATE_UCI)

# Use any() to handle multiple rows safely
in_z1 <- any(in_z1, na.rm = TRUE)
in_z2 <- any(in_z2, na.rm = TRUE)

if (in_z1 && in_z2) {
  # Inside both → do nothing
  log_debug("flag within both intervals")

} else if (xor(in_z1, in_z2)) {
  # Inside one, outside the other → paint the one that's outside
  if (!in_z1) {
    dot <- paint_nodes(dot, Z1, "'#FFC107'")
    log_debug("flag outside Z1 only")
  }
  if (!in_z2) {
    dot <- paint_nodes(dot, Z2, "'#FFC107'")
    log_debug("flag outside Z2 only")
  }

} else if (!in_z1 && !in_z2) {
  # Outside both → paint both in red
  dot <- paint_nodes(dot, Z1, "'#C00000'")
  dot <- paint_nodes(dot, Z2, "'#C00000'")
  log_debug("flag outside both intervals")
}
  dot <- sub("(digraph[^\\{]*\\{)", "\\1\n  graph [tooltip=\"Final Causal graph\"]", 
           dot, perl = TRUE)
  return(dot)
}

graphviz_to_png_data <- function(dot_code) {
  svg_str <- DiagrammeRsvg::export_svg(grViz(dot_code))
  svg_obj <- xml2::read_xml(svg_str)

  # Write PNG to raw vector
  raw_conn <- rawConnection(raw(0), open = "wb")
  on.exit(close(raw_conn), add = TRUE)

  rsvg::rsvg_png(svg_obj, file = raw_conn)

  rawConnectionValue(raw_conn)
}
# --- END INLINE: scripts/AIR_functions.R ---

disconnectMessage2()
actionButton("disconnect", "Disconnect the app")

if (!dir.exists(paste0(AIRHome, "/data/"))) {
  dir.create(paste0(AIRHome, "/data/"), recursive = TRUE)
}

if (!dir.exists(paste0(AIRHome, "/input/"))) {
  dir.create(paste0(AIRHome, "/input/"), recursive = TRUE)
}

if (Sys.info()["sysname"] == "Linux") {
  # --- BEGIN INLINE: scripts/tetrad_utils.R ---
# ----- Utility Functions -----

# ----- Graph Visualization -----
visualize_graph <- function(graph) {
  if (requireNamespace("rstudioapi", quietly = TRUE) && rstudioapi::isAvailable()) {
    b <- TRUE
  } else if (Sys.getenv("RSTUDIO") == "1") {
    b <- TRUE
  } else {
    b <- FALSE
  }
  
  if (b) {
    if (!is.null(graph)) {
      dot <- graph_to_dot(graph)
      dot <- sub("(digraph[^\\{]*\\{)", "\\1\n  graph [tooltip=\"Causal graph\"]", dot, perl = TRUE)
      grViz(dot)
    } else {
      log_debug("No graph generated. Please check the BOSS execution.")
    }
  }
}
  # --- END INLINE: scripts/tetrad_utils.R ---
# --- BEGIN INLINE: scripts/TetradSearch.R ---
# This class translates some select methods from TetradSearch.py in py-tetrad
# for use in R using rJava.
#
# This is a temporary class, as a much better effort at translating these
# methods is underway by another group.

TetradSearch <- setRefClass(
  "TetradSearch",
  
  fields = list(
    data = "data.frame",          # Input dataset
    sample_size = "numeric",      # Sample size
    data_model = "ANY",           # Data Model (Tabular data or Covariance Matrix)
    score = "ANY",                # Score object
    test = "ANY",                 # IndependenceTest object
    mc_test = "ANY",              # IndependenceTest for the Markov Checker
    mc_ind_results = "ANY",       # Markov Checker independence test results
    knowledge = "ANY",            # Background knowledge object
    graph = "ANY",                # Resulting graph
    search = "ANY",               # Search object
    params = "ANY"                # Parameters object
  ),
  
  methods = list(
    
    # Initialize the TetradSearch object
    #
    # @param data A data frame containing the dataset to be analyzed.
    # @return A TetradSearch object.
    initialize = function(data) {
      log_debug("Initializing TetradSearch object...")
      
      if (!is.data.frame(data)) {
        stop("Data must be a data.frame")
      }
      
      .self$data <- data
      .self$sample_size <- nrow(data)
      log_debug("Data frame dimensions: {msg}", msg = dim(data))
      log_debug("Sample size set to: {msg}", msg = .self$sample_size)
      
      .self$data_model <- .self$data_frame_to_tetrad_dataset(data)
      .self$data_model <- .jcast(.self$data_model, "edu.cmu.tetrad.data.DataModel")
      
      log_debug("Tetrad DataSet created.")
      
      .self$params <- .jnew("edu.cmu.tetrad.util.Parameters")
      
      .self$knowledge <- .jnew("edu/cmu/tetrad/data/Knowledge")
      log_debug("Knowledge instance created.")
      log_debug("TetradSearch object initialized successfully.")
    },
    
    # Make sure the score object is initialized
    .check_score = function() {
      if (is.null(.self$score)) {
        stop("Error: The 'score' field has not been initialized yet. Please \
                 set a score before running the algorithm.")
      }
    },
    
    .setParam = function(key, value) {
      safe_jcall(.self$params, "V", "set", key, .jcast(.jnew("java/lang/Boolean", value), "java/lang/Object"))
    },
    
    .setParamInt = function(key, value) {
      safe_jcall(.self$params, "V", "set", key, .jcast(.jnew("java/lang/Integer", as.integer(value)), "java/lang/Object"))
    },
    
    .set_knowledge = function() {
      safe_jcall(.self$search, "V", "setKnowledge", .jcast(.self$knowledge, "edu.cmu.tetrad.data.Knowledge"))
    },
    
    # Run the search algorithm, for the typical case
    .run_search = function() {
      .self$.set_knowledge()
      .self$graph <- .jcast(.self$search$search(), "edu.cmu.tetrad.graph.Graph")
    },
    
    # Make sure the test object is initialized
    .check_test = function() {
      if (is.null(.self$test)) {
        stop("Error: The 'test' field has not been initialized yet. Please \
                 set a test before running the algorithm.")
      }
    },
    
    # Add a variable to a specific tier in the knowledge
    #
    # @param tier The tier to which the variable should be added.
    # @param var_name The name of the variable to add.
    add_to_tier = function(tier, var_name) {
      log_debug("Adding variable {var_name} to tier {tier}...", var_name = var_name, tier = tier)
      tryCatch({
        tier <- as.integer(tier)
        var_name <- as.character(var_name)
        safe_jcall(.self$knowledge, "V", "addToTier", tier, var_name)
        log_debug("Variable {var_name} added to tier {tier}", var_name = var_name, tier = tier)
      }, error = function(e) {
        log_debug("Error adding variable to tier: {msg}", msg = e$message)
      })
    },
    
    # Set the verbose flag
    #
    # @param verbose TRUE or FALSE
    set_verbose = function(verbose) {
      .self$.setParam("verbose", verbose)
    },
    
    # Set the score to the SEM BIC.
    #
    # @param penalty_discount The penalty discount to use in the SemBicScore calculation.
    use_sem_bic = function(penalty_discount = 2) {
      .self$.setParamDouble("penaltyDiscount", penalty_discount)
      .self$score <- .jnew("edu.cmu.tetrad.algcomparison.score.SemBicScore")
      .self$score <- .jcast(.self$score, "edu.cmu.tetrad.algcomparison.score.ScoreWrapper")
      log_debug("SemBicScore object created with penalty discount set.")
    },
    
    # Set the test to Fisher Z
    # # @param alpha The significance cutoff.
    use_fisher_z = function(alpha = 0.01, use_for_mc = FALSE) {
      .self$.setParamDouble("alpha", alpha)
      
      if (use_for_mc) {
        .self$mc_test <- .jnew("edu.cmu.tetrad.algcomparison.independence.FisherZ")
        .self$mc_test <- .jcast(.self$mc_test, "edu.cmu.tetrad.algcomparison.independence.IndependenceWrapper")
      } else {
        .self$test <- .jnew("edu.cmu.tetrad.algcomparison.independence.FisherZ")
        .self$test <- .jcast(.self$test, "edu.cmu.tetrad.algcomparison.independence.IndependenceWrapper")
      }
      
      log_debug("Fisher Z object created with alpha set.")
    },
    
    # --- Internal parameter helpers ---
    
    .setParamDouble = function(key, value) {
      safe_jcall(.self$params, "V", "set", key, .jcast(.jnew("java/lang/Double", as.double(value)), "java/lang/Object"))
    },
    
    # Run the BOSS algorithm
    #
    # @param num_starts The number of random restarts to do; the model with the best BIC score overall is returned.
    # @param use_bes TRUE if the algorithm should finish up with a call to BES (Backward Equivalence Search from
    #   the FGES algorithm) to guarantee correctness under Faithfulness.
    # @param time_lag Default 0; if > 1, a time lag model of this order is constructed.
    # @param use_data_order TRUE if the original data order should be used for the initial permutation. If
    #   num_starts > 1, random permuatations are used for subsequent restarts.
    # @param output_cpdag TRUE if a CPDAG should be output, FALSE if a DAG should be output.
    # @return The estimated graph.
    run_boss = function(num_starts = 30, use_bes = FALSE, time_lag = 0, use_data_order = TRUE, output_cpdag = TRUE) {
      log_debug("Running BOSS algorithm...")
      
      .self$.setParam("useBes", use_bes)
      .self$.setParamInt("numStarts", num_starts)
      .self$.setParamInt("timeLag", time_lag)
      .self$.setParam("useDataOrder", use_data_order)
      .self$.setParam("outputCpdag", output_cpdag)
      
      dataModel <- .jcast(.self$data_model, "edu.cmu.tetrad.data.DataModel")
      
      boss <- .jnew("edu.cmu.tetrad.algcomparison.algorithm.oracle.cpdag.Boss", .self$score)
      safe_jcall(boss, "V", "setKnowledge", .self$knowledge)
      
      graph <<- safe_jcall(boss, "Ledu/cmu/tetrad/graph/Graph;", "search", dataModel, .self$params)
      .self$graph <- graph
      
      log_debug("BOSS search completed.")
    },

    get_java = function() {
      return(.self$graph)
    },
    
    # This method prints the structure of the graph estimated by the most recent algorithm call.
    print_graph = function() {
      log_debug("Attempting to print the graph...")
      if (is.null(.self$graph)) {
        log_debug("No graph generated yet. Please run an algorithm first.")
      } else {
        log_debug("Graph structure:")
        log_debug("{msg}", msg = .self$graph)
      }
      invisible(.self$graph)
    },
    
    # An adjustment set for a pair of nodes <source, target> for a CPDAG is a set of nodes that blocks
    # all paths from the source to the target that cannot contribute to a calculation for the total effect
    # of the source on the target in any DAG in a CPDAG while not blocking any path from the source to the target
    # that could be causal. In typical causal graphs, multiple adjustment sets may exist for a given pair of
    # nodes. This method returns up to maxNumSets adjustment sets for the pair of nodes <source, target>
    # fitting a certain description.
    #
    # The description is as follows. We look for adjustment sets of variables that are close to either the
    # source or the target (or either) in the graph. We take all possibly causal paths from the source to the
    # target into account but only consider other paths up to a certain specified length. (This maximum length
    # can be unlimited for small graphs.)
    #
    # Within this description, we list adjustment sets in order or increasing size. Hopefully, these parameters
    # along with the size ordering can help to give guidance for the user to choose the best adjustment set for
    # their purposes when multiple adjustment sets are possible.
    #
    # @param source                  The source node whose sets will be used for adjustment.
    # @param target                  The target node whose sets will be adjusted to match the source node.
    # @param maxNumSets              The maximum number of sets to be adjusted. If this value is less than or equal to
    #                                0, all sets in the target node will be adjusted to match the source node.
    # @param maxDistanceFromEndpoint The maximum distance from the endpoint of the trek to consider for adjustment.
    # @param nearWhichEndpoint       The endpoint(s) to consider for adjustment; 1 = near the source, 2 = near the
    #                                target, 3 = near either.
    # @param maxPathLength           The maximum length of the path to consider for backdoor paths. If a value of -1 is
    #                                given, all paths will be considered.
    # @return A list of adjustment sets for the pair of nodes &lt;source, target&gt;. Return an smpty
    # list if source == target or there is no amenable path from source to target.
    get_adjustment_sets = function(graph, source, target, max_num_sets = 10, max_distance_from_point = 5,
                                   near_which_endpoint = 1, max_path_length = 20) {
      log_debug("Getting adjustment sets for: {source} -> {target}", source = source, target = target)
      
      # Look up Node objects by name
      source_node <- safe_jcall(graph, "Ledu/cmu/tetrad/graph/Node;", "getNode", source)
      target_node <- safe_jcall(graph, "Ledu/cmu/tetrad/graph/Node;", "getNode", target)
      
      if (is.jnull(source_node)) stop(paste("Source node", source, "not found in the graph."))
      if (is.jnull(target_node)) stop(paste("Target node", target, "not found in the graph."))
      
      # Get Paths object from Graph
      paths <- safe_jcall(graph, "Ledu/cmu/tetrad/graph/Paths;", "paths")
      
      # Java List<Set<Node>>
      sets_list <- safe_jcall(paths,
                          "Ljava/util/List;",
                          "adjustmentSets",
                          source_node,
                          target_node,
                          "dag",
                          as.integer(max_num_sets),
                          as.integer(max_distance_from_point),
                          as.integer(near_which_endpoint),
                          as.integer(max_path_length))
      
      
      size <- safe_jcall(sets_list, "I", "size")
      log_debug("Number of adjustment sets: {size}", size = size)
      
      # Convert Java List<Set<Node>> to R list of character vectors
      size <- safe_jcall(sets_list, "I", "size")
      result <- vector("list", size)
      
      for (i in seq_len(size)) {
        jset <- safe_jcall(sets_list, "Ljava/lang/Object;", "get", as.integer(i - 1))
        jarray <- safe_jcall(jset, "[Ljava/lang/Object;", "toArray")
        result[[i]] <- sapply(jarray, function(n) safe_jcall(n, "S", "getName"))
      }
      
      return(result)
    },
    
    print_adjustment_sets = function(adjustment_sets) {
      if (length(adjustment_sets) == 0) {
        log_debug("No adjustment sets found.")
        return()
      }
      
      for (i in seq_along(adjustment_sets)) {
        set <- adjustment_sets[[i]]
      
        # headline
        log_debug("Adjustment set {i}:", i = i)
      
        if (length(set) == 0) {
          log_debug("Adjustment set is empty.")
        } else {
          # collapse first, then log
          log_debug(paste(set, collapse = ", "))
        }
      }
    },
    
    # Performs a Markov check on a graph with respect to the supplied dataset and returns statistics
    # showing performance on that check.
    #
    # @param graph The graph to perform the Markov check on. This may be a DAG, CPDAG, MAG or PAG.
    # @param percent_resample Tests are done using random subsamples of the data per test, if this is
    #   less than 1, or all of the data, if it is equal to 1.
    # @param condition_set_type The type of conditioning set to use for the Markov check, one of:
    #   GLOBAL_MARKOV, LOCAL_MARKOV, PARENTS_AND_NEIGHBORS, MARKOV_BLANKET, RECURSIVE_MSEP, NONCOLLIDERS_ONLY,
    #   ORDERED_LOCAL_MARKOV, or ORDERED_LOCAL_MARKOV_MAG
    # @param find_smallest_subset Whether to find the smallest subset for a given set that yields independence.
    # @param parallelized TRUE if conditional independencies should be checked in parallel.
    # @effective_sample_size The effective sample size to use for calculations, or -1 if the actual sample size.
    # @return Marov checker statistics as a named list.
    markov_check = function(graph, percent_resample = 1, condition_set_type = "ORDERED_LOCAL_MARKOV",
                            find_smallest_subset = FALSE, parallelized = TRUE, effective_sample_size = -1) {
      log_debug("Running Markov check...")
      
      if (is.null(.self$mc_test)) {
        stop("A test for the Markov Checker has not been set. Please call a `use_*` method with `use_for_mc = TRUE`.")
      }
      
      condition_set_type_ <- .jfield("edu.cmu.tetrad.search.ConditioningSetType",
                                     name = condition_set_type,
                                     sig = "Ledu/cmu/tetrad/search/ConditioningSetType;")
      
      dataModel <- .jcast(.self$data_model, "edu.cmu.tetrad.data.DataModel")
      
      test_ <- safe_jcall(.self$mc_test, "Ledu/cmu/tetrad/search/test/IndependenceTest;",
                      "getTest", dataModel, .self$params)
      
      mc <- .jnew("edu.cmu.tetrad.search.MarkovCheck", graph, test_, condition_set_type_)
      
      # Configure it
      safe_jcall(mc, "V", "setFractionResample", as.double(percent_resample))
      safe_jcall(mc, "V", "setFindSmallestSubset", find_smallest_subset)
      safe_jcall(mc, "V", "setParallelized", parallelized)
      
      # Generate results
      safe_jcall(mc, "V", "generateAllResults")
      .self$mc_ind_results <- safe_jcall(mc, "Ljava/util/List;", "getResults", TRUE)
      
      # Set sample size if specified
      ### temporary fix ((mdk)) ??
      ### originally: if (sample_size != -1) {
      if(effective_sample_size != -1) {
        ### temporary fix ((mdk)) ??   
        ### originally: safe_jcall(mc, "V", "setSampleSize", as.integer(sample_size))
        safe_jcall(mc, "V", "setSampleSize", as.integer(effective_sample_size))
      }
      
      # Extract statistics
      ad_ind <- safe_jcall(mc, "D", "getAndersonDarlingP", TRUE)
      ad_dep <- safe_jcall(mc, "D", "getAndersonDarlingP", FALSE)
      ks_ind <- safe_jcall(mc, "D", "getKsPValue", TRUE)
      ks_dep <- safe_jcall(mc, "D", "getKsPValue", FALSE)
      bin_indep <- safe_jcall(mc, "D", "getBinomialPValue_", TRUE)
      bin_dep <- safe_jcall(mc, "D", "getBinomialPValue_", FALSE)
      frac_dep_ind <- safe_jcall(mc, "D", "getFractionDependent", TRUE)
      frac_dep_dep <- safe_jcall(mc, "D", "getFractionDependent", FALSE)
      num_tests_ind <- safe_jcall(mc, "I", "getNumTests", TRUE)
      num_tests_dep <- safe_jcall(mc, "I", "getNumTests", FALSE)
      
      # Return as a named list
      return(list(
        ad_ind = ad_ind,
        ad_dep = ad_dep,
        ks_ind = ks_ind,
        ks_dep = ks_dep,
        bin_indep = bin_indep,
        bin_dep = bin_dep,
        frac_dep_ind = frac_dep_ind,
        frac_dep_dep = frac_dep_dep,
        num_tests_ind = num_tests_ind,
        num_tests_dep = num_tests_dep,
        mc = mc
      ))
    },

    # Converts the given R data frame to a (possibly mixed) Tetrad DataSet.
    #
    # @param df The R data frame to translate. Continuous columns should be of type 'numeric' and the
    #   discrete columns of type 'integer'.
    data_frame_to_tetrad_dataset = function(df) {
      stopifnot(require(rJava))
      
      nrows <- nrow(df)
      ncols <- ncol(df)
      
      # Create Java ArrayList<Node>
      var_list <- .jnew("java/util/ArrayList")
      
      # Prepare empty double[][] and int[][] (as Java arrays)
      cont_data <- vector("list", ncols)
      disc_data <- vector("list", ncols)
      
      for (j in seq_len(ncols)) {
        name <- colnames(df)[j]
        col <- df[[j]]
        
        if (is.numeric(col)) {
          variable <- .jnew("edu/cmu/tetrad/data/ContinuousVariable", name)
          node <- .jcast(variable, "edu/cmu/tetrad/graph/Node")
          safe_jcall(var_list, "Z", "add", .jcast(node, "java/lang/Object"))
          cont_data[[j]] <- .jarray(as.numeric(col), dispatch = TRUE)
          disc_data[[j]] <- .jnull("[I")  # null int[] for discrete
        } else if (is.integer(col) || is.factor(col)) {
          num_categories <- length(unique(na.omit(col)))
          variable <- .jnew("edu/cmu/tetrad/data/DiscreteVariable", name, as.integer(num_categories))
          node <- .jcast(variable, "edu/cmu/tetrad/graph/Node")
          safe_jcall(var_list, "Z", "add", .jcast(node, "java/lang/Object"))
          cont_data[[j]] <- .jnull("[D")  # null double[] for continuous
          disc_data[[j]] <- .jarray(as.integer(col), dispatch = TRUE)
        } else {
          stop(paste("Unsupported column type:", name))
        }
      }
      
      # Convert R lists of arrays to Java double[][] and int[][]
      j_cont_data <- .jarray(cont_data, dispatch = TRUE)
      j_disc_data <- .jarray(disc_data, dispatch = TRUE)
      
      # Call static Java helper method
      ds <- safe_jcall("edu.cmu.tetrad.util.DataSetHelper",
                   "Ledu/cmu/tetrad/data/DataSet;",
                   "fromR",
                   .jcast(var_list, "java.util.List"),
                   as.integer(nrows),
                   .jcast(j_cont_data, "[[D"),
                   .jcast(j_disc_data, "[[I"))
      
      return(ds)
    }
  )
)
# --- END INLINE: scripts/TetradSearch.R ---
  
  if (!dir.exists(paste0(AIRHome, "/data/"))) {
    dir.create(paste0(AIRHome, "/data/"), recursive = TRUE)
  }
  
  if (!dir.exists(paste0(AIRHome, "/input/"))) {
    dir.create(paste0(AIRHome, "/input/"), recursive = TRUE)
  }
} 

```

```{r shinyjs-init, context="ui", echo=FALSE}
useShinyjs()
```
# Analysis of Bias

##  {.sidebar}

```{r sidebar}
br()
h5("Step 1- Upload data file:")
uiOutput("ui_file1")
uiOutput('ui_file2')
uiOutput('ui_buildButton')
uiOutput('step3')
uiOutput('xvar')
uiOutput('ui_threshold_x')
uiOutput('yvar')
uiOutput('ui_threshold_y')
uiOutput('step4')
uiOutput('ui_model_exist')
uiOutput('ui_model_upload')
uiOutput('ui_ate_upload')
uiOutput('ui_goButton')
br()
uiOutput('ui_dl_btn')

br()

tagList(
  # 1) Normal Shiny button lives in the sidebar
  uiOutput("ui_reset_btn"),

  # 2–3) JS: switch between fixed-at-bottom (no scroll) and sticky (when scroll shows up)
  tags$script(HTML("
    (function () {
      const slot = document.getElementById('reload_slot');
      if (!slot) return;

      // placeholder preserves original position for 'sticky' mode
      const ph = document.createElement('span');
      ph.id = 'reload_slot_ph';
      slot.parentNode.insertBefore(ph, slot);

      function getSidebar() {
        // Quarto dashboard sidebar container
        return document.querySelector('.quarto-dashboard-sidebar')
            || document.querySelector('.sidebar')
            || document.querySelector('[class*=sidebar]');
      }

      function getScrollEl(sidebar) {
        // find the actual scrolling element inside the sidebar
        const nodes = [sidebar, ...sidebar.querySelectorAll('*')];
        for (const n of nodes) {
          const cs = getComputedStyle(n);
          if ((cs.overflowY === 'auto' || cs.overflowY === 'scroll') &&
              n.scrollHeight > n.clientHeight + 1) return n;
        }
        return sidebar;
      }

      function placeFixed() {
        const sb = getSidebar();
        if (!sb) return;
        const rect = sb.getBoundingClientRect();
        // move slot to body so overflow clipping won't affect it
        if (slot.parentNode !== document.body) document.body.appendChild(slot);
        Object.assign(slot.style, {
          position: 'fixed',
          bottom: '12px',
          left: (rect.left + 12) + 'px',
          width: (rect.width - 24) + 'px',
          zIndex: 10000,
          background: 'transparent'
        });
      }

      function placeSticky() {
        // move slot back under the placeholder, full-width of sidebar content
        if (slot.parentNode !== ph.parentNode) {
          ph.parentNode.insertBefore(slot, ph.nextSibling);
        }
        Object.assign(slot.style, {
          position: 'sticky',
          bottom: '12px',
          left: '',
          width: '100%',
          zIndex: 1,
          background: 'transparent'
        });
      }

      function refreshPlacement() {
        const sb = getSidebar();
        if (!sb) return;
        const scroller = getScrollEl(sb);
        const needsScroll = scroller.scrollHeight > scroller.clientHeight + 1;
        if (needsScroll) placeSticky(); else placeFixed();
      }

      // run now and on likely changes
      setTimeout(refreshPlacement, 0);
      window.addEventListener('resize', refreshPlacement);
      document.addEventListener('shiny:value', refreshPlacement);
      document.addEventListener('shiny:outputinvalidated', refreshPlacement);

      const sb = getSidebar();
      if (sb) {
        new MutationObserver(refreshPlacement)
          .observe(sb, { childList: true, subtree: true });
      }
    })();
  "))
)

```

## Column {width = "40%"}

```{r graph-pane}
# grVizOutput('blankGraph')
uiOutput('ui_graph_pane')
```

## Column {width = "60%"}

```{r right-of-graph}
# plotOutput('histogram_x', height = "50%")
uiOutput('second_column_content')
# uiOutput('ui_top_right_pane')
```
# Docs

```{r info-tab, echo=FALSE, message=FALSE, warning=FALSE, results='asis'}
# Read the file and output its content as-is
knitr::opts_knit$set(root.dir = AIRHome)
knitr::opts_chunk$set(progress = FALSE)
options(knitr.progress = FALSE)

cat(knitr::knit_child("info.qmd", quiet = TRUE, envir = knitr::knit_global()))
```

# Key

-   Graph Colors
    -   [Blue]{style="color: #5C9AFF; font-weight: bold;"}: X Variable of interest, selected by the user in step 2. 
    -   [Violet]{style="color: #5C6CFF; font-weight: bold;"}: Y Variable of interest, selected by the user in step 2.
    -   [Dark Gray]{style="color: #9394A2; font-weight: bold;"}: Any nodes with this color have been identified as belonging to the first identified adjustment set, used to calculate causal effect estimates.
    -   [Light Gray]{style="color: #D4C7C7; font-weight: bold;"}: Any nodes with this color have been identified as belonging to the second identified adjustment set, used to calculate causal effect estimates.
    -   [Red]{style="color: #C00000; font-weight: bold;"}: Any nodes with this color have been flagged as introducing bias into the results of the input model. Nodes will change from (dark/light) gray to red if the model ATE falls outside the 95% confidence interval for a given adjustment set.
-   Histogram Colors
    -   [Blue]{style="color: #5C9AFF; font-weight: bold;"}: This is the 'treated' or 'success' portion of the data. Data falling within this range are categorized as 1. This represents the presence of a treatment, success, category of interest, etc...X Variable of interest, selected by the user in step 2. 
    -   [Violet]{style="color: #5C6CFF; font-weight: bold;"}: Y Variable of interest, selected by the user in step 2.
    -   [Gray]{style="color: #EAE1D7; font-weight: bold;"}: This is the 'untreated' or 'fail' portion of the data. Data falling within this range are categorized as 0. This represents the absence of a treatment, success, category of interest, etc...
-   Ribbon Plot Colors
    -   [Red]{style="color: #C00000; font-weight: bold;"}: The range of effect sizes that are outside the 95% confidence intervals of both adjustment sets. Values falling in this range are considered non-significant.
    -   [Yellow]{style="color: #FCB514; font-weight: bold;"}: The range of effect sizes that are outside of one 95% confidence interval for one effect size, but inside for the other. Values falling in this range are considered suspect, and should be closely monitored for signs of bias.
    -   [Green]{style="color: #378855; font-weight: bold;"}: The range of effect sizes that are inside both adjustment set's 95% confidence intervals. Values falling in this range are consistent with those calculated by AIR's causal estimation and are considered bias-free.
    -   [Dark Gray]{style="color: #9394A2; font-weight: bold;"}: The dot (causal estimate) and lines (95% confidence interval) drawn in this color correspond to the first identified adjustment set in your data. 
    -   [Light Gray]{style="color: #D4C7C7; font-weight: bold;"}: The dot (causal estimate) and lines (95% confidence interval) drawn in this color correspond to the second identified adjustment set in your data. 
 
# Logs {layout-ncol=2}
:::: {.columns}
::: {.column width=100%}
```{r running-log, echo=FALSE}
uiOutput("live_log_tail", 
         style = "height: 100%; min-height: 400px; display: flex; flex-direction: column;")

```
:::


::: {.column width=100%}
```{r logs-dl, echo=FALSE}
uiOutput("selected_run")
downloadButton("download_selected", "Download selected log")
# Full-width table below
uiOutput("log_table")

```
:::
::::

```{r backend-compute}
#| context: server


output$ui_reset_btn <- renderUI({
  req(graph_complete())
  actionButton(
    "reload_btn", "↻ Reset Selections",
    class = "btn btn-light btn-sm",
    style = "width:100%; color:#212529;"
  )
})
# Live log tail tracker ----
output$live_log_tail <- renderUI({
  logfile <- tryCatch(
    get(".airtool_log_file", envir = .GlobalEnv),
    error = function(e) NULL
  )

  invalidateLater(1000)

  if (is.null(logfile) || !file.exists(logfile)) {
    return(pre("No active log file found."))
  }

  log_tail <- tryCatch(
    tail(readLines(logfile, warn = FALSE), n = 25),
    error = function(e) "<unable to read log file>"
  )

  msg <- paste(log_tail, collapse = "\n")
  msg <- gsub("ERROR", "<b style='color:red'>ERROR</b>", msg, fixed = TRUE)
  msg <- gsub("WARN",  "<b style='color:orange'>WARN</b>", msg, fixed = TRUE)
  msg <- gsub("INFO",  "<b style='color:green'>INFO</b>", msg, fixed = TRUE)

  HTML(paste0(
    "<div style='height:100%; min-height:100%; display:flex; flex-direction:column;'>",
      "<div style='flex:1; height:100%; overflow-y:auto; ",
                  "background-color:#ddd; color:#111; ",
                  "font-family:monospace; font-size:0.9em; padding:10px; border-radius:6px;'>",
        "<strong>Live log:</strong><br><br>",
        "<pre style='white-space:pre-wrap; margin:0; height:100%;'>",
           msg,
        "</pre>",
      "</div>",
    "</div>"
  ))
})


### file input text --------------------------------------
output$file1_txt <- renderUI({
  req(input$file1)
  HTML(paste0("<i>Found ", nrow(df()), " instances of ", ncol(df()), " variables</i>"))
})

output$file2_txt <- renderUI({
  req(input$file2)
  HTML(paste0("<i>Found ", nrow(knowledge()), " variables with ", length(unique(knowledge()$level)), " knowledge tiers</i>"))
})
### log table stuff -----------------------------------
  log_root   <- "/workspace/logs"          # parent dir of run-stamped folders
  display_tz <- "America/Pittsburgh"          # adjust if desired

  # Build a data.frame of runs with times (newest first)
  get_runs_df <- function() {
    entries <- dir(log_root, full.names = TRUE, recursive = FALSE, all.files = FALSE)
    dirs <- entries[file.info(entries)$isdir]
    if (!length(dirs)) {
      return(data.frame(run=character(), created=as.POSIXct(character()),
                        modified=as.POSIXct(character()), size_mb=numeric()))
    }

    runs <- basename(dirs)

    # Parse "created" from folder name like "YYYY-mm-dd_HH-MM-SS"
    created <- suppressWarnings(as.POSIXct(runs, format = "%Y-%m-%d_%H-%M-%S", tz = display_tz))

    # Look for a log file inside each folder (use airtool.log, else newest *.log)
    get_log_info <- function(d) {
      p1 <- file.path(d, "airtool.log")
      cand <- if (file.exists(p1)) p1 else {
        lf <- list.files(d, pattern="\\.log$", full.names=TRUE)
        if (length(lf)) lf[which.max(file.info(lf)$mtime)] else NA_character_
      }
      if (is.na(cand)) return(list(mtime=as.POSIXct(NA), size_mb=NA_real_))
      fi <- file.info(cand)
      list(mtime = fi$mtime, size_mb = round(fi$size / 1024^2, 2))
    }

    info <- lapply(dirs, get_log_info)
    modified <- as.POSIXct(sapply(info, `[[`, "mtime"), origin="1970-01-01", tz=display_tz)
    size_mb  <- as.numeric(sapply(info, `[[`, "size_mb"))

    df <- data.frame(run = runs, created = created, modified = modified, size_mb = size_mb,
                     stringsAsFactors = FALSE)

    # Order: newest first by created (fallback to modified)
    ord_key <- ifelse(is.na(df$created), as.numeric(df$modified), as.numeric(df$created))
    df[order(df$run, decreasing = TRUE), , drop = FALSE]
  }

  # Render the table with a radio column
output$log_table <- renderUI({
  df <- get_runs_df()
  if (!nrow(df)) return(tags$p("No logs found."))

  # Determine which run should be selected:
  # - current selection if any
  # - otherwise default to the first (newest)
  selected <- isolate(input$run_choice)
  if (is.null(selected)) selected <- df$run[1]

  rows <- lapply(seq_len(nrow(df)), function(i) {
    run <- df$run[i]
    tags$tr(
      tags$td(class = "radio",
        tags$input(
          type = "radio",
          name = "run_choice",          # same group name for all
          value = run,
          checked = if (identical(run, selected)) TRUE else NULL,
          onclick = sprintf(
            "Shiny.setInputValue('run_choice', %s, {priority:'event'});",
            jsonlite::toJSON(run, auto_unbox = TRUE)
          )
        )
      ),
      tags$td(tags$code(run)),
      tags$td(ifelse(is.na(df$created[i]), "", format(df$created[i], tz = display_tz, usetz = TRUE))),
      tags$td(ifelse(is.na(df$modified[i]), "", format(df$modified[i], tz = display_tz, usetz = TRUE))),
      tags$td(ifelse(is.na(df$size_mb[i]), "", sprintf('%.2f', df$size_mb[i])))
    )
  })

  # Only initialize the Shiny value once (on first render when there isn't one)
  init_script <- if (is.null(isolate(input$run_choice))) {
    tags$script(HTML(sprintf(
      "Shiny.setInputValue('run_choice', %s, {priority:'event'});",
      jsonlite::toJSON(selected, auto_unbox = TRUE)
    )))
  } else NULL

  tagList(
    init_script,
    tags$table(
      class = "loglist",
      tags$thead(
        tags$tr(
          tags$th("Select"),
          tags$th("Run"),
          tags$th("Created"),
          tags$th("Last modified"),
          tags$th("Size (MB)")
        )
      ),
      tags$tbody(rows)
    )
  )
})
  # Optional: show which run is selected
  output$selected_run <- renderUI({
    req(input$run_choice)
    tags$p(HTML(sprintf("Selected run: <code>%s</code>", htmltools::htmlEscape(input$run_choice))))
  })

  # Download the selected run's log
output$download_selected <- downloadHandler(
  filename = function() {
    shiny::req(input$run_choice)
    paste0("airtool-", input$run_choice, ".log")
  },
  contentType = "text/plain; charset=utf-8",
  content = function(file) {
    shiny::req(input$run_choice)

    log_root <- "/workspace/logs"
    run_dir  <- file.path(log_root, input$run_choice)

    shiny::validate(shiny::need(
      grepl("^\\d{4}-\\d{2}-\\d{2}_\\d{2}-\\d{2}-\\d{2}$", basename(run_dir)),
      "Invalid run folder selected."
    ))

    # Prefer airtool.log; else newest *.log
    src <- file.path(run_dir, "airtool.log")
    if (!file.exists(src)) {
      cand <- list.files(run_dir, pattern = "\\.log$", full.names = TRUE)
      shiny::validate(shiny::need(length(cand) > 0, "No .log file found in the selected folder."))
      mt   <- file.info(cand)$mtime
      src  <- cand[which.max(mt)]
    }
    shiny::validate(shiny::need(file.exists(src), "Log file not found."))

    # --- retry wrapper to avoid races/overlayfs quirks ---
    copy_log_with_retries <- function(from, to, tries = 5, base_pause = 0.12) {
      for (i in seq_len(tries)) {
        ok <- try({
          # Snapshot the file as text (safer for logs than raw-bytes sizing)
          lines <- readLines(from, warn = FALSE, encoding = "UTF-8")

          con <- file(to, open = "w", encoding = "UTF-8")
          on.exit(try(close(con), silent = TRUE), add = TRUE)
          writeLines(lines, con, sep = "\n", useBytes = TRUE)
          flush(con)

          TRUE
        }, silent = TRUE)

        if (isTRUE(ok)) return(TRUE)
        Sys.sleep(base_pause * i)  # backoff
      }
      FALSE
    }

    ok <- copy_log_with_retries(src, file)
    shiny::validate(shiny::need(ok, "Could not stage log for download after several attempts."))
  }
) 
### variable declarations --------------------------------------------------------------
# --- Reactive value containers --- #
calc_complete <- reactiveVal(FALSE)
graph_complete <- reactiveVal(FALSE)
graph_ready <- reactive({
  vars_ready()
})
graph_updated <- reactiveVal(FALSE)
file_check <- reactiveVal(FALSE)


### input reactives -------------------------------------------------------------------
## ------------------------------------------------------------
## Read either .rds or .rda / .RData and return the first object
## ------------------------------------------------------------
model_in <- reactive({
  req(input$model_in)

  shiny::validate(
    shiny::need(tools::file_ext(input$model_in$name) %in% c("rds","rda","RData"),
                "Model must be .rds or .rda/.RData"),
    shiny::need(file.size(input$model_in$datapath) < 1000*1024^2,
                "Model file too large (100 MB limit)")
  )

  log_debug("Upload received: {input$model_in$name}")

  ## use the original filename to preserve the extension
  ext <- tolower(tools::file_ext(input$model_in$name))

  log_debug("Upload received: {b} → {c}  (ext = {d})",
            b=input$model_in$name, c=input$model_in$datapath, d=ext)

  tryCatch({

    mdl <- switch(
      ext,

      ## ----------------------------------------------------------------
      ## 1)  readRDS – single-object files
      ## ----------------------------------------------------------------
      "rds" = base::readRDS(input$model_in$datapath),

      ## ----------------------------------------------------------------
      ## 2)  load() – standard .rda / .RData
      ##     – if that fails we *try* readRDS as a fall-back
      ## ----------------------------------------------------------------
      "rda" = , "rdata" = {
        env <- new.env()
        obj_names <- tryCatch(
          base::load(input$model_in$datapath, envir = env),
          error = function(e) {
            message("load() failed, trying readRDS() …")
            NULL
          }
        )

        if (is.null(obj_names)) {
          ## ‘.rda’ file was actually an .rds – give it one more shot
          base::readRDS(input$model_in$datapath)

        } else if (length(obj_names) == 0) {
          stop("No objects found inside the .rda/.RData file")

        } else {
          env[[obj_names[1]]]
        }
      },

      ## ----------------------------------------------------------------
      ## 3)  anything else → hard-fail
      ## ----------------------------------------------------------------
      stop(sprintf("Unsupported model file extension: %s", ext))
    )

    log_debug("Loaded object of class {cls}: {nm}  (ext = {ext})",
              cls = paste(class(mdl), collapse = ", "),
              nm  = input$model_in$name,
              ext = ext)

    ## ------- basic sanity-check: does it respond to predict() ? ----------
    safe_data <- tryCatch({
      d <- head(df() %||% data.frame())
      if (!is.null(input$yvar) && input$yvar %in% names(d)) {
        d <- d[setdiff(names(d), input$yvar)]
      }
      d
    }, error = function(e) data.frame())

    log_debug("Testing predict() using columns: {cols}",
              cols = paste(names(safe_data), collapse = ", "))

    pred_test <- tryCatch(predict(mdl, safe_data), error = function(e) e)

    if (inherits(pred_test, "error")) {
      log_debug("Predict() failed during validation: {msg}",
                msg = conditionMessage(pred_test))
      shiny::validate(
        shiny::need(FALSE,
                    "Loaded model could not generate predictions with the uploaded data.")
      )
    } else if (is.factor(pred_test) || is.character(pred_test)) {
      log_debug("Predict() returned factor/character output; treating as classification model.")
    } else if (!is.numeric(pred_test)) {
      log_debug("Predict() returned unsupported output of class: {cls}",
                cls = paste(class(pred_test), collapse = ", "))
      shiny::validate(
        shiny::need(FALSE,
                "Model predictions must be numeric, factor, or character.")
      )
    } else {
      log_debug("Predict() validation succeeded; numeric output confirmed.")
    }

    log_debug("Model loaded ({a}): {b}", a=ext, b=input$model_in$name)
    mdl

  }, error = function(e) {
    log_debug("Error reading model: {a}", a=conditionMessage(e))
    NULL
  })
})

df <- reactive({
  req(input$file1)
  ext <- tolower(tools::file_ext(input$file1$name))

  shiny::validate(
    shiny::need(ext %in% c("csv", "xlsx"), "Data file must be .csv or .xlsx"),
    shiny::need(file.size(input$file1$datapath) < 5000 * 1024^2, "File too large (500 MB limit)")
  )

  tryCatch({
    if (ext == "csv") {
      first_line <- readLines(input$file1$datapath, n = 1)
      delim <- if (grepl("\t", first_line)) "\t" else ","

      # ← IMPORTANT: don’t force .default = col_number()
      df <- readr::read_delim(
        input$file1$datapath,
        delim = delim,
        guess_max = 10000,         # improve type guessing on wide/long files
        show_col_types = FALSE
      )

      log_debug("Read in CSV file: {msg}", msg = input$file1$datapath)
      log_debug("Detected delimiter: '{delim}'", delim = ifelse(delim == "\t", "tab", delim))
    } else {
      shinyalert::sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Unsupported file type",
        text = "Only .csv files are supported.",
        type = "error",
        btn_labels = "Continue",
        closeOnClickOutside = FALSE
      )
      return(NULL)
    }

    log_debug("Column names: {names}", names = names(df))

    # If you truly need to coerce only specific columns, do it TARGETED, not blanket:
    # df[c(1, ncol(df))] <- lapply(df[c(1, ncol(df))], readr::parse_number)

    df
  }, error = function(e) {
    log_debug("Error reading data file: {msg}", msg = conditionMessage(e))
    NULL
  })
})

high_corr_pairs <- reactive({
  d <- df()
  req(!is.null(d))

  # keep only truly numeric columns
  num <- d[, vapply(d, is.numeric, logical(1)), drop = FALSE]

  # need at least 2 numeric columns
  if (ncol(num) < 2) {
    return(data.frame(var1 = character(), var2 = character(), r = numeric()))
  }

  # drop constant / all-NA columns (cor() would return NA)
  keep <- vapply(num, function(x) {
    x2 <- stats::na.omit(x)
    length(unique(x2)) > 1
  }, logical(1))
  num <- num[, keep, drop = FALSE]
  if (ncol(num) < 2) {
    return(data.frame(var1 = character(), var2 = character(), r = numeric()))
  }

  # compute correlation
  cm <- stats::cor(num, use = "pairwise.complete.obs", method = "pearson")
  # keep only upper triangle (no dupes, no diagonal)
  cm[lower.tri(cm, diag = TRUE)] <- NA

  thr <- 0.9  # make this an input if you want
  idx <- which(abs(cm) >= thr, arr.ind = TRUE)
  if (!nrow(idx)) {
    return(data.frame(var1 = character(), var2 = character(), r = numeric()))
  }

  data.frame(
    var1 = rownames(cm)[idx[, 1]],
    var2 = colnames(cm)[idx[, 2]],
    r    = cm[idx],
    row.names = NULL,
    stringsAsFactors = FALSE
  )[order(-abs(cm[idx])), ]
})

observe({
  pairs <- high_corr_pairs()
  req(nrow(pairs) > 0)

  # Nice HTML list for readability inside SweetAlert
  items <- paste0(
    "<ul>",
    paste(sprintf("<li><b>%s</b> &amp; <b>%s</b> (r = %.2f)</li>", pairs$var1, pairs$var2, pairs$r),
          collapse = ""),
    "</ul>"
  )

if (nrow(pairs) > 0) {
  # Sort and format nicely
  pairs <- pairs[order(-abs(pairs$r)), ]
  
  fmt_num <- function(x) formatC(x, format = "f", digits = 4)
  esc      <- function(x) htmltools::htmlEscape(x)
  
  items <- paste0(
    "<ul style='text-align:left;
                list-style-type:disc;
                font-family:ui-monospace, SFMono-Regular, Menlo, Monaco, Consolas, 
                             \"Liberation Mono\", \"Courier New\", monospace;
                line-height:1.3;'>",
    paste(
      sprintf(
        "<li><span style='display:inline-block; width:4.5em; text-align:right;'><b>%s</b></span>  [%s -- %s]</li>",
        fmt_num(pairs$r),
        esc(pairs$var1),
        esc(pairs$var2)
      ),
      collapse = ""
    ),
    "</ul>"
  )
  shinyWidgets::sendSweetAlert(
    session = shiny::getDefaultReactiveDomain(),
    title   = "⚠️ High Correlation Detected",
    text    = htmltools::HTML(paste0(
      "<p><b>Some variables in your dataset are strongly correlated (|r| ≥ 0.90).</b></p>",
      "<p>This may indicate redundant information and can affect downstream analyses such as regression or causal estimation.</p>",
      "<hr style='margin:8px 0;'>",
      "<p><b>Highly correlated pairs:</b></p>",
      items
    )),
    type = "warning",
    btn_labels = "OK",
    html = TRUE
  )
}



})


knowledge <- reactive({
  req(input$file2)

  shiny::validate(
    shiny::need(tolower(tools::file_ext(input$file2$name)) %in% c("csv", "xlsx"),
                "Knowledge file must be .csv or .xlsx"),
    shiny::need(file.size(input$file2$datapath) < 5 * 1024^2,
                "Knowledge file too large (5 MB limit)")
  )

  tryCatch({
    path2 <- input$file2$datapath
    ext <- tolower(tools::file_ext(path2))
    log_debug("Read in knowledge file: {p}", p = path2)

    # --- Detect delimiter if CSV ---
    if (ext == "csv") {
      first_line <- readLines(path2, n = 1)
      delim <- if (grepl("\t", first_line)) "\t" else ","

      knowledge_in <- readr::read_delim(
        path2,
        delim = delim,
        col_names = TRUE,
        show_col_types = FALSE
      )

      log_debug("Detected delimiter: '{delim}'", delim = ifelse(delim == "\t", "tab", delim))
    } else  {
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Unsupported file type",
        text = "Only .csv files are supported.",
        type = "error",
        btn_labels = "Continue",
        closeOnClickOutside = FALSE
      )
      return(NULL)
    }

    # --- Column validation and auto-fix ---
    if (!all(colnames(knowledge_in) %in% c("level", "variable"))) {
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Warning: knowledge file not formatted properly",
        text = "Attempting to fix header and assign columns.
                Please fix knowledge file formatting to ensure proper functionality.",
        type = "warning",
        btn_labels = "Continue",
        closeOnClickOutside = FALSE
      )
    }

    fixed <- fix_knowledge(knowledge_in)

    shiny::validate(
      shiny::need(all(c("level", "variable") %in% names(fixed)),
                  "Knowledge file must have ‘level’ and ‘variable’ columns")
    )

    if (is.character(fixed)) {
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Upload Failed",
        text = "An error occurred while fixing the knowledge file.",
        type = "error",
        btn_labels = "Continue",
        closeOnClickOutside = FALSE
      )
      return(NULL)
    }

    log_debug("Knowledge file names: {msg}", msg = fixed$variable)
    log_debug("Knowledge file levels: {lvl}", lvl = fixed$level)
    file_check(TRUE)

    fixed
  }, error = function(e) {
    log_debug("Error in knowledge(): {msg}", msg = conditionMessage(e))
    sendSweetAlert(
      session = shiny::getDefaultReactiveDomain(),
      title = "Upload Failed",
      text = "An error occurred while reading the knowledge file.",
      type = "error",
      btn_labels = "Continue",
      closeOnClickOutside = FALSE
    )
    NULL
  })
})

cols_mismatch <- reactive({
  req(input$file2)

  matches <- !all(knowledge()$variable %in% names(df()))
  return(matches)
})

vars_ready <- reactive({
  nzchar(input$xvar) && nzchar(input$yvar)
})

x_ready   <- reactive({ nzchar(vars_sel()$x) })         # X chosen & valid
y_ready   <- reactive({ nzchar(vars_sel()$y) })         # Y chosen & valid
both_ready <- reactive({ x_ready() && y_ready() })      # X + Y ready

tv_dir <- reactive({
  input$tv_dir_in
})

ov_dir <- reactive({
  input$ov_dir_in
})

tv_val <- reactive({
  req(x_ready(), !is.null(input$tv_threshold))          # ← new guard
  safe_thresh(input$tv_threshold,
              df()[[vars_sel()$x]], "X threshold")
})

ov_val <- reactive({
  req(y_ready(), !is.null(input$ov_threshold))          # ← new guard
  safe_thresh(input$ov_threshold,
              df()[[vars_sel()$y]], "Y threshold")
})

vars_sel <- reactive({
  req(df())                           # data present

  list(
    x = if (!is.null(input$xvar) && nzchar(input$xvar))
           check_column_exists(df(), input$xvar, "xvar")
         else "",
    y = if (!is.null(input$yvar) && nzchar(input$yvar))
           check_column_exists(df(), input$yvar, "yvar")
         else ""
  )
})

## ───────────────────────────────────────────────────────────
## 1. highlighted_graph  – always returns a dot *with* fills
##    whenever x / y / Z change
## ───────────────────────────────────────────────────────────
highlighted_graph <- reactive({
  req(graph_complete())                # we *do* need the cpdag
  req(both_ready())

  dot <- graph_to_dot(rv$best_cpdag_seen_so_far)

  Z <- if (!is.null(rv$Zvars)) rv$Zvars
    else if (exists("Zvars", inherits = TRUE)) get("Zvars", inherits = TRUE)
    else NULL

  
  ## — colour the experimental variable if one is chosen
  if (!is.null(vars_sel()$x) && nzchar(vars_sel()$x))
    dot <- change_node_color(dot, vars_sel()$x, "'#5C9AFF'")

  ## — colour the outcome variable if one is chosen
  if (!is.null(vars_sel()$y) && nzchar(vars_sel()$y))
    dot <- change_node_color(dot, vars_sel()$y, "'#5C6CFF'")


  if (!is.null(Z) && nrow(Z)) {
      if (any(Z$grp == "Z1")) dot <- change_node_color(dot, Z[Z$grp == "Z1", "Z"], "'#9394A2'")
      if (any(Z$grp == "Z2")) dot <- change_node_color(dot, Z[Z$grp == "Z2", "Z"], "'#D4C7C7'")
  }  

  dot
})

### event observations ----------------------------------------------------------------
options(shiny.error = function() {
  tb <- traceback()  # Capture the call stack
  showModal(modalDialog(
    title = "An error occurred",
    paste0("AIR Tool crashed. Refresh your browser and see log file, ", log_file, ", for details"),
    easyClose = TRUE
  ))
})

# Disconnect observer
observeEvent(input$disconnect, {
  tryCatch({
    session$close()
  }, error = function(e) {
    log_debug("Error in disconnect observer: {msg}", msg =conditionMessage(e))
  })
})

observeEvent(input$reload_btn, {
  ## ── "Reset Selections" – soft-reset back to Step 2 ─────────────────────
  ## Preserves: uploaded files, parsed data, knowledge, causal graph
  ## Resets:    variable selections, thresholds, model config, all results
    req(graph_complete())        # only meaningful if a graph exists

    # ── 1. Clear all result / estimation buffers ──────────────────────────

    results_buffer(data.frame())
    results_out_buffer(data.frame())
    datafile_buffer(NULL)
    superlearner_output_buffer(list())
    tmle_buffer(data.frame())
    ml_buffer(data.frame())
    combined_buffer(data.frame())
    classifier_buffer(NULL)

    # ── 2. Reset control flags ────────────────────────────────────────────
    locked(FALSE)
    calc_complete(FALSE)
    graph_updated(FALSE)
    # NOTE: step1lock and graph_complete intentionally left alone

    # ── 3. Clear per-run reactive values (keep graph artefacts) ───────────
    rv$Zvars     <- NULL
    rv$model_yn  <- NULL
    rv$model_ate <- NULL

    # ── 4. Tidy up global-env leftovers from the previous run ─────────────
    for (nm in c("Zvars", "df_vars", "results_detailed")) {
      if (exists(nm, envir = .GlobalEnv)) rm(list = nm, envir = .GlobalEnv)
    }

    # ── 5. Reset UI inputs back to blanks ─────────────────────────────────
    updateSelectInput(session, "xvar", selected = "")
    updateSelectInput(session, "yvar", selected = "")
    # Threshold / direction inputs will vanish automatically because their

    # renderUI calls req(vars_sel()$x) / req(vars_sel()$y), which fail
    # when the selects are empty.

    # ── 6. Re-render the plain (uncoloured) causal graph ──────────────────
    output$blankGraph <- renderGrViz({
      req(graph_complete(), rv$dot)
      dot <- sub(
        "(digraph[^\\{]*\\{)",
        "\\1\n  graph [tooltip=\"Causal graph\"]",
        rv$dot, perl = TRUE
      )
      grViz(dot)
    })

    log_info("Reset Selections: state rolled back to Step 2.")

    sendSweetAlert(
      session = shiny::getDefaultReactiveDomain(),
      title   = "Selections Reset",
      text    = "Your causal graph has been preserved. Please select new X and Y variables to continue.",
      type    = "info",
      btn_labels        = "OK",
      closeOnClickOutside = TRUE
    )
})

# File input triggers reactive evaluation
observeEvent(input$file2, {
  tryCatch({
    result <- knowledge()
    # Optionally log result
    log_debug("{line}", line = paste(capture.output(str(result)), collapse = "\n"))
  }, error = function(e) {
    log_debug("Error in file2 observer: {msg}", msg = conditionMessage(e))
  })
})

# Handle the user’s choice from the SweetAlert
observeEvent(input$confirm_mismatch, ignoreInit = TRUE, {
    if (isTRUE(input$confirm_mismatch)) {
      # User clicked "Continue"
      run_calc()
    } else {
      # User clicked "Cancel" — do nothing (i.e., break out)
      showNotification("Calculation cancelled.", type = "message")
    }
  })


# Causal graph build button
observeEvent(input$buildButton, {
  if (isTRUE(cols_mismatch())) {
      bad <- htmltools::htmlEscape(setdiff(knowledge()$variable, names(df())))
      confirmSweetAlert(
        session,
        inputId = "confirm_mismatch",
        title = "Warning: Knowledge file variables don't match data headers",
        type = "warning",
        btn_labels = c("Cancel", "Continue"),
        danger_mode = TRUE,
        closeOnClickOutside = FALSE,
        html = TRUE,
        text = htmltools::HTML(
          paste0(
            "<b>Mismatches:</b> ", paste(bad, collapse = ", "),
            "<br><br>Continue anyway?"
          )
        )
      )
      # IMPORTANT: return here so we don't fall through and run
      return(invisible(NULL))
    }

    run_calc()
})

# run the graph calculation 
run_calc <- function() {
  tryCatch({
    withProgress(message = "Building Causal Graph",
                 style = "notification",
                 value = 0.1, {
      incProgress(0.2, message = "Generating cpDAG: ", detail = "Calling TETRAD API")

      if (Sys.info()[["sysname"]] != "Linux") {
        stop("Non-Linux execution path no longer supported: remove or re-implement without file I/O")
      }

      graphlist <- AIR_getGraph(df(), knowledge(), graphtxt_buffer = graphtxt_buffer)

      incProgress(0.6, message = "Generating cpDAG: ", detail = "Formatting output")
      graphtxt_raw <- graph_to_dot(graphlist[[1]])
      graphtxt_stripped <- gsub("(digraph g \\{)|\\}", "", graphtxt_raw)

      log_debug("graphtxt_stripped: {msg}", msg = graphtxt_stripped)
      if (is.null(trimws(graphtxt_stripped)) || !nzchar(trimws(graphtxt_stripped))) {
        sendSweetAlert(
          session = shiny::getDefaultReactiveDomain(),
          title = "Tetrad Returned an Empty Graph",
          text = "Please review your data to ensure that a causal relationship exists in your data",
          type = "warning",
          btn_labels = "Continue",
          btn_colors = "#3085d6",
          closeOnClickOutside = FALSE,
          showCloseButton = TRUE
        )
        stop("Tetrad returned an empty graph string")
      }

      graphtxt <- tryCatch(
        gsub("(?s)Graph Attributes:.*", "", graphtxt_raw, perl = TRUE),
        error = function(e) {
          log_debug("Regex error on graph string: {msg}", msg = conditionMessage(e))
          graphtxt_raw  # fallback to raw if regex fails
        }
      )

      count_components_from_dot <- function(dot) {
        # pattern that matches either "quoted tokens" or bare A-Za-z0-9_ tokens
        node_pat <- '(?:"[^"]+"|[A-Za-z0-9_]+)'
        edge_pat <- paste0(node_pat, "[[:space:]]*(?:->|--)[[:space:]]*", node_pat)

        # find A -> B or A -- B (allowing quotes); needs PCRE
        edges <- unlist(regmatches(dot, gregexpr(edge_pat, dot, perl = TRUE)))
        if (length(edges) == 0) return(0L)

        # normalize "A -> B" / "A -- B" into "A B"
        edges_norm <- gsub("[[:space:]]*(?:->|--)[[:space:]]*", " ", edges, perl = TRUE)

        # split into pairs and drop the quotes around node names
        edge_pairs <- do.call(rbind, strsplit(edges_norm, " "))
        edge_pairs <- gsub('^"|"$', "", edge_pairs)  # remove surrounding quotes

        # build graph directly from edgelist matrix
        g <- igraph::graph_from_edgelist(edge_pairs, directed = TRUE)
        igraph::components(igraph::as.undirected(g, mode = "collapse"))$no
      }

      subgraphs <- count_components_from_dot(graphtxt)
      log_debug("Number of subgraphs: {no}", no = subgraphs)

      rv$ts <- graphlist[[2]]
      rv$MC_passing_cpdag_already_found <- graphlist[[3]]
      rv$best_cpdag_seen_so_far <- graphlist[[4]]

      dot_raw <- graph_to_dot(graphlist[[1]])
      rv$dot <- gsub("(?m)^\\s*\\n", "", dot_raw %||% "", perl = TRUE)

      log_debug("Graph: {msg}", msg = rv$dot)
      graph_complete(TRUE)
    }) # <- closes withProgress

    output$blankGraph <- renderGrViz({
      req(graph_complete(), rv$dot)
      log_debug("Graph Built.")
      dot <- sub("(digraph[^\\{]*\\{)", "\\1\n  graph [tooltip=\"Causal graph\"]", rv$dot, perl = TRUE)
      grViz(dot)
    }) # <- closes renderGrViz

    if (subgraphs > 1) {
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Tetrad Returned Multiple Sub-graphs",
        text = "The discovered graph has more than one disconnected component. This can happen when certain variables have no causal relationship with others, or when key mediators are missing. If you expected a single connected graph, consider whether important variables are absent, or whether some variables should be treated as non-causal.",
        type = "warning",
        btn_labels = "Continue",
        btn_colors = "#3085d6",
        closeOnClickOutside = TRUE,
        showCloseButton = FALSE
      )
    } else {
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Your Causal Graph is Ready",
        text = NULL,
        type = "success",
        btn_labels = "Continue",
        btn_colors = "#3085d6",
        closeOnClickOutside = TRUE,
        showCloseButton = FALSE
      )
    } # <- closes if/else
  }, # <- ends tryCatch expr arg
  error = function(e) {
    err_summary <- tryCatch({
      paste0("Error in buildButton observer: ", conditionMessage(e), "\n")
    }, error = function(inner) {
      paste0("Unknown error in buildButton observer: ", capture.output(str(e)), collapse = "\n")
    })
    log_debug(err_summary)
  }) # <- closes tryCatch
} # <- closes run_calc()

observeEvent( {input$xvar; input$yvar}, ignoreInit = TRUE, {
  req(both_ready())
  req(graph_complete())

# ── NEW: Y must be reachable from X to keep Tetrad happy

shiny::validate({
  v   <- vars_sel()
  buf <- graphtxt_buffer()
  orp <- get_orphans(graphtxt_buffer)

  shiny::need(!is.null(v$x) && !is.null(v$y) && nzchar(v$x) && nzchar(v$y),
              "Please select both X and Y")

  # X and Y must differ
  shiny::need(isTRUE(v$x != v$y),
              "X must be a different column than Y")

  # X must NOT be an orphan (i.e., must have at least one parent)
  shiny::need(!(v$x %in% orp),
              "X must have parents")

  # Y must be a descendant of X
  shiny::need(v$y %in% get_X_descendents(v$x, graphtxt_buffer),
              "Y must be a descendant of X")
})

## ---------- 1. maintain the global helpers -------------------------------
  df_vars <- data.frame(
    var = c("TV", "OV"),
    val = c(vars_sel()$x, vars_sel()$y),
    stringsAsFactors = FALSE
  )
  assign("df_vars", df_vars, envir = .GlobalEnv)

  ### adjustment-set calculation (unchanged) ---------------------------------
  adj_list <- AIR_getAdjSets(
    rv$ts,
    vars_sel()$x,
    vars_sel()$y,
    rv$MC_passing_cpdag_already_found,
    rv$best_cpdag_seen_so_far
  )

  if (identical(adj_list, NULL)){
    sendSweetAlert(
      session = shiny::getDefaultReactiveDomain(),
      title = "Graph Error",
      text = "Adjustment sets could not be calculated because the graph estimated from your data does not satisfy the statistical conditions needed for this step.",
      type = "error",
      closeOnClickOutside = FALSE
      )
    
  } else if (adj_list[1] == "error"){
    if (adj_list[2] == "edges"){
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "No adjustment set found",
        text  = "Undirected edges detected in graph. Adjustment sets may or may not be missing
                due to undirected edges. If undirected edges exist near your variables of interest,
                consider revising knowledge file so that search result has no undirected edges
                and try again. Ensure that nodes on either side of an undirected edge are on
                separate levels in the knowledge file to ensure directionality is enforced.
                If this does not resolve the issue, there may simply be no valid adjustment
                sets (in this case, the undirected edges are likely a red herring).
                Automatic resolution of this issue is not yet supported by the AIR
                Tool, so users must resolve this issue manually.",
        type  = "error",
        closeOnClickOutside = FALSE
      )
    } else if (adj_list[2] == "structure"){
      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "No adjustment set found",
        text  = "The structure of this data is not amenable to causal inference as no adjustment sets can be found.",
        type  = "error",
        closeOnClickOutside = FALSE
      )
    }
  } else {
    Z1 <- strsplit(gsub("[{}']", "", toString(adj_list[[1]])), ", ")[[1]]
    Z2 <- strsplit(gsub("[{}']", "", toString(adj_list[[2]])), ", ")[[1]]
  
    Zvars_loc <- rbind(
      data.frame(name = vars_sel()$x, grp = "Z1", Z = Z1),
      data.frame(name = vars_sel()$x, grp = "Z2", Z = Z2)
    )
  
    rv$Zvars <- Zvars_loc
    assign("Zvars", Zvars_loc, envir = .GlobalEnv)
  
    ## ---------- 2. logging (unchanged) ---------------------------------------
    log_debug("Graph Updated.")
    log_debug("Selected X variable: {msg}", msg = vars_sel()$x)
    log_debug("Selected Y variable: {msg}", msg = vars_sel()$y)
    if (!is.null(input$tv_threshold))
      log_debug("Selected X threshold: {msg}", msg = tv_val())
    if (!is.null(input$ov_threshold))
      log_debug("Selected Y threshold: {msg}", msg = ov_val())
  
    ## ---------- 3. re-render the graph ---------------------------------------
    output$blankGraph <- renderGrViz({
      req(graph_complete())
      #input$xvar; input$yvar
      #highlighted_graph()
      grViz(get_updated_graph(
                              rv$dot,
                              vars_sel()$x,
                              vars_sel()$y,
                              rv$Zvars))    # highlighted_graph() already reacts                       
    })
    graph_updated(TRUE)    
  }
})

observeEvent(input$goButton, {
  req(both_ready())
  req(graph_complete())


shiny::validate({
  v   <- vars_sel()
  buf <- graphtxt_buffer()
  orp <- get_orphans(graphtxt_buffer)

  shiny::need(!is.null(v$x) && !is.null(v$y) && nzchar(v$x) && nzchar(v$y),
              "Please select both X and Y")

  # X and Y must differ
  shiny::need(isTRUE(v$x != v$y),
              "X must be a different column than Y")

  # X must NOT be an orphan (i.e., must have at least one parent)
  shiny::need(!(v$x %in% orp),
              "X must have parents")

  # Y must be a descendant of X
  shiny::need(v$y %in% get_X_descendents(v$x, graphtxt_buffer),
              "Y must be a descendant of X")
})

    tryCatch({
    locked(TRUE)

    print(ls(pattern = "df", envir = .GlobalEnv))

    log_debug("Selected model_exist: {msg}", msg = input$model_exist)
    log_debug("Selected ate_in: {msg}", msg = input$ate_in)

    withProgress(message = 'Building Causal Graph', style = "notification", value = 0.1, {
      incProgress(0.1, message = "Calculating Adjustment Sets", detail = "Creating Compute Environment")

      datafile_buffer(df())
      rv$model_yn <- input$model_exist
      rv$model_ate <- input$ate_in

      results_out_buffer(data.frame(
        Row = numeric(),
        Treatment = character(),
        Group = character(),
        Mean = numeric(),
        LCI = numeric(),
        UCI = numeric(),
        stringsAsFactors = FALSE
      ))

      log_debug("Calculating Adjustment Sets")

      tryCatch({
        combos <- expand.grid(i = unique(Zvars$name), j = unique(Zvars$grp), stringsAsFactors = FALSE)
        tv_threshold_val <- tv_val()
        ov_threshold_val <- ov_val()
        tv_dir_val <- tv_dir()
        ov_dir_val <- ov_dir()

        n_cores <- min(nrow(combos), parallel::detectCores() - 1L)
        cluster <- makeCluster(n_cores)


        on.exit(stopCluster(cluster), add = TRUE)

        clusterEvalQ(cluster, {
          library(tmle3)
          library(sl3)
        })

        ## Only truly-static objects go through the automatic
        ## foreach `.export` mechanism
        vars_to_export <- c(
          "Zvars",
          "runSuperLearner",
          "log_file"
        )

        ## while we *explicitly* push the big lookup table once
        ## so every worker sees it

        # new – export both objects in one vector
        parallel::clusterExport(
          cluster,
          c("df_vars",            # the lookup table
            "safe_jcall",         # shielded .jcall
            "log_debug", "log_info", # NEW
            "log_warn",  "log_error",
            "graph_to_dot",
            "get_tetrad_path"),       #     (anything else the workers touch)
          envir = .GlobalEnv
        )

        registerDoParallel(cluster)

        incProgress(0.1, detail = "Calculating ATE for multiple Adjustment Sets")

        # ---- take immutable copies while we're still in the main thread ----
        df_snapshot     <- isolate(datafile_buffer())   # <-- plain data.frame now
        graph_snapshot  <- isolate(rv$dot)              # whatever rv$dot is

        # pass the snapshots to the workers; nothing reactive inside anymore
        results_list <- foreach(k = 1:nrow(combos),
                .packages = c("AIPW", "dplyr", "e1071", "earth", "ggplot2",
                              "hal9001", "nnet", "randomForest", "readr",
                              "scales", "sl3", "tidyr", "tmle3", "xgboost",
                              "foreach", "doParallel", "shiny"),
                .export   = vars_to_export) %dopar% {
          i <- combos[k, "i"]
          j <- combos[k, "j"]
          settings <- data.frame(
            doc_title = paste0(i, "-", j),
            nfold = 20,
            Z_level = j,
            varName = i,
            confounders = paste0(Zvars[Zvars$name == i & Zvars$grp == j, ]$Z, collapse = " ")
          )

          res <- tryCatch(
            runSuperLearner(settings,
                            graph_snapshot,
                            df_snapshot,
                            tv_dir_val,
                            ov_dir_val,
                            tv_threshold_val,
                            ov_threshold_val,
                            log_file),
            error = function(e) {
              log_debug("SuperLearner: Error i={i}, j={j}: {msg}", msg = conditionMessage(e), i = i, j = j)
              NULL
            }
          )
          res
        }
  

        ## Main thread again.
        for (res in results_list) {
          if (!is.null(res)) {
              results_buffer(
              dplyr::bind_rows(results_buffer(), res$tmle)
            )
              superlearner_output_buffer(
              append(superlearner_output_buffer(), res$diag)
            )
          }
        }

        # -----------------------------------------------------------------
        # Build the final combined data-frame so the UI can see it
        # -----------------------------------------------------------------
        tryCatch({
          # “settings_stub” supplies the few fields processResults() still reads
          settings_stub <- data.frame(
            doc_title  = "Combined",
            nfold      = 20,
            Z_level    = "Z?",
            varName    = vars_sel()$x,
            confounders = ""
          )
        
          processResults(settings_stub,
                         AIRHome       = graph_snapshot,
                         tv_dir        = tv_dir_val,
                         ov_dir        = ov_dir_val,
                         tv_threshold  = tv_threshold_val,
                         ov_threshold  = ov_threshold_val,
                         model_in      = if (identical(input$model_exist, "Yes"))
                                            model_in() else NULL,
                         model_yn      = input$model_exist,
                         model_ate     = input$ate_in,
                         log_file      = log_file)
        }, error = function(e) {
          log_debug("processResults() failed: {msg}", msg = conditionMessage(e))
        })

        incProgress(0.2, message = "Processing Results")
        log_debug("Successfully closed parallel cluster")

      }, error = function(e) {
        log_debug("Fatal error during ATE estimation: {msg}", msg = conditionMessage(e))
      })

     calc_complete(TRUE)

      output$blankGraph <- renderGrViz({
        dot <- get_final_graph(rv$dot, vars_sel()$x, vars_sel()$y, rv$Zvars)
        grViz(dot)
      })

      sendSweetAlert(
        session = shiny::getDefaultReactiveDomain(),
        title = "Causal Estimates Successfully Calculated",
        text = "",
        type = "success",
        btn_labels = "Continue",
        btn_colors = "#3085d6",
        closeOnClickOutside = TRUE,
        showCloseButton = FALSE
      )
    })

  }, error = function(e) {
    log_debug("Top-level observeEvent(goButton) failure: {msg}", msg = conditionMessage(e))
  })
})

### output definitions ------------------------------------------------------------------

output$ui_graph_pane = renderUI({
  grVizOutput('blankGraph')
})

output$second_column_content <- renderUI({
    if (calc_complete()) {
      tags$div(
        style = "display: flex; flex-direction: column; height: 100%;",
        tags$div(
          style = "flex: 0 0 25%;",
          uiOutput("ui_ribbon_plot")
        ),
        tags$div(
          style = "flex: 0 0 25%;",
          uiOutput("ui_figurecaption")
        ),
        tags$hr(
          style = "border: none; border-top: 1px solid #ccc; margin: 5px 0;"
        ),
        tags$div(
          style = "flex: 1;",
          uiOutput("ui_interpretation")
        )
      )
    } else {
      tags$div(
        tags$div(
          style = "flex: 0 0 25%",
          plotOutput('histogram_x')
        ),
        tags$div(
          style = "flex: 0 0 25%",
          plotOutput('histogram_y')
        )
      )
    }
  })

output$ui_dl_btn <- renderUI({
  req(calc_complete())
  downloadBttn(
              outputId = "download_report",
              style = "simple",
              color = "primary"
            )
})

output$step3 = renderUI({
  req(graph_complete())
  tagList(
    hr(),
    h5("Step 2- Select variables:")
  )
})



# --- X -----------------------------------------------------------------------
output$xvar <- renderUI({
  req(graph_complete(), df())

  all_cols   <- names(df())
  orphans    <- get_orphans(graphtxt_buffer)
  spinsters  <- get_spinsters(graphtxt_buffer)
  base_choices <- sort(unique(setdiff(all_cols, union(orphans, spinsters))))
  blocked      <- sort(setdiff(all_cols, base_choices))

  log_debug("orphans: {orphans}", orphans = orphans)
  log_debug("spinsters: {spinsters}", spinsters = spinsters)
  log_debug("base_choices: {base_choices}", base_choices = base_choices)

  grouped_x_choices <- list(
    "Available (estimable)" = c("— Select X —" = "", setNames(base_choices, base_choices)),
    "Not estimable"         = setNames(blocked, blocked)
  )

  selected_val <- {
    prev <- isolate(input$xvar)
    if (!is.null(prev) && prev %in% base_choices) prev else ""  # nothing selected on first load
  }

  tagList(
    selectInput(
      inputId   = "xvar",
      label     = "Experimental (X) variable:",
      choices   = grouped_x_choices,  # <-- no nested optgroups
      selected  = selected_val,
      selectize = FALSE
    ),
    tags$style(HTML("#xvar option:disabled { color:#9aa0a6; }")),
    tags$script(HTML(sprintf("
      (function(){
        $(document).off('.xvarBlock');

        var blocked = %s;

        function applyXDisabled(){
          var sel = document.getElementById('xvar');
          if(!sel) return;

          $('#xvar option').each(function(){
            this.disabled = false;
            this.title = '';
          });

          if (sel.value && blocked.indexOf(sel.value) !== -1) {
            sel.value = '';
            $(sel).trigger('change');
          }

          $('#xvar option').each(function(){
            if (this.value && blocked.indexOf(this.value) !== -1){
              this.disabled = true;
              this.title = 'Not selectable (orphan/spinster)';
            }
          });
        }

        var tries=0, t=setInterval(function(){
          applyXDisabled();
          if(++tries>20) clearInterval(t);
        }, 100);

        $(document).on('shiny:value.xvarBlock shiny:visualchange.xvarBlock', applyXDisabled);
      })();
    ", jsonlite::toJSON(blocked, auto_unbox = TRUE))))
  )
})

# --- Y -----------------------------------------------------------------------
output$yvar <- renderUI({
  req(x_ready(), graph_complete(), df())

  step1lock(TRUE)

  # freeze X for this render; don't re-run on Y changes
  xcol <- isolate(vars_sel()$x)
  req(!is.null(xcol), nzchar(xcol))

  # candidates for Y (descendants of X, not X)
  descs <- get_X_descendents(xcol, graphtxt_buffer)
  descs <- sort(setdiff(unique(descs), xcol))

  if (!isTRUE(locked())) {
    shiny::validate(
      shiny::need(length(descs) > 0, "Graph not ready yet – build the causal graph first")
    )
  }

  all_cols   <- names(df())
  blocked_y  <- sort(setdiff(all_cols, descs))  # visible but not selectable

  if (isTRUE(locked())) {
    # LOCKED: show only canonical Y and disable control
    selected_val <- isolate(vars_sel()$y)
    if (is.null(selected_val)) selected_val <- ""
    choices <- if (nzchar(selected_val)) selected_val else ""

    tagList(
      selectInput(
        inputId   = "yvar",
        label     = "Outcome (Y) variable:",
        choices   = choices,
        selected  = selected_val,
        selectize = FALSE
      ),
      tags$script(HTML("
        var tries=0,t=setInterval(function(){
          var $sel=$('#yvar');
          if($sel.length){ $sel.prop('disabled',true); clearInterval(t); }
          if(++tries>50) clearInterval(t);
        },100);
      "))
    )

  } else {
    # UNLOCKED
    grouped_y_choices <- list(
      "Available (descendants of X)" = c("— Select Y —" = "", setNames(descs, descs)),
      "Not selectable"               = setNames(blocked_y, blocked_y)
    )

    selected_y <- {
      prev <- isolate(input$yvar)
      if (!is.null(prev) && prev %in% descs) prev else ""  # nothing selected on first load
    }

    tagList(
      selectInput(
        inputId   = "yvar",
        label     = "Outcome (Y) variable:",
        choices   = grouped_y_choices,
        selected  = selected_y,
        selectize = FALSE
      ),
      tags$style(HTML("#yvar option:disabled { color:#9aa0a6; }")),
      tags$script(HTML(sprintf("
        (function(){
          // Remove previous handler to prevent accumulation
          $(document).off('.yvarBlock');

          var blocked = %s;

          function applyYDisabled(){
            var sel = document.getElementById('yvar');
            if(!sel) return;

            // Re-enable ALL first, then disable current blocked set
            $('#yvar option').each(function(){
              this.disabled = false;
              this.title = '';
            });

            if (sel.value && blocked.indexOf(sel.value) !== -1) {
              sel.value = '';
              $(sel).trigger('change');
            }

            $('#yvar option').each(function(){
              if (this.value && blocked.indexOf(this.value) !== -1){
                this.disabled = true;
                this.title = 'Not selectable (not a descendant of X)';
              }
            });
          }

          var tries=0, t=setInterval(function(){
            applyYDisabled();
            if(++tries>20) clearInterval(t);
          }, 100);

          $(document).on('shiny:value.yvarBlock shiny:visualchange.yvarBlock', applyYDisabled);
        })();
      ", jsonlite::toJSON(blocked_y, auto_unbox = TRUE))))
    )
  }
})

# ui: useShinyjs()
#observeEvent(locked(), {
#  session$onFlushed(function() {
#    if (isTRUE(locked())) {
#      shinyjs::disable("xvar"); shinyjs::disable("yvar")
#      runjs("$('#xvar,#yvar').closest('.form-group').css('opacity',0.5)")
#    } else {
#      shinyjs::enable("xvar"); shinyjs::enable("yvar")
#      runjs("$('#xvar,#yvar').closest('.form-group').css('opacity',1)")
#    }
#  }, once = FALSE)
#}, ignoreInit = TRUE, priority = 1000)

# Snap X/Y back if they change while locked
#observeEvent(input$xvar, {
#  if (isTRUE(locked())) {
#    sel <- isolate(vars_sel()$x); if (is.null(sel)) sel <- ""
#    isolate(updateSelectInput(session, "xvar", selected = sel))
#  }
#}, ignoreInit = TRUE, priority = 1000)

# ── Reset Y and downstream state when X changes ──────────────────
observeEvent(input$xvar, ignoreInit = TRUE, {
  # Only act if graph exists and X is valid
  req(graph_complete(), nzchar(input$xvar))
  
  # Reset Y selection to force a fresh pick
  updateSelectInput(session, "yvar", selected = "")
  
  # Clear stale adjustment set state
  rv$Zvars <- NULL
  graph_updated(FALSE)
})

fine_step <- function(vec, n_bins = 30) {
  rng <- range(vec, na.rm = TRUE)
  span <- rng[2] - rng[1]
  if (span <= 0 || is.na(span) || is.infinite(span)) return(1)
  # Step = half a bin width
  step <- span / (2 * n_bins)
  signif(step, 1)
}

 # vec = df()[[vars_sel()$x]], step = fine_step(vec)
default_split <- function(vec) {
  rng <- range(vec, na.rm = TRUE)
  step <- fine_step(vec)
  raw  <- round(mean(rng, na.rm = TRUE) / step) * step
  split_threshold(vec, raw)
}

split_threshold <- function(v, candidate) {
  # If candidate would give all-same class, nudge half a bin
  too_low  <- sum(v >= candidate) == 0
  too_high <- sum(v <  candidate) == 0

  if (too_low)  candidate <- candidate + .Machine$double.eps
  if (too_high) candidate <- candidate - .Machine$double.eps

  candidate
}

make_threshold_ui <- function(id, vec, n_bins = 30) {
  rng   <- range(vec, na.rm = TRUE)
  step  <- fine_step(vec, n_bins)
  # Default threshold is mean of the range, snapped to a valid step
  raw   <- round(mean(rng, na.rm = TRUE) / step) * step
  safe  <- split_threshold(vec, raw)          # nudge if all-same split

  numericInput(id, NULL,
               value = safe,
               step  = step,
               width = "90px")
}

output$ui_threshold_x <- renderUI({
  req(vars_sel()$x)
  vec <- df()[[vars_sel()$x]]
  val <- {
          prev <- isolate(input$tv_threshold)
          if (!is.null(prev) && !is.na(prev)) prev else default_split(vec)
          }
  dir_val <- {
          prev <- isolate(input$tv_dir_in)
          if (!is.null(prev) && !is.na(prev)) prev else TRUE
          }
 
  step <- fine_step(vec)  # half a bin, n_bins = 30
  tagList(
    tags$div(
  style = "display:flex;align-items:center;gap:.5rem;flex-wrap:wrap;",
  tags$span(paste("Test the effect when ", vars_sel()$x, " is ")),

  tags$div(
    style = "display:inline-block;",
    switchInput(
      inputId  = "tv_dir_in",
      label    = NULL,
      value    = dir_val,
      onLabel  = HTML("<span class='op-mono'>&ge;</span>"),
      offLabel = HTML("<span class='op-mono'>&le;</span>"),
      inline   = TRUE
    )
  ),

  tags$div(
    style = "display:inline-block;",
    numericInput("tv_threshold", NULL, value = val, step = step, width = "90px")
  )
),
    if (locked()) tags$script(HTML("
      var tvDisableTries = 0;
      var disableTv = setInterval(function() {
        var $el = $('#tv_threshold');
        if ($el.length) {
          $el.prop('disabled', true);
          $el.css({'pointer-events':'none', 'opacity':0.5});
          clearInterval(disableTv);
        }
        tvDisableTries++;
        if (tvDisableTries > 10) clearInterval(disableTv);
      }, 100);
    "))
  )
})

output$ui_threshold_y <- renderUI({
  req(vars_sel()$y)
  vec <- df()[[vars_sel()$y]]
  val <- if (locked() && !is.null(input$ov_threshold))
           isolate(ov_val()) else default_split(vec)
  step <- fine_step(vec)  # half a bin, n_bins = 30
  tagList(
  tags$div(
  style = "display:flex;align-items:center;gap:.5rem;flex-wrap:wrap;",
  tags$span(paste("Where success is defined as ", vars_sel()$y, " being ")),

  tags$div(
    style = "display:inline-block;",
    switchInput(
      inputId  = "ov_dir_in",
      label    = NULL,
      value    = TRUE,
      onLabel  = HTML("<span class='op-mono'>&ge;</span>"),
      offLabel = HTML("<span class='op-mono'>&le;</span>"),
      inline   = TRUE
    )
  ),

  tags$div(
    style = "display:inline-block;",
    numericInput("ov_threshold", NULL, value = val, step = step, width = "90px")
  )
),
    if (locked()) tags$script(HTML("
      var ovDisableTries = 0;
      var disableOv = setInterval(function() {
        var $el = $('#ov_threshold');
        if ($el.length) {
          $el.prop('disabled', true);
          $el.css({'pointer-events':'none', 'opacity':0.5});
          clearInterval(disableOv);
        }
        ovDisableTries++;
        if (ovDisableTries > 10) clearInterval(disableOv);
      }, 100);
    "))
  )
})

output$ui_file1 <- renderUI({
  if (!locked()) {
    # Show fileInput before locking
  div(
    style = "
      border: 1px solid #ccc;
      border-radius: 8px;
      padding: 15px;
      background-color: #f9f9f9;
      margin-bottom: 10px;
    ",
    fileInput("file1", "Upload data file:", accept = ".csv"),
    uiOutput("file1_txt")  # this will show your italic text
  )
  } else {
    # Show locked/greyed label after "Conduct Experiment"
    fname <- if (!is.null(input$file1)) input$file1$name else "<none selected>"
    tags$div(
      style = "margin-bottom: 10px; color: #aaa; font-style: italic;",
      tags$span(icon("file"), "Uploaded data file: "),
      tags$b(fname)
    )
  }
})

output$ui_file2 <- renderUI({
  req(input$file1)  # Only show after file1 is uploaded
  if (!locked()) {
    div(
      style = "
      border: 1px solid #ccc;
      border-radius: 8px;
      padding: 15px;
      background-color: #f9f9f9;
      margin-bottom: 10px;
      ",
      fileInput("file2", "Upload knowledge file:", accept = ".csv"),
      uiOutput("file2_txt")  # this will show your italic text
    )
  } else {
    fname <- if (!is.null(input$file2)) input$file2$name else "<none selected>"
    tags$div(
      style = "margin-bottom: 10px; color: #aaa; font-style: italic;",
      tags$span(icon("file"), "Knowledge file: "),
      tags$b(fname)
    )
  }
})

output$histogram_x <- renderPlot({
  req(vars_sel()$x, !is.null(input$tv_threshold))   # wait for slider
  p <- get_histogram_x(df(), vars_sel()$x, tv_dir(), tv_val())
  p
  }, bg = "transparent")

output$histogram_y <- renderPlot({
  req(vars_sel()$y, !is.null(input$ov_threshold))   # wait for slider
  p <- get_histogram_y(df(), vars_sel()$y, ov_dir(), ov_val())
  p
  }, bg = "transparent")

output$ui_model_exist <- renderUI({
  req(graph_ready())
  req(graph_updated())
  selected_val <- if (locked()) isolate(input$model_exist) else "null"
  tagList(
    selectInput(
      inputId  = "model_exist",
      label    = "Do you have an existing model?",
      choices  = c(
        "Yes: I can upload it"      = "Yes",
        "Yes: I can provide an ATE" = "ATE",
        "No: Do it all for me"      = "No",
        "No: Just show me the effect" = "null"
      ),
      selected = selected_val,
      width    = "100%"
    ),
    if (locked()) tags$script(HTML("
      var modelExistDisableTries = 0;
      var disableModelExist = setInterval(function() {
        var $sel = $('#model_exist');
        var $inp = $('#model_exist-selectized');
        if ($sel.length && $inp.length && $inp.parent().hasClass('selectize-input')) {
          $sel.prop('disabled', true);
          $inp.prop('disabled', true);
          $inp.parent().css({'pointer-events':'none', 'opacity':0.5});
          clearInterval(disableModelExist);
        }
        modelExistDisableTries++;
        if (modelExistDisableTries > 10) clearInterval(disableModelExist);
      }, 100);
    "))
  )
})

output$ui_ate_upload = renderUI({
  req(graph_ready())
  req(input$model_exist)
  if (input$model_exist == "ATE") {
    val <- if (locked()) isolate(input$ate_in) else 0
    if (locked()) {
      # After calculation: show value as static text
      tags$div(
        style = "color: #999; font-style: italic; margin-top: 0.5em;",
        paste("ATE provided:", val)
      )
    } else {
      numericInput("ate_in", label = "ATE: ", value = val, min = -1, max = 1, step = 0.1)
    }
  } else {
    return(NULL)
  }
})


observeEvent(input$ate_in, ignoreInit = TRUE, {
  v <- input$ate_in
  if (is.null(v) || is.na(v)) return()

  v2 <- max(min(v, 1), -1)   # bounds
  if (!identical(v, v2)) {
    updateNumericInput(session, "ate_in", value = v2)
  }
})
output$ui_model_upload = renderUI({
  req(input$model_exist)
  if (input$model_exist == "Yes") {
    if (locked()) {
      # After calculation: show label, not file picker
      file_label <- if (!is.null(input$model_in)) {
        paste0("Uploaded model: ", input$model_in$name)
      } else {
        "No model file selected."
      }
      tags$div(
        style = "color: #999; font-style: italic; margin-top: 0.5em;",
        file_label
      )
    } else {
      # Before calculation: show actual fileInput
      fileInput("model_in", "Upload model file", accept = c(".rds", ".rda", ".model"))
    }
  } else {
    return(NULL)
  }
})

output$step4 = renderUI({
  req(graph_ready())
  req(graph_updated())
  tagList(
    hr(),
    h5("Step 3- About the model:")
  )
})


output$ui_graphViz = renderUI({
  req(calc_complete())
  h4("Causal graph")
  grVizOutput("graphViz")
})

output$ui_ci_plot = renderUI({
  req(calc_complete())  
  h4("Comparison of ATE for AIR and ML")
  plotOutput("ci_plot")
})

output$ui_ribbon_plot <- renderUI({
  req(calc_complete())
  tagList(
    h4("Comparison of ATE for AIR and ML"),
    imageOutput("ribbon_plot")
  )
})

output$ui_ribbon_plot <- renderUI({
  req(calc_complete())
  tags$div(
    style = "max-width: 600px;
            margin: 0 auto;",
    h4("Comparison of ATE for AIR and ML"),
    # make the <img> responsive
    imageOutput("ribbon_plot",
                width  = "100%",
                height = "auto")
  )
})

output$ribbon_plot <- renderImage(
  {
    req(calc_complete())
    tmp <- tempfile(fileext = ".png")
    plt <- get_ribbon_plot(AIRHome)

    # bail out early if we didn't get an actual plot
    if (!inherits(plt, "gg")) {
      return(list(           # let Shiny show the <div> instead of an image
        src  = NULL,
        alt  = "No ribbon plot to display"
      ))
    }

    ggsave(tmp, plot = plt,
           width = 6.5, height = 4, units = "in", bg = "transparent")
           list(src         = tmp,
                contentType = "image/png",
                alt         = "ATE ribbon plot",
                width       = "100%",   # let CSS handle final size
                height      = NULL)     # keep aspect ratio
  },
  deleteFile = TRUE
)

output$ui_figurecaption = renderUI({
  req(calc_complete())
  uiOutput("figurecaption")
})

output$ui_goButton = renderUI({
  req(graph_ready())
  req(graph_updated())
  tagList(
    actionBttn(
      inputId = "goButton",
      label = "Analyze Data",
      style = "simple", 
      color = "primary"
    ),
    if (locked()) tags$script(HTML("
      var tries = 0;
      var disableGoButton = setInterval(function() {
        var $inp = $('#goButton');
        if ($inp.length) {
          $inp.prop('disabled', true);
          $inp.css({'pointer-events':'none', 'opacity':0.5});
          clearInterval(disableGoButton);
        }
        tries++;
        if (tries > 10) clearInterval(disableGoButton);
      }, 100);
    "))
  )
})

output$ui_buildButton = renderUI({
  req(file_check())
  if (isTRUE(step1lock())) return(invisible(NULL))
  tagList(
    actionBttn(
      inputId = "buildButton",
      label = "Build Graph",
      style = "simple", 
      color = "primary"
    ),
    if (locked()) tags$script(HTML("
      var tries = 0;
      var disableBuildButton = setInterval(function() {
        var $inp = $('#buildButton');
        if ($inp.length) {
          $inp.prop('disabled', true);
          $inp.css({'pointer-events':'none', 'opacity':0.5});
          clearInterval(disableBuildButton);
        }
        tries++;
        if (tries > 10) clearInterval(disableBuildButton);
      }, 100);
    "))
  )
})

output$figurecaption = renderUI({
  req(calc_complete())
  caption <- get_figure_caption(AIRHome, df_vars)
  tags$div(style = "font-size:12px;", caption)
  })

output$ui_interpretation = renderUI({
  req(calc_complete())  

  result_text <- get_ui_interpretation(AIRHome, df_vars, Zvars,
                                       model_yn = rv$model_yn,
                                       n_obs    = nrow(datafile_buffer()))

  tagList(
    tags$h3("Interpreting your results:"),
    tags$div(style = "font-size:24px;", HTML(result_text))
  )
})

# -----------------------------------------------------------------
# 2) DiagrammeR graph → PDF   (SVG → PDF via librsvg)
# -----------------------------------------------------------------
save_dot_to_pdf <- function(gr,
                            path,
                            width  = 8,   # inches
                            height = 6) {

  svg_raw <- charToRaw(DiagrammeRsvg::export_svg(gr))

  # rsvg >= 2.6 has rsvg_pdf(); earlier versions only rsvg_svg() + system()
  if (isTRUE(utils::packageVersion("rsvg") >= "2.6")) {
    rsvg::rsvg_pdf(svg_raw, file = path,
                   width  = width * 72,   # points
                   height = height * 72)
  } else {
    tmp_svg <- tempfile(fileext = ".svg")
    writeBin(svg_raw, tmp_svg)
    rsvg::rsvg_pdf(tmp_svg, file = path,
                   width  = width * 72,
                   height = height * 72)
    unlink(tmp_svg)
  }
  invisible(path)
}

output$download_report <- downloadHandler(
  filename = function() paste0("AIRTool-Report_", Sys.Date(), ".pdf"),
  content = function(file) {
    req(both_ready())

    tryCatch({
      tmpdir <- tempfile("airtool_report_")
      dir.create(tmpdir, recursive = TRUE)
      on.exit(unlink(tmpdir, recursive = TRUE), add = TRUE)
  
      # ── freeze reactive values ─────────────────────────────────────────────
      vals <- isolate(list(
        xvar          = vars_sel()$x,
        yvar          = vars_sel()$y,
        df            = df(),
        tv_dir        = tv_dir(),
        ov_dir        = ov_dir(),
        tv_threshold  = tv_val(),
        ov_threshold  = ov_val(),
        Zvars         = Zvars,
        graph_update  = graph_ready()
      ))

      #highlighted_dot <- isolate(highlighted_graph())
      graph_snapshot  <- isolate(rv$dot)              # whatever rv$dot is

      ## --- rebuild the little helper table *locally* ------------------
      df_vars <- data.frame(
        var = c("TV", "OV"),
        val = c(vals$xvar, vals$yvar),
        stringsAsFactors = FALSE
      )
  
      figure_cap   <- get_figure_caption(AIRHome, df_vars)
      result_text  <- get_ui_interpretation(AIRHome, df_vars, vals$Zvars,
                                            model_yn = isolate(rv$model_yn),
                                            n_obs    = nrow(vals$df))

  
      # -----------------------------------------------------------------
      # 1) ggplot → PDF  (histograms, ribbon)
      # -----------------------------------------------------------------
      save_plot_pdf <- function(plot_obj,
                              path,
                              width  = 8,   # inches
                              height = 6) {

      grDevices::cairo_pdf(               # <-- NOTE:  file=  not filename=
          file   = path,
          width  = width,
          height = height,
          family = "DejaVu Sans",
          onefile = TRUE
      )
      on.exit(dev.off(), add = TRUE)

      print(plot_obj)
      invisible(path)
      }
  
      # ── generate images ────────────────────────────────────────────────────
      save_plot_pdf(
        get_histogram_x(vals$df, vals$xvar, vals$tv_dir, vals$tv_threshold),
        file.path(tmpdir, "xhist.pdf")
      )
      save_plot_pdf(
        get_histogram_y(vals$df, vals$yvar, vals$ov_dir, vals$ov_threshold),
        file.path(tmpdir, "yhist.pdf")
      )
      save_plot_pdf(get_ribbon_plot(AIRHome), file.path(tmpdir, "ribbon.pdf"))

      stopifnot(
        is.data.frame(vals$Zvars),
        all(c("name", "grp", "Z") %in% names(vals$Zvars))
      )


      # ---- updated ribbon / final graph PNGs ------------------------------------
      save_dot_to_pdf(
        DiagrammeR::grViz(                       # << wrap the DOT string
          #highlighted_dot
          get_updated_graph(
                            graph_snapshot,
                            vals$xvar,
                            vals$yvar,
                            vals$Zvars)
        ),
        file.path(tmpdir, "updatedgraph.pdf")
      )

      save_dot_to_pdf(
        DiagrammeR::grViz(                       # << same here
          get_final_graph(graph_snapshot,
                          vals$xvar, vals$yvar, vals$Zvars)
        ),
        file.path(tmpdir, "finalgraph.pdf")
      )
        
      # ── render report ──────────────────────────────────────────────────────
      file.copy("report.qmd", file.path(tmpdir, "report.qmd"))
  
      withr::with_dir(tmpdir, {
        quarto::quarto_render(
          input         = "report.qmd",
          output_format = "pdf",
          output_file   = "AIRReport.pdf",
          execute_params = list(
            AIRHome      = AIRHome,
            xvar         = vals$xvar,
            yvar         = vals$yvar,
            figure_cap   = figure_cap,
            result_text  = result_text
          )
        )
      })
  
      final_path <- file.path(tmpdir, "AIRReport.pdf")
      if (!file.exists(final_path))
        stop("Failed to generate report")
  
      # deliver to browser
      file.copy(final_path, file, overwrite = TRUE)
    },
    error = function(e) {
      message("PDF-build failed: ", e$message)
      shiny::showNotification(
        paste("Report generation failed:", e$message), type = "error", duration = NULL
      )
      stop(e)   # re-throw for Shiny’s 500 handler
    })
  }
  ,contentType = "application/pdf"
)
```